基于激光雷达探测的晴空条件下飞机尾流特征参数反演方法

技术领域

本发明属于航空安全技术领域，涉及一种常见的晴空条件下飞机尾流特征参数(即环量和涡心位置)的估计方法。

背景技术

飞机尾流是飞机飞行时在其后方产生的一对反向旋转的强烈漩涡，其对后续飞机的飞行安全具有重大影响，可能使其产生抖动、颠簸，甚至失去控制。导致265人死亡的2001年美航587飞机在纽约肯尼迪机场坠毁即是因尾流所致，这也是美国至今除911之外最重大的航空事故。为了减小此类事故发生的概率，国际民航组织(ICAO)制定了两架飞机起飞或降落的最小安全距离准则。该准则虽然可以在很大程度上保证飞机免受尾流的影响，但往往过于保守，限制了机场跑道的飞机起降密度，这也是导致飞机晚点的一个重要原因。

因此，为在保证飞机飞行安全的前提下有效提升机场容量，亟需发展飞机尾流的实时探测方法，通过估计飞机尾流特征参数(环量和涡心位置)对飞机尾流的危害进行评估，并据此实时调整飞机的起降间隔。当前欧洲的SESAR(欧洲单一天空空中交通管理研究)计划和美国的NextGen(美国下一代空中交通管理系统)计划均设了专题对尾流的实时探测问题进行研究。

激光雷达是晴空条件(干燥大气，大气不含雨、雾、雪等降水粒子)下飞机尾流探测最常用的方式。现有的基于激光雷达探测的尾流特征参数反演方法总体来说可以分为两大类：

第一类属于速度包络匹配算法，通常首先根据速度包络估计涡心位置，再基于速度模型估计环量值。比如一种方法是在频谱中设定阈值得到速度包络，从中估计涡心位置，最后通过比对现有的尾流速度模型得到左右漩涡环量。还有方法就是从波束视线方向(LOS)最大速度估计环量的方法，但在每个波束视线上只使用了径向速度最大的那一个距离单元处的数据，而舍弃了其他距离单元携带的信息，方法的鲁棒性和数据的利用率还可进一步提升。另外，涡心位置的估计对这类方法非常重要，其可能会给环量估计的鲁棒性和正确率造成很大影响。

第二类属于模板匹配算法，这类算法通常建立尾流特征参数与测量变量之间的数学模型，并通过大量探测单元的数据组合求解尾流的特征参数。比如一种方法是通过比较频谱数据和漩涡模型，分别建立了对称和非对称飞机尾流的特征参数与测量变量的两个极大似然函数模型。还有方法就是首先利用包络速度得到涡心位置粗略的估计，接下来基于环量是自由参数的假设精确地计算涡心位置，方法的鲁棒性得到提高。但是总体来说，这一类算法要求的计算量较大。

发明内容

本发明的目的在于提出一种基于激光雷达探测的晴空条件下飞机尾流特征参数反演方法，该方法能够避免涡心位置估计给下一步环量估计造成的影响，同时提升了数据利用率。同时本发明采用了最优化方法，将参数反演问题转变为全局优化和参数甄别的问题，减少了对计算量的需求。

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

基于激光雷达探测的晴空条件下飞机尾流特征参数反演方法，包括以下步骤：

步骤一、利用激光雷达探测系统，获取激光雷达扫描区域内不同激光雷达探测单元处的多普勒速度。

在机场跑道的一侧设置激光雷达，以激光雷达所在位置为坐标原点，建立三维直角坐标系XYZ，其中X轴表示到激光雷达的水平距离，Y轴表示到激光雷达的竖直高度差，Z轴为机场跑道方向，激光雷达在与机场跑道正交的XOY平面上由下至上，再由上至下进行来回扫描，激光雷达扫描的区域为激光雷达扫描区域，激光雷达在x0y平面内的扫描的角度范围为[α_min,α_max]，α_min,α_max为预设值。如2014年香港国际机场开展了一系列尾流探测实验，用于实验的激光雷达波束仰角的边界分别为：α_min＝0.83°，α_max＝10.77°。

激光雷达扫描区域内设有L个驻留的角度。激光雷达在来回扫描过程中，在每个驻留的角度上发射一束用于探测的主振激光，在被主振激光照亮的窄扇形区域中顺着主振激光的传播方向划分为多个连续的长度尺寸相同的小单元格，这些小单元格称为激光雷达探测单元，通过激光雷达探测单元内气溶胶粒子的后向散射的回波光得到不同激光雷达探测单元处的多普勒速度。

