0 引　言

每年中国有上千万人因自然灾害造成生命财产损失，地震所造成的损失占比约为50%以上。为了保证人民的生命财产安全，寻找准确分析地震前兆的方法仍然是急需解决的问题。因为地壳运动会将CO₂等气体排出到地表，在地震前夕断裂带附近的CO₂地表浓度会出现异常，对地震带附近CO₂气体进行高精度测量，可以为分析地震前兆提供重要的帮助^[1-2]。

目前，主要测量CO₂浓度的方法有非分散性红外线技术(NDIR)分析技术、电化学技术、色谱分析技术等，但上述方法普遍存在易受其背景气体干扰，精度较低、无法实现实时监测等缺点^[3-5]。可调谐半导体激光吸收光谱技术(TDLAS)技术由于其具有不受背景气体干扰、精度高、可实时监测等优点，近年来已经成为国内外的研究热点，在气体检测领域得到广泛应用^[6]。2008年，Schaeffer等利用TDLAS技术对大气中CO₂浓度进行测量，现场测量得到的最低探测下限(MDL)为4.5 ppmv。2020年，王彪采用新型VCSEL激光器，研制了一套CO₂气体检测系统。系统最低探测下限达到9 ppmv，测量浓度为100 ppmv的CO₂系统误差为7.3%。2021年，吕淑媛利用分布反馈式激光器(DFB) 作为光源搭建了基于TDLAS技术的CO₂检测系统，其最低探测下限达到15 ppmv^[7-9]。然而，以上研究均没有对于TDLAS系统中存在的噪声进行分析并利用算法对系统噪声进行抑制处理，仪器的精度存在较大的提升空间。

文中首先介绍了TDLAS原理，然后介绍了CO₂检测系统的各部分组成，最后通过实验验证了采用卡尔曼小波算法的TDLAS系统具有高准确度、高稳定性、低探测下限的优秀性能，为高准确度CO₂系统预测地震前兆提供良好前景。

1 系统原理及设计

1.1 TDLAS原理

根据Beer-Lambert定律，对于单一频率激光，其通过气体吸收后的光强可表示为：

$ {I_\lambda } = {I_0}\exp [ - \alpha (\upsilon ,P,T)cL] $ (1)

式中： $ {I_0} $为入射光强； $ {I_\lambda } $为出射光强； $ \alpha $为 $ \upsilon ,P,T $的函数； $ \upsilon $为单色光频率； $ P $, $ T $为单色光所通过介质的压强与温度；c为气室中分子数浓度； $ L $为光程。

波长调制技术(WMS)在TDLAS技术基础上，引入高频正弦波作为调制信号，加入到激光器驱动信号中，产生的光信号通过多通池后经由锁相放大器对吸收谱线的二次谐波进行解调。当 $ \alpha (\upsilon ,P,T)cL $远小于0.05时，公式(1)变为：

$ {I_\lambda } = {I_0}[ - \alpha (\upsilon ,P,T)cL] $ (2)

经过简省后，二次谐波系数与浓度关系式可表示为：

$ {I_{2f}} \propto {I_0}{\alpha _0}cL $ (3)

式中： $ {\alpha _0} $是 $ \upsilon ,P,T $为定值时 $ \alpha (\upsilon ,P,T) $得到的常数。当 $ L $为定值时，由公式(3)可得到 $ {I_{2 f}} $信号与气体浓度成正相关关系^[10]。

1.2 系统设计

CO₂吸收谱线的选取如图1所示。

图 1. 波数为4978.202 cm⁻¹的CO₂吸收谱线

Fig. 1. CO₂ absorption spectra with wave numbers of 4978.202 cm⁻¹

下载图片 查看所有图片

为了实现痕量CO₂的高精度测量，同时避免与水蒸气吸收谱线重叠与交叉干扰，系统选取了较强的CO₂吸收谱线（4978.202 cm⁻¹）。系统通过调节注入DFB激光器的电流使得其输出的中心波长在2 008 nm并作为CO₂的检测光源。

TDLAS痕量CO₂检测系统原理图如图2所示。首先，将信号发生器产生的低频锯齿波信号和高频正弦波信号通过加法器叠加，提供给电流源用以扫描和调制DFB激光器波长。激光器输出的光束经过准直器后射入长光程多通池，出射光由光电探测器接收。

为了提高CO₂气体浓度的探测下限，该系统采用有效光程为40 m的自主研发柱面镜多通池。该多通池可在温度范围0~40 ℃及压强范围1.333~101.325 kPa的情况下稳定工作，以确保系统在现场测量过程中的可靠性。

多通池温度控制方面，温度采集电路通过传感器实时采集多通池的工作温度，通过16位A/D模块AD7606将多通池温度模拟信号转变为数字信号发送给STM32 F103 ZET6芯片，芯片中采集的温度值与设定值进行比较，运用PID算法对温度进行控制，后通过单片机发送给D/A模块，并输出相应的PWM控制 TEC驱动芯片，从而控制TEC实现对被控对象的进行温度控制。温度控制系统实验室内控制精密度可达到0. 01 ℃，调节时间为20 min。

