張靜，郭思?jí)簦n迎春，雷亞平，邢芳芳，杜文麗，李亞兵，*，馮璐，*

（1.鄭州大學(xué)農(nóng)學(xué)院，棉花生物學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室鄭州大學(xué)研究基地，鄭州 450000；2.中國(guó)農(nóng)業(yè)科學(xué)院棉花研究所，棉花生物學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，河南 安陽(yáng) 455000）

棉花是我國(guó)最重要的經(jīng)濟(jì)作物之一，棉花產(chǎn)量的快速準(zhǔn)確估算對(duì)其精準(zhǔn)管理和生產(chǎn)決策十分重要。傳統(tǒng)的作物產(chǎn)量估算主要依靠田間測(cè)產(chǎn)調(diào)查，不僅費(fèi)時(shí)費(fèi)工，而且準(zhǔn)確性較低[1]。近年來(lái)，遙感技術(shù)在農(nóng)業(yè)生產(chǎn)監(jiān)測(cè)中應(yīng)用越來(lái)越廣泛。衛(wèi)星遙感技術(shù)由于易受天氣影響不能滿足田塊尺度監(jiān)測(cè)的需要[2-3]，而低空無(wú)人機(jī)遙感技術(shù)具有便于攜帶、低成本、高分辨等優(yōu)點(diǎn)[1,4]，彌補(bǔ)了這一缺點(diǎn)，已成為精準(zhǔn)農(nóng)業(yè)研究的熱門(mén)方向。

目前，國(guó)內(nèi)外學(xué)者基于無(wú)人機(jī)搭載不同圖像采集器（可見(jiàn)光相機(jī)、熱紅外相機(jī)、多光譜相機(jī)）對(duì)作物表型信息開(kāi)展了大量研究[5-6]，實(shí)現(xiàn)了對(duì)作物生長(zhǎng)參數(shù)，如株高[7]、生物量[8]、葉面積指數(shù)[9-10]等的有效評(píng)估。近年來(lái)，作物產(chǎn)量估算成為無(wú)人機(jī)生產(chǎn)監(jiān)測(cè)的重要研究?jī)?nèi)容。王來(lái)剛等[11]基于多光譜圖像提取植被指數(shù)建立玉米估產(chǎn)模型，發(fā)現(xiàn)與單生育期相比，多生育期估算效果更好。申洋洋等[12]基于無(wú)人機(jī)的多光譜圖像利用結(jié)合機(jī)器學(xué)習(xí)算法提取植被指數(shù)，建立了小麥產(chǎn)量估算模型，在抽穗期估算效果最好。此外，基于可見(jiàn)光相機(jī)采集的RGB圖像提取顏色、紋理等[13]特征，在小麥產(chǎn)量估測(cè)[14-15]中也得到了較好的驗(yàn)證，且多特征融合比單一特征估算精度高；多類(lèi)型圖像特征融合（如多光譜和熱紅外）較單一類(lèi)型圖像特征估算精度高[16]。因此，無(wú)人機(jī)遙感技術(shù)在作物產(chǎn)量估算中具有良好的應(yīng)用前景。然而，現(xiàn)有的圖像獲取大多是在同一飛行高度進(jìn)行，研究發(fā)現(xiàn)3 m高度下的水稻估算產(chǎn)量模型優(yōu)于6和9 m[17]，但是不同高度的比較研究仍然較少，此外，基于無(wú)人機(jī)估算作物產(chǎn)量的研究多集中于小麥、玉米等糧食作物，對(duì)棉花的研究相對(duì)較少，且已有關(guān)于棉花產(chǎn)量估算方面的研究多是基于多年氣象和棉花物候信息對(duì)全省范圍的產(chǎn)量建立模型[18-19]，或者是通過(guò)光譜圖像對(duì)中型田塊進(jìn)行棉花產(chǎn)量估算[1,20]，少見(jiàn)基于無(wú)人機(jī)RGB圖像對(duì)田塊尺度的棉花產(chǎn)量估算。

