CN115265599B - 双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法 - Google Patents

Abstract

本发明公开了一种双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法，步骤为：S1、基于IMU坐标系和地理坐标系，构建旋转调制后的平均地球物理场相关陀螺零偏模型；S2、构建地球物理场相关加速度计零偏模型；S3、构建地球物理场相关陀螺和加速度零偏导致的双轴旋转惯导长航时导航位置误差模型；S4、标定地球物理场相关陀螺零偏并进行误差补偿；该快速标定方法不仅能够精确标定出地球物理场相关陀螺零偏，提高双轴旋转惯导的长航时导航精度，且将标定时间缩短至在2～6h之间，使标定效率提高十数倍；同时与现有标定方法相比，只需要三个采样点的数据；即兼具补偿精度高、标定试验时间短、效率高、对数据采集要求低的多重优势，实用性佳。

Description

双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法

技术领域

本发明涉及双轴旋转惯导误差补偿技术领域，特别涉及双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法。

背景技术

双轴旋转惯导由惯性测量单元(Initial Measurement Unit，简称为IMU)和转位机构组成，IMU安装于转位机构上，通过IMU的对称转位将惯性器件的常值误差调制成均值为零的周期变化量，实现在现有的惯性器件精度水平上大幅度提高双轴旋转惯导长航时导航精度的目的。因此，双轴旋转惯导被广泛地应用在以远洋船舶为代表的需要高精度长航时导航信息的应用场合中。

目前，双轴旋转惯导中惯性器件的常值误差补偿技术较为成熟，例如，授权发明专利CN109459063B提供了一种双轴旋转惯导的十六位置误差调制方法，发明专利CN113465631A公开了一种基于对角线旋转的双轴旋转惯导误差调制方法。在通过旋转调制消除常值误差后，无法通过旋转调制消除的非常值误差成为了影响双轴旋转惯导长航时导航精度最主要的误差源。其中，地球物理场相关陀螺零偏是双轴旋转惯导非常值误差中最主要的误差源。地球物理场相关陀螺零偏包括地磁场相关陀螺零偏和重力场相关陀螺零偏，地磁场相关陀螺零偏和重力场相关陀螺零偏分别是地球磁场和重力场导致的陀螺零偏误差。由于地球磁场和重力场的方向均固定，因此，双轴旋转惯导中每个陀螺的地球物理场相关陀螺零偏随着转位机构的旋转而不断变化，是非常值误差，严重影响了双轴旋转惯导的长航时导航精度。

为了解决上述问题，申请号为CN 2022108702755的发明专利提供了一种双轴旋转惯导地球物理场相关陀螺零偏的补偿方法，该方法能准确的标定出地球物理场相关陀螺零偏，然而，该方法在制定标定方案时忽略地球振荡，为了保证标定精度，该方法的限定条件为采集的总导航时间需要大于48小时，且需要采集标定实验过程全程的数据，存在着标定时间长、对数据采集要求高的问题。

发明内容

本发明的目的是提供一种解决目前现有的双轴旋转惯导地球物理场相关陀螺零偏的补偿方法存在着的标定时间长、对数据采集要求高的问题的双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法。

为此，本发明技术方案如下：

一种双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法，步骤如下：

S1、构建IMU坐标系，即m系；构建地理坐标系，即n系；基于两个坐标系，构建地球物理场相关陀螺零偏模型，进而构建旋转调制后的平均地球物理场相关陀螺零偏模型；其中，

地球物理场相关陀螺零偏模型为：

式中，为在m系中地球物理场相关陀螺零偏向量，/>为在m系中地磁场相关陀螺零偏向量，/>为在m系中重力场相关陀螺零偏向量；

旋转调制后的平均地球物理场相关陀螺零偏模型为一个三维向量，其三个分量均为非零常数，其表达式为：

式中，为旋转调制后的平均地球物理场相关陀螺零偏在X_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Z_n轴的分量；/>和/>均为非零常数；

S2、构建地球物理场相关加速度计零偏模型，其表达式为：

式中，旋转调制后的平均地球物理场相关加速度计零偏向量，/>为旋转调制后的平均地球物理场相关加速度计零偏在X_n轴的分量，/>为旋转调制后的平均地球物理场相关加速度计零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关加速度计零偏在Z_n轴的分量；/>和/>均为非零常数；

S3、构建地球物理场相关陀螺和加速度零偏导致的双轴旋转惯导长航时导航位置误差模型，其表达式为：

式中，δL(t)、δλ(t)分别为导航时间为t时双轴旋转惯导的纬度误差和经度误差；L(t)、λ(t)分别为导航时间为t时双轴旋转惯导输出的纬度和经度；ω_ie为地球自转角速率；δα₀为初始横摇角误差，δβ₀为初始纵摇角误差，δγ₀为初始航向角误差；为旋转调制后的平均地球物理场相关陀螺零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Z_n轴的分量；

S4、标定地球物理场相关陀螺零偏并进行误差补偿；

S401、构建标定方程中的观测向量的分块向量为：

z_j＝[δL(t_j)δλ(t_j)]^T，

式中，j为标定实验数据采集的第j个采样点，z_j为第j个采样点对应的标定方程中的观测向量的第j个分块向量，t_j为第j个采样点对应的导航时间，δL(t_j)和δλ(t_j)分别为第j个采样点对应的纬度误差和经度误差；

