張根偉，彭思龍，曹樹亞*，趙 將，楊 柳， 楊 杰，楊俊超，黃啟斌*

1. 軍事科學院防化研究院，北京 102205 2. 國民核生化災害防護國家重點實驗室，北京 102205 3. 中國科學院自動化研究所，北京 100190 4. 中國科學院大學，北京 100190

引 言

離子遷移譜(ion mobility spectrometry, IMS)技術是目前最為廣泛使用的痕量化學物質檢測技術之一。離子遷移譜檢測儀具有常壓下檢測、分析靈敏度高、響應時間短、成本低等優(yōu)點，廣泛應用于化學戰(zhàn)劑、爆炸物和毒品等檢測領域[1-5]。

IMS檢測儀產(chǎn)生的原始信號具有如下的特點： 信號?。?干擾信號多，其中主要有來自漂移管的高壓電場(約300 V·cm-1)、多路氣泵的PWM控制脈沖信號及其工作時的振動，此外還有來自控制電路等的電磁干擾。因此，采集到的原始離子遷移譜譜圖信號噪聲較大，微弱信號容易淹沒在噪聲中而無法檢測。此外，由于儀器自身及外部環(huán)境干擾，采集到的譜圖基線基本都不在零基線上，不可避免的存在基線漂移現(xiàn)象。Jiang等[6-7]研究了滑動時間數(shù)據(jù)平均、多項式數(shù)據(jù)平滑、傅里葉變換降噪濾波、小波分析降噪方法對實際采集到的離子遷移譜信號數(shù)據(jù)進行處理的降噪效果?；€一般被看作是緩變的背景，經(jīng)過譜圖的下方但是不能越過波峰，它無任何明確的數(shù)學定義。通過對基線賦予不同的先驗知識，能夠建立不同的基線校正模型。目前主要有以下幾類基線校正方法： 導數(shù)與光滑性方法、曲線擬合方法，非對稱最小二乘方法、背景估計方法及稀疏性約束方法等[7-9]。

稀疏表示是過去近十年來信號處理界一個非常引人關注的研究領域，稀疏表示的目的就是在給定的超完備字典中用盡可能少的原子來表示信號，可以獲得信號更為簡潔的表示方式，從而使我們更容易地獲取信號中所蘊含的信息，更方便進一步對信號進行處理。傳統(tǒng)的離子遷移譜譜圖信號去噪和基線校正預處理分別采用不同的算法實現(xiàn)，沒有將兩者結合起來考慮，導致需單獨評價各自預處理效果。本文通過增加基線光滑度約束以及基于Gaussian字典，提出一種基于稀疏表示的可同時實現(xiàn)離子遷移譜譜圖去噪和基線校正的快速算法。

1 方法原理

1.1 稀疏表示模型

離子遷移譜信號x一般認為由三部分構成： 純譜s、基線b和噪聲n組成，x可以寫成

x=s+b+n

(1)

基于稀疏表示，本文提出了一種可同時實現(xiàn)離子遷移譜譜圖去噪及基線校正模型如式(2)所示

(2)

純譜s由超完備字典Φ組成，α為相應的表達系數(shù)，即s=Φα。這是一個標準的優(yōu)化問題，為了得到稀疏解，采用l1范數(shù)作為線性懲罰項。λ為正則化參數(shù)用來調節(jié)懲罰項在優(yōu)化過程中的比例，同時增加基線平滑度約束條件，D為差分矩陣。通過稀疏表示模型，用最少的字典原子重建離子遷移譜信號，將噪聲項系數(shù)去除實現(xiàn)去噪平滑的效果。

1.2 求解算法

離子遷移譜在弱電場下形成的反應離子峰一般認為是Gaussian分布，Gaussian分布函數(shù)表達式如式(3)

(3)

采用Gaussian分布函數(shù)擬合離子遷移譜反應離子峰，離子遷移譜譜圖可以看成由不同的Gaussian分布函數(shù)擬合而成。超完備字典的構造方法為： 根據(jù)Gaussian分布函數(shù)表達式(3)，將反應峰峰位參數(shù)μ，峰寬參數(shù)σ作為調節(jié)變量，在峰位及峰寬變化區(qū)間內構造超完備字典。

稀疏系數(shù)α的求解有幾種常用的方法，如交替方向乘子算法(ADMM)[10]、最小角回歸法(LARS)[11]、可分離替代函數(shù)法(SSF)[12]，迭代重加權最小二乘算法(IRLS)[13]等。綜合考慮算法運行速度及收斂性，本文采用迭代重加權最小二乘算法求解稀疏表示模型，同時實現(xiàn)去噪及基線校正功能，算法流程如下：

算法: 迭代重加權最小二乘法(IRLS)

步驟1 輸入離子遷移譜譜圖信號x, 字典Φ，正則化參數(shù)λ1,λ2，差分矩陣d;

步驟2 參數(shù)初始化權重因子W=I, 迭代次數(shù)k=0;

步驟3 求解迭代更新基線、稀疏表示系數(shù)和權重因子:

(1) 基線校正:b(k+1)=(I+λ1DTD)-1(x-Φα(k)+λ1DTDb(k))

(3) 權重因子:W=diag(α(k+1))

步驟4 檢查迭代終止條件: 如果滿足求解精度，則結束迭代; 否則，返回步驟3。

