大田玉米作物系數無人機多光譜遙感估算方法

韓文霆 邵國敏 馬代健 ZHANG Huihui 王 毅 牛亞曉

(1.西北農林科技大學機械與電子工程學院，陜西楊凌 712100； 2.西北農林科技大學水土保持研究所，陜西楊凌 712100；3.美國農業部農業研究服務屬，柯林斯堡 CO 80526)

0 引言

作物蒸散量(Evapotranspiration，ET)主要由土壤蒸發(Evaporation)和作物蒸騰(Transpiration)兩部分組成[1]，是連接生態與水文過程的重要紐帶[2]，其快速監測對準確制定和管理大田灌溉制度及提高大田灌溉用水效率有著非常關鍵的作用。常用的估算蒸散量的方法包括：渦度相關法、蒸滲儀法、土壤水量平衡法、 波文比能量平衡法、Penman-Monteith模型等[3-4]。其中PM模型基于能量平衡、空氣動力學原理和冠層阻力來估算ET，具有明確的物理依據，可以清晰地分析ET變化過程及影響機制[2]。但PM模型計算復雜，因此聯合國糧農組織(Food and agricultural organization，FAO)提出了參考作物蒸散量(Reference evapotranspiration，ET0)，由ET0和作物系數估算ET[5]。FAO-56作物系數法[6-7]是世界公認的估算作物蒸散量方法，具有操作簡便、精度可靠、實用性強等特點，在世界范圍內被廣泛應用。

目前，國內外研究人員[1, 3, 8-15]利用蒸滲儀數據、渦度相關數據或者水量平衡法獲得作物的實際蒸散量ETc，并采用氣象數據獲得參考蒸散量ET0，結果已證實上述方法具有很好的精度和普適性，其中雙作物系數法估算蒸散量的精度更好[4]，馮禹等[9]的黃土高原旱作玉米蒸散估算結果表明，引入葉面積指數計算覆蓋度并校正雙作物系數的方法能夠較好地估算旱作玉米ET(R2=0.871)。

由于作物生長狀況、氣象條件和水分脅迫等因素的影響，作物系數必須進行校正才能估算實際蒸散量[7]，但及時發現這些外界因素造成的影響并快速調整作物系數曲線較難，并且在時空尺度上，田間作物用水存在差異性，使用推薦的作物系數計算作物蒸散量誤差較大。作物系數Kc和作物冠層光譜植被指數均受作物覆蓋度、葉面積指數LAI、綠度等因素的影響[16-18]，此關系使遙感技術估算作物系數成為可能。李賀麗等[4]采用地面光譜技術研究了華北平原地區小麥作物系數與植被指數(NDVI、差值植被指數DVI、RVI、SAVI、綠波段葉綠素指數CIgreen、EVI、修正型土壤調整植被指數MASVI和重歸一化植被指數RDVI)的關系以及水分和氮素脅迫對其的影響，結果表明不同施氮水平下Kc無明顯規律性差異，且植被指數和Kc的相關性較弱。HUNSAKER等[8]采用地面光譜技術研究了美國西南部沙漠地區棉花充分覆蓋前和充分覆蓋后的基礎作物系數Kcb與NDVI和積溫GDD的關系。ER-RAKI等[19]基于地面光譜技術研究了采用摩洛哥中心干旱地區小麥的NDVI估算其基礎作物系數Kcb和覆蓋度fc的方法，結果顯示植被指數與Kcb和fc有很高的相關性。GONTIA等[20]基于衛星遙感技術研究了印度西孟加拉邦小麥植被指數(NDVI和SAVI)與其生育期內每月的作物系數Kc的關系。MASELLI等[21]采用意大利中部地區250 m精度的MODIS-NDVI衛星影像估算其植被覆蓋度并結合不同下墊面的氣象站觀測數據估算區域實際蒸散量。