构建标定系数向量x为：

构建观测矩阵的分块矩阵为：

式中，j为标定实验数据采集的第j个采样点，H_j为第j个采样点对应的标定方程中的观测矩阵中的第j个分块矩阵，t_j为第j个采样点对应的导航时间，S_j＝sinω_iet_j，C_j＝cosω_iet_j；

进而，构建标定方程为：

z＝Hx+n，

式中，z为标定方程的观测向量，z＝[z₁ z₂ z₃]^T，z₁、z₂和z₃分别为第1个采样点、第2个采样点和第3个采样点获得的观测向量，即：z_j＝[δL(t_j)δλ(t_j)]^T；n为观测误差；H为标定方程的观测矩阵，H＝[H₁ H₂ H₃]^T，H₁、H₂和H₃分别为第1个采样点、第2个采样点和第3个采样点获得的观测矩阵；

S402、设计标定实验方案为：1)将双轴旋转惯导静止放置在纬度L₀、经度λ₀的实验室内并将纬度L₀和经度λ₀装订到双轴旋转惯导中；2)待双轴旋转惯导对准后，在导航过程中实施现有的十六位置旋转调制方法用以消除所有的常值误差；3)对准完成后，进入导航状态的第一个时刻记为第一个采样点时刻t₁，获得第一个采样点时刻导航时间t₁和t₁时刻下双轴旋转惯导输出的经度导航结果λ(t₁)和纬度导航结果L(t₁)，而后通过公式：δλ(t₁)＝λ(t₁)-λ₀和δL(t₁)＝L(t₁)-L₀，计算得到t₁时刻对应的双轴旋转惯导的经度误差为δλ(t₁)和纬度误差δL(t₁)；4)基于采样周期T_S，将t₂＝t₁+T_S时刻作为第二个采样点时刻为t₂，获得第二个采样点时刻导航时间t₂和t₂时刻下双轴旋转惯导输出的经度导航结果λ(t₂)和纬度导航结果L(t₂)，而后通过公式：δλ(t₂)＝λ(t₂)-λ(t₁)和δL(t₂)＝L(t₂)-L(t₁)，计算得到t₂时刻对应的双轴旋转惯导的经度误差为δλ(t₂)和纬度误差δL(t₂)；5)将t₃＝t₂+T_S时刻作为第二个采样点时刻为t₃，获得第二个采样点时刻导航时间t₃和t₃时刻下双轴旋转惯导输出的经度导航结果λ(t₃)和纬度导航结果L(t₃)，而后通过公式：δλ(t₃)＝λ(t₃)-λ(t₂)和δL(t₃)＝L(t₃)-L(t₂)，计算得到t₃时刻对应的双轴旋转惯导的经度误差为δλ(t₃)和纬度误差δL(t₃)；6)将由步骤3)至步骤5)获得的三个采样点的数据代入至标定方程：z＝Hx+n中，利用最小二乘法求解标定方程，以得到标定系数向量；

S403、标定结果处理：将步骤S402获得的标定系数向量，以及 代入至由步骤S302得到的时域下的双轴旋转惯导长航时导航位置误差模型中，获得实时的纬度误差δL(t)和经度误差δλ(t)，进而将双轴旋转惯导输出的实时纬度L(t)和经度λ(t)分别减去δL(t)和δλ(t)，即得到补偿后的纬度和经度。

进一步地，在步骤S1中，

在m系中地磁场相关陀螺零偏向量的表达式为：

式中，M为陀螺地磁场感应系数矩阵，其表达式为：M_xx为X_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_xy为Y_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_xz为Z_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_yx为X_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_yy为Y_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_yz为Z_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_zx为X_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量，M_zy为Y_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量，M_zz为Z_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量；Hⁿ为在n系中地磁场强度向量，其表达式为：Hⁿ＝[0 H_G 0]^T，H_G为地磁场强度大小；/>为IMU在转位机构旋转下的姿态矩阵；

在m系中重力场相关陀螺零偏向量的表达式为：

式中，G为陀螺重力场感应系数矩阵，G_xx为X_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_xy为Y_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_xz为Z_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_yx为X_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_yy为Y_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_yz为Z_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_zx为X_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量，G_zy为Y_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量，G_zz为Z_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量；gⁿ为在n系中重力加速度向量；gⁿ为在n系中重力加速度向量，gⁿ＝[0 0 -g]^T，g为重力加速度大小。

进一步地，步骤S3的具体实施过程如下：

S301、将步骤S1和步骤S2构建的旋转调制后的平均地球物理场相关陀螺零偏模型：和加速度计零偏模型：/>代入惯导误差方程，得到地球物理场相关的陀螺零偏和加速度计零偏导致的双轴旋转惯导导航误差模型：

