[发明专利]一种大气颗粒物PM2.5数密度检测的方法

权利要求书

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(n)=Σi=0M-1wi(n)y(n-i)=WT(n)Y(n)=YT(n)W(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₂的平均值：

X1‾=(Σk=1Lx(k))/L]]>

X2‾=(Σk=1Lx(L+k))/L]]>

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

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

σ12=Σk=1L[x(k)-X1‾]2]]>

σ22=Σk=1L[x(L+k)-X2‾]2]]>

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

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

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

A0C1C2=b(0),]]>其中，A0=2X1‾σ12-2σ122X2‾σ22-2σ22,]]>b(0)=11,]]>C₁＝a²,C2=a4N]]>

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

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

C1C2=P0A0Tb(0),]]>其中，P₀＝(A₀^TA₀)^-1，A₀^T为矩阵A₀的转置，(A₀^TA₀)^-1为A₀^TA₀的逆矩阵，令C1(0)C2(0)=C1C2,]]>将在[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₃的平均值及方差σ₃²：

X3‾=(Σk=1Lx(n0+k))/L]]>

σ32=Σk=1L[x(n0+k)-X3‾]2,]]>

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

dT=2X3‾σ32-2σ32,]]>b⁽¹⁾＝1，

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

dTC1C2=b(1),]]>b⁽¹⁾＝1；

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

P1=P0-P0ddTP01+dTP0d]]>

C1(1)C2(1)=C1(0)C2(0)+(b(1)-dTC1(0)C2(0))P1d,]]>

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

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

a²＝C₁,N=C12C2;]]>

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

ρ=NV]]>

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

P₀＝P₁

C1(0)C2(0)=C1(1)C2(1);]]>

步骤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,.......。