基于自適應(yīng)閾值約束的無監(jiān)督聚類智能速度拾取

王迪， 袁三一*， 袁煥， 曾華會， 王尚旭

1 中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室， 北京 102249 2 中國石油天然氣股份有限公司勘探開發(fā)研究院西北分院， 蘭州 730000

0 引言

速度分析是常規(guī)地震資料處理中的重要環(huán)節(jié)之一，主要包括疊加速度分析、偏移速度分析和速度層析反演等幾種主要方法.其中，疊加速度作為速度分析中的重要內(nèi)容，不但是地震資料常規(guī)疊加處理和時間偏移的必要過程，也是層析反演、疊前偏移、阻抗反演和全波形反演等初始速度模型建立的基礎(chǔ)(杜啟振等，2011；Zhang et al.，2015；Du et al.，2017；徐文君等，2017；Yuan et al.，2019).為了獲得疊加速度場，需要在每個共中心點(CMP)道集對應(yīng)的速度譜上進行拾取，現(xiàn)在普遍應(yīng)用的方法是人工拾取.但是，人工拾取需要花費大量的時間，成本高昂，拾取效率較低，且拾取精度易受到人為經(jīng)驗的影響.這些問題在寬方位、高密度地震資料處理時顯得更為突出.因此，在保證速度拾取精度的前提下，有必要研究一種快速且有效的自動拾取速度譜的方法，以減少地震資料處理人員的人工操作，提高速度拾取效率.

隨著地震勘探程度的不斷深入和技術(shù)條件的不斷改善，自動拾取速度的研究有了一定的進展.目前發(fā)展的方法大部分可以分為優(yōu)化搜尋和機器學(xué)習(xí)這兩類.優(yōu)化搜尋速度自動拾取方法本質(zhì)上是一種反演類方法，主要是利用最優(yōu)化算法和最大相似度量準則，設(shè)定一定的速度約束條件，對初始速度模型加以擾動，自動尋找速度譜中疊加能量的最優(yōu)解.在自動拾取發(fā)展的早期，優(yōu)化搜尋類方法發(fā)展迅速.Toldi(1989)最先提出了一種基于速度譜疊加能量的自動速度拾取方法，以最大化疊加能量和作為目標函數(shù)，反演得到速度場，但該方法假設(shè)模型線性化，也很難很好處理信噪比低的地區(qū)，需要進一步的約束.Zhang和Claerbout(1990)將非線性優(yōu)化方法應(yīng)用于速度拾取中；Zhang(1991)后續(xù)又對此方法做了進一步的優(yōu)化.Lumley(1997)提出了一種利用蒙特卡洛隨機擾動方法進行自動層速度的拾取.Meshbey等(2002)和Fomel(2009)將速度拾取視為一個變分問題，Almarzoug和Ahmed(2012)在其基礎(chǔ)上發(fā)展了自動地震速度拾取方法.這些方法在獲取相對簡單地質(zhì)體的速度時已經(jīng)取得了一定的效果，但大部分方法依賴于初始層速度模型，并對于橫向速度突變的地區(qū)，拾取準確性相對不高.還有些方法從振幅隨偏移距變化(AVO)分析中挖掘速度信息，Swan(2001)在疊前道集中通過AVO分析提取特定屬性特征，構(gòu)建剩余速度指示器，從而迭代更新速度場.由于其方法的穩(wěn)定性欠佳，適用范圍有限，Ratcliffe和Roberts(2003)在此基礎(chǔ)上增添了監(jiān)控收斂準則，從而更能適應(yīng)實際生產(chǎn)要求.基于疊前AVO分析的這類方法需要部分道準確速度作為初始模型，只能部分解決速度自動拾取問題.

