基于补偿阈值的曲波变换地面微地震弱信号检测方法

引用本文

何柯, 周丽萍, 于宝利, 邓勇, 王丽丽, 王茜. 基于补偿阈值的曲波变换地面微地震弱信号检测方法[J]. 物探与化探, 2016,40(1): 55-60.
HE Ke, ZHOU Li-Ping, YU Bao-Li, DENG Yong, WANG Li-Li, WANG Xi. The ground microseismic weak signal detection method based on compensation threshold of curvelet transform[J]. Geophysical and Geochemical Exploration, 2016,40(1): 55-60.
Doi:10.11720/wtyht.2016.1.10 复制到剪切板

Permissions

基于补偿阈值的曲波变换地面微地震弱信号检测方法

何柯, 周丽萍, 于宝利, 邓勇, 王丽丽, 王茜

中国石油集团东方地球物理公司研究院,河北涿州 072751

作者简介: 何柯(1987-),男,硕士研究生,毕业于中国石油大学(华东),从事微地震资料处理与地震资料解释工作。

收稿日期: 2014-12-25

基金: 中国石油天然气股份有限公司“新疆大庆”重大科技专项(2012E-34-03)资助

摘要

地面微地震弱信号在经过曲波变换后得到的曲波系数有小的幅值,这部分小幅值系数与噪声的曲波系数产生重叠,采用常规阈值函数,不能区分该部分有效弱信号系数,会被当作噪声系数去除。对于地面微地震资料弱信号的处理,应该在尽少损伤弱信号的基础上进行噪声的压制。因此笔者在常规阈值函数的基础上,提出一种补偿阈值函数算法,在弱信号和噪声之间加入一个新的阈值,这样就可对之前损伤的与噪声曲波系数重叠的弱信号进行补偿。通过对模型的试算和实际资料的处理,补偿阈值函数算法能较好地对微地震资料中的随机噪声和相干噪声进行压制,提高资料的信噪比,达到弱信号检测的目的。

关键词: 地面微地震; 弱信号检测; 曲波变换; 补偿阈值; 信噪比

中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2016)01-0055-06 doi: 10.11720/wtyht.2016.1.10

The ground microseismic weak signal detection method based on compensation threshold of curvelet transform

HE Ke, ZHOU Li-Ping, YU Bao-Li, DENG Yong, WANG Li-Li, WANG Xi

BGP INC.,China National Petroleum Corporation,Zhuozhou 072751,China

Abstract

The ground microseismic weak signal curvelet coefficient has a small amplitude after curvelet transform;this part of the weak signal will be covered by the noise of curvelet coefficient.The conventional threshold function fails to differentiate effective weak signal from ineffective ones,and hence the effective weak signal will be removed as noise.For the ground weak signal in microseismic data processing,the weak signal should be tackled with less damage in noise suppression.In this paper,the authors put forward a compensation threshold algorithm on the basis of the conventional threshold function.The authors add a new threshold between the weak signal and noise,so that it can compensate the damaged weak signal overlapped with the noise of curvelet coefficient.Model test and actual data processing show that compensation threshold function algorithm is better in suppressing random and coherent noise of the microseismic data than the conventional threshold function.It can also improve the signal-to-noise ratio of data and achieve the purpose of weak signal detection.

Keyword: ground microseismic; weak signal detection; curvelet transform; compensation threshold; signal-to-noise ratio

Show Figures

在微地震资料中, 微地震弱信号指的是相对于噪声来说振幅和能量都很小, 被噪声部分或者完全淹没, 在剖面上形成较弱的同相轴, 甚至完全看不见轴的有效信号^{[1, 2]}。对于弱信号检测一般方法为在时域和频域进行噪声和信号的分离, 将噪声完全去除而只保留信号, 这样就达到了目的。而在实际处理中, 只将噪声完全去除的同时而不损伤信号这点是很难的, 有些方法虽然能去除很多噪声, 但同时也一定程度地损伤了有效信号。

小波变换适用于识别一维或点信号的奇异性, 但对于二维信号处理后效果不是很好。近几年提出的Curvelet变换方法理论, 比小波变换更具有优势, 不仅包括尺度信息及位置参数, 还有方位参数, 曲波对曲线的逼近方式是以支撑区间长方形为基的, 而小波变换是以正方形为基的, 因此曲波比小波更能逼近地震同相轴的趋势, 曲波变换在地震数据处理中有很多的应用^{[3, 4, 5, 6, 7]}。

