随钻电磁波测井技术得到众多石油公司的认可和使用^[1,2,3,4,5]，国内外学者也对随钻电磁波数值模拟做了相关研究。1991年，L.C.Shen^[6]首次对随钻电磁波测井仪器响应进行了理论研究,得到了轴旋转对称地层的解析解；Q.Zhou等人^[7]研究了水平多层介质中随钻电磁波仪器响应的解析解及数值计算解；J.R.Lovell等人^[8]和T.Wang等人^[9]分别采用有限元法和时域有限差分法进行随钻电磁波仪器响应数值研究；Liu Qinghuo等人^[10]使用共轭梯度法来求解轴对称条件下不均匀介质层的散射解；孙向阳等人^[11]采用矢量有限元法模拟随钻测井仪的电磁响应；魏宝君等人^[12]采用混合法和递推矩阵算法模拟层状介质中的随钻电磁波仪器响应。这些研究已取得了一定的成果，但如何达到计算精度高、收敛速度快以及实现动态数值模拟，是目前各种数值模拟算法有待解决的问题。

自适应hp-FEM算法是一种研究工程结构问题的有效计算方法，适用于具有复杂边界形状或边界条件、含有复杂媒质的问题，在分析井周复杂地层时具有明显的优势^{[13,14,15,16,17,18,19]}。hp-FEM算法能对不同区域采取自适应细化策略，通过减小局部误差来提高精度，降低计算次数，可以显著地提高有限元计算分析的效率及可靠性。笔者利用hp-FEM算法结合HERMERS软件进行了建模和运算，根据电磁场理论建立数学模型，考虑边界条件，对所建立的二维旋转对称性地层模型进行了数值模拟，得到了随钻电磁波测井的幅值比和相位差响应曲线，刻度成更直观的电阻率曲线，并分析了仪器响应随频率和源距的变化规律，以期为仪器参数设计优化提供参考。

1 随钻电磁波测井原理和方法 1.1 时谐电磁场理论

在随钻电磁波测井中，基于时谐电磁场理论，利用麦克斯韦方程组进行求解，微分形式的麦克斯韦方程组可写成：

式中：H为磁场强度，A/m；E为电场强度，V/m；ω为发射源角频率(ω≠0)，rad/s；σ为媒介的电导率，S/m；ε为介电常数，F/m；μ为磁导率，H/m；ρ为空间电荷密度，C/m³；J^imp为发射线圈上外加电流密度，A/m²；i=－1。

可得到时谐麦克斯韦方程在求解域Ω上的波动方程：

又因为

则式(2)化简后得到时谐麦克斯韦方程在均质介质求解域Ω上满足的波动方程：

其中

令

式中：k为传播常数，也称波数；β为相位常数；α为衰减常数。

1.2 边界条件处理

在利用有限元方法进行随钻电磁波正演计算时，根据所处理问题的边界条件进而求解。

1.2.1 理想导体边界条件

由理论电磁场可知，任意导体分界面上，电场强度的切向分量和磁感应强度的法向分量都是连续的。换句话说，理想导电体表面不可能存在电场切向分量和磁场法向分量。因此，在理想导体边界求解域Γ₁上满足的边界条件为：

式中：B为磁感应强度，T；n为在发射线圈边界求解域Γ₂所满足边界条件下的任意试函数沿法线方向的矢量。

1.2.2 发射线圈边界条件

由于发射线圈模型中通入外加电流J^imp，为方便计算，根据场的等效原理，可将外加电流等效为面电流J^imp_s，来建立发射线圈模型。因此，在求解域Γ₂上满足的边界条件为：

1.2.3 非理想导体边界条件

与理想导体不同，非理想导体中存在电阻、损耗，即σ≠0，求解域Γ₃为非理想导体边界，故其边界条件为：

2 随钻电磁波数值模拟方法求解 2.1 时谐麦克斯韦方程的变分方程

随钻电磁波测井响应模拟过程是对高频情况下的矢量电场进行分析，因此采用矢量有限元法对H(curl)高频电磁场进行求解，从而确保各单元相邻表面之间电场切向的连续性和场的旋度零空间的正确表示。在H(curl)空间求解电场可以保证求解过程中不出现伪解，正演仿真结果准确可靠。H(curl)空间定义为：

式(9)、(10)、(11)结合其相应的封闭区域边界条件作为约束条件代入式(4)中，然后在方程中引入矢量试函数F，并在求解域Ω上对各项进行积分运算，即可得到在H(curl)空间下有限元变分问题的标准方程：

式中：M^imp为外加磁场强度，A/m;Γ₂为发射线圈边界条件；F为任意的一个试函数，且F∈H={E∈H(curl,Ω):n·E=0 on Ω}；F_t为在Γ₂边界条件下的任意试函数沿切线方向的矢量。

当外加磁场M^imp为0时，式(13)可写成一个双线性函数和一个线函数的组合，即：

利用有限元法，对式(13)进行离散化处理，建立了大型稀疏线性方程组：

