汪六三，魯翠萍，王儒敬，黃　偉，郭紅燕，汪玉冰，林志丹，王　鍵，蔣　慶，宋良圖

(1.中國(guó)科學(xué)院合肥智能機(jī)械研究所，安徽合肥　230031;　2.合肥電子工程學(xué)院，安徽合肥　230037)

1　引　言

氮素是作物最重要的營(yíng)養(yǎng)元素，對(duì)絕大部分作物來(lái)說(shuō)，生長(zhǎng)所需的氮素90%以上來(lái)源于土壤［1］。分析測(cè)定土壤中各種形態(tài)的氮含量，對(duì)了解土壤氮的供給水平和指導(dǎo)施肥具有重要意義;在精準(zhǔn)農(nóng)業(yè)生產(chǎn)中，土壤不同形態(tài)氮素水平的測(cè)定是精準(zhǔn)施肥或平衡施肥的重要依據(jù)。土壤堿解氮是土壤氮素的重要指標(biāo)之一，堿解氮含量的高低能反映出近期土壤氮素豐缺。

傳統(tǒng)土壤堿解氮的測(cè)定是基于濕化學(xué)法［2］，操作復(fù)雜、測(cè)量時(shí)間長(zhǎng)、成本高，難以滿(mǎn)足快速檢測(cè)土壤堿解氮含量的需求。可見(jiàn)/近紅外光譜技術(shù)具有快速、無(wú)損、無(wú)污染等特點(diǎn)，它可以在短時(shí)間內(nèi)分析大量土壤樣品，實(shí)現(xiàn)土壤參數(shù)的快速測(cè)量，在土壤養(yǎng)分檢測(cè)中得到大量應(yīng)用［3-9］。通?？梢?jiàn)/近紅外光譜包括百千個(gè)波長(zhǎng)變量，而這些變量中含有與待測(cè)屬性無(wú)關(guān)的冗余噪聲［10］。因此，建模前有必要對(duì)波長(zhǎng)變量進(jìn)行選擇。這樣可以剔除不相關(guān)的變量，簡(jiǎn)化模型，更重要的是可以提高校正模型的預(yù)測(cè)能力。常用的變量選擇方法包括:無(wú)信息變量消除法(Uninformative variable elimination，UVE)［11］、遺 傳 算 法 (Genetic algorithms，GA)［12］、連續(xù)投影算法(Successive projections algorithm，SPA)［13］、競(jìng)爭(zhēng)性自適應(yīng)重加權(quán)算法(Competitiveadaptivereweightedsampling，CARS)［14］、隨機(jī)蛙跳(Random frog，RF) 等［15］。劉燕德等［16］以江西臍橙果園土壤為研究對(duì)象，通過(guò)遺傳算法、連續(xù)投影算法、競(jìng)爭(zhēng)性自適應(yīng)加權(quán)算法篩選出有機(jī)質(zhì)敏感波段，結(jié)果顯示競(jìng)爭(zhēng)性自適應(yīng)加權(quán)算法更加有效。賈生堯等［17］利用變量投影重要性(Variable importance in the projection，VIP)、無(wú)信息變量消除法選擇波長(zhǎng)變量，并分別用遞歸偏最小二乘回歸法(Recursive partial least squares regression，RPLS)、偏最小二乘回歸法(Partial least squares regression，PLS)建立土壤全氮和有機(jī)質(zhì)含量預(yù)測(cè)模型，結(jié)果表明VIP-RPLS對(duì)土壤全氮和有機(jī)質(zhì)含量具有更高的預(yù)測(cè)精度。高洪智等［18］將連續(xù)投影算法與貢獻(xiàn)值結(jié)合篩選土壤總氮的特征波長(zhǎng)，所建模型的預(yù)測(cè)精度優(yōu)于全譜偏最小二乘回歸結(jié)果。林志丹等［19］針對(duì)土壤有機(jī)質(zhì)含量預(yù)測(cè)，使用連續(xù)投影算法和遺傳算法進(jìn)行波長(zhǎng)優(yōu)選，主成分回歸構(gòu)建預(yù)測(cè)模型。結(jié)果顯示連續(xù)投影算法和遺傳算法都能夠有效減少參與建模的波長(zhǎng)數(shù)且提高模型預(yù)測(cè)精度。

