基于迭代反演的層間多次波壓制方法

包培楠， 王孝， 謝俊法， 王維紅， 石穎， 徐嘉亮， 楊育臣

1 東北石油大學地球科學學院， 大慶 163318 2 中國石油天然氣股份有限公司勘探開發研究院西北分院， 蘭州 730070 3 東北石油大學環渤海能源研究院， 秦皇島 066004 4 東北石油大學非常規油氣研究院， 大慶 163318 5 “陸相頁巖油氣成藏及高效開發”教育部重點實驗室， 大慶 163318 6 密蘇里科技大學， 羅拉 65401， 美國

0 引言

層間多次波是地震資料中的規則干擾波，其動校正時差和有效波差別較小，難以識別和有效壓制，對有效波的成像和識別帶來很大的困難(Berkhout, 2006)，同時層間多次波的存在和其他噪聲一樣，對后續地震資料的解釋、反演和應用都會帶來不利影響，導致構造和油氣識別的精度下降，鉆探成功率降低，對油氣的勘探開發帶來極大的風險(孫小東等，2020).因此，有效壓制層間多次波，進而提高地震資料成像的精度，是油氣地震勘探領域的一個重要研究方向.

多次波壓制方法主要分為兩類，一類是濾波方法，另一類是預測相減方法.濾波方法主要利用一次波與多次波的時差關系和周期特征進行多次波的識別和壓制(Taner, 1980; Lokshtanov, 1999)，如拋物Radon濾波(Hampson, 1986)和雙曲Radon濾波(Foster and Mosher, 1992; 石穎和王維紅, 2012; Shi et al., 2017)等，該類方法計算效率高，算法容易實現，當有效波和多次波之間的動校正時差較大時，可取得滿意的多次波壓制效果.但對于復雜介質，如速度梯度較小(或速度反轉)的情況，或構造變化劇烈的介質，應用濾波法難以有效識別有效波和多次波，往往得不到理想的多次波壓制結果(Berkhout and Verschuur, 1997).預測相減方法基于波動理論，能更好的適應于復雜介質的情況，諸多地球物理工作者都對該類方法進行了系統深入的研究.該類多次波壓制方法包括逆散射級數法(Araujo et al., 1994; Weglein et al., 1997, 2003; Otnes et al., 2004; Kristopher, 2017)、共聚焦點技術(Berkhout and Verschuur, 2006)以及近年來發展的基于虛同相軸的層間多次波壓制方法(吳靜等, 2013)和基于Marchenko的層間多次波壓制方法(Behura et al., 2014; Verschuur and Berkhout, 2015；匡偉康等，2018).上述方法中，基于散射理論的逆散射級數法無需地下的先驗信息，預測一次可以得到與所有界面相關的同階層間多次波，在沒有有效手段區分有效波和多次波時，是層間多次波壓制較為有效的方法，但是該方法存在計算量大，預測遠偏移距多次波效果差的缺點.共聚焦點方法適合復雜介質條件的多次波壓制，但該方法需要全波場數據，當地震資料的近偏移距缺失時需要相應的波場重建方法進行數據重建(王錦妍等，2020)，且一次只能預測得到與某一界面相關的層間多次波，在一定程度上依賴初始速度模型以求取準確的聚焦算子.虛同相軸方法能夠較準確地預測層間多次波，但對觀測系統要求較高，尚未擺脫對人工操作的依賴.Marchenko層間多次波壓制方法在模型數據上取得了很好的效果，也可應用于實際數據(Zhang and Slob, 2020)，但是對于層間多次波發育且信噪比較低的陸上地震資料多次波壓制還存在計算不穩定的情況.

Jakubowicz(1998)提出了一種數據驅動構造層間多次波的新方法，劉戰等 (2019) 對其進行了詳細推導，巧妙地將地下散射點移到表面，但是對于實際數據無法進行精確的層間多次波壓制.Ikelle(2006)、Ikelle等(2009)等通過引入虛源點的概念較為有效地解決了該方法應用于實際地震資料層間多次波壓制的問題.Ikelle認為通過一步預測相減就可有效壓制某一界面產生的所有層間多次波，但是該方法對多次波振幅預測精度不夠.提高多次波壓制效果的方法主要有兩種方法，一是增大自適應相減的濾波算子長度以校正多次波的振幅和相位，另一種是提高多次波的預測精度.第一種方法對有效波的振幅產生較大傷害，針對提高多次波預測精度，在數據驅動和CFP層間多次波壓制方法的基礎上，提出一種基于迭代反演的層間多次波壓制方法(MSI)，該方法直接利用地震數據本身進行層間多次波預測，將多次波壓制看做優化問題，避免構建CFP道集所需的聚焦運算和地下速度的估計，該方法的優點有兩個：一是計算效率高，二是多次波預測精度高.理論模型數據測試表明該方法具有高精度的特點.

