购买
下载掌阅APP,畅读海量书库
立即打开
畅读海量书库
扫码下载掌阅APP

1.4 作物需水量与耗水量点面尺度转换研究

作物需水量与耗水量跨越了不同的空间尺度,小到水分子在植株器官中的运动,大到整个国家甚至全球的作物需水状况,其跨度达十几个数量级。同时,由于空间尺度范围大,气候、地形、植被或土地利用、土壤水分状况等影响作物需水量与耗水量的因子具有随机性和不确定性,作物需水量与耗水量也表现出较强的空间异质性,致使大尺度作物需水与耗水研究成果不能是小尺度作物需水与耗水研究成果的简单叠加,而小尺度研究成果亦不是大尺度研究成果的简单插值或分解。因此很多学者一直在探索研究作物需水量与耗水量的空间分布特征、作物需水量与耗水量的采样策略及由点到面的尺度转换技术。由于技术设备、人力、财力等的限制,作物需水量与耗水量及其影响因子的测定通常只能在短时间、小范围测定,同时影响作物需水量与耗水量的主要因子(如区域的土壤类型、气候等)存在时空异质性,使得某一尺度下作物需水量与耗水量具有较高的尺度依赖性,在某一尺度上建立的模型或方法一般不能移植到高一级或低一级时空问题中求解,致使小尺度(田间站点数据)的各类作物需水量与耗水量计算方法运用于大尺度、长时间上将产生较大误差。因此研究如何以较小的代价获取较多的作物需水量与耗水量信息、探索不同尺度之间的作物生命需水转换方法对于提高无资料或站点稀疏地区的作物需水量与耗水量精度、优化区域灌溉制度、实现区域农业水资源可持续发展有着非常重要的理论和实际意义。无论进行作物需水量与耗水量时间尺度的研究还是空间尺度的研究,传统的方法都需以区域内大量试验数据为基础,但受生产条件、技术水平等影响,实测资料十分有限,推算得出的公式大多为经验公式,具有较强的区域限制性。因此,引入新的理论和计算方法研究作物生命需水的时空尺度特征显得尤为迫切。

由叶片到个体的尺度转换建立在单叶的蒸腾速率和不同层次的叶面积的精确测定,个体到群体的尺度转换建立在个体蒸腾量准确获得的基础上 [86] 。应用热技术直接测定个体蒸腾耗水量简化了空间尺度转换过程。根据叶片或个体尺度的测定结果,结合环境因子及作物冠层结构特征的尺度转换可得到作物单株或群体的需水量或耗水量。实际应用中,常借助液流和作物茎秆直径、茎秆截面积、叶面积及叶量和冠幅及冠层投影面积和蒸腾有关系的生物学因子,结合作物取样调查数据,推算研究区内每一个个体的耗水量和群体耗水量 [87]

不同中间因子的尺度转换结果存在差异,如Hatton等认为选择与个体大小直接有关的因子——茎秆直径作为中间转换因子效果较好 [88] 。而Cermák等认为尺度转换因子中最佳的指标为叶量,因为叶量综合了叶面积和叶干重的信息,其中叶干重是叶片辐射吸收量的线性函数 [89] 。尺度外推时的取样方法及样本数的确定也引起了研究者的注意,即怎样合理安排有限资源使之可以最大限度地反应群体耗水信息。样株选取可按完全随机或总体特征分布函数的分位数的方法确定。许多研究者认为,在液流个体间变异性最主要的驱动因子的分布函数的分位数处选取样株进行测定,可使尺度外推效果较好,尤其是在研究对象液流的个体差异较大时,随机选取样株的方法结果可能不可靠 [89] 。Bethenod等在测定玉米田耗水量时,通过随机量取研究地34株样本测定其直径,分析其频率分布,测定直径出现频率最高的植株的液流,按照分析结果,他们选择直径范围为20mm左右的8株玉米进行液流测定,测定结果与波文比测定结果较吻合 [90]

Miller等研究认为,根据研究对象自身因子和环境因子聚类的方法来选取代表株的方法,可以提高尺度转换的精度 [91] 。在研究区域内按照网格化布置样点,获取研究地土壤样品测定其质地及水分含量,并利用地统计学空间插值的方法获取空间每点的土壤含水量数据;利用LIDAR成像来提取叶面积、直径、树高和地表及冠层的海拔高度等信息,根据这些信息对研究区域进行聚类分析,将研究区域分成几个亚区,每个亚区的个体有相似的直径和土壤水分状况,假定其液流量也相似。这样只需在每个亚区选择一个代表株进行测定便可推求群体的耗水量。但对研究区域分区的信息不易取得,某种程度上限制了其应用。

