格爾木河流域土壤濕度時空變化及其影響因素研究

張 京，金曉媚，張緒財，朱曉倩

(中國地質大學(北京)水資源與環境學院, 北京 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。