格尔木河流域土壤湿度时空变化及其影响因素研究

张 京，金晓媚，张绪财，朱晓倩

(中国地质大学(北京)水资源与环境学院, 北京 100083)

土壤湿度是气候、植被、地形及土壤因素等自然条件的综合反应，是地下水—地表水—大气水循环系统的核心和纽带，是干旱、半干旱地区生态系统和植被建设的基础，对于维持干旱地区内靠天然降水和地下水支撑的天然、人工植被的稳定性具有重要作用。土壤湿度作为地表特征参数之一，对作物的生长具有决定性的作用，也是进行土地退化、植被覆盖及干旱等生态环境研究的重要指标之一，是地下水资源评价的关键部分。区域尺度甚至全球范围的土壤湿度反演是陆面过程模式研究中必不可少的参数，在改善区域及全球气候、预测区域干湿情况研究中意义重大[1-2]。

土壤湿度的研究方法可分为传统方法和遥感方法两大类。传统方法主要有重量水分法、土壤湿度计法、蒸渗法和电阻法等，但因为样地系数代表范围有限，因而宏观性不够，数据得不到及时更新[3]。20世纪80年代，遥感监测方法飞速发展，遥感监测土壤水分克服了以上的缺陷，其具有快速、实时、长期、动态和大区域监测以及良好的时间、空间分辨率等优点，大大提高了工作效率[4-5]。辛景峰[6]利用1981—1994年的NOAA/AVHRR数据，利用全国102个固定农业观测站20 cm的土壤湿度数据加以验证，研究得出土壤湿度与地表温度/植被指数(Ts/NDVI)的定量关系。Doninck[7]等采用 Aqua 星和 Terra 星采集的交叉时相的MODIS数据改进了表观热惯量模型中地表温度的算法，对非洲南部地区进行了土壤湿度的反演，结果接近NASA提供的基于AMSR-E数据的土壤湿度，可靠性较高。杨树聪[8]设计实验研究表明在植被覆盖率较低(NDVI<0.35)情况下，表观热惯量与土壤含水量之间相关性很高，相关系数大于0.7，而植被覆盖较好(NDVI>0.35)的情况下，表观热惯量不再适用。吴黎[9]等尝试用改进的热惯量模型计算不同植被覆盖下、不同土壤含水量实验区的热惯量值，并与实测土壤含水量进行相关性分析，发现NDVI≤0.35时热惯量方法反演的土壤含水量与实测土壤湿度值匹配度较高。

本文研究基于中等分辨率的MODIS遥感数据，利用表观热惯量(ATI)法反演得到格尔木河流域的土壤湿度时空分布特征，并分析了其影响因素，对当地的生态环境保护及水资源评价具有重要的科学意义，也为该区地下水资源开发利用提供了合理依据。

1 研究区概况

格尔木河流域地处青海省柴达木盆地南缘中部，南倚昆仑山脉，北临达布逊湖，平均海拔2 780 m,面积4 566 km2。该地区年均降水量42.7 mm，降雨集中在6、7、8三月，多年平均蒸发量1 495 mm，年均气温4.7 ℃，相对湿度32%。流域地理位置遥感图见图1。

图1 格尔木河流域地理位置遥感图Fig.1 Remote sensing map of the geographical location of the Golmud River Basin

格尔木地区具有典型的高原大陆性气候特征，干旱、多风少雨，蒸发强烈。流域内格尔木河是柴达木盆地第二大内陆河流，以地下水补给为主，补给量占年径流量的65%以上，其余为冰雪融水和雨水，分别占年径流量的22.76%和12.09%[10]。格尔木河上游为山区，主要为高寒干草原、高寒荒漠和高寒草地，其植物群落的生长发育主要依靠大气降水和地表水滋补，基本不受地下水影响[11]，主要植被类型有沙蒿。河水出山后在山前戈壁带大量入渗补给地下水，至细土带地下水又大量溢出形成多条泉集河，此处土质厚，水质较好，生长有白刺、芦苇、芨芨草等植被；至冲湖积平原区，土壤含盐量增高，地下水质变差，则主要生长红柳、矮芦苇等旱生耐盐植被；湖积平原主要为盐沼地，除河岸两带有少量植被呈带状生长外，其余基本无植被生长。

格尔木河径流年内分配比较均匀，径流年际变化较小，是一条水量变化小而稳定的河流。终端的东达布逊湖则是格尔木河地表、地下水排汇中心。

2 数据选取与研究方法

2.1 数据选取与预处理

