CN106053386A - 一种用于吸收光谱测量的无近似Voigt线型拟合方法 - Google Patents

Abstract

本发明针对可调谐半导体激光吸收光谱技术中的线型拟合的关键环节，提供了一种没有任何近似步骤的Voigt线型拟合方法，从根本上消除了近似运算对线型拟合以及吸收光谱测量精度的影响。该方法包括Voigt线型计算步骤、Voigt线型函数偏导数计算步骤、最小二乘拟合步骤。第一步，利用初始化光谱参数和FFT方法计算Voigt线型函数；第二步，将Voigt线型函数对待拟合参数的偏导数转换为Gauss或Lorentz线型函数对该参数的偏导数，然后再利用FFT方法运算获得Voigt线型函数对该参数的偏导数；第三步，利用实际测量的Voigt线型曲线和最小二乘拟合算法进行Voigt线型拟合，根据优值函数判断最小二乘拟合运算结束，还是进入下一次迭代运算。

Description

一种用于吸收光谱测量的无近似Voigt线型拟合方法

技术领域

本发明属于光学流场诊断领域，涉及半导体激光吸收光谱技术和信号处理技术，可用于提高激光吸收光谱进行温度、组分浓度、压强测量的精度。

背景技术

可调谐二极管激光吸收光谱技术(TDLAS)是目前进行燃烧流场温度、组分浓度、速度或压强等参数测量的非常有效的方法。它的基本原理是Beer-Lambert定律，利用气体分子对激光光谱吸收特性，测量穿过吸收气体的激光光强，可以得到具有特定吸收线型光强吸收曲线。吸收线型函数与场的参数直接相关，通过对测量吸收曲线的线型拟合就可以得到这些参数的数值。

对于气体吸收环境，吸收线型主要由两种物理机制决定，分子热运动引起的多普勒展宽和分子间碰撞引起的碰撞展宽。这两种展宽分别产生两种相对应的线型函数，Gauss线型和Lorentz线型。当气体压力较小时，多普勒展宽占主导地位，吸收线型用Gauss线型描述；而在高压条件下，分子间碰撞加强，碰撞展宽占主导，吸收线型用Lorentz线型。而实际上，在大多数情况下两种展宽特性的程度相差不大，这时需要用Voigt线型来描述吸收线型。Voigt线型表述为Gauss线型和Lorentz线型的卷积。这个卷积表达式没有明确的解析形式。这为Voigt线型函数的应用带来了两个麻烦：第一是采用积分形式进行Voigt线型的数值计算需要消耗大量的计算时间；第二是积分形式的Voigt线型函数无法使用现有的线型拟合工具进行曲线拟合，因为一般的线型拟合工具都要求一个解析表达式作为输入参数。

目前对于解决以上两个问题的通用方法是对Voigt线型进行解析表达式近似，用来代替其原有的积分表达式。然而近似解析表达式也有其固有的缺陷——近似必然会引入误差。若希望减小近似误差就必须采用复杂程度很高近似解析表达式，甚至是复数表达式。复杂的近似表达式又会引起较长的计算时间，以及在线型拟合中不收敛情况。因此对于吸收光谱测量而言，一种能够快速计算的无近似的Voigt线型拟合方法具有非常重要的意义。

下述文献报道涉及激光吸收光谱测量中的Voigt线型函数的计算和拟合相关内容。

1.美国耶鲁大学S.M.Abrarov等人在论文“Rapidly convergent series forhigh-accuracy calculation of the Voigt function”(Journal of QuantitativeSpectroscopy&Radiative Transfer 111(2010)372–375)提出用一种基于傅里叶展开的指数函数序列来实现对Voigt函数的高精度近似。这种近似方法的计算准确度能够到达10^-9。虽然其计算速度要明显快于积分方法，但近似表达式是一个序列之和的形式，非常复杂。

2.加拿大皇后大学A.B.McLean在论文“Implementation of an efficientanalytical approximation to the Voigt function for photoemission lineshapeanalysis”(Journal of Electron Spectroscopy and Related Phenomena 64(1994)125-132)提出了一种表达式较简单的Voigt函数近似表达式。该表达式由4个形式相同的多项式之和构成，每个多项式包含4个数值固定的参数值。该近似表达式形式简单，计算速度很快，也适用于线型拟合，但近似的准确度不高。

