基于CEEMDAN的光声光谱CO检测技术

李磊，汤利苹，马秋阳，高子健，高杨，乔莹莹

（郑州大学 物理（微电子）学院，郑州 450001）

0 引言

一氧化碳（CO）是一种无色、无味的有毒气体，主要由含碳材料的不完全燃烧产生［1-2］。CO被吸入人体后，使血红蛋白丧失携氧能力，导致人体各组织器官缺氧甚至死亡［3］。鉴于其危害性，世界卫生组织（World Health Organization，WHO）现行空气质量指南中规定CO在15 min内暴露量上限为87.3×10-6，1 h内暴露量上限为30.6×10-6［3］。因此，实现对CO的实时、高灵敏度检测对保障日常生活和工作场所的安全极为重要。

目前，光声光谱（Photoacoustics Spectroscopy，PAS）检测技术因其高灵敏度、实时、高动态范围等优点，在痕量气体检测领域受到广泛关注［4-6］。PAS是一种基于光谱吸收的技术，调制的光信号照射待测气体，待测气体吸收光能后产生与调制光同频的声波，通过声传感器检测的声信号可反演待测气体浓度［7］。根据声传感器的类型可将光声光谱分为传统共振式PAS、悬臂梁增强型光声光谱（Cantilever Enhanced PAS，CEPAS）［8］和石英音叉增强型光声光谱（Quartz Enhanced PAS，QEPAS）［9］。传统共振式PAS利用共振式光声池放大声信号，并用商用麦克风检测声波信号，技术成熟，性价比高。CEPAS利用悬臂梁检测声信号，灵敏度高，响应频带宽，但成本较高。QEPAS利用石英音叉探测声信号，体积小，品质因数高。但石英音叉工作频率固定，使得该技术的响应频带较窄。

鉴于PAS的优异性能，近年来基于PAS的CO检测技术也得到了广泛研究。2013年，MA Yufei等采用987 mW的4.61 μm激光器结合QEPAS，获得了1.5×10-9（体积分数，下同）的最小检测限［10］；2016年，MAO Xuefeng等采用法布里-珀罗光纤声传感器结合近红外光源，最小检测极限达到4.6×10-6［11］；2018年，CUI Ruyue等使用2.33 μm激光器和14.5 m多反射池，最小检测极限达到1×10-6［12］；2019年，YIN Xukun等报道了一种采用1.56 μm激光器结合10 W光纤放大器的CO传感器，获得了10-9级别的最小检测极限［13］；2020年，LIU Xin等采用2.3 μm激光器作为激发光源，结合水汽增强光声信号的方法，实现了1.88×10-7的最小检测限［14］。PAS技术中可通过选择中红外光源［8-9］、增大光功率［8，13］、增大吸收光程［11］、提高弛豫率［8，14］等方式提高系统的信噪比。此外，通过算法降噪实现信噪比的提升是另一种常用的改善系统性能的方法，能在不增加系统复杂度的基础上提高系统性能，对于提高PAS痕量气体检测性能具有重要的研究意义。

在PAS系统中，气体吸收峰与波长调制的激光相互作用形成光声信号，是一个非线性非平稳的过程，光声信号同时受到环境中广泛存在的乘性噪声、池壁和窗口吸收产生的背景噪声等随机和非随机噪声的干扰［15-17］。为降低噪声干扰，信号平均算法［18］、卡尔曼滤波［4］、回归算法［19］、小波变换［20-21］等方法被提出且成功应用于PAS领域。但上述算法主要针对线性非平稳或者非线性平稳信号，在处理非线性非平稳信号时有局限性，因此，PAS系统需要引入非线性非平稳算法对光声信号进行处理。经验模态分解（Empirical Mode Decomposition，EMD）因其优异的非线性非平稳信号处理能力，自提出以来就得到广泛应用［22］。为了降低EMD中模态混叠效应，TORRES M E等提出了基于自适应噪声的完备经验模态分解（Complete Ensemble EMD with Adaptive Noise，CEEMDAN）［23］。

