基于仿生優化算法的水稻葉綠素含量反演模型

李曉凱，于海業，于 躍，王洪健，張 蕾，張 昕，隋媛媛

吉林大學生物與農業工程學院，吉林 長春 130022

引 言

水稻為我國主要糧食作物，其農業生產面積在全國耕地面積中占比高，約為25%，居世界第2位。隨著人口規模的增長和生活質量的日益改善，對糧食需求量也越來越高，水稻產業得到了快速發展。據國家統計局2020年度數據顯示，水稻產量在我國三大農作物中占比約為34%, 基本能保證60%的人民對糧食的需要，所以確保水稻的質量和產量，直接關系到我國糧食安全問題。葉綠素是植物進行光合作用的主要色素，在作物體內許多重要物質和蛋白質的合成中起著重要作用，其含量與植被的光合能力、生長發育以及營養狀況密切相關。當植物受外界環境脅迫時，葉綠素含量會降低，從而間接反映植物脅迫、生長和衰老等狀況[1-2]。有研究證明其含量可用于水稻葉片氮素營養診斷及優化、稻田氮素追肥系統開發與優化、監測與評估水稻病蟲害等多方面研究[3]，具有重要的意義。葉綠素含量直接影響水稻的光譜反射率，尤其是可見光與近紅外波段[4-5]，因此基于光譜反射率反演葉綠素含量具有夯實的理論基礎。

目前基于光譜信息預測水稻葉片葉綠素含量多以單波段、多波段或全光譜信息為主。單波段光譜常結合與葉綠素含量有較好相關性的植被指數或改進后的植被指數進行建模，如比值指數、歸一化指數、三角植被指數等，但常會導致光譜有效信息無法充分利用，模型精度偏低且適用性差[6-7]。選取全波段光譜信息建模時，又由于存有復雜的冗余信息，同樣會使模型精度受到一定的影響。此外，光譜數據易受光、噪音和基線漂移等因素的干擾，影響后續分析處理和建模效果，因此需在建模進行預處理。常采用預處理手段是根據光譜信號的特征依照經驗選擇預處理方法，雖無需建立模型，但要求選擇者有敏銳的判斷力和豐富的經驗；為克服上述缺點，增強模型準確性，可采用依據預處理目的將預處理方法分類，通過全排列組合后的實際建模效果選擇最佳預處理方法，但目前預處理方法較多，耗時較長。極限學習機(extreme learning machine, ELM)與核極限學習機(kernel function extreme learning machine，KELM)在風電功率預測[8]、作物蒸散量預測[9]、利用光譜對水源類型識別[10]等領域均取得較高精度，但存有計算繁瑣且易陷入局部最優解等問題，并且在反演預測植物相關生理信息方面尚缺乏應用。

隨著人工智能技術的發展，最優算法在各個研究領域被廣泛應用。然而在求解多峰值、多變量或目標函數不可微分等復雜優化問題時，傳統優化方法常導致求解過程計算繁瑣或難尋得全局最優解等問題。仿生群智能算法憑借對目標函數與約束條件的限制較少，能高效地尋求最優解，且與其他模型結合性好，使得諸多復雜的問題得以更好地解決。為解決單純采用機器學習模型反演葉綠素含量時精確性和穩定性差且易陷入局部最小值等的問題，引入由Heidari，Mirjalili等在2019年提出的一種基于群體優化的仿生群智能算法——哈里斯鷹算法(Harris Hawks optimizer, HHO)[11]。該算法具有較強的全局收斂性及魯棒性，能夠很好地解決非線性規劃、組合優化等問題，同時可根據不同情況選擇合適的控制變量以獲得更高的尋優精度，因此備受研究人員的重視并迅速應用在各種領域。Hossein等[12]將其應用于土木工程領域，通過HHO算法優化了土質邊坡穩定條件因子計算權重，提高了預測準確度，有效地解決了土木工程中邊坡穩定性計算的難題；Du等[13]將HHO算法應用于空氣質量檢測領域，利用改進后的HHO算法對極限學習機參數進行優化，提出了一種預測精度較高的PM2.5和PM10極限學習機混合模型。HHO算法還在水火電調度優化、電力預測、圖像分割、物聯網領域、新冠病毒感染人數預測等多領域取得進展[14-16]。然而，基于HHO算法在植物生理信息預測方面的應用還有待深入研究。

