閻守國，謝馥勵，龔丹，章成廣，張碧星

1中國科學院聲學研究所聲場聲信息國家重點實驗室，北京 100190

2長江大學油氣資源與勘探技術教育部重點實驗室，荊州 434102

1 引言

隨著世界油氣勘探的發(fā)展，裂縫性儲集層已成為油氣勘探的主要對象之一.在聲波測井中，由于地層中的裂縫寬度一般很小，裂縫與井壁之間的位置關系也不確定，且往往是多條裂縫一起構成裂縫帶，因此裂縫地層對陣列聲波測井波形的影響十分復雜，這使得對裂縫地層的識別與評價成為當今測井所面臨的難題之一.通過數(shù)值模擬計算，對含裂縫地層的井孔聲場進行研究，掌握裂縫對井內(nèi)聲場的影響規(guī)律，對于聲波測井數(shù)據(jù)處理和解釋具有重要意義.

針對裂縫性地層的井孔聲場問題，早期的研究工作采用平板狀裂縫模型（Hornby et al.，1989），利用低頻解析公式分析計算了斯通利波通過單一裂縫時的傳播及反射特性，之后類似的研究（Tang，1990；Tang and Cheng，1993）指出，聲波在經(jīng)過滲透帶處也會產(chǎn)生反射的斯通利波，后來也有學者通過不同的方法（Derov et al.，2009；Bakku et al.，2013）研究了裂縫存在時的井孔聲場，但這些方法大都假定聲場具有軸對稱性；積分方程法（Spring and Dudley，1992）、有限差分法（Kostek et al.，1998；陳德華等，2004；叢健生，2004）和有限元方法（Matuszyk，2013）陸續(xù)被應用于裂縫性地層井孔聲波場的研究中，但他們的計算都針對二維情況，并且水平層或裂縫寬度較大.由于有限元方法在地震波模擬問題上比有限差分法更耗時，因此對于井孔聲場的模擬大都采用有限差分方法，三維交錯網(wǎng)格有限差分法近年來被應用于大斜度井、傾斜分層地層及各向異性介質(zhì)地層等非軸對稱問題的數(shù)值模擬研究（Leslie and Randall，1992；Cheng et al.，1995；Sinha，2006；林偉軍等，2006；閻守國等，2011），得到了較好的效果，但采用該方法對非軸對稱裂縫地層的井孔聲場研究還未見報道.這主要是由于當裂縫寬度較小時，需要劃分細小的網(wǎng)格來滿足薄裂縫的計算需求，而這將大大增加三維有限差分的計算量，導致模擬計算無法順利進行.本文針對這一問題，對含有傾斜薄裂縫孔隙地層中的井孔聲波場采用三維不規(guī)則交錯網(wǎng)格有限差分法進行數(shù)值模擬計算.考察不同裂縫及地層參數(shù)情況下的井孔聲場特性，為進一步研究利用井孔聲場信息反演裂縫參數(shù)奠定基礎.

2 物理模型及計算方法

2.1 物理模型

本文采用三維應力－速度交錯網(wǎng)格有限差分法數(shù)值模擬計算含有傾斜裂縫地層的井孔聲場問題.圖1為含有傾斜裂縫地層的井孔模型示意圖.假定聲源位于笛卡爾坐標系的原點，井軸與z軸重合，井外裂縫與xz軸構成的平面垂直，與xy軸構成的平面傾角為α，裂縫下界面與z軸的交點與原點的距離為dz，裂縫的垂直厚度為H，裂縫間隔為h.

圖1 傾斜裂縫地層中的井孔模型Fig.1 The configuration of borehole in porous formation with tilted thin fracture

2.2 變網(wǎng)格有限差分處理