多通池压力控制方面，待测气体经比例阀进入多通池，其次经过三通由压力传感器得到多通池中实时的压强值，再通过比例阀与气泵相接。压力传感器得到的实时压强值，以模拟量的形式传给压力控制电路。压力控制电路将实时压强值与目标压强值对比，改变驱动比例阀的PWM信号的占空比，改变比例阀的开口程度，从而实现对多通池内压强的控制。

图 2. TDLAS痕量CO₂检测系统原理图

Fig. 2. Schematic diagram of TDLAS trace amount of CO₂ detection system

下载图片 查看所有图片

电源模块采用外部24 V直流电源输入，可输出两路12 V/100 mA、一路24 V/1250 mA以及两路12 V/2 000 mA对系统进行供电，电源纹波均小于50 mV，实现了系统电源供给一体化。

采集和信号处理部分，系统测量采用STM32 F103 ZET6芯片进行控制与数据采集。在信号处理端，使用锁相放大器(LIA)解调光电探测器输出的信号，输出二次谐波信号由AD7606芯片转换为数字信号供STM32 F103 ZET6处理器采集与进行算法处理。

1.3 TDLAS系统噪声分析

探测器噪声和光学条纹干涉是影响TDLAS系统实际检测结果的重要因素。探测器噪声以热噪声和散粒噪声为主，均属于白噪声。系统的光学条纹干涉噪声受标准具长度、激光器温度、激光器驱动电流频率和调制幅度的影响，其在信号中的表达式推导如公式(4)。背景信号的二次函数可以展开为：

$ \begin{gathered} {S_{bg}}(\upsilon ,x) = \eta {I_1}\left(1 - \frac{F}{2}\right) + \eta \delta I\left(1 - \frac{F}{2}\right)\cos (\omega t - \phi ) + \\ \frac{{\eta {I_1}F}}{2}\cos (\varphi (\upsilon ,x)) + \frac{{\eta \delta IF}}{2}\cos (\omega t - \phi )\cos (\varphi (\upsilon ,x)) \\ \end{gathered} $ (4)

式中： $ {I_1} $为直流输出强度， $ \eta $为仪器响应参数； $ \delta I $为电流调制强度；x为标准具长度； $ \varphi (\upsilon ,x) $= $ 4\pi \upsilon {\text{x}} $； $ F $为精细度系数； $ \phi $为频率调制和强度调制的相位差。二次谐波信号主要来源于公式(4)中的第三项：

$ \frac{{\eta {I_1}F}}{2}\cos (\varphi (\upsilon ,x)) = \frac{{\eta {I_1}F}}{2}\cos (4\pi \upsilon x + 4\pi {\delta _\upsilon }x\cos \omega t) $ (5)

式中： $ {\delta _\upsilon } $为频率调制幅度；定义关于调制电流归一化的频率调制幅度 $ m = 2 x{\delta _\upsilon } $。将 $ m $代入公式(5)并进行贝塞尔展开，最终背景信号二次谐波 $ \;\mu (t) $可以表示为：

$ \mu (t) = \eta {I_1}F\cos (4\pi x{\upsilon _1}){J_2}(2\pi m) $ (6)

式中： $ {J_2} $为贝塞尔函数。通过上式可知，二次谐波背景信号中的光学条纹干涉会引起二次谐波波形以频率为 $ {\upsilon _1} $的低频余弦形式波动^[11]。这将在很大程度上影响系统精度。

2 实 验

2.1 卡尔曼-小波分析去除TDLAS噪声原理

为了解决TDLAS系统中的探测器噪声与光学条纹干涉的影响，文中采用了卡尔曼-小波去噪的方法。卡尔曼滤波技术是利用公式(6)更新估计和公式(7)观测估计两个基本理论方程来对待采集数据完成估值推测的：

$ {x_{k + 1}} = {x_k} + {\varepsilon _k} $ (7)

$ {z_k} = {x_k} + {\eta _k} $ (8)

式中： $ {x_{k + 1}} $为系统状态矩阵； $ {z_k} $为状态矩阵的实际观测量； $ {\varepsilon _k} $和 $ {\eta _k} $均属于系统探测器噪声引起的相互独立的随机变量。噪声观测估计方程将新测量的数据代入到预先的数据估计中从而得到更加接近于实际值的修正数据估计。利用卡尔曼滤波对于数据的预测性，可以有效的去除系统中的探测器噪声。二次谐波信号 $ {\mu _{(t)}} $经过卡尔曼滤波后，对其做小波变换并表示为：

$ {W_{\psi u}} = (a,b) = \frac{1}{{\sqrt a }}\int_{ - \infty }^{ + \infty } {u(t)} \psi \left(\frac{{t - b}}{a}\right){\rm{d}}t $ (9)

