任永甜，胡 儀，陳 駿，陳 鈞

中國工程物理研究院表面物理與化學(xué)重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621908

引 言

氣體分析應(yīng)用廣泛，關(guān)系到國計(jì)民生。無論是在科學(xué)研究中還是生產(chǎn)實(shí)踐中，氣體分析都有重要作用。氣體成分分析的手段有氣相色譜法、質(zhì)譜法、紅外光譜法、吸收光譜法和激光拉曼光譜法等。在光譜技術(shù)中，激光拉曼光譜具有非接觸、無損等優(yōu)點(diǎn)[1]，但是，拉曼散射強(qiáng)度較低，實(shí)際測得的拉曼信號會受到射線、熒光、電荷耦合元件暗電流以及光學(xué)系統(tǒng)雜散光等干擾。為了提高信號質(zhì)量，除了增加曝光時(shí)間或者提升硬件性能外，還可以利用信號處理技術(shù)對信號進(jìn)行降噪處理。

目前用于光譜平滑降噪處理的方法有小波變換(wavelet transform)[2]、經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)[3]、Sawitzky-Golay(SG)濾波等[4]。但是在一些實(shí)際應(yīng)用場景中，為了觀察動(dòng)態(tài)過程，信號測量時(shí)間較短，獲得的信號的信噪比也比較低，而目前的一些降噪算法處理微弱信號的效果不是很好。除了小波去噪，常用于微弱信號檢測的方法有盲源分離、相關(guān)檢測、混沌檢測和壓縮感知(compressed sensing,CS)等。近年來，研究人員將壓縮感知用于光譜的重建，取得不錯(cuò)的效果[5-6]。

氫氣與我們的生活息息相關(guān)，氫氣分析測量在國際熱核聚變實(shí)驗(yàn)堆(international thermonuclear experimental reactor,ITER)、氫能等領(lǐng)域有重要意義。此外，氫氣分子結(jié)構(gòu)簡單，是理想的研究模型。因此，本研究選取氫氣和氘氣作為測量對象，獲取在采集時(shí)間較短情況下低信噪比的氫同位素氣體拉曼光譜，研究壓縮感知理論用于解析氫氣拉曼光譜的可行性，并與其他降噪方法比對，結(jié)果表明壓縮感知理論為氫氣拉曼光譜降噪提供了一種有效的途徑。

1 壓縮感知重構(gòu)與算法

2006年，Donoho，Candès，Romberg和Tao[7-9]等在研究信號稀疏性的時(shí)候發(fā)現(xiàn)只要一個(gè)信號是稀疏的或可壓縮的，就可以通過合適的基或者框架對高維信號進(jìn)行稀疏表示，然后使用高效算法從高度不完整的線性測量中恢復(fù)原始信號。可以通過更少的測量來感知稀疏信號，因此得名壓縮感知。壓縮感知不僅可以壓縮采樣，縮短采樣時(shí)間，而且可以降噪，提高信噪比，以更好地實(shí)現(xiàn)原始信號的恢復(fù)和重建。壓縮感知一般是由三個(gè)部分組成：信號的稀疏表示，信號測量和信號重構(gòu)。

1.1 信號的稀疏表示

一個(gè)有限長一維離散信號x是RN空間N×1維的列向量，能夠由N×N的稀疏矩陣Ψ表示為

x=Ψs

(1)

式(1)中，N×1的列向量s為信號x的稀疏系數(shù)，x如果僅有少數(shù)非零項(xiàng)，Ψ則是單位矩陣。除了傅里葉變換、小波變換等基函數(shù)，還可以用超完備的冗余字典稀疏表示數(shù)據(jù)。

本研究中兩種壓縮感知分析分別使用多組正交基構(gòu)建的正交基字典和根據(jù)氫同位素氣體拉曼光譜特征設(shè)計(jì)的字典。根據(jù)拉曼光譜中峰的強(qiáng)度和面積，設(shè)計(jì)洛倫茲函數(shù)線性組合F(x)，F(xiàn)(x)作為純凈無噪聲的拉曼信號，‖s‖0=K，由F(x)=Ψs可以得到稀疏矩陣Ψ。其中K表示稀疏度。