現有研究中，地面光譜技術可以較好地估算基礎作物系數Kcb，但難以估算作物系數Kc；衛星遙感技術可以較好地估算區域作物的作物系數，但衛星的拍攝頻率和影像分辨率難以滿足大田作物的日蒸散量估算要求。無人機遙感技術由于其平臺構建容易、運行維護簡便、分辨率高、作業周期短等特點，廣泛應用于植被信息監測研究中[22-24]，為估算大田作物系數提供了新的解決方案[25]。目前基于無人機遙感技術估算大田作物的作物系數研究較少，且準確估算水分脅迫大田作物的作物系數對制定大田灌溉制度有重要意義。基于此，本文以2017年內蒙古達拉特旗昭君鎮實驗站大田玉米、土壤、氣象等數據為基礎，采用經氣象因子和作物覆蓋度校正后的雙作物系數法計算不同水分脅迫與不同生長時期玉米的作物系數，并使用自主研發的無人機多光譜系統航拍玉米的冠層多光譜影像，提取6種植被指數(SR、EVI、VARI、GNDVI、SAVI和NDVI)，研究不同生長時期(快速生長期、生長中期和生長后期)玉米Kc與6種常用植被指數的關系及水分脅迫對其影響，分析無人機多光譜遙感估算玉米Kc的可行性和適用性，為大田玉米作物系數的遙感估算研究提供技術支持。

1 材料與方法

1.1 實驗區概況

實驗區域位于內蒙古自治區鄂爾多斯市達拉特旗昭君鎮(40°26′0.29″N，109°36′25.99″E)，海拔1 010 m，屬于典型溫帶大陸性氣候。實驗區域玉米(鈞凱918)播種時間為2017年5月20日，出苗時間為2017年6月1日，抽穗時間為2017年7月20日，收獲時間為2017年9月7日(青儲)，生育期(Day after planting，DAP)共110 d。實驗數據采集時間為2017年6月26日—8月29日，對應玉米的生長階段見表1。實驗地塊(圖1中圓形區域)總面積為1.13 hm2，播種深度約5 cm，行距58 cm，株距25 cm，播種時秸稈還田，并按照當地種植習慣施底肥，出苗后根據出苗情況進行補苗，并在玉米拔節期、大喇叭口期分別追加氮肥(240 kg/hm2)。實驗階段總降雨量為44 mm，日平均氣溫22.9℃，日平均風速0.5 m/s，日平均相對濕度為57.1%，日平均太陽凈輻射為145.1 W/m2。生育期內自然降雨難以滿足作物需水要求，主要的水分供給方式為噴灌。實驗區土壤為砂壤土，其中砂粒(0.05～2 mm)占80.7%、粉粒(0.002～0.05 mm)占13.7%、黏粒(0～0.002 mm)占5.6%，耕層土壤有機質質量比為47.15 g/kg，田間持水率為29%(體積含水率)，土壤容重為1.56 g/cm3。實驗地中選取2個扇形區域(圖1)，在區域1和區域2內分別選取一塊10 m×10 m的方形樣地(圖1中A和B區域)，并對樣地A和樣地B的玉米進行長期監測(2017年6月26日—8月29日)。

表1 2017年玉米生長階段劃分Tab.1 Maize growth stage division in 2017

注：玉米生長階段劃分依據參考FAO-56指南[7]。

圖1 實驗地俯視圖(2017年8月15日，無人機RGB影像，5 cm/pixel)Fig.1 Experimental ground view (2017-08-15,unmanned aerial vehicle RGB image, 5 cm/pixel)

1.2 水分脅迫處理

實驗地灌溉方式和灌溉量測量方法如圖2所示，采用中心支軸噴灌機(圖2a)灌溉，在玉米不同生育期進行不同的水分脅迫處理，如表2所示，田間持水量的50%為充分灌溉，采用雨量筒(圖2b)測量樣地A和樣地B的灌溉量和降雨量，每個樣地內均勻放置3個雨量筒，取其平均值。

圖2 灌溉方式和灌溉量測量方法Fig.2 Irrigation method and irrigation depth measuring method

1.3 無人機系統及影像數據采集方法

