基于無人機遙感影像的多品種水稻抽穗時間估算

莫佳才，彭 漪，劉小娟，王 靖

（武漢大學遙感信息工程學院，湖北 武漢 430079）

水稻是我國三大糧食作物之一，近年來我國每年的水稻總產(chǎn)量穩(wěn)定在2億t以上，約占全國糧食總產(chǎn)量的30%。水稻產(chǎn)量維系著國家糧食安全，提高水稻總產(chǎn)量仍是我國水稻產(chǎn)業(yè)今后的工作重點[1]，而水稻的育種選種是實現(xiàn)這一目標的基礎工作之一。傳統(tǒng)的水稻育種選種工作依賴于水稻工作者人工觀察和記錄，費時費力，而遙感技術(shù)能夠便捷快速地捕捉信息，提高工作效率，在水稻育種栽培上正得到越來越廣泛的應用。水稻生育期監(jiān)測是對水稻各生育時期的生長狀態(tài)信息進行記載的過程，對水稻的田間管理和產(chǎn)量預估來說至關(guān)重要[2]。在移栽稻中，水稻移栽后的生育期可分為返青期、分蘗期、拔節(jié)期、孕穗期、抽穗期、乳熟期和成熟期[3]。其中抽穗期是育種者選擇優(yōu)良親本和進行授粉操作的關(guān)鍵時期，準確地監(jiān)測抽穗期對于水稻育種和產(chǎn)量估算具有重要意義[4]。

研究表明，在水稻生長發(fā)育過程中，其冠層光譜反射率的變化呈現(xiàn)出一定的規(guī)律性，由反射率計算得到的遙感植被指數(shù)（Vegetation index，VI）的規(guī)律性則更為明顯[5]。水稻的植被指數(shù)通常在營養(yǎng)生長期逐漸上升，到抽穗期達到最高值，隨后轉(zhuǎn)為逐漸下降，直至收割[6-8]。正是這種規(guī)律性，使研究者們能夠從水稻植被指數(shù)的時間序列曲線中識別出水稻的生育時期。由于氣象條件、觀測儀器擾動等內(nèi)外因素的影響，觀測到的原始曲線通常存在噪聲，需要進行平滑處理，常用的方法有濾波法、動態(tài)閾值法和函數(shù)擬合法等[9]。平滑后的時序植被指數(shù)曲線存在特征節(jié)點，研究者們常從這些節(jié)點直接提取水稻抽穗期。如WANG 等[10]在浙江省德清縣、楊浩等[11]在湖南省沅江市用國產(chǎn)HJ衛(wèi)星數(shù)據(jù)，從時序植被指數(shù)最大值識別水稻抽穗期，誤差在5 d 以內(nèi)；MA 等[12]以多個研究區(qū)的地面和無人機遙感數(shù)據(jù)為基礎，將擬合的植被指數(shù)時序曲線的曲率最大值點對應的時期認為是水稻的抽穗期，取得了滿意的效果。在小范圍、水稻品種較為單一的情況下，這種方法可以達到一定的精確度，但當水稻品種復雜化時，誤差會變大。如SUN 等[13]用MODIS EVI 時序曲線最大值提取全中國不同地區(qū)站點的水稻抽穗期，絕大多數(shù)結(jié)果與實際相差小于16 d；WANG 等[14]和LI 等[15]分別從MODIS NDVI 序列最大值中識別湖南省和江西省水稻抽穗期，均方根誤差（RMSE）分別為20 d 與10 d。一般而言，水稻抽穗期時長為10～15 d[16]，因而，上述3 個研究試驗提取的抽穗期結(jié)果誤差過大。相關(guān)研究在以下2個方面存在不足。一方面是多基于衛(wèi)星遙感數(shù)據(jù)對水稻抽穗期進行識別。與衛(wèi)星相比，無人機具有方便、靈活和受大氣影響較小的特點，使用無人機拍攝遙感影像可以獲得更高的空間和時間分辨率，更適用于田間尺度的精細農(nóng)業(yè)，因此需要深入研究基于無人機遙感數(shù)據(jù)的水稻抽穗期估算方法。另一方面，多從植被指數(shù)時序曲線的特征點上直接提取水稻抽穗期，忽略了特征點與水稻抽穗期之間可能存在的數(shù)學關(guān)系以及水稻品種的多樣性。因此，在已有研究的基礎上，利用3個獨立試驗區(qū)的1 102個水稻品種的多時相無人機遙感數(shù)據(jù)，構(gòu)建多時相植被指數(shù)曲線，提取關(guān)鍵特征點，探索其與水稻抽穗期的數(shù)學聯(lián)系，以此實現(xiàn)多品種水稻的抽穗期估算，為多品種水稻育種和大范圍農(nóng)業(yè)監(jiān)測提供參考和指導。

