CN102680425A - 多组分气体傅立叶变换光谱分析的多分析模型信息融合方法 - Google Patents

Abstract

本发明公开了一种多组分气体傅立叶变换红外光谱在线分析的信息融合方法，该方法首先采用常规多组分气体的特征变量提取方法、分析模型建立方法，为每种目标气体建立多个分析模型，然后建立每个模型的噪声水平估计模型。对于给定的光谱图，计算每个模型的结果，估计其噪声水平，并根据各模型的噪声水平估计其权值系数。最后对同一种目标气体的分析模型进行加权平均，计算该组分气体的最终分析结果，以提高多组分气体光谱分析结果的准确度。本方法可用于石油、天然气勘探的气测录井，矿井安全，产品质量检查与故障诊断、五金、化工、环境保护等涉及多组分气体光谱定量分析的应用领域中。

Description

多组分气体傅立叶变换光谱分析的多分析模型信息融合方法

技术领域：

本发明涉及光谱在线定量分析领域，涉及一种多组分气体傅立叶变换红外光谱多分析模型的信息融合方法。 

背景技术：

傅立叶变换红外光谱分析是多组分气体在线定量分析的重要方法之一。气体分析涉及到科学研究、矿井安全、环境保护、设备的故障诊断、产品质量的检测、天然气与石油勘探录井等多个领域。在这些领域中，组分气体的在线定量分析发挥着非常重要的作用。如：变压器的故障诊断方法中，有一种就是通过检测变压器油中乙炔等气体的含量来识别的，天然气的燃烧值是通过检测天然气中甲烷、乙烷等组分的含量来确定的。虽然气体吸收光谱的光谱分析是气体定量分析的一种重要方法，但这种方法还存在一些没有解决的问题，其中之一就是光谱基线的漂移，待分析气体中可能存的干扰气体，谱线的噪声会给分析结果带来偏差。在存在多种具有相同分子基团的组分情况下，如丙烷、异丁烷、正丁烷等烷烃气体都具有-CH₂-基团，由于其吸收光谱严重交叠，光谱所形成的矩阵条件数大，光谱中一点小的干扰就可能引入大的偏差。为了尽可能得到更为准确的分析结果，需要想办法降低这些偏差的水平。此外，光谱线的灵敏度与气体的大浓度分析范围总是一对矛盾的问题。如何在大浓度范围与灵敏度之间进行无缝连接，在小浓度范围时能实用高吸光度谱线进行气体分析，而在大浓度范围时，考虑到高吸光度谱线已经吸收饱和，能用合适吸光度的谱线进行分析。 

发明内容：

本发明的目的在于，提供一种多组分气体傅立叶变换红外光谱分析的多模型信息融合方法，提高分析结果准确性。 

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

一种多组分气体傅立叶变换红外光谱分析的多模型信息融合方法，该方法为每种气体的提取多个特征变量，然后根据特征变量之间的不同组合建立多个分析模型，并估计各模型在不同气体浓度时的噪声水平，在实际的多组分气体光谱在线分析过程中，根据各模型的分析结果及噪声水平，对其分析结果进行加权平均，以获得更为准确的分析结果。 

特征变量提取的方法是Tikhonov正则化方法或者是前向选择法； 

对于以吸光度为输出的光谱图，Tikhonov正则化法是基于式(1)的线性模型的： 

Y_i＝XB_i+E    (1a) 

式中X包含m个标定样本在w条谱线上的光谱；Y_i是m×1的向量，它包含了第i种被分析气体的浓度信息；B_i是w×1的回归向量，表示第i种气体的特征变量系数；E也是m×1的向量，表示随机误差；对于以透射率为输出的光谱图，模型修改为： 

Y_i＝-ln(X)B_i+E    (1b) 

式中ln(·)是表示自然对数运算。是将(1)式转化为式(2)最优化问题求解： 

$\min ({\left|\right|\mathrm{XB}-y\left|\right|}_a^a+\mathrm{\λ}{\left|\right|\mathrm{LB}\left|\right|}_b^b)---\left(2\right)$

式中‖·‖_p表示回归向量的p-模，p＝2时就是欧几里模；1≤a，b＜∞，L表示一个正则化算子，它迫使B的估计值属于相应的特性良好的函数的子空间；λ表示控制第2项相对于第1项的权值；min(·)表示最小值求取算子；明显地，式(3)中的第一项是最小化式(2)中的随机误差E，而第2项则是为了最小化回归向量B，以获得灵敏度高、信噪比高的特征谱线，以提高标定模型分析结果的准确性； 

对于前向选择法，其方法是观察并比对相同浓度的各组分气体的吸收光谱，估计各组分气体吸收光谱的交叠程度，并按照交叠程度将这些光谱图进行排序，首先为交叠度最低的气体选择特征变量，然后依次为每种气体提取特征变量；其步骤如下： 

(1)为第一种气体寻找两条或多条谱线的组合，使得其对该种气体的灵敏度最高，而对其它气体的灵敏度尽可能得低； 

(2)为第二种气体的寻找两条或多条谱线的组合，使得它对该种气体的灵敏度大于第三种及以后的气体的灵敏度； 

(3)为第三种气体的寻找两条或多条谱线的组合，使得它对该种气体的灵敏度大于第四种及以后的气体的灵敏度。 

第四种及以后的气体，以此类推。 

各组分气体的分析模型采用多项式偏最小二乘法，或者采用前向神经网络法来建立；每个模型的输入集合中，包含每种气体的一个特征变量；特征变量组合中存在某个空间，所有的特征变量所包含的谱线均没有吸收饱和，即透射率在0.1以上； 

采用偏最小二乘法建模的模型时，其形式为： 

$y_{\mathrm{im}}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{a}_j{f}_{\mathrm{im}}+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{p}_{\mathrm{kj}}{f}_{\mathrm{km}}{f}_{\mathrm{im}}+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{q}_{\mathrm{kj}}{f}_{\mathrm{km}}{f}_{\mathrm{jm}}^2+e---\left(3\right)$

