基于自寻优小波降噪算法的海洋重力数据滤波

李冬毅，覃方君，黄春福，李 安

（海军工程大学 电气工程学院，武汉 430033）

海洋重力场信息是地球固有的物理特性之一，实现海洋重力场的精密测量，对地球物理研究、地震火山监测、惯性导航等都具有非常重要的意义[1]，船载重力测量是测量沿海和近海地区海洋重力场的基本方法，与其他技术相比，它在沿海地区具有能实现高精度和高分辨率测量的特殊优势[2]。船载重力测量过程会受到温度、湿度、海浪、季风等因素的影响，测得的原始信号会包含复杂的噪声信号。文献[3]指出，动态条件下干扰信号的幅值可能比重力异常信息高出百倍甚至几千倍。为了消弱干扰噪声的影响，必须对重力传感器原始观测数据进行滤波处理。国内外目前已开展了大量的滤波方法研究，主要包括数字滤波器[4]、Kalman 滤波[5]、小波降噪[6]等。

为了解决干扰噪声与有用信号的频谱混叠问题，数字滤波器设计时需要具有良好的窄带低通以及尖锐的高频衰减特性[3]。为了满足这一条件，通常需要结合具体测量数据有针对性地调整参数，过程比较繁琐，导致数字滤波器工程实用性不强；另外，数字滤波器处理结果存在较大的边沿效应和相位延迟，严重影响了有效数据的个数[7]。Kalman 滤波是一种基于时域模型的最优估计方法，能有效地进行重力测量原始数据降噪，但其结果依赖于精确的系统模型和测量噪声模型，在重力测量领域的处理效果还有待进一步提升[7]。小波降噪是对小波分析的典型应用，与数字滤波器相比，它具有多分辨分析的特点，能够有效地区分信号中的突变和噪声，在非平稳信号降噪方面有着良好的性能[8]。刘万国模拟了重力异常数据，利用Daubechies 小波函数分析发现，由于重力信息信噪比非常小，阈值去噪法取得的降噪效果与强制去噪法差别不大[9]。赵立业将其应用于海洋重力测量数据处理中，取得了优于数字滤波器和Kalman滤波方法的滤波效果[10]。然而，小波降噪本质上是带通滤波，无法对指定频段外的噪声进行很好的处理[11]，滤波效果受噪声复杂程度影响较大，需要开展相关的改进研究。

针对船测重力数据的滤波需求，对传统小波降噪方法的改进主要有两种方案，一种是基于小波分析原理，拓展更多的滤波方法，找出更适合用于滤波的方法[12,13]；另一种方案是将小波降噪理论与其他算法结合，通过其他算法弥补小波降噪理论的固有缺陷，得到更优的滤波方法[14,15]。显然，第二种方案可以有针对性地提升传统小波降噪方法的滤波效果，进一步改善滤波质量。遗传算法起源于对生物遗传和进化开展的模拟研究，是一种基于生物种群发展机制的自适应概率优化技术。该算法直接处理的对象不是特定的参数，而是根据一定规律编码得到的基因个体，因此没有算法模型限制，可行解的表示形式也非常广泛[16]。利用遗传算法的全局搜索能力，可以有效拓展小波分析方法的适用性，具有广泛的应用前景。肖渝在开展水质光谱数据的去噪研究中，选用信噪比作为适应度函数，利用遗传算法对小波降噪的阈值开展全局寻优，解决了传统小波去噪方法在该领域存在的信号失真问题，较好地保留了光谱的完整性[17]。然而，根据定义可知，信噪比需要在有用信号已知的情况下才能计算，但在重力数据处理过程中，有用信号是希望得到的处理结果，事先并不确定，因此该改进方法不能直接应用于重力数据处理，需要对适应度函数进行调整。

针对船载重力信号的缓变特性，本文将标准差作为适应度函数，利用遗传算法搜索小波降噪中阈值参数的全局最优解，得到自寻优小波降噪算法，并将新算法与常用滤波算法进行对比，验证新算法的性能与优越性。

1 传统小波降噪方法的局限

传统小波降噪方法的主要步骤如下[18]：

1)小波分解。选定合适的小波及需要分解的层次，对原始信号进行分解计算，离散小波变换会将信号分解成低通子带A 和高通子带D，如图1 所示。图中↓2 表示样点数减少一半。经过下采样后得到的Wn为细节系数，Vn为近似系数。

图1 小波分解示意图Fig.1 Schematic diagram of wavelet decomposition

2)阈值量化。对分解得到的一系列细节系数和近似系数利用选择的阈值进行量化处理。

