基于理論干濕邊與改進(jìn)TVDI的麥田土壤水分估算研究

蔡慶空 陶亮亮 蔣瑞波 蔣金豹

(1.河南工程學(xué)院土木工程學(xué)院， 鄭州 451191； 2.南京信息工程大學(xué)地理科學(xué)學(xué)院， 南京 210044；3.中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院， 北京 100083)

0 引言

土壤水分是地球系統(tǒng)的重要組成部分，影響著水圈、生物圈、大氣圈的水熱平衡，是農(nóng)業(yè)、水文、生態(tài)、氣候等領(lǐng)域的關(guān)鍵參數(shù)，對旱情監(jiān)測及農(nóng)作物生長發(fā)育有著十分重要的作用。

傳統(tǒng)基于點(diǎn)的測量方法無法滿足目前用于制作大面積、高精度土壤水分產(chǎn)品的應(yīng)用需求，而遙感技術(shù)具有快速、大面積以及實(shí)時等優(yōu)點(diǎn)，為估算大面積土壤水分時空信息、實(shí)時監(jiān)測土壤水分動態(tài)變化提供了便利。光學(xué)遙感作為目前最常用、最成熟的土壤水分遙感監(jiān)測方法而被廣泛使用，其波段豐富、成像質(zhì)量好、空間拓?fù)潢P(guān)系清晰。目前，利用光學(xué)遙感反演地表土壤水分含量最常用的監(jiān)測方法包括歸一化植被指數(shù) (Normalized difference vegetation index，NDVI)[1]、表觀熱慣量 (Apparent thermal inertia，ATI)[2-3]、溫度植被干旱指數(shù) (Temperature vegetation drought index，TVDI)[4]、作物缺水指數(shù) (Crop water stress index，CWSI)[5-6]以及高光譜遙感等。其中，TVDI方法通過將可見光波段與熱紅外波段結(jié)合構(gòu)建特征空間來估算土壤水分含量，在土壤水分的監(jiān)測中取得了較多的研究成果[7-12]。PATEL等[13]利用TVDI方法對印度地區(qū)的土壤水分進(jìn)行估算，結(jié)果表明，TVDI與土壤水分之間存在非常好的負(fù)相關(guān)關(guān)系，尤其是在植被稀疏的情況下。ZHANG等[14]采用MODIS遙感產(chǎn)品數(shù)據(jù)構(gòu)建地表溫度及NDVI特征空間，對關(guān)中地區(qū)土壤水分空間變化進(jìn)行研究，結(jié)果表明，TVDI與地表0～10 cm和10～20 cm土壤水分均有較高的相關(guān)性，因此TVDI可以反映并指示地表0～20 cm深度的土壤水分狀況，可以用于關(guān)中地區(qū)干旱動態(tài)監(jiān)測。

基于地表溫度(Land surface temperature，LST)和植被指數(shù)(Vegetation index，VI)特征空間的TVDI在選取干濕邊以及用于土壤水分監(jiān)測時存在不足，理論上TVDI特征空間應(yīng)包含3種植被覆蓋類型，即裸土、部分植被覆蓋以及全植被覆蓋類型，但實(shí)際應(yīng)用過程中，研究區(qū)域的覆蓋狀況往往無法完全滿足植被覆蓋條件，得到的干濕邊實(shí)質(zhì)上只是理論特征空間內(nèi)部存在的邊，而非理論干濕邊。本文從地表能量平衡方程出發(fā)，根據(jù)能量交換過程中的極端狀況，獲取理論干濕邊的端點(diǎn)方程，從而得到理論意義的特征空間，進(jìn)而對試驗(yàn)區(qū)域內(nèi)地表土壤水分進(jìn)行反演研究，以期獲得更高精度的土壤含水率估算值。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于關(guān)中平原中部地區(qū)，主要包括楊凌區(qū)、扶風(fēng)縣以及武功縣(北緯34.1°～34.5°，東經(jīng)107.8°～108.3°，圖1)，地勢平坦，氣候溫和，屬于暖溫帶大陸性季風(fēng)氣候，年平均植被蒸發(fā)量為993.2 mm，平均日照時數(shù)2 163.8 h，主要種植作物為小麥和玉米。研究區(qū)內(nèi)以水澆地(即農(nóng)田)為主，部分區(qū)域有城鎮(zhèn)用地、河流、山地等地表覆蓋類型。