1 Curvelet变换原理

1.1 脊波变换

定义双变量函数f(x, y)的Ridgelet变换为

$\begin{matrix} Φ_{a, b, θ} (x, y) = a^{- 1 / 2} Φ [(xcosθ + ysinθ - b) / a], \end{matrix}$ (1)

式中:Φ 为小波函数, a为尺度因子, b 为偏移截距, θ 为变换角度。

对二维数组f(x₁, x₂)而言, Ridgelet变换重构后公式为^[8]

$\begin{matrix} f (x) = \overset{2 π}{\int_{0}} \overset{\infty}{\int_{- \infty}} \overset{\infty}{\int_{0}} R_{f} (a, b, θ) Φ_{a, b, θ} (x) \frac{da}{a^{3}} db \frac{dθ}{4 π} 。 \end{matrix}$ (2)

1.2 拉东变换^[9]

设函数f(x, y)∈ L²(D), 在平面D区域中是平方可积的, 对它作Radon变换得

$\begin{matrix} R_{f} (θ, t) = \underset{D}{\int \int} f (x, y) δ [t - (xcosθ + ysinθ)] dxdy, \end{matrix}$ (3)

式中:δ 是冲激函数, t-xcosθ -ysinθ =0表示极坐标下的任意一条直线方程。

通过对Radon变换的系数进行Wavelet变换, 并进行分析, 可以得到Ridgelet变换系数

$\begin{matrix} R_{f} (a, b, θ) = \int R_{f} (θ, t) a^{- 1 / 2} Φ [(t - b) / a] dt, \end{matrix}$ (4)

式中:θ 为方位参数。

1.3 曲波变换的步骤^{[10, 11, 12, 13]}

1)子带分解

先定义子带波器P₀, 当二维目标函数f(x₁, x₂)经过P₀滤波后, 可以分解一系列子带f→ (P₀f, Δ ₁f, Δ ₂f, …), 其中P₀是低通滤波, Δ ₁, Δ ₂, …是带通滤波。

2)平滑分割

$\begin{matrix} Δ_{n} f \to (w_{Q} Δ_{n} {f)}_{Q \in Q_{n}}, \end{matrix}$ (5)

式中:w_Q为不同子带Q=[k₁/2ⁿ, (k₁+1)/2ⁿ]× [k₂/2ⁿ, (k₂+1)/2ⁿ]上的窗函数集, 分割每一个平滑子带。

3)正规化

$\begin{matrix} g_{Q} = 2^{- n} (T_{Q})^{- 1} (w_{Q} Δ_{n} {f)}_{Q \in Q_{n}}, \end{matrix}$ (6)

其中:(T_Qf)(x₁, x₂)=f(2ⁿx₁-k₁, 2ⁿx₂-k₂)。

4)Ridgelet 分析

$\begin{matrix} α_{u} < g_{Q}, ρ_{λ} >, u = (Q, λ) \end{matrix}$ 先对每个数据进行局部脊波重构, 再整合分割的部分, 最后重新组合成新数据。

1.4 阈值函数

信号和噪声经过Curvelet变换后得到的Curvelet系数是不同的, 所以可以选择合适的阈值和阈值函数来消除噪声, 保留信号。Chang等提出了一种基于 Bayes准则的阈值, 其依据为小波图像的特点^[14]

其中: $\begin{matrix} {\tilde{σ}}_{noise}^{2} \end{matrix}$ 和 $\begin{matrix} {\tilde{σ}}_{signal} \end{matrix}$ 分别为噪声方差 $\begin{matrix} σ_{noise}^{2} \end{matrix}$ 和信号标准差σ _signal的估计值。

Gnanadurai和 Sadasivam提出了依赖于分解层次的阈值^[15]

$\begin{matrix} T = 2^{L - K} σ^{2} - (| A_{M} - G_{M} |), \end{matrix}$ (8)

其中:L是信号分解, K是去噪层次, |A_M-G_M|为子带系数的算术平均和几何平均之差的绝对值。笔者选用的阈值为

$\begin{matrix} T_{h} = {\tilde{σ}}_{noise}^{2} / {\tilde{σ}}_{signal} - | A_{M} - G_{M} | 。 \end{matrix}$ (9)

虽然常规阈值函数去噪有一定的效果, 但都有不足之处。在受到强噪声干扰时, 硬阈值函数去噪后的信号可能会有动荡和突变, 且不连续。当Curvelet系数较大时, 软阈值函数去噪后的信号虽然有连续性, 但是重构后的信号与真实信号的逼近程度会有所欠缺。

