0 引言
在飞机飞行阶段，大气湍流是一个影响航空安全的重大隐患。通常会导致飞机突然摇晃或颠簸，严重的湍流可能会导致飞机失去控制而偏离正常的飞行状态，造成结构损坏或伤亡事件[1]。2019年7月11日，加拿大航空公司的AC33客机在温哥华飞往悉尼途中遭遇湍流，短时间内迅速下降数百米，舱内设备严重损坏，多人受伤，最终迫降在檀香山机场。因此，为避免严重损失，探究大气湍流的检测方法尤为重要。
为保障飞行安全，民航飞机通常利用机载气象雷达探测前方路径是否存在湍流。目前，国际民航界统一标准规定速度谱宽大于5 m/s的气象目标为湍流[2]。但是，湍流对飞机的危害程度与飞机特有属性密切相关，谱宽为5 m/s的湍流，对翼载荷大的飞机一般不会构成威胁，若对飞机告警会造成不必要的绕飞，降低飞行效率。但是翼载荷小的飞机操控性差，若不告警会导致不必要的伤亡事故[3]。
2016年3月，美国航空无线电技术委员会(RTCA)修订了机载气象雷达最低运行性能标准DO-220A，要求在湍流检测时考虑飞机本身的因素[4]，因此引入垂直载荷因子量化湍流对不同类别飞机的影响。有学者提出可从ADS-B(automatic dependent surveillance-broadcast，简称广播式自动相关监视)数据中提取飞机航行中的高度数据和地速信息[5]，基于此数据可估算飞机的垂直载荷因子。ADS-B数据无需机组干预可自动传输，易于获取且精确度高、更新频率快[6]，因此使用ADS-B数据估算垂直载荷因子成本低、效率高、精度准，有利于研究大气湍流对飞机的影响。
为提升我国机载雷达湍流告警技术水平，基于DO-220A标准，本文研究基于垂直载荷因子的湍流检测方法，并结合飞机ADS-B数据进行验证和分析。首先构建湍流模型，根据飞机对湍流的频响函数求得飞机比例因子，然后结合湍流谱宽，估算飞机垂直载荷因子。在满足DO-220A检测标准的前提下，根据贝叶斯准则，确定基于统计特性的湍流检测门限。最后将应用ADS-B数据估算的飞机垂直载荷因子与检测门限作比较，判断是否进行告警处理。本文详细研究梳理了基于垂直载荷因子的湍流告警方法的实现步骤，准确、可靠地量化了湍流对不同机型的危害程度。
1 湍流检测方法的流程
湍流是一种由漩涡和垂直气流引起的不规则的空气运动，会产生垂直载荷作用于飞机上，导致飞机轻微颠簸或失去控制造成结构破坏。因此提前检测湍流并告警对保证飞行安全至关重要。
NASA TRAWS早期计划确定了量化湍流对飞机影响的度量标准——垂直载荷因子，即采用 5 s内飞机垂直加速度的均方根来量化湍流的危险性，记为σΔn[7]。这种基于雷达观测数据预测飞机垂直载荷的算法结构[7]近似如下：
(1)
式中，σΔn/unitσw表示飞机比例因子，
表示雷达湍流谱宽相关量，
是理论上雷达脉冲体积的补偿系数。
为后文叙述方便，引入以下定义[7]：
z=xy
(2)
式中：z表示飞机垂直载荷因子的估计值，量化了湍流对飞机的危害程度；x表示飞机比例因子的估计值，与飞机特有属性相关;y表示经过脉冲体积补偿的雷达湍流回波的速度谱宽估计值[7]，一般采用脉冲对处理(PPP)方法估算谱宽[3,8]。
以下将结合飞机ADS-B数据，分析飞机对湍流的响应，求解飞机比例因子x,根据式(2)求得飞机的垂直载荷因子。继而根据DO-220A检测标准，基于统计数据和贝叶斯准则确定湍流检测的最佳门限。最终实现了结合ADS-B数据的基于垂直载荷因子的湍流检测，检测流程如图1所示。
图1 湍流检测流程图
2 基于ADS-B数据求解垂直载荷因子
求解飞机垂直载荷因子主要是求解飞机比例因子和回波谱宽。飞机比例因子依赖于飞机对湍流的响应，因此需建立湍流风场模型和飞机模型[3]，应用气动、飞机力学模型等相关知识求解飞机对湍流的频响函数[3,9]，继而求得飞机比例因子，同时考虑谱宽因素，计算出飞机垂直载荷因子。
2.1 湍流模型及其功率谱密度函数
量化湍流对飞机的影响时，首要考虑因素是飞机系统的输入，即建立三维湍流场。对飞行中的飞机来说，湍流是一种方向和强度均有明显变化的阵风[3]，为简化分析，单独考虑机翼对对称垂直阵风分量的响应，将三维空间的问题简化成一维的形式。Von Karman模型连续阵风的功率谱密度(PSD)函数表示如下[9]：
(3)
式中，σ表示湍流强度，L表示湍流尺度，V表示风速。
2.2 垂直载荷因子
飞机对湍流的响应十分复杂，一般在频域采用功率谱法求解。假设质量为m的飞机在遭遇湍流前处于升力等于重力的配平状态，遇湍流后只作沉浮运动[3]。飞机主升力面进入阵风的瞬间有效攻角发生变化，飞机沉浮运动响应和阵风速度引起的升力也随之改变[7]。根据牛顿第二定律建立飞机浮沉运动方程，得到任意频率点处飞机垂直加速度响应和垂直阵风速度之间的频响函数[4]为
(4)
式中，ρ为空气密度，VTAS表示飞机真空速，Sw表示机翼面积，CL为相对机翼面积的全机升力线斜率。
因此，单位垂直阵风响应的均方根载荷，即飞机比例因子[3]可表示为
(5)
当湍流回波的谱宽已知时，根据式(2)可估算出飞机的垂直载荷因子。垂直载荷因子的计算流程如图2所示。
图2 垂直载荷因子的计算流程
根据垂直载荷因子取值范围对湍流强度进行如表1所示的分类[7]。
表1 湍流强度的分类
z取值范围湍流强度z≤0.1g轻度湍流0.1g
2.3 应用ADS-B数据
ADS-B是一种全新的技术，装有ADS-B系统的飞机可以实时提供自身的精确位置和飞行高度、地速、航向等信息，从中提取高度和地速数据求解飞机的实时重量W和升力系数CL，代入式(4)，结合式(5)，可估算飞机实时比例因子。
利用ADS-B数据中的高度数据h计算飞机的爬升速率Vroc，结合地速Vgr计算飞机真空速VTAS和飞行航迹角γ[5,10]:
(6)
(7)
(8)
飞机实时重量是非公开数据，一般的估算方法复杂且仅使用于起飞阶段，本文采用一种简单并适用于整个航程的估算方法[10]。此方法将总重量分为3个部分：空重Wempty、燃料重量Wfuel和有效载荷重量Wpayload。即
W=Wempty+Wfuel+Wpayload
(9)
空重已知，其他两项计算如下:
(10)
(11)
式中,mcr表示每秒钟油耗重量，Stotal表示飞机航行距离，Vaverage表示飞机巡航时的平均速度，Sto表示飞机起飞时的测量距离，Stomax表示飞机最大起飞距离，Wmplw表示飞机最大有效载荷。
升力系数CL也会影响飞机对湍流的响应，按式(12)进行估算，式中需将飞机重量W的单位由kg转化为N。
(12)
根据上述参数可估算出飞机实时垂直载荷因子。
3 基于统计特性确定检测门限
湍流检测时，需将飞机垂直载荷因子与检测门限比较，若大于此时的门限则告警。确定湍流检测门限一般采用“重量输入”法和“通用”法[7]。“重量输入”法依赖飞机实时翼载荷，但实时翼载荷是不断变化的，会增加系统运算复杂度，所以在求解门限时更倾向于使用“通用”法。“通用”法以飞机起降重量的统计数据为基础，在一定的误差范围内可消除对飞机实时翼载荷的依赖。在特定的飞行条件下，“通用”法求解步骤如下。
1)求解飞机比例因子x的统计特性
根据飞机起降实际重量的统计数据，求解x的概率密度函数(PDF)fx(x)[7]。fx(x)服从高斯分布，均值用μx表示，标准差用σx表示。
2)求解湍流谱宽y的统计特性
估算湍流的谱宽时通常采用时域的PPP方法。可根据文献[8]求解PPP法估计谱宽时的性能。若一定范围内估计的湍流平均多普勒速度谱宽表示为μy，标准差表示为σy，y的PDF也服从高斯分布，一般表示为[3]
(13)
3)求解飞机垂直载荷因子z的统计特性
已知fx(x)、fy(y)服从高斯分布，由概率论相关知识可知z的PDF也服从高斯分布[11]。x，y，z的统计特性具有如下的关系[7]：
μz=μxμy
(14)
(15)
若fx(x)、fy(y)已知，由式(14)、式(15)可计算z的PDF。
4)确定检测门限
根据假设检验知识，假设H0为湍流不存在，H1为湍流存在。根据垂直载荷因子的统计特性，在满足DO-220A检测条件和其他先验知识的前提下，基于贝叶斯准则[12]，可计算得出最佳检测门限。也就是说在已知假设H0和假设H1的先验概率以及各种判决代价因子给定的情况下，寻找使得判决所付出的平均代价最小的判决门限[3]。确定检测门限的流程如图3所示。
图3 确定基于统计特性的湍流检测门限流程
4 验证与分析
DO-220A中基于翼载荷标准定义了A、B、C三类飞机，本节以机型为B777-200LR的A类飞机为例，对湍流检测方法进行仿真分析。首先根据2.3节内容，利用飞机ADS-B数据计算飞机真空速VTAS、实时重量W和升力系数CL，代入式(4)求解飞机对湍流的频响函数，结合Von Karman阵风模型，根据式(5)求解飞机实时比例因子，仿真所需参数如表2所示。假设已知雷达测得湍流目标的谱宽为5 m/s，求解垂直载荷因子。
表2 比例因子仿真参数
仿真参数参数值湍流强度σ/(m·s-1)1湍流尺度L/m762机翼面积Sw/m2427.8空重Wempty/kg148181油耗系数mcr/(kg·s-1)-2巡航平均速度Vaverage/(m·s-1)250起飞测量距离Sto/m2970最大起飞距离Stomax/m3536最大有效载荷Wmplw/kg201912
若已获得该飞机比例因子x的统计数据，且假设漏警是虚警代价的2.5倍，根据DO-220A检测性能标准，采用“通用”法求“最佳”湍流检测门限。
4.1 估算垂直载荷因子
在flightradar24网站下载B777-200LR的ADS-B数据如图4所示。由图4(a)可见，飞行途中t1=29 590 s到t2=29 673 s的83 s内飞机高度迅速下降了约600 m，这段时间前后，图4(b)中飞机地速在252.9 m/s到263.7 m/s之间波动。
(a)飞机实际飞行高度
(b)飞机实际地速
图4 飞机ADS-B数据
根据式(6)、式(7)计算飞机的爬升速率和真空速如图5所示。由图5(a)可见，忽略飞机的起降情况，t=29 660 s时飞机的爬升速率约为-8 m/s,有明显的变化且为负值，表明飞机此时遭遇湍流迅速下降。
(a)飞机爬升速率
(b)飞机真空速
图5 飞机爬升速率和真空速
根据式(9)、式(12)计算飞机实时重量和升力系数如图6所示，飞行中燃料会损耗，因而飞机重量随之减小，遭遇湍流时，升力系数也减小了约0.07。
(a)飞机重量
(b)飞机升力系数
图6 飞机重量和升力系数
图7、图8分别显示了根据式(3)、式(5)得到的Von Karman模型连续阵风的PSD和此模型下飞机的比例因子。由图8可知飞机在遇到湍流时，比例因子迅速增大到0.039 02 g/(m/s)。
图7 Von Karman阵风功率谱密度
图8 飞机比例因子
由式(2)计算t=29 660 s时飞机的垂直载荷因子，z=xy=0.039 02*5=0.195 1g。根据表1的分类标准，飞机遭遇的是中度湍流。
4.2 确定湍流检测门限
1)比例因子的统计特性
根据不同飞行高度下飞机比例因子统计数据的均值和标准差表格[7]，选取飞机遭遇湍流时的飞行高度，可确定比例因子x的均值μx= 0.041 26,标准差σx=0.003 716。
2)谱宽的统计特性
根据DO-220A标准可知，对于A类飞机，湍流不存在和存在时的垂直载荷因子z的均值分别为0.1和0.3。由步骤1)已知x的均值和标准差，那么y的均值可由式(14)计算得出。假设雷达仿真参数设置如表3所示，根据PPP法可估计谱宽的标准差。
表3 雷达仿真参数
主要参数数值脉冲数/个8脉冲重复频率/Hz3000雷达波长/m0.032信噪比/dB10
3)垂直载荷因子的统计特性及门限的确定
x,y的统计特性已知，由式(15)确定z的标准差，最终可得出两种假设条件下z的概率密度函数。表4列出了H0和H1假设下x,y,z的统计特性。
表4 x,y,z的统计特性
假设H0假设H1μz(×g)0.10.3σz(×g)0.04410.0950μx(×g/(m·s-1))0.04120.0412σx(×g/(m·s-1))0.003710.00371μy/(m·s-1)2.42407.2710σy/(m·s-1)1.04252.1997
H0和H1假设下垂直载荷因子的概率密度函数如图9所示。
图9 假设H0和假设H1条件下垂直载荷因子的PDF
在检测概率不低于85%、虚警概率小于20%的前提下，如图10(a)所示，可以确定检测门限zt的范围是0.138g～0.201g。当H0和H1的先验概率相等、漏检代价是虚警概率的2.5倍时，即 P(H0)=P(H1)=0.5，代价因子c01=2.5c10，c00=c11=0时，根据贝叶斯准则，如图10(b)所示，可知当平均代价最小时，可以确定最佳检测门限为0.160 g。
(a)不同检测门限下虚警概率和检测概率关系曲线
(b)不同检测门限下的平均代价
图10 确定基于统计特性的检测门限
由4.1节的仿真结果可知飞机的垂直载荷因子z=0.195 1 g，大于湍流的检测门限0.160 g，因此需向飞行员告警，与报道中的实际情况相符，从而验证了利用ADS-B数据的基于垂直载荷因子的湍流检测方法的合理性。flightradar24网站提供大部分航班的ADS-B数据，应用这些数据能够精确可靠地估算飞机垂直载荷因子，对湍流检测的研究有重要意义。
5 结束语
本文主要研究了基于垂直载荷因子的湍流检测方法。首先构建了湍流模型，基于飞机对湍流的响应求解比例因子，然后结合雷达回波的谱宽，得到垂直载荷因子，量化了湍流对飞机的影响。接着根据DO-220A，推导基于统计特性的检测门限范围,结合贝叶斯准则，确定湍流检测的最佳门限。仿真算例选取遭遇湍流飞机的ADS-B数据求解垂直载荷因子，所求结果大于湍流检测的门限，应向飞行员告警，这与实际情况相符，说明本文研究的湍流检测方法可以实现对湍流的告警。在今后的研究中, 将尝试其他湍流模型，进一步详细分析和比较验证湍流告警性能，并考虑其他飞机模型，更充分地验证湍流检测方法的可靠性。
[1]GOLDING W L.Turbulence and Its Impact on Commercial Aviation[J].Journal of Aviation/Aerospace Education & Research，2002，11(2):19-29.
[2]卢晓光，夏冬.基于统计置信度的湍流检测门限确定方法[J].中国民航大学学报，2011，29(4):27-30.
[3]卢晓光,范源丹,李海,等.基于垂直载荷因子的增强型湍流检测方法仿真[J].现代雷达,2019,41(9):31-36.
[4]RTCA.Minimum Operational Performance Standards(MOPS)for Airborne Weather Radar System[S].DO-220A Washington D.C: RTCA Inc，2016.
J M, KWLATKOWSKI K, HAAN S D, et al.Retrieving Clear-Air Turbulence Information from Regular Commercial Aircraft Using Mode-S and ADS-B Broadcast[J].Atmospheric Measurement Techniques Discussions，2015(11):11817-11852.
[6]范源丹.晴空湍流对飞机的影响分析[J].交通技术,2019, 8(1):1-10.
[7]BOWLES R L , BUCK B K.A Methodology for Determining Statistical Performance Compliance for Airborne Doppler Radar with Forward-Looking Turbulence Detection Capability[R].NASA，NASA/CR-2009-215769，2009.
[8]DOVIAK R J , ZRNIC D S.Doppler Radar and Weather Observations[M].New York: Dover Pubns，2006.
[9]WRIGHT J R , COOPER J E.Introduction to Aircraft Aeroelasticity and Loads[M].New York: John Wiley & Sons，2008.
[10]Gloudemans T.Aircraft Performance Parameter Estimation Using Global ADS-B and Open Data[D].Netherlands:Delft University of Technology, 2016.
[11]盛 骤，谢世千，潘承毅，等.概率论与数理统计[M].北京:高等教育出版社，2001.
[12]赵树杰，赵建勋.信号检测与估计理论[M].北京:清华大学出版社，2005.
