CN104297115B - 一种大气颗粒物pm2.5数密度检测的方法 - Google Patents

Abstract

本发明提供了一种大气颗粒物PM2.5数密度的检测方法，其采用激光光源输出的红外光，照射含有PM2.5的气体，因散射效应，形成散射光，用光学部件，将各散射光汇聚，形成干涉光，采用光电探测器等测量干涉光强，为消除噪声污染，用自适应FIR滤波器进行信号处理，并计算干涉光强平均值、方差，依据该平均值、方差与散射光的总光强、形成散射光的PM2.5总数目的相关性原理，按***参数估计算法，完成对散射光强、PM2.5总数目的估计，根据已知测试室容积和PM2.5总数目的估计值，计算被测气体的PM2.5数密度，实现大气PM2.5的检测，其检测原理新颖、独特，其算法基于递推最小二乘RLS原理，无需矩阵求逆运算，运算量低，性能稳定，其***结构简洁，成本低，可靠性高。

Description

一种大气颗粒物PM2.5数密度检测的方法

技术领域

本发明涉及大气颗粒物检测技术，特别涉及一种大气颗粒物PM2.5数密度检测的方法及其装置。

背景技术

PM是颗粒物(Particulate Matter)的英文缩写，用于指代大气中包括固体和液体在内的颗粒物。通常根据空气动力学粒径的不同，PM又分为PM10和PM2.5，分别对应空气动力学粒径10μm及2.5μm以下的大气颗粒物。PM2.5相对于PM10，粒径更小、毒性更大、传输距离更远、停留时间更长；PM2.5也称为可入肺颗粒物，携带各种有害物质，通过呼吸道进入肺泡，经过血液循环到达其他器官，对人体呼吸***和其他功能***的损害极大。

目前，广泛采用的PM2.5测定方法是：滤膜称重法、β射线吸收法、微量振荡天平法。

滤膜称重法：该方法在一定的流量采样,将空气中的颗粒物捕集于高性能滤膜上,称量滤膜采样前后的质量,由其质量差求得捕集的颗粒物质量,再根据采样前后滤膜的重量差与采样空气体积,计算出颗粒物的质量数密度。它的缺点是需要手工进行操作,影响测量精度的因素较多且操作程序繁杂、占用房间和设备较多、采样时间较长、仪器维修量大、花费成本较高等。因该方法自动化程度低,只能定期进行监测,很难用于在线监测。

β射线吸收法：用玻璃纤维滤纸收集大气中的微颗粒子，在颗粒物采样前，首先采用β射线穿过清洁的玻璃纤维滤纸，记下射线衰减值；再让样本气体通过玻璃纤维滤纸，完成颗粒物采集；在颗粒物采集后，再次用β射线穿过同一滤纸；由于β射线穿过有尘滤纸时，射线强度将衰减，其强度的衰减仅与尘的质量有关，根据两次β射线被吸收的变化量，来求得收集在滤纸上的颗粒物的质量。

微量振荡天平法：基于锥形元件振荡微量天平原理，天平的振荡频率会随着滤膜上收集的颗粒物的质量变化发生改变，通过准确测量频率的变化得到采集到的颗粒物的质量，然后根据收集这些颗粒物时采集的样品体积计算得出样品的数密度。

以上三种主流方法，监测精度较好，为各国环保部门更多的采用。采用过滤膜收集PM2.5，过滤膜污浊，需要定期对其清洗、烘干、更换等处理。采样分析过程繁琐，不能实现长期监控。

另外，光散射法基于光散射原理，空气中的微粒在光的照射下会发生散射现象，光散射与微粒大小、光波波长、微粒折射率、颗粒粒径分布和数密度等因素有关。光散射法是通过测量颗粒物受光照射后所发出的散射光信号，来计算颗粒物的粒径分布和质量数密度。近年来，光散射法因具有测量速度快、精度高、重复性好，且适用于在线测量等优点，而受到了广泛关注。

中国发明专利“基于广角傅立叶变换的PM2.5检测装置”(专利申请号：201310519249.9),采用激光器、准直器、负压产生单元、复变透镜、衍射屏和摄像机等部件组成。利用准直器将激光器输出光，准直为所需要的直径光束，颗粒物在激光光束照射下，在复变透镜的后焦面形成衍射，通过衍射屏形成明暗相间的同心圆，该同心圆的太小随颗粒粒径变化而变化。由摄像机获得同心圆图像，通过对该同心圆图像进行处理，可以得到颗粒物的大小值。因此，该发明仅涉及颗粒物直径的检测方法，且***比较复杂。

发明内容

为了解决现有技术中问题，本发明提供了一种大气颗粒物PM2.5数密度检测方法，其包括以下步骤：

步骤1：含有多种粒径颗粒物的被测空气，经过不同孔径的过滤器分离后，形成一股含有PM2.5及更小颗粒物的气流，并进入气体测试室；

步骤2：在气体测试室中，用激光光源输出的红外光照射被测气流，一部分红外光与气流中颗粒物发生散射，形成散射光，另一部分红外光不与颗粒物发生散射，直接穿透被测气流，形成穿透光，采用光学部件，将散射光和穿透光汇聚在一起；

步骤3：汇聚在一起的散射光和穿透光，发生干涉效应，并输出干涉光；

步骤4：采用光电探测器测量以上干涉光的强度信号，其包含各散射光的自相关强度、穿透光的自相关强度、散射光与穿透光之间的互相关强度、各散射光之间的互相关强度四项内容，该干涉光的强度信号是随机信号，其概率分布密度函数与形成散射光的颗粒子总数目N相关；

步骤5：用A/D转换器对以上的干涉光总强度信号，以采样周期为T，进行离散化抽样，得到抽样信号序列y(n)，用y(n)表示在nT时刻的抽样信号，n为正整数；

步骤6：从被噪声污染的抽样信号中，恢复出原始的干涉光强信号，采

用递推最小二乘RLS算法，设计一个M阶的自适应FIR滤波器，对抽样

信号y(n)，按下式进行滤波：

其中，w_i(n)为nT时刻的M阶FIR滤波器滤波系数，i＝0,1,...,M-1，W(n)

为nT时刻滤波系数向量：W(n)＝[w₀(n),...,w_M-1(n)]^T，而x(n)为滤波后的信

