CN103955005B - 一种火箭橇轨道重力实时测量方法 - Google Patents

Abstract

一种火箭橇轨道重力实时测量方法(1)在当前导航周期，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的火箭橇橇体高度信息h_B，进而得到高度误差信号δh＝h‑h_B；(2)将高度误差信号δh经过积分控制得到控制分量x，高度误差信号δh经过PID控制，得到高度控制量u_h；(3)将分量x经过低通滤波得到重力偏差x′，再将重力偏差x′与地球重力模型值g_h做和，得到实时的测量重力值g₀；(4)将u_h及g₀反馈至惯性导航系统，由惯性导航系统解算得到修正后的火箭橇橇体的高度h；(5)进入下一导航周期，惯性导航系统解算火箭橇橇体的纬度信息，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的h_B，重新得到高度误差信号δh；转步骤(2)循环执行。

Description

一种火箭橇轨道重力实时测量方法

技术领域

本发明涉及一种组合导航系统设计方法，尤其一种火箭橇轨道重力实时测量方法，可用于具有高精度导航精度的场合，例如航空航天、测绘等领域。

背景技术

在惯性导航系统中，高度通道的特征方程具有一个正根在高度通道的误差传播中具有与成比例的按时间随指数规律增长的分量。显然，如果天向通道存在初始高度误差Δh₀、初始天向速度误差、天向加速度计的标度因子误差以及常值误差，则高度通道随时间增长其误差是分散的。为了解决这一问题，需要在高度通道中引入外部信息进行阻尼。常用的外部高度信息包括雷达高度表、气压高度表以及大气数据中心等提供的数据，这些高度信息的误差不会发散。因此，在长航行时的精确导航中可以认为外部高度信息的精度应高于惯性导航系统提供的高度精度，利用高度计对系统做高度和天向速度补偿，这就是所谓的高度计/惯性组合导航系统或高度阻尼办法克服高度误差发散问题。引起纯惯性高度通道发散的根本原因是系统无阻尼，引入阻尼的两种常用途径为：(1)采用回路反馈法，此时系统阻尼系数是固定不变的；(2)采用卡尔曼滤波法，此时系统阻尼系数是时变得，根据外来高度参考信息误差特征实时确定。

另外，目前使用的重力测量方法是根据大地测量系统进行定点测量，测量时间长，测量点数多，且不能连续测量。因此，为了能够减小大地重力测量的难度和连续性，且为了能够实时的测量火箭橇轨道重力，需要提出一种火箭橇轨道重力实时测量方法。

发明内容

本发明的技术解决问题：克服现有技术的不足，提供一种火箭橇轨道重力实时测量方法，实现了对惯导高度通道的修正，不仅使得导航结果精确，且能够用于大地重力连续测量。

本发明的技术解决方案：一种火箭橇轨道重力实时测量方法，步骤如下：

(1)在当前导航周期，惯性导航系统解算火箭橇撬体的高度h以及火箭橇橇体的纬度信息，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的火箭橇橇体高度信息h_B，进而得到高度误差信号δh＝h-h_B；所述的惯性导航系统置于火箭橇橇体中；

(2)将高度误差信号δh经过积分控制得到控制分量x，高度误差信号δh经过PID控制，得到高度控制量u_h；

(3)将上述控制分量x经过低通滤波得到重力偏差x′，再将重力偏差x′与地球重力模型值g_h做和，得到实时的测量重力值g₀；

(4)将步骤(2)中得到的高度控制量u_h以及步骤(3)中得到实时的测量重力值g₀反馈至惯性导航系统，由惯性导航系统解算得到修正后的火箭橇撬体的高度h；

(5)进入下一导航周期，惯性导航系统解算火箭橇橇体的纬度信息，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的火箭橇橇体高度信息h_B，将步骤(4)中修正后的高度h与h_B做差，得到高度误差信号δh；转步骤(2)循环执行。

所述步骤(4)中的修正后的火箭橇撬体的高度h采用将下式积分后得到，

$\stackrel{..}{h}=a_z+a_{\mathrm{zB}}-g_0+\frac{2g_0}{R}h-u_h$