采用自主研發的無人機多光譜影像采集系統(圖3)，該系統采用開源飛控Pixhawk控制，經緯 M600型機架，最大載量5 kg，最大續航時間30 min。搭載的多光譜相機(圖3b)為RedEdge(MicaSense，USA)，該相機有5個波段(表3)，相機焦距為5.5 mm，視場角為47.2°，圖像分辨率為1 280像素×960像素。該相機配備了兩塊3 m×3 m的灰板(GroupVIII，USA)(圖3c)及光強傳感器(圖3d)。光強傳感器可校正航拍過程中外界光線的變化對光譜影像造成的影響，灰板具有固定的反射率，如表3所示，可對航拍影像進行反射率校正，生成反射率影像圖，并提取植被指數。

表2 實驗區實際灌溉量與降雨量Tab.2 Actual irrigation amount and rainfall in experimental area mm

圖3 無人機多光譜影像采集系統Fig.3 UAV multispectral images acquisition system

表3 RedEdge多光譜相機參數及灰板對其中心波長的反射率Tab.3 Multispectral camera parameters and reflectivity of gray plate to its center wavelength

選取晴朗無云天氣，于11:00—13:00拍攝玉米冠層多光譜影像，拍攝時鏡頭垂直向下，如圖4所示，每隔3 d拍攝一次，飛行高度70 m，地面分辨率為5 cm/pixel，航向和旁向重疊度均為80%，每次航拍采用固定航線。無人機作業時，首先飛至作業高度并拍攝地面兩塊灰板，然后進行航拍作業。整個生育期共采集20組多光譜影像數據。

圖4 無人機航拍作業現場圖Fig.4 UAV aerial photography

1.4 地面數據采集方法

(1)氣象站數據采集：農業氣象站由中國河北清易電子科技有限公司組裝和調試，位于實驗區西南100 m處(40°25′55.53″N，109°36′22.69″E)，其下墊面為青草，草平面高度保持在12 cm左右。監測參數包括空氣溫度、相對濕度、風速、太陽凈輻射和降雨量，每30 min采集一次數據。

(2)土壤含水率測量：采用干燥法測量實驗區玉米的土壤含水率變化。在A和B樣地中心點附近隨機選取3組樣點，測量180 cm土層含水率變化。測量時，按照30 cm間隔對土壤進行分層采樣，每3 d采集1次，并在灌溉前、灌溉后和降雨后進行加測，取180 cm厚土層各采樣層土壤含水率的平均值。采用線性插值法對土壤含水率進行處理，得到樣地A和樣地B土壤含水率的逐日數據序列。

(3)玉米生態指標觀測：采用隨機采樣法對玉米的葉面積指數(LAI)、株高等指標進行采集，每7 d測量一次。在樣地A和樣地B內隨機選取10個采樣點，采用LAI-2200C型冠層分析儀測量葉面積指數(LAI)，取其平均測量值；在樣地A和樣地B內隨機選取10株具有代表性的玉米，采用米尺測量株高，取其平均測量值。

1.5 作物系數與植被指數計算方法

1.5.1FAO雙作物系數法

FAO-56[7]中指出雙作物系數法的計算流程。首先將作物的生育期分成4個生長階段：生長初期、快速生長期、生長中期和生長后期；然后在FAO-56中查找該作物類型在標準情況下(無任何脅迫，田間管理良好，最小相對濕度為45%，冠層上方2 m高度風速為2 m/s)不同生長階段的基礎作物系數推薦值Kcb,tab，可得到標準情況下的基礎作物系數逐日曲線；標況下玉米生長初期、生長中期和生長末期的基礎作物系數分別取0.15、1.15、0.15。根據實驗區氣候條件、作物生理參數(作物高度、LAI等)、土壤水分、田間管理等信息對標況下的基礎作物系數曲線進行校正，并按照FAO-56給出的算法[7]計算水分脅迫系數Ks和土壤蒸發系數Ke，計算公式為

(1)

Ke=Kc-KsKcb,a

(2)

其中

Kcb,a=Kc,min+(Kcb,adjust-Kc,min)[1-exp (-0.7LAI)]

(3)

(4)

式中Kcb,adjust——不同生長階段經氣象因子校正后的基礎作物系數

Kcb,tab——標準情況下的基礎作物系數

