基于构造导向滤波的高斯束层析速度建模方法及其应用

doi:10.11720/wtyht.2020.1300

基于构造导向滤波的高斯束层析速度建模方法及其应用

郑浩¹, 蔡杰雄¹, 王静波²

1. 中国石化石油物探技术研究院,江苏南京 211103

2. 中国石化勘探分公司,四川成都 610041

Gaussian beam tomography with structure-filtering and its applications

ZHENG Hao¹, CAI Jie-Xiong¹, WANG Jing-Bo²

1. Sinopec Geophysical Research Institute,Nanjing 211103,China

2. Sinopec Exploration Company,Chengdu 610041,China

责任编辑: 叶佩

收稿日期: 2019-06-3 修回日期: 2020-01-6 网络出版日期: 2020-04-20

基金资助:

国家重点研发计划项目. 2017ZX05036005-006

Received: 2019-06-3 Revised: 2020-01-6 Online: 2020-04-20

作者简介 About authors

郑浩(1991-),男,中国石化石油物探技术研究院,主要从事地球物理反演与Q值建模、智能化地震资料处理等研究工作。Email:zhenghao.swty@sinopec.com 。

摘要

高斯束层析速度建模方法利用高斯束传播算子计算层析核函数,能够提高速度模型反演的精度与稳定性,但尚未引入地质先验信息有时会导致反演结果缺乏地质意义。基于此,通过引入预条件模型正则化算子,推导得到构造导向滤波约束下的高斯束层析速度建模方法,该方法的关键是构造导向滤波算子的构建。本文通过地震剖面计算构造的结构张量,利用各向异性光滑算子加入地质构造特征作为反演的正则化约束条件。该方法完全基于数据驱动,实现了对层析反演的“软约束”,既缓解层析反演求解的多解性,同时又提高反演分辨率,得到更加符合地质认识的高精度速度模型。数值模拟实验及实际数据应用验证了构造导向滤波下的高斯束层析速度建模的有效性与实用性。

关键词： 高斯束层析 ; 速度建模 ; 构造导向滤波 ; 结构张量 ; 正则化约束

Abstract

Gaussian beam tomography uses the Gaussian beam operator to compute the travel-time tomographic kernel function instead of the conventional ray tomographic kernels,which can improve the accuracy and stability of velocity building.However,the absence of prior geological information sometimes leads to the lack of geological understanding of the results.On such a basis,by introducing the preconditioned model regularization operator,a method of Gaussian beam tomography velocity building under the constraint of structure-guided filter is derived,whose key is to construct the structure-guided filter operator.In this paper,structure tensors were calculated by seismic profiles,and geological structure features were introduced by anisotropic diffusion smoothing operators as regularization constraints for inversion.This method is entirely data-driven and achieves the "soft constraint" of tomographic inversion,which not only reduces the issue of multiple solutions of inversion but also improves the inversion resolution,and hence can obtain a reasonable velocity model with geological recognition.The validity and practicability of the Gaussian beam tomography velocity building with structure-guided filtering were verified by its application to synthetics and practical data.

Keywords： Gaussian beam tomography ; velocity building ; structure-guided filter ; structure tensors ; regularization

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

本文引用格式

郑浩, 蔡杰雄, 王静波. 基于构造导向滤波的高斯束层析速度建模方法及其应用. 物探与化探[J], 2020, 44(2): 372-380 doi:10.11720/wtyht.2020.1300

ZHENG Hao, CAI Jie-Xiong, WANG Jing-Bo. Gaussian beam tomography with structure-filtering and its applications. Geophysical and Geochemical Exploration[J], 2020, 44(2): 372-380 doi:10.11720/wtyht.2020.1300

0 引言