圖1 研究區(qū)概況Fig.1 Overview of study area

以陜西楊凌區(qū)作為核心區(qū)，開展多次野外試驗(yàn)，主要獲取地面土壤水分?jǐn)?shù)據(jù)以及植株參數(shù)等，共設(shè)置3個小麥試驗(yàn)站點(diǎn)，每個試驗(yàn)站點(diǎn)根據(jù)農(nóng)田觀測面積設(shè)置不同的固定樣點(diǎn)，同時以標(biāo)桿標(biāo)記并進(jìn)行數(shù)據(jù)記錄。觀測時間為2014年3月29日和4月22日，觀測數(shù)據(jù)主要包括地表土壤水分含量、地表溫度、植株高度、葉面積指數(shù)以及光譜數(shù)據(jù)等。

所使用的氣象科學(xué)數(shù)據(jù)從中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)下載，主要包括中國地面氣候資料日值數(shù)據(jù)集和中國農(nóng)作物生長發(fā)育和農(nóng)田土壤濕度旬值數(shù)據(jù)集等。

1.2 衛(wèi)星數(shù)據(jù)與預(yù)處理

中分辨率成像光譜輻射計(jì)(Moderate resolution imaging spectroradiometer，MODIS)共有36個離散波段，空間分辨率為250、500、1 000 m，掃描寬度為2 330 km。目前MODIS提供多種業(yè)務(wù)化運(yùn)行產(chǎn)品，包括地表參量、海洋參量以及大氣參量等從原始數(shù)據(jù)0級到應(yīng)用模型開發(fā)5級及以上級別產(chǎn)品。本研究主要使用2014年3月29日(DOY088)和4月22日(DOY112)MYD09GA每日地表反射率數(shù)據(jù)以及MYD11A1每日地表溫度數(shù)據(jù)(可從https:∥ladsweb.modaps.eosdis.nasa.gov/下載)。其中MYD09GA地表反射率數(shù)據(jù)共7個波段，空間分辨率為500 m。該產(chǎn)品為經(jīng)過大氣校正的2級地表反射率產(chǎn)品數(shù)據(jù)，可以得到星下點(diǎn)和標(biāo)準(zhǔn)太陽高度角，用于計(jì)算每日反射率及歸一化植被指數(shù)(NDVI)。MYD11A1每日地表溫度數(shù)據(jù)空間分辨率為1 000 m，利用MRT(MODIS reprojection tool)軟件進(jìn)行處理，投影類型為Albers Equal Area，并將影像分辨率從1 000 m重采樣至500 m，用于計(jì)算空氣溫度、理論干濕邊及觀測干濕邊。

1.3 研究方法

1.3.1溫度植被干旱指數(shù)及Ts-Fv特征空間

研究發(fā)現(xiàn)地表溫度與植被指數(shù)(VI)所構(gòu)建的特征空間內(nèi)存在多條土壤水分含量等值線，并且對于每條等值線上的點(diǎn)所對應(yīng)的土壤水分含量均相等，該特征空間將地表覆蓋區(qū)域分為裸土區(qū)、部分植被覆蓋區(qū)以及全植被覆蓋區(qū)，其所構(gòu)建的干濕邊與植被指數(shù)呈現(xiàn)顯著的線性關(guān)系。基于此，SANDHOLT等[4]提出溫度植被干旱指數(shù)(TVDI)，該指數(shù)與地表土壤水分密切相關(guān)，是表征植被受水分脅迫以及地表干濕狀況的重要指標(biāo)，其表達(dá)式為

TVDI=(Ts-Tmin)/(Tmax-Tmin)

(1)

其中

Tmax=a+bVI

(2)

