徐逸鶴， 徐濤， 王敏玲， 白志滕吉文

1 中國科學(xué)院地質(zhì)與地球物理研究所，巖石圈演化國家重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 中國科學(xué)院大學(xué)， 北京 100049 3 中國科學(xué)院青藏高原地球科學(xué)卓越創(chuàng)新中心， 北京 100101

?

井中震源的遠(yuǎn)場(chǎng)波場(chǎng)特征研究

1 中國科學(xué)院地質(zhì)與地球物理研究所，巖石圈演化國家重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 中國科學(xué)院大學(xué)， 北京 100049 3 中國科學(xué)院青藏高原地球科學(xué)卓越創(chuàng)新中心， 北京 100101

井中震源在逆VSP、隨鉆地震和采礦地球物理研究中都有廣泛應(yīng)用.滿足“小井孔”(井孔半徑遠(yuǎn)小于特征波長(zhǎng))及“遠(yuǎn)場(chǎng)”(炮檢距大于特征波長(zhǎng))假設(shè)時(shí)，井中震源的遠(yuǎn)場(chǎng)波場(chǎng)存在解析解.為了檢驗(yàn)解析解在不同情況下的適用性，本文使用最速下降積分計(jì)算了不滿足上述假設(shè)時(shí)井中震源遠(yuǎn)場(chǎng)波場(chǎng)的合成地震記錄，即半解析解.模型試驗(yàn)表明，解析解只能在同時(shí)滿足“小井孔”和“遠(yuǎn)場(chǎng)”假設(shè)時(shí)使用；當(dāng)這兩個(gè)假設(shè)條件不滿足時(shí)，解析解的振幅和波形相對(duì)于半解析解會(huì)有明顯的偏差.隨著假設(shè)不滿足程度的增加，偏差會(huì)逐漸增加，并會(huì)逐漸影響走時(shí)的準(zhǔn)確拾??；這種條件下，采用半解析解才能獲得準(zhǔn)確的井中震源波場(chǎng).

井中震源; 遠(yuǎn)場(chǎng)波場(chǎng); 解析解; 最速下降積分; 最速下降法

1 引言

井中震源是指具有圓柱形結(jié)構(gòu)的震源，常用于人工震源的地震探測(cè)和勘探中，例如深地震測(cè)深、地震勘探和采礦地球物理等領(lǐng)域.常見的井中震源包括炸藥震源，井下的空氣槍、水槍、射孔槍、以及鉆頭震源等(Chen et al., 1990).鉆井的存在給井中震源問題帶來了較為復(fù)雜的邊界條件，使得井中震源波場(chǎng)體現(xiàn)出與常規(guī)震源不同的特性(Lee and Balch, 1982; Meredith et al., 1993).

以震源類型劃分，井中震源可以大致分為三種：即徑向應(yīng)力源、軸向應(yīng)力源和旋轉(zhuǎn)應(yīng)力源，其他震源可以由這三種震源的組合而成.采礦地球物理中的炸藥震源可以簡(jiǎn)化為徑向應(yīng)力源(Blair, 2007, 2010)，隨鉆地震勘探中的鉆頭震源可以簡(jiǎn)化成軸向應(yīng)力源和旋轉(zhuǎn)應(yīng)力源的組合(Rector and Hardage, 1992).以研究區(qū)域劃分，井中震源的波場(chǎng)可以分為井內(nèi)波場(chǎng)和井外波場(chǎng).井內(nèi)波場(chǎng)作為聲波測(cè)井的主要研究區(qū)域，前人已經(jīng)對(duì)其有大量而詳實(shí)的研究(Tsang and Rader, 1979; Cheng and Toks?z, 1981; Tubman et al., 1984; Cheng, 1994; 沈建國和張海瀾，2000)；而井外波場(chǎng)的研究還相對(duì)較少(Heelan, 1953; Lee and Balch, 1982; Meredith, 1991; 劉銀斌等，1993a, b; 張釙等，1995; Blair, 2007).但是隨著采礦地球物理和隨鉆地震的發(fā)展，井外波場(chǎng)研究的重要性日益明顯(Rector and Marion, 1991; Rector and Hardage, 1992; Haldorsen et al., 1995; Poletto, 2005; Vasconcelos and Snieder, 2008; 陸斌等，2009；王鵬等，2009；黃偉傳等2010；吳何珍等，2010).

井中震源的波場(chǎng)研究方法主要分為數(shù)值模擬方法(Cheng, 1994; Dong and Toks?z, 1995)和解析法(Heelan, 1953; Tsang and Rader, 1979; Lee, 1986; Meredith et al., 1993; De Hoop et al., 1994; Blair, 2007)兩類.而井外波場(chǎng)問題的計(jì)算區(qū)域(如炮檢距為1000m)和研究對(duì)象(如鉆井半徑為0.1m)的尺度相差很大，因此在使用數(shù)值模擬方法時(shí)，如果要求網(wǎng)格尺寸逼近鉆井半徑，波場(chǎng)模擬效果較好，但計(jì)算量會(huì)過大；如果網(wǎng)格尺寸過大則不能很好地研究井中震源的地震波場(chǎng)特征(王鵬等，2009).所以，解析法成為目前研究井外波場(chǎng)的主要手段.

