程介虹，陳爭(zhēng)光,2*

1. 黑龍江八一農(nóng)墾大學(xué)電氣與信息學(xué)院，黑龍江 大慶 163319 2. 黑龍江省水稻生態(tài)育秧裝置及全程機(jī)械化工程技術(shù)中心，黑龍江 大慶 163319

引 言

近紅外光譜區(qū)(800～2 500 nm)的含氫基團(tuán)的倍頻和合頻吸收峰組成的吸收強(qiáng)度較弱靈敏度較低，吸收帶較寬且嚴(yán)重重疊。若采用全譜建模，不僅會(huì)存在某些光譜區(qū)域與待測(cè)組分相關(guān)性弱，而且相鄰的波長(zhǎng)高度相關(guān)，包含了大量的冗余信息，這都會(huì)影響模型的精度和穩(wěn)健性?？朔@些問題的有效途徑是對(duì)所測(cè)得的光譜進(jìn)行波長(zhǎng)選擇，減少建模所需的波長(zhǎng)點(diǎn)和計(jì)算工作量，進(jìn)而得到預(yù)測(cè)能力強(qiáng)、魯棒性高的模型。在眾多特征波長(zhǎng)選擇算法中，隨機(jī)蛙跳(random frog, RF)[1]是近年來(lái)提出的一種新型特征波長(zhǎng)選擇算法。其依據(jù)不同的變量具有不同的被選擇可能性，通過多次迭代，計(jì)算每個(gè)變量被選擇的概率，選擇概率高的變量為特征波長(zhǎng)。

陳立旦等[2]通過RF選出特征波長(zhǎng)后，建立最小二乘支持向量機(jī)(least squares support vector machine, LS-SVM)模型，對(duì)生物柴油的含水量進(jìn)行預(yù)測(cè)，發(fā)現(xiàn)RF-LS-SVM模型的相關(guān)系數(shù)大于0.95，可以準(zhǔn)確地預(yù)測(cè)生物柴油的含水量。胡孟晗等[3]通過RF對(duì)特征波長(zhǎng)進(jìn)行提取，建立LS-SVM模型預(yù)測(cè)藍(lán)莓硬度和彈性模量，與全譜模型對(duì)比，RF算法可以有效地去除冗余信息，提升模型預(yù)測(cè)準(zhǔn)確率。孫紅等[4]采用相關(guān)系數(shù)法(correlation coefficient, CC)和RF算法篩選對(duì)葉綠素含量敏感的波長(zhǎng)，建立偏最小二乘回歸(partial least squares regression, PLSR)模型對(duì)馬鈴薯作物的葉綠素含量進(jìn)行預(yù)測(cè)，結(jié)果表明RF-PLSR模型預(yù)測(cè)精度優(yōu)于CC-PLSR，可實(shí)現(xiàn)馬鈴薯不同葉位葉綠素含量的無(wú)損檢測(cè)。此外，Yu等[5]采用RF和PLSR建立校準(zhǔn)模型，發(fā)現(xiàn)通過380～1 030 nm區(qū)域的波長(zhǎng)可實(shí)現(xiàn)辣椒植物的總氮含量的預(yù)測(cè)。Zhao等[6]通過RF算法選擇特征波長(zhǎng)，建立RF-PLSR和RF-LS-SVM模型預(yù)測(cè)桑葚果實(shí)的總可溶性固體值含量，兩個(gè)模型皆具有良好的性能。以上結(jié)果表明RF算法在數(shù)據(jù)降維方面是有效的。

盡管RF算法在特征波長(zhǎng)選擇方面具有一定優(yōu)勢(shì)，但存在兩方面的不足：其一是，初始變量集V0的產(chǎn)生是隨機(jī)的，難以保證初始信息的有效性；算法為保證運(yùn)行過程中遍歷整個(gè)數(shù)據(jù)集，要求迭代次數(shù)N需足夠大，從而導(dǎo)致算法的運(yùn)行時(shí)間長(zhǎng)、收斂速度慢。其二是，RF在選擇特征波長(zhǎng)時(shí)，選擇被選概率值大于閾值的變量為特征波長(zhǎng)，但對(duì)閾值的設(shè)定無(wú)理論依據(jù)，易受人為因素影響。

