改進(jìn)隨機(jī)蛙跳算法在大豆品種快速鑒別中的應(yīng)用

李 偉, 譚 峰, 張 偉, 高陸思, 李金山

1. 黑龍江八一農(nóng)墾大學(xué)工程學(xué)院, 黑龍江 大慶 163319 2. 黑龍江八一農(nóng)墾大學(xué)信息與電氣工程學(xué)院, 黑龍江 大慶 163319 3. 黑龍江省農(nóng)業(yè)科學(xué)院綏化分院, 黑龍江 綏化 152052 4. 大慶市綠色農(nóng)產(chǎn)品監(jiān)測中心, 黑龍江 大慶 163311

引 言

大豆是中國重要的油料和糧食作物, 在國家的農(nóng)業(yè)生產(chǎn)、 糧食安全以及社會發(fā)展等領(lǐng)域擁有重要的戰(zhàn)略地位[1]。 黑龍江省作為中國大豆的主產(chǎn)區(qū), 具有得天獨(dú)厚的地理環(huán)境和生態(tài)氣候, 區(qū)位優(yōu)勢突出, 肩負(fù)著推進(jìn)大豆產(chǎn)業(yè)振興的重任。 隨著大豆雜交技術(shù)的廣泛應(yīng)用, 種子品類繁多, 在大豆的種植、 加工、 儲存以及流通的過程中極易發(fā)生混雜, 直接導(dǎo)致種子質(zhì)量下降, 體現(xiàn)不出特性優(yōu)勢, 給農(nóng)民和國家利益造成很大損失。 因此, 快速準(zhǔn)確地鑒別大豆品種對于鑒定種子品質(zhì), 凈化種業(yè)市場以及保障糧食安全具有重要意義。

目前, 農(nóng)作物品種鑒別的傳統(tǒng)方法有形態(tài)分析法[2]、 化學(xué)鑒定法[3]和分子標(biāo)記法[4]等。 在這些方法中, 形態(tài)分析法最為簡便, 但需要鑒定者對農(nóng)作物不同品種的外部形態(tài)和內(nèi)部結(jié)構(gòu)有豐富的認(rèn)知, 往往受主觀因素影響使得鑒別準(zhǔn)確度不高。 化學(xué)鑒定法、 分子標(biāo)記法鑒別準(zhǔn)確度能夠滿足需求, 但是存在檢測周期長、 過程繁瑣以及技術(shù)成本高等缺點(diǎn), 不適宜對樣本進(jìn)行快速批量分析。 而拉曼光譜是一種基于物質(zhì)分子的散射光譜, 具有“指紋”識別特性, 可以通過譜峰位移、 強(qiáng)度和譜線走勢分析物質(zhì)的組分及其結(jié)構(gòu)變化信息[5]。 與傳統(tǒng)鑒別方法相比, 拉曼光譜具有分析速度快、 準(zhǔn)確性高、 無需試樣制備以及綠色環(huán)保等優(yōu)點(diǎn), 在農(nóng)作物品種的快速鑒別中已經(jīng)得到廣泛應(yīng)用。 Dibs等[6]采集了玉米轉(zhuǎn)基因和非轉(zhuǎn)基因的拉曼光譜數(shù)據(jù), 采用線性判別分析(LDA)方法建立品種鑒別模型。 實(shí)驗(yàn)結(jié)果表明, Whittaker濾波器、 標(biāo)準(zhǔn)正態(tài)變換(SNV)和卷積平滑組合(savitzky-golay, SG)卷積平滑組合預(yù)處理方法效果最好, 再結(jié)合遺傳算法(GA)進(jìn)行特征變量, 所建模型預(yù)測準(zhǔn)確率最高已達(dá)到87.5%。 沙敏等[7]采集糯米、 秈米以及粳米的拉曼光譜數(shù)據(jù)共計72份, 采用小波變換以及歸一化對光譜數(shù)據(jù)進(jìn)行預(yù)處理, 再結(jié)合主成分分析(PCA)、 層次聚類分析(HCA)以及采用支持向量機(jī)(SVM)三種不同模式識別方法進(jìn)行建模分析。 實(shí)驗(yàn)結(jié)果表明, 采用SVM方法大米品種平均識別率高達(dá)98.86%。 Liu等[8]采集了8個糯小麥品種的拉曼光譜數(shù)據(jù), 通過MSC、 SG平滑和二階導(dǎo)數(shù)組合方法預(yù)處理光譜數(shù)據(jù), 再采用主成分分析法(PCA)進(jìn)行分類判別, 其中校正集預(yù)測準(zhǔn)確率為94.4%, 驗(yàn)證集預(yù)測準(zhǔn)確率為94.6%。 然而使用拉曼光譜技術(shù)鑒別大豆品種的研究并不多見。