Heelan(1953)在遠(yuǎn)場(chǎng)近似下得到了有限長(zhǎng)度的三種基本井中震源的遠(yuǎn)場(chǎng)解析解.Jordan(1962)和Abo-Zena(1977)使用不同的方法也獲得了井中震源的遠(yuǎn)場(chǎng)解析解.Lee和Balch(1982)考慮了鉆井中流體的存在，得到含液井的遠(yuǎn)場(chǎng)解析解.劉銀斌等(1993a，b)將多個(gè)震源的解疊加起來，得到了井中震源陣列的遠(yuǎn)場(chǎng)解析解.上述研究均基于兩個(gè)假設(shè)條件：(1)井孔半徑遠(yuǎn)小于特征波長(zhǎng)；(2)炮檢距大于特征波長(zhǎng).本文中我們分別稱之為“小井孔”和“遠(yuǎn)場(chǎng)”假設(shè).Meredith(1991)和Blair(2007)使用離散波數(shù)法對(duì)頻率波數(shù)域的解進(jìn)行數(shù)值積分，得到了精確的遠(yuǎn)場(chǎng)波場(chǎng)，并分別應(yīng)用于井間地震和采礦地球物理的研究.離散波數(shù)法通過增加等間距的虛擬震源，將積分轉(zhuǎn)化為求和，是計(jì)算波數(shù)域積分的常用方法(Bouchon and Aki, 1977; Bouchon, 1978, 2003).但是波數(shù)域積分中極點(diǎn)(pole)和支線(branch cut)所代表的豐富物理意義也隨之消失(Lapwood, 1949).為了保留波數(shù)域積分的優(yōu)勢(shì)，井內(nèi)波場(chǎng)的研究通常采用實(shí)軸積分法(Tsang and Rader, 1979; 沈建國和張海瀾，2000).但是當(dāng)研究井外波場(chǎng)時(shí)，由于炮檢距一般遠(yuǎn)大于鉆孔半徑，實(shí)軸積分的積分函數(shù)會(huì)出現(xiàn)高頻振蕩的現(xiàn)象，嚴(yán)重影響數(shù)值積分的精度和計(jì)算量，不適用于井外波場(chǎng)的計(jì)算.

為了克服波數(shù)域?qū)嵼S積分中出現(xiàn)的高頻振蕩現(xiàn)象，本文提出一種基于最速下降路徑的數(shù)值積分方法，并通過比較數(shù)值積分和解析解之間的誤差，探討遠(yuǎn)場(chǎng)解析解的適用性條件.

2 井中震源波場(chǎng)解析解

2.1 頻率波數(shù)域的解析解

假設(shè)在各向同性彈性全空間中存在一個(gè)無限長(zhǎng)的圓柱形空洞(半徑為a)，用來模擬鉆井或者炮眼.井內(nèi)的震源通過圓柱形的邊界對(duì)井外空間產(chǎn)生的影響，通?？梢杂镁趦蓚?cè)位移或者應(yīng)力的連續(xù)性來傳遞.本文中我們考慮井中為真空，井壁為自由邊界條件的情況，因此只需考慮應(yīng)力源對(duì)井外空間的作用.

為了便于邊界條件的描述，我們采用柱坐標(biāo)系，并將柱坐標(biāo)系的z軸與圓柱形空洞的對(duì)稱軸重合(圖1a).在柱坐標(biāo)系下，邊界上的應(yīng)力可以用σrr，σr θ，σrz三個(gè)分量來表示，它們分別對(duì)應(yīng)于徑向、旋轉(zhuǎn)和軸向應(yīng)力源(圖1b).

柱坐標(biāo)系下的平衡方程(忽略體力項(xiàng))為

圖1 井中震源示意圖(a)震源結(jié)構(gòu).震源S呈軸對(duì)稱狀分布于半徑為a的井壁上，檢波器位于P(r,z)點(diǎn).(b)三種基本的井中震源.Fig.1 Diagram of downhole seismic sources(a) Source geometry. The source S is distributed axisymmetrically around the borehole wall of radius a. The receiver is placed at point P(r, z). (b) Three basic types of downhole seismic sources.

(1)

(2)其中λ,μ是Lamé常數(shù).幾何關(guān)系是

(3)

使用Helmholtz分解和極環(huán)分解，我們可以將位移場(chǎng)u=(ur,uθ,uz)T分成三個(gè)標(biāo)量場(chǎng)的組合(Meredith, 1991)：

(4)

其中φ,ψ,為標(biāo)量函數(shù)，稱為勢(shì)函數(shù)，分別對(duì)應(yīng)P波，SV波和SH波.將(2)—(5)式代入(1)式，整理可得，

(5)

(6)

其中b1,b2,b3分別對(duì)應(yīng)于徑向、旋轉(zhuǎn)和軸向應(yīng)力源.將波動(dòng)方程(5)和邊界條件(6)變換到頻率波數(shù)域，求解這個(gè)邊值問題，可得頻率波數(shù)域的井中震源解：

(7)

(8)

而d,Lij分別為

(10)

前人關(guān)于井中震源的解析解研究中，盡管采用的具體方法不盡相同，但最終都會(huì)得到與(9)式等價(jià)的解(Heelan, 1953; Abo-Zena, 1977; Lee and Balch, 1982; Meredith, 1991; 劉銀斌等，1993a, b; Blair, 2007).而對(duì)(10)式中雙重積分的計(jì)算方法則主要分為兩類.一類是在遠(yuǎn)場(chǎng)的假設(shè)下，解析地計(jì)算兩個(gè)積分，得到遠(yuǎn)場(chǎng)解析解(Heelan, 1953; Abo-Zena, 1977; Lee and Balch, 1982)；另一類是采用數(shù)值積分法計(jì)算，得到半解析解(Meredith, 1991; Blair, 2007).

2.2 近似條件下的遠(yuǎn)場(chǎng)解析解

(11)

其中γ是Euler常數(shù)，Γ(n)是gamma函數(shù).并且可以進(jìn)一步推知，

(12)

(13)

假設(shè)三種應(yīng)力源的形式都是δ(z)G(t)，其中G(t)是震源時(shí)間函數(shù).那么對(duì)(13)式作反Fourier變換，可以得到這三類震源在物理空間域的解.