1 材料和方法

1.1 試驗區(qū)域

本研究試驗區(qū)域由三部分組成，分別為湖北省鄂州市的試驗區(qū)1 和試驗區(qū)2、海南省陵水縣的試驗區(qū)3。鄂州市的2 個試驗區(qū)（30°22′32″N，114°44′46″E）均位于鄂州市杜山鎮(zhèn)東港村，其中試驗區(qū)1 的試驗田分成1 000 余個面積約1.6 m2的小區(qū)，共種植1 014個水稻品種，于2018年夏季進行試驗；試驗區(qū)2 的試驗田由48 個36 m2的小區(qū)組成，種植了48 個水稻品種，于2019 年夏季進行試驗。試驗區(qū)3（18°31′51″N，110°3′18″E）坐落在陵水縣椰林鎮(zhèn)文官村，包含種植在40 個60 m2小區(qū)的40 個水稻品種，于2018 年冬春季進行試驗。3 個試驗區(qū)的位置和田塊如圖1 所示。3 個試驗區(qū)之間相互獨立且存在明顯差異：鄂州市與陵水縣分別位于亞熱帶與熱帶，氣候差異明顯，試驗期間光溫差異大；種植的水稻品種差異大，僅在試驗區(qū)2 和試驗區(qū)3 之間存在少量相同品種；水稻品種的抽穗時間差異大。

圖1 試驗區(qū)示意圖Fig.1 Location of the experiment

1.2 數(shù)據(jù)獲取與預處理

1.2.1 實際抽穗日期觀測 水稻抽穗是水稻在分蘗拔節(jié)結(jié)束后，完全發(fā)育的稻穗逐漸伸出頂部的現(xiàn)象[17]。當?shù)咎镏腥后w抽穗率達到10%時，為水稻的始穗期[18]。在水稻生長發(fā)育期間，實地持續(xù)觀測水稻長勢情況，判斷并記錄水稻所處生育時期，由此獲得各試驗區(qū)每個品種的實際始穗時間。3個試驗區(qū)的品種、氣象和觀測情況等信息如表1所示。

表1 各試驗區(qū)數(shù)據(jù)概況Tab.1 Data overview of the experiment areas

1.2.2 無人機數(shù)據(jù)采集與預處理 從水稻返青期至成熟期，每隔一定的時間，使用無人機執(zhí)行飛行任務，搭載多光譜相機對試驗區(qū)拍攝遙感影像，獲取各小區(qū)水稻的遙感數(shù)據(jù)。試驗區(qū)1、2、3分別獲取了10、13、15 期的無人機數(shù)據(jù)。本研究使用的是美國Tetracam 公司生產(chǎn)的Micro-MCA 12 通道多光譜相機，內(nèi)含從可見光到近紅外的12 個波段通道，中心波長分別為490、520、550、570、670、680、700、720、800、850、900、950 nm，其中前10個波段的寬度為10 nm，第11、12 個波段的寬度則分別為20、40 nm。為保證影像空間分辨率的一致性與輻射信息的準確性，3 個試驗區(qū)的無人機采用了基本相同的飛行高度且盡量選擇在晴朗無云或少云的天氣狀況下飛行。對拍攝得到的原始影像進行波段間幾何配準，使12 個波段的影像處于同一坐標系下。隨后，以影像拍攝期間鋪設在水稻田旁的平整路面上的標準定標毯為基準，采用線性模型對影像進行輻射定標，將原始影像各波段的像元值轉(zhuǎn)化為光譜反射率。為獲取各品種水稻的反射率，在每個小區(qū)內(nèi)部靠中心位置勾畫合適大小的感興趣區(qū)（Region of interest，ROI），計算ROI 內(nèi)像元的平均反射率，作為該小區(qū)水稻光譜反射率的代表值。