式中X_m表示第m个样本的光谱值向量；y_im表示第i种气体的第m个估计浓度；f_im＝X_mB_i表示第i种气体的第m个样本的特征变量值；e表示模型的残差；N表示气体的种类数；a_j、p_kj和q_kj是模型的待定系数，其估计公式为： 

P_i＝inv(F_iΛF_i ^T)F_i ^TΛY_i(4) 

式中P_i＝[a_j，p_kj，q_kj]，表示第i种气体多项式偏最小二乘法模型的系数向量，F_i＝[f₁，f₂，…，f_N]，Λ为用户设定的权值，Y_i＝[y_i1，y_i2，…，y_iM]^T，上标T表示转置算子，M表示样本的数量。 

采用前向神经网络法进行分析模型建模时，网络的隐层节点响应函数为S型函数 

$s_i=\frac{2}{1+e^{-t_i}}-1,$ $t_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{ji}}{f}_j,(i=\mathrm{1,2},...,S)---\left(5\right)$

式中s_i表示隐层第i个节点的输出，t_i表示隐层第i个节点的总输入，W_ji表示隐层第i个节点与输入层第j个节点的连接权值，隐层一共有S个节点，S＞N；输出层节点响应函数为线性函数： 

y_i＝r_i+d_i， $r_i=\underset{j=1}{\overset{s}{\mathrm{\Σ}}}{V}_{\mathrm{ji}}{s}_j,(i=\mathrm{1,2},...,N)---\left(6\right)$

式中r_i表示输出层第i个节点的总输入，V_ji表示隐层第j个节点与输出层第i个节点的连接权值，y_i表示输出层第i个节点的输出，d_i表示输出层第i个节点截长，网络的参数W_ji、V_ji和d_i学习算法采用动量法： 

W_ji(k+1)＝W_ji(k)-α(i)[(1-λ)ΔW_ji(k)+λΔW_ji(k-1)](7a) 

V_ji(k+1)＝V_ji(k)-β(i)[(1-λ)ΔV_ji(k)+λΔV_ji(k-1)](7b) 

d_i(k+1)＝d_i(k)-γ(i)[(1-λ)Δd_i(k)+λΔd_i(k-1)](7c) 

式中α(i)、β(i)和γ(i)分别为气体权值W_ji、V_ji和截长d_i的学习率，λ为动量因子，且0≤λ≤1；对于ΔW_ji、ΔV_ji和Δd_i分别用下式： 

$\mathrm{\Δ}{W}_{\mathrm{ji}}=\frac{g(W_{j1},...,W_{\mathrm{ji}}+e,...,W_{\mathrm{jS}},V_{\mathrm{ji}},d_i)-g(W_{j1},...,W_{\mathrm{ji}},...,W_{\mathrm{jS}},V_{\mathrm{ji}},d_i)}{e}---\left(8a\right)$

$\mathrm{\Δ}{V}_{\mathrm{ji}}=\frac{g(V_{j1},...,V_{\mathrm{ji}}+e,...,V_{\mathrm{jN}},W_{\mathrm{ji}},d_i)-g(V_{j1},...,V_{\mathrm{ji}},...,V_{\mathrm{jN}},W_{\mathrm{ji}},d_i)}{e}---\left(8b\right)$

$\mathrm{\Δ}{d}_{\mathrm{ji}}=\frac{g(W_{\mathrm{ji}},V_{\mathrm{ji}},d_1,...,d_i+e,...,d_N)-g(W_{\mathrm{ji}},V_{\mathrm{ji}},d_1,...,d_i,...,d_N)}{e}---\left(8c\right)$

式中e是一个足够小的量，可取e=0.001，函数f(·)是由附图2的前向神经网络结构在式(5)和(6)响应函数情况下形成的函数。 

各组分气体的每个模型的信噪比采用如下步骤估计： 

1)估计每条谱线的噪声水平；在气室中充满氮气，并连续扫描光谱8～15次左右，求每条谱线的均方根。对于以透射率为输出的光谱，按式(9)计算： 

$\mathrm{\Δ}{x}_i=\sqrt{\frac{\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{(x_{\mathrm{ij}}-{\stackrel{\‾}{x}}_i)}^2}{L-1}}---\left(9\right)$

式中L表示扫描的次数，Δx_i表示第i条谱线的均方差，x_ij表示第j次扫描第i条谱线的谱线值， 

表示L次光谱扫描结果中第i条谱线的均值；而对于以吸光度为输出的光谱，每条谱线值进行指数运算exp(-v_ij)后代入式(9)进行计算； 

2)估计特征变量f_i的噪声水平，对于以透射率和吸光度为输出的光谱图，给定光谱数据X以后，特征变量的噪声分别由式(10a)和(10b)估计，或者分别由式(11a)和(11b)估计： 

$\mathrm{\Δ}{f}_i=2\underset{j=1}{\overset{R}{\mathrm{\Σ}}}\left|\right|b_j\left|\right|\mathrm{\Δ}{x}_j/x_j---\left(10a\right)$

$\mathrm{\Δ}{f}_i=2\underset{j=1}{\overset{R}{\mathrm{\Σ}}}\left|\right|b_j\left|\right|\mathrm{\Δ}{x}_j/\exp (-x_j)---\left(10b\right)$

$\mathrm{\Δ}{f}_i=2\sqrt{\underset{j=1}{\overset{R}{\mathrm{\Σ}}}{(b_j\mathrm{\Δ}{x}_j/x_j)}^2}---\left(11a\right)$

$\mathrm{\Δ}{f}_i=2\sqrt{\underset{j=1}{\overset{R}{\mathrm{\Σ}}}{(b_j\mathrm{\Δ}{x}_j/\exp (-x_j))}^2}---\left(11b\right)$

式中b_j是由特征变量提取算法确定的向量B_i的第j个元素，这里R表述B_i的长度； 

3)估计每个分析模型的噪声水平。首先估计每组分气体分析模型对每个输入f_i的微分： 

$\mathrm{\Δ}{g}_{\mathrm{jki}}=\frac{g_{\mathrm{jk}}(f_1,...,f_i+e,...,f_N)-g_{\mathrm{jk}}(f_1,...,f_i,...,f_N)}{e}---\left(12\right)$