式中：A为容量矩阵;C为施加条件数;x为需要求解的未知变量。

hp-FEM是以细分单元网格或提高插值函数阶次为手段的一种自适应算法，在求解区域较大的网格剖分问题上具有计算能力强、计算精度高、收敛速度好的特点，因此采用hp-FEM有限元法并结合UMFPACK求解器对式(15)进行求解。

2.2 hp-FEM算法求解过程

用初始网格T_h,p将Ω划分为多个互不重叠的单元k₁,k₂,…,k_n，其中每个单元的尺寸h₁,h₂,…,h_n＞0，对应的多项式阶次p₁,p₂,…,p_n≥1。那么，在求解空间H(curl)下，有限元求解子空间定义为：

式中：E_h,p为k_h,p中依据初始网格得到的电场E的近似解。

由H(curl)空间特性可知，此时空间H_h,p中的矢量场分布是不连续的，但是该矢量场的切向分量在整个求解区域各单元表面上保持连续变化，满足高频情况下求解随钻电磁波测井过程中电场分布的要求。通常情况下，自适应hp-FEM算法的求解步骤如图 1所示。

图 1中，T_h,p指初始网格分布，由多个互不重叠的单元组成；TOL是求解过程中最大容许误差；E_h,p是初始网格得到的电场的近似解；E_ref是计算电场的参考解；H_h,p是求解空间范围内的电场近似解的集合；K是选择的一个初始网格，K_i(i=1,2,…，n)是由K细化为多个子单元；ERR_i(i=1,2,…，n)指每个子单元误差，ERR是根据子单元误差求解的全局误差。

3 模拟算例 3.1 无限大均匀地层验证

随钻电磁波测井仪测量的幅值比和相位差数据与仪器附近的地层电导率存在着对应关系，因此用 幅值比和相位差数据来表征地层的电导率。

随钻电磁波测井仪器如图 2所示，发射线圈T的匝数为n_T，半径为a_T，其中通以交变电流I_T=I₀e^iwt；接收线圈R₁和R₂的匝数和半径相同，均为n_R和a_R，R₁和T的距离是L₁，R₂和T的距离是L₂。由于L₁a_T，为简化起见，把发射线圈看作磁偶极子处理。发射线圈T发射高频电磁波，经过地层到达接收线圈R₁和R₂，并记录接收线圈间感应电动势的幅值比E_ATT和相位差ΔΦ。这2个参数的优点为：一是幅值比和相位差为相对测量值，可以降低井眼和线圈尺寸的影响；二是无需去掉直耦信号，简化了仪器结构，降低了仪器实现的复杂程度。

根据式(4)，计算得到均匀各向同性介质中接收线圈间幅度值和相位差的解析式：

地层模型1是不考虑井眼和钻铤存在条件下的无限大均质各向同性地层，参数设置如下：圆柱状地层半径为300 m，高为500 m(将该模型假设为无限大地层模型)，发射频率2 MHz，源距L₁为1.0 m，L₂为1.2 m，地层电阻率为1 Ω·m，相对磁导率设为1.0。均质单一地层利用hp-FEM算法的计算结果如图 3 所示，显示了hp-FEM方法的网格划分和网格多项式阶次分布。从图 3可以看出，均质地层中的电场分布对称，网格细化对称，一定程度上验证了该方法和程序的正确性和稳定性。

图 4为均匀地层hp-FEM数值模拟结果相位差计算值与解析值^[13,18]的对比，两者能很好地吻合，进一步验证了算法的正确性。

由式(18)得到的解析值和由h-FEM、hp-FEM正演模拟得到的相位差计算值见表 1，可以看出2种方法得到的计算值和解析值很接近，hp-FEM的相对误差比h-FEM小一个数量级，进一步验证了自适应 hp-FEM 正演算法的准确性和可靠性。误差主要源于地层模型的尺寸(虽然很大，但仍达不到无限大情况)。对比2种方法可以看出，hp-FEM具有迭代次数更少、计算精度高、网格剖分少和计算速度快等优点，因此hp-FEM可以在随钻电磁波测井正演模拟过程中应用。

3.2 视电导率曲线转化

采用式(17)和式(18)，对无限大均质地层计算随钻电磁波测井仪器响应幅值比和相位差电导率的转换关系。考虑发射频率400 kHz和2 MHz，幅值比和相位差主要反映地层电导率的变化，故此处假设地层介电常数为常数。

根据hp-FEM算法数值正演结果，在2 MHz和400 kHz发射频率下，地层无倾斜，源距L₁为1.0 m，L₂为1.2 m，相对磁导率为1.0的条件下，相位差和幅值比在均匀介质中随电导率的变化关系如图 5所示。从图 5可以看出，在较低的电导率情况下，相位差和幅度比数值相应较小；相同的电导率条件下，2 MHz频率下的相位差比400 kHz的数值大，2条曲线整体变化趋势相同。另外，相位差和地层电导率在双对数坐标下具有比较好的线性关系，幅值比刻度曲线也有类似特点，因此，可用线性插值的方法对相位差和幅值比进行视电导率的转换。

