石油地球物理勘探  2020, Vol. 55 Issue (6): 1305-1311  DOI: 10.13810/j.cnki.issn.1000-7210.2020.06.016
0
文章快速检索     高级检索

引用本文 

段心标, 王华忠, 邓光校. 应用反射理论的最小二乘逆时偏移成像. 石油地球物理勘探, 2020, 55(6): 1305-1311. DOI: 10.13810/j.cnki.issn.1000-7210.2020.06.016.
DUAN Xinbiao, WANG Huazhong, DENG Guangxiao. Least-squares reverse-time migration based on reflection theory. Oil Geophysical Prospecting, 2020, 55(6): 1305-1311. DOI: 10.13810/j.cnki.issn.1000-7210.2020.06.016.

本项研究受国家重点研发计划“深层含盐层系复杂波场高精度地震成像”(2017YFC0602804-02)和国家科技重大专项“缝洞储集体高精度地震成像技术研究”(2016ZX05014-001-002)联合资助

作者简介

段心标  副研究员, 在读博士, 1982年生; 2005年获中国石油大学(华东)勘查技术与工程专业学士学位, 2008年获中国科学院地质与地球物理所固体地球物理学硕士学位; 现就职于中国石化石油物探技术研究院, 主要从事地震偏移成像技术研发及推广应用

段心标, 江苏省南京市江宁区上高路219号中国石化石油物探技术研究院软件所, 211103。Email:duanxb.swty@sinopec.com

文章历史

本文于2020年2月14日收到,最终修改稿于同年9月2日收到
应用反射理论的最小二乘逆时偏移成像
段心标①② , 王华忠 , 邓光校     
1 同济大学海洋与地球科学学院波现象与智能反演成像研究组, 上海 200092;
2 中国石化石油物探技术研究院, 江苏南京 211103;
3 中国石化西北油田分公司, 新疆乌鲁木齐 830011
摘要:最小二乘偏移是在线性反演理论下求解模型空间的精确解,相对于常规偏移具有更高成像分辨率、保幅性,且可减少成像假象。经典最小二乘偏移以散射波线性化表达为基础,成像实质是估计介质的散射强度。由于地下实际地层以层状介质为主,反射成像是实际应用中最小二乘偏移技术的目标。文中重点讨论散射理论与反射理论两种线性化表达方法的不同,根据实际应用中偏移成像的需求,选择忽略角度信息的反射理论线性化表达反偏移方法,建立了一套估计平均反射系数的最小二乘逆时偏移流程。Sigsbee2a模型验证了最小二乘逆时偏移方法的成像效果;实际三维探区应用结果表明,与常规逆时偏移方法相比,该方法“串珠”成像收敛效果更好,对深层大断裂、层间小断层的刻画更清晰。
关键词最小二乘偏移    散射理论    反射理论    反射系数    串珠成像    
Least-squares reverse-time migration based on reflection theory
DUAN Xinbiao①② , WANG Huazhong , DENG Guangxiao     
1 Wave Phenomena and Intelligent Inversion Imaging Group(WPI), School of Ocean and Earth Sciences, Tongji University Shanghai 200092, China;
2 Sinopec Geophysical Research Institute, Nanjing, Jiangsu 211103, China;
3 Sinopec Northwest Company, Urumqi, Xiujiang 830011, China
Abstract: Least-squares migration can provide an accurate solution to the model space based on the linear inversion theory.Compared to conventional migration methods, least-squares migration can increase imaging resolution, improve amplitude preservation and reduce migration artifacts.The classical least-squares migration is based on the linearized expression of the scattered wave, and its goal is essentially estimating the scattering intensity of underground medium.While real underground medium is mainly layered, the target of least-squares migration in production is reflection-coefficient imaging.This paper focuses on the difference between linearized expression method of scattering theory and reflection theory.According to the requirement of migration imaging in actual production, then the demigration method based on linearized expression of reflection theory and the least-squares RTM process for estimating the reflection coefficient are established.Sigsbee2a model test shows a good imaging effect using the proposed method.An application of three-dimensional field data shows that the least-squares RTM method is better than RTM.The string beads imaging has better convergence effect, and the deep large faults and small interbed faults are more clearly described.
Keywords: least-squares migration    scattering theory    reflection theory    reflection coefficient    string beads imaging    
0 引言

