基于Psi角误差修正模型的极地双惯导协同标定方法

基于psi角误差修正模型的极地双惯导协同标定方法
技术领域
1.本发明属于导航技术领域，涉及惯导系统的外场标定方法，特别涉及基于psi角误差修正模型的极地双惯导协同标定方法，适用于两套及以上带有双轴或三轴转位机构的惯导系统间在极地的协同标定。


背景技术：

2.极地地区自然资源丰富，战略地位重要，具有很高的科学研究和战略价值，北极问题已超出北极区域的范畴，攸关人类生存和发展的共同命运。为确保舰船等载体平台在极地的安全航行，载体精确的定位导航技术是需要解决的关键问题之一。由于经线收敛于地理极点、地磁线收敛于极点附近，此外，极地环境复杂，磁暴、太阳风暴频发，这导致很多低纬度地区常用的导航手段并不适用于极地范围内的导航定位。惯性导航能够为载体连续提供姿态、速度、位置信息，是极地环境中最主要导航方式。然而惯性导航在极区也面临着计算溢出导致误差增大、缺少航向参考等问题，需要在进入极区后切换导航坐标系。
3.由于惯导系统自身特性以及外界环境的影响，惯导系统在一段时间后其误差参数会发生变化，这会严重影响惯导系统的定位精度。此外，在极地环境中，恶劣的环境会导致惯导系统的元器件发生故障的可能性增大。这些情况下需要对惯导系统重新进行标定以保证惯导系统的定位精度。在中低纬度下载体一般具备外界准确参考信息作为观测，利用卡尔曼滤波即可对惯导系统实现系统级标定。然而在极区特殊的环境下，能够接收到的外界参考信息极其有限，外场标定技术的使用会受到限制。这对于需要定期标定或者故障后修复的惯导系统而言其导航精度会受到严重影响。具备极地航行能力的舰船航时长，对系统可靠性的要求高，通常会搭载多套带有转位机构的惯导系统。利用两套惯导系统的冗余信息，以两套惯导系统间的相对速度、位置作为约束观测，构建联合状态卡尔曼滤波器能够使惯导系统的系统性误差得到标定估计。传统的phi角误差模型定义在真实坐标系，然而真实坐标系是未知的，通常将计算坐标系近似为真实坐标系使用，存在近似误差，而psi角误差模型定义在计算坐标系，与位置误差解耦，更适合极区长航时航行的舰船。此外，若在舰船运动过程中进行标定，传统的速度误差方程中的比力项需微分获得，在动态环境下误差较大，也会影响标定的精度。
4.本发明针对目前存在的问题，提出基于psi角误差修正模型的极地双惯导协同标定方法，适用于装备多套带有转位机构惯导系统的载体在极地环境下的标定。本发明以两套惯导系统在横平台坐标系下的相对速度、横计算坐标系下的相对位置作为约束观测，建立了地球椭球模型横地球坐标系下基于psi角误差修正模型的双惯导系统联合状态卡尔曼滤波器对待标定惯导系统的误差进行估计。该方法不受载体运动状态的影响，在静基座、动基座下均能对待标定惯导系统的全误差参数进行在线标定，解决了极地无外界基准信息情况下惯导系统的外场标定难题；以两套惯导系统间的相对误差为观测量，标定精度不受惯导系统绝对误差的影响；通过误差修正模型消除了模型中的比力项，提高了在运动状态下的标定精度；psi角误差模型定义在计算坐标系，与位置误差解耦，更适合极地航行的长航
时舰船惯导系统的标定。


技术实现要素：

