0 引言

地震旅行时成像是以地震波在不同介质中传播规律为基础，在地表或井中观测地震信号，通过旅行时反演重构地质体内部速度分布的一种地球物理勘探方法。该方法广泛应用于工程检测、油气勘察、地球内部结构研究等^[1-3]。地震初至波旅行时成像能识别、追踪目标地质体，稳定性较好，对浅层成像精度较高。

地震旅行时成像的过程主要包括正演和反演两部分。早期的地震旅行时正演采用经典射线追踪方法，如打靶法^[4]、高斯射线束法^[5]、弯曲射线法^[6]，但可能存在多条射线路径，反演极易出现多解。为此，后续发展了基于程函方程的地震旅行时算法，如有限差分算法^[7]、快速推进算法（Fast Marching Method，FMM）^[8]、最短路径算法^[9]等。其中FMM是一种网格数值算法，通过有限差分或有限元方法求解程函方程追踪地震波前界面，具有较高的精度和效率，因而被广泛应用^[10]。由于地震资料观测的局限性、观测数据中存在误差等原因，地震旅行时反演是一个不适定问题，反演算法主要有阻尼最小二乘法^[11]、基于正交分解的最小二乘法（LSQR）^[12]、子空间反演法^[13]、代数重建技术（ART）^[14]等，这些算法在不同程度上能降低反演问题的不适定性。此外，正则化技术^[15]通过引入地球物理先验信息减少了反问题的多解性、提高了解的稳定性，因此被广泛用于地球物理反演^[16]。最小结构模型约束是将最小模型稳定泛函与最大光滑稳定泛函加入到正则化项中，可有效降低反演问题的不适定性，但地质体边界刻画模糊，尤其是对射线分布稀疏的地质体分辨率较低。

近年来，模糊C均值（Fuzzy C-means，FCM）聚类算法已用于地球物理反演^[17-20]。该算法通过反复修改聚类中心及每个网格单元对聚类中心的隶属度，实现对反演结果的自动分类^[17]。此外，通过将先验信息作为参考聚类中心的方式，将FCM聚类算法从无监督学习变成监督学习，引导每个网格单元向先验信息的方向修正，可有效提高反演结果的精度^[19-20]。

本文在地震初至波旅行时成像中，采用非结构三角形进行网格剖分，采用Sethian等^[8]提出的FMM算法进行正演计算，在地震旅行时正演过程中直接计算灵敏度但不存储射线路径^[16]。采用以模型灵敏度信息为依据的渐进反演网格优化技术^[21-22]，在最小结构模型约束正则化反演基础上，引入FCM聚类模型约束，实现了基于FCM聚类模型约束的二维地震初至波旅行时渐进反演；对比了无监督学习与先验信息监督学习反演结果，并应用简单模型讨论了FCM聚类模型约束的权重、先验信息引导项权重等参数选取策略；最后，通过实测数据的反演验证了方法的有效性和实用性。

1 方法原理 1.1 地震初至旅行时反演

在最小结构模型约束正则化技术^[23]的基础上，引入FCM聚类模型约束项，构建地震初至波旅行时反演的目标函数

$ \begin{aligned} \varphi(\boldsymbol{m}) & =\varphi_{\mathrm{d}}(\boldsymbol{m})+\lambda \varphi_{\mathrm{m}}(\boldsymbol{m})+\beta \varphi_{\mathrm{FCM}}(\boldsymbol{m}) \\ & =\left\|\boldsymbol{W}_{\mathrm{d}}\left[\boldsymbol{d}^{\mathrm{obs}}-\boldsymbol{A}(\boldsymbol{m})\right]\right\|_2^2+ \\ & \lambda\left\|\boldsymbol{W}_{\mathrm{m}}\left(\boldsymbol{m}-\boldsymbol{m}^{\mathrm{ref}}\right)\right\|_2^2+\beta \varphi_{\mathrm{FCM}}(\boldsymbol{m}) \end{aligned} $

