基于相关分析与经验模分解的捷联式重力仪误差矫正方法

摘要

本发明属于重力测量领域，公开了基于相关分析与经验模分解的捷联式重力仪误差矫正方法。通过等效带通滤波器研究相应频段内的重力仪误差以提高方法的精度，利用经验模分解得到测量结果的本征模态分量，利用本征模态分量与伪比力的相关性来进行误差分离，并且利用迭代的方式避免经验模分解中的模态混叠，本发明在不影响捷联式航空重力测量分辨率的情况下提高了捷联式航空重力测量的精度，对于捷联式重力仪的精度提升具有较大意义。

说明书

技术领域：

本发明涉及一种捷联式航空重力仪误差矫正方法，特别涉及了基于相关分析与经验模分解的捷联式重力仪误差矫正方法。

背景技术：

地球重力场的信息在地球物理、大地测量以及大地动力学等领域具有非常重要的作用。在所有测量地球重力场的手段当中，航空重力测量是一种获取高精度大范围地球重力场的高效率手段。航空重力仪分为平台式与捷联式重力仪，捷联式重力仪同平台式重力仪最大的不同在于捷联式重力仪使用数学方式而不是实际的稳定平台来获取重力场信息。捷联式重力仪具有结构简单，体积较小，成本低以及功耗低的优势，并且前期研究工作表明捷联式航空重力仪在获取高精度高分辨率重力场方面具有同平台式重力仪相当甚至更高的表现，所以捷联式重力仪已经成为航空重力测量的趋势所在。捷联式航空重力仪面对的一个问题是，它容易受到气流以及飞行器运动特性的影响。这是由于捷联式重力仪直接固定在舱室内部没有物理平台用于隔离气流以及载体运动对捷联式重力仪的影响，这类由于气流以及载体动态特性导致的误差被称为动态性误差。

经验模分解是一种自适应的信号分解手段，通过经验模分解可以得到信号的本征模态分量。经验模分解可以用于海洋与地震数据分析、机械振动信号分析、声音处理、纹理分析、SAR影像滤波以及气象等领域，并取得了较好的结果。从目前情况来看，直接针对捷联式重力仪动态性误差的研究工作以及利用经验模分解对重力测量信号进行处理的文献较少。文献(Zhao L.et al.Airborne Gravity Data Denoising Based on Empirical Mode Decomposition:A Case Study for SGA-WZ Greenland Test Data.ISPRS Int.J.Geo-Inf.2015,4,2205-2218)应用了经验模分解对实际测量得到的捷联式重力仪数据进行了处理，得到了同航空重力测量中通常使用的低通滤波器相接近的效果，但是对测量结果精度的提高并没有显著帮助。从重力测量理论来看，动态性误差的源头不能通过卡尔曼估计的方式得到，但是动态性误差与比力是强相关的。因此，实现一种基于相关性分析与经验模分解的捷联式航空重力仪动态性误差矫正方法具有很高的实际应用价值。

发明内容:

本发明针对捷联式航空重力仪受外界气流以及载体动态特性影响，进而产生相关的动态性误差的问题，提出了一种基于相关分析与经验模分解的捷联式重力仪误差矫正方法，进而减小动态性对航空重力测量的影响，提升航空重力测量的精度。

主要的技术方案如下：

基于相关分析与经验模分解的捷联式重力仪误差矫正方法，包括以下步骤：

步骤一，对实际测量得到的重力数据进行重力解算并且经过低通滤波得到特定分辨率下的测量结果：

捷联式重力测量的模型为

其中，和为由GPS测量得到的相对地球的加速度与速度，f^b是捷联式惯导系统测量得到的在载体系下的比力信息，是由载体系到导航系的姿态变化矩阵，是地球相对于惯性系的转动角速度，是导航系相对于地球的转动角速度，γⁿ是在导航系下的正常重力值，dgⁿ是重力异常结果在导航系下的表示；

步骤二，利用等效带通滤波器对测量结果进行滤波并得到对应频段测量结果，其中等效带通滤波器的结果相当于两个低通滤波器的差值；

步骤三，利用经验模分解对重力测量的带通滤波结果进行分解，得到其本征模态分量；

步骤四，利用GPS信息计算伪比力，之后同样利用步骤二中的带通滤波器对其进行带通滤波；

伪比力的计算公式为

其中，IF是所求的不包含惯性测量误差的伪比力；

步骤五，计算步骤三中得到的带通滤波结果的本征模态分量与伪比力带通滤波结果的相关性系数，其中相关性系数最大的本征模态分量为由动态性误差构成的本征模态分量，将其从测量结果中减去以得到中间结果；

步骤六，由于经验模分解中可能存在模态混淆问题，对中间结果再次进行带通滤波，并对结果进行经验模分解，得到中间结果经带通滤波后的本征模态分量；

步骤七，计算新的本征模态分量同伪比力带通滤波结果的相关性系数，依旧认为相关性系数最大的本征模态分量为由动态性误差构成的本征模态分量应从测量结果中减去，得到最终结果；

在本发明中，通过以上七个步骤，便可以实现基于相关分析与经验模分解的捷联式重力仪误差矫正方法，在不改变分辨率的情况下分离了动态性误差。

本发明具有以下优点：

1)直接针对捷联式重力仪中的动态性误差进行矫正，实现误差分离；

