CN113919194A - 一种基于滤波误差法的非线性飞行动力学系统辨识方法 - Google Patents

Abstract

本发明提供一种基于滤波误差法的非线性飞行动力学系统辨识方法，包括如下步骤：步骤S100，将飞行动力学方程转换为含有加性过程噪声和测量噪声的非线性动力学系统，待辨识的模型参数θ由系统参数、初始状态和滤波增益参数构成；步骤S200，对于给定的测量噪声协方差矩阵R，采用Gauss‑Newton法最小化负对数似然函数，获得参数θ的最大似然估计，其中状态估计采用线性化Kalman滤波器；步骤S300，估计测量噪声协方差矩阵R；步骤S400，重复执行步骤S200～S300，直至收敛，得到非线性飞行动力学系统辨识结果；步骤S500，计算非线性飞行动力学系统辨识结果的不确定度。本发明适用范围更广、实用性更强。

Description

一种基于滤波误差法的非线性飞行动力学系统辨识方法

技术领域

本发明涉及航天技术领域，具体而言，涉及一种基于滤波误差法的非线性飞行动力学系统辨识方法。

背景技术

对于动力学系统辨识，最大似然法(Maximum Likelihood,ML)具有卓越的理论特性和经实践证明的实用性。根据其对随机噪声的不同处理方法，最大似然法分为三类：只考虑测量噪声的输出误差法(Output Error，OE)、只考虑过程噪声的方程误差法(EquationError，EE)、同时考虑过程噪声和测量噪声的滤波误差法(Filter Error，FE)。其中：

输出误差法是当前飞行器飞行动力学系统辨识领域最广泛使用的方法。其技术方案是：在动力学系统模型中只考虑测量噪声而不考虑过程噪声，待辨识参数由系统参数和初始状态构成，采用松弛迭代法交替估计未知参数和测量噪声协方差矩阵，其中，状态预测直接由状态方程数值积分得到。该方法的主要缺点是：只考虑测量噪声而不考虑过程噪声，当飞行试验遭遇湍流或存在其它过程噪声时，辨识结果散布显著增大，并且根据Cramer-Rao界给出的辨识结果不确定度不可信。

滤波误差法虽然理论提出得很早，但实际应用较少，主要原因是，理论上的滤波误差法需要预知系统过程噪声和测量噪声的统计特性，并且算法复杂，存在收敛性和可辨识性问题。较少的滤波误差法应用也主要局限于线性系统。对于非线性飞行动力学系统，仅R.V.Jategaonkar和E.Plaetschke在文献中提出了一种滤波误差法形式，并将其应用于飞机飞行试验气动参数辨识。其技术方案是，对于含加性过程噪声和测量噪声的中等非线性飞行动力学系统，将过程噪声协方差矩阵作为未知参数，与其他模型参数共同构成待辨识参数向量，采用松弛迭代法交替估计未知参数和测量噪声协方差，其中，状态估计采用一个基于恒定增益矩阵K的近似滤波器，该矩阵通过求解连续时间Riccati方程获得。Riccati方程中状态矩阵A和观测矩阵C采用系统关于初始状态的一阶近似，而不考虑他们随状态的变化。该方法的主要缺点是：1)辨识算法仅适用于稳定飞行器，而不适用于中立稳定和不稳定飞行器。这是因为，在该方法中，滤波增益矩阵K通过求解连续时间Riccati方程来得到。理论已经证明，Riccati方程存在唯一正定解的充要条件是，由状态矩阵A和观测矩阵C构成的线性系统是稳定且可观测的。对于中立稳定和不稳定飞行器，不满足稳定性条件。2)辨识算法仅适用于弱到中等非线性动力学系统，而不适用于强非线性动力学系统。这是因为，Riccati方程中状态矩阵A和观测矩阵C采用初始状态近似计算，而不考虑他们随状态的变化。对于强非线性动力学系统，状态矩阵A和观测矩阵C随状态变化较大，由上述Riccati方程计算的滤波增益矩阵K不能反映系统的真实情况，从而导致辨识结果偏差，甚至可能导致辨识算法不收敛。

发明内容

