基于無人機數碼影像的冬小麥氮素營養診斷研究

楊福芹，李天馳，馮海寬 ，謝 瑞，肖天豪，陳 超

（1. 河南工程學院土木工程學院，河南 鄭州 451191；2. 國家農業信息化工程技術研究中心，北京 100097；3. 河南省測繪工程院，河南 鄭州 451191）

0 引言

【研究意義】氮素是冬小麥生長所必需的營養元素之一。現有研究表明，氮素能調控冬小麥不同生育期的生長狀況，對冬小麥的產量、蛋白質含量等有著顯著影響[1]。因此及時監測作物的氮素營養狀況對于其生長有著重要意義[2?3]。【前人研究進展】目前測定冬小麥氮含量的方法有SPAD-502葉綠素測量儀[4]、Auto Analyzer3型連續流動分析儀[5]、Dualex氮素指數測量儀[6]等。但是這些方法費時費力，而且測定成本較高。近來發展迅速的測繪無人機以體積小、易操作等優點廣泛應用于各方面研究[7?9]。在林業研究方面，楊全月等[10]、金忠明等[11]利用無人機影像進行樹冠信息的研究工作。周小成等[12]利用兩期無人機影像對伐區蓄積量進行估算研究。在水環境研究方面，侍昊等[13]利用微型無人機數據對城市水環境信息進行提取研究，陳橋驛等[14]把無人機航測技術用于河道水環境治理中，馮磊等[15]利用無人機遙感技術監測海域陸源排污口。在農業方面，通過無人機搭載不同的傳感器，可進行農作物生長監測[16]、產量預測[17]、作物表型參數提取[18]等工作。李長春等[19]利用無人機數碼影像建立大豆育種材料葉面積指數預測模型，有效地估測了大豆育種材料不同生育期的葉面積指數。陳鵬等[20]利用無人機搭載的多光譜相機獲取了馬鈴薯的影像信息與紋理信息，通過二者融合估算了馬鈴薯葉片葉綠素含量。田明璐等[21]將成像光譜儀搭載于低空無人機，使用多種模型對棉花高光譜影像進行反演，估測了棉花葉面積指數。Mengmeng D等[22]通過無人機遙感影像獲得了9種植被指數，利用這些植被指數采用逐步回歸法構建了冬小麥產量估算模型。近年來，相關學者提出了用于診斷作物氮營養狀況的氮營養指數（Nitrogen nutrition index, NNI）。該指數是基于作物實際氮濃度值與臨界氮濃度值建立的數理模型[23]，有效地結合了作物生物量與氮含量，診斷效果要優于使用單一作物農學參數[24]。國內外學者將氮營養指數廣泛應用于各類作物的氮營養診斷，王新等[25]建立了加工番茄氮素營養指數模型，驗證結果與種植區實測診斷結果一致。王仁紅等[26]建立了冬小麥關鍵生育期的氮營養指數經驗模型，模型決定系數（R2）高于0.85。CHRISTOS A等[27]對3個亞麻籽品種施加不同濃度的氮肥，發現植株氮素利用效率與NNI之間呈負相關關系。【本研究切入點】目前基于無人機開展氮素營養診斷的研究，多是基于氮含量或氮累積量進行研究，從氮營養指數角度進行氮素營養狀況的研究比較缺乏。以往研究多是直接使用估算模型進行圖像指數的篩選，采用方差膨脹因子結合相關性篩選相關性好、圖像指數間多重共線性弱的圖像指數方法也較少見。【擬解決的關鍵問題】本研究以冬小麥為研究對象，利用無人機高分辨率數碼影像提取圖像指數，然后結合相關性和方差膨脹因子，綜合選取相關性高、多重共線性小的圖像指數，探討利用無人機數碼影像估算冬小麥氮營養診斷模型的可行性及效果，以期為實時監測冬小麥氮素營養狀況及生長狀況提供有力保障。

1 材料與方法

1.1 試驗設計

研究區位于北京市昌平區小湯山國家精準農業研究示范基地，所選取的試驗田有48個小區，每個小區面積是56 m2。試驗品種為京9843（J9843）和中麥175（ZM175）；4個氮素水平分別為未施尿素（0，N1）、195 km·hm?2尿素（1/2正常，N2）、390 kg·hm?2尿素（正常，N3）、585 kg·hm?2尿素（3/2正常，N4）；3個水分水平分別為雨養（W1）、正常水（W2）、2倍正常水（W3）。16個小區為一組，每種處理重復3次，重復1（編號1-1到13-1）和重復2（編號1-2到13-2）用于建模，重復3（編號1-3到13-3）用于 驗證。試驗方案如圖1所示。

