地震倒譜特征參數譜聚類地震相分析方法

桑凱恒 張繁昌* 李傳輝

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580；②中國地質大學(北京)地球物理與信息技術學院，北京 100083)

0 引言

“相”是一定巖層生成時的古地理環境及其物質表現的總和[1]，而地震相是沉積相的地震響應。由于沉積環境對油氣形成具有重要影響，所以準確分析地震相對油氣勘探和油藏評價具有指導意義[2-4]。傳統的地震相劃分主要由人工直接在地震剖面上解釋，俗稱“相面法”，人為因素較多，具有極大的主觀性和不確定性，同時受限于地震資料分辨率，無法滿足精細勘探的精度需求[5]。

隨著地震數字處理技術和計算機技術的發展，地震相分析進入數據分析時代，目前常用的地震相分析方法有隨機模擬、神經網絡[6-7]、聚類算法[8-10]和深度學習[11-12]等。隨機模擬以地質統計學為基礎，以已知信息為約束條件，通過隨機函數產生等概率地震相模型[13-14]。但隨機模擬結果易受隨機模型影響，而且在地質結構復雜地區難以準確劃分地震相。神經網絡和深度學習具有較強的容錯性和泛化能力，但需要海量訓練樣本數據，同時訓練網絡的計算開銷巨大。聚類算法作為最常用的無監督學習方法之一，通過對無標記樣本的學習揭示數據的內在性質和規律[15-16]。K均值聚類、C模糊聚類等經典聚類算法在簡單數據集上均獲得了理想的聚類結果，但對于非凸數據集并不能實現全局最優。以圖論為基礎的譜聚類方法能夠在任意形狀的數據空間實現全局最優劃分[17-18]。因此本文選擇譜聚類方法劃分地震相，通過與井資料的標定，建立地震相與地質體間的對應關系。

地震數據中包含豐富的地質信息，但直接用地震數據劃分地震相的計算量巨大，還容易造成維數災難。常規地震屬性可以提取地震數據中的特征信息，但種類繁多的地震屬性往往對應不同的地質特征，人為選擇地震屬性劃分地震相往往導致預測結果的多解性[19-20]。反射系數是地質特征的直觀反映，但很難直接準確求取地層反射系數。

本文提出一種基于地震數據倒譜特征參數的譜聚類地震相分析方法。該方法以地震倒譜特征參數為譜聚類的輸入變量劃分地震相，然后通過井標定，建立地震相與地質體間的對應關系。其優點是：一方面能夠提取反射系數信息作為地震相劃分依據；另一方面可以將數據集的聚類轉化為圖的最優分割問題，能夠適應復雜數據空間的分類。通過模型和實際數據驗證了方法的可行性。

1 方法原理

1.1 倒譜特征參數

地震勘探中，地層反射系數直觀反映了地質特征，因此可以利用反射系數作為輸入信息分析地震相。目前從地震數據中提取反射系數信息的常用方法有反褶積[21-23]、譜反演[24-25]等，但仍很難準確求取地層反射系數。Xie等[26]認為，語音信號處理中的聲學模型與地震勘探中的褶積模型相近，可以借鑒語音信號的處理方法提取地震反射系數信息。因此本文從地震數據中提取反射系數的表征參數——地震倒譜特征參數，代替反射系數作為譜聚類的輸入特征向量分析地震相。

假設s(n)為長度為N的一段地震信號，根據褶積理論，在時間域s(n)為子波x(n)與反射系數h(n)的褶積

s(n)=x(n)*h(n)

(1)

式中n表示時間域采樣點。將式(1)變換到倒譜域

Z-1[lnS(z)]=Z-1[lnX(z)]+Z-1[lnH(z)]

(2)

式中：S(z)、X(z)、H(z)分別為s(n)、x(n)、h(n)的Z域表示，z表示Z域采樣點；Z-1{·}表示逆Z變換。

式(2)表明，對信號進行Z變換后取對數，再進行逆Z變換就得到信號的倒譜域表示。在倒譜域褶積地震信號轉化為加性信號，容易分離地震子波與反射系數。利用模型數據驗證式(2)的正確性(圖1)。可見，倒譜域地震記錄的振幅(圖1c右)為反射系數振幅(圖1a右)與子波振幅(圖1b右)之和，實驗結果驗證了式(2)。根據式(2)，有

