一种微震资料去噪方法、系统、电子设备及存储介质

本发明属于微震监测资料处理技术领域，涉及一种微震资料去噪方法、系统、电子设备及存储介质，具体涉及一种基于同步挤压小波变换时频谱图像处理的微震资料去噪方法、系统、电子设备及存储介质。

背景技术：

微震监测技术已广泛应用于非常规油气勘探、矿山安全监测等资源开发中，微震监测资料处理的关键步骤是对微震事件进行精确定位，从而有效地提高微震监测水平。但是，微震数据的可靠性和参数提取的准确性，如微震事件定位所需的到达时间、极性、振幅等属性参数，都直接受到各种噪声的干扰，影响微震事件的定位精度，导致资源开发过程中的不确定性因素增加，降低资源开采水平。

带通滤波是在特定频带范围内对噪声进行去噪的一种常用方法，但对于信号频率相同的噪声去噪效果不佳，且需预先知道有效的滤波频带。在此基础上，其他更复杂、有效的滤波方法也随之被提出，如主成分分析法(principalcomponentanalysis,pca)、奇异值分解法(singularvaluedecomposition,svd)和基于经验模态分解法(empiricalmodedecomposition,emd)等滤波方法；上述滤波方法均是在时间域(一维)进行，可操作性不强，且计算过程较为繁琐。

时频变换(timefrequencytransform，tft)将一维时间域的波形变换到二维时间-频率域内，可以在时频谱上较为清楚的识别出微震信号的位置。常用的时频变换方法有短时傅里叶变换(shorttimefouriertransform,stft)、小波变换(wavelettransform,wt)和s变换(stransform,st)等变换方法。近年来，同步挤压(synchrosqueezed,ss)作为一种新的重新分配技术被引入，其具有更高时频分辨率特征而被广泛应用于时频域滤波，其中时频域中的小波系数阈值方法因原理简单、去噪效果显著而受到广泛关注。但是，常用的硬、软阈值时频域滤波方法分别存在背景噪音滤除不完全、信号失真等问题，而导致滤波效果不理想。

因此，有必要研究合适的去噪方法，提高微震数据的信噪比,从而拾取到更加精确的震相到时，为微震事件精确定位提供可靠的到时数据。

技术实现要素：

针对现有技术中存在的上述问题，本发明提供一种基于同步挤压小波变换时频谱图像处理的微震资料去噪方法、系统、电子设备及存储介质。

本发明公开了一种基于同步挤压小波变换时频谱图像处理的微震资料去噪方法，包括：

根据微震信号的方差比曲线极小值所对应的时刻，估计纯背景噪音范围；

利用同步挤压小波变换将微震记录变换到时频域；

根据所述时频域内各个频点的纯背景噪音小波系数的经验累积分布函数，估计各个频点的小波系数阈值；

根据各个频点的所述小波系数阈值对所述时频域内各个频点的小波系数进行阈值滤波；

计算经小波系数阈值滤波后时频谱图像的连通分量面积以及连通分量面积阈值；

根据所述连通分量面积阈值对所述时频谱图像再次进行阈值滤波；

将经连通分量面积阈值滤波后的时频谱图像进行同步挤压小波反变换，得到去噪后的微震记录。

作为本发明的进一步改进，所述方差比曲线为：

式中，i为微震信号i时刻，s(t)为微震信号t时刻的振幅，n为微震信号长度；

根据微震信号的方差比曲线极小值所对应的时刻t0＝min(rov)，估计得到纯背景噪音范围r＝[1,t0]。

作为本发明的进一步改进，所述估计各个频点的小波系数阈值的计算公式为：

β(a)＝ecdfa^-1(p)

式中，β(a)为频点为a的小波系数阈值，ecdfa为频点为a的纯背景噪音小波系数经验累积分布函数，ecdfa^-1为经验累积分布函数的反函数，p为经验累积分布概率阈值。

作为本发明的进一步改进，对所述时频域内各个频点的小波系数进行阈值滤波的计算公式为：

式中，为频点为a，t时刻滤波后的小波系数阈值，t(a,t)为频点为a，t时刻原始未滤波的小波系数；