近年来随着“两宽一高”地震采集和层析反演方法取得突破进展,地震成像分辨率得到有效拓宽^[1,2], 高分辨率层析反演速度建模方法取得了长足的进步。传统的深度域射线层析反演速度建模方法的分辨率存在中波数段信息缺失问题^[3],其可以较好地恢复地下速度的背景低波数分量,但难以准确反演中波数段的高精度速度模型^[4,5]。学者们提出波动方程层析反演方法^[6],将速度反演分辨率向中高波数端扩展。该类方法基于波动理论,避免了高频近似假设引入的误差,反演结果精度高,但该类方法存在计算效率问题,且反演精度非常依赖于数据品质,当数据较差时存在反演稳定性问题^[7,8],在实际应用中推广较少。为了降低射线类层析矩阵的稀疏性,同时克服波动方程层析的计算量与稳定性问题,诸多学者提出介于常规射线理论和波动方程理论之间的方法。早在1995年,Vasco^[9]便提出一种介于射线和波动方程之间的胖射线层析方法,该方法通过将射线加宽,提高射线路径的投影范围,从而提高层析精度,但该方法缺乏严格的理论推导,仅是对波动理论感性认识,应用较为局限。菲涅尔体层析^[10]及高斯束层析^[11,12]的出现实现了对波动方程线性化近似的“束”层析。菲涅尔体层析利用射线周围的波动性质来拟合走时和波场振幅,相比于射线层析的高频近似假设,该方法引入了菲涅尔体的有限频效应,更接近真实的波场传播,这降低常规射线层析方法中矩阵的稀疏性,提高计算稳定性^[13]。高斯束层析^[14]实际上是将临近区域空间内的多条射线合并为一条高斯束射线从而建立层析方程,对于观测系统中的一个炮检对,通过提高射线覆盖范围来提高反演精度及分辨率,降低层析方程的病态性,且该方法相对于波动方程层析,计算速度快、计算结果稳定。但早期的高斯束层析理论推导不够严格,蔡杰雄^[15]从高斯束偏移角度道集出发,在波动方程一阶Born近似和Rytov近似下,推导出了基于高斯束算子的成像域走时层析方程及其敏感度核函数,完善了高斯束层析反演的理论基础,推动了高斯束层析反演的实用化。

层析反演速度建模的另一个重要研究内容是各种约束条件的利用,即在反演过程中加入包括层位、构造产状、井速度等多种已知信息约束以实现正则化,提高反演稳定性及分辨率。目前,层析反演中的构造约束方式主要有两类,一类是在速度层析过程中引入地质层位信息,通过层约束提高反演精度,但这类方法需要解释人员进行地质构造的精细化解释,对解释结果有较高的精度要求,若解释结果出现错误,会使得层析结果误差较大。该类方法大都基于主观认知强行加入反演过程中,比如层位约束层析反演就是在本轮的反演中不修改层位而仅更新层位内速度值,待完成本轮反演后重新偏移成像,在新剖面上重新进行层位拾取,之后再加入到新一轮的层析速度更新中,通常我们称这种约束为“强”约束^[16,17]。强约束的反演具有双刃性,它一方面可能快速促进收敛,但另一方面可能引入解释人员的主观错误信息,导致反演结果出现误差;另一类方法通过引入含有构造信息的正则化算子,用于导引速度更新方向,避免了解释人员的主观错误信息,但构造约束精度取决于正则化算子的精度。我们称这种约束为“软”约束。目前有多种预条件正则化算子,包括Sobel算子、Reborts算子、基于斜率的平滑滤波算子、基于扩散方程解析解的构造导向滤波算子等。这些方法通过提取地震偏移图像的构造特征,从而构建预条件算子,实现构造约束层析反演。

本文从高斯束层析核函数出发,通过引入含有构造信息的预条件模型正则化算子,推导得到构造导向滤波约束下的高斯束层析速度建模方法。在预条件模型正则化算子的设计上,本文引入结构张量算法,用于提取偏移图像中的平坦、边缘及角点区域,通过扩散方程解析解实现构造导向平滑滤波,建立层析反演正则化算子,实现反演过程的自动“软”约束,提高高斯束速度层析反演的精度与稳定性。

