張新平,權(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)提供借鑒。
黃河內(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 呈現(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ù)源。
水體與建設(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è)用地。
本研究借助估算模型精度評(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 中,遙感反演的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
圖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
圖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
本研究首先提出了基于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)。