一种捷联式航空重力仪误差分离方法

摘要

本发明公开了一种捷联式航空重力仪误差分离方法，其步骤为：（1）、捷联式航空重力仪原始数据的预处理：将捷联式航空重力仪得到的加速度计和陀螺的原始脉冲信号通过当量转换成正常的测量信息；（2）、确定捷联式航空重力仪在飞行测量过程中产生的姿态误差；（3）、通过步骤（2）确定的姿态误差对捷联式重力仪的测量结果进行修正，实现误差的分离。本发明具有原理简单、操作简便、分离精度高、易推广使用等优点。

说明书

技术领域

本发明主要涉及到捷联式航空重力仪领域，特指一种捷联式航空重力仪误差分离方法。

背景技术

在捷联式航空重力仪中加速度计和陀螺是关键器件。由于加速度计和陀螺存在偏差和漂 移，加速度计和陀螺的精度就受到了限制，尤其是陀螺的精度会影响到最终估计出来的重力 仪的姿态，使得估计出来的姿态中存在较大的姿态误差，这样就会在捷联式航空重力仪最终 的测量结果中造成测量误差。

目前，对于姿态误差的分离并没有一个完全成熟的办法，国外有学者提出了小波降噪、 波数相关滤波、比例尺校正等方法进行误差分离。这些方法都是事后在数值处理上采用的一 些数据处理方面的技巧，并没有从根本上对姿态误差进行估计，所能获得的效果也不是很理 想。所以，要想解决捷联式航空重力仪的误差分离问题就必须设计一种全新的姿态误差的估 计方法。

发明内容

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种原理 简单、操作简便、分离精度高、易推广使用的捷联式航空重力仪误差分离方法。

为解决上述技术问题，本发明采用以下技术方案：

一种捷联式航空重力仪误差分离方法，其步骤为：

（1）、捷联式航空重力仪原始数据的预处理：将捷联式航空重力仪得到的加速度计和陀 螺的原始脉冲信号通过当量转换成正常的测量信息；

（2）、确定捷联式航空重力仪在飞行测量过程中产生的姿态误差；

（3）、通过步骤（2）确定的姿态误差对捷联式重力仪的测量结果进行修正，实现误差的 分离。

作为本发明的进一步改进：所述步骤（1）为：根据捷联式航空重力仪标定得到的标定参 数对重力仪得到的原始脉冲数据进行当量转换，以获取正常的测量信息；

f=K_a·N_a+B_a   (1)

w=K_g·(N_g-B_g)    (2)

其中，K_a和B_a分别表示加速度计的刻度因子和零偏；K_g和B_g分别表示陀螺的刻度因 子和零偏；N_a和N_g分别表示加速度计和陀螺的原始脉冲数据；f和w分别表示当量转换后 的测量信息。

作为本发明的进一步改进：所述步骤（2）中，通过步骤（1）获得捷联式航空重力仪的 正常测量信息之后，分以下三步确定捷联式航空重力仪的姿态误差：

（2.1）、在当地地理坐标系中捷联式航空重力仪的误差动态方程为：

$\delta \stackrel{.}{p}=\mathrm{\delta v}---\left(3\right)$

$\delta \stackrel{.}{v}=[f^n\times ]\psi +C_b^n\delta f^b-({2\omega }_{\mathrm{ie}}^n+{\omega }_{\mathrm{en}}^n)\times \mathrm{\delta v}$

     (4)

$-(2\delta {\omega }_{\mathrm{ie}}^n+\delta {\omega }_{\mathrm{en}}^n)\times v-\delta g^n$

$\stackrel{.}{\psi }=-{\omega }_{\mathrm{in}}^n\times \psi +\delta {\omega }_{\mathrm{in}}^n-C_b^n\delta {\omega }_{\mathrm{ib}}^b---\left(5\right)$

其中，δp、δv和ψ分别表示捷联式航空重力仪的位置误差、速度误差、姿态误差；fⁿ表示在当地地理坐标系下的加速度计比力测量值；表示从坐标系到当地地理坐标系的转移 矩阵；δf^b和分别表示加速度计和陀螺的测量误差；表示地球坐标系相对于惯性坐标 系的转动速度在当地地理坐标系中的投影；表示当地地理坐标系相对于地球坐标系的转 动速度在当地地理坐标系下的投影；δgⁿ表示当地地理坐标系下的重力扰动；

选取位置误差、速度误差和姿态误差为状态变量，可以得到卡尔曼滤波的状态方程：