本文提出了一种基于PAS痕量气体检测技术的CO传感器。该传感器采用1.56 μm分布反馈式激光器和掺铒光纤放大器（Erbium-Doped Optical Fiber Amplifier，EDFA）作为激发光源系统，设计了一种差分光声池（Differential Photoacoustic Cell，DPAC）来消除气流噪声和环境噪声。在信号处理过程中，采用以CEEMDAN为核心的降噪重构算法处理光声信号以获得较高的检测灵敏度。

1 基本原理

EMD是一种自适应的时域信号分解方法，可以将非线性、非平稳信号分解成有限个本征模态函数（Intrinsic Mode Function，IMF），每个IMF分量代表不同时间尺度特征。然而，EMD分解常存在模态混叠的问题，CEEMDAN方法通过自适应加入白噪声解决了此问题［24］。

设x(t)为待分解时域信号，vj(t)为满足标准正态分布的高斯白噪声，Ek(·)为EMD分解得到的第k个模态分量，Ik为CEEMDAN分解得到的第k个模态分量。CEEMDAN算法的具体分解过程为［25-26］

1）首先，向原始信号x(t)添加高斯白噪声vj(t)，获得第j次添加高斯白噪声后的时间序列xj(t)。对xj(t)进行EMD分解，并对产生的N个模态分量进行总体平均得到CEEMDAN的第一个本征模态分量［22］，即

式中，j为加入白噪声的次数，j=1，2，…，N，ε0为噪声系数。

2）一阶残差r1(t)为原始信号x(t)和一阶模态分量I1(t)之差为

3）在一阶残差r1(t)的基础上继续添加高斯白噪声vj(t)，获得新信号r1j(t)，进行EMD分解。重复计算N次，对N个EMD分解得到的模态分量集成平均可得CEEMDAN的第二个本征模态分量，即

4）二阶残差r2(t)为新的原始信号r1(t)和二阶本征模态分量I2(t)之差，即

5）重复步骤3）～4），直到获得的残差信号rK(t)为单调信号不能再分解为止。此时得到的本征模态分量的数量为K，原始信号可表示为

在CEEMDAN的基础上结合Savitxky-Golay（S-G）滤波算法，对CEEMDAN分解的模态分量进行进一步降噪处理，在去除噪声的同时最大程度地保留原始信号的波形特征［27］。S-G滤波器利用最小二乘拟合系数作为滤波器响应函数，是一种具有卷积平滑的高频噪声全自适应滤波方法，已经被成功应用于光谱降噪［20，28-29］。计算降噪后各分量和原始信号之间的相关系数，对于各次本征模态分量，若相关系数大于0.4，则认为此模态分量由信号成分主导，选择用来重构信号，否则，认为此分量由噪声成分主导，直接舍弃。残差分量选择用来重构信号。本文中，将CEEMDAN分解，S-G滤波再结合相关系数重构的降噪算法称为以CEEMDAN为核心的降噪重构算法。

为了说明以CEEMDAN为核心的降噪重构算法的性能优势，将其与光声光谱领域常用的降噪方法——小波降噪进行性能比较。基于光声光谱基本原理，利用MATLAB仿真CO吸收光能产生的光声信号的二次谐波并归一化。选择最常见的高斯白噪声添加到纯二次谐波波形上，分别计算含噪信号、小波降噪信号和以CEEMDAN为核心的降噪重构算法降噪信号的信噪比（Signal to Noise Ratio，SNR），以及上述信号与纯二次谐波信号之间的结构相似度（Structural Similarity，SSIM）［20，30］。信号的SNR表达式为

式中，Ps为信号的有效功率，Pn为噪声的有效功率。

结构相似度SSIM是衡量两幅图片相似度的指标，被用作描述含噪信号、小波降噪信号、以CEEMDAN为核心的降噪重构算法降噪信号与纯二次谐波信号之间的相似程度。SSIM的取值范围为0～1，数值越接近1，说明两个信号之间的相似程度越高。SSIM的计算公式为

