李 勤，沈鴻雁，王 鑫，趙 靜，李 萌

(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è)精度。

1 方法原理

1.1 路徑積分速度延拓

速度延拓(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)繞射波的偏移成像。

1.2 成像流程

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ù)矩陣)。

2 模型數(shù)據(jù)處理與分析

2.1 單點(diǎn)繞射模型

(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)繞射波成像。

2.2 多點(diǎ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

3 實(shí)際資料處理與分析

實(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

4 結(jié) 論

(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ì)量。