冷佳宣,喻志超,冯方方,张逸伦,何 川

(1.北京大学地球与空间科学学院,北京100871;2.国家超级计算深圳中心(深圳云计算中心),广东深圳518055;3.中国科学院地质与地球物理研究所,北京100029)

微地震监测技术在非常规油气储层改造中得到广泛应用,通过观测和分析水力压裂诱发微地震事件可以实时监测压裂施工中缝网的发育过程,评价压裂改造效果[1-2]。在水力压裂作业的同时,开展微地震监测数据采集,传感器除了接收到数量有限、发震时刻随机和能量震级较小的微地震信号外,还接收到因压裂施工产生的强能量背景干扰信号及随机噪声[3]。这些干扰信号严重降低了微地震监测数据的质量,影响实时监测过程中低信噪比微地震信号的识别及震相初至的拾取,因此,研究微地震监测记录背景干扰信号压制方法具有非常重要的意义。

受储层岩性、压裂施工方案以及监测条件等因素综合影响,连续压裂监测记录中的有效微地震信号特征及数量并不明确,记录中绝大部分的能量来源于背景干扰信号。为了准确识别微地震事件,在微地震监测数据处理前通常需要对监测记录进行预处理以提高微地震信号的信噪比,常见的预处理方法包括带通滤波、极化滤波[4]和时频域滤波[5-6]。压裂监测现场存在的背景干扰信号持续时间长、能量强且具有时变非平稳的特征,常规的预处理方法难以在去除这些信号的同时最大程度保留有效微地震信号的能量,对微地震监测记录的信噪比提升有限,在一定程度上影响了后续微地震事件识别与初至拾取的效果。因此,在微地震信号特征不完全明确的情况下,背景干扰信号的识别和压制是一种更合理的降低监测记录噪声水平的预处理方式。

本文以井中微地震监测三分量信号之间的同步与相关分析为基础,提出一种基于多维经验模态分解的背景干扰信号识别与压制方法,能够在实际微地震资料处理过程中更好地消除背景干扰对后续处理效果的影响。首先,基于微地震监测记录的表达式提出背景干扰信号压制的基本思路,然后,介绍多维经验模态分解方法的原理,并根据不同阶次本征模态函数分量的能量强弱、偏振特征实现与背景干扰信号对应的能量成分识别及有效去除,最后,以实际资料的应用检验了方法的有效性。

1 方法原理

1.1 微地震监测记录的表达

微地震监测记录可以表示为:

X=AS

(1)

式中:S=[s1,s2,,sK]T表示K个源信号;X=[x1,x2,,xM]T表示M个观测信号;A为源信号与观测信号间的关系矩阵。源信号可以为微地震或者干扰信号,观测信号可以表示不同分量、不同空间测位的记录。

连续监测记录中的微地震信号通常发震时刻随机、数量有限、持续时间短(<1s),能量占比远小于背景干扰信号。图1为实际一次水力压裂施工曲线(油压、排量和砂量)和井下一级检波器垂直分量记录的时频分析结果,其中时频分析采用短时傅里叶变换

图1 压裂施工曲线与连续监测记录垂直分量的时频分析结果

方法,选择1s的时间窗口和50%窗口重叠进行分析,从图中可以明显看出,连续监测记录中包含一些能量较强,持续时间较长的成分,同时其强弱变化与压裂作业阶段相关联,具有时变非平稳的特征。

常规时域或频域滤波的预处理方式在微地震信号特征不完全明确的情况下无法达到最优的处理效果,即欠处理的情况下难以最大程度压制干扰信号,过度滤波则将对有效微地震信号的能量造成损失。图2和图3分别为一段1s监测信号的两个不同带通参数滤波的处理效果,在微地震信号特征不明确情况

图2 带通滤波去噪(50～400Hz)方法的处理效果a 原始记录; b 处理后的记录; c 去除的信号成分

图3 带通滤波去噪方法(100～200Hz)的处理效果a 原始记录; b 处理后的记录; c 去除的信号成分

下,常规预处理过程采用带通滤波的方法,当滤波参数选择过大或者过小时可能无法有效压制噪声且会造成有效微地震信号部分能量损失。从图2c和图3c中可以看出,均存在着噪声压制效果不理想且对微地震有效信号能量(特别是S波)造成损失的问题,这将会对后续的微地震事件震相初至拾取产生很大的影响。

