全波形反演(Full Waveform Inversion，简称FWI, Tarantola，1984；Pratt，1999)利用全波形信息，对地震炮记录与模型正演炮记录的差值建立目标泛函，利用最优化思想计算梯度方向，通过不断迭代更新地下介质参数使两者残差越来越小，从而实现优化介质参数的目的.全波形反演不仅考虑了地震波传播的走时等运动学信息，而且加入了振幅、相位等动力学信息，能够适应强横向变速介质和各向异性介质的速度反演.它是目前精度最高的一种速度反演方法，通过迭代反演得到高精度的地下构造，能够满足目前勘探开发日益复杂的需求，可以为叠前成像技术提供更准确的速度模型，最终提高油气勘探的成功率.

全波形反演方法需要大量的模拟地震数据，其计算量正比于炮点和检波点的数量，对一个大型的地震数据进行波形反演的计算量很大.为了提高计算效率，Trezeguet等(1989)将偏移和波形反演结合，迭代中采用成像条件和波形拟合双重标准，提高速度估计的效率.Gong等(2008)提出了基于逆时偏的波形反演方法，减少了炮点和检波点数量，提高了计算效率.Krebs等(2009)利用相位编码技术对所有的炮进行编码，减少了反演过程中计算量.Wang和Goo(2010)提出了多炮地震数据的快速全波形反演方法，这种方法通过生成几个超级炮点，使得炮的数量减少，在不影响反演质量的情况下提高计算效率.Mao和Wu(2012)在波形反演中引入了GPU加速方法，在很大程度上提高了FWI的计算效率.

FWI是一个高度非线性欠定问题(董良国等，2013；杨午阳等，2013)，观测数据和模型数据之间是非线性对应关系，如果初始速度模型和实际的速度模型相差太多，就会出现局部极值以及周期移位现象，造成迭代反演结果不收敛.针对反演中存在的这些问题，Bunks等(1995)采用空间网格和频率多尺度方法进行波形反演，在反演过程中从低频开始不断地增加高频成分，并对应对模型网格进行加密，很好地解决了各个尺度上解的收敛问题.Luo和Schuster(1990)利用走时互相关作为目标函数进行反演，减弱了对初始模型的依赖；Shin和Cha(2009)将FWI拓展到拉普拉斯域来反演速度的低频信息；Xu等(2012)提出反射波波形反演方法，利用反射波与低波数的速度结构之间的线性关系反演速度场中的背景速度；Chi等(2013)利用基于能量包络的目标函数得到观测数据中的超低频信息进行波形反演，消除了局部极小值的影响，克服了因缺失低频成分无法迭代收敛的缺陷.利用走时层析或者偏移速度分析也可以给FWI提供比较准确的初始速度模型，从而提高反演解的迭代收敛速度(De Hoop, 1960; Shen and Symes, 2008).Weibull等(2012)把目标函数改成联合偏移速度和全波形反演的形式，解决了初始模型差引起的周期移位问题.陈生昌和陈国新(2016)提出利用基于时间二阶积分的波场反演背景速度.为了提高反演的整体收敛效果，有些研究者还提出了对地震数据加窗的反演思路(Yoon et al., 2003, 2012)，还有学者提出了利用波场能量加权梯度的方法(Zhang et al., 2011, 2012).另外加入一些额外的先验信息也可以约束解收敛到全局最优，Guitton等(2010)通过得到的地质倾角信息作为约束条件，Asnaashari等(2012)通过已知的井数据插值成先验速度信息作为反演过程中的约束，得到了很好的迭代反演结果.

van Leeuwen和Herrmann(2013)针对计算效率和局部极值问题，提出了一种改进的波形反演方法：波场重构反演(Wavefield Reconstruction Inversion，简称WRI).该方法通过将波动方程作为惩罚项加入到目标函数中，同时将波场和模型参数作为反演参数进行反演.该方法不需要存储正演波场和计算伴随波场，因此计算效率得到了很大提高，通过拓宽解的寻找空间，WRI的解更加稳定，受局部极值的影响较弱.波场重构反演方法通过惩罚因子将波形匹配项与波动方程惩罚项连接在一起构成目标函数，该常数用来增强或减弱波动方程对目标函数的约束，对反演的精度和稳定性至关重要.van Leeuwen等(2014)优化了WRI的相关理论，给出了重构波场的物理含义并对目标函数中惩罚因子的取值进行了经验性分析；Fang等(2015)在Bayesian理论下对WRI的不确定性进行了定量分析.van Leeuwen和Herrmann(2016)给出了惩罚因子的选取策略：通过使波形匹配项和波动方程项达到同一数量级来求取惩罚因子初始值，然后在反演过程中逐步增大.该策略与da Silva(2015)提出的利用积性函数求取惩罚因子的方法类似，但精度较低，难以保证反演的稳定性.