因此，本研究以棉花為研究對(duì)象，通過(guò)無(wú)人機(jī)搭載可見(jiàn)光相機(jī)采集不同生育時(shí)期和不同高度的RGB圖像，為解決單一因素的預(yù)測(cè)精度不高的問(wèn)題[4,6]，提取顏色指數(shù)特征和紋理特征并將其結(jié)合，通過(guò)逐步回歸分析和因子分析2種方法處理數(shù)據(jù)，建立多元回歸模型估算棉花產(chǎn)量，對(duì)比分析不同生育時(shí)期和不同高度建立的棉花產(chǎn)量估算模型，篩選出建模的最佳生育時(shí)期和最優(yōu)高度。旨在為利用無(wú)人機(jī)快速估算棉花產(chǎn)量提供參考。

1 材料與方法

1.1 研究區(qū)域與試驗(yàn)設(shè)計(jì)

本試驗(yàn)于2020年7月至2020年9月在河南省安陽(yáng)市白璧鎮(zhèn)中國(guó)農(nóng)業(yè)科學(xué)院棉花研究所東場(chǎng)試驗(yàn)基地（36°06′N(xiāo),114°21′E）開(kāi)展。供試棉花品種為魯棉研28(SCRC28)。試驗(yàn)田共設(shè)24個(gè)小區(qū)，各小區(qū)面積為64 m2（8 m×8 m），種植行間距0.8 m，共10行。每組6個(gè)小區(qū)，均設(shè)置6個(gè)密度處理（1.5、3.3、5.1、6.9、8.7、10.5萬(wàn)株·hm-2），共4次重復(fù)（重復(fù)Ⅰ~Ⅳ）。棉花田間管理采用常規(guī)高產(chǎn)管理方式。

1.2 棉花產(chǎn)量測(cè)定

在收獲期按小區(qū)實(shí)收測(cè)產(chǎn)，2020年分別在9月28日和11月2日收取2次棉絮，10月29日收獲青鈴，將其相加計(jì)算實(shí)收籽棉產(chǎn)量，公式如下。

式中，Y1為收取棉絮后計(jì)算的產(chǎn)量，kg·hm-2；Y2為收取的青鈴產(chǎn)量，Y為總籽棉產(chǎn)量，W1為籽棉毛重，W2為口袋重，W3為小區(qū)青鈴重，單位均為kg；S指小區(qū)面積，單位為hm2。

1.3 圖像獲取與處理

本研究采用大疆精靈四旋翼無(wú)人機(jī)平臺(tái)，搭載2 048萬(wàn)像素CMOS傳感器，F(xiàn)OV84°鏡頭，自動(dòng)對(duì)焦。分別在棉花花鈴期和吐絮期采集冠層圖像，采集時(shí)間分別為7月10日至8月21日和8月25日至9月25日，各采集10次，圖像采集當(dāng)天時(shí)間為9∶00—10∶00，選擇晴朗無(wú)風(fēng)天氣，陰雨天氣順延。飛行高度分別設(shè)置為20、30和40 m，每隔4 d采集1次，各小區(qū)在每個(gè)高度下各采集1幅圖像，總計(jì)72幅圖像。鏡頭垂直向下，獲取圖像像素大小為5 472×3 078，同步存儲(chǔ)于Micro SD卡，以JEPG格式導(dǎo)入計(jì)算機(jī)中。

借助Photoshop CC2019和PyCharm 2020.2.3軟件對(duì)RGB棉花冠層圖像進(jìn)行預(yù)處理和特征提取。預(yù)處理包括圖像裁剪、分割、灰度化處理等。裁剪是按照小區(qū)劃分裁剪成單一密度圖像，并將裁剪后的圖像保存為JPG格式；通過(guò)HSV閾值分割方法將棉株和背景（土壤、傳感器等）分離；采用視覺(jué)法對(duì)分割后的圖像進(jìn)行灰度化進(jìn)而提取紋理特征。經(jīng)預(yù)處理后，棉花花鈴期和吐絮期各飛行高度分別獲得240幅冠層圖像。

1.4 圖像特征參數(shù)計(jì)算