近年來，由于計算機性能的顯著提高，機器學(xué)習(xí)也被應(yīng)用于速度自動拾取，主要有深度學(xué)習(xí)和聚類分析方法.根據(jù)不同特征的神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)能力，Ma等(2018)及Park和Sacchi(2020)分別提出了基于卷積神經(jīng)網(wǎng)絡(luò)的疊加速度拾取方法；Biswas等(2018)及Fabien-Ouellet和Sarkar(2020)發(fā)展了基于循環(huán)神經(jīng)網(wǎng)絡(luò)的速度拾取方法；Zhang等(2019)利用長短期記憶網(wǎng)絡(luò)(LSTM)算法提出了一種深度學(xué)習(xí)拾取速度方法.現(xiàn)有的深度學(xué)習(xí)拾取速度方法主要思路是通過學(xué)習(xí)大數(shù)據(jù)，用人工神經(jīng)網(wǎng)絡(luò)來建立CMP道集或速度譜數(shù)據(jù)(輸入)和速度模型(輸出)之間的非線性映射關(guān)系，模擬CMP道集或速度譜數(shù)據(jù)到速度的逆過程.在部分地區(qū)中能獲得較高精度的速度場，但要學(xué)習(xí)大量且豐富的樣本標簽數(shù)據(jù)，對預(yù)測的輸入數(shù)據(jù)質(zhì)量和采樣間隔等要求較高，且其預(yù)測過程和結(jié)果目前仍很難可解釋.

聚類分析是無監(jiān)督機器學(xué)習(xí)算法中的一種，可以挖掘數(shù)據(jù)中的隱含結(jié)構(gòu)，對具有相同特征的數(shù)據(jù)進行分組，從而得到數(shù)據(jù)的大致分布.相比于深度學(xué)習(xí)方法，聚類算法不需要標簽及訓(xùn)練網(wǎng)絡(luò)，算法簡潔，可解釋性強.Smith(2017)利用多種屬性特征劃分出速度范圍，并首次將聚類方法應(yīng)用于速度譜自動拾取.K均值聚類算法是聚類分析方法中一種基本的劃分式方法(MacQueen，1967；Hung et al.，2013；Ahmad and Hashmi，2016)，該算法簡便易實現(xiàn)，且計算速度較快，通常被作為大樣本聚類分析的首選算法，近年來也常被用于地震勘探中(Galvis et al.，2017；Araya-Polo et al.，2017).Song等(2018)提出了一種固定K值的智能速度拾取方法，以鄰近道人工拾取速度作為篩選聚類范圍的約束，并應(yīng)用K均值算法進行聚類，能減少人工挑選的時間，從而降低處理數(shù)據(jù)集的成本.針對難以選取固定K值的問題，密度聚類等其他無監(jiān)督方法也被應(yīng)用到自動拾取速度方法中(Bin Waheed et al.，2019；Zhang and Lu，2016).而Chen和Schuster(2018)則提出了一種自底向上的策略來實現(xiàn)K均值迭代過程，一定程度上解決了不同CMP道集選擇不同K值的問題.但仍有兩個問題需要解決：一是對于不同深度反射能量差異大的道集，如煤系地層沉積時，往往造成煤層和目標儲層存在上下疊置關(guān)系，從而導(dǎo)致有效砂體反射信息淹沒于煤層反射(田忠斌等，2016；李勤等，2019)，容易會忽視弱反射的速度拾取；二是聚類范圍的篩選至關(guān)重要，經(jīng)典的K均值聚類算法容易受到噪聲和離群點的影響，使得聚類中心偏離主要數(shù)據(jù)區(qū)域，從而使得速度拾取精度受到影響.

因此，本文提出了一種基于自適應(yīng)閾值約束的無監(jiān)督聚類智能速度拾取方法，減少人工勞動，提高速度拾取效率并確保其拾取精度.根據(jù)速度譜中不同深度的能量團聚焦特征，采用時窗方法在速度譜中計算自適應(yīng)閾值，從而識別出一次反射波速度能量團作為速度拾取的候選區(qū)域.然后，根據(jù)K均值算法將不同的速度能量團聚類，并將聚類中心作為拾取的疊加速度，實現(xiàn)速度的自動拾取.同時，根據(jù)人工拾取經(jīng)驗設(shè)計了離群速度處理環(huán)節(jié)，以獲得更為光滑的速度場.最后，數(shù)值模型和實際地震數(shù)據(jù)驗證了該方法的可行性，并獲得了較高的拾取精度.