Tmin=c+dVI

(3)

式中Ts——影像相應(yīng)像元的地表溫度

Tmin——Ts-VI特征空間的最低地表溫度，即特征空間的“濕邊”

Tmax——Ts-VI特征空間的最高地表溫度，即特征空間的“干邊”

VI——植被指數(shù)

a、b、c、d——經(jīng)驗(yàn)參數(shù)，可線性擬合得到

目前最常用的植被指數(shù)為歸一化植被指數(shù)，但是在植被高覆蓋區(qū)域，NDVI易于飽和，同時其對于植被冠層背景影響較為敏感，使土壤水分估算不確定，因此本文將NDVI替換為改進(jìn)植被覆蓋度Fv[15]，其計(jì)算公式為

(4)

式中NDVIi——影像任一像元的NDVI

NDVImax——全植被覆蓋對應(yīng)的NDVI

NDVImin——全裸土覆蓋對應(yīng)的NDVI

構(gòu)建的Ts-Fv特征空間(圖2)不僅可以表示植被的生長狀況，同時能夠一定程度上避免TVDI特征空間估算時對于地表覆蓋類型的限制，應(yīng)用范圍更廣泛。

圖2 理論干邊與觀測干邊在梯形特征空間位置圖Fig.2 Schematic of theoretical and observed dry edges in trapezoidal feature space

1.3.2地表能量平衡方程

無論是觀測干濕邊還是理論干濕邊，實(shí)質(zhì)上都受到地表能量平衡方程影響，而求導(dǎo)理論干邊及理論濕邊4個端點(diǎn)的方法均需要從地表能量平衡方程出發(fā)。地表能量平衡方程表示了太陽輻射能量除去被大氣吸收以及被大氣頂層反射回太空以外到達(dá)地表的能量守恒情況[16]。其可以表示為

Rn-G=H+LE

(5)

(6)

αs=0.160α1+0.291α2+0.243α3+0.116α4+
0.112α5+0.081α7-0.001 5

(7)

(8)

(9)

εs=0.261+0.314ε31+0.411ε32

(10)

G=Rn[Γv+(Γs-Γv)(1-fv)]

(11)

(12)

(13)

(14)

d0=0.65h

(15)

z0=0.13h

(16)

(17)

(18)

(19)

(20)

γ=0.665×10-3P

(21)

式中Rn——地表凈輻射通量

G——地表土壤熱通量

H——地表顯熱通量

LE——地表潛熱通量

αs——地表反射率[17]

S0——太陽總輻射，利用6S輻射傳輸模型得到

σ——斯-玻常數(shù)，取5.67×10-8W/(m2·K4)

α1、α2、α3、α4、α5、α7——MODIS09反射率產(chǎn)品的1、2、3、4、5、7波段

εa——無云時大氣有效發(fā)射率[18]

ea、Ta——參考溫度水氣壓和近地表空氣溫度，可由半經(jīng)驗(yàn)方程估算得到[19]

Tsky——天空等效溫度[20]

εs——地表發(fā)射率[21]

ε31、ε32——MODIS產(chǎn)品第31、32波段發(fā)射率

Γv——植被覆蓋下G和Rn比率，取0.05[22]

Γs——裸土覆蓋下G和Rn比率，取0.315[22]

fv——植被覆蓋度

ρ——空氣密度，取1.29 kg/m3

CP——空氣定壓比熱，取1.005 kJ/(kg·K)

ra——空氣動力學(xué)阻抗[23]

z——地表參考高度，取2 m

k——卡曼常數(shù)，取0.41

u——參考高度2 m處的風(fēng)速

d0——零平面位置高度[22]

z0——表面粗糙度[16]

h——植株平均高度

Δ——飽和水氣壓曲線斜率

e0(Ta)——飽和水氣壓[24]

dz——參考高度水汽飽和差

Tdmax、Tdmin——日最高氣溫、日最低氣溫

γ——干濕球常數(shù)P——大氣壓

1.3.3理論干濕邊端點(diǎn)選取及方程推算