从地震数据中准确提取并利用复杂的地下岩性与储层信息,是油气勘探的关键环节。地震波成像是连接地震观测数据与储层解释的桥梁,常规偏移成像技术[1-5]主要目的是获得能精细描述地质构造的地震成果资料,用于部署和落实井位。

随着油气勘探程度的日趋深入和地震成像技术的不断进步,直接估算地下弹性参数的反演成像方法受到越来越多的关注。早在20世纪60年代,Backus等[6-7]从数学角度研究了地球物理反问题,奠定了早期地球物理反演的理论基础。Tarantola[8-9]基于Bayes估计理论,建立了一套地震波全波形反演的理论方法体系,推动了地震反演成像技术的快速发展。然而,全波形反演具有很强的非线性特征,难以由观测数据直接反演得到较准确的介质弹性参数,无法满足油藏描述的精度需求。实际中应用的依然是经典的地震波成像处理和储层描述流程,而如何获得更高偏移成像质量是其中重要环节。

最小二乘偏移[10-18]是常规偏移的后续发展,其核心思想是在线性反演理论下求解模型空间的精确解,相当于在常规偏移构造成像基础上,考虑了振幅补偿和波形校正,从而得到更保幅和更高分辨率的成像结果。

经典最小二乘偏移以散射波线性化表达为基础,认为模型参数是速度(或慢度)的高波数扰动,其成像实质是估计介质的散射强度。Zhang等[19-20]提出基于逆时偏移叠加剖面的反偏移方法,并形成相应的最小二乘偏移技术,但并未给出反偏移方法的详细理论依据。陈生昌等[21-22]提出地震反射波数据偏移成像方法和相应的反射波动方程最小二乘偏移方法,并引入“反射率”概念,即最小二乘偏移的实质是估计反射面上的反射率。本文认为:因地下介质以层状为主,获取高保真反射系数是地震成像走向储层刻画的关键,估算不同角度的反射系数是最小二乘偏移的目标。然而,考虑到反演不同角度反射系数的计算量巨大,且高质量偏移叠加成像依然是当今实际应用的主要目标,因此以最小二乘偏移获得忽略角度信息的平均反射系数具有现实意义。

本文首先简介线性最小二乘偏移方法基础理论,重点讨论散射理论线性化表达方法和反射理论线性化表达方法,并分别建立散射强度、反射系数与地震数据的线性关系式;然后从偏移成像实际应用需求出发,选择基于反射理论并忽略角度信息的Kirchhoff正演反偏移表达方法,实现与实际应用需求相一致的最小二乘逆时偏移技术;最后,通过模型数据和实际资料验证本文方法的效果。

1 方法原理 1.1 线性最小二乘偏移方法

基于Bayes估计理论,地震反演成像可表示为后验概率密度最大化问题;在地震反演通常假设条件下,进而可转化为如下目标函数最小化问题

$ E(\mathit{\boldsymbol{m}}) = \frac{1}{2}{[L(\mathit{\boldsymbol{m}}) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}]^{\rm{T}}}[L(\mathit{\boldsymbol{m}}) - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}] $ (1)

式中:m为地下介质参数,如速度、散射强度、反射系数或波阻抗等;dobs为观测数据;L(·)表示波场传播算子。该式非线性反演问题求解难度大,需对非线性反演问题逐步进行线性化,依次反演参数的低、中、高波数信息。最小二乘偏移是一种线性化的地震反演成像方法,用于估计介质参数的高波数成分,其目标函数可表示为