1 方法原理

首先，我們簡單介紹地震資料處理中人工拾取疊加速度的基本過程：根據(jù)CMP道集計算速度譜，識別出不同時間位置的速度能量團，拾取能量團中心點(中淺層的速度能量團一般為峰值點)，對拾取的時間速度對進行插值即得到疊加速度曲線.不難發(fā)現(xiàn)識別速度能量團和拾取中心點是從速度譜中拾取疊加速度的關(guān)鍵.因此，無監(jiān)督聚類智能速度拾取的思路主要是圍繞這兩項內(nèi)容類比并模擬人工速度拾取過程.一是在速度譜中，根據(jù)一定的規(guī)則識別出一次反射速度能量團作為速度拾取的候選區(qū)域，二是應(yīng)用聚類算法對此區(qū)域樣本進行聚類，并拾取中心點.為了避免噪聲等因素的影響，根據(jù)疊加速度的基本變化規(guī)律和人工拾取經(jīng)驗，加入了離群速度處理環(huán)節(jié).

1.1 自適應(yīng)確定速度拾取候選區(qū)

任何地下地層之間存在的波阻抗差異均會引起地震反射，而不同巖性的地層界面引起的反射能量可能會存在較大差異(馬水平等，2010)，導(dǎo)致在同一CMP道集中強弱反射相間存在，進而不同深度的速度能量團能量有較大的變化.而目前的聚類速度自動拾取方法(Chen and Schuster，2018)是在整個速度譜的基礎(chǔ)上，設(shè)定統(tǒng)一的閾值確定速度拾取候選區(qū)，固定的閾值使得其拾取精度不僅一定程度上受到速度能量團聚焦性的制約，且拾取時容易忽視弱反射速度能量團.為了盡可能識別到所有有效的一次反射速度能量團，同時保留強反射和弱反射，我們設(shè)計了一種滑動窗方法對速度譜進行處理.先定義速度譜矩陣S(從上至下對應(yīng)時間增加，而從左至右對應(yīng)速度增大)，即

(1)

(2)

(3)

(4)

根據(jù)反射波雙曲線的幾何特征，速度和深度是影響其雙曲線形態(tài)的主要因素.常規(guī)速度譜的生成本質(zhì)上是一個應(yīng)用不同時間速度對掃描進行動校正并計算疊加能量的過程(Neidell and Taner, 1971).然而不同深度的雙曲線走時對速度變化敏感度不一，在相同的速度變化中位于中淺層的雙曲線走時變化幅度更劇烈，從而導(dǎo)致其不同時窗間的疊加能量差異大，最終計算的該位置速度能量團聚焦性高，即速度譜中淺層對速度變化敏感.因此在深度和速度的共同影響下，實際數(shù)據(jù)中速度能量團聚焦性逐漸下降，橫向?qū)挾仍龃?，速度分辨率降?根據(jù)這一特點，實際應(yīng)用中淺層速度的拾取精度要求高，而深層速度能量團聚焦性差，拾取精度要求相對較低.我們結(jié)合這一特征設(shè)計閾值百分比p(i)為

p(i)=pmin+(pmax-pmin)e-i，

(5)

其中，pmin為最小百分比閾值，pmax為最大百分比閾值.pmin和pmax的取值遵循的基本原則是p=pmax時需要能篩選出較為聚焦的淺層能量團，而p=pmin需要能篩選出有一定寬度且有遞增趨勢的深層能量團，具體取值可以根據(jù)速度譜能量團的分布試驗而定.從上式可以看出，p(i)與深度呈負相關(guān)，深度增加，閾值百分比減小，從而降低了閾值門檻，擴大了深層速度拾取候選區(qū)范圍，這基本符合實際速度譜中能量團的變化特征.同時，在淺層動校正對速度誤差敏感，較小的閾值百分比p可以提高淺層速度拾取精度；而在深層人工拾取的速度并不一定在能量團峰值處，較大的閾值百分比p擴大了拾取候選區(qū)范圍，更符合人工拾取經(jīng)驗.

(6)

