融合气象观测数据的高精度大气数据估计算法研究

2022-12-26蒋保睿肖地波王志强

计算机测量与控制 2022年12期

蒋保睿，肖地波，张勇，王志强，林茜

(1.成都信息工程大学自动化学院，成都 610225；2.南京航空航天大学无人机研究院，南京 211106；3.湖南华南光电(集团)有限责任公司，湖南常德 415005)

0 引言

飞行器运动方向与其周围空气的运动方向间的夹角是重要的飞行数据之一，机体与空气的相对速度大小体现在空速与马赫数上，与空气的相对速度方向体现在攻角与侧滑角上[1]，攻角、侧滑角、马赫数、动压和静压统称为大气数据。大气数据系统(ADS, air data system)可以对上述大气数据进行测量，目前普遍使用的探针式大气数据系统由三部分组成：安装在飞机前端的(空速管、角度传感器等)、进行数据解算的计算机、显示与数据输出装置。ADS技术成熟[2]，当前多种民航飞机、战斗机和运输机均广泛采用此系统，如空客与波音的客机、中国的运八运输机，美国的F-35战斗机，俄罗斯的苏-30战斗机上均安装了以空速管为标志的传统大气数据系统[3]。

空速管又名皮托管，可以直接测量总压和静压，动压与静压的比值与马赫数两个之间存在确定的、单调的非线性函数关系，因此通过空速管可以较为简便的求出马赫数。角度传感器安装于飞行器表面，运用风向标原理测量攻角和侧滑角。但空速管和角度传感器突出安装于飞行器顶端的结构存在易受外力损坏、增加雷达反射面积、影响飞行器本身气动性能等问题[4]。并且风标式传感器缺少成熟的冗余解决方案，即使设置了多个传感器互为备份，依然发生过惨痛的事故。如2018年和2019年，PT Lion Mentari航空公司和埃塞俄比亚航空公司发生了两起重大空难，导致数百人死亡。经过漫长的调查，一个根本原因是波音公司的737MAX客机在飞行过程中对攻角传感器的测量异常，这反过来影响了飞行控制系统，并最终导致坠毁[5]。综上所述，确保大气数据的准确测量具有重要意义，直接关系到整个飞行安全。但传统的探针式大气数据系统诸多问题使得人们开始寻求更加可靠的估计方法。

近年来，随着航空事业的发展，各种先进飞行器，如高超声速、隐身等飞行器的出现，使得传统的大气数据系统越来越难以满足飞行器对大气数据的测量要求，多种新型的大气数据系统被提出从数据解算原理进行区分，这些新型大气数据系统可分为嵌入式大气数据系统(FADS, flush air data sensing system)、光学/多普勒大气数据系统、虚拟大气数据系统。

嵌入式大气数据系统依靠嵌入在飞行器前端的测压孔将压力引入到内部的压力传感器上，并通过大气数据计算机从压力传感器上敏感到的压力来计算出攻角、侧滑角、动压和静压等大气数据[6]。FADS的关键技术主要包括测压孔布局、压力分布模型及解算、冗余及故障管理、与其他系统融合等，得到了广泛研究[7-11]。压力分布模型用以描述压力分布于大气数据之间的关系，通过求解该模型可以得到攻角, 侧滑角, 动压, 静压四个参数，剩余大气参数可通过上述四个量计算得到[7]。FADS是典型的多传感器系统，加上现代飞行器工作环境可能越来越恶劣，发生故障的可能性增加，FADS的冗余设计和故障管理也是FADS研究的主要内容之一，FADS的结构使得其冗余备份和故障管理较为方便[8]。由于FADS技术中的压力感受装置内嵌于飞行器内与飞行器表面平齐，因此不会影响气动外形，适用于大马赫数、大攻角的飞行状态，对帮助隐形也有积极作用；同时，由于FADS系统的压力传感器一般置于飞行器仓内，这使其更能适应高超声速飞行和恶劣的飞行环境，因此FADS技术代表了大气数据传感技术未来的发展方向，将在未来各类新型飞行器中得到更为广泛的应用[9-11]。FADS主要由取气装置、引气管路、压力测量单元、数据预处理单元以及一系列软件算法组成，工作过程包括压力传感器数据采集与数据转换、大气数据计算机计算与输出两个部分。经解算得到攻角、侧滑角、动静压与马赫数，通过飞行器表面特定区域的压力分布反推得到飞行参数[12]。

