CN101059428A - 基于修正卡尔曼滤波理论的气体浓度定量分析仪 - Google Patents

Abstract

本发明公开了一种基于修正卡尔曼滤波理论的气体浓度定量分析仪，硬件部分由光源、接收系统、光谱仪、数据传输线及计算机等组成。软件部分的修正卡尔曼滤波算法由标准卡尔曼滤波算法和噪声参数自适应估计算法两部分组成。光源射出的光被被测气体吸收进接收系统后，进入光谱仪；光谱仪获得波长相对应的透射率关系给出光谱数据经过数据传输线进入计算机。进入计算机内的光谱数据建模成测量方程，浓度分布不随时间变化，由此建立气体浓度状态方程，然后利用修正卡尔曼滤波算法进行处理，鲁棒反演气体浓度。

Description

基于修正卡尔曼滤波理论的气体浓度定量分析仪

技术领域

本发明涉及一种气体浓度定量分析装置，尤其是指一种基于修正卡尔曼滤波理论的气体浓度定量分析仪。

背景技术

环境污染是世界范围普遍关心的问题。对各类燃烧器、工业及商用锅炉产生的气体、机动车尾气、各类化工厂有毒气体(如HCl、NH3、HF)等进行环境污染物的监测和预报，是环境保护、控制治理污染的必要前提。

传统的气体污染监测是以湿式化学技术和吸收取样后的实验分析为基础，这些仪器通常只限于单点测量，响应时间比较慢。相比而言，光学和光谱学遥感技术以其大范围、连续实时监测方式而成为环境气体监测、分析的理想工具。其中差分吸收光谱技术DOAS(Differential OpticalAbsorption Spectomery)以其出众的测试方法及技术特点成为大气污染模式研究和大气污染监测的常用方法之一。

DOAS技术中一般采用最小二乘原理对气体浓度进行反演，但是这种方法的缺点的在反演求解过程中假设系统数据(如压强、温度等)精确已知。而实际情况是系统往往会受到各种因素的扰动，因此造成该反演算法的鲁棒性不够，而且在数据含有大噪声的情况下反演精度也不够理想。为了解决上述问题，针对实际应用中可能遇到的各种扰动因素，把气体浓度随时间变化的方程建模成状态方程，从DOAS测量原理出发推导出测量方程，从而建立浓度气体反演的状态空间表达，随后充分考虑系统的不确定性及其测量数据的噪声问题，采用卡尔曼滤波改进原最小二乘算法的不足。

DOAS是德国Heidelberg大学环境物理研究所的U Platt和D Perner在20世纪70年代末提出的，随后在大气实用探测方面迅速发展。光谱测量基本原理：已知从光源发出的光束经过某一待测气体，由于不同气体分子对光的吸收不同，使得吸收光谱的强度和结构都发生了相应的改变。提取这些特征结构，然后利用最小二乘法，就可以求得相应气体的浓度。这就是DOAS方法的主要思想。光强分布服从Beer-Lambert公式：

       I(λ)＝I₀(λ)exp[-(Lcσ(λ)+ε_R(λ)+ε_M(λ))]+N        (1)

(1)式中，λ为波长，I(λ)为吸收光谱的光强测量值，I₀(λ)为出射光光强初始值，L为光程，c为待测气体浓度，σ(λ)为气体吸收截面，ε_R(λ)为Reileigh散射系数，ε_M(λ)为Mie散射系数，N为各种噪声总和。对(1)式两边取对数得： $\ln \left[\frac{I_0\left(\mathrm{\λ}\right)}{I\left(\mathrm{\λ}\right)}\right]=\mathrm{\σ}\left(\mathrm{\λ}\right)\mathrm{cL}+{\mathrm{\ϵ}}_R+{\mathrm{\ϵ}}_N+N^{\′}---\left(2\right)$

通过光谱仪测量，我们可以测得对应波长λ_mk(mk＝1，2，……，n)的强度I(λ_mk)和I₀(λ_mk)，然后我们采用最小二乘法对这n组超定方程求解浓度c。

最小二乘法只对量测值

和相应的估计值

之间的误差的平方和求最小，而且在求解过程中气体吸收界面数据来自先验的实验室测量，而没有考虑到实际系统环境中温度及其压强的变化，这个是最小二乘方法的局限之一。其次，在噪声较大的情况下，这就会使得浓度估计值大大的偏离了真实浓度值，造成反演精度下降。因此通常的做法是首先对方程(2)中的噪声进行滤除以降低噪声的干扰，然后再用最小二乘法求解，这样便可进一步提高反演的精确度。

