井孔與彈性波的相互作用：理論、方法和應用

唐曉明， 古希浩， 李楊虎， 蘇遠大*

1 中國石油大學(華東)地球科學與技術學院， 青島 266580 2 青島海洋科學與技術國家實驗室海洋礦產資源評價與探測技術功能實驗室， 青島 266071 3 中國石油大學(華東)深層油氣重點實驗室， 青島 266580

0 引言

地殼巖石中鉆井與彈性波的相互作用是一個經久不衰的科學問題，從早期各種井中的地震勘探測量(如垂直地震反射剖面、井間地震等)，聲波測井，到近期發展起來的聲波遠探測技術，都與該問題密切相連.地震測量主要關心井孔對入射的低頻彈性波場的調制和響應(參見White，1953；Schoenberg，1986；Peng et al.，1993)，而聲波遠探測則要考慮中高頻段的井中聲源對井外地層的輻射波在地層中的反射，以及在井中對反射聲波接收的全過程 (唐曉明和魏周拓，2012; Tang et al.，2014; Xu et al.，2019).對于偶極橫波遠探測技術，還必須考慮聲源發射和接收的方位效應.除了波的輻射和接收，另一個重要的問題是井對彈性波的反射和散射，這個問題由于遠探測技術在鄰井探測(蘇遠大等，2014；Tang et al.，2016)和井叢防碰中的應用變得重要起來.

綜上所述，我們將井中的輻射和接收以及井對波的散射統一歸結為井孔與彈性波的相互作用問題，并結合遠探測技術的應用對此進行理論分析、求解.遠探測波場的特點是波的井外輻射及對井的入射都具有球面波的傳播特征，而很多早期，甚至包括近期 (如Hirabayashi et al.，2017)對井孔的彈性波入射分析大都采用平面波入射的假設.許家旗和胡恒山(2019)結合虛源法和互易定理求解井孔的球面波入射問題，將輻射到地層界面、反射回井孔的入射波等效為來自虛源的球面波，在此基礎上利用虛源與井中聲源的互易關系得到該球面波在井中的激發響應.本文采用將入射球面波展開為柱面波的方法(Tang and Cheng，2004)，可以像求解井孔的輻射問題那樣，方便地得到對井孔入射問題的解.

遠探測波場的分析結果為近期發展起來的聲波遠探測技術奠定了理論基礎并獲得了良好的應用效果.鑒于該課題的重要性和挑戰性, 許多學者在各個方面做了大量工作.這些工作和結果在相關的文獻中都已有介紹；本文的重點在于介紹和討論在實際應用中被證明是行之有效的理論分析和結果，這些應用包括運用井孔散射理論的鄰井探測和具有方位特征的偶極聲波遠探測.

1 井中聲源輻射的理論分析

圖1給出了利用井中聲波測井儀對井旁地質反射體進行探測的示意圖，儀器上的聲源向井外輻射彈性波，經反射后入射到井中，被儀器上的探測器接收.采用柱坐標系，井中流體中位于坐標(r0,φ0,z0)處的點聲源向場點(r,φ,z)輻射的球面波可寫成柱面波展開的形式(Tang and Cheng，2004):

+sin(nφ)sin(nφ0)]eik(z-z0)dk，

(1)

將(1)式的聲源置于井中，將在井內流體和井外彈性地層中激發出彈性波動，井內外的位移場可分別表示為

(2)

其中φf和φ分別是流體和地層中縱(P)波的位移勢，和Γ分別是地層中SH和SV橫波的位移勢，為z軸方向的單位向量.頻率-波數域中各個位移勢函數的通解一般表達為

(3)

(4)

式中的上標d和f分別表示聲源輻射的直達聲場和井壁在流體中產生的反射聲場.由(1)、(2)和(3)式計算出直達波和井壁在流體和地層中的位移和應力，代入(4)式得到兩個矩陣方程：

M×[AnBnCnDn]T=b,

M×[A′nB′nC′nD′n]T=b′，

(5)

式中的M為一4×4矩陣；b和b′為4×1向量，它們的表達式由附錄A給出.(3)式中的波幅系數確定后，聲源在井內流體和井外地層中激發的波場就完全確定.舉例來說，由(3)式得到地層輻射的SH位移勢的頻譜函數為

×Kn(sr)eik(z-z0)dk.

(6)

與遠探測有關的波場輻射范圍，一般都在遠大于波長的距離之外，滿足|sr|?1的條件，由此可對(6)式中的波數積分采用最速下降方法求解(郭敦仁，1978；Tang and Patterson, 2009)，得到

