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

中國工程物理研究院表面物理與化學重點實驗室，四川 綿陽 621908

引 言

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

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

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

1 壓縮感知重構與算法

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

1.1 信號的稀疏表示

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

x=Ψs

(1)

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

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

1.2 信號測量

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

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

(2)

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

1.3 信號的重構及重構算法

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

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

(3)

本研究選取正交匹配追蹤(orthogonal matching pursuit,OMP)[10-11]算法對氫同位素氣體拉曼光譜信號進行重構，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內積中的最大值，即找索引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(氘)組成的氫氣作為研究對象，基于壓縮感知理論，①對拉曼光譜傅里葉變換，帶阻濾波后反變換回到時域，并結合多組正交基構建正交基字典，用OMP算法(下文簡稱FTOMP)對拉曼光譜進行重構；②對拉曼光譜進行尋峰和多項式擬合基線[12]，然后通過洛倫茲函數(shù)設計原子，構造字典，用OMP算法(下文簡稱LoOMP)對拉曼光譜進行重構。并與小波軟閾值、小波硬閾值和SG濾波處理結果比對，分析幾種方法處理后的信噪比和均方根誤差。

2 實驗部分

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

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

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

圖1 (a)無噪聲的仿真光譜；(b)加入噪聲的仿真光譜；(c)FTOMP重構的光譜；(d)LoOMP重構的光譜；(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高斯噪聲。然后使用幾種方法對含噪聲譜圖進行處理，圖1(b)為加入5 db高斯白噪聲后的仿真拉曼光譜，圖1(c)為FTOMP算法的重構光譜,圖1(d)為LoOMP算法的重構光譜，圖1(e)為小波軟閾值處理得到的降噪光譜，圖1(f)為小波硬閾值處理得到的降噪光譜，圖1(g)為經(jīng)SG濾波獲得的降噪光譜。其中使用LoOMP算法重構的拉曼光譜的參數(shù)信息如表2所示。

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

并通過計算降噪后的信噪比(SNR)和均方根誤差(RMSE)以及譜峰強度偏差(Intensity error)這三個指標對降噪效果進行評價，表達式分別如式(4)—式(6)

(4)

(5)

(6)

圖2 壓縮感知方法重構信號(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ù)據(jù)處理

實驗采用自制拉曼光譜儀對氫氣同位素氣體進行測量，激光波長為532 nm。樣品池中充入氫氣和氘氣的混合氣體，光譜采集時間分別為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支振動信號。

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

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

圖4 (a)0.1 s采集的H2和D2混合氣體的拉曼光譜；(b)FTOMP重構的光譜；(c)LoOMP重構的光譜； (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信號的譜峰強度偏差Table 4 The difference in the D2 signal intensity by the five methods

圖5 兩種壓縮感知方法處理后光譜中D2信號的強度偏差(a)：信號強度偏差與采樣率M/N的關系曲線；(b)：信號強度偏差與稀疏度K的關系曲線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

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

圖6 五種方法處理的拉曼光譜中D2信號(a)：譜峰強度與采集時間的關系曲線；(b)：譜峰強度偏差與采集時間的關系曲線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 結果與討論

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

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

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

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

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

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

4 結 論

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