1 引言

立体层析是一种极具特色的层析反演方法，该方法将局部相干同相轴在炮道集与检波点道集内的射线参数(或P参数)水平分量、以及炮点坐标和检波点坐标都纳入到反演的数据空间之中，重新定义了层析反演的数据分量和模型分量，使得立体层析成为射线类层析反演方法中唯一一种可以同时反演速度结构、反射点位置与反射层局部形态的方法(Billette and Lambare, 1998).国内外学者在立体层析反演的理论和实践方面有很多工作(Prieux et al., 2013; 倪瑶等，2013；倪瑶等，2014；任丽娟等，2014；李振伟等，2014).Li等(2015)将块状约束加入立体层析FRECHET导数矩阵的建立过程中，使得反演结果更具有地质一致性.王宇翔等(2016)将结构张量算法用于P参数水平分量的搜索，实现了高密度的二维立体层析反演.杨锴等(2016)基于非降阶Hamilton算子实现了三维直角坐标系下的立体层析反演方法.长期以来，改进传统立体层析中过于光滑的正则化项一直是学术界和工业界共同关注的问题.实践表明，引入地层格架信息并追求地质一致性的反演结果是最强有力的正则化手段.这种策略可以极大压缩立体层析反演层析矩阵的规模，使得反演的精度和可靠性大为提高.

从射线方法的特点出发，若引入地层格架信息应该使用更优化的模型描述方式.三角剖分建模方式正是这样一种可以自然体现地层格架信息的模型描述方式.为保证反演得到的速度模型具有地质一致性，本文使用了地质界面约束下的慢度平方模型三角网格剖分方式(Hale and Cohen, 1991).二维情形下的模型剖分对应于三角网格剖分，三维情形下的模型剖分则对应于四面体网格剖分.若不加特殊说明，本文提到三角网格剖分皆为二维模型的描述.作者首次提出并完成了慢度平方三角剖分模型下的二维立体层析反演方法.不仅强波阻抗界面在三角剖分模型内可以得到自然描述，而且一个地质块体内的速度采用慢度平方描述使得射线追踪存在解析解，相比传统的Runge-Kutta步进算法更为简洁高效.基于考虑界面的射线扰动理论(Farra and Madariaga, 1989)，我们完成了二维立体层析数据空间对慢度平方三角剖分模型中的所有模型分量的偏导数求取.本文提供的理论算例表明，在慢度平方三角剖分模型下建立的二维立体层析FRECHET矩阵规模得到了大幅压缩，迭代次数得到降低的同时求解精度也得到了很好的保证.由于慢度平方三角剖分模型是现有已知的最为稀疏的模型描述方式，因此在慢度平方三角剖分模型下实现立体层析可以认为是是目前已知的最为高效的二维立体层析层析算法.

2 三角网格剖分模型中二维立体层析数据空间与模型空间简介

首先以二维立体层析为例对涉及到的数据空间和模型空间各个分量给予介绍.图 1显示了二维立体层析数据空间的各个分量与模型空间的各个分量.对于二维空间的一个射线对而言，假设地表水平，那么地表能观测到的数据为d=(S_x, R_x, T_sr, P_sx, P_rx), 其中S_x, R_x为炮检点的X坐标；T_sr为这个射线对的总走时；P_sx, P_rx为炮检点处的射线参数水平分量.这些数据分量与地下的模型空间m=(x₀, z₀, θ_s, θ_r, V)有关，其中x₀, z₀为地下反射点坐标, θ_s, θ_r代表从炮点与检波点一侧出射的地下散射角, V代表地下介质速度信息.众所周知，共炮点道集和共检波点道集内的局部相干同相轴的斜率等同于检波点处和炮点处慢度矢量的水平分量(Billette and Lambare, 1998).事实上，立体层析反演中“立体”这个称谓的由来，即与炮检点处的局部传播方向被引入数据空间，可以用于约束层析反演密切相关.

