基于粒子群優化和保幅成像條件的廣義屏偏移

何 潤 楊午陽 王恩利 魏新建 謝春輝 閆國亮

(中國石油勘探開發研究院西北分院，甘肅蘭州 730020)

0 引言

現今地震勘探已不僅僅關注構造形態的準確性，更力求得到可靠的巖性、物性等信息，這就要求偏移成像技術在得到足夠高的成像精度的同時還具有一定振幅保真能力。疊前深度偏移成像技術已由最初的Kirchhoff算法發展到單程波偏移算法和逆時偏移算法。其中的單程波算法對復雜構造具有一定成像能力，且運算效率較高等特點，因此近年來有較快發展。

自Stoffa等[1]提出裂步傅里葉法(Split-step Fourier，SSF)和Ristow等[2]提出傅里葉有限差分法(Fourier finite-difference，FFD)這兩種對橫向速度變化介質有一定適應能力的偏移算法之后，頻率—波數域單程波偏移算法逐漸受到關注。朱遂偉等[3]采用模擬退火算法全局優化FFD算法，提高了對陡傾地層的成像精度。梅金順等[4]通過引入ω循環型邊界條件，結合正則化方法，克服了FFD算法的邊界效應并加快收斂速度。黃建平等[5]和朱峰等[6]通過推導黏聲VTI介質中的FFD波場延拓算子，進一步擴大了FFD算法的適用范圍。另外，de Hoop等[7]和Le Rousseau等[8]通過對相位屏理論的研究，提出了廣義屏偏移算子(Generalized screen propagator，GSP)，它也能準確描述地震波在陡傾地層或復雜構造中的傳播情形。此后，Shin等[9]將廣義屏算子應用于TTI介質構造成像，Kim等[10]將該偏移算子引入彈性介質多波多分量地震勘探，均拓展了廣義屏算子的實際應用范圍。

另一方面，隨著勘探目標從過去的構造油氣藏轉向更加隱蔽的巖性、地層油氣藏，甚至是致密油氣或頁巖氣領域，地震波的振幅信息在地震資料精細解釋和儲層預測工作中也越來越重要，使波動方程疊前深度偏移方法也向振幅保真的方向發展，力求能更精確地定量描述儲層巖性、物性、流體性質等。張關泉[11]通過分別推導上、下行波波動方程，得出在非均勻介質中符合地震波傳播動力學的上、下行波耦合方程組。張宇[12]根據從雙程波方程中分裂得出的上、下行單程波方程組，在波場延拓中求解真振幅共炮集波動方程。在此基礎上，劉定進等[13-14]分別對傅里葉有限差分算子和廣義屏算子進行真振幅推導，并且修改邊界條件，從而得到保幅的疊前偏移結果。呂彬等[15]探討了保幅型SSF疊前深度偏移方法的理論及應用；崔興福等[16]實現了非均勻介質中波動方程混合法保幅偏移；葉月明等[17]通過合成平面波的方法顯著提高了單程波偏移算法的計算效率。

波場延拓和成像條件是偏移成像的兩個部分，二者需相互匹配才能得到準確的地質構造圖像。對于真振幅的偏移成像，不僅需要在波場延拓時保持能量的穩定，還需要在保幅的成像條件下進行，這樣才能保留更完整的巖性、物性以及流體性質等儲層信息。目前對保幅的成像條件研究較少，在常用的成像條件互相關和反褶積中，僅反褶積成像條件具備部分保幅能力。但由于它為除法運算，所以存在成像不穩定的現象。針對反褶積成像的不穩定性，學者們利用不同的方法提高了該算法的實用性。如Guitton等[18-19]通過對算式中分母項添加阻尼因子以確保結果的穩定，而葉月明等[20]則采用平滑算法來保持計算穩定。

本文從成像精度和保幅能力兩個方面對常規保幅高階廣義屏算法進行改進： 首先利用粒子群優化算法對廣義屏算子中近似展開式的常系數項做最優化求解，提高廣義屏算法的成像角度(精度)；再利用高斯窗函數反褶積成像條件進行成像，有效提高了偏移成像的保幅效果。

1 高階廣義屏保幅偏移算法