由群体或农田到区域尺度的作物需水量与耗水量的尺度转换较多采用空间插值方法,或用分布式水文模型方法与遥感方法估算区域作物需水量与耗水量的空间分布。空间插值可以分为几何方法、统计方法、空间统计方法、函数方法、随机模拟方法、物理模型模拟方法和综合方法。最常用的几何方法有泰森多边形和反距离加权方法 [92] ;常用的统计方法有趋势面方法和多元回归方法。空间统计又称地质统计学,是建立在空间相关的先验模型之上,假定空间随机变量具有二阶平稳性,或者是服从空间统计的本征假设(intrinsic hypothesis),空间统计方法以Kriging及其各种变种(Cokriging等)为代表;常用的函数方法有傅里叶级数、样条函数、双线性内插、立方卷积法等。常用的随机模拟方法有高斯过程、马尔科夫过程、蒙特卡罗方法、人工神经网络方法等。每种插值方法都有优缺点,主要依赖插值点数据的特性,对某些数据拟合很好的方法可能不适合其他一些数据 [93] 。如牛振国等根据公式对气压、辐射进行坡度、坡向和高度校正,利用GIS的空间分析功能,建立了基于DEM的内蒙古半干旱鄂尔多斯高原考考赖沟流域ET 0 的分布式模型 [94] 。史海滨等根据118万km 2 面积内135个非规则分布气象站30年气象资料,通过对非规则采样信息的正则化处理,在传统最优内插估计方法基础上,进行了误差控制下的适当外推,构成了内插与外推的合理结合,采用地质统计学原理,充分利用ET 0 区域信息的特征,采用Kriging最优无偏估计方法对区域信息进行估计,绘制各月的ET 0 最优等值线图及估计精度 等值线图 [95] 。佟玲等利用石羊河流域及周边17个站的参考作物蒸发蒸腾量与海拔高程建立回归方程,基于DEM借助GIS软件得到石羊河流域参考作物蒸发蒸腾量的空间分布图 [96] 。常秀华利用辽宁抚顺地区3个气象站的年平均ET 0 与海拔高程建立的回归方程,基于DEM借助ArcView3.2软件得到该地区参考作物蒸发蒸腾量的空间分布图 [97] 。Phillips等采用Kriging方法对哥伦比亚河流域700~1000个温度、湿度、风速数据进行插值,解包含温度、湿度、风速三个变量的非线性方程组,模拟1990年三个时段(1月9—13日、4月4—8日、8月2—6日)的潜在蒸发蒸腾量 [98] 。Martínez-Cob采用普通克里格Kriging、协克里格Cokriging及改进的残差克里格方法插值西班牙东北部山区多年平均ET 0 ,总体上,三种方法在验证站的估值与观测值一致性较好,改进的残差克里格插值要比其他两个方法差些 [99] 。Ashraf等采用美国内布拉斯加州、堪萨斯州、科罗拉多州1989—1990两年17个站的日气象数据,通过Kriging、反距离平方、反距离、Cokriging方法插值ET 0 ,对比发现Kriging方法RMSE最小,且先插值气象变量后计算ET 0 比先计算ET 0 后插值的RMSE低 [100] 。Dalezios等将带有两个偏移系数的指数半变异模型作为最适合的拟合理论模型,在希腊380个50km×50km格网内进行了ET 0 的Kriging插值,并绘制了月ET 0 与年ET 0 等值线图。对比结果表明,空间统计学可用于绘制复杂地形大区域ET 0 等值线图 [101] 。Harcum等基于美国科罗拉多州17个气象站1982年5月30日—9月1日气象资料与内布拉斯加州15个气象站1983年和1984年5月30日—9月1日气象资料,采用卡尔曼滤波(kalman filter)插值法获得区域日ET 0 ,把ET 0 看作一个随机过程,考虑了时间与空间相关,并考虑了测量误差与模型误差。而月ET 0 可看作是时间上不相关的,因此可用Kriging插值月ET 0 等值线,经检验卡尔曼滤波是可行的空间插值法 [102] 。倪广恒等基于全国范围210个气象站1976—2000年逐日气象资料,应用FAO-56 Penman-Monteith公式计算各站ET 0 ,利用反距离空间加权插值方法得到全国ET 0 分布图 [103]

