李 勤， 李慶春， 張 林

（1.西安科技大學地質與環(huán)境學院，陜西西安 710054；2.長安大學地質工程與測繪學院，陜西西安 710054 3.陜西省地質調查中心，陜西西安 710016）

基于射線理論的各向異性煤層地震響應

李 勤1， 李慶春2， 張 林3

（1.西安科技大學地質與環(huán)境學院，陜西西安 710054；2.長安大學地質工程與測繪學院，陜西西安 710054 3.陜西省地質調查中心，陜西西安 710016）

以VTI型煤層為研究對象，從介質中的地震波的射線方向、走時計算等多個方面展開討論，分析VTI型煤層的試射法射線追蹤技術，充分利用VC＋＋的面向對象類框架泛型編程的特點和Fortran數值計算方面的優(yōu)勢，通過混合編程實現(xiàn)VTI型煤層中的地震波試射法射線追蹤；通過各向同性煤層和VTI型煤層的射線追蹤試算對比，論證對VTI型煤層進行各向異性射線追蹤的必要性，同時驗證本文VTI型煤層射線追蹤方法的正確性.最后，通過計算時間成本比較，表明該方法計算速度較快，能滿足地震正演模擬的要求.該方法適用于各向異性煤層的多波射線追蹤.

VTI型煤層；射線追蹤；各向異性

0 引言

煤層裂隙是國內外煤礦地震勘探研究的重點和難點，裂隙發(fā)育復雜［1－3］，其中一種較為典型的裂隙系統(tǒng)可近似為平行于地表的裂隙系統(tǒng)，這種裂隙系統(tǒng)稱為VTI型，而這種煤層就稱之為VTI型煤層［4］.VTI型煤層作為各向異性煤層的一種常見模型，裂隙發(fā)育導致的各向異性會對地震波的記錄產生較大影響［5－6］.但是，由于煤層相對地震子波來說較薄，各向異性煤層地震波正演模擬也相對復雜.目前，對于煤層進行各向異性的研究主要集中在實驗室分析［7－8］，波動方程法地震波模擬［9－14］及AVO特征分析［15－17］.射線追蹤法和波動方程法各有優(yōu)缺點.波動方程法相對來說，所生成的地震記錄波形豐富，信息量大，比較接近實際地震記錄；而射線追蹤法適合于追蹤單一波形，所得模擬地震記錄可直接用于地震波各向異性響應特征研究.射線追蹤法地震波模擬由于其具有顯示直觀、高效計算的特點，在石油地震勘探方面取得了顯著成效［18－21］，然而卻一直未見用于各向異性煤層中.試射法是一種基于射線的算法，也是一種經典的射線追蹤算法，有很大的發(fā)展空間［22］.目前的地震記錄正演模擬軟件很多，如Tesseral等，它們大多是對層厚較大的模型進行模擬，還沒有對VTI型煤層的正演模擬模塊.

本文以各向異性煤層為研究對象，采用試射法對煤層中地震波進行射線追蹤，分析VTI型煤層中走時及射線路徑的傳播特征；并將VC＋＋與Fortran語言有機結合，實現(xiàn)混合編程，研發(fā)基于VTI型煤層的地震記錄正演模擬模塊，將其作為一個開放式模塊，可加載到其它軟件系統(tǒng)中，提高正演模擬的精度和效率.

1 各向異性煤層中試射法射線追蹤

相對于一般VTI型介質來說，VTI型煤層的結構更為復雜，射線方向的計算精度要求更高.本文將VTI型煤層的試射法射線追蹤分為以下幾個步驟進行：給定群角進行射線試射，進行群相關系轉換，尋找射線出射角度，計算走時.

1.1 群相關系

各向同性介質中，相速度與群速度傳播方向相同，大小相等，而在各向異性介質中，波的相速度指波前傳播的速度，沿著相角方向傳播；群速度是指能量的傳播速度，沿著群角方向傳播，決定地震射線的分布.群、相速度關系式為

式中v（θ）為相速度，v（φ）為群速度.