式中，ψ为姿态误差向量，ψ＝[δα δβ δγ]^T，其中，δα为横摇角误差，δβ为纵摇角误差，δγ为航向角误差；vⁿ为速度向量， 为东向速度，/>为北向速度，为天向速度；δvⁿ为速度误差向量，/> 为东向速度误差，/>为北向速度误差，/>为天向速度误差；L、λ和h分别为纬度、经度和高度；δL、δλ和δh分别为纬度误差、经度误差和高度误差；/>为n系相对于惯性系的角速率向量；/>为n系相对于惯性系的角速率向量误差；/>为m系中陀螺组件测量的角速率向量；/>为m系中陀螺组件测量的角速率向量误差；fⁿ为在n系中加速度组件测量的加速度计向量；f^m为在m系中加速度组件测量的加速度计向量；δf^m为在m系中加速度组件测量的加速度计向量误差；/>为在n系中地球自转角速率；/>为在n系中n系相对于地球的角速率向量；R_N和R_E分别为当地地球的子午圈半径和卯酉圈半径；

S302、将步骤S301构建的双轴旋转惯导导航误差模型进行拉普拉斯变换，同时忽略傅科振荡和舒勒振荡这两个对长航时位置精度影响很小的误差项，得到时域下的双轴旋转惯导长航时导航位置误差模型。

进一步地，在步骤S402中，采样周期T_S满足：1h＜T_S＜3h。优选地，在步骤S402中，采样周期T_S为2h。

与现有技术相比，该双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法的有益效果在于：

1)本申请的方法能够精确标定出地球物理场相关陀螺零偏，并通过补偿能较好地提高双轴旋转惯导的长航时导航精度，经过实验验证，本申请的方法在静态条件下将双轴惯导的经度解算精度提高了68.52％，在海上动态条件下将双轴惯导的经度解算精度提高了53.33％；

2)相对于现有的双轴旋转惯导地球物理场相关陀螺零偏的补偿方法对标定试验时间需要大于48小时的要求对比，本申请的方法将标定时间缩短至在2～6h之间，使标定效率提高十数倍；

3)相对于现有的双轴旋转惯导地球物理场相关陀螺零偏的补偿方法在标定试验中需要采集标定试验全程的数据，对数据采集有较高的要求，本申请的方法只需要采集3个采样点的数据，对数据采集要求低；

综上，本申请兼具补偿精度高、标定试验时间短、效率高、对数据采集要求低的多重优势，实用性佳。

附图说明

图1为本发明的双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法的流程示意图；

图2(a)为本发明对比例1、对比例2和实施例1在静态条件下双轴旋转惯导的位置导航精度中的纬度精度的对比示意图；

图2(b)为本发明对比例1、对比例2和实施例1在静态条件下双轴旋转惯导的位置导航精度中的经度精度的对比示意图；

图3(a)为本发明对比例1、对比例2和实施例1在海上动态条件下双轴旋转惯导的位置导航精度中的纬度精度的对比示意图；

图3(b)为本发明对比例1、对比例2和实施例1在海上动态条件下双轴旋转惯导的位置导航精度中的经度精度的对比示意图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

如图1所示，该双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法的具体实施步骤如下：

S1、构建地球物理场相关陀螺零偏模型；

具体地，该步骤S1的实施步骤如下：

S101、构建IMU坐标系(m系)和地理坐标系(n系)；其中，

1)IMU坐标系，即m系，其表达式为o-X_mY_mZ_m系，该坐标系为正交坐标系，其原点o点为双轴旋转惯导中IMU内陀螺与加速度计的测量中心点，其X_m轴与IMU内的X向陀螺敏感轴方向一致，其Y_m轴在IMU内X向陀螺敏感轴和Y向陀螺敏感轴构成的平面中且垂直于X_m轴，其Z_m轴垂直于X_m轴与Y_m轴构成的平面；

2)地理坐标系，即n系，其表达式为o-X_nY_nZ_n系，该坐标系为正交坐标系，其原点o点为m系的原点，其X_n轴指向地理东向，其Y_n轴指向地理北向，其Z_n轴指向地理天向；

S102、构建地球物理场相关陀螺零偏中地磁场相关陀螺零偏模型；

在m系中，构建地磁场相关陀螺零偏模型为：

其中，为在m系中地磁场相关陀螺零偏向量，/>ε_Mx为地磁场相关陀螺零偏向量在X_m轴上的分量，ε_My为地磁场相关陀螺零偏向量在Y_m轴上的分量，ε_Mz为地磁场相关陀螺零偏向量在Z_m轴上的分量；M为陀螺地磁场感应系数矩阵，其中，M_xx为X_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_xy为Y_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_xz为Z_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_yx为X_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_yy为Y_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_yz为Z_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_zx为X_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量，M_zy为Y_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量，M_zz为Z_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量；H^m为在m系中地磁场强度向量，H^m＝[H_xH_y H_z]^T，其中，H_x为X_m轴上的地磁场强度分量，H_y为Y_m轴上的地磁场强度分量，H_z为Z_m轴上的地磁场强度分量；

由于地磁场强度在地理坐标系中的方向为从南指向北，因此，在n系中地磁场强度向量表示为：

Hⁿ＝[0 H_G 0]^T，

式中，Hⁿ为在n系中地磁场强度向量；H_G为地磁场强度大小；

考虑到双轴旋转惯导中转位机构的旋转，在m系中地磁场相关陀螺零偏模型的表达式为：

式中，为IMU在转位机构旋转下的姿态矩阵；

S103、构建地球物理场相关陀螺零偏中重力场相关陀螺零偏模型；

在m系中，构建重力场相关陀螺零偏模型为：

