CN115169002A - 直气复合飞行控制系统时变参数辨识方法 - Google Patents

Abstract

直气复合飞行控制系统时变参数辨识方法，它属于飞行器飞行参数辨识领域。本发明解决了采用小扰动线性化方法得到的常值气动参数无法满足飞行器设计要求的问题。本发明方法采取的主要技术方案为：建立飞行器体坐标系和飞行器速度坐标系，基于坐标系建立俯仰通道姿态动力学方程，基于俯仰通道姿态动力学方程设计时变参数滤波器；再基于设计的滤波器系统对气动参数进行在线辨识。本发明依据魏尔斯特拉斯逼近定理，通过多项式拟合，将时变参数估计转化为非时变参数估计。本发明方法可以应用于飞行器飞行参数辨识。

Description

直气复合飞行控制系统时变参数辨识方法

技术领域

本发明属于飞行器飞行参数辨识领域，具体涉及一种直气复合飞行控制系统时变参数辨识方法。

背景技术

飞行器模型参数辨识是设计和研制新型号飞行器的重要基础和依据，建立精确的飞行器动力学模型是飞行器控制系统设计、飞行品质评价、飞行器自适应及自学习控制以及地面仿真的重要一环。随着需求的提升和飞行器气动设计的发展，飞行器的机动能力在不断提升，经常会做出大攻角，高过载的机动动作。尤其是采用直气复合控制策略的飞行器，当直接力参与时会对飞行器运动状态产生较大影响，随之而来的飞行器气动特性的非线性特征愈发明显，由小扰动线性化方法得到的常值气动参数已经无法满足飞行器的设计要求。在其时变气动参数的整个变化范围内，如何对其进行有效辨识已经成为飞行控制器设计人员迫切需要解决的一个问题。

发明内容

本发明的目的是为解决采用小扰动线性化方法得到的常值气动参数无法满足飞行器设计要求的问题，而提出的一种直气复合飞行控制系统时变参数辨识方法。

本发明为解决上述技术问题所采取的技术方案是：

直气复合飞行控制系统时变参数辨识方法，所述方法具体包括以下步骤：

步骤一、直气复合飞行控制时变参数滤波器系统设计

步骤一一、建立飞行器体坐标系ox₁y₁z₁和飞行器速度坐标系ox₃y₃z₃；

步骤一二、基于步骤一一的坐标系建立俯仰通道姿态动力学方程；所述俯仰通道姿态动力学方程为：

其中，ω_z为俯仰角速率，

为ω_z的一阶导数，δ_z为舵偏角，F_z为直接力，α为攻角，气动参数a₁、a₂和a₃分别定义为：

直接力参数l_z定义为：l_z＝-l/J_z，J_z为飞行器绕坐标轴oz₃的转动惯量，

为俯仰力矩对俯仰角速率的偏导数，

为俯仰力矩对攻角的偏导数，

为俯仰力矩对舵偏角的偏导数，l是侧推发动机到飞行器质心的距离；

步骤一三、基于步骤一二的俯仰通道姿态动力学方程设计时变参数滤波器；

对于气动参数a₂，利用关于α的一次多项式近似a₂(α)，即：

a₂(α)≈a₂₀+k₂₁α (3)

其中，a₂₀和k₂₁为常数；

将式(3)代入俯仰通道姿态动力学方程(1)得到：

滤波器系统状态变量x₁、x₂、x₃、x₄、x₅和x₆分别定义为：x₁＝ω_z，x₂＝a₁，x₃＝a₂₀，x₄＝k₂₁，x₅＝l_z，x₆＝a₃；

则滤波器系统状态向量X为：

X＝[x₁ x₂ x₃ x₄ x₅ x₆]^T (5)

定义滤波器系统状态方程为：

其中，

为X的一阶导数，f(X)为滤波器系统状态转移函数；

滤波器系统量测方程h为：

h＝[1 0 0 0 0 0] (8)

步骤二、基于步骤一中设计的滤波器系统对a₂₀和k₂₁进行在线辨识；

同理，对气动参数a₁和a₃进行辨识。

进一步地，所述步骤一一中，建立飞行器体坐标系ox₁y₁z₁和飞行器速度坐标系ox₃y₃z₃的具体过程为：