5.本发明提出基于psi角误差修正模型的极地双惯导协同标定方法，实现了在极地环境无外界基准信息时，对具备自标定能力的惯导系统的陀螺标度因数、加速度计标度因数、安装误差角的标定，psi角误差模型更适合极地航行的长航时舰船的惯导系统，对速度误差模型进行修正提高了运动状态下的标定精度，本发明的标定精度能够满足导航级惯导系统的需求，具有重要工程实用价值。
6.为解决上述技术问题，本发明提出的解决方案为：
7.基于psi角误差修正模型的极地双惯导协同标定方法，所述方法包括以下步骤：
8.(1)构建基于地球椭球模型下的横地球坐标系和横计算坐标系；
9.以北纬0
°
东经90
°
点为横地球坐标系下的北极点，定义为横北极点，北纬0
°
西经90
°
点为横地球坐标系下的南极点，定义为横南极点，0
°
经线和180
°
经线包围的椭圆面为横赤道面，取横北极、横南极和北极点组成的半个大椭圆为0
°
横经线，所在平面为横本初子午线，将地球坐标系与新定义的横地球坐标系之间的转换关系表示为：
[0010][0011]
式中，e系表示地球坐标系，e
′
系表示横地球坐标系，表示地球坐标系与横地球坐标系转换的方向余弦矩阵；
[0012]
基于横经纬网定义横计算坐标系，横北向指向横北极，所在位置的法线向上为天向，按照右手坐标系定义横东向，将横计算坐标系与计算坐标系之间的转换关系表示为：
[0013][0014]
式中，c系表示计算坐标系，c
′
系表示横计算坐标系，β表示计算坐标系与横计算坐标系之间的旋转角度；
[0015]
确定β与经纬度、横经纬度之间的转换关系：
[0016][0017][0018]
将横计算坐标系与横平台坐标系之间的转换关系表示为：
[0019][0020]
式中，p
′
系表示横平台坐标系，i表示单位矩阵，[ψ
×
]表示横计算坐标系下漂移误差角的反对称矩阵；
[0021]
将载体所在位置的法线与横赤道面的夹角定义为横纬度，与横本初子午面的夹角
定义为横经度，将地球坐标系下定义的经纬度与横经纬度之间的转换关系表示为：
[0022][0023][0024][0025][0026]
式中，l
t
表示横纬度，λ
t
表示横经度，l表示纬度，λ表示经度；
[0027]
(2)构建两套惯导系统的误差模型；
[0028]
定义正常工作的双轴旋转调制惯导系统为惯导1，其体坐标系b1定义为“右-前-上”，待标定的惯导系统为惯导2，其体坐标系b2定义为“右-前-上”；
[0029]
惯导1的标度因数误差及安装误差很小进而忽略，将惯导1的误差模型定义为：
[0030][0031][0032]
其中，
[0033][0034][0035]
式中，表示惯导1的陀螺组件误差，表示惯导1的加速度计组件误差，表示惯导1的x轴陀螺漂移，表示惯导1的y轴陀螺漂移，表示惯导1的z轴陀螺漂移，表示惯导1的x轴加速度计零偏，表示惯导1的y轴加速度计零偏，表示惯导1的z轴加速度计零偏，表示惯导1的陀螺漂移，表示惯导1的加速度计零偏，为惯导1的陀螺噪声，为惯导1加速度计噪声；
[0036]
考虑标度因数误差、安装角误差及零偏误差，将惯导2的误差模型定义为：
[0037][0038][0039]
其中，
[0040][0041][0042]
式中，表示惯导2的陀螺组件误差，表示惯导2的加速度计组件误差，表示惯导2陀螺组件输出的理论角速度矢量，表示惯导2加速度计组件测得的理论比力矢量，表示惯导2的x轴陀螺漂移，表示惯导2的y轴陀螺漂移，表示惯导2的z轴陀螺漂
移，表示惯导2的x轴加速度计零偏，表示惯导2的y轴加速度计零偏，表示惯导2的z轴加速度计零偏，表示惯导2的陀螺漂移，表示惯导2的加速度计零偏，为惯导2的陀螺噪声，为惯导2加速度计噪声；δκg和δμg表示陀螺的标度因数误差矩阵和安装误差矩阵，δκa和δμa表示加速度计的标度因数误差矩阵和安装误差矩阵；
[0043]
确定δκg和δκa：
[0044][0045][0046]
式中，δκ
gx
、δκ
gy
和δκ
gz
分别表示x轴陀螺、y轴陀螺和z轴陀螺的标度因数误差，δκ
ax
、δκ
ay
和δκ
az
分别表示x轴加速度计、y轴加速度计和z轴加速度计的标度因数误差；
[0047]
确定δμg和δμa：
[0048][0049][0050]
式中，δμ
gyx
、δμ
gzx
和δμ
gzy
表示陀螺组件的三个安装误差角，δμ
ayx
、δμ
azx
、δμ
azy
、δμ
axy
、δμ
axz
和δμ
ayz
表示加速度计组件的六个安装误差角；
[0051]
(3)利用两套惯导系统输出的姿态、速度、位置相关信息，联合状态卡尔曼滤波器，具体步骤为：
[0052]
(3.1)确定系统联合误差方程：
[0053][0054][0055][0056][0057][0058][0059]
其中，
[0060]
[0061][0062][0063]
式中，ψ1表示横计算坐标系下惯导1的漂移误差角，表示横平台坐标系下误差修正后惯导1的速度误差矢量，表示横计算坐标系下惯导1的位置误差，表示横计算坐标系下惯导1的地球自转角速度，表示横计算坐标系下惯导1的转移角速度，表示惯导1体坐标系至横平台坐标系的方向余弦矩阵，表示在横计算坐标系下惯导1的重力矢量，ψ2表示横计算坐标系下惯导2的漂移误差角，表示横平台坐标系下误差修正后惯导2的速度误差矢量，表示横计算坐标系下惯导2的位置误差，表示惯导2横计算坐标系下的地球自转角速度，表示横计算坐标系下惯导2的转移角速度，表示惯导2体坐标系至横平台坐标系的方向余弦矩阵，表示在横计算坐标系下惯导2的重力矢量，v
p
′
表示载体在横平台坐标系下的速度，h表示载体所在位置的高度，r
x
为载体所在位置的横东向的曲率半径，ry为载体所在位置的横北向的曲率半径，为载体所在位置的扭曲率，re和rn分别为载体所在位置的卯酉圈半径和子午圈半径；
[0064]
(3.2)确定联合状态方程：
[0065][0066]
其中，
[0067][0068]
[0069][0070][0071][0072][0073][0074][0075][0076]
[0077][0078][0079][0080][0081][0082]
式中，0i×j表示i行j列的零矩阵，和为惯导1输出的横纬度和横经度，和为惯导2输出的横纬度和横经度，分别表示横平台坐标系下载体在横东向、横北向、天向的速度，ω
ie
表示地球自转角速度大小，表示惯导1输出位置处的重力加速度的值，表示惯导2输出位置处的重力加速度的值，c
23
表示矩阵的第二、第三列，c3表示矩阵的第三列，表示惯导2的x轴陀螺输出值，表示惯导2的y轴陀螺输出值，表示惯导2的z轴陀螺输出值，表示惯导2的x轴加速度计输出值，表示惯导2的y轴加速度计输出值，表示惯导2的z轴加速度计输出值；
[0083]
将状态向量x(t)表示为：
[0084][0085]
式中，ψ
e1
、ψ
n1
、ψ
u1
分别表示惯导1在横东向、横北向、天向的漂移误差角，分别表示惯导1在横东向、横北向、天向的漂移误差角，分别表示在横平台坐标系下惯导1在横东向、横北向、天向的速度误差，表示惯导1在横东向的位置误差，表示惯导1在横北向的位置误差，表示惯导1的天向误差，ψ
e2
、ψ
n2
、ψ
u2
表示惯导2在横东向、横北向、天向的漂移误差角，分别表示横平台坐标系下惯导2在横东向、横北向、天向的速度误差，表示惯导2的在横东向
的位置误差，表示惯导2的在横北向的位置误差，表示惯导2在天向的位置误差；
[0086]
将噪声分布矩阵及噪声矩阵表示为：
[0087][0088][0089]
(3.3)确定状态约束观测方程：
[0090]
将惯导1、惯导2系统输出的速度、位置分别表示为：
[0091][0092][0093]
式中，和分别表示惯导1和惯导2输出的横平台坐标系下的速度矢量，惯导1与惯导2之间的外杆臂参数在系统安装完毕后通过标定补偿，表示惯导1输出的横计算坐标系下的位置信息，表示惯导2输出的横计算坐标系下的位置信息，表示惯导2体坐标系相对于横平台坐标系的旋转角速度，表示两套惯导间的外杆臂，表示两套惯导间的外杆臂在惯导2体坐标系下的投影，rc′
表示公共点的位置在横计算坐标系下的真值；
[0094]
由于两套系统反映的是同一载体的速度信息、位置信息，观测量实质上构成了惯导1、惯导2各自速度误差、位置误差的约束，将观测量表示为：
[0095][0096][0097]
式中，υv、υr为相应的速度观测噪声、位置观测噪声；
[0098]
增加外界高度信息的观测
[0099][0100]
式中，表示惯导1输出的高度信息，δh1表示惯导1的高度误差，υh为高度观测的噪声；将观测方程表示为：
[0101]
z(t)＝h(t)x(t)+υ(t)
[0102]
其中，
[0103][0104]
[0105]
h1＝[0 0 1]
[0106]
υ(t)＝[(υv)
t (υr)
t υh]
t
[0107]
式中，i3×3表示3行3列的单位矩阵；
[0108]
(4)确定两套惯导系统的转位次序：
[0109]
惯导1的转位次序为双轴16次序，具体转位流程如下：
[0110]
次序1：z轴以9
°
/s正向旋转180
°
，转停100s；
[0111]
次序2：y轴以9
°
/s反向旋转180
°
，转停100s；
[0112]
次序3：z轴以9
°
/s正向旋转180
°
，转停100s；
[0113]
次序4：y轴以9
°
/s反向旋转180
°
，转停100s；
[0114]
次序5：y轴以9
°
/s反向旋转180
°
，转停100s；
[0115]
次序6：z轴以9
°
/s正向旋转180
°
，转停100s；
[0116]
次序7：y轴以9
°
/s反向旋转180
°
，转停100s；
[0117]
次序8：z轴以9
°
/s正向旋转180
°
，转停100s；
[0118]
次序9：z轴以9
°
/s反向旋转180
°
，转停100s；
[0119]
次序10：y轴以9
°
/s正向旋转180
°
，转停100s；
[0120]
次序11：z轴以9
°
/s反向旋转180
°
，转停100s；
[0121]
次序12：y轴以9
°
/s正向旋转180
°
，转停100s；
[0122]
次序13：y轴以9
°
/s正向旋转180
°
，转停100s；
[0123]
次序14：z轴以9
°
/s反向旋转180
°
，转停100s；
[0124]
次序15：y轴以9
°
/s正向旋转180
°
，转停100s；
[0125]
次序16：z轴以9
°
/s反向旋转180
°
，转停100s；
[0126]
惯导2的转位次序为18次序，具体转位流程如下：
[0127]
次序1：y轴以9
°
/s正向旋转90
°
，转停180s；
[0128]
次序2：y轴以9
°
/s正向旋转180
°
，转停180s；
[0129]
次序3：y轴以9
°
/s正向旋转180
°
，转停180s；
[0130]
次序4：z轴以9
°
/s正向旋转90
°
，转停180s；
[0131]
次序5：z轴以9
°
/s正向旋转180
°
，转停180s；
[0132]
次序6：z轴以9
°
/s正向旋转180
°
，转停180s；
[0133]
次序7：x轴以9
°
/s正向旋转90
°
，转停180s；
[0134]
次序8：x轴以9
°
/s正向旋转180
°
，转停180s；
[0135]
次序9：x轴以9
°
/s正向旋转180
°
，转停180s；
[0136]
次序10：x轴以9
°
/s正向旋转90
°
，转停180s；
[0137]
次序11：x轴以9
°
/s正向旋转90
°
，转停180s；
[0138]
次序12：x轴以9
°
/s正向旋转90
°
，转停180s；
[0139]
次序13：z轴以9
°
/s正向旋转90
°
，转停180s；
[0140]
次序14：z轴以9
°
/s正向旋转90
°
，转停180s；
[0141]
次序15：z轴以9
°
/s正向旋转90
°
，转停180s；
[0142]
次序16：y轴以9
°
/s正向旋转90
°
，转停180s；
[0143]
次序17：y轴以9
°
/s正向旋转90
°
，转停180s；
[0144]
次序18：y轴以9
°
/s正向旋转90
°
，转停180s；
[0145]
基于联合转位方式，惯导1处于双轴旋转调制导航状态，惯导2处于标定状态，其陀螺标度因数误差、加速度计标度因数误差、安装误差均得到激励，根据步骤(3)所述方案，建立联合状态卡尔曼滤波器即实现惯导2的极区外场标定。
[0146]
进一步的，本发明方法对载体的运动状态无要求，载体处于系泊状态或运动状态均能实现在线标定。
[0147]
进一步的，本发明方法不仅仅适用于高纬度地区，在非横北极、横南极的中纬度地区同样适用。
[0148]
进一步的，惯导1与惯导2处于零位时的相对姿态在安装完毕后标定得到，惯导2在标定初始时刻的姿态基于两套惯导的相对姿态通过与惯导1传递对准获得。
[0149]
进一步的，步骤(4)所述的联合转位次序适用于两套及以上具有双轴转位机构的惯导系统间的在线标定，对于双轴和三轴惯导系统之间、多套三轴惯导系统间的在线标定也适用。
[0150]
进一步的，步骤(4)所述的联合转位次序仅为基于两套具有双轴转位机构的惯导系统的优选方案，对于其他旋转调制次序与标定次序之间的联合转位方案，也属于本发明的范畴。
[0151]
综上所述，本发明的优点及积极效果为：本发明建立了地球椭球模型横地球坐标系下基于psi角误差修正模型的联合状态卡尔曼滤波器，解决了传统的以当地水平坐标系为导航坐标系的惯性导航系统在极区失效的问题，通过两套惯导系统协同转位，利用两套惯导系统的冗余信息实现了在极地环境下的标定，psi角误差模型定义在计算坐标系，与位置误差解耦，更适合极地航行的长航时舰船惯导系统的标定，通过对速度误差的修正，消除比力项，提高了动基座下的标定精度，该标定方法不受极地特殊地理位置和恶劣自然环境的影响，具有重要工程意义。
附图说明
[0152]
图1是本发明实施例提供的流程图。
具体实施方式
[0153]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0154]
由于经线在极地快速收敛，以当地水平地理坐标系为导航坐标系的惯导系统会存在很大误差。由于磁暴、电离层的干扰等，在极区缺少可靠的外界参考信息，在极区的导航定位主要依赖于惯性导航系统。当惯导系统在极地环境下需要定期标定或者故障后更换了部分组件需要重新标定，传统的外场标定技术不能满足这一需求，这会影响惯导系统的精度。此外，传统的phi角误差模型定义在真实坐标系，然而真实坐标系是未知的，通常将计算坐标系近似为真实坐标系使用，存在近似误差，传统速度误差方程中的比力项也无法直接测得，在动基座下比力计算存在误差，这些因素会影响标定的精度。针对这些问题，本发明提出基于psi角误差修正模型的极地双惯导协同标定方法，所述的标定方法如图1所示。具
体实施方式如下：
[0155]
(1)构建基于地球椭球模型下的横地球坐标系和横计算坐标系；
[0156]
以北纬0
°
东经90
°
点为横地球坐标系下的北极点，定义为横北极点，北纬0
°
西经90
°
点为横地球坐标系下的南极点，定义为横南极点，0
°
经线和180
°
经线包围的椭圆面为横赤道面，取横北极、横南极和北极点组成的半个大椭圆为0
°
横经线，所在平面为横本初子午线，将地球坐标系与新定义的横地球坐标系之间的转换关系表示为：
[0157][0158]
式中，e系表示地球坐标系，e
′
系表示横地球坐标系，表示地球坐标系与横地球坐标系转换的方向余弦矩阵；
[0159]
基于横经纬网定义横计算坐标系，横北向指向横北极，所在位置的法线向上为天向，按照右手坐标系定义横东向，将横计算坐标系与计算坐标系之间的转换关系表示为：
[0160][0161]
式中，c系表示计算坐标系，c
′
系表示横计算坐标系，β表示计算坐标系与横计算坐标系之间的旋转角度；
[0162]
确定β与经纬度、横经纬度之间的转换关系：
[0163][0164][0165]
将横计算坐标系与横平台坐标系之间的转换关系表示为：
[0166][0167]
式中，p
′
系表示横平台坐标系，i表示单位矩阵，[ψ
×
]表示横计算坐标系下漂移误差角的反对称矩阵；
[0168]
将载体所在位置的法线与横赤道面的夹角定义为横纬度，与横本初子午面的夹角定义为横经度，将地球坐标系下定义的经纬度与横经纬度之间的转换关系表示为：
[0169][0170][0171]
[0172][0173]
式中，l
t
表示横纬度，λ
t
表示横经度，l表示纬度，λ表示经度；
[0174]
(2)构建两套惯导系统的误差模型；
[0175]
定义正常工作的双轴旋转调制惯导系统为惯导1，其体坐标系b1定义为“右-前-上”，待标定的惯导系统为惯导2，其体坐标系b2定义为“右-前-上”；
[0176]
惯导1的标度因数误差及安装误差很小进而忽略，将惯导1的误差模型定义为：
[0177][0178][0179]
其中，
[0180][0181][0182]
式中，表示惯导1的陀螺组件误差，表示惯导1的加速度计组件误差，表示惯导1的x轴陀螺漂移，表示惯导1的y轴陀螺漂移，表示惯导1的z轴陀螺漂移，表示惯导1的x轴加速度计零偏，表示惯导1的y轴加速度计零偏，表示惯导1的z轴加速度计零偏，表示惯导1的陀螺漂移，表示惯导1的加速度计零偏，为惯导1的陀螺噪声，为惯导1加速度计噪声；
[0183]
考虑标度因数误差、安装角误差及零偏误差，将惯导2的误差模型定义为：
[0184][0185][0186]
其中，
[0187][0188][0189]
式中，表示惯导2的陀螺组件误差，表示惯导2的加速度计组件误差，表示惯导2陀螺组件输出的理论角速度矢量，表示惯导2加速度计组件测得的理论比力矢量，表示惯导2的x轴陀螺漂移，表示惯导2的y轴陀螺漂移，表示惯导2的z轴陀螺漂移，表示惯导2的x轴加速度计零偏，表示惯导2的y轴加速度计零偏，表示惯导2的z轴加速度计零偏，表示惯导2的陀螺漂移，表示惯导2的加速度计零偏，为惯导2的陀螺噪声，为惯导2加速度计噪声；δκg和δμg表示陀螺的标度因数误差矩阵和安装误差矩阵，δκa和δμa表示加速度计的标度因数误差矩阵和安装误差矩阵；
[0190]
确定δκg和δκa：
[0191]
[0192][0193]
式中，δκ
gx
、δκ
gy
和δκ
gz
分别表示x轴陀螺、y轴陀螺和z轴陀螺的标度因数误差，δκ
ax
、δκ
ay
和δκ
az
分别表示x轴加速度计、y轴加速度计和z轴加速度计的标度因数误差；
[0194]
确定δμg和δμa：
[0195][0196][0197]
式中，δμ
gyx
、δμ
gzx
和δμ
gzy
表示陀螺组件的三个安装误差角，δμ
ayx
、δμ
azx
、δμ
azy
、δμ
axy
、δμ
axz
和δμ
ayz
表示加速度计组件的六个安装误差角；
[0198]
(3)利用两套惯导系统输出的姿态、速度、位置相关信息，联合状态卡尔曼滤波器，具体步骤为：
[0199]
(3.1)确定系统联合误差方程：
[0200][0201][0202][0203][0204][0205][0206]
其中，
[0207][0208][0209][0210]
式中，ψ1表示横计算坐标系下惯导1的漂移误差角，表示横平台坐标系下误差修正后惯导1的速度误差矢量，表示横计算坐标系下惯导1的位置误差，表示横计算坐标系下惯导1的地球自转角速度，表示横计算坐标系下惯导1的转移角速度，表示惯导1体坐标系至横平台坐标系的方向余弦矩阵，表示在横计算坐标系下惯导1的重力矢
量，ψ2表示横计算坐标系下惯导2的漂移误差角，表示横平台坐标系下误差修正后惯导2的速度误差矢量，表示横计算坐标系下惯导2的位置误差，表示惯导2横计算坐标系下的地球自转角速度，表示横计算坐标系下惯导2的转移角速度，表示惯导2体坐标系至横平台坐标系的方向余弦矩阵，表示在横计算坐标系下惯导2的重力矢量，v
p
′
表示载体在横平台坐标系下的速度，h表示载体所在位置的高度，r
x
为载体所在位置的横东向的曲率半径，ry为载体所在位置的横北向的曲率半径，为载体所在位置的扭曲率，re和rn分别为载体所在位置的卯酉圈半径和子午圈半径；
[0211]
(3.2)确定联合状态方程：
[0212][0213]
其中，
[0214][0215][0216][0217]
[0218][0219][0220][0221][0222][0223][0224]
[0225][0226][0227][0228][0229]
式中，0i×j表示i行j列的零矩阵，和为惯导1输出的横纬度和横经度，和为惯导2输出的横纬度和横经度，分别表示横平台坐标系下载体在横东向、横北向、天向的速度，ω
ie
表示地球自转角速度大小，表示惯导1输出位置处的重力加速度的值，表示惯导2输出位置处的重力加速度的值，c
23
表示矩阵的第二、第三列，c3表示矩阵的第三列，表示惯导2的x轴陀螺输出值，表示惯导2的y轴陀螺输出值，表示惯导2的z轴陀螺输出值，表示惯导2的x轴加速度计输出值，表示惯导2的y轴加速度计输出值，表示惯导2的z轴加速度计输出值；
[0230]
将状态向量x(t)表示为：
[0231][0232]
式中，ψ
e1
、ψ
n1
、ψ
u1
分别表示惯导1在横东向、横北向、天向的漂移误差角，分别表示惯导1在横东向、横北向、天向的漂移误差角，分别表示在横平台坐标系下惯导1在横东向、横北向、天向的速度误差，表示惯导1在横东向的位置误差，表示惯导1在横北向的位置误差，表示惯导1的天向误差，ψ
e2
、ψ
n2
、ψ
u2
表示惯导2在横东向、横北向、天向的漂移误差角，分别表示横平台坐标系下惯导2在横东向、横北向、天向的速度误差，表示惯导2的在横东向的位置误差，表示惯导2的在横北向的位置误差，表示惯导2在天向的位置误差；
[0233]
将噪声分布矩阵及噪声矩阵表示为：
[0234][0235][0236]
(3.3)确定状态约束观测方程：
[0237]
将惯导1、惯导2系统输出的速度、位置分别表示为：
[0238][0239][0240]
式中，和分别表示惯导1和惯导2输出的横平台坐标系下的速度矢量，惯导1与惯导2之间的外杆臂参数在系统安装完毕后通过标定补偿，表示惯导1输出的横计算坐标系下的位置信息，表示惯导2输出的横计算坐标系下的位置信息，表示惯导2体坐标系相对于横平台坐标系的旋转角速度，表示两套惯导间的外杆臂，表示两套惯导间的外杆臂在惯导2体坐标系下的投影，rc′
表示公共点的位置在横计算坐标系下的真值；
[0241]
由于两套系统反映的是同一载体的速度信息、位置信息，观测量实质上构成了惯导1、惯导2各自速度误差、位置误差的约束，将观测量表示为：
[0242][0243][0244]
式中，υv、υr为相应的速度观测噪声、位置观测噪声；
[0245]
增加外界高度信息的观测
[0246][0247]
式中，表示惯导1输出的高度信息，δh1表示惯导1的高度误差，υh为高度观测的噪声；
[0248]
将观测方程表示为：
[0249]
z(t)＝h(t)x(t)+υ(t)
[0250]
其中，
[0251][0252][0253]
h1＝[0 0 1]
[0254]
υ(t)＝[(υv)
t (υr)
t υh]
t
[0255]
式中，i3×3表示3行3列的单位矩阵；
[0256]
(4)确定两套惯导系统的转位次序：
[0257]
惯导1的转位次序为双轴16次序，具体转位流程如下：
[0258]
次序1：z轴以9
°
/s正向旋转180
°
，转停100s；
[0259]
次序2：y轴以9
°
/s反向旋转180
°
，转停100s；
[0260]
次序3：z轴以9
°
/s正向旋转180
°
，转停100s；
[0261]
次序4：y轴以9
°
/s反向旋转180
°
，转停100s；
[0262]
次序5：y轴以9
°
/s反向旋转180
°
，转停100s；
[0263]
次序6：z轴以9
°
/s正向旋转180
°
，转停100s；
[0264]
次序7：y轴以9
°
/s反向旋转180
°
，转停100s；
[0265]
次序8：z轴以9
°
/s正向旋转180
°
，转停100s；
[0266]
次序9：z轴以9
°
/s反向旋转180
°
，转停100s；
[0267]
次序10：y轴以9
°
/s正向旋转180
°
，转停100s；
[0268]
次序11：z轴以9
°
/s反向旋转180
°
，转停100s；
[0269]
次序12：y轴以9
°
/s正向旋转180
°
，转停100s；
[0270]
次序13：y轴以9
°
/s正向旋转180
°
，转停100s；
[0271]
次序14：z轴以9
°
/s反向旋转180
°
，转停100s；
[0272]
次序15：y轴以9
°
/s正向旋转180
°
，转停100s；
[0273]
次序16：z轴以9
°
/s反向旋转180
°
，转停100s；
[0274]
惯导2的转位次序为18次序，具体转位流程如下：
[0275]
次序1：y轴以9
°
/s正向旋转90
°
，转停180s；
[0276]
次序2：y轴以9
°
/s正向旋转180
°
，转停180s；
[0277]
次序3：y轴以9
°
/s正向旋转180
°
，转停180s；
[0278]
次序4：z轴以9
°
/s正向旋转90
°
，转停180s；
[0279]
次序5：z轴以9
°
/s正向旋转180
°
，转停180s；
[0280]
次序6：z轴以9
°
/s正向旋转180
°
，转停180s；
[0281]
次序7：x轴以9
°
/s正向旋转90
°
，转停180s；
[0282]
次序8：x轴以9
°
/s正向旋转180
°
，转停180s；
[0283]
次序9：x轴以9
°
/s正向旋转180
°
，转停180s；
[0284]
次序10：x轴以9
°
/s正向旋转90
°
，转停180s；
[0285]
次序11：x轴以9
°
/s正向旋转90
°
，转停180s；
[0286]
次序12：x轴以9
°
/s正向旋转90
°
，转停180s；
[0287]
次序13：z轴以9
°
/s正向旋转90
°
，转停180s；
[0288]
次序14：z轴以9
°
/s正向旋转90
°
，转停180s；
[0289]
次序15：z轴以9
°
/s正向旋转90
°
，转停180s；
[0290]
次序16：y轴以9
°
/s正向旋转90
°
，转停180s；
[0291]
次序17：y轴以9
°
/s正向旋转90
°
，转停180s；
[0292]
次序18：y轴以9
°
/s正向旋转90
°
，转停180s；
[0293]
基于联合转位方式，惯导1处于双轴旋转调制导航状态，惯导2处于标定状态，其陀螺标度因数误差、加速度计标度因数误差、安装误差均得到激励，根据步骤(3)所述方案，建立联合状态卡尔曼滤波器即实现惯导2的极区外场标定。
[0294]
本发明方法对载体的运动状态无要求，载体处于系泊状态或运动状态均能实现在线标定。
[0295]
本发明方法不仅仅适用于高纬度地区，在非横北极、横南极的中纬度地区同样适用。
[0296]
惯导1与惯导2处于零位时的相对姿态在安装完毕后标定得到，惯导2在标定初始时刻的姿态基于两套惯导的相对姿态通过与惯导1传递对准获得。
[0297]
步骤(4)所述中的联合转位次序适用于两套及以上具有双轴转位机构的惯导系统间的在线标定，对于双轴和三轴惯导系统之间、多套三轴惯导系统间的在线标定也适用。
[0298]
步骤(4)所述中的联合转位次序仅为基于两套具有双轴转位机构的惯导系统的优选方案，对于其他旋转调制次序与标定次序之间的联合转位方案，也属于本发明的范畴。
[0299]
以上所述仅是本发明的优选实施方式，并不用以限制本发明，凡属于本发明思路下的技术方案均属于本发明的保护范围。在不脱离本发明原理前提下的若干改进和润饰等，这些改进和润饰也应视为本发明的保护范围。

