黏性介質(zhì)疊前時間偏移方法

李雪英，文慧儉，張 桐，于國楠，王金麗

（1．東北石油大學地球科學學院，黑龍江大慶 163318； 2．大慶油田博士后科研工作站，黑龍江大慶 163458； 3．東北石油大學博士后科研流動站，黑龍江大慶 163318； 4．大慶油田有限責任公司第三采油廠，黑龍江大慶163113； 5．大慶油田有限責任公司第二采油廠，黑龍江大慶 163416）

黏性介質(zhì)疊前時間偏移方法

李雪英1，2，3，文慧儉1，張 桐4，于國楠5，王金麗5

（1．東北石油大學地球科學學院，黑龍江大慶 163318； 2．大慶油田博士后科研工作站，黑龍江大慶 163458； 3．東北石油大學博士后科研流動站，黑龍江大慶 163318； 4．大慶油田有限責任公司第三采油廠，黑龍江大慶163113； 5．大慶油田有限責任公司第二采油廠，黑龍江大慶 163416）

根據(jù)衡Q模型，將吸收系數(shù)引入黏性介質(zhì)的相速度中，以相移法為基礎，利用穩(wěn)相點原理，推導基于疊加速度和等效Q值的黏性介質(zhì)疊前時間偏移的走時和振幅計算公式，研究適合于黏性介質(zhì)的疊前時間偏移方法．理論數(shù)據(jù)處理結(jié)果表明，該方法具有精確的復雜構(gòu)造成像能力，能夠?qū)ば匝a償結(jié)合到偏移過程中，有效恢復深層界面的高頻信息；在準確構(gòu)造成像的同時，提高地震成像的分辨率．

疊前時間偏移；黏性吸收補償；等效Q值；穩(wěn)定性控制；黏性介質(zhì)

0 引言

在常規(guī)疊前時間偏移過程中，忽略地震波能量的吸收與頻散效應是不合理的．因此，人們提出應該在剔除地震波傳播效應的同時，完成對地震波能量補償?shù)酿ば越橘|(zhì)偏移．黏性補償技術(shù)包括兩類：一類是反Q濾波算法；另一類是黏性介質(zhì)深度偏移算法．在早期的黏性補償技術(shù)中，反Q濾波算法得到快速發(fā)展．Hargreaves N D等提出一種與Stolt偏移算法類似的反Q濾波算法，通過引入類似于Stolt偏移的坐標變換，利用快速Fourier變換，可以有效地校正時變的相位扭曲［1—2］．基于地震波場延拓理論，Wang Y H提出一種適用于水平層狀和Q值隨時間或深度連續(xù)變化，并且可以同時完成幅度補償和相位校正的反Q濾波算法［3—5］．反Q濾波技術(shù)無法處理Q值橫向變化的疊前地震資料，因此在地震資料處理中的應用價值受到限制．

為了能夠?qū)ば匝a償技術(shù)推向疊前地震資料處理中，人們提出黏性介質(zhì)疊前深度偏移算法．Mitte R等根據(jù)黏性介質(zhì)疊前深度偏移的基本思想，建立黏性介質(zhì)的聲波波動方程及其波場延拓方法，給出適用于黏性介質(zhì)偏移的成像條件，并進行相應的數(shù)值算例驗證［6—7］．Causse E等實現(xiàn)黏性介質(zhì)的逆時偏移成像［8—9］．Emmerich H、Carcione J M和Stekl I等在地震正演模擬中，以及Keers H等在反演中考慮地震波能量被吸收的實際效果，模擬地震波在地下真實傳播情況［10—13］．Zhang J F等采用加權(quán)最小平方法，設計最優(yōu)空間褶積短算子適應介質(zhì)橫向變化，通過限制波場最大傳播角度的策略，實現(xiàn)算法穩(wěn)定性，給出黏性非均勻介質(zhì)疊前深度偏移的精確成像，提高算法的精度與效率［14］．這些研究主要集中在對黏性介質(zhì)疊前深度偏移的論述上，對于黏性介質(zhì)的疊前時間偏移方法鮮有論及．

筆者根據(jù)穩(wěn)相點原理，推導黏性介質(zhì)疊前時間偏移的高頻漸進解公式，引入等效Q值，使得黏性介質(zhì)疊前時間偏移方法可以應用到疊前地震資料的高分辨率成像處理中，解決反Q濾波不能處理Q值非均勻的疊前資料及黏性介質(zhì)疊前深度偏移層狀Q值建模問題．

1 理論推導

1．1 走時及振幅

基于單程波方程和穩(wěn)相點理論，參考深度偏移方法研究黏性介質(zhì)疊前時間偏移中走時及振幅計算方法［15—19］．對于水平層狀介質(zhì)，炮和單道數(shù)據(jù)可以統(tǒng)一看作空間脈沖函數(shù)，在頻率—波數(shù)域用f（ω）表達，基于深度偏移相移法，對于水平層狀黏性介質(zhì)波場的深度延拓［15，20］可以表示為

