基于無人機遙感的青貯夏玉米水分虧缺指數反演研究

李星恕 程雙飛 薛 志 熊秀芳 韓文霆,2 張立元

(1.西北農林科技大學機械與電子工程學院, 陜西楊凌 712100；2.西北農林科技大學中國旱區農業節水研究院， 陜西楊凌 712100)

0 引言

我國是嚴重缺水國家，農業用水超過全國用水量的60%。玉米作為我國主要的糧食作物，其種植面積達到了0.2億hm2[1]。在干旱半干旱地區，玉米種植主要依靠灌溉，而傳統灌溉方式水資源利用率低、浪費嚴重。為了降低玉米生產過程中的灌溉用水，提高水資源利用效率，實時、準確監測大田玉米旱情非常必要。傳統旱情監測方法一般通過地面傳感器監測土壤含水率和植被生理參數，效率低、成本較高，監測較大面積的玉米旱情信息非常困難[2-4]。遙感技術提供了一種實時、無損的廣域范圍旱情監測方法，但是傳統的衛星遙感時空精度較低，無法滿足大田尺度的旱情信息監測[5]。而無人機遙感由于方便、快捷、精度高、成本低等優點，彌補了衛星遙感的不足，近年來在農業領域廣泛應用[6-8]。

監測旱情的指標很多，其中冠層溫度和以冠層溫度為基礎建立的多種干旱指數應用較多[9-12]，其中水分脅迫指數(Crop water stress index，CWSI)應用最為廣泛。但CWSI干旱指數的計算需要精確獲取作物的冠層溫度[13-15]。遙感監測時，由于熱紅外相機分辨率較低，作物覆蓋度較低的遙感影像中會存在大量的土壤-植被混合像元，從而影響作物冠層溫度的提取[16]。雖然有學者采用基于RGB濾波等預處理方法消除混合像元對冠層溫度提取的影響，但是這些方法需要主觀選擇閾值，存在耗時長、穩定性差等缺點，在實際應用中難以獲得滿意的水分脅迫監測效果[17]。因此，有學者在CWSI的基礎上進一步提出了水分虧缺指數(Water deficit index，WDI)的概念[18]。WDI根據以植物的陸氣溫差(陸地表面混合溫度與空氣溫度之差)和植被指數為橫、縱坐標而形成的植被指數-溫度梯形空間(Vegetation index-temperature trapezoid,VIT)的方法計算得到[19]；通過定義該梯形空間的4個邊界點，在理論上將植物所有干旱情況包含在梯形空間內。與CWSI相比，WDI不需確定作物冠層溫度，因此在植被覆蓋度低的區域仍有較好的反演效果。近年來，有學者利用WDI監測矮化綠豆，發現WDI可監測產量和季節性作物水分變化狀況[20]；利用衛星遙感反演WDI結合潛在蒸散量計算實際蒸散量，發現WDI可用于輔助計算實際蒸散量[21]；利用無人機遙感反演WDI研究水分脅迫對不同基因型蘋果樹的影響，結果表明，WDI作為干旱指標與果樹莖水勢密切相關，適用于不同覆蓋度和冠層結構的蘋果樹旱情監測，可為蘋果育種提供一種新的實時監測方法[22]；對不同生育期的大麥進行無人機遙感WDI反演，表明能夠準確監測不同生育期的大麥旱情變化，且能用于及時發現大麥早熟[23]。近年來雖然利用無人機遙感、采用WDI指數對多種作物進行了干旱脅迫研究，但是利用無人機遙感技術實時監測拔節期青貯夏玉米旱情，研究不同時間尺度下不同水分脅迫對WDI干旱指數監測效果的影響，迄今為止研究較少。

為了研究水分脅迫對不同時間尺度下拔節期夏玉米WDI和陸氣溫差監測效果的影響，本文利用無人機遙感數據和地面采集數據，建立植被指數-溫度梯形空間，計算WDI干旱指數并生成大田尺度的青貯夏玉米旱情分布圖，對基于WDI與陸氣溫差的不同時間尺度下青貯夏玉米氣孔導度和土壤含水率的相關性進行分析，并研究水分脅迫對其影響，為大田作物干旱信息的實時監測提供參考。

1 實驗材料與方法

1.1 實驗區概況

