地球物理学报  2020, Vol. 63 Issue (2): 666-686   PDF    
基于走时约束的分频分步零偏VSP波动方程多参数反演
韩璇颖1,2, 印兴耀2, 曹丹平2     
1. 中国石油化工股份有限公司石油物探技术研究院, 南京 211103;
2. 中国石油大学(华东)地球科学与技术学院, 青岛 266580
摘要:在地震勘探中,描述复杂介质的正演和反演问题通常包含许多反映介质不同特性的参数.同时获得这些参数对进行更准确的岩性描述和油藏预测具有重要的理论和现实意义.为了提高频率域黏弹性波动方程的零偏VSP多参数反演的精度,本文对多参数反演的可行性进行分析,明确了目标函数的敏感程度及参数之间的耦合情况,提出了一种基于走时约束的分频分步多参数反演策略.首先利用零偏VSP资料构建先验信息,然后分别利用高、低频数据进行两步反演,也就是"三个参数反演+五个参数反演"的过程,以提高反演的稳健性和精度.利用此方法可同时得到零偏VSP数据可靠的弹性波速度、密度和品质因子,为精确的时-深关系及含油气的解释和预测奠定基础,同时也可以为地面地震叠前反演提供可靠有效的约束,增强地面地震反演精度.
关键词: 频率域      黏弹性波动方程      五参数反演      零偏VSP     
Multi-parameter inversion of frequency-divided zero-offset VSP wave equation based on traveltime constraint
HAN XuanYing1,2, YIN XingYao2, CAO DanPing2     
1. Sinopec Geophysical Research Institute, Nanjing 211103, China;
2. China University of Petroleum(East China), Qingdao 266580, China
Abstract: In seismic exploration, forward modeling and inverse problems of describing complex media usually contains many parameters which reflect the different properties. In order to improve the accuracy of zero offset VSP multi-parameter inversion in the frequency domain, the feasibility of such inversion is analyzed. The sensitivity of the objective function and the coupling between the parameters are clarified, and a strategy for frequency-divided wave equation multi-parameter inversion based on traveltime constraint is proposed. First, the prior information is constructed by using the zero offset VSP data, then a two-step inversion is carried out on low-and high-frequency data, which is a process of "three parameter inversion + five parameters inversion", in order to improve the robustness and accuracy of the inversion. By using this method, the reliable elastic wave velocities, density and quality factors of the zero offset VSP data can be obtained simultaneously. It lays a foundation for accurate time depth relationship and interpretation and prediction of oil and gas. It can also provide a reliable constraint for the ground seismic pre-stack inversion and enhance accuracy of ground seismic inversion.
Keywords: Frequency domain    Viscoelastic wave equation    Five parameter inversion    Zero offset VSP    
0 引言

中国东部老油田石油勘探的主要地质目标已经转为隐蔽油气藏和复杂油气藏,对勘探精度和解释结果提出了更高的要求.速度和密度是介质基本的物性参数,同时,实际地层介质为黏弹性介质,描述黏弹性介质固有衰减的品质因子Q与介质内部的渗透率、饱和度等参数有关(武银婷和朱光明,2017).因此,在地震勘探中,对介质的纵横波速度、密度和Q值等多种参数同时反演的研究具有重要的理论意义和实际应用价值(王小杰等,2011印兴耀等,2015宗兆云等, 2012a, b).

作为一种井中地震勘探方法,VSP在油气田勘探开发研究中,发挥着连接地面地震、钻井、地质资料的桥梁作用(未晛等,2019).Li和Li(2018)在巴黎盆地利用VSP数据和地面地震数据进行裂缝性储层描述,结果表明了剪切波各向异性解释在揭示裂缝性储层的横向变化方面具有优势.利用VSP方法获取的准确的纵横波速度、密度和Q值,也可以为之后地震地质层位的标定,成果解释及油气预测奠定基础,同时为叠前反演提供可靠的约束.人们通常先对VSP资料进行波场分离(Rao and Wang, 2016),然后利用走时的信息获取参数.Stewart(1984)Lines等(1984)对走时反演VSP速度做了系统的研究;邹强等(2003)提出同时利用阻抗和走时反演层速度;姚忠瑞等(2007)以北101井VSP测井为例,反演零偏VSP资料的横波速度.但在这个过程中往往会受到走时拾取误差、地震波能量强度等多种因素的影响,得到的速度往往不够准确(张继国和牟风明,2006).同时,仅仅利用走时信息不能直接反演出与时间无关的密度的信息,也不利于品质因子Q的反演.因此,要想同时获得精确的纵横波速度、密度和Q值,不仅仅需要走时,还需要波形等其他波场信息.波动方程反演可以充分地利用走时及波场的其他信息如振幅、频率和相位等,具有高分辨率的优势,可以提高反演精度(杨勤勇等,2014张广智等,2017廖建平等,2018李冀蜀等,2019梁展源等,2019).Owusu等(2016)根据阿拉伯湾VSP数据信息,利用全波形反演恢复了合理的速度模型;Egorov等(2017)利用全波形反演可以得到VSP纵波的速度同时量化其变化情况.