針對(duì)上述兩點(diǎn)，對(duì)RF算法進(jìn)行了改進(jìn)，提出一種聯(lián)合區(qū)間隨機(jī)蛙跳(synergy interval-random frog, Si-RF)算法，以一組公開的土壤樣本近紅外光譜數(shù)據(jù)為例，分別利用RF和改進(jìn)的Si-RF進(jìn)行特征波長(zhǎng)選擇，建立多元線性回歸(multiple linear regression, MLR)模型，比較預(yù)測(cè)精度，并與全譜的PLSR模型進(jìn)行對(duì)比，以證明改進(jìn)的Si-RF算法的有效性。

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

1.1 樣本數(shù)據(jù)

所用數(shù)據(jù)為一組土壤樣本近紅外光譜數(shù)據(jù)，來(lái)自于網(wǎng)站Quality & Technology。該數(shù)據(jù)集包含108個(gè)土壤樣本。樣本光譜的波長(zhǎng)范圍為400～2 500 nm，采樣間隔為2 nm，共計(jì)1 050個(gè)波長(zhǎng)點(diǎn)。本文以土壤有機(jī)質(zhì)(soil organic matter, SOM)的含量作為因變量進(jìn)行波長(zhǎng)選擇及近紅外光譜數(shù)據(jù)建模預(yù)測(cè)分析。

1.2 隨機(jī)蛙跳算法

1.2.1 算法步驟

RF是Li[1]提出的一種類似于可逆跳躍馬爾可夫鏈蒙特卡羅(reversible jump Markov Chain Monte Carlo, RJMCMC)的算法，它以迭代的方式進(jìn)行，計(jì)算每個(gè)變量在每次迭代中被選擇的概率，概率越高變量重要性越大，優(yōu)選概率高的變量為特征變量。

隨機(jī)蛙跳的主要步驟包括以下三步[1]：

(1)初始化：參數(shù)設(shè)置，隨機(jī)選擇一個(gè)包含Q個(gè)變量的變量子集V0；

(2)概率引導(dǎo)模型搜索：基于V0，選擇包含Q*(隨機(jī)產(chǎn)生)個(gè)變量的候選變量子集V*，以一定概率接受V*作為V1，并用V1代替V0，循環(huán)此步驟直至N次迭代完成；

(3)變量評(píng)估：計(jì)算每個(gè)變量被選擇的概率，概率越高變量重要性越大。

其中概率引導(dǎo)模型搜索和變量評(píng)估具體方法如下。

1.2.2 概率引導(dǎo)模型搜索

首先，從均值為Q、方差為0.3Q的正態(tài)分布中隨機(jī)選擇一個(gè)整數(shù)Q*，之后通過以下三種方式之一產(chǎn)生一個(gè)包含Q*個(gè)變量的候選變量子集V*：

(1)如果Q*=Q，則令V*=V0。

(2)如果Q*

(3)如果Q*>Q，則從V-V0(V代表包含全部p個(gè)變量的集合)中隨機(jī)抽取ω(Q*-Q)個(gè)變量，ω默認(rèn)值為3，生成一個(gè)變量子集T，通過V0和T的組合建立PLS模型，保留模型中回歸系數(shù)最大的Q*個(gè)變量，并將其設(shè)為候選子集V*。

簡(jiǎn)而言之，利用所提出的正態(tài)分布控制變量數(shù)，實(shí)現(xiàn)變量的增、刪操作。在得到候選變量子集V*后，下一步是確定V*是否可以被接受。分別對(duì)V0和V*建立PLS模型，計(jì)算交叉驗(yàn)證均方根誤差(cross-validation root mean square error, RMSECV)，得到RMSECV和RMSECV*。如果RMSECV*≤RMSECV，接受V*為V1，否則接受V*為V1概率為0.1RMSECV/RMSECV*。最后，使用V1中的變量更新V0，并重復(fù)N次迭代，直至循環(huán)結(jié)束。

1.2.3 變量評(píng)估

N次迭代之后，總共獲得N個(gè)變量子集。對(duì)于每個(gè)變量，可以使用式(1)計(jì)算其被選擇的概率。

(1)

式(1)中，Nj為第j個(gè)變量在N次迭代中被選擇的次數(shù)，變量越重要，被這N個(gè)變量子集選擇的機(jī)會(huì)就越多。因此，該選擇概率可以用作變量重要性的度量，可以用作變量選擇的標(biāo)準(zhǔn)。

1.3 對(duì)RF算法的改進(jìn)

1.3.1V0子集的初選

在RF算法中，初始變量集V0的產(chǎn)生是隨機(jī)的，具有較大的不確定性，可能會(huì)產(chǎn)生無(wú)信息變量或干擾信息，從而導(dǎo)致算法的迭代次數(shù)大，運(yùn)行時(shí)間長(zhǎng)。為了提高初始集V0變量的有效性，減少迭代次數(shù)，對(duì)V0子集的產(chǎn)生進(jìn)行改進(jìn)。

