樂(lè)友喜 楊杰飛 姜 蕾 陳學(xué)國(guó) 吳佳偉 陳藝都

(①中國(guó)石油大學(xué)(華東)，山東青島 266580；②山東省深層油氣重點(diǎn)實(shí)驗(yàn)室，山東青島 266580；③中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北武漢 430074；④中國(guó)石化勝利油田分公司物探研究院，山東東營(yíng) 257022；⑤中國(guó)石化勝利油田分公司勘探開(kāi)發(fā)研究院，山東東營(yíng) 257015)

0 引言

火成巖是一種特殊的巖性體，具有高速度、非均質(zhì)、原生孔隙度大等特點(diǎn)?；鸪蓭r與圍巖有多種接觸關(guān)系，火成巖反射波場(chǎng)和轉(zhuǎn)換波場(chǎng)較復(fù)雜，致使火成巖發(fā)育區(qū)的地震波場(chǎng)混亂，對(duì)于下伏圍巖的地震響應(yīng)特征具有很強(qiáng)屏蔽作用，導(dǎo)致火成巖下伏儲(chǔ)層反射信號(hào)能量較弱，嚴(yán)重制約了火成巖發(fā)育區(qū)的儲(chǔ)層預(yù)測(cè)。

目前去除地震強(qiáng)屏蔽方法多集中在多子波分解與重構(gòu)技術(shù)。對(duì)于信號(hào)分解，匹配追蹤(Matching Pursuit, MP)算法在信號(hào)的稀疏表示方面效果很好，是多子波分解的核心算法。Mallat等[1]利用基于Gabor原子形成的完備原子庫(kù)，通過(guò)不斷迭代，求取原信號(hào)與原子庫(kù)子波的最佳稀疏表示——MP算法；Chakraborty等[2]首次把MP算法應(yīng)用于地震信號(hào)分析，其中MP原子庫(kù)(即子波類型和子波分辨率)的選擇十分重要。該項(xiàng)技術(shù)主要是利用不同頻率、振幅、相位的各類原子與地震道匹配[3-8]，分解成一個(gè)原子庫(kù)，其中最大能量原子對(duì)應(yīng)于強(qiáng)背景，在重構(gòu)中舍去，從而達(dá)到去屏蔽的效果。

多子波分解與重構(gòu)技術(shù)往往存在以下兩個(gè)問(wèn)題[9]：首先，對(duì)于分解的子波實(shí)際意義不明確，地震剖面上最強(qiáng)子波不一定都代表火成巖，也可能是其他巖體；其次，地震道子波分解是對(duì)單道處理，無(wú)法考慮地層橫向關(guān)系，對(duì)于復(fù)雜的火成巖發(fā)育區(qū)，火成巖存在形態(tài)和厚度變化，反射特征并不一致，容易造成匹配過(guò)度和重構(gòu)剖面連續(xù)性差。

金成志等[9]在2017年首次結(jié)合Wheeler變換和主成分分析(Principal Component Analysis, PCA)技術(shù)去除地震強(qiáng)屏蔽，可了解強(qiáng)屏蔽層之下的砂體展布特征，為去除強(qiáng)屏蔽提供了新思路。長(zhǎng)、短旋回波形特征的主要區(qū)別在于相位和頻率，因此，采用PCA提取反射波形的空間一致性相位信息分離長(zhǎng)、短旋回，可去除地震強(qiáng)屏蔽的影響。由于“長(zhǎng)旋回反射相位特征一致”的認(rèn)識(shí)僅在相對(duì)地質(zhì)年代域中成立，所以在利用PCA分解長(zhǎng)、短旋回波形時(shí)，需要進(jìn)行Wheeler變換。

