李 勤,沈鴻雁,王 鑫,趙 靜,李 萌
(1.西安石油大學(xué) 地球科學(xué)與工程學(xué)院,陜西 西安 710065;2.陜西省油氣成藏地質(zhì)學(xué)重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710065)
煤田地球物理勘探的任務(wù)之一是準(zhǔn)確識(shí)別和追蹤斷層、陷落柱、富水層等異質(zhì)體(區(qū)),以避免(或減輕)開采過程中的地質(zhì)災(zāi)害[1-9]。地震繞射波由于具有較高的分辨力和與小尺度異質(zhì)體的準(zhǔn)確對(duì)應(yīng)關(guān)系,近年來備受勘探人員重視。其實(shí),早在20世紀(jì)50年代,KREY[10]等就已經(jīng)發(fā)現(xiàn)來自斷塊、尖滅等處的繞射波,之后TROREY[11],BERRYHILL[12]和Klem-Musatov[13]等對(duì)其進(jìn)行了研究并建立了基本的繞射波理論體系。隨著研究的不斷深入,人們發(fā)現(xiàn)除了能有效檢測(cè)異質(zhì)體的存在,繞射波還可用于描述異質(zhì)體的方位及邊界等特征[14],充分利用繞射波對(duì)地下異質(zhì)體(區(qū))的精細(xì)勘探具有十分重要的意義[15-18]。
繞射波利用的關(guān)鍵是對(duì)其進(jìn)行準(zhǔn)確成像,但獨(dú)立的時(shí)差規(guī)律使得它在與反射波同時(shí)成像時(shí)不僅得不到正確收斂,而且還會(huì)因能量較弱而被掩蓋。此外,繞射成像對(duì)速度擾動(dòng)較為“敏感”[19],很小的速度變化都可能使其離開最佳的收斂狀態(tài),因此偏移時(shí)需要提供精確的速度模型,這也導(dǎo)致了繞射成像實(shí)現(xiàn)起來比較困難。雖然近年發(fā)展出許多專門針對(duì)繞射波的成像方法和技術(shù),也取得了一定的成像效果[20-25],但大多還是對(duì)速度模型有著較強(qiáng)的依賴[26]。
路徑積分法是一種較新的地震成像技術(shù),因無需預(yù)先構(gòu)建成像速度模型就可實(shí)現(xiàn)成像,所以引起了研究者的極大關(guān)注[27-30]。目前,該方法在直接成像中的應(yīng)用基于稱為“速度延拓”的偏移算法。CLAERBOUT[31]提出了速度延拓偏微分方程,將其應(yīng)用于速度分析;FOMEL[32-33]指出在各向同性情況下,零偏移距速度延拓是一個(gè)連續(xù)和解析的過程,并給出了方程的f-k域解;BURNETT等[34]認(rèn)為路徑積分和速度延拓的結(jié)合拓展出一種不依賴速度模型的成像方法;LANDA等[35]由共成像點(diǎn)道集的平面度計(jì)算了權(quán)值,說明加權(quán)因子的引入可增強(qiáng)路徑積分成像的協(xié)調(diào)性;MERZLIKIN等[36]進(jìn)行了路徑積分法繞射波成像,并且推導(dǎo)了偏移濾波器的解析形式。上述研究驗(yàn)證了路徑積分成像的可行性,探索了其在地震勘探中的應(yīng)用潛力。
筆者在已有理論的基礎(chǔ)上,研究針對(duì)疊后繞射波的f-k域路徑積分成像方法,以進(jìn)一步改善繞射波的成像質(zhì)量,提高破碎帶、巖溶、采空區(qū)等小尺度異質(zhì)體(區(qū))的探測(cè)精度。
速度延拓(Velocity Continuation,VC)是一種以速度為傳播變量的波形過程[33],描述了時(shí)間偏移同相軸在不同速度下的形態(tài)變化。零偏移距各向同性時(shí),該過程可表示為一個(gè)線性雙曲形式的偏微分方程[31,33,36]:
(1)
式中,P(x,t,v)為地震波場(chǎng),x,t,v分別為空間坐標(biāo)、單程走時(shí)和半偏移速度。
當(dāng)v=C時(shí)(C為非負(fù)常數(shù)),求解式(1)的過程稱為等速延拓,解得的P(x,t,C)即該速度下的恒速時(shí)間偏移剖面。一般可引入變量τ來替換單程走時(shí)t(τ=t2),通過對(duì)時(shí)間軸的“拉伸”(重采樣)處理將方程簡化為傅里葉域中的常微分方程并得到其解析解:
(2)
基于速度延拓的連續(xù)性以及繞射同相軸頂點(diǎn)在該延拓過程中的“平穩(wěn)性”,可利用連續(xù)疊加等速延拓剖面的方式來實(shí)現(xiàn)繞射波成像[37-39],繼而得到基于路徑積分的速度延拓偏移濾波器:
(3)
其中,F(xiàn)PI(Ω,k,va,vb)為f-k域路徑積分偏移濾波器,其性能取決于Ω,k以及速度積分限va和vb。通過式(3)得到的成像結(jié)果可認(rèn)為是由給定區(qū)間內(nèi)的所有等速延拓剖面等權(quán)重疊加而成。
為了進(jìn)一步增加成像的協(xié)調(diào)性,提高成像質(zhì)量,可在上述濾波器中加入高斯因子:
(4)
濾波器的解析表達(dá)式(式(4)第2個(gè)等號(hào))基于復(fù)指數(shù)積分函數(shù)與虛誤差函數(shù)成正比的性質(zhì)得出,通過求取虛誤差函數(shù)端點(diǎn)值的方式,就可快速求出相應(yīng)的偏移濾波器。濾波結(jié)果經(jīng)二維逆傅里葉變換后,可進(jìn)一步輸出為成像剖面。與疊后t-x域的偏移不同,f-k域路徑積分法主要利用濾波運(yùn)算進(jìn)行幅值和相位控制,通過增強(qiáng)繞射波場(chǎng)中的水平(近水平)分量、弱化具有一定傾角的繞射弧側(cè)翼來實(shí)現(xiàn)繞射波的偏移成像。
f-k域路徑積分法疊后繞射波偏移成像流程如圖1所示,具體實(shí)施步驟:
圖1 路徑積分法成像流程Fig.1 Flowchart of path integral method
(1)輸入原始地震記錄(水平疊加剖面),對(duì)其進(jìn)行預(yù)處理(能量補(bǔ)償、干擾波壓制等),然后選擇合適的波場(chǎng)分離方法提取繞射波場(chǎng)Dx×y(x行y列)??蛇x用SVD(Singular Value Decomposition)法提取繞射波[40-44],基于水平疊加剖面中反射波和繞射波橫向相干性和連續(xù)性方面的明顯差異,通過SVD帶通濾波實(shí)現(xiàn)2者的分離。
(3)按照拉伸填充后的時(shí)空域剖面的網(wǎng)格參數(shù)生成坐標(biāo)矩陣Ωm×1和k1×n,再選取積分限va和vb,由式(4)構(gòu)造路徑積分偏移濾波器Fm×n(復(fù)矩陣)。
(1)模型構(gòu)建。建立了一個(gè)單繞射點(diǎn)地質(zhì)模型來檢驗(yàn)方法的正確性(圖2(a))。如圖2(a)所示,速度為1 500 m/s的均勻地層中,存在1個(gè)半徑為15 m的球形異質(zhì)體,球心距地表375 m,其速度為4 500 m/s。采用有限差分解波動(dòng)方程的方法進(jìn)行正演模擬,自激自收觀測(cè)系統(tǒng),震源子波為Ricker子波,主頻為40 Hz,道距為2 m,道數(shù)為1 000道,采樣率為1 ms,采樣點(diǎn)數(shù)為1 000點(diǎn)。
正演獲得的地震剖面如圖2(b)所示,該剖面可等效為水平疊加剖面。從該地震剖面可看出,異質(zhì)體引起的地震響應(yīng)為一條具有雙曲線特征的繞射波同相軸,頂點(diǎn)在500 ms處,跨度約2 km,能量隨著遠(yuǎn)離繞射源呈現(xiàn)出衰減特性。為了貼近實(shí)際情況,在地震剖面中加入了隨機(jī)噪聲(圖2(c)),其信噪比約為10∶1。
圖2 單點(diǎn)繞射模型及地震剖面Fig.2 Single-point diffraction model and seismic section
(2)等速延拓成像。由于路徑積分速度延拓成像基于等速延拓原理,因此需要先檢驗(yàn)等速延拓的有效性。根據(jù)偏移流程,首先對(duì)原始數(shù)據(jù)(圖2(c))進(jìn)行時(shí)間軸拉伸處理,獲得地震剖面如圖2(d)所示。在拉伸后的剖面中,繞射同相軸頂點(diǎn)“上移”,其附近變得“平緩”,側(cè)翼則更加“陡峭”。
根據(jù)式(2),在3個(gè)不同的速度下對(duì)拉伸剖面(圖2(d))進(jìn)行了等速延拓實(shí)驗(yàn),成像結(jié)果在還原時(shí)間軸后如圖3所示。當(dāng)延拓速度vvc遠(yuǎn)小于最佳成像速度v0時(shí),繞射同相軸兩翼變短內(nèi)收(圖3(a));vvc接近v0時(shí),收斂于頂點(diǎn)(圖3(b));隨著vvc的繼續(xù)增大,繞射弧朝反向重新展開(圖3(c))。延拓速度的連續(xù)變化使繞射同相軸兩翼持續(xù)產(chǎn)生相移,而在該過程中,其頂點(diǎn)位置保持靜止。
圖3 不同速度下的等速延拓結(jié)果比較Fig.3 Comparison of constant velocity continuation results with different imaging velocities
上述實(shí)驗(yàn)說明等速延拓可有效進(jìn)行繞射波偏移,其成像的關(guān)鍵在于延拓速度的選取。只有在最佳成像速度附近進(jìn)行延拓才能使繞射波準(zhǔn)確成像,否則將導(dǎo)致其欠偏移(速度太小)或過偏移(速度太大)。值得注意的是,對(duì)于等速延拓而言,最佳的延拓速度并非精確對(duì)應(yīng)于真實(shí)速度,原因或與剖面時(shí)間軸拉伸有關(guān)。
(3)路徑積分速度延拓成像。早先的速度延拓成像是通過一系列離散的速度進(jìn)行等速延拓,然后優(yōu)選高聚焦量的部分拼接形成最終的成像剖面。這種做法不僅運(yùn)算成本高,而且成像質(zhì)量還會(huì)受速度離散化程度的影響,路徑積分為該問題的解決創(chuàng)造了條件,通過積分解可實(shí)現(xiàn)速度區(qū)間內(nèi)等速延拓剖面的連續(xù)疊加。
通過式(3)構(gòu)造的路徑積分偏移濾波器以及相應(yīng)的成像結(jié)果如圖4所示,速度區(qū)間取700~2 500 m/s。圖4(a),(b)分別為其對(duì)應(yīng)的f-k幅值譜、相位譜,濾波器幅值譜中的黑色區(qū)域稱為“主瓣”,兩側(cè)伴有若干條淺灰色細(xì)長“旁瓣”,它們關(guān)于k=0對(duì)稱,且波數(shù)寬度隨著頻率的升高呈非線性增加;圖4(c)為繞射波成像剖面,經(jīng)過偏移后繞射弧頂點(diǎn)附近的同相軸能量得到了顯著增強(qiáng),側(cè)翼也得到一定程度的壓制,但其對(duì)側(cè)卻出現(xiàn)了一條與之共享繞射頂點(diǎn)的“反向弧”,它們形成的“X”狀同相軸對(duì)應(yīng)4條未完全收斂的“尾線”;圖4(d)給出了相應(yīng)的幅值譜響應(yīng)(Ω=375 Hz),主、旁瓣光滑可辨,其幅度和波數(shù)寬度隨著遠(yuǎn)離零波數(shù)軸衰減強(qiáng)烈;圖4(e)為濾波器相位譜響應(yīng),其值從k=0時(shí)的零相位向兩側(cè)逐漸減小到 -π/2,并且在|k|持續(xù)增大的過程中發(fā)生若干次正負(fù)“跳躍”,整體呈“鋸齒”狀,換言之,主瓣區(qū)相位穩(wěn)定,而旁瓣的相位變化則比較劇烈。
圖4 路徑積分濾波器和繞射波成像結(jié)果Fig.4 Path integral filter and diffraction imaging result
圖4(c)中的“反向弧”(或稱為尾線)是路徑積分法固有的偏移“剩余”,它與濾波器旁瓣的存在密切相關(guān)。圖4(f)為地震剖面濾波前后的幅值譜響應(yīng)(Ω=375 Hz),相較于原始曲線(灰色線),濾波后曲線的幅值(黑色線)在零波數(shù)附近有一個(gè)明顯的增強(qiáng),遠(yuǎn)離k=0處則被大幅削弱。然而,由于旁瓣(特別是圖4(d)中箭頭處的“第1旁瓣”)的“泄漏”效應(yīng),導(dǎo)致曲線主峰兩側(cè)出現(xiàn)了明顯的“旁峰”(圖4(f)中箭頭所指),這些旁峰在變換回t-x域時(shí),可能共同產(chǎn)生了偏移剖面中的尾線。從疊加角度來說,給定速度范圍的開區(qū)間內(nèi),繞射雙曲線的側(cè)翼在連續(xù)等速延拓偏移求和的過程中相互抵消,而區(qū)間端點(diǎn)處的偏移卻因無法抵消而被保留,從而導(dǎo)致成像剖面中產(chǎn)生2條對(duì)應(yīng)積分上、下限的同相軸。
(4)加權(quán)路徑積分速度延拓成像。在積分區(qū)間較大的情況下,經(jīng)過路徑積分法(式(3))偏移的繞射弧頂點(diǎn)和尾線存在較大的幅值差異(圖4(f)),因此后者產(chǎn)生的影響基本可以忽略;但區(qū)間較小時(shí),尾線可能與其他同相軸發(fā)生相長干涉而產(chǎn)生虛假特征,需要采取一定的措施來壓制這類干擾。為此,構(gòu)造了高斯加權(quán)路徑積分偏移濾波器(式(4))來削弱尾線的影響。取速度為700~2 500 m/s,vbias=1 600 m/s,σ=200 m/s,構(gòu)建的濾波器f-k幅值譜、相位譜及相應(yīng)譜的響應(yīng)(Ω=375 Hz)分別如圖5(a),(b),(d),(e)所示。與未加權(quán)的濾波器(圖4(a),(b),(d),(e))相比較,高斯加權(quán)有效壓制了濾波器的幅值譜旁瓣,并使其相位在遠(yuǎn)離零波數(shù)軸的區(qū)域“失穩(wěn)”,降低了大傾角分量貢獻(xiàn)的“穩(wěn)定疊加”。圖5(c)為高斯加權(quán)路徑積分法繞射波成像剖面,不僅有效弱化了尾線,還使繞射點(diǎn)得到了進(jìn)一步的聚焦。從濾波前、后剖面的幅值譜響應(yīng)(Ω=375 Hz)(圖5(f))也可看出,兩側(cè)的幅值被削弱,而中間的小波數(shù)區(qū)域得到了保留。此外,在原始繞射剖面中加入的隨機(jī)噪聲,偏移后表現(xiàn)為細(xì)微的水平同相軸,對(duì)成像結(jié)果并無明顯影響,說明該方法具有一定的抗噪能力。
圖5 高斯加權(quán)路徑積分濾波器和繞射波成像結(jié)果Fig.5 Gaussian weighted path integral filter and diffraction imaging result
(5)成像參數(shù)控制。高斯加權(quán)路徑積分法成像的質(zhì)量受控于積分區(qū)間[va,vb]、速度偏量vbias和標(biāo)準(zhǔn)差σ等參數(shù)。為此,基于單點(diǎn)繞射模型就成像參數(shù)的優(yōu)選開展了實(shí)驗(yàn)研究。
無論加權(quán)與否,路徑積分法的成像效果均受濾波器主瓣區(qū)的控制。理論上,增大積分范圍(速度區(qū)間)會(huì)減小主瓣的波數(shù)寬度,使其變窄,這意味著成像剖面中對(duì)應(yīng)區(qū)間端點(diǎn)的欠、過偏同相軸的形態(tài)將發(fā)生改變;同時(shí),適當(dāng)增大速度區(qū)間還可以提高主瓣幅值,進(jìn)而凸顯繞射頂點(diǎn)。然而,過大增加積分范圍將以犧牲運(yùn)算成本為代價(jià)。
當(dāng)σ限定時(shí),最佳成像速度與積分區(qū)間[va,vb]的包含關(guān)系將直接影響成像效果,區(qū)間整體低于或高于最佳速度時(shí)均無法正確成像。此外,積分下限越遠(yuǎn)離最佳速度,成像剖面中繞射點(diǎn)兩側(cè)的同相軸欠偏越嚴(yán)重(圖6(a));反之,則會(huì)導(dǎo)致嚴(yán)重過偏(圖6(c))。獲得良好成像結(jié)果的前提是必須使積分區(qū)間包含最佳成像速度(圖6(b)),此時(shí)尾線的長短受積分上、下限影響,而繞射點(diǎn)的幅值強(qiáng)弱則與速度區(qū)間的大小有關(guān)。
圖6 不同積分限的路徑積分成像結(jié)果比較Fig.6 Comparison of path integral imaging results with different velocity ranges
一般而言,對(duì)于確定的地震剖面,其成像速度往往分布于一個(gè)較大的區(qū)間內(nèi),但受目標(biāo)層的限制,必然存在一個(gè)更小的、分布更集中的子區(qū)間,可將這個(gè)含有主要成像速度的區(qū)間稱為“優(yōu)勢(shì)速度區(qū)間”。速度偏量vbias指定了優(yōu)勢(shì)速度區(qū)間所處的位置,在以該速度為中心的鄰域內(nèi),貢獻(xiàn)疊加的權(quán)值將被增強(qiáng)。在缺乏先驗(yàn)信息的情況下,vbias可取為積分區(qū)間的中點(diǎn),認(rèn)為優(yōu)勢(shì)速度區(qū)間分布在以該點(diǎn)為中心的一定范圍內(nèi)。
在進(jìn)行高斯加權(quán)路徑積分法成像時(shí),標(biāo)準(zhǔn)差σ
規(guī)定了優(yōu)勢(shì)速度區(qū)間內(nèi)貢獻(xiàn)的增強(qiáng)程度。由圖7可看出,在積分限va,vb和速度偏量vbias限定的情況下,標(biāo)準(zhǔn)差σ越小,vbias附近貢獻(xiàn)的權(quán)值就越高,繞射收斂越集中(圖7(b)),但太小的σ值則容易放大錯(cuò)誤速度引起的貢獻(xiàn),使尾線干擾明顯增強(qiáng)(圖7(c));當(dāng)σ值越大時(shí),剖面則越趨向于未加權(quán)路徑積分法的成像結(jié)果(圖7(a))。因此,適當(dāng)選擇σ值可有效提高繞射點(diǎn)的分辨率。
圖7 不同標(biāo)準(zhǔn)差σ的路徑積分成像結(jié)果比較(va=1 000 m/s,vb=2 200 m/s,vbias=1 600 m/s)Fig.7 Comparison of path integral imaging results with different standard deviation σ(va=1 000 m/s,vb=2 200 m/s,vbias=1 600 m/s)
當(dāng)路徑積分法應(yīng)用于實(shí)際時(shí),可先大致選取一個(gè)合理的速度范圍進(jìn)行偏移,良好的成像剖面應(yīng)具有均衡的欠、過偏成分,若存在明顯的欠偏同相軸,則需要提高積分下限va,而大范圍的過偏弧則說明積分上限vb過高。通過不斷調(diào)整積分限,就可獲得一個(gè)既包含主要成像速度、又能較好地突出繞射點(diǎn)的速度區(qū)間,然后以此為基礎(chǔ)確定速度偏量vbias和標(biāo)準(zhǔn)差σ,最后再利用高斯加權(quán)路徑積分法實(shí)現(xiàn)繞射波成像。
為了進(jìn)一步檢驗(yàn)對(duì)多繞射點(diǎn)剖面的成像效果,建立了一個(gè)含有5個(gè)繞射源(標(biāo)識(shí)為a,b,c,d,e)的地質(zhì)模型(圖8(a)),其中均勻地層的速度為3 000 m/s,異質(zhì)體速度為4 500 m/s,半徑為8 m,采用與單繞射點(diǎn)模型相同的模擬方法和數(shù)據(jù)采集參數(shù)獲得如圖8(b)所示的地震記錄(剖面中加入了隨機(jī)噪聲,SNR=10),利用路徑積分速度延拓對(duì)該剖面進(jìn)行偏移成像,取va=1 500 m/s,vb=4 600 m/s,vbias=3 050 m/s、σ=200 m/s。從成像結(jié)果來看(圖8(c)),5個(gè)異質(zhì)體均得到了準(zhǔn)確成像,說明該方法適用于多繞射源成像。
圖8 多點(diǎn)繞射模型成像Fig.8 Seismic imaging of multi-point diffraction model
實(shí)際地震資料的原始水平疊加剖面如圖9(a)所示,該剖面構(gòu)造較復(fù)雜,0~0.5 km存在一個(gè)尖滅,由此形成多條繞射波同相軸;在整條剖面(100~300 ms)中存在一套不整合面,由于剝蝕程度的差異,導(dǎo)致地震剖面上也存在數(shù)條繞射同相軸,能量較弱,幾乎淹沒于反射波之下;在剖面的右側(cè),存在一組空間分布較寬(1.2~1.6 km)、深度跨度較長(100~800 ms)的繞射波,僅靠水平疊加剖面還難以判斷該繞射源的形態(tài)。反射波疊后偏移成像剖面(圖9(b))中繞射波得到了有效收斂,尤其是右側(cè)的繞射波組收斂成一條近乎連續(xù)的波阻抗界面,結(jié)合該界面兩側(cè)錯(cuò)斷的同相軸可初步判斷該處發(fā)育一條較大型的逆斷層,但受強(qiáng)反射能量的影響,尖滅點(diǎn)、斷點(diǎn)、不整合面的凸點(diǎn)等信息幾乎被掩蓋在反射同相軸之下,縱、橫向分辨率較低。
為了更好地揭示異質(zhì)體(區(qū))的特征,提取了地震資料中的繞射波場(chǎng)(圖10(a)),并分別采用Kirchhoff法和路徑積分法對(duì)其進(jìn)行疊后偏移成像處理。與疊后Kirchhoff偏移成像結(jié)果(圖10(b))相比較,路徑積分法成像剖面(圖10(c))中,破碎帶、風(fēng)化殼等異質(zhì)區(qū)引起的繞射波得到了有效收斂,清晰地刻畫出構(gòu)造細(xì)節(jié),且縱、橫向分辨率以及風(fēng)化殼界面波的連續(xù)性均優(yōu)于疊后Kirchhoff偏移剖面。基于紊亂的地震響應(yīng)特征,可在剖面中圈劃出2個(gè)地質(zhì)異質(zhì)區(qū):300 ms以淺,清晰展示了由風(fēng)化殼引起的尖滅點(diǎn)、凸點(diǎn)信息;剖面右側(cè)則揭示出斷層及斷層引起的破碎異質(zhì)區(qū),而相關(guān)信息在反射波疊后偏移成像剖面(圖9(b))中則幾乎看不到。
圖9 實(shí)際資料全波場(chǎng)成像Fig.9 Full wave field imaging of real seismic data sets
(1)f-k域路徑積分法是濾波形式的疊加類成像方法,可在給定區(qū)間內(nèi)自動(dòng)累積穩(wěn)定貢獻(xiàn),無需預(yù)先提供成像速度模型,減少了速度建模時(shí)的主觀影響,具有流程簡單、精度高、穩(wěn)定性好、運(yùn)算效率高等優(yōu)點(diǎn)。
(2)以疊后繞射剖面作為f-k域路徑積分法的輸入時(shí),強(qiáng)反射波場(chǎng)對(duì)成像結(jié)果的影響較大,因此在提取繞射波時(shí)應(yīng)盡量保證其干凈和完整。
(3)高斯加權(quán)路徑積分法受控于積分區(qū)間[va,vb]、速度偏量vbias和標(biāo)準(zhǔn)差σ等參數(shù),需恰當(dāng)選擇以確保繞射波的成像質(zhì)量。