聯(lián)合區(qū)間偏最小二乘法(synergy interval partial least squares, SiPLS)是Norgaard提出的一種波長(zhǎng)選擇算法。該方法將光譜劃分為等寬的n個(gè)子區(qū)間，對(duì)其中m個(gè)子區(qū)間任意組合為聯(lián)合區(qū)間。基于聯(lián)合區(qū)間建立PLS模型，比較各PLS模型的RMSECV的值，將最小RMSECV值所對(duì)應(yīng)的聯(lián)合區(qū)間的波長(zhǎng)設(shè)為初始變量集V0，開始迭代，可以消除V0的隨機(jī)性，避免無(wú)信息變量及噪聲的干擾，從而減少迭代次數(shù)。

1.3.2 建模波長(zhǎng)的優(yōu)選

在RF算法中，一般選擇概率值較大的前10或15個(gè)變量，或者通過人為設(shè)定概率的閾值，取概率值大于閾值的變量來(lái)選擇符合要求的特征波長(zhǎng)，建模波長(zhǎng)數(shù)量選擇存在不確定性。

本文的改進(jìn)是：對(duì)排序后的變量從第一個(gè)波長(zhǎng)開始，每次增加一個(gè)波長(zhǎng)，建立光譜數(shù)據(jù)和有機(jī)質(zhì)含量數(shù)據(jù)之間的MLR模型。計(jì)算每個(gè)模型的驗(yàn)證均方根誤差(root mean square error of validation, RMSEV)值，其中最小RMSEV值所對(duì)應(yīng)的變量子集即為特征波長(zhǎng)。RMSEV可以使用式(2)計(jì)算

(2)

這樣可以找到預(yù)測(cè)精度最優(yōu)所包含的波長(zhǎng)數(shù)，提高預(yù)測(cè)精度。

1.4 建模方法

現(xiàn)有研究大多對(duì)RF所選特征波長(zhǎng)建立PLSR模型。而MLR是一種常規(guī)的校正方法，直觀簡(jiǎn)單，且具有良好的統(tǒng)計(jì)特性，應(yīng)用非常普遍，其優(yōu)點(diǎn)是產(chǎn)生的模型比主成分回歸(principal components regression, PCR)和PLSR模型更簡(jiǎn)單，更易于解釋。

本工作建立三種模型：基于全譜的PLSR模型、基于RF波長(zhǎng)選擇的MLR模型和基于Si-RF改進(jìn)的波長(zhǎng)選擇的MLR模型。通過三種模型預(yù)測(cè)能力的比較驗(yàn)證本法的有效性。模型的預(yù)測(cè)能力主要通過校正相關(guān)系數(shù)(Rc)、校正均方根誤差(RMSEC)、預(yù)測(cè)相關(guān)系數(shù)(Rp)、預(yù)測(cè)均方根誤差(RMSEP)指標(biāo)來(lái)評(píng)價(jià)。其中，R取值越接近1，RMSEC和RMSEP越接近0，模型的擬合性越好，預(yù)測(cè)精度越高。

1.5 數(shù)據(jù)分析

軟件采用MATLAB R2015b及The Unscrambler X 10.3 (64-bit)，光譜數(shù)據(jù)的預(yù)處理、建模分析及預(yù)測(cè)在Unscrambler軟件中實(shí)現(xiàn)，特征波長(zhǎng)提取、圖形的繪制在MATLAB中實(shí)現(xiàn)。計(jì)算機(jī)硬件的配置為Intel(R)Core(TM)i5-3450CPU@3.50GHz處理器，8GB內(nèi)存，操作系統(tǒng)為windows10。

2 結(jié)果與討論

2.1 光譜數(shù)據(jù)特征

土壤樣本的原始近紅外光譜圖如圖1(a)所示。為校正光譜基線，消除其他背景的干擾，提高光譜分辨率，并且在一定程度上減少各變量間的線性相關(guān)性，利用Savitzky-Golay窗口寬度為11的一階求導(dǎo)法對(duì)原始光譜數(shù)據(jù)進(jìn)行預(yù)處理，預(yù)處理后的近紅外光譜圖如圖1(b)所示，可以發(fā)現(xiàn)通過預(yù)處理后的近紅外光譜曲線，能更精確地確定吸收峰的位置。

圖1 原始光譜圖及預(yù)處理后的光譜圖(a): 原始光譜圖；(b): S-G一階導(dǎo)處理后的光譜圖Fig.1 Original and pre-processed spectra(a): Original; (b): S-G first derivative

