基于声波测井资料的波动方程走时层析反演速度建模

李卿卿,符力耘*,杜启振,冉亚楠,胡勇

1 深层油气全国重点实验室(中国石油大学(华东)),山东青岛 266580 2 山东省油藏地质重点实验室,中国石油大学(华东),山东青岛 266580 3 青岛海洋科学与技术国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛 266071 4 中国矿业大学资源与地球科学学院,江苏徐州 221116

0 引言

地震层析反演方法由于主要利用地震波的走时信息,其非线性较全波形反演弱,因此常用于建立运动学偏移速度模型或为全波形反演提供初始参数模型.地震层析反演方法主要分为射线与波动方程两类.

目前,工业上主要采用基于射线的层析反演方法,该类方法主要分为两种:一种是将角道集(CIG道集)的剩余时差(RMO)转换为走时差(Stork,1992),通过射线追踪计算射线路径并构造相应的层析线性方程组,而层析线性方程组的求解过程可以解释为将剩余时差反投影到模型空间实现速度反演(Meng et al.,2004; Woodward et al.,2008; Jones,2010);第二种是利用程函方程(井西利等,2007;刘玉柱等,2009,2011;张新彦等,2017;蔡杰雄等,2017;张兵和王华忠,2019;董良国等,2021),通过求取初至波走时差信息,实现初至波层析反演.然而,射线走时层析技术基于高频近似假设,而地震波是一种包含多频率成分的信号.因此,射线路径往往只在高速区域产生优势采样,而其他区域反演精度不足.此外,射线层析方法无法处理多路径问题,导致层析剖面存在焦散、阴影区等问题(Luo and Schuster,1991).综上所述,虽然基于射线理论的层析方法计算效率高且操作灵活,但对于较为复杂的介质该类方法往往难以获得令人满意的反演效果.

针对传统走时层析方法高频假设存在的问题,Luo和Schuster(1991)提出了波动方程层析成像方法,即以波路径替代传统走时层析成像中的射线路径.由于波动方程走时层析成像方法不存在高频近似假设,且能够处理多路径问题.因此,基于波动方程走时层析成像方法在处理复杂介质时层析精度远优于基于射线理论的层析反演方法.根据所用地震波的类型,波动方程走时层析反演可分为初至波和反射波两类.

Luo 和 Schuster (1991)最早提出了初至波波动方程走时层析反演方法,该方法通过计算观测数据与模拟数据互相关的最大值来获取走时差(de Hoop and van der Hilst,2005),然后利用走时差伴随场计算目标函数的梯度(Tarantola,1987),避免了射线追踪和Fréchet导数的计算.由于仅利用初至波走时信息,初至波走时层析反演方法稳定性高,可实现浅地表、中深层大尺度或井间速度建模.然而,初至波走时层析反演方法仅利用了地震记录中的初至波信息,往往反演精度较低,无法对地下中深层实现较高精度速度建模.

反射波走时层析反演由于其对中深层速度进行准确建模能力而受到广泛关注.该方法也基本分为两类,即数据域反射波走时层析和成像域反射波走时层析反演.在数据域,由于反射波走时差沿时间方向变化,且反射波层析反演初始速度模型通常为均匀模型,无法直接从实际模拟地震数据中提取走时信息,继而无法获得走时层析反演最为重要的走时差伴随场.因此,数据域反射波走时层析反演往往将Born一阶近似的反偏移记录作为模拟地震反射波信息,继而利用Gaussian窗函数或者动态编码方法(DIW)获得走时差伴随场(Hale,2013; Ma and Hale,2013),最优化更新速度信息,实现反射波层析反演速度建模.在成像域,则需在成像过程中提取角道集(CIG),并依据CIG中同相轴曲率计算走时差伴随场 (Jin and Beydoun,2000; Meng et al.,2004; Zhang et al.,2011,2015a,b; Yang and Sava,2013; Alkhalifah and Choi,2014; Chi et al.,2015),再通过计算目标函数的梯度,实现成像域反射波层析反演速度建模.

