基于Hilbert變換的全波場分離逆時偏移成像

王一博，鄭憶康，薛清峰，常旭，Tong W.Fei，Yi Luo

1 中國科學院地質與地球物理研究所，北京 1000292 EXPEC ARC,Saudi Aramco，Dhahran 31311，Saudi Arabia

?

基于Hilbert變換的全波場分離逆時偏移成像

王一博1，鄭憶康1，薛清峰1，常旭1，Tong W.Fei2，Yi Luo2

1 中國科學院地質與地球物理研究所，北京 1000292 EXPEC ARC,Saudi Aramco，Dhahran 31311，Saudi Arabia

逆時偏移方法利用雙程波算子模擬波場的正向和反向傳播，通常采用互相關成像條件獲得偏移剖面，是一種高精度的成像方法.但是傳統的互相關成像條件會在偏移結果中產生低頻噪聲；此外，如果偏移速度中存在劇烈速度變化還可能進一步產生偏移假象.為了提高逆時偏移的成像質量，可在成像過程中先對震源波場和檢波點波場分別進行波場分離，然后選擇合適的波場成分進行互相關成像.本文基于Hilbert變換，推導了可在偏移過程中進行上下行和左右行波場分離的高效波場分離公式以及相應的成像條件，結合Sigsbee 2B合成數據，給出了不同波場成分的互相關成像結果.數值算例結果表明，采用本文提出的高效波場分離算法以及合理的波場成分互相關成像條件可以獲得高信噪比的成像結果.

逆時偏移；互相關成像條件；波場分離；Hilbert變換

1 引言

深度域偏移成像方法是油氣勘探領域的一項關鍵技術，傳統的偏移成像方法，如Kirchhoff偏移、單程波偏移等，存在一定局限性，無法滿足復雜地區，復雜地表情況下的高精度成像要求.基于射線理論的Kirchhoff積分偏移法(Wiggins，1984；Keho and Beydoun，1988)，求取的是波動方程高頻近似解，對復雜構造中的高陡傾角地層不能很好地成像.基于單程波算子的波動方程偏移法，例如相移法(Gazdag，1978)、分步傅里葉法(Stoffa et al.,1990)、傅里葉有限差分法(Ristow and Rühl,1994)等，盡管有算法穩定、背景噪聲小的優點，但在橫向速度變化特別大的時候成像精度有限.在20世紀80年代，一些學者提出了逆時偏移成像方法及實現策略(Baysal et al.，1983；McMechan,1983；Whitmore，1983)，但是受當時的計算能力限制，該方法并沒有得到廣泛應用.近年來，隨著計算能力的飛速發展，逆時偏移以其能對高陡傾角成像，并能夠處理地震數據中多種成分(反射波，繞射波，棱柱波等)的優點，在學術界和工業界都成為了研究熱點.

在逆時偏移中，首先需要求解雙程波動方程，構建震源波場和檢波點波場，然后利用合適的成像條件生成成像結果.一般我們采用的是有限差分法來求解常密度聲波方程，在時間空間域進行波場的正向和反向傳播.最常用的成像條件是經典的互相關成像條件(Claerbout，1971)，即震源波場和檢波點波場在所有時間采樣點上進行零延遲互相關后再進行疊加.這種成像方法在計算和存儲能力滿足需求的情況下比較容易實現，但其主要問題是成像結果中存在高振幅的低頻噪聲.這種低頻噪聲是由波動方程雙向傳播的特點所引起的，會嚴重干擾成像質量.為了避免在反射界面處產生強噪聲，可以在波場模擬中進行針對性處理.Baysal等(1984)提出在波場模擬中采用無反射的雙程波動方程，Loewenthal等(1987)提出對偏移速度進行充分光滑以降低反射波能量.此外，Zhang等(2009)提出采用濾波方法來消除低頻噪聲.根據低頻噪聲的產生原理，如果能在成像時分離不同方向的波場能量，則可以從本質上解決低頻噪聲問題.但常規的基于波場分離的逆時偏移成像方法，通常需要存儲所有時刻的波場值并進行傅里葉變換，需要較大的存儲量和計算量.Yoon和Marfurt(2006)利用Poynting矢量確定波場傳播方向并據此來進行波場分離成像.Liu等(2011)提出了一種基于Hilbert變換的波場分離成像條件，可以有效去除高振幅的低頻噪聲.Fei等(2015)進一步指出，在偏移速度模型存在速度劇烈變化區域時，常規逆時偏移方法會產成成像假象.為此他們利用Hilbert變換構建了一個波場分離算子，來實現下行震源波場和上行檢波點波場的互相關成像，并有效去除了成像假象，但他們只是探討了成像過程中的上下行波場分離，針對的成像問題有一定局限性.