其中，xij表示大于閾值時，速度譜中能量團sij對應(yīng)的時間速度對的位置坐標.經(jīng)過閾值門檻的篩選，將速度譜矩陣降維，僅保留非零時間速度對xij作為速度拾取候選區(qū).

1.2 K均值聚類速度拾取

(7)

其中，dk表示歐氏距離.按照距離最小原則進行劃分，將數(shù)據(jù)點劃分到與其距離最近的聚類中心所在的類.所有數(shù)據(jù)點分類后，重新計算k個聚類中每個聚類的全部數(shù)據(jù)點的平均值，該平均值所在的數(shù)據(jù)點成為新的聚類中心.其中，新的聚類中心可寫為

(8)

其中，ck為第k類的所有數(shù)據(jù)點的合集，Nk為第k類中的總數(shù)據(jù)點數(shù)，上標l=1,2,…,L為迭代次數(shù).在確定新的聚類中心后，計算每個類別中數(shù)據(jù)點與聚類中心的距離和，其目標函數(shù)J為

(9)

然而在實際數(shù)據(jù)中，不同疊前道集中的速度譜所需拾取的能量團個數(shù)不一致，針對這一問題，我們選取一種自底向上的迭代策略來求取目標函數(shù)最佳解.即先給定一個初始類別數(shù)K0，按照初始聚類中心的分布，將數(shù)據(jù)點初步分為K0類，并逐步迭代找到聚類類別最少的聚類分布，聚類中心個數(shù)在迭代過程中逐步減小.其具體迭代過程如下：

(2)由式(7)計算距離，按照最小距離原則將數(shù)據(jù)點進行重新分類.對于未分配到數(shù)據(jù)點xij的類別ck，剔除所在類別對應(yīng)的聚類中心，不參與下次迭代，同時總聚類數(shù)相應(yīng)減?。?/p>

(3)根據(jù)式(8)逐個類別重新計算新的聚類中心；

(4)重復(fù)步驟(2)和步驟(3)，直到目標函數(shù)滿足收斂要求.

1.3 離群速度處理

K均值聚類算法中每個聚類的中心都由每個聚類里所有數(shù)據(jù)點求均值得到.當存在一定噪聲(如多次波、野值等)時，聚類中心易受干擾，偏離速度能量團峰值，影響聚類效果.因此，為了獲得更接近人工拾取的速度場，在實際應(yīng)用時我們對離群速度進行了后處理.人工拾取速度經(jīng)驗指出，一般情況下疊加速度會隨著時間的增加而增大，即使在速度倒轉(zhuǎn)的情況下，速度倒轉(zhuǎn)也有一定范圍(Han et al.，2007；Yan et al.，2016；Sun and Zhang，2020).基于速度譜的以上兩個特征，判斷速度是否倒轉(zhuǎn)的準則qk如下：

(10)

其中，ymk為拾取速度mk處的位置縱坐標，xmk為拾取速度mk處的位置橫坐標，將拾取的速度點與鄰近點相比較，判斷出是否為速度倒轉(zhuǎn)點.若mk為速度倒轉(zhuǎn)點且速度倒轉(zhuǎn)變化范圍大于判斷值，則去除速度異常點.

2 例子

2.1 一維數(shù)值模型測試

為了測試提出方法的可行性，首先使用一維數(shù)值模型進行測試.如圖1a為一個預(yù)定義速度模型創(chuàng)建的CMP道集.該道集由主頻為30 Hz的零相位雷克子波經(jīng)射線追蹤方法數(shù)值模擬獲得，共包含了1300個時間采樣點和31道，采樣間隔為1 ms，道間距為50 m.為了測試本方法對弱反射的速度拾取能力，模型中設(shè)計了3個弱反射同相軸(0.6 s、0.75 s和1.2 s)和4個強反射同相軸相間分布.在 0.9 s及0.91 s處有兩個位置相近的同相軸，用于測試該方法在子波干涉時的速度拾取準確性.對圖1a中的CMP道集進行速度分析，得到速度譜(圖1b).從圖1b可以看出，弱反射的速度能量團能量遠小于其他同相軸，而兩個干涉同相軸(0.9 s和0.91 s)的速度能量團混疊在了一起.此外，由于不同深度速度敏感性不同，相比于淺層的能量團，深層的能量團速度方向分辨率更低.

