程 千 徐洪剛 曹引波 段福義 陳 震

(中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)田灌溉研究所，新鄉(xiāng) 453002)

0 引言

收獲前，及時(shí)準(zhǔn)確地估測(cè)農(nóng)作物產(chǎn)量有助于相關(guān)糧食政策的調(diào)整及制定。遙感技術(shù)是監(jiān)測(cè)作物生長(zhǎng)狀況的有效手段，遙感估產(chǎn)則是將遙感信息作為輸入變量，通過(guò)建立遙感估產(chǎn)模型來(lái)估測(cè)作物產(chǎn)量。隨著科技的發(fā)展，無(wú)人機(jī)能夠以低成本和更實(shí)用的方式獲取高時(shí)空分辨率的遙感數(shù)據(jù)，并將其應(yīng)用于作物生長(zhǎng)監(jiān)測(cè)[1-6]和產(chǎn)量估測(cè)中[7-9]，彌補(bǔ)了衛(wèi)星遙感運(yùn)行周期長(zhǎng)、空間分辨率低和受氣象條件影響大等不足。

多光譜傳感器包含的近紅外波段能夠較好地捕捉作物整個(gè)生育期內(nèi)葉片結(jié)構(gòu)變化產(chǎn)生的光譜反射特征，其與其他波段構(gòu)建的植被指數(shù)包含更多的作物生長(zhǎng)信息，有助于產(chǎn)量估測(cè)。通過(guò)建立植被指數(shù)與產(chǎn)量間的統(tǒng)計(jì)關(guān)系進(jìn)行遙感估產(chǎn)的方法，簡(jiǎn)單易行，目前已有很多學(xué)者利用無(wú)人機(jī)平臺(tái)搭載多光譜傳感器對(duì)不同作物進(jìn)行產(chǎn)量估測(cè)[6,10-13]。

在冬小麥生長(zhǎng)階段，冠層反射率在不斷發(fā)生變化，冬小麥生殖生長(zhǎng)階段的光譜指數(shù)與產(chǎn)量具有更高的相關(guān)性，可更好地用來(lái)估測(cè)產(chǎn)量，如開(kāi)花期[14]或灌漿初期[15]的光譜植被指數(shù)等。隨著生長(zhǎng)過(guò)程的延續(xù)，產(chǎn)量與單一時(shí)期光譜指數(shù)之間的相關(guān)性達(dá)到非常顯著的水平，但是無(wú)法在單個(gè)時(shí)期內(nèi)確定產(chǎn)量形成的過(guò)程。相反，由多個(gè)生育期生長(zhǎng)信息構(gòu)成的數(shù)據(jù)集能更好地反映作物生長(zhǎng)變化過(guò)程，可以提供更多有用的信息，有利于提高算法的估測(cè)精度。文獻(xiàn)[16-17]利用多個(gè)生育期的植被指數(shù)的累積值進(jìn)行作物產(chǎn)量估測(cè)，提高了估測(cè)精度。本文嘗試將新生育期的植被指數(shù)作為樣本數(shù)據(jù)新的屬性值添加到上一個(gè)生育期的數(shù)據(jù)集中，形成新的多時(shí)相植被指數(shù)數(shù)據(jù)集，從而對(duì)冬小麥產(chǎn)量進(jìn)行估測(cè)。

本文采用無(wú)人機(jī)搭載多光譜傳感器，對(duì)冬小麥多個(gè)生育期進(jìn)行遙感監(jiān)測(cè)，提取多種光譜植被指數(shù)，并利用多時(shí)相的植被指數(shù)構(gòu)建數(shù)據(jù)集，訓(xùn)練機(jī)器學(xué)習(xí)算法，進(jìn)行冬小麥產(chǎn)量的估測(cè)分析。

1 材料與方法

1.1 試驗(yàn)區(qū)概況