傳統(tǒng)選取干濕邊的方法即為選擇植被指數(shù)最大和最小時地表溫度的最大與最小值，構(gòu)建的Ts-NDVI特征空間一般將濕邊處理為與NDVI軸相平行的直線，干邊處理為與NDVI呈線性關(guān)系，但是多數(shù)情況下，干邊的最高溫度往往隨著NDVI的增大呈現(xiàn)先增大后減小的趨勢，并且隨著時間、季節(jié)以及緯度的變化而有所差異。因此利用地表能量平衡方程獲取理論干濕邊的端點(diǎn)顯得尤為重要，理論濕邊是由蒸騰和蒸發(fā)最大時裸土及植被覆蓋下最低地表溫度確定，理論干邊則是由無蒸發(fā)和蒸騰時裸土和植被覆蓋下最高地表溫度確定，由此得到的干濕邊具有明確的物理意義，更加適合植被覆蓋下的水分估算研究。

在確定理論干邊的端點(diǎn)時，地表處于沒有土壤水分蒸發(fā)以及沒有植被蒸騰的狀態(tài)，此時的熱量交換只有顯熱交換，能量平衡方程變?yōu)?/p>

Rn-G=H

(22)

對于裸土或植被分布較少的區(qū)域，如果僅僅選擇NDVI最小時所對應(yīng)的地表溫度為理論最高溫度，可能會出現(xiàn)理論最高溫度低于實(shí)際溫度的情況，因此選擇地表溫度隨著NDVI增大而升到最高的點(diǎn)作為干邊的端點(diǎn)，此時地表溫度記為Tsd。植被完全覆蓋情況下，隨著NDVI升高，地表溫度表現(xiàn)為穩(wěn)定下降趨勢，因此選擇地表覆蓋度fv趨近于1時地表溫度最高值的點(diǎn)作為干邊的另一個端點(diǎn)，地表溫度記為Tvd，其計(jì)算公式分別為

Tsd=

(23)

Tvd=

(24)

式中Tave——地表平均溫度

理論濕邊上的點(diǎn)具有最大蒸發(fā)和蒸騰的特點(diǎn)，熱交換形式為潛熱交換，因此能量平衡方程可表示為

Rn-G=LE

(25)

在選擇濕邊端點(diǎn)時，對于裸土最大蒸發(fā)的點(diǎn)選擇fv接近0時所對應(yīng)的最低地表溫度，記為Tsw，對于全植被覆蓋下最大蒸騰的點(diǎn)則選擇fv接近1時所對應(yīng)的最低地表溫度，記為Tvw，其計(jì)算公式分別為

(26)

(27)

1.4 評價指標(biāo)

主要采用相關(guān)系數(shù)r和均方根誤差RMSE 2個評價參數(shù)進(jìn)行結(jié)果分析和驗(yàn)證，r越大，RMSE越小，表明土壤水分的估算精度越高，地表干濕狀況監(jiān)測

更精確。

另外為了更好地驗(yàn)證與評價結(jié)果，將計(jì)算得到的溫度植被干旱指數(shù)轉(zhuǎn)換為土壤水分估算值，其轉(zhuǎn)換公式為

SMestimated=(1-TVDI)(SMmax-SMmin)+SMmin

(28)

式中SMmax——實(shí)測土壤水分含量的最大值

SMmin——實(shí)測土壤水分含量的最小值

2 結(jié)果與分析

2.1 特征空間構(gòu)建及理論干濕邊端點(diǎn)選取

結(jié)合MODIS遙感影像數(shù)據(jù)及地面觀測數(shù)據(jù)，對計(jì)算得到的NDVI以0.01的步長提取對應(yīng)地表溫度影像上的最大值和最小值，同時對所得到的散點(diǎn)圖進(jìn)行線性擬合得到觀測干濕邊，如圖3所示。

圖3 觀測、理論干濕邊及其特征空間分布圖Fig.3 Observation, theoretical wet and dry edges and their feature space

