基于经验模态分解的航空重力测量动态误差分离

引用本文

邹欣蕾, 蔡劭琨, 吴美平, 曹聚亮, 张开东. 基于经验模态分解的航空重力测量动态误差分离[J]. 物探与化探, 2016,40(6): 1217-1221.
ZOU Xin-Lei, CAI Shao-Kun, WU Mei-Ping, CAO Ju-Liang, ZHANG Kai-Dong. Dynamic errors separation of airborne gravimetry based on empirical mode decomposition[J]. Geophysical and Geochemical Exploration, 2016,40(6): 1217-1221.
Doi:10.11720/wtyht.2016.6.25 复制到剪切板

Permissions

基于经验模态分解的航空重力测量动态误差分离

邹欣蕾, 蔡劭琨, 吴美平, 曹聚亮, 张开东

国防科学技术大学机电工程与自动化学院,湖南长沙 410073

作者简介: 邹欣蕾(1990-),女,国防科学技术大学机电工程与自动化学院导航技术实验室硕士研究生,主要从事航空重力测量数据处理方法研究工作。

收稿日期: 2016-01-25

基金: 国家高技术研究发展计划(“863”计划)项目(2013AA063902)

摘要

航空重力测量数据在测量的过程中很容易受到飞机动态运动的影响,飞机动态运动强时,航空重力测量数据的质量也随之降低,虽然经过数字滤波后可以将航空重力异常数据中的大部分噪声去除,但一部分与信号重叠的噪声还是无法被消除。利用经验模态分解对受飞机动态影响较大的一组实测数据进行深化处理。结果表明,在不影响航空重力异常数据分辨率的情况下,该方法很好地消除了由飞机动态运动引起的动态误差,使得航空重力测量的内符合精度由1.43 mGal提升到了1.27 mGal。

关键词: 航空重力测量; 经验模态分解; 重力异常提取; 相关分析

中图分类号:P631 文献标志码:A 文章编号:1000-8918(2016)06-1217-05 doi: 10.11720/wtyht.2016.6.25

Dynamic errors separation of airborne gravimetry based on empirical mode decomposition

ZOU Xin-Lei, CAI Shao-Kun, WU Mei-Ping, CAO Ju-Liang, ZHANG Kai-Dong

Mechatronies and Automation College,National University of Defense Technology,Changsha 410073,China

Abstract

The data of airborne gravimetry measurement is easily interfered by the dynamic movement of planes,and the quality of the data of airborne gravimetry will also decline when the dynamic movement of planes is strong.Through digital filtering, most of the noises in gravity anomaly data can be removed.However,a part of noises overlapped with useful signals cannot be wiped out.Using empirical mode decomposition,we can process a group of measured data in depth which have great influence on dynamic movement of planes.The result indicates that the dynamic errors caused by the dynamic movement of planes could be eliminated without affecting the data resolution of airborne gravity anomaly,thus enhancing the internal accord accuracy from 1.43 to 1.27 mGal.

Keyword: airborne gravimetry; empirical mode decomposition; extraction of gravity anomaly; correlation analysis

Show Figures

在航空重力测量中, 由于各种外界因素的影响, 使得观测数据的有效信息往往淹没在大量的噪声之中, 如何降低甚至消除噪声, 是航空重力数据处理中一项关键技术。消除动态环境下的强噪声影响并保留有用的重力信息, 是滤波设计中应当重点考虑的问题^{[1, 2, 3, 4, 5]}。

