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

(1.西安石油大學 地球科學與工程學院，陜西 西安 710065；2.陜西省油氣成藏地質學重點實驗室，陜西 西安 710065)

煤田地球物理勘探的任務之一是準確識別和追蹤斷層、陷落柱、富水層等異質體(區(qū))，以避免(或減輕)開采過程中的地質災害[1-9]。地震繞射波由于具有較高的分辨力和與小尺度異質體的準確對應關系，近年來備受勘探人員重視。其實，早在20世紀50年代，KREY[10]等就已經發(fā)現(xiàn)來自斷塊、尖滅等處的繞射波，之后TROREY[11],BERRYHILL[12]和Klem-Musatov[13]等對其進行了研究并建立了基本的繞射波理論體系。隨著研究的不斷深入，人們發(fā)現(xiàn)除了能有效檢測異質體的存在，繞射波還可用于描述異質體的方位及邊界等特征[14]，充分利用繞射波對地下異質體(區(qū))的精細勘探具有十分重要的意義[15-18]。

繞射波利用的關鍵是對其進行準確成像，但獨立的時差規(guī)律使得它在與反射波同時成像時不僅得不到正確收斂，而且還會因能量較弱而被掩蓋。此外，繞射成像對速度擾動較為“敏感”[19]，很小的速度變化都可能使其離開最佳的收斂狀態(tài)，因此偏移時需要提供精確的速度模型，這也導致了繞射成像實現(xiàn)起來比較困難。雖然近年發(fā)展出許多專門針對繞射波的成像方法和技術，也取得了一定的成像效果[20-25]，但大多還是對速度模型有著較強的依賴[26]。

路徑積分法是一種較新的地震成像技術，因無需預先構建成像速度模型就可實現(xiàn)成像，所以引起了研究者的極大關注[27-30]。目前，該方法在直接成像中的應用基于稱為“速度延拓”的偏移算法。CLAERBOUT[31]提出了速度延拓偏微分方程，將其應用于速度分析；FOMEL[32-33]指出在各向同性情況下，零偏移距速度延拓是一個連續(xù)和解析的過程，并給出了方程的f-k域解；BURNETT等[34]認為路徑積分和速度延拓的結合拓展出一種不依賴速度模型的成像方法；LANDA等[35]由共成像點道集的平面度計算了權值，說明加權因子的引入可增強路徑積分成像的協(xié)調性；MERZLIKIN等[36]進行了路徑積分法繞射波成像，并且推導了偏移濾波器的解析形式。上述研究驗證了路徑積分成像的可行性，探索了其在地震勘探中的應用潛力。

筆者在已有理論的基礎上，研究針對疊后繞射波的f-k域路徑積分成像方法，以進一步改善繞射波的成像質量，提高破碎帶、巖溶、采空區(qū)等小尺度異質體(區(qū))的探測精度。

1 方法原理

1.1 路徑積分速度延拓

速度延拓(Velocity Continuation，VC)是一種以速度為傳播變量的波形過程[33]，描述了時間偏移同相軸在不同速度下的形態(tài)變化。零偏移距各向同性時，該過程可表示為一個線性雙曲形式的偏微分方程[31,33,36]：

(1)

式中，P(x,t,v)為地震波場，x，t，v分別為空間坐標、單程走時和半偏移速度。

當v=C時(C為非負常數(shù))，求解式(1)的過程稱為等速延拓，解得的P(x,t,C)即該速度下的恒速時間偏移剖面。一般可引入變量τ來替換單程走時t(τ=t2)，通過對時間軸的“拉伸”(重采樣)處理將方程簡化為傅里葉域中的常微分方程并得到其解析解：

(2)

基于速度延拓的連續(xù)性以及繞射同相軸頂點在該延拓過程中的“平穩(wěn)性”，可利用連續(xù)疊加等速延拓剖面的方式來實現(xiàn)繞射波成像[37-39]，繼而得到基于路徑積分的速度延拓偏移濾波器：

(3)