号序列，x(n)可视为在nT时刻的原始的干涉光强信号估计值，Y(n)为nT

时刻所对应的信号向量，该向量由M个抽样信号组成，其如下：

Y(n)＝[y(n),y(n-1)...,y(n-M+1)]^T

W^T(n)、Y^T(n)分别为W(n)和Y(n)的转置：

W^T(n)＝[w₀(n),...,w_M-1(n)]，Y^T(n)＝[y(n),y(n-1),...,y(n-M+1)]；

步骤7：设滤波后的信号序列x(n)的序号指针为n₀，从n₀＝0开始，连续取最前面的二段长度为L的信号，L为正整数，用信号向量X₁、X₂表示：

X₁＝[x(1),...,x(L)]^T

X₂＝[x(L+1),...,x(2L)]^T

X₁、X₂分别对应在[0-LT]、[LT-2LT]时间内的干涉光强度信号， 分别是信号X₁、X₂的转置形式：

X₁ ^T＝[x(1),...,x(L)]

X₂ ^T＝[x(L+1),...,x(2L)]

按下式，计算X₁、X₂的平均值：

L,M均为正整数，且L＞M；

步骤8：按下式，计算X₁、X₂的方差：

σ₁ ²和σ₂ ²分别对应在[0-LT]和[LT-2LT]时间内的干涉光强度的涨落:

步骤9：根据干涉光强度在[0-LT]、[LT-2LT]时间内的平均值和方差：σ₁ ²、σ₂ ²，建立以下的参数C₁、C₂的矩阵方程：

其中，C₁＝a²,

且a²为粒子散射光的总光强，N为形成散射光的粒子总数目，将C₁、C₂称为***参数，这里假设，a²和N在[0-LT]和[LT-2LT]时间内保持不变，因此，在此时间内，***参数C₁、C₂近似是一组常量；

步骤10：通过以上方程的最小二乘解，得到[0-LT]和[LT-2LT]时间内的***参数C₁和C₂的估计值：

其中，P₀＝(A₀ ^TA₀)^-1，A₀ ^T为矩阵A₀的转置，(A₀ ^TA₀)^-1为A₀ ^TA₀的逆矩阵，令将在[0-LT]和[LT-2LT]时间内的参数C₁、C₂，分别用表示；

步骤11：调整信号序列x(n)的序号指针n₀，令n₀＝n₀+2L，使序号指针，指向下一段长度为L的滤波后的信号；

步骤12：从当前的序号指针n₀开始，取一段长度为L滤波后的信号：X₃＝[x(n₀+1),...,x(n₀+L)]^T，X₃对应于在[n₀T-(n₀+L)T]时间内的干涉光强度信号，按下式，计算X₃的平均值及方差σ₃ ²：

用新生成的平均值及方差σ₃ ²，构成信号更新向量d^T及b⁽¹⁾：

b⁽¹⁾＝1，

建立在[n₀T-(n₀+L)T]时间内的参数C₁、C₂的矩阵方程：

b⁽¹⁾＝1；

步骤13：根据所构成的信号更新向量d^T及b⁽¹⁾，采用递推最小二乘RLS算法，按下式，对在[n₀T-(n₀+L)T]时间内的对***参数C₁、C₂的进行估计，并分别用表示之：

其中，向量d与d^T互为转置关系，与分别对应于本次迭代运算以前与之后的***参数C₁、C₂；

步骤14：根据以上[n₀T-(n₀+L)T]时间内的***参数C₁和C₂，按下式，计算气体测试室在该时间内，散射光的总光强a²及形成散射光的粒子总数目N:

a²＝C₁,

步骤15：根据已知的气体测试室的容积V和以上求得的粒子总数目N，能计算出在[n₀T-(n₀+L)T]时间内，被测气体的颗粒物数密度ρ，即单位体积的气体中PM2.5粒子数量：

步骤16：根据递推最小二乘RLS算法要求，按下式，重置迭代运算各

单元的初值：

步骤17：调整信号序列X的序号指针n₀，令n₀＝n₀+L，使序号指针指向下一段长度为L的滤波后的信号；

步骤18：重复步骤12-17，如此循环，完成在起始时刻n₀＝0之后的，任意时间[n₀T-(n₀+L)T]内的***参数C₁、C₂的估计，进而，得到任意时间[n₀T-(n₀+L)T]内的散射光的总光强a²及形成散射光的粒子总数目N，实现被测气体在起始时刻n₀＝0之后的，任意时间[n₀T-(n₀+L)T]内的颗粒物数密度的检测，其中，n₀＝0,2L,3L,4L,5L,.......。

作为本发明的进一步改进，步骤1中，被测空气在负压发生单元的作用下，经过二级不同孔径的过滤器分离。

作为本发明的进一步改进，步骤2中，所述光学部件为凸透镜或凸透镜组。

作为本发明的进一步改进，步骤6中，自适应FIR滤波器最小二乘RLS算法，其特征在于，其包括以下步骤：

步骤1：初始化：

1)定义n＝-1时，W^T(-1)＝[0,0,...,0]^T，R(-1)＝δ₁ eye，即R^-1(-1)＝δ₁ ^-1 eye，δ₁是很小的正数，eye为单位对角矩阵；

2)定义n＝0时的抽样信号向量Y(0)：Y(0)＝[y(0),0..,0]，Y(0)为M维向量；

3)定义n＝0时的原始信号估计值d(0)，d(0)＝δ₂，δ₂是很小的正数；

4)选常数λ：0＜λ＜1。

步骤2：由Y(0)、R^-1(-1)及λ，按下式,计算n＝0时的向量K(0)：

步骤3：按下式,计算n＝0时的自适应FIR滤波系数向量W(0)：

e(0)＝d(0)-W^T(-1)Y(0)＝d(0)

W(0)＝W(-1)+K(0)e(0)＝K(0)e(0)

步骤4：按下式，计算n＝0时的矩阵R^-1(0)：

步骤5：令n＝n+1，利用当前nT时刻的抽样信号y(n)，更新信号向量Y(n)：

Y(n)＝[y(n),y(n-1),...,y(n-M+1)]^T

步骤6：由当前nT时刻信号向量Y(n)和(n-1)T时刻矩阵R^-1(n-1)，按下式，计算当前nT时刻的向量K(n)、矩阵R^-1(n)：

