飞行器上的IMU包括三轴加速度计和三轴陀螺仪。三轴加速度计测量除重力之外的合力引起的飞行器沿机体系三轴方向上的加速度,除重力之外的力有气动力和发动机的推力。三轴陀螺仪测量飞行器绕机体系三轴的角速率。IMU是飞行器的核心传感器,它的故障如果没有被及时诊断将会导致严重的飞行事故。1995年,美国Perseus飞行器由于俯仰角速度传感器故障导致飞行速度超出正常范围限制,机翼承受过大气动力而发生折断[1]。2005年,马航MH124客机从珀斯出发前往吉隆坡,在珀斯西北240 km印度洋上空发生异常情况。遭遇主飞行仪表数据显示矛盾、自动右转向及自动超高度飞行等一系列问题,事故原因为大气数据惯性基准组件中的加速度传感器故障[1]。
现阶段,IMU的故障诊断主要通过硬件余度和表决检测来实现。但是硬件余度存在成本高、质量重和体积大等缺点,且相似余度的传感器配置易发生共性故障,对特定类型的故障难以进行有效的诊断。解析余度是基于软件的一种故障诊断方法,它有效地规避了硬件余度的上述缺点。但是解析余度的故障诊断算法依赖于准确的系统建模,模型不确定性是阻碍解析余度应用的一个主要问题。故障诊断算法的设计需要尽可能降低对模型不确定性的敏感性甚至不敏感,但同时要求不损失对故障的敏感性,这是解析余度故障诊断算法设计的难点[2]。
在故障诊断算法设计中,为了降低对模型不确定性的敏感性,通常使用卡尔曼滤波器来进行故障估计。噪声是一种典型的模型不确定性,卡尔曼滤波假设传感器噪声服从零均值的高斯分布。由于高斯分布的线性传播封闭性,通过均值和协方差2个参数足以表征噪声通过线性系统传递的特性。对于线性系统,在噪声高斯分布假设下,卡尔曼滤波可以得到最优的状态估计,基于卡尔曼滤波器的故障诊断算法可以大幅降低对传感器噪声的敏感性。针对非线性系统而言,扩展卡尔曼滤波[3]通过在上一时刻估计处进行泰勒级数展开,保留线性部分,进而将非线性系统局部线性化,将卡尔曼滤波推广到了非线性系统。但是扩展卡尔曼滤波存在截断误差积累的问题,另一方面,雅克比矩阵计算较麻烦,在强非线性情况下,数值稳定性较低。无迹卡尔曼滤波[4]同样假设噪声服从高斯分布,用一组确定的Sigma点表征高斯分布的均值和协方差,Sigma点通过非线性系统传播了高斯分布的均值和协方差。经过证明,无迹卡尔曼滤波至少可以达到二阶扩展卡尔曼滤波的精度[5],是适用于非线性高斯系统的一种较为优秀的滤波方法。
卡尔曼滤波框架下的故障诊断算法,为了不损失对故障的敏感性,通常将故障建模成随机游走模型[6-8]。参考文献[9]中,将飞行器的IMU故障建模成随机游走模型,并且利用OTSUKF技术解耦IMU故障和系统状态,但是随机游走模型的协方差矩阵是按照经验给定的,当该协方差矩阵偏离经验值时,故障诊断的效果将大幅度降低甚至导致结果不可用。针对于滤波器的数学模型与实际情况不符合的问题,Mehra提出了新息协方差匹配自适应滤波技术[10],保证了卡尔曼滤波器的鲁棒性。随机游走模型中的协方差反应了故障的动态特性,本文将新息协方差匹配技术与OTSUKF技术相结合,提出了ATSUKF。该方法在线自适应调整随机游走模型的协方差矩阵,不需要按照经验给定常值的协方差矩阵,提升了OTSUKF的故障诊断性能和鲁棒性。
1 问题描述 飞行器主动容错控制[11]的结构如图 1所示,本文研究的问题是设计一种故障诊断算法,可以融合全球定位卫星系统(GNSS)和故障的IMU信号,实现了IMU的故障诊断。
2 系统建模 2.1 非线性系统模型 基于解析余度的故障诊断算法设计依赖于飞行器的非线性系统模型,文献[12]推导了适用于飞行器IMU故障诊断的非线性系统模型。系统方程为:
|
(1) |
式中:u, v, w表示飞行器速度在机体坐标系三轴的投影; φ, θ, ψ表示姿态角; axm, aym, azm表示加速度计的测量值; pm, qm, rm表示飞行器绕机体坐标系三轴旋转的角速率测量值; fax, fay, faz表示加速度计发生的故障幅值; fp, fq, fr表示陀螺仪发生的故障幅值; ωax, ωay, ωaz表示加速度计的测量噪声; ωp, ωq, ωr表示陀螺仪的测量噪声。观测方程为
|
(2) |
式中:uGSm, vGSm, wGSm表示GNSS的速度测量值; 多天线GNSS可以测量飞行器的姿态角[13], φm, θm, ψm表示GNSS的姿态角测量值; νu, νv, νw表示GNSS的速度测量噪声; νφ, νθ, νψ表示GNSS的姿态角测量噪声; RDCM表示方向余弦矩阵, 它反映了地球坐标系到机体坐标系的坐标转换关系, 它的定义为
|
(3) |
2.2 连续时间状态滤波模型 定义系统模型的状态向量、输入向量、IMU故障向量、观测向量、IMU噪声向量和观测噪声向量分别为
|
(4) |
|
(5) |
|
(6) |
|
(7) |
|
(8) |
|
(9) |
非线性系统模型可以写成如下形式
|
(10) |
系统方程中的噪声为非线性噪声, 对于包含非线性系统噪声的滤波模型, 无迹卡尔曼滤波有2种处理方法, 一种是将非线性噪声增广为状态[14], 另外一种方式将非线性噪声通过线性化方法转换为加性噪声。线性化方法是
|
(11) |
|
(12) |
式中:矩阵G(t)为噪声分布矩阵;矩阵F(t)为故障分布矩阵。通过对噪声和故障的线性化处理, 得到连续时间状态滤波模型
|
(13) |
2.3 离散时间状态滤波模型 无迹卡尔曼滤波假设连续时间的滤波模型(13)的噪声服从零均值的高斯分布
|
(14) |
|
(15) |
|
(16) |
式中:q表示系统噪声强度矩阵; R表示测量噪声强度矩阵; δij表示狄利克雷函数。对连续的滤波模型(13)的系统方程按欧拉方法进行离散化
|
(17) |
离散的滤波模型观测方程可以写成
|
(18) |
定义离散系统噪声和测量噪声为
|
(19) |
|
(20) |
离散后的系统噪声强度矩阵
|
(21) |
|
(22) |
则离散滤波模型(17)式和(18)式的噪声强度矩阵
|
(23) |
|
(24) |
高斯分布经过线性变换仍然为高斯分布, 所以离散的滤波模型系统噪声仍然服从零均值的高斯分布
|
(25) |
|
(26) |
定义
|
(27) |
|
(28) |
可以得到离散时间状态滤波模型
|
(29) |
3 自适应无迹卡尔曼滤波 3.1 最优二步无迹卡尔曼滤波及其缺点 无迹卡尔曼滤波的关键技术是UT变换, UT变换依赖一组确定性的Sigma点, 这组Sigma点反映了高斯分布噪声的均值和协方差。考虑如下包含加性噪声的离散时间状态滤波模型
|
(30) |
文献[4-5]推导并给出了无迹卡尔曼滤波算法。
离散时间状态滤波模型(30)式没有考虑到传感器故障的影响, 无迹卡尔曼滤波不能直接应用于IMU的故障诊断。针对包含传感器故障的离散时间状态滤波模型(29)式, 文献[9]提出了适用于飞行器IMU故障诊断的OTSUKF。在文献[15]中, 给出了二步卡尔曼滤波的最优解以及证明。
OTSUKF并列运行一个无偏差滤波器和一个偏差滤波器, 通过耦合矩阵融合2个滤波器的估计值得到最终的状态估计和故障估计。其中无偏差滤波器是基于忽略故障信息的滤波模型设计的, 偏差滤波器是基于故障的滤波模型设计的。偏差滤波器的输入是无偏差滤波器的新息, 而故障的滤波模型建模成随机游走模型
|
(31) |
|
(32) |
|
(33) |
式中,Qkf是随机游走模型的协方差矩阵。
随机游走模型协方差矩阵对角线元素的值反映故障变化的快慢, 对角线元素的值越大, 表征故障变化越快。该协方差矩阵非对角线元素的值反映了不同传感器故障变化的线性相关程度。应用OTSUKF进行飞行的IMU故障诊断存在的问题是:随机游走模型的协方差矩阵的值是按照设计者的经验给定的一个常值矩阵, 随机游走模型的协方差矩阵很难与实际故障的动态特性相匹配, 当协方差矩阵与实际的故障特性匹配较差的时候, OTSUKF的故障诊断性能会大幅下降。
3.2 新息协方差匹配 ATSUKF的自适应方法设计依据于新息协方差匹配技术。传感器的故障建模成随机游走模型, 该模型是线性模型。针对线性系统
|
(34) |
文献[3, 5]推导并给出了卡尔曼滤波算法。卡尔曼滤波中, 新息协方差按照公式(35)计算
|
(35) |
真实的新息协方差可以根据最近的m步的新息进行辨识
|
(36) |
根据新息协方差匹配, 得到以下关系
|
(37) |
则噪声ωk的协方差矩阵Qk可以通过下式估计
|
(38) |
式中:Ak是系统转移矩阵, Hk+1是观测矩阵。
针对于包含传感器故障的滤波模型(29)式而设计的自适应故障滤波器(41)中, 故障的系统转移矩阵为单位矩阵, 故障的观测矩阵为Sk+1且该矩阵为非奇异矩阵, 得到随机游走模型的协方差矩阵的估计公式
|
(39) |
3.3 自适应二步无迹卡尔曼滤波 ATSUKF在OTSUKF的基础上, 将故障滤波器改进为自适应故障滤波器。ATSUKF通过新息协方差匹配技术, 自适应地获取随机游走模型的协方差矩阵, 解决了OTSUKF中存在的随机游走模型的协方差矩阵与故障动态特性不匹配的问题。
伪代码(40)至(42)式给出了ATSUKF算法。无故障滤波器如(40)式所示, 自适应故障滤波器如(41)式所示, 融合算法如(42)式所示。其中sigma(·)表示Sigma点的生成, f-UT(·)和h-UT(·)分别表示Sigma点经过滤波模型的系统方程和观测方程的UT变换。无故障滤波器的输入参数xkff和Pkff是上一次迭代的输出; QkRk+1根据公式(23)、(24)计算; uk是系统的输入; yk+1是系统的输出; ukopt和Qkopt是上一次迭代的自适应故障滤波器的输出。故障的估计
是通过自适应故障滤波器得到的, 它的输入参数fkiPkf和
是上一次迭代的输出; Uk+1是融合算法的输出; yk+1f是无故障滤波器的输出; m是估计随机游走模型协方差矩阵所使用新息的长度。系统的状态估计xk+1est和Pk+1est是通过融合算法得到的。融合算法的输入参数Vk是上一次迭代的输出; Fk是故障分布矩阵, 按照公式(27)计算; Qkxf是IMU噪声与随机游走模型噪声的协方差, IMU噪声和随机游走模型噪声是相互独立的, 它的值取零[15];
是随机游走模型的协方差矩阵估计, 它是自适应故障滤波器的上一次迭代的输出;
和Kk+1ff是无故障滤波器的输出;
和Sk+1是自适应故障滤波器的输出。
|
(40) |
|
(41) |
|
(42) |
4 仿真实验与结果分析 4.1 仿真实验 仿真实验分别用二组不同的值初始化随机游走模型的协方差矩阵, 对比验证了OTSUKF方法和ATSUKF的IMU故障诊断效果。
MATLAB仿真中, 采样时间取0.01 s, 分别在10~20 s, 25~35 s, 40~50 s给加速度计和陀螺仪注入阶梯故障F1、斜坡故障F2和正弦故障F3。故障注入类型和时间如表 1所示, 其中加速度计故障的单位是m/s, 陀螺仪故障的单位是rad/s。
表 1 加速度计和陀螺仪典型故障注入
f |
t |
10~20 s |
25~35 s |
40~50 s |
fax |
F1 |
-F2 |
F3 |
fay |
F2 |
-F3 |
F1 |
faz |
F3 |
-F1 |
F2 |
fp |
F1 |
-F2 |
F3 |
fq |
F2 |
-F3 |
F1 |
fr |
F3 |
-F1 |
F2 |
注入故障的数学表达式为
|
(43) |
当随机游走模型的协方差矩阵Qkf初始化为
|
(44) |
基于OTSUKF和ATSUKF的IMU故障诊断结果分别如图 2至3所示。
当随机游走模型的协方差矩阵Qkf初始化为
|
(45) |
基于OTSUKF和ATSUKF的IMU故障诊断结果分别如图 4和图 5所示。
4.2 结果分析 本文引入均方误差作为评价故障诊断算法准确度的指标, 均方误差的定义为
|
(46) |
IMU故障诊断算法的均方误差如表 2所示, 表 2前3列是OTSUKF和ATSUKF分别在随机游走模型的协方差矩阵初始化为公式(44)情况下的结果, 后3列是OTSUKF和ATSUKF分别在随机游走模型的协方差矩阵初始化为公式(45)情况下的结果。
表 2 故障诊断的均方误差
RMSE1 |
OTSUKF |
ATSUKF |
RMSE2 |
OTSUKF |
ATSUKF |
fax |
0.249 0 |
0.092 5 |
fax |
0.333 9 |
0.092 5 |
fay |
0.241 5 |
0.086 2 |
fay |
0.332 0 |
0.086 2 |
faz |
0.241 1 |
0.087 7 |
faz |
0.331 9 |
0.087 7 |
fp |
0.247 9 |
0.088 3 |
fp |
0.331 7 |
0.088 3 |
fq |
0.246 6 |
0.087 1 |
fq |
0.332 1 |
0.087 1 |
fr |
0.243 8 |
0.087 6 |
fr |
0.331 2 |
0.087 6 |
图 2和图 4的结果反映了基于OTSUKF方法的IMU故障诊断的性能对于随机游走模型协方差初始值的鲁棒性较差。在工程应用中, 为了得到满意的故障诊断性能, 参数整定的难度较大。图 3和图 5反映了ATSUKF在随机游走模型的协方差取不同的初始值情况下, IMU故障诊断均具有较好的性能。表 2的第3列和第6列反映了二次实验的ATSUKF故障诊断算法的均方误差。二次实验的均方误差一致表明了ATSUKF方法对于随机游走模型协方差初始值的鲁棒性较好。对比表 2的第2列和第3列、表 2的第3列和第6列表明基于ATSUKF的IMU故障算法较OTSUKF算法具有更好的性能。
OTSUKF方法中的随机游走模型协方差是一个常值矩阵, 值为设计者给定的初始值。而ATSUKF方法中随机游走模型协方差是一个变化的矩阵, 其取值是根据协方差匹配得到。故障建模成随机游走模型(31), 其噪声项协方差反映了故障变化特性。协方差矩阵对角线元素的值反映故障变化的快慢, 对角线元素值越大, 表征故障变化越快。协方差矩阵的非对角线元素的值反映不同传感器故障变化的线性相关程度。在IMU故障诊断中, 各通道的故障是不相关的, 所以协方差矩阵的非对角元素为零。
传感器故障的动态特性是不可预知的, 可以将故障建模成带未知协方差矩阵的随机游走模型, 随机游走模型的协方差矩阵反映了故障的动态特性。图 6给出了ATSUKF在初始化参数1情况下对角线元素值变化的情况。当故障发生的时候, ATSUKF会根据故障的变化快慢自适应的调整对角线元素的值, 而OTSUKF的协方差矩阵是一个常值矩阵。相较于OTSUKF, ATSUKF方法的滤波模型可以更好地匹配实际系统, 所以具有更好的故障估计效果。另一方面, 在ATSUKF方法中, 随机游走模型的协方差矩阵是基于协方差匹配技术在线计算的, 它对于协方差矩阵的初始值具有较强的鲁棒性。
4.3 算法实时性分析 IMU的更新频率为1 000 Hz, GNSS测量的速度和姿态更新频率为100 Hz。仿真实验中, OTSUKF和ATSUKF算法的更新频率设置为100 Hz。为了避免飞行器仿真器对故障诊断算法运行速度的影响, 提取飞行器正常飞行60 s的传感器数据。在IMU的数据中, 按照表 1的方式注入故障。基于故障数据, 进行故障诊断的蒙特卡洛仿真实验。仿真环境和算法消耗时间如表 3所示。
表 3 仿真实验环境与算法消耗时间
硬件 |
型号 |
软件 |
名称 |
算法 |
时间/s |
CPU |
i7-7700HQ |
仿真 |
MATLAB |
ATSUKF |
1.36 |
内存 |
16GB |
系统 |
macOS |
OTSUKF |
1.17 |
实际飞行60 s的数据, 基于ATSUKF的故障诊断算法需要处理1.36 s, 基于OTSUKF的故障诊断算法需要处理1.17 s。自适应的随机游走模型协方差矩阵的辨识涉及矩阵求逆操作, 是造成ATSUKF运行较慢的原因。根据实验结果, 所提出的IMU故障诊断算法满足实时性的要求。
5 结论 本文提出了一种自适应二步无迹卡尔曼滤波方法,融合全球定位卫星系统和故障的惯性测量单元传感信号,实现了惯性测量单元的故障诊断。一方面,自适应二步无迹卡尔曼滤波对传感器噪声进行优化,得到状态和故障的估计,降低了故障诊断算法对于噪声的敏感性。另一方面,得益于随机游走模型和新息协方差匹配自适应技术的引入,提升了故障诊断算法对于故障的敏感性。仿真实验结果表明,针对于飞行器惯性测量单元的故障诊断,所提出的方法具有较好的性能和鲁棒性。
Aircraft Inertial Measurement Unit Fault Diagnosis Based on Adaptive Two-Stage UKF