由于需要計算的裂縫很薄，并受井孔模型大小的影響，如果采用很小的網(wǎng)格進行數(shù)值計算會導致計算速度過慢，甚至出現(xiàn)無法計算的情況.所以我們采用不均勻網(wǎng)格的辦法在裂縫處采用小網(wǎng)格，在非裂縫處采用大網(wǎng)格，這樣可以大大的減小計算量，在保證計算精度的同時提高計算速度.在使用變網(wǎng)格時，網(wǎng)格步長的變化可能會導致數(shù)值反射的出現(xiàn).其原因是因為波場離散后，相速度是網(wǎng)格步長的函數(shù)，當相速度梯度較大時，即使速度和密度都沒有變化，入射波的能量也會部分反射回來，導致數(shù)值反射現(xiàn)象.在網(wǎng)格變化的區(qū)域?qū)Σ▓鲞M行插值計算雖可以在某種程度上壓制數(shù)值反射，但是插值算法計算量大，且效果并不令人滿意，所以本文沒有采用插值計算.本文采用的算法是根據(jù)各點所對應的步長的變化計算其差分系數(shù)，此方法不會產(chǎn)生數(shù)值反射，并且減小了計算量（孫衛(wèi)濤和楊慧珠，2004）.下面給出了具有四階空間精度的不規(guī)則網(wǎng)格差分算子（圖2）.

圖2 變網(wǎng)格差分示意圖Fig.2 The configuration of nonuniform finite－difference grid

節(jié)點i處的場量ψ對空間的一階導數(shù)?ψ／?x可以表示為相鄰節(jié)點的線性組合形式：

其中，

（1a）即為網(wǎng)格發(fā)生變化時的四階空間精度差分算子.而當s1＝r1＝r2＝Δr時，即在均勻網(wǎng)格情況下，（1a）式可以退化為均勻網(wǎng)格四階空間精度差分算子，此時

2.3 流體－孔隙介質(zhì)界面處的差分方程統(tǒng)一

由于實際地層中往往是裂縫與孔隙同時存在，為了能夠更加接近實際情況，同時也為了考察孔隙度滲透率等地層參數(shù)對裂縫產(chǎn)生的影響，本文在數(shù)值模擬時采用流體飽和孔隙介質(zhì)作為井外地層介質(zhì)，地層中的彈性波傳播采用Biot孔隙彈性波方程組描述（Biot，1955；1956；1962a；1962b；王秀明等，2003）.孔隙介質(zhì)中場量和參數(shù)更多，波動方程形式復雜，與流體波動方程形式不統(tǒng)一，不能直接應用參數(shù)平均法獲得流體－孔隙介質(zhì)交界面上場量的離散表達式，因此無法直接進行差分計算.關威等（Guan et al.，2009）曾利用井壁邊界上的應力和位移連續(xù)性條件推導出了井壁網(wǎng)格點的差分離散表達式，實現(xiàn)了對孔隙地層中聲波測井的有限差分模擬算法，但這種處理方法在算法的實現(xiàn)上相對麻煩，尤其是地層中存在裂縫時，需要利用界面處的邊界條件重新推導差分公式.本文采用將裂縫中的流體利用孔隙介質(zhì)方程參數(shù)取流體極限的辦法，統(tǒng)一了彈性波動方程，簡化了介質(zhì)交界面處的處理過程.

依據(jù)Biot理論，孔隙介質(zhì)的本構關系為

其中，Tij為總應力張量，S為流體應力，Pf為流體壓強，分別為顆粒、骨架以及流體的體積模量，u為固體位移矢量，w為固體相對流體位移矢量，φ為孔隙度，μb是骨架剪切模量，λc是拉梅系數(shù)，λc＝H－2μb，H＝λb＋2μb＋α2M，M＝

將方程（2）兩邊對時間t求導，并進行分量展開，可寫為

其中，Vx，Vy，Vz分別表示固體介質(zhì)在x，y，z方向上的質(zhì)點速度分量，Wx，Wy，Wz分別表示固體相對流體介質(zhì)在x，y，z方向上的質(zhì)點速度分量.

孔隙介質(zhì)的運動方程為

式中，σij＝Tij＋φPfδij，ρ11＝ （1－φ）ρs－ （1－E）φρf，ρ12＝ （1－E）φρf，ρ22＝Eφρf，ρs，ρf分別為固體顆粒與孔隙流體密度，E為孔隙迂曲度，可近似表示為，K為滲透率，η為黏滯系數(shù).

如果令：Q1＝－ρ11ρ22，R1＝ρ12／Q1，R2＝ρ22／Q1，R3＝ （ρ12＋ρ22）／Q1，P1＝ （ρ11＋ρ12）／Q1，P2＝ （ρ12＋ρ22）／Q1，P3＝ （ρ11＋2ρ12＋ρ22）／Q1，則孔隙介質(zhì)中的速度分量可以寫為