為此，本文首先分析了常規逆時偏移過程中噪聲和假象的產生原理，然后在上下行波場分離成像條件的基礎上，構建了基于Hilbert變換的全波場分離算子，推導了上左、上右、下左以及下右四個方向的波場分離算子、相應的互相關成像條件以及具體實現策略，最后針對Sigsbee 2B模型給出了16種波場分離成像條件對應的逆時偏移成像結果.數值計算結果表明本文方法不僅可以有效地壓制常規逆時偏移成像方法的噪聲和假象，還可以根據成像目標的地質情況獲取最佳成像結果.

2 上下行波場分離及成像條件

常規逆時偏移采用的是不進行波場分離的互相關成像條件(Claerbout，1971)，其表達式為

(1)

這里，x是地下某一點的空間位置，t表示時間，Tmax是最大波場傳播時間，一般選擇為地震數據的最大記錄時間.S(x,t)是正向傳播的震源波場，R(x,t)是反向傳播的檢波點波場.這兩個波場可以通過求解雙程波動方程獲得，例如采用有限差分法求解常密度聲波方程.為了避免存儲整個正傳或反傳波場，逆時偏移可以采用以下數值步驟實現：(1) 對于每一炮，正演震源波場，記錄所有邊界點的波場值和最后時刻的波場快照；(2) 反向傳播檢波點波場，同時根據存儲的邊界點波場值和最后時刻的波場快照按時間逆推震源波場；(3) 在每個時間采樣點應用零延遲互相關成像條件，即對檢波點波場和震源波場進行點乘，然后疊加所有時間采樣點和所有炮的結果，獲得偏移剖面.

逆時偏移采用雙程波動方程模擬波場傳播過程，可以準確地生成沿各個方向傳播的波組成分.即使是在其他成像方法中難以處理的回折波或棱柱波，逆時偏移也可以對他們正確成像.但是常規互相關成像條件會在成像結果中生成強振幅的低頻噪聲，尤其當偏移速度中存在強反射界面時，噪聲更為明顯.為了提高成像精度，有必要改進常規成像條件.

假設地下介質存在水平反射層，當波場遇到該反射層時，由于雙程波動方程的特點，入射波場會分為向下透射和向上反射兩部分，可以用公式表示為

(2)

(3)

(4)

Liu等(2011)指出正是最后兩項，SdRd和SuRu，生成了逆時偏移中強振幅的低頻噪聲.這兩部分表示的是波場中沿相同方向傳播的能量，不應該包括在互相關成像條件中.為此，他們提出應當在公式(4)中去掉最后兩項.Fei等(2015)進一步指出，在偏移速度中存在速度變化特別劇烈的區域時，公式(4)中的第二項SuRd將產生假象并嚴重干擾成像質量，應該在互相關成像條件中去除，則公式(4)所示的成像條件可進一步表示為

(5)

對于公式(5)所示的波場分離成像條件，其實現難點在于如何快速有效地進行波場分離.以二維地震數據的逆時偏移為例，傳統的波場分離方法是將震源與檢波點波場全部存儲，然后變換到頻率波數域，采用公式(6)至公式(9)實現波場分離，獲得上行波和下行波：

(6)

(7)

(8)