試驗(yàn)區(qū)位于河南省新鄉(xiāng)市七里營(yíng)鎮(zhèn)中國(guó)農(nóng)業(yè)科學(xué)院新鄉(xiāng)綜合試驗(yàn)基地(113°45′38″N，35°8′10″E)，該試驗(yàn)區(qū)地形為黃河北部沖積扇平原，海拔73 m，土壤類(lèi)型為粘壤土。氣候類(lèi)型屬于溫帶季風(fēng)氣候，四季分明，雨熱同期，多年年平均氣溫14.5℃，年平均降水量560 mm左右，主要集中在6—9月，平均無(wú)霜期199 d，試驗(yàn)研究區(qū)位置見(jiàn)圖1。本次研究共有3種灌溉水量處理(Irrigation treatment，IT)：240 mm(IT1)、190 mm(IT2)和145 mm(IT3)，每種灌水處理包含60個(gè)小區(qū)；試驗(yàn)選取了30個(gè)小麥品種，每個(gè)小區(qū)1個(gè)品種，隨機(jī)分布種植，則每個(gè)灌水處理對(duì)應(yīng)30個(gè)品種，重復(fù)兩次，試驗(yàn)設(shè)計(jì)如圖1所示。

1.2 無(wú)人機(jī)平臺(tái)與遙感數(shù)據(jù)獲取

采用大疆經(jīng)緯M210型四旋翼無(wú)人機(jī)平臺(tái)，無(wú)人機(jī)凈質(zhì)量4.69 kg，最大承載質(zhì)量1.45 kg,續(xù)航時(shí)間20 min左右。無(wú)人機(jī)搭載MicaSense RedEdge-MX型多光譜相機(jī)拍攝多光譜影像。多光譜相機(jī)質(zhì)量232 g,圖像尺寸1 280像素×960像素,具有藍(lán)(Blue，B，(475±20)nm)、綠(Green，G，(560±20)nm)、紅(Red，R，(668±10)nm)、紅邊(Red-edge，RE，(717±10)nm)和近紅外(Near infrared，NIR，(840±40)nm)共5個(gè)通道。

冬小麥抽穗之后，進(jìn)入生殖生長(zhǎng)階段，各器官營(yíng)養(yǎng)逐步轉(zhuǎn)化為籽粒的形成，小麥葉片開(kāi)始衰落，此時(shí)雖然冠層光譜反射率不斷降低，但與產(chǎn)量的相關(guān)性逐步升高[12]。因此，從冬小麥開(kāi)花期到灌漿期，選取多次采集的多光譜影像，進(jìn)行產(chǎn)量預(yù)測(cè)。具體采集日期為：開(kāi)花期(S1，2020年4月14日)、灌漿初期(S2，2020年4月23日)、灌漿中期(S3，2020年4月30日)和灌漿末期(S4，2020年5月10日)。

影像采集時(shí)使用大疆GSPro地面站對(duì)試驗(yàn)區(qū)覆蓋范圍規(guī)劃航線(xiàn)自主飛行。為保證影像質(zhì)量，減少大氣影響，選擇天氣晴好、無(wú)云低風(fēng)速且11:00—14:00時(shí)間段飛行。飛行高度30 m，影像重疊率為航向85%和旁向80%,地面空間分辨率為2 cm。

1.3 數(shù)據(jù)處理

1.3.1影像預(yù)處理

采用Pix4Dmapper軟件對(duì)無(wú)人機(jī)的多光譜影像進(jìn)行拼接處理，并導(dǎo)入每次飛行前后拍攝的校準(zhǔn)板的影像進(jìn)行輻射校正，生成該試驗(yàn)區(qū)地物反射率的正射影像圖。將正射影像圖導(dǎo)入開(kāi)源軟件QGIS，進(jìn)行幾何校正、試驗(yàn)小區(qū)矢量文件制作、多光譜植被指數(shù)計(jì)算、各小區(qū)多光譜植被指數(shù)統(tǒng)計(jì)等，得到各小區(qū)多光譜植被指數(shù)均值作為產(chǎn)量估測(cè)的輸入變量。

1.3.2多光譜植被指數(shù)選取