现有的气体浓度分析装置都是基于气体光谱吸收技术，一般采用最小二乘法对气体浓度进行求解。其数学模型可以表达为：Y＝HC+v(v代表噪声，C是待求的气体浓度，Y为气体浓度分析仪器测量所得的光谱数据)，至今为止，最小二乘算法在计算过程中把H看成固定不变的。然而实际测量的光谱的漂移、拉伸或压缩或者因为测量条件的变化使得H发生变化，因此把H看成固定不变的可能会带来误差。

由此可见，最小二乘方法比较简便，但是在外界噪声干扰较大的情况下，探测的精度不高。尤其在实际应用中，由于环境因素变化，气体的吸收截面(数据H)跟先验吸收截面的有区别，如果若仍旧采用先验吸收截面数据，那反演的结果误差就较大了。

最佳线性滤波理论起源于40年代美国科学家Wiener和前苏联科学家Kолмогоров等人的研究工作，后人统称为维纳滤波理论。从理论上说，维纳滤波的最大缺点是必须用到无限过去的数据，不适用于实时处理。

为了克服这一缺点，60年代卡尔曼(Kalman)把状态空间模型引入滤波理论，并导出了一套递推估计算法，后人称之为卡尔曼滤波理论。卡尔曼滤波是基于最小均方误差为估计的最佳准则的递推估计的算法，其基本思想是：利用前一时刻的估计值和现时刻的观测值来更新对状态变量的估计，求出最优估计。它适合于实时处理和计算机运算。

而标准卡尔曼滤波模型要求系统的动态噪声和测量噪声为零均值且统计特性已知的白噪声，但在实际的气体浓度测量中，存在着Reileigh散射、Mie散射以及光电探测器噪声等噪声，不可能准确地测得动态噪声的统计特性。为此，应估计测量方程协方差阵和状态误差协方差阵。

自适应估计开窗逼近法是Mehra于1970年提出的，其基本原理是利用前N步历史信息的信息和测量残差来估计当前状态噪声协方差矩阵和测量噪声协方差矩阵，从而使状态协方差矩阵和测量协方差矩阵趋向于实际噪声协方差，并且能自适应于当前状态信息和测量信息，再应用标准卡尔曼滤波模型进行计算，获得系统状态的最优估值。来估计测量方程协方差阵及状态误差协方差阵。

结合标准的卡尔曼滤波和噪声自适应估计算法构成修正卡尔曼滤波理论，可以有效的解决测量噪声或者系统带来的不确定性问题，从而鲁棒反演气体浓度。鲁棒性的概念最初出现在统计决策理论中.用于贝叶斯规则的灵敏度分析，指的是如果决策规则对于某个函数的变化不灵敏，那么我们说它对于这个函数是鲁棒的。从反演的角度看，鲁棒性是指反演过程的稳定性与健壮性，即在受到干扰的情况下，反演结果继续维持正常情况的程度。

发明内容

本发明在气体光谱吸收技术的基础上，从噪声统计和考虑系统不确定性的角度，利用测量值不断对浓度状态值进行估计修正，提供了一种以气体的吸收光谱为测量数据对气体浓度进行定量分析的基于修正卡尔曼滤波理论的气体浓度定量分析仪。

一种基于修正卡尔曼滤波理论的气体浓度定量分析仪，由光源、接收系统、光谱仪和计算机组成，光谱仪接入计算机，所述的计算机基于修正卡尔曼滤波理论的气体浓度分析软件对气体浓度进行定量分析，所述的气

体浓度分析软件包括：

卡尔曼滤波迭代模块，用于结合测量数据对气体浓度进行在线估计

噪声参数自适应估计模块，用于结合测量数据对卡尔曼计算过程中涉及的状态噪声和测量噪声进行在线实时修正；

光源发出光波经被测气体后经接收系统进入光谱仪，光谱仪获得光波中每个波长相对应的透射率关系，给出对应的光的强度光谱数据，发送给计算机，计算机内气体浓度分析软件将光谱数据建模成测量方程，建立气体浓度状态方程，利用修正卡尔曼滤波算法进行处理，鲁棒反演气体浓度。

气体浓度分析软件对所获得的测量数据进行分析的过程具体如下：

(1)建立气体吸收截面数据