步骤7：由以下递推公式，计算当前nT时刻的滤波系数向量W(n)：

d(n)＝W^T(n-1)Y(n-1)

e(n)＝d(n)-W^T(n-1)Y(n)

W(n)＝W(n-1)+K(n)e(n)

并对抽样信号y(n)，按下式进行自适应FIR滤波：

步骤8：重复以上步骤5-7，如此循环，获得任意时刻的滤波后的信号序列x(n)。

一种大气颗粒物PM2.5数密度检测装置，其包括激光光源、光隔离器、光电探测器、A/D转换器、信号处理单元、控制单元、气体测试室、进气罩、负压发生单元，所述的激光光源依次通过第一光纤、光隔离器、第二光纤与气体测试室相连，所述的光电探测器通过第三光纤与气体测试室相连，光电探测器的电信号输出端与A/D转换器相连，A/D转换器与信号处理单元相连，所述控制单元分别与信号处理单元、负压发生单元相连，所述进气罩通过进气管与气体测试室相连，所述的负压发生单元通过出气管与气体测试室相连。

作为本发明的进一步改进，进气罩内设置两级不同口径的过滤器，第一级过滤器能滤除较大尺寸的颗粒，第二级滤除粒径大于2.5微米的颗粒。

作为本发明的进一步改进，所述的气体测试室内置两组光纤式自聚焦透镜、凸透镜，其中，第一组自聚焦透镜通过第二光纤、光隔离器与激光光源相连，第一组自聚焦透镜、第一凸透镜负责发射光子，第二组自聚焦透镜通过第三光纤与光电探测器相连，第二组自聚焦透镜、第二凸透镜负责接收光子，该两组自聚焦透镜、凸透镜的中心位于同一光轴上，第一组自聚焦透镜的端面，位于第一凸透镜的前焦点处，第二组自聚焦透镜的端面，位于第二凸透镜的后焦点处，第一凸透镜、第二凸透镜分别起到光线扩束、光线汇聚的作用；

作为本发明的进一步改进，进气管与出气管之间的连线方向，与气体测试室内置两组光纤式自聚焦透镜、凸透镜的光轴方向相垂直，以保证被测气体的流动方向与光子发射方向相垂直。

作为本发明的进一步改进，进气管与出气管之间的连线，位于第一凸透镜与第二凸透镜之间，以保证被测气体的流动区域，主要集中在该两组凸透镜之间的区域。

有益效果如下：

本发明提供了一种大气颗粒物PM2.5数密度检测的方法及其装置，其采用激光光源输出的红外光，照射含有PM2.5的气体，一部分红外光与气流中PM2.5颗粒物发生散射，形成散射光。另一部分红外光不与颗粒物发生散射，直接穿透被测气流，形成穿透光。采用光学部件，将散射光和穿透光汇聚在一起。汇聚在一起的散射光和穿透光，发生干涉效应，并输出干涉光。采用光电探测器、A/D转换器等测量以上干涉光的强度信号，其包含各散射光的自相关强度、穿透光的自相关强度、散射光与穿透光之间的互相关强度、各散射光之间的互相关强度四项内容。该干涉光强度信号是一随机信号，其概率分布密度函数与PM2.5气室的颗粒物数目N相关。关于本发明基本原理的理论推导表明(参见，后面的具体实施方式部分)：该随机信号的平均值、方差与散射光的总光强、形成散射光的PM2.5总数目是相关的。该结论是本发明的基本原理所在。

依据本发明的基本原理，首先用自适应FIR滤波器，对测量数据进行滤波处理，从被噪声污染的抽样信号中，恢复出原始的干涉光强信号滤波。即用滤波后的测量数据，作为原始干涉光强信号的估计。该自适应FIR滤波器，基于递推最小二乘RLS算法。算法的核心是，以指数加权的估计误差平方和构成目标函数，通过该目标函数最小化，来确定FIR滤波系数向量：W(n)＝[w₀(n),...,w_M-1(n)]^T。该系数向量是由自适应算法实现更新和时变。自适应滤波器的最重要特征就在于它能够在未知环境中自动地、有效工作，并能够跟踪输入信号的时变特征，适应***环境的变化，以获得期望响应，使目标函数最小化，完成对未知的原始信号的最优估计。当被测气体的PM2.5颗粒物浓度变化或外界干扰的情形下，其干涉光强度信号为非平稳或时变状态的随机信号。

而传统的FIR滤波器，其滤波系数为固定常数，适合处理平稳环境下的随机信号，无法处理非平稳环境或时变状态的随机信号。

本发明在完成测量数据的滤波处理后，并将滤波后的测量数据，按时间分段，由滤波后的测量数据(即干涉光强信号)，计算出每个时间段所对应的干涉光强平均值、方差。再依据干涉光强的平均值、方差与散射光的总光强、形成散射光的PM2.5总数目是相关性原理，采用***参数估计算法，对该时间段内的散射光的总光强、PM2.5总数目进行估计，并根据PM2.5总数目的估计值和已知的测试室容积，计算出被测气体的该时间段所对应的PM2.5数密度，即单位体积的气体中PM2.5粒子数量。从而实现大气颗粒物PM2.5数密度的检测。

该***参数估计算法也是基于递推最小二乘RLS原理，其算法的核心是，以估计误差平方和构成目标函数，通过该目标函数最小化，来确定***参数C₁和C₂。首先利用最初始的两个时间段对应的干涉光强的平均值和方差σ₁ ²、σ₂ ²，作为信号向量，按最小二乘原理，得到***参数C₁和C₂的估计初值；利用下一个时间段对应的干涉光强的平均值及方差σ₃ ²，对信号向量进行更新；按由最小二乘原理导出的递推公式，在估计初值的基础上，计算出新时间段所对应的***参数C₁和C₂。如此，逐段逐段地完成从初始之后的任意时间段内的参数C₁和C₂估计。

再根据参数C₁和C₂的定义式：C₁＝a²,

计算出从初始之后的任意时间段内的散射光的总光强a²和PM2.5总数目N的估计值。

因此，该参数估计算法，其每一个时间段的参数估计值，是由当前时间段的干涉光强的平均值、方差及前一时间参数估计值，按递推公式来共同决定。其算法无需矩阵求逆运算，具有运算量低，性能稳定的特点。

