地球物理学进展  2016, Vol. 31 Issue (4): 1678-1687   PDF    
全波形反演初始模型建立策略研究综述
王连坤1,2, 方伍宝1, 段心标1, 曹文俊2, 李振春2     
1. 中国石油化工股份有限公司石油物探技术研究院, 南京 211103
2. 中国石油大学(华东)地球科学与技术学院, 青岛 266580
摘要: 全波形反演充分利用叠前地震波场中的运动学和动力学信息估计地球内部介质的弹性参数,是一种高精度的速度建模方法和油藏描述的手段.但是对于实际资料,尤其是陆地资料,全波形反演很难收敛到正确的结果上,往往无法实现其建立高精度速度模型的目标,需要较准确的初始模型来降低误差泛函的非线性性.研究与调研认为,建立全波形反演初始模型的方法主要有三种:走时层析、偏移速度分析和Laplace域反演.本文首先分析了全波形反演对初始模型的严重依赖性,其次对全波形反演初始模型建立策略进行简要的介绍与分析,最后给出了建立全波形反演初始模型的认识及建议.
关键词全波形反演     初始模型     走时层析反演     Laplace域反演     特征波反演    
Review of full waveform inversion initial model building strategy
WANG Lian-kun1,2 , FANG Wu-bao1 , DUAN Xin-biao1 , CAO Wen-jun2 , LI Zhen-chun2     
1. SINOPEC Geophysical Research Institute, Nanjing 211103, China
2. School of Geoscience, China University of Petroleum (East China), Qingdao 266580, China
Abstract: Full waveform inversion makes full use of the kinematic and dynamic information,which are in the prestack seismic data, to estimate the elasticity parameters of the earth's interior medium, and is a kind of high-accuracy velocity model building method and means of reservoir description. But for the actual data, especially the land data, full waveform inversion is difficult to converge to the correct results, it often can not reach the goal of building the high-accuracy velocity model, this is because the relationship between the seismic wave field (especially, reflection wave field) and the geophysical parameters is strongly nonlinear, the error function of the FWI exists a huge number of local extreme value point, and this will seriously affect the convergence results of the inversion, so we need a high-accuracy initial model to reduce the nonlinear of the error function in the actual application. From the inversion strategic of the FWI, a prerequisite condition for the FWI successful implementation is to estimate a macro velocity model with higher accuracy. Of course, how to establish an accurate initial model is also a research hotspot in the field of exploration. After a long research and investigation, we think that the building methods of the full waveform inversion initial model basically has the following three kinds: traveltime tomography inversion, migration velocity analysis and Laplace domain inversion. This paper first briefly analyses the full waveform inversion serious dependence on initial model, then introduces and analyses briefly the full waveform inversion initial model building strategy, finally gives the suggestions to build the full waveform inversion initial model.
Key words: Full waveform inversion     Initial model     Travel time tomography inversion     Laplace domain inversion     Characteristic wave inversion    
0 引言

全波形反演(full waveform inversion, FWI)充分利用叠前地震波场中的运动学和动力学信息估计地球内部介质的弹性参数, 是一种高精度速度建模的方法(Virieux and Operto, 2009; 刘国峰等, 2012)和油藏描述的手段(王西文等, 2013).从地球物理的角度看, FWI涉及模型参数化、正演模拟、误差泛函构建、数据预处理、优化方法的选取、正则化约束、迭代终止条件、子波估计等多项研究内容(胡光辉等, 2013).从某种程度上讲, FWI是一套完美的理论体系, 它考虑了极其丰富的波现象(如一次反射、多次反射、透射、绕射、散射等), 其数据空间是地震全波场信息, 在理论上是具有最高反演分辨率的速度估计方法.

但是, 从数学角度来看, FWI却是一个高度病态的非线性问题(Tarantola, 1984), 需要解决其收敛性和多解性问题(杨午阳等, 2013).理论模型测试结果表明, FWI可以建立高精度的速度模型, 但是对于实际资料, 尤其是陆地资料, FWI很难收敛到正确的结果上, 往往无法实现其建立高精度速度模型的目标.这是因为地震波波场(尤其是反射波场)与地球物理参数之间的关系是强非线性的, 误差泛函存在非常多的局部极值点, 这会严重影响FWI的收敛结果(Gauthier et al., 1986), 所以在实际应用时往往需要一个较好的宏观速度模型(速度模型的中低波数部分)来降低误差泛函的非线性性(Virieux and Operto, 2009).

反演策略上认为, FWI成功实施的前提条件是估计一个精度较高的宏观速度模型.当然, 如何建立一个精确的初始模型也是FWI领域的一个研究热点.目前, 建立初始模型的方法很多, 在油气勘探中, 可以使用反射波层析成像和基于偏移的速度分析方法(migration velocity analysis, MVA)来建立FWI的初始模型(Vireux and Operto, 2009).经过长时间的研究与调研认为, 建立FWI初始模型的方法主要有以下三种:走时层析(Bishop et al., 1985; Lambaré, 2008; 成谷等, 2007; 刘玉柱等, 2014)、偏移速度分析(Al-yahya, 1989; Symes and Carazzone, 1991; 李振春, 2002; 刘守伟等, 2008)和Laplace域反演(Shin and Cha, 2008; Shin and Ha, 2008; Shin et al, 2013; Ha et al, 2015; 胡英等, 2015).本文首先分析了FWI对初始模型的严重依赖性, 其次对FWI初始模型建立策略进行简要的介绍与分析, 最后给出了认识及研究建议.

1 全波形反演对初始模型的严重依赖性

在当前计算能力飞速发展的环境下, FWI已经从2D走向了3D, 从模型测试验证逐步走向了实用化阶段. 但是, 就目前来说, FWI的实际应用还仅仅局限在针对海上资料, 以挪威Valhall油田为代表的海上资料应用陆续出现, 并取得了很好的建模效果, 但是陆地资料还很难实现. 究其原因, 相对于海上地震资料, 陆地资料很难实现全方位、大偏移距的观测; 低频信息缺失; 震源稳定性较差, 子波未知; 随机干扰严重; 正演模拟波场与实际波场的差异过大; 观测数据的照明角度不足等(王华忠等, 2012).其中, 低频信息与初始速度模型的耦合是全波形反演在实际资料应用中的最大瓶颈(Sirgue, 2006; 杨勤勇等, 2014).