对于第m个驻留的角度上的第i个激光雷达探测单元，(x_i,y_i)为第i个激光雷达探测单元其中心位置处坐标，R_L、R_R分别为第i个激光雷达探测单元其中心位置处坐标与尾流的左右涡心之间的距离。

步骤二利用激光雷达扫描区域上属于非尾流区域的各激光雷达探测单元的多普勒速度求解背景风速度。

将激光雷达在XOY平面内的扫描区域初步划分为尾流区域和远离尾流的非尾流区域：以机场跑道中心线为参考，在激光雷达扫描区域上，在机场跑道中心线两侧且距离机场跑道中心线的水平距离225m以内的区域认为是尾流区域，尾流区域外前后200m的区域认为是非尾流区域。

利用激光雷达扫描区域上属于非尾流区域的各激光雷达探测单元的多普勒速度，按下式计算求解背景风速度：

其中为在非尾流区域内的一主振激光的路径上的第j个激光雷达探测单元处的多普勒速度，j＝1,2,3…N；U为矩阵α_j∈[α_min,α_max]为激光雷达发出的主振激光穿过第j个激光雷达探测单元时的仰角，N为非尾流区域中激光雷达探测单元的总数。

步骤三：利用属于尾流区域的各激光雷达探测单元的多普勒速度求解涡心位置初值；

利用激光雷达扫描区域上属于尾流区域中的各激光雷达探测单元的多普勒速度求得尾流区域中竖直方向(即Y轴方向)上所有激光雷达探测单元的速度梯度的分布，其中x，y分别表示尾流区域中各个激光雷达探测单元其中心位置的横、纵坐标值，V_D(x,y)表示尾流区域内中心位置坐标为(x,y)的激光雷达探测单元处测得的多普勒速度；

取尾流区域中速度梯度为正值且最大(即)的两个激光雷达探测单元其中心位置所在的点，设为点A和点B；取尾流区域中速度梯度为负值且最小(即)的两个激光雷达探测单元其中心位置所在的点，设为点C和点D；将这4个点按组合中一个点梯度正最大，另一个点梯度负最小的形式组成4种组合，即组合AC,AD,BC和BD。

引入方向约束：一个组合中速度梯度正最大的点(A点或B点)应比负最小的点(C点或D点)更远离激光雷达。

引入距离约束：一个组合中两个点的间距在20m到90m范围内，两个点的水平距离在 20m到90m范围内，且两个点的高度差小于30m。

通过引入上述两个约束来判别这4个组合中是否存在有效的左右涡心位置初值的组合：

如果没有组合同时满足这两个约束，则判定该尾流区域中没有尾流；反之，如果有不少于一个组合同时满足上面两个约束条件，则比较这些组合中两个点对应的速度梯度值，将满足约束且梯度最大的组合作为涡心位置初值。

步骤四：利用最优化方法求解飞机尾流特征参数与多普勒速度、背景风速度以及飞机尾流理论速度之间的目标方程，得到尾流特征参数解的集合。

其中，上述约束函数最小的值所对应的解称为最优解，M为尾流区域中激光雷达探测单元的总数；为在尾流区域内第i个激光雷达探测单元处的多普勒速度，i＝1,2,3…M；Γ_L和Γ_R分别代表尾流的左右漩涡的涡环量，O_L和O_R分别代表尾流的左右漩涡的涡心位置，Γ_L、Γ_R、>L和O_R为待求的尾流特征参数；(x_i,y_i)为尾流区域中第i个激光雷达探测单元中心位置处坐标。F(·)为表征背景风的径向速度分量尾流理论速度在径向上的投影多普勒速度与待求的尾流特征参数的约束关系的目标方程：

其中，为第i个激光雷达探测单元处的多普勒速度；为背景风投影到穿过第i个激光雷达探测单元的激光波束上的速度分量；为尾流速度投影到穿过第i个激光雷达探测单元的激光波束上的分量；α_i＝arctan(y_i/x_i)为激光雷达扫描波束穿过第i个激光雷达探测单元时的仰角，V_wind为步骤二中由式(1)得到的背景风速度，为Hallock-Burnham(哈洛克-伯纳姆)模型对应的飞机尾流理论速度：