本文以寧夏吳忠地區(qū)75個(gè)水稻土樣為研究對(duì)象，利用可見(jiàn)/近紅外光譜技術(shù)采集土壤樣品光譜，比較不同預(yù)處理方法和原始光譜對(duì)預(yù)測(cè)模型的效果，探討GA、SPA、CARS、RF 4種波長(zhǎng)選擇算法對(duì)土壤堿解氮含量預(yù)測(cè)模型的影響，為土壤堿解氮含量快速獲取和儀器開(kāi)發(fā)提供技術(shù)支持。

2　材料與方法

2.1　樣本采集、制備和劃分

在寧夏吳忠市葉盛貢米種植基地2號(hào)地(954畝)采集土壤樣品75個(gè)。土樣的采樣范圍為北緯38°07'23.65″～ 北緯 38°07'39.64″，東經(jīng)106°11'27.66″～ 東經(jīng) 106°12'03.28″;采集的土壤去除雜物后，讓其自然風(fēng)干。風(fēng)干后的土壤樣品先經(jīng)過(guò)人工碾磨，然后過(guò)篩。根據(jù)分析項(xiàng)目的要求，過(guò)60目數(shù)的篩子。每個(gè)樣品分兩份，一份用于實(shí)驗(yàn)室化學(xué)分析，另一份用于光譜分析，用自封袋對(duì)土樣進(jìn)行密封保存，并統(tǒng)一編號(hào)。土壤樣本共75個(gè)，使用SPXY(Sample set partitioning based on joint x-y distance)方法對(duì)樣本按4∶1進(jìn)行劃分，其中建模集60個(gè)，預(yù)測(cè)集15個(gè)樣本，如表1所示。

表1　土壤樣品堿解氮測(cè)量值分布Tab.1　Measurement distribution of soil available nitrogen mg·kg－1

2.2　土壤樣本光譜采集

本實(shí)驗(yàn)采用美國(guó)Veris公司的可見(jiàn)/近紅外光譜土壤檢測(cè)系統(tǒng)。該系統(tǒng)將光源和漫反射收集集成于一個(gè)犁頭內(nèi)。犁頭底部固定有藍(lán)寶石窗片。這種窗片在近紅外波段透過(guò)率大約為90%且足夠耐用，能夠承受與土壤的連續(xù)接觸。系統(tǒng)使用鹵素?zé)粽丈渫寥?，漫反射光通過(guò)光纖傳導(dǎo)進(jìn)入光譜儀。置于犁頭內(nèi)的光學(xué)快門(mén)每隔5 min自動(dòng)獲取暗光譜。光譜儀每秒大約采集20個(gè)光譜，同時(shí)通過(guò)USB接口傳輸給電腦保存。采集程序使用 Labview(National Instruments，Austin，TX，USA)編寫(xiě)。光纖傳導(dǎo)漫反射光進(jìn)入兩個(gè)分立的光譜儀。一個(gè)光譜儀(Ocean optics，USB4000)使用硅CCD作為探測(cè)器測(cè)量342～1 050 nm波段光強(qiáng)，光譜分辨率為3 nm;另一個(gè)光譜儀(Hamamatsu，C9914GB)使用InGaAs陣列探測(cè)器測(cè)量1 000～2 221.5 nm波段光強(qiáng)，光譜分辨率為8 nm。犁頭安裝于一個(gè)支架上，支架上安裝有光譜儀、控制器等光譜采集與控制處理系統(tǒng)。整個(gè)系統(tǒng)包含兩種測(cè)量模式——田間動(dòng)態(tài)測(cè)量模式和實(shí)驗(yàn)室靜態(tài)測(cè)量模式。本研究中，我們將犁頭從支架上拆卸下來(lái)，在實(shí)驗(yàn)室條件下采集土壤光譜。