1.4.1 顏色指數(shù)特征提取圖像有多種顏色空間，無(wú)人機(jī)獲取的原始圖像是基于RGB彩色空間的[21]。本研究基于分割后的棉花冠層圖像提取RGB顏色特征值，為了減少陰影的影響，首先對(duì)其進(jìn)行歸一化處理并計(jì)算出r、g、b值，然后根據(jù)不同組合計(jì)算12個(gè)可見(jiàn)光顏色指數(shù)（color index,CI）[22-26]，分別是可見(jiàn)光大氣阻抗植被指數(shù)（visible atmospherically resistant index,VARI）、超綠植被指數(shù)（excess green vegetation index,EXG）、超紅植被指數(shù)（excess red vegetation index,EXR）、超藍(lán)植被指數(shù)（excess blue vegetation index,EXB）、綠紅差值指數(shù)（excess green minus excess red vegetation index,EXGR）、歸一化差分指數(shù)（normalized greenred difference index,NDI）、改良綠紅植被指數(shù)（modified green red vegetation index,MGRVI）、Woebbecke指數(shù)（woebbecke index,WI）、川島指數(shù)（kawashima index,IKAW）、綠葉指數(shù)（green leaf algorithm,GLI）、紅綠藍(lán)植被指數(shù)（red green blue vegetation index,RGBVI）、植物指數(shù)（vegetativen,VEG）。

式中，r、g、b分別為色度坐標(biāo)，R、G、B分別為各分量的最大值。

1.4.2 紋理特征提取紋理反映了物體本身的屬性，整體上表現(xiàn)出某種規(guī)律性，其灰度分布具有周期性[21]。數(shù)字圖像中的紋理是相鄰像素的空間位置變化的視覺(jué)表現(xiàn)[21]。本研究基于灰度共生矩陣（gray-level co-occurrence matrix,GLCM）[14,27]和Tamura紋理特征[28-29]提取灰度圖像的紋理特征（textural feature,TF）。其中，灰度共生矩陣選擇了4個(gè)常用參數(shù)，分別是能量（energy）、對(duì)比度（contrast）、熵（entropy）、逆差分矩（inverse different moment,IDM）；而Tamura紋理特征選擇的3種常見(jiàn)指標(biāo)是粗糙度（coarseness）、對(duì)比度（contrast）、方向度（directionality）。

1.5 產(chǎn)量估算

本文主要使用逐步回歸分析（stepwise regression,SWR）和因子分析（factor analysis,FN）方法對(duì)RGB圖像特征數(shù)據(jù)進(jìn)行處理并建立多元回歸模型，主要步驟為：①各生育時(shí)期采集的10次圖像提取特征值計(jì)算平均值；②結(jié)合圖像特征參數(shù)（12個(gè)顏色指數(shù)特征、7個(gè)紋理特征），利用逐步回歸分析方法對(duì)所有特征值進(jìn)行篩選，同時(shí)利用因子分析方法處理特征值，得到方差貢獻(xiàn)率高于90%的因子；③分別探討圖像特征參數(shù)與產(chǎn)量建立的模型，確定估算棉花產(chǎn)量的最佳生育時(shí)期和最佳高度；④選擇重復(fù)Ⅰ~Ⅲ為建模數(shù)據(jù)集，對(duì)建立的模型進(jìn)行分析，用未參與建模的重復(fù)Ⅳ樣本驗(yàn)證產(chǎn)量估算模型的精度。

逐步回歸分析原理是依次不斷反復(fù)引入自變量直到將所有顯著變量引入，采用F檢驗(yàn)和t檢驗(yàn)進(jìn)行顯著性檢驗(yàn)，最終將不顯著自變量刪除后得到最優(yōu)解釋變量集[30]；因子回歸是利用降維方法將多變量歸結(jié)為少數(shù)幾個(gè)新變量，即因子的統(tǒng)計(jì)[31-32]。

1.6 評(píng)價(jià)指標(biāo)