其中，r_c是涡核半径(一般取r_c＝0.052b₀，其中b₀为左右漩涡涡心间的涡间距)，(x_i,y_i)为尾流区域中第i个激光雷达探测单元的坐标，r_L,i，r_R,i是尾流区域中第i个激光雷达探测单元到左右涡心的距离，O_L和O_R为激光雷达在与机场跑道正交的扫描区域上由下至上扫描一次(仰角变化范围为α_min～α_max)的时间内的实时涡心位置(考虑背景风和下沉风)：

其中，t是涡心的下沉时间,O_L0,O_R0为第三步求解得到的涡心位置初值，上标x,y分别表示左右涡心在以激光雷达为原点的XOY坐标系中横坐标和纵坐标的值。

步骤五：对尾流特征参数解的集合融合甄别得到尾流特征参数解；

对激光雷达在与机场跑道正交的扫描区域上由下至上或由上至下扫描一次所扫描的激光雷达扫描区域内的尾流区域求解式(2)，对应可得到包括最优解在内的若干组优化解。

由于噪声等非理想因素的存在，最优解可能并非最适于真实情况的解，需要从解的集合中甄选最佳参数。

激光雷达在与机场跑道正交的扫描区域上由下至上后再由上至下这样来回扫描多次，对于每次扫描的激光雷达扫描区域，利用Hough变换分别找到最多最优解对应的左右涡心位置所在的直线Line_L和Line_R，其中Line_L是最多左涡心所在的直线，Line_R是最多右涡心所在的直线，Line_L，Line_R可看作左右涡心位置演化比较合理的直线轨迹。

然后在每次扫描的激光雷达扫描区域内的各优化解中，按左右涡心位置到对应直线的距离之和(即左涡心到Line_L的距离与右涡心到Line_R的距离的距离和)进行排序，取距离和最小的优化解中的左右涡环量，涡心位置作为该扫描区域上反演得到的尾流特征参数。

本发明具有以下技术效果：

本发明提供的该方法能够避免涡心位置估计给下一步环量估计造成的影响，同时提升了数据利用率。同时本发明采用了最优化方法，将参数反演问题转变为全局优化和参数甄别的问题，减少了对计算量的需求。

附图说明

图1是激光雷达扫描的三维示意图；

图2是激光雷达、机场跑道和尾流之间的几何关系图；

图3是本发明流程图；

图4是在一具体实施例中激光雷达测得的多普勒速度在扫描区域的分布图；

图5为在一具体实施例中利用激光雷达测得的多普勒数据反演得到的环量结果图；

图6为在一具体实施例中利用激光雷达测得的多普勒数据反演得到的涡心位置随时间的变化图。

具体实施方式

下面将结合本发明实施例图中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，作进一步详细说明。

参照图1，图1是激光雷达扫描的三维示意图。在机场跑道的一侧设置激光雷达，激光雷达在与机场跑道正交的平面上来回扫描，来回扫描的区域为激光雷达扫描区域。以激光雷达所在位置为坐标原点，建立三维直角坐标系XYZ，其中X轴表示激光雷达的水平方向，Y 轴表示激光雷达的竖直方向，Z轴为跑道方向，激光雷达在XOY平面上由下至上，再由上至下进行来回扫描，激光雷达由下至上完成角度范围为α_min～α_max的扇形区域扫描后，转而由上至下扫描角度范围为α_max～α_min的扇形区域，由此循环，得到多次激光雷达扫描数据。

激光雷达由下至上或由上之下的每次扫描过程中均设有L个驻留的角度。激光雷达扫描时，在每个驻留的角度上发射一束用于探测的主振激光(如图1中画出的五个驻留角上的波束)，在主振激光的路径上依次分布着多个激光雷达探测单元(图1中扫描区域上的小方形表示激光雷达的探测单元，激光雷达探测单元是激光雷达最小可分辨单元，如图2中的p^th，q^th和i^th探测单元；激光雷达在每个驻留角上发射一束用于探测的主振激光，在被主振激光照亮的窄扇形区域中，顺着主振激光的传播方向可以划分为一个一个连续的长度尺寸相同的小单元格。激光雷达可以得到这些小单元格内的多普勒速度信息，因此将这些小单元格称为激光雷达探测单元。因为扫描区域在XOY平面上，所以激光雷达探测单元也位于XOY平面)。通过激光雷达探测单元内气溶胶粒子的后向散射的回波光得到不同激光雷达探测单元处的多普勒速度。