本发明采用的光电探测器，是光纤式点探测元件，属非面成像元件。其成本低廉，无需衍射屏、CCD摄像机等部件。相对传统技术，***结构简洁，成本低，可靠性更高。

本发明的气体测试室内置两组光纤式自聚焦透镜、凸透镜。其中，第一组自聚焦透镜、凸透镜，负责发射光子和光线扩束；第二组自聚焦透镜、凸透镜，负责光线汇聚和接收光子。

来自激光光源的红外光，经光纤进入第一组自聚焦透镜后，由该自聚焦透镜的端面出光，该端面的出光与光轴成一夹角，其全部光线照射在第一凸透镜上，因该端面正好位于第一凸透镜的前焦点处，光线经第一凸透镜后，形成扩束的平行于光轴的平行光。

由于进气管与出气管之间的连线方向，与气体测试室内置两组光纤式自聚焦透镜、凸透镜的光轴方向相垂直，能使该扩束的平行光与被测气体的流动方向相垂直。同时，由于进气管与出气管之间的连线，位于第一凸透镜与第二凸透镜之间，使被测气体的流动区域，主要集中在该两组凸透镜之间的区域。该平行光在有效地照射被测气体后，因散射效应，形成的散射光线，能有效地到达第二凸透镜上，第二凸透镜将该散射光线，汇聚于其后焦点处。由于该第二凸透镜的后焦点位置，与第二组自聚焦透镜的端面重合。这样，汇聚之后的散射光线，由第二组自聚焦透镜的端面耦合，进入第三光纤，并经该第三光纤被光电探测器接收。

上述设计原理与结构，既能保证红外光充分照射被测气体，也能保证红外光与被测气体所有PM2.5粒子发生散射，其散射效应更充分。同时，保证散射光线，被光电探测器有效的接收。对提高***的测量性能有利。

附图说明

图1是本发明一种大气颗粒物PM2.5数密度检测装置结构示意图。

具体实施方式

下面结合附图对本发明做进一步说明。

一种大气颗粒物PM2.5数密度检测装置，其包括激光光源、光隔离器、光电探测器、A/D转换器、信号处理单元、控制单元、气体测试室、进气罩、负压发生单元，所述的激光光源依次通过第一光纤、光隔离器、第二光纤与气体测试室相连，所述的光电探测器通过第三光纤与气体测试室相连，光电探测器的电信号输出端与A/D转换器相连，A/D转换器与信号处理单元相连，所述控制单元分别与信号处理单元、负压发生单元相连，所述进气罩通过进气管与气体测试室相连，进气罩内设置两级不同口径的过滤器，第一级过滤器能滤除较大尺寸的颗粒，第二级滤除粒径大于2.5微米的颗粒。负压发生单元通过出气管与气体测试室相连。

气体测试室内置两组光纤式自聚焦透镜、凸透镜，其中，第一组自聚焦透镜通过第二光纤、光隔离器与激光光源相连，第一组自聚焦透镜、第一凸透镜负责发射光子，第二组自聚焦透镜通过第三光纤与光电探测器相连，第二组自聚焦透镜、第二凸透镜负责接收光子，该两组自聚焦透镜、凸透镜的中心位于同一光轴上，第一组自聚焦透镜的端面，位于第一凸透镜的前焦点处，第二组自聚焦透镜的端面，位于第二凸透镜的后焦点处，第一凸透镜、第二凸透镜分别起到光线扩束、光线汇聚的作用。

1、本发明的大气颗粒物数密度检测方法的基本原理

含有多种粒径颗粒物的被测空气，在负压发生单元的作用下，经过二级不同孔径的过滤器分离后，形成一股含有PM2.5及更小颗粒物的气流，并进入气体测试室。在气体测试室中，用激光光源输出的红外光照射被测气流，一部分红外光与气流中颗粒物发生散射，形成所谓的散射光。另一部分红外光不与颗粒物发生散射，直接穿透被测气流，形成所谓的穿透光。采用光学部件，将散射光和穿透光汇聚在一起。汇聚在一起的散射光和穿透光，发生干涉效应，并输出干涉光。

如果假设气体测试室中的当前颗粒物数目为N，忽略散射过程的能量损耗，且假设每个粒子形成散射光的偏振方向一致，幅值也一样的，所有粒子散射光的相位是在(-π,π)均匀分布，并且相互独立的。则N个粒子形成的总散射光的电场为

其中，为第i个粒子的散射光电场向量，忽略其向量特性，用标量表示，为第i个粒子的散射光电场幅度，φ_i为第i个粒子的散射光电场相位。

这里，对每个粒子的散射光电场向量E_i，引入标度因子是为了在粒子数目趋于很大时，其所有散射光的能量总和保持有限。即：

a²为粒子散射光的总光强。

用表示透光的电场向量，用E_T和φ_T表示穿透光电场向量的幅度和相位。当散射光和穿透光汇聚在一起的，发生干涉效应，形成干涉光，该干涉光的总电场向量

如果假设散射光与穿透光的偏振方向一致，则干涉光的总电场写成标量形式：

采用光电探测器测量干涉光的总强度I，由于光电探测器是平方律器件，故I＝E_I*E_I ^* (3)

其中，E_I ^*为干涉光的总电场E_I的复共轭。将式(2)代入(3)：

展开上式，得到(4)式:

因为

将式(5)代入(4)：

φ_T-φ_i穿透光与散射光的相位差，φ_i-φ_k为不同粒子的散射光之间的相位差。这里a²、E_T、φ_T可以近似为常量。

从式(6)可以看出：

第一项是穿透光的自相关强度，第二项为散射光的自相关强度；第三项是散射光与穿透光的互相关强度、第四项是各散射光之间的互相关强度；干涉光的总强度I的统计平均值

E()表示对括号符()内的随机变量做统计平均。

假设所有粒子散射光的相位φ_i是在(-π,π)均匀分布，并且相互独立的。且a²、E_T、φ_T近似为常量，故：

当i≠k:

将式(8a)、(8b)代入(7)，得I的统计平均值：

利用式(6)和(9)，计算

根据式(10)，对进行统计平均计算，得到干涉光的总强度I的方差σ²：

对式(10)做统计平均，得到：