不同于二维常规立体层析的模型空间，三角网格慢度平方介质下的模型空间为m=(x₀, z₀, θ_s, θ_r, x_n, z_n, u²).其中x_n和z_n是控制界面节点的横、纵坐标；u²是可以线性变化的慢度平方即u²=u₀²+k_xx +k_zz, 其中k_x, k_z是慢度平方的空间梯度.数据空间的残差Δd与模型空间残差Δm之间的关系由射线扰动理论建立，这里简写为

(1)

(1) 式代表了一个矩阵方程，其系数矩阵F即为立体层析FRECHET偏导数矩阵，在慢度平方介质下其完整形式可写为

(2)

在矩阵中，显然有

(3)

其次，根据费马原理，走时T在速度结构、激发点与接收点位置都固定时已经是最小走时.可以认为散射角对走时T不存在一阶扰动影响，即有

(4)

将上述关系代入(2) 式，可以得到简化后的三角网格立体层析矩阵(5)：

(5)

(5) 中σ的含义为传统二维立体层析中针对各类数据的尺度加权因子，用于平衡不同类型的数据对于立体层析反演的误差泛函的贡献比例，避免由于物理量纲不同造成数据类型量级之间差别过大的问题，参数的大小需要通过做具体的测试得到.对理论数据而言，无论传统立体层析还是本文的三角网格立体层析中这个量级差别都是比较固定的.李振伟(2014)和王宇翔(2016)处理实际数据时所用的σ参数与倪瑶等(2013)用于理论数据的σ参数并无不同.作者预计三角网格立体层析也应该可以类似的将σ参数推广到实际数据情形.即便到了三维多参数下，这个策略应该依然奏效，当然，用自适应的方式获取σ参数应该完全可行.只要在立体层析矩阵建立之后，按照五类数据分量进行FRECHET导数的量级评估，然后根据评估结果即可确定最优的σ参数，在反演实际数据时自适应获取σ参数应该是一个更好的选择.

立体层析反演的理论基础是射线扰动理论(Farra and Madariaga, 1987).它的重要性在于，在傍轴射线追踪和Hamilton传播算子基础上建立了速度扰动与描述射线传播的相空间各参数扰动量之间的一阶线性关系.三角网格剖分由于在模型描述中引入了界面，因此需要在扰动理论中考虑界面扰动对立体层析数据空间的影响.Farra等(1989)提出的考虑界面扰动的射线扰动理论恰恰是三角网格立体层析算法中建立FRECHET敏感度矩阵所需要的.考虑到文章的自洽性，在第3节中我们首先给出三角网格剖分下的二维射线追踪解析式推导和傍轴射线追踪系统；第4，5节给出基于界面扰动的射线扰动理论建立的二维情形下、三角网格立体层析数据空间分量对模型空间分量之间的一阶线性扰动关系.

3 三角网格慢度平方介质中二维射线追踪及傍轴射线追踪系统 3.1 三角网格慢度平方介质的二维射线追踪

(1) 慢度平方梯度介质中的运动学射线追踪

慢度平方介质中二维程函方程的Hamilton算子形式为

(6)

其中u²(x, z)是随空间位置而变化的慢度平方函数，根据上述Hamilton形式可推导出二维运动学射线方程为

(7)

为了简化表达，(7) 式使用了数字下标形式表示射线P参数与空间坐标的两个分量，其中P₁对应于P_x，P₂对应于P_z，X₁对应于水平x方向，X₂对应于垂直z方向.当u²不是关于空间位置线性变化时，上述微分方程组可用Runge-Kutta法求数值解.但是当u²关于空间位置线性变化时，那么u²关于坐标的偏导数就是常数，这时(7) 式便有了解析解.线性变化的慢度平方可以写为

(8)

其中，u₀²为一个三角块内的参考点处的慢度平方值，k₁(或k_x), k₂(或k_z)为慢度平方在x, z方向上的梯度.从(7)、(8) 式出发容易得到慢度向量的解析解(Virieux et al., 1988)为

(9)

射线坐标的解析解为

(10)

走时的解析解为

(11)

其中积分参量σ的计算公式为