层析反演已是一种较为成熟的速度建模方法.然而,传统初至波走时层析反演方法对地下介质分辨能力有限,无法获得准确的中深层速度信息.相比之下,反射波层析反演方法能明显提高中深层速度建模的精度,但数据域反射波层析反演速度建模需利用反偏移模拟数据,这无疑大大增加了计算量,降低了计算效率;同时,基于Born近似的反偏移过程是一种线性近似,其模拟数据与实际地震资料也存在一定差距.此外,基于成像域反射波层析反演速度建模方法需从CIG中提取走时差信息,也同样面临计算量较大、实用性较差等问题.为此,一些学者提出了融合初至波与反射波的走时层析反演方法(张兵和王华忠,2019),但如何确定两者之间的权重也是一个不可忽视的难题.

声波测井资料是地球物理勘探中最为重要且精确的数据之一,其含有丰富的低频与高频成分,具有垂向分辨率高的特点.将测井信息引入初始速度模型建立可以在一定程度上弥补地震频带缺失、反射数据精度低、信噪比低的不足.由于声波测井资料具有高精度的速度及可靠的低频信息,可从目标泛函约束的角度引入先验模型信息,但该过程的操作难度较大,常用于AVO反演.吴国忱和田军(2014)通过在AVO反演过程中加入尺度约束项,提高了反演的稳定性及可靠性.冯国峰等(2003)针对二维波动方程,将大范围收敛的同伦摄动方法引入到算子参数识别过程,并利用声波测井资料进行约束反演.Han等(2007)和傅红笋等(2012)研究了测井约束下地震波形反演的同伦摄动求解方法.Asnaashari等(2013)提出通过测井、地质解释建立先验模型,并使用正则化项和先验模型项扩展反演目标函数,将先验模型用于反演优化过程,以此约束最小化目标泛函,该方法是将约束条件引入全波形反演中的经典性代表.杜泽源等(2017)研究了基于测井约束的全波形反演方法,对反演过程进行约束,提高了反演精度.然而,在反演过程中引入声波测井资料约束项的方法存在两个缺陷:测井项梯度与主梯度更新方向不一致,导致反演错误;声波测井资料约束项的权重系数无法由理论确定,而往往采用经验的做法.

本文区别于将传统声波测井资料作为目标函数约束项的方法,提出了一种基于声波测井资料的三步法波动方程全波走时层析反演速度建模方法,旨在高效地实现浅层及中深层速度建模.本文架构如下:首先,利用波动方程早至波走时层析反演方法获得初始低波数速度模型,并通过一次偏移成像获得初始成像剖面.其次,采用平面波分解从一次成像剖面中提取构造角度,并结合声波测井信息沿构造角度进行稀疏反演,建立更高精度的中波数速度剖面.最后,通过求取全波走时差,采用最优化方法更新速度,实现全波高波数层析反演速度建模.Marmousi模型数值结果表明,本文提出的层析反演速度建模方法能够有效提高浅层及中深层速度建模精度.

1 方法原理

1.1 波动方程早至波走时层析反演速度建模方法简介

波动方程早至波走时层析与射线走时层析的原理类似,均基于速度扰动仅对地震波走时产生影响,而不影响其振幅的假设.因此,在波动方程早至波走时层析中可定义基于走时信息的目标函数:

(1)

其中,Δτ为观测地震数据与模拟地震数据中早至波走时差;s,r分别表示炮和检波点.作为走时层析反演方法中关键步骤,走时差的求取可通过观测与模拟地震数据早至波的互相关运算获得,即

(2)

其中,x表示空间位置,τ为观测与模拟地震数据的时移量,pobs和pcal分别代表观测与模拟波场.当式(2)取得最大值时,即可得到两种地震记录的走时差:

Δτ=argmax(f(x,τ)),τ∈[-T,T],

(3)

其中,T为最大走时差.此时,式(2)中互相关函数f(x,τ)对时移差τ的一阶偏导数应为零,即:

(4)

由公式(1)可知,目标函数J的求取是一个最小二乘问题,可通过梯度法来进行迭代求解(本文中采用共轭梯度法).基于梯度法求解最小二乘问题的关键在于如何计算目标函数梯度,即目标函数J对速度的一阶导数:

(5)

其中,

(6)

将式(4)代入式(6),可得:

(7)

根据Tarantola(1987),式(7)中模拟波场对速度的一阶导数可表示为