為了使光源保持穩(wěn)定，鹵素?zé)艄庠粗辽兕A(yù)熱20 min左右。將過(guò)篩的土壤樣品裝入樣品杯，用樣品杯蓋將樣品杯內(nèi)的樣品壓實(shí)壓平。樣品杯貼近探頭窗口，采集土壤樣本光譜，轉(zhuǎn)動(dòng)樣品杯，每個(gè)樣品重復(fù)掃描3次，求其平均光譜曲線(xiàn)。

圖1　可見(jiàn)/近紅外土壤屬性測(cè)定儀結(jié)構(gòu)示意圖。其中:①藍(lán)寶石窗;②鹵素?zé)?③漫反射光收集探頭;④光纖;⑤光譜儀;⑥電源。

2.3　光譜預(yù)處理

樣本表面不均引起的散射現(xiàn)象、裝樣量不一致引起的光程差、暗電流和儀器隨機(jī)噪聲引起的光譜曲線(xiàn)不重復(fù)現(xiàn)象和基線(xiàn)漂移現(xiàn)象，以及樣本不同成分自檢相互干擾引起的背景因素和多重共線(xiàn)性等無(wú)用信息對(duì)光譜曲線(xiàn)均有影響［20］。為了達(dá)到較好的預(yù)測(cè)效果，通常需要對(duì)光譜數(shù)據(jù)進(jìn)行預(yù)處理。本研究分別對(duì)土壤樣本光譜數(shù)據(jù)進(jìn)行了Savitzky Golay平滑(SG smoothing)、多元散射校正(Multiple scatter correction，MSC)、標(biāo)準(zhǔn)正態(tài)變量變換(Standard normal variate，SNV)等預(yù)處理，并對(duì)各種預(yù)處理效果和原始光譜數(shù)據(jù)效果進(jìn)行了比較，具體見(jiàn)下文。光譜預(yù)處理使用Uncrambler X 10.4(CAMO，Norway)和 Matlab2012b(Math-Works，USA)軟件。

2.4　波長(zhǎng)篩選

原始光譜數(shù)據(jù)量大，同時(shí)存在冗余噪聲和大量共線(xiàn)信息，對(duì)光譜建模存在干擾，需要對(duì)波長(zhǎng)進(jìn)行篩選。本文對(duì)文中使用的遺傳算法(GA)、連續(xù)投影算法(SPA)、競(jìng)爭(zhēng)性自適應(yīng)重加權(quán)算法(CARS)和隨機(jī)蛙跳(RF)算法作一個(gè)簡(jiǎn)單介紹，具體的可以參見(jiàn)文獻(xiàn)［12-15］。

遺傳算法是模擬達(dá)爾文生物進(jìn)化論的自然選擇和遺傳學(xué)機(jī)理的生物進(jìn)化過(guò)程的計(jì)算模型，是一種通過(guò)模擬生物進(jìn)化隨機(jī)尋優(yōu)求解的常用算法。它利用選擇、交換和突變等算子的操作，隨著不斷的遺傳迭代，使目標(biāo)函數(shù)值較優(yōu)的變量被保留，較差的變量被淘汰，最終達(dá)到最優(yōu)結(jié)果。

連續(xù)投影算法是一種減少變量共線(xiàn)性問(wèn)題的前向循環(huán)選擇方法。利用向量的投影分析，該方法從一個(gè)波長(zhǎng)開(kāi)始，每次循環(huán)都計(jì)算它在未選入的波長(zhǎng)上的投影，將投影向量最大的波長(zhǎng)引入到波長(zhǎng)組合。每一個(gè)新選入的波長(zhǎng)，都與前一個(gè)線(xiàn)性關(guān)系最小。

