0 引言

地震资料处理中，精确的速度场不仅是高质量成像的基础，也是储层预测的基础^[1-2]。目前获得地震速度场的常规方法是在速度谱上拾取速度，但人工拾取存在主观误差大、效率低下等问题^[3-4]。随着三维地震数据的规模日益增加，处理员一般不会对工区内每个道集进行速度拾取，而是间隔一定距离进行速度拾取，对于速度变化较快的区域，这种做法会降低速度场的精度，最终影响地震剖面的质量^[5]。为此，许多学者不断努力探索地震速度拾取的自动化方法，避免处理员不必要的重复劳动，提高速度拾取效率。经过多年的发展，地震速度自动拾取方法较多，根据实现原理，可分为优化搜寻与基于机器学习2种，优化搜寻方法认为幅值最大的位置是能量团中心，采取最优化搜索的方法去找最大幅值，然后获取对应的时间与速度信息。Toldi^[6] 利用共轭梯度算法寻找速度谱中的最大幅值，并获得了对应位置的时间-速度对。张建彬等^[7]将速度谱当作二维矩阵，利用非线性函数寻优法搜索能量团中心，进而获得地震速度。优化搜寻法通常要建立目标函数，利用给定的初始值、约束条件搜索寻优，因此，拾取精度对初始值和约束条件的依赖程度较高。机器学习分为有监督与无监督2种算法，近年来有监督算法在油气勘探领域逐渐兴起^[8-9]，通过标记过的训练样本，将输入的速度谱和输出的时间-速度对建立映射关系，实现速度谱的拾取，主要方法包含YOLO（You Only Look Once）网络与长短期记忆模型相结合^[10]、回归神经网络^[11]、卷积神经网络（CNN）^[12-13]、FOCS神经网络^[14]、多信息融合网络^[15]、深度神经网络^[16-17]、Faster-RCNN网络^[18]等。这些研究提高了计算效率和拾取精度，然而，基于机器学习的有监督算法需要制作大量标签、设计和训练网络，且泛化能力差，整体上计算效率可能不高。无监督学习通常把分析对象分组，发掘数据内在结构，也叫做聚类分析，省去了网络设计与训练的过程，具有算法简单、计算效率高的优势^[19-20]。K均值聚类方法是聚类分析的一种，Galvis等^[21]使用K均值聚类算法实现了面波的自动检测与分类。王迪等^[2]采用常规K均值聚类算法，通过先验速度约束和自适应阈值剔除能量团边缘的速度点，加快了计算速度。Wang等^[22]以速度趋势带作为初始时间条件并应用于Chan-Vese模型，对速度谱进行分割，得到速度候选区域，并采用mean-shift聚类方法进行聚类。然而，常规K均值聚类算法也存在着缺陷，一是速度拾取前要给K和聚类中心赋初值^[23]；二是初值会影响迭代次数和拾取结果；三是最终的速度值与能量团的形状及噪声能量有关，与真实速度有一定的误差。Chen^[24]采取从下往上的拾取方式，部分解决了提前给出K值的难题。Xie等^[25]利用先验速度与3个阈值产生矩形信号，不仅得到了K值，还得到了聚类中心的初始值，初始值与最终拾取的结果较为接近，计算效率和精度都得到较大的提高。周竹生等^[26]在常规K均值聚类算法中引入了权值，提高了速度拾取的精度和抗干扰能力，但该方法需要先验速度，增加了处理员的工作量，且使用幅值阈值剔除速度点，导致弱幅值能量团的拾取能力较低。

基于应用速度点的幅值，改进常规K均值聚类算法，提出一种加权K均值聚类方法，通过设置幅值阈值和比例系数确定剔除不必要的速度点，在时间方向上开时窗，获取幅值最大的若干个速度点，并计算这些速度点的平均速度与平均时间，完成初始值的计算；根据初值剔除或者合并部分聚类中心，迭代时逐步剔除每个中心边缘的速度点；拾取完成后，检查拾取结果，若存在速度反转，则利用前2个拾取结果计算速度反转的相对大小，进而实现异常速度的判断与清除。将该方法与常规K均值聚类方法在模型和实际工区进行测试，以期提高拾取精度的同时，提高计算效率。