吸收截面反映了气体吸收光强能力的大小，它跟波长、温度、压强有关。气体的吸收截面是气体浓度反演的基础，现有的吸收截面数据可以选择哈佛大学的HITRAN公开数据库，它包含了绝大多数常见气体的吸收截面数据，如O₃，CO₂，SO₂，NO₂等，其中也包含了气体在各个温度条件下的吸收截面。结合选用的光谱仪的响应函数得到仪器专用的气体吸收截面数据。

(2)测量光谱预处理

在DOAS中，会获得经气体吸收的测量光谱和未经气体衰减的光源光谱。有时也会获取背景光谱。测量光谱和光源光谱减去背景光谱，消除背景光干扰。然后，用减去背景光谱的测量光谱除以做了同样处理的光源光谱，对其结果取对数，就得到我们所需要的光谱测量数据。

(3)计算机反演

首先假设气体浓度达到稳态的时候，不随时间变化，有：

$\frac{\mathrm{dc}\left({\mathrm{\λ}}_k\right)}{\mathrm{dt}}=0---\left(3\right)$

对其进行离散化，得到浓度状态方程为：

            c_k＝c_k-1+w_k-1                    (4)

其中c_k为λ_k下的气体浓度，w_k为状态噪声，表示对系统建模的不确定性。

根据(2)式对其进行简化，令 $y_k=\ln \left[\frac{I_0\left({\mathrm{\λ}}_{\mathrm{mk}}\right)}{I\left({\mathrm{\λ}}_{\mathrm{mk}}\right)}\right],$ H_k＝σ(λ_mk)L，v_k＝ε_R+ε_M+N从而得到测量方程：

            y_k＝H_kc_k+v_k，k＝1，2，…，n         (5)

其中y_k为测量值，H_k为气体吸收截面矩阵与总光程的乘积，v_k为测量噪声。

则式(4)、(5)构成了气体浓度测试的状态方程，即气体浓度反演的状态空间表达。