微地震信号预处理的核心在于最大程度提高有效微地震信号与背景干扰信号的信噪比。常规技术侧重于通过提取有效微地震信号的时频域特征及道间相关性等特征,从而最大程度辨识微地震信号的能量成分,提升识别能力。基于多维经验模态分解的背景干扰信号识别与压制方法则试图通过识别与压制连续监测记录中能量较强的背景干扰信号来达到提高有效微地震信号信噪比的目标,因此,提取具有时变非平稳特征的干扰信号是该方法的关键。

1.2 多维经验模态分解

基于经验模态分解(empirical mode decomposition,EMD)的时频分析方法是分析时变和非平稳信号的有力工具,被广泛应用于地震勘探、机械故障诊断以及生物医学信号分析等科学研究和工程应用领域[7-10]。这类方法具有多分辨率和自适应的特点,能够将信号自适应分解成不同尺度的本征模态函数(intrinsic mode function,IMF),与傅里叶分析和小波分析方法相比在非线性非平稳信号分析中具有一定的优势。经验模态分解方法不需要选择基函数,其原理是通过逐步减去时序信号上下包络的平均值得到有限个数的有效本征模态函数分量,这些本征模态函数分量能够描述原信号中不同时间尺度局部特征信息且能够精准重构原始信号[7]。

在实际应用中,抑制或消除算法中存在的模态混叠现象(Mode Mixing)是经验模态分解类方法面临的重要问题,其改进方法如集合经验模态分解和完备经验模态分解等方法的发展在地震信号提取处理中取得了较好的应用效果[11-13]。但是,以往的方法存在两方面的不足,一方面是在分解三分量监测记录时,经验模态分解对每个维度信号分别进行分解的处理方式会出现不同通道尺度排列不确定和同阶尺度频率不一致等问题;另一方面,微地震信号具有“异常事件”的特征,记录中微地震信号的存在会引起间歇现象,使得分解结果中表现出模态混叠,且在微地震信号特征不完全明确的情况下,对有效信号的重构可能存在欠处理或者过处理的现象。

(2)

虽然多维经验模态分解能够实现多通道信号的同步与相关分析,但依旧存在模态混叠的问题。噪声辅助多维经验模态分解算法(noise assisted MEMD,NA-MEMD)通过引入额外的噪声通道信号,利用白噪声表现出的滤波器组特性,可以减小多维本征模态函数中的模态混叠及模式对齐问题对后续信号特征提取的影响[16]。噪声辅助多维经验模态分解算法首先构造一个包含p维原始信号和q维高斯白噪声的(p+q)维复合信号,然后对信号进行多维经验模态分解,在所得的(p+q)维分解结果中删除q维噪声即可得到原p维信号的分解结果[16]。算法步骤为:

为了说明噪声辅助多维经验模态分解算法对含强能量时变非平稳信号的低信噪比微地震监测记录的处理效果,此处构建由三分量的时变非平稳信号U、实际微地震监测记录V以及随机噪声W组成的合成记录S,如公式(3)所示:

S=U+V+W

(3)

图4a为模拟的非平稳调频信号,其各分量的公式为:

(4)

非平稳信号的中心频率f0为35Hz,f0的调频频率f1为5Hz,式中:ρ=0.2为频率变化的幅度。图4b为实际微地震监测记录的三个分量,图4c为包含随机噪声的合成记录,图中随机噪声的振幅方差为0.55,图4d为合成记录各个分量对应的时频分析结果。

图4 三分量合成记录及时频分析a 含随机噪声的时变非平稳信号成分; b 微地震信号成分; c 合成记录; d 三分量对应的时频分析结果

图5至图7分别为经验模态分解、多维经验模态分解和噪声辅助多维经验模态分解对三分量合成记录分解的前7阶本征模态函数分量,图中红线为合成记录中的时变非平稳信号。图5中c3对应时变非平稳信号,c1～c2包含随机噪声和微地震信号;图6中c5对应时变非平稳信号,c1对应高频随机噪声,c2～c4中包含微地震信号和部分随机噪声;图7中c5对应时变非平稳信号,c1～c2对应高频随机噪声,c3～c4中包含微地震信号成分。可以看出,经验模态分解这类方法能够自适应地将信号分解成从高频到低频的若干个信号成分,对比分解的结果可以看出,多维经验模态分解方法所得多维本征模态函数相比于经验模态分解方法同阶尺度频率能够保持一致,可以实现不同通道信号之间的同步与相关分析。此外,噪声辅助的多维经验模态分解方法能够减小多维本征模态函数中的模态混叠现象,更为准确地提取时变非平稳信号。