1 方法原理

基于迭代反演的層間多次波壓制方法，可以看作是自由表面多次波壓制方法(SRME，Surface-Related Multiple Elimination)的擴展.在CFP層間多次波壓制方法的基礎上，首先借助于數據矩陣表示法和反饋迭代模型(WRW)，建立反射波的數學模型.然后利用“反饋迭代模型”表示出一次波與多次波之間存在的關系，預測出層間多次波.最終利用自適應匹配濾波法將預測的多次波從地震數據中減去.

1.1 CFP層間多次波產生及壓制

假定在地表z0處激發震源產生地震波，如圖1a所示，地震波在地下第zm層處發生反射，置于地表的檢波器接收一次反射波.依據圖1b所示的WRW模型，檢波器接收到的一次反射波地震記錄可表示為

ΔP(z0,z0)=D(z0)ΔX(z0,z0)S(z0)

=D(z0)ΔP-(z0,z0)，

(1)

其中:

ΔP-(z0,z0)=ΔX(z0,z0)S(z0)，

(2)

(3)

式中，ΔP(z0,z0)為只包含一次反射波的地震數據矩陣，D(z0)為檢波點矩陣，ΔX(z0,z0)為一次反射波的傳遞矩陣，S(z0)為震源矩陣，ΔP-(z0,z0)為反射波場，W-(z0,zm)為上行傳播矩陣，W+(zm,z0)為下行傳播矩陣，R(zm,zm)為地下第zm層反射系數矩陣，m為地下反射層序號，(z0,z0)表示地表激發地表接收.

若地震波在地表z0處發生至少一次下行反射后，被置于地表的檢波器所接收，其生成表面多次波的射線路徑如圖2a所示，依據圖2b所示包含表面多次波的WRW模型，則地震記錄可表示為

P(z0,z0)=D(z0)ΔX(z0,z0)[S(z0)+R(z0,z0)

×P-(z0,z0)]，

(4)

式中，P-(z0,z0)表示在不考慮檢波點對地震波場影響情況下的含表面多次波的反射波場，R(z0,z0)表示向下延拓算子.則有:

P-(z0,z0)=ΔP-(z0,z0)+ΔX(z0,z0)[R(z0,z0)

×P-(z0,z0)]，

(5)

P(z0,z0)=D(z0)P-(z0,z0)，

(6)

由于總波場P(z0,z0)含一次波和表面多次波兩部分，根據公式(4)，表面多次波{δM(z0,z0)}0可表示為

{δM(z0,z0)}0=D(z0)ΔX(z0,z0)[R(z0,z0)

×P-(z0,z0)].

(7)

引入界面算子A(z0,z0)，假定其與震源矩陣、檢波點矩陣以及向下延拓算子有關，且存在關系式：

A(z0,z0)=S-1(z0)R(z0,z0)D-1(z0)，

(8)

假定在自由表面處，向下延拓算子R(z0,z0)可近似為-I(I為單位矩陣)，則有：

A(z0,z0)≈-[D(z0)S(z0)]-1，

(9)

若等式(7)右邊乘以一個單位矩陣I=-A(z0,z0)[D(z0)S(z0)]，結合公式(1)和公式(5)，表面多次波又可表示為

{δM(z0,z0)}0=ΔP(z0,z0)A(z0,z0)P(z0,z0)，

(10)

從總波場中去除表面多次波，可得到壓制多次波后的地震波場，即：

圖1 一次波傳播路徑及其WRW模型(a) 一次波傳播路徑； (b) 產生一次波的WRW模型.Fig.1 Propagation path of primary waves and WRW model(a) Propagation path of primary waves; (b) WRW model of primary wave generation.

ΔP(z0,z0)=P(z0,z0)-{δM(z0,z0)}0.

(11)

若令X(z0,z0)為含一次波和多次波的傳遞矩陣，根據公式(1)，包含一次波和多次波的地震記錄P(z0,z0)可表示為

P(z0,z0)=D(z0)X(z0,z0)S(z0).

(12)

假定地震波在地面以下的zn界面發生至少一次下行反射后，被置于地表的檢波器所接收，產生層間多次波的射線路徑如圖3a所示，WRW模型如圖3b所示，則只包含一次反射波和層間多次反射波的地震記錄可表示為

