謝德紅，李俊鋒，劉 菂，萬曉霞，葉 藝

1. 南京林業(yè)大學(xué)輕工與食品學(xué)院，江蘇 南京 210037 2. 河南牧業(yè)經(jīng)濟(jì)學(xué)院包裝與印刷工程學(xué)院，河南 鄭州 450046 3. 北京工商大學(xué)食品安全大數(shù)據(jù)技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100048 4. 武漢大學(xué)湖北省文物顏色信息數(shù)字化與虛擬再現(xiàn)工程研究中心，湖北 武漢 430079

引 言

農(nóng)藥的高效防治功能，對促進(jìn)農(nóng)作物增產(chǎn)、緩解農(nóng)產(chǎn)品供需矛盾有舉足輕重的作用。然而，農(nóng)藥大量不合理使用，導(dǎo)致農(nóng)藥殘留超過食品安全標(biāo)準(zhǔn)，引發(fā)諸多食品安全問題[1]。當(dāng)前一些化學(xué)、生物等檢測方法精度雖高，但存在耗時費(fèi)力、破損被測物、污染環(huán)境等缺點(diǎn)，難以滿足現(xiàn)代農(nóng)業(yè)快速、無損、批量及實(shí)時檢測的需求。近紅外光譜因其無損、快速、無污染的特點(diǎn)，在農(nóng)藥殘留定性和定量檢測應(yīng)用中優(yōu)勢明顯。但是，近紅外光譜是典型的弱信號，被測物化學(xué)組分的譜峰相對微弱且重疊度高。此外，由于近紅外光譜采集儀器和采集環(huán)境等因素，實(shí)測近紅外光譜常受噪聲干擾。在后續(xù)近紅外光譜數(shù)據(jù)分類中，無法有效區(qū)別此光譜中被測物化學(xué)組分信息和噪聲，進(jìn)而影響分類精度。有效地去除光譜數(shù)據(jù)中的噪聲，對建立穩(wěn)定性好、分類精度高的農(nóng)藥殘留檢測的分類模型十分重要。

小波分解[2-3]和經(jīng)驗(yàn)?zāi)Ｌ纸?empirical mode decomposition, EMD)[4-5]可較好刻畫光譜信號的特征，近些年常被用于近紅外等各類光譜信號的降噪中，并在信噪比小時有一定抑制噪聲的效果。小波分解方法是在小波分解的基礎(chǔ)上通過閾值手段去除噪聲，因此小波分解的基函數(shù)(即小波基)和閾值選擇與染噪光譜特征的適應(yīng)性對降噪至關(guān)重要。有研究[3]提出利用Stein無偏風(fēng)險估計(jì)(Stein’s unbiased estimate of risk, SURE)有效地解決了閾值自適應(yīng)問題，但小波基與染噪光譜特征的匹配問題卻鮮有研究關(guān)注。據(jù)文獻(xiàn)[6]研究發(fā)現(xiàn)，當(dāng)小波基與被分解信號特征不匹配時，重構(gòu)信號易產(chǎn)生Gibbs震蕩現(xiàn)象。此震蕩波破壞了光譜的幾何特征，導(dǎo)致在后續(xù)近紅外光譜數(shù)據(jù)分類中易被誤認(rèn)為被測物化學(xué)組分的譜峰，惡化分類精度、影響微弱的農(nóng)藥殘留成分的正確分析。EMD克服了小波基與被分解信號相適應(yīng)的要求，可將任意特征的光譜信號分解成不同的震蕩成分，即本征模態(tài)函數(shù)(intrinsic mode functions, IMFs)，但在信號重構(gòu)時易產(chǎn)生混疊模態(tài)問題?；殳B模態(tài)是由經(jīng)驗(yàn)?zāi)７纸獾倪^程中信號因局部極值在很短的時間(或很短的波段)間隔內(nèi)發(fā)生多次跳變導(dǎo)致的。在光譜信號染噪后，混疊模態(tài)中的信號狀態(tài)很難完全避免。為此，有學(xué)者提出了集合模態(tài)分解(ensemble empirical mode decomposition，EEMD)[4]。此方法是一種噪聲協(xié)助分析方法，因而當(dāng)分解出的IMFs數(shù)量比較有限時，各IMFs不可避免的存在噪聲，導(dǎo)致重構(gòu)信號的信噪比降低?；パa(bǔ)集合模態(tài)分解(complete ensemble empirical mode decomposition, CEEMD)[5]可抑制IMFs里的殘余噪聲，但當(dāng)參數(shù)選擇不當(dāng)時，會產(chǎn)生錯誤成分使最后獲得的IMFs不能真正滿足IMFs的定義，導(dǎo)致重構(gòu)信號仍會存在一些小震蕩，在近紅外光譜中呈現(xiàn)為額外的譜峰。為了提高近紅外光譜的信噪比、減少光譜信號的失真以及額外譜峰的產(chǎn)生或殘留、滿足果蔬農(nóng)藥殘留定性檢測的應(yīng)用需求，本研究提出了L曲線與Hodrick-Prescott分解模型相結(jié)合的自適應(yīng)降噪方法，復(fù)原被測物表面原始近紅外光譜，為后續(xù)近紅外光譜數(shù)據(jù)分類、農(nóng)藥殘留檢測、及農(nóng)藥殘留化學(xué)組分信息分析和挖掘奠定基礎(chǔ)。

