利用时频分析法精化GNSS/INS组合导航先验滤波模型*

甘 雨 隋立芬 张清华 桑 渤

1)信息工程大学地理空间信息学院，郑州 450052

2)72946 部队，淄博255000

1 引言

GNSS/INS 组合导航应用Kalman 滤波进行数据融合，经典Kalman 滤波要求有可靠的函数模型及随机模型。先验滤波模型的确定不仅困难，而且有时并不能完全合理地描述整个导航过程，因此许多学者研究了多种自适应滤波理论［5－6］，并对经典Kalman 滤波进行了扩展，解决了许多实际应用问题，特别是抗差自适应滤波，可以同时控制状态扰动和观测异常的影响。

然而，应用Kalman 滤波时，先验信息的确定往往具有一定的随意性和较强的经验性，缺乏理论基础。如果先验模型与真实模型相差较大，不仅影响经典Kalman 滤波，也会降低自适应滤波的可靠性。在组合导航系统中，状态随机模型的确定是难点。现有文献极少涉及先验模型确定问题，影响了Kalman 滤波在组合导航中的实际应用。

组合导航系统本质上是连续线性系统，具有功率谱密度形式的频域随机模型。提出利用功率谱密度分析、自相关函数分析、Allan 方差分析等时频特性分析方法，由惯性元件数据计算各状态噪声成分的功率谱密度，为组合导航提供较为准确的先验信息，避免实际应用Kalman 滤波时繁琐的调试。算例验证了时频分析精化先验滤波模型的有效性。

2 GNSS/INS 组合导航的随机模型

假设线性系统用微分方程可表示为

式中，X(t)为系统状态，u(t)为系统驱动白噪声过程，F(t)和G(t)为系数阵。u(t)满足:

q 为u(t)的功率谱密度阵，又称方差强度阵。问题的关键是求取状态噪声谱密度阵q。

陀螺误差可表示为［7，8］:

其中，εc为非白噪声部分，wε为白噪声部分。εc可以描述为随机常值、高斯马尔可夫过程、角速率随机游走或者它们的组合形式。

同理，加速度误差▽可表示为:

其中，▽c为非白噪声部分，w▽为白噪声部分。

设εc和▽c的驱动白噪声分别为ηε和η▽，则系统的状态噪声u 为:

由此可知，精化组合导航的状态随机模型，就是由陀螺和加速度计数据计算wε、w▽、ηε、η▽的功率谱密度，然后形成谱密度阵q。

3 基于时频特性分析的先验滤波模型精化

以X 轴陀螺为例说明利用时频特性分析提取噪声成分功率谱密度信息的方法。

3.1 角速率白噪声的模型

设角速率白噪声对应于陀螺误差中的wεx(wε的X 轴分量)。对于该误差项，只需研究其随机模型，即确定wεx的功率谱密度qωx共有三种方法。

1)功率谱密度法

白噪声的谱密度为常数，对应于双对数功率谱密度(logSω(f)－logf)斜率为0 的部分。先计算静态陀螺数据的功率谱密度Sω(f)，再对双对数谱密度(logSω(f)－logf)拟合出斜率为0 的直线，与纵轴交点为logqωx，qωx即为角速率白噪声的谱密度。若计算的是单边谱密度［9］而非双边谱密度Sω(f)，双对数单边谱拟合直线与纵轴交点为log2qωx。

图1 角速率白噪声谱密度拟合Fig.1 Fitting the PSD of angular rate white noise

2)Allan 方差法

角度随机游走系数与功率谱密度的转换关系为［10，11］:

3)方差法

功率谱密度对频率的积分即为方差，对于某一纯粹的白噪声过程x(t)，Sx(f)为常数，则有:

其中，Δf 为频率带宽。功率谱密度为

如果Δf 未知，可取其上限Fs，Fs 为采样频率。则有qwx的计算公式为

式中，ΔT 为采样间隔。

式(9)是纯粹白噪声过程的功率谱密度求取方法。对于许多陀螺数据，我们无法求取其中白噪声的方差，只能求取陀螺数据整体的方差，而这一方差还包括其他有色噪声成分的影响，因此，该方法在白噪声起主导作用的陀螺上比较有效。

