2024, Vol. 41 
2. 中国科学院大学 北京燕山地球关键带国家野外科学观测研究站，北京 101408;
3. 中国农业科学院农业资源与农业区划研究所，北京 100081
2. Beijing Yanshan Earth Critical Zone National Research Station, University of Chinese Academy of Sciences, Beijing 101408, China;
3. Institute of Agricultural Resources and Regional Planning, Chinese Academy of Agricultural Sciences, Beijing 100081, China
黄河流域是我国重要的生态安全屏障，生态环境十分脆弱，2019年9月18日，习近平总书记在黄河流域生态保护和高质量发展座谈会上指出，黄河流域生态保护和高质量发展是我国现阶段的重大国家战略，同时也是事关中华民族伟大复兴的千秋大计[1]。黄河及其支流作为北方地区最大的供水源，承载着全国15%的耕地灌溉以及12%的人口供水需求[2]。然而黄河流域大部分位于干旱、半干旱气候区，年平均降水量低于500 mm，干旱的气候条件使黄河流域水资源条件先天不足，土壤水资源匮乏。所以，土壤水分时空分布状况是制约黄河流域农业生产生活等经济发展的关键因素之一，同时也是影响黄河流域水资源循环和分配的重要环节。因此，研究黄河流域长时间尺度的土壤水分分布状况及时空变化格局，探索土壤水分对生态环境的响应机制，对黄河流域可持续发展具有极大的科学价值和现实意义[2]。
近些年来，许多专家对黄河流域诸多环境要素的时空变化进行了研究[3-7]，但大多集中在气候、生态变化等方面。袁丽华等[3]利用MOD13Q1数据研究黄河流域2000—2010年归一化植被指数(normalized difference vegetation index，NDVI)的时空变化状况，发现自2005年以来NDVI表现出快速增长的趋势，同时近年来植被覆盖得到极大的改善；曹云等[4]利用2000—2020年的黄河流域气象站长时间观测数据和中分辨率成像光谱仪(moderate resolution imaging spectroradiometer，MODIS)数据，分析黄河流域植被生态质量的时空变化特征，探讨气候变化背景下极端气温指标和极端降水指标与植被生态质量指数的变化关系；王有恒等[6]基于黄河流域142个气象观测站1961—2020年的气候观测资料分析流域气温、降水的变化趋势；巩杰等[5]利用黄河流域气象站点资料分析近61年黄河流域干湿状况时空演变格局。也有相关研究对黄河流域的水文状况时空分布特征进行分析。邵晓梅和严昌荣[8]利用土壤湿度站点资料分析黄河流域半湿润偏旱区的土壤水分动态变化规律，发现土壤水分变化范围在12.69%~18.84%;史建国等[9]利用气象站点资料计算水分亏缺量来研究黄河流域水分变化特点及时空格局; 王毅荣[10]利用玛曲1985—2003年的土壤湿度站点资料分析甘南草原地区的水热配置状况和土壤水分变化特征，发现黄河上游牧草区土壤水分存在显著下降趋势，且在1995年发生了突变; 高科等[11]利用黄河源区25 km的多源卫星数据融合全球地表土壤水分遥感数据分析不同季节影响土壤水分分布和空间变异的主要影响因子; 王晓峰等[12]利用全球陆地蒸散发阿姆斯特丹模型产品的月尺度土壤水分数据探讨黄河中游地区生态恢复前后的水文时空变化特征，其中在生态恢复区土壤水分呈下降趋势，下降速率为每年0.001 3%。
已有研究多使用土壤水分站点监测数据[8, 10]，然而土壤水分实测资料获取困难，资料时间长度有限，并且在空间上代表性差，不能体现土壤水分的区域变化特征。遥感技术能够实现大范围数据的获取，且具有良好的时空分辨率，为长时间序列、高空间分辨率土壤水分数据的获取提供了可能。目前的长时序土壤水分产品主要是被动微波反演结果或者再分析数据，
空间分辨率较为粗糙(≥25 km)[13-14]，难以满足流域尺度时空变化研究的需求[11-12]，而现有的高空间分辨率(如1 km)土壤水分数据集多采用机器学习模型[15-16]或降尺度方法生产[17]，缺乏严谨的理论基础，往往依赖于多种辅助数据，具有较大的不确定性。鉴于黄河流域土壤水分对生态环境的重要性以及土壤水分在空间上存在强烈的变异性，有必要生产高空间分辨率的土壤水分数据实现黄河流域土壤水分的变化研究。相较于微波遥感，光学、热红外遥感可以提供更精细的观测数据，这使得高空间分辨率土壤水分的反演成为可能。Leng等[13]基于植被指数-地表温度特征空间方法实现了中国大陆土壤水分的全天候计算，平均均方根误差为0.063 m3/m3，具有良好的精度，可满足土壤水分生产和应用的需求，因此可用于黄河流域高分辨率及高精度土壤水分的反演研究。由于夏季能量循环最为活跃，土壤水分对于能量变化的响应更强烈，故选取夏季研究土壤水分变化的驱动因素更能体现土壤水分的能量响应机制。因此，本研究基于“植被指数-地表温度”特征空间模型反演土壤水分来研究黄河流域夏季的土壤水分空间分布特征和时空变化趋势，并进一步分析土壤水分时空变化的驱动因素，以期为黄河流域水源涵养与生态工程建设提供科学依据。
1 研究区与数据源 1.1 研究区研究区为黄河流域(图 1)。黄河是我国的第2大河流，发源于青海省巴颜喀拉山脉北麓，在山东省东营市流入渤海，横跨9个省区，全长约5 400 km，流域总面积约8×105 km2。黄河流域是指受黄河水系影响的地理生态区[18]，从东部沿海到西部源区地表生态环境差异较大，依次分布着落叶阔叶林带、草原带、荒漠带、青藏高原带等多种典型生态环境区。黄河流域地形起伏较大，自西向东有青藏高原、内蒙古高原、黄土高原、华北平原4个地形区，其中黄河源区位于高原地区，平均海拔较高；中部地区属于典型的黄土地貌，由于植被稀疏，水土流失严重；东部以黄河冲积平原为主，海拔较低。研究区年平均降水量约为440 mm，多年平均气温为7 ℃[5]，气候复杂多样，东部为温带大陆性气候，中部地区属半干旱气候，西部为干旱气候[3]。
|
Download:
|
| 红线内即为黄河流域。 图 1 黄河流域位置示意图 Fig. 1 Schematic diagram of the location of Yellow River basin | |
本研究使用的遥感数据为2001—2020年6—8月黄河流域的MODIS影像数据产品，主要产品参数见表 1。为使数据的空间参考和分辨率一致，利用MRT工具对MODIS数据进行影像镶嵌、投影与格式转换等预处理，并重新采样到相同的空间分辨率0.01°。为保持时相一致，将MYD13A2的NDVI产品、MYD09A1的反射率及方位角产品以及MYD15A2H的叶面积指数(leaf area index, LAI)产品利用Savitzky Golay滤波插值成逐日的栅格数据。
|
表 1 MODIS产品参数 Table 1 MODIS product parameters |
全天候地表温度数据采用中国陆域及周边逐日1 km全天候地表温度数据集(TRIMS LST)产品，该数据来源于“国家青藏高原科学数据中心”平台(http://data.tpdc.ac.cn)。TRIMS LST数据集的空间分辨率为1 km，时间分辨率为2次/d，时间范围为2000—2021年；空间范围包括我国大部分区域(72°~135°E，19°~55°N)，可用于长时间尺度土壤水分计算研究[19-20]。
土壤属性数据包括田间持水量和永久萎蔫点的水分含量，其中永久萎蔫点数据由田间持水量和植物可利用水含量差值计算可得，相关数据来源于天津大学Zhang等[21]发表的产品，空间分辨率为1 km，下载地址为(https://dataverse.harvard.edu/dataset.xhtml?PersistentId=doi:10.7910/DVN/UI5LCE)。
土地利用类型数据采用的是MODIS的MCD12Q1产品，该数据空间分辨率为500 m，时间分辨率为每年。由于国际地圈生物圈计划的全球植被分类方案数据类型最为丰富，因此本研究采用MCD12Q1的第1个数据集分类方案(LC_Type1)的数据，并进行拼接、投影转换、重采样、裁剪等处理。
气象数据采用2001—2020年6—8月的美国宇航局全球陆地数据同化系统(global land data assimilation system, GLDAS)产品数据集，主要包括短波净辐射通量、长波净辐射通量、风速、气温、相对湿度5个参数，空间分辨率为0.25°，时间分辨率为3 h。由于使用的MODIS产品为Aqua卫星数据，其过境时间为13：30，而TRIMS LST产品也是由MODIS的Aqua卫星数据作为输入而计算得到，因此为保持与上述产品的时间一致性，选取与其时间较为接近的15:00的GLDAS气象数据，并进行数据格式转换、重投影、裁剪以及重采样等预处理。为实现黄河流域土壤水分变化驱动因素分析，还采用1 km分辨率的平均月温度数据和月总降水量数据产品[22]，该数据来自于黄土高原科学数据中心(http://loess.geodata.cn/)，涵盖中国的主要陆地地区(18.2°~53.5° N，73.5°~135.0° E)，空间分辨率为0.5′(约1 km)，其中月温度数据单位为0.1 ℃，月总降水量数据单位为0.1 mm。
为对土壤水分反演结果进行精度评价，本研究搜集了黄河流域及其周边60个土壤水分观测站点2017年的土壤水分逐日数据以及土壤湿度产品。其中站点数据来源于中国气象局，使用观测深度为10 cm的土壤水分数据。土壤湿度产品采用Li等[15]提供的中国范围内1 km高质量土壤湿度数据集SMCI 1.0，数据来源于“国家青藏高原科学数据中心”平台(http://data.tpdc.ac.cn)。SMCI数据集包含2000—2020年以10 cm为间隔10层深度(10~100 cm)的日尺度土壤湿度数据，采用10 cm处的土壤湿度数据，并利用python进行数据读取、重采样、裁剪等预处理，用于土壤水分反演结果的验证分析。
2 研究方法 2.1 梯形特征空间模型本研究在“植被指数-地表温度”梯形特征空间模型的基础上实现了黄河流域的土壤水分数据集反演，具体反演方法如下[13, 23]：
在梯形特征空间模型(图 2)中，若给定像元的植被覆盖度(fraction vegetation coverage, FVC)和地表温度(land surface temperature, LST)即可确定任意地表观测对象在梯形特征空间中对应的位置，并可反演出其土壤水分。如P点处的土壤水分可由下式进行计算：
| $ \theta=\theta_{\mathrm{wp}}+\left(\theta_{\mathrm{fc}}-\theta_{\mathrm{wp}}\right) M_{0}, $ | (1) |
|
Download:
|
| 图 2 梯形特征空间示意图 Fig. 2 Schematic diagram of the trapezoidal feature space | |
式中: θfc，θwp分别为土壤田间持水量(field capacity，fc)和凋萎点(wilting point，wp)含水量；M0为土壤水分可利用率。
根据Moran等[24]的研究成果，梯形特征空间4个极值点(干燥裸土温度Ts, max，湿润裸土温度Ts, min，受水分胁迫的全植被覆盖温度Tv, max和水分状况良好的全植被覆盖温度Tv, min)对应的地表温度可以通过能量守恒方程确定：
| $ \left\{\begin{align*} T_{\mathrm{s}, \max }= & \frac{r_{\mathrm{a}}\left(R_{\mathrm{n}, \mathrm{s}}-G\right)}{C_{\mathrm{v}}}+T_{\mathrm{a}}, \\ T_{\mathrm{s}, \min }= & \frac{r_{\mathrm{a}}\left(R_{\mathrm{n}, \mathrm{s}}-G\right)}{C_{\mathrm{v}}} \frac{\gamma}{\mathit{\Delta }+\gamma}-\frac{\mathrm{VPD}}{\mathit{\Delta }+\gamma}+T_{\mathrm{a}}, \\ T_{\mathrm{v}, \max }= & \frac{r_{\mathrm{a}}\left(R_{\mathrm{n}, \mathrm{v}}-G\right)}{C_{\mathrm{v}}} \frac{\gamma\left(1+\frac{r_{\mathrm{cx}}}{r_{\mathrm{a}}}\right)}{\mathit{\Delta }+\gamma\left(1+\frac{r_{\mathrm{cx}}}{r_{\mathrm{a}}}\right)}- \\ & \frac{\mathrm{VPD}}{\mathit{\Delta }+\gamma\left(1+\frac{r_{\mathrm{cx}}}{r_{\mathrm{a}}}\right)}+T_{\mathrm{a}}, \\ T_{\mathrm{v}, \min }= & \frac{r_{\mathrm{a}}\left(R_{\mathrm{n}, \mathrm{v}}-G\right)}{C_{\mathrm{v}}} \frac{\gamma\left(1+\frac{r_{\mathrm{cp}}}{r_{\mathrm{a}}}\right)}{\mathit{\Delta }+\gamma\left(1+\frac{r_{\mathrm{cp}}}{r_{\mathrm{a}}}\right)}- \\ & \frac{\mathrm{VPD}}{\mathit{\Delta }+\gamma\left(1+\frac{r_{\mathrm{cp}}}{r_{\mathrm{a}}}\right)}+T_{\mathrm{a}}, \end{align*}\right. $ | (2) |
式中: Rn, s和Rn, v分别为土壤和植被组分的表面净辐射(W/m2)；G为土壤热通量(W/m2)；Ta为空气温度(K)；Cv为空气体积热容(J/(℃·m3))；γ为干湿表常数(kPa/℃)；VPD为饱和水汽压差(kPa)；Δ为饱和水汽压随气温变化的斜率(kPa/℃)；rcx为植被受水分胁迫状态下的冠层阻力，可以设定为2 000 s/m；rcp为植被生长水分供应充足时的冠层阻力，可以通过LAI计算；ra为空气动力学阻抗[25]。
在确定了极端位置的地表温度之后，可得到梯形特征空间模型的干湿边方程：
| $ T_{\mathrm{dry}}=\left(T_{\mathrm{v}, \text { max }}-T_{\mathrm{s}, \text { max }}\right) f+T_{\mathrm{s}, \text { max }}, $ | (3) |
| $ T_{\mathrm{wet}}=\left(T_{\mathrm{v}, \text { min }}-T_{\mathrm{s}, \text { min }}\right) f+T_{\mathrm{s}, \text { min }}, $ | (4) |
式中: Tdry表示干边的温度，Twet表示湿边的温度，f为植被覆盖度。
对P位置处，可得土壤水分可利用率M0
| $ \begin{array}{*{20}{c}} M_{0}=\frac{d_{\mathrm{MP}}}{d_{\mathrm{MN}}}= \\ \frac{\left(T_{\mathrm{v}, \text { max }}-T_{\mathrm{s}, \text { max }}\right) f+T_{\mathrm{s}, \text { max }}-T_{\mathrm{s}}}{\left(T_{\mathrm{v}, \text { max }}-T_{\mathrm{s}, \text { max }}\right) f+T_{\mathrm{s}, \text { max }}-\left[\left(T_{\mathrm{v}, \text { min }}-T_{\mathrm{s}, \text { min }}\right) f+T_{\mathrm{s}, \text { min }}\right]}, \end{array} $ | (5) |
式中：dMP和dMN分别为点M和P与点M和N之间的距离。结合式(1)和式(5)，即可算出任一像元的土壤含水量。
2.2 Sen斜率与Mann-Kendall趋势检验分析Theil-Sen Median方法又称为Sen斜率估计，是一种广泛使用的非参数统计趋势计算方法。该方法通过计算时间序列的中值可以较好地避免数据中异常值对分析结果的干扰，适用于长时间序列数据的变化趋势分析。Sen斜率估计的原理表达式[26-27]如下
| $ \beta=\text { Median }\left(\frac{x_{j}-x_{i}}{j-i}\right), \quad \forall j>i , $ | (6) |
式中：i和j为以年为单位的序列(2001≤i < j≤2020)；xi和xj分别为第i和j年的土壤水分值。若β大于零，则表明土壤水分呈现出增长的变化趋势，反之则表现为随时间降低的变化趋势。
Mann-Kendall检验(简称M-K检验)是目前应用较为广泛的一种非参数统计时间序列变化趋势检验方法[3, 28]，该方法采用双边趋势检验，根据标准正态分布特征的Z值实现变化趋势的显著性检验。在给定显著性水平α下，当|Z|>Z1-α/2时，表示时间序列变化趋势显著，反之变化不显著。Z1-α/2为在置信水平α下，标准正态函数分布表对应的值[27]。
该方法的适用范围较广，不需要样本服从特定的统计分布，能够克服异常值的影响，可以实现对研究变量时间序列数据的趋势显著性检验。因此，本研究结合Sen斜率估计与M-K趋势检验以增强方法的抗噪性，提高检验结果的稳定性。研究中显著性水平α取值0.05，对照标准正态分布表，当|Z|大于1.65、1.96和2.58时，表明趋势分别通过了信度为90%、95%和99%的显著性检验，将研究变量的变化趋势分为9个等级[27]，见表 2。
|
|
表 2 变化趋势分级表 Table 2 Trend grading table |
本文采用王劲峰和徐成东[29]提出的地理探测器方法来探究土壤水分变化的驱动因子及其影响机制的差异。空间分异性是地理学的重要特征之一，该特征是指层内方差小于层间方差的现象，地理探测器通过分析因子的空间分异性特征来获取驱动其变化的影响要素[30]。主要包括以下2个探测器：
研究利用分异及因子探测器探究某影响因子对目标变量的解释程度，模型用q值度量，其计算原理如下：
| $ q=1-\frac{\sum\nolimits_{h=1}^{L} N_{h} \sigma_{h}^{2}}{N \sigma^{2}}, $ | (7) |
式中: h为目标变量Y或影响因子X的分类或分层；Nh为要素在类型h内的样本数；N为研究区内的总样本数；σ2为全区Y值的方差。q为衡量因子影响程度的指标，其中，q值越大，表示变量X对Y的影响越大，解释程度越高，该变量的值域范围为0~1。
同时，研究利用交互作用探测器识别不同影响因子之间的交互作用对目标变量的影响程度，模型通过对比2个因子分别对目标变量单独作用时的q值之和与2个因子共同作用时的q值，判断各个因素相互作用是否会增强或减弱对目标变量Y的解释能力，即当2个因子交互作用时对目标变量的影响能力增大还是减小。
以每年夏季平均土壤水分数据作为因变量，选取气温、降水、NDVI、DEM、土地利用类型5个要素对土壤水分变化进行归因分析，各要素分级情况如图 3。
|
Download:
|
| 图 3 黄河流域驱动要素分级图 Fig. 3 Classification diagram of driving factors in Yellow River basin | |
本文基于“植被指数-地表温度”梯形特征空间模型计算得到黄河流域土壤水分数据集。为确保反演结果的可靠性，共收集黄河流域60个土壤水分观测站点数据并选取2017年具有代表性的3天(6月15日，7月15日，8月15日)对结果进行精度检验。图 4(a)为土壤水分反演结果与实测站点数据的对比图。结果表明，两者表现为较好的一致性，决定系数(R2)均高于0.5，均方根误差(root mean squared error, RMSE)范围为0.045~0.051 m3/m3，偏差(bias)为0.029~0.031 m3/m3，表明该模型在黄河流域具有良好的应用潜力。
|
Download:
|
| 图 4 土壤水分反演结果精度验证图 Fig. 4 Precision verification diagram of soil moisture inversion results | |
为进一步评估土壤水分反演结果的有效性，采用随机采样方法在黄河流域范围内生成800个像素点作为样本将反演结果与已发表的1 km SMCI土壤水分产品进行对比，如图 4(b)~4(d)所示。对3天的散点图统计分析发现，散点均较为均匀地分布在1∶1线附近，表明反演结果与SMCI产品具有较好的相关性。此外，3天的平均RMSE为0.060 m3/m3，平均偏差为0.046 m3/m3，这与Leng等[13, 23, 31]的验证结果相当。通过上述验证结果表明，本文反演的土壤水分能够提供可靠的土壤水分信息，可用于黄河流域的应用分析。
综上分析，由于能够获取的数据有限，土壤水分的验证结果可能存在一些偏差，但总体上在合理的误差范围内。本研究仅搜集到土壤下0~10 cm处的平均土壤水分验证数据，与遥感探测深度可能存在不匹配问题，但目前对土壤水分遥感反演的最佳深度问题并无统一的看法，有学者认为20 cm左右深度土壤水分与遥感资料具有较好的相关性[32]，也有部分学者认为10 cm左右或30 cm左右为最佳遥感反演土壤水分深度[33-34]，因此本研究采用的验证数据仅能为反演结果提供一定的参考。针对这一问题，在未来的研究中，可以搜集不同深度的土壤水分实测数据以及更高空间分辨率的产品进行进一步验证。
3.2 黄河流域土壤水分时空变化分析研究利用梯形特征空间模型对黄河流域2001—2020年夏季(6—8月)的土壤水分进行反演与统计分析，结果如图 5所示。从图 5(a)可以看出在2001—2020年间，黄河流域夏季年均土壤水分在0.08~0.27 m3/m3之间变化，且土壤水分具有较大的空间异质性，在黄河源区和下游地区土壤较为湿润，而在黄河中游特别是黄土高原地区相对比较干旱。
|
Download:
|
| 图 5 黄河流域土壤水分时空变化图 Fig. 5 Temporal and spatial variations of soil moisture in Yellow River basin | |
通过对黄河各流域土壤水分变化趋势分析可知(图 5(b))，整体上黄河流域在近20年来夏季平均土壤水分呈现波动上升的变化趋势，平均以每年0.000 5 m3/m3的速度增长，2017年土壤水分含量达到最高值(0.123 m3/m3)，2001年土壤水分含量最低(0.108 m3/m3)。在2014年以后，土壤水分含量显著增高，增长速度相比2014年之前的每年0.000 4 m3/m3变为每年0.000 6 m3/m3，且土壤水分呈现先减少后增加的趋势，这可能与2014年开始实施新一轮退耕还林还草工程有关。在植被恢复期间，植被生长需要吸收水分，使得土壤水分减少，而随着植被的恢复，植被的水源涵养作用增强，土壤水分会呈现增加趋势[12]。根据《退耕还林综合效益监测国家报告》，截至2014年底，长江、黄河流域中上游退耕还林工程涵养水源总价值量由2014年的3 102.14亿元/a增长为2019年的4 630.22亿元/a，说明退耕还林工程对水源涵养具有重要作用。
此外，黄河各流域夏季年均土壤水分表现为不同的变化趋势，上游地区增长速度最大，土壤水分以每年0.000 8 m3/m3的速度增长；中游地区土壤水分呈现增长趋势，但变化较为缓慢(每年0.000 1 m3/m3)；下游地区土壤水分年际变化波动较大，变化趋势不显著(R2为0.008)。王晓峰等[12]研究发现1998—2018年黄河中游土壤水分有所减少，这与我们的结果存在差异，出现该问题的原因可能与研究年限和研究季节不同有关，夏季多降雨，往往夏季土壤水分含量高于其他季节且变化波动大。此外，各流域土壤水分变化出现差异的原因一方面与上游降水增多、中下游降水减少有一定关系[6]，另一方面冰川积雪融水是黄河流域水资源的重要来源，每年对黄河径流的补给占径流的1.9%，并且随着全球气候变暖，黄河源区的冰川加速消融[6]，造成地表径流量增加，同时地表径流渗透进土壤可以增加土壤含水量。黄河中游是主要的农业种植区，且黄土高原地区水土流失严重，尽管生态恢复在一定程度上缓解了水土流失，但气候的暖干化[1]增大了地面蒸发和植被蒸腾，同时人类活动用水的加剧使得中游土壤水分增长趋势缓慢。黄河下游流域面积小且流经华北平原地区，主要为作物种植区，人类活动是土壤水分变化的主要干扰因素。
为探索土壤水分在空间上的变化趋势，研究采用Sen斜率与Mann-Kendall趋势检验结合的方法得到黄河流域土壤水分变化趋势图。从图 5(c)可以看出黄河上游地区土壤水分主要呈现增长趋势，Sen斜率最高为0.016，黄河中游北部地区特别是内蒙古高原地区土壤水分表现为减少趋势。图 5(d)结果显示黄河流域土壤水分空间变化趋势并不显著，主要表现为不显著增加和不显著减少趋势，分别占区域总面积的39.54%和58.01%，显著变化区域仅占全区的2.44%。
3.3 土壤水分变化驱动因素分析为探索黄河流域土壤水分时间变化的驱动因素，选取气温、降水、NDVI这3个因子并分别获取各因子每年夏季的平均值，根据土壤水分变化分析结果，这里将2014年作为分界点探索2014年前后各因子对土壤水分变化的影响，各因子间的关系如图 6(a)所示。结果表明，2014年之前3个因子对土壤水分时间变化的影响程度：降水量>气温>NDVI，其中降水量是土壤水分年际变化的主导因素[5]，相关系数为0.71，气温和NDVI与土壤水分的相关性次之，但两者均与土壤水分表现为负相关关系，这可能与气温升高增大地面蒸发和植被蒸腾有关。2014年后3个因子对土壤水分时间变化的影响程度：降水量=NDVI >气温，降水和NDVI对土壤水分变化起到主导作用，相关系数均为0.94，而气温与土壤水分的相关系数为0.07，对土壤水分变化的影响并不显著。其中NDVI和气温与土壤水分均呈现出正相关关系，这可能与2014年实施新一轮植树种草有关，植被恢复增强了水土保持，减少水土流失，使土壤水分增长。综上，降水和NDVI是黄河流域土壤水分时间变化的主要驱动因子，气温对土壤水分变化的影响相对较小。
|
Download:
|
| 图 6 黄河流域因子相关性分析及NDVI变化图 Fig. 6 Factor correlation analysis and NDVI variation map in Yellow River basin | |
黄河流域地域广阔，而土壤水分具有较强的空间异质性，因此有必要探索黄河流域土壤水分空间变化的影响因素。由于黄河下游各因子地理分异较小，因此利用地理探测器仅对黄河上游和中游进行变化驱动要素分析。根据因子探测结果(表 3)可知，黄河上游和下游的土壤水分空间变化受到不同要素的影响，并且不同影响要素对各流域土壤水分的影响模式和程度存在差异。根据p值，除土地利用的p值接近1，表示结果不显著外，其他影响因子均具有显著性，表明黄河流域土壤水分的分布具有空间异质性。对于黄河上游，气温和DEM对土壤水分的空间异质性影响最为显著，NDVI和降水影响相对较小，而由于黄河上游土地利用类型主要为草地，其空间变异小，因此土地利用类型对土壤水分的影响并不显著。
|
|
表 3 土壤水分的驱动因子探测 Table 3 Detection of driving factors for soil moisture |
对于黄河中游，研究中所选择的各因子对土壤水分空间变化均表现出了显著影响(p < 0.05)，其中NDVI、降水和土地利用类型是影响土壤水分空间变化的关键因子，气温和DEM的影响相对较小。一方面，黄河中游特别是黄土高原地区水土流失严重，而植被可以起到涵养水源、水土保持的作用。自1999年实施退耕还草还林工程以来，黄河流域的植被状况发生了显著变化，图 6(b)显示，黄河流域植被NDVI在2001年以来呈现显著的增长趋势，增长区域占黄河流域的87%左右，其中黄土高原地区生态恢复工程实施面积达3.33×107 hm2，植被覆盖度也提高至63.6%[35]。尽管有关研究发现黄河中游地区生态恢复区土壤水分呈下降趋势[12]，但植被恢复的水土保持作用在长时间阶段还是起到了积极的作用。根据《黄河流域水土保持公报(2020)》，相比于1999年，黄河流域水土保持率提高20.61%，而黄土高原是水土保持的主要地区(水土保持率为63.44%)，说明植被恢复对于黄河中游土壤水分变化具有重要影响。另一方面，降水对黄河中游土壤水分的变化的影响也不容忽视，黄河中游主要位于干旱半干旱地区，降水能够渗入地面，对土壤水分变化具有直接影响。
研究利用交互探测器分析不同影响因子之间的交互作用对土壤水分的影响，结果如表 4所示，任意2种影响因子组合对土壤水分变化的影响均表现为增强效果(具体表现为非线性增强或双因子增强)，并且不同因子的交互作用对上游和中游土壤水分的影响模式不同。从q值来看，各因子交互后对黄河中游的作用比上游更加显著。对于黄河上游，温度与DEM交互后对土壤水分的影响表现为双因子增强作用(非线性增强效果要大于双因子增强)，这是由于温度与DEM是2个相互关联的量，一般而言海拔越高，气温越低，因此两者的作用是一致的。其他影响因子交互后均表现为非线性增强，同时，不同影响因子交互后的非线性增强的贡献率大于单个因子独自作用时的2倍。这一现象说明各因子的交互作用均加强了对土壤水分空间分布的解释力，上游土壤水分变化是各因子综合作用的结果，但从各因子交互后对土壤水分的贡献来看，气温与NDVI交互后的影响最大，而降水与土地利用类型之间交互的影响最小。对于黄河中游，降水与NDVI或土地利用类型、NDVI与土地利用类型交互后对土壤水分的影响均为双因子增强作用，其他各因子均为非线性增强作用，其中降水与土地利用类型交互后对土壤水分的影响最大，而气温与DEM交互后的影响最小，这也与因子探测结果一致。由于NDVI和土地利用类型易受人类活动的影响，因此各因子交互后的作用减弱，表明中游人为干扰对土壤水分空间变化具有一定的影响。
|
|
表 4 土壤水分的交互因子探测结果 Table 4 Detection results of interaction factors for soil moisture |
根据地理探测器因子探测结果，黄河流域土壤水分空间分布表现出明显的空间异质性，同时影响流域上、中游土壤水分时空变化的因素有较大差异。黄河上游和中游生态以及社会环境的差异是造成这种现象的主要原因。这一定程度上说明由于土壤水分的空间变异性较大，在实际研究中仅仅利用站点实测数据进行插值来分析区域的土壤水分时空变化缺乏说服力，因此，利用遥感反演土壤水分研究黄河流域土壤水分变化和驱动要素分析具有较大的优势。未来应考虑如何提高土壤水分的反演精度，反演得到时空分辨率更高的土壤水分数据集对于区域尺度的生态水文以及气候变化研究等具有重要的应用价值。
4 结论本研究以黄河流域这一典型的生态环境区为研究对象，利用2001—2020年夏季的MODIS数据产品以及气象数据为基础，反演土壤水分含量数据，并研究该区域的土壤水分时空变化及其驱动因子。研究结果表明：
1) 基于梯形特征空间模型计算黄河流域土壤水分数据集，与土壤水分实测站点数据验证对比RMSE范围为0.045~0.051 m3/m3，与SMCI产品对比平均RMSE为0.060 m3/m3，具有良好的反演精度，能够有效地反映黄河流域的土壤水分状况。
2) 黄河流域的土壤水分含量呈现明显的空间异质性，在黄河源区和下游地区土壤水分含量较高，较为湿润；而黄河中游地区土壤含水量较低，相对比较干旱。
3) 黄河流域夏季平均土壤水分含量在2001—2020年呈波动增长趋势，平均以每年0.000 5 m3/m3的速度增长。各流域土壤水分变化不同，上游地区增长速度最快，每年增加0.000 8 m3/m3，中游地区土壤水分增长缓慢，下游地区土壤水分变化波动较大。整体上，黄河流域土壤水分含量空间变化主要表现为不显著增加和不显著减少趋势，分别占区域总面积的39.54%和58.01%，未来需要进一步加强黄河流域的水土保持工作。
4) 降水是黄河流域土壤水分含量时间变化的主导要素，2014年以来NDVI对土壤水分含量的影响增大，气温影响最小。基于地理探测器的土壤水分空间变化驱动要素研究表明，黄河上游和下游土壤水分主要受不同要素的影响。其中，气温和DEM是影响黄河上游土壤水分含量空间分布的主要因子；中游土壤水分含量的时空变化主要受到NDVI和降水量的影响，且近20年来NDVI增长区域占全流域的87%左右，为黄河流域水土保持提供了良好条件。各因子之间的交互作用对流域不同区域土壤水分含量的影响不同，各因子交互后对黄河中游的作用要强于上游，且中游土壤水分的变化主要受人为活动的影响。
本研究以特征空间土壤水分反演模型为基础，探讨黄河流域近20年来夏季土壤水分的变化情况及其影响因素，为黄河流域生态环境保护与水资源管理提供了参考。然而，本研究仍存在一些不足之处，特征空间土壤水分反演模型的实现依赖于气象数据的输入，由于数据资源有限，本研究采用的气象数据产品空间分辨率较低，可能使土壤水分反演结果产生误差，未来通过反演高精度的气象数据来提高特征空间模型的反演精度是一种可行的方案。本研究对土壤水分变化的分析局限在夏季尺度以及区域尺度，可能无法全面地反映黄河流域土壤水分的变化状况，后续将扩大研究时序和增加分区，对黄河流域土壤水分变化进行更加细致的研究。此外，在驱动因素分析过程中本研究仅探究了部分自然因素对土壤水分的影响力，然而社会经济因子如人口数量、生态环境等是否也会对土壤水分产生作用？未来需结合相关数据进行进一步探索。
| [1] |
马柱国, 符淙斌, 周天军, 等. 黄河流域气候与水文变化的现状及思考[J]. 中国科学院院刊, 2020, 35(1): 52-60. Doi:10.16418/j.issn.1000-3045.20191223002 |
| [2] |
杨胜天, 刘昌明, 王鹏新. 黄河流域土壤水分遥感估算[J]. 地理科学进展, 2003, 22(5): 454-462. Doi:10.11820/dlkxjz.2003.05.003 |
| [3] |
袁丽华, 蒋卫国, 申文明, 等. 2000—2010年黄河流域植被覆盖的时空变化[J]. 生态学报, 2013, 33(24): 7798-7806. Doi:10.5846/stxb201305281212 |
| [4] |
曹云, 孙应龙, 陈紫璇, 等. 2000—2020年黄河流域植被生态质量变化及其对极端气候的响应[J]. 生态学报, 2022, 42(11): 4524-4535. Doi:10.5846/stxb202105151273 |
| [5] |
巩杰, 高秉丽, 李焱, 等. 1960—2020年黄河流域气候干湿状况时空分异及变化趋势[J]. 中国农业气象, 2022, 43(3): 165-176. Doi:10.3969/j.issn.1000-6362.2022.03.001 |
| [6] |
王有恒, 谭丹, 韩兰英, 等. 黄河流域气候变化研究综述[J]. 中国沙漠, 2021, 41(4): 235-246. Doi:10.7522/j.issn.1000-694X.2021.00086 |
| [7] |
黄建平, 张国龙, 于海鹏, 等. 黄河流域近40年气候变化的时空特征[J]. 水利学报, 2020, 51(9): 1048-1058. Doi:10.13243/j.cnki.slxb.20200603 |
| [8] |
邵晓梅, 严昌荣. 黄河流域半湿润偏旱区土壤水分动态变化规律研究[J]. 干旱地区农业研究, 2006, 24(2): 91-95, 109. Doi:10.3321/j.issn:1000-7601.2006.02.020 |
| [9] |
史建国, 严昌荣, 何文清, 等. 黄河流域水分亏缺时空格局变化研究[J]. 自然资源学报, 2008, 23(1): 113-119. Doi:10.3321/j.issn:1000-3037.2008.01.013 |
| [10] |
王毅荣. 黄河上游牧草区水热变化及土壤水分响应[J]. 自然资源学报, 2007, 22(2): 260-267. Doi:10.3321/j.issn:1000-3037.2007.02.013 |
| [11] |
高科, 盖艾鸿, 潘韬, 等. 黄河源区土壤水分空间变异及其主控因子[J]. 甘肃农业大学学报, 2019, 54(1): 166-174. Doi:10.13432/j.cnki.jgsau.2019.01.023 |
| [12] |
王晓峰, 杨丹, 冯晓明, 等. 黄河中游地区生态恢复对流域水资源的影响[J]. 水土保持通报, 2020, 40(6): 205-212. Doi:10.13961/j.cnki.stbctb.2020.06.030 |
| [13] |
Leng P, Li Z L, Liao Q Y, et al. Determination of all-sky surface soil moisture at fine spatial resolution synergistically using optical/thermal infrared and microwave measurements[J]. Journal of Hydrology, 2019, 579: 124167. Doi:10.1016/j.jhydrol.2019.124167 |
| [14] |
孙亚勇, 马建威, 黄诗峰, 等. 被动微波遥感土壤水分反演研究进展[J]. 中国防汛抗旱, 2021, 31(3): 8-13. Doi:10.16867/j.issn.1673-9264.2021012 |
| [15] |
Li Q L, Shi G S, Shang Guan W, et al. A 1 km daily soil moisture dataset over China using in situ measurement and machine learning[J]. Earth System Science Data, 2022, 14(12): 5267-5286. Doi:10.5194/essd-14-5267-2022 |
| [16] |
Han Q Q, Zeng Y J, Zhang L J, et al. Global long term daily 1km surface soil moisture dataset with physics informed machine learning[J]. Scientific Data, 2023, 10(1): 1-12. Doi:10.1038/s41597-023-02011-7 |
| [17] |
Song P L, Zhang Y Q, Guo J P, et al. A 1km daily surface soil moisture dataset of enhanced coverage under all-weather conditions over China in 2003-2019[J]. Earth System Science Data, 2022, 14(6): 2613-2637. Doi:10.5194/essd-14-2613-2022 |
| [18] |
徐严, 刘程军, 王睿, 等. 黄河流域地级市绿色化空间格局演化及其影响因素研究[J]. 资源开发与市场, 2021, 37(11): 1301-1308. Doi:10.3969/j.issn.1005-8141.2021.11.004 |
| [19] |
Zhang X D, Zhou J, Liang S L, et al. A practical reanalysis data and thermal infrared remote sensing data merging (RTM) method for reconstruction of a 1 km all-weather land surface temperature[J]. Remote Sensing of Environment, 2021, 260: 112437. Doi:10.1016/j.rse.2021.112437 |
| [20] |
Zhou J, Zhang X D, Zhan W F, et al. A thermal sampling depth correction method for land surface temperature estimation from satellite passive microwave observation over barren land[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(8): 4743-4756. Doi:10.1109/TGRS.2017.2698828 |
| [21] |
Zhang Y G, Schaap M G, Zha Y Y. A high-resolution global map of soil hydraulic properties produced by a hierarchical parameterization of a physically based water retention model[J]. Water Resources Research, 2018, 54(12): 9774-9790. Doi:10.1029/2018WR023539 |
| [22] |
Peng S Z, Ding Y X, Liu W Z, et al. 1 km monthly temperature and precipitation dataset for China from 1901 to 2017[J]. Earth System Science Data, 2019, 11(4): 1931-1946. Doi:10.5194/essd-11-1931-2019 |
| [23] |
Leng P, Li Z L, Duan S B, et al. A practical approach for deriving all-weather soil moisture content using combined satellite and meteorological data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2017, 131: 40-51. Doi:10.1016/j.isprsjprs.2017.07.013 |
| [24] |
Moran M S, Clarke T R, Inoue Y, et al. Estimating crop water deficit using the relation between surface-air temperature and spectral vegetation index[J]. Remote Sensing of Environment, 1994, 49(3): 246-263. Doi:10.1016/0034-4257(94)90020-5 |
| [25] |
廖前瑜, 任超, 冷佩, 等. 全天候区域地表蒸散发反演: 以黑河流域为例[J]. 中国农业信息, 2019, 31(1): 35-47. Doi:10.12105/j.issn.1672-0423.20190104 |
| [26] |
乐荣武, 张娜, 王晶杰, 等. 2000—2019年内蒙古草地地上生物量的时空变化特征[J]. 中国科学院大学学报, 2022, 39(1): 21-33. Doi:10.7523/j.ucas.2020.0047 |
| [27] |
王佃来, 刘文萍, 黄心渊. 基于Sen+Mann-Kendall的北京植被变化趋势分析[J]. 计算机工程与应用, 2013, 49(5): 13-17. Doi:10.3778/j.issn.1002-8331.1206-0282 |
| [28] |
莫崇勋, 祝灿, 阮俞理, 等. 广西水资源压力时空特征分析及未来状态预测研究[J]. 中国农村水利水电, 2021(12): 99-104, 110. Doi:10.3969/j.issn.1007-2284.2021.12.016 |
| [29] |
王劲峰, 徐成东. 地理探测器: 原理与展望[J]. 地理学报, 2017, 72(1): 116-134. Doi:10.11821/dlxb201701010 |
| [30] |
周健, 靳诚, 李平星. 基于高速公路流的江苏省交通网络地域分异及影响因素[J]. 中国科学院大学学报, 2022, 39(1): 64-73. Doi:10.7523/j.ucas.2020.0029 |
| [31] |
Leng P, Li Z L, Duan S B, et al. First results of all-weather soil moisture retrieval from an optical/thermal infrared remote-sensing-based operational system in China[J]. International Journal of Remote Sensing, 2019, 40(5/6): 2069-2086. Doi:10.1080/01431161.2018.1468119 |
| [32] |
郭铌, 陈添宇, 雷建勤, 等. 用NOAA卫星可见光和红外资料估算甘肃省东部农田区土壤湿度[J]. 应用气象学报, 1997, 8(2): 212-218. Doi:10.1007/s00376-997-0061-6 |
| [33] |
罗秀陵, 薛勤, 张长虹, 等. 应用NOAA-AVHRR资料监测四川干旱[J]. 气象, 1996, 22(5): 35-38. Doi:10.7519/j.issn.1000-0526.1996.5.008 |
| [34] |
Engman E T, Chauhan N. Status of microwave soil moisture measurements with remote sensing[J]. Remote Sensing of Environment, 1995, 51(1): 189-198. Doi:10.1016/0034-4257(94)00074-w |
| [35] |
山仑, 徐炳成. 新时期黄土高原退耕还林(草)有关问题探讨[J]. 水土保持通报, 2019, 39(6): 295-297. Doi:10.13961/j.cnki.stbctb.2019.06.043 |