现代大型船舶、飞机对导航系统的性能要求越来越高，单一导航系统已不能满足大型船舶、飞机对导航精度和可靠性的要求，由惯性导航系统(inertial navigation system，INS)、全球导航卫星系统(global navigation satellite system，GNSS)、天文导航系统(celestial navigation system，CNS)等若干子系统相结合构成的组合导航系统，可充分利用各导航子系统特点，实现优势互补，极大地提高了导航系统的精度和可靠性，已经成为大型船舶、飞机广泛采用的导航系统^[1]。

组合导航系统通常采用联邦滤波作为基本的信息融合算法^[2]。但是，联邦滤波算法的自身特性决定了当子系统中的一个或几个发生故障时，故障数据会污染整个组合导航系统，从而导致组合导航系统的性能严重下降。为了避免这种情况，需要对组合导航系统及其各子系统的故障进行实时有效的检测及诊断^[3-4]。根据故障的时间特性，组合导航系统的故障可分为突变故障和缓变故障^[5]。组合导航系统常用的故障检测方法主要有χ²检验法^[6-12]、广义似然比检测^[13-14]、小波分析法^[15-16]、基于神经网络的方法等^{[7, 17-18]}。针对基于联邦滤波的组合导航系统的特点, 常采用的方法是χ²检验法，该方法不需要确定造成故障的原因。χ²检验法又可分为状态χ²检验法和残差χ²检验法^[2]两类。

状态χ²检验法中，卡尔曼滤波器提供的状态与量测相关，会受到故障量测影响；“状态递推器”^[2]的状态只与设定的初值相关，不受故障量测影响。通过对二者的差异进行奇偶校验，便可以实现故障的检测和隔离。但“状态递推器”不引入量测更新，系统误差将随着时间积累，进而导致状态χ²检验的检测灵敏度下降直至失效。文献[4]通过2个辅助滤波器定期交换校正，解决了这一问题，但增加了计算机负担，不适用于实时性的组合导航系统。残差χ²检测法中，残差定义为系统实际测量值与预测值差值。通过构造残差的检测函数，对系统故障进行实施有效的检测。与状态χ²检验法比较，残差χ²检测法算法简单，实时性更好，更适用于动态系统。但是，残差χ²检验法进行状态更新时引入了量测，当系统故障时，状态会受到故障量测影响，无法对缓变故障进行实时有效的检测。针对这一问题，文献[5]引入模糊逻辑和加权平均计算检测阈值，通过突变信号大小自动调节检测时间长短，改善对突变故障的检测灵敏度；文献[6-7]将自适应滤波和残差检测法相结合，改善了对缓变故障的检测灵敏度，但要求系统量测噪声阵确切已知；文献[8]联合残差χ²算法和改进序贯比算法，弥补了残差χ²算法故障检测灵敏度低和改进序贯比算法无法判断故障结束时间的缺点；文献[9]利用状态χ²检验法和残差χ²检验法两者相互配合，提高对系统软硬故障的检测精度，但难以很好的兼顾实时性。

本文从组合导航系统的物理特征出发，采用时间序列分析的方法，通过建立组合导航系统残差自回归模型，采用实时数据进行残差变化分析，提高残差χ²检验法对组合导航系统缓变故障检测灵敏度和准确性。

1 组合导航系统残差的时间序列分析 1.1 组合导航系统的残差

组合导航系统常用带故障的离散系统动态模型^[10]：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{X}}\left( k \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {k,k - 1} \right)\mathit{\boldsymbol{X}}\left( {k - 1} \right) + \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {k - 1} \right)\mathit{\boldsymbol{W}}\left( {k - 1} \right)\\ \mathit{\boldsymbol{Z}}\left( k \right) = \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{X}}\left( k \right) + \mathit{\boldsymbol{V}}\left( k \right) + f\left( {k,\varphi } \right)\mathit{\boldsymbol{\gamma }},k \ge 0 \end{array} \right. $

(1)