根據張關泉[11]和張宇[12]的推導，二維保幅單程聲波方程為

(1)

(2)

式中：pU、pD分別為上行、下行波場；ω為圓頻率；Λ和Γ在頻率—波數域中的表達式分別為

式中：kx為(沿x方向的)水平波數；kz為(沿z方向的)垂直波數。

在保幅單程波方程中，地震波既能保持其運動學特征，又能保留其動力學特征。對Λ的求解情況決定構造形態的準確性；對Γ的求解情況決定地震波振幅的保真效果。目前對Λ和Γ的求解均是對其中的變量進行泰勒展開近似計算，其具體表達式[14，21-22]分別為

(3)

(4)

2 粒子群優化算法

2.1 廣義屏算法成像精度的探討

(5)

將該多項式系數項設為待定系數做重新擬合逼近，其展開式可表示為如下形式

(6)

式中aj為待定系數。

根據函數最佳平方逼近理論，求解多元函數方程的最優值，即可寫為

(7)

其中

由于垂直波數算子kz具有實際物理意義，表示在波場延拓中隨延拓深度的遞增而反復迭代運算得到的波場值，所以求取全局最優解才能得到更準確的波場值，也即取得更好逼近效果。

2.2 粒子群優化算法原理

粒子群優化算法(Particle Swarm Optimization,PSO)是一種原理簡單易懂、算法易于實現的智能尋優算法，自Kennedy等[23-24]提出后便廣泛應用于眾多領域。目前在地球物理相關領域的應用大多為反演和數據處理方面，而在成像偏移上涉及較少。

其基本原理簡述為：在第k次迭代時，m維的第i個粒子具有速度vi(t)和位置xi(t)，而該粒子下一時刻的速度vi(t+1)和位置xi(t+1)將由該粒子當前的速度和位置、當前最優位置Pi(t)和全局最優位置Pg(t)共同決定，通過反復迭代最終得到全局最優解。用公式表述為

vi(t+1)=ωvi(t)+c1r1[Pi(t)-xi(t)]+

c2r2[Pg(t)-xi(t)]

(8)

xi(t+1)=xi(t)+vi(t+1)

(9)

式中：ω為慣性權值，調節算法收斂速度；c1和c2為學習因子，用于控制算法不陷入局部極值；r1和r2為取值在[0，1]區間的隨機數。

2.3 PSO-GSP應用計算

將粒子群優化算法應用于高階廣義屏偏移算子的計算流程如下：

(1)設置相關參數，如確定求解的維度m，粒子個數n(即求解待定常系數)，以及每個粒子初始化位置和飛行速度等；

(2)調整慣性權值，計算每個粒子的歷史最好位置和歷史全局最好位置；

(3)根據式(7)做全局最優化求解；

(4)判斷是否為全局最優值，若不是則返回(2)；

(5)將求得的待定常系數代入式(3)中進行計算。

為了直觀驗證PSO-GSP算法對成像角度的提高，現將四階PSO-GSP與常規四階GSP以及常用的兩種單程波偏移算法(SSF和FFD)在脈沖試驗中進行對比分析，如圖1～圖4所示。從單程波偏移算子的推導中可知，波場延拓過程中均會參考v0用以計算，v0即同一深度中所有點所對應實際速度v的最小值。當v0遠小于實際速度v時，即地層橫向速度變化劇烈時，速度擾動的增大就會影響波場傳遞的結果，使地震波無法以嚴格的球形傳播規律向下傳播，成像角度便會大幅度減小，以至于成像結果不準確，而不同的偏移算法對地震波保持球形傳播的能力有所不同。那么，在脈沖試驗的均勻介質中，分別選取v0=v/2和v0=v/3，檢驗四種單程波偏移算法對角度成像的優劣。

由圖1～圖4可知，對于不同的v0，SSF的成像角度均最小，四階GSP適中，FFD和PSO-GSP成像角度均較高；對比FFD、四階GSP和PSO-GSP發現，當v0=v/2時，PSO-GSP優于FFD和四階GSP，當v0=v/3時，四階GSP成像角度較小，其余兩者成像角度相當，但FFD算法可看到明顯的背景噪聲，這對偏移成像質量有影響。由此可見，四階PSO-GSP算法對橫向速度變化較為劇烈的地層成像效果更好；在地層相對平緩的地區，與FFD的成像效果相當。

