面向高超声速飞行器双重不确定性的自适应状态估计

冯肖雪，刘 萌，李笑宇，潘 峰,2

(1. 北京理工大学自动化学院，北京 100081；2. 北京理工大学昆明产业技术研究院，昆明 650106)

0 引 言

高超声速飞行器一般指临近空间内飞行速度超过5倍声速的飞行器，具有飞行速度快、反应迅速、突防成功率高等特点[1]，是现今世界各国的研究热点。高超声速飞行器工作在大空域、宽速域等复杂多变的飞行环境中，相对于亚声速及超声速飞行器其飞行特性有着明显差异性，具有更强的不确定性，具体体现在：(1)高超声速飞行器的飞行环境复杂多变，获得精确的气动特性非常困难，通过风洞试验来获取准确气动参数的方法并不适用于高超声速飞行器[2];(2)临近空间的大气运动复杂，各种难以预测的扰动对高超声速飞行器有明显影响[3];(3)高超声速飞行器一般机体较轻，机体易产生弹性形变，且超高速飞行下的气动加热会降低机体的刚性，造成高超声速飞行器动力学模型的强不确定性[4];(4)传感器、执行器等飞行控制系统硬件自身误差[5]。

由外部复杂的飞行环境以及内部动力学飞行特性带来的不确定性因素对高超声速飞行器的飞行状态有巨大影响，而现今全球范围内高超声速飞行器的技术跨度大、成熟度较低且地面试验和飞行试验匮乏，所以通过理论研究来提高内外部双重不确定性因素影响下的系统状态估计精度，对高超声速飞行器的故障诊断、轨迹预测[6]、目标跟踪、抗扰动等领域，以及提升高超声速飞行器控制系统的安全性和可靠性都具有十分重要的研究意义。高超声速飞行器的内外部双重不确定性因素可采用系统状态演化方程和量测方程中的未知干扰项进行建模[7]，高超声速飞行器双重不确定性状态估计问题则转化为了含双重未知干扰项的系统状态估计问题。目前解决含未知干扰项的状态估计问题的算法大致可分为以下几类：

1)状态演化方程中含有未知干扰项：文献[8]针对含有状态未知干扰的系统提出了最小上界滤波器，该滤波器在卡尔曼滤波器的基础上引入自适应调整因子刻画状态未知干扰。但是该方法得到的估计结果并非传统最小均方误差意义下最优。文献[9]针对含未知过程输入的离散系统提出了基于方程式的未知输入滤波器，该方法通过矩阵函数估计未知干扰，并且利用未知干扰的估计值计算最优滤波增益。但该方法假设未知干扰和滤波误差满足一定的约束条件。文献[10]针对含未知干扰输入的多传感器离散系统，基于典范型分解，给出了不依赖于未知干扰的线性无偏最小方差最优状态滤波器。但该方法假设系统变换后的量测矩阵和未知干扰分布矩阵满足解耦约束条件。文献[11]针对含有状态未知干扰和执行器故障的线性离散不确定系统，给出了一种基于受限系统模型的状态估计算法。但是该算法认为未知干扰是范数有界的。

2)量测方程中含有未知干扰项：文献[12]针对含量测未知干扰的马尔可夫跳跃线性系统，将系统建模为广泛意义下具有随机权重的多面体结构，给出了该结构相应的上界滤波器。文献[13]针对量测方程含有结构化未知量测干扰的系统提出一种迭代最小上界滤波器，但该方法将结构化未知干扰的方差建模为对角阵。上述方法均使用了基于鲁棒估计准则的上界滤波器解决估计问题，但估计结果并非传统最小均方误差意义下最优。文献[14]针对含量测未知干扰系统提出了一种滤波器，基于线性无偏最小方差估计准则，通过构造含拉格朗日因子的性能指标并对其求导给出了滤波器参数设计。但是该方法要求未知干扰分布矩阵的维数不能超过量测的维数。

3)状态演化方程和量测方程中均含有未知干扰项：文献[15]针对幅值未知但波形已知的干扰作用于系统的情况，通过构造一个系统输出并为其设计观测器进而同时估计状态变量和未知干扰。文献[16]考虑状态方程和量测方程同时受到未知干扰的系统，基于线性无偏最小方差估计准则设计了未知输入T-S观测器。但上述两种方法均假设状态干扰和量测干扰为相同的干扰。文献[17]针对状态方程和量测方程同时受到未知干扰的系统，采用自校正外部模型，给出了基于加权最小二乘递推算法和卡尔曼滤波器两种新的状态估计算法。但是二者均假设未知干扰分布函数已知。文献[18]针对未知干扰与状态不独立情况下的含双重未知干扰的离散随机系统，提出了两阶期望最大化算法来辨识未知干扰。但该算法假设未知干扰建模为均值和方差未知的高斯项。文献[19]针对含有未知干扰的线性时变离散随机系统，提出了状态变量和故障同时估计的扩维鲁棒三阶卡尔曼滤波器。但该方法将未知干扰建模为宽平稳过程的随机游走噪声。