针对波场重构反演中惩罚因子精度较低等问题，本文将波场重构反演拓展到时间域，并利用梯度法进行波场重构.根据时间域增广波动方程，提出了一种新的基于约束优化理论的惩罚因子算法.根据本文提出的算法，针对波形反演在应用时可能存在的噪音干扰、子波错误和低频信息缺失等问题，利用部分Sigsbee2A模型合成数据分别进行了测试，并取得了较好的反演结果.

全波形反演方法利用波场的振幅、走时以及相位信息，根据波形匹配准则通过求取梯度来得到地下介质参数.传统的FWI目标函数表达式如下：

(1)

上式中m为速度，d为观测数据；P是炮点算子，以此来确定炮点对应的模拟数据.A为正演算子，具体表达式：，t为时间，Δ是离散的拉普拉斯算子.

FWI的目标函数是L₂范数下的泛函，可以用伴随状态法(Plessix, 2006)求解:

(2)

求解上式需要计算并存储所有炮点的正演波场：A(m)u=q和伴随波场ψ：A(m)^Tψ=P^T(Pu－d)，对内存的要求很高，而且计算效率较低；在没有准确的初始模型的情况下反演结果很容易受局部极小值影响.

为解决上述问题，van Leeuwen和Herrmann(2013)在频率域通过修改目标函数，提出一种基于波场重构的反演方法：将状态方程加入到目标函数中，波场和模型都作为参数进行反演.WRI目标函数的一般数学表示为

(3)

上式第一项为传统全波形反演中的模拟数据与观测数据的匹配项；第二项为惩罚项，即波动方程，λ > 0为惩罚因子.

WRI的目标函数包含两个未知量，无法直接求解，但在给定模型参数m的情况下，可以得到一个扩展的状态变量，即重构波场：

(4)

其中T代表矩阵的转置，u_λ为重构波场，该式又称为增广波动方程.可以看出重构波场与惩罚因子的取值有直接关系，而我们无法根据上式得到一个准确的惩罚因子算法.

下面将WRI拓展到时间域.首先将WRI的目标函数重写为如下形式：

(5)

在最小二乘的理论框架下，对目标函数关于波场求导，整理可得重构波场的更新表达式：

(6a)

(6b)

其中U为重构波场，g_u为波场更新梯度；但该式由于难以判断收敛性而无法进行求解，因此将重构波场的表达式重写为

(7)

可以看出，重构波场是观测波场与模型波场通过平衡因子的组合波场，该波场满足已知波动方程，同时包含一定的观测数据信息.新的惩罚因子λ′的取值需要根据初始波场以及观测数据的准确性综合给出，具体算法见下节.

在重构波场之后，波场与模型参数不再耦合，可以直接求取模型更新的梯度：

(8)

WRI引入惩罚因子的目的是用来调节对远离波动方程约束项的惩罚，但在时间域中，惩罚因子转化为权重项来组合模拟波场与观测波场，以期得到真实的波场或者弥补两者各自的不足从而得到更加准确的反演结果.惩罚因子的准确性直接影响到反演结果的精度，为此我们将其求解过程总结如下等式约束优化问题：

(9)

其中u_t为真实波场.基于以下原因：1)波动方程准确；2)重构波场与模型参数不再耦合；3)观测数据品质的缺陷可由正演算子进行修正；我们拟采用波动方程作为目标函数：

(10)

因此可以得到下列关于惩罚因子的方程：

(11)

同样采用梯度法计算惩罚因子的更新梯度：

(12)

在本文中取中间值λ=0.5作为初始值进行更新.注意到惩罚因子的求取是在当前时刻的初始波场和观测波场下计算的，是关于一个常数的更新，计算量较小，易于在反演中实现.

惩罚因子的计算过程本质上是在波动方程的约束下，对初始波场和观测波场互相匹配优化的过程.当炮记录包含噪音或者缺少低频信息，或者地震子波不准确时，WRI可以通过准确的惩罚因子来过滤干扰信息从而准确地重构波场进行反演.