(12)

(2) 射线穿越界面时的处理

当射线从当前三角网格穿过分界面进入另一个三角网格时，需要用到Snell定律.其穿越界面前后的射线参数满足以下关系：

(13)

其中p，为过网格前后的射线参数，n₁为三角网格边的法向量, ×为向量积运算.

3.2 慢度平方介质中的二维傍轴射线追踪

(1) 层内的傍轴射线追踪

假设背景介质中慢度分布为u(x), 在背景介质中传播的射线称为中心射线.定义中心射线上的位置矢量和慢度矢量分别为x_c(σ)和p_c(σ).通过一阶射线扰动理论即可算出邻近的傍轴射线，其位置矢量和慢度矢量由下式给出:

(14)

位置和慢度向量的扰动统一写为δy(σ)=(δx, δp)，其扰动关系可用傍轴射线方程表示如下：

(15)

根据式(7)，(15) 式中的系数矩阵可以进一步简化为

(16)

其中I为单位矩阵，V为慢度平方的二阶导：

(17)

在慢度平方介质中，由于块内慢度平方为线性变化，所以其二阶导V为0，此时傍轴射线追踪亦有非常简洁的解析解表示形式：

(18)

将(18) 式整理为与传播矩阵的乘积形式为

(19)

我们将此传播矩阵简记为P₀矩阵.

(2) 射线穿越界面时的傍轴射线追踪

Farra等(1989)将射线扰动理论推广到了考虑界面的情形下，其中最重要的贡献是详细讨论了当射线通过界面时的傍轴射线追踪系统.如图 2所示，中心射线与界面相交于O点，此时O点对应的傍轴扰动为δx.而相应的傍轴射线上的对应点则位于Q处.一般情况下傍轴射线与中心射线不可能同时打到界面，所以Q点不可能在界面上，即Q点与I点不重合.那么如何根据O点与Q点信息计算出I点处的相空间扰动量？首先通过图 3中的矢量关系可以看出O点与I点之间的矢量差如下：

(20)

这里dσ代表d从Q处至I处对应的积分参量变化，∂H₀/∂P, ∂H₀/∂x代表Hamilton算子在O点处的变化率.这里再加入两个定义：矢量a、b的张量积为|a〉〈b|，定义二者点积为〈a|b〉，把界面定义成f₀(x)=0，▽f₀定义成O点表面处的法向量.很明显I点属于界面(见图 2)，因此I点的位置满足界面的定义：f₀[x₀(σ₀)+dx]=0，▽f₀·dx=0，基于上述定义，即可计算出积分参量的变化

(21)

有了dσ的表达式，那么I点处关于O点的傍轴扰动量dy=(dx, dp)即可写为传播矩阵形式：

(22)

其中

(23)

(23) 式中

(24)

(3) 傍轴射线的连续性条件

Farra等(1989)基于傍轴射线过界面以及慢度矢量的连续性条件，导出了界面扰动关于相空间扰动的表达式.很明显，在界面两侧射线入射点和出射点一致，可以写为

(25)

慢度矢量的连续性条件由界面两侧的扰动Hamilton算子相等导出.显然在界面两侧必须满足，将dH与的表达式代入得到

(26)

这里dp为过界面之前的慢度矢量扰动，为过界面后的慢度矢量扰动，根据Snell定律有

(27)

这里▽f代表界面I点处的法向量(界面的方程定义为f=0), 当界面比较光滑时，存在二阶导

(28)

这里▽▽f₀代表界面x₀(σ_i)处的二阶导数.将式(28) 代入式(27) 中可得：

(29)

联合(25)、(26) 和(29) 式，即可通过射线达到界面之前的dx和dp表达射线通过界面后的：

(30)

经过一系列代数操作，最终可以将过界面后相空间扰动量表达为过界面前相空间的扰动量Δy的函数，注意与Δy同时包括了坐标扰动量和慢度矢量扰动量.具体推导细节请参考(Farra et al., 1989)：

(31)

其中的转换矩阵T由下式定义：

(32)

