鐘翔君, 楊 麗*, 張東興, 崔 濤, 和賢桃, 杜兆輝

1. 中國(guó)農(nóng)業(yè)大學(xué)工學(xué)院， 北京 100083 2. 農(nóng)業(yè)農(nóng)村部土壤-機(jī)器-植物系統(tǒng)技術(shù)重點(diǎn)實(shí)驗(yàn)室， 北京 100083

引 言

土壤有機(jī)質(zhì)(soil organic matter， SOM)是影響播量的土壤關(guān)鍵參數(shù)， 對(duì)作物的生長(zhǎng)發(fā)育起至關(guān)重要的作用[1-2]。 根據(jù)田間SOM信息對(duì)播量進(jìn)行實(shí)時(shí)調(diào)控， 可以充分挖掘土壤潛力、 節(jié)約良種用量， 對(duì)作物的提質(zhì)增效具有重要意義[3-6]。 傳統(tǒng)SOM信息的獲取多以實(shí)驗(yàn)室化學(xué)分析為主[7]， 雖然應(yīng)用較廣， 但分析過(guò)程繁瑣、 時(shí)效性差、 成本高且采樣的密度難以滿足大面積檢測(cè)需求。 近年來(lái)可見(jiàn)-近紅外光譜分析因其具有操作方便、 采樣速率快等優(yōu)勢(shì)， 還可提供高分辨率和豐富的土壤光譜信息， 成為SOM快速獲取的熱門(mén)途徑。

國(guó)內(nèi)外許多學(xué)者對(duì)SOM含量的光譜預(yù)測(cè)已開(kāi)展了大量研究[8-12]， 其中， 光譜特征篩選方法[13]有效解決光譜信息量大、 數(shù)據(jù)冗雜等造成預(yù)測(cè)模型效率低的問(wèn)題， 是光譜分析過(guò)程的重要環(huán)節(jié)。 Vohland等[14]對(duì)德國(guó)不同類型的土樣進(jìn)行光譜分析， 通過(guò)競(jìng)爭(zhēng)性自適應(yīng)重加權(quán)算法(competitive adaptive reweighted sampling， CARS)結(jié)合偏最小二乘回歸(partial least squares regression， PLSR)建立了SOM含量預(yù)測(cè)模型。 Viscarra Rossel等[15]對(duì)澳大利亞不同類型土壤有機(jī)碳的組成進(jìn)行研究， 基于決策樹(shù)算法推導(dǎo)傳遞函數(shù)并構(gòu)建預(yù)測(cè)模型來(lái)預(yù)測(cè)土壤總有機(jī)碳組分。 Shi等[16]對(duì)不同省份的土類數(shù)據(jù)進(jìn)行可見(jiàn)-近紅外光譜分析， 通過(guò)空間約束局部-偏最小二乘方法建立了SOM預(yù)測(cè)模型。 張智濤等[17]基于分?jǐn)?shù)階微分結(jié)合支持向量機(jī)分類-隨機(jī)森林構(gòu)建荒漠土SOM含量預(yù)測(cè)模型。 張娟娟等[18]分析了5種砂姜黑土樣本的光譜特征， 通過(guò)遺傳算法篩選特征波長(zhǎng)并結(jié)合支持向量機(jī)建立了預(yù)測(cè)模型。 于雷等[19]采用不同變量篩選方法對(duì)漢江平原土樣進(jìn)行特征提取， 并構(gòu)建了SOM含量預(yù)測(cè)模型。 Hong等[20]通過(guò)分?jǐn)?shù)階微分結(jié)合不同的變量篩選方法， 分析了華中地區(qū)土樣的光譜特征并構(gòu)建了SOM含量預(yù)測(cè)模型。 綜上可以看出， 利用特征變量篩選方法可以有效優(yōu)化模型， 但是不同類型土壤差異較大， 構(gòu)建的模型大多僅針對(duì)某種特定類型的土壤， 對(duì)不同土壤類型的估測(cè)精度和適用性難以估測(cè)[21]。