前人的研究可知,利用波动方程进行纵波速度、横波速度和密度等多参数的同时反演,速度的反演结果较好,但密度的结果会出现偏差(Liu et al., 1995).为了提高地面地震多参数反演的精度,国内外学者提出了多种反演策略.Jeong等(2012)为了提高多参数反演中密度的反演效果,提出可以先将密度设为常数反演出弹性模量,然后换算成速度再进行速度和密度的同时反演.Prieux等(2013)为了提高多参数反演的稳定性,利用了大偏移距和全偏移距的数据进行反演.杨积忠等(2014)提出了一种有利于速度、密度同时反演的两步反演策略,可以同时得到比较精确的速度、密度.Yang等(2016)研究了基于散射积分的截断高斯-牛顿法求解高斯-牛顿方程,可以同时反演密度和纵波速度.Sun等(2017)将该方法进一步改进,应用到弹性介质中.关于品质因子的多参数反演,Malinowski等(2011)在频率域同时反演了速度和品质因子;Operto等(2013)总结了多参数同时反演的基本原则及如何选择合适的参数组合等问题,并给出了数据驱动和模型驱动的黏弹性波动方程成像策略.

但前人波动方程多参数反演多是基于地面地震的情况(Brossier et al., 2009Brossier,2011Guitton et al., 2012),并没有对VSP资料尤其是含有透射波信息的情况进行研究,同时针对纵横波速度、密度和品质因子QPQS五个参数的同时反演分析较少.针对这几个方面,本论文首先对零偏VSP资料的波动方程的多参数反演进行可行性分析,然后提出一种基于走时约束的分频分步多参数反演策略.首先利用零偏VSP数据的地震资料得到纵、横波速度和密度的先验信息,然后再进行分频分步反演.在分频分步反演过程中先利用地震记录的低频信息再利用高频信息依次反演,一方面可以提高反演的稳健性,另一方面在低频反演过程中可以进行合理的假设,减少反演参数,为高频精确反演提供约束,提高了反演的稳健性和精度.

1 方法基础 1.1 波动方程正演

本文反演的基础是频率域二维波动方程有限差分正演.当体力为零时,弹性波波动方程为

(1)

这里ρ为密度,λμ为拉梅系数.本文使用优化9点有限差分算法(Štekl and Pratt, 1998)进行正演模拟.

地震波的衰减作为黏弹性介质的一种重要的属性,逐渐成为人们研究的热点.介质本身所固有的衰减特性通常用品质因子Q值来描述.在频率域通过使用与频率和Q值相关的复速度,将吸收衰减项引入到地震波传播方程中,可以模拟出黏弹性介质中地震波能量衰减的过程,进行黏弹性介质的数值试验.地震波复波数k(ω)的公式为

(2)

其中,α为衰减系数,v(ω)为相速度.根据波数和速度的关系可以得到复速度c(ω)为

(3)

又因为:

(4)

ω>0时,复速度的表达形式又可以表示为

(5)

再结合Futterman Q值衰减模型的公式,推导出λμ的表达式,进行频率域黏弹性波动方程正演模拟.

1.2 波动方程反演

在地震全波形反演中,首先需要确定反演的目标函数,用来衡量观测数据和模拟数据之间的匹配程度.使用L-2范数作为目标函数是最简单和最常用的方法,具有计算效率高的特点.其形式为

(6)

其中,d为观测的地震数据,f为模拟的结果,F为目标函数. fF与纵横波速度、密度和品质因子有关.考虑到本文应用的零偏VSP资料数据量小,同时为了避免局部优化算法的不足(易陷入局部极值等),这里利用快速模拟退火算法进行反演.

2 反演可行性分析 2.1 纵、横波速度和密度三个参数的反演

本文首先研究三个参数的反演.分析目标函数对不同参数的敏感程度,研究不同参数之间的耦合程度,是建立多参数反演策略的有效途径.

2.1.1 三参数的敏感性分析

杨积忠等(2014)已对地面地震数据目标函数的敏感度进行分析,提出无论是大偏移距,中偏移距还是小偏移距,L-2范数的目标函数对密度较小尺度摄动较敏感,但利用反射波无法反演较大尺度的密度摄动成份.本文则对含有透射波的零偏VSP数据进行目标函数的敏感度分析.从基于零偏VSP资料波动方程出发,详细研究数据随纵、横波速度和密度的不同空间摄动尺度的变化关系,并系统分析这些变化关系在目标函数下的敏感程度,为反演提供理论指导.试验首先选择一个参考模型,并对参考模型作不同尺度的正弦摄动,然后基于参考模型和摄动后模型分别进行正演模拟.利用模拟的地震数据分别绘制出目标函数随参数摄动尺度变化的曲线.

(1) 参考模型

参考模型选为MarmousiⅡ模型的纵、横波速度和密度的一部分,模型深度为200 m,纵、横波速度和密度的模型随深度的变化如图 1所示.

图 1 模型参数 (a)纵波速度; (b)横波速度; (c)密度. Fig. 1 Model parameters (a) P-wave velocity; (b) S-wave velocity; (c) Density.

(2) 模型的摄动方案

对各个参数作不同尺度的正弦摄动.参数沿深度的摄动方式为(杨积忠等,2014):

(7)

这里,a(z)代表不同的参数(纵、横波速度和密度),a0(z)表示参数在深度为z时的值;Δa表示摄动量;λ为摄动波长,也就是摄动的尺度;amax表示参数的最大值.摄动幅度为最大值的10%.为了比较完整地再现整个摄动尺度范围内目标函数随尺度的变化,本文在小尺度采样较密,在大尺度采样较稀,共选取了196个尺度样点.

(3) 敏感性分析

