基于压缩感知的复杂地表菲涅尔束偏移

doi:10.11720/wtyht.2019.0077

基于压缩感知的复杂地表菲涅尔束偏移

杨飞龙¹^,², 俞岱³, 孙渊³

1. 西安石油大学地球科学与工程学院,陕西西安 710065

2. 陕西省油气成藏地质学重点实验室,陕西西安 710065

3. 长安大学地质工程与测绘学院,陕西西安 710054

Fresnel beam migration method based on compressed sensing under complex topographic conditions

YANG Fei-Long¹^,², YU Dai³, SUN Yuan³

1. College of the Geoscience and Engineering,Xi'an Shiyou University,Xi'an 710065,China

2. Shaanxi Key Laboratory of Petroleum Accumulation Geology,Xi'an 710065,China

3. College of Geology Engineering and Geomatics,Chang'an University,Xi'an 710054,China

责任编辑: 叶佩

收稿日期: 2019-02-13 修回日期: 2019-09-12 网络出版日期: 2019-12-20

基金资助:

陕西省教育厅科研计划项目“基于压缩感知的起伏地表Kirchhoff型动态聚焦束偏移方法”. 19JK0668

Received: 2019-02-13 Revised: 2019-09-12 Online: 2019-12-20

作者简介 About authors

杨飞龙(1988-),男,西安石油大学地球科学与工程学院讲师,主要从事地震数据处理及偏移成像方法研究工作。

摘要

高斯束偏移是一种束类地震数据成像方法,该方法可以处理多波至走时,具有良好的成像精度。但在复杂地表条件下应用存在两个问题:一是原有的局部平面波分解方法计算精度不足,影响低信噪比数据成像质量;二是偏移精度取决于射线束传播算子,高斯束传播算子随着传播距离增大发散过快,难以同时保证浅层和中深层构造成像精度,并且高斯束射线传播算子在模型浅部覆盖范围不足,对复杂地表模型成像精度较低。为了解决上述问题,本文将基于压缩感知理论的平面波分解技术应用到复杂地表菲涅尔束偏移中,不仅有效提高了局部平面波分解精度,同时菲涅尔束传播算子有效解决了高斯束算子随传播距离增加束宽发散过大的问题,典型的模型算例验证了本文方法的有效性和稳定性。

关键词： 高斯束 ; 压缩感知 ; 菲涅尔束偏移 ; 复杂地表 ; 局部平面波分解 ; 束算子

Abstract

Gaussian beam migration (GBM) is a ray-based seismic imaging method,which can handle multi-arrivals and has good imaging accuracy.However,two problems should be taken into consideration when this method is used in complex topographic conditions:First,the original local plane wave decomposition method has insufficient calculation accuracy,which affects the imaging quality of low signal-noise ratio data.Second,the ray beam propagator is the most important factor for migration accuracy.With the enhancement of the propagation distance,the Gaussian beam width increases rapidly and cannot ensure the imaging accuracy of the near surface and depth structures at the same time.In addition,the insufficient coverage of the Gaussian beam in the shallow part of the model might affect the imaging quality of this region.In order to solve the above problem,the authors apply wave decomposition technique based on compressed sensing theory to the complex surface in Fresnel beam migration.It not only effectively improves the precision of local plane wave decomposition,but also solves the problem that the width of Gaussian beam operator is increased quickly with the increase of propagation distance according to the Fresnel beam operator.Typical numerical examples prove the validity and stability of this method.

Keywords： Gaussian beam migration (GBM) ; compressed sensing(CS) ; Fresnel beam migration(FBM) ; complex topography ; local plane wave decomposition ; beam propagator