$ E(\mathit{\boldsymbol{m}}) = \frac{1}{2}{\left\| {\mathit{\boldsymbol{Lm}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|^2} $ (2)

式中L为矩阵形式的线性化波场传播算子。

令目标函数梯度为零,则有如下的法方程

$ {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{Lm}} = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $ (3)

最小二乘偏移成像结果表示为

$ \mathit{\boldsymbol{\hat m}} = {({\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}})^{ - 1}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $ (4)

式中:LTL为Hessian矩阵,可用H表示;LTdobs是对观测数据进行常规偏移成像。最小二乘偏移是Hessian矩阵的逆H-1作用于常规偏移成像,能够消除Hessian矩阵导致的成像模糊化效应,提高成像分辨率和保幅性,减少成像假象。

由于“目标函数梯度为零”仅是一种理想假设,且Hessian矩阵非常庞大,其计算和存储都很困难,Hessian矩阵的求逆更难以完成,因而无法直接求解法方程以实现最小二乘偏移。为此,通常采用迭代法求解最小二乘偏移问题,假设目标函数在初始值附近满足二次型,从初始解逐步逼近真实解。

1.2 地震波场的线性化表达 1.2.1 散射理论的Born正演线性化表达

根据散射理论,地下空间位置x处实际介质速度v(x)分解为两部分:背景参考速度v0(x)和变化的扰动速度Δv(x)。相应地,总波场$\hat{u}$(x, ω)分解为参考波场(或入射波场)$\hat{u}_{0}$(x, ω)和散射波场Δ$\hat{u}$(x, ω)。

一阶Born近似假设下,散射波场可写成

$ \begin{array}{*{20}{l}} {\Delta \hat u(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega )}\\ {\quad = {\omega ^2}\int {{G_0}} (\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega )m({\mathit{\boldsymbol{x}}^\prime }){{\hat u}_0}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ){\rm{d}}{\mathit{\boldsymbol{x}}^\prime }} \end{array} $ (5)

式中:xs表示震源点;G0(x, x′, ω)表示背景介质下的Green函数;m(x′)=$\frac{1}{{{v^2}\left({{\mathit{\boldsymbol{x}}^\prime }} \right)}} - \frac{1}{{v_0^2\left({{\mathit{\boldsymbol{x}}^\prime }} \right)}}$表示空间某点的慢度平方的扰动。

对于带限子波$\hat{f}$(ω),相应入射地震波场为

$ {\hat u_0}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) = \hat f(\omega ){G_0}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) $ (6)

进而,散射波场可表示为

$ \begin{array}{l} \Delta \hat u(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) = {\omega ^2}\hat f(\omega )\int {{G_0}} (\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega )m({\mathit{\boldsymbol{x}}^\prime }) \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {G_0}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ){\rm{d}}{\mathit{\boldsymbol{x}}^\prime } \end{array} $ (7)

式(7)定义了一个Fredholm第一类积分方程,它描述了一个线性问题。其积分核可写成

$ \mathit{\boldsymbol{L}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) = {\omega ^2}\hat f(\omega )\sum\limits_{{\mathit{\boldsymbol{x}}^\prime }} {{G_0}} (\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega ){G_0}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) $ (8)

在散射假设下,可认为地表接收到的携带地下信息的信号就是散射波场记录,则可将式(7)写成如下线性算子表达式

$ \mathit{\boldsymbol{Lm}} = {\mathit{\boldsymbol{d}}_{{\rm{obs}}}} $ (9)

此式中的m设定为慢度平方扰动,则该式表达了散射强度与地震数据之间的线性关系。

1.2.2 反射理论的Kirchoff正演线性化表达

反射波理论假设下,地下界面Ω上入射波场引起的二次场可表示为

$ {\hat u_{\rm{s}}} = R\hat f(\omega ){G_{\rm{s}}} $ (10)

式中:R是反射系数;Gs是震源点xs到反射面的Green函数;$\hat{f}$(ω)Gs是入射波场。式(10)的描述更符合物理实际。

地下界面Ω上波场的法向导数为

$ \frac{{\partial {{\hat u}_{\rm{s}}}}}{{\partial n}} = - R\hat f(\omega )\frac{{\partial {G_{\rm{s}}}}}{{\partial n}} $ (11)

利用WKBJ近似方法[23],有

$ \frac{{\partial {G_{\rm{s}}}}}{{\partial n}} \approx {\rm{i}}\omega \mathit{\boldsymbol{n}} \cdot \nabla {\tau _{\rm{s}}} \cdot {G_{\rm{s}}} $ (12)

式中:τs是WKBJ近似下的波的走时;n·τs表达了波传播方向与地面法线方向之间的夹角关系,控制了不同方向上波的幅值。

由Green定理,反射面上方地面上任意一点x的反射波场可写成如下积分式

$ \begin{array}{l} {{\hat u}_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega } \right) = \hat f(\omega )\int_\varOmega {\left[ {{{\hat u}_{\rm{s}}}\left( {{\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega } \right)\frac{{\partial {G_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega } \right)}}{{\partial n}} - } \right.} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {{G_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega } \right)\frac{{\partial {{\hat u}_{\rm{s}}}\left( {{\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega } \right)}}{{\partial n}}} \right]{\rm{d}}\varOmega \end{array} $ (13)

