高光譜和NSGA2-ELM算法的粳稻葉片氮素含量反演

馮 帥，曹英麗,2*，許童羽,2，于豐華,2，陳春玲,2，趙冬雪，金 彥

1. 沈陽農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院，遼寧 沈陽 110161 2. 沈陽農(nóng)業(yè)大學(xué)遼寧省農(nóng)業(yè)信息化工程技術(shù)中心，遼寧 沈陽 110161

引 言

氮素是農(nóng)作物生長發(fā)育過程中重要的營養(yǎng)成分，實時監(jiān)測和評估農(nóng)作物的氮素含量對于農(nóng)作物田間精準管理和長勢預(yù)測等均具有十分重要的意義[1]。

目前，采用高光譜檢測法對果蔬[2-3]和糧食作物[4-6]的氮素營養(yǎng)診斷已成為國內(nèi)外學(xué)者研究的主要內(nèi)容。王樹文[7]等研究表明基于主成分分析和相關(guān)分析結(jié)合多元回歸分析模型的差值指數(shù)、 多變量單波段指數(shù)等模型反演效果較好，預(yù)測集R2為0.869，RMSE為0.085。劉明博[8]等采用連續(xù)投影法(SPA)篩選的有效波段、 光譜指數(shù)RVI、 NDVI以及全光譜波段構(gòu)建多種水稻葉片氮素含量反演模型。對比發(fā)現(xiàn)，基于SPA有效波段構(gòu)建的模型的估測效果明顯優(yōu)于光譜指數(shù)所建，但略差于全光譜波段所建模型。Tian[9]等通過分析多種高光譜植被指數(shù)與水稻葉片氮素含量的定量關(guān)系，得出采用綠色比率指數(shù)SR(R553,R537)反演葉片氮素含量具有最佳估測精度。Du[10]等采用高光譜激光雷達(HSL)技術(shù)構(gòu)建兩種積分指數(shù)NOAC和RII反演水稻葉片全氮含量(LNC)。方美紅[11]等采用小波系數(shù)構(gòu)建水稻葉片氮含量反演模型，研究表明該模型有較高估測精度，預(yù)測值與估測值的復(fù)相關(guān)系數(shù)高達0.99，顯著優(yōu)于傳統(tǒng)光譜指數(shù)反演模型。

粳稻生長過程中受氮肥影響較為明顯，不同施氮水平下的粳稻生長趨勢有一定的差異，葉片對自然光的吸收和反射也表現(xiàn)出不同的變化，進而對粳稻產(chǎn)量產(chǎn)生較大影響。鑒于此，本研究為實現(xiàn)粳稻葉片氮素含量的高效、 快速和精準反演，分別采用特征波段組合和光譜植被指數(shù)組合作為輸入，構(gòu)建3種氮素含量反演模型并對比分析，進而提出基于高光譜的粳稻葉片氮素含量反演方法，以期為粳稻氮素營養(yǎng)診斷和田間精準管理提供科學(xué)依據(jù)和理論支持。

1 實驗部分

1.1 試驗設(shè)計

粳稻小區(qū)試驗于2018年6月至8月在遼寧省沈陽市沈河區(qū)沈陽農(nóng)業(yè)大學(xué)南區(qū)水稻試驗基地(118°53′E，38°43′N，平均海拔40 m)開展，供試品種為沈稻9816。共劃分12個小區(qū), 每個試驗小區(qū)面積為30 m2(7.61 m×4.20 m), 如圖1所示。人為設(shè)置小區(qū)試驗無肥、 低氮、 正常和高氮4種情況，共設(shè)4個施氮梯度，分別為N0(不含氮)、 N1(150 kg·hm-2)、 N2(240 kg·hm-2)、 N3(330 kg·hm-2)，每個水平3次重復(fù)。同時各試驗小區(qū)之間采取隔離措施，保證小區(qū)之間水肥不互相滲透，其他田間管理水平均按當(dāng)?shù)卣Ｋ竭M行。

圖1 12個粳稻小區(qū)分布圖Fig.1 Distribution of 12 rice plots

1.2 葉片高光譜數(shù)據(jù)采集