其中，F(xiàn)PI(Ω,k,va,vb)為f-k域路徑積分偏移濾波器，其性能取決于Ω，k以及速度積分限va和vb。通過式(3)得到的成像結果可認為是由給定區(qū)間內的所有等速延拓剖面等權重疊加而成。

為了進一步增加成像的協(xié)調性，提高成像質量，可在上述濾波器中加入高斯因子：

(4)

濾波器的解析表達式(式(4)第2個等號)基于復指數(shù)積分函數(shù)與虛誤差函數(shù)成正比的性質得出，通過求取虛誤差函數(shù)端點值的方式，就可快速求出相應的偏移濾波器。濾波結果經二維逆傅里葉變換后，可進一步輸出為成像剖面。與疊后t-x域的偏移不同，f-k域路徑積分法主要利用濾波運算進行幅值和相位控制，通過增強繞射波場中的水平(近水平)分量、弱化具有一定傾角的繞射弧側翼來實現(xiàn)繞射波的偏移成像。

1.2 成像流程

f-k域路徑積分法疊后繞射波偏移成像流程如圖1所示，具體實施步驟：

圖1 路徑積分法成像流程Fig.1 Flowchart of path integral method

(1)輸入原始地震記錄(水平疊加剖面)，對其進行預處理(能量補償、干擾波壓制等)，然后選擇合適的波場分離方法提取繞射波場Dx×y(x行y列)?？蛇x用SVD(Singular Value Decomposition)法提取繞射波[40-44]，基于水平疊加剖面中反射波和繞射波橫向相干性和連續(xù)性方面的明顯差異，通過SVD帶通濾波實現(xiàn)2者的分離。

(3)按照拉伸填充后的時空域剖面的網格參數(shù)生成坐標矩陣Ωm×1和k1×n，再選取積分限va和vb，由式(4)構造路徑積分偏移濾波器Fm×n(復矩陣)。

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

2.1 單點繞射模型

(1)模型構建。建立了一個單繞射點地質模型來檢驗方法的正確性(圖2(a))。如圖2(a)所示，速度為1 500 m/s的均勻地層中，存在1個半徑為15 m的球形異質體，球心距地表375 m，其速度為4 500 m/s。采用有限差分解波動方程的方法進行正演模擬，自激自收觀測系統(tǒng)，震源子波為Ricker子波，主頻為40 Hz，道距為2 m，道數(shù)為1 000道，采樣率為1 ms，采樣點數(shù)為1 000點。

正演獲得的地震剖面如圖2(b)所示，該剖面可等效為水平疊加剖面。從該地震剖面可看出，異質體引起的地震響應為一條具有雙曲線特征的繞射波同相軸，頂點在500 ms處，跨度約2 km，能量隨著遠離繞射源呈現(xiàn)出衰減特性。為了貼近實際情況，在地震剖面中加入了隨機噪聲(圖2(c))，其信噪比約為10∶1。

圖2 單點繞射模型及地震剖面Fig.2 Single-point diffraction model and seismic section

(2)等速延拓成像。由于路徑積分速度延拓成像基于等速延拓原理，因此需要先檢驗等速延拓的有效性。根據(jù)偏移流程，首先對原始數(shù)據(jù)(圖2(c))進行時間軸拉伸處理，獲得地震剖面如圖2(d)所示。在拉伸后的剖面中，繞射同相軸頂點“上移”，其附近變得“平緩”，側翼則更加“陡峭”。

根據(jù)式(2)，在3個不同的速度下對拉伸剖面(圖2(d))進行了等速延拓實驗，成像結果在還原時間軸后如圖3所示。當延拓速度vvc遠小于最佳成像速度v0時，繞射同相軸兩翼變短內收(圖3(a))；vvc接近v0時，收斂于頂點(圖3(b))；隨著vvc的繼續(xù)增大，繞射弧朝反向重新展開(圖3(c))。延拓速度的連續(xù)變化使繞射同相軸兩翼持續(xù)產生相移，而在該過程中，其頂點位置保持靜止。

圖3 不同速度下的等速延拓結果比較Fig.3 Comparison of constant velocity continuation results with different imaging velocities