(32) 式中的子矩阵T_Ⅰ和T_Ⅱ通过下式给出：

(33)

式(31) 代表了当射线穿过一个界面时如何根据矩阵Π和矩阵T求得过界面之后的相空间扰动量.在计算下一个界面的相空间扰动量时，便可把此次扰动量作为下一层的初始扰动量进行计算.图 4显示了利用上述傍轴扰动方程在一个慢度平方介质中计算得到的中心射线与傍轴射线的轨迹.

4 三角网格慢度平方介质中二维立体层析数据空间中x，P_x与模型空间各分量的一阶扰动关系

根据上述Farra等(1989)提出的界面扰动的射线扰动理论即可建立地表坐标x，地表水平慢度P_x与模型空间各分量的一阶扰动关系.扰动由两部分组成，第一部分是相空间初始扰动，即射线出射位置处的坐标x₀，z₀以及出射方向p_x0，p_z0；第二部分是模型参数扰动，其中界面扰动可由P，T，Π这三个传播矩阵进行计算.这里所说的初始相空间扰动既可以理解为直接改变相空间的初始值，也可以理解为射线到达当前地质块之前的模型参数变化而引起的相空间扰动.

假设在背景速度上有一个光滑扰动，定义这种模型的慢度变化为从u₀至u=u₀+Δu，慢度扰动产生相应的Hamilton的扰动响应H=H₀+ΔH，其中H₀为参考慢度u₀所对应的Hamilton算子，且

(34)

其中设Δk_x与Δk_z为慢度梯度扰动.根据Farra和Madariaga(1987)，总的射线相空间信息可以写为y(τ)=y₀(τ)+Δy(τ)，考虑到总扰动量有相空间扰动和背景速度扰动两方面的影响，对相空间进行全微分计算后可以得到

(35)

第一项为初始相空间扰动所引起，第二项为背景速度扰动所引起，且

(36)

其中.根据Gilbert和Backus (1966)，(35) 的传播矩阵解为

(37)

其中传播矩阵P₀的形式已经由(19) 式给出.

4.1 x，P_x与慢度梯度的一阶扰动关系

从以上推导可以看出，在一个地质块内x，P_x与背景慢度的一阶扰动关系为0，但是与慢度梯度k_x, k_z的关系还是存在的.在一个地质块之内，(37) 式可以进一步展开为(38) 式

(38)

(38) 式与直接从解析式(9) 式和(10) 式出发求关于k_x, k_z的偏导数等价.从另一个角度证明了其本身的正确性.但是，地下地质构造一般是由多个地质块组成的.必须根据(35) 式，将(38) 式乘以传播矩阵后得到更一般的表达式(39).该式即为x，P_x与慢度梯度的一阶扰动关系.

(39)

其中.

4.2 x，P_x与反射点，两个散射角的一阶扰动关系

令反射点位置扰动为(Δx₀, Δz₀), 射线的散射角扰动量为Δθ，散射角扰动量不能直接用传播矩阵进行运算，需要先计算出其对应的对初始相空间的扰动量，即：

(40)

再利用传播矩阵即可得到一阶扰动量关系：

(41)

4.3 x，P_x与界面节点坐标的一阶扰动关系

设n_j, j=1, 2为界面左右两点的坐标，界面扰动为Δf, Δf是一个关于节点坐标的函数，这里以线性的界面扰动为例，给出Δf的形式.用界面左右端点(x₁, z₁)，(x₂, z₂)表示界面，其解析式为：

(42)

对左端点在z方向上做Δz₁的扰动，扰动后的界面表示为

(43)

前后两个界面函数相减，即可得到Δf表达形式：

(44)

根据Farra等(1989)，将界面扰动量转换为相空间初始扰动量，再利用传播矩阵即可算得数据空间与界面扰动的一阶扰动关系：

(45)

其中

(46)

5 三角网格慢度平方介质中走时T_sr与模型空间各分量之间的一阶扰动关系

走时T_sr为炮检点接收到的走时之和，即T_sr=t_s+t_r, 显然有：