式中，为在m系中重力场相关陀螺零偏向量，/>ε_Gx为重力场相关陀螺零偏向量在X_m轴上的分量，ε_Gy为重力场相关陀螺零偏向量在Y_m轴上的分量，ε_Gz为重力场相关陀螺零偏向量在Z_m轴上的分量；G为陀螺重力场感应系数矩阵，其中，G_xx为X_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_xy为Y_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_xz为Z_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_yx为X_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_yy为Y_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_yz为Z_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_zx为X_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量，G_zy为Y_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量，G_zz为Z_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量；g^m为在m系中重力加速度向量，g^m＝[g_x g_yg_z]^T，g_x为X_m轴上的重力加速度分量，g_y为Y_m轴上的重力加速度分量，g_z为Z_m轴上的重力加速度分量；

由于重力加速度在地理坐标系中的方向为从天指向地，因此，在n系中重力加速度向量表示为：

gⁿ＝[0 0 -g]^T，

其中，gⁿ为在n系中重力加速度向量，g为重力加速度大小；

考虑到双轴旋转惯导中转位机构的旋转，在m系中重力场相关陀螺零偏模型写为：

其中，为IMU在转位机构旋转下的姿态矩阵；

S104、综合地磁场相关陀螺零偏和重力场相关陀螺零偏，得到地球物理场相关陀螺零偏模型；

由于地球物理场相关陀螺零偏由地磁场相关陀螺零偏和重力场相关陀螺零偏组成，地球物理场相关陀螺零偏由步骤S102构建的地磁场相关陀螺零偏和步骤S103构建的重力场相关陀螺零偏相加得到，其表达式为：

式中，为在m系中地球物理场相关陀螺零偏向量，/>ε_Tx为地球物理场相关陀螺零偏向量在X_m轴上的分量，ε_Ty为地球物理场相关陀螺零偏向量在Y_m轴上的分量，ε_Tz为地球物理场相关陀螺零偏向量在Z_m轴上的分量；

S105、构建旋转调制后的平均地球物理场相关陀螺零偏模型；

由旋转调制原理，在旋转调制后的在n系中的平均地球物理场相关陀螺零偏模型由下式计算得到：

式中，为在n系中旋转调制后的平均地球物理场相关陀螺零偏，/>为转位机构旋转角度，k为倍数，由于在双轴旋转惯导的旋转调制中，转位机构按照180°的整数倍进行对称均匀旋转，所以k为正整数，其取值为1,2,3,4……；

由步骤S102～S104中的各式，的计算结果为三个分量均为非零常数的三维向量，/>记为：

式中，为旋转调制后的平均地球物理场相关陀螺零偏在X_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Z_n轴的分量；/>和/>均为非零常数；

S2、构建地球物理场相关加速度计零偏模型；

具体地，该步骤S2的实施步骤如下：

按照与步骤S1相同的思路，构建旋转调制后的平均地球物理场相关加速度计零偏模型，其表达式为：

式中，旋转调制后的平均地球物理场相关加速度计零偏向量，/>为旋转调制后的平均地球物理场相关加速度计零偏在X_n轴的分量，/>为旋转调制后的平均地球物理场相关加速度计零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关加速度计零偏在Z_n轴的分量；/>和/>均为非零常数；

S3、构建地球物理场相关陀螺和加速度零偏导致的双轴旋转惯导长航时导航位置误差模型；

具体地，该步骤S3的实施步骤如下：

S301、将步骤S1和步骤S2构建的旋转调制后的平均地球物理场相关陀螺零偏模型：和加速度计零偏模型：/>代入惯导误差方程，得到地球物理场相关的陀螺零偏和加速度计零偏导致的双轴旋转惯导导航误差模型：

式中，ψ为姿态误差向量，ψ＝[δα δβ δγ]^T，其中，δα为横摇角误差，δβ为纵摇角误差，δγ为航向角误差；vⁿ为速度向量， 为东向速度，/>为北向速度，为天向速度；δvⁿ为速度误差向量，/> 为东向速度误差，/>为北向速度误差，/>为天向速度误差；L、λ和h分别为纬度、经度和高度；δL、δλ和δh分别为纬度误差、经度误差和高度误差；/>为n系相对于惯性系的角速率向量；/>为n系相对于惯性系的角速率向量误差；/>为m系中陀螺组件测量的角速率向量；/>为m系中陀螺组件测量的角速率向量误差；fⁿ为在n系中加速度组件测量的加速度计向量；f^m为在m系中加速度组件测量的加速度计向量；δf^m为在m系中加速度组件测量的加速度计向量误差；/>为在n系中地球自转角速率；/>为在n系中n系相对于地球的角速率向量；R_N和R_E分别为当地地球的子午圈半径和卯酉圈半径；

S302、将步骤S301构建的双轴旋转惯导导航误差模型进行拉普拉斯变换，同时忽略傅科振荡和舒勒振荡这两个对长航时位置精度影响很小的误差项，得到时域下的双轴旋转惯导长航时导航位置误差模型，其表达式为：