(1)徑向應(yīng)力源：

(14)

(2)軸向應(yīng)力源：

(15)

(3)旋轉(zhuǎn)應(yīng)力源：

φG′(t-R/β).

(16)

其中G′(t)為震源時(shí)間函數(shù)的導(dǎo)數(shù)，是徑向和旋轉(zhuǎn)應(yīng)力源產(chǎn)生的遠(yuǎn)場(chǎng)波形，而軸向壓力源產(chǎn)生的遠(yuǎn)場(chǎng)波形則正比于震源時(shí)間函數(shù)G(t).

在徑向和軸向應(yīng)力源的位移表示為P波和SV波的組合.為了分離P和SV的貢獻(xiàn)，我們引入坐標(biāo)系(R,φ,θ)，其中

r=Rcosφ,z=Rsinφ,

(17)

那么R,φ方向的位移可以通過r,z方向的位移旋轉(zhuǎn)得到

uR=urcosφ+uzsinφ,uφ=-ursinφ+uzcosφ,

(18)

在新的坐標(biāo)系下，徑向和軸向應(yīng)力源的位移分別為

和

(20)

考慮到本文中φ的定義與前人使用的φL的關(guān)系φ=φL-π/2，上述公式與前人的結(jié)果一致(Heelan, 1953;LeeandBalch, 1982).從(19)、(20)式中可以看出，兩種震源激發(fā)的SV波能量都比P波能量大，并且能量差會(huì)隨著兩種波速差的增大而增大(圖2).

圖2 井中震源的遠(yuǎn)場(chǎng)輻射圖樣(修改自Lee and Balch, 1982；參數(shù)相同)震源分別是(a)施加在裸眼井井壁上的徑向壓力源，(b)施加在充液井井壁上的徑向壓力源，和(c)位于充液井對(duì)稱軸上的單極聲波源.Fig.2 Far-field radiation patterns for downhole seismic sources (modified after Lee and Balch, 1982; same parameters are used)The sources are respectively (a) a radial stress source applied on an empty borehole, (b) a radial stress source applied on a fluid-filled borehole, and (c) a monopole acoustic source placed at the center of a fluid-filled borehole.

3 井中震源波場(chǎng)半解析解

為了計(jì)算最速下降路徑積分，我們首先求解最速下降路徑的解析表達(dá)式，然后采用數(shù)值積分方法，沿著最速下降路徑計(jì)算頻率波數(shù)域解的積分.

理想情況下，計(jì)算F(h)的最速下降積分路徑應(yīng)遵循如下步驟：

(21)

其中積分函數(shù)exp(lnF(h))的振蕩是由lnF(h)虛部的變化導(dǎo)致的.而沿lnF(h)的虛部等值線進(jìn)行積分時(shí)，lnF(h)虛部保持不變，積分函數(shù)的振蕩性就會(huì)完全消失.從復(fù)平面上的大多數(shù)點(diǎn)出發(fā)沿虛部等值線積分，兩個(gè)方向的積分函數(shù)分別是指數(shù)增加和指數(shù)減少.但是對(duì)于某些特殊的點(diǎn)來說，兩個(gè)方向的積分函數(shù)都以指數(shù)形式減少.我們稱這些點(diǎn)為鞍點(diǎn)(saddlepoint)，記為hs.它們滿足如下條件：

(22)

從鞍點(diǎn)出發(fā)的所有方向中，沿虛部等值線兩個(gè)方向，積分函數(shù)的減小是最快的.因此，這條經(jīng)過鞍點(diǎn)的lnF(h)虛部等值線被稱為最速下降路徑(SteepestDescentPath，簡(jiǎn)稱SDP).SDP上任意一點(diǎn)h都滿足

lnF(h)=lnF(hs)-X2,

(23)

其中X是任意正實(shí)數(shù).已知hs就可利用上式計(jì)算出整條最速下降路徑.

(24)

該函數(shù)在數(shù)值計(jì)算中被稱為歸一化Hankel函數(shù)(scaledHankelfunction).代入(7)式可得

記

(26)

其中

取Imq>0∪(Imq=0∩ωReq>0)，

(27)

其中c=α,β.與近似條件下的解析解中使用的f(h)=i(qr+hz)不同，我們?cè)谧钏傧陆德窂降挠?jì)算中考慮了井孔半徑a的影響.令f′(h)=0，可以得到

圖3 最速下降積分和實(shí)軸積分的對(duì)比圖(a)(c)分別是實(shí)軸積分和最速下降積分在復(fù)波數(shù)平面內(nèi)的積分路徑，圖(b)(d)則是兩個(gè)路徑上的積分函數(shù).波數(shù)h的單位是m-1.Fig.3 Comparison of Steepest Descent Integration Method and Real-axis Integration Method Panels (a) (c) show two different integration paths for Real-axis Integration Method and Steepest Descent Integration Method in complex wavenumber plane, while (b) (d) are the corresponding integrands along the paths.

(28)

]=0,

(29)

方程解為

(30)

其中Δ是(29)式的判別式.當(dāng)X=0時(shí)，

(31)

(32)

因此，h1所代表的分支在鞍點(diǎn)hs的右側(cè)，我們稱之為SDP的右支，記為C1；同理h2是SDP的左支，記為C2.C1,2的正方向?yàn)閄2增加的方向，所以最速下降路徑CSDP=-C2+C1.結(jié)合復(fù)變函數(shù)中的Jordan引理，可以將實(shí)軸上的積分完全變成沿SDP的積分，如下式所示:

(33)

由于積分函數(shù)以exp(-X2)衰減，根據(jù)計(jì)算精度要求，可以很容易地確定上式的積分截?cái)嗌舷?同時(shí)，由于積分函數(shù)非常光滑，較低的采樣率也可以得到相當(dāng)高的精度.