传统的FIR(finite impulse response, 有限冲激响应)低通滤波器已经在航空重力测量数据处理中得到了较为广泛的应用, 其中窗函数法是最常见的滤波器设计方法^{[6, 7, 8, 9]}。FIR低通滤波器可以将数据中的高频噪声去掉, 对于噪声和有用信号在频域有明显区别的信号来说, FIR低通滤波器可以很有效地将噪声去除。而航空重力异常数据的特性除了噪声的幅值很大、有用信号主要集中在低频段的很窄部分之外, 噪声与有用信号之间还存在重叠, 若要使用FIR低通滤波器将这部分误差消除, 就必须要减小截止频率, 这无疑会降低分辨率, 并不是一个理想的结果。在测量的过程中, 飞机最佳的状态为匀速、等高的直线运动, 但由于多种内、外界因素的影响, 使得在飞行状态上会有一定的波动, 而且波动越剧烈, 对所采集到的数据影响越大^{[10, 11, 12]}。这些由飞机的动态运动所引起的误差主要集中在低频段, 仅通过FIR低通滤波器无法在保证分辨率不降低的情况下将其有效消除。经验模态分解方法认为任何复杂的时间序列都可以分解为一组从高频到低频的本征模态函数, 这些本征模态函数可以很好地反映信号在任何时间局部的频率特征。笔者将经验模态分解应用到航空重力测量数据处理中, 对已经过FIR低通滤波的航空重力测量数据做深化处理, 并使用实测数据对这一方法进行验证。

1 航空重力测量的数学模型

经典力学中, 根据牛顿第二定律, 在惯性坐标系i中, 质点的运动方程^[13]:

$\begin{matrix} {\ddot{x}}^{i} = f^{i} + g^{i}, (1) \end{matrix}$

其中:fⁱ为惯性系i系下的比力测量值, 由加速度计测量得到; gⁱ为惯性系i系下的质点所受的引力加速度; $\begin{matrix} {\ddot{x}}^{i} \end{matrix}$ 为惯性系i系下的质点运动加速度。

将式(1)转换为重力测量方程:

$\begin{matrix} g^{i} = {\ddot{x}}^{i} - f^{i}, (2) \end{matrix}$

将式(2)转换到当地地理坐标系n系下可得方程:

$\begin{matrix} g^{n} = {\dot{v}}^{n} + (2 ω_{ie}^{n} + ω_{en}^{n}) v^{n} - f^{n}, (3) \end{matrix}$

其中: $\begin{matrix} {\dot{v}}^{n} \end{matrix}$ 为n系载体运动加速度; vⁿ为载体运动速度; $\begin{matrix} ω_{ie}^{n} \end{matrix}$ 为地球自转角速度在n系下的投影; $\begin{matrix} ω_{en}^{n} \end{matrix}$ 为n系相对地球坐标系e系的角速度在n系下的投影; (2 $\begin{matrix} ω_{ie}^{n} \end{matrix}$ + $\begin{matrix} ω_{en}^{n} \end{matrix}$ )vⁿ为Coriolis加速度; fⁿ为n系下的比力测量值; gⁿ为n系下的重力加速度矢量。

式(3)中的重力加速度矢量gⁿ又可以表示为正常重力矢量γ ⁿ和重力扰动矢量δ gⁿ之和, 因此可得出航空重力矢量测量的基本模型为:

$\begin{matrix} δ g^{n} = {\dot{v}}^{n} - C_{b}^{n} f^{b} + (2 ω_{ie}^{n} + ω_{en}^{n}) v^{n} - γ^{n}, (4) \end{matrix}$

其中: $\begin{matrix} C_{b}^{n} \end{matrix}$ 为载体坐标系b系到当地地理坐标系n系的方向余弦转换矩阵; f^b为载体系下的比力测量值。

可将式(4)写为重力扰动值在当地地理坐标系3个方向上的分量形式:

$\begin{matrix} \begin{matrix} δ g_{N} = {\dot{v}}_{N} - f_{N} + 2 ω_{ie} \cdot v_{E} \cdot sinL - \\ \frac{v_{N} \cdot v_{D}}{R_{M} + h} + \frac{v_{E}^{2} \cdot tanL}{R_{N} + h}, (5) \\ δ g_{E} = {\dot{v}}_{E} - f_{E} + 2 ω_{ie} \cdot (v_{N} \cdot sinL + v_{D} \cdot cosL) - \\ \frac{v_{E} \cdot v_{N} \cdot tanL}{R_{N} + h} - \frac{v_{E} \cdot v_{D}}{R_{N} + h}, (6) \\ δ g_{D} = {\dot{v}}_{D} - f_{D} + (2 ω_{ie} \cdot cosL + \frac{v_{E}}{R_{N} + h}) \cdot v_{E} + \\ \frac{v_{N}^{2}}{R_{M} + h} - γ, (7) \end{matrix} \end{matrix}$