受限于传感器的测量误差、传感器阵列的建模偏差和压力传递延迟，通常FADS计算的马赫数等大气数据误差较大。而广泛运用在导航与控制系统的惯性导航系统(INS, inertial navigation system)恰好具有瞬时精度高、采样周期短的特点，大量研究表明INS与全球导航卫星系统(GNSS, global navigation satellite system)的数据融合能得到高精度的位置和速度信息[13]。但其无法估计有风时的攻角和侧滑角，也无法通过有效的方式获得空速。因此，FADS与INS融合可以用来估计大气数据。肖地波等人将惯性测量元件输出的数据与FADS数据进行融合，建立高维度飞行状态模型，获得长期可靠且准确的飞行数据[14]。Lugo等融合FADS和INS估计火星降落器的大气数据[15]。Prabhu 等提出了大气数据冗余分析方法，在大气数据传感器故障时利用最优扩展卡尔曼滤波算法依靠惯性元件提供大气传感器故障时的测量[16]。蒋保睿等使用FADS与INS的数据，经非线性滤波算法计算，在无风或风速可以忽略的情况下可以得到较为精确的大气数据[17]。贾乾磊等在建立FADS气动模型后，使用模糊理论处理数据，利用聚合运算进行数据融合，提高了FADS的测量精度[18]。但是当风速过大以至不可忽略时，INS根据地速输出数据，与FADS测量的空速间存在差值，导致地速与空速实为不同的速度参数，此时简单的融合已经失去意义。为此，Rhudy等人考虑了风扰动、GPS速度分量和本体速度分量，即考虑了空速、地速和它们的差值，提出了一种无人机姿态、航向和风速的估计算法，结合GPS和ADS、INS的数据，提高实际大气环境下获取大气数据的准确性[19]。Q.Z.He等将定常风场作为弹道速度测量中的未知定常偏差，采用两级EKF进行传感器融合。当风速、攻角、侧滑角传感器工作时，两级EKF估计风速、攻角、侧滑角、风速。当传感器不工作时，实现风速、攻角、侧滑角的信息重建。风场数据的添加完美地实现了飞控系统中关键传感器的信息融合和重构，该系统通常被称为VADS，比传统的ADS更可靠[20]。程胡华等研究了风场随时间的变化规律，用大量实验证明了3.5小时内的观测数据最有利于修正风场模型[21]。

与风速同样属于大气扰动的还有气压等压面的漂移现象，在通常情况下，FADS只能通过当地国际标准气压(ISA, international standard atmosphere)表计算高度，这样的气压表存有高度与气压的对应关系。然而这种关系是在假设空气为理想标准大气条件下得到的，受到气流、温度等条件影响，精度较低。如果飞行器飞行地点的实际大气状况不符合标准的条件，FADS测量的气压便不能正确反映所在地的真实高度，存在“原理误差”。据测量，同为10 000 m处的气压，7月与12月的气压最大相差20 hPa，相当于500米的高度差[22]。而通过观测可以获得各海拔高度对应的气压值，进而补偿原理误差。鲍雪通过实验，补偿了气压高度原理误差的惯性/气压高度计/GPS 组合导航系统，结果显示气压补偿使原有组合导航系统的高度误差明显减小[23]。

本文在上述分析的基础上，提出了新的飞行器大气数据估计方法，该方法为了充分获取可用数据，在飞行器起飞前收集近期的气象观测信息，预测出当地风速和气压高度数据并导入系统，依靠气象、FADS、INS三者信息估计大气数据。在不使用其他设备辅助的条件下，在风速不可忽略、气压高度漂移的情况下估计气流角和马赫数的大小。

