激光雷达探测中基于贝叶斯网络的飞机尾流反演

尾流是飞机飞行时的必然产物，对航空安全有重大威胁，还限制着航空效率和容量的提升，飞机尾流涡核的精准辨识是动态缩减尾流间隔的前提，目前晴空尾涡探测的主要工具是相干多普勒激光雷达（CDL）。针对使用激光雷达进行飞机尾流探测中受限于雷达时空分辨率和背景风场的影响导致的识别和反演尾流关键参数误差较大这一问题，提出一种在激光雷达探测数据基础上基于贝叶斯网络（BN）和均方误差（MSE）构建的尾涡参数反演模型。搭建大气背景风和湍流环境，并将其叠加到模拟的尾涡速度场上，得到用于训练模型的仿真数据集。实验结果表明：所提算法能够得到误差较小的参数反演结果（仿真算例中涡核位置偏差在2 m以内，环量偏差在5%以内）；在实际算例中，所提算法与传统算法相比，反演速度场均方误差显著降低（平均超过50%）。本研究可用于机场实时尾涡监测，对尾涡间隔标准制定有重要意义。

Abstract

Aircraft wake, an inevitable byproduct of aircraft flight, poses a major threat to aviation safety and limits the improvement of aviation efficiency and capacity. Accurate identification of aircraft wake vortex nuclei is a prerequisite for dynamically reducing wake intervals, and coherent Doppler Lidar (CDL) is the main tool for clear-air wake detection. To address the significant errors in identifying and inverting key parameters of aircraft wake turbulence caused by the limitations of CDL spatiotemporal resolution and background wind field effects, this study proposes a wake vortex parameter inversion model based on Bayesian network (BN) and mean squared error (MSE) using CDL detection data. An atmospheric background wind and turbulence environment are built and superimposed onto the simulated wake velocity field to obtain a simulation dataset for training the model. The results show that the proposed model can obtain parameter inversion results with small errors (within 2 meters deviation of the vortex core position and within 5% deviation of the ring volume in the simulated case) at an acceptable computational level. In actual cases, the mean squared error of the inversion velocity field is significantly reduced (more than 50% on average) compared with the conventional algorithm. This research can be used for real-time monitoring of wake vortices at airports and is of great significance for the development of wake interval standards.

1　引言

尾流是飞机升力的副产物，若处置不当会对尾随其后的飞机飞行安全产生威胁，现行间隔标准过于保守，使得跑道和空域容量成为限制飞行流量增长的瓶颈^［1］。未来尾涡间隔标准将向动态化方向发展，在保障飞行安全的前提下，对尾涡进行探测并识别涡核具体位置有助于缩短间隔并提高运行效率^［2-5］。

在多普勒激光雷达探测数据的基础上对飞机尾涡的精确识别已成为这一领域的研究重点^［6］，近年来在尾涡雷达探测及其数据处理方面国内外主要进行了以下研究。2003年，Holzäpfel等^［7］应用机场多普勒激光雷达开展了多次尾涡探测实验，提出一种基于激光雷达测量数据估计尾涡环量的方法。Köpp等^［8-9］应用2 μm多普勒雷达探测尾涡，提出一种基于雷达回波数据的多普勒频谱、径向速度等数据的处理算法。日本电子导航研究所的Komatsubara等^［10］和法国Thales公司的Barbaresco等^［11］分别在机场应用多普勒激光雷达进行不同民航机型的尾涡探测实验。2008年，德国宇航中心在实验中使用发烟装置，利用机载2 μm连续波激光雷达获得了较好的尾涡探测结果^［12］。Smalikho等^［13］利用雷达探测出风场和尾涡相关参数，研究了风和大气湍流对尾涡的影响。在此基础上，2015年他们又基于1.55 μm多普勒激光雷达实测数据提出了径向速度法（RV）反演尾涡特征参数，该方法在提取较低信噪比的大飞机回波数据方面性能较好，适用于稳定大气和弱湍流条件^［14］。Darracq等^［15］利用机载雷达探测尾涡设计了利用轴向速度探测尾涡特征参数的算法。Ehlers等^［16］同样在机载雷达探测的基础上，提出利用尾涡特征进行尾涡识别的算法。Hallermeyer等^［17］提出一种将原有的速度包络法和参数化方法相结合的新的尾涡特征参数反演算法，对不同机型和不同天气状态下产生的尾涡进行了模拟，最终得到的环量反演均方根误差小于20%。Smalikho等^［18-20］分析不同飞机机型的尾涡演化过程，考虑近地面效应的影响，并据此进一步确定了不同机型和不同湍流强度条件下利用RV法进行尾涡特征参数反演的局限性。王筱晔等^［21］应用多普勒激光雷达在四川开展了尾涡观测实验，提出一种基于谱宽的尾涡识别算法。赵丽雅等^［22-23］考虑尾涡下沉及消耗，提出一种基于雷达回波的动态尾涡特征参数计算方法。香港天文台的Hon等^［24-25］在香港机场开展了一系列激光雷达探测进近阶段飞机尾涡的实验，并取得了初步成果。Gao等^［26-30］提出一种基于单部激光雷达的三维风场局地/全局变分反演算法，结合风场线性切变和最小二乘拟合法构建参数化尾流行为预测模型，解决了经典尾流预测模型气象环境参数未随时间演化实时调整的问题。

目前许多对尾涡结构的研究建立在反演所得速度场的基础上，而在实际探测中用激光雷达对尾涡流场进行扫描时，时空分辨率受到限制，雷达距离门很难总是恰好与涡核位置重合。本文研究旨在基于可实现的雷达时空分辨率，进一步提高尾流参数反演的准确性，提出一种基于贝叶斯网络和均方误差寻优法的尾涡参数反演模型，结合仿真分析评估模型实用性，并展示了该模型在实际算例中的效果。结果表明，所提尾涡识别模型与传统的梯度法和极差法等方法相比，能得到更准确的结果，具有较高的实用价值。

2　尾涡定位及强度反演模型

2.1　基于贝叶斯网络的涡核位置识别方法