技术特征：
1.基于psi角误差修正模型的极地双惯导协同标定方法，其特征在于，所述方法包括以下步骤：(1)构建基于地球椭球模型下的横地球坐标系和横计算坐标系；以北纬0
°
东经90
°
点为横地球坐标系下的北极点，定义为横北极点，北纬0
°
西经90
°
点为横地球坐标系下的南极点，定义为横南极点，0
°
经线和180
°
经线包围的椭圆面为横赤道面，取横北极、横南极和北极点组成的半个大椭圆为0
°
横经线，所在平面为横本初子午线，将地球坐标系与新定义的横地球坐标系之间的转换关系表示为：式中，e系表示地球坐标系，e
′
系表示横地球坐标系，表示地球坐标系与横地球坐标系转换的方向余弦矩阵；基于横经纬网定义横计算坐标系，横北向指向横北极，所在位置的法线向上为天向，按照右手坐标系定义横东向，将横计算坐标系与计算坐标系之间的转换关系表示为：式中，c系表示计算坐标系，c
′
系表示横计算坐标系，β表示计算坐标系与横计算坐标系之间的旋转角度；确定β与经纬度、横经纬度之间的转换关系：确定β与经纬度、横经纬度之间的转换关系：将横计算坐标系与横平台坐标系之间的转换关系表示为：式中，p
′
系表示横平台坐标系，i表示单位矩阵，[ψ
×
]表示横计算坐标系下漂移误差角的反对称矩阵；将载体所在位置的法线与横赤道面的夹角定义为横纬度，与横本初子午面的夹角定义为横经度，将地球坐标系下定义的经纬度与横经纬度之间的转换关系表示为：为横经度，将地球坐标系下定义的经纬度与横经纬度之间的转换关系表示为：为横经度，将地球坐标系下定义的经纬度与横经纬度之间的转换关系表示为：
式中，l
t
表示横纬度，λ
t
表示横经度，l表示纬度，λ表示经度；(2)构建两套惯导系统的误差模型；定义正常工作的双轴旋转调制惯导系统为惯导1，其体坐标系b1定义为“右-前-上”，待标定的惯导系统为惯导2，其体坐标系b2定义为“右-前-上”；惯导1的标度因数误差及安装误差很小进而忽略，将惯导1的误差模型定义为：惯导1的标度因数误差及安装误差很小进而忽略，将惯导1的误差模型定义为：其中，其中，式中，表示惯导1的陀螺组件误差，表示惯导1的加速度计组件误差，表示惯导1的x轴陀螺漂移，表示惯导1的y轴陀螺漂移，表示惯导1的z轴陀螺漂移，表示惯导1的x轴加速度计零偏，表示惯导1的y轴加速度计零偏，表示惯导1的z轴加速度计零偏，表示惯导1的陀螺漂移，表示惯导1的加速度计零偏，为惯导1的陀螺噪声，为惯导1加速度计噪声；考虑标度因数误差、安装角误差及零偏误差，将惯导2的误差模型定义为：考虑标度因数误差、安装角误差及零偏误差，将惯导2的误差模型定义为：其中，其中，式中，表示惯导2的陀螺组件误差，表示惯导2的加速度计组件误差，表示惯导2陀螺组件输出的理论角速度矢量，表示惯导2加速度计组件测得的理论比力矢量，表示惯导2的x轴陀螺漂移，表示惯导2的y轴陀螺漂移，表示惯导2的z轴陀螺漂移，表示惯导2的x轴加速度计零偏，表示惯导2的y轴加速度计零偏，表示惯导2的z轴加速度计零偏，表示惯导2的陀螺漂移，表示惯导2的加速度计零偏，为惯导2的陀螺噪声，为惯导2加速度计噪声；δκ
g
和δμ
g
表示陀螺的标度因数误差矩阵和安装误差矩阵，δκ
a
和δμ
a
表示加速度计的标度因数误差矩阵和安装误差矩阵；确定δκ
g
和δκ
a
：
表示横计算坐标系下惯导2的漂移误差角，表示横平台坐标系下误差修正后惯导2的速度误差矢量，表示横计算坐标系下惯导2的位置误差，表示惯导2横计算坐标系下的地球自转角速度，表示横计算坐标系下惯导2的转移角速度，表示惯导2体坐标系至横平台坐标系的方向余弦矩阵，表示在横计算坐标系下惯导2的重力矢量，v
p
′
表示载体在横平台坐标系下的速度，h表示载体所在位置的高度，r
x
为载体所在位置的横东向的曲率半径，r
y
为载体所在位置的横北向的曲率半径，为载体所在位置的扭曲率，r
e
和r
n
分别为载体所在位置的卯酉圈半径和子午圈半径；(3.2)确定联合状态方程：其中，其中，其中，其中，
式中，0
i
×
j
表示i行j列的零矩阵，和为惯导1输出的横纬度和横经度，和为惯导2输出的横纬度和横经度，分别表示横平台坐标系下载体在横东向、横北向、天向的速度，ω
ie
表示地球自转角速度大小，表示惯导1输出位置处的重力加速度的值，表示惯导2输出位置处的重力加速度的值，c
23
表示矩阵的第二、第三列，c3表示矩阵的第三列，表示惯导2的x轴陀螺输出值，表示惯导2的y轴陀螺输出值，表示惯导2的z轴陀螺输出值，表示惯导2的x轴加速度计输出值，表示惯导2的y轴加速度计输出值，表示惯导2的z轴加速度计输出值；将状态向量x(t)表示为：式中，ψ
e1
、ψ
n1
、ψ
u1
分别表示惯导1在横东向、横北向、天向的漂移误差角，分别表示惯导1在横东向、横北向、天向的漂移误差角，分别表示在横平台坐标系下惯导1在横东向、横北向、天向的速度误差，表示惯导1在横东向的位置误差，表示惯导1在横北向的位置误差，表示惯导1的天向误差，ψ
e2
、ψ
n2
、ψ
u2
表示惯导2在横东向、横北向、天向的漂移误差角，分别表示横平台坐标系下惯导2在横东向、横北向、天向的速度误差，表示惯导2的在横东向的位置误差，表示惯导2的在横北向的位置误差，表示惯导2在天向的位置误差；将噪声分布矩阵及噪声矩阵表示为：
(3.3)确定状态约束观测方程：将惯导1、惯导2系统输出的速度、位置分别表示为：将惯导1、惯导2系统输出的速度、位置分别表示为：式中，和分别表示惯导1和惯导2输出的横平台坐标系下的速度矢量，惯导1与惯导2之间的外杆臂参数在系统安装完毕后通过标定补偿，表示惯导1输出的横计算坐标系下的位置信息，表示惯导2输出的横计算坐标系下的位置信息，表示惯导2体坐标系相对于横平台坐标系的旋转角速度，表示两套惯导间的外杆臂，表示两套惯导间的外杆臂在惯导2体坐标系下的投影，r
c
′
表示公共点的位置在横计算坐标系下的真值；由于两套系统反映的是同一载体的速度信息、位置信息，观测量实质上构成了惯导1、惯导2各自速度误差、位置误差的约束，将观测量表示为：惯导2各自速度误差、位置误差的约束，将观测量表示为：式中，υ
v
、υ
r
为相应的速度观测噪声、位置观测噪声；增加外界高度信息的观测增加外界高度信息的观测式中，表示惯导1输出的高度信息，δh1表示惯导1的高度误差，υ
h
为高度观测的噪声；将观测方程表示为：z(f)＝h(f)x(f)+υ(f)其中，其中，h1＝[0 0 1]
υ(t)＝[(υ
v
)
t (υ
r
)
t υ
h
]
t
式中，i3×3表示3行3列的单位矩阵；(4)确定两套惯导系统的转位次序：惯导1的转位次序为双轴16次序，具体转位流程如下：次序1：z轴以9
°
/s正向旋转180
°
，转停100s；次序2：y轴以9
°
/s反向旋转180
°
，转停100s；次序3：z轴以9
°
/s正向旋转180
°
，转停100s；次序4：y轴以9
°
/s反向旋转180
°
，转停100s；次序5：y轴以9
°
/s反向旋转180
°
，转停100s；次序6：z轴以9
°
/s正向旋转180
°
，转停100s；次序7：y轴以9
°
/s反向旋转180
°
，转停100s；次序8：z轴以9
°
/s正向旋转180
°
，转停100s；次序9：z轴以9
°
/s反向旋转180
°
，转停100s；次序10：y轴以9
°
/s正向旋转180
°
，转停100s；次序11：z轴以9
°
/s反向旋转180
°
，转停100s；次序12：y轴以9
°
/s正向旋转180
°
，转停100s；次序13：y轴以9
°
/s正向旋转180
°
，转停100s；次序14：z轴以9
°
/s反向旋转180
°
，转停100s；次序15：y轴以9
°
/s正向旋转180
°
，转停100s；次序16：z轴以9
°
/s反向旋转180
°
，转停100s；惯导2的转位次序为18次序，具体转位流程如下：次序1：y轴以9
°
/s正向旋转90
°
，转停180s；次序2：y轴以9
°
/s正向旋转180
°
，转停180s；次序3：y轴以9
°
/s正向旋转180
°
，转停180s；次序4：z轴以9
°
/s正向旋转90
°
，转停180s；次序5：z轴以9
°
/s正向旋转180
°
，转停180s；次序6：z轴以9
°
/s正向旋转180
°
，转停180s；次序7：x轴以9
°
/s正向旋转90
°
，转停180s；次序8：x轴以9
°
/s正向旋转180
°
，转停180s；次序9：x轴以9
°
/s正向旋转180
°
，转停180s；次序10：x轴以9
°
/s正向旋转90
°
，转停180s；次序11：x轴以9
°
/s正向旋转90
°
，转停180s；次序12：x轴以9
°
/s正向旋转90
°
，转停180s；次序13：z轴以9
°
/s正向旋转90
°
，转停180s；次序14：z轴以9
°
/s正向旋转90
°
，转停180s；次序15：z轴以9
°
/s正向旋转90
°
，转停180s；次序16：y轴以9
°
/s正向旋转90
°
，转停180s；次序17：y轴以9
°
/s正向旋转90
°
，转停180s；次序18：y轴以9
°
/s正向旋转90
°
，转停180s；
基于联合转位方式，惯导1处于双轴旋转调制导航状态，惯导2处于标定状态，其陀螺标度因数误差、加速度计标度因数误差、安装误差均得到激励，根据步骤(3)所述方案，建立联合状态卡尔曼滤波器即实现惯导2的极区外场标定。2.如权利要求1所述的基于psi角误差修正模型的极地双惯导协同标定方法，其特征在于，所述方法对载体的运动状态无要求，载体处于系泊状态或运动状态均能实现在线标定。3.如权利要求1所述的基于psi角误差修正模型的极地双惯导协同标定方法，其特征在于，所述方法不仅仅适用于高纬度地区，在非横北极、横南极的中纬度地区同样适用。4.如权利要求1所述的基于psi角误差修正模型的极地双惯导协同标定方法，其特征在于，惯导1与惯导2处于零位时的相对姿态在安装完毕后标定得到，惯导2在标定初始时刻的姿态基于两套惯导的相对姿态通过与惯导1传递对准获得。5.如权利要求1所述的基于psi角误差修正模型的极地双惯导协同标定方法，其特征在于，所述步骤(4)中的联合转位次序适用于两套及以上具有双轴转位机构的惯导系统间的在线标定，对于双轴和三轴惯导系统之间、多套三轴惯导系统间的在线标定也适用。

技术总结
本发明属于导航技术领域，公开了基于Psi角误差修正模型的极地双惯导协同标定方法，适用于极地环境无外界参考信息时双惯导系统的协同标定。本发明在椭球模型横地球坐标系的基础上，将两套惯导系统在横平台坐标系下的相对速度、横计算坐标系下的相对位置作为约束观测，建立了Psi角误差修正模型下的双惯导系统联合状态卡尔曼滤波器，对待标定惯导系统的误差参数进行标定估计。本发明方法标定精度不受正常工作惯导系统绝对误差的影响，在运动状态下也能进行标定；通过对速度误差模型的修正，避免动态情况下比力不准影响标定精度；Psi角误差模型定义在计算坐标系，与位置误差解耦，更适合极地航行的长航时舰船惯导系统的标定。更适合极地航行的长航时舰船惯导系统的标定。更适合极地航行的长航时舰船惯导系统的标定。


技术研发人员：王林 梁钟泓 廖志坤 郭鸿刚 王元涵 郝艳 陈宏翔 宋东阳
受保护的技术使用者：中国人民解放军国防科技大学
技术研发日：2023.03.11
技术公布日：2023/7/26