常海濤， 蔡 靜， 溫 悅， 朱玉梅

(北京長城計量測試技術(shù)研究所，北京 100095)

1 引 言

可調(diào)諧半導(dǎo)體激光吸收光譜技術(shù)(tunable diode laser absorption spectroscopy, TDLAS)具有非侵入、高靈敏度和快速時間響應(yīng)等優(yōu)點(diǎn)，廣泛用于航空發(fā)動機(jī)燃?xì)鉁囟仍诰€測量[1～5]。但在實際測量中，由于待測氣體濃度常為毫克每立方米至納克每立方米量級，容易受到探測器噪聲、激光額外噪聲的影響[6～8]。Lines等提出去除噪聲的方法通常從硬件和軟件數(shù)字濾波兩個角度去實施[9]。

硬件去噪方法通常使用高穩(wěn)定性恒流源、液體制冷劑精確控制激光器電流和溫度，通入零氣扣除背景吸收光等方式，但去噪效果受硬件本身性能限制較大，且需要加入元件或增加檢測步驟，增加成本。

常用軟件數(shù)字濾波方法有：算數(shù)平均、滑動平均、非線性最小二乘、小波變換、最大相關(guān)等。Li等經(jīng)分析比較提出軟件濾波不同程度地存在對有用信號影響大、不能同時滿足去除各種噪聲的要求、計算耗時長、實際應(yīng)用時需要更換選擇合適線型、信號噪聲頻率混疊難以分辨等問題[10]。

主成分分析(principle component analysis, PCA)是一種非常有用的統(tǒng)計技術(shù)，它已經(jīng)應(yīng)用于人臉識別和圖像壓縮等領(lǐng)域中，并且是高維數(shù)據(jù)計算模型的常用技術(shù)[11]。本文采用基于PCA的吸收光譜數(shù)據(jù)濾波方法，利用PCA將高維數(shù)據(jù)降成低維數(shù)據(jù)，在保證濾波效果的同時保證了良好的計算效率。

2 原理介紹

2.1 TDLAS雙線測溫原理

TDLAS技術(shù)是建立在Beer-Lambert吸收定律的基礎(chǔ)之上的，當(dāng)一束頻率為ν激光束穿過長度為L的被測氣體時，一部分光強(qiáng)將被吸收，吸收率可以表示為：

(1)

式中：I0和It分別為入射激光強(qiáng)度和透射激光強(qiáng)度；p為總壓力，Pa；T(x)為x處的溫度，K；Xabs(x)為x處的氣體摩爾百分比(氣體濃度)；φ為線型函數(shù)，cm；S(T(x))為躍遷時吸收譜線的線強(qiáng)度，cm-2atm-1，可以表示為：

(2)

式中：h為普朗克常數(shù)，J·s；c為光速，cm/s；k是玻爾茲曼常數(shù)，J/K；Q(T)是配分函數(shù)，它反映了在所處溫度T下，在吸收對應(yīng)低能級上的粒子數(shù)占總粒子數(shù)的比值。Q(T)采用多項式擬合的辦法得到近似值：

Q(T)=a+bT+cT2+dT3

(3)

對于不同的氣體和溫度范圍，系數(shù)a，b，c，d的取值不同，可以在HITRAN光譜數(shù)據(jù)庫中查到。T0是參考溫度，通常為296 K。ν0為吸收譜線中心頻率，cm-1；E″為吸收躍遷的低能級能量，cm-1。可見，吸收譜線的線強(qiáng)度S(T)只是以溫度為變量的函數(shù)。

(4)

本文采用雙線測溫方法，如圖1所示，在掃描波長方案中，在相同壓力、相同摩爾分?jǐn)?shù)和相同路徑長度下同時測量兩個躍遷的積分吸光度(面積)。兩吸收面積之比可消除路徑長度、濃度、壓強(qiáng)和光源波段等對測量的影響，使用雙線直接吸收光譜方法，兩條吸收線的積分吸收率之比R僅僅是溫度的函數(shù)：

(5)

因此，可求得激光路徑上平均溫度的表達(dá)式：

(6)

圖1 雙線測溫原理圖Fig.1 Illustration of TDLAS two-line thermometry

2.2 PCA濾波重建方法

PCA是圖像處理中經(jīng)典的降維和特征提取方法。依據(jù)特征值大小表征光譜變化信息多少原則，對吸收光譜信號和噪聲信號進(jìn)行分離。PCA分析是光譜矩陣，而TDLAS測量得到的吸收光譜是一維矩陣，因此，如何構(gòu)造光譜矩陣是PCA濾波法重要的一步。

本文采用構(gòu)造Hankel矩陣兩種方式對TDLAS原始吸收光譜進(jìn)行處理。令離散的TDLAS吸收光譜數(shù)據(jù)為Aspec={a(1),a(2),…,a(N)}，將其轉(zhuǎn)換成Hankeln×m矩陣A：

(7)

矩陣中每一行代表一個完成周期的吸收光譜信號，橫向代表掃描頻率的變化，縱向代表時間的變化。將光譜矩陣A經(jīng)奇異值分解拆分成標(biāo)準(zhǔn)列正交得分矩陣U、對角陣S和標(biāo)準(zhǔn)行正交載荷矩陣VT乘積的形式：

