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表示检波点的位置，检波点位于区域Ω的边界上，T^obs(r) 表示r处的实际观测走时，c(x) 表示空间位置x处的介质速度. t(r) 表示在当前速度模型下空间位置r处计算的理论走时，在高频假设下，走时场t(r) 满足程函方程：

(2)

其中，t(x_source)=0表示震源点处的地震波走时为零，符号∇ 为梯度算子.

利用伴随状态法，可以确定上述目标函数(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方法中，伴随场的计算依赖于模型表面的法向量，这是有悖射线理论实际的.

在传统AST方法中，在有检波点的边界处，伴随场λ(x) 满足(4)式的边界条件n·λ(x)∇t(x)=t(r)－T^obs(r)，显然与检波点位置处的模型边界法向量有关，这是不合理的.因为这样就忽视了地震波走时只依赖于地震波所经过的介质、与观测点处的边界曲率无关这一客观事实.这样，就会引起计算的伴随场不合理，从而导致梯度的畸变，造成传统伴随状态法走时层析成像反演结果的不正确，我们在2019年就已经发现了这个现象(Han et al., 2019).

举一个单点激发、单点接收的例子.走时场t(x)和唯一检波点处的时差t(r)－T^obs(r) 确定且不为零.如果在这个唯一的检波点处的地表法向量与走时场在该处的梯度∇t(x)|_x=r垂直(这种情况一定存在)，这时，传统方法的伴随方程的边界条件n·λ(x)∇t(x)=t(r)－T^obs(r) 永远不成立！因此，根据传统AST方法的边界条件无法得到正确的伴随场(特别是在地表起伏剧烈的情况下)，最终导致计算的梯度不正确.

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

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

Waheed等(2016)在采用传统AST方法反演的文章附录中指出，在计算伴随场时，井中接收时采用∇·[λ(x)∇t(x)]=0的处理方式(他们的文章中没有说明原因，也没有进行井中观测的试验)，但在地面检波点处仍然还是使用了传统AST方法中伴随方程，结果仍然依赖于地表的法向量.地面和井中检波器采用了不同的处理方式，造成井中和地表的检波器对梯度的贡献不是一个数量级，致使地表走时数据基本不起作用，这是不合理的.因为，尽管地表与井中观测数据的照明范围不同，但走时观测数据的地位应该是等同的.

2 改进的AST方法

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

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

(7)

其中，r_i表示第i个检波点的位置(共有N个检波点)，x表示计算区域Ω′内任意位置，使用δ函数控制区域Ω′内检波点处目标函数的计算. T^obs(x) 表示炮点至空间位置x处的实际观测初至波走时，t(x) 表示在当前模型下炮点至空间位置x处计算的理论走时，它满足(2)式的程函方程.

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

(8)

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

(9)

(10)

(11)

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

(12)

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

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

(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方法中新的伴随方程的正确性.

3 模型试验

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

3.1 简单起伏地表模型

对于传统AST方法，从(4)式计算伴随场的边界条件n·λ(x)∇t(x)=t(r)－T^obs(r) 可以看出：在边界上存在检波点的位置，伴随场初值与理论走时的梯度在该点地表的法向的投影的乘积，等于该点的初至波走时残差.因此，传统AST方法的伴随场以及梯度计算结果依赖于地表法向量.对于地表平缓的模型，地表的法向差别不大，地表法向的存在对反演结果影响并不大；但对于地表起伏剧烈、地表的法向变化较大的模型，这个问题就很突出.

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

模型离散为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方法的伴随场与地表法向量无关，它可以有效修正梯度，从而最终能得到正确的反演结果.

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

3.2 复杂起伏地表模型

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

在地表 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).

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

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

4 结论

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

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

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

对目标函数(8)式中的理论走时t做的扰动，并舍去高阶项，有：

即：

由于扰动方向的任意性，目标函数L对变量t的梯度可以写为:

(A1)

根据关系式，并令，则(A1)式可以写为：

(A2)

对公式(A2)的第二个积分项使用高斯散度定理，得到：

(A3)

其中，n为扩展的边界的外法向单位矢量.

当目标函数达到极小值时，由于目标函数对走时场变量t的导数为零，即(10)式的. 再考虑到扰动方向的任意性，根据(A3)式，可以得到边界上的边界条件：

(A4)

和区域Ω′内部的反传播方程：

(A5)

由于区域Ω及其边界均位于扩展的区域Ω′之内，因此，无论在区域Ω内还是在边界上，伴随场λ(x) 均满足(A5)式，即伴随方程为：

该式即为改进后的AST方法的伴随方程(13)式.

可以看到，在实际需要计算的区域范围内，无论是在模型边界上还是在模型内部，改进后的伴随方程具有统一的形式，而且与地表的法向量无关.