當場量和介質(zhì)參數(shù)取極限的狀態(tài)下，描述孔隙介質(zhì)中彈性波傳播的Biot方程組可以退化為單相的固體或流體介質(zhì)中的彈性波方程組.本文僅考慮退化為流體的情況，當孔隙介質(zhì)參數(shù)選取如下極限狀態(tài)時，退化為流體介質(zhì)：φ＝1，Kb＝0，Ks＝Kf，H＝M＝Kf，μ＝0，a＝1，b＝0，由于此時ρ11＝0，ρ12＝0，ρ22＝ρf，使得Q＝0，P，R的計算不再有意義，可根據(jù)流體中的波動方程直接定義R1＝1／ρf，R2＝0，R3＝0，流體中不存在相對滲流位移，可直接令W＝0.在流體與孔隙介質(zhì)交界面處進行參數(shù)平均處理，密度取算數(shù)平均值，彈性系數(shù)取調(diào)和平均值.經(jīng)過以上的處理，在計算中遇到流體和孔隙介質(zhì)交界面時可以直接進行差分計算，非常適合處理含有裂縫的地層模型.

2.4 差分算法的實現(xiàn)

本文采用三維應力－速度交錯網(wǎng)格有限差分，各物理量在網(wǎng)格上的離散規(guī)則如下：

其中，下標i，j，k表示為整網(wǎng)格點，下標i＋1／2，j＋1／2，k＋1／2為半網(wǎng)格點，聲源類型為點源，時域脈沖選擇變形的Ricker子波，表達式為f（t）＝，時域波形及頻譜如圖3所示.

圖3 聲源時域波形及頻譜Fig.3 The waveform in time domain and spectrum in frequency domain of source pulse

3 計算結果及分析

為符合裂縫性儲層的實際特點，本文選擇低孔隙度、低滲透率地層進行數(shù)值模擬計算.所采用的介質(zhì)參數(shù)如表1所示，其中孔隙介質(zhì)參數(shù)Ks，ρs為孔隙固體顆粒的體積模量及密度，VrP，VrS為骨架縱、橫波速度，η，φ，K分別為黏滯系數(shù)、孔隙度及滲透率，Vf，ρf為孔隙內(nèi)流體的速度及密度.井孔半徑a＝0.1m，計算過程中當聲源頻率不大于5kHz時，空間網(wǎng)格最大可選為2cm，變網(wǎng)格差分中可任意更改小網(wǎng)格大小，根據(jù)裂縫寬度不同，本文計算過程中最小網(wǎng)格選取為20μm.

表1 模型的介質(zhì)參數(shù)Table 1 The parameters of the media

3.1 算法正確性驗證

當井外為無裂縫的孔隙介質(zhì)地層時，將有限差分結果與實軸積分法結果進行對比，對比結果如圖4所示.聲源頻率為2.5kHz，接收器與聲源間距為1m.圖中黑線為有限差分法計算結果，紅線為實軸積分法計算結果.兩種方法計算結果基本吻合，證明了有限差分方法的正確性.圖中幅度上的差異是由歸一化造成的，不影響本文結果.

圖4 有限差分與實軸積分法的結果對比圖Fig.4 The comparison results of real axis integrate method and finite－difference method

采用完全匹配層（PML）來處理各方向吸收邊界（張魯新等，2010）.圖5為井外為均勻孔隙介質(zhì)情況下，xz平面上z方向應力（Tzz）的聲場快照，圖中所示時刻為1.35ms，x，y，z方向上的計算距離均為±2.0m，從圖5可以看出各方向均無反射波存在，邊界吸收效果良好.

圖5 井外無裂縫時xz平面上Tzz的聲場快照圖Fig.5 The Tzzacoustical fields snapshot of xz plane when there is no fracture outside borehole

3.2 裂縫對井內(nèi)聲場的影響

