一种基于集成框架的微地震数据去噪方法

技术领域

本发明属于微地震数据处理领域，尤其涉及一种基于集成框架的微地震数据去噪方法。

背景技术

在油气田开发阶段，微地震接收地震压裂产生的振动，由此判断地下结构、实时监测压裂情况，及时调整压裂参数，从而降低开发成本。由于震源能量低，传播路径复杂，通常微地震采集得到的数据中包含大量噪声，包括呈随机分布的蓝色噪声、谐波噪声等。国内外大量学者对随机噪声的压制方法进行大量研究，基本包括三大类：单道信号的噪声压制方法、具有多道结构特征的随机噪声压制算法和深度学习类算法。单道信号的噪声压制方法（Huang et al.， 2020）需要考虑信号自身的特征（如频率、振幅等），如，基于滤波器的双边滤波去噪（Zhang and Guo， 2017）、局部投影去噪（Chen et al.， 2018）、基于递归滤波器的随机噪声压制算法（Mousavi and Langston， 2016）等，也包括利用自适应信号频率分解算法进行噪声衰减，包括：EMD（Empirical Mode Decomposition， Gómez and Velis，2016; Li et al.， 2020）、EEMD（Ensemble Empirical Mode Decomposition， Han andBaan， 2015; Jia et al.， 2016; Zuo et al.， 2019）和VMD（Variational modedecomposition， Chen et al.， 2021）等以及时频域去噪算法（Vera Rodriguez et al.，2012;Mousavi et al.， 2016; Akram et al.， 2016）。具有多道结构特征的随机噪声压制算法考虑数据的横向连续性，选择合适的稀疏变换，在变换域中完成噪声衰减，包括Radon变换（Sabbione et al.， 2015）、 变换（Forghani-Arani et al.， 2013）、Shearlet变换（Liang et al.， 2018）等。深度学习类包括：自编码器卷积神经网络（Dong et al.，2019）、自适应稀疏字典学习及去噪（Wang et al.， 2020）等。总的来说，这些算法基本针对数据的随机噪声，并且大部分根据数据的可预测性和频率上的差异完成降噪。

而在图像处理领域，以非局部块相似性为正则化算法，可完成数据的去噪、修补、去模糊等操作，由于叠后地震数据的非局部块相似性特征，众多学者将BM3D (Dabov etal. 2007)、NLM (Non-local mean， Buades et al. 2008)算法引入地震数据中，同压缩感知（Donoho， 2006）结合，实现地震数据的去噪（Buades et al.， 2005; Wang et al.，2019）、重建（Chen et al.， 2021）、高分辨率（Chen et al.， 2022）等工作。地震信号作为反射系数的低通滤波，其波形存在大量的相似，这也是反演多解性的原因。因此地下地层结构不能完全相似，就像不可能存在两片完全相同的树叶。因此，基于块匹配的非局部相似性或多或少会存在一些问题。除随机噪声外，微地震中还可能含有振幅不稳定的规则谐波噪声，目前少有文献对微地震这种噪声的产生机理深入研究。

发明内容

为解决上述问题，本发明提供了如下方案：一种基于集成框架的微地震数据去噪方法，包括：

基于周期性的谐波噪声，将有效信号视为非周期，通过奇异谱分析，提取数据的周期状谐波信号，获得去除谐波的信号；

基于所述去除谐波的信号，根据频率的差异，在时频域划分有效信号和随机噪声；然后利用非局部波形相似性，计算具有波形相似结构的信号均值，实现微地震数据中的随机噪声压制。

优选地，通过奇异谱分析，提取数据的周期状谐波信号，获得去除谐波的信号的过程包括，基于构造在时间序列上的特定矩阵进行奇异值分解，从所述时间序列中分解获得趋势、振荡分量和周期规则的谐波噪声。

优选地，基于构造在时间序列上的特定矩阵进行奇异值分解，从所述时间序列中分解获得趋势、振荡分量和周期规则的谐波噪声的过程包括，

