基于超虛干涉法約束的模糊C均值聚類地震初至自動提取

譚家煒 李 靜* 李飛達 曾昭發

（①吉林大學地球探測科學與技術學院，吉林長春130021；②吉林省勘查地球物理研究院，吉林長春130062）

0 引言

在地震數據處理中，初至拾取的精度和效率是一個基本且重要的問題。初至拾取尤其在靜校正、層析成像等地震數據處理過程中具有重要作用，進而直接影響后期綜合地質解釋。

早期人工地震初至拾取方法雖然精度高，但效率低下，限制了地震數據處理的速度；另外，因拾取精度受人為因素操控，穩定性難以保障［1］。針對現今海量的地震數據，選用自動初至拾取方法，則可大幅度減少人為因素干擾并提高工作效率。

早在20世紀70年代初，Peraldi等［2］就率先提出基于相鄰地震道互相關獲得初至到達與信號峰值之間的時差，但當相鄰道波形相差很大或存在缺失道時，該方法難以奏效。

Coppens［3］于80 年 代 創 立 的 能 量 比 值 法 已 成為迄今應用最普遍的自動初至拾取方法。即通過計算地震信號的長短時窗特征函數的比值（STA／LTA），設置觸發閾值，比值大于閾值的第一個時刻即為地震初至波到達時刻［4］。但該方法易受初至波到達前的隨機干擾影響［5-7］。

Boschetti等［8］曾提出基于分形維方法的初至拾取方法，但因需做插值處理，導致計算速度慢，且其結果取決于插值的準確性［9］。

李輝峰等［10］提出基于邊緣檢測的初至拾取方法，通過把待處理的地震數據轉變成灰二值化灰度圖，再進行邊緣檢測。該方法雖能有效地拾取地震初至，但需大量訓練數據以提高拾取的準確性，耗時長，效率低。在低信噪比情況下難以達到預期效果，需人工干預。

譚玉陽等［11］設計了一種SLPEA 算法用于微地震自動初至拾取，利用地震信號與環境噪聲之間的差異，根據信號信噪比構造不同檢測函數拾取初至。該方法克服了傳統方法抗噪性弱的缺點，但計算量大、耗時長，且時窗的選擇依據不明確，依賴經驗或反復調試［12］。

隨著地震勘探的不斷深入，上述方法已無法滿足日益增長的龐大數據量要求，所以需要一種人工干預小的準確地震初至拾取方法，尤其是針對低信噪比的、海量的地震數據情形。

基于上述成果調研，本文提出一種基于超虛干涉（Super-virtual interferometry，SVI）約束的模糊C均值（Fuzzy C-means，FCM）聚類地震自動初至拾取方法。FCM 聚類分析是一種非監督的機器學習方法，只使用數據本身對數據進行分類［13］。與嚴重依賴于輸入訓練數據的監督機器學習方法不同，聚類分析僅依賴于數據本身，因此更加靈活，更容易應用于實際的地震初至拾取。除此以外，由于球面擴散、衰減和環境噪聲的影響，遠炮檢距初至波的信號弱，噪聲強，信噪比低，特征值非常混亂，常常無法直接拾取。本文利用SVI法加強初至波的能量，提高地震數據的信噪比［14-16］，在此基礎上開展FCM 聚類分析，并自動拾取低信噪比地震資料的初至。

1 方法原理

1.1 FCM 聚類算法

地震記錄由有效信號和噪聲兩部分構成。有效信號部分的第一個時刻被視為地震的發生，即初至波的到達時刻。因此，在給定一組數據的情況下，初至拾取的本質就變成了分類問題。當一組訓練數據與預定義的數據特征一起給出時，這種分類問題可被視為監督分類問題。若只使用數據本身的屬性特征對地震記錄進行分類，該問題就變成了典型的聚類分析問題，即非監督機器學習［17］。

模糊聚類分析是一種利用模糊數學語言對事物按照一定標準進行分類的數學方法［18］，通過利用聚類對象的屬性特征建立模糊矩陣，并在其之上按適當隸屬度進行聚類分析。因此每個點可屬于具有不同隸屬度的兩個或更多聚類，即不再像傳統分類一樣只屬于某一類，而是可通過多個隸屬度矩陣表示數據點屬于不同分類的概率［19］。FCM 是一種目標函數迭代優化基礎上的非監督聚類分析方法。相比于監督機器學習方法，非監督學習方法只對數據本身進行模糊聚類分析且在實現過程中不需人為干預；另外，無需大量訓練樣本，相比于監督機器學習方法，其計算速度更快。