第一项：

因为当i≠k时:

E[cos(φ_T-φ_i)cos(φ_T-φ_l)]＝E[cos(φ_T-φ_i)]*E[cos(φ_T-φ_l)]＝0；

故，第一项：

第二项：

第三项：

第四项：

因为当l≠i,k≠i,l≠m,m≠k时：

当l＝i,m＝k,k≠i,l≠m时：

故，第四项：

将(12)、(13)、(14)和(15)代入(11)，得到干涉光的总强度I的方差σ²：

将前面的式(9)：代入式(16)，故方差σ²写成：

(17)式表明，干涉光的总强度I的方差σ²、平均值与散射光的总光强、形成散射光的PM2.5总数目是相关的。

因为干涉光的总强度I为测量数据，I的方差σ²和平均值可由测量数据计算出来。因此，在获得方差σ²和平均值后，可利用(17)式，对形成散射光的粒子总数N和散射光的总光强a²进行估计。再根据已知测量气室的容积，可计算出被测气流的颗粒物数密度，即单位体积的气体中所含颗粒数量。

以上结论，即为本发明的大气颗粒物数密度检测方法的基本原理。

2、本发明的大气颗粒物检测的参数估计算法

用A/D转换器对以上的干涉光总强度信号，以采样周期为T，进行离散化抽样，得的抽样信号序列Y，用信号向量表示：Y＝[y(1),...,y(n),...]^T,Y^T为信号Y的转置形式，其中y(n)是在nT时刻的抽样信号，n正整数。

为了消除外界噪声的影响，从被噪声污染的抽样信号中，恢复出原始的干涉光强信号，采用一个自适应的M阶FIR滤波器，对测量信号Y进行滤波处理：滤波后的信号为X：X＝[x(1),...,x(n),...]^T；

设起始时刻为n₀T，从n₀＝0开始，连续取最前面的二段长度为L的信号：

X₁＝[x(1),...,x(L)]^T

X₂＝[x(L+1),...,x(2L)]^T

X₁、X₂分别对应在[0-LT]、[LT-2LT]时间内的干涉光强度信号，按下式，计算X₁,X₂的平均值、方差：

n₀,L均为正整数；

由(17)式：

其中a₁ ²、N₁和a₂ ²、N₂分别对应在[0-LT]和[LT-2LT]时间的散射光的总功率、形成散射光的粒子总数目。

将(18)、(19)式等号两边，分别除以σ₁ ²、σ₂ ²：

如果假设，在两相邻时间段内的散射光的总功率、形成散射光的粒子总数近似相等，即：

a₁ ²＝a₂ ²＝a²、N₁＝N₂＝N

则(20)式写成矩阵方程：

令

定义参数C₁和C₂：

则(21)式可写成：

用分别表示为在[0-LT]和[LT-2LT]时间内的参数C₁、C₂：

通过以上方程的最小二乘解，得到

其中，A₀ ^T为矩阵A₀的转置，P₀＝(A₀ ^TA₀)^-1,(A₀ ^TA₀)^-1为A₀ ^TA₀的逆矩阵。

通过调整信号序列X的序号指针n₀，令n₀＝n₀+2L，使序号指针指向下一段长度为L的滤波后的信号；从当前的序号指针n₀开始，取一段长度为L滤波后的信号：X₃＝[x(n₀+1),...,x(n₀+L)]^T。X₃对应于在[n₀T-(n₀+L)T]时间内的干涉光强度信号，按下式，并计算X₃的平均值及方差σ₃ ²：

同理，根据(17)式：

其中，a₃ ²、N₃分别对应在[n₀T-(n₀+L)T]时间内的散射光的总功率、形成散射光的粒子总数目。

将(28)式等号两边，分别除以σ₃ ²：

将(28)式写成矩阵形式：

定义信号更新向量d^T及b⁽¹⁾：

相对于前面的数据向量d^T及b⁽¹⁾来自最新的采样数据，因此，将d^T及b⁽¹⁾称为信号更新向量。

根据参数C₁和C₂的定义，式(23)，知道在[n₀T-(n₀+L)T]时间内C₁和C₂：

则(29)可写成：

因为矩阵dd^T：

dd^T为奇异方阵，不可求逆，无法按式(25)，求式(32)的最小二乘解。故将式(32)与式(25)组合成下列矩阵方程：

将在[n₀T-(n₀+L)T]时间内的参数C₁、C₂，分别用表示之，通过求解(33)的最小二乘解：

令

故(34)式写成：

由式(25)，有：

将上式等号两边，同乘以P₀ ^-1：

将式(37)代入式(36)：

根据矩阵求逆引理：

[A+UV]^-1＝A^-1-A^-1U[eye+VA^-1U]^-1VA^-1 (39)

其中，eye为单位对角矩阵，矩阵A为非奇异方阵，且矩阵U,V满足，使eye+VA^-1U为非奇异矩阵。

利用以上矩阵求逆引理，并代入式(35)：

利用P₀＝[A₀ ^T A₀]^-1，有：

综合以上(41)式、(38)式，构成递推式最小二乘算法(Recursive Least-SquaresAlgorithm,RLS)：

综上所述，按RLS递推式最小二乘算法，实现参数C₁和C₂估计的过程总结如下：

一、根据最初始的两段测量数据X₁,X₂，计算其平均值和方差σ₁ ²,σ₂ ²；

二、计算矩阵A₀,P₀,b⁽⁰⁾

P₀＝(A₀ ^TA₀)^-1；

三、计算X₁,X₂对应的参数C₁和C₂：

四、利用新到来的一段测量数据X₃，计算其平均值和方差σ₃ ²，并计算更新向量d^T及b⁽¹⁾：

b⁽¹⁾＝1

五、计算新到的测量数据X₃对应的参数C₁和C₂：

六、按下式，重置迭代运算各单元的初值：

七、连续获取更新的测量数据X₃，重复以上四至六，进行新一轮迭代，逐段递推，直至实现所有测量数据对应的参数C₁和C₂的估计。

八、根据参数C₁和C₂，按下式，计算PM2.5气室在该时间内，散射光的总光强a²及形成散射光的粒子总数目N:

a²＝C₁,

九、根据已知的气体测试室的容积V和粒子总数目N，进而计算出被测气体的颗粒物数密度ρ，即单位体积的气体中PM2.5粒子总数量：