(7)

圖1 井外地層反射體聲波遠探測示意圖Fig.1 A schematic model of a borehole with a reflector in the formation

將(1)式中的點聲源加以組合，可以構造出各種形態的聲源，最為典型的是測井常用的多極子聲源.例如，將(1)式中的點源居中放置(取r0=0), 這時(3)式的級數求和中只有n=0項的貢獻，即輻射是軸對稱的，不產生SH波，這時該式地層的位移勢函數中=0，而φ和Γ分別對應于居中單極聲源輻射的P波和SV波.這個問題在Meredith(1990)和唐曉明和魏周拓(2012)論文中討論過.

對于偶極聲源，可以把(1)式所示的兩個點源放在(r0,φ0,z0)和(r0,φ0+π,z0)的位置上；二者的激發聲源函數S(ω)大小相等，符號相反，在兩源足夠靠近(r0→0)的條件下形成偶極輻射，這時(3)式的級數求和中只有n=1項的貢獻，在類似于(7)式的遠場漸近條件下，可以得到φ、和Γ的表達式，分別求出P、SH和SV波的遠場解及其輻射指向性.偶極輻射的問題也有多人討論過(Tang et al.，2014；Xu et al.，2019).類似地，把四個正負相間的點源相對于井心對稱放置，可組成四極聲源，不復贅述.

對于分布式的聲源結構，可將(1)式的點聲源沿結構的表面積分，點源的激發強度由其在結構表面的分布函數決定，由此可得任意聲源在井外的輻射波場.一個很好的例子是隨鉆聲波儀器上的環狀聲源(唐曉明和鄭傳漢，2004；Xu et al.，2018)，將(1)式按強度分布為cos(nφ0)的圓環積分，可得到環狀多極子聲源及其激發的波場.

2 井對入射彈性波的響應

以上的討論給出了井中聲源在井外地層產生的輻射波場的理論分析和求解方法.采用相同的理論和方法，也可以分析井對彈性波入射的響應，這種響應包括了兩個重要方面：一是在井內流體中對入射波場進行接收，二是井對入射波場的反射或散射.

井中聲源輻射的球面波(參考(7)式)被地層中的反射體反射后，在反射體尺寸大于波長的條件下，反射波也以球面波的形式傳播(Aki and Richards，1980)，其形式如下：

(8)

(9)

(1)式對(8)式的應用將界面反射波轉化為來自虛源的球面入射波；而上述貝塞爾函數的調整使得(9)式的波場通解可以用來求解球面波對井孔的入射問題，避免了對此問題常用的平面波入射假設(例如White，1953；Schoenberg，1986；Peng et al.，1993；Hirabayashi et al.，2017).確定(9)式中波幅系數的邊界條件方程(4)也要做相應的改變，即把井中聲源在井壁處產生的反射和向外輻射變為井外虛源在井壁處產生的散射和透射.

(10)

式中的上標i，s和f分別表示(虛源)入射、散射和井中流體的透射波場.(10)式中的入射波即為井中聲源輻射到地層，遇到反射體后的反射波，其形式與(8)式類似.以(7)式中的SH波勢函數為例，將式中的輻射距離R用(8)式中的傳播距離D置換，再運用(1)式, 就得到其入射波的表達式.根據這一推導，可以得到入射P、SH和SV波勢函數的一般形式

(11)

式中的下標P、SH、SV分別表示P、SH、SV波，RD為井中聲源的輻射函數(其幅值即為隨角度(θ,φ)的變化的輻射指向因子)，由求解(5)式中的波幅系數得到；ρ為地層密度； RF為反射體的反射系數(參見(8)式).

與從(2)、(3)和(4)式推導出(5)式的過程類似，可從(2)、(9)、(10)和(11)式推導出與(5)式形式相同的矩陣方程組.矩陣M與(5)式的相同(在附錄A中給出).但是，代表聲源激發貢獻的向量b和b′與井中流體激發時很不相同，因為它們代表著來自虛源的P、SH、SV波的貢獻，可根據(2)式和(11)式計算得到，具體形式由附錄B給出.

向量b和b′確定之后，便可從(5)式中求解出(9)式中的波幅系數，從而得到井中的透射和井外的散射波場.散射波位移勢的譜函數可寫為

(12)

(13)

其中ρf為流體密度.與(6)式中的輻射分析一樣，(12)和(13)式中的波數積分可用最速下降法準確求解.兩式貝塞爾函數Kn(kvr′0)中的虛源的徑向距離遠大于波長，滿足|kvr′0|?1的條件.以(13)式中的流體徑向位移為例，其最速下降解為

