抑制地下流体观测数据干扰的期望最大化算法∗

周 洋 王 俊

(1.中国地震局地震研究所(地震大地测量重点实验室)，湖北 武汉 430071；2.湖北省地震局，湖北 武汉 430071；3.中南民族大学计算机科学学院，湖北 武汉 430071)

湖北区域地球物理台网目前共计23 个台站83套仪器，其中地下流体仪器41 套，SWY－II 型数字式水位仪9 套、SZW－1A 及SZW－II 数字式温度计15套、RTP－II 型气温气压雨量综合观测仪3 套、WYY－1型气象三要素观测仪14 套[1]。地下流体仪器种类繁多、在网运行数量也是在各观测手段(重力、形变、电磁、流体)中居首，同时受干扰因素众多且难以区分。故寻求一种应用前景更为广阔的数据干扰抑制技术显得尤为重要。

近年来兴起的小波分析方法以其良好的时频分析能力迅速成为非平稳信号处理的有力工具，基于小波分析的去噪方法更是大量涌现，且被证明其具有传统傅里叶变换去噪方法所不能比拟的优越性[2－7]。去除地震信号中随机噪声技术虽然已经有了很大发展，但仍然不能满足人们对地震数据去噪效果的要求，因此改进现有去噪技术，寻求更加精确、有效、快速、便捷的去噪算法，是一项十分有意义、尚待研究的课题[8－12]。目前国内在数字信号处理如降噪方面的技术已日趋成熟，如文献[13－16]介绍了基于时空域数字信号处理的地震资料降噪方法，而文献[17－21]则研究了基于深度学习的自适应学习去噪方法。

本文提出一种基于神经网络的期望最大化算法，用于地震地下流体观测数据噪声抑制。该算法较其他小波算法及神经网络算法去噪在保留原始数据信号能量前提下，更能快速、准确定位并抑制噪声干扰。与常规幅度阈值相比，所得到的概率阈值更适合于数据。它通过概率阈值为用户提供了量化是否将大幅度异常视为噪声的置信度的可能性。该方法显示出比优化的常规方法稍好的性能，且参数测试和变化少得多。

1 期望最大化算法

1.1 简介

地震数据总是会因各种噪声而损坏，从而降低了数据质量。因此，噪声衰减是有助于解释地震数据处理的重要步骤。可以通过基于模型的信号处理来分离信号和噪声。

该方法假定信号和噪声具有可以由特定数学模型捕获的不同特性。该模型可以基于诸如波动理论之类的物理学定律，该定律被用于许多分析方法中。或者，模型可以利用确定性特征，例如不同的迹线偏移或统计属性(例如独立性)以区分信号和噪声。

我们关注高振幅噪声衰减的问题。考虑一个简单的噪声模型，其中所有大于特定阈值的数据样本都可能是噪声。阈值可以以多种形式指定，但通常用阈值因子乘以某些数据统计量(例如，分析窗口内计算的平均值，中位数或均方根值)来表示。阈值因子和数据统计信息是用户定义的参数。通常由于其简单性而调用此模型。但是，其主要缺点是阈值的频繁变化，导致整个数据噪声功率的变化。

我们提出了一种自动阈值确定技术，用于检测选择域中的大振幅噪声。以频率偏移域中的衰减噪声衰减为例，但结果可以推广到其他类型的噪声，例如衍射的多个噪声和尖峰，其中在某个时间或时间偏移中会出现大幅度采样窗口被检测到，然后通过插值删除。我们的技术与常规方法相比，它具有更好的数据适应性。

首先，我们简要回顾一下涌浪噪声的问题。然后，我们概述了自动阈值确定的新方法。最后，我们结合湖北地震台地下流体数据给出一些真实的数据示例。

1.2 降噪衰减

涌浪噪声是由恶劣的天气条件引起的，并且是获取地震数据时经常遇到的问题。它对地震数据质量有不利影响，甚至可能导致暂停采集。它的特点是振幅大，频率低。

讨论了产生涌浪噪声的几种可能机制，得出结论，对于现代的泡沫填充拖缆，最可能的原因是强烈的海浪或拖缆表面上的动态压力变化引起海洋垂直运动，从而引起静水压力波动，这是由于拖缆周围存在湍流层而导致的。