本研究利用拉曼光譜儀采集了黑龍江省4個高蛋白大豆品種(黑農(nóng)88、 黑農(nóng)98、 綏農(nóng)71以及綏農(nóng)76)的樣本數(shù)據(jù)。 為了去除光譜中冗余和無效波長點(diǎn)對模型預(yù)測精度的影響, 提出一種改進(jìn)隨機(jī)蛙跳(modified random frog, MRF)算法, 實(shí)現(xiàn)大豆品種分類特征信息的有效提取, 再通過偏最小二乘判別分析法(parital least squares discriminant analysis, PLS-DA)構(gòu)建大豆品種鑒別模型, 以期達(dá)到提高品種分類精度的目的。

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

1.1 儀器和軟件

拉曼光譜數(shù)據(jù)由奧譜天成(廈門)光電有限公司制造的便攜式拉曼光譜儀進(jìn)行采集, 儀器激發(fā)波長為785 nm, 測量范圍在200～3 400 cm-1區(qū)間, 分辨率為6.58 cm-1。 光譜強(qiáng)度穩(wěn)定性低于5%, 符合穩(wěn)定性標(biāo)準(zhǔn)要求。 數(shù)據(jù)處理軟件采用MathWorks 公司的Matlab R2017b 實(shí)現(xiàn)。

1.2 樣品獲取和光譜采集

實(shí)驗(yàn)選用黑龍江省2020年主要種植的4種高蛋白大豆種子為研究對象, 品種分別是黑農(nóng)88、 黑農(nóng)98、 綏農(nóng)71以及綏農(nóng)76, 由黑龍江省農(nóng)業(yè)科學(xué)院提供。 從每個品種中隨機(jī)選取40 粒成熟、 無損傷的種子作為該品種樣本集, 4個品種共計160個樣本用于光譜數(shù)據(jù)采集。 拉曼光譜儀采集數(shù)據(jù)時, 激光功率設(shè)為500 mW, 積分時間設(shè)為5 000 ms。 對于每個樣本連續(xù)進(jìn)行3次采集, 將平均光譜作為該樣本的原始光譜數(shù)據(jù)。

1.3 方法

1.3.1 隨機(jī)蛙跳算法

隨機(jī)蛙跳(random frog, RF)算法[9]借鑒可逆跳躍馬爾可夫鏈蒙特卡羅(reversible jump Markov Chain Monte Carlo, RJMCMC)算法的思想, 建立了一種數(shù)學(xué)上簡單并且計算效率高的特征波長選擇算法。 RF通過在模型空間中實(shí)現(xiàn)固定維度和跨維度移動搜索獲得一個偽MCMC鏈, 用于計算每個特征變量的被選概率, 作為該特征變量選擇的標(biāo)準(zhǔn), 主要包括以下三個步驟[9]:

(1) 參數(shù)初始化:K0為初始變量集F0的個數(shù), 在1～P之間隨機(jī)設(shè)置(P為全部變量個數(shù)); 迭代次數(shù)N為10 000;θ為正態(tài)分布控制參數(shù), 取默認(rèn)值0.3;ω為大于1的候選變量個數(shù)控制參數(shù), 取默認(rèn)值3;η為接受概率, 取默認(rèn)值0.1。

(2) 候選變量子集構(gòu)造: 通過N(K0, 0.3K0)產(chǎn)生一個隨機(jī)數(shù), 就近取整確定為候選子集變量個數(shù)K*, 進(jìn)而構(gòu)造一個包含K*個變量的候選子集F*,F為包含全部p個變量的集合。

① 若K*=K0, 則F*=F0。

② 若K*

③ 若K*>K0, 則先從F-F0中隨機(jī)選取3(K*-K0)個變量, 并與F0構(gòu)成新變量子集T后建立PLS-DA模型, 保留模型中回歸系數(shù)絕對值最大的K*個變量構(gòu)成候選子集F*。

候選變量子集F*初選后, 判斷其是否被接受進(jìn)行下一次迭代。 分別用F0和F*構(gòu)建偏最小二乘分析模型, Err0和Err*為兩個模型的預(yù)測誤差。 如果Err*≤Err0, 則接受F*作為F1, 否則以0.1Err0/Err*的概率接受F*作為F1。 用F1代替F0, 進(jìn)行下一次迭代計算, 直到結(jié)束N次循環(huán)。