总结来看，现有的针对状态演化方程和量测方程中均含有未知干扰项的状态估计算法存在以下局限性：(1)要求状态干扰和量测干扰为相同的干扰。(2)要求未知干扰分布具有一定的先验信息。对于高超声速飞行器控制系统，复杂飞行环境下的未知扰动和动力学飞行特性带来的不确定性从产生机理以及数学模型上都是完全不相同的，也没有任何先验信息可取，针对此类问题现有的研究方法束手无策。因此，本文将受到内外部双重不确定性影响的高超声速飞行器建模为状态演化方程和量测方程含有不同未知干扰输入项的随机系统，设计了基于自适应方差极小化的递推状态估计器。首先建立状态估计递推模型，然后解耦滤波误差中的量测未知干扰，之后引入自适应调整因子刻画状态未知干扰，最后利用最小方差估计准则设计了滤波器中的待设计参数矩阵。以外部突风和内部传感器故障为例，受内外部双重不确定性因素影响下的高超声速飞行器仿真实验验证了本文算法的有效性。

1 问题提出

考虑线性离散系统

xk+1=Akxk+Bkuk+qk+dk

(1)

yk+1=Ck+1xk+1+ηk+1+Dk+1bk+1

(2)

式中：xk+1∈Rn为系统的状态向量，yk+1∈Rm为量测向量，uk为已知的控制输入，qk∈Rn和ηk+1∈Rm分别为满足假设1的系统噪声和量测噪声。dk∈Rn和bk+1∈Rp分别为未知状态干扰输入矢量和未知量测干扰输入矢量。Ak、Bk、Ck+1和Dk+1为适当维数的矩阵。

假设1.过程噪声qk∈Rn和量测噪声ηk+1∈Rm分别为零均值的高斯白噪声并且满足

(3)

式中：Qk≥0为系统噪声方差阵，Rk+1≥0为量测噪声方差阵，δkj为Kronecker符号。

假设2.rank(Dk+1)=p,m>p，rank(*)代表矩阵的秩。

注1.假设1采取的高斯白噪声是滤波估计、信号处理、通信系统、深度学习等诸多领域理论研究中最常采用的噪声模型，具有公式简洁、便于分析计算、算法设计等优势。依据中心极限定理，高斯白噪声可以很好地模拟真实系统中的加性噪声。过程噪声和量测噪声不相关的假设可以避免算法推导过程中交叉项的引入，降低公式复杂度。假设2是实现量测未知干扰解耦的基本条件，要求量测向量维数m必须大于未知量测干扰维数p，否则便没有足够的信息来解耦未知量测干扰。

状态未知干扰dk项除了可以描述加性干扰外，还可以描述各种系统建模的不确定性，例如系统建模过程中的非线性项、模型降阶简化以及关键参数扰动等引入的状态未知干扰输入。量测未知干扰Dk+1bk+1可以描述传感器标定误差、外界环境对传感器的干扰以及信息传输过程中的信息丢失等引入的量测未知干扰。当dk=0，Dk+1bk+1不为0时，式(1)和(2)表示的系统等同于只受到量测干扰的系统[20]。当Dk+1bk+1=0，dk不为0时，式(1)和式(2)表示的系统等同于只受到状态干扰的系统[21]。此外，若将状态未知干扰建模为dk=Ekεk，则系统模型等同于状态干扰分布矩阵已知的模型[19]，若要求dk=bk+1，则系统模型等同于状态干扰与量测干扰相同的模型[15-16]。由此可知，式(1)和式(2)可以表征更为普遍意义下含双重未知干扰的随机不确定离散系统数学模型。

本文的目的是针对状态方程和量测方程含不同未知干扰(且无任何先验信息)的高超声速飞行器系统，基于自适应估计的思想，设计一种基于自适应方差极小化的递推状态估计器。

2 基于自适应方差极小化的递推状态估计器设计