1 模型建立

各高度下大气风场包括四个部分：平均风、大气紊流、风切变和突风。平均风是风速的基准值，大气紊流现象的形成和出现与很多因素有关，热交换、地形因素甚至飞机尾流都会形成紊流。风切变是指平均风的时间或空间变化。低空风切变是影响飞机飞行安全的一大因素。突风即俗称的阵风，可以叠加在平均风或紊流上面进行分析。风速观测分为风向观测和风速观测，可以结合上述风速模型，通过风向标、风速计、多普勒测风雷达等设备进行观测[24]。

在周围环境无风或风速已知时，可从INS获得攻角和侧滑角：

(1)

FADS测压孔处测得的压力与当地表面外形、攻角、侧滑角、动压和静压有关，可以描述为[25]：

(2)

其中：pi为第i个孔的压力数据，θi为第i个测压孔处的气流入射角，定义为来流方向与当地表面法线方向的夹角，qc与p∞分别为自由流动压与静压，τ(·)为某一动静压之比对应的修正系数。因为qc/P∞(即动静压之比)与马赫数呈确定的单调函数关系，所以τ(·)可以理解为马赫数的函数τ(Ma)，经牛顿理论可以确定τ(Ma)为一个单调增函数。

设置位于机头附近的五个测压孔组成本文所用的FADS，其中一个孔位于机头中心，其余4个孔均匀分布在其四周，如图1所示。

图1 测压孔位置示意图

分别选择竖直对称轴和水平对称轴上的三个测压孔，使用“三点法”得到局部攻角与局部侧滑角，如下式[26]：

(3)

式中，αe和βe分别为测压孔处的当地攻角和当地侧滑角，往往通过校准得到自由流攻角和侧滑角[26]；A,B,C和A′,B′,C′均可由测压孔位置和实时测量的压力值求得，前者与竖直方向的测压孔有关，后者与水平方向的测压孔有关，其具体表达式见参考文献[26]。

马赫数Ma可由动静压之比Pt/P∞迭代求解[27]。气压高度Hp是静压P∞的单调函数，随P∞的增加呈单调递减的关系，故可以根据任何一值求解另一值；结合温度，可以得到当地音速。获得这些值后，三轴空速可由下式求解：

(4)

(5)

f(·)的具体形式为：

(6)

h是一个线性函数，用于取其自变量的第4-6项，因其自变量X为一个9行的向量，所以此处h(x)为X左乘3行×9列的矩阵，观测矩阵h的具体形式为：

(7)

2 算法设计

因为INS测量的速度为飞行器相对于地面的速度，FADS测量的为空速，根据实时气象观测信息的风速数据推测实际风速的大小和方向，进而使FADS与INS在计算上优势互补。INS测量的速度误差来源于惯性测量元件(加速度计与角速度计)的零偏、初始对准过程中的误差和离散化累加造成的算法误差等；FADS对攻角和侧滑角的测量误差主要来源于压力传感器的精度、压力传感器阵列的建模偏差；对高度、马赫数的测量误差除了源自于传感器测量精度问题外还受到国际标准气压与当地实际气压偏差的影响。风速、大气压强观测误差主要来源于气象观测精度。

用飞行器模型与INS的数据对飞行状态建模，建立系统的状态方程；由FADS和气象观测融合的数据作为观测量，构建滤波模型如式所示，输出高精度的空速矢量与飞行高度，即可计算出攻角、侧滑角、马赫数等大气参数。

算法主要过程如下：

1)FADS大气数据解算：FADS根据压力测量值独立解算出攻角、侧滑角、马赫数和当地音速，攻角和侧滑角采用式所示的三点法公式计算，马赫数和静压采用迭代求解[28]，气压高度根据静压得到，当地空速根据气压高度获取；

2)计算地面坐标系下飞行器相对于空气的速度：根据FADS输出的攻角α、侧滑角β、马赫数Ma、当地音速Vair和INS系统输出的姿态角(φ,ψ,θ)，可以得到飞行器在地面坐标系下的三轴空速分量：