作物需水量与耗水量空间尺度转换方面,Zhao等根据1961—2000年逐月气象数据计算黑河流域及周边15个县的ET 0 ,基于30m分辨率的DEM并借助GIS估算黑河流域ET 0 的空间分布,采用《中国农业百科全书》春小麦不同生育阶段的平均作物系数K c ,并由ET 0 ×K c 得到春小麦的需水量ET c 的分布 [104] 。Hashmi等首先将落基山脉东边科罗拉多州6个气象站的气象数据在流域内进行插值,由于Penman-Monteith公式所需的气象资料不能完全满足,此研究采用FAO-24 Blaney-Criddle公式计算ET 0 ,然后依据DEM调整ET 0 ,依据土地利用图生成K c 分布图,两者相乘得到研究区作物需水量ET c 分布,发现依赖距离的插值方法、泰森多边形与反距离加权方法不适合此研究,Kriging法效果好一些 [105] 。Ray等把遥感与GIS结合起来估算印度中西部半干旱地区MRBC灌区的区域作物需水量,由IRS-1C遥感影像获取土壤调整植被指数(SAVI),建立作物系数与SVAI的回归方程,计算每个栅格的作物系数,用FAO-24修正的Blaney-Criddle公式计算 3个站的月ET 0 值,借助GIS软件采用反距离平方插值法得到ET 0 分布图,ET 0 与K c 相乘得到区域作物需水量分布图 [106] 。王景雷等把GIS的空间数据管理功能和地统计学的空间分析功能有机结合,用地统计学、反距离加权、多项式与基于DEM的趋势面分析法等多种插值方法进行山东省冬小麦需水量等值线图的绘制,所得结果与当地冬小麦实际用水分布状况较为吻合 [107,108] 。佟玲利用石羊河流域春小麦需水量与海拔高程建立回归方程,基于DEM借助GIS软件得到石羊河流域春小麦需水量的空间分布图,并与反距离加权插值法和样条插值法的结果比较,认为基于DEM的回归方法效果最好 [109]

FAO-56推荐的Penman-Monteith方法主要适用于单点的单一作物耗水量计算。对于区域多种作物组合的耗水量计算,首先要根据不同代表点的气象观测资料计算代表点的耗水量,然后用插值法绘制区域耗水量的分布图。但这种方法很难克服气象因素和作物耗水的空间变异所产生的较大计算误差,没有考虑多种作物组合中作物与作物间的交互作用。区域作物耗水量的空间分布受气候、地形、植被或土地利用、土壤水分状况等因素影响,不少学者采用不同的方法进行了研究,但这些研究所考虑的影响因素较少。如果研究区的高差和范围都相对较大,气象站点的分布又比较稀疏,地形所带来的气温的影响,坡度和坡向对辐射的影响,海拔高度对气压的影响都不能忽略 [93]

20世纪60年代后期,遥感技术的应用为用能量平衡法计算区域作物耗水量提供了可能。20世纪80年代以后,利用遥感作物冠层温度估算区域耗水量分布的研究变得十分活跃,并在一些发达国家得到了一定的应用。结合大尺度区域的遥感和地面气象资料,许多学者对估算作物耗水量的各种方法进行了研究。地表能量平衡方程和Penman-Monteith阻力模型是物理基础较坚实且应用最广泛的两种方法。基于地球表面能量平衡原理估算作物耗水量的遥感方法已被证明是有效的、成本相对较低的方法。1990年以来,这一方法已取得较大的发展。其原理是:在不考虑由平流引起的水平能量传输的情况下,地表单位面积的垂向净收入能量的分配形式主要包括用于大气升温的感热通量,用于水在物态转换时(如蒸发、凝结、升华、融化等)所需的潜热通量,用于地表加热的土壤热通量,还有一部分消耗于植被光合作用、新陈代谢活动引起的能量转换和植物组织内部及植冠空间的热量储存,这一部分能量通常比测量主要成分的误差还要小,常常忽略不计。利用遥感数据计算作物耗水量时,一般在分别计算出净辐射量、土壤热通量及感热通量后,将潜热通量作为余项求出,净辐射量可由遥感方法得出,土壤热通量通常是找出与净辐射量的经验关系,但求解区域的感热通量较复杂,需要区域分布的气温、地表粗糙长度及风速等地面资料,一般较难获得。因而,感热通量的精确反演,一直为遥感蒸发蒸腾模型研究领域的热点 [110]