1.2 VTI型煤層中的Snell定律

由惠更斯原理可知，介質中的反射與透射均遵循Snell定律.在VTI型煤層中，介質對稱軸一般不與界面法線方向一致，因此，其反射與透射關系并不能用簡單的Snell定律來描述.根據圖1并結合Snell定律，將VTI型煤層中的Snell定律表述為

圖1 VTI介質反射、透射示意圖Fig.1 Reflection and transmission in VTImedia

式中α1，α2，α3分別為 P波的入射角、反射角、透射角，v（θ1），v（θ2），v（θ3）分別為入射、反射、透射時的相速度.

1.3 射線角度和方向計算

1）初始角度

試射法射線追蹤首先需要給定射線的初始角度.因為群角是與射線走向相關的，射線追蹤時初始群角是可以設定的，給定初始群角φ的范圍為：arctan（xmin／2／z）～arctan（xmax／2／z），其中xmin，xmax分別為最小炮檢距和最大炮檢距，z為反射面的深度.

2）根據群角計算相角

根據上式作f（φ）的近似式，通過遞推得到相角θ.

3）反射相角、透射相角計算

由式（2）可知，反射、透射角度不能直接通過Snell定律的計算得出，因此，首先計算射線參數P，并給定一個射線收斂的精度范圍，然后通過多次修改初始射線參數進行多次試射，使出射點與接收點接近.本文采用非線性方程求解方法——對分法來求反射相角、透射相角.

根據Snell定律（式（2）），假設f（θ）＝P－sinα／v（θ）＝0，則對于非線性方程f（θ）＝0，設其搜索區(qū)間為［0，π／2］，從端點θ0＝0開始，以步長λ逐步向后進行搜索．對于在搜索區(qū)間范圍內的［θi，θi＋1］，分f（θi）＝0，f（θi）f（θi＋1） ＝0，f（θi）f（θi－1） ＞0，f（θi）f（θi－1） ＜0四種情形搜索方程的實根，從區(qū)間的左端點開始，直到區(qū)間的右端點結束.采用對分法求取方程的實根時，步長的選擇不應過大或過小，選擇過大的步長，會導致一些存在的實根的丟失，而選擇過小的步長，則在一定程度上增加計算量，因此以較大步長但不丟失實根為宜.

1.4 走時計算

對于給定的射線參數，先由Snell定律通過計算得出射線在第i層的相角θi及與相角相對應的相速度vi（θ），再求出與相角相對應的群速度vi（φ）以及群角φi，則可以得到射線經過n層介質的PP波的射線走時為

式中，hi為第i層的厚度，v（φ）i為φi對應的群速度.

2 模型試算

表1為模型參數，第2層和第4層為VTI型煤層，其中第2層的頂、底界面起伏，但厚度保持不變，第4層界面起伏，厚度變化.各層均表現(xiàn)各向異性.模型大小為1 000 m×700 m.震源位于地表.

表1 地質模型參數Table 1 Parameters in geologicalmodel

采用觀測系統(tǒng)：單邊放炮，每炮60道，道間距10 m，最小偏移距150 m.為了在圖形顯示時，射線能被清晰地分辨，對射線按5∶1的比例抽稀.通過對模型進行射線追蹤，所得射線路徑如圖2所示，圖中煤層的頂、底界面反射、透射清晰.如果將各層的各向異性參數ε，δ均設為0，則煤層為各向同性煤層，射線追蹤結果如圖3所示.將上述兩圖對比可知：各向同性煤層與VTI型煤層中的射線路徑有明顯區(qū)別；沒有出現(xiàn)射線盲區(qū).圖4將各向同性模型與各向異性模型追蹤的走時進行對比，可以看出：地震走時記錄中有四條同相軸，能有效分辨薄層；各向同性煤層與各向異性煤層的走時有較大差別，更進一步說明在煤層地震波的研究中，考慮煤層的各向異性給地震波傳播帶來的影響是必要的.