(8)

3)当前风速初步估计：根据飞行前输入系统的气象观测数据，采用外插的方法预测出当前的风速。

4)地面坐标系下飞行器相对于地面的速度估计：由气象观测数据估计的风速(Vwindx,Vwindy,Vwindz)加上地面坐标系下的三轴空速可以得到地面坐标系下飞行器相对于地面的三轴速度。

5)计算飞行器高度和马赫数：将实际测量的海拔高度与大气静压的对应关系带入算法中，FADS实时测量气压的大小，并由气象预报信息构建高度与气压的对应关系，达到对高度值的校准功能。算法结构如图 2所示。

图2 高度校正示意图

图中FADS计算得到的静压数据P∞可通过ISA，由静压与高度的关系换算出高度HFADS，此高度精度较低。另一方面，气象观测数据中提取出观测的高度与P∞的对应关系，与静压P∞进行对比，得到气象预报与FADS数据结合的HWEA。

(9)

(2)计算增益K并更新方差P0：

(10)

图3 融合气象观测估计算法结构图

根据INS测量值、FADS测量值和气象观测数据构建式(5)所示的状态方程和观测方程，通过一个扩展卡尔曼滤波的过程来获取飞行器相对于地面三轴速度的精确估计值，具体的EKF滤波过程如下：

(1)设定滤波初值。令k=1，初始时刻的状态滤波误差协方差矩阵P1取随机的9维对角矩阵，状态量X1的各项取值使用已知的系统状态初值。

(2)对于k-1时刻的状态滤波误差协方差矩阵Pk-1，经下式得到预测协方差矩阵：

(11)

(3)由k-1时刻的输入向量、状态向量Xk-1和状态转移方程估计k时刻的状态Xk∣k-1：

Xk|k-1=f(Xk-1,uk-1)

(12)

估计观测值：

Zk=HkXk|k-1

(13)

(4)求卡尔曼增益：

(14)

(5)更新状态与协方差矩阵：

Pk=(In-KHk)Pk|k-1

(15)

(6)重复循环上述2～4步的过程。

(16)

上述过程在FADS和INS信息融合估计大气数据的过程中，融入了气象观测数据中包含的大气风、气压等大气信息，在不增加额外的机载设备的情况下，利用了更多的有用信息，理论上能够获得更加精确的大气数据估计值。

3 仿真验证

为了验证上述算法的有效性，采用数值仿真对算法进行验证。假设飞行过程包括上升、巡航和下降过程，飞行高度0.01～20 km(0.01 km以下未采集数据)，马赫数0.1～7(0.1马赫以下未采集数据)，飞行时间1 000 s。

风速考虑平均风和大气紊流，平均风在x/y/z轴上的分量在-30～30 m/s中变化，紊流由文献[29]的方法进行建模。压力传感器噪声和INS系统噪声均建模为一阶马尔科夫过程[30]。在同一水平面气压观测噪声为15 Pa(3σ)，风速观测数据采用飞行前3.5小时观测的数据，观测噪声为高斯白噪声。

因气象观测数据分为风速数据与气压数据，在仿真实验中先分析风速数据融合后的滤波效果，再探究完整的气象观测数据融合后的效果。

3.1 风速融合仿真实验

如图 3所示，在上述仿真条件下，融合风速观测与未融合任何气象观测数据的EKF滤波估计方法[31]估计得到的攻角与侧滑角估计曲线如图 4所示，对应的误差曲线如图 5所示，马赫数及其误差曲线如图 6所示。

图4 融合气象观测风速信息的角度估算曲线

图5 融合气象观测风速信息的角度误差曲线图

图6 融合气象观测风速信息的马赫数误差曲线

从图4～图6中可以看出：考虑风速情况下，在融合了气象观测数据中的风速后，本文提出的融合估计算法能够在整个飞行过程中获取攻角、侧滑角和马赫数的估计值；并且，相对于仅融合FADS/INS的大气数据估计算法，本文所提出的融合了气象观测数据的风速后的算法在攻角、侧滑角和马赫数的估计上，精度都有了明显提高。