(9)

這里的S(kz,ω)對應震源波場的二維傅里葉變換，R(kz,ω)對應檢波點波場的二維傅里葉變換.公式(6)至公式(9)所示的分離方法需要較大的數據存儲量，在實際生產中存在一定問題.為此Fei等(2015)提出了一種利用Hilbert變換的高效波場分離算法.

對于一個給定的函數f(t)，Hilbert變換具有如下性質：

(10)

這里的F表示傅里葉變換算子，H表示Hilbert變換算子，下標t表示變換作用的變量，sgn表示符號函數，定義為

(11)

可以構建波場分離算子U：

(12)

該算子的傅里葉變換具有如下性質：

(13)

(14)

-Hz(s)Ht(r)}dt,

(15)

這里的s代表震源波場，r代表檢波點波場，r滿足r=g*r0，其中g是基于波動方程計算獲得的格林函數，r0表示實際地震記錄，可證明Ht(r)=g*Ht(r0)，即Ht(r)表示的是一個地震記錄經過Hilbert變換之后正演模擬得到的波場.

3 全波場分離及成像條件

在某些復雜情況下，僅使用上下行波的波場分離成像條件是不夠的，需要進一步獲取精細成像結果，例如在地下存在高陡傾角構造時，就需要將上下行波場進一步分離為左行與右行波場.

震源與檢波點波場按照左上、左下、右上與右下四個方向分解后可以表示為

S(x,t)=Slu(x,t)+Sld(x,t)+Sru(x,t)+Srd(x,t),

(16)

R(x,t)=Rlu(x,t)+Rld(x,t)+Rru(x,t)+Rrd(x,t),

(17)這里的下標l和r分別表示左行波和右行波.把公式(16)和公式(17)代入互相關成像條件(1)中，將獲得16個成像結果.可以根據地層傾角特點對16個成像結果進行選擇性組合，以獲取相對合理的成像結果，例如對位于震源和檢波點下方左側的垂直反射層成像時，相對合理的成像結果可能是

(18)

以下推導全波場分離公式，首先定義一個新的算子V，其表達式為

(19)

將該算子作用于空間變量x,可得

(20)

=s-HzHt(s)-HxHt(s)-HxHz(s),

(21)

=s+HzHt(s)-HxHt(s)+HxHz(s),

(22)

=s-HzHt(s)+HxHt(s)+HxHz(s),

(23)

=s+HzHt(s)+HxHt(s)-HxHz(s).

(24)

實際計算時，只需在常規逆時偏移基礎上額外計算一個Hilbert變換后的地震記錄的正演波場，對其沿不同坐標方向進行Hilbert變換并加以組合，就可以獲得四個方向的波場分離結果.整個計算過程不需要存儲波場快照和進行傅里葉變換，具有很高的計算效率.

根據上述波場分離公式，可進一步推導相應的波場分離互相關成像條件，例如左下行震源波場和右上行檢波點波場互相關成像條件為

(25)

值得注意的是由于算子Ut作用的震源波場或者檢波點波場在頻率域內對于負的ω均為0，我們只需在時間變量上對震源波場或者檢波點波場作用一次算子Ut，則公式(25)可以簡化為

(26)

將波場分離算子的定義公式(12)和(19)代入公式(26)，可得左下行震源波場與右上行檢波點波場互相關成像條件的最終表達式：

-Hz(s)Ht(r)-Hz(s)Hz(r)-Hz(s)Hx(r)+Hz(s)HxHzHt(r)

-Hx(s)Ht(r)-Hx(s)Hz(r)-Hx(s)Hx(r)+Hx(s)HxHzHt(r)

-HxHz(s)r+HxHz(s)HzHt(r)+HxHz(s)HzHt(r)+HxHz(s)HzHz(r)]dt.

(27)

在對波場進行左右行和上下行分離后，我們一共可以得到16種互相關成像結果，每一種成像公式的表達形式均與公式(25)類似.我們認為其中有四種成像結果是Sigsbee 2B模型相對合理的成像結果，第一種即為公式(27)所示，其余三種成像公式經推導具有如下表達式.