本发明旨在提供一种基于滤波误差法的非线性飞行动力学系统辨识方法，以解决现有滤波误差法因求解Riccati方程带来的适用范围受到极大限制的的问题。

本发明提供的一种基于滤波误差法的非线性飞行动力学系统辨识方法，包括如下步骤：

步骤S100，将飞行动力学方程转换为含有加性过程噪声和测量噪声的非线性动力学系统，构建待辨识的模型参数θ，所述待辨识的模型参数θ由系统参数、初始状态和滤波增益参数构成；

步骤S200，对于给定的测量噪声协方差矩阵R，采用Gauss-Newton法最小化负对数似然函数，获得所述待辨识的模型参数θ的最大似然估计；

步骤S300，估计测量噪声协方差矩阵R；

步骤S400，重复步骤S200和S300直至收敛，得到非线性飞行动力学系统辨识结果。

进一步的，步骤S100包括如下子步骤：

步骤S101，将飞行动力学方程写成如下含有加性过程噪声和测量噪声的非线性动力学系统：

x(t)＝f[x(t)，u(t)，β]+w(t) (1)

y(t)＝g[x(t)，u(t)，β] (2)

z(t_k)＝y(t_k)+v(t_k)，k＝1，2，…，N (3)

式中，x为n维状态向量；y为m维观测向量；u为l维控制向量；β为p维待辨识的系统参数；z为观测向量在N个离散时间点上的测量值，t_k为第k个离散时间点；系统函数f和g分别为n维和m维非线性实数值向量函数；w为过程噪声、v为测量噪声；

步骤S102，构建待辨识的模型参数θ，表示为：

式中，β为系统参数；x₀为非线性飞行动力学系统初始状态；λ为滤波增益参数；上标“T”表示向量或矩阵的转置。

进一步的，步骤S200包括如下子步骤：

步骤S201，采用负对数似然函数作为非线性飞行动力学系统辨识的判据函数J(θ)：

式中，R为测量噪声防方差矩阵；

步骤S202，采用一阶近似的线性化Kalman滤波器对非线性飞行动力学系统进行状态估计：

状态预测方程为：

状态校正方程为：

观测方程为：

式中，K为Kalman滤波器的滤波增益矩阵；

步骤S203，采用中心差分计算观测量对待辨识的模型参数θ的灵敏度：

其中，δθ_j为待辨识的模型参数θ中的未知变量θ_j上的小扰动；

和

分别为在小扰动+δθ_j和-δθ_j作用下的响应变量。对于待辨识的模型参数θ中的每一个未知变量上的小扰动±δθ_j，可以采用与步骤S202中的状态估计类似的扰动系统方程计算得到扰动响应变量

和

步骤S204，采用Gauss-Newton法修正待辨识的模型参数θ：

式中，上标“(l)”表示第l步迭代，

表示第l步迭代的非线性飞行动力学系统辨识结果，Δθ表示待辨识的模型参数θ的修正量；

步骤S205，重复执行步骤S202～S204直至收敛，获得待辨识的模型参数θ的最大似然估计。

进一步的，步骤S300中的测量噪声协方差矩阵R采用下式进行估计：

进一步的，所述基于滤波误差法的非线性飞行动力学系统辨识方法还包括：

步骤S500，计算非线性飞行动力学系统辨识结果的不确定度。

进一步的，所述计算非线性飞行动力学系统辨识结果的不确定度的方法包括如下子步骤：

步骤S501，在有色残差下，非线性飞行动力学系统参数辨识结果

的修正协方差矩阵为：

式中，M为信息矩阵；S(i)为观测灵敏度；R_υυ(i-j)为输出残差的自相关函数；υ(i)为输出残差；分别采用下式计算：

υ(i)＝z(t_i)-y(t_i) (17)

步骤S502，非线性飞行动力学系统参数辨识结果

的标准差

为修正协方差矩阵

中对角线元素的平方根：

步骤S503，置信水平0.95的非线性飞行动力学系统辨识结果

的不确定度为：

式中，

为置信水平0.95的非线性飞行动力学系统辨识结果

的不确定度。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

