地震波在地下介质中传播的速度是地震勘探领域最重要的参数之一，它对于高精度的地震成像以及岩性与储层刻画都是必不可少的.传统的偏移速度分析方法(Biondi et al.,2004)和走时层析反演方法(Schuster et al.,1993)都只能得到宏观速度场，即速度的低频成分.全波形反演(FWI)以地震数据的波形信息为依据，直接基于波动方程，可以反演得到地下速度场的中高频信息(Sheng et al.,2006)， 能够满足复杂构造成像方法对速度参数的精度要求，具备岩性反演与储层预测的潜力，应用前景非常广泛.

“全波形反演”是Pan在1988年明确提出的(Pan et al.,1988)，与早期文献中提到的“波形反演”(Tarantola,1984; Gauthier et al.,1986; Mora,1986)本质上是一样的，它起源于Laily(1983)和Tarantola(1984)，他们将地震成像表述为一个局部最优化反演问题——寻求模拟数据与实际地震记录之间最小二乘意义下的最佳拟合，且给出了以震源正传波场与数据残差反传波场相关的方式来计算反演梯度的方法，避免了Fréchet导数的直接计算，掀起了全波形反演的研究热潮(Mora,1988; Pratt et al.,1990; Pratt et al.,1996; Operto et al.,2004).全波形反演的每次迭代都需要若干次地震波正演计算，且反演过程是非线性的，因此提高反演的计算效率和避免反演落入局部极小值是全波形反演研究的两大难题.

地震波正演是全波形反演的计算引擎，它是决定全波形反演规模与效率的最核心因素.地震波正演可以在频率域(Liao et al.,1996; Pratt et al.,1998; 刘国峰等，2012)或时间域(Bunks et al.,1995; Zhou et al.,1995;丁继才等，2007)进行.Pratt等(1990)将频率域正演应用到了全波形反演，二维(2D)频率域全波形反演常使用LU分解来模拟波场，在多炮模拟时具有优势(Operto et al.,2004).LU分解需要存储分解矩阵，对内存的需求量大，不能适应大模型尤其是三维(3D)模型的计算问题，如Operto等(2007)给出的三维黏声波正演方法，对于边长10 km的3D立方体模型，在网格间隔50 m及频率10 Hz时，需要约320GB的内存来存储单精度的LU分解结果.这迫使3D频率域正演要采用迭代法(Warner et al.,2008)，而迭代法的效率依赖预条件选取的好坏，预条件选的不好，迭代次数会较多，导致正演效率很低，另外迭代法正演需要对每炮分别计算，这使得频域全波形反演失去了多炮情况下的优势.

时间域正演方法通过在时间方向上的迭代来求解各个时间步的波场值，内存需求相对较少，能够适应3D问题.目前用于全波形反演的时间域正演方法多是基于规则网格(Bunks et al.,1995;丁继才等，2007)，规则网格的大小取决于模型最小速度参数与波场最大频率成分的比值，而时间步长又取决于网格大小与最大速度参数的比值.依据介质速度参数调整网格尺寸，则能显著减少网格单元数和加大时间步长，进而减少计算量和内存需求量.复杂介质条件下，速度沿深度方向和横向上的变化很大，而非规则网格局部单元的大小能与介质局部速度值相匹配，速度大的地方，网格大而疏，速度小的地方，网格小而密；而规则网格的网格单元大小只与最小速度值相匹配，其在模型全区域都是小而密的.显然，采用规则网格的地震波正演需要更多的内存量进行波场值的计算和存储，计算量也会增大，这会给大模型的全波形反演造成困难.另外，全波形反演的每次迭代都会对速度场进行更新，非规则网格能通过自身的调整更好地适应速度的变化，以避免出现数值频散影响反演结果.目前在地震勘探领域，将非规则网格地震波正演应用到时间域全波形反演的方法还未见报道.

非线性优化算法也是影响全波形反演效率的一个重要因素，好的优化算法能带来更快的收敛.全波形反演在数学上一般被抽象为一个二次型优化问题，对于二次型优化问题牛顿法的收敛速度最快，但是牛顿法需要计算Hessian矩阵，即目标函数的二阶导数，该矩阵维度是梯度维度的平方，因此计算和存储开销太大(Pratt et al.,1998).最速下降法直接用目标函数梯度的反方向作为搜索方向，无需其他处理，简单易实施，但其收敛速度较慢.寻求收敛速度快而又不显著增加计算量的优化算法一直是全波形反演研究的一个热点，在过去几十年里，多数学者采用了共轭梯度法(Mora,1986;Crase et al.,1990)，该方法采用共轭方向代替梯度反方向作为反演的搜索方向，每次迭代所需的共轭方向由当前梯度方向与上次迭代的共轭方向组合产生，理论上共轭梯度法是n步收敛的，n等于梯度的维度，可见其收敛效果仍不理想.将近似Hessian矩阵用于优化算法中可以得到接近牛顿法的收敛效果，拟牛顿法即为这类方法，实际应用中为减少存储量和计算量多采用L-BFGS算法(Nocedal,1980)，Brossier等对比了共轭梯度法和L-BFGS法，显示了L-BFGS在全波形反演中的优越性(Brossier et al.,2009).

