東北寒地水稻莖稈纖維素含量近紅外光譜反演

徐 博， 許童羽, 2*， 于豐華, 2， 張國圣， 馮 帥， 郭忠輝， 周長獻(xiàn)

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

引 言

倒伏是影響農(nóng)作物產(chǎn)量的一個(gè)重要因素。 水稻莖稈細(xì)胞壁組成(木質(zhì)纖維素)極大地影響莖稈的機(jī)械強(qiáng)度并決定其抗倒伏能力[1]。 在水稻抗倒伏育種中， 莖稈纖維素含量檢測一般采用傳統(tǒng)的化學(xué)方法， 費(fèi)時(shí)費(fèi)力， 成本較高[2]， 而近紅外光譜檢測具有高效、 無損等優(yōu)點(diǎn)， 已廣泛應(yīng)用于農(nóng)業(yè)領(lǐng)域的各項(xiàng)指標(biāo)檢測。 在利用近紅外光譜對(duì)水稻莖稈纖維素含量進(jìn)行定量檢測的研究中， Diallo[3]利用近紅外光譜結(jié)合競爭性自適應(yīng)加權(quán)(CARS)篩選特征波段建立偏最小二乘回歸模型(PLSR)對(duì)水稻莖稈纖維素等元素含量進(jìn)行預(yù)測， 模型測試集RMSE為25.4 mg·g-1。 Huang[4]等在利用近紅外光譜法檢測轉(zhuǎn)基因水稻生物質(zhì)糖化壁聚合物， 建立了水稻莖稈纖維素含量的偏最小二乘(PLSR)回歸預(yù)測模型， 模型的RPD為2.79。 Huang[5]等利用近紅外光譜結(jié)合標(biāo)準(zhǔn)多元散射校正(SMSC)與偏最小二乘回歸模型(PLSR)對(duì)水稻莖稈纖維素含量快速估計(jì)， 模型的RPD為2.36。 上述研究表明利用近紅外光譜對(duì)水稻莖稈纖維素含量定量檢測的可行性， 但所采用的預(yù)處理和光譜特征選擇方法得到的光譜特征變量中仍存在一定的噪音和冗余信息； 鑒于此， 本研究在光譜預(yù)處理和特征波長選擇中進(jìn)行了一定的研究， 采用SNV-CWT方法對(duì)光譜進(jìn)行預(yù)處理， 降低由散射所導(dǎo)致的光譜噪聲， 采用SiPLS-IRIV特征波段選擇方法有效的去除了光譜中存在的冗余信息， 并建立了支持向量機(jī)回歸(ε-support vectormachine regression, εSVR)和核極限學(xué)習(xí)機(jī)(kernel-based extreme learning machine， KELM)的水稻莖稈纖維素含量反演模型， 以期為水稻莖稈纖維素含量反演提供一定的科學(xué)依據(jù)。

1 實(shí)驗(yàn)部分

1.1 樣本獲取與農(nóng)學(xué)參數(shù)測定

試驗(yàn)樣本于2018年8月—9月在遼寧省沈陽市沈北新區(qū)清水臺(tái)鎮(zhèn)柳條河村(123°63′E， 42°01′N)試驗(yàn)田獲取。 試驗(yàn)田總面積為1 333 m2， 等分4個(gè)試驗(yàn)小區(qū)， 進(jìn)行4個(gè)梯度的氮肥處理， 分別為0， 50， 100和150 kg·hm-2， 試驗(yàn)水稻品種為秋光， 種植密度為13.3 cm×30 cm。 在水稻灌漿期至成熟期共進(jìn)行15次田間采樣， 每次對(duì)各個(gè)小區(qū)進(jìn)行4次重復(fù)采樣， 截取水稻莖稈基部倒數(shù)2和3節(jié)進(jìn)行殺青烘干處理， 經(jīng)粉碎機(jī)粉碎成可通過篩孔為1 mm(18目)的樣本， 樣本的農(nóng)學(xué)參數(shù)為水稻莖稈纖維素含量值， 采用VanSoest干重法測定， 剔除異常含量值后， 共計(jì)234個(gè)樣本。 采用Kennard-Stone(KS)方法， 按照7∶3的比例將樣本集劃分為訓(xùn)練集(163)和測試集(71)， 如表1所示， 訓(xùn)練集和測試集的莖稈纖維素含量范圍分別為141.86～399.39和155.46～351.16 mg·g-1， 測試集莖稈纖維素含量在訓(xùn)練集范圍內(nèi)， 標(biāo)準(zhǔn)差相差不大， 樣本分散度基本一致， 訓(xùn)練集和測試集的樣本具有較好的代表性。

表1 訓(xùn)練集和測試集統(tǒng)計(jì)分析

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