De Groot等[10]在2006年正式提出三維 Wheeler 域變換技術(shù)；De Bruin 等[11]于2007年提出Wheeler域與時(shí)間域地震數(shù)據(jù)交互解釋的技術(shù)思路。至此，Wheeler域變換技術(shù)被正式推廣至油氣勘探領(lǐng)域。目前，應(yīng)用地震數(shù)據(jù)構(gòu)建Wheeler數(shù)據(jù)體的方法主要有三種[12]：第一種是由De Groot等[10]、Nordlund等[13]提出的基于地震層位解釋的Wheeler數(shù)據(jù)體生成方法，通過(guò)追蹤盡可能多的層位，解釋層序邊界構(gòu)建Wheeler數(shù)據(jù)體。但利用人工拾取同相軸構(gòu)建Wheeler數(shù)據(jù)體效率低、工作量大。第二種方法是Lomask[14]、Fomel[15]提出的將數(shù)據(jù)體拉平進(jìn)而生成Wheeler數(shù)據(jù)體的方法，其優(yōu)勢(shì)是不需要層位拾取，但是在存在不整合面和斷層等的復(fù)雜地質(zhì)情況下，無(wú)法正確拉平同相軸。第三種方法由Stark[16-17]提出，首先通過(guò)地震瞬時(shí)相位解纏生成相對(duì)地質(zhì)時(shí)間(RGT)體，然后在RGT體的基礎(chǔ)上構(gòu)建Wheeler數(shù)據(jù)體。這種方法的精度取決于瞬時(shí)相位解纏的精度，難于取得理想的應(yīng)用效果。

由于火成巖相帶類型多、地震波振幅橫向變化快、反射特征復(fù)雜，因此在利用長(zhǎng)、短旋回波形分析法去除火成巖強(qiáng)屏蔽時(shí)，往往存在以下問(wèn)題[9，12]：①需要對(duì)火成巖發(fā)育區(qū)多個(gè)層序界面進(jìn)行精細(xì)層位解釋，以形成準(zhǔn)確的地層切片體[12]；②若Wheeler域中火成巖反射同相軸無(wú)法拉平，通過(guò)PCA提取的長(zhǎng)旋回波形存在一定偏差，從而影響長(zhǎng)、短旋回分離后的短旋回波形的有效精度；③局部區(qū)域發(fā)育的火成巖爆發(fā)相、噴溢相的界面往往與長(zhǎng)旋回界面不平行，難以準(zhǔn)確剝離和剔除火成巖強(qiáng)反射；④在地下斷層發(fā)育、地層出現(xiàn)尖滅、同相軸連續(xù)性差等復(fù)雜情況下，經(jīng)Wheeler正、反變換后，重構(gòu)的地震信號(hào)可能出現(xiàn)局部失真，并且重構(gòu)過(guò)程較復(fù)雜。

為了剝離火成巖強(qiáng)反射，可根據(jù)解釋的火成巖層位信息設(shè)置相應(yīng)的時(shí)窗寬度，通過(guò)搜索縱、橫向地震波形特征準(zhǔn)確獲得火成巖層位，實(shí)現(xiàn)局部層拉平，再利用子體波形PCA剝離火成巖強(qiáng)反射?；诖耍疚牟捎萌S多項(xiàng)式曲面擬合代替Wheeler變換實(shí)現(xiàn)局部層拉平，形成了基于三維多項(xiàng)式曲面擬合的子體波形PCA火成巖強(qiáng)屏蔽剝離技術(shù)，避開(kāi)了Wheeler變換帶來(lái)的問(wèn)題。同時(shí)，在剝離過(guò)程中無(wú)需引入振幅閾值控制，對(duì)于火成巖反射振幅橫向變化快的問(wèn)題，只要在設(shè)置的子體窗口范圍內(nèi)火成巖反射振幅變化不大，就可以提取該子體窗口的火成巖強(qiáng)反射，并通過(guò)滑動(dòng)子體窗口，在橫向剝離火成巖強(qiáng)反射，提高了方法的實(shí)用性。

1 三維多項(xiàng)式曲面擬合局部層拉平

首先對(duì)解釋的火成巖層位進(jìn)行三維多項(xiàng)式曲面擬合，通過(guò)空間搜索獲取精確的層位，從而有效解決Wheeler變換中的問(wèn)題①；針對(duì)問(wèn)題②和③，如果存在局部發(fā)育的火成巖相帶，無(wú)論其界面是否與本地區(qū)長(zhǎng)旋回界面平行，根據(jù)解釋的火成巖層位信息，設(shè)置相應(yīng)的時(shí)窗寬度，通過(guò)搜索縱、橫向波形特征擬合火成巖的準(zhǔn)確層位，從而實(shí)現(xiàn)局部層拉平，最終基于子體波形PCA分解方法剝離火成巖強(qiáng)反射，從而有效避開(kāi)了問(wèn)題④。具體過(guò)程如下。