全波形反演的局部极小值问题是非线性反演所 固有的难题，当初始模型足够接近真实模型(Gauthier et al., 1986)，或有长偏移距数据和透射波数据来重建反演模型中的低频和中频信息时(Mora,1986,1988; Pratt et al.,1990; Pratt et al.,1996)，可使反演正确收敛到全局极小值.学者们还提出了一些提高反演线性度的方法来最大程度地避免全波形反演落入局部极小值：Sheng等(2006)提出了加时间窗选取早至数据以提高反演的线性度；Bunks等(1995)提出了由低频到高频来重建反演模型的时间域多尺度反演策略，低频数据的全波形反演要比高频时更线性化； Sirgue和Pratt(2004)将Bunks的时域多尺度反演策略推广到了频率域，且从波数覆盖的角度对频率域全波形反演最优频率的选择进行了研究.加窗及多尺度反演本质上都是通过改变目标函数的特性，以拓宽其超曲面“最深谷”的范围来保证收敛性(张剑锋，1993).

为了满足全波形反演在内存和计算效率上的需求，本文采用了基于非规则三角网格的地震波正演方法，应用随速度参数自适应最优化剖分技术来离散模型空间，该方法能根据迭代更新后的模型参数自适应地调整和优化网格，以实现反演分辨率与模型参数的最佳匹配，且能避免反演中可能出现的频散问题；结合Wiener滤波给出了时间域全波形反演最优频带选取方法，采用多尺度反演策略来避免反演落入局部极小值，应用L-BFGS优化算法来获取更快的收敛.为验证本文这套时间域全波形反演方法的有效性，文中给出了Overthrust数值算例.

常见的基于非规则网格的地震波正演方法有谱元法和有限元法，它们都要利用由变分法得到的刚度矩阵，存储这个矩阵需要内存较多，若每次重新计算则计算量又很大.本文使用的地震波正演方法是格子法(Zhang et al.,1999; 徐义等，2008)，该方法综合了有限元与有限差分法的优点，一般采用非规则、非结构化三角网格，便于精细刻画复杂构造，内存需求小，且计算效率与差分法相当.

2.1 正演方法

本小节简略介绍一下格子法的基本原理，模型边界处采用非规则网格PML方法来处理.为便于处理介质的非均匀性，使用如下声波方程(徐义等，2008)：

方程中p(x ,t)为波场值，ρ(x)和v(x)分别为地下介质密度和速度，s(x _s,t)为震源函数.

格子法的核心是积分平衡的微分方程弱形式，如图 1所示，在k点邻域对式(1)作面积分并应用格林定理可得到

式中m是绕节点k的三角单元数，S_kl是节点k周围第l个三角形单元中的虚线段，α与β是虚线段包络线的外法线方向余弦,是震源函数的面积分.利用动力学计算中的集中质量模型及三角单元线性插值方法，可得到式(2)的空间离散形式：

其中式中下标r指代三角单元中的 节点，A是三角单元面积，b_k和a_k是三角单元的几 何参数，以图 1中三角单元ijk为例，b_k=(z_i－z_j)/2,a_k=(x_i－x_j)/2. 用二阶中心差分来离散(3)式左端项中的时间导数，即可实现波场值的迭代更新.

图1 非规则三角网格局部示意图 Fig.1 Local irregular triangle-grid

非规则网格PML衰减带内的计算引入了随边界点位置变化而改变的局部坐标系，建立了基于局部坐标系的PML分裂方程， 类似地，可得到其积分平衡弱形式，继而推导得到其离散公式，具体可参考相关文献(徐义等，2008)，在此不再赘述.

大规模的非线性反演一般采用局部求导类优化算法加线性搜索来实现，待反演的参数太多，全局搜索方法效率太低，在目前的计算条件下很难应用.局部求导类优化算法的收敛性和线性搜索的效率是我们选择时考虑的首要因素，下面介绍一下本文所采用的L-BFGS优化算法和步长搜索策略.

3.1 数值优化算法

BFGS方法是以其发明人Broyden、Fletcher、Goldfarb和Shanno命名的，它是拟牛顿法的一种，具有快速的收敛性，在实际应用时多采用L-BFGS 算法.BFGS的核心是利用梯度来构建近似逆Hessian 阵，其公式表示如下(Nocedal,1980)：