本发明通过含有过程噪声和测量噪声的非线性飞行动力学系统，将滤波增益矩阵中的元素作为未知参数，与系统参数和初始状态共同构成待辨识的模型参数，并采用松弛迭代法交替估计待辨识的模型参数和测量噪声协方差矩阵，从而具有以下优点：(1)适用范围更广、实用性更强。不仅适用于稳定飞行器，而且适用于中立稳定和不稳定飞行器；不仅适用于弱到中等非线性动力学系统，而且适用于强非线性动力学系统。(2)实现更便捷。无需求解Riccati方程，在现有的输出误差算法框架下，只需将滤波增益增补为待辨识的模型参数，并在状态预测中增加校正步骤即可实现。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的流程图。

图2为本发明实施例的基于滤波误差法的非线性飞行动力学系统辨识方法应用的详细流程图。

图3为示例一中的有噪声和无噪声下的飞机响应展示图。

图4为示例一种的Monte-Carlo仿真辨识结果分布图。

图5a为示例一中的输出误差法的仿真辨识拟合结果展示图。

图5b为示例一中的本发明基于滤波误差法的非线性飞行动力学系统辨识方法的仿真辨识拟合结果展示图。

图6为示例二中的一次偏航激励机动中的参数展示图。

图7为示例二中的模型飞行试验气动参数辨识结果展示图。

图8a为示例二中机动3的辨识拟合结果展示图。

图8b为示例二中机动11的辨识拟合结果展示图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例

如图1所示，本实施例提出一种基于滤波误差法的非线性飞行动力学系统辨识方法，包括：

步骤S100，将飞行动力学方程转换为含有加性过程噪声和测量噪声的非线性动力学系统，构建待辨识的模型参数θ，所述待辨识的模型参数θ由系统参数、初始状态和滤波增益参数构成。具体包括如下子步骤：

步骤S101，将飞行动力学方程写成如下含有加性过程噪声和测量噪声的非线性动力学系统：

x(t)＝f[x(t)，u(t)，β]+w(t) (1)

y(t)＝g[x(t)，u(t)，β] (2)

z(t_k)＝y(t_k)+v(t_k)，k＝1，2，…，N (3)

式中，x为n维状态向量；y为m维观测向量；u为l维控制向量；β为p维待辨识的系统参数；z为观测向量在N个离散时间点上的测量值，t_k为第k个离散时间点；系统函数f和g分别为n维和m维非线性实数值向量函数；w为过程噪声、v为测量噪声。

步骤S102，对于飞行器飞行试验，用于气动参数辨识的激励机动时间通常较短，因此可以合理地假设滤波增益矩阵K是不变的。本步骤将滤波增益矩阵K中的元素(即滤波增益参数)记为向量λ，于是，可以构建待辨识的模型参数θ，表示为：

式中，x₀为非线性飞行动力学系统初始状态；上标“T”表示向量或矩阵的转置，即β^T为系统参数β的转置、

为的初始状态x₀转置、λ^T为滤波增益参数λ的转置。换言之，所述待辨识的模型参数θ由系统参数β、初始状态x₀、滤波增益参数λ构成。

步骤S200，对于给定的测量噪声协方差矩阵R，采用Gauss-Newton法最小化负对数似然函数，获得所述待辨识的模型参数θ的最大似然估计。具体包括如下子步骤：

步骤S201，采用负对数似然函数作为非线性飞行动力学系统辨识的判据函数J(θ)：

式中，R为测量噪声协方差矩阵；

步骤S202，由于需要考虑过程噪声，因此必须采用合适的状态估计器。对于线性系统，Kalman滤波器是一种最优的状态估计器。对于本实施例在步骤S100中提出的非线性动力学系统(1)～(3)，则采用一阶近似的线性化Kalman滤波器对非线性飞行动力学系统进行状态估计：

状态预测方程为：

状态校正方程为：

观测方程为：

式中，K为Kalman滤波器的滤波增益矩阵；

步骤S203，观测灵敏度