3、自适应FIR滤器的递推最小二乘算法

传统的FIR滤波器，其滤波系数为固定常数，适合处理平稳环境下的随机信号。处理非平稳或时变情形下的随机信号，常采用所谓的自适应FIR滤波器。本发明涉及的自适应FIR滤波器，基于递推最小二乘算法(Recursive Least-Squares Algorithm,RLS)。

假设y(n)是在nT时刻的抽样信号，n为正整数，y(n)是外界噪声污染后的信号。为从被噪声污染的抽样信号中，恢复出的原始信号，采用FIR滤波器，对y(n)进行滤波处理：

x(n)为nT时刻滤波后的信号。

引入nT时刻滤波系数向量W(n)：

W(n)＝[w₀(n),...,w_M-1(n)]^T (44-1)

引入nT时刻抽样信号向量Y(n)：

则(43)是写成：

x(n)＝W^T(n)Y(n)＝Y^T(n)W(n) (44-3)

其中，W^T(n)和Y^T(n)Y^T分别为W(n)和Y(n)的转置：

W^T(n)＝[w₀(n),...,w_M-1(n)],Y^T(n)＝[y(n),...,y(n-M+1)]

如果设d(n)为nT时刻的未受污染的原始信号的预期值，把滤波后的信号x(n)视为nT时刻d(n)的估计值。将原始信号d(n)与其估计x(n)之间的差值ε(n)，定义为估计误差：

ε(n)＝d(n)-x(n)＝d(n)-W^T(n)Y(n)＝d(n)-Y^T(n)W(n) (45)

自适应FIR滤器的递推最小二乘算法的核心是，以指数加权的估计误差平方和构成目标函数，通过该目标函数最小化，来确定FIR滤波系数向量W(n)：

式中，指数加权因子0＜λ＜1称作遗忘因子，其作用是对离n时刻越近的误差加比较大的权重，而对离n时刻越远的误差加比较小的权重，即赋予老数据与新数据以不同的权重值。

最小化准则要求J(n)对W(n)的导数等于零，即

$\frac{\∂J\left(n\right)}{\∂W}=-2\underset{i=0}{\overset{n}{\Σ}}{\mathrm{\λ}}^{n-i}\[d\left(n\right)-Y^T\left(n\right)W\left(n\right)\]Y\left(n\right)=0---\left(47\right)$

经整理后，得到标准方程

定义，

则标准方程，式(48)写成，R(n)W(n)＝q(n) (50)

该方程的解为，W(n)＝R^-1(n)q(n) (51)

这里R^-1(n)是方阵R(n)的逆矩阵。

根据定义式(49)，有

则,以上二式，写成以下递推形式：

R(n)＝λR(n-1)+Y(n)Y^T(n) (52-c)

q(n)＝λq(n-1)+d(n)Y(n) (52-d)

将矩阵求逆引理，[A+UV]^-1＝A^-1-A^-1U[I+VA^-1U]^-1VA^-1 (39)作用在(52-c)上：

定义向量K(n)：

将式(54)代入式(53)，则

将式(55)代入式(51)，有

因为K(n)的定义式，(54)式，则

将式(57)代入式(56)，得到自适应滤波系数向量的递推公式：

W(n)＝W(n-1)+K(n)e(n) (59)

在(59)式中，定义e(n)：

e(n)＝d(n)-W^T(n-1)Y(n)＝d(n)-Y^T(n)W(n-1) (60)

e(n)称为先验估计误差

根据前面的估计误差ε(n)的定义式，即式(45)：

ε(n)＝d(n)-x(n)＝d(n)-W^T(n)Y(n)＝d(n)-Y^T(n)W(n) (45)

对比式(60)与式(45)可知，e(n)与估计误差ε(n)不同，具体如下：

(60)式中，W^T(n-1)Y(n)是nT时刻自适应FIR滤器输出，利用的是(n-1)T时刻滤波系数向量W(n-1)，用W^T(n-1)Y(n)估计nT时刻的原始信号d(n)。其产生的估计误差e(n)，称为先验估计误差。

而在式(45)中，W^T(n)Y(n)是nT时刻自适应FIR滤器输出，但利用的是nT时刻滤波系数向量W(n)。用W^T(n)Y(n)估计nT时刻的原始信号d(n)，产生的估计误差e(n)，也被称为后验估计误差。

由于无法确切地知道nT时刻的原始信号d(n)，故用W^T(n-1)Y(n-1)来估计d(n)，用d(n)表示之：

d(n)＝W^T(n-1)Y(n-1) (61)

将式(60)中的d(n)用d(n)代替，则自适应滤波系数向量的递推公式：

综上所述，自适应FIR滤器算法如下：

一、初始化：

(1)定义n＝-1时，W^T(-1)＝[0,0,...,0]^T，R(-1)＝δ₁ eye，即R^-1(-1)＝δ₁ ^-1 eye，δ₁是很小的正数，eye为单位对角矩阵；

(2)定义n＝0时的抽样信号向量Y(0)：Y(0)＝[y(0),0..,0]，Y(0)为M维向量；

(3)定义n＝0时的原始信号估计值d(0)，d(0)＝δ₂，δ₂是很小的正数；

(4)选常数λ：0＜λ＜1。

二、由Y(0)、R^-1(-1)及λ，根据式(45),计算n＝0时的向量K(0)：

三、根据式(62),计算n＝0时的自适应FIR滤波系数向量W(0)：

e(0)＝d(0)-W^T(-1)Y(0)＝d(0)

W(0)＝W(-1)+K(0)e(0)＝K(0)e(0)

四、根据式(55)，计算n＝0时的矩阵R^-1(0)：

五、令n＝n+1，利用当前nT时刻的抽样信号y(n)，更新信号向量Y(n)：

Y(n)＝[y(n),y(n-1),...,y(n-M+1)]^T

六、由当前nT时刻信号向量Y(n)和(n-1)T时刻矩阵R^-1(n-1)，根据式(54)、(55)，计算当前nT时刻的向量K(n)、矩阵R^-1(n)：

七、由递推公式(62)式，计算当前nT时刻的滤波系数向量W(n)：

d(n)＝W^T(n-1)Y(n-1)