改写为乘积形式：

其中 分别是第k次反演迭代时的模型参数、反演梯度和近似逆Hessian阵， H ₀选择为一正定对角阵即可.当只保留有限项更新时即为L-BFGS方法，如保留m项，则当k+1>m时逆Hessian矩阵有如下更新形式：

3.2 步长搜索

步长搜索是个一维搜索问题，即沿着优化搜索方向搜索迭代步长，目的是加速收敛.数学上有黄金分割法、斐波那契序列法等，这些方法适用于目标函数值计算量较小的数学问题，对于全波形反演而言，其目标函数的每一次计算都需要一次波场正演，因此寻求如何用较少的搜索次数找到合适的步长非常重要.本文采用了抛物线拟合的方式来估计步长(Vigh et al，2009),以减少搜索次数，如图 2所示.该方法的思想是：找到两个步长α₁和α₂，使其 对应的目标函数满足条件S(α₁)0)和S(α₁)2)， 其中α₀=0，由三个步长所对应的点拟合一条抛物线，抛物线极小点处对应的步长α_opt即为最终估计步长.

图2 抛物线拟合法步长搜索 Fig.2 Step search by parabolic fitting

多尺度全波形反演最早由Bunks等(1995)在时间域实现，他们利用低通滤波器和多网格法由低频到高频逐步反演，以避免反演落入局部极小值.与频域多尺度全波形反演相比较，时域方法利用了数据中的冗余频率信息，具有更好的抗噪性(Sirgue et al.,2004),而实际数据总是含噪的，故时域方法具有更好的实际应用价值；参考Sirgue等(2004)的频率优选方案，时间域多尺度全波形反演 也可以通过类似的方式选择最优频带(Boonyasiriwat et al.,2009)，通常两到三个频带即可，这极大提高了反演效率.

4.3 变尺度三角网格

目前多尺度全波形反演一般采用矩形网格来 离散计算区域(Mora, 1988; Bunks et al.,1995; Boonyasiriwat et al.,2009)，2D地震波场模拟时网 格步长和时间步长要满足稳定性条件min(dx,dz) >a·dt·v_max和频散条件max(dx,dz)min/(bf_max)， dx和dz分别是离散网格沿水平和垂直方向的采样步长(允许dx=dz)，dt是数值计算时间步长，v_min,v_max分别是介质速度的最小值和最大值，a,b是常数(不同的计算方法其取值不同).

如果能让网格空间步长与网格单元的局部速度关联起来，而不是和最小速度相关联，则能大幅较少网格单元数，且能采用更大的时间步长，从而显著减少正演内存需求和计算量.为此，我们采用了与速度参数相关的变尺度三角网格来离散模型空间，相应地需要满足条件：dt<min[h_local/(av_local)]，h_locallocal/(bf_max)，h_local代表局部三角单元的边长尺度，v_local是相应的局部介质速度.

将与局部速度参数相关的网格剖分方法应用到全波形反演上还有两个显著优点：一是实现了反演分辨率在空间上的自适应变化，速度小的地方地震波长小，分辨率高，需要细粒度的网格，反之需要粗粒度的网格；二是能更好地克服反演过程中出现的频散问题，我们的网格尺度是与局部速度相匹配的，对反演中出现的速度更新适应性更好，当局部速度变化过大时能调整网格尺度重新适应新速度.图 3是Overthrust模型，图 4是与该模型相适应的非结构化三角网格，该图是取f_max=1.5 Hz时剖分的示意图，实际计算中的网格要比这个加密很多，由此图我们可以直观地看到网格尺度随速度的变化.

图3 二维Overthrust速度模型 Fig.3 2D Overthrust velocity model

图4 基于二维Overthrust速度模型的非规则三角网格 Fig.4 Irregular-grid based on 2D Overthrust velocity model

二维Overthrust模型含有一个大的逆冲推覆构造，速度变化较为复杂，高速层和低速层交错叠放，浅部和深度分别有两处小尺度速度异常体(如图 3所示)，常规基于射线追踪的层析成像方法对其难以识别，该模型可以用来测试本文方法的有效性.以此模型作为真实模型，震源采用25 Hz的Ricker子波，共模拟68炮地震记录，炮间距260 m，最大偏移距8000 m，道间距30 m，采样间隔4 ms，记录长度为4 s.

根据本文所述频带选择方法，选择了两个频带，主频分别是3 Hz和10 Hz.先用wiener低通滤波器对数据和震源子波做滤波，使其以3 Hz为主频，以图 5所示的大尺度平滑速度场作为初始速度模型，剖分得到三角单元数为222390.在CPU集群上使用18个节点进行并行运算，反演迭代了30次，耗 时4小时38分，得到如图 7所示的反演结果，在地表 5000 m、7000 m和13000 m处的三个CMP点(下面分别简称其为A、B、C点)切出三条速度垂向变化曲线，用来对比检验反演效果，如图 8所示.