1 改進(jìn)Hodrick-Prescott降噪的原理與方法

Hodrick-Prescott分解模型[7]是一種信號分離方法，廣泛用于經(jīng)濟(jì)數(shù)據(jù)中經(jīng)濟(jì)趨勢和各時間點(diǎn)震蕩波的分離。在果蔬農(nóng)藥殘留的近紅外光譜中，用以表征被測物化學(xué)組分的譜峰相對較弱且時常發(fā)生重疊，因而其光譜波形隨波長的變化相對緩慢，如同經(jīng)濟(jì)趨勢的特征； 通常，近紅外光譜采集過程中所產(chǎn)生的噪聲隨波長變化相對劇烈，如同經(jīng)濟(jì)數(shù)據(jù)中隨時間點(diǎn)變化的震蕩波特征。因此將依據(jù)染噪近紅外光譜的特征，建立基于Hodrick-Prescott分解模型的近紅外光譜降噪的方法，并與bior6.8小波分解方法、sym8小波分解方法、互補(bǔ)集合模態(tài)分解方法進(jìn)行比較，驗(yàn)證本方法的有效性。

1.1 Hodrick-Prescott降噪原理與方法

設(shè)被測物的原始近紅外光譜(即未染噪聲的近紅外光譜)為f=(fi，…，fj)∈Rn，其近紅外光譜在采集過程中受廣義高斯噪聲的污染獲得染噪近紅外光譜為s=(si，…，sj)T∈Rn，它們與噪聲e=(ei，…，ej)T∈Rn，關(guān)系如下

s=f+e

(1)

下標(biāo)i和j表示波長，n為近紅外波長范圍內(nèi)的波長點(diǎn)的數(shù)量。

利用Hodrick-Prescott分解方法降噪是通過引入一定約束條件，將光譜信號復(fù)原為具有存在性、唯一性且適定的結(jié)果，且此結(jié)果為受噪聲干擾小的過程。已知染噪近紅外光譜信號s，求原始近紅外光譜f的Hodrick-Prescott分解方法可轉(zhuǎn)化為目標(biāo)函數(shù)J(f)的最小化問題

J(f)=(1/2)‖s-f‖2+κ‖Df‖2

(2)

此目標(biāo)函數(shù)由殘差項(xiàng)和正則化項(xiàng)兩部分構(gòu)成。其中，‖s-f‖2為殘差項(xiàng)，保證復(fù)原信號的保真度； ‖Df‖2為正則化項(xiàng)，約束近紅外光譜f的波形特征；D(n-2)×n為二階差分矩陣

