CN113408057B - 一种基于em算法的飞机气流角高精度离线估计方法 - Google Patents

Abstract

本发明公开了一种基于EM算法的飞机气流角高精度离线估计方法，首先对飞行器进行受力分析，利用运动学方程及动力学方程对飞行器进行建模，得到适配的状态空间模型；接下来利用已经建立的飞行器模型，加入不同的状态噪声及过程噪声，得到不同噪声环境下的飞行数据；然后载入仿真数据，运行EM算法，得到不同噪声背景下的气流角估计值，并将估计值与真实值进行比较，通过蒙特卡洛仿真验证EM算法的可行性；最后对EM算法进行封装，使得使用者只需设置好飞行器各项参数即可得到相应的气流角估计值。本发明方法的估计结果精度较高，可以满足飞行器测试者对精确气流角的需求。

Description

一种基于EM算法的飞机气流角高精度离线估计方法

技术领域

本发明属于飞行器技术领域，具体涉及一种飞机气流角高精度离线估计方法。

背景技术

迎角和侧滑角是飞行力学的两个重要的飞行状态参数，也是飞行控制及导航系统所需要的两个关键参数。由于有效载荷和预算有限，许多飞机系统，特别是无人机，没有配备气流角度传感器，使得利用其他相关传感器(例如GPS和INS)的数据来估计两个气流角成为一个亟待解决的问题。

一般情况下，气流角是由安装在飞机上的风标传感器、压差式传感器和零压差式传感器等直接测量得到的，但这类传感器不仅使用成本较高，其由于受到结冰及与飞行状态有关的局部环流的影响，几乎不可避免地会造成很大的零点偏差。迎角的容限偏差对于飞行器而言将对控制律的计算造成很大的误差，并且有可能在某种迎角情况下导致飞机不可控。另外，迎角也是飞机飞控系统和失速告警所需要的重要参数，而侧滑角的偏差将会导致飞机产生不必要的能量损失。

飞行器动态系统存在系统噪声时，气动参数辨识成为了一个典型的状态与参数的联合估计问题。其中状态估计与参数估计之间深度耦合、相互影响、互为因果；状态估计误差会引起参数辨识误差，辨识误差反过来会带来状态估计误差，从而使联合估计问题呈现出较高的复杂性。特别是当系统噪声和量测噪声统计特性未知的情况下，需要从试验数据中提取大量信息，从而进一步使问题复杂化。随着最优估计理论的不断发展，期望最大化算法(EM algorithm)与滤波/平滑算法的结合逐渐趋于成熟，气流角估计问题则可以转化为自适应滤波问题，EM算法的数值稳定性较好，并且近期运用EM迭代优化思想实现估计与辨识的一体化处理也取得了重要的研究进展，因此考虑采用EM算法对系统的过程噪声和量测噪声进行估计，能够有效的减小噪声对辨识结果的影响，大幅提高辨识精度。

发明内容

为了克服现有技术的不足，本发明提供了一种基于EM算法的飞机气流角高精度离线估计方法，首先对飞行器进行受力分析，利用运动学方程及动力学方程对飞行器进行建模，得到适配的状态空间模型；接下来利用已经建立的飞行器模型，加入不同的状态噪声及过程噪声，得到不同噪声环境下的飞行数据；然后载入仿真数据，运行EM算法，得到不同噪声背景下的气流角估计值，并将估计值与真实值进行比较，通过蒙特卡洛仿真验证EM算法的可行性；最后对EM算法进行封装，使得使用者只需设置好飞行器各项参数即可得到相应的气流角估计值。本发明方法的估计结果精度较高，可以满足飞行器测试者对精确气流角的需求。

本发明解决其技术问题所采用的技术方案包括如下步骤：

步骤1：对飞行器进行受力分析，利用运动学方程及动力学方程对飞行器进行建模，得到适配的状态空间模型；

步骤1-1：通过离散体轴动力学方程和旋转运动学方程，得到状态方程如下：

将式(1)改写为：

x_k+1＝f(x_k,u_k)+ω_k (2)

式中