第二種是左下行震源波場和左上行檢波點波場的互相關：

-Hz(s)Ht(r)-Hz(s)Hz(r)+Hz(s)Hx(r)-Hz(s)HxHzHt(r)

-Hx(s)Ht(r)-Hx(s)Hz(r)+Hx(s)Hx(r)-Hx(s)HxHzHt(r)

-HxHz(s)r+HxHz(s)HzHt(r)-HxHz(s)HzHt(r)-HxHz(s)HzHz(r)]dt.

(28)

第三種是右下行震源波場和右上行檢波點波場的互相關：

-Hz(s)Ht(r)-Hz(s)Hz(r)-Hz(s)Hx(r)+Hz(s)HxHzHt(r)

+Hx(s)Ht(r)+Hx(s)Hz(r)+Hx(s)Hx(r)-Hx(s)HxHzHt(r)

+HxHz(s)r+HxHz(s)HzHt(r)-HxHz(s)HzHt(r)-HxHz(s)HzHz(r)]dt.

(29)

第四種是右下行震源波場和左上行檢波點波場互相關：

-Hz(s)Ht(r)-Hz(s)Hz(r)+Hz(s)Hx(r)-Hz(s)HxHzHt(r)

+Hx(s)Ht(r)+Hx(s)Hz(r)-Hx(s)Hx(r)+Hx(s)HxHzHt(r)

+HxHz(s)r-HxHz(s)HzHt(r)+HxHz(s)HzHt(r)+HxHz(s)HzHz(r)]dt.

(30)

當炮點和檢波點都在地表，地下只有水平反射層時，僅有下行震源波場和上行檢波點波場的互相關是正確的成像結果.但如果地下反射層的傾角接近90°時，就應該采用左行震源波場與右行檢波點波場的互相關成像結果，或者采用右行震源波場與左行檢波點波場的互相關成像結果.我們需要針對待處理工區地質結構的特點，對各種成像結果進行選擇性組合，以獲取最佳成像效果.

4 數值算例

我們采用一個簡單模型來驗證前文推導的波場分離公式的有效性.圖1是一個雙層速度模型，該模型的水平方向有300個網格點，垂直方向有150個網格點，網格點間距為10 m.震源位置在地表水平坐標1500 m處.圖2是0.6 s時刻的震源波場快照.圖3a至圖3d是震源波場的四個方向分離結果，如圖所示，波場的左下、左上、右下與右上四個不同方向的波場成分均得到了有效分離.

我們采用Sigsbee 2B數據來驗證波場分離成像條件的有效性.如圖4a所示，該模型的水平方向有3201個網格點，垂直方向有1201個網格點，網格點間距為7.62 m.該數據有496炮，炮點間距是45.72 m，單炮最大道數為348道，炮點間距為45.72 m，檢波點間距為22.86 m，炮點和檢波點深度均為7.62 m，炮記錄的時間長度為12 s，時間采樣間隔為8 ms.我們采用自由表面多次波去除之后的數據進行逆時偏移.

圖1 雙層速度模型Fig.1 A two-layer velocity model

圖2 0.6 s的震源波場快照Fig.2 The source wavefield snapshot at 0.6 s

圖3 四個方向的波場分離結果：(a) 左下；(b) 左上；(c) 右下；(d) 右上Fig.3 Four directional wavefield after decomposition.(a) Left-downgoing;(b) Left-upgoing;(c) Right-downgoing;(d) Right-upgoing

圖4 (a) Sigsbee 2B 速度模型；(b) Sigsbee 2B偏移速度模型Fig.4 (a) Sigsbee 2B stratigraphic velocity model;(b) Sigsbee 2B migration velocity model

圖5 (a) Sigsbee 2B數據的常規逆時偏移剖面；(b) 對圖(a)所示的常規逆時偏移剖面進行Laplace濾波的結果Fig.5 (a) The conventional RTM result of Sigsbee 2B dataset;(b) Laplace filtering result of the conventional RTM image as shown in (a)