圖1 一維數(shù)值模型(a) CMP道集； (b) 速度譜.Fig.1 One dimensional model(a) CMP gather； (b) Velocity spectrum.

應(yīng)用自適應(yīng)閾值和固定閾值的無監(jiān)督聚類方法(Chen and Schuster，2018)進行智能速度拾取，并將拾取的速度進行動校正.根據(jù)速度譜能量團的大致分布，兩種方法均選定初始類別K0為12，由于速度譜中能量團分布較為稀疏，在自適應(yīng)閾值方法中選取較大的滑動窗窗長dwin為100 ms，而為了進一步說明常規(guī)聚類方法的局限性，比較了兩種不同固定閾值(p=90%和p=95%)對聚類效果的影響.圖2(a—f)為常規(guī)聚類方法中兩種閾值情況下的速度拾取候選區(qū)、速度拾取結(jié)果和動校正后的道集.在常規(guī)聚類方法中，為了拾取到弱反射速度能量團，先測試了較小的固定閾值百分比(p=90%)進行速度拾取，較小的p其對應(yīng)閾值偏小，增大了深層能量團樣本量(圖2a中紅色部分)，從而增加了后續(xù)智能速度拾取的運算量.同時，K均值算法易于將密集的、球狀或團狀的數(shù)據(jù)聚類，篩選的速度拾取候選區(qū)也需要符合這一特征，而p偏小時篩選的深層候選區(qū)范圍并不聚焦，由此易拾取到多個聚類中心，顯著降低深層速度拾取準確性(圖2b)，經(jīng)動校正后0.9 s及0.91 s處同相軸顯著向上彎曲.當選取較大的固定閾值百分比(p=95%)時，較大的p其對應(yīng)閾值偏大，只能兼顧能量較強的速度能量團，篩選速度拾取候選區(qū)時忽視了弱反射能量團(圖2d)，淺層的2個弱反射和深層的1個弱反射都未能拾取到速度(圖2e)，不利于速度曲線的精細刻畫，也降低了拾取精度.從而0.8s以上的淺層同相軸明顯校正不足，同相軸向下彎曲嚴重(圖2f).圖2(g—i)為自適應(yīng)閾值方法的速度拾取及動校正結(jié)果，自適應(yīng)閾值在識別能量團時，強度不同的速度能量團均能被識別(圖2g)，所有反射的速度能量團均被包含在候選區(qū)內(nèi).因此，本文方法仍能拾取到正確速度(圖2h)，即使在0.9 s處左右兩個干涉同相軸情況下，且動校正后同相軸基本被拉平(圖2i).需要注意的是，由于K均值算法其劃分準則是歐氏距離，本質(zhì)上更適用于密集的、球狀或團狀的數(shù)據(jù)聚類，而速度能量團并不是標準的圓，處于深層的能量團形狀更接近于橢圓，因此同一能量團可能被分為兩類.而在0.9 s處左右兩個干涉同相軸的速度能量團在此情況下，綜合初始模型影響拾取到兩個聚類中心.

圖2 不同方法拾取速度比較(a) 常規(guī)聚類方法確定的速度拾取候選區(qū)(p=90%)； (b) 在(a)中拾取的疊加速度； (c) 應(yīng)用(b)中速度動校正后道集； (d) 常規(guī)聚類方法確定的速度拾取候選區(qū)(p=95%)； (e) 在(d)中拾取的疊加速度； (f) 應(yīng)用(e)中速度動校正后道集； (g) 自適應(yīng)閾值方法確定的速度拾取候選區(qū)； (h) 自適應(yīng)閾值方法拾取的疊加速度； (i) 應(yīng)用(h)中速度動校正后道集.Fig.2 Comparison of velocities picked by different methods(a) Candidate set selected by the conventional method (p=90%)； (b) Stacking velocity curve picked on (a)； (c) NMO gather corrected by the velocity picked in (b)； (d) Candidate set selected by the conventional method (p=95%)； (e) Stacking velocity curve picked on (d)； (f) NMO gather corrected by the velocity picked in (e)； (g) Candidate set selected by the proposed method； (h) Stacking velocity curve picked by the proposed method；(i) NMO gather corrected by the velocity picked in (h).

