0 引　言

姿态基准传递技术通过准确测量船体局部点位与基准惯导所在点位之间的变形角，使船载系统（如雷达系统、武器发射系统及导航系统等）工作在统一的坐标基准之下，从而保证各系统的协同工作。在诸多姿态基准传递方法中，基于惯性量匹配的方法凭借其测量精度高、受环境影响小等优势，已成为广受研究的主流方法^{[1 – 4]}。惯性量匹配法的原理是利用主惯导和待测点位捷联惯导的惯性量输出建立匹配关系、借助Kalman滤波实现对船体变形角的估计，保证其估计精度的前提之一是确保2组惯性量之间严格的数据同步。为解决数据同步问题，文献[5]与文献[6]分别研究了时间延迟对于角速度匹配算法与惯性姿态匹配算法的精度影响，并设计了时延在线估计补偿算法，但均未考虑从采集方式上解决时延问题；文献[7]提出通过时钟同步解决数据同步采集问题的方案，但缺少对变形角测量结果的试验验证。

本文在角速度匹配算法的基础上，首先通过仿真对时延误差的影响进行研究。为消除时延对姿态基准传递过程的影响，一方面提出基于晶振校正的时间同步采集方案，以实现从采集机制上控制时延，同时设计了引入时延的扩维Kalman滤波器，进一步从算法角度完成对时延误差的补偿，达到高精度测量船体变形角的目的。通过仿真试验，验证了所提基于时间同步的船体姿态基准传递方案的准确性与可靠性，对于该领域的研究具有一定的参考价值。

1 角速度匹配算法 1.1 角速度匹配原理

基于角速度匹配的船体姿态基准传递方法利用2套光纤陀螺捷连惯导INS1和INS2，它们的安装位置及坐标定义如图1所示。其中，INS1作为船体主惯导，提供全船设备的姿态基准信息，定义其坐标系为 ${b_1}$ 。INS2安装于待测位置，定义其坐标系为 ${b_2}$ ，INS1与INS2大致对准，但由于船体变形角的存在， ${b_1}$ 、 ${b_2}$ 之间存在一个小角度偏差。

如图2所示，定义总变形角 ${{\varphi }} = {\left[ {{{{\varphi }}_x},{{{\varphi }}_y},{{{\varphi }}_z}} \right]^{\rm{T}}}$ 。由文献[8]研究工作可知，总变形角中同时包含静态成分与动态成分，其中静态角由温度变化、日光照射等慢变因素导致，动态角由海浪冲击、随机风力等快变因素导致，即有

${{\varphi }} = {{\varPhi }} + {{\theta}}(t)\text{。}$

(1)

此时INS1敏感 ${b_1}$ 系相对 $i$ 系的角速度在 ${b_1}$ 系下的投影，表示为 ${{\omega}}_{i{b_1}}^{{b_1}}$ ；INS2敏感 ${b_2}$ 系相对于 $i$ 系的角速度在 ${b_2}$ 系下的投影，表示为 ${{\omega}}_{i{b_2}}^{{b_2}}$ 。可推导出如下关系：

${{\omega}}_{i{b_2}}^{{b_2}} = {{C}}{{_{{b_2}}^{{b_1}}}^{\rm{T}}}{{\omega}}_{i{b_1}}^{{b_1}} + {{\omega}}_{{b_1}{b_2}}^{{b_2}} = {{C}}{{_{{b_2}}^{{b_1}}}^{\rm{T}}}{{\omega}}_{i{b_1}}^{{b_1}} + {\dot{{\theta}}} \text{，}$

(2)

在变形角为小角度的前提下，由小角度近似得到转换矩阵 ${{C}}_{{b_2}}^{{b_1}}$ 形式如下：

$C_{{b_2}}^{{b_1}} \approx \left[ {\begin{array}{*{20}{c}} 1&{ - {\varphi _z}}&{{\varphi _y}} \\ {{\varphi _z}}&1&{ - {\varphi _x}} \\ { - {\varphi _y}}&{{\varphi _x}}&1 \end{array}} \right] = {{I}} + {\rm{(}}{{\varphi }} \times {\rm{)}}\text{，}$

(3)