e(n)＝d(n)-W^T(n-1)Y(n)

W(n)＝W(n-1)+K(n)e(n)

并对抽样信号y(n)，按下式进行自适应FIR滤波：

重复以上五至七，如此循环，获得任意时刻的滤波后的信号序列x(n)。

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

Claims (4)

1.一种大气颗粒物PM2.5数密度检测方法，其特征在于，其包括以下步骤：

步骤1：含有多种粒径颗粒物的被测空气，经过不同孔径的过滤器分离后，形成一股含有PM2.5及更小颗粒物的气流，并进入气体测试室；

步骤2：在气体测试室中，用激光光源输出的红外光照射被测气流，一部分红外光与气流中颗粒物发生散射，形成散射光，另一部分红外光不与颗粒物发生散射，直接穿透被测气流，形成穿透光，采用光学部件，将散射光和穿透光汇聚在一起；

步骤3：汇聚在一起的散射光和穿透光，发生干涉效应，并输出干涉光；

步骤4：采用光电探测器测量以上干涉光的强度信号，其包含各散射光的自相关强度、穿透光的自相关强度、散射光与穿透光之间的互相关强度、各散射光之间的互相关强度四项内容，该干涉光的强度信号是随机信号，其概率分布密度函数与形成散射光的颗粒子总数目N相关；

步骤5：用A/D转换器对以上的干涉光总强度信号，以采样周期为T，进行离散化抽样，得到抽样信号序列y(n)，用y(n)表示在nT时刻的抽样信号，n为正整数；

步骤6：为了从被噪声污染的抽样信号中，恢复出原始的干涉光强信号，采用递推最小二乘RLS算法，设计一个M阶的自适应FIR滤波器，对抽样信号y(n)，按下式进行滤波：

$x\left(n\right)=\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}{w}_i\left(n\right)y(n-i)=W^T\left(n\right)Y\left(n\right)=Y^T\left(n\right)W\left(n\right)$

其中，w_i(n)为nT时刻的M阶FIR滤波器滤波系数，i＝0,1,...,M-1，W(n)为nT时刻滤波系数向量：W(n)＝[w₀(n),...,w_M-1(n)]^T，而x(n)为滤波后的信号序列，x(n)可视为在nT时刻的原始的干涉光强信号估计值，Y(n)为nT时刻所对应的信号向量，该向量由M个抽样信号组成，其如下：

Y(n)＝[y(n),y(n-1)...,y(n-M+1)]^T

W^T(n)、Y^T(n)分别为W(n)和Y(n)的转置：

W^T(n)＝[w₀(n),...,w_M-1(n)]，Y^T(n)＝[y(n),y(n-1),...,y(n-M+1)]；

步骤7：设滤波后的信号序列x(n)的序号指针为n₀，从n₀＝0开始，连续取最前面的二段长度为L的信号，L为正整数，用信号向量X₁、X₂表示：

X₁＝[x(1),...,x(L)]^T

X₂＝[x(L+1),...,x(2L)]^T

X₁、X₂分别对应在[0-LT]、[LT-2LT]时间内的干涉光强度信号， 分别是信号X₁、X₂的转置形式：

X₁ ^T＝[x(1),...,x(L)]

X₂ ^T＝[x(L+1),...,x(2L)]

按下式，计算X₁、X₂的平均值：

$\stackrel{\‾}{X_1}=\left(\underset{k=1}{\overset{L}{\Σ}}x\right(k\left)\right)/L$

$\stackrel{\‾}{X_2}=\left(\underset{k=1}{\overset{L}{\mathrm{\Σ}}}x\right(L+k\left)\right)/L$

L,M均为正整数，且L＞M；

步骤8：按下式，计算X₁、X₂的方差：

${{\mathrm{\σ}}_1}^2=\underset{k=1}{\overset{L}{\Σ}}{\[x\left(k\right)-\stackrel{\‾}{X_1}\]}^2$

${{\mathrm{\σ}}_2}^2=\underset{k=1}{\overset{L}{\Σ}}{\[x(L+k)-\stackrel{\‾}{X_2}\]}^2$

σ₁ ²和σ₂ ²分别对应在[0-LT]和[LT-2LT]时间内的干涉光强度的涨落:

步骤9：根据干涉光强度在[0-LT]、[LT-2LT]时间内的平均值和方差：

σ₁ ²、σ₂ ²，建立以下的参数C₁、C₂的矩阵方程：

其中，

a²为粒子散射光的总光强，N为形成散射光的粒子总数目，将C₁、C₂称为***参数，这里假设，a²和N在[0-LT]和[LT-2LT]时间内保持不变，因此，在此时间内，***参数C₁、C₂近似是一组常量；

步骤10：通过以上方程的最小二乘解，得到[0-LT]和[LT-2LT]时间内的***参数C₁和C₂的估计值：

其中，P₀＝(A₀ ^TA₀)^-1，A₀ ^T为矩阵A₀的转置，(A₀ ^TA₀)^-1为A₀ ^TA₀的逆矩阵，令将在[0-LT]和[LT-2LT]时间内的***参数C₁、C₂，分别用表示；

步骤11：调整信号序列x(n)的序号指针n₀，令n₀＝n₀+2L，使序号指针，指向下一段长度为L的滤波后的信号；

步骤12：从当前的序号指针n₀开始，取一段长度为L滤波后的信号：X₃＝[x(n₀+1),...,x(n₀+L)]^T，X₃对应于在[n₀T-(n₀+L)T]时间内的干涉光强度信号，按下式，计算X₃的平均值及方差σ₃ ²：

$\begin{array}{c}\stackrel{\‾}{X_3}=\left(\underset{k=1}{\overset{L}{\mathrm{\Σ}}}x\right(n_0+k\left)\right)/L\\ {{\mathrm{\σ}}_3}^2=\underset{k=1}{\overset{L}{\Σ}}{\[x\left(n_0+k\right)-\stackrel{\‾}{X_3}\]}^2\end{array},$

用新生成的平均值及方差构成信号更新向量d^T及b⁽¹⁾：

$d^T=\left[\begin{array}{cc}\frac{2\stackrel{\‾}{X_3}}{{{\mathrm{\σ}}_3}^2}& \frac{-2}{{{\mathrm{\σ}}_3}^2}\end{array}\right],b^{\left(1\right)}=1,$