1 成像域高斯束层析方法

频率域高斯束成像条件可以表示为:

(1)

I_{GBM} (x, θ, φ, w) = S (x, p_{S}, w; x_{S}) R^{*} (x, p_{R}, w; x_{S})

式中:I_GBM(x,θ,φ,w)表示共成像点道集,S(x,p_S,w;x_S)表示来自激发点的高斯束下行波场,R(x,p_R,w;x_S)表示检波点接收的高斯束上行波场,x=(x,y,z)表示共成像点位置的坐标,p_S和p_R分别代表激发点和接收点位置的慢度矢量,θ表示反射张角,φ表示成像点的方位角,w=2πf表示角频率,*表示共轭。

从角度域高斯束偏移成像条件(1)出发,在波动方程的一阶Born近似下,可得到剩余时差与速度修正量的表达式,再引入Rytov近似可进一步简化得到剩余时差与速度修正量的线性关系式^[18]:

(2)

\begin{array}{l} Δ t_{GBM} (x, θ, φ) = \int_{V_{S}} K_{S}^{T} (y; x, p, x_{S}) Δ v (y) d y + \\ \int_{V_{R}} K_{R}^{T} (y; x, p, x_{S}) Δ v (y) d y \end{array}

其中,Δt(x,θ,φ,ω)为高斯束偏移的剩余时差,K^T表示高斯束层析核函数,其可以表示为两部分:

(3a)

K_{S}^{T} (y; x, p_{S}, x_{S}) = \int W (w) \cdot Im [\frac{2 w}{v_{0}^{3} (y)} \frac{G_{S} (x; p_{S}, w, y) G_{0} (y; p_{S}, w, x_{S})}{G_{0} (x; p_{S}, w, x_{S})}] d w

(3b)

K_{R}^{T} (y; x, p_{R}, x_{S}) = \int W (w) \cdot Im [\frac{2 w}{v_{0}^{3} (y)} \frac{G_{R} (x; p_{R}, w, y) G_{0}^{*} (y; p_{R}, w, x_{S})}{G_{0}^{*} (x; p_{R}, w, x_{S})}] d w

式中: $K_{S}^{T}$ 为 $K_{R}^{T}$ 分别表示激发点到成像点和成像点到接收点的层析核函数;W(w)表示归一化的加权函数,可以有多种取法,刘玉柱^[19]采用归一化的高斯函数,王华忠^[8]采用模拟数据的归一化能谱作为权函数,本文采用高斯函数作为加权函数;v₀表示初始速度;G_S(x;p_S,ω,y)表示震源波场中某点y到成像点x的格林函数;G_R(x;p_R,ω,y)表示检波点波场中某点y到成像点x的格林函数;*和Im分别表示共轭和虚部。由上面可知,在波动方程线性化近似下推导出的高斯束层析核函数包含炮点到成像点以及成像点到检波点两部分,见图1。显然,相比于常规射线层析,高斯束层析核函数模拟的波传播路径更加接近波的实际传播方式,这缓解了射线方法中矩阵的稀疏性问题,反演精度更高。

基于式(2)表示的层析核函数,就可以建立类似于射线层析的方程:

(4)

A Δ m = Δ t 。

图1

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

图1 深度域高斯束层析反演核函数

Fig.1 Gaussian beam tomography kernel function in imaging domain

A是利用式(2)层析核函数计算得到的敏感度核函数矩阵,Δm是待反演的速度更新量,Δt是共成像点道集(common image gathers,CIG)的走时残差,在反射点局部,可以假定炮点和检波点的慢度p_S和p_R是常数,那么该点的慢度可以表示为:

(5)

Δ p_{z} = Average (p_{S} + p_{R}) 。

通过在CIG道集上拾取剩余深度差Δz(RMO),则可利用式(5)得到走时残差Δt:

(6)

Δ t = Δ p_{z} Δ z 。

通常,式(4)中的矩阵A是极度稀疏的,直接求解稳定性较差,且多解性较强,往往需要引入正则化项,加入一些先验约束,使解在约束条件下收敛。这里采用预条件模型正则化的方式加入先验信息,待反演的模型扰动量可表示为:

(7)

Δ m = Du 。

其中算子D表示预条件算子,即引入的先验信息,这样可以将Δm的求解问题转化为算子u的求解问题,将式(7)代入到式(4)中,可得:

(8)

ADu = Δ t 。

式(8)即为构建的预条件层析方程,其阻尼最小二乘方程可表示为:

(9)

D^{T} A^{T} ADu + εu = D^{T} A^{T} Δ t,

求解上式得到u,代入到式(7)中即可得到速度更新量Δm。

通过上述推导,实现了波动方程线性化的深度域走时层析反演方法,再与叠前深度偏移相结合进行迭代:速度误差带来成像不聚焦误差,表现为成像道集走时差;再利用成像误差反过来修正速度,实现两者的迭代收敛,获得更精确的速度场和成像结果。

2 基于结构张量构造导向滤波正则化技术

从式(9)中可以看出,预条件模型正则化的高斯束层析反演关键是算子D的构建,传统方法常选择空间平滑算子,主要是通过减小层析系数矩阵的条件数,缓解反问题病态性,在矩阵求解的稳定性方面有一定作用,但该方法容易模糊构造边界,其反演结果对于边界的刻画较差。本文采用各向异性扩散方程解析解^[20]构建预条件算子D,其隐式表达式为如下的偏微分方程:

(10)

g (x) - α \nabla \cdot T (x) \nabla g (x) = f (x),

其中:f(x)是待滤波的地震数据,g(x)是滤波后的地震数据,α是表示平滑力度强弱的常数;T(x)是从地震剖面中计算得到的结构张量。结构张量实际上相当于沿构造梯度的一个平滑解,它可以用于构造的倾角扫描,且平滑过程沿构造进行,达到保边界的效果。对于一个三维数据体而言,其某点的结构张量T(x)可以通过一个3×3的对称正定矩阵表示:

(11)

T (x) = [\begin{array}{l} d x^{2} & d x d y & d x d z \\ d x d y & d y^{2} & d y d z \\ d x d z & d y d z & d z^{2} \end{array}],

其中:dx、dy、dz分别为数据沿x,y,z三个方向的梯度。通过矩阵的特征值分解,上式可以表示为:

(12)

T (x) = λ_{u} u u^{T} + λ_{v} v v^{T} + λ_{w} w w^{T},

其中:λ_u、λ_v和λ_w为分布在0~1的正实数,表示矩阵T(x)的特征值,u、v和w是对应的特征向量。通常定义λ_u≥λ_v≥λ_w≥0,这种情况下,特征向量u表示梯度最大的方向,即表示垂直于构造倾角的方向。特征向量v和w分别表示inline和crossline平行于构造的方向,λ越大意味沿对应特征方向的平滑尺度越大。

这里以二维数据为例,对特征向量进行说明。二维情况下,式(12)只有前两项,图2a、b中黄色的椭圆分别是地震切片与剖面上结构张量计算结果, 实际上在无构造(各向同性)区域,该点的结构张量形态是圆形,如图3a所示,此时u和v的长度相等;当存在构造形态(各向异性)时,就会出现椭圆形态的结构张量,如图3b所示,其中v表示平行于构造方向的特征向量,u表示垂直于构造方向的特征向量,这样我们就可以计算出所有点的结构张量。

图2

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

图2 地震切片(a)及剖面(b)计算得到的结构张量

Fig.2 Structural tensor calculated on seismic slice (a) and section (b)

图3

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

图3 各向同性(a)和各向异性(b)的结构张量算子

Fig.3 Isotropic(a) and anisotropic(b) structural tensor operators