式中：Z(k)∈R^m是系统的量测值；X(k)∈Rⁿ是系统状态；Φ(k, k-1)∈R^m×n是系统状态转移矩阵；Γ(k)∈R^n×r是系统噪声矩阵；W(k)∈R^r和W(k)∈R^m是相互独立的高斯白噪声序列，且:

$ \left\{ \begin{array}{l} {\rm{E}}\left[ {\mathit{\boldsymbol{W}}\left( k \right){\mathit{\boldsymbol{W}}^{\rm{T}}}\left( j \right)} \right] = \mathit{\boldsymbol{Q}}\left( k \right)\delta \left( {k,j} \right)\\ {\rm{E}}\left[ {\mathit{\boldsymbol{V}}\left( k \right){\mathit{\boldsymbol{V}}^{\rm{T}}}\left( j \right)} \right] = \mathit{\boldsymbol{R}}\left( k \right)\delta \left( {k,j} \right)\\ {\rm{E}}\left[ {{\mathit{\boldsymbol{W}}_k}} \right] = 0\\ {\rm{E}}\left[ {{\mathit{\boldsymbol{V}}_k}} \right] = 0 \end{array} \right. $

式中：δ(k, j)为克朗尼克函数；γ是随机向量，它表示故障的大小；Q(k)是系统噪声序列的方差阵，假定为非负定；R(k)是量测噪声序列的方差阵，假定为非负定；f(k, φ)是分段函数：

$ f\left( {k,\varphi } \right) = \left\{ \begin{array}{l} 1,k \ge \varphi \\ 0,k < \varphi \end{array} \right. $

式中φ是故障发生的时间。

k时刻系统状态递推值${\hat{\mathit{\boldsymbol{X}}}}\left( {k,\;k - 1} \right)$和量测预报值${\hat{\mathit{\boldsymbol{Z}}}}\left( {k,\;k - 1} \right)$计算公式分别如下^[9]

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\hat X}}\left( {k,k - 1} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {k,k - 1} \right)\mathit{\boldsymbol{\hat X}}\left( {k - 1} \right)\\ \mathit{\boldsymbol{\hat Z}}\left( {k,k - 1} \right) = \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{\hat X}}\left( {k,k - 1} \right)\\ \mathit{\boldsymbol{\gamma }}\left( k \right) = \mathit{\boldsymbol{Z}}\left( k \right) - \mathit{\boldsymbol{\hat Z}}\left( {k,k - 1} \right) \end{array} \right. $

(2)

式中γ(k)在卡尔曼滤波器中称为新息，即残差。

1.2 残差的AR建模^[21]

组合导航系统的残差信号中的加性误差会随着时间随机变化，当前时刻残差信号与该时刻之前的残差信号并不完全独立。通过分析观测数据得到组合导航系统残差特性，并建立组合导航系统残差模型。

以SINS(strapdown inertial navigation system，SINS)/GNSS组合导航系统为例，SINS/GNSS组合导航系统选取SINS输出的位置信息与GNSS输出的对应位置信息相减作为量测量，SINS/GNSS组合导航系统残差为系统量测与系统预测值的差值。

连续采集SINS/GNSS组合导航系统静止和运动时的残差数据，高斯投影后转换为平面直角坐标系下的残差数据(x, y, z)，证明无论是静态还是动态条件下，残差序列均具有明显的自相关性。

时间序列是随时间改变而随机变化的信号序列，时间序列分析的意义在于研究某一时间序列在长期变动过程中所存在的统计规律性，建立序列近似的、简化的数学模型，并将其应用于系统动态特性的描述、预测分析和误差补偿等方面^[11]。根据SINS/GNSS组合导航系统残差特性，采用自回归模型(auto regressive model，AR)来描述和预测SINS/GNSS组合导航系统残差信号的变化规律。

对于时间序列{X(t), t=0, ±1, ±2, …}而言，AR(p)模型为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_t} = {\varphi _1}{\mathit{\boldsymbol{x}}_{t - 1}} + {\varphi _2}{\mathit{\boldsymbol{x}}_{t - 2}} + \cdots + {\varphi _p}{\mathit{\boldsymbol{x}}_{t - p}} + {\varepsilon _t}}\\ {\left( {t = p + 1,p + 2, \cdots ,p + n} \right)} \end{array} $