2.2 二維數(shù)值模型測試

為了進一步證明本文方法復(fù)雜模型中的有效性，我們對二維的Marmousi模型(Martin et al.，2006)進行測試.在Marmousi模型中等間隔抽取了74道作為本次測試的真實層速度場(圖3)，進而由速度模型和Gardner密度公式正演得到一系列CMP道集.根據(jù)速度譜能量團的大致分布規(guī)律(能量團個數(shù)，單個能量團延續(xù)時長等)，兩種方法均選定初始類別K0為15，并通過多次試驗，在自適應(yīng)閾值方法中選取滑動窗窗長dwin為40ms，最小閾值百分比pmin=85%，最大閾值百分比pmax=95%，而在常規(guī)聚類方法選取了固定閾值百分比p=90%.

圖3 層速度場模型Fig.3 Interval velocity model

圖4(a,b)分別展示了Marmousi模型中第10道的CMP道集及常規(guī)聚類方法拾取的疊加速度(綠線)、本文方法拾取的疊加速度(紅線)和真實均方根速度(黑線)的對比.在CMP道集中，淺層和深層反射波能量差異大，速度譜中不同深度的能量團間幅值相差較大，深層反射波能量強，0.5 s以上存在若干弱反射.在此情況下，本文方法拾取的速度仍與真實值基本一致，且相較于常規(guī)聚類方法，在深層(約1.2 s)和弱反射(約0.4 s)處拾取的速度均準確性更高.應(yīng)用兩種方法所拾取的速度場進行動校正，如圖4(c,d)所示.從動校正后的道集中可以看出，本文拾取的速度拉平了所有同相軸，證明了算法的有效性，也驗證了速度拾取的正確性.而常規(guī)方法拾取的速度偏離真實值，其動校正后道集在深層約1.2 s處同相軸出現(xiàn)上翹現(xiàn)象.

圖4 Marmousi模型中的單道CMP拾取結(jié)果(a) CMP 道集； (b) 兩種方法自動拾取的速度曲線； (c) 常規(guī)聚類方法后動校正； (d) 本文方法拾取速度動校正.Fig.4 CMP gathers and the corresponding velocity curves picked by two different methods(a) CMP gather； (b) Automatic estimated stacking velocity curves corresponding to two methods； (c) NMO gather corrected by the conventional automatic estimated stacking velocity； (d) NMO gather corrected by the proposed automatic estimated stacking velocity.

對整個Marmousi模型重復(fù)上述單道的處理流程，得到了兩種不同方法自動拾取的疊加速度場.圖5展示了Marmousi模型真實均方根速度場(圖5a)、常規(guī)聚類方法拾取的速度場(圖5b)和自適應(yīng)閾值方法拾取的速度場(圖5c).總體上看，常規(guī)方法在兩側(cè)楔狀體等構(gòu)造位置拾取精度較低(如藍色箭頭所示)，且速度場橫向上連續(xù)性差，如黑色箭頭指示的CMP道拾取到速度異常值，與相鄰道差異大.而自適應(yīng)閾值方法拾取的速度場在構(gòu)造走勢上基本與真實值一致，在大部分道上比常規(guī)聚類方法的拾取精度更高，對于楔狀體等構(gòu)造細節(jié)刻畫地更為準確，速度場變化較為連續(xù).為了更直觀看到兩種方法速度拾取精度上的差異，我們應(yīng)用圖5中的三個速度場進行疊加，其疊加剖面如圖6所示.與常規(guī)聚類方法(圖6b)相比，本文方法的疊加剖面(圖6c)能更清晰地反映出主要地質(zhì)構(gòu)造，尤其是在速度拾取精度更高的位置(如紅色橢圓內(nèi))，動校正后同相軸被拉平同相疊加，清晰地成像出了在1.3 s處的楔狀體.