4 近似條件下遠(yuǎn)場(chǎng)解析解的適用性特征4.1 最速下降積分法試驗(yàn)

我們首先驗(yàn)證最速下降積分法在計(jì)算遠(yuǎn)場(chǎng)波場(chǎng)時(shí)的適用性，采用如下參數(shù)模型：井孔半徑為0.1 m，檢波器放置在距離鉆井1000 m處；井外地層為常見的Pierre頁巖，其P波波速為2074 m·s-1，S波為869 m·s-1，密度為2250 kg·m-3(Meredith, 1991)；震源為徑向應(yīng)力震源，震源時(shí)間函數(shù)是主頻為30 Hz的Ricker子波(圖4a).該主頻的震源在這種地層中的特征波長(zhǎng)λm約為346 m，既遠(yuǎn)大于井孔半徑2 m，又遠(yuǎn)小于1000 m，同時(shí)滿足“小井孔”和“遠(yuǎn)場(chǎng)”假設(shè)，近似條件下的解析解可以準(zhǔn)確地得到位移解.

從徑向應(yīng)力源的遠(yuǎn)場(chǎng)解析解(14)式可知，遠(yuǎn)場(chǎng)的波形正比于Ricker子波的導(dǎo)數(shù)(圖4b).圖4c是地震記錄的徑向分量，其中實(shí)線代表半解析解，虛線代表解析解，倒三角形標(biāo)注了P波的理論到時(shí).數(shù)值試驗(yàn)結(jié)果表明，用最速下降路徑積分得到的數(shù)值解與近似條件下的解析解吻合很好(圖4c).由于在本例條件下，解析解對(duì)位移刻畫較為準(zhǔn)確，所以這一結(jié)果進(jìn)一步驗(yàn)證了最速下降積分法的正確性.

4.2 徑向應(yīng)力震源

4.2.1 高頻情況

徑向應(yīng)力震源通常用來模擬炮眼中的爆炸震源，典型的炮眼半徑約為0.1 m.為了避免震源激發(fā)時(shí)地面的劇烈振蕩破壞檢波器和近井的井壁面波的影響，檢波器一般放置在離井一段距離之外，放置在r=1000 m,z=0 m處，其他參數(shù)也同4.1節(jié).

我們首先檢驗(yàn)“小井孔”假設(shè)失效時(shí)的情況.在井孔絕對(duì)半徑不變的情況下，通過減小特征波長(zhǎng)同樣可以使“小井孔”假設(shè)失效.在地層速度不變的前提下，提高震源的主頻可以減小特征波長(zhǎng).計(jì)算結(jié)果如圖5所示.該模型中，“小井孔”的條件為特征波場(chǎng)λ?a=0.1 m，圖5a中可以看出，該假設(shè)條件完全滿

圖4 徑向應(yīng)力源激發(fā)波場(chǎng)的半解析解和解析解對(duì)比(a)震源時(shí)間函數(shù)(Ricker子波)；(b)解析解遠(yuǎn)場(chǎng)波形(Ricker子波的導(dǎo)數(shù))；(c)實(shí)線表示最速下降積分得到ur的半解析解，虛線為解析解.倒三角形為P波的到時(shí).Fig.4 Comparison of the semi-analytical solution and the analytical solution of wave field excited by a radial stress source (a) Source time function (Ricker wavelet); (b) Far-field waveform of analytical solution (derivative of Ricker wavelet); (c) Solid line denotes the semi-analytical solution of ur obtained by Steepest Descent Integration Method and dashed line is the analytical solution. The inverted triangle denotes the theoretical arrival time of P waves.

圖5 高頻情況下徑向應(yīng)力源激發(fā)波場(chǎng)的半解析解(實(shí)線)和近似條件下的解析解(虛線)對(duì)比圖(a—f)分別是Ricker子波主頻為30, 50, 100, 300, 500, 1000 Hz時(shí)，位于r=1000 m,z=0 m處的檢波器記錄的徑向分量.倒三角形為P波理論到時(shí).Fig.5 Comparison of the semi-analytical solution (solid line) and the analytical solution (dashed line) of wave field excited by a radial stress source for high frequency Panels (a—f) are radial components of seismograms recorded by the receiver placed at r=1000 m,z=0 m with the peak frequency of Ricker wavelet being 30, 50, 100, 300, 500, 1000 Hz respectively. The inverted triangle denotes the theoretical arrival time of P waves.

足時(shí)，兩者之間差異很小.當(dāng)震源主頻逐漸增加時(shí)，子波波長(zhǎng)逐漸減小，圖5b和圖5c的假設(shè)條件接近失效，兩者之間的差異增加，但波形差異仍然較小.當(dāng)子波波長(zhǎng)進(jìn)一步減小時(shí)，“小井孔”假設(shè)失效(圖5d—5f)，雖然數(shù)值解和近似條件下的解析解的P波到時(shí)一直保持在理論到時(shí)0.4821 s處，但兩者在振幅和相位方面的差異開始凸顯.其中，半解析解的振幅略高于解析解的振幅，最大振動(dòng)的到時(shí)相對(duì)滯后.

解析解將波場(chǎng)對(duì)圓柱形井壁這一特殊的邊界條件的復(fù)雜響應(yīng)簡(jiǎn)化為一個(gè)尺度因子a2(見式(14)).滿足“小井孔”假設(shè)時(shí)，由于井孔半徑遠(yuǎn)小于特征波長(zhǎng)，鉆井對(duì)遠(yuǎn)場(chǎng)波場(chǎng)影響不大，解析解的簡(jiǎn)化較為合理.而當(dāng)“小井孔”假設(shè)失效時(shí)，井孔對(duì)波場(chǎng)的作用開始體現(xiàn)，即使在遠(yuǎn)場(chǎng)的地震記錄上也會(huì)有響應(yīng).