建立在[n₀T-(n₀+L)T]时间内的***参数C₁、C₂的矩阵方程：

$d^T\left[\begin{array}{c}{C}_1\\ C_2\end{array}\right]=b^{\left(1\right)},b^{\left(1\right)}=1;$

步骤13：根据所构成的信号更新向量d^T及b⁽¹⁾，采用递推最小二乘RLS算法，按下式，完成在[n₀T-(n₀+L)T]时间内的对***参数C₁、C₂的进行估计，并分别用表示之：

$\begin{array}{c}{P}_1=P_0-\frac{P_0{\mathrm{dd}}^T{P}_0}{1+d^T{P}_{0}d}\\ \left[\begin{array}{c}{C}_1^{\left(1\right)}\\ C_2^{\left(1\right)}\end{array}\right]=\left[\begin{array}{c}{C}_1^{\left(0\right)}\\ C_2^{\left(0\right)}\end{array}\right]+\left(b^{\left(1\right)}-d^T\left[\begin{array}{c}{C}_1^{\left(0\right)}\\ C_2^{\left(0\right)}\end{array}\right]\right)P_{1}d\end{array},$

其中，向量d与d^T互为转置关系，与分别对应于本次迭代运算以前与之后的参数C₁、C₂；

步骤14：根据以上[n₀T-(n₀+L)T]时间内的***参数C₁和C₂，按下式，计算气体测试室在该时间内，散射光的总光强a²及形成散射光的粒子总数目N:

$a^2=C_1,N=\frac{{C_1}^2}{C_2};$

步骤15：根据已知的气体测试室的容积V和以上求得的粒子总数目N，能计算出在[n₀T-(n₀+L)T]时间内，被测气体的颗粒物数密度ρ，即单位体积的气体中PM2.5粒子数量：

$\mathrm{\ρ}=\frac{N}{V}$

步骤16：根据递推最小二乘RLS算法要求，按下式，重置迭代运算各单元的初值：

$\begin{array}{c}{P}_0=P_1\\ \left[\begin{array}{c}{C}_1^{\left(0\right)}\\ C_2^{\left(0\right)}\end{array}\right]=\left[\begin{array}{c}{C}_1^{\left(1\right)}\\ C_2^{\left(1\right)}\end{array}\right]\end{array};$

步骤17：调整信号序列X的序号指针n₀，令n₀＝n₀+L，使序号指针指向下一段长度为L的滤波后的信号；

步骤18：重复步骤12-17，如此循环，完成在起始时刻n₀＝0之后的，任意时间[n₀T-(n₀+L)T]内的***参数C₁、C₂的估计，进而，得到任意时间[n₀T-(n₀+L)T]内的散射光的总光强a²及形成散射光的粒子总数目N，实现被测气体在起始时刻n₀＝0之后的，任意时间[n₀T-(n₀+L)T]内的颗粒物数密度的检测，其中，n₀＝0,2L,3L,4L,5L,.......。

2.根据权利要求1所述的一种大气颗粒物PM2.5数密度检测方法，其特征在于：步骤1中，被测空气在负压发生单元的作用下，经过二级不同孔径的过滤器分离。

3.根据权利要求1所述的一种大气颗粒物PM2.5数密度检测方法，其特征在于：步骤2中，所述光学部件为凸透镜或凸透镜组。

4.根据权利要求1所述的一种大气颗粒物PM2.5数密度检测方法，其特征在于：步骤6中，自适应FIR滤波器最小二乘RLS算法，其特征在于，其包括以下步骤：

步骤1：初始化：

1)定义n＝-1时，W^T(-1)＝[0,0,...,0]^T，R(-1)＝δ₁eye，即R^-1(-1)＝δ₁ ^-1eye，δ₁是很小的正数，eye为单位对角矩阵；

2)定义n＝0时的抽样信号向量Y(0)：Y(0)＝[y(0),0..,0]，Y(0)为M维向量；

3)定义n＝0时的原始信号估计值d(0)，d(0)＝δ₂，δ₂是很小的正数；

4)选常数λ：0＜λ＜1；

步骤2：由Y(0)、R^-1(-1)及λ，按下式,计算n＝0时的向量K(0)：

$K\left(0\right)=\frac{R^{-1}(-1)Y\left(0\right)}{\mathrm{\λ}+Y^T\left(0\right)R^{-1}(-1)Y\left(0\right)}$

步骤3：按下式,计算n＝0时的自适应FIR滤波系数向量W(0)：

e(0)＝d(0)-W^T(-1)Y(0)＝d(0)

W(0)＝W(-1)+K(0)e(0)＝K(0)e(0)

步骤4：按下式，计算n＝0时的矩阵R^-1(0)：

$R^{-1}\left(0\right)=\frac{1}{\mathrm{\λ}}\[R^{-1}(-1)-K\left(0\right)Y^T{R}^{-1}(-1)\]$

步骤5：令n＝n+1，利用当前nT时刻的抽样信号y(n)，更新信号向量Y(n)：

Y(n)＝[y(n),y(n-1),...,y(n-M+1)]^T

步骤6：由当前nT时刻信号向量Y(n)和(n-1)T时刻矩阵R^-1(n-1)，按下式，计算当前nT时刻的向量K(n)、矩阵R^-1(n)：

$K\left(n\right)=\frac{R^{-1}(n-1)Y\left(n\right)}{\mathrm{\λ}+Y^T\left(n\right)R^{-1}(n-1)Y\left(n\right)}$

$R^{-1}\left(n\right)=\frac{1}{\mathrm{\λ}}\[R^{-1}(n-1)-K\left(n\right)Y^T{R}^{-1}(n-1)\]$

步骤7：由以下递推公式，计算当前nT时刻的滤波系数向量W(n)：

d(n)＝W^T(n-1)Y(n-1)

e(n)＝d(n)-W^T(n-1)Y(n)

W(n)＝W(n-1)+K(n)e(n)

并对抽样信号y(n)，按下式进行自适应FIR滤波：

$x\left(n\right)=\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}{w}_i\left(n\right)y(n-i)=W^T\left(n\right)Y\left(n\right)$

步骤8：重复以上步骤5-7，如此循环，获得任意时刻的滤波后的信号序列x(n)。