根据上述试验方案,利用图 1的模型得到利用零偏VSP数据的地震记录计算的L-2范数目标函数随不同参数的不同摄动尺度的敏感曲线,如图 2所示.目标函数随VP摄动的曲线数值相对较大,说明VP是影响目标函数变化的重要参数.目标函数随VS的摄动曲线在较小尺度处呈现出峰值,而在较大尺度处变化平缓,值较小,说明目标函数对VS的较小尺度的摄动具有一定的敏感度,而对长波长成分的敏感度较低.因此仅仅利用零偏VSP资料数据构造的常规目标函数,在没有有效约束条件的情况下对VS的较小波长成分具有较好的反演能力,但难以恢复中长波长的成分.在较小的摄动尺度下,密度的扰动对目标函数的敏感曲线有两个峰值,说明目标函数对密度的较小尺度的摄动具有较好的敏感度.而对大中尺度的密度摄动,虽然较低,但并不为零,这也为利用零偏VSP资料反演密度提供了可能.当同时扰动纵波速度和横波速度时,目标函数的变化曲线同仅扰动纵波速度时的几乎相同,同样,同时扰动密度和纵波速度及同时扰动密度和横波速度的目标函数变化曲线,同仅扰动敏感度较高的参数的变化曲线很接近.敏感度较低的参数的变化,往往不易在目标函数的变化曲线中体现出来.这也给多参数反演提供了难度,因此需要有效的约束信息再进行反演.

图 2 L-2范数的目标函数随不同参数不同摄动尺度的敏感曲线 Fig. 2 Sensitive curves of L-2 norm objective function to different perturbation scales with different parameters
2.1.2 三参数耦合性的分析

为了进一步分析在波动方程反演过程中纵、横波速度和密度之间的相互耦合关系,对突变模型进行如下数值实验:突变三个参数中的某个参数,再同时反演三个参数,观察该参数突变对其他参数的影响.本节数值试验采用的模型为500 m×100 m的层状模型,网格离散间距为10 m,震源采用主频20 Hz的Ricker子波.反演过程中纵、横波速度和密度的初值都为背景值(2000 m·s-1, 1500 m·s-1, 2000 kg·m-3). 图 3显示的为横波速度有突变,纵波速度和密度都为背景值模型的三参数反演结果.很显然,在反演过程中,由于横波速度的异常同时影响了三个参数的反演结果.图 4为只对纵波速度模型进行突变时的三参数反演结果,结果同图 3类似,但在速度突变量相同的情况下,纵波速度突变比横波速度突变对密度影响更大.图 5为突变密度时的反演结果,密度的突变,也同时影响了纵横波速度的反演结果,但影响较小.综合图 345,可以知道速度对密度的反演的影响较大,密度对速度的反演的影响相对较小,这同地面地震反演的结论一致.因此,在进行三参数反演之前,如果能够先得到速度较小的变化范围区间,然后再进行波动方程反演,此时,速度的变化较小,对密度的影响也较小,三个参数之间的耦合效应较不明显,便可以得到较为理想的结果.

图 3 横波速度突变情况下纵波速度(a)、横波速度(b)和密度(c)的反演结果 虚线表示真实值,实线表示反演值. Fig. 3 Inversion results of P-wave velocity (a) and S-wave velocity (b) and density (c) when S-wave velocity has a sudden change Dot lines represent real values and solid line represent inversion values.
图 4 纵波速度突变情况下纵波速度(a)、横波速度(b)和密度(c)的反演结果 虚线表示真实值,实线表示反演值. Fig. 4 Inversion results of P-wave (a) and S-wave velocity (b) and density (c) when P-wave velocity has a sudden change Dot lines represent real values and solid lines represent inversion values.
图 5 密度突变情况下纵波速度(a)、横波速度(b)和密度(c)的反演结果 虚线表示真实值,实线表示反演值. Fig. 5 Inversion results of P-wave (a) and S-wave velocity (b) and density (c) when density has a sudden change Dot lines represent real values and solid lines represent inversion values.

从上面对目标函数的敏感性分析及各个参数的耦合性分析可以看出,零偏VSP资料L-2范数目标函数受纵波速度的影响较大,受横波速度和密度的中长尺度的扰动影响较小,且参数之间相互耦合相互影响.如果直接利用零偏VSP数据对纵横波速度和密度同时反演,得到的结果不够精确.因此,在三个参数反演的时候,可以先利用走时获得有效的先验信息再进行波动方程反演,目的是在进行弹性波波动方程反演之前有效的约束纵横波速度的取值范围,以减小在反演中纵横波速度和密度之间的相互影响,从而得到更加精确的结果.具体过程为:先根据零偏VSP数据的垂直分量和水平分量走时信息,分别得到纵、横波速度;再将走时得到的速度作为波动方程反演的初值,并进行上下浮动,作为反演中速度的上下边界值;然后进行弹性波三参数波动方程反演,得到更加精确的纵、横波速度和密度.

2.2 Q值的敏感性分析

在进行频率域含衰减的地震正演模拟的过程中,引入了两个参:QPQS,分别表示对纵波和横波的衰减.这里采用同上文一致的试验过程分析QPQS对高低频率成分地震记录的敏感程度.QPQS的参考模型深度为200 m,模型随深度的变化如图 6所示.图 7是利用零偏VSP数据高、低频成分计算的目标函数随QPQS的不同摄动尺度的敏感曲线.横坐标为扰动尺度λ,纵坐标为目标函数,也就是利用摄动后的模型进行正演模拟的地震记录和利用参考模型正演模拟得到的地震记录的差值.实线表示只利用地震记录高频成分(16~30 Hz)的目标函数随QPQS不同摄动尺度的敏感曲线,虚线则表示只利用地震记录低频部分(1~15 Hz)的结果.两个图中实线都比虚线的值大很多,因此,无论是QP还是QS,同低频部分相比,地震记录高频成分的目标函数值很大,Q值对地震记录的高频成分影响较大,而对地震记录的低频成分影响较小.由于QPQS对地震记录中不同频率段影响不同,可以考虑在反演中进行分频反演,以提高反演的准确性.

