1 引 言

对于砂泥岩薄互储层和裂缝型油气储层这两类非常规油气储层的准确评价远比常规油气储层复杂和困难，因为这类储层表现出电各向异性特征，即其水平电阻率和垂直电阻率不同，在水平分层垂直井中，常规的感应测井只对水平电阻率敏感，其测量结果主要受地层中低电阻率部分的影响，由此得出较低的含油饱和度，从而漏失这类储层^[1].目前，三分量感应测井和定向测井两种仪器在评价这类储层上较原有的感应测井仪器有明显的进步.前者主要采用三对互相垂直的发射接收线圈组来测量评价井周油气层，这类测井仪器在电缆测井中得到了较快的发展^[1-9]，后者主要采用不同法向角度的线圈组合(以下称为倾斜线圈系)来识别评价地层的电各向异性特征，主要应用于随钻测井环境中^[10-11].倾斜线圈系除了能够在导向井中识别地层边界，修正井眼轨迹外，也可以识别和评价这类各向异性储层，但有关定量求取电阻率各向异性参数方面的研究相对较少^[12]，本文通过采用数值模式匹配法高效模拟倾斜线圈系在横向各向同性地层的响应，为未来定量评价各向异性参数提供了坚实的基础.

为了分析此类倾斜线圈系在各向异性地层中响应特征，需要借助数值模拟的方法.2004 年，Hue等^[10]采用三维有限差分法(Finite Difference Time Domain method，FDTD)模拟倾斜线圈系在地层中的响应，由于三维有限差分法的运算速度较慢，2007年，采用数值模式匹配法模拟倾斜线圈系在地层中的响应，并在吸收边界处采用了完全匹配层(Perfectly Matched Layers，PML)^[11].国内，杨锦舟等^[13]将线圈系简化成磁偶极子，利用并矢Green函数计算倾斜线圈系在水平层状介质中的响应，重点研究了接收线圈法向角度对接收信号的影响.

本文详细给出了模拟倾斜线圈系在横向各向同性地层中响应的数值模式匹配法(Numerical Mode-Matching，NMM)，该方法是一种半有限元半解析解，它的效率比其它方法(如有限元素法，有限差分法)要高很多，而且非轴对称的倾斜线圈系是一种三维源，即便是在轴对称地层模型的条件下，如果采用有限差分法或是有限元素法必须考虑三维模型^[10].在数值模式匹配法中，采用垂向本征模式描述非轴对称源的积分形式，因而能在二维条件下求解仪器响应特征，由此进一步降低了运算量，提高效率.数值模式匹配方法在感应测井模拟应用起源于20 世纪80 年代^[14].Tsang等^[15]采用数值模式匹配法对普通电阻率测井的数值模拟表明，其效率是有限元素法的数倍.在前人工作的基础上，Pai^[16]发展了垂向数值模式匹配法.Liu^[17]采用数值模式匹配法模拟偏轴条件下的电磁波测井响应，从而使该方法不再局限于二维轴对称模型.通过对椭圆井眼和非轴对称的侵入条件下的侧向测井响应模拟，Guo^[18]将该方法扩展到三维模型.国内，张庚骥等^[19-21]最早将此方法应用于电法测井响应的数值模拟中，并且通过基函数的改进，极大地提高了该方法的计算效率和精度.沈金松^[22]利用垂向数值模式匹配法研究了多频电磁波测井响应，并计算了水淹层和低电阻率储层的高频电磁波测井响应.

本文在借鉴了Hue的工作基础上，详细推导了利用垂向数值模式匹配方法模拟倾斜线圈系在横向各向同性地层的表达式，并通过与解析解和有限差分解对比，验证了该方法的准确性.由于采用了幅度和斜度基函数，运算效率得到了有效提高.

2 基本原理

考虑地层模型由N个水平层状横向各向同性介质组成(如图 1)，各个地层的复介电常数张量、磁导率张量和电导率张量表示为ε=diag(ε_H，ε_H，ε_V)，μ = diag(μ_H，μ_H，μ_V)，σ = diag(σ_H，σ_H，σ_V)，其中ε_H = [ε_rε₀ - (σ_H/iω)]S_z，ε_V = [ε_rε₀-(σ_V/iω)]/S_z，μ_H =μ₀S_z，μ_V =μ₀/S_z，在z方向上采用完全匹配层^[23]，截断无限求解域(Sz表示衰减因子).σ_H 和σ_V 分别表示水平电导率和垂直电导率，ε_r 表示相对介电常数，ε₀ 为真空中的介电常数，μ₀ 为真空中的磁导率常数，i=.模拟时将模型在径向上划分多个同轴环形域，在每个环形域中利用一维有限元方法求解本征模式，然后在径向界面处通过电场和磁场的连续性条件确定各本征模式之间的反射系数矩阵和透射系数矩阵.