(3) 變量選擇概率計算: 經(jīng)過N次迭代計算共獲取N個變量子集。Nj表示第j個變量在N次迭代中被選擇的次數(shù)。 則每個變量被選擇概率根據(jù)式(1)進(jìn)行計算

(1)

由式(1)可知, 變量越重要被隨機(jī)蛙跳算法選擇的概率就越大。 因此, 該被選概率可以作為衡量變量重要性的指標(biāo)。 通過設(shè)定不同閾值, 選取被選概率大于閾值的變量作為特征波長變量, 對比所建模型性能從而確定最優(yōu)特征波長。

1.3.2 改進(jìn)隨機(jī)蛙跳算法

在采用RF算法提取特征變量時, 由于初始變量集F0是從原始數(shù)據(jù)集中隨機(jī)抽取產(chǎn)生的, 可能會引入無用或干擾波長變量, 使得算法的預(yù)測能力和收斂速度降低。 因此, 為了提高初始變量集的有效性, 提出一種基于最小絕對收斂和選擇算子(least absolute shrinkage and selection operator, LASSO)的改進(jìn)隨機(jī)蛙跳算法。 LASSO[10]是一種基于普通偏最小二乘法的特征波長選擇算法, 通過增加L1范數(shù)懲罰, 在預(yù)測殘差平方和最小準(zhǔn)則下, 將模型的回歸系數(shù)的絕對值之和限定在一個常數(shù)范圍內(nèi), 從而將一些不重要變量的回歸系數(shù)嚴(yán)格收縮為0。 多元線性回歸分析模型如式(2)

Y=Xβ+ε

(2)

式(2)中,X=[x1,x2, …,xn]T∈Rn×m為大豆樣本組分變量,Y=[y1,y2, …,yn]T∈Rn×1為樣本對應(yīng)品種屬性, 其中,m為樣本的特征變量數(shù),n為大豆樣本數(shù)量。β=[β1,β2, …,βm]T∈Rm×1是回歸系數(shù)向量。 LASSO準(zhǔn)則的目標(biāo)函數(shù)如式(3)

(3)

式(3)中,λ1為L1正則化系數(shù), 當(dāng)λ1=0時為常規(guī)最小二乘法, 當(dāng)λ1逐漸增大時L1正則化影響也越大, 越來越多的系數(shù)將接近于0或者等于0,λ1的取值將根據(jù)RMSECV取最小值時確定。

LASSO算法可以篩選與屬性變量最相關(guān)特征變量, 采用該方法提取大豆重要特征波長點(diǎn)作為RF初始變量集F0, 可以有效消除初始變量的隨機(jī)性, 在此基礎(chǔ)上開始迭代計算, 可以改善RF算法所需迭代次數(shù)大、 算法收斂慢等問題。 同時, RF算法通常選取概率排序1～10的變量, 或者選取概率大于某一閾值的變量, 因此提取的特征波長往往具有不確定性。 改進(jìn)方法如下: 首先去除選擇概率為0的變量, 對于排序后變量以10個波長點(diǎn)為間隔, 每次增加1個間隔構(gòu)建波長變量與大豆品種屬性的偏最小二乘分析模型。 然后計算每個模型的RMSECV找到最優(yōu)特征波長以提高模型的預(yù)測性能。

1.3.3 ElaticNet特征選擇方法

ElaticNet方法結(jié)合了LASSO回歸和嶺回歸的思想, 在普通最小二乘回歸基礎(chǔ)上增加L1懲罰和L2懲罰, 可以同時對變量系數(shù)的絕對值和平方項(xiàng)進(jìn)行壓縮。 LASSO回歸對所有波長變量的回歸系數(shù)進(jìn)行同等程度的壓縮, 則可能由于過度壓縮導(dǎo)致模型的預(yù)測準(zhǔn)確度下降。 嶺回歸[11]是通過增加L2懲罰, 收縮不重要變量回歸系數(shù)以達(dá)到提高預(yù)測準(zhǔn)確度的目的, 但不具有變量選擇的能力。 Elastic Net結(jié)合這兩種懲罰函數(shù)的特點(diǎn), 在高維數(shù)據(jù)的情況下剔除冗余變量, 且一般不會過度壓縮回歸系數(shù), 從而達(dá)到提高模型預(yù)測精度的目的。 Elastic Net準(zhǔn)則的目標(biāo)函數(shù)如式(4)

(4)