(8)

(9)

将公式(9)代入公式(5)中,可得波动方程早至波层析反演的梯度与伴随源表达式:

(10)

其中,δ(x,t)为走时差伴随源.通过最优化方法,可实现波动方程初至波层析反演速度建模.

1.2 基于声波测井数据稀疏反演速度建模方法

平面波分解(PWD)是一种获取地下构造信息的有效手段,该方法类似于预测误差滤波器,即:

rix,ix+1=dix+1-Fix,ix+1(α)dix,

(11)

其中,rix,ix+1代表地震道dix+1与预测地震道dix之间的剩余误差.Fix,ix+1(α)表示与地下构造倾角α相关的预测算子.以早至波层析反演结果为偏移速度对其实施逆时偏移可获得初步成像剖面.通过对成像剖面进行平面波分解,可获得相应地下构造的倾角剖面α.由于地下倾角剖面是源自于成像剖面,因此,该倾角剖面较早至波层析反演速度结果能够更好地反映地下构造信息.

速度剖面是一种全网格数据体,而声波测井数据是一种稀疏数据.在不考虑速度频散效应的前提下,声波测井数据可以看作为一种速度剖面的抽样表达,即:

w=Mv,

(12)

其中,w和v分别表示声波测井数据和地下速度剖面,M表示采样矩阵.由于公式(12)中采样矩阵M具有不确定性,因而无法直接利用矩阵变换求取速度场v.可将公式(12)转化为最小二乘最优化问题进行求解:

(13)

式中速度场v为待求速度场.为了在速度场反演中引入平面波分解所得倾角信息,速度场可变换为

v=Sq,

(14)

其中,S是与构造角度相关的平滑算子,q为反演系数基(如傅里叶或小波分解系数).将公式(14)代入公式(13)可得:

(15)

式中ε为正则化因子.通过最优化方法(如共轭梯度)迭代求解公式(15),可得基于声波测井数据的速度场:

(16)

1.3 波动方程全波层析反演速度建模

波动方程层析反演速度建模方法往往仅利用早至波或反射波信息,而无法利用全波信息,究其原因主要在于:地震波类型越全,反演非线性越强;层析反演速度建模初始模型往往为均匀或光滑速度模型,无法提取模拟地震数据与观测地震数据之间的走时差.而本文提出的利用早至波走时层析结合声波测井数据稀疏反演进行初始速度建模方法能够较好地解决初始模型建立问题,因此可利用全波信息进行层析反演速度建模.

全波层析速度反演速度建模方法最重要的是走时差的计算.本文利用局部域互相关方法计算模拟数据与观测数据的全波时变走时差.利用高斯窗函数将模拟及观测地震记录由数据域转换为局部数据域:

(17)

其中,Dobs(x,k,t)与Dcal(x,k,t)分别代表局部数据域的观测与模拟地震数据,k为高斯窗函数的长度.对Dobs(x,k,t)与Dcal(x,k,t)进行互相关运算:

(18)

并搜索互相关函数C(x,τ,t)的最大值,可获得局部数据域观测地震数据Dobs(x,k,t)与模拟地震数据Dcal(x,k,t)的走时差:

(19)

值得说明的是,局部域互相关算法中高斯窗函数长度的选择影响着走时差的准确性,本文采用最优化方法反演获得高斯窗函数长度,即:

(20)

由式(20)可知,当∂Δτ(x,k,t)/∂k=0时,可获得最佳高斯窗长度.

由于全波走时层析反演方法的初始速度模型采用的是1.1—1.2节中反演的速度结果,其对应的模拟数据已包含各种类型地震波信息.因此,无需利用反偏移算法获得模拟地震记录,避免了传统Born反偏移所导致的计算量大问题.基于此,全波走时层析反演伴随源与梯度仅需将初至波走时层析速度层析反演中初至波走时差改为全波走时差,即:

(21)

1.4 反演步骤及策略

如图1所示,以声波测井数据及地震记录为出发点,以一维线性速度模型为初始模型,通过利用地震波全波走时及声波测井信息,提出了一种基于声波测井数据的全波走时层析反演速度建模方法,具体反演步骤及策略为:

图1 基于声波测井数据全波走时层析反演速度建模技术流程图