u2——不同生長階段冠層上方2 m處的日平均風速，m/s

RH,min——不同生長階段日平均最小相對濕度，%

Kcb——不同生長階段經水分脅迫系數校正后的基礎作物系數

h——不同生長階段平均株高，m

Kcb,a——根據作物實際覆蓋度校正的基礎作物系數

Kc——作物系數

Kc,min——裸土最小作物系數，取值為0.15～0.20

LAI——實測的冠層葉面積指數

1.5.2玉米植被指數提取方法

植被指數可以對地表植被狀況進行簡單、有效地度量，前人研究中也使用了植被指數建立作物系數的估算模型。常用的植被指數中，NDVI可以估算植被的生理生態參數(覆蓋度、葉面積指數、光合有效輻射等)，SR可監測植被覆蓋變化且在植被覆蓋濃密的情況下效果最好，SAVI可降低土壤背景的影響，EVI可解決植被指數易飽和的問題，VARI可減小光照差異和大氣效應的影響，GNDVI對作物的色素變化敏感[26-27]。因為作物系數Kc綜合考慮了作物蒸騰和土壤蒸發兩方面的影響因素，所以可監測作物系數變化情況的植被指數應能較好地反映植被和土壤的變化情況。

因此本文選取了6種植被指數(NDVI、SR、SAVI、EVI、VARI和GNDVI)，其計算公式見表4。采用Pix4DMapper軟件平臺對RedEdge相機拍攝的多光譜影像進行幾何校正、高斯均值濾波和拼接處理，生成正射影像圖，并利用灰板和光照傳感器對光譜影像進行反射率校正，生成反射率影像。采用ENVI軟件平臺裁剪樣地A和樣地B的反射率影像，并提取上述6種植被指數。計算時，灰板對多光譜相機各波段中心波長的平均反射率分別作為藍、綠、紅、近紅外波段的反射率。

表4 植被指數計算公式Tab.4 Vegetation indices

注：RBlue、RGreen、RRed、RNir分別為灰板對RedEdge相機藍、綠、紅、近紅外波段的平均反射率。

2 結果與分析

2.1 玉米生理指標

快速生長期至生長后期樣地A和樣地B玉米株高和LAI如圖5和圖6所示。

圖5 樣地A和樣地B玉米株高變化Fig.5 Variation of maize height in sample fields A and B

圖6 樣地A和樣地B玉米LAI變化Fig.6 Variation of maize LAI in sample fields A and B

圖5中，樣地A和樣地B在快速生長期玉米高度逐漸增大，且兩樣地的株高差異較小，生長中期達到最大值，最大高度分別為273 cm和268 cm。圖6中，樣地A和樣地B在快速生長期玉米LAI逐漸增大；生長中期兩樣地的LAI都先增大再減小，樣地A和樣地B玉米LAI最大值分別為3.07和2.70，兩樣地LAI的整體變化規律符合作物的生長規律；生長中期和生長后期由于樣地B的灌溉頻率和灌溉量小于樣地A的值，樣地B出現水分脅迫，所以樣地B的LAI小于樣地A的值。

2.2 植被指數

快速生長期至生長后期樣地A和樣地B玉米植被指數如圖7和圖8所示。

圖7 樣地A植被指數變化曲線Fig.7 Curves of vegetation indices in sample field A

圖8 樣地B植被指數變化曲線Fig.8 Curves of vegetation indices in sample field B

由圖7和圖8可知，樣地A和樣地B在快速生長期植被指數逐漸增大，生長中期和后期植被指數NDVI、EVI、SAVI、GNDVI、VARI有緩慢降低趨勢，植被指數SR有較大的減小趨勢；生長中期和后期，水分脅迫使得樣地B的葉面積指數小于樣地A，植被指數SR對作物的覆蓋度變化比較敏感，因此樣地B的植被指數SR比樣地A減小幅度大。

2.3 作物系數

快速生長期至生長后期樣地A和樣地B玉米的基礎作物系數Kcb、土壤蒸發系數Ke、水分脅迫系數Ks和作物系數Kc的逐日變化曲線如圖9所示。