目前已有多種多光譜植被指數(shù)用于作物產(chǎn)量估測(cè)研究，如表1所示，大致分為以下幾類(lèi)[18]：①基于紅波段和近紅外波段運(yùn)算的多光譜指數(shù)：差異植被指數(shù)(Difference vegetation index，DVI)、增強(qiáng)植被指數(shù)(Enhanced vegetation index，EVI)、比率植被指數(shù)(Ratio vegetation index，RVI)、改善簡(jiǎn)單比率指數(shù)(Modified simple ratio，MSR)、歸一化植被指數(shù)(Normalized difference vegetation index，NDVI)、植被色素比率(Plant pigment ratio，PPR)、結(jié)構(gòu)不敏感色素指數(shù)(Structure insensitive pigment index，SIPI)和改進(jìn)的非線(xiàn)性植被指數(shù)(Modified nonlinear vegetation index，MNVI)。②優(yōu)化土壤背景的紅波段和近紅外波段運(yùn)算的多光譜指數(shù)：土壤調(diào)節(jié)植被指數(shù)(Soil-adjusted vegetation index，SAVI)和優(yōu)化土壤調(diào)節(jié)植被指數(shù)(Optimization of soil-adjusted vegetation index，OSAVI)。③利用紅邊波段特征的多光譜指數(shù)：紅邊葉綠素指數(shù)(Red-edge chlorophyll index，CIRE)、紅邊歸一化植被指數(shù)(Red-edge NDVI，NDVIRE)、改良葉綠素吸收率指數(shù)(Modified chlorophyll absorption ratio index，MCARI)和轉(zhuǎn)化葉綠素吸收反射指數(shù)(Transformed chlorophyll absorption reflectance index，TCARI)。④基于從綠波段到近紅外波段光譜反射率曲線(xiàn)形狀的指數(shù)：三角植被指數(shù)(Triangular vegetation index，TVI)和改善三角植被指數(shù)2(Modified triangular vegetation index 2，MTVI2)等。

表1 植被指數(shù)計(jì)算公式Tab.1 Vegetation index calculation formula

1.4 數(shù)據(jù)分析方法

水分虧缺條件下冬小麥生長(zhǎng)特征出現(xiàn)差異，遙感影像的像元灰度變化能夠反映生長(zhǎng)特征的變化。通過(guò)多光譜植被指數(shù)評(píng)價(jià)水分虧缺條件下，冬小麥生育期的生長(zhǎng)差異，并采用皮爾遜相關(guān)系數(shù)和決定系數(shù)評(píng)價(jià)不同生育期的多光譜植被指數(shù)與冬小麥產(chǎn)量的線(xiàn)性相關(guān)程度,選擇合適的生育期和植被指數(shù)進(jìn)行產(chǎn)量估測(cè)。

分別采用偏最小二乘回歸(Partial least squares regression, PLSR)、支持向量機(jī)回歸(Support vector regression, SVR)、隨機(jī)森林回歸(Random forest regression, RFR)等統(tǒng)計(jì)學(xué)習(xí)算法進(jìn)行冬小麥產(chǎn)量估測(cè)；并采用決定系數(shù)R2、均方根誤差(Root mean squared error, RMSE)和平均絕對(duì)誤差(Mean absolute error, MAE)等指標(biāo)評(píng)價(jià)不同統(tǒng)計(jì)學(xué)習(xí)算法的估測(cè)精度。以上數(shù)據(jù)分析均使用Python語(yǔ)言中的Scikit-learn[34]機(jī)器學(xué)習(xí)包完成。

2 結(jié)果與分析

2.1 植被指數(shù)與產(chǎn)量相關(guān)性分析

不同生育期多光譜植被指數(shù)與冬小麥產(chǎn)量的相關(guān)系數(shù)如圖2所示，可以看出，進(jìn)入開(kāi)花期后，光譜指數(shù)與產(chǎn)量出現(xiàn)相關(guān)性，隨著冬小麥生長(zhǎng)，相關(guān)性逐漸增強(qiáng)，灌漿末期多個(gè)植被指數(shù)與產(chǎn)量相關(guān)系數(shù)達(dá)到0.7。灌漿期后大部分植被指數(shù)與產(chǎn)量呈正相關(guān)關(guān)系，且植被指數(shù)之間相關(guān)性較強(qiáng)，而TCARI和SIPI與產(chǎn)量及其他植被指數(shù)呈較強(qiáng)的負(fù)相關(guān)關(guān)系。

