基于SGA-RF算法的農(nóng)業(yè)土壤鎘濃度反演研究

王軒慧 陳建毅 鄭西來(lái) 朱 成 王軒力 單春芝

(1.中國(guó)海洋大學(xué)海洋環(huán)境與生態(tài)教育部重點(diǎn)實(shí)驗(yàn)室， 青島 266100; 2.青島農(nóng)業(yè)大學(xué)理學(xué)與信息科學(xué)學(xué)院， 青島 266109; 3.中國(guó)聯(lián)合網(wǎng)絡(luò)通信有限公司濟(jì)南軟件研究院項(xiàng)目管理部， 濟(jì)南 250000; 4.山西工程技術(shù)學(xué)院信息工程與自動(dòng)化系， 陽(yáng)泉 045000; 5.國(guó)家海洋局北海環(huán)境監(jiān)測(cè)中心， 青島 266033)

0 引言

農(nóng)業(yè)土壤重金屬污染已成為我國(guó)主要的環(huán)境問(wèn)題之一[1-2]。土壤中的鎘會(huì)在農(nóng)作物中富集并進(jìn)入食物鏈，從而對(duì)人類健康造成嚴(yán)重威脅[3-4]。因此，密切監(jiān)測(cè)農(nóng)業(yè)土壤中的重金屬鎘超標(biāo)情況對(duì)于防止農(nóng)產(chǎn)品污染，保障人體健康有著重要意義。

目前，常規(guī)化學(xué)方法對(duì)土壤重金屬濃度的測(cè)定不僅周期長(zhǎng)，成本高，而且會(huì)產(chǎn)生大量危害環(huán)境的廢棄物[5]。相比之下，可見(jiàn)/近紅外光譜技術(shù)具有成本低、效率高、環(huán)境污染少等優(yōu)點(diǎn)[6-7]，逐漸成為有效監(jiān)測(cè)土壤重金屬污染情況的重要手段[8-9]。CHEN等[10]采用正交信號(hào)校正預(yù)處理方法和反向傳播神經(jīng)網(wǎng)絡(luò)建立污水灌溉區(qū)的重金屬鎘反演模型。GHOLIZADEH等[11]結(jié)合各種預(yù)處理方法與支持向量機(jī)回歸建立土壤重金屬含量與反射光譜之間的對(duì)應(yīng)關(guān)系。TAYEBI等[12]使用可見(jiàn)/近紅外反射光譜和偏最小二乘(PLSR)模型來(lái)預(yù)測(cè)農(nóng)業(yè)土壤中重金屬鐵的濃度，指出二階微分預(yù)處理與PLSR結(jié)合建立的模型反演精度較高。

盡管上述方法成功地將可見(jiàn)/近紅外光譜技術(shù)應(yīng)用于土壤重金屬高光譜反演領(lǐng)域，但由于原始光譜數(shù)據(jù)中存在大量冗余和不相關(guān)特征，從而嚴(yán)重影響了反演模型的性能。特別是農(nóng)業(yè)土壤中重金屬鎘所對(duì)應(yīng)的可見(jiàn)/近紅外光譜波段非常少，因而直接使用原始光譜構(gòu)建回歸反演模型不僅準(zhǔn)確率和穩(wěn)定性偏低，而且運(yùn)算效率較差。因此，只有通過(guò)使用特征選擇方法才能選出具有最小冗余度且最具代表性的敏感波段，進(jìn)一步有效地提高農(nóng)業(yè)土壤鎘回歸反演模型的預(yù)測(cè)能力[13]。JIANG等[14]利用遺傳算法結(jié)合PLSR建立土壤鎘含量反演模型。雖然上述方法取得了一定成果，但是該方法對(duì)相鄰光譜波段的特征值進(jìn)行了融合，這會(huì)導(dǎo)致部分潛在重要特征的丟失。XIA等[15]采用相關(guān)系數(shù)法結(jié)合PLSR成功反演長(zhǎng)江流域沉積物中的鎘含量。但基于相關(guān)系數(shù)法的特征選擇方法只適用于線性統(tǒng)計(jì)方法建立的模型，對(duì)于非線性相關(guān)及樣本分布不均勻的情況，反演精度較低。

