馬明明，陳 浩，何 曉，王秀明

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

1 引 言

地層橫波慢度信息對石油勘探開發(fā)有著重要的意義，它與地層巖石性質(zhì)密切相關(guān).井壁周圍地層性質(zhì)常常因為井壁應(yīng)力集中、鉆井機械損壞和泥漿入侵等外力因素發(fā)生變化，進而導(dǎo)致井壁周圍地層出現(xiàn)橫波慢度的徑向變化帶.反之，如果已知徑向變化帶橫波慢度分布情況，就可以了解井壁周圍地層性質(zhì)的改變.對于慢地層（地層橫波慢度大于井孔流體縱波慢度），傳統(tǒng)的單極子儀器探測不到橫波首波，原因是井外地層的橫波慢度大于井孔流體縱波慢度時，全波列中無法探測到橫波首波信號.因此，人們通常采用偶極彎曲波來反演地層橫波的信息［1－2］.

1996—2006年Sinha等發(fā)展了基于微擾法［3］和BG理論［4］的反演算法［5－8］，來反演地層橫波慢度的徑向分布，在處理實際數(shù)據(jù)時，此方法不能同時滿足高分辨率和高精度需求.2000年，Braunisch等又提出一種利用多頻、多模式波的信息聯(lián)合反演地層橫波慢度的方法［9］，但此方法在慢地層應(yīng)用受限.2010年，Tang等建立了參量化的、高頻約束條件下的反演方程［10］，方程假定地層橫波慢度是徑向半徑r的指數(shù)函數(shù)，并且在反演的過程中對彎曲波進行高頻約束，但高頻彎曲波的激發(fā)強度較小，易受噪聲等的干擾，為反演結(jié)果的準(zhǔn)確性帶來了挑戰(zhàn).2011—2012年，Yang等建立了套管井模型下的參量化的反演方程［11－12］，他們也假定地層橫波慢度是徑向半徑r的指數(shù)函數(shù)，指數(shù)函數(shù)的變化形式由3個參數(shù)控制，反演地層橫波慢度的徑向分布，轉(zhuǎn)換為對這3個參數(shù)的反演，但是反演計算方法的收斂性受初始值影響嚴(yán)重.

另外，國內(nèi)的陳雪蓮和王瑞甲研究了TI地層模式波對井孔流體參數(shù)和地層參數(shù)的靈敏度，分析了各模式波的徑向探測深度，其中彎曲波探測深度約為4個井孔半徑［13］.這也說明可通過彎曲波來獲取井壁周圍地層的橫波慢度，進而了解地層的性質(zhì).

前人的工作都說明了，利用彎曲波頻散可以求解井壁周圍地層橫波慢度的徑向分布，進而了解井壁周圍的地層性質(zhì).但是，前人對所采用的反演方法的特性分析較少，并很少涉及如何提高計算效率.本文的工作則主要側(cè)重于這兩個方面.本文在前人的工作基礎(chǔ)上，首先計算了地層橫波慢度徑向均勻（以下簡稱均勻模型）和非均勻情況下（以下簡稱非均勻模型，非均勻模型的產(chǎn)生是由于均勻模型的橫波慢度在徑向上發(fā)生了變化）的時域全波列波形.其次利用頻域方法［14］提取了波形的彎曲波頻散曲線，并對彎曲波頻散曲線進行了擬合校正［15］.再次，討論了如何結(jié)合微擾法與BG理論，基于非均勻模型和均勻模型之間彎曲波的差異，進行地層橫波慢度徑向分布的反演計算.最后，我們分析了如何在反演的迭代過程中提高計算效率，并且用多種計算實例，驗證了所采用的反演算法的魯棒性，為它的現(xiàn)場應(yīng)用打下基礎(chǔ).

2 偶極井孔聲場理論

在三維柱坐標(biāo)系（r，θ，z）下，以偶極源的中心為原點，并使z軸與井軸方向一致.根據(jù)文獻［16］，井孔內(nèi)位移勢函數(shù)φf在頻率波數(shù)內(nèi)的穩(wěn)態(tài)聲場表達式如下：

井外地層聲場縱波位移勢函數(shù)，SH橫波位移勢函數(shù)和SV橫波位移勢函數(shù)在頻率波數(shù)內(nèi)的穩(wěn)態(tài)聲場表達式如式：