為進一步提高水稻葉綠素含量預測精度與穩定性，以東北粳稻吉粳88為研究對象，采用光譜檢測技術，首先依據多波段建模對水稻葉片光譜數據的預處理方法展開研究?；?類預處理方法全排列組合的優勢，在HHO模型基礎上，結合連續投影法將各預處理組合所篩選的特征波段輸入模型中，依據實際建模效果選擇最佳預處理方法。其次，利用HHO仿生優化算法對KELM訓練過程進行優化，提出一種HHO-KELM融合模型，進一步提高水稻葉片葉綠素含量反演模型的準確性與泛化性。旨在探究光譜預處理的必要性及尋求最優的預處理方法，拓寬群智能算法在植物生理信息預測方向的研究，建立了更精準穩定的水稻葉片葉綠素含量診斷模型，為作物葉片葉綠素含量診斷提供新的思路與理論依據。

1 實驗部分

1.1 試驗設計

實驗于吉林大學農業實驗基地內進行，水稻品種選用吉粳88。于2021年5月20日移栽，實驗共設5個氮肥梯度，分別為N0(0)，N1(120)，N2(160)，N3(200)，N4(240)，單位為kg·hm-2，氮肥梯度采用網格種植，各網格面積為10 m×10 m，各水平設置3次重復。每小區選取20株長勢旺盛、健康無病蟲害的水稻葉片作為實驗待測樣本。

1.2 數據采集與處理

實驗選定在陽光明媚、無風的情況下進行，實驗時間為2021年7月10日(分蘗盛期)10：00—14：00，對300株樣本進行數據采集，每株測量倒二、三完全展開葉，采集葉片的高光譜和葉綠素含量數據，樣本量共為600個。高光譜數據采用美國Analytical Spectral Devices公司生產的手持式地物光譜儀FieldSpec HandHeld 2測定，測量范圍325～1 075 nm，采樣間隔1.4 nm，分辨率3 nm，采用葉夾式對每片葉取避開葉脈的1/3，1/2和2/3處進行測量。葉綠素含量采用SPAD-502葉綠素儀測定(日本柯尼卡美能達)，同上述位置3次重復測量。SPAD-502讀數與葉綠素含量密切相關[17]，因此以其值表示葉綠素含量。以上2種數據均取平均值作為該葉片高光譜與葉綠素含量數據。數據處理與分析軟件為ViewSpecPro，Matlab R2018b和 Origin 2018。

在光譜預數據處理前，首先去除噪聲，保留401～1 000 nm波段，采用蒙特卡洛交叉驗證法(MCCV)在篩選光譜值與化學值的同時進行異常樣本剔除，保留樣本561個。然后應用不同預處理方法結合連續投影算法(successive projections algorithm, SPA)進行特征波段選擇，降低光譜數據維度。樣本采用3∶1隨機劃分訓練集(420個)和測試集(141個)，分別建立KELM模型以及優化后的HHO-KELM模型。

1.3 光譜預處理方法

表1 預處理方法種類

1.4 模型的建立

1.4.1 核極限學習機

極限學習機(ELM)是由南洋理工大學的Huang等[18]于2004年提出的基于單隱藏層前饋神經網絡的算法，最初是為改進反向傳播算法(BP)以提升學習效率和簡化學習參數的設定。相較于傳統前饋神經網絡訓練速度慢、易陷入局部極小值點、訓練過程不穩定等缺點，ELM算法可以在輸入層與隱含層之間隨機產生連接權值(w)和隱藏層神經元的閾值(b)，且后續訓練過程中無需調整，只需要設置隱含層神經元的個數即可獲得最優解。

為進一步增強ELM的泛化能力以及模型穩定性，Huang等[19]通過深度分析ELM與支持向量機(SVM)的原理，提出將ELM與SVM中的核函數結合的改進算法KELM，能夠在保留ELM優點的基礎上提高模型的適應性與穩定性，使模型的預測性能得到有效提升。其最終表達式為

(1)

式(1)中：x為輸入向量，(x1,x2, …,xn)為訓練樣本輸入向量，K()為核函數，n為樣本數量，I為單位矩陣，C為正則化系數，ΩELM為核矩陣，L為期望輸出。

1.5 哈里斯鷹優化算法

1.5.1 算法仿生原理