采用美國海洋光學(xué)公司生產(chǎn)的光纖光譜儀HR2000+，其光譜范圍為400～1 000 nm，光譜分辨率為0.45 nm。按標號將磨碎的葉片粉末放置于操作平臺，光譜探頭緊壓在葉片粉末上；測量前通過黑白板校正。通過自帶OceanView軟件完成粳稻葉片高光譜數(shù)據(jù)的采集。

1.3 葉片氮素含量測定

從每個試驗小區(qū)采樣點獲取粳稻不同部位葉片20片左右，分別裝入自封袋中并標注小區(qū)名稱和編號，立即帶回實驗室。在室內(nèi)，首先對葉片進行洗滌，去除葉片表面灰塵等無用物質(zhì)，然后在105 ℃條件下殺青30 min，并在70 ℃的烘箱中烘干至恒量，稱量粉碎。最后采用凱氏定氮法測定粳稻葉片氮素含量。所采集的數(shù)據(jù)樣本中去除粗大誤差，最終共得到粳稻葉片氮素含量有效數(shù)據(jù)280組，其概率密度函數(shù)如圖2所示。

圖2 粳稻280組葉片氮素含量的概率密度函數(shù)Fig.2 Probabilistic density function of nitrogen content in 280 groups of japonica rice leaves

由圖2可知，280組粳稻葉片氮素含量數(shù)據(jù)呈正態(tài)分布，均值為2.860 mg·g-1，最大值為4.530 mg·g-1，最小值為1.060 mg·g-1，標準差為0.825 mg·g-1，變異系數(shù)為28.846%，滿足氮素含量反演要求。同時采用Kennard-Stone算法(KS)將樣本按照訓(xùn)練集與驗證集3∶1的比例進行劃分，其氮素含量統(tǒng)計表如表1所示。

表1 訓(xùn)練集與驗證集粳稻葉片氮素含量數(shù)據(jù)統(tǒng)計表Table 1 Statistical table of nitrogen content in leaves of japonica rice in training set and validation set

1.4 光譜特征參數(shù)的篩選

1.4.1 光譜特征波段篩選方法

高光譜數(shù)據(jù)具有較高的分辨率，因此在光譜數(shù)據(jù)中包含了大量的數(shù)據(jù)信息，但在龐大的數(shù)據(jù)信息中仍存在一定的冗余信息，給后續(xù)的研究帶來了一定困擾。篩選特征波段進可減少數(shù)據(jù)冗余，獲得對模型起關(guān)鍵性作用的特征波段是不可或缺的。分別采用隨機青蛙法(random frog)與迭代和保留信息變量算法(iteratively retaining informative variables，IRIV)兩種方法篩選特征波段。隨機青蛙法(random frog)采用蒙特卡洛法通過維度轉(zhuǎn)換進行采樣，以每個波段變量的選擇頻率作為剔除冗余變量的依據(jù)，最終得到最佳特征波段。迭代和保留信息變量算法是基于二進制矩陣重排過濾器提出的特征變量選擇算法，通過多次迭代，計算包含某一波段和不包含該波段變量時的RMSECV平均值，得到兩個均值之差DMEAN和曼-惠特尼檢驗的P值，從而判斷波段的重要性，其判斷規(guī)則如表2所示[12]，并依據(jù)該規(guī)則消除無信息和干擾波段變量，保留強和弱波段變量，最后通過反向消除獲得最佳特征波段。

表2 迭代和保留信息變量算法判斷規(guī)則Table 2 The judgment rule of iteratively retaining informative variables

1.4.2 光譜植被指數(shù)的構(gòu)建

植被指數(shù)常常被用來精準監(jiān)測植被生長狀況和生物量信息等植被信息，同時具有增強植被信息，減弱非植被信息的作用，因此選取3種對粳稻葉片氮素含量具有較好反演效果的高光譜植被指數(shù)，即比值植被指數(shù)(ratio spectral index，RSI)、 差值植被指數(shù)(difference spectral index, DSI)以及歸一化植被指數(shù)(normalized difference spectral index, NDSI)。其計算公式如式(1)—式(3)