除了振幅和相位相對(duì)的差異外，兩者的絕對(duì)振幅都會(huì)隨著頻率的變大而迅速地增加.這是因?yàn)閮烧叩奈灰贫即笾抡扔赗icker子波的導(dǎo)數(shù).

4.2.2 低頻情況

當(dāng)震源頻率變低時(shí)，特征波長(zhǎng)會(huì)隨之增加.雖然這時(shí)“小井孔”假設(shè)不會(huì)受到影響，但是有可能會(huì)導(dǎo)致“遠(yuǎn)場(chǎng)”假設(shè)失效.低頻情況下的計(jì)算結(jié)果如圖6所示.該模型中“遠(yuǎn)場(chǎng)”假設(shè)的條件是λ?r=1000 m.從圖6a中可以看出，該假設(shè)完全滿足，兩者差異很小.但是當(dāng)震源主頻逐漸到5 Hz時(shí)，“遠(yuǎn)場(chǎng)”假設(shè)接近失效，兩者差異已經(jīng)開始體現(xiàn)(圖6b).隨著震源主頻的進(jìn)一步增加，兩者之間的差異越來越明顯(圖6c—6f).在低頻情況下，半解析解和解析解同樣出現(xiàn)了振幅和相位上的偏差.數(shù)值解的振幅大于解析解，并且比高頻情況下更加明顯(圖6d—6f).與高頻情況不同的是，低頻情況下數(shù)值解的波形與解析解也有比較明顯的差別，逐漸從Ricker子波的導(dǎo)數(shù)變成了Ricker子波.

一般情況下，遠(yuǎn)場(chǎng)是指炮檢距的尺度遠(yuǎn)大于震源尺度；而在井中震源問題中，“遠(yuǎn)場(chǎng)”假設(shè)比較的是炮檢距的尺度和特征波長(zhǎng)的尺度.所以雖然本例中炮檢距(1000 m)遠(yuǎn)大于震源的尺度(0.1 m)，但是并不一定滿足“遠(yuǎn)場(chǎng)”假設(shè).如圖6c—6f中，特征波長(zhǎng)大于或等于炮檢距時(shí)，解析解與數(shù)值解出現(xiàn)明顯的差異.4.3 軸向應(yīng)力震源

軸向應(yīng)力震源會(huì)用來模擬一些附著在井壁上的

圖6 低頻情況下徑向應(yīng)力源激發(fā)波場(chǎng)的半解析解(實(shí)線)和解析解(虛線)的對(duì)比圖(a—f)分別是Ricker子波主頻為30, 5, 2, 1, 0.5, 0.1 Hz時(shí)檢波器記錄的徑向分量.倒三角形為P波理論到時(shí).Fig.6 Comparison of the semi-analytical solution (solid line) and the analytical solution (dashed line) of wave field excited by a radial stress source for low frequency Panels (a—f) are radial components of seismograms recorded by the receiver placed at r=1000 m,z=0 m with the peak frequency of Ricker wavelet being 30, 5, 2, 1, 0.5, 0.1 Hz respectively. The inverted triangle denotes the theoretical arrival time of P waves.

圖7 軸向應(yīng)力源激發(fā)波場(chǎng)的半解析解(實(shí)線)和解析解(虛線)的對(duì)比圖(a—c)分別是Ricker子波主頻為30, 0.5, 500 Hz時(shí)檢波器的垂向分量.倒三角形是SV波的理論到時(shí).Fig.7 Comparison of the semi-analytical solution (solid line) and the analytical solution (dashed line) of wave field excited by a axial stress source Panels (a)—(c) are vertical components of seismograms with the peak frequency of Ricker wavelet being 30, 0.5, 500 Hz respectively. The inverted triangle denotes the theoretical arrival time of SV waves.

圖8 旋轉(zhuǎn)應(yīng)力源激發(fā)波場(chǎng)半解析解(實(shí)線)和解析解(虛線)的對(duì)比圖(a—c)分別是Ricker子波主頻為30, 0.5, 500 Hz時(shí)檢波器的橫向分量.倒三角形是SH波的理論到時(shí).Fig.8 Comparison of the semi-analytical solution (solid line) and the analytical solution (dashed line) of wave field excited by a torsional stress source Panels (a—c) are transverse components of seismograms with the peak frequency of Ricker wavelet being 30, 0.5, 500 Hz respectively. The inverted triangle denotes the theoretical arrival time of SH waves.

振蕩器震源、隨鉆震源，或者與徑向震源來共同模擬一些較為復(fù)雜的震源.計(jì)算結(jié)果如圖7所示.因?yàn)檩S向震源在檢波器處的位移場(chǎng)只有垂向分量，所以圖中對(duì)比的是解析解和數(shù)值解的垂向分量.軸向震源的遠(yuǎn)場(chǎng)波形是Ricker子波(圖7)，因此兩者絕對(duì)振幅都無明顯變化.當(dāng)滿足“小井孔”假設(shè)和“遠(yuǎn)場(chǎng)”假設(shè)時(shí)，兩者差異很小(圖7a)，而當(dāng)“遠(yuǎn)場(chǎng)”假設(shè)失效(圖7b)或“小井孔”假設(shè)失效(圖7c)時(shí)，兩者之間會(huì)有較為明顯的差異.

4.4 旋轉(zhuǎn)應(yīng)力震源