可以使用解析微分或有限差分近似进行计算。对于线性系统，通常采用解析微分法来计算灵敏度。然而，这需要通过对动力学系统中方程(1)和(2)求偏微分得到灵敏度方程的显式解，此外，还需要滤波增益矩阵的梯度。对于非线性动力学系统，针对每个新的非线性模型形式，都需要根据系统模型修改计算程序，这是十分繁琐甚至不可行的。而通过灵敏度的有限差分近似，则可以消除这些实际困难。因此本发明采用中心差分来计算观测量对待辨识的模型参数θ的灵敏度：

其中，δθ_j为待辨识的模型参数θ中的未知变量θ_j上的小扰动；

和

分别为在小扰动+δθ_j和-δθ_j作用下的响应变量。对于待辨识的模型参数θ中的每一个未知变量上的小扰动±δθ_j，采用与步骤S202中的状态估计(方程(6)～(8))类似的扰动系统方程计算得到扰动响应变量

和

对于系统参数β的每一个元素，扰动状态变量和扰动观测变量可以按照下式进行演化：

δβ为系统参数β上的小扰动，式(10)中扰动状态变量预测值

通过数值积分进行计算，式(11)中扰动状态变量校正值

和式(12)中扰动观测变量y_p的计算非常简单。对于待辨识的模型参数θ的其他元素，即初始状态x₀、滤波增益参数λ，考虑相应的扰动，采用与方程(11)～(12)类似的方法进行演化。

因此，将有限差分近似推广到滤波算法中时，针对待辨识的模型参数θ的每个元素，需要对扰动状态方程进行数值积分。这一变化带来的计算量的增加不会很大，即使是使用PC机，通常也能在数十秒内完成这些计算。

步骤S204，采用Gauss-Newton法修正待辨识的模型参数θ：

式中，上标“(l)”表示第l步迭代，

表示第l步迭代的非线性飞行动力学系统辨识结果，Δθ表示待辨识的模型参数θ的修正量；

步骤S205，重复执行步骤S202～S204直至收敛，获得待辨识的模型参数θ的最大似然估计。换言之，在给定测量噪声协方差矩阵R的情况下，采用Gauss-Newton法优化判据函数(5)，辨识得到模型参数θ的最大似然估计。

步骤S300，估计测量噪声协方差矩阵R，具体可以采用下式对测量噪声协方差矩阵R进行估计：

步骤S400，重复执行步骤S200～S300直至收敛，得到非线性飞行动力学系统辨识结果。换言之，待辨识的模型参数θ和测量噪声协方差矩阵R的估计采用两步松弛迭代技术，第一步，在给定测量噪声协方差矩阵R的情况下，采用Gauss-Newton法优化判据函数(5)，辨识出模型参数θ，第二步，由式(15)估计测量噪声协方差矩阵R。两步估计交替进行，直至收敛即可得到非线性飞行动力学系统辨识结果。

进一步的，本实施例还提出了辨识准度估计方法，即上述的基于滤波误差法的非线性飞行动力学系统辨识方法还包括：

步骤S500，计算非线性飞行动力学系统辨识结果的不确定度。具体包括如下子步骤：

步骤S501，在有色残差下，非线性飞行动力学系统参数辨识结果

的修正协方差矩阵为：

式中，M为信息矩阵；S(i)为观测灵敏度；R_υυ(i-j)为输出残差的自相关函数；υ(i)为输出残差。分别采用下式计算：

υ(i)＝z(t_i)-y(t_i) (20)

步骤S502，非线性飞行动力学系统参数辨识结果

的标准差

为修正协方差矩阵

中对角线元素的平方根：

步骤S503，置信水平0.95的非线性飞行动力学系统辨识结果

的不确定度为：

式中，

为置信水平0.95的非线性飞行动力学系统辨识结果

的不确定度。需要说明的是，本实施例只是优选置信水平为0.95，当置信水平为其他值时，可以相应地调整式(21)的系数。

以下通过对上述的基于滤波误差法的非线性飞行动力学系统辨识方法的步骤S100～S500的具体应用进一步详述：

对于飞行器非线性动力学系统辨识，采用刚体六自由度动力学方程作为状态方程，其中线位移动力学方程建立在速度系下，角位移动力学方程通常建立在体轴系下。这种处理方法的优点是，所有状态变量均有测量数据，换言之，状态变量集是观测变量集的一个子集，滤波增益矩阵K的主元与次元分明，在大多数情况下可以只辨识主元，而将次元设置为0。