假定动态噪声{w_k}和量测噪声{v_k}为无不相关的零均值白噪声序列，并且符合： $\left\{\begin{array}{c}E\left(w_k\right)=0,\mathrm{cov}\left(w_k{w}_k^T\right)=Q_k{\mathrm{\δ}}_{\mathrm{kj}}\\ E\left(v_k\right)0,\mathrm{cov}\left(v_k{v}_k^T\right)=R_k{\mathrm{\δ}}_{\mathrm{kj}}\\ \mathrm{cov}\left(w_k{v}_j^T\right)=0\end{array}\right.---\left(6\right)$

通过利用标准卡尔曼滤波，得如下解：

${\hat{c}}_{k/k-1}={\hat{c}}_{k-1}---\left(7\right)$

${\hat{c}}_k={\hat{c}}_{k/k-1}+K_k(y_k-H_k{\hat{c}}_{k/k-1})---\left(8\right)$

$K_k=P_{k/k-1}{H}_k^T{[H_k{P}_{k/k-1}{H}_k^T+R_k]}^{-1}---\left(9\right)$

             P_k＝(I-K_kH_k)P_k/k-1          (10)

             P_k/k-1＝P_k-1+Q_k-1           (11)

输入为测量值y_k，假设浓度初始值为

，浓度误差协方差为P₀，迭代从初始值c₀，P₀，Q₀，R₀出发，通过每次得到的测量值y_k，不断修正估计出来的浓度值

，最终给出估计浓度。在稳态条件下，状态方程中的状态矩阵A_k，k-1为单位矩阵。假定实际噪声一般噪声，其模型为： $\left\{\begin{array}{c}E\left(w_k\right)=q_k,E\left((w_k-q_k){(w_k-q_k)}^T\right)=Q_k\\ E\left(v_k\right)=r_k,E\left((v_k-r_k){(v_k-r_k)}^T\right)=R_k\end{array}\right.---\left(12\right)$

设定门限长度N值，在卡尔曼迭代进行N步后，进入噪声参数自适应模块。输入为前N项浓度估计值和测量值，输出为噪声的均值和协方差。

状态信息为： $L_k={\hat{c}}_k-A_{k,k-1}{\hat{c}}_{k-1}---\left(13\right)$

测量信息为： $M_k=H_k{\hat{c}}_k-y_k---\left(14\right)$

利用前N次的历史信息，则可得噪声的协方差：

${\hat{q}}_k=\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{L}_{k-i}---\left(15\right)$

$Q_k=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(L_{k-1}-{\hat{q}}_k){(L_{k-1}-{\hat{q}}_k)}^T---\left(16\right)$

${\hat{r}}_k=\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{M}_{k-1}---\left(17\right)$

$R_k=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(M_{k-i}-{\hat{r}}_k){(M_{k-1}-{\hat{r}}_k)}^T+\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H}_{k-1}{P}_{k-1}{H}_{k-1}^T---\left(18\right)$

将方程(13)～(18)结合到方程(7)～(11)中，就可对未知噪声进行修正卡尔曼滤波计算。

所述的光源选用光谱范围宽的氙弧灯、半导体红外激光器等。

所述的接收系统选用望远系统、光纤等。

所述的光谱仪选用紫外光谱仪、红外光谱仪或其他波段的光谱仪。

所述的光谱仪和计算机通过数据传输线连接，数据传输线选用RS-232串口、485总线、USB形式或其它用于数据传输的连接线。

所述的计算机选用个人计算机、笔记本计算机、嵌入式计算机或工控计算机。

本发明的优点为：

鲁棒性高，能消除由于环境因素变化带来的影响，提高结果的可重复性；

测量精度高，由于卡尔曼本身具有滤波和估计的双重特点，使得估计精度提高。

由于卡尔曼滤波的基本方程是时间域内的递推形式，其计算过程是一个不断地“预测一修正”过程，在求解时不要求存储大量数据，并且一旦观测到了新的数据，随时可以算得新的滤波值，因此非常便于实时处理和计算机实现。

附图说明

图1为本发明装置的硬件结构框图；

图2为本发明装置气体浓度分析软件的操作流程图。

图3为本发明装置所测试的气体一组参考光谱和测量光谱图

具体实施方式

如图1所示，一种基于卡尔曼滤波理论的气体浓度定量分析仪，包括光源、接收系统、光谱仪、计算机以及连接光谱仪和计算机的数据传输线。

光源可以根据实际情况确定，可以采用光谱范围宽的氙弧灯，也可以采用半导体红外激光器。

接收系统可以根据实际情况确定，可以采用望远系统。

光谱仪，用于测量光波中每个波长对应的光的强度数据，可以是紫外光谱仪，也可以是红外光谱仪，也可以是其他波段的光谱仪。光谱仪的基本参数包括波长间隔、扫描波长范围、扫描速度等参数，可以根据实际情况进行选用。

数据传输线，可以采用RS-232串口形式，485总线、USB形式或者其他能与计算机进行数据传输的信号形式与计算机相连。

计算机可以是个人计算机、笔记本计算机、嵌入式计算机、工控计算机等计算机。

气体浓度分析软件对测量数据的分析处理过程如图2所示：

(1)先验气体吸收截面数据库建立

吸收截面反映了气体吸收光强能力的大小，它跟波长、温度、压强有关。气体的吸收截面是气体浓度反演的基础，现有的吸收截面数据可以选择哈佛大学的HITRAN公开数据库，它包含了绝大多数常见气体的吸收截面数据，如O₃，CO₂，SO₂，NO₂等，其中也包含了气体在各个温度条件下的吸收截面；也可以选用自己测定的数据。然后与仪器响应函数卷积，再与光程相乘计算得到该仪器的气体吸收截面数据库。

(2)测量光谱预处理

在DOAS中，获得经气体吸收的测量光谱和未经气体衰减的光源光谱。有时也会获取背景光谱。测量光谱和光源光谱减去背景光谱，消除背景光干扰。然后，用减去背景光谱的测量光谱除以做了同样处理的光源光谱，对其结果取对数，就得到我们所需要的光谱测量数据。

(3)计算机反演

建立状态空间方程，结合经过预处理后的测量数据，通过修正卡尔曼滤波模块(公式(7)～(18))对气体浓度进行反演，输出结果。

实施例：

测量的气体成分为臭氧，使用的光谱仪器为CCD光纤光谱仪，覆盖180nm-800nm的光谱范围，通过计算机采集到的光强数据，再通过修正卡尔曼滤波对气体浓度实现反演。

(1)建立先验气体吸收截面数据库：采用哈佛大学的HITRAN公开数据库做为所需系统参数的先验数据，把气体吸收截面矩阵与仪器函数卷积后再与光程的乘积得到矩阵H值。

(2)实际浓度已知的气体，加入混合气体，形成不同的测量噪声，利用光谱仪得到测量光谱和未经气体衰减的光源光谱。然后，测量光谱除以光源光谱，所得的结果取对数，就得到我们所需要的光谱测量数据，即Y值。图3给出了一组测量数据。

(3)设置任意浓度初始值为

，浓度误差协方差为P₀，噪声协方差矩阵Q₀，R₀，利用公式(7)～(18)进行迭代计算。做为对比，也使用最小二乘算法进行了反演计算。测量中，对于每次实验，都测量了多组数据，然后对于每组数据利用修正卡尔曼滤波和最小二乘方法进行反演，最后比较计算其平均误差、均方误差和最大误差。

定义信噪比为：

$\mathrm{snr}=10{\log }_{10}^{\frac{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(H_j*c)}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_i^2}}$

