井孔与弹性波的相互作用：理论、方法和应用

唐晓明， 古希浩， 李杨虎， 苏远大*

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相同，在此不复赘述.