若原始未滤波的小波系数大于或等于小波系数阈值，则保留，否则将该原始未滤波的小波系数设置为0。

作为本发明的进一步改进，所述计算经小波系数阈值滤波后时频谱图像的连通分量面积，包括：

对经小波系数阈值滤波后的时频谱图像进行二值像素处理，得到二值像素时频谱图像；

计算所述二值像素时频谱图像中各个连通性分量面积。

作为本发明的进一步改进，所述连通分量面积阈值的计算公式为：

式中，sthre为连通性分量面积阈值，m为连通性分量的数量，sk为第k个连通性分量面积，median表示为中位数。

作为本发明的进一步改进，对所述时频谱图像再次进行阈值滤波的计算公式为：

其中，为连通性分量面积阈值滤波后的时频谱图像，s为时频谱图像的连通性分量面积；

若时频谱图像的连通性分量面积s大于或等于连通性分量面积阈值sthre，则保留，否则将该连通性分量去除。

本发明还公开了一种用于实现上述微震资料去噪方法的系统，包括：

第一估计模块，用于根据微震信号的方差比曲线极小值所对应的时刻，估计纯背景噪音范围；

变换模块，用于利用同步挤压小波变换将微震记录变换到时频域；

第二估计模块，用于根据所述时频域内各个频点的纯背景噪音小波系数的经验累积分布函数，估计各个频点的小波系数阈值；

第一滤波模块，用于根据各个频点的所述小波系数阈值对所述时频域内各个频点的小波系数进行阈值滤波；

计算模块，用于计算经小波系数阈值滤波后时频谱图像的连通分量面积以及连通分量面积阈值；

第二滤波模块，用于基于所述连通分量面积阈值对所述时频谱图像再次进行阈值滤波；

反变换模块，用于将经连通分量面积阈值滤波后的时频谱图像进行同步挤压小波反变换，得到去噪后的微震记录。

本发明还公开了一种电子设备，包括：处理器、通信接口、存储器和通信总线，其中，处理器、通信接口、存储器通过通信总线完成相互间的通信；

所述存储器中存储有计算机程序，当所述程序被所述处理器执行时，使得所述处理器执行上述微震资料去噪方法的步骤。

本发明还公开了一种存储介质，其存储有可由电子设备执行的计算机程序，当所述程序在所述电子设备上运行时，使得所述电子设备执行上述微震资料去噪方法的步骤。

与现有技术相比，本发明的有益效果为：

本发明通过结合时频域内的小波系数和图像连通分量面积阈值方法，实现了对背景噪音更高效的去除，极大地提高了微震记录信噪比，为后续微震事件精确定位提供可靠的初至信息。

附图说明

图1为本发明一种实施例公开的微震资料去噪方法的流程图；

图2为本发明一种实施例公开的微震资料去噪系统的框架图；

图3为本发明一种实施例公开的未加噪音的模拟微震记录和时频谱图像；

图4为对图3的微震记录加入一定程度噪音后的低信噪比微震记录和时频谱图像；

图5为利用方差比曲线极小值估计得到的纯背景噪音范围；

图6为利用经验累积分布函数(ecdf)小波系数阈值滤波后的微震记录及其时频谱图像；

图7为运用时频谱图像的连通性分量面积阈值滤波后的微震记录和时频谱图像；

图8(a)为滤波后的微震记录与图3中未加噪音的模拟微震记录的对比图；

图8(b)为滤除的噪音记录和所对应的时频谱图像；

图9为低信噪比高采样率(6khz)金属矿山微震记录滤波前后的波形和时频谱对比图；

图10为水力压裂(1khz)地面微震监测记录滤波前后的波形和时频谱对比图；

图11为天然地震记录(20hz)滤波前后的波形和时频谱对比图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提供一种基于同步挤压小波变换时频谱图像处理的微震资料去噪方法及系统，其利用信号的方差比曲线极小值估计纯背景噪音范围，为小波系数阈值的精确设置提供较为准确的纯背景噪音范围；并利用更符合纯背景噪音小波系数统计分布特征的经验累积分布函数(ecdf)设置小波系数阈值，将大部分背景噪音小波系数去除；再利用小波系数阈值滤波后的时频谱图像的连通性分量面积阈值对剩余的噪音小波系数进行去除，进一步提高微震信号信噪比，为微震事件精确定位提供可靠的震相信息。

