基于VSWI 和TVDI 差異的河套灌區(qū)沈?yàn)豕嘤?耕地灌溉面積遙感監(jiān)測(cè)

田 鑫，李瑞平，王思楠，范雷雷，于明輝

（1.內(nèi)蒙古農(nóng)業(yè)大學(xué) 水利與土木建筑工程學(xué)院，呼和浩特 010018；2.威海市城市規(guī)劃技術(shù) 服務(wù)中心有限公司，山東 威海 264200；3.內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)大數(shù)據(jù)研究與應(yīng)用重點(diǎn)實(shí)驗(yàn)室， 呼和浩特 010018；4.威海市國(guó)土資源地理信息中心，山東 威海 264200）

0 引 言

【研究意義】我國(guó)是個(gè)農(nóng)業(yè)發(fā)達(dá)但水資源匱乏的國(guó)家，伴隨著我國(guó)農(nóng)業(yè)的發(fā)展，水資源合理的分配已成為我國(guó)農(nóng)業(yè)得以健康發(fā)展的基礎(chǔ)，而灌溉面積監(jiān)測(cè)作為水資源合理配置的重要依據(jù)也越來越受到重視[1-2]。沈?yàn)豕嘤蜃鳛楹犹坠鄥^(qū)的重要組成部分，區(qū)域內(nèi)干旱少雨，蒸發(fā)量遠(yuǎn)大于降水量，水資源緊缺、有效灌溉面積增速減緩緩慢的現(xiàn)狀日益嚴(yán)重，因此通過對(duì)沈?yàn)豕嘤虻墓喔让娣e進(jìn)行監(jiān)測(cè)，為沈?yàn)豕嘤蛩Y源合理配置提供理論依據(jù)十分必要。

【研究進(jìn)展】20 世紀(jì)的90 年代初期，杜文才等[3]利用衛(wèi)星遙感技術(shù)對(duì)我國(guó)寧夏中寧地區(qū)的灌溉情況展開了調(diào)查和分析，并繪制出了中寧地區(qū)的水資源、農(nóng)業(yè)灌溉設(shè)施分布圖等。隨后，水利部將河南省的灌區(qū)作為研究對(duì)象，結(jié)合美國(guó)陸地衛(wèi)星4～5 號(hào)專題制圖儀獲取的多波段掃描影像以及河南省現(xiàn)有的數(shù)據(jù)資料等，進(jìn)行實(shí)地的勘察，基本實(shí)現(xiàn)了利用遙感技術(shù)實(shí)現(xiàn)對(duì)有效灌溉面積的測(cè)量[4]。許迪等[5]對(duì)我國(guó)山東省簸箕李灌區(qū)的農(nóng)業(yè)灌溉系統(tǒng)展開了分析和調(diào)查，建立了識(shí)別系統(tǒng)，對(duì)于灌溉面積測(cè)量方法的開發(fā)起到了積極的作用。21 世紀(jì)初期，我國(guó)開展了寧夏引黃灌溉面積及作物種植結(jié)構(gòu)遙感調(diào)查項(xiàng)目研究，結(jié)合遙感技術(shù)、地理信息系統(tǒng)和全球定位系統(tǒng)，實(shí)現(xiàn)了對(duì)寧夏引黃灌區(qū)的灌溉面積、農(nóng)作物分布、土地結(jié)構(gòu)等數(shù)據(jù)的收集和分析。易珍言等[6]、沈靜[7]分別采用對(duì)PDI、MPDI 差值的閾值進(jìn)行確定的方法，結(jié)合少量實(shí)測(cè)數(shù)據(jù)，對(duì)灌溉面積進(jìn)行了提取，取得了較好的成果，但均需提前知曉灌溉行為發(fā)生的準(zhǔn)確時(shí)間。