其中，n＝1代表偶極源，ra為偶極距，r為徑向半徑，εn＝2，ω為角頻率，In和Kn為n階虛宗量貝塞爾函數(shù)，A′n、An、Bn、Cn、Dn、En和Fn是待定系數(shù)，由井壁處的邊界條件確定分別為井孔流體縱波、井外地層縱波和井外地層橫波的徑向波數(shù)，sf、sp和ss分別為井孔流體縱波、井外地層縱波和井外地層橫波慢度，k是軸向波數(shù).

在相鄰2個層的界面上聲場應(yīng)滿足邊界條件.如果在半徑rj（下標(biāo)j＝0，…，N表示徑向上的層的序號，N表示徑向上的層數(shù)）處的界面兩側(cè)都是固體，在界面上膠結(jié)良好，則界面上的位移和應(yīng)力分量應(yīng)該連續(xù)，即

式中，ur為法向位移分量，uθ和uz為切向位移分量，σrr為法向應(yīng)力分量，σrz和σrθ為切向應(yīng)力分量.

在井壁r0處，其界面的內(nèi)側(cè)是液體，外側(cè)是固體，此時邊界條件變成法向位移ur和法向應(yīng)力σrr連續(xù)，切向應(yīng)力σrz和σrθ為零，切向位移分量uθ和uz無需考慮，即

層與層之間的聲場由湯姆森—哈斯克（Thomson－Hanskell）傳播矩陣［17］來連接，反復(fù)利用各層介質(zhì)的傳遞矩陣和邊界條件，就可以得求的徑向分層地層的聲場.對于偶極聲源，測量的是徑向位移，是瞬態(tài)聲場，需要考慮聲源的影響，其在空間－時間域的具體表達式為［16］

2.1 時域波形

計算分析時，本文給出了兩個模型，模型1的井外是無限大均勻地層，模型2的井外地層橫波慢度徑向分層，它是由模型1的地層橫波慢度徑向擾動所致.兩個模型的井孔半徑r0＝0.1m，井孔流體密度ρf＝1000kg／m3，井孔流體縱波慢度sf＝666.7μs／m.模型1的地層參數(shù)見表1，模型2中，地層橫波慢度徑向是非均勻的，導(dǎo)致剪切模型模量徑向非均勻，進而引起縱波慢度徑向非均勻，然而縱波慢度對彎曲波頻散影響很?。?］，所以為了提高計算效率，計算過程中假設(shè)模型2的地層縱波慢度與模型1相同，此外，地層密度對彎曲波頻散影響也很小，在計算時也同樣假設(shè)模型2的地層密度與模型1相同，模型2中地層橫波慢度ss的徑向分布數(shù)值見表2.圖1為模型2的示意圖，其中徑向半徑b和井壁r0之間為橫波慢度的變化帶，b＝0.3m.

表1 地層參數(shù)Table 1 Parameters of formation

表2 變化帶橫波慢度值Table 2 The values of shear wave slowness of altered zone

圖1 井孔與地層模型Fig.1 The model of borehole and formation

文中計算用到的聲源函數(shù)為F（ω）＝聲源中心頻率f0＝3kHz，聲源帶寬fb＝0.4f0，聲源與第一個接收器之間的距離d＝3.0m，共八個接收器，相鄰接收器之間的距離dz＝0.1524m.

圖2中黑色和紅色曲線分別表示模型1和模型2的偶極全波波列，從圖2可以看出模型2的全波波列明顯延后于模型1的全波波列，這是由于模型2的井壁周圍地層橫波慢度增加所致.

2.2 頻散分析

在圖2的全波列中，占主要地位的是彎曲波成分.為了分析彎曲波的頻散特征，可以采用頻域法從時域波形中提取各頻率下彎曲波的慢度.本文利用頻域加權(quán)相似法對圖2中兩個模型的全波波列進行處理，提取了偶極全波波列的彎曲波.

圖3a中黑色和紅色曲線分別表示模型1和模型2中第一道波列的頻譜.圖3b中，黑色和紅色圓圈狀曲線分別為對模型1和模型2的波列提取的彎曲波頻散，對它們利用正切（或反正切）函數(shù)進行擬合，分別得到圖中黑色和紅色實線，利用頻散方程［19］計算兩個模型的彎曲波頻散見三角狀曲線，這三種方法求解的頻散吻合性較好.擬合的好處是可以得到寬頻帶范圍內(nèi)的彎曲波，有效地消除了頻散曲線中異常數(shù)據(jù)點對反演結(jié)果的影響，提高了反演計算的準(zhǔn)確性.此外，從圖3b中可以看出頻率在1kHz以下時，紅色和黑色實線幾乎重合，并且隨著頻率向低頻移動，二者的彎曲波慢度都趨向于模型1中地層橫波慢度，這也就幫助確定發(fā)生徑向變化前地層橫波慢度的數(shù)值，在反演計算過程中，此數(shù)值為參考模型的地層橫波慢度.