裂縫引起的介質(zhì)突變會使聲波在界面處產(chǎn)生反射，并使穿過裂縫的聲波能量減弱.圖6是井外存在裂縫情況下xz平面上場量Tzz的快照圖，計算范圍為±2.0m，水平裂縫位于聲源上方0.9m處，裂縫寬度為2mm.圖中P，S，ST分別代表縱波、橫波和斯通利波，下標d代表直達波，r代表由裂縫產(chǎn)生的反射波.由圖6a可見，此時陣列波形圖中可看到直達的縱波、橫波及斯通利波，斯通利波只在井孔內(nèi)沿井軸方向傳播.縱波通過裂縫后直達場減弱，此時從圖中還觀察不到反射波；聲波經(jīng)過一段時間的傳播，在圖6b中可以看到反射縱波，并且由于孔隙介質(zhì)的耗散作用，井內(nèi)的斯通利波在傳播過程中強度逐漸減弱；圖6c中可以看到反射橫波的出現(xiàn)，此時由于各波列之間的相互干擾，以及頭波的影響，井外聲場顯得比較復雜.雖然從陣列波形圖中可以明顯地看到反射波的存在，裂縫的位置也可以很直觀的反映出來，但是由于實際測井中只能記錄井內(nèi)的聲場，因此需要針對陣列波形圖進行研究.

圖6 不同時刻，井外存在裂縫時xz平面上的聲場快照圖（Tzz）（a）0.51ms；（b）0.68ms；（c）0.85ms.Fig.6 The Tzzacoustical fields snapshot of xz plane when there are fractures outside borehole

圖7為井外存在裂縫時井軸上的陣列波形圖，橫軸是接收時間，縱軸為接收點距離聲源的距離，波列的幅度經(jīng)過歸一化處理.從圖7中可以看到，井軸上的陣列波形圖存在縱波（P）、橫波（S）以及斯通利波（ST），由于聲源激發(fā)頻率及地層參數(shù)的影響，此時縱、橫波激發(fā)幅度較小，斯通利波占主要成分，當正常傳播的各種波型遇到裂縫時，就會產(chǎn)生相應的反射波，圖7中可以觀察到的反射波列有兩條，按照其傳播速度判斷分別為橫波及斯通利波反射.

圖7 井軸上的陣列波形圖Fig.7 The array waveform neceived on the borehole axis

如果能夠在陣列波形圖中觀察到反射波列，即可發(fā)現(xiàn)裂縫的存在.但是由于全波列聲波測井儀器的聲源頻率一般在15kHz以內(nèi)，比較常用的頻率在2.5～10kHz之間，而地層中裂縫的寬度一般都比較小，介于20～2000μm之間，裂縫產(chǎn)生的反射波強度往往很弱，因此在聲波測井的頻段范圍內(nèi)能否觀察到細小裂縫的存在是我們首先要考察的問題.

圖8是數(shù)值模擬計算得到的聲源頻率為2.5kHz時不同裂縫寬度情況下，井軸上接收到的陣列全波圖.此時全波圖的主要成分是斯通利波，縱、橫波的激發(fā)幅度很小，從波列圖中可以得到斯通利波的速度約為1459m·s－1.圖中a、b分別是裂縫寬度為2cm和2mm的情況，從圖8a中可以觀察到裂縫反射的斯通利波，但隨著裂縫寬度減小，在圖8b中已經(jīng)無法觀察到反射波的存在.可以看出低頻聲源對單條裂縫的檢測能力較弱.

圖8 數(shù)值模擬計算得到的聲源頻率為2.5kHz時不同裂縫寬度情況下，井軸上接收到的陣列全波圖（a）H＝2cm；（b）H＝2mm.Fig.8 The simulation results of full waveform array received on the borehole axis with different fracture width，the source frequency is 2.5kHz

圖9 數(shù)值模擬計算得到的聲源頻率為5.0kHz時不同裂縫寬度情況下，井軸上接收到的陣列全波圖（a）H＝2cm；（b）H＝2mm；（c）H＝0.2mm；（d）H＝0.02mm.Fig.9 The simulation results of full waveform array received on the borehole axis with different fracture width，the source frequency is 5.0kHz

圖9是數(shù)值模擬計算得到的聲源頻率為5.0kHz時不同裂縫寬度情況下，井軸上接收到的陣列全波圖.與2.5kHz時的情況相同，此時波列中占主要成分的依然是斯通利波，其速度約為1490m·s－1.從圖中可以看到，當裂縫寬度較大時（圖9a），存在明顯的兩條反射波列，按照其傳播速度判斷分別為地層中沿井壁滑行的橫波及井內(nèi)流體中斯通利波，當裂縫寬度減小后，反射波列中的斯通利波明顯減弱，在裂縫寬度為2mm時（圖9b）就已經(jīng)很難觀察到反射斯通利波波列，但反射波列中的橫波卻并沒有減弱，并且當裂縫寬度減小到0.02mm時依然可以觀察到明顯的反射橫波.這一現(xiàn)象有助于我們探測水平薄裂縫的存在.