圖8 左下類型的震源波場與分離后的兩個檢波點波場進行互相關成像的結果.(a) SldRld；(b) SldRrd

圖6 Sigsbee 2B模型對應的四種相對合理的成像結果(a) SldRru；(b) SldRlu；(c) SrdRru；(d) SrdRlu.Fig.6 Four relatively correct RTM images related with Sigsbee 2B model

圖7 左上類型的震源波場與分離后的四個檢波點波場進行互相關成像的結果(a) SluRlu;(b) SluRld;(c) SluRru;(d) SluRrd.Fig.7 The RTM images constructed by the cross-correlation of the left-upgoing type source wavefield with decomposed four receiver wavefields

圖9 右上類型的震源波場與分離后的四個檢波點波場進行互相關成像的結果(a) SruRlu；(b) SruRld；(c) SruRru；(d) SruRrd.Fig.9 The RTM images constructed by the cross-correlation of the right-upgoing type source wavefield with decomposed four receiver wavefields

圖10 右下類型的震源波場與分離后的兩個檢波點波場進行互相關成像的結果(a) SrdRld；(b) SrdRrd.Fig.10 The RTM images constructed by the cross-correlation of the right-downgoing type source wavefield with decomposed two receiver wavefields

圖4b是Sigsbee 2B數據對應的偏移速度模型，圖5a是采用常規互相關成像條件的逆時偏移結果，圖5b是圖5a經過Laplace濾波之后的結果.圖6a至圖6d分別是采用公式(27)至公式(30)所示的四種波場分離互相關成像條件獲得的逆時偏移成像結果.圖6a至圖6d是Sigsbee 2B模型對應的四個相對合理的成像結果，但圖6b與圖6c的鹽丘上部均存在一些噪聲，這是由于鹽丘上邊界為傾斜邊界，而圖6b與圖6c對應的成像公式中存在水平方向同向傳播的震源與檢波點波場，這兩個同向傳播波場的互相關導致了鹽丘上部出現低頻噪聲.圖7a至圖7d分別是左上類型的震源波場與分離后的四種檢波點波場采用互相關成像條件的逆時偏移結果.圖8a與圖8b分別是左下類型的震源波場與左下、右下兩種檢波點波場采用互相關成像條件的逆時偏移結果.圖9a至圖9d分別是右上類型的震源波場與分離后的四種檢波點波場采用互相關成像條件的逆時偏移結果.圖10a與圖10b分別是右下類型的震源波場與左下、右下兩種檢波點波場采用互相關成像條件的逆時偏移結果.圖7至圖10所示的12種成像結果大都是偏移噪聲和假象，主要是由同方向傳播的震源與檢波點波場互相關產生的，如SldRld,SluRlu,SrdRrd,SruRru等.

5 結論

逆時偏移是地震勘探領域一種重要的成像方法，它成像精度高并且適用于復雜地質構造成像，但常規逆時偏移結果存在低頻噪聲，而且在偏移速度梯度較大時還可能存在偏移假象.本文在上下行波場分離互相關成像條件的基礎上，推導了基于Hilbert變換的全波場分離公式以及相應的互相關成像條件，并針對Sigsbee 2B模型給出了16種采用全波場分離互相關成像條件的偏移結果，可以有效地分離常規逆時偏移方法產生的噪聲和假象.本文建立的全波場分離逆時偏移成像方法計算效率高，易于實現，在實際應用中可以利用本文方法根據目標地層的特點高效地獲取最優成像結果，這對于逆時偏移的實際應用有一定意義.

Baysal E,Kosloff D D,Sherwood J W C.1983.Reverse time migration.Geophysics,48(11):1514-1524.

Baysal E,Kosloff D D,Sherwood J W C.1984.A two-way nonreflecting wave equation.Geophysics,49(2):132-141.