其中，表示高度信息h的二阶导数，a_z为加速度计敏感到的加速度在地理坐标系中的投影分量，v_e为火箭橇橇体东向速度，v_n为火箭橇橇体北向速度，ω_ie为地球自转角速度，为火箭橇橇体所在纬度，R为地球半径，M为地球子午圈曲率半径，N为地球子午圈曲率半径。

所述步骤(5)中的惯性导航系统解算火箭橇橇体的纬度信息通过下述步骤实现：

(5.1)设置火箭橇轨道坐标系(OX_lY_lZ_l)，该坐标系的原点为火箭橇轨道起始点，OX_l轴指向火箭橇橇体运动前进方向，OZ_l轴朝上垂直于轨道，OY_l轴在水平面内垂直于轨道，且三者满足右手坐标系；

(5.2)惯性测量系统进行自对准或者进行传递对准，得到火箭橇橇体的轨道坐标系下的三个姿态角初值；

(5.3)计算地球转速以及重力加速度在火箭橇轨道坐标系下的分量，结合惯性测量系统中捷联陀螺的输出，更新火箭橇橇体在火箭橇轨道坐标系下的三个姿态角；

(5.4)利用步骤(5.3)中更新后的姿态角计算火箭橇轨道坐标系到捷联本体坐标系的姿态变换矩阵；

(5.5)利用步骤(5.4)中的姿态变换矩阵、步骤(5.3)中的重力加速度在火箭橇轨道坐标系下的分量，结合惯性测量系统中捷联加表的输出，得到火箭橇橇体在火箭橇轨道坐标系的加速度，进而得到火箭橇橇体在火箭橇轨道坐标系的速度以及位置；

(5.6)将步骤(5.5)中得到的火箭橇橇体在火箭橇轨道坐标系下的位置信息转换为在地理坐标系下的火箭橇橇体位置信息。

本发明与现有技术相比的优点如下：

(1)现有的惯性导航算法高度通道不稳定，会随着时间增长误差不断增大。本发明能够使得高度通道稳定，能够实时地获得高度信息值，进行导航运算。

(2)现有的大地重力测量方法采用非连续的，定点式测量方法，无法连续实时的测量大地重力值，工作量大，且重复度高。

(3)与现有的大地重力测量方法相比，本发明不仅能够精确快速的测量出大地的重力值，而且能够连续的不间断的测量大地重力值，自动化水平高且人工参与度小，比现有的方法更进了一步。

(4)基于火箭撬发射点的轨道坐标系的导航算法中，对于直线轨道来说，基于火箭撬发射点的轨道坐标系OX轴与轨道平行，因此，可以直接得到橇体沿轨道运行的距离信息；对于与地球水平面平行的弯曲轨道，就可在基于火箭撬发射点的轨道坐标系OXZ平面内描述橇体沿轨道运行的距离信息。该计算方法简单，且物理意义明确。

(5)基于火箭撬发射点的轨道坐标系的导航算法姿态角更新时只需考虑陀螺仪组合的测量值即可，不需要考虑速度的影响，算法简单。而基于地理坐标系的导航算法在姿态更新时需要考虑速度的影响。

附图说明

图1为本发明的高度通道回路和重力测量流程图；

图2为本发明的高度通道回路和重力测量原理图；

图3为使用本发明方法设计的组合导航系统的高度通道输出结果图；

图4为使用本发明方法设计的重力偏差和重力测量曲线。

具体实施方式

由于惯性导航系统重力通道不稳定，惯性导航系统高度信息工作时导航误差随时间增长而增长，且真实重力值与重力模型值有误差。为此，需要引入外部高度信以对惯性导航高度通道进行修正，同时测量得到大地的真实重力值。

火箭橇轨道是一种特殊的跑车环境，根据理论设计严格施工，轨道的高度信息与纬度信息存在唯一的一一对应关系，这个关系在轨道设计中已经确定。所以在本专利中通过确定火箭橇橇体纬度得到火箭橇橇体实时的真实高度信息。