令α=λ1/(λ1+λ2),λ=λ1+λ2, 則式(4)為

(5)

從式(5)中可以看出, 當(dāng)α=0時為嶺回歸, 當(dāng)α=1為LASSO回歸, 當(dāng)α∈(0, 1)時為Elastic Net回歸, 它是L1和L2懲罰的特殊線性組合, 采用α的大小來調(diào)節(jié)它們之間的權(quán)重。 當(dāng)λ=0時為普通最小二乘法, 當(dāng)λ逐漸增大時正則化的影響也將越大,λ的取值將根據(jù)RMSECV取最小值時確定。

1.3.4 模型構(gòu)建方法與評價

PLS是在主成分回歸基礎(chǔ)上提出的多元校正方法, 在計算主成分時, 考慮主成分方差盡量最大以提取有用信息, 同時還保證主成分與樣本屬性變量盡可能地相關(guān), 使得所建模型具有較高的預(yù)測精度和較強(qiáng)的穩(wěn)健性。 因此以改進(jìn)隨機(jī)蛙跳法選擇后的特征波長為輸入變量, 使用PLS-DA方法建立大豆品種鑒別模型。 并采用均方根誤差(root mean square error, RMSE)、 決定系數(shù)(coefficient of determination,R2)和準(zhǔn)確率(accuracy)三個指標(biāo)對模型效果進(jìn)行評價。 RMSE值越小, 樣本預(yù)測值與實(shí)際值之間的差異越小,R2越接近1, 樣本預(yù)測值與實(shí)際值之間的相關(guān)度越高, 表明模型的預(yù)測性能越好。 對于識別正確率, 分別將黑農(nóng)88、 黑農(nóng)98、 綏農(nóng)71以及綏農(nóng)76四個大豆品種樣本賦值 1、 2、 3、 4, 根據(jù)文獻(xiàn)[12], 將PLS-DA鑒別模型的誤差閾值設(shè)置為±0.5, 結(jié)果分別在1±0.5、 2±0.5、 3±0.5和4±0.5之間時, 大豆樣本識別正確, 否則識別錯誤。

2 結(jié)果與討論

2.1 樣本光譜數(shù)據(jù)分析

大豆樣本的原始拉曼光譜如圖1(a)所示。 標(biāo)準(zhǔn)正態(tài)變換(standard normal variate transformation, SNV)[13]是一種有效的預(yù)處理方法, 可以校正樣本不均勻帶來的散射效應(yīng)對光譜數(shù)據(jù)的影響, 預(yù)處理效果如圖1(b)所示。 再利用KS (Kennard-Stone)算法按照3∶1的比例劃分大豆樣本集, 則校正集樣本為120個、 預(yù)測集樣本為40個。 采用PCA方法對校正集和預(yù)測集的共160個大豆拉曼光譜數(shù)據(jù)進(jìn)行降維處理, 第一、 第二和第三主成分的貢獻(xiàn)率分別為88.56%、 4.65%和3.45%, 前三個主成分的累積貢獻(xiàn)率達(dá)到96.66%, 說明可以反映原始樣本所提供的絕大多數(shù)信息。 采用前三個主成分繪制大豆樣本校正集和預(yù)測集的得分圖如圖2所示, 可知校正集和預(yù)測集的樣本分布比較均勻, 校正集能夠覆蓋預(yù)測集, 說明該樣本劃分方法有效并可以用于建模分析。

圖1 大豆原始及預(yù)處理后拉曼光譜圖(a): 原始拉曼光譜; (b): SNV預(yù)處理后拉曼光譜Fig.1 Raw and SNV pre-processed Raman spectra of soybean(a): Raw Raman spectra; (b): SNV pre-processed Raman spectra

圖2 大豆樣本前三主成分得分圖Fig.2 Score plot of the top three PCA for soybean samples

2.2 拉曼光譜特征波長優(yōu)選

2.2.1 RF特征波長優(yōu)選

表1 不同閾值的優(yōu)選波長Table 1 Optimized wavelengths with different thresholds

圖3 RF算法的波長被選概率Fig.3 Selected probability of each wavelength by RF

2.2.2 MRF特征波長優(yōu)選