式中g_jk(·)表示由步骤2)的前向神经网络或多项式确定的第j组分气体的第k个分析模型，e是一个很小的量，直接设定为0.0005～0.005之间的数；然后按式(13)进行该模型的噪声的估计： 

${\mathrm{\δ}}_{\mathrm{jk}}=\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left(\mathrm{\Δ}{g}_{\mathrm{jki}}\mathrm{\Δ}{f}_i\right)}^2}---\left(13\right)$

式中δ_jk表示第j组分气体的第k个分析模型的噪声水平。 

各组分气体的最终分析结果的信息融合按式(14)进行： 

${\tilde{y}}_j=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{y}_{\mathrm{jk}}\frac{{\mathrm{\τ}}_k}{{\mathrm{\κ}}_j}---\left(14\right)$

式中y_jk表示第j组分气体的第k个分析模型的分析结果， 

表示第j组分气体的融合结果，K表示第j组分气体的分析模型个数，κ_j和τ_k是由分模型噪声水平估计的系数，由式(15)估计： 

${\mathrm{\κ}}_j=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\τ}}_k,$ ${\mathrm{\τ}}_k=\frac{\underset{l=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\δ}}_{\mathrm{jl}}}{{\mathrm{\δ}}_{\mathrm{jk}}}---\left(15\right).$

该方法为每种气体的提取多个特征变量，然后根据特征变量之间的不同组合建立多个分析模型，并估计各模型在不同气体浓度时的噪声水平，在实际的多组分气体光谱在线分析过程中，根据各模型的分析结果及噪声水平，对其分析结果进行加权平均，以获得更为准确的分析结果。本方法可用于石油、天然气勘探的气测录井，矿井安全，产品质量检查与故障诊断、五金、化工、环境保护等涉及多组分气体光谱定量分析的应用领域中。 

附图说明：

图1(a)是浓度均为1％的甲烷、乙烷、丙烷、异丁烷和正丁烷五组分气体的傅立叶变换中红外吸收光谱图； 

图1(b)是浓度均为1000ppm的甲烷、乙烷、丙烷、异丁烷和正丁烷五组分气体在2800～3100波数段的吸收光谱图； 

图1(c)是浓度均为1％的甲烷、乙烷、丙烷、异丁烷和正丁烷五组分气体的低波数段傅立叶变换中红外吸收光谱图； 

图2(a)多输出两层前向神经网络结构图； 

图2(b)单输出两层前向神经网络结构图。 

具体实施方式：

以下以甲烷、乙烷、丙烷、异丁烷、正丁烷、异戊烷和正戊烷的傅立叶变换中红外吸收光谱在线分析为例，来说明本发明的具体实施方式： 

不失一般性，本发明以石油、天然气勘探气测录井过程中所需要的分析的气体的光谱 分析为例，结合本发明的原理对本发明的实施方式作进一步的详细说明。气测录井所要分析的气体包括甲烷、乙烷、丙烷、异丁烷和正丁烷，其傅立叶变换中红外吸收光谱图如图1所示。 

1)在不同的波数段范围内为每种气体提取不同特征变量。提取的方法可以是Tikhonov正则化方法，也可以使是前向选择法。 

对于以吸光度为输出的光谱图，Tikhonov正则化法是基于式(1)的线性模型的： 

Y_i＝XB_i+E    (1a) 

式中X包含m个标定样本在w条谱线上的光谱；Y_i是m×1的向量，它包含了第i种被分析气体的浓度信息；B_i是w×1的回归向量，表示第i种气体的特征变量系数；E也是m×1的向量，表示随机误差。对于以透射率为输出的光谱图，模型修改为： 

Y_i＝-ln(X)B_i+E    (1b) 

式中ln(·)是表示自然对数运算。是将(1)式转化为式(2)最优化问题求解： 

$\min ({\left|\right|\mathrm{XB}-y\left|\right|}_a^a+\mathrm{\λ}{\left|\right|\mathrm{LB}\left|\right|}_b^b)---\left(2\right)$

式中‖·‖_p表示回归向量的p-模，p＝2时就是欧几里模；1≤a，b＜∞，L表示一个正则化算子，它迫使B的估计值属于相应的特性良好的函数的子空间；λ表示控制第2项相对于第1项的权值；min(·)表示最小值求取算子。明显地，式(3)中的第一项是最小化式(2)中的随机误差E，而第2项则是为了最小化回归向量B，以获得灵敏度高、信噪比高的特征谱线，以提高标定模型分析结果的准确性； 

对于前向选择法，其方法是观察并比对相同浓度的各组分气体的吸收光谱，估计各组分气体吸收光谱的交叠程度，并按照交叠程度将这些光谱图进行排序，首先为交叠度最低的气体选择特征变量，然后依次为每种气体提取特征变量。其步骤如下： 

(1)为第一种气体寻找两条或多条谱线的组合，使得其对该种气体的灵敏度最高，而对其它气体的灵敏度尽可能得低； 

(2)为第二种气体的寻找两条或多条谱线的组合，使得它对该种气体的灵敏度大于第三种及以后的气体的灵敏度； 

(3)为第三种气体的寻找两条或多条谱线的组合，使得它对该种气体的灵敏度大于第四种及以后的气体的灵敏度。 

其它依次类推。 

2)采用多项式偏最小二乘法、前向神经网络法为每种目标气体在不同的浓度范围段内采用在不同波数段提取的特征变量组合进行分析模型建模，每个模型的输入集合中，包含 每种气体的一个特征变量。特征变量组合中存在某个空间，所有的特征变量所包含的谱线均没有吸收饱和，即透射率在0.1以上； 

采用偏最小二乘法建模的模型时，其形式为： 

$y_{\mathrm{im}}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{a}_j{f}_m+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{p}_{\mathrm{kj}}{f}_{\mathrm{km}}{f}_{\mathrm{jm}}+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{q}_{\mathrm{kj}}{f}_{\mathrm{km}}{f}_{\mathrm{jm}}^2+e---\left(3\right)$