(1) 利用互相关方法,计算模拟数据与观测数据的早至波走时差,并利用公式(10)更新初始线性速度模型,获得初至波层析反演速度模型v1;

(2) 基于步骤(1)所得速度模型及地震记录,通过逆时偏移算法得到初始成像剖面,并采用平面波分解方法获得地下构造倾角剖面α;

(3) 利用地下构造倾角剖面及声波测井数据,通过公式(15)进行稀疏速度反演,获得精细速度模型v2.考虑到声波测井速度通常为非全井段数据,因此,在该步骤中需利用步骤(1)中建立的速度模型v1对其进行补全,即公式(12)中测井数据可表示为

(22)

式中void表示空缺.

(4) 以步骤(3)所得速度模型v2为初始速度模型,利用公式(18)—(20)反演获得高斯窗函数的长度及走时差伴随场,通过共轭梯度法更新速度模型,最终获得基于声波测井数据的波动方程全波走时层析反演速度模型.

2 数值模拟

在本节中,将选用较为复杂的Marmousi模型对本文方法的有效性进行验证.首先,利用早至波波动方程层析反演获得速度模型v1;其次,以速度v1为初始速度模型,基于声波测井速度稀疏反演更新获得速度模型v2;最后,以速度模型v2为初始模型,利用波动方程全波层析反演方法获得最终偏移速度模型v3.

如图2a所示为Marmousi速度模型,该模型网格大小为663×234,空间采样间隔为10 m.设置炮点均匀分布于地表,共40炮,采用全排列接收,最大偏移距为6.63 km.震源为主频15 Hz的雷克子波,其时间采样间隔为0.5 ms,共模拟7000个时间步长.为使本文方法更能符合野外实际情况,所有地震记录均滤除了5 Hz以下的频率信号.

图2 Marmousi速度模型

2.1 早至波波动方程层析反演速度建模

采用如图2b所示一维线性速度模型作为初始速度模型,其速度范围介于1.95～4.7 km·s-1之间.图3a为滤除5 Hz以下低频且去除直达波后单炮地震记录,炮点位于(3300.0 m,0.0 m).由于早至波波动方程层析速度反演方法仅利用地震记录中的早至波信息,因此,本文采用图3a地震记录中直达波以下200 ms内的地震数据作为早至波地震记录(如图3b).利用公式(2)互相关算法,可获得如图4a所示模拟早至波与观测早至波地震记录的走时差场,继而根据公式(10)中走时差伴随源算法,可得早至波波动方程层析反演速度建模走时差伴随源(图4b).

图3 Marmousi模型地震记录

图4 Marmousi模型

根据图4b所示走时差伴随源及公式(10),采用共轭梯度法,通过30次迭代更新初始速度模型(图2b),获得如图5所示早至波波动方程层析反演速度剖面.从反演结果可以看出,早至波波动方程层析反演速度建模方法能较好的刻画浅层及大尺度速度信息,而对中深层及小尺度速度信息的反演效果有限.而从图6所示逆时偏移成像剖面中也可进一步看出,浅层及小倾角地层能够进行较好的成像,而深层及大角度细节构造产生了偏移假象.图7所示为x=1500 m、x=3500 m及x=5500 m处角道集剖面,其中左侧与右侧分别为早至波层析反演速度及真实速度结果,从图中对比可清晰看出角道集同向轴无法拉平.由此可知,仅依靠早至波波动方程层析反演速度建模方法,无法满足目前对中深层及陡倾角等细节信息的成像需求.

图5 早至波波动方程层析反演速度模型(迭代30次)

图6 基于图5速度的成像剖面

图7 Marmousi模型角道集

2.2 基于声波测井数据稀疏反演速度建模方法

由2.1节可以看出,传统早至波层析反演速度建模方法仅能实现浅层或大尺度速度建模,而无法实现小尺度及中深层速度建模.事实上,目前仅有全波形反演方法(FWI)能够实现中深层及小尺度复杂构造速度建模.对比图5及图6可知,成像剖面相比于速度剖面更能体现地下的小尺度及细节信息.因此,基于图6所示Marmousi模型成像剖面,采用平面波分解方法提取包含地下更多细节信息的倾角信息.事实上,从成像剖面中提取角度信息类似于通过反偏移提取一次反射波地震记录的思路,均利用了成像剖面的高分辨率特性.然而,在平面波分解求取地下构造倾角信息过程中存在奇点影响.本文采用滑动窗口平滑方法提取地下倾角信息,以避免成像剖面存在的奇点影响.从图8所示倾角剖面可以看出,Marmousi模型的陡倾角及中深层倾角信息均可正确体现.

