一种旋转捷联惯导系统误差系统级标定方法

1.本发明涉及捷联惯性导航技术领域，具体涉及一种旋转捷联惯导系统误差系统级标定方法。


背景技术：

2.惯性导航系统的主要组成部分为惯性测量单元，在生产制造时不可避免的存在相关器件误差，如常值误差、刻度系数误差、安装误差等。为了抑制这些误差对惯导系统的影响，除了旋转调制技术外，常常通过误差标定技术对误差进行估计补偿，来达到高精度导航的目的。标定主要分为分立式标定和系统级标定。分立式标定技术将惯导系统安装在转动机构上，通过不同位置编排对陀螺仪和加速度计进行标定，将标定出的误差对惯导系统进行补偿，但分立式标定具有拆装繁琐，并且对转动机构要求精度高，无法外场标定的缺点。旋转捷联惯性导航系统可以利用自身的转动机构进行系统级标定，可以避免从载体上拆装的复杂性并且可以为长航时精确导航提供了现实基础。
3.系统级标定通过设计合理的位置编排来激励惯性测量单元的误差，建立惯性测量单元输入输出数学模型，以导航误差(速度误差、位置误差)作为观测量，通过最小二乘法或卡尔曼滤波将误差参数估计出来。sagem公司设计了一种18位置标定方案(camberlein l,mazzanti f.calibration technique for laser gyro strapdown inertial navigation systems[j].ortung und navigation,1985:5.0-5.13.)用来标定激光陀螺惯导系统，经过标定可以达到系统级精度要求，并且该方案有着丰富的工程实践经验，是一种传统的误差标定方法。文献(谢波,秦永元,万彦辉.激光陀螺捷联惯导系统多位置标定方法[j].中国惯性技术学报,2011,19(02):157-162+169.)提出了一种19位置标定方法，通过多次初始对准、位置翻转和静态导航过程激励惯性测量单元误差参数，采用最小二乘方法可以在较短时间内可以标定出21个误差参数。系统级标定技术中旋转标定方案是关键因素。


技术实现要素：