FCM 方法的實現過程：首先是構建目標函數；然后不斷迭代更新隸屬度矩陣和聚類中心，讓目標函數取到極小值；最終根據最小誤差原則獲取較好聚類分析結果。該聚類結果并非每個數據點屬于哪個分類，而是每個數據點對聚類中心的隸屬概率，它允許一個數據屬于兩個或多個聚類［20］。

FCM 是基于以下目標函數的最小化［21］

式中：m 為模糊指標；xi表示nx維數據中的第i個數據點（nx是xi的維數）；ui，j表示分類j 中xi的隸屬度；cj表示nc維的聚類中心；‖·‖表示任意數據點與聚類中心之間的相似性的某種范數（如L2）；C 表示聚類數目；N 表示輸入數據的個數。

xi和cj是相同維度向量，其維數由特征函數的數量決定，本文用均值M、能量E 和長短時窗均值（STA、LTA）之比R 作為初至拾取的三個特征函數，每個xi由對應采樣點i處所有特征函數值組成。

三個特征函數定義如下

式中：Tw為時窗長度；d（i）為地震記錄；NLTA、NSTA表示長、短時窗長度。

通過式（1）的迭代優化進行模糊分類，其中隸屬度ui，j和聚類中心cj的更新表達式分別為

當maxi，j｛|u（k＋1）i，j－u（k）i，j|｝＜ε 時，停止 迭代，其中ε為0～1之間的終止閾值，k 為迭代次數。此過程最后收斂于目標函數J 的局部極小值或鞍點。

為了驗證FCM 方法的有效性，合成了簡單的無噪聲（圖1a）和含噪聲（圖1b，信噪比為1d B）的單道數據，得到對應數據的三個特征函數（圖2）。

圖1 無噪（a）和含噪（b，1dB）簡單合成數據

圖2 無噪（a）和含噪（b，1d B）數據聚類算法特征函數

圖3為簡單合成數據計算得到的拾取結果和隸屬度，可見利用三個特征函數（圖3a、圖3b為無噪聲，圖3c、圖3d為含噪聲）同時進行聚類分析能得到較好結果，可拾取有效信號，拾取位置大致相同，即證明該方法能有效地拾取初至。利用三個特征函數對計算互相約束，可提高聚類分析的穩定性，增強方法的抗干擾能力。

1.2 SVI方法

SVI是一種完全數據驅動的方法，不依賴地震數據的頻率或視速度信息［22］。傳統的SVI方法使用互相關和卷積疊加的方法，通過大幅提高信噪比增強初至信號［23］。但經過處理后，地震子波旁瓣數量增加，旁瓣能量也同時提升。Nakata等［24］提出了基于互相干代替互相關的地震數據干涉重構方法，很好地從本質上避免了SVI處理后初至波產生的虛假旁瓣。Place等［25］將該方法應用到淺層主動源地震數據中，并對比與互相關應用效果的差別。

首先，對檢波點R1與R2處記錄的信號做相互干（圖4a），得到虛擬記錄，它表示由虛擬源R′1產生的信號（包括負旅行時間，圖4a中的虛線）在檢波點R2處被接收。對位于檢波點R1臨界炮檢距（初至波由直達波變為折射波所對應的炮檢距）后的任何震源S 重復該過程，可產生相同效果的虛擬記錄。對這些記錄疊加求和，增強由R′1產生的虛擬信號，并衰減非相干噪聲，以提高數據信噪比。當源位于位置S 時，將疊加的虛擬記錄與檢波點R1的原始記錄進行卷積，產生SVI記錄，就像由震源S產生在檢波點R2接收的記錄一樣。對于S 與R2之間超出臨界距離的所有檢波點R1，可獲得相同效果的SVI記錄，疊加求和進一步提高數據信噪比［26］。

圖3 聚類分析的計算結果

圖4 SVI方法原理示意圖

將兩個不同點R1和R2處記錄的兩條地震數據的互相關在傅里葉域中定義為

式中：u（R1，S）表示震源S 在檢波點R1記錄的地震波場；*表示復共軛。

互相干定義為頻率歸一化的互相關，表示為

式中|·|表示振幅譜。

互相干表達式（式（8））的分子與互相關表達式（式（7））相同，分母為波的振幅譜的乘積。與互相關相比，互相干忽略振幅信息的影響，只使用相位信息，有效地從提取的地震信號中去除震源的特征，即使在強噪聲情況下也能得到穩定的結果，因此互相干比互相關更穩健。

忽略噪聲的存在，可將從震源S 處到R 處記錄的地震波場寫為傅里葉域中的乘積

式中：W（S）是信號源在S 處產生的信號；G（R，S）是R 與S 之間的格林函數。