圖9 快速生長期至生長后期作物系數逐日變化曲線Fig.9 Curves of daily crop coefficient from rapid growth stage to late growth stage

圖9a中，樣地A和樣地B玉米基礎作物系數Kcb在快速生長期逐漸增大，生長中期Kcb達到穩定的較大值，生長后期Kcb逐漸降低，該變化規律主要由于作物覆蓋度變化所致；生長中期和生長后期，由于水分脅迫的影響，樣地B的作物蒸騰受限，所以其Kcb小于樣地A的值；生長中期和后期的Kcb變化是土壤水分條件和作物長勢狀況等因素綜合影響的結果。圖9b中，快速生長期和生長中期土壤蒸發系數Ke逐漸減小，主要由于作物覆蓋度增大限制了土壤蒸發所致；生長后期樣地A的Ke有增大趨勢，樣地B的Ke逐漸減小，主要原因在于生長后期作物覆蓋度逐漸降低，充分灌溉條件下有利于土壤蒸發，水分脅迫會限制土壤蒸發；生長后期的Ke變化也是土壤水分條件和作物長勢狀況等因素綜合影響的結果；Ke出現的波動變化主要因為灌溉調整土壤含水率所致。圖9c中，快速生長期至生長后期，由于樣地A無水分脅迫，所以其水分脅迫系數Ks保持不變，樣地B有水分脅迫，隨著作物蒸騰和土壤蒸發，其水分脅迫系數逐漸降低。圖9d中，樣地A的作物系數Kc在快速生長期逐漸增大，生長中期保持相對穩定的較大值，生長后期逐漸減小，此變化規律和FAO-56中標況下的Kc值變化情況基本相符。樣地B的Kc在快速生長期有一定的波動變化，是水分脅迫和土壤蒸發綜合作用結果；生長中期和生長后期作物系數Kc開始快速下降，主要由于較長時間的水分脅迫(表2)使得水分脅迫系數Ks持續下降所致。上述結果與馮禹等[9]在黃土高原東部地區春玉米試驗結果較為一致。

2.4 玉米植被指數與Kc關系模型

快速生長期樣地A和樣地B玉米植被指數與作物系數Kc的關系如圖10、11所示。

圖10 快速生長期樣地A植被指數與作物系數Kc的關系(n=25)Fig.10 Relationship between vegetation indices and Kc of maize in sample field A in rapid growing stage (n=25)

圖11 快速生長期樣地B植被指數與作物系數Kc的關系(n=25)Fig.11 Relationship between vegetation indices and Kc of maize in sample field B in rapid growing stage (n=25)

由圖10、11可知，快速生長期樣地A玉米植被指數與作物系數Kc的相關性較強(R2為0.731 2～0.940 1，p<0.05，RMSE為0.025 0～0.053 0，n=25)，見表5，其中比值植被指數SR與Kc的相關性最好。主要原因在于快速生長期玉米的植被指數(圖7)逐漸增大，灌溉充足玉米的Kc也逐漸增大，因此植被指數和作物系數具有較好的相關性。樣地B玉米植被指數和作物系數Kc相關性較弱(R2為0.000 2～0.083 0)，因為快速生長期樣地B在水分脅迫和作物覆蓋度變化的影響下其Ke波動較大，使得Kc波動也較大，而植被指數變化趨勢較平滑(圖8)，所以兩者相關性較弱。

表5 快速生長期樣地A玉米植被指數與作物系數Kc的關系Tab.5 Relationship between vegetation indices and Kc of maize in sample field A in rapid growing stage

生長中期至生長后期樣地A和樣地B玉米植被指數與作物系數Kc的關系如圖12、13所示。

圖12 生長中期至生長后期樣地A植被指數與作物系數Kc的關系(n=40)Fig.12 Relationship between vegetation indices and Kc of maize in sample field A from middle growing stage to later growing stage (n=40)

圖13 生長中期至生長后期樣地B植被指數與作物系數Kc的關系(n=40)Fig.13 Relationship between vegetation indices and Kc of maize in sample field B from middle growing stage to later growing stage (n=40)