得到每点的结构张量后,通过合理调整λ_u、λ_v、λ_w,从而实现沿层位和断层的构造导向滤波处理,平滑效果如图4所示。

图4

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

图4 构造导向滤波前(a)后(b)对比,可同时加强层位与断层特征

Fig.4 The comparison of (a) before and (b) after structural-guided filtering,which can enhance the structure and fault characteristics

引入有限差分近似,预条件算子D可写作:

(13)

D = (I + \nabla^{T} {T \nabla)}^{- 1},

其中T即为式(12)所求的结构张量,单位矩阵I的作用是在λ_u、λ_v、λ_w均为零时保持图像不改变,将式(13)所得结果代入式(9)中即可实现基于构造导向滤波的高斯束层析反演。此外,相比于其他的构造约束算法,该方法无需指定平滑角度,完全自动提取,计算灵活;且方法的抗噪性较好,即便是在信噪比极低的情况下,构造导向滤波算法也可以得到相对精确的结果。

3 模型及实际资料测试

3.1 二维层状模型数据测试

设计二维层状地质模型验证本文方法的精度,模型中填充速度属性值。所设计模型为速度沿纵向递增的层状模型, 其中模型中部有一高速隆起构造(图5a),纵向采样点为601,采样间隔为10 m,横向共有1201道,道间距为10 m。利用声波波动方程正演得到单炮记录(图6),观测系统采用中心放炮,两边接收,每炮801道。初始模型采用常梯度模型,见图5,其中第一层速度给成准确值。图5c、5d展示了常规高斯束层析与构造约束预条件高斯束层析的反演结果,显然,常规的高斯束层析对于浅层的水平层状构造反演结果尚可,但对于中部隆起构造反演结果较差,分辨率明显不足,且中深部反演结果误差较大;相比较而言,基于构造导向滤波的高斯束层析算法所得结果精度较高,与真实模型基本吻合,更有利于精确成像。

图5

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

图5 理论地质模型(速度值)及反演结果

a—理论模型;b—初始模型;c—常规高斯束层析反演结果;d—构造导向滤波高斯束层析反演结果

Fig.5 Velocity model and the inverted results

a—velocity model;b—initial model;c—updated velocity model by conventional Gaussian beam tomography;d—updated velocity model by Gaussian beam tomography based on structural-guided filtering

图6

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

图6 模型正演的单炮记录

Fig.6 The modelling common shot gathers

为了更加直观地看出造构造导向滤波高斯束层析的优势,在层析结果竖直抽取某一道进行对比,如图7a所示。图7b表示了理论速度值、平滑速度值和反演的速度值的对比曲线。,显然,基于构造导向滤波的高斯束层析反演方法不仅可以反演得到中波数分量的结果,提供更加丰富的速度信息,且能够达到保边界的效果,构造边界刻画更为清晰,这说明基于构造导向滤波的高斯束层析对于构造的刻画明显优于常规高斯束层析方法,且反演精度更高,减少了由于人为干预带来的误差,结果更加合理。

图7

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

图7 反演结果单道对比分析

a—抽取模型中心道;b—真实值、常规平滑正则化约束层析反演结果及构造约束正则化层析反演结果单道对比

Fig.7 Single trace comparison

a—the location of the extracted trace in the model;b—comparison of theoretical result,conventional smoothing regularization tomography inversion result and structural-guided regularization tomography result

深度域高斯束层析反演技术可以缓解层析矩阵的病态问题,反演得到速度模型的中波束分量,提高地震勘探的成像精度。基于此,通过引入基于结构张量的构造导向滤波算法用于约束速度更新量的计算,能反演出高精度的地下速度模型,增强断层突变位置的反演效果,弥补地震成像的中波数段缺失,提高地震成像精度。

3.2 实际地震资料应用