水稻莖稈近紅外光譜數(shù)據(jù)采用OceanOptics公司生產(chǎn)的NIRQuest512型號(hào)非成像高光譜儀于實(shí)驗(yàn)室內(nèi)獲取。 光譜采集系統(tǒng)由HL-2000系列鹵鎢光源， 適用于可見光近紅外波段(360～2 000 nm)， SLIT-10光柵分光系統(tǒng)， 濱松銦鎵砷化物(InGaAs)陣列探測器和QR400-7-VIS-NIR反射探頭組成。 光譜波段范圍為900～1 700 nm， 采樣間隔為1.5 nm， 共512個(gè)波段。 光譜數(shù)據(jù)采集時(shí)， 光源預(yù)熱15 min達(dá)到熱平衡和穩(wěn)定后， 關(guān)閉快門， 進(jìn)行暗噪聲測量； 將反射探頭垂直固定在試驗(yàn)臺(tái)架上， 使用標(biāo)準(zhǔn)漫反射板測量反射光譜作為標(biāo)準(zhǔn)漫反射， 然后將反射探頭與水稻莖稈貼合； 設(shè)定光譜掃描次數(shù)為10次， 每個(gè)樣本重復(fù)測量3次取平均光譜作為該樣本的光譜反射率。

盡管使用多次掃描的平均光譜作為樣本光譜， 但由于反射探頭與光譜儀采用6繞1光纖束連接， 測量過程中， 難以避免光纖位置變化， 導(dǎo)致部分光纖內(nèi)的光會(huì)因散射損失帶來一定的光譜誤差； 同時(shí)由于水稻莖稈為柱形結(jié)構(gòu)且相對(duì)粗糙， 存在一定的表面散射， 光譜噪聲仍不可避免。 為了提高模型預(yù)測能力， 采用標(biāo)準(zhǔn)變量正態(tài)變換(standard normal variate， SNV)、 連續(xù)小波變換(continuous wavelet transform, CWT)[6]及兩種方法結(jié)合進(jìn)行光譜去噪。 為進(jìn)一步減少模型輸入量， 去除光譜中的冗余信息， 采用聯(lián)合區(qū)間偏最小二乘法(synergy interval PLS， SiPLS)、 迭代保留信息變量法(iteratively retaining informative variables， IRIV)[7]及兩種方法結(jié)合進(jìn)行光譜特征提取， 提取的光譜特征變量作為模型輸入。

1.3 模型參數(shù)優(yōu)化與模型建立

選用引入不敏感函數(shù)的支持向量機(jī)(εSVR)和核極限學(xué)習(xí)機(jī)(KELM)兩種方法進(jìn)行建模， 模型均采用核函數(shù)作為非線性特征映射函數(shù)， 核函數(shù)的選擇和模型參數(shù)確定直接影響模型的精度， 因此選擇適用性和穩(wěn)定性較好的徑向基函數(shù)(radial basis function， RBF)作為核函數(shù)， 采用灰狼算法(grey wolf optimizer， GWO)[8-9]， 差分進(jìn)化灰狼算法(differential evolution grey wolf optimizer， DEGWO)[10]， 自適應(yīng)差分進(jìn)化灰狼算法(self-adaptive differential evolution grey wolf optimizer， SaDEGWO)三種優(yōu)化算法確定模型的懲罰系數(shù)C， 核函數(shù)系數(shù)γ和不敏感參數(shù)ε(εSVR)。

1.4 模型評(píng)價(jià)標(biāo)準(zhǔn)

2 結(jié)果與討論

2.1 光譜預(yù)處理對(duì)比分析

