井中偶极声源在横向各向同性地层中的辐射波场

古希浩， 唐晓明,3*， 苏远大,3

1 中国石油大学(华东)地球科学与技术学院， 青岛 266580 2 中国石油大学(华东)深层油气重点实验室， 青岛 266580 3 青岛海洋科学与技术国家实验室海洋矿产资源评价与探测技术功能实验室， 青岛 266071

0 引言

各向异性地层中声源的辐射与波传播问题一直是地球物理领域内十分重要的方向(Levin，1979；Thomsen，1986；Wang，2011；周欣等，2020).特别是在井孔地球物理模型中，由于各向异性地层与充液井孔的共同作用，井中声源激发的弹性波场变得更为复杂(Fang et al.，2014).随着非常规页岩油气勘探开发的进一步深入，研究井中声源在各向异性(尤其是横向各向同性)地层中的辐射特性(刘金霞等，2012；孙伟家等，2013；李希强等，2013)，对井间地震走时的准确拾取，VSP剖面的正确追踪，特别是声波远探测技术(唐晓明和魏周拓，2012；Wang et al.，2015；Li and Yue，2017；Li et al.，2021；刘汇鑫等，2021)中井外地质体的精确刻画具有重要意义.明确各向异性和主频对声源辐射特性的影响对井下仪器的研制以及成像方法的改进具有指导作用.

多年来，国内外学者对井中声源的辐射特性已经开展了一系列研究工作.Meredith(1990)采用解析和数值相结合的方法研究了轴对称声源在井外近场与远场中的辐射指向性，为井间地震现场数据的理解提供了基本的模拟工具.随着偶极声波远探测技术逐渐成为研究热点，Tang等(2014)基于最速下降法，给出了井中偶极声源的远场辐射波形，通过与严格解所得结果对比，验证了该方法在辐射波场计算中的有效性.许家旗和胡恒山(2020)通过远场近似方法获得了充液井中偏心膨胀源的辐射波场，并阐明了声源偏心距与波场方位性之间的关系.但是上述研究皆设置井外地层为各向同性均匀介质，对井外各向异性地层的研究涉略较少.目前仅有少量学者进行了研究，如White(1982)考虑到非常规储层广泛存在的各向异性特征，推导出井中低频震源激发的弹性波远场位移的严格解，并讨论了波前面与各向异性的关系.Dong和Toksöz(1995)考察了井中无指向性的爆炸源在中等各向异性地层中的辐射指向性，但因受到当时计算机的计算精度限制，导致耦合分量的辐射指向性出现毛刺状尖峰，且其未考虑特殊各向异性地层中的三叉区.近年来，具有方位识别能力的多极声源在井下探测中广泛使用(王瑞甲和乔文孝，2015)，但是对于井中多极尤其是偶极声源在各向异性地层中激发的波场尚未见详细研究.而且与各向同性介质不同的是，各向异性介质中的弹性波速度与传播方向有关，因此无法根据最速下降法求取远场位移表达式的显式解(Dong and Toksöz，1995)，这也为各向异性地层井外波场模拟的实现增加了的难度.

本文首先对井中偶极声源在各向异性(主要考虑横向各向同性)地层中激发弹性波的远场位移表达式进行严格解析求解，然后运用最速下降法，借助群速度面与相速度面之间的转换关系，得到各向异性地层中准P(简称qP)、准SV(简称qSV)和SH波的远场位移渐近解，对比二者傅里叶变换后所得远场波形，验证了远场渐近解的精确性.最后根据渐近解表达式计算得到辐射指向性与辐射波形，着重研究了声源频率和各向异性对井中偶极声源辐射特性的影响.本文的结果可以为井外全波场的观测解释以及偶极声波远探测技术在各向异性地层中的应用提供理论支持.

1 VTI地层井外远场位移严格解