華北平原是全國(guó)重要的糧食和經(jīng)濟(jì)作物區(qū)， 同時(shí)是我國(guó)玉米主產(chǎn)區(qū)之一， 通過(guò)研究該區(qū)域SOM信息指導(dǎo)播種、 施肥及其他土壤改良作業(yè)， 可有效降低生產(chǎn)投入、 提高肥料利用率。 基于此， 以該區(qū)域北部的砂壤潮土為研究對(duì)象， 以高靈敏度微型可見(jiàn)-近紅外光譜儀采集并分析300～2 500 nm波長(zhǎng)范圍的光譜反射率， 以多種波長(zhǎng)選擇方法篩選出特征波長(zhǎng)， 在對(duì)不同特征波長(zhǎng)進(jìn)行建模分析的基礎(chǔ)上， 找出反演SOM的優(yōu)選方法， 為該區(qū)域SOM的快速獲取設(shè)備的設(shè)計(jì)方法和模型選擇提供參考。

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

1.1 土樣采集與處理

研究區(qū)位于河北省廊坊市(39°19′N， 116°17′E)中部平原地帶， 地處華北平原北部， 是我國(guó)玉米生產(chǎn)主區(qū)之一。 地勢(shì)平坦， 土壤類型以砂壤質(zhì)為主， 占土壤總面積90%以上， 光照充足， 溫差較大， 這些獨(dú)特的土壤及氣候條件， 使得該地區(qū)以種植玉米、 花生、 甘薯等作物為主。

在常年耕作的地塊上以五點(diǎn)采樣法采集0～20 cm耕作層的土壤樣本， 采集時(shí)去除地表殘茬及礫石， 并將采集的土樣密封帶回實(shí)驗(yàn)室進(jìn)行處理。 共采集了60份土樣， 每份大約3 kg。 為不破壞其內(nèi)部成分， 將取回的土樣分別置于恒溫干燥箱(DHG-9123A型， 上海)并在40 ℃下烘干24 h至恒重， 然后將烘干后的土壤研磨并過(guò)1 mm篩網(wǎng)后備用， 分別供實(shí)驗(yàn)室分析及光譜測(cè)試用。

1.2 土樣實(shí)驗(yàn)室檢測(cè)

樣品的SOM含量采用TOC元素分析儀(Elementar vario TOC cube， 德國(guó))進(jìn)行測(cè)定。 首先分別用萬(wàn)分之一電子天平(FA324型， 上海)稱取研磨后的土樣15～20 mg， 并置于準(zhǔn)備好的直徑4 mm、 高6 mm開(kāi)口銀囊中， 隨后在每個(gè)銀囊滴入1 mol·L-1HCl將土樣完全浸潤(rùn)， 靜置30 min后轉(zhuǎn)移至恒溫干燥箱中干燥至恒重。 將烘干后的銀囊封口并用錫紙包裹、 壓實(shí)， 隨后依次投放于TOC元素分析儀中測(cè)量其SOM含量。 為保證數(shù)據(jù)的有效性， 每個(gè)樣本準(zhǔn)備5個(gè)重復(fù)并求均值， 得到SOM含量統(tǒng)計(jì)結(jié)果如表1所示。

表1 SOM含量統(tǒng)計(jì)

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

