现在的大型民用航空飞机一般采用多套惯性导航系统和多套大气数据系统，利用冗余信息对导航系统的完好性进行监测，提高导航系统的可靠性和安全性。现在民用导航系统常见故障可分为硬故障和软故障^[1]。目前针对多余度系统的故障检测方法有等价空间法，包括广义似然比法^[2](generalized likelihood ratio test, GLT)、最优奇偶向量法^[3](optimal parity vector test, OPT)和奇异值分解法^[4]。GLT是基于残差量的最大似然估计故障检验方法，该方法计算量小、实时性强，常用于余度系统的故障检测与隔离，但该方法对软故障的检测性能较低^[5-7]。这是由于GLT方法的故障检测函数只与当前时刻的噪声和可能的故障有关^[8-9]，且当有故障发生时，其故障检测函数近似与奇偶残差中的故障量平方成正比，故当发生软故障时，由于在故障刚发生的一段时间内故障的幅值较小，导致故障检测函数值持续低于故障检测门限，这就导致该方法对软故障的检测不灵敏，存在较大的延时。

在故障检测中，残差信息经常被用来反映系统故障，序贯概率比检验(sequential probability ratio test, SPRT)就是一种用残差信息进行故障判决的常用方法。文献[10]针对服从正态分布的残差，提出了一种改进序贯概率比(improved sequential probability ratio test, IM-SPRT)方法，该方法采用迭代递推的方法计算故障检测值，充分利用了之前时刻的故障测量值，可有效提高对软故障的检测灵敏度，但文献[11-12]中指出，IM-SPRT方法无法判断故障结束时间，会造成很大的虚警。

本文结合GLT和IM-SPRT方法各自的优势，提出了一种基于GLT/IM-SPRT的多余度系统高灵敏度故障检测方法。设计了多余度大气/惯性系统高灵敏度故障检测与隔离方法框架，实现了灵敏度高，漏检率低的故障检测与隔离。对比该方法与GLT方法的故障检测结果，验证了该方法的有效性。

1 GLT与IM-SPRT性能分析 1.1 GLT方法

有故障发生时, m个余度传感器的量测方程为

(1)

式中：X∈Rⁿ是n个待测状态变量；Z∈R^m是m个传感器的量测值(m≥n); H是传感器安装矩阵, f是故障向量；ξ是m维零均值、协方差为σ²的高斯白噪声序列, 即E(ξ)=0, E(ξξ^T)=σ²I_m, I_m为单位矩阵。

通过构造满足条件V′H=0与条件V′V=I_m-n(V′为V的转置)的奇偶矩阵V使得

(2)

式中：P是对状态X解耦的奇偶残差, 仅与噪声或可能的故障有关, 当传感器无故障时；P仅是噪声的函数, 即

(3)

1) 故障检测

GLT方法故障检测决策的假设检验:

无故障状态H₀:E(P)=0, E(PP′)=σ²V′V;

有故障状态H₁：

E(P)=μ≠0, E((P-μ)(P-μ)′)=σ²V′V; 其中μ=V′f。

利用奇偶残差P构造故障检测判决函数

(4)

式中，F_{D, GLT}~χ²(m-n)。

则故障检测过程可描述为

T_d是先验故障检测门限。在已知故障检测虚警率及m-n下, 通过χ²分布表获得。当故障检测门限T_d增大时, 故障漏检率会增大, 软故障检测平均延迟时间会增加; 反之, 故障虚警率增大, 软故障检测平均延迟时间会减小。

2) 故障隔离

H_i表示假设第i(i=1, 2, …m)个传感器故障, 对第i个传感器构造故障隔离判决函数

(5)

V_i^*为V的第i列。故障传感器的隔离函数值是所有传感器隔离函数值中的最大者。

由公式(2)和(4)可知, GLT方法故障检测函数只和当前时刻的噪声和可能发生的故障有关, 当有故障发生时, 其故障检测函数值近似与奇偶残差中的故障信息的平方成正比。如图 1所示, 假设在时刻200~400有硬故障发生, 故障幅值会在时刻200迅速增大到较大值, 对应GLT方法故障检测函数值在时刻200处也迅速增大而超过故障检测门限, 对硬故障的检测灵敏度高; 假设在时刻200~300有软故障发生, 在故障刚发生的一段时间内, 奇偶残差中的故障特征不明显, 导致故障检测函数值很小, 不能检测出故障, 随故障幅值的持续增大, 故障检测函数值也增大, 直到增大到一定程度, 才可检测出故障, 这就导致GLT方法对软故障的检测灵敏度低。