實驗田位于內蒙古自治區鄂爾多斯市達拉特旗昭君鎮(40°26′0.29″ N，109°36′25.99″ E)，海拔1 010 m，面積1.13 hm2。實驗田土壤為砂壤土，田間持水量16%。種植作物為青貯夏玉米(品種為鈞凱918)；播種時間為2018年5月8日，收獲時間為9月8日，全生育期共120 d，拔節期為6月26日—7月26日。播種深度5 cm，行距58 cm，株距25 cm，按照當地種植習慣進行田間管理。實驗田采用自走中心支軸式噴灌機灌溉。

實驗田劃分出3個扇形灌溉區域。為了研究不同水分脅迫程度對青貯夏玉米旱情監測的影響，控制噴灌機灌水量使不同灌溉區域處于100%充分灌溉(田間持水量的95%，1區)、65%充分灌溉(2區)、40%充分灌溉(3區)，如圖1所示。為了方便采集地面數據，在距離實驗田中心點25 m的圓弧上均勻選擇3個采樣區、距離中心點47.5 m的圓弧上均勻選擇6個采樣區(圖1中矩形)，每個采樣區面積6 m×6 m；這樣每個灌溉區內有3個均勻分布的采樣區。為了提高遙感影像精度，采用 RTK x5plus測量系統(科利達，中國)在實驗田內均勻布設5個像控點用于后期遙感影像拼接過程中的圖像校正(圖1中圓形)。在每個矩形采樣區選擇對角線中點和端點3個點作為玉米生理參數的測量點。

圖1 實驗地樣區概況Fig.1 Experimental area

1.2 數據采集

數據采集主要包括地面數據和無人機遙感數據獲取。地面數據采集用于建立WDI梯形空間邊界。無人機獲取遙感數據提取青貯夏玉米的NDVI和陸地表面混合溫度用于計算WDI。

1.2.1無人機數據獲取

采用自主研發的無人機多光譜影像采集系統和無人機熱紅外影像采集系統分別獲取青貯夏玉米的多光譜和熱紅外影像。兩套系統及相機具體參數見表1和表2。

表1 無人機多光譜影像采集系統主要參數Tab.1 Main parameters of UAV multispectral image acquisition system

表2 無人機熱紅外影像采集系統主要參數Tab.2 Main parameters of UAV thermal infrared image acquisition system

在青貯夏玉米拔節期間選擇晴朗無風的天氣進行無人機遙感監測。飛行時間為12：00—13:00，采用固定航線飛行，飛行高度70 m，航向和旁向重疊度為85%，多光譜影像分辨率為4.6 cm，熱紅外影像分辨率為7.8 cm。在后期遙感影像處理過程中，以遙感影像中的像控點為基準拼接多光譜和熱紅外影像，多光譜影像利用黑色和白色漫反射板進行反射率校正；熱紅外影像利用雷泰ST60+型手持測溫儀(RAYTEK，美國)測量黑色漫反射板、白色漫反射板和水體溫度對熱紅外影像進行溫度校準；再利用重采樣的方法將兩類影像像素統一降至8 cm；然后從多光譜影像中提取歸一化植被指數(Normalized difference vegetation index，NDVI)，從熱紅外影像中提取陸地表面混合溫度。歸一化植被指數計算式為

NDVI=(IR-R)/(IR+R)

式中IR——近紅外波段R——紅光波段

1.2.2地面數據獲取

為了獲取建立WDI模型和表征干旱的土壤含水率(Soil water content，SWC)、氣孔導度(Stomatal conductance，Gs)等數據，需要采集以下地面數據：

(1)氣象站數據：為了減小氣象站點和實驗田的氣象數據差異，根據FAO-56建立標準的下墊面和氣象站[24]，從而能夠利用一個氣象站監測更廣范圍內的氣象數據。農業氣象站(中國河北清易電子科技有限公司，中國)位于實驗區南1 km處，下墊面為紫花苜蓿。可采集太陽凈輻射、空氣溫度、相對濕度、風速數據，采集頻率30 min，選擇12:00—14:00的氣象數據計算NDVI-陸氣溫差梯形空間邊界。

(2)土壤含水率：采用干燥法測量實驗田土壤含水率。在采樣區中心點附近采集10 cm深度的土樣，將干燥箱溫度設置105℃干燥8 h至恒質量。

(3)氣孔導度：在采樣區的3個測量點附近選擇有代表性的玉米植株，用AP4型氣孔計(Delta-T公司，英國)測量被太陽直射的玉米第1片全展葉的葉片中間部分的氣孔導度；為了確定凋萎狀態下和充分灌溉下的玉米氣孔導度，在玉米生育期內根據種植經驗選擇實驗田接近凋萎狀態和充分灌溉的玉米植株各10株，測量其氣孔導度并取平均值。選擇最小值作為凋萎狀態的玉米氣孔導度，選擇最大值作為充分灌溉狀態的玉米氣孔導度。

