CN113009437B - 一种根据实时气溶胶分布求解激光雷达动态雷达比的方法 - Google Patents

Abstract

本发明公开了一种根据实时气溶胶分布求解激光雷达动态雷达比的方法，适用于晴空、少云的条件下。本发明包括以下步骤：1)气溶胶消光系数传感器测量得到地面实测气溶胶消光系数α_a(r₀)；2)计算地面处雷达信号P(r₀)；3)给定适合本地的预估雷达比；4)拟合地面及混合均匀处的消光系数；5)反代入Fernald后向积分方程求解动态雷达；6)将反代算法进行迭代；7)反演获得气溶胶消光系数的分布廓线。本发明在Fernald法的基础上，提出一种根据实时气溶胶分布求解动态雷达比的方法，能够实现雷达比本地化的目的，用于提高激光雷达反演算法的精度，从而可以获得更精确的气溶胶消光系数的分布廓线。

Description

一种根据实时气溶胶分布求解激光雷达动态雷达比的方法

技术领域

本发明涉及一种求解雷达比的方法，具体涉及一种根据实时气溶胶分布求解激光雷达动态雷达比的方法。

背景技术

激光雷达的原始方程为

其中，r表示高度；P(r)表示激光雷达归一化后向散射信号；C表示雷达常数，E表示脉冲能量，在雷达稳定工作的情况下，这两个值被视为常数。值得注意的是，(1)式中包含两个未知变量，显然无法直接求解。Fernald在求解(1)时将大气看成大气分子与气溶胶两部分组成，即认为大气消光系数α是气溶胶消光系数α_a与大气分子的消光系数α_m之和(标注a、m分别代表气溶胶和大气分子，后同)；后向散射系数β也看成由气溶胶后向散射系数β_a和大气分子后向散射系数β_m两部分组成。则(1)式化为：

通常假设消光系数与后向散射系数比值(雷达比)为常数，对于大气分子其消光α_m(r)和后向散射系数β_m(r)可以根据标准大气廓线计算得理论值，故有：对于气溶胶，有S_a＝α_a/β_a，其数值范围一般在10到100sr之间。据此Fernald给出了激光雷达方程的解析解：

其中，参考位置r_c取6km，并且有α_a(r_c)+α_m(r_c)＝1.05α_m(r_c)。由后向积分公式(3)可知，根据某一参考位置的气溶胶消光系数α_a(r_c)和该时刻的气溶胶消光-后向散射比S_a，即可反演出气溶胶消光系数的垂直廓线。另外推导过程中基本没有省略与近似处理，雷达比的选取是最直接的影响反演精度的因子。

雷达比(其中β为大气后向散射系数，α为大气消光系数)是时间和空间的复杂函数，它不仅与气溶胶粒子的尺度谱分布有关，还和气溶胶的组成紧密相关，其数值范围一般在10到100sr之间。雷达比的精确性直接决定了气溶胶消光廓线的反演精度，在气溶胶分布变化较大的地区不精确的雷达比对反演结果的影响更大。在实际应用中，由于气溶胶或云雾粒子的光学特性非常复杂，准确的激光雷达比取决于它们的光学性质要实时获取气溶胶的属性参数进而计算动态雷达比并不容易，因此在实际应用中往往将雷达比取固定值参与雷达反演。但是，固定值的使用必然会给雷达反演带来误差。

目前雷达比的计算方法，主要有两大类：独立测量法和联合反演法。独立测量分为多角度测量法、高光谱分辨率雷达法和拉曼雷达法，目前广泛使用拉曼雷达法。联合反演主要有利用卫星测得的光学厚度限制激光雷达测得的气溶胶消光系数的垂直积分法。尽管拉曼雷达和高光谱分辨率激光雷达能够实现不同高度气溶胶激光雷达比的有效探测，但二者分光系统复杂、回波信号微弱、探测距离短、信噪比低，这限制了二者的实际应用，而通过卫星数据限制激光雷达的联合法，因对单站激光雷达来说卫星数据空间分辨率过低，故两者数据不能完全对应从而会造成较大误差。

发明内容

为了解决现有反演方法存在的缺陷，进一步提高雷达比的反演精度，本发明在Fernald法的基础上，提出一种根据实时气溶胶分布求解激光雷达动态雷达比的方法，能够实现雷达比本地化，动态化的目的，用于提高激光雷达反演算法的精度，从而可以获得更精确的气溶胶消光系数的分布廓线。本发明提出的求解雷达比的方法适用于晴空，少云的条件下。本发明的方法包括以下步骤：

