改进的伴随状态法初至波走时层析成像方法

2021-03-08董良国张建明韩佩恩

地球物理学报 2021年3期

董良国，张建明，韩佩恩

同济大学海洋地质国家重点实验室，上海 200092

0 引言

反演地下介质速度分布是地震学的核心问题之一.鉴于初至波最先到达、能量强、信息可靠、易于拾取等优点，目前被大量研究和应用，并产生了利用初至波不同信息的多种反演方法，如初至波走时反演(Olson, 1989)、初至波包络反演(敖瑞德等，2015)、初至波相位反演(Choi and Alkhalifah, 2013)、初至波波形反演(Sheng et al., 2006)以及初至波多信息联合反演(Liu and Zhang, 2017; 张建明等，2020)等方法.

从理论上讲，利用初至波动力学信息可以取得更高的反演精度，但影响初至波动力学信息的因素比较多，信息的稳定性差，导致反演的非线性强，对初始模型要求更高，也更容易陷入局部极值.而初至波走时信息基本上只受速度的宏观分布影响，信息可靠性高，因此初至波走时反演方法稳定性更强，在近地表以及井间速度建模等方面得到了更广泛的应用.

使用不同的地震波正演引擎，产生了不同的初至波走时反演方法，例如，基于射线追踪(Bishop et al., 1985)、基于程函方程走时计算(Taillandier et al., 2009)以及基于波动方程模拟(Luo and Schuster, 1991)的初至波走时反演方法.在这些方法中，不需要射线追踪、基于伴随状态法的初至波走时层析成像方法(Adjoint State Tomography method，以下简称AST方法)近年来得到了很大的关注.这种方法基于伴随状态法、通过程函方程计算的走时场和伴随方程计算的伴随场来计算目标函数的梯度，避免了射线追踪和Fréchet导数矩阵的计算.该方法对内存的要求主要依赖于地下模型的网格数，与观测的地震数据量无关，内存需求极低，特别适合于并行计算(Noble et al.，2010；Benaichouche et al.，2015).至于伴随状态法，自从它首次在优化控制领域被提出以来(Lions, 1971)，目前已经被广泛应用于包括地震学在内的许多领域中(Plessix，2006；Fichtner and Trampert, 2011).

基于程函方程的AST方法最早由Sei和Symes(1994)所提出，自此之后，几乎所有的关于AST方法的理论方法及应用的文献，都是基于面积分来定义目标函数，由此得到的伴随方程都依赖于地表的法向量，本文称这种方法为传统AST方法.

随后，这种基于面积分来定义目标函数的AST方法得到了不断的发展和完善.Leung和Qian(2006)重新推导得到了与观测系统地表法向量有关的伴随状态法走时层析成像方法，构造了伴随方程的差分格式.Taillandier等(2009)以管量场概念对反传播方程的计算格式和意义作了进一步说明，但从其中的地表形态不对称的理论模型试验结果中，可明显看到伴随场的不对称现象，进而导致了梯度形态畸变的不合理现象.Noble等(2010)指出了伴随状态法层析成像技术在并行计算时所具有的巨大效率优势.Bretaudeau等(2014)进一步发展了基于高斯-牛顿优化方法的二阶伴随状态法走时层析方程.为了在走时反演中有效刻画地下模型的界面信息，在AST方法中引入了水平集技术(Li and Leung, 2013; Li et al., 2014)，可以考虑多传播路径问题.这种依赖于地表边界法向量的传统AST方法，目前已经发展到透射和反射走时的同时反演(Huang and Bellefleur，2012；Li et al., 2014)以及各向异性介质模型参数的反演(Waheed et al., 2014,2016)中.

在国内，谢春等(2014)讨论了基于射线理论的伴随状态法初至波走时层析成像技术，并扩展至基于有限频的AST走时层析成像(谢春等, 2015).为优化反演效率和提高反演效果，李勇德等(2017)提出了使用近似Hessian矩阵进行预条件的AST走时反演层析成像方法.

