首先设定如下坐标系:地心地固坐标系(e系)、载体坐标系(b系)、初始时刻的与e系重合并凝滞的坐标系e0系以及初始时刻与b系重合并凝滞的坐标系b0系,则b系到e系的方向余弦阵可做如下链式分解:
C b e = C e 0 e ∗ C b 0 e 0 ∗ C b b 0 C_b^e = C_{e0}^{e} * C_{b0}^{e0} * C_b^{b0} Cbe=Ce0e∗Cb0e0∗Cbb0
因为:
C ˙ e e 0 = C e e 0 ( ω i e e × ) = C e e 0 Ω i e e \dot C_{e}^{e0} = C_{e}^{e0} (\omega_{ie}^e ×) = C_{e}^{e0} \Omega_{ie}^e C˙ee0=Cee0(ωiee×)=Cee0Ωiee
可得:
C e e 0 = I + s i n ( ω i e t ) ω i e t Ω i e e t + 1 − c o s ( ω i e t ) ( ω i e t ) 2 ( Ω i e e t ) 2 C_{e}^{e0} = I + \frac{sin(\omega_{ie}t)}{\omega_{ie}t}\Omega_{ie}^e t + \frac{1-cos(\omega_{ie}t)}{(\omega_{ie}t)^2}(\Omega_{ie}^e t)^2 Cee0=I+ωietsin(ωiet)Ωieet+(ωiet)21−cos(ωiet)(Ωieet)2
其中:
ω
i
e
=
∣
ω
i
e
e
∣
\omega_{ie} = |\omega_{ie}^e|
ωie=∣ωiee∣,t为从开始到现在的时间长度。
因为:
C ˙ b b 0 = C b b 0 ( ω b 0 b b × ) = C b b 0 ( ω i b b × ) = C b b 0 Ω i b b \dot C_{b}^{b0} = C_{b}^{b0} (\omega_{b0b}^b ×) = C_{b}^{b0} (\omega_{ib}^b ×) = C_{b}^{b0} \Omega_{ib}^b C˙bb0=Cbb0(ωb0bb×)=Cbb0(ωibb×)=Cbb0Ωibb
可得:
C b b 0 = I + s i n ( ω i b t ) ω i b t Ω i b b t + 1 − c o s ( ω i b t ) ( ω i b t ) 2 ( Ω i b b t ) 2 C_{b}^{b0} = I + \frac{sin(\omega_{ib}t)}{\omega_{ib}t}\Omega_{ib}^b t + \frac{1-cos(\omega_{ib}t)}{(\omega_{ib}t)^2}(\Omega_{ib}^b t)^2 Cbb0=I+ωibtsin(ωibt)Ωibbt+(ωibt)21−cos(ωibt)(Ωibbt)2
其中: ω i b = ∣ ω i b b ∣ \omega_{ib} = |\omega_{ib}^b| ωib=∣ωibb∣,t为从开始到现在的时间长度。但是应用中:
C b b 0 ( t ) = C b b 0 ( t − 1 ) [ I + s i n ( ω i b Δ t ) ω i b Δ t Ω i b b Δ t + 1 − c o s ( ω i b Δ t ) ( ω i b Δ t ) 2 ( Ω i b b Δ t ) 2 ] C_{b}^{b0} (t) = C_{b}^{b0} (t-1) [I + \frac{sin(\omega_{ib} \Delta t)}{\omega_{ib} \Delta t}\Omega_{ib}^b \Delta t + \frac{1-cos(\omega_{ib} \Delta t)}{(\omega_{ib} \Delta t)^2}(\Omega_{ib}^b \Delta t)^2] Cbb0(t)=Cbb0(t−1)[I+ωibΔtsin(ωibΔt)ΩibbΔt+(ωibΔt)21−cos(ωibΔt)(ΩibbΔt)2]
此时等效旋转矢量也为 Δ t \Delta t Δt时间内的。不需限制 ω i b b \omega_{ib}^{b} ωibb大小,具有很强的抗角运动的能力;
下面计算
C
b
0
e
0
C_{b0}^{e0}
Cb0e0
将e系中的重力矢量投影到e0l系:
g e 0 = C e e 0 g e g^{e0} = C_{e}^{e0} g^{e} ge0=Cee0ge
将b系中的比力投影到b0系:
f s f b 0 = C b b 0 f s f b f_{sf}^{b0} = C_{b}^{b0} f_{sf}^b fsfb0=Cbb0fsfb
容易得到:
C b 0 e 0 f s f b 0 = − g e 0 C_{b0}^{e0} f_{sf}^{b0} = -g^{e0} Cb0e0fsfb0=−ge0
即:
C b 0 e 0 C b b 0 f s f b = − C e e 0 g e C_{b0}^{e0} C_{b}^{b0} f_{sf}^b = -C_{e}^{e0} g^{e} Cb0e0Cbb0fsfb=−Cee0ge
实际上加计测量的比力总是含有误差,对准过程中多多少少会有载体线加速度误差,所以上式应为:
C b 0 e 0 ( C b b 0 f s f b − ∇ b 0 ) = − C e e 0 g e C_{b0}^{e0} (C_{b}^{b0} f_{sf}^b - \nabla^{b0} ) = -C_{e}^{e0} g^{e} Cb0e0(Cbb0fsfb−∇b0)=−Cee0ge
为了消除上述误差的影响,考虑对其进行积分,于是一次积分:
(1) C b 0 e 0 F b 0 = G e 0 \tag{1} C_{b0}^{e0} F^{b0} = G^{e0} Cb0e0Fb0=Ge0(1)
其中:
F
b
0
=
∫
0
Δ
t
(
C
b
b
0
f
s
f
b
−
∇
b
0
)
d
t
F^{b0} = \int_0^{\Delta t} (C_{b}^{b0} f_{sf}^b - \nabla^{b0}) dt
Fb0=∫0Δt(Cbb0fsfb−∇b0)dt,
G
e
0
=
−
∫
0
Δ
t
C
e
e
0
g
e
d
t
G^{e0} = -\int_0^{\Delta t} C_{e}^{e0} g^{e} dt
Ge0=−∫0ΔtCee0gedt
从开始到现在经过的时间为t,选取0.5t为另一矢量:
F
0.5
t
b
0
,
G
0.5
t
e
0
F^{b0}_{0.5t}, G^{e0}_{0.5t}
F0.5tb0,G0.5te0,得到基于速度的双矢量定姿:
C b 0 e 0 [ F b 0 ∣ F b 0 ∣ F b 0 × F 0.5 t b 0 ∣ F b 0 × F 0.5 t b 0 ∣ F b 0 × F 0.5 t b 0 × F b 0 ∣ F b 0 × F 0.5 t b 0 × F b 0 ∣ ] = [ G e 0 ∣ G e 0 ∣ G e 0 × G 0.5 t e 0 ∣ G e 0 × G 0.5 t e 0 ∣ G e 0 × G 0.5 t e 0 × G e 0 ∣ G e 0 × G 0.5 t e 0 × G e 0 ∣ ] C_{b0}^{e0} \begin{bmatrix} \frac{F^{b0}}{|F^{b0}|} & \frac{F^{b0}×F^{b0}_{0.5t}}{|F^{b0}×F^{b0}_{0.5t}|} & \frac{F^{b0}×F^{b0}_{0.5t}×F^{b0}}{|F^{b0}×F^{b0}_{0.5t}×F^{b0}|} \end{bmatrix} = \begin{bmatrix} \frac{G^{e0}}{|G^{e0}|} & \frac{G^{e0}×G^{e0}_{0.5t}}{|G^{e0}×G^{e0}_{0.5t}|} & \frac{G^{e0}×G^{e0}_{0.5t}×G^{e0}}{|G^{e0}×G^{e0}_{0.5t}×G^{e0}|} \end{bmatrix} Cb0e0[∣Fb0∣Fb0∣Fb0×F0.5tb0∣Fb0×F0.5tb0∣Fb0×F0.5tb0×Fb0∣Fb0×F0.5tb0×Fb0]=[∣Ge0∣Ge0∣Ge0×G0.5te0∣Ge0×G0.5te0∣Ge0×G0.5te0×Ge0∣Ge0×G0.5te0×Ge0]
得到:
C b 0 e 0 = [ G e 0 ∣ G e 0 ∣ G e 0 × G 0.5 t e 0 ∣ G e 0 × G 0.5 t e 0 ∣ G e 0 × G 0.5 t e 0 × G e 0 ∣ G e 0 × G 0.5 t e 0 × G e 0 ∣ ] [ F b 0 ∣ F b 0 ∣ F b 0 × F 0.5 t b 0 ∣ F b 0 × F 0.5 t b 0 ∣ F b 0 × F 0.5 t b 0 × F b 0 ∣ F b 0 × F 0.5 t b 0 × F b 0 ∣ ] − 1 C_{b0}^{e0} = \begin{bmatrix} \frac{G^{e0}}{|G^{e0}|} & \frac{G^{e0}×G^{e0}_{0.5t}}{|G^{e0}×G^{e0}_{0.5t}|} & \frac{G^{e0}×G^{e0}_{0.5t}×G^{e0}}{|G^{e0}×G^{e0}_{0.5t}×G^{e0}|} \end{bmatrix} \begin{bmatrix} \frac{F^{b0}}{|F^{b0}|} & \frac{F^{b0}×F^{b0}_{0.5t}}{|F^{b0}×F^{b0}_{0.5t}|} & \frac{F^{b0}×F^{b0}_{0.5t}×F^{b0}}{|F^{b0}×F^{b0}_{0.5t}×F^{b0}|} \end{bmatrix} ^{-1} Cb0e0=[∣Ge0∣Ge0∣Ge0×G0.5te0∣Ge0×G0.5te0∣Ge0×G0.5te0×Ge0∣Ge0×G0.5te0×Ge0][∣Fb0∣Fb0∣Fb0×F0.5tb0∣Fb0×F0.5tb0∣Fb0×F0.5tb0×Fb0∣Fb0×F0.5tb0×Fb0]−1
另外,还可以对公式(1)再进行一次积分,得到基于位置的双矢量定姿,在此不再过多介绍。可以看出,两次积分使得误差项 ∇ b 0 \nabla^{b0} ∇b0几乎可以忽略,即增强了抗线晃动的能力。
测试
另将对准期间的IMU数据作平均处理后,选取e系的重力与地球自转角速度作为双矢量,选取b系的比力与角速度作为双矢量,进行对准。载体处于准静止状态,即没有进行晃动模拟,对准时长为20分钟,对准结束后继续静止250s,进行惯导解算。本方案的对准结果如下:

基于惯性系的对准方案解算结果如下:

对比两幅图可以看出,基于惯性系的对准方案并没有取得很好的对准效果,甚至在X方向的漂移误差很大,所以对于准静态的载体而言,基于惯性系的对准方案并不合适,另外该方案计算量较大,降低了可用性。但是对于晃动的载体而言,基于惯性系的这种抗晃动的对准方案应该会具有很大的优势,但是没有实践过便没有发言权。
另外,这两种对准方法只是粗对准的方案,对于精度要求高的惯导,还是要进行精对准。看过一些文献,基于速度量测的卡尔曼滤波算法会有比较高的精度,后面会进行尝试。
参考:
[1]严恭敏,翁浚.捷联惯导算法与组合导航原理讲义.
[2]严恭敏,白亮,翁浚,秦永元.基于频域分离算子的SINS抗晃动干扰初始对准算法[J].宇航学报,2011,32(07):1486-1490.