式中Gr(x, x′, ω)是反射波传播的Green函数(图 1)。

图 1 波场传播示意图

将式(10)~式(12)代入式(13),可得

$ \left\{ \begin{array}{l} {{\hat u}_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) \approx {\rm{i}}\omega \hat f(\omega )\int_\varOmega R \mathit{\boldsymbol{n}}(\nabla {\tau _{\rm{s}}} + \nabla {\tau _{\rm{r}}}) \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {G_{\rm{s}}}({x^\prime },{x_{\rm{s}}},\omega ){G_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega ){\rm{d}}\varOmega \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {\rm{i}}\omega \hat f(\omega )\int_\varOmega R \mathit{\boldsymbol{n}} \cdot \nabla \tau \cdot {G_{\rm{s}}}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {G_{\rm{r}}}({\mathit{\boldsymbol{x}}_r},{\mathit{\boldsymbol{x}}^\prime },\omega ){\rm{d}}\varOmega \\ \tau = {\tau _{\rm{s}}} + {\tau _{\rm{r}}} \end{array} \right. $ (14)

根据程函方程关系,式(14)可进一步写成

$ \begin{array}{*{20}{c}} {{{\hat u}_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) \approx - {\rm{i}}\omega \hat f(\omega )\int_\varOmega R \mathit{\boldsymbol{n}} \cdot \mathit{\boldsymbol{s}}({\mathit{\boldsymbol{x}}^\prime }) \times }\\ {{G_{\rm{s}}}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ){G_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega ){\rm{d}}\varOmega } \end{array} $ (15)

式中s(x′)为慢度矢量。

式(15)建立了地震数据与不同方向反射系数之间的联系[24-25]。然而,在缓变速度结构中,不同角度的反射波能量较均衡;在速度剧变结构中,不同角度的反射系数难以估计。同时,考虑到直接基于式(15)进行波场模拟的计算量巨大,因此建立地震数据与平均反射系数之间的关系就具有现实意义。若不考虑背景速度对振幅的影响,式(15)就退化为

$ \begin{array}{*{20}{l}} {{{\hat u}_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) \approx - {\rm{i}}\omega \hat f(\omega ) \times }\\ {\quad \int_D {{R^\prime }} {G_{\rm{s}}}({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ){G_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega ){\rm{d}}{\mathit{\boldsymbol{x}}^\prime }} \end{array} $ (16)

式中:D表示散射体;R′表示平均反射系数。

在反射假设下,认为地表接收到的携带地下信息的信号就是反射波场记录。式(16)同样可写成线性算子(式(9))表达形式,此时m表示忽略地下角度信息的平均反射系数,其积分核为

$ \mathit{\boldsymbol{L}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ) = - {\rm{i}}\omega \hat f(\omega )\sum\limits_{{\mathit{\boldsymbol{x}}^\prime }} {{G_{\rm{s}}}} ({\mathit{\boldsymbol{x}}^\prime },{\mathit{\boldsymbol{x}}_{\rm{s}}},\omega ){G_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}^\prime },\omega ) $ (17)
1.3 最小二乘逆时偏移实现方法

逆时偏移技术现已在实际处理中得到广泛应用,其成像基础是Claerbout[26]给出的成像方法:波场外推+成像条件。成像条件为:成像点处,入射波的到达时等于反射波的出发时。显然,常规逆时偏移是估计地下介质的反射系数。由于地下介质以层状为主,以定位反射界面的位置或产生角度反射系数道集为目标的反射系数成像是当前地震成像的主流需求。

如前所述,由于角度反射系数反演的计算量巨大,因而实际应用中,最小二乘偏移的目的是在常规逆时偏移基础上进一步提高偏移叠加成像质量,使成像结果分辨率更高,成像假象更少。

在反演成像方法中,反偏移采用Kirchhoff正演思路,具体迭代反演步骤如下。

(1) 以设定的偏移速度场对观测数据dobs作逆时偏移成像,作为最小二乘偏移的初始值m0

(2) 对第i次迭代后偏移成像值mi(i=0, 1, 2, …)做反偏移,正演模拟预测数据dcal,反偏移公式为