式中，δL(t)、δλ(t)分别为导航时间为t时双轴旋转惯导的纬度误差和经度误差；L(t)、λ(t)分别为导航时间为t时双轴旋转惯导输出的纬度和经度；ω_ie为地球自转角速率；δα₀为初始横摇角误差，δβ₀为初始纵摇角误差，δγ₀为初始航向角误差；为旋转调制后的平均地球物理场相关陀螺零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Z_n轴的分量；从时域下的双轴旋转惯导长航时导航位置误差模型可以看出，地球物理场相关陀螺零偏/>中的/>和/>对双轴旋转惯导长航时导航位置误差有累加影响，而/>对双轴旋转惯导长航时导航误差无累加影响，同时，地球物理场相关加速度计零偏对双轴旋转惯导长航时导航误差无累加影响；

S4、标定地球物理场相关陀螺零偏并进行误差补偿；

具体地，步骤S3的实施步骤如下：

S401、构建标定方程：

S4011、基于步骤S302构建的时域下的双轴旋转惯导长航时导航位置误差模型，构建标定方程中的观测向量的分块向量，其表达式为：

z_j＝[δL(t_j) δλ(t_j)]^T，

式中，j为标定实验数据采集的第j个采样点，z_j为第j个采样点对应的标定方程中的观测向量的第j个分块向量，t_j为第j个采样点对应的导航时间，δL(t_j)和δλ(t_j)分别为第j个采样点对应的纬度误差和经度误差；

S4012、构建标定系数向量：

由于地球物理场相关陀螺零偏中的/>和/>对双轴旋转惯导长航时导航位置误差有累加影响，而/>对双轴旋转惯导长航时导航误差无累加影响，因此，构建的标定系数向量需要包含/>和/>以标定出对双轴旋转惯导长航时导航位置误差有累加影响的误差源；同时，地球物理场相关加速度计零偏对双轴旋转惯导长航时导航误差无累加影响并考虑标定方程求解的可行性，

进而，基于步骤S302构建的时域下的双轴旋转惯导长航时导航位置误差模型，以及地球物理场相关陀螺零偏中的/>和/>构建标定方程中的标定系数向量x，其表达式为：

S4013、基于S4012构建的标定方程中的标定系数向量和步骤S302构建的时域下的双轴旋转惯导长航时导航位置误差模型，构建标定方程中的观测矩阵的分块矩阵，其表达式为：

式中，j为标定实验数据采集的第j个采样点，H_j为第j个采样点对应的标定方程中的观测矩阵中的第j个分块矩阵，t_j为第j个采样点对应的导航时间，S_j＝sinω_iet_j，C_j＝cosω_iet_j；

S4014、根据步骤S4011～S4013构建的标定方程中的观测向量的分块向量、标定系数向量和观测矩阵的分块矩阵；并考虑标定系数向量为六维向量，且每个采样点能获得一个二维的观测向量的分块向量以及两行的观测矩阵的分块矩阵，因此，要完成标定系数向量的求解，需要三个采样点的数据构成六维的观测向量和六行的观测矩阵；因此，标定方程构建为：

z＝Hx+n，

式中，z为标定方程的观测向量，z＝[z₁ z₂ z₃]^T，z₁、z₂和z₃分别为第1个采样点、第2个采样点和第3个采样点获得的观测向量，即：z_j＝[δL(t_j) δλ(t_j)]^T；n为观测误差；H为标定方程的观测矩阵，H＝[H₁ H₂ H₃]^T，H₁、H₂和H₃分别为第1个采样点、第2个采样点和第3个采样点获得的观测矩阵，即：

S402、设计标定实验方案：

1)双轴旋转惯导静止放置在纬度L₀、经度λ₀的实验室内，开机，在双轴旋转惯导中装订实验室准确的位置基准信息，即装订的纬度为L₀，经度为λ₀；

2)双轴旋转惯导在对准后进入导航状态，旋转方案采用现有的十六位置误差调制方法，以旋转调制消除所有的常值误差；其中，上述十六位置误差调制方法可参见已授权专利CN109459063B的一种双轴旋转惯导的十六位置误差调制方法；

3)对准完成后，进入导航状态的第一个时刻记为第一个采样点时刻t₁，获得第一个采样点时刻导航时间t₁和t₁时刻下双轴旋转惯导输出的经度导航结果λ(t₁)和纬度导航结果L(t₁)，而后通过公式：δλ(t₁)＝λ(t₁)-λ₀和δL(t₁)＝L(t₁)-L₀，计算得到t₁时刻对应的双轴旋转惯导的经度误差为δλ(t₁)和纬度误差δL(t₁)；

4)基于设定的采样周期T_S，将t₂＝t₁+T_S时刻作为第二个采样点时刻为t₂，获得第二个采样点时刻导航时间t₂和t₂时刻下双轴旋转惯导输出的经度导航结果λ(t₂)和纬度导航结果L(t₂)，而后通过公式：δλ(t₂)＝λ(t₂)-λ(t₁)和δL(t₂)＝L(t₂)-L(t₁)，计算得到t₂时刻对应的双轴旋转惯导的经度误差为δλ(t₂)和纬度误差δL(t₂)；

5)同样地，基于设定的采样周期T_S，将t₃＝t₂+T_S时刻作为第二个采样点时刻为t₃，获得第二个采样点时刻导航时间t₃和t₃时刻下双轴旋转惯导输出的经度导航结果λ(t₃)和纬度导航结果L(t₃)，而后通过公式：δλ(t₃)＝λ(t₃)-λ(t₂)和δL(t₃)＝L(t₃)-L(t₂)，计算得到t₃时刻对应的双轴旋转惯导的经度误差为δλ(t₃)和纬度误差δL(t₃)；