(3)

用于懲罰復(fù)原光譜信號的震蕩波動，迫使其趨向梯度減少的方向?！ぁ?表示l2范數(shù)。κ≥0為正則化參數(shù)，用以調(diào)整目標(biāo)函數(shù)兩部分的權(quán)重。若κ=0，正則化項(xiàng)的懲罰作用失效，求得的復(fù)原光譜信號f趨于染噪光譜信號; 反之，若κ→+∞, 則正則化項(xiàng)懲罰起主導(dǎo)作用，求得的復(fù)原光譜信號f只會盡可能滿足‖Df‖2最小，因而其復(fù)原光譜信號逼真程度比較差、偏離原始光譜信號，甚至無法呈現(xiàn)原始光譜信號的基本波形。由此可知，選擇合適的正則化參數(shù)κ，對獲得最優(yōu)降噪效果至關(guān)重要。

1.2 正則化參數(shù)自適應(yīng)獲取方法

依上所析，正則化參數(shù)κ越大、懲罰越嚴(yán)重，復(fù)原出的光譜信號越平滑、含震蕩波動程度越弱。雖易去除噪聲，但也存在平滑掉有利于數(shù)據(jù)分類的小譜峰的危險； 反之，正則化參數(shù)κ越小，正則化項(xiàng)的懲罰作用越小，復(fù)原光譜信號則越易殘留噪聲。此外，原始Hodrick-Prescott方法中的正則化參數(shù)κ的選擇是非自適應(yīng)的，選擇一個固定值，難以適應(yīng)實(shí)測近紅外光譜非穩(wěn)定特征及實(shí)測近紅外光譜信噪比有大有小的狀況。基于此，為了使正則化參數(shù)與實(shí)測染噪近紅外光譜能自適應(yīng)，提出利用L-曲線方法[8]獲得Hodrick-Prescott方法中的自適應(yīng)值。

L-曲線是指各正則化參數(shù)下對應(yīng)復(fù)原光譜信號的殘差項(xiàng)‖s-fκ‖2與正則化項(xiàng)‖Df‖2所構(gòu)成的log-log尺度曲線，其形狀接近L形。

l:={(log‖s-fκ‖2, log‖Dfκ‖2):κ≥0}

(4)

其中，“l(fā)”表示L-曲線，fκ為每個κ值下的復(fù)原光譜信號。由式可知，每個κ值均存在一個最優(yōu)fκ使目標(biāo)函數(shù)J(f)趨于最小。由此，L-曲線上每個點(diǎn)則表示每個κ值下的最小J(f)值，其拐點(diǎn)表示所有最小J(f)值中最小值，且拐點(diǎn)處對應(yīng)的κ值則為與當(dāng)前染噪光譜信號最優(yōu)適應(yīng)值。

2 結(jié)果與討論

實(shí)驗(yàn)采集341個上海青(一種青菜)近紅外光譜樣本，并在采集的近紅外光譜上人工添加廣義高斯噪聲作為降噪對象。為驗(yàn)證本方法的有效性，選用bior6.8和sym8小波基的SURE閾值方法[3]及CEEMD方法[5]作對比，并采用信噪比(signal noise ratios, SNR)、均方根誤差(root mean square error, RMSE)及降噪光譜數(shù)據(jù)訓(xùn)練所得分類模型識別率評價降噪效果。

2.1 實(shí)驗(yàn)數(shù)據(jù)