充液井中声源激发的弹性波场主要分为井内与井外两部分，各向异性地层情况下的井内波场基于声波导理论已经得到详实地模拟与分析(Ellefsen，1990；He and Hu，2009)，本文主要关注井外的辐射波场.采用图1所示的坐标系进行描述，设半径为a的无限长圆柱井孔中充满密度为ρf的流体，井外是以井轴为对称轴的横向各向同性地层，其各向异性特征由五个独立的弹性常数表征(唐晓明和郑传汉，2004).井中的偶极声源位于坐标原点，指向x轴正方向.θ为辐射波射线与z轴正方向的夹角，φ为辐射波射线和井轴所在矢面与声源指向之间的夹角，R为声源到辐射场点之间的距离，uR、uθ和uφ分别表示各向异性地层中qP、qSV和SH波产生的位移，其偏振方向与传播方向并不平行或垂直.

图1 描述井中偶极声源在横向各向同性地层中的辐射声场所采用的坐标示意图

根据霍姆赫兹定理，各向异性地层中弹性波的位移矢量场表示为：

(1)

(2)

其中，r0为偶极声源的分布半径，r和z分别为场点的径向和轴向距离，ω为圆频率，S(ω)为声源函数，k是轴向波数，f、qP、qSV和qSH分别表示流体、qP、qSV和SH波的径向波数(具体表达式参见(唐晓明和郑传汉，2004)).Kn是第二类变形贝塞尔函数，考虑偶极声源取n=1.由于各向异性造成qP波与qSV波的耦合，式(2)中qP和qSV波的位移势函数φ和Γ都包含qp、qsv两种本征波动，qP波中的qsv分量由b′控制，qSV波中的qp分量由a′控制，其中：

(3)

此外SH波可独立传播，其势函数只包含单纯的环向分量.

式(2)中井外弹性波分量的振幅系数可由矩阵方程(4)求取：

(4)

其中，An为井内偶极声波的振幅系数，表征着井内导波声场的强弱；Bn、Dn和Fn分别为地层中qp、sh和qsv波动的振幅系数，是井外辐射波场的重点考虑对象.方程左侧为Q矩阵与振幅系数的乘积，表示井内导波与地层中弹性波在井壁上产生的应力和法向位移，右侧b向量表示井中偶极声源在井壁处的直接贡献，井壁两侧的应力和法向位移连续，由此边界条件建立矩阵方程.其中Q矩阵和b向量的元素表达式由附录A给出.

将式(2)中的位移势函数代入式(1)，可得到各向异性地层中qP、qSV和SH波的位移分量的积分解：

(5)

井中偶极声源辐射到在VTI地层中的弹性波以球面波的形式向外传播，辐射波的位移表达式也可用球坐标系来表示，对于地层中的qP波：

(6)

式中，uR-p代表qP波中的qp分量，由柱坐标系下qP波的径向和轴向位移中的qp分量投影到球坐标系而得；uR-sv代表qP波中的qsv分量，由柱坐标系下qP波的径向和轴向位移中的qsv分量经坐标系转换得来.

对于地层中的qSV波：

(7)

式中，uθ-sv代表垂向夹角为θ方向上的qSV波中的qsv分量，由柱坐标系下qSV波的径向和轴向位移中的qsv分量投影到球坐标系而得；uθ-p代表qSV波中的qp分量，由柱坐标系下qSV波的径向和轴向位移中的qp分量经坐标系转换得来.此外，各向异性地层中SH波位移uφ在球坐标系与柱坐标系中并无二致，表达式相同.

2 VTI地层井外远场位移渐近解

井间地震、远探测测井和VSP剖面所涉及的波场辐射范围一般都远大于波长，满足远场条件|ξr|≫1(ξ为井外地层中弹性波的径向波数，r为声源与场点的径向距离).将最速下降法(Aki and Richards，1980)用于上述严格解，忽略O(1/R2)的高阶项，可得各向异性地层中qP波位移的远场渐近解表达式：

(8)

qSV波位移的远场渐近解表达式：

(9)

SH波位移的远场渐近解表达式：

(10)

最速下降法的求解步骤为：定义ψ(k)=ikcosθ-qsinθ，确定满足Ψ′(k0)=0的鞍点k0，使用不同振型(qp、qsv、sh)的轴向和径向波数的表达式可分别求取不同振型的最速下降解，其中kP0、kSV0和kSH0分别为qp、qsv分量和SH波轴向波数的最速下降解，qP0、qSV0和qSH0分别为qp、qsv分量和SH波径向波数的最速下降解，然后代入式(8)—(10)中相应的位移分量中可得各向异性地层中弹性波的位移远场渐近解，ψ″(k)的具体表达式由附录B给出.