圖1 v0=v/2 (a)和v0=v/3 (b)時的SSF脈沖響應

圖2 v0=v/2 (a)和v0=v/3 (b)時的FFD脈沖響應

圖3 v0=v/2 (a)和v0=v/3 (b)時的四階GSP脈沖響應

圖4 v0=v/2 (a)和v0=v/3 (b)時的四階PSO-GSP脈沖響應

圖5 當v0=v/2(a)和v0=v/3 (b)時四階PSO-GSP與各階GSP隨入射角度變化曲線圖

為進一步量化四階PSO-GSP算法在成像精度上的提高幅度，對四階PSO-GSP與各階GSP分別在不同程度的陡傾地層中成像精度隨入射角度增加的變化情況進行對比分析。如圖5a所示，以v0=v/2模擬橫向速度變化劇烈的地層情況，隨著入射角的增大，四階GSP在入射角大于30°后便逐漸偏離準確成像值；而四階PSO-GSP在入射角度小于60°時的成像精度與10階GSP相當，在小于80°時的成像精度也能保持與20階GSP相似。以v0=v/3模擬更為陡傾的地層情況(圖5b)，可見各種算法的成像誤差都有所增大，四階GSP在入射角度大于20°后誤差即大幅增加；但PSO-GSP依然保持與20階GSP相當的成像精度，在入射角大于60°后誤差才逐漸增加。因此，在相同的近似展開階數情況下，PSO-GSP的成像精度相較于GSP大幅提升，對陡傾地層的成像能力更強。

3 保幅成像條件

地震資料偏移成像的過程，不僅需對波場外推以獲得地震波在不同深度點的波場值，還需要合適的成像條件對所求取的上、下行波場進行相關成像，得到直觀的地下構造特征。而根據式(4)，在保幅的波場延拓條件下，對成像條件的要求也隨之變高。成像條件不僅需要提供準確的相位信息，更要在保幅波場延拓計算中使振幅保真，得到真實且完整的地下巖性、物性以及流體性質等儲層信息。

目前，應用最為廣泛與普遍的成像條件有兩種，分別為互相關和反褶積，其公式如下

(10)

(11)

式中σ表示阻尼因子。

這兩種成像條件均有各自的優缺點： ①互相關成像條件只有乘法運算，算法相對穩定，但其只能保證成像中相位正確，無法滿足振幅保真的需求；②反褶積成像條件加入除法運算，具有一定的振幅保真能力，但當分母項趨近于零值時，會出現成像不穩定現象；為了避免除法運算的不穩定性，需額外加入阻尼因子σ，但卻帶來了阻尼因子選取的問題，且大幅提高了成像過程中的人為因素，進而導致振幅的失真或不可控的噪聲干擾。

為了既保證在成像過程中的振幅保真，又避免成像條件的不穩定或干擾，葉月明等[20]首次提出了尋找一種窗函數平滑反褶積成像中的分母項，消除成像的奇異值影響。本文利用高斯窗函數對反褶積成像中的分母項進行平滑，其表達式如下

(12)

式中λ為平滑參數，即高斯窗的長度。該參數類似高斯分布函數的均方差，值越大表示高斯窗中的散點越多，則參與平滑運算的點數也變多，那么相對應的成像效果就會越好，其穩定性就會越高，當然計算量也會隨之增加，所以對于平滑參數λ的選取需根據實際情況而定。

利用高斯窗函數對反褶積成像條件中的分母項進行平滑，就相當于褶積濾波的過程，即將高斯窗函數當做是濾波器，反褶積的分母項當作是信號輸入，則平滑處理的表達式如下

(13)

將高斯窗函數平滑后的下行波場值代入式(11)中，即可得到基于高斯窗函數平滑的反褶積成像條件

(14)