A=USVT

(8)

對光譜矩陣進(jìn)行特征值分解的本質(zhì)是將原有信號分解到時間和頻率兩個不同的向量子空間，不同奇異值大小代表不同頻率和時間下信號分量投影的大小，特征值越大代表其攜帶光譜變化信息越多。取前n個最大特征值和特征向量作為A的主成分，代表了吸收光譜信號的形狀變化，剩余特征值對應(yīng)的光譜信息均為噪聲。令S*為S矩陣前n個最大特征值組成的對角陣，取前n個對應(yīng)的特征向量的得分矩陣U為U*，載荷矩陣V為V*，則有A的主成分矩陣A0為：

A0=U*S*VT*

(9)

A0為代表吸收光譜主要特征的光譜數(shù)據(jù)，即為濾除噪聲的吸收光譜數(shù)據(jù)。

3 實驗系統(tǒng)搭建

為了驗證本文提出的濾波方法的有效性，搭建TDLAS測溫系統(tǒng)，采用基于主成分濾波方法對實測H2O吸收光譜曲線進(jìn)行濾波處理，并將濾波前后曲線進(jìn)行比較。此外，為了量化濾波效果，本文利用MATLAB軟件模擬理想吸收光譜曲線，附加1倍和5倍隨機(jī)噪聲，并對本文提出的PCA濾波方法與傳統(tǒng)的Savitzky-Golay濾波方法的濾波效果進(jìn)行了比較。

圖2 實驗裝置整體示意圖Fig.2 Schematic of experimental device

整個實驗系統(tǒng)如圖2所示，由光源、分束器、石英氣室、標(biāo)準(zhǔn)具、管式高溫爐、探測器組成。兩個激光器分時復(fù)用，中心波長7 444 cm-1和7 185 cm-1，通過光纖分束器一分為二，一路通過標(biāo)準(zhǔn)具直接由探測器接收，記錄激光波長信息，另一路經(jīng)準(zhǔn)直后通過石英氣室內(nèi)H2O氣體吸收，由探測器獲得吸收光譜。熱源有高溫管式爐提供，將石英管放入管式爐中，同時插入作為溫度測量基準(zhǔn)的熱電偶；將準(zhǔn)直鏡與探測器鏡頭各自對準(zhǔn)石英管的管口。測量過程中石英管兩端通N2避免空氣中H2O吸收干擾。通過電腦控制激光器，通過示波器查看探測器的波形，根據(jù)需要調(diào)整其放大倍數(shù)，令信號在0～2 V內(nèi)。啟動管式爐并設(shè)定溫度(150～1 000 ℃，每隔100 ℃取1點(diǎn))待穩(wěn)定后由高速數(shù)采卡采集吸收光譜數(shù)據(jù)，保存至上位機(jī)。

圖3 實驗裝置實物圖Fig.3 Experimental device

4 實驗驗證與仿真結(jié)果分析

4.1 實驗驗證

基于上一節(jié)搭建的實驗系統(tǒng)，對500 ℃的TDLAS吸收光譜進(jìn)行采集，并采用本文提出的PCA濾波方進(jìn)行濾波處理，對濾波結(jié)果進(jìn)行分析。試驗中激光器采用NEL的DFB激光器，中心波長7 185 cm-1，信號源產(chǎn)生100 Hz鋸齒波信號，輸入給激光器控制器，實現(xiàn)光譜掃描。對采集到吸收光譜信號進(jìn)行預(yù)處理，三次多項式擬合后誤差平方和10-3量級，提高階次誤差減小并不明顯，因此本文采用三次多項式擬合零吸收基線。吸收光譜預(yù)處理如圖4所示，其中圖4(a)給出了5個周期掃描光譜信號，圖4(b)對多周期掃描光譜信號進(jìn)行提取，圖4(c)為三次多項式擬合基線與吸收光譜曲線對比，圖4(d)為去除基線后的TDLAS吸收光譜曲線。

圖4 吸收光譜數(shù)據(jù)預(yù)處理Fig.4 Pretreatment of absorption spectrum data

提取整周期吸收光譜信號，并調(diào)整數(shù)據(jù)長度保持一致，取10組吸收光譜信號，構(gòu)造Hankel矩陣，每一行代表1個完成周期的吸收光譜信號。對光譜進(jìn)行進(jìn)行PCA濾波重建，圖5給出了不同主成分?jǐn)?shù)濾波重建后曲線與原始吸收光譜曲線對比，圖5(a)～圖5(d)分別對應(yīng)1～4個主成分?jǐn)?shù)。如圖所示，主成分?jǐn)?shù)為1時，濾波效果最好，隨著主成分?jǐn)?shù)的增加，冗余的噪聲信息逐漸滲入吸收光譜信號中，濾波后曲線明顯攜帶了高頻噪聲，并且隨著主成分?jǐn)?shù)的增加高頻噪聲攜帶量逐漸增大。

圖5 不同主成分?jǐn)?shù)濾波效果比較Fig.5 Comparison of filtering effects of different principal components