1.2.3 植被指數(shù)提取 選擇了4 種常用植被指數(shù)，分別為歸一化植被指數(shù)（Normalized difference vegetation index，NDVI）、歸一化差分紅邊植被指數(shù)（Normalized difference red edge index，NDRE）、簡單比值植被指數(shù)（Simple ratio index，SR）、二波段增強型植被指數(shù)（Two-band enhanced vegetation index，EVI2），對應計算公式見表2。其中，ρNIR、ρred和ρrededge分別表示近紅外、紅和紅邊波段的反射率，在本研究中對應采用800、670、720 nm波段。

表2 所用植被指數(shù)Tab.2 Vegetation index used

1.3 研究方法

1.3.1 多時相VI 曲線擬合 原始的多時相VI 曲線雖然存在明顯的變化規(guī)律，但由于無人機數(shù)據(jù)采集時的氣象條件變化、無人機及內(nèi)部傳感器擾動等內(nèi)外因素的影響，曲線會存在一定的噪音，使多時相曲線不夠平滑甚至有突變現(xiàn)象，因此需對原始曲線進行平滑處理。采用函數(shù)擬合法來重構(gòu)多時相VI曲線。常用的擬合函數(shù)包括雙logistic 函數(shù)、非對稱高斯函數(shù)和對稱高斯函數(shù)等，其中非對稱高斯函數(shù)常采用雙高斯函數(shù)形式，對稱高斯函數(shù)即單高斯函數(shù)。前2 種函數(shù)的表達式中，需要解算的未知參數(shù)個數(shù)均至少有6 個[23]，而本研究試驗區(qū)的影像期數(shù)較少，用這2 種函數(shù)易產(chǎn)生過擬合現(xiàn)象。對稱高斯函數(shù)，即單高斯函數(shù)的參數(shù)個數(shù)最少，為4個。高斯函數(shù)的基本公式為：

式中，a是高斯曲線的峰值，b是峰值對應的橫坐標，c是標準差，決定高斯鐘形曲線的寬度，d為曲線的常數(shù)項?？紤]到水稻種植前和收割后水稻田中土壤的各波段反射率無明顯差別，VI接近于0，故公式中常數(shù)項d設置為0，實際為3 個參數(shù)，不易過擬合。

1.3.2 擬合曲線特征提取與抽穗時間估算模型構(gòu)建 水稻抽穗期是水稻從營養(yǎng)生長轉(zhuǎn)為生殖生長的過渡期，植株從只有綠色的葉片發(fā)展為綠色葉片和黃綠色稻穗并存，這一變化反映在冠層光譜上。對各試驗區(qū)全部品種的多時相VI 曲線進行初步探究后發(fā)現(xiàn)，實測始穗日期都集中在擬合曲線峰值右側(cè)，實測始穗DAS與曲線峰值DAS之間存在很強的線性相關(guān)性。據(jù)此，提出本研究的抽穗時間估算模型：

其中，IHDDAS為initial heading date 的DAS，即始穗日期的DAS（簡稱為始穗DAS），PEAKDAS為峰值的DAS，a1和a2分別為線性函數(shù)關(guān)系式的斜率和截距。