步驟5 輸出基線和稀疏表示求解系數(shù)。

2 實驗部分

2.1 儀器

采用Ni63離子源便攜式離子遷移譜檢測儀。該檢測儀工作模式分為檢測模式和反吹清潔模式，檢測模式下通過流速為470 mL·min-1氣泵進行采樣，氣體通過隔水薄膜后進入遷移管電離區(qū)； 帶電離子在通過采樣周期為20 ms，開門時間為200 μs的離子門后進入電場遷移區(qū)，電場遷移區(qū)電場強度為300 V·cm-1。帶電離子通過法拉第盤采集電流信號，經(jīng)過電流-電壓轉換后通過串口上傳至采集控制軟件，形成橫坐標為時間，縱坐標為電壓幅值的二維離子遷移譜譜圖。

2.2 數(shù)據(jù)

實驗數(shù)據(jù)由仿真數(shù)據(jù)和真實樣本數(shù)據(jù)組成。為了驗證算法，仿真數(shù)據(jù)由三個不同Gaussian分布函數(shù)字典原子、Gaussian白噪聲和基線(分別為正弦波函數(shù)、指數(shù)函數(shù)和線性函數(shù))組成。Gaussian白噪聲標準差為0，幅值為譜圖信號最大幅值的0.5%。真實樣本為水楊酸甲酯(Methylis salicylas, 國藥集團化學試劑有限公司，97%)，水楊酸甲酯因其易于與電子結合且峰位穩(wěn)定常被用作離子遷移譜檢測儀校準物。本實驗通過動態(tài)配氣儀(CPR001, 北京康爾興)將水楊酸甲酯與氮氣在氣袋中混合，樣本濃度為1 mg·m-3。

2.3 系統(tǒng)配置

計算機配置和編程語言如下：

CPU: i7-4790S, 3.20 GHz; 內存: 8 GB。

操作系統(tǒng)： 64bit-Windows 7 Home Edition。

編程語言： MATLAB R2014a (MathWorks, Inc.)。

3 結果與討論

3.1 仿真樣本

仿真數(shù)據(jù)分別采用正弦波函數(shù)、指數(shù)函數(shù)和線性函數(shù)作為基線，采用迭代重加權最小二乘算法對稀疏表示模型求解，迭代求解后基線及基線校正后的譜圖分別如圖1—圖3所示。從圖中可以看出，IRLS算法得到的基線很好的擬合了仿真疊加的三種不同類型基線； Gaussian白噪聲毛刺得到了有效去除，譜圖變得平滑； 經(jīng)過基線校正后的譜圖信號，基線基本在零基線上。采用IRLS算法獲得的仿真譜圖稀疏表示系數(shù)如圖4所示，實驗結果表明，采用稀疏表示模型時通過將噪聲表達系數(shù)舍棄，用最少的字典原子實現(xiàn)了譜圖重建。

圖1 仿真譜圖采用正弦波函數(shù)基線后的結果

圖2 仿真譜圖采用指數(shù)函數(shù)基線后的結果

圖3 仿真譜圖采用線性函數(shù)基線后的結果

為了評估本文提出的可同時實現(xiàn)離子遷移譜譜圖去噪及基線校正的稀疏表示模型及IRLS算法的性能，采用非對稱

圖4 仿真譜圖采用IRLS算法獲得的稀疏表示系數(shù)

最小二乘方法(AsLS)基線校正方法，外加Savitzky-Golay(SG)及小波變換(Wavelet)平滑去噪方法同IRLS算法進行了對比。對采用不同去噪及基線校正算法后的仿真譜圖采用均方根誤差(RMSE)和信噪比(SNR)兩種指標來評估方法性能[計算公式分別為式(4)和式(5)]，計算結果如表1所示。從表1中可以看出IRLS算法在RMSE和SNR指標上都優(yōu)于其他兩種算法，尤其當基線是正弦波基線時，性能有了顯著的提升。

(4)

(5)

表1 采用不同算法后RMSE和SNR指標對比

3.2 真實樣本

水楊酸甲酯在9.84 ms處出現(xiàn)反應離子特征峰，約化遷移率K0為1.62。從采集到的原始水楊酸甲酯離子遷移譜譜圖中可以看出，信號中噪聲波動大、噪聲幅值約為最大幅值的2%； 譜圖基線并不是線性基線。經(jīng)過稀疏表示模型求解后，得到的譜圖基線在零基線上且噪聲得的了很好的去除，實驗結果如圖5所示，不同算法RMSE和SNR指標結果如表2所示。從表2可以看出IRLS算法在RMSE和SNR指標上都顯著優(yōu)于其他兩種算法。

圖5 水楊酸甲酯譜圖稀疏表示模型求解結果

表2 采用不同算法后RMSE和SNR指標對比

4 結 論

應用稀疏表示對離子遷移譜譜圖去噪及基線校正算法進行了研究，提出了一種可同時實現(xiàn)離子遷移譜譜圖去噪和基線校正的稀疏表示模型。仿真和真實樣本實驗表明，采用稀疏表示模型及求解算法得到的實驗結果有了顯著性能提升。離子遷移譜譜圖經(jīng)過去噪和基線校正預處理后，可為下一步化學物質種類的準確識別和定量研究提供參考。