3.2 高斯马尔可夫过程的模型

常用一阶高斯马尔可夫(G-M)过程对陀螺的相关噪声成分进行建模，设一阶G-M 模型表示的X 轴陀螺误差成分为εcgx，则有:

其中，Tc为相关时间，ηεx即为该G-M 误差成分的驱动白噪声，作为组合系统状态噪声的一部分，其谱密度为

其中，σ2为εcgx的方差。

将εcgx引入到组合导航的状态方程中，需要求得其方差σ2及相关时间Tc。虽然通过Allan 方差可以得到一些G-M 过程的信息，但是Allan 方差对G-M 过程的表达不够精确［12］，更理想的分析方法是自相关函数法。

利用自相关函数经过式(12)可以求取σ2及Tc，即可求得驱动白噪声ηεx的谱密度qηx。相关时间Tc=300 s，方差σ2=1.0 的一阶G-M 过程的自相关函数如图2 所示。

图2 一阶G-M 过程的自相关函数Fig.2 Autocorrelation function of 1st order G-M process

对于高阶G-M 过程，其模型比一阶G-M 过程略微复杂，建模后状态参数的个数更多，基本分析方法与一阶G-M 过程类似。

许多文献使用AR 模型描述惯性元件误差，实际上，任意阶G-M 过程都可以用一定阶数的AR 过程表示。这里侧重于用连续的随机过程表达系统的物理特征。

3.3 角速率随机游走的模型

随机游走可以看成是相关时间很长的一阶G-M过程。当Tc→∞时，一阶G-M 过程变成，即随机游走。设角速率随机游走表示的X 轴陀螺误差成分为εcrx，则:

式中，ξηx为该过程的驱动白噪声，谱密度为qξx。易知，随机游走的方差随时间而增长，而其本质上为非平稳过程，因此无法对其进行自相关分析或功率谱分析，需要使用其它方法求取ξηx的功率谱密度qξx。

3.3.1 Allan 方差法

角速率随机游走的功率谱密度为［10］:

K 为速率随机游走系数。根据线性系统传递函数的理论可由驱动白噪声ξηx的谱密度qξx求得随机游走εcrx功率谱密度为［13］:

显然有:

对Allan 方差进行最小二乘拟合或分段最小二乘拟合，可以求取K 的值，如果谱密度的单位取deg2/h4/Hz，可将式(16)用下式代替:

3.3.2 方差法

Maybeck 给出了随机游走过程的驱动白噪声谱密度所满足的关系式［14］，对ξηx而言即:

变形可得:

εcrx(t+ΔT)－εcrx(t)/ΔT 实际上近似为随机游走过程的微分——白噪声ξηx。显然，式(19)本质上与式(9)相同。同理，该方法在角速率随机游走为主要噪声源的陀螺上比较有效。

3.4 应用中的模型与方法选择

若用εcbx表示X 轴陀螺随机常值分量，则有:

在实际应用中，如果对各轴陀螺和加速度计建模，并不一定将上述四种模型都引入到滤波的状态方程中。仍以X 轴陀螺为例，如果同时引入wεx、εcgx、εcrx、εcbx，除白噪声wεx外，其他三种属于非白噪声εc的模型都将增加状态向量的维数，加重滤波解算的负担。一般情况下，白噪声wεx要引入到状态模型，它不会增加滤波负担。εcbx反映一些非随机信息，一般要引入，但在需要保证滤波效率的情况下也可不引入。εcgx和εcrx的作用都是描述噪声的时相关性，实用中只选其一。有文献认为随机游走过程并不是描述角速率误差的现实模型，因为随机游走的方差无限增大，而角速率输出中一般会含有指数相关的一阶马尔可夫过程和宽带白噪声［15］。

在模型精化的各种方法中，Allan 方差法和自相关函数法都需要对长时间的静态惯性数据进行分析才能达到满意的精度，方差法在噪声成分比较单一的条件下才满足其理论条件，否则计算的方差不能反映该类噪声的特点。可见，方法的选择与数据的长度与特性有关，应根据实际的分析结果进行判断取舍。

4 计算与分析