本文在重金屬鎘高光譜反演領(lǐng)域，針對(duì)上述特征選擇技術(shù)的不足，提出一種基于斯皮爾曼等級(jí)相關(guān)分析的遺傳隨機(jī)森林(Spearman’s rank correlation coefficient-based genetic algorithm using random forest，SGA-RF)特征選擇算法?；诖蠊梁恿饔虮韺愚r(nóng)業(yè)土壤可見(jiàn)/近紅外高光譜數(shù)據(jù)，通過(guò)與傳統(tǒng)遺傳算法和基于K近鄰的遺傳算法進(jìn)行比較，探討SGA-RF算法在農(nóng)業(yè)土壤重金屬高光譜反演領(lǐng)域應(yīng)用的可能性，以期為利用高光譜技術(shù)檢測(cè)農(nóng)業(yè)土壤中低含量重金屬污染狀況提供一定的理論依據(jù)。

1 SGA-RF特征選擇算法

圖1 SGA-RF特征選擇流程圖Fig.1 Flow chart of SGA-RF feature selection

SGA-RF特征選擇算法流程如圖1所示。該算法由2個(gè)子算法組成，是一種基于Wrapper[16]型的特征選擇方法。首先，為了消除冗余和無(wú)關(guān)信息，采用斯皮爾曼相關(guān)分析特征選擇方法(算法1)對(duì)2 051個(gè)預(yù)處理后的初始特征集合進(jìn)行預(yù)選，得到了彼此無(wú)關(guān)聯(lián)的特征波段子集作為算法2的輸入。其次，在遺傳隨機(jī)森林算法(算法2)中采用基于隨機(jī)森林袋外誤差的適應(yīng)度函數(shù)來(lái)評(píng)價(jià)個(gè)體的優(yōu)劣，最終選出適應(yīng)度函數(shù)最小值所對(duì)應(yīng)的特征子集為最優(yōu)特征波段子集。最后，利用最優(yōu)特征波段子集構(gòu)建隨機(jī)森林回歸(Random forest regression, RFR)反演模型。

以斯皮爾曼等級(jí)相關(guān)分析作為SGA-RF算法的特征預(yù)選方法，主要目的是有效減少初始特征集維數(shù)并剔除原始光譜特征波段之間的冗余和無(wú)關(guān)信息。該算法包括：①數(shù)據(jù)準(zhǔn)備階段(步驟1～步驟2)，通過(guò)計(jì)算兩兩特征波段之間的斯皮爾曼等級(jí)相關(guān)系數(shù)將強(qiáng)相關(guān)的特征波段存入二維矩陣X中。②算法的核心組成部分(步驟3～步驟8)，通過(guò)遍歷矩陣X構(gòu)造所有相關(guān)特征的分組Y。③算法輸出部分(步驟9～步驟11)，通過(guò)對(duì)特征波段集合Q與Y中每組代表波段求并集得到最優(yōu)光譜特征波段子集。算法的具體描述如算法1所示。

(1)算法1——基于斯皮爾曼等級(jí)相關(guān)分析的特征選擇算法

輸入：①預(yù)處理后的光譜特征集合M。②循環(huán)變量i1。③循環(huán)變量j1。④循環(huán)變量b1。⑤循環(huán)變量c1。⑥最大波長(zhǎng)數(shù)L2 051。

輸出：所有相關(guān)光譜特征分組集合Y(所有相關(guān)的分組)；最優(yōu)光譜特征子集F。

步驟1: 計(jì)算Mi與Mj兩個(gè)特征波段之間的斯皮爾曼等級(jí)相關(guān)系數(shù)，即

(1)

式中N——樣本數(shù)量

d——兩列成對(duì)變量的等級(jí)差分集合

在式(1)中,對(duì)Mi與Mj進(jìn)行排序(同時(shí)為升序或降序)，得到兩個(gè)元素等級(jí)集合m和n，其中mi為mi在Mi中的等級(jí)，ni為ni在Mj中的等級(jí)。將集合m、n中的元素對(duì)應(yīng)相減得到一個(gè)排行差分集合d，其中di=mi-ni，1≤i≤N。如果ρ的絕對(duì)值(|ρ|)大于或等于0.8且標(biāo)準(zhǔn)偏差p小于10%，則特征值Mi與Mj被認(rèn)為是顯著相關(guān)的，將Mi與Mj存入矩陣X中，X為二維矩陣。

步驟2:重復(fù)執(zhí)行步驟1，直至j達(dá)到最大波長(zhǎng)數(shù)L。

步驟3:把i賦值給j，重復(fù)執(zhí)行步驟1～步驟2，直至i達(dá)到最大波長(zhǎng)數(shù)L。

步驟4：把X中的第一行兩個(gè)元素賦值給Yb。設(shè)置ismatch為0，a為1。