(1)气溶胶消光系数传感器测量得到地面实测气溶胶消光系数α_a(r₀)

(2)计算地面处雷达信号P(r₀)；

(3)给定适合本地的预估雷达比；

(4)拟合地面及混合均匀处的消光系数；

(5)反代入Fernald后向积分方程求解动态雷达比；

(6)将反代算法进行迭代；

(7)利用动态雷达比反演获得气溶胶消光系数的分布廓线。

进一步地，所述激光雷达在地面r₀处的信号P(r₀)，由激光雷达在最低有效位置r_m处的信号P(r_m)和地面实测气溶胶消光系数α_a(r₀)求得：

其中α_m为大气分子的消光系数。

进一步地，调整Fernald后向积分方程中Sa的值使得反演得到的气溶胶消光系数α_a与地面实测气溶胶消光系数α_a(r₀)大小相近，此时Sa的值为适合本地的预估雷达比；

由预估雷达比及步骤(2)得到的地面信号P(r₀)根据Fernald后向积分方程拟合地面消光系数σ_a(r′₀)，如下式所示：

式中：r表示高度，P(r)表示激光雷达归一化后向散射信号，C表示雷达常数，E表示脉冲能量，在雷达稳定工作的情况下，这两个值被视为常数。α_a为气溶胶消光系数；α_m为大气分子的消光系数，Sa为雷达比；Sm为大气分子的雷达比；

矫正系数H为：

进一步地，激光雷达原始方程为：

对上式两边求导，得：

利用最小二乘法拟合ln(P(r))与高度r，并计算求得这段高度上所有点偏离拟合直线的残差平方和，残差平方和最小处为所求混合均匀的高度。

根据步骤(3)给出的预估雷达比，反演得气溶胶消光廓线，从而根据所述气溶胶消光廓线的对应高度得到混合均匀处的消光系数σ_a(r′_u)；σ_a(r′_u)乘上相应的矫正系数H，即为所求混合均匀处的消光系数α_a(r_u)。

进一步地，所述Fernald给出的激光雷达方程的解析解中共有3个未知量分别为：参考点处的消光系数、待求点处的消光系数和雷达比S_a；

使用所述步骤(4)中得到的混合均匀处的消光系数α_a(r_u)作为参考点处消光；地面实测处消光系数α_a(r₀)作为待求点处消光系数。

进一步地，使用步骤(2)中计算得到的地面处雷达信号P(r₀)作为待求点的回波信号，与所述的混合均匀处的消光系数α_a(r_c)和地面处实测消光系数α_a(r₀)一同带入Fernald后向积分方程，用于反推出雷达比。

进一步地，所述步骤(6)中将所述的算法进行迭代，选择进行两次迭代。

进一步地，在反演气溶胶消光系数分布廓线时，取参考位置r_c为6km，并且在此高度有：

α_a(r_c)+α_m(r_c)＝1.05α_m(r_c)

根据后向积分公式：

由参考点位置处的气溶胶消光系数α_a(r_c)和反演得到的动态雷达比S_a，得到气溶胶消光系数的垂直廓线。

本发明相比现有技术具有以下优点：

本发明可以得到本地化的动态雷达比，即给定一组回波信号与气溶胶消光系数，就可以得到与之对应的雷达比，并且动态雷达比的反演精度远高于固定雷达比。本发明为分析激光雷达的探测数据提供了更实用和精确的方法，并且具有本地化的意义，这更加符合对数据分析的要求。

附图说明

图1为本发明的动态雷达比反演方法的流程图。

图2为2017年7月7日15：10时的气溶胶消光系数的分垂直布廓线图。

图3为不同动态雷达比Sa反演的气溶胶消光系数与实测值的相对误差图。

图4为本发明的动态雷达比及固定雷达比(雷达比Sa＝30sr)的反演值与实测值的相对误差的概率分布图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明中的技术方案的具体实施进行清楚、完整地描述。

(1)气溶胶消光系数传感器测量得到地面实测气溶胶消光系数α_a(r₀)，计算地面处雷达信号P(r₀)

由于激光雷达存在发射光束与接收望远镜视场角的不完全几何重叠区，造成近场回波信号与真实大气回波信号不一致，从而导致激光雷达近场信号存在高度为0.2698km的盲区。在高度在0-800m内的近地层大气，除一些极端天气状况，该区域大气对流较为充分，气溶胶粒子混合均匀，气溶胶消光系数基本满足随高度按指数规律递减，可表示为：