(3)

式中：φ₁, φ₂, …, φ_p为自回归系数；p为AR模型的阶数; ε_t为均值为0，方差为σ²的白噪声。自回归系数φ由AR(p)序列的自协方差函数[γ₀, γ₁, …, γ_p]通过Yule-Walker方程唯一决定：

$ \left[ {\begin{array}{*{20}{c}} {{\gamma _1}}\\ {{\gamma _2}}\\ \cdots \\ {{\gamma _p}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\gamma _0}}&{{\gamma _1}}& \cdots &{{\gamma _{p - 1}}}\\ {{\gamma _1}}&{{\gamma _0}}& \cdots &{{\gamma _{p - 2}}}\\ \cdots&\cdots&\cdots&\cdots \\ {{\gamma _{p - 1}}}&{{\gamma _{p - 2}}}& \cdots &{{\gamma _0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\varphi _1}}\\ {{\varphi _2}}\\ \cdots \\ {{\varphi _p}} \end{array}} \right] $

(4)

白噪声的方差σ²由下式求解：

$ {\sigma ^2} = {\gamma _0} - \left( {{\varphi _1}{\gamma _1} + {\varphi _2}{\gamma _2} + \cdots + {\varphi _p}{\gamma _p}} \right) $

(5)

参数估计是在给定阶次的情况下进行的。由于事先无法判断模型的阶次，因此在建模过程中先给定模型的某个阶次，然后按照上述估计方法，估计出AR模型的参数，得到各阶次模型，最后通过Akaike^[19]信息准则定阶来确定AR模型。

1.3 AR模型的检验

当安装SINS/GNSS组合导航系统的载体静止时，选择p阶模型，采用一步预测的方法对误差进行预测

$ {{\mathit{\boldsymbol{\hat x}}}_t} = {\varphi _1}{\mathit{\boldsymbol{x}}_{t - 1}} + {\varphi _2}{\mathit{\boldsymbol{x}}_{t - 2}} + \cdots + {\varphi _p}{\mathit{\boldsymbol{x}}_{t - p}}\left( {t = p,11, \cdots ,n} \right) $

$ {\mathit{\boldsymbol{\omega }}_t} = {\mathit{\boldsymbol{x}}_t} - {{\mathit{\boldsymbol{\hat x}}}_t} = {\mathit{\boldsymbol{x}}_t} - \left( {{\varphi _1}{\mathit{\boldsymbol{x}}_{t - 1}} + {\varphi _2}{\mathit{\boldsymbol{x}}_{t - 2}} + \cdots + {\varphi _p}{\mathit{\boldsymbol{x}}_{t - p}}} \right) $

式中：${{\hat{\mathit{\boldsymbol{X}}}}_t}$为误差的一步预测值, ω_t为预测误差和实际误差的差值。将实际量测与预测误差相减所得到的系统残差信息将更为准确。

定义真实值为载体静止时在系统多次测量的算术平均值，建立系统残差AR模型。通过模型得到的预测值对实际测量的残差数据进行修正，并分析其自相关性。分析表明，修正后的残差数据测量噪声的自相关性明显变弱，即所建立的SINS/GNSS组合导航系残差AR模型能较好的描述系统动态特性，具有适用性。

2 基于残差的时间序列建模的χ²检验 2.1 组合导航系统模型

SINS/GNSS组合导航系统状态方程：

$ \mathit{\boldsymbol{\dot X}} = \mathit{\boldsymbol{F}}\left( t \right)\mathit{\boldsymbol{X}}\left( t \right) + \mathit{\boldsymbol{G}}\left( t \right)\mathit{\boldsymbol{W}}\left( t \right) $

(6)

其中

$ \begin{array}{l} \mathit{\boldsymbol{X}}\left( t \right) = \left[ {{\phi _E}\;{\phi _N}\;{\phi _J}\;\phi {\nu _E}\;\phi {\nu _N}\;\phi {\nu _U}\;\delta L\;\delta \lambda \;\delta h\;{\varepsilon _x}\;{\varepsilon _y}} \right.\\ \;\;\;\;\;\;\;\;\;\;\;{\left. {{\varepsilon _z}\;{\nabla _x}\;{\nabla _y}\;{\nabla _z}} \right]^{\rm{T}}}, \end{array} $