步驟5:取出矩陣Xa中第a行的第1個(gè)元素Xa1。如果矩陣Yb中不含有Xa1，轉(zhuǎn)入步驟6；否則設(shè)置ismatch為1，再判斷矩陣Y中是否含有Xa2，如果矩陣Y中不含Xa2，將Xa2存入Yb中,轉(zhuǎn)入步驟7。

步驟6:取出矩陣Xa中第a行的第2個(gè)元素Xa2。如果矩陣Yb中含有Xa2，設(shè)置ismatch為1，并判斷矩陣Yb中是否含Xa1，如果矩陣Yb中不含有Xa1，將Xa1存入Yb中。

步驟7:如果ismatch為1，在X中剔除Xa，得到新的矩陣Xa與Yb。

步驟8:重復(fù)執(zhí)行步驟5～步驟7，直到遍歷完Xa的所有元素。Yb中存儲(chǔ)一組顯著相關(guān)的特征波段，Xa為去掉與Yb相關(guān)特征波長(zhǎng)之后的剩余特征波長(zhǎng)變量。若Yb中含有不相鄰波段，對(duì)該組進(jìn)行標(biāo)注。

步驟9:重復(fù)執(zhí)行步驟4～步驟8，直到遍歷完X的所有元素，輸出所有相關(guān)特征的分組Y，通過(guò)計(jì)算MY(Y在M中的余集)得到特征波段集合Q。

步驟10:計(jì)算Yc這一組內(nèi)部所有斯皮爾曼相關(guān)系數(shù)的平均值和相應(yīng)的標(biāo)準(zhǔn)偏差，保留與待測(cè)組分相關(guān)性最大的特征波段作為這一組的代表性特征波段。

步驟11:重復(fù)執(zhí)行步驟10，直至處理完Y的所有分組，每組代表波段組成集合P。

步驟12:通過(guò)計(jì)算集合Q∪P，輸出最優(yōu)光譜特征子集F。

采用算法1的特征選擇方法，不僅能夠剔除原始全光譜的所有相關(guān)特征值，而且大幅度地縮減了下一步特征精選(算法2)的篩選范圍，能夠達(dá)到提高算法反演精度和執(zhí)行效率的雙重目的。對(duì)于波長(zhǎng)特征選擇，遺傳算法是一種非常有效的方法。適應(yīng)度函數(shù)[17]的設(shè)計(jì)在遺傳特征選擇算法中起著至關(guān)重要的作用，但是常規(guī)的適應(yīng)度函數(shù)與具體的解決問(wèn)題領(lǐng)域聯(lián)系不夠緊密，導(dǎo)致最終篩選出的特征子集建立的模型預(yù)測(cè)精度較低。為此提出了一種新型的適應(yīng)度函數(shù)，綜合考慮了特征子集的回歸性能和算法的運(yùn)行效率，該算法的具體描述如算法2所示。

(2)算法2——遺傳隨機(jī)森林的特征選擇算法

輸入：①光譜特征集合F。②每一個(gè)體的染色體長(zhǎng)度(基因)m集合F中元素的數(shù)量。③染色體的基因位數(shù)p1。④初始群體包含的個(gè)體數(shù)f。⑤最大迭代次數(shù)Generations200。⑥停滯代數(shù)StallGenLimit100。⑦交叉概率pn0.8。⑧變異算子pm0.5。

輸出：最優(yōu)光譜特征子集E。

步驟1(參數(shù)編碼): 用一串含有mp個(gè)0/1字符(基因)的字符串(染色體串)來(lái)表示每種區(qū)間組合。若基因?yàn)?表示該特征被選中，基因?yàn)?表示該特征未被選中。

步驟2(種群初始化):隨機(jī)產(chǎn)生fm個(gè)字符作為初始群體。

步驟3(新的適應(yīng)度函數(shù))：本算法中的適應(yīng)度函數(shù)被定義為隨機(jī)森林袋外誤差和特征維度的加權(quán)和，即

(2)

式中α——RFR袋外誤差，是衡量RFR模型性能的重要指標(biāo)

β——特征子集含有的特征值數(shù)量，直接決定算法的運(yùn)行效率

ω——回歸精度在適應(yīng)度函數(shù)中所占的比重

m——染色體長(zhǎng)度(原始特征值數(shù)量)

步驟4(復(fù)制)：采用精英保留策略。設(shè)置精英個(gè)數(shù)c1為2，按照式(2)計(jì)算所有個(gè)體的適應(yīng)度，選出適應(yīng)度最低的前兩名個(gè)體不進(jìn)行配對(duì)交叉自動(dòng)推到下一代。