利用三維地震數(shù)據(jù)體擬合三維多項(xiàng)式曲線時(shí)，需要選取一個(gè)分析區(qū)域，首先劃分一個(gè)初始窗口，然后假設(shè)某一時(shí)窗W的時(shí)間取值范圍為t∈[-l,l]，橫向窗口范圍為x∈[-m,m]、y∈[-n,n]，對(duì)該時(shí)窗沿x、y方向建立離散坐標(biāo)系

D={(x,y)|x∈[-m,m],

y∈[-n,n]；x,y∈Z}

(1)

式中：x、y分別為Crossline、Inline方向的坐標(biāo)；Z為整數(shù)集合；D為研究區(qū)的二維空間。

時(shí)窗W的中心點(diǎn)時(shí)間擬合多項(xiàng)式為

(2)

式中：T(x，y)為由擬合得到的點(diǎn)(x，y)的強(qiáng)反射同相軸時(shí)間；Ci為第i個(gè)正交多項(xiàng)式的擬合系數(shù)；Gi為D的二元二次正交多項(xiàng)式，即

(3)

根據(jù)多道互相關(guān)性最強(qiáng)的原則確定最佳Ci。對(duì)W內(nèi)數(shù)據(jù)S(x，y，t)可以計(jì)算歸一化多道互相關(guān)系數(shù)，即

(4)

其中

(5)

為了不影響擬合效果，求取其他系數(shù)時(shí)，需要假設(shè)C0的值為0，可以確保擬合的信號(hào)時(shí)間固定在窗口中心。

以C1的掃描為例，選取初始窗口，以初始窗口邊界的樣點(diǎn)時(shí)間作為掃描范圍，通過(guò)不同的樣點(diǎn)時(shí)間得到不同的C1；由式(4)求出相應(yīng)的相關(guān)系數(shù)，其中最大相關(guān)系數(shù)對(duì)應(yīng)的C1即為最終系數(shù)。由于Gi為二元二次正交多項(xiàng)式，故可以用獨(dú)立掃描的方式求得Ci，從而簡(jiǎn)化計(jì)算。具體過(guò)程為：①固定C0(可取為0)，先掃描C1，此時(shí)C2～C5均為0，比較多道互相關(guān)值確定C1；②保持已確定的多道互相關(guān)系數(shù)，再掃描C2，依此類推。掃描過(guò)程中窗口的形狀隨多項(xiàng)式系數(shù)的變化而變化，求取系數(shù)的過(guò)程就是搜索信號(hào)的過(guò)程，由掃描確定的窗口就是期望信號(hào)的窗口。

2 子體波形PCA分解剝離火成巖強(qiáng)屏蔽

PCA的核心思想是將數(shù)據(jù)降維，它能夠有效地提取數(shù)據(jù)中的“主要”信息，降低了復(fù)雜數(shù)據(jù)的維度，從而去除冗余成分，挖掘數(shù)據(jù)中的隱藏信息。PCA方法原理簡(jiǎn)單，已用于眾多領(lǐng)域。因?yàn)镻CA能將一組相關(guān)性較強(qiáng)的變量轉(zhuǎn)化為一組相關(guān)性較弱的矩陣，最大特征值對(duì)應(yīng)的即為主要成分，而小特征值對(duì)應(yīng)的往往是相關(guān)性較差的噪聲。因而PCA常用于去除隨機(jī)噪聲和信息聚類，本文采用基于三維多項(xiàng)式曲面擬合局部層拉平數(shù)據(jù)體，火成巖反射層位經(jīng)過(guò)層拉平，在橫向上波形特征基本一致，通過(guò)PCA將其作為背景信息提取。