(14)

3 理論的應用

3.1 基于井孔散射分析的鄰井探測和應用實例

井孔對彈性波的散射問題在地震測量中一直不受重視，原因是缺少應用需求.隨著井中遠探測技術的發展和需要，其重要性日益突出，主要的應用是鄰井探測：對兩個相鄰的鉆井，把其中一個作為目標井，另一個作為探井放置測量儀器，用測得的聲波數據對目標井成像.鄰井探測的一個重要應用是驗證偶極橫波遠探測技術.在中國石油大學(華東)鉆了兩口相距10 m的井(地層巖性為火山凝灰巖)，SH橫波的探測結果獲得了目標井的良好成像.相關的理論分析和結果可參見Tang 等(2016).鄰井探測的一個新近的重要應用是井叢的防碰問題.在海上油田密集的井叢中鉆井，井之間的防碰問題對油田的安全生產至關重要.密集的井叢一般位于淺海松軟的地層中，對此采用偶極聲源產生的低頻縱波，而不是橫波，會獲得較好的應用效果.

鄰井探測的理論分析如圖2所示，包括了波與目標井和探測井的兩次入射/散射過程.波從探井出發，經目標井散射后回到探井，波在探井井壁處被再次散射，其透射部分被井中儀器接收后進行成像處理.

圖2 鄰井探測的理論分析模型Fig.2 Analysis model for delineating a target borehole using elastic waves from a nearby well

參見前面對偶極聲源輻射的描述及遠場漸近解，得到其輻射到目標井的縱波位移勢函數為

(15)

式中R為探井至目標井的傳播距離，θ為波的出射角，RDP是偶極聲源的縱波輻射因子.

-pKn+1(pr)]Kn(kvr′0)eik(z-z0)dk，

(16)

式中的波幅系數Bn和B′n需在(15)式中的縱波對目標井入射的條件下求解(5)式得到.在井間距離r′0遠大于波長及(入射方向)波長在井軸方向投影的條件下(|kvr′0|?1，|pr|?1)得到(16)式在探井周圍的最速下降解為：

+B′nsin(nφ)]eiω(R′+r′0sinθ)/α-πi/4，

(17)

在探井中接收(17)式的散射波，從原理上可以再次求解(5)式的矩陣方程得到井中流體的透射波，也可以用大家熟知的源與接收之間的互易性來簡化這個問題.對于(17)式的入射縱波位移，其在井中接收到的流體位移可表示為

(18)

式中RC稱為井的接收響應函數(Peng et al.，1993).在源和接收器都位于井軸的條件下，運用彈性波的互易原理，可以證明(Tang et al.，2014；Xu et al.，2019)

RC(ω)=RD(ω)，

(19)

即井的聲源輻射函數與井對同一聲源產生的波場的接收響應函數是等同的.這樣一來，遠探測模擬分析中波從井中的出射和入射只需使用同一個RD函數，大大提高了計算效率.

我們用一個計算實例來驗證以上分析的正確性.表1給出了流體和地層(軟地層)的聲學參數，井徑為0.2 m，探測井與目標井的距離為5 m.井中放置偶極聲源和8個接收器，其至源的井軸距離為3到4.05 m，聲源的激發函數是中心頻率為3 kHz的Kelly子波.為了模擬探井中聲源激發后來自目標井的縱波信號，分別采用了解析方法((15)到(18)式)和有限差分數值方法，差分法采用的空間網格和時間步長分別是4.5 mm和0.42 μs.圖3給出了解析法(實線)和差分法(虛線)的比較，二者的結果幾乎完全一致，但計算效率差別巨大，差分法需耗時十幾小時，而解析法僅需數十秒便得到圖中的結果.

表1 模型參數Table 1 Model parameters

圖3 分別利用有限差分方法和漸近解析方法計算的鄰井散射縱波的結果對比Fig.3 Comparison of compressional waves scattered from a nearby well. Solid and dashed curves are calculated using analytical and finite difference modeling, respectively

