基于修正NDVI時間序列的大區(qū)域冬小麥全生育期墑情監(jiān)測

王金鑫，趙光成，張廣周，于百順，羅蔚然，張成才，李 穎

(1.鄭州大學水利與環(huán)境學院，鄭州 450001；2. 中國氣象局·河南省農業(yè)氣象保障與應用技術重點實驗室，鄭州 450003；3. 河南省氣象科學研究所，鄭州 450003)

0 引 言

遙感墑情監(jiān)測是目前遙感技術應用的前沿領域[1,2]，也是公認的世界性研究難題之一[1-4]。究其原因，影響墑情的因素包括土壤類型、降水、太陽高度、大氣、植物長勢、地表狀況以及田間管理等等，使土壤墑情遙感反演異常困難[1,5,6]。基于植被指數的墑情監(jiān)測或反演是遙感墑情監(jiān)測的常用方法。在一定區(qū)域、一定時間內，氣溫、土壤、地形等因素相對穩(wěn)定，只降雨量變化會對作物生長造成較顯著影響[7]，因而，植被指數的變化與土壤墑情之間有較顯著相關性。但植被受缺水脅迫在短期內不會在生長狀態(tài)上明顯表現出來，具有一定的時滯性[8]，無法及時地反映出土壤水分狀況(后面有詳細討論)，所以植被指數是一個靈敏度較差的水分脅迫指標[9]，但在農作物的整個生育期內，墑情的波動與植被指數的變化頻率與趨勢應該是一致的。除早期的基于歸一化植被指數(NDVI)的距平植被指數(AVI)、條件植被指數(VCI)[10]外，考慮地表溫度影響(其是十分復雜的[11])，人們又提出了植被供水指數(VSWI)[12,13]、溫度植被旱情指數(TVDI)[14]、條件植被溫度指數(VTCI)[15,16]以及特征空間法[14]等新的概念和方法用于土壤墑情監(jiān)測，相應于研究案例都取得了較好的效果。經過多年的研究積累，墑情遙感監(jiān)測也達成一些基本共識：①從土壤深度角度，以10～20 cm左右效果最好[1-3,17]，表層和深層土壤水分與遙感資料的相關關系往往不佳[3]；②從地表覆蓋角度，裸土和低植被以熱慣量模型和微波遙感效果最好[3,4,18]；完全植被覆蓋以作物缺水指數(CWSI)、距平植被指數、條件植被指數、特征空間模型、能量指數法和植被供水指數法等效果較好[1-3,19]；部分植被覆蓋采用農田蒸散雙層模型(但模型復雜)、改進求解的熱慣量模型和溫度調節(jié)指數效果較好[3,18]。概而言之，裸土時，土壤的熱特性起決定作用，基于溫度的方法較好；植被茂密時，植被的生長狀態(tài)是重要的指示因子，指數法占優(yōu)；半植被覆蓋時，要同時考慮土壤和植被熱特性差異的雙重效應，綜合指數或方法最好。充分證明了地表溫度和植被指數是農業(yè)墑情光學遙感監(jiān)測中的主要指示因子[4]。

國內外已有研究表明，任何單一的指數或方法都有自己的特點和適用條件，地表生態(tài)的復雜性和遙感技術的特點也決定了不存在一種普適的墑情監(jiān)測指數和方法。大范圍農田墑情遙感監(jiān)測更是存在遙感數據的適用性選擇以及光熱條件的區(qū)域分異問題。本文擬提出一種基于高時間分辨率、中空間分辨率的遙感時序數據，考慮地域和物候期差異，利用距平植被指數，實現大區(qū)域冬小麥全生育期土壤相對濕度的監(jiān)測方法。以河南省中東部黃淮海平原冬小麥主產區(qū)為研究對象，以MODIS NDVI為基礎墑情指數，針對不同地區(qū)、不同生長階段引進相應因子進行自適應改進，實現對農田墑情及其時空分布規(guī)律的監(jiān)測。

1 研究區(qū)域與數據源

河南省黃淮海平原位于北緯32°08′～36°21′和東經112°51′～116°35′之間，北起豫冀交界，南至淮河一線，西起海拔100 m等高線和豫西北丘陵邊緣，東至豫魯、豫皖分界線(為方便數據處理，本文以該區(qū)域的行政界線為準)，面積約8.7 萬km2，占黃淮海平原總面積的25.6%，河南省總面積的52.4%，屬于南暖溫帶半濕潤大陸性季風氣候區(qū)[20]。該地區(qū)是河南省的主要農業(yè)區(qū)，冬小麥主產區(qū)，年農業(yè)總產值占全省農業(yè)總產值的70%以上。