图2所示为激光雷达、机场跑道和尾流之间的几何关系。尾流的左右涡环量分别为Γ_L和Γ_R，两涡心位置分别为O_L和O_R，涡心间距b₀。对于第m个驻留的角度上的第i个激光雷达探测单元，(x_i,y_i)为第i个激光雷达探测单元其中心位置处坐标，R_L、R_R分别为第i个激光雷达探测单元其中心位置处坐标与尾流的左右涡心之间的距离。

激光雷达在与机场跑道正交的平面上由下而上完整扫描一次，由下而上完整扫描一次的仰角变化范围为α_min～α_max，获得扫描过程中所有激光雷达探测单元的多普勒数据后，得到多普勒速度在激光雷达扫描平面上的分布。完成本次扫描后，激光雷达在扫描区域上转而由上至下完成下一次扫描，由上至下完整扫描一次的仰角变化范围为α_max～α_min，得到此次扫描过程中多普勒速度在激光雷达扫描区域上的分布。依次规律，通过激光雷达不断地来回，上下扫描，可以得到多次扫描数据，即可得到各次扫描过程中不同时刻多普勒速度在激光雷达扫描区域上的分布。

图3为本发明利用激光雷达测量数据对尾流特征参数进行反演，融合，甄别的流程图。具体实施步骤如下：

第一步：利用激光雷达扫描区域上属于非尾流区域的各激光雷达探测单元的多普勒速度求解背景风速度。

在晴空条件下，采用位于跑道一侧且扫描平面与跑道垂直的高分辨率激光雷达上下扫描飞机尾流。

将激光雷达在XOY平面内的扫描区域初步划分为尾流区域和远离尾流的非尾流区域：以机场跑道中心线为参考，在激光雷达扫描区域上，在机场跑道中心线两侧且距离机场跑道中心线的水平距离225m以内的区域认为是尾流区域(如图3中Region①所示)，尾流区域外前后200m的区域认为是非尾流区域(如图3中Region②所示)。

激光雷达在与机场跑道正交的平面上由下而上完整扫描一次，获得扫描过程中所有激光雷达探测单元的多普勒数据后，得到多普勒速度在扫描平面上的分布。利用激光雷达扫描区域上属于非尾流区域的各激光雷达探测单元的多普勒速度，按下式计算求解背景风速度：

其中为在非尾流区域内的一主振激光的路径上的第j个激光雷达探测单元处的多普勒速度，j＝1,2,3…N；U为矩阵α_j∈[α_min,α_max]为激光雷达发出的主振激光穿过第j个激光雷达探测单元时的仰角，N为非尾流区域中激光雷达探测单元的总数。

第二步、利用属于尾流区域的各激光雷达探测单元的多普勒速度求解涡心位置初值。

利用激光雷达扫描区域上属于尾流区域中的各激光雷达探测单元的多普勒速度求得尾流区域中竖直方向(即Y轴方向)上所有激光雷达探测单元的速度梯度的分布，其中x，y分别表示尾流区域中各个激光雷达探测单元其中心位置的横、纵坐标值，V_D(x,y)表示尾流区域内中心位置坐标为(x,y)的激光雷达探测单元处测得的多普勒速度；

取尾流区域中速度梯度为正值且最大(即)的两个激光雷达探测单元其中心位置所在的点，设为点A和点B；取尾流区域中速度梯度为负值且最小(即)的两个激光雷达探测单元其中心位置所在的点，设为点C和点D；将这4个点按组合中一个点梯度正最大，另一个点梯度负最小的形式组成4种组合，即组合AC,AD,BC和BD。

引入方向约束：一个组合中速度梯度正最大的点(A点或B点)应比负最小的点(C点或D点)更远离激光雷达。

引入距离约束：一个组合中两个点的间距在20m到90m范围内，两个点的水平距离在 20m到90m范围内，且两个点的高度差小于30m。

通过引入上述两个约束来判别这4个组合中是否存在有效的左右涡心位置初值的组合：

如果没有组合同时满足这两个约束，则判定该尾流区域中没有尾流；反之，如果有不少于一个组合同时满足上面两个约束条件，则比较这些组合中两个点对应的速度梯度值，将满足约束且梯度最大的组合作为涡心位置初值。

第三步、利用最优化方法求解飞机尾流特征参数与多普勒速度、背景风速度以及飞机尾流理论速度之间的目标方程，得到尾流特征参数解的集合。