为待估计的未知状态量，其中θ_k,φ_k分别为俯仰角和滚转角，u_k,v_k,w_k分别为空速在体轴系三个方向上的投影；u_k＝[a_xk a_yk a_zk p_k q_k r_k]^T为已知输入量，其中p_k,q_k,r_k分别为滚转角速率，俯仰角速率和偏航角速率，a_xk,a_yk,a_zk分别为加速度传感器测得的体轴系三个方向上的加速度大小；ω_k＝[ω_k(1) ω_k(2) ω_k(3) ω_k(4)ω_k(5)]^T为服从零均值且方差未知高斯分布的状态噪声，即ω_k～N(0,Q)，且Q为正定矩阵；

步骤1-2：假设飞机配备了空速传感器和平台惯性导航系统，那么无风条件下的测量方程建立如式(16)：

将上式改写为：

z_k＝h(x_k,u_k)+ζ_k (4)

其中

为k时刻的量测量，且V为飞行器的空速；ζ_k＝[ζ_k(1) ζ_k(2)ζ_k(3) ζ_k(4) ζ_k(5) ζ_k(6)]^T为服从零均值且方差未知高斯分布的量测噪声，即ζ_k～N(0,R)，且R为正定矩阵；

由式(2)及式(4)共同组成飞行器非线性离散系统状态空间模型；

如果速度分量u_k,v_k,w_k已知，飞机迎角α和侧滑角β由式(5)计算得到：

步骤2：利用步骤1建立的飞行器状态空间模型，加入不同的状态噪声及过程噪声，得到不同噪声环境下的飞行数据，生成测试用仿真数据；

步骤2-1：首先给定一组状态初值，通过龙格库塔法得到其余时刻的状态量，再加上人为设定的噪声，得到含噪声干扰的状态量；

步骤2-2：再通过龙格库塔法计算出一组无噪声干扰的状态量，通过式(5)计算得到飞机迎角和侧滑角的真值；

步骤2-3：将步骤2-2得到的状态量代入式(2)，再加上人为设定的量测噪声即得到带有噪声的观测量；

步骤2-4：在每一次仿真时改变状态噪声和量测噪声，即改变Q和R的大小，得到不同噪声水平下的仿真值；

步骤3：载入步骤2生成的仿真数据，运行EM算法，得到不同噪声背景下的气流角估计值，并将估计值与真实值进行比较，通过蒙特卡洛仿真验证EM算法的可行性；

步骤3-1：将EM算法应用于非线性离散系统状态空间模型中，将状态向量X_N＝{x₁,x₂,…,x_N}作为EM中的不可观测数据，量测向量Z_N＝{z₁,z₂,…,z_N}作为EM中的可观测数据，

为未知统计量；P₁为初始时刻状态量协方差矩阵，

为初始时刻状态量预测均值；

由EM算法，其代价函数为：

其中，

为未知统计量的估计值，

为在给定量测数据Z_N和当前未知统计量估计值

下隐变量数据X_N的条件概率分布，p_μ(X_N,Z_N)为以μ为模型参数建立的状态向量X_N和量测向量Z_N的联合概率分布；

步骤3-2：在给定参数μ的条件下，系统状态和观测值的对数联合概率分布为：

将式(7)代入式(6)，得到：

其中

步骤3-3：通过最大化成本函数推导出未知统计量的最优估计，在高斯分布的假设条件下，采用三阶球面逼近算法对式(9)-式(11)进行计算，统计量μ有解析解，具体如下：

P₁＝P_1∣N (13)

式中

i＝1,2,…4n

其中，N为样本总数；ξⁱ为依据球容积准则选择的西格玛点；n为状态向量的维数；m为量测量的维数；

步骤4：设置飞行器各项参数即能得到相应的气流角估计值。

本发明的有益效果如下：

使用传统气流角辨识算法进行迎角及侧滑角估计时需要状态噪声及量测噪声的方差统计特性，而不准确的甚至错误的噪声假设往往对估计结果的精度影响很大，进而影响最终的估计结果。本方法中用期望最大化算法对噪声的统计特性进行估计，避免了噪声统计特性不准确而导致的估计结果精度降低。相对于滤波，平滑方法因为利用更多的量测信息而获得更高精度的估计，另外，由于容积卡尔曼平滑器具有估计精度高、计算量相对较小等优点，非常适合用于高维状态量的估计问题，本方法摒弃了传统的仅仅使用在卡尔曼滤波的技术路线，转而在滤波的基础上增加了容积平滑器，使得最终估计结果准确度较传统算法有所提升。通过式(2)及式(4)可知，本算法在构建飞行器非线性离散系统状态空间模型时并未使用气动模型和气动导数等相关信息，使得状态空间模型独立于飞行器的气动布局，即在特殊布局下也可正常进行气流角的估计工作，拓展了算法的应用范围，并且为下一步飞行器气动模型参数的辨识打好基础。

