微惯性传感器高过载对准方法

High-Overload Alignment Method Based on MEMS

摘要: 针对高速、过载条件下的微惯性传感器对准问题，本文设计了一种加速度计冗余对称安装的配置方案，提出了一种基于卫导信息 + 自旋信号分解的对准方法。根据卫导信息建立姿态推算算法快速解算弹体航向和俯仰初始角度。针对弹体自旋情况，利用对称双加速度计配置方案，根据双加速度计量测信号的相位辨识及误差补偿技术，结合微惯性陀螺转速信号，运用递推最小二乘算法，突破基于周期信号的横滚姿态快速获取技术，实现制导炮弹空中高动态快速对准。通过仿真分析在弹体高速旋转条件下，水平对准误差小于0.2˚，航向误差小于0.5˚的精度水平。

Abstract: Aiming at the alignment problem of Micro-Electro-Mechanical System (MEMS) under high speed and overload conditions, this paper designed a configuration scheme of redundant and symmetrical installation of accelerometer, and proposed an alignment method based on the satellite navigation information + spin signal decomposition. Based on the satellite navigation information, the attitude algorithm is established to quickly solved the projectile heading and inertial pitching angle. In view of the spin condition of the projectile body, the phase identification and error compensation technology of measured signal are used according to the configuration scheme of symmetrical double accelerometer. Combined with gyro signal and recursive least square algorithm, this method breaks through the technology of rapid acquisition of roll attitude based on periodic signal, and realizes high dynamic and rapid aerial alignment of guided projectiles. Through simulation analysis, the horizontal alignment error is less than 0.2˚ and the course error is less than 0.5˚ under the condition of high speed rotation of projectile body.

1. 引言

微机电惯性传感器在高过载环境中面临两大问题，一是高过载过程不易存活，即陀螺在过载作用后失效，丧失陀螺效应；二是高过载后陀螺性能退化，即过载前后陀螺相关参数发生较大变化，其中以陀螺零偏的变化表现明显。因此抗高过载技术研究，已成为微惯性传感器领域的重点研究方向，美国高级研究计划局(DARPA)更是把抗高过载陀螺仪作为关键器件之一列在了“μPNT”计划中。

随着微机电惯性传感器工程化程度不断提高，对相关元件级、系统级加固和过载防护技术研究的深入开展，使微机电惯性传感器能够承受的高过载量级逐步提升，进一步促进了以微机电惯性测量单元为基础的弹药制导技术得到长足发展，形成一条低成本、高效益的弹药精确制导化发展道路 [1] [2]。

美国Honeywell公司研制的陀螺结构，20,000 g过载后存活，但零偏稳定性从18˚/h下降到180˚/h，性能退化近10倍；法国THALES公司研制的微机电陀螺仪，20,000 g过载后虽然存活，但多项参数均明显恶化，如：零偏稳定性退化587%，标度因数稳定性退化218%。

在基于制导炮弹的MEMS惯性导航系统的对准中，需要对制导炮弹自旋角度进行的较为准确的分解，但传统的三轴正交加速度计在信号分解求解横滚角度时，受到弹体高速旋转和飞行运动的影响，旋转炮弹特有的运动特征在飞行过程中由于轴向存在自旋角速率而产生锥形运动。弹体将绕速度轴作周期性划圆运动，表现为进动和章动的耦合运动。一般情况下，由于制导炮弹处于高旋和高马赫飞行状态，弹体也会产生锥形运动。此时的偏航角速率和俯仰角速率呈现不规则的正弦振荡形式，这将导致测量出的轴向角速率值并非真值，耦合了锥形运动引起的附加信息 [3]。

当采用弹载加速度计测试法，为了避免炮弹在高动态飞行过程中存在进动和章动，其法向加速度会受到空气升力的影响发生变化，影响转速测量精度，在自旋轴敏感方向配置了对称加速度计，用以消除弹体高速自旋和线运动的影响，达到准确估计横滚角的目的 [4] [5]。

MEMS惯性组件的配置方案本质上是在空间纬度上的数学平台的构建，传统的微惯性组件在结构设计时如图1所示，把加速度计的敏感点设计在载体的旋转轴上，避免加速度计由于角运动带来的尺寸效应。但是由于制导炮弹内部空间的约束条件，使得物理上难以保证，因此加速度计根据弹体实际情况分布在弹体的不同位置，加速度计敏感点不在载体的旋转中心上，将引起尺寸效应，导致加速度计的输出存在向心加速度和切向加速度带来的加速度误差干扰。

