CN110161532B - 一种基于多波长激光雷达反演气溶胶微物理特性的方法 - Google Patents

Abstract

本发明公开一种基于多波长激光雷达反演气溶胶微物理特性的方法，包括步骤：输入光学特性参数，对光学特性参数进行归一化处理得到与气溶胶负载无关的光学特性参数量；建立气溶胶的微物理特性参量与多波长光学特性参数之间相关联的查找表；计算输入光学特性参数与查找表中各个个体在光学特性域的距离，并仅保留K个具有最小距离的个体作为可行解的范围；根据Monte Carlo随机采样理论，按照概率采样生成N_MC个不同的搜索执行顺序，并依次求解得到N_MC个备选解；平均N_MC个备选解，得到气溶胶微物理特性参数的最终反演结果。与现有技术相比，本发明无需额外的输入辅助判断信息，且执行效率高，可以满足激光雷达海量数据自动化处理的应用需求。

Description

一种基于多波长激光雷达反演气溶胶微物理特性的方法

技术领域

本发明涉及激光雷达大气探测反演技术，特别是一种基于多波长激光雷达反演气溶胶微物理特性的方法。

背景技术

气溶胶-云相互作用是目前气候变化预测中最大的不确定性来源之一，其一方面是由于气溶胶-云相互作用本身是比较复杂的系统，另一方面是由于缺少气溶胶-云的有效监测数据。激光雷达因其较高的时间和空间分辨率，被广泛应用于大气气溶胶的探测研究。高光谱分辨率激光雷达和Raman激光雷达的应用，在理论上解决了普通Mie散射激光雷达探测气溶胶光学特性参数(后向散射系数和消光系数)时“一个方程两个未知数”的数学欠定问题，提高了气溶胶光学特性的探测精度。但是，气溶胶的后向散射系数和消光系数不仅与气溶胶数浓度相关，还与其类型、粒径大小等密切相关，这就使得采用气溶胶的后向散射系数和消光系数直接表征气溶胶浓度，去研究气溶胶-云的相互作用时会存在较大的不确定性。如果能够探测得到气溶胶的复折射率、有效粒径、数浓度等微物理特性参数，则能更加有效的表征气溶胶的本质属性，能够提高对气溶胶-云相互作用研究的认知，进而降低气候变化预测的不确定性。

气溶胶的后向散射系数和消光系数是由其复折射率和粒径分布参数共同决定的，不同波段的后向散射系数和消光系数对同一气溶胶表现出不同的响应，也就是说多波长的后向散射系数和消光系数蕴含了气溶胶粒子的复折射率和粒径等微物理特性信息，这就为激光雷达探测气溶胶的微物理特性提供了必要的理论基础。目前已经开发了多种基于多波长激光雷达反演气溶胶微物理特性参数的算法，但是激光雷达探测反演气溶胶微物理特性参数依旧是比较有挑战性的工作。

基于正则化、混合正则化的气溶胶微物理特性反演方法是目前比较常用的方法，在已知气溶胶复折射率和大致粒径分布范围的条件下，该方法可以较好的反演气溶胶的粒子谱分布等微物理特性。但是实际情况下，气溶胶的复折射率和粒径分布的先验知识较难获得，这就需要选择不同的复折射率和粒径分布范围组合反复求解，并比较得到的反演结果，最终选择合理的结果。这一过程不仅耗时，而且不够自动化，可能会受到数据分析专家的主观经验影响。基于最小化偏差的反演方法不仅受到查找表本身的影响较大，而且对噪声较为敏感。NASA的科学家E.Chemyakin等提出了一种基于排列、平均的气溶胶微物理特性反演方法，该方法可以自动化执行，且降低了对数据查找表本身取值的依赖，对数据噪声具有一定的抑制能力，但代价是该方法执行效率比较低，难以处理具有较多维光学特性输入的情况。

发明内容