式中：j為虛數(shù)單位；ω為頻率；kx為水平波數(shù)；P（kx，ω，z）為深度為z處的頻率—波數(shù)域波場；Δzi為各層介質(zhì)的厚度；n為目的層以上的介質(zhì)層數(shù)；f（ω）為炮點或檢波點的時域信號的傅立葉變換；vi（ω）為各層復速度，且

式中：vi為各層的實相速度；Qi為各層的品質(zhì)因子．用時間表示波場延拓步長，對于黏性介質(zhì)，在地面波場延拓到T深度的波場為

式中：T為單程旅行時表達深度；Δτi為各層的時間厚度，為深度T處頻率—波數(shù)域波場．

令px＝kx／ω，則式（3）的指數(shù)函數(shù)的指數(shù)項可近似表達為

令

式中：vrms為均方根速度；Qeff為等效Q值；vQrms為與Q有關(guān)的速度，在實際應用時可以近似等于均方根速度．

將式（5）和式（6）代入式（3），做空間Fourier反變換，得到空間—頻率域波場：

式中：x為炮（檢波）點到成像點之間的水平距離；px為射線參數(shù)，與頻率無關(guān)，令

對式（9）求取一階和二階導數(shù)得

應用穩(wěn)相點原理，對式（7）求得黏性吸收介質(zhì)相移公式的漸進解：

式中：p0x為φ′（px）的零點，為成像點所對應的射線參數(shù)，即穩(wěn)相點．令φ′（px）＝0，解得：

將p0x代入式（11）得：

將p0x，φ″（p0x）代入式（9）、式（11）、式（8），得到地震波從源點（0，0）到（x，T）的走時t、振幅A和吸收Θ 為

當ω大于2Hz時，式（12）是式（7）的很好近似，對于地震勘探是滿足的，因此可以將式（12）作為頻率域—空間波場的解析表達式．

由式（12）定義下行波在吸收介質(zhì)的正向傳播算子為

其中下標B代表吸收介質(zhì)．

為了能夠推導吸收介質(zhì)中上行波的逆向傳播算子，假設存在一種與吸收介質(zhì)相伴隨的補償介質(zhì)，體積模量定義為

式中：＊為復共軛符號；M′（ω）為補償介質(zhì)的體積模量，與吸收介質(zhì)的體積模量互為復共軛．補償介質(zhì)只是為了計算方便而引入的一種與吸收介質(zhì)相伴隨的表達方式（下標A表示）．

Zhang J F給出吸收介質(zhì)和補償介質(zhì)上、下行波的正向傳播算子與逆?zhèn)鞑ニ阕又g的關(guān)系［14］：

由互易定理得：

同理，可得：

式中：x′為炮點或者檢波點坐標；x為成像點坐標；F+A、F+B為下行波在補償介質(zhì)和吸收介質(zhì)中的逆向傳播算子；F—A和F—B為上行波在補償介質(zhì)和吸收介質(zhì)中的逆向傳播算子；W+A、W+B為下行波在補償介質(zhì)和吸收介質(zhì)中的正向傳播算子；W—A、W—B為上行波在補償介質(zhì)和吸收介質(zhì)中的正向傳播算子．

由式（20）和式（22）或式（23）和式（24）給出吸收介質(zhì)逆向傳播算子和補償介質(zhì)正向傳播算子之間的轉(zhuǎn)換關(guān)系，即互為復共軛；反之，亦然．對于黏性介質(zhì)疊前時間偏移，吸收介質(zhì)中的上行波場向下延拓可以通過對波場乘以逆向傳播算子完成；下行波場的向下延拓根據(jù)不同的成像條件選擇不同的正向傳播算子；如果是在吸收介質(zhì)中延拓，則選取W+B（x，x′）；在補償介質(zhì)中延拓，則選擇W+A（x，x′）．

1．2 成像條件

有2種成像條件可以用于黏性介質(zhì)的偏移成像，第一種是反褶積成像條件：

式中：PD（x，ω）為吸收介質(zhì)中正向延拓的下行波場；PU（x，ω）為吸收介質(zhì)中反向延拓的上行波場；ε為一個很小穩(wěn)定性常數(shù)．董春暉分別考察單極源和偶極源采用褶積成像條件，可以實現(xiàn)疊前時間偏移的保幅處理［］．

第二種成像條件是在補償介質(zhì)中設計一個正向外推算子，以增強波場．該正向算子與吸收介質(zhì)的逆向傳播算子互為復共軛（見式（23）），因而很容易獲取算子．如果設［PD（x，ω）］′為補償介質(zhì)中的正向延拓下行波場，則第二種成像條件（相關(guān)成像條件）可以表達為