上述實驗說明等速延拓可有效進行繞射波偏移，其成像的關鍵在于延拓速度的選取。只有在最佳成像速度附近進行延拓才能使繞射波準確成像，否則將導致其欠偏移(速度太小)或過偏移(速度太大)。值得注意的是，對于等速延拓而言，最佳的延拓速度并非精確對應于真實速度，原因或與剖面時間軸拉伸有關。

(3)路徑積分速度延拓成像。早先的速度延拓成像是通過一系列離散的速度進行等速延拓，然后優(yōu)選高聚焦量的部分拼接形成最終的成像剖面。這種做法不僅運算成本高，而且成像質量還會受速度離散化程度的影響，路徑積分為該問題的解決創(chuàng)造了條件，通過積分解可實現(xiàn)速度區(qū)間內等速延拓剖面的連續(xù)疊加。

通過式(3)構造的路徑積分偏移濾波器以及相應的成像結果如圖4所示，速度區(qū)間取700～2 500 m/s。圖4(a)，(b)分別為其對應的f-k幅值譜、相位譜，濾波器幅值譜中的黑色區(qū)域稱為“主瓣”，兩側伴有若干條淺灰色細長“旁瓣”，它們關于k=0對稱，且波數(shù)寬度隨著頻率的升高呈非線性增加；圖4(c)為繞射波成像剖面，經過偏移后繞射弧頂點附近的同相軸能量得到了顯著增強，側翼也得到一定程度的壓制，但其對側卻出現(xiàn)了一條與之共享繞射頂點的“反向弧”，它們形成的“X”狀同相軸對應4條未完全收斂的“尾線”；圖4(d)給出了相應的幅值譜響應(Ω=375 Hz)，主、旁瓣光滑可辨，其幅度和波數(shù)寬度隨著遠離零波數(shù)軸衰減強烈；圖4(e)為濾波器相位譜響應，其值從k=0時的零相位向兩側逐漸減小到 -π/2，并且在|k|持續(xù)增大的過程中發(fā)生若干次正負“跳躍”，整體呈“鋸齒”狀，換言之，主瓣區(qū)相位穩(wěn)定，而旁瓣的相位變化則比較劇烈。

圖4 路徑積分濾波器和繞射波成像結果Fig.4 Path integral filter and diffraction imaging result

圖4(c)中的“反向弧”(或稱為尾線)是路徑積分法固有的偏移“剩余”，它與濾波器旁瓣的存在密切相關。圖4(f)為地震剖面濾波前后的幅值譜響應(Ω=375 Hz)，相較于原始曲線(灰色線)，濾波后曲線的幅值(黑色線)在零波數(shù)附近有一個明顯的增強，遠離k=0處則被大幅削弱。然而，由于旁瓣(特別是圖4(d)中箭頭處的“第1旁瓣”)的“泄漏”效應，導致曲線主峰兩側出現(xiàn)了明顯的“旁峰”(圖4(f)中箭頭所指)，這些旁峰在變換回t-x域時，可能共同產生了偏移剖面中的尾線。從疊加角度來說，給定速度范圍的開區(qū)間內，繞射雙曲線的側翼在連續(xù)等速延拓偏移求和的過程中相互抵消，而區(qū)間端點處的偏移卻因無法抵消而被保留，從而導致成像剖面中產生2條對應積分上、下限的同相軸。