1.2 信號測量

通過M×N維的測量矩陣Φ與信號x做內(nèi)積，得到測量值y，因?yàn)镸遠(yuǎn)遠(yuǎn)小于N，所以可以減少測量數(shù)據(jù)量，達(dá)到降維的目的。測量值y可以寫成

y=Φx+e=ΦΨs+e=Θs+e

(2)

式(2)中，y為M×1的列向量，Θ=ΦΨ為M×N的傳感矩陣，e表示有界能量的噪聲項(xiàng)[8]，‖e‖2≤ε。

1.3 信號的重構(gòu)及重構(gòu)算法

信號的重構(gòu)可以轉(zhuǎn)換為l1最小范數(shù)下的最優(yōu)化問題

min‖x‖1s.t.‖y-Φx‖2≤ε

(3)

本研究選取正交匹配追蹤(orthogonal matching pursuit,OMP)[10-11]算法對氫同位素氣體拉曼光譜信號進(jìn)行重構(gòu)，OMP算法步驟如下：

輸入：M×N維的傳感矩陣Θ=ΦΨ，M×1的測量值y，信號稀疏度K。

在算法中，用?表示空集，Λt表示t次迭代的索引的集合，Jt表示t次迭代的索引號。

初始化：殘差r0=y，索引集Λt=?，Θ0=?，迭代次數(shù)t=1。

循環(huán)執(zhí)行步驟(1)—步驟(5)：

(1)找出殘差rt-1和傳感矩陣Θj內(nèi)積中的最大值，即找索引Jt，使得Jt=argmaxj=1，2，…，n|〈rt-1，Θj〉|，其中，Θj是傳感矩陣第j列；

(2)更新索引集Λt=Λt-1∪{Jt}，記錄找到的傳感矩陣中的重建原子集合Dt=[Dt-1，ΘJt]；

(5)令t=t+1，若t

本研究選用H2(氕)和D2(氘)組成的氫氣作為研究對象，基于壓縮感知理論，①對拉曼光譜傅里葉變換，帶阻濾波后反變換回到時(shí)域，并結(jié)合多組正交基構(gòu)建正交基字典，用OMP算法(下文簡稱FTOMP)對拉曼光譜進(jìn)行重構(gòu)；②對拉曼光譜進(jìn)行尋峰和多項(xiàng)式擬合基線[12]，然后通過洛倫茲函數(shù)設(shè)計(jì)原子，構(gòu)造字典，用OMP算法(下文簡稱LoOMP)對拉曼光譜進(jìn)行重構(gòu)。并與小波軟閾值、小波硬閾值和SG濾波處理結(jié)果比對，分析幾種方法處理后的信噪比和均方根誤差。

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

2.1 仿真數(shù)據(jù)處理

首先進(jìn)行了仿真實(shí)驗(yàn)，使用洛倫茲函數(shù)來仿真純凈的拉曼光譜信號。為了仿真弱拉曼光譜信號，仿真信號設(shè)計(jì)的強(qiáng)度較低，每個(gè)峰的信號強(qiáng)度均不超過20。仿真信號峰的位置、峰的強(qiáng)度和峰的半高寬等相關(guān)參數(shù)如表1所示。圖1(a)為仿真的無噪聲拉曼光譜。

表1 仿真信號參數(shù)表Table 1 Parameters of only simulated Raman signal

圖1 (a)無噪聲的仿真光譜；(b)加入噪聲的仿真光譜；(c)FTOMP重構(gòu)的光譜；(d)LoOMP重構(gòu)的光譜；(e)小波軟閾值降噪的光譜；(f)小波硬閾值降噪的光譜；(g)經(jīng)SG濾波的光譜Fig.1 (a)Simulated spectrum;(b)Simulated spectrum with noise;(c)Reconstructed spectrum by FTOMP;(d)Reconstructed spectrum by LoOMP;(e)Denoising spectrum with Wavelet soft threshold;(f)Denoising spectrum with Wavelet hard threshold;(g)SG filtered spectrum