1 方法原理

加权K均值聚类方法使用权值改进常规K均值聚类算法以实现全自动速度拾取，无需先验速度信息，根据速度谱的特征，可以自动计算候选区域、K值和聚类中心的初值。

1.1 候选区域的确定

实际生产中，一般采用相关速度谱进行速度拾取，相关速度谱的公式如下^[27]：

$ S=\frac{\sum\limits_{j=-\lambda / 2}^{\lambda / 2}\left[\sum\limits_{i=1}^{N} u\left(t_{i}+j, x_{i}\right)\right]^{2}}{N \sum\limits_{j=-\lambda / 2}^{\lambda / 2} \sum\limits_{i=1}^{N} u\left(t_{i}+j, x_{i}\right)^{2}} $

(1)

$ t_{i}=\sqrt{t_{0}^{2}+\frac{x_{i}^{2}}{v^{2}}} $

(2)

式中：$S$ 为相关值；$u\left(t_{i}+j, x_{i}\right)$ 为振幅；$x_{i}$ 为第$i$ 道数据的炮检距，$\mathrm{m}; ~ \lambda$ 为数据时窗的大小；$N$ 为地震道数；$j$ 为时间方向的顺序号；$v$ 为扫描速度，$\mathrm{m} / \mathrm{s}; t_{0}$ 为双程旅行时，$\mathrm{s}; ~ t_{i}$ 为旅行时间，s。

式（1）以t_i为中心，在时间上通过变量j选取时窗±λ/2内的数据参与计算，当扫描速度与真实速度相等时，计算的相关值最大^[27]。相关速度谱中，具有0 ≤ S ≤ 1的特征，选定的幅值阈值一般能够适应较多的工区。

根据速度谱中能量团中心及其附近区域的幅值更大，而其他区域的幅值较小的特点，剔除部分弱幅值速度点。选定幅值阈值thre1，只允许幅值高于该阈值的速度点参与计算，从而剔除能量团边缘的速度点及其他弱幅值速度点。

