無人機光譜分析在水稻產量預測中的應用

隋麗娜,房 建,郭立峰

(河北民族師范學院 數學與計算機科學學院,河北 承德 067000)

0 引言

水稻在我國廣泛種植,是我國主要的糧食作物,水稻的穩產高產直接關系到國家糧食安全[1]。由于種植條件的特殊性,在東北等主要水稻產地,水稻主要為大規模的集中式種植[2]。該種植模式有助于管理,但當某個生長時期出現病災和蟲災時,會對水稻產量產生嚴重影響[3]。因此,急需建立一套水稻全生長周期的產量預測模型,當一個生長時期出現問題時,通過其他生長周期進行調整,保證水稻穩產。目前,水稻產量模型主要包括單株圖像采集法[4]、飛行器圖像分析[5]和衛星光譜分析[6]。單株圖像分析對一株水稻不同生長時期進行圖像處理,分析其結節數、植株高度和植株寬度等參數,然后對其產量進行建模分析。該方法建模精度高、針對性強,但受限于樣本數量及樣品選擇,對于大范圍的水稻產量預測容易出現偏差[7]。飛行器圖像分析方法對水稻冠層圖像進行分析,背景噪聲處理難度較大,且受到稻田環境因素影響,算法精度有待提高[8]。衛星光譜分析采樣范圍過大,光譜信息龐雜,建模難度大,且獲得光譜數據成本高[9]。無人機技術的發展,使獲得低成本光譜數據成為可能。以無人機搭載光譜儀進行40m高度均速光譜采集,通過分析水稻不同生長時期的光譜,得到水稻畝產預測模型。同時,進行多組“留一法”實驗,結果表明:該方法具有良好的預測精度,性能穩定,預測精度高。該模型數據獲得方便,預測精度較高,適于大規模一線推廣。

1 系統組成

系統包括無人機光譜系統、水稻產量模型系統、模型建立依據和產量模型建立4部分,如圖1所示。

圖1 系統結構圖Fig.1 Structure of system。

無人機光譜系統保持無人40m高度勻速飛行,采用機載反射光譜傳感器檢測水稻冠層葉綠素光譜。模型建立依據,葉綠素可以表征植物生長狀態,冠層葉綠素含量可以預測水稻產量。通過無人機機載光譜儀檢測水稻冠層葉綠素光譜,由于葉綠素在紅邊區域、近紅外區域均有特征吸收波長,因此采用差值植被指數綜合考慮光譜包含信息。計算水稻全生命周期中的差值植被指數及其對應的畝產量,對其進行共線性分析,確定模型建立因子,并采用逐步線性擬合的方式建立模型。同時,采用“留一法”進行畝產預測,對模型進行線性擬合評價和統計指標評價,測試模型可靠性。

2 建模原理

太陽與地球間距離遙遠,照射到地球上的太陽光是理想的白色平行光,當白色平行光照射到水稻上時,葉綠素會吸收特定波長的光,然后將白光反射到定高飛行的無人機光譜傳觀器上。由于特定波長的光被吸收,在光譜傳感器輸出光譜曲線上形成波谷,通過分析光譜即可確定葉綠素含量。葉綠素吸收光能,將水和CO2轉化為有機物和氧氣,實現光合作用,可以作為衡量水稻長勢及產量的標志。因此,通過無人機光譜遙感計算葉綠素含量可以預測水稻最終產量。

2.1 葉綠素含量與水稻產量

水稻不同生長周期葉綠素SPAD值與畝產關系如圖2所示。

圖2 水稻不同生長周期葉綠素SPAD值與畝產關系Fig.2 Relationship between chlorophyll SPAD and rice yield in different growth period。

葉綠素吸收光能并離子化,能量被儲存在三磷酸腺苷中,最終將水和CO2轉化為有機物和氧氣,實現光合作用。由于無人機反射光譜檢測重點針對冠層葉綠素含量,現采用SPAD-502型葉綠素檢測儀對不同生長時期的水稻冠層葉片進行檢測,每葉片3點取樣,計算水稻冠層葉片SPAD值,即

(1)

其中,i為單株水稻冠層葉片數量,測試結果如圖2(a)所示。

葉綠素SPAD值從拔節期開始顯著增長,拔節期到孕穗期增長速率最大,到抽穗期達到最大,隨后含量開始下降。這是由于在抽穗期之前光合作用主要集中在植物生長上;從抽穗期開始,轉向生殖生長,葉片中的氮開始向稻穗集中轉移,葉綠素含量逐漸降低,此時水稻處于生殖生長階段。設冠層葉片葉綠素含量A為

(2)

其中,LAI為冠層葉面指數。不同生長時期冠層葉片葉綠素含量A與畝產擬合直線,相關系數如圖2(b)所示。其在孕穗期達到最高,乳熟期次之。這是由于該時期葉片發育茂盛,葉綠素含量能較好反映水稻生長狀況,抽穗期與乳熟期主要進行稻穗生殖生長,葉片變黃,葉綠素含量降低。總體上葉綠素含量與畝產相關系數均達到0.72以上,擬合精度較高,冠層葉綠素含量可以作為水稻產量的評價對象。