步驟5(選擇): 采用錦標(biāo)賽選擇策略。當(dāng)精英子代從當(dāng)前種群移除后，對(duì)剩余m-2個(gè)個(gè)體執(zhí)行錦標(biāo)賽選擇算法。設(shè)置參加錦標(biāo)賽的個(gè)體個(gè)數(shù)為2，依據(jù)式(2)評(píng)估所有個(gè)體的適應(yīng)度，從種群中選出2個(gè)最優(yōu)個(gè)體作為父代個(gè)體。

步驟6(交叉): 對(duì)錦標(biāo)賽策略生成的父代個(gè)體進(jìn)行交叉操作。選用的交叉函數(shù)是算術(shù)交叉類型的異或運(yùn)算，對(duì)兩個(gè)父類染色體進(jìn)行異或操作產(chǎn)生交叉子代。產(chǎn)生交叉子代的數(shù)量為c2=round((m-c1)pn)。

步驟7(變異): 變異操作采用均勻變異的方式，最終產(chǎn)生變異子代的數(shù)量c3=m-c1-c2。

步驟8(生成子代個(gè)體): 精英子代，交叉子代和變異子代的所有個(gè)體組成了新一代種群。利用式(2)的適應(yīng)度函數(shù)評(píng)價(jià)新一代種群的適應(yīng)性。

步驟9:迭代步驟4～步驟8，直到滿足停止條件。停止條件包含如下兩類：①最大迭代次數(shù)達(dá)到設(shè)定值Generations。 ②停滯代數(shù)達(dá)到設(shè)定值StallGenLimit。

步驟10:輸出最優(yōu)光譜特征子集E。

傳統(tǒng)的適應(yīng)度函數(shù)一般只考慮回歸模型中的交叉驗(yàn)證均方根誤差和相關(guān)系數(shù)，提出的適應(yīng)度函數(shù)綜合考慮特征子集的回歸性能和算法的運(yùn)行效率。由式(2)可以得出，當(dāng)隨機(jī)森林的袋外誤差取得最小值，并且選出的特征值數(shù)量最少時(shí)，本文算法設(shè)計(jì)的適應(yīng)度函數(shù)取得最小值。

在變異操作算法中，首先生成符合某一范圍內(nèi)均勻分布的隨機(jī)數(shù)集合(RD)，每個(gè)隨機(jī)數(shù)的值與染色體中每個(gè)基因座上的原有基因相關(guān)聯(lián)。然后從左到右掃描染色體，將每個(gè)RD的值與變異算子pm進(jìn)行比較。如果位置i處的RD小于pm，則位置i處的基因被翻轉(zhuǎn)；否則，該基因不被翻轉(zhuǎn)。

2 實(shí)驗(yàn)

2.1 數(shù)據(jù)來(lái)源

在青島市大沽河流域采集具有代表性的農(nóng)田表層土壤樣品124個(gè)，攤開并置于陰涼處，待自風(fēng)干后，去除土樣中的石子和動(dòng)植物殘?bào)w等明顯雜物，用研缽研磨，再過(guò)100目篩子保存于自封袋中待測(cè)。準(zhǔn)確稱取0.250 0 g樣品于150℃下以“HNO3-HF-H2O2”消解體系和密封高壓釜消解罐法進(jìn)行消解，利用電感耦合等離子體發(fā)射光譜儀(ICP-OES，Optima 8000，PerkinElmer，USA)來(lái)檢測(cè)124個(gè)樣品的鎘含量，并保留每份土樣剩余部分用于高光譜數(shù)據(jù)采集。

2.2 光譜采集及預(yù)處理

本次實(shí)驗(yàn)選取ASD FieldSpec 3型便攜式光譜儀測(cè)量土壤樣品的光譜反射率。將土壤樣品放入直徑10 cm、厚度約為2 cm的培養(yǎng)皿中。在50 W鹵素?zé)糇鳛楣庠吹陌凳抑袦y(cè)量土壤樣品的光譜，光源距樣品35 cm，天頂角為30°，光譜儀探頭安裝在樣品垂直上方15 cm處。每次實(shí)驗(yàn)開始之前和每測(cè)量10組樣品之后使用具有100%反射率的標(biāo)準(zhǔn)化白板進(jìn)行校正。每個(gè)土樣測(cè)量時(shí)轉(zhuǎn)動(dòng)培養(yǎng)皿3次，每次轉(zhuǎn)動(dòng)90°，每個(gè)方向上取樣10次，4個(gè)方向上共得到40條光譜測(cè)量值，算術(shù)平均后作為該樣品的最終光譜反射率。