将式（3）代入式（2）中，并与 ${{\omega}}_{i{b_1}}^{{b_1}}$ 作差，进一步推导得到角速度匹配方程：

$ \Delta {{{\omega}} = {{\omega}} }_{i{b_1}}^{{b_1}} - {{\omega}}_{i{b_2}}^{{b_2}} = ({{\varphi }} \times ){{\omega}}_{i{b_1}}^{{b_1}} - {\dot{{\theta}}} = ( - {{\omega}}_{i{b_1}}^{{b_1}} \times ){{\varphi }} - {\dot{{\theta}}} \text{，} $

(4)

其中，

$ - {{\omega}}_{i{b_1}}^{{b_1}} \times = \left[ {\begin{array}{*{20}{c}} 0&{{{\omega}}_{i{b_1}z}^{{b_1}}}&{ - {{\omega}}_{i{b_1}y}^{{b_1}}} \\ { - {{\omega}}_{i{b_1}z}^{{b_1}}}&0&{{{\omega}}_{i{b_1}x}^{{b_1}}} \\ {{{\omega}}_{i{b_1}y}^{{b_1}}}&{ - {{\omega}}_{i{b_1}x}^{{b_1}}}&0 \end{array}} \right]\text{。} $

由角速度匹配方程式（4）可知，角速度差值与总变形角、动态变形角速度之间存在线性关系，在此基础上合理设计Kalman滤波器，可实现对变形角的最优估计，实现船体姿态基准的传递。

1.2 卡尔曼滤波方程

在建立Kalman滤波方程之前，首先对静态变形角、动态变形角进行建模。在本文研究中，将静态变形角建模为短时间内的常值模型，即有 ${\dot{{\varPhi}}} {\rm{ = 0}}$ ；对于动态变形角，建模为白噪声驱动的二阶Markov过程，其滤波器方程可表示为：

$ {\ddot \theta _i} + 2{\mu _i}{\dot \theta _i} + b_i^2{\theta _i} = 2{b_i}{\sigma _i}\sqrt {{\mu _i}} \left[ {w(t)} \right],{\rm{ }}i = x,y,z\text{。} $

(5)

其中， $b_i^2 = \mu _i^2 + \lambda _i^2$ ， ${\mu _i}$ 为衰减因子， ${\lambda _i}$ 为支配频率； ${\sigma _i}$ 为幅度标准差； $w(t)$ 为均值为0、方差为1的高斯白噪。

对角速度匹配方程做进一步整理，在式（4）基础上，考虑光纤陀螺漂移影响，得到：

$\Delta {{\omega}} = - {{\omega}}_{i{b_1}}^{{b_1}} \times {{\varphi }} - {\dot{{\theta}}} + {{\tilde{{\varepsilon}}}_1} - {{\tilde{{\varepsilon}}}_2}\text{，}$

(6)

式中： ${{\tilde{{\varepsilon}}}_1}$ 和 ${{\tilde{{\varepsilon}}}_2}$ 分别为2套INS的陀螺漂移，包含常值漂移与随机漂移。记INS1的常值漂移与随机漂移分别为 ${{{\varepsilon}}_{1c}}$ 与 ${{{\varepsilon}}_{1r}}$ ，INS2的常值漂移与随机漂移分别为 ${{{\varepsilon}}_{2c}}$ 与 ${{{\varepsilon}}_{2r}}$ 。认为常值漂移为恒定值，并提出随机漂移误差模型：