由圖3可知，DOY112的觀測干濕邊上的點(diǎn)分布較為分散，尤其是在[0.3，0.7]區(qū)間內(nèi)的點(diǎn)，而DOY088上的點(diǎn)分布較為集中均勻，同時研究區(qū)內(nèi)當(dāng)日的溫度較高，最高地表溫度達(dá)到307.44 K，最低地表溫度為295.22 K，而DOY112的當(dāng)日溫度相對較低，最高地表溫度為298.96 K，最低地表溫度為290.16 K，另外DOY088及DOY112上觀測干濕邊線性擬合結(jié)果的相關(guān)系數(shù)均較高(表1)，同時圖3也表示出研究區(qū)內(nèi)理論干濕邊的端點(diǎn)位置(黑色方形實(shí)點(diǎn))及特征空間分布情況，其中綠色實(shí)線為理論干邊，綠色虛線為理論濕邊。由此可見，觀測干濕邊形成的特征空間僅為理論干濕邊特征空間的一部分，根據(jù)理論干濕邊進(jìn)行干旱指數(shù)構(gòu)建以及土壤水分估算將更加貼近實(shí)際情況。

表1 觀測及理論干濕邊相關(guān)參數(shù)Tab.1 Observation and theoretical parameters of wet and dry edges

2.2 土壤含水率估算結(jié)果與分析

根據(jù)式(1)和表1各干濕邊方程分別計(jì)算DOY088和DOY112以觀測干濕邊以及理論干濕邊為特征空間的溫度植被干旱指數(shù)，分別記為TVDIob和TVDIth。本次試驗(yàn)DOY088和DOY112的樣點(diǎn)數(shù)據(jù)分別為37、35個，利用其中的23、22個樣點(diǎn)數(shù)據(jù)(約2/3)分別對計(jì)算得到的溫度植被干旱指數(shù)進(jìn)行擬合分析，其擬合結(jié)果如圖4和表2所示。

圖4 TVDIob和TVDIth與實(shí)測土壤含水率擬合結(jié)果Fig.4 Fitting results between TVDIob, TVDIth and measured soil moisture content

表2 TVDIob和TVDIth與實(shí)測土壤含水率擬合參數(shù)Tab.2 Fitting parameters between TVDIob, TVDIth and measured soil moisture content

圖4a、4b為地表溫度Ts與改進(jìn)植被覆蓋度Fv構(gòu)建的特征空間所計(jì)算得到的溫度植被干旱指數(shù)與實(shí)測土壤含水率擬合結(jié)果，而圖4c、4d為地表溫度Ts與歸一化植被指數(shù)(NDVI)構(gòu)建的特征空間所計(jì)算得到的溫度植被干旱指數(shù)與實(shí)測土壤含水率擬合結(jié)果。從圖中可以看出，不論是DOY088還是DOY112，利用觀測干濕邊及理論干濕邊計(jì)算得到的溫度植被干旱指數(shù)TVDIob和TVDIth與實(shí)測土壤含水率均有比較好的負(fù)相關(guān)性，同時TVDIth與實(shí)測土壤含水率之間的相關(guān)性比TVDIob高，相關(guān)系數(shù)均在-0.700左右，同時均方根誤差也比TVDIob小，均不大于0.060 cm3/cm3。說明利用理論干濕邊得到的溫度植被干旱指數(shù)對土壤水分有著更高的敏感性，更適合植被覆蓋區(qū)域的土壤水分監(jiān)測，主要是因?yàn)槔碚摳蓾襁吀荏w現(xiàn)水熱平衡的實(shí)際狀況。同時在DOY088和DOY112，Ts和Fv所計(jì)算得到的溫度植被干旱指數(shù)均比由Ts和NDVI的計(jì)算值與實(shí)測土壤含水率的擬合效果更好，說明改進(jìn)植被覆蓋度的引入和使用能夠拓寬干旱指數(shù)的應(yīng)用范圍，避免了傳統(tǒng)特征空間干濕邊估算時必須包含裸土、部分植被及全植被覆蓋的地表植被覆蓋狀況的限制，為不完全包含所有的地表覆蓋類型的區(qū)域及影像應(yīng)用提供了技術(shù)支持。