使用Kennard-Stone[18]算法從原始光譜數(shù)據(jù)集(包含124個(gè)樣本)中選出100個(gè)樣本作為校驗(yàn)集，剩余的24個(gè)樣本作為預(yù)測(cè)集用于模型驗(yàn)證。剔除具有低信噪比的350～399 nm、2 451～2 500 nm兩個(gè)邊緣波段，在400～2 450 nm波段上對(duì)原始光譜數(shù)據(jù)采用多種預(yù)處理方法。基于隨機(jī)森林回歸模型，不同預(yù)處理方法所建模型的預(yù)測(cè)性能如表1所示。

由表1可得，一階微分處理方法所建立的RFR模型結(jié)果最優(yōu)，此時(shí)預(yù)測(cè)集相關(guān)系數(shù)、相對(duì)分析誤差最大，分別為0.897 3和1.26，均方根誤差達(dá)到最小值0.079 4?？梢?jiàn)，一階微分預(yù)處理最適于獲取含量較少組分的光譜信息，與文獻(xiàn)[19]研究結(jié)果一致。因此，后續(xù)特征選擇與建模均在一階微分分析基礎(chǔ)上進(jìn)行。

大沽河流域124個(gè)土壤樣本的原始可見(jiàn)/近紅外光譜曲線如圖2a所示。從圖2a可以看出，所有樣本的反射率曲線特征基本相似。與文獻(xiàn)[20]報(bào)道的一致，大多數(shù)土壤光譜在波長(zhǎng)1 400 nm、1 800～2 000 nm和2 200～2 400 nm附近都會(huì)出現(xiàn)明顯的水分吸收峰。圖2b為最優(yōu)預(yù)處理方式(一階微分)后的光譜曲線。從圖2b可以看出，一階微分預(yù)處理不僅能夠加強(qiáng)原始光譜的3個(gè)強(qiáng)吸收峰，而且能夠明顯增強(qiáng)原始光譜中1 000 nm和高于2 200 nm的弱吸收峰。這表明一階微分預(yù)處理方法能夠增強(qiáng)樣品之間的光譜特征差異，較適應(yīng)于農(nóng)業(yè)土壤鎘元素的高光譜響應(yīng)。

表1 不同光譜預(yù)處理方法的隨機(jī)森林建模結(jié)果Tab.1 Prediction results of RFR modeling by using different pre-processing methods

注：D1表示對(duì)光譜矩陣求一階導(dǎo)，D2表示對(duì)光譜矩陣求二階導(dǎo)，SNV(Standard normal variate transformation)表示對(duì)光譜矩陣進(jìn)行標(biāo)準(zhǔn)正態(tài)變量變換，S-G(Savitzky-Golay)表示對(duì)光譜矩陣進(jìn)行Savitzky-Golay卷積平滑。

3 結(jié)果與討論

3.1 基于斯皮爾曼等級(jí)相關(guān)分析的特征選擇算法結(jié)果分析

一般來(lái)說(shuō)，特征值之間的相關(guān)性是通過(guò)相關(guān)系數(shù)來(lái)衡量的。由于斯皮爾曼等級(jí)相關(guān)系數(shù)對(duì)數(shù)據(jù)條件的要求沒(méi)有皮爾遜相關(guān)系數(shù)嚴(yán)格，并且可以使用單調(diào)函數(shù)來(lái)描述變量之間的相關(guān)性，因此SGA-RF算法選用斯皮爾曼等級(jí)相關(guān)分析作為特征預(yù)選方法剔除冗余特征。對(duì)于初始的2 051個(gè)特征波段應(yīng)用斯皮爾曼等級(jí)相關(guān)系數(shù)分析，在校驗(yàn)集上篩選的特征波長(zhǎng)變量分布如圖3所示。

圖2 124個(gè)土壤反射光譜Fig.2 Soil reflectance spectra of 124 soil samples

圖3 算法1篩選的波長(zhǎng)變量分布圖(不相鄰波段之間的相關(guān)性用紅色豎線表示)Fig.3 Distribution diagrams of spectral features selected by algorithm 1 (there was also a correlation between nonadjacent bands, which were shown on red long string)