Penman-Menteith模型综合了可以解释辐射和感热及平流传输的能量平衡与空气动力学传输方程,故有着坚实的物理基础。但该模型仍然只描述了一维垂向通量。遥感技术的应用为Penman-Menteith模型直接估算区域作物耗水量提供了一种有效途径。其基本步骤为:首先由遥感数据和具有空间代表性的地面气象资料结合Penman-Menteith模型计算得到作物的潜在耗水量;然后通过遥感信息,如地表温度、植被指数等,计算得到实际耗水量与潜在耗水量的比值系数;最后将该系数与之前所得的潜在耗水量相乘得到实际耗水量 [111-113]

常用的区域作物耗水量遥感模型的基本方法大致分类如下:

(1)单层模型。把土壤和植被作为一个整体、一个边界层来研究其传输过程,假设所有传输发生在地表有效粗糙长度。此类模型首先将感热通量用一个一维通量梯度表达式模拟,然后由能量平衡方程用“余项法”计算耗水量。

(2)附加阻力模型。由于用遥感地表温度代替空气动力学温度计算感热通量会带来误差,尤其在半干旱区与部分植被覆盖区,将得到过高的感热通量估计值 [114] 。因而,许多学者引入“附加阻力”来改善其计算准确度(需要指出,这里的“附加阻力”不同于微气象领域中的“空气动力学附加阻力” [116] )。

(3)双层模型。单层模型通过单一的表面温度、表面阻力及空气动力学阻力计算能量通量,未考虑土壤表面和植被冠层各自的能量平衡、温度及水汽压系统之间的区别;在将非均匀下垫面较为复杂的动力传送和热量传输过程进行简化的同时,也牺牲了计算精度。

为了更精确地表达作物耗水过程,利用Shuttleworth-Wallace [64] 双层模型(简称S-W模型),将非均匀下垫面的土壤表面和植被冠层作为两个边界层分别进行能量平衡计算。同时,根据能量守恒原理,从参考高度处散发的总能量通量是各层能量通量之和。以感热通量为例,从土壤表面输送到源汇层界面的感热通量加上从冠层表面输送到该界面的感热通量应该等于从该界面输送到参考高度的感热通量,在将该模型综合遥感信息进行计算时,一种思路是:遥感源汇层界面向上的通量,即在确定净辐射量和土壤热通量后,通过联解参考高度处的感热平衡方程及总能量平衡方程求得。从而实现了用一层界面向上的通量代替土壤、植被冠层向界面的二层通量。S-W模型将土壤和植被叶片看作连续的湍流输送源,被称为系列模型,但应用仍局限于农业气象的田间尺度上。Lhomme等 [116] 基于S-W模型,并假设热红外表面温度为土壤与植被冠层表面温度的加权平均值,权重因子分别为土壤和植被的覆盖率,推导出一种感热通量双层模式。Norman等 [115] 对系列模型进行了简化,提出了应用遥感数据的平行模型(简称N95模型),该模型假设土壤通量和冠层通量互相平行,土壤表面和植被冠层分别与上层大气进行独立的能量和水汽交换。Norman等认为,这种假设在半干旱地区、较低或中等叶面积指数及中等风速的情况下是成立的。研究表明,在植被稀疏且分布不均匀时,中等风速情况下,棵间土壤表面蒸发与植被冠层蒸腾只有微弱的耦合关系。这种简化的平行模型综合利用遥感数据和地面数据,易于求解,适用于半干旱区较为常见的稀疏植被条件,因而促使双层模型在区域尺度上的应用向前迈进了一步。该模型利用遥感所获多角度或单角度亮度温度反演或进行迭代计算得出表面组分温度,在将经典双层模型进行简化的同时,保持了一定的精度;而且,模型所需输入的参数较少,因此实际操作性较强。但在植被较为稀疏的情况下,土壤热辐射对净辐射的贡献取决于土壤表面温度,这时将净辐射通量在土壤和冠层间采用比尔定律进行分配,会产生较大的系统误差。Kustas和Norman [117] 对模型进行了改进,建立了机理性更强的净辐射量分配的算法;另外,当植被叶片呈如行播作物的聚集簇生时,叶面只能截留稀疏散布叶片所接受太阳辐射量的70%~80%,因而改进的模型引入一个聚集因子进行修正,同时,该因子也可修正簇生带来的冠层内部风速的改变。上述系列模式属于“分层模型”(layer model) [118] :土壤层在植被冠层之下,各源通量在冠层顶部汇合,相互耦合;而将土壤蒸发和植被蒸腾分开考虑的还有一种补丁模型(patch model),即植被呈斑块状镶嵌在裸露的土壤表面,各源通量只有与空气的垂直作用,而无相互作用,蒸腾和蒸发并列放置,不存在耦合关系,因而被称为“分块模型” [119]