2.2 冠層葉綠素無人機光譜

水稻處于不同生長時期的無人機光譜如圖3所示。其中,總體上拔節期、孕穗期、抽穗期和乳熟期光譜曲線走勢一致,葉綠素對于綠光波段吸收能力很弱,故在530～600nm形成波峰,對于藍光和紅光吸收能力較強,在400～500nm形成較小波谷;在620～680nm形成較深波谷,同時在700～800nm紅光波段與近紅外線波段交接區域,反射率急劇增長,形成紅邊區域;在900～960nm近紅外波段形成波谷。由圖3可知:葉綠素具有3個特征吸收波段,其中紅邊區域與近紅外波段吸收效果顯著,因此采用采用紅邊區域與近紅外波段組合參數模式,可提高對于葉綠素含量檢測的穩定性與敏感性。

圖3 不同生長時期冠層葉綠素光譜Fig.3 The rice canopy chlorophyll spectral in different growth period。

由于水稻產量與水稻4個生長時期的冠層葉綠素含量擬合,相關系數達到0.7以上,通過冠層葉綠素含量可以有效預測水稻產量;同時,冠層葉綠素含量的無人機光譜數據具有良好的特征吸收波長,通過紅邊區域與近紅外波段組合參數模式可以有效預測其含量,因此通過分析無人機光譜可以有效的預測水稻產量。

3 建模建立

葉綠素具有3個特征吸收波段,在紅邊區域和近紅外波段吸收效果顯著,采用單一波段不能概括光譜曲線全部信息,因此采用多波段組合進行組合計算,得到水稻植被組合因子。在水稻生長全生命周期中,拔節期到孕穗期過程中,水稻冠層葉片處于高速生長狀態;當處于抽穗期和乳熟期時,冠層葉片由綠色向黃色發展,要充分考量冠層葉綠素含量的變化,因此選用對于背景變化敏感的差值植被指數[10]作為模型建立因子,計算公式如式(3)所示。其中,DNNIR為近紅外光譜反射率;DNR為紅邊區域光譜反射率。

DVI(NIR,R)=DNNIR-DNR

(3)

由于水稻生長分為拔節期、孕穗期、抽穗期和乳熟期,為了綜合不同時期對最終畝產的影響,將4個時期的差值植被指數作為模型建立的因子。

3.1 不同生長時期差值植被指數與畝產關系

水稻不同生長時期差值植被指數與畝產之間的相關系數R灰度云圖,如圖4所示。

圖4 不同時期水稻DVI(NIR,R)與產量相關系數云圖Fig. 4 Cloud chart of linear correlation coefficient between DVI(NIR,R)and rice yield in different growth period。

圖4中,紅邊波長為橫坐標,近紅外特征波長為縱坐標,每個坐標相關系數R數值用該點的灰度表示,即坐標點越亮,該點對于的相關系數R越大。具體數值參照右側灰度強度條。

圖4(a)為水稻拔節期差值植被指數與畝產相關系數云圖,該時期擬合相關系數R普遍在0.4以下,最大值為0.631 6。水稻處于拔穗期時,冠層葉片發育不成熟,葉片發育程度較低,土壤背景雜波處于主要因素。圖4(b)為孕穗期相關系數R云圖,總體上較拔節期有所提高,大部分處于0.5以上,最大值為0.771 3。這是由于冠層葉片在孕穗期顯著生長,葉綠素含量顯著增加,紅邊區域吸收能力增強。圖4(c)為抽穗期相關系數R云圖,大部分處于0.6以下,最大值為0.634 1,出現在1 080nm和1 089nm處。由于兩特征吸收波段相近,造成紅邊向量與近紅外向量線性相關較大,因此擬合關系系數R較低。圖4(d)為乳熟期相關系數R云圖,最大值出現在900nm及1 100nm處,該波段為淀粉特征吸收波段,紅光和藍光波段相關系數強度下降,葉綠素含量降低。在乳熟期水稻生長以生殖生長為主,稻穗中淀粉開始生成。不同生長時期差值植被指數與畝產相關系數最大值及出現波長如表1所示。

表1 不同生長時期DVI(NIR,R)特征波長及與畝產相關系數Table 1 Characteristic absorption wavelength and correlation coefficient between DVI(NIR,R)and rice yield in different growth period。

3.2 共線性判定

水稻在生長過程中分為4個時期,現對4個時期的差值植被指數與最終畝產量進行線性擬合,得到水稻產量模型,如式(4)所示。其中,X1～X4分別為拔節期、孕穗期、抽穗期和乳熟期差值植被指數DVI(NIR,R)向量。

Y=aX1+bX2+cX3+dX4+m