(1)

式中：m为模型速度向量；φ_d（m）为初至波旅行时观测数据d^obs与正演响应A（m）间的数据拟合项；φ_m（m）为模型拟合项；φ_FCM（m）为FCM模型约束项；λ为正则化因子；β为FCM聚类模型约束项的权重；W_m为模型误差矩阵；m^ref为参考模型向量；W_d为数据方差矩阵，形式为

$ {\boldsymbol{W}}_{\mathrm{d}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\frac{1}{{\sigma }_{1}+\zeta }, \frac{1}{{\sigma }_{2}+\zeta }, \cdots , \frac{1}{{\sigma }_{N}+\zeta }\right) $

(2)

式中：σ_i为的数据方差；N为观测的初至波旅行时数据个数；$ \zeta $为极小的正有理数，保证分母不为0。

本文采用最小结构稳定因子^[24]作为模型拟合项

$ \begin{array}{l}{\varphi }_{\mathrm{m}}\left(\boldsymbol{m}\right)=\int \left[{\alpha }_{\mathrm{s}}{\left(\boldsymbol{m}-{\boldsymbol{m}}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right)}^{2}+\right.\\ \begin{array}{ccc}& & \end{array}\left.{\alpha }_{x}{\left(\frac{\partial \boldsymbol{m}}{\partial x}\right)}^{2}+{\alpha }_{z}{\left(\frac{\partial \boldsymbol{m}}{\partial z}\right)}^{2}\right]\mathrm{d}s\end{array} $

(3)

式中：α_s为最小模型稳定泛函的比例系数；α_x、α_z分别为最大光滑稳定泛函在x和z方向的比例系数；∂m/∂x和∂m/∂z分别表示模型向量在x和z方向的偏导数；积分是对任意单元的面积分。

1.2 FCM聚类模型约束

FCM聚类算法的思想是在迭代过程中，反复修改聚类中心和聚类中心的隶属度，是一种无监督学习。FCM聚类的优化目标函数为

$ {\varphi }_{\mathrm{F}\mathrm{C}\mathrm{M}}\left(\boldsymbol{m}\right)=\sum\limits _{l=1}^{C}\sum\limits _{j=1}^{M}{u}_{j, l}^{q}{\left({m}_{j}-{v}_{l}\right)}^{2} $

(4)

式中：C为聚类中心个数；M为模型单元总数；m_j为第j个模型单元速度；隶属度u_j，l表示第j个模型单元属于第l个参考簇的程度；模糊化参数q控制隶属度u_j，l模糊化程度，本文取q=2；v_l为第l个簇中心的速度值。u_j，l为非负数，满足

$ \sum\limits _{l=1}^{C}{u}_{j, l}=1 \quad(0\le j\le M) $

(5)

为了改善反演效果，将先验信息引入FCM聚类目标函数φ_FCM（m）^[17]，将聚类算法从无监督学习升级为监督学习，即

$ \begin{array}{l}{\varphi }_{\mathrm{F}\mathrm{C}\mathrm{M}}\left(\boldsymbol{m}\right)=\sum\limits _{l=1}^{C}\sum\limits _{j=1}^{M}{u}_{j, l}^{q}{\left({m}_{j}-{v}_{l}\right)}^{2}+\\ \quad\begin{array}{cccc}& & & \sum\limits _{l=1}^{C}{\kappa }_{l}{\left({v}_{l}-{t}_{l}\right)}^{2}\end{array}\end{array} $

(6)

式中：t_l为第l个参考聚类中心；κ_l为第l个参考聚类中心的权重因子。

φ_FCM（m）的第一项可表示为

$ \sum\limits _{l=1}^{C}\sum\limits _{j=1}^{M}{u}_{j, l}^{q}{\left({m}_{j}-{v}_{l}\right)}^{2}=\sum\limits _{l=1}^{C}(\boldsymbol{m}-{\boldsymbol{v}}_{l}{)}^{T}{\boldsymbol{U}}_{l}(\boldsymbol{m}-{\boldsymbol{v}}_{l}) $