我们已经发现，这种传统AST方法无论在理论方法上还是在应用过程中还存在一定的问题(Han et al., 2019)，本文的目的是对这些问题进行进一步的讨论.首先对这种传统AST方法简要介绍之后，从理论模型试验以及地震波传播物理规律的角度指出了这些问题的存在，并提出了改进的方法，即不依赖于地表法向量的改进的AST走时反演层析成像方法，并通过两个模型反演试验证明了改进方法的正确性和有效性.至于传统AST方法存在的问题的严格数学证明，是我们下一步的工作目标.

需要指出的是，本文采用快速扫描法(Zhao, 2005)求解程函方程和伴随方程，确定走时场和伴随场，从而计算AST方法的迭代梯度.

1 目前传统AST方法存在的问题

1.1 传统AST方法简介

观测和理论走时达到最佳匹配是走时反演的基本原则.目前传统伴随状态法初至波走时层析普遍采用面积分来定义目标函数(Leung and Qian，2006；Taillandier et al.，2009；Noble et al.，2010；Huang and Bellefleur，2012; Li and Leung, 2013; Li et al., 2014；Bretaudeau et al.，2014；Waheed et al., 2014; 谢春等，2014，2015；Benaichouche et al.，2015；Waheed et al., 2016；李勇德等，2017)：

(1)

其中，r表示检波点的位置，检波点位于区域Ω的边界∂Ω上，Tobs(r)表示r处的实际观测走时，c(x)表示空间位置x处的介质速度.t(r)表示在当前速度模型下空间位置r处计算的理论走时，在高频假设下，走时场t(r)满足程函方程：

(2)

利用伴随状态法，可以确定上述目标函数(1)的梯度(Leung and Qian，2006；Taillandier et al.，2009；Noble et al.，2010；Li and Leung, 2013; Li et al., 2014；Waheed et al., 2014, 2016)：

(3)

其中，伴随场λ(x)满足：

(4)

(5)

(4)式中的n为边界∂Ω上位置r处的边界外法向单位矢量.

利用(4)式，根据检波点处的走时残差设置检波点处的伴随场λ的初值后，再利用方程(5)将走时残差反向传播即可获得伴随场λ(x)，进而通过公式(3)就可以得到目标函数对速度的梯度，然后通过线性搜索算法可以很方便地获得一个优化的步长αn(其中n为迭代次数)，这样就可以利用下面的局部优化算法进行速度迭代，得到最终的速度反演结果.

(6)

可见，这种AST方法的核心是利用伴随方程计算伴随场λ(x)，关键是伴随方程.自从Sei和Symes(1994)提出AST走时反演方法之后，有关AST方法的所有文献中，所用的伴随场方程的边界条件均是(4)式，该式依赖于边界的法向量.

1.2 传统AST方法存在的问题

下面具体分析上述传统AST方法存在的问题.

问题一：传统AST方法采用面积分来定义目标函数，决定了其无法处理检波器在模型内部的情况.

方程(1)采用面积分来定义目标函数，就预先假定了所有检波器均位于模型表面，这就造成了传统AST方法只适用于检波点在边界上、无法考虑检波器在模型内部的情况.当进行井中观测以及井间观测时，检波器在模型内部，这时方程(1)的面积分定义目标函数的方式就不适用，伴随方程(4)中的所谓的边界法向当然无从谈起.可见，对于检波点在地下的情况(如井中观测)，这种传统的AST方法是不适用的.

问题二：传统AST方法中，伴随场的计算依赖于模型表面的法向量，这是有悖射线理论实际的.

设计一个除地表形状外中间对称的高速异常体模型(图1)，初始速度也是两边对称的纵向梯度模型(图2).在该模型的中心地表激发，在两侧具有相同横向距离和同等深度、但地表形态不同的两个位置(图1中两个R点处)进行接收.由于炮点两侧真实和初始速度模型以及两个观测点都是完全对称的，只要地表起伏不是极端剧烈，从地震波传播的客观现实来讲，第一轮迭代的梯度在激发点两侧也应该是对称的.但是，在首次迭代中，尽管左右两个检波点具有相同的走时差t(r)-Tobs(r)，但由于左右两个检波点处的地表法向量n不同，因此，利用依赖于地表法向量的传统AST方法得到的伴随场两侧并不对称(图3)，由此得到的梯度场当然也不对称(图4).从这个简单的例子可以看出，传统AST方法的伴随方程存在明显的不合理之处.