圖1 反射系數(a)、Ricker子波(b)、地震數據(c)的時間域(左)和倒譜域(右)表示

(3)

對于倒譜特征參數的計算，參照聲紋識別技術[27]，建立H(Z)的全極點模型

(4)

式中：ai為線性預測系數(Linear Predictive Coefficients，LPC)；p為線性預測階數；G為預測模型增益。

(5)

由自相關方法和Levinson-Durbin遞推算法[28]計算地震數據得到ai。

為了分析提取倒譜特征參數的有效性，建立楔形體模型(圖2a)，用30Hz雷克子波正演得到合成地震記錄，抽取其中三分之一(圖2b)進行分析；然后以單道地震記錄為輸入，提取地震數據的1～15階倒譜特征參數(圖3)。可見，與正演地震記錄30、60道位置反映的速度體邊界(圖2b)相比，倒譜特征參數反映的邊界(圖3)更清楚。說明地震倒譜特征參數可以充分提取反射系數信息，可作為譜聚類地震相分析的輸入變量。

圖2 楔形體模型(a)及其正演地震記錄(b)30 、60道位置對應楔形體中不同速度體邊界

圖3 圖2b倒譜特征參數紅、綠、藍色曲線對應圖2a中不同速度的地質體

1.2 譜聚類

聚類作為一種不需要特定標簽數據進行訓練的無監督學習方法，已成為研究最多、應用最廣的機器學習方法之一[15]。譜聚類方法能夠將數據集的聚類轉化為圖的劃分問題。圖中的一個結點對應一個數據點，結點間的邊權值對應數據間的相似度。圖的切割就是找到一條或幾條邊，以這些邊為切割點，將圖劃分為兩個或多個子圖，使同一個子圖數據點間相似度高，不同子圖數據點間相似度低[29-30]。

首先定義最小割集準則，將全部樣本劃分為集合A和B，目標函數Cut(A,B)為

(6)

根據拉普拉斯矩陣的性質，目標函數轉化為

(7)

因此，可以將圖的劃分問題轉化為尋找“全局最優”問題，即

(8)

式中m為聚類樣本數。

根據瑞利商性質，式(8)可以轉換為求解L的特征值問題，即使fT(D-W)f最小的解為L最小特征值對應的特征向量。將二分類問題推廣到k類分類，譜聚類可表示為以下優化問題

(9)

在實際應用中發現，上述經典譜聚類方法存在相似度矩陣維數過大的問題。如對于10000個樣本點的數據集，相似度矩陣大小為10000×10000，難以實現矩陣的存儲和計算。因此本文優化了相似度矩陣計算方法，構建了稀疏相似度矩陣。首先設定稀疏系數ξ，對于某一樣本點的相似度向量Wi=[wi1,wi2,…，wim]，只保留相似度最大的前ξ個相似系數，其余以0值存儲。通過對所有樣本相似度向量的稀疏化，可以得到稀疏相似度矩陣。然后利用稀疏矩陣特征值分解的方法求取相似度矩陣的特征值和特征向量。

本文利用Three-circle 數據集分析、對比K均值與譜聚類的效果(圖4)。可見：當相同標記的樣本為同一種顏色，不同標記的樣本為不同顏色時，分類精度更高；K均值的分類精度較低，這是由于其基于原型數據歐式距離進行分類，只考慮數據之間的歐式距離，因此對復雜結構數據集無法得到滿意的分類結果(圖4左)；譜聚類將原始數據映射到特征向量的數據空間中進行劃分，因此對復雜結構數據集也能以極高的精度實現樣本聚類(圖4右)。

圖4 Three-circle 數據集K均值聚類(左)、譜聚類(右)效果不同的標記點代表真實不同分類，顏色代表分類結果

1.3 方法技術路線

根據前文的方法原理，建立方法的技術路線：

(1)根據目標要求，按照一定大小的時窗提取每道地震記錄；

(2)以地震記錄為輸入，利用自相關方法和Levinson-Durbin遞推算法計算線性預測系數ai；

(5)根據W求出相應正則化或非正則化L，并通過對L進行特征分解求出特征值及其對應特征向量；

(6)選擇特征向量構成特征向量空間，在特征向量空間采用常規聚類算法進行聚類，得到地震相自動分析結果。