式中X_m表示第m个样本的光谱值向量；y_im表示第i种气体的第m个估计浓度；f_im＝X_mB_i表示第i种气体的第m个样本的特征变量值；e表示模型的残差；N表示气体的种类数；a_j、p_kj和q_kj是模型的待定系数，其估计公式为： 

P_i＝inv(F_iΛF_i ^T)F_i ^TΛY_i(4) 

式中P_i＝[a_j，p_kj，q_kj]，表示第i种气体多项式偏最小二乘法模型的系数向量，F_i＝[f₁，f₂，…，f_N]，Λ为用户设定的权值，Y_i＝[y_i1，y_i2，…，y_iM]^T，上标T表示转置算子，M表示样本的数量。 

采用前向神经网络法进行分析模型建模时，其模型如附图2所示，网络隐层节点响应函数为S型函数 

$s_i=\frac{2}{1+e^{-t_i}}-1,$ $t_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{ji}}{f}_j,(i=\mathrm{1,2},...,S)---\left(5\right)$

式中s_i表示隐层第i个节点的输出，t_i表示隐层第i个节点的总输入，W_ji表示隐层第i个节点与输入层第j个节点的连接权值，隐层一共有S个节点，S＞N。输出层节点响应函数为线性函数： 

y_i＝r_i+d_i， $r_i=\underset{j=1}{\overset{S}{\mathrm{\Σ}}}{V}_{\mathrm{ji}}{s}_j,(i=\mathrm{1,2},...,N)---\left(6\right)$

式中r_i表示输出层第i个节点的总输入，V_ji表示隐层第j个节点与输出层第i个节点的连接权值，y_i表示输出层第i个节点的输出，d_i表示输出层第i个节点截长，网络的参数W_ji、V_ji和d_i学习算法采用动量法： 

W_ji(k+1)＝W_ji(k)-α(i)[(1-λ)ΔW_ji(k)+λΔW_ji(k-1)](7a) 

V_ji(k+1)＝V_ji(k)-β(i)[(1-λ)ΔV_ji(k)+λΔV_ji(k-1)](7b) 

d_i(k+1)＝d_i(k)-γ(i)[(1-λ)Δd_i(k)+λΔd_i(k-1)](7c) 

式中α(i)、β(i)和γ(i)分别为气体权值W_ji、V_ji和截长d_i的学习率，λ为动量因子，且0≤λ≤1。对于ΔW_ji、ΔV_ji和Δd_i分别用下式： 

$\mathrm{\Δ}{W}_{\mathrm{ji}}=\frac{g(W_{j1},..,W_{\mathrm{ji}}+e,...,W_{\mathrm{jS}},V_{\mathrm{ji}},d_i)-g(W_{j1},...,W_{\mathrm{ji}},...,W_{\mathrm{jS}},V_{\mathrm{ji}},d_i)}{e}---\left(8a\right)$

$\mathrm{\Δ}{V}_{\mathrm{ji}}=\frac{g(V_{j1},...,V_{\mathrm{ji}}+e,...,V_{\mathrm{jN}},W_{\mathrm{ji}},d_i)-g(V_{j1},...,V_{\mathrm{ji}},...,V_{\mathrm{jN}},W_{\mathrm{ji}},d_i)}{e}---\left(8b\right)$

$\mathrm{\Δ}{d}_{\mathrm{ji}}=\frac{g(W_{\mathrm{ji}},V_{\mathrm{ji}},d_1,...,d_i+e,...,d_N)-g(W_{\mathrm{ji}},V_{\mathrm{ji}},d_1,...,d_i,...,d_N)}{e}---\left(8c\right)$

式中e是一个足够小的量，可取e=0.001，函数f(·)是由附图2的前向神经网络结构在式(5)和(6)响应函数情况下形成的函数； 

3)估算每组分气体的每个模型的信噪比 

(1)估计每条谱线的噪声水平。在气室中充满氮气，并连续扫描光谱8～15次左右，求每条谱线的均方根。对于以透射率为输出的光谱，按式(9)计算： 

$\mathrm{\Δ}{x}_i=\sqrt{\frac{\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{(x_{\mathrm{ij}}-{\stackrel{\‾}{x}}_i)}^2}{L-1}}---\left(9\right)$

式中L表示扫描的次数，Δx_i表示第i条谱线的均方差，x_ij表示第j次扫描第i条谱线的谱线值， 

表示L次光谱扫描结果中第i条谱线的均值。而对于以吸光度为输出的光谱，每条谱线值进行指数运算exp(-v_ij)后代入式(9)进行计算； 

(2)估计特征变量f_i的噪声水平。对于以透射率和吸光度为输出的光谱图，给定光谱数据X以后，特征变量的噪声分别由式(10a)和(10b)估计，或者分别由式(11a)和(11b)估计： 

$\mathrm{\Δ}{f}_i=2\underset{j=1}{\overset{R}{\mathrm{\Σ}}}\left|\right|b_j\left|\right|\mathrm{\Δ}{x}_j/x_j---\left(10a\right)$

$\mathrm{\Δ}{f}_i=2\underset{j=1}{\overset{R}{\mathrm{\Σ}}}\left|\right|b_j\left|\right|\mathrm{\Δ}{x}_j/\exp (-x_j)---\left(10b\right)$

$\mathrm{\Δ}{f}_i=2\sqrt{\underset{j=1}{\overset{R}{\mathrm{\Σ}}}{(b_j\mathrm{\Δ}{x}_j/x_j)}^2}---\left(11a\right)$

$\mathrm{\Δ}{f}_i=2\sqrt{\underset{j=1}{\overset{R}{\mathrm{\Σ}}}{(b_j\mathrm{\Δ}{x}_j/\exp (-x_j))}^2}---\left(11b\right)$

式中b_j是由特征变量提取算法确定的向量B_i的第j个元素，这里R表述B_i的长度。 

(3)估计每个分析模型的噪声水平。首先估计每组分气体分析模型对每个输入f_i的微分： 

