一种基于吸收光谱图重构的多组分气体浓度定量分析方法

摘要

本发明公开了一种基于光谱图重构的多组分气体定量分析方法。首先预估待分析气体中各组分的浓度和折算吸光率，根据Lambert-Beer定理重构光谱图。然后以重构光谱图与实际光谱图的谱线值之差作为输入，构建一个勾型损失函数。输入以零为界，输入小于零时，函数值递减，减幅较小，输入大于零时，函数递增，增幅较大。然后采用动量法或遗传算法对待分析气体中各组分浓度进行迭代运算，直到相邻两次迭代结果中各组分浓度差小于预设的容许值，结束迭代运算。本发明公开的方法分析速度快、结果准确、干扰气对分析结果影响小，可用于气录井、产品质量检查与故障诊断，以及五金、化工、环保等领域的多组分气体定量分析中。

说明书

技术领域

本发明涉及光谱定量分析领域，涉及一种多组分气体浓度定量分析的方法，特别涉及一种基于光谱图重构的多组分气体定量分析方法。

背景技术

光谱分析常用于混合物的成分分析，也是多组分气体在线定量分析的重要方法之一。气体分析涉及到科学研究、环境保护、设备的故障诊断、产品质量的检测、天然气与石油勘探录井等多个领域，在这些领域中，组分气体的定量分析发挥着非常重要的作用。如：变压器的故障诊断方法中，有一种就是通过检测变压器油中乙炔等气体的含量来识别的，天然气的燃烧值是通过检测天然气中甲烷、乙烷等组分的含量来确定的。虽然气体吸收光谱的光谱分析是气体定量分析的一种重要方法，但这种方法还存在一些没有解决的问题，其中之一就是由于被分析背景比较复杂，除了目标气之外，可能还存在某些未知的干扰气体。例如，家具以及五金器材等的质量检测指标中，有一项是单位时间内甲醛的挥发量，但家具以及五金器材等产品所挥发的气体中，除了甲醛外，可能还有乙醛等其它有机化合物，而且其中有一些成分的吸收光谱可能与甲醛吸收光谱有交叠。于是，如何减小这些干扰气的存在对目标气体分析的影响，就成了多组分气体光谱定量分析的重要内容。另外，在科学研究中，如化学反应，也可以通过光谱在线分析来识别是否有新物质的产生，或者中间过程是否有新的物质成分。

发明内容

本发明的目的在于，提供一种基于吸收光谱图重构的多组分气体定量分析方法，该方法在Lambert-Beer定理基础上，通过构建一个称之为勾型损失函数的函数，并采用动量法神经网络学习算法来对各种目标气的浓度进行学习，最终达到减小未知气体对目标气定量分析带来的干扰的目的。

为了实现上述任务，本发明采用如下技术解决方案：

一种基于吸收光谱图重构的多组分气体浓度定量分析方法，首先预估各组分目标气体的浓度，估计各组分气体在光谱所感兴趣的波数段范围内的折算吸光率，获得初始重构光谱图；选择所感兴趣的波数段范围内的重构光谱图与实际扫描得到的光谱图之差为重构偏差，以这个偏差作为输入，构建勾型函数作为光谱图重构的损失函数，并采用神经网络学习算法中的动量法对重构光谱进行迭代运算，以更新每种目标气体浓度值，直到相邻两次迭代结果中，所有目标气体浓度值之间的偏差小于预定的容许偏差；

所述勾型函数具有如下式形式：

$f\left({ϵ}_i\right)=k_1{ϵ}_i^2+{(\exp \left(k_2{ϵ}_i\right)-1)}^2$

式中ε_i＝u_i-v_i，表示第i条谱线实际光谱值u_i与重构光谱值v_i之差；k₁、k₂称之为勾型损失函数的形状系数。

各组分气体的折算吸光率是根据其当前浓度估计值或预设值，以及相应单组分气体样本气体光谱通过曲线拟合法来估计的；曲线拟合法既是分段线性化、多项式逼近法或样条插值法；

如果光谱仪输出光谱图是以透光率作为输出的，则对于单组份气体A_j的扫描光谱图的第i个光谱值v_i的自然对数除以单组份气体A_j的浓度C_Aj，即为目标气A_j的折算吸光率：

δ_ij(k)＝log(v_i(k))/C_Aj(k)，k＝1，2，…，S

式中log(·)表示自然对数运算，S为单组份气体A_j的样本数量；