土壤樣品的光譜數(shù)據(jù)用美國(guó)海洋光學(xué)公司的QE Pro高性能光譜儀及NIR Quest系列近紅外光譜儀同步采集。 其中， NIR Quest512-2.5近紅外光譜儀采用穩(wěn)定性高的濱松銦鎵砷化物(InGaAs)陣列探測(cè)器， 可測(cè)量900～2 500 nm波長(zhǎng)范圍的光譜數(shù)據(jù)， 光學(xué)分辨率為9.0 nm。 QE Pro高性能光纖光譜儀采用低噪音的電子部分與18位A/D轉(zhuǎn)換器， 同時(shí)配備高容量的板存緩沖區(qū)， 具有高靈敏度與寬動(dòng)態(tài)范圍特性， 可大大提高光譜檢測(cè)的準(zhǔn)確度， 同時(shí)具有很高的信噪比(大于1 000∶1)和穩(wěn)定性， 可測(cè)量185～1 100 nm可見(jiàn)-近紅外波長(zhǎng)范圍的光譜數(shù)據(jù)， 滿足高速及寬濃度范圍的快速高精度的光譜測(cè)量， 光學(xué)分辨率為1.7 nm。

圖1為光譜采集裝置實(shí)物圖， 其中， 光源為5W HL-2000-FHSA型鹵鎢燈光源(Ocean Optics, Inc., 美國(guó))， 其內(nèi)部集成風(fēng)扇冷卻、 快門(mén)和手動(dòng)衰減器功能， 可以保證持續(xù)穩(wěn)定的光源輸出； 光源配合實(shí)驗(yàn)試級(jí)QR200-12-MIXED型全光譜一分三光纖(Ocean Optics, Inc., 美國(guó))進(jìn)行試驗(yàn)， 該光纖主要包括1個(gè)入射光纖、 2個(gè)反射光纖(UV-Vis和Vis-NIR)和光纖探頭組成； 光纖探頭固定在Stage-RTL-T型多功能檢測(cè)臺(tái)(Ocean Optics, Inc., 美國(guó))光具座上， 裝有土樣的培養(yǎng)皿置于檢測(cè)臺(tái)下方的樣品支座上； 通過(guò)筆記本電腦的Ocean View軟件采集樣本的反射光譜。

圖1 光譜采集裝置圖

為降低環(huán)境及儀器噪聲的影響， 獲取高精度的光譜反射率數(shù)據(jù)， 樣品測(cè)量前用美國(guó)海洋光學(xué)公司99%漫反射標(biāo)準(zhǔn)白板進(jìn)行校正， 分別獲取開(kāi)啟光源及關(guān)閉光源后得到的亮、 暗光譜數(shù)據(jù)后， 根據(jù)式(1)運(yùn)算得到校正后的反射率數(shù)據(jù)。

(1)

式(1)中：WS為開(kāi)啟光源得到的校正亮光譜，DS為關(guān)閉光源的得到的校正暗光譜，RS為樣品初始反射率光譜，Rf為校正后的樣品反射率光譜。

經(jīng)白板校準(zhǔn)后， 將不同SOM含量的土壤樣本置于直徑3.5 mm的培養(yǎng)皿中， 通過(guò)調(diào)節(jié)檢測(cè)臺(tái)滑軌使光纖探頭位于樣品上表面， 試驗(yàn)時(shí)每采集5個(gè)樣本， 用標(biāo)準(zhǔn)白板校正1次。 其中， 試驗(yàn)時(shí)光纖探頭距標(biāo)準(zhǔn)白板及樣品的上表面高度均為3 mm。 采用五點(diǎn)法選取樣本5個(gè)位置采集光譜， 每個(gè)位置連續(xù)采集5次的均值作為該位置的反射光譜， 每個(gè)土壤樣本準(zhǔn)備3個(gè)重復(fù)， 試驗(yàn)共得到900條光譜數(shù)據(jù)。

1.4 光譜數(shù)據(jù)處理

由于低于380 nm和高于2 400 nm波長(zhǎng)的數(shù)據(jù)噪聲較大， 因此將上述波段從每組光譜數(shù)據(jù)中去除， 只保留380～2 400 nm范圍的光譜數(shù)據(jù)用于后續(xù)分析。 為降低因儀器噪聲、 測(cè)量環(huán)境及土樣表面粗糙度等因素對(duì)采樣的影響， 采用蒙特卡洛交叉驗(yàn)證法(Monte Carlo cross validation， MCCV)篩選異常數(shù)據(jù)并剔除。 對(duì)剔除異常樣本后的光譜數(shù)據(jù)采用Savitzky-Golay(SG)平滑法進(jìn)行預(yù)處理， 并用作后續(xù)分析。