不同生育期多光譜植被指數(shù)與冬小麥產(chǎn)量的線(xiàn)性回歸決定系數(shù)如表2所示，可以看出，隨著冬小麥生殖階段不斷發(fā)展，決定系數(shù)不斷增大，在灌漿末期達(dá)到最大。除PPR外，其他植被指數(shù)均能在一定程度上反映產(chǎn)量的變化特征?？梢?jiàn)本文選取的多光譜植被指數(shù)與冬小麥產(chǎn)量有較強(qiáng)的相關(guān)性，能夠用來(lái)估測(cè)冬小麥產(chǎn)量。

表2 不同生育期植被指數(shù)與產(chǎn)量決定系數(shù)Tab.2 Coefficient of determination between vegetation index and yield in different growth stages

2.2 多時(shí)相植被指數(shù)估測(cè)產(chǎn)量

選取不同生育期多光譜植被指數(shù)，采用PLSR、SVR和RFR等統(tǒng)計(jì)學(xué)習(xí)算法估測(cè)冬小麥產(chǎn)量，結(jié)果如圖3和表3所示。單時(shí)相時(shí)(如S1或S2等)，將文中選取的16種植被指數(shù)分別作為自變量的1個(gè)維度，形成自變量向量，即每個(gè)樣本數(shù)據(jù)有16個(gè)屬性值。當(dāng)構(gòu)建多時(shí)相植被指數(shù)數(shù)據(jù)集時(shí)(如S1+S2等)，將新生育期(S2)的16個(gè)植被指數(shù)合并到前1個(gè)生育期(S1)的自變量向量中，形成新的多時(shí)相植被指數(shù)數(shù)據(jù)集，每增加1個(gè)時(shí)相，自變量向量增加16個(gè)維度，即每個(gè)樣本數(shù)據(jù)增加16個(gè)屬性值。

在訓(xùn)練統(tǒng)計(jì)學(xué)習(xí)算法時(shí)，為綜合提高算法精度和泛化能力，常采用交叉驗(yàn)證的方法，由于試驗(yàn)有3種灌水處理，為減少灌水處理對(duì)算法的影響，保證數(shù)據(jù)分布的均勻性，本文采用10次5折交叉驗(yàn)證，取均值作為每種算法的評(píng)測(cè)指標(biāo)值，圖3和表3中各評(píng)測(cè)指標(biāo)值針對(duì)測(cè)試集的數(shù)據(jù)，更能代表模型的估測(cè)能力。

由圖3和表3可知，當(dāng)使用單個(gè)生育期的植被指數(shù)估測(cè)冬小麥產(chǎn)量時(shí)，由開(kāi)花期到灌漿初期、灌漿中期和灌漿末期，不同算法的預(yù)測(cè)精度不斷提高，這與植被指數(shù)與產(chǎn)量相關(guān)性變化規(guī)律一致。

表3 不同回歸模型估測(cè)產(chǎn)量統(tǒng)計(jì)分析Tab.3 Validation statistics of different regression models for grain yield prediction

使用生育期(S1+S2)的多時(shí)相植被指數(shù)比單個(gè)時(shí)相(S2)的植被指數(shù)估測(cè)產(chǎn)量的精度高，PLSR模型估測(cè)精度R2提高約0.021；SVR模型R2提高約0.015；RFR模型R2提高約0.051。同理，生育期(S1+S2+S3)的多時(shí)相植被指數(shù)比單個(gè)時(shí)相(S1或S2或S3)的植被指數(shù)估測(cè)產(chǎn)量的精度高，說(shuō)明在灌漿末期(S4)到來(lái)前，這種數(shù)據(jù)集構(gòu)建方式能夠提高產(chǎn)量估測(cè)精度，對(duì)產(chǎn)量作出較好的提前估測(cè)。灌漿末期(S4)的植被指數(shù)與產(chǎn)量相關(guān)性更好，估測(cè)精度更高，此時(shí)構(gòu)建的多時(shí)相植被指數(shù)，不同模型均有較高的估測(cè)精度。