圖5 拾取的速度場與真實速度場比較(a) 真實速度場； (b) 常規(guī)聚類方法拾取的速度場； (c) 自適應(yīng)閾值方法拾取的速度場.Fig.5 Comparisons of picked stacking velocity structures and the true reference(a) True stacking velocity structure； (b) Velocity structure picked by the conventional method； (c) Velocity structure picked by the adaptive threshold method.

圖6 Marmousi模型疊加剖面(a) 真實速度場疊加； (b) 常規(guī)聚類方法速度場疊加； (c) 自適應(yīng)閾值方法速度場疊加.Fig.6 Stacked sections of picked stacking velocity structures and the true reference(a) Stacked section obtained by true velocity field； (b) Stacked section obtained by the conventional method； (c) Stacked section obtained by the adaptive threshold method.

2.3 實際數(shù)據(jù)應(yīng)用

最后，采用了一條陸上的二維實際測線對本文方法進行測試，該測線共包含500個CMP道集.每個CMP道集的偏移距及道數(shù)不等，數(shù)據(jù)的時間總長度2.5 s，采樣間隔為1 ms.根據(jù)速度譜能量團的大致分布，兩種方法均選取初始類別K0為12，自適應(yīng)閾值方法中選取滑動窗窗長dwin為300 ms，而在常規(guī)聚類方法選取了固定閾值百分比為90%.由于實際數(shù)據(jù)噪聲較多，為了保證速度模型具有一定光滑性，在速度拾取后對離群速度進行了處理，去除了大于5%的速度倒轉(zhuǎn)異常值.

圖7 實際數(shù)據(jù)CMP道集速度拾取及動校正結(jié)果(a) CMP 100道集； (b) CMP 100道集的速度拾取結(jié)果； (c) 本文方法拾取速度對(a)進行動校正后道集； (d) 常規(guī)聚類方法拾取速度對(a)進行動校正后道集； (e) CMP 225道集； (f) CMP 225道集的速度拾取結(jié)果； (g) d本文方法拾取速度對(e)進行動校正后道集； (h) 常規(guī)聚類方法拾取速度對(e)進行動校正后道集.Fig.7 Velocity picking and NMO results of field data(a) CMP 100 gather； (b) Estimated stacking velocity of CMP 100； (c) NMO result of (a) using velocity picked by the proposed method； (d) NMO result of (a) using velocity picked by the conventional method； (e) CMP 225 gather； (f) Estimated stacking velocity of CMP 225； (g) NMO result of (e) using velocity picked by the proposed method； (h) NMO result of (e) using velocity picked by the conventional method.

圖7展示了CMP 100和225道集的速度拾取結(jié)果及動校正道集.圖7(a—d)展示了在CMP 100道集中，三種方法拾取速度曲線比較以及常規(guī)聚類方法和自適應(yīng)閾值方法的動校正后道集.從疊前CMP道集可以看出，淺層信噪比高，深層信噪比較低，因而本文方法在淺層拾取精度更高.從速度譜可以看出，在 0.5 s處均存在高速速度異常，常規(guī)方法速度拾取范圍包含了部分高速異常值，而最終的聚類的中心是本類中所有數(shù)據(jù)點的均值，受部分高速值影響，導(dǎo)致拾取的聚類中心向異常值靠攏，降低了拾取精度.而本文中的閾值是隨時間深度變化而變化，在同一時間位置能對異常速度能量有一定的抑制作用，避開了高速速度異常.動校正后，本文中的方法(圖7c)在0.6 s左右，反射同相軸被拉平，而常規(guī)聚類方法0.6 s處校正不足，同相軸向下彎曲.圖7(e—h)分別為三種方法對CMP 225道集進行拾取的速度曲線比較，以及常規(guī)聚類方法和自適應(yīng)閾值方法的動校正后道集比較.相比于CMP 100道集，CMP 225道集(圖7e)的反射波能量從淺到深差異小，常規(guī)聚類方法與本文方法在主要反射處拾取的速度精度相近(圖7f).但常規(guī)聚類方法在0.2 s處拾取的速度遠偏離人工拾取值(圖7f)，其對應(yīng)的同相軸也未被拉平(圖7h).而本文方法淺層未拾取到速度異常值，且在弱反射(2s左右)處拾取的速度更接近人工拾取值，基本拉平了同相軸(圖7g).