6)将由步骤3)至步骤5)获得的三个采样点的数据代入至标定方程：z＝Hx+n中，利用最小二乘法求解标定方程，以得到标定系数向量；

其中，为了保证标定方程求解的非奇异性，采样周期T_S应满足小于地球振荡周期的1/8，因此，基于地球振荡周期为24h，采样周期T_S应小于3h的条件；同时，考虑到导航误差随时间累加，因而为了提高标定方程的可观测性，采样周期T_S应大于1h；进而，采样周期T_S应满足：1h＜T_S＜3h，相应地，总标定时间应满足：2h＜总标定时间＜6h；优选地，采样周期T_S为2h，总标定时间为4h；

S403、将步骤S402获得的标定系数向量代入步骤S202获得的时域下的双轴旋转惯导长航时导航位置误差模型，并进行误差补偿；

具体地，将步骤S402获得的标定系数向量，以及代入至由步骤S302得到的时域下的双轴旋转惯导长航时导航位置误差模型中，获得实时的纬度误差δL(t)和经度误差δλ(t)，进而将双轴旋转惯导输出的实时纬度L(t)和经度λ(t)分别减去δL(t)和δλ(t)，即得到补偿后的纬度和经度。

为了验证本发明提出方法的有效性和正确性，分别构建对比例1、对比例2和实施例1，并分别在静态实验和海上动态实验中获得三种情况下的导航位置，以实现基于现有技术(即对比例1、对比例2)，对本申请实现的导航位置精度提升效果进行对比。

对比例1

选用某双轴旋转惯导作为对比例1的实验对象，选用的双轴旋转惯导中IMU由三个精度为0.008°/h的激光陀螺仪和三个精度为10ug的加速度计组成，转位机构的定位精度为5″，对选用的双轴旋转惯导只采用发明专利CN109459063B提供的一种双轴旋转惯导的十六位置误差调制方法，以旋转调制消除所有的常值误差，而不对双轴旋转惯导地球物理场相关陀螺零偏进行补偿。

对比例2

采用与对比例1相同的双轴旋转惯导作为对比例2的实验对象，选用的双轴旋转惯导首先采用授权发明专利CN109459063B提供的一种双轴旋转惯导的十六位置误差调制方法，以旋转调制消除所有的常值误差；然后采用背景技术中提及的申请号为CN2022108702755的发明专利提供的一种双轴旋转惯导地球物理场相关陀螺零偏的补偿方法对双轴旋转惯导地球物理场相关陀螺零偏进行标定和补偿；在双轴旋转惯导地球物理场相关陀螺零偏的标定中，标定实验时间为125h。

实施例1

采用与对比例1相同的双轴旋转惯导作为实施例1的实验对象，选用的双轴旋转惯导首先采用授权发明专利CN109459063B提供的一种双轴旋转惯导的十六位置误差调制方法，以旋转调制消除所有的常值误差；然后采用本发明提出的方法对双轴旋转惯导地球物理场相关陀螺零偏进行标定和补偿；在双轴旋转惯导地球物理场相关陀螺零偏的标定中，采样周期T_S设置为2h，标定实验时间为4h。

(一)静态试验：

利用GPS的信息作为位置基准，进而分别在对比例1、对比例2和实施例1中获得双轴旋转惯导在实验室静态条件下的位置精度。

如图2(a)所示为对比例1、对比例2和实施例1分别在实验室静态条件下的位置精度中的纬度精度对比示意图；从图中可知，纬度误差由于不随导航时间增大而线性增加，最大纬度误差没有明显变化；附图2(b)为对比例1、对比例2和实施例1在实验室静态条件下的位置精度中的经度精度对比示意图，从图中可知，利用本发明方法将双轴惯导的静态情况下最大经度误差从2.7海里减小到0.85海里，精度提高了68.52％，与对比例2实现的精度提高程度相当，但从标定时间考量，实施例1的标定时间相比对比例2相比节约了96.8％，且在标定实验过程中只需要采集三个采样点的数据，可见实施例1相对于对比例2更具高效性，操作上也更具简便性。

(二)海上动态试验：

利用GPS的信息作为位置基准，进而分别在对比例1、对比例2和实施例1中获得安装于海上船舶的双轴旋转惯导在海上动态条件下的位置精度。

如图3(a)所示为对比例1、对比例2和实施例1分别在海上动态条件下的位置精度中的纬度精度对比示意图；从图中可知，纬度误差由于不随导航时间增大而线性增加，最大纬度误差没有明显变化；如图3(b)所示为对比例1、对比例2和实施例1分别在海上动态条件下的位置精度中的经度精度对比示意图；从图中可知，利用本发明方法将双轴惯导的海上动态条件下最大经度误差从3海里减小到1.4海里，精度提高了53.33％，与对比例2的精度提高程度相当，但从标定时间考量，实施例1的标定时间相比对比例2相比节约了96.8％，且在标定实验过程中只需要采集三个采样点的数据，可见实施例1相对于对比例2更具高效性，操作上也更具简便性，同时实施例1的方法与对比例2相同，均可适用于静止环境中和动态环境中。