【切入點(diǎn)】基于以上可以發(fā)現(xiàn)，遙感技術(shù)在灌區(qū)灌溉面積監(jiān)測(cè)上具有廣泛的應(yīng)用，基于指數(shù)差異閾值的灌溉面積監(jiān)測(cè)方法也取得了一定的進(jìn)展，但植被供水指數(shù)VSWI、溫度植被干旱指數(shù)TVDI 等在灌溉面積提取方面的應(yīng)用較少，多用于研究區(qū)的旱情監(jiān)測(cè)等，因此，本文嘗試采用基于Landsat-8 衛(wèi)星數(shù)據(jù)的植被供水指數(shù)VSWI、溫度植被干旱指數(shù)TVDI 結(jié)合地表溫度LST 變化對(duì)研究區(qū)的灌溉面積進(jìn)行提取。

【擬解決的關(guān)鍵問題】探索VSWI、TVDI 在灌溉面積遙感監(jiān)測(cè)模型創(chuàng)建中作為主要參數(shù)的可行性。

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

1.1 研究區(qū)概況

沈?yàn)豕嘤蛭挥诎蛷┠谞柺械奈髂喜浚⒐珮屑~西北部，106°08′E—107°10′E，40°08′N—40°57′N，沈?yàn)豕嘤蚴莾?nèi)蒙古河套灌區(qū)的重要組成部分。灌域東南至黃河與鄂爾多斯市隔河相望，西南與阿拉善盟毗鄰，東北部與解放閘灌域接壤，西北靠陰山山脈，東西長(zhǎng)92 km，南北寬65 km。承擔(dān)著中國(guó)林業(yè)科學(xué)研究院沙林中心，巴彥淖爾市6 個(gè)國(guó)營(yíng)農(nóng)場(chǎng)，磴口縣5個(gè)蘇木、鎮(zhèn)、辦事處，杭錦后旗1 個(gè)鎮(zhèn)，阿拉善盟1個(gè)國(guó)營(yíng)農(nóng)場(chǎng)及鄂爾多斯市杭錦旗巴拉亥鎮(zhèn)部分耕地的灌排任務(wù)，灌域內(nèi)土壤類型以壤質(zhì)砂土和砂質(zhì)壤土居多，并含有少部分的砂土，沈?yàn)豕嘤蚓哂忻黠@的大陸性氣候特點(diǎn)，冬季嚴(yán)寒，夏季短促。氣候干旱少雨，蒸發(fā)強(qiáng)烈，年平均降雨量在130～215 mm 之間，年平均蒸發(fā)量在2 100～2 400 mm 之間。蒸發(fā)量與降水量的嚴(yán)重失衡造成沈?yàn)豕嘤驅(qū)σS灌溉的絕對(duì)依賴。

1.2 數(shù)據(jù)來源

本次研究選用的數(shù)據(jù)為L(zhǎng)andsat-8 衛(wèi)星影像數(shù)據(jù)，共選用了3 景影像，根據(jù)往年灌溉時(shí)間，2016 年7—9月，選取每月云量較少的影像，用于研究區(qū)8 月底第5 次灌溉面積的提取，詳細(xì)情況見表1。

表1 研究所用Landsat OLI/TIRS 影像數(shù)據(jù) Table1 Landsat OLI/TIRS image data for research institute

收集到耕地野外調(diào)查樣本379 份，用于模型構(gòu)建及土地利用類型圖的繪制，其中玉米145 個(gè)，葵花177 個(gè)，瓜類57 個(gè)。非耕地采樣點(diǎn)采用隨機(jī)采樣的方式，共采集沙地樣本點(diǎn)34 個(gè)，建設(shè)用地采樣點(diǎn)25個(gè)，林地采樣點(diǎn)26 個(gè)，同時(shí)于8 月初對(duì)研究區(qū)已灌溉區(qū)域進(jìn)行采樣，共計(jì)182 個(gè)采樣點(diǎn)，用于灌溉面積反演結(jié)果的檢驗(yàn)。

2 模型構(gòu)建

