張新平，權(quán) 全

(1.西安理工大學(xué) 藝術(shù)與設(shè)計(jì)學(xué)院,陜西 西安 710054；2.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710048)

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

1 研究區(qū)概況與數(shù)據(jù)來(lái)源

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

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

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

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

2 研究方法與思路

2.1 TVDI 計(jì)算方法

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

圖2 被指數(shù)與地表溫度特征空間和非線(xiàn)性干濕邊界遙感反演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

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

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

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

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

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

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

式中：v為NDVI；ρ3、ρ4分別為L(zhǎng)andsat-5/8 的紅、近紅外波段的地表反射率。

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

2.2 水體與建設(shè)用地的剔除

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

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

2.3 TVDI 估算結(jié)果驗(yàn)證

本研究借助估算模型精度評(píng)價(jià)指標(biāo)[9,22]，即平均誤差（ME，簡(jiǎn)記為Λ）、平均相對(duì)誤差（MRE，簡(jiǎn)記為Δ）和均方根誤差（RMSE，簡(jiǎn)記為Ω）來(lái)評(píng)價(jià)TVDI 的遙感反演結(jié)果的準(zhǔn)確性。Λ=0，表示無(wú)偏差；Λ＞0，表示高估；Λ＜0，表示低估。Δ與Ω是對(duì)模型估算精度的度量，理論上，其值越低越好。

3 結(jié)果與分析

3.1 非線(xiàn)性干濕邊界擬合結(jié)果

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

圖3 1986—2020 年黃河內(nèi)蒙古段4 類(lèi)土地利用下LST-NDVI 特征空間干燥邊界與濕潤(rùn)邊界非線(xiàn)性擬合結(jié)果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 時(shí)空變化

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

圖4 1986—2020 年黃河內(nèi)蒙古段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 估算結(jié)果驗(yàn)證

圖5（a）呈現(xiàn)了8 期54 對(duì)遙感反演TVDI 與數(shù)據(jù)集（M2TUNXLND）的干燥度（1-GWETTOP）2D 散點(diǎn)圖對(duì)比結(jié)果。本研究反演的TVD I 與（1-GWETTOP）二維散點(diǎn)圖相對(duì)均衡地分布在1∶1 參考線(xiàn)（藍(lán)線(xiàn)）兩側(cè)，整體表明TVDI 反演方法可行、結(jié)果準(zhǔn)確。從圖5（b）中模型估算精度評(píng)價(jià)指數(shù)折線(xiàn)圖可知，4 類(lèi)土地利用/覆被下的TVDI 在近35 年間呈先增后減變化，變化的轉(zhuǎn)折時(shí)間出現(xiàn)在2015 年，TVDI 值大小順序?yàn)椋何蠢玫兀巨r(nóng)田＞草地＞林地，林地和草地在正常與干旱之間變化，未利用地與農(nóng)田在干旱和極干旱之間變化。TVDI 標(biāo)準(zhǔn)偏差呈增減波動(dòng)變化，表層土壤旱情的空間異質(zhì)性也呈現(xiàn)出高低交錯(cuò)變化。8 個(gè)研究期平均誤差分別為：0.030 6、0.004 3、0.060 9、-0.002 6、0.078 1、-0.003 7、-0.002 1 和-0.049 4，平均相對(duì)誤差（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%，相對(duì)準(zhǔn)確。這證實(shí)了本研究引言部分提出的兩個(gè)科學(xué)假設(shè)，即考慮LULC、非線(xiàn)性干濕邊界的多景Landsat 影像反演TVDI 方法合理、結(jié)果準(zhǔn)確。

圖5 TVDI 遙感反演精度評(píng)價(jià)結(jié)果Fig.5 Evaluation results of TVDI remote sensing inversion accuracy

4 結(jié)語(yǔ)

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

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

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

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