综上所述，对比例1、对比例2和实施例1在静态条件和海上动态条件的位置精度对比表明了本发明能够精确标定出地球物理场相关陀螺零偏，通过补偿能较好地提高双轴旋转惯导的长航时导航精度，补偿精度与现有的双轴旋转惯导地球物理场相关陀螺零偏的补偿方法相当，但与现有的双轴旋转惯导地球物理场相关陀螺零偏的补偿方法(对比例2)对标定试验时间需要大于48小时的要求对比，本申请的方法最大优势与提升是标定时间，以实施例1为例的标定时间仅需要4h，效率相对现有技术提高12倍；同时，现有的双轴旋转惯导地球物理场相关陀螺零偏的补偿方法在标定试验中需要采集标定试验全程的数据，对数据采集有较高的要求，而本申请的方法只需要采集3个采样点的数据，对数据采集要求低；因此，与现有技术相比，本申请的方法兼具补偿精度高、标定试验时间短、效率高、以及对数据采集要求低的多重优势，实用性好。

本发明未详细公开的部分属于本领域的公知技术。尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化时显而易见的，一切利用本发明构思的发明创造均为保护之列。

Claims (5)

1.一种双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法，其特征在于，步骤如下：

S1、构建IMU坐标系，即m系；构建地理坐标系，即n系；基于两个坐标系，构建地球物理场相关陀螺零偏模型，进而构建旋转调制后的平均地球物理场相关陀螺零偏模型；其中，

地球物理场相关陀螺零偏模型为：

式中，为在m系中地球物理场相关陀螺零偏向量，/>为在m系中地磁场相关陀螺零偏向量，/>为在m系中重力场相关陀螺零偏向量；

旋转调制后的平均地球物理场相关陀螺零偏模型为一个三维向量，其三个分量均为非零常数，其表达式为：

式中，为旋转调制后的平均地球物理场相关陀螺零偏在X_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Z_n轴的分量；/>和/>均为非零常数；

S2、构建地球物理场相关加速度计零偏模型，其表达式为：

式中，旋转调制后的平均地球物理场相关加速度计零偏向量，/>为旋转调制后的平均地球物理场相关加速度计零偏在X_n轴的分量，/>为旋转调制后的平均地球物理场相关加速度计零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关加速度计零偏在Z_n轴的分量；/>和/>均为非零常数；

S3、构建地球物理场相关陀螺和加速度零偏导致的双轴旋转惯导长航时导航位置误差模型，其表达式为：

式中，δL(t)、δλ(t)分别为导航时间为t时双轴旋转惯导的纬度误差和经度误差；L(t)、λ(t)分别为导航时间为t时双轴旋转惯导输出的纬度和经度；ω_ie为地球自转角速率；δα₀为初始横摇角误差，δβ₀为初始纵摇角误差，δγ₀为初始航向角误差；为旋转调制后的平均地球物理场相关陀螺零偏在Y_n轴的分量，/>为旋转调制后的平均地球物理场相关陀螺零偏在Z_n轴的分量；

S4、标定地球物理场相关陀螺零偏并进行误差补偿；

S401、构建标定方程中的观测向量的分块向量为：

z_j＝[δL(t_j) δλ(t_j)]^T，

式中，j为标定实验数据采集的第j个采样点，z_j为第j个采样点对应的标定方程中的观测向量的第j个分块向量，t_j为第j个采样点对应的导航时间，δL(t_j)和δλ(t_j)分别为第j个采样点对应的纬度误差和经度误差；

构建标定系数向量x为：

构建观测矩阵的分块矩阵为：

式中，j为标定实验数据采集的第j个采样点，H_j为第j个采样点对应的标定方程中的观测矩阵中的第j个分块矩阵，t_j为第j个采样点对应的导航时间，S_j＝sinω_iet_j，C_j＝cosω_iet_j；

进而，构建标定方程为：

z＝Hx+n，

式中，z为标定方程的观测向量，z＝[z₁ z₂ z₃]^T，z₁、z₂和z₃分别为第1个采样点、第2个采样点和第3个采样点获得的观测向量，即：z_j＝[δL(t_j) δλ(t_j)]^T；n为观测误差；H为标定方程的观测矩阵，H＝[H₁ H₂ H₃]^T，H₁、H₂和H₃分别为第1个采样点、第2个采样点和第3个采样点获得的观测矩阵；