以飞行器质心为原点o，以飞行器的中心对称轴为ox₁轴，以飞行器速度方向为ox₃轴，oy₁轴和oy₃轴均在飞行器的铅垂面上，且oy₁轴与oy₃轴的关系通过攻角和偏航角决定，oz₁轴与ox₁轴、oy₁轴满足右手定则，oz₃轴与ox₃轴、oy₃轴满足右手定则。

进一步地，所述步骤二的具体过程为：

对式(7)进行离散化，得到离散系统：

其中，f(X_k)为第k步的滤波器系统状态转移函数，T为采样间隔，X_k为第k步的滤波器系统状态向量；

再基于离散系统对参数a₂₀和k₂₁进行在线辨识。

进一步地，所述对式(7)进行离散化采用的是差分法。

进一步地，所述基于离散系统对参数a₂₀和k₂₁进行在线辨识，采用的是EKF算法。

更进一步地，所述EKF算法的流程为：

步骤1、预测状态和误差协方差矩阵：

其中，

为状态的一步预测结果，

为第k-1步的滤波器系统状态转移函数估计结果，P_k|k-1为误差协方差矩阵的一步预测结果，P_k-1|k-1为第k-1步的误差协方差矩阵估计结果，F_k为第k步的雅克比矩阵，上角标T代表矩阵的转置，Q_k为第k步的过程噪声；

步骤2、计算卡尔曼增益：

其中，z_k为第k步的实际测量值，

是使用

计算出来的预测测量值，

为残差，

为第k步的状态转移矩阵，R_k为第k步的量测噪声协方差矩阵，S_k为中间变量矩阵，上角标-1代表矩阵的逆，K_k为第k步的增益矩阵；

步骤3、更新状态和误差协方差矩阵：

P_k|k＝(I-K_kH_k)P_k|k-1 (16)

其中，

为第k步的状态估计结果，P_k|k为第k步的误差协方差矩阵估计结果，I为单位矩阵。

本发明的有益效果是：

本发明针对带有时变参数的系统提出了一种滤波器建模新方法，依据魏尔斯特拉斯逼近定理，通过多项式拟合，将时变参数估计问题转化为非时变参数估计问题，并在仿真中引入直接力作用，通过EKF联合滤波，在线辨识时变气动参数与直接力参数。从仿真结果看出，本发明方法能有效地估计时变气动参数与直接力参数，可以满足飞行器设计要求。

附图说明

图1是飞行器速度坐标系与飞行器体坐标系关系的示意图；

图中，M_z是俯仰力矩，β是侧滑角；

图2是飞行器攻角变化情况图；

图3是侧推发动机点火情况图；

图4是俯仰角速率ω_z估计值、真实值和测量值的对比图；

图5是a₂估计值与真实值对比图；

图6是对参数a₂₀的估计结果图；

图7是对参数k₂₁的估计结果图；

图8是对a₁的估计结果图；

图9是对a₃的估计结果图；

图10是对l_z的估计结果图；

图11是舵偏角变化情况图；

图12是俯仰角与俯仰角指令变化情况图。

具体实施方式

具体实施方式一、本实施方式所述的直气复合飞行控制系统时变参数辨识方法，所述方法具体包括以下步骤：

步骤一、直气复合飞行控制时变参数滤波器系统设计

步骤一一、建立飞行器体坐标系ox₁y₁z₁和飞行器速度坐标系ox₃y₃z₃；

步骤一二、基于步骤一一的坐标系建立俯仰通道姿态动力学方程；考虑气动非线性特性的情况下，所述俯仰通道姿态动力学方程为：

其中，ω_z为俯仰角速率，俯仰角速率的估计值、真实值和测量值的对比情况如图4所示，

为ω_z的一阶导数，δ_z为舵偏角，舵偏角的变化情况如图11所示，F_z为直接力，α为攻角，气动参数a₁、a₂和a₃分别定义为：

直接力参数l_z定义为：l_z＝-l/J_z，J_z为飞行器绕坐标轴oz₃的转动惯量，

为俯仰力矩对俯仰角速率的偏导数，

为俯仰力矩对攻角的偏导数，