其中:δ g_N为重力扰动的北向分量; δ g_E为重力扰动的东向分量; δ g_D为重力扰动的垂向分量; v_N、v_E、v_D为载体速度在北向、东向、地向上的分量; $\begin{matrix} {\dot{v}}_{N} \end{matrix}$ 、 $\begin{matrix} {\dot{v}}_{E} \end{matrix}$ 、 $\begin{matrix} {\dot{v}}_{D} \end{matrix}$ 为载体加速度在北向、东向、地向上的分量; f_N、 f_E、 f_D为比力在北向、东向、地向上的分量; L为纬度; h为高度; R_M子午圈半径; R_N为卯酉圈半径。通常所说的航空重力标量测量就是获取重力扰动的垂向分量δ g_D, 又称为重力异常。

在某次航空重力重复测线测量中3条重复测线 (L61501、L61502、L61508)上所获得的重力异常结果如图1所示。这3条重复测线的重力异常结果均已使用同一FIR低通滤波器将高频噪声滤除。

	Figure Option View Download New Window
	图1 三条测线重力异常数据

由图1可以看出, 经过FIR低通滤波器处理后的测线L61501和测线L61502相对平滑, 效果较好; 但测线L61508的振荡相比起前两条测线要大很多, 说明在测线L61508的重力异常结果中还含有一定量的噪声。通过对3条测线的飞行高度变化这一表征飞机飞行动态性的数据进行分析, 认为测线L61508中的振荡可能是由飞机飞行中相对剧烈的动态运动引起的, 属于动态误差。

3条测线的飞行高度如图2所示:由图2可以看出, 与测线L61501和L61502相比, 测线L61508的高度起伏变化最剧烈, 尤其是中后段, 上下起伏的幅度比较大, 从直观上来看, 这与测线L61508上重力异常结果的振荡具有很强的相关性。

	Figure Option View Download New Window
	图2 三条测线的飞行高度

为了消除测线L61508中的动态误差, 需要进一步对重力异常数据做深化处理, 此时如果再使用FIR低通滤波器, 想要滤掉剩余噪声(动态误差)就只能降低截止频率, 这样一来会降低分辨率。因此, 笔者使用经验模态分解的方法对重力异常数据做进一步的处理。

2 经验模态分解的原理

经验模态分解(empirical mode decomposition, EMD)是一种自适应的信号时频分解方法, 能够根据数据特性把数据分解成一系列的模态函数^[14]。EMD分解方法假设任何复杂的信号都是由一些具有不同带宽、反映信号局部频率特征的本征模态函数(intrinsic mode function, IMF)组成, 因此可以基于信号本身特性自适应地将原信号逐次分解, 获得一组J个从高频到低频的基本序列IMF以及一个残余信号r(t) 。其中, 分解得到的IMF是自适应的, 各个IMF所包含的频率范围, 即分辨率也是自适应的。IMF的自适应性是因为信号的EMD分解过程是依据信号自身特性进行的, 从信号分析的角度来看, IMF就是一组频率和幅度都可变的基函数, 是完全自适应的。分解结果如式(8)所示:

EMD分解步骤(图3):

	Figure Option View Download New Window
	图3 经验模态分解算法流程

1) 确定信号中x(t)中的极小值点和极大值点。

2)采用三次样条法分别对极大值点和极小值点进行插值, 构造信号x(t)的上下包络线x_max(t)和x_min(t)并计算上下包络线的均值:

$\begin{matrix} m (t) = \frac{(x_{\max} (t) + x_{\min} (t))}{2} 。 (9) \end{matrix}$

3) 计算信号x(t)与包络线均值m(t)的差值:

$\begin{matrix} h (t) = x (t) - m (t) 。 (10) \end{matrix}$

4) 判断差值h(t)是否满足IMF的条件。若满足, 则将h(t)作为信号x(t)的一个IMF分量, 并求出信号x(t)与IMF分量的差值:

$\begin{matrix} r (t) = x (t) - h (t) 。 (11) \end{matrix}$

继续执行步骤5); 若不满足, 则将h(t)作为新的信号重新返回执行1)~4)。

5) 将r(t)作为新的信号重新返回执行1)~4),

循环执行提取出所有的IMF分量。直到提取的最后一个IMF分量或者余项小于预先设定的阈值或已成为单调函数即可结束算法。

其中每个IMF分量必须满足的条件:①整个信号上, 极值点数目与零点数目相差1; ②极大值与极小值分别构成的包络线均值为0。

3 基于经验模态分解的动态误差分离

采用EMD方法对航空重力测量数据进行去噪处理, 首先对数据进行EMD分解, 然后根据相关性分析对IMF进行筛选, 最后进行重构得到去噪后的信号, 基本步骤为(图4):

	Figure Option View Download New Window
	图4 经验模态分解的动态误差分离流程

1)将动态误差较大的重力数据进行EMD分解, 得到了一组本征模态函数IMF和残余量。

2)采用特定的IMF筛选准则, 将认为是噪声的分量从IMF中去除。首先将各个IMF分量与表征飞机动态的物理量(如飞行高度的变化)进行相关分析, 相关性强的IMF分量认为是由飞机运动引起的动态误差, 最后得到一组新的IMF。

3)将新的IMF与原有的残余量进行叠加, 叠加的结果就是去除动态误差后的重力异常结果。

3.1 EMD分解

在前面的部分中已经介绍了EMD分解的原理, 下面就对测线L61508的重力异常结果进行EMD分解, 得到如图5所示的分解结果。

	Figure Option View Download New Window
	图5 测线L61508重力异常经验模态分解的IMF分量(mGal)

图中IMF1和IMF2为分解得到的两个IMF分量, 这两个IMF分量是通过前面所述的EMD分解方法自动得到的, residual为分解后的残余信号。

3.2 IMF筛选

下面对分解得到的两个IMF分量做相关性分析, 分析结果将作为筛选的依据。相关分析包括两个IMF分量与飞机飞行高度变化的相关分析以及两个IMF分量与测线L61501的重力异常结果的相关分析。与飞机飞行高度变化进行相关分析是为了判断两个IMF分量中哪个是噪声分量, 与测线L61501作相关分析是为了判断哪个分量与重力异常有用信号相关性强, 选择测线L61501的重力异常结果是因为其动态误差较小, 数据真实性好。在实际应用中, 仅需要做IMF分量与飞机飞行高度变化的相关性分析, 在此做IMF分量与测线L61501的相关性分析是为了更加透彻说明筛选原因。分析结果如表1所示。

表1 相关性分析结果

从上表可以看出, IMF1和飞行高度变化的相关性要比IMF2与高度误差的相关性大很多, 而IMF1与测线L61501的相关性则比IMF2与测线L61501的相关性小了近一个数量级。由此可以说明, IMF1中包含的主要是由飞机的动态运动所引起的动态误差, 在筛选中可以将其去除。

3.3 数据重构与结果分析

对经过筛选后保留下的部分进行叠加, 得到深化处理后的航空重力异常结果, 再与处理前的结果进行对比(图6), 可以看出经过EMD分解去噪后, 重力异常结果中的振荡得到了很好的分离。

	Figure Option View Download New Window
	图6 测线L61508重力异常数据深化处理前后对比