為了說明鄰井探測的重要應用，圖4給出了淺海井叢防碰的一個應用實例.第一道給出了30～60 m井段地層的自然伽馬(GR)和縱波慢度(DTP)的測井曲線.對于這種淺海超軟慢速地層，縱波較橫波易于激發且所需記錄時間較短，適宜于采用上述的縱波鄰井探測的分析方法.第二道給出了偶極測井的全波波形的變密度圖，聲源的中心頻率為3 kHz.圖中高波幅的振相為沿井筒傳播的直達縱波(對軟地層稱之為泄漏式縱波).對全波信號進行分波處理后提取的鄰井散射波的變密度圖見第三道，處理后的直達波相對于原始數據被大幅壓制，較弱的鄰井散射波得以顯示.將其用第一道中的縱波慢度偏移成像后，給出了井旁的一個高角度反射體，如第四道中的近豎直條帶所示，距井僅5～6 m.經地表井位分布圖驗證后，證實這是井旁的一個高角度鄰井.對于間距僅數米的鄰井，鉆井繼續下鉆需要不斷檢測兩井之間的距離，防止碰撞發生.

圖4 偶極縱波鄰井探測在淺海井叢防碰中的應用實例Fig.4 Field example of nearby borehole target delineation in a shallow-marine densely drilled area using dipole P-wave data

3.2 四分量偶極橫波遠探測數據的模擬分析和成像實例

偶極橫波遠探測采用四分量數據(xx，xy，yx，yy)可以確定反射體距井的位置和走向(唐曉明和魏周拓，2012).根據(14)式給出的井中流體位移，可以對四分量數據進行精確的模擬和分析.實際四分量正交偶極儀器的數據接收的組合如圖5所示，x(或y)向換能器偏心放置.將位于+x和-x(或+y和-y)軸兩個換能器的聲壓相減，便得到x(或y)向的接收位移.例如，取(13)式中的聲壓徑向導數?P/?r，再除以ρfω2，便得到(14)式中的徑向位移.在r很小(r→0)的條件下，級數求和中只有A1cosφ+A′1sinφ這一項，對應于儀器居中時的偶極響應.可見，對稱偏心換能器的數據組合，是居中測量流體位移的一個很好的近似.在居中接收的條件下，可以采用前述的源和接收的互易性(如(19)式)來模擬四分量數據.以交叉分量為例，這樣模擬的數據滿足xy=yx，因為根據該互易性原理，y-向接收的由x-向偶極源發射的數據必須等于x-向接收的由y-向偶極源發射的數據.

圖5 井中x和y軸指向的接收器接收入射波場的示意圖Fig.5 Sensor configuration in x and y coordinates for receiving incident waves

但是，如果接收的位移不居中，例如在偏心距為r的+x和+y軸上接收yx和xy位移，這樣的數據就不滿足互易性，即xy≠yx.由此產生的一個重要問題是: 實際測量中儀器偏心對四分量偶極數據的采集和處理會有多大的影響? 我們用(14)式來模擬分析這種情況.采用圖1中的模型，反射體為表1硬地層a和b之間的界面，傾角為60°，走向為150°(從x軸方向順時針計算).井中儀器與圖3相同，x和y換能器的偏心距r=0.04 m.聲源到反射體一側的虛源徑向距離為30 m，源的中心頻率為3 kHz，由此模型計算得到的反射橫波數據如圖6a所示.其中下、中、上圖分別顯示了計算得到的SV、SH和SH+SV(即模擬數據)的四分量(xx，xy，yx，yy)波形圖.顯然xy與yx的幅度出現了差別(xy>yx)，說明x和y接收器不在同一位置上時，上述的互易性就不嚴格成立.從模擬結果還可以看到, SV分量比SH分量的幅度小很多，說明偏心接收的數據SH+SV中仍以SH的貢獻為主.

圖6 (a) 圖1模型計算的偏心接收的四分量偶極橫波遠探測數據(上圖)，其中的SH和SV波的貢獻由中、下圖分別給出； (b) 四分量數據旋轉得到的波幅方位分布圖(實線，居中接收結果由圖中虛線給出)，其峰值對應著反射體走向Fig.6 (a) Simulation of four-component dipole shear-wave data (xx, xy, yx, yy) acquisition for the eccentric reception (upper figure) for the model in Fig.1. The middle (lower) figure shows the SH (SV) component of the modeled data. (b) The azimuthal rotation of wave amplitude (solid line, dotted line is the centered reception result) for the 4C data, with the peak amplitude corresponding to reflector strike

將圖6a上圖中的四分量數據，按實際偶極數據處理那樣，做四分量旋轉，組成一個新的XX分量(Li et al.，2019; 李盛清等，2020)

XX=xxcos2φ-sinφcosφ(xy+yx)+yysin2φ,

(20)