$ S_1(v, t)=\left\{\begin{array}{cc} S(v, t), & S(v, t) \geqslant \text { thre } 1 \\ 0, & \text { else } \end{array}\right. $

(3)

式中：t为时间，s。

幅值阈值thre1的选取较为重要，通常为能量团边缘的数值。在实际速度谱中，可能存在噪声或速度谱收敛性差等情况，在某些特殊情况下，如二维地震资料的速度谱中，当噪声的能量比一次波强，或者速度谱的能量团几乎没有规律时，应当采取人工的方式拾取速度。在信噪比相对较高的情况下，一次波的能量团包含的速度点通常数量较多且幅值较大，为了进一步剔除噪声及不收敛的离群速度点，以某幅值大于零的速度点（v_m，t_m）为中心开二维数据窗口，设定二维数据时窗在速度和时间方向的半窗口长度，则该二维数据时窗内的最大速度点数为

$ { num }_{\text {all }}=\left(2 v_{\text {win }}+1\right)\left(2 t_{\text {win }}+1\right) $

(4)

式中：v_win为二维数据时窗在速度方向的半窗口长度，m/s；t_win为二维数据时窗在时间方向的半窗口长度，s。

二维数据时窗的速度区间为（v_m - v_win，v_m + v_win），时间区间为（t_m - t_win，t_m + t_win），若速度点数为num_win，且num_win < thre2 · num_all（thre2为设置的比例阈值，满足0 < thre2 < 1），则将速度点（v_m，t_m）的幅值置0，通过二维时窗内统计速度点个数的方法，能够剔除能量团边缘等离群速度点，提高速度拾取的计算效率与精度。

1.2 初始值的获取

根据一次波的能量团包含的速度点通常数量较多且幅值较大的特点，计算初始值时从浅至深划定时窗，若每个时窗宽度为λ_win，则聚类中心的初始个数K₀为

$ K_{0}=\frac{\left(t_{\text {end }}-t_{\text {first }}\right)}{\lambda_{\text {win }}} $

(5)

式中：$t_{\text {first }}$ 和$t_{\text {end }}$ 分别为速度拾取的起始、截止采样点。

速度点在每个时窗内的集合为$X=\left\{x_{i, j}\right\}, x_{i, j}$ 表示$S_{1}(v, t)$的速度点，下标i与j分别表示行与列。初始值计算中应尽量避免噪声和能量团不规则的不利影响，在集合X内，按照速度点的幅值大小，选取num_max个最大幅值的速度点，在信噪比不低的情况下，位于能量团边缘的噪声速度点幅值较低，而大幅值速度点在能量团中心的可能性很大，且受噪声和能量团不规则的影响最小，将幅值最大的若干个速度点的平均速度与平均时间作为聚类中心$m_{k}^{0}$ 的初值。这种方法获得的初值通常与真实速度较为接近，速度拾取时迭代次数也大幅减少，聚类中心的初始值的计算式为

$ m_{k}^{0}=\frac{1}{n u m_{\max }} \sum\limits_{i=1}^{n u m_{\max }} x_{i, j} $

(6)

式中：$m_{k}^{0}$的上标为迭代次数，0为第0次迭代，即初始值；k为顺序号。

由于式（4）剔除了较多速度点，因此根据式（6）计算的初始值可能存在数值为0或者相邻聚类中心较近的情况，需要剔除或者合并。具体做法是剔除速度为0的聚类中心，若2个或者多个中心的时间距离小于λ_win/2，则将之合并为一个，然后对每个聚类中心采用式（6）计算初始时间和初始速度。

1.3 加权K均值聚类的全自动速度拾取

对初始聚类中心进行剔除与合并之后，能让聚类中心和能量团的对应关系更加合理，若聚类中心最终保留K个，每个速度点到这K个聚类中心的距离分别为

$ d_{k}=\left\|x_{i, j}-m_{k}\right\|^{2} $

(7)

依据距离最小的规则进行聚类或分组，对一个聚类而言，每个速度点均距其中心最近。聚类后，每个聚类有若干个速度点，这些速度点的中心可能与初始中心不一致，需获得新的中心：

$ m_{k}^{l}=\frac{1}{M_{k}} \sum\limits_{X \in c_{k}} x_{i, j} $

(8)

式中：l为迭代次数；c_k和M_k分别第k个聚类的速度点集合和速度点的个数。

由式（8）可知，计算聚类中心时所有的速度点一视同仁，没有能量团中心与边缘的分别，因此得到的聚类中心受限于该聚类中速度点组成的几何形状，拾取的是几何形状的中心。若能量团具有规则的对称形状，这种做法拾取的速度较为精确，而当能量团不规则或者某个时窗内含有多次波等噪声的速度点时，拾取速度的精度将会降低甚至出现较大误差。因此，对不同的速度点赋予不同的权值^[28]，使聚类中心逼近能量团中心，权值的计算式为

$ W_{i, j}=k_{\mathrm{a}} A_{i, j}^{n} $

(9)

式中：k_a为正实数；A_i，j为速度点的幅值；n是幂次方，通常不小于1。由式（9）可知，速度点的幅值和对应的权值成正比，因此聚类中心采用如下公式更新：

$ m_{k}^{l}=\frac{1}{N_{k}} \sum\limits_{X \in c_{k}} x_{i, j} W_{i, j} $

(10)

式中：$N_{k}=\sum\limits_{X \in c_{k}} W_{i, j}$。一般情况下，速度点距离能量团中心越近，其幅值越大。由于权值和幅值呈正相关，采用式（10）所示的聚类中心更新方法，得到的结果更加逼近能量团中心。

为了避免无限次迭代，需建立终止迭代的目标函数及设置最大迭代次数，目标函数如下：

$ J^{l}=\sum\limits_{k=1}^{K} \sum\limits_{X \in c_{k}}\left\|x_{i, j}-m_{k}^{l}\right\|^{2} $

(11)

若$\left|J^{l}-J^{l-1}\right| \leqslant$ error则终止迭代，其中error为很小的正数（如0.01）。为了加快计算速度、避免不规则能量团对拾取精度的干扰，在更新速度聚类中心之后，需剔除每个聚类中最边缘的小部分速度点。因此，随着迭代的进行，速度点的数量越来越少且更加逼近能量团中心，能量团的形状也将越来越规则，在提高计算效率的同时，也提高了速度拾取的精度。

1.4 异常速度的判断与剔除

在实际地震速度拾取过程中，由于噪声等影响，拾取结果可能存在异常，为了保证拾取结果的准确性，应对拾取结果进行处理。首先，剔除与整体速度轨迹偏离较大的速度点，通过对拾取结果进行平滑，比较拾取结果与平滑结果的速度差异，若拾取速度大于给定的阈值，则剔除该速度点。其次，剔除可能存在的多次波，在速度谱中，多次波与速度反转有一定的相似性，但速度反转时，速度减小的范围有限^{[2, 29]}，而多次波速度减小的数值较大，根据这一特点可以进行多次波剔除。

若第k个聚类中心的速度v_k小于第k - 1个聚类中心的速度v_{k - 1}，先计算与相邻拾取结果的斜率值：

$ q_{k, k-1} =\frac{v_{k}-v_{k-1}}{t_{k}-t_{k-1}} $

(12)

$ q_{k-1, k-2} =\frac{v_{k-1}-v_{k-2}}{t_{k-1}-t_{k-2}} $

(13)

式中：t_k为第k个聚类中心的时间，s。

2条直线的斜率与夹角（θ）存在如下关系：

$ \tan \theta=\frac{q_{k, k-1}-q_{k-1, k-2}}{1+q_{k, k-1} q_{k-1, k-2}} $

(14)

得出夹角的计算式：

$ \theta=\arctan \left(\frac{q_{k, k-1}-q_{k-1, k-2}}{1+q_{k, k-1} q_{k-1, k-2}}\right) $

(15)

按照式（12）—（15），若计算出的夹角偏小，则当作正常的速度反转而保留；反之，则为多次波，应该剔除（图 1）。

2 模型数据测试

模型数据测试流程为

（1）建立包含5个层介质的水平层状速度模型。模型从上至下5个地层厚度依次为900 m、500 m、10 m、550 m和540 m，速度分别为2 900 m/s、3 700 m/s、2 800 m/s、3 800 m/s和4 100 m/s；其中第3层为薄层，每个界面的反射系数较大，能够形成多次波（图 2）。

（2）对速度模型正演。获得含3个一次波同相轴和多个多次波同相轴的CMP道集，该道集的道距10 m，最小炮检距为0，总道数为300（图 3）。

（3）计算CMP道集的速度谱。可得到含3个一次波能量团（如图 4a中字母P所示）和5个多次波能量团的相关速度谱，多次波能量团中，3个幅值较强，其余2个较弱。

（4）剔除弱幅值速度点。选定速度谱幅值阈值，剔除时间为2.0~2.5 s时的弱能量团以及0.5~ 2.0 s能量团外侧的弱幅值速度点，得到仅保留幅值高于阈值的速度点的速度谱（图 4b）。

（5）剔除离群速度点。对剔除了弱幅值速度点的速度谱，以幅值大于0的速度点为中心开二维数据窗口，统计二维数据窗口内幅值大于0的速度点，若统计的速度点数小于一定的比例，则剔除该速度点。以该方法对全部速度点进行处理，将离群速度点剔除，得到的速度谱中能量团边缘的部分速度点被剔除，能提高计算效率和拾取精度（图 4c）。

（6）根据式（5）和式（6）计算初始值。将大幅值速度点的平均速度与平均时间作为聚类中心的初值，因在1.2 s附近的一次波能量团和多次波能量团在时间方向上无法区分，将这2个能量团当成一个处理，即速度谱中只有5个聚类中心。

（7）自动速度拾取。分别采用常规、加权K均值聚类方法对速度谱自动拾取，并进行效果对比。拾取结果（表 1）显示，2种自动拾取方法的结果相差不大，本文方法提取的5个聚类中心速度整体略高，在1.2 s附近，常规K均值聚类法的拾取速度明显偏小，分析认为是该时窗包含多次波的速度点，常规方法将多次波速度点当作一次波，且迭代过程中无法剔除导致。

（8）根据计算夹角[式（15）]，进行异常速度初步判断与剔除。本文方法利用3个相邻拾取结果计算夹角，最终剔除了第4、第5个拾取结果。在速度谱中展示3个拾取结果和线性插值的速度拟合曲线，结果（图 5a）显示，拾取结果较常规K均值聚类方法的精度更高，常规K均值聚类方法在1.2 s附近聚类能量团的几何形状向低速移动，拾取的速度低于真实值。利用2种自动拾取方法获得的3组时间-速度对，结合通过内插和外推手段获得其他时间的速度，对得到的速度对CMP道集进行动校正，2种方法拾取的第3个能量团的速度虽然相差60 m/s，但相对误差仅为1.7%，因此，同相轴动校正后的差异也很小（图 5b，5c），2种方法得到的动校正道集相似度很高。

综合分析可知，在模型测试时，加权K均值聚类方法虽然具有权值的优势，且迭代过程中能够剔除速度点，但模型数据速度点与聚类中心均较少，加上能量团的形状较为规则，其计算效率仍然与常规K均值聚类方法相当。

3 实际应用效果

以塔里木盆地TD区块为例，选取了一条含750个道集的测线测试加权K均值聚类方法自动拾取的效果。图 6a为该测线信号反动校正后的CRP道集，共65道，道间距为100 m，最小和最大炮检距分别为100 m和6 500 m，记录时长为7.0 s，计算得到其速度谱如图 6b所示，剔除其中弱幅值速度点和离群速度点，最终得到的速度谱中能量团更加聚焦，且速度点分布更加规则（图 7）。

分别采用常规K均值聚类、加权K均值聚类2种方法对速度谱进行自动速度拾取。结果（图 8）显示：当时间为2.0~4.5 s时，能量团的收敛性较好，2种自动拾取方法的结果基本一致；当时间为4.5~6.0 s时，能量团更加发散，常规K均值聚类方法拾取结果与真实速度存在少量偏差，加权K均值聚类方法拾取结果精度更高。2种方法计算的速度差异不明显，动校正结果的差异非常小；常规K均值聚类每个道集迭代了44次，而加权K均值聚类每个道集仅迭代6次，计算效率约提升了7倍。

对整条测线，人工拾取是每间隔50个道集进行拾取，未拾取的道集采用插值方法得到结果，因此速度场较为平滑（图 9a）。常规、加权K均值聚类方法是自动对全部地震道进行拾取，得到的速度场（图 9b，9c）基本相同，精度比人工拾取明显更高，在距离为3 km，2.5~3.0 s时存在明显的速度变化。

在自动拾取速度过程中可能出现少数聚类中心受到噪声干扰等异常，为了消除这一影响，从时间与空间2个维度对自动拾取的速度场进行平滑。分别采用人工拾取速度与2种自动拾取速度对整条测线的CRP道集进行动校正，然后采用外切方法将浅层远道出现拉伸畸变的数据切除，最后进行叠加得到叠加剖面。结果（图 10）显示：3个叠加剖面的同相轴倾角、断层位置等基本相同，但常规、加权K均值聚类方法得到的叠加图精度更高，表现为断裂更加清晰（如图 10中红色椭圆框所示）；叠加能量更强，能清晰地看到背斜构造，且波组特征更加合理（如图 10中的黑色方框所示）。

4 结论

（1）加权K均值聚类法通过剔除弱值和离群速度点，提高了计算效率和速度拾取的精度；在某个时窗内选取若干个幅值最大的速度点，以其均值作为初值，与最终的拾取结果较接近，完成拾取所需的迭代次数也大幅减少，加快了计算速度。

（2）加权K均值聚类法对不同的速度点赋予不同的权值，使聚类中心逼近能量团中心，能量团更加规则，相较于常规K均值聚类法，计算效率和拾取精度均更高，在模型测试及塔里木盆地TD区块的实际资料的应用结果显示，其计算效率约提升了7倍，精度提高了1.7%。

（3）在某些特殊情况下，当噪声的能量比一次波强，或者速度谱的能量团没有规律时，无法使用加权K均值聚类法进行自动速度拾取。