[0004]
本发明的目的是解决惯导系统误差标定问题，提供了一种新的旋转捷联惯导系统误差标定方法，能够同时标定出陀螺仪和加速度计的常值误差、刻度系数误差和安装误差，并对惯导系统进行补偿，进而达到提高系统精度的目的。
[0005]
为实现上述目的，本发明所设计的旋转捷联惯导系统误差标定方法，包括如下步骤：
[0006]
步骤s11：定义坐标系和定义惯性测量单元中陀螺仪和加速度计的常值误差参数、刻度系数误差参数和安装误差参数
[0007]
确定旋转惯导系统惯性测量单元所处的惯性坐标系(设为i系)，惯性坐标系的原点oi位于地球的中心，由原点出发指向北极为o
i-zi轴，由原点出发指向平均春分点为o
i-xi轴，o
i-yi轴与o
i-zi轴构成右手直角坐标系；
[0008]
确定旋转惯导系统惯性测量单元所处的导航坐标系(设为n系)，该坐标系为当地
地理坐标系，其原点on位于载体质心，其o
n-xn、o
n-yn、o
n-zn轴由原点分别指向东、北、天方向；
[0009]
确定旋转惯导系统惯性测量单元所处的载体坐标系(设为b系)，载体的质心为载体坐标系的原点ob，载体系坐标系的o
b-xb、o
b-yb和o
b-zb轴分别从载体的质心指向载体的右侧、前侧和上侧；
[0010]
确定惯性测量单元坐标系(设为p系)，惯性测量单元由三组正交安装的陀螺仪和加速度计组成，惯性测量单元坐标系的原点o
p
位于惯性测量单元质心，其o
p-x
p
、o
p-y
p
、o
p-z
p
轴指向由原点出发分别与三组正交安装的陀螺仪和加速度计指向平行；
[0011]
确定安装坐标系坐标系(设为m系)，在制作惯性测量单元时，无法保证三组陀螺仪和加速度计是理想正交安装的，设实际安装时惯性测量单元坐标系为安装坐标系，安装坐标系的原点om与p系相同；
[0012]
陀螺仪和加速度计常值误差分别表示为ε和具体表示如下：
[0013][0014]
其中，ε
x
、εy、εz分别为陀螺仪o
p-x
p
、o
p-y
p
、o
p-z
p
轴的常值误差，分别为加速度计常值误差在o
p-x
p
、o
p-y
p
、o
p-z
p
轴的分量。
[0015]
陀螺仪和加速度计刻度系数误差分别表示为δkg和δka，具体表示如下：
[0016][0017]
其中，δk
gx
、δk
gy
、δk
gz
分别为陀螺仪刻度系数误差在o
p-x
p
、o
p-y
p
、o
p-z
p
轴的分量，δk
ax
、δk
ay
、δk
az
分别为加速度计刻度系数误差在o
p-x
p
、o
p-y
p
、o
p-z
p
轴的分量。
[0018]
陀螺仪和加速度计安装误差分别表示为δag和δaa，具体表示如下：
[0019][0020]
其中，δa
gyx
为o
m-ym轴陀螺仪与o
p-x
p-y
p
平面的安装误差角，δa
gzx
为o
m-zm轴陀螺仪与o
p-x
p-z
p
平面的安装误差角，δa
gzy
为o
m-zm轴陀螺仪与o
p-y
p-z
p
平面的安装误差角，δa
gxy
为o
m-xm轴陀螺仪与o
p-x
p-y
p
平面的安装误差角，δa
gxz
为o
m-xm轴陀螺仪与o
p-x
p-z
p
平面的安装误差角，δa
gyz
为o
m-ym轴陀螺仪与o
p-y
p-z
p
平面的安装误差角，δa
ayx
为o
m-ym轴加速度计与o
p-x
p-y
p
平面的安装误差角，δa
azx
为o
m-zm轴加速度计与o
p-x
p-z
p
平面的安装误差角，δa
azy
为o
m-zm轴加速度计与o
p-y
p-z
p
平面的安装误差角，δa
axy
为o
m-xm轴加速度计与o
p-x
p-y
p
平面的安装误差角，δa
axz
为o
m-xm轴加速度计与o
p-x
p-z
p
平面的安装误差角，δa
ayz
为o
m-ym轴加速度计与o
p-y
p-z
p
平面的安装误差角。
[0021]
步骤s21:构建陀螺仪和加速度计误差模型
[0022]
为了减少标定参数的数量，以便于提高运算速率，缩短标定时间，现假定陀螺仪正交系的o
p-x
p
轴与安装坐标系下的o
m-xm轴重合，并且o
p-y
p
轴在o
m-x
m-ym平面内。则陀螺仪误差模型如下：
[0023][0024]
其中，分别为陀螺仪角增量误差在o
b-xb、o
b-yb、o
b-zb轴方向的分量，分别为陀螺仪角增量在o
b-xb、o
b-yb、o
b-zb轴方向的分量。
[0025]
加速度计误差模型：
[0026][0027]
其中，分别为加速度计比力增量误差在o
b-xb、o
b-yb、o
b-zb轴方向的分量，分别为加速度计比力增量在o
b-xb、o
b-yb、o
b-zb轴方向的分量。
[0028]
步骤s31:通过步骤21陀螺仪和加速度计误差模型经过翻转激励出步骤11定义的误差参数。
[0029]
初始时惯性测量单元的三个轴向o
p-x
p
、o
p-y
p
、o
p-z
p
分别指向天-东-北，转动机构的中轴与惯性测量单元的o
p-y
p
轴重合，转动机构的内轴与惯性测量单元的o
p-z
p
轴重合，转动机构的外轴与惯性测量单元的o
p-x
p
轴重合。设计30位置标定路径如下：1.绕转动机构中框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；2.绕转动机构中框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；3.绕转动机构中框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；4.绕转动机构中框轴反向旋转180度，旋转角速度为ω度/秒，转动完成后停止t秒；5.绕转动机构中框轴反向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；6.绕转动机构内框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；7.绕转动机构中框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；8.绕转动机构中框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；9.绕转动机构中框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；10.绕转动机构中框轴反向旋转180度，旋转角速度为ω度/秒，转动完成后停止t秒；11.绕转动机构中框轴反向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；12.绕转动机构内框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；13.绕转动机构内框轴正向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；14.绕转动机构内框轴反向旋转180度，旋转角速度为ω度/秒，转动完成后停止t秒；15.绕转动机构内框轴反向旋转90度，旋转角速度为ω度/秒，转动完成后停止t秒；16-30位置与1-15位置转轴次序相同，转动方向相反。
[0030]
步骤s41:根据步骤31中的误差参数建立kalman滤波模型
[0031]
经过步骤31的一系列翻转，可以激励陀螺仪和加速度计产生误差。然后将其带入惯导姿态误差、速度误差和位置误差方程即可得到系统状态误差模型。
[0032]
建立状态向量包含姿态误差、速度误差、位置误差、陀螺仪和加速度计常值误差、刻度系数误差和安装误差的30维kalman状态方程：
[0033][0034]
其中，30维状态向量x(t)为
[0035][0036]
为状态向量x的微分，f为状态转移矩阵，表达式如下：
[0037][0038]
其中各分块矩阵具体表达式如下：
[0039][0040]
其中，ω
ie
为地球自转角速率，rm、rn分别为地球子午圈、卯酉圈曲率半径，l为载体所在纬度，h为载体所在高度。
[0041]
[0042][0043]
其中，ve、vn分别为北向、北向速度。
[0044][0045]
其中，分别为导航坐标系下加速度计输出在东向、北向和天向的分量。
[0046][0047]
其中，vu为天向速度。
[0048][0049]
[0050][0051]
从载体坐标系到导航坐标系的姿态变换矩阵如下：
[0052][0053][0054][0055]
式(3)中g为系统噪声驱动矩阵，表达式如下：
[0056][0057]
w(t)为系统噪声矩阵，由陀螺仪噪声wg(t)和加速度计噪声wa(t)组成，均为零均值高斯白噪声向量序列，其表达式如下：
[0058][0059]
建立如下以速度误差和位置误差为观测量的kalman滤波器观测方程：
[0060]
z(t)＝h(t)x(t)+v(t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0061]
其中，z(t)为观测量，由速度误差和位置误差组成，表达式如下：
[0062]
z(t)＝[δv
e δv
n δv
u δl δλ δh]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0063]
h为观测矩阵，表达式如下：
[0064][0065]
v(t)为观测噪声矩阵，代表能够获取观测量的系统噪声，包括速度观测噪声和位置观测噪声，均为零均值高斯白噪声向量序列，且上述w(t)和v(t)互不相关，表达式如下：
[0066]
[0067]
步骤s51：kalman滤波模型离散化获得陀螺仪和加速度计的常值误差、刻度系数误差和安装误差估计值
[0068]
给定系统状态空间模型：
[0069][0070]
式中，xk为离散后的状态向量，zk为离散后的量测向量，f
k/k-1
、hk和γ
k/k-1
分别为状态转移矩阵、量测矩阵和系统噪声驱动矩阵，w
k-1
和vk分别为系统噪声向量和量测噪声向量。
[0071]
从k-1时刻到k时刻的离散状态转移矩阵f
k/k-1
具体表达式如下：
[0072][0073]
其中，i为单位矩阵，ts为滤波周期。
[0074]
kalman滤波全套算法如下：
[0075]
状态一步预测
[0076][0077]
状态一步预测均方误差阵
[0078][0079]
滤波增益
[0080][0081]
状态估计
[0082][0083]
状态估计均方误差
[0084]
pk＝(i-k
khk
)p
k/k-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(28)
[0085]
式(25)中q为陀螺仪和加速度计输出白噪声的均方误差矩阵，r为观测噪声方差矩阵。
[0086]
综上所述，结合步骤31一系列翻转运动，将惯性测量单元输出数据进行惯导解算，与此同时进行kalman滤波，在给定系统状态向量初始值x0、状态估计均方差初始值p0，陀螺仪和加速度计输出白噪声的均方误差矩阵q0和观测噪声方差矩阵初始值r0的条件下，以及确定观测量z、状态转移矩阵f和观测矩阵h的选取，由式(24)至式(28)可以递推估计出任意时刻的状态量x，从而可以得到陀螺仪和加速度计的常值误差、刻度系数误差和安装误差估计值。
[0087]
步骤s61:将步骤51中得到的陀螺仪和加速度计的常值误差、刻度系数误差和安装误差估计值对系统进行校正。误差修正模型如下所示：
[0088]
陀螺仪和加速度计常值误差反馈：
[0089][0090][0091]
其中，ε和为反馈后系统常值误差，和为反馈前系统常值误差，和常值误差估计值。
[0092]
陀螺仪刻度系数误差和安装误差反馈：
[0093][0094]
其中kg和ag为陀螺仪反馈后刻度系数和安装角矩阵，和为陀螺仪反馈前刻度系数和安装角矩阵，δkg和δag为陀螺仪刻度系数误差和安装误差估计值。
[0095]
陀螺仪和加速度安装误差反馈：
[0096][0097]
其中ka和aa为加速度计反馈后刻度系数和安装角矩阵，和为加速度计反馈前刻度系数和安装角矩阵，δka和δaa为加速度计刻度系数误差和安装误差估计值。
[0098]
与现有技术相比，本发明的有益效果为：本发明利用旋转惯导系统自身转动机构优势以及设计的30位置旋转标定方案，可以在外场环境下充分激励出惯导相关器件误差，然后通过建立30维kalman滤波模型，以速度误差和位置误差作为观测量可以较高精度的估计出陀螺仪和加速度计的常值误差、刻度系数误差和安装误差。
附图说明
[0099]
图1为转动机构示意图；
[0100]
图2为旋转标定策略图；
[0101]
图3为仿真实验误差参数估计曲线对比
[0102]
图4为估计的惯性器件误差补偿后导航定位误差对比
具体实施方式
[0103]
下面结合附图和具体实施例对本发明作进一步的详细说明。
[0104]
本发明可行性可通过仿真实验验证，如图1惯性测量单元安装在旋转机构上，b系下标定旋转方案如图2所示，仿真实验采取对比实验，以sagem公司设计的传统18位置方案为参照，比对本发明方案的标定效果：
[0105]
实验条件：初始纬度30.6度，初始经度114.2度，初始高度25米；初始姿态为(0
°
,-90
°
,-90
°
)；仿真采取静态仿真，所以初始速度为0；初始姿态误差(0.5
′
,-0.5
′
,0.5
′
)；初始速度误差(0.1m/s,0.1m/s,0.1m/s,)；初始位置误差(1m,1m,3m)；仿真时长4小时。
[0106]
kalman滤波模型初始参数设置：
[0107]
系统状态向量初始值x0＝[0
30
×1]；
[0108]
状态估计均方差初始值p0：
[0109]
p0＝diag{(0.5
′
)2,(-0.5
′
)2,(0.5
′
)2,(0.1m/s)2,(0.1m/s)2,(0.1m/s)2,(1m)2,(1m)2,(3m)2,(0.05
°
/h)2,(0.05
°
/h)2,(0.05
°
/h)2,(80ug)2,(90ug)2,(100ug)2,(40ppm)2,(10
″
)2,(10
″
)2,(40ppm)2,(10
″
)2,(40ppm)2,(20ppm)2,(10
″
)2,(10
″
)2,(10
″
)2,(20ppm)2,(10
″
)2,(10
″
)2,(10
″
)2,(20ppm)2}
[0110]
白噪声的均方误差矩阵q0：
[0111]
q0＝diag{(0.001
°
/h)2,(0.001
°
/h)2,(0.001
°
/h)2,(1ug)2,(1ug)2,(1ug)2}
[0112]
观测噪声方差矩阵初始值r0：
[0113]
r0＝diag{(0.1m/s)2,(0.1m/s)2,(0.1m/s)2,(1m)2,(1m)2,(1m)2}
[0114]
仿真实验中陀螺仪和加速度计常值误差、安装误差、刻度系数误差参数设定值与估计值如下表1所示。
[0115]
表1标定误差参数设定值与估计值结果
[0116][0117]
如附图3中3a陀螺仪x轴常值误差、3b陀螺仪y轴常值误差、3c陀螺仪z轴常值误差、3d加速度计x轴常值误差、3e加速度计y轴常值误差、3f加速度计z轴常值误差、3g陀螺仪x轴刻度系数误差、3h陀螺仪y轴刻度系数误差、3i陀螺仪z轴刻度系数误差、3j加速度计x轴刻度系数误差、3k加速度计y轴刻度系数误差、3l加速度计z轴刻度系数误差、3m陀螺仪安装误差yx、3n陀螺仪安装误差zx、3o陀螺仪安装误差zy、3p加速度计安装误差yx、3q加速度计安装误差zx、3r加速度计安装误差xy、3s加速度计安装误差zy、3t加速度计安装误差xz、3u加速度计安装误差yz所示为陀螺仪和加速度计21个相关误差参数估计曲线。从图中可以看出，设计的30位置方案各项误差参数在标定旋转策略以及kalman滤波的作用下和传统18位置方案一样得到了收敛，且大部分误差参数收敛效果更好。
[0118]
结合表1结果，设计的30位置标定方案所估计出的21项误差参数中16项残余误差比18位置方案小。将两种标定方案估计的惯性器件误差对imu补偿后进行48小时的静基座
导航，两种方案的导航定位误差曲线如附图4，从图中可以看出，设计的30位置标定方案的最大定位误差比传统18位置方案小3639米，可以验证本发明所提出的标定方案具有更高的标定精度。