图1 不同地表法向量的对称接收点试验的真实模型Fig.1 The true velocity model with two symmetrical receivers but different surface normal vectors

图2 初始速度模型Fig.2 The initial velocity model

图3 传统AST方法得到的归一化后的伴随场分布Fig.3 The normalized adjoint field distribution by the conventional AST method

图4 传统AST方法计算的梯度Fig.4 The gradient obtained by the conventional AST method

问题三：当地表和井中同时观测时，传统AST方法无法正确计算梯度.

2 改进的AST方法

若使用传统AST方法中的方程(1)来定义走时层析目标函数，检波点只能位于区域Ω的边界∂Ω上.我们将∂Ω′面定义为∂Ω面向外的一个拓展面，拓展的空间部分充填空气，速度为340 m·s-1.这样，原来地表上以及井中的检波点均变为区域Ω′的内部点，如图5.由于实际地表之上本来就是空气层，并且在区域Ω的边界∂Ω为凸边界的条件下，这样的拓展并不会改变区域Ω内的地震波走时.

图5 计算区域的拓展示意图Fig.5 The extended sketch map of the computational domain

现在基于拓展的新空间Ω′及其边界∂Ω′，利用体积分定义伴随状态法初至波走时层析目标函数：

(7)

将控制方程(2)与伴随变量λ(即乘子)引入目标函数(7)中，形成新的目标函数：

(8)

在最优化过程中，走时场变量t、伴随变量λ和速度c相互独立，所以当目标函数达到极小值时，目标函数仍满足以下三个关系：

(9)

(10)

(11)

由(9)式可以得到控制方程，即程函方程(2)式.而由(11)式可以得到目标函数(7)的梯度：

(12)

可见，由于每一次迭代中的当前速度模型c(x)已知，伴随场λ(x)的求取是确定目标函数梯度的关键.而由(10)式就可以导出控制伴随场λ(x)的伴随方程.

我们使用扰动法求解(10)式，得到了改进的AST方法的伴随方程(推导过程见附录A)：

x∈Ωand ∂Ω.

(13)

由于每一次迭代中的理论走时场t(x)可以通过求解程函方程(2)式确定，这样，求解(13)式的伴随方程即可确定伴随场λ(x)，再利用(12)式即可获得目标函数的梯度，进而利用(6)式就可以进行速度迭代反演.

需要指出的是，我们这种改进的AST方法与传统AST方法的唯一差别，就是确定伴随场λ(x)的伴随方程不同.传统AST方法的伴随方程是(4)和(5)式，利用(4)式设置边界上检波点处的伴随场λ的初值，再利用(5)式确定伴随场λ(x).而我们这种改进的AST方法的伴随方程是统一的(13)式，地表和地下检波点具有相同的处理方式，并且伴随场的分布与模型边界的法向量无关，自然克服了前面指出的传统AST方法存在的三个缺陷.

也就是说，改进的AST方法适用于各种地表和井中同时观测方式，地下和井中观测数据的处理方式得到了统一，从理论方法上避免了Waheed等(2016)处理井中观测时的井中和地表观测数据对梯度贡献不一致的问题.无论是内部点还是边界上，改进方法的伴随方程都遵从统一形式，伴随场的计算与模型边界的法向量无关，不需要区分检波点位置处的地表法向量，这样更符合实际物理意义，保证了梯度场的正确性.

仍然采用图1所示的起伏地表高速异常体模型，初始模型仍然为图2所示的纵向常梯度模型，使用上述改进的AST方法计算得到的平滑后的首轮迭代梯度如图6所示.可以发现，炮点两侧的梯度场是对称的，这是符合理论预期的，证明了我们得到的改进AST方法中新的伴随方程的正确性.