对于长度为

基于所述向量组建Hankel矩阵，对所述Hankel矩阵进行奇异值分解，然后进行分组操作，根据下标集合划分获得互不相交的子集；其中，所述子集代表不同的趋势成分；

根据信号特征选择部分子集参与重建，进行对角线平均，得到分解后的序列。

优选地，所述向量的表达式为：

所述Hankel矩阵的表达式为：

。

优选地，对所述Hankel矩阵进行奇异值分解，然后进行分组操作的过程包括，

令

令

；其中，/>

然后进行分组操作，将下标集合

优选地，所述根据信号特征选择部分子集参与重建，进行对角线平均，得到分解后的序列的过程包括，

令

将每个矩阵

令

优选地，所述对角线平均公式的表达式为：

。

优选地，基于所述去除谐波的信号，根据频率的差异，在时频域划分有效信号和随机噪声的过程包括，

使用所有尺度计算小波变换系数的叠加振幅，利用小波变换将数据转换到时频域后，所述数据在低频时间分辨率低，频率分辨率高，高频时间分辨率高，频率分辨率低；

基于在数据时频系数中的部分频带内，微地震数据的随机噪声频带较低，且分布连续，能量均衡；有效信号呈间断式能量分布，分布在较高频带的特性将高频从微地震事件的频率范围中分离出来；

所述叠加振幅的计算公式为：

其中，

优选地，利用非局部波形相似性，计算具有波形相似结构的信号均值的过程包括，

根据信号的非局部均值算法，单道含加性噪声的微地震数据

对于给定的位置

其中，

式中，

在纯噪声位置计算自适应参数

其中，

优选地，所述微地震数据去噪方法还包括基于相对误差和余弦相似度对去噪能力进行检测；

其中，所述相对误差表示去噪后数据和干净数据的差异，相对误差越小，去噪效果越好；所述相对误差定义为：

所述余弦相似度表示去噪后数据和干净数据的相似程度，相似度越接近于100%，表示去噪后数据和干净数据越相似，定义为：

。

与现有技术相比，本发明具有如下优点和技术效果：

本发明将奇异谱分析方法与非局部波形相似性原理结合，这两种方法的联合使用可以压制微地震数据的大部分噪声，具有较大实际应用价值。

本发明利用奇异谱分析算法压制数据中的谐波噪声，基于非局部波形相似性提取有效信号。通过奇异谱分析算法考虑了数据中的谐波影响，并且充分利用有效信号和随机噪声的频率差异及两者的波形相似性差异完成降噪。

附图说明

构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1为本发明实施例的基于集成框架的微地震数据去噪处理流程图；

图2为本发明实施例的时频域分割数据示意图；

图3为本发明实施例的非局部波形相似性的计算原理图；

图4为本发明实施例的合成地震记录及加噪后数据图；

图5为本发明实施例的算法去噪后数据及误差图；

其中，（a）为算法去噪后的数据图；

（b）为去噪后数据和无噪数据的差图；

图6为本发明实施例的去噪前后时间域对比图；、

其中，（a）表示TF-Thre方法去噪前的时间域对比图；

（b）表示TF-Thre方法去噪后的时间域对比图；

（c）表示SSA-TF-Thre方法去噪前的时间域对比图；

（d）表示SSA-TF-Thre方法去噪后的时间域对比图；

（e）表示VMD方法去噪前的时间域对比图；

（f）表示VMD方法去噪后的时间域对比图；

（g）表示SSA-NLWM方法去噪前的时间域对比图；

（h）表示SSA-NLWM方法去噪后的时间域对比图；

图7为本发明实施例的算法鲁棒性测试图；

图8为本发明实施例的算法抗噪性测试图；

图9为本发明实施例的TF-Thre算法和NLWM算法去噪结果图；