贝叶斯网络是一种可以将先验知识和已有数据结合起来学习因果关系的概率模型，是目前不确定知识和概率推理领域中最有效的方法之一，它以图形化的方式直观地表达各变量的联合概率分布，可以应对噪声和数据缺失等情况，能在大量样本学习过程中不断优化网络参数以得到更好的推理结果，具有很强的环境适应性^［31］。目前已有研究对尾涡的演变机制尚不完全了解，且受复杂大气和雷达时空分辨率影响，仅通过尾流速度场呈现出的特征难以准确识别尾涡位置，贝叶斯网络模型可以结合尾涡已知特征构建网络，通过输入数据集学习网络参数，从而解决尾涡定位这一复杂的不确定性推理问题。

2.1.1　贝叶斯原理介绍及数据处理

贝叶斯定理是经典的概率模型之一，基于先验信息和观测数据得到目标变量的后验分布是贝叶斯的核心理论。贝叶斯网络由一个有向无环图和若干条件概率表组成，在本研究中，通过有向无环图定性描述探测点速度特征与涡核相对位置等变量之间的依赖和独立关系，通过条件概率表对这些关系进行定量描述^［32］。以要构建的模型为例，记中间节点 $X$ 的状态输出父节点为 $P a (X)$ ，观测输入子节点为 $C h (X)$ ，贝叶斯网络应满足：存在一个表示中间节点状态的变量集 $V = \{X_{i}\}$ ，其中， $i = 1,2, \dots, n$ ，变量对应节点之间有向边集合为 $E$ ；变量取值可以是离散或连续的，由变量对应的节点之间的有向边构成一个有向无环图 $G = 〈V, E〉$ ，其中， $V$ 为节点集， $E$ 为有向边集；对每个节点 $X_{i}$ 和它父节点集合 $P a (X_{i})$ 都对应一个条件概率分布表 $P [X_{i} | P a (X_{i})]$ ，且满足

P (X_{1}, X_{2}, \dots, X_{n}) = \prod_{i = 1}^{n} P [X_{i} | P a (X_{i})]

。(1)

所提贝叶斯网络模型主要分为3个阶段：准备工作阶段、模型训练阶段和应用阶段，如图1所示。

图 1. 贝叶斯网络实现技术路线图

Fig. 1. Technology roadmap of Bayesian network

下载图片查看所有图片

1）准备工作阶段：收集训练数据集和测试数据集，用于之后的训练和应用阶段，选择能反映涡核位置的特征，构建合适的贝叶斯网络模型，并决定模型中每个节点代表的特征。

2）模型训练阶段：将原始数据处理后生成数据集，使用训练数据集对贝叶斯网络模型进行训练，训练过程包括对每个节点概率分布的估计，以确定每个节点的状态与其他节点的状态之间的联系。

3）应用阶段：将测试数据集输入训练好的贝叶斯网络模型中，通过已经学习好的节点概率分布进行推理，根据给定数据和贝叶斯网络模型计算出每个类别的概率，选择具有最高概率的类别作为结果。

2.1.2　雷达探测数据预处理

应用激光雷达探测尾涡可以得到整个探测范围内所有探测点的回波参数（包括径向速度、谱宽、谱强、载波噪声比等），这通常是一个庞大的数据集，直接应用于尾涡参数反演会导致计算量过大，难以满足实际运行中的效率需求，因此需要通过以下方法对数据集进行处理，提取关键的特征参数后对其离散化，得到适用于贝叶斯网络模型参数训练或推理的数据集，其中，在模型训练阶段使用的规范数据集应包括推理结果（即涡核与最近探测点相对位置范围的类别代码）。从理论上来说，涡心所在的径向距离（相对于激光雷达）上多普勒速度的波动范围较大，而涡核上下探测单元应该具有反向的多普勒速度，因此涡核的仰角应该在正最大速度所在的仰角和负最小速度所在的仰角中间位置。分别计算各距离门上径向速度标准差和不同探测仰角上的径向速度，提取以下数据：

1）左右涡径向速度达到标准差峰值和次大值对应的距离门，左右涡处对应的峰值距离门分别记为 $P_{l e f t}$ 和 $P_{r i g h t}$ ，次大值距离门分别记为 $S_{l e f t}$ 和 $S_{r i g h t}$ ；

2）左右涡径向速度达到极大值和极小值时对应探测仰角的中值，分别记为 $a_{l e f t}$ 和 $a_{r i g h t}$ ；

3）距离门i上的径向速度标准差记为 $σ (i)$ ，计算左右涡处峰值与次大值比值 $R_{l e f t}$ 、 $R_{r i g h t}$ 。

R_{l e f t} = \frac{σ (P_{l e f t})}{σ (S_{l e f t})}

，（2）

R_{r i g h t} = \frac{σ (P_{r i g h t})}{σ (S_{r i g h t})}

。（3）

基于已知的尾涡参数特征，可以构建以下决策树（图2）作为知识库：

图 2. 基于特征参数的涡核范围决策树

Fig. 2. Decision tree of vortex core range based on characteristic parameters

下载图片查看所有图片

图2中，扫描点坐标 $(i, j)$ 中 $i$ 表示距离门数， $j$ 表示探测角度。图2表达了通过对探测点所在距离门和仰角上的速度特征进行判断，可以输出该探测点与涡核的相对位置关系的分类结果。

鉴于决策树分类器过于简单，采用精确门限的硬划分，在表达和推理过程中无法体现出复杂的尾涡特征，识别准确率与其他方法相比没有显著提高。基于其能一定程度上反映涡核特征但不够灵活且容错性较低的特点，将决策树数据库作为先验知识，在此基础上利用贝叶斯网络和基于均方误差的寻优算法来设计涡核识别模型，该方法能在充分利用已有先验知识的前提下提高模型的柔和性和容错性，提高处理不确定知识和推理分类的能力。

2.1.3　网络构建及节点定义

基于实际尾涡速度场特征和第2.1.2节中的决策树，可以提取出关键特征作为节点（观测节点、中间节点、输出节点）并构建贝叶斯网络，网络模型如图3所示，其中，需要输入的A、B、R为表示速度场特征的观测节点，M_A、M_B为中间节点，M_R为输出节点。表示探测角、距离门、径向速度标准差的观测节点A、B、R分别可能存在3、5、5种状态，如表1所示。表示中间节点的M_A、M_B和表示输出节点的M_R分别可能存在3、7、8种状态，分别定义如表2和表3所示。

图 3. 贝叶斯网络模型

Fig. 3. Bayesian network model

下载图片查看所有图片

表 1. 观测节点状态集定义及编号

Table 1. Definition and numbering of state set for observation nodes

