中國東北黑土區地表土壤含水量時空分異特征研究

李雪冬,紅 英,費 龍,劉昀昊,王 浩,薛 澤

(1.長春師范大學地理科學學院,吉林 長春 130032;2.內蒙古師范大學地理科學學院,內蒙古 呼和浩特 010000)

0 前言

土壤含水量(Soil Water Content,SWC)強烈控制著顯熱和潛熱通量之間的能量再分配,是地表與大氣界面的重要狀態參量,直接影響土壤的性質和植物生長狀況,間接控制著植物分布以及生態系統小氣候的變化[1]。土壤水分是氣候、水文、生態、農業等領域衡量土壤干旱程度的重要指標,更是陸面過程模型、水文模型與生態過程模型的一個必要變量[2-3]。土壤含水量的時空分布與變化研究有助于進一步明確植被與土壤的關系,便于理解陸氣間熱量平衡與溫度變化,同時也可為土地覆被空間格局、作物估產、抗旱救災以及水土資源利用等提供決策支持。

SWC時空格局特征與旱情監測通常涉及較大區域,具有持續性和偶然性的特點,同時受外界環境作用導致空間上具有較強的變異性。傳統的土壤水分測量方法[4],如烘干法、土壤時域反射儀(TDR)測量法等,雖然能夠非常精確地獲得SWC信息,但在短時間內獲取大范圍的土壤水分信息時,會受到人力與財力方面的限制。遙感技術具有時效性強、監測范圍廣且成本較低等優點,憑借較高的時空分辨率,可以實現長時間大區域的動態監測。眾多學者基于土壤水分與植被生長狀況和地表溫度的關系理論,尋求基于遙感的土壤水分監測和旱情判別方法,溫度植被干旱指數(Temperature Vegetation Dryness Index,TVDI)是最常用方法之一[5-6],基于MODIS數據的TVDI干旱監測方法在國際上得到了廣泛的應用和認可,也是我國氣象部門干旱監測業務化應用的主要方法之一。PRICE、CARLSON、GOETZ和SANDHOLT等認為NDVI和Ts特征空間呈三角形關系,高中靈等[7]在此基礎上對植被指數進行對數變換,建立了ATVDI方法;HOLZMAN等[8]建立了LST與其他植被指數的關系,發現植被指數對土壤水分反演結果影響較大;劉立文等[9]引入多種植被指數,并對比了各種指數之間的關系,證實了上述結論,同時還發現利用DEM對Ts進行地形校正,反演TVDI具有更高的精度。

合成孔徑雷達(Synthetic Aperture Radar,SAR)憑借穿透性強的特點,使其獲取數據不受云層和光線的影響,并且對土壤表面水分高度敏感,在SWC反演中被廣泛使用。豐富的模型為微波遙感反演土壤水分信息提供了物理基礎和多樣化的工具,但在植被覆蓋條件下,微波信號的組成則變得十分復雜[10],反演植被覆蓋下的SWC重點在于如何有效地分離出植被對微波散射信號的相對貢獻[11]。借助植被的生物、物理特性挖掘植被散射信號的特征,并借助光學遙感數據定量反演植被參數信息是近年來研究熱點[12]。BAGHDADI等[13]利用植被指數來減少植被層的散射貢獻,降低土壤水分反演的不確定性;BAO等[14]使用歸一化水體指數校正水云模型,減少了植被對估算SWC的影響,在植被特征參數表達后向散射信號的能力評價、模型參數的識別以及整體求解方案等方面的研究取得了一定的成果。微波數據雖然在反演SWC時表現出較強的能力,但數據獲取時間和資源上還不及可見光數據。

我國東北黑土區被稱為世界三大黑土區之一,主要分布在松嫩平原以及周圍低臺地區域,具有土層厚、有機質含量高、肥力高等特點,素有“土中之王”的美稱,對于我國糧食安全有著舉足輕重的意義。本研究從土壤水分含量對植被生長過程健康程度與農作物產量影響的相對重要性出發,協同Sentinel-1的合成孔徑雷達數據與Sentinel-2多光譜數據,結合修正的水云模型,實現東北黑土區的SWC定量反演,分析SWC的空間分異特征;為了彌補Sentinel數據量大、獲取時間序列數據困難的不足,借助時間序列MOD11A2和MOD13A3產品計算TVDI,構建黑土區SWC時間序列數據集,分析不同土地利用類型的SWC變化趨勢,進而分析土壤水分的時空分布特征與氣象因子的關系,挖掘土壤水分變化的主導因子。研究結果可為黑土地保護、抗旱監測、作物估產以及研究土地覆被格局特征等提供數據支撐。