常用的遥感作物耗水模型包括:

(1)SEBAL(surface energy balance algorithm for land)。SEBAL是由Bastiaansen等 [120] 提出的多步骤求取地面特征参数进而得到区域耗水量的模型。该模型利用遥感数据反演得到的地表温度,半球反照率和归一化植被指数NDVI及其相互关系得出不同地表类型的宽带地表通量(包括感热通量和土壤热通量)后,用余项法逐像元地计算区域的分布。SEBAL为基于能量平衡原理的单层模型,它无论作为计算模型还是验证模型已在许多研究中得到了成功的应用 [121] 。其主要优点是物理基础较为坚实,适合于不同的气候条件区;另外,模型可以利用各种具有可见光、近红外和热红外波段的卫星遥感数据,并结合常规地面资料(如风速、气温、太阳净辐射等)计算能量平衡的各分量,因而可得出不同时空分辨率的分布图。运用空气和地表温度的实验方程,消去难以获取的空气温度参数。但模型在计算感热通量时,需要在遥感图像中选择干(热)点和湿(冷)点,某些情况下可能不易选择。在参数化时,对表面粗糙度的物理过程描述不充分,例如,模型只将下垫面粗糙度表示为叶面积指数的函数,未考虑地表结构,如植被行距、形状及分布状态等对其的影响;模型中仍应用了一些经验关系,在使用时需要进行率定,如宽带行星反照率与大气顶反照率的关系;模型仅适合在晴朗无云的天气和植被茂密的平原地区应用 [92]

(2)植被指数-温度梯形模型。由于作物耗水直接影响了作物与土壤的热量平衡,因此也影响了作物和土壤本身的温度。作物耗水与表面温度的这种依赖关系,为利用热红外遥感信息推求作物耗水量提供了一种途径。Jackson等利用Penman-Monteith公式和能量平衡单层模型估算得出能代表最小和最大作物耗水速率的叶面最高和最低温度。然后将这些值与实际叶面温度进行比较,得出作物实际耗水与潜在耗水之比。基于此,提出了作物水分胁迫指数(crop water stress index,CWSI),该指数在灌溉制度的制定、作物估产及植物病虫害监测等领域有着广泛的应用 [94] 。但计算该指数需要叶温,而大多数航天航空遥感器所能测得的是地表的混合温度。因此CWSI仅限于在农田等植被密集区应用。为克服CWSI的弱点,Moran等 [111] 基于能量平衡双层模型,引入水分亏缺指数(water deficit index,WDI),提出植被指数-温度梯形模型(vegetation index/temperature trapezoid,VITT),成功地将Jackson等提出的作物水分胁迫模型扩展到部分植被覆盖的区域。Yang等 [122] 在所选甘蔗试验田,通过测定不同土壤水分和作物状况下的光谱反射率及地表温度确立VITT梯形。定义土壤水分有效性指数(Ma=1-WDI),从而将水量平衡方程和遥感能量平衡方程联系起来。试验表明,两者所得Ma具有很好的相关性。Yang等认为,LANDSAT TM数据结合VITT模型是局地尺度作物耗水估算的一种较为实用的方法。姜杰等 [113] 采用LANDSAT7 ETM+遥感数据,与地面同步观测相结合,利用VITT模型,通过Penman-Monteith模型定义地气温差与植被覆盖度的理论边界,并运用修正的Jackson正弦关系,计算得到华北平原太行山山前平原农区的日蒸发蒸腾量频率分布。研究者将设在农区的试验站所在30m×30m像元的日耗水量与站内的蒸渗仪观测值比较后,发现模拟值与实测值具有较好的一致性。虽然VITT模型在非均匀下垫面条件下的区域作物耗水分布图计算中已显示出一定的优势,但在实际应用中仍存在一些不足:一方面,依据散点图进行梯形包络线范围的确定时,人为因素对此影响较大,应用时需要借助Penman-Monteith模型与特定日期、特定时刻的气象资料确立梯形的理论边界;另外,需要了解研究区域的土地利用类型;另一方面,Penman-Monteith模型的输入参数,如饱和水汽压、空气温度、风速等,通常不易获取,也很难从测定“点”扩展到与遥感数据对应像元尺度的“面”。