以無錫迅杰光遠(yuǎn)科技有限公司的IAS-2000近紅外光譜儀為近紅外光譜采集設(shè)備，光譜儀采集參數(shù)設(shè)置如下： 平均采樣次數(shù)為30次，采集波段為900～1 700 nm，測量間隔為1 nm。光譜采集環(huán)境： 裝有空調(diào)的恒溫環(huán)境(23 ℃)內(nèi)。采集的實(shí)驗(yàn)樣品為無農(nóng)藥殘留上海青和噴灑了樂果(農(nóng)藥，C5H12NO3PS2)的上海青。其中，無農(nóng)藥殘留上海青購自農(nóng)貿(mào)市場，并通過小蘇打和食鹽浸泡、自來沖、蒸餾水先后沖洗晾干； 殘留農(nóng)藥上海青則由無農(nóng)藥殘留上海青噴曬不同濃度配比的樂果溶液并晾干所得。其中，樂果為河南省周口市德貝爾生物化學(xué)品工程有限生產(chǎn)的40%樂果乳油； 不同濃度配比的樂果溶液由蒸餾水稀釋所得，且稀釋的最小和最大倍數(shù)分別為50倍和1 500倍。實(shí)驗(yàn)分別采集殘留農(nóng)藥樣本291個、未殘留農(nóng)藥樣本50個。從農(nóng)藥殘留樣本和未殘留農(nóng)藥樣本中分別選取146個農(nóng)藥殘留樣本和25個未殘留農(nóng)藥樣本構(gòu)成訓(xùn)練樣本集，剩余為測試樣本集。其中，146個農(nóng)藥殘留樣本上的樂果濃度橫跨最小與最大配比。

為了便于驗(yàn)證降噪效果，設(shè)定上述采集的近紅外吸收光譜為原始純光譜信號，模擬疊加不同程度高斯噪聲后的光譜為實(shí)驗(yàn)待處理的染噪光譜。圖1上排圖為稀釋50倍的農(nóng)藥殘留樣本，疊加噪聲后的近紅外光譜信噪比分別為19.07 dB[如圖1(a)]和18.79 dB[如圖1(b)]。

圖1 染噪近紅外光譜和其L-曲線(a)： SNR=19.07 dB； (b)： SNR=18.79 dBFig.1 Noise-contaminated NIR spectroscopy and their L-curve(a)： SNR=19.07 dB； (b)： SNR=18.79 dB

2.2 降噪效果與分析

對任一組染噪近紅外光譜，可利用L-曲線方法獲得此光譜的自適應(yīng)κ值，其步驟執(zhí)行如下： 先設(shè)定κ的范圍和步長，設(shè)定κ=0.001, 0.002, …, 29.999, 30； 然后分別計(jì)算各κ值下殘差項(xiàng)‖s-fκ‖2與正則化項(xiàng)‖Df‖2的值，在獲得L-曲線的基礎(chǔ)上，獲取曲線最大曲率點(diǎn)處的κ值； 最后，在此κ值下，利用Hodrick-Prescott方法復(fù)原的近紅外光譜fκ即為降噪后的近紅外光譜。依上述步驟，圖1中兩組染噪近紅外光譜的自適應(yīng)正則化參數(shù)分別為0.91和0.11，如圖1(a,b)下排圖所示。

為了衡量、比較降噪效果，首先采用了SNR和RMSE兩個指標(biāo)。其中，SNR值大小與降噪效果成正比，RMSE值大小則與降噪效果成反比。表1利用SNR和RMSE兩個指標(biāo)比較了bior6.8和sym8小波的SURE閾值方法(均為五層小波分解)[3]、CEEMD方法[5]、以及本方法。首先，通過SNR和RMSE值的橫向比較發(fā)現(xiàn)，本方法的SNR值最大、RMSE值最小，說明相對其他三種方法，本方法的降噪效果最優(yōu)。其次，結(jié)合兩組染噪近紅外光譜的縱向數(shù)據(jù)比較發(fā)現(xiàn)，降噪方法優(yōu)劣順序未隨噪聲不同而變化，說明本文方法對不同噪聲的降噪效果相對穩(wěn)定。SNR和RMSE均是從全局的角度考察降噪前后光譜相似程度，未能反映局部缺陷，而局部缺陷恰恰容易惡化后續(xù)數(shù)據(jù)的分類精度。例如，當(dāng)降噪后光譜出現(xiàn)沿未染噪光譜上下波動的小波峰波谷，且波峰和波谷構(gòu)成的增加和減少的面積能相互抵消時，SNR和RMSE則不能從數(shù)據(jù)上體現(xiàn)額外小波峰波谷的局部缺陷存在，但此額外產(chǎn)生的小波峰波谷在后續(xù)數(shù)據(jù)分類中又極易誤認(rèn)為是物質(zhì)化學(xué)組分的譜峰，影響分類的精確度。因此，本研究將繼續(xù)從局部和后續(xù)數(shù)據(jù)分類精度進(jìn)一步驗(yàn)證降噪效果。