圖2 VTI型煤層中PP波射線路徑Fig.2 Ray path of PPwave in VTI coal

圖3 各向同性煤層中PP波射線路徑Fig.3 Ray path of PPwave in isotropicmedia

在INTEL P4 2.8 GHz的CPU和1GB的內存配置下，采用Fortran和VC＋＋混合編程自主開發(fā)VTI型煤層的射線追蹤軟件，并比較了4個反射界面各向異性和各向同性射線追蹤的計算效率.如圖5所示，VTI型煤層中的計算時間略高于各向同性介質的射線追蹤時間，反射界面多的模型射線追蹤計算時間略高于簡單模型（大約7.6% ～12%），但是總耗時不超過30 s，表明計算速度較快.

3 結語

分析了VTI型煤層的射線追蹤，采用VC＋＋和Fortran語言混合編程技術，實現(xiàn)了VTI型煤層的PP波射線追蹤，并自主研發(fā)出實用的地震記錄正演模擬軟件（已申請軟件著作權保護），只需改變速度和角度的相應參數，該方法還可適用于PSV轉換波的記錄模擬，實現(xiàn)多種波形的射線追蹤，并可將射線追蹤法推廣到HTI型煤層，以及更為復雜的且與實際煤層更相近的TTI型煤層中；通過模型試算，對各向同性煤層與VTI型煤層走時進行對比，論證了煤層中的各向異性對地震記錄有明顯影響，模型模擬結果與理論分析一致，為煤礦地震勘探中各向異性分析方法的驗證提供了重要依據.

圖4 各界面反射到地表的走時對比（實線表示煤層為VTI型，虛線表示煤層為各向同性）Fig.4 Travel time of reflection waves（Solid lines∶in VTImedia；Dotted lines∶in isotropic media）

圖5 計算時間的比較（無填充柱體表示煤層為各向同性，填充柱體表示煤層為各向異性）Fig.5 Comparison of calculation time（Column∶in isotropic coal；Filling column∶in VTI coal）

［1］David G.Seismic anisotropy in coal beds［J］.SEG Expanded Abstracts，2005：142－146.

［2］Arjun T，Shilpi G，Archana S，et al.Seismic azimuthal anisotropy：An important tool for detection of coal bed methane，an unconventional source of energy—A review［J］.Geohorizons，2010，12：23－27.

［3］Liu E R，Crampin S，Queen J，et al.Velocity and attenuation anisotropy caused by microcracks and macrofractures in multiazimuth reverse VSP［J］.Canadian Journal of Exploration Geophysics，1993，29（1）：177－188.

［4］Chen T J，Cui R F，Liu E R.AVO forward modeling for VTI coal［J］.Journal of China Coal Society，2009，34（4）：438－442.

［5］Wang J，Chen D H，Zhang H L，etal.Studies on phase and group velocities from acoustic logging［J］.Applied Geophysics，2012，9（1）：108－113.

［6］Wei Z T，F(xiàn)an Y R，Chen X L，et al.Application of shear wave anisotropy in fractures and in-situ stress analysis［J］. Progress in Geophys，2012，27（1）：217－224.

［7］Dong SH.Test on elastic anisotropic coefficients of gas coal［J］.Chinese JGeophys，2008，51930：947－952.

［8］Liu Y Y，Wang Y，Zhang MG.Discussion on“Teston elastic anisotropy coefficients ofgas coal”［J］.Progress in Geophys，2012，27（4）：1832－1836.

［9］Yang D Y，Peng S P，Wang Y，et al.Seismic record simulation for coal-seam containing fractures［J］.Coal Geology＆Exploration，2007，35（6）：61－66.

［10］Pei Z L，Dong Y S，Peng SP.Seismic record simulation azimuthally anisotropy for coal-seam containing fractures［J］.Oil Geophysical Prospecting，2007，42（6）：665－672.

［11］Li D H，Dong S H，Zhao X C，et al.Seismic wave simulation of anisotropy in coal-bed media［J］.Geophysical＆Geochemical Exploration，2010，34（5）：604－609.