一种火箭橇轨道重力实时测量方法，如图1所示，具体步骤如下：

(1)在当前导航周期，惯性导航系统解算火箭橇撬体的高度h以及火箭橇橇体的纬度信息，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的火箭橇橇体高度信息h_B，进而得到高度误差信号δh＝h-h_B；所述的惯性导航系统置于火箭橇橇体中；

(2)将高度误差信号δh经过积分控制得到控制分量x，高度误差信号δh经过PID控制，得到高度控制量u_h；

如图2所示，将导航高度信息h和真实高度信息h_B做差，得到高度误差信号δh，高度误差信息经过积分控制器得到控制分量x，再将高度误差信号δh通过PID控制，得到高度控制量u_h。

PID控制器的形式为：

$C\left(s\right)=K_{1}s+K_2+\frac{K_3}{s}$

其中，K₁,K₂,K₃是PID控制器的三个系数，高度误差信号经过积分控制器得到控制分量x。

(3)将上述控制分量x经过低通滤波得到重力偏差x′，再将重力偏差x′与地球重力模型值g_h做和，得到实时的测量重力值g₀；

g₀＝g_h+x′

其中式中系数a₁,a₂,a₃,a₄,a₅,a₆的值取决于所用旋转椭球的参数值。

例：采用WGS84的重力模型值，PID控制器的三个控制量为 其中τ＝0.5，这里采用的低通滤波器为其中f＝0.001HZ。

(4)将步骤(2)中得到的高度控制量u_h以及步骤(3)中得到实时的测量重力值g₀反馈至惯性导航系统，由惯性导航系统解算得到修正后的火箭橇撬体的高度h；

修正后的火箭橇撬体的高度h采用将下式积分后得到，

$\stackrel{..}{h}=a_z+a_{\mathrm{zB}}-g_0+\frac{2g_0}{R}h-u_h$

其中，表示高度信息h的二阶导数，a_z为加速度计敏感到的加速度在地理坐标系中的投影分量，v_e为火箭橇橇体东向速度，v_n为火箭橇橇体北向速度，ω_ie为地球自转角速度，为火箭橇橇体所在纬度，R为地球半径，M为地球子午圈曲率半径，N为地球子午圈曲率半径。

(5)进入下一导航周期，惯性导航系统解算火箭橇橇体的纬度信息，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的火箭橇橇体高度信息h_B，将步骤(4)中修正后的高度h与h_B做差，得到高度误差信号δh；转步骤(2)循环执行。

上述惯性导航系统解算得到火箭橇橇体高度和纬度信息可以采用目前常规的导航解算方式，该方法在地理坐标系下进行，最终结果为直接获得火箭橇橇体的高度和纬度信息，此部分为本领域公知技术，不进行赘述。

因为火箭橇轨道的特殊性，还可以采用另一种导航解算方法，该方法在火箭橇轨道坐标系下进行，该坐标系的原点在火箭橇轨道起始点，该方法直接获得的为火箭橇橇体在火箭橇轨道坐标系下的三维坐标信息(x,y,z)，可以通过坐标变换将火箭橇橇体在火箭橇轨道坐标系l下的三维坐标信息(x,y,z)转换为在地理坐标系L下的火箭橇橇体位置信息转换方式如下所示：

其中，为火箭橇橇体在地理坐标系L下的位置信息，(x,y,z)火箭橇橇体在火箭橇轨道坐标系l下的三维坐标信息，为火箭橇轨道起始点在地理坐标系L下的位置信息，为从火箭橇轨道坐标系l到地理坐标系L的转换矩阵，可由轨道的设计图中直接获得。

下面介绍上面提及的新的导航解算方法，具体步骤如下：

(1)设置火箭橇轨道坐标系(OX_lY_lZ_l)，该坐标系的原点为火箭橇轨道起始点，OX_l轴指向火箭橇橇体运动前进方向，OZ_l轴朝上垂直于轨道，OY_l轴在水平面内垂直于轨道，且三者满足右手坐标系；

(2)惯性测量系统进行自对准或者进行传递对准，得到火箭橇橇体的轨道坐标系下的三个姿态角初值；