以中點放炮、雙邊接收的簡單水平層模型為例(圖6a)，對比以上三種成像條件的結果： ①互相關(圖6b)和反褶積(圖6c)的成像結果隨著地震波入射角的增大在水平層兩端的同相軸不同程度地出現拉伸、變粗現象，而高斯窗函數反褶積(圖6d)的成像結果在不同角度入射時可保持同相軸形態的穩定；②提取水平反射層附近的最大振幅值(圖6e)發現，三種成像條件的結果有相同規律，即隨入射角增大而振幅值增大，當入射角超過可成像角度后振幅值迅速衰減。其中，高斯窗函數反褶積在可成像角度范圍內的振幅值最穩定，保持相近振幅值，且在超出可成像角度范圍后具有與互相關相近的振幅衰減速率，而反褶積的振幅衰減較慢，這是導致成像結果中有背景干擾噪聲的主要原因。對比該水平層模型三種成像條件結果可知，高斯窗函數在振幅保真方面優于另兩種傳統成像條件。

圖6 水平層模型及其成像結果

4 Marmousi模型試算

利用Marmousi模型(圖7)分別檢驗高斯窗函數對偏移成像的保幅效果及粒子群優化算法相對于傳統廣義屏算子的成像精度對比。該模型的橫向采樣點數為300，橫向采樣間隔dx=12.5m，縱向采樣點數為560，縱向采樣間隔dz=6m。用于偏移的炮集記錄共240個，120道接收，炮間距為25m，道間距為12.5m，記錄時間長度為3.5s，時間采樣間隔為4ms。

首先，應用互相關、反褶積和高斯窗函數反褶積三種不同的成像條件，對基于四階PSO-GSP波場延拓算子偏移得到波場值進行成像，得到的剖面分別如圖8～圖10所示。圖8顯示互相關成像條件的淺層成像效果較好，能量較強，但深層能量明顯減弱(白色箭頭所指)；此外，位于約2200m處的中部背斜構造的內幕細節構造成像差，難以分辨(黑色箭頭所指)。圖9是利用反褶積成像條件進行成像后的結果，選取的阻尼因子σ=0.01，圖中振幅保真較互相關成像條件有所改善，但地層傾角較大時成像干擾較嚴重(黑色方框所指)，且陡傾地層以下成像較為模糊；另外，中部背斜構造內幕的細節特征依然模糊，無太大改善。圖10是采用高斯窗函數反褶積成像條件得到的成像結果，選取的平滑參數λ=150，與前兩種方法對比，高斯窗函數反褶積成像條件的振幅保真效果較好，3000m處的古地層構造清晰可見，在地層傾角較大處也能準確描繪同相軸。在設置平滑參數后，整幅圖的噪聲干擾明顯減弱，位于約2200m處的中部背斜構造的內幕反射特征清晰可辨，同相軸的描繪與模型基本一致。上述對比分析結果充分說明高斯窗函數具有良好保幅和抗干擾能力。

圖7 Marmousi速度模型

圖8 基于保幅四階PSO-GSP算子利用互相關成像條件的成像剖面

圖9 基于保幅四階PSO-GSP算子利用阻尼反褶積成像條件的成像剖面(阻尼因子σ=0.01)

圖10 基于保幅四階PSO-GSP算子利用高斯窗函數反褶積成像條件的成像剖面(平滑參數λ=150)

其次，將四階PSO-GSP在Marmousi模型中的成像效果與四階GSP、單程波偏移中最常見的SSF和FFD對比，主要從成像精度和抗噪能力等方面探討分析。對SSF、FFD、四階GSP三種單程波偏移算子均選用與四階PSO-GSP(圖10)一致的平滑參數(λ=150)的高斯窗反褶積成像條件進行成像，所得Marmousi成像剖面如圖11～圖13所示。