其中最低有效探测距离zm为0.2898km，气溶胶标高h₀为能见度的函数，如表1所示能见度越高h₀越大，由于考虑的为晴朗、少云时的情况，能见度均高于25km，此时对应的标高h₀为1.45km。所以最低有效位置处的消光系数为气溶胶消光系数的83.02％，即整个盲区内的气溶胶消光系数与地面处的消光系数的误差均在17％以内，并且高度越低两者相对误差越小。另外使用下午12点到下午17点对流较强时的数据，以减小盲区和地面的差别。

Table 1不同地面能见度条件下的气溶胶标高

由此可以提出一个合理的假设：在激光雷达盲区内气溶胶的消光系数保持一致。另外使用下午12点到下午17点对流较强时的数据，以减小盲区和地面的差别。

以实测气溶胶消光系数作为参考点的消光系数α_a(r₀)，但是由于激光雷达存在盲区(0.2698km)并不能直接得到与地面处消光系数相对应的回波信号。雷达方程在地面处(r₀＝0)的解为

P(r₀)＝CE(β_a(r₀)+β_m(r₀)) (5)

P(r₀)表示地面处雷达信号。由于激光雷达盲区的存在，激光雷达在最低有效位置r_m处的归一化后向散射信号为

基于激光雷达盲区内气溶胶的消光系数保持一致的假定，则上式化为

因此，可由激光雷达在r_m处的信号P(r_m)和实测气溶胶消光系数α_a(r₀)求得激光雷达在地面处的信号P(r₀)。

2017年7月7日15：10时激光雷达在最低有效探测高度r_m处的信号P(r_m)为1.5103，则此时有：

P(r₀)＝1.6332

(2)给定适合本地的预估雷达比

通过调整Fernald后向积分方程中Sa的值使得反演得到的气溶胶消光系数与实测气溶胶消光系数大小基本相同，那么此时Sa的值为适合本地的预估雷达比，此时Sa＝7sr。将此雷达比带入后向积分Fernald后向积分方程中，得到预估消光廓线。

由预估雷达比及步骤(1)得到的地面信号P(r₀)根据式(3)Fernald后向积分方程可拟合地面消光σ_a(r′₀)。

此时有：

σ_a(r′₀)＝0.0254km^-1

另外有矫正系数H为：

此时

(3)拟合地面及混合均匀处的消光系数

对(2)式两边取对数，则(2)式化为：

若存在一段混合均匀的高度，则在此范围内气溶胶后向散射系数和消光系数不随高度改变，此时ln(P(r))与高度r呈线性关系，如(10)式所示

若ln(P(r))与r高度相关，则可认为该时刻在这段距离上气溶胶均匀混合的假定是合理的。利用最小二乘法拟合ln(P(r))与高度r，并计算所有点偏离拟合直线的残差平方和，根据残差平方和的最小值确定混合均匀的高度r_u。根据步骤(2)给出的预估雷达比结合(3)式可反演得到气溶胶预估消光廓线，从而可以得到混合均匀的高度所对应的消光系数α_a(r_u)′。为了减小预估雷达造成的误差，α_a(r_u)′应乘上对应的矫正系数H，记为α_a(r_u)。

则此时有：

α_a(r_u)′＝0.0053km^-1

α_a(r_u)＝α_a(r_u)′·H＝0.0037km^-1

(4)反代求解动态雷达比

Fernald给出的激光雷达方程的解析解(3)式中共有3个未知量：α_a(r_c)——参考点的消光系数、α_a(r)——待求点处的消光系数和雷达比S_a。在此方法中，以地面(r₀)处作为待求点，混合均匀(r_u)处作为参考点。将混合均匀处的消光系数α_a(r_u)、地面处消光系数α_a(r₀)及其对应的回波信号P(r₀)带入(3)式中，即可反推出此时对应的雷达比S_a。

则此时有：

Sa＝4.55sr

(5)将反代算法进行迭代