具體實(shí)現(xiàn)過(guò)程為：首先以一個(gè)地震道為中心，取若干相鄰道形成一個(gè)PCA單元，再通過(guò)層位時(shí)間信息時(shí)窗約束進(jìn)行層拉平；其次通過(guò)PCA波形分析提取局部層拉平后的相似波形，將此波形從PCA單元中剔除，達(dá)到剝離火成巖強(qiáng)屏蔽的目的。上述基于子體波形的PCA方法稱為子體波形PCA分解方法。具體原理闡述如下。

設(shè)X=(s1，s2，…，sq)為一個(gè)p×q的子窗口內(nèi)的地震數(shù)據(jù)矩陣，其中si為某道地震數(shù)據(jù)，子窗口內(nèi)的道數(shù)為q，每道的時(shí)間樣點(diǎn)數(shù)為p。經(jīng)過(guò)線性變換將原始數(shù)據(jù)X轉(zhuǎn)換到另一個(gè)變量空間Y，Y中的各變量彼此不相關(guān)，即

Yr×q=Kr×pXp×q

(6)

其中r

(7)

(8)

3 模型測(cè)試及實(shí)際資料應(yīng)用

圖1為惠民凹陷YHM地區(qū)火成巖地質(zhì)模型及其正演模擬記錄。由圖可見(jiàn)：火山巖溢流相厚度小，橫向分布范圍廣，產(chǎn)狀近似水平(圖1a)，頂、底形成強(qiáng)反射，波形、振幅基本保持不變(圖1b)；火山巖通道相表現(xiàn)為空白反射或雜亂反射(圖1b)；火山巖爆發(fā)相、噴溢相的橫向厚度變化較大，中間厚，向兩邊逐漸變薄(圖1a)，多個(gè)界面的反射疊合在一起形成復(fù)波，同相軸呈凸起狀，并且淹沒(méi)了下伏地層反射(圖1b)。

圖2為圖1b的Wheeler正變換結(jié)果。由圖可見(jiàn)，Wheeler正變換的時(shí)間段為100～600ms，主要包含5組同相軸。在Wheeler域通過(guò)基于PCA分解方法提取強(qiáng)振幅、相位特征相同的并與長(zhǎng)旋回對(duì)應(yīng)的子體波形，在此基礎(chǔ)上再進(jìn)行Wheeler反變換，得到火成巖強(qiáng)反射(圖3a)，最終得到剔除兩套火成巖強(qiáng)反射的結(jié)果(圖3b)?？梢?jiàn)：①在Wheeler域火山巖溢流相和火山巖沉積相界面下方同相軸被拉平，但火山巖爆發(fā)相、噴溢相反射同相軸并沒(méi)有完全拉平，因此剝離的火山巖強(qiáng)反射存在一定偏差，影響剝離火成巖強(qiáng)屏蔽效果(圖3a)。②火山巖溢流相橫向分布范圍廣，產(chǎn)狀近似水平，因此被很好地剔除，只是在斷層附近有一定殘余；火成巖爆發(fā)相和噴溢相的位置仍存留一定剩余波形，同時(shí)在斷層左側(cè)火成巖不發(fā)育位置提取了假火成巖反射，說(shuō)明僅通過(guò)振幅閾值約束往往難以判定火成巖反射的可靠性，需要解釋人員結(jié)合研究區(qū)地質(zhì)條件分析進(jìn)行人工干預(yù)(圖3b)。

圖1 惠民凹陷YHM地區(qū)火成巖地質(zhì)模型(a)及其正演模擬記錄(b)

圖2 圖1b的Wheeler正變換結(jié)果

圖3 基于長(zhǎng)、短旋回PCA分解法剝離的火成巖強(qiáng)反射(a)及剔除火成巖強(qiáng)反射(b)的結(jié)果

圖4為三維多項(xiàng)式曲面擬合得到的火成巖強(qiáng)反射層在時(shí)窗中心點(diǎn)處的時(shí)間剖面。由圖可見(jiàn)，獲得的火成巖反射層橫向連續(xù)性強(qiáng)，并且在斷層位置處較精準(zhǔn)地提取了火成巖空間位置，在此基礎(chǔ)上可對(duì)火成巖反射層進(jìn)行局部層拉平。