一般来讲，对于高精度惯性测量系统，可采用自对准方案，而对于精度相对较低的惯性测量系统，采用外部水平和方位传递对准的方式。

从地理坐标系到轨道坐标系的转换矩阵由大地测量得到，可以认为没有误差，地理坐标系原点为橇体中心，其中x、y、z三轴满足东北天准则。

$R_{L_0}^l=\left[\begin{array}{ccc}\cos {\mathrm{\β}}_0& 0& \sin {\mathrm{\β}}_0\\ 0& 1& 0\\ -\sin {\mathrm{\β}}_0& 0& \cos {\mathrm{\β}}_0\end{array}\right]\left[\begin{array}{ccc}\cos {\mathrm{\γ}}_0& \sin {\mathrm{\γ}}_0& 0\\ -\sin {\mathrm{\γ}}_0& \cos {\mathrm{\γ}}_0& 0\\ 0& 0& 1\end{array}\right]$

其中，β₀和γ₀分别为轨道相对发射点地理坐标系的俯仰角和偏航角，其示意图如图3所示。

在初始对准时，保持惯组静止一段时间，记录这段时间内的在橇体坐标系三个轴向测量到的平均加速度为平均角速度为因此可以得到橇体坐标系b和地理坐标系L₀之间的姿态角为：

$\left\{\begin{array}{c}{\mathrm{\φ}}_x^{b,L_0}=\arcsin \frac{A_y^b}{g_0}\\ {\mathrm{\φ}}_y^{b,L_0}=-\arcsin \frac{A_x^b}{A_z^b}\\ {\mathrm{\φ}}_z^{b,L_0}=\arcsin \frac{G_x^b\cos {\mathrm{\φ}}_y^{b,L}+G_z^b\sin {\mathrm{\φ}}_y^{b,L}}{{\mathrm{\ω}}_{\mathrm{ie}}\cos L}\end{array}\right.$

其中，g₀为测试点重力加速度，ω_ie为地球自转角速度，L为测试点纬度。橇体坐标系b与橇体固连，原点为橇体中心，x轴指向运动方向，z轴指天，y轴分别与x、z轴垂直，且满足右手准则。

地理坐标系L₀到橇体坐标系b的姿态变换矩阵为

$R_{L_0}^b=\left[\begin{array}{ccc}\cos {\mathrm{\φ}}_y^{b,L_0}& 0& -\sin {\mathrm{\φ}}_y^{b,L_0}\\ 0& 1& 0\\ -\sin {\mathrm{\φ}}_y^{b,L_0}& 0& \cos {\mathrm{\φ}}_y^{b,L_0}\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\mathrm{\φ}}_x^{b,L_0}& \sin {\mathrm{\φ}}_x^{b,L_0}\\ 0& -{\sin \mathrm{\φ}}_x^{b,L_0}& \cos {\mathrm{\φ}}_x^{b,L_0}\end{array}\right]\left[\begin{array}{ccc}\cos {\mathrm{\φ}}_z^{b,L_0}& \sin {\mathrm{\φ}}_z^{b{,L}_0}& 0\\ -\sin {\mathrm{\φ}}_z^{b,L_0}& \cos {\mathrm{\φ}}_z^{b,L_0}& 0\\ 0& 0& 1\end{array}\right]$

且

$R_b^{L_0}={\left(R_{L_0}^b\right)}^T$

由此可得从橇体坐标系到轨道坐标系的姿态变换矩阵初值以及 $R_l^b={\left(R_b^l\right)}^{-1}.$ 因此，有

