張小鳴　冒智康　李紹穩(wěn)　金秀　朱正偉

摘要：光譜變量多、光譜信息易重疊、數(shù)據(jù)存在冗余且存在大量噪聲，導(dǎo)致近紅外光譜分析建模的時(shí)間長、模型的泛化能力不強(qiáng)、預(yù)測精度不高等問題。為此，提出了一種基于變量有效性精英蟻群系統(tǒng)（EAS）的波長選擇算法，通過蟻群搜索選出特征波長。將變量有效性作為蟻群的初始信息素，在信息素更新過程中引入變量重要性投影方法，加快算法的收斂速度，避免陷入局部最優(yōu)。為了驗(yàn)證算法的有效性，選取193個(gè)沙姜黑土土壤樣本的近紅外光譜數(shù)據(jù)為研究對(duì)象，采用變量有效性精英蟻群算法選擇特征波長，構(gòu)建土壤速效磷含量偏最小二乘回歸預(yù)測模型，對(duì)全譜偏最小二乘模型、蟻群算法偏最小二乘模型進(jìn)行比較研究。結(jié)果表明，采用變量有效性精英蟻群算法選出的特征波長建立的回歸模型精度比普通蟻群算法選出的特征波長和全光譜高，建立回歸模型更加簡單。

關(guān)鍵詞：土壤速效磷;特征變量選擇;蟻群算法;變量有效性;精英策略

中圖分類號(hào)： TP391文獻(xiàn)標(biāo)志碼： A

文章編號(hào)：1002-1302（2019）19-00 -

收稿日期：2018-07-18

基金項(xiàng)目：農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)物聯(lián)網(wǎng)技術(shù)集成與應(yīng)用重點(diǎn)實(shí)驗(yàn)室開放基金（編號(hào)：2016KL07）。

作者簡介：張小鳴（1958—），男，安徽合肥人，博士，教授，主要從事嵌入式系統(tǒng)應(yīng)用研究。E-mail：mdd2942942@163.com。

近年來，過度耕種、施肥已經(jīng)嚴(yán)重影響了土壤氮、磷、鉀循環(huán)及其平衡[1]?？焖儆行y定土壤氮、磷、鉀含量，對(duì)指導(dǎo)農(nóng)業(yè)耕種中作物種植和施水施肥有一定的指導(dǎo)意義[2]。光譜檢測技術(shù)具有分析速度快、可多組分同時(shí)測定、非破壞性、低成本和操作簡單等特點(diǎn)，已被廣泛應(yīng)用于土壤養(yǎng)分含量檢測[3]。但是近紅外光譜存在嚴(yán)重的共線性，信噪比低，并且包含大量無關(guān)信息，因此特征波長選擇成為近紅外光譜分析中不可或缺的重要步驟，是光譜分析技術(shù)領(lǐng)域的研究熱點(diǎn)[4]。近年來，國內(nèi)外學(xué)者提出了大量特征變量選擇算法，主要包括基于偏最小二乘回歸（PLSR）模型參數(shù)的算法，基于連續(xù)投影策略的算法，基于波長區(qū)間選擇的算法，基于群智能優(yōu)化算法和基于模型集群分析策略等[5]，其中群智能算法成為研究熱點(diǎn)[6]。郭志明等使用普通蟻群算法優(yōu)選預(yù)測蘋果可溶性固形物含量的近紅外光譜特征波長，建立偏最小二乘模型，使得波長變量數(shù)減少，模型預(yù)測性能提高，決定系數(shù)和預(yù)測均方根誤差分別為 0.97 和0.51[7]。但是這種方法將蟻群初始信息素濃度都設(shè)為1，隨機(jī)性強(qiáng)，容易出現(xiàn)搜索停滯現(xiàn)象。Shamsipur等提出一種蟻群算法進(jìn)行特征波長選擇，但是波長需要人為設(shè)定，具有一定的不確定性[8]。基于以上問題，本研究提出基于變量有效性的精英蟻群算法，同時(shí)結(jié)合前向選擇方法優(yōu)選特征波長，采用精英蟻群搜索，并將變量有效性作為蟻群初始信息素，以避免蟻群陷入搜索停滯，搜索時(shí)間過長。最后，通過前向選擇實(shí)現(xiàn)特征波長自主選擇，以期提高土壤速效磷含量分析模型的準(zhǔn)確性和穩(wěn)定性。

