基于Landsat 時序影像和非線性邊界監(jiān)測土壤干旱研究
——以1986—2020 年黃河內蒙古段生長季為例

張新平，權 全

(1.西安理工大學 藝術與設計學院,陜西 西安 710054；2.西安理工大學 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室,陜西 西安 710048)

土壤水分（soil moisture,SM）是干旱脅迫生態(tài)系統(tǒng)生產的主要驅動因素[1]。目前，基于遙感反演SM 的方法主要有：基于反射率的光學方法（ORM）、被動微波（MWP）、主動微波（MWA）和協同方法（光學+熱紅外，MWP+MWA，微波+光學）。實際工作中，通常依據SM 與土壤溫度植被干旱指數（TVDI）間的關系，通過遙感反演表層土壤的干旱程度來間接地表征土壤的相對含水量[2-3]。目前，監(jiān)測土壤水分與旱情的遙感指數主要有：水分虧缺指數（WDI）、溫度植被干旱指數（TVDI）、土壤水分虧缺指數（SWDI）、土壤水分指數（SMI）、垂直干旱指數（PDI）、改進的垂直干旱指數（MPDI）、第二種改進的垂直干旱指數（MPDI1）和溫度植被土壤水分干燥指數（TVMDI）。其中，TVDI 僅需通過多光譜衛(wèi)星影像定量計算得出的地表溫度（LST）和植被指數（NDVI），因沒有其他限制條件，成為應用最普遍的干旱指數。關于長時間、廣空間高精度遙感反演TVDI 的研究，主要集中在多源數據融合[4-5]、反演模型改進[6]和下墊面分類考慮[7]。三者同時兼顧的研究較少，僅見于基于Landsat 和Sentinel 在單景影像空間尺度上和LST-NDVI 和STR-NDVI 特征空間中開展的土壤水分與干燥度的遙感反演研究[4]。已有研究認為，土壤含水量等值線與模擬干邊呈非線性相交，特別是高NDVI 地區(qū)的蒸發(fā)修正量隨NDVI 增大而增長更迅速，而低NDVI 則增長緩慢，所以遙感反演TVDI 的過程中需要進行非線性干濕邊界的擬合修訂[8]。前期研究中發(fā)現在地形復雜、地表覆蓋物異質性高的區(qū)域，梯形或三角模型的干濕邊界通常表現為曲線，且非線性邊界的預測精度高于線性邊界[9]。基于Landsat 衛(wèi)星的多景影像拼接后，在不同的土地利用/覆被類型下，通過非線性干濕邊界遙感反演土壤水分與干燥度的研究，尚未見報道。鑒此，本研究提出以下2 個假設并進行證實：（1）基于非線性干濕邊界反演的TVDI 數據是否可靠？（2）長時間、廣空間的TVDI 遙感反演是否可以采用如下流程：同年度植被生長季內不同DOY（day of year）同源影像的LST、NDVI 各自像元水平上的均值合成→邊緣非數值（NaN）區(qū)域裁剪、多景空間拼接→生成各個土地利用/覆被類型的LST-NDVI 特征空間的干濕邊界→反演出各地類的TVDI→合并成整幅TVDI。本研究以黃河內蒙古段為例，以1986—2020 年的步長為4～5 a 的8 期Landsat-5/8 時序影像為數據源，開展了TVDI 估算、精度評價、空間分異等方面的研究，以期為干旱半干旱地區(qū)的土壤干旱監(jiān)測與預報提供借鑒。

1 研究區(qū)概況與數據來源

黃河內蒙古河段地處黃河流域最北端，起自寧夏的石嘴山，止于內蒙古伊克昭盟準格爾旗的馬柵鄉(xiāng)，全長約820 km，為典型的水沙異源、徑流量季節(jié)性波動河段[10]，其中巴彥高勒至頭道拐約520 km 河道為強沖積性河道（圖1）。行政區(qū)劃上涵蓋了呼和浩特、包頭、鄂爾多斯、巴彥淖爾和烏海5 座城市，占地面積約198 671 km2，年均降雨量介于80～500 mm，為干旱半干旱區(qū)域。從區(qū)位經濟體量上看，該區(qū)域是黃河“幾”字灣都市圈的核心地段，區(qū)內有中國“三大一首”自流引水灌溉的河套灌區(qū)，鹽堿化嚴重。據統(tǒng)計，該區(qū)域2019 年底總人口達1 038.12 萬人，建成區(qū)面積為597.79 km2，市區(qū)平均人口密度5 310 人/km2，人均年收入為5 438.3 元，水資源總量為10 198 760 t，相對匱乏。王國慶等研究發(fā)現黃河流域未來30 年水資源量將減少，全流域經濟社會的可持續(xù)發(fā)展將受到威脅[11]。因此，亟待開展高分辨率長時間序列的地表干旱監(jiān)測研究，尋找氣候適應對策，指導該區(qū)域水資源的精準管控和永續(xù)利用。