图6 改进的AST方法计算的归一化后的梯度分布Fig.6 The normalized gradient obtained by the improved AST method

3 模型试验

下面，通过两个理论模型试验来说明改进的AST方法的正确性及反演效果.

3.1 简单起伏地表模型

为此，设计一个地表起伏比较剧烈的理论模型(见图7a)，在2000 m·s-1的均匀速度模型内部存在一个速度为3000 m·s-1高速异常体.在地表等间距激发的201炮产生的地表处的初至波走时数据作为“观测数据”，相邻炮点的水平间距为20 m.每炮都在地表均匀设置401个检波器，相邻检波点水平间距为10 m.

图7 简单起伏地表速度模型及改进前后AST反演结果(a) 简单起伏地表速度模型； (b) AST-IV反演结果； (c) AST-DV反演结果.Fig.7 The simple topographic velocity model and the AST inversion results(a) The simple topographic velocity model; (b) The AST-IV result; (c) The AST-DV result.

模型离散为401×101网格点，纵横向网格间距均为10 m，采用速度值为2000 m·s-1的常速模型作为反演的初始模型.分别使用改进的不依赖于地表法向量的AST方法(AST-IV)和依赖于地表法向量的传统AST方法(AST-DV)进行试验，反演结果分别见图7b和图7c.其中，AST-IV和AST-DV分别迭代了693代和575代，它们的目标函数都不再下降.

可以看到，改进方法比传统方法能更精确地反演出近地表附近的局部高速异常，而传统方法的反演结果非常差.由于二者的唯一差别，就是确定伴随场λ(x)的伴随方程不同，说明传统AST方法在伴随方程的求解过程中，检波点处的走时残差反向传播是不正确的，导致其计算的梯度和模型修改量不正确，说明了传统AST方法的伴随场依赖于地表法向量是不合理的.

图8展示了改进AST方法与传统AST方法不同迭代次数(第1、20、100、300次)的归一化梯度值.可以发现，随着迭代的进行，改进AST方法的梯度形态逐渐变得与高速异常体形态一致，而传统AST方法的梯度没有得到本质的改善，反演过早就陷入了局部极值.这进一步证实了传统AST方法在伴随方程的求解过程中，检波点处的走时残差反向传播是不正确的，导致其梯度和模型修改量形态异常，造成该问题的根源就是来自于传统AST方法的伴随场错误地依赖于地表法向量，当地表起伏剧烈时，这个问题就很突出.而改进的AST方法的伴随场与地表法向量无关，它可以有效修正梯度，从而最终能得到正确的反演结果.

图8 传统与改进的AST方法不同迭代次数的归一化梯度从上到下分别为第1、20、100、300代的梯度(左侧为传统的AST方法，右侧为改进的AST方法).Fig.8 The normalized gradients of different iterations for the conventional and the improved AST methodsFrom top to bottom show the gradients corresponding to 1, 20, 100 and 300 iterations. The AST-DV and the AST-IV results are shown on the left and right respectively.

从本试验的目标函数收敛曲线(图9)上可以发现，改进AST方法不仅收敛速度明显快于传统AST方法，且反演结果的走时残差远小于传统AST方法的走时残差.这主要是由于传统AST方法的伴随场错误地依赖于地表法向量，对于这个地表起伏剧烈的模型，造成了梯度不正确，使反演过早地陷入了局部极值.

图9 归一化后的目标函数收敛曲线虚线为传统AST方法，实线为改进的AST方法.Fig.9 The normalized objective function iteration curvesThe dashed and the solid lines represent the conventional and the improved AST methods respectively.

3.2 复杂起伏地表模型

下面使用一个相对复杂的模型来验证改进方法的反演效果.图10显示了该模型5～35 km范围内的近地表附近的速度变化.该模型深部结构比较简单，近地表速度横向变化很大，还存在一些低速和高速异常体，30 km的横向范围内地表最大起伏只有0.5 km，起伏还是比较平缓的.采用的纵向常梯度初始速度模型如图11所示，模型离散为4001×151网格点，纵横向网格间距均为10 m.