1.2 農學參數獲取

試驗分別在冬小麥挑旗期（2015年4月26日）和開花期（2015年5月13日）獲取小麥葉片氮、植株 氮以及生物量等農學參數。

1.2.1 地上干生物量 在室外調查固定樣方2行×1 m（行距15 cm）的每個小區內，破壞性選取具有代表性的20株小麥作為樣本，器官分離后迅速稱其質量后分別置于紙袋中，置于烘箱105 ℃殺青30 min，然后設置溫度到75 ℃烘干至恒定質量（48 h以上）。計算含水量，根據已獲得的地上鮮生物量計算地上干 生物量。

1.2.2 植株氮含量 將烘干后的各器官稱重后進行粉 碎，采用凱式定氮法測定不同植株器官的氮含量。

圖 1 試驗設計Fig. 1 Experimental design

1.2.3 氮營養指數 基于冬小麥植株氮濃度和地上部生物量和JUSTES等[28]建立的冬小麥臨界氮濃度稀釋模型[式（1）]，來計算冬小麥氮營養指數[式（2）]。

式中， Nct表示植株臨界氮質量分數（%）；DM表示地上部生物量（t ·hm?2）。

式中， Nct表示植株臨界氮質量分數（%）；N為小 麥地上部分實際植株氮質量分數（%）。

1.3 無人機數碼數據獲取及處理

采用八旋翼無人機為近地飛行平臺搭載數碼相機，數碼相機型號是DSC-QX100，質量0.179 kg，像素達到2 090萬，選擇晴朗無云天氣，在12:00～13:00對冬小麥試驗區進行航拍。將獲取的數碼影像數據采用俄羅斯Agisoft LLC公司開發的Agisoft Photo Scan Professional軟件處理，主要步驟為：（1）數碼影像和POS點數據輸入；（2）點云生成；（3）網格構建；（4）紋理生成；（5）生成無人機數字正射影像（Digital Orthophoto Map，DOM）。最后獲取的無 人機影像如圖1所示。

1.4 數碼圖像指數的選取

可見光波段是植物葉片的強吸收波段，由于植物色素的吸收，特別是葉綠素a和b的強吸收，在可見光波段形成兩個吸收谷（藍光和紅光）和一個反射峰（綠光）。圖像指數是將不同波段的DN值按照一定的數學運算規則進行運算。氮素是葉綠素的主要組成成分之一，葉綠素含量直接影響著氮素的豐缺。根據文獻和資料中已有的可見光植被指數[29?36]，以及定義的r、g、b圖像指數，共選25個圖像指數對 氮營養指數進行反演，如表1所示。

1.5 數據分析方法

1.5.1 方差膨脹因子 在衡量多元線性回歸模型中變量間多重共線性嚴重程度時通常使用方差膨脹因子（Variance Inflation Factor，VIF），它表示解釋變量之間的相關性[37]。其定義為：

式中， Ri2表示第i個變量與其他變量的決定系數。一般來說，當0≤VIF＜10時，回歸模型不存在多重共線性；當10≤VIF≤20時，回歸模型存在一定的多重共線性；當VIF＞20時，回歸模型存在嚴重 的多重共線性[38]。

1.5.2 偏最小二乘法 偏最小二乘法（Partial Least Square Regression，PLSR）其實質是將典型關聯分析、主成分分析及多元線性回歸分析于一身，其表達式為：

式中，yi為氮營養指數，i=1，2, ···, n，a0、a1、a2、 ···、an為回歸系數，x1、x2、 ···、xn為圖像指數。

1.5.3 精度評價 為了評價模型的效果，選取決定系數（Correlation Coefficient, R2）、均方根誤差（Root Mean Square Error, RMSE）及平均絕對誤差（Mean Absolute Error, MAE）作為結果精度評價指標。計算公式為：

表 1 研究用到的圖像指數Table 1 Image indices applied

2 結果與分析

2.1 圖像指數敏感性的篩選

式中，n可表示為樣本總數， Xi、 Yi可分別表示為第i個樣本氮營養指數實測值與預測值， X 、 Y可表 示為樣本氮營養指數實測平均值與預測平均值。