式中，x为原始信号，即纯二次谐波波形；y为降噪后重构信号；μx、μy分别为x和y的平均；σxy为x的协方差；代表方差；c1和c2接近于零，是正则化常数，用来避免分母的不稳定性。

降噪效果如图1所示，以CEEMDAN为核心的降噪重构算法降噪后的波形更加贴合纯二次谐波信号波形。计算结果如表1所示，以CEEMDAN为核心的降噪重构算法处理后信号的SNR和SSIM分别为18.17和0.66，大于小波降噪后信号的SNR和SSIM（10.98，0.53）。因此，在处理光声信号时，以CEEMDAN为核心的降噪重构算法性能明显优于小波降噪方法。

图1 仿真降噪性能Fig.1 Simulation of noise reduction performance

表1 小波和以CEEMDAN为核心的降噪重构算法的降噪效果对比Table 1 Performance comparison of wavelet algorithm and the noise reduction and reconstruction algorithm based on CEEMDAN

2 实验系统

2.1 吸收峰的选择

根据HITRAN光谱数据库，CO在红外波长范围内有三个吸收谱带，吸收谱线的位置和对应的吸收系数如图2（a）所示［31］。可以看出，CO近红外的吸收系数比中红外的吸收系数低3～4个数量级。但是，近红外1.58 μm处的吸收带正好位于光通信波段，得益于光通信技术的发展，该波段的激光器以及相关的光放大技术比较成熟且性能稳定，成熟的近红外光放大技术能够补偿由弱吸收线导致的灵敏度降低问题，因此本系统选择该波段执行气体检测。为避免交叉干扰影响检测性能，在选择吸收峰时，需避开空气中的常见气体如H2O、CO2、H2S在近红外波段的吸收谱线。依据HITRAN数据库对纯气在296 K和1 atm（标准大气压，1 atm=101325 Pa）条件下的吸收系数进行仿真计算，考虑到空气中H2O含量大，为表征真实情况，将H2O的吸收系数放大100倍后与相应位置的CO吸收谱线的吸收系数进行比较。如图2（b）所示，位于1568.04 nm的CO吸收峰吸收系数大且无交叉干扰，因此被用于CO气体测量。在PAS系统中，光声信号的幅值正比于系统激励光源的光功率，可以通过放大光功率提高系统信噪比。系统选用的分布反馈式（Distributed Feedback，DFB）激光器（G&H E0067929）的出射激光功率有限，使用EDFA放大光功率至300 mW［11，31-33］。

图2 CO吸收系数图谱Fig.2 CO absorption coefficient diagram

2.2 光源调制频率确定

光声池的共振频率是谐振式光声气体传感器的重要参数，光源调制频率应调谐到光声池的一个本征频率处，使谐振腔内形成驻波，从而实现光声信号的放大［34］。为了准确地测量光声池共振频率，利用质量流量控制器（MFC，MF-200C，INHA）将1% CO/N2标准气体以恒定流速充入DPAC，在室温常压的条件下测量光声池的共振频率。正弦调制电流的频率从900 Hz持续调节到1800 Hz，记录产生的二次谐波信号的幅值，用洛伦兹曲线拟合实验数据得到光声池频率响应曲线，如图3（a）所示。可知，二次谐波信号的幅值随着正弦频率的增加先增大后减小，在1512 Hz达到峰值，因此推知光声池的一阶纵向共振频率为调制频率的2倍，即3024 Hz。同时，充分考虑进出气口和麦克风与谐振腔之间的空气间隙等因素，在COMSOL中对光声池结构做一比一的建模。基于一维谐振腔理论和热粘滞声学理论，利用有限元分析方法（Finite Element Method，FEM）仿真光声池的一阶纵向共振模态如图3（b）所示。光声信号在谐振腔内形成驻波，谐振腔中部的声压最大，并且两个谐振腔内声压相位相反。由于光声池内两个谐振腔结构完全一样，内部产生一样的气流噪声和环境噪声。对两个谐振腔上麦克风测得的信号使用差分模式能够有效去除噪声、放大光声信号，从而提高信噪比。

