应用于扬尘颗粒物在线监测的mie散射快速计算方法与流程

本发明涉及扬尘颗粒物浓度测量技术领域，具体涉及一种应用于扬尘颗粒物在线监测的mie散射快速计算方法。

背景技术：

大气中颗粒物来源复杂，研究结果表明，扬尘是造成大气颗粒物污染的重要因素。扬尘粒径范围大、浓度变化快速，因此监测扬尘变化需要适用范围广、测量速度快、可在线实时测量的方法。

光散射法测量颗粒物浓度具有测量速度快、装置简单使用方便、能够做到实时测量等优点。光散射的原理是指光照射到介质中的颗粒物时，与介质相互作用会偏离传播方向，沿各个方向散射导致光的衰减。光沿各个方向的散射情况取决于光的波长、颗粒物的粒径以及颗粒物的折射率等参数。Mie散射理论是德国科学家Gustav Mie用经典波动光学理论的麦克斯韦方程推导得到的任意成分均匀球型粒子的散射光强角分布的严格数学解。通过100多年的研究与发展，尤其是随着计算机科学的发展，Mie散射理论已经趋于成熟并应用于生产实践。

目前实际应用于大气颗粒物浓度测量的光散射方法主要是角散射法、衍射散射法和全散射法(又称透射法、消光法或浊度法)。常用的数据处理方法是函数限制法和独立模式解法。本发明用于函数限制解法，即假设颗粒的尺寸分布符合某个函数关系，通过优化算法不断反演还原颗粒的尺寸分布，进而获得颗粒物浓度信息。但目前数据处理过程中因为初值选取不当和反复调用Mie散射函数导致速度降低，一般每次数据处理时间达到90s以上，极大地影响了 仪器的响应速度，不能满足扬尘颗粒物浓度的跟踪测量需求。

技术实现要素：

为解决上述问题，本发明提供了一种应用于扬尘颗粒物在线监测的mie散射快速计算方法。

为实现上述目的，本发明采取的技术方案为：

应用于扬尘颗粒物在线监测的mie散射快速计算方法，包括如下步骤：

S1、计算粒子的散射系数，形成散射系数表格：

根据扬尘的折射率m与入射光波长λ_i(1≤i≤N)基于Mie散射理论计算出扬尘空间存在的各种粒径D_j(1≤j≤M)下粒子的理论散射系数k_ext，形成N*M粒子散射系数表格：

$k_{ext}=\frac{2}{{\mathrm{\α}}^2}\underset{n=1}{\overset{\mathrm{\∞}}{\Σ}}(2n+1)\mathrm{Re}(a_n+b_n)$

式中，α＝πD_j/λ_i为尺寸参数，a_n、b_n为m、D_j、λ_i的函数，由Mie散射理论获得，这里不再赘述；

S2、建立粒子粒径分布数学模型，选定优化算法中模型参量的合理初始值：

S21、根据扬尘样本数据或者理论数据拟定多峰分布数学模型f(a₁,a₂,...,a_N,x)，其中a₁,a₂,...,a_N为数学模型的模型参量，参量的个数与入射光波长数量相等，可根据模型参量个数调整入射光波长的个数；x为粒子粒径；

S22、判断是否为首次测量，如果是，则在扬尘浓度分布合理范围内使用网格法寻找目标函数参量的合适初始值；如果该次测量不是首次测量，则将上一次数据处理获得的最优模型参量作为下一次数据处理的模型参量的初始值，

即${(a_1,a_2,...,a_N)}_n^0={(a_1,a_2,...,a_N)}_{n-1}^{best};$

S3、反演颗粒物粒径分布，计算颗粒物浓度：

S31、基于函数限制解法，利用粒子浓度分布模型结合粒子理论散射系数计算出各个波长下的一组浊度值，与实测浊度值对比计算残差平方和，通过 L-M优化算法的反复迭代，找到与实际粒子分布最接近的粒子模型最优参量；

S32、将最优参量代入到粒子浓度分布模型，通过积分得到固定粒径范围内粒子的浓度信息。

优选的，所述步骤S1所得的表格仅在所测扬尘颗粒物折射率发生较大变化或者使用的单色光波长变化时做更改，数据处理时仅调用表格中的数值。

优选的，所述步骤S22的具体步骤为：根据扬尘分布模型参量的合理范围，在每个参量的合理范围均等取出L个数值形成N维网格，共有L^N个网格节点；分别以各个节点处的数学模型参数值为初始值，计算各光路浊度值与实际值的最小残差平方和ε，取其中最小的ε所对应的模型参量初始值作为本次浓度计算的初始值。

优选的，所述步骤S32的具体步骤为：

1)将粒子粒径分布模型参量与浊度τ的关系设为函数表达式

τ_i＝g(a₁,a₂,...,a_N,λ_i)；

式中，τ_i为第i个光波长下的计算出的浊度值，a₁,a₂,...,a_N为模型参量，λ_i为第i个光路的波长；

设置L-M算法的迭代计算的最大次数K、和迭代结果允许残差平方和最大值ε_max；

2)计算以模型参量值为基础的各光路浊度值与实际值的残差平方和ε；

$\mathrm{\ϵ}=\underset{i=1}{\overset{N}{\Σ}}|{\mathrm{\τ}}_i-{\mathrm{A\τ}}_i{|}^2$

式中，ε为残差平方和，τ_i为计算第i个光波长下的浊度值，Aτ_i为实测第i个光波长下的浊度值；