针对式(1)和式(2)建模系统的状态估计问题，最直观的解决办法是对状态未知干扰和量测未知干扰进行双重解耦，获取自适应方差极小化意义下的最优估计。但由于本文考虑的系统中状态未知干扰的干扰分布矩阵未知，且这是采用干扰解耦方法的根本条件，因此本文建模的系统无法采用文献[23]中干扰解耦的方法实现状态干扰解耦。另外，采用干扰解耦的办法对量测干扰实施解耦会牺牲待设计矩阵参数的部分自由度，而剩余的自由度能否对状态干扰实施二次解耦也是值得商榷的问题，需要针对系统模型具体问题进行具体分析。此外，针对本文所考虑的系统，扩维卡尔曼滤波[19]也是一种较为直接的方法，该方法通过将未知状态干扰以及量测干扰引入扩维状态变量，较好地解决了状态干扰和量测干扰之间的耦合问题。但由于双重干扰的强不确定性，扩维卡尔曼滤波的计算量和滤波误差会随着扩维状态变量的维数上升而急剧增加。文献[8]中鲁棒估计方法很好地规避了上述所提的参数设计自由度不足以及滤波误差发散的问题，鲁棒估计方法通过引入自适应调整因子来刻画未知干扰，无需牺牲待设计滤波增益矩阵的自由度。本文针对状态方程和量测方程含双重未知干扰的随机不确定系统，拟采用干扰解耦的办法实现量测未知干扰解耦，而针对状态干扰则采用引入自适应调整因子求取自适应状态估计的方法。

本文针对含双重干扰系统的状态估计问题，提出了一种基于自适应方差极小化的递推状态估计器。首先建立状态估计递推滤波器，实现滤波误差中的量测未知干扰解耦，之后引入自适应调整因子刻画状态未知干扰并得到了最小上界滤波误差协方差矩阵，最后利用最小方差估计准则设计了滤波器中的量测增益反馈矩阵。

2.1 量测干扰解耦递推滤波器设计

针对系统(1)～(2)，设计具有如下形式的状态估计递推滤波器

(4)

ek)+TkBkuk+Lk+1(Ck+1xk+1+ηk+1+

Dk+1bk+1)=Fk(xk-ek)+TkBkuk+Lk+1·

Ck+1(Akxk+Bkuk+qk+dk)+Lk+1ηk+1+

Lk+1Dk+1bk+1=(Fk+Lk+1Ck+1Ak)xk-

Fkek+(Tk+Lk+1Ck+1)Bkuk+Lk+1Ck+1·

(qk+dk)+Lk+1ηk+1+Lk+1Dk+1bk+1

(5)

联立式(1)和式(5)得到k+1时刻滤波误差为

(I-Lk+1Ck+1-Tk)Bkuk+(I-Lk+1Ck+1)·

(dk+qk)-Lk+1ηk+1-Lk+1Dk+1bk+1

(6)

对任意的未知输入bk+1和dk,为了实现量测干扰和滤波误差解耦，引入如下的解耦约束条件。

(7)

将式(7)代入式(6)，状态估计误差化简为

ek+1=Fkek-Lk+1ηk+1+(I-Lk+1Ck+1)(dk+qk)=

(I-Lk+1Ck+1)Akek-Lk+1ηk+1+Tk(dk+

qk)=Tk(Akek+dk)+Tkqk-Lk+1ηk+1

(8)

根据滤波误差方程(8)，可得滤波误差方差阵为

(9)

式(9)中的状态未知干扰dk阻碍了滤波协方差的计算，需要引入自适应调整因子αk+1(αk+1≥1)刻画状态未知干扰并定义滤波误差协方差的上界为

(10)

2.2 自适应调整因子设计

自适应调整因子可以通过滤波残差协方差矩阵求取。首先计算滤波残差：

(11)

对上式左乘Lk并计算滤波残差协方差Vk+1，利用约束条件(7)将量测干扰从滤波残差中解耦可得：

(12)

(13)

其中自适应调整因子的求解方法由定理1给出。

定理1.如果以下三个条件均满足：

2)Ck+1是满秩的，即rank{Ck+1}=n

(14)

(15)

(16)

集合Ωk+1为

(17)

证.

1) 最优自适应调整因子的存在性：

定义中间变量ΔΨk为

(18)

由式(12)～(13)得：

(19)

由Lk+1Ck+1满秩可得：

ΔΨk≥0

(20)

由Tk+1满秩以及式(9)～(10)可得：

(21)

故可得到下列不等式：

(22)

(23)

2) 最优自适应调整因子的求解：

(24)

(25)

(26)

由式(25)～(26)有

(27)

(28)

结合式(22)～(23)有

即式(14)～(15)成立。

(29)

(30)

(31)

2.3 量测增益反馈矩阵设计

(32)

式中:Λk+1，H和G为中间变量。

(33)

(34)

(35)

证.

应用最小方差估计准则，并由约束条件(7)可引出如下辅助方程：

(36)

将式(31)代入式(36)得：

tr{Qk+Lk+1Ck+1Qk(Lk+1Ck+1)T}-

(37)

(38)

定义中间变量H和G如式(33)～(34)，则式(38)可以表示为

(39)

将式(39)和约束条件(7)联立得到矩阵方程组：

(40)

写为分块矩阵形式：

(41)

由假设2可知分块矩阵方程的系数矩阵的逆存在，解分块矩阵方程(41)得：