(4)加權路徑積分速度延拓成像。在積分區(qū)間較大的情況下，經過路徑積分法(式(3))偏移的繞射弧頂點和尾線存在較大的幅值差異(圖4(f))，因此后者產生的影響基本可以忽略；但區(qū)間較小時，尾線可能與其他同相軸發(fā)生相長干涉而產生虛假特征，需要采取一定的措施來壓制這類干擾。為此，構造了高斯加權路徑積分偏移濾波器(式(4))來削弱尾線的影響。取速度為700～2 500 m/s,vbias=1 600 m/s，σ=200 m/s，構建的濾波器f-k幅值譜、相位譜及相應譜的響應(Ω=375 Hz)分別如圖5(a)，(b)，(d)，(e)所示。與未加權的濾波器(圖4(a)，(b)，(d)，(e))相比較，高斯加權有效壓制了濾波器的幅值譜旁瓣，并使其相位在遠離零波數(shù)軸的區(qū)域“失穩(wěn)”，降低了大傾角分量貢獻的“穩(wěn)定疊加”。圖5(c)為高斯加權路徑積分法繞射波成像剖面，不僅有效弱化了尾線，還使繞射點得到了進一步的聚焦。從濾波前、后剖面的幅值譜響應(Ω=375 Hz)(圖5(f))也可看出，兩側的幅值被削弱，而中間的小波數(shù)區(qū)域得到了保留。此外，在原始繞射剖面中加入的隨機噪聲，偏移后表現(xiàn)為細微的水平同相軸，對成像結果并無明顯影響，說明該方法具有一定的抗噪能力。

圖5 高斯加權路徑積分濾波器和繞射波成像結果Fig.5 Gaussian weighted path integral filter and diffraction imaging result

(5)成像參數(shù)控制。高斯加權路徑積分法成像的質量受控于積分區(qū)間[va,vb]、速度偏量vbias和標準差σ等參數(shù)。為此，基于單點繞射模型就成像參數(shù)的優(yōu)選開展了實驗研究。

無論加權與否，路徑積分法的成像效果均受濾波器主瓣區(qū)的控制。理論上，增大積分范圍(速度區(qū)間)會減小主瓣的波數(shù)寬度，使其變窄，這意味著成像剖面中對應區(qū)間端點的欠、過偏同相軸的形態(tài)將發(fā)生改變；同時，適當增大速度區(qū)間還可以提高主瓣幅值，進而凸顯繞射頂點。然而，過大增加積分范圍將以犧牲運算成本為代價。

當σ限定時，最佳成像速度與積分區(qū)間[va,vb]的包含關系將直接影響成像效果，區(qū)間整體低于或高于最佳速度時均無法正確成像。此外，積分下限越遠離最佳速度，成像剖面中繞射點兩側的同相軸欠偏越嚴重(圖6(a))；反之，則會導致嚴重過偏(圖6(c))。獲得良好成像結果的前提是必須使積分區(qū)間包含最佳成像速度(圖6(b))，此時尾線的長短受積分上、下限影響，而繞射點的幅值強弱則與速度區(qū)間的大小有關。

圖6 不同積分限的路徑積分成像結果比較Fig.6 Comparison of path integral imaging results with different velocity ranges

一般而言，對于確定的地震剖面，其成像速度往往分布于一個較大的區(qū)間內，但受目標層的限制，必然存在一個更小的、分布更集中的子區(qū)間，可將這個含有主要成像速度的區(qū)間稱為“優(yōu)勢速度區(qū)間”。速度偏量vbias指定了優(yōu)勢速度區(qū)間所處的位置，在以該速度為中心的鄰域內，貢獻疊加的權值將被增強。在缺乏先驗信息的情況下，vbias可取為積分區(qū)間的中點，認為優(yōu)勢速度區(qū)間分布在以該點為中心的一定范圍內。

在進行高斯加權路徑積分法成像時，標準差σ

規(guī)定了優(yōu)勢速度區(qū)間內貢獻的增強程度。由圖7可看出，在積分限va,vb和速度偏量vbias限定的情況下，標準差σ越小，vbias附近貢獻的權值就越高，繞射收斂越集中(圖7(b))，但太小的σ值則容易放大錯誤速度引起的貢獻，使尾線干擾明顯增強(圖7(c))；當σ值越大時，剖面則越趨向于未加權路徑積分法的成像結果(圖7(a))。因此，適當選擇σ值可有效提高繞射點的分辨率。

圖7 不同標準差σ的路徑積分成像結果比較(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)