在各向同性介质中，轴向波数的最速下降解为k0=ωcosθ/v，其中v为地层的弹性波速.在各向异性介质中，由于相速度与群速度具有差异，显式类型的最速下降解不复存在，此时可以借助相速度面与群速度面之间的转换关系来获得隐式的最速下降解(Dong and Toksöz，1995).如图2所示，井外为各向异性地层(参数见表1)，根据Kelvin-Christoffel方程(Carcione，2014)可求得弹性波相速度在特定方向的解析表达式：

(11)

式中可以看出各向异性介质中弹性波的相速度vph具有角散性质，速度大小随着相角θph改变.相速度vph与群速度vg可通过几何关系互相转换：

(12)

群角θg与相角θph也符合下述转换公式：

(13)

针对式(8)—(10)中给定的观测角度θ(观测方向与井孔的夹角，即群角θg)，得出qp、qsv和sh波动在群速度面上的对应点(见图2a)，然后根据式(12)和(13)可以反推找到给定的观测角度θ在相速度面上所对应的点(见图2b)，以确定vph和θph.各向异性介质中轴向波数的最速下降解为k0=ωcosθph/vph(对应式(8)—(10)中的kP0、kSV0和kSH0)，径向波数的最速下降解为q0=-iωsinθph/vph(对应式(8)—(10)中的qP0、qSV0和qSH0).值得注意的是在特殊各向异性地层(如石膏与土质夹层)中，各向异性参数满足关系式(Payton，1983)

图2 借助qp、qsv和sh波动的(a)群速度面和(b)相速度面来确定中等各向异性介质中辐射波波数的隐式最速下降解

(c13+c44)2-c11(c33-c44)>0，

(14)

或：

此时qsv波动的群角在某些范围内会出现三个群速度与之对应的情况(Auld，1973)，如图3所示.此时针对特定观测角度θg对应的三个群速度分别找到其相应的相速度，获取三个最速下降解，分别代入qsv分量的位移表达式可得其位移在特殊各向异性地层中的远场渐近解.

图3 借助qp、qsv和sh波动的(a)群速度面和(b)相速度面来确定特殊各向异性介质中辐射波波数的隐式最速下降解.这种情况下qsv波动的群速度面会出现三叉区

为考察井孔与各向异性对偶极声源的影响，基于式(8)、(9)和(10)分别定义qP、qSV和SH波的辐射指向性为：

(16)

将最速下降解代入式(8)、(9)和(10)可以计算任意频率下各向异性介质中qP、qSV和SH波的远场辐射波形，利用式(16)便可以得到三种波的辐射指向性.

3 数值计算结果及分析

本文从各向异性地层中偶极声场的远场渐近解出发，着重研究声源频率和各向异性对井中偶极声源的远场辐射特性的影响.

3.1 算法对比验证

为了证明该方法的可靠性，首先考察在砂岩条件下各向异性远场近似方法与Tang等(2014)论文中各向同性方法的自洽性，然后将奥斯汀白垩层中远场渐近解(式(8))模拟的波形结果与式(6)中的严格解计算的波形结果进行对比，模型中流体与地层的弹性参数见表1.

表1 模型参数

图4a是在砂岩地层情况下，使用各向异性远场近似方法和Tang等(2014)中各向同性方法计算的P波远场波形的对比结果，计算采用中心频率为3 kHz的Kelly声源，井孔半径为0.1 m；偶极声源指向x轴正方向，在xoz面内存在由10个间距为1 m的接收器组成的接收阵列，并且在距离声源5m处垂直于x轴放置.其中黑色实线为各向同性方法计算的结果，红色虚线为本文各向异性方法计算的结果.在砂岩地层中各向异性为零，群速度面与相速度面重合，式(3)中的耦合因子取值为零，P波位移中不包含qsv分量.两种方法计算所得波形的到时和幅度都十分吻合，表明各向异性远场渐近结果退化为各向同性时，与已发表的各向同性结果自洽.