本文以內(nèi)蒙古河套灌區(qū)沈?yàn)豕嘤驗(yàn)檠芯繀^(qū)，基于Landsat-8 數(shù)據(jù)產(chǎn)品對(duì)灌溉面積進(jìn)行監(jiān)測(cè)，構(gòu)建灌溉面積遙感監(jiān)測(cè)模型。通過對(duì)研究區(qū)歸一化植被指數(shù)NDVI、地表溫度LST、植被供水指數(shù)VSWI、溫度植被干旱指數(shù)TVDI 的反演，計(jì)算相鄰影像的VSWI、TVDI 差值，進(jìn)而計(jì)算耕地實(shí)測(cè)樣本點(diǎn)的VSWI、TVDI差值的均值，通過差值的均值來對(duì)每幅影像的灌溉行為進(jìn)行分析，構(gòu)建基于VSWI 差異和TVDI 差異的灌溉面積監(jiān)測(cè)模型，結(jié)合基于SEE5.0 決策樹分類法獲取的耕地邊界矢量圖，對(duì)耕地的灌溉面積進(jìn)行提取。

2.1 歸一化植被指數(shù)NDVI 反演

遙感光譜數(shù)據(jù)在計(jì)算分析生物量、植被生長(zhǎng)狀況上都需要借助植被指數(shù)，是一項(xiàng)無量綱參數(shù)，在光譜數(shù)據(jù)分析中起到了至關(guān)重要的作用。沒有任何條件的輔助以及約束，植被狀態(tài)都可以通過光譜信息表達(dá)出來，同時(shí)能夠?qū)χ脖坏纳L(zhǎng)活力、生物量、覆蓋面積等等這些信息進(jìn)行定性、定量的評(píng)價(jià)和分析[8]。1978年，Deering[9]首次提出了RVI（簡(jiǎn)化比值植被指數(shù)），通過非線性的歸一化處理后，能夠得到NDVI（歸一化植被指數(shù)），并將其約束在[-1，1]范圍內(nèi)，那么可以得到公式：

式中：ρNIR近紅外波段的真實(shí)地表反射率；ρR為紅光波段的真實(shí)地表反射率；NDVI 就是RVI 通過非線性歸一化處理后的結(jié)果。

2.2 地表溫度反演

地表溫度是當(dāng)前我們?cè)u(píng)估大氣層與陸地表層之間的水分以及能量平衡的關(guān)鍵性指標(biāo)之一，同時(shí)也能夠反映出土壤墑值、評(píng)估作物的產(chǎn)量等，因此對(duì)于遙感模型的建立是非常關(guān)鍵的，此外對(duì)于土壤的鹽堿化、保水能力、氣候變化等研究都提供了重要的參考價(jià)值。遙感技術(shù)的發(fā)展推動(dòng)了地表溫度遙感反演算法的精確度和準(zhǔn)確度。目前，大范圍的地表溫度進(jìn)行監(jiān)控常常會(huì)用到熱紅外遙感。地表溫度反演的算法主要有3種：?jiǎn)未八惴ā瓮ǖ浪惴芭八惴╗10-12]，本研究采用劈窗算法，如式（2）所示：

根據(jù)研究區(qū)往年的灌溉時(shí)間，8 月初的第5 次灌溉和10 月底的秋灌，灌溉作物種類基本涵蓋研究區(qū)主要農(nóng)作物，但10 月底前后衛(wèi)星影像云量普遍較高，因此選取7 月24 日，8 月25 日，9 月16 日3 景影像對(duì)地表溫度分布情況進(jìn)行反演，反演結(jié)果如圖1 所示根據(jù)耕地、林地、建設(shè)用地、沙地的實(shí)測(cè)點(diǎn)坐標(biāo)，計(jì)算7—8 月、8—9 月的地表溫差LST7-LST8、LST8-LST9，根據(jù)表2 可以看出，林地、建設(shè)用地、沙地的地表溫度7—9 月始終呈下降趨勢(shì)，而耕地7—8 月地表溫度驟降，降幅較其余3 種土地類型更為明顯，8—9 月耕地地表溫度小幅回升。