图5 经验模态分解对合成记录分解的结果

图6 多维经验模态分解对合成记录分解的结果

图7 噪声辅助多维经验模态分解对合成记录分解的结果

1.3 背景干扰信号能量识别与去除

完成经验模态分解之后通常需要对不同尺度信号进行分离和重构,前人方法中有利用原始信号与各本征模态函数之间的互相关系数[17]、求取相邻本征模态函数分量之间的互信息熵[8]和定义自适应间隔阈值[13]等方法来辨识出有效信号,但是这些处理方式在信号特征不完全明确的情况下应用效果并不理想。我们处理的目标信号具有偏振特征,从偏振分析的角度进行识别是最为有效的方式。

由于不同类型的波通常具有各自不同的偏振特征,对于噪声源的位置、激发方式和能量大小的范围等属性稳定的信号,通过极化分析求取波的极化属性可以作为背景干扰信号类型的判别依据。极化分析方法通过三分量记录构成一个三阶协方差矩阵,此协方差矩阵的特征值和特征向量定义了一个椭球体,该椭球是协方差时窗内记录到质点运动的最小平方近似,特征值和其对应的特征向量则表征了质点振动的主要特征[18]。所构建的协方差矩阵为:

(5)

(6)

极化度η取值范围为0～1,η=0表示质点的振动轨迹为圆,η=1表示质点线性振动。

基于多维经验模态分解的背景干扰信号识别与压制方法关注连续监测记录中的背景干扰信号成分,对极化度低的高频信号可以设置阈值η0进行筛除。由于受随机信号或者信号频段相近的影响,经验模态分解可能会出现模态混叠的现象,为了提取长持续时间的强干扰信号同时尽可能保留随机微地震信号的能量,对可能包含微地震信号的本征模态函数分量作进一步偏振滤波处理。

由于微地震信号与背景干扰信号的偏振方向不同,且连续记录中背景干扰信号持续时间长,可以选择合适的滑动时窗长度,计算时间窗口内质点振动轨迹的偏振主方向,并求出时窗内偏振主方向和时窗内平均偏振主方向之间的夹角β。根据信号偏振方向的差异,定义一个自适应偏振滤波器。

x′(t)=x(t)f(t)

(7)

f(t)=cosp[β(t)]

(8)

式中:x(t),x′(t)为目标信号滤波前、后数据;f(t)为滤波系数;p为矢量方向夹角的影响控制系数;p取值为1～2。

1.4 方法流程

由于不同干扰信号的能量和偏振特征存在差异,其对有效微地震信号的影响程度不同,为了尽可能地去除强能量的干扰信号,同时避免去除干扰信号的过程中对微地震信号造成能量损失,依据干扰信号的能量强弱以及偏振属性提出强背景干扰信号的压制方法,技术流程见图8,主要步骤为:

图8 强背景干扰信号压制方法流程

1) 对三分量监测信号进行噪声辅助多维经验模态分解实现信号的分解,获得各阶次多维本征模态函数分量;

2) 求取各多维本征模态函数分量信号的能量和偏振属性,并将信号按照能量大小进行排序;

3) 对于能量较弱的多维本征模态函数分量根据信号的极化度属性判断随机噪声,对于能量较强的多维本征模态函数分量利用自适应偏振滤波提取背景干扰信号对应的成分;

4) 去除随机噪声以及背景干扰信号,将剩余的信号成分重构,得到压制干扰信号后的监测记录。

2 应用实例

将基于多维经验模态分解的背景干扰信号识别与压制方法应用于实际的井中微地震监测资料以证明方法的效果。以图3中的微地震事件记录为例,图9 为三分量记录经过多维经验模态分解的结果,图10 为各多维本征模态函数分量能量占比,图11为按照能量从大到小排序的各多维本征模态函数分量信号的质点振动轨迹及信号的极化度。从图中可以看出,信号中绝大部分能量位于低频段,高频随机噪声(多维本征模态函数第1分量)能量占比很小,如果仅去除高频随机噪声无法达到提高信噪比的目标。由于模态混叠的影响,微地震信号的能量分布在几个连续的多维本征模态函数分量中(多维本征模态函数第2～6分量),在低频成分中多维本征模态函数第4～6分量包含较强能量的背景干扰信号同时包含部分微地震信号能量,此时将这几个分量不作额外处理的保留或者去除都将无法达到理想的处理效果。