$\frac{d}{\mathrm{dt}}\left(\begin{array}{c}\mathrm{\delta p}\\ \mathrm{\delta v}\\ \psi \end{array}\right)=F\left(\begin{array}{c}\mathrm{\delta p}\\ \mathrm{\delta v}\\ \psi \end{array}\right)+\mathrm{Gu}---\left(6\right)$

（2.2）、利用卡尔曼滤波方法进行状态变量X的前向估计，一个完成的卡尔曼滤波过程 可以表示如下：

X_k+1/k=φ_kX_k/k   (10)

$P_{k+1/k}={\phi }_k{P}_{k/k}{\phi }_k^T+Q_k---\left(11\right)$

X_k+1/k+1=X_k+1/k-K_k+1(H_k+1X_k+1/k-z_k+1)   (12)

$P_{k+1/k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{(I-K_{k+1}{H}_{k+1})}^T-K_{k+1}{\mathrm{PK}}_{k+1}^T---\left(13\right)$

其中，

$K_{k+1}=P_{k+1/k}{H}_{k+1}^T(H_{k+1}{P}_{k+1/k}{H}_{k+1}^T+R)---\left(14\right)$

（2.3）、利用一个周期内前向估计得到的状态变量进行反向恢复得到该周期起始点处的 状态变量的最终估计值，其中的姿态误差就是最终估计出来的姿态误差；反向恢复的计算 如下：

X_0/0=φ_0,NX_N/N   (15)

其中，X_N/N表示一个前向估计—反向恢复计算周期终点处的状态估计值；X_0/0表示一 个前向估计—反向恢复计算周期起点处的状态估计值；φ_0,N表示一个前向估计—反向恢复计 算周期终点到起点的转移矩阵。

作为本发明的进一步改进：所述步骤（3）通过下式来进行误差的分离：

$d\tilde{g}=\mathrm{dg}+\delta f^n---\left(16\right)$

$\delta f^n=f^n\times \tilde{\psi }---\left(17\right)$

其中，表示误差分离后得到的重力扰动；dg表示误差分离前重力扰动测量结果；δfⁿ表示由于姿态误差引起的测量误差的大小。

与现有技术相比，本发明的优点在于：

(1)本发明的捷联式航空重力仪误差分离方法理论依据强。本发明以引起捷联式航空重 力仪测量误差的原因为出发点，估计出了姿态误差这一捷联式航空重力仪测量误差的误差源， 误差分离的结果更加可以精确，更具有指导意义。

(2)本发明的捷联式航空重力仪误差分离方法分离精度高。本发明的针对性很强，只针 对引起捷联式航空重力仪测量误差的源进行误差分离，避免了有效信息的误分离。与原来的 方法相比，这样的误差分离更贴近于实际的情况，误差分离后的结果精度也更高。

(3)本发明的捷联式航空重力仪误差分离方法实现简便。本发明只需要在卡尔曼滤波的 前向估计过程中根据滤波收敛情况设计一个反向恢复过程即可，算法实现所需的额外工作量 较小。

附图说明

图1是本发明方法的流程示意图。

图2是本发明中采用前向估计—反向恢复方法的原理示意图。

图3是本发明中采用的前向估计—反向恢复方法在应用实例中的流程示意图。

具体实施方式

以下将结合说明书附图和具体实施例对本发明做进一步详细说明。

本发明的捷联式航空重力仪误差分离方法，为一种基于卡尔曼滤波的捷联式航空重力仪 误差分离的前向估计—反向恢复方法，如图1、图2和图3所示，其具体步骤为：

（1）、捷联式航空重力仪原始数据的预处理：捷联式航空重力仪得到的加速度计和陀螺 的原始数据都是脉冲信号，需要通过当量转换才能转换成正常的测量信息。

根据捷联式航空重力仪标定得到的标定参数对重力仪得到的原始脉冲数据进行当量转 换，以获取正常的测量信息；

f=K_a·N_a+B_a   (1)

w=K_g·(N_g-B_g)   (2)

其中，K_a和B_a分别表示加速度计的刻度因子和零偏；K_g和B_g分别表示陀螺的刻度因 子和零偏；N_a和N_g分别表示加速度计和陀螺的原始脉冲数据；f和w分别表示当量转换后 的测量信息。

（2）、确定捷联式航空重力仪在飞行测量过程中产生的姿态误差：通过步骤(1)获得捷联 式航空重力仪的正常测量信息之后，分以下三步确定捷联式航空重力仪的姿态误差：

（2.1）、在当地地理坐标系中捷联式航空重力仪的误差动态方程为：

$\delta \stackrel{.}{p}=\mathrm{\delta v}---\left(3\right)$

$\delta \stackrel{.}{v}=[f^n\times ]\psi +C_b^n\delta f^b-({2\omega }_{\mathrm{ie}}^n+{\omega }_{\mathrm{en}}^n)\times \mathrm{\delta v}$

    (4)

$-(2\delta {\omega }_{\mathrm{ie}}^n+\delta {\omega }_{\mathrm{en}}^n)\times v-\delta g^n$

$\stackrel{.}{\psi }=-{\omega }_{\mathrm{in}}^n\times \psi +\delta {\omega }_{\mathrm{in}}^n-C_b^n\delta {\omega }_{\mathrm{ib}}^b---\left(5\right)$

其中，δp、δv和ψ分别表示捷联式航空重力仪的位置误差、速度误差、姿态误差；fⁿ表示在当地地理坐标系下的加速度计比力测量值；表示从坐标系到当地地理坐标系的转移 矩阵；δf^b和分别表示加速度计和陀螺的测量误差；表示地球坐标系相对于惯性坐标 系的转动速度在当地地理坐标系中的投影；表示当地地理坐标系相对于地球坐标系的转 动速度在当地地理坐标系下的投影；δgⁿ表示当地地理坐标系下的重力扰动。

选取位置误差、速度误差和姿态误差为状态变量，可以得到卡尔曼滤波的状态方程：

$\frac{d}{\mathrm{dt}}\left(\begin{array}{c}\mathrm{\delta p}\\ \mathrm{\delta v}\\ \psi \end{array}\right)=F\left(\begin{array}{c}\mathrm{\delta p}\\ \mathrm{\delta v}\\ \psi \end{array}\right)+\mathrm{Gu}---\left(6\right)$

观测方程为：

$y=H\left(\begin{array}{c}\mathrm{\delta p}\\ \mathrm{\delta v}\\ \psi \end{array}\right)---\left(7\right)$

其中：

$H=\left(\begin{array}{ccc}I& 0& 0\\ 0& I& 0\end{array}\right)---\left(8\right)$

观测量为：

$Z=\left(\begin{array}{c}\mathrm{\delta p}\\ \mathrm{\delta v}\end{array}\right)=\left(\begin{array}{c}{p}_{\mathrm{GPS}}^n-p_{\mathrm{INS}}^n\\ v_{\mathrm{GPS}}^n-v_{\mathrm{INS}}^n\end{array}\right)·---\left(9\right)$

其中，表示由GPS测得的当地地理坐标系中的重力仪位置；表示由重力仪中 的加速度计和陀螺得到的重力仪位置；表示由GPS测得的当地地理坐标系中的重力仪速 度；表示由重力仪中的加速度计和陀螺得到的重力仪速度。

（2.2）、利用卡尔曼滤波方法进行状态变量X的前向估计，一个完成的卡尔曼滤波过程 可以表示如下：

X_k+1/k=φ_kX_k/k   (10)

$P_{k+1/k}={\phi }_k{P}_{k/k}{\phi }_k^T+Q_k---\left(11\right)$

X_k+1/k+1=X_k+1/k-K_k+1(H_k+1X_k+1/k-z_k+1)   (12)

$P_{k+1/k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{(I-K_{k+1}{H}_{k+1})}^T-K_{k+1}{\mathrm{PK}}_{k+1}^T---\left(13\right)$

其中，

$K_{k+1}=P_{k+1/k}{H}_{k+1}^T(H_{k+1}{P}_{k+1/k}{H}_{k+1}^T+R)---\left(14\right)$

（2.3）、利用一个周期内前向估计得到的状态变量进行反向恢复得到该周期起始点处的 状态变量的最终估计值，其中的姿态误差就是最终估计出来的姿态误差。反向恢复的计算 如下：

X_0/0=φ_0,NX_N/N   (15)

其中，X_N/N表示一个前向估计—反向恢复计算周期终点处的状态估计值；X_0/0表示一 个前向估计—反向恢复计算周期起点处的状态估计值；φ_0,N表示一个前向估计—反向恢复计 算周期终点到起点的转移矩阵。

（3）、通过步骤（2）确定的姿态误差对捷联式重力仪的测量结果进行修正，从而实现 误差的分离：

$d\tilde{g}=\mathrm{dg}+\delta f^n---\left(16\right)$

$\delta f^n=f^n\times \tilde{\psi }---\left(17\right)$

其中，表示误差分离后得到的重力扰动；dg表示误差分离前重力扰动测量结果；δfⁿ表示由于姿态误差引起的测量误差的大小。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于 本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术 人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。