梯度类的全波形反演算法使用基于局部梯度的最优化方法实现观测波场和模拟波场的残差最小, 该算法最关键的制约就是FWI对初始模型的严重依赖性(Virieux and Operto, 2009).低频信息的缺失, 使得常规的速度建模方法难以满足FWI对初始模型精度的要求(Symes, 2008), 正演合成的数据与实测数据中对应的同相轴到达时差异过大(超过所有频段的半个波长), 就会引起周波跃迁现象(Hu, 2012), 造成地震反演对长波长和中间波长信息的低敏感性(Jannane et al., 1989), 使得反演陷入局部极值, 造成模型更新迭代过程中的错误收敛(Bunks et al., 1995). 只有初始模型能够较为准确地反映地震波走时, FWI才能在迭代过程中收敛到真实解附近. 因此, 一个精确的初始模型显得尤为重要(Virieux and Operto, 2009; Prieux et al., 2013).

20世纪80年代, Gauthier等(1986) 等学者率先实现了2D地震资料的时间域全波形反演.随后, Pratt(1990) Geller和Hara(1993) Song等(1995) Ravaut等(2004) Gao等(2006) Operto等(2006) Bleibinhaus等(2007) Sirgue等(2010) Fichtner和Trampert(2011) 也实现了FWI的应用, 这些应用表明, FWI是一种高精度的速度建模方法, 它具有精确刻画地下构造细节的潜力, 同时也指出FWI严重依赖初始速度模型.

同济大学的任浩然(2012) 以经典的洼陷模型(图 1)为例, 测试了FWI对三种初始模型(光滑速度模型, 常速模型, 常梯度模型)的依赖性, 其中常速模型的速度值为3500 m/s, 它在第二层的速度偏离了原值的12.5%; 常梯度模型的速度从浅到深从3500 m/s到4500 m/s按照常梯度变化, 在浅层的速度值错误, 垂直入射至凹陷底部的走时差约为7.1%.图 2给出了这三种模型的测试结果.从图 2中可以看到, 当初始速度偏离真实速度较多时, FWI难以收敛到真实解附近. 在本次测试中, 尤其是常速模型在更新最下层位置时显得无能为力; 尽管经过151次迭代, 下层界面的位置仍不能收敛到正确位置, 而是在一个假的位置上“徘徊”, 即反演过程错误收敛, 陷入了局部极小值.常梯度模型在凹陷部位同样出现错误收敛.这再次说明, FWI对初始模型准确度的要求是非常高的.

图 1 洼陷速度模型(任浩然, 2012) (a)真实速度模型; (b)平滑速度模型. Figure 1 Sag velocity model(Ren, 2012) (a)The real velocity model; (b)The smooth velocity model.
图 2 全波形反演对初始模型的依赖性测试结果(任浩然, 2012) (a)为平滑速度模型反演结果; (b)为常速模型反演结果; (c)为常梯度模型反演结果. Figure 2 FWI dependence on initial model test results(Ren, 2012) (a)The smooth velocity model FWI result;(b)The constant velocity model FWI result;(c)The constant gradient velocity model FWI result.
2 走时层析

地震层析成像是一种常用的速度建模的方法, 它主要利用以地震波走时为代表的运动学信息进行速度建模, 属于数据拟合类的速度估计方法(刘福田等, 1989; 成谷等, 2007; 刘玉柱等, 2014).目前, 勘探地球物理领域常用的层析成像方法主要有反射层析、初至走时层析、立体层析等, 下面分别对这三种方法的原理以及在FWI领域的应用进行简单地介绍.

2.1 反射层析

反射层析(Bishop et al., 1985; Farra and Madariaga, 1988; Yang et al., 2005; Woodward et al., 2008)是最经典的层析成像方法, 其模型空间定义为速度分布和一些连续分布的反射界面, 利用反射射线追踪的模拟走时和拾取走时的残差建立误差泛函, 误差泛函的最优化问题可以归结为线性最小二乘优化问题.Bishop等(1985) 最先利用反射地震数据进行层析研究, Farra和Madariaga(1988) 通过补偿函数的方式加入了先验信息, 并用一个衰减因子加速收敛.后来又有很多学者进行了研究(Stork, 1992; Delprat-Jannaud and Lailly, 1995; Clapp et al., 2004; Belkebir et al., 2005; 秦宁等, 2012), 但是反射层析需要在地震波场中沿着连续的同相轴拾取反射同相轴, 但是这在具有波形复杂和低信噪比的实际叠前地震数据中是很难办到的!此外, 反射层析的正问题为两点反射射线追踪, 在3D情况下实施极为困难, 这在一定程度上影响了反射层系的稳定性.

反射层析作为一种经典的速度建模方法, 它在全波形反演方面的应用还是很多的. Sirgue等(2010) 在挪威Valhall油田进行了3D频率域全波形反演的研究和应用, 其初始模型是由基于射线的反射层析得到, 反演过程中实现了150 m处古河道和1050 m处气云的准确成像, 如图 3, 图 3a3c是反射层析方法的结果, 图 3b3d是FWI的结果, 3a、3b为150 m的深度切片, 3c、3d为1050 m的深度切片.仅仅通过比较两种方法(基于射线的反射层析和全波形反演)的速度深度切片, 可以看出全波形反演在构造细节方面的优势, 明显提高了分辨率.Zelt等(2005) Operto等(2006) Brenders and Pratt(2006) Ma和Hale(2013) 等也进行了相关的研究.Prieux等(2013) 指出, 虽然反射层析充分利用反射波走时构建高精度的宏观速度模型, 但是它不适用于直达波, 而直达波对全波形反演来说却是至关重要的!

图 3 Valhall地区FWI结果(Sirgue et al., 2010) (a)为反射层析得到的速度切片(150 m); (b)为FWI得到的速度切片(150 m); (c)为反射层析得到的速度切片(1050 m); (d)为FWI得到的速度切片(1050 m). Figure 3 Imaging of the Valhall field by 3D FWI(Sirgue et al., 2010) (a)The velocity slice form reflection tomography(150 m);(b)The velocity slice form FWI(150 m);(c)The velocity slice form reflection tomography(1050 m);(d)The velocity slice form FWI(1050 m).
2.2 立体层析