${\dot \varepsilon _{ri}} + {\kappa _i}{\varepsilon _{ri}} = {\sigma '_i}\sqrt {2{\kappa _i}} w(t),{\rm{ }}i = x,y,z\text{。}$

(7)

式中： ${\kappa _i}$ ， ${\sigma '_i}$ 分别为不规则系数与表征漂移不规则程度的均方差。将状态向量设为如下21维形式：

$ {{x}} = {\left[ {\begin{array}{*{20}{c}} {{\Phi }}&{{\theta}}&{{\dot{{\theta}}} }&{{{{\varepsilon}}_{1c}}}&{{{{\varepsilon}}_{2c}}}&{{{{\varepsilon}}_{1r}}}&{{{{\varepsilon}}_{2r}}} \end{array}} \right]^{\rm{T}}}\text{，} $

(8)

进一步写出系统状态方程：

${{\dot{{x}}} = {{Fx}} + {{Gw}}}\text{，}$

(9)

其中：

$\begin{array}{l} {{F}} = \left[ {\begin{array}{*{20}{c}} {{{{F}}_1}}&{{{{0}}_{9 \times 12}}} \\ {{{{0}}_{12 \times 9}}}&{{{{F}}_2}} \end{array}} \right]\text{，}\\ {{{F}}_1} = \left[ {\begin{array}{*{20}{c}} {{{{0}}_{3 \times 3}}}&{{{{0}}_{3 \times 3}}}&{{{{0}}_{3 \times 3}}} \\ {{{{0}}_{3 \times 3}}}&{{{{0}}_{3 \times 3}}}&{{{{I}}_{3 \times 3}}} \\ {{{{0}}_{3 \times 3}}}&{{A_1}}&{{A_2}} \end{array}} \right]\text{，} \end{array}$

${A_1} = diag\{ \left[ {\begin{array}{*{20}{c}} { - b_x^2,}&{ - b_y^2,}&{ - b_z^2} \end{array}} \right]\}\text{，} $

${A_2} = diag\{ \left[ {\begin{array}{*{20}{c}} { - 2{\mu _x},}&{ - 2{\mu _y},}&{ - 2{\mu _z}} \end{array}} \right]\} \text{，}$

$ {F_2} = diag\{ \left[ {{0_{{\rm{1}} \times {\rm{6}}}}, - {\kappa _{1x}}, - {\kappa _{1y}}, - {\kappa _{1z}}, - {\kappa _{2x}}, - {\kappa _{2y}}, - {\kappa _{2z}}} \right]\}\text{，} $

$ {{G}} = {\left[ {\begin{array}{*{20}{c}} {{{{0}}_{3 \times 3}}}&{{{{0}}_{3 \times 3}}}&{{{{G}}_1}}&{{{{0}}_{3 \times 3}}}&{{{{0}}_{3 \times 3}}}&{{{{G}}_2}}&{{{{G}}_3}} \end{array}} \right]^{\rm{T}}}\text{，} $

$ {{{G}}_1} = diag\{ \left[ {\begin{array}{*{20}{c}} {2{b_x}{\sigma _x}\sqrt {{\mu _x}} ,}&{2{b_y}{\sigma _y}\sqrt {{\mu _y}} ,}&{2{b_z}{\sigma _z}\sqrt {{\mu _z}} } \end{array}} \right]\}\text{，} $

$ {{{G}}_2} = diag\{ \left[ {\begin{array}{*{20}{c}} {{{\sigma '}_{1x}}\sqrt {2{\kappa _{1x}}} ,}&{{{\sigma '}_{1y}}\sqrt {2{\kappa _{1y}}} ,}&{{{\sigma '}_{1z}}\sqrt {2{\kappa _{1z}}} } \end{array}} \right]\}\text{，} $

${{{G}}_3} = diag\{ \left[ {\begin{array}{*{20}{c}} {{{\sigma '}_{2x}}\sqrt {2{\kappa _{2x}}} }&{{{\sigma '}_{2y}}\sqrt {2{\kappa _{2y}}} }&{{{\sigma '}_{2z}}\sqrt {2{\kappa _{2z}}} } \end{array}} \right]\} \text{。}$

系统量测方程为：

${{z}} = {{Hx}} + {{v}}\text{，}$

(10)

其中，量测量 ${{z}} = \Delta {{\omega}} = {{\omega}}_{ib}^b - {{\omega}}_{ib}^s$ ， ${{v}}$ 为量测噪声，量测矩阵 ${{H}}$ 为：

${{H}} = \left[ {\begin{array}{*{20}{c}} {{{{H}}_1}}&{{{{H}}_1}}&{ - {{{I}}_{3 \times 3}}}&{{{{I}}_{3 \times 3}}}&{ - {{{I}}_{3 \times 3}}}&{{{{I}}_{3 \times 3}}}&{ - {{{I}}_{3 \times 3}}} \end{array}} \right]\text{。}$

其中， $ {{{H}}_1} = - {{\omega}}_{ib}^b \times \varphi$ 。

2 时间同步方法设计 2.1 时延产生的原因

对于实际测量条件来说，由于2套系统信息采集与信息输出时间的差异，将会导致二者数据存在延迟。如果时间延迟量过大，将给测量值引入较大的误差，影响算法的收敛速度，并最终导致船体变形角估计值的准确性严重降低，无法完成姿态基准的精确传递。

时延产生的原因可以归结为两点：一是由于2套惯导系统内部晶振不同步，导致数据采集过程存在延时，该延时量大概在40～120 ms之间；二是由于当2套惯导系统相隔较大距离，且传输线路的材质密度不同时，会使得信息传输速度存在差异从而导致信息传输过程存在延时。

2.2 基于晶振校正的同步采集

经2.1节分析，导致时延的很大方面原因在于2套系统的晶振不同步，因此直接采集的2组数据将存在时延偏差。为从采集机制设计上解决时延问题，本论文设计基于晶振校正的同步采集方案，即将GPS的1PPS作为外部基准源，利用1PPS信号驯服2套惯导系统的晶振^{[9 – 10]}。

由于GPS提供的1PPS信号具有很高的精度，一般接收机输出的1PPS信号的上升沿与UTC的统计误差精度为1 μm，并且不存在累计误差，所以可用1PPS信号对2套系统陀螺晶振进行驯服校正，从而保证2套系统的时钟同步。另外，惯导系统在接收到1PPS信号的同时也接收UTC时间信息，可以为惯导数据添加时间标记，即“打时戳”，便于在数据同步采集单元内部进行时间取齐，同步采集方案如图3所示。

基于上述同步采集原理，理论上能将时延量控制在10 ms以内。但仅靠该方法无法彻底解决时延影响，还需从算法角度对残余时延误差进行补偿。

2.3 引入时延量的扩维滤波器

1.2节对于Kalman滤波器的设计中未考虑2组陀螺的数据延迟，在此条件下，将角速度匹配方程重写为：

$ \Delta {{\omega}} = {{{\omega}}_{1{t_1}}} - {{{\omega}}_{2{t_1}}} = - {{\omega}}_{i{b_1}}^{{b_1}} \times {{\varphi }} - {\dot{{\theta}}} + {{\tilde{{\varepsilon}}}_1} - {{\tilde{{\varepsilon}}}_2}\text{，} $

(11)

实际情况下，假设INS1在 ${t_2}$ 时刻的输出与INS2在 ${t_1}$ 时刻的输出相对应，即时延为 $\Delta t = {t_1} - {t_2}$ 。在此情况下，对于坐标系 $o{x_2}{y_2}{z_2}$ 而言，其在 ${t_1}$ 与 ${t_2}$ 时刻必定存在角速度差，即有

${{{\omega}}_{2{t_2}}} - {{{\omega}}_{2{t_1}}} = {{\dot{{\omega}}}_2} \cdot \Delta\text{，} $

(12)

将式（12）代入式（11）中，得到考虑时延之后的角速度匹配方程：

$ \begin{split} \Delta {{\omega}} =& {{{\omega}}_{1{t_2}}} - {{{\omega}}_{2{t_1}}} = {{{\omega}}_{1{t_2}}} - ({{{\omega}}_{2{t_2}}} - {{{\dot{{\omega}}}}_2} \cdot \Delta t) = \\ & - {{\omega}}_{i{b_1}}^{{b_1}} \times {{\varphi }} - {\dot{{\theta}}} + {{{\tilde{{\varepsilon}}}}_1} - {{{\tilde{{\varepsilon}}}}_2} + {{{\dot{{\omega}}}}_2} \cdot \Delta t\text{。} \end{split} $

(13)

将时延 $\Delta t$ 扩充为滤波状态量，即状态量取如下22维形式：

$ {{x}}=\left[ {\begin{array}{*{20}{c}} {{{\Phi}}}&{{{\theta}}}&{{\dot{{\theta}}}}&{{{{\varepsilon}}_{{1c}}}}&{{{{\varepsilon}}_{{2c}}}}&{{{{\varepsilon}}_{{1r}}}}&{{{{\varepsilon}}_{{2r}}}}&{\Delta t}^{\rm{T}} \end{array}} \right]\text{，}$

(14)

新的状态转移阵 ${{F'}}$ 与噪声驱动阵 ${{G'}}$ ，如下式：

${{F'}} = \left[ {\begin{array}{*{20}{c}} {{{{F}}_1}}&{{{{0}}_{9 \times 12}}}&{{{{0}}_{9 \times 1}}} \\ {{{{0}}_{12 \times 9}}}&{{{{F}}_2}}&{{{{0}}_{12 \times 1}}} \\ {{{{0}}_{1 \times 9}}}&{{{{0}}_{1 \times 12}}}&0 \end{array}} \right]\text{，}$

(15)

${{G'}} = {\left[ {\begin{array}{*{20}{c}} {{{{0}}_{3 \times 3}}}&{{{{0}}_{3 \times 3}}}&{{{{G}}_1}}&{{{{0}}_{3 \times 3}}}&{{{{0}}_{3 \times 3}}}&{{{{G}}_2}}&{{{{G}}_3}}&{{{{0}}_{3 \times 1}}} \end{array}} \right]^{\rm{T}}}\text{，}$

(16)

量测阵 ${{H'}}$ 为：

${{H'}} = \left[ {\begin{array}{*{20}{c}} {{{{H}}_1}}\;\;{{{{H}}_1}}\;\;{ - {{{I}}_{3 \times 3}}}\;\;{{{{I}}_{3 \times 3}}}\;\;{ - {{{I}}_{3 \times 3}}}\;\;{{{{I}}_{3 \times 3}}}\;\;{ - {{{I}}_{3 \times 3}}}\;\;{{{{H}}_2}} \end{array}} \right]\text{，}$

(17)

其中，

${{{H}}_2} = \left[ {\begin{array}{*{20}{c}} {{{{\dot{{\omega}}}}_{2x}}} \\ {{{{\dot{{\omega}}}}_{2x}}} \\ {{{{\dot{{\omega}}}}_{2x}}} \end{array}} \right]\text{。}$

3 仿真试验与分析

为验证本文所提基于时间同步算法的姿态基准传递方法的有效性，设置仿真试验。仿真中假设船体在系统时处于摇摆状态，分别绕纵摇轴、横摇轴和方位轴以正弦形式进行三轴模拟运动；2套惯导系统陀螺精度分别为0.01°/h和0.05°/h；假设经数据同步采集之后，残余时延长度分别为10 ms（理想情况）与40 ms（非理想情况），对2种情况船体变形角度进行仿真解算，仿真时间为10 min，得到结果如图4 ~ 图9所示。

整理后的估计误差见表1，分析图表结果可得出结论：引入时延的扩维滤波器能够准确估计出时延量，时延补偿后的估计曲线可以较好地跟随真实值的变化，且10 ms时延情况下的解算情况优于40 ms情况。由表1可知，时延补偿算法能够达到10 ms左右的估计精度，较未补偿情况提升明显。另外，鉴于时延增大时算法的估计精度有所下降，若配合本文所提基于晶振校正的时间同步采集方案将时延量控制在10 ms以内，再结合时延补偿算法进行变形角估计，理论上可得到较为理想的解算精度，从而实现高精度的船体姿态基准传递。

4 结　语

本文研究基于角速度匹配算法的船体高精度姿态基准传递技术，并重点针对其中2套惯导系统数据同步问题开展方案设计，从2个角度入手解决时延问题：一是从数据采集机制上，利用GPS提供的1PPS作为高精度基准时钟源，设计基于晶振校正的同步采集方案；二是从算法设计上，研究引入时延量的扩维Kalman滤波算法，通过2种方法相结合以实现对时延偏差的控制与补偿，从而准确解算出船体变形角，达到船体姿态基准高精度传递的目的。

从仿真试验结果来看，本文所提基于时间同步的姿态基准传递方案，能够将测量误差控制在10′左右，具有一定的工程价值。下一步工作应针对该方案在实船环境下的解算效果做进一步研究。