采用的GPS/INS 实验数据共约30 分钟，其中初始静止状态约5 分钟。IMU 采样频率100 Hz，GPS 采样频率1 Hz。IMU 的XYZ 三轴分别指向下、左、前方向，标称陀螺漂移1 deg/h，标称加速度偏置10－4g。对GPS/INS 数据进行松组合导航。采用双差C/A 码解算的位置和Doppller 解算的速度作为GPS 输出值与INS 输出值构成观测量，以载波相位差分获得的位置作为参考解。

组合导航滤波状态变量中，加速度计误差中可忽略时相关误差，这是由于该分量相对较小，同时也为了使滤波器的维数尽量低些［8］，因而这里加速度计误差用随机常值和白噪声进行建模。陀螺误差也包括了随机常值和白噪声，但还使用高斯马尔可夫模型描述时相关误差。整个系统的状态噪声为u=［wεηεw▽］T，其谱密度矩阵为q，有

滤波初始参数及观测值精度设置如下:初始姿态失准角100.0 s、100.0 s、500.0 s，陀螺随机常值和高斯马尔可夫过程初始标准差1.0 deg/h，加速度计随机常值的初始标准差10－4g;GPS 水平方向位置标准差1.0 m，速度标准差0.05 m/s，垂直方向位置标准差2.0 m，速度标准差0.10 m/s。

对于状态随机模型相关的滤波参数q，用时频特性分析方法进行计算。

由于条件限制，只有5 分钟的静止数据。因为数据太短，Allan 方差和自相关函数的估计精度不高，因而这里主要用功率谱密度法提取白噪声谱密度qωx和q▽x，对于高斯马尔可夫噪声，这里的相关时间的估计精度也受数据长度限制，只能根据元件规格取经验值Tc=600 s，再由噪声方差通过式(12)计算高斯马尔可夫驱动白噪声的谱密度qηx。分析计算的结果如表1 所示，将其代入式(24)即可得到频域谱密度矩阵q。

表1 功率谱密度计算结果Tab.1 Calculating results of PSD

对GPS/INS 组合数据分别进行INS 解算、双差C/A 码GPS 解算和GPS/INS 组合解算。图3 所示为INS 纬度及参考纬度，虚线为参考值，点线为INS值;图4 所示为GPS/INS 的纬度误差和GPS 的纬度误差，“+”为GPS 结果，粗实线为GPS/INS 结果，RMS 及最大误差(MAX)如表2 所示;图5 所示为283560.0－283566.0 s 内的GPS 及GPS/INS 纬度值，“+”为GPS 结果，实线为GPS/INS 结果。

表2 GPS 和GPS/INS 的RMS 和MAXTab.2 RMS and MAX of GPS and GPS/INS

图3 INS 纬度与参考纬度Fig.3 INS latitude and the reference latitude

图4 GPS 与GPS/INS 的纬度误差Fig.4 Latitude error of GPS and GPS/INS

由上述结果可以看出:

1)受到惯性元件误差、初始误差等因素的影响，单独INS 的误差随时间积累，需要借助其他导航手段如GNSS 对其误差进行抑制。但是当GNSS 失锁时，只能依赖于INS，此时可通过对INS 误差进行消噪处理来削弱误差。

2)使用时频分析方法对实测数据分析计算得到的滤波参数能够合理地平衡INS 和GPS 这两大系统的关系，可以在不需要过多经验信息和复杂调试的条件下精化滤波先验模型。当先验模型较为准确时，GPS/INS 对GPS 具有类似于拟合的效果，可以平滑掉部分GPS 误差。

3)组合导航系统在保证精度的前提下具有更高的结果输出频率，在本例中，输出频率为100 Hz，结果的精度优于单独的GPS 系统或INS 系统。

4)在少数历元内，组合导航系统的误差大于GPS 的误差，这是因为运动中存在一些状态突变，导致了与先验滤波模型之间的偏差，需要用自适应滤波理论加以处理。

图5 GPS 与GPS/INS 纬度结果对比Fig.5 Comparison between the latitude results of GPS and GPS/INS

5 结语