圖1 研究區(qū)地理信息和所用Landsat 影像行列號Fig.1 Geographical information of study area and the path and row numbers of Landsat images

研究數據來源于美國NASA Landsat 5/8多光譜與熱紅外影像（https://glovis.usgs.gov）。影像分帶信息見圖1，時間跨度為35 a（1986—2020 年），步長為4～5 a，8 期Landsat 系列衛(wèi)星數據，共816 景無云影像（TM 612 景、OLI&TIRS 204 景），在影像云量覆蓋度較高的情況下，用相鄰年份的年積累日（DOY）接近的影像替代，NDVI 與LST 詳細的反演過程見“2.1 TVDI 計算方法”部分。

與研究期一致的8 期（1980s,1990,1995,2000,2005,2010,2015 和2020 年）30 m 土地利用柵格數據集來自中國資源環(huán)境科學與數據中心（https://www.resdc.cn），1980—2020 年每月0.500°×0.625°土壤濕度數據MERRA-2(M2T1NXLND 5.12.4)，獲取自：https://disc.gsfc.nasa.gov/datasets/M2TMNXLND_5.12.4/summary?keywords=https:%2F%2Fdisc.gsfc.nasa.gov；該表層土壤濕度數據為netCDF 格式，其中soil moisture L4，GWETTOP 變量用于本研究TVDI 估算結果的驗證。采用ArcGIS10.2，ENVI5.3，OriginPro2015 處理數據。

2 研究方法與思路

2.1 TVDI 計算方法

圖2 呈現了依據Sandholt 等[12]提出的植被指數-地表溫度特征空間反演TVDI 的光學原理，數學計算過程見式(1)，在三角形UVW中，U點為干燥裸土，V點為濕潤裸土，W點為濕潤密閉冠層；m1表示（t-tmin），代表某一像元與相同植被覆蓋情況下最濕像元的溫度差；m2表示（tmax-tmin），代表在一定植被覆蓋條件下最大溫度差。特征空間內的斜線可看作θ相同的等值線，θ自下而上升高，斜率絕對值大的相對于斜率絕對值小的較為干旱，因此，在一定面積區(qū)域擬合出其特征空間的干邊（LSTmax）、濕邊（LSTmin），即可得到每個像元的干旱指數。

圖2 被指數與地表溫度特征空間和非線性干濕邊界遙感反演TVDI 原理示意Fig.2 Vegetation index and surface temperature space and schematic diagram of TVDI remote sensing version based on the nonlinear dry and wet edges

式中：θ為溫度植被干旱指數（TVDI）；t為地表溫度；tmin、tmax分別為給定NDVI 值下的最小地表溫度、最大地表溫度。本研究中溫度變量的單位均為℃。研究區(qū)海拔較高，平均1 309 m，地表起伏明顯，氣溫差異大，地表溫度受高程影響顯著，需用式（2）對LST 數據進行高程訂正[3,13]。

式中：tdem為修正后地表溫度；t0為遙感反演的地表溫度；訂正系數α=0.006 ℃/m；h為30 m DEM。

本研究參照覃藝等[3]提出的自變量等間隔區(qū)間的LSTmin-NDVI 與LSTmax-NDVI 的矩陣散點數據獲取方法，在ArcGIS 10.2 獲取點對數據，然后在Origin Pro 2015 中繪制成2D 散點圖，并根據擬合模型的決定系數（R2）顯著性和實際符合情況，在3～9 次冪多項式函數中選擇最佳的擬合方程，分別獲取LSTmin與LSTmax的關于NDVI 的非線性方程。

地表溫度通過Landsat 衛(wèi)星的熱紅外波段反演得到，其中TM-6 采用單通道算法[14]；TIRS-10/11 采用廣義單通道[14]與劈裂窗協方差-方差比率（SWCVR）[15]相結合的算法[16]；依據式（3）和（4）獲得地表溫度。

式中：γ與δ為中間反演參數[12,17]；ε為地表比輻射率；ψ1、ψ2和ψ3為大氣參數；d為Landsat 衛(wèi)星的熱紅外波段的數字量化值；G和O為增益與偏置值（頭文件中獲得）。

NDVI 是植被覆蓋的一種表征，在生態(tài)和環(huán)境領域的研究中得到了廣泛應用，其計算式為：