1.3.3 模型應用與評價 將水稻抽穗時間估算模型首先在建模區(qū)應用，隨后將模型推廣至完全獨立的驗證區(qū)進行檢驗，并使用Pearson相關(guān)系數(shù)r、均方根誤差（Root mean square error，RMSE）和相對均方根誤差（Relative root mean square error，RRMSE）來綜合評價此模型的適用性。

2 結(jié)果與分析

2.1 VI與多時相曲線擬合方法選擇

在計算各時期的植被指數(shù)后，可獲得各品種的多時相VI序列。圖2展示了試驗區(qū)2部分水稻品種的多時相VI 散點圖。由此可見，NDVI 在觀測初期就迅速上升到全生育期內(nèi)的最大值，隨后基本保持平緩，到后期略有下降，在中后期的數(shù)值很接近，無法以此區(qū)分生育期；NDRE 則在整個生育期呈現(xiàn)先上升后下降的明顯變化，且曲線基本上左右對稱，該曲線形態(tài)對于抽穗期的識別十分有利；SR序列與NDVI 相似，且由于水稻封壟后紅波段反射率極小，SR 容易出現(xiàn)異常大的值；EVI2 表現(xiàn)優(yōu)于NDVI 和SR，但波動比NDRE 大。因此，后續(xù)研究將選用NDRE繼續(xù)進行。四者對比可發(fā)現(xiàn)由紅邊波段構(gòu)建的植被指數(shù)的優(yōu)越性，這與前人的研究結(jié)果相符[24]。

圖2 試驗區(qū)2部分水稻品種的多時相VI散點圖Fig.2 Scatter diagram of multi-time VIs of some cultivars in experiment area 2

以試驗區(qū)2 的珞優(yōu)9311 品種為例，使用3 種擬合函數(shù)對該品種的多時相NDRE 曲線進行擬合，結(jié)果如圖3 所示。原始數(shù)據(jù)在前期快速上升，在中間時期有小幅波動，但不影響總體趨勢為先上升后下降。雙高斯函數(shù)曲線雖然整體上對原始數(shù)據(jù)點擬合最貼近，但明顯出現(xiàn)雙峰現(xiàn)象，與單季水稻的NDRE 曲線只有單峰的事實不相符，這是曲線的過擬合導致的。雙logistic 函數(shù)曲線形態(tài)整體偏向于營養(yǎng)生長期，頂峰與原始數(shù)據(jù)的頂峰不重合，不能反映真實的NDRE變化情況。單高斯函數(shù)曲線雖然貼合的點最少，但與原始數(shù)據(jù)的曲線形態(tài)最為相似，頂峰位置也最為接近，是比另2種函數(shù)更合適的選擇。因此，本研究將用單高斯函數(shù)進行所有品種的多時相NDRE曲線的擬合。

圖3 3種曲線擬合方法示例Fig.3 Examples of three curve fitting methods

2.2 模型構(gòu)建與應用

為探究峰值日期與實測始穗日期之間的函數(shù)關(guān)系，由于試驗區(qū)1 的水稻品種最多且實際抽穗時間跨度最大，為保證建立的模型具有通用性，故將其作為建模區(qū)。圖4 描繪了試驗區(qū)1 部分品種的NDRE擬合曲線與實測始穗日期的相對位置。對試驗區(qū)1 的1 014 個水稻品種的多時相NDRE 曲線，用單高斯函數(shù)進行擬合，解算參數(shù)b，提取峰值DAS，隨后與這些品種的實測始穗DAS作比較（圖5）?？梢园l(fā)現(xiàn)，實測始穗DAS 與峰值DAS 高度線性相關(guān)，兩者的Pearson相關(guān)系數(shù)達到0.906。根據(jù)試驗區(qū)1的1 014 種水稻的實測始穗DAS 與峰值DAS，以最小二乘準則，解算公式（2）的參數(shù)a1和a2，得a1=1.3和a2=-18，代入式（2）中，最終可得本研究的水稻抽穗時間估算模型：

圖4 部分水稻品種的實測始穗日期與擬合曲線相對位置Fig.4 Relative position of measured initial heading date and fitting curve of some cultivars