附图说明

图1为本发明实施例生成仿真数据所使用的输入量变化图像。

图2为本发明实施例进行飞行仿真得到的状态量及量测量变化图像。

图3为本发明实施例使用单次仿真数据进行气流角重构的结果。

图4为本发明实施例进行100次蒙特卡洛仿真并使用均方误差对每一次的重构结果进行评估得到的结果。

图5为本发明实施例与传统EKF算法进行对比得到的结果。

图6为本发明实施例的噪声辨识结果(以Q_u,Q_θ,R_V和R_θ为例)。

图7为本发明实施例使用单次真实飞行数据进行气流角重构的结果。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

本发明是利用EM算法较好的数值稳定性及容积卡尔曼平滑器(CKS)对高维状态下的估计问题具有辨识精度高、稳定性好的特点，从而建立和推导一种基于EM和CKS联合估计(EM-CKS)的气流角重构新算法，能够同时解决动态系统中存在的系统方程形式复杂并且非线性程度高、噪声总是存在且无法消除等问题，为飞行器迎角及侧滑角估计提供有效的工具。

本发明目的在于为飞行器提供一种气流角估计的新思路，传统的估计方法是基于非线性运动学和测量模型的扩展卡尔曼滤波器(EKF)，而卡尔曼滤波法的本质特征是采用状态空间法描述系统模型，它与实际动态系统的误差包含系统噪声和量测噪声，即过程噪声和量测噪声的统计特性要作为先验信息参与估计。而这些统计特性通常是通过传感器地面测试或专家经验获得的，存在一定的偏差，导致最后的估计结果精度较低。而本发明在噪声未知的情况下很好地对气流角进行估计，且结果精度较高。

一种基于EM算法的飞机气流角高精度离线估计方法，包括如下步骤：

步骤1：对飞行器进行受力分析，利用运动学方程及动力学方程对飞行器进行建模，得到适配的状态空间模型；

步骤1-1：通过离散体轴动力学方程和旋转运动学方程，得到状态方程如下：

将式(1)改写为简明形式：

x_k+1＝f(x_k,u_k)+ω_k (2)

式中

为待估计的未知状态量，其中θ_k,φ_k分别为俯仰角和滚转角，u_k,v_k,w_k分别为空速在体轴系三个方向上的投影；u_k＝[a_xk a_yk a_zk p_k q_k r_k]^T为已知输入量，其中p_k,q_k,r_k分别为滚转角速率，俯仰角速率和偏航角速率，a_xk,a_yk,a_zk分别为加速度传感器测得的体轴系三个方向上的加速度大小；ω_k＝[ω_k(1) ω_k(2) ω_k(3) ω_k(4)ω_k(5)]^T为服从零均值且方差未知高斯分布的状态噪声，即ω_k～N(0,Q)，且Q为正定矩阵；

步骤1-2：假设飞机配备了空速传感器和平台惯性导航系统，那么无风条件下的测量方程建立如式(16)：

将上式改写为简洁形式：

z_k＝h(x_k,u_k)+ζ_k (4)

其中

为k时刻的量测量，且V为飞行器的空速；ζ_k＝[ζ_k(1) ζ_k(2)ζ_k(3) ζ_k(4) ζ_k(5) ζ_k(6)]^T为服从零均值且方差未知高斯分布的量测噪声，即ζ_k～N(0,R)，且R为正定矩阵；

由式(2)及式(4)共同组成飞行器非线性离散系统状态空间模型；

如果速度分量u_k,v_k,w_k已知，飞机迎角α和侧滑角β由式(5)计算得到：

步骤2：利用步骤1建立的飞行器状态空间模型，加入不同的状态噪声及过程噪声，得到不同噪声环境下的飞行数据，生成测试用仿真数据；