1.5 SOM含量特征篩選方法

1.5.1 CARS算法

CARS方法首先抽取部分樣本作為校正集， 利用MCCV方法及PLSR構(gòu)建模型， 以模型中回歸系數(shù)絕對(duì)值權(quán)重作為基準(zhǔn)， 保留模型中權(quán)重值大的特征波長(zhǎng)并建立新的模型， 經(jīng)過(guò)多次計(jì)算， 結(jié)合交叉驗(yàn)證確定交叉驗(yàn)證均方根誤差(root mean square error of cross validation， RMSECV)小的波長(zhǎng)集合為最優(yōu)特征組合[19]。 該方法可以降低冗余數(shù)據(jù)的干擾， 從而選出優(yōu)化后的變量組合， 提高模型的穩(wěn)定性及預(yù)測(cè)效果。

1.5.2 連續(xù)投影算法

連續(xù)投影算法(successive projections algorithm， SPA)首先將校正集波長(zhǎng)矩陣投影到其他波長(zhǎng)上， 計(jì)算出每個(gè)波長(zhǎng)點(diǎn)對(duì)應(yīng)的投影值， 以投影值為基準(zhǔn)， 篩選并保留最大投影值所在的波長(zhǎng)， 通過(guò)不斷計(jì)算篩選出最優(yōu)的波長(zhǎng)組合。 通過(guò)SPA方法選擇的是冗余信息低及共線性少的變量組合， 可以在一定程度上避免光譜信息重疊， 有利于簡(jiǎn)化模型結(jié)構(gòu)、 提高運(yùn)算效率。

1.5.3 其他特征提取算法

無(wú)信息變量消除(uninformative variables elimination， UVE)方法通過(guò)噪聲信息加入到光譜數(shù)據(jù)中， 通過(guò)交叉驗(yàn)證剔除無(wú)效信息變量并建立PLSR模型， 通過(guò)對(duì)比系數(shù)矩陣的絕對(duì)值大小， 確定出特征變量組合。 變量組合集群分析法(variable combination population analysis， VCPA)采用二進(jìn)制矩陣采樣策略， 利用指數(shù)衰減函數(shù)篩選無(wú)效變量， 并依據(jù)交叉驗(yàn)證均方根誤差最終選擇出特征變量組合。

1.6 模型構(gòu)建及檢驗(yàn)

利用光譜-理化值共生距離法(sample set partitioning based on joint x-y distance， SPXY)將樣本集按7∶3劃分為建模集和預(yù)測(cè)集。 分別以全波長(zhǎng)及CARS， SPA， UVE， VCPA及CARS-SPA等不同方法篩選的特征波長(zhǎng)為自變量， SOM含量為因變量， 基于PLSR結(jié)合交叉驗(yàn)證構(gòu)建SOM含量預(yù)測(cè)模型。 分別以決定系數(shù)(R2)、 校正均方根誤差(root mean square error of calibration， RMSEC)、 預(yù)測(cè)均方根誤差(root mean square error of prediction， RMSEP)及剩余預(yù)測(cè)偏差(residual prediction deviation， RPD)等作為模型的評(píng)價(jià)指標(biāo)[16]。 其中， RPD越大、R2越接近1、 RMSEC與RMSEP越小表明模型效果越好。

2 結(jié)果與討論

2.1 預(yù)處理結(jié)果分析