實際儲層中的裂縫往往不是單一出現(xiàn)的，而是由多條薄裂縫一起構成一條裂縫帶，我們定義裂縫帶密度為d＝H／（H＋h），它相當于裂縫在整個裂縫帶范圍內(nèi)的占空比.圖10給出了由單條寬度為2mm的裂縫構成的裂縫帶在不同的裂縫帶密度情況下的陣列波形圖，圖中裂縫帶寬度均為6cm，圖10（a，b）是頻率為2.5kHz，裂縫帶密度分別為0.1和0.25時的波形圖，圖10（c，d）是頻率為5.0kHz時的對應情況.圖8b及圖9b是單條2mm的裂縫在這兩個頻率下的陣列波形圖，通過對比可以看出聲波通過裂縫帶時會產(chǎn)生相比于單條裂縫更強的反射斯通利波，且隨著裂縫帶密度的增加反射波也隨之增強，最終會形成相當于同等寬度的寬裂縫的反射效果.同時，對比2.5kHz與5.0kHz時的結果，發(fā)現(xiàn)高頻聲源在裂縫帶密度相對較低時的反射波列更加明顯.因此，在聲波測井常用的頻率范圍內(nèi)，高頻聲源對低密度裂縫帶的檢測能力更強，更適合裂縫帶的檢測.

在實際地層中裂縫與井孔的關系往往很復雜，可能出現(xiàn)與井孔相交或者不相交等不同情況，經(jīng)過模擬計算發(fā)現(xiàn)，在本文所使用的聲源類型及頻率范圍內(nèi)，未與井孔相交的裂縫或裂縫帶對井內(nèi)聲場產(chǎn)生的影響很難通過陣列波形圖加以辨別，因此本文只考慮裂縫或裂縫帶與井孔相交存在一定傾斜角度的情況.圖11是聲源頻率為5.0kHz時，裂縫帶相對井孔的傾斜角度分別為15°、30°和45°時的陣列波形圖，裂縫帶由單條寬度為2mm的裂縫構成，裂縫帶密度為0.25，總寬度為6cm.通過與相同狀態(tài)下水平裂縫的對比（圖10d），可以發(fā)現(xiàn)此時裂縫帶的反射橫波消失了，即使傾斜角度較小時（圖11a）也無法觀察到；但裂縫傾斜角度對反射斯通利波的影響并不大，圖中用符號ST標注出了裂縫產(chǎn)生的反射斯通利波波列，對比圖11（a—c）可發(fā)現(xiàn)，即使在裂縫傾斜角度較大時依然可以觀察到明顯的反射斯通利波波列.

圖10 單條寬度為2mm的裂縫構成的裂縫帶在不同的裂縫帶密度情況下的陣列波形圖，裂縫帶寬度均為6cm（a）f＝2.5kHz，d＝0.1；（b）f＝2.5kHz，d＝0.25；（c）f＝5.0kHz，d＝0.1；（d）f＝5.0kHz，d＝0.25.Fig.10 The simulation results of full waveform array with different fracture zone density，the fracture zone constituted by the fractures 2mm width，and the fracture zone width is always 6cm.

圖11 不同裂縫傾斜角度時的陣列波形圖（5.0kHz）（a）α＝15°；（b）α＝30°；（c）α＝45°.Fig.11 The simulation results of full waveform array with different tilted angle of the fracture zone and borehole axis，the source frequency is 5.0kHz.

圖12 不同滲透率條件下的陣列波形圖（a）K＝0.002μm2；（b）K＝0.005μm2.Fig.12 The simulation results of full waveform array with different permeability of formation

圖12顯示了孔隙介質(zhì)滲透率改變時陣列波形圖的變化情況，聲源激發(fā)頻率為5kHz，水平裂縫，圖中標出了各波列的位置，P，S，ST分別對應縱波、橫波及斯通利波.從圖12可以看到，從聲源發(fā)出的聲波在未遇到裂縫時有縱波、橫波及斯通利波三種成分，由于此時縱、橫波激發(fā)強度較弱，很難從陣列波形圖中觀察到它們的反射波，而斯通利波幅度相對較大，當滲透率相對較小時，斯通利波衰減較弱，當其受到界面的影響時，產(chǎn)生較強的反射橫波及斯通利波；隨著滲透率的增大，斯通利波的衰減明顯增強，使得斯通利波在到達裂縫時能量減弱，此時反射波還可以觀察到反射的橫波但反射的斯通利波則明顯減弱，很難從波形圖中觀察到.