其中，（a）表示TF-Thre算法去噪结果图；

（b）表示NLWM算法去噪结果图；

图10为本发明实施例的含噪多道压裂数据图；

图11为本发明实施例的微震数据去噪后剖面图；

其中，（a）表示采用TF-Thre算法的微震数据去噪后剖面图；

（b）表示采用VMD算法的微震数据去噪后剖面图；

（c）表示采用SSA-NLWM算法的微震数据去噪后剖面图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。

需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

如图1所示，本发明所提供的一种基于集成框架的微地震数据去噪方法，包括，奇异谱分析（SSA）和非局部波形均值（NLWM）。对于周期性的谐波噪声，将有效信号视为非周期，通过奇异谱分析，提取数据的周期状谐波信号。此外，对于去除谐波的信号，将只包含随机噪声和有效信号；然后根据频率的差异，首先在时频域大致划分信号和随机噪声，利用非局部波形相似性，计算具有波形相似结构的信号均值，达到压制随机噪声的目的。这两种方法的串联可以压制绝大部分微地震数据的噪声。

利用非局部均值算法压制随机噪声会忽略频率域特征。为此，本实施例通过时频域划分方案分离部分噪声和有效信号。最后，由于非局部波形的相似度寻找相似的信号结构，而对周期性较强的谐波无效，本实施例遵循“先规则后随机”的去噪规律，引入奇异谱分析压制谐波。

进一步地，本实施例的奇异谱分析是一种处理非线性时间序列数据的方法，可以对时间序列进行分析和预测。基于构造在时间序列上的特定矩阵的奇异值分解（SVD），从一个时间序列中分解出趋势、振荡分量和噪声。SSA具有非常广泛的适用性，对于时间序列，既不需要假设参数模型，也不需要假设平稳性条件。

具体地，对于长度为

这些向量组成Hankel矩阵：

对Hankel矩阵进行奇异值分解，令

令

；其中，/>

然后进行分组操作，将下标集合

为了重建信号，需要进行对角线平均，将上式的每个矩阵

通过上述方法即可实现含谐波微地震数据的去噪。

进一步地，针对时频域信号和噪声的划分，本实施例基于时频域信号噪声分割策略进行信号的预处理，算法假设：信号的能量集中在时频系数的高频部分，而噪声则由值较小，且频率较低的系数表示，且分布面积广。算法根据时频系数的峰值分布，大致区分有效信号和随机噪声的分布频率范围。

利用小波变换（Continuous wavelet transformation，CWT）将数据转换到时频域后，其低频时间分辨率低，频率分辨率高；高频时间分辨率高，频率分辨率低。在数据时频系数中的部分频带内，微震数据的随机噪声频带较低（较高尺度），且分布连续，能量均衡；而有效信号呈间断式能量分布，分布在较高频带（较低尺度）。利用这一特性将高频从微地震事件的频率范围中分离出来。为此，使用所有尺度计算CWT系数的叠加振幅：

其中，

进一步地，根据信号的非局部均值算法，单道含加性噪声的微地震数据

对于给定的位置

其中，

式中，

其中，

本实施例的非局部波形均值的权重

测试一：

进一步地，本实施例的微地震数据去噪方法还包括基于相对误差和余弦相似度对去噪能力进行检测；

其中，所述相对误差表示去噪后数据和干净数据的差异，相对误差越小，去噪效果越好；所述相对误差定义为：

所述余弦相似度表示去噪后数据和干净数据的相似程度，相似度越接近于100%，表示去噪后数据和干净数据越相似，定义为：

首先使用一个简单模型测试非局部波形均值（NLWM）的去噪能力，基于随机稀疏反射系数同Ricker子波褶积合成地震记录如图4所示，其中，线条①为添加全频带随机噪声后的数据，线条②为干净数据，在去噪测试中，选择的对比算法为：小波阈值去噪（WTD）、经验小波变换（EWT）和变分模态分解（VMD）算法。