Figure 1. The scheme of MEMS configuration

图1. MEMS惯组配置方案

2. 卫导/微惯性组合对准应用

针对弹体自旋情况，利用对称双加速度计配置方案，运用双加速度计量测信号的相位辨识及误差补偿技术，结合微惯性陀螺转速信号，获取精确的横滚姿态。假设弹体出炮口后只受空气动力的影响而无其它干扰，弹体旋转角速度也趋于平稳。弹体的攻角较小、并在垂直平面内飞行，则无控段时弹体的俯仰角和偏航角可以用弹道倾角和弹道偏角替代。利用卫导测得的速度信息分量可以计算弹体的弹道倾角和弹道偏角 [3]，计算公式(1)如下：

$\begin{array}{l}{\theta }_v=\arctan \left(v_y/\sqrt{v_x^2+v_z^2}\right)\\ {\Psi }_v=\arctan \left(-v_z/v_x\right)\end{array}$ (1)

弹体在无控阶段、在气动力的作用下弹道基本稳定，无控阶段的攻角和侧滑角都很小，姿态角都在弹道角附近摆动，因此可以采取以下近似：

$\begin{array}{l}\theta ={\theta }_v\\ \stackrel{\dot{}}{\theta }={\stackrel{\dot{}}{\theta }}_v\\ \stackrel{\dot{}}{\Psi }={\stackrel{\dot{}}{\Psi }}_v\end{array}$ (2)

利用卫导得到 $\stackrel{\dot{}}{\Psi }$ 和 $\stackrel{\dot{}}{\theta }$，等式右边的数值都可以求得。

弹体绕质心转动的运动学方程(3)。

$\begin{array}{l}\frac{\text{d}\theta }{\text{d}t}={\omega }_x\sin \gamma +{\omega }_z\cos \gamma \\ \frac{\text{d}\Psi }{\text{d}t}=\left({\omega }_y\cos \gamma -{\omega }_z\sin \gamma \right)/\cos \theta \\ \frac{\text{d}\gamma }{\text{d}t}={\omega }_x-\tan \theta \left({\omega }_y\cos \gamma -{\omega }_z\sin \gamma \right)\end{array}$ (3)

(3)式中， ${\omega }_x,{\omega }_y,{\omega }_z$ 为弹体转动角度 $\omega$ 在弹体坐标系三个轴上的分量，由于角速率 $\omega$ 可由三轴速率陀螺测量得到，因此为了求解横滚角速率可以将上式变形为 [6]：

$\begin{array}{l}{\omega }_y=\stackrel{\dot{}}{\theta }\sin \gamma +\stackrel{\dot{}}{\Psi }\cos \theta \cos \gamma \\ {\omega }_z=-\stackrel{\dot{}}{\Psi }\cos \theta \sin \gamma +\stackrel{\dot{}}{\theta }\cos \gamma \end{array}$ (4)

在卫导信息得到 $\stackrel{\dot{}}{\Psi }$ 和 $\stackrel{\dot{}}{\theta }$ 情况下，可以求得滚转角 $\gamma$。

3. 对称加速度计自旋信号分解

虽然在2.2中，卫导与微惯性的组合对准中可以获得滚转角，但是该方法对于弹体的运行轨迹要求理想化，当考虑气动干扰时该方法受弹体运行限制条件，得到滚转角精度低。因此，本文采用对称加速度计自旋信号分解的方式来进一步更为准确的横滚角。

对称加速度计脉冲输出可表示如下：

