康智清 (中石化河南油田地球物理勘探公司物探研究所,河南 南陽473132)
劉恩 (中石化中原油田地質(zhì)錄井公司,河南 濮陽457000)
司杰戈 (中石化勝利油田物探公司,山東 東營257000)
陳康 (中國石油大學(xué) (華東)地球科學(xué)與技術(shù)學(xué)院,山東 青島266555)
逆時偏移是地震偏移方法的重要發(fā)展。與傳統(tǒng)偏移方法不同,逆時偏移是在時間軸上實現(xiàn)外推,可以看作是沿時間反方向的正演模擬過程。傳統(tǒng)的沿深度方向的偏移方法基于單程波方程,而逆時偏移則基于全波方程,允許地震波在全方位傳播,因而不存在傾角限制。隨著計算機技術(shù)的發(fā)展和對復(fù)雜構(gòu)造成像的更高要求,逆時偏移技術(shù)的研究也漸漸深入,由以前的二維發(fā)展到三維、聲波發(fā)展到彈性波、各向同性發(fā)展到各向異性。但是偏移假象仍然是影響逆時偏移剖面效果的重要方面[1~5]。
整體看來,逆時偏移需要解決的問題主要包括幾方面。①正演方面的問題:包括波動方程的改進;正演模擬的精度;邊界吸收處理;數(shù)值頻散;穩(wěn)定性問題以及彈性波的波場分離問題。②逆時延拓過程中的問題:逆時延拓中產(chǎn)生的次生干擾;反射和透射損失問題。③成像條件的問題:如何獲得正確的保幅的成像條件;如何通過成像條件消除各種偏移假象。④偏移假象問題:假象產(chǎn)生的原因和機制;直達波、邊界反射、回轉(zhuǎn)波等的影響造成的偏移假象。⑤算法的效率問題:通過改進算法提高運算效率;通過GPU和并行運算的方法提高運算效率。筆者主要研究第④方面的問題,關(guān)于逆時偏移去噪的研究,主要包括修改波動方程、修改成像條件、成像后進行濾波處理等。具體地說,Valenciano等[6]根據(jù)反演理論提出了反褶積成像條件。Mulder等[5]運用一個空間域的低通濾波器來消除噪聲。Yoon等[3]提出了在零延遲互相關(guān)條件中加入Poynting矢量來消除成像噪聲。Liu等[4]把全波場分解成單程波分量,并運用成像條件去結(jié)合這些波場分量達到消除成像噪聲的目的。還有方向性衰減、速度平滑、波場分離成像、無反射波動方程逆時偏移成像等方法[7~12]。
筆者從偏移假象產(chǎn)生的機制出發(fā),重點分析了逆時偏移假象產(chǎn)生的各種原因,針對這些原因提出和總結(jié)了解決的方法。特別針對成像條件引入的低頻噪聲,筆者采用了拉普拉斯濾波算法,通過模型測試結(jié)果可知該算法效果較好。
逆時偏移主要包括波動方程的正演模擬、逆時外推和成像條件的確定3個步驟[13],在進行正演模擬和逆時外推過程中分別保存其波場,然后運用成像條件進行求和,得到局部成像數(shù)據(jù)體。最后將所有炮集的逆時偏移結(jié)果進行疊加,得到最終的疊前深度偏移成像結(jié)果。三維逆時偏移的定解問題可以描述為:
利用微分和差分關(guān)系,對式(1)進行差分離散,得到用于正演模擬和逆時深度偏移的高階差分方程的初始方程:
式中:u為不同時刻的波長值,m;x、y、z為三維坐標(biāo);v為速度,m/s;t為時間,s;Δt為時間采樣間隔,s;O為高階項。式 (2)代表正演過程,式 (3)代表逆時外推過程。將得到的正演波場和逆時外推波場進行互相關(guān)就可以得到成像數(shù)據(jù)體,最后進行疊加就可以得到成像剖面。
在逆時偏移第一個階段,炮點波場在正向延拓過程中會產(chǎn)生直達波,如果將含有直達波的波場進行互相關(guān)就會在成像剖面頂端產(chǎn)生很嚴(yán)重的偏移假象。因此,必須在正演記錄中將直達波切除,消除直達波的影響。圖1(a)為速度模型,模型大小為200(道)×200(采樣點),網(wǎng)格大小為10m×10m,在第100道激發(fā),各道接收,可以得到圖1(b),此時為沒有去除直達波的正演記錄。用圖1(b)的正演記錄進行逆時外推,最終進行互相關(guān)成像就可以得到圖1(c),可以看到,在成像結(jié)果頂端產(chǎn)生偏移假象。圖1(d)為去除了直達波的正演記錄。圖1(e)為對應(yīng)的單炮逆時偏移結(jié)果,可以看出,去除直達波以后,頂端的偏移假象得到了消除。
圖1 直達波對逆時偏移成像剖面的影響
在正演的過程中,邊界的處理是很重要的,如果邊界問題處理不好,將會引入邊界反射,邊界反射會將一些干擾波引入到正演記錄中,最終直接導(dǎo)致偏移剖面中出現(xiàn)偏移假象問題。圖2(a)為層狀速度模型,圖2(b)代表含有邊界反射的正演記錄。圖2(c)代表含邊界反射的單炮偏移剖面,可以明顯地看到,在偏移剖面中,存在另外一些同相軸,即是邊界反射帶來的偏移假象。筆者采用了PML(perfect matched layer) 吸 收 邊 界 條件消除邊界反射,通過邊界反射的消除,圖2(d)中的正演記錄不含邊界反射,圖2(e)為消除邊界反射后對應(yīng)的單炮偏移剖面。對比圖2(e)和(c)可知,通過邊界條件的使用,偏移剖面中由于邊界反射引起的偏移假象得到了很好的壓制。
圖2 直達波對成像剖面的影響
目前,常用的成像條件有零時刻成像條件、互相關(guān)成像條件、波阻抗成像條件[14]。在逆時偏移成像中,通常使用的是互相關(guān)成像條件[15]。由于一些不正確的互相關(guān),導(dǎo)致了成像剖面中存在著偏移噪聲,這些噪聲主要分布在成像剖面淺部,遮蓋地下真實形態(tài),從而降低了成像質(zhì)量。下面具體以層狀模型為例分析噪聲的形成機制。
圖3(a)為t=0.2s時刻的炮點波場,波是向下傳播的。圖3(b)為對應(yīng)于圖3(a)的逆時波場,在對應(yīng)的點波傳播方向是向上的,圖3(a)和 (b)的互相關(guān)就形成了圖3(c)中的一道弧線,即為成像噪聲。圖3(d)為t=0.4s時刻的炮點波場,波場存在向下傳播的透射波和向上傳播的反射波;圖3(e)為同時刻的逆時外推波場,波場中也存在逆時入射波和逆時反射波,圖3(d)和 (e)的互相關(guān)就構(gòu)成了圖3(f)的A、B兩點和一條弧線,其中A、B兩點是正確的像點,而弧形是互相關(guān)引入的噪聲。
對于成像過程中的噪聲,目前已經(jīng)由很多的學(xué)者提出很多去除的方法,主要有方向性衰減、速度平滑、波場分離成像、無反射波動方程逆時偏移成像等。但是方向性衰減、波場分離成像、波因廷矢量成像條件等在實現(xiàn)上有很大的難度。無反射波動方程效果不是太好,特別是對于大角度入射時反射還是比較明顯。速度平滑是一種比較簡單的方法,效果也不錯,可是對精確的速度模型來說不是一個好的方法,因為速度模型的平滑引入了速度誤差。因此,筆者采用拉普拉斯濾波算法消除成像過程中的噪聲,它在去噪的過程中既起到了高通濾波的作用,也有速度平滑的效果。
對于常規(guī)的二維拉普拉斯濾波算子:
式中:為拉普拉斯算子;I為波場值。從式(4)中可以看出,拉普拉斯算子可以表示為一個二階微分的形式。
圖3 成像條件引入噪聲的形成機制
圖4 為不同階數(shù)拉普拉斯算子對低頻的壓制效果,可以看到,微分有提升高頻的作用,同時可以壓制低頻,且隨著階數(shù)越高,微分作用對低頻的削弱越明顯,對高頻成分也是非線性的提升,2階微分對高頻的提升和對低頻的壓制作用明顯強于1階微分。通過圖4可以得出下面的結(jié)論:拉普拉斯算子是基于2階微分的組合算子,具有很好的壓制低頻和提升高頻的作用,把一個信號或者圖像通過拉普拉斯算子,其過程相當(dāng)于一個高通濾波器,相比高通濾波器,拉普拉斯算子濾波還能起到平滑的效果。
圖5(a)是一個合成地震記錄的深度域表示,由一個傾斜軸和一個水平軸組成,Kz表示Z方向的波數(shù)域,Kx表示X方向的波數(shù)域。圖5(b)是將圖5(a)通過拉普拉斯算子后的結(jié)果,可以發(fā)現(xiàn),通過拉普拉斯算子,低頻成分得到了很好的壓制,同時提升和保留了高頻成分。
圖4 不同階數(shù)拉普拉斯算子對低頻的壓制效果
圖5 濾波前、后波數(shù)域分析
圖6 、7是對拉普拉斯濾波的效果測試。圖6(a)和圖7(a)分別是速度模型和sigbee模型,圖6(b)是30炮逆時偏移剖面,可以看到成像噪聲分布在剖面各個地方,嚴(yán)重影響了剖面的質(zhì)量;同理可以看到圖7(b)剖面中的噪聲也是非常嚴(yán)重。圖6(c)和圖7(c)是將含噪剖面進行拉普拉斯算子濾波的結(jié)果,可以看出,通過濾波,成像剖面的質(zhì)量都得到了很好的改善,成像噪聲得到了很好的壓制。
圖6 速度模型拉普拉斯算子濾波效果
圖7 sigbee模型拉普拉斯算子濾波效果
筆者分成像前和成像過程兩個階段分析了逆時偏移剖面中的偏移假象的原因,成像前主要有直達波和邊界反射的影響,成像過程主要是成像條件的影響。針對這些問題,筆者提出了有效的解決方法,在成像前分別用切除直達波和采用PML邊界條件的方法達到目的;對于成像條件引入的噪聲,筆者采用了拉普拉斯算法進行壓制,通過模型測試結(jié)果可以看出,拉普拉斯算法能夠有效地對成像條件引入的噪聲進行壓制。
[1]Baysal E,Kosloff D D,Sherwood J W C.Reverse-time migration [J].Geophysics,1984,48 (4):1514~1524.
[2]Chang W F,McMechan G A.Elastic reverse-time migration [J].Geophysics,1987,52 (4):1365~1378.
[3]Yoon K,Marfurt K J.Reverse-time migration using the Poynting vector.Exploration Geophysics,2006,37 (1):102~107.
[4]Liu Faqi,Zhang,Guanquan,Morton S A,et al.An effective imaging condition for reverse-time migration using wavefield decomposition [J].Geophysics,2011,76 (1):S29~S39.
[5]Mulder W A,Plessix R E.A comparison between one way and two way wave equation migtation [J].Geophysics,2004,69 (6):1491~1504.
[6]Valenciano A A,Biondi B.Deconvolution imaging condition for reverse-time migration [J].Stanford Exploration Project,2002,Report 112:83~96.
[7]Liu Faqi,Zhang G,Morton S A,et al.Reverse-time migration using one-way wavefield imaging condition [A].2007SEG Annual Meeting [C].San Antonio,2007-09-23~28.
[8]Guitton A,Kaelin B.Least-squares attenuation of reverse-time migration artifact [J].Geophysics,2004,69 (6):1491~1504.
[9]Rickett J,Sava P.Offset and angle-domain common image-point gathersfor shot-profile migration [J].Geophysics,2002,67 (7):883~889.
[10]Fletcher R P,F(xiàn)owler P J.Suppressing unwanted internal reflections in prestack reverse-time migration [J].Geophysics,2006,71(6):79~82.
[11]紅偉,劉洪,鄒振 .地震疊前逆時偏移中的去噪與存儲 [J].地球物理學(xué)報,2010,53(9):2171~2180.
[12]張軍華,呂寧,田連玉,等 .地震資料去噪方法技術(shù)綜述評述 [J].地球物理學(xué)進展,2006,21(2):546~553.
[13]楊仁虎,常旭,劉伊克 .疊前逆時偏移影響因素分析 [J].地球物理學(xué)報,2010,53(8):1902~1913.
[14]陳可洋 .基于高階有限差分的波動方程疊前逆時偏移方法 [J].石油物探,2010,49(2):153~157.
[15]丁亮,劉洋 .逆時偏移成像技術(shù)研究進展 [J].地球物理學(xué)進展,2011,26(3):1085~1100.