旋轉(zhuǎn)應(yīng)力震源可能由鉆頭旋轉(zhuǎn)產(chǎn)生，也有可能由鉆纜與井壁的摩擦產(chǎn)生.通常情況下旋轉(zhuǎn)應(yīng)力源的強(qiáng)度較小，產(chǎn)生的地震波能量也較弱.計(jì)算結(jié)果如圖8所示.與前兩種震源產(chǎn)生P-SV波不同，旋轉(zhuǎn)應(yīng)力震源只產(chǎn)生SH波，所以圖8中比較的是位移的橫向分量.由于旋轉(zhuǎn)震源遠(yuǎn)場(chǎng)波形是Ricker子波的導(dǎo)數(shù)，所以兩者的絕對(duì)振幅與頻率呈現(xiàn)正相關(guān).在“遠(yuǎn)場(chǎng)”假設(shè)(圖8b)或“小井孔”假設(shè)(圖8c)失效時(shí)，兩者同樣會(huì)出現(xiàn)明顯的差異.

5 結(jié)論

本文采用沿最速下降路徑的數(shù)值積分計(jì)算了井中震源的遠(yuǎn)場(chǎng)波場(chǎng).由于沿該路徑的積分函數(shù)不存在任何振蕩，因此獲得了高精度的數(shù)值積分.最速下降積分在計(jì)算時(shí)可以將P波和SV波分離，避免了波形之間的互相干擾.該方法保留了對(duì)波數(shù)域極點(diǎn)和支線進(jìn)行分析的可能性，可以解析地得到面波，折射波出現(xiàn)的位置和波形.模型試驗(yàn)表明，當(dāng)“小井孔”和“遠(yuǎn)場(chǎng)”假設(shè)不滿足時(shí)，近似條件下解析解的振幅和波形相對(duì)于半解析解都會(huì)有明顯的偏差，使用半解析解能獲得更準(zhǔn)確的波場(chǎng)信息.

附錄A 井中震源的頻率波數(shù)域解

柱坐標(biāo)系下的波動(dòng)方程為

(A1)

(A2)

(A3)

(A4)

其中

(A5)

是柱面波，也可以看作是頻率波數(shù)域的解.將(A4)、(A5)代入(A3)式中，并考慮到Hankel函數(shù)的定義

(A6)

可以得到

(A7)

(A7)式是關(guān)于f1,f2,f3的線性方程組，如要得到非零解，必須滿足

(A8)

(A9)

可以看到，P，SV，SH波在各向同性介質(zhì)中都是解耦的.從(A9)式可以得到q的三對(duì)共軛解，根據(jù)Sommerfeld輻射條件，我們選擇Imq>0的三個(gè)解來保證波場(chǎng)解在無窮遠(yuǎn)處的能量有限.當(dāng)Imq=0時(shí)，我們選擇ωReq>0的解來保證波場(chǎng)是向外傳播的.所以，勢(shì)函數(shù)在頻率波數(shù)域的解為

(A10)

(θ,z,t).

(A11)

將邊界條件也變換到頻率波數(shù)域，可以得到關(guān)于f1,2,3的線性方程組

(A12)

(A13)

(A14)

從(A12)、(A14)可以看出，任意類型的井中震源都可以產(chǎn)生P，SV和SH波.

軸對(duì)稱情況對(duì)應(yīng)于n=0的情況.這時(shí)Dij的表達(dá)式退化成

(A15)

解這個(gè)方程組得到f1,2,3后，將(A10)代入勢(shì)函數(shù)的表達(dá)式，得到頻率波數(shù)域的位移解表達(dá)式：

(A16)

記(A16)式中的矩陣為K，則

(A17)

其中detD是D的行列式，adjD是D的伴隨矩陣.

Abo-Zena A M. 1977. Radiation from a finite cylindrical explosive source.Geophysics, 42(7): 1384-1393, doi:10.1190/1.1440799.Aki K, Richards P G. 2002. Quantitative seismology. Sausalito, California: University Science Books.

Blair D P. 2007. A comparison of Heelan and exact solutions for seismic radiation from a short cylindrical charge.Geophysics, 72(2): E33-E41, doi:10.1190/1.2424543.

Blair D P. 2010. Seismic radiation from an explosive column.Geophysics, 75(1): E55-E65, doi:10.1190/1.3294860.

Bouchon M, Aki K. 1977. Discrete wave-number representation of seismic-source wave fields.Bull.Seismol.Soc.Am., 67(2): 259-277.Bouchon M. 1978. The importance of the surface or interface P wave in near-earthquake studies.Bull.Seismol.Soc.Am., 68(5): 1293-1311.Bouchon M. 2003. A review of the discrete wavenumber method.PureAppl.Geophys., 160(3-4): 445-465, doi:10.1007/PL00012545.

Chen S T, Eriksen E A, Miller M A. 1990. Experimental studies on downhole seismic sources.Geophysics, 55(12): 1645-1651, doi:10.1190/1.1442818.

Cheng C H, Toks?z M N. 1981. Elastic wave propagation in a fluid-filled borehole and synthetic acoustic Logs.Geophysics, 46(7): 1042-1053, doi:10.1190/1.1441242.

Cheng N Y. 1994. Borehole wave propagation in isotropic and anisotropic media: three-dimensional finite difference approach [Ph. D. thesis]. Massachusetts: Massachusetts Institute of Technology.

De Hoop A T, De Hon B P, Kurkjian A L. 1994. Calculation of transient tube-wave signals in cross-borehole acoustics.J.Acoust.Soc.Am., 95(4): 1773-1789, doi:10.1121/1.408697.Dong W J, Toks?z M N. 1995. Borehole seismic-source radiation-pattern in transversely isotropic media.Geophysics, 60(1): 29-42, doi:10.1190/1.1443759.

Haldorsen J B U, Miller D E, Walsh J J. 1995. Walk-away VSP using drill noise as a source.Geophysics, 60(4): 978-997, doi:10.1190/1.1443863.

