CN101520515A - 基于激光雷达的混合层高度的自动反演算法 - Google Patents

Abstract

本发明涉及一种基于激光雷达的混合层高度的自动反演算法，首先对激光雷达回波信号进行反演，获得气溶胶后向散射廓线图，然后设定拟合初始值，通过比较由不同初始值而得出的不同混合层高度来确定合适初始值，并使用改进后的梯度法，用误差函数拟合原始信号来获取混合层高度信息。本发明利用云高仪的气溶胶后向散射回波信号图，得出气溶胶消光系数廓线分布廓线图，通过对传统梯度法的改进用误差函数去拟合与原始信号最接近的曲线，从而能够准确求得混合层高度，具有误差小，鉴别力高，能够很好的识别所需信息的优点。

Description

基于激光雷达的混合层高度的自动反演算法

技术领域

本发明涉及激光雷达技术领域，尤其是一种基于激光雷达的混合层高度的自动反演算法

背景技术

大气混合层是对流层下部直接受地面强烈影响的一层大气。因为受非均匀下垫面的影响，不仅会对辐射收支、气象场等产生影响，而且对污染物的沉降、输送及排出等都有重要作用，如对空气污染物的分布具有支配性作用。同时，由于人类的生命活动几乎都发生在这一层内，对它的研究与天气预报、气候预测、大气物理研究、工农业生产、环境保护和航空安全等密切相关。混合层高度是评价空气质量受城市污染源、长距离输送或沉降影响的关键参数，也是检验空气污染数值模拟效果的重要参数，如何方便有效地确定混合层高度并准确监测其变化过程，对空气污染物的扩散、传输模式以及污染物预报模式都具有十分重要的意义。目前还没有直接测量混合层高度的方法，比较常用的技术有无线电探空，远程探测和参数化法。在最新的远程探测法中，利用云高仪的气溶胶后向散射廓线图推断混合层高度是一种理想的方法。确定混合层的高度有多种，如目测法、阈值法和梯度法。目测法分析有较强的主观性，视觉误差不可避免，只可以用来大致估计混合层和夹带层。阈值法是以后向散射信号低于一个预先设定的阈值作为判断基准，或者以后向散射信号超过清洁大气信号一定值为判断基准，但由于回波信号起伏很大，因此采用的廓线图一般都是经过平均处理的，这就有可能丢失一些瞬时的有用信息。梯度法的理论基础是混合层相对不变的浓度一般高于混合层之上大气的浓度，因此混合层到混合层之上就会看到一个变化，从混合层内相对较强的散射变为混合层上较小的散射。梯度方法就是通过选择后散射系数的负梯度最大值(即—dβ/dR)的地方作为混合层顶。但是因为有限差分导致噪声被放大，在实际过程中很难区分何时出现后向散射系数的最大值，为解决这一间题，通常选择时间或空间的一段连续区域进行平均，并对平均后的廓线图进行垂直平滑处理，但即使进行时空平均处理后，该方法仍有可能会把混合层之上的散射层或混合层之下的一些卷云当作混合层。

发明内容

本发明的目的是提供一种基于激光雷达的混合层高度的自动反演算法，以解决传统方法的鉴别力不高，误差大，不能较好地识别出所需信息，需要平滑处理的问题。

为了到达上述目的，本发明所采用的技术方案为：

1、一种基于激光雷达的混合层高度的自动反演算法，其特征在于：首先对激光雷达回波信号进行反演，获得气溶胶后向散射廓线图，然后设定拟合初始值，通过比较由不同初始值而得出的不同混合层高度来确定合适初始值，并使用改进后的梯度法，用误差函数拟合原始信号来获取混合层高度信息；所述算法步骤为：

(1)列出激光雷达方程，选取激光波长为905nm，对应的激光雷达方程为：P(R)＝P_tCη(R)R^-2[β_α(R)+β_m(R)]T_a ²(R)T_m ²(R)；

