利用三維多項式曲面擬合的子體波形主成分分析剝離火成巖強屏蔽

樂友喜 楊杰飛 姜 蕾 陳學國 吳佳偉 陳藝都

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

0 引言

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

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

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

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

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

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

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

1 三維多項式曲面擬合局部層拉平

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

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

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

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

(1)

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

時窗W的中心點時間擬合多項式為

(2)

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

(3)

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

(4)

其中

(5)

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

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

2 子體波形PCA分解剝離火成巖強屏蔽

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

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

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

Yr×q=Kr×pXp×q

(6)

其中r

(7)

(8)

3 模型測試及實際資料應用

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

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

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

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

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

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

圖4 三維多項式曲面擬合得到的火成巖強反射層在時窗中心點處的時間剖面

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

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

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

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

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

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

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

圖8 對圖7剝離火成巖強屏蔽的結(jié)果

4 結(jié)論

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

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

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

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