競(jìng)爭(zhēng)性自適應(yīng)重加權(quán)算法是一種基于回歸系數(shù)進(jìn)行波長(zhǎng)選擇的方法。該方法模擬達(dá)爾文進(jìn)化理論的“適者生存”原則，將每個(gè)波長(zhǎng)看作一個(gè)體，對(duì)波長(zhǎng)進(jìn)行逐步淘汰，每次采樣過(guò)程中利用自適應(yīng)重加權(quán)采樣技術(shù)和指數(shù)衰減函數(shù)結(jié)合的方法優(yōu)選出PLS模型中回歸系數(shù)絕對(duì)值大的波長(zhǎng)點(diǎn)，去除PLS模型中回歸系數(shù)值權(quán)重較小的波長(zhǎng)，并計(jì)算變量子集RMSECV值，最后選擇最小RMSECV子集作為最優(yōu)變量子集。

隨機(jī)蛙跳是一種迭代方式工作的波長(zhǎng)選擇方法。首先一個(gè)變量子集被隨機(jī)選擇，基于這個(gè)選擇的變量子集按照一定的幾率進(jìn)行迭代，每迭代一次，變量子集更新一次。經(jīng)過(guò)N次迭代后，每個(gè)變量選擇的幾率被計(jì)算出來(lái)。然后幾率高的變量用來(lái)建立PLS模型，RMSECV最小的PLS模型對(duì)應(yīng)的變量為最優(yōu)變量。

本研究中，GA和SPA算法分別使用GAPLS工具箱、SPA工具箱實(shí)現(xiàn)［21］;CARS和RF算法使用LibPLS工具箱［22］實(shí)現(xiàn)。GAPLS工具箱要求初始變量小于200。因此首先對(duì)光譜數(shù)據(jù)使用三點(diǎn)平均后波長(zhǎng)數(shù)由380減為127，再使用GA進(jìn)行波長(zhǎng)篩選。

2.5　模型驗(yàn)證與評(píng)價(jià)

模型分別采用完全交互驗(yàn)證和外部驗(yàn)證對(duì)其性能進(jìn)行評(píng)價(jià)，由決定系數(shù)(R2)、建模均方根誤差(RMSECV)、預(yù)測(cè)均方根誤差(RMSEP)和預(yù)測(cè)相對(duì)分析誤差(RPD)來(lái)評(píng)價(jià)。在建模分析中，R2和RPD值高為好，同時(shí)RMSECV和RMSEP值低為好。當(dāng)RPD值位于1.5～2之間表示所建立的模型對(duì)分析指標(biāo)具有一定的預(yù)測(cè)能力，當(dāng)位于2～2.5之間表示所建立的模型對(duì)分析指標(biāo)定量預(yù)測(cè)是可行的，2.5～3之間表示所建立的模型具有較好的預(yù)測(cè)精度［3，23］。

3　結(jié)果與討論

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

75份土壤樣品的可見(jiàn)/近紅外光譜吸收?qǐng)D譜如圖2所示，在1 400 nm和1 900 nm附近有明顯的水分子吸收峰。

圖2　75份土壤樣品的可見(jiàn)/近紅外光譜

去除340～366 nm噪聲大的波段，選取366～2 221.5 nm波段光譜數(shù)據(jù)作為后續(xù)數(shù)據(jù)處理。對(duì)土壤光譜數(shù)據(jù)分別采用SG smoothing、SNV、MSC等預(yù)處理方法進(jìn)行光譜預(yù)處理，分別建立PLS模型。本文建立的所有預(yù)測(cè)模型由決定系數(shù)(R2)、建模均方根誤差(RMSECV)、均方根預(yù)測(cè)誤差(RMSEP)和相對(duì)預(yù)測(cè)偏差(RPD)評(píng)價(jià)，其中R2和RPD值越大，RMSECV、RMSEP值越小，模型性能越好。

表2　不同預(yù)處理方法的PLS建模結(jié)果Tab.2　Results of PLSmodels obtained with different spectral pretreatmentmethods

表2給出了不同預(yù)處理方法PLS建模結(jié)果。本次實(shí)驗(yàn)發(fā)現(xiàn)采集的光譜不經(jīng)過(guò)任何預(yù)處理效果最好，是由于光譜采集時(shí)儀器性能穩(wěn)定，制備的樣品過(guò)了60目篩，樣品的顆粒度比較小和均勻。