1 土壤樣本采集與處理

試驗(yàn)所用的193個(gè)土壤樣本取自安徽省蒙城縣、宿州市埇橋區(qū)和懷遠(yuǎn)縣，這些地區(qū)的主要作物是小麥和玉米，土壤類型為沙姜黑土。根據(jù)作物種植密度、施肥濃度梯度的不同選取采樣點(diǎn)，以保證樣本速效磷含量存在一定差異。將采得的土壤樣本封存到實(shí)驗(yàn)室進(jìn)行風(fēng)干后碾磨，過20目篩。將處理后的土壤樣品等分為2份，一份用于標(biāo)準(zhǔn)定量測試，另一份用于光譜數(shù)據(jù)采集。速效磷含量測定采用碳酸氫鈉浸提-鉬銻抗分光光度法[9]。使用OFS1700地物光譜儀（海洋光學(xué)亞洲公司）和50 W鹵鎢燈接觸式反射探頭采集土壤近紅外光譜數(shù)據(jù)，測量時(shí)每個(gè)樣本旋轉(zhuǎn)90°測量6次取平均值作為光譜數(shù)據(jù);光譜范圍為350～1 655 nm，光譜分辨率為1 nm，共采集波長數(shù)1 306個(gè)。去除檢驗(yàn)值相差較大的異常樣本，最終選出188個(gè)樣本用于建模分析。

2 光譜數(shù)據(jù)預(yù)處理

188個(gè)土壤樣本在波長為350～1 655 nm處的原始近紅外光譜反射率如圖1-a所示，可以發(fā)現(xiàn)，光譜在400～500 nm處反射率較低，在500～700 nm處明顯上升;在 930 nm 處有明顯的波動(dòng)，可能是土壤中鐵氧化物（針鐵礦、赤鐵礦）對(duì)光譜吸收引起的[10];在1 400 nm處有明顯吸收峰，可能是由于土壤樣品中殘留水分的影響[11]。光譜圖像出現(xiàn)較多毛刺，伴有較多噪聲，可能是由于土壤顆粒大小、均勻度不同;光譜儀受高頻噪聲干擾、基線漂移的影響等，因此需要對(duì)光譜進(jìn)行預(yù)處理。Savitzky-Golay（S-G）平滑算法可消除基線漂移、傾斜等引起的噪聲。標(biāo)準(zhǔn)正態(tài)變量變換（SNV）和多元散射校正（MSC）可減少固體顆粒大小不同、光散射造成的影響。矢量歸一化可使數(shù)據(jù)保持相同標(biāo)度，減少數(shù)值差異對(duì)建模結(jié)果的影響[11]。根據(jù)比爾定律，土壤某養(yǎng)分的近紅外光譜吸光度與該養(yǎng)分的濃度成正比，因此將原始反射率光譜圖轉(zhuǎn)換成吸光度光譜圖，具體方式為將原始近紅外光譜反射率數(shù)據(jù)取倒數(shù)再取對(duì)數(shù)轉(zhuǎn)換為吸光度數(shù)據(jù)，然后采用S-G平滑算法+SNV對(duì)吸光度數(shù)據(jù)進(jìn)行預(yù)處理，得到188個(gè)土壤樣本在波長為350～1 655 nm近紅外光譜波段內(nèi)的吸光度光譜圖（圖1-b）。

3 光譜特征變量選擇與建模方法

3.1 校正集與驗(yàn)證集的劃分