为了进一步提高算法的反演精度，将算法进行迭代，迭代流程图如图1所示。通过对比相邻两次迭代雷达比的绝对误差，发现在第18次迭代后迭代结果趋于收敛，相邻两次迭代结果的绝对误差小于3×10^-5；在150次迭代后，相邻两次迭代结果的绝对误差小于7×10^-6，故认为迭代是收敛的。另外发现，反代法得到的气溶胶消光系数的反演值与实测值的平均相对误差，随着迭代次数出现了一个震荡减小的趋势，在第18次迭代后平均相对误差开始趋于稳定的状态，平均相对误差稳定在13.4％附近，而第2次迭代得到反演值与实测值的平均相对误差最低，为8.39％。通过分析不同迭代次数下反演值与实测值之间的误差，发现进行两次迭代的平均相对误差最小，明显优于单次反代算法(反演值与实测值的相对误差为18％)和迭代收敛后的精度，故选择进行两次迭代。

(6)利用动态雷达比反演获得气溶胶消光系数的分布廓线。

在反演气溶胶消光系数分布廓线时，我们取参考位置r_c为6km，并且在此高度有：

α_a(r_c)+α_m(r_c)＝1.05α_m(r_c) (11)

则此时有：

α_a(r_c)＝2.7×10^-4km^-1

将参考点处的气溶胶消光系数α_a(r_c)和反演得到的动态雷达比S_a，带入后向积分公式(3)中，r取不同高度即可得到气溶胶消光系数的垂直廓线。

若此时取r为最低有效探测高度0.2689km，则此时有：

α_a(0.2689)＝0.0151km^-1

2017年7月7日15：10时消光系数的分垂直布如图2所示。

(7)动态雷达比与传统固定雷达比反演结果的对比

根据《环境空气质量标准》GB 3095-2012文件中对PM2.5日均浓度限值的规定，将2016-2018年间的11000组数据分为：PM2.5＜35μg/m³(包含7400组数据)、35μg/m³＜PM2.5＜75μg/m³(包含3400组数据)，以及PM2.5＞75μg/m³(包含200组数据)，三组。雷达比取30sr为本地采取的固定雷达比，固定雷达比、动态雷达比的反演值与实测值的误差分析如表2所示。

Table 2固定雷达比、动态雷达比的反演结果与实测值的误差分析表

由表2可知本发明的反演精度明显优于固定雷达比。在不同PM2.5浓度下，动态雷达比的反演结果和实测值的绝对误差与固定雷达比的绝对误差相比均小一个量级。在PM2.5＜75μg/m³时，动态雷达比的反演结果与实测值的平均相对误差均低于15％，较固定雷达比不低于136％的平均相对误差来说，本发明的动态雷达比使反演精度得到大幅度提升。尽管随着PM2.5浓度的提高动态雷达比的反演精度有所降低，但在本发明反演精度最低的情况下(PM2.5＞75μg/m³)其反演结果与实测值的平均相对误差为38.41％，

依然较固定雷达比91.20％的相对误差低50％左右。

利用本发明对激光雷达信号进行反演，得到动态雷达比，并利用动态雷达比反演消光系数廓线，可以得到以下结论：

(1)雷达比的取值对消光系数的影响主要体现在反演值与实测值的误差上。取雷达比从1-100sr(间隔为1)反演气溶胶消光系数，反演值与实测值的平均相对误差随雷达比的变化如图3所示。可以发现反演值与实测值的相对误差随着雷达比的增加迅速增加，当雷达比从10sr增加到100sr时，相对误差由20％增加到500％。其中雷达比取30sr为本地采取的固定雷达比，本发明的动态雷达比及固定雷达比(雷达比＝30sr)的反演值与实测值的误差分析如表2所示。由表2可知动态雷达比的反演精度明显优于固定雷达比。在不同PM2.5浓度下，动态雷达比的反演结果和实测值的绝对误差与固定雷达比的绝对误差相比均小一个量级。在PM2.5＜75μg/m³时，动态雷达比的反演结果与实测值的平均相对误差均低于15％，较固定雷达比不低于136％的平均相对误差来说，本发明使反演精度得到大幅度提升。尽管随着PM2.5浓度的提高动态雷达比的反演精度有所降低，但在动态雷达比反演精度最低的情况下(PM2.5＞75μg/m³)其反演结果与实测值的平均相对误差为38.41％，依然较固定雷达比91.20％的相对误差低50％左右。

(2)由图3可知当雷达比在6-10sr之间时反演值与实测值之间的相对误差较小，在±20％之间，这与利用反代法反演得到的雷达比(平均值9.52sr)是较为接近的。可知本地雷达比的取值要略小于一般认为的雷达比的取值范围(10-100sr)，一般认为的雷达比的取值范围并不适合本地，由此可见雷达比本地化的重要性。