图4b是井外为奥斯汀白垩层的情况下，qP波远场位移的严格解和远场渐近解计算的对比结果，黑色实线为严格解计算的辐射波形，红色虚线为渐近解计算的辐射波形.与砂岩地层不同的是，此时qP波辐射波场中存在qp和qsv两种波动形式，qp分量到时较早并且幅度较大，由于各向异性的存在而耦合出的qsv分量到时稍晚并且幅度较小.无论是qp分量还是qsv分量，两种方法计算所得波形都几乎完全重合，验证了各向异性条件下渐近解计算所得结果的正确性和可靠性.

图4 (a)砂岩地层情况下，P波位移的各向同性与各向异性远场渐近结果的自洽验证；(b)奥斯汀白垩层中，qP波位移的严格解与远场渐近结果对比

接下来利用各向异性条件下的远场渐近解，详细讨论声源频率和地层各向异性对充液井中偶极声源的辐射指向性的影响.

3.2 不同声源频率下的偶极辐射特性

在井孔地球物理勘探中，不同类型的勘探方法所使用声源的激发主频通常也在不同的区间，VSP和井间地震方法通常使用低频震源(100 Hz左右)，声波远探测中偶极声源激发频带涵盖2～5 kHz，相对其他地震成像方法具有分辨率高的优势(古希浩等，2020)，同时也存在更远距离成像的需求(董经利等，2020).在各向异性条件下，考察频率对井中偶极声源辐射特性的影响，优选出适宜的激发主频，对井下探测仪器的设计与改进具有重要的意义.

偶极声波远探测技术通常使用SH波对井外裂缝进行探测，因此我们首先设置式(16)中φ=90°，此时场点在yoz平面内，通过计算RSH来关注不同主频情况下奥斯汀白垩层中偶极SH波的辐射指向性(如图5所示)，其中环向刻度代表声源辐射方向与井轴的夹角，径向刻度为单位强度的声源辐射到地层中弹性波的相对幅值.从图中可以看出偶极声源向井外辐射的SH波较为纯净，不存在其他耦合分量.当声源主频为低频0.1 kHz时，由于波长远大于井径，声源的辐射并不受井孔的调制作用(Tang et al.，2014)，但与各向同性介质不同的是，各向异性地层中SH波的辐射指向性并非均匀辐射的圆形，而是轴向稍扁、径向略鼓的椭圆形(见图5a)，这说明水平偏振的SH波更倾向于沿着层理面传播，各向异性的存在使得SH波的辐射能量向径向压缩，这样对井外竖直裂缝等反射体的探测敏感度更高.随着频率的增加(图5a—f)，SH波的辐射指向性的幅度先增大后减小，在3 kHz和4 kHz之间激发强度最大，相比于砂岩地层中辐射能流峰值位于4.2 kHz(曹景记等，2016)，奥斯汀白垩层中优势激发频率向低频移动.而且在优势频率(3.5 kHz)附近SH波辐射指向性的覆盖性良好，径向探测能力突出，更有利于识别与井平行的反射体.

图5 不同声源主频条件下，井外各向异性地层中偶极SH波的辐射指向性

接下来设置式(16)中φ=0°，此时场点在xoz平面内，通过计算RqSV可得出不同主频情况下偶极qSV波的辐射指向性.与各向同性地层中SV波类似，qSV波中qsv分量(简称qSV-sv波，用空心圆圈表示)的辐射指向性呈现蝶形.受地层的各向异性影响，qSV波还存在式(3)中耦合项a′所贡献的qp分量(简称qSV-p波，用实线表示，并在右侧放大显示)，这是各向异性地层充液井中偶极声源辐射的重要特征.如图6a所示，当声源主频为低频0.1 kHz时，qSV-sv波的辐射指向性呈现上下对称的两瓣，相对于各向同性中的圆形(Tang et al.，2014)，各向异性的存在使得辐射能量向45°汇聚，这是因为垂直层理结构阻碍了能量的传播，而qSV-sv波在径向存在零点，导致能量只能从0°至90°的中间角度向外辐射.值得注意的是，qSV-p波的辐射指向性也是蝶形，在径向和垂向为零，说明层理面或者垂直层理面方向qSV波还原为严格意义上的SV波.随着激发频率的增加(图6a—f)，qSV-sv波的幅度变化不大，但辐射能量更加集中于45°方向，辐射覆盖范围也更加狭窄；qSV-p波的波幅随着频率的增加先增大后减小，在5 kHz左右激发强度最大，但是频率对其辐射形状影响不大.由于qSV波无径向探测能力，声波远探测中常不予考虑，但是井外若存在45°倾斜的反射体，该振型或许有微弱响应.