(4)株高：在采樣區3個測量點附近隨機選擇3株玉米測量株高，取平均值。

(5)葉面積指數：在采樣區內側，選擇距離采樣區邊界50 cm的玉米植株，用LAI-2200C型(LICOR，美國)冠層分析儀測量葉面積指數。

(6)溫度：在采樣區3個測量點附近，利用雷泰ST60+型手持測溫儀在距離地表20 cm處垂直測量土壤表面溫度，在玉米株高80%高度處沿水平方向測量冠層溫度，取平均值。

1.3 基于無人機遙感的青貯夏玉米WDI反演方法

采用NDVI-陸氣溫差建立WDI梯形空間，WDI反演技術路線見圖2。通過定義梯形空間的4個邊界點：1點為水分充足的郁閉冠層、2點為凋萎狀態下的郁閉冠層、3點為水分充足的裸土、4點為嚴重水分脅迫的裸土，再結合地面數據計算和確定青貯夏玉米的WDI梯形空間邊界；線段13為WDI模型的濕邊，線段24為WDI模型干邊。通過多光譜影像中提取的NDVI和熱紅外影像中提取的陸地表面混合溫度(與空氣溫度之差為陸氣溫差)，確定青貯夏玉米在梯形空間中的具體位置，見圖3中B點。

圖2 WDI反演技術路線圖Fig.2 Technology road map for inversion WDI

圖3 梯形空間示意圖Fig.3 VIT trapezoid space

水分虧缺指數WDI計算公式為[15]

(1)

式中Ts——陸地表面混合溫度，℃

Ta——空氣溫度，℃

WDI梯形空間的建立是以能量平衡為基礎[19]。作物陸氣溫差的計算公式為[14]

Ts-Ta=

(2)

式中Rn——太陽凈輻射，W/m2

G——土壤熱通量，W/m2

ra——空氣動力學阻力，s/m

rc——冠層阻力，s/m

Cv——空氣體積熱容，J/(m3·℃)

Δ——飽和蒸氣壓與空氣溫度的斜率，kPa/℃

γ——計量常數，Pa/K

VPD——飽和蒸氣壓，kPa

其中Δ、γ、VPD參考FAO-56得到[24]。

利用式(2)計算WDI梯形空間4個頂點時，1、2點分別表示水分充足和凋萎狀態下的郁閉冠層,因此冠層阻力rc分別用潛在蒸散量狀態下的冠層阻力rcp和凋萎狀態下的冠層阻力rex代替；3、4點分別為水分充足和水分脅迫時的裸土，因此可以定義3點和4點的冠層阻力rc分別為零和無窮大。裸土時(3、4點)，熱量傳遞過程主要受到空氣阻力ra和土壤阻力rs的影響，因此計算陸氣溫差時將ra替換為ra+rs。

確定梯形空間所需的參數計算公式如下：

冠層阻力rc[25]

(3)

式中Gs——葉面氣孔導度，mol/(m2·s)

LAI——葉面積指數

空氣動力學阻力ra[26]

(4)

其中

zo≈0.13hcd0=0.65hc

式中z——風速計高度，m

hc——植株最大高度，m

u——風速，m/s

zo——表面粗糙度，m

d0——零點高度，m

土壤阻力rs[27]

(5)

式中rs——土壤阻力，s/m

Tsoil——裸土表面溫度，℃

Tc——冠層溫度，℃b、c——常數

us——土壤表面上方0.2 m處的風速，m/s

其中c=0.002 5,b=0.012,us參考FAO-56利用2 m處風速計算得到[24]。

2 結果與分析

2.1 青貯夏玉米拔節期WDI梯形空間

青貯夏玉米拔節期建立WDI梯形空間邊界所需數據隨生長期的變化如表3所示。表中數據均利用地面測量的數據計算得到。

表3 青貯夏玉米拔節期梯形空間所需數據Tab.3 Silage summer maize jointing stage data of VIT trapezoids

根據表3數據計算生育期第49、57、65天(拔節期3 d)邊界點，根據邊界點的定義，1點和3點為郁閉冠層，NDVI值為1；2點和4點為裸土，NDVI值為0。邊界點陸氣溫差如表4所示。