RSI=Ri/Rj

(1)

DSI=Ri-Rj

(2)

NDSI=Ri-Rj/Ri+Rj

(3)

式中，Ri與Rj分別表示在400～1 000 nm范圍內(nèi)的光譜反射率值。

1.5 反演模型的構(gòu)建與評價標準

分別將特征波段組合和植被指數(shù)作為模型輸入，實測粳稻葉片氮素含量作為輸出，構(gòu)建BP神經(jīng)網(wǎng)絡(luò)、 支持向量機和NSGA2-ELM三種反演模型。由于ELM模型采用隨機生成輸入權(quán)重和隱含層偏向值，這將導(dǎo)致ELM輸出波動大，模型不穩(wěn)定。因此提出一種將NSGA2算法與ELM模型相結(jié)合的混合模型算法。具體NSGA2優(yōu)化過程如下：

(1)根據(jù)樣本數(shù)據(jù)集確定ELM的網(wǎng)絡(luò)拓撲結(jié)構(gòu)，將神經(jīng)元之間的權(quán)重和偏向值構(gòu)成實數(shù)向量，用以表示種群M中的個體。同時隨機生成實數(shù)向量的初始值構(gòu)成大小為N的第一代父代種群P。

(2)對父代種群進行非支配排序，并采用傳統(tǒng)的遺傳算法對父代種群進行選擇、 交叉和變異操作產(chǎn)生大小為N的子代種群P1。將種群P和P1合并為大小為2N的種群B。

(3)對種群進行非支配排序，獲得非支配解的前端Ft，即為非支配面的F1，F(xiàn)2和F3，并計算擁擠度。之后采用精英保留策略篩選最優(yōu)個體，即由于子代和父代種群個體均包含在種群B中，則采用非支配排序后的F1中的個體為種群B中最佳的，因此現(xiàn)將F1全部個體放入新父代種群P3中。若P3小于N，則繼續(xù)將F2中的個體加入種群P3中。若P3仍小于N，則對F3進行擁擠度排序，取N-|P3|個種群個體添加至P3，直至P3種群大小為N。然后采用遺傳算法對種群P3進行選擇、 交叉和變異產(chǎn)生新的種群P4。

(4)重復(fù)n次上述計算過程，達到設(shè)定的最大迭代次數(shù)則停止迭代，得到最佳ELM的最佳權(quán)重和偏向值，完成優(yōu)化。

同時為了檢測粳稻葉片氮素含量反演模型的準確性和可靠性，采用決定系數(shù)R2和均方根誤差RMSE兩個評價指標檢測模型的擬合效果和估測能力。其計算公式如式(4)和式(5)

(4)

(5)

整體技術(shù)路線如圖3所示。

圖3 研究技術(shù)路線圖Fig.3 Research technology roadmap

2 結(jié)果與討論

2.1 粳稻葉片氮素含量的光譜分析

由圖4可知，在不同施氮水平下的粳稻葉片光譜反射率曲線變化趨勢大致相同。在可見光波段范圍(400～770 nm)光譜反射率均呈現(xiàn)出“先升后降”的變化規(guī)律。同時，在綠波段(550 nm附近)出現(xiàn)明顯的反射峰，在紅波段(680 nm附近)出現(xiàn)明顯的吸收谷。在680～770 nm范圍內(nèi)不同施氮水平的光譜反射率曲線變化基本一致，反射率急劇升高。在近紅外波段范圍(>770 nm)形成較高的反射平臺。然而，在可見光和近紅外波段范圍，各個施氮水平的葉片光譜反射率表現(xiàn)出相反的變化規(guī)律。在可見光波段范圍，光譜反射率隨著施氮水平的升高而減小，尤其在“紅峰”處極其明顯，在近紅外波段范圍，光譜反射率隨著施氮水平的升高而增大，尤其在“近紅外反射平臺”處較為明顯。

圖4 不同施氮水平下的粳稻葉片平均光譜反射率N0：不含氮；N1：150 kg·hm-2；N2：240 kg·hm-2；N3：330 kg·hm-2Fig.4 Average spectral reflectance of japonica rice leaves under different nitrogen levelsN0: Nitrogen-free; N1: 150 kg·hm-2;N2: 240 kg·hm-2; N3: 330 kg·hm-2