3)小波重构。重构是分解的逆过程，通过上采样和滤波由细节系数和近似系数恢复出去噪后的信号，如图2 所示。

图2 小波重构示意图Fig.2 Schematic diagram of wavelet reconstruction

上述步骤中，决定小波降噪滤波效果的要素主要包括小波种类、小波分解层数、阈值的量化方法[19]。关于船测重力数据滤波过程中，如何选择小波降噪各项要素的问题，已有许多学者开展了相关的研究，并形成了一定的成果。关于小波种类、小波分解层数的选择，韦建成采用捷联式海洋重力测量系统的实测数据，对多种小波种类、多种分解层数的小波降噪滤波效果开展了研究，试验结果表明，为处理海洋重力测量数据，可以选用的小波种类包括db6、db7、db8、db9、db10、sym6、sym7、sym8、sym9、sym10、coif3、coif4、coif5，可以选用的小波分解层数为8、9[20]。关于阈值量化方法的选择，除了阈值计算方法外，其他因素在原理上都有比较明显的优劣之分，在实际重力数据处理的过程中也得到了验证。例如，阈值去噪与强制去噪相比，在区分有用信号和噪声信号时设置了一定的门限，能够更大限度地保留有用信号的细节，可以弱化去噪结果的边界效应，减少测线端部的无效数据量[20]；软阈值与硬阈值相比，具有更好的数学特性，在二维重力数据滤波处理过程中效果更优[21]。

综合上述研究内容，可以看出，在船测重力数据滤波过程中应用的小波降噪方法，其小波种类、分解层数及阈值量化方法中的大多数因素都已经有了相对成熟的选择指标，合适的阈值计算方法是影响小波降噪效果进一步提升的关键因素。常用的阈值计算方法如表1 所示。

表1 常用阈值计算方法Tab.1 Common threshold calculation methods

将基于这四种阈值计算方法的小波降噪应用于实测重力数据的处理中，得到的结果如图3 所示，统计滤波后信号的标准差、信噪比及均方根误差如表2 所示。

表2 基于不同阈值的小波降噪效果分析Tab.2 Effect analysis of wavelet denoising based on different thresholds

图3 基于不同阈值的小波降噪效果Fig.3 Wavelet denoising effect based on different thresholds

从上述图表可以看出，针对重力测量数据处理领域，应用常用阈值计算方法设计的小波降噪可以取得一定的滤波效果，但仍存在波形畸变，无法达到最优的滤波效果，且针对更为复杂的噪声环境，有可能存在无法滤除噪声的情况。为了适应高精度海洋重力测量的需求，有必要结合海洋重力测量的特点，参考常用阈值计算方法，研究出更适用于海洋重力数据处理的小波降噪方法。

2 自寻优小波降噪算法

由上一节可知，常用阈值计算方法无法适应复杂多变的信号特征，尽管采用自适应阈值可以解决不同尺度下噪声特性的差异问题，但仍无法确保计算出的阈值能够达到最优滤波效果，需要结合待处理的数据特点来优化阈值计算方法。针对此问题，考虑到船载重力测量过程中航速低，航行区域有限，测得的重力数据中的有用信号大多属于缓变信息的特点，本文提出了自寻优小波降噪算法，算法流程图如图4 所示。

图4 自寻优小波降噪算法流程图Fig.4 Self-optimizing wavelet denoising algorithm flowchart

自寻优小波降噪算法借鉴了遗传算法的思想，对传统小波降噪方法的基本步骤进行了重新组合，具有能自动修正阈值参数的优点。在传统小波降噪方法的基础上，计算出常用阈值并代入到遗传算法中，编码后生成初始种群，以小波重构后信号的标准差作为适应度函数确定出最优个体，在保留最优个体的条件下对种群进行选择、交叉、变异操作形成新种群，如此往复，直至满足终止条件，译码并输出最优阈值，利用最优阈值进行阈值量化和小波重构即可得到滤波后的有用信号。

适应度函数的定义式如下

式(1)表征了小波重构后信号的标准差，标准差是表示样本中所有个体变异大小的指标，可以反映出整个样本相对于样本平均数的离散程度。实际重力信息与地球质量分布有关，变化极为缓慢，且船载重力测量过程时间有限，系泊状态下的有用信号可以看作常值，航行状态下的测得的有用信号也受制于低船速、短测线的影响，通常变化幅度也比较小；另外，船测重力数据的噪声幅值远大于有用信号变化幅值，理论上通过滤波后信号的离散程度来评价滤波质量是完全可行的。因此，选用标准差作为适应度函数能够满足海洋重力数据滤波的需求，适应度函数越小，表明重构后信号离散程度小，对应的阈值越好。

3 试验分析