式中：v為NDVI；ρ3、ρ4分別為Landsat-5/8 的紅、近紅外波段的地表反射率。

依據均值合成原理，取各研究期的6—9 月份3 景LST 與NDVI 進行像元水平上的均值計算，將各自的均值合成結果作為各研究期的TVDI 估算的數據源。

2.2 水體與建設用地的剔除

水體與建設用地是地表土壤水分含量接近1 與0 的兩種特殊地類，由于反演TVDI 指數的梯形或三角模型對濕度過度飽和與完全干燥的像元均比較敏感，所以待評估區(qū)域中不應含有滯留的水域和連片的不透水面[18]。借鑒Xu[19]提出的修正的歸一化水體指數（MNDWI，簡記為β）與Feyisa 等[20]提出的自動的水體提取指數（AWEI，簡記為η），同時依據建筑指數（IBI，簡記為λ）和土壤指數（SI，簡記為ξ）的均值合成的歸一化差異建筑與土壤指數（NDBSI，簡記為μ）[21]：

式中：β為MNDWI；η為AWEI；λ為IBI；ξ為SI；μ為NDBSI；ρ1、ρ2、ρ3、ρ4、ρ5和ρ6分別為Landsat-5/8 的藍、綠、紅、近紅外、短波紅外1 波段和短波紅外2 波段的地表反射率。計算AWEI 和NDBSI，以同期高分辨率Google Earth 和土地利用數據設定為閾值，剝除各個時期的水域與建設用地。

2.3 TVDI 估算結果驗證

本研究借助估算模型精度評價指標[9,22]，即平均誤差（ME，簡記為Λ）、平均相對誤差（MRE，簡記為Δ）和均方根誤差（RMSE，簡記為Ω）來評價TVDI 的遙感反演結果的準確性。Λ=0，表示無偏差；Λ＞0，表示高估；Λ＜0，表示低估。Δ與Ω是對模型估算精度的度量，理論上，其值越低越好。

3 結果與分析

3.1 非線性干濕邊界擬合結果

圖3 中，遙感反演的8 期4 類土地利用/覆被類型下的研究區(qū)LST-NDVI 特征空間的非線性干、濕邊方程結果如圖所示?？梢?，干濕邊界不是完美的線性，而是高次多項式曲線，包絡呈不規(guī)則的近似梯形形狀。黃河內蒙古段植物生長季（6—9 月）的LST-NDVI 特征空間的干燥邊界（tmax）與濕潤邊界（tmin），除個別呈9 次多項式特例外，都呈7 次多項式曲線。干燥邊界的R2顯著高于濕潤邊界，且兩者均達到0.05 水平上的顯著，這表明擬合模型是合理的。

圖3 1986—2020 年黃河內蒙古段4 類土地利用下LST-NDVI 特征空間干燥邊界與濕潤邊界非線性擬合結果Fig.3 Non-linear fitting results of dry edge and wet edge of LST-NDVI space in Inner Mongolia section of the Yellow River under four land use types from 1986 to 2020

3.2 TVDI 時空變化

圖4 呈現了8 期6—9 月平均TVDI 的空間分布格局。從空間分布規(guī)律看，依據TVDI 土壤濕度分級標準：0＜θ≤0.2，極濕潤；0.2＜θ≤0.4，濕潤；0.4＜θ≤0.6，正常；0.6＜θ≤0.8，干旱；0.8＜θ≤1.0，極干旱。自1986 年以來，干旱等級經歷了如下變化：干旱主導→正常、干旱共存→正常兼有干旱→干旱兼有正?！珊?、極干旱平分→干旱主導、正常與濕潤鑲嵌。2015 年土壤干旱程度達到峰值，極干旱范圍由研究區(qū)西北部的連片的未利用地（烏拉特后旗）向整個研究區(qū)擴散。2000 年出現了濕潤半濕潤帶，主要分布在研究區(qū)東北部武川縣大青山林區(qū)和黃河北岸河套灌區(qū)一帶。由圖4（i）可知，黃河內蒙古段土地利用/覆被以草地和未利用地為主體，近35 年來，建設農用地、林地面積呈穩(wěn)步增加變化，草地面積呈先增后減變化，水體、農田和未利用地面積呈輕微波動變化?？梢姡蠢玫兀ㄉ车?、鹽堿地、沼澤地、裸土地、裸巖等）是黃河內蒙古段土地干旱化的主要地類，需要加強監(jiān)測和修復治理。另外，因受地形及資源環(huán)境的限制，農田的干旱化傾向，加重了對農業(yè)經濟的影響，農業(yè)用水需求應考慮統(tǒng)籌調配。