為降低原始光譜中的噪聲， 提升光譜對(duì)莖稈纖維素含量的解析性， 分別采用了SNV， CWT和SNV-CWT對(duì)原始光譜進(jìn)行處理。 如圖1所示， 所有樣本的原始光譜曲線趨勢基本一致， 光譜反射特征具有一定差異， 與原始光譜對(duì)比明顯看出經(jīng)SNV處理后的光譜反射特征差異更為明顯， 光譜分布更為集中。 CWT處理光譜時(shí)， 選擇雙正交， 緊支撐和對(duì)稱性， 且具有線性相位特性的bior1.1小波函數(shù)作為基函數(shù)， 分別在21， 22， 23， 24， 25， 26， 27， 28， 29， 210尺度下對(duì)光譜進(jìn)行分解， 得到對(duì)應(yīng)尺度的小波系數(shù)， 尺度因子越小， 小波系數(shù)越小， 對(duì)應(yīng)的頻率越高， 具有較好的時(shí)間(波長)分辨率但頻率分辨率較差， 尺度因子越大則相反。 通常噪聲為高頻信號(hào)， 分解尺度越大， 小波系數(shù)越大， 頻率越低， 降噪性能越好， 但時(shí)間(波長)分辨率降低， 使有效信息特性被消減， 因此需要確定合適的分解尺度， 保證較高的頻率分辨率同時(shí)具有較好的時(shí)間(波長)分辨率。 如圖1(c)和(d)分別為對(duì)單個(gè)樣本的原始光譜和經(jīng)SNV處理后的光譜進(jìn)行不同尺度分解得到的小波系數(shù)曲線， 8～10尺度的光譜特征曲線起伏較大， 具有較好的頻率分辨率， 但時(shí)間(波長)分辨率低， 1～4尺度的光譜特征曲線起伏變化較小， 具有較好的時(shí)間(波長)分辨率， 但頻率分辨率較差， 而5～7尺度的光譜特征曲線起伏相似， 具有較好的頻率分辨率和時(shí)間(波長)分辨率， 為確定最佳分解尺度， 對(duì)不同分解尺度得到的小波系數(shù)與水稻莖稈纖維素含量進(jìn)行了相關(guān)性分析， 如圖2(a)和(b)所示， 從7尺度到10尺度相關(guān)性較低， 呈遞減趨勢， 而1尺度到6尺度相關(guān)性呈現(xiàn)遞增趨勢且相關(guān)性系數(shù)大于0.7的光譜區(qū)域逐漸增寬， 表明隨著分解尺度的增大， 有效去除高頻噪聲， 適當(dāng)增大分解尺度可有效的發(fā)掘光譜中隱含信息， 當(dāng)分解尺度過大， 導(dǎo)致有效信息被消除。 因此確定CWT和SNV-CWT的最佳分解尺度為6。

圖1 樣本的光譜曲線

圖2 小波能量系數(shù)與纖維素的相關(guān)系數(shù)

為了進(jìn)一步探究預(yù)處理方法的有效性， 分別對(duì)原始全光譜、 SNV處理光譜、 CWT處理光譜和SNV-CWT處理光譜采用10折偏最小二乘法建立水稻莖稈纖維素含量的定量分析模型。 結(jié)果如表2所示。 采用經(jīng)過SNV-CWT-bior1.1-26預(yù)處理后的全光譜建立的模型的測試集RPD為2.71， 高于未處理建立的模型(RPD=2.34)、 SNV處理的模型(RPD=2.60)及CWT-bior1.1-26處理建立的模型(RPD=2.55)。 表明采用SNV-CWT-bior1.1-26預(yù)處理方法能有效的減少光譜中的噪聲， 提升了光譜對(duì)莖稈纖維素含量的解析能力， 因此確定該方法作為光譜預(yù)處理方法。

表2 基于不同預(yù)處理方法的結(jié)果比較

2.2 光譜特征選擇方法對(duì)比與分析

由于全波段光譜數(shù)據(jù)中存在一定的無用和冗余信息， 影響建模的精度， 因此首先采用SiPLS方法將全波段512個(gè)光譜變量按16個(gè)波段等分成32個(gè)區(qū)間， 設(shè)定4個(gè)區(qū)間組合， 依據(jù)交叉驗(yàn)證選取RMSE最小值確定最優(yōu)聯(lián)合區(qū)間， 最終第12， 17， 23和24區(qū)間組合為最優(yōu)聯(lián)合區(qū)間， 區(qū)間對(duì)應(yīng)波長范圍分別為1 182～1 207， 1 310～1 334， 1 462～1 486和1 487～1 511 nm， 所選區(qū)間在全波段光譜中的位置如圖3所示。 SiPLS方法對(duì)全波段光譜按照預(yù)設(shè)區(qū)間值進(jìn)行劃分， 以各區(qū)間排列組合方式， 確定特征波長的范圍， 一定程度上去除了無用信息和干擾信息， 然而受區(qū)間值的限制， 難以避免優(yōu)選區(qū)間中仍然存在冗余信息。 因此， 采用IRIV算法對(duì)優(yōu)選區(qū)間進(jìn)一步進(jìn)行特征波段篩選， 全波段光譜經(jīng)過SiPLS方法將光譜變量數(shù)由512個(gè)降至64個(gè)， 再經(jīng)過IRIV進(jìn)行3次迭代， 變量消減至17個(gè)， 反向消除11個(gè)， 最終得到6個(gè)特征波長， 特征波長為1 200, 1 207, 1 325, 1 470, 1 482和1 492 nm， 如圖3所示， 依據(jù)前人研究纖維素的貢獻(xiàn)波段范圍主要分布在1 140～1 200 nm的C—H伸縮振動(dòng)的第二倍頻帶， 1 350～1 470 nm的C—H和O—H的第一倍頻帶， 1 490～1 600 nm的O—H的第一倍頻帶[12]。 本研究的結(jié)論基本一致， 但因受氫鍵效應(yīng)、 振動(dòng)耦合因素影響， 使得纖維素在近紅外吸收峰的位置存在部分位移。 為了進(jìn)一步說明SiPLS和IRIV方法組合的合理性， 采用IRIV方法對(duì)全波段光譜進(jìn)行特征波長選擇， 結(jié)果如圖3所示， 512個(gè)光譜變量經(jīng)4次迭代后， 變量縮減至30個(gè)， 再經(jīng)過反向消除14個(gè)， 得到最終特征波長16個(gè)(987， 1 033， 1 138， 1 141， 1 157， 1 217，1 285， 1 287， 1 296， 1 330， 1 336， 1 480， 1 481， 1 483， 1 485和1 508 nm)， 說明IRIV方法能夠有效的去除全波段光譜數(shù)據(jù)中的無用信息和干擾信息。 對(duì)比SiPLS-IRIV和IRIV方法發(fā)現(xiàn)， SiPLS-IRIV選取的特征波長更少， 波長間的共線性程度相對(duì)更低。