为俯仰力矩对舵偏角的偏导数，l是侧推发动机到飞行器质心的距离；

步骤一三、基于步骤一二的俯仰通道姿态动力学方程设计时变参数滤波器；

以攻角气动参数a₂为例，考虑a₂为关于攻角α的函数，即a₂＝a₂(α)。依据魏尔斯特拉斯逼近定理第一定律(Stone–Weierstrass Theorem)，对于任何一个连续函数，一定可以找到最高项次数不大于n的多项式，使之与该连续函数差为零。即，闭区间上的连续函数可用多项式级数一致逼近，将该函数通过多项式表示，可写成：

a₂＝a₂(α)＝a₂₀+k₁α+k₂α²+k₃α³+... (2)

其中，a₂₀、k₁、k₂、...为常数；

对于气动参数a₂，当攻角在临界攻角范围内变化时，可以利用关于α的一次多项式近似a₂(α)，即：

a₂(α)≈a₂₀+k₂₁α (3)

其中，a₂₀和k₂₁为常数；

将式(3)代入俯仰通道姿态动力学方程(1)得到：

滤波器系统状态变量x₁、x₂、x₃、x₄、x₅和x₆分别定义为：x₁＝ω_z，x₂＝a₁，x₃＝a₂₀，x₄＝k₂₁，x₅＝l_z，x₆＝a₃；

则滤波器系统状态向量X为：

X＝[x₁ x₂ x₃ x₄ x₅ x₆]^T (5)

定义滤波器系统状态方程为：

其中，

为X的一阶导数，f(X)为滤波器系统状态转移函数；

滤波器系统量测方程h为：

h＝[1 0 0 0 0 0] (8)

步骤二、基于步骤一中设计的滤波器系统对a₂₀和k₂₁进行在线辨识；

同理，对气动参数a₁和a₃进行辨识(辨识方法同气动参数a₂)。

对于直气复合控制飞行器，飞行器飞行状态由于侧推发动机的点火会产生剧烈变化，其气动参数不再像传统飞行器一样可以考虑为常数，而是时变的，时变气动参辨识问题一直是设计飞行器控制器的主要障碍。本发明针对直气复合飞行控制系统时变参数估计问题，提出一种新型的工程解决算法，通过多项式拟合，将时变参数拟合为若干由定常参数构成的多项式进行辨识，从而将时变参数辨识问题转化为定常参数辨识问题，降低实现难度。本发明的飞行器可以为导弹。

具体实施方式二：结合图1说明本实施方式。本实施方式与具体实施方式一不同的是，所述步骤一一中，建立飞行器体坐标系ox₁y₁z₁和飞行器速度坐标系ox₃y₃z₃的具体过程为：

以飞行器质心为原点o，以飞行器的中心对称轴为ox₁轴，以飞行器速度方向为ox₃轴，oy₁轴和oy₃轴均在飞行器的铅垂面上，且oy₁轴与oy₃轴的关系通过攻角和偏航角决定，oz₁轴与ox₁轴、oy₁轴满足右手定则，oz₃轴与ox₃轴、oy₃轴满足右手定则。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤二的具体过程为：

对式(7)进行离散化，得到离散系统：

其中，f(X_k)为第k步的滤波器系统状态转移函数，T为采样间隔，X_k为第k步的滤波器系统状态向量；

再基于离散系统对参数a₂₀和k₂₁进行在线辨识。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述对式(7)进行离散化采用的是差分法。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述基于离散系统对参数a₂₀和k₂₁进行在线辨识，采用的是EKF算法。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述EKF算法的流程为：

步骤1、预测状态和误差协方差矩阵：

其中，

为状态的一步预测结果，

为第k-1步的滤波器系统状态转移函数估计结果，P_k|k-1为误差协方差矩阵的一步预测结果，P_k-1|k-1为第k-1步的误差协方差矩阵估计结果，F_k为第k步的雅克比矩阵，上角标T代表矩阵的转置，Q_k为第k步的过程噪声；

步骤2、计算卡尔曼增益：

其中，z_k为第k步的实际测量值，

是使用

计算出来的预测测量值，

为残差，