如果光谱仪输出光谱图是以吸光率作为输出的，则对于单组份气体A_j的扫描光谱图的第i个光谱值v_i除以单组份气体A_j的浓度，即为目标气A_j的折算吸光率：

δ_ij(k)＝v_i(k)/C_Aj(k)，k＝1，2，…，S

于是，如果单组份气体A_j的估算浓度为C_Aj，若采用分段线性化估算δ_ij，有：

δ_ij＝(δ_ij(k-1)-δ_ij(k))(C_Aj-C_Aj(k))/(C_Aj(k+1)-C_Aj(k))

式中C_Aj大于C_Aj(k)，且小于C_Aj(k+1)。

若采用多项式逼近法估算δ_ij，有：

${\delta }_{\mathrm{ij}}=\underset{s=0}{\overset{Q}{\Sigma }}{b}_s{C}_{\mathrm{Aj}}^T$

式中Q≤S-1，系数b_s(s＝1，2，…，Q)由下式确定：

B＝inv(H^TH)H^TC

式中B＝[b₀，b₁，…，b_Q]^T，C＝[C_Aj(1)，C_Aj(2)，…，C_Aj(S)]^T，

$H=\left(\begin{array}{cccc}1& C_{\mathrm{Aj}}\left(1\right)& L& C_{\mathrm{Aj}}{\left(1\right)}^Q\\ 1& C_{\mathrm{Aj}}\left(2\right)& L& C_{\mathrm{Aj}}{\left(2\right)}^Q\\ M& M& O& M\\ 1& C_{\mathrm{Aj}}\left(S\right)& L& C_{\mathrm{Aj}}{\left(S\right)}^Q\end{array}\right).$

所述迭代运算方法采用神经网络的学习算法中的动量法、最速下降法、遗传算法或模拟退火法来实现，并对每种气体浓度的增量设定了上限和下限，以避免学习收敛过程中的来回跳变，增快收敛速度；具体包括下列步骤：

(1)预估N组分气体的目标气A₁、A₂、…、A_N的浓度分别为：C_A1、C_A2、…、C_AN；预估的方法是预设浓度值法、特征谱线的最小二乘法、主分量回归法、加权最小二乘估计法或反最小二乘估计法；

令光谱分析的M组标定样本中Y＝[C₁ C₂ …C_N]表示N组分气体的目标气的期望浓度矩阵，E＝[E₁ E₂ …E_N]表示N组分气体的残差矩阵，H＝[V₁ V₂ … V_L]表示所感兴趣的L条谱线值形成的矩阵，其中C_i＝[C_i1 C_i2 …C_NM]^T，E_i＝[E_i1 E_i2 …E_NM]^T，i＝1，2，…，N，V_j＝[V_j1 V_j2 …V_jM]^T，j＝1，2，…，L。

如果光谱仪输出光谱图是以透光率作为输出的，则最小二乘法、加权最小二乘估计法可按(1a)式计算：

[C_A1，C_A2，…，C_AN]＝log(V)W         (1a)

对于最小二乘法，式中W＝inv(Q^TQ)Q^TY，对于加权最小二乘法，式中W＝inv(Q^TΛQ)Q^TΛY，其中inv(·)表示求逆运算，Q_ij＝log(H_ij)，Λ是一个L×L的加权矩阵；

如果光谱仪输出光谱图是以吸光率作为输出的，则最小二乘法、加权最小二乘估计法可按(1a)式计算：

[C_A1，C_A2，…，C_AN]＝VW    (1c)

对于最小二乘法，式中W＝inv(H^TH)H^TY，对于加权最小二乘法，式中W＝inv(H^TΛH)H^TΛY；

(2)按权利要求3，根据各组分的浓度，计算各组分折算吸光率δ_ij。δ_ij表示在所感兴趣光谱段中第i谱线的光谱值对应的目标气A_j的折算吸光率，它是实际吸光率与光程的乘积；

如果光谱仪输出光谱图是以透光率作为输出的，则对于单组份气体A_j的扫描光谱图的第i个光谱值v_i的自然对数除以单组份气体A_j的浓度C_Aj，即为目标气A_j的折算吸光率：

δ_ij＝log(v_i)/C_Aj    (2a)

式中log(·)表示自然对数运算；

如果光谱仪输出光谱图是以吸光率作为输出的，则对于单组份气体A_j的扫描光谱图的第i个光谱值v_i除以单组份气体A_j的浓度，即为目标气A_j的折算吸光率：