圖11為SSF偏移剖面，由于SSF對強橫向速度變化的介質適應性較差，導致其對Marmousi模型中部的陡傾地層繞射波無法歸位，盡管淺部及兩側的地層構造能較好地成像，基本無法分辨斷面，位于陡傾地層下方的深部斷層由于淺層構造歸位不準確，導致深層構造形態十分模糊(白色箭頭所指)。圖12為FFD偏移剖面，相較于SSF的成像結果，FFD算子的斷面歸位清晰，地層層狀特征明顯，但對陡傾地層的成像出現假象(黑色方框所指)，且整幅剖面依然在成像平滑后有噪聲，這或許是由圖2所示的FFD在大角度時波形恢復有假象造成；另外在2200m處中部背斜地層內部的細節構造也能成像(黑色箭頭所指)，但存在干擾。圖13為四階GSP偏移剖面，該算法的成像效果與FFD較為相似，地層和構造特征基本一致，且對陡傾地層的成像均有部分歸位不完全的現象，但其噪聲干擾明顯減少。四階PSO-GSP偏移剖面(圖10))與FFD和四階GSP偏移剖面相比，可見其對陡傾地層(黑色方框所指)的成像效果最佳，斷點分明，繞射波基本完全歸位。

圖11 基于SSF真振幅偏移剖面

圖12 基于FFD真振幅偏移剖面

圖13 基于四階GSP真振幅偏移剖面

5 實際資料分析

選取一條二維地震數據進行偏移成像處理。該實際資料取自中國南海某區塊的數據，在進行疊前偏移之前，對該資料炮集記錄做了一些常規的地震資料預處理，包括定義觀測系統、線性干擾壓制、預測反褶積、疊加速度分析以及壓制海上多次波(拉東變換和SRME)等。在獲得了較為準確的層速度(圖14)之后，分別利用SSF、FFD、四階GSP和四階PSO-GSP對其進行偏移成像。

該實際地震資料速度模型尺寸為5250m×4250m，縱向采樣間隔dz=6m，橫向采樣間隔dx=25m。采集炮集記錄共計742炮，每炮360道接收，炮檢距和道間距均為25m，接收時長為5s，時間采樣間隔為2ms?；谝陨纤姆N單程波波場延拓算子偏移得到的成像剖面如圖15～圖18所示，其成像條件均為高斯窗函數的反褶積成像條件。

由于該實際資料橫向太長，圖15～圖18只截取其中第200～第1200個CDP共計1001個采樣點做對比分析。對比四幅圖發現，由于該海上資料淺層沉積較平整，四種偏移算子的成像結果差別甚微；但對于深部的背斜，四種偏移算子的成像剖面有所差異，如圖15的黑色方框所示，SSF偏移得到的地層同相軸的連續性很差，受到的干擾很嚴重，與圖16和圖17相比，其地層傾角也不一致，SSF剖面中地層為水平但不連續；而FFD和四階GSP的偏移剖面中地層有一定傾角且連續性較好，這兩種方法對此實際資料的成像能力相當。圖18為四階PSO-GSP偏移剖面，在黑色方框內分辨率比FFD和四階GSP更高，連續性更強。放大四種算法剖面的中部基巖黑色方框處(圖19)可發現，四階PSO-GSP紅色方框內的同相軸連續性更好，FFD次之，其余兩種算法同相軸連續性較差；而豎直紅線右側部分，四階PSO-GSP算法更能保持地層能量信息，優于其余三種算法。

圖14 海上實際數據層速度模型

圖16 基于FFD真振幅偏移剖面

圖17 基于四階GSP真振幅偏移剖面

圖18 基于四階PSO-GSP真振幅偏移剖面

圖19 各真振幅偏移剖面的局部放大

6 結論

本文從成像精度和保幅能力兩個方面對傳統保幅高階廣義屏疊前深度偏移算法進行改進。在成像精度方面，利用粒子群優化算法對高階廣義屏算子中垂直波數的常系數項進行全局最優化求解，進一步提高了高階廣義屏算法的成像角度；在保幅能力方面，利用高斯窗反褶積成像條件對保幅波場延拓的高階廣義屏算法的振幅保真效果進一步加強。通過脈沖試驗和Marmousi模型對四種算法的對比分析，基于粒子群優化和保幅成像條件的高階廣義屏算法相較于其他三種傳統單程波保幅偏移算法不僅有更高的成像精度、能夠更好地適應陡傾地層的成像，還能有效提高地震波振幅保真效果，為下一步解釋提供有力依據。通過對實際資料的應用效果分析，證明基于粒子群優化和保幅成像條件的高階廣義屏算法能有效提高保幅單程波偏移的成像質量，具有實際應用價值。