图9 三分量记录经噪声辅助多维经验模态分解结果(为了更清楚地显示信号,不同子图间未采用统一的振幅增益)

图10 各多维本征模态函数分量能量占比

图11 按照能量排序的多维本征模态函数分量质点振动轨迹及极化度(为了更清楚地显示轨迹,图中未采用统一标尺)

为了保留有效微地震信号成分的能量,我们对多维本征模态函数第5和第6分量进行自适应偏振滤波。图12和图13分别为多维本征模态函数第5分量和第6分量经过自适应偏振滤波处理的结果,图中红、绿、蓝线分别代表X、Y、Z三分量,从滤波结果可以看出,微地震S波信号的能量从背景干扰信号中分离出来,偏振轨迹的对比也说明了滤波后的记录中S波信号成分减少。通过去除分离后的强能量背景干扰信号以及高频随机噪声,可以得到处理后的记录,如图14所示。

图12 多维本征模态函数第5分量经过自适应偏振滤波后的结果a 多维本征模态函数第5分量; b 自适应偏振滤波系数; c 滤波后结果; d 滤波前(黑色)、后(红色)偏振轨迹对比

图13 多维本征模态函数第6分量经过自适应偏振滤波后的结果a 多维本征模态函数第6分量; b 自适应偏振滤波系数; c 滤波后结果; d 滤波前(黑色)、后(红色)偏振轨迹对比

图14 压制强背景干扰信号后的去噪效果a 原始记录; b 强背景干扰信号压制处理结果; c 强背景干扰信号压制去除的背景干扰信号

图15显示了实际微地震事件记录的背景干扰噪声的压制效果。图15a为原始记录,图15b为基于多维经验模态分解背景干扰信号识别与去除方法处理结果,图15c为该方法去除的背景干扰信号,图15d为带通滤波(50～200Hz)的处理结果,图15e为带通滤波去除的噪声成分,图中:1～12道为x分量,13～24道为y分量,25～36道为z分量。从图15可以看出,各级检波器噪声成分存在差异,此方法能够自适应地实现对背景干扰信号压制且压制效果比常规带通滤波方法更好,有效保留了微地震S波信号能量。

图15 实际微地震事件记录的去噪效果对比a 原始记录; b 强背景干扰信号压制方法处理结果; c 强背景干扰信号压制方法去除的背景干扰信号; d 带通滤波(50～200Hz)处理结果; e 带通滤波去除的噪声成分

长短时窗能量比方法(short-term average/long-term average,STA/LTA)是常用的微地震事件识别以及初至拾取的方法。图16为使用带通滤波和强背景干扰信号压制方法预处理后的P波、S波震相的长短时窗能量比比值对比,计算能量比值所选择的长、短时窗大小分别为0.050s和0.002s。图16a为图15中第6级检波器记录的长短时窗能量比曲线对比,经强背景干扰信号压制方法处理后的P波、S波震相初至处的长短时窗能量比比值相比于常规带通滤波处理结果更为突出,图16b为对该压裂段中46个微地震事件处理结果统计分析,从图中可以看出,除了极少数可能由于初至拾取不准确的原因导致能量比降低外,经强背景干扰信号压制方法处理后的初至位置的能量比值绝大多数明显增强,说明该方法可以提高微地震记录的信噪比,有助于微地震信号的检测。

图16 带通滤波与强背景干扰信号压制方法预处理后的P波、S波震相的长短时窗能量比比值对比a 图15中第6级检波器记录的长短时窗能量比曲线对比; b 该压裂段46个微地震事件震相初至处的长短时窗能量比对比

3 结论

微地震监测记录中存在的长持续时间和强能量的背景干扰信号,常规预处理方法对这些干扰信号压制效果并不理想,且过度滤波对微地震信号造成损害。基于三分量监测信号的同步与相关分析,提出了一种基于多维经验模态分解的背景干扰信号识别与去除方法。与常规方法的不同之处在于此方法以背景干扰信号为处理目标,尽可能去除背景干扰信号的同时避免微地震信号的损失。考虑到背景干扰信号时变非平稳的特征,该方法通过多维经验模态分解得到不同阶次的本征模态函数分量,并根据偏振分析判断多通道信号成分实现干扰信号的识别、分离和压制。与常规的预处理方法相比,此方法在高效去除监测记录的强能量背景干扰的同时,较好地保留了微地震信号的有效能量成分,因而提高了微地震监测记录的信噪比,有利于后续事件识别与初至拾取。