对于完整的六自由度动力学方程，状态向量x^T＝(V，α，β，p，q，r，θ，ψ，φ)，观测向量y^T＝(V，α，β，p，q，r，a_x，a_y，a_z，θ，ψ，φ)。

状态方程：

式中，V为飞行速度；α为迎角；β为侧滑角；p，q，r为体轴系角速率分量；θ，ψ，φ为俯仰、偏航、滚转姿态角；

为动压；ρ为大气密度；c，l分别为纵向和横侧向参考长度；S为参考面积；m为飞行器质量；J_x，J_y，J_z，J_xz为相对于体轴的转动惯量和惯性积；

C_D，C_L，C_C为速度系气动力系数；C_l，C_m，C_n为体轴系气动力矩系数；F_e为发动机推力，假设沿飞行器纵轴；w_i为过程噪声。

观测方程：

V_m＝V+v_V (24a)

α_m＝α+v_α (24b)

β_m＝β+v_β (24c)

p_m＝p+v_p (24d)

q_m＝q+v_q (24e)

r_m＝r+v_r (24f)

θ_m＝θ+v_θ (24j)

ψ_m＝ψ+v_ψ (24k)

φ_m＝φ+v_φ (24l)

式中，下标“m”表示测量值；n_x，n_y，n_z为体轴系过载(即用重力加速度g无因次化的线加速度，其中n_x沿体轴x_b为正，n_y沿体轴y_b为正，n_z沿体轴z_b反方向为正)；C_A，C_N，C_Y为体轴系气动力系数；v_i为观测噪声。

速度系与体轴系气动力系数的关系式为：

C_D＝C_A cosαcosβ-C_Ysinβ+C_Nsinαcosβ (25a)

C_L＝-C_A sinα+C_Ncosα (25b)

C_C＝C_A cosαsinβ+C_Ycosβ+C_Nsinαsinβ (25c)

体轴系气动力和力矩系数采用线性模型：

式中，δ_e，δ_r，δ_a分别为升降舵、方向舵、副翼偏角；M为马赫数；a为声速。

对于式(23)～(26)所描述的非线性飞行动力学系统，应用上述的基于滤波误差法的非线性飞行动力学系统辨识方法的步骤S100～S500，就可以进行气动参数辨识，并给出辨识结果的不确定度。如图2所示，具体如下：

Box 1：输入初始数据。包括：1)输入飞行器主要物理、几何参数；2)输入飞行试验实测数据；3)指定状态变量和观测变量；4)指定待辨识参数。

Box 2：确定待辨识模型参数的初值θ⁽⁰⁾。包括：气动参数(对应系统参数向量)、初始状态、滤波增益矩阵的辨识起始值。其中，气动参数取地面预测结果或方程误差法辨识结果，初始状态取实测结果，滤波增益矩阵主元(K_V，V，K_α，α，K_β，β，K_p，p，K_q，q，K_r，r，K_θ，θ，K_ψ，ψ，K_φ，φ)取[0，1]区间的小正值(例如0.2)，其他元素取为0。

Box 3：初步估计测量噪声协方差矩阵R。采用现有方法(如Morelli E.Estimatingnoise characteristics from flight test data using optimal Fouriersmoothing.Journal ofAircraft，1995，32(4)：689-695)先初步估计测量噪声协方差矩阵。

Box 4：计算判据函数初值J⁽⁰⁾。利用松弛迭代得到的θ和R，由式(6)～(8)进行状态预测和校正以及观测量预测，得到y(t_k)，再由式(5)计算J⁽⁰⁾。

Box 5：状态预测与校正，计算观测量y(t_k)。利用第l-1次迭代得到的参数θ^(l-1)和R，由式(6)～(8)进行状态预测和校正以及观测量预测，得到y(t_k)。

Box 6：计算观测灵敏度

采用中心差分法，由式(9)计算观测灵敏度

其中扰动量δθ_j取为max(10^-6|θ_j|，10^-8)。