$\mathrm{\Δ}{g}_{\mathrm{jki}}=\frac{g_{\mathrm{jk}}(f_1,...,f_i+e,...,f_N)-g_{\mathrm{jk}}(f_1,...,f_i,...,f_N)}{e}---\left(12\right)$

式中g_jk(·)表示由步骤2)的前向神经网络或多项式确定的第j组分气体的第k个分析模型， e是一个很小的量，可以直接设定为0.0005～0.005之间的数。然后按式(13)进行该模型的噪声的估计： 

${\mathrm{\δ}}_{\mathrm{jk}}=\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left(\mathrm{\Δ}{g}_{\mathrm{jki}}\mathrm{\Δ}{f}_i\right)}^2}---\left(13\right)$

式中δ_jk表示第j组分气体的第k个分析模型的噪声水平； 

4)根据各分析模型的噪声水平进行最终分析结果的信息融合，可按式(14)进行： 

${\tilde{y}}_j=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{y}_{\mathrm{jk}}\frac{{\mathrm{\τ}}_k}{{\mathrm{\κ}}_j}---\left(14\right)$

式中y_jk表示第j组分气体的第k个分析模型的分析结果， 

表示第j组分气体的融合结果，K表示第j组分气体的分析模型个数，κ_j和τ_k是由分模型噪声水平估计的系数，由式(15)估计： 

${\mathrm{\κ}}_j=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\τ}}_k,$ ${\mathrm{\τ}}_k=\frac{\underset{l=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\δ}}_{\mathrm{jl}}}{{\mathrm{\δ}}_{\mathrm{jk}}}---\left(15\right)$

实施例 

不失一般性，本实例以甲烷、乙烷、丙烷、异丁烷和正戊烷的分析为例来说明本发明的实施方式。这五组分气体的红外吸收光谱图如附图1所示，其中甲烷、乙烷为目标气体，甲烷的分析溶度范围为0～100％，乙烷的分析溶度范围为0～50％，其它为干扰气体，浓度范围为0～1％。从图1中可以看出，低波数段光谱交叠小，但分辨率也低，高波数段光谱分辨率高，但交叠严重。本实施例可按如下步骤完成： 

1)分别在波数段2850～3100cm^-1范围内，以及650～1400cm^-1范围内为每组分气体提取特征变量。 

观察图1(b)可直接看出，3016.5cm^-1处，甲烷有较大的钉状吸收峰，其它气体在该处吸收较小，且随着波数的减小，各组分气体的吸光度有所增强，于是，取该钉状吸收峰谱线值及其两端的谱线值的差值，甲烷的特征值明显地比其它组分气体要高很多，因此是甲烷一个比较好的特征变量。由于图1是以透光率为输出的谱图，各谱线值取对数以后再进行上述组合，该特征变量的线性度将更好一些，于是该特征变量可化为式(16) 

f_m1＝ln(v_3025.1)+ln(v_3010.8)-2ln(v_3016.5)(16) 

式中v_c表示波数为c的谱线值，f_m1表示甲烷的第一个特征变量； 

观察图1(b)还可直接看出，3023～3080cm^-1范围内，随着波数的减小，乙烷的吸光度增大，甲烷振荡性略有减小，其它气体只有轻微的增大，因此，这段光谱范围内取两条谱线的差值，将是乙烷的一个备选特征变量 

f_e1＝ln(v_3079.7)-ln(v_3023.7)(17) 

式中f_e1表示乙烷的第一个特征变量； 

观察图1(b)还可直接看出，甲烷的钉状吸收峰两端，随着波数的减小，丙烷的吸光度增大最快，甲烷基本不变，因此，这两条谱线的差值，将是丙烷的一个良好备选特征变量 

f_p1＝ln(v_3025.1)-ln(v_3010.8)(18) 

式中f_p1表示丙烷的第一个特征变量； 

观察图1(b)还可直接看出，甲烷的钉状吸收峰两端，随着波数的减小，丙烷的吸光度增大最快，甲烷基本不变，因此，这两条谱线的差值，将是丙烷的一个良好备选特征变量 

f_p1＝ln(v_3025.1)-ln(v_3010.8)(19) 

式中f_p1表示丙烷的第一个特征变量； 

观察图1(b)还可直接看出，2943～2953cm^-1范围内，异丁烷谱图曲线的斜率最大，也就是说变化最快，因此在这段波数范围内，取两条谱线的差值，将是异丁烷的一个良好备选特征变量： 

f_p1＝ln(v_2943.2)-ln(v_2953.3)(20) 

式中f_p1表示异丁烷的第一个特征变量； 

观察图1(b)还可直接看出，2943～2953cm^-1范围内，异丁烷谱图曲线的斜率最大，也就是说变化最快，因此在这段波数范围内，取两条谱线的差值，将是异丁烷的一个良好备选特征变量： 

f_i1＝ln(v_2943.2)-ln(v_2953.3)(21) 

式中f_i1表示异丁烷的第一个特征变量； 

观察图1(b)还可直接看出，2898～2919cm^-1范围内，正丁烷谱图曲线的斜率最大，也就是说变化最快，因此在这段波数范围内，取两条谱线的差值，将是正丁烷的一个良好备选特征变量： 

f_n1＝ln(v_2898.7)-ln(v2_918.8)(22) 

式中f_n1表示正丁烷的第一个特征变量； 

观察图1(c)可直接看出，1210～1250cm^-1范围内，随着波数的增大，甲烷的吸光度增大，其它气体只有轻微的增大，因此，这段光谱范围内取两条谱线的差值，将是甲烷的一 个备选特征变量。考虑到甲烷的分析浓度范围大，吸光度高的谱线在甲烷浓度高时容易吸收饱和，这里分别选择 

f_m2＝ln(v_1207.2)-ln(v_1231.6)(23) 

f_m3＝ln(v_1207.2)-ln(v_1221.5)(24) 

式中f_m2和f_m3分别表示甲烷的第二和第三个特征变量； 

观察图1(c)还可直接看出，825～1105cm^-1范围内，随着波数的减小，乙烷的吸光度增大，其它气体吸光度有轻微的减小，因此，这段光谱范围内取两条谱线的差值，将是乙烷的一个备选特征变量 