$\{\begin{array}{l}{A}_{x1}=a\times \sin \left(\omega t+{\gamma }_0\right)\\ A_{x2}=a^{\prime }\times \sin \left(\omega t+{\gamma }_0+\text{π}\right)\end{array}$ (5)

其中， ${\gamma }_0$ 为弹体初始状态的滚转角。

在对称加速度计处理时频信息提取初始滚转角 ${\gamma }_0$，采用基于加权递推最小二乘法来进行求解，根据递推回归原理，通过加权处理环节，实现对处于不同空间位置的双加速度计的观测量权值调节，使得即使不在同一基准面配置的加速度计能够准确估计滚转初始角 [7]。

最小二乘法的测量方程如下：

$\{\begin{array}{l}{f}_1\left(x_1,x_2,x_3,\cdots ,x_k\right)=y_1-e_1\hfill \\ f_2\left(x_1,x_2,x_3,\cdots ,x_k\right)=y_2-e_2\hfill \\ ⋮\hfill \\ f_n\left(x_1,x_2,x_3,\cdots ,x_k\right)=y_n-e_n\hfill \end{array}$ (6)

其中， $x_1,x_2,\cdots ,x_k$ 是要估计的“ $k$ ”个状态； $y_1,y_2,\cdots ,y_n$ 是 $n$ 次观测； $e_1,e_2,\cdots ,e_n$ 是相应的 $n$ 次观测误差。要实现这些状态“最佳”估计，则选用 $x_1,x_2,\cdots ,x_k$ 的值，使得 $Q_n={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[f_i\left(x_1,x_2,\cdots ,x_k\right)-y_i\right]}^2}$ 最小。

假设公式(6)中函数关系 $f_i\left(x_1,x_2,\cdots ,x_k\right)-y_i\left(i=1,2,\cdots ,n\right)$ 的偏导数是连续的，利用雅克比线性化法进行线性化，得：

$\{\begin{array}{l}{f}_1\left(x_1,x_2,x_3,\cdots ,x_k\right)=b_{11}{x}_1+b_{12}{x}_2+\cdots +b_{1k}{x}_k\hfill \\ f_2\left(x_1,x_2,x_3,\cdots ,x_k\right)=b_{21}{x}_1+b_{22}{x}_2+\cdots +b_{2k}{x}_k\hfill \\ ⋮\hfill \\ f_n\left(x_1,x_2,x_3,\cdots ,x_k\right)=b_{n1}{x}_1+b_{n2}{x}_2+\cdots +b_{nk}{x}_k\hfill \end{array}$ (7)

由(6)、(7)可得：

$\{\begin{array}{l}{b}_{11}{x}_1+b_{12}{x}_2+\cdots +b_{1k}{x}_k=y_1-e_1\hfill \\ b_{21}{x}_1+b_{22}{x}_2+\cdots +b_{2k}{x}_k=y_2-e_2\hfill \\ ⋮\hfill \\ b_{n1}{x}_1+b_{n2}{x}_2+\cdots +b_{nk}{x}_k=y_n-e_n\hfill \end{array}$ (8)

其中 $n\ge k$。

式(8)改为矩阵形式为：

$B_{n}X=Y_n-{\epsilon }_n$ (9)

其中： $B_n=\left(\begin{array}{ccc}{b}_{11}& \cdots & b_{1k}\\ ⋮& \ddots & ⋮\\ b_{n1}& \cdots & b_{nk}\end{array}\right)$， $X=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_k\end{array}\right]$， $Y=\left[\begin{array}{c}{Y}_1\\ ⋮\\ Y_n\end{array}\right]$， ${\epsilon }_n=\left[\begin{array}{c}{e}_1\\ ⋮\\ e_n\end{array}\right]$。

令：

$d{Q}_n={\left(B_{n}X-Y_n\right)}^{\text{T}}{B}_{n}dX+d{X}^{\text{T}}{B}_n^{\text{T}}\left(B_{n}X-Y_n\right)=0$ (10)

假设再进行一次观测，可得：

$M_{n+1}X=y_{n+1}-e_{n+1}$ (11)

状态量 $X$ 的估计值 ${\hat{X}}_{n+1}$ 表达式：

$\left(B_{n+1}^{\text{T}}{B}_{n+1}\right){\hat{X}}_{n+1}=B_{n+1}^{\text{T}}{Y}_{n+1}$ (12)

令 ${\hat{X}}_{n+1}={\hat{X}}_n+\left({\hat{X}}_{n+1}-{\hat{X}}_n\right)$，利用式(12)可得：

$\left[B_n^{\text{T}}{B}_n+M_{n+1}^{\text{T}}{M}_{n+1}\right]\left({\hat{X}}_{n+1}-{\hat{X}}_n\right)=M_{n+1}^{\text{T}}\left(y_{n+1}-M_{n+1}{\hat{X}}_n\right)$ (13)

为避免在引入新观测量时出现 $k\times k$ 的矩阵求逆，应用Sherman-Morisson和Bartlett公式，令 $P_n={\left(B_n^{\text{T}}{B}_n\right)}^{-1}$， $P_{n+1}={\left(B_{n+1}^{\text{T}}{B}_{n+1}\right)}^{-1}={\left(B_n^{\text{T}}{B}_n+M_{n+1}^{\text{T}}{M}_{n+1}\right)}^{-1}$ 可得：

$\begin{array}{l}{P}_{n+1}=P_n-P_n{M}_{n+1}^{\text{T}}{\left(1+M_{n+1}{P}_n{M}_{n+1}^{\text{T}}\right)}^{-1}{M}_{n+1}{P}_n\hfill \\ {\hat{X}}_{n+1}={\hat{X}}_n+P_{n+1}{M}_{n+1}^{\text{T}}\left(y_{n+1}-M_{n+1}{\hat{X}}_n\right)\hfill \end{array}$ (14)

对测量信息进行加权处理，使噪声期望值 $Q_n$ 最小，避免星光/惯性组合导航系统在动态干扰条件下量测噪声会带来滤波震荡 [8] [9]。

噪声期望值 $Q_n$ 表达式：

$Q_n={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[\frac{f_i\left(x_1,x_2,\cdots ,x_k\right)-y_i}{{\sigma }_i}\right]}^2}$ (15)

由式(15)可得：

$Q_n={\left(B_{n}X-Y_n\right)}^{\text{T}}{\Sigma }_n^{-1}\left(B_{n}X-Y_n\right)$ (16)

其中 ${\Sigma }_n=\left(\begin{array}{ccc}{\sigma }_1^2& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\sigma }_n^2\end{array}\right)$。