(7)

式中

$ \left\{\begin{array}{l}{\boldsymbol{v}}_{l}={\left(\begin{array}{c}{v}_{l}\\ {v}_{l}\\ ⋮\\ {v}_{l}\end{array}\right)}_{M\times 1}\\ {\boldsymbol{U}}_{l}={\left(\begin{array}{cccc}{u}_{1, l}^{q}& & & \\ & {u}_{2, l}^{q}& & \\ & & \ddots & \\ & & & {u}_{M, l}^{q}\end{array}\right)}_{M\times M}\end{array}\right. $

(8)

1.3 目标函数优化求解

本文反演采用高斯—牛顿法求解模型参数向量$ \boldsymbol{m} $。假设第n+1次迭代，对正演算子A（m）做一阶泰勒展开，并将第n+1次迭代目标函数φ^（n+1）（m）对Δm^（n）进行求导并取0，可得高斯—牛顿方程

$ \begin{array}{l}\left\{{\left[{\boldsymbol{J}}^{\left(n\right)}\right]}^{\mathrm{T}}{\boldsymbol{W}}_{\mathrm{d}}^{\mathrm{T}}{\boldsymbol{W}}_{\mathrm{d}}{\boldsymbol{J}}^{\left(n\right)}+\lambda {\boldsymbol{W}}_{\mathrm{m}}^{\mathrm{T}}{\boldsymbol{W}}_{m}+\beta \sum\limits _{l=1}^{C}{\boldsymbol{U}}_{l}^{\left(n\right)}\right\}\mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)}\\ ={\left[{\boldsymbol{J}}^{\left(n\right)}\right]}^{\mathrm{T}}{\boldsymbol{W}}_{\mathrm{d}}^{\mathrm{T}}{\boldsymbol{W}}_{\mathrm{d}}\left[{\boldsymbol{d}}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{A}}^{\left(n\right)}\left(\boldsymbol{m}\right)\right]-\\ \begin{array}{c}\end{array}\lambda {\boldsymbol{W}}_{\mathrm{m}}^{\mathrm{T}}{\boldsymbol{W}}_{m}\left[{\boldsymbol{m}}^{\left(n\right)}-{\boldsymbol{m}}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right]-\beta \sum\limits _{l=1}^{C}{\boldsymbol{U}}_{l}^{\left(n\right)}\left[{\boldsymbol{m}}^{\left(n\right)}-{\boldsymbol{v}}_{l}^{\left(n\right)}\right]\left(9\right)\end{array} $

(9)

式中：$ {\boldsymbol{W}}_{\mathrm{m}}^{\mathrm{T}}{\boldsymbol{W}}_{\mathrm{m}}={\alpha }_{\mathrm{s}}\boldsymbol{I}+{\alpha }_{x}{\boldsymbol{W}}_{x}^{\mathrm{T}}{\boldsymbol{W}}_{x}+{\alpha }_{z}{\boldsymbol{W}}_{z}^{\mathrm{T}}{\boldsymbol{W}}_{z} $，其中I为单位矩阵，W_x、W_z为水平和垂直方向粗糙度矩阵^[24]；J^（n）为第n次迭代的灵敏度矩阵，是正演算子A^（n）（m）对m的偏导数矩阵；第n+1次迭代的簇中心$ {v}_{l}^{(n+1)} $和隶属度$ {u}_{j, l}^{(n+1)} $分别为

$ {v}_{l}^{(n+1)}=\frac{\sum\limits _{j=1}^{M}{\left({u}_{j, l}^{q}\right)}^{\left(n\right)}{m}_{j}+{\kappa }_{l}{t}_{l}}{\sum\limits _{j=1}^{M}{\left({u}_{j, l}^{q}\right)}^{\left(n\right)}+{\kappa }_{l}} $

(10)