2.2 光譜特征波段的篩選

鑒于高光譜的波段數(shù)較多，且存在無用和干擾波段，對建模的效率和精度產(chǎn)生較大影響。因此首先采用隨機青蛙算法(random-frog)初選特征波段，去除冗余信息。確定random-frog算法的潛在變量的最大數(shù)為6，初始采樣變量數(shù)為1 000，但由于random-frog算法以蒙特卡洛法為篩選原理，每次篩選出的特征波段略有差異。因此共運行random-frog算法50次，最后取結(jié)果的平均值作為特征波段的判斷依據(jù)，其結(jié)果如圖5所示。由圖5可知，random-frog算法篩選的結(jié)果較好，大部分光譜波段的選樣概率較小，小部分光譜波段的選樣概率較大。同時，設(shè)定篩選閾值為0.1，選樣概率大于0.1的波段被認定為特征波段，最終篩選出23個波段作為特征波段。

圖5 隨機青蛙算法篩選結(jié)果圖Fig.5 The screening result graph of random frog

random-frog方法篩選得到23個特征波段，波段數(shù)量較大，且random-frog方法在篩選過程具有隨機性，每次的篩選結(jié)果略有差異，可能仍存在無信息波段。因此，為提高模型的穩(wěn)健性和準確性，采用迭代和保留信息變量算法(IRIV)對波段進行進一步篩選，消除冗余信息和無用信息。

通過多次測試確定IRIV算法的最大主因子數(shù)為10，交叉驗證次數(shù)為5。IRIV算法篩選特征波段共進行3輪，如圖6(a)所示。前2輪波段變量個數(shù)逐漸減少，從23個波段減少到10個，其中在第2輪迭代過程中基本剔除了無用信息波段和干擾波段，余下波段的DMEAN和P值如圖6(b)所示，在第3輪采用反向消除，最終得到8個特征波段，其中弱波段變量僅有1個為439.6 nm，強波段變量僅有7個，分別為414.2，430.9，447.9，682.7，685.4，686.3和999.1 nm。此時，大部分波段集中于400～700 nm之間，即在可見光波段區(qū)域之間有7個特征波段，僅有1個特征波段位于近紅外波段區(qū)域。從篩選結(jié)果來看，雖然各個波段間仍存在少量的共線性，但相比于初次篩選，已較好地降低了波段間的共線性和冗余信息。

圖6 IRIV算法篩選過程圖(a)：IRIV迭代保留波段個數(shù)；(b)：第2輪波段的DMEAN和P值Fig.6 IRIV algorithm screeniog process diagram(a)：Number of retained variables in iterative round of IRIV；(b)：Dmean and P values in second iterat

2.3 植被指數(shù)的篩選

圖7為400～1 000 nm中任意兩個光譜波段隨機組合構(gòu)成的DSI(Ri,Rj)，RSI(Ri,Rj)和NDSI(Ri,Rj)與粳稻葉片氮素含量的決定系數(shù)等勢圖。從決定系數(shù)等勢圖中可知與葉片氮素含量相關(guān)性最優(yōu)的光譜植被指數(shù)和相關(guān)性最好的波段范圍。在圖7(a)中，對于DSI(Ri,Rj)而言，在707.3～852.5與625.5～692 nm的光譜波段組合與葉片氮素含量較好，決定系數(shù)R2均大于0.7，其中DSI(R648.1,R738.1)與葉片氮素含量的相關(guān)性最好，R2為0.811 4。相比DSI(Ri,Rj)而言，RSI(Ri,Rj)與葉片氮素含量相關(guān)性較好的波段范圍較窄，多集中于可見光波段范圍[圖7(b)]。其中在596.1～686.2和507.2～607.4 nm的波段組合相關(guān)性最佳，R2均大于0.7，最好的波段組合構(gòu)建的RSI植被指數(shù)為RSI(R532.8,R677.3)，R2為0.829 7。

