常海涛， 蔡 静， 温 悦， 朱玉梅

(北京长城计量测试技术研究所，北京 100095)

1 引 言

可调谐半导体激光吸收光谱技术(tunable diode laser absorption spectroscopy, TDLAS)具有非侵入、高灵敏度和快速时间响应等优点，广泛用于航空发动机燃气温度在线测量[1～5]。但在实际测量中，由于待测气体浓度常为毫克每立方米至纳克每立方米量级，容易受到探测器噪声、激光额外噪声的影响[6～8]。Lines等提出去除噪声的方法通常从硬件和软件数字滤波两个角度去实施[9]。

硬件去噪方法通常使用高稳定性恒流源、液体制冷剂精确控制激光器电流和温度，通入零气扣除背景吸收光等方式，但去噪效果受硬件本身性能限制较大，且需要加入元件或增加检测步骤，增加成本。

常用软件数字滤波方法有：算数平均、滑动平均、非线性最小二乘、小波变换、最大相关等。Li等经分析比较提出软件滤波不同程度地存在对有用信号影响大、不能同时满足去除各种噪声的要求、计算耗时长、实际应用时需要更换选择合适线型、信号噪声频率混叠难以分辨等问题[10]。

主成分分析(principle component analysis, PCA)是一种非常有用的统计技术，它已经应用于人脸识别和图像压缩等领域中，并且是高维数据计算模型的常用技术[11]。本文采用基于PCA的吸收光谱数据滤波方法，利用PCA将高维数据降成低维数据，在保证滤波效果的同时保证了良好的计算效率。

2 原理介绍

2.1 TDLAS双线测温原理

TDLAS技术是建立在Beer-Lambert吸收定律的基础之上的，当一束频率为ν激光束穿过长度为L的被测气体时，一部分光强将被吸收，吸收率可以表示为：

(1)

式中：I0和It分别为入射激光强度和透射激光强度；p为总压力，Pa；T(x)为x处的温度，K；Xabs(x)为x处的气体摩尔百分比(气体浓度)；φ为线型函数，cm；S(T(x))为跃迁时吸收谱线的线强度，cm-2atm-1，可以表示为：

(2)

式中：h为普朗克常数，J·s；c为光速，cm/s；k是玻尔兹曼常数，J/K；Q(T)是配分函数，它反映了在所处温度T下，在吸收对应低能级上的粒子数占总粒子数的比值。Q(T)采用多项式拟合的办法得到近似值：

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

(3)

对于不同的气体和温度范围，系数a，b，c，d的取值不同，可以在HITRAN光谱数据库中查到。T0是参考温度，通常为296 K。ν0为吸收谱线中心频率，cm-1；E″为吸收跃迁的低能级能量，cm-1。可见，吸收谱线的线强度S(T)只是以温度为变量的函数。

(4)

本文采用双线测温方法，如图1所示，在扫描波长方案中，在相同压力、相同摩尔分数和相同路径长度下同时测量两个跃迁的积分吸光度(面积)。两吸收面积之比可消除路径长度、浓度、压强和光源波段等对测量的影响，使用双线直接吸收光谱方法，两条吸收线的积分吸收率之比R仅仅是温度的函数：

(5)

因此，可求得激光路径上平均温度的表达式：

(6)

图1 双线测温原理图Fig.1 Illustration of TDLAS two-line thermometry

2.2 PCA滤波重建方法

PCA是图像处理中经典的降维和特征提取方法。依据特征值大小表征光谱变化信息多少原则，对吸收光谱信号和噪声信号进行分离。PCA分析是光谱矩阵，而TDLAS测量得到的吸收光谱是一维矩阵，因此，如何构造光谱矩阵是PCA滤波法重要的一步。

本文采用构造Hankel矩阵两种方式对TDLAS原始吸收光谱进行处理。令离散的TDLAS吸收光谱数据为Aspec={a(1),a(2),…,a(N)}，将其转换成Hankeln×m矩阵A：

(7)

矩阵中每一行代表一个完成周期的吸收光谱信号，横向代表扫描频率的变化，纵向代表时间的变化。将光谱矩阵A经奇异值分解拆分成标准列正交得分矩阵U、对角阵S和标准行正交载荷矩阵VT乘积的形式：