圖5 建模區(qū)峰值DAS與實測始穗DAS散點圖Fig.5 Scatter plot of PEAKDAS and measured IHDDAS in the modeling experiment area

將此模型應用到建模區(qū)即試驗區(qū)1 的1 014 個水稻品種，可以得到這些品種的始穗時間估算值。如圖6所示，始穗DAS估算值多與實測值接近，二者構(gòu)成的二維散點都分布在1∶1 線附近。圖6 中有個別點偏離1∶1線較遠，經(jīng)統(tǒng)計，始穗DAS估算值與實測值之差在10 d 以上的只有6 個品種，占比為0.6%。絕大多數(shù)品種的差值在5 d 以內(nèi)，品種個數(shù)為910，占比達到89.7%。整體上估算值與實測值的RMSE僅為3.15 d，說明抽穗時間估算模型能夠很好地估算建模區(qū)水稻的始穗時間。

圖6 建模區(qū)模型估算結(jié)果散點圖Fig.6 Scatter plot of estimated results by the model in the modeling experiment

2.3 模型推廣與驗證

為探究水稻抽穗時間估算模型的泛化性能，將模型IHDDAS= 1.3×PEAKDAS- 18 推廣至2 個驗證區(qū)進行驗證，其中驗證區(qū)1 為試驗區(qū)2，驗證區(qū)2 為試驗區(qū)3。如圖7所示，該模型在驗證區(qū)2的應用效果非常好，40 個品種的始穗時間估算RMSE僅為2.99 d，這個誤差小于建模區(qū)；同時，經(jīng)統(tǒng)計僅有3個品種的誤差大于5 d。而在驗證區(qū)1，始穗時間估算值小于實測值的現(xiàn)象較為明顯，但誤差并不很大，48 個品種的RMSE為5.02 d，誤差超過10 d 的也僅有1個品種，整體誤差可以接受。綜合來看，由建模區(qū)建立的水稻抽穗時間估算模型在驗證區(qū)表現(xiàn)良好，能以較小的誤差估算2 個驗證區(qū)內(nèi)水稻的始穗時間。表3是模型在3個試驗區(qū)應用的統(tǒng)計結(jié)果。

圖7 驗證區(qū)模型估算結(jié)果散點圖Fig.7 Scatter plot of estimated results by the model in the verification experiment areas

表3 3個試驗區(qū)的模型應用結(jié)果Tab.3 Statistics of model application results in the three experiment areas

3 結(jié)論與討論

3.1 模型適用性討論

在本研究的3 個試驗區(qū)中，建模區(qū)和驗證區(qū)進行的試驗有很大的差別。從時間和地點上看，建模區(qū)和驗證區(qū)1 處于同一個試驗地點，但試驗時間不一致，間隔了1 a；建模區(qū)和驗證區(qū)2 的時間、地點均不一致：前者位于華中地區(qū)湖北省的亞熱帶季風氣候區(qū)，在2018 年夏季進行試驗，6—8 月平均氣溫為25～33 ℃，月均日照長度約為140 h，后者位于華南地區(qū)海南省的熱帶季風氣候區(qū)，在2018年冬春季進行試驗，1—4 月的平均氣溫為20～27 ℃，月均日照長度約為180 h。從品種上來看，建模區(qū)的試驗材料有1 014 個水稻品種，涵蓋了中國大部分地區(qū)典型的雜交水稻品系，但與驗證區(qū)1 的48 個品種和驗證區(qū)2的40個品種均不相同，僅驗證區(qū)1和2中有8個相同品種。