式中： ${\psi _{a,b}} = \dfrac{1}{{\sqrt a }}\psi \left(\dfrac{{t - b}}{a}\right)$为母小波时间平移b和尺度伸缩a的结果。由于小波变换可以分别在时域层面和频域层面进行对信号进行分解，可以根据光学条纹干涉在频域上频率为 $ \upsilon $与时域上近似余弦波波形的特征，对其进行分离，以达到去除光学条纹干涉的目的。综上，利用卡尔曼-小波去噪方法可以有效去除系统中存在的探测器噪声与光学条纹干涉的影响，从而提升系统的稳定性与准确度^[12-14]。

2.2 二次谐波的降噪处理

系统采用卡尔曼-小波分析方法消除探测器噪声以及光学条纹干涉的影响。图3为系统对浓度为50 ppmv CO₂样本吸收线的二次谐波采用卡尔曼-小波分析去噪方法前后效果对比图。

图 3. 利用卡尔曼-小波分析算法滤波前后效果对比

Fig. 3. Comparison of the second harmonic collected by the system before and after filtering by Kalman-wavelet analysis algorithm

下载图片 查看所有图片

图3上图中蓝色波形是滤波前的气体浓度二次谐波波形,红色波形是利用滤波后的波形。SNR由7.5提升到了15.5,是原信噪比的2.06倍。拟合残差1 $ \sigma $为0.38946 V。

2.3 气体浓度标定实验

图4展示了气体浓度标定实验的结果。

图 4. 50~8000 ppmv CO₂浓度标定图

Fig. 4. Diagram of 50-8000 ppmv CO₂ concentration calibration

下载图片 查看所有图片

实验室配气系统配备了五种不同浓度的CO₂气体(50、300、1000、4000、8000 ppmv)，每个浓度值均采样500点以上并取平均值，在实验室25 ℃的环境下对该仪器的气体浓度进行了标定实验。如图4所示，CO₂浓度标定线性度达到了0.997 7。

2.4 气体检测准确度对比实验

实验针对五种不同浓度的CO₂气体(50、300、1000、4000、8000 ppmv)进行测量，对该仪器应用卡尔曼-小波分析滤波前后的气体浓度检测准确度进行了测量与对比。

图5为卡尔曼-小波分析算法滤波前CO₂气体浓度误差对比图，黑色部分为滤波前的五种不同浓度CO₂气体测量误差曲线，误差为4.47%~4.83%。红色部分为滤波后的误差曲线，误差为2.57%~2.66%。可以看出，使用卡尔曼-小波分析滤波算法可以明显减小系统误差，提高检测准确度。

图 5. 卡尔曼-小波分析算法滤波前后CO₂误差

Fig. 5. Error result of Kalman-wavelet analysis algorithm

下载图片 查看所有图片

2.5 气体检测稳定性对比实验

系统在4000 ppmv浓度条件下进行了稳定性对比实验，结果如图6所示。

图 6. 系统在4000 ppmv浓度下CO₂ 检测稳定性对比实验

Fig. 6. Comparison experiment of CO₂ detection stability of the system at 4000 ppmv concentration

下载图片 查看所有图片

图6中红色部分为1000 s内采用卡尔曼-小波分析滤波前得到的浓度结果，波动范围为3 886.3~4 122.1 ppmv，标准差为35 ppmv。蓝色部分为系统滤波后探测到的浓度结果，波动范围为3 952.3~4 062.5 ppmv，标准差为20.9 ppmv。由图可知，系统采用卡尔曼-小波分析滤波后测得的浓度数据波动范围更小，得到数据的标准差由35 ppmv下降为20.9 ppmv。系统稳定性得到了较大提升。

2.6 系统探测下限对比实验

系统对浓度为300 ppmv的CO₂气体样本进行测量，然后利用Allan方差来分析探测下限。

如图7所示，图7(a)与图7(b)中呈下降趋势的红色实线表示系统在白噪声占主导的区域呈现的预期的响应。图7(a)表示系统直接测量的浓度值得到的艾伦方差与积分时间 $ \tau $的关系。由图可知，最优的积分时间约72 s时所对应的最低探测下限为7.2 ppmv。图7(b)表示系统通过卡尔曼-小波分析滤波后的浓度值得到的艾伦方差与积分时间t的关系。由图可知，最优的积分时间约61 s时所对应的最低探测下限为5.5 ppmv。实验利用卡尔曼-小波分析算法提升了系统的探测下限，使得最低探测下限由7.2 ppmv提升到5.5 ppmv。

图 7. Allan方差探测下限图

Fig. 7. Allan deviation method detection limit

下载图片 查看所有图片

3 结论与展望

文中研制了面向地震监测的高精度CO₂检测系统，并针对系统噪声利用卡尔曼-小波分析算法进行滤波，实验表明滤波后系统准确度、稳定性、最低探测下线指标均得到了提升。结合该方法的系统可以使地球化学气体测量具有更广泛的应用前景，为实现准确分析地震前兆提供重要的帮助。下一步将在目前的高精度CO₂检测系统研制基础上加入主成分分析(PCA)算法，通过对二次谐波主特性进行提取的方式进一步提升系统准确度与稳定性。