图3 工作在一阶纵向共振模式的差分光声池Fig.3 Differential photoacoustic cell operating in a first-order longitudinal resonance mode

2.3 调制深度优化

对于PAS系统，激光器驱动电流的调制幅度严重影响系统性能［35-36］。为了确定激光器驱动电流的最佳调制幅度，室温常压下将1% CO/N2气体混合物作为样品气体注入光声池内。图4（a）中绘制了二次谐波信号幅值随激光器驱动电流调制幅度变化的趋势，说明随着调制幅度的增加，二次谐波信号的幅值会趋于不变。图4（b）展示了不同调制幅度下，二次谐波的波形。可知，当调制电流幅度增大到1.6 V以上时，虽然二次谐波信号幅值没有随着调制电流幅度的增大而下降，但信号波形由于残余强度调制（Residual Amplitude Modulation，RAM）效应导致了严重扭曲［37］。因此，为了保证在信号不失真的前提下获得最大的二次谐波信号幅值，最终选定调制幅度为1.6 V，并将此最优调制幅度用于后续实验。

图4优化调制幅度Fig.4 Optimized modulation amplitude

2.4 系统装置

图5 （a）为PAS痕量气体检测系统原理，主要功能模块包括：光源及驱动，EDFA，光纤准直器，光声池，锁相放大模块，数据处理模块。在实验过程中，为了降低由窗口和池壁吸收引入的相干背景噪声的影响，使用波长调制和二次谐波检测技术。信号发生器（Fluke 294）产生低频锯齿波和高频正弦波的叠加信号调制激光波长，使激光波长缓慢地扫过CO吸收峰，正弦波的频率为光声池一阶纵向共振频率的一半。此叠加信号送入激光器驱动电流源（ILX Lightwave LDX-3232），并与激光器温度控制器（ILX Lightwave LDT-5525B）一起调控激光器的输出波长。激光器产生的光信号送入定制的EDFA中并将功率泵浦到～300 mW。为了减少池壁吸收，EDFA出射的激光光束经光纤准直器耦合进入DPAC。利用嵌入在光声池共振腔中部的麦克风（BSWA MPA416570074）探测产生的光声信号，信号经前置放大器后，送入锁相放大器（Stanford Research Systems，USA，Model SR830）进行解调。

图5 PAS系统及差分光声池Fig.5 The PAS system and differential photoacoustic cell

双通道DPAC的结构可有效降低系统中气流噪声和环境噪声的影响，改善系统性能。实验所用的光声池为树脂材料3D打印的双通道DPAC，谐振腔直径为6 mm，长度为30 mm，谐振腔两端连接的缓冲室长度为15 mm，直径为30 mm，用于减轻气流噪声。信号通道和参考通道谐振腔的中间部分嵌入相同型号的麦克风用来收集光声信号。低噪声DPAC具体结构如图5（b）所示。

3 实验结果与讨论

传感器系统中存在气体流速噪声、池窗口及池壁强吸收产生的背景噪声以及电子设备自身产生的其他干扰噪声等，为了提高系统信噪比，使用以CEEMDAN为核心的降噪重构算法抑制噪声。测量记录了1500×10-6的CO/N2混合气体通入光声池时产生的二次谐波信号，对其进行CEEMDAN分解，加入的白噪声标准差为0.2，噪声添加次数为50。经CEEMDAN分解得到13个IMF分量和一个残余分量如图6（a）所示。在图6（a）中，各次IMF分量Ii的频率从高到低排列，最后一个分量r为表征波形趋势的残余分量。对分解得到的IMF分量和残余分量进行S-G滤波处理得到Di和Dr，如图6（b）。表2中列出了降噪后的IMF分量和残余分量与原始二次谐波信号之间的相关系数Ci/Cr，用降噪后分量和原始信号之间的相关系数自适应选取重构信号的分量。