3.中科院安徽光学精密机械研究所张庆礼等人在论文“X射线衍射双峰的Voigt峰形函数拟合”(人工晶体学报,2009年38卷第2期)中适用M个节点的高斯-厄米数值积分公式对Voigt线型函数进行近似，并利用该近似表达式对双峰的Voigt线型曲线进行拟合。拟合结果表明其收敛速度和稳定性都较好，但文章同时指出节点的增多同时会增加计算量，在多次数值迭代情形下回明显见面计算速度。

4.瑞典Umea大学Jonas Westberg在论文“Fast and non-approximatemethodology for calculation of wavelength-modulated Voigt lineshape functionssuitable for real-time curve fitting”(JJournal of Quantitative Spectroscopy&Radiative Transfer 113(2012)2049-2057)提供了一种针对激光吸收光谱测量中波长调制Voigt线型函数的快速无近似计算方法。该方法的应用目标实时曲线拟合。文中利用快速傅里叶变换方法计算与卷积相关的表达式，详细介绍了如何求解波长调制Voigt线型函数进行傅里叶展宽之后的n阶调制谐波系数。整个计算过程没有任何表达式近似，而且计算速度较快。但文章没有介绍如何利用该方法进行Voigt线型拟合。

上述文献证明了精确计算Voigt线型函数的复杂性，以及无近似计算的可行性，但关于Voigt线型函数的拟合的研究仍然是基于近似表达式。无近似Voigt线型函数的计算为无近似Voigt线型拟合提供了可能性。本发明在此基础上实现一种用于吸收光谱测量的无近似Voigt线型拟合方法。

发明内容

本发明的目的是为了提高TDLAS方法进行气体参数测量的精度，提供一种无近似Voigt线型拟合方法。该方法以非线性最小二乘拟合方法为基础，使用FFT方法计算Voigt线型函数以及该函数对待拟合参数的一阶偏导数，将计算结果应用到最小二乘拟合算法内部，实现整个拟合过程的无近似。该方法从理论上解决了长期以来Voigt线型的近似计算所引起的各类问题。该方法适用于所有可以利用TDLAS进行流场参数测量的领域，如发动机羽流、超燃冲压发动机隔离段/燃烧室和内燃机等温度和组分浓度测量，以及火力发电、燃煤锅炉等领域的燃烧效率优化、污染物监控等。

本发明描述了无近似Voigt线型函数计算方法，实现方案如下：

(1)根据所给出的中心波长、Gauss线宽、Lorentz线宽等线型参数，分别计算归一化的Gauss线型函数和Lorentz线型函数；

(2)分别计算Gauss线型函数和Lorentz线型函数的傅里叶变换；

(3)利用卷积的傅里叶变换特性，得到Voigt线型函数的傅里叶变换等于Gauss线型函数和Lorentz线型函数的傅里叶变换之积；

(3)利用逆傅里叶变换得到归一化的Voigt线型函数。

本发明提供了无近似地计算Voigt线型函数线型参数的一阶偏导数的方法，实现方案如下：

(1)利用卷积的微分特性，将Voigt线型函数对线型参数的一阶偏导数转换为Gauss线型函数或Lorentz线型函数对这些线型参数的一阶偏导数的形式；

(2)计算Gauss线型函数或Lorentz线型函数对这些线型参数的一阶偏导数；

(3)利用卷积的傅里叶变换性质，计算Voigt线型函数对这些线型参数的一阶偏导数；

本发明改进了非线性最小二乘拟合算法，使该算法适用于没有明确的解析表达式的曲线拟合，实现方案如下：

(1)设置待拟合线型参数的初始值，并利用该初始值计算Voigt线型函数和及其一阶偏导数；

(2)利用计算的Voigt线型函数和实测的线型曲线计算优值因子SSE；

(3)利用计算的Voigt线型函数的一阶偏导数构建黑塞矩阵和梯度方程，求解待拟合参数的增量；

(4)利用增量之后的拟合参数重新计算Voigt线型函数和优值因子SSE。