改进的软硬阈值折中方法为

$\begin{matrix} c' (i, j) = \{\begin{matrix} sgn [c (i, j)] [1 - {|\frac{λ}{c (i, j)}|}^{m}], | c (i, j) | \geq λ \\ 0, | c (i, j) | < λ \end{matrix} \end{matrix}$ (10)

其中:c(i, j)为变换后的Curvelet系数, c'(i, j)为处理后的Curvelet系数, λ 为阈值, sgn()是判断正负的函数。

2 改进的补偿阈值函数

由于地面微地震监测资料中噪声较强, 信号部分或者全部淹没在噪声之中, 看到的有效信号能量也较弱, 地面微地震弱信号一般具有弱振幅、弱连续性的特点, 因此有效信号经过曲波变换后得到的系数会有小的幅值, 与噪声部分相重合。在式(10)中, 当曲波系数小于阈值时, 就设为0, 那么弱信号的曲波系数就被作为噪声系数而滤除, 这样在去除噪声的同时也把微地震弱信号所去除。那么对式(10)进行改进, 引入一个新的阈值T₀

$\begin{matrix} c' (i, j) = \{\begin{matrix} (1 - α) c (i, j) + αsgn [c (i, j)] [| c (i, j) | - T], | c (i, j) | \geq T \\ λsgn [c (i, j)] [| c (i, j) | - (T - T_{0})], T_{0} \leq | c (i, j) | < T \\ 0, | c (i, j) | < T_{0} \end{matrix} \end{matrix}$ (11)

其中:-1< α < 1, T₀=λ T, 0< λ < 1。这样, T₀和T之间所对应的噪声和信号的Curvelet 系数也能够被分离, 能防止有效信号的剔除, 保护了微地震弱信号。当T₀=1时, 该函数就是软硬阈值折中函数, 只要取合适的补偿系数λ 和折中系数α , 就能在尽量不损伤有效信号的程度上尽可能多的去除噪声, 这种改进的阈值方法为阈值补偿法^{[16, 17, 18, 19, 20]}。改进的曲波变换方法步骤为:

1)对地面微地震资料进行去直流分离、带通、强能量压制等一系列预处理, 得到预处理后的资料。

2)选择合适的分解层数, 对预处理后的微地震资料进行曲波正变换, 得到相应的Curvelet 系数。

3)求出噪声和信号方差, 运算阈值函数得到T和T₀。

4)运用改进的阈值补偿法, 选择合适的补偿系数和折中系数, 对Curvelet 系数进行阈值处理, 得到新的曲波系数。

5)将处理后的Curvelet 系数进行反变换重构得到去噪后的地面微地震信号。

3 理论模型和实际资料的试算

为了验证方法的可行性, 设计了微地震模型。地面微地震模型共有126道, 每级检波器间隔距离为15 m, 数据采样间隔为2 ms, 1 000个采样点。模型中共包含3个能量不同的有效事件, 三个微震源相对能量比为1:2:4, 如图1所示。对原始模型加入随机噪声后得到了信噪比小于1的微地震记录(图2), 能量最弱的一个有效事件的同相轴已被淹没, 可以称为弱信号。设原始信号为x, 降噪信号为y, 都是M× N二维信号, 则信噪比(R_SN)公式为

$\begin{matrix} R_{SN} = 10 \lg \{[\overset{M}{\sum_{i = 1}} \overset{N}{\sum_{j}} x_{ij}] / [\overset{M}{\sum_{i = 1}} \overset{N}{\sum_{j}} (y_{ij} - x_{ij})]\} 。 \end{matrix}$ (12)

	Figure Option View Download New Window
	图1 原始模型

	Figure Option View Download New Window
	图2 加噪模型

常规阈值算法虽然可以将弱信号凸显出来, 但是处理后剖面中还残留着少量的噪声, 处理后的微地震资料的信噪比为5.221(图3), 经过补偿阈值算法后, 地面微地震资料的信噪比提高到7.125(图4), 噪声得到了较好地压制, 资料的信噪比大幅度得到了提高。

	Figure Option View Download New Window
	图3 常规阈值去噪结果

	Figure Option View Download New Window
	图4 补偿阈值去噪结果