4 結語

本文采用有限差分法數(shù)值模擬了有裂縫及裂縫帶存在時孔隙介質(zhì)地層中的井孔聲場問題，采用非均勻網(wǎng)格有限差分，所計算的最低裂縫寬度達到了20μm，利用統(tǒng)一形式的聲波波動方程，通過將孔隙介質(zhì)退化為流體的辦法，簡化了流體與孔隙介質(zhì)交界面處的處理過程.通過對裂縫及裂縫帶寬度、裂縫傾斜角度及孔隙介質(zhì)滲透率等因素的考察總結出以下結論：

（1）當存在水平裂縫時，使用低頻聲源（2.5kHz）只能觀察到反射斯通利波，而使用高頻聲源（5kHz以上）還會觀察到反射橫波波列，隨著裂縫寬度的減小，反射斯通利波逐漸減弱，但反射橫波強度變化不大，可以用來檢測單條水平薄裂縫的存在.

（2）當多條裂縫構成裂縫帶時，裂縫帶的密度和寬度會對反射波的強度產(chǎn)生較大影響，當密度較小時，裂縫帶無法形成整體作用，反射波較弱，但隨著裂縫帶密度的增大，反射波逐漸增強，最終會形成相當于同等寬度的寬裂縫的反射效果.

（3）與井孔相交的裂縫隨著傾斜角度的增加，其反射橫波會消失，但對斯通利波的影響不大，傾斜角度很大時依然存在反射斯通利波，這對傾斜裂縫的檢測是有利的.

（4）斯通利波在傳播過程中會向孔隙介質(zhì)輻射能量，產(chǎn)生衰減，隨著滲透率的增大，斯通利波的衰減也隨之增大，會導致界面處斯通利波幅度減小，可能觀察不到界面的反射斯通利波.

本文通過數(shù)值模擬的辦法，定性研究了孔隙地層井孔中聲波通過裂縫及裂縫帶時的反射規(guī)律，更加深入的工作是找到聲波的衰減及反射透射規(guī)律與裂縫參數(shù)之間更深層次的定量關系，最終實現(xiàn)利用陣列聲波信息對地層中裂縫相關參數(shù)的反演，這是我們下一步將要進行的主要工作.

Bakku S K，F(xiàn)ehler M，Burns D.2013.Fracture compliance estimation using borehole tube waves.Geophysics，78（4）：D249－D260.

Biot M A.1955.Theory of elasticity and consolidation for a porous anisotropic solid.J.Apply.Physica，26（2），doi：10.1063／1.1721956.

Biot M A.1956.Theory of deformation of a porous viscoelastic anisotropic solid.JournalofAppliedPhysics，27（5），doi：10.1063／1.1722402.

Biot M A.1962a.Mechanics of deformation and acoustic propagation in porous media.J.Apply.Physica，33（4），doi：10.1063／1.1728759.

Biot M A.1962b.Generalized theory of acoustic propagation in porous dissipative media.TheJournaloftheAcoustical SocietyofAmerica，34（9）：1254－1264.

Chen D H，Cong J S，Xu D L，et al.2004.The simulation of borehole field in fractural formation.JournalofDaqing PetroleumInstitute，28（3）：4－6，13，116.

Cheng N Y，Cheng C H，Toks?z M N.1995.Borehole wave propagation in three dimensions.J.Acoust.Soc.Am，97（6）：3483－3493.

Cong J S.2004.Simulated for the acoustic field of the layer and fractures in and outside of a borehole using finite difference method［Master thesis］.Daqing：Daqing Petroleum Institute.Derov A，Maximov G，Lazarkov M，et al.2009.Characterizing hydraulic fractures using slow waves in the fracture and tube waves in the borehole.2009SEG Annual Meeting，4115－4119

Guan W，Hu H S，He X.2009.Finite－difference modeling of the monopole acoustic logs in a horizontally stratified porous formation.J.Acoust.Soc.Am，125（4）：1942－1950.