S402、设计标定实验方案为：1)将双轴旋转惯导静止放置在纬度L₀、经度λ₀的实验室内并将纬度L₀和经度λ₀装订到双轴旋转惯导中；2)待双轴旋转惯导对准后，在导航过程中实施现有的十六位置旋转调制方法用以消除所有的常值误差；3)对准完成后，进入导航状态的第一个时刻记为第一个采样点时刻t₁，获得第一个采样点时刻导航时间t₁和t₁时刻下双轴旋转惯导输出的经度导航结果λ(t₁)和纬度导航结果L(t₁)，而后通过公式：δλ(t₁)＝λ(t₁)-λ₀和δL(t₁)＝L(t₁)-L₀，计算得到t₁时刻对应的双轴旋转惯导的经度误差为δλ(t₁)和纬度误差δL(t₁)；4)基于采样周期T_S，将t₂＝t₁+T_S时刻作为第二个采样点时刻为t₂，获得第二个采样点时刻导航时间t₂和t₂时刻下双轴旋转惯导输出的经度导航结果λ(t₂)和纬度导航结果L(t₂)，而后通过公式：δλ(t₂)＝λ(t₂)-λ(t₁)和δL(t₂)＝L(t₂)-L(t₁)，计算得到t₂时刻对应的双轴旋转惯导的经度误差为δλ(t₂)和纬度误差δL(t₂)；5)将t₃＝t₂+T_S时刻作为第二个采样点时刻为t₃，获得第二个采样点时刻导航时间t₃和t₃时刻下双轴旋转惯导输出的经度导航结果λ(t₃)和纬度导航结果L(t₃)，而后通过公式：δλ(t₃)＝λ(t₃)-λ(t₂)和δL(t₃)＝L(t₃)-L(t₂)，计算得到t₃时刻对应的双轴旋转惯导的经度误差为δλ(t₃)和纬度误差δL(t₃)；6)将由步骤3)至步骤5)获得的三个采样点的数据代入至标定方程：z＝Hx+n中，利用最小二乘法求解标定方程，以得到标定系数向量；

S403、标定结果处理：将步骤S402获得的标定系数向量，以及 代入至由步骤S302得到的时域下的双轴旋转惯导长航时导航位置误差模型中，获得实时的纬度误差δL(t)和经度误差δλ(t)，进而将双轴旋转惯导输出的实时纬度L(t)和经度λ(t)分别减去δL(t)和δλ(t)，即得到补偿后的纬度和经度。

2.根据权利要求1所述的双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法，其特征在于，在步骤S1中，

在m系中地磁场相关陀螺零偏向量的表达式为：

式中，M为陀螺地磁场感应系数矩阵，其表达式为：M_xx为X_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_xy为Y_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_xz为Z_m轴上的地磁场强度分量导致的X_m轴上的陀螺零偏分量，M_yx为X_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_yy为Y_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_yz为Z_m轴上的地磁场强度分量导致的Y_m轴上的陀螺零偏分量，M_zx为X_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量，M_zy为Y_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量，M_zz为Z_m轴上的地磁场强度分量导致的Z_m轴上的陀螺零偏分量；Hⁿ为在n系中地磁场强度向量，其表达式为：Hⁿ＝[0 H_G 0]^T，H_G为地磁场强度大小；/>为IMU在转位机构旋转下的姿态矩阵；

在m系中重力场相关陀螺零偏向量的表达式为：

式中，G为陀螺重力场感应系数矩阵，G_xx为X_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_xy为Y_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_xz为Z_m轴上的重力加速度分量导致的X_m轴上的陀螺零偏分量，G_yx为X_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_yy为Y_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_yz为Z_m轴上的重力加速度分量导致的Y_m轴上的陀螺零偏分量，G_zx为X_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量，G_zy为Y_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量，G_zz为Z_m轴上的重力加速度分量导致的Z_m轴上的陀螺零偏分量；gⁿ为在n系中重力加速度向量；gⁿ为在n系中重力加速度向量，gⁿ＝[0 0 -g]^T，g为重力加速度大小。

3.根据权利要求1所述的双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法，其特征在于，步骤S3的具体实施过程如下：

S301、将步骤S1和步骤S2构建的旋转调制后的平均地球物理场相关陀螺零偏模型：和加速度计零偏模型：/>代入惯导误差方程，得到地球物理场相关的陀螺零偏和加速度计零偏导致的双轴旋转惯导导航误差模型：

式中，ψ为姿态误差向量，ψ＝[δα δβ δγ]^T，其中，δα为横摇角误差，δβ为纵摇角误差，δγ为航向角误差；vⁿ为速度向量， 为东向速度，/>为北向速度，/>为天向速度；δvⁿ为速度误差向量，/> 为东向速度误差，/>为北向速度误差，/>为天向速度误差；L、λ和h分别为纬度、经度和高度；δL、δλ和δh分别为纬度误差、经度误差和高度误差；/>为n系相对于惯性系的角速率向量；/>为n系相对于惯性系的角速率向量误差；/>为m系中陀螺组件测量的角速率向量；/>为m系中陀螺组件测量的角速率向量误差；fⁿ为在n系中加速度组件测量的加速度计向量；f^m为在m系中加速度组件测量的加速度计向量；δf^m为在m系中加速度组件测量的加速度计向量误差；/>为在n系中地球自转角速率；/>为在n系中n系相对于地球的角速率向量；R_N和R_E分别为当地地球的子午圈半径和卯酉圈半径；

S302、将步骤S301构建的双轴旋转惯导导航误差模型进行拉普拉斯变换，同时忽略傅科振荡和舒勒振荡这两个对长航时位置精度影响很小的误差项，得到时域下的双轴旋转惯导长航时导航位置误差模型。

4.根据权利要求1所述的双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法，其特征在于，在步骤S402中，采样周期T_S满足：1h＜T_S＜3h。

5.根据权利要求4所述的双轴旋转惯导地球物理场相关陀螺零偏的快速标定方法，其特征在于，在步骤S402中，采样周期T_S为2h。