图8 图6提取的倾角信息

图9所示为假定的测井数据,共5口井,分别位于x=1300 m,x=2300 m,x=3300 m,x=4300 m,x=5300 m.由于实际测井数据往往仅包含目的层段,且井位置具有空间随机性.因此,为符合实际测井数据,本文所假定5口井数据均缺失一定地层段的测井数据.基于图8所示地下构造倾角信息,并结合公式(15)—(16)所示稀疏反演方法,可获得基于声波测井数据稀疏反演速度(图10).通过图10与波动方程早至波层析反演速度剖面(图5)对比可以看出,该方法能够较为准确获得地下大尺度速度信息,尤其是深层速度信息.然而,基于声波测井数据稀疏反演方法的实质是声波测井速度在滑动窗口内沿倾角方向进行横向插值,虽然可反演得到大尺度及深层速度信息,但其速度值可能偏离真实速度,且在插值过程中存在奇点.因此,利用声波测井速度稀疏反演速度剖面无法直接用于地震成像,需通过层析反演方法对速度剖面进行进一步修正.

图9 Marmousi模型声波测井数据(5口井)

图10 基于图9所示声波测井数据稀疏反演速度剖面

2.3 波动方程全波层析反演速度建模

图10所示速度模型已包含较为丰富的浅层及中深层速度信息,因此以图10所示速度模型作为初始速度模型,可获得地震波数较为丰富的模拟地震数据,而无需额外利用反偏移算法获得对中深层速度建模更为重要的反射波信息.由1.3节可知,波动方程全波层析反演速度建模方法的关键在于获取准确的全波走时差Δτ(x,t),而在全波走时差计算过程中选择合适的高斯窗长度是个难题.图11所示为仅含单个地震子波走时差随高斯窗大小变化曲线,从图中可以看出,该曲线随时窗长度的增加呈现先减小后增大的趋势;当走时曲线取最小值时,可得到最佳走时差.

图11 走时差随高斯窗大小变化曲线

如图12a所示,为验证公式(18)—(20)高斯窗大小选取方法的有效性,本文采用包含多子波的单道地震数据来进行验证,其中红色与蓝色曲线分别代表模拟与观测地震道.图12b—e分别代表采用高斯窗长度为240 ms、636 ms、700 ms及1000 ms时获得的走时差曲线,其中实线与虚线分别代表计算与真实走时差.从图中可以看出,选择较短或较长的高斯窗长度均无法得到准确走时差.而由公式(18)—(20)求取的走时差曲线图12c可知,其计算走时差(实线)与真实走时差(虚线)高度吻合,证实了本文提出的动态计算高斯窗的长度的可行性,继而在高斯窗内求取走时差的方法是有效可行的.

图12 (a) 单道模拟与观测地震记录; (b)—(e) 不同高斯窗大小走时差信息,其中(b) 240 ms,(c) 636 ms,(d) 700 ms,(e) 1000 ms

本节以图10所示速度模型为初始速度模型,因而模拟地震数据(图13a)中已包含丰富的地震波信息.利用公式(18)—(20)获得走时差(图13b)及走时差伴随源(图13c),采用共轭梯度法对该初始速度模型进行更新.图14展示了经过30次迭代后的层析反演结果,对比图5可知,本文方法的反演精度明显优于传统早至波层析反演方法,能够有效地获得中深层及小尺度地下速度信息.

图13 Marmousi 模型

图14 全波走时层析反演速度剖面

以图14所示速度模型为成像速度,经逆时偏移得到成像剖面(图15),从成像剖面可以看出,高陡构造与深层构造均能够实现高分辨率成像,且与真实模型高度吻合.为更直观体现本文方法有效性,相同地提取了x=1500 m、x=3500 m及x=5500 m处角道集剖面(图16),其中左侧与右侧分别为本文层析反演速度及真实速度结果.从图中可以看出,深层角道集同向轴能够较好的拉平,且与真实速度模型角道集吻合度较高.因此,本文提出的利用声波测井资料进行层析速度建模的方法是有效可行的.