$ \left\{ \begin{array}{l} \frac{1}{{v_0^2}}\frac{{{\partial ^2}{u_{\rm{s}}}}}{{\partial {t^2}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},t) - {\nabla ^2}{u_{\rm{s}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},t)\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = f(t)\delta (\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}})\\ \frac{1}{{v_0^2}}\frac{{{\partial ^2}{u_{\rm{r}}}}}{{\partial {t^2}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},t) - {\nabla ^2}{u_{\rm{r}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},t)\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = m(\mathit{\boldsymbol{x}})\frac{{\partial {u_{\rm{s}}}(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_{\rm{s}}},t)}}{{\partial t}}\\ {\mathit{\boldsymbol{d}}_{{\rm{cal}}}}({\mathit{\boldsymbol{x}}_r},{\mathit{\boldsymbol{x}}_{\rm{s}}},t) = {u_{\rm{r}}}({\mathit{\boldsymbol{x}}_r},{\mathit{\boldsymbol{x}}_{\rm{s}}},t) \end{array} \right. $ (18)

(3) 计算观测数据dobs与预测数据dcal的残差波场δd,当残差满足设定条件时迭代停止。

(4) 对δd进行逆时偏移成像,得到目标泛函的梯度gi

(5) 计算共轭梯度方向

$ {\mathit{\boldsymbol{h}}_k} = \left\{ {\begin{array}{*{20}{c}} { - {\mathit{\boldsymbol{g}}_k}}&{k = 1}\\ { - {\mathit{\boldsymbol{g}}_k} + {\alpha _k}{\mathit{\boldsymbol{h}}_{k - 1}}}&{k \ge 2} \end{array}} \right. $ (19)

式中αk=$\frac{\left(\boldsymbol{g}_{k}-\boldsymbol{g}_{k-1}\right)^{\mathrm{T}} \boldsymbol{g}_{k}}{\boldsymbol{g}_{k-1}^{\mathrm{T}} \boldsymbol{g}_{k-1}}$。

(6) 计算迭代步长

$ {\mu _{k + 1}} = \frac{{\mathit{\boldsymbol{h}}_k^{\rm{T}}{\mathit{\boldsymbol{g}}_k}}}{{{{(\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{h}}_k})}^{\rm{T}}}\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{h}}_k}}} $ (20)

(7) 更新偏移成像值

$ {\mathit{\boldsymbol{m}}_{k + 1}} = {\mathit{\boldsymbol{m}}_k} - {\mu _{k + 1}}{\mathit{\boldsymbol{h}}_k} $ (21)

(8) 转到步骤(2),展开新一轮迭代。

由于迭代算法计算量巨大,三维最小二乘逆时偏移中采用了GPU加速计算技术,在此不再赘述。

2 数值模型试验及实际应用

首先利用二维Sigsbee2a模型(图 2a)做数值试验,验证本文方法效果;进而将三维最小二乘逆时偏移技术应用于实际地震资料成像处理。

图 2 Sigsbee2a模型及其成像结果 (a)速度模型;(b)常规逆时偏移;(c)最小二乘逆时偏移
2.1 Sigsbee2a模型数据

为了便于计算,把Sigsbee2a经典模型速度值单位转为m/s,并重新定义速度网格,新网格的尺寸为3200×1200,间隔为5m×5m。对该模型正演模拟120炮数据,炮间距为100m,检波点间隔为20m,炮检距范围是-2000~2000m,记录长度为7s。

利用本文所提最小二乘逆时偏移方法进行反演迭代,得到第20次迭代结果(图 2c),并与常规逆时偏移剖面(图 2b)进行对比;图 3是目标泛函迭代误差分析曲线,可见反演逐步收敛。从整体剖面及盐丘左下局部放大图(图 4)的对比可见,相对于常规逆时偏移方法,本文方法的成像效果(图 4b)显著改善,分辨率变高,深部成像更清晰。

图 3 迭代误差曲线

图 4 Sigsbee2a模型的常规逆时偏移(a)和最小二乘逆时偏移(b)成像的局部放大
2.2 实际资料

本次试验所用实际数据取自中国西北M油田,测试处理的目的是提高储层的成像精度,力求对奥陶系缝洞体及小断裂精确成像。所选数据满覆盖面积为401km2,共48000炮,成像面元是30m×30m。利用100GPU节点(每节点4块K10卡)3次迭代,共用时59天。