Box 7：计算待辨识参数修正量Δθ。利用松弛迭代得到的R、观测量实测值z(t_k)、Box 5计算得到的y(t_k)、以及Box 6计算得到的

由式(14)计算待辨识参数修正量Δθ。

Box 8：计算判据函数J^(l)。利用松弛迭代得到的R、观测量实测值z(t_k)、以及Box 5计算得到的y(t_k)，由式(5)计算判据函数J^(l)。若J^(l)＜J^(l-1)且|[J^(l)-J^(l-1)]/J^(l-1)|＜ε，则Gauss-Newton迭代收敛。ε为设定阈值，根据需求设定。

Box 9：修正Δθ。若J^(l)＞J^(l-1)，将减小步长，重新计算J^(l)，直至J^(l)＜J^(l-1)。

Box 10：估计测量噪声协方差矩阵R。Gauss-Newton迭代收敛后，比较当前准则函数值J^(l)与前一步松弛迭代的判据函数J_pre，若|[J^(l)-J_pre]/J_pre|＜ε，则松弛迭代收敛，否则，由式(15)估计测量噪声协方差矩阵R，并令J_pre＝J^(l)，进行下一次Gauss-Newton迭代。

Box 11：计算参数估计的不确定度

利用松弛迭代得到的R、观测量实测值z(t_k)、Box 5计算得到的y(t_k)、以及Box 6计算得到的

由式(16)～(22)计算参数估计的不确定度

Box 12：输出辨识结果。包括但不限于：1)输出待辨识模型参数的估计值

及其不确定度

2)输出观测量预测值y(t_k)；3)输出观测噪声协方差矩阵R。

在实际的气动参数辨识飞行试验中，纵向激励和横侧向激励往往是分别施加的，此时需要根据激励机动的具体情况，选择部分变量作为状态量和观测量，相应地选择部分气动参数、初始状态、滤波增益矩阵元素作为待辨识参数。

以下通过具体示例对本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的步骤S100～S500的有益效果进行说明。

示例一：飞行仿真气动参数辨识

采用VFW-614ATTAS纵向飞行仿真数据，研究本发明的基于滤波误差法的非线性飞行动力学系统辨识方法在受控环境中的性能。飞行仿真初始状态为：高度h₀＝2517.2m，速度V₀＝130m/s，迎角α₀＝1.9°，俯仰角速率q₀＝0°/s，俯仰角θ₀＝1.75°。

仿真采用的气动力模型为：

气动参数基准值采用飞行试验数据辨识结果。

采用升降舵“3211”输入，以激发飞机的纵向运动模态。仿真计算过程中加入过程噪声，在仿真输出上加载测量噪声，过程噪声和测量噪声均为Gauss白噪声，其标准差列于表1。采样时间步长为40ms。

表1，飞行仿真过程噪声和测量噪声标准差：

共进行300次仿真，图3给出了仿真所采用的升降舵偏输入以及其中一次的仿真结果与无噪声结果的比较。可见，在过程噪声的影响下，飞行状态参数存在偏离无噪声结果的情况。

针对300次仿真结果，分别采用输出误差法和本发明的基于滤波误差法的非线性飞行动力学系统辨识方法开展气动参数辨识。辨识模型的状态方程为：

观测方程为：

V_m＝V+v_V (29a)

α_m＝α+v_α (29b)

q_m＝q+v_r (29c)

θ_m＝θ+v_θ (29e)

气动力模型为：

速度系气动力系数(C_D，C_L)与体轴系气动力系数(C_A，C_N)的关系式为：

C_D＝C_A cos+C_N sinα (31a)

C_L＝-C_A sinα+C_Ncosα (31b)

待辨识的模型参数为：

图4显示了俯仰控制导数

和俯仰阻尼导数C_mq辨识结果的分布情况，可见Monte-Carlo仿真辨识结果近似服从正态分布。从两种方法对比来看，辨识结果均值相近，但本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的散布范围明显较小。

表2列出了Monte-Carlo仿真辨识结果的统计特性，其中“均值”表示参数估计值的统计平均值，“2σ”表示参数估计值统计标准差的2倍，“平均U₉₅”表示不确定度U₉₅的统计平均值。表中还给出了气动参数的基准值。可见，两种方法辨识结果的均值都接近于基准值，但本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的2σ显著小于输出误差法，这与图2所示分布情况是一致的。本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的平均U₉₅也小于输出误差法。