算法1將原始全光譜變量從2 051個(gè)減少到108組特征波段,每組內(nèi)部選出與鎘濃度相關(guān)性最強(qiáng)的特征值為該組代表波段，所有代表波段在400～2 450 nm的分布如圖3a所示。為了更加清晰地觀察不相鄰波段之間是否存在相關(guān)性，將圖3a中藍(lán)色框內(nèi)的局部波長(zhǎng)變量放大于圖3b中，將圖3a中紫色框內(nèi)的局部波長(zhǎng)變量放大于圖3c中。其中，對(duì)特征值的標(biāo)注已在算法1中步驟8完成。在圖3c中以2 446 nm為代表波段的這一組被標(biāo)識(shí)為黑色，表明該組內(nèi)部所有成員均為相鄰波段。在圖3a中以999 nm為代表波段的這一組被標(biāo)識(shí)為紅色，表明該組內(nèi)部包含不相鄰波段。以特征波段999 nm與1 829 nm的相關(guān)性為例，算法1判定二者具有相關(guān)性；但是采用傳統(tǒng)遺傳算法的區(qū)間分割法，二者不具有相關(guān)性，原因是999 nm和1 829 nm肯定不會(huì)被劃分到同一個(gè)子區(qū)間內(nèi)部。因此，傳統(tǒng)遺傳算法按照固定區(qū)間大小均分原始光譜，會(huì)導(dǎo)致部分潛在重要特征值的丟失。

從圖3a可以看出，預(yù)選出的108個(gè)特征波段主要位于1 550～1 850 nm和2 300～2 450 nm的范圍內(nèi)。分析圖3b 可得，1 550～1 720 nm范圍內(nèi)共有9組特征波段內(nèi)部含有不相鄰波段。由圖3c可得，1 800～2 450 nm范圍內(nèi)共有7組特征波段內(nèi)部含有不相鄰波段。108組相關(guān)特征中共有17組內(nèi)部存在不相鄰波段之間具有相關(guān)性的現(xiàn)象?？梢?jiàn)，相鄰波段之間存在高相關(guān)性的概率較大，但是不相鄰波段之間也有可能存在相關(guān)性。因此，基于斯皮爾曼等級(jí)相關(guān)分析的特征選擇方法克服了傳統(tǒng)遺傳算法只對(duì)相鄰波段特征值進(jìn)行融合的局限性，采用一種更加合理的方法縮減遺傳算法的自變量個(gè)數(shù)，使得遺傳算法能夠在冗余度最小且最具代表性的敏感波段中進(jìn)行全局搜索，最終達(dá)到優(yōu)化反演模型，提高預(yù)測(cè)能力的目的。

3.2 系數(shù)對(duì)適應(yīng)度函數(shù)的影響

適應(yīng)度函數(shù)的表達(dá)式(式(2))中包含一個(gè)權(quán)重參數(shù)ω，取值范圍為[0,1]，該參數(shù)表示袋外誤差在整個(gè)適應(yīng)度函數(shù)中所占比重。為了找到適應(yīng)度函數(shù)取最小值時(shí)所對(duì)應(yīng)的權(quán)重參數(shù)ω，需要考察不同權(quán)重對(duì)反演結(jié)果的影響。表2給出了不同權(quán)重參數(shù)ω所對(duì)應(yīng)的適應(yīng)度以及相對(duì)應(yīng)的相關(guān)系數(shù)、均方根誤差和預(yù)測(cè)相對(duì)分析誤差。

從表2可以看出，當(dāng)權(quán)重參數(shù)ω為0.3，適應(yīng)度函數(shù)取得最小值。因此，最優(yōu)權(quán)重參數(shù)ω為0.3。

3.3 特征選擇結(jié)果分析

在算法2中，遺傳算法種群規(guī)模被定義為100，運(yùn)行總代數(shù)為200，基因位數(shù)設(shè)為1。

表2 不同權(quán)重參數(shù)ω所對(duì)應(yīng)的適應(yīng)度以及相對(duì)應(yīng)的相關(guān)系數(shù)、均方根誤差和預(yù)測(cè)相對(duì)分析誤差Tab.2 Values of fitness function (fi) and corresponding R, RMSE and RPD according to several ω values

圖4 SGA-RF算法特征選擇過(guò)程適應(yīng)度函數(shù)進(jìn)化曲線Fig.4 Evolution curves of fitness function during feature selection process of SGA-RF

由圖4可得，從第109代到第200代，適應(yīng)度在0.014～0.016 5之間，變動(dòng)幅度很小。最佳個(gè)體的適應(yīng)度為0.015 397 7，平均值約為0.019 616 8。因此，可以判斷SGA-RF算法沒(méi)有出現(xiàn)過(guò)早收斂現(xiàn)象，原因是SGA-RF算法設(shè)置了正確有效的適應(yīng)度，同時(shí)采用精英保留策略和錦標(biāo)賽選擇算法確保了該算法的必然收斂。