步骤2-1：首先给定一组状态初值，通过龙格库塔法得到其余时刻的状态量，再加上人为设定的噪声，得到含噪声干扰的状态量；

步骤2-2：再通过龙格库塔法计算出一组无噪声干扰的状态量，通过式(5)计算得到飞机迎角和侧滑角的真值，，作为衡量通过本估计方法得到的气流角重构结果优劣的真值；

步骤2-3：将步骤2-2得到的状态量代入式(2)，再加上人为设定的量测噪声即得到带有噪声的观测量；

步骤2-4：在每一次仿真时改变状态噪声和量测噪声，即改变Q和R的大小，得到不同噪声水平下的仿真值；

通过状态方程及量测方程可以看出，本实施例是在无风、系统噪声及量测噪声未知且飞行器未配备气流角传感器的情况下进行的；

步骤3：载入步骤2生成的仿真数据，运行EM算法，得到不同噪声背景下的气流角估计值，并将估计值与真实值进行比较，通过蒙特卡洛仿真验证EM算法的可行性；

步骤3-1：将EM算法应用于非线性离散系统状态空间模型中，将状态向量X_N＝{x₁,x₂,…,x_N}作为EM中的不可观测数据，量测向量Z_N＝{z₁,z₂,…,z_N}作为EM中的可观测数据，

为未知统计量；

由EM算法，其代价函数为：

其中，

为未知统计量的估计值；

步骤3-2：在给定参数μ的条件下，系统状态和观测值的对数联合概率分布为：

将式(7)代入式(6)，得到：

其中

步骤3-3：通过最大化成本函数推导出未知统计量的最优估计，在高斯分布的假设条件下，采用三阶球面逼近算法对式(9)-式(11)进行计算，统计量μ有解析解，具体如下：

P₁＝P_1∣N (13)

式中

i＝1,2,…4n

其中，N为样本总数；ξⁱ为依据球容积准则选择的西格玛点；n为状态向量的维数；m为量测量的维数；

为了更加简洁明了，本文提出的算法步骤归纳如下：

①给出

P₁,Q和R的初值；

②执行容积卡尔曼滤波(CKF)和容积卡尔曼平滑(CRTSS)过程，计算

和

③执行EM算法，更新

P₁,Q和R；

④重复步骤①到步骤③，直到达到收敛条件；

⑤通过式(5)计算迎角及侧滑角的值；

步骤4：设置飞行器各项参数即能得到相应的气流角估计值。

如图1所示，为本发明实施例生成仿真数据所使用的输入量变化图像，图2为本发明实施例进行飞行仿真得到的状态量及量测量变化图像。

由图3可以看出，使用随机生成的单次飞行数据对算法进行验证，可以得到较为理想的迎角及侧滑角估计结果。但是，单次试验具有偶然性，不能充分说明CKS-EM算法对于气流角重构这一问题具有良好的适用性。因此，为了消除单次试验的偶然性，增加本算法的可靠性和可信度，本算例中进行了100次蒙特卡洛仿真，并使用MSE(均方误差)对每一次的传感器误差参数估计结果进行评估。通过图4可以看出，在一百次蒙特卡洛仿真中，无论是迎角还是侧滑角，每一次的均方误差的量级均在小数点后九位，说明对于任何一组随机生成的气流角，CKS-EM算法都可以很好地重构出十分接近真值的迎角及侧滑角，从而不仅验证了通过本算法可以得到理想的气流角重构值，同时也证明了本算法具有一定的通用性。另外，由于，及对于任何外形的飞行器都是成立的，因此本文所阐述的使用CKS-EM算法进行气流角重构这一方法是具有一定可推广性的。

图5展示了分别将系统噪声方差阵和量测噪声方差阵扩大10倍及20倍，代入EKF算法进行滤波，得到接近真值的状态量，再通过得到重构的迎角及侧滑角。将上述过程重复100次(即进行100次蒙特卡洛仿真)，得到每一次的均方误差，再与CKS-EM算法得到的结果进行对比。图6为本发明实施例的噪声辨识结果(以Q_u,Q_θ,R_V和R_θ为例)。

从图7可以看出本算法针对真实飞行数据也具有良好的估计性能。