为进一步定量分析误差情况，计算出攻角、侧滑角、马赫数估计误差的最大值(误差绝对值的最大值)和平均值(误差绝对值的平均值)，分别如表1～3所示。

表1 攻角误差统计表(风速融合)

从表1和表2中的数据可以得到入下结论：相对于仅仅融合FADS和INS数据，融合了气象观测数据后，最大误差降低了四分之一，平均误差降低了四分之三，融合气象观测数据对于提高估计精度有明显作用。

表2 侧滑角误差统计表(风速融合)

通过表3的数据，可以看出融合风速观测对马赫数的估计精度提升明显。

表3 马赫数误差统计表(风速融合)

以上为攻角，侧滑角，马赫数与气压高度的统计特性。通常，飞行器大气数据数据系统获取的大气数据还包括静压、动压和气压高度，其中气压高度和静压存在一一对应关系，而马赫数、动压和静压三项可以在获取其中两项的同时求得最后一项。

3.2 气压融合仿真实验

图7 融合气压观测信息前后马赫数及误差

图8 融合气压观测信息前后气压高度及误差

从图7和图8中可以看出，融合了气象观测数据的气压信息后，马赫数误差变化不大，但气压高度估计误差有了明显降低。

为进一步定量分析气象数据中的气压高度信息的对马赫数和高度估计的影响，统计飞行过程中马赫数和气压高度估计误差的绝对值的最大值和平均值，分别如表 4和表 5所示。

表4 马赫数误差统计表(气压融合)

表5 气压高度误差统计表(气压融合)

从表4和表5中可以看出，融合气息观测数据中的气压信息后，马赫数估计最大误差和平均误差分别降低了17%、10%，说明气象观测数据中的气压高度信息是影响马赫数估计的因素之一；气压高度估计误差最大值降低了62%，平均值降低了87%，说明气象观测数据中的气压信息对于提高高度估计精度具有重要意义。

3.3 完整气象观测数据融合仿真实验

在相同的仿真条件下，同时融合气象观测数据的大气风和气压高度信息后，本文提出的FADS/INS/气象观测数据融合算法估计得到的攻角、侧滑角与仅融合大气风时估计得到的攻角和侧滑角相同(如和所示)，融合完整的气象观测数据(气压、风)与仅仅融合了气象观测数据中的大风信息，攻角和侧滑角误差差别不大，说明大气风是影响融合算法攻角和侧滑角估计的主要因素。

马赫数估计误差如表 6所示。

表6 马赫数误差统计表(气压与风速融合)

由表 6可以看出：融合了完整的气象观测数据后，马赫数估计误差相对于未融合气象观测数据有明显降低。与表 3和表 4比较可以看出，融合了完整的气象观测数据之后，马赫数估计误差相对于仅融合风速或气压观测，误差均有明显降低，说明大气风和大气压都会影响融合算法马赫数估计。

马赫数估计误差如图 9所示。

图9 融合气象观测信息的气压高度及误差曲线

从图9中可以看出，融合了完整的气象观测数据后，气压高度估计误差相对于未融合气象观测数据有明显降低。与图 8相比可以看出，融合完整气象观测数据的气压高度估计误差相比仅仅融合气压观测数据的高度估计误差有明显降低，说明大气风也是影响气压高度估计的因素。

为进一步测试算法的有效性和稳定性，采用蒙特-卡诺方法对算法进行验证。在相同的输入时，未融合气象观测数据、仅融合气象观测数据中的大气风信息、仅融合气象观测数据中的气压信息、融合完整的气象观测数据时的将算法各运行100次，并统计各个算法获得的攻角、侧滑角、马赫数和气压高度最大误差(绝对值)，如表 7所示。表中，仅融合气压观测数据时，仅仅使用图 2所示的流程估计气压高度和马赫数，未计算攻角和侧滑角。仅融合风速观测数据时，未计算气压高度。