其中φ為x軸方向順時針計算的旋轉角.圖6b的極坐標圖給出幅值|XX(φ)|隨旋轉角φ的變化曲線，曲線的極大值(極小值)正好對應反射體的走向(傾向), 這正是四分量偶極數據求取反射體方位的理論依據(這是因為換能器居中時xy=yx, 四分量數據組成的對稱矩陣存在最大和最小兩個特征方向, 如圖中的虛線圖形所示).圖中換能器居中(虛線)與偏心(實線)的結果對比還表明：測井時接收換能器的偏心并不太影響反射體方位的確定.但是，按(20)式得到的反射體(走向)的方位存在180°不確定性，因為圖6b圖形的峰值大小相等、方向相反(即差180°).這時，考慮偏心和居中接收的差別(如圖6b所示)有助于解決180°不確定性，值得進一步研究(許家旗和胡恒山，2020).

四分量偶極橫波遠探測技術在油氣勘探和生產中已經取得了良好的應用效果.我們用一個遠探測方位成像的可重復性實例來說明這種方法的方位特征和可靠性.在西北某油田六千多米的深井井段做了四分量偶極橫波的重復測井.圖7給出了40 m井段內兩次測量數據的成像處理結果.如第一道所示，用Run 1和Run 2表示的兩次測量過程中儀器都在井中旋轉，且儀器的方位(即x或y極板的指向)很不相同.將兩套四分量數據分別用(20)式做旋轉變換，再用圖6b所示的方式求取反射體的走向，得出的結果十分一致，都在北北西/南南東(NNW/SSE)方向, 如圖7右邊的方位玫瑰圖所示.沿此方位所做的兩套數據在井旁20 m范圍內的遠探測成像分別在第二道和第三道給出.兩圖中數米到數十米尺度的反射體成像都十分一致(個別除外).特別是圖中箭頭所指的高角度反射體，被解釋為NNW/SSE走向的裂縫/斷層.該裂縫在六千多米深處的地層壓力下還能清晰成像，說明其被高壓油氣充注呈開啟狀態，是油氣運移的重要通道.對遠探測技術來說，在儀器聲源發射和接收任意指向的條件下，對地層反射的小振幅信號(相對于井筒直達波而言)成像的可重復性說明了技術的可靠性和結果的準確性.

圖7 四分量偶極橫波遠探測方位成像重復性驗證實例. 盡管是在兩次測量(Run 1和Run 2)完全不同的儀器方位(第一道)條件下，但是反射體走向的識別結果(第四道)在成像圖(第二道和第三道)中具有相當好的可重復性Fig.7 Repeatability test of four-component dipole shear-wave reflection imaging. Despite the very different tool rotation (Panel 1) for two logging passes (Run 1 and Run 2), repeatable imaging results (Panels 2 and 3) are obtained along the reflector strike (Panel 4)

4 結語和展望

隨著井孔彈性波測量技術的發展和需求，井孔與彈性波相互作用這一經典問題變得日益重要起來，這種相互作用主要包括井中聲源向外的輻射以及波從井外反射體向井的入射、散射和透射.將井孔波場的柱面波展開和波數積分通解與波的遠場漸近解及井孔輻射和接收的互易性相結合，得到了求解這一問題的快速有效方法.根據這一方法，可以準確模擬該問題中各種有關的波動現象并分析其波場特征，這些分析結果為近期發展起來的聲波遠探測技術奠定了理論基礎并獲得了良好的應用效果.具體而言，對偶極聲源的輻射特性及其遠探測波場的方位特征進行四分量數據采集能夠對井外反射體準確成像并確定其走向; 超軟慢速地層鉆井中縱波的輻射和接收響應特征分析促進了鄰井探測技術在淺海井叢防碰工程中的應用.遠探測技術的進展和應用也帶來了更多的需求和亟待解決的問題.例如，該技術在深部油氣藏中探測到許多斷溶構造，油氣開發需要確定斷溶體的方位，為此必須克服偶極探測方位的180°不確定性.這是一個挑戰性的課題，不僅需要更深入的理論探討和分析，還需要在實際應用中驗證分析結果的有效性.

致謝謹此祝賀陳颙先生從事地球物理教學科研工作60周年.

附錄A

矩陣方程(5)式的完整形式如下：

(A1)

和

(A2)

矩陣M中各元素的具體表達式如下：

M11=-nIn(fa)/a-fIn+1(fa),

M12=nKn(pa)/a-pKn+1(pa),

M13=nKn(sa)/a,

M14=ik[nKn(sa)/a-sKn+1(sa)],

M21=ρfω2In(fa),