根据上节提出的惩罚因子算法，首先给出声波近似下一个模型试算结果来验证该方法的正确性，然后在该算法的基础上分别从噪音干扰，子波错误以及低频信息缺失三个问题分析WRI的适应性.

这里本文以部分Sigsbee2A模型为例进行试算，检验WRI以及惩罚因子算法的正确性.图 1为真实模型，模型大小为横向7000 m，纵向2400 m，纵横向采样点分别为700、300，网格大小横向10 m，纵向8 m；地表放炮，地面接收，起始炮点在500 m处，共30炮，炮间距200 m；震源子波是主频为10 Hz的Ricker子波，最大接收时间为2.5 ms，时间采样率为1 ms；图 2为初始平滑模型，以下测试基本参数与本节相同；图 3为迭代10次后的WRI反演结果.

图 1

Fig. 1

图 1 真实模型 Fig. 1 True model

图 2

Fig. 2

图 2 初始模型 Fig. 2 Initial model

图 3

Fig. 3

图 3 10次迭代后反演结果 Fig. 3 Inversion result after 10 iterations

由图 3可以看出，在10次迭代之后，WRI可以给出相对准确的反演结果.分别抽取3000 m(图 4a)和4000 m(图 4b)处的反演速度与初始速度和真实速度进行比较可以发现，针对浅层位置的速度，WRI可以得到准确的速度；但由于来自深层的反射波较弱，难以准确重构真实波场，所以反演结果精度相对较低，但基本可以得到深层的基础构造；针对深层异常体，在迭代次数较少时只能得到其准确位置.

图 4

Fig. 4

图 4 (a) 3000 m处速度曲线对比；(b) 4000 m处速度曲线对比 Fig. 4 (a) Comparison of velocity at 3000 m; (b) Comparison of velocity at 4000 m

在实际生产中，地震记录中会存在大量噪音，反演方法对噪音的敏感性至关重要.下面我们将炮记录加入一定的噪音进行反演，测试该方法对噪音的适应性.

本节所用噪音类型为高斯噪音，信噪比0.41.图 5a为普通正演炮记录，图 5b为包含噪音干扰的炮记录.图 6为加入噪音后FWI与WRI的反演结果.由于噪音能量较强，所以在反演时对表层能量进行切除，对表层以下速度进行更新.传统的FWI反演方法直接通过观测数据构造伴随波场，当观测数据包含较多噪音时，其反演结果会受其影响而导致速度不准确，反演结果包含很多假象，难以辨别各处构造(图 6a).而WRI在进行模型更新之前首先计算惩罚因子，利用优化的惩罚因子重构波场，可以减弱噪音的影响，使反演结果更加准确(图 6b).抽取3000 m(图 7a)和4000 m(图 7b)处的真实速度与FWI反演结果和WRI反演结果进行比较：FWI的速度波动较大，偏离真实速度，而WRI反演结果较接近真实速度.

图 5

Fig. 5

图 5 (a) 模拟炮记录；(b)加噪音后炮记录 Fig. 5 (a) Seismograms; (b) Seismograms with noise interference

图 6

Fig. 6

图 6 (a) 加噪音后FWI反演结果；(b)加噪音后WRI反演结果 Fig. 6 (a) FWI inversion result with noise; (b) WRI inversion result with noise

图 7

Fig. 7

图 7 (a) 3000 m处速度曲线比较；(b) 4000 m处速度曲线比较 Fig. 7 (a) Comparison of velocity at 3000 m; (b) Comparison of velocity at 4000 m

在实际生产中，激发地震的子波是未知的，一般需要首先进行子波评定才可以进行反演.而FWI是基于波形匹配准则进行反演的，当地震子波不准确时，FWI易产生周期移位现象，导致反演结果不准确.

频率域的WRI目标函数和增广波动方程对子波的依赖性较强，子波的不准确直接导致重构波场以及模型梯度的错误(Fang et al., 2015).但通过将WRI拓展到时间域可以发现，波场重构过程以及模型梯度是利用波场进行计算的，不涉及波形匹配，对地震子波要求较低.而且基于新的惩罚因子算法，可以尽可能准确地重构真实波场，减弱观测波场和模拟波场子波不匹配的影响.下面我们利用不同的子波进行正演和反演，测试基于新的惩罚因子的WRI对地震子波的敏感性.