即式(32)和式(35)成立。

2.4 算法实现步骤

基于2.1、2.2和2.3节的推导和证明，针对高超声速飞行器双重不确定性系统，本文所提的基于自适应方差极小化的递推状态估计器的计算步骤如下：

1)k=1，设置状态初值x0，误差协方差初值P0；

2)由式(11)计算残差γk+1；

6)将参数代入式(33)和式(34)计算中间变量H,G；

7)将中间变量H和G代入式(35)计算Λk+1；

11)k+1，返回步骤2)。

3 仿真校验及结果分析

3.1 仿真环境

采用文献[24]给出的高超声速飞行器纵向小扰动线性模型进行仿真验证。其中：状态变量为x=[v,γ,α,q,h]T(v,γ,α,q,h分别为飞行速度，航迹倾角，迎角，俯仰角速度和高度)，控制输入为uk=[δe,η]T(δe,η分别为升降舵偏角，油门设定)，状态方程中Ak，Bk和过程噪声协方差矩阵Qk分别为

未知状态扰动dk设定为外部突风[21]，设置如下：前10 s未知干扰为0，第10 s到第20 s受到幅值为1的阶跃响应，第20 s到第30 s受到幅值为-1的阶跃响应，从第30 s一直到第100 s未知干扰始终为0。

控制输入为uk=[-0.2389, 0.1566]T，初始状态设置x0=[4525.6, 0, 0.978, 0, 30000]T，状态协方差P0=0.12I5×5。

3.2 系统状态估计结果及分析

本节对所提出的基于自适应方差极小化的递推状态估计器(AVMRE)的有效性进行验证，并与KF和最小上界滤波器[8](Minimumupper-bound filter, MUBF)进行对比。图1和图5分别为不同算法对第一维和第五维状态变量的估计效果对比，这两维度都不受突风和传感器故障的影响，可以看出本文所提算法、MUBF与KF算法估计效果相差不大，都可较好地跟踪真实值。

图2为第二维状态变量估计效果对比，由状态干扰dk和量测干扰Dk+1bk+1可知该变量只受到传感器故障影响，当量测变量在50～70 s受到传感器故障影响时，本文所提算法对故障干扰进行了解耦设计，故估计效果较好，而未考虑故障干扰解耦的KF和MUBF估计效果不佳。

图3为第三维状态变量估计效果对比，由状态干扰dk和量测干扰Dk+1bk+1可知该维度变量只受到突风干扰影响，可以看出本文所提算法和MUBF由于引入自适应调整因子刻画突风干扰，进而在0～50 s的估计值较好地跟踪上了真实值，而KF由于忽略了突风的影响导致估计效果不佳。同时可以看出由于传感器故障的影响，MUBF和KF的量测值已经产生了偏差，导致相应维度的状态估计效果也偏离了真实值。

图4为第四维状态变量估计效果对比，可以看出本文算法和MUBF算法受到突风干扰估计效果相差不大，而KF算法估计效果不佳。系统在50～70 s受到故障干扰， KF和MUBF由于忽略了故障干扰的影响导致估计效果较差，本文算法通过对故障干扰进行了解耦设计，故估计效果依然较好。

为了更清晰地对比三种算法的估计性能，下面给出了100次蒙特卡洛仿真下状态估计RMSE对比。为了节省篇幅，只给出了第三维和第四维的结果图，可以看出本文所提算法的RMSE较小，而KF和MUBF在50～70 s由于传感器故障的存在，两种算法的RMSE发生了突变。总结来说，从图6可以看出，本文所提算法对于单独受到突风干扰或传感器故障影响的高超声速飞行器的状态估计效果较好；从图7可以看出，对同时受到突风和传感器故障双重干扰影响的高超声速飞行器，本文所提算法依

图1 第一维状态估计对比

图2 第二维状态估计对比

图3 第三维状态估计对比

图4 第四维状态估计对比

图5 第五维状态估计对比

图6 第三维的状态估计RMSE对比

图7 第四维的状态估计RMSE对比

然可以获得令人满意的状态估计结果。

4 结 论

本文将高超声速飞行器中的不确定性因素建模为未知干扰输入项，针对状态方程和量测方程含有不同未知干扰输入的高超声速飞行器系统状态估计问题开展了研究，提出了一种基于自适应方差极小化的递推状态估计器。首先建立了状态估计递推滤波器，实现状态估计误差中的量测未知干扰解耦，之后引入自适应调整因子刻画状态未知干扰并得到了最小上界估计误差协方差矩阵，最后，利用最小方差估计准则设计了滤波器中的量测增益反馈矩阵，实现了自适应方差极小化意义下的状态估计。以外部突风和内部传感器故障为例，受内外部双重不确定性因素影响下的高超声速飞行器仿真实验验证了本文算法的有效性。