［12］Li D H，Dong SH，Zhao X C，et al.Seismic wave simulation of bipbase EDA media in coal-bed media［J］.Progress in Geophys，2011，26（2）：654－663.

［13］Qian Jin，Cui Ruofei，Chen Tongjun.Anisotropic numerical simulation of coal-bearing strata with finite-difference［J］.Coal Geology＆Exploration，2010，38（2）：63－67.

［14］Sun ZQ，Sun JG，Han FX.Traveltime computation using fastmarchingmethod under complex topographical conditions［J］.Chinese Journal of Computational Physics，2010，27（2）：281－286.

［15］Chen S C，Ma Z T，Wu R S.Wave equation method for wavefield propagation path tracing［J］.Chinese Journal of Computational Physics，2008，25（2）：186－194.

［16］Peng SP，Gao Y F，Yang R S，et al.Theory and application of AVO for detection of coalbed methane-A case from the Huainan coalfield［J］.Chinese JGeophys，2005，48（6）：1475－1486.

［17］Shi Y，Wang Y，Lu J.Analysis of azimuthal anisotropy in coalmeasures by 3CVSP data［J］.Journal of China Coal Society，2007，32（8）：813－817.

［18］Gao EG，Xu GM，Jiang X Y，etal.Segment ray tracing in 3Dmedia［J］.Oil Geophysical Prospecting，2002，37（1）：11－16.

［19］Zhao A H，Ding Z F.New approximate expressions of seismic group velocities for weakly anisotropy media［J］.Progress in Geophys，2005，20（4）：916－919.

［20］Li L.Phase velocity and group velocity in TImedia［J］.Geophysical Prospecting for Petroleum，2008，47（4）：334－338.

［21］Li JG，Li Y P，Guo X L.Shooting ray tracing in VTImedia［J］.Oil Geophysical Prospecting，2010，45（1）：491－496.

［22］Li Q，Bai C Y.Review on seismic wavefront and raytracing in complex media［J］.Progress in Geophys，2012，27（1）：0092－0104.

［23］Thomsen L.Weak elastic anisotropy［J］.Geophysics，1986，51（10）：1954－1966.

［24］Li Juan，Zhang Qin.Mixed programming of Visual C＋＋ based on dynamic link library［J］.Computer Engineering and Design，2010，31（9）：2140－2144.

［25］Liu JX，JiW，Liu Y，et al.Simulation and visualization of 1D fCSEMbased on Matlab and Fortran mix-programming［J］. Computing Techniques for Geophysical and Geochemical Exploration，2012，34（2）：224－228.

Seism ic Response of Anisotropic Coal Based on Ray Theory

LIQin1，LIQingchun2，ZHANG Lin3
（1.Xi'an University ofScience and Technology，College ofGeological and Environmental Sciences，Xi'an 710054，China；2.Chang'an University，College ofGeology Engineering and Geomatics，Xi'an 710054，China；3.Shaanxi Center ofGeological Survey，Xi'an 710016，China）

We illustrate a technique of ray tracing by shooting method and analyze rules of ray in VTI coal.With mixed-language programming of VC＋＋and Fortran，ray tracing in VTI coal with a complex velocitymodel is realized.It demonstrates necessaries of anisotropic ray tracing in VTI coal by comparison of ray tracing in VTI coal and isotropic media.By contrast of calculation time，it indicates that thismethod is efficient tomeet needs of seismic records forward.The technique can be used inmulti-wave ray tracing in anisotropic coal as giving responding formulae and parameters.

VTI coal；ray tracing；anisotropy

date：2013－07－27；Revised date：2013－11－25

P631

A

1001-246X（2014）04-0444-05

2013－07－27；

2013－11－25

國家自然科學基金（41304105，41374145）；高等學校博士學科點專項科研基金（20120205130002）；陜西省教育廳專項科研計劃（12JK0470）資助項目

李勤（1979－），女，博士，講師，從事應用地球物理信息技術工作，E-mail：eriliqin＠126.com