图5 反演Overthrust模型的初始速度模型 Fig.5 The initial velocity model used to invert Overthrust model

图6 基于初始速度模型的非规则三角网格 Fig.6 Irregular-grid based on the initial velocity model

图7 主频3 Hz数据的速度反演结果 Fig.7 The velocity model obtained after inversion using 3 Hz peak-frequency data

从图 7可以直观地看到速度模型的基本结构形态均已反演出来，但分辨率不高，断层的断面和断点不清晰，浅部和深部的两处小尺度速度异常体没有反演出来.如图 8所示，厚的速度层反演效果较好，如B点深度1000 m和3000 m处，C点深度2000 m和3000 m处；而薄层反演效果较差，如A点2800 m处，位于一个厚的高速层中的薄低速层没有反演出来，这是因为主频3 Hz的数据包含的是速度体的低波数信息，对高波数信息不敏感.

图8 主频3 Hz反演结果在5000 m (a)、7000 m (b)和13000 m (c)处的垂向速度曲线对比 Fig.8 Comparision of vertical velocity graphs at 5000 m (a),7000 m (b) and 13000 m (c) after inversion using 3 Hz peak-frequency data

用wiener滤波器将数据和震源子波做低通滤波，使其主频为10 Hz，以3 Hz反演得到的最终速度场作为初始速度模型(图 7)，剖分得到三角单元数为1970634.使用CPU集群上的18个节点进行并行运算，迭代12次，用时38小时46分.反演结果见图 9;同样地，图 10显示了A、B、C三点处速度垂向变化曲线的对比情况.

图9 主频10 Hz数据的速度反演结果 Fig.9 The velocity model obtained after inversion using 10 Hz peak-frequency data

图10 主频10 Hz反演结果在5000 m (a)、7000 m (b)和13000 m (c)处的垂向速度曲线对比 Fig.10 Comparision of vertical velocity graphs at 5000 m (a),7000 m (b) and 13000 m (c) after inversion using 10 Hz peak-frequency data

从图 9可以看到，经过高频带反演，速度模型的细节显现出来，逆冲推覆体的断面、断点清晰可见，浅层和深层的两处小尺度异常体也反演出来了.由速度曲线对比图(图 10)可以看到，反演结果和真实模型匹配得非常好，之前低频带反演没有呈现出来的薄层也得到了很好的识别.

多尺度反演方法不仅能很好地引导全波形反演收敛到全局极小，而且还能使反演更加高效.进行低频带反演效率非常高，可以快速得到反演模型的大部分信息，然后采用较高频带进行很少的迭代即可得到高精度的结果.如果不分频带，以上面的模型为例，25 Hz时剖分单元数为10553018，为3 Hz时剖分单元数的47倍，这么多的三角单元参与计算会使得反演效率很低，且极可能因为速度初始模型中低波数成分的缺失导致反演收敛到局部极小值.经过两个频带的反演，最终结果与真实模型之间只有微 小的差异，如果要进一步提高反演精度可以采用更高频带的数据.

图 11是目标函数值随反演迭代次数的变化趋势图，该图显示出反演过程收敛很快，低频带反演在第十次迭代时即得到了很好的收敛，高频带部分迭代4次后即得到了很好的收敛.高频带反演的快速收敛是因为低频带反演给其提供了较好的初始模型.总的来说，L-BFGS优化方法的使用提供了二阶导数信息，使得反演收敛加速.

图11 归一化数据残差随迭代次数变化曲线图 Fig.11 Normalized residual versus iteration number

本文结合自适应变尺度三角网格剖分和格子法正演，实现了二维时间域多尺度全波形反演计算.非规则三角网格的使用既减少了内存需求又提高了正演效率；自适应最优化网格剖分技术能让网格随着模型参数的更新而自适应变化，使网格粒度与反演分辨率相匹配，且避免了可能出现的数值频散而影响反演结果；采用的频带选择方案配合多尺度反演策略在减少数据频率信息冗余度的同时，只需二到三个频带即可使得反演收敛；L-BFGS的使用为反演过程补充了目标函数二阶导数信息，减少了反演迭代次数.二维Overthrust数值算例的成功实施验证了本文方法的有效性.

本文这套全波形反演方法虽然是在二维下实施的，但较容易扩展到三维，应用前景较好.今后我们会深入研究密度及震源对波形反演的影响，且进一步通过炮集编码和GPU加速提升反演效率，以推动本文方法面向工业应用.