$ {u}_{j, l}^{(n+1)}=\frac{1}{\sum\limits _{k=1}^{C}{\left[\frac{‖{m}_{j}^{\left(n\right)}-{v}_{l}^{\left(n\right)}‖}{‖{m}_{j}^{\left(n\right)}-{v}_{k}^{\left(n\right)}‖}\right]}^{\frac{2}{q-1}}} $

(11)

式（9）可写为

$ {\boldsymbol{H}}^{\left(n\right)}\mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)}={\boldsymbol{P}}^{\left(n\right)} $

(12)

采用稳定双共轭梯度法（BICGSTAB）^[21]求解$ \mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)} $

$ \mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)}={\left[{\boldsymbol{H}}^{\left(n\right)}\right]}^{-1}{\boldsymbol{P}}^{\left(n\right)} $

(13)

再根据线性搜索步长更新模型参数

$ {\boldsymbol{m}}^{(n+1)}={\boldsymbol{m}}^{\left(n\right)}+\eta \mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)} $

(14)

式中η为搜索步长，可表示为

$ \eta =\frac{{\left[\mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)}\right]}^{\mathrm{T}}{\boldsymbol{P}}^{\left(\mathrm{n}\right)}}{{\left[\mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)}\right]}^{\mathrm{T}}{\boldsymbol{H}}^{\left(n\right)}\mathrm{\Delta }{\boldsymbol{m}}^{\left(n\right)}+\xi } $

(15)

式中$ \xi $为极小的正实数，确保分母不为0。

1.4 基于模型灵敏度的自适应反演

本文采用变网格技术实现地震初至波旅行时反演。在反演的初期采用较粗的网格，通过减少单元数降低反演的多解性，再利用模型灵敏度信息进行自适应网格加密^[22]。以上一轮反演结果作为网格细化后反演的初始和参考模型，确保反演的稳定性，同时改善反演效果。定义模型灵敏度为

$ {S}_{j}=\sqrt{\frac{1}{N}\sum\limits _{i=1}^{N}{J}_{ij}} $

(16)

式中J_ij为灵敏度矩阵的元素。模型灵敏度S_j表示单元的模型参数m_j改变时对全部观测数据集的影响。对于测点处单元的模型灵敏度较大，网格的细化会从测点逐渐向反演区域的深部进行，因此具有对反演区域整体优化的特点。基于模型灵敏度的自适应反演过程如下。

（1）读入观测数据，设置每重网格最大的迭代次数、最小单元面积、网格优化比例（为反演单元总数的2%）。

（2）从数据信息生成正演模型的几何描述，设置反演范围，并做非结构三角形网格剖分。

（3）对当前网格下使用BICGSTAB最优化算法求解高斯—牛顿优化的反演目标函数，直至满足BICGSTAB的迭代次数，并求取模型改变量。

（4）通过模型灵敏度网格优化方案计算单元模型灵敏度信息，按照网格优化比例选取需要细化的网格单元，进行网格的自适应加密。

（5）将上一轮反演结果作为细化后反演的初始和参考模型，并判断是否达到反演终止的条件。满足则终止反演，否则返回第三步开始网格细化后的反演。

1.5 参数选取

地震初至波旅行时反演的目标函数φ（m）中的正则化因子λ是平衡数据拟合项与模型拟合项的关键参数。本文采用李曼等^[25]提出的改进“L曲线”算法求取λ。该算法无需利用导数信息，可以自动舍弃错误的曲率点，能提高正则化因子的求取精度。