由于发射源是倾斜的，将激发出非轴对称的交流电:J (ρ，φ，z)=I_Tδ(ρ$\hat{\rho }$-ρ_T)${\hat{\delta }}$(z-ζ_T)(${\hat{\phi }}$+${\hat{z}}$ξ_T)，其中$\hat{\rho }$，${\hat{\phi }}$，${\hat{z}}$表示单位矢量，I_T 为电流，ζ_T(φ)=z_T+ρ_Ttan(θ_T)cos(φ-φ_T)，ξ_T(φ)=tan(θ_T)sin(φ-φ_T)，ρ_T 表示电流源在水平面上的投影半径，z_T 表示电流源的z向坐标位置，θ_T 表示法向角度，φ_T 表示方位角度(本文设为0°)，电磁场的时间变化因子为e^-iωt，ω 为发射源角频率，E表示电场强度，H表示磁场强度，倾斜电流源激发的电场和磁场可以表示成如下矢量波动方程^[11]:

(1)

(2)

以下给出利用垂向数值模式匹配方法求解上述方程的原理及方法.

2.1 电磁波的轴向分量

在轴对称条件下，由于电场E的z方向分量e_z可以表达为^[23]:

(3a)

(3b)

上式中，J_ν(k_ρρ)表示ν阶贝塞尔函数.同理，对z方向的磁场也可以进行谱分解.

将式(3)代入(1)式，可得z方向的ν 阶谱分量e_zν(ρ，z):

(4)

式中，d_zν 为电通密度，其与电场强度e_zν 的关系为:d_zν =ε_Ve_zν，j_z表示z方向的电流密度，j_z(ρ，z)j=

同理，可以得出z方向的ν 阶谱分量b_zν(ρ，z):

(5)

式中，b_zν 为磁通密度，其与磁场强度h_zν 的关系为:b_zν =μ_Vh_zν，j_φ 表示φ 方向的电流密度，j_φ(ρ，z)=

在无源区域，采用(3b)对e_zν(ρ，z)进行汉科尔变换，(4)式可以化成^[11]:

(6)

在z方向上采用基函数将d_zν 展开，

(7)

式中，Λ_α 为第α 个区间的基函数(α，β=1，2，…，M，M表示z方向展开电通密度的区间总数)，将(7)式代入(6)式，得到TM 波z向分量的本征方程:

(8)

式中，

上撇′表示求导数，〈〉表示内积(〈f，g〉=f(z)g(z)dz)，a_εq表示特征向量，k_ρεq表示特征值，ε表示TM 模式，将(8)式中的下角标μ 替换为ε可得TE 模式的本征方程.

假设在介质的电参数仅随z方向变化，与ρ 方向无关，所以电通密度可以写成:

(9)

式中，b_εq(ρ)表示第q个本征模式ψ_εq(z)的传播幅度.

将(9)式代入(4)式化简成贝塞尔方程形式^[11]:

(10)

(11)

通过电磁对偶性，同理化简(5)式可以得出TE 模式磁场的贝塞尔方程形式:

(12)

(13)

式中，a_μq表示TE 模式的特征向量，k_ρμq表示TE模式的特征值，b_μq(ρ)表示第q个本征模式ψ_μq(z)的传播幅度.

考虑均匀介质条件下，传播幅度在径向上可以写成贝塞尔函数或者汉科尔函数形式^[11]，

(14)

H_ν⁽¹⁾ 表示ν 阶第一类汉科尔函数.

通过(10)式和(12)式，径向上满足源项的连续性条件可以求得X_ενq^J，X_ενq^H，S_μνq^J和S_μνq^H，

由此可以得到z方向上的电通密度d_zν 和磁通密度b_zν:

(15)

式中，

T 表示矩阵转置.

2.2 电磁波的横向分量

由麦克斯韦方程式，可得φ 方向电磁场与z方向电磁场的关系(因为接收线圈不受ρ 方向上的电磁场的影响，所以本文仅考虑横向分量中φ 方向上的电磁场)^[24-25]:

(16)

(17)

电场的φ 向分量e_φν 和磁场的φ 向分量h_φν 可以表示:

(18)

式中，

2.3 反射矩阵和透射矩阵

由于电磁波在径向界面上反射和透射，源区内的电磁波z向分量和φ 向分量可表示^[11]:

(19)

(20)

下角标“0"表示发射层(发射层和接收层在同一径向分层中)，为从0 层(源层)向-1 层的广义反射系数矩阵，为从0层(源层)向1层的广义反射系数矩阵，可以通过边界条件递归来求得和，边界条件是在最内层边界(钻铤表面)ρ_min 处存在全反射，由此计算电场的反射矩阵元素为R_min，min+1^- =，磁场的反射矩阵元素为R_min，min+1^- =.在最外层边界ρmax 处不存在反射波，由此计算反射矩阵元素为=0^[25-26].广义反射系数矩阵表示为

(21)

(22)

此处下角标n表示第n个径向分层，R_n，n+1^-，R_n+1，n⁺表示界面ρn上狭义反射系数矩阵，T_n+1，n⁺，T_n，n+1^- 表示界面ρ_n上狭义透射系数矩阵，I表示2M×2M的单位矩阵，如图 2.

(23)

(24)

(25)

(26)

式中，Ψ_n，n+1 =a_n^T·p_n+1·a_n+1，“+/-"分别表示外行波和内行波.

W₀^- 和W₀⁺ 为源层的幅度系数，可以通过径向上源项的连续性条件可得^[11]，

(27)

(28)

2.4 接收线圈的感应电压及视电阻率

倾斜接收线圈的等效电阻表示为(具体推导公式见附录)^[11]，

(29)

式中，

ρ_R 表示接收线圈在水平面上的投影半径，z_R 表示接收线圈的z向坐标位置，θ_R 表示法向角度，φ_R 表示方位角度(本文均取0°，如图 1)，下角标“r"表示接收线圈所在层(一般考虑与源层在同一径向分层中，即“0"层)，*是指仅与右项矩阵的前半部分相乘，表示仅考虑电场贡献.

2.5 基函数

基函数的选取很大程度上影响了数值模式匹配法的速度和精度，本文采用幅度基函数和斜度基函数，考虑z方向的单元(i)(单元的两个结点zi和z_i+1 )(见图 3)，基函数形式表示如下^[20]:

其中，.

因为单元内的4个基函数Λ₁，Λ₂，Λ₃ 和Λ₄ 满足如下条件:

所以幅度基函数和斜度基函数保证了z方向上的σ_zE_z和μ_zH_z连续及其导数连续，而样条插值函数仅能够满足z方向上的σ_zE_z和μ_zH_z连续，但是不能满足σ_zE_z和μ_zH_z的导数连续，为了达到一定的精度需要更多的基函数^[11].采用本文的基函数，可以降低矩阵的维数，因此运算速度也相应提高了(大约可以提高1倍).

3 数值算例及其结果分析

本节采用改进的数值模式匹配法模拟倾斜线圈系在横向各向同性介质中的响应，并对比分析常规水平线圈系和倾斜线圈系在横向各向同性介质中的响应特征.

3.1 算法验证

为了验证数值模式匹配法的精度，对无限厚各向同性地层，分别用模式匹配方法和解析解方法计算了不同法向角度的线圈系响应，倾斜线圈系的参数结构见图 4.无限厚地层电导率为0.001S/m 到10S/m.仪器发射频率为2 MHz，发射电流(I_T )为1A，仪器其它参数见表 1.

(30)

由图 5(a，b)可以看出，无论在有铤还是无铤条件下，采用数值模式匹配法模拟倾斜线圈系在无限厚均匀地层响应值与解析解结果吻合较好.

图 6为两层横向各向同性介质地层模型，上层介质水平电导率和垂直电导率分别为1.0S/m 和5.0S/m，下层介质水平电导率和垂直电导率分别为5.0S/m 和1.0S/m，井半径为0.1270 m，井内泥浆电导率为0.0005S/m，线圈法向角度为45°，线圈距为0.762m(30in)，其它仪器参数同表 1.分别采用数值模式匹配法和有限差分方法模拟倾斜线圈系在该地层模型中的响应，对比结果见图 7(a，b)，有限差分结果取自文献^[11].

由图 7(a，b)可以看出，在考虑横向各向同性介质时，数值模式匹配法与有限差分法结果吻合较好，但由于数值模式匹配法是一种半解析半数值解，在运算过程中仅考虑一个方向上的有限元求解，其运算速度比有限差分方法要提高很多.

3.2 钻铤对倾斜线圈系响应的影响

为了考察钻铤对倾斜线圈系响应的影响，本文分别计算两种不同倾角的仪器在有铤和无铤条件下响应信号的相位角之差和幅度值之差，仪器参数见表 2.