圖7 任意兩波段組合構(gòu)成的DSI (a), RSI (b)和NDSI (c)與葉片氮素含量的決定系數(shù)等勢圖>Fig.7 Contour of R2 between DSI (a), RSI (b), NDSI (c) and leaf nitrogen content in any two bands

圖7(c)為任意兩個波段組合構(gòu)建的NDSI(Ri,Rj)與葉片氮素含量的決定系數(shù)等勢圖。在601.5～682.7與514.1～609.2 nm的光譜波段組合與氮素含量的相關(guān)性較好，R2大于0.7，其中NDSI(R654.8,R532.9)與葉片氮素含量相關(guān)性最好，R2為0.816 9。

2.4 粳稻葉片氮素含量反演建模與分析

2.4.1 BP神經(jīng)網(wǎng)絡(luò)反演建模

采用random-frog和IRIV算法篩選的8個特征波段和全生長期的3種較好的植被指數(shù)(DSI(R648.1,R738.1)，RSI(R532.8,R677.3)和NDSI(R654.8,R532.9))分別作為模型輸入，粳稻葉片氮素含量作為模型輸出，構(gòu)建BP神經(jīng)網(wǎng)絡(luò)反演模型。將Tansig和Purelin法分別設(shè)置為隱含層和輸出層的傳遞函數(shù)，將Trainlm法設(shè)置為訓(xùn)練函數(shù)，訓(xùn)練最大迭代次數(shù)為1 000次，學(xué)習(xí)速率lr和訓(xùn)練精度goal分別均為0.1和0.01。同時，通過逐步試驗得出隱含層神經(jīng)節(jié)點個數(shù)為6時，模型估測效果最佳。建模結(jié)果如圖8所示。

圖8 BP神經(jīng)網(wǎng)絡(luò)建模結(jié)果(a)：特征波段組合；(b)：植被指數(shù)組合Fig.8 Modeling results of BP neural network(a)：Characteristic band combination；(b)：Vegetation index combination

由圖8可知，以8個特征波段組合作為輸入構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)反演模型，其建模效果要優(yōu)于植被指數(shù)組合。模型訓(xùn)練集的決定系數(shù)R2和均方根誤差RMSE分別為0.741 5和0.430 1，驗證集的R2和RMSE分別為0.782 2和0.365 6。由此可見，采用8個特征波段組合作為BP神經(jīng)網(wǎng)絡(luò)的輸入進行建模，提高了反演粳稻葉片氮素含量的預(yù)測能力。

2.4.2 支持向量機反演建模

分別采用特征波段組合和植被指數(shù)組合作為支持向量機(SVR)模型的輸入，構(gòu)建粳稻葉片氮素含量反演模型。將SVR的核函數(shù)設(shè)置為Sigmoid，特征波段組合和植被指數(shù)組合的懲罰因子c分別為1.3和2，Sigmoid核函數(shù)的系數(shù)g分別為10.1和9，損失函數(shù)值分別為0.17和0.1，SVR建模結(jié)果如圖9所示。

圖9 支持向量機建模結(jié)果(a)：特征波段組合；(b)：植被指數(shù)組合Fig.9 Modeling results of SVR(a)：Characteristic band combination；(b)：Vegetation index combination

由圖9比較可知，無論是以特征波段組合還是以植被指數(shù)組合作為輸入所構(gòu)建的SVR模型，可靠性和估測能力具有一致性。從驗證集的建模結(jié)果來看，特征波段組合的SVR模型的R2略高于植被指數(shù)組合，但其RMSE略低于植被指數(shù)組合。特征波段組合和植被指數(shù)組合的訓(xùn)練集R2和RMSE分別為：0.720 4，0.445 9和0.747 6和0.420 6，驗證集R2和RMSE分別為：0.741 2，0.569 6和0.739 7和0.416 7。

2.4.3 非支配的精英策略遺傳算法優(yōu)化極限學(xué)習(xí)機反演建模