由于地面干扰源的存在, 地面微地震资料中会含有一些相干性的噪声, 一般出现在局部道内, 呈局部线性干扰。为了更好地模拟地面微地震记录, 在原始模型(图1)中增加一组能量较强的视速度低于有效信号的干扰噪声(图5a)。对比图5b和图5c, 常规阈值曲波变换去噪后剖面中相干噪声能量较强, 随机噪声基本去除完毕, 补偿阈值曲波变换去噪后剖面中噪声基本消除。通过式(12)的计算, 常规阈值算法后得到的资料的信噪比为2.672, 补偿阈值算法为5.832, 对于实际地面微地震资料中的噪声, 经过改进后的补偿阈值曲波变换算法能更好地提高资料的信噪比, 而且对于微地震资料中同时出现的相干噪声和随机噪声也能更好地进行压制。

	Figure Option View Download New Window
	图5 加噪模型常规阈值与补偿阈值去噪结果对比a— 加噪模型; b— 常规阈值去噪结果; c— 补偿阈值去噪结果

将算法运用到实际资料的处理中。选取地面微地震事件中一条测线的微地震记录, 如图6a所示, 采样间隔为2 ms, 共126道, 500个采样点, 原始数据中红框显示部分能量较弱, 资料信噪比很低, 剖面中存在大量随机噪声, 前51道有效信号能量相对于后几十道有效信号能量较弱, 同相轴连续性也没后面的好, 可以看作弱信号。从去噪后处理结果来看(图6b~图6f), 常规阈值曲波变换方法有了很大的改善, 虽然处理后资料的剖面含有的随机噪声较少, 弱信号已检测出来, 强能量的同相轴更连续, 常规阈值曲波变换方法在去除噪声时也去除了一些有效信号。补偿阈值曲波变换去噪结果是最佳的, 弱信号的同相轴较为连续, 很大程度上提高了信噪比, 信号保真上也更进一步, 易于更好地进行初至拾取, 为下一步的反演提高较精确的初至时间。

	Figure Option View Download New Window
	图6 地面微地震事件1中某测线微地震记录去噪效果对比a— 原始资料; b— 图a结果局部放大; c— 常规阈值去噪结果; d— 图c结果局部放大; e— 补偿阈值去噪结果; f— 图e结果局部放大

再对另一微地震事件进行处理。该原始资料采样间隔为2ms, 126道, 1 000个采样点, 如图7a所示, 剖面中还有大量随机噪声, 有效信号的同相轴不连续。对比图7b和图7c的去噪结果可以看出, 补偿阈值曲波变换方法处理后的效果较常规阈值算法有了改善, 较大的提高了信噪比, 同相轴清晰且具有较好的连续性, 信号恢复特征明显。

	Figure Option View Download New Window
	图7 地面微地震事件2中某测线微地震记录去噪效果对比a— 原始资料; b— 常规阈值去噪结果; c— 补偿阈值去噪结果

4 结论