图 8a和图 9a为有铤或者无铤条件下，仪器a，b在无限厚地层中的相位值和幅度值信号.总体上钻铤没有改变相位和幅度衰减的变化趋势，但对二者的影响有比较大的差异，为了更清晰地表示钻铤对不同线圈系接收信号的影响，对同一种仪器将有钻挺和无钻挺时对应接收器接收的信号差显示在图 8b和图 9b中.由图 8b可以看出钻铤对响应信号的相位值存在影响，对仪器a而言，钻铤影响随着地层电导率增大而增大，当地层电导率等于10S/m时，两者差异大约10°.对仪器b 而言，钻铤对相位影响总体较小，但二者的差异是非单调的.由图 9b可以看出钻铤的存在使得响应信号的幅度值有较大的影响，总体降低了大约30dB，而且影响随着地层电导率增大而有所增大.总之，钻铤对响应信号的影响比较复杂，不仅随地层电导率变化，并且和线圈的法向角度有关.

3.3 不同法向角度的倾斜线圈系在横向各向同性地层中的响应

为了对比倾斜线圈系与水平线圈系在横向各向同性地层中的探测特性，在图 10中的三层地层模型分别模拟了两种仪器的响应.倾斜线圈系的参数结构见图 4，两种仪器的源距，线圈间距，钻铤半径及频率相同.模型中间层为横向各向同性介质，水平电导率和垂直电导率分别为1.0S/m 和0.2S/m，层厚为0.5m，上下两层分别为各向同性介质，其电导率为1S/m，井半径为0.1270 m，井内泥浆电导率为1S/m.图 11(a，b)分别为两种线圈系相位差和幅度比的对比结果.

由图 11(a，b)可看出，在水平分层垂直井条件下，即使在横向各向同性介质中，水平线圈系响应信号的相位差和幅度比仍然呈一条直线，仅反映出水平电导率的变化，对垂直电导率无响应，而倾斜线圈系能有效地反映出地层各向异性特征，即其响应不仅受到水平电导率影响，也受到垂直电导率影响.

为了进一步考察不同法向角度的线圈组合在各向异性地层中的响应特征，分别给出9 种不同倾斜角的线圈系在上述地层中的响应，模拟时两个接收线圈与发射线圈中心距离分别为0.6096 m(24in)和0.762m(30in).地层模型如图 11a所示.

图 12，13中括号的数值依次表示发射线圈、短源距接收线圈和长源距接收线圈的法向角度.由图 12，13可以看出，在垂直井中，无论相位值差信号还是幅度值比信号，都只有在线圈组里存在一对倾斜的发射接收线圈才能对横向各向同性地层的电各向异性敏感，同时看出，不同的组合对地层各向异性响应也不同，组合不同线圈系组合有可能反演出地层垂直和水平电阻率值.由于线圈组合的种类非常之多，相应的探测特性也会不尽相同，限于篇幅，不再详述.

4 结论和讨论

本文采用垂向数值模式匹配法计算了倾斜线圈系在横向各向同性地层中的响应特征，计算过程中通过引入幅度基函数和斜度基函数，明显提高了计算速度(大约1 倍计算速度).通过计算结果发现:钻铤的存在会给响应结果带来一定的影响，且此类影响并非恒定的，是与地层电导率信息和线圈法向角度有关的，所以对于倾斜线圈系而言，当考虑不同法向角度的线圈组合时，必须对钻铤的影响进行校正.本文还在水平分层垂直井的条件下，对比了不同倾斜线圈系在各向同性和横向各向同性介质中的响应，可以看出不同线圈系的探测特性和对各向异性的敏感性都有较大的差异，更进一步的模拟表明发射线圈、两个接收线圈的法向角度分别为45°、0°和45°时有更好的适应性.然而，当前倾斜线圈系主要应用在地质导向中，用以探测地层边界，在地层的定量评价中相对较少，这可能是未来倾斜线圈系研究和应用的重要发展方向，本文的研究还不完善，也不深入，未来还需要开展更多的正演模拟和参数反演研究.

附 录

这里推导接收线圈的感应电动势，得出视电阻率，文章所涉及到的问题仅考虑接收线圈的电场变化，所以采用“*”表示乘以右侧矩阵的前半部分.

发射线圈的电流也可以通过类似的推导进行化简.

致谢

作者对两位匿名评审专家认真细致的审阅和富有建设性的意见和建议表示感谢，同时第一作者也要对张海澜研究员和何晓博士在数值模拟方面的指导表示感谢，感谢首席科学家肖加奇老师在文章修改方面给予的建议.