在生育期長短方面，水稻從抽穗到成熟經(jīng)歷的時間約為40 d[25]，因此可用抽穗期的早晚來反映全生育期的長短。建模區(qū)水稻的始穗期從播種后59 d 到103 d 不等，其中大部分在80 d 前，生育期較短；驗證區(qū)2 水稻的始穗期晚于播種后100 d，基本都在110 d 以上，生育期很長；驗證區(qū)1 水稻的始穗期位于兩者之間，生育期也比較長?？偟膩碚f，建模區(qū)的水稻屬于短生育期，驗證區(qū)的水稻屬于長生育期。值得注意的是，水稻的生育期長短與溫度、光照等外界環(huán)境條件有很大的關(guān)系[26]，陵水冬春季較低的溫度和長日照在一定程度上延長了水稻的生育期。

建模區(qū)的模型在2個與之差異大的驗證區(qū)都有良好的表現(xiàn)，這表明本研究建立的水稻抽穗時間估算模型對氣候、品種、生育期的長短均不敏感，模型具有很好的適用性。未來試驗可以在更多不同的時間、地點、品種上應用，以更廣泛地探索模型的適用性。

3.2 抽穗時間與多時相VI曲線峰值時間的關(guān)系

本研究根據(jù)始穗時間與多時相VI 曲線的峰值時間之間的線性函數(shù)關(guān)系，建立了估算效果良好的水稻抽穗時間估算模型。然而，在已有文獻中，最常用的是直接用峰值時間作為抽穗時間的估算值。為探究2 種方法的估算結(jié)果差異，用3 個試驗區(qū)的數(shù)據(jù)，對比各品種的多時相NDRE 曲線峰值時間與實際始穗時間，發(fā)現(xiàn)峰值DAS 整體上普遍小于實測始穗DAS，且實際抽穗時間越晚，這種差異越明顯。經(jīng)計算，用峰值時間直接估算實際抽穗時間的誤差遠大于本研究提出的抽穗日期估算模型的估算誤差。

這一結(jié)果與水稻生長發(fā)育期間冠層反射率的變化情況密切相關(guān)。從分蘗期到拔節(jié)孕穗期，水稻植株的葉片不斷增多增大，葉片層數(shù)的增加使得附加反射率增加，從而引起近紅外波段反射率增大，而紅邊波段反射率沒有明顯變化，NDRE 在此期間不斷增大；進入幼穗分化期，葉片不再增多并需要向穗部供給養(yǎng)分以維持其生長，此時近紅外波段反射率開始下降，紅邊波段依舊變化不大，故NDRE在此期間開始轉(zhuǎn)為減??；到抽穗期，水稻冠層結(jié)構(gòu)因為綠色稻穗從葉鞘中抽出而發(fā)生明顯變化，對冠層結(jié)構(gòu)變化敏感的紅邊波段開始上升，近紅外波段變化不大，因此抽穗期內(nèi)NDRE 減??；之后，稻穗逐漸成熟并彎曲，葉片變黃變小，使得紅邊波段反射率不斷增大，近紅外波段反射率減小，從而NDRE繼續(xù)減小[27-28]。因此，從整個生育期來看，NDRE 在幼穗分化階段已經(jīng)達到最大值，達到峰值的時間早于抽穗時間，使得峰值時間不能成為抽穗時間的理想估計值，而本研究根據(jù)峰值時間與抽穗時間的線性關(guān)系建立的模型能更準確地估算水稻的抽穗時間。

本研究提出了一種以無人機遙感為手段的多品種水稻抽穗時間估算方法。以不同氣候區(qū)域的上千品種水稻的多時相無人機多光譜遙感數(shù)據(jù)為基礎，根據(jù)多時相植被指數(shù)曲線的峰值時間與實際始穗時間的強線性關(guān)系，建立了統(tǒng)一適用的多品種水稻抽穗時間遙感估算模型。模型在建模區(qū)1 014種水稻的始穗期估算誤差RMSE達到3.15 d，在品種數(shù)為48、40 個的2 個驗證區(qū)的RMSE分別為5.02、2.99 d，整體估算精確度良好。該模型可以不區(qū)分氣候、地點、品種等差異，為多品種水稻抽穗時間監(jiān)測提供了有力途徑，是育種栽培和田間管理的理想手段。