偶极qP波与qSV波在相同平面内，通过计算式(16)中的RqP得到其辐射指向性.从图中可以看出偶极声源向井外辐射的qP波中存在qp和qsv两种成分，其中qp分量(简称qP-p波，用空心圆圈表示)的辐射指向性与各向同性地层中的纯P波相似，以井孔为轴呈现出对称的两瓣，具有良好的辐射覆盖性.其耦合分量(简称qP-sv波，用实线表示)的辐射指向性的形状与上文中的qSV-p波相似，呈现出蝶形.从图7可以看出，qP波中耦合分量与主分量的幅度在同一量级，在低频(图7a)尤为明显，相比之下图6中qSV波中耦合分量要比主分量微弱的多，说明横波耦合分量对各向异性的敏感度要强于纵波.随着激发频率的增加(图7a—f)，qP-p波的幅度先增大后减小，在4～5kHz之间最易激发.值得注意的是，qP波的辐射指向性的幅度远小于比SH波，而qSV波无径向探测能力，因此SH波在辐射波场中依然占据主导地位，在各向异性地层依然需要使用SH波进行远探测成像.

图6 不同声源主频条件下，各向异性地层中偶极qSV波的辐射指向性

图7 不同声源主频条件下，各向异性地层中偶极qP波的辐射指向性

3.3 各向异性对偶极辐射特性的影响

接下来利用式(8)、(9)和(10)计算井中偶极声源激发，通过井外竖直平面内各个方向接收器所记录的波场，来观测辐射波场随地层各向异性的变化规律.声源主频设置为偶极测井常用的3 kHz，井外接收器以5m等源距环绕在声源周围，分别考虑地层为各向同性(选用奥斯汀白垩层中的c11和c66为各向同性地层的压缩和剪切模量)、中等各向异性(奥斯汀白垩层)和特殊各向异性(石膏-土质夹层)的情况，SH、qSV和qP波的辐射波场分别以各向同性地层中该类型辐射波的最大幅度为标准进行归一化，地层具体弹性参数见表1(数据来源于(White，1982)).

图8给出不同各向异性情况下，井外偶极SH波的辐射波场.在各向同性地层中(图8a)，SH辐射波良好地覆盖井外各个方向，受到井孔的调制作用，径向(90°)辐射波形幅度最大；此时弹性波传播速度与角度无关，到时面也呈现出规则的圆形.随着地层各向异性增强(图8b、c)，SH波的辐射能量向径向集中，其他方向的能量减弱，在特殊各向异性的石膏-土质夹层中，能量甚至极端汇聚到径向，大部分能量沿着水平层理方向辐射传播，这与上文辐射指向性的分析一致.受到弹性波速度角散性质的影响，不同接收器所记录SH波的传播速度不同，到时面呈现椭圆形甚至蝶形.此时若使用井中测量波速进行速度建模会造成反射体位置误差，因此在远探测成像时考虑各向异性校正(Tang et al.，2008)非常有必要.

图8 (a)各向同性、(b)中等各向异性和(c)特殊各向异性地层中，井外环向接收到的偶极SH波的辐射波场，声源主频为3 kHz

图9给出不同各向异性情况下，井外偶极qSV波的辐射波场.在各向同性地层中(图9a)，纯SV辐射波以声源为圆心在平面内均匀分布，但是波幅在径向和轴向为零.在中等各向同性地层中(图9b)，qSV波由主分量qSV-sv波和微弱的耦合分量qSV-p波两部分构成，qSV-sv波幅度在30°～60°较强，qSV-p波在50°～80°较为明显；两种波分别以qsv波动与qp波动的速度传播，由于受到速度角散的影响，它们的到时面也呈现出椭圆形.在特殊各向同性的石膏-土质夹层中(图9c)，qSV-sv波场在15°～75°出现三叉区(Triplication or cusps).从图中可以看到，在三叉区的尖点附近还会呈现出振幅奇异性，这是因为常规球面波的几何衰减为R-1，但是尖点处qSV-sv波传播的几何衰减为R-0.8(White，1982)，所以造成尖点方向的波形振幅的异常.