式中：ϕ_E、ϕ_N、ϕ_U为数学平台失准角，δν_E、ϕδν_N、δν_U分别为载体的东向、北向和天向速度误差，δL、δλ、δh分别为纬度误差、经度误差、高度误差，ε_x、ε_y、ε_z、∇_x、∇_y、∇_z分别为陀螺随机常值漂移和加速度随机常值零偏。

系统的噪声转移矩阵G(t)为

$ \mathit{\boldsymbol{G}}\left( t \right) = {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_b^m}&{{{\bf{0}}_{3 \times 3}}}\\ {{{\bf{0}}_{3 \times 3}}}&{\mathit{\boldsymbol{C}}_b^n}\\ {{{\bf{0}}_{9 \times 3}}}&{{{\bf{0}}_{9 \times 3}}} \end{array}} \right]_{15 \times 6}} $

系统噪声矢量由陀螺仪和加速度计的随机误差组成，表达式为

$ \mathit{\boldsymbol{w}} = {\left[ {\begin{array}{*{20}{c}} {{\omega _{{\varepsilon _x}}}}&{{\omega _{{\varepsilon _y}}}}&{{\omega _{{\varepsilon _z}}}}&{{\omega _{{\nabla _x}}}}&{{\omega _{{\nabla _y}}}}&{{\omega _{{\nabla _z}}}} \end{array}} \right]^{\rm{T}}} $

系统的状态转移矩阵F(t)为

$ \mathit{\boldsymbol{F}}\left( t \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{F}}_N}}&{{\mathit{\boldsymbol{F}}_S}}\\ {{{\bf{0}}_{6 \times 9}}}&{{\mathit{\boldsymbol{F}}_M}} \end{array}} \right] $

式中：F_N为对应的9维基本导航参数系统阵^[21]，F_S和F_M分别为

$ {\mathit{\boldsymbol{F}}_S} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_b^n}&{{{\bf{0}}_{3 \times 3}}}\\ {{{\bf{0}}_{3 \times 3}}}&{\mathit{\boldsymbol{C}}_b^n}\\ {{{\bf{0}}_{3 \times 3}}}&{{{\bf{0}}_{3 \times 3}}} \end{array}} \right] $

$ {\mathit{\boldsymbol{F}}_M} = \left[ {{\mathit{\boldsymbol{0}}_{6 \times 6}}} \right] $

SINS/GNSS组合导航系统采用SINS与GNSS的位置之差作为卡尔曼滤波的量测信息，得到量测方程为

$ \mathit{\boldsymbol{Z}}\left( t \right) = \left[ \begin{array}{l} {L_{IE}} - {L_{GE}}\\ {\lambda _{IN}} - {\lambda _{GN}}\\ {h_{IU}} - {h_{GU}} \end{array} \right] = \mathit{\boldsymbol{HX}}\left( t \right) + \mathit{\boldsymbol{V}}\left( t \right) $

(7)

其中

$ \mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {{{\bf{0}}_{3 \times 6}}}&{{\rm{diag}}\left( {{R_{\rm{M}}},{R_{\rm{N}}},\cos L,1} \right)}&{{{\bf{0}}_{3 \times 6}}} \end{array}} \right] $

$ {\mathit{\boldsymbol{V}}_p} = \left[ {\begin{array}{*{20}{c}} {{N_{{\rm{GE}}}}}&{{N_{{\rm{GN}}}}}&{{N_{{\rm{GU}}}}} \end{array}} \right] $

式中：下标I表示SINS，下标G表示GNSS；N_GE、N_GN、N_GU为GNSS沿东北天方向的位置误差；R_M为地球子午平面内的主曲率半径，R_N为与子午面垂直且共法线的椭圆主曲率半径，L为地心纬度。

2.2 组合导航系统残差的χ²检验