采用濃度梯度法將樣本集按照速效磷含量的參照值進(jìn)行排序，以3 ∶1的比例劃分為校正集與驗(yàn)證集。由表1可知，校正集與驗(yàn)證集有相似的數(shù)據(jù)分布特征結(jié)構(gòu)，說明選出的校正集與驗(yàn)證集可以代表整體數(shù)據(jù)集的分布特征。

3.2 精英策略蟻群算法

3.2.1 蟻群算法波長選擇原理

螞蟻在覓食過程中，會(huì)在途經(jīng)線路上分泌信息素，后面的螞蟻會(huì)根據(jù)信息素判斷覓食方向。蟻群算法通過模擬螞蟻覓食行為，根據(jù)各路徑上的信息素沉積，蟻群最終收斂到最優(yōu)解，此最優(yōu)解即為特征波長變量。

3.2.2 精英策略蟻群算法

波長變量有效性可以表征變量對(duì)于模型預(yù)測指標(biāo)的重要性，其值越大，重要性越強(qiáng)。變量重要性計(jì)算方法如下：

采用蒙特卡羅法（MCS）采樣，從n個(gè)樣本中隨機(jī)抽取部分樣本建立偏最小二乘回歸（PLSR）模型，得到相應(yīng)的回歸系數(shù)向量（βi）、決定系數(shù)（R2），經(jīng)過m次采樣，第i個(gè)波長變量的有效性計(jì)算公式為

從公式中可以看出，回歸模型的決定系數(shù)越大、均方根誤差（RMSE）越小，該變量的有效性越大，重要性越強(qiáng)。

在模擬蟻群過程中，若沒有信息素，螞蟻就會(huì)隨機(jī)選擇一個(gè)方向，而且不會(huì)向著同一方向運(yùn)動(dòng)。引入精英蟻群思想并將變量有效性作為螞蟻初始信息素的濃度選擇依據(jù)，為蟻群提供初始方向。在信息素更新時(shí)引入變量重要性投影（VIP）系數(shù)，加快算法的收斂，最后通過前向選擇方法選出最終特征波長數(shù)量。

傳統(tǒng)蟻群算法計(jì)算時(shí)間長，且容易出現(xiàn)所有螞蟻求出的解完全一致的現(xiàn)象[12]?；诰⒉呗缘南伻盒畔⑺馗聶C(jī)制是根據(jù)重要性分配不同權(quán)重，將螞蟻按權(quán)重排序，避免算法陷入局部最優(yōu)，通過迭代實(shí)現(xiàn)全局尋優(yōu)[13]（圖2）。每次迭代將螞蟻所選波長變量代入偏最小二乘回歸模型，計(jì)算RMSE，RMSE越小，說明所選的波長變量精英度越高。

基于精英策略的蟻群信息素更新計(jì)算公式為

τi（t+1）=（1-ρ）×τi（t）+ρ×Δτ。（2）

式中：τi（t+1）為第i個(gè)螞蟻（第i個(gè)波長變量）在t+1時(shí)刻更新后的信息素含量;τi（t）為第i個(gè)螞蟻在t時(shí)刻的信息素含量;ρ為信息素衰減率常數(shù)，取值在0～1之間;Δτ為蟻群信息素增量。

蟻群信息素增量Δτ的計(jì)算公式為

Δτ=∑[DD（]mi=2[DD）]ωi×τi+λ×F。（3）

式中：ωi為第i個(gè)螞蟻信息素更新所占的權(quán)重;τi（t）為第i個(gè)螞蟻在t時(shí)刻的信息素含量;λ是信息素更新系數(shù)常量，取值為0～1;F為回歸模型的RMSE極小值的倒數(shù)。

F和ωi計(jì)算公式如下：

F=Q1+minRMSEωi=（1-λ）i。（4）