Heelan P A. 1953. Radiation from a cylindrical source of finite length.Geophysics, 18(3): 685-696, doi:10.1190/1.1437923.Huang W C, Ge H K, Wang B S, et al. 2010. Deconvolution interferometry and its application to seismic while drilling data processing.ProgressinGeophysics(in Chinese), 25(3): 951-956, doi:10.3969/j.issn.1004-2903.2010.03.032.

Jordan D W. 1962. The stress wave from a finite, cylindrical explosive source.J.Math.Mech., 11(4): 503-551.

Lapwood E R. 1949. The disturbance due to a line source in a semi-infinite elastic medium.Philos.Trans.Roy.Soc.Lond.AMath.Phys.Sci., 242(841): 63-100, doi:10.1098/rsta.1949.0005.Lee M W, Balch A H. 1982. Theoretical seismic-wave radiation from a fluid-filled borehole.Geophysics, 47(9): 1308-1314, doi:10.1190/1.1441391.

Lee M W. 1986. Low-frequency radiation from point sources in a fluid-filled borehole.Geophysics, 51(9): 1801-1807, doi:10.1190/1.1442226.

Liu Y B, Li Y M, Wu R S, et al. 1993a. Far-field radiation pattern of dowhole source and array source.OGP(in Chinese), 28(4): 379-388.

Liu Y B, Li Y M, Wu R S, et al. 1993b. Radiative energies from downhole source and array source.OGP(in Chinese), 28(4): 389-395.

Lu B, Ge H K, Wu H Z, et al. 2009. SWD data preprocessing using wavelet transform of correlation domain.ChineseJournalGeophysics(in Chinese), 52(9): 2349-2356, doi:10.3969/j.issn.0001-5733.2009.09.020 Meredith J A. 1991. Numerical and analytical modelling of downhole seismic sources, the near and far field [Ph.D.thesis]. Massachusetts: Massachusetts Institute of Technology.Meredith J A, Toks?z M N, Cheng C H. 1993. Secondary shear waves from source boreholes.Geophys.Prospect., 41(3): 287-312, doi:10.1111/j.1365-2478.1993.tb00571.x.

Poletto F. 2005. Energy balance of a drill-bit seismic source, part 1: Rotary energy and radiation properties.Geophysics, 70(2): T13-T28, doi:10.1190/1.1897038.

Rector J W, Marion B P. 1991. The use of drill-bit energy as a downhole seismic source.Geophysics, 56(5): 628-634, doi:10.1190/1.1443079.

Rector J W, Hardage B A. 1992. Radiation pattern and seismic waves generated by a working roller-cone drill bit.Geophysics, 57(10): 1319-1333, doi:10.1190/1.1443199.

Shen J G, Zhang H L. 2000. Numerical study on 3D acoustic field generated by eccentric sources in borehole.ChineseJournalGeophysics(in Chinese), 43(2): 279-286.

Tsang L, Rader D. 1979. Numerical evaluation of the transient acoustic waveform due to a point source in a fluid-filled borehole.Geophysics, 44(10): 1706-1720, doi:10.1190/1.1440932.

Tubman K M, Cheng C H, Toks?z M N. 1984. Synthetic full waveform acoustic logs in cased boreholes.Geophysics, 49(7): 1051-1059, doi:10.1190/1.1441720.

Vasconcelos I, Snieder R. 2008. Interferometry by deconvolution: Part 2-Theory for elastic waves and application to drill-bit seismic imaging.Geophysics, 73(3): S129-S141, doi:10.1190/1.2904985.

Wang P, Ge H K, Lu B, et al. 2009. Numerical simulation on SWD seismic wave propagation and data processing.PetroleumDrillingTechniques(in Chinese), 37(2): 5-9.

Wu H Z, Ge H K, Yang D H, et al. 2010. A research of cepstrum analysis of drill string vibration and extracting the bit source signals.ChineseJournalGeophysics(in Chinese), 53(8): 1968-1975, doi:10.3969/j.issn.0001-5733.2010.08.023.

Zhang B, Li Y M, Liu Y B. 1995. The calculation of cross-well seismic in isotropic porous layered media.ActaGeophysicaSinica(in Chinese), 38(4): 507-518.

附中文參考文獻(xiàn)

黃偉傳, 葛洪魁, 王寶善等. 2010. 反褶積干涉成像及其在隨鉆地震數(shù)據(jù)處理中的應(yīng)用. 地球物理學(xué)進(jìn)展, 25(3): 951-956, doi:10.3969/j.issn.1004-2903.2010.03.032.

劉銀斌, 李幼銘, 吳如山等. 1993a. 井下震源和陣列震源的遠(yuǎn)場(chǎng)輻射花樣. 石油地球物理勘探, 28(4): 379-388.

劉銀斌, 李幼銘, 吳如山等. 1993b. 井下震源和陣列震源的輻射能量. 石油地球物理勘探, 28(4): 389-395.

陸斌, 葛洪魁, 吳何珍等. 2009. 利用相關(guān)域小波變換進(jìn)行SWD資料預(yù)處理. 地球物理學(xué)報(bào), 52(9): 2349-2356, doi:10.3969/j.issn.0001-5733.2009.09.020.

沈建國, 張海瀾. 2000. 井內(nèi)偏心聲源激發(fā)的三維聲場(chǎng)的數(shù)值研究. 地球物理學(xué)報(bào), 43(2): 279-286.

王鵬, 葛洪魁, 陸斌等. 2009. 隨鉆地震波場(chǎng)傳播與數(shù)據(jù)處理方法的數(shù)值實(shí)驗(yàn). 石油鉆探技術(shù), 37(2): 5-9.