HHO算法是依據美國亞利桑那州南部猛禽哈里斯鷹在捕食時的合作行為和追逐策略。根據場景動態性和獵物逃跑模式，對哈里斯鷹群進行了數學模擬，并結合萊維飛行(Levy-Flight)實現對復雜多維問題優化求解。在合作行為中，群鷹各自分工不同，尋找目標、驅趕、狩獵各有所司。在追逐策略中，主要為“突襲圍捕”，突襲指群鷹從不同方向合作突襲獵物，同時向獵物周圍匯聚，一次性捕獲；圍捕指初次失敗致使獵物外逃時，在靠近獵物的地方進行多次快速短距離的突襲，并根據獵物的反應和躲避方向做出不同的追逐策略。這兩種追逐方式可看作一種動態博弈過程。由此建立了新的數學模型來處理各種優化問題[8]。在HHO中，哈里斯鷹是候選解，獵物是隨迭代逼近的最優解，分為全局探索階段、探索與開發轉換階段和局部開發階段，局部開發階段又分為軟包圍、硬包圍、漸進式快速俯沖的軟包圍以及漸進式快速俯沖的硬包圍四種策略。

1.5.2 具體算法與流程

(1)探索階段

哈里斯鷹群體分散度很高，其個體隨機潛伏在某些位置，利用敏銳的雙眼偵查環境、追蹤獵物，并根據兩種策略等待并發現獵物：當q<0.5時，哈里斯鷹會根據其他成員和獵物的位置進行棲息；當q≥0.5時，哈里斯鷹會在種群活動范圍內隨機棲息在某棵樹上。

X(t+1)=

(2)

式(2)中，X(t)和X(t+1)分別為當前和下一次迭代時個體的位置，t為迭代次數，Xrand(t)和Xrabbit(t)分別為隨機選出和具有最優適應度的個體位置，r1，r2，r3，r4和q都是[0, 1]區間的隨機數，ub和lb為搜索空間變量的上下界。根據q來隨機選擇所采用的策略，Xm(t)為個體平均位置，其表達式為

(3)

式(3)中，N為種群規模,Xi(t)為種群中第i個個體的位置。

(2)探索與開發階段的轉換

根據獵物的逃逸跑量在探索和開發行為之間轉換，獵物逃逸能量(E)定義為

E=2E0(1-t/T)

(4)

式(4)中，E0是獵物的初始能量，為[-1, 1]之間的隨機數，隨迭代自動更新；t為迭代次數；T為最大迭代次數。當|E|≥1時進入全局探索階段，|E|<1當時進入局部開發階段。

(3)局部開發階段

在此階段中，哈里斯鷹開始對獵物進行突襲追捕。根據逃逸能量E及區間(0, 1)的隨機數r的取值大小來選擇四種攻擊策略捕獲獵物，具體體現為：

①軟包圍

當0.5≤|E|<1且r≥0.5時，獵物有足夠的體力去逃跑，而鷹群通過環繞追逐獵物，使其疲勞。此階段數學模型為式(6)

Xi(t+1)=Xrabbit(t)-Xi(t)-E|JXrabbit(t)-Xi(t)|

(5)

式(5)中，J=2(1-r5)為獵物的隨機跳躍強度，r5為區間(0, 1)上的隨機數。

②硬包圍

當|E|<0.5且r≥0.5時，獵物筋疲力盡已不足以逃脫追捕，此階段的模型為

Xi(t+1)=Xrabbit(t)+E|Xrabbit(t)-Xi(t)|

(6)

③漸進式快速俯沖式軟包圍

當r<0.5，|E|≥0.5時，獵物仍有能量去逃跑，而鷹群將會產生一個環繞軟包圍來追捕獵物。為了模擬獵物的逃跑模式和跳躍動作，將Lévy函數LF集成在HHO算法中，此階段模型為

(7)

式(7)中，F(·)表示最小化問題的適應度值，Y和Z分別為

Y=Xrabbit(t)-E|JXrabbit(t)-Xi(t)|

(8)

Z=Y+S×LF(D)

(9)

式(8)和式(9)中，S為大小為D×1的(0, 1)上的隨機向量，D為目標函數涉及的維度，LF(·)為服從Levy-Flight 分布的向量，其一維計算方法為式(10)

(10)

(11)