f_e2＝ln(v_1105.1)-ln(v_827.7)(25) 

式中f_e2表示乙烷的第二个特征变量； 

观察图1(c)还可直接看出，748cm^-1附近丙烷有钉状吸收峰，因此，做类是f_m1的运算，将是丙烷的一个良好备选特征变量 

f_p2＝ln(v_744.4)+ln(v_754.5)-2ln(v_748.8)(26) 

式中f_p2表示丙烷的第二个特征变量； 

观察图1(c)还可直接看出，1177cm^-1附近，异丁烷谱图有钉状吸收峰，也就是说变化最快，因此在这段波数范围内，取两条谱线的差值，将是异丁烷的一个良好备选特征变量： 

f_i2＝ln(v_2943.2)-ln(v_2953.3)(28) 

式中f_i2表示异丁烷的第二个特征变量； 

观察图1(c)还可直接看出，970cm^-1附近，正丁烷有一个较大的吸收峰，在该区域附近的低波数端，正丁烷吸收峰还没完全衰减完，丙烷和异丁烷的吸收峰已明显出现，因此，选择该吸收峰高波数端的一条谱线与峰值谱线做差值运算，将是正丁烷的一个良好备选特征变量： 

f_n2＝ln(v_1029.0)-ln(v_970.0)(29) 

式中f_n2表示正丁烷的第二个特征变量； 

2)为每种目标气体建立多个分析模型 

(1)所有目标气体的第一个特征变量是从高波数段2850～3100cm^-1提取的，目标气体在该波数段的吸光率比较高，因此，这些特征变量组合的一个集合作为输入建立的分析模型，可以作为低浓度情况下各组分气体的分析模型。选择多项式作为目标气体的分析模型： 

$c_{m1}=a_1{f}_{m1}+a_2{f}_{e1}+a_3{f}_{p1}+a_4{f}_{i1}+a_5{f}_{n1}+b_1{f}_{m1}^2+b_2{f}_{e1}^2$

$=1.56f_{m1}-0.018f_{e1}+0.00062f_{i1}+0.0043f_{n1}+1.21{f_{m1}}^2+0.0028{f_{e1}}^2---\left(30a\right)$

$c_{e1}=p_1{f}_{m1}+p_2{f}_{e1}+p_3{f}_{p1}+p_4{f}_{i1}+p_5{f}_{n1}+q_1{f}_{m1}^2+q_2{f}_{e1}^2$

$=0.06f_{m1}+604f_{e1}+0.0022f_{p1}+0.0025f_{n1}+0.012{f_{m1}}^2+0.12{f_{e1}}^2---\left(30b\right)$

式中c_m1表示甲烷第一个分析模型的浓度值；c_e1表示乙烷第一个分析模型的浓度值；a₁～a₅，以及b₁和b₂为甲烷第一个分析模型的系数，求得每个样本光谱的特征变量值后，由式(4)确定；p₁～p₅，以及q₁和q₂为乙烷第一个分析模型的系数，求得每个样本光谱的特征变量值后，由式(4)确定； 

(2)按照式(30a)的形式选择f_m2和其它组分气体的第2个特征变量组成输入集合，建立甲烷的第2个分析模型，模型输出用c_m2表示；按照式(30b)的形式选择f_e2和其它组分气体的第2个特征变量组成输入集合，建立乙烷的第2个分析模型，模型输出用c_e2表示； 

(3)按照式(30a)的形式选择f_m3和其它组分气体的第2个特征变量组成输入集合，建立甲烷的第3个分析模型，模型输出用c_m3表示； 

3)估算每组分气体的每个模型的信噪比 

(1)估计每条谱线的噪声水平。在气室中充满氮气，并连续扫描光谱15次，把15个谱图中所要用到的谱线的值代入式(9)，求该谱线的噪声水平。即把式(16)～(29)中的v_x(x表示波长)对应的15个谱图值代入式(9)，求该谱线的噪声水平。例如，把式(17)中的v_3079.7和v_3023.6分别代入式(9)计算得到： 

$\mathrm{\Δ}{v}_{3079.7}=\sqrt{\frac{\underset{j=1}{\overset{15}{\mathrm{\Σ}}}{(x_{\mathrm{ij}}-{\stackrel{\‾}{x}}_i)}^2}{15-1}}=0.00046$

$\mathrm{\Δ}{v}_{3023.6}=\sqrt{\frac{\underset{j=1}{\overset{15}{\mathrm{\Σ}}}{(x_{\mathrm{ij}}-{\stackrel{\‾}{x}}_i)}^2}{15-1}}=0.00045$

(2)估计特征变量f_i的噪声水平。光谱是以透射率为输出的光谱图，于是，在所有气体浓度为0时，特征变量f_e1的噪声由式(10a)估计可得： 

$\mathrm{\Δ}{f}_{e1}=2\underset{j=1}{\overset{R}{\mathrm{\Σ}}}\left|\right|b_j\left|\right|\mathrm{\Δ}{x}_j/x_j=2(\mathrm{\Δ}{v}_{3097.7}/v_{3097.7}+\mathrm{\Δ}{v}_{3023.6}/v_{3023.6})=2(0.00046/1+0.00045/1)=0.00182$

采用同样的方法可以估计其它特征变量的噪声水平。 

(3)估计每个分析模型的噪声水平。首先估计每组分气体分析模型对每个输入f_i的微分，根据式(12)，对于如式(30b)所示乙烷的分析模型可得： 

Δg_e11＝0.06；Δg_e12＝6.4；Δg_e13＝0.0022；Δg_e14＝0；Δg_e15＝0.0025 

然后按式(13)进行该模型的噪声的估计可得乙烷的第一个分析模型的噪声水平为： 

${\mathrm{\δ}}_{e1}=\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left(\mathrm{\Δ}{g}_{\mathrm{jki}}\mathrm{\Δ}{f}_i\right)}^2}=0.01165$

按照同样的方法可求得其它分析模型的噪声水平。 

