一种基于惯性匹配的船体姿态基准传递方法

袁二凯，杨功流，于 沛，唐 改

(1．北京航空航天大学 仪器科学与光电工程学院，北京100191;2．中国人民解放军93716 部队，天津301700)

0 引 言

现代舰船都配备有雷达、导弹发射架、飞机等攻防设备，为保证它们的正常运行，需要舰船为其提供一定精度的姿态、位置等运动参数。舰船上的中心航姿系统是高精度惯导系统或平台罗经，它可以测量并提供这些参数［1］。但是由于船体变形角的存在，设备安装点和中心航姿系统之间存在角偏差，直接影响到设备的使用精度。消除船体变形的影响，获取船体局部高精度姿态信息成为近年来国内外研究的一个热点问题。目前部分舰船采用光学方法测量主惯导与设备之间的船体变形，在中心主惯导高精度姿态基准上补偿变形量以获取战位点处的高精度姿态信息，但光学测量系统成本高昂，结构复杂，应用并不广泛。在各用户设备处安装高精度航姿系统成本很高，而且俄罗斯的早期实践证明这一方法并不完善［2－4］。

本文提出一种基于惯性姿态匹配法实时估计船体变形以实现姿态基准传递的方法。在对光学设备测得的船体变形角数据进行频域分析的基础上，对船体变形角高精度建模。通过在用户设备旁边安装低精度的捷联惯导系统，利用中心主惯导与捷联惯导2 个惯性测量单元(IMU)输出的姿态信息作为观测量，建立卡尔曼滤波方程，实现船体变形角的实时最优估计，在主惯导输出姿态信息的基础上补偿变形量，将高精度姿态基准传递给用户设备。

1 姿态匹配算法

1.1 坐标定义

高精度姿态基准传递方法中各系统的安装位置及坐标定义如图1所示。IMU1 为船体的中心主惯导，是高精度惯导系统或平台罗经，全船的姿态信息均以IMU1 的姿态信息为基准，oy 轴沿着船体的纵轴，oz 轴与甲板面相垂直并指向上方，ox 轴与oy、oz 轴构成右手直角坐标系，即为右前上坐标系。IMU2 安装在用户设备附近，与IMU1 近似对准，为较低精度捷联惯导，坐标系为ox′y′z′。姿态基准传递方案利用IMU1 和IMU2 的姿态输出，估计二者之间的实时船体变形角，将IMU1 的高精度姿态信息传递到IMU2 所在战位点。定义ox 方向的船体变形为纵挠角，oy 方向的船体变形为横挠角，oz 方向的船体变形为首挠角。

图1 系统的安装位置及坐标定义Fig.1 The position of the system and the coordinate definition

对算法中各坐标系定义如下:

定义IMU1 确立的载体坐标系为b1，IMU2 确立的载体坐标系为b2。IMU2 为低精度惯导，由于陀螺漂移等因素的影响，t 时刻IMU2 的计算坐标系会慢慢偏离b2，定义t 时刻IMU2 的计算坐标系b′2。

1.2 姿态匹配算法原理

根据坐标系的定义和坐标变换原理可得:

式(1)中左右两边状态转移矩阵对应的欧拉角相等。记t 时刻2 个IMU 之间的船体变形角为φ对应的姿态欧拉角为θε，则对应的欧拉角可近似为φ+θε。

记

可得

由式(5)可见，2 个IMU 输出姿态信息，实时变形角与IMU2 姿态误差角之间存在线性关系，可以此为基础设计卡尔曼滤波器实现对船体变形角的实时最优估计，从而实现高精度姿态基准的传递。

2 船体变形角的高精度建模

研究表明［1］，船体变形的频率集中在0.3 Hz 以内，由变化较快的动态变形θ 和变化缓慢的静态变形Φ 两部分组成。记t 时刻2 个IMU 之间的船体变形角为φ，可得:

为了实现对船体变形角的准确建模，本文选取一段光学设备采集的船体变形实测数据(理论精度达到1″，可近似视为船体的实际变形)进行了频域分析。

纵挠角的频谱如图2所示。图中频率低于0.3 Hz 的2 个峰值体现的是船体变形角的频率特性，高于0.4 Hz 的峰值反映的是由于设备震动产生的噪声，可以在不影响变形角成分的情况下将其与高频白噪声一起滤除。而且从频率上看，静态变形和动态变形分布界限明显:静态变形频率集中在0.05 Hz 以内，动态变形集中在0.1 ～0.3 Hz，因此在频域内将2 种变形分离并建模可行。滤除高频噪声和分离2 类变形都可以通过小波分析法实现［7］。

图2 纵挠角的频谱图Fig.2 The spectrum curve of the pitching deformation

其中，σ 为静态变形角速率ωΦi的标准差，作为白噪声的驱动系数，可以在分离静态变形角后求得。

动态变形θ 可视为二阶Markov 过程，其对应的变形滤波器方程可表示为:

其中，μ 为动态变形不规则系数 (相关时间的倒数);λ 为动态变形的支配频率;D 为动态变形的方差;b2=μ2+λ2;w(t)为均值为0、方差为1的高斯白噪声。动态变形的支配频率λ 即为图2中0.1 ～0.3 Hz 范围内峰值对应的频率。相关时间和方差可以在分离动态变形后通过数学方法解算得到。

通过在频域内分离静态变形和动态变形，式(7)和式(8)可以实现对船体变形角的精确建模，并作为状态量在卡尔曼滤波器中实现最优估计。

3 卡尔曼滤波器模型的建立

3.1 状态方程的建立

IMU2 的姿态误差角可表示为:

卡尔曼滤波器状态向量取如下18 维形式:

其中:Φ 为静态变形角;ωΦ为静态变形角速度;θ为动态变形角;˙θ 为动态变形角速度;θε为IMU2 陀螺漂移造成的姿态失准角;ε 为IMU2 的陀螺漂移。

系统状态方程的矩阵形式如下:

其中F 为18 ×18 维的稀疏矩阵，由式(8) ～式(10)可得:

G 为系统噪声系数矩阵，是18 ×18 维的对角阵，其中非零元素为:

3.2 量测方程的建立

卡尔曼滤波量测方程的矩阵形式为:

Z 的定义如式(4)，由式(6)、式(7)和式(10)得:

4 仿真验证与结果分析

4.1 仿真条件

1)船体绕x，y，z 三轴以正弦规律做摇摆运动，绕x，y，z 三轴的摇摆周期分别为8 s，7 s，6 s;摇摆幅度分别为4°，5°，3°;初始水平姿态角为0;航向角为30°;仿真时间50 h。

2)IMU1 作为姿态基准，姿态输出不含陀螺漂移等误差信息;IMU2 的姿态输出是在船体姿态角的基础上，加入陀螺漂移和船体变形角。IMU2 的x，y，z 三轴陀螺漂移分别为0.005°/h，0.005°/h，0.01°/h;船体变形角将实测船体变形角速度加到IMU2 陀螺的角速度输出中参与姿态解算。

4.2 仿真结果分析

利用小波分析法滤除高频噪声分离静态和动态变形结果如图3所示。计算可得动态变形模型参数为:μ=［0.1，0.0625，0.0743］，1/s，λ=［0.13，0.1，0.13］，Hz，D=［9.6 ×10－7，2.25 ×10－5，2.2 ×10－7］，deg;静态变形模型参数为:σ=［3.7×10－6，2.7 ×10－6，7.3 ×10－7］，deg。

图3 小波分析法分离静态和动态变形结果Fig.3 Separation result of static and dynamic deformation

卡尔曼滤波估计结果如图4所示，算法实现了对船体静态变形和动态变形的分离估计，IMU2 的x，y，z 三个轴向陀螺漂移分别收敛于0.005°/h，0.005°/h，0.01°/h，与仿真设定值一致。总体变形角的估计结果与理论值对比曲线和估计误差曲线如图5所示，船体变形角的估计结果与理论值一致，x，y，z 三个轴向的变形角估计误差均方差分别为2.3512″，4.565″，3.1242″。姿态匹配算法有效抑制了捷联惯导陀螺漂移的影响，实现了对船体变形角的高精度实时估计。捷联惯导所在局部位置的理想姿态为船体中心姿态和局部变形量之和，实现了船体变形角的高精度估计，在中心主惯导高精度测量结果的基础上补偿变形量就能得到设备所在位置的姿态信息，所以局部位置姿态信息的精度取决于中心主惯导的姿态测量精度和船体变形角的估计精度。

捷联惯导的姿态测量误差和姿态基准传递法得到的姿态误差对比如图6所示。捷联惯导由于陀螺漂移的影响，姿态测量误差发散较快，姿态匹配算法可以实现船体变形角的实时高精度估计，当中心主惯导的姿态测量结果为理想值时，由于变形角估计均方误差在5″以内，所以基准传递精度可以达到5″以内。

图4 卡尔曼滤波估计结果Fig.4 Estimation result of the Kalman filter

图5 变形估计结果与理论值对比曲线及误差曲线Fig.5 Comparison of simulation result and the initial measurement data and the estimation error curve

图6 姿态测量误差和姿态传递误差曲线Fig.6 Curves of attitude measurement error and the attitude reference transfer error

5 结 语

本文提出一种基于惯性姿态匹配法实时估计船体变形以实现姿态基准传递的方法。姿态匹配法利用中心主惯导和用户设备处捷联惯导的姿态输出构建卡尔曼滤波方程，在对船体变形高精度建模的基础上，可以实现船体变形角的实时高精度估计，将主惯导高精度的姿态基准传递给用户设备。仿真验证表明，算法能够有效抑制低精度捷联惯导陀螺漂移因素对船体变形估计结果的影响，中心主惯导的姿态测量结果为理想值时，姿态基准传递精度在5″以内，说明船体变形角建模方法和姿态匹配算法的可行性，为船体局部高精度姿态信息获取提供了理论参考。

［1］MOCHALOV A V．A system for measuring deformation of large-sized objects［C］．RTO/ NATO，France，1999，15:1－15．

［2］MOCHALOV A V，KAZANTSEV．Use of the ring laser units for measurement of the moving object deformation［C］．Proceedings of SPIE，2002(4680):85－92．

［3］朱昀炤，汪顺亭，缪玲娟，等．船体变形测量技术综述［J］．船舶工程，2007，29(6):58－61．

ZHU Yun-zhao，WANG Shun-ting，MIAO Ling-juan，et al．Review of measuring technique for ship deformation［J］．Ship Engineering，2007，29(6):58－61．

［4］TITTERTON D H，WESTON J L．Strapdown inertial navigation technology(Second edition)［M］．The Institution of Electrical Engineers，2004．

［5］郑佳兴，秦石乔，王省书．基于姿态匹配的船体形变测量方法［J］．中国惯性技术学报，2010，18(2):175－180．

ZHENG Jia-xing，QIN Shi-qiao，WANG Xing-shu．Ship hull angular deformation measurement taking slow-varying quasi-static component into account［J］．Journal of Chinese Inertial Technology，2010，18(2):175－180．

［6］SCHNEIDER A M，Kalman filter formulations for transfer alignment of strapdown inertial units［J］．Journal of the Institute of Navigation，1983，30(1):72－89．

［7］WEN M．Spectrum denoise by wavelet analyze and its realize in Matlab［J］．Modern Electronics Technique，2003，26(24):47－49．