式中：Q是一個(gè)顯著因子常量，取值為0～1。若RMSE越小，則F值越大，Δτ值也越大，表示模型預(yù)測能力越好，螞蟻精英度越高。若ωi值越大，則Δτ值也越大，代表精英螞蟻分泌的信息素越多。Δτ值越大，表示蟻群算法更新速度越快。

為了進(jìn)一步加快蟻群算法收斂速度，減少計(jì)算時(shí)間，在基于精英策略的蟻群信息素基礎(chǔ)上，再引入VIP，計(jì)算公式如下：

VIPj=p∑w2k（SSYa/SSYt）SSY=∑[DD（]ni=1[DD）]（yi-y）2。（5）

式中：VIPj表示第j個(gè)波長變量的VIP得分;p表示波長變量數(shù)量;wk表示對(duì)應(yīng)的權(quán)重系數(shù);SSYa表示用a個(gè)潛變量建模時(shí)的SSY值;SSYt表示用全部潛變量建模時(shí)的SSY值。SSY表示每個(gè)土壤樣本速效磷含量預(yù)測值（yi）與所有土壤樣本速效磷含量預(yù)測值的平均值 y 之差的平方和;蟻群信息素改進(jìn)計(jì)算公式如下：

τ1=τ×VIP。（6）

式中：τ1為算法改進(jìn)后的蟻群信息素濃度;τ是蟻群原始信息素濃度;VIP表征波長變量對(duì)預(yù)測的重要性[14]。當(dāng)VIP>1時(shí)，τ1>τ，信息素被放大，波長變量被選中的概率變大。反之，信息素被縮小，波長變量被選中的概率變小。

最后，將通過精英策略蟻群算法選出的特征波長變量按選擇次數(shù)從大到小排列，并將排序后的特征波長變量依次代入偏最小二乘交互檢驗(yàn)?zāi)Ｐ陀?jì)算RMSE，根據(jù)采樣前向選擇（FS）規(guī)則進(jìn)行特征波長變量選擇，如果RMSE變小，則選擇此變量作為特征波長變量，否則刪除。

4 試驗(yàn)結(jié)果與分析

4.1 普通蟻群算法+PLSR

使用MATLAB R2014a進(jìn)行光譜預(yù)處理和PLS回歸建模分析。以1 306個(gè)（3 50～1 655 nm）全波長變量作為選擇對(duì)象，根據(jù)經(jīng)驗(yàn)值與優(yōu)選策略，設(shè)置普通蟻群算法的初始種群數(shù)量為100，單個(gè)螞蟻尋找最大變量數(shù)為40，顯著性因子為0.01，衰減系數(shù)為0.65。對(duì)每次迭代選出的最優(yōu)波長變量組中波長變量的相應(yīng)貢獻(xiàn)加1，迭代完成后，根據(jù)波長變量的貢獻(xiàn)率，按照前向選擇規(guī)則選出最佳特征波長變量組合。將普通蟻群算法（ACO）的迭代次數(shù)分別設(shè)置為30、50、100次，對(duì)產(chǎn)生的特征波長變量數(shù)及其建立的PLSR模型與全波PLSR模型的評(píng)價(jià)參數(shù)進(jìn)行比較，如表2所示。

從表2可以看出，使用普通蟻群選出的特征波長可以提高模型精度。隨著迭代次數(shù)的增加，模型的校正集均方根誤差（RMSECV）逐漸降低，這是因?yàn)樗惴ㄖ袑⒛繕?biāo)函數(shù)設(shè)置為模型的均方根誤差。但驗(yàn)證集均方根誤差（RMSEP）有先變小后增大的趨勢，綜合各項(xiàng)指標(biāo)以及算法時(shí)間復(fù)雜度，選擇迭代次數(shù)為50次。運(yùn)行MATLAB R2014a軟件，得到基本蟻群算法運(yùn)行后波長變量累計(jì)貢獻(xiàn)率圖（圖3）和特征波長點(diǎn)篩選圖（圖4）。圖3中縱向直線代表每個(gè)波長被選擇的次數(shù)即貢獻(xiàn)率，波長被選擇次數(shù)越多說明該波長貢獻(xiàn)越大，代表該波長被選為特征波長的概率越大。以被選中的波長貢獻(xiàn)率為依據(jù)，通過前向選擇選出最終的特征波長;圖4中圓圈代表對(duì)應(yīng)波長被選擇為特征波長，可以清楚地看出最終選擇出的特征波長分布情況，同時(shí)可以看出，該方法選出的特征波長點(diǎn)數(shù)較多。

