CN105866068A - 一种基于递推最小二乘法的激光吸收光谱谐波信号峰值点提取方法 - Google Patents

Abstract

本发明提供一种基于递推最小二乘法的激光吸收光谱谐波信号峰值点提取方法，分为以下步骤：(1)以激光吸收光谱的谐波信号采样点和采样点的序列数作为输入，使用三次多项式拟合谐波信号的主峰部分；(2)使用第一组输入数据计算初始值；(3)增加测量数据代入递推公式进行计算；(4)使用的递推结果计算拟合函数峰值点；(5)对计算得的峰值点进行判断，若符合条件，拟合函数的峰值点即为谐波信号的峰值点，若不符合条件，则返回(3)。本发明用于激光吸收光谱监测系统中谐波信号峰值点的查找，使监测系统在准确测量的基础上大幅提高测量速度，具有广泛的应用前景。

Description

一种基于递推最小二乘法的激光吸收光谱谐波信号峰值点提取方法

技术领域

本发明涉及一种基于递推最小二乘法的激光吸收光谱谐波信号峰值点提取方法，具体涉及在激光吸收光谱测量技术中，使用递推最小二乘法快速、准确地提取谐波信号的峰值点。

背景技术

激光吸收光谱测量技术广泛应用于工业领域。例如，污染物的排放检测，燃烧产物测量，航空发动机性能监测等，可以测量温度、压力以及多种气体浓度。专利“一种车载式机动车尾气中氧气浓度的激光检测系统”(专利号：CN201510680434.5)使用该技术测量了机动车尾气中的氧气浓度。专利“基于紫外激光吸收光谱的非接触式火焰温度及OH基浓度测量装置及测量方法”(专利号：201410225315.6)使用该技术测量了火焰温度及OH基浓度。Beer-Lambert定律是激光吸收光谱测量技术的基础，该技术多应用在近红外光谱区域，在理想情况下，假设吸收介质为均匀介质，沿光程的气流参数相同，且入射激光为单一频率的前提下，关系式表示为：

$I_v=I_{v,0}\·\exp (-P\·X_{H_{2}O}\·S(T)\·\mathrm{\φ}(v)\·L),---\left(1\right)$

其中，v[cm^-1]为入射光频率，I_v为透射光强，I_v，0为入射光强，L[cm]为均匀待测气流长度，即吸收长度，P[atm]为静压，X为吸收组分摩尔浓度，S(T)[cm^-2atm^-1]为吸收线的线强度，φ(ν)[cm]为线型函数。其中线型函数满足归一化条件，

激光吸收光谱测量技术中，通过在激光器上加载一个较慢频率的斜坡信号与高频正弦波信号对激光器进行驱动，激光器的输出光强为 为未调制时的激光器输出光强，i₀为光强调制幅值。接着使用锁相放大器分析通过吸收气体后的激光信号，可以获得调制频率多次谐波信号。因为二次谐波信号为对称信号，在线型中心处有最大值，且为偶次谐波信号中幅值最大的谐波信号，所以分析中最常用的为调制频率的二次谐波信号。根据论文“Second-Harmonic Detection withTunable Diode Lasers Comparison of Experiment and Theory”(Applied Physics B,1981(26):pp.203-210.)，理想二次谐波H₂表达式为：

$H_2(\stackrel{\‾}{v},a)=-\frac{S\left(T\right)\·P\·X_{H_{2}O}\·L}{\mathrm{\π}}{\∫}_{\mathrm{\π}}^{+\mathrm{\π}}\mathrm{\φ}(\stackrel{\‾}{v}+acos\mathrm{\θ})\·cos\left(2\mathrm{\θ}\right)\·d\mathrm{\θ},---\left(2\right)$