图16 Marmousi模型角道集

在实际应用中,测井数据往往存在空间分布不均且数量较少问题.为进一步验证本文提出的基于声波测井数据的波动方程全波层析反演速度建模方法的有效性,如图17所示,减少Marmousi模型中测井数量为3口井,其井位置分别位于x=1000 m,x=4000 m,x=5500 m,且均缺失一定地层段的测井数据.相同地,基于图8所示地下构造倾角信息,利用本文中稀疏反演方法,可获得如图18所示速度剖面.通过对比图18与图10可以看出,当测井数量减少且分布不均时,仍可反演得到较为准确的地下大尺度及深层速度信息,充分证明了本文方法的稳定性和有效性.

图17 Marmousi模型声波测井数据(3口井)

图18 基于图17所示声波测井数据稀疏反演速度

以图18所示速度模型为初始速度模型,利用本文中全波层析反演速度建模方法,经过30次迭代(与包含5口井时相同迭代次数),得到如图19所示速度反演结果.从图19中可清晰看出,中深层及小尺度地下速度信息得到了较好的恢复,且反演结果与图14所示反演剖面吻合度较高.因此,当模型中测井数量较少且分布不均时,本文方法依然具有较高的反演精度.

图19 全波走时层析反演速度剖面

相同地,以图19所示速度模型为成像速度,经逆时偏移得到成像剖面(图20),从成像剖面可以看出,高陡构造与深层构造均能够实现高分辨率成像,且从图21所示角道集剖面(x=1500 m、x=3500 m及x=5500 m,其中左侧与右侧分别为仅含3口井时层析反演速度及真实速度结果)可以看出,利用本文反演方法获得的角道集同向轴能够较好的拉平,且与真实速度模型角道集吻合度较高,可以充分说明,当模型中仅包含3口声波测井时,本文速度建模精度能够满足逆时偏移成像对速度精度的要求.综上所述,本文提出的利用声波测井资料进行层析速度建模的方法对测井数据的数量依赖性较小,具有较强的稳定性,能够有效的提高中浅层及深层的速度建模精度.

图20 基于图19的成像剖面

3 讨论

本文在走时仅受速度控制的假设下实现了中浅层及深层高分辨率速度建模,而未考虑其他介质参数对走时的影响.但当地下构造较为复杂时,走时将受其他多种地下介质参数影响,例如各向异性及衰减(Q)参数.

由于各向异性参数和衰减参数与速度均存在耦合关系,因而理论上需对速度、各向异性及衰减参数进行联合层析,但由于各向异性参数和衰减参数对走时不敏感,导致联合层析效果往往较差,且存在多参数反问题的串扰问题.因此,在实际层析反演中,往往采用先对速度进行层析建模,然后对其他参数进行层析建模的两步法策略.由于声波测井数据能提供准确的速度信息,因此,在多参数联合层析中,可利用本文中提出的稀疏反演策略建立较为准确的初始速度模型,有望降低速度与其他参数的耦合性,可在一定程度上提高多参数反演的精度.

4 结论

本文针对中深层走时层析速度建模方法面临的难题,提出了一种基于声波测井资料的三步法波动方程全波走时层析反演速度建模方法.以一维线性速度模型为初始速度模型,采用早至波波动方程层析反演实现对地下大尺度及浅层构造低波数速度建模,继而利用声波测井速度及初始成像剖面结合稀疏反演获得深层及较高分辨率的中波数速度场;以稀疏反演速度场为初始模型,利用波动方程全波层析反演方法获得最终偏移速度场剖面.

Marmousi模型算例表明,本文所提出的速度建模方法不仅能够实现对浅层及中深层的较高精度速度建模,而且计算效率较高,克服了传统方法反演精度低或计算代价高的缺陷.但值得注意的是,本文所提出的层析反演方法及策略需利用声波测井数据,如何有效避免声波测井数据的空间分布及测井数据段对最终层析反演结果产生的影响有待进一步研究.