2018, Vol. 61
随着油气资源需求的日益增长，常规油气资源不断消耗, 非常规油气资源的开采已迫在眉睫.微地震监测技术作为一种极为有效的非常规油气勘探手段，在油气资源开采中发挥着越来越重要的作用(Maxwell and Urbanic, 2001; Maxwell et al., 2010).地面微地震监测是微地震压裂监测领域中十分重要的一项技术，其特殊的震源产生和接收方式使得所采集到的数据，有效信号能量微弱，给数据处理带来了较大的困难(Shemeta and Anderson, 2010; Kushnir et al., 2014).如何对地面微震数据中所包含的随机噪声进行有效的压制，较好的保留有效信号的能量，提高信噪比，是地面微震数据要解决的难题之一.
过去的数年时间里，国内外优秀的地震勘探和信号处理领域的专家就如何压制微地震数据中的随机噪声进行了不懈的的探索.极化滤波作为一种空间滤波手段，通过对微地震数据中的有效信号和噪声之间的偏振差异进行研究，达到压制噪声干扰的目的(Du et al., 2000).但由于其单一的方向性，使其无法完成复杂情况下的信噪分离.F-k滤波是一种二维的滤波技术，其主要通过利用微地震数据中有效信号与噪声之间的视速度差异来进行噪声的压制(宋维琪等，2008).但在其窄带二维滤波的过程中容易发生信号畸变的现象，导致假同相轴的产生.2012年基于多尺度分析的小波变换被应用于微地震数据的噪声压制当中(徐宏斌等，2012).该方法是通过对微地震信号进行多尺度分解，对分解到各个尺度层的小波系数设置阈值，以达到噪声压制的效果.此方法在微地震数据噪声压制中起到一定的作用，但由于小波变换在二维及高维数据的处理中存在方向局限性，只能在有限的方向上表示有效信号的信息，进而对微地震数据处理的效果造成一定的影响(林洪波等, 2011; Mousavi and Langston, 2016).
近年来多尺度几何分析得到了长足的发展，一些多尺度变换如Ridgelet(脊波)变换(Candès，1998)，Curvelet(曲波)变换(Donoho and Duncan, 2000)，Contourlet(轮廓波)变换(Do and Vetterli, 2005)以及Shearlet(剪切波)变换(Guo and Labate, 2007a，2007b)等，基于它们的多尺度性多方向性以及各向异性等特点，使得有效信号的分解稀疏性增强，为微地震数据的噪声压制提供了有利条件.其中，具有合成膨胀仿射系统构造的Shearlet变换是一种接近最优的多维函数稀疏表示方法，它可以通过高维操作算子对一基本函数进行伸缩平移及旋转操作来完成多尺度多方向的稀疏分解(Easley et al., 2009).相较于其他多尺度几何分析方法，Shearlet变换具有更敏感的方向性，最优逼近的稀疏表示性，比起Curvelet和Contourlet变换等拥有更为简洁的数学结构，运算速率较高，更易实现离散化形式；Curvelet变换的频域划分方式为隔层细分的，相较于Shearlet变换的逐层细分来说更加繁琐，不利于连续数据与离散数据的处理，且在一定程度上削弱了它的稀疏表示性能；Shearlet变换与Contourlet变换的时域实现方式比较相似，但是Shearlet变换划分方向的滤波器比起Contourlet变换的方向滤波器在方向数目，及支撑集大小上都没有限制；同时，Shearlet变换也是多尺度领域中唯一可以统一处理连续与离散系统的变换方法.
在对微地震数据处理中，Shearlet变换可以很好的将数据进行稀疏化，根据有效信号与噪声在变换域的方向、频率与位置差异，Shearlet阈值收缩算法可以准确的提取有效信号，压制随机噪声，提高信噪比.传统的阈值方法是采用统一的阈值对Shearlet系数进行处理，该方法过度削减了Shearlet系数, 造成部分有效信息的丢失，之后虽有基于尺度阈值的改进算法，但其依旧不能很好地保留有效信号幅度.本文提出一种基于Context模型的Shearlet变换去噪算法, 本算法考虑不同尺度不同方向上有效信号与噪声的分布不同，设计自适应的Context模型，一定程度上增强信号，加强有效信号的识别，便于Shearlet系数的分组设计，考虑系数位置阈值优化问题，以便可以更好的保留有效信号幅值，提高微震数据信噪比.
1 Shearlet变换基本原理Shearlet变换是一种可以通过对基函数进行膨胀，平移及剪切等仿射操作，来生成具有不同特征函数的多尺度几何分析工具(Kittipoom et al., 2012; Guo and Labate, 2013).其可以自适应的跟踪二维信号奇异曲线的方向，并且可以通过尺度参数的变化较准确的描述信号奇异性特征.
1.1 连续Shearlet变换有一连续二维函数φ∈L2(R2)，对其分别进行膨胀，平移操作，得
|
(1) |
其中，DM是膨胀操作，定义为DMf(x)=|detM|－1/2× φ(M－1x)；Tt是平移操作，定义为Ttf(x)=f(x－t)；G为一个2×2的矩阵，
|
(2) |
其中，膨胀矩阵
|
(3) |
那么，函数f∈ L2(R2)的连续Shearlet变换定义如下:
|
(4) |
其中SHψ(·)表示为连续Shearlet变换，〈, 〉表示内积.当一连续Shearlet函数ψ∈ L2(R2)满足条件：
|
(5) |
其中，ξ1和ξ2表示二维频域坐标，
|
(6) |
其中SHψ－1(·)表示连续Shearlet逆变换.
1.2 离散Shearlet变换当选择{ (a, s, t)=(2－j, －l, k):j, l∈ Z, k∈Z2 }时，获得矩阵
|
(7) |
当
|
(8) |
根据公式(4)和(6)，函数F∈ L2(R2)的离散Shearlet变换定义如下：
|
(9) |
|
(10) |
其中〈, 〉表示卷积，F代表信号函数，ψj, l, k代表离散剪切波原函数，离散参数j, l, k分别对应着相应的参数a, s, t.从图 1中我们可以看出每一个ψj, l, k的频域形式都紧支撑着一组梯形对，这些梯形对将频域结构进行划分.
|
图 1 剪切波频域支撑 Fig. 1 Shearlet support tiling in frequency domain |
Shearlet变换可以把数据分解在不同的尺度上，噪声在不同尺度上的能量分布基本相同，其变换域系数会随着尺度的增加而相应减小被分配到不同的尺度上.因而，我们通过分辨原始信号与噪声在Shearlet变换域上的分布不同来区分有效信号与噪声.基于剪切波(Shearlet)多方向的特性，当基函数与地震信号方向一致时信号的剪切波系数值较大，而当两者之间方向差别较大时，其得到的剪切波系数就较小.因而，应用合理的阈值函数，滤除较小系数，再进行逆变换，便可实现去噪.然而，同一尺度层的Shearlet系数数值区间较广，统一设定的阈值对大部分系数来说并不是最佳的，基于方向的阈值在去噪结果的表现中也是有限的，因而我们考虑了引入Context模型，实现基于系数位置的阈值去噪方法.
Context模型原是应用在信号压缩编码领域中的一种方法，之后被Chang等(2000)引入在噪声的消除方面.Context模型是用来优化不同尺度不同方向的系数阈值，提高去噪结果的.对同一尺度的变换域系数来说，它的数值变化区间比较广，设置同一阈值并不能达到最佳的去噪效果.当将这一尺度系数细化到各个方向上来看时，由于不同的方向，有效信号与噪声的分布存在差异，造成它的各方向系数分布也不尽相同，因而设置方向阈值会更有利于信号与噪声的分离，达到更良好的噪声压制效果.但是对于每一方向的变换域系数而言，它的系数数值依旧会存在较大的变化波动，为了能够更加细化的利用各方向系数的分布特点，完善信号的细节效果，本文结合Context模型对各尺度各方向系数进行再划分，以一种接近位置阈值的方法来完成变换域系数的处理，以求能达到更高的细节处理效果.
图 2中展示了利用基于Context模型的Shearlet变换位置阈值算法相较于传统的Shearlet变换方向阈值算法的优越性，(a)是含噪地面微地震记录经过Shearlet分解后第三尺度第五方向的系数图，(b)是该系数图基于方向阈值收缩的结果图，(c)是基于位置阈值收缩的结果图.从两组阈值收缩结果的比较中看出，基于方向的阈值收缩方法在收缩系数的时候，去除噪声系数时, 削减了部分有效信号的系数.基于位置阈值的系数收缩方法，在去除噪声的同时也尽可能保留有效信号系数的完整性，因而基于位置阈值的去噪方法更有利于噪声消减及微弱有效信号的恢复.
|
图 2 第三尺度第五方向Shearlet分解系数图 (a)含噪地面微震记录系数图; (b)方向阈值去噪结果; (c)位置阈值去噪结果. Fig. 2 The coefficients map for the third scale of fifth direction after a shearlet decomposition (a) The coefficients map of a noisy surface microseismic record; (b) The result of direction threshold denoising; (c) The result of position threshold denoising. |
Context模型属于一种能量划分模型，针对各尺度各方向的变换域系数，利用系数内部存在的关联性将能量相近的系数划分成同一组，则存在每一能量组的系数方差近似等于该组有效信号能量与噪声能量之和，针对每一能量组系数，分别设置合理的阈值并进行系数收缩处理，则可实现信噪分离.Context模型打破了同一尺度空间系数的约束，对同一尺度不同方向的系数在不破坏其相似性的前提下重新排列分类，更好的利用了细节系数的强相关性，增强了信号降噪过程中部分被忽视的局部特征，使得信号的处理效果更佳清晰准确.
结合Shearlet变换多尺度多方向等性能在地震信号处理上的优势以及基于Context模型的位置阈值在有效地震信号细节保留上的特点，本文建立了基于自适应Context模型的Shearlet变换随机噪声压制算法.
将含噪的地面微地震数据进行Shearlet稀疏分解，记分解后的Shearlet系数为I(x, y)，x, y为坐标.定义Context模型为3×3方块，如图 3所示，Shearlet系数I(x, y)相对应的Context值C(x, y)计算为
|
(11) |
|
图 3 Context模型 Fig. 3 The Context Model |
Context系数考虑了每一位置系数的临近位置系数对该位置的影响，在不破坏原有系数特征的前提下，综合减少了噪声对于有效信号系数的影响，增强有效信号的强度，使得能量组划分能更加准确.由于不同尺度层不同方向之间所包含的有效信号的能量不同，设置相同的能量划分模型会影响计算的准确度和效率，因此根据不同尺度不同方向的系数能量及方差不同，建立自适应的Context模型，减少计算量，以提高计算精度.
在有效信息较多，噪声较少的方向，其系数能量相对较大，方差相对较小，对于这样的系数方向我们尽可能的更加细致的划分系数能量组，使得之后阈值收缩的结果精确度更高.在有效信息包含较少，噪声较多的系数方向，我们可以减少划分组数来提高计算效率.所以，我们给出自适应能量组划分为
|
(12) |
其中Ek (Ik(x, y))代表各尺度各方向的Shearlet系数能量，Vk (Ik(x, y))代表各尺度各方向的Shearlet系数方差，k代表各尺度层的方向数目，m为一常数，|_ _|表示为向下取整数，pj, k即为每一尺度每一方向的划分组数.
根据得到的能量组数pj, k，统计每一尺度每一方向的Context系数Cj, k(x, y)的能量范围，根据对应得到的能量组数pj, k，将这一能量范围均分为p组，并判断每一Context系数所属的能量组，将Context系数分组.根据划分得到的Context系数能量组，得到对应的Shearlet系数能量组，计算每一组的贝叶斯阈值并进行系数收缩.
针对各尺度方向的各能量组，设置不同的阈值来进行系数筛选，本文采用贝叶斯阈值来估计各能量组阈值，阈值为
|
(13) |
根据Shearlet系数则有
|
(14) |
其中噪声n的标准差σn估计为
|
(15) |
那么，则有有效信号s方差σs2估计为
|
(16) |
采用硬阈值函数：
|
(17) |
基于Context模型的Shearlet变换算法流程图如图 4.
|
图 4 CMST算法流程图 Fig. 4 The algorithm flow diagram of CMST |
综上所述，基于Context模型的Shearlet变换地面微震随机噪声压制的算法步骤如下.
(1) 将含噪的地面微地震记录进行Shearlet稀疏分解；
(2) 计算每一尺度及方向的Shearlet系数Context值，系数能量及方差，实现自适应能量分组的目的；
(3) 通过得到的Context模型自适应分组，得到相应的Shearlet系数能量分组，估计各个能量组的贝叶斯阈值，并对Shearlet系数进行阈值收缩处理；
(4) 对收缩后的Shearlet系数进行逆变换重构，得到噪声压制后的地面微地震信号记录.
3 实验结果 3.1 合成地震记录为了验证本文提出的去噪方法的可行性和适用性，首先将其应用于仿真模型上.利用Ricker子波合成一幅21道的地面微地震模拟记录，每道512个采样点，采样频率为1000 Hz，地上主频设为30 Hz，道间距为50 m，各层传播速度分别为1500 m·s-1, 2400 m·s-1, 和3600 m·s-1，而各层岩石密度则分别设置为1.8×103 kg·m-3，2.4×103 kg·m-3，3.0×103 kg·m-3，各层深度分别为100 m，200 m及100 m，如图 5a所示，对该合成的模拟地面微地震记录加入高斯白噪声使其信噪比达到-10 dB，如图 5b所示，图中可以明显看到第三轴被噪声污染比较严重.分别应用基于传统硬阈值的Shearlet变换方法和本文提出的基于自适应Context模型的Shearlet变换方法对该含噪模拟记录进行噪声压制处理，得到处理后结果图分别为图 5c及图 5d.
|
图 5 模拟地面微震处理结果 (a)纯净记录; (b)含噪记录; (c) TST处理结果; (d) CMST处理结果. Fig. 5 Results of a synthetic surface microseismic record (a) Pure record; (b) Noisy record; (c) Result of TST; (d) Result of CMST. |
从图 5c及5d的对比中我们可以看出，两种方法都对噪声起到了较好的压制效果，但本文方法对同相轴的恢复效果更好，噪声压制更彻底，对比方法的噪声残留较多并且对于第三层微弱同相轴的重建基本无效.
为了更清晰的比较两种方法的去噪效果差异，分别从图 5c及5d中选取同一道数据(第14道)进行频域与时域的对比.从频谱图 6a的对比中可看出，本文方法的能量保持较对比方法的效果好，同时噪声压制效果也较强，100 Hz以上幅值基本为零，而对比方法还有一定的幅度值.在时域单道波形对比图 6b中，可以看出在噪声较强条件下，本文方法的幅度保持效果较好，噪声压制也较彻底.为了方便观察，对该道信号每一同相轴的波形进行了放大，如图 6(c，d，e)所示，从图 6e中可以明显看出本算法在微弱同相轴恢复上的优势.
|
图 6 两类算法在第14道的频域图和波形图 (a)频域结果图；(b)时域单道波形图；(c)第一同相轴波形放大结果；(d)第二同相轴波形放大结果；(e)第三同相轴波形放大结果. Fig. 6 The amplitude spectrum and waveform of the two methods in the 14th channel (a) Results of amplitude spectrum; (b) Results of waveform; (c) Results of enlarged drawing in first event; (d) Results of enlarged drawing in second event; (e) Results of enlarged drawing in third event. |
表 1给出了两种方法的量化比较结果，从这些数据中可以看出，应用对比方法进行噪声压制后信噪比提高了15.0377 dB，本文方法提高了18.3741 dB，本文方法的均方根误差比起对比方法降低了0.0315，更加说明了本文算法的优越性.
|
表 1 信噪比和均方根误差定量比较 Table 1 The SNR and RMSE comparison |
为了验证本文算法的实用性，将其应用于一实际地面微地震数据记录中，对其进行噪声压制.图 7a为某地一幅道数为71道，每道1536个采样点的实际地面微地震数据记录图，其中部分区域存在大量的噪声.分别采用对比方法以及本文方法对其进行噪声压制，结果图分别为图 7b及7c.
|
图 7 实际地面微震数据比较 (a)实际地面微震数据；(b) TST处理结果；(c) CMST处理结果. Fig. 7 The comparison of field surface microseismic data (a) Field surface microseismic noisy record; (b) Result of the TST; (c) Result of CMST. |
对比实验中可以看出，两种方法都去除了大量的随机噪声，但本文算法的噪声压制后记录要优于对比方法，噪声压制更为彻底，且同相轴连续性更好，对比方法在强噪声区域消噪效果并不理想.为了进行更清晰的比较，将图中(红)方框的局部区域进行放大，如图 8(a，b，c)所示，从局部放大图中可以看出，在噪声较强的情况下，本文的方法可以较为清晰连续的恢复能量相对较弱的同相轴，而对比方法则无法，因而可以说明本文所提算法在处理实际地面微地震数据中具有一定的实用性和优势性.
|
图 8 图 7(a—c)局部放大区域比较图 (a)局部实际微震记录；(b) TST处理局部放大结果；(c) CMST处理局部放大结果. Fig. 8 The comparison of the zoomed parts in the rectangles of the Fig. 7(a—c) (a) Zoomed noisy record; (b) Zoomed result of the TST; (c) Zoomed result of the CMST. |
本文为了解决在较低信噪比条件下的噪声压制问题，提出了一种基于Context模型的Shearlet变换算法，并应用于地面微震信号噪声压制中，它通过利用Context模型进行信号增强，根据Shearlet系数能量方差进行系数分组，并进一步考虑系数的位置阈值，分别对每组系数设置合理阈值，改进了传统阈值方法的不足.本文算法在低信噪比条件下，可以有效的压制噪声，保持信号幅值完整，较好的恢复微弱有效信号，使得同相轴更加清晰连贯，信噪比有所提高.仿真实验及实际数据记录处理结果证明了本文算法在地面微地震噪声压制及微弱有效信号重构中的优势.
Candès E J. 1998. Ridgelets: Theory and applications[Ph. D. thesis]. Stanford University.
|
Chang S G, Yu B, Vetterli M. 2000. Spatially adaptive wavelet thresholding with context modeling for image denoising. IEEE Transactions on Image Processing, 9(9): 1522-1531. DOI:10.1109/83.862630 |
Do M N, Vetterli M. 2005. The contourlet transform:An efficient directional multiresolution image representation. IEEE Transactions on Image Processing, 14(12): 2091-2106. DOI:10.1109/TIP.2005.859376 |
Donoho D L, Duncan M R. 2000. Digital Curvelet transform:Strategy, implementation, and experiments. Proceedings of SPIE, 4056: 12-30. DOI:10.1117/12.381679 |
Du Z J, Foulger G R, Mao W J. 2000. Noise reduction for broad-band, three-component seismograms using data-adaptive polarization filters. Geophysical Journal International, 141(3): 820-828. DOI:10.1046/j.1365-246x.2000.00156.x |
Easley G R, Labate D, Colonna F. 2009. Shearlet-based total variation diffusion for denoising. IEEE Transactions on Image Processing, 18(2): 260-268. DOI:10.1109/TIP.2008.2008070 |
Guo K H, Labate D. 2007a. Optimally sparse multidimensional representation using Shearlets. SIAM Journal on Mathematical Analysis, 39(1): 298-318. |
Guo K H, Labate D. 2007b. Sparse shearlet representation of Fourier integral operators. Electronic Research Announcements in Mathematical Sciences, 14(29): 7-19. |
Guo K H, Labate D. 2013. The construction of smooth parseval frames of Shearlets. Mathematical Modelling of Natural Phenomena, 8(1): 82-105. DOI:10.1051/mmnp/20138106 |
Kittipoom P, Kutyniok G, Lim W Q. 2012. Construction of compactly supported shearlet frames. Constructive Approximation, 35(1): 21-72. DOI:10.1007/s00365-011-9142-y |
Kushnir A, Varypaev A, Dricker I, et al. 2014. Passive surface microseismic monitoring as a statistical problem:Location of weak microseismic signals in the presence of strongly correlated noise. Geophysical Prospecting, 62(4): 819-833. DOI:10.1111/gpr.2014.62.issue-4 |
Lin H B, Li Y, Xu X C. 2011. Segmenting time-frequency peak filtering method to attenuation of seismic random noise. Chinese Journal of Geophysics (in Chinese), 54(5): 1358-1366. DOI:10.3969/j.issn.0001-5733.2011.05.025 |
Maxwell S C, Urbancic T I. 2001. The role of passive microseismic monitoring in the instrumented oil field. The Leading Edge, 20(6): 636-639. DOI:10.1190/1.1439012 |
Maxwell S C, Rutledge J, Jones R, et al. 2010. Petroleum reservoir characterization using downhole microseismic monitoring. Geophysics, 75(5): A129-A137. |
Mousavi S M, Langston C A. 2016. Hybrid seismic denoising using higher-order statistics and improved wavelet block thresholding. Bulletin of the Seismological Society of America, 106(4): 1380-1393. DOI:10.1785/0120150345 |
Shemeta J, Anderson P. 2010. It's a matter of size:Magnitude and moment estimates for microseismic data. The Leading Edge, 29(3): 296-302. DOI:10.1190/1.3353726 |
Song W Q, Chen Z D, Mao Z H. 2008. Hydro-Fracturing Break Microseismic Monitoring Technology (in Chinese). Dongying: China University of Petroleum Press.
|
Xu H B, Li S L, Chen J J. 2012. A study on method of signal denoising based on wavelet transform for micro-seismicity monitoring in large-scale rockmass structures. Acta Seismologica Sinica (in Chinese), 34(1): 85-96. |
林红波, 李月, 徐学纯. 2011. 压制地震勘探随机噪声的分段时频峰值滤波方法. 地球物理学报, 54(5): 1358-1366. DOI:10.3969/j.issn.0001-5733.2011.05.025 |
宋维琪, 陈泽东, 毛中华. 2008. 水力压裂裂缝微地震监测技术. 东营: 中国石油大学出版社.
|
徐宏斌, 李庶林, 陈际经. 2012. 基于小波变换的大尺度岩体结构微震监测信号去噪方法研究. 地震学报, 34(1): 85-96. DOI:10.3969/j.issn.0253-3782.2012.01.008 |