为试验自寻优小波降噪算法的性能，本文选取常用的IIR 滤波器、FIR 滤波器、Kalman 滤波算法、传统小波降噪方法作为参考，比较不同降噪方法的处理效果。其中，IIR 滤波器选用4 阶Bessel 滤波器，FIR 滤波器选用切比雪夫最佳一致逼近法进行设计。

3.1 评价指标

本文选用滤波后信号的标准差、信噪比、均方根误差作为评价指标来衡量不同方法的滤波质量。标准差越小、信噪比越大，有用信号与滤波后信号的均方根误差越小，则降噪效果越好。记有用信号为x，降噪后的信号为，的均值为μ，信号数据长度为N，则标准差(STD)定义如下：

信噪比(SNR)定义如下

均方根误差(RMSE)定义如下

3.2 仿真数据处理

船载重力测量时，重力仪会同时受到环境因素和测量船运动带来的干扰，使得测得的重力数据中同时包括有用信息、背景噪声及外部干扰。相对于有用信号而言，“干扰”和“噪声”均属于高频信号。因此，本文将外部干扰与背景噪声统一视作“噪声”，用滤波方法进行处理。

仿真噪声信号幅值设置在有用信号的10 倍以上，利用缓变的正弦波模拟待分离重力信息，即有用信号，用高频正弦波和满足正态分布的随机信号模拟噪声信号，得到原始信号如图5 所示。

图5 仿真原始信号Fig.5 Simulated original signal

利用常用滤波方法和自寻优小波降噪算法对原始信号进行处理，得到逼近误差如图6 所示，滤波后参数统计如表3 所示。

表3 仿真噪声条件下不同降噪方法分析Tab.3 Analysis of different denoising methods under simulated noise

图6 不同降噪方法对仿真信号的滤波效果Fig.6 Filtering effects of different denoising methods on simulated signals

通过对仿真数据的处理可以看出，这几种滤波方法都能抑制噪声，滤波效果从差到好依次为：IIR 滤波器、FIR 滤波器、Kalman 滤波、传统小波降噪、自寻优小波降噪。从图6 可以看出，数字滤波器滤波后结果存在较大的边沿效应，与理论分析一致，需要舍弃两端的部分数据才能达到较好的滤波效果；Kalman 滤波器的滤波效果较好，但不能很好地复现缓变有用信号的变化趋势，滤波结果波动较大；自寻优小波降噪算法保留了小波降噪本身的优势，与数字滤波器、Kalman 滤波算法相比，滤波结果也更为平滑，更贴合实际重力数据的变化规律。自寻优小波降噪算法相比传统小波降噪算法，滤波后信号精度提升31.9%。

3.3 系泊状态重力测量滤波

2021 年3 月，本团队联合中科院精测院、哈尔滨工业大学完成了基于原子重力仪的湖上重力测量试验。试验中选用的重力仪主要有三种：提供起始测点静态重力基准的CG-5 型重力仪、提供动态重力基准的dgShip型重力仪、被试的冷原子重力仪。CG-5 型高精度静态相对重力仪观测精度优于5μGal，静态长期漂移小于1.0 mGal/天。dgShip 型相对重力仪动态重复测量精度优于0.3 mGal。被试的冷原子重力仪测量灵敏度可达，长期稳定性为5.5μGal。

接下来将利用该次试验数据验证不同降噪方法的滤波效果。为了凸显滤波方法的效果，本文所选用的试验数据已对相对重力仪、原子重力仪的测量数据进行了厄特弗斯改正。

选取系泊状态下测得的原子重力仪数据作为原始信号，如图7 所示。在动态条件下，原子重力仪噪声强度更大，主要原因是在动态条件下载体运动会导致原子重力仪中拉曼光反射镜的振动，从而加大原子干涉的相位噪声。由于此次试验缺少系泊状态下相对重力仪测量数据，本文选取航行状态下的初始段相对重力仪测量结果的平均值作为重力基准，所选数据测定时测量船并未启动，且与原始信号测点位置相同，依据该平均值扩充数据长度，得到有用信号。

图7 系泊状态下重力测量原始信号Fig.7 Original signals of gravity measurement in moored state

利用常用滤波方法和自寻优小波降噪算法对原始信号进行处理，得到逼近误差如图8 所示，滤波后参数统计如表4 所示。

表4 系泊状态下不同降噪方法分析Tab.4 Analysis of different denoising methods in mooring state

图8 系泊状态下不同降噪方法的滤波效果Fig.8 Filtering effects of different denoising methods in moored state