對(duì)比3種統(tǒng)計(jì)學(xué)習(xí)算法，針對(duì)本文的數(shù)據(jù)集，當(dāng)數(shù)據(jù)屬性值增多時(shí)，PLSR估測(cè)能力低于SVR和RFR。PLSR模型針對(duì)生育期(S1+S2+S4)數(shù)據(jù)測(cè)試集的估測(cè)精度最高，平均R2為0.459，RMSE為1 822.746 kg/hm2；SVR模型針對(duì)生育期S4數(shù)據(jù)測(cè)試集的估測(cè)精度最高，平均R2為0.540，RMSE為1 676.520 kg/hm2；RFR模型針對(duì)生育期(S3+S4)數(shù)據(jù)測(cè)試集的估測(cè)精度最高，平均R2為0.560，RMSE為1 633.896 kg/hm2。而且，RFR模型針對(duì)生育期S4后不同時(shí)相植被指數(shù)組合的數(shù)據(jù)測(cè)試集，均表現(xiàn)出較高的估測(cè)精度。

圖4為不同灌水處理時(shí)，不同品種冬小麥產(chǎn)量估測(cè)值與實(shí)測(cè)值對(duì)比，圖中估測(cè)值針對(duì)某一生育期的全部數(shù)據(jù)樣本，且估測(cè)模型分別選取估測(cè)精度最高的訓(xùn)練模型，PLSR模型針對(duì)生育期(S1+S2+S4)，SVR模型針對(duì)生育期S4，RFR模型針對(duì)生育期(S3+S4)。可以看出，3種模型的估測(cè)值均反映不同灌溉處理?xiàng)l件下產(chǎn)量變化規(guī)律，并且RFR模型的估測(cè)值與實(shí)測(cè)值更為接近。

正常灌溉處理(IT1)時(shí)，冬小麥產(chǎn)量平均達(dá)到8 700 kg/hm2，水分虧缺處理(IT2和IT3)時(shí)，產(chǎn)量逐步降低，且水分虧缺嚴(yán)重時(shí)，產(chǎn)量下降較多，相比正常灌溉，冬小麥減產(chǎn)將近1/2。對(duì)比冬小麥產(chǎn)量估測(cè)值與實(shí)測(cè)值，IT1灌水處理時(shí)，估測(cè)值偏??；IT2和IT3灌水處理時(shí)，估測(cè)值與實(shí)測(cè)值較為接近；但每個(gè)灌水處理?xiàng)l件下估測(cè)值的分布范圍均比實(shí)測(cè)值的小。由于試驗(yàn)區(qū)冬小麥品種較多，同一灌水處理?xiàng)l件下，特別是正常灌水處理時(shí)，產(chǎn)量相差較大，使得數(shù)據(jù)分布較為分散。

3 討論

本文將新生育期的植被指數(shù)作為新的屬性值，添加到上一個(gè)生育期的數(shù)據(jù)集中，形成新的多時(shí)相植被指數(shù)數(shù)據(jù)集，提高了冬小麥產(chǎn)量估測(cè)精度，這與文獻(xiàn)[11,13,17]的結(jié)論一致。而且文獻(xiàn)[11]研究中同時(shí)使用了株高參數(shù)，并發(fā)現(xiàn)灌漿期株高比植被指數(shù)對(duì)冬小麥產(chǎn)量有更高的估測(cè)精度，說(shuō)明不同類(lèi)型傳感器獲取的冠層參數(shù)，經(jīng)過(guò)數(shù)據(jù)融合能進(jìn)一步提高作物產(chǎn)量估測(cè)精度?？梢?jiàn)，無(wú)論是多時(shí)相數(shù)據(jù)或多源數(shù)據(jù)，均提高了樣本數(shù)據(jù)屬性值的維度，有助于產(chǎn)量估測(cè)模型的訓(xùn)練。