4.2 精英策略蟻群算法+PLSR

在變量有效性精英策略蟻群算法（EEAS）的波長變量選擇方法研究中，參數(shù)設(shè)置與普通蟻群算法相同，初始種群大小為100，最大變量數(shù)為40，顯著性因子為0.01，衰減系數(shù)為

0.65，迭代次數(shù)分別為30、50、100次。引入前向選擇規(guī)則，根據(jù)波長點(diǎn)貢獻(xiàn)率大小進(jìn)行再次篩選，最終得到最佳特征波長，圖5為通過MATLAB R2014a畫出的基于變量有效性精英蟻群算法的所有波長點(diǎn)累計(jì)貢獻(xiàn)率圖，圖中縱向直線表示精英蟻群算法運(yùn)行過程中波長被選擇的次數(shù)，即貢獻(xiàn)，可以清楚地看出哪些波長點(diǎn)被選中的頻次高，即該波長點(diǎn)被選為特征波長的概率大。圖6為通過MATLAB R2014a畫出的基于變量有效性精英蟻群算法的特征波長點(diǎn)篩選圖，圖中圓圈代表最終選擇出的特征波長，可以清楚了解最終選擇出的特征波長分布，與圖4相比，本算法所選出的特征波長明顯減少。其對(duì)應(yīng)偏最小二乘回歸（PLSR）模型評(píng)價(jià)參數(shù)如表3所示。與普通蟻群算法相比，迭代次數(shù)為30、50、100次時(shí)，模型的校正集均方根誤差均有降低，且模型的決定系數(shù)提高，說明基于變量有效性的精英蟻群算法可以更準(zhǔn)確選出與建模相關(guān)的特征變量。另外，從所選擇的波長點(diǎn)圖（圖6）可以看出，基于變量有效性的精英蟻群算法選擇出的特征波長明顯減少，波長點(diǎn)分布在 400～480 nm、550～750 nm之間以及950、1 400、1 600 nm 附近，與文獻(xiàn)[15-16]對(duì)土壤速效磷含量預(yù)測有重要作用的波長大致相同，證明了本算法的正確性。

由圖7可以看出，基于變量有效性的精英蟻群算法收斂速度明顯加快，而且建立的回歸模型的RMSECV更小。在相同的參數(shù)下通過EEAS特征波長算法建立的PLSR模型效果更好，選出的波長點(diǎn)少，可降低模型的復(fù)雜度。

從圖8可以發(fā)現(xiàn)，EEAS-PLS模型的點(diǎn)集中在y=x附近，其次是ACO-PLS模型，而全譜模型點(diǎn)比較分散。說明采用特征波長建立的速效磷含量預(yù)測模型優(yōu)于全譜PLS模型，EEAS選出特征波長建立模型效果優(yōu)于ACO選出特征波長建立的模型。所以使用基于變量有效性的精英蟻群算法選擇特征波長，可以有效去除近紅外光譜中的光譜信息重疊、冗余，降低模型復(fù)雜度、節(jié)省時(shí)間，同時(shí)提高模型預(yù)測精度。

5 結(jié)論