其中，a[cm^-1]为光强调制幅值i₀所对应的激光频率调制幅值。实际经过锁相放大器所得信号并不是标准的二次谐波信号，且谐波信号中受激光噪声，激光散射等影响，包含一定的背景噪声。根据论文“Wavelength modulation absorption spectroscopy with 2f detection using multiplexed diode lasers forrapid temperature measurements in gaseous flows”(Applied Physics B:Lasers and Optics,2004.78(3-4):pp.503-511)，实际测得二次谐波信号与理想的谐波信号峰值相等，将二次谐波除以对应的一次谐波可以去除背景噪声对谐波峰值的影响。

根据表达式(2)可以知道，在温度和压力已知的情况下，气体浓度与二次谐波峰值成正比。在使用已知浓度的气体对二次谐波峰值进行标定后，被测气体浓度可以根据以下公式算出：

$X_{mex}=X_{calib}\frac{H_2{(\stackrel{\‾}{v},a)}_{meas}}{H_2{(\stackrel{\‾}{v},a)}_{calib}},---\left(3\right)$

其中，下标“meas”表示被测值，下标“calib”表示标定值

在一般的实验条件下，需要将谐波信号上传到上位机，经过遍历查找一个周期内的峰值，例如专利“用可调谐激光吸收光谱的气相H₂O₂浓度的测量方法及装置”(专利号：CN201310561962.X)，“一种基于激光吸收光谱的气体浓度在线测量方法”(专利号：CN201210012470.0)，这样大大降低了激光吸收光谱的测量速度。而递推最小二乘法是实时参数辨识的有效方法，如专利“基于递推最小二乘法的空调所属建筑物一阶模型实时参数辨识方法”(专利号：CN201510571379.6)和“基于递推最小二乘法RLS的结晶器ARX模型辨识方法”(专利号：CN201210066222.4)。因此，研究激光吸收光谱技术和递推最小二乘法相结合的方法，可以在有效提高测量速度的基础上，进一步提高测量的精度。

基于以上背景，本文发明了一种基于递推最小二乘法的激光吸收光谱谐波信号峰值点提取方法，通过对谐波信号一个周期中的部分采样点进行在线拟合获得实时峰值点信息。在减少噪声对测量结果影响的同时，提高了激光吸收光谱技术的测量速度。

发明内容

针对激光吸收光谱技术，为了能在线高速高准确的提取谐波信号的峰值点，本发明提供一种基于递推最小二乘法的激光吸收光谱谐波信号峰值提取方法

本发明所采用的技术方案如下：

步骤一，以激光光谱吸收信号的谐波信号采样点s(n)和采样点在一个周期中的序列数n为输入，写出对谐波信号的主峰部分进行拟合的三次多项式：

在白噪声影响下，所述三次多项拟合函数可表示为

s(n)＝u(n)·x+e(n)， (1)

其中，e(n)为白噪声，

$\begin{array}{c}u\left(n\right)=\left[\begin{array}{cccc}1& n& n^2& n^3\end{array}\right],\\ x={\left[\begin{array}{cccc}{x}_0& x_1& x_2& x_3\end{array}\right]}^T,\end{array}---\left(2\right)$

将所述N(N≥1)个相邻输入s(n)和n表示为式(3)(4)的矩阵形式

S_N＝[s(1) s(2)…s(N)]^T， (3)

$U_N=\left[\begin{array}{c}u\left(1\right)\\ u\left(2\right)\\ .\\ .\\ .\\ u\left(N\right)\end{array}\right]=\left[\begin{array}{cccc}1& 1^1& 1^2& 1^3\\ 1& 2^1& 2^2& 2^3\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& N& N^2& N^3\end{array}\right],---\left(4\right)$

则待求参数x_N可以通过式(5)计算出

$x_N=J_N{U}_N^T{S}_N,---\left(5\right)$

其中

$J_N={\left(U_N^T{U}_N\right)}^{-1},---\left(6\right)$

步骤二，计算递推初始值：

使用第一组输入数据(N＝1)计算递推初始值，