(3)半经验模型。半经验模型一般是用从遥感信息容易获得的参数得出瞬时或日蒸发蒸腾量。其中,应用最广泛的是由瞬时地表-空气温差得出日耗水量的简化模型(simplified method) [123] 。该模型最早由Jackson等于 1977 年提出,以后又有许多学者(如Seguin、Caselles等)对其进行研究和改进。尽管模型所需输入的参数较少,但可以很好地模拟有不同植被覆盖度的各种下垫面的复杂耗水机理。该模型的不足之处在于:在自然界,一天之内不时有云层穿过,风速变化,此时各气象因子会随之发生改变,使得地表温度的日变化幅度较大。因此用地方时13:00时的瞬时地表温度得出的日通量,其精度会受到影响 [124]

遥感作物耗水量的估算主要是利用可见光、近红外及热红外波段的反射和辐射信息及其变化规律进行相关地表参数的反演后,结合近地层大气的风速、温度和湿度等信息估算而来。近些年来,虽然在非均匀及稀疏植被下垫面能量传输机制的研究方面取得了较大的进展,但在遥感信息与作物耗水机理模型的链接中仍存在如下问题:

(1)地表温度的反演问题。热红外传感器探测的是地表辐射温度,又称为地球表面的“皮肤”温度(skin temperature)。然而,地表远非“皮肤”状或均一的二维实体,各样的组分及其各异的几何结构均增加了地表真实温度的反演难度。作物耗水模型中利用遥感地表温度或代替较难获得的空气动力学温度计算感热通量,或进行一些参数的计算(如WDI)。因而,地表温度的反演准确度直接影响着蒸发蒸腾量估算的精确度。

(2)尺度问题。包括时间延拓和空间延拓两个方面。在将瞬时蒸发蒸腾量进行日蒸发蒸腾量的扩展时所要求的“绝对晴天”在现实中出现的几率不会很大 [125] 。空间延拓主要指作物耗水模型中所需的气象参数由点测资料标定遥感像元面的数据,进而再从像元面扩展到“区域”甚至“全球”;另外,用来进行模型结果比较的局地观测数据与计算时所利用的遥感数据的尺度也存在差异 [126] 。然而,不同尺度信息之间往往是非线性、不确定的,时空尺度的延拓应是未来的研究重点。

(3)阻力问题。“面”上的气孔阻力、表面阻力(对于植被下垫面,常称为冠层阻力)及空气动力学阻力等对于区域蒸发蒸腾量估算关键的参数仍然需要依靠冠层高度及风速等“点”上资料来推算得到平均信息。如何充分利用遥感数据建立机理性较强的辅助性的阻力模型,也是需要进一步探讨的问题。

目前,遥感方法是区域作物耗水估算及其时空分布规律研究的一种有力工具,并且随着遥感技术的不断发展,特别是多光谱、多角度、多分辨率的遥感影像的应用以及非均匀下垫面区域气象场结构研究的继续深入,地表参数反演的精度在不断提高,区域作物耗水遥感模型也在逐步完善。作物耗水遥感模型在挖掘现有遥感信息精确模拟作物耗水机理的同时,也在朝着简单、便于实际应用的方向发展。在美国爱达荷州,利用SEBAL模型绘制的季节性作物耗水图,已被用来预测对地下水系统的回补量和灌溉对熊河及蛇河流域上游流量减小的影响 [127] 。遥感技术将为区域水资源的有效利用与合理配置提供一种更加有效的途径。 iX8oVsS0X4pzORO/LM3PE4D4GDAk48/B3X5XL/vY8sr8XbDV22EiEyB2xhQ3gLaU

点击中间区域
呼出菜单
上一章
目录
下一章
×