3.3 简单模型算例

地层模型2由6层层状水平地层组成，每个地层的厚度均为5 m，电阻率依次为2，10，1，25，1和15 Ω·m。井眼半径0.1 m，钻井液电阻率为0.1 Ω·m。随钻测井仪器参数设定如下：接收线圈间距0.2 m，相对磁导率为1，仪器源距分别为0.5 m和1.0 m，发射频率分别为400 kHz和2 MHz。按照上述地层模型正演,得到其对应的相位差曲线和幅值比曲线，及其刻度对应的电阻率曲线，如图 6所示。

从图 6可以看出，相对于源距因素，频率对相位差的影响较大；2 MHz下的相位差的变化范围明显大于400 kHz下的相位差变化范围，而且在同一深度时，2 MHz对应的2种相位差的数值比较大。从图 6(c)可以看出，仪器源距越小，幅值比越大；在两地层的交界面处，相位差曲线存在明显的极化角现象；而且极化角处对应的地层深度与所设模型的地层界面相近；在高阻层的分界面深度为5，15和25 m处，仪器源距越小，尖角现象越明显。幅值比曲线(见图 6(c))与相位差曲线有相似的结果，在由低阻层进入高阻层时,幅值比曲线极化角现象十分突出，说明幅值比曲线在确定高阻含油地层具有较好的参考价值。

4 影响因素分析

采用控制变量法研究源距和频率对随钻电磁波测井仪器响应的影响。模型参数为：地层厚度12 m(下围岩厚度为4 m，目的层厚度为3 m，上围岩厚度为5 m)，均为水平层状地层，且电阻率依次为1，10和1 Ω·m，相对磁导率为1；井眼半径为0.1 m， 设钻井液电阻率为0.1 Ω·m；接收线圈间距为0.2 m，发射电流为0.005 A；迭代误差为2%。

4.1 仪器源距因素

设仪器发射频率2 MHz，源距依次设为0.25，0.50，0.75，1.00，1.25和1.50 m，考察不同源距条件下数值模拟得到的测井响应曲线，分析源距对随钻电阻率电磁波仪器响应的影响。

在保持仪器间距和发射频率固定的情况下，不同仪器源距条件下得到的相位差和幅值比曲线如图 7所示。

从图 7(b)可以看出：源距对幅值比曲线影响较大；源距越小，幅值比越大；当仪器穿越不同地层界面时，源距越小，极化角越明显；在高阻地层界面分界面处，仪器的源距越小，幅值比曲线的分辨能力越强。

然而，相对于幅值比曲线，相位差曲线能更好地反映地层界面的响应效果，更能反映地层模型的分层情况。从图 7(a)可以看出，源距越小，相位差越大；在两地层交界面处，仪器源距越小，极化角越明显。由低阻层过渡到高阻层时，相位差曲

4.2 仪器发射频率因素

设线圈源距为1.0 m，发射频率依次为400 kHz、1 MHz和2 MHz，考察不同发射频率条件下数值模拟得到的测井响应曲线，分析源距对于随钻电磁波仪器响应的影响。

在仪器源距保持固定的情况下，不同发射频率条件下得到的相位差和幅值比曲线如图 8所示。从图 8可以看出，相对于幅值比曲线，相位差曲线在反映地层界面信息方面的效果明显较好，能更清楚地反映地层模型的分层情况。

从图 8(a)可以看出，仪器频率越大，相位差变化幅度越大，得到的相位差也越大；在两地层交界面处，相位差曲线存在明显的极化角，极化角对应的地层深度与所设模型的地层界面相近；在高阻层的分界面处，2 MHz对应的相位差曲线越大，极化角现象也越明显。由低阻层过渡到高阻层时，相位差曲线中出现了一个明显的极化角，而由高阻层过渡到低阻层时，极化角不大，因此仪器在由低阻地层穿过高阻地层时，地层分辨能力更强。图 8(b)的幅值比曲线与相位差曲线有相似的结果。

5 结 论

1) 采用自适应hp-FEM算法建立随钻电磁波测井仪器响应数值模拟模型，在保证计算精度的前提下对自由度、计算时间以及迭代次数等参数进行了计算，得到了幅值比和相位差响应曲线。

2 ) 在不同的仪器结构和地质条件下，测井响应值会有很大的变化；相位差曲线能更好地反映地层界面的响应效果，更趋向于地层模型的分层情况；仪器源距越小，发射频率越高，地层分界面处极值角越明显，分辨能力更强；仪器在由低阻穿过高阻时，极化角更为明显，可以为地质导向决策提供信息。

3) 基于hp-FEM的数值模拟能够根据地层模型的实际情况和误差指示自动选择合适的细化方式和计算策略，通过减小局部误差来提高计算精度，降低了计算次数，提高了有限元计算分析的效率及可靠性。