对处理前后的数据进行进一步的频谱分析, 图7所示的是测线L61508重力异常结果处理前后的功率谱密度, 从图中可以看出, 经过EMD分解去噪后的重力异常结果的分辨率并没有降低, 分辨率均为4 km左右, 说明在去除动态误差的同时有效地保留了有用信号。

	Figure Option View Download New Window
	图7 测线L61508重力异常数据深化处理前后功率谱密度对比

最后, 对深化处理前后测线L61501、测线L61502和测线L61508的内符合精度进行计算, 得深化处理前3条测线的内符合精度为1.43 mGal, 深化处理后3条测线的内符合精度为1.27 mGal, 相对处理前, 内符合精度有一定程度上的提高。

4 结论

航空重力测量在人力、效率、作业便利等方面都具有无可比拟的优越性, 有效地弥补了传统地面重力测量、海洋测量以及卫星重力测量的诸多不足。由于测量数据混有高频噪声, 使得有效的滤波算法成为提升航空重力测量数据质量的关键。通过分析经FIR低通滤波器预处理后的航空重力异常数据, 最终确定使用经验模态分解去噪的方法来对航空重力异常数据进行深化处理, 并通过分析相关性确定了本征模态函数的筛选依据, 重构后得到了去除动态误差后的航空重力异常结果, 通过分析数据的功率谱密度和内符合精度可知, 在没有降低数据分辨率的情况下将数据的内符合精度由1.43 mGal提升至1.27 mGal, 说明了处理是有效的。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	张开东. 基于SINS/DGPS的航空重力测量方法[D]. 长沙: 中国人民解放军国防科学技术大学, 2007. [本文引用:1]
[2]	张开东, 吴美平, 胡小平. 基于捷联惯导的航空重力测量滤波算法[J]. 中国惯性技术学报, 2007, 15(1): 5-8. [本文引用:1]
[3]	孙中苗. 航空重力测量理论、方法及应用研究[D]. 郑州: 中国人民解放军信息工程大学, 2004. [本文引用:1]
[4]	Youcef H. A Comparison of Filtering Techniques for Airborne Gravimetry[D]. Calgary: University of Calgary, 1996. [本文引用:1]
[5]	Bruton A M. Improving the accuracy and tesolution of SINS/DGPS airborne gravimetry[D]. Calgary: University of Calgary, 2000. [本文引用:1]
[6]	蔡劭琨, 吴美平, 张开东. 航空重力测量中FIR低通滤波器的比较[J]. 物探与化探, 2010, 34(1): 74-78. [本文引用:1]
[7]	孙中苗, 夏哲仁. 航空重力测量中低通滤波器参数的优选[J]. 测绘学院学报, 2003, 19(1): 11-14, 18. [本文引用:1]
[8]	田颜锋. 航空重力测量数据滤波处理算法研究[D]. 郑州: 中国人民解放军信息工程大学, 2010. [本文引用:1]
[9]	李姗姗. 航空重力测量中数字滤波器的设计与应用[J]. 解放军测绘学院学报, 1999, 16(2): 83-85, 89. [本文引用:1]
[10]	蔡劭琨. 航空重力矢量测量及误差分离方法研究[D]. 长沙: 中国人民解放军国防科学技术大学, 2014. [本文引用:1]
[11]	孙中苗, 夏哲仁, 石磐, 等. 航空重力测量数据的滤波与处理[J]. 地球物理学进展, 2004, 19(1): 119-124. [本文引用:1]
[12]	黄杨明. 高精度捷联式航空重力仪误差估计方法研究[D]. 长沙: 中国人民解放军国防科学技术大学, 2013. [本文引用:1]
[13]	熊盛青. 航空重力勘探理论方法及应用[M]. 北京: 地质出版社, 2010. [本文引用:1]
[14]	赵磊. 捷联式航空重力测量数据滤波技术研究[D]. 长沙: 中国人民解放军国防科学技术大学, 2015. [本文引用:1]

2007

0.0