恢復震源特征是SVI方法的重要步驟［25］，將式（9）分別代入式（7）和式（8），可得

虛擬記錄與實際原始記錄卷積形成SVI地震記錄，重新恢復震源的特征（圖4b），表示為

式中uSVI表示SVI處理后的波場，第二個總和項表示疊加的互相干函數（圖4a）。

2 初至拾取流程

基于SVI方法約束的FCM 地震初至拾取流程如圖5所示，具體實現過程如下。

（1）加窗處理獲取包含初至信號的地震資料——加窗可避免對全波場進行運算，降低計算量并減少虛假同相軸的產生；

（2）選取某炮，將檢波點R1與R2的波場做互相干處理，得到虛擬記錄；

（3）將每炮互相干處理結果做疊加求和，增強虛擬記錄的能量，提高信噪比，得到檢波點R1與R2之間的虛擬記錄；

（4）選取震源S，將檢波點R1波場與檢波點R1和R2之間虛擬記錄做卷積處理，得到SVI記錄；

（5）對每道卷積結果疊加求和，增強SVI記錄，再次提高信噪比，得到最終經SVI處理的信號；

（6）將處理后的記錄作為FCM 算法的輸入數據進行計算，給定算法參數，如聚類數目C、模糊指標m、迭代停止條件ε；

（7）初始化隸屬度矩陣U，U（0）＝［ui，j］；

（8）利用U（k）根據式（6）計算聚類中心c（k）＝［cj］；

（9）根據式（5）更新U（k），得到U（k＋1）；

（10）如果‖U（k＋1）－U（k）‖＜ε，則停止迭代，輸出拾取結果；否則k＝k＋1，返回第3步繼續計算。

圖5 基于SVI方法約束的FCM 地震初至拾取流程

3 理論模型測試

3.1 水平層狀模型

為了驗證本文方法的有效性和穩定性，建立了一個三層水平層狀縱波速度模型（圖6a）并進行聲波有限差分正演模擬（圖6b）。第一層（0～20m）、第二層（20～40m）、第三層（40～60m）的速度依次為800、1600、2000m／s，采用主頻為30Hz的雷克子波，時間采樣間隔為0.5ms，炮間距為5m，道間距為5m。在正演數據中加入不同程度的高斯隨機噪聲，模擬不同信噪比時的初至拾取。

圖7為當原始記錄SNR＝25dB 時的FCM 方法和STA／LTA 方法的單炮初至拾取結果（藍色圓點）。圖8為使用基于射線追蹤方法的層析成像方法對兩種方法拾取結果進行縱波速度層析成像，圖中黑色虛線為真實模型的速度分界面。可見當噪聲不太強時，FCM 方法能較準確地拾取初至，反演結果也與實際模型較相符；而STA／LTA 方法的拾取雖大致準確，但由于參數選擇固定，存在一定誤差；縱波速度成像結果與真實模型不能較好匹配。

當SNR＝15d B 時（圖9），FCM 和STA／LTA方法初至拾取都不準確，FCM 方法只能拾取信噪比較高的近炮檢距信號，而STA／LTA 方法受噪聲影響更嚴重，錯誤拾取非常多。

使用基于SVI的FCM 方法進行拾取（圖10a），經過SVI處理后，折射波能量明顯加強，噪聲受到壓制，遠炮檢距信噪比明顯提高，可有效判定并識別真實的初至信號。從其反演結果（圖10b）看出，與真實速度模型較為匹配，進一步說明拾取的準確性。圖11a和圖11b分別為對原始數據使用FCM 拾取和使用基于SVI的FCM 拾取的初至時間二維圖。在原始數據中只能拾取近炮檢距信號，遠炮檢距無法拾取，但經SVI處理后，遠炮檢距拾取效果明顯改善。

圖6 水平層狀模型（a）及其正演模擬數據（b）

圖7 SNR＝25dB時FCM（a）和STA／LTA（b）方法地震記錄初至拾取結果

圖8 FCM（a）和STA／LTA（b）方法縱波速度層析成像結果

圖9 SNR＝15dB時FCM（a）和STA／LTA（b）方法地震記錄初至拾取結果

圖10 SNR＝15dB時基于SVI的FCM 初至拾取結果（a）和縱波速度層析成像結果（b）

圖11 針對原始數據單獨FCM（a）和基于SVI的FCM（b）方法初至時間二維圖

3.2 臺階速度模型

為進一步驗證本文方法的有效性，建立如圖12a所示的臺階速度模型。其上層速度為800m／s，左側深度為15m，右側深度為30m，臺階在150m 處，下層速度為2000m／s，時間采樣間隔為0.5ms，采用主頻為30Hz的雷克子波，炮間距為4m，道間距為4m。正演得到的地震記錄如圖12b所示。