测试结果如图5所示，可以发现NLWM和VMD的去噪结果类似，去噪后信号和无噪信号的相似度差距较小（图5（a），NLWM和VMD的相似度分别为97%和96%），且去噪后信号和无噪信号的差剖面上也证明了这个结论（图5a，NLWM和VMD的相对误差分别为23%和27%）。总体上，NLWM和VMD去噪效果优于WTD和EWT算法。但这种数值上的微弱优势并不能证明NLWM在随机噪声压制方面的优势。但是需注意，在如箭头所示的位置，无波形显示，NLWM在这部分的去噪能力较强，而微震信号绝大部分为稀疏分布，这使得NLWM在微震数据的随机噪声压制中具有得天独厚的优势。

进一步地，使用稀疏模型继续测试算法的去噪能力，本模型测试部分通过合成采样间隔为1ms的地震数据，有效信号使用主频50Hz的Ricker子波和随机反射系数褶积形成（图6（a）、6（c）、6（e）和6（g）黑色实线），同时加入正弦波噪声和蓝色噪声模拟微震数据中的谐波和随机噪声。合成模型如图6（b）、6（d）、6（f）和6（h）黑色实线所示。各算法的临时简称及算法来源如表1所示。

表1

对比方案选择时频域随机噪声阈值压制算法（下称TF-Thre）、时频域随机噪声压制算法同SSA的结合（下称SSA-TF-Thre），同时具备压制谐波和随机噪声能力的变分模态分解算法（VMD）和我们建议的SSA同非局部波形均值的结合（SSA-NLWM），上述四种方法的预期效果及算法来源见表1。去噪结果如图6所示，建议的SSA-NLWM表现出良好的去噪性能，其它三种方法存在不同程度的残余噪声（图6（g）和6（h））；同时，比较图6（a）和图6（c）可以发现，SSA具有相当可靠压制周期噪声能力，具体的，其算法性能在Sacchi（2009）中有详细说明，这体现了算法在本模型的适用性。

为了测试算法的鲁棒性和抗噪性，本实施例取消了随机数种子的设置，这使得每次执行算法生成的噪声、随机反射系数各不相同，基于该模型，使用TF-Thre、SSA-TF-Thre和SSA-NLWM进行算法鲁棒性测试（VMD依赖分解的参数，需要人工选择固有模态，未参与大量测试）。共测试1000次，相对误差曲线（高斯平滑后）如图7所示，建议的SSA-NLWM误差一直最小（线条①），同时SSA-TF-Thre（线条②）和SSA-NLWM（线条①）相差“随机噪声”。

然后固定随机数种子，改变信噪比，测试算法的去噪能力，如图8所示，可以发现，随信噪比的降低，几种方法的误差均随之增大，但是建议的SSA-NLWM误差一直最小（线条①）。并且，SSA-TF的稳定性要优于其它两种算法，这是因为算法对谐波无压制作用。

测试二：

基于S. Mostafa Mousavi (2015)文章中应用的开源数据进行测试，如图9所示，建议的NLWM去噪后的赋值和原始信号的匹配度更高，特别对于弱信号（箭头）；并且如图中虚线框所示，在强振幅处，算法也有着准确匹配振幅的优势，去噪效果略优于TF-Thre算法。

另外，使用一个多道的实际压裂监测微震数据进行去噪测试，原始数据如图10所示，数据共26道，其中第1道的信噪比最低，第3-24道有不同程度的谐波噪声干扰，分别应用TF-Thre、VMD和SSA-NLWM进行噪声压制测试，测试结果如图11所示。TF-Thre对谐波能量强的数据处理效果较差，且存在部分残余噪声（如图11（a）第1道）；VMD算法剩余噪声较多，且损伤有效信号（图11（b））。建议的方法虽然也损失部分有效信号，但是仍体现出较高的去噪能力（图11（c））。

以上，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。