图 8中蓝色为正演使用的地震子波，红色为反演时采用的子波.图 9a为子波错误时WRI的反演结果，图 9b为对应FWI反演结果.抽取3000 m和4000 m处速度进行比较(图 10)，可以看出当使用错误的子波进行反演时，由于波形不匹配，FWI在1000 m处和1500 m处存在明显的周期移位现象，反演结果偏离真实速度，但WRI通过重构波场计算模型更新梯度，对子波的依赖性较小，部分地方(1700 m和2000 m)的速度由于波场的不准确会有所波动，但整体速度基本准确，接近真实速度.

图 8

Fig. 8

图 8 不同地震子波 Fig. 8 Different seismic wavelet

图 9

Fig. 9

图 9 (a) WRI反演结果；(b) FWI反演结果 Fig. 9 (a) WRI inversion result; (b) FWI inversion result

图 10

Fig. 10

图 10 (a) 3000 m处速度曲线比较；(b) 4000 m处速度曲线比较 Fig. 10 (a) Comparison of velocity at 3000 m; (b) Comparison of velocity at 4000 m

在实际地震勘探中，低频信息是缺失的.而没有充足的低频信息，FWI易陷入局部极值，导致无法收敛.但WRI同时从模型空间和数据空间寻求最优解，对低频信息的依赖性较小，可以得到一个比较好的反演结果.本节通过对地震子波进行滤波处理来检测WRI对低频信息的依赖性.

图 11a中蓝色为时间域正常地震子波，红色为滤除5 Hz以下频率后的子波，图 11b为频率域示意图.图 12为不同子波的炮记录频谱，可以发现在对子波进行滤波之后，炮记录中几乎不包含任何低频信息，比较符合实际地震勘探情况.为测试本文算法对低频信息的敏感性，本节采用梯度模型(图 13)作为初始速度场.图 14a为FWI针对缺少低频信息的炮记录反演得到的速度场，图 14b为WRI反演结果.两者在浅层的速度准确，但针对深层的异常体，FWI反演结果并未收敛，而且由于缺失低频信息，FWI在深层的速度偏高，远离真实速度，而且深层的异常体并未收敛.WRI通过准确重构波场可以补充部分低频信息，而且同时从模型空间和波场空间进行求解，受局部极值影响较弱，因此反演结果准确.同样抽取3000 m(图 15a)和4000 m(图 15b)处的速度进行比较，可以看出在缺少低频信息的情况下，基于同样的迭代次数，FWI只能反演处出浅层的速度，在1000 m之后反演结果远偏离真实速度，而WRI的反演结果从浅层至深层一直相对比较准确，深层的异常体已经收敛，而且位置比较准确.

图 11

Fig. 11

图 11 (a) 时间域滤波前后子波；(b)频率域滤波前后子波 Fig. 11 (a) Wavelet before and after filter in time domain; (b) Wavelet before and after filter in frequency domain

图 12

Fig. 12

图 12 滤波前后炮记录频谱 Fig. 12 Spectrum of seismograms before and after filter

图 13

Fig. 13

图 13 初始速度场 Fig. 13 Initial Velocity

图 14

Fig. 14

图 14 (a) 无低频信息下的FWI反演结果；(b)无低频信息下的WRI反演结果 Fig. 14 (a) FWI inversion result without low frequency; (b) WRI inversion result without low frequency

图 15

Fig. 15

图 15 (a) 3000 m处速度曲线比较；(b) 4000 m处速度曲线比较 Fig. 15 (a) Comparison of velocity at 3000 m; (b) Comparison of velocity at 4000 m

本文在时间域对波场重构反演方法进行了研究，推导了时间域增广波动方程和模型梯度，提出了基于约束优化的惩罚因子算法.基于准确的惩罚因子，并结合波场重构反演理论的优点，从噪音、地震子波错误和低频信息缺失等三个方面对该理论的适用性进行了测试.数值实验表明波场重构反演具体很高的精度，而且通过对惩罚因子的准确计算，可以实现对模拟波场和观测波场的准确组合，在无法准确预测地震波传播特点的情况下(子波错误，噪音等)，仍然可以给出较好的反演结果.针对地震勘探中低频信息不足的问题，WRI通过同时从模型空间和数据空间寻求最优解，能够有效地避免局部极值，给出准确的反演速度.