圖6給出了不同主成分的貢獻(xiàn)值，其中第1主成分貢獻(xiàn)值最大，表示其承載的光譜曲線特征變化信息最多，從第2個組成貢獻(xiàn)值開始，后面貢獻(xiàn)值基本相等，接近于0，表明其表征信息的基本上均為高頻噪聲信息。這與圖5中給出的1個主成分進(jìn)行重建后的光譜曲線濾波效果最好結(jié)論一致。

圖6 各主成分貢獻(xiàn)值分布Fig.6 The magnitude of the eigenvalues in decreasing order

4.2 仿真驗證及濾波效果評估

為了進(jìn)一步量化濾波效果，本文用MATLAB軟件模擬理想吸收光譜曲線，并附加1倍和5倍隨機(jī)噪聲，形成噪聲水平不同的吸收光譜矩陣。分別對1倍和5倍隨機(jī)噪聲吸收光譜曲線進(jìn)行PCA濾波處理，結(jié)果如圖7～圖8所示。

圖7(a)為加載1倍隨機(jī)噪聲前后吸收光譜曲線；圖7(b)給出了特征值隨主成分?jǐn)?shù)變化，從第4個主成分開始，特征值接近于0，僅表征噪聲信息；圖7(c)是濾波重建曲線與理想曲線比較，濾波后曲線光滑平整，且與理想曲線的重合度較高，說明濾波效果良好。由圖7(d)可知濾波后曲線與理想曲線偏差值，基本在10-3V量級，相對偏差0.1%左右。

圖7 1倍隨機(jī)噪聲吸收光譜PCA濾波效果比較Fig.7 The filtering process of PCA algorithm to contaminated curve with 1-times noise

圖8給出了5倍隨機(jī)噪聲吸收光譜PCA濾波效果比較，噪聲明顯更大，第3個主成分開始特征接近為0，重建曲線與理想曲線相對偏差相對偏差0.2%左右，濾波效果良好。

圖8 5倍隨機(jī)噪聲吸收光譜PCA濾波效果比較Fig.8 The filtering process of PCA algorithm to contaminated curve with 5-times noise

最后將本文提出的PCA濾波方法與傳統(tǒng)的Savitzky-Golay濾波方法進(jìn)行比較，用理想吸收光譜曲線與濾波后曲線的誤差平方和(SSR)表征濾波效果，表1給出了2種濾波方法濾波結(jié)果比較。由表1可知，對于1倍隨機(jī)噪聲吸收光譜信號和5倍噪聲吸收光譜信號，PCA濾波后SSR值分別從2.21和55.252降低到了0.002和0.018，噪聲去除率達(dá)81%，且優(yōu)于Savitzky-Golay濾波的0.006和0.006 2，充分證明本文所提出的PCA濾波方法對于不同噪聲水平的吸收光譜曲線高頻噪聲濾除明顯，且噪聲去除率優(yōu)于傳統(tǒng)Savitzky-Golay濾波方法。

表1 不同濾波方法濾波效果比較Tab.1 SSR of different filters

采用圖3所示實驗裝置，對溫度范圍300～900 ℃，步長100 ℃，7個溫度點(diǎn)H2O氣體吸收光譜進(jìn)行測量，采用直接吸收雙線測溫法進(jìn)行溫度解算，以熱電偶溫度傳感器讀數(shù)值為真值。溫度解算過程中分別采用主成分濾波和Savitzky-Golay濾波兩種方式進(jìn)行光譜數(shù)據(jù)的預(yù)處理，最終溫度解算結(jié)果如表2所示。主成分濾波后解算出的溫度值與熱電偶讀數(shù)之間偏差明顯更小，標(biāo)準(zhǔn)差從8.9降到5.3，進(jìn)一步證明本文所提出的濾波方法濾波效果優(yōu)于傳統(tǒng)方法。

表2 基于不同濾波方法氣體溫度解算結(jié)果比較Tab.2 Comparison of gas temperature calculated by different filtering methods ℃

5 結(jié) 論

本文介紹了一種基于主成分分析的TDLAS吸收光譜高頻噪聲濾波技術(shù)，對TDALS雙線測溫原理、PCA濾波重建方法的原理進(jìn)行了詳細(xì)闡述。搭建了TDLAS溫度測量裝置，對實測吸收光譜曲線進(jìn)行PCA濾波，實驗結(jié)果表明，該方法濾波效果明顯。為進(jìn)一步對濾波效果進(jìn)行量化和評估，對理想高斯吸收譜線疊加1倍隨機(jī)噪聲和5倍隨機(jī)噪聲，形成兩種不同噪聲水平的吸收光譜信號，用本文提出的PCA濾波方法和傳統(tǒng)的Savitzky-Golay濾波進(jìn)行高頻噪聲濾除，并通過理想吸收光譜曲線與濾波后曲線的誤差平方和表征濾波效果進(jìn)行比較。從試驗結(jié)果可以看出，PCA濾波后曲線與理想曲線重合度較高，誤差平方和分別降低到了0.002和0.018，噪聲去除率達(dá)81%，且優(yōu)于Savitzky-Golay濾波的0.006和0.006 2。