δ_ij＝v_i/C_Aj   (2b)

(3)在光谱图中所感兴趣的光谱段内，获得初始重构吸收光谱图：

如果光谱仪输出光谱图是以透光率作为输出的，则重构光谱图中第i个光谱值为：

$v_i=\exp (-\underset{j=1}{\overset{N}{\Sigma }}{\delta }_{\mathrm{ij}}{C}_{\mathrm{Aj}})---\left(3a\right)$

如果光谱仪输出光谱图是以吸光率作为输出的，则重构光谱图中第i个光谱值为：

$v_i=\underset{j=1}{\overset{N}{\Sigma }}{\delta }_{\mathrm{ij}}{C}_{\mathrm{Aj}}---\left(3b\right)$

式中exp(·)表示自然指数运算，N表示目标气体的种数；

(4)如果光谱仪输出光谱图是以透光率作为输出的，则令重构偏差为：

ε_i＝u_i-v_i    (4a)

如果光谱仪输出光谱图是以吸光率作为输出的，则令重构偏差为：

ε_i＝(u_i-v_i)exp(-v_i)    (4b)

式中u_i表示光谱扫描得到的吸收光谱图中第i个光谱值。选择勾型函数作为损失函数：

$f(C_{A1},C_{A2},\Lambda ,C_{\mathrm{AN}})=\underset{i}{\overset{M}{\Sigma }}(k_1{ϵ}_i^2+{(\exp \left(k_2{ϵ}_i\right)-1)}^2)---\left(5\right)$

式中M表示所重构的吸收光谱段的光谱值个数，k₁、k₂称之为勾型损失函数的形状系数；

(5)用神经网络学习方法中的动量法对目标气体的浓度进行迭代运算：

C_Aj(t)＝C_Aj(t-1)-α(j)[(1-λ)ΔC_Aj(t-1)+γΔC_Aj(t-2)]      (6)

式中α(j)为气体A_j浓度定量分析的学习率，λ为动量因子，且0≤λ≤1，t为学习迭代次数。对于ΔC_Aj，可采用式(7)来估计：

${\mathrm{\Delta C}}_{\mathrm{Aj}}=\frac{\partial f\left(E\right)}{{\partial C}_{\mathrm{Aj}}}\approx \frac{f(C_1,C_2,\Lambda ,C_{\mathrm{Aj}}+e,\Lambda ,C_{\mathrm{AN}})}{e}---\left(7\right)$

式中e是一个足够小的量，可取e＝0.01；

如果ΔC_Aj＞Max_Aj，则令ΔC_Aj＝Max_Aj；如果ΔC_Aj＜Min_Aj，则令ΔC_Aj＝Min_Aj，其中Max_Aj和Min_Aj分别表示所设定的气体组分A_j浓度的增量上限和下限，且Max_Aj＞0，Min_Aj＜0；

(6)用迭代结果取代(1)中的预估值，依次重复上述步骤(2)、(3)、(4)、(5)，直到每种目标气的浓度与上一次迭代结果之差小于设定的容许偏差。

该方法首先预估各组分目标气体的浓度，估计各组分气体在不同浓度情况下在不同光谱数据点的折算吸光率，即各组分气体在不同光谱数据点的实际吸光率与光程的乘积，然后根据Lambert-Beer定理，获得初始重构光谱图，选择重构光谱图与实际扫描得到的光谱图之差为重构偏差，构建一种称之为勾型函数作为光谱图重构的损失函数，并采用神经网络学习算法中的动量法、最速下降法、遗传算法、模拟退火算法等对重构光谱进行迭代运算，以更新每种目标气体浓度值，直到各组分气体浓度收敛到相邻两次迭代运算结果变化小于所设定的容许偏差。

本发明实现的多组分气体浓度定量分析速度快、结果准确，特别有助于减小待分析混合气中存在的干扰气对目标组分分析结果的影响，可用于石油、天然气勘探的气录井，产品质量检查与故障诊断、五金、化工、环境保护等领域的多组分气体光谱定量分析应用中。

附图说明

图1(a)是甲烷、乙烷、丙烷、异丁烷、正丁烷、异戊烷和正戊烷七组分气体浓度分别为6.09％、1.58％、0.935％、0.444％、0.477％、0.233％和0.234％标准混合气400-1400波数段傅立叶变换中红外吸收光谱图；