$\left\{\begin{array}{c}{J}_1={\left(U_1^T{U}_1\right)}^{-1}\\ x_1=J_1{U}_1^T{S}_1\end{array}\right.,---\left(7\right)$

步骤三，增加测量数据代入递推公式进行计算：

将所述S_N+1和U_N+1表示为递推形式

$S_{N+1}=\left[\begin{array}{c}{S}_N\\ s\left(N\right)\end{array}\right],U_{N+1}=\left[\begin{array}{c}{U}_N\\ u\left(N\right)\end{array}\right],---\left(8\right)$

则

$\begin{array}{c}{J}_{N+1}={\left(U_{N+1}^T{U}_{N+1}\right)}^{-1}={\left({\left[\begin{array}{c}{U}_N\\ u\left(N\right)\end{array}\right]}^T\left[\begin{array}{c}{U}_N\\ u\left(N\right)\end{array}\right]\right)}^{-1}\\ ={\[U_N^T{U}_N+u^T\left(N\right)u\left(N\right)\]}^{-1}\\ ={\[J_N^{-1}+u^T\left(N\right)u\left(N\right)\]}^{-1}\end{array},---\left(9\right)$

也可写成

$J_N^{-1}=J_{N+1}^{-1}-u^T\left(N\right)u\left(N\right),---\left(10\right)$

根据矩阵求逆定理

(B+C·D·E)^-1＝B^-1-B^-1C(D^-1+E·B^-1·C)^-1E·B^-1， (11)

令

$B=J_N^{-1},C=u^T\left(N\right),D=I,E=u\left(N\right),---\left(12\right)$

其中，I为单位矩阵，所以式(9)可以表示为式(13)所示的递推公式

$\begin{array}{c}{J}_{N+1}={\[J_N^{-1}+u^T\left(N\right)\·I\·u\left(N\right)\]}^{-1}\\ =J_N-J_N{u}^T\left(N\right){\[I+u\left(N\right)J_N{u}^T\left(N\right)\]}^{-1}u\left(N\right)J_N\end{array},---\left(13\right)$

结合式(5)，式(8)和式(13)，x_N+1可以表示成

$\begin{array}{c}{x}_{N+1}={\left(U_{N+1}^T{U}_{N+1}\right)}^{-1}{U}_{N+1}^T{S}_{N+1}=J_{N+1}{\left[\begin{array}{c}{U}_N\\ u\left(N\right)\end{array}\right]}^T\left[\begin{array}{c}{S}_N\\ s\left(N\right)\end{array}\right]\\ =J_{N+1}{J}_N^{-1}\left(J_N{U}_N^T{S}_N\right)+J_{N+1}{u}^T\left(N\right)s\left(N\right)\\ =J_{N+1}{J}_N^{-1}{x}_N+J_{N+1}{u}^T\left(N\right)s\left(N\right).\end{array},---\left(14\right)$

将式(10)代入式(14)，得到第二个递推公式

$\begin{array}{c}{x}_{N+1}=J_{N+1}\[J_{N+1}^{-1}-u^T\left(N\right)u\left(N\right)\]x_N+J_{N+1}{u}^T\left(N\right)s\left(N\right)\\ =x_N+J_{N+1}{u}^T\left(N\right)\[s\left(N\right)-u\left(N\right)x_N\]\end{array},---\left(15\right)$

经过以上计算过程，得到递推公式为，

$\left\{\begin{array}{c}{J}_{N+1}=J_N-J_N{u}^T\left(N\right){\[I+u\left(N\right)J_N{u}^T\left(N\right)\]}^{-1}u\left(N\right)J_N\\ x_{N+1}=x_N+J_{N+1}{u}^T\left(N\right)\[s\left(N\right)-u\left(N\right)x_N\]\end{array}\right.,---\left(16\right)$

步骤四，使用递推结果计算拟合函数的峰值点：

可根据递推结果x获得拟合函数的各项系数x_k(k＝0,1,2,3)，所述拟合函数的导数为

s′(n)＝x₁+2x₂n+3x₃n²， (17)

令s′(n)＝0，求得峰值点的位置

$N_{max}=\frac{-2x_2-\sqrt{4x_2^2-6x_1{x}_3}}{6x_3},---\left(18\right)$

求得峰值点的大小

P＝s(N_max)， (19)

步骤五，对计算得的峰值点进行判断：

当N≥20，对N-19，N-18，…，N这20组结果设定如下判断条件

$\begin{array}{c}\frac{1}{20}\sqrt{\underset{i=N-19}{\overset{N}{\Σ}}{\[N_{\max }\left(i\right)-{\stackrel{\‾}{N}}_{\max }\]}^2}\≤\frac{1}{1000}{N}_{\max }\left(N\right)\\ \frac{1}{20}\sqrt{\underset{i=N-19}{\overset{N}{\Σ}}{\[P\left(i\right)-\stackrel{\‾}{P}\]}^2}\≤\frac{1}{1000}P\left(N\right)\end{array},---\left(20\right)$

其中，N_max(i)，P(i)表示计算的峰值点位置与峰值点大小，为20组结果的平均值，若符合条件，则拟合函数的峰值点(N_max(N)，P(N))即为谐波信号的峰值点，若不符合条件，则返回步骤三。

本发明的效果：通过使用递推最小二乘法，可以在线获得谐波信号的峰值点信息，有效地提高测量速度。同时，基于递推最小二乘法使用拟合函数对谐波信号进行拟合，使得计算的峰值点信息更加准确，从而提高激光吸收光谱技术的测量速度和准确性。

附图说明

图1是基于激光吸收光谱技术的测量系统结构原理图。

图2是递推最小二乘法流程图。

图3是测量系统测得的谐波信号示意图，矩形框中为需要递推最小二乘法进行拟合的部分。

图4是谐波信号与通过递推最小二乘法(Recursive Least Square,RLS)计算的拟合信号之间的对比图，N_max和P为测量所需的峰值点信息。

图5是设定CO₂浓度与使用递推最小二乘法(RLS)的CO₂浓度测量结果的对比图。

具体实施方式

在本实施案例中，通过预设气体浓度，然后利用该条件下的激光吸收光谱测量数据进行实验，证明了该方法的有效性，并在最后将测得的气体浓度与设定浓度进行比较。

下面结合附图对本发明作进一步的说明：

步骤一，将一定比例的CO₂、N₂气体充入增强腔中，选取CO₂吸收谱线进行实验研究，中心频率为6357.312cm^-1。增强腔中的吸收路径为50m，所有的实验都在常温常压下进行。图1为基于激光吸收光谱技术的测量系统结构原理图；

步骤二，使用激光控制器对激光器进行温度和电流的控制，实现对激光器输出频率的调谐；

步骤三，使用光纤分路器将激光分为两束，其中一束通过增强腔获得带有气体浓度信息的激光吸收光谱，另外一束通过标准具获得激光输出频率与时间的关系，每路都由光电二极管组成的光电探测器所检测；

步骤四，将通过被测气体的激光检测信号输入锁相放大器获得一次谐波信号和二次谐波信号，并将二次谐波信号除以对应的一次谐波信号实现标准化处理；

步骤五，标准化处理过的二次谐波就是递推最小二乘法中输入信号，对输入信号进行采样获得s(n)与采样点的序列数n，写出对谐波信号的主峰部分进行拟合的三次多项式：

在白噪声影响下，所述三次多项拟合函数可表示为

s(n)＝u(n)·x+e(n)， (1)

其中，e(n)为白噪声，

$\begin{array}{c}u\left(n\right)=\left[\begin{array}{cccc}1& 0& n^2& n^3\end{array}\right],\\ x={\left[\begin{array}{cccc}{x}_0& x_1& x_2& x_3\end{array}\right]}^T,\end{array}---\left(2\right)$

将所述N(N≥1)个相邻输入s(n)和n表示为式(3)(4)的矩阵形式

S_N＝[s(1) s(2)…s(N)]^T， (3)

$U_N=\left[\begin{array}{c}u\left(1\right)\\ u\left(2\right)\\ .\\ .\\ .\\ u\left(N\right)\end{array}\right]=\left[\begin{array}{cccc}1& 1^1& 1^2& 1^3\\ 1& 2^1& 2^2& 2^3\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& N& N^2& N^3\end{array}\right],---\left(4\right)$

则待求参数x_N可以通过式(5)计算出

$x_N=J_N{U}_N^T{S}_N,---\left(5\right)$

其中

$J_N={\left(U_N^T{U}_N\right)}^{-1},---\left(6\right)$

步骤六，计算递推初始值：

使用第一组输入数据(N＝1)计算递推初始值，

$\left\{\begin{array}{c}{J}_1={\left(U_1^T{U}_1\right)}^{-1}\\ x_1=J_1{U}_1^T{S}_1\end{array}\right.,---\left(7\right)$

步骤七，增加测量数据代入递推公式进行计算：

所述递推公式为：

$\left\{\begin{array}{c}{J}_{N+1}=J_N-J_N{u}^T\left(N\right){\[I+u\left(N\right)J_N{u}^T\left(N\right)\]}^{-1}u\left(N\right)J_N\\ x_{N+1}=x_N+J_{N+1}{u}^T\left(N\right)\[s\left(N\right)-u\left(N\right)x_N\]\end{array}\right.,---\left(8\right)$

步骤八，使用递推结果计算拟合函数的峰值点：

可根据递推结果x获得拟合函数的各项系数x_k(k＝0,1,2,3)，所述拟合函数的导数为

s′(n)＝x₁+2x₂n+3x₃n²， (9)

令s′(n)＝0，求得峰值点的位置

$N_{max}=\frac{-2x_2-\sqrt{4x_2^2-6x_1{x}_3}}{6x_3},---\left(10\right)$

求得峰值点的大小

P＝s(N_max)， (11)

步骤九，对计算得的峰值点进行判断：

当N≥20，对N-19，N-18，…，N这20组结果设定如下判断条件

$\begin{array}{c}\frac{1}{20}\sqrt{\underset{i=N-19}{\overset{N}{\Σ}}{\[N_{\max }\left(i\right)-{\stackrel{\‾}{N}}_{\max }\]}^2}\≤\frac{1}{1000}{N}_{\max }\left(N\right)\\ \frac{1}{20}\sqrt{\underset{i=N-19}{\overset{N}{\Σ}}{\[P\left(i\right)-\stackrel{\‾}{P}\]}^2}\≤\frac{1}{1000}P\left(N\right)\end{array},---\left(12\right)$

其中，N_max(i)，P(i)表示计算的峰值点位置与峰值点大小，为20组结果的平均值，若符合条件，则拟合函数的峰值点(N_max(N)，P(N))即为谐波信号的峰值点，若不符合条件，则返回步骤七。在测量中，只需谐波信号的一部分采样点就可以计算得准确的峰值点信息，参与计算部分如图3中矩形区域内所示。当结果满足条件时，拟合函数与谐波信号重合，如图4所示。

步骤十，计算气体浓度：在温度和压力已知的情况下，气体浓度与谐波峰值成正比，根据激光吸收光谱技术中峰值点信息计算出CO₂气体浓度，测量结果如图5所示；

实验结果表明，使用部分采样点就可以计算得准确的峰值信息，在保证了测量准确性的同时提高测量速度。

以上对本发明及其实施方式的描述，并不局限于此，附图中所示仅是本发明的实施方式之一。在不脱离本发明创造宗旨的情况下，不经创造地设计出与该技术方案类似的结构或实施例，均属本发明保护范围。

Claims (1)

1.一种基于递推最小二乘法的激光吸收光谱谐波信号峰值点提取方法，该方法包括以下步骤：

步骤一，以激光光谱吸收信号的谐波信号采样点s(n)和采样点在一个周期中的序列数n为输入，写出对谐波信号的主峰部分进行拟合的三次多项式：

在白噪声影响下，所述三次多项拟合函数可表示为

s(n)＝u(n)·x+e(n)， (1)

其中，e(n)为白噪声，

$\begin{array}{c}u\left(n\right)=\left[\begin{array}{cccc}1& n& n^2& n^3\end{array}\right],\\ x={\left[\begin{array}{cccc}{x}_0& x_1& x_2& x_3\end{array}\right]}^T,\end{array}---\left(2\right)$

将所述N(N≥1)个相邻输入s(n)和n表示为式(3)(4)的矩阵形式

$S_N={\left[\begin{array}{cccc}s\left(1\right)& s\left(2\right)& ...& s\left(N\right)\end{array}\right]}^T,---\left(3\right)$

$U_N=\left[\begin{array}{c}u\left(1\right)\\ u\left(2\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ u\left(N\right)\end{array}\right]=\left[\begin{array}{cccc}1& 1^1& 1^2& 1^3\\ 1& 2^1& 2^2& 2^3\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 1& N& N^2& N^3\end{array}\right],---\left(4\right)$

则待求参数x_N可以通过式(5)计算出

$x_N=J_N{U}_N^T{S}_N,---\left(5\right)$

其中

$J_N={\left(U_N^T{U}_N\right)}^{-1},---\left(6\right)$

步骤二，计算递推初始值：

使用第一组输入数据(N＝1)计算递推初始值，

$\left\{\begin{array}{c}{J}_1={\left(U_1^T{U}_1\right)}^{-1}\\ x_1=J_1{U}_1^T{S}_1\end{array}\right.,---\left(7\right)$

步骤三，增加测量数据代入递推公式进行计算：

所述递推公式为：

$\left\{\begin{array}{c}{J}_{N+1}=J_N-J_N{u}^T\left(N\right){\[I+u\left(N\right)J_N{u}^T\left(N\right)\]}^{-1}u\left(N\right)J_N\\ x_{N+1}=x_N+J_{N+1}{u}^T\left(N\right)\[s\left(N\right)-u\left(N\right)x_N\]\end{array}\right.,---\left(8\right)$

步骤四，使用递推结果计算拟合函数的峰值点：

可根据递推结果x获得拟合函数的各项系数x_k(k＝0,1,2,3)，所述拟合函数的导数为

s′(n)＝x₁+2x₂n+3x₃n²， (9)

令s′(n)＝0，求得峰值点的位置

$N_{max}=\frac{-2x_2-\sqrt{4x_2^2-6x_1{x}_3}}{6x_3},---\left(10\right)$

求得峰值点的大小

P＝s(N_max)， (11)

步骤五，对计算得的峰值点进行判断：

当N≥20，对N-19，N-18，…，N这20组结果设定如下判断条件

$\begin{array}{c}\frac{1}{20}\sqrt{\underset{i=N-19}{\overset{N}{\mathrm{\Σ}}}{\[N_{\max }\left(i\right)-{\stackrel{\‾}{N}}_{\max }\]}^2}\≤\frac{1}{1000}{N}_{\max }\left(N\right)\\ \frac{1}{20}\sqrt{\underset{i=N-19}{\overset{N}{\mathrm{\Σ}}}{\[P\left(i\right)-\stackrel{\‾}{P}\]}^2}\≤\frac{1}{1000}P\left(N\right)\end{array},---\left(12\right)$

其中，N_max(i)，P(i)表示计算的峰值点位置与峰值点大小， 为20组结果的平均值，

若符合条件，则拟合函数的峰值点(N_max(N)，P(N))即为谐波信号的峰值点，若不符合条件，则返回步骤三。