在該數據中加入不同程度的高斯隨機噪聲。圖13a和圖13b分別為當原始記錄SNR＝25dB時的FCM 方法和STA／LTA 方法的單炮初至拾取結果（藍色圓點）。圖14a和圖14b為使用基于射線追蹤方法的初至層析對兩種方法的拾取結果進行層析成像，黑色虛線為真實模型的速度分界面。可知當噪聲強度不太強時，FCM方法能較準確拾取初至，反演結果也與實際模型較為匹配，只在右側分界面稍上移，原因是200m 炮檢距后的拾取不夠準確。而STA／LTA 方法的拾取在近炮檢距也大致準確，但遠炮檢距的拾取結果不準確，基本變成直線，無細節變化，從反演結果（圖14b）可見，導致臺階右側的反演結果與實際模型不能匹配。

圖12 臺階理論模型（a）及其正演模擬數據（b）

圖13 SNR＝25dB時地震記錄初至拾取結果

圖14 縱波速度層析成像結果

使用基于SVI的FCM 方法進行初至拾取，結果見圖13c。與只使用FCM 方法拾取相比，遠炮檢距拾取更準確，從反演結果（圖14c）可知，右側高速層反演效果明顯改善，與實際模型較吻合。

當SNR＝15dB時（圖15），遠炮檢距的初至信號已經無法分辨。FCM 方法和STA／LTA 方法初至拾取都不準確，遠炮檢距都只能錯誤地拾取到直達波上。使用基于SVI的FCM 方法進行拾取，經過SVI處理后（圖15c），初至信號能量明顯加強，噪聲受到壓制，信噪比明顯提高，可準確地拾取初至信號。由反演結果圖（圖15d）看出，反演模型較為匹配，進一步說明拾取結果的準確性。

圖15 SNR＝15dB時地震記錄初至拾取結果及基于SVI的FCM 拾取后層析成像結果

4 實際數據應用

選取沙特吉達地區實際地震資料進行應用測試。數據剖面長度為1200m，采樣間隔為1ms，共包含212炮記錄，每炮240道接收，道間距和炮間距均為5m。圖16為該資料的單炮記錄，該資料遠炮檢距初至信號極弱，信噪比極低，噪聲干擾非常嚴重，初至難以準確拾取。

原始地震數據經過SVI處理后（圖17），遠炮檢距的噪聲得到有效壓制，初至同相軸清晰可辨。由于信噪比的提高，能夠拾取遠炮檢距的初至信號。圖18a和圖18b分別為使用基于SVI的FCM 拾取的和只使用FCM 拾取的初至二維圖，對比看出SVI處理后遠炮檢距信息拾取效果明顯提高。

圖19a為利用SVI處理后的初至拾取結果進行層析成像的速度成像結果。根據前期在該區域開展的地質調查和相關地球物理研究成果［27］，測線橫跨一條南北向的Qademah斷層，該斷層緣于紅海海域隆起形成的不連續速度結構。根據前人研究成果，在測線300m 附近是斷層位置。圖19a所示的速度成像結果中的低速下凹與地質情況吻合。為了驗證反演結果的可靠性，與共炮檢距剖面（圖19b）進行對比。在缺少鉆井資料情況下，可認為共炮檢距剖面能夠反映重要的地下結構分布特征［28］。圖19b所示的共炮檢距剖面水平方向不連續分布與速度成像高低速變化帶有很好的一致性，特別是在橫向300m 和700m 附近的低速不連續帶吻合很好。

圖16 實際地震數據單炮記錄

圖17 基于SVI的FCM 初至拾取結果

圖18 針對實際資料單獨FCM（a）和基于SVI的FCM（b）方法初至時間二維圖

圖19 實際地震縱波速度層析成像（a）及共炮檢距剖面與速度成像剖面擬合圖（b）

5 結論

本文針對低信噪比地震資料中，特別是遠炮檢距初至信號弱、噪聲強，無法準確地自動拾取初至信號的問題，提出了利用基于SVI約束的FCM 聚類分析初至自動拾取方法。其基本思路是首先通過SVI方法提高地震記錄遠炮檢距信號的信噪比，在此基礎上，利用FCM 算法的初至拾取方法對加強后的信號進行自動拾取。理論模型正演數據和實際地震資料對方法的有效性和適用性進行了驗證，即實現了在低信噪比數據條件下，自動、快速拾取地震記錄初至，達到了滿意的計算效率和計算精度。通過層析成像速度反演結果也驗證了拾取結果的可靠性。本文研究方法為地震數據處理提供了一種穩健、可靠的技術手段。