鉴于反射层析方法数据存在拾取困难, Billette和Lambaré(1998)Billette等(2003) Duveneck(2004) Lambaré(2008) 等发展了立体层析方法.它不仅利用了地震波走时, 还充分考虑了共炮点和共检波点道集的局部相干同相轴的斜率以及炮检点位置范围, 其模型空间则被定义为速度分布、反射点位置、反射角与入射角以及反射点到炮点和检波点的单程走时(倪瑶等, 2013), 其拾取对象为地震记录中具有局部相干性的反射波同相轴, 并将其解释为一个反射/绕射过程, 利用地震波走时和地震波在地表的局部传播方向约束速度模型, 较反射层析具有更强的稳定性.另外, 此方法还在数据空间中引入了炮检点处走时的一阶导数(分别对应于炮检波点处慢度矢量的水平分量), 这样, 在进行数据拾取时就不需要沿着连续的反射界面进行, 这相对于反射层析的连续拾取是更加合理的, 而数据空间中新引入的炮检点处旅行时的一阶导数可以有效地约束炮检点处射线的局部传播方向, 很大程度上避免了射线多路径现象(胡光辉等, 2014).

立体层析方法可以认为是射线理论意义下比较完善的一种速度建模方法, 很多学者对立体层析进行了系统的研究, Chauris和Noble(2001) 认识到了立体层析与DSQ(differential semblance optimization)的联系, Lambaré等(2004) 发展了一种半自动的初至拾取方法, 倪瑶等(2013) 利用射线中心坐标系下的射线扰动理论与射线传播矩阵的性质推导出了立体层析Frechet导数的简洁形式, 促进了立体层析速度建模的进一步发展.同时, 也有很多学者(Billette和Lambaré., 1998; Billette et al., 2003; Lambaré and Alerini, 2005; Prieux et al., 2009, 2013)利用立体层析来建立FWI的初始模型, 进行全波形反演的应用研究.

2.3 初至走时层析

初至走时层析(first-arrival traveltime tomography, FATT)利用直达波、折射波、回折波、透射波等不同类型的初至信息进行非线性反演, 以此来产生平滑的模型(Hole, 1992; Zelt and Barton, 1998; 秦宁等, 2012, 2014; 桑运云等, 2013).走时残差沿着射线路径反向传播, 以此来计算灵敏度矩阵, 采用LSQR(Paige and Saunders, 1982)求解, 也可以采用伴随状态法求解, 该方法避开了灵敏度矩阵的直接求解, 就像在FWI中一样(Taillandier et al, 2009).

将FATT得到的初始模型用于FWI的例子很多, 例如, Ravaut等(2004) Jaiswal等(2008, 2009)、Malinowski和Operto(2008) Prieux等(2009, 2011)将其应用于地面观测的资料中; Dessa和Pascal(2003) 将其应用与实验室超声实验中; Pratt和Goulty(1991) 将其应用于井间资料中; Gao等(2006) 将其应用到了VSP资料中.

Brenders和Pratt(2006, 2007)研究发现, 当使用FATT获得的初始模型进行FWI时, 为了得到可靠的反演结果, 必须要有低频信息和大偏移距资料.例如, 只有当Brenders和Pratt(2007) 使用初始频率为0.5 Hz和最大偏移距为16 km进行全波形反演时, BP基准模型上部区域才可能准确成像.此外, 当地下介质中存在低速体时, 可靠的初至拾取也是很困难的, FATT是不适用的(Prieux et al., 2011). Ellefsen(2009) 认为, 在使用复数频率值时, 用一个频率域波形反演算法可以把FATT转化为初至相位反演, 与基于高频近似的FATT相比, 该方法有助于解释在层析灵敏度核函数中的有效频率效应.

Prieux等(2009) 分别用FATT和立体层析建立了FWI的初始模型, 然后用FWI重新对合成的Valhall模型进行测试, 其结果如图 4, 图 4a为真实的Valhall速度模型, 图 4b为真实模型平滑后进行FWI的结果, 图 4c为FATT方法建立初始模型后进行FWI的结果, 图 4d为立体层析方法建立初始模型后FWI的结果.测试结果可以看出, 立体层析方法成功恢复了速度的长波长信息, 而利用FATT得到的初始模型进行FWI之后, 模型的浅层部分得到了精确的恢复.

图 4 FWI对初始模型的灵敏性测试(Prieux et al., 2009) (a)为真实Valhall模型; (b)为利用平滑速度模型得到的FWI结果; (c)为利用FATT得到的速度模型进行FWI的结果; (d)为利用立体层析得到的速度模型进行FWI的结果. Figure 4 FWI sensitivity testing of initial model(Prieux et al., 2009) (a)The real Valhall velocity model;(b)The FWI result by smooth velocity model; (c)The FWI result by velocity from FATT;(d)The FWI result by velocity from stereo-tomography.
3 偏移速度分析

一般来说, 偏移后的成像数据体比起叠前数据来说, 信噪比更高, 绕射波得到收敛, 同相轴的相干性更好, 波形更简单.这些优点使得偏移速度分析(migration velocity analysis, MVA)在反演数据空间的提取方面较层析成像有较大的优势(李振春, 2002; 陈志德等, 2002; 刘守伟等, 2008). Al-Yahya(1989) 是比较早的一位将剩余曲率分析引入到叠前偏移速度分析的地球物理学家, 他基于常速介质假设推导出了一种迭代的速度分析方法, 该方法的实现过程是:首先用初始速度模型进行叠前深度偏移, 然后估计速度误差, 并且更新初始速度模型, 重复这一操作, 直至误差小于某一条件, 该方法实现起来比较稳定, 但需要在每次更新模型后实施偏移并拾取共成像点道集上的深度差, 并且当速度横向变化剧烈时, 这种方法精度很差、计算效率非常低.Liu(1997) 基于射线理论, 推导出了成像深度关于速度的导数的解析形式, 该导函数揭示了剩余时差和剩余速度的关系, 然后Liu利用该导函数和扰动理论, 导出了一个速度分析的新方法, 该方法可以MVA的灵敏度和误差估计, 有助于对速度估计的可靠性进行定量分析.该方法虽然在一定程度上缓解了MVA方法对速度模型复杂度的严格限制, 但依然需要进行反射界面拾取.

