1 层析算法

基于射线理论的地震走时层析成像在工业界获得了成功应用, 现在仍是商业软件的核心建模模块。其基本原理建立在对正问题的一阶Taylor级数展开的基础上:

$ t=f(s) \approx f\left(s_{0}\right)+f^{\prime}\left(s_{0}\right) \Delta s $

(1)

式中:t为射线的走时; f(s)为表达走时与慢度s之间非线性关系的函数; s₀为背景慢度。在射线理论框架下, 正问题可以用矩阵公式表达为:

$ f(\mathit{\boldsymbol{s}}) = \mathit{\boldsymbol{L}}(\mathit{\boldsymbol{s}}) \cdot \mathit{\boldsymbol{s}} $

(2)

式中:L(s)表示慢度模型s下的射线路径长度矩阵。根据Fermat原理, 射线路径是相对稳定的, 因此:

$ {f^\prime }\left( {{\mathit{\boldsymbol{s}}_0}} \right) = \mathit{\boldsymbol{L}}\left( {{s_0}} \right) $

(3)

将(3)式代入(1)式即得传统的射线走时层析方程组:

$ \mathit{\boldsymbol{L}}\left( {{\mathit{\boldsymbol{s}}_0}} \right)\Delta \mathit{\boldsymbol{s}} = \Delta \mathit{\boldsymbol{t}} $

(4)

式中:Δs为慢度扰动向量; Δt为真实走时与理论合成走时之间的走时残差向量。传统层析成像通过求解上述方程组实现局部慢度向量的更新。如此, 利用局部线性化+迭代方式实现非线性走时层析反演。

受计算能力的限制, 20世纪70—80年代以来, 地震层析成像一直采用经验性的求解大规模线性方程组(4)的方式进行模型的更新^[1]。典型的算法包括反投影法(BPT)、代数重建法(ART)和联合代数重建法(SIRT)等。这类方法基于加权平均的方式对模型进行更新, 计算效率高但反演精度低。另一类反演类算法包括最小二乘QR分解(LSQR)和截断奇异值分解(SVD)等, 基于奇异矩阵的广义逆算法, 反演效果更好, 但计算量比较大, 算法不稳定^[2]。但无论经验性解法还是广义逆解法都需要存储层析矩阵, 对计算机内存要求很高。

借鉴全波形反演的实现过程^[3], 在目标函数的局部最优化基础上^[4]建立了一种走时层析成像实现方式, 该方式通过寻找某一慢度向量s^*, 使得如下目标函数极小:

$ \phi (\mathit{\boldsymbol{s}}) = \frac{1}{2}\left\| {{\mathit{\boldsymbol{t}}_0} - {\mathit{\boldsymbol{t}}_c}} \right\|_2^2 \to \min $

(5)

式中:t₀为观测走时; t_c为满足(2)式的理论合成走时。结合(3)式得到上述目标函数的梯度^[5]:

$ \mathit{\boldsymbol{g}} = \frac{{\partial \phi }}{{\partial \mathit{\boldsymbol{s}}}} = - {\mathit{\boldsymbol{L}}^{\rm{T}}}\Delta \mathit{\boldsymbol{t}} $

(6)

不难推而广之, 只要采用(5)式的目标函数形式, 任何反演应用都可以得到如(6)式的目标函数梯度表达式, 区别在于L取其一般形式为依赖于正问题的敏感核函数K, Δt取其一般形式为数据残差Δd。因此, 按照下列公式沿着最速下降方向(也可以采用其它更优的更新方向)即可实现模型的更新。

$ {\mathit{\boldsymbol{s}}^{k + 1}} = {\mathit{\boldsymbol{s}}^k} - \alpha \mathit{\boldsymbol{g}} $

(7)

式中:k为反演中的迭代次数; α为步长, 可以利用线性搜索法计算得到^[6]。不难看出, 目标函数梯度的计算是关键, 理论上可以根据(6)式算得, 但需要计算并存储核函数矩阵, 在实际应用中不现实。

随着伴随状态法在全波形反演(FWI)中的成功应用^[7], 越来越多的学者将伴随状态法应用到层析成像中^{[4, 8-10]}。伴随状态法走时层析将程函方程作为传统目标函数的约束项, 利用拉格朗日乘子法进行求解, 进而间接得到目标函数梯度的表达式:

$ \frac{{\partial \phi }}{{\partial s}} = \int_\Omega {\frac{{\lambda (x)}}{{v(x)}}} {\rm{d}}x $

(8)

式中:v为地下介质速度分布; λ为伴随旅行时场, 在地表处满足:

$ \lambda \mathit{\boldsymbol{n}} \cdot \nabla t = {t_0} - {t_c} $

(9)

在介质内部满足:

$ \nabla \lambda \cdot \nabla t = 0 $

(10)

式中:n为地表法向量。层析成像无需射线追踪即可直接求得目标函数的梯度, 配合快速扫描正演算法, 计算效率高, 内存占用小。但无法实施预条件, 为此给出了一种等效的预条件算子^[10]。

LIU等^[11]在全波形反演研究中提出了一种改进的散射积分算法, 之后又将该算法应用于地震走时层析成像^[12-13]。该算法与伴随状态法皆用于局部最优化算法中梯度的求取, 但改进的散射积分法将梯度的计算表达为核函数—残差乘的累积求和形式:

$ \begin{aligned} &\boldsymbol{p}_{\mathrm{sD}}=-\boldsymbol{g}=\boldsymbol{L}^{\mathrm{T}} \Delta t=\\ &\left[\begin{array}{ccccc} {l_{11}} & {\cdots} & {l_{i 1}} & {\cdots} & {l_{m 1}} \\ {\vdots} & {\ddots} & {\vdots} & {\ddots} & {\vdots} \\ {l_{1 j}} & {\cdots} & {l_{i j}} & {\cdots} & {l_{m j}} \\ {\vdots} & {\ddots} & {\vdots} & {\ddots} & {\vdots} \\ {l_{1 n}} & {\cdots} & {l_{i n}} & {\cdots} & {l_{m n}} \end{array}\right]\left[\begin{array}{c} {\Delta t_{1}} \\ {\vdots} \\ {\Delta t_{i}} \\ {\vdots} \\ {\Delta t_{m}} \end{array}\right]=\sum\limits_{i=1}^{m}\left[\begin{array}{c} {l_{i 1}} \\ {\vdots} \\ {l_{i j}} \\ {\vdots} \\ {l_{i n}} \end{array}\right] \Delta t_{i} \end{aligned} $

(11)

式中:p_SD为最速下降方向, 即负梯度方向; l_ij为敏感核函数矩阵L中的元素, 表示第i条射线在第j个模型网格里的长度; Δt_i为Δt向量中的元素, 表示第i个数据的走时残差。可见, 该方法尽管需要计算核函数, 但无需存储核函数, 内存占用极小、计算效率高、物理意义明确。更重要的是, 由于计算了核函数, 预条件与Hessian矩阵的利用非常方便。如利用公式(12), 获得梯度的同时就可以得到准确的预条件算子, 基本不需额外的计算量。这是伴随状态法层析无法实现的。

$ \boldsymbol{H}_{0}=\operatorname{diag}\left\{\boldsymbol{L}^{\mathrm{T}} \boldsymbol{L}\right\}=\sum\limits_{i=1}^{m}\left|\begin{array}{ccccc} {l_{i 1}^{2}} & {\cdots} & {0} & {\cdots} & {0} \\ {\vdots} & {\ddots} & {\vdots} & {\ddots} & {\vdots} \\ {0} & {\cdots} & {l_{i j}^{2}} & {\cdots} & {0} \\ {\vdots} & {\ddots} & {\vdots} & {\ddots} & {\vdots} \\ {0} & {\cdots} & {0} & {\cdots} & {l_{i n}^{2}} \end{array}\right| $

(12)

式中:H₀为Hessian矩阵的对角线元素, 即准确的预条件算子。图 1展示了预条件算子对层析成像反演的影响, 尤其当初始模型梯度严重小于真实模型梯度时, 无预条件的层析在浅层形成了高速屏蔽层, 导致深部的速度结构很难被反演出来。

2 基于不同理论的地震层析成像方法

根据(4)式, 地震层析成像包含3个要素:核函数、物性参数与数据。这3个要素相互关联, 某些情况下又可以相互独立。核函数依赖于表达物性参数与数据之间关系的方程(正问题), 因此根据所基于的正演理论不同, 层析成像可分为射线层析、菲涅尔体层析与波动层析。传统的地震层析成像基于射线理论^[14-16], 射线本身即为层析的核函数。由于其简单、高效, 在层析静校正^[17]、背景模型建立^[18-19]中应用广泛。但受其高频假设的限制, 其反演结果分辨率不高^[20]。

