古诗三首教学反思第28卷第6期中国惯性技术学报V ol.28 No.6 2020年12月Journal of Chine Inertial Technology Dec. 2020 文章编号：1005-6734(2020)06-0834-07 doi.10.13695/jki.12-1222/o3.2020.06.021
渐消SPRT在SINS/CNS/SRS导航系统软故障诊断中的应用
高广乐1,2，高社生1,2，彭旭3，胡高歌1,2
（1. 西北工业大学自动化学院，西安 710072；2. 西北工业大学深圳研究院，深圳 518057；
3. 南瑞集团（国网电力科学研究院）有限公司，南京 211106）
摘要：为提高导航系统的可靠性，实时的故障诊断与隔离十分必要。序贯概率比检验(SPRT)对传感器中缓慢增长的软故障具有较高的灵敏度，但该方法却存在故障检测延迟甚至无法检测故障结束的缺陷，造成故障的漏警和误警。对SPRT在故障检测中存在的缺陷进行了分析并提出了一种渐消SPRT 方法。渐消SPRT方法通过引入渐消因子，降低了历史信息对故障时刻统计量变化率与统计量的影响，进一步增强SPRT对故障的灵敏度，从而实现降低故障检测延迟并避免SPRT无法检测故障结束的问题。最后基于捷联惯导/天文/光谱红移(SINS/CNS/SRS)组合导航系统，进行了导航系统实施故障诊断仿真验证。仿真结果表明相比于SPRT，渐消SPRT能够检测到缓变故障的结束并将故障开始时刻的检测延迟降低了41%，
极大地提高了SINS/CNS/SRS组合导航系统的实时估计精度与可靠性。在实时导航系统中，渐消SPRT能够起到良好的软故障检测与隔离作用，保证了系统的稳定性。
关键词：故障检测；序贯概率比检验；渐消因子；光谱红移导航；组合导航
中图分类号：V249.3 文献标志码：A
Fading SPRT method for soft fault diagnosis in SINS/CNS/SRS
integrated navigation system
GAO Guangle1,2, GAO Shesheng1,2, PENG Xu3, HU Gaoge1,2
(1. School of Automation, Northwestern Polytechnical University, Xi’an 710072, China; 2. Rearch &
Development Institute of Northwestern Polytechnical University in Shenzhen, Shenzhen 518057, China;长沙的气候
3. NARI Group Corporation (State Grid Electric Power Rearch Institute), Nanjing 211106, China)
Abstract: In order to improve the stability of the navigation system, real-time fault diagnosis and isolation are esntial. The quential probability ratio test (SPRT) has nsitivity to the detection of slow-growing soft faults, but SPRT has an obviously defect that has time delay in fault detection and even cannot detect the fault end. The defect of SPRT is analyzed and further a fading SPRT method is propod. The Fading SPRT method introduces the Fading factor in the calculation of SPRT statistics, which reduces the influence of historical measurement on the statistics at the time with failure, thus avoids the problem that SPRT cannot detect the end of failure and reduce the delay in fault detection. Finally, bad on the propod SINS / CNS / SRS integrated navigation system, the fading SPRT method is verified. As the simulation results shown, fading SPRT can detect the fault end and reduced the time delay in fault start detection by 41%, which greatly improves the real-time estimation accuracy and reliability of SINS/CNS/SRS integrated navigation system under slow-growing fault.
Key words: fault detection; SPRT; fading factor; spectral redshift navigation; integrated navigation
收稿日期：2020-07-04；修回日期：2020-11-14
基金项目：国家自然科学基金项目（42004021）
作者简介：高广乐（1993—），男，博士研究生，从事导航，制导以及控制方向研究。E-mail：*****************：
联 系 人：高社生（1956—），男，教授，博士生导师，从事导航，制导以及控制方向研究。E-mail：*****************
第6期高广乐等：渐消SPRT在SINS/CNS/SRS导航系统软故障诊断中的应用835
为构建一个长期稳定且可靠的组合导航系统，采用多传感器的组合导航方式受到了广泛关注[1-3]。在众多导航方法中，光谱红移导航系统(SRS)是一种具有前瞻性的自主导航方法。光谱红移导航只利用天体光谱的多普勒频移效应实现飞行器导航，具有无时间延迟、选星容易以及自主性强等优点[4,5]。而将SRS与捷联式惯性导航系统(SINS)以及天文导航系统(CNS)进行组合的SINS/CNS/SRS组合导航系统，则不仅能够在天文导航选星困难时，通过SRS辅助纠正捷联惯导系统误差，同时也保证了导航系统的高度自主性。
但是，在恶劣的环境之中，CNS或SRS系统传感器经常会出现故障导致导航系统不可靠，因此在组合导航系统中进行故障的隔离与诊断则十分必要。Carlson等首先提出了联邦卡尔曼滤波理论，使各导航子系统相互分离，降低系统之间相互干扰，一定程度提高了系统容错性，为导航系统故障隔离的实现奠定了基础[6,7]。关于故障诊断方法，卡方检验法目前是工程中最常用的故障诊断方法，对突变故障
有较好的检测效果，但对缓慢增加的软故障，卡方检验法却不够灵敏[8,9]。进而，Joerger M提出了一种累积的卡方检验法，将加权的残差平方累积作为统计量，能够较快地检测出系统的变化，但随着时间的推移，累积量会不断变大且难以下降[10]。为解决累积量难以下降的问题，提出了基于卡尔曼滤波的序贯概率比检验方法(SPRT)，用于检测传感器中存在的缓变故障与持续偏差[11,12]。但由于在传统SPRT检验中，统计量的计算涉及了较多的历史量测信息，因此也不能快速变化，及时反映传感器故障的发生与结束，甚至还会引起无法检测故障结束的问题，造成了系统漏警与误警[13]。在组合导航系统中，长时间的误警会导致正常的CNS 或SRS观测信息被长期隔离而失效，使得SINS系统发散，导航系统精度与可靠性降低。
为弥补SPRT方法的缺陷，本文提出了一种渐消SPRT方法。该方法引入渐消因子，降低历史量测信息对新时刻统计量以及统计量变化率的影响，避免了传统SPRT方法无法检测故障结束的缺陷并降低了故障开始时刻的检测延时，从而减少检测系统的误警与漏警。最后基于SINS/CNS/SRS组合导航系统，对渐消SPRT方法诊断缓慢故障进行了仿真，验证了该方法能够及时检测与隔离导航中存在的传感器软故障，保证了导航系统的估计精度与可靠性。
1 联邦卡尔曼滤波算法
为方便故障传感器隔离，提高导航容错性，在多传感器组合导航系统中可采用联邦卡尔曼滤波器（FKF）[7]。
假设系统方程为
1
k k-
X=FX+W(1)
,,
+
i k i k k i
有亭翼然
=
Z H X V(2) 式中，下标k代表采样时间，X为系统状态向量，F 为系统状态转移矩阵，下标i表示提供量测信息的不同传感器，Z表示量测向量，H为量测矩阵，~(0,)
N
W Q和~(0,)
N
V R分别代表状态噪声矩阵和量测噪声矩阵。
基于以上系统方程，FKF的主要步骤如下：
1）状态预测
/11
ˆˆ
k k k k
--
=
X F X(3)
T
/
11
k k k k k
--
=+
P F P F Q(4)
式中，
1
ˆ
k-
X为k-1时刻状态估计，
/1
ˆ
k k-
X为k时刻状态
预测，/1
k k-
P为预测协方差阵。
2）子滤波器的量测更新
T T1
,/1,,1,
()
i k k k i k i k k i k i
-
--
=+
K P H H P H R (5) ,/1,,,/1
ˆˆˆ
()
i k k k i k i k i k k k
--
=+-
X X K Z H X (6)
,,,/1
()
i k i k i k k k-
=-
P K H P
I(7) 式中，,i k
K为i子滤波器的增益阵；
,
ˆ
i k
X为i子滤波器下的状态估计；,i k
P为i子滤波器下状态协方差阵估计。
3）多个子滤波器信息融合
11
,
1
()
n
k i k
i
--
=
=∑
P P(8)
1
,,
1
ˆˆ
()
n
k k i k i k
i
-
=
=∑
X P P X(9)
式中，k P为融合后的协方差阵估计；ˆ
k
X为融合后的
状态估计。
2 基于卡尔曼滤波的渐消序贯概率比检验
2.1 基于卡尔曼滤波的SPRT故障检测的基本原理
设故障下的导航子滤波器的量测方程可表示为
,
++
i k i k k
=
Z H X Vρ(10) 其中
,
,
=
0,
i k
i k
⎧
⎨
⎩
ρ
ρ
，有故障
无故障
(11) 因此导航子系统的新息值可表示为
,,,
ˆ
()=+
i i k i k/k-1i k i k
k=-
νZ H X Vρ(12) 选取k时刻内的导航子滤波器下的所有新息值组成新
息序列{()|1,2}
i
j j k
=
ν 。设系统状态方程无扰动，且子系统量测量中无故障时，该子滤波器中的新息序列应相互独立且服从正态分布(0,)
NΣ[13]。此时新息期望值和方差可表示为
T
E[]0E[]
i i i
==
νννΣ
陈三愿
，(13)
反问和设问有什么区别836 中国惯性技术学报 第28卷
而当系统有故障时，可得新息期望与方差为
T
ˆˆˆE[]0E[()()]i i i i i i =≠--=νρ
νρνρΣ， (14) 根据式(13)(14)，可设置如下假设检验：
0H :ˆ0i =ρ
，系统正常稳定；1H :ˆ0i ≠ρ，则系统出现未知故障。且知两种假设下新息样本的概率密度为
T 11
()()()20/21/2
1
(()|H )2()
i i k k i n p k e π--=νΣννΣ (15) T 1
1ˆˆ()()())
2
1/21/2
1
(()|H )2()i i i i k k i n p k e
π----=
νρΣνρνΣ（）（）(16)
其中
1
1
ˆ()k
i i i j j k =≈=∑ρ
νν (17) T
,/1,/1,=zz i k k i k k k i k i --=+ΣP H P H R (18)
进一步，根据两类假设的概率密度函数，则可计算两类假设下的似然比为
11100((1),(2),,()|H )(()|H )
()((1),(2),,()|H )(()|H )
k
i i i i j i i i i p k p k L k p k p k ==
=∏ννννννννΛΛ(19)
再对()L k 求对数，所得变量即为基于卡尔曼滤波的
SPRT 检验方法的统计量
110T 1T 1,/1/1j 1
(()|H )柜子英文
()In()(()|H )(())()(())()()()
2
k
i
j i
zz
zz k i i i k k i i i k k i p k k p k k k k k λ=----===---∏∑
ννννP νννP ν(20)
设每时刻的统计量变化率可表示为
T 1T 1
,/1,/1()()()(())()(())()=22
zz zz i i j j i i i i j j i i j j j j j λ------∆-
νP νννP νν(21) 则可获取SPRT 统计量的迭代形式为
()(-1)+()k k k λλλ=∆ (22)
传统的Wald 提出的SPRT 方法中，常采用双阈值进行故障判断[12]，其双阈值设置如下
up =In(
1m f P T P -，down 1 =In(
m
f
P T P - (23) 其中，f P 表示误警率；m P 表示漏警率。
最后，根据阈值和统计量，可写出基于卡尔曼滤波的SPRT 故障诊断方法的检验函数为
up down up down ,
() (),(),()k T e k T k T k T λλλ>⎧⎪
=⎨⎪<⎩故障继续检测正常 (24)
由式(24)中的判断函数可以看出，双阈值法在统计量居于两个阈值中间时，则该时刻系统不做出判断。而实时性系统，如导航系统，则需要在每个时刻都知道传感器量测信息是否可被引入式(5)~(7)，进行量测更新。因此双阈值法完全不适用实时导航系统。为避免双阈值法导致无法得到实时检测结果的问题，可对SPRT 法进行改进，只采用单阈值进行检测。
除此之外，在卡尔曼滤波过程中，误警会导致量测失效，使系统只进行状态预测。而短时间进行状态预测不会造成过多的状态量偏差。并且若存在其他未发生故障传感器及时纠正偏差，则系统仍可获取准确的导航信息。相反，漏警会造成存在故障的量测值被引入系统导致误差更大的状态估计值。随着
卡尔曼滤波过程迭代，故障引起的偏差会在状态量中累积，使得状态量跟随故障，导致系统更无法检测出故障而产生更多的漏警。因此在采用单阈值法时，系统采用较小阈值即
down T ，进行单阈值检测来避免系统产生过多漏警。相
应的，SPRT 故障诊断检验函数修改为
down down ,
() (),()k T e k k T λλ⎧=⎨<⎩故障正常 (25)
2.2 基于卡尔曼滤波的SPRT 故障检测的缺陷分析
阈值改进后的传统SPRT 方法对缓慢增加故障的检测灵敏度较高，但仍存在延迟和无法检测故障结束的缺陷。为解决此问题，首先需要对缺陷出现的原因进行分析，具体如下：
SPRT 统计量可表示为
=1
()()k
j k j λλ=∆∑ (26)
且()j λ∆可变换为
T 1
,/1T 1
,/1T 1T 1
,/1,/1T 1T 1
,/1,/1T ()()()
()=
2
(())()(())
=
2
()()()()()()1trace =()+2()()21trace 2()2
zz i i j j i zz i i i j j i i zz zz i i i j j i i i j j zz zz i i i j j i i i j j i i i j j j j j j j j j j j λ------------∆---⎡⎤-⎢⎥-⎢⎥⎣⎦
-νP νννP ννννP ννP ννP ννP νννT 1
,/1T 1
,/1()=1trace ()(2())2
zz i i j j zz i i j j i i j ----⎡⎤⎣⎦⎡⎤-⎣⎦νP νP νν(27) 同时将i ν也写作迭代形式为
1111
()()=(1)+()k i i i j k k j k v k k k k
=-=-∑ννν (28)
当量测系统长时间正常时，0i →ν。进而由式 (27)(28)可知，此时也存在()0j λ∆→,()0k λ→。
假设在t 1时刻系统发生故障。当t 1较大时，由式 (27)(28)可看出，在故障初期(1)i t ν仍接近于零，导致()j λ∆增长缓慢，使得统计量需要累积较长时间才能够检测到故障开始，造成较长时间漏警。
进一步，假设已经被检测到的故障在t 2时刻结束。因此，(21)t λ-应已超过阈值，且故障越大，(21)t λ-越大。而当故障结束后，统计量需要较长的时间累计此时为负值的统计变化量()j λ∆，使自身下降至阈值以下，实现对故障结束的检测。但是，在故障结束之后的时间里，量测系统也已恢复正常，使得i ν逐渐趋于0，
中统计量未下降至阈值之下时，趋于零的统计量变化率则很难再使统计量快速降低，从而导致在故障结束后统计量一直大于阈值，出现无法检测故障结束的问题，进而引起系统持续的误警，导致在卡尔曼滤波过程中，提供量测量的传感器彻底失效。
2.3 基于卡尔曼滤波的渐消SPRT 检验方法基本原理
基于以上分析可看出，SPRT 无法及时检测故障开始而导致漏警的原因是由于故障开始时，正常时刻
幸福课的新息仍对i ν的计算产生较大的影响；而无法检测故障结束的主要原因是由于故障时刻的历史统计量在故障结束后，仍然对统计量计算产生影响。
为解决SPRT 的缺陷，本文对SPRT 方法进行了改进：可在计算新息均值i ν和统计量λ时引入渐消因子，对SPRT 中历史统计量以及历史新息采取渐消处理，在故障时刻降低历史信息的影响，提高统计量对故障的敏感度。
首先，为增强均值在故障发生时刻的变化率，减少故障发生时刻检测延迟，可在均值计算中引人渐消技术，降低历史时刻新息在均值计算中的权重，从而增强当前时刻新息样本的影响。改进后公式如下： 11
()=(1)+()i i k k
k k k k
ααα
--ννv (29)
式中，α为渐消因子，文中选取0.85。
由式(29)可看出，当量测系统无故障时，仍存在0i →ν。而当量测系统发生故障时，渐消因子的引入
加强了新时刻的新息量的占比，使得在故障开始时信息均值对故障的发生更为敏感，从而可减少故障开始检测的延迟并降低故障的漏警。贫困
同时，为及时检测故障结束，可改进统计量累计方式为：
-1
1
() ()()(1)()k j k j k k k λαλλαλλ==∆+∆=⋅-+∆∑ (30)
由式(27)可看出，采用渐消因子后，随着时间的推移，时间越久的历史时刻的λ∆对统计量()k λ的影响会越来越小。因此在故障结束后，随着时间推移，即使新时刻的()0j λ∆→，受渐消因子影响故障时刻
的历史统计量变化率对故障结束后统计量的影响也会逐渐消失，使得统计量快速下降为零，从而可避免统计量持续大于阈值导致故障结束无法被检测的情况出现，进一步也消除了此情况所带来的持续误警。
3 SINS/CNS/SRS 组合导航系统
3.1 光谱红移导航原理
根据光谱红移原理与多普勒公式可知以下关系
(31)
式中，z 表示飞行器上实拍光谱中计算出的红移值，p v 表示惯性系中飞行器与光源之间相对速度矢量，c 表示
真空中的光速，
θ表示惯性坐标系中光源和飞行器之间连线与视向速度矢量r v 的绝对夹角。
假设飞行器在同一时间可接收光源A ，光源B 和光源C 的光谱信息，可计算出每时刻三条光谱相对于静止光谱的红移值123,,z z z 。将红移值带入式(31)后，再结合以下运算关系
1p 11
2p 223p 33
()()()r r r v v v ⎧=-⋅⎪
=-⋅⎨⎪
=-⋅⎩v v u v v u v v u (32) 可得
p 11
1p 222p 333()(10()(10()(10z z z ⎧-⋅-+=⎪⎪⎪-⋅-+=⎨⎪
⎪-⋅-+=⎪⎩
v v u v v u v v u (33)
式中，123,,v v v 分别表示三个参考天体在惯性系中的速
度矢量，可以通过查询相关天体星历获取。123,,u u u 分别表示为惯性系中各天体指向飞行器位置矢量的单位矢量，可由星敏感器测得。p =cos r v θv 为惯性飞行器沿着光源方向的径向速率。
为简便计算可将式(33)线性化。设
p p ()()(1ce ce ce ce Z z =-⋅-+v v v u (34)
式中，ce 表示不同的观测天体。
再将式(34)在=0p v 处一阶泰勒展开后得
p p p p
p p =0=0=0=0()=()|
|||ce ce ce ce px py px py ce
pz Z
pz Z Z Z Z
v v v v Z v v +∂∂⋅+⋅+∂∂∂⋅+∆∂v v v v v v (35) 式中，Z ∆代表高阶项。
最后，省略式(35)中高阶项后，联立方程组可得
p p p p p p p p p p p p 1111p =0=0=0=03222p =0=0=0=03333p =0=0=0=0()||||0()||||0()||||0
px py pz px py pz px py pz px py pz
px py pz px py pz Z Z Z
Z v v v v v v Z Z Z Z v v v v v v Z Z Z
Z v v v v v v ⎧∂∂∂+⋅+⋅+⋅=⎪∂∂∂⎪
⎪∂∂∂⎪
+⋅+⋅+⋅=⎨∂∂∂⎪∂∂∂+⋅+⋅+⋅=∂∂∂⎩
v v v v v v v v v v v v v v v ⎪⎪⎪(36) 最后通过求解式(36)即可获取p v 值。
3.2 惯性/天文/光谱红移组合导航系统
考虑到光谱红移导航优点，本文提出了一种惯性/
838 中国惯性技术学报
第28卷
天文/光谱(SINS/CNS/SRS)红移组合导航系统。在天文导航失效时，惯性/天文/光谱红移组合导航系统可及时通过SRS 求解飞行器速度信息，及时纠正SINS 的误差，保证组合导航系统的自主性与可靠性。SINS/CNS/SRS 组合导航系统结构如图1所示。
图1 SINS/CNS/SRS 组合导航系统结构图 Fig.1 Structure of SINS/CNS/SRS integrated
navigation system
3.2.1 SINS/CNS/SRS 导航系统状态方程
根据惯性导航系统误差模型，我们可给出组合导航系统状态方程为[4]
1k k k -X =F X +W (37)
式中，F 为惯导系统的状态转移矩阵，X 为惯导系统的状态向量，具体为 E N H E
N H b
b b
T
X
Y Z X Y Z =[ ]v v v L h
φφφδδδδδλδεεε∇∇∇X (38) 式中，E N H (,,)φφφ表示导航系下平台失准角；E N H (,,)v v v δδδ表示导航系下速度误差；(,,)L h δδλδ表
示导航系下位置误差；b b b X
Y Z (,,)εεε和X Y Z (,,)∇∇∇分别表示陀螺随机游走以及加速度计常值误差。
W 为系统噪声矩阵，具体为
a T E
N
H
E
N
H []g
g g a a w w w w w w =W (39)
式中，g g g
E N H w w w （，，）分别表示陀螺仪的随机漂移。a a a
E N H w w w （，，）表示加速度计误差。
3.2.2 SINS/CNS/SRS 组合导航系统量测方程
SINS/SRS 导航子系统的量测值为捷联惯性导航系统解算的速度信息与光谱红移导航系统解算速度之间的差值。
由于红移导航速度计算都是在惯性系下完成，因此建立量测方程前首先需要进行坐标转换，转换公式如下
i i e
SRS g p e g p SRS +==v C v C C v V (40)
式中，i e C 表示地心惯性系到地球坐标系的转换矩
阵，
e g
C 表示地球坐标系到地理坐标系的转换矩阵，SRS v 表示地理坐标系下光谱红移速度量测值,SRS V
表示量测噪声。
SINS/SRS 导航子系统测量方程可表示为
v SINS SRS SRS =E N H v v v δδδ⎡⎤
⎢⎥=-+⎢⎥⎣⎦
Z v v V (41)
SINS/CNS 导航子系统经纬度量测方程为SINS 系统输出的经纬度位置信息与CNS 系统算出的经纬度信息的差值，具体如下
SINS CNS p p CNS SINS CNS =k L L λλ-⎡⎤
=+⎢⎥-⎣⎦
Z H X V (42)
式中，CNS CNS (,)L λ为天文导航系统输出的经纬度信息；SINS SINS (,)L λ为惯性导航系统输出的经纬度信息；CNS
V 为天文测距量测噪声，p H 为量测矩阵。
天文导航系统无法获取高度信息，为了防止SINS 高度通道发散，引入气压高度表测量高度信息，由气压
高度表和SINS 输出的高度之差作为量测量，得高度差的量测方程为
[]SINS BA BA h h k h h =-=+Z H X V (43)
式中，
SINS h 和BA h 分别为惯导和气压高度表输出的高度信息，h H 为高度量测阵，BA V 为高度量测噪声阵。
4 仿真验证
基于SINS/CNS/SRS 组合导航系统对渐消SPRT 方法进行软故障诊断仿真验证。导航系选取东北天坐标系，飞行器飞行轨迹如图2所示。
110.5
34.246经度/(°)
纬度/(°)
高度 /m
图2 飞行器轨迹图
Fig.2 Flight trajectory
设置飞行总时长为1200 s ，飞行器初始位置为
(108.90934.24610000m)︒︒，，。常规情况下传感器参数如表1所示。
表1 传感器仿真参数设置
Tab.1 Parameter of nsors in simulation 传感器 误差参数 工作周期陀螺仪 0.5(/h)︒(Constant drift) 0.05 s 0.5(/h)︒(STD of white noi)
加速度计 0.1m g (Constant bias) 0.05 s 0.1m g (STD of white noi)
CNS 位置 20 m 1 s SRS 速度 0.5 (m/s) 1 s 气压高度计
10 m
0.1 s
同时设置缓变故障模型数学模型为