图1(b)是一种混合油气的傅立叶变换中红外实际吸收谱图的400-1400波数段。该混合油气的气相色谱分析结果为：甲烷：1.0652％、乙烷：0.0253％、丙烷：0.0199％、异丁烷：0.0105％、正丁烷：0.0113％、异戊烷：0.0041％和正戊烷0.0031％；

图2是浓度均为1％的甲烷、乙烷、丙烷、异丁烷、正丁烷、异戊烷和正戊烷七组分气体的400-1400波数段傅立叶变换中红外吸收光谱图；

图3是勾型损失函数的曲线图；

图4(a)是图1(a)所示光谱图及相应的逼近光谱，光谱逼近分析结果分别为：甲烷：6.02％、乙烷：1.62％、丙烷：0.95％、异丁烷：0.46％、正丁烷：0.48％、异戊烷：0.24％、正戊烷：0.24％。

图4(b)图1(b)所示光谱图及其逼近光谱图。逼近光谱图分析结果分别为：甲烷：1.08％、乙烷：0.03％、丙烷：0.029％、异丁烷：0.0085％、正丁烷：0.02％、异戊烷：0.01％和正戊烷0.008％；

具体实施方式

以下以甲烷、乙烷、丙烷、异丁烷、正丁烷、异戊烷和正戊烷的中红外吸收光谱定量分析为例，结合实际仪扫描的光谱与光谱重构结果等附图，以及本发明的原理对本发明的实施方式作进一步的详细说明。

本发明通过构建一个勾型损失函数：

P＝k₁ε²+(exp(k₂ε)-1)²     (1)

式中，exp(·)表示自然对数函数；ε是自变量；P表示损失函数；参数k₁和k₂控制函数的曲线形状。明显地，该函数中，权重在光谱透射率小于重构光谱相应波数透射率时较小，反之则大，从而达到了减小被分析混合气中存在的未知气体对目标气体分析结果的影响。

具体包括下列步骤：

1)预估N组分气体的目标气A₁、A₂、…、A_N的浓度分别为：C_A1、C_A2、…、C_AN，估计的方法可以采用预设浓度值，也可以采用特征谱线的最小二乘法估计、加权最小二乘估计法、反最小二乘估计法等方法求得；

2)估计各目标气体在所感兴趣的光谱段各组分气体在各自浓度情况下各数据点的折算吸光率δ；

3)在光谱图中感兴趣的光谱段内，根据Lambert-Beer定理，重构光谱；

4)求取重构偏差与损失函数值；

5)用神经网络学习方法中的最速下降法、动量法、遗传算法、模拟退火算法等对目标气体的浓度进行迭代运算。如果气体浓度增量超出最大增量范围，则令其为增量范围边界值；

6)用迭代结果取代1)中的预估值，重复上述步骤2)、3)、4)、5)，直到每种目标气的浓度与上一次迭代结果之差小于设定的容许偏差。

下面以甲烷、乙烷、丙烷、异丁烷、正丁烷、异戊烷和正戊烷的计算为例，说明本发明专利的具体实施方式。

(1)预估目标气体中各组分的浓度

对于图1(a)所示的混合气体的光谱图，预设7组分气体的目标气的浓度分别为：5％、1.5％、1％、0.5％、0.5％、0.2％和0.2％；对于图1(b)所示的混合气体的光谱图，预设7组分气体的目标气的浓度分别为：1％、0、0、0、0、0和0。

(2)估计目标气体各组分在所感兴趣的光谱段各数据点的折算吸光率

浓度均为1％的单组份甲烷、乙烷、丙烷、异丁烷、正丁烷、异戊烷和正戊烷气体在中红外400到1400波数段的吸收光谱如附图2所示，该光谱图是采用Bruker公司Alpha光谱仪获得，分辨率为4波数，输出值为透光率，因此各气体组分的折算吸光率为

δ_ij＝log(v_ij)/C_Aj＝log(v_ij)(i＝1，2，…，522；A_j∈{甲烷、乙烷、丙烷、异丁烷、正丁烷、

                                                   异戊烷和正戊烷})               (2)