Claerbout J F.1971.Toward a unified theory of reflector mapping.Geophysics,36(3):467-481.

Fei T W,Luo Y,Yang J R,et al.2015.Removing false images in reverse time migration:The concept of de-primary.Geophysics,80(6):S237-S244.

Gazdag J.1978.Wave equation migration with the phase-shift method.Geophysics,43(7):1342-1351.

Keho T H,Beydoun W B.1988.Paraxial ray Kirchhoff migration.Geophysics,53(12):1540-1546.

Liu F Q,Zhang G Q,Morton S A,et al.2011.An effective imaging condition for reverse-time migration using wavefield decomposition.Geophysics,76(1):S29-S39.

Loewenthal D,Stoffa P L,Faria E L.1987.Suppressing the unwanted reflections of the full wave equation.Geophysics,52(7):1007-1012.

McMechan G A.1983.Migration by extrapolation of time-dependent boundary values.Geophysical Prospecting,31(3):413-420.Ristow D,Rühl T.1994.Fourier finite-difference migration.Geophysics,59(12):1882-1893.

Stoffa P L,Fokkema J T,de Luna Freire R M,et al.1990.Split-step Fourier migration.Geophysics,55(4):410-421.

Whitmore N D.1983.Iterative depth migration by backward time propagation.∥53rd Annual International Meeting.SEG Expanded Abstracts,382-385.

Wiggins J W.1984.Kirchhoff integral extrapolation and migration of nonplanar data.Geophysics,49(8):1239-1248.

Yoon K,Marfurt K J.2006.Reverse-time migration using the Poynting vector.Exploration Geophysics,37(1):102-107,doi:10.1071/EG06102.

Zhang Y,Sun J.2009.Practical issues in reverse time migration:True amplitude gathers,noise removal and harmonic source encoding.First Break,26:29-35.

(本文編輯 胡素芳)

Reverse time migration with Hilbert transform based full wavefield decomposition

WANG Yi-Bo1,ZHENG Yi-Kang1,XUE Qing-Feng1,CHANG Xu1,FEI W.Tong2,LUO Yi2

1 Institute of Geology and Geophysics,Chinese Academy of Sciences,Beijing 100029,China2 EXPEC Advanced Research Center,Saudi Aramco,Dhahran 31311,Saudi Arabia

Reverse time migration (RTM) is an effective technique for complex subsurface imaging.It uses two-way wave-equation in wavefield forward and backward simulation,and employs cross-correlation imaging condition.It can offer subsurface images with high accuracy and high resolution.However,the conventional imaging condition generates low-frequency noises and may form false images around the strong velocity gradients or velocity interfaces in the migration velocity model.The wavefield components which propagate along different directions should be decomposed and selectively cross-correlated to improve the final image quality.In this manuscript,we use an efficient operator which based on Hilbert transform to decompose down-up and left-right going wavefields.The decomposed wavefields are then selectively cross-correlated according to the acquisition geometry and target location.In the numerical example section,the proposed RTM with wavefield decomposition is applied to Sigsbee 2B synthetic dataset.The results show that the effective wavefield decomposition operator combined with appropriate imaging condition can generate subsurface images with high signal-to-noise ratio.

Reverse time migration;Cross-correlation imaging condition;Hilbert transform;Wavefield decomposition

王一博，鄭憶康，薛清峰等.2016.基于Hilbert變換的全波場分離逆時偏移成像.地球物理學報，59(11):4200-4211,

10.6038/cjg20161122.

Wang Y B,Zheng Y K,Xue Q F.2016.Reverse time migration with Hilbert transform based full wavefield decomposition.Chinese J.Geophys.(in Chinese),59(11):4200-4211,doi:10.6038/cjg20161122.

國家自然科學基金項目(41422403)資助.

王一博，研究員，博士生導師，主要從事勘探地球物理和儲層地質力學方面的研究工作.E-mail:wangyibo@mail.iggcas.ac.cn

10.6038/cjg20161122

P631

2016-05-17，2016-07-02收修定稿