本研究采用決定系數(shù)（R2）、平均絕對(duì)百分比誤差（mean absolute percentage error，MAPE）、均方根誤差（root mean square error，RMSE）和標(biāo)準(zhǔn)均方根誤差（normal root mean square error，NRMSE）4個(gè)指標(biāo)來(lái)評(píng)價(jià)預(yù)測(cè)值與實(shí)測(cè)值擬合性。R2越大模型擬合效果越好；MAPE、RMSE越小模型精度越高；NRMSE越小，預(yù)測(cè)值和實(shí)測(cè)值一致性越好。

式中，xi為實(shí)測(cè)的棉花產(chǎn)量，為實(shí)測(cè)的棉花產(chǎn)量平均值；yi為預(yù)測(cè)的棉花產(chǎn)量；n為樣本數(shù)。

2 結(jié)果與分析

2.1 特征值篩選和因子選擇結(jié)果

將顏色指數(shù)特征和紋理特征結(jié)合起來(lái)，分別進(jìn)行逐步回歸分析和因子分析。逐步回歸結(jié)果如表1所示，可以看出，WI、方向度和IKAW、方向度分別是高度20 m花鈴期和吐絮期篩選后的特征值；WI、ENG和WI分別是高度30 m花鈴期和吐絮期的影響因素；而高度40 m花鈴期和吐絮期篩選的特征值分別是IDM和粗糙度。表2中FA1、FA2、FA3分別表示不同高度、不同生育時(shí)期經(jīng)因子分析后貢獻(xiàn)率超90%的3個(gè)因子，經(jīng)計(jì)算得出各因子的系數(shù)值如表2所示。

表1 逐步回歸分析結(jié)果Table 1 Results of stepwise regression analysis

表2 因子分析結(jié)果Table 2 Results of factor analysis

2.2 不同產(chǎn)量模型擬合效果分析

分別將逐步回歸篩選出的特征值以及因子分析得出的因子與產(chǎn)量建立多元回歸模型，對(duì)同一高度下花鈴期和吐絮期的產(chǎn)量估算模型進(jìn)行對(duì)比分析。由表3可知，在20 m高度下，花鈴期的R2比吐絮期高，RMSE、NRMSE比吐絮期低；在30 m高度下，花鈴期的R2均在0.7以上（=0.784 9，=0.757 4），而吐絮期的R2僅為0.415 9和0.505 2，RMSE、NRMSE均比花鈴期低；在40 m高度下，花鈴期和吐絮期的R2都在0.5左右，RMSE、NRMSE值較高，模型擬合效果一般。綜合來(lái)看，20和30 m的高度下，花鈴期產(chǎn)量估算模型擬合度均優(yōu)于吐絮期且擬合效果較好，而40 m高度花鈴期和吐絮期擬合效果一般。

表3 不同產(chǎn)量模型對(duì)比分析Table 3 Comparative analysis of different yield models

2.3 不同生育時(shí)期產(chǎn)量模型驗(yàn)證分析

2.3.1 高度20 m模型分析對(duì)建立的產(chǎn)量估算模型進(jìn)行驗(yàn)證分析，結(jié)果表明，對(duì)于高度20 m（圖1）建立的模型，在花鈴期通過(guò)逐步回歸方法建立的模型驗(yàn)證結(jié)果顯著（P<0.05）而通過(guò)因子分析方法建立的估產(chǎn)模型不顯著（P>0.05）；在吐絮期兩種方法建立的估產(chǎn)模型驗(yàn)證結(jié)果都顯著（P<0.05），其中，通過(guò)因子分析方法的R2更高、擬合效果更好，但通過(guò)逐步回歸方法驗(yàn)證的R2與之僅相差約0.05，而MAPE、RMSE和NRMSE更小，其模型具有更高的精度，與實(shí)測(cè)產(chǎn)量的一致性更好。綜合對(duì)比逐步回歸方法在花鈴期和吐絮期的擬合模型和驗(yàn)證結(jié)果，花鈴期的產(chǎn)量估算模型擬合效果更好，模型精度更高。