圖2 偶極子全波列波形圖，模型1（黑色）和模型2（紅色）Fig.2 Dipole waveforms of model one（black）and model two（red）

3 反演方法

由圖2和圖3可以看出，地層橫波慢度在徑向上的變化分布影響了偶極子彎曲波傳播特征，因此彎曲波傳播特征的改變可以被用來確定地層橫波慢度在徑向上的變化情況，也即地層性質(zhì)的改變.

本文選擇的反演方法結(jié)合了微擾法和BG理論，根據(jù)此反演方法，待求模型的井外有一個橫波慢度的徑向變化帶，此變化帶是由于參考模型（參考模型的井外是無限大均勻地層）的橫波慢度在近井壁處發(fā)生徑向擾動所致.處理待求模型的偶極全波波列得到的彎曲波慢度在低頻處的趨近值，可以認(rèn)為是參考模型的地層橫波慢度srefer.利用牛頓－拉夫森（Newton－Raphson）方法［20］計算頻散方程：

可以求解出參考模型的彎曲波頻散vrefer＝ω／k，進而求解出擾動部分貢獻的彎曲波頻散?v＝vrealvrefer，其中vreal表示提取的彎曲波頻散，符號?及其后面符號代表小量的意思.利用?v，得到地層橫波慢度發(fā)生的擾動?ss，它與參考模型的橫波慢度之和，即為實際模型的橫波慢度徑向分布.

3.1 微擾法

在柱坐標(biāo)系下，根據(jù)微擾法，當(dāng)參考模型空間中有擾動發(fā)生后，擾動后的模型的彎曲波相速度也相應(yīng)的發(fā)生了改變：

圖3 （a）兩個模型的第一道波列的頻譜，（b）兩種模型的彎曲波頻散：提取（圓圈狀）、擬合（線狀）和利用頻散方程計算（三角形）Fig.3 （a）The first spectrum of each model；（b）Flexural dispersion curves of each model：the extractive ones（circle），the fitted ones（line）and the calculated ones（triangle）based on dispersion equation

因此只要已知沒有擾動時參考模型的相速度，就可求解發(fā)生擾動后的模型的相速度：

因此根據(jù)方程（9），當(dāng)?shù)貙訖M波慢度從井壁r0處到徑向半徑b處，發(fā)生了擾動可以計算地層橫波慢度從井壁r0處到b發(fā)生徑向擾動后的彎曲波頻散，只是這里的積分區(qū)域變成了r0到b，即

3.2 BG理論

Backus和Gilbert提出由于反問題解的非唯一性，不去求反問題解本身而去求它的某種平均乃是消除反問題非唯一性的一種途徑，這種平均稱為解估計.這里結(jié)合微擾法和BG理論，給出了利用彎曲波頻散求解地層橫波慢度的主要公式.

假設(shè)已知在頻率點ωi處（i＝1，…，N），彎曲波相速度的相對變化?vi／vi，是由r′處剪切模量發(fā)生的相對擾動?μ（r′）／μ造成的，即有：

根據(jù)BG理論，在反演?μ（r′）／μ時，不去求解?μ（r′）／μ本身，而是求解?μ（r′）／μ的估計?μ＾（r′）／μ具有更好的穩(wěn)定性.如下：

因此只要求解出待定系數(shù)組｛ai｝i＝1，…，N，就能計算出r′處的剪切模量的相對變化的估計.Backus和Gilbert的研究給出了系數(shù)組的一系列準(zhǔn)則和方法，具體的過程見參考文獻［5］.

下面介紹如何提高所采用的結(jié)合微擾法和BG理論的反演方法的計算效率，這也是本文的重點.

根據(jù)BG理論，正確評價解估計的準(zhǔn)則是在展布函數(shù)和方差之間取合理的折衷.Backus和Gilbert建議采用以下目標(biāo)函數(shù)：

其中，aTSa是平均核函數(shù)A（r，r′）的展布函數(shù)，aTEa是解估計的方差，β為折衷參數(shù)且β∈（0，π／2）.