则 ${\hat{X}}_n={\left(B_n^{\text{T}}{\Sigma }_n^{-1}{B}_n\right)}^{-1}{B}_n^{\text{T}}{\Sigma }_n^{-1}{Y}_n$，进行加权可得额外的信息时，则有：

$\left[B_n^{\text{T}}{\Sigma }_n^{-1}{B}_n+\frac{M_{n+1}^{\text{T}}{M}_{n+1}}{{\sigma }_{n+1}^2}\right]{\hat{X}}_{n+1}=B_n^{\text{T}}{\Sigma }_n^{-1}{Y}_n+\frac{M_{n+1}^{\text{T}}{Y}_{n+1}}{{\sigma }_{n+1}^2}$ (17)

利用式(17)得到：

$\{\begin{array}{l}{P}_{n+1}=P_n-P_n{M}_{n+1}^{\text{T}}{\left(1+\frac{M_{n+1}{P}_n{M}_{n+1}^{\text{T}}}{{\sigma }_{n+1}^2}\right)}^{-1}\frac{1}{{\sigma }_{n+1}^2}{M}_{n+1}{P}_n\\ {\hat{X}}_{n+1}={\hat{X}}_n+P_{n+1}{M}_{n+1}^{\text{T}}\left(\frac{y_{n+1}-M_{n+1}{\hat{X}}_n}{{\sigma }_{n+1}^2}\right)\end{array}$ (18)

其中 $P_n={\left(B_n^{\text{T}}{\Sigma }_n^{-1}{B}_n\right)}^{-1}$， $P_n={\left(B_n^{\text{T}}{\Sigma }_n^{-1}{B}_n+\frac{M_{n+1}^{\text{T}}{M}_{n+1}}{{\sigma }_{n+1}^2}\right)}^{-1}$。

在高动态条件下，卫导信息快速解算弹体的俯仰和偏航角度时，其精度受卫导的采样频率影响明显，在仿真中卫导以10 Hz采样，俯仰角32˚，偏航角2˚约在5秒内可以收敛到0.2˚的误差范围内，仿真结果如图2所示。

在双加速度计的弹道滚转量测仿真，预设弹丸从膛内发射后飞行稳定性，对称配置的加速度计抵消了绝大部分飞行过程中线运动引入的干扰因素，在弹丸转速满足高速旋转10 r/s的动态条件时，采用递推最小二乘可以获取滚转角 ${\gamma }_0$ 初始值。在高速旋转条件下， ${\gamma }_0$ 的精度可以在0.5˚以内，仿真结果如图3所示。

Figure 2. The error of pitch and yaw alignment based on satellite navigation information

图2. 基于卫导信息的俯仰及偏航对准误差

Figure 3. Initial error of roll angle based on recursive square spin decomposition

图3. 基于递推二乘自旋分解的横滚角初始误差

4. 结论

基于卫导信息 + 自旋信号分解的对准技术。根据卫导信息，建立姿态推算算法快速解算弹体航向和俯仰初始角度；针对弹体自旋情况，利用对称双加速度计配置方案，根据双加速度计量测信号的相位辨识及误差补偿技术，结合MEMS陀螺转速(角频率)信号，运用递推最小二乘算法，突破基于周期信号的横滚姿态快速获取技术，解决了弹体横滚角的精度测量问题，实现制导炮弹空中高动态快速对准。通过仿真分析在弹体高速旋转条件下，水平对准误差小于0.2˚，航向误差小于0.5˚的精度水平。