引入FCM聚类模型约束的同时，也引入了聚类中心的个数C、参考聚类中心t_l、参考聚类中心的权重因子κ_l、FCM聚类模型约束项权重β等众多参数。这些参数势必会影响反演的稳定性及效果。因此如何选取合适参数是引入FCM聚类模型约束反演的一个关键问题。其中，聚类中心的个数$ C $及参考聚类中心t_l可通过获取有效且精确的先验信息确定。参考聚类中心的权重因子κ_l的作用是引导聚类中心向真值靠近，FCM聚类模型约束项权重β的作用是控制反演中FCM聚类模型约束项的比重。为此，对κ_l和β参数不同组合的反演效果进行了对比分析，讨论了两个参数在反演过程中的作用；结合渐进反演网格优化技术的特点，提出了一种FCM聚类模型约束参数的自动选取方案。该方案在选取一个适当的聚类中心权重因子κ=100后，通过减少反演前期FCM聚类模型约束比重，让正则化项占主导，后期增大FCM聚类模型约束比重实现FCM聚类模型约束。将每重网格的FCM聚类模型约束项权重修改为

$ \beta =100\times K $

(17)

式中K为每重网格的迭代次数。

2 算例分析 2.1 模型算例1

为了验证算法的有效性，设计了图 1左所示的理论模型。在背景速度为2000 m/s的均匀半空间中，设计了两个尺寸为30 m×30 m的高速异常体，速度均为3000 m/s，顶面埋深分别为10和90 m。

采用双井激发、井中和地面同时接收的观测系统。井中深度2.5 m处开始设置炮点和检波器各60个，地表均匀布设60个检波器，炮点距和检波点距都为2.5 m，井间距为120 m。基于上述参数进行地震初至旅行时正演，将理论旅行时加入5%的随机噪声作为观测数据（共14400个）。

反演的初始速度模型是速度为2000 m/s的均匀半空间，初始正则化因子λ设置为10000。最小结构模型约束反演结果如图 1右所示，可以大致判断两个高速异常体的位置，浅部的异常体的规模和速度与实际值接近，但深部异常体规模较大，速度与实际值相差亦较大。

为了分析基于模型灵敏度渐进反演过程，对于图 1右所示反演模型，设计了四重网格，第一重网格的迭代次数为5，其他的最大迭代次数为10，反演过程的网格变化如图 2所示，从测点处逐渐向反演区域内部细化。表 1给出了每重网格的反演单元个数、迭代次数及反演开始与结束的均方根误差（RMS）。单元数由1961逐渐增加到3883；初始网格和第一次细化网格反演的RMS下降很快，后两次的RMS下降量减少并趋于稳定。

首先根据精确的先验信息，设置了FCM聚类模型约束反演的两个参考聚类中心：v=2000和3000 m/s。其次，为分析参考聚类中心的权重因子κ（本文算例中不同参考聚类中心的权重因子设为相同）和FCM聚类模型约束项的权重因子β对反演结果的影响，设置κ和β分别为1、10、100、1000、10000，对比不同参数组合的反演结果（图 3）。

由图 3可见：

（1）当聚类中心权重因子κ=1或10时，随β增大，反演结果的聚类特征更明显，但速度值较真值差异变大（图 3的第1、第2行）。

（2）当FCM聚类模型约束项权重β=1、10、100时，随聚类中心权重因子κ的增大，异常体的特征改善不明显，也没有明显的聚类特征（图 3的第1~第3列）。原因是聚类中心权重因子κ较大时引导每个单元的速度向真实值靠近，但控制聚类程度的FCM聚类模型约束项权重β较小，聚类效果不明显。

（3）适当的参数组合可以得到理想的反演结果，如κ=100、β=1000或10000时的反演结果异常体尺寸与真实模型接近（图 3第3行的第4和第5列），而β=10000时反演的异常体速度与真实更接近（图 3第3行的第5列）。

（4）越大的κ和β组合的反演结果的聚类特征越明显（图 3第5列的第4和第5行），但异常体的尺寸偏小。