圖3　土壤光譜和GA變量篩選結(jié)果

表3　不同波長(zhǎng)篩選方法的PLS建模結(jié)果Tab.3　Results of PLSmodes obtained with differentwavelength selection

圖4　SPA運(yùn)行結(jié)果

3.2　波長(zhǎng)篩選

圖3顯示了基于GA算法的波長(zhǎng)篩選結(jié)果。GA算法操作過(guò)程中依賴(lài)于隨機(jī)數(shù)發(fā)生器，每次運(yùn)行的解是不同的，多次運(yùn)行可以獲得廣泛可能的解，對(duì)GA算法進(jìn)行了5次運(yùn)行。圖中上半部分為建模集樣品可見(jiàn)/近紅外光譜，圖中下半部分橫線(xiàn)為5次運(yùn)行的波長(zhǎng)篩選結(jié)果。從圖中可以看出5次運(yùn)行篩選的波長(zhǎng)主要集中在400～700 nm、1 100～1 500 nm波段?；谶x擇的波長(zhǎng)建立GA-PLS土壤堿解氮定量預(yù)測(cè)模型，5次預(yù)測(cè)結(jié)果的均值如表3所示。

圖4表示使用SPA算法進(jìn)行波長(zhǎng)選擇的結(jié)果。從圖中可以看出選取了4個(gè)變量，對(duì)應(yīng)波長(zhǎng)為535.85，619.66，1 091.49，1 817.58 nm?；谶@4個(gè)變量建立了SPA-PLS土壤堿解氮定量預(yù)測(cè)模型，預(yù)測(cè)結(jié)果如表3所示。雖然SPA-PLS模型參與建模的變量數(shù)僅獲得4個(gè)變量，占原始變量的 1%，SPA-PLS模型的性能適當(dāng)，R2為0.726，RPD值為1.906，結(jié)果表明PSA-PLS模型對(duì)土壤堿解氮有一定的預(yù)測(cè)能力。

圖5表示運(yùn)行一次CARS算法，隨著采樣數(shù)的增加，采樣變量數(shù)(圖5(a))、十折交叉驗(yàn)證RMSECV值(圖5(b))、各變量回歸系數(shù)(圖5(c))的變化趨勢(shì)。從圖5中可以看出，CARS變量選擇過(guò)程分兩個(gè)階段進(jìn)行。第一階段由于指數(shù)衰減函數(shù)的作用，采樣的變量數(shù)減少速度很快，稱(chēng)為粗選;第二階段采樣的變量數(shù)減少速度減慢，稱(chēng)為細(xì)選。在1～25采樣區(qū)間，RMSECV值隨著采樣數(shù)的增加減少快速，表明消除了與堿解氮無(wú)關(guān)變量;采樣超過(guò)25后，隨著采樣次數(shù)的繼續(xù)增加，RMSECV值又逐漸遞增，表明光譜中與堿解氮相關(guān)的重要信息被剔除了。圖中顯示采樣25時(shí)，RMSECV值最小(圖5(c)中星號(hào)垂線(xiàn)標(biāo)示)。因此，第25次蒙特卡羅采樣獲得的變量為預(yù)測(cè)土壤堿解氮含量的關(guān)鍵變量，共計(jì)29個(gè)特征波長(zhǎng)，依次為455.18，484.26，501.56，530.17，558.47，597.58，608.65，679.35，778.83，789.03，995.24，999.74，1 013.12，1 075.28，1 080.7，1 091.49，1 096.86，1 323.96，1 328.95，1 769.18，1 890.49，1 898.91，1 907.31，1 977.34，1 981.38，1 985.42，2 168.37，2 198.73，2 221.53 nm?；贑ARS算法獲得的29個(gè)特征波長(zhǎng)建立PLS模型對(duì)土壤堿解氮進(jìn)行定量預(yù)測(cè)，預(yù)測(cè)性能如表3所示。

