地震倒譜特征參數(shù)譜聚類地震相分析方法

桑凱恒 張繁昌* 李傳輝

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

0 引言

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

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

地震數(shù)據(jù)中包含豐富的地質(zhì)信息，但直接用地震數(shù)據(jù)劃分地震相的計算量巨大，還容易造成維數(shù)災(zāi)難。常規(guī)地震屬性可以提取地震數(shù)據(jù)中的特征信息，但種類繁多的地震屬性往往對應(yīng)不同的地質(zhì)特征，人為選擇地震屬性劃分地震相往往導(dǎo)致預(yù)測結(jié)果的多解性[19-20]。反射系數(shù)是地質(zhì)特征的直觀反映，但很難直接準(zhǔn)確求取地層反射系數(shù)。

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

1 方法原理

1.1 倒譜特征參數(shù)

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

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

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

(1)

式中n表示時間域采樣點(diǎ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域采樣點(diǎn)；Z-1{·}表示逆Z變換。

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

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

(3)

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

(4)

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

(5)

由自相關(guān)方法和Levinson-Durbin遞推算法[28]計算地震數(shù)據(jù)得到ai。

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

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

圖3 圖2b倒譜特征參數(shù)紅、綠、藍(lán)色曲線對應(yīng)圖2a中不同速度的地質(zhì)體

1.2 譜聚類

聚類作為一種不需要特定標(biāo)簽數(shù)據(jù)進(jìn)行訓(xùn)練的無監(jiān)督學(xué)習(xí)方法，已成為研究最多、應(yīng)用最廣的機(jī)器學(xué)習(xí)方法之一[15]。譜聚類方法能夠?qū)?shù)據(jù)集的聚類轉(zhuǎn)化為圖的劃分問題。圖中的一個結(jié)點(diǎn)對應(yīng)一個數(shù)據(jù)點(diǎn)，結(jié)點(diǎn)間的邊權(quán)值對應(yīng)數(shù)據(jù)間的相似度。圖的切割就是找到一條或幾條邊，以這些邊為切割點(diǎn)，將圖劃分為兩個或多個子圖，使同一個子圖數(shù)據(jù)點(diǎn)間相似度高，不同子圖數(shù)據(jù)點(diǎn)間相似度低[29-30]。

首先定義最小割集準(zhǔn)則，將全部樣本劃分為集合A和B，目標(biāo)函數(shù)Cut(A,B)為

(6)

根據(jù)拉普拉斯矩陣的性質(zhì)，目標(biāo)函數(shù)轉(zhuǎn)化為

(7)

因此，可以將圖的劃分問題轉(zhuǎn)化為尋找“全局最優(yōu)”問題，即

(8)

式中m為聚類樣本數(shù)。

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

(9)

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

本文利用Three-circle 數(shù)據(jù)集分析、對比K均值與譜聚類的效果(圖4)?？梢姡寒?dāng)相同標(biāo)記的樣本為同一種顏色，不同標(biāo)記的樣本為不同顏色時，分類精度更高；K均值的分類精度較低，這是由于其基于原型數(shù)據(jù)歐式距離進(jìn)行分類，只考慮數(shù)據(jù)之間的歐式距離，因此對復(fù)雜結(jié)構(gòu)數(shù)據(jù)集無法得到滿意的分類結(jié)果(圖4左)；譜聚類將原始數(shù)據(jù)映射到特征向量的數(shù)據(jù)空間中進(jìn)行劃分，因此對復(fù)雜結(jié)構(gòu)數(shù)據(jù)集也能以極高的精度實(shí)現(xiàn)樣本聚類(圖4右)。

圖4 Three-circle 數(shù)據(jù)集K均值聚類(左)、譜聚類(右)效果不同的標(biāo)記點(diǎn)代表真實(shí)不同分類，顏色代表分類結(jié)果

1.3 方法技術(shù)路線

根據(jù)前文的方法原理，建立方法的技術(shù)路線：

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

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

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

(6)選擇特征向量構(gòu)成特征向量空間，在特征向量空間采用常規(guī)聚類算法進(jìn)行聚類，得到地震相自動分析結(jié)果。

2 模型與實(shí)際數(shù)據(jù)測試

2.1 模型試驗(yàn)

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

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

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

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

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

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

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

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

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

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

(10)

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

(11)

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

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

2.2 實(shí)際數(shù)據(jù)測試

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

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

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

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

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

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

圖12 目標(biāo)層段鉆井資料

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

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

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

3 討論

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

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

4 結(jié)論與認(rèn)識

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

(1)以圖論為基礎(chǔ)的譜聚類方法將數(shù)據(jù)的聚類轉(zhuǎn)化為圖的分割問題，通過圖的最優(yōu)分割實(shí)現(xiàn)數(shù)據(jù)的精確聚類。相較于傳統(tǒng)以原型數(shù)據(jù)聚類的方法，譜聚類對復(fù)雜高維數(shù)據(jù)的劃分能力更強(qiáng)，劃分精度更高。通過優(yōu)化相似度矩陣計算方法，構(gòu)建稀疏相似度矩陣，可以降低矩陣維度過大引起的存儲和計算量大的問題，使譜聚類更適用于劃分三維空間地震相。

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

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