沈 淳 高 航 王雪松 李健兵
(国防科技大学电子科学学院 长沙 410073)
飞机尾流是飞机飞行时机翼上下表面压力差而在其后方形成的一对反向旋转的湍流漩涡,具有空间尺度大、持续时间长、旋转强烈等特点[1]。作为伴随飞机出现的大气湍流中的一个新类别,飞机尾流对航空安全有重大威胁,当后机进入前机的尾流时,可能发生颠簸、横滚,乃至失去控制。据美国国家交通安全委员会的统计,1993年—2000年间美国境内约1/3的航空事故与尾流相关[2];2001年,一架美联航飞机在纽约肯尼迪机场起飞后因遭遇前机尾流而坠毁,导致265人死亡,这是美国自“911”后最大的航空事故。在机场和航空母舰的进港/离港航线上,由于飞机起降密集,飞机尾流的影响更大,是制约机场吞吐容量、航母舰载机安全起降的主要因素之一。
为了避免飞机尾流的危害,国际民航组织(International Civil Aviation Organization,ICAO)于上世纪70年代起陆续制定了飞机之间的间隔规则以规避尾流的危害[3],但这些规则没有考虑天气条件、机场跑道布局等差异,总体比较保守,过大的间隔很大程度上限制了机场的起降容量,亦容易导致飞机晚点等现象。
随着经济发展全球化进程的不断提升,空中交通在近几十年来得到了快速发展(大约每15年增长1倍),这使得机场和空域运行资源日趋紧张,亟需发展新的尾流间隔标准,特别是基于雷达探测的动态尾流间隔。欧洲统一天空(Single European Sky Air traffic management Research,SESAR)[4]和美国下一代航空运输系统(Next Generation air transportation system,NextGen)[5]这两大空管计划以及国际民航组织的航空系统模块升级(Aviation System Block Upgrade,ASBU)计划都将飞机尾流间隔优化作为其中的重要内容。中国民航亦在逐步引入飞机尾流间隔优化的概念,并于2019年12月5日起依托中南空管局在广州白云机场与深圳宝安机场开展航空器尾流间隔重新分类(RE-CATegorization of aircrafts,RECAT)的试验运行工作,验证其对中国民航的可行性。
飞机尾流间隔优化的前提是对飞机尾流进行有效的探测和监视,在晴空条件下,飞机尾流内部包括浮尘、气溶胶、花粉等微粒,其尺寸非常小,与激光波长量级相当。同时激光雷达具有角分辨率高等优点,所以是当前探测晴空飞机尾流广泛使用的传感器。激光雷达对晴空飞机尾流探测的主要研究团队包括德国宇航中心(Deutsches zentrum für Luft-und Raumfahrt,DLR)[6]、俄罗斯科学院[7]、Leosphere公司[8]、日本宇航探索局(Japan Aerospace eXploration Agency,JAXA)[9]等;国防科技大学王雪松教授课题组自2005年开始,对飞机尾流的动力学特性和电磁散射特性进行了大量研究积累[10—21],在此基础上于2019年8月开始在长沙黄花机场开展实地探测,通过实验发现激光雷达可以对晴空条件下的飞机尾流在近程上实现很好的探测。
基于激光雷达探测构建飞机尾流特征参数反演系统,可得到飞机尾流的演化轨迹和强度信息,为后续利用尾流特征参数进行飞机尾流的短时行为预测和危害评估提供数据支撑,为制定适合中国各地机场地形条件和实际气象环境变化的安全起降间隔标准提供系统仿真和实测验证的技术支持。对于我国突破国际垄断,为军、民航以及航母舰载机快速起降积累核心技术,具有重大的应用价值和科学意义。
基于激光雷达探测的飞机尾流特征参数反演系统以激光雷达探测器和高性能计算机为平台,构建软硬件结合的数据处理和分析系统,实现雷达实测数据采集、回波数据杂波滤除和修正、飞机尾流定位和强度估计的完整处理流程,同时实现了飞机尾流回波的模拟仿真,可进行算法分析和验证。系统主要由4个模块构成:激光雷达回波仿真模块、雷达实测数据接收与显示模块、飞机尾流涡心定位模块和飞机尾流速度环量估计模块。系统流程如图1所示。
系统数据分析和处理流程可概述为:首先通过激光雷达回波模块产生回波信号,或者通过雷达实测数据接收模块接收回波信号;然后通过飞机尾流涡心定位模块对飞机尾流进行定位得到飞机尾流的涡心位置、涡心间距和涡心轨迹等特征参数;最后,对定位得到的飞机尾流进行速度环量的估计,得到飞机尾流的强度信息。
图1 飞机尾流特征参数反演系统流程图Fig.1 Flow chart of aircraft wake vortex parameter-retrieval system
3.1.1 复杂大气环境模拟模型
飞机尾流所处的背景风场往往包含复杂的湍流、风切变等情况,因此为模拟激光雷达探测飞机尾流回波模型更加贴近实际情况,需模拟出复杂风场背景。
其中,平均风是指风速的基准值,通常用来表明特定时间的风速平均值。风切变是指空间中任意两点间风向或风速的突然变化,低空情况(距离地面高度小于600 m)下可划分为水平风的水平切变、水平风的垂直切变和垂直风的切变[22]。大气湍流是介质的不规则运动,其参数极其复杂且无规律变化,可采用基于冯·卡门(Von Kármán)频谱的拟小波方法来模拟随机湍流场。拟小波方法采用成群的小漩涡结构(拟小波)来模拟大气湍流,这些拟小波在空间中有随机的涡心位置和结构,以此代表湍流的随机性[23,24]。
3.1.2 飞机尾流速度模型
飞机尾流稳定段左右漩涡消散速度基本一致,尾流速度可以简化为左右两个漩涡产生的速度场叠加的合速度,如图2所示。
通过已有的尾流漩涡速度模型可以得到单一漩涡在某一位置处的切向速度受两个漩涡共同影响的探测单元的合速度为
常见的尾流速度模型有Rankine速度模型、Hallock-Burnham速度模型和Lamb-Oseen 速度模型等[25]。Hallock-Burnham 速度模型与实测数据吻合度较高[26],而Rankine模型和Lamb-Oseen模型存在对环量的低估[27],故采用Hallock-Burnham速度模型来计算尾流漩涡的切向速度V α。模型的表达式为
图2 飞机尾流速度矢量示意图Fig.2 Velocity components of aircraft wake vortices
其中,Γ和rc分别为漩涡的环量和漩涡核半径,r是(xi,yi)到漩涡涡心的距离。对于飞机尾流来说,rc一般取两个漩涡涡心间隔的0.052倍。
3.1.3 飞机尾流演化模型
飞机尾流的演化情况和生命长度与大气温度、湍流能量耗散率、风切变以及尾流的初始强度、初始位置等都有关系[28]。系统采用Frank Holzäpfel[29]提出的两阶段(Probabilistic Two-Phase,P2P)尾流耗散模型来模拟飞机尾流的强度和位置在复杂大气环境中的演化过程。这个模型充分考虑了近地面效应、侧风、湍流和稳定分层结构等情况,模型可以输出在一定置信概率区间内尾流的环量,涡心位置随观测时间的变化。P2P方法已被广泛应用于对二维情况下飞机尾流演化情况的模拟。
3.1.4 激光雷达回波模型
激光雷达接收的是探测单元内大量随机分布的气溶胶的所有后向散射信号,忽略粒子之间的互耦,可认为每个探测单元的后向散射信号为所有气溶胶后向散射信号的相干叠加[30]。根据Zrnic的多普勒激光雷达大气信号模型[31],激光雷达接收到的探测单元内第k个气溶胶的回波信号是后向散射信号和噪声的叠加:
其中,Ts=1/B为采样时间间隔,B为激光雷达采样率,λ为波长,j 为虚数单位,nk为随机噪声信号,vk为探测单元内第k个浮尘粒子的速度(由尾流速度和背景风速度共同决定),ρk为第k个粒子的激光散射强度,sk为激光雷达的发射信号。
激光雷达信号处理把相邻N个采样回波信号补零到M点后再作M点的快速傅里叶变换(Fast Fourier Transform,FFT)得到多普勒频谱。由于空气中浮尘等的散射较弱,激光雷达单次探测的信噪比SNR有限。为提升回波的信噪比,可以对多个脉冲重复周期(PRT)之间的多普勒频谱进行累积,以提高激光雷达的探测能力[7],累积后的频谱为
进一步用矩估计方法可从频谱SD的1阶矩和2阶矩得到回波信号的多普勒速度和多普勒谱宽等信息。
激光雷达探测飞机尾流一般采用侧视方式,对飞机飞行经过的垂直平面进行机械扫描得到尾流径向速度分布图。具体探测场景如图3所示,激光雷达放置于飞机跑道一侧,雷达波束垂直于飞机跑道,根据飞机到达雷达扫描切面位置时的高度设置一定的机械扫描角度范围(例如0~30°),上下对飞机尾流进行切片扫描,形成连续的RHI(Range Height Indicator)速度分布图(如图4所示)。实测数据接收模块的基本功能是:根据不同型号激光雷达的数据格式,设置相应的接收模式,对数据进行分析判读,完整无损反演成逐帧RHI格式的数据,同时根据需要接收相应的气象条件参数和飞机飞行数据,以便后续的飞机尾流特征参数反演。
从激光雷达探测得到的RHI速度分布图中可以发现,飞机尾流稳定段左右漩涡的速度场与背景风场速度有明显区别,尾流涡心上下单元具有相反的多普勒速度(如图5所示),涡心位置正好处于正负多普勒速度之间,可以根据此差别来判断是否存在飞机尾流并进行定位。
图3 激光雷达探测飞机尾流示意图Fig.3 Geometry setup of wake vortex Lidar detection
图4 激光雷达探测多普勒速度RHI分布示意图Fig.4 RHI distribution of Doppler velocity by Lidar
图5 飞机尾流左右涡心回波数据示意图Fig.5 Velocity distribution of the left and right wake vortices
对涡心定位的方法有许多种,如基于Gabor滤波器的涡心定位法、多普勒速度极差法、多普勒速度绝对值之和法、多普勒速度平方和法等,可单独使用其中的一种方法进行定位,也可以结合多种方法进行定位,例如首先采用Gabor滤波进行粗定位,再在粗定位划定的范围内利用多普勒速度极差等方法进行精细定位。
3.3.1 Gabor滤波定位模块
二维Gabor滤波器由Daugman[32]首次提出,是一个经过复指数函数调制的高斯核函数,可以获得较高的时域分辨率和频域分辨率,与人类视觉系统中细胞的视觉刺激响应相似,可以应用于图像的特征提取。
飞机尾涡在垂直方向(y方向)上具有一正一负的速度起伏,故只需在垂直方向上对RHI速度分布进行处理,用于飞机尾流涡心定位的Gabor滤波核函数为
其中,x和y分别代表水平距离和垂直距离,j为虚数单位,σ表示Gabor函数的高斯因子沿x和y方向的标准差,一般取值为0.5倍Gabor滤波器尺寸,滤波器尺寸大小约为15 m(稍大于涡心上下极大和极小多普勒速度位置之间的距离);µ是与σ和Gabor函数波长λ有关的系数,Gabor函数波长λ=µσ,µ一般取值为0.5σ。
图6 Gabor滤波后幅度二维分布图Fig.6 Two-dimensional amplitude distribution after Gabor filter
Gabor滤波涡心定位具体方法如下:首先利用Gabor滤波器对RHI回波数据图进行滤波,可以得到处理后的Gabor滤波值二维分布图(如图6所示),其中分布图中存在的极大和极小值可能分别对应尾流左右涡心;然后对极值进行判断,因存在噪声和湍流等因素的影响,需移除RHI回波数据图中的边界以消除边界突变,再将极大和极小值进行一一配对,利用距离判断法对极值配对进行筛选,得到满足条件的极值配对,再选择其中滤波值绝对值相乘最大的极值对作为尾流的左右涡心位置。
其中距离判断法是指在极值点距离地面垂直距离>1.5倍飞机翼展时,极值对需满足极大值和极小值的水平方向距离≤1.5倍的飞机翼展,垂直方向距离≤1倍的飞机翼展;在极值点距离地面垂直距离≤1.5倍飞机翼展(考虑近地面效应可能使左右涡距离变得较大)时,极值对需满足极大值和极小值的水平方向距离≤2倍的飞机翼展,垂直方向距离≤1倍的飞机翼展。
3.3.2 多普勒速度信息定位模块
理论上来说,漩涡涡心所在的径向距离(相对于激光雷达)上多普勒速度的波动范围较大,可以利用这个特性来对涡心位置进行定位。
(1) 径向距离确定
(a) 多普勒速度极差方法:该方法利用距离激光雷达为Rk的所有探测单元的最大多普勒速度和最小多普勒速度之差来表征一个径向距离上多普勒速度的波动大小,称为“极差”
其中,φ为激光雷达扫描波束的仰角。
得到的极差在不同径向距离(相对于激光雷达)上的值如图7,可以发现在不同距离上极差的两个峰值,这两个峰值分别对应着左右漩涡涡心,其中Rc1,Rc2分别表示左右漩涡到激光雷达的径向距离。
(b) 多普勒速度绝对值之和方法:该方法利用距离激光雷达为Rk的所有探测单元的多普勒速度绝对值之和来表征波动大小,与多普勒速度极差方法相似,可以得到在不同距离上多普勒速度之和的两个峰值,这两个峰值分别对应着左右漩涡涡心
图7 多普勒速度极差随径向距离的变化Fig.7 Variation of Doppler velocity range along radial distance
(c) 多普勒速度平方和方法:该方法利用距离激光雷达为Rk的所有探测单元的多普勒速度平方之和来表征波动大小,同样可以得到在不同距离上多普勒速度平方之和的两个峰值,这两个峰值分别对应着左右漩涡涡心
(2) 俯仰角确定
漩涡核上下探测单元具有反向的多普勒速度,因此上述3种多普勒速度信息定位方法中,漩涡核的仰角应该在正最大速度所在的仰角和负最小速度所在的仰角中间位置[7],如图8所示,即通过下式分别计算得到两个漩涡涡心的仰角
图8 定位漩涡涡心仰角的说明图Fig.8 Determination of wake vortex cores’ elevation angles
其中,i=1表示左漩涡,i=2表示右漩涡。
飞机尾流速度环量代表尾流漩涡强度,是后续飞机遭遇尾流危害评估的关键指标。系统中实现了3种速度环量估计算法,包括基于路径积分、基于最优化方法和基于多普勒速度极差的速度环量估计,可有效反演出尾流环量特征参数。
3.4.1 基于路径积分的速度环量估计
基于路径积分的速度环量估计[33]是根据飞机尾流内部空气浮尘和胶质体的动力学特性,在前期涡心定位获取的左右涡心位置上下提取多个LOS(Line Of Sight)上的多普勒分辨单元,进行闭合积分得到尾流速度环量强度。
具体方法如下:首先在激光雷达获取的RHI二维多普勒速度分布图中去除背景风分量;然后根据前期涡心定位获得的涡心位置,在涡心位置的上下提取多个不同长度的分辨单元(如图9(a)所示),采用沿着点 An到点 Bn,再到点 Cn,最后回到 An的闭合路径进行积分的方法进行尾流速度环量的反演(如图9(b)所示);最后去除激光雷达RHI机械扫描引起的尾流压缩/扩展效应,采用最小二乘法,利用计算得到的速度环量对不同时刻的涡心位置不断进行迭代修正,直到得到稳定的速度环量,环量计算公式为
3.4.2 基于最优化方法的速度环量估计
最优化方法[34]基于大气中悬浮粒子的弱惯性,建立了表征多普勒速度和特征参数之间关系的目标方程;然后利用最优化方法求解该目标方程,得到特征参数。
第i个探测单元的径向速度(与多普勒速度大小相等,速度相反表征该探测单元处流体速度(Vflow,(i))在激光雷达波束上的投影
图9 尾流环量积分示意图Fig.9 Path integration of wake vortex circulation
其中,Vwind,(i)和Vwake,(i)分别为该探测单元处的背景风速度和尾流速度,X|r=K·X为速度矢量X在激光雷达波束上的投影,K为仰角为α的波束的方向矢量
根据式(12),得到由背景风速度、尾流速度、多普勒速度构成的飞机尾流特征参数反演的目标方程
其中,包含6个未知数(OL,OR中各含有两个未知坐标),因此当采用M ≥6个探测单元的数据(i ≤M)时,可以通过最优化方法(如非线性最小二乘法)求解目标方程,得到待求的目标参数[ΓL,ΓR,OL,OR]
3.4.3 基于多普勒速度极差的速度环量估计
极差方法[35]利用涡心定位得到的涡心位置,通过涡心上下的多普勒速度之差求出左右尾流漩涡的环量Γ1,Γ2。
图10中圆圈为尾流漩涡的抽象示意,虚线为激光雷达发射的波束,点划线给出了探测单元到漩涡涡心位置的距离。
图10中详细画出了左涡心上下两个探测单元(Rci,φk,i)(i=1,2,k=±1)的速度,φ1,1和φ−1,1是这两个探测单元的仰角。下标k,i表示探测单元的位置,其中编号k为正表示探测单元在涡心上方,为负表示探测单元在涡心下方;i=1表示探测单元到激光雷达的径向距离与左涡心的径向距离相等,而i=2表示探测单元到激光雷达的径向距离与右涡心的径向距离相等。这两个探测单元的径向速度分别
涡心径向距离上两个速度极值点之间的距离一般为20~50 m,而激光雷达到尾流的距离比较远,因此可以假定Vwind·cosφ1,1≈Vwind·cosφ−1,1。通过把式(15)中的两个公式相减,可以抵消掉背景风的影响,得到的速度差可以表示为与Γ1,Γ2相关的量。
并有rc=0.052b0为漩涡的涡心半径,b0是两个漩涡涡心的间距。从式(16)和式(17)可以求出左右漩涡环量Γ1,Γ2
图10 左涡心径向距离上的两个探测单元速度分解Fig.10 Velocity analysis of two detection units above and below the left vortex core
飞机尾流特征参数反演系统对飞机尾流演化规律进行探索研究,包含数据输入、尾流仿真、数据显示、参数反演等模块,可实现雷达回波模拟仿真、数据信息接收和显示、特征参数反演、结果对比分析全过程,系统整体界面如图11所示。下面从仿真数据和实测数据两个方面对系统运行和反演效果进行分析。
飞机尾流仿真可根据需要设置不同的雷达工作和气象参数,以模拟不同的激光雷达型号和不同的天气条件,仿真参数设置界面如图12。
为验证系统性能,飞机尾流左右涡心初始位置设置为(550 m,107 m)和(610 m,105 m),初始环量大小均为400 m2/s,切向风速度设置为2 m/s,湍流耗散率为0.05 m2/s3。雷达模拟仿真参数如表1所示,雷达机械扫描速度为1.4°/s,俯仰角范围为1~15°。
飞机尾流仿真结果如图13所示,对飞机尾流和背景条件进行了仿真模拟,得到设置背景条件下的飞机尾流演化轨迹、理论环量演化和雷达回波数据等,可进一步对飞机尾流特征参数反演性能进行误差分析和对比。
利用仿真模拟得到的飞机尾流探测回波信号进行涡心定位,首先使用Gabor滤波进行粗定位得到初步范围,再利用多普勒信息进行精细定位,结果如图14所示,其中,‘×’ 和‘□’ 代表单独使用Gabor滤波器的左右涡心定位结果,‘△’ 和‘▽’ 代表Gabor和多普勒速度极差联合定位的左右涡心定位结果,‘◦’ 和‘+’ 代表Gabor和多普勒速度平方和联合定位的左右涡心定位结果,‘▷’ 和‘◁’ 代表Gabor和多普勒速度绝对值之和联合定位的左右涡心定位结果。可以发现Gabor滤波、多普勒速度极差、多普勒速度平方之和和多普勒速度平方和等方法,在相对干净的背景湍流和噪声环境条件下可以对涡心进行准确的定位。
为比较不同涡心定位方法的定位效果,将涡心位置与真实位置的误差与飞机翼展进行对比,计算公式为
利用尾流回波仿真数据(仿真参数设置如4.1节所述),对上述4种涡心定位方法进行误差对比计算,单独使用Gabor滤波器、Gabor和多普勒速度极差联合定位、Gabor和多普勒速度平方和联合定位、Gabor和多普勒速度绝对值之和联合定位的左涡心定位误差分别为0.09b0,0.07b0,0.11b0和0.32b0;相应的右涡心定位误差分别为0.10b0,0.07b0,0.13b0和0.24b0,结果证明使用联合方法即先用Gabor滤波器进行粗定位,再用多普勒速度极差进行精细定位的效果最好。
图11 飞机尾流特征参数反演系统界面Fig.11 Interface of wake vortex parameter-retrieval system
图12 飞机尾流仿真参数设置界面Fig.12 Interface of wake vortex simulation parameter setup
进一步可利用涡心定位结果对飞机尾流的速度环量进行估计,采用3.4.1节中的路径积分算法进行计算,并将反演值与理论真实值对比,结果如图15所示,可以发现反演环量与理论真实值吻合很好。
为比较不同速度环量计算方法的反演效果,引入相对误差,计算公式为
路径积分方法、最优化方法和多普勒速度极差方法得到的左涡环量值估计误差分别为11.1%,8.65%和17.32%;相应的右涡环量值估计误差分别为8.88%,7.27%和15.45%,结果证明使用路径积分和最优化方法可以得到较好的估计精度,多普勒速度极差方法因受到复杂背景风场演化的影响,涡心定位发生偏离导致环量估计误差较大。但是最优化方法计算量较大,不利于环量的实时反演。综合分析,路径积分算法可以同时兼顾准确率和实时性的要求,是一种值得推荐的环量反演方法。
表1 激光雷达仿真参数Tab.1 Simulation parameters of the Lidar
图13 飞机尾流多普勒速度和涡心轨迹的时间演化Fig.13 Evolution of wake vortex Doppler velocity and vortex-core trajectory
图14 飞机尾流涡心定位结果Fig.14 Results of wake vortex core location
图15 飞机尾流速度环量的理论与估计值对比Fig.15 Comparison of wake vortex theory and estimated circulation
飞机尾流特征参数反演系统可对不同型号的激光雷达数据进行接收并进行涡心定位和速度环量估计,目前已在香港国际机场、青岛流亭国际机场、长沙黄花国际机场等地开展了实地探测并进行验证。以香港国际机场探测实验为例,采用Wind-Cube200s多普勒激光雷达(雷达参数如表2所示),雷达布置于亚洲国际博览馆屋顶,与跑道垂直距离400 m左右,如图16所示。
系统可对飞机尾流探测实验数据应用不同的方法进行涡心位置和速度环量的估计,系统操作界面如图17所示。
以最优化方法为例,对飞机尾流涡心位置和速度环量进行估计,图18给出飞机尾流涡心定位轨迹,图19给出飞机尾流速度环量演化,可以发现,飞机尾涡在向下飘降的过程中,受到复杂风场的作用影响,左涡演化速度较快,同一RHI扫描反演得到的左右尾涡速度环量会存在差异;同时左涡环量会有一个增大后再减小的趋势,主要由于复杂背景导致的涡心位置和环量的估计误差,但是环量的起伏不剧烈,说明反演结果总体比较平稳。反演结果显示系统可对尾流进行有效的涡心定位和环量估计,为后续尾流的演化规律认知研究提供了数据支撑。
图16 香港机场飞机尾流激光雷达探测示意图Fig.16 Lidar detection scene at Hongkong international airport
本文基于激光雷达探测构建了飞机尾流特征参数反演系统,对代表飞机尾流位置和强度的特征参数进行反演,为机场和航空安全管制提供了手段,主要结论如下:
图17 飞机尾流速度环量估计方法设置Fig.17 Interface of wake vortex circulation estimation algorithm
(1) 对飞机尾流回波信号进行仿真,应用仿真回波信号验证了系统涡心定位和环量估计算法的可行性和有效性;
图18 飞机尾流实测数据涡心位置演化Fig.18 Retrieval results of wake vortex core trajectory from detected data
图19 飞机尾流实测数据速度环量估计Fig.19 Retrieval results of wake vortex circulations from detected data
(2) 飞机尾流特征参数反演系统可高效实现飞机尾流的探测数据处理,包括激光雷达探测数据的导入、飞机尾流涡心位置的定位和速度环量的估计;
(3) 得到的飞机尾流特征参数可进一步开展尾流演化规律的认知研究和危害评估,为飞机尾流动态起降间隔标准制定提供支撑。