图 5图 6是两条不同测线的常规逆时偏移与最小二乘逆时偏移成像剖面对比,可见最小二乘逆时偏移剖面中“串珠”成像收敛更好,对深层大断裂、层间小断层的刻画更清晰,T78(奥陶系蓬莱坝组顶界面)以下成像能量更强,成像分辨率得到提高。对比奥陶系目的层的成像频谱(图 7),可见最小二乘逆时偏移成像拓宽了频谱(图 7b,频带展宽5Hz);T74(一间房组顶界面)下40~120ms均方根振幅,也证实本文方法成像结果(图 8b)对“串珠”刻画更清晰。

图 5 L1测线的常规逆时偏移(a)和最小二乘逆时偏移(b)成像剖面

图 6 L2测线的常规逆时偏移(a)和最小二乘逆时偏移(b)成像剖面

图 7 常规逆时偏移(蓝线)和最小二乘逆时偏移(红线)结果的目的层频谱

图 8 常规逆时偏移(a)与最小二乘逆时偏移(b)的均方根振幅平面属性对比

以上结果均展示了最小二乘逆时偏移成像在提高成像分辨率、振幅均衡性及成像精度等方面的效果。

然而也须指出,最小二乘逆时偏移实际应用效果与其理论优势尚有一段距离。近年国外地球物理公司虽进行了大量的最小二乘逆时偏移应用尝试,但整体而言现今该方法的应用仍处于试验探索阶段。

3 结论与讨论

最小二乘偏移是一种线性化的地震反演成像方法,线性近似可以是“背景速度+速度扰动”模式,也可以是“背景速度+反射系数”模式。由于地下介质以层状为主,现今业界常规偏移成像的主要目标是反射系数。本文分析了两种线性化模式的差异,选择基于忽略角度信息的反射系数的反偏移正问题表达方法,提出一种估计平均反射系数的最小二乘逆时偏移方法。理论模型测试结果表明本方法能提高成像分辨率,深部成像更清晰。将研发的三维最小二乘逆时偏移技术应用于中国西部探区资料处理,所得“串珠”成像比常规逆时偏移结果收敛更好,深层大断裂、层间小断层的刻画更清晰,成像分辨率获得提高,充分体现了最小二乘逆时偏移成像在提高成像分辨率和振幅均衡性及提高成像精度方面的优势。

然而,目前最小二乘逆时偏移实际应用效果仍不甚理想,主要原因包括反偏移得到的无噪声一次波模拟数据与实际观测记录差异较大、实际地震子波未知、现有建模精度达不到最小二乘逆时偏移要求及三维技术计算量大等。为了降低最小二乘逆时偏移对速度模型和地震子波的依赖,近年也发展了动态拉伸时差校正方法和空变子波估计方法,一定程度上减弱了反演的病态性。但考虑到反问题的复杂性,这些应用瓶颈问题并未得到有效解决,它们会持续伴随最小二乘逆时偏移的发展过程,有待逐步克服和解决。