采用變量有效性精英蟻群算法選擇的特征波長變量數(shù)少于普通蟻群算法。試驗(yàn)結(jié)果表明，利用特征波長變量建立的回歸模型預(yù)測效果優(yōu)于全譜波長變量?；谧兞坑行跃⑾伻核惴ǖ腜LS模型可以實(shí)現(xiàn)土壤速效磷含量的較高精度預(yù)測。

參考文獻(xiàn)：

[1]楊愛霞，丁建麗，李艷紅，等. 基于可見-近紅外光譜變量選擇的荒漠土壤全磷含量估測研究[J]. 光譜學(xué)與光譜分析，2016，36（3）：691-696.

[2]馬慧敏，陳 亮，潘康兵，等. 防磷固定劑對(duì)土壤有效性的影響[J]. 西南農(nóng)業(yè)學(xué)報(bào)，2017，30（11）：2533-2536.

[3]王樸杰，王世東，張合兵，等. 基于高光譜的復(fù)墾農(nóng)田土壤有機(jī)質(zhì)含量估測[J]. 土壤，2018，50（3）：558-565.

[4]林志丹，汪玉冰，王儒敬，等. 波長優(yōu)選對(duì)土壤有機(jī)質(zhì)含量可見光/近紅外光譜模型的優(yōu)化[J]. 發(fā)光學(xué)報(bào)，2016，37（11）：1428-1435.

[5]宋相中，唐 果，張錄達(dá)，等. 近紅外光譜分析中的變量選擇算法研究進(jìn)展[J]. 光譜學(xué)與光譜分析，2017，37（4）：1048-1052.

[6]賓 俊，范 偉，周冀衡，等. 智能優(yōu)化算法應(yīng)用于近紅外光譜波長選擇的比較研究[J]. 光譜學(xué)與光譜分析，2017，37（1）：95-102.

[7]郭志明，黃文倩，彭彥昆，等. 自適應(yīng)蟻群優(yōu)化算法的近紅外光譜特征波長選擇方法[J]. 分析化學(xué)，2014，303（40）：513-518.

[8]Shamsipur M，Zare-Shahabadi V，Hemmateenejad B，et al. Ant colony optimisation：a powerful tool for wavelength selection[J]. Journal of Chemometrics，2006，20：146-157.

[9]Bray R H，Kurtz L T. Determination of total，organic，and available forms of phosphorus in soils[J]. Soil Science，1945，59（1）：39-45.

[10]Viscarra Rossel R A，Bui E N，de Caritat P，et al. Mapping iron oxides and the color of Australian soil using visible-near-infrared reflectance spectra[J]. Journal of Geophysical Research：Earth Surface，2010，115：F4.

[11]Rossel R A V，Behrens T. Using data mining to model and interpret soil diffuse reflectance spectra[J]. Geoderma，2010，158（1/2）：46-54.

[12]江明珠，韓邦興，顏 暉，等. 便攜式近紅外光譜儀快速無損鑒別霍山石斛楓斗和河南石斛楓斗[J]. 林產(chǎn)化學(xué)與工業(yè)，2017，37（5）：102-106.

[13]吳碧霞. 一種基于精英策略的蟻群優(yōu)化QoS單播路由算法[J]. 莆田學(xué)院學(xué)報(bào)，2011，18（5）：73-76.

[14]邸 億，龍 飛，李卓越，等. 一種基于改進(jìn)蟻群算法的多目標(biāo)跟蹤數(shù)據(jù)關(guān)聯(lián)方法[J]. 計(jì)算機(jī)應(yīng)用與軟件，2013，30（4）：306-309.

[15]Mouazen A M，Maleki M R，De Baerdemaeker J，et al.On-line measurement of some selected soil properties using a VIS-NIR sensor[J]. Soil and Tillage Research，2007，93（1）：13-27.

[16]Daniel K W，Tripathi N K，Honda K. Artificial neural network analysis of laboratory and in situ spectra for the estimation of macronutrients in soils of Lop Buri （Thailand）[J]. Australian Journal of Soil Research，2003，41（1）：47-59.