二维VTI介质qP波斜率层析方法

谷巍巍,杨 锴,王 潇

(同济大学海洋与地球科学学院反射地震学科组,上海200092)

立体层析是将一个局部相干同相轴在炮道集与检波点道集内的射线参数水平分量(psx、prx,以下简称地表观测px参数)纳入到数据空间之中,使得数据空间的信息相对传统反射层析更为丰富。除了旅行时t外,地表观测px参数与炮检点坐标也都纳入到立体层析数据空间中,使得立体层析成为层析成像方法中唯一一种可以同时反演速度、反射点位置、反射张角与构造倾角的方法[1-4]。

斜率层析属于立体层析的一种特殊情况,在利用运动学反偏移重建数据空间的过程中,会出现炮检点位置和旅行时信息已经匹配而地表px参数无法匹配的情况,这时可利用px残差完成对速度模型、反射点坐标和散射角的更新[5-6]。斜率层析算法也因此而得名。本文将上述情况下的各向同性介质斜率层析方法推广到VTI介质。与前人研究工作的不同之处在于,当VTI介质中的运动学反偏移完成之后,我们利用px残差仅反演三个各向异性参数,而不反演反射点坐标和散射角这些与地下成像点有关的局部运动学信息。这是由于各向异性参数更新之后,会立即执行VTI介质中的运动学射线追踪,以找到更新后介质中的成像点位置(即运动学偏移),运动学偏移完成后相当于更新了成像点的局部运动学信息。这样做的优点是缓解了算法实际应用中需要联合反演多种物理参数的压力。二维典型理论数据试验结果表明,仅利用地表观测的px残差确实可以完成VTI介质中三个各向异性参数的估计,将其与VTI介质中的运动学偏移/反偏移相结合,是一种非常实用的各向异性参数估算策略。本文考虑到qP波依然是大部分实际地震采集中需要面对的主要波动现象,基于文献[7]提出的VTI介质拟声波程函方程的汉密尔顿形式,推导了qP波斜率层析所需的Fréchet偏导数公式,给出了VTI介质中的运动学反偏移与VTI介质qP波斜率层析相结合的处理流程。需要强调的是,本文所讨论的斜率层析只用到了px参数关于各向异性三参数的偏导数,因此地表观测位置(sx,rx)和旅行时t关于地下各向异性介质的偏导数将被略去。

1 二维VTI介质qP波斜率层析偏导数推导

任何一种层析反演算法都必须建立数据空间残差与模型空间残差之间的线性关系,即:

Δd=FΔm

(1)

基于方程(1),利用观测到的数据残差就可以计算出模型空间的更新量,达到更新初始模型的目的。其关键步骤是建立方程(1)所示的Fréchet偏导数矩阵(以下简称F矩阵)。注意,在二维VTI介质qP波立体层析中,数据空间为d=[sxrxtpsxprx],其中sx、rx为炮检点坐标;t为走时;psx、prx为炮检点处的射线参数水平分量。地下模型空间为m=[x0z0θsθrvvεδ],其中x0、z0为地下反射点x的坐标;θs、θr分别代表从炮点一侧与检波点一侧出射的散射相角;vv、ε、δ代表VTI介质中的三个各向异性参数。

(2)

方程(2)展示了二维VTI立体层析Fréchet偏导数矩阵,其中第一行、第二行、第三行分别为旅行时t、炮点坐标sx、检波点坐标rx关于立体层析模型空间各分量的偏导数,本文不予讨论。考虑到本文讨论的斜率层析仅限于反演各向异性参数,这里仅关注炮点射线参数水平分量psx、检波点射线参数水平分量prx关于三个各向异性参数的偏导数,关于散射相角和反射点位置的偏导数将不予讨论。因此,本文使用的VTI介质qP波斜率层析Fréchet偏导数矩阵退化为:

(3)

注意,公式(2)、公式(3)中的κ为针对矩阵中每一行的数量级不同设置的均衡系数,目的是保证每一行的数值范围基本一致。

WANG等[8]基于重新参数化后的二维VTI介质拟声波程函方程和射线扰动理论[9],推导了二维qP波立体层析Fréchet偏导数矩阵。本文采用类似的方式推导出适用于二维qP波的斜率层析Fréchet偏导数矩阵。二维VTI介质拟声波程函方程如下[7]:

(4)

式中,vv是介质垂向速度;vnmo是NMO速度。根据文献[7]可知各向异性参数之间有如下关系:η=(ε-δ)/(1+2δ)和vnmo=vv(1+2δ)1/2。为了保证变量之间的独立性,不妨采用一种等价的程函方程,其形式如下:

(5)

因此,二维VTI介质下的qP波汉密尔顿算子可以写为:

(6)

相应地,与射线参数水平分量和垂直分量有关的射线追踪一阶微分方程组为:

(7)

式中:σ为沿射线方向的滑动积分变量。

对方程(7)做全微分,可得到由初始相空间扰动(Δpx,Δpz)及各向异性参数扰动(Δvv,Δε,Δδ)所引起的射线轨迹上的相空间扰动:

(8)

其中,

(9a)

根据射线扰动理论[9],关于模型空间3个各向异性参数的扰动汉密尔顿形式可以写为:

(9b)

(8)式中后三项形式如下:

(9c)

结合(9a)、(9b)、(9c)式,利用射线传播矩阵理论[10]即可求得方程(8)的传播矩阵解:

(10)

式中:σ0代表射线出射位置,σ1代表射线传播的最终位置,Π(σ1,σ0)表示从射线出射位置到最终位置的传播矩阵。由于本文仅利用px作为数据残差反演各向异性参数,而不反演出射位置处的散射角信息,因此(10)式可进一步简化为:

(11)

利用(11)式即可获得二维VTI各向异性介质qP波斜率层析所需的数据分量与模型空间各分量之间的一阶扰动关系,也意味着斜率层析Fréchet偏导数矩阵F得到建立,后续通过求解线性方程组(1),就可以反演各向异性参数。

2 二维VTI介质qP波斜率层析中运动学反偏移及实现流程

首先介绍文献[5]和文献[6]提出的运动学反偏移概念以及由此导出的各向同性介质中的斜率层析方法。如图1所示,v是地下成像点(x,z)处的速度,ts、tr分别是地下成像点(x,z)到炮点s、检波点r的单程走时,h是半偏移距,m是炮点s与检波点r的中点。图1a表示地下模型信息与地表数据信息的几何关系;图1b为共偏移距偏移成像剖面,在该剖面上拾取构造倾角ξ;图1c为偏移距域共成像点道集(CIG),在该道集上拾取剩余曲率(RMO)的斜率信息tanφ。基于偏移后共偏移距成像剖面上搜索到的构造倾角ξ,向地表发射线找到对应的炮检点坐标,对应的射线水平参数分量设为psx和prx。

图1 运动学反偏移原理(据文献[5])a 地下模型信息与地表数据信息的几何关系; b 偏移后共偏移距剖面; c 偏移后共成像点道集

(12a)

(12b)

式中,α=2vcosβcosξ,β是反射角,它与两个散射角θs、θr以及构造倾角ξ的关系如下:

cosθs+cosθr=2cosβcosξ

(12c)

将上述各向同性介质中的运动学反偏移推广到VTI介质非常简单,只需要将各向同性介质的运动学射线追踪更换为VTI介质的运动学射线追踪,将各向同性克希霍夫叠前深度偏移更换为VTI介质克希霍夫叠前深度偏移即可。

图2 VTI介质的运动学反偏移及斜率层析流程

3 理论数据测试

利用二维VTI多层背斜模型数据,验证了本文计算公式的正确性和实现流程的可靠性。图3为6层背斜VTI参数模型,基于该模型进行VTI波动方程正演模拟。正演数据共601炮,每炮801道,炮间距20m,道间距10m,最大偏移距4000m,记录时间为4s。图4展示了正演数据的3个单炮记录。大量实际资料分析表明,δ参数通常是地震响应较弱的一个参数,目前工业界主要通过利用井资料对各向同性深度偏移成像剖面进行标定后,再利用真实深度与成像深度的差转换得到δ剖面,这是一个比较成熟的技术[13]。因此,本文重点讨论如何利用提取出的地表px参数残差反演垂直速度vv和ε,将δ剖面设置为正确值,不对其进行反演。

图5展示了初始vv剖面和初始ε剖面(常梯度)。为了检验本文算法联合反演这两个参数的能力,这里使用的初始vv各向同性模型为经过倾角时差校正(DMO)速度分析后转到深度域、并做了大尺度平滑后的速度模型。