3)判断当前的残差平方和ε是否小于允许残差平方和最大值ε_max,如果是则跳转执行步骤4)；否则判断此时迭代次数是否小于迭代计算的最大次数N，如果是则这用L-M算法计算合理的迭代步长，更新模型参量a₁,a₂,...,a_N的预测值，跳转执行步骤2)；否则跳转执行步骤4)；L-M优化算法的具体实现过程 在此不再赘述；

4)将迭代计算结束后的最后一组模型参量a₁,a₂,...,a_N作为最优模型参量；将最优模型参量代入到粒子浓度分布模型，通过积分得到固定粒径范围内粒子的浓度信息。

本发明具有以下有益效果：

本发明提前计算粒子的理论散射系数，采用查表的方法用于模型参量优化计算，避免多次循环运算，极大缩减计算时间；算法中利用网格法寻求进行第一次优化算法的合理模型参量初始值，以此准确获得目标函数的全局最小值；非第一次测量的模型参量初始值设定为上一次计算出的模型最优参量，减少迭代次数，增加运算速度；本发明可应用于扬尘颗粒物浓度测量，能够适应扬尘颗粒物粒径范围大、浓度变化快速的特性，用于测量系统做到对扬尘颗粒物浓度的实时测量。

附图说明

图1为本发明实施例应用于扬尘颗粒物在线监测的mie散射快速计算方法中步骤S1的流程图；

图2为本发明实施例应用于扬尘颗粒物在线监测的mie散射快速计算方法中步骤S2的流程图；

图3为本发明实施例应用于扬尘颗粒物在线监测的mie散射快速计算方法中步骤S3的流程图。

具体实施方式

为了使本发明的目的及优点更加清楚明白，以下结合实施例对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1-3所示，本发明实施例提供了一种应用于扬尘颗粒物在线监测的 mie散射快速计算方法，包括如下步骤：

S1、计算粒子的散射系数，形成散射系数表格：

根据扬尘的折射率m与入射光波长λ_i(1≤i≤N)基于Mie散射理论计算出扬尘空间存在的各种粒径D_j(1≤j≤M)下粒子的理论散射系数k_ext，形成N*M粒子散射系数表格：

$k_{ext}=\frac{2}{{\mathrm{\α}}^2}\underset{n=1}{\overset{\mathrm{\∞}}{\Σ}}(2n+1)\mathrm{Re}(a_n+b_n)---\left(1\right)$

式中，α＝πD_j/λ_i为尺寸参数，a_n、b_n是m、D_j、λ_i的函数，由Mie散射理论获得，这里不再赘述；该表格仅在所测扬尘颗粒物折射率发生较大变化或者使用的单色光波长变化时做更改，数据处理时仅调用表格中的数值。

S2、建立粒子粒径分布数学模型，选定优化算法中模型参量的合理初始值：

S21、根据扬尘样本数据或者理论数据拟定多峰分布数学模型f(a₁,a₂,...,a_N,x)，其中a₁,a₂,...,a_N为数学模型的模型参量，参量的个数与入射光波长数量相等，可根据模型参量个数调整入射光波长的个数；x为粒子粒径。

S22、判断是否为首次测量，如果是，则在扬尘浓度分布合理范围内使用网格法寻找目标函数参量的合适初始值。根据扬尘分布模型参量的合理范围，在每个参量的合理范围均等取出L个数值形成N维网格，共有L^N个网格节点。分别以各个节点处的数学模型参数值为初始值，计算各光路浊度值与实际值的最小残差平方和ε(具体计算方法见第三部分)，取其中最小的ε所对应的模型参量初始值作为本次浓度计算的初始值。如果该次测量不是首次测量，则将上一次数据处理获得的最优模型参量作为下一次数据处理的模型参量的初始值。

即${(a_1,a_2,...,a_N)}_n^0={(a_1,a_2,...,a_N)}_{n-1}^{best}---\left(2\right)$

S3、反演颗粒物粒径分布，计算颗粒物浓度：

S31、基于函数限制解法，利用粒子浓度分布模型结合粒子理论散射系数 计算出各个波长下的一组浊度值，与实测浊度值对比计算残差平方和，通过L-M优化算法的反复迭代，找到与实际粒子分布最接近的粒子模型最优参量；

S32、将最优参量代入到粒子浓度分布模型，通过积分得到固定粒径范围内粒子的浓度信息。具体实现方式如下：

1)将粒子粒径分布模型参量与浊度τ的关系设为函数表达式(3)，设置L-M算法的迭代计算的最大次数K、和迭代结果允许残差平方和最大值ε_max；

τ_i＝g(a₁,a₂,...,a_N,λ_i) (3)

式中，τ_i为第i个光波长下的计算出的浊度值，a₁,a₂,...,a_N为模型参量，λ_i为第i个光路的波长；

2)计算以模型参量值为基础的各光路浊度值与实际值的残差平方和ε；

$\mathrm{\ϵ}=\underset{i=1}{\overset{N}{\Σ}}|{\mathrm{\τ}}_i-{\mathrm{A\τ}}_i{|}^2---\left(4\right)$

式中，ε为残差平方和，τ_i为计算第i个光波长下的浊度值，Aτ_i为实测第i个光波长下的浊度值；

3)判断当前的残差平方和ε是否小于允许残差平方和最大值ε_max,如果是则跳转执行步骤4)；否则判断此时迭代次数是否小于迭代计算的最大次数N，如果是则这用L-M算法计算合理的迭代步长，更新模型参量a₁,a₂,...,a_N的预测值，跳转执行步骤2)；否则跳转执行步骤4)。L-M优化算法的具体实现过程在此不再赘述；

4)将迭代计算结束后的最后一组模型参量a₁,a₂,...,a_N作为最优模型参量；将最优模型参量代入到粒子浓度分布模型，通过积分得到固定粒径范围内粒子的浓度信息。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。