該成像條件首先由Mittet R等提出［6］，后由Causse E等給出嚴格的數(shù)學證明［9］．Zhang J F通過數(shù)值算例證實，褶積成像條件與相關(guān)成像條件相比可以很好地壓制噪聲，從而給出更加清晰的成像而實現(xiàn)保幅處理［14］，文中選用反褶積成像條件．

1．3 穩(wěn)定性控制

黏性介質(zhì)疊前時間偏移成像是通過設計補償算子完成地震數(shù)據(jù)的高頻補償，但是在高頻恢復過程中也會增強高頻噪聲，使算法表現(xiàn)固有的不穩(wěn)定性．因此，穩(wěn)定性控制成為黏性介質(zhì)偏移成像和反Q濾波最為關(guān)鍵的問題．

為了使噪聲在偏移過程中不出現(xiàn)不必要的擴大，設計一個適用于高頻的補償算子．重新構(gòu)造一個以自變量為χ（χ＝ωT）的補償函數(shù)Γ（χ），設置截止值χc，當χ≤χc時，補償函數(shù)為原有精確的補償算子；當χ＞χc時，補償算子經(jīng)過一個一階光滑可導函數(shù)改造后而漸變?yōu)樽畲蟮恼穹T限值Glim，經(jīng)過改造后的補償函數(shù)Γ（χ）隨著χ變化與精確補償算子對比結(jié)果見圖1（實線為原有的補償算子的模，虛線為改造后的模）．該方法可以保持高頻的振幅恢復，對深層和遠偏移距的很弱有效信息進行最大限度恢復與補償，雖然在改造過程中不可避免地增強背景噪聲，但是噪聲在疊前偏移后的共反射點疊加中得到有效壓制．

圖1 補償算子與經(jīng)過改造后的補償算子譜Fig．1 The comparison of compensation operator spec—trum with spectrum after smoothing

式中：a（χ）為一階光滑可導函數(shù)．

2 模型驗證

采用二維頻域有限差分方法雙程聲波波動方程，建立黏性介質(zhì)的衰減正演模擬理論模型．正演模擬的宏觀介質(zhì)參數(shù)模型見圖2．模擬后的數(shù)據(jù)包含51個共炮點道集，震源采用的是最大頻率為55Hz的雷克子波，炮點分布位置為447．75～3 147．75m，炮間距為54m；每炮為800道，檢波點間距為4．5m；記錄長度為1．68s，采樣間隔為4ms；在模型兩側(cè)每邊共設20個網(wǎng)格點的衰減區(qū)，以減少邊界反射．為了保證模擬數(shù)據(jù)更接近真實數(shù)據(jù)，在模擬中保留表面多次波和層間多次波．

圖2 正演模擬的宏觀介質(zhì)參數(shù)模型Fig．2 Subsurface real phase velocity and Q model

采用共反射點疊加方案，對模擬的理論數(shù)據(jù)進行常規(guī)疊前時間偏移和黏性介質(zhì)疊前時間偏移（見圖3和圖4）．由圖3和圖4可以看出，2個成像剖面上所有反射軸正確歸位且精確成像；在模型右側(cè)的楔狀體處，波從高速介質(zhì)入射到低速介質(zhì)，波形出現(xiàn)極性反轉(zhuǎn)．在常規(guī)疊前時間偏移成像剖面上，由于地震波在傳播過程中存在衰減和頻散，地震波的能量隨著深度增加而逐漸減弱；使得深層反射界面成像模糊，振幅減弱，分辨率較差．為了更好地顯示分辨率改善效果，在增益后的圖像中任意選擇5處設定標尺，然后在黏性補償剖面中按照相同的長度畫出，在對吸收給出正確補償后（見圖4），所有界面得到正確構(gòu)造歸位，深層記錄頻帶得到拓寬（見圖5），不考慮黏性補償?shù)臅r頻譜見圖5（a），考慮黏性補償?shù)臅r頻譜見圖5（b）．由圖5可以看出，空間子波被壓縮，振幅增強，分辨率提高．

圖3 常規(guī)疊前時間偏移結(jié)果Fig．3 The profile of conventional prestack time migration

圖4 考慮黏性補償?shù)寞B前時間偏移結(jié)果Fig．4 The prestack time migration profile of compensation for anelastic absorption

圖5 廣義S變換前后時頻譜Fig．5 The comparison time frequency spectrum gained by generalized S transform

從2個成像剖面同一位置處抽取一道數(shù)據(jù)，利用廣義S變換進行時頻分析（見圖5）．在不考慮黏性補償?shù)臅r頻譜上，隨著深度的增加，地層子波的頻帶逐漸變窄，峰值頻率向低頻移動．經(jīng)過黏性補償后，深層的子波頻帶得到恢復，峰值頻率向高頻移動，深淺層的子波頻帶寬度保持一致．