1.2 IM-SPRT方法

设p为奇偶残差P的任意一维分量, 它的k次序贯独立样本为{p₁, p₂, …, p_k}, 由概率论与数理统计原理近似可得p~N(p_k, σ_k²), p_k为样本均值, σ_k²为样本方差, 对无故障和有故障2种状态下的概率密度函数有似然比

(6)

两边同时取对数, 得对数似然比为

(7)

式中, p₀表示无故障时的真实值。

当系统正常时, 由切比雪夫大数定理可知, 样本均值p_k依概率收敛于真值p₀, 此时, λ(k)会趋向零, 即λ(k)→0;当系统发生故障时, 随着k的增加, 样本均值p_k会逼近故障时的真值, 该值与p₀不同, 可得λ(k)的递推计算公式, 见(8)式。

(8)

式中:若变量p的先验方差σ²已知, 可直接代替σ_k²; 若变量p的先验方差σ²未知, 可实时计算样本方差代替σ_k²。

文献[10]对该方法设定一个故障检测门限T_sprt, 该门限由虚警率P_F和漏检率P_M唯一确定, 为

(9)

IM-SPRT方法故障检测门限T_sprt的变化对故障检测性能的影响分析与GLT方法相同。

则故障检测过程可描述为

由公式(8)可知, IM-SPRT方法采用迭代递推的方法来计算故障检验函数, 充分利用之前时刻的检测值, 如图 2所示, 假设在时刻200~400有硬故障发生, 故障的幅值在时刻200达到较大值, 对应IM-SPRT方法的故障检测函数值也会迅速增大到故障检测门限之上, 对硬故障的检测灵敏度高; 假设在时刻200~300有软故障发生, 该方法的故障检测函数值也是累加增大, 这样可以提高软故障的检测灵敏度。当故障结束时, 该方法的故障检测函数仍会在故障检测门限之上, 会造成很大的虚警。

此外, 对IM-SPRT方法, 将其故障检测函数λ(k)展开, 有

(10)

进一步可得

(11)

即当k值减小时, λ(k)的值会增大, 这样通过设置重置周期T周期性将k值重置为1, 可提高λ(k)值的累加速度, 进而进一步提高IM-SPRT方法对软故障的检测灵敏度。

2 基于GLT/IM-SPRT的高灵敏度故障检测方法

由第1节的分析可知, GLT方法虚警率低, 能准确的判断故障结束时间, 但对软故障的检测灵敏度低。IM-SPRT方法能快速检测出软故障, 但无法准确判断故障结束时间, 会造成很大的虚警。本文针对这些问题, 提出结合这2种方法优势的软、硬故障检测与隔离方法, 该方法用IM-SPRT方法检测故障, 通过GLT方法判断故障结束时刻来重置IM-SPRT方法的检测值, 并用GLT方法隔离故障, 方法框架见图 3。

在对奇偶残差P的每一维都用IM-SPRT方法进行故障检测的同时也用GLT方法进行故障检测。当奇偶残差P的每一维的IM-SPRT方法的检测值都没有超过给定阈值, 则系统正常; 当奇偶残差P的IM-SPRT方法检测值中有超过给定阈值的值, 则系统故障。在检测出系统故障时, 用GLT故障隔离方法隔离故障, 在GLT方法判断出故障结束时, 将IM-SPRT方法的检测值置为零, 这样当故障消失时, IM-SPRT的故障检测值可以迅速减小到检测门限以下, 并以当前时刻为起始时刻开始重新迭代递推, 这样当非第一次发生的故障出现时, 还可以继续检测故障。设计的故障检测与隔离算法流程如下:

1) 计算奇偶矩阵V, 设置重置周期T, 将k值置为零;

2) 计算奇偶残差P;

3) 判断k值是否超过重置周期T, 若超过, 则重置为1, 若未超过则继续累加, 并对奇偶残差用2种方法进行故障检测;

4) 若有检测值超出了IM-SPRT方法的故障检测门限, 则同时用GLT方法的故障隔离策略隔离故障, 并返回到步骤2);

5) 若GLT方法判断出此次故障结束, 则将IM-SPRT的故障检测值置为零, 并返回到步骤2)进行下一个故障的检测与隔离。

3 仿真分析

为了验证本文提出算法的有效性, 首先完成惯性导航系统和大气数据系统的导航解算, 在此基础上, 给获得的大气或惯性数据加均值为零特定方差大小的高斯白噪声, 产生3套带有随机误差的传感器数据, 假设其中仅有1套数据出现故障, 在实验数据中添加一定大小的硬故障或软故障信号, 分别用本文方法和GLT方法进行1 000次蒙特卡洛仿真, 计算故障检测平均延迟时间与漏检率。