(1)MODIS数据。本次反演研究区土壤湿度选用的遥感数据是Aqua卫星采集的中等分辨率MODIS数据，005数据集版本(2017、2018年为006数据集版本)，主要有MYD09A1(地表反射率产品、500 m空间分辨率)和MYD11A2(地表温度产品、1 000 m空间分辨率)，时间分辨率均为8天，时间为2002—2018年每年6—9月，共计488景。数据来源于美国NASA网站(https://ladsweb.modaps.eosdis.nasa.gov/)。

为了研究数据的一致性，将所有数据在MRT软件中进行投影转换，重采样为500 m，地理坐标系为WGS84，投影坐标系为Geographic Lat/Lon，然后在ENVI 5.1软件中进行波段的计算和叠加处理。

(2)CLDAS数据。“CLDAS”全称中国气象局陆面同化系统。“CLDAS土壤体积含水量(m3·m-3)分析产品V2.0”为覆盖亚洲区域(0°～65°N，60°～160°E)，逐小时、垂直分为5层(0～5、0～10、10～40、40～100、100～200 cm)，0.0625°×0.0625°空间分辨率的经纬度网格土壤湿度产品。该数据集研制技术和精度与国际同类产品(如GLDAS、NLDAS产品)相当，在中国区域质量优于国际同类产品，且时空分辨率更高。本次研究采用2017年6—9月CLDAS数据(共计20景)对MODIS数据反演的土壤湿度结果进行验证。数据来源于中国气象数据网(http://data.cma.cn/)。首先需要在NetCDF4Excel中对CLDAS数据进行读取，然后在ArcGIS中进行插值处理。

2.2 研究方法

(1)表观热惯量法

土壤热惯量是度量土壤热惰性大小的物理量。是土壤的热学特性，表征了土壤阻止其自身温度变化能力的大小：

(1)

式中：P——热惯量；

K——土壤热传导率/(J·cm-1·s-1·K-1)；

ρ——土壤密度/(g·cm-3)；

c——土壤比热容/(J·g-1·K-1)。

对于确定的土壤类型，土壤热惯量随土壤水分增加而增大。但由于常规土壤热惯量的计算模型涉及很多物理参数，如表面粗糙度、空气湿度、风速等物理参数获取及计算较为困难复杂，Price根据能量平衡原理提出表观热惯量进而估算土壤含水量[17]。热惯量简化后的形式为：

(2)

式中：ATI——表观热惯量；

A——地表反照率；

Td、Tn——表示一天中最高、最低地表温度。

地表反照率A的计算采用Liang S L提出的用宽波段反照率代替全波段反照率算法，α代表各波段地表反射率[18]，计算公式为：

A=0.16α1+0.291α2+0.243α3+0.116α4+

0.112α5+0.081α7-0.0015 (3)

表观热惯量的计算简化了参数，式中一天中最高地温Td、最低地温Tn及计算地表反照率的各参数容易从遥感数据中获得，计算较为简单，经过MRT预处理之后在ENVI5.1中进行波段计算，比较容易推广使用。

经MODIS数据计算，多年来格尔木河流域地区年平均最高NDVI值为0.097(2016年)，年内最高NDVI值为0.135(2014年8月)，且全区NDVI小于0.35的区域占全区的95.48%，属于低植被干旱半干旱地区，所以本次研究采用表观热惯量法。

(2)Mann-Kendall趋势检验法[19-23]

采用Mann-Kendall方法度量研究区ATI变化趋势的显著性，其趋势检验法过程如下：

对于序列Xt=(x1,x2,…,xn)，∀i,j≤n且i≠j，做如下假设：H0：序列中的数据随机排列，即无显著趋势;H1：序列存在上升或下降单调趋势。定义统计量S：

(4)

当n≥10 时，统计量S近似服从标准正态分布，使用检验统计量Ζ进行趋势检验，Ζ值由式(6)计算:

(6)

(3) Mann-Kendall法突变检测[19-23]

采用Mann-Kendall方法度量研究区ATI变化过程中的突变点，其突变检测过程如下：

给定序列为xn,xn-1,…,x1，Sk表示第i个样本xi>xj(1≤j≤i≤n)的累计数，定义统计量：

(7)

(j=1,2,…,i;k=1,2,…,n)

在时间序列随机独立的假定下，Sk的均值和方差分别为：

E[Sk]=k(k-1)/4,

Var[Sk]=k(k-1)(2k+5)/72 ，1≤k≤n(8)

将Sk标准化：

(9)

UBk=-UFk,i′=n+1-ii,i′=1,2,…,n(10)

3 结果与分析

3.1 土壤湿度时空分布特征

计算表观热惯量ATI的MODIS数据时间分辨率为8天，每月3或4幅。首先在MRT中对MODIS数据进行预处理，接下来在ENVI5.1中进行波段叠加计算，统计分析，利用表观热惯量法计算日幅ATI空间分布图像，并用每月内多幅日值ATI影像均值作为月值ATI的反映。由于格尔木河流域地处柴达木盆地内部，每年6—9月气温较高，其他月份气温较低，不适合用表观热惯量法计算，因此将每年6—9月份的平均ATI作为此地区年平均ATI值对土壤湿度进行反演，得到格尔木河流域表观热惯量分布图(图2)，选取2003年、2008年、2013年、2016年图为例。

图2 格尔木河流域特征年份表观热惯量(ATI)空间分布图Fig.2 Spatial distribution of the apparent thermal inertia (ATI) in characteristic years of the Golmud River Basin

表观热惯量在一定程度上能够反应土壤湿度的分布特征。一般来说，表观热惯量越大，土壤湿度越大。根据图2可知，格尔木河流域大部分为荒漠区，ATI值比较低，在0.02～0.03之间；由于植被主要发育在研究区中南部，植被对土壤湿度有一定的响应作用，所以在植被发育较好的地区ATI值相对较高。主要表现在以格尔木市为中心，西侧的格尔木西农场植被以耕地为主，主要种植枸杞、藜麦等农作物，需水量较大，ATI值比较高，在0.035～0.045之间；在格尔木市东部主要发育芦苇、芨芨草、沙蒿、红柳等灌木种类，ATI值在0.03～0.04之间；研究区南部主要为荒漠，基本无植被发育，ATI值最小，在0.025左右；向北延伸，盐渍化不断加重，地表含盐量不断升高，最北部的察尔汗盐湖，盐场面积逐年扩大，ATI值最高在0.05以上。

2002—2016年间，研究区年平均ATI值总体较小，土壤湿度有限，生态环境脆弱，值在0.020～0.030之间波动，呈波动上升趋势(图3)。

图3 格尔木河流域2002—2016年ATI年际变化趋势Fig.3 Interannual variation trend of ATI in the Golmud River Basin from 2002 to 2016

3.1.1Mann-Kendall趋势检验法

图4 格尔木河流域2002—2016年ATI变化趋势图Fig.4 Variation trend of ATI in the Golmud River Basin from 2002 to 2016

经过修正后，临界值为-1.27和2.65。大于0的范围表观热惯量呈上升趋势，小于0的部分表观热惯量呈下降趋势。由图中可知，小于-1.27属于显著下降，主要分布在中部偏东的一小部分土壤盐渍化较为严重地带。大于2.65属于显著上升，主要集中在研究区北部盐池以及格尔木市东西两侧农田及人工防护林。总体来说，ATI上升范围明显大于下降范围，研究区ATI总体呈上升趋势。

表1 95%显著性检验参数Table 1 Parameters of 95% significance test

3.1.2Mann-Kendall法突变检测

在Excel中实现Mann-Kendall法突变检测，根据UFk和UBk曲线走势，格尔木河流域的ATI在2006年以前呈下降趋势，自2006年以后开始呈上升趋势，且从2011年起这种上升趋势越发显著超过临界水平，表明研究区年平均ATI上升的趋势是十分显著的。根据UFk和UBk曲线交点的位置，研究区年平均ATI在2008年开始突变(图5)。突变检验结果与图3所示结果基本符合，证明研究结果具有一定的可靠性。

图5 格尔木河流域ATI突变检验曲线图Fig.5 Mutation test curve of ATI in the Golmud River Basin

3.2 不同用地类型的表观热惯量分布

对2015年ATI空间分布图按照北—南(AB)和西—东(CD)两向分别进行剖面提取(图1)，分析同一剖面下不同用地类型的表观热惯量分布及关系(图6)。表2通过提取不同用地类型，利用ArcGIS统计分析了单一用地类型的ATI平均值。

图6 格尔木河流域不同用地类型的ATI变化Fig.6 ATI changes in different types of landuse in the Golmud River Basin

用地类型盐池自然植被耕地、林地建筑用地河流、湖泊裸土ATI均值0.0880.0340.0330.0310.0300.028

通过图6、表2分析可得，研究区最北部盐池含水量最高，ATI均值最高，达到0.088；其次是自然植被(0.034)和耕林地(0.033)，均分布于格尔木市东西侧周围地带；河流、湖泊分布较少且稀疏，ATI均值较植物略低，但由于存在混合像元并且近年来河道较干，所以在剖面图中像元较少；中北部和南部荒漠裸土的ATI均值最低，为0.028。图、表所得结果大致一致，证明计算结果具有一定的可靠性。

3.3 CLDAS土壤湿度产品验证

由于受CLDAS数据获取时间的限制，采用2017年6-9月中的7月份MODIS数据和CLDAS数据进行验证。用NetCDF4Excel获取CLDAS土壤含水量数据，在ArcGIS中将ATI图像重采样为对应的CLDAS产品像元分辨率大小0.0625°，以0.01步长进行重分类；对CLDAS数据进行插值分析，最后以表格显示分区统计，得到ATI与CLDAS各层土壤含水量的相关性图(图7)。如图所示，CLDAS土壤湿度数据各层之间分布趋势一致，但10～40 cm、40～100 cm、100～200 cm层的土壤体积含水量比0～10 cm层的高。ATI与CLDAS各层土壤含水量数据拟合都较好，拟合度最高的是10～40 cm层，相关性系数达到0.85，说明MODIS数据反演研究区土壤湿度分布特征是可靠的。

图7 ATI与CLDAS各层数据相关性分析Fig.7 Correlation analysis of each layer of ATI and CLDAS

3.4 野外实测土壤含水量数据验证

为了进一步验证反演结果的可靠性，2018年8月进行研究区野外实地数据测量与采集(图8)，所用仪器为ML3 Theta Probe 便携式土壤湿度探头配合HH2 读数表，仪器基于FDR(Frequency Domain Reflectometry)频域反射原理来测量土壤的表观介电常数，从而得到土壤容积含水量，测量精度±1%。本次野外定点共55个，每个点在开挖剖面的10 cm、20 cm、30 cm层处不同方位分别进行3次测量，并取平均值作为该点该层的土壤含水量值。

图8 土壤湿度实测点位及野外路线图Fig.8 Soil moisture measurement point and field road map

采用2018年8月MODIS遥感数据反演的结果与2018年8月野外实测不同深度平均土壤含水量数据进行相关性分析(图9)，两者相关系数达到0.59，证明遥感反演土壤湿度结果具有一定可靠性。

图9 反演土壤湿度与不同深度平均实测土壤湿度相关性分析Fig.9 Correlation analysis of inversion of soil moisture and average measured soil moisture at different depths

3.5 影响因素分析

3.5.1植被与土壤湿度的相关性

在干旱半干旱地区，植被对土壤湿度有一定的指示作用。植被可以涵养水源，改善干旱地区的生态环境；而土壤相对比较湿润的地区又可以促进植物的生长。图10(a)是利用ENVI5.1计算得到的研究区2015年的植被指数(NDVI)空间分布图，结合图10(b)2015年的表观热惯量(ATI)空间分布图，不计北部的盐池，两者空间分布具有一致性：NDVI值高的地区，ATI相对也高，植被长势较好。

图10 格尔木河流域2015年NDVI(a)与ATI(b)空间分布图Fig.10 Spatial distribution of NDVI (a) and ATI (b) in 2015 in the Golmud River Basin

为了进一步分析ATI与NDVI的定量关系，利用ArcGIS中的区域统计模块，建立了ATI与NDVI(不包括水体)的相关性分析图(图11)。结果表明，2015年ATI与NDVI(不包括水体)两者之间的相关性系数达到0.92，呈正相关，表明植被与土壤湿度有良好的相关性。

图11 格尔木河流域2015年ATI与NDVI(不包括水体)相关性分析Fig.11 Correlation analysis between ATI and NDVI (excluding water bodies) in 2015 in the Golmud River Basin

3.5.2蒸发与土壤湿度的相关性

蒸发是水量平衡和地表能量平衡的主要组成部分。分析蒸发与土壤湿度的相关性对于当地的水资源合理开发利用具有重要意义。图12(a)是利用GLDAS数据通过SEBS模型[24]计算得到的格尔木河流域2016年的年蒸散发(ET)空间分布图，结合图12(b)2016年的表观热惯量(ATI)空间分布图，发现两者空间分布具有一致性，ET值高的地区，ATI相对也高，两者具有良好的响应关系。

图12 格尔木河流域2016年ET(a)与ATI(b)空间分布图Fig.12 Spatial distribution of ET(a) and ATI (b) in 2016 in the Golmud River Basin

进一步分析ATI与ET的定量关系。从图13(a)可知2002—2016年ET与ATI变化趋势大致一致，呈上升趋势；从图13(b)可知，两者的相关系数达到0.79，蒸散发和表观热惯量之间呈正相关。

图13 格尔木河流域2002—2016年6—9月间ATI与ET相关性分析Fig.13 Correlation analysis of ATI and ET between June to Sepetember from 2002 to 2016 in the Golmud River Basin

4 结论

(1)2002—2016年间，研究区的ATI总体呈上升趋势，并在2008年发生突变。

(2)研究区北部盐池平均ATI值最高，为0.088；中部及南部大面积裸土平均ATI值最低，为0.028；自然植被与耕、林地平均ATI值分别为0.034和0.033。

(3)利用CLDAS土壤湿度产品和野外实测土壤含水量数据对MODIS数据反演结果进行验证，相关系数分别为0.85和0.59，验证结果均证明了反演结果的可靠性。

(4)植被和蒸散发对土壤湿度都有积极的响应作用，呈正相关，相关性系数分别为0.92、0.79。