基于上述两个剖面以及正确的δ剖面,实施第1轮VTI介质克希霍夫深度偏移,获得初始小偏移距叠加成像剖面(偏移距范围0～500m)和若干初始共成像点道集(图6)。基于图6a所示的初始成像剖面,拾取了700多个数据点(图6b),然后从每个数据点出发,利用结构张量算法从零偏移距到最大偏移距每隔5个偏移距向地表发出一对射线,这样得到的数据空间一共有1.4×104个射线对。图6d 展示了1500m偏移距的偏移成像剖面局部放大结果,以及CDP1430处的CIG,图中水平方向红线对应深度为820m。两条红线的交点为A,利用结构张量算法可以快速地在1500m偏移距成像剖面内搜索出关于A点的构造倾角ξ,同时在CDP1430处的CIG内搜索出RMO曲线的斜率信息tanφ。

图3 6层背斜VTI模型参数a 垂直速度vv剖面; b ε剖面; c δ剖面

图4 二维VTI介质波动方程模拟的3个单炮记录

图5 初始垂直速度vv(a)与初始ε(b)剖面

图6 VTI介质中通过运动学反偏移获取数据空间的实现过程a 初始小偏移距叠加成像剖面(偏移距范围0～500m); b 基于图6a拾取的数据点位置; c 利用图5所示初始vv、初始ε以及正确的δ剖面进行VTI克希霍夫叠前深度偏移获得的共成像点道集; d 1500m共偏移距成像剖面局部放大结果及CDP1430处的CIG道集; e 坐标为7100m的共炮点记录; f 坐标为8600m的共检波点记录

图9、图10分别为基于初始模型和最终反演模型的CIG道集和VTI叠前深度偏移结果。图9a对比了5km处基于初始模型和最终反演模型得到的CIG道集,图9b 对比了10km处基于初始模型和最终反演模型得到的CIG道集,可见基于最终反演模型获得的CIG均得到了有效拉平。图10对比了基于初始模型和最终模型的PSDM偏移剖面,可以看出中深部的成像品质大为提升,说明本文方法的反演结果能够满足叠前深度偏移成像的要求。

图7 联合反演结果a vv剖面; b ε剖面; c 误差泛函下降曲线

图8 不同水平位置vv和ε抽线对比a 5km处vv抽线对比; b 10km处vv抽线对比; c 5km处ε抽线对比; d 10km处ε抽线对比

图9 5km(a)和10km(b)处的参数模型更新前后CIG道集对比

图10 基于初始模型(a)和最终模型(b)PSDM偏移剖面对比

4 结论和讨论

本文将各向同性介质中的斜率层析方法推广到二维VTI介质qP波的情况。基于射线扰动理论推导了地表px参数关于地下三个各向异性参数的偏导数,建立了适合于二维VTI介质qP波斜率层析的线性方程组,并将各向同性介质中的运动学偏移/反偏移策略推广到VTI介质中,最终实现了二维VTI介质qP波斜率层析方法。研究结果表明,VTI介质中的运动学偏移/反偏移与二维VTI介质qP波斜率层析线性方程组相结合,可以很好地实现二维VTI介质中的各向异性参数重建。该方法也可以方便地推广到三维或者TTI介质。

需要指出的是,无论是各向同性介质还是各向异性介质中的斜率层析方法,对初始模型都有一定要求:初始成像体的品质能够保证CIG上提取的RMO斜率tanφ和共偏移距成像剖面上提取的构造倾角ξ都是可靠的。如果这个要求无法满足,那么后续通过运动学反偏移获取的地表px参数精度将难以保证。

其次,在斜率层析中忽略旅行时对反射点位置和散射角的导数项不会降低算法的适用性。这是因为,如果在立体层析反演中能够将数据残差通过合理的方式转移到某一个数据分量上,则会使得数据空间得到合理的压缩,这对提高反演的稳定性是有益的。这一点已经在前人发展的斜率层析和控制定向层析方法中得到了证实。同时,本文方法没有将反射点位置和散射角纳入到模型空间内,实际上又压缩了传统斜率层析的模型空间,进一步提高了反演的稳定性,有利于本文方法的推广应用。