Hornby B E，Johnson D L，Winkler K W，et al.1989.Fracture evaluation using reflected Stoneley－wave arrivals.Geophysics，54（10）：1274－1288.

Kostek S，Johnson D L，Randall C J.1998.The interaction of tube waves with borehole fractures.PartⅠ：Numerical models.Geophysics，63（3）：800－808.

Leslie H D，Randall C J.1992.Multipole sources in boreholes penetrating anisotropic formations：Numerical and experimental results.J.Acoust.Soc.Am.，91（1）：12－17.

Lin W J，Wang X M，Zhang H L.2006.Acoustic wave propagation in a borehole penetrating an inclined layered formation.Chinese J.Geophys.（in Chinese），49（1）：284－294，doi：10.3321／j.issn：0001－5733.2006.01.036.

Matuszyk P J，Torres－Verdín C，Pardo D.2013.Frequency－domain finite－element simulations of 2Dsonic wireline borehole measurements acquired in fractured and thinly bedded formations.Geophysics，78（4）：D193－D207.

Sinha B K，?im?ek E，Liu Q H.2006.Elastic－wave propagation in deviated wells in anisotropic formations.Geophysics，71（6）：D191－D202.

Spring C T，Dudley D G.1992.Acoustic－wave propagation in a cylindrical borehole with fractures.J.Acoust.Soc.Am.，91（2）：658－669.

Sun W T，Yang H Z.2004.A3－D finite difference method using irregular grids for elastic wave propagation in anisotropic media.ChineseJ.Geophys.（in Chinese），47（2）：332－337，doi：10.3321／j.issn：0001－5733.2004.02.022.

Tang X M.1990.Acoustic logging in fractured and porous formations［Ph.D.thesis］.Cambridge：Massachusetts Institute of Technology.

Tang X M，Cheng C H.1993.Borehole stoneley wave propagation across permeable structures 1.GeophysicalProspecting，41（2）：165－187.

Wang X M，Zhang H L，Wang D.2003.Modelling of seismic wave propagation in heterogeneous poroelastic media using a highorder staggered finite－difference method.ChineseJ.Geophys.（in Chinese），46（6）：842－849.

Yan S G，Song R L，LüW G，et al.2011.Numerical simulation of acoustic fields excited by cross－dipole source in deviated wells in transversely isotropic formation.ChineseJ.Geophys.（in Chinese），54（9）：2412－2418，doi：10.3969／j.issn.0001－5733.2011.09.025.

Zhang L X，F(xiàn)u L Y，Pei Z L.2010.Finite difference modeling of Biot′s Poroelastic equations with unsplit convolutional PML and rotated staggered grid.ChineseJ.Geophys.（in Chinese），53（10）：2470－2483，doi：10.3969／j.issn.0001－5733.2010.10.021.

附中文參考文獻

陳德華，叢健生，徐德龍等.2004.裂縫性地層中的井孔聲場模擬.大慶石油學院學報，28（3）：4－6，13，116.

叢健生.2004.利用有限差分法模擬計算具有分層和裂縫地層井內(nèi)外聲場［碩士論文］.大慶：大慶石油學院.

林偉軍，王秀明，張海瀾.2006.傾斜地層中的井孔聲場研究.地球物理學報，49（1）：284－294，doi：10.3321／j.issn：0001－5733.2006.01.036.

孫衛(wèi)濤，楊慧珠.2004.各向異性介質(zhì)彈性波傳播的三維不規(guī)則網(wǎng)格有限差分方法.地球物理學報，47（2）：332－337，doi：10.3321／j.issn：0001－5733.2004.02.022.

王秀明，張海瀾，王東.2003.利用高階交錯網(wǎng)格有限差分法模擬地震波在非均勻孔隙介質(zhì)中的傳播.地球物理學報，46（6）：842－849.

閻守國，宋若龍，呂偉國等.2011.橫向各向同性地層斜井中正交偶極子激發(fā)聲場的數(shù)值模擬.地球物理學報，54（9）：2412－2418，doi：10.3969／j.issn.0001－5733.2011.09.025.

張魯新，符力耘，裴正林.2010.不分裂卷積完全匹配層與旋轉(zhuǎn)交錯網(wǎng)格有限差分在孔隙彈性介質(zhì)模擬中的應用.地球物理學報，53（10）：2470－2483，doi：10.3969／j.issn.0001－5733.2010.10.021.