采用SGA-RF特征波段選擇方法最終選出37個(gè)敏感波段，它們?cè)?00～2 450 nm波段的分布圖如圖5所示。圖中黑色曲線為某個(gè)土壤樣本的原始光譜曲線，紅色空心圓圈表示特征波段所在位置。圖5b和圖5c分別放大了圖5a中矩形框和正方形框內(nèi)的敏感波段密集區(qū)域。

圖5 基于SGA-RF算法的鎘可見(jiàn)/近紅外光譜特征波段分布圖Fig.5 Visible/NIR spectra distribution diagrams for soil Cr based on variables selected by SGA-RF (shown in red hollow circle markers)

由圖5可知，重金屬鎘的敏感波段主要位于1 600～1 800 nm和2 350～2 450 nm的范圍內(nèi)。已有研究表明，位于這2個(gè)范圍的吸收峰主要與有機(jī)化合物中的O—H、 N—H和C—H基團(tuán)有關(guān)[21]。這可能與重金屬鎘的有機(jī)化合物吸附有關(guān)，吸附機(jī)理有待進(jìn)一步深入研究。

3.4 與其他特征選擇算法的比較與分析

為了驗(yàn)證算法1，分別對(duì)不采用算法1的遺傳隨機(jī)森林算法(GA-RF)和采用算法1做預(yù)選的SGA-RF算法進(jìn)行建模性能比較，結(jié)果如表3所示。兩種算法均采用RFR建立回歸模型。

表3 評(píng)估算法1對(duì)本文算法預(yù)測(cè)效果的影響Tab.3 Impact of algorithm 1 on predictive effect by using SGA-RF

由表3可以看出，與GA-RF算法相比，SGA-RF算法不僅建模時(shí)間短，特征值數(shù)量少，而且模型性能明顯較高。因?yàn)樗惴?在有效剔除初始特征集冗余信息的同時(shí)能夠保留與鎘元素相關(guān)性最強(qiáng)的重要特征，大幅度減少特征精選階段的輸入特征值維數(shù)，不僅減少了計(jì)算復(fù)雜度，而且提高了模型反演能力。

為了進(jìn)一步驗(yàn)證SGA-RF特征選擇算法的有效性，本文算法與其他3類特征選擇方法進(jìn)行對(duì)比，結(jié)果如表4所示。這3類特征選擇方法分別為：全光譜不做特征波段選擇方法、單一特征選擇方法和以K最近鄰法構(gòu)建適應(yīng)度函數(shù)的遺傳算法 (SGA-KNN)。單一特征選擇方法包含兩種方法：原始遺傳算法與偏最小二乘算法結(jié)合的特征波段選擇方法(GA-PLSR)和斯皮爾曼等級(jí)相關(guān)分析(Spearman Rank)方法。這3類特征選擇算法均采用RFR作為回歸模型，并以選出特征波段的數(shù)量、建模時(shí)間、相關(guān)系數(shù)、均方根誤差和預(yù)測(cè)相對(duì)分析誤差作為模型的評(píng)價(jià)指標(biāo)。

依據(jù)文獻(xiàn)[21]，GA-PLSR算法采用2～20 nm不同波段間隔對(duì)原始光譜進(jìn)行重采樣，每個(gè)子區(qū)間取平均數(shù)作為新的光譜自變量，以交叉驗(yàn)證均方根誤差(RMSECV)最小值作為最優(yōu)特征譜區(qū)篩選的標(biāo)準(zhǔn),最終得到12 nm為最優(yōu)光譜區(qū)間劃分間隔。GA-PLSR算法使用PLS-Toolbox[22]重復(fù)運(yùn)行50次，以PLSR算法的RMSECV作為適應(yīng)度函數(shù)。SGA-KNN算法采用文獻(xiàn)[23]中的基于K最近鄰分類法(KNN)的適應(yīng)度函數(shù)，其中參數(shù)K選取3。上述2種算法的種群規(guī)模、運(yùn)行總代數(shù)和SGA-RF算法保持一致。

表4 本文特征選擇算法與其他特征選擇算法性能比較Tab.4 Comparison of performance between proposed and other feature selection methods