式中P(R)为激光雷达回波功率，P_t为激光发射功率， $C=\frac{\mathrm{c\τ}}{2}{\mathrm{TA}}_r$ 为大气探测激光雷达系统常数，c为光速，τ为激光器脉宽，T为雷达接收系统的光学效率，A_r是雷达接收单元的有效孔径，β_a(R)和β_m(R)分别为高度R处大气气溶胶和空气分子的后向散射系数； $T_a\left(R\right)=\exp [-{\∫}_0^r{\mathrm{\α}}_a(r^{,},\mathrm{\λ})d{r}^{,}]$ 为大气探测激光雷达至对应高度处大气气溶胶透过率， $T_m\left(R\right)=\exp [-{\∫}_0^r{\mathrm{\α}}_m(r^{,},\mathrm{\λ})d{r}^{,}]$ 是相应的空气分子透过率，α_a(r’，λ)和α_m(r’，λ)分别为高度R处大气气溶胶和空气分子的消光系数，η(R)是激光雷达的几何重叠因子，所述几何重叠因子总为1；

(2)对所述激光雷达方程进行距离修正，方程两边同时乘以地面到被测气溶胶粒子群高度的平方R²：

X(R)＝P(R)R²＝P_tCη(R)[β_α(R)+β_m(R)]T_a ²(R)T_m ²(R)；

(3)确定标定高度R_c、大气分子消光系数α_m(R)、气溶胶消光系数标定值α_a(R_c)、大气分子消光后向比S₁、气溶胶消光后向比S₂；

所述标定高度R_c的确定：所述R_c为近乎不含气溶胶的清洁大气层所在对流层顶附近的高度，根据资料选取R_c＝20km；

所述大气分子消光后向比S₁的确定：对于905nm波长的激光，选取S₁＝40Sr；

所述气溶胶消光后向比S₂的确定：根据瑞利散射理论得到S₂＝α_m(R)/β_m(R)＝8π/3；