$\left\{\begin{array}{c}{\mathrm{\φ}}_x=\arcsin \left[R_l^b\left(\mathrm{2,3}\right)\right]\\ {\mathrm{\φ}}_y=-\mathrm{arctam}\left[\frac{R_l^b\left(\mathrm{1,3}\right)}{R_l^b\left(\mathrm{3,3}\right)}\right]\\ {\mathrm{\φ}}_z=-\arctan \left[\frac{R_l^b\left(\mathrm{2,1}\right)}{R_l^b\left(\mathrm{2,2}\right)}\right]\end{array}\right.$

(3)计算地球转速以及重力加速度在火箭橇轨道坐标系下的分量，结合惯性测量系统中捷联陀螺的输出，更新火箭橇橇体在火箭橇轨道坐标系下的三个姿态角；

由于轨道坐标系与地球固联，则地球转速在轨道坐标系的分量为常值向量

其中，为发射点的地理纬度，ω_ie为地球自转速度。

在捷联陀螺测量得到后，就可利用下面微分方程对三个姿态角进行更新：

其中，为第n次更新时的时间tn对应的捷联陀螺的输出；为地球转速在轨道坐标系中的投影分量；

将第n次更新时的时间t_n对应的三个姿态角根据上述积分即可得到下一更新时间t_n+1对应的姿态角其中，当n＝0时，三个姿态角的取值为步骤(2)的结果。

(4)利用步骤(3)中更新后的姿态角计算火箭橇轨道坐标系到捷联本体坐标系的姿态变换矩阵；

(5)利用步骤(4)中的姿态变换矩阵、步骤(3)中的重力加速度在火箭橇轨道坐标系下的分量，结合惯性测量系统中捷联加表的输出，得到火箭橇橇体在火箭橇轨道的加速度，进而得到火箭橇橇体在火箭橇轨道的速度以及位置。

发射点轨道坐标系的重力加速度分量为

$\left[\begin{array}{c}{g}_x^l\\ g_y^l\\ g_z^l\end{array}\right]=\left[\begin{array}{c}g\left(r_x\right)\sin \mathrm{\β}\left(r_x\right)\\ 0\\ g\left(r_x\right)\cos \mathrm{\β}\left(r_x\right)\end{array}\right]=\left[\begin{array}{c}-(g_0+b_1{r}_x+b_2{r}_x^2)\sin \frac{\left|P_{0}P\right|-r_x}{N+h_p}\\ 0\\ (g_0+b_1{r}_x+b_2{r}_x^2)\cos \frac{\left|P_{0}P\right|-r_x}{N+h_p}\end{array}\right]$

其中，地球长半轴a＝6378137±2m，地球偏心率e²，为载体所在点的纬度值；h_p为橇体相对水准面的高度；|P₀P|＝(N+h_p)β′，β′为矢量OP₀和OP的夹角，其中矢量OP₀为地球中心O到轨道初始点P₀的矢量，OP为地球中心O到P点的矢量，P为轨道与地球表面切点，为火箭橇轨道重力加速度模型，其中r_x为橇体沿轨道运行的距离；g₀为发射点位置的重力加速度，b₁和b₂为常值，由轨道重力测量值拟合确定，具体拟合即采用现有方法每隔固定距离(例如60m)测量重力加速度，将所有测量值采用最小二乘法进行拟合后即可得到b₁和b₂值。

惯性导航系统中的捷联加表测量值为f^b(t_n)，就可利用下述速度微分方程

${\stackrel{.}{V}}^l\left(t_n\right)=R_b^l\left(t_{n+1}\right)f^b\left(t_n\right)-2{\mathrm{\Ω}}_{\mathrm{ie}}^l{V}^l\left(t_n\right)+g^l$

对速度进行更新，就可得到一个新的在轨道坐标系下的速度矢量V^l(t_n+1)。其中g^l为重力矢量在轨道坐标系中的投影，为地球转速在轨道坐标系投影分量的反对称矩阵，

在得到速度矢量V^l(t_n+1)后，就可利用位置微分方程

${\stackrel{.}{r}}^l\left(t_n\right)=V^l\left(t_{n+1}\right)$

对位置进行更新，就可得到一个新的在轨道坐标系下的位置矢量r^l(t_n+1)即(x,y,z)；在导航时间内，重复步骤(3)-(5)，完成惯性测量系统基于火箭橇轨道坐标系的导航解算。

附图3为惯性导航系统导航解算结果中的高度与根据纬度获得的轨道真实高度进行比较，从图中可以看出导航结果高度与轨道真实高度吻合。附图4为重力偏差(本发明方法)和重力偏差与理论重力值之和即为本发明测得的实时重力值，可以看出理论重力值与测得的实时重力值之间存在一定的重力偏差，可以通过本方明将其测量出来，得到真实的重力值。将测量得到的真实重力值反馈到导航解算中，得到更精确的导航结果，获得更高的导航定位精度，进而可以完成基于外部高度信息的惯性导航系统高度通道修正方法组合导航系统设计，并完成了火箭橇轨道重力的实时测量。

本发明未详细描述内容为本领域技术人员公知技术。

Claims (3)

1.一种火箭橇轨道重力实时测量方法，其特征在于步骤如下：

(1)在当前导航周期，惯性导航系统解算火箭橇橇体的高度h以及火箭橇橇体的纬度信息，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的火箭橇橇体高度信息h_B，进而得到高度误差信号δh＝h-h_B；所述的惯性导航系统置于火箭橇橇体中；

(2)将高度误差信号δh经过积分控制得到控制分量x，高度误差信号δh经过PID控制，得到高度控制量u_h；

(3)将上述控制分量x经过低通滤波得到重力偏差x′，再将重力偏差x′与地球重力模型值g_h做和，得到实时的测量重力值g₀；

(4)将步骤(2)中得到的高度控制量u_h以及步骤(3)中得到实时的测量重力值g₀反馈至惯性导航系统，由惯性导航系统解算得到修正后的火箭橇橇体的高度h；

(5)进入下一导航周期，惯性导航系统解算火箭橇橇体的纬度信息，根据火箭橇轨道中的纬度与高度信息的关系，得到当前导航周期下真实的火箭橇橇体高度信息h_B，将步骤(4)中修正后的高度h与h_B做差，得到高度误差信号δh；转步骤(2)循环执行。

2.根据权利要求1所述的一种火箭橇轨道重力实时测量方法，其特征在于：所述步骤(4)中的修正后的火箭橇橇体的高度h采用将下式积分后得到，

$\stackrel{\·\·}{h}=a_z+a_{zB}-g_0+\frac{2g_0}{R}h-u_h$

其中，表示高度信息h的二阶导数，a_z为加速度计敏感到的加速度在地理坐标系中的投影分量，v_e为火箭橇橇体东向速度，v_n为火箭橇橇体北向速度，ω_ie为地球自转角速度，为火箭橇橇体所在纬度，R为地球半径，M为地球子午圈曲率半径，N为地球子午圈曲率半径。

3.根据权利要求1所述的一种火箭橇轨道重力实时测量方法，其特征在于：所述步骤(5)中的惯性导航系统解算火箭橇橇体的纬度信息通过下述步骤实现：

(5.1)设置火箭橇轨道坐标系OX_lY_lZ_l，该坐标系的原点为火箭橇轨道起始点，OX_l轴指向火箭橇橇体运动前进方向，OZ_l轴朝上垂直于轨道，OY_l轴在水平面内垂直于轨道，且三者满足右手坐标系；

(5.2)惯性测量系统进行自对准或者进行传递对准，得到火箭橇橇体的轨道坐标系下的三个姿态角初值；

(5.3)计算地球转速以及重力加速度在火箭橇轨道坐标系下的分量，结合惯性测量系统中捷联陀螺的输出，更新火箭橇橇体在火箭橇轨道坐标系下的三个姿态角；

(5.4)利用步骤(5.3)中更新后的姿态角计算火箭橇轨道坐标系到捷联本体坐标系的姿态变换矩阵；

(5.5)利用步骤(5.4)中的姿态变换矩阵、步骤(5.3)中的重力加速度在火箭橇轨道坐标系下的分量，结合惯性测量系统中捷联加表的输出，得到火箭橇橇体在火箭橇轨道坐标系的加速度，进而得到火箭橇橇体在火箭橇轨道坐标系的速度以及位置；

(5.6)将步骤(5.5)中得到的火箭橇橇体在火箭橇轨道坐标系下的位置信息转换为在地理坐标系下的火箭橇橇体位置信息。