求解方程（15）的過程是個對β不斷尋優(yōu)的過程，采用傳統(tǒng)的迭代方式，計算效率較低，本文采用了1995年楊文采提出的改進措施［21］.在迭代開始時，先求低分辨率小方差的解估計，隨非線性迭代逐步提高分辨率.具體地，在開始時取折衷因子β0＝π／2，隨迭代逐步減小.即：

迭代次數(shù)q增加時，βq減小，這樣的選取方式避免了內(nèi)循環(huán)中解線性方程組的無謂計算，大大降低了計算成本.

3.3 反演流程

結(jié)合微擾法與BG理論，利用彎曲波頻散反演地層橫波慢度徑向分布的基本流程可以歸納為：

（1）井孔半徑、井孔流體參數(shù)和地層參數(shù)（除橫波慢度外）已知，從波形數(shù)據(jù)中提取彎曲波頻散vreal，對其進行擬合得到了寬頻帶范圍的彎曲波慢度，其在低頻處趨近的數(shù)值即為參考模型的橫波慢度；

（2）對于選定的參考模型，計算偶極頻散方程，得到參考模型的彎曲波頻散vrefer；

（3）比較vreal與vrefer，若二者之間的差別小于某一小值min，則認(rèn)為地層橫波慢度沒有發(fā)生徑向擾動，否則認(rèn)為地層橫波慢度相對參考模型發(fā)生了徑向擾動，則開始反演計算：

首先計算核函數(shù)（10）式，其次聯(lián)合（10）—（15）式計算出待定系數(shù)組｛ai｝，i＝1，2，…，N，則得到了徑向半徑r′處的地層橫波慢度的解估計：

改變r′值就可得到井外地層橫波慢度的徑向分布，它與地層橫波慢度真實值之間相對誤差為：

4 計算實例

本文假設(shè)模型2的偶極全波列為實際測得數(shù)據(jù).為了分析結(jié)合微擾法和BG理論反演方法的特性，本文對偶極全波列存在以下5種情況時的數(shù)據(jù)進行了處理，得到了各情況下地層橫波慢度的徑向分布.這5種情況是：（1）無噪聲和參數(shù)誤差；（2）井孔流體慢度sf存在10%以內(nèi)誤差；（3）井外地層縱波慢度sp存在10%以內(nèi)誤差；（4）全波波列存在SNR＝20dB的噪聲；（5）全波波列存在SNR＝10dB的噪聲.

4.1 無噪聲和參數(shù)誤差

假設(shè)圖2中的紅色波列數(shù)據(jù)已知，其所對應(yīng)的模型井外地層橫波慢度徑向分布情況未知.根據(jù)第2章介紹的反演方法，得到的參考模型的彎曲波頻散見圖4a中黑色曲線，處理圖2中的紅色波列得到的未知模型的彎曲波頻散見圖4a中紅色曲線.本文選擇6個頻率點下兩個模型的彎曲波慢度之差作為反演計算的輸入?yún)?shù)，相鄰點之間頻率間隔為Δf＝1.0kHz，每個頻率點下兩個模型的彎曲波慢度用直線連接，見圖中圓圈所標(biāo)示.

圖4b中紅色曲線為反演計算得到的地層橫波慢度的徑向分布，黑色階梯狀曲線為地層橫波慢度真實的徑向分布.反演結(jié)果相對誤差的徑向分布見圖4c，從圖中可以看到反演結(jié)果相對誤差基本控制在10%以內(nèi)，并且離井壁越近，誤差相對越大.

4.2 其他參數(shù)存在誤差

實際數(shù)據(jù)比理論數(shù)據(jù)要更復(fù)雜，可能導(dǎo)致獲取的其他地層參數(shù)或者井孔參數(shù)不準(zhǔn)確，為此，我們還分析了這些參數(shù)存在誤差時對反演結(jié)果的影響.本文首先研究了井孔流體慢度sf分別存在－10%、－5%、5%和10%的誤差時對反演結(jié)果的影響.反演的與實際的地層橫波慢度的徑向分布見圖5a，反演結(jié)果的相對誤差的徑向分布見圖5b.我們注意到，sf存在誤差與sf無誤差的反演結(jié)果保持了很好的一致性.