图6 IMF分量Fig.6 The IMF components

表2 分量与原始信号之间的相关系数Table 2 The correlation coefficients of the components and original signal

原始信号与经以CEEMDAN为核心的降噪重构算法重构后的信号对比如图7（a）所示。重构信号在保留了峰值和宽度等波形信息的基础上去除了噪声成分。原始噪声信号和经小波降噪、以CEEMDAN为核心的降噪重构算法处理后的噪声信号对比如图7（b）所示。光声光谱中信号的信噪比为二次谐波幅值与噪声标准差（σ）之比，取二次谐波旁边无信号区域作为噪声信 号［21］。系 统 的 最 低 检 测 极 限（Minimum Detection Limit，MDL）可由特定气体浓度和信噪比给出，即cmin=c/SNR，式中，cmin为最低可测量浓度，c为特定气体浓度，SNR为信噪比，具体数据在表3中列出。经小波降噪后信号的信噪比由128.2提高到213.6，是原来的1.7倍，最低可检测浓度由11.7×10-6降至7.0×10-6。经以CEEMDAN为核心的降噪重构算法处理后信号的信噪比由128.2提高到585.8，是原来的4.6倍，最低可检测浓度由11.7×10-6降至2.6×10-6。与仿真结果一致，实验结果也证实了以CEEMDAN为核心的降噪重构算法应用于PAS领域的可行性与有效性，若要进一步降低PAS传感器对CO气体的检测极限，可通过加大激发光源功率等手段实现［38-39］。

表3 光声光谱系统性能参数Table 3 Performance parameters of PAS system

为了标定系统性能，室温常压下将不同浓度的CO/N2混合气体充入光声池中，测定产生的光声信号。将激光器的电流和温度分别设定为280 mA和39.3℃，使DFB激光器的输出波长锁定CO吸收峰1568.04 nm。用由0.1 Hz的锯齿波和1512 Hz的正弦波组成的叠加信号控制激光器的电流驱动，对输出激光波长进行调制。利用双通道配气系统把不同浓度的CO/N2混合气体以200 sccm的流速充入光声池内。谐振腔内产生的声波被麦克风拾取后，通过前置放大器，送入锁相放大器中进行解调放大，记录锁相放大器积分时间为100 ms时输出的光声信号。图8绘制了二次谐波信号的峰值与CO浓度之间的关系。可以看出，降噪前后二次谐波信号的峰值与气体浓度的关系呈现良好的线性关系，与理论相符合。图8中的插图为不使用/使用以CEEMDAN为核心的降噪重构算法后长时间的噪声，可以看到采用以CEEMDAN为核心的降噪重构算法后噪声标准差大大降低，意味着系统的精度和稳定性都得到了明显提高。

图8 气体浓度和二次谐波幅值之间的关系（插图为纯N2背景时测量的噪声）Fig.8 The relation between gas concentrations and second harmonic signal amplitudes（inset shows the noise measured at the pure N2 background）

4 结论

本文设计了一种基于PAS技术的痕量CO气体传感器。使用通信波段性能稳定的近红外分布反馈式激光器和EDFA作为系统激发光源，采用双通道结构的DPAC降低系统噪声。利用以CEEMDAN为核心的降噪重构算法处理采集的非线性光声信号，结合S-G滤波器对CEEMDAN分解产生的分量进一步降噪，计算降噪后的分量和原始信号的相关系数自适应地选取重构信号的有效分量，从而提高信噪比，优化系统性能。实验结果表明，在积分时间为100 ms时，以CEEMDAN为核心的降噪重构算法将系统信噪比提高到原来的4.6倍，将系统的最小检测极限从11.7×10-6降到了2.6×10-6，该算法可应用于PAS痕量气体检测领域。