其中，上述约束函数最小的值所对应的解称为最优解，M为尾流区域中激光雷达探测单元的总数；为在尾流区域内第i个激光雷达探测单元处的多普勒速度，i＝1,2,3…M；Γ_L和Γ_R分别代表尾流的左右漩涡的涡环量，O_L和O_R分别代表尾流的左右漩涡的涡心位置，Γ_L、Γ_R、>L和O_R为待求的尾流特征参数；(x_i,y_i)为尾流区域中第i个激光雷达探测单元中心位置处坐标。F(·)为表征背景风的径向速度分量尾流理论速度在径向上的投影多普勒速度与待求的尾流特征参数的约束关系的目标方程：

其中，为第i个激光雷达探测单元处的多普勒速度；为背景风投影到穿过第i个激光雷达探测单元的激光波束上的速度分量；为尾流速度投影到穿过第i个激光雷达探测单元的激光波束上的分量；α_i＝arctan(yi/x_i)为激光雷达扫描波束穿过第i个激光雷达探测单元时的仰角，V_wind为步骤二中由式(1)得到的背景风速度，为Hallock-Burnham(哈洛克-伯纳姆)模型对应的飞机尾流理论速度：

其中，r_c是涡核半径(一般取r_c＝0.052b₀，其中b₀为左右漩涡涡心间的涡间距)，(x_i,y_i)为尾流区域中第i个激光雷达探测单元的坐标，r_L,i，r_R,i是尾流区域中第i个激光雷达探测单元到左右涡心的距离，O_L和O_R为激光雷达在与机场跑道正交的扫描区域上由下至上扫描一次(仰角变化范围为α_min～α_max)的时间内的实时涡心位置(考虑背景风和下沉风)：

其中，t是涡心的下沉时间,O_L0,O_R0为第二步求解得到的涡心位置初值，上标x,y分别表示左右涡心在以激光雷达为原点的XOY坐标系中横坐标和纵坐标的值。

第四步：对尾流特征参数解的集合融合甄别得到尾流特征参数解

对激光雷达每一次由上至下或由下至上扫描过程中所扫描的扫描区域内的尾流区域求解式(2)，对应可得到包括最优解在内的若干组优化解。

由于噪声等非理想因素的存在，最优解可能并非最适于真实情况的解，需要从解的集合中甄选最佳参数。

首先对于激光雷达每一次由上至下或由下至上扫描过程中所扫描的扫描区域，利用Hough 变换分别找到最多最优解对应的左右涡心位置所在的直线Line_L和Line_R，其中Line_L是最多左涡心所在的直线，Line_R是最多右涡心所在的直线，Line_L，Line_R可看作左右涡心位置演化比较合理的直线轨迹。然后在每次扫描的激光雷达扫描区域内的各优化解中，按左右涡心位置到对应直线的距离之和(即左涡心到Line_L的距离与右涡心到Line_R的距离的距离和)进行排序，取距离和最小的优化解中的左右涡环量，涡心位置作为该扫描区域上反演得到的尾流特征参数。

本发明所提供的方法在香港国际机场的飞机尾流激光雷达探测数据上得到成功应用，图4 示出2014年8月3日07:17:25激光雷达测得的多普勒速度在扫描区域的分布图。其中，横轴是激光雷达探测单元到激光雷达的水平距离，纵轴是激光雷达探测单元和激光雷达的高度差，颜色分布代表多普勒速度大小，颜色越黑表示该探测单元处的多普勒速度绝对值越大。两条垂直虚线之间的区域为尾流区域，黑色垂直虚线和黑色垂直实线之间的区域为非尾流区域。空心三角形符号(△)和圆形符号(○)分别为按第二步估计得到的左右涡心位置初值O_L0,O_R0。可以看出尾流特殊的双对顶结构，这有助于较好地实现尾流涡心位置估计。

图5为本发明利用2014年8月3日07:17:25-07:18:34的激光雷达测量数据反演得到的环量结果图。其中，横轴为从尾流产生到扫描的时间，纵轴为环量强度。三角形代表左涡环量，圆形代表右涡环量。从这两个图可以看出，反演环量虽然呈现一定的振荡，但是随着尾流演化时间的增长，整体呈现一个降低的趋势，与实际情况吻合。

图6为本发明利用2014年8月3日07:17:25-07:18:34的激光雷达测量数据反演得到的涡心位置随时间的变化图。其中横轴为探测单元到激光雷达的水平距离，纵轴为探测单元离地面的高度(已考虑了激光雷达的高度)。三角形代表左涡位置，圆形代表右涡位置。从图中可以看出随着尾流演化时间的增长，左右涡心在高度上都有所下降，但是左涡心下降更快，说明在尾流演变过程中，两个涡心不再处于同一个高度。