参考文献
[1]
Baysal E, Kosloff D D, Sherwood J W C. Reverse time migration[J]. Geophysics, 1983, 48(11): 1514-1524. DOI:10.1190/1.1441434
[2]
Whitmore N D.Iterative depth migration by backward time propagation[C].SEG Technical Program Expanded Abstracts, 1983, 2: 382-385.
[3]
Lumley D, Claerbout J, Bevc D.Anti-aliased Kirchhoff 3D migration[C].SEG Technical Program Expanded Abstracts, 1994, 13: 1282-1285.
[4]
Ristow D, Ruhl T. Fourier finite-difference migration[J]. Geophysics, 1994, 59(12): 1882-1893. DOI:10.1190/1.1443575
[5]
Sun J, Abma R, Bernitsas N.Anti-aliasing in Kirchhoff migration[C].SEG Technical Program Expanded Abstracts, 1999, 18: 1134-1137.
[6]
Backus G E, Gilbert F. Numerical applications of a formalism for geophysical inverse problems[J]. Geophysical Journal of the Royal Astronomical Society, 1967, 13(2): 247-276.
[7]
Backus G E, Gilbert J F. Uniqueness in the inversion of inaccurate gross earth data[J]. Philosophical Transactions of the Royal Society B:Biological Sciences, 1970, 266(1): 123-192.
[8]
Tarantola A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754
[9]
Tarantola A.Inverse Problem Theory and Methods for Model Parameter Estimation[M].Society for Industrial and Applied Mathematics, 2005.
[10]
Cole Sand Karrenbach M.Least-Squares Kirchhoff Migration[R].Stanford Exploration Project Report, 1992.
[11]
Nemeth T, Wu C J and Schuster T. Least squares migration incomplete reflection data[J]. Geophysics, 1999, 64(1): 208-221.
[12]
Tang Y.Wave-equation Hessian by phase encoding[C].SEG Technical Program Expanded Abstracts, 2008, 27: 2201-2205.
[13]
Dai W, Wang X and Schuster G T. Least-squares migration of multisource data with a deblurring filter[J]. Geophysics, 2011, 76(5): R135-R146. DOI:10.1190/geo2010-0159.1
[14]
郭振波, 李振春. 最小平方逆时偏移真振幅成像[J]. 石油地球物理勘探, 2014, 49(1): 113-120.
GUO Zhenbo, LI Zhenchun. Least squares reverse time migration for true amplitude imaging[J]. Oil Geophysical Prospecting, 2014, 49(1): 113-120.
[15]
李闯, 黄建平, 李振春, 等. 预条件最小二乘逆时偏移方法[J]. 石油地球物理勘探, 2016, 51(3): 513-520.
LI Chuang, HUANG Jianping, LI Zhenchun, et al. Preconditioned least squares reverse time migration[J]. Oil Geophysical Prospecting, 2016, 51(3): 513-520.
[16]
匡斌, 唐祥功, 张猛, 等. 三维自适应最小二乘逆时偏移技术[J]. 石油地球物理勘探, 2018, 53(1): 73-79.
KUANG Bin, TANG Xianggong, ZHANG Meng, et al. 3D adaptive least-squares reverse-time migration[J]. Oil Geophysical Prospecting, 2018, 53(1): 73-79.
[17]
王华忠, 胡江涛, 郭颂. 最小二乘叠前深度偏移成像理论与方法[J]. 石油物探, 2017, 56(2): 159-170.
WANG Huazhong, HU Jiangtao, GUO Song. Theory and method of least square prestack depth migration and imaging[J]. Geophysical Prospecting for Petro-leum, 2017, 56(2): 159-170.
[18]
姚振岸, 孙成禹, 喻志超, 等. 融合点扩散函数的预条件黏声最小二乘逆时偏移[J]. 石油地球物理勘探, 2019, 54(1): 73-83.
YAO Zhen'an, SUN Chengyu, YU Zhichao, et al. Preconditioned visco-acoustic least squares reverse time migration integrated with point spread function[J]. Oil Geophysical Prospecting, 2019, 54(1): 73-83.
[19]
Zhang Y, Duan L, Xie Y.A stable and practical implementation of least-squares reverse time migration[C].SEG Technical Program Expanded Abstracts, 2013, 32: 3716-3720.
[20]
Zhang Y, Duan L, Xie Y. A stable and practical implementation of least-squares reverse time migration[J]. Geophysics, 2015, 80(1): V23-V31.
[21]
陈生昌, 周华敏. 再论地震数据偏移成像[J]. 地球物理学报, 2016, 59(2): 643-654.
CHEN Shengchang, ZHOU Huamin. Re-exploration to migration of seismic data[J]. Chinese Journal of Geophysics, 2016, 59(2): 643-654.
[22]
陈生昌, 周华敏. 地震数据的反射波动方程最小二乘偏移[J]. 地球物理学报, 2018, 61(4): 1413-1420.
CHEN Shengchang, ZHOU Huamin. Least squares migration of seismic data[J]. Chinese Journal of Geophysics, 2018, 61(4): 1413-1420.
[23]
Clayton R W, Stolt R H. A Born-WKBJ inversion me-thod for acoustic reflection data[J]. Geophysics, 1981, 46(11): 1559-1567. DOI:10.1190/1.1441162
[24]
Jaramillo H H, Bleistein N. The link of Kirchhoff migration and demigration to Kirchhoff and Born mode-ling[J]. Geophysics, 1999, 64(6): 1793-1805. DOI:10.1190/1.1444685
[25]
Bleistein N, Zhang Y, Xu S. Migration/inversion:think image point coordinates, process in acquisition surface coordinates[J]. Inverse Problems, 2005, 21(5): 1715-1744. DOI:10.1088/0266-5611/21/5/013
[26]
Claerbout J F. Towards a unified theory of reflector mapping[J]. Geophyscis, 1971, 36(3): 467-481.