抑制涌浪噪声的常规技术首先计算滑动窗口内数据的幅度或功率谱。窗口内给定频率上超过某个阈值的所有频谱值均被视为噪声。然后对这些噪声样本进行衰减或内插。但是，在此过程中必须定义一个阈值。

对于给定的频率f，我们用rk＝｜Dk(f)｜2定义Sn＝{r1，r2，…，rn}，其中Dk(f)是数据窗口中第k条迹线的傅里叶变换。因此，Sn代表给定频率下的一组迹线的功率谱值。阈值通常计算为

1.3 离群值检测及阈值计算

本节介绍了一种通过检测异常值来确定适当阈值的自动方法。

考虑一个数据集Sn＝{r1，r2，…，rn}，其中样本rk被假定为独立且均匀分布的，由概率密度函数(PDF)g(r)生成。在应用统计中，在集合Sn中查找异常样本称为异常值检测。目的是找到那些与其余数据显示不同统计特性的样本。在我们的问题中，离群值(噪声＋信号)的幅度大于常规数据(仅信号)的幅度。因此，异常值与常规数据的总体不同之处在于，可以通过使用PDF 建模来捕获一些明显的统计量。

便于标记，离群点的PDF 用p(r｜θ1)表示，常规数据的PDF 用p(r｜θ0)表示。为简单起见，假设两个分布都属于同一系列的参数化PDF，但具有不同的参数值，即θ1≠θ0。让标量表示集合Sn中离群值的分数。该参数的统计含义是:从Sn随机抽取的样本ri是离群值的先验概率(即前数据建模)。现在，数据PDF 包含两种模型的混合结构:

可以从Sn估算方程(2)中定义的模型的所有参数。估算完成后，可以使用贝叶斯规则计算给定离群值样本r＝ri的后验概率(即后数据建模)。在贝叶斯规则中，假设事件B已发生，则A和B是两个事件，而P(A｜B)是事件A发生的概率。它表示两个事件的条件似然，例如P(A｜B)P(B)＝P(B｜A)P(A)。这里，r是离群值。

如果离群值ri的后验概率大于概率阈值β，即

例如，对于β＝0.5，式(4)表示仅当模型确定至少有50%是异常值时才选择异常值。式(4)中概率阈值β的选择比式(1)中任意阈值因子α客观得多，并且反映了我们可能需要将任何数据分类为异常值的统计置信度。

计算阈值时，涉及三个步骤:选择一个PDF，识别算法并解释结果。

为了使用式(3)给出的检测标准，需要选择式(2)中p(r｜θ)的参数形式。当r增加时，p(r｜θ)的形式应提供不同的衰减率。需要这样做，以确保与常规数据相比，离群值人群中大振幅值的可能性更大。使用傅里叶变换的实部和虚部是具有等方差的零均值高斯随机变量的近似值，可以证明功率谱样本的PDF 具有指数形式。

指数分布具有一个参数，θ＝λ，它是分布的平均值:

然后，式(2)中的模型变为

用λ1>λ0表示一个事实，即异常值的分布更可能产生较大的振幅。

在使用方程式进行推论之前，人们可能会质疑模型在式(6)中的拟合优度。这个问题的答案很重要，因为它需要对原假设进行统计检验:数据遵循式(6)中的模型。数学模型有助于我们理解和描述现实世界，并且公认的是，真正的数据生成模型比任何可能的假设模型都要复杂得多。因此，模型的正确性问题尚未正式解决，但本文稍后将讨论模型敏感性测试。

式(6)中的模型具有三个未知参数(λ0，λ1，ε)，将从Sn进行估计。我们建议在以下优化问题中使用最大似然估计器(MLE)，以获得其所需的统计属性，例如一致性和效率:

此问题没有封闭形式的解决方案，因此参数使用迭代程序估算值。此过程属于一类称为期望最大化(EM)算法的技术，用于在数据丢失或不完整时查找模型参数的MLE。迭代EM 算法由以下步骤给出:

(2)如下更新参数

在该特定应用中，EM 算法的收敛速度很快，并且在很大程度上与初始条件无关。设置初始条件，使得等于[ε(0)n]个最大振幅的平均值等于其余数据的平均值。

在将未知参数λ0，λ1和ε替换为其估计值的同时，结合式(3)、式(4)和式(6)会得出基于幅度的阈值检测标准。幅度r是概率为β(如果r>rexp)的离群值，其中exp 表示表达式中的指数。