1 數據源與方法

1.1 研究區

由于對“黑土”土壤類型的看法不同而導致東北黑土區的界限劃定存在不同方案。本研究選取的“黑土區”為劉寶元等以土壤發生學分類的黑土、黑鈣土、栗鈣土和灰色森林土作為黑土區標志、采用“中心引力集聚法”反復迭代產生的東北典型黑土區[15],該區域面積為33.32萬km2,共涉及138個縣級行政區(表1),分為松嫩黑土亞區、蒙東黑土亞區以及三江黑土亞區。松嫩黑土亞區主要分布在松嫩平原及其四周臺地低丘區,蒙東黑土亞區主要分布在內蒙古東部地區,三江黑土亞區為向東延伸出的部分(圖1)。本研究選取的黑土區土壤結構良好,土層疏松軟綿,呈中性或微酸性,理化性能好,素有“土中之王”的美稱,對我國糧食安全起著關鍵作用,是重要的商品糧基地。該黑土區氣候冬季寒冷,1月平均氣溫-6～-32 ℃,春秋兩季干旱多風,年降雨量一般在400～1 200 mm,夏季雨量充沛,主要集中于6月至9月(全年降水量的70%以上)。地貌分布規律性強且類型多樣,自南向北跨越了暖溫帶、中溫帶、寒溫帶3個不同的自然地帶,從東到西橫穿濕潤區、亞濕潤區、亞干旱區等不同的氣候區,自然植被由東南向西北依次為落葉闊葉林、落葉闊葉及針葉混交林、針葉林和草原,植被生長茂盛,根系發達,在地上與地下積累了大量的有機質。

圖1 研究區位置圖

表1 東北典型黑土區面積統計

1.2 數據源