表2，Monte-Carlo仿真辨识结果统计特性：

图5a、图5b给出了其中一次的仿真辨识拟合结果，即辨识模型输出与仿真响应数据的比较。本发明的本发明的基于滤波误差法的非线性飞行动力学系统辨识方法由于增加了状态校正，辨识拟合结果明显好于输出误差法。

上述仿真辨识结果表明，在受过程噪声影响的情况下，本发明的本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的辨识结果比输出误差法更可信。

示例二：

中国空气动力研究与发展中心利用缩比模型飞行试验，开展了一飞机布局的气动参数辨识飞行试验。飞行高度为1500～2000m，飞行速度为45～60m/s。测量参数包括体轴系线加速度和角速率、姿态角、操纵面偏转角、GPS位置和速度、气流角、静压、动压以及发动机转子的转速等，采样频率为100Hz。

为了辨识气动稳定和控制导数，分别在俯仰、偏航、滚转通道加载连续三个偶极方波激励，以激发试验模型的相关运动模态。进行了4个架次的飞行试验，共39个激励机动。从飞行试验数据来看，试验模型存在程度不同的持续“颠簸”现象，初步判断飞行空域有轻度到中度湍流。

本示例使用13个偏航激励机动飞行数据辨识航向气动参数。辨识模型的状态方程为：

观测方程为：

β_m＝β+v_β (34a)

r_m＝r+v_r (34b)

气动力模型为：

速度系气动力系数C_C与体轴系气动力系数C_Y的关系式为：

C_C＝C_A cosαsinβ+C_Y cosβ+C_N sinαsinβ (36)

待辨识参数为：

由于大气湍流和阵风等的影响，飞行试验中三通道之间的耦合较严重。图6给出了一次偏航激励机动中的舵偏角和角速率，可见机动过程中，湍流和阵风等产生俯仰和滚转扰动，飞控系统需要偏转升降舵和副翼以抑制扰动。因此，在辨识模型中需要考虑俯仰和滚转对偏航运动的影响，一种简单而实用的方法是，在方程(32)～(35)中，除状态参数β和r外，其他参数(V，α，p，q，θ，φ)取为实测数据。

分别采用输出误差法和本发明的基于滤波误差法的非线性飞行动力学系统辨识方法进行气动参数辨识。图7给出了主要气动参数

的辨识结果及其置信水平0.95不确定度。可见，本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的辨识结果的散布范围总体上小于输出误差法。图中还给出了本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的辨识结果的均值，用虚线表示。值得注意的是，对于本发明的基于滤波误差法的非线性飞行动力学系统辨识方法，辨识结果均值基本位于各次机动辨识结果的不确定度范围内，表明不确定度估计结果是合理的。很显然，输出误差法不具备这一特性。尽管输出误差法不确定度估计也是采用有色残差情况下的修正协方差方法，但由于辨识模型对实测数据的拟合较差，严重背离了参数辨识方法的适用条件，导致不确定度估计结果不可信。

图8a、图8b给出了机动3和机动11的辨识拟合结果。本发明的基于滤波误差法的非线性飞行动力学系统辨识方法对测量数据的拟合都很好。从输出误差法的辨识拟合结果来看，机动3的拟合较差，相应地，辨识结果的不确定度较大，而机动11的拟合较好，辨识结果不确定度也相应地较小，初步判断这是由不同湍流强度造成的。

表3列出了滤波增益矩阵主元K_β，β和K_r，r的辨识结果及其不确定度，其他元素取为0。与估计值相比，不确定度较小，说明判据函数对滤波增益较敏感，辨识结果是可靠的。

表3，飞行试验滤波增益辨识结果：

上述辨识结果和分析表明，对于湍流中的飞行试验，本发明的基于滤波误差法的非线性飞行动力学系统辨识方法的辨识结果优于输出误差法。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于滤波误差法的非线性飞行动力学系统辨识方法，其特征在于，包括如下步骤：