將108個(gè)土壤樣本通過SPXY(sample set portioning based on joint x-y distance)算法分為75%訓(xùn)練集和25%預(yù)測(cè)集，建模集包含81個(gè)樣本，預(yù)測(cè)集包含27個(gè)樣本，土壤有機(jī)質(zhì)含量統(tǒng)計(jì)數(shù)據(jù)結(jié)果如表1所示。劃分后的建模集的SOM含量范圍涵蓋預(yù)測(cè)集的SOM含量，建模集具有代表性。

表1 土壤有機(jī)質(zhì)含量統(tǒng)計(jì)數(shù)據(jù)結(jié)果Table 1 Statistical data of soil organic matter content

2.2 特征波長(zhǎng)選取

2.2.1 RF變量選擇結(jié)果

如前所述，首先對(duì)RF進(jìn)行初始化參數(shù)設(shè)置，N設(shè)定為10 000，Q設(shè)定為10，開始運(yùn)行。每個(gè)變量被選擇的概率結(jié)果如圖2所示，選擇概率大于0.2的變量為最終特征波長(zhǎng)，得到滿足條件的有10個(gè)波長(zhǎng)點(diǎn)分別為1 420，1 390，1 392，1 394，1 388，1 422，2 318，1 424，1 396和1 922 nm。

圖2 RF運(yùn)行結(jié)果Fig.2 The result of random frog

2.2.2 Si-RF變量選擇結(jié)果

表2 SiPLS子區(qū)間優(yōu)選結(jié)果Table 2 Sub-interval optimization results of SiPLS

由表2可以發(fā)現(xiàn)，將全譜等分為30個(gè)區(qū)間，組合數(shù)設(shè)置為3時(shí)，RMSECV最小，此時(shí)所選的特征波長(zhǎng)點(diǎn)為104個(gè)，將這三個(gè)波段1 182～1 250，1 392～1 460和2 288～2 354 nm，共計(jì)104個(gè)波長(zhǎng)點(diǎn)作為初始變量子集V0，RF算法的迭代次數(shù)分別設(shè)置為500，1 000，1 500和2 000次，得到結(jié)果如表3所示。

表3 不同迭代次數(shù)的優(yōu)選結(jié)果Table 3 Optimal results of different iteration times

由表3可知，當(dāng)N設(shè)置為1 000次時(shí)，RMSEV值最小。該情況下Si-RF運(yùn)行結(jié)果如圖3所示，每個(gè)變量被選擇的概率結(jié)果如圖3(a)所示。將每個(gè)變量被選擇的概率值進(jìn)行降序排列，從第一個(gè)波長(zhǎng)開始，逐次增加一個(gè)波長(zhǎng)建立MLR模型。各模型的RMSEV值如圖3(b)所示，正方形標(biāo)記所示為最低RMSEV值，為0.818 4，此時(shí)選擇的特征波長(zhǎng)數(shù)為17個(gè)，分別為1 392，1 394，1 420，2 332，2 330，1 418，1 440，1 348，1 920，1 402，2 000，1 424，2 312，1 442，1 426，1 444和2 364 nm。

圖3 Si-RF運(yùn)行結(jié)果(a): 各變量被選概率；(b): 各模型RMSEV值Fig.3 The result of Si-RF(a): Selection probability of each variable;(b): RMSEV values of each model

2.3 模型建立與比較

將全譜、RF以及 Si-RF選擇的特征波長(zhǎng)，建立回歸模型比較預(yù)測(cè)能力，得到模型的校正、預(yù)測(cè)相關(guān)系數(shù)和校正、預(yù)測(cè)均方根誤差的值如表4所示。

表4 不同波長(zhǎng)選擇方法下模型的結(jié)果Table 4 Results of model with different wavelength selection methods

從表4可以看出，RF和Si-RF模型的各項(xiàng)參數(shù)均優(yōu)于全譜，改進(jìn)的Si-RF算法模型的各項(xiàng)參數(shù)均優(yōu)于RF。基于RF所選特征波長(zhǎng)的MLR模型的Rp為0.9354，RMSEP為1.627 6，而改進(jìn)后Si-RF選擇的特征波長(zhǎng)MLR模型的Rp為0.984 8，RMSEP減小到0.818 4，大大提升了預(yù)測(cè)精度。

