桑凱恒 張繁昌* 李傳輝

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

0 引言

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

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

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

本文提出一種基于地震數(shù)據(jù)倒譜特征參數(shù)的譜聚類地震相分析方法。該方法以地震倒譜特征參數(shù)為譜聚類的輸入變量劃分地震相，然后通過井標(biāo)定，建立地震相與地質(zhì)體間的對(duì)應(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)為，語(yǔ)音信號(hào)處理中的聲學(xué)模型與地震勘探中的褶積模型相近，可以借鑒語(yǔ)音信號(hào)的處理方法提取地震反射系數(shù)信息。因此本文從地震數(shù)據(jù)中提取反射系數(shù)的表征參數(shù)——地震倒譜特征參數(shù)，代替反射系數(shù)作為譜聚類的輸入特征向量分析地震相。

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

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

(1)

式中n表示時(shí)間域采樣點(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)表明，對(duì)信號(hào)進(jìn)行Z變換后取對(duì)數(shù)，再進(jìn)行逆Z變換就得到信號(hào)的倒譜域表示。在倒譜域褶積地震信號(hào)轉(zhuǎn)化為加性信號(hào)，容易分離地震子波與反射系數(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)的時(shí)間域(左)和倒譜域(右)表示

(3)

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

(4)

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

(5)

由自相關(guān)方法和Levinson-Durbin遞推算法[28]計(jì)算地震數(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道位置對(duì)應(yīng)楔形體中不同速度體邊界

圖3 圖2b倒譜特征參數(shù)紅、綠、藍(lán)色曲線對(duì)應(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)化為圖的劃分問題。圖中的一個(gè)結(jié)點(diǎn)對(duì)應(yīng)一個(gè)數(shù)據(jù)點(diǎn)，結(jié)點(diǎn)間的邊權(quán)值對(duì)應(yīng)數(shù)據(jù)間的相似度。圖的切割就是找到一條或幾條邊，以這些邊為切割點(diǎn)，將圖劃分為兩個(gè)或多個(gè)子圖，使同一個(gè)子圖數(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(zhǎng)最小特征值對(duì)應(yīng)的特征向量。將二分類問題推廣到k類分類，譜聚類可表示為以下優(yōu)化問題

(9)

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

本文利用Three-circle 數(shù)據(jù)集分析、對(duì)比K均值與譜聚類的效果(圖4)。可見：當(dāng)相同標(biāo)記的樣本為同一種顏色，不同標(biāo)記的樣本為不同顏色時(shí)，分類精度更高；K均值的分類精度較低，這是由于其基于原型數(shù)據(jù)歐式距離進(jìn)行分類，只考慮數(shù)據(jù)之間的歐式距離，因此對(duì)復(fù)雜結(jié)構(gòu)數(shù)據(jù)集無法得到滿意的分類結(jié)果(圖4左)；譜聚類將原始數(shù)據(jù)映射到特征向量的數(shù)據(jù)空間中進(jìn)行劃分，因此對(duì)復(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)要求，按照一定大小的時(shí)窗提取每道地震記錄；

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

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

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

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

2.1 模型試驗(yàn)

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

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

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

通過直觀對(duì)比、分析發(fā)現(xiàn)：①K均值聚類方法在原始數(shù)據(jù)空間計(jì)算距離，對(duì)于復(fù)雜數(shù)據(jù)劃分能力有限，地震相劃分結(jié)果(圖7c、圖7d、圖8c、圖8d、圖9c、圖9d)與真實(shí)地震相(圖7a、圖8a、圖9a)差距較大，而且隨著聚類簇?cái)?shù)增加，劃分結(jié)果的隨機(jī)性也逐漸增大，表現(xiàn)為：以地震記錄為輸入的K均值聚類結(jié)果(圖7d、圖8d、圖9d)存在很大的隨機(jī)性，可解釋性差，即使聚類簇?cái)?shù)較低仍存在很大誤差；以地震倒譜特征參數(shù)為輸入的K均值聚類結(jié)果(圖7 c、圖8 c、圖9 c)可以消除地震波形的影響，充分利用反射系數(shù)劃分地震相，可提高劃分精度。②譜聚類方法將數(shù)據(jù)從原始數(shù)據(jù)空間映射到特征向量空間，避免直接計(jì)算數(shù)據(jù)距離，在復(fù)雜數(shù)據(jù)空間中仍有較高的劃分精度，地震相劃分結(jié)果(圖7b、圖8b、圖9b)與真實(shí)地震相(圖7a、圖8a、圖9a)吻合良好。為了說明聚類簇?cái)?shù)為3、6和10的波阻抗的數(shù)值范圍，統(tǒng)計(jì)了每一分類結(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，縱向上分為三層，每層時(shí)間厚度均為100ms，其中上、下兩層為均勻?qū)樱俣葹?000m/s，密度為2.2g/cm3，中間層由Marmousi模型(圖6)縱向疊置組成為了更直觀地定量評(píng)價(jià)模型劃分精度，定義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ù)測(cè)試

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

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

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

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

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

圖13為P區(qū)地震相劃分結(jié)果。由圖可見：①與地震瞬時(shí)振幅地震相劃分結(jié)果(圖13b)相比，譜聚類地震相劃分結(jié)果(圖13a)的地震相邊界更清晰，可解釋性更強(qiáng)。②與地震瞬時(shí)振幅地震相劃分結(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下古地貌圖。對(duì)比圖13a與圖14可見，前者的扇體發(fā)育區(qū)對(duì)應(yīng)于后者的湖底區(qū)域，兩者劃分的邊界吻合良好，驗(yàn)證了譜聚類地震相分析的準(zhǔn)確性。

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

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

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

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

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

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

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