所述大气分子消光系数α_m(R)的确定：先由公式 $\left\{\begin{array}{c}{\mathrm{\β}}_m\left(R\right)=1.54\×{10}^{-3}\exp (-\frac{R}{7}){\left(\frac{532}{\text{\λ}}\right)}^4\\ {\mathrm{\α}}_m\left(R\right)={\mathrm{\β}}_m\left(R\right)\×\frac{8}{3}\mathrm{\π}\end{array}\right.$ 得出905nm波长的大气分子后向散射系数β_m(R)，再由公式S₂＝α_m(R)/β_m(R)＝8π/3，求得大气分子消光系数α_m(R)；

所述气溶胶消光系数标定值α_a(R_c)的确定：由气溶胶散射比和波长之间的公式 $R\′({\mathrm{\λ}}_2,R)=1+{\left(\frac{{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1}\right)}^{4-n}[R\′({\mathrm{\λ}}_1,R)-1],$ 取n＝1，由标准值R₅₃₂＝1.01得到905nm波长的气溶胶散射比 $R_{905}^{\′}=1.05$ ，再将R_c＝20km代入公式 $\left\{\begin{array}{c}{\mathrm{\β}}_m\left(R\right)=1.54\×{10}^{-3}\exp (-\frac{R}{7}){\left(\frac{532}{\text{\λ}}\right)}^4\\ {\mathrm{\α}}_m\left(R\right)={\mathrm{\β}}_m\left(R\right)\×\frac{8}{3}\mathrm{\π}\end{array}\right.$ 得到大气分子标定系数消光值α_m(R_C)和大气分子后向散射系数标定值β_m(R_C)，将和β_m(R_C)带入气溶胶散射比公式R′＝1+β_a(R_C)/β_m(R_C)求得气溶胶后向散射系数标定值β_a(R_c)，最后由公式

${\mathrm{\α}}_a\left(R\right)=-\frac{S_1}{S_2}{\mathrm{\α}}_m\left(R\right)$

$+\frac{X\left(R\right)\·\exp \left[2(\frac{S_1}{S_2}-1){\∫}_R^{R_c}{\mathrm{\α}}_m\left(R\right)\mathrm{dR}\right]}{\frac{X\left(R_c\right)}{{\mathrm{\α}}_a\left(R_c\right)+\frac{S_1}{S_2}{\mathrm{\α}}_m\left(R_c\right)}+2{\∫}_R^{R_c}X\left(R\right)\exp \left[2(\frac{S_1}{S_2}-1){\∫}_R^{R_c}{\mathrm{\α}}_m\left(R\right)\mathrm{dR}\right]\mathrm{dR}}$ 获得消光系数标定值α_a(R_c)；

(4)获得气溶胶消光系数后，根据大气分子消光后向比S₁，获得气溶胶后向散射系数β_a(R)，绘出后向散射廓线图；

(5)用误差函数 $\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\mathrm{\π}}}\underset{0}{\overset{x}{\∫}}{e}^{-t^2}\mathrm{dt}$ 拟合原始信号，得到拟合后理想曲线表达方程 $B\left(R\right)=\frac{B_m+B_u}{2}-\frac{B_m-B_u}{2}\mathrm{erf}\left(\frac{R-H}{s}\right),$ 方程中B_m是混合层后向散射系数平均值，B_u是混合层上空后向散射系数平均值，R是高度，H是混合层平均厚度，s与夹带层的厚度有关；

(6)设 $M_1=\frac{B_m+B_u}{2},$ $M_2=\frac{B_m-B_u}{2},$ 从后向散射廓线图中进行观察，如有两个明显的梯度变化，不能马上确定M₁初始值，选择第一个梯度出现的地方p1，给p1处的M₁设定一个初始值；

(7)选择第二个梯度出现的地方p2，给p2处的M₁设定另一个初始值；

(8)比较M₁为不同初始值时取得的混合层高度，结合测量地点、该地近期气象参数等因素，选择符合该地区上空混合层高度的特征数值。

本发明利用云高仪的气溶胶后向散射回波信号图，得出气溶胶消光系数廓线分布廓线图，通过对传统梯度法的改进用误差函数去拟合与原始信号最接近的曲线，从而能够准确求得混合层高度，具有误差小，鉴别力高，能够很好的识别所需信息的优点。

附图说明

图1为气溶胶后向散射廓线图。

图2为气溶胶后向散射廓线图和拟合后的曲线。

图3为近地面处气溶胶后向散射廓线图和拟合后的曲线。

图4为混合层高度为500米处气溶胶后向散射廓线图和拟合后的曲线。

图5为实际情况与理论情况比较的混合层高度为500米处气溶胶后向散射廓线图和拟合后的曲线。

图6一天时间内测得的混合层高度变化图。

图7为算法流程图。

具体实施方式

该算法实施过程主要分为两大部分：第一步先获得气溶胶后向散射廓线图，第二步是用改进后的梯度法去拟合廓线图从而获得混合层高度等信息。

第一部分：如何获得气溶胶后向散射廓线图：

步骤1：从激光雷达方程出发，对于905nm波长，激光雷达方程为：P(R)＝P_tCη(R)R^-2[β_α(R)+β_m(R)]T_a ²(R)T_m ²(R)     (1)

步骤2：对雷达方程进行距离修正，方程两边同时乘以R²得：

X(R)＝P(R)R²＝P_tCη(R)[β_α(R)+β_m(R)]T_a ²(R)T_m ²(R)

步骤3：确定R_c，α_m(R)，α_a(R_c)，S₁，S₂五参数。

欲求气溶胶的后向散射系数，先从求气溶胶消光系数着手。

事先已知某一高度R_c处(标定高度)该波长大气气溶胶和空气分子的消光系数(标定值)，则在标定高度R_c以下该波长在各高度的气溶胶消光系数为(后向积分)：

${\mathrm{\α}}_a\left(R\right)=-\frac{S_1}{S_2}{\mathrm{\α}}_m\left(R\right)$

$+\frac{X\left(R\right)\·\exp \left[2(\frac{S_1}{S_2}-1){\∫}_R^{R_c}{\mathrm{\α}}_m\left(R\right)\mathrm{dR}\right]}{\frac{X\left(R_c\right)}{{\mathrm{\α}}_a\left(R_c\right)+\frac{S_1}{S_2}{\mathrm{\α}}_m\left(R_c\right)}+2{\∫}_R^{R_c}X\left(R\right)\exp \left[2(\frac{S_1}{S_2}-1){\∫}_R^{R_c}{\mathrm{\α}}_m\left(R\right)\mathrm{dR}\right]\mathrm{dR}}---\left(2\right)$

由上面公式可以看出，若要从激光雷达测量的回波信号P(R)中得到大气气溶胶的消光系数α_a(R)，必须事先知道标定高度R_c、大气分子消光系数α_m(R)、气溶胶消光系数标定值α_a(R_c)、大气分子消光后向比S₁和气溶胶消光后向比S₂这五个参数。

R_c：对于标定高度R_c是通过选取近乎不含气溶胶的清洁大气层所在的高度来确定，这个高度通常在对流层顶附近，结合国内外相关探空资料，在这里我们选取R_c＝20km。

S₁：S₁＝α_a(R)/β_a(R)是气溶胶消光后向散射比，它依赖于发射的激光波长、气溶胶的尺度谱分布和折射指数，对于905nm波长，取S₁＝40Sr。

S₂：S₂是大气分子消光后向散射比，根据瑞利散射理论可得

S₂＝α_m(R)/β_m(R)＝8π/3                (3)

α_m(R)：先由如下公式得出905nm波长的大气分子消光系数β_m(R)：

$\left\{\begin{array}{c}{\mathrm{\β}}_m\left(R\right)=1.54\×{10}^{-3}\exp (-\frac{R}{7}){\left(\frac{532}{\text{\λ}}\right)}^4\\ {\mathrm{\α}}_m\left(R\right)={\mathrm{\β}}_m\left(R\right)\×\frac{8}{3}\mathrm{\π}\end{array}\right.---\left(4\right)$

再由公式(3)，通过大气分子的后向散射系数β_m(R)可以求得大气分子的消光系数α_m(R)。

α_a(R_c)：先求气溶胶散射比R′(λ₂，R)。由于气溶胶散射比和波长满足下列关系：

$R\′({\mathrm{\λ}}_2,R)=1+{\left(\frac{{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1}\right)}^{4-n}[R\′({\mathrm{\λ}}_1,R)-1]---\left(5\right)$

一般取n＝1，则由标准值R₅₃₂＝1.01计算得到905nm波长的气溶胶散射比 $R_{905}^{\′}=1.05,$ 然后根据公式(4)得到大气分子标定系数消光值，将和β_m(R_C)带入气溶胶散射比公式

R′＝1+β_a(R_C)/β_m(R_C)                (6)

求得气溶胶后向散射系数标定值β_a(R_c)，由气溶胶消光后向散射比S₁获得消光系数标定值α_a(R_c)。

步骤4：获得气溶胶消光系数后，根据S₁，获得气溶胶后向散射系数，从而绘出后向散射廓线图，见图1。

第二部分：如何用改进后的梯度法得出混合层高度：

获得雷达信号的后向散射廓线图后，就可以利用改进后的梯度法(误差拟合法)求得混合层高度。这种方法的特点不是从雷达回波信号得出混合层高度，而是选择一条与原始回波信号最接近的理想曲线，从中得出混合层高度。拟合的理想曲线表达方程如下：

$B\left(R\right)=\frac{B_m+B_u}{2}-\frac{B_m-B_u}{2}\mathrm{erf}\left(\frac{R-H}{s}\right)---\left(7\right)$

其中，B_m是混合层后向散射系数平均值，B_u是混合层上空后向散射系数平均值，erf是误差函数，表达式为 $\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\mathrm{\π}}}\underset{0}{\overset{x}{\∫}}{e}^{-t^2}\mathrm{dt},$ R是高度，H是混合层平均厚度，s与夹带层的厚度有关。

步骤5：设定 $M_1=\frac{B_m+B_u}{2},$ $M_2=\frac{B_m-B_u}{2}.$

步骤6：从后向散射廓线图中观察是否存在一个以上的明显梯度变化，图2中有两个明显的梯度，不能马上确定M₁初始值，先选择第一个梯度出现的地方(p1)，则设M₁＝230，拟合后的混合层高度为325m。

步骤7：选择第二个梯度出现的地方(p2)，则设M1₁＝100，拟合后的混合层高度为575m。

步骤8：比较根据M₁不同初始值取得的混合层高度，结合测量地点、该地近期气象参数等因素，选择符合该地区上空混合层高度的特征数值575m，舍弃325m。

采用本方法在多种回波信号下进行拟合的结果分析：

图3—5是云高仪测得的三种情况的气溶胶后向散射廓线图和拟合后的曲线，图3中在近地面处的后向散射系数值较小，可能是因为测量时间为早晨6点，无污染源，大气较清洁，尽管这种方法在底部廓线和原始信号的拟合并不完美，但不影响混合层高度的确定，混合层高度较高；图4的混合层在500米左右，偏低，但利用仍能较好地判断出混合层，在约1400米处后向散射系数突然增大，拟合分析结果为一375米厚的云层；从图5可以看出，除一些噪声外，廓线拟合与实际情况较吻合。图6为在一天内测得的混合层高度变化图，从下午约一点左右混合层高度开始下降。我们对此给出的解释是：太阳未升起前，由近地面热对流引起的混合层高度较低，太阳升起后，在晴朗干净的天气情况下，地面开始存储大量的热能，地面上空能保持的热对流不断积累，而当大气层面处于不稳定状态时，热对流又与空气的垂直浮动有关，剧烈的热对流加快空气垂直浮动，混合层高度变高；又因为热对流造成自由大气向下浮动，同时云的反馈作用使地面接收的热辐射通量减少，地面温度偏低，大气湿度较大，因此混合层高度在某一时间开始下降，夹带层厚度增加。总体来说，由近地层对流混合形成的城市混合层高度具有大气混合层典型的日变化特点，即早晚比较低，白天逐渐升高，在某一时间又开始下降，混合层高度变化一般在0.5—2km范围之内。

研究表明，利用本发明的方法进行混合层的探测是有效可行的，即使在一些较困难的情况下，如混合层高度较低，混合层和其上层后向散射信号较相近、难以区分转折点，在混合层中存在强散射气流团等时，这种方法仍然有效。

Claims (1)

1、一种基于激光雷达的混合层高度的自动反演算法，其特征在于：首先对激光雷达回波信号进行反演，获得气溶胶后向散射廓线图，然后设定拟合初始值，通过比较由不同初始值而得出的不同混合层高度来确定合适初始值，并使用改进后的梯度法，用误差函数拟合原始信号来获取混合层高度信息；所述算法步骤为：

(1)列出激光雷达方程，选取激光波长为905nm，对应的激光雷达方程为：P(R)＝P_tCη(R)R^-2[β_α(R)+β_m(R)]T_a ²(R)T_m ²(R)；

式中P(R)为激光雷达回波功率，P_t为激光发射功率， $C=\frac{\mathrm{c\τ}}{2}{\mathrm{TA}}_r$ 为大气探测激光雷达系统常数，c为光速，τ为激光器脉宽，T为雷达接收系统的光学效率，A_r是雷达接收单元的有效孔径，β_a(R)和β_m(R)分别为高度R处大气气溶胶和空气分子的后向散射系数； $T_a\left(R\right)=\exp [-{\∫}_0^r{\mathrm{\α}}_a(r^{,},\mathrm{\λ})d{r}^{,}]$ 为大气探测激光雷达至对应高度处大气气溶胶透过率， $T_m\left(R\right)=\exp [-{\∫}_0^r{\mathrm{\α}}_m(r^{,},\mathrm{\λ})d{r}^{,}]$ 是相应的空气分子透过率，α_a(r’，λ)和α_m(r’，λ)分别为高度R处大气气溶胶和空气分子的消光系数，η(R)是激光雷达的几何重叠因子，所述几何重叠因子总为1；

(2)对所述激光雷达方程进行距离修正，方程两边同时乘以地面到被测气溶胶粒子群高度的平方R²：

X(R)＝P(R)R²＝P_tCη(R)[β_α(R)+β_m(R)]T_a ²(R)T_m ²(R)；

(3)确定标定高度R_c、大气分子消光系数α_m(R)、气溶胶消光系数标定值α_a(R_c)、大气分子消光后向比S₁、气溶胶消光后向比S₂；

所述标定高度R_c的确定：所述R_c为近乎不含气溶胶的清洁大气层所在对流层顶附近的高度，根据资料选取R_c＝20km；

所述大气分子消光后向比S₁的确定：对于905nm波长的激光，选取S₁＝40Sr；

所述气溶胶消光后向比S₂的确定：根据瑞利散射理论得到S₂＝α_m(R)/β_m(R)＝8π/3；

所述大气分子消光系数α_m(R)的确定：先由公式 $\left\{\begin{array}{c}{\mathrm{\β}}_m\left(R\right)=1.54\×{10}^{-3}\exp (-\frac{R}{7}){\left(\frac{532}{\mathrm{\λ}}\right)}^4\\ {\mathrm{\α}}_m\left(R\right)={\mathrm{\β}}_m\left(R\right)\×\frac{8}{3}\mathrm{\π}\end{array}\right.$ 得出905nm波长的大气分子后向散射系数β_m(R)，再由公式S₂＝α_m(R)/β_m(R)＝8π/3，求得大气分子消光系数α_m(R)；

所述气溶胶消光系数标定值α_a(R_c)的确定：由气溶胶散射比和波长之间的公式 $R\′({\mathrm{\λ}}_2,R)=1+{\left(\frac{{\mathrm{\λ}}_2}{{\mathrm{\λ}}_1}\right)}^{4-n}[R\′({\mathrm{\λ}}_1,R)-1],$ 取n＝1，由标准值R₅₃₂＝1.01得到905nm波长的气溶胶散射比 $R_{905}^{\′}=1.05,$ 再将R_c＝20km代入公式 $\left\{\begin{array}{c}{\mathrm{\β}}_m\left(R\right)=1.54\×{10}^{-3}\exp (-\frac{R}{7}){\left(\frac{532}{\mathrm{\λ}}\right)}^4\\ {\mathrm{\α}}_m\left(R\right)={\mathrm{\β}}_m\left(R\right)\×\frac{8}{3}\mathrm{\π}\end{array}\right.$ 得到大气分子标定系数消光值α_m(R_C)和大气分子后向散射系数标定值β_m(R_C)，将

和β_m(R_C)带入气溶胶散射比公式R′＝1+β_a(R_C)/β_m(R_C)求得气溶胶后向散射系数标定值β_a(R_c)，最后由公式 ${\mathrm{\α}}_a\left(R\right)=-\frac{S_1}{S_2}{\mathrm{\α}}_m\left(R\right)$

$+\frac{X\left(R\right)\·\exp \left[2(\frac{S_1}{S_2}-1){\∫}_R^{R_c}{\mathrm{\α}}_m\left(R\right)\mathrm{dR}\right]}{\frac{X\left(R_c\right)}{{\mathrm{\α}}_a\left(R_c\right)+\frac{S_1}{S_2}{\mathrm{\α}}_m\left(R_c\right)}+2{\∫}_R^{R_c}X\left(R\right)\exp \left[2(\frac{S_1}{S_2}-1){\∫}_R^{R_c}{\mathrm{\α}}_m\left(R\right)\mathrm{dR}\right]\mathrm{dR}}$ 获得消光系数标定值α_a(R_c)；

(4)获得气溶胶消光系数后，根据大气分子消光后向比S₁，获得气溶胶后向散射系数β_a(R)，绘出后向散射廓线图；

(5)用误差函数 $\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\mathrm{\π}}}\underset{0}{\overset{x}{\∫}}{e}^{-t^2\mathrm{dt}}$ 拟合原始信号，得到拟合后理想曲线表达方程 $B\left(R\right)=\frac{B_m+B_u}{2}-\frac{B_m-B_u}{2}\mathrm{erf}\left(\frac{R-H}{s}\right),$ 方程中B_m是混合层后向散射系数平均值，B_u是混合层上空后向散射系数平均值，R是高度，H是混合层平均厚度，s与夹带层的厚度有关；

(6)设 $M_1=\frac{B_m+B_u}{2},$ $M_2=\frac{B_m-B_u}{2},$ 从后向散射廓线图中进行观察，如有两个明显的梯度变化，不能马上确定M₁初始值，选择第一个梯度出现的地方p1，给p1处的M₁设定一个初始值；

(7)选择第二个梯度出现的地方p2，给p2处的M₁设定另一个初始值；

(8)比较M₁为不同初始值时取得的混合层高度，结合测量地点、该地近期气象参数等因素，选择符合该地区上空混合层高度的特征数值。

Images