表4 梯形空間邊界點陸氣溫差

Tab.4 Boundary pointTs-Taof VIT trapezoid℃

生育期1點2點3點4點第49天-7.31.93.068.3第57天-4.11.43.851.9第65天-3.51.43.747.9

根據梯形空間邊界點建立WDI梯形空間邊界，從多光譜遙感影像中提取NDVI，從熱紅外遙感影像中提取地表溫度。青貯夏玉米旱情的WDI梯形空間如圖4所示。

圖4 青貯夏玉米拔節期梯形空間Fig.4 Silage summer maize jointing stage VIT trapezoid

由圖4可以看出，梯形空間的右側干邊可以很好地包含右側所有的像素點，但是左側濕邊只有生育期第49天的梯形空間很好地包含了左側所有像素點，原因是由于生育期第49天空氣濕度明顯低于另外兩天。在3 d的梯形空間中，生育期第57天的梯形空間的值整體偏左，這是由于生育期第55天發生了降雨事件，表明梯形空間對短期降雨反應敏感,可以反映農田整體的干旱情況。

2.2 陸氣溫差和WDI分布圖

陸氣溫差為陸地表面混合溫度與空氣溫度之差，也反映作物的干旱程度。圖5、6為拔節期青貯夏玉米的陸氣溫差和WDI分布圖。由圖5、6可知，生育期第49、57、65天的陸氣溫差和WDI分布圖中均有5條明顯的分界線，是為了方便田間管理在實驗田種植過程中預留的裸地。在分布圖中可以明顯地區分裸地和玉米，且隨著青貯夏玉米的生長，裸地與玉米的差異更加明顯。由于生育期第55天下雨的原因，可以明顯看出生育期第57天陸氣溫差和WDI分布圖整體偏暗，證明陸氣溫差和WDI分布圖對短期降水事件反應敏感。生育期第65天的干旱梯度最明顯，可以看出1區為充分灌溉，WDI值最小。由圖5、6可知，1區和3區有明顯偏暗的區域，這是由于灌溉不均勻導致局部土壤含水率過高，在陸氣溫差和WDI分布圖中形成了明顯偏暗的區域。

圖5 青貯夏玉米拔節期陸氣溫差Ts-Ta分布圖Fig.5 Silage summer maize jointing stage Ts-Ta maps

圖6 青貯夏玉米拔節期 WDI分布圖Fig.6 Silage summer maize jointing stage WDI maps

2.3 日間尺度下WDI與陸氣溫差對土壤含水率的相關性分析

圖7為拔節期青貯夏玉米WDI和陸氣溫差與土壤含水率的線性擬合，關系式如表5所示。由圖7和表5可知，隨著玉米生長，WDI和陸氣溫差與土壤含水率擬合關系式的斜率整體呈減小趨勢；日間尺度下，WDI和陸氣溫差與土壤含水率的相關性無明顯差異。生育期第65天，WDI和陸氣溫差與土壤含水率的相關性明顯高于其他兩天。原因是生育期第65天干旱嚴重，干旱梯度明顯。當干旱越嚴重干旱梯度越明顯，則WDI和陸氣溫差與土壤含水率的相關性越顯著。

圖7 WDI和Ts-Ta與土壤含水率(SWC)的關系Fig.7 Relationship between WDI or Ts-Ta and soil water content

表5 WDI和Ts-Ta與土壤含水率(SWC)的關系Tab.5 Relationship between WDI or Ts-Ta and SWC

2.4 日間尺度下WDI與陸氣溫差對氣孔導度的相關性分析

圖8為拔節期間青貯夏玉米的WDI和陸氣溫差與氣孔導度的線性擬合，關系式如表6所示。由圖8和表6可知，隨著玉米生長，WDI和陸氣溫差與氣孔導度擬合關系式的斜率呈減小趨勢。與陸氣溫差相比，WDI與氣孔導度擬合關系式的斜率變化更加明顯。由表6可知，日間尺度下，WDI和陸氣溫差與氣孔導度的相關性沒有明顯差異；對比拔節期3d的數據可知，干旱越嚴重，干旱梯度越明顯，則WDI和陸氣溫差與氣孔導度相關性越顯著。

2.5 旬間尺度下WDI與陸氣溫差對土壤含水率、氣孔導度的相關性分析

