技术领域
[0001] 本
发明涉及捷联式航空重力仪水平分量测量结果中误差的修正方法，特指利用重力场模型对捷联式航空重力仪水平分量测量结果中的误差进行修正的方法，属于重力测量领域。
背景技术
[0002] 传统的动态重力仪是利用双轴稳定平台来保证重力
传感器敏感轴的垂直指向。捷联式航空重力仪是利用数学平台来取代物理平台以实现重力传感器的垂直指向和水平指向。
加速度计和陀螺是捷联式航空重力仪的核心器件，由于惯性导航的固有误差特性，利用加速度计和陀螺数据进行
姿态解算会存在姿态误差，由此会在捷联式重力仪测量结果的水平分量中带来低频误差。
[0003] 目前对于姿态误差的修正方法尚未成熟，国外一些研究机构提出了诸如
波数相关滤波，端点匹配等误差分离方法。其中波数相关滤波需要对每条测线进行重复测量或者要求测线之间间距比较近；端点匹配的方法要求提供测线两端的准确的重力信息，这在实际测量过程中是很难满足的。因此，需要设计一种实用的误差修正方法来实现捷联式重力仪测量水平分量误差的有效修正。
发明内容
[0004] 本发明要解决的问题就在于：针对现有的误差修正方法需要借助本条测线之外的其他信息，实用性不强这一问题，根据捷联式航空重力仪的姿态解算误差特性，利用重力场模型设计一种新的误差修正方法。
[0005] 为了实现上述发明目的，本发明思路为：基于捷联式航空重力仪姿态解算的误差特性，并利用地球重力场模型对捷联式航空重力仪水平分量的测量误差进行修正。
[0006] 本发明的技术方案为：利用重力场模型的捷联式重力仪水平分量误差修正方法，包含以下步骤：
[0007] (1)、对捷联式航空重力仪原始数据进行预处理：将捷联式航空重力仪得到的加速度计、陀螺的原始脉冲数通过当量转换成正常的测量信息；
[0008] (2)、对GNSS原始数据进行处理得到载体的
位置、速度和加速度信息：分以下两步确定载体的位置、速度和加速度信息：
[0009] ①、通过差分
定位方式获得载体的位置信息；
[0010] ②、对载体的位置信息进行两次差分分别得到速度信息和加速度信息；
[0011] (3)、以步骤(2)得到的GNSS载体位置信息和速度信息作为外部观测，对通过步骤(1)获得的当量转换完后的加速度计和陀螺数据进行卡尔曼滤波组合导航解算，得到比力估计结果fn；
[0012] (4)、对比力估计结果fn和GNSS载体加速度信息进行求差得到重力矢量估计结果δgn；
[0013] (5)、对步骤(4)求得的重力矢量估计结果δgn进行低通滤波，将其包含的高频噪声滤除；
[0014] (6)、通过步骤(1)～(5)求得的重力矢量估计结果δgn的水平分量 和 中存在由姿态误差引起的低频测量误差，分以下五步去除水平分量估计结果中存在的由姿态误差引起的低频测量误差：
[0015] ①、以姿态误差的舒勒振荡周期84.4分钟对应的
频率为截止频率设计高通
滤波器，将水平分量 和 中包含低频误差的频段整体去掉，得到不含低频部分的水平分量估计结果；
[0016] ②、将姿态误差的舒勒周期(T＝84.4分钟)与航测飞机的飞行速度v相结合，确定步骤①中去掉的频段所对应的
波长λ；
[0017] ③、通过步骤②得到的波长λ确定重力场模型计算时所需要的阶次n；
[0018] ④、根据步骤③确定的阶次n利用地球重力场模型计算步骤①中被整体去除的低频部分；
[0019] ⑤、将步骤①得到的不含低频部分的水平分量估计结果和步骤④通过地球重力场模型计算得到的低频部分相融合，得到最终消除了姿态误差的水平分量估计结果。
[0020] 进一步地，步骤(6)中利用地球重力场模型的计算值来代替原始测量结果中具有姿态误差的低频部分。
[0022] (1)理论
基础好。本发明根据水平分量测量结果中所存在的姿态误差的
频谱特性对水平分量测量结果进行修正，针对性强，理论基础好。
[0023] (2)可操作性强。本发明不需要测区提供额外的重力场信息对测量结果进行修正，不需要增加额外的工作量，可操作性强。
[0024] (3)实现简便。本发明只需要增加一个重力场模型的计算过程，实现简便，
算法所需的额外工作量较小。
附图说明
[0025] 图1为本发明方法的流程示意图；
[0026] 图2为卡尔曼滤波组合导航解算得到比力估计结果的流程示意图；
[0027] 图3为利用地球重力场模型计算低频部分流程示意图。
具体实施方式
[0028] 以下结合附图对本发明做进一步详细说明。
[0029] 如图1所示，本发明所述的利用重力场模型的捷联式重力仪水平分量误差修正方法，其具体流程为：
[0030] 1、对捷联式航空重力仪原始数据进行预处理：将捷联式航空重力仪得到的加速度计、陀螺的原始脉冲数Na和Ng通过当量转换成正常的比力f和
角速度w，当量转换所需的标定参数加速度计刻度因子Ka、加速度计零偏Ba、陀螺刻度因子Kg和陀螺零偏Bg是通过对捷联式航空重力仪进行标定而得到。
[0031] f＝Ka·Na+Ba (1)
[0032] w＝Kg·(Ng-Bg) (2)
[0033] 2、对GNSS原始数据进行处理得到载体位置、速度和加速度信息：分以下两步确定载体的位置、速度和加速度信息：
[0034] (1)、通过差分定位方式获得载体的位置信息；
[0035] (2)、对载体的位置信息进行两次差分分别得到速度信息和加速度信息：
[0036]
[0037]
[0038] 其中，r(k-1)和v(k-1)分别表示t(k-1)时刻载体的位置和速度，r(k+1)和v(k+1)分别表示t(k+1)时刻载体的位置和速度；v(k)和a(k)分别表示t(k)时刻载体的速度和加速度，分别通过位置信息对时间的一次差分和二次差分得到。
[0039] 3、以步骤2得到的GNSS载体位置信息和速度信息作为外部观测，对通过步骤1获得的当量转换完后的加速度计比力和陀螺角速度数据进行卡尔曼滤波组合导航解算，得到加速度计的比力测量结果在当地地理
坐标系中的投影fn。
[0040] 在当地地理坐标系中捷联式航空重力仪的误差动态方程为：
[0041]
[0042] 其中，δp、δv和ψ分别表示捷联式航空重力仪的加速度计和陀螺数据进行惯性导航解算所存在的位置误差、速度误差和姿态误差；fn表示加速度计的比力测量结果在当地地b理坐标系下的投影； 表示从体坐标系到当地地理坐标系的姿态转移矩阵；δf 和 分别表示加速度计和陀螺的测量误差； 表示地球坐标系相对于惯性坐标系的转动角速度在当地地理坐标系中的投影； 表示当地地理坐标系相对于地球坐标系的转动角速度在当地地理坐标系下的投影；δgn表示当地地理坐标系下的重力扰动； 表示当地地理坐标系相对于惯性坐标系的转动角速度在当地地理坐标系中的投影，计算如下：
[0043]
[0044] 选取位置误差δp、速度误差δv、姿态误差ψ、加速度计零偏ba和陀螺零偏bg为状态变量，可以得到卡尔曼滤波的
状态方程：
[0045]
[0046] 以从GNSS中得到的位置和速度信息作为卡尔曼滤波组合导航解算的外部观测量。观测方程为：
[0047]
[0048] 其中，观测量表示如下：
[0049]
[0050] 其中， 和 分别表示由GNSS得到的位置和速度信息； 和 分别表示由捷联式航空重力仪中的加速度计和陀螺进行惯性导航解算得到的位置和速度。
[0051] 通过卡尔曼滤波估计出系统的误差状态之后(如图2所示)，一方面可以用位置误差δp、速度误差δv和姿态误差ψ的估计值对惯性导航解算中所存在的位置误差、速度误差和姿态误差进行修正，另一方面还可以用加速度计零偏ba和陀螺零偏bg的估计值对器件误差进行修正，由此得到更为精确的比力估计结果fn。
[0052] 4、对比力估计结果fn和GNSS载体加速度信息进行求差得到重力矢量估计结果δgn：
[0053]
[0054] 其中， 表示载体的运动加速度； 和 可以通过GNSS的位置信息和速度信息进行计算；γn表示正常重力，可以通过GNSS的位置信息进行计算。
[0055] 5、对步骤4求得的重力矢量估计结果δgn进行低通滤波，将其包含的高频噪声滤除。
[0056] 6、通过步骤1～5求得的重力矢量估计结果δgn的水平分量 和 中存在由姿态误差引起的低频测量误差。
[0057] 其中， 表示重力矢量水平分量中的北向分量； 表示重力矢量水平分量中的东向分量。
[0058] 分以下五步去除水平分量估计结果中存在的由姿态误差引起的低频测量误差：
[0059] (1)、姿态误差以84.4分钟为周期进行舒勒振荡，因此在水平分量 和 中，周期大于84.4分钟的成分里含有姿态误差引起的低频误差，以该周期对应的频率为截止频率设计高通滤波，由此将水平分量测量结果中包含低频误差的频段整体去除，得到不含低频部分的水平分量估计结果。
[0060] (2)、将姿态误差的舒勒周期(T＝84.4分钟)与航测飞机的飞行速度v相结合，确定步骤(1)中去掉的频段所对应的波长λ：
[0061] λ＝T·v (11)
[0062] (3)、通过步骤(2)得到的波长λ确定重力场模型计算时所需要的阶次n：
[0063]
[0064] (4)、根据步骤(3)确定的阶次n利用地球重力场模型计算步骤(1)中被整体去除的低频部分：
[0065]
[0066]
[0067] 其中， 和 分别表示利用地球重力场模型计算出的重力矢量的北向分量和东向分量；GM表示地球引力常数和地球
质量的乘积；a表示平均地球半径；r表示计算点的向径；θ表示计算点的余纬；λ表示计算点的经度； 表示规范化的勒让德多项式； 和表示完全正常化的地球扰动引力位系数；n表示重力场模型的阶数；m表示重力场模型的次数。
[0068] (5)、将步骤(1)得到的不含低频部分的水平分量估计结果和步骤(4)通过地球重力场模型计算得到的低频部分(如图3所示)相融合，得到最终消除了姿态误差的水平分量估计结果。
[0069] 以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述
实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