PDF (2799KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

杨飞龙, 俞岱, 孙渊. 基于压缩感知的复杂地表菲涅尔束偏移. 物探与化探[J], 2019, 43(6): 1285-1290 doi:10.11720/wtyht.2019.0077

YANG Fei-Long, YU Dai, SUN Yuan. Fresnel beam migration method based on compressed sensing under complex topographic conditions. Geophysical and Geochemical Exploration[J], 2019, 43(6): 1285-1290 doi:10.11720/wtyht.2019.0077

0 引言

随着油气勘探的深入,地震勘探的工作重心逐渐向沙漠、戈壁、黄土塬等复杂地区转移^[1]。在这些勘探地区,受近地表高程变化造成的静校正问题难以获得高质量的地震剖面^[2]。高斯束偏移方法克服了Kirchhoff偏移不能处理多波至和单程波动方程偏移不能对陡倾构造精确成像的缺点,并且保留了射线类偏移方法对复杂地表的适应性。近年来束偏移在复杂地表条件下的应用得到了广泛关注。Ceveny将平面波分解方法引入到高斯束当中,对束偏移中平面波的分解奠定了基础^[3]。孙建国等在地震波传播理论的基础上考虑了地震波传播时射线束宽度与菲涅尔带半径的关系^[4]。Gray通过局部高程静校正的方法将衰减窗内接收点的高程校正到各个窗中心点所在的高度,解决了起伏地形对偏移方法的影响问题,在复杂地表条件下实现了共炮域高斯波束偏移^[5]。Hu和Stoffa提出慢度驱动的高斯束叠前深度偏移方法,在这方法中考虑了局部窗内数据的相似性,并根据这一性质定义了成像公式中的加权函数从而加强有效信号对偏移结果的贡献^[6]。Yue等在Gray方法的基础上,考虑了高程的影响并对局部倾斜叠加公式进行修正,提高了复杂地表下高斯束偏移的精度^[7]。岳玉波等把真振幅高斯波束偏移扩展到了复杂地表条件^[8]。孙辉等通过分析子波特征改进射线束偏移精度,提出了基于子波的波束偏移方法^[9]。Wu等从反演角度对平面波进行分解,将τ-p变换方法引入到束偏移中提高偏移成像精度^[10]。杨继东等将菲涅尔束偏移方法扩展到了复杂地表条件下^[11]。高成等将高斯束偏移方法扩展至时间域,提出了二维时间域高斯束偏移,改善了偏移算法对速度模型的依赖^[12]。王华忠等在局部平面波分解过程中引入压缩感知概念,提出了特征波成像方法^[13]。孙辉等将压缩感知思想引入到波束偏移成像中,提高了局部平面波分解精度,改善了复杂地表下成像质量^[14-15]。

文中通过对压缩感知理论及束偏移方法的研究,将基于压缩感知的平面波分解技术应用到复杂地表菲涅尔束偏移中。通过基于压缩感知的平面波分解方法来提高τ-p域下数据体的计算精度;同时使用菲涅尔束传播算子有效控制射线束宽度解决其随传播距离发散过快等问题,保证了复杂构造浅层及中深层构造同时精细成像,并使用典型的模型算例及不同方法对比验证了本文方法的有效性。

1 方法原理

1.1 局部平面波分解

常规束偏移方法使用局部倾斜叠加(即线性拉东变换LRT)来对某一个窗中心地震数据的局部平面波进行分解。将压缩感知(CS)引入到局部平面波分解中,从反演角度考虑构造如下形式的最小化目标函数:

(1)

J = ‖ D (ω) - L (ω) P (ω) ‖^{2} + ‖ W_{P} (ω) P (ω) ‖^{2},

其中,W_P(ω)为空间加权对角矩阵。上述公式对P(ω)求导,并令其为0可以得到下面的解:

(2)

P (ω) = [L^{H} (ω) L (ω) + W_{P}^{- 1} {(ω)]}^{- 1} L^{H} (ω) D (ω) 。

使用L₀范数来表示平面波分解结果的稀疏性,上述平面波分解问题可以表示为:

(3)

\begin{array}{l} \min ‖ P (ω) ‖_{0}^{0}, \\ S.T. ‖ D (ω) - L (ω) P (ω) ‖_{2}^{2} < ε 。 \end{array}

其中L₀范数表示方程非0解的个数极小,ε表示模型空间的噪声大小。如此便成功地利用有效信号的稀疏性作为约束条件,在压缩感知理论框架下将常规的LRT转化为一个压缩感知问题。

图1a为多个平面波的合成数据,图1b为该数据的常规LRT结果,图1c为基于压缩感知的LRT结果。从图中可以看出基于压缩感知的LRT可以获得分辨率更高、质量更好的平面波数据。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 不同平面波分解方法对比

a—平面波数据;b—LRT结果;c—CS-LRT结果

Fig.1 Comparison of different plane wave decomposition methods

a—plane wave data;b—LRT results;c—CS-LRT results

1.2 菲涅尔束算子

菲涅尔束传播算子本质上仍然是高斯射线束传播算子,高斯射线束的理论基础是波动方程的高频近似解:

(4)

u (s, n, ω) = \sqrt[]{\frac{V (s)}{Q (s)}} \exp \{i ωτ (s) + \frac{i ω}{2} \frac{P (s)}{Q (s)} n^{2}\},

其中s、n为中心射线坐标系的坐标,ω为角频率,u为波场值,V为中心射线上的速度值, $τ (s) = \int_{s_{0}}^{s} \frac{d s}{V (s)}$ 为射线束沿中心射线的旅行时,P(s)、Q(s)分别为动力学射线参数,可以表示为:

(5)

\frac{P (s)}{Q (s)} = \frac{ε p_{1} (s) + p_{2} (s)}{ε q_{1} (s) + q_{2} (s)} 。

在高斯射线束传播算子中ε是一个常数,但在菲涅尔束传播算子中,ε是一个与菲涅尔带半径相关的变量,其表达式为:

(6)

ε (s) = - i \frac{ω_{ref} r_{F}^{2} / 2 + \sqrt[]{| (ω_{ref} r_{F}^{2} {/ 2)}^{2} - 4 q_{1}^{2} (s) q_{2}^{2} (s) |}}{2 q_{1}^{2} (s)}

其中ω_ref为参考频率,r_F对应的表达式如下:

(7)

r_{F} (s) = \sqrt[]{\frac{π}{ω_{ref}} \frac{q_{2} (s)}{p_{2} (s)}} 。

地震波传播到检波点处的有效半宽度Ls可以表示为:

(8)

\begin{array}{r} Ls & = {\{\frac{ω}{2} Im [\frac{p (s)}{q (s)}]\}}^{- \frac{1}{2}} \\ = {\{\frac{ω}{2} Im [\frac{ε p_{1} (s) + p_{2} (s)}{ε q_{1} (s) + q_{2} (s)}]\}}^{- \frac{1}{2}} 。 \end{array}

根据上式可知,地震波传播是射线束的宽度与ε有关。图2为各向同性介质中不同射线束形态对比图,图2a、图2b分别对应高斯射线束、菲涅尔束,其中射线束起始点位于(500 m,0 m)处,射线束初射方向为0°,射线束的频率为20 Hz,参考频率为30 Hz,传播速度为1 500 m/s,高斯波束和菲涅尔束的初始波束宽度均为40 m。可以发现:高斯波束宽度随着距离的增加迅速增大,菲涅尔束则可以很好地控制射线束宽度。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 各向同性介质中不同射线束振幅实部示意

a—高斯射线束;b—菲涅尔束

Fig.2 Comparison diagram of beam geometry in constant velocity model

a—Geometry of Gaussian beam;b—Geometry of Fresnel beam

1.3 走时振幅计算

如图3所示,基于Cerveny有效邻域波场近似理论,起伏地表条件下菲涅尔束偏移中地震波传播到检波点R处的波场大小可以表示为:

(9)

\begin{array}{l} U (R, ω) = A (L) \exp \{i ω [τ (L) + p_{x}^{L} (x - x^{L}) + p_{s}^{L} (z - z^{L}) + \frac{1}{2} X^{T} Re (W) X]\} \\ \exp [- \frac{1}{2} |\frac{ω}{ω_{ref}}| \frac{(x - x^{L})^{2}}{l_{0}^{2}} + |\frac{ω}{ω_{ref}}| \frac{(x - x^{L}) (z - z^{L}) \sin θ^{L}}{l_{0}^{2} \cos θ^{L}} - \frac{1}{2} |\frac{ω}{ω_{ref}}| \frac{(z - z^{L})^{2} si n^{2} θ^{L}}{l_{0}^{2} co s^{2} θ^{L}}] \end{array}

上式中,第一个指数项为接收点波场相对于束中心波场的相位校正项,第二个指数项为相应的振幅校正项。振幅校正项中的- $\frac{1}{2} |\frac{ω}{ω_{ref}}| \frac{(x - x^{L})^{2}}{l_{0}^{2}}$ 是由于接收点和束中心水平坐标不一致引起的,在菲涅尔束偏移中表现为局部倾斜叠加的指数衰减项;- $\frac{1}{2} |\frac{ω}{ω_{ref}}| \frac{(z - z^{L})^{2} si n^{2} θ^{L}}{l_{0}^{2} co s^{2} θ^{L}}$ 和 $|\frac{ω}{ω_{ref}}| \frac{(x - x^{L}) (z - z^{L}) \sin θ^{L}}{l_{0}^{2} \cos θ^{L}}$ 是由于地表高程差异引起的振幅衰减项。相位校正项中τ(L)为束中心的实值走时,在偏移中作为时间的初始值,一般取为零; $p_{x}^{L}$ (x-x^L)是由接收点和束中心水平坐标不一致引起的旅行时差,在菲涅尔束偏移中体现为局部倾斜叠加的线性相位校正项; $p_{s}^{L}$ (z-z^L)是由于地表高程差异引起的线性相位校正项,该项对束偏移保幅成像起到重要作用,在一定程度上改善束偏移成像质量; $\frac{1}{2}$ X^TRe(W)X为由起伏地表和近地表速度变化引起的二次时差校正项,当地表高程和近地表速度剧烈变化时此项的影响不可忽略。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 起伏地表条件下坐标变换关系

Fig.3 The Coordinate transformation relation in topographic surface conditions

2 模型试算

2.1 起伏地表层状模型试算

图4a为各向同性介质起伏地表层状模型,模型网格大小为401×401,网格间距为2.5 m×2.5 m,模型中从上到下层速度分别为700、1 300、1 800 m/s。图4b为对该模型进行波动方程正演得到的地震波场记录,炮点位于横向位置为500 m处的起伏地表上,第一个检波点的横向位置为 0 m,共有101个接收道,道间距为10 m,道与道之间横向位置均匀排列,每一个接收道时间采样间隔为1 ms,记录长度为3 s。图4c为复杂地表条件下高斯束偏移结果,图4d为基于压缩感知的起伏地表菲涅尔束偏移结果。从偏移结果中可以看出,基于压缩感知的菲涅尔束偏移面对复杂地表条件平面波分解精度更高,成像效果更好。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 起伏地表层状模型偏移结果

a—速度模型;b—地震记录;c—高斯束偏移结果;d—基于压缩感知的菲涅尔束偏移结果

Fig.4 Migration result of layer model under complex topographic conditions

a—velocity model;b—seismic record;c—migration of GBM;d—migration result of CS-FBM

2.2 起伏地表Marmousi模型试算

图5a为起伏地表Marmousi速度模型,其横向有737个网格节点,横向网格间距为12.5 m,纵向上有750个网格节点,纵向网格间距为4 m。地震记录由240炮组成,炮间距为90 m,每炮有101个接收道,道间距为25 m,放炮方式为单边放炮,偏移距范围从0~2 500 m不等,每道有800个采样点,采样间隔为4 ms。图5b与图5c中偏移结果分别对应高斯束束偏移方法和文中的新方法,两种方法的射线角度范围均为-50°~50°,射线数量为25根。可以发现两者在模型中深部都可以得到较好的偏移成像效果,但在模型浅部,两种方法的偏移结果有明显差异。首先在起伏地表处,原方法对应的结果在起伏地表处偏移能量断断续续,无规律性可言,而本文方法的偏移结果在起伏地表处的偏移能量规律性更强,其次本文方法的偏移结果信噪比更高,反映的地质构造更加清晰。对比图中画白圈的断层构造可以发现,高斯束偏移成像剖面中较难识别出断层,而从本文方法对应的偏移结果中,可以较为清晰地识别出断层构造。造成这一现象的主要原因在于高斯射线束宽度在模型浅部较窄,在上述射线参数前提下,难以保证模型浅层的射线覆盖率,而文中新方法拥有初始波束宽度,从而避免了这一问题。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 起伏地表Marmousi模型偏移结果对比

a—速度模型;b—高斯束偏移结果;c—基于压缩感知的菲涅尔束偏移结果

Fig.5 Migration result of Marmousi model under complex topographic conditions

a—velocity model;b—migration of GBM;c—migration result of CS-FBM

3 结论

笔者在高斯束偏移算法研究的基础上,针对该方法在局部平面波分解步骤及射线束传播算子中存在的问题,分别使用基于压缩感知的局部倾斜叠加方法和菲涅尔束传播算子改进高斯束偏移算法,理论模型试算结果验证了对于复杂地表条件下的陡倾构造偏移成像效果,优于常规高斯束偏移算法。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

赵文智, 胡素云, 董大忠 , 等.

“十五”期间中国油气勘探进展及未来重点勘探领域

[J]. 石油勘探与开发, 2007,34(5):513-520.