本文采用的遙感數據包括：研究區(qū)域2011-2015小麥生產年度的b1、b2和b7波段的Terra MODIS數據、MODIS的每日陸面溫度產品等，從共享網站(https:∥ladsweb.nascom.nasa.gov/)下載。 b1、b2和b7波段數據的數據等級是L2級的陸地表面反射率數據，已經過大氣糾正和幾何糾正。b1和b2波段的空間分辨率為250 m，共下載1367幅；b7波段的分辨率為500 m，根據需要下載了75幅；MODIS每日陸面溫度數據(簡稱遙感地表溫度數據，包括每天白天和夜晚兩幅數據，這里取其平均值)是從共享網站下載的成品數據，空間分辨率是1 000 m，數量和使用情況同b1、b2波段。主要問題是有一些日期的數據“漏洞”很多(沒有溫度反演數據)。MODIS源數據為Sinusoidal projection，利用MRT軟件統一轉換為WGS84基準的Albers等積投影，并將地溫數據和b7波段數據重采樣為250 m，以與b1、b2波段匹配，進行像元級計算。

本研究所涉及的墑情及其他數據：研究區(qū)域2014- 2015年度10和20 cm土壤水分含量實測數據、2011-2014年2月上旬和11月下旬的10和20 cm土壤水分含量實測數據、豫南、中、北地區(qū)的多個監(jiān)測站1981-2010年(30年)逐月平均氣溫數據、2014-2015年河南省冬小麥種植區(qū)域數據以及河南省行政界線數據等，從河南省氣象局調研獲取。墑情實測數據是研究區(qū)域均勻分布的100多個站點的10和20 cm每日土壤體積含水量的平均值。以各觀測站30年的逐月平均氣溫進行內插(本文的內插方法均采用基于球面函數的普通克里金插值)，分別得到各區(qū)域每個像元的常年逐月平均氣溫。將2014- 2015年冬小麥種植區(qū)域數據、河南省界線數據也轉換為與上述遙感數據統一的基準與投影，然后對影像數據進行剪裁和植被指數計算。由于冬小麥是河南省的主要糧食作物，種植范圍年際變化不大，所以本文假設2014-2015年的小麥種植范圍同前三年的種植范圍相同，并用之對MODIS數據做掩膜，得到最終的研究區(qū)域。

2 研究方法與過程

2.1 研究思想闡述

中原地區(qū)的冬小麥生長周期較長，從頭年十月到來年六月，跨深秋、冬季、春季和初夏幾個季節(jié)，長達八個月，溫度起伏大，生態(tài)因素多變。要進行大尺度麥田墑情監(jiān)測，既要考慮地域差異，又要考慮冬小麥的物候特征。根據相關研究成果，黃河和駐(馬店)漯(河)交界是河南省重要的綜合地理分區(qū)和小麥生態(tài)分區(qū)界線[21,22]。本文以33°N緯線(駐馬店市與漯河市交界處)和黃河為界，把河南省黃淮海平原冬小麥主產區(qū)劃分為豫南、豫中、豫北三個地區(qū)(圖1)；大致以返青和乳熟期為基點(咨詢河南農業(yè)科學院和國家小麥工程中心專家)把小麥全生育期劃分為：前期(從播種到返青，為裸土到半植被覆蓋期)、中期(從返青到乳熟，全植被覆蓋期)和后期(從乳熟到收割，為麥子變黃時期)，見表1。

表1 研究區(qū)冬小麥生長期劃分Tab.1 Division of growing period of winter wheat in research

圖1 河南省冬小麥主產區(qū)區(qū)域劃分Fig.1 Regional division of main wheat producing areas in Henan Province

在小麥生長前期，麥苗尚不能完全覆蓋裸土，NDVI不能完全反映墑情狀況，必須考慮土壤溫度因子。根據土壤的熱特性，在同一(季節(jié))時間(光熱條件下)，若土壤含水量高，其溫度要降低；若土壤含水量低，其溫度必升高?？紤]到NDVI為一個無量綱的量，這里對原始NDVI進行以下修正：

(1)