其次研究了井外地層縱波慢度sp分別存在－10%、－5%、5%和10%的誤差時對反演結(jié)果的影響.反演的與實際的地層橫波慢度的徑向分布見圖6a，反演結(jié)果的相對誤差的徑向分布見圖6b.我們注意到，sp存在誤差與sp無誤差的反演結(jié)果也保持了很好的一致性，這也是計算時假設(shè)地層縱波慢度為一固定值的原因.

4.3 存在噪聲

實際數(shù)據(jù)中往往存在噪聲，為此還應(yīng)針對波列數(shù)據(jù)存在噪聲的情況做進一步的分析.這部分的工作是通過對圖2中的紅色波列分別添加SNR＝20dB和SNR＝10dB的高斯白噪聲各4次來展開的.

圖7顯示了分別添加了一次SNR＝20dB和SNR＝10dB的高斯白噪聲后的波形曲線和頻譜曲線.圖7a中黑色、紅色和藍色曲線分別表示：無噪聲時、存在SNR＝20dB和SNR＝10dB的噪聲時的偶極全波列波形；圖7b的曲線為相應(yīng)的第一道波形，圖7c為相應(yīng)的第一道波形的頻譜（圖中標(biāo)示了果的準(zhǔn)確性大大降低.擬合后的彎曲波頻帶范圍變寬，可利用的彎曲波數(shù)據(jù)增加，進而降低了噪聲的干擾.加入4次SNR＝20dB噪聲后反演的與實際的地層橫波慢度的徑向分布見圖8e，反演結(jié)果相對誤差的徑向分布見圖8f.加入4次SNR＝10dB噪聲后反演的與實際的地層橫波慢度的徑向分布見圖9e，反演結(jié)果相對誤差的徑向分布見圖9f.從圖8f和9f可見：存在SNR＝20dB和SNR＝10dB噪聲時，相對誤差也基本控制在10%以內(nèi)，并且與沒有噪聲情況下的誤差吻合度較高，這說明采用的反演方法對噪聲不敏感.

圖7 （a）偶極全波列波形；（b）第一道波列；（c）第一道波列的頻譜；（d）SNR＝20dB時八道波列的頻譜；（e）SNR＝10dB時八道波列的頻譜Fig.7 （a）Dipole full waveforms；（b）The first waveform of each case；（c）The first spectrum of each case；（d）Eight spectrums when SNR＝20dB；（e）Eight spectrums when SNR＝10dB

5 結(jié)論與討論

本文以井外地層橫波慢度徑向分層模型為例，計算了模型的偶極全波波列，提取了模型的彎曲波頻散，利用此彎曲波慢度與井外地層均勻時彎曲波慢度在多頻率點下的差異，采用結(jié)合微擾法和BG理論的反演方法，得到了徑向分層模型井外地層橫波慢度的徑向分布，并且討論了如何提高反演計算的效率，以適應(yīng)現(xiàn)場對數(shù)據(jù)處理速度的要求.

此外，為了分析采用的反演方法的特性，本文分別對以下5種情況的偶極全波波列數(shù)據(jù)進行處理，得到了各情況下地層橫波慢度的徑向分布.這5種情況是：（1）無噪聲和參數(shù)誤差時；（2）井孔流體慢度sf存在10%以內(nèi)誤差時；（3）井外地層縱波慢度sp存在10%以內(nèi)誤差時；（4）全波波列存在SNR＝20dB的噪聲時；（5）全波波列存在SNR＝10dB的噪聲時.結(jié)果表明：（2）和（3）的反演結(jié)果與（1）的反演結(jié)果保持了很好的一致性，這是因為sf和sp對彎曲波的影響較小的緣故；在（4）和（5）兩種情況下，通過對提取的彎曲波頻散進行擬合，得到了較寬頻帶范圍內(nèi)的彎曲波數(shù)據(jù)，有效地消除了頻散曲線中異常數(shù)據(jù)點對反演結(jié)果的影響，使得（4）和（5）的反演結(jié)果相對于（1）的反演結(jié)果都沒有發(fā)生明顯的改變，這說明采用的反演方法對噪聲不敏感.

圖9 （a—d）SNR＝10dB時提?。ㄈ睿┖蛿M合（線狀）的彎曲波頻散；（e）反演的與實際的地層橫波慢度的徑向分布；（f）反演的與實際的地層橫波慢度之間相對誤差的徑向分布Fig.9 （a—d）Extractive（circle）and fitted（line）flexural dispersion curves when SNR＝20dB；（e）Inverted and actual radial profiles of shear wave slowness；（f）Radial profile of relative error between inverted and actual shear wave slowness