步骤S100，将飞行动力学方程转换为含有加性过程噪声和测量噪声的非线性动力学系统，构建待辨识的模型参数θ，所述待辨识的模型参数θ由系统参数、初始状态和滤波增益参数构成；

步骤S200，对于给定的测量噪声协方差矩阵R，采用Gauss-Newton法最小化负对数似然函数，获得所述待辨识的模型参数θ的最大似然估计；

步骤S300，估计测量噪声协方差矩阵R；

步骤S400，重复执行步骤S200～S300，直至收敛，得到非线性飞行动力学系统辨识结果。

2.根据权利要求1所述的基于滤波误差法的非线性飞行动力学系统辨识方法，其特征在于，步骤S100包括如下子步骤：

步骤S101，将飞行动力学方程写成如下含有加性过程噪声和测量噪声的非线性动力学系统：

x(t)＝f[x(t)，u(t)，β]+w(t) (1)

y(t)＝g[x(t)，u(t)，β] (2)

z(t_k)＝y(t_k)+v(t_k)，k＝1，2，…，N (3)

式中，x为n维状态向量；y为m维观测向量；u为l维控制向量；β为p维待辨识的系统参数；z为观测向量在N个离散时间点上的测量值，t_k为第k个离散时间点；系统函数f和g分别为n维和m维非线性实数值向量函数；w为过程噪声、v为测量噪声；

步骤S102，构建待辨识的模型参数θ，表示为：

式中，β为系统参数；x₀为非线性飞行动力学系统初始状态；λ为滤波增益参数；上标“T”表示向量或矩阵的转置。

3.根据权利要求2所述的基于滤波误差法的非线性飞行动力学系统辨识方法，其特征在于，步骤S200中包括如下子步骤：

步骤S201，采用负对数似然函数作为非线性飞行动力学系统辨识的判据函数J(θ)：

式中，R为测量噪声协方差矩阵；

步骤S202，采用一阶近似的线性化Kalman滤波器对非线性飞行动力学系统进行状态估计：

状态预测方程为：

状态校正方程为：

观测方程为：

式中，K为Kalman滤波器的滤波增益矩阵；

步骤S203，采用中心差分计算观测量对待辨识的模型参数θ的灵敏度：

其中，δθ_j为待辨识的模型参数θ中的未知变量θ_j上的小扰动；y_i，p+和y_i，p-分别为在小扰动+δθ_j和-8θ_j作用下的响应变量；

步骤S204，采用Gauss-Newton法修正待辨识的模型参数θ：

式中，上标“(l)”表示第l步迭代，

表示第l步迭代的非线性飞行动力学系统辨识结果，Δθ表示待辨识的模型参数θ的修正量；

步骤S205，重复执行步骤S202～S204直至收敛，获得待辨识的模型参数θ的最大似然估计。

4.根据权利要求3所述的基于滤波误差法的非线性飞行动力学系统辨识方法，其特征在于，步骤S300中的测量噪声协方差矩阵R采用下式进行估计：

5.根据权利要求1-4任一项所述的基于滤波误差法的非线性飞行动力学系统辨识方法，其特征在于，所述基于滤波误差法的非线性飞行动力学系统辨识方法还包括：

步骤S500，计算非线性飞行动力学系统辨识结果的不确定度。

6.根据权利要求5所述的基于滤波误差法的非线性飞行动力学系统辨识方法，其特征在于，所述计算非线性飞行动力学系统辨识结果的不确定度的方法包括如下子步骤：

步骤S501，在有色残差下，非线性飞行动力学系统参数辨识结果

的修正协方差矩阵为：

式中，M为信息矩阵；S(i)为观测灵敏度；R_υυ(i-j)为输出残差的自相关函数；υ(i)为输出残差；分别采用下式计算：

υ(i)＝z(t_i)-y(t_i) (17)

步骤S502，非线性飞行动力学系统参数辨识结果

的标准差

为修正协方差矩阵

中对角线元素的平方根：

步骤S503，置信水平0.95的非线性飞行动力学系统辨识结果

的不确定度为：

式中，

为置信水平0.95的非线性飞行动力学系统辨识结果

的不确定度。

Images