根據(jù)以上分析，造成耕地溫度變化的原因是灌溉，進(jìn)一步分析，耕地7—8 月的溫度降低若是降雨引起，首先其溫度降幅應(yīng)低于或持平于其他土地利用類型，其次其地表溫度變化趨勢(shì)應(yīng)與其余3 種土地利用類型基本相同。

表2 不同土地類型地表溫度時(shí)間變化 Table 2 Time variation of surface temperature in different land types ℃

圖1 沈?yàn)豕嘤虻乇頊囟确囱萁Y(jié)果 Fig.1 Inversion results of surface temperature in Shenwu irrigation area

2.3 VSWI、TVDI 差異模型構(gòu)建

Carlson 等[13]首次提出了植被供水指數(shù)（Vegetation Supply Water Index，VSWI）的概念，其是在地表溫度和植被指數(shù)的基礎(chǔ)上建立的一種監(jiān)控模型，能夠?qū)崿F(xiàn)對(duì)植被的干旱情況進(jìn)行觀察。具體的表達(dá)式為：

式中：地表溫度以TS表示，單位為k；歸一化植被指數(shù)為NDVI；圖像增強(qiáng)系數(shù)為B。

Sandholt 等[14]首次提出了溫度植被干旱指數(shù)(Temperature Vegetation Dryness Index，TVDI)的概念，其是在Ts-NDVI 特征空間基礎(chǔ)上發(fā)展而來的。這種干旱指數(shù)能夠計(jì)算出植被指數(shù)以及地表溫度，其數(shù)據(jù)來源于遙感影像，為此可以有效地避免出現(xiàn)土壤水分狀態(tài)監(jiān)控中的問題，大大降低了植被覆蓋度的影響程度，進(jìn)而提高監(jiān)測(cè)數(shù)據(jù)的準(zhǔn)確度，為此對(duì)于地形復(fù)雜、面積較大區(qū)域的干旱監(jiān)測(cè)具有很好的適應(yīng)性。計(jì)算式（4）如示：

式中：濕邊方程以Ts-min表示；干邊方程以Ts-max表示；地表溫度以Ts表示。TVDI 的值域?yàn)閇0，1]。

植被供水指數(shù)和溫度植被干旱指數(shù)可以描述植被的供水情況和土壤的含水率，植被供水指數(shù)越大，植被的供水情況越好，溫度植被干旱指數(shù)越小，土壤含水率越高，這2 種情況都能證明灌溉行為的發(fā)生，因此可以利用以上2 種指數(shù)對(duì)灌溉面積進(jìn)行監(jiān)測(cè)，對(duì)于灌溉面積的提取，關(guān)鍵在于確定灌溉后的影像，基于以上分析可以發(fā)現(xiàn)，灌溉后影像的植被供水指數(shù)VSWI 會(huì)高于灌溉前的影像。而TVDI 會(huì)低于灌溉前的影像，以此可以得出灌溉行為與上述2 種指數(shù)的關(guān)系，如式（5）和式（6）：

式中：Iv、IT分別表示研究區(qū)某個(gè)像元處的灌溉程度；VSWIt1、TVDIt1為前一幅影像的植被供水指數(shù)和溫度植被干旱指數(shù)；VSWIt2、TVDIt2為后一幅影像的植被供水指數(shù)和溫度植被干旱指數(shù)。

在排除降水影響的情況下，對(duì)于植被供水指數(shù)，如果后一幅影像某一處像元的植被供水指數(shù)大于前一幅影像，即Iv小于0，可以認(rèn)為在該像元處，發(fā)生了灌溉行為，且Iv值越小，灌溉行為越明顯。對(duì)于溫度植被干旱指數(shù)來說，如果后一幅影像某一處像元的溫度植被干旱指數(shù)小于前一幅影像，即IT大于零，則可以認(rèn)為在該像元處發(fā)生了灌溉行為。且IT值越大，灌溉行為越明顯。