针对现有多波长激光雷达反演气溶胶其微物理特性算法存在的不足，本发明的目的在于提供一种基于多波长激光雷达反演气溶胶微物理特性参数的方法。该方法仅需要输入待反演个例的光学特性参数，无需额外的输入辅助信息，即可反演得到气溶胶的复折射率、有效粒径、数浓度、表面积浓度及体积浓度等微物理特性反演结果。该方法操作简单，运行效率高，可自动化执行，能够适用于海量激光雷达数据处理的要求。

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

一种基于多波长激光雷达反演气溶胶微物理特性的方法，包括步骤：

S1：输入多波长激光雷达探测气溶胶后得到的光学特性参数，包括后向散射系数和消光系数，对所述光学特性参数进行归一化处理，得到与气溶胶负载无关的N个归一化光学特性参数；

S2：建立一个气溶胶微物理特性参数与光学特性参数关联的查找表，存储有对应数浓度为1cm^-3时的气溶胶个体的多波长后向散射系数、消光系数、归一化后的光学特性参数及对应的气溶胶复折射率、有效粒径、表面积浓度、体积浓度等微物理特性参数；

S3：计算所述光学特性参数与查找表中各气溶胶个体#j在光学特性域的距离d^#j，并仅保留K个具有最小距离的气溶胶个体作为可行解的范围；

S4：根据Monte Carlo随机采样原理，按概率采样生成N_MC个搜索执行顺序，在每个搜索执行顺序中，依次计算待反演个体与查找表中可行解的范围内气溶胶个体的距离，排除筛选后，取对于最后一个光学特性量有最小距离的气溶胶个体作为单次求解的备选解，共得到N_MC个备选解；

S5：计算各备选解对应的物理特性参数，其中数浓度、表面积浓度和体积浓度的公式如下：

式中，

和

分别为第#p个备选解在查找表中对应数浓度为1cm^-3时的表面积浓度和体积浓度；

和

分别为第#p个备选解的数浓度、表面积浓度和体积浓度，λ_i对应不同激光雷达探测波长，n和m分别为后向散射系数和消光系数的探测波长数量，β^*和α^*分别为输入个例对应的后向散射系数和消光系数，β^#p和α^#p分别为第#p个备选解在查找表中对应的后向散射系数和消光系数。

平均N_MC个备选解对应的物理特性参数，最终反演得到气溶胶微物理特性参数，包括气溶胶复折射率、有效粒径、数浓度、表面积浓度、体积浓度等。

在步骤S1中，按下述方式对光学特性参数进行归一化处理得：

S_i,j＝α(λ_j)/β(λ_i)

式中，α为消光系数，β为后向散射系数，i＝1,...,n,j＝1,...,m,λ为激光雷达的探测波长。

所述步骤S3中，距离d^#j的计算公式如下：

d^#j＝∑_G∈Ω|(G^*-G^#j)/G^*|

式中，Ω为步骤S1中得到的N个光学特性参数集合，G*表示输入光学特性参数对应的归一化后的光学特性参数，G^#j表示查找表中第#j个个体对应的归一化光学特性参数。

所述的查找表中存储有气溶胶个体的多波长后向散射系数、消光系数及归一化后的光学特性参数。

所述步骤S4中，Monte Carlo随机采样的次数N_MC小于全排列

且生成1～N的随机排列的采样方法可以不同。

所述步骤S4中，所述的排除筛选为：按照距离

由小到大排列，排除距离较大的(1-ω)部分[ω∈(0,1)]，依次缩小备选解的可行解空间，以便最终确定该次求解的备选解。系数ω的取值与输入光学特性参数β和α的个数、查找表的构成有关，可以根据仿真数据分析，选取使得气溶胶微物理特性反演误差最小的ω取值，再将其应用到实际的激光雷达数据处理。

所述步骤S5中，每个备选解#p，直接在所述的查找表获得对应的气溶胶复折射率

有效粒径

求平均后得到最终的气溶胶复折射率和有效粒径。