圖4 1986—2020 年黃河內蒙古段TVDI 空間分布和土地利用/覆被面積變化Fig.4 Spatial distribution of TVDI and changes of land use/cove ratio in Inner Mongolia section of the Yellow River

3.3 TVDI 估算結果驗證

圖5（a）呈現了8 期54 對遙感反演TVDI 與數據集（M2TUNXLND）的干燥度（1-GWETTOP）2D 散點圖對比結果。本研究反演的TVD I 與（1-GWETTOP）二維散點圖相對均衡地分布在1∶1 參考線（藍線）兩側，整體表明TVDI 反演方法可行、結果準確。從圖5（b）中模型估算精度評價指數折線圖可知，4 類土地利用/覆被下的TVDI 在近35 年間呈先增后減變化，變化的轉折時間出現在2015 年，TVDI 值大小順序為：未利用地＞農田＞草地＞林地，林地和草地在正常與干旱之間變化，未利用地與農田在干旱和極干旱之間變化。TVDI 標準偏差呈增減波動變化，表層土壤旱情的空間異質性也呈現出高低交錯變化。8 個研究期平均誤差分別為：0.030 6、0.004 3、0.060 9、-0.002 6、0.078 1、-0.003 7、-0.002 1 和-0.049 4，平均相對誤差（MRE）依次為：0.075 1、0.021 2、0.192 0、0.482 2、0.232 3、0.018 6、0.020 7 和0.094 2，均方根誤差分別為：0.064 1、0.021 5、0.127 3、0.234 4、0.144 3、0.018 2、0.018 3 和0.104 6，總體誤差為2%～8%，相對準確。這證實了本研究引言部分提出的兩個科學假設，即考慮LULC、非線性干濕邊界的多景Landsat 影像反演TVDI 方法合理、結果準確。

圖5 TVDI 遙感反演精度評價結果Fig.5 Evaluation results of TVDI remote sensing inversion accuracy

4 結語

本研究首先提出了基于Landsat 時序數據、土地利用/覆被數據和非線性干濕邊界生成長時間、廣空間的TVDI 的理論框架與技術流程，然后以黃河內蒙古段為案例進行了實證研究，并將遙感反演TVDI 與NASA 發(fā)布的0.500°×0.625°月尺度的地表診斷數據集（M2TUNXLND）衍生的表層土壤干燥度（1-SM）進行比對評估，證明了本研究所提出的兩個科學假設，即（1）“季節(jié)內均值合成→分用地類型擬合高次多項式干濕邊界→合并生成全區(qū)域TVDI”，該方法流程是可行的；（2）基于上述方法遙感反演的TVDI 數值是可靠的。

傳統(tǒng)的TVDI 反演模型的假設前提是土壤干燥度（1-SM）與LST 呈線性關系，且干燥、濕潤邊界均為線性的。本研究則是通過更高次冪多項式（7 次或9 次）研究了三者之間的關系，干濕邊界的64 個擬合曲線決定系數為0.497 6～0.920 9，均達到了p＜0.001 水平上的顯著性要求，并獲得了可靠的預測結果。這表明，考慮LST-NDVI 特征空間的干濕邊界的非線性關系，在一定程度上可以改善區(qū)域性TVDI 的預測精度。

在利用TVDI 或改進TVDI 進行干旱分析時，大多數研究人員將研究區(qū)域作為一個整體考慮，較少考慮到不同的土地利用/覆被的影響。由于林地、草地、農田和未利用地之間的植被高度和蓋度存在較大差異，相應的地表最高和最低溫度也存在較大差異，導致TVDI 計算結果存在較大差異。同種類型土地利用/覆被提供了相對均質的地表環(huán)境，有利于準確表達TVDI 與LST 和NDVI 之間的關系，有效降低了TVDI 遙感反演的環(huán)境差異誤差。本研究所得結論與SHI 等[7]的結果一致。

在后續(xù)的土壤旱情預測預報研究中，需要深入地開展干燥、濕潤邊界的線性與非線擬合效果評估研究，同時兼顧地貌、氣候和土壤質地等因素。此外，應考慮衛(wèi)星遙感數據的分辨率與下墊面的異質性相匹配。為了貫徹黃河流域生態(tài)保護與高質量發(fā)展重大國家戰(zhàn)略，需要以水定人、以水定產，精準核算黃河“幾”字灣區(qū)域的水資源承載力，助力該區(qū)域都市群高質量發(fā)展，實現黃河大保護生態(tài)目標。