另外，比較DOY088和DOY112的評價結(jié)果可以發(fā)現(xiàn)，DOY088時的擬合結(jié)果均比DOY112相應(yīng)的結(jié)果好，DOY088時TVDIth與實(shí)測土壤含水率擬合結(jié)果相關(guān)性最高，相關(guān)系數(shù)達(dá)到-0.715，均方根誤差也最小，為0.029 cm3/cm3。同時DOY112時散點(diǎn)分布相比較于DOY088分散得多，主要是因?yàn)樾←溤?月處于生長狀態(tài)，其長勢旺盛，同時降水量遠(yuǎn)遠(yuǎn)達(dá)不到小麥的生理需求量，因此當(dāng)時大部分小麥種植區(qū)域進(jìn)行農(nóng)業(yè)灌溉，而有些區(qū)域未能及時進(jìn)行灌溉，使得DOY112時的土壤含水率的分布范圍較大，其擬合精度小于DOY088。

利用另外獨(dú)立的DOY088和DOY112時采集的14、13個樣點(diǎn)數(shù)據(jù)分別對土壤含水率估算結(jié)果進(jìn)行驗(yàn)證評價。利用式(28)將DOY088和DOY112時計(jì)算得到的TVDIth均轉(zhuǎn)換為土壤含水率估算值。其驗(yàn)證結(jié)果如圖5所示。從圖5可以看出，由TVDIth計(jì)算得到的土壤含水率估算值與實(shí)測土壤含水率均有著較好的線性相關(guān)性，其相關(guān)系數(shù)r在DOY088時為0.805，在DOY112時為0.637，均方根誤差分別為0.002、0.018 cm3/cm3。而對于DOY088，Ts-Fv的擬合結(jié)果比Ts-NDVI模式下的好，但是在DOY112時期其差距較小，主要是因?yàn)樵贒OY088時地表覆蓋類型未能完全達(dá)到具備完全植被覆蓋的情況，而DOY112時的地表覆蓋情況對于溫度植被干旱指數(shù)的計(jì)算更加適合，使得不管利用Fv還是NDVI，對土壤水分的估算結(jié)果影響不大。

圖5 基于TVDIth計(jì)算的土壤含水率估算值與實(shí)測值驗(yàn)證結(jié)果Fig.5 Validation results between estimated soil moisture from TVDIth and measured soil moisture content

3 結(jié)論

(1)根據(jù)地表能量平衡方程，提出一種理論干濕邊端點(diǎn)獲取方法及以Ts-Fv特征空間構(gòu)建的改進(jìn)TVDI模型，拓寬了TVDI在干旱監(jiān)測及土壤水分估算中的應(yīng)用范圍，引入改進(jìn)植被覆蓋度參數(shù)，在一定程度上避免了TVDI參數(shù)對地表覆蓋類型的限制。

(2)理論干濕邊所構(gòu)建的特征空間更能反映自然現(xiàn)實(shí)的水熱轉(zhuǎn)換狀況，更適合真實(shí)土壤水分的遙感估算。

(3)改進(jìn)TVDIth與實(shí)測土壤含水率有著較好的線性相關(guān)性，相關(guān)系數(shù)在-0.700左右，均方根誤差不大于0.060 cm3/cm3。DOY088和DOY112的土壤含水率估算值均與實(shí)測土壤含水率有著較好的擬合關(guān)系，尤其是DOY088的反演結(jié)果更加貼近實(shí)際地表干濕狀況。

(4)理論干濕邊端點(diǎn)的選取方法需要多種地表參數(shù)，使得該方法更加復(fù)雜，對地面觀測數(shù)據(jù)的完備性要求較高。因此，在后續(xù)的工作中，將針對小麥不同階段的物候狀態(tài)及不同地表覆蓋類型的研究區(qū)域開展研究和應(yīng)用，進(jìn)一步驗(yàn)證該方法的有效性。