LUO等^[21]提出了基于波动方程的走时层析成像方法。该方法利用波形互相关构造走时残差, 然后利用伴随状态法求解最小二乘目标函数, 迭代更新模型。随后, MARQUERING等^[22]提出了有限频地震层析成像方法。该方法基于波动理论构建核函数, 同样利用波形互相关构造走时残差, 在天然地震学研究中得到了成功应用。但这些方法需要波动方程的正演模拟计算, 计算量庞大。

为此, SPETZLER等^[23]和LIU等^[24]对走时菲涅尔体的性质进行了研究, 刘玉柱等^[25]进而提出了带限菲涅尔体地震层析成像方法。该方法基于波动理论显式构建核函数((13)式), 但只考虑对单一震相贡献最大的菲涅尔体范围内的波路径(图 2)。

$ K(r, \omega | g, s)=\frac{2 \omega}{v_{0}^{3}(r)} \cdot \operatorname{Im}\left[\frac{G(g, r) u(r, s)}{u(g, s)}\right] $

(13)

式中:ω为圆频率; Im表示取复数的虚部; G(g, r)为无扰动速度场v₀(r)中r点在检波点g处的格林函数; u(r, s)为炮点s点在空间任意一点r处的格林函数。

后续LIU等^[13]又将改进的散射积分法引入梯度的计算, 将射线走时残差替代为波形互相关残差, 计算效率显著提高, 内存占用量大幅降低, 反演精度有所提高。菲涅尔体层析与射线层析的理论模型反演效果如图 3所示, 对比图 3b和图 3c可见, 反演中考虑波动理论核函数的确可以提高反演精度。图 4和图 5展示了西部山区某实际资料的散射积分法菲涅尔体层析与射线层析反演结果及其对应的层析静校正叠加剖面^[13], 可以看到菲涅尔体层析的优势明显。

3 利用不同初至信息的地震层析成像方法

传统的走时层析成像利用的是射线到达时, 即初至走时信息。该走时与射线路径核函数相匹配。随着有限频理论和波动理论在层析成像中的应用, 波形互相关走时替代了射线走时, 被用于构建走时残差。但无论是射线走时还是互相关走时, 都是单一走时, 反演精度和分辨率提高有限。

为了提高反演分辨率, SHENG等^[26]提出了初至波形反演, 由于相对于单一走时信息, 初至波形包含的信息量更多, 因此反演精度和分辨率显著提高。但该方法需要进行波动方程的正演模拟, 计算量与传统的全波形反演相当, 因此难以应用到实际当中。刘玉柱等^[27]提出了利用高斯束合成初至波形, 以提高初至波形反演计算效率的方法, 取得了一定的效果。

为了避开时间域波动方程的正演模拟, 同时利用初至波形中更多的信息, CHOI等^[28]提出了瞬时走时层析成像方法, 刘玉柱等^[29]提出了初至波相位走时层析方法。这些方法基于波动理论, 利用初至中频率依赖的相位信息实现初至波多频率走时反演, 反演精度介于射线层析与波形层析之间, 但计算效率较波形反演明显降低。图 6给出了利用射线走时、瞬时走时与初至波形的层析反演结果, 可以看出随着利用的信息量的增加, 反演精度逐步提升。相对于射线层析, 波形层析的反演精度有很大的提升。但值得一提的是, 随着信息量的增加, 层析反演的非线性程度也在增强, 即对初始模型的依赖性也逐次增强。因此, 实际应用中可以通过逐步增加信息量来逐级提高反演的精度。

图 7a展示了某海上OBS实际资料的初至波形层析反演结果。图 7b是与该反演速度无关的商业软件叠前时间偏移(PSTM)剖面, T23轴和T25轴与图 7a中深度为2.0km和4.7km的界面匹配得非常好^[30]。

4 从单参数到多参数

传统的层析成像方法以反演P波速度为主, 但这已经不能满足成像、储层识别等对地下介质多参数反演的需求。因此, 多参数层析成像研究逐渐多起来, 包括各向异性层析^[31-33]、Q层析^[34-36]、纵横波速度与密度层析^[37-38]等。但多参数层析的主要难点在于多参数耦合, 即不同的参数对观测数据都有影响。虽然理论上可以导出并分析不同参数的敏感核函数, 但同时反演多参数仍然很困难, 尤其存在强弱参数之分的情况下, 弱参数往往很难被反演出来。研究表明, 除利用必要的先验信息作为约束外, 还需要考虑以下3个方面。