相较于现有技术，本发明的反演方法，实现起来十分简单，无需数据处理人员的经验判断作为额外输入，可以实现算法的自动化运行。另外，该方法通过K-近邻约束可以减小求解结果的方差，提高求解的准确度，通过MC采样可以提升算法的执行效率，这样就可以实现海量激光雷达数据的实时自动化处理分析的需求。

附图说明

图1是本发明实施例提供的一种基于多波长激光雷达反演气溶胶微物理特性参数方法流程图。

图2是本发明实施例提供的基于多波长激光雷达反演气溶胶微物理特性参数方法的单个备选解求解流程图。

图3是本发明实施例提供的一种基于多波长激光雷达反演气溶胶微物理特性参数方法气溶胶有效粒径反演统计结果的示意图。

具体实施方式

下面结合实施例和附图来详细说明本发明，但本发明并不仅限于此。

图1是本发明提出的气溶胶微物理特性反演算法的总体流程图，包括：

步骤S1：假定输入的光学特性参数(后向散射系数β和消光系数α)来自三波长(采用Nd:YAG激光器作为激光光源，探测波长为355nm，532nm和1064nm)高光谱分辨率激光雷达/Raman激光雷达，即n＝m＝3。那么可以得到N＝15个归一化光学特性参数，分别为：

B₁＝β(λ₁)/β₃,B₂＝β(λ₂)/β₃,B₃＝β(λ₃)/β₃，

A₁＝α(λ₁)/α₃,A₂＝α(λ₂)/α₃,A₃＝α(λ₃)/α₃，

S_1,1＝α(λ₁)/β(λ₁),S_1,2＝α(λ₂)/β(λ₁),S_1,3＝α(λ₃)/β(λ₁)，

S_2,1＝α(λ₁)/β(λ₂),S_2,2＝α(λ₂)/β(λ₂),S_2,3＝α(λ₃)/β(λ₂)，

S_3,1＝α(λ₁)/β(λ₃),S_3,2＝α(λ₂)/β(λ₃),S_3,3＝α(λ₃)/β(λ₃)，

式中，λ₁＝355nm，λ₂＝532nm，λ₁₀₆₄＝1064nm，

那么，也就可以得到归一化光学特性集合为Ω＝{B₁,B₂,B₃,A₁,A₂,A₃,S_1,1,S_1,2,S_1,3,S_2,1,S_2,2,S_2,3,S_3,1,S_3,2,S_3,3}。

步骤S2：建立气溶胶微物理特性参数与光学特性参数关联的查找表。查找表中存储了气溶胶个体的多波长后向散射系数、消光系数及归一化光学特性参数，其微物理特性主要由粒子的复折射率和粒子谱分布确定，在本实施例中，查找表的构建考虑气溶胶粒子的粒径分布为单峰对数正太分布，即

式中，r为气溶胶的粒径，r_med为粒径分布模态的均值，σ为粒径分布模态的模式宽度。气溶胶的数浓度均固定为1cm^-3。本实施中，查找表中气溶胶复折射率和粒径分布参数的范围为复折射率实部m_r的取值范围为1.29～1.71，取值间隔为0.02；复折射率虚部m_i的取值范围为0和0.25×10^-3～50×10^-3，取值间隔为10^-3；粒径参数r_med取值范围为15～305nm，取值间隔为10nm；粒径分布宽度参数σ的取值范围为1.475～2.525，取值间隔为0.05。查找表中共有755040个气溶胶个例，其对应的有效粒径、表面积浓度及体积浓度等可以按照下式计算，

有效粒径：r_eff＝∫r³n(r)dr/∫r²n(r)dr＝r_medexp(-2.5ln²σ)，

表面积浓度：s_t＝4π∫r²n(r)dr＝4πr_medexp(-2ln²σ)，

体积浓度：v_t＝4π/3∫r³n(r)dr＝4π/3r_medexp(-4.5ln²σ)。