下面结合附图对本发明做进一步的详细描述：

如图1所示，本发明提供一种基于同步挤压小波变换时频谱图像处理的微震资料去噪方法，包括：

步骤1、根据微震信号的方差比曲线极小值所对应的时刻，估计纯背景噪音范围；其中，

方差比曲线为：

式中，i为微震信号i时刻，s(t)为微震信号t时刻的振幅，n为微震信号长度；

根据微震信号的方差比曲线极小值所对应的时刻t0＝min(rov)，估计得到纯背景噪音范围r＝[1,t0]。

步骤2、利用高时频分辨率的同步挤压小波变换方法将一维时间域微震记录变换到二维时频域。

β(a)＝ecdfa^-1(p)

步骤3、根据时频域内各个频点的纯背景噪音小波系数的经验累积分布函数，估计各个频点的小波系数阈值；其中，

估计各个频点的小波系数阈值的计算公式为：

式中，β(a)为频点为a的小波系数阈值，ecdfa为频点为a的纯背景噪音小波系数经验累积分布函数，ecdfa^-1为经验累积分布函数的反函数，p为经验累积分布概率阈值。

步骤4、根据各个频点的小波系数阈值对时频域内各个频点的小波系数进行阈值滤波；其中，

上述计算公式为：

式中，为频点为a，t时刻滤波后的小波系数阈值，t(a,t)为频点为a，t时刻原始未滤波的小波系数；

若原始未滤波的小波系数大于或等于小波系数阈值，则保留t(a,t)，否则将该原始未滤波的小波系数设置为0。

步骤5、计算经小波系数阈值滤波后时频谱图像的连通分量面积以及连通分量面积阈值；其中，

计算经小波系数阈值滤波后时频谱图像的连通分量面积，包括：对经小波系数阈值滤波后的时频谱图像进行二值像素处理，即小波系数绝对值大于零则赋值为1，否则赋值为0；计算二值像素时频谱图像中各个连通性分量面积；

连通分量面积阈值的计算公式为：

式中，sthre为连通性分量面积阈值，m为连通性分量的数量，sk为第k个连通性分量面积，median表示为中位数。

步骤6、根据连通分量面积阈值对时频谱图像再次进行阈值滤波；其中，

上述计算公式为：

其中，为连通性分量面积阈值滤波后的时频谱图像，s为时频谱图像的连通性分量面积；

若时频谱图像的连通性分量面积s大于或等于连通性分量面积阈值sthre，则保留该连通性分量，否则将该连通性分量去除。

步骤7、将经连通分量面积阈值滤波后的时频谱图像进行同步挤压小波反变换，得到去噪后的微震记录。

如图2所示，本发明提供一种用于实现上述微震资料去噪方法的系统，包括：

第一估计模块，用于实现上述步骤1；

变换模块，用于实现上述步骤2；

第二估计模块，用于实现上述步骤3；

第一滤波模块，用于实现上述步骤4；

计算模块，用于实现上述步骤5；

第二滤波模块，用于实现上述步骤6；

反变换模块，用于实现上述步骤7。

本发明提供一种电子设备，包括：处理器、通信接口、存储器和通信总线，其中，处理器、通信接口、存储器通过通信总线完成相互间的通信；存储器中存储有计算机程序，当程序被处理器执行时，使得处理器执行上述微震资料去噪方法的步骤。

该电子设备的技术方案与上述去噪方法的技术方案属于同一构思，电子设备的技术方案未详细描述的细节内容，均可以参见上述去噪方法的技术方案的描述。

电子设备可以是任何类型的静止或移动电子设备，包括移动计算机或移动电子设备(例如，平板计算机、个人数字助理、膝上型计算机、笔记本计算机、上网本等)、移动电话(例如，智能手机)、可佩戴的电子设备(例如，智能手表、智能眼镜等)或其他类型的移动设备，或者诸如台式计算机或pc的静止电子设备；电子设备还可以是移动式或静止式的服务器。