4)根据各分析模型的噪声水平进行最终分析结果的信息融合。首先根据各组分气体的分析模型的噪声水平按照式(15)估计系数τ_k和κ_j。乙烷的有两个分析模型，在各组分气体浓度为0时，其噪声水平分别为0.01165和0.23156，于是可得： 

τ₁＝(0.01165+0.23156)/0.01165＝20.8764；τ₂＝(0.01165+0.23156)/0.23156＝1.0507 

κ_j＝τ₁+τ₂＝21.9271 

然后根据求得的系数、各分析模型计算结果按照式(14)计算融合结果。在各组分气体浓度为0时，乙烷的两个分析模型的分析结果分别为0.001％8和0.0352％，于是其融合结果为： 

${\tilde{y}}_j=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{y}_{\mathrm{jk}}\frac{{\mathrm{\τ}}_k}{{\mathrm{\κ}}_j}=(20.8764\×0.0018\%+1.0507\×0.0352\%)/21.9271=0.0034\%.$

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的方法及技术内容作出些许的更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，仍属于本发明技术方案的范围内。 

Claims (5)

1.一种多组分气体傅立叶变换红外光谱分析的多模型信息融合方法，其特征在于，该方法为每种气体的提取多个特征变量，然后根据特征变量之间的不同组合建立多个分析模型，并估计各模型在不同气体浓度时的噪声水平，在实际的多组分气体光谱在线分析过程中，根据各模型的分析结果及噪声水平，对其分析结果进行加权平均，以获得更为准确的分析结果。

2.根据权利要求1所述的多模型信息融合方法，其特征在于，特征变量提取的方法是Tikhonov正则化方法或者是前向选择法；

对于以吸光度为输出的光谱图，Tikhonov正则化法是基于式(1)的线性模型的：

Y_i＝XB_i+E    (1a)

式中X包含m个标定样本在w条谱线上的光谱；Y_i是m×1的向量，它包含了第i种被分析气体的浓度信息；B_i是w×1的回归向量，表示第i种气体的特征变量系数；E也是m×1的向量，表示随机误差；对于以透射率为输出的光谱图，模型修改为：

Y_i＝-ln(X)B_i+E    (1b)

式中ln(·)是表示自然对数运算。是将(1)式转化为式(2)最优化问题求解：

$\min ({\left|\right|\mathrm{XB}-y\left|\right|}_a^a+\mathrm{\λ}{\left|\right|\mathrm{LB}\left|\right|}_b^b)---\left(2\right)$

式中‖·‖_p表示回归向量的p-模，p＝2时就是欧几里模；1≤a，b＜∞，L表示一个正则化算子，它迫使B的估计值属于相应的特性良好的函数的子空间；λ表示控制第2项相对于第1项的权值；min(·)表示最小值求取算子；式(3)中的第一项是最小化式(2)中的随机误差E，而第2项则是为了最小化回归向量B，以获得灵敏度高、信噪比高的特征谱线，以提高标定模型分析结果的准确性；

对于前向选择法，其方法是观察并比对相同浓度的各组分气体的吸收光谱，估计各组分气体吸收光谱的交叠程度，并按照交叠程度将这些光谱图进行排序，首先为交叠度最低的气体选择特征变量，然后依次为每种气体提取特征变量；其步骤如下：

(1)为第一种气体寻找两条或多条谱线的组合，使得其对该种气体的灵敏度最高，而对其它气体的灵敏度尽可能得低；

(2)为第二种气体的寻找两条或多条谱线的组合，使得它对该种气体的灵敏度大于第三种及以后的气体的灵敏度；

(3)为第三种气体的寻找两条或多条谱线的组合，使得它对该种气体的灵敏度大于第四种及以后的气体的灵敏度；

第四种及以后的气体，以此类推。

3.根据权利要求1所述的多模型信息融合方法，其特征在于，各组分气体的分析模型采用多项式偏最小二乘法，或者采用前向神经网络法来建立；每个模型的输入集合中，包含每种气体的一个特征变量；特征变量组合中存在某个空间，所有的特征变量所包含的谱线均没有吸收饱和，即透射率在0.1以上；

采用偏最小二乘法建模的模型时，其形式为：

$y_{\mathrm{im}}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{a}_j{f}_{\mathrm{jm}}+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{p}_{\mathrm{kj}}{f}_{\mathrm{km}}{f}_{\mathrm{jm}}+\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{q}_{\mathrm{kj}}{f}_{\mathrm{km}}{f}_{\mathrm{jm}}^2+e---\left(3\right)$

式中X_m表示第m个样本的光谱值向量；y_im表示第i种气体的第m个估计浓度；f_im＝X_mB_i表示第i种气体的第m个样本的特征变量值；e表示模型的残差；N表示气体的种类数；a_j、p_kj和q_kj是模型的待定系数，其估计公式为：

P_i＝inv(F_iΛF_i ^T)F_i ^TΛY_i(4)

式中P_i＝[a_j，p_kj，q_kj]，表示第i种气体多项式偏最小二乘法模型的系数向量，F_i＝[f₁，f₂，…，f_N]，Λ为用户设定的权值，Y_i＝[y_i1，y_i2，…，y_iM]^T，上标T表示转置算子，M表示样本的数量。

采用前向神经网络法进行分析模型建模时，网络的隐层节点响应函数为S型函数

$s_i=\frac{2}{1+e^{-t_i}}-1,$ $t_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{ji}}{f}_j,(i=\mathrm{1,2},...,S)---\left(5\right)$

式中s_i表示隐层第i个节点的输出，t_i表示隐层第i个节点的总输入，W_ji表示隐层第i个节点与输入层第j个节点的连接权值，隐层一共有S个节点，S＞N；输出层节点响应函数为线性函数：

y_i＝r_i+d_i， $r_i=\underset{j=1}{\overset{S}{\mathrm{\Σ}}}{V}_{\mathrm{ji}}{s}_j,(i=\mathrm{1,2},...,N)---\left(6\right)$

