摘要
近年来，基于光学陀螺惯导系统的动态重力测量技术能高效获取地球重力场数据，得到快速发展；但现有滤波方法对重力信号与噪声的频域混叠处理仍具有局限性，严重限制了测量的精度与分辨率。为此，创新性地提出基于经验模态分解（empirical mode decomposition，EMD）与波数相关滤波（wavenumber correlation filtering，WCF）的联合降噪方法，先对原始重力结果进行低通滤波，再对滤波结果进行经验模态分解并设置阈值保留低阶本征模态函数，对重复测线保留结果进行相关性滤波，实现信号的重构与降噪。仿真结果表明，所提方法较传统EMD信噪比提高了44%、均方误差减小了48%，证明了EMD-WCF联合降噪方法的有效性，能兼顾重力测量精度与分辨率的同时，为高频重力信号降噪提供了新的解决方案。
Abstract
In recent years, the dynamic gravity measurement technology utilizing optical-gyro inertial navigation systems has advanced rapidly, owing to its high efficiency capabilities in acquiring Earth′s gravity field data. However, existing filtering methods exhibit limitations in addressing the frequency-domain aliasing of gravity signals and noise, which severely constrains the accuracy and resolution of measurements. To address these challenges, a novel joint noise reduction method integrating EMD (empirical mode decomposition) and WCF (wavenumber correlation filtering) was innovatively proposed. The raw gravity results were low-pass filtered initially. The filtered results then underwent EMD, and a threshold was applied to retain low-order intrinsic mode functions. Subsequently, correlation filtering was implemented on repeated line retention results to reconstruct signals and suppress noise. Simulation experiments demonstrate that the proposed method achieves a 44% improvement in signal-to-noise ratio and a 48% reduction in mean square error compared to traditional EMD processing. These results confirm the effectiveness of the EMD-WCF method in balancing measurement accuracy and resolution, offering a new strategy for high-frequency gravity signal denoising.
近年来，动态重力测量技术发展迅猛，能够快速获得大范围、高分辨率的地球重力场信息，特别是采用激光陀螺/光纤陀螺和石英挠性加速度计的重力仪在海空重力测量方面得到广泛的应用，填补了静态重力测量和卫星重力测量之间的空白[1-5]。在实施动态测量过程中，飞机、车辆等载体的发动机振动或自身运动都会带来误差，造成原始观测结果中背景噪声的幅值远高于目标信号[6]，如何在强噪声干扰下提取出微弱重力信息是进一步发展动态重力测量亟须解决的问题。
由于重力测量对实时性要求不高、允许事后处理，已有学者对动态重力测量数据处理方法展开研究。针对部分频率远高于真实信号的噪声，通常采用低通滤波的方式去除，工程中使用的数字滤波器主要分为有限冲激响应（finite impulse response，FIR）低通滤波器和无限冲激响应（infinite impulse response，IIR）低通滤波器两种。文献[7-8]描述IIR滤波器处理重力数据时，由于其通带平坦、阻带衰减特性快，能够获得较为平滑的结果，在一些航空重力测量作业中得到推广[6]。但IIR滤波器不具有线性相位，无法保证滤波前后波形一致性，降低了滤波处理的精度[9]，因此，重力测量中普遍使用具有良好线性相位特性的FIR低通滤波器[10-13]。除了FIR和IIR，柳林涛等[14-15]提出一种函数滤波器，属于改进的小波滤波方法，实现了截止尺度参数与过渡带参数的分离，使滤波器设计和滤波实施变得简明灵活。郎骏健等[16]提出了傅里叶基追踪低通滤波方法，在精度相当条件下对比传统FIR处理过程，边缘数据丢失减少，提高了数据利用率。
针对部分频率与真实信号混叠的噪声，低通滤波方法处理则存在缺陷。如果想进一步提高测量精度，就需要降低滤波器截止频率，增大滤波带宽，将信号与噪声混合的部分滤除；或采用经验模态分解（empirical mode decomposition，EMD）方法[17]，直接去除一些本征模态函数（intrinsic mode function，IMF）。Cheng等[18]分别采用200 s和300 s低通滤波处理相同测线的原始数据，得到不同空间分辨率的精度结果。邹欣蕾等[19]利用EMD方法剔除了一些与飞机运动相关性弱的IMF，保留剩余分量重构出降噪后的信号。两种方法本质上都是通过舍弃重力信号的高频部分来提升重复测线内符合精度，但包含有用信号特别是反映小尺度地质变化的高频信息会丢失，造成测量结果失真，不利于构建高分辨率地球重力场数据库。
为了兼顾动态重力测量的精度与分辨率，本文研究了一种基于EMD的改进重力数据处理方法，以期实现对传统低通滤波后的信号进一步降噪，在不损失高频信息的条件下提高测量精度。
1 联合降噪方法实现
EMD作为一种时域处理方法，具有高度自适应性，它不需要预设基函数或窗口，而是根据信号的局部特征进行分解，能够准确提取信号的瞬时频率和局部特征。因此，EMD适合处理非平稳重力信号。波数相关滤波（wavenumber correlation filtering，WCF）则是统计学中基于相关性分析的一种方法，同样适合处理重复测线重力结果[20-22]。下文将简要介绍EMD与WCF技术的特点及实现方法，并提出一种基于EMD和WCF的联合降噪方法解决高频重力信号误差耦合问题。
1.1 经验模态分解
20 世纪90年代，学者黄锷为了应对经典傅里叶变换无法同时分辨信号的时间和频率问题，提出了按照数据自身时间尺度特征分解信号的方法，称为经验模态分解，该方法将复杂信号分解为有限个IMF和残差（residual，Res），各IMF代表原信号不同时间尺度的局部特征，Res是分解出所有IMF后的残余信号，能够描述信号的整体趋势[23-24]。其中，单一IMF反映信号某一时间尺度的频率特性，多个IMF共同组成了整个信号的频率特性，所以EMD适合处理频率随时间变化的非平稳信号。
根据文献[25]，经验模态分解流程如图1所示，具体步骤如下：
步骤1：确定待处理信号x（t）的所有极大值和极小值点。
步骤2：通过样条曲线拟合对局部极大、极小值分别拟合出上、下包络线，并计算出均值m（t）。
步骤3：计算原始信号与上一步求出的包络均值的差，形成新信号h（t）=x（t）-m（t）。
步骤4：判断差值后的信号h（t）是否满足本征模态函数条件，如果满足，h（t）就作为x（t）其中一个IMF，并在x（t）中扣掉该分量，将剩余信号r（t）=x（t）-h（t）代替x（t）返回步骤1继续执行，不断迭代直到剩余分量单调或小于设定的阈值，形成残差rn；如果不满足条件，则将h（t）信号替换原始待处理信号x（t），从步骤1开始执行。
图1EMD分解流程
Fig.1Flow chart of EMD decomposition
步骤5：经过数次迭代形成了若干IMF和残差，最终信号可表示为，表示EMD完成，分解得到的IMF与残差组合能够复现原始信号。
1.2 波数相关滤波
实施重力测量作业时，通常对同一条测线进行多次重复测量，采用计算内符合精度的方法评估系统的稳定性[26]。如果没有噪声干扰，多次测量相同测线敏感到的重力信号应该是不变的，而随机噪声和误差会随着测次和载体运动方向改变等存在差别。基于位置不变、重力不变原理，衍生出利用相关性分析处理重复测线数据的方法，将多测线结果相关性强的部分认为是重力信号，相关性弱的部分认为是随机噪声，通过设置阈值剔除相关性弱的分量，只保留相关性强的分量，实现信号与噪声的分离[20-22]。
假设两条重复测线重力结果分别为X（t）和Y（t），通过傅里叶变换将信号转换到频域，得到它们在频域中的不同波数分量，可以将两个信号的第k级波数分别表示为：
(1)
(2)
其中：FX（k）、FY（k）分别表示X（t）、Y（t）对应的频域第k级频率分量；|FX（k）|、|FY（k）|表示对应频率下的幅值；和分别代表对应频率的相位。通过计算相应波数下两组数据的相关系数Ck来构建它们的波数相关频谱，计算方法为：
(3)
其中，Ck表示k波数下的相关系数，为FX（k）、FY（k）之间的相位差，×表示矢量积。
基于计算出的k级相关系数，给定一个阈值，如果对应波数下的相关系数小于这个阈值，则将该波数分量认为是噪声并剔除。在处理重力数据时，阈值符号的选取同时要参照载体的运动方向。如果待处理数据源于运动方向相同的测线，例如都是从起点到终点，阈值的选择在正相关系数范围，如T＞0.7；如果是连续两条重复测线测量，从起点到终点后再由终点返回至起点，则阈值的选择在负相关系数范围，如T≤-0.7。以正相关系数为例，剔除弱相关分量如式（4）所示：
(4)
此外，阈值幅值的选取也会对最终结果产生影响：为了保留较高相关性的波数构成信号，实际应用中需要设置阈值下限；同时，为了符合真实重力测量特性，阈值过高可能导致部分具有弱相关特性的真实信号被误判为噪声而滤除，需要设置阈值上限。文献[22]提出一种阈值确定方法：在仿真中通过对所有可能选择的阈值集合（如0.3≤T≤0.8，ΔT=0.1）的降噪效果进行评估，比较不同相关滤波结果与真实信号的均方根误差，选择最小均方根误差对应的阈值。实测数据则在阈值上、下限范围内选择重构后重力内符合精度最高的参数进行滤波。
WCF属于频域信号处理，两组数据的所有波数分量筛选后，利用快速傅里叶反变换（inverse fast Fourier transform，IFFT）将滤波后的重力信号恢复到时域中，得到这条重复测线上的重力处理结果：
(5)
利用WCF处理的流程如图2所示。
图2WCF处理重力信号流程
Fig.2Flow chart of WCF processing gravity signal
1.3 EMD-WCF联合降噪处理
为了兼顾测量精度与分辨率，本文提出EMD-WCF联合降噪方法，对重力数据进行处理，具体步骤如下：
步骤1：原始重力测量结果首先进行FIR低通滤波，得到滤波后的重力信号。
步骤2：对滤波后的结果进行EMD处理，选择部分低阶IMF构成高频待降噪分量。
步骤3：对筛选出的高频待降噪分量进行波数相关滤波，得到降噪后的低阶IMF。
步骤4：将降噪后的低阶IMF与筛选前剩下的高阶IMF重构出降噪后的重力结果。
处理实测数据时，由于不同测线测量时载体运动状态不同以及系统偏差影响，会存在分解出的IMF数目不同的情况，步骤2中待降噪IMF分量的选择影响了EMD-WCF算法的精度。文献[27]提出地球重力场模型的长波结果相对准确，可以利用地球重力场模型补偿重力低频分量测量误差。因此，通过导出测线位置的模型重力值作为低频基准，选取与基准处于不同频段的IMF作为待WCF处理的高频信号。
传统的EMD降噪方法会丢失高频重力信息，因此，创新性地提出基于EMD与WCF组合的高频重力信号降噪方法，对EMD得到的低阶IMF进行相关性滤波，解耦高频信号与噪声，提取真实重力，理论上能够在不降低分辨率的情况下提高测量精度。该方法具体流程如图3所示。
图3EMD-WCF降噪方法流程
Fig.3Flow chart of EMD-WCF noise reduction method
2 仿真数据分析
为了方便对比降噪处理前后信号的信噪比（signal to noise ratio，SNR）和均方误差（mean squared error，MSE），利用仿真数据验证所提方法的有效性。仿真数据基于已发表论文中动态重力测量获取的真实重力特性设置：一方面，实际测量中载体的运动噪声频率较高，需要采用截止频率较低的滤波器进行预处理，采用300 s低通滤波已成为一种常用的参数选择[28-29]。大量相关研究都采用该参数并取得了较好的成果，选择300 s滤波也能方便实验数据处理结果与已发表成果进行对比。因此，重力的最高频率设置为1/300 Hz，对应原始测量值经过300 s低通滤波预处理的结果。另一方面，不同位置的重力信号存在差异，其物理特性在空间域中属于多波长信号的叠加，载体匀速作业时测量结果反映在时域中为调频信号，即频率随着时间变化。因此，确定最高频为1/300 Hz后，生成一段频率逐渐降低的信号模拟重力。
首先生成一段1 000 s的啁啾信号表示原始低通滤波后的结果O（t），重力测量中低通滤波器截止频率fs选择1/300 Hz，设置生成信号的频率f变化范围为，同时在1/300 Hz处添加额外的固定频率信号作为高频信息，如图4所示。
图4生成的原始低通滤波后的重力信号
Fig.4The generated original low-pass filtered gravity signal
为了模拟真实测量值，对生成的重复测线重力信号X（t）和Y（t）分别添加不同幅值的固定频率噪声，如式（6）所示：
(6)
其中：E1-high（t）表示最高频耦合噪声，频率设置为O（t）的最高频率，模拟截止频率处的信号失真和边界噪声；E2-high（t）表示次高频噪声，其频率略低于O（t），与E1-high（t）组合形成宽频带噪声；KX和KY表示不同的幅值；RX（t）和RY（t）表示白噪声，仿真测量过程其他因素的影响。添加误差后的两个信号如图5所示。
对信号X和Y进行EMD处理，得到7个IMF分量和一个Res残差（IMF8），分解后的结果如图6和图7所示。仿真数据中两个信号分解得到的图像能够直观表现出相似的频率分布特性，前5个IMF分量幅值和振荡频率比较大，说明信号高频分量集中在这个部分。而后2个IMF分量（不包括残差）形状几乎相同且振荡频率较低，对应信号的低频分量，也是信号能量的主要集中部分。因此，选择X和Y的前5个IMF构成新的信号进行WCF处理，将处理后的信号与剩余IMF和残差组合得到降噪后的结果。
图5添加误差后的模拟重力信号
Fig.5Simulated gravity signal after adding error
图6X信号EMD处理后的结果
Fig.6Results of X signal after EMD processing
提取低阶IMF分量重构后的信号如图8所示。通过图8可以看出，经过EMD筛选重构出的高频信号起始幅值偏小，说明这部分信号包含的能量较小、长波重力信息较少，这也是大多数学者在研究重力数据处理时选择直接丢弃低阶IMF的原因。这对于处理垂向重力扰动数据来说尚可以接受，但是水平重力扰动的高频段反映垂线偏差信息[30]，如果直接丢弃，不利于进一步发展矢量重力测量技术。采用本文提出的方法对高频段耦合噪声部分的重力信号进行降噪处理，利用波数相关滤波提取真实重力，降噪处理后恢复的高频段结果如图9所示。
图7Y信号EMD处理后的结果
Fig.7Results of Y signal after EMD processing
图8经过EMD筛选后的高频重力信号
Fig.8The high frequency gravity signal after EMD screening
图9EMD-WCF方法提取出的高频重力信号
Fig.9The high frequency gravity signal extracted by EMD-WCF method
图9表明，提取后的信号相比于图8中提取前的结果，噪声波动明显减小，将提取后的信号与EMD处理剩余的IMF和残差组合，重构出降噪后的重力结果，如图10所示。
图10EMD-WCF方法降噪后的重力信号
Fig.10The gravity signal after denoising by EMD-WCF method
将EMD-WCF联合降噪算法与传统的EMD处理的结果以及原始低通滤波结果进行对比，并计算三组信号的信噪比和均方误差。如图11所示。起始阶段受限于数据长度，信号的频率特性没有完整地反映出来，EMD-WCF保留信号中部分噪声未被完全滤除，导致前期精度稍低；而传统EMD舍弃了所有高频量，在前期获得了更平滑和更准确的结果。整体来看，传统EMD方法处理后的误差大于所提EMD-WCF方法处理后的结果。通过计算，EMD-WCF处理后的信号信噪比为8.60 dB，传统EMD处理后的信号信噪比为5.97 dB，模拟的原始低通滤波信号信噪比为0.24 dB，EMD和EMD-WCF处理结果的信噪比相比于原信号均有较大提升，而EMD-WCF处理后的信号均方误差更小、信噪比更高。
图11EMD-WCF方法对比原始滤波结果和传统EMD处理后的误差
Fig.11The EMD-WCF method compares the error between the original filtering result and the traditional EMD processing
因此，EMD-WCF算法对重力信号降噪具有积极作用，能够在不失分辨率的情况下提高测量精度，同时也防止了重力数据受到低频规律噪声的污染。
3 结论
本文针对当前动态重力测量数据处理时难以平衡测量结果的精度与分辨率问题，提出基于经验模态分解与波数相关滤波的高频重力降噪方法，对低通滤波后的信号进行进一步处理，降低滤波器边界效应等原因造成的高频信号耦合噪声的影响。仿真结果表明，EMD-WCF降噪方法相比于传统EMD处理信噪比提高了44%、均方误差减小了48%，证明了所提方法的有效性，为真实重力数据特别是矢量重力的进一步处理提供了新的思路；同时，基于分频降噪的误差处理方法也为后续开展低频重力误差补偿研究打下基础。