2 模型與實際數據測試

2.1 模型試驗

為了驗證上述地震相分析方法的有效性，建立三維地質模型(圖5a)，假設同一地震相具有相同或相近波阻抗數值。用30Hz雷克子波正演得到合成地震記錄(圖5b)。

圖5 三維地質模型 (a)及其正演地震記錄(b)

按照本文方法利用正演地震記錄得到聚類簇數k分別為3、6和10的地震相劃分結果(圖7、圖8和圖9)，倒譜特征參數階數為1～15。

通過直觀對比、分析發現：①K均值聚類方法在原始數據空間計算距離，對于復雜數據劃分能力有限，地震相劃分結果(圖7c、圖7d、圖8c、圖8d、圖9c、圖9d)與真實地震相(圖7a、圖8a、圖9a)差距較大，而且隨著聚類簇數增加，劃分結果的隨機性也逐漸增大，表現為：以地震記錄為輸入的K均值聚類結果(圖7d、圖8d、圖9d)存在很大的隨機性，可解釋性差，即使聚類簇數較低仍存在很大誤差；以地震倒譜特征參數為輸入的K均值聚類結果(圖7 c、圖8 c、圖9 c)可以消除地震波形的影響，充分利用反射系數劃分地震相，可提高劃分精度。②譜聚類方法將數據從原始數據空間映射到特征向量空間，避免直接計算數據距離，在復雜數據空間中仍有較高的劃分精度，地震相劃分結果(圖7b、圖8b、圖9b)與真實地震相(圖7a、圖8a、圖9a)吻合良好。為了說明聚類簇數為3、6和10的波阻抗的數值范圍，統計了每一分類結果的類內波阻抗數值范圍(表1)。可見，與K均值聚類地震相劃分結果(c、d)相比，譜聚類地震相劃分結果(b)不同類間波阻抗數據重疊最小，與真實地震相(a)的波阻抗范圍基本一致，吻合度更高。

表1 k=3、6、10的類內波阻抗數值范圍 單位：×106 kg·m-2·s-1

圖6 圖5a中間層阻抗平面圖

圖7 地震相劃分結果(k=3)(a)真實地震相(根據波阻抗數據進行均分)；(b)譜聚類(以地震倒譜特征參數為輸入)；(c)K均值聚類(以地震倒譜特征參數為輸入)；(d)K均值聚類(以地震記錄為輸入)

圖8 地震相劃分結果(k=6)(a)真實地震相(根據波阻抗數據進行均分)；(b)譜聚類(以地震倒譜特征參數為輸入)；(c)K均值聚類(以地震倒譜特征參數為輸入)；(d)K均值聚類(以地震記錄為輸入)

圖9 地震相劃分結果(k=10)(a)真實地震相(根據波阻抗數據進行均分)；(b)譜聚類(以地震倒譜特征參數為輸入)；(c)K均值聚類(以地震倒譜特征參數為輸入)；(d)K均值聚類(以地震記錄為輸入)

模型尺寸為750×737×300，縱向上分為三層，每層時間厚度均為100ms，其中上、下兩層為均勻層，速度為3000m/s，密度為2.2g/cm3，中間層由Marmousi模型(圖6)縱向疊置組成為了更直觀地定量評價模型劃分精度，定義Rand指數

(10)

式中S、N為樣本集合，即

(11)

圖10為不同方法模型劃分的RI。由圖可見：以地震記錄為輸入的K均值聚類地震相劃分精度較低(RI約為0.5)；以地震倒譜特征參數為輸入的K均值聚類地震相劃分精度明顯提高(RI約為0.7)；以地震倒譜特征參數為輸入的譜聚類地震相劃分精度最高(RI均大于0.9)，再次驗證了本文方法的有效性。

圖10 不同方法模型劃分的RIb、c、d分別代表譜聚類(以地震倒譜特征參數為輸入)、K均值聚類(以地震倒譜特征參數為輸入)、K均值聚類(以地震記錄為輸入)

2.2 實際數據測試

利用實際地震資料測試本文方法的技術路線。圖11為P區Inline和Crossline地震剖面，目的層段為Es3下砂礫巖儲層，層位S3為Es3下的頂。圖12為目標層段鉆井資料。利用疊后三維地震數據分析地震相宏觀展布，詳細流程如下。