GNSS/INS 组合导航的状态噪声由陀螺仪和加速度计的误差所构成，建立先验滤波模型的前提是对陀螺和加速度计中的误差建模，尤其是以功率谱密度为代表的频域随机模型的计算。

由于连续系统的频域随机模型和离散系统的时域随机模型可以互相转换，也可以由数据分析计算出几种常见误差源的时域随机模型，但是其过程更为复杂，而且频域的描述更贴近于系统的本质，因此这里主要使用频域模型。

利用时频分析精化先验滤波模型的方法，既减少了确定滤波模型时的主观性和经验性，也避免了实际滤波应用前的繁琐调试。它基于实际的数据进行分析，该方法不仅可以用于组合导航，对其他的Kalman 滤波相关应用领域都具有一定参考价值。

1 Jazwinski A H.Stochastic processes and filtering theory［M］.New York:Mathematics in Science and Engineering，Academic Press，1970.

2 Mohamed A H and Schwarz K P.Adaptive Kalman Filtering for INS/GPS［J］.Journal of Geodesy，1999，73(4):193－203.

3 夏启军，孙优贤，周春晖.渐消卡尔曼滤波器的最佳自适应算法及其应用［J］.自动化学报，1990，16(3)，210－216.(Xia Qijun，Sun Youxian and Zhou Chunhui.An optimal adaptive algorithm for fading Kalman filter and its application［J］.Acta Automatica Sinaca，1990，16(3)，210－216)

4 欧吉坤，柴艳菊，袁运斌.自适应选权滤波［A］.大地测量与地球动力学进展［C］.武汉:湖北科学技术出版社，2004，816－823.(Ou Jikun，Chai Yanju and Yuan Yunbin.Adaptive flitering by selecting the parameter weight factor［A］.Progress in geodesy and geodynamics［C］.Wuhan:Hubei Science and Technology Press，2004，816－823 )

5 Yang Yuanxi，He Haibo and Xu Guochang.Adaptively robust filtering for kinematic geodetic positioning［J］.Journal of Geodesy，2001，75(2/3):109－116.

6 杨元喜，何海波，徐天河.论动态自适应滤波［J］.测绘学报，2001，30(4):293－298.(Yang Yuanxi，He Haibo and Xu Tianhe.Adaptive robust filtering for kinematic GPS positioning［J］.Acta Geodaetica et Cartographica Sinica，30(4):293－298)

7 Gelb A.Applied optimal estimation［M］.USA:The Analytic Science Corporation，1974.

8 秦永元，张洪钺，汪叔华.卡尔曼滤波与组合导航原理［M］.西 安:西 北 工 业 大 学 出 版 社，1998.(Qin Yongyuan，Zhang Hongyue and Wang Shuhua.Kalman filter and the principle of integrated navigation［M］.Xian:Northwest Industrial Press，1998)

9 何正嘉，等.现代信号处理及工程应用［M］.西安:西安交通大学出版社，2007.(He Zhengjia，et al，Modern signal processing and its application［M］.Xi’an:Xi’an Jiaotong University Press，2007)

10 毛奔，林玉荣.惯性器件测试与建模［M］.哈尔滨:哈尔滨工程大学出版社，2008.(Mao Ben and Lin Yurong.Inertial sensor testing and modeling［M］.Harbin:Harbin Engineering University Press，2008)

11 IEEE Std 952-1997.IEEE standard specification format guide and test procedure for single-axis interferometric fiber optic gyros［S］.New York:IEEE Standard Board，1997.

12 Flenniken W S.Modeling inertial measurement units and analyzing the effect of their errors in navigation applications［D］.Graduate Faculty of Auburn University，2005.

13 熊凯，雷拥军，曾海波.基于Allan 方差法的光纤陀螺建模与仿真［J］.空间控制技术与应用，2010，36(3)，7－13.(Xiong Kai，Lei Yongjun and Zeng Haibo.Modeling and simulation of fiber optic gyros based on Allan variance method［J］.Aerospace Control and Application，2010，36(3):7－13)

14 Maybeck P S.Stichastic models，estimation and control，volume 1［M］.New York:Academic Press，1979.

15 Gebre-Egziabher D.Design and performance analysis of a low-cost aided dead reckoning navigator［D］.USA:Stanford University，2004.