无故障发生时系统卡尔曼滤波器的残差γ(k)是零均值高斯白噪声，其方差为

$ \mathit{\boldsymbol{A}}\left( k \right) = \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{P}}\left( {k,k - 1} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right) + \mathit{\boldsymbol{R}}\left( k \right) $

(8)

定义系统故障检测函数为

$ \mathit{\boldsymbol{\lambda }}\left( k \right) = {\mathit{\boldsymbol{\gamma }}^{\rm{T}}}\left( k \right){\mathit{\boldsymbol{A}}^{ - 1}}\left( k \right)\mathit{\boldsymbol{\gamma }}\left( k \right) $

(9)

λ(k)服从自由度为n的χ²分布，即λ(k)~χ²(n)。因此判断准则为

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\lambda }}\left( k \right) \ge {T_D},有故障\\ \mathit{\boldsymbol{\lambda }}\left( k \right) < {T_D},无故障 \end{array} \right. $

(10)

式中T_D是预先设置的门限，它决定了故障检测的性能。由奈曼皮尔逊准则可知，当限定误警率P_f=α时，则P_f=P[λ_K > T_D/H₀]解出的门限T_D可使漏检率P_f=P[λ_K≤T_D/H₁]最小。因此T_D可由误警率P_f确定。

2.3 基于时间序列模型的残差χ²检验

根据组合导航系统的物理特征，由式(3)中组合导航系统残差AR模型得到k时刻系统基于AR模型的残差，其方差为

$ {\mathit{\boldsymbol{A}}_a}\left( k \right) = \mathit{\boldsymbol{H}}\left( k \right){\mathit{\boldsymbol{P}}_a}\left( {k,k - 1} \right){\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right) + \mathit{\boldsymbol{R}}\left( k \right) $

(11)

定义系统故障检测函数为

$ {\mathit{\boldsymbol{\lambda }}_a}\left( k \right) = \mathit{\boldsymbol{\gamma }}_a^{\rm{T}}\left( k \right)\mathit{\boldsymbol{A}}_a^{ - 1}\left( k \right){\mathit{\boldsymbol{\gamma }}_a}\left( k \right) $

(12)

λ_a(k)服从自由度为n的χ²分布，即λ_a(k)~χ²(n)。其判断准则与门限T_D取值方法与式(10)所示相同，在此不加累述。

3 基于时间序列模型残差χ²检验的仿真分析

根据时间序列模型残差χ²检验模型对系统进行仿真，仿真条件设置捷联惯导SINS的输出频率为50 Hz，GPS的更新频率为1 Hz。陀螺的零偏稳定性为0.1 (°)/h，加速度计零偏稳定性为10^-4 g。载体初始位置为北纬32°，东经118°。初始姿态角为0°，初始速度30 m/s，载体进行匀速航行。仿真中分别利用残差χ²检验和基于时间序列模型残差χ²检验方法对组合导航系统进行数据处理，时间序列建模点迹数目选择为3 000个，仿真时间为300 s，在100~200 s加入故障。故障的设置条件为100~200 s，在经度量测上加入0.05(t-100) m的缓变故障。

图 1、2中是故障模式下传统残差χ²检验和基于时间序列模型的残差χ²检验的残差数据比较。从图中可以看出两者在无故障时均接近于0。但当GPS在100~200 s出现失锁现象后，传统残差χ²检验的残差数据因受到量测污染，变化很小，始终无法达到检测阈值，因而无法及时检测出GPS故障，出现了漏检的现象，当故障结束时，近似于产生了一个突变，此时会认为有故障发生。而基于时间序列模型的残差χ²检验方法的残差数据则随着故障引入逐渐增大，在故障发生45 s后达到检测阈值，对缓变故障做出告警。由此可见，对比传统残差χ²检验方法，基于时间序列模型的残差χ²检验方法，对于组合导航系统软故障检测具有更高的灵敏度和准确性。

4 结论

1) 基于时间序列模型的残差χ²检验法，能够实时有效地检测组合导航缓变故障。

2) 提高了故障检测的灵敏度和准确性。改善了组合导航系统的整体性能。