(3)由图4可以发现，固定雷达比反演得到的气溶胶消光系数与实测值间的相对误差更为离散，其方差为0.531，可知反演值与实测值间的相对误差存在明显的波动，即固定雷达比反演得到的气溶胶消光系数并不能很好的反应实际情况的变化，而动态雷达比反演得到的气溶胶消光系数与实测值间的相对误差的方差为0.0124，远小于固定雷达比，反演值与实测值间的相对误差明显更为固定与集中，既动态雷达比可以反应实际消光系数的改变，从而使反演值与实测值同步变化，从而达到精确反演的目的。

Claims (4)

1.一种根据实时气溶胶分布求解动态雷达比的方法，其特征在于，包括以下步骤：

步骤(1)，气溶胶消光系数传感器测量得到地面实测气溶胶消光系数α_a(r₀)；

步骤(2)，计算激光雷达在地面r₀处信号P(r₀)，包括：

激光雷达原始方程为：

其中，r表示高度；P(r)表示激光雷达归一化后向散射信号；C表示雷达常数，E表示脉冲能量，在雷达稳定工作的情况下，这两个值被视为常数；α为气溶胶消光系数；β为后向散射系数；

所述激光雷达在地面r₀处的信号P(r₀)，由激光雷达在最低有效位置r_m处的信号P(r_m)和地面实测气溶胶消光系数α_a(r₀)求得：

其中α_m为大气分子的消光系数；

步骤(3)，给定适合本地的预估雷达比，包括：

调整Fernald后向积分方程中Sa的值使得反演得到的气溶胶消光系数α_a与地面实测气溶胶消光系数α_a(r₀)大小相同，此时Sa的值为适合本地的预估雷达比；

由预估雷达比及步骤(2)得到的地面信号P(r₀)根据Fernald后向积分方程拟合地面消光系数σ_a(r′₀)，如下式所示：

式中：r₀表示地面；r_c为参考点高度，P(r_c)表示高度r_c处的激光雷达归一化后向散射信号；α_a为气溶胶消光系数；α_m为大气分子的消光系数；α_a(r_c)为r_c处气溶胶消光系数；α_m(r_c)为r_c处大气分子的消光系数；Sa为雷达比；Sm为大气分子的雷达比；

矫正系数H为：

步骤(4)，拟合地面及混合均匀处的消光系数，包括：

对式(1)两边求导，得：

利用拟合算法拟合ln(P(r))与高度r的线性关系，并计算求得这段高度上所有点偏离拟合直线的残差平方和，残差平方和最小处为所求混合均匀的高度；

根据步骤(3)给出的预估雷达比，反演得气溶胶消光廓线，从而根据所述气溶胶消光廓线的对应高度得到混合均匀处的消光系数σ_a(r′_u)；σ_a(r′_u)乘上相应的矫正系数H，即为所求混合均匀处的消光系数α_a(r_u)；

步骤(5)，反代入Fernald后向积分方程求解动态雷达比，包括：

所述Fernald给出的激光雷达方程的解析解中共有3个未知量分别为：参考点处的消光系数、待求点处的消光系数和雷达比S_a；

使用所述步骤(4)中得到的混合均匀处的消光系数α_a(r_u)作为参考点处消光；地面实测处消光系数α_a(r₀)作为待求点处消光系数，将α_a(r_u)，α_a(r₀)及待求点处对应的回波信号P(r₀)带入Fernald后向积分方程，反推出其对应的雷达比Sa；

步骤(6)，将由步骤(5)得到的雷达比Sa作为新的预估雷达比输入步骤(3)，从而将步骤(3)-(5)进行迭代；

步骤(7)，利用动态雷达比反演获得气溶胶消光系数的分布廓线。

2.根据权利要求1所述的求解动态雷达比的方法，其特征在于，所述步骤(4)包括：

所述拟合算法采用最小二乘法。

3.根据权利要求1所述的求解动态雷达比的方法，其特征在于：

所述步骤(6)中进行两次迭代。

4.根据权利要求1至3任一项所述的求解动态雷达比的方法，其特征在于，所述步骤(7)还包括：

在反演气溶胶消光系数分布廓线时，取参考位置r_c为6km，并且在此高度有：

α_a(r_c)+α_m(r_c)＝1.05α_m(r_c) (5)

根据后向积分公式：

由参考点位置处的气溶胶消光系数α_a(r_c)和反演得到的动态雷达比S_a，得到气溶胶消光系数的垂直廓线。