在小麥生長中期，麥子生長茂盛，基本全部覆蓋地表，NDVI可以指示根部墑情。但正如前面所述，NDVI相對于墑情波動具有滯后性。圖2是豫中龍城站中期NDVI與10和20cm土壤含水量的曲線圖。NDVI與10、20cm土壤含水量的相關系數分別為-0.218，-0.213；將植被指數曲線前移5d，NDVI與10、20cm土壤含水量的相關系數分別為0.664，0.467?？梢?，NDVI滯后于墑情波動約4～6d。此外，眾所周知，在高植被覆蓋期，NDVI具有飽和性。需注意的是，這里的NDVI是利用遙感波段計算的原始值，由于受到太陽光照角度、觀測角度、云、氣溶膠等因素的影響，其往往包含著很多噪聲，并不都是植被指數的真實值(一般比實際值要小)，尤其是低值部分。

圖2 豫中龍城監(jiān)測站NDVI相對于墑情變化的滯后性Fig.2 The lag of NDVI relative to soil moisture change in Longcheng monitoring station

基于上述考慮，對于中期NDVI修改如下：

(2)

式中：KMNDVI指中期修正NDVI(MediumNDVI)；β為修正加系數；Tn為某像元某天的實際地表溫度；Tn-1為相應像元該天的前一天的實際地表溫度；其他符號同前。

這里修改的前半部分(分子第一項)主要作用是將NDVI的滯后性提前。其理由是如果當天下雨或大面積灌溉，地表溫度肯定低于前一天的溫度，溫差較大時，修改效果明顯；如果不下雨或沒有灌溉，那么二者的溫度接近，修改效果可以忽略。后一部分(分子第二項)的作用與前期類似，這里主要用于改善NDVI的飽和性。

在小麥生長后期，麥子將逐漸變黃，NDVI的值要逐漸變小，但并不一定意味著土壤濕度低。因此也要對后期NDVI進行修正。杜曉等[23]依據對水的吸收率曲線和植被與土壤反射率曲線的分析，指出MODIS的b6、b7波段反映了含水土壤和植被的波譜信息，并構造了地表含水量指數(SWCI)模型。在此基礎上，張紅衛(wèi)等[24]等考慮綠色植被葉綠素對光譜的吸收作用，引入b1波段信息，以b1/b7為修正系數，構造了增強型土壤表層水分含水量指數(ESWCI)模型，并應用于植被茂盛期的墑情監(jiān)測。小麥生長成熟期是一個葉綠素含量逐漸減少，作物由茂盛過渡到黃枯的階段，所以，與上述修正系數相反，本文對該階段的植被指數修改如下：

(3)

式中：KLNDVI指后期修改NDVI(LateNDVI)；γ為修正乘系數；b1、b2、b7為MODIS相應波段光譜反射率。

將修正NDVI與實測墑情數據進行相關分析，驗證其有效性。在通過驗證的基礎上，構建修改NDVI時間序列，并與前三年的平均值進行距平分析，得到監(jiān)測年與前三年平均相比的墑情分布時空規(guī)律。

2.2 技術路線與過程

綜上所述，本文的技術路線如圖3所示。

圖3 研究技術路線Fig.3 Research technical route

本文通過旬最大值合成(MVC)和改進S-G濾波法[25]，去除噪聲數據，構建反映冬小麥生長實際規(guī)律的、較光滑連續(xù)的修正NDVI生長曲線。

3 結果與分析

3.1 修正前后像元生長曲線對比

圖4表示豫中淮陽站2014-2015年植被指數修改前后的情況。

圖4 豫中淮陽監(jiān)測站2014-2015年NDVI修正前后的效果Fig.4 The effect comparisons of NDVI before and after modification of Huaiyang observation station at central Henan province in 2014-2015

可以發(fā)現，修正前后，植被指數曲線的走勢基本是一樣的，但修正后植被指數的變化幅度一般較小。

3.2 修正NDVI與實測墑情的相關分析

為了驗證修正NDVI的有效性，分別在3個地區(qū)均勻選擇8個監(jiān)測站，將修正前后的植被指數與10和20cm實測墑情進行了相關分析。其結果如表2所示。

由表2可知，原始NDVI與實測墑情的相關情況，后期比前期和中期要好，中期比前期要好。修正后NDVI與墑情相關的規(guī)律與原始NDVI基本類似，但修正NDVI在各時期的相關系數都有不同程度提高，說明了修正系數的有效性，修正后的NDVI更適宜作為農田的墑情指數。

表2 NDVI修正前后與實測墑情的相關分析對比Tab.2 Comparison of correlation between NDVI before and after modified and measured soil moisture

續(xù)表2 NDVI修正前后與實測墑情的相關分析對比

注：*表示在0.05水平(雙側)上顯著相關，**表示在0.01水平(雙側)上顯著相關。