{P(z0,z0)}0=D(z0){X(z0,z0)}0S(z0)，

(13)

式中，{}0表示與界面z0相關的多次波都已被壓制，{P(z0,z0)}0表示與表面相關的多次波壓制后的地震記錄，{X(z0,z0)}0表示與表面相關的多次波壓制后的傳遞矩陣.對于地表以下的界面，做類似定義：

{P(z0,z0)}n=D(z0){X(z0,z0)}nS(z0)，

(14)

式中，n=0,1,2,…,∞，當n=0時，表示與自由表面相關的情況.{}n表示關于界面z≤zn的多次波都已被壓制，只有關于界面z>zn的層間多次波存在.

CFP方法是將表面多次波壓制方法引入到層間多次波壓制中，該方法假設與層z≤zn-1相關的所有多次波在之前都已被壓制，得到向下外推炮記錄{P(zn,z0)}n-1作為壓制層間多次波程序的輸入數據，表達式為

{P(zn,z0)}n-1=Γ(zn,z0){P(z0,z0)}n-1，

(15)

式中，{P(zn,z0)}n-1表示震源在地表z0處，檢波器在地下zn處的地震記錄.Γ(zn,z0)表示檢波器一側的向下延拓算子.根據公式(3)，關于界面zn的傳遞矩陣為

(16)

根據公式(7)，與界面zn相關的層間多次波為

{δM(z0,z0)}n=

(17)

圖2 表面多次波傳播路徑及其WRW模型(a) 表面多次波傳播路徑； (b) 產生表面多次波的WRW模型.Fig.2 Propagation path of surface multiples and WRW model(a) Propagation path of surface multiples; (b) WRW model of generation of surface multiples.

圖3 層間多次波傳播路徑及其WRW模型(a) 邊界zn處產生層間多次波的波傳播路徑； (b) 一次波(m>n)和層間多次波正演的WRW反饋模型.Fig.3 Propagation path of internal multiples and WRW model(a) Propagation path of internal multiples at boundary zn; (b) WRW feedback model of forward modeling for primary waves (m>n) and internal multiples.

{δM(z0,z0)}n=

(18)

A(zn,zn)=S-1(zn)R(zn,zn)D-1(zn)

≈-[D(zn)S(zn)]-1，

(19)

則壓制與該界面相關的層間多次波公式為

{P(z0,z0)}n={P(z0,z0)}n-1-{δM(z0,z0)}n

×{P(zn,z0)}n-1.

(20)

通過以上分析可知，層間多次波壓制算法與表面多次波壓制算法相似，只是對于層間多次波而言，需要由檢波器在深度zn時的炮記錄來代替檢波器位于表面時的炮記錄.

1.2 MSI層間多次波壓制算法

CFP層間多次波壓制方法需要求取地震數據的延拓算子，在很大程度上增加了層間多次波預測過程的計算量.為了降低計算成本，由前面闡述的CFP方法入手，假設經過k次迭代壓制層間多次波后的地震記錄{P(z0,z0)}′k可以由上一次迭代結果{P(z0,z0)}′k-1得到，引入一個類似于界面算子A(zn,zn)的卷積因子T(k-1)，則經過k次迭代的層間多次波{M(z0,z0)}′k表示為

{M(z0,z0)}′k=T(k-1)[P(z0,z0)-{P(z0,z0)}′k-1]，

(21)

由公式(20)可知，界面算子A(zn,zn)為

(22)

將公式(22)代入公式(20)，得到：

×({P(z0,z0)}n-1-{P(z0,z0)}n).

(23)

在實際應用中，矩陣的逆可以根據最小平方的形式進行求取(Berkhout，2006)，則公式(23)又可表示為

