王瑞甲,喬文孝*,鞠曉東
1 中國(guó)石油大學(xué)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 102249
2 北京市地球探測(cè)與信息技術(shù)重點(diǎn)實(shí)驗(yàn)室,北京 102249
隨鉆聲波測(cè)井在節(jié)省井架占用時(shí)間、利用測(cè)得的聲波速度模型與地震勘探數(shù)據(jù)相結(jié)合實(shí)時(shí)確定地層界面的位置、估算地層孔隙壓力等方面有著電纜測(cè)井無(wú)法比擬的優(yōu)勢(shì)[1].關(guān)于隨鉆聲波測(cè)井的研究,國(guó)內(nèi)外已做了大量的工作.Minear和Legget成功實(shí)現(xiàn)了地層隨鉆縱波測(cè)量[2-3];Tang等認(rèn)為采用四極子聲源進(jìn)行隨鉆橫波測(cè)量有著偶極子聲波測(cè)井無(wú)法比擬的優(yōu)勢(shì)[4];Sinha等也研究了隨鉆測(cè)井模型下各向同性地層井孔內(nèi)導(dǎo)波的基本響應(yīng)特征[5].目前,隨鉆聲波測(cè)井儀已基本實(shí)現(xiàn)了地層縱、橫波測(cè)量的功能,下一步所面臨的挑戰(zhàn)是對(duì)地層的聲學(xué)各向異性進(jìn)行測(cè)量.
各向異性測(cè)量主要包括快橫波面方位的確定和橫波各向異性值的測(cè)量?jī)蓚€(gè)方面.盡管四極子聲源在隨鉆地層橫波測(cè)量方面取得了成功,但是限于其方位特性,很難利用四極子聲源實(shí)現(xiàn)隨鉆地層各向異性測(cè)量.雖然部分學(xué)者已經(jīng)在此方面開(kāi)展了一些工作[6],但是至今未見(jiàn)成功利用四極子聲源實(shí)現(xiàn)地層各向異性測(cè)量的報(bào)道.采用正交偶極子聲波測(cè)井方式評(píng)價(jià)地層各向異性的方法已經(jīng)在電纜測(cè)井中得到了廣泛的應(yīng)用[7].因?yàn)殡S鉆四極子聲波測(cè)井儀換能器的安裝位置同正交偶極子聲波測(cè)井儀器相同,且其接收站兼具備正交偶極子接收功能,通過(guò)合理的電路設(shè)計(jì),可以較為方便地實(shí)現(xiàn)隨鉆正交偶極子聲波測(cè)井,所以采用正交偶極子聲源進(jìn)行地層各向異性評(píng)價(jià)的方法為隨鉆地層各向異性測(cè)量的首選方式.研究各向異性地層隨鉆正交偶極子聲波測(cè)井的響應(yīng)特征,對(duì)偶極子聲源在含鉆鋌各向異性地層井孔內(nèi)激發(fā)的聲場(chǎng)進(jìn)行分析,可以幫助理解在隨鉆條件下各向異性地層井孔內(nèi)沿井軸方向傳播的彎曲波的特征,為新一代隨鉆聲波各向異性測(cè)量?jī)x器的設(shè)計(jì)及測(cè)量方案的設(shè)計(jì)提供理論指導(dǎo).
有關(guān)各向異性地層井孔聲場(chǎng)的研究,國(guó)內(nèi)外已經(jīng)做了大量的工作.Cheng采用三維直角坐標(biāo)系有限差分方法模擬了正交各向異性地層包圍的井孔內(nèi)多極子聲源激發(fā)的聲場(chǎng)[8].Schmitt研究了介質(zhì)對(duì)稱軸同井軸平行情況下,多極子聲源激發(fā)的模式波的頻散曲線及衰減曲線,并分析了各地層參數(shù)對(duì)于井孔內(nèi)導(dǎo)波的影響[9].Sinha采用三維柱坐標(biāo)系有限差分方法模擬研究了典型的硬地層和軟地層條件下,TI地層斜井情況下多極子聲源激發(fā)的聲場(chǎng)以及儀器的存在對(duì)于井內(nèi)模式波頻散特征的影響,他認(rèn)為,各向異性地層中,彎曲波在低頻下的傳播速度為對(duì)應(yīng)地層橫波的相速度[10].王秀明采用三維直角坐標(biāo)系有限差分方法計(jì)算了TI地層斜井中的單極子聲源和偶極子聲源激發(fā)的聲場(chǎng),他的模擬結(jié)果表明聲波測(cè)井所測(cè)得的彎曲波的速度同各向異性地層體波的群速度一致[11].張碧星分別采用實(shí)軸積分和攝動(dòng)積分的方法研究了TI地層中模式波的頻散特性和激發(fā)強(qiáng)度[11].陳雪蓮和王瑞甲采用實(shí)軸積分法模擬了徑向分層TI孔隙介質(zhì)井孔內(nèi)多極子聲源激發(fā)的聲場(chǎng),并著重研究了滲透率對(duì)模式波衰減和幅度的影響以及井孔模式波的探測(cè)深度問(wèn)題[12-13].He和Hu等從理論上推導(dǎo)了井孔彎曲波的低頻極限速度公式,并采用三維柱坐標(biāo)系有限差分算法模擬了TI介質(zhì)斜井中的彎曲波,他們的研究結(jié)果表明,大多數(shù)情況下,快、慢彎曲波的慢度近似等于沿井軸方向傳播的地層快、慢橫波的慢度[14-15].閆守國(guó)和宋若龍等也模擬了橫向各向同性斜井中偶極子聲源激發(fā)的聲場(chǎng),并提出了采用守恒積分的方法解決柱坐標(biāo)系波動(dòng)方程在井軸上出現(xiàn)的奇異點(diǎn)的問(wèn)題[16].上述的研究均為電纜測(cè)井情況下各向異性地層聲波測(cè)井模擬,鮮見(jiàn)有關(guān)在隨鉆條件下各向異性地層偶極子聲源激發(fā)聲場(chǎng)研究的報(bào)道.
即使在地層為各向同性的情況下,由于鉆鋌占據(jù)了井內(nèi)的大部分空間,隨鉆條件下的彎曲波的頻散特性、激發(fā)特征均與電纜測(cè)井不同[1,17].在地層為各向異性的情況下,偶極子聲源激發(fā)的聲場(chǎng)將更為復(fù)雜,無(wú)法采用解析的方法進(jìn)行模擬.地層介質(zhì)最為廣泛存在的一種各向異性介質(zhì)模型為橫向各向同性(TI)介質(zhì).本文采用三維有限差分方法模擬研究了橫向各向同性(TI)地層隨鉆正交偶極子聲波測(cè)井,對(duì)地層的聲學(xué)各向異性在隨鉆正交偶極子聲波測(cè)井中的響應(yīng)特征進(jìn)行了分析,對(duì)采用隨鉆正交偶極子聲波測(cè)井方式進(jìn)行各向異性測(cè)量的可行性進(jìn)行了評(píng)價(jià).
與電纜測(cè)井不同,在隨鉆聲波測(cè)井中,鉆挺占據(jù)了井孔內(nèi)的大部分空間.由于鉆鋌的存在,井孔內(nèi)沿井軸方向傳播的各種模式波的性質(zhì)同電纜測(cè)井不同.圖1a為T(mén)I地層隨鉆測(cè)井聲學(xué)模型示意圖,S方向?yàn)門(mén)I介質(zhì)的對(duì)稱軸方向,它與井軸的夾角為α.如圖1b所示,TI地層隨鉆測(cè)井聲學(xué)模型可以簡(jiǎn)化為柱狀徑向分層聲學(xué)模型,沿井徑方向從內(nèi)向外的介質(zhì)依次為水、鋼(鉆鋌)、水、地層,各介質(zhì)的外徑分別為r0、r1、r2和無(wú)窮大.井孔內(nèi)充滿流體.井孔外地層為無(wú)限大TI介質(zhì).鉆鋌位于井孔中央,鉆鋌中間的水眼中充滿水.在實(shí)際測(cè)井中,尤其是在鉆進(jìn)過(guò)程中,鉆鋌并非完全居中,此時(shí)井中的聲場(chǎng)將更為復(fù)雜.為了突出本文所關(guān)心的問(wèn)題,本文的模型假設(shè)鉆鋌在井孔中完全居中.
由于地層為各向異性介質(zhì),該問(wèn)題不存在解析解,必須采用數(shù)值方法來(lái)模擬地層中的聲傳播.三維有限差分方法是模擬復(fù)雜介質(zhì)中聲傳播問(wèn)題的常用方法[6,8,10-11,15-16,18].本文采用三維有限差分方法來(lái)模擬隨鉆情況下的各向異性地層井孔中的聲傳播.任意各向異性介質(zhì)中的運(yùn)動(dòng)方程和本構(gòu)方程分別為式(1—3)和式(4):
圖1 隨鉆測(cè)井聲學(xué)模型示意圖,包括(a)TI地層井孔隨鉆測(cè)井聲學(xué)模型和(b)井孔橫截面示意圖Fig.1 Schematic of LWD acoustic model,including(a)acoustic model of borehole surrounded by TI formation in LWD conditions and(b)the cross section of the borehole
其中vx、vy、vz分別為x、y、z方向上質(zhì)點(diǎn)振動(dòng)速度分量;τxx、τyy、τzz分別為x、y、z方向上的正應(yīng)力;τxy、τyz、τxz為剪切應(yīng)力;ρ為介質(zhì)的密度;cab(a=1~6,b=1~6)是各向異性介質(zhì)的剛性系數(shù).特別地,對(duì)于TI介質(zhì),當(dāng)介質(zhì)對(duì)稱軸同z軸平行時(shí),式(4)中僅c11、c12、c13、c22、c23、c33、c44、c55和c66不為零,且滿足c12=c11-2c66、c23=c13、c22=c11和c44=c55,其它元素為零,這樣采用c11、c13、c33、c44、c66五個(gè)參數(shù)即可描述TI介質(zhì)中的波傳播現(xiàn)象.通過(guò)Bond變換可以獲得當(dāng)介質(zhì)對(duì)稱軸同z軸呈一定夾角時(shí)介質(zhì)的剛性系數(shù)矩陣[19].對(duì)于TI介質(zhì),當(dāng)介質(zhì)對(duì)稱軸在x-z平面內(nèi)且同z軸呈一定夾角時(shí),除上述幾個(gè)參數(shù)之外,c15、c25、c35和 c46也不為零.右側(cè)的gab(a、b=x~z)表示力變化速度的體積源,和體力源fi(i=x~z)組合使用可以模擬各種聲源.
圖2 交錯(cuò)網(wǎng)格1/8元胞示意圖Fig.2 Schematic for 1/8cell of staggered grid
我們采用了交錯(cuò)網(wǎng)格的方式來(lái)實(shí)現(xiàn)差分的顯式迭代過(guò)程.圖2為采用的交錯(cuò)網(wǎng)格1/8元胞示意圖.式(5)為速度和應(yīng)力各分量在空間和時(shí)間上的位置,其中l(wèi)表示網(wǎng)格的空間位置.正應(yīng)力各分量τxx,τyy,τzz均位于整數(shù)網(wǎng)格節(jié)點(diǎn)上,切應(yīng)力各分量和速度分量分別位于各自對(duì)應(yīng)的半整數(shù)網(wǎng)格節(jié)點(diǎn)上.對(duì)于本文研究的TI地層與井軸斜交的情況,由于剛性系數(shù)矩陣元素c15、c25、c35和c46不為0,根據(jù)網(wǎng)格上各個(gè)物理量之間的位置關(guān)系,僅采用交錯(cuò)網(wǎng)格無(wú)法進(jìn)行差分近似處理.如圖2所示,在采用式(4)計(jì)算時(shí),由于c15、c25、c35和c46不為零,需要網(wǎng)格 (lx,ly,lz+1/2)處的速度值以及網(wǎng)格 (lx+1/2,ly+1/2,lz+1/2)的速度值,網(wǎng) 格(lx+1/2,ly+1/2,lz+1/2)的速度值針對(duì)此問(wèn)題,一種解決方法是采用對(duì)速度場(chǎng)進(jìn)行插值的方法獲取上述點(diǎn)的速度值,另外一種方法是采用輔助交錯(cuò)網(wǎng)格的方法.本文采用了對(duì)速度場(chǎng)進(jìn)行插值的方法.
如式(6—9)所示,在計(jì)算時(shí),首先計(jì)算該網(wǎng)格點(diǎn)上的應(yīng)力值,應(yīng)力值τxx的計(jì)算方法如式(6)所示,其他應(yīng)力值的計(jì)算與式(6)類似,此處不做贅述.普通交錯(cuò)網(wǎng)格處的速度值仍舊按照式(1—3)進(jìn)行計(jì)算,網(wǎng)格上 (lx,ly,lz+1/2)處的速度值網(wǎng)格 (lx+1/2,ly+1/2,lz+1/2)的速度值以及網(wǎng)格 (lx+1/2,ly,lz)的速度值可以通過(guò)對(duì)速度場(chǎng)進(jìn)行插值的方法獲得,如式(7—9)所示.據(jù)此,可以完成差分算法的顯示迭代.
其中,式(7—9)中,I表示平移算子,其兩個(gè)下標(biāo)分別表述平移算子的空間階數(shù)和平移算子的方向.N為差分算子所采用的空間階數(shù)的一半.
采用本文的網(wǎng)格劃分方法對(duì)式(1—4,6—9)進(jìn)行離散化處理,得到顯式的差分迭代格式.式(10)和式(11)分別為速度分量vx和應(yīng)力分量τxx離散差分格式,其他分量及輔助交錯(cuò)網(wǎng)格各分量的迭代格式形式類似.
式(10~11)中,δx、δy、δz分別代表物理量在x、y、z方向的差分,Δt為計(jì)算采用的時(shí)間步長(zhǎng).
在直角坐標(biāo)系下,對(duì)于一般的各向異性介質(zhì),有限差分計(jì)算方法的穩(wěn)定性條件為[11]
式(12~13)中,vmax和vmin代表計(jì)算模型速度的最大值和最小值,am為采用的差分系數(shù),Δx、Δy、Δz分別代表x、y、z方向的空間步長(zhǎng),fmax為聲源覆蓋的最高頻率.
本文采用兩個(gè)緊貼鉆鋌外側(cè)振動(dòng)相位相反的點(diǎn)聲源來(lái)模擬偶極子聲源.點(diǎn)聲源的加載方法為式(14)所述形式[20]:
式(14)中,I3為三維直角坐標(biāo)系中網(wǎng)格的脈沖響應(yīng),sn為第n次迭代時(shí)加載的聲源值,δ為單位脈沖函數(shù),xs、ys和zs分別為聲源在x、y、z方向的坐標(biāo).聲源函數(shù)采用了雷克子波函數(shù),如式(15)所示.
在計(jì)算中,介質(zhì)的剛性系數(shù)和密度均賦在整數(shù)節(jié)點(diǎn)上,對(duì)于非整數(shù)網(wǎng)格點(diǎn)的物理量,通過(guò)臨近網(wǎng)格的物理量的平均得到.非整數(shù)網(wǎng)格點(diǎn)處的密度,通過(guò)式(16~18)所示的平均的方法獲得.非整數(shù)網(wǎng)格點(diǎn)處的剛性系數(shù)通過(guò)如式(19~21)所示的計(jì)算方法獲得.這樣對(duì)于固液界面(鉆鋌-流體邊界和流體-地層邊界),通過(guò)式(16~21)所給出的平均的方法,邊界條件自動(dòng)滿足.
為了模擬無(wú)限大的地層,采用了完全匹配層(PML)技術(shù)來(lái)吸收向地層內(nèi)傳播的波[18].PML層厚度選為地層縱波波長(zhǎng)的一半.
由于模型計(jì)算量較大,采用傳統(tǒng)的串行計(jì)算方法無(wú)法滿足計(jì)算需求.我們采用了OpenMP和MPI混合編程技術(shù),將有限差分算法在集群上實(shí)現(xiàn).MPI是目前在集群上應(yīng)用最為廣泛的并行計(jì)算技術(shù).OpenMP雖然僅適用于單機(jī)多核計(jì)算,但是其計(jì)算效率高,易于編程實(shí)現(xiàn),且目前大部分編譯器都已經(jīng)支持OpenMP技術(shù).本文通過(guò)采用OpenMP和MPI混合編程技術(shù),簡(jiǎn)化了并行算法的復(fù)雜性,提高了程序的執(zhí)行效率.如圖3所示,x和z分別代表直角坐標(biāo)系的x方向和z方向,m和n分別代表采用的進(jìn)程數(shù)和線程數(shù),雙向箭頭表示相鄰進(jìn)程之間的通信.通過(guò)合理的計(jì)算區(qū)域劃分,將計(jì)算任務(wù)分配到每個(gè)參與計(jì)算的節(jié)點(diǎn)上.通過(guò)MPI技術(shù),在每個(gè)計(jì)算節(jié)點(diǎn)上開(kāi)辟一個(gè)進(jìn)程,通過(guò)進(jìn)程間的通信和協(xié)作,實(shí)現(xiàn)計(jì)算的并行.在一個(gè)節(jié)點(diǎn)上,運(yùn)用OpenMP技術(shù),開(kāi)辟多個(gè)線程,利用多核協(xié)同工作,加快計(jì)算的速度.
圖3 并行計(jì)算方案Fig.3 Parallel implementation of the algorithm
對(duì)于240×240×300個(gè)網(wǎng)格,20000個(gè)時(shí)間步長(zhǎng)的數(shù)值模型,采用5個(gè)CPU核心數(shù)為12的節(jié)點(diǎn)進(jìn)行計(jì)算,每個(gè)節(jié)點(diǎn)開(kāi)辟的線程數(shù)目為11,采用雙精度進(jìn)行計(jì)算時(shí)完成計(jì)算所需的時(shí)間大約為20h.
圖4為數(shù)值計(jì)算模型示意圖,包括(a)模型主計(jì)算區(qū)域和(b)井孔橫截面示意圖.模型主計(jì)算區(qū)域的尺寸為1m×1m×4.8m,x、y、z方向的空間采樣間隔分別為0.0075m、0.0075m和0.0125m.井孔位于模型中央,井軸與z軸平行.介質(zhì)對(duì)稱軸S位于x-z平面內(nèi),介質(zhì)對(duì)稱軸與井軸夾角為α.特別地,當(dāng)α=0°時(shí),井軸與介質(zhì)對(duì)稱軸平行,相當(dāng)于豎直井井孔沿對(duì)稱軸穿過(guò)VTI地層,當(dāng)α=90°時(shí),井軸與介質(zhì)對(duì)稱軸垂直,相當(dāng)于豎直井井孔沿垂直于介質(zhì)對(duì)稱軸的方向穿過(guò)HTI地層.數(shù)值模擬時(shí)采用的地層參數(shù)為實(shí)驗(yàn)室內(nèi)測(cè)量的各向異性介質(zhì)的參數(shù),該介質(zhì)在TI對(duì)稱軸與z軸平行情況下的剛性參數(shù)如表1所示.井孔內(nèi)流體及鉆鋌參數(shù)見(jiàn)表2,鉆鋌內(nèi)徑、外徑及井眼直徑分別為0.054m、0.180m和0.240m.聲源加載在距離底界面0.8m處,采用在鉆鋌外徑處加載兩個(gè)震動(dòng)相位相反的點(diǎn)聲源的方法來(lái)模擬偶極子聲源.偶極子接收器同樣放置于鉆鋌外徑處,接收器源距為2.0~3.5m,間距為0.15m.由于本文重點(diǎn)研究的對(duì)象為地層彎曲波,不涉及隔聲及鉆鋌波問(wèn)題的研究,為壓制鉆鋌波,在發(fā)射器到源距最小的接收器之間將鉆鋌截?cái)?
表1 地層參數(shù)Table1 Formation parameters
表2 鉆鋌及鉆鋌內(nèi)外流體的參數(shù)Table 2 Parameters of the collar and the fluid in and out of the collar
圖4b為x-y平面內(nèi)井孔橫截面示意圖.為描述方便,定義地層橫向同性面和井軸垂直面之交線與偶極子聲源偏振方向的夾角為β,發(fā)射探頭和接收探頭所對(duì)應(yīng)的夾角分別為βT和βR.特別地,對(duì)于本文α=90°的井孔模型,β為偶極子聲源偏振方向同快橫波面的夾角.當(dāng)βT=βR時(shí),接收器與發(fā)射器偏振方向相同,測(cè)得波形為同向分量波形;當(dāng)βR=βT+90°時(shí),接收器與發(fā)射器偏振方向相差90°,測(cè)得波形為正交分量波形.特別地,當(dāng)βT=0°時(shí),聲源的偏振方向與地層中傳播的SH波偏振方向一致;當(dāng)βT=90°時(shí),聲源的偏振方向與地層中傳播的準(zhǔn)SV波偏振方向一致.
圖4 數(shù)值模擬采用的模型示意圖(a)及x-y截面示意圖(b)Fig.4 Schematic diagram of numerical simulation model(a)and diagram of x-ycross section(b)
本文首先模擬了α=90°,βT=0°、21.25°、45°、68.75°和90°情況下正交偶極子聲波測(cè)井,借以研究在井軸與TI介質(zhì)對(duì)稱軸垂直的情況下,隨鉆正交偶極子聲波測(cè)井對(duì)地層各向異性的評(píng)價(jià)能力,然后計(jì)算了α=0°、15°、30°、45°、60°、75°、90°,βT=0°、90°情況下偶極子聲源激發(fā)的聲場(chǎng),研究不同井斜情況下的隨鉆正交偶極子聲波測(cè)井的響應(yīng)特征.
圖5 α=90°時(shí),βT=0°、90°偶極子在井孔中激勵(lì)的偶極子波形,源距為2.0~3.5mFig.5 Dipole waveforms excited in the borehole by dipole source of source-receiver space 2.0~3.5m withα=90°,βT=0°、90°
數(shù)值模擬了α=90°,βT=0°、90°情況下偶極子聲源在井孔中激發(fā)的模式波.圖5為數(shù)值模擬結(jié)果,其中實(shí)線為聲源和接收器的方向βT,βR=0°時(shí)測(cè)得的偶極子波形,虛線為聲源和接收器的方向βT,βR=90°時(shí)測(cè)得的偶極子波形.可以看到,βT=0°時(shí)偶極子聲源激勵(lì)的彎曲波的傳播速度大于βT=90°時(shí)偶極子激勵(lì)的彎曲波的傳播速度,同對(duì)應(yīng)的地層體波SH波和SV波(在α=90°時(shí),縱波和SV波不耦合)的速度相一致.這表明,在鉆鋌存在的情況下,不同方向的偶極子聲源激發(fā)的彎曲波速度的差異同地層橫波速度的各向異性有關(guān).
數(shù)值模擬了α=90°,βT=0°、21.25°、45°和68.75°和90°情況下,偶極子聲源在含鉆鋌井孔內(nèi)激發(fā)的聲場(chǎng).圖6為不同角度下偶極子聲源在井孔中激發(fā)的同向分量波形(a)和正交分量波形(b).從圖6中可見(jiàn),當(dāng)βT=0°或者90°時(shí),幾乎接收不到正交分量波形信號(hào),當(dāng)βT=21.25°、45°、68.75°時(shí),正交分量能量較強(qiáng),且當(dāng)βT=45°時(shí)正交分量能量最強(qiáng);而同向分量波形能量在βT=0°或者90°時(shí)幅度較強(qiáng),在βT=45°時(shí)幅度相對(duì)較弱.同向分量和正交分量的能量變化表明,在βT=0°或者90°時(shí),彎曲波未發(fā)生分裂現(xiàn)象;βT=21.25°、45°、68.75°時(shí),彎曲波發(fā)生了分裂現(xiàn)象.圖6的模擬結(jié)果證實(shí)了,對(duì)于井軸同介質(zhì)對(duì)稱軸垂直的TI地層井孔,同電纜測(cè)井一致,在隨鉆條件下,彎曲波分裂現(xiàn)象仍然存在.
圖6 α=90°,βT=0°、21.25°、45°、68.75°和90°時(shí),模擬得到的(a)同向分量波形和(b)正交分量信號(hào),源距為2mFig.6 Simulated inline component waveforms(a)and cross-line component waveforms(b)of source-receiver space 2mwithα=90°,βT=0°、21.25°、45°、68.75°and 90°
圖7 α=90°,βT=68.75°時(shí),數(shù)值模擬得到的四分量偶極子波形,包括同向分量波形(a)XX 和(c)YY,以及正交分量波形(b)XY和(d)YX,源距為2~3.5mFig.7 Simulated four-component dipole waveforms including incline component waveforms(a)XXand(c)YY,and cross-line component waveforms(b)XYand(d)YXof source-receiver space 2~3.5mwithα=90°,βT=68.75°
假設(shè)正交偶極子聲源的兩個(gè)正交的方向分別標(biāo)記為X和Y.模擬了X方向同快橫波面的夾角分別為βT=0°、21.25°、45°、68.75°、90°幾種情況下正交偶極子聲源在含鉆鋌井孔內(nèi)激發(fā)的四分量偶極子波形.圖7為βT=68.75°情況下模擬的四分量偶極子波形,包括同向分量波形(a)(c)和正交分量波形(b)(d).對(duì)于各向同性地層,由于各方向地層聲學(xué)參數(shù)均相同,不會(huì)接收到正交分量信號(hào);對(duì)于各向異性地層,當(dāng)聲源偏振方向同介質(zhì)對(duì)稱軸呈一定夾角時(shí),由于各向異性地層的耦合作用,會(huì)接收到正交分量.從圖7中可以看到,正交分量信號(hào)XY和YX 均有較強(qiáng)的幅度,說(shuō)明發(fā)生了彎曲波分裂現(xiàn)象.綜上所述,在含鉆鋌TI地層井孔中,βT=0°、90°時(shí)偶極子聲源激勵(lì)的彎曲波未發(fā)生分裂現(xiàn)象,分別以較快、較慢的速度沿井軸傳播;當(dāng)聲源偏振方向同介質(zhì)對(duì)稱軸呈一定夾角的情況下,彎曲波分裂成以快、慢速度傳播的兩種模式波,能夠接收到較強(qiáng)幅度的正交分量信號(hào).
以上分析表明,同電纜測(cè)井條件下各向異性地層井孔中偶極子聲源激發(fā)的聲場(chǎng)類似,隨鉆條件下偶極子聲源激發(fā)的彎曲波也存在分裂現(xiàn)象,且βT=0°時(shí)偶極子聲源激勵(lì)的彎曲波的速度大于βT=90°時(shí)偶極子聲源激勵(lì)的彎曲波的速度.
快橫波面定義為快橫波偏振方向與井軸確定的平面.采用各向異性分析方法,通過(guò)四分量偶極子波形的旋轉(zhuǎn),從模擬得到的陣列波形中提取了快橫波面的方位,通過(guò)對(duì)比反演得到的方位角同正演模型采用的方位角,分析采用隨鉆正交偶極子測(cè)井進(jìn)行快橫波面方位測(cè)量的可行性.
圖8為α=90°井內(nèi),在快橫波面同聲源偏振方向的夾角βT=68.75°的情況下,通過(guò)Alford四分量波形旋轉(zhuǎn)方法[12]得到的正交分量相對(duì)能量隨儀器旋轉(zhuǎn)角度的變化圖,其極小值對(duì)應(yīng)著目的快橫波面方位.正交分量的相對(duì)能量定義為正交分量波形能量占四分量波形總能量的比例.從圖8中反演得到的快橫波面同聲源的偏振方向的夾角β′T為69.71°.
2.2.1 設(shè)置數(shù)據(jù)驅(qū)動(dòng)頁(yè)面 數(shù)據(jù)驅(qū)動(dòng)頁(yè)面(Data Driven Pages)是通過(guò)設(shè)置索引圖層生成多個(gè)輸出頁(yè)面的方法。索引圖層用于通過(guò)單個(gè)布局生成多個(gè)輸出頁(yè)面,每個(gè)頁(yè)面顯示不同范圍的數(shù)據(jù),范圍由索引圖層中的要素定義[2]。
圖9為采用圖8所示方法從模擬得到的陣列波形中反演得到的快橫波面的方位同模型實(shí)際采用的方位的對(duì)比圖,其中實(shí)線為模擬時(shí)實(shí)際的快橫波面同聲源的偏振方向的夾角βT,空心圓圈為從數(shù)值模擬的波形中反演得到的聲源的偏振方向同快橫波面的夾角β′T.從圖9中可見(jiàn),反演得到的快橫波面方位同模型實(shí)際的方位一致性非常好.數(shù)值模擬結(jié)果證明,在隨鉆條件下采用正交偶極子聲波測(cè)量方式能夠?qū)Φ貙拥目鞕M波方位進(jìn)行評(píng)價(jià).在這一點(diǎn)上,隨鉆條件下的正交偶極子聲波測(cè)井同電纜測(cè)井情況一致.
圖8 α=90°井內(nèi),βT=68.75°情況下,正交分量相對(duì)能量隨儀器旋轉(zhuǎn)角度變化圖Fig.8 The relative energy of cross-line component waveforms with the changes of rotation angle of tools,whenα=90°,βT=68.75°
圖9 從數(shù)值模擬得到的波形中提取的快橫波面方位同模型實(shí)際的方位的對(duì)比Fig.9 Comparison of the fast shear wave direction obtained from simulated waveforms with the actual direction
采用矩陣束方法[21]從α=90°,βT=0°、90°時(shí)偶極子聲源激發(fā)的波形中提取了各模式波的頻散曲線,提取的結(jié)果如圖10所示,其中黑色實(shí)線和虛線分別為βT=0°、90°時(shí)偶極子聲源激勵(lì)的彎曲波的頻散曲線,白色實(shí)線和虛線分別為該角度下地層體波SH波和SV波的慢度.圖10a和圖10b分別為βT=0°時(shí)偶極子聲源激勵(lì)的聲場(chǎng)的頻散圖和βT=90°時(shí)偶極子聲源激勵(lì)的聲場(chǎng)的頻散圖.人們?cè)贖TI地層井孔中的數(shù)值模擬結(jié)果表明,對(duì)于電纜測(cè)井,在大多數(shù)情況下,彎曲波的低頻傳播速度接近于對(duì)應(yīng)的快、慢橫波的傳播速度[10,15].從圖10可見(jiàn),與電纜測(cè)井相比,隨鉆條件下的彎曲波的頻散規(guī)律有兩點(diǎn)不同:在隨鉆條件下,彎曲波的低頻速度并非趨近于地層的橫波速度;彎曲波隨著頻率的變化并非單調(diào)變化.對(duì)于低頻段(0~1kHz)的彎曲波,其速度隨著頻率增加而增加,且βT=0°時(shí)偶極子聲源和βT=90°時(shí)偶極子聲源激發(fā)的彎曲波的速度差異不大;在頻率1.5~8kHz下,βT=0°時(shí)偶極子聲源激發(fā)的彎曲波的速度大于βT=90°時(shí)偶極子聲源激勵(lì)的彎曲波的傳播速度,二者差異較大,同電纜測(cè)井情況一致.數(shù)值模擬的結(jié)果表明,對(duì)于我們所研究的地層,由于在隨鉆條件下彎曲波的低頻速度不再趨近于地層的橫波速度,無(wú)法通過(guò)彎曲波的測(cè)量直接獲得地層橫波速度;不過(guò),在某些頻段內(nèi)βT=0°時(shí)偶極子聲源激發(fā)的彎曲波的速度大于βT=90°時(shí)偶極子聲源激發(fā)的彎曲波的速度,地層彎曲波的各向異性仍然能夠反應(yīng)地層橫波速度的各向異性.
在圖10中,可以觀察到隨頻率增加慢度變小的鉆鋌波和隨頻率降低慢度逐漸趨近于地層橫波慢度的六極子波.之所以能夠接收到六極子波,是因?yàn)樵诒疚牡哪M中,偶極子聲源裝在鉆鋌的外側(cè),極距較大,對(duì)于近場(chǎng)而言并非理想的偶極子源,激發(fā)的模式波中含有六極子和更高極性的成分[17].另外一點(diǎn),在圖10所示的頻散圖中,我們觀察到了以地層橫波的速度傳播的模式(圖中圓圈標(biāo)注區(qū)域).該模式在整個(gè)計(jì)算的頻率段內(nèi)都能夠觀察到,在低頻段和鉆鋌波混合在一起,難以區(qū)分.該模式的速度同地層橫波速度一致,當(dāng)βT=0°時(shí),該模式的速度同地層的SH波一致;當(dāng)βT=90°時(shí),該模式的速度同地層的SV波速度一致.SV波的速度為1399.3m/s,小于本文模型中的井內(nèi)流體速度1500m/s,這說(shuō)明,該模式能夠在軟地層中激勵(lì).通過(guò)對(duì)于該模式波速度的提取,可以直接獲得地層橫波時(shí)差及橫波速度的各向異性值.在偶極子聲源激發(fā)的六極子模式波中,未觀察到明顯的與地層各向異性相關(guān)的信息.
數(shù)值模擬了α=0°、15°、30°、45°、60°、75°、90°,βT=0°、90°情況下偶極子聲源激發(fā)的聲場(chǎng),采用慢度時(shí)間相關(guān)法(STC)提取了同向分量彎曲波的速度.圖11a為提取的彎曲波的速度與介質(zhì)的體波速度的對(duì)比圖,其中實(shí)線為對(duì)應(yīng)偏振方向的體波的速度,實(shí)心圓圈和空心圓圈標(biāo)注的曲線為彎曲波的速度.對(duì)于α=90°的井,βT=0°和βT=90°時(shí)模擬得到同向分量彎曲波的速度為1098.1m/s和1056.3m/s,均小于對(duì)應(yīng)的體波的速度1500.6m/s和1399.3m/s.從圖11a可見(jiàn),隨α變化,彎曲波的速度同體波速度變化趨勢(shì)不完全一致.βT=0°時(shí)偶極子聲源激勵(lì)的彎曲波的速度隨α的增加而增加,同地層SH體波的速度變換趨勢(shì)基本一致,βT=90°時(shí)偶極子聲源激勵(lì)的彎曲波的速度隨α的增加先減小后增加,而地層準(zhǔn)SV體波的速度隨α先增加后減小,二者變化趨勢(shì)不一致.
He等的研究結(jié)果表明,在電纜測(cè)井情況下,對(duì)于大多數(shù)橫向各向同性地層,彎曲波的低頻傳播速度接近于對(duì)應(yīng)的地層橫波的速度,且隨井斜傾角的變化同地層體波的變化基本一致[15].而在隨鉆條件下,對(duì)于本文所研究的地層模型,彎曲波的速度遠(yuǎn)小于地層橫波速度,且二者的變化趨勢(shì)也不一致.
通過(guò)式(22)計(jì)算了不同角度下通過(guò)彎曲波測(cè)得的各向異性值,其中時(shí)偶極子聲源激勵(lì)的彎曲波速度時(shí)偶極子聲源激勵(lì)的彎曲波的速度.圖11b中的空心圓圈和實(shí)線分別為通過(guò)彎曲波計(jì)算的各向異性值和地層真實(shí)的橫波速度的各向異性值.從圖11b中可見(jiàn),在斜井隨鉆條件下采用彎曲波獲得的各向異性值同地層的真實(shí)各向異性值差別很大:當(dāng)井軸方向同介質(zhì)對(duì)稱軸方向的夾角α小于45°時(shí),采用彎曲波獲得的各向異性值遠(yuǎn)大于地層的真實(shí)的各向異性值;當(dāng)α大于45°時(shí),隨著α的增加,二者逐漸接近,在α為65°左右時(shí),二者的值相當(dāng),然后隨α的增加,二者差異變大.當(dāng)α大于60°時(shí),雖然二者隨α的變化趨勢(shì)不一致,但是彎曲波的各向異性值基本能夠反映地層橫波速度的各向異性.此種情況下,可以考慮采用合適的反演方法來(lái)獲取地層的各向異性信息.
模式波關(guān)于某參數(shù)p的靈敏度定義為歸一化后的波的相速度對(duì)該參數(shù)的偏微分,見(jiàn)式(23)[1].通過(guò)靈敏度分析可以幫助理解各參數(shù)對(duì)模式波速度的影響程度.由于各向異性地層井孔聲場(chǎng)不存在解析解,本文選取各向同性地層參數(shù)來(lái)分析隨鉆測(cè)井聲學(xué)模型各參數(shù)對(duì)彎曲波速度的影響程度.選取的地層縱橫波速度分別為3227.2m/s和1399.3m/s,其它模型參數(shù)同上文中的各向異性地層模型的一致.圖12a和圖12b分別為理論計(jì)算的相速度曲線、群速度曲線和各參數(shù)(地層縱波速度、地層橫波速度、鉆鋌縱波速度、鉆鋌橫波速度、井內(nèi)流體速度)的靈敏度曲線.從圖12a可見(jiàn),理論相速度曲線的變化規(guī)律與從數(shù)值模擬得到的頻散曲線的變化規(guī)律一致.從圖12b靈敏度曲線可見(jiàn),在隨鉆條件下,地層橫波對(duì)彎曲波的控制作用相對(duì)電纜測(cè)井減小,彎曲波除了受地層橫波控制之外,還受到鉆鋌和其它地層參數(shù)的影響.對(duì)于本文的模型,在低于1kHz的頻率內(nèi),鉆鋌橫波速度為彎曲波的最主要的控制因素;在高于1.5kHz內(nèi),地層橫波速度是影響彎曲波的最主要因素.這解釋了圖10中在低頻段(0~1kHz)βT=0°時(shí)偶極子聲源和βT=90°時(shí)偶極子聲源激發(fā)的彎曲波速度基本無(wú)差異,而在高頻段(1.5kHz)二者差異較大的現(xiàn)象.這說(shuō)明,盡管彎曲波的低頻速度并非接近地層橫波速度,但是仍舊可以考慮結(jié)合理論的頻散曲線從彎曲波頻散信息中提取地層的橫波速度.對(duì)于各向異性地層,可以考慮首先選擇合適的頻率區(qū)間,然后分別從快、慢彎曲波波形的頻散信息中提取出地層的快、慢橫波速度,進(jìn)一步地可以計(jì)算地層橫波速度的各向異性值.
在隨鉆條件下,由于鉆鋌的存在,充液井孔中傳播的模式波的特征與電纜測(cè)井非常不同[1,4,17],地層聲學(xué)各向異性在正交偶極子聲波測(cè)井中的響應(yīng)特征也與電纜測(cè)井不同.本文的模擬結(jié)果驗(yàn)證了井孔與TI介質(zhì)對(duì)稱軸垂直情況下,采用隨鉆正交偶極子聲波測(cè)井方式進(jìn)行地層快橫波方位角評(píng)價(jià)和地層橫波速度的各向異性評(píng)價(jià)的可行性.另外,在以下幾個(gè)方面,還需要進(jìn)行進(jìn)一步的研究工作.
同電纜測(cè)井相比,隨鉆條件下的彎曲波的傳播速度在低頻范圍內(nèi)并非趨近于地層的橫波速度,無(wú)法實(shí)現(xiàn)直接橫波測(cè)井.盡管上文中證實(shí)了地層橫波速度在一定的頻率區(qū)間內(nèi)仍是影響彎曲波的最主要因素,可以考慮采用結(jié)合正演模型的反演方法來(lái)提取地層的快、慢橫波速度及各向異性值,但是正演模型的誤差,正演模型輸入?yún)?shù)的誤差以及測(cè)量的彎曲波速度的誤差等都可能導(dǎo)致反演得到的橫波速度及各向異性值存在較大的誤差.可以考慮結(jié)合四極子聲源的測(cè)量結(jié)果或者偶極子波形中的六極子模式的處理結(jié)果,對(duì)地層的快、慢橫波速度及各向異性值進(jìn)行聯(lián)合反演,這會(huì)在一定程度上提高反演結(jié)果的可靠性.另外,建立適用于反演的隨鉆條件下的TI地層井孔多極子聲源激發(fā)模式波的頻散模型是實(shí)現(xiàn)上述反演算法的關(guān)鍵.
當(dāng)井孔與TI介質(zhì)對(duì)稱軸存在一定夾角時(shí),情況變得非常復(fù)雜,很難應(yīng)用于現(xiàn)場(chǎng)資料的解釋.雖然,當(dāng)α大于60°時(shí),彎曲波的各向異性值能夠基本反映該角度下地層橫波速度的各向異性,但是目前的隨鉆測(cè)井技術(shù)很難準(zhǔn)確測(cè)得α,仍舊不能夠得到地層真實(shí)的快、慢橫波速度及橫波速度的各向異性值.另外,當(dāng)α大于60°時(shí),彎曲波的各向異性值同地層真實(shí)橫波速度的各向異性值隨α的變化趨勢(shì)也不一致,這也給反演結(jié)果帶來(lái)了一定的不確定性.
在圖10的頻散分析結(jié)果中,我們觀察到了以地層橫波速度傳播的模式波,而且對(duì)于βT=0°的偶極子聲源,該模式波速度與地層SH波一致,對(duì)于βT=90°的偶極子聲源,該模式波速度同地層SV波一致.這意味著,可以通過(guò)偶極子聲源實(shí)現(xiàn)直接橫波測(cè)量和各向異性測(cè)量.不過(guò)該模式幅度相對(duì)較低,受鉆鋌波以及鉆進(jìn)過(guò)程中的噪聲影響較大,如果想實(shí)現(xiàn)工程化應(yīng)用,需要在隔聲、采集以及信號(hào)處理方面做進(jìn)一步的工作.
本文的研究中未涉及鉆進(jìn)過(guò)程中的噪聲、儀器的偏心以及鉆鋌波對(duì)測(cè)量結(jié)果的影響,實(shí)際測(cè)井情況下,這些因素對(duì)測(cè)量結(jié)果有著較大的影響[1].針對(duì)這些問(wèn)題,需要在儀器設(shè)計(jì)及信號(hào)處理等方面進(jìn)行進(jìn)一步的研究工作.
我們采用三維有限差分方法對(duì)TI地層井孔中隨鉆正交偶極子聲波測(cè)井進(jìn)行了數(shù)值模擬,通過(guò)慢度時(shí)間相關(guān)法、頻散分析方法及各向異性分析方法對(duì)模擬得到的波形進(jìn)行了分析,驗(yàn)證了在隨鉆條件下的HTI地層井孔中進(jìn)行橫波速度的各向異性評(píng)價(jià)以及快橫波方位測(cè)量的可靠性,并對(duì)介質(zhì)對(duì)稱軸同井軸斜交情況下偶極子聲源激發(fā)的聲場(chǎng)進(jìn)行了計(jì)算,主要得到以下認(rèn)識(shí):
(1)在隨鉆測(cè)井條件下,對(duì)于井軸與TI地層對(duì)稱軸垂直的情況,當(dāng)偶極子聲源偏振方向同介質(zhì)對(duì)稱軸呈一定夾角時(shí),井孔中的彎曲波能夠分裂成以快慢兩種速度傳播的模式波,證實(shí)了鉆鋌存在情況下各向異性地層井孔內(nèi)的彎曲波分裂現(xiàn)象.通過(guò)對(duì)模擬得到的四分量偶極子波形的處理證實(shí)了正交偶極子測(cè)量方式評(píng)價(jià)地層快橫波方位的可行性.
(2)在介質(zhì)對(duì)稱軸同井軸斜交情況下的TI地層井孔內(nèi),偶極子聲源激發(fā)的聲場(chǎng)非常復(fù)雜.不過(guò),當(dāng)α大于60°時(shí)彎曲波的各向異性值基本能夠反映地層橫波速度的各向異性.
(3)隨鉆條件下井孔內(nèi)彎曲波頻散特征同電纜測(cè)井情況下的頻散特征非常不同,但是在某些頻率區(qū)間內(nèi),地層橫波速度仍然是控制彎曲波的最主要的因素.實(shí)際測(cè)量時(shí),由于鉆鋌模型固定,可以考慮選擇與正演理論相結(jié)合的反演算法來(lái)提取地層的快、慢橫波速度及各向異性值.
(4)頻散分析結(jié)果顯示,偶極子波形中存在以地層橫波速度傳播的模式,而且該模式的速度與聲源的偏振方向有關(guān).對(duì)于βT=0°的偶極子聲源,該模式波速度與地層SH波一致,對(duì)于βT=90°的偶極子聲源,該模式波速度同地層SV波一致.這意味著,可以通過(guò)偶極子聲源實(shí)現(xiàn)直接橫波測(cè)量和各向異性測(cè)量.
致 謝 感謝匿名審稿專家提出的寶貴意見(jiàn)和建設(shè)性建議.
(References)
[1]Tang X M,Cheng C H A.Quantitative Borehole Acoustic Methods.San Diego:Elsevier Science Publishing Co.Inc.,2004.
[2]Minear J,Birchak R,Robbins C,et al.Compressional slowness measurements while drilling.// 36thAnnual Logging Symposium Transactions,Society of Professional Well Log Analysis,1995.
[3]Leggett J VIII,Dubinsky V,Patterson D,et al.Field test result demonstrating improved resl-time data quality in an advanced LWD acoustic system.SPE Annual conference and Exhibition,New Orleans,Louisiana,USA,2001.
[4]Tang X M,Wang T,Patterson D.Multipole acoustic logging-while-drilling.//72ndAnn.Internat.Mtg.,Soc.Explor.Geophys.,Expanded Abstracts,2002:364-368.
[5]Sinha B K,Simsek E,Asvadurov S.Influence of a pipe tool on borehole modes.Geophysics,2009,74(3):E111-E123.
[6]Wang T,Tang X M.LWD quadrupole shear measurement in anisotropic formations.//73ndAnn.Internat.Mtg.,Soc.Explor.Geophys.,Expanded Abstracts.
[7]喬文孝,王瑞甲,車小花等.利用國(guó)產(chǎn)多極子聲波測(cè)井儀評(píng)價(jià)地層各向異性的實(shí)例分析.聲學(xué)技術(shù),2008,27(5):270-271.Qiao W X,Wang R J.Che X H,et al.Case study of anisotropy evaluation by homemade multipole acoustic logging tool.Technical Acoustics (in Chinese),2008,27(5):270-271.
[8]Chen N Y,Toks?z M N,Cheng C H.Borehole wave propagation in isotropic and anisotropic media:threedimensional finite difference approach.Department of Earth,Atmospheric,and Planetary Sciences,MIT,1994.
[9]Schmitt D P.Acoustic multipole logging in transversely isotropic poroelastic formations.JASA.1989,86(6):2397-2421.
[10]Sinha B K,Simsek E,Liu Q.Elastic-wave propagation in deviated wells in anisotropic formations.Geophysics,2006,71(6):D191-D202
[11]張海瀾,王秀明,張碧星.井孔中的聲場(chǎng)與波.北京:科學(xué)出版社,2004.Zhang H L,Wang X M,Zhang B X.Acoustic Field and Waves in Borehole (in Chinese).Beijing:Science Press,2004.
[12]陳雪蓮,王瑞甲.徑向分層TI孔隙介質(zhì)井孔中激發(fā)的模式波的數(shù)值研究,地球物理學(xué)報(bào),2008,51(4):1270-1277.Chen X L,Wang R J.A numerical study on the mode waves excited by multipole sources in the fluid-filled borehole in radially layered transversely isotropic porous medium.Chinese J.Geophys.(in Chinese),2008,51(4):1270-1277.
[13]陳雪蓮,王瑞甲.橫向各向同性彈性地層井孔中模式波的探測(cè)深度.吉林大學(xué)學(xué)報(bào) (地球科學(xué)版),2008,38(3):502-507.Chen X L,Wang R J.Investigating depth of mode waves in the borehole surrounded by transversely isotropic elastic formation.Journal of Jilin University(Earth Science Edition)(in Chinese),2008,38(3):502-507.
[14]He X,Hu H S,Guan W.Fast and slow flexural waves in a deviated borehole in a homogeneous or layered anisotropic formation.Geophys.J.Int.,2010,181:417-426.
[15]He X,Hu H S.Borehole flexural modes in anisotropic formations:The low-frequency asymptotic velocity.Geophysics,2009,74(4):E149-E158.
[16]閆守國(guó),宋若龍,呂偉國(guó)等.橫向各向同性地層斜井中正交偶極子激發(fā)的聲場(chǎng)的數(shù)值模擬.地球物理學(xué)報(bào),2011,54(9):2412-2418.Yan S G,Song R L,LüW G,et al.Numerical simulation of acoustic fields excited by cross-dipole source in deviated wells in transversely isotropic formation.Chinese J.Geophys(in Chinese),2011,54(9):2412-2418.
[17]崔志文.多孔介質(zhì)聲學(xué)模型與多極源聲電效應(yīng)測(cè)井和多極隨鉆聲測(cè)井的理論與數(shù)值研究[博士論文].長(zhǎng)春:吉林大學(xué),2004.Cui Z W.Theoretical and numerical study of modified Biot′s models,acoustoelectric well logging and acoustic logging while drilling excited by multipole acoustic sources[Ph.D.thesis](in Chinese).Changchun:Jilin University,2004.
[18]Chew W C, Liu Q H.Perfectly matched layers for elastodynamics:A new absorbing boundary conditions.J.Computat.Acoustics,1996,4:341-359.
[19]梁鍇.TI介質(zhì)地震波傳播特征與正演方法研究[博士論文].東營(yíng):中國(guó)石油大學(xué) (華東),2009.Liang K.The study on propagation feature and forward modeling of seismic wave in TI media[Ph.D.thesis](in Chinese).Dongying:China University of Petroleum (East China),2009.
[20]Schneider J B,Wagner C L,Broschat S L.Implementation of transparent sources embedded in acoustic finite-difference time-domain grids.J.Acoust.Soc.Am.,1998,103(1):136-142.
[21]Ekstrom M P.Dispersion estimation from borehole acoustic arrays using a modified matrix pencil algorithm.//Proceedings of the 29th Asilomar Conference on Signals,Systems and Computers.IEEE Computer Society,1995:449-453.