自动拾取的MVA方法充分考虑了CIG道集的相干性, 通过构造目标泛函的方式来刻画CIG道集是否拉平.Jin和Madariaga(1994) 发展了一种两步Monte Carlo非线性速度迭代反演算法, 由于其相干性的计算是在模型空间而不是数据空间进行的, 故该算法成功避开了需要大量计算的正演模拟, 在很大程度上提高了计算效率.Jervis等(1996) 使用了模拟退火法和遗传算法来进行偏移速度估计, 并将其应用到了墨西哥湾实测数据中.但是这些方法目标泛函都是强非线性的, 具有许多次级极值, 如果采用全局最优化的方式寻优, 那么将面对非常巨大的计算量, 需要高昂的计算成本; 但是如果以局部最优化的方式寻优, 则容易陷入局部极小值, 难以取得成功.因此, 这类优化问题陷入了一个两难的局面.

Symes和Carazzone(1991) 通过数值试验提出了DSQ泛函的概念, 并进行了理论性的研究.Chauris和Noble(2001) 发展了基于DSQ泛函的反射地震数据宏观速度估计方法, 他们指出, DSQ泛函较其他的泛函而言具有良好的凸性, 是一种自动化的方法, 既不需要拾取反射界面, 也不需要解释, 并且可以通过局部(梯度)优化的方式收敛, 从而更加适合以局部最优化方式进行寻优.Ling等(2014) 针对山前带地区信噪比较低的现状, 指出在复杂地表情况下很难实现常规的速度反演方法, 并在速度模型空间发展了一种基于马尔科夫链的蒙特卡洛背景速度反演方法, 在实用化过程中, 他们使用了Kirchhoff叠前深度偏移确定误差泛函(比如DSQ), 但是仍然需要进一步提高反演效率.

近年来, 随着高斯束偏移、逆时偏移的发展, MVA方法也得到了进一步发展(李振春, 2011), 但是国内外研究学者很少利用MVA方法来建立全波形反演的初始模型, MVA方法还是主要用于偏移成像领域.SEP的Biondi和Almomin(2012, 2013)联合MVA和FWI, 提出了TFWI(tomographic full waveform inversion), 并对Marmousi模型进行了测试, 其结果如图 5, 图 5a为真实速度模型, 图 5b为初始模型, 图 5c为FWI的结果, 图 5d为TFWI的结果, 数值测试表明TFWI对初始模型的精度要求不是很高, 但也能很好地全局收敛, 并展示出了很好精确性.

图 5 Marmousi模型FWI测试结果(Biondi and Almomin, 2012) (a)为真实速度模型; (b)为初始模型; (c)为FWI的结果; (d)为TFWI的结果. Figure 5 FWI test results of Marmousi model(Biondi and Almomin, 2012) (a)the real velocity model;(b)the initial velocity model;(c)the FWI result;(d)the TFWI result.
4 Laplace域波形反演

Laplace域波形反演(Shin and Cha, 2008; Shin and Ha, 2008; Shin et al., 2013; Ha et al., 2015)充分利用了Laplace域对频率不敏感的特性, 可以利用简单的初始模型从缺失低频信息的地震资料中得到长波长速度模型(曹书红和陈景波, 2014).该算法的实现过程是:利用衰减波形记录的直流分量进行反演, 产生一个对初始模型不敏感的宏观模型, 以此作为初始模型参与频率域全波形反演(Shin and Cha, 2008; 胡英等, 2015).它可以被视为使用复数频率值的频率域波形反演, 其实部为0, 虚部控制着地震波场在时间域上的衰减, 它对初始模型敏感性差一些.Laplace域反演建立了一个平滑的速度模型(长波长速度信息).在Shin和Cha(2008) 的研究中, 该算法得到的速度模型可以作为频率域全波形反演的低频速度成分.Laplace域反演比起时间域反演和频率域反演, 在有限的计算资源下, 具有高效求解速度模型的能力(Pyun et al., 2011), Laplace域正演模拟可以选用大网格间隔, 这使得阻抗矩阵“变小”了, 而且频散可以忽略(Shin et al., 2014).

虽然Laplace域波形反演可以恢复速度的长波长信息, 避开低频缺失带来的影响, 但是Laplace域反演深度依赖于偏移距和衰减参数(Shin and Cha, 2009).在分析了Laplace域反演这一缺点后, Shin和Cha(2009) 发展了Laplace-Fourier波形反演算法, 它充分利用了阻尼波场的低频分量.Laplace-Fourier域相当于对沿时间方向衰减的地震图进行一次反演.Shin和Cha(2009)的研究表明, 该方法可以应用于频率小于地震带宽的最小频率的情况下, 它的这种能力被Shin和Cha(2009)称为是低频的“幻影复活”(mirage resurrection). 图 6是在BP标准模型中联合应用Laplace域、Laplace-Fourier域和Fourier域FWI的案例(Shin和Cha, 2009).初始模型是一个简单的速度梯度模型(图 6b), 然后用Laplace域反演得到第一个长波长的速度模型(图 6c), 接下来用该模型作为初始模型进行Laplace-Fourier域反演(图 6d).在这里, 在衰减数据的反演中所使用的起始频率低至0.01Hz.最后进行频率域FWI, 其结果见图 6e, 模型中所有的结构均被成功成像.

图 6 BP模型FWI测试结果(Shin and Cha, 2009) (a)为真实速度模型; (b)为初始速度模型; (c)由Laplace域反演得到的速度模型; (d)为对(c)速度模型进行Laplace-Fourier域反演的结果; (e)为对(d)进行频率域FWI的结果. Figure 6 FWI test results of BP model(Shin and Cha, 2009) (a)The real velocity model;(b)The initial velocity model;(c)The Laplace domain inversion result by the initial velocity model;(d)The Laplace-Fourier domain inversion result by the velocity model form(c); (e)The frequency domain FWI result by the velocity model form(c).

虽然Shin等人开展了一些关于Laplace-Fourier 域反演的研究(Shin and Cha, 2009), 但总体上, Laplace-Fourier域全波形反演的理论和方法研究尚不充分.Laplace-Fourier域全波形反演需要选择不同的衰减常数与多个频率点组合进行反演计算, 其计算效率低, 占用内存大, 运算时间长.基于此, 胡英等(2015) 在前人研究的基础上, 结合多尺度网格和多尺度计算区域两种方法, 发展了一种Laplace-Fourier域多尺度高效全波形反演方法, 它通过调整网格间距大小和模型的计算区域深度减少反演网格数, 频率点的高低决定反演网格间距大小, Laplace衰减常数大小决定反演模型计算深度, 提高反演效率, 在反演精度保持不变的情况下, 有效提高了Laplace-Fourier域的计算效率.