“哨兵”系列衛星是歐洲哥白尼空間部分(GSC)的專用衛星系列,Sentinel-1搭載了C波段SAR傳感器,涉及波模式、干涉寬幅模式、條帶模式和超寬幅模式等4種數據獲取模式,時間分辨率為12 d,極化方式為VV和VH。利用歐空局提供的SNAP軟件對Sentinel-1A的SAR數據進行斑點濾波、輻射定標和地理編碼,提取相應的微波后向散射系數。Sentinel-2多光譜數據的時間分辨率為5 d,包含13個波段,本研究采用空間分辨率為20 m的植被紅邊波段和短波紅外波段計算NDWI用于反演植被含水量。土地利用類型數據為Sentinel-2衛星圖像構建的10 m分辨率數據。本研究還用到了MODIS產品數據,其中包括地表溫度產品(MODIS11A2,LST)和植被指數產品(MODIS13A3,NDVI),空間分辨率1 km,數據來源于美國國家航空航天局(https://modis.gsfc.nasa.gov),研究區內產品數據的可信度均在95%以上。DEM數據來源于地理空間數據云(http://www.gscloud.cn),圖像預處理包括拼接、裁剪以及壞道處理等,空間分辨率重采樣至1 km,將高程與坡度利用等間隔進行分割,便于后期分析。氣象數據(月平均氣溫和降水數據)來源于中國氣象局氣象數據中心(http://data.cma.cn),氣溫數據對海拔變化較為敏感,插值考慮到了海拔因素的影響,采用ANUSPLIN氣象插值獲得,這樣在很大程度上消除山地區域海拔等因素對溫度的影響,因此空間插值的準確率得到了極大提高。降水采用克里金插值,插值空間分辨率與MODIS產品數據一致。

1.2.2 地面實測數據

地面采樣時間為2020年8月1～5日(圖1),由于Sentinel-1A衛星8月2日從研究區過境,盡量保障了地面實測與衛星過境時間的同步性,樣點主要分布在吉林省和黑龍江省境內。在研究區內共設置29個樣方,每個樣方在2 m×2 m的正方形范圍內,選擇對角線與中心5個位置進行取樣,使用GPS記錄坐標,每次取樣選擇土壤表面0～5 cm左右的土壤,充分混合后取雜質較少的土壤放入鋁盒中稱重,之后在105 ℃的烘干箱中進行充分烘干,待土壤恒重后進行稱重記錄,測得土壤含水量[16]。

1.3 研究方法

1.3.1 土壤含水量反演模型

1.3.1.1 TVDI計算方法

研究表明植被指數與地表溫度構建的特征空間呈現三角形或者梯形分布。SANDHOLT等利用簡化的LST-NDVI特征空間提出了TVDI[17],其中TVDI的表達式為:

在地質遺跡調查中發現的冰臼大多屬于側洞冰臼，就是冰臼的開口朝向旁側，與山頂部位開口向上的冰臼相區別。朗鄉花崗巖石林地質公園是側洞冰臼分布最多的地區。

(1)

式中,T代表TVDI,L代表的是地表溫度(℃);Lmax代表的是給定植被指數下的最大地表溫度(℃);Lmin為給定植被指數下的地表溫度的最小值(℃)。Lmax和Lmin的表達式分別為:

Lmax=a+b×ND,

(2)

Lmin=c+d×ND.

(3)

式(1)為干、濕邊的擬合方程,T數值在0～1之間,其中T值越大,則代表土壤含水量越低,旱情越嚴重,在0.4之內為地表濕潤狀態,無干旱現象,當T在0.4～0.6之間為輕旱,而T大于0.6為較干旱現象,ND表示植被指數NDVI。

1.3.1.2 修正水云模型

1978年ATTEMA和ULABY[18]在估算農作物覆蓋下土壤水分時提出了水云模型(Water-Cloud model),該模型假設植被覆蓋層為均值的散射體,不考慮植被與地表的多次散射,因此將獲取的散射信號看作是農作物的直接反射以及農作物雙次衰減后的地表散射之和,因此微波后向散射可以表示為公式(4):

σC=σV+γ2σS,

(4)

式中,σC為植被覆蓋地表獲取的總的微波后向散射系數,σV與σS分別為植被層的后向散射系數和土壤層直接后向散射系數,γ2為微波穿透植被層的雙層衰減因子。它們的表達式分別為:

σV=Amvcosθ(1-γ2),

(5)

γ2=exp(-2Bmvsecθ),

(6)

式中,mv(kg/m2)為植被含水量,θ為入射角,γ2為植被觀測透過率。A和B為模型的經驗參數,水云模型參數采用BINDLISHI等[19]在不同土地覆被方式中的參數結果。將式(5)、式(6)代入式(4),并變形可獲得土壤層直接后向散射系數為:

(7)

植被含水量是水云模型的重要輸入參數,在植被覆蓋區域的土壤水分反演中具有重要的作用,研究表明植被含水量與光譜指數之間具有良好的擬合關系,本研究采用經驗模型法對不同的植被含水量進行反演,其中模型參數借助JACKSON等[20]的研究結果。

1.3.1.3 模型評價

本研究反演的土壤水分數據用野外實測的數據進行校驗,反演結果精度用決定系數(R2)和均方根誤差(Re)進行評價,決定系數(R2)用于表征實測值和模擬值之間的擬合優度,值越接近于1,代表擬合性越好,Re反映了模擬值與實測值間的偏離程度,值越小模型的反演精度越高。

(8)

式中,Qi為實測值,Pi為模型反演值,n為驗證樣本數量。

1.3.2 分析方法

1.3.2.1 相關分析

本研究中的相關系數采用的是Pearson相關系數,主要用于分析TVDI與溫度和降水間的相關關系,相關系數可用來表示兩個變量之間的緊密程度,公式如下:

(9)

式中,x表示TVDI,y表示溫度或降水,n為樣本數量。rxy為所求變量相關系數,絕對值越大代表相關程度越大。采用t檢驗判斷rxy是否顯著。

1.3.2.2 趨勢分析

本文應用一元線性回歸法,分析黑土區2000—2020年SWC變化趨勢,其變化斜率采用最小二乘法計算,公式為:

(10)

式中,n為研究時段年數;Ti為某像元第i年的TVDI均值;θs為該像元TVDI年際變化斜率,θs大于零時TVDI為增大趨勢,則旱情加重,反之表示旱情減輕。

2 結果分析

2.1 模型適用性評價

本研究分別用到了基于Sentinel微波數據的修正水云模型和MODIS產品數據的TVDI兩種方式反演研究區土壤水分含量,并且對比了Sentinel微波數據VH/VV不同極化方式下的后向散射系數,反演結果用實測的29個采樣點對其進行驗證,結果如圖2所示。圖2(a)、圖2(b)分別為VH和VV極化方式下,總后向散射系數和考慮植被后的后向散射系數對比;圖3為MODIS的TVDI反演結果精度。通過圖2(a)、圖2(b)可以發現,實測的數據與借助Sentinel反演的數據均在擬合值附近,說明Sentinel在反演土壤水分時具有一定的潛力,分別分析圖2(a)、圖2(b)可以發現,在考慮到植被對后向散射系數的影響后,修正的水云模型在VV/VH不同極化方式下的地表后向散射反演結果精度均較高,肯定了植被與地表的交叉散射在總后向散射中的貢獻,在考慮植被影響后的水云模型均較原始的總后向散射系數值反演的結果有明顯提升。對比圖2(a)與圖2(b)可以發現,無論是否考慮到植被影響的條件,VV極化條件下土壤后向散射系數的反演精度均比VH極化條件下得到的土壤后向散射系數的結果要好,其中VV在考慮植被影響的前提下,反演精度最高,R2和Re分別為0.75和0.015。基于Sentinel反演的土壤水分數據雖然可以得到較為理想的效果,但該數據的應用受數據獲取時間的限制,因此本研究采用MODIS計算的TVDI,進行時間序列土壤水分數據反演。由圖3分析可知,基于MODIS的TVDI反演出的土壤水分數據與實測數據也具有較強的相關關系,其中R2和Re分別為0.71和0.053。

圖2 基于Sentinel不同極化方式反演土壤水分(SWC)精度驗證

圖3 TVDI反演土壤水分(SWC)精度驗證

2.2 黑土區地表水分的空間變化特征

為了分析研究區內土壤水分的空間分異情況,本研究基于MODIS產品2000—2020年生長季TVDI的均值以及Sentinel的VV極化方式下基于修正水云模型反演的2020年土壤含水量分別生成了黑土區土壤水分的空間分布圖(圖4)。由圖4可知基于Sentinel的反演結果與MODIS的TVDI反演結果具有相似的空間分布特征,黑土區的土壤水分情況呈現較強的空間異質性。從空間特征而言,內蒙古黑土亞區的土壤水分含量相對較低,干旱現象較為明顯,在松嫩平原區腹地干旱現象較為突出,低值區域主要分布在該區域的耕地區域,在松嫩平原的周邊區域,土壤濕度相對較高。整個黑土區的土壤水分含量相對較低,呈現干旱現象。

圖4 黑土區土壤水分空間分布特征

2.3 黑土區土壤含水量的時間變化特征

為了討論黑土區多年來土壤水分的變化情況以及時間序列特征,本研究分別獲取黑土區2000—2020年各土地利用類型的土壤水分變化趨勢,如圖5所示,(a)表示的是整個黑土區,(b)～(f)分別代表的土地利用類型依次為水田、旱地、有林地、灌木林地、草地。由圖5(a)可知黑土區的TVDI值隨著年份的變化呈現下降趨勢并且具有較強的顯著性特征,在各個土地利用類型中TVDI值均隨著年份均呈現下降趨勢,其中旱地的TVDI值最高為0.84,灌木林地和草地的最大TVDI也均大于0.8,而水田的TVDI值雖然最高值小于0.75,但是在變化趨勢上水田的下降趨勢卻最為明顯,其次為草地和灌木林地,有林地相對變化較弱。不同植被類型區TVDI的平均值由大到小排序為旱地、灌木林地、草地、有林地、水田,這不僅反應了不同土地利用類型TVDI的差異,也強調了空間尺度的重要性。

圖5 不同土地利用類型土壤水分變化情況

由于年均的TVDI值并不能全面反映出研究區內的土壤水分的時空演變特征,因此本研究借助時間序列趨勢分析方法,獲得了研究區內2000—2020年的逐象元的土壤水分變化趨勢的空間分布,得到的趨勢分布和P值(顯著性檢驗)如圖6所示。由圖可知TVDI值呈現升高趨勢的區域面積較小,而呈現下降趨勢的區域則相對較大,研究區內的部分區域趨向干旱,其中變化最為明顯的為松嫩平原的南部以及內蒙古亞黑土區的北部,而三江平原亞區則呈現波動下降的趨勢。

圖6 2000—2020年黑土區土壤水分變化趨勢圖

2.4 影響因素

為了探討氣候因素對土壤水分的影響,本研究將土壤水分數據與溫度和降雨數據進行了相關分析(圖7),結果表明,研究區內的TVDI 與降水量呈負相關,有50.2%的地區TVDI與降水量的負相關關系通過顯著性檢驗,尤其是在內蒙古亞區,土壤水分與降雨的相關性顯著特征較為明顯。研究期間整體黑土區的土壤水分情況與氣溫的相關性特征并不顯著,其中有80.2%的區域沒有通過顯著性檢驗,說明隨著氣溫的升高,TVDI 的變化并不明顯。就整個區域而言,在松嫩平原尤其是南部區域,土地利用類型相對復雜,并且耕地相對較多,人為干擾較為嚴重,從而導致干旱情況不受單一因素控制,是多種因素共同作用的結果,降水量只在一定程度上影響這些地區的干旱情況。

圖7 TVDI變化影響因素分析

3 結論

本文從黑土區退化嚴重以及該區域對我國糧食安全的重要性出發,基于Sentinel-1微波數據,結合Sentinel-2數據,借助水云模型實現了研究區內的土壤水分反演,基于實測數據驗證了Sentinel-1不同極化方式的反演精度,探討了Sentinel在土壤水分反演過程中的潛力。因受Sentinel數據獲取時間限制,本研究進一步借助MOD11A2 LST、MOD13A2 NDVI數據集,計算了TVDI干旱監測指數,并分析了黑土區的時空分布格局,討論了不同土地利用類型長時間序列的土壤水分變化情況。研究的主要結論如下:

第一,Sentinel-1數據和MODIS數據集在監測黑土區土壤水分過程中均具有一定的潛力。在Sentinel-1微波反演土壤水分過程中,植被對后向散射系數影響較大,在去除植被影響后,Sentinel的VV/VH極化方式反演精度均有所提升,且VV極化方式表現出更大的潛力。對比Sentinel-1數據和MODIS數據可以發現,Sentinel反演的土壤水分值要低于MODIS的TVDI反演結果,主要原因可能是MODIS數據并不能有效去除植被的影響。

第二,黑土區內植被生長季的TVDI平均值為0.6,具有較強的空間異質性,其中干旱現象主要分布在松嫩平原腹地的旱地區域。分析2000—2020年的黑土區土壤水分含量的時間變化趨勢可以發現,研究區內的干旱情況較為顯著,但干旱有緩解的趨勢,其中旱地、草地區域旱情較為明顯。

第三,黑土區內的土壤水分情況受到氣溫和降雨的共同作用,研究區內TVDI與氣溫呈現正相關與降雨呈現負相關,分析研究區內的逐象元的控制因子可知,降雨相對于氣溫對研究區內的干旱影響更為直接,其中受降雨影響的區域較大,而受溫度影響的區域較小。整體而言在松嫩平原中部受降雨的影響較為嚴重,而在內蒙古黑土亞區受溫度的影響較為嚴重。

黑土SWC的時空變化特征研究對于黑土地保護、資源合理利用、農業的可持續發展、農業災害的防治等具重要意義。借助微波的土壤水分反演可以提供較高的反演精度,但反演的過程復雜,同時受到地表粗糙度、植被覆蓋度以及土壤質地等因素影響嚴重。微波數據的獲取受時間的限制難以實現長時間序列的推演,因此如何實現快速、便捷的時間序列土壤含水量反演依舊是未來研究的熱點。相對于微波數據基于MODIS的TVDI反演方法更加簡單,但是基于MODIS反演的數據空間分辨率較低。本研究反演結果的驗證主要是借助實測數據,就研究區的面積而言,驗證點的數量略顯不足,并且驗證點主要分布在黑龍江和吉林省境內,因此在驗證的過程中存在一定的局限性。另外,在黑土區的土壤含水量反演結果與氣象因素的相關分析中,氣象因子的插值過程考慮到了地形的影響,但并未單獨分析地形因素對土壤含水量的作用。