圖4 三維多項(xiàng)式曲面擬合得到的火成巖強(qiáng)反射層在時(shí)窗中心點(diǎn)處的時(shí)間剖面

通過(guò)分析子體波形PCA分解的各分量發(fā)現(xiàn)，在火成巖反射層局部拉平的子體窗口內(nèi)，主要能量集中在最大特征值對(duì)應(yīng)的分量上，且橫向波形特征一致，該分量即為局部拉平后的火成巖強(qiáng)反射，而被火成巖淹沒(méi)的其他反射層信號(hào)分布在其余特征值對(duì)應(yīng)的分量中。因此，選取最大特征值對(duì)應(yīng)的分量進(jìn)行重構(gòu)，可以剝離火成巖強(qiáng)反射(圖5a)，并得到剔除火成巖強(qiáng)反射的剖面(圖5b)?？梢?jiàn)，本文方法剝離的火成巖強(qiáng)反射橫向連續(xù)性強(qiáng)(圖5a)，能夠較準(zhǔn)確地剔除火成巖強(qiáng)反射界面的綜合響應(yīng)，保留火成巖下方的地層反射信息，從而更好地削弱火成巖對(duì)下伏地層反射波的屏蔽作用(圖5b)。

圖5 基于多項(xiàng)式曲面擬合子體波形PCA 分解方法剝離的火成巖強(qiáng)反射(a)及剔除火成巖強(qiáng)反射(b)的結(jié)果

為了驗(yàn)證基于多項(xiàng)式曲面擬合子體波形PCA 分解方法的實(shí)際應(yīng)用效果，對(duì)惠民凹陷YHM地區(qū)的實(shí)際地震資料進(jìn)行處理。該地區(qū)火成巖十分發(fā)育，火成巖頂、底界面較大的波阻抗差異導(dǎo)致出現(xiàn)強(qiáng)反射同相軸，強(qiáng)反射同相軸嚴(yán)重屏蔽了火成巖下方儲(chǔ)層的弱信號(hào)，降低了儲(chǔ)層預(yù)測(cè)精度。

W1井綜合錄井圖(圖6a)顯示：W1井T3界面上方火成巖為玄武巖，T3界面下方沙三段3砂組儲(chǔ)層發(fā)育。W3井綜合錄井圖(圖6b)顯示：W3井T3界面上方火成巖為凝灰?guī)r，T3界面下方沙三段3砂組儲(chǔ)層不發(fā)育。由W1-W2-W3井連井地震剖面(圖7)可見(jiàn)，T3反射層(圖7綠色虛線)上方火成巖發(fā)育，為中強(qiáng)連續(xù)反射(圖7藍(lán)框位置)，主要目的層沙三段3砂組(圖7黃色虛線)在T3反射層下方50ms處。

圖6 W1井(a)、W3井(b)綜合錄井圖