由圖12、13可知，生長中期至生長后期樣地A玉米植被指數與作物系數Kc相關性(R2為0.276 5～0.373 2)較弱，因為生長中期至生長后期樣地A玉米的覆蓋度較高，植被指數已達到飽和，對作物系數變化不敏感，所以兩者的相關性較弱。樣地B部分植被指數與作物系數Kc的相關性(R2為0.366 2～0.848 7，p<0.05，RMSE為0.114 2～0.233 8，n=40)較強，見表6，主要因為該生長中期至生長后期樣地B有水分脅迫，玉米覆蓋度逐漸減小，部分植被指數減小趨勢較明顯，所以兩者相關性較高。

表6 生長中期至生長后期樣地B植被指數與作物系數Kc的關系Tab.6 Relationship between vegetation indices and crop coefficient Kc of maize in sample field B from middle growing stage to later growing stage

由以上可知，不同生長時期玉米的比值植被指數(SR)和作物系數相關性最好(R2最大為0.940 1，p<0.05，RMSE=0.025 0，n=25)，本研究相對于前人研究取得了較好的結果，李賀麗等[4]基于地面光譜技術建立小麥的VIs-Kc模型，決定系數R2最大為0.15，p<0.01，n=195，PARK等[34]基于衛星遙感技術建立NDVI、LAI和土壤含水率與Kc模型，決定系數R2最大為0.64。本文采用的水分脅迫程度較大，模型的適用性有待進一步驗證，后期實驗應劃分多個水分脅迫水平，對多種水分脅迫玉米的Kc估算精度進行評價。時間上非連續采集土壤含水率，對模型也存在一定的影響。目前水分脅迫的監測手段有熱紅外影像[35]、地面傳感器網絡等。后期實驗可同時采集多光譜影像和熱紅外影像數據，并搭建地面傳感器網絡連續監測作物、土壤等相關參數，對大田玉米作物系數的無人機遙感技術估算模型進一步修正。

3 結論

(1)采用無人機多光譜技術估算Kc具有一定的可行性。在快速生長期充分灌溉和生長后期水分脅迫條件下，均可以采用無人機多光譜技術估算玉米作物系數Kc。

(2)不同時期玉米的植被指數和作物系數的相關性不同。快速生長期充分灌溉玉米的VIs-Kc模型的相關性(R2為0.731 2～0.940 1，p<0.05，RMSE為0.025 0～0.053 0，n=25)與生長中期至生長后期充分灌溉玉米的VIs-Kc模型的相關性(R2為0.276 5～0.373 2，p<0.05，RMSE為0.073 8～0.079 3，n=40)不同；快速生長期水分脅迫玉米的VIs-Kc模型的相關性(R2為0.000 2～0.083 0，p<0.05，RMSE為0.099 4～0.114 4，n=25)與生長中期至生長后期水分脅迫玉米VIs-Kc模型的相關性(R2為0.366 2～0.848 7，p<0.05，RMSE為0.114 2～0.233 8，n=40)不同。

(3)不同水分脅迫玉米的植被指數和作物系數的相關性不同。快速生長期充分灌溉玉米的VIs-Kc模型相關性(R2為0.731 2～0.940 1，p<0.05，n=25)比水分脅迫玉米的VIs-Kc模型的相關性(R2為0.000 2～0.083 0，p<0.05，n=25)強；生長中期至生長后期充分灌溉玉米的VIs-Kc模型相關性(R2為0.276 5～0.373 2，p<0.05，n=40)比水分脅迫玉米的VIs-Kc模型相關性(R2為0.366 2～0.848 7，p<0.05，n=40)弱。

(4)不同植被指數和作物系數相關性不同，比值植被指數SR與作物系數Kc的相關性最好。快速生長期充分灌溉玉米的VIs-Kc模型的相關性由大到小依次為：SR、EVI、VARI、GNDVI、SAVI、NDVI；生長中期至生長后期水分脅迫玉米的VIs-Kc模型的相關性由大到小依次為：SR、GNDVI、VARI、NDVI、SAVI、EVI；其中比值植被指數SR與作物系數Kc的相關性最好。