2.1.1 圖像指數與氮營養指數的相關性分析 將無人機數碼影像的圖像指數與冬小麥氮營養指數進行相關性分析，結果如表2所示。從表2可以看出，冬小麥挑旗期、開花期的圖像指數大部分都與氮營養指數存在極顯著相關（P＜0.01）。在挑旗期，圖像指數r、r/b、g/b、r?b、g?b、IKAW、（r?g?b）/（r+g）、ExR、ExGR、r/g、CIVE與氮營養指數呈現極顯著負相關（P＜0.01），其中極顯著負相關最大的圖像指數是r/b，其值是?0.855 2；圖像指數b、g+b、GRVI、MGRVI、NDI、VARI、WI、GLA、ExGR與氮營養指數呈現極顯著正相關（P＜0.01），其中極顯著正相關最大的圖像指數是b，其值為0.855 7。按照相關性絕對值從大到小排列的前10個圖像指數依次為：b、r/b、IKAW、r?b、（r?g?b）/（r+g）、r、g+b、GLA、g/b、VARI。在開花期，圖像指數r、g/b、r?b、IKAW、（r?g?b）/（r+g）、ExR、NDI、r/g、r+b、CIVE與氮營養指數呈現極顯著負相關（P＜0.01），其中極顯著負相關最大的圖像指數是IKAW，其值是?0.818 4；圖像指數g、b、g+b、ExG、GRVI、MGRVI、VARI、EXGR、GLA、GLI、ExGR與氮營養指數呈現極顯著正相關（P＜0.01），其中極顯著正相關的圖像指數與挑旗期的圖像指數一樣，都是b，其值為0.783 4。按照相關性絕對值從大到小排列的前10個圖像指數依次為：IKAW、r/b、r?b、b、（r?g?b）/（r+g）、g+b、r、GLA、V ARI、ExR。

表 2 不同生育期圖像指數與氮營養指數的相關性Table 2 Correlation coefficients between image indices and NNIs at two growth stages of winter wheat crop

2.1.2 圖像指數間多重共線性篩選 為防止偏最小二乘回歸模型入選參量之間的高相關，采用方差膨脹因子對圖像指數進行多重共線性分析，如圖2所示。從圖2可以看出，當0≤VIF＜10時，圖像指數間不存在多重共線性，如g/b與GRVI 、g與r、g?b與r等；當10≤VIF≤20時，圖像指數間存在一定的共 線 性，如GRVI與g+b、（r?g?b）/（r+g）與ExR等；當VIF＞20時，圖像指數間存在嚴重的共線性，如r/b與r?b、ExR與VARI等。為了有效降低模型多重共線性，同時保證對氮營養指數反演敏感的圖像指數，將相關系數與方差膨脹因子進行整合，選擇圖像指數與氮營養指數間相關性強、圖像指數間共線性弱的圖像指數進行建模。結合表2和圖2，選取挑旗期圖像指數b、g/b、（r?g?b）/（r+g）、NDI、WI，開花期圖像指數b、r/b、（r?g?b）/（r+g）、V ARI作為多元變量入選偏最小二乘回歸估算模型。

2.2 冬小麥氮營養指數的估算及填圖

結合相關系數與方差膨脹因子獲得挑旗期和開花期對氮營養指數較為敏感的圖像指數，通過偏最小二乘回歸建立挑旗期和開花期的氮營養指數反演模型。估算和驗證結果見表3和圖3。挑旗期的估算模型為：

開花期估算模型為：

從表3和圖3可以看出，挑旗期圖像指數反演氮營養指數的 R2為0.711 6，RMSE為0.107 5，MAE為0.086 0；驗證模型的RMSE為0.135 8，MAE為0.105 2，實測值與預測值分布在1∶1線附近，表明該模型具有較高的精度和可靠性。開花期圖像指數反演氮營養指數的 R2為0.720 4，RMSE為0.085 8，MAE為0.071 7；驗證模型的RMSE為0.187 1，MAE為0.131 0，實測值與預測值分布在1∶1線附近偏下，一些氮營養指數值被嚴重低估，建模效果較好，但驗證效果較差。

圖 2 各生育期圖像指數之間的方差膨脹因子Fig. 2 Variance inflation factor of image indices at flagging and flowering stages

表 3 建模與驗證Table 3 Modeling and validation