(5)比较新的优值因子与之前的优值因子的差异，判断是结束运算还是重复步骤(1)-(4)直到拟合差异收敛到一个可接受的较小值。

附图说明

图1为实现本发明所述方法的原理流程图。

图2为利用本发明所述方法计算的Voigt线型函数及其对线型参数的一阶偏导数。

图3为本发明所述方法对仿真数据拟合的结果。

图4为本发明所述方法对实测数据拟合的结果。

图5为本发明所述方法计算消耗时间随数据点数增加的变化曲线。

具体实施方式

现结合附图和实施例对本发明用于吸收光谱测量的无近似Voigt线型拟合方法做进一步详细描述。图1给出了本发明的无近似Voigt线型拟合方法的原理流程图。

本发明的实现方案如下：本发明充分利用卷积运算的特性以及Voigt线型函数的特征，结合实际激光吸收光谱测量的应用条件，提出了无近似Voigt线型拟合方法。发明分为三个步骤，即无近似Voigt线型函数计算步骤、无近似Voigt线型函数一阶参数偏导数计算步骤和Voigt线型的无近似最小二乘拟合步骤。具体实施方案如下：

步骤一、无近似Voigt线型函数计算

Voigt线型函数在形式上Gauss线型函数与Lorentz线型函数的卷积，共包括4个线型参数，中心波长ν₀，Gauss线宽Δν_G，Lorentz线宽Δν_L和线型幅度K。Voigt线宽可以通过Gauss线宽和Lorentz线宽求解得到。在计算Voigt线型函数时，这4个线型参数是已知量。具体计算步骤如下：

1)计算归一化的Gauss线型函数φ_G(ν)和归一化的Lorentz线型函数φ_L(ν)，这两个线型函数都有明确的解析表达式，表示如下。

${\mathrm{\φ}}_G\left(v\right)=\frac{2}{{\mathrm{\Δv}}_G}\sqrt{\frac{ln2}{\mathrm{\π}}}\exp \[-4ln2{\left(\frac{v}{{\mathrm{\Δv}}_G}\right)}^2\]---\left(1\right)$

${\mathrm{\φ}}_L\left(v\right)=\frac{2}{\mathrm{\π}}\frac{{\mathrm{\Δv}}_L}{4v^2+{\mathrm{\Δv}}_L^2}---\left(2\right)$

相应地，峰值位于中心波长ν₀处的Voigt线型函数的表达式表示为：

${\mathrm{\φ}}_V(v-v_0)=K\·{\mathrm{\φ}}_G(v-v_0)\⊗{\mathrm{\φ}}_L\left(v\right)---\left(3\right)$

2)计算归一化Gauss线型函数和归一化Lorentz线型函数的傅里叶变换,H_G(ω)和H_L(ω)；

3)利用卷积的傅里叶变换特性，得到Voigt线型函数的傅里叶变换H_V(ω)等于Gauss线型函数和Lorentz线型函数的傅里叶变换之积，即

H_V(ω)＝H_G(ω)H_L(ω) (4)

4)进行傅里叶逆变换得到归一化的Voigt线型函数，再乘以线型幅度A就得到实际的Voigt线型函数φ_V(ν)。

步骤二、无近似Voigt线型函数一阶参数偏导数

Voigt线型函数有4个线型参数，因此它对线型参数的一阶偏导数也有4个。其中对线型幅度的一阶偏导数就是归一化的Voigt线型函数，非常简单，下面主要介绍其他三个参数的一阶偏导数的具体计算步骤：

1)以参数中心波长ν₀为例，从方程(3)可以得到，Voigt线型函数对Gauss线宽的一阶偏导数表示为：

$\frac{\∂{\mathrm{\φ}}_V(v-v_0)}{\∂v_0}=K\frac{\∂}{\∂v_0}{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}{\mathrm{\φ}}_G(u-v_0)\·{\mathrm{\φ}}_L(v-u)du---\left(5\right)$

利用卷积的微分特性，方程(5)可表示为：

$\begin{array}{c}\frac{\∂{\mathrm{\φ}}_V(v-v_0)}{\∂v_0}=K{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}\frac{\∂{\mathrm{\φ}}_G(u-v_0)}{\∂v_0}\·{\mathrm{\φ}}_L(v-u)du\\ =K\·\frac{\∂{\mathrm{\φ}}_G(v-v_0)}{\∂v_0}\⊗{\mathrm{\φ}}_L\left(v\right)\end{array}---\left(6\right)$