式(11)中，μ和ν為(0, 1)上的隨機數，β設置為1.5。

④漸進式快速俯沖式硬包圍

當r<0.5，|E|<0.5時，獵物沒有足夠的能量去逃跑，鷹群構建了一個環繞硬包圍來捕獲獵物，此階段的數學模型為

(12)

式(12)中，Y和Z分別為式(13)和式(14)

Y=Xrabbit(t)-E|JXrabbit(t)-Xm(t)|

(13)

Z=Y+S×LF(D)

(14)

算法流程如圖1所示。

圖1 HHO算法流程示意圖

2 結果與討論

2.1 優化預處理參數

在建模之前需對KELM的正則化系數C和核函數參數S進行預處理，預處理方法中SG平滑的窗口參數、CWT的小波函數和分解尺度參數，均需結合SPA進行優化。SG平滑優化選取窗口數從3～15進行篩選，間隔為2，以模型的RMSE值為評價標準，RMSE最小值對應的窗口即為最佳窗口數。CWT優化方法為選取Daubechies(db2, db3, …, db15)，Coiflets(coif1, coif2, …, coif5)，Symmlets(sym2, sym3, …, sym7)共25個小波函數，分解尺度從1～10，同樣以模型的RMSE值為評價標準，最小值為最佳小波函數和分解尺度。最終尋優結果如表2所示。

表2 預處理方法及KELM參數優化結果

2.2 KELM的建模效果

通過對預處理方法優化后，對樣本集進行KELM建模，其結果如圖2所示。由圖可知，采用基線校正與尺度縮放以及兩類方法結合可顯著提高建模效果，當引入CWT時，發現在大部分情況下可顯著提高預測精度，但當加入散射校正后，除結合CWT外的其余任意種類預處理效果反而下降，尤其是引入MSC后顯著下降。最終選取較優預處理組合為CWT+MMS，CWT+MSC+SG+SS和CWT+SS。

圖2 預處理方法的KELM模型測試集決定系數

2.3 HHO優化KELM模型

基礎KELM中C和S取尋優后的固定值，選取高斯核函數(RBF)為核函數。利用HHO算法優化時，依據其尋優特性對KELM的C和S選取進行優化，同樣選取核函數為RBF，進而形成對比以驗證模型預測精度的變化情況。HHO中的適應度函數選擇訓練集誤差的均方誤差(MSE)，其值越小表明預測的數據與實際葉綠素數據擬合性越好。適應度值定義如式(15)

fitness=argmin(MSEpridect)

(15)

HHO-KELM算法流程如圖3所示。

注：點畫線框內為HHO優化KELM

2.4 KELM與HHO-KELM建模

表3 三種預處理組合所對應的KELM和HHO-KELM模型結果

兩模型所對應的測試集訓練結果和預測誤差如圖4所示。由圖4(a)可看出HHO-KELM比KELM模型預測準確度更高，觀察圖4(b)和(c)得KELM在訓練集和測試集中預測存在較大偏差且穩定性差，說明基礎KELM有時無法充分利用數據信息，所建模型預測性能較差。而利用HHO算法優化后可以較好地解決，表明利用HHO算法可以對KELM訓練過程進行優化，自動最優調節KELM模型參數，使統計分析結果更準確且更有說服力。圖4(a)中測試集第59，62和76個樣本點的模型偏差顯著，經反復實驗發現此3點為MCCV未篩選出的異常樣本。由此證明了利用HHO算法優化數學模型用于反演植物生理信息的可行性，為今后尋求更優群智能算法用于研究精準、無損和穩定的植物生理信息檢測技術提供了理論依據。

圖4 CWT+MMS所對應的測試集預測結果對比圖與散點圖

基于HHO-KELM模型建立植被生理參數反演模型的研究中仍有一些問題有待解決，如本研究中以水稻的某一生育期為研究對象，將其研究結果應用于其他植物或不同生育期的可靠性有待進一步考察。是否可嘗試結合多波段植被指數以減弱誤差影響，再提升模型的準確性與適用性。其次是否所有的群智能算法都可應用于機器學習模型，從而提升植物生理參數的預測性能，還需更多驗證。最后在反復實驗中發現，因樣本數量的關系，HHO算法易出現過擬合現象。對于植物生理參數反演模型的可靠性和普適性仍需更深入的研究與實踐。

3 結 論