圖1 高度20 m花鈴期和吐絮期產(chǎn)量估算模型驗(yàn)證結(jié)果Fig.1 Verification results of flowering and boll development stage and boll-opening stage yield estimation model at the height of 20 m

2.3.2 高度30 m模型分析 在30 m高度（圖2）下觀察到了類(lèi)似的結(jié)果，無(wú)論是花鈴期還是吐絮期，驗(yàn)證結(jié)果均顯著（P<0.05）。所有驗(yàn)證的R2均達(dá)到了0.9以上，通過(guò)逐步回歸驗(yàn)證的4個(gè)評(píng)價(jià)指標(biāo)（R2、MAPE、RMSE和NRMSE）均優(yōu)于因子分析方法的驗(yàn)證結(jié)果；對(duì)比通過(guò)逐步回歸方法在花鈴期和吐絮期的模型驗(yàn)證結(jié)果，花鈴期的R2更高、擬合效果更好，RMSE和NRMSE更小，說(shuō)明其模型具有更高的精度、與實(shí)測(cè)產(chǎn)量的一致性更好。因此，綜合對(duì)比產(chǎn)量估算模型和驗(yàn)證結(jié)果，通過(guò)逐步回歸方法在花鈴期建立的產(chǎn)量估算模型擬合效果更好，模型精度更高。

圖2 高度30 m花鈴期和吐絮期產(chǎn)量估算模型驗(yàn)證結(jié)果Fig.2 Verification results of flowering and boll development stage and boll-opening stage yield estimation model at the height of 30 m

2.3.3 高度40 m模型分析40 m建立的產(chǎn)量估算模型驗(yàn)證結(jié)果如圖3所示，可以看出，吐絮期驗(yàn)證的R2在0.8以上（=0.801 2，=0.946 9），但與擬合模型相差較大（=0.506 0，=0.448 1）（表3），盡管通過(guò)逐步回歸方法驗(yàn)證的R2比因子分析方法低，但MAPE、RMSE和NRMSE均比因子分析小，因此，與通過(guò)因子分析方法預(yù)測(cè)的結(jié)果相比，逐步回歸方法預(yù)測(cè)值更接近實(shí)測(cè)值。

圖3 高度40 m吐絮期產(chǎn)量估算模型驗(yàn)證結(jié)果Fig.3 Verification results of boll-opening stage yield estimation model at the height of 40 m

綜合來(lái)看，3個(gè)高度的2種方法對(duì)比結(jié)果均顯示，通過(guò)逐步回歸方法建立的產(chǎn)量估算模型擬合度和精度優(yōu)于因子分析方法。在高度20 m花鈴期建立的產(chǎn)量估算模型預(yù)測(cè)效果更好、精度更高；高度30 m在花鈴期建立的產(chǎn)量估算模型預(yù)測(cè)效果和精度相當(dāng)；而40 m高度下，與花鈴期相比，在吐絮期建立的產(chǎn)量估算模型效果更好，但擬合效果一般，穩(wěn)定性較差。

2.4 不同高度模型對(duì)比

進(jìn)一步比較通過(guò)SWR方法在20（花鈴期）、30（花鈴期）和40 m（吐絮期）3個(gè)高度下建立的3個(gè)模型并篩選出產(chǎn)量估算模型的最佳高度和最佳生育時(shí)期。由表3可知，高度30 m在花鈴期建立的模型R2（0.784 9）最高，20 m高度下在花鈴期建立的模型的R2（0.733 3）次之，高度40 m在吐絮期的R2（0.506 0）最低，而RMSE和NRMSE的值從小到大依次是高度30、20、40 m。綜合來(lái)看，在花鈴期、高度為30 m建立的產(chǎn)量估算模型擬合度優(yōu)于20和40 m。

分析3個(gè)擬合模型的驗(yàn)證結(jié)果可知（圖1～3），以花鈴期、高度30 m圖像建立的模型驗(yàn)證結(jié)果R2（0.977 6）最高，MAPE（2.31%）、RMSE（93.7 kg·hm-2）和NRMSE（2.81%）最小，與擬合的指標(biāo)相差最小。因此，無(wú)論是從模型擬合度還是驗(yàn)證集的評(píng)估指標(biāo)來(lái)看，高度30 m下的模型擬合度最好、穩(wěn)定性最強(qiáng)、精度最高。