图 6 QP(a)和QS(b)的参考模型 Fig. 6 Reference models of QP (a) and QS (b)
图 7 低频和高频成分地震记录的目标函数随QP(a)和QS(b)的不同摄动尺度的敏感曲线 虚线表示低频成分,实线表示高频成分. Fig. 7 Sensitive curves of the objective function of low frequency and high frequency components in seismic records with different perturbations of QP (a) and QS (b) Dot lines represent low frequency component and solid lines represent high frequency component.
2.3 先验信息的构建

零偏VSP地震记录可以分别利用垂直分量和水平分量的走时信息得到纵、横波速度,作为反演的一个重要的先验信息.但由于真实的地震记录含有衰减,根据时深关系得到的并不是真实的弹性波速度.因此,需要进一步分析含衰减地震记录中走时信息得到的纵、横波速度和真实的弹性波速度的差别.

本文选取三种模型进行正演模拟,分别为弹性模型(VP=2000 m·s-1)、QP=50的衰减模型和QP=10的衰减模型,这三种模型的横波速度和密度均相同(VS=1000 m·s-1ρ=1000 g·cm-3).衰减模型中当QP趋于无穷大时纵波速度趋于2000 m·s-1.对正演得到的地震记录分别进行初至拾取,通过时深关系得到纵波速度并对比分析.图 8为弹性模型的结果,得到的平均速度为1998 m·s-1图 9图 10分别为QP=50和QP=10的衰减模型的结果,得到的平均速度分别为1976.1 m·s-1和1927.7 m·s-1.可以看出QP可以引起地震记录能量的衰减,纵波速度的减小.但速度的差异并不大,即使衰减很大(QP=10),在含衰减的地震记录中根据走时信息得到的衰减后速度值,不小于真实速度的5%.综上所述,在含衰减的地震记录中根据走时信息得到的速度依然可以作为反演弹性速度的先验信息.

图 8 VP=2000 m·s-1弹性模型的地震记录(a),相应的时距曲线(b)及得到的纵波速度(c) 虚线表示参考值(2000 m·s-1),实线表示根据走时信息得到的速度. Fig. 8 The seismic records of elastic model (VP=2000 m·s-1) (a), corresponding time-depth curve (b) and obtained P-wave velocity (c) Dot line represents reference value (2000 m·s-1) and solid line represents velocity based on traveltime information.
图 9 (a) QP=50的衰减模型的地震记录(a),相应的时距曲线(b)及得到的纵波速度(c) 虚线表示参考值(2000 m·s-1),实线表示根据走时信息得到的速度. Fig. 9 The seismic record of the attenuation model (QP=50) (a), corresponding time-depth curve (b) and the obtained P-wave velocity (c) Dot line represent reference value (2000 m·s-1) and solid line represent velocity based on travel time information.
图 10 (a) QP=10的衰减模型的地震记录(a),相应的时距曲线(b)及得到的纵波速度(c) 虚线表示参考值(2000 m·s-1),实线表示根据走时信息得到的速度. Fig. 10 Seismic record of attenuation model (QP=10) (a), corresponding time-depth curve (b) and obtained P-wave velocity (c) Dot line represents reference value (2000 m·s-1) and solid line represents velocity based on traveltime information.
3 基于走时约束的分频分步VSP波动方程多参数反演策略

为了避免局部优化算法的劣势,同时考虑到零偏VSP数据量小,本文采用模拟退火算法进行波动方程纵横波速度、密度和QPQS五个参数的反演.直接使用含衰减的地震记录反演这五个参数,反演结果的准确性和稳定性都不够,因此,本文提出了一种基于走时约束的分频分步反演策略来提高反演的可靠性.反演过程如图 11所示.

图 11 基于走时约束的分频分步VSP波动方程多参数反演流程图 Fig. 11 Flowchart of frequency-divided zero offset VSP wave equation multi-parameter inversion based on travel time constraint

整个反演分为两个过程:先验信息的构建和分频分步反演.先验信息的构建是利用地震记录中的走时信息,根据时深关系得到纵、横波速度的估计值并利用经验公式估计密度值,为接下来的分频分步反演提供速度和密度的初值及速度小范围的约束信息.由于Q值对地震记录的高频成分影响较大,而对地震记录的低频成分影响较小,因此,第二阶段的反演可以分为两步,第一步将QPQS值设为固定的常数,仅利用与Q值关系较小的地震记录低频成分(1~15 Hz,间隔为2 Hz)同时进行反演,这个过程相当于三参数反演,结合上一阶段的先验信息,可以得到较为可靠的纵、横波速度和密度的估计值.将该结果作为下一步反演的约束信息,第二步仅利用地震记录的高频成分(16~30 Hz,间隔为2 Hz)同时进行五参数反演,可以得到精确的QPQS、纵横波速度和密度,Q因子反演的稳定性也可以得到提高.

4 模型测试

正演模型如图 12所示,深度为900 m,偏移距为100 m,空间采样10 m,时间采样间隔1 ms,时间采样点数800,频率采样间隔1 Hz,最大频率30 Hz,震源为主频为30 Hz的雷克子波,采用PML吸收边界条件.正演模拟得到的含衰减的地震记录如图 13所示.

图 12 正演模型 (a)纵波速度; (b)横波速度; (c)密度;(d) QP;(e) QS. Fig. 12 Forward model (a) P-wave velocity; (b) S-wave velocity; (c) Density; (d) QP; (e) QS.
图 13 黏弹性模型正演的地震记录 (a)垂直分量;(b)水平分量. Fig. 13 Seismic records of forward modeling on viscoelastic model (a) Vertical component; (b) Horizontal component.