圖9為拔節期青貯夏玉米在旬間尺度下WDI和陸氣溫差與土壤含水率和氣孔導度的線性擬合。由圖9可知，旬間尺度下監測拔節期青貯夏玉米旱情的連續變化時， WDI和土壤含水率、氣孔導度的相關性優于陸氣溫差。原因是玉米在拔節期覆蓋度和生理狀態發生明顯改變，由表3可知，3 d的氣象因素也有明顯的差異，WDI在計算過程中，考慮到了植物生理參數和氣象因素的變化。證明在旬間尺度下，對拔節期青貯夏玉米的旱情進行連續監測時，WDI比陸氣溫差更具優勢。

圖8 WDI和Ts-Ta與氣孔導度的關系Fig.8 Relationship between WDI or Ts-Ta and stomatal conductance

表6 WDI和Ts-Ta與葉片氣孔導度(Gs)的關系Tab.6 Relationship between WDI or Ts-Ta and Gs

圖9 WDI和陸氣溫差與土壤含水率(SWC)和氣孔導度(Gs)的關系Fig.9 Relationships between WDI or Ts-Ta and soil water content or Gs

2.6 旬間尺度下不同水分脅迫的WDI和陸氣溫差與土壤含水率和氣孔導度的相關分析

圖10為拔節期青貯夏玉米在旬間尺度不同水分脅迫下數據的擬合，關系式如表7所示。充分灌溉下(1區100%充分灌溉)，WDI和陸氣溫差與氣孔導度和土壤含水率的相關性都不高。而在不同水分脅迫下(2區65%充分灌溉，3區40%充分灌溉)WDI和陸氣溫差與氣孔導度和土壤含水率都有較強的相關性。對比2區和3區，發現在水分脅迫不同的情況下，WDI與土壤含水率和氣孔導度的相關性無明顯差異，而陸氣溫差與土壤含水率和氣孔導度的相關性卻出現了較大差異。在不同的水分脅迫下，WDI能較穩定地監測青貯夏玉米的旱情，而陸氣溫差對青貯夏玉米的旱情監測效果卻出現了較大的波動。

3 結論

(1)植被指數-溫度梯形空間、WDI分布圖和陸氣溫差分布圖對短期降雨反應敏感，WDI和陸氣溫差分布圖可明顯區分青貯夏玉米拔節期的裸地和玉米，以及灌溉不均勻導致的含水率異常點。

(2)日間尺度下WDI和陸氣溫差與土壤含水率和氣孔導度的相關性無明顯差異。日間尺度下監測青貯夏玉米拔節期旱情在空間上的分布時，易于獲取的陸氣溫差比WDI更具優勢。

圖10 不同水分脅迫下WDI和陸氣溫差(Ts-Ta)與土壤含水率(SWC)和氣孔導度(Gs)的關系Fig.10 Relationship of WDI or Ts-Ta and soil water content or stomatal conductance under different water stresses

水分脅迫自變量x因變量y擬合關系式決定系數R2100%充分灌溉WDITs-TaSWCy=-0.0300x+0.07780.0889Gsy=-0.4889x+0.65590.1139SWCy=-0.0007x+0.07790.0978Gsy=-0.0112x+0.65520.108965%充分灌溉WDITs-TaSWCy=-0.0593x+0.07150.7593Gsy=-1.0105x+0.63460.8200SWCy=-0.0023x+0.07770.6939Gsy=-0.0411x+0.75450.807440%充分灌溉WDITs-TaSWCy=-0.0976x+0.07910.7283Gsy=-1.4127x+0.72200.7725SWCy=-0.0260x+0.07650.3566Gsy=-0.0412x+0.71520.4662

(3)旬間尺度下WDI與土壤含水率和氣孔導度的相關性明顯優于陸氣溫差。在監測旬間尺度下青貯夏玉米拔節期旱情的連續變化時，WDI的監測效果更好；旬間尺度下對不同水分脅迫梯度的數據進行相關性分析發現，充分灌溉下WDI和陸氣溫差都不能很好地監測青貯夏玉米拔節期旱情。WDI在不同的水分脅迫下可以較為穩定地監測青貯夏玉米拔節期旱情，而陸氣溫差的監測效果卻出現了較大的波動，這表明與陸氣溫差相比，WDI監測不同水分脅迫的青貯夏玉米拔節期旱情有著更好的穩定性。

(4)衛星遙感數據、無人機遙感數據、地面監測數據3種數據融合進行監測，可提高監測范圍和監測精度。采用無人機遙感反演WDI監測青貯夏玉米旱情，可為后期“星-機-地”數據融合的新型遙感監測方法提供了研究基礎。