吳何珍, 葛洪魁, 楊頂輝等. 2010. 鉆柱振動(dòng)倒譜分析及其鉆頭源信號(hào)提取方法研究. 地球物理學(xué)報(bào), 53(8): 1968-1975, doi:10.3969/j.issn.0001-5733.2010.08.023.

張釙, 李幼銘, 劉銀斌. 1995. 層狀各向同性多孔介質(zhì)中井間地震的數(shù)值計(jì)算. 地球物理學(xué)報(bào), 38(4): 507-518.

(本文編輯 胡素芳)

1StateKeyLaboratoryofLithosphericEvolution,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China3CASCenterforExcellenceinTibetanPlateauEarthSciences,Beijing100101,China

Borehole sources, whose scope goes far beyond sources in boreholes, are of extreme importance in research with active seismic sources, including deep seismic sounding, reverse vertical seismic profiling (RVSP), seismic while drilling, mining geophysics, etc. Sources used in these studies are all of cylindrical structures, which is the reason why they are called borehole sources and why their wave fields has unique characteristics. Previous studies on borehole sources are mostly based on analytical solutions obtained when small-borehole assumption (the borehole radius is significantly smaller than the characteristic wave length) and far-field assumption (the offset is greater than the characteristic wave length) are satisfied. It is still an open question whether the analytical solutions are applicable to cases that violate the two assumptions.This study is based on the synthetic seismograms computed by both analytical solutions and semi-analytical solutions. The analytical solutions used in previous studies are obtained through asymptotic analysis, while the semi-analytical solutions are computed by numerical integration. The semi-analytical solutions are of higher accuracy and therefore regarded as “true solutions”. Synthetic seismograms from the analytical solutions are compared to true solutions to validate whether the analytical solutions are applicable to certain cases or not. Accuracy is crucial to the comparison. Yet the high oscillation of solutions in frequency-wavenumber domain brings out a great challenge. We developed a brand-new numerical method called Steepest Descent Integration Method (SDIM). The new method is inspired by the Method of Steepest Descent (SDM) in asymptotic analysis that is specially designed for highly oscillatory integral and is the very method used to obtain the analytical solutions. Replacing approximate integration path and approximate integrand in SDM with accurate ones, SDIM breaks the restraints of small borehole and far field and can compute seismograms at arbitrary offset and arbitrary source frequency with extremely high accuracy efficiently. We calculate the seismograms by both SDIM and SDM for a large offset (1000 m, significantly large compared to borehole radius of 0.1 m) and varied source frequency (0.1～1000 Hz). The assumption of small-borehole is violated in high frequency cases, while far-field assumption fails when the frequency is low. The same experiment is conducted for all three basic borehole sources.The works presented in the paper can be categorized into two parts, namely the new SDIM and comparison of seismograms. The study of SDIM shows that: (1) The solutions of borehole sources problem in frequency-wavenumber domain are highly oscillatory. The oscillation depends on source frequency and offsets. High frequency sources result in severe oscillation, so as large offsets. (2) The oscillation is attributed to Hankel functions in the solutions whose exponential part account for most of it. Hence, exponential functions are used in the derivation of SDIM instead of Hankel functions, making the work much easier. (3) The only difference between SDIM and SDM is the accuracy of the steepest descent path and the integrand. SDIM uses the accurate path and integrand, while SDM uses approximate ones. In addition, four numerical examples are presented in the paper. Each is designed specifically. They demonstrate that: (1) Results from SDIM are identical to ones from SDM when small-borehole assumption and far-field assumption are satisfied, which supports the validity of SDIM. (2) When small-borehole assumption is violated, the SDM results differ much from the SDIM ones that are considered as true results. It infers that the influence from borehole might not be ignored even for far-field wave field. (3) When far-field assumption fails, the results from SDM are inaccurate as well, which means the absolute value of the offset cannot guarantee far-field. The ratio of the offset to the characteristic wave length matters. (4) The same phenomenon occurs in the wave field of all the three basic borehole sources.Obtaining accurate far-field seismograms is the key problem of borehole sources research. Yet it is challenging because of highly oscillatory integral involved. By taking advantage of the special form of analytical solutions, we developed a brand-new method for computing highly oscillatory wavenumber integration. It completely avoids the oscillation and results in numerical integration of a fully smooth function, leading to synthetic seismograms with high precision. It also allows us to compute P, S and surface waves separately, reducing their mutual interference. Numerical experiments demonstrate that the results from SDM are considerably different from ones from SDIM, in both amplitude and phase when the small-borehole assumption or the far-field assumption fails. Therefore, the SDM has its constraint and SDIM is a better alternative if accurate wave-field information is needed.

Downhole seismic sources; Far-field wave field; Analytical solution; Steepest Descent Integration Method; Method of steepest descent

中國地震局公益性行業(yè)科研專項(xiàng)(20140823)和國家自然科學(xué)基金(41174075，41274070，41374062,41474068)聯(lián)合資助.

徐逸鶴，男，1990年生，博士生，主要從事井中震源和隨鉆震源波場(chǎng)研究. E-mail: xuyihe@mail.iggcas.ac.cn

*通訊作者徐濤，男，1978年生，副研究員，主要從事地震射線理論與殼幔結(jié)構(gòu)成像研究.E-mail: xutao@mail.iggcas.ac.cn

10.6038/cjg20150824.

10.6038/cjg20150824

P631

2015-02-09，2015-06-23收修定稿

徐逸鶴，徐濤，王敏玲等. 2015. 井中震源的遠(yuǎn)場(chǎng)波場(chǎng)特征研究.地球物理學(xué)報(bào)，58(8):2912-2926,

Xu Y H, Xu T, Wang M L, et al. 2015. Far-field wavefield characteristics of downhole seismic sources.ChineseJ.Geophys. (in Chinese),58(8):2912-2926,doi:10.6038/cjg20150824.