綜合不同生育時(shí)期、不同高度的分析結(jié)果，產(chǎn)量估算模型的最佳生育時(shí)期為花鈴期，最優(yōu)高度為30 m。

3 討論

多項(xiàng)研究表明，作物產(chǎn)量和圖像顏色特征、紋理特征存在相關(guān)性并能夠擬合模型進(jìn)而估測(cè)產(chǎn)量，如龔紅菊[32]驗(yàn)證了基于紋理特征估算水稻產(chǎn)量的可行性。但是通過(guò)無(wú)人機(jī)獲取圖像時(shí)，光照條件不同可能導(dǎo)致基于單一特征指數(shù)的預(yù)測(cè)效果隨機(jī)性較大，將顏色特征和紋理特征結(jié)合，其結(jié)果比單一顏色特征建立的小麥產(chǎn)量預(yù)測(cè)模型精度高[14]。本研究選擇結(jié)合棉花冠層圖像的顏色指數(shù)特征和紋理特征對(duì)產(chǎn)量進(jìn)行估算并得到較穩(wěn)定的模型，與萬(wàn)亮等[6]研究結(jié)果一致。

本研究探討了不同棉花生育時(shí)期建立的產(chǎn)量估算模型，但在不同的生育時(shí)期預(yù)測(cè)效果不同，高度20和30 m下，花鈴期建模效果均優(yōu)于吐絮期，模型擬合程度更高、驗(yàn)證結(jié)果更顯著。劉煥軍等[33]基于時(shí)間序列Landsat影像建立棉花最優(yōu)產(chǎn)量估算模型，發(fā)現(xiàn)其隨著生長(zhǎng)階段表現(xiàn)出差異，得出花鈴期更適合棉花產(chǎn)量估測(cè)的結(jié)果，與本研究結(jié)果一致。而在40 m高度下，吐絮期產(chǎn)量估算模型擬合度優(yōu)于花鈴期，但擬合模型的R2與驗(yàn)證數(shù)據(jù)集相差較大，盡管棉絮占產(chǎn)量的大部分，但吐絮期葉片大量脫落、棉株倒伏[34]等原因可能導(dǎo)致提取特征不精確，也可能與樣本數(shù)量較少有關(guān)。本研究用一年數(shù)據(jù)進(jìn)行建模，而田間試驗(yàn)受氣溫、降水等自然環(huán)境條件、管理措施等影響較大[35]，擬進(jìn)一步增加樣本數(shù)量以及用不同品種、不同年限的棉花數(shù)據(jù)進(jìn)行分析，得到更普遍的產(chǎn)量估算模型。

本文探討了不同高度的棉花產(chǎn)量模型對(duì)比分析，得到30 m的模型預(yù)測(cè)結(jié)果最好。一般來(lái)說(shuō)，不同高度對(duì)試驗(yàn)地的圖像分辨率不同，飛行高度越低精度越高[36]，但本文結(jié)果不同，可能與高度設(shè)置較小有關(guān)，有研究者在10[14]、30[37]、50[1]、80 m[11]都得到了較好的預(yù)測(cè)結(jié)果；也有可能與采集圖像時(shí)的風(fēng)力、光照變化[38]及預(yù)處理時(shí)使用的圖像分割方法相關(guān)。因此，今后可以設(shè)置更多高度梯度以找出更合適的飛行高度進(jìn)一步優(yōu)化棉花產(chǎn)量估算模型，也可以基于聚類(lèi)和深度學(xué)習(xí)[39]等方法對(duì)圖像進(jìn)行分割處理找出更合適的分割算法。本研究結(jié)果表明，通過(guò)提取圖像特征建立產(chǎn)量估算模型是可行的，但仍需要進(jìn)行更多的研究，例如氣候、農(nóng)藝性狀等信息與圖像結(jié)合進(jìn)一步建立更復(fù)雜的模型。