當(dāng)被測物含量有限時，其近紅外光譜(或其對應(yīng)化學(xué)組分的近紅外光譜譜峰)也比較弱。例如，在本實(shí)驗(yàn)中，上海青本身的化學(xué)組分含量一般遠(yuǎn)高于殘留農(nóng)藥的化學(xué)組分，因而農(nóng)藥殘留上海青的近紅外光譜中除較弱的、反映農(nóng)藥化學(xué)組分的譜峰外，主要為寬波段特征、且信號較強(qiáng)的上海青化學(xué)組分的近紅外光譜。此特點(diǎn)使得殘留樣本的近紅外光譜與未殘留樣本的近紅外光譜識別度小。為了增強(qiáng)農(nóng)藥譜峰的可識別度，在近紅外光譜數(shù)據(jù)分類、判斷農(nóng)藥殘留與否的過程中，常用導(dǎo)數(shù)處理方式抑制背景光譜信號(如上海青的近紅外光譜)、增強(qiáng)目標(biāo)識別物化學(xué)組分的譜峰(如農(nóng)藥化學(xué)組分的譜峰)。一階導(dǎo)數(shù)的處理方式是雙刃劍，會增強(qiáng)農(nóng)藥化學(xué)組分譜峰也會放大噪聲等非光譜成份，但均有利于從局部判斷降噪后光譜與原未染噪光譜的相似度。因此，利用光譜一階導(dǎo)數(shù)譜圖進(jìn)一步評判各方法的降噪效果，如圖2所示。

圖2 降噪后光譜和未染噪聲光譜的一階導(dǎo)數(shù)譜圖(a): 小波法(bior6.8小波核)； (b): 小波法(sym8)； (c): CEEMD法； (d): 本方法Fig.2 First derivative spectrum of the denoised NIR spectroscopy and its original(a): Wavelet method (with bior6.8 wavelet kernel); (b): Wavelet method (with sym8 wavelet kernel);(c): CEEMD method; (d): The proposed method

圖2各分圖均為圖1(a)對應(yīng)染噪近紅外光譜經(jīng)過降噪和一階求導(dǎo)后的結(jié)果。圖中實(shí)線代表降噪后光譜的一階求導(dǎo)結(jié)果，虛線則代表原未染噪聲光譜的一階求導(dǎo)結(jié)果。在圖2(a)中，實(shí)線圖譜大約于980，1 020，1 100，1 290，1 410，1 500和1 620 nm等波段處呈現(xiàn)明顯的震蕩波，而此震蕩波在虛線的譜圖上均不可見。由此推測，此額外震蕩波為bior6.8小波方法降噪所產(chǎn)生、非被測物化學(xué)組分，會干擾后續(xù)光譜數(shù)據(jù)的分類、惡化分類精度； 圖2(b)顯示sym8小波方法也會產(chǎn)生較強(qiáng)的震蕩波； 圖2(c)實(shí)線圖譜中雖無明顯的震蕩波，但相對于虛線所示譜圖，仍然呈現(xiàn)額外的小波峰波谷，此小波峰波谷與某些物質(zhì)化學(xué)組分的近紅外二級和三級倍頻波[9]難以區(qū)分，也極易被誤當(dāng)作上海青或農(nóng)藥的化學(xué)組分的小譜峰； 圖2(d)中實(shí)線與虛線表示的譜圖比較接近，且無明顯額外的震蕩波和小譜峰。由此分析可知，從局部特征判斷，相對上述其他三種降噪方法，本方法降噪后的效果優(yōu)，也更有益于后續(xù)化學(xué)組分分析以及農(nóng)藥殘留檢測為目的的數(shù)據(jù)分類。