直接利用含衰减的地震记录反演得到的结果如图 14所示,图 15为能量随迭代次数变化的曲线.可以看出,虽然最终目标函数值已经很小,但反演结果较差.这是多参数反演的强烈非线性造成的,下面将采用如上所述的走时约束分频分步反演策略.

图 14 直接反演的结果 (a)纵波速度;(b)横波速度;(c)密度;(d) QP;(e) QS.虚线表示真实值,实线表示反演值. Fig. 14 Results of direct inversion (a) P-wave velocity; (b) S-wave velocity; (c) QP; (d) QS. Dot lines represent real values and solid lines represent inversion values.
图 15 能量随迭代次数的变化曲线 Fig. 15 Curve of energy varying with iteration number
4.1 先验信息的建立

由于从含衰减的地震记录中获得的速度值,一般不会小于真实速度的5%.因此,从零偏VSP初至走时获得的速度值可以作为下一步反演的初始模型,并且向上向下浮动5%可以作为反演的速度约束.

图 16为VSP地震记录的垂直分量及水平分量.图 17a为垂直分量利用走时得到的纵波速度和真实纵波速度的对比.可以看出利用走时得到的纵波速度还不够精确.将初至拾取的纵波速度上下浮动5%作为反演时的约束信息,如图 18a所示.地震记录的水平分量如图 16b所示,选取两条能量较强的下行波A和B,B波无疑是直达纵波.分析这两条波传播到同一深度的时间,粗略估计出它们的平均速度,发现A波的速度是B波速度的0.6左右,初步判断A波是转换横波.结合图 16a,A波在水平分量上能量较强,在垂直分量上则能量很小,说明是振动极化方向与波的传播方向相互垂直的波,也就是横波.因此,可以拾取A波(图 16b红点线)来得到横波速度.160 m之前的A波不明显,无法直接得到相应位置的横波速度,本文通过延伸得到了这部分浅层的速度.图 17b为利用走时信息得到整个模型的横波速度和真实横波速度的对比.由于利用转换横波获得的横波速度较利用直达波获得的纵波速度相比精度低,本文在将该横波速度作为下一步约束信息的时候,上下浮动了10%,如图 18b所示.对于密度的估计值,本文根据Gardner公式(马中高和解吉高,2005)与先验的纵波速度和横波速度求得,如图 18c所示.利用经验公式得到的密度有一定的误差,因此本文只将此结果作为下一步反演的初值.

图 16 (a) 真实记录垂直分量;(b)真实记录水平分量的下行波(160 m及以下) Fig. 16 (a) Real vertical component record; (b) Down-going wave of real horizontal component record(below 160 m)
图 17 速度曲线的对比 (a)纵波速度;(b)横波速度.虚线表示真实值,实线表示根据走时信息得到的速度. Fig. 17 Comparison of velocity curves (a) P-wave velocity; (b) S-wave velocity. Dot lines represent real values and solid lines represent velocity based on travel time information.
图 18 约束信息 (a)纵波速度;(b)横波速度;(c)密度.虚线表示真实值,实线表示反演的初始值,点划线表示约束. Fig. 18 Constraint information (a) P-wave velocity; (b) S-wave velocity; (c) Density. Dot lines represent real values, solid line represent initial value for inversion and chain lines represent constraints.
4.2 分频分步反演

QPQS对地震记录中的高频成分更加敏感.因此,本文采用分频分步反演策略:第一步只利用低频地震记录进行三参数反演,将Q值设为某个常数,作为已知量;第二步反演仅利用地震记录的高频部分,目的是为了得到精确的纵横波速度、密度和QPQS.在第二步反演过程中,将第一步反演得到的纵、横波速度和密度作为初始值,并上下浮动作为约束信息;Q值的初始值为第一次反演中的固定值.

本文在第一步反演过程中QPQS分别设为60和30.纵、横波速度和密度的初始模型和约束信息如图 18所示.第一步模拟退火算法中初始温度为200 ℃,Markov链长为180,迭代次数为628次.图 19为第一步反演之后得到的三个参数,能量随迭代次数的变化如图 20所示.第二步反演过程中速度和密度的初始值为第一步反演的结果,同时将该结果上下浮动5%作为这次反演过程中的最大值和最小值.QPQS的初始值同样为60和30,反演过程中QP约束在的(100,10)范围内,QS约束在(80,10)范围内.第二步反演中模拟退火算法的初始温度为200 ℃,Markov链长为180,迭代次数为628次.图 21给出了最后反演得到的五个参数.图 22显示的是深度为200 m处真实的地震记录和反演结果合成地震记录的对比.能量随迭代次数变化的曲线如图 23所示.最后进行了加噪测试,地震记录信噪比为4时得到的反演结果如图 24所示.