式中log(·)表示自然对数运算；v_i表示目标气第j组分A_j谱图中所感兴趣波数段中第i谱线的光谱值；δ_ij表示目标气第j组分A_j谱图中所感兴趣波数段中第i谱线的的折算吸光率，它是实际吸光率与光程的乘积；C_Aj表示目标气第j组分A_j的浓度。考虑到波数1300～1500的中红外光谱段受水汽影响较大，波数400～750的中红外光谱段包含上述七种气体浓度信息少，且受二氧化碳影响较大，这里所感兴趣波数段只考虑波数750～1300的中红外光谱段。

对于图1中所示的光谱图，除甲烷外，其它气体浓度较小，各自的吸光率在较小浓度变化范围内不会有较大差异，可以选取为1％时的吸光率。而图1(a)甲烷的浓度较高，其折算吸光率与1％时相比相对要小一些，因此需要修正。修正的方法可以是分段线性化，也可以是样条插值法、多项式回归法。为方便起见，本实例采用分段线性化。如果甲烷单组分样本中，其浓度分别为：1％，3％，10％，30％，60％和100％，甲烷预设浓度为5％，那么，采用分段线性化可计算甲烷第i条谱线的吸光率：

${\delta }_{i1}=\mathrm{rate}\times {\delta }_{i1}^{10}+(1-\mathrm{rate})\times {\delta }_{i1}^3---\left(3\right)$

式中rate＝(5-3)/(10-3)＝2/7；和表示甲烷浓度分别为3％和10％时，甲烷第i条谱线的吸光率。

(3)在光谱图750到1300波数段，根据Lambert-Beer定理，获得初始重构吸收光谱图：

$v_i=\exp (-\underset{j=1}{\overset{N}{\Sigma }}{\delta }_{\mathrm{ij}}{C}_{\mathrm{Aj}})---\left(4\right)$

式中exp(·)表示自然指数运算，N表示目标气体的种数，这里N＝7；

(4)计算重构偏差与损失函数值

ε_i＝u_i-v_i    (5)

选择勾型函数作为损失函数：

$f(C_{A1},C_{A2},\Lambda ,C_{\mathrm{AN}})=\underset{i}{\overset{M}{\Sigma }}(k_1{ϵ}_i^2+{(\exp \left(k_2{ϵ}_i\right)-1)}^2)---\left(6\right)$

式中u_i表示实际扫描光谱图中所感兴趣波数段的第i条谱线的值；M表示所重构的吸收光谱段的光谱值个数；k₁、k₂为式(1)中的勾型损失函数的形状系数。这里选择k₁＝1、k₂＝100，对应于式(1)表示的勾型函数的曲线如图3所示。

(5)用神经网络学习方法中的动量法对目标气体的浓度进行迭代运算

C_Aj(t)＝C_Aj(t-1)-α(j)[(1-λ)ΔC_Aj(t-1)+γΔC_Aj(t-2)]    (7)

式中α(j)为待分析气体第j组分A_j浓度定量分析的学习率，λ为动量因子，且0≤λ≤1，t为学习迭代次。对于ΔC_Aj，可采用式(8)来估计：

${\mathrm{\Delta C}}_{\mathrm{Aj}}=\frac{\partial f\left(E\right)}{{\partial C}_{\mathrm{Aj}}}\approx \frac{f(C_1,C_2,\Lambda ,C_{\mathrm{Aj}}+e,\Lambda ,C_{\mathrm{AN}})}{e}---\left(8\right)$

式中e是一个足够小量，可取e＝0.01。如果ΔC_Aj＞Max_Aj，则令ΔC_Aj＝Max_Aj，如果ΔC_Aj＜Min_Aj，则令ΔC_Aj＝Min_Aj。其中Max_Aj和Min_Aj分别表示设定的气体组分A_j浓度的增量上限和下限，且Max_Aj＞0，Min_Aj＜0。本实例中设定七组分烷烃气体对应的Max_Aj分别为：0.2，0.1，0.1，0.05，0.05，0.02和0.02，Min_Aj分别为：-0.2，-0.1，-0.1，-0.05，-0.05，-0.02和-0.02；

(6)用迭代结果取代(2)中的预估值，重复上述步骤(3)、(4)、(5)，直到每种目标气的浓度与上一次迭代结果之差小于设定的容许偏差E。对于图1(a)，这里选择E＝[0.10.050.050.020.020.010.01]，对于图1(b)，这里选择E＝[0.050.010.010.010.010.010.01]。最后得到的逼近光谱图如图4所示。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施方式仅限于此，对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单的推演或替换，都应当视为属于本发明由所提交的权利要求书确定专利保护范围。