Variable	State classification	Numbering
Elevation angle $A$	$a_{l e f t}$	1
	$a_{r i g h t}$	2
	else	3
Range gate $B$	$[1, P_{l e f t} - 1]$	1
	$[P_{l e f t} + 1, P_{r i g h t} - 1]$	2
	$[P_{r i g h t} + 1, \infty)$	3
	$P_{l e f t}$	4
	$P_{r i g h t}$	5
Standard deviation ratio $R$	［2.8，∞）	1
	［2.67，2.8）	2
	［1.42，2.67）	3
	［1，1.42）	4
	［0，1）	5

查看所有表

表 2. 中间节点状态集定义及编号

Table 2. Definition and numbering of state set for intermediate nodes

Variable	State classification	Numbering
M_A	Elevation angle of the left vortex core	1
	Elevation angle of the right vortex core	2
	Else	3
M_B	Left side of the left vortex	1
	Between the two vortices	2
	Right side of the right vortex	3
	Closest point to the left vortex core	4
	Range gate closest to the left vortex core，which not be the closest point	5
	Closest point to the right vortex core	6
	Range gate closest to the right vortex core，which not be the closest point	7

查看所有表

表 3. 输出节点状态集定义及编号

Table 3. Definition and numbering of state set for output nodes

Variable	State classification	Numbering
M_R	Left side of the left vortex	1
	Between the two vortices	2
	Right side of the right vortex	3
	Range gate closest to the vortex core（non-nearest point）	4
	Coinciding with vortex core	5
	Within（0，2.5］m from the vortex core	6
	Within（2.5，5］m from the vortex core	7
	Within（5，7.5］m from the vortex core	8

查看所有表

2.1.4　模型参数学习

常用的激光探测雷达通过发射激光束在一定角度范围内上下来回扫描，扫描到不同角度的时刻也是不同的，因此在仿真数据集时需要考虑到尾涡随时间的演变，德国宇航中心提出的P2P模型^［33-34］和Burnham等^［35］提出的B-H模型分别在尾涡演变过程和速度场反演方面拟合得较好，传统的研究中一般假定左右涡高度和强度相同。但由于实际探测中飞行情况复杂，还会受到大气湍流的影响，故左右涡可能会出现强度或高度不同的情况，因此研究将基于可实现的雷达时空分辨率，以常见的A330-300型号飞机为例，仿真1000组实际中可能出现的不同强度和高度尾涡组合作为数据集进行实验，用以训练网络模型参数，进而得到更准确的尾涡识别模型。相关机型参量及假设大气参数等计算条件如表4所示。

表 4. A330-300机型及环境相关参数

Table 4. Related parameters of A330-300 model and environment

Parameter	Value
Take-off weight $m / k g$	230000
Span $p / m$	60.30
Load coefficient $S$	π/4
Take-off speed $V / (m \cdot s^{- 1})$	85
Air density $ρ / (k g \cdot m^{- 3})$	1.16
Gravity acceleration $g / (m \cdot s^{- 2})$	9.81
Normalized eddy dissipation rate $ε^{*}$	0.07
Normalized Brunt-Väisälä frequency $N^{*}$	0

查看所有表

利用加拿大Norsys公司开发的Netica软件构建贝叶斯网络模型。对1000组数据进行预处理，根据前文方法分别提取其特征，并且离散化，按75%训练集、15%测试集、15%验证集随机分组，训练集数据结合期望最大（EM）算法进行参数优化构建条件概率表（CPT），实现贝叶斯网络的参数学习。将CPT代入后，可以得贝叶斯网络模型参数，如图4所示。

图 4. 训练得到的贝叶斯网络模型参数

Fig. 4. The parameters of Bayesian network model obtained by training

下载图片查看所有图片

之后将总量的15%数据作为测试样本，作为贝叶斯网络特征节点的输入。根据贝叶斯网络学习过程确定的类节点的边缘概率及特征节点的条件概率进行贝叶斯推理，以相对位置节点概率的大小作为判断的依据，即取概率最大处对应的状态作为识别结果。

2.2　基于均方误差对涡核精确位置和环量寻优

基于前文方法计算得到的左右涡核及其最接近的探测点位置参数（Y轴坐标、Z轴坐标、径向距离、探测仰角）可分别表示为：左涡核 $(y_{c l}, z_{c l}, r_{c l}, θ_{c l})$ ；右涡核 $(y_{c r}, z_{c r}, r_{c r}, θ_{c r})$ ；左涡最近探测点 $(y_{d l}, z_{d l}, r_{d l}, θ_{d l})$ ；右涡最近探测点 $(y_{d r}, z_{d r}, r_{d r}, θ_{d r})$ 。涡核和该探测点的距离 $r_{S L}$ ， $r_{S R}$ 范围分别为 $[T_{l 1}, T_{l 2}]$ 和 $[T_{r 1}, T_{r 2}]$ ，以0.5 m为单位（即 $r_{S \cdot} = T_{\cdot 1} ∶ 0.5 ∶ T_{\cdot 2}$ 的一组径向距离）在已知探测仰角上对尾涡对的精确位置和环量进行双重寻优。

首先可以根据径向距离和仰角推算出需要寻优的涡核对位置组合，以 $T_{\cdot 2} - T_{\cdot 1} = 2.5 m$ 为例，该组合应为 $6 \times 6$ 的矩阵，同一行表示左涡位置相同，同一列表示右涡位置相同，则第 $i$ 行第 $j$ 列处的左右涡坐标计算公式为