采用MCCV方法分別對(duì)不同樣本的反射率數(shù)據(jù)進(jìn)行異常篩選， 其中每個(gè)樣本的光譜數(shù)據(jù)作為一個(gè)獨(dú)立的數(shù)據(jù)點(diǎn)， 分別以樣本的標(biāo)準(zhǔn)偏差作為y軸， 平均預(yù)測(cè)誤差為x軸， 對(duì)所有樣本光譜數(shù)據(jù)(數(shù)據(jù)點(diǎn))進(jìn)行篩選， 不同樣本的數(shù)據(jù)集分布結(jié)果如圖2所示。 從圖中可以看出， 不同土樣光譜的數(shù)據(jù)集離散程度不一樣， 但大部分?jǐn)?shù)據(jù)點(diǎn)在某范圍內(nèi)呈現(xiàn)集中分布。 將遠(yuǎn)離大部分?jǐn)?shù)據(jù)集分布的數(shù)據(jù)點(diǎn)(即平均誤差和標(biāo)準(zhǔn)偏差越大)視為異常樣本并予以剔除， 留下的樣本數(shù)據(jù)作為有效數(shù)據(jù)， 用于后續(xù)分析與運(yùn)算。 經(jīng)過(guò)異常值的篩選剔除， 最終共保留了809個(gè)有效數(shù)據(jù)。

圖2 MCCV異常值篩選結(jié)果

對(duì)剔除異常數(shù)據(jù)后的光譜進(jìn)行SG平滑， 得到平滑后的光譜曲線如圖3所示。 從圖中可以看出， 不同SOM含量的光譜反射率曲線總體變化趨勢(shì)類似， 隨著波長(zhǎng)的增加， 光譜反射率呈現(xiàn)先增加后減小的趨勢(shì)。 同時(shí)， 所有光譜曲線均在1 410， 1 910和2 200 nm附近出現(xiàn)明顯的水分吸收谷， 這與Laamrani等[13]得到的光譜曲線特征結(jié)論類似。 另外， 由于兩臺(tái)光譜儀在Ocean View軟件中進(jìn)行拼接， 所以在970 nm附近的反射率出現(xiàn)明顯波動(dòng)。

圖3 光譜反射率曲線

2.2 SOM含量特征變量篩選

經(jīng)過(guò)CARS， SPA， CARS-SPA， UVE及VCPA方法篩選變量結(jié)果如圖4所示， 從圖中可以看出， 不同篩選方法篩選出的波長(zhǎng)數(shù)目及波長(zhǎng)所在位置存在顯著差異。 從圖4(a)中可以看出， CARS算法在采樣次數(shù)增加至200次的過(guò)程中， 特征變量的個(gè)數(shù)逐漸減少， 其趨勢(shì)由快速下降逐漸變?yōu)槠骄彛?而RMSECV的值呈現(xiàn)先減小后增加的趨勢(shì)， 這與Hong等[20]對(duì)漢江平原土樣進(jìn)行光譜數(shù)據(jù)處理得到的結(jié)論類似。 如圖4(a)中黑色豎直線標(biāo)注， 當(dāng)采樣次數(shù)為46次時(shí)RMSECV取得最小值， 該采樣次數(shù)對(duì)應(yīng)篩選出的特征波長(zhǎng)個(gè)數(shù)為288個(gè)， 使得波段數(shù)目壓縮至全波段數(shù)目的23.4%， 波長(zhǎng)的分布如圖4(b)所示。 將基于SPXY方法劃分好的建模集和預(yù)測(cè)集數(shù)據(jù)通過(guò)SPA算法進(jìn)行計(jì)算， 結(jié)合圖4(c)可以看出， 隨著變量個(gè)數(shù)的增加， RMSECV的值大致呈現(xiàn)快速減小然后趨于穩(wěn)定的趨勢(shì)， 而當(dāng)變量個(gè)數(shù)為138個(gè)時(shí)， 其值達(dá)到最小， 篩選出的特征變量分布如圖4(d)所示， 波段數(shù)目壓縮至全波段的11.2%。 相較于CARS方法， SPA法篩選的變量共線性達(dá)到最小， 極大地減少了建模所需的波長(zhǎng)個(gè)數(shù)， 而經(jīng)過(guò)CARS方法篩選的變量個(gè)數(shù)雖然相較于全波長(zhǎng)有所降低， 但是波長(zhǎng)數(shù)量仍然較多， 在全波長(zhǎng)范圍內(nèi)均有分布， 所以采用SPA算法對(duì)CARS篩選后的變量進(jìn)行二次篩選， 進(jìn)一步優(yōu)化變量的結(jié)構(gòu)， 結(jié)果如圖4(e)和(f)所示， 共篩選出了185組特征波長(zhǎng)， 波段數(shù)目壓縮至全波段的15.0%。 通過(guò)比較UVE方法運(yùn)算得到的系數(shù)矩陣， 篩選出248組特征波長(zhǎng)， 波段數(shù)目壓縮至全波段的20.1%， 如圖4(g)所示， 該方法篩選出的波段較為集中。 經(jīng)過(guò)對(duì)比RMSECV的值， 基于VCPA方法最終篩選出100組特征波長(zhǎng)， 波段數(shù)目壓縮至全波段的8.1%， 波長(zhǎng)分布如圖4(h)所示。