采用統(tǒng)計(jì)學(xué)習(xí)模型進(jìn)行回歸預(yù)測(cè)時(shí)，樣本數(shù)據(jù)的選擇對(duì)模型的訓(xùn)練效果有較大影響，特別是樣本數(shù)據(jù)量較小時(shí)，訓(xùn)練數(shù)據(jù)集的劃分更為重要。交叉驗(yàn)證法適用于樣本數(shù)量不大的情況，而且對(duì)模型偏差和方差能夠兼顧。本文采用P次K折交叉驗(yàn)證法，減小訓(xùn)練數(shù)據(jù)隨機(jī)選取對(duì)模型預(yù)測(cè)精度帶來(lái)的影響。P次為重復(fù)K折交叉驗(yàn)證的次數(shù)，通過(guò)P次重復(fù)減小隨機(jī)誤差，本文P選取10。K折即把數(shù)據(jù)集劃分為K等份，其中K-1份作為訓(xùn)練集，剩余1份作為測(cè)試集，由于樣本總量為180個(gè)，單生育期時(shí)數(shù)據(jù)屬性為16個(gè)，多時(shí)相時(shí)數(shù)據(jù)屬性最多有64個(gè)，為保證不同生育期訓(xùn)練模型的偏差和方差均較低，K值不宜過(guò)大或過(guò)小，本文K選取5作為數(shù)據(jù)集劃分標(biāo)準(zhǔn)。

冬小麥產(chǎn)量構(gòu)成因素較多，與多光譜植被指數(shù)之間存在較為復(fù)雜的非線(xiàn)性關(guān)系，PLSR模型雖然包含了主成分分析過(guò)程，可以有效解決自變量之間的多個(gè)共線(xiàn)性問(wèn)題，但當(dāng)自變量與因變量之間存在復(fù)雜的關(guān)系時(shí)，其仍有較大的局限性。SVR模型和RFR模型對(duì)復(fù)雜的非線(xiàn)性關(guān)系均有一定的回歸能力，并且RFR模型具有更好的抗噪能力，本文冬小麥試驗(yàn)中有30個(gè)小麥品種，數(shù)據(jù)分布較為離散，針對(duì)灌漿末期的不同時(shí)相植被指數(shù)，RFR模型均表現(xiàn)出較高的估測(cè)精度。但是當(dāng)樣本數(shù)據(jù)量過(guò)大，或者樣本數(shù)據(jù)特征過(guò)多時(shí)，傳統(tǒng)統(tǒng)計(jì)學(xué)習(xí)算法的預(yù)測(cè)能力將下降，需要訓(xùn)練更優(yōu)秀的模型，如深度神經(jīng)網(wǎng)絡(luò)模型等。

4 結(jié)論

(1)隨著冬小麥的不斷生長(zhǎng)，植被指數(shù)與產(chǎn)量的相關(guān)性不斷增強(qiáng)，灌漿末期多個(gè)植被指數(shù)與產(chǎn)量相關(guān)系數(shù)達(dá)到0.7，基于此生育期植被指數(shù)估測(cè)冬小麥產(chǎn)量的精度也較高。

(2)采用多時(shí)相植被指數(shù)能進(jìn)一步提高冬小麥產(chǎn)量的估測(cè)精度，采用冬小麥開(kāi)花期和灌漿初期的多時(shí)相植被指數(shù)進(jìn)行估產(chǎn)比采用單個(gè)生育期的植被指數(shù)估測(cè)產(chǎn)量的精度高，PLSR模型估測(cè)決定系數(shù)R2提高約0.021，SVR模型R2提高約0.015，RFR模型R2提高約0.051。灌漿末期的多時(shí)相植被指數(shù)，3種模型均有較高的估測(cè)精度。PLSR模型估測(cè)精度最高時(shí)R2、RMSE分別為0.459、1 822.746 kg/hm2，SVR模型估測(cè)精度最高時(shí)R2、RMSE分別為0.540、1 676.520 kg/hm2，RFR模型估測(cè)精度最高時(shí)R2、RMSE分別為0.560、1 633.896 kg/hm2，相比之下，RFR模型具有更好的穩(wěn)定性。

(3)長(zhǎng)時(shí)間的水分虧缺導(dǎo)致冬小麥產(chǎn)量下降，與正常灌溉水量相比，長(zhǎng)時(shí)間的水分虧缺將造成冬小麥減產(chǎn)約1/2，此時(shí)需及時(shí)制定灌溉計(jì)劃。