M22=2ρβ2[(n2-n)Kn(pa)/a2+pKn+1(pa)/a]

+ρ(2k2β2-ω2)Kn(pa),

M23=2nρβ2[(n-1)Kn(sa)/a2-sKn+1(sa)/a],

M24=2ikρβ2{[(n2-n)/a2+s2]Kn(sa)

+sKn+1(sa)/a},

M31=0,

M32=2nρβ2[(1-n)Kn(pa)/a2+pKn+1(pa)/a],

M33=-ρβ2{[2(n2-n)/a2+s2]Kn(sa)

+2sKn+1(sa)/a},

M34=2iknρβ2[(1-n)Kn(sa)/a2+sKn+1(sa)/a],

M41=0,

M42=2ikρβ2[nKn(pa)/a-pKn+1(pa)],

M43=iknρβ2Kn(sa)/a,

M44=-(s2+k2)ρβ2[nKn(sa)/a-sKn+1(sa)].

(A3)

聲源位于井內流體中時，向量b和b′中各元素的具體表達式如下：

b1=εn[nKn(fa)/a-fKn+1(fa)]cos(nφ0),

b2=-εnρfω2Kn(fa)cos(nφ0),

b3=0,b4=0,

b′1=εn[nKn(fa)/a-fKn+1(fa)]sin(nφ0),

b′2=-εnρfω2Kn(fa)sin(nφ0),

b′3=0,b′4=0.

(A4)

附錄B

對于位于井外的虛源，入射波類型不同，其向量b和b′中的元素也不一樣.

b1=-εn[nIn(pa)/a+pIn+1(pa)]cos(nφ′0),

b2=-εnρ{2β2[(n2-n)In(pa)/a2-pIn+1(pa)/a]

+(2k2β2-ω2)In(pa)}cos(nφ′0),

b3=-2nεnρβ2[(1-n)In(pa)/a2-pIn+1(pa)/a]

×cos(nφ′0),

b4=-2εnikρβ2[nIn(pa)/a+pIn+1(pa)]cos(nφ′0),

b′1=-εn[nIn(pa)/a+pIn+1(pa)]sin(nφ′0),

b′2=-εnρ{2β2[(n2-n)In(pa)/a2-pIn+1(pa)/a]

+(2k2β2-ω2)In(pa)}sin(nφ′0),

b′3=-2nεnρβ2[(1-n)In(pa)/a2-pIn+1(pa)/a]

×sin(nφ′0),

b′4=-2εnikρβ2[nIn(pa)/a+pIn+1(pa)]sin(nφ′0).

(B1)

b2=-2nεnρβ2[(n-1)In(sa)/a2+sIn+1(sa)/a]

×sin(nφ′0),

b3=εnρβ2{[2(n2-n)/a2+s2]In(sa)-2sIn+1(sa)/a}

×sin(nφ′0),

b4=-iknεnρβ2In(sa)sin(nφ′0)/a,

b′1=nεnIn(sa)cos(nφ′0)/a,

b′2=2nεnρβ2[(n-1)In(sa)/a2+sIn+1(sa)/a]

×cos(nφ′0),

b′3=-εnρβ2{[2(n2-n)/a2+s2]In(sa)-2sIn+1(sa)/a}

×cos(nφ′0),

b′4=iknεnρβ2In(sa)cos(nφ′0)/a.

(B2)

b1=-ikεn[nIn(sa)/a+sIn+1(sa)]cos(nφ′0),

b2=-2ikεnρβ2{[(n2-n)/a2+s2]In(sa)

-sIn+1(sa)/a}×cos(nφ′0),

b3=-2iknεnρβ2[(1-n)In(sa)/a2

-sIn+1(sa)/a]cos(nφ′0),

b4=εn(s2+k2)ρβ2[nIn(sa)/a+sIn+1(sa)]

×cos(nφ′0),

b′1=-ikεn[nIn(sa)/a+sIn+1(sa)]sin(nφ′0),

b′2=-2ikεnρβ2{[(n2-n)/a2+s2]In(sa)

-sIn+1(sa)/a}sin(nφ′0),

b′3=-2iknεnρβ2[(1-n)In(sa)/a2-sIn+1(sa)/a]

×sin(nφ′0),

b′4=εn(s2+k2)ρβ2[nIn(sa)/a+sIn+1(sa)]

×sin(nφ′0).

(B3)

其中，S為聲源函數；RD為輻射函數；RF為反射系數；φ′0為虛源的方位；其他符號定義與附錄A相同，在此不復贅述.