圖8為對(duì)圖7剝離火成巖強(qiáng)屏蔽的結(jié)果。由圖可見(jiàn)：①多子波分解法剝離的火成巖強(qiáng)反射(圖8a)橫向連續(xù)性較差，其使用的子波庫(kù)為雷克子波庫(kù)，剔除火成巖強(qiáng)反射剖面(圖8b)殘留部分火成巖信息。②基于Wheeler正、反變換的長(zhǎng)、短旋回PCA分解法在斷層位置處剝離了火成巖強(qiáng)反射(圖8c)，并錯(cuò)誤地識(shí)別了火成巖范圍(藍(lán)框區(qū)域)，導(dǎo)致剔除的火成巖強(qiáng)反射不準(zhǔn)確(圖8d)；基于長(zhǎng)、短旋回PCA分解法無(wú)法準(zhǔn)確提取局部火成巖同相軸(紅框區(qū)域)，如果擴(kuò)大提取的時(shí)間范圍，會(huì)損失有效信息。③三維多項(xiàng)式曲面擬合得到的火成巖強(qiáng)反射層在時(shí)窗中心點(diǎn)處的時(shí)間剖面(圖8e)在橫向上光滑、連續(xù)，對(duì)應(yīng)層位準(zhǔn)確。④采用基于三維多項(xiàng)式曲面擬合的子體波形PCA分解法剝離的火成巖強(qiáng)反射考慮了火成巖發(fā)育區(qū)的地震反射特征，實(shí)現(xiàn)了局部層拉平，再進(jìn)行子體波形PCA分解，準(zhǔn)確地提取了火成巖強(qiáng)反射(圖8f)，剔除火成巖強(qiáng)反射剖面(圖8g)的橫向連續(xù)性強(qiáng)，較準(zhǔn)確地剔除了強(qiáng)反射界面的影響，從而在一定程度上增強(qiáng)了下伏地層反射。⑤利用基于時(shí)頻分析的地震波能量衰減補(bǔ)償方法[18]對(duì)圖8g進(jìn)行高頻能量衰減補(bǔ)償處理(圖8h)，可提高對(duì)地震弱反射信息的刻畫能力。圖9為圖7藍(lán)框區(qū)域及圖8h藍(lán)框區(qū)域的局部放大，可見(jiàn)：W1井T3界面上方的火成巖強(qiáng)反射(圖9a紅圈處)淹沒(méi)了沙二段1砂組以及1砂組與火成巖之間的一套薄互層砂組(圖9a紅色虛線框)，導(dǎo)致在原始地震剖面上無(wú)法識(shí)別這兩套砂組；在剔除火成巖強(qiáng)反射及能量衰減補(bǔ)償處理后，這兩套砂組對(duì)應(yīng)的同相軸均呈弱反射信號(hào)特征，能夠進(jìn)行橫向追蹤和識(shí)別(圖9b紅色虛線框)，同時(shí)準(zhǔn)確地刻畫了3砂組，砂體邊界向右延伸了175m(圖9中藍(lán)色箭頭處)。實(shí)際應(yīng)用表明：本文方法能夠削弱火成巖對(duì)下伏地層反射的屏蔽作用，提高對(duì)地震弱反射信息的刻畫能力。

圖7 W1—W2—W3井連井地震剖面

圖8 對(duì)圖7剝離火成巖強(qiáng)屏蔽的結(jié)果

4 結(jié)論

本文采用三維多項(xiàng)式曲面擬合實(shí)現(xiàn)火成巖強(qiáng)反射局部層拉平，結(jié)合子波相位特性，利用子體波形PCA分解方法剝離火成巖強(qiáng)反射，避開(kāi)了Wheeler變換帶來(lái)的問(wèn)題，形成了一種非常適用的火成巖強(qiáng)屏蔽剝離技術(shù)。

(1)該項(xiàng)技術(shù)在剝離火成巖強(qiáng)反射過(guò)程中不需要引入振幅閾值的控制，通過(guò)分析研究區(qū)火成巖反射特征，合理設(shè)置子體窗口，有效解決了火成巖反射振幅橫向變化快的問(wèn)題；同時(shí)，通過(guò)子體窗口的連續(xù)滑動(dòng)，解決了火成巖強(qiáng)反射橫向剝離問(wèn)題，克服了多子波分解方法和基于Wheeler變換的長(zhǎng)、短旋回PCA分解方法橫向連續(xù)性差的缺陷。

(2)該項(xiàng)技術(shù)在實(shí)際應(yīng)用過(guò)程中，通過(guò)對(duì)解釋的火成巖層位進(jìn)行三維多項(xiàng)式曲面擬合，利用空間搜索獲取精確層位，自動(dòng)生成局部層拉平數(shù)據(jù)體，為子體波形PCA分解提供了可靠的基礎(chǔ)數(shù)據(jù)，簡(jiǎn)化了實(shí)現(xiàn)步驟，實(shí)用性大大增強(qiáng)。

(3)模型測(cè)試和實(shí)際資料應(yīng)用表明，該方法剝離的火成巖強(qiáng)屏蔽橫向連續(xù)性強(qiáng)，能夠較準(zhǔn)確地剔除火成巖強(qiáng)反射界面的綜合響應(yīng)，在一定程度上削弱了火成巖對(duì)下伏地層反射波的屏蔽作用，提高了對(duì)地震弱反射的識(shí)別能力，為弱信號(hào)的能量補(bǔ)償和后續(xù)處理奠定了基礎(chǔ)。