表1：测量数据含有噪声时，修正卡尔曼滤波和最小二乘算法的对比

信噪比		平均误差	均方误差	最大误差
					-52.17	卡尔曼	-1.001	1.002	1.120
最小二乘算法	-43.775	282.98	673.620
				-35.58		卡尔曼	-0.869	0.969	1.978
最小二乘算法	1.792	22.092	50.858
					-16.13	卡尔曼	0.171	1.133	3.050
最小二乘算法	0.313	2.263	4.853
				7.79		卡尔曼	0.010	0.252	0.555
最小二乘算法	0.024	0.227	0.825

(4)实际浓度已知的气体，在不同的状态环境下，利用光谱仪得到测量光谱和未经气体衰减的光源光谱。然后，测量光谱除以光源光谱，所得的结果取对数，就得到我们所需要的光谱测量数据，即Y值。设置任意浓度初始值为

，浓度误差协方差为P₀，噪声协方差矩阵Q₀，R₀，利用公式(7)～(18)进行迭代计算。做为对比，使用最小二乘算法进行了反演计算。测量中，对于每次实验，都测量了多组数据，然后对于每组数据利用修正卡尔曼滤波和最小二乘方法进行反演，最后比较计算其平均误差、均方误差和最大误差。

表2：外界压强温度变化时，修正卡尔曼滤波和最小二乘算法的对比

信噪比		平均误差	均方误差	最大误差
					-43.38	卡尔曼	-0.971	0.977	1.187
最小二乘算法	10.555	96.783	217.231
				-23.64		卡尔曼	-0.459	0.947	2.230
最小二乘算法	2.501	9.936	23.826
					-3.21	卡尔曼	0.014	0.702	1.922
最小二乘算法	0.033	0.851	2.369
				15.86		卡尔曼	0.008	0.118	0.344
最小二乘算法	0.017	0.107	0.2489

Claims (6)

1.一种基于修正卡尔曼滤波理论的气体浓度定量分析仪，由光源、接收系统、光谱仪和计算机组成，光谱仪接入计算机，其特征在于：所述的计算机基于修正卡尔曼滤波理论的气体浓度分析软件对气体浓度进行定量分析，所述的气体浓度分析软件包括：

卡尔曼滤波迭代模块，用于结合测量数据对气体浓度进行在线估计；

噪声参数自适应估计模块，用于结合测量数据对卡尔曼计算过程中涉及的状态噪声和测量噪声进行在线实时修正；

光源发出光波经被测气体后经接收系统进入光谱仪，光谱仪获得光波中每个波长相对应的透射率关系，给出对应的光的强度光谱数据，发送给计算机，计算机内气体浓度分析软件将光谱数据建模成测量方程，建立气体浓度状态方程，利用修正卡尔曼滤波算法进行处理，鲁棒反演气体浓度。

2.如权利要求1所述的气体浓度定量分析仪，其特征在于：所述的光源选用光谱范围宽的氙弧灯、半导体红外激光器或其他类型的光源。

3.如权利要求1所述的气体浓度定量分析仪，其特征在于：所述的接收系统选用望远系统或光纤。

4.如权利要求1所述的气体浓度定量分析仪，其特征在于：所述的光谱仪选用紫外光谱仪、红外光谱仪或其他波段的光谱仪。

5.如权利要求1所述的气体浓度定量分析仪，其特征在于：所述的光谱仪和计算机通过数据传输线连接，数据传输线选用RS-232串口、485总线、USB形式或其它用于数据传输的连接线。

6.如权利要求1所述的气体浓度定量分析仪，其特征在于：所述的计算机选用个人计算机、笔记本计算机、嵌入式计算机或工控计算机。

Images