2.3 模型建立與檢驗(yàn)

分別基于全波長(zhǎng)及不同變量篩選方法得到的特征波長(zhǎng)為自變量， SOM含量為因變量， 采用SPXY法將光譜數(shù)據(jù)按7∶3分為建模集和預(yù)測(cè)集， 結(jié)合留一法交叉驗(yàn)證， 構(gòu)建PLSR預(yù)測(cè)模型， 得到不同模型的預(yù)測(cè)效果如圖5所示。

圖5 不同波長(zhǎng)PLSR建模結(jié)果

利用光譜可以實(shí)現(xiàn)SOM的預(yù)測(cè)， 但是光譜波段多、 數(shù)據(jù)信息冗雜， 且土壤光譜反射率易受土壤質(zhì)地、 顏色及外部工作環(huán)境等多種因素的影響， 均為SOM的快速預(yù)測(cè)及儀器設(shè)計(jì)增加了難度。 本研究針對(duì)玉米主產(chǎn)區(qū)之一華北平原地帶的砂壤潮土進(jìn)行一致的處理以后， 對(duì)比不同的波長(zhǎng)篩選方法提取有效變量， 降低了無(wú)效信息對(duì)預(yù)測(cè)效果的干擾， 實(shí)現(xiàn)SOM含量預(yù)測(cè)。 在后續(xù)研究中， 需要考慮其他影響因素如光照、 溫度、 土壤類型等對(duì)預(yù)測(cè)效果的影響， 優(yōu)化數(shù)據(jù)處理及建模方法， 以進(jìn)一步提高SOM的預(yù)測(cè)精度， 實(shí)現(xiàn)田間SOM快速高精度檢測(cè)。

3 結(jié) 論

以玉米主產(chǎn)區(qū)之一華北平原為研究區(qū)域， 對(duì)該區(qū)域砂壤潮土進(jìn)行可見(jiàn)-近紅外光譜采集， 通過(guò)不同的波長(zhǎng)篩選方法提取有效變量并進(jìn)行SOM含量預(yù)測(cè)， 得到主要結(jié)論如下：

(1)不同方法篩選的波長(zhǎng)數(shù)目及波長(zhǎng)位置存在顯著差異， CARS和SPA算法選擇的光譜特征在整個(gè)光譜范圍都有分布， UVE和VCPA篩選的波段較為集中， 且基于CARS-SPA方法可以進(jìn)一步優(yōu)選特征變量， 其特征波長(zhǎng)僅為全波長(zhǎng)數(shù)量的15%。

(2)通過(guò)對(duì)比不同模型的建模及預(yù)測(cè)效果， 除UVE和VCPA算法外， 其余算法構(gòu)建的模型均能實(shí)現(xiàn)SOM含量的有效預(yù)測(cè)， 其RPD值均大于2.0。