方程(6)仍然是一个卷积形式，这个卷积包含两个部分，第一部分是Gauss线型函数中心波长ν₀的一阶偏导数，第二部分是Lorentz线型函数。由于Gauss线型函数有解析表达式方程(1)，因此该一阶偏导数可以用解析表达式表示：

$\frac{\∂{\mathrm{\φ}}_G(v-v_0)}{\∂v_0}=8ln2\frac{v-v_0}{{\mathrm{\Δv}}_G^2}{\mathrm{\φ}}_G(v-v_0)---\left(7\right)$

因此方程(6)表示的Voigt线型函数对ν₀的一阶偏导数是两个解析表达式的卷积，可以利用计算Voigt函数类似的傅里叶变换方法进行计算。

2)根据方程(7)计算Gauss线型函数对ν₀的一阶偏导数

3)计算Gauss线型函数对ν₀一阶偏导数的傅里叶变换；

4)利用卷积的傅里叶变换特性得到Voigt线型函数对中心波长的一阶偏导数的傅里叶变换值；

5)使用逆傅里叶变换计算得到Voigt线型函数ν₀的一阶偏导数

6)使用(1)-(5)同样的步骤可以计算Voigt线型函数对Gauss线宽和Lorentz线宽的一阶偏导数，和二者的卷积表达式如下。

$\frac{\∂{\mathrm{\φ}}_V(v-v_0)}{\∂{\mathrm{\Δv}}_G}=K\·\frac{\∂{\mathrm{\φ}}_G(v-v_0)}{\∂{\mathrm{\Δv}}_G}\⊗{\mathrm{\φ}}_L\left(v\right)---\left(8\right)$

$\frac{\∂{\mathrm{\φ}}_V(v-v_0)}{\∂{\mathrm{\Δv}}_L}=K\·{\mathrm{\φ}}_G(v-v_0)\⊗\frac{\∂{\mathrm{\φ}}_L\left(v\right)}{\∂{\mathrm{\Δv}}_L}---\left(9\right)$

图2是利用无近似计算方法所画出的在激光吸收光谱测量中比较典型的Voigt线型曲线以及其对每个线型参数的一阶偏导数。

步骤三、Voigt线型的无近似最小二乘拟合

Voigt线型的无近似最小二乘拟合是针对Voigt线型函数对经典的Levenberg-Marquardt线型拟合方法的一种改进。最小二乘拟合是一种迭代算法，每一迭代周期内的具体步骤如下：

1)设置待拟合线型参数向量p＝[K,ν₀,Δν_G,Δν_L]的初始值，并利用该初始值计算Voigt线型函数及其对每个线型参数的偏导数；

2)利用实测的线型曲线和计算的Voigt线型函数来计算优值因子χ²(p)，用公式表示为

${\mathrm{\χ}}^2\left(p\right)=\underset{i=0}{\overset{N-1}{\Σ}}{\[\frac{{\mathrm{\φ}}_{V,i}-{\mathrm{\φ}}_V(v_i;p)}{{\mathrm{\σ}}_i}\]}^2---\left(10\right)$

其中N为实测数据点数；φ_V,i表示第i个测量数据；σ_i表示第i个测量数据的方差，在不知道具体数值的的情况下可以将其设置为1。

3)计算优值因子χ²(p)对待拟合参数的一阶偏导数；

$\frac{\∂{\mathrm{\χ}}^2\left(p\right)}{\∂p_j}=2{\mathrm{\β}}_j=-2\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\frac{\[{\mathrm{\φ}}_{V,i}-{\mathrm{\φ}}_V(v_i;p)\]}{{\mathrm{\σ}}_i^2}\frac{\∂{\mathrm{\φ}}_V(v_i;p)}{p_j}---\left(11\right)$

(11)表示的是优值因子对第i个拟合参数的偏导数，为了后续描述的方便，用β_j表示该式中的求和项。

4)利用计算的Voigt线型函数对每个线型参数的偏导数计算并构建Hessian矩阵，该矩阵的元素α_j,k的计算公式表示为：