以上反演結(jié)果都證實了，結(jié)合微擾法和BG理論的反演方法具有很好的魯棒性，可以被用于反演井壁周圍地層橫波慢度的徑向分布，進而了解井壁周圍地層的性質(zhì).

（

）

［1］ Schmitt D P.Shear wave logging in elastic formations.JournaloftheAcousticalSocietyofAmerica，1988，84（6）：2215－2229.

［2］ Winbow G A.A theoretical study of acoustic S－wave and Pwave velocity logging with conventional and dipole sources in soft formations.Geophysics，1988，53（10）：1334－1342.

［3］ Ellefsen K J，Cheng C H，Toks?z M N.Applications of perturbation theory to acoustic logging.JournalofGeophysical Research，1991，96（B1）：537－549.

［4］ Backus G，Gilbert F.Uniqueness in the inversion of inaccurate gross earth data.Philosophical Transactions of the Royal Society A：Mathematical.PhysicalandEngineering Sciences，1970，266（1173）：123－192.

［5］ Burridge R，Sinha B K.Inversion for formation shear modulus and radial depth of investigation using borehole flexural waves.SEG Expanded Abstracts，1996，158－161.

［6］ Sinha B K，Burridge R.Radial profiling of formation shear velocity from borehole flexural dispersions.IEEE UltrasonicsSymposiumProceedings，2001，391－396.

［7］ Sinha B K，Asvadurov S.Dispersion and radial depth of investigation of borehole modes.GeophysicalProspecting，2004，52（4）：271－286.

［8］ Sinha B K，Vissapragada B，Renlie L，et al.Radial profiling of the three formation shear moduli and its application to well completions.Geophysics，2006，71（6）：E65－E77.

［9］ Braunisch H，Habashy T M，Sinha B K，et al.Inversion of borehole dispersions for formation elastic moduli.IEEE UltrasonicsSymposium，2000，551－556.

［10］ Tang X M，Patterson D J.Mapping formation radial shearwave velocity variation by a constrained inversion of borehole flexural－wave dispersion data.Geophysics，2010，75（6）：E183－E190.

［11］ Yang J，Sinha B K，Habashy T M.A Parameterized－modelbased radial profiling for formation shear slowness in cased boreholes.SEG San Antonio 2011Annual Meeting，2011，449－453.

［12］ Yang J Q，Sinha B K，Habashy T M.A theoretical study on formation shear radial profiling in well－bonded cased boreholes using sonic dispersion data based on a parameterized perturbation model.Geophysics，2012，77（3）：WA 197－WA 210.

［13］ 陳雪蓮，王瑞甲.橫向各向同性彈性地層井孔中模式波的探測深度.吉林大學(xué)學(xué)報（地球科學(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.JournalofJilinUniversity（EarthScience Edition）（in Chinese），2008，38（3）：502－507.

［14］ Nolte B，Huang X J.Dispersion analysis of split flexural waves：Borehole acoustics and logging and reservoir delineation consortia annual report.Massachusetts Institute of Technology，1997.

［15］ Tang X M，Li C，Patterson D J.A curve－fitting technique for determining dispersion characteristics of guided elastic waves.Geophysics，2010，75（3）：E153－E160.

［16］ Paillet F L，Cheng C H.Acoustic waves in boreholes.Boca Raton／Ann Arbor／Boston London：CRC Press，1991.

［17］ Schmitt D P.Effects of radial layering when logging in saturated porous formations.JournaloftheAcoustical SocietyofAmerica，1988，84（6）：2200－2214.

［18］ 張海瀾，王秀明，張碧星.井孔的聲場和波.北京：科學(xué)出版社，2004.

Zhang H L，Wang X M，Zhang B X.Acoustic wavefield and waves in the borehole（in Chinese）.Beijing：Science Press，2004.

［19］ Tang X，Wang T，Patterson D.Multipole acoustic loggingwhile－drilling.／／72nd Annual International Meeting，SEG Expanded Abstracts，2002，364－368.

［20］ Press W H，Teukolsky S A，Vetterling W T，et al.Numerical recipes in Fortran 77：The art of scientific computing.Cambridge：Cambridge University Press，1992.

［21］ 楊文采.地震波場反演的BG一逆散射方法.地球物理學(xué)報，1995，38（3）：358－366.

Yang W C.BG－inverse scattering method for inversion of seismic wavefield.ChineseJ.Geophys.（in Chinese），1995，38（3）：358－366.