为第k步的状态转移矩阵，R_k为第k步的量测噪声协方差矩阵，S_k为中间变量矩阵，上角标-1代表矩阵的逆，K_k为第k步的增益矩阵；

步骤3、更新状态和误差协方差矩阵：

P_k|k＝(I-K_kH_k)P_k|k-1 (16)

其中，

为第k步的状态估计结果，P_k|k为第k步的误差协方差矩阵估计结果，I为单位矩阵。

其它步骤及参数与具体实施方式一至五之一相同。

本发明对时变参数a₂(α)进行辨识之前，需要分析非线性滤波器可观性，具体过程如下：

由非线性系统弱可观理论可知，系统有6个状态变量，1个可观测变量，则需要求取5阶李导数：

依据上式所构成非线性系统

零阶李导数为：

d₀＝H＝x₁ (17)

一阶李导数为：

二阶李导数为：

三阶李导数为：

四阶李导数为：

五阶李导数为：

令：

则观测矩阵可写做：

经过系统可观性分析，通过本发明的滤波器设计方法，将时变参数a₂(α)的辨识问题转化为两个非时变参数a₂₀和k₂₁的辨识问题，由此可以对时变气动参数进行在线辨识。

仿真部分

本发明kalman滤波通过威尔施特拉斯逼近定理，将时变参数逼近为n阶多项式，从而实现将时变参数估计问题转化为定常参数估计问题，下面通过一个仿真算例讲解本发明具体实施方案：

仿真实验中，设过程噪声矩阵Q＝0.001·I_6*6，测量噪声R＝0.01，真实气动参数a₁为0.075，a₃为24.5，真实直接力参数l_z为0.05。依据以往实验经验，气动参数a₂随攻角α增大而增大，并在攻角从0°变化至30°时，a₂的值增大为初值的二倍，攻角变化如图2所示，故设计a₂₀真实值为10，k₂₁真实值为18。为了更好地体现滤波效果，令飞行器弹俯仰角指令θ_c呈式(25)变化，俯仰角指令变化如图12所示，即：

θ_c＝p₁cos(p₂t)+p₃t (25)

其中，p₁，p₂，p₃为可调节参数。为估计直接力参数l_z且符合飞行器实际飞行工况，对飞行器施加开环脉冲侧向力信号，侧推发动机点火情况如图3所示，其信号幅值为200N，点火持续时间为4ms，两次点火间隔为36ms，脉冲发动机开机时间为仿真开始后1s～5s。

滤波器状态变量定义为x₁＝ω_z，x₂＝a₁，x₃＝a₂₀，x₄＝k₂₁，x₅＝l_z，x₆＝a₃，则系统状态向量为：

X＝[x₁ x₂ x₃ x₄ x₅ x₆]^T (26)

定义状态方程为：

其中，系统状态转移函数为：

量测方程为：

h＝[1 0 0 0 0 0] (29)

依据差分法得到离散系统X_k+1＝f(X_k)；

其中：

则EKF算法流程如下：

1)预测状态和误差协方差矩阵：

2)计算卡尔曼增益：

其中

为状态转移矩阵，R_k为量测噪声协方差矩阵。

3)更新状态和误差协方差矩阵：

P_k|k＝(I-K_kH_k)P_k|k-1 (37)

通过上述滤波器设计方法，将时变参数a₂(α)的辨识问题转化为两个非时变参数a₂₀和k₂₁的辨识问题，由此可以对时变气动参数进行在线辨识。a₂的估计值与真实值的对比情况如图5所示，参数a₂₀的估计结果如图6所示，参数k₂₁的估计结果如图7所示。同理，对a₁、a₃和l_z进行辨识，得到a₁的估计结果如图8所示，得到a₃的估计结果如图9所示，得到l_z的估计结果如图10所示。

本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (6)

1.直气复合飞行控制系统时变参数辨识方法，其特征在于，所述方法具体包括以下步骤：

步骤一、直气复合飞行控制时变参数滤波器系统设计

步骤一一、建立飞行器体坐标系ox₁y₁z₁和飞行器速度坐标系ox₃y₃z₃；

步骤一二、基于步骤一一的坐标系建立俯仰通道姿态动力学方程；所述俯仰通道姿态动力学方程为：