圖6表示基于Random frog算法運(yùn)行10 000次結(jié)果的均值，由于Random frog算法每次運(yùn)行的結(jié)果略有差異，所以需要多次運(yùn)行并且結(jié)果取平均值。越重要的變量，被選擇的概率越大。依據(jù)經(jīng)驗(yàn)，以0.3為閾值，選擇概率大于閾值的波長(zhǎng)將被作為特征波長(zhǎng)。因此，獲得了13個(gè)超過(guò)點(diǎn)線(xiàn)(0.3)的波長(zhǎng)(501.56，789.03，912，949.56，1 283.78，1 323.96，1 719.88，1 769.18，1 773.6，1 839.26，2 095.77，2 118.83 nm)作為特征波長(zhǎng)，作為 Random frog算法運(yùn)行10 000次后的結(jié)果?；?3個(gè)特征波長(zhǎng)建立RF-PLS土壤堿解氮定量預(yù)測(cè)模型，預(yù)測(cè)結(jié)果如表3所示。

從表3可知，采用SPA提取的特征波長(zhǎng)建立的PLS模型的效果最好，預(yù)測(cè)集的決定系數(shù)為0.726，均方根預(yù)測(cè)誤差 3.616，相對(duì)分析誤差為1.906。采用CARS選擇的特征波長(zhǎng)建立的模型預(yù)測(cè)效果最差，預(yù)測(cè)集的決定系數(shù)為0.536，均方根預(yù)測(cè)誤差8.005，相對(duì)分析誤差為0.861。比較基于GA、Random frog算法提取的特征波長(zhǎng)建立的PLS模型可知，其預(yù)測(cè)能力均低于全譜PLS模型，可能是GA、Random frog選取的特征波長(zhǎng)中包含無(wú)用信息，沒(méi)有達(dá)到最優(yōu)選作用。依據(jù)模型評(píng)價(jià)標(biāo)準(zhǔn)，SPA、GA和Random frog算法及全譜PLS模型均對(duì)土壤堿解氮具備一定的預(yù)測(cè)能力，相對(duì)于全譜PLS，SPA、GA和Random frog算法簡(jiǎn)化了模型;CARS算法模型未能體現(xiàn)其對(duì)土壤堿解氮的預(yù)測(cè)能力，分析可能是由于CARS算法在波長(zhǎng)篩選過(guò)程中剔除了重要波長(zhǎng)變量。

圖5　CARS運(yùn)行結(jié)果。(a)篩選過(guò)程中采樣變量數(shù)的變化趨勢(shì);(b)10折交叉驗(yàn)證所得殘差RMSECV;(c)各變量回歸系數(shù)隨著采樣次數(shù)增加的變化趨勢(shì)。

圖6　Random frog選擇波長(zhǎng)的概率

4　結(jié)　論

利用可見(jiàn)/近紅外光譜技術(shù)對(duì)土壤堿解氮含量進(jìn)行了檢測(cè)。用SG、MSC、SNV對(duì)原始光譜數(shù)據(jù)進(jìn)行了預(yù)處理，基于全光譜建立了土壤堿解氮含量PLS預(yù)測(cè)模型，比較了不同預(yù)處理方法和原始光譜的預(yù)測(cè)效果。結(jié)果表明原始光譜數(shù)據(jù)建立的模型最優(yōu)，分析是由于光譜采集時(shí)儀器性能穩(wěn)定，樣品的顆粒度比較小和均勻(過(guò)60目篩)。在原始光譜數(shù)據(jù)的基礎(chǔ)上，采用GA算法、SPA算法、CARS算法和Random frog算法提取特征波長(zhǎng)，并基于特征波長(zhǎng)建立PLS預(yù)測(cè)模型。結(jié)果顯示，與全譜PLS模型相比，SPA算法建立的預(yù)測(cè)模型效果最好，GA算法、CARS算法和Random frog算法建立的模型的預(yù)測(cè)效果均低于全譜PLS模型，說(shuō)明SPA算法對(duì)于土壤堿解氮是一種有效的波長(zhǎng)篩選方法。本文的研究為土壤堿解氮含量的快速檢測(cè)和儀器開(kāi)發(fā)奠定了基礎(chǔ)。