由表4可知，與其他3類特征選擇算法相比，本文算法的建模時(shí)間最短，相關(guān)系數(shù)和預(yù)測(cè)相對(duì)分析誤差最高，取得最佳預(yù)測(cè)效果。全光譜和斯皮爾曼等級(jí)相關(guān)分析建模的預(yù)測(cè)相對(duì)分析誤差均在1.4以下，說(shuō)明這2種方法建立土壤鎘含量預(yù)測(cè)模型的性能較差，原因是全光譜中含有大量的冗余和無(wú)關(guān)信息；斯皮爾曼等級(jí)相關(guān)分析只能去除特征波段之間的冗余，無(wú)法剔除那些與待測(cè)組分無(wú)關(guān)的特征值，導(dǎo)致模型精度較低。GA-PLSR算法和SGA-KNN算法建模的預(yù)測(cè)相對(duì)分析誤差均在1.0以下，說(shuō)明以GA-PLSR算法和SGA-KNN算法作為特征波段選擇方法無(wú)法建立土壤鎘含量預(yù)測(cè)模型。原因?yàn)椋孩?GA-PLSR算法不僅忽略了不相鄰波段之間的相關(guān)性，同時(shí)也沒(méi)有考慮特征波段與待測(cè)組分之間的相關(guān)性，而且該算法用于評(píng)估每個(gè)染色體的PLSR模型為線性模型。②由于KNN算法對(duì)參數(shù)K的取值，距離的度量和多重共線性都非常敏感，導(dǎo)致該算法對(duì)特征的回歸能力較差，并不適用于構(gòu)造遺傳算法的適應(yīng)度函數(shù)。③SGA-RF算法首先采用斯皮爾曼等級(jí)相關(guān)分析進(jìn)行特征預(yù)選，去掉所有特征波段之間的冗余信息；其次，結(jié)合遺傳算法的全局搜索能力和RFR的較高反演能力，優(yōu)選出適應(yīng)度最低的特征子集。其中，在適應(yīng)性函數(shù)中采用非線性模型RFR評(píng)估每個(gè)個(gè)體的性能，而RFR對(duì)于參數(shù)設(shè)置與多重共線性問(wèn)題均不敏感。

在SGA-RF算法特征篩選結(jié)果的基礎(chǔ)上，以RFR作為回歸模型，模型校驗(yàn)集和預(yù)測(cè)集樣本的估測(cè)值和實(shí)測(cè)值之間的關(guān)系如圖6所示。

圖6 SGA-RF算法結(jié)合RFR模型建模樣本、檢驗(yàn)樣本的實(shí)測(cè)值與估測(cè)值比較Fig.6 Comparison of measured Cd contents and estimated values of modeling and testing samples through SGA-RF algorithm with RFR model

由圖6可得，采用SGA-RF特征選擇算法在校驗(yàn)集和預(yù)測(cè)集均取得了較好的反演效果，驗(yàn)證了該算法與隨機(jī)森林回歸模型相結(jié)合能夠有效解決農(nóng)業(yè)土壤中痕量級(jí)重金屬含量的反演問(wèn)題。

4 結(jié)論

(1)基于大沽河流域農(nóng)業(yè)土壤樣本真實(shí)數(shù)據(jù)集，SGA-RF算法能夠?qū)? 051個(gè)原始特征波段優(yōu)選至37個(gè)。以選出的37個(gè)敏感波段作為自變量建立大沽河流域農(nóng)業(yè)土壤鎘高光譜反演模型，該模型具有較低的預(yù)測(cè)均方根誤差(0.060 1)，較高的相關(guān)系數(shù)(0.950 2)和預(yù)測(cè)相對(duì)分析誤差(2.03)。建模結(jié)果證實(shí)了該算法在降低計(jì)算復(fù)雜度的基礎(chǔ)上提高了土壤鎘含量反演模型的精度。

(2)斯皮爾曼等級(jí)相關(guān)分析將原始2 051波段縮減到108組特征子集，通過(guò)分析每組特征子集的組內(nèi)相關(guān)性，發(fā)現(xiàn)在不相鄰波段之間也存在相關(guān)性。說(shuō)明基于斯皮爾曼等級(jí)的特征預(yù)選方法能夠克服傳統(tǒng)遺傳算法的局限性，最大限度地縮減遺傳算法的自變量個(gè)數(shù)，是SGA-RF算法的重要組成部分。

(3)以相關(guān)系數(shù)、均方根誤差和預(yù)測(cè)分析相對(duì)誤差作為評(píng)價(jià)標(biāo)準(zhǔn)，本文算法分別與其他4類特征選擇方法進(jìn)行了對(duì)比。比較結(jié)果說(shuō)明，SGA-RF算法可以有效剔除冗余波長(zhǎng)變量，提取最具代表性的特征波長(zhǎng)變量，減少建模時(shí)間，提高預(yù)測(cè)效果。

(4)SGA-RF算法有效解決了土壤重金屬鎘高光譜反演領(lǐng)域建模時(shí)間長(zhǎng)，預(yù)測(cè)精度差的問(wèn)題，可用于其他農(nóng)業(yè)土壤的重金屬鎘污染監(jiān)測(cè)與識(shí)別。