针对地面微地震资料进过曲波变换后, 常规阈值算法在去噪的同时会一定程度地损伤有效信号, 而地面微地震信号本身就比较微弱, 而笔者提出的改进的补偿阈值算法正是对这部分有效信号进行一定的补偿, 通过对模型和实际资料的试算, 提高了地面微地震资料的信噪比, 且尽可能地不损伤有效信号, 达到微地震弱信号检测的目的。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	宋维琪, 陈泽东, 毛中华. 水力压裂裂缝微地震监测技术[M]. 东营: 中国石油大学出版社, 2008: 1-5. [本文引用:1]
[2]	Forghani-Arani F, Willis M, Haines S, et al. Noise suppression in surface microseismic date by τ-p transform[J]. The Leading edge, 2012, 31(12): 1496-2501. [本文引用:1]
[3]	姜宇东, 杨勤勇, 何柯, 等. 基于曲波变换的地面微地震资料去噪方法研究[J]. 石油物探, 2012, 51(6): 620-624 [本文引用:1]
[4]	包乾宗, 高静怀, 陈文超. Curvelet域垂直地震剖面波场分离[J]. 西安交通大学学报, 2007, 41(6): 650-654. [本文引用:1]
[5]	张恒磊, 张云翠, 宋双, 等. 基于Curvelet域的叠前地震资料去噪方法[J]. 石油地球物理勘探, 2008, 43(5): 508-513. [本文引用:1]
[6]	张恒磊, 刘天佑, 张云翠. 基于高阶相关的Curvelet域和空间域的倾角扫描噪声压制方法[J]. 石油地球物理勘探, 2010, 45(2): 208-214. [本文引用:1]
[7]	张军华, 梁晓腾, 傅金荣, 等. Curvelet变换及其在地震资料去噪中的应用[C]//CPS/SEG Beijing, 2009. [本文引用:1]
[8]	Ramesh N, AnatoLly I B, Domiinique G G, et al. Coherent and rand om noise attenuation using the Curvelet transform[J]. The Leading Edge, 2008, 27(2): 240-248. [本文引用:1]
[9]	Cand ès E J, Donoho D L. New tight frames of curvelets and optimal representations of objects with smooth singularities[R]. USA: Department of Statistics, Stanford University, 2002. [本文引用:1]
[10]	Cand es E J, Demanet L, Donoho D L, et al. Fast discrete curvelet transforms[J]. Multiscale Modeling Simulation, 2005, 5(3): 861-899. [本文引用:1]
[11]	Cand es E J, Demanet L. Fast discrete curvelet transforms, applied and computational mathematics[R]. California Institute of Technology, 2005: 1-43. [本文引用:1]
[12]	Herrmann F J, EOS-UBC, Wang D L, et al. Seismic date processing with Curvelet: A multiscale and nonlinear approach[C]//Expand ed Abstracts of 77^th SEG Annual Internet Meeting, 2007: 2220-2224. [本文引用:1]
[13]	Felix J H, Gilles H. Non-parametric seismic data recovery with curvelet frames[J]. Geophysics, 2008, 173: 233-248. [本文引用:1]
[14]	Chang S G, Yu B, Marti N V. Adaptive wavelet thresholding for image denoising and compression[J]. IEEE Transactions on Image Processing, 2000, 9(9): 1532-1546. [本文引用:1]
[15]	Gnanadurai D, Sadasivam V. Image denoising using double density wavelet transform based adaptive thresholding technique[J]. International Journal of Wavelets Multiresolution and Information Processing, 2005, 3(1): 141-152. [本文引用:1]
[16]	包乾宗, 陈文超, 高静怀. 基于第二代Curvelet变换的地震资料随机噪声衰减[J]. 煤田地质与勘探, 2010, 38(1): 66-70. [本文引用:1]
[17]	吴爱弟, 赵秀玲. 基于曲波变换的地震信号去噪方法[J]. 中国石油大学学报: 自然科学版, 2010, 34(3): 30-33. [本文引用:1]
[18]	张军华. 地震资料去噪方法——原理、算法、编程及应用[M]. 北京: 中国石油大学出版社, 2011: 180-189. [本文引用:1]
[19]	何柯. 微地震弱信号检测方法研究[D]. 青岛: 中国石油大学, 2013: 58-77. [本文引用:1]
[20]	彭才, 常智, 朱仕军. 基于曲波变换的地震数据去噪方法[J]. 石油物探, 2008, 47(5): 461-464. [本文引用:1]

2008

0.0

... 在微地震资料中,微地震弱信号指的是相对于噪声来说振幅和能量都很小,被噪声部分或者完全淹没,在剖面上形成较弱的同相轴,甚至完全看不见轴的有效信号^[1,2] ...

2012

0.0

2012

0.0

... 近几年提出的Curvelet变换方法理论,比小波变换更具有优势,不仅包括尺度信息及位置参数,还有方位参数,曲波对曲线的逼近方式是以支撑区间长方形为基的,而小波变换是以正方形为基的,因此曲波比小波更能逼近地震同相轴的趋势,曲波变换在地震数据处理中有很多的应用^[3,4,5,6,7] ...

2007

0.0

2008

0.0

2010

0.0

2009

0.0

2008

0.0

... 对二维数组f(x₁,x₂)而言, Ridgelet变换重构后公式为^[8] ...

2002

0.0

... 2 拉东变换^[9]设函数f(x,y)#cod#x02208 ...

2005

0.0

... 3 曲波变换的步骤^{[10,11,12,13]}1)子带分解 ...

2005

0.0

... 3 曲波变换的步骤^{[10,11,12,13]}1)子带分解 ...

2007

0.0

... 3 曲波变换的步骤^{[10,11,12,13]}1)子带分解 ...

2008

0.0

... 3 曲波变换的步骤^{[10,11,12,13]}1)子带分解 ...

2000

0.0

... Chang等提出了一种基于 Bayes准则的阈值,其依据为小波图像的特点^[14] ...

2005

0.0

... Gnanadurai和 Sadasivam提出了依赖于分解层次的阈值^[15] ...

2010

0.0