其中，ω_z为俯仰角速率，

为ω_z的一阶导数，δ_z为舵偏角，F_z为直接力，α为攻角，气动参数a₁、a₂和a₃分别定义为：

直接力参数l_z定义为：l_z＝-l/J_z，J_z为飞行器绕坐标轴oz₃的转动惯量，

为俯仰力矩对俯仰角速率的偏导数，

为俯仰力矩对攻角的偏导数，

为俯仰力矩对舵偏角的偏导数，l是侧推发动机到飞行器质心的距离；

步骤一三、基于步骤一二的俯仰通道姿态动力学方程设计时变参数滤波器；

对于气动参数a₂，利用关于α的一次多项式近似a₂(α)，即：

a₂(α)≈a₂₀+k₂₁α (3)其中，a₂₀和k₂₁为常数；

将式(3)代入俯仰通道姿态动力学方程(1)得到：

滤波器系统状态变量x₁、x₂、x₃、x₄、x₅和x₆分别定义为：x₁＝ω_z，x₂＝a₁，x₃＝a₂₀，x₄＝k₂₁，x₅＝l_z，x₆＝a₃；

则滤波器系统状态向量X为：

X＝[x₁ x₂ x₃ x₄ x₅ x₆]^T (5)

定义滤波器系统状态方程为：

其中，

为X的一阶导数，f(X)为滤波器系统状态转移函数；

滤波器系统量测方程h为：

h＝[1 0 0 0 0 0] (8)

步骤二、基于步骤一中设计的滤波器系统对a₂₀和k₂₁进行在线辨识；

同理，对气动参数a₁和a₃进行辨识。

2.根据权利要求1所述的直气复合飞行控制系统时变参数辨识方法，其特征在于，所述步骤一一中，建立飞行器体坐标系ox₁y₁z₁和飞行器速度坐标系ox₃y₃z₃的具体过程为：

以飞行器质心为原点o，以飞行器的中心对称轴为ox₁轴，以飞行器速度方向为ox₃轴，oy₁轴和oy₃轴均在飞行器的铅垂面上，且oy₁轴与oy₃轴的关系通过攻角和偏航角决定，oz₁轴与ox₁轴、oy₁轴满足右手定则，oz₃轴与ox₃轴、oy₃轴满足右手定则。

3.根据权利要求2所述的直气复合飞行控制系统时变参数辨识方法，其特征在于，所述步骤二的具体过程为：

对式(7)进行离散化，得到离散系统：

其中，f(X_k)为第k步的滤波器系统状态转移函数，T为采样间隔，X_k为第k步的滤波器系统状态向量；

再基于离散系统对参数a₂₀和k₂₁进行在线辨识。

4.根据权利要求3所述的直气复合飞行控制系统时变参数辨识方法，其特征在于，所述对式(7)进行离散化采用的是差分法。

5.根据权利要求4所述的直气复合飞行控制系统时变参数辨识方法，其特征在于，所述基于离散系统对参数a₂₀和k₂₁进行在线辨识，采用的是EKF算法。

6.根据权利要求5所述的直气复合飞行控制系统时变参数辨识方法，其特征在于，所述EKF算法的流程为：

步骤1、预测状态和误差协方差矩阵：

其中，

为状态的一步预测结果，

为第k-1步的滤波器系统状态转移函数估计结果，P_kk-1为误差协方差矩阵的一步预测结果，P_k-1k-1为第k-1步的误差协方差矩阵估计结果，F_k为第k步的雅克比矩阵，上角标T代表矩阵的转置，Q_k为第k步的过程噪声；

步骤2、计算卡尔曼增益：

其中，z_k为第k步的实际测量值，

是使用

计算出来的预测测量值，

为残差，

为第k步的状态转移矩阵，R_k为第k步的量测噪声协方差矩阵，S_k为中间变量矩阵，上角标-1代表矩阵的逆，K_k为第k步的增益矩阵；

步骤3、更新状态和误差协方差矩阵：

P_kk＝(I-K_kH_k)P_kk-1 (16)

其中，

为第k步的状态估计结果，P_kk为第k步的误差协方差矩阵估计结果，I为单位矩阵。

Images