计算机指令包括计算机程序代码，计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。

本发明提供一种存储介质，其存储有可由电子设备执行的计算机程序，当程序在电子设备上运行时，使得电子设备执行上述微震资料去噪方法的步骤。

该存储介质的技术方案与上述去噪方法的技术方案属于同一构思，存储介质的技术方案未详细描述的细节内容，均可以参见上述去噪方法的技术方案的描述。

存储介质可以包括：能够携带计算机程序代码的任何实体或装置、记录介质、u盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(rom，read-onlymemory)、随机存取存储器(ram，randomaccessmemory)、电载波信号、电信信号以及软件分发介质等。

本发明的优点为：

1、提高了微震资料信噪比：针对地面微震资料信噪比低，震相难以识别和初至不易拾取的问题，利用同步挤压小波时频变换方法，得到更高时频分辨率的时频谱。将时频谱看做为一副二维图像，结合小波系数阈值和时频谱图像的连通性分量面积阈值方法，将背景噪音完全去除，得到更加干净的微震信号能量，极大地提高了微震资料信噪比，为震源参数反演提供可信的波形信息；

2、为时频域滤波方法的发展提供了一个新的方向：将一维时间域信号变换到二维时频域，则可以将时频域滤波当作成二维图像处理问题，运用成熟的图像处理方法把所需的图像细节信息保留下来，完成微震信号滤波。

实施例1：模拟数据

先合成未加噪音的模拟微震记录，并利用同步挤压小波变换得到原始微震记录的时频谱如图3所示。

对图3中的微震记录加入一定程度的噪音，使得微震信号初至被噪音淹没，震相信息识别困难；如图4所示，从加噪后的微震记录时频谱图可以看出，背景噪音能量分布于整个时频域，难以利用常规的带通滤波方法将背景噪音完全滤除。

如图5所示，在微震信号到来时刻，方差比曲线有一个较为明显的极小值，从该极小值所对应的时刻可以较为准确估计出纯背景噪声的范围，为后续处理提供合理的纯噪声数据。

如图6所示，背景噪声的小波系数绝大多数被滤除，微震信号几乎完全保留，微震信号的信噪比显著提高。但是，从滤波后的波形和时频谱中可以看到，较少的低振幅、孤立的背景噪声仍然保留了下来，这也是小波系数阈值滤波方法的固有问题。因此，在后续的进一步处理步骤中，利用连通性分量阈值方法将前一步遗留下来的孤立背景噪声的小波系数几乎完全被压制，最终的去噪微震记录信噪比进一步提高，如图7所示。

如图8a、8b所示，去噪后的信号与原始信号在整个波形上非常吻合，初至波的到时、极性等信息十分明显，为震源参数反演提供准确的初至信息。

实施例2：实际资料

金属矿山微震监测由于微震事件定位精度要求较高，往往需要高采样率记录。使得采矿环境下记录的地震信号十分复杂，微震记录中包括微震事件、爆炸和机械噪声等多种信号的混杂，信噪比较低，难以直接提取到所需的微震事件信号。如图9所示，某金属矿山的低信噪比微震记录，初至信息被噪音完全淹没，利用本滤波方法后所得到的微震记录和时频谱，背景噪音几乎完全被滤除，微震信号初至清晰可见，对比滤波前后时频谱可以看出，微震信号能量几乎被完全保留了下来，背景噪音能量被完全去除，表明本方法的滤波效果明显。图10为水力压裂地面微震监测数据滤波前后波形对比，可以看到滤波效果十分明显，对比其中第1、2、9、11和18道微震记录的滤波效果，被噪音完全淹没的微震信号经过本方法滤波后，极大地提高了微震记录信噪比，各道的微震事件信号信息完全显现出来。图11为天然地震信号的滤波效果对比，从原始地震记录的时频谱可以看出，原始地震记录主要被低频噪音所污染，经过本方法滤波后，低频噪音被完全去除，且同时将有效频带内的噪音也去除了，震相信息(p、s波)清晰可见，表明本方法的滤波效果十分明显。

以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。