(4)

在進行多元線性擬合前必須確保X1-X4線性無關,即4向量不共線。判定依據為容忍度T和主成分分析特征值。其中,容忍度T計算公式為

(4)

其中,Rij為Xj與Xi的線性相關系數。容忍度越小,該向量與其他向量的線性先關性越高。當T<0.1時,認為線性相關性非常嚴重;主成分分析特征值E越小,線性相關性越高。現計算拔節期、孕穗期、抽穗期和乳熟期差值植被指數DVI(NIR,R)向量的容忍度和主成分特征值,如圖5所示。其中,拔節期容忍度為0.097,主成分特征值為0.14;抽穗期容忍度為0.051,主成分特征值為0.09。差值植被指數DVI(NIR,R)向量在拔節期與抽穗期高度共線,應給予舍棄。因此,以孕穗期和乳熟期差值植被指數DVI(NIR,R)向量作為模型建立因子。

圖5 共線性判定Fig.5 The collinear test。

3.3 逐步回歸分析

由共線性分析可知:拔節期與抽穗期容忍度和主成分特征值很小,線性相關性高,在模型中剔除。因此,水稻產量模型簡化為

Y=aX1+bX2+m

(5)

現采用逐步線性擬合的方法進行求解,主要步驟如下:

1)將孕穗期差值指數、乳熟期差值指數與畝產量進行線性擬合,對回歸系數βi預測值與實際值進行F檢驗。

y=Xiβi+εi(i=1,2)

(6)

孕穗期與乳熟期對應的檢驗結果為F1、F2,取F(1)=max(F1,F2),查表的置信區間為1-a情況下,臨界值為F1,當F(1)>F1時,引入第1個變量。

2)引入第2個量,建立畝產量y與X1和X2的二元線性回歸方程,對回歸系數預測值與實際值進行F檢驗,計算結果為F(1)。查表的置信區間為1-a情況下,臨界值為F2,當F(2)>F2引入第2個變量,完成產量預測模型。

現取8組樣本光譜數據,采用“留一法”進行交叉預測,以其中7組樣品光譜作為建模樣本,另一組作為預測樣本,進行8次模型預測,結果表2所示。表2中,8組預測值和實際值平均相對誤差值為7.102%,表明模型具有較高的預測精度。

表2 水稻產量模型系數及其預測值Table 2 Coefficient and predictive value of production prediction of rice。

4 模型評價

4.1 線性擬合評價

多元線性擬合的模型評價標準為線性決定系數R2和均方根差RMSE。其中,線性決定系數R2表示模型預測值與實際值的符合程度,計算公式為

(7)

均方根差RMSE為預測值與真實值的誤差平方根的均值,計算公式為

(8)

現計算采用采用“留一法”進行交叉預測,計算每一組模型R2和RMSE值,如圖6所示。

圖6 模型線性性能評價Fig. 6 Linear fitting evaluation。

其中,R2值分布在0.607～0.648之間,均值為0.630;RMSE值分布在21.612～25.145之間,均值為24.300。這說明,模型擬合線性性能良好。

4.2 統計指標評價

統計學中常將總體相對誤差Rs和預估精度P作為模型預測精度的評價標準,兩標準的計算公式為

(9)

采用式(9)計算得到8組樣品的Rs為0.0142,預估精度P為93.71。測試結果表明,該模型在統計學上具有良好的預測精度。

5 結論

為了實現大范圍的水稻產量預測,基于無人機水稻反射光譜分析建立了水稻產量模型。由于冠層葉綠素濃度可以有效表征植物生長情況,預估水稻產量。通過分析無人機測量的水稻冠層葉綠素光譜發現:在620～680nm形成較深波谷;在700～800nm紅光波段與近紅外線波段交接區域,反射率急劇增長,形成紅邊區域;在900～960nm近紅外波段形成波谷,包含豐富信息。因此,采取多波段組合進行綜合計算,以完整地提取光譜信息;選用包含紅邊區域與近紅外區域的水稻差值植被指數,作為建模因子。同時,分析差值植被指數在拔節期、孕穗期、抽穗期和乳熟期的全生長周期中的特征吸收波長,將其和畝產進行線性擬合,得到線性相關系數。對其在不同生長時期進行共線性判定,發現差值植被指數在拔節期與抽穗期的容忍度和主成分特征值均很低,存在嚴重的線性相關,因此給予舍棄,并對孕穗期與乳熟期差值植被指數與水稻畝產進行逐次線性擬合。對8組樣本進行“留一法”預測表明:模型的線性決定系數R2值分布在0.607～0.648之間,均值為0.630;RMSE值分布在21.612～25.145之間,均值為24.300。由此表明,模型精度較高。計算統計量總體相對誤差Rs=0.0142和預估精度P=93.71,表明該模型在統計學上具有較高精度。本方法數據采集方便,光譜處理簡單,適用于大規模推廣。