通过对系泊状态重力测量数据的处理结果分析可以看出，几种方法的滤波效果从差到好依次为：传统小波降噪、Kalman 滤波、IIR 滤波器、FIR 滤波器、自寻优小波降噪。尽管传统小波降噪方法在处理仿真数据时取得了较好的滤波效果，但针对具有更复杂噪声特性的原始信号，其滤波效果的稳定性难以得到保证，甚至可能无法滤除噪声，而自寻优小波降噪算法仍能有效地滤除噪声，且滤波效果优于常用滤波方法，滤波后信号精度至少提升50%。

3.4 航行状态重力测量滤波

选取相对重力仪测得的重力数据作为有用信号，将相同时段、相同航迹下测得的原子重力仪数据作为原始信号，如图9 所示。

图9 航行状态下重力测量原始信号Fig.9 Gravity measurement original signal in sailing state

利用常用滤波方法和自寻优小波降噪算法对原始信号进行处理，得到逼近误差如图10 所示，滤波后参数统计如表5 所示。

表5 航行状态下不同降噪方法分析Tab.5 Analysis of different denoising methods in sailing state

图10 航行状态下不同降噪方法的滤波效果Fig.10 Filtering effect of different denoising methods in sailing state

针对航行状态下重力测量数据的滤波处理，几种方法的滤波效果从差到好依次为：FIR 滤波器、IIR 滤波器、Kalman 滤波、传统小波降噪、自寻优小波降噪。从图10 可以看出，针对缓变重力信号的滤波处理，数字滤波器不可避免地会出现边沿效应，需要对数据进行截短处理，影响有效数据的长度；Kalman 滤波方法的滤波结果存在频繁突变，不能很好地重构缓变信号；传统小波降噪方法能够在一定程度上恢复信号，但仍存在部分噪声无法完全滤除。与上述方法相比，自寻优小波降噪算法滤波后信号精度至少提升17.6%，逼近误差最小，且能更大程度地复现缓变重力信号的变化规律。

3.5 算法处理速度分析

面对海量重力数据处理时，算法的处理速度、功耗是一个必须面对的问题。目前重力数据处理方法均为事后处理，一般基于台式计算机平台进行运算，因此功耗问题可转化为处理方法运行时间问题。

本文试验采用的硬件平台为RedmiBook Pro15 计算机，CPU 型号为 Intel(R) Core(TM) i5-11300H，采用220V 交流电供电。软件平台为Windows10 操作系统和Matlab R2021a。将3.2 至3.4 节中不同降噪方法的处理时间进行统计，如表6 所示。

表6 不同降噪方法的处理时间Tab.6 Data processing time of different noise reduction methods

从表中可以看出，相比于常用滤波方法，自寻优小波降噪算法的数据处理时间较长，这是由于自寻优小波降噪算法的复杂性导致的。自寻优小波降噪算法需要将遗传算法的全局搜索功能融合到传统小波降噪算法中，不可避免地会涉及到多子代、多种群中不同个体的适应度函数计算，计算量会大幅增加，因此信号处理速度会远低过传统小波降噪算法，与其他滤波算法相比也没有显著优势。但目前船载原子重力测量领域需要解决的问题主要是测量噪声导致的数据精度不高，可以考虑牺牲一定的实时性来换取更高的精度，因此本文所提的自寻优小波降噪算法仍有较大的应用价值。

4 结论与展望

本文围绕船载重力测量数据滤波需求，依据重力测量数据特点，提出了自寻优小波降噪算法，主要结论如下：

1) 常用重力测量滤波算法都存在一定的缺陷，不能适应含有复杂噪声信号的消噪需求。数字滤波器的滤波结果存在边沿效应，需要对数据长度进行截短才能得到较好的滤波效果；Kalman 滤波处理后的数据存在频繁突变，不适用于提取缓变重力信息；传统小波降噪方法的滤波效果受噪声信号复杂程度影响较大，泛化能力差。

2) 自寻优小波降噪算法摆脱了传统小波降噪方法的局限，在具有复杂噪声特性的原始信号消噪中具有显著优势，其滤波效果优于传统小波降噪、Kalman 滤波方法、FIR 滤波器、IIR 滤波器等常用滤波方法，为从含噪海洋重力测量数据中提取高精度重力信息提供了可能。

目前，本文所提方法的有效性仅在重力数据后处理的过程中进行了验证。后续应当增加更为丰富的实测试验，可以考虑将自寻优小波降噪算法应用在实际重力测量过程中，进行实时处理，并与其他数据处理方法对比内存占用量、信息处理速度、功耗等参数，从多角度去评估新方法对重力测量过程的影响。此外，本文所提方法仅适用于重力缓变区域数据，针对重力变化较大的海沟等区域数据还需要探索其它降噪方法。