式中r_i表示输出层第i个节点的总输入，V_ji表示隐层第j个节点与输出层第i个节点的连接权值，y_i表示输出层第i个节点的输出，d_i表示输出层第i个节点截长，网络的参数W_ji、V_ji和d_i学习算法采用动量法：

W_ji(k+1)＝W_ji(k)-α(i)[(1-λ)ΔW_ji(k)+λΔW_ji(k-1)](7a)

V_ji(k+1)＝V_ji(k)-β(i)[(1-λ)ΔV_ji(k)+λΔV_ji(k-1)](7b)

d_i(k+1)＝d_i(k)-γ(i)[(1-λ)Δd_i(k)+λΔd_i(k-1)](7c)

式中α(i)、β(i)和γ(i)分别为气体权值W_ji、V_ji和截长d_i的学习率，λ为动量因子，且0≤λ≤1；对于ΔW_ji、ΔV_ji和Δd_i分别用下式：

$\mathrm{\Δ}{W}_{\mathrm{ji}}=\frac{g(W_{j1},...,W_{\mathrm{ji}}+e,...,W_{\mathrm{jS}},V_{\mathrm{ji}},d_i)-g(W_{j1},...,W_{\mathrm{ji}},...,W_{\mathrm{jS}},V_{\mathrm{ji}},d_i)}{e}---\left(8a\right)$

$\mathrm{\Δ}{V}_{\mathrm{ji}}=\frac{g(V_{j1},...,V_{\mathrm{ji}}+e,...,V_{\mathrm{jN}},W_{\mathrm{ji}},d_i)-g(V_{j1},...,V_{\mathrm{ji}},...,V_{\mathrm{jN}},W_{\mathrm{ji}},d_i)}{e}---\left(8b\right)$

$\mathrm{\Δ}{d}_{\mathrm{ji}}=\frac{g(W_{\mathrm{ji}},V_{\mathrm{ji}},d_1,...,d_i+e,...,d_N)-g(W_{\mathrm{ji}},V_{\mathrm{ji}},d_1,...,d_i,...,d_N)}{e}---\left(8c\right)$

式中e是一个足够小的量，可取e＝0.001，函数f(·)是由附图3的前向神经网络结构在式(5)和(6)响应函数情况下形成的函数。

4.根据权利要求1所述的多模型信息融合方法，其特征在于，各组分气体的每个模型的信噪比采用如下步骤估计：

1)估计每条谱线的噪声水平；在气室中充满氮气，并连续扫描光谱8～15次左右，求每条谱线的均方根；对于以透射率为输出的光谱，按式(9)计算：

$\mathrm{\Δ}{x}_i=\sqrt{\frac{\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{(x_{\mathrm{ij}}-{\stackrel{\‾}{x}}_i)}^2}{L-1}---\left(9\right)}$

式中L表示扫描的次数，Δx_i表示第i条谱线的均方差，x_ij表示第j次扫描第i条谱线的谱线值，

表示L次光谱扫描结果中第i条谱线的均值；而对于以吸光度为输出的光谱，每条谱线值进行指数运算exp(-v_ij)后代入式(9)进行计算；

2)估计特征变量f_i的噪声水平，对于以透射率和吸光度为输出的光谱图，给定光谱数据X以后，特征变量的噪声分别由式(10a)和(10b)估计，或者分别由式(11a)和(11b)估计：

$\mathrm{\Δ}{f}_i=2\underset{j=1}{\overset{R}{\mathrm{\Σ}}}\left|\right|b_j\left|\right|\mathrm{\Δ}{x}_j/x_j---\left(10a\right)$

$\mathrm{\Δ}{f}_i=2\underset{j=1}{\overset{R}{\mathrm{\Σ}}}\left|\right|b_j\left|\right|\mathrm{\Δ}{x}_j/\exp (-x_j)---\left(10b\right)$

$\mathrm{\Δ}{f}_i=2\sqrt{\underset{j=1}{\overset{R}{\mathrm{\Σ}}}{(b_j\mathrm{\Δ}{x}_j/x_j)}^2}---\left(11a\right)$

$\mathrm{\Δ}{f}_i=2\sqrt{\underset{j=1}{\overset{R}{\mathrm{\Σ}}}{(b_j\mathrm{\Δ}{x}_j/\exp (-x_j))}^2}---\left(11b\right)$

式中b_j是由特征变量提取算法确定的向量B_i的第j个元素，这里R表述B_i的长度；

3)估计每个分析模型的噪声水平；首先估计每组分气体分析模型对每个输入f_i的微分：

$\mathrm{\Δ}{g}_{\mathrm{jki}}=\frac{g_{\mathrm{jk}}(f_1,...,f_i+e,...,f_n)-g_{\mathrm{jk}}(f_1,...,f_i,...,f_N)}{e}---\left(12\right)$

式中g_jk(·)表示由步骤2)的前向神经网络或多项式确定的第j组分气体的第k个分析模型，e是一个很小的量，直接设定为0.0005～0.005之间的数；然后按式(13)进行该模型的噪声的估计：

${\mathrm{\δ}}_{\mathrm{jk}}=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\mathrm{\Δ}{g}_{\mathrm{jki}}\mathrm{\Δ}{f}_i\right)}^2}---\left(13\right)$

式中δ_jk表示第j组分气体的第k个分析模型的噪声水平。

5.根据权利要求1所述的多模型信息融合方法，其特征在于，各组分气体的最终分析结果的信息融合按式(14)进行：

${\tilde{y}}_j=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{y}_{\mathrm{jk}}\frac{{\mathrm{\τ}}_k}{{\mathrm{\κ}}_j}---\left(14\right)$

式中y_jk表示第j组分气体的第k个分析模型的分析结果，

表示第j组分气体的融合结果，K表示第j组分气体的分析模型个数，κ_j和τ_k是由分模型噪声水平估计的系数，由式(15)估计：

${\mathrm{\κ}}_j=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\τ}}_k,$ ${\mathrm{\τ}}_k=\frac{\underset{l=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\δ}}_{\mathrm{jl}}}{{\mathrm{\δ}}_{\mathrm{jk}}}---\left(15\right).$

Images