A=USVT

(8)

对光谱矩阵进行特征值分解的本质是将原有信号分解到时间和频率两个不同的向量子空间，不同奇异值大小代表不同频率和时间下信号分量投影的大小，特征值越大代表其携带光谱变化信息越多。取前n个最大特征值和特征向量作为A的主成分，代表了吸收光谱信号的形状变化，剩余特征值对应的光谱信息均为噪声。令S*为S矩阵前n个最大特征值组成的对角阵，取前n个对应的特征向量的得分矩阵U为U*，载荷矩阵V为V*，则有A的主成分矩阵A0为：

A0=U*S*VT*

(9)

A0为代表吸收光谱主要特征的光谱数据，即为滤除噪声的吸收光谱数据。

3 实验系统搭建

为了验证本文提出的滤波方法的有效性，搭建TDLAS测温系统，采用基于主成分滤波方法对实测H2O吸收光谱曲线进行滤波处理，并将滤波前后曲线进行比较。此外，为了量化滤波效果，本文利用MATLAB软件模拟理想吸收光谱曲线，附加1倍和5倍随机噪声，并对本文提出的PCA滤波方法与传统的Savitzky-Golay滤波方法的滤波效果进行了比较。

图2 实验装置整体示意图Fig.2 Schematic of experimental device

整个实验系统如图2所示，由光源、分束器、石英气室、标准具、管式高温炉、探测器组成。两个激光器分时复用，中心波长7 444 cm-1和7 185 cm-1，通过光纤分束器一分为二，一路通过标准具直接由探测器接收，记录激光波长信息，另一路经准直后通过石英气室内H2O气体吸收，由探测器获得吸收光谱。热源有高温管式炉提供，将石英管放入管式炉中，同时插入作为温度测量基准的热电偶；将准直镜与探测器镜头各自对准石英管的管口。测量过程中石英管两端通N2避免空气中H2O吸收干扰。通过电脑控制激光器，通过示波器查看探测器的波形，根据需要调整其放大倍数，令信号在0～2 V内。启动管式炉并设定温度(150～1 000 ℃，每隔100 ℃取1点)待稳定后由高速数采卡采集吸收光谱数据，保存至上位机。

图3 实验装置实物图Fig.3 Experimental device

4 实验验证与仿真结果分析

4.1 实验验证

基于上一节搭建的实验系统，对500 ℃的TDLAS吸收光谱进行采集，并采用本文提出的PCA滤波方进行滤波处理，对滤波结果进行分析。试验中激光器采用NEL的DFB激光器，中心波长7 185 cm-1，信号源产生100 Hz锯齿波信号，输入给激光器控制器，实现光谱扫描。对采集到吸收光谱信号进行预处理，三次多项式拟合后误差平方和10-3量级，提高阶次误差减小并不明显，因此本文采用三次多项式拟合零吸收基线。吸收光谱预处理如图4所示，其中图4(a)给出了5个周期扫描光谱信号，图4(b)对多周期扫描光谱信号进行提取，图4(c)为三次多项式拟合基线与吸收光谱曲线对比，图4(d)为去除基线后的TDLAS吸收光谱曲线。

图4 吸收光谱数据预处理Fig.4 Pretreatment of absorption spectrum data

提取整周期吸收光谱信号，并调整数据长度保持一致，取10组吸收光谱信号，构造Hankel矩阵，每一行代表1个完成周期的吸收光谱信号。对光谱进行进行PCA滤波重建，图5给出了不同主成分数滤波重建后曲线与原始吸收光谱曲线对比，图5(a)～图5(d)分别对应1～4个主成分数。如图所示，主成分数为1时，滤波效果最好，随着主成分数的增加，冗余的噪声信息逐渐渗入吸收光谱信号中，滤波后曲线明显携带了高频噪声，并且随着主成分数的增加高频噪声携带量逐渐增大。

图5 不同主成分数滤波效果比较Fig.5 Comparison of filtering effects of different principal components