(47)

分析知与的求法完全一致，两者分别对应着地下反射点到炮点和检波点方向的透射过程.

5.1 走时T与慢度平方，与其梯度的一阶扰动关系

设某层的慢度平方为u₀²，k_x和k_z为慢度平方在x，z方向上的梯度，通过(11) 式即可得到走时T与慢度平方，与其梯度的一阶扰动关系：

(48)

5.2 走时T与反射点坐标的一阶扰动关系

以图 5a为例，当反射点横坐标有Δx₀的扰动时，扰动后射线路径与原射线路径之差可近似为图中AC的距离，由三角关系可以很快地得到AC长度为Δx₀sinθ，除以速度即得到时间差.反射点z方向扰动情况亦可如此推导.最终的走时与反射点坐标位置的扰动关系如下：

(49)

5.3 走时T与界面节点的一阶扰动关系

走时与界面节点的一阶扰动关系在Trinks等(2005)中已经有详细推导.这里直接给出其表达式：

(50)

(51)

其中，θ_i, θ_t分别是入射角与透射角.

至此，三角网格立体层析所需的所有数据空间分量对模型空间分量的偏导数定量关系都已经得到.

6 理论数据敏感度测试

在得到三角网格立体层析所需的数据空间各分量与模型空间各分量之间的一阶线性扰动关系之后，以下基于三角网格理论模型中的一个射线对信息，来测试不同数据分量与不同模型分量之间的敏感度.图 6a展示了一个二维三层慢度平方梯度模型，射线穿过2层界面，3个地层，共5个节点，通过射线追踪得到的射线对，反射点在(3 km，2 km)处，在z=0 km处可观测到数据空间(S_x, R_x, T_sr, P_sx, P_rx).考虑到物理量之间的对称性，这里仅展示数据空间分量(S_x, P_sx, T_sr)对模型空间分量(u₀², k_x, k_z, θ_s, z_n)的偏导数，且其中P_sx/u₀²，S_x/u₀²，T/θ_s为0.

图 6c，6e展示了S_x对模型参数θ_s, z_n的偏导数测试.以图 6e为例，偏导数∂S_x(z_n)/∂z_n通过傍轴射线追踪得出，蓝色曲线S_x(z_n+δz_n)通过射线追踪得出，可以清晰地看出绿线S_x(z_n) + (∂S_x(z_n)/∂z_n)δz_n为蓝色曲线的切向量，由此可以验证偏导数的正确性.同理其他几幅图也同样证明了其偏导数推导的正确性.

7 三角网格立体层析理论数据反演测试 7.1 层状背斜模型理论数据测试

本实验在如图 7a所示的盐丘模型中进行，横向12 km，纵向6 km，中间有6个强波阻抗界面, 分为7个地质块体(块体的编号见图 7b)，每个块体内设为常慢度即梯度系数k_x, k_z均为0，从上到下的6个块体内的真实慢度平方信息以及最终反演得到的慢度平方信息如表 1所示.

从图 7c可以看到实验用的射线信息(仅使用108个射线对).从图 7d所示的初始模型可以看到，不仅两侧的反射界面比真实模型明显呈上翘趋势，界面的中间形态也与真实模型也有差距.初始模型中的慢度平方信息为真实模型的0.8倍.图 7e，7f，7g分别为迭代1，2，4次后的模型.可以看到通过4次迭代后的模型(如图 7g所示)，盐丘模型的形状与慢度信息都得到较好的恢复.无论从图 7h所示的目标泛函下降曲线还是表 1中的慢度反演结果都可以看到迭代到第4轮时候残差就已经很小，证明了反演的正确性和高效性.王宇翔等(2016)曾经基于B样条描述反演过同样的模型，需使用多达1800个B样条系数以及6000个射线对，还要施加很强的沿层光滑的正则化，需要20轮以上的迭代次数才能够反演出比较可靠的结果.相比之下，本文的反演只用到了108个射线对，而且没有加入任何正则化，表明三角网格剖分加慢度平方介质在描述模型的简约性方面具有明显优势.这个优势对于射线照明稀疏的深部构造反演显然是一个利好消息.