為了測試壓縮感知算法的降噪能力，對圖1(a)中無噪聲的仿真光譜信號加入5 db高斯噪聲。然后使用幾種方法對含噪聲譜圖進(jìn)行處理，圖1(b)為加入5 db高斯白噪聲后的仿真拉曼光譜，圖1(c)為FTOMP算法的重構(gòu)光譜,圖1(d)為LoOMP算法的重構(gòu)光譜，圖1(e)為小波軟閾值處理得到的降噪光譜，圖1(f)為小波硬閾值處理得到的降噪光譜，圖1(g)為經(jīng)SG濾波獲得的降噪光譜。其中使用LoOMP算法重構(gòu)的拉曼光譜的參數(shù)信息如表2所示。

表2 LoOMP重構(gòu)信號的參數(shù)Table 2 Parameters of reconstructed Raman signals by LoOMP

并通過計(jì)算降噪后的信噪比(SNR)和均方根誤差(RMSE)以及譜峰強(qiáng)度偏差(Intensity error)這三個(gè)指標(biāo)對降噪效果進(jìn)行評價(jià)，表達(dá)式分別如式(4)—式(6)

(4)

(5)

(6)

圖2 壓縮感知方法重構(gòu)信號(a)：M/N-SNR曲線；(b)：M/N-RMSE曲線；(c)：K-SNR曲線；(d)：K-RMSE曲線Fig.2 Compressed sensing method to reconstruct the signal(a):M/N-SNR curves;(b):M/N-RMSE curves;(c):K-SNR curves;(d):K-RMSE curves

表3 含噪聲的仿真光譜處理前后的SNR和RMSETable 3 SNR and RMSE before and after processing the treatments of the simulated noisy Raman spectrum

2.2 實(shí)際測量數(shù)據(jù)處理

實(shí)驗(yàn)采用自制拉曼光譜儀對氫氣同位素氣體進(jìn)行測量，激光波長為532 nm。樣品池中充入氫氣和氘氣的混合氣體，光譜采集時(shí)間分別為1，0.15，0.10，0.08，0.06，0.05和0.04 s。其中1 s采集的拉曼光譜如圖3所示。圖中4 163和2 993 cm-1處的譜峰分別為H2和D2的Q1支振動(dòng)信號。

圖3 1 s采集的H2和D2混合氣體的拉曼光譜Fig.3 Raman spectrum of H2 and D2 mixed gas with 1 s collection time

以0.1 s拉曼光譜作為處理對象，并對其進(jìn)行歸一化處理，比較五種方法降噪處理后的譜圖，如圖4所示，圖4(a)是0.1 s采集的H2和D2混合氣體的拉曼光譜，圖4(b)是FTOMP算法重構(gòu)的光譜，圖4(c)是LoOMP算法重構(gòu)的光譜，圖4(d)是小波軟閾值降噪后的光譜，圖4(e)是小波硬閾值降噪后的光譜，圖4(f)是經(jīng)SG濾波降噪后的光譜。從譜圖上可以看出，LoOMP重構(gòu)處理的譜圖曲線最光滑。五種方法處理后，拉曼光譜中D2信號的譜峰強(qiáng)度與同一采集時(shí)間下連續(xù)采集100次拉曼光譜得到的平均譜作比較，得到D2信號強(qiáng)度偏差如表4所示。然后計(jì)算了FTOMP重構(gòu)后的光譜與LoOMP重構(gòu)后的光譜在不同采樣率M/N和稀疏度K下D2信號的譜峰強(qiáng)度偏差，結(jié)果如圖5所示，M表示測量值，N表示信號的長度。