实际资料选自中国西南部复杂山地页岩气某工区。该地区地表起伏大,悬崖峭壁密布,山高谷深,数据信噪比非常低,且地下结构复杂,速度横向变化剧烈,速度建模难度大,成像困难,具有典型的复杂地表和复杂构造的“双复杂”特性,其对地震资料成像精度要求远高于常规勘探。

首先对该区域的资料进行预处理、叠前时间偏移、深度域初始建模及偏移成像等基础工作,得到高品质的CIG道集,在此基础上,通过基于构造导向滤波的高斯束层析反演得到深度域速度模型,结合克希霍夫叠前深度偏移得到最终的成像剖面。图8是某实际数据深度域常规高斯束层析建模和构造导向滤波高斯束层析速度建模结果对比,从两张对比图上可以看出,构造导向滤波高斯束层析速度建模结果表现了速度模型更多细节,清晰地刻画了速度纵横向空间变化,提高了速度建模的精度。

图8

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

图8 某实际数据深度域常规高斯束层析建模(a)及构造导向滤波高斯束层析建模(b)结果对比

Fig.8 Comparison of conventional Gaussian beam tomography (a) and structure-guided filter Gaussian beam tomography (b) in depth domain

速度模型的更新和剩余速度分析采用速度更新量和剩余速度谱联合质量控制方法,图9是质控图,图10是更新前后的CIG道集。从更新前的道集上看,存在同相轴未校平现象,更新后这些同相轴变平,图上前头所指。由于存在各向异性,还有一些同相轴未校平。另外,该区域存在多次波的影响(如图中红框区域所示),因此仍有部分道集下拉现象。

图9

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

图9 速度更新量(a)与剩余谱(b)联合质控

Fig.9 Joint quality control of velocity update (a) and residual spectrum (b)

图10

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

图10 速度更新前(a)后(b)CIG道集

Fig.10 The CIGs before (a) and after (b) velocity update

图11是对应于图8速度模型范围的偏移结果,显然,采用相同的成像方法,构造导向滤波下的高斯束层析反演对应的偏移结果中断点更加连续,复杂构造区信噪比更高,整体成像质量明显优于常规方法所得结果,这对后续的构造落实有较好的指导意义。通过该实际资料处理过程可以认识到,基于构造导向滤波的高斯束层析反演速度建模能够反演得到速度的中高波数分量,所得成像结果细节更加丰富,信噪比高,成像质量好。

图11

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

图11 某实际数据PSDM成像结果

a—常规高斯束层析反演对应的偏移结果;b—构造导向滤波高斯束层析反演对应的偏移结果

Fig.11 PSDM imaging results

a—the imaging results corresponding to conventional Gausssian beam tomography;b—the imaging results corresponding to structural-guided filter Gaussian beam tomography

4 结论

基于波动方程线性化近似的深度域高斯束层析速度建模方法具有较高的反演精度与实用性。在此基础上,本文引入结构张量算法,推导了构造导向滤波正则化技术。该技术通过计算地震数据的结构张量,并利用各向异性光滑算子将其应用于高斯束层析反演中,作为正则化约束项,实现了构造导向滤波约束下的高斯束层析速度建模方法。相比于常规的高斯束层析速度建模技术,该方法完全基于数据驱动,实现了对层析反演的“软约束”,既缓解了层析反演求解的多解性,同时提高反演分辨率与稳定性,得到更加符合地质认识的高精度速度模型,有效改善成像效果,提高了高斯束层析反演对复杂模型的反演分辨率,为高斯束实用化推广铺平了道路。

进一步可考虑将基于结构张量的构造导向滤波算法应用于各向异性层析建模,通过更高精度且自动化地计算构造的方位角及倾角推进各向异性高斯束层析方法的实用化。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Brillatz

, Vigée

, Coléou

, et al.

Getting closer to the geology contrasts by starting stratigraphic inversion with velocity models from new tomographic methods

[C]//SEG Technical Program Expanded Abstracts, 2014: 4763-4767.