图 19 第一步反演的结果 (a)纵波速度;(b)横波速度;(c)密度.虚线表示真实值,实线表示反演值. Fig. 19 Results of the first step inversion (a) P-wave velocity; (b) S-wave velocity; (c) Density. Dot lines represent real values and solid lines represent inversion values.
图 20 第一步反演中能量随迭代次数的变化 Fig. 20 Energy varying with iteration number in the first step of the inversion
图 21 第二步反演的结果 (a)纵波速度;(b)横波速度;(c)密度;(d) Q;(e) QS.虚线表示真实值,实线表示反演值. Fig. 21 Results of the second step inversion (a) P-wave velocity; (b) S-wave velocity; (c) Density; (d) QP; (e) QS. Dot lines represent real values and solid lines represent inversion values.
图 22 地震记录的对比(h=200 m) 虚线表示真实记录,实线表示合成记录. Fig. 22 Comparison of seismic records (h=200 m) Dot line represents real record and solid line represents synthetic record.
图 23 第二步反演中能量随迭代次数的变化 Fig. 23 Energy varying with iteration number in the second step of the inversion
图 24 SNR=4时基于走时约束的分频分步反演结果 (a)纵波速度;(b)横波速度;(c)密度;(d) Q;(e) QS.虚线表示真实值,实线表示反演值. Fig. 24 Results of frequency-divided inversion based on travel time constraints (a) P-wave velocity; (b) S-wave velocity; (c) Density; (d) QP; (e) QS. Dot lines represent real values and solid lines represent inversion values.
5 实际零偏VSP资料的应用

在反演之前需要对某工区实际资料进行预处理,包括旋转定向和几何扩散修正.在实际三维三分量VSP数据采集过程中,井下三分量检波器随着深度的变化,沿着井径的垂直分量的方向不发生改变,但两个水平分量是随机取向的.经旋转定向可将三个分量旋转成正确的垂直、径向和切向分量.由于本文进行正演模拟的过程中采用的是二维模型,而实际资料为三维,还需要将实际数据转化为近二维数据.道集需要乘以T1/2将球面几何扩散转化为圆柱形几何扩散模式(Pica et al., 1990),修正之后再利用P波和SV波所在平面内的垂向和径向分量进行反演.

该零偏纵波偏移距为24 m,高程为1.83 m.本论文研究的深度范围为100~1000 m.零偏VSP资料的垂直和径向分量如图 25所示.

图 25 实际零偏VSP资料垂直分量(a)和径向分量(b) Fig. 25 Vertical component (a) and radial component (b) of real zero offset VSP data
5.1 先验信息的建立

按照第4节的内容,分别利用实际零偏VSP数据的垂直分量和径向分量进行初至拾取(图 26),得到的纵、横波速度如图 27ab所示.根据得到的纵、横波速度,依据经验公式求得的密度如图 27c所示.

图 26 (a) 实际零偏VSP资料垂直分量及(b)径向分量(下行波)初至拾取结果 Fig. 26 First arrival pick up results of vertical component (a) and radial component of real zero offset VSP data (down-going wave)
图 27 走时信息得到(a)纵波速度、(b)横波速度及经验公式得到的(c)密度 Fig. 27 P-wave velocity and (a) S-wave velocity (b) obtained from travel time information and density (c) obtained by empirical formula
5.2 五个参数的波动方程反演

本文研究的反演基础为频率域波动方程正演模拟.模型深度为900 m,偏移距为24 m.边界条件为PML.其他参数:空间采样10 m,时间采样间隔1 ms,时间采样点数2800,频率采样间隔1 Hz,最大频率30 Hz.震源采用爆炸震源,震源子波先通过叠加垂直分量下行波场的多个道集估计,然后再乘以一个波场与合成记录的能量比得到(Egorov et al., 2017).

整个反演过程采用前文基于走时约束的“三参数+五参数”反演策略.首先将利用走时得到的数据作为初始值进行第一步反演.得到的纵、横波速度和密度的结果如图 28所示.将该结果作为先验条件进行第二步反演.第二步纵、横波速度和密度的反演结果与相应的测井速度对比如图 29所示,由于该资料没有横波速度测井,这里仅对比了声波测井和密度测井.图 30Q值的反演结果,(a)为QP与该资料利用谱比法得到的品质因子的对比.图 31为深度600 m处利用反演结果的合成记录与实际资料的对比.以上反演结果的对比表明,纵波速度的反演结果与测井曲线拟合较好,密度和品质因子反演出现部分偏差.原因可能有以下几点:(1)采用的3-D转为2-D振幅校正不够精确,导致对密度及品质因子的反演造成了影响;(2)仅使用一个源进行反演,也可导致反演结果的不准确性;(3)层状模型假设造成了误差等.下一步可以从这几个方面进行完善.

图 28 零偏资料第一步反演的结果 (a)纵波速度;(b)横波速度;(c)密度.虚线表示初始值,实线表示反演值. Fig. 28 Results of first step inversion of the zero offset data (a) P-wave velocity; (b) S-wave velocity; (c) Density. Dot lines represent initial values and solid lines represent inversion values.
图 29 零偏资料第二步反演的结果 (a)纵波速度;(b)横波速度;(c)密度.虚线表示初始值,实线表示反演值,红线表示测井值. Fig. 29 Results of second step inversion of the zero offset data (a) P-wave velocity; (b) S-wave velocity; (c) Density. Dot lines represent initial values, solid lines represent inversion values and red lines represent logging values.
图 30 零偏资料第二步反演结果 (a) QP;(b) QS.虚线表示初始值,实线表示反演值,红线表示谱比法值. Fig. 30 Results of second step inversion of the zero offset data (a) QP; (b) QS. Dot lines represent initial values, solid lines represent inversion values, and red line represents the result of spectral ratio method.
图 31 深度为600 m零偏资料垂直分量地震记录的对比 虚线表示真实值,实线表示反演值. Fig. 31 Comparison of vertical component in zero-offset seismic records at depth 600m Dot line represents real values and solid line represents inversion values.
6 结论