表1 降噪前后近紅外光譜信噪比與均方根誤差對比

Table 1 Comparison of signal-to-noise ratio and the root mean square error of the denoised NIR spectroscopy

染噪光譜SNR評價指標(biāo)sym8bior6.8CEEMD本方法19.07SNR32.6332.5633.0735.75RMSE0.011 70.011 80.011 10.008 218.79SNR31.6931.5031.7133.35RMSE0.013 00.013 30.013 00.010 8

表2 降噪前后近紅外數(shù)據(jù)的SVM分類結(jié)果Table 2 SVM classification results baseddifferent denoising methods

研究近紅外光譜降噪的主要目的是減少噪聲對后續(xù)農(nóng)藥殘留檢測中數(shù)據(jù)分類的干擾，從而提高訓(xùn)練所得分類模型的擬合能力和泛化能力。由此，對染噪訓(xùn)練樣本數(shù)據(jù)和測試樣本數(shù)據(jù)降噪處理后，再經(jīng)過一階導(dǎo)數(shù)預(yù)處理、主成分分析(PCA)選取十維特征的前提下，建立檢測農(nóng)藥殘留與否的支持向量機(jī)(SVM)[10]分類模型，并分別通過訓(xùn)練樣本數(shù)據(jù)和測試樣本數(shù)據(jù)獲得此分類模型的識別率，以進(jìn)一步評價降噪方法優(yōu)劣。表2比較了原未染噪光譜及四種方法降噪后光譜的分類效果。由分類原理可知，上述額外震蕩波和小譜峰的存在，會導(dǎo)致分類數(shù)據(jù)在特征空間的分布結(jié)構(gòu)發(fā)生變化，不僅影響訓(xùn)練所得分類模型擬合數(shù)據(jù)的能力，更會惡化分類模型的泛化能力(即正確識別非訓(xùn)練近紅外光譜數(shù)據(jù)所屬類別的能力)。表2中的數(shù)據(jù)顯示了與理論推導(dǎo)一致的結(jié)論。由表2中的分類模型的識別率還可發(fā)現(xiàn)，本方法降噪后的光譜數(shù)據(jù)訓(xùn)練所得的SVM分類模型，其訓(xùn)練集和測試集的識別率遠(yuǎn)高于其他三種方法的結(jié)果，接近原始無噪聲近紅外光譜數(shù)據(jù)的結(jié)果。此結(jié)果表明，本方法降噪效果優(yōu)良，且對分類模型的準(zhǔn)確率影響最優(yōu)。

3 結(jié) 論

提出了一種改進(jìn)Hodrick-Prescott分解模型的自適應(yīng)近紅外光譜降噪方法，充分利用Hodrick-Prescott分解模型中正則化項(xiàng)對復(fù)原光譜的懲罰作用，迫使復(fù)原光譜傾向于梯度減少的方向，并結(jié)合L-曲線方法自適應(yīng)地獲取了Hodrick-Prescott分解模型中正則化參數(shù)，從而實(shí)現(xiàn)了染噪近紅外光譜自適應(yīng)降噪的目的。實(shí)驗(yàn)從信噪比、均方根誤差及分類模型識別率數(shù)據(jù)證明了本降噪方法具有一定的優(yōu)越性，可為近紅外光譜定性檢測提供可靠的噪聲預(yù)處理方法。