×({P(z0,z0)n-1-{P(z0,z0)}n),

(24)

公式(24)中{δM(z0,z0)}n表示關于界面zn的層間多次波，若考慮所有的層間多次波，公式(24)可寫為

{M(z0,z0)}′k={P(z0,zn)}′k({P(z0,z0)}′k-1)H

×[{P(z0,z0)}′k-1({P(z0,z0)}′k-1)H]-1

×(P(z0,z0)-{P(z0,z0)}′k),

(25)

式中，H表示共軛轉置，則卷積T(k-1)因子為

T(k-1)={P(z0,z0)}′k-1({P(z0,z0)}′k-2)H

×[{P(z0,z0)}′k-2({P(z0,z0)}′k-2)H]-1，

(26)

最終壓制的多次波結果為

{P(z0,z0)}′k=P(z0,z0)-{M(z0,z0)}′k.

(27)

由以上可知，在MSI算法中，不需要顯式的表面算子和顯式的震源矩陣，表面算子被原始數據P(z0,z0)及先前兩步的多次波壓制結果{P(z0,z0)}′k-1和{P(z0,z0)}′k-2所替代.MSI算法在多次波預測中隱含的考慮了表面算子的空間變化，并在多次波預測不斷迭代反演更新中實現層間算子的空間變化，因此在迭代反演過程中，未被預測或者預測不完整的層間多次波將隨迭代次數的增加逐漸完善，同時在一定程度上解決了后續多次波自適應相減方法中的非線性問題.

在實際應用中，將去除直達波的地震數據輸入到程序中進行迭代以求取卷積因子并對多次波模型進行更新，根據數據的具體情況及計算效率的綜合考慮選擇合適的迭代次數(一般迭代2～3次)得到最終的多次波模型，然后利用自適應匹配濾波進行相減.其實現流程如圖4所示.

圖4 MSI層間多次波壓制算法實現流程圖Fig.4 Implementation flowchart of MSI internal multiple suppression algorithm

2 模型試算

為了證明上文所述方法的實用性和有效性，采用兩個不同的地質模型進行測試.模型一為含楔狀構造的地質模型，模型二是針對我國西北某地區的地質特征，通過測井資料插值的方式生成速度模型正演的數模數據.

2.1 模型一

首先對圖5所示模型進行測試，其中圖5a為速度模型，包含三個強阻抗界面，第三層為高速層，地震波遇到層界面將產生較強的層間多次波.圖5b為密度模型，第三層為高密度層.模型大小為4500 m×1400 m，采用的觀測系統為中間放炮兩邊接收，炮點和檢波點均勻分布于地表，炮間距和道間距均為20 m，震源激發201炮，每炮201個檢波器接收.激發震源采用主頻為25 Hz的雷克子波，采樣間隔為4 ms，地震反射記錄為2.5 s.為了避免表面多次波的影響，地震波場正演計算時模型四周都采用了吸收邊界.

圖6為第101炮地震波場正演模擬記錄，黑色箭頭指示的同相軸分別是三個反射界面的一次波，其余均為層間多次波.圖7為預測的層間多次波，相比迭代1次，經過迭代3次預測的層間多次波信息更完整，其相位和能量都更與實際的層間多次波相符，尤其在1.3～1.5 s之間，經過3次迭代后，這個時間段的層間多次波基本上都被預測出.注意在層間多次波預測過程中，為了避免產生噪聲而傷害有效波，會設置時窗來控制多次波預測范圍.圖8為經過3次迭代壓制層間多次波后地震記錄，其中圖8a為迭代1次壓制層間多次波后的地震記錄，圖8b為迭代3次壓制層間多次波后的地震記錄，觀察可知，迭代1次后，層間多次波得到部分壓制，一次反射波的原始波場特征完好；迭代3次即獲得較為理想的層間多次波壓制效果，多次波殘余能量很少，同時一次波的振幅和相位信息得到較好的保護.圖9為3次迭代層間多次波壓制后的零偏移距剖面，與原始模型地質特征相符，層位特征明顯，無虛假構造出現.

圖5 (a) 速度模型； (b) 密度模型Fig.5 (a) Velocity model； (b) Density model

圖6 含層間多次波的地震記錄Fig.6 Seismic record with internal multiples

圖10為單道數據進行層間多次波壓制前后的對比圖，黑色虛曲線為原始單道數據，綠色實曲線為壓制層間多次波后的，紅色實曲線為預測的層間多次波.從圖中可以看到預測的層間多次波與實際的層間多次波在到時、相位上都有很好的一致性，只在振幅上有微弱差異，此差異可利用后續的自適應匹配相減來消除.黑色與紅色實曲線的對比也可以看出自適應匹配相減后層間多次波得到了有效的壓制，且不傷害有效波，說明迭代反演的層間多次波壓制方法對模型數據的壓制效果較好.

圖7 預測的層間多次波(a) 迭代1次； (b) 迭代3次.Fig.7 Predicted internal multiples(a) By one iteration; (b) By three iterations.

2.2 模型二

圖11為通過層位數據和測井資料控制的方式得到西北地區研究區的速度模型，速度模型大小為9300 m×3200 m，共有17個波阻抗分界面，圖中白色箭頭所指向的四個層為低速薄煤層，煤層速度在2800 m·s-1左右，厚度約10 m，與上下高速地層(4000 m·s-1左右)形成強反射界面，層間多次波較發育，因此層間多次波與一次波的差異很小且能量較強，壓制時很容易傷害有效波，壓制難度較大.

依據速度模型，采用有限差分正演方法模擬地震記錄，圖12為炮數據層間多次波壓制前后對比圖，每炮400道，道間距為10 m，炮檢距為10 m，采樣間隔為0.2 ms，模型的網格間距在水平方向和垂直方向上均為2 m，由于MSI方法是一個迭代反演的過程，直達波的存在會影響層間多次波的壓制效果，在壓制層間多次波之前首先需去除直達波.圖12a、b為壓制層間多次波前后的單炮記錄對比，觀察圖12a發現，層間多次波主要分布在2.7～4.0 s之間.從圖12b易知，在反射時間3.0～4.0 s之間層間多次波得到明顯壓制，而有效波得到較好地保護.對層間多次波壓制前后的數據進行頻譜分析，如圖12c、d所示，主頻范圍稍有拓寬，拓寬5 Hz左右，能量整體有所抬升，這在一定程度上提高了分辨率，說明該方法未傷及有效波，只有層間多次波被有效壓制(圖12中白色箭頭所示都為層間多次波).

圖13為圖12a、b實際數模單炮記錄的局部放大圖，對比發現，2.7～4.0 s之間的層間多次波得到有效壓制.

圖14為層間多次波壓制前后疊加剖面對比圖，從圖中可以看出層間多次波壓制效果較為明顯，為了能夠更清楚的看出層間多次波壓制效果，將整個地震剖面(圖14a、b)上白色框所示部分放大，分別得到圖14c、d.從整個疊加剖面及局部放大圖中可以進一步看出，3.0～4.0 s之間的層間多次波基本被壓制，壓制效果較好.

3 結論

在共聚焦點層間多次波壓制方法基礎上，引入迭代反演策略，將層間多次波壓制的計算轉化為反演過程，用卷積因子代替CFP算法中的聚焦算子，形成基于迭代反演的層間多次波壓制方法，進而應用該方法實現數據驅動的高精度層間多次波預測.MSI方法通過將地下散射點移至地表，使得該方法的實用性大大增強，而且是完全數據驅動，不需要地下介質的任何先驗信息，層間多次波預測精度較高，MSI方法的試算和應用表明，一般通過三次迭代即可得到高精度的多次波預測結果，具有很高的計算效率.采用自適應匹配濾波方法將預測得到的層間多次波進行壓制，可得到多次波壓制結果.理論模擬算例表明MSI方法數據適應性強、計算精度高、計算成本低，具有良好的實際地震資料層間多次波壓制的應用前景.

圖8 層間多次波壓制(a) 迭代1次； (b) 迭代3次.Fig.8 Internal multiples suppression(a) By one iteration; (b) By three iterations.

圖9 經3次迭代壓制層間多次波的零偏移距剖面(a) 原始零偏移距剖面； (b) 迭代3次.Fig.9 Zero-offset profiles after internal multiples suppression by three iterations(a) Original zero-offset profile; (b) After three iterations.

圖10 單道數據層間多次波壓制前后對比圖Fig.10 Comparison of single-trace data before and after internal multiples suppression

圖11 速度模型Fig.11 Velocity model

圖12 層間多次波壓制前后對比圖(a) 層間多次波壓制前的單炮記錄； (b) 層間多次波壓制后的單炮記錄； (c) 層間多次波壓制前頻譜； (d) 層間多次波壓制后頻譜.Fig.12 Comparison before and after internal multiples suppression(a) Shot gather before internal multiples suppression; (b) Shot gather after the internal multiples suppression; (c) Spectrum before internal multiple suppression; (d) Spectrum after internal multiple suppression.

圖13 單炮部分數據層間多次波壓制前后對比圖(a) 層間多次波壓制前； (b) 層間多次波壓制后.Fig.13 Comparison of part shot gather before and after internal multiples suppression(a) Before multiples suppression; (b) After internal multiple suppression.

圖14 疊加數據層間多次波壓制前后對比圖(a) 層間多次波壓制前； (b) 層間多次波壓制后； (c) 部分層間多次波壓制前； (d) 部分層間多次波壓制后.Fig.14 Comparison of stacked data before and after multiples suppression(a) Stacked data before internal multiples suppression; (b) Stacked data after multiples suppression; (c) Part stacked data before multiples suppression; (d) Part stacked data after internal multiple suppression.