本文详细研究了零偏VSP资料多参数波动方程反演.首先对反演目标函数的敏感性及参数的耦合情况等进行了可行性分析,为多参数反演提供了理论依据;然后提出一种基于走时约束的分频分步VSP波动方程五参数反演策略.本文充分利用地震记录中的走时的信息,得到纵、横波速度和密度的约束条件,作为反演的先验信息.接下来的分频分步反演一方面先用低频得到大尺度的构造信息,再进行高频精细反演,提高了多参数反演的稳健性;另一方面Q值对低频信息影响小,对高频信息影响大,先用低频反演可以适当忽略Q值的影响,减少了反演的参数,然后再利用高频信息得到精确的纵横波速度、密度和Q值.最后进行了模型测试和实际资料的应用,为之后地震地质层位的标定,成果解释及油气预测奠定了基础,同时也可以为叠前反演提供可靠的约束.受多种因素的影响,实际资料的应用性还有进一步的改进空间.

References
Bian A F, Yu W H, Zhou H W. 2010. Progress in the frequency-domain full waveform inversion method. Progress in Geophysics (in Chinese), 25(3): 982-993.
Brossier R, Operto S, Virieux J. 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion. Geophysics, 74(6): WCC105-WCC118. DOI:10.1190/1.3215771
Brossier R. 2011. Two-dimensional frequency-domain visco-elastic full waveform inversion:Parallel algorithms, optimization and performance. Computers & Geosciences, 37(4): 444-455.
Egorov A, Pevzner R, Bóna A, et al. 2017. Time-lapse full waveform inversion of vertical seismic profile data:Workflow and application to the CO2CRC Otway project. Geophysical Research Letters, 44(14): 7211-7218. DOI:10.1002/2017GL074122
Guitton A, Ayeni G, Díaz E. 2012. Constrained full-waveform inversion by model reparameterization. Geophysics, 77(2): R117-R127. DOI:10.1190/geo2011-0196.1
Jeong W, Lee H Y, Min D J. 2012. Full waveform inversion strategy for density in the frequency domain. Geophysical Journal International, 188(3): 1221-1242. DOI:10.1111/j.1365-246X.2011.05314.x
Li H Q, Li X Y. 2018. Fractured reservoir characterization using VSP and surface data:a case study from the Paris basin. Journal of Geophysics and Engineering, 15(4): 1470-1483. DOI:10.1088/1742-2140/aab288
Li J S, Li Z G, Wu J L. 2019. Establishment of initial velocity field for full waveform inversion based on reflective energy constraints. Progress in Geophysics (in Chinese), 34(4): 1453-1459.
Liang Z Y, Wu G C, Zhang X Y. 2019. Envelope inversion method based on frequency-shifted objective function. Progress in Geophysics (in Chinese), 34(4): 1481-1488.
Liao J P, Liu H X, Dai S X, et al. 2018. Research on 2D elastic wave full waveform inversion in time-space domain-theoretical model numerical experiments. Progress in Geophysics (in Chinese), 33(2): 0671-0678.
Lines L R, Bourgeois A, Covey J D. 1984. Traveltime inversion of vertical seismic profiles-a feasibility study. Geophysics, 49(3): 250-264. DOI:10.1190/1.1441657
Liu P C, Hartzell S, Stephenson W. 1995. Non-linear multiparameter inversion using a hybrid global search algorithm:applications in reflection seismology. Geophysical Journal International, 122(3): 991-1000. DOI:10.1111/j.1365-246X.1995.tb06851.x
Ma Z G, Xie J G. 2005. Relationship among conpressional wave, shear wave velocities and density in rocks. Progress on Geophysics (in Chinese), 20(4): 905-910.
Malinowski M, Operto S, Ribodetti A. 2011. High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full-waveform inversion. Geophysical Journal International, 186(3): 1179-1204. DOI:10.1111/j.1365-246X.2011.05098.x
Operto S, Gholami Y, Brossier R, et al. 2013. A guided tour of multiparameter full-waveform inversion with multicomponent data:from theory to practice. The Leading Edge, 32(9): 1040-1054. DOI:10.1190/tle32091040.1
Owusu J C, Podgornova O, Charara M, et al. 2016. Anisotropic elastic full-waveform inversion of walkaway vertical seismic profiling data from the Arabian Gulf. Geophysical Prospecting, 64(1): 38-53. DOI:10.1111/1365-2478.12227
Pica A, Diet J P, Tarantola A. 1990. Nonlinear inversion of seismic reflection data in a laterally invariant medium. Geophysics, 55(3): 284-292. DOI:10.1190/1.1442836
Prieux V, Brossier R, Operto S, et al. 2013. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field. Part 1: imaging compressional wave speed, density and attenuation. Geophysical Journal International, 194(3): 1640-1664.
Rao Y, Wang Y H. 2016. VSP wave separation by adaptive masking filters. Journal of Geophysics and Engineering, 13(3): 412-421. DOI:10.1088/1742-2132/13/3/412
Štekl I, Pratt R G. 1998. Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators. Geophysics, 63(5): 1779-1794. DOI:10.1190/1.1444472
Stewart R R. 1984. VSP interval velocities from traveltime inversion. Geophysical Prospecting, 32(4): 608-628. DOI:10.1111/j.1365-2478.1984.tb01709.x
Sun M A, Yang J Z, Dong L G, et al. 2017. Density reconstruction in multiparameter elastic full-waveform inversion. Journal of Geophysics and Engineering, 14(6): 1445-1462. DOI:10.1088/1742-2140/aa93b0
Wang X J, Yin X Y, Wu G C. 2011. Estimation of stratigraphic quality factors on pre-stack seismic data. Oil Geophysical Prospecting (in Chinese), 46(3): 423-428. DOI:10.1007/s12182-011-0118-0
Wei X, Yang Z F, Yan X F. 2019. Multi-scale issue in seismic exploration and its research progress. Progress in Geophysics (in Chinese), 34(6): 2353-2360.
Wu Y T, Zhu G M. 2017. Comparison and application of three different methods for inversing Q value with zero-offset VSP data. Progress in Geophysics (in Chinese), 32(5): 2107-2114.
Yang J Z, Liu Y Z, Dong L G. 2014. A multi-parameter full waveform inversion strategy for acoustic media with variable density. Chinese Journal of Geophysics (in Chinese), 57(2): 628-643.
Yang J Z, Liu Y Z, Dong L G. 2016. Simultaneous estimation of velocity and density in acoustic multiparameter full-waveform inversion using an improved scattering-integral approach. Geophysics, 81(6): R399-R415. DOI:10.1190/geo2015-0707.1
Yang Q Y, Hu G H, Wang L X. 2014. Research status and development trend of full waveform inversion. Geophysical Prospecting for Petroleum (in Chinese), 53(1): 77-83.
Yao Z Y, He X H. 2007. Calculation of S-wave velocity from zero-offset VSP data:Taking the VSP data of well 101 in Qibei survey area as an example. Progress in Exploration Geophysics (in Chinese), 30(2): 100-103.
Yin X Y, Zong Z Y, Wu G C. 2015. Research on Seismic fluid identification driven by rock physics. Science China Earth Sciences, 58(2): 159-171. DOI:10.1007/s11430-014-4992-3
Zhang G Z, Jiang L J, Sun C L, et al. 2017. The stepped multi-parameter FWI of acoustic media in time-domain by L-BFGS method with illumination analysis. Geophysical Prospecting for Petroleum (in Chinese), 56(1): 31-37, 74.
Zhang J G, Mu F M. 2006. Study of practicality of VSP S-wave velocity inversion. Oil Geophysical Prospecting (in Chinese), 41(6): 697-701.
Zong Z Y, Yin X Y, Wu G C. 2012a. Fluid identification method based on compressional and shear modulus direct inversion. Chinese Journal of Geophysics (in Chinese), 55(1): 284-292.
Zong Z Y, Yin X Y, Zhang F, et al. 2012b. Reflection coefficient equation and pre-stack seismic inversion with Young's modulus and Poisson ratio. Chinese Journal of Geophysics (in Chinese), 55(11): 3786-3794.
Zou Q, Zhou X X, Zhong B S. 2003. Acquiring interval velocity by joint inversion of travel-time and wave impedance. Oil Geophysical Prospecting (in Chinese), 38(4): 396-399.
李冀蜀, 李志刚, 吴建鲁. 2019. 基于反射能量约束全波形反演初始速度场建立研究. 地球物理学进展, 34(4): 1453-1459.
梁展源, 吴国忱, 张晓语. 2019. 基于频移目标函数的包络反演方法. 地球物理学进展, 34(4): 1481-1488.
廖建平, 刘和秀, 戴世鑫, 等. 2018. 二维时间空间域弹性波全波形反演-理论模型数值试验. 地球物理学进展, 33(2): 0671-0678.
马中高, 解吉高. 2005. 岩石的纵、横波速度与密度的规律研究. 地球物理学进展, 20(4): 905-910. DOI:10.3969/j.issn.1004-2903.2005.04.004
王小杰, 印兴耀, 吴国忱. 2011. 基于叠前地震数据的地层Q值估计. 石油地球物理勘探, 46(3): 423-428.
未晛, 杨志芳, 晏信飞. 2019. 地震勘探中的多尺度问题及其研究进展. 地球物理学进展, 34(6): 2353-2360.
武银婷, 朱光明. 2017. 基于零偏VSP三种Q值反演方法对比分析及应用. 地球物理学进展, 32(5): 2017-2114.
杨积忠, 刘玉柱, 董良国. 2014. 变密度声波方程多参数全波形反演策略. 地球物理学报, 57(2): 628-643.
杨勤勇, 胡光辉, 王立歆. 2014. 全波形反演研究现状及发展趋势. 石油物探, 53(1): 77-83. DOI:10.3969/j.issn.1000-1441.2014.01.011
姚忠瑞, 何惺华. 2007. 利用零偏VSP资料求取横波速度-以七北101井VSP测井为例. 勘探地球物理进展, 30(2): 100-103.
印兴耀, 宗兆云, 吴国忱. 2015. 岩石物理驱动下地震流体识别研究. 中国科学:地球科学, 45(1): 8-21.
张广智, 姜岚杰, 孙昌路, 等. 2017. 基于照明预处理的分步多参数时间域声波全波形反演方法研究. 石油物探, 56(1): 31-37, 74. DOI:10.3969/j.issn.1000-1441.2017.01.004
张继国, 牟风明. 2006. VSP横波速度反演实用性研究. 石油地球物理勘探, 41(6): 697-701. DOI:10.3321/j.issn:1000-7210.2006.06.017
宗兆云, 印兴耀, 吴国忱. 2012a. 基于叠前地震纵横波模量直接反演的流体检测方法. 地球物理学报, 55(1): 284-292.
宗兆云, 印兴耀, 张峰, 等. 2012b. 杨氏模量和泊松比反射系数近似方程及叠前地震反演. 地球物理学报, 55(11): 3786-3794.
邹强, 周熙襄, 钟本善. 2003. 走时与波阻抗联合反演求取层速度. 石油地球物理勘探, 38(4): 396-399. DOI:10.3321/j.issn:1000-7210.2003.04.011