3.3 經MVC和S-G濾波后像元生長曲線情況

圖5是豫北某像元修正NDVI曲線情況。

圖5 豫北某像元經MVG合成和改進S-G濾波后的修正NDVI生長曲線Fig.5 The growth curve of a Pixel modified NDVI after MVC synthesis and improved S-G filtering in the north of Henan province

3.4 基于修正NDVI距平分析的冬小麥墑情監(jiān)測

基于生長曲線得到距平值系列(圖6)，并進行距平值等級劃分，最終得到監(jiān)測年(2014-2015年)與前三年均值相比，研究區(qū)域冬小麥田相對濕度的時空分布規(guī)律(圖7)。

圖6 研究區(qū)像元的距平值時間序列Fig. 6 The time series of the anomalies of the study area

圖7 2014-2015年河南省冬小麥主產區(qū)土壤濕度相對于前三年均值偏離情況的時空分布Fig.7 Temporal and spatial distribution of soil moisture relative to the mean of the previous three years in main winter wheat production areas in Henan province from 2014 to 2015

這里的關鍵是距平值(距平植被指數)的等級劃分。實驗數據(圖6)表明，距平值的分布區(qū)間約為[-1,1]，集中在[-0.2,0.2]之內。如果等級劃分過粗，那么就看不出墑情的空間分布結構；如果劃分過細，則會引起數據失真，不能反映真實情況。綜合相關研究文獻[18,26]和實驗結果，本文把距平值劃分為7個等級：<-0.5，[-0.5,-0.15), [-0.15,-0.07)，[-0.07,0.07)，[0.07,0.15)，[0.15,0.5)，>0.5，分別表示相對濕度從干到濕的分布序列。其中，[-0.07,0.07)為與往年持平。需要說明的是，圖7所示的情況是監(jiān)測年相對于前三年平均的濕度偏離情況，并不是絕對的墑情等級。

從圖7中可以看出，2014-2015年的10月下旬、3月中旬、4月中旬至5月上旬，研究區(qū)冬小麥田的墑情基本與前三年平均持平；3月下旬、4月上旬和5月中旬至下旬稍偏干；11月中旬至12月上旬偏干；12月下旬至3月上旬偏濕。

為了驗證上述結果的正確性，我們選取了明顯偏干的2014年11月下旬和明顯偏濕的2015年2月上旬作為實證。將上述兩個時期前三年的10和20 cm實測墑情(取最大值)的平均值與2015年相同時期實測墑情(取最大值)的值做比較(均值為被減數)，其結果如圖8所示(等級劃分的相對比例與圖7一致)。從圖8(a)、8(c)與8(b)，8(d)、8(f)與8(e)的比較可以看出，監(jiān)測值與實測值分布規(guī)律基本是一致的，誤差基本都在一個量級內；比較而言，監(jiān)測值與10 cm實測土壤水分符合度更好。需要說明的是，圖8的色標適用于除8(b)、8(e)以外的其他圖。

圖8 相對濕度驗證實例Fig.8 Comparison of relative moisture verification

4 結 論

土壤墑情是作物旱情診斷和灌溉管理的基礎參數，是農作物產量的限制因子，監(jiān)測土壤墑情對提高水分利用效率和效益具有重要意義[27]。本文充分利用MODIS遙感的高時間分辨率特性，構建基于區(qū)域和作物生育期自適應修改的植被指數生長曲線，通過距平植被指數，實現大區(qū)域全生育期的冬小麥墑情監(jiān)測。研究表明：

(1)與NDVI值相比，基于地域和物候期分異的修正NDVI值與土壤水分的相關性更好，更適合作為農田墑情的標識指數；

(2)基于自適應修正NDVI作物生長曲線的遙感墑情監(jiān)測是一種成本低、精度較好、適宜性較廣的有效方法；

(3)農作物生長曲線是一條隱含諸多農情信息的特征曲線，對其進行農情信息挖掘是一個很有價值的研究方向。

本文提出對中期NDVI修正的初衷之一是將其滯后性提前，但實驗中我們發(fā)現，恰是在降雨的時候，往往也是云層遮擋較為嚴重的時候，MODIS的每日陸面溫度就存在較為嚴重的漏洞。為保證結果的科學性，這時我們采取的措施是不修改。從圖4可以看出，滯后性修改的效果不明顯，這在一定程度上影響了實驗結果。此外，距平植被指數法雖然不能準確地劃定墑情的等級，但如果時間序列積累足夠長，則可直接得到監(jiān)測年相對于常年的偏離情況，具有一定的實用參考價值。