采用LASSO算法對大豆全光譜數(shù)據(jù)進(jìn)行初選。 正則化參數(shù)λ1通過10折交叉驗(yàn)證法確定, 首先將預(yù)處理后大豆樣本數(shù)據(jù)劃分為10份, 每次選取其中9份作為校正集, 剩余1份作為預(yù)測集, 采用LASSO方法選擇特征波長建模并計算每個模型均方根誤差, 經(jīng)過10次迭代計算后取其平均值作為最終RSMECV值。 當(dāng)λ1=0.003 0時, RMSECV最小為0.195 8, 此時共提取了106個與大豆品種屬性最相關(guān)的特征變量, 將其作為RF算法的初選變量集F0, 開始進(jìn)行迭代計算。 將MRF算法的迭代次數(shù)分別設(shè)置為1 000、 1 500、 2 000和2 500次, 運(yùn)算結(jié)果如表2所示。

表2 不同迭代次數(shù)的優(yōu)選波長Table 2 Optimized wavelengths with different iteration times

對于不同迭代次數(shù)選取的特征波長, 首先去除被選概率為0的波長點(diǎn), 對剩余波長進(jìn)行排序后以10個波長點(diǎn)為間隔, 每次增加1個間隔建立PLS-DA模型, 選擇RMSECV取最小值時的波長組合為優(yōu)選特征波長。 由表2可知, 當(dāng)?shù)螖?shù)為2 000次時, 不同波長組合的RMSECV計算結(jié)果如圖4(a)所示, RMSECV值最小時為0.111 0, 最終選取300個有效特征波長點(diǎn), 如圖4(b)所示。

圖4 MRF特征波長優(yōu)選結(jié)果(a): 特征波長組合的RMSECV值; (b): 優(yōu)選特征波長分布Fig.4 Results of MRF characteristic wavelength optimization(a): RMSECV of characteristic wavelength combinations; (b): Distributionof optimized characteristic wavelength

2.2.3 ElasticNet特征波長優(yōu)選

由式(5)可知, ElasticNet的優(yōu)化函數(shù)L(β)包含參數(shù)α(0<α<1)和正則化系數(shù)λ(λ>0)。 通過10折交叉驗(yàn)證法確定調(diào)節(jié)參數(shù)α, 由表3可知, 當(dāng)參數(shù)α取0.016時, RMSECV值最小為0.145 6, 此時最優(yōu)正則化系數(shù)為λ=0.031 7。 基于所選最優(yōu)參數(shù)(α=0.016,λ=0.031 7), ElasticNet法通過坐標(biāo)下降法(coordinate descent)迭代計算3201個大豆拉曼光譜數(shù)據(jù)的稀疏系數(shù), 非零稀疏系數(shù)對應(yīng)成分變量即為所選特征波長點(diǎn), 通過ElasticNet法共選擇了1 345個大豆特征波長點(diǎn), 波長分布如圖5所示。

表3 調(diào)節(jié)參數(shù)α選擇Table 3 Selection of adjustment parameter α

圖5 ElasticNet選擇的特征波長分布Fig.5 Characteristic wavelengths distribution by ElasticNet

2.3 模型構(gòu)建結(jié)果評價與分析

將預(yù)處理后大豆光譜數(shù)據(jù), 分別采用RF、 MRF、 LASSO和ElasticNet算法進(jìn)行特征波長提取, 將這4種方法得到的特征波長點(diǎn)以及全光譜數(shù)據(jù)作為預(yù)測模型輸入變量, 大豆品種屬性作為輸出變量, 建立大豆品種鑒別的PLS-DA模型, 各模型鑒別結(jié)果如表4所示。

表4 不同波長選擇方法的模型對比Table 4 Model comparisons of different wavelength selection methods

3 結(jié) 論

提出的MRF特征波長優(yōu)選方法具有高效的全局搜索能力, 可以有效改善RF算法初始變量集隨機(jī)性、 所需迭代次數(shù)大、 閾值選取不確定的問題。 利用該方法結(jié)合拉曼光譜建立黑龍江省高蛋白大豆品種的PLS鑒別模型, 并與全光譜以及RF、 LASSO和ElasticNet算法提取的特征波長建模結(jié)果進(jìn)行對比分析。 通過MRF波長選擇方法, 可以有效地剔除無關(guān)波長點(diǎn), 簡化模型結(jié)構(gòu), 而且模型的預(yù)測性能得到一定程度地提升, 識別準(zhǔn)確率達(dá)到100%, 為黑龍江省大豆在種植、 加工、 儲存與流通過程中的品種快速鑒別提供了一種新思路。 在后續(xù)研究中, 應(yīng)進(jìn)一步增加不同產(chǎn)地、 不同品種的大豆樣本數(shù)量, 以提高鑒別模型的穩(wěn)健性和泛化能力, 以期在實(shí)際應(yīng)用中取得良好效果。