在構(gòu)建速度場時，由于拾取工作量較大，因此僅人工拾取了20道，CMP間隔為25，而智能方法可以逐點進行拾取，考慮到速度場光滑性問題，最終選擇以5個CMP點為間隔進行拾取.圖7展示了人工拾取速度場(圖8a)、常規(guī)聚類方法拾取的速度場(圖8b)和本文方法拾取的速度場(圖8c)的比較.從圖可看出，受噪聲等因素的影響，常規(guī)聚類方法拾取速度場(圖8b)橫向上連續(xù)性差，且在0.5～1 s處部分道拾取到了高速異常值，準確度較低.圖8c為本文方法拾取的速度場，速度變化趨勢與人工拾取的速度場基本一致，與人工拾取結(jié)果(圖8a)基本吻合，但總用時是人工拾取的約0.1倍，平均單道拾取用時是人工拾取的0.02倍，拾取效率高.

圖8 實際數(shù)據(jù)速度場對比(a) 人工拾取速度場； (b) 常規(guī)聚類方法拾取速度場； (c) 自適應(yīng)閾值方法拾取速度場.Fig.8 Comparisons of picked stacking velocity structures from field data(a) Velocity structure picked by the manual picking； (b) Velocity structure picked by the conventional method； (c) Velocity structure picked by the adaptive threshold method.

為了進一步檢驗拾取結(jié)果的正確性，應(yīng)用圖7中三種方法拾取的速度場進行動校正和疊加處理.從疊加剖面(圖9)可以看出，本文方法的疊加剖面(圖9c)與人工結(jié)果(圖9a)相近，在0.7 s及1 s處的同相軸明顯比常規(guī)聚類方法(圖9b)更連續(xù)，成像更清晰.

圖9 三種方法實際數(shù)據(jù)疊加剖面(a) 人工拾取速度的疊加剖面； (b) 常規(guī)聚類方法拾取速度的疊加剖面； (c) 自適應(yīng)閾值方法拾取速度的疊加剖面.Fig.9 Stacked sections corresponding to three methods(a) Stacked section obtained by the manual picking； (b) Stacked section obtained by the conventional method； (c) Stacked section obtained by the adaptive threshold method.

3 結(jié)論

本文通過類比人工速度譜拾取過程，提出了一種基于自適應(yīng)閾值約束的無監(jiān)督聚類智能速度拾取方法.該方法利用時窗方法在速度譜中計算自適應(yīng)閾值，從而識別出一次反射速度能量團作為速度拾取的候選區(qū)域.隨后K均值聚類方法將不同的速度能量團聚類，并拾取疊加速度.最后，依據(jù)人工拾取速度的經(jīng)驗，對離群速度點進行后處理，為解決地震資料速度譜自動拾取問題提供了一種新的方案.

例子表明本文方法整體與人工拾取的精度相當，對速度變化的跟蹤較為準確，能夠快速拾取疊加速度曲線，比人工拾取快約50倍數(shù)量級，提高了拾取效率，降低了人工成本，為多維數(shù)據(jù)體逐點拾取速度場提供了可能.同比常規(guī)聚類智能拾取方法，自適應(yīng)閾值的引入提高了識別速度能量團精度，主要體現(xiàn)在識別弱反射，適用于拾取同一CMP道集內(nèi)能量強弱相間的反射波速度.離群速度處理這一步驟能規(guī)避部分速度異常值，以獲得更光滑的速度場，更符合人工解釋常識，避免加劇動校正子波拉伸現(xiàn)象.但本文方法在一定程度上仍受到速度譜精度以及逐個CMP道集處理方式的限制，下一步工作是從CMP道集角度以及多個CMP同時處理角度出發(fā)建立速度場，從而進一步提高實際速度拾取精度.