由以上分析可知，适当的聚类中心权重因子κ能有效地引导聚类中心向真值靠近，而FCM聚类模型约束项权重β可以控制聚类程度。因此，本文采用前文提出的FCM聚类模型约束参数的选取方案，即通过减小反演前期FCM聚类模型约束项比重β，让模型正则化项占主导，后期增大FCM聚类模型约束项的权重实现FCM模型约束。新参数选取方案的反演结果如图 4所示，与未加入FCM聚类模型约束的反演结果（图 1右）相比，背景更干净，两个高速异常体尺寸、位置及速度值与实际更吻合，边界刻画清晰。与图 3第3行第5列的反演结果接近。

为进一步分析FCM聚类模型约束的地震初至波旅行时反演效果，对最小结构模型约束和加入FCM模型约束反演结果的速度频数进行统计（图 5）。可见，FCM聚类模型约束通过引入先验信息，反演结果中出现两个峰值与实际给定的两个聚类中心接近，形成了有效的聚类效果。

2.2 模型算例2

为进一步验证基于FCM聚类模型约束的地震初至波旅行时反演能提高成像分辨率，设计了如图 6左所示的层状模型，区域1、区域2和区域3是速度分别为1000、2000、3000 m/s的均匀介质。采用与模型算例1相同的观测方法和参数进行正演计算。

反演的初始模型是速度为2000 m/s的均匀半空间，初始正则化因子λ为10000，设置了v=1000、2000、3000 m/s三个参考聚类中心；采用模型算例1中FCM聚类模型约束参数的选取方案。

图 6中为最小结构模型约束的反演结果，可以大致判断层状速度带，但各层边界模糊，界面附近的速度值与实际相差较大；图 6右为加入FCM聚类模型约束的反演结果，与图 6中相比，能明显区分3个区域的层状速度带，界面位置较准确，刻画清晰，速度值与真值更吻合。

图 6两种方法反演结果的速度频数统计如图 7所示，加入FCM聚类模型约束后，三个峰值的速度与引入的三个先验信息吻合，形成了有效的聚类特征。

2.3 实测数据反演

使用强风化层中的孤石（微风化或未风化的花岗岩）探测的井中地震初至波旅行时数据验证本文方法的有效性和实用性。实测数据由单井激发、单井接收观测系统采集。在孤石两侧各布置一个深度为20 m的钻孔，钻孔间距为9.4 m，共布置17个炮点和36个检波器，左侧从3.0 m深度开始布设炮点，炮点距为1 m；右侧从深度2.5 m开始布置检波器，道间距为0.5 m。共采集有效旅行时观测数据612个。

在测线5 m处的验证钻孔揭示：0~2.6 m为人工堆积黏性土；2.6~7.0 m为坡积而成的粉质黏土；7.0~9.8 m为花岗岩风化残积而成的砂质黏土；9.8~12.2 m为花岗岩微风化体（孤石）；12.2~16.5 m为花岗岩风化残积而成的砂质黏性土；16.5~33.6 m为全风化层，局部夹强风化花岗岩；稳定的潜水面深度为2.5 m。

图 8左为未加入FCM聚类模型约束的反演速度剖面，可以大致判断孤石（高速异常体）位置，与钻孔信息吻合，但尺寸较小，且速度与背景差异不够明显。图 8右为加入FCM聚类模型约束后的反演速度剖面，为了凸显高速异常体的明显特征，设置了v=2000和3000 m/s两个参考聚类中心，采用模型算例1中FCM聚类模型约束参数的选取方案。由图 8右能清晰地判断孤石的位置及边界，异常体的尺寸和速度值都与实际钻孔信息吻合。

3 结束语

为了提高地震初至波旅行时对浅层结构的成像精度，本文提出了基于FCM聚类模型约束的地震初至波旅行时反演方法。通过理论研究、模型数据和实测数据试算，取得以下几点认识：

（1）FCM聚类模型约束易与传统正则化反演结合，可有效提高反演结果的分辨率；

（2）FCM聚类模型约束通过引入精确的速度先验信息，可有效提高速度反演精度；

（3）FCM聚类模型约束权重、先验信息引导项权重直接影响反演结果，本文的自动选择策略有效改善了反演效果，降低了参数选择难度。