2)不改变分辨率的情况下提高测量精度。

附图说明：

图1为某次重力测量实验中的高度变化情况；

图2为飞行实验中测线1、5、6的结果图；

图3为带通滤波器过后测线1、5、6的结果图；

图4为经验模分解的流程图；

图5为测线1的各个本征模态分量；

图6为矫正后飞行实验中测线1、5、6的结果图；

图7为整个发明的流程图。

具体实施方式：

下面结合附图和某次实际航空重力测量实验对本发明中的方法做进一步详细阐述，实验的组织实施情况如表1，在该次飞行实验中包含6条重复测线(以下统称飞行实验)，其中第5、6条飞行测线为起伏飞行以验证动态性误差，飞行实验的高度变化如图1所示。

表1

步骤一，对实际测量得到的重力数据进行重力解算并且经过低通滤波得到3km分辨率下的测量结果；

捷联式重力测量的模型为

其中，和为由GPS测量得到的相对地球的加速度与速度，f^b是捷联式惯导系统测量得到的在载体系下的比力信息，是由载体系到导航系的姿态变化矩阵，是地球相对于惯性系的转动角速度，是导航系相对于地球的转动角速度，γⁿ是在导航系下的正常重力值，dgⁿ是重力异常结果在导航系下的表示；

对飞行实验的数据进行重力解算，得到6条重复测线的测量精度为1.73mGal/3km(1mGal＝1*10^-5m/s²)，其中测线1、5、6的精度对比如图2。

步骤二，利用等效带通滤波器对测量结果进行滤波并得到对应频段测量结果，其中等效带通滤波器的结果相当于两个低通滤波器的差值；

将数据经过带通滤波器，图3给出了经过带通滤波器后测线1、5、6的结果。

步骤三，利用经验模分解对重力测量的带通滤波结果进行分解，得到其本征模态分量：

经验模态分解的步骤如下所示(图4为其流程图)：

(1)：x(t)为将要进行经验模分解的信号，确定x(t)的所有极大值和极小值，分别对极大值点和极小值点进行插值(例如采用三次样条进行插值计算)，构造的上下包络线x_up(t)和x_low(t)，计算上下包络线的均值m₁(t)＝(x_up(t)+x_low(t))/2；

(2)：把原信号序列x(t)减去该均值m₁(t)可得到一个去掉低频的新数据序列g₁(t)，即g₁(t)＝x(t)-m₁(t)；

(3)：判断g₁(t)是否是一个本征模态分量(intrinsic>

如果g₁(t)符合IMF的定义条件，是一个IMF，则提取g₁(t)作为第一个IMF分量，令imf₁(t)＝g₁(t)，并求原信号与imf₁之间的差值r₁(t)，r₁(t)＝x(t)-imf₁(t)。

如果g₁(t)不是一个IMF，则将g₁(t)视为一个新的信号序列，重复第一步和第二步，求其包络均值及g₁(t)与m₁₁(t)间的差值g₁(t)；对g₁(t)重复上述过程n次，直到符合IMF的定义条件，则令为x(t)的第一个IMF分量，并求原信号与imf₁之间的差值r₁(t)，即r₁(t)＝x(t)-imf₁(t)；

(4)：将r₁(t)作为一个新的“原始”信号，重复步骤1至3，提取第二个本征模态函数分量imf₂，令r₂(t)＝r₁(t)-imf₂(t)；将r₂(t)作为一个新的“原始”信号，抽取第三个本征模态函数分量imf₃；以此类推，直到第K次的余项r_K＝r_(K-1)-imf_k满足终止条件，则停止迭代，x(t)的经验模态分解完成。

经以上步骤，此时原始信号x(t)可以被表示为以上分解得到的所有本征模态函数和一个残余趋势项之和：

式中imf_k(t)表示第k个本征模态函数，每个本征模态函数imf₁，imf₂，……，imf_k表示原始信号从小尺度到大尺度的时间特征。

对各条测线进行经验模分解，得到其本征模态分量，图5给出了测线1的各个本征模态分量。

步骤四，利用式(2)计算伪比力，之后利用步骤二中相同的带通滤波器对其进行带通滤波；

步骤五，计算步骤三中得到的本征模态分量与伪比力带通滤波结果的相关性系数，其中相关性系数最大的本征模态分量为由动态性误差构成的本征模态分量，将其从测量结果中减去以得到中间结果。

测线1中本征模态分量与伪比力带通滤波结果的相关性系数如表2所示，其中第2个本征模态分量的相关性系数最大，所以认为第2个本征模态分量为由动态性误差构成的本征模态分量，在测量结果中去除,得到中间结果；

表2

本征模态分量编号123456相关性系数-0.2-0.56-0.13-0.090.040.01

步骤六，由于经验模分解中可能存在模态混淆问题，对中间结果再次利用步骤二中带通滤波器进行带通滤波，并对结果进行经验模分解，得到新的本征模态分量。

步骤七，计算新的本征模态分量与伪比力带通滤波结果的相关性系数，依旧认为相关性系数最大的本征模态分量为由动态性误差构成的本征模态分量应从测量结果中减去，得到最终结果，图7为整个算法流程图。

经过矫正后，飞行实验的精度从1.73mGal/3km提升到1.27mGal/3km，图6给出了提升后测线1、5、6的精度。以上结果充分说明了本发明中算法的突出效果。

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