\{\begin{matrix} r_{S L} (i) = 0.5 \times (i - 1) \\ r_{S R} (j) = 0.5 \times (j - 1) \\ r_{c l} = T_{l 1} + r_{S L} (i) \\ r_{c r} = T_{r 1} + r_{S R} (j) \\ θ_{c \cdot} = θ_{d \cdot} \\ y_{c l} = r_{c l} \times c o s θ_{c l} \\ z_{c l} = r_{c l} \times s i n θ_{c l} \\ y_{c r} = r_{c r} \times c o s θ_{c r} \\ z_{c r} = r_{c r} \times s i n θ_{c r} \end{matrix}

。（4）

因为一般探测时扫描角分辨率较高，可以在估算时认为涡核与探测点在同一扫描角上，即 $θ_{c \cdot} = θ_{d \cdot}$ ，而当角度分辨率 $α$ 不够高（如 $α > 0.5$ ）时， $θ_{c \cdot}$ 可以取 $θ_{d \cdot} \pm 0.5 α$ 。

在已知距离涡核5~15 m内所有探测点（左右涡分别有 $n_{1}$ ， $n_{2}$ 个）的径向速度基础上，对尾涡环量 $(Γ_{L}, Γ_{R})$ 进行寻优。

目标函数的表达式为

\begin{matrix} m i n E_{M S} = \\ \frac{\sum_{1}^{n_{1}} {[V_{s p e c} (a) - V_{r e a l} (a)]}^{2} + \sum_{1}^{n_{2}} {[V_{s p e c} (b) - V_{r e a l} (b)]}^{2}}{n_{1} + n_{2}} \end{matrix}

，（5）

式中： $E_{M S}$ 表示左右涡探测点推测速度和实际探测速度均方误差； $V_{s p e c}$ 和 $V_{r e a l}$ 分别表示探测点在假设环量下的径向速度和实际探测得到的径向速度。当 $E_{M S}$ 最小时，认为寻找到与实际情况最接近的环量推测结果。

V_{s p e c} = f (O_{L}, O_{R}, Γ_{L}, Γ_{R})

，（6）

\{\begin{array}{l} O_{L} = (y_{c l}, z_{c l}) \\ O_{R} = (y_{c r}, z_{c r}) \end{array}

，（7）

f = \{\begin{array}{l} V_{s p e c} = V_{y} \times c o s (a r c t a n \frac{z}{y}) + V_{z} \times s i n (a r c t a n \frac{z}{y}) \\ V_{y} = V_{y L} + V_{y R} = \frac{Γ_{L} (z - z_{c l})}{2 π [{(y - y_{c l})}^{2} + {(z - z_{c l})}^{2} + r_{c L}^{2}]} + \frac{Γ_{R} (z - z_{c r})}{2 π [{(y - y_{c r})}^{2} + {(z - z_{c r})}^{2} + r_{c R}^{2}]} \\ V_{z} = V_{z L} + V_{z R} = \frac{Γ_{L} (y - y_{c l})}{2 π [{(y - y_{c l})}^{2} + {(z - z_{c l})}^{2} + r_{c L}^{2}]} - \frac{Γ_{R} (y - y_{c r})}{2 π [{(y - y_{c r})}^{2} + {(z - z_{c r})}^{2} + r_{c R}^{2}]} \end{array}

，（8）

式中： $O_{L}$ 和 $O_{R}$ 分别表示左右涡核坐标； $Γ_{L}$ 和 $Γ_{R}$ 表示左右涡环量； $(y, z)$ 表示尾涡流场中任意一点坐标； $r_{c L}$ 和 $r_{c R}$ 分别表示左右涡涡核半径。

2.3　背景风估计及剔除

在有风廓线观测条件的机场可以通过观测数据估算大气条件，包括用于尾涡反演的湍流耗散率 $ε$ ^［36-37］。而在没有侧风廓线的情况下，可以在径向探测风速的基础上在大致尾涡位置左右两侧划定侧风取样区域来估计侧风，两涡核间距 $D$ 与机型翼展有关，以 $D$ 约为50 m为例，以 $\frac{D}{2}$ 为单位确定尾涡影响区域，设左右涡核坐标分别为 $(X_{l}, Y_{l})$ 和 $(X_{r}, Y_{r})$ ，则此区域边界顶点为 $(X_{l} - 25, Y_{l} - 25)$ 、 $(X_{l} - 25, Y_{l} + 25)$ 、 $(X_{r} + 25, Y_{r} + 25)$ 、 $(X_{r} + 25, Y_{r} - 25)$ 。同时选取涡核两侧一定范围（本研究中选取距涡核左右25 m以外、水平方向长100 m、竖直方向高50 m的区域）作为侧风提取区域，如图5所示。

图 5. 尾涡影响区域及侧风取值区域

Fig. 5. Wake vortex influence area and crosswind value area

下载图片查看所有图片

分别计算左右侧区域内平均侧风，对尾涡影响区域按探测点到两边侧风区域中心点的距离加权实现侧风剔除，计算过程如下：

1）左侧风区域中心点横坐标 $a = X_{l} - 50$ ，右侧风区域中心点横坐标 $b = X_{r} + 50$ ；

2）将两边侧风区域内探测点处的径向风速转化到水平方向后，分别求其左右两侧的平均风，记为 $w_{l e f}$ 和 $w_{r i g}$ 。

探测点 $(x x, y y)$ 处探测仰角为 $a r c t a n (\frac{y y}{x x})$ ，此处侧风为

w_{(x x, y y)} = \frac{|x x - b|}{|x x - a| + |x x - b|} \times w_{l e f} \times a r c t a n (\frac{y y}{x x}) + \frac{|x x - a|}{|x x - a| + |x x - b|} \times w_{r i g} \times a r c t a n (\frac{y y}{x x})

。（9）

探测得到的信噪比参数 $R_{S N}$ 和谱宽参数 $b_{w}$ 也与涡核位置呈现相关性^［21］，因此在实际应用中，当涡核位置无法准确确定时，可以根据 $R_{S N}$ 、 $b_{w}$ 和雷达跑道相对位置确定尾涡的大致范围后按距离加权剔除侧风，以谱宽为例，以激光雷达为原点建立坐标系，用横坐标表示探测点与雷达的水平距离，纵坐标表示高度，谱宽探测结果如图6所示。可以看出，谱宽最大值所在探测点基本位于跑道中线的延长线上。

图 6. 探测到尾涡时的频谱宽度图

Fig. 6. Spectral bandwith diagram when wake vortex is detected

下载图片查看所有图片

以广州白云机场安装在西跑道的某激光雷达为例，雷达在跑道东侧距跑道中线320 m，距跑道入口1500 m处，按3°下滑角计算，当飞机经过激光雷达探测切面时高度约为80 m，据此划定有效尾涡范围和侧风区域：

1）以激光雷达为原点建立坐标轴，设跑道中线对应的横坐标为 $x_{r c l}$ ，当 $b_{w}$ 最大值探测点 $(x x_{S}, y y_{S})$ 横坐标 $x x_{S}$ 在 $[(x_{r c l} - 150) m, (x_{r c l} + 150) m]$ 范围内，纵坐标 $y y_{S}$ 在120 m以下，则认为在尾涡在跑道范围内；

2）以探测点 $(x x_{S}, y y_{S})$ 为核心划定涡核计算范围，是以 $(x x_{S} - 100, y y_{S} - 50)$ 、 $(x x_{S} - 100, y y_{S} + 50)$ 、 $(x x_{S} + 100, y y_{S} + 50)$ 、 $(x x_{S} + 100, y y_{S} - 50)$ 为顶点的矩形；

3）在此基础上向左右两侧扩展100 m作为侧风计算区域，即左侧是以 $(x x_{S} - 200, y y_{S} - 50)$ 、 $(x x_{S} - 200,$ $y y_{S} + 50)$ 、 $(x x_{S} - 100, y y_{S} + 50)$ 、 $(x x_{S} - 100, y y_{S} - 50)$ 为顶点的矩形，右侧是以 $(x x_{S} + 200, y y_{S} - 50)$ 、 $(x x_{S} + 200, y y_{S} + 50)$ 、 $(x x_{S} + 100, y y_{S} + 50)$ 、 $(x x_{S} + 100, y y_{S} - 50)$ 为顶点的矩形。

3　数据仿真及结果验证

为验证算法性能，设计了一组包括尾涡参数和背景风参数的实验环境，首先根据此实验条件对尾涡流场进行模拟，依据激光雷达参数反演得到径向扫描风速数据，之后用所提算法对尾涡参数进行反演，与最初假设参数对比，以检验其效果。

提取出距离门上探测得到的径向速度，根据第2.1.2节中决策树对所需特征进行离散化，得到各项特征所属类别代码。

1）标准差峰值距离门。计算每个距离门上径向速度的标准差，因为在实际探测数据中可能出现两个以上峰值，通过设定阈值寻找符合条件的峰值及其对应距离门，即在判断是否为峰值时加入了一个条件——峰值必须大于等于阈值（最大峰值的60%），该阈值用于过滤掉不明显的峰值。返回的结果是两个数组，分别是标准差数列中的明显的峰值和它们所对应的距离门。

表 5. 仿真尾涡参数

Table 5. Simulation vortex parameters

Parameter		Left vortex	Right vortex
Position coordinate /m		（450，65）	（510，65）
Circulation		600	600
Background wind	Mean wind /（m·s^-1）	-1
Background wind	Turbulence intensity /（m·s^-1）	1

查看所有表

表 6. BN-MSE寻优推理参数及结果

Table 6. BN-MSE optimization inference parameters and results

Calculated parameter	Left vortex	Right vortex
$E_{M S}$	0.137
BN inference result	7	7
Standard deviation	（4.931，2.217）	（4.185，2.285）
Range gate	（25，26）	（29，30）
Radial distance range	［452.5，455］	［512.5，515］
Estimated angle	8.25	7.25
Estimated coordinates	（448.58，64.80）	（508.68，64.80）
Estimated circulation	593	591

查看所有表

表 7. 所提算法与传统算法推理结果对比

Table 7. Comparison of inference results between the proposed algorithm and traditional algorithms

Parameter	Range method	Gradient method	Proposed algorithm（without side-wind exclusion）	Proposed algorithm（with side-wind exclusion）
Simulation vortex core position（left/right）	（450，65）/（510，65）
Estimated result（left/right）	（445.89，60.68）/506.46，59.94）	（445.62，62.63）/（506.20，62.15）	（448.27，64.75）/（508.07，64.72）	（448.58，64.80）/（508.68，64.80）
Relative deviation of position（left/right）	（6.02 m，7.14 m）	（5.54 m，5.88 m）	（1.75 m，1.95 m）	（1.42 m，1.33 m）
Relative deviation of circulation（left/right）	（-17，-33）	（-18，-21）	（-13，-10）	（-7，-9）

查看所有表

2）梯度筛选。如果步骤1）中得到超过两个峰值，则通过距离门上的速度梯度进一步筛选，因为尾涡在涡核上下发生速度突变，故速度梯度场中左右涡应分别对应正最大-负最小所在的距离门，且左右涡核之间距离一般与飞机翼展相当，在此基础上可进一步通过左右涡相对位置判断选择合理的坐标组合。

3）计算峰值和次大值比值、仰角等各项参数，并按第2.1.3节中的节点定义进行离散化后输入贝叶斯网络，得到与涡核位置最接近的探测点及两者之间的相对位置关系。

4）在贝叶斯推断结果范围（2.5 m）内进行以 $E_{M S}$ 最小为目标函数的参数寻优（径向距离取值间隔为0.5 m，环量取值范围为［1，600］，间隔为1），对左右涡参数进行嵌套循环，代入H-B速度场模型可以得到范围内不同坐标和强度下的探测点径向速度 $V_{s p e c}$ ，结合该点实际探测速度 $V_{r e a l}$ 可以得到本组尾涡参数的均方误差 $E_{M S}$ ，认为 $E_{M S}$ 最小时的参数组合是与实际情况最接近的尾涡坐标和环量值。

相关机型及假设大气参数如表1所示，其他包括涡核坐标、环量及背景风场的计算条件如表5所示。

基于H-B速度场模型模拟以上实验环境下的尾涡流场，进而根据雷达探测原理仿真其探测结果，以激光雷达为原点建立坐标系，用横坐标表示探测点与雷达的水平距离，纵坐标表示高度，仿真速度场探测结果如图7所示，从中可以看出明显的尾涡特征。

图 7. 激光雷达探测风场仿真

Fig. 7. Lidar detection wind field simulation

下载图片查看所有图片

得到实验所需的径向速度场数据并预处理后，将提取出的特征参数代入已构建的BN模型，可以得到距涡核最近的雷达探测点及两者相对位置范围，在本例中，通过Netica软件可以将观测节点结果输入贝叶斯网络，状态固定的节点会变为灰色，最后可以得到输出节点概率分布，取概率最大处对应的状态作为识别结果，计算过程参数如图8和表6所示。

图 8. 仿真样本贝叶斯网络推理结果。（a）左涡最近探测点；（b）右涡最近探测点

Fig. 8. Bayesian network inference results of simulation sample. (a) Left vortex nearest detection point; (b) right vortex nearest detection point

下载图片查看所有图片

由图8可以看出，本例中的输出结果均为编号为7的状态，即最近的探测点与涡核距离在（2.5 m，5 m］范围内。基于第2.2节中的方法在该范围内反演寻优，可以得到尾涡参数和对应的均方误差，具体如表6所示。

传统方法一般基于涡核附近速度发生突变这一特性来推测涡核位置，如梯度法和极差法。这类方法计算较为简便，能大致推测出涡核所在范围，但缺点是在探测时空分辨率有限或背景风场复杂的情况下难以得到精确的涡核位置。将前文介绍的仿真算例分别代入传统算法（极差法、梯度法）和所提算法，得到的尾涡坐标、环量推理结果与仿真参数偏差如表7所示。

4　实际算例

4.1　实验环境

尾涡探测实验地点为广州白云国际机场，该机场西跑道旁共装有两部脉冲相干多普勒激光雷达，该型近年来多次用于航空领域探测和安全保障^［38-40］，雷达选址和参数如表8所示，其中，PPI、DBS、RHI分别为平面位置显示、多普勒波束锐化、距离高度显示的缩写。

表 8. 广州白云国际机场雷达扫描参数

Table 8. Lidar scanning parameters at Guangzhou Baiyun international airport

Parameter	Lidar 1	Lidar 2
Lidar model	Wind3D 6000	Wind3D 6000
Location	North end of west runway（19 end）	South of the west runway（01 end）
Scanning strategy	PPI，DBS，RHI	RHI
Range resolution /m	150	15
Azimuth angle /（°）	RHI：16；	286
Azimuth angle /（°）	PPI：0‒359；	286
Elevation angle /（°）	RHI：0‒180	0‒180
	PPI：3；
	DBS：71.38

查看所有表

其中，激光雷达2的扫描策略和距离分辨率更适用于尾涡定位和参数反演，故采用此雷达的探测数据作为实验算例，该激光雷达被放置在西跑道南端（01端）以南的地面上，主要用来扫描机场西跑道上飞机起降阶段产生的尾流。激光雷达在垂直于机场跑道的平面上循环扫描，相应的三维示意图如图9所示。其中，用五角星表示的激光雷达与雷达中线的距离约为320 m，在距离跑道入口约1500 m处，按3°下滑角计算，飞机到达雷达扫描截面的高度约为80 m。

图 9. 激光雷达与跑道相对位置示意图

Fig. 9. The relative position diagram of lidar and runway

下载图片查看所有图片

激光雷达输出扫描平面上各探测单元的多普勒速度，实验中采用的激光雷达为中国海洋大学与青岛镭测创芯科技有限公司联合研制的Wind3D 6000三维扫描型雷达。该激光雷达的主要技术指标如表9所示^［21］。

表 9. 激光雷达技术指标

Table 9. Technical specifications of Lidar

Parameter	Specification
Wavelength /nm	1550
Pulse repetition rate /kHz	10
Pulse energy /μJ	160
Pulse width /ns	100‒200
Power consumption /W	$<$ 300
Radial velocity measurement range /（m $∙$ s^-1）	-37.5‒37.5
Velocity measurement uncertainty /（m $∙$ s^-1）	$\leq$ 0.1
Measurement range /m	40‒6000
Range resolution /m	15‒30

查看所有表

除此之外，RHI可以根据实际情况设定，本研究中，RHI为161°~179°，一个RHI周期为8 s。

4.2　探测实验结果验证

以广州白云机场2022年11月1日的探测数据为例，与当日航班信息匹配后，本小节选取了8组不同机型的尾涡扫描数据来展示所提算法在实测数据上的应用情况，验证其是否能基于不同机型有效反演出尾涡参数，表10介绍了选取的8个实例应用所提算法得到的反演结果，并与传统算法反演结果的均方误差进行了比较。

表 10. 所提算法与传统算法反演结果均方误差比较

Table 10. Mean squared error comparison of inversion results between the proposed model and traditional algorithms

Case number	Aircraft type	Estimated value		$E_{M S}$	Improvement compared to traditional method /%
Case number	Aircraft type	Position coordinate	circulation	$E_{M S}$	Gradient method	Range method
1	B763	（-350.20，64.02）（-308.98，62.05）	（157，160）	0.438	48.83	36.34
2	A359	（-359.79，48.13）（-322.11，61.85）	（144，129）	0.249	59.58	80.82
3	B744	（-355.85，47.42）（-300.18，48.03）	（143，110）	0.207	81.19	17.20
4	B77L	（-356.74，55.29）（-323.36，48.67）	（161，145）	0.316	65.95	74.23
5	B77L	（-376.17，53.81）（-327.44，41.02）	（125，138）	0.428	63.52	49.13
6	A346	（-459.39，24.12）（-401.66，16.28）	（181，137）	0.545	48.33	34.16
7	A333	（-376.59，69.80）（-320.53，69.59）	（176，135）	0.256	37.11	53.24
8	B752	（-361.41，33.91）（-334.93，37.33）	（111，113）	0.212	61.03	55.04

查看所有表

以上8个实例的雷达探测速度场如图10所示，图中用黑色实心圆标注出了反演的涡核位置，从图中可以看出，反演结果基本与已知的尾涡特征吻合。

图 10. 8组算例涡核反演结果。（a）算例1；（b）算例2；（c）算例3；（d）算例4；（e）算例5；（f）算例6；（g）算例7；（h）算例8

Fig. 10. Eight sets of vortex core inversion results. (a) Example 1; (b) example 2; (c) example 3; (d) example 4; (e) example 5; (f) example 6; (g) example 7; (h) example 8

下载图片查看所有图片

从图10可以看出，所提算法可以对涡核在高度10~80 m范围内的尾涡位置和环量进行反演，并得到均方误差较小的结果。

5　结论

在给出研究中常用的多普勒激光雷达探测原理和基于探测结果的尾涡特征参数计算方法的基础上，对飞机尾涡的探测和识别进行了研究，建立了基于BN和均方误差寻优的尾涡参数反演模型，通过对有背景风条件下的尾涡速度场雷达探测结果进行模拟并代入模型，验证了所提算法的有效性。与传统算法中基于极差和梯度特点推测涡核位置相比，所提算法能够得到更精确的结果。

通过雷达探测结果可以计算得到距离门上的极差、梯度和标准差等参数，基于已知的尾涡特征可以建立决策树模型判断雷达探测点与涡核的相对位置，考虑到参数缺失和探测结果存在误差的情况，可以在决策树的基础上对参数和结果进行离散化以构建BN模型，模拟1000组尾涡探测数据作为模型的训练集，在输入探测参数后模型可以输出涡核坐标的范围，径向距离上分辨率为2.5 m。

与传统方法能得到的15 m范围相比，上述步骤已经将涡核位置缩小至径向2.5 m范围内，基于最小均方误差在此范围内进行涡核坐标和环量的循环寻优，可以极大减少计算量，以较高效率得到与模拟情况较为接近的反演结果，在代入的算例验证中，径向距离估算结果偏差不超过2 m，且基于估计的涡核位置反演得到的环量偏差在5%以内。

在与航班时刻对比后选取8组不同机型的尾涡探测数据，预处理后应用所提算法与传统算法进行对比，反演速度场的均方误差显著降低，相比梯度法和极差法分别平均降低了63.89%和52.15%。

本研究目前建立在15 m径向分辨率激光雷达探测数据的基础上，之后的研究可以讨论对不同分辨率雷达数据的适用情况。且在实例验证阶段发现，所提算法在尾涡消散后段的反演性能有所下降，这是地面效应和尾涡强度减小导致的，后续希望加入对地面效应的研究考虑，并考虑尾涡演变的时序性。

基于激光雷达探测结果的尾涡研究可以深化对飞机涡流的认识，进一步可以在保证安全的前提下缩短间隔，提升跑道和空域容量，这在民航业的快速发展，空中飞行流量增长的背景下具有重大意义。

参考文献

[1] HolzäpfelF, KauertzS, KonopkaJ, et al. Aircraft wake vortex state-of-the-art & research needs[R]. Oberpfaffenhofen: WakeNet3-Europe, 2015.

[2] Barbaresco F. Wake vortex detection, prediction and decision support tools: new challenge for airports to increase capacity and safety[J]. Revue REE, 2013, 3: 15-25.

[3] MutuelL H, BarbarescoF, JugeP, et al. ATM decision support tool for wake vortex hazard management combining sensors and modeling[C]∥Proceedings of the 6th AIAA Atmospheric and Space Environments Conference, June 13-17, 2014, Atlanta, GA. Reston, Virginia: AIAA, 2014: 2332.

[4] 潘卫军, 罗玉明, 韩帅, 等. ARJ21飞机遭遇不同前机尾流的响应和安全性研究[J]. 飞行力学, 2022, 40(1): 13-18, 25.

Pan W J, Luo Y M, Han S, et al. Research on the response and safety of ARJ21 aircraft encounters wake vortex of different front aircrafts[J]. Flight Dynamics, 2022, 40(1): 13-18, 25.

[5] ChengJ, HoffA, TittsworthJ, et al. The development of wake turbulence re-categorization in the United States (invited)[C]∥Proceedings of the 8th AIAA Atmospheric and Space Environments Conference, June 13-17, 2016, Washington, D.C. Virginia: AIAA, 2016: 3434.

[6] 周艳宗, 王冲, 刘燕平, 等. 相干测风激光雷达研究进展和应用[J]. 激光与光电子学进展, 2019, 56(2): 020001.

Zhou Y Z, Wang C, Liu Y P, et al. Research progress and application of coherent wind lidar[J]. Laser & Optoelectronics Progress, 2019, 56(2): 020001.

[7] Holzäpfel F, Gerz T, Köpp F, et al. Strategies for circulation evaluation of aircraft wake vortices measured by lidar[J]. Journal of Atmospheric and Oceanic Technology, 2003, 20(8): 1183-1195.

[8] Köpp F, Rahm S, Smalikho I. Characterization of aircraft wake vortices by 2-μm pulsed Doppler lidar[J]. Journal of Atmospheric and Oceanic Technology, 2004, 21(2): 194-206.

[9] Köpp F, Rahm S, Smalikho I, et al. Comparison of wake-vortex parameters measured by pulsed and continuous-wave lidars[J]. Journal of Aircraft, 2005, 42(4): 916-923.

[10] KomatsubaraT, KakuN. Wake vortex detection of departure aircraft(WSANE2006)[R]. Online: Technical Report of Ieice Sane, 2006.

[11] BarbarescoF, JeantetA, MeierU. Wake vortex detection & monitoring by X-band Doppler radar: Paris Orly radar campaign[C]∥IET International Conference on Radar Systems 2007, October 15-18, 2007, Edinburgh, UK. London: IEE, 2007: 1-5.

[12] Rahm S, Smalikho I. Aircraft wake vortex measurement with airborne coherent Doppler lidar[J]. Journal of Aircraft, 2008, 45(4): 1148-1155.

[13] Smalikho I N, Rahm S. Lidar investigations of the effects of wind and atmospheric turbulence on an aircraft wake vortex[J]. Atmospheric and Oceanic Optics, 2010, 23(2): 137-146.

[14] Smalikho I N, Banakh V A, Holzäpfel F, et al. Method of radial velocities for the estimation of aircraft wake vortex parameters from data measured by coherent Doppler lidar[J]. Optics Express, 2015, 23(19): A1194-A1207.

[15] Darracq D, Corjon A, Ducros F, et al. Simulation of wake vortex detection with airborne Doppler lidar[J]. Journal of Aircraft, 2000, 37(6): 984-993.

[16] EhlersJ, FezansN. Airborne Doppler LiDAR sensor parameter analysis for wake vortex impact alleviation purposes[M]∥Bordeneuve-Guibé J, Drouin A, Roos C. Advances in aerospace guidance, navigation and control. Cham: Springer, 2015: 433-453.

[17] HallermeyerA, Dolfi-BouteyreA, VallaM, et al. Development and assessment of a Wake Vortex characterization algorithm based on a hybrid LIDAR signal processing[C]∥Proceedings of the 8th AIAA Atmospheric and Space Environments Conference, June 13-17, 2016, Washington, D.C. Reston, Virginia: AIAA, 2016: 3272.

[18] Smalikho I N, Banakh V A, Falits A V. Measurements of aircraft wake vortex parameters by a Stream Line Doppler lidar[J]. Atmospheric and Oceanic Optics, 2017, 30(6): 588-595.

[19] Smalikho I N. Taking into account the ground effect on aircraft wake vortices when estimating their circulation from lidar measurements[J]. Atmospheric and Oceanic Optics, 2019, 32(6): 686-700.

[20] Smalikho I N, Banakh V A, Falits A V, et al. Experimental study of aircraft wake vortices on the airfield of tolmachevo airport in 2018[J]. Atmospheric and Oceanic Optics, 2020, 33(2): 124-133.

[21] 王筱晔, 吴松华, 刘晓英, 等. 基于相干多普勒激光雷达的飞机尾涡观测[J]. 光学学报, 2021, 41(9): 0901001.

Wang X Y, Wu S H, Liu X Y, et al. Observation of aircraft wake vortex based on coherent Doppler lidar[J]. Acta Optica Sinica, 2021, 41(9): 0901001.

[22] 赵丽雅, 谷润平, 魏志强. 基于激光雷达回波的动态尾涡特征参数计算[J]. 武汉科技大学学报, 2018, 41(5): 388-394.

Zhao L Y, Gu R P, Wei Z Q. Calculation of characteristic parameters of dynarnic wake vortex based on lidar echo[J]. Journal of Wuhan University of Science and Technology, 2018, 41(5): 388-394.

[23] 谷润平, 赵丽雅, 魏志强. 飞机尾涡流场特征参数估算方法研究[J]. 航空计算技术, 2017, 47(6): 14-17, 23.

Gu R P, Zhao L Y, Wei Z Q. Study on estimation method of characteristic parameters of aircraft wake vortex[J]. Aeronautical Computing Technique, 2017, 47(6): 14-17, 23.

[24] Hon K, Chan P. Aircraft wake vortex observations in Hong Kong[J]. Journal of Radars, 2017, 6(6): 709-718.

[25] Shen C, Li J B, Zhang F L, et al. Two-step locating method for aircraft wake vortices based on Gabor filter and velocity range distribution[J]. IET Radar, 2020, 14(12): 1958-1967.

[26] Gao H, Li J B, Chan P W, et al. Parameter-retrieval of dry-air wake vortices with a scanning Doppler Lidar[J]. Optics Express, 2018, 26(13): 16377-16392.

[27] Li J B, Shen C, Gao H, et al. Path integration (PI) method for the parameter-retrieval of aircraft wake vortex by lidar[J]. Optics Express, 2020, 28(3): 4286-4306.

[28] 沈淳, 高航, 王雪松, 等. 基于激光雷达探测的飞机尾流特征参数反演系统[J]. 雷达学报, 2020, 9(6): 1032-1044.

Shen C, Gao H, Wang X S, et al. Aircraft wake vortex parameter-retrieval system based on lidar[J]. Journal of Radars, 2020, 9(6): 1032-1044.

[29] 高航, 周洁, 胡健, 等. 基于激光雷达探测的三维风场变分反演算法[J]. 光学学报, 2021, 41(20): 2028002.

Gao H, Zhou J, Hu J, et al. Variational retrieval algorithm for three-dimensional wind field based on lidar detection[J]. Acta Optica Sinica, 2021, 41(20): 2028002.

[30] 沈淳, 李健兵, 高航, 等. 基于数据同化的飞机尾流行为预测[J]. 雷达学报, 2021, 10(4): 632-645.

Shen C, Li J B, Gao H, et al. Aircraft wake vortex behavior prediction based on data assimilation[J]. Journal of Radars, 2021, 10(4): 632-645.

[31] 张连文, 郭海鹏. 贝叶斯网引论[M]. 北京: 科学出版社, 2006.

ZhangL W, GuoH P. Introduction to Bayesian networks[M]. Beijing: Science Press, 2006.

[32] 李硕豪, 张军. 贝叶斯网络结构学习综述[J]. 计算机应用研究, 2015, 32(3): 641-646.

Li S H, Zhang J. Review of Bayesian networks structure learning[J]. Application Research of Computers, 2015, 32(3): 641-646.

[33] Holzäpfel F. Probabilistic two-phase wake vortex decay and transport model[J]. Journal of Aircraft, 2003, 40(2): 323-331.

[34] Holzäpfel F. Probabilistic two-phase aircraft wake-vortex model: further development and assessment[J]. Journal of Aircraft, 2006, 43(3): 700-708.

[35] BurnhamD C, HallockJ N. Chicago monostatic acoustic vortex sensing system[R]. New York: U.S. Department of Transportation, 1982.

[36] 靳翔, 宋小全, 刘佳鑫, 等. 基于多普勒激光雷达的边界层内湍流参数估算[J]. 中国激光, 2021, 48(11): 1110001.

Jin X, Song X Q, Liu J X, et al. Estimation of turbulence parameters in atmospheric boundary layer based on Doppler lidar[J]. Chinese Journal of Lasers, 2021, 48(11): 1110001.

[37] 陈晓敏, 张洪玮, 孙康闻, 等. 基于相干多普勒激光雷达的斜程湍流参数反演方法研究[J]. 大气与环境光学学报, 2023, 18(1): 1-13.

Chen X M, Zhang H W, Sun K W, et al. Inversion methods of slant turbulence parameters based on coherent Doppler lidar[J]. Journal of Atmospheric and Environmental Optics, 2023, 18(1): 1-13.

[38] Wu S H, Zhai X C, Liu B Y. Aircraft wake vortex and turbulence measurement under near-ground effect using coherent Doppler lidar[J]. Optics Express, 2019, 27(2): 1142-1163.

[39] 张洪玮, 吴松华, 尹嘉萍, 等. 基于短距相干测风激光雷达的机场低空风切变观测[J]. 红外与毫米波学报, 2018, 37(4): 468-476.

Zhang H W, Wu S H, Yin J P, et al. Airport low-level wind shear observation based on short-range CDL[J]. Journal of Infrared and Millimeter Waves, 2018, 37(4): 468-476.

[40] Zhang H W, Wu S H, Wang Q C, et al. Airport low-level wind shear lidar observation at Beijing Capital International Airport[J]. Infrared Physics & Technology, 2019, 96: 113-122.

2.2　基于均方误差对涡核精确位置和环量寻优

谷润平, 鹿彤, 魏志强. 激光雷达探测中基于贝叶斯网络的飞机尾流反演[J]. 激光与光电子学进展, 2024, 61(4): 0428006. Runping Gu, Tong Lu, Zhiqiang Wei. Aircraft Wake Inversion Based on Bayesian Network in Lidar Detection[J]. Laser & Optoelectronics Progress, 2024, 61(4): 0428006.