5 认识与展望

5.1 能否建立一个符合FWI精度要求的初始模型是FWI成功实施的关键因素, 上文分别从走时层析、偏移速度分析、Laplace域反演等方面阐述了FWI初始模型的建立策略. 但在实际资料应用中也很难做到初始模型与初始最低频率的耦合, 因为在常规的地震勘探资料采集中往往缺少可靠的低频成分, 应采用初至波走时层析为速度场提供可靠的浅层低频低波数成分, 中深层的中低频中低波数成分可通过测井资料等先验信息约束的反演提供, 二者综合建模后再利用网格层析技术进一步校正速度场, 为FWI提供一个较为可靠的初始速度模型(胡光辉等, 2014).

5.2 此外, 特征波波形反演在反演效果上要优于全波形反演(王华忠等, 2012). 特征波波形反演使用某一或某一些特征波的走时信息和波形信息进行匹配, 而不是匹配所有波场, 改善了反演的结果, 提高了反演的精度, 是解决陆上资料FWI问题的有效途径.基于偏移/反偏移的反射波全波形反演(Xu et al., 2012)降低了对低频信息的依赖, 实际上是特征波波形反演的一个特例.

5.3 FWI是一个系统工程, 它的成功实施必须依赖于可靠的初始速度模型来降低其误差泛函的非线性性.而现阶段常用的FWI初始速度建模方法(走时层析、偏移速度分析、Laplace域反演)都有其固有的缺陷和适用范围, 仅仅使用一种方法很难得到满足精度要求的初始速度模型, 因此就需要综合多种建模手段, 充分发挥每种建模手段的优势. 此外, 特征波波场反演和构建新的目标泛函(董良国等, 2013)也是非常重要的思路和研究方向. 总之, FWI实用化的进程漫长而曲折, 但相信有一天FWI会发挥其高精度成像、多参数建模、进行储层描述的优势, 促进地球物理勘探技术的继续发展.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持!
参考文献
[1] Al-Yahya K.1989. Velocity analysis by iterative profile migration[J]. Geophysics, 54 (6) : 718–729. DOI:10.1190/1.1442699
[2] Belkebir K, Chaumet P C, Sentenac A.2005. Superresolution in total internal reflection tomography[J]. Journal of the Optical Society of America A, 22 (9) : 1889–1897. DOI:10.1364/JOSAA.22.001889
[3] Billette F, Lambaré G.1998. Velocity macro-model estimation from seismic reflection data by stereotomography[J]. Geophysical Journal International, 135 (2) : 671–690. DOI:10.1046/j.1365-246X.1998.00632.x
[4] Billette F, Le Bégat S, Podvin P, et al.2003. Practical aspects and applications of 2D stereotomography[J]. Geophysics, 68 (3) : 1008–1021. DOI:10.1190/1.1581072
[5] Biondi B, Almomin A. 2012. Tomographic full waveform inversion (TFWI) by combining full waveform inversion with wave-equation migration velocity analysis[C].// 82nd Annual International Meeting, SEG. Expanded Abstracts, 547-552.
[6] Biondi B, Almomin A.2013. Tomographic full-waveform inversion (TFWI) by combining FWI and wave-equation migration velocity analysis[J]. The Leading Edge, 32 (9) : 1074–1080. DOI:10.1190/tle32091074.1
[7] Bishop T N, Bube K P, Cutler R T, et al.1985. Tomographic determination of velocity and depth in laterally varying media[J]. Geophysics, 50 (6) : 903–923. DOI:10.1190/1.1441970
[8] Bleibinhaus F, Hole J A, Ryberg T, et al.2007. Structure of the California coast ranges and San Andreas Fault at SAFOD from seismic waveform inversion and reflection imaging[J]. Journal of Geophysical Research: Solid Earth, 112 (B6) : B06315.
[9] Brenders A J, Pratt R G.2006. Full waveform tomography for lithospheric imaging: Results from a blind test in a realistic crustal model[J]. Geophysical Journal International, 168 (1) : 133–151.
[10] Brenders A J, Pratt R G.2007. Efficient waveform tomography for lithospheric imaging: Implications for realistic, two-dimensional acquisition geometries and low-frequency data[J]. Geophysical Journal International, 168 (1) : 152–170. DOI:10.1111/gji.2007.168.issue-1
[11] Bunks C, Saleck F M, Zaleski S, et al.1995. Multiscale seismic waveform inversion[J]. Geophysics, 60 (5) : 1457–1473. DOI:10.1190/1.1443880
[12] Cao S H, Chen J B.2014. Studies on complex frequencies in frequency domain full waveform inversion[J]. Chinese Journal of Geophysics (in Chinese), 57 (7) : 2302–2313. DOI:10.6038/cjg20140724
[13] Chauris H, Noble A M.2001. Two-dimensional velocity macro model estimation from seismic reflection data by local differential semblance optimization: Applications to synthetic and real data sets[J]. Geophysical Journal International, 144 (1) : 14–26. DOI:10.1046/j.1365-246x.2001.00279.x
[14] Chen Z D, Liu Z K, Li C B.2002. 3-D pre-stack depth migration velocity analysis and automatic Monte Carlo velocity picking in depth[J]. Chinese Journal of Geophysics (in Chinese), 45 (2) : 246–254. DOI:10.3321/j.issn:0001-5733.2002.02.011
[15] Cheng G, Zhang B J, Ma Z T.2007. Regularization in traveltime tomography for reflection seismic data[J]. Progress in Geophysics (in Chinese), 22 (6) : 1715–1721. DOI:10.3969/j.issn.1004-2903.2007.06.005
[16] Clapp R G, Biondi B L, Claerbout J F.2004. Incorporating geologic information into reflection tomography[J]. Geophysics, 69 (2) : 533–546. DOI:10.1190/1.1707073
[17] Delprat-Jannaud F, Lailly P.1995. Reflection tomography: How to handle multiple arrivals?[J]. Journal of Geophysical Research, 100 (B1) : 703–715. DOI:10.1029/94JB02461
[18] Dessa J X, Pascal G.2003. Combined traveltime and frequency-domain seismic waveform inversion: A case study on multi-offset ultrasonic data[J]. Geophysical Journal International, 154 (1) : 117–133. DOI:10.1046/j.1365-246X.2003.01956.x
[19] Dong L G, Chi B X, Tao J X, et al.2013. Objective function behavior in acoustic full-waveform inversion[J]. Chinese Journal of Geophysics (in Chinese), 56 (10) : 3445–3460. DOI:10.6038/cjg20131020
[20] Duveneck E.2004. Velocity model estimation with data-derived wavefront attributes[J]. Geophysics, 69 (1) : 265–274. DOI:10.1190/1.1649394
[21] Ellefsen K J.2009. A comparison of phase inversion and traveltime tomography for processing near-surface refraction traveltimes[J]. Geophysics, 74 (6) : WCB11–WCB24. DOI:10.1190/1.3196857
[22] Farra V, Madariaga R.1988. Non-linear reflection tomography[J]. Geophysical Journal International, 95 (1) : 135–147. DOI:10.1111/j.1365-246X.1988.tb00456.x
[23] Fichtner A, Trampert J.2011. Resolution analysis in full waveform inversion[J]. Geophysical Journal International, 187 (3) : 1604–1624. DOI:10.1111/gji.2011.187.issue-3
[24] Gao F C, Levander A R, Pratt R G, et al.2006. Waveform tomography at a groundwater contamination site: VSP-surface data set[J]. Geophysics, 71 (1) : H1–H11. DOI:10.1190/1.2159049
[25] Gauthier O, Virieux J, Tarantola A.1986. Two-dimensional nonlinear inversion of seismic waveforms: Numerical results[J]. Geophysics, 51 (7) : 1387–1403. DOI:10.1190/1.1442188
[26] Geller R J, Hara T.1993. Two efficient algorithms for iterative linearized inversion of seismic waveform data[J]. Geophysical Journal International, 115 (3) : 699–710. DOI:10.1111/gji.1993.115.issue-3
[27] Ha W, Kang S G, Shin C.2015. 3D Laplace-domain waveform inversion using a low-frequency time-domain modeling algorithm[J]. Geophysics, 80 (1) : R1–R13. DOI:10.1190/geo2013-0332.1
[28] Hole J A.1992. Nonlinear high-resolution three-dimensional seismic travel time tomography[J]. Journal of Geophysical Research, 97 (B5) : 6553–6562. DOI:10.1029/92JB00235
[29] Hu G H. 2012. Three-dimensional acoustic Full Waveform Inversion: Method, algorithms and application to the Valhall petroleum field[Ph. D. thesis]. Paris: Université De Grenoble.
[30] Hu Guanghui, Jia Chunmei, Xia Hongrui, et al.2013. Implementation and validation of 3D acoustic full waveform inversion[J]. Geophysical Prospecting for Petroleum, 52 (4) : 417–425. DOI:10.3969/j.issn.1000-1441.2013.04.012
[31] Hu Guanghui, Wang Lixin, Fang Wubao, et al.2014. Full Waveform Inversion Method and Application[M]. Beijing: Petroleum Industry Press .
[32] Hu Guanghui, Wang Lixin, Wang Jie, et al.2015. Characteristics waveform inversion based on early arrival waves[J]. Geophysical Prospecting for Petroleum, 54 (1) : 71–76. DOI:10.3969/j.issn.1000-1441.2015.01.010
[33] Hu Ying, Zhang Dong, Yuan Jianzheng, et al.2015. An efficient multi-scale waveform inversion method in Laplace-Fourier domain[J]. Petroleum Exploration and Development, 42 (3) : 338–346. DOI:10.11698/PED.2015.03.10
[34] Jannane M, Beydoun W B, Crase E, et al.1989. Wavelengths of earth structures that can be resolved from seismic reflection data[J]. Geophysics, 54 (7) : 906–910. DOI:10.1190/1.1442719
[35] Jaiswal P, Zelt C A, Bally A W, et al.2008. 2-D traveltime and waveform inversion for improved seismic imaging: Naga Thrust and Fold Belt, India[J]. Geophysical Journal International, 173 (2) : 642–658. DOI:10.1111/gji.2008.173.issue-2
[36] Jaiswal P, Zelt C A, Dasgupta R, et al.2009. Seismic imaging of the Naga thrust using multiscale waveform inversion[J]. Geophysics, 74 (6) .
[37] Jervis M, Sen M K, Stoffa P L.1996. Prestack migration velocity estimation using nonlinear methods[J]. Geophysics, 61 (1) : 138–150. DOI:10.1190/1.1443934
[38] Jin S, Madariaga R.1994. Nonlinear velocity inversion by a two-step Monte Carlo method[J]. Geophysics, 59 (4) : 577–590. DOI:10.1190/1.1443618
[39] Lambaré G, Alerini M, Baina R, et al.2004. Stereotomography: A semi-automatic approach for velocity macromodel estimation[J]. Geophysical Prospecting, 52 (6) : 671–681. DOI:10.1111/gpr.2004.52.issue-6
[40] Lambaré G, Alerini M. 2005. Semi automatic PP-PS Stereotomography: Application to the synthetic Valhall dataset[C].// SEG Technical Program. Expanded Abstracts, 943-946.
[41] Lambaré G.2008. Stereotomography[J]. Geophysics, 73 (5) .
[42] Li Zhenchun.2002.The Study of Multi-gather Migration Velocity Modeling[D].Shanghai: School of ocean and earth science, Tongji university.
[43] Li Zhenchun.2011. Seismic prestack imaging theory and method[M]. Dongying: China university of petroleum press .
[44] Ling Y, Wang H Z, Liu S Y. 2014. Monte Carlo background velocity inversion method[C].//Beijing 2014 International Geophysical Conference & Exposition. Beijing, China, 735-738, doi: 10.1190/IGCBeijing2014-188.
[45] Liu F T, Qu K X, Wu H, et al.1989. Seismic tomography of the Chinese continent and adjacent region[J]. Acta Geophysica Sinica (in Chinese), 32 (3) : 281–291.
[46] Liu G F, Liu H, Meng X H, et al.2012. Frequency-related factors analysis in frequency domain waveform inversion[J]. Chinese Journal of Geophysics (in Chinese), 55 (4) : 1345–1353. DOI:10.6038/j.issn.0001-5733.2012.04.030
[47] Liu S W, Wang H Z, Cheng J B, et al.2008. Space-time-shift imaging condition and migration velocity analysis[J]. Chinese Journal of Geophysics (in Chinese), 51 (6) : 1883–1891. DOI:10.3321/j.issn:0001-5733.2008.06.031
[48] Liu Y Z, Wang G Y, Dong L G, et al.2014. Joint inversion of VTI parameters using nonlinear traveltime tomography[J]. Chinese Journal of Geophysics (in Chinese), 57 (10) : 3402–3410. DOI:10.6038/cjg20141026
[49] Liu Z Y.1997. An analytical approach to migration velocity analysis[J]. Geophysics, 62 (4) : 1238–1249. DOI:10.1190/1.1444225
[50] Ma Y, Hale D.2013. Wave-equation reflection traveltime inversion with dynamic warping and full-waveform inversion[J]. Geophysics, 78 (6) : R223–R233. DOI:10.1190/geo2013-0004.1
[51] Malinowski M, Operto S.2008. Quantitative imaging of the Permo-Mesozoic complex and its basement by frequency domain waveform tomography of wide-aperture seismic data from the Polish Basin[J]. Geophysical Prospecting, 56 (6) : 805–825. DOI:10.1111/gpr.2008.56.issue-6
[52] Ni Yao, Yang Kai, Chen Baoshu.2013. Stereotomography inversion method: theory and application testing[J]. Geophysical Prospecting for Petroleum, 52 (2) : 121–130. DOI:10.3969/j.issn.1000-1441.2013.02.002
[53] Operto S, Virieux J, Dessa J X, et al.2006. Crustal seismic imaging from multifold ocean bottom seismometer data by frequency domain full waveform tomography: Application to the eastern Nankai trough[J]. Journal of Geophysical Research: Solid Earth, 111 (B9) : B09306.
[54] Pratt R G. 1990. Inverse theory applied to multi-source cross-hole tomography. Part 1: Acoustic wave-equation method[J]. Geophysical Prospecting, 38(3): 287-310.
[55] Pratt R G, Goulty N R.1991. Combining wave-equation imaging with traveltime tomography to form high-resolution images from crosshole data[J]. Geophysics, 56 (2) : 208–224. DOI:10.1190/1.1443033
[56] Prieux V, Operto S, Brossier R, et al. 2009. Application of acoustic full waveform inversion to the synthetic Valhall velocity model[C].// SEG Technical Program. Expanded Abstracts, 2268-2272.
[57] Prieux V, Brossier R, Gholami Y, et al.2011. On the footprint of anisotropy on isotropic full waveform inversion: The Valhall case study[J]. Geophysical Journal International, 187 (3) : 1495–1515. DOI:10.1111/gji.2011.187.issue-3
[58] Prieux V, Lambaré G, Operto S, et al.2013. Building starting models for full waveform inversion from wide-aperture data by stereotomography[J]. Geophysical Prospecting, 61 (S1) : 109–137.
[59] Pyun S, Son W, Shin C.2011. 3D acoustic waveform inversion in the Laplace domain using an iterative solver[J]. Geophysical Prospecting, 59 (3) : 386–399. DOI:10.1111/gpr.2011.59.issue-3
[60] Qin Ning, Li Zhenchun, Yang Xiaodong, et al.2012. Image domain travel-time tomography velocity inversion based on automatic picking[J]. OGP, 47 (3) : 392–398.
[61] Qin N, Li Z C, Sang Y Y, et al.2014. The research of travel time tomographic velocity modeling method on first break[J]. Progress in Geophysics (in Chinese), 29 (1) : 255–260. DOI:10.6038/pg20140136
[62] Ravaut C, Operto S, Improta L, et al.2004. Multiscale imaging of complex structures from multifold wide-aperture seismic data by frequency-domain full-waveform tomography: application to a thrust belt[J]. Geophysical Journal International, 159 (3) : 1032–1056. DOI:10.1111/gji.2004.159.issue-3
[63] Ren Haoran.2012.A study of the seismic wave inversion method for imaging in the acoustic media[D].Shanghai: School of ocean and earth science, Tongji university.
[64] Sang Yunyun, Li Zhenchun, Zhang Kai.2013. Shortest path ray tracing based on parabolic travel-time interpolation[J]. OGP, 48 (3) : 403–409.
[65] Shin C, Ha W.2008. A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains[J]. Geophysics, 73 (5) : VE119–VE133. DOI:10.1190/1.2953978
[66] Shin C, Ha W, Kim Y.2013. Subsurface model estimation using Laplace-domain inversion methods[J]. The Leading Edge, 32 (9) : 1094–1099. DOI:10.1190/tle32091094.1
[67] Shin C, Ha W, Jun H, et al.2014. 3D Laplace-domain full waveform inversion using a single GPU card[J]. Computers & Geosciences, 67 : 1–13.
[68] Sirgue L. 2006. The importance of low frequency and large offset in waveform inversion[C].//68th Annual International Meeting. London: EAGE, A037.
[69] Sirgue L, Barkved O I, Dellinger J, et al.2010. Full waveform inversion: The next leap forward in imaging at Valhall[J]. First Break, 28 (4) : 65–70.
[70] Song Z M, Williamson P R, Pratt R G.1995. Frequency-domain acoustic-wave modeling and inversion of crosshole data; Part II, inversion method, synthetic experiments and real-data results[J]. Geophysics, 60 (3) : 796–809. DOI:10.1190/1.1443818
[71] Stork C.1992. Reflection tomography in the postmigrated domain[J]. Geophysics, 57 (5) : 680–692. DOI:10.1190/1.1443282
[72] Symes W W, Carazzone J J.1991. Velocity inversion by differential semblance optimization[J]. Geophysics, 56 (5) : 654–663. DOI:10.1190/1.1443082
[73] Symes W W.2008. Migration velocity analysis and waveform inversion[J]. Geophysical Prospecting, 56 (6) : 765–790. DOI:10.1111/gpr.2008.56.issue-6
[74] Taillandier C, Noble M, Chauris H, et al.2009. First-arrival traveltime tomography based on the adjoint-state method[J]. Geophysics, 74 (6) : WCB1–WCB10. DOI:10.1190/1.3250266
[75] Tarantola A.1984. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 49 (8) : 1259–1266. DOI:10.1190/1.1441754
[76] Virieux J, Operto S.2009. An overview of full-waveform inversion in exploration geophysics[J]. Geophysics, 74 (6) : WCC1–WCC26. DOI:10.1190/1.3238367
[77] Wang Huazhong, Wang Xiongwen, Wang Xiwen.2012. Analysis of the basic problems of seismic wave inversion[J]. Lithologic Reservoirs, 24 (6) : 1–9.
[78] Wang X W, Yang W Y, Lu B, et al.2013. Grasp of the world geophysical technology development status, to promote seismic data processing and interpretation of technological progress[J]. Progress in Geophysics (in Chinese), 28 (1) : 224–239. DOI:10.6038/pg20130124
[79] Woodward M J, Nichols D, Zdraveva O, et al.2008. A decade of tomography[J]. Geophysics, 73 (5) : VE5–VE11. DOI:10.1190/1.2969907
[80] Xu S, Wang D L, Chen F, et al. 2012. Inversion on reflected seismic wave[C].//SEG Annual Meeting. Las Vegas, Nevada: SEG, 1-7.
[81] Yang K, Wang H Z, Ma Z T.2005. An output imaging scheme of the common reflection surface stack[J]. Journal of Seismic Exploration, 14 (3) : 131–150.
[82] Yang Qinyong, Hu Guanghui, Wang Lixin.2014. Research status and development trend of full waveform inversion[J]. Geophysical Prospecting for Petroleum, 53 (1) : 77–83. DOI:10.3969/j.issn.1000-1441.2014.01.011
[83] Yang W Y, Wang X W, Yong X S, et al.2013. The review of seismic Full waveform inversion method[J]. Progress in Geophysics (in Chinese), 28 (2) : 766–776. DOI:10.6038/pg20130225
[84] Zelt C A, Barton P J.1998. Three-dimensional seismic refraction tomography: A comparison of two methods applied to data from the Faeroe Basin[J]. Journal of Geophysical Research, 103 (B4) : 7187–7210. DOI:10.1029/97JB03536
[85] Zelt C A, Pratt G, Brenders A, et al. 2005. Advancements in long-offset seismic imaging: A blind test of Traveltime and waveform tomography[C].//American Geophysical Union, Spring Meeting. New Orleans: American Geophysical Union.
[86] 曹书红, 陈景波.2014. 频率域全波形反演中关于复频率的研究[J]. 地球物理学报, 57 (7) : 2302–2313. DOI:10.6038/cjg20140724
[87] 陈志德, 刘桭宽, 李成斌.2002. 三维叠前深度偏移速度分析及蒙特卡洛自动层速度拾取[J]. 地球物理学报, 45 (2) : 246–254. DOI:10.3321/j.issn:0001-5733.2002.02.011
[88] 成谷, 张宝金, 马在田.2007. 浅谈反射地震走时层析中的正则化[J]. 地球物理学进展, 22 (6) : 1715–1721. DOI:10.3969/j.issn.1004-2903.2007.06.005
[89] 董良国, 迟本鑫, 陶纪霞, 等.2013. 声波全波形反演目标函数性态[J]. 地球物理学报, 56 (10) : 3445–3460. DOI:10.6038/cjg20131020
[90] 胡光辉, 贾春梅, 夏洪瑞, 等.2013. 三维声波全波形反演的实现与验证[J]. 石油物探, 52 (4) : 417–425.
[91] 胡光辉, 王立歆, 方伍宝, 等.2014. 全波形反演方法及应用[M]. 北京: 石油工业出版社 .
[92] 胡光辉, 王立歆, 王杰, 等.2015. 基于早至波的特征波波形反演建模方法[J]. 石油物探, 54 (1) : 71–76.
[93] 胡英, 张东, 袁建征, 等.2015. Laplace-Fourier域多尺度高效全波形反演方法[J]. 石油勘探与开发, 42 (3) : 338–346.
[94] 李振春. 2002. 多道集偏移速度建模方法研究[博士论文]. 上海: 同济大学.
[95] 李振春.2011. 地震叠前成像理论与方法[M]. 东营: 中国石油大学出版社 .
[96] 刘福田, 曲克信, 吴华, 等.1989. 中国大陆及其邻近地区的地震层析成象[J]. 地球物理学报, 32 (3) : 281–291.
[97] 刘国峰, 刘洪, 孟小红, 等.2012. 频率域波形反演中与频率相关的影响因素分析[J]. 地球物理学报, 55 (4) : 1345–1353. DOI:10.6038/j.issn.0001-5733.2012.04.030
[98] 刘守伟, 王华忠, 程玖兵, 等.2008. 时空移动成像条件及偏移速度分析[J]. 地球物理学报, 51 (6) : 1883–1891. DOI:10.3321/j.issn:0001-5733.2008.06.031
[99] 刘玉柱, 王光银, 董良国, 等.2014. VTI介质多参数联合走时层析成像方法[J]. 地球物理学报, 57 (10) : 3402–3410. DOI:10.6038/cjg20141026
[100] 倪瑶, 杨锴, 陈宝书.2013. 立体层析反演方法理论分析与应用测试[J]. 石油物探, 52 (2) : 121–130.
[101] 秦宁, 李振春, 杨晓东, 等.2012. 自动拾取的成像空间域走时层析速度反演[J]. 石油地球物理勘探, 47 (3) : 392–398.
[102] 秦宁, 李振春, 桑运云, 等.2014. 初至波走时层析速度建模方法研究[J]. 地球物理学进展, 29 (1) : 255–260. DOI:10.6038/pg20140136
[103] 任浩然. 2012. 声介质地震波反演成像方法研究[博士论文]. 上海: 同济大学.
[104] 桑运云, 李振春, 张凯.2013. 抛物旅行时插值最短路径射线追踪[J]. 石油地球物理勘探, 48 (3) : 403–409.
[105] 王华忠, 王雄文, 王西文.2012. 地震波反演的基本问题分析[J]. 岩性油气藏, 24 (6) : 1–9.
[106] 王西文, 杨午阳, 吕彬, 等.2013. 把握世界物探技术发展现状, 促进地震资料处理、解释技术进步[J]. 地球物理学进展, 28 (1) : 224–239. DOI:10.6038/pg20130124
[107] 杨勤勇, 胡光辉, 王立歆.2014. 全波形反演研究现状及发展趋势[J]. 石油物探, 53 (1) : 77–83.
[108] 杨午阳, 王西文, 雍学善, 等.2013. 地震全波形反演方法研究综述[J]. 地球物理学进展, 28 (2) : 766–776. DOI:10.6038/pg20130225