7.2 复杂断陷模型理论数据测试

本实验在如图 8a所示的中海油断陷模型中进行，该模型横向24 km，纵向8 km，其有18个界面，分为17个地质块体，对块体编号如图 8b所示，大部分块体慢度平方都带有梯度.真实模型与最终反演结果的背景慢度及纵横向梯度信息如表 2所示.

从图 8c可以看到实验所用的射线信息(仅使用394个射线对).图 8d为初始模型，其几何拓扑结构要求与真实模型一致.所谓几何拓扑一致，简单的说就是初始模型的分块数和真实模型相同.一般而言只要建模时的初始成像剖面品质不很差，这个条件在实际应用中依然是容易满足的.可以看出初始模型的界面比较平坦，相对真实模型有明显差异.另外除表层慢度信息正确外(可以理解为海水或近地表速度已知)，其余区块的初始慢度信息均为真实模型的1.2倍.

图 8e显示了初始模型与真实模型相减之后的残差，可以看到二者的差距相当可观.图 8f，8g，8h分别为迭代1，2，4次后的模型.可以看到4次迭代后的模型，断陷模型的形状与慢度信息都得到较好的恢复.从图 8i所示的目标泛函下降曲线也可以看到迭代到第4轮时候残差已经很小.图 8j显示了最后一轮反演模型与真实模型相减之后的残差，可以看到二者的差距相比初始残差有明显缩小.反演模型和真实模型的每块内部的参考慢度平方慢度差已经到了非常接近的程度，只是因为节点坐标不可能反演到与真实模型完全一致，所以在边界处才有比较大的残差.具体差值见表 2.在每块内横向变化系数k_x都能反演出很好的精度，然而纵向变化系数k_z精度略差.这是因为反射射线能保证反演结果的横向分辨率却无法保证纵向分辨率所致.如果想要在反演结果中得到很好的纵向分辨率，在反射观测时必须有足够多的大偏移距射线或甚至井间观测射线才能够获得，因此这并非立体层析方法的问题，而是观测局限所致.总之上述各个指标都体现了反演算法的正确性.

8 结论与讨论

本文首次提出并完成了慢度平方三角剖分模型下的二维立体层析反演.在慢度平方三角剖分模型下实施立体层析是一种自然的得到地质一致性反演模型的层析反演算法.该方法存在诸多优点.首先，界面在三角剖分模型内可以得到自然的描述，界面和慢度将同时得到有效的更新.其次，由于三角形内部引入慢度平方描述，无论是运动学或动力学射线追踪均存在解析解，无需使用Runge-Kutta算法，计算成本得到大幅削减.最重要的一点，相比常规的B样条描述或者大网格描述模型的方式，由于模型空间的稀疏度大幅提高导致未知数的个数急剧减少，无论对于求解精度的提高还是计算成本的下降都具有极大的帮助.反演获得的地质一致性模型可以用于逆时偏移成像或作为全波形反演的初始模型，也非常有利于后续的地质解释与岩石物理分析.

本文实现的三角网格二维立体层析可以顺利扩展到三维四面体网格剖分立体层析，其中三角网格剖分下的二维FRECHET导数推广到三维是非常直接的，理论方面没有任何困难.三角剖分立体层析在实用中的挑战其实不在于反演算法本身，而在于如何实现高效的三角剖分建模.在三维实际应用中这个挑战尤为明显.因为在三维四面体剖分建模之前，必须先在初始三维成像体上解释出反射界面，而后在反射界面上选择出若干节点后再进行四面体剖分建模.如何在三维初始成像数据体上实现自动的三维地质层序精细解释，实现高效准确的三维四面体建模，是日后面向应用的研究需要重点攻关的一个问题.在当前自动层序解释软件与相关算法迅猛发展的大背景下，基于二维、三维自动层序解释的三角剖分与四面体剖分将是下一步实用化研究要重点攻关的部分.