式(12)中给出的阈值考虑了数据的统计信息(通过式(6)中模型参数的估计)和β值的置信度要求。当log(β/(1－β))＝0，β＝0.5，用户对阈值的影响最小。这与以下事实相吻合:对于要做出的任何二元决策，每个概率为0.5 的结果表示信息量不足的情况下，该概率的先验知识不会影响决策。因此，对于建议的概率阈值，β＝0.5 将是更好的默认值。

对式(12)中定义的新阈值标准的分析揭示了几点。首先，当用户增加值β时，接受异常值所需的统计置信度就会增加，从而会增加阈值水平。因此，β具有与常规方法的阈值因子α类似的作用。其次，在统计模型下，较小的表示存在少量异常值，并且增加阈值水平以反映该异常值。另一方面，当时，模型会推断出大多数数据都是异常值，因此阈值水平降低了。最后，当时，方程(6)中两个混合模型的假设没有数据支撑。因此，的值没有统计意义。在这种情况下，我们假定数据仅由常规样本组成，因为排除了数据仅为异常值样本的情况。阈值会自动设置为较大的值，因此不会选择异常值。

为了研究本技术对所选统计模型的敏感性，我们考虑在数据建模中使用瑞利分布而不是指数分布。

一旦通过式(12)中的阈值标准识别出噪声幅度，就可以将它们从Sn中或从相邻样本中删除(例如，通过设置ri＝0)。但是，我们通常用恒定因子对噪声样本进行重新缩放以选择更保守的选项，使得噪声样本的新均值等于指数模型的常规数据的均值，以及噪声样本的新均方根值等于瑞利模型的常规数据的均方根值。

2 实验及结果分析

本文在上述模型算法的基础上，结合湖北省地震地下流体数据实例进行仿真实验，并与传统神经网络及小波去噪算法比较。

2.1 信号去噪

选取湖北省房县三海村台、钟祥马岭台、秭归台水温仪观测数据做仿真实验，时间跨度200 d。如图1～图3。

图1 房县三海村台水温数据

图2 钟祥马岭台水温数据

图3 秭归台水温数据

从图中得到的结果来看，期望最大化算法对噪声有很好的抑制作用，保留信号能量的同时，去噪效果明显，为实际工程应用及地震数据分析预报提供了一个新的研究思路。实现仿真时间段可根据实际需要实时调整，系统可移植性、灵活性强。

2.2 与传统算法的比较

同传统的数据处理方法相比，小波滤波方法具有独特的优势，即能够在去除噪声的同时，很好保留信号的突变部分或图像的边缘。而本文提出的一种新的期望最大化算法则比传统的小波滤波算法更具优越性。

以2020 年5 月荆门台1 测点水温数据为例，给出数据信号在4 种不同的去噪算法处理下的结果，如图4。

由图4 可知本文提出的算法较传统神经网络及小波算法(如模极大值法、空域相关法、小波阈值法)去噪相比，效果更理想。表1 为上述几种算法的检测性能比较，表2 给出几种不同算法的定性比较。

图4 几种去噪方法比较

表1 几种算法检测的结果比较

表2 几种滤波方法的定性比较

从SNR 和RMSE 的定义可知，SNR 值越大，检测效果越好；RMSE 反之。能量成分比越小，检测效果越好，标准残差反之。而局部线性模型法很好地保留了信号发展初期的高频特性，最大限度地反映了原信号本身的性质，且性能参数优于其他几种算法。

3 结论

本文提出一种基于神经网络的期望最大化算法，用于地震地下流体观测数据干扰抑制的新方法。结合湖北地震台站地下流体数据实例做仿真实验，并与其他神经网络及小波分析方法做出比较。得出以下结论:

(1)使用统计建模可以自动检测大幅度的噪声。最终的统计检测方法等效于优化的常规方法，但参数测试和变化少得多。噪声衰减测试显示，在整个数据范围内各种变化的噪声水平下，结果都是一致的。该方法是通用的，可以推广到其他类型的异常高振幅噪声抑制，例如尖峰和衍射倍频噪声抑制。

(2)该算法对地下流体观测数据信号消噪效果好，能有效抑制各种因素干扰，为台站数据分析预报提供一种新的研究方向和方法，为仪器维护人员提供一种新的处理手段和途径。

(3)该算法不依赖于具体的对象的数学模型，可移植性强。该算法可推广应用至其他观测手段的数据干扰抑制，如形变、重力、电磁等。