图10 起伏地表真实速度模型Fig.10 The true complex topographic velocity model

在地表5～35 km范围内共激发751炮，相邻炮点的水平间距为40 m.每炮都在地表设置500个检波器，均采用双边接收方式，相邻检波点水平间距20 m，最大偏移距5 km.基于准确速度模型、利用最短路径法(Moser，1991)计算的走时作为“观测”的初至波到达时.分别使用依赖于地表法向量的传统AST方法和改进的不依赖于地表法向量的AST方法进行试验，两种方法在反演中均使用快速扫描法(Zhao, 2005)计算走时场以及伴随场.

两种方法的反演结果分别见图12和图13.由于本模型的地表起伏相对比较平缓(30 km的横向范围内地表最大起伏只有0.5 km)，传统和改进的AST方法都反演出了近地表速度的基本变化趋势，但改进方法比传统方法能更精确地反演出近地表附近的高速和低速局部异常.另外，在传统方法中有较严重的能量聚焦现象(图12)，在同一速度层中，不同的速度参数呈条带状分布，横向速度波动较大，说明在传统AST方法伴随方程的求解过程中，检波点处的走时残差反向传播是不正确的，导致其梯度和模型修改量形态异常，而这些不合理现象在改进的AST方法中得到了较好的改善(图13).

图11 初始速度模型Fig.11 The initial velocity model

图12 传统AST方法反演结果Fig.12 The inversion result by the conventional AST method

图13 改进AST方法反演结果Fig.13 The inversion result by the improved AST method

对两种反演方法第80代的反演结果，抽取地表以下20、30、40、50 m深度的速度剖面(见图14)，可以发现，无论是什么深度，改进后AST方法(红线)的反演精度均高于传统方法(蓝线)，改进方法的反演结果更加接近真实模型.

图14 第80代反演结果地表下不同深度的速度剖面对比真实模型(黑实线)、初始模型(黑虚线)、传统AST方法反演结果(蓝线)、改进AST方法反演结果(红线).Fig.14 The inverted lateral velocity slices with depth of (a) 20 m, (b) 30 m, (c) 40 m and (d) 50 m below the surface for the 80th iterationThe black solid line, the black dashed line, the blue line and the red line represent the true model, the initial model, the AST-DV inversion result and the AST-IV inversion result respectively.

从炮点在15、20、25 km处的三个单炮的走时残差上也可以明显看出(图15)，改进AST方法的反演结果的初至波走时残差更小，说明改进后AST方法与传统AST方法相比，反演精度得到了提高.

图15 方法改进前后反演结果的剩余走时对比(a) 炮点在15 km处的单炮; (b) 炮点在20 km处的单炮; (c) 炮点在25 km处的单炮. 真实模型(黑线)、传统AST方法(蓝线)、改进AST方法(红线).Fig.15 The residual first-arrival traveltime of the final inversion results for different shots with sources located at (a) 15 km, (b) 20 km and (c) 25 kmThe black line, the blue line and the red line represent the true model, the AST-DV inversion result and the AST-IV inversion result respectively.

4 结论

本文分析了目前传统的伴随状态法初至波走时层析成像方法的不足：一是伴随方程依赖于地表法向量不合理，二是无法合理处理井中观测问题.在此基础上，论文进一步提出了不依赖地表法向量的改进的伴随状态法初至波走时层析成像方法，从理论上得到了不依赖于地表法向量的伴随方程，克服了传统方法中伴随场依赖于地表法向量的缺陷，使得检波点处的走时残差可以正确地反传播至地下，进而可以得到更加合理的速度修正方向，提高了速度反演的精度.同时，提出的改进方法还可以适应地表和井中同时观测的任意观测系统.两个理论模型反演试验结果，说明了改进的AST方法的正确性和有效性.

附录A 改进的AST方法伴随方程的推导

下面使用扰动法(Nocedal and Wright, 2006)求解方程(10)，从而得到改进的AST方法的伴随方程.