而其对应的后向散射系数和消光系数可以按照T-矩阵或者Mie散射理论计算得到，归一化光学特性参数则由步骤S1所述方式计算。在本实施例中，假定粒子为球形等效粒子，采用Mie散射计算其光学特性

式中，K_β和K_α分别为后向散射系数和消光系数的核函数，r为气溶胶粒径，λ为激光波长，m_r+im_i为气溶胶复折射率，m_r和m_i分别为复折射率的实部和虚部。K_β(r,λ,m_r+im_i)＝πr²Q_bsc(r,λ,m_r+im_i)，K_α(r,λ,m_r+im_i)＝πr²Q_ext(r,λ,m_r+im_i)，Q_bsc和Q_ext分别为Mie散射计算得到的后向散射系数和消光效率因子。

步骤S3：计算输入个例与查找表中每一个个体的在归一化光学特性域Ω的距离

并按照从小到的顺序排序，选取查找表中距离最小的K＝50000个个体作为可行解的范围。

步骤S4：图2给出了步骤S4的具体查找备选解的流程，按照光学特性构建的搜索路径共有

种等可能的搜索路径，这里按照Monte Carlo采样的思路生成一次随机的搜索路径，在本实施例中介绍一种生成随机搜索路径的方式，但是MC采样生成随机路径的方法不局限于此方法。首先生成15个[0,1]均匀分布的随机数，

1	2	3	4	5	6	7	8
								0.1270	0.9134	0.6324	0.0975	0.2785	0.5469	0.9575	0.9649
9	10	11	12	13	14	15
								0.1576	0.9706	0.9572	0.4854	0.8003	0.1419	0.4218

将这15个随机数按照从小到大排序：

4	1	14	9	5	15	12	6
								0.0975	0.1270	0.1419	0.1576	0.2785	0.4218	0.4854	0.5469
3	13	2	11	7	8	10
								0.6324	0.8003	0.9134	0.9572	0.9575	0.9649	0.9706

则可以得到采样的随机序列为4→1→14→…→8→10，也即该次备选解的搜索次序为G₄→G₁→G₁₄→…→G₈→G₁₀，也即A₁→B₁→S_3,2→…→S_1,2→S_2,1。

从n＝1～14,也即依次按照光学变量A₁→B₁→S_3，2→…→S_1，2计算

并排序，每次仅保留距离

最小的ω＝55％作为可能备选解的搜索范围，即可将备选解范围缩小到K*ω¹⁴＝50000*0.55¹⁴≈12个个体，对于最后一维光学变量(n＝15)S_2,1计算距离参量

并选择

最小的个体作为该次求得的备选解。重复该过程，得到N_MC＝500个备选解。值得注意的是，每一个备选解的求解过程都是相互独立的，因此很容易进行并行化算法设计。

步骤S5：计算每一个备选解#p对应的

数浓度：

表面积浓度：

体积浓度：

备选解对应的气溶胶复折射率

及有效粒径

可以从查找表中直接得到。进而平均500个备选解，得到最终的气溶胶复折射率、有效粒径、数浓度、表面积浓度、体积浓度等微物理特性参数的反演结果。图3给出了对2880个不同的输入光学特性参数反演气溶胶有效粒径的相对误差的统计结果，直方图表示统计的个例数，对应右边的y坐标，曲线表示反演相对误差的积分概率密度，对应左边的y坐标。

以上所述仅为本发明的较佳实施举例，并不用于限制本发明，凡在本发明精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于多波长激光雷达反演气溶胶微物理特性的方法，其特征在于，包括步骤：

S1：输入多波长激光雷达探测气溶胶后得到的光学特性参数，包括后向散射系数和消光系数，对所述光学特性参数进行归一化处理，得到与气溶胶负载无关的N个光学特性参数；

S2：建立气溶胶微物理特性参数与光学特性参数关联的查找表；