理論數(shù)據(jù)處理結(jié)果表明：文中給出的方法在對非均勻吸收介質(zhì)進行正確構(gòu)造成像的同時，可以完成保幅、提高分辨率的地震處理要求，在成像過程中體現(xiàn)考慮黏性補償?shù)闹匾裕?/p>

3 結(jié)束語

文中給出的方法將黏性補償結(jié)合到偏移過程中，能夠很好地完成對吸收與頻散的補償；深層界面信息得到有效恢復，頻帶拓寬，振幅增強；地震剖面分辨率得到明顯提高．在精確構(gòu)造成像的基礎上，能夠豐富地震信息，為后續(xù)的地震解釋奠定基礎．

［1］ Hargreaves N D，Calvert A J．Inverse Q filtering by Fourier transform［J］．Geophysics，1991，56（1）：519—527．

［2］ Hargreaves N D．Similarity and the inverse Q filter：some simple algorithms for inverse Q filtering［J］．Geophysics，1992，57（2）：944 —947．

［3］ Wang Y H．A stable and efficient approach of inverse Q filtering［J］．Geophysics，2002，67（2）：657—663．

［4］ Wang Y H．Quantifying the effectiveness of stabilized inverse Q filtering［J］．Geophysics，2003，68（1）：337—345．

［5］ Wang Y H．Inverse Q—filter for seismic resolution enhancement［J］．Geophysics，2006，71（3）：V51—V60．

［6］ Mittet R，Sollie R，Hokstad K．Prestack depth migration with compensation for absorption and dispersion［J］．Geophysics，1995，60 （5）：1485—1494．

［7］ Mittet R．A simple design procedure for depth extrapolation operators that compensate for absorption and dispersion［J］．Geophysics，2007，72（2）：S105—S102．

［8］ Causse E，Ursin B．Asymptotic error analysis of constant—velocity viscoacoustic migration［J］．Geophysics，1999，64（4）：1036—1045．

［9］ Causse E，Ursin B．Viscoacoustic reverse—time migration［J］．Journal of Seismic Exploration，2000，9（1）：165—184．

［10］ Emmerich H，Korm M．Incorporation of attenuation into time—domain computations of seismic wavefields［J］．Geophysics，1987，52 （9）：1252—1264．

［11］ Carcione J M，Kosloff D，kosloff R．Wave propagation simulation in a linear viscoelastic medium［J］．Geophysics Journal Interna—tional，1999，95（3）：597—611．

［12］ Stekl I，Pratt R G．Accurate viscoelastic modeling by frequency—domain finite differences using rotated operators［J］．Geophysics，1998，63（5）：1779—1794．

［13］ Keers H，Vasco D W，Johnson L R．Viscoacoustic cross—well imaging using asymptotic waveforms［J］．Geophysics，2001，66（3）：861—870．

［14］ Zhang J F，Wapenaar K．Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media［J］．Geophysical prospecting，2002，50（6）：629—643．

［15］ Gazdag J．Wave equation migration with the phase—shift method［J］．Geophysics，1978，43（7）：1342—1351．

［16］ Bleistein N，Cohen J K，Stockwell J W．Mathematics of multidimensional seismic inversion［M］．Berlin：Springer，2001．

［17］ Zhang J F，Wapenaar C P A．High—resolution depth imaging with sparseness—constrained inversion［J］．Geophysical Prospecting，2006，54（1）：49—62．

［18］ Zhang Jianfeng，Liu Linong．Optimum split—step Fourier 3Ddepth migration：Developments and practical aspects［J］．Geophysics，2007，72（3）：S167—S175．

［19］ Claerbout J F．Imaging the earth＇s interior［M］．Oxford：Blackwell Scientific Publications，1985．

［20］ 董春暉，張劍鋒．起伏地表下的直接疊前時間偏移［J］．地球物理學報，2009，52（1）：239—244．Dong Chunhui，Zhang Jianfeng．Prestack time migration including surface topography［J］．Chinese Journal of Geophysics，2009，52 （1）：239—244．

TE132．1

A

2095—4107（2013）04—0098—07

DOI 10．3969／j．issn．2095—4107．2013．04．015

2013—07—05；編輯：任志平

國家科技重大專項（2011ZX05008—006—42）；黑龍江省博士后基金項目（LBH—Z12002）；黑龍江省教育廳科技項目（12513069）

李雪英（1972—），男，博士，教授，主要從事地震波傳播與成像方面的研究．

文慧儉（1972—），女，博士，副教授，主要從事油田開發(fā)地質(zhì)方面的研究，whjdqpi＠163．com．