本次仿真以三余度大气数据系统子系统级气压高度数据、三余度惯性导航系统传感器级东向加速度数据的故障检测与隔离为例, 给定虚警率为0.001, 本文所给方法的故障检测阈值为3, GLT方法的故障检测阈值为13.8155。本文中给定k值重置周期为100 s。

1) 气压高度

给获得的气压高度数据添加零均值, 标准差为σ=1 m的高斯白噪声, 并取其中1 000个时刻数据。

情况一  给生成的第一套数据在第200~400时间段内添加幅值10的硬故障, 在第600~700时间段内添加斜率为0.4的软故障。

本文方法和GLT方法在一次仿真中的结果如下所示, 仿真结果中故障检测标志“1”表示有故障发生, “0”表示无故障发生。

对比图 4a)至4b)和图 5a)至5b)可知, 本文方法和GLT方法均在时刻200~400检测出了故障, 对所设置的硬故障检测均无延迟, 而对于软故障的检测, 从图 4b)可以看出, 本文给出的方法在第607个时刻时检测出故障, 存在7个周期的延时, 从图 5b)可以看出, GLT方法在第612个时刻时检测出故障, 存在12个周期的延时, 相比之下, 本文所给方法的灵敏度更高。对比图 4c)和图 5c), 因为2种方法都采用GLT方法的隔离策略, 故在故障隔离上没有明显差异。

进一步给出2种方法1 000次蒙特卡洛仿真后计算出的漏检率和平均延迟时间, 见表 1, 从表中可以看出, 对于硬故障检测2种方法无明显差异, 而对于软故障, GLT方法漏检率约为本文方法的2.5倍, 平均延迟时间比本文多出5个周期, 可见本文方法对软故障的检测灵敏度更高。

情况2  给生成的第一套数据在第200~400时间段内添加幅值11的硬故障, 在第600~700时间段内添加斜率为0.3的软故障。

由于篇幅有限, 仅给出1 000次蒙特卡洛仿真后2种方法计算出的漏检率和平均延迟时间, 结果如表 2所示，可见本文对软故障的检测灵敏度更高。

2) 东向加速度

给获得的东向加速度数据添加零均值, 标准差为σ=0.02 m/s²的高斯白噪声, 取其中1 000个时刻数据。

情况一  给生成的第一套数据在第200~400时间段内添加幅值0.2的硬故障, 在第600~700时间段内添加斜率为0.005的软故障。

进一步给出2种方法1 000次蒙特卡洛仿真后计算出的漏检率和平均延迟时间, 见表 3。

从表中可以看出，对于硬故障检测2种方法无明显差异, 而对于软故障，GLT方法漏检率约为本文方法的1.6倍，平均延迟时间比本文多出9个周期，可见本文方法对软故障的检测灵敏度更高。

情况2  给生成的第一套数据在第200~400时间段内添加幅值0.22的硬故障, 在第600~700时间段内添加斜率为0.004的软故障。

同样仅给出1 000次蒙特卡洛仿真之后2种方法计算的故障检测漏检率和平均延迟时间, 结果见表 4。可知, 在此仿真条件下, 对于硬故障的检测量中方法无明显差异, 而对于软故障, GLT方法漏检率约为本文方法的1.7倍, 平均延迟时间比本文多出11个周期, 可见本文对软故障的检测灵敏度更高。

由对气压高度数据和东向加速度数据的故障检测与隔离的仿真分析结果可知, 在假设同一时刻只有一个故障发生的情况下, 本文所给的方法能有效检测软、硬故障情况, 且相比于GLT方法, 本文所给方法对软故障的检测灵敏度更高。

4 结论

为了提高GLT方法对软故障的检测灵敏度, 本文提出了一种基于GLT/IM-SPRT的多余度大气/惯性高灵敏度故障检测方法, 结合GLT方法与IM-SPRT各自的优势, 将2种方法结合使用, 用IM-SPRT方法检测故障, 用GLT方法隔离故障, 并通过GLT方法判断故障结束时刻来重置IM-SPRT方法检测值; 设计了多余度大气/惯性系统高灵敏度故障检测与隔离框架, 仿真通过对比本文所提方法与GLT方法的检测性能, 验证了本文所提算法的有效性, 结果表明, 本文所提出的算法相比于GLT方法对软故障的检测灵敏度更高, 能有效提高对民机多余度大气/惯性系统的故障检测能力。噪声对该方法性能的影响还有待进一步研究。