采用random-frog和IRIV算法篩選的8個特征波段和全生長期的3種較好的植被指數(shù)(DSI(R648.1,R738.1)，RSI(R532.8,R677.3)和NDSI(R654.8,R532.9))分別作為非支配的精英策略遺傳算法優(yōu)化極限學(xué)習(xí)機(NSGA2-ELM)模型的輸入，構(gòu)建粳稻葉片氮素含量反演模型。將特征波段組合和植被指數(shù)組合的NSGA2-ELM模型的種群個數(shù)均設(shè)置為30，交叉概率和變異概率分別設(shè)置為0.96，0.97和0.001，0.001，最大擁擠距離均為10 000，保留最優(yōu)個體比例均為0.1，迭代次數(shù)均為50次。建模結(jié)果如圖10所示。

圖10 NSGA2-ELM建模結(jié)果(a)：特征波段組合；(b)：植被指數(shù)組合Fig.10 Modeling results of NSGA2-ELM(a)：Characteristic band combination；(b)：Vegetation index combination

由圖10可知，以特征波段組合作為NSGA2-ELM模型的輸入構(gòu)建模型的估測效果明顯優(yōu)于植被指數(shù)組合，模型訓(xùn)練集的決定系數(shù)R2和均方根誤差RMSE分別為0.817 2和0.355 5，驗證集的R2和RMSE分別為0.849 7和0.301 1。相比BP神經(jīng)網(wǎng)絡(luò)模型和支持向量機模型，NSGA2-ELM模型無論是以特征波段組合還是以植被指數(shù)組合作為輸入，在模型精度和估測能力上都有顯著提高，說明采用NSGA2算法優(yōu)化ELM模型對粳稻葉片氮素含量的預(yù)測具有較大的優(yōu)勢。分析原因在于NSGA2算法能夠保持種群的多樣性，并引入精英策略，不僅能增大樣本采集空間，而且能較好地防止最優(yōu)個體丟失，也避免了類似BP神經(jīng)網(wǎng)絡(luò)等模型陷入過擬合和局部最優(yōu)現(xiàn)象，從而使模型具有更好的非線性映射能力和魯棒性。同時從不同輸入類型來看，3種反演模型多以特征波段組合作為輸入，模型估測效果較好。說明使用random-frog和IRIV篩選出的特征波段具有較好的代表性，包含了粳稻葉片高光譜的大部分特征。因此采用random-frog和IRIV方法獲取的特征波段所建立的反演模型，具有更好的穩(wěn)定性和估測能力。

3 結(jié) 論

以粳稻葉片高光譜數(shù)據(jù)為數(shù)據(jù)源，先后采用隨機青蛙算法與迭代和保留信息變量算法篩選特征波段，并構(gòu)建3種較優(yōu)植被指數(shù)，分別以特征波段組合和植被指數(shù)組合作為模型輸入，構(gòu)建BP神經(jīng)網(wǎng)絡(luò)、 支持向量機和非支配的精英策略遺傳算法優(yōu)化極限學(xué)習(xí)機粳稻葉片氮素含量反演模型并比較其優(yōu)劣。主要結(jié)論如下：

(1)針對高光譜較多的光譜波段，通過采用random-frog和IRIV算法相結(jié)合的方式篩選光譜特征波段，共得到8個特征波段，其中可見光波段有7個，分別為414.2, 430.9, 439.6, 447.9, 682.7, 685.4和686.3 nm，近紅外波段僅有1個為999.1 nm。結(jié)果表明能夠較好地剔除無用和干擾光譜波段，大大降低了波段間的共線性和冗余信息，提高了建模效率和精度。

(2)分別采用特征波段組合和植被指數(shù)組合作為模型輸入，構(gòu)建3種粳稻葉片氮素含量反演模型。對比分析得出，NSGA2-ELM模型的估測效果要優(yōu)于BP神經(jīng)網(wǎng)絡(luò)模型和SVR模型，這是由于NSGA2算法能保持種群的多樣性，增大了樣本空間，防止最優(yōu)個體丟失，提高了ELM模型的估測精度和魯棒性。同時，對比不同輸入類型建?？芍?，以特征波段組合作為輸入，模型估測效果較好，驗證集決定系數(shù)R2為0.849 7，均方根誤差RMSE為0.301 1。