圖3 基于SiPLS-IRIV算法的光譜特征變量優(yōu)選

為了探究特征選擇方法的有效性， 分別對(duì)SiPLS， IRIV和SiPLS-IRIV特征選擇方法采用10折偏最小二乘法建立水稻莖稈纖維素含量的定量分析模型。 結(jié)果如表3所示。 采用SiPMS-IRIV特征選擇方法選取的6個(gè)特征波長建立的模型的測試集RPD為2.79， 相對(duì)全光譜(RPD=2.71)、 SiPLS特征選擇方法(RPD=2.72)及IRIV特征選擇方法(RPD=2.75)選取的特征波長建立的模型測試集PRD均有所提升， 同時(shí)變量壓縮率高達(dá)98.83%， 表明在保證模型預(yù)測精度前提下， SiPLS-IRIV特征選擇方法有效簡化了建模的復(fù)雜性。

表3 基于不同特征變量選擇方法的結(jié)果比較

2.3 模型對(duì)比與分析

表4 基于不同優(yōu)化算法的模型參數(shù)選擇結(jié)果比較

圖4 GWO, DEGWO和SaDEGWO優(yōu)化εSVR(a, b, c)和KELM(d, e, f)模型的訓(xùn)練和測試結(jié)果

在構(gòu)建基于高光譜水稻莖稈纖維素含量的定量反演模型過程中發(fā)現(xiàn)SNV-CWTbior1.1-26光譜預(yù)處理方法使得模型精度提升效果較好， 因光譜采集過程中， 受環(huán)境等不可避免因素影響， 造成光譜中含有一定噪聲， 經(jīng)SNV處理消除噪聲的光譜通過連續(xù)小波變換在時(shí)域和頻域進(jìn)行不同尺度展開， 從光譜信號(hào)中有效提取重要信息， 進(jìn)而使得模型精度提升。 采用未去噪聲的光譜進(jìn)行連續(xù)小波變換， 分解得到光譜信息中可能會(huì)保留一定噪聲信息。

在光譜特征篩選過程中， 發(fā)現(xiàn)采用IRIV方法提取光譜特征變量作為模型輸入， 精度高于SiPLS方法， 由于IRIV方法能夠有效的去除光譜中無信息變量和干擾信息， 而SiPLS方法提取的聯(lián)合區(qū)間仍然包含部分無用信息和干擾信息； 同時(shí)發(fā)現(xiàn)SiPLS-IRIV優(yōu)于IRIV， 原因在于IRIV方法在全波段范圍進(jìn)行特征篩選時(shí)， 相對(duì)變量數(shù)明顯增多， 存在一定的變量辨識(shí)度降低的問題， 同時(shí)表明采用IRIV對(duì)SiPLS聯(lián)合區(qū)間再進(jìn)行光譜特征提取方法的可行性， 進(jìn)一步表明適當(dāng)縮小IRIV變量篩選范圍一定程度上提升該方法的有效性。

最后從模型對(duì)比發(fā)現(xiàn)， 采用GWO， DEGWO及SaDEGWO優(yōu)化KELM模型的精度基本一致， 且略高于PLSR線性回歸模型， 采用SaDEGWO優(yōu)化εSVR模型精度均高于其他模型， 原因在于在解決本研究問題中， εSVR模型對(duì)徑向基函數(shù)的高維映射解析力更強(qiáng)， 適用性更好。

此外， 在采用SaDEGWO優(yōu)化建立εSVR模型時(shí)， 收斂速度相對(duì)較慢， 未來將在變異策略進(jìn)行改進(jìn)， 并考慮變異因子和交叉概率因子整體分布， 及個(gè)體之間的變異因子擾動(dòng)。

3 結(jié) 論