圖4分別為對(duì)建模集、預(yù)測(cè)集樣本的全譜-PLS、RF-MLR和Si-RF-MLR模型的SOM的實(shí)測(cè)值和預(yù)測(cè)值相關(guān)圖。從圖中可以更加直觀的看出，基于Si-RF波長(zhǎng)選擇算法的MLR模型優(yōu)于全譜模型及RF算法的MLR模型。

圖4 不同模型下土壤有機(jī)質(zhì)的實(shí)測(cè)值和預(yù)測(cè)值相關(guān)圖(a): 全譜-PLS；(b): RF-MLR；(c): Si-RF-MLRFig.4 Correlation between measured and predicted values of SOM obtained from different models(a): Full spectrum PLS; (b): RF-MLR; (c): Si-RF-MLR

由于RF算法對(duì)初始變量集的產(chǎn)生是隨機(jī)的，有較大的不確定性，可能會(huì)包含無(wú)信息變量或干擾信息，從而導(dǎo)致算法的迭代次數(shù)大、運(yùn)行時(shí)間長(zhǎng)。而通過SiPLS特征波長(zhǎng)初選，得到的波長(zhǎng)對(duì)于目標(biāo)變量變化最為敏感，同時(shí)避免了其他光譜無(wú)信息變量與噪聲的影響。所以首先對(duì)全譜通過SiPLS進(jìn)行特征波長(zhǎng)初選，將其初選結(jié)果作為RF的初始變量子集V0，這樣可以改善RF收斂速度慢的問題，減少RF算法的迭代次數(shù)，大大節(jié)省運(yùn)行時(shí)間，并且由于初始變量子集是針對(duì)于有效信息的波長(zhǎng)，有利于RF每次迭代中V*所包含的波長(zhǎng)的選擇，可以提高預(yù)測(cè)精度。在運(yùn)行中，迭代次數(shù)也由10 000次減少至1 000次，提高運(yùn)行效率。

通過Si-RF選出的特征波長(zhǎng)點(diǎn)的范圍在1 348～1 444，1 920～2 364 nm之間，這與許多前人研究所選波長(zhǎng)點(diǎn)范圍基本一致。如：白婷等[7]針對(duì)艾比湖60個(gè)表層土樣，基于CARS算法提取的SOM特征波段主要集中在1 970和2 340 nm附近；朱亞星等[8]通過UVE-CARS優(yōu)選出84個(gè)變量做為預(yù)測(cè)SOM含量的特征波長(zhǎng)，分布于561～721和1 920～2 280 nm波段；于雷等[9]通過CARS-SPA優(yōu)選出的37個(gè)特征波長(zhǎng)，集中在近紅外區(qū)域1 800～2 400 nm，而且基于波長(zhǎng)選擇建立的SOM含量的PLSR模型預(yù)測(cè)精度最優(yōu)。本工作Si-RF優(yōu)選出的波段與圖2B近紅外光譜曲線吸收峰的位置也基本一致，符合高志海等[10]的論點(diǎn)，即光譜曲線上的凸起區(qū)可能對(duì)提取土壤有機(jī)質(zhì)信息有實(shí)際意義。

對(duì)比RF及Si-RF所選波長(zhǎng)點(diǎn)范圍，RF的范圍在1 388～1 424和1 922～2 318 nm之間，Si-RF的范圍在1 348～1 444和1 920～2 364 nm之間，可以發(fā)現(xiàn)Si-RF已經(jīng)基本涵蓋RF所選波長(zhǎng)的大部分，這也在一定程度上說(shuō)明可以減少算法迭代次數(shù)。

3 結(jié) 論

提出了一種近紅外光譜分析中特征波長(zhǎng)選擇的Si-RF算法，該方法通過對(duì)全譜進(jìn)行SiPLS特征波長(zhǎng)初選，將所得的波長(zhǎng)做為初始變量子集，使得初始變量子集涵蓋有效信息，以解決RF中迭代次數(shù)過多，運(yùn)行效率較低的問題。將RF和改進(jìn)的Si-RF應(yīng)用于一組土壤樣本近紅外光譜數(shù)據(jù)集，將由RF選擇的特征波長(zhǎng)和改進(jìn)的Si-RF選擇的特征波長(zhǎng)提取出來(lái)，建立MLR模型，發(fā)現(xiàn)Si-RF-MLR模型的預(yù)測(cè)精度優(yōu)于RF-MLR，并且在運(yùn)行時(shí)間上也大大降低，提高運(yùn)行效率；相較于全譜的PLSR模型，也極大的提高了預(yù)測(cè)精度，簡(jiǎn)化模型的復(fù)雜度。證明改進(jìn)的Si-RF是一種有效的特征波長(zhǎng)選擇算法。