圖11 P區Inline(a)和Crossline(b)地震剖面

(1)首先以S3為頂、向下50ms為底，截取時窗內的地震數據；

(2)對截取的單道地震數據提取1～15階地震倒譜特征參數，并做歸一化處理；

(3)以提取的地震倒譜特征參數為輸入，采用譜聚類方法對全部地震道劃分地震相。

圖13為P區地震相劃分結果。由圖可見：①與地震瞬時振幅地震相劃分結果(圖13b)相比，譜聚類地震相劃分結果(圖13a)的地震相邊界更清晰，可解釋性更強。②與地震瞬時振幅地震相劃分結果(圖13b)相比，多地震屬性地震相劃分結果(圖13c)的湖底扇體邊界更清楚；但在湖底扇區域外，A5、A6井以南的區域仍顯示大片扇體發育區(圖13c)，與古地貌(圖14)不符，因此地震相劃分精度低于譜聚類。鉆井資料(圖12)揭示，A3、A4、A7、A8、A9、A10井位于扇體發育區，A1、A2、A5、A6、A11井位于扇體邊緣或扇體不發育區。

圖12 目標層段鉆井資料

圖14為P區Es3下古地貌圖。對比圖13a與圖14可見，前者的扇體發育區對應于后者的湖底區域，兩者劃分的邊界吻合良好，驗證了譜聚類地震相分析的準確性。

圖13 P區地震相劃分結果(a)譜聚類；(b)地震瞬時振幅；(c)多地震屬性聚類結果中的紅色和黃色區域標定為扇體發育區，綠色和藍色區域標定為扇體欠發育區

圖14 P區Es 3下古地貌圖綠色為古地貌中的湖底部分

3 討論

本文提出一種基于地震數據倒譜特征參數的譜聚類地震相分析方法，其中倒譜特征參數階數的選擇對地震相分析結果具有重要影響。李韻[27]指出，在一定范圍內，線性預測階數與倒譜特征參數包含的信息量呈成正比，但是超過了某個值，高階參數包含的信息量很少，而且會增加計算難度和量度，因此線性預測階數的最佳取值范圍為8～32。鐘林鵬[31]通過實驗論證了倒譜特征參數有效信息主要集中在前幾階。分析式(5)發現：當倒譜特征參數階數小于線性預測階數時，每一階倒譜特征參數都會有新的線性預測系數加入；當倒譜特征參數的階數大于線性預測階數時，每一階倒譜特征參數不再有新的線性預測系數加入，只是前幾階倒譜特征參數值的累加，因此包含的信息量有限，與文獻[27,31]結果一致。

為了在實際地震數據測試中選擇合適階數的倒譜特征參數，要在理論分析的基礎上根據實際情況決定。對于大尺度目標，如目標層段為數百毫秒的地震數據，可以取較大的線性預測階數，更利于挖掘宏觀信息，相應地可取較高的倒譜特征參數階數。對于小尺度目標，都應取較小的線性預測階數和倒譜特征參數階數。

4 結論與認識

本文以地震倒譜特征參數為輸入變量，研究了利用譜聚類分析地震相的方法。通過模型測試和實際數據應用，取得如下認識：

(1)以圖論為基礎的譜聚類方法將數據的聚類轉化為圖的分割問題，通過圖的最優分割實現數據的精確聚類。相較于傳統以原型數據聚類的方法，譜聚類對復雜高維數據的劃分能力更強，劃分精度更高。通過優化相似度矩陣計算方法，構建稀疏相似度矩陣，可以降低矩陣維度過大引起的存儲和計算量大的問題，使譜聚類更適用于劃分三維空間地震相。

(2)以地震數據的倒譜特征參數代替原始地震記錄作為地震相劃分的輸入變量，一方面能減少數據的維數，減少計算復雜度和數據量；另一方面能消除波形的影響，提取直接反映地質情況的反射系數信息，提高劃分精度。

(3)實際數據應用表明，與地震瞬時振幅、多地震屬性地震相劃分結果相比，文中方法劃分的地震相帶與古地貌吻合更好，邊界更清晰，可解釋性也更強。將地震相劃分結果與井資料進行標定可以劃分不同的沉積單元，為油氣勘探和油藏評價提供數據支撐。