圖4 (a)0.1 s采集的H2和D2混合氣體的拉曼光譜；(b)FTOMP重構(gòu)的光譜；(c)LoOMP重構(gòu)的光譜； (d)小波軟閾值降噪的光譜；(e)小波硬閾值降噪的光譜；(f)經(jīng)SG濾波的光譜Fig.4 (a)Raman spectrum of H2 and D2 mixed gas with 1 s collection time;(b)Reconstructed spectrum by FTOMP;(c)Reconstructed spectrum by LoOMP;(d)Denoising spectrum by Wavelet soft threshold;(e)Denoising spectrum by Wavelet hard threshold;(f)SG filtered spectrum

表4 五種方法處理后拉曼光譜中D2信號的譜峰強(qiáng)度偏差Table 4 The difference in the D2 signal intensity by the five methods

圖5 兩種壓縮感知方法處理后光譜中D2信號的強(qiáng)度偏差(a)：信號強(qiáng)度偏差與采樣率M/N的關(guān)系曲線；(b)：信號強(qiáng)度偏差與稀疏度K的關(guān)系曲線Fig.5 The difference in the D2 signal intensity by LoOMP and FTOMP(a)：Intensity error of D2 signal vs M/N value;(b)：Intensity error of D2 signal vs K value

六組采集時(shí)間下的拉曼光譜經(jīng)過幾種方法處理得到的譜圖在D2信號處的強(qiáng)度如圖6(a)所示，圖6(b)是五種方法處理后D2信號處的強(qiáng)度偏差。

圖6 五種方法處理的拉曼光譜中D2信號(a)：譜峰強(qiáng)度與采集時(shí)間的關(guān)系曲線；(b)：譜峰強(qiáng)度偏差與采集時(shí)間的關(guān)系曲線Fig.6 The difference in the D2 signal intensity by LoOMP and FTOMP(a)：Intensity of D2 signal vs collection time;(b)：Intensity error of D2 signal vs collection time

3 結(jié)果與討論

從表1—表3以及圖1可以看出，幾種方法處理后光譜的信噪比均有所提高，均方根誤差減小。本研究提出的LoOMP算法能夠很好地重構(gòu)拉曼光譜的拉曼峰，獲得的譜圖的信噪比最大，均方根誤差最小，其次是小波軟閾值處理后的光譜，F(xiàn)TOMP算法重構(gòu)光譜與SG濾波降噪光譜的信噪比最低，均方根誤差最大。當(dāng)M/N=0.8，K=50時(shí)，F(xiàn)TOMP算法重構(gòu)譜圖的信噪比為10.274 8，均方根誤差為0.877 1，比LoOMP重構(gòu)效果差。通過仿真實(shí)驗(yàn)可知LoOMP算法處理拉曼光譜的降噪效果相對最優(yōu)。

圖2探究了采樣率和稀疏度對壓縮感知算法重構(gòu)信號效果的影響。隨著采樣率的增大，LoOMP算法重構(gòu)光譜的SNR和RMSE變化較小；當(dāng)M/N達(dá)到0.4～0.5時(shí)，F(xiàn)TOMP的重構(gòu)效果達(dá)到最佳，SNR和RMSE隨M/N變化趨于穩(wěn)定；隨著稀疏度的增大，兩種壓縮感知算法重構(gòu)光譜的SNR和RMSE增大，當(dāng)K=50時(shí)，此時(shí)SNR和RMSE最佳，并且隨著K的增加趨于穩(wěn)定。比較LoOMP算法和FTOMP算法重構(gòu)光譜的SNR和RMSE，LoOMP算法重構(gòu)效果要好。

結(jié)合圖4和表4，可以看出，實(shí)驗(yàn)數(shù)據(jù)和仿真結(jié)果一致，LoOMP算法在拉曼峰處重構(gòu)效果最佳，譜峰強(qiáng)度偏差最小，降噪效果最好；小波軟閾值降噪效果較好，但是在拉曼峰處強(qiáng)度偏差較大；小波硬閾值和SG濾波處理結(jié)果在拉曼峰處強(qiáng)度偏差居中，但是降噪效果較差；FTOMP處理效果最差。