Claims (1)

1.一种基于EM算法的飞机气流角高精度离线估计方法，其特征在于，包括以下步骤：

步骤1：对飞行器进行受力分析，利用运动学方程及动力学方程对飞行器进行建模，得到适配的状态空间模型；

步骤1-1：通过离散体轴动力学方程和旋转运动学方程，得到状态方程如下：

将式(1)改写为：

x_k+1＝f(x_k,u_k)+ω_k (2)

式中

为待估计的未知状态量，其中θ_k,φ_k分别为俯仰角和滚转角，u_k,v_k,w_k分别为空速在体轴系三个方向上的投影；u_k＝[a_xk a_yk a_zk p_k q_k r_k]^T为已知输入量，其中p_k,q_k,r_k分别为滚转角速率，俯仰角速率和偏航角速率，a_xk,a_yk,a_zk分别为加速度传感器测得的体轴系三个方向上的加速度大小；ω_k＝[ω_k(1) ω_k(2) ω_k(3) ω_k(4) ω_k(5)]^T为服从零均值且方差未知高斯分布的状态噪声，即ω_k～N(0,Q)，且Q为正定矩阵；

步骤1-2：假设飞机配备了空速传感器和平台惯性导航系统，那么无风条件下的测量方程建立如式(3)：

将上式改写为：

z_k＝h(x_k,u_k)+ζ_k (4)

其中

为k时刻的量测量，且V为飞行器的空速；ζ_k＝[ζ_k(1) ζ_k(2) ζ_k(3)ζ_k(4) ζ_k(5) ζ_k(6)]^T为服从零均值且方差未知高斯分布的量测噪声，即ζ_k～N(0,R)，且R为正定矩阵；

由式(2)及式(4)共同组成飞行器非线性离散系统状态空间模型；

如果速度分量u_k,v_k,w_k已知，飞机迎角α和侧滑角β由式(5)计算得到：

步骤2：利用步骤1建立的飞行器状态空间模型，加入不同的状态噪声及过程噪声，得到不同噪声环境下的飞行数据，生成测试用仿真数据；

步骤2-1：首先给定一组状态初值，通过龙格库塔法得到其余时刻的状态量，再加上人为设定的噪声，得到含噪声干扰的状态量；

步骤2-2：再通过龙格库塔法计算出一组无噪声干扰的状态量，通过式(5)计算得到飞机迎角和侧滑角的真值；

步骤2-3：将步骤2-2得到的状态量代入式(2)，再加上人为设定的量测噪声即得到带有噪声的观测量；

步骤2-4：在每一次仿真时改变状态噪声和量测噪声，即改变Q和R的大小，得到不同噪声水平下的仿真值；

步骤3：载入步骤2生成的仿真数据，运行EM算法，得到不同噪声背景下的气流角估计值，并将估计值与真实值进行比较，通过蒙特卡洛仿真验证EM算法的可行性；

步骤3-1：将EM算法应用于非线性离散系统状态空间模型中，将状态向量X_N＝{x₁,x₂,…,x_N}作为EM中的不可观测数据，量测向量Z_N＝{z₁,z₂,…,z_N}作为EM中的可观测数据，

为未知统计量；P₁为初始时刻状态量协方差矩阵，

为初始时刻状态量预测均值；

由EM算法，其代价函数为：

其中，

为未知统计量的估计值，

为在给定量测数据Z_N和当前未知统计量估计值

下状态向量X_N的条件概率分布，p_μ(X_N,Z_N)为以μ为模型参数建立的状态向量X_N和量测向量Z_N的联合概率分布；

步骤3-2：在给定参数μ的条件下，系统状态和观测值的对数联合概率分布为：

将式(7)代入式(6)，得到：

其中

步骤3-3：通过最大化成本函数推导出未知统计量的最优估计，在高斯分布的假设条件下，采用三阶球面逼近算法对式(9)-式(11)进行计算，统计量μ有解析解，具体如下：

式中

i＝1,2,…4n

其中，N为样本总数；ξⁱ为依据球容积准则选择的西格玛点；n为状态向量的维数；m为量测量的维数；

步骤4：设置飞行器各项参数即能得到相应的气流角估计值。

Images