按照所收獲到的地表溫度的反演結(jié)果以及歸一化植被指數(shù)的運(yùn)算數(shù)據(jù)，以VSWI 運(yùn)算公式為基準(zhǔn)，對(duì)B 進(jìn)行取值，即B=10 000，在運(yùn)算之后得出2016年7 月24 日、8 月25 日、9 月16 日的植被供水指數(shù)，根據(jù)式（3），對(duì)各幅影像的VSWI 進(jìn)行計(jì)算，并根據(jù)式（5），計(jì)算時(shí)間相鄰影像的VSWI 差值，根據(jù)耕地采樣點(diǎn)經(jīng)緯度數(shù)據(jù)，得到耕地采樣點(diǎn)VSWI 差值的均值，結(jié)果見表3。從表3 可以看出，8 月耕地采樣點(diǎn)的植被供水指數(shù)比7 月同比高12.5，相比9 月同比高8.7，基于以上分析，結(jié)合本文構(gòu)建的基于VSWI差異的灌溉面積監(jiān)測(cè)模型可以看出，8 月的灌溉行為十分明顯，可以作為灌溉面積提取的依據(jù)。

表3 耕地采樣點(diǎn)VSWI 差值均值 Table3 Mean value of VSWI difference of cultivated land measured points

故而，按照國(guó)家氣象局氣候司所發(fā)布的相應(yīng)規(guī)范對(duì)干旱等級(jí)進(jìn)行劃分，即重旱、中旱、輕旱、正常、濕潤(rùn)依次為10 cm 土壤相對(duì)濕度<40%、40%～50%、50%～60%、60%～80%、>80%；而基于VSWI 則依次為0～5、5～10、10～15、15～25、25～35。獲取了研究區(qū)7—9 月的旱情分布情況，見圖2。從圖2 同樣可以看出，8 月耕地的植被供水指數(shù)情況確實(shí)相比7月和9 月明顯提高。在獲取沈?yàn)豕嘤蜓芯繀^(qū)地表溫度反演結(jié)果及NDVI 運(yùn)算結(jié)果之后，在對(duì)研究區(qū)TVDI進(jìn)行運(yùn)算的時(shí)候選用了以ENVI 感圖像處理軟件為基礎(chǔ)二次開發(fā)的TVDI 功能模塊。按照NDVI 對(duì)應(yīng)的地表溫度，從起最高、最小值點(diǎn)對(duì)干、濕邊進(jìn)行擬合[15-17]，從而得出該區(qū)域Landsat-8 數(shù)據(jù)7 月24 日、8 月25日、9 月16 日的相應(yīng)擬合曲線，見圖3。

圖2 基于VSWI 的旱情分布 Fig.2 Drought distribution map based on VSWI

圖3 沈?yàn)豕嘤騎S-NDVI 特征空間擬合曲線Fig.3 Ts-NDVI characteristic space fitting curve in Shenwu irrigation area

基于干旱等級(jí)劃分標(biāo)準(zhǔn)把研究區(qū)TVDI 分成了濕潤(rùn)、正常、輕旱、中旱、重旱5 個(gè)等級(jí)，即：0～0.2、0.2～0.4、0.4～0.6、0.6～0.8、0.8～1.0[18]，圖4 為詳盡的結(jié)果。利用式（4）對(duì)3 幅影像的溫度植被干旱指數(shù)進(jìn)行計(jì)算，并根據(jù)野外耕地實(shí)測(cè)點(diǎn)，計(jì)算相鄰影像實(shí)測(cè)點(diǎn)的TVDI 差值的均值。從表4 可以看出，8 月實(shí)測(cè)點(diǎn)的TVDI 均值小于7 月和9 月，根據(jù)基于TVDI差異的灌溉面積監(jiān)測(cè)模型可以看出，8 月的灌溉行為較為明顯，因此選取8 月的影像對(duì)灌溉面積進(jìn)行提取。