圖5展示了隨著采樣率和稀疏度的變化，拉曼峰強(qiáng)度的重構(gòu)偏差結(jié)果。隨著采樣率的增大，強(qiáng)度偏差變化趨勢較小，LoOMP的強(qiáng)度偏差要比FTOMP小得多；隨著稀疏度的增大，強(qiáng)度偏差在減小，LoOMP算法重構(gòu)光譜的譜峰強(qiáng)度偏差比FTOMP算法小很多，可以得出LoOMP算法相對FTOMP算法重構(gòu)的信號效果要好。圖6是不同采集時(shí)間下的拉曼光譜經(jīng)過幾種方法重構(gòu)后D2信號的強(qiáng)度的變化以及偏差?？梢钥闯?，LoOMP重構(gòu)譜峰強(qiáng)度和原始光譜譜峰強(qiáng)度偏差最小，且相對穩(wěn)定；小波硬閾值和SG濾波兩種方法在譜峰處重構(gòu)效果較好；小波軟閾值和FTOMP在譜峰處重構(gòu)強(qiáng)度偏差較大。但是在采集時(shí)間較長情況下小波軟閾值重構(gòu)效果好于FTOMP。

LoOMP和FTOMP兩種壓縮感知算法的處理結(jié)果差別較大，主要原因是LoOMP根據(jù)拉曼光譜特征，使用洛倫茲函數(shù)構(gòu)建字典的原子，這樣可以更好地保留光譜特性，該方法適用于拉曼光譜的降噪，可以恢復(fù)拉曼光譜的特定主峰，其效率在某些方面優(yōu)于其他方法，對評價(jià)拉曼光譜和檢測被測物具有重要意義，在安檢等一些需要定性分析場景中具有一定的應(yīng)用價(jià)值；FTOMP使用正交基構(gòu)建的字典對拉曼光譜稀疏表示效果較差，而且經(jīng)過傅里葉濾波處理，譜峰強(qiáng)度會有所降低。

將壓縮感知用于拉曼光譜重構(gòu)的工作，重構(gòu)拉曼光譜信號峰的頻率、峰寬和強(qiáng)度[13]，可以掃描幾個(gè)光譜基[14]，也可以減少信號采集時(shí)間[5]。本研究與現(xiàn)有的降噪方法相比，主要是將壓縮感知算法應(yīng)用于氫氣拉曼光譜降噪處理，結(jié)合氫氣拉曼光譜信號特征，使用洛倫茲函數(shù)設(shè)計(jì)稀疏矩陣的原子，可以更好地稀疏表示原始拉曼光譜數(shù)據(jù)，以達(dá)到在重構(gòu)過程中對氫氣拉曼光譜降噪的目的。而且LoOMP算法處理氫氣拉曼光譜的時(shí)候不需要對拉曼光譜進(jìn)行分段處理，可實(shí)現(xiàn)直接對幾個(gè)拉曼峰同時(shí)重構(gòu)，降噪的同時(shí)還保留了拉曼信號的強(qiáng)度和譜峰基本形狀信息。

4 結(jié) 論

運(yùn)用壓縮感知算法分析氫同位素氣體的拉曼光譜，研究兩種不同的稀疏矩陣構(gòu)建方式對重構(gòu)結(jié)果的影響。通過對仿真光譜和實(shí)際測量拉曼光譜進(jìn)行處理比對，研究表明：提出的LoOMP壓縮感知算法可以有效地降低噪聲，同時(shí)信號強(qiáng)度偏差很小。通過比對，LoOMP壓縮感知算法重構(gòu)的氫同位素氣體拉曼光譜的信噪比和均方根誤差優(yōu)于FTOMP壓縮感知算法、小波軟閾值法、小波硬閾值法和SG濾波方法，且拉曼峰的強(qiáng)度重構(gòu)偏差最小。鑒于氣體拉曼信號特征相似，本研究的LoOMP壓縮感知算法也適用于其他氣體的拉曼光譜的降噪處理。