图6给出了不同主成分的贡献值，其中第1主成分贡献值最大，表示其承载的光谱曲线特征变化信息最多，从第2个组成贡献值开始，后面贡献值基本相等，接近于0，表明其表征信息的基本上均为高频噪声信息。这与图5中给出的1个主成分进行重建后的光谱曲线滤波效果最好结论一致。

图6 各主成分贡献值分布Fig.6 The magnitude of the eigenvalues in decreasing order

4.2 仿真验证及滤波效果评估

为了进一步量化滤波效果，本文用MATLAB软件模拟理想吸收光谱曲线，并附加1倍和5倍随机噪声，形成噪声水平不同的吸收光谱矩阵。分别对1倍和5倍随机噪声吸收光谱曲线进行PCA滤波处理，结果如图7～图8所示。

图7(a)为加载1倍随机噪声前后吸收光谱曲线；图7(b)给出了特征值随主成分数变化，从第4个主成分开始，特征值接近于0，仅表征噪声信息；图7(c)是滤波重建曲线与理想曲线比较，滤波后曲线光滑平整，且与理想曲线的重合度较高，说明滤波效果良好。由图7(d)可知滤波后曲线与理想曲线偏差值，基本在10-3V量级，相对偏差0.1%左右。

图7 1倍随机噪声吸收光谱PCA滤波效果比较Fig.7 The filtering process of PCA algorithm to contaminated curve with 1-times noise

图8给出了5倍随机噪声吸收光谱PCA滤波效果比较，噪声明显更大，第3个主成分开始特征接近为0，重建曲线与理想曲线相对偏差相对偏差0.2%左右，滤波效果良好。

图8 5倍随机噪声吸收光谱PCA滤波效果比较Fig.8 The filtering process of PCA algorithm to contaminated curve with 5-times noise

最后将本文提出的PCA滤波方法与传统的Savitzky-Golay滤波方法进行比较，用理想吸收光谱曲线与滤波后曲线的误差平方和(SSR)表征滤波效果，表1给出了2种滤波方法滤波结果比较。由表1可知，对于1倍随机噪声吸收光谱信号和5倍噪声吸收光谱信号，PCA滤波后SSR值分别从2.21和55.252降低到了0.002和0.018，噪声去除率达81%，且优于Savitzky-Golay滤波的0.006和0.006 2，充分证明本文所提出的PCA滤波方法对于不同噪声水平的吸收光谱曲线高频噪声滤除明显，且噪声去除率优于传统Savitzky-Golay滤波方法。

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

采用图3所示实验装置，对温度范围300～900 ℃，步长100 ℃，7个温度点H2O气体吸收光谱进行测量，采用直接吸收双线测温法进行温度解算，以热电偶温度传感器读数值为真值。温度解算过程中分别采用主成分滤波和Savitzky-Golay滤波两种方式进行光谱数据的预处理，最终温度解算结果如表2所示。主成分滤波后解算出的温度值与热电偶读数之间偏差明显更小，标准差从8.9降到5.3，进一步证明本文所提出的滤波方法滤波效果优于传统方法。

表2 基于不同滤波方法气体温度解算结果比较Tab.2 Comparison of gas temperature calculated by different filtering methods ℃

5 结 论

本文介绍了一种基于主成分分析的TDLAS吸收光谱高频噪声滤波技术，对TDALS双线测温原理、PCA滤波重建方法的原理进行了详细阐述。搭建了TDLAS温度测量装置，对实测吸收光谱曲线进行PCA滤波，实验结果表明，该方法滤波效果明显。为进一步对滤波效果进行量化和评估，对理想高斯吸收谱线叠加1倍随机噪声和5倍随机噪声，形成两种不同噪声水平的吸收光谱信号，用本文提出的PCA滤波方法和传统的Savitzky-Golay滤波进行高频噪声滤除，并通过理想吸收光谱曲线与滤波后曲线的误差平方和表征滤波效果进行比较。从试验结果可以看出，PCA滤波后曲线与理想曲线重合度较高，误差平方和分别降低到了0.002和0.018，噪声去除率达81%，且优于Savitzky-Golay滤波的0.006和0.006 2。