當路徑積分法應用于實際時，可先大致選取一個合理的速度范圍進行偏移，良好的成像剖面應具有均衡的欠、過偏成分，若存在明顯的欠偏同相軸，則需要提高積分下限va，而大范圍的過偏弧則說明積分上限vb過高。通過不斷調整積分限，就可獲得一個既包含主要成像速度、又能較好地突出繞射點的速度區(qū)間，然后以此為基礎確定速度偏量vbias和標準差σ，最后再利用高斯加權路徑積分法實現(xiàn)繞射波成像。

2.2 多點繞射模型

為了進一步檢驗對多繞射點剖面的成像效果，建立了一個含有5個繞射源(標識為a，b，c，d，e)的地質模型(圖8(a))，其中均勻地層的速度為3 000 m/s，異質體速度為4 500 m/s，半徑為8 m，采用與單繞射點模型相同的模擬方法和數(shù)據(jù)采集參數(shù)獲得如圖8(b)所示的地震記錄(剖面中加入了隨機噪聲，SNR=10)，利用路徑積分速度延拓對該剖面進行偏移成像，取va=1 500 m/s,vb=4 600 m/s，vbias=3 050 m/s、σ=200 m/s。從成像結果來看(圖8(c))，5個異質體均得到了準確成像，說明該方法適用于多繞射源成像。

圖8 多點繞射模型成像Fig.8 Seismic imaging of multi-point diffraction model

3 實際資料處理與分析

實際地震資料的原始水平疊加剖面如圖9(a)所示，該剖面構造較復雜，0～0.5 km存在一個尖滅，由此形成多條繞射波同相軸；在整條剖面(100～300 ms)中存在一套不整合面，由于剝蝕程度的差異，導致地震剖面上也存在數(shù)條繞射同相軸，能量較弱，幾乎淹沒于反射波之下；在剖面的右側，存在一組空間分布較寬(1.2～1.6 km)、深度跨度較長(100～800 ms)的繞射波，僅靠水平疊加剖面還難以判斷該繞射源的形態(tài)。反射波疊后偏移成像剖面(圖9(b))中繞射波得到了有效收斂，尤其是右側的繞射波組收斂成一條近乎連續(xù)的波阻抗界面，結合該界面兩側錯斷的同相軸可初步判斷該處發(fā)育一條較大型的逆斷層，但受強反射能量的影響，尖滅點、斷點、不整合面的凸點等信息幾乎被掩蓋在反射同相軸之下，縱、橫向分辨率較低。

為了更好地揭示異質體(區(qū))的特征，提取了地震資料中的繞射波場(圖10(a))，并分別采用Kirchhoff法和路徑積分法對其進行疊后偏移成像處理。與疊后Kirchhoff偏移成像結果(圖10(b))相比較，路徑積分法成像剖面(圖10(c))中，破碎帶、風化殼等異質區(qū)引起的繞射波得到了有效收斂，清晰地刻畫出構造細節(jié)，且縱、橫向分辨率以及風化殼界面波的連續(xù)性均優(yōu)于疊后Kirchhoff偏移剖面。基于紊亂的地震響應特征，可在剖面中圈劃出2個地質異質區(qū)：300 ms以淺，清晰展示了由風化殼引起的尖滅點、凸點信息；剖面右側則揭示出斷層及斷層引起的破碎異質區(qū)，而相關信息在反射波疊后偏移成像剖面(圖9(b))中則幾乎看不到。

圖9 實際資料全波場成像Fig.9 Full wave field imaging of real seismic data sets

4 結 論

(1)f-k域路徑積分法是濾波形式的疊加類成像方法，可在給定區(qū)間內自動累積穩(wěn)定貢獻，無需預先提供成像速度模型，減少了速度建模時的主觀影響，具有流程簡單、精度高、穩(wěn)定性好、運算效率高等優(yōu)點。

(2)以疊后繞射剖面作為f-k域路徑積分法的輸入時，強反射波場對成像結果的影響較大，因此在提取繞射波時應盡量保證其干凈和完整。

(3)高斯加權路徑積分法受控于積分區(qū)間[va,vb]、速度偏量vbias和標準差σ等參數(shù)，需恰當選擇以確保繞射波的成像質量。