1) 参数化方法^[39-42]。采用合适的参数化方式可以有效降低多参数之间的强弱差异, 有望实现多参数的同时反演。如VTI介质走时层析成像中, 速度属于强参数, 各向异性参数属于弱参数, 但如果采用三速度参数化方式则可以使得对应的核函数在相同的数量级上, 这对反演很有利^[5]。

2) 反演策略^[43-45]。尽管多个参数对地震数据都有影响, 但不同的参数对不同的数据子集可能有不同的影响权重, 采用针对性的反演策略往往可以收到事半功倍的效果。如刘玉柱等^[5]在各向异性走时层析反演中采用了两步法反演策略, 即第一步同时反演速度和各向异性参数, 然后将反演的速度和最开始的各向异性参数作为初始模型再进行第二轮反演, 得到可靠的各向异性参数(如图 8所示)。

3) Hessian矩阵的利用^{[11, 46-49]}。Hessian矩阵是目标函数对模型参数的二阶偏导数构成的系数矩阵。在单参数反演中, Hessian矩阵的对角线元素反映了波场的几何扩散效应, 非主对角线元素反映了空间不同散射体对地震波传播的干涉效应; 在多参数反演中, 如果将Hessian矩阵视为块状矩阵, 那么Hessian矩阵的主对角块元素反映的是单参数对波场的以上影响, 非主对角块元素反映的则是包含几何扩散和干涉效应在内的多参数耦合作用。如在变密度声波方程全波形反演中, Hessian矩阵可以表达为^[50]:

$ \begin{aligned} &\boldsymbol{H}=\left(\begin{array}{l} {\boldsymbol{K}_{v}^{\mathrm{T}}} \\ {\boldsymbol{K}_{\rho}^{\mathrm{T}}} \end{array}\right)\left(\boldsymbol{K}_{v} \quad \boldsymbol{K}_{\rho}\right)=\left(\begin{array}{cc} {\boldsymbol{K}_{v}^{\mathrm{T}} \boldsymbol{K}_{v}} & {\boldsymbol{K}_{v}^{\mathrm{T}} \boldsymbol{K}_{\rho}} \\ {\boldsymbol{K}_{\rho}^{\mathrm{T}} \boldsymbol{K}_{v}} & {\boldsymbol{K}_{\rho}^{\mathrm{T}} \boldsymbol{K}_{\rho}} \end{array}\right)=\\ &\left(\begin{array}{ll} {\boldsymbol{H}_{v v}} & {\boldsymbol{H} v \rho} \\ {\boldsymbol{H} \rho v} & {\boldsymbol{H} \rho \rho} \end{array}\right) \end{aligned} $

(14)

式中:K_v, K_ρ分别表示速度和密度的核函数矩阵; T表示复矩阵的共轭转置。因此, 在多参数反演中考虑Hessian矩阵可以实现多参数的同时反演, 提高多参数的反演精度。但遗憾的是, 目前在近地表物性参数层析反演(尤其是多参数反演)中较少使用完整的Hessian矩阵。

5 结论与展望

关于地震层析成像的研究持续了三十多年, 从射线理论到波动理论, 从走时反演到波形反演, 从单参数到多参数, 从方程求解到局部最优化, 对其理论、方法与应用研究从简单到复杂, 反演精度与分辨率逐步提高。从以上综述不难看出, 影响反演精度最主要的因素是利用的信息, 波形反演要远比走时反演获得的反演分辨率高; 其次是理论基础, 由于地震波传播的有限频特性, 相较于射线层析, 波动理论反演可以获得更高的反演精度; 最后是反演算法、反演策略、参数化等对反演质量也具有一定的影响。

多参数反演是目前层析成像的研究热点, 除Hessian矩阵的有效利用外, 反演策略的制定也很关键, 值得深入研究。另外, 本文提出的改进的散射积分算法方便预条件、Hessian矩阵利用充分、内存占用小、反演稳定等诸多优点使其尤其适合于求解地震层析成像反演问题。相信它会有更加广泛的实际应用价值。