图9 (a)各向同性、(b)中等各向异性和(c)特殊各向异性地层中，井外环向接收到的偶极qSV波的辐射波场，声源主频为3 kHz

图10给出在不同各向异性情况下，井外偶极qP波的辐射波场.各向同性地层中的偶极P波向井外均匀辐射，到时面为圆形，波幅在径向略高.随着各向异性增加，在奥斯汀白垩层中井外波场除了qP-p波主分量外，在45°附近还耦合出具有相当可观的qP-sv波，受各向异性的影响，两种波的到时面为椭圆形，qP-p波在径向波速较快.在特殊各向异性的石膏-土质夹层中，qP-p波的到时面更是向径向集中压缩，其幅度在径向也是极为强烈.此外qP-sv波在特殊各向异性地层中也出现三叉区并产生振幅奇异性，与White(1982)文章中所得结论吻合，证实了基于远场近似理论的井外偶极波场模拟方法在各向异性地层中的有效适用性.

图10 (a)各向同性、(b)中等各向异性和(c)特殊各向异性地层中，井外环向接收到的偶极qP波的辐射波场，声源主频为3 kHz

4 结论

本文基于最速下降法，借助相速度与群速度之间的转换关系，给出了井中偶极声源在各向异性地层中辐射波场的快速模拟方法，并深入研究了声源激发主频与各向异性对偶极辐射特性的影响，得到以下几点认识：

(1)SH波在辐射波场中依然占据主导地位，各向异性的存在使得其辐射能量向径向压缩，其优势辐射频率也发生改变，由于到时面呈现椭圆形甚至蝶形，因此有必要在远探测成像时考虑各向异性校正.

(2)受地层的各向异性影响，qSV与qP波的辐射波场中均存在qp和qsv两种分量，qSV波的辐射能量集中于45°方向，井外若存在45°倾斜的反射体，该振型或许有微弱响应.横波耦合分量对各向异性的敏感度要强于纵波.

(3)在各向异性参数满足一定条件的特殊各向异性的地层中，qSV-sv和qP-sv波的辐射波场都会出现三叉区(Triplication or cusps)和振幅奇异性，这是特殊各向异性介质的重要特征.

附录A

Q矩阵和b向量中各元素的具体表达式为：

Q11=-nIn(fa)/a-fIn+1(fa)，

Q12=(1+ika′)[nKn(qPa)/a-qPKn+1(qPa)]，

Q13=nKn(qSHa)/a,

Q14=(b′+ik)[nKn(qSVa)/a-qSVKn+1(qSVa)]，

Q21=ρfω2In(fa)，

+2c66(1+ika′)[(n2-n)Kn(qPa)/a2

+qPKn+1(qPa)/a]，

Q23=2nc66[(n-1)Kn(qSHa)/a2-qSHKn+1(qSHa)/a],

+2c66(b′+ik)[(n2-n)Kn(qSVa)/a2

+qSVKn+1(qSVa)/a]，

Q31=0,

Q32=2nc66(1+ika′)[(1-n)Kn(qPa)/a2

+qPKn+1(qPa)/a],

+2qSHKn+1(qSHa)/a},

Q34=2nc66(b′+ik)[(1-n)Kn(qSVa)/a2

+qSVKn+1(qSVa)/a],

Q41=0,

-qPKn+1(qPa)],

Q43=iknc44Kn(qSHa)/a,

-qSVKn+1(qSVa)],

b1=εn[nKn(fa)/a-fKn+1(fa)]，

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

b3=0，b4=0，b4=0,

其中a是井孔半径，εn是纽曼因子，当n=0时εn=1，当n>0时εn=2.

附录B

计算qp分量时ψ″(k0)的表达式为：

(B1)

其中：

(B2)

ψ″SV0(kSV0)的表达式与式(B2)类似，只需将式中径向波数与轴向波数进行替换.

计算SH波时ψ″(k0)的表达式较为简便：

(B3)

式中参数的含义在正文中已给出.