${\mathrm{\α}}_{j,k}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\frac{1}{{\mathrm{\σ}}_i^2}\[\frac{\∂{\mathrm{\φ}}_V(v_i;p)}{p_j}\frac{\∂{\mathrm{\φ}}_V(v_i;p)}{p_k}\]---\left(11\right)$

5)利用一个小的最小二乘算法因子λ给Hessian矩阵的对角线元素施加一个的偏置，形成一个新的矩阵，这个矩阵的元素表示为：

$\begin{array}{c}{\mathrm{\α}}_{j,j}^{\′}\≡{\mathrm{\α}}_{j,j}(1+\mathrm{\λ})\\ {\mathrm{\α}}_{j,k}^{\′}\≡{\mathrm{\α}}_{j,k}(j\≠k)\end{array}---\left(12\right)$

6)求解下列线型方程组得到待拟合参数的增量Δp；

$\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{\α}}_{j,k}^{\′}{\mathrm{\δp}}_k={\mathrm{\β}}_j---\left(13\right)$

7)待拟合线型参数加上求解得到的增量，利用新的线型参数重复步骤(1)和(2)计算新的优值因子χ²(p+Δp)；

8)比较新优值因子χ²(p+Δp)和旧优值因子χ²(p)判断是否结束迭代循环。如果χ²(p+Δp)>χ²(p)，说明待拟合参数增量使计算Voigt线型和实测的线型曲线之间的SSE更差，需要增大算法因子λ；如果χ²(p+Δp)<χ²(p)，首先判断两者差值χ²(p)-χ²(p+Δp)是否小于收敛阈值ε，如果χ²(p)-χ²(p+Δp)<ε,则结束迭代循环，否则，用p+Δp代替p，并减小算法因子λ，重复迭代过程。

图3给出的是利用本发明提出的无近似Voigt线型函数拟合方法，对一个信噪比为6dB的仿真吸收光谱数据进行Voigt线型拟合的结果。上图是拟合曲线和仿真数据，下图为拟合误差，可以看出拟合误差在3％以内。

图4给出了利用本发明提出的无近似Voigt线型函数拟合方法，对实际吸收光谱测量数据的Voigt线型拟合结果。

图5给出了本发明提出的无近似Voigt线型函数拟合算法所消耗时间随拟合数据点数的变化关系，可以看出对于1000点的数据，整个拟合过程只需要10ms。

Claims (1)

1.一种用于吸收光谱测量的无近似Voigt线型拟合方法，包含以下步骤：

步骤一、无近似Voigt线型函数计算

(1)根据所给出的中心波长、Gauss线宽、Lorentz线宽等线型参数，分别计算归一化的Gauss线型函数和Lorentz线型函数；

(2)分别计算Gauss线型函数和Lorentz线型函数的傅里叶变换；

(3)利用卷积的傅里叶变换特性，得到Voigt线型函数的傅里叶变换等于Gauss线型函数和Lorentz线型函数的傅里叶变换之积；

(3)利用逆傅里叶变换得到归一化的Voigt线型函数；

步骤二、无近似Voigt线型函数一阶参数偏导数

(1)利用卷积的微分特性，将Voigt线型函数对线型参数的一阶偏导数转换为Gauss线型函数或Lorentz线型函数对这些线型参数的一阶偏导数的形式；

(2)计算Gauss线型函数或Lorentz线型函数对这些线型参数的一阶偏导数；

(3)利用卷积的傅里叶变换性质，计算Voigt线型函数对这些线型参数的一阶偏导数；

步骤三、Voigt线型的无近似最小二乘拟合

(1)设置待拟合线型参数的初始值，并利用该初始值计算Voigt线型函数和及其一阶偏导数；

(2)利用计算的Voigt线型函数和实测的线型曲线计算优值因子SSE；

(3)利用计算的Voigt线型函数的一阶偏导数构建黑塞矩阵和梯度方程，求解待拟合参数的增量；

(4)利用增量之后的拟合参数重新计算Voigt线型函数和优值因子SSE。

(5)比较新的优值因子与之前的优值因子的差异，判断是结束运算还是重复步骤(1)-(4)直到拟合差异收敛到一个可接受的较小值。