結合相關性和方差膨脹因子篩選的變量，采用偏最小二乘回歸方法構建冬小麥氮營養指數模型，對模型進行填圖及可視化，結果如圖4所示。從圖中可以看出，顏色由藍色過渡到紅色表示冬小麥氮素營養狀況由氮肥不足到氮肥過量，小區顏色越偏向藍色，表明該小區氮肥越不足，越接近紅色，表明該小區氮肥越過量。結合圖1分析可知，在圖4a冬小麥挑旗期，建立的冬小麥氮營養診斷模型總體上能反映該生育期冬小麥氮素營養狀況，有一部分試驗區被低估，但與實測值相差較小。在圖4b冬小麥開花期，建立的冬小麥氮營養診斷模型能夠較為準確地反映冬小麥氮素營養狀況，只有極少部分試驗區內作物氮素被低估。總之，采用遙感監測圖可以直觀地反映冬小麥氮素營養空間分布狀況，為缺氮小區或過氮小區及時調整氮肥施加量提供了有效依據。

圖 3 冬小麥各生育期NNI實測值與偏最小二乘模型估測值關系Fig. 3 Correlation between measured and estimated values of NNIs of winter wheat at each growth stage by partial least squares model

圖 4 冬小麥各生育期NNI預測值Fig. 4 Estimated values of NNIs of winter wheat at each growth stage

3 討論

氮營養指數是一個可以相對準確地判定作物氮營養狀況的指標，當氮營養指數小于1時，處于缺氮狀態，作物生長緩慢，莖生長細且很短；當氮營養指數大于1時，處于過氮狀態，作物貪青徒長，推遲成熟期；當氮營養指數等于1時，作物植株氮素營養適量。本研究首先結合相關性和方差膨脹因子篩選了對氮營養指數敏感的圖像指數，然后采用偏最小二乘法構建了挑旗期和開花期的氮營養指數模型。在挑旗期利用篩選的圖像指數預測氮營養指數變化范圍為0.50～1.04，在相同施氮水平下預測的氮營養指數值與實測氮營養指數值接近，在開花期利用篩選的圖像指數預測氮營養指數變化范圍為0.45～0.96，在相同施氮水平下預測的氮營養指數一部分與實測的氮營養指數值相近，一部分被低估。魏鵬飛等[39]用多光譜傳感器獲取多光譜影像反演夏玉米葉片氮含量狀況，其模型的決定系數的范圍為0.56～0.64，均方根誤差為20.50%～31.12%。陳鵬飛等[40]用多光譜傳感器獲取多光譜影像反演棉花氮素營養狀況，其模型的決定系數范圍為0.19～0.57，均方根誤差的范圍為0.17%～0.23%。劉帥兵等[41]利用數碼影像反演冬小麥氮含量，其模型的決定系數范圍0.689 7～0.815 5，均方根誤差為0.106 1～0.375 8。秦占飛等[42]通過UHD185機載成像光譜儀獲取高光譜影像反演氮素營養狀況，其模型的決定系數是0.673，均方根誤差是0.329。本研究使用的數碼相機傳感器所能提供的光譜信息很少，通過結合相關系數-方差膨脹因子，其估算作物氮素營養狀況的效果與這些研究相當或高于上述研究結果，為冬小麥氮素營養診斷反演提供了一種有效的思路。

另外，本研究僅利用了冬小麥挑旗期和開花期的數碼影像及地面實測數據，而且是分生育期進行建模，對構建的模型的通用性有一定的影響，應增加不同年限及不同區域的數據，以提高模型的穩定性和普適性。另外，數碼影像僅包含紅綠藍三個波段的信息，接下來的研究是提取數碼影像的圖像頻率，利用圖像頻率結合圖像指數對作物氮素營養狀況進行研究是今后努力的方向。

4 結論

本研究從數碼影像紅綠藍三個波段形成的圖像指數著手，結合相關性和方差膨脹因子篩選對氮營養指數敏感的圖像指數，采用偏最小二乘回歸方法探討篩選的敏感指數估算氮營養指數的能力，并利用該模型對數碼影像進行填圖和空間可視化分析，結果表明：

（1）結合相關系數-方差膨脹因子，篩選出與氮營養指數相關性好、指標與指標間多重共線性小的圖像指數。其中，挑旗期篩選的圖像指數是b、g/b、（r?g?b）/（r+g）、NDI、WI，開花期篩選的圖像指數是b、r/b、（r?g?b）/（r+g）、VARI。

（2）將挑旗期和開花期構建的氮營養指數模型填圖于無人機數碼影像中，兩個生育期都能反映出各自生育期的氮素營養狀況，但開花期氮營養指數模型能更加準確地反映冬小麥氮素營養狀況。本研究結果可為無人機數碼影像在精準施肥中的應用提供理論依據。