S3：计算所述光学特性参数与查找表中各气溶胶个体#j在光学特性域的距离d^#j，并仅保留K个具有最小距离的气溶胶个体作为可行解的范围；

S4：根据Monte Carlo随机采样原理，按概率采样生成N_MC个搜索执行顺序，在每个搜索执行顺序中，依次计算待反演个体与查找表中可行解的范围内气溶胶个体的距离，排除筛选后，取对于最后一个光学特性量有最小距离的气溶胶个体作为单次求解的备选解，共得到N_MC个备选解；

S5：计算各备选解对应的物理特性参数，其中数浓度、表面积浓度和体积浓度的公式如下：

式中，

和

分别为第#p个备选解在查找表中对应数浓度为1cm^-3时的表面积浓度和体积浓度；

和

分别为第#p个备选解的数浓度、表面积浓度和体积浓度，λ_i对应不同激光雷达探测波长，n和m分别为后向散射系数和消光系数的探测波长数量，β^*和α^*分别为输入个例对应的后向散射系数和消光系数，β^#p和α^#p分别为第#p个备选解在查找表中对应的后向散射系数和消光系数；

平均N_MC个备选解对应的物理特性参数，最终反演得到气溶胶微物理特性参数。

2.如权利要求1所述的基于多波长激光雷达反演气溶胶微物理特性的方法，其特征在于，在步骤S1中，按下述方式对光学特性参数进行归一化处理得：

S_i,j＝α(λ_j)/β(λ_i)

式中，i＝1,...,n,j＝1,...,m，β为后向散射系数和α为消光系数，n和m分别为后向散射系数和消光系数的探测波长数量。

3.如权利要求1所述的基于多波长激光雷达反演气溶胶微物理特性的方法，其特征在于，所述步骤S2中，查找表中存储有对应数浓度为1cm^-3时的气溶胶个体的多波长后向散射系数、消光系数、归一化后的光学特性参数及对应的气溶胶复折射率、有效粒径、表面积浓度、体积浓度；

后向散射系数：β(λ)＝∫K_β(r,λ,m_r+im_i)n(r)dr

消光系数：α(λ)＝∫K_α(r,λ,m_r+im_i)n(r)dr

有效粒径：r_eff＝∫r³n(r)dr/∫r²n(r)dr

表面积浓度：s_t＝4π∫r²n(r)dr

体积浓度：v_t＝4π/3∫r³n(r)dr

式中，λ为对应的激光探测波长，r为气溶胶的粒径，n(r)为气溶胶粒子的粒径分布，K_β和K_α分别为后向散射系数和消光系数的核函数，m_r+im_i为气溶胶复折射率，i为虚数单位，m_r和m_i分别为复折射率的实部和虚部。

4.如权利要求1所述的基于多波长激光雷达反演气溶胶微物理特性的方法，其特征在于，所述步骤S3中，距离d^#j的计算公式如下：

d^#j＝∑_G∈Ω|(G^*-G^#j)/G^*|

式中，Ω为步骤S1中得到的N个光学特性参数集合，G^*表示输入光学特性参数对应的归一化后的光学特性参数，G^#j表示步骤S2中所构建查找表中第#j个个体的归一化后的光学特性参数。

5.如权利要求1所述的基于多波长激光雷达反演气溶胶微物理特性的方法，其特征在于，所述步骤S4中，Monte Carlo随机采样的次数N_MC小于全排列

6.如权利要求1所述的基于多波长激光雷达反演气溶胶微物理特性的方法，其特征在于，所述步骤S4中，所述的排除筛选为：按照距离

由小到大排列，排除距离较大的(1-ω)部分，ω∈(0,1)。

7.如权利要求1所述的基于多波长激光雷达反演气溶胶微物理特性的方法，其特征在于，所述步骤S5中，根据每个备选解#p，直接在所述的查找表获得对应的气溶胶复折射率

有效粒径

求平均后得到最终的气溶胶复折射率和有效粒径。

Images