表4 耕地實(shí)測(cè)點(diǎn)TVDI 差值均值 Table4 Mean value of TVDI difference between the measured points of cultivated land

3 結(jié)果與分析

3.1 土地利用分類

綜合Landsat-8 影像的光譜特征，植被指數(shù)，ISODATA 分類結(jié)果作為樣本集的屬性值，將研究區(qū)5種主要的地物類型作為訓(xùn)練數(shù)據(jù)集，基于SEE5.0 的決策樹分類方法進(jìn)行土地利用分類，其對(duì)于耕地、林地、沙地、建設(shè)用地的分類精度分別為90.0%、88.5%、88.2%、88.0%，總精度為88.7%，其土地利用類型分類結(jié)果如圖5 所示，進(jìn)而得到沈?yàn)豕嘤蚋厥噶窟吔纭?/p>

圖5 沈?yàn)豕嘤蛲恋乩妙愋?Fig.5 Land use type map of Shenwu irrigation area

圖6 基于耕地矢量邊界的沈?yàn)豕嘤蚬喔让娣e分布 Fig.6 Distribution map of irrigation area in Shenwu irrigation area based on cultivated land vector boundary

表5 灌溉面積監(jiān)測(cè)結(jié)果精度評(píng)定 Table5 Accuracy evaluation of irrigation area monitoring results

3.2 耕地灌溉面積提取

根據(jù)耕地矢量邊界，將基于VSWI 差異模型、TVDI 差異模型所獲取的研究區(qū)第5 次灌溉面積進(jìn)行提取，結(jié)果如圖6 所示可以看出，二者具有較高的一致性，提取的灌溉面積區(qū)域基本相同，同時(shí)根據(jù)統(tǒng)計(jì)數(shù)據(jù)對(duì)2 種模型的精度進(jìn)行對(duì)比，結(jié)果如表5。根據(jù)反演已灌溉點(diǎn)位占灌溉采樣點(diǎn)的百分比，VSWI 和TVDI 的模型精度分別為85.3%和89.7%。

4 討 論

通過指數(shù)閾值選取對(duì)灌溉面積進(jìn)行提取，關(guān)鍵在于對(duì)于灌溉行為發(fā)生時(shí)間的確認(rèn)及閾值的確定，相較之前的成果，本次研究通過地表溫度的變化確定灌溉行為發(fā)生的時(shí)間，針對(duì)河套灌區(qū)等存在大量沙地的區(qū)域可取的較好的成果。今后的TVDI 模型創(chuàng)設(shè)，會(huì)用到一些參數(shù)，譬如地表溫度、植被信息等[19-22]，在遙感參數(shù)繁雜的推算期間或許會(huì)對(duì)一定的經(jīng)驗(yàn)和半經(jīng)驗(yàn)公式加以運(yùn)用，關(guān)鍵參數(shù)反演結(jié)果在地表復(fù)雜度較高的區(qū)域內(nèi)會(huì)出現(xiàn)較大的偏差，從而使TVDI 監(jiān)測(cè)旱情情況的精準(zhǔn)度無法得到有效地確保。因而，在后期探究中，需要將各個(gè)模型參數(shù)的適應(yīng)范疇和運(yùn)算條件理清，結(jié)合多年數(shù)據(jù)，對(duì)模型參數(shù)進(jìn)行進(jìn)一步優(yōu)化，為定量化遙感分析打好基礎(chǔ)。

5 結(jié) 論

2 種模型反演灌溉面積與灌溉面積統(tǒng)計(jì)數(shù)據(jù)的比值分別為90.2%和91.3%，TVDI 的反演結(jié)果更接近于統(tǒng)計(jì)數(shù)據(jù)，根據(jù)反演正確點(diǎn)位與野外采樣點(diǎn)的比值分別為85.3%、89.7%，證明了2 種模型的可行性，基于TVDI 差異的監(jiān)測(cè)模型精度高于基于VSWI 差異的監(jiān)測(cè)模型。