0 引言

井间电磁波层析成像是一种高分辨率的井中地球物理探测技术，它通过在钻孔中发射和接收电磁波，能够获得井间剖面的电磁参数分布.自1923年由前苏联科学家А.Д.佩特罗夫斯基提出以后，在矿产勘探、工程勘查及油气开发等领域取得了显著的应用效果.井间电磁波层析成像依据地层中不同介质对电磁波具有不同的吸收作用，通过测量电磁场强度得到吸收系数的分布，由此来判断目标地质体的位置与形态.随着钻探成本的提高，人们需要从钻孔中提取更多的信息来了解钻孔间的地质情况，井间电磁波层析成像改变了横向探测能力不足的弱点，为井间的地质探测提供了重要的技术支撑.

考虑到电磁波正演和反演计算的复杂性，人们常会采用一些近似的方法对数据进行解释.井间电磁波层析成像通常使用的高频电磁波(频率从0.1 MHz到50 MHz)，可以认为电磁波是沿着射线传播的，其原理同医学CT类似，自成功应用于地学领域以后(Dines and Lytle, 1979)，这种基于射线理论的成像方法成为了电磁波层析成像数据解释的主要手段(Thomson and Hinde, 1993; Yu et al., 1998; Korpisalo and Heikkinen, 2015)，然而井间电磁波层析成像与医学CT在尺度和场源的可预测性等方面存在着差异，使得数据的解释更为复杂.

井间电磁波层析成像包含了初始辐射场强E₀和方向性因子两个重要的辐射参数.在医学CT领域，场源的强度在计算过程中可作为一个参考值，因此可以将其视为常数，而在地球物理领域，初始辐射场强不仅与发射功率、天线的结构相关，而且还随着岩石的电磁参数而变化(孙建国，2002).因此确定初始辐射场强E₀是井间电磁波层析成像的关键步骤之一.在反演中一般将E₀设为常量，直接凭经验取值或者由“迭代线性拟合”(刘立振，1986)得到；也有采用“相对衰减层析成像”的方法(刘立振，1990)，回避直接求E₀.根据经验设置E₀值通常很难符合实际情况，而用“迭代线性拟合”求取E₀也常常得到不合理的结果.“相对衰减层析成像”回避了直接求E₀的问题，但是需要对背景场进行估算，很难保证方程之间相容.冯锐等介绍了一种计算感应波动场的正演公式，提出了计算原地辐射常数E₀值的交替追赶方法(冯锐等，1997)，魏明果提出了一种“电磁波井间场强幅值的相对归一化”方法，用以解决成像控制方程之间的相容性问题(魏明果和刘润泽，1999).曹俊兴研究了双频透射电磁波电导率层析成像方法，建立了包含两个频率电场强度的层析成像方程，克服了电磁波吸收系数层析成像的初始辐射场的计算问题，应用实例表明该方法可获得更精确的层析图像(曹俊兴，2001；Cao et al., 2003；王文娟等, 2009, 2011)，但该方法的只适用于电导率大于0.01 S·m^-1的高导介质，限制了其应用范围.Liu等(1998)在地震质心频移衰减层析成像(Quan and Harris, 1997)的基础上，提出了基于质心频移的电磁波层析成像方法，虽然该方法对几何扩散、仪器响应、天线耦合等影响不敏感，但为了能够得到频移及方差，必须使用宽频带的数据.Day-Lewis在2002年提出一种差分吸收系数成像的方法，通过测量同一位置流体盐度变化前后的场强变化值，求解得到吸收系数的变化值(Day-Lewis et al., 2002)，能够避免初始辐射场和天线效应的影响，但只能应用于物性随时间变化的区域.综合分析这些解决问题的方法，在实际应用中仍然存在局限性.

不仅如此，计算天线的方向性因子也是一个非常困难的问题，在以往工作中，通常采用的是均匀介质中半波偶极天线的方向性因子，但实际地质条件很难满足介质均匀的条件，而且钻孔附近的介质也会影响天线的辐射模式(Holliger et al., 2001)，使得半波偶极天线的方向性因子不再适合实际情况.此外，传统成像方法采用的代数重建、联合迭代、最小二乘等算法通常是不带约束的，往往会导致成像的分辨率降低.为此，本文首先采用数值模拟分析初始辐射强度和方向性因子的变化特征，建立估计辐射参数的井间电磁波层析成像方程，实现考虑钻孔约束的正则化反演，以期提高反演结果的准确性.

1 射线层析成像理论

根据电磁波辐射理论(吴以仁，1982)，假设电磁波沿直线传播，那么在无限均匀介质中偶极天线在远区的辐射电场强度E(r)可表示为

(1)

其中E₀是有效初始辐射场强(单位：dB)；β是电磁波的吸收系数(单位：dB·m^-1)；r表示发射天线中点到观测点的距离；f(θ)是偶极天线的方向性因子，它是同发射天线与观测射线的夹角θ有关.

对(1)式取对数，并规定参量，那么可以将(1)写成关于吸收系数的线性方程，即：

(2)

如果确定了初始辐射场强E₀及方向性因子f(θ)，通过求解方程(2)就可以得到吸收系数.通常假设E₀为常量，由线性拟合方式得到其估计值.井间电磁波测量采用的是半波偶极天线，因此方向性因子通常按照均匀介质中的公式计算.

2 数值模拟与分析 2.1 柱坐标下的时域有限差分法

为了评价上述射线成像理论有效性，并分析发射天线初始辐射强度和辐射方向性因子的变化特征，计算了部分典型模型的井间电磁响应.采用的数值模拟方法是柱坐标系下的时域有限差分法，将发射天线设置在轴对称中心线上，并认为周围介质具有轴对称的分布特征，把三维计算转化为二维，提高计算效率，同时正确反映偶极天线的辐射特征和电磁波的几何扩散特性(Holliger and Bergmann, 2002).

将麦克斯韦方程从三维笛卡儿坐标变换到二维圆柱坐标可以得到两组偏微分方程，TM(Transverse Magnetic)模式和TE(Transverse Electric)模式.二维TM模式只含有电场E_r、E_z和磁场H_φ三个分量，采用电偶极子作为发射和接收天线时，对应TM模式下的偏微分方程为

(3)

其中μ为磁导率，ε为介电常数，σ为电导率.

在柱坐标系下，发射天线位于对称轴上，因此可以被准确地模拟，而接收天线只能近似为无限小的电偶极子天线.将与天线表面相切的电场分量设置为零，模拟有限长度的天线，可用来分析不同长度天线的辐射特征.为了避免边界截断处非物理的电磁波反射，在模型空间的上、下、右边界采用Mur二阶吸收边界条件，只需要在边界按照近似的差分公式计算H_φ分量(姚军等, 2009).

2.2 均匀介质模型

为了验证柱坐标系下时域有限差分法计算的正确性，设计的模型是一个高80 m、半径为35 m的均匀介质模型(图 1a).模型介电常数ε=9ε₀，磁导率μ=μ₀，电导率σ=0.1 mS·m^-1.发射天线分别位于高程-20 m、-40 m和-60 m处，接收点位于距发射天线30 m的垂直测线上，接收点距为1 m.发射天线的工作频率为30 MHz，可以计算得到电磁波在介质中的波长λ=3.33 m，介质的吸收系数β=0.0545 dB·m^-1.采用半波偶极天线(长度2l=0.5λ)作为发射源，计算观测点处的E_z分量，然后同解析公式计算的结果进行比较(吴以仁，1982).

对比图 1b、图 1c和图 1d可以看出，发射天线在不同点时，解析法(图 1中蓝色实线)和时域有限差分法(图 1中红色虚线)的计算结果均非常吻合，说明了时域有限差分法计算的正确性.

为了分析不同长度发射天线的方向图特征，将发射天线置于图 1a模型中-40 m处，接收点位于以发射点为圆心，半径为30 m半圆上(图 2a所示).分别采用长度2l=0.5λ、1λ、1.5λ、2λ的偶极天线作为发射源，计算观测点处的辐射场强，然后对观测点处的场强归一化后得到偶极天线的辐射方向图(图 2b).

对比图 2b中的结果可见，天线的方向性因子和天线长度与波长的比值2l/λ有关，选择不同长度的发射天线，其辐射方向图也会随之改变.在相同的观测条件下，图 2c中由于2l/λ不同导致方向性因子与半波偶极天线的偏差最高可达20 dB.实际工作中很难保证观测条件满足半波偶极天线的假设，因此直接使用半波偶极天线的方向性因子对数据进行解释可能会出现偏差.

2.3 充填钻孔模型

浅层的钻孔通常是干孔，部分或者完全被井液充填.电磁波在空气中的衰减非常微弱，电磁波从空气进入介质中会因为反射等因素的影响，导致只有部分能量进入到介质中，但电磁波的反射对于空气充填钻孔中天线的方向图几乎没有影响；然而被井液充填的钻孔情况就大不相同了，电磁波在水中的传播速度是普通地层的1/3左右，使得电磁波的反射更加明显并以导波的形式沿钻孔方向传播(Holliger and Bergmann, 2002).

为了分析充填空气和充填井液的钻孔对天线辐射场的影响，在图 2a模型的观测条件下设计了钻孔被充填的模型.图 3a中围绕发射钻孔设计了宽10 cm的空气，设置其介电常数ε=1.0006ε₀，磁导率μ=μ₀，电导率σ=10^-20S·m^-1；图 3b中则模拟钻孔被井液充填的情况，将钻孔附近设置为高介电常数介质，其介电常数ε=80ε₀，磁导率μ=μ₀，电导率σ=100 mS·m^-1，发射天线的长度和工作频率和图 1a一致.

图 3c对比了均匀介质、空气充填钻孔和井液充填钻孔中偶极天线的方向图，图 3d中对比了它们三者中辐射场强的大小.可以看到，空气充填钻孔中的方向图和均匀介质中的基本一致，但图 3d显示空气充填钻孔中的辐射场强值比均匀介质中小20 dB以上，说明空气充填不会影响方向性因子，而会影响初始辐射场强的大小；井液充填钻孔中方向图与均匀介质中的相比，整体形态已经发生了改变，并且辐射场强幅值也低于均匀介质中的，说明井液充填会同时影响方向性因子和初始辐射场强.

2.4 水平层状模型

上述模型只考虑了钻孔周围介质是均匀的情况，然而多数情况下，地层的物性都是变化的.当电偶极子位于物性变化的界面附近时，辐射能量耦合到界面上，导致天线辐射的模式与均匀介质中的偶极天线完全不同(Holliger et al., 2001).通常，井间电磁波测量的工作模式是在钻孔内某一深度发射电磁波，而在另一钻孔不同点上接收，然后移动发射天线的位置，继续在不同点上接收电磁波场，直到观测射线将研究区域全部覆盖.在这种模式下进行等步长的测量，那么观测数据中常存在着一些具有相同发射角度的射线.

为了分析层状模型中这些观测点处场强的特征，设计了图 4a所示含两层介质的水平层状模型，上层介质的介电常数ε=4ε₀，磁导率μ=μ₀，电导率σ=2 mS·m^-1，下层介质的参数与图 1a中的模型保持一致.设置4个发射点的高程分别为-70 m、-60 m、-35 m和-25 m，接收点位置从-75~-5 m，点距1 m，发射点与接收点的径向距离30 m，发射天线参数与图 1a中一致，那么上层介质的吸收系数β=1.618 dB·m^-1.

图 4b和图 4c显示了不同发射点的模拟结果.按深度显示的结果中(图 4b)，4种观测方式在-40 m处均有明显的异常变化，指示了界面处物性变化对电磁波的影响，但这些异常之间没有明显的联系.将四个发射点的异常数据按照发射天线与射线夹角θ排序得到图 4c，为了避免其他因素的影响，只绘制了发射点和接收点在同一介质中的异常.图 4c中在-25 m定点发射的曲线(DF-25 m)和在-35 m定点发射的曲线(DF-35 m)，在-60 m定点发射的曲线(DF-60 m)和在-70 m定点发射的曲线(DF-70 m)均基本重合.可以看出，物性的变化改变了电磁波的辐射场，在不同介质中初始辐射场强不同，而在同一介质中发射和接收电磁波，初始辐射场强和方向性因子没有发生变化.

3 估计辐射参数的层析成像 3.1 层析成像方程的建立

通过以上分析，天线长度与波长的比值、钻孔的充填情况、钻孔周围介质的物性均会影响初始辐射场强或者方向性因子，导致初始辐射场强不再是常量，方向性因子不再和均匀介质中半波偶极天线的一致，这样给反演计算带来一些困难.所幸的是，井间电磁波测量在不同发射点存在着相同的发射角度，同时也较易由钻探编录获取钻孔周围介质的分层情况，结合上述分析，在同一介质中传播的电磁波初始辐射场强和方向性因子不会变化，若只考虑发射天线对初始辐射场强和方向性因子的影响，设发射天线在同一地层段，初始辐射场强和方向性因子均保持不变，并将它们作为未知数参与反演.那么根据介质的分层情况，设在第k段地层内发射电磁波对应的初始辐射场强为E₀(k)，方向性因子为f_k(θ)，在第k段地层内初始辐射场强不变，而方向因子是随发射角度θ变化的，那么可以将(2)式改写为

(4)

其中A′=－ln(E(r)r)，对于合理的离散化模型，可以将(4)式写为

(5)

其中E_k=－ln(E₀(k))，F_kⁿ=－ln(f_k(θ_n))，β_j是第j的模型单元的吸收系数，θ_n表示第n个发射角度，矩阵系数r_ij为第i条射线在第j个单元中的长度. A′_i和r_ij都可由已知的数据计算得到，那么可将(5)写成A=Gm形式的线性方程组，m=(β₁, β₂, …, E₁, E₂, …, F₁¹, F₁², …)^T，G是由系数r_ij和1构成的矩阵.

3.2 钻孔约束下的正则化反演

由于求解辐射参数增加了未知量的个数，因此采用带约束的正则化方法求解方程(5).钻孔的分层情况不仅在上述方程建立中起到作用，还能为反演添加物性约束条件，使得反演结果更加符合实际地质情况.假设钻孔附近第k段(k=1, 2, …, l)深度范围内的岩性不变，在反演中这段范围对应的模型吸收系数为(m₁, m₂, m₃, …, m_K)，那么这段范围内吸收系数应该全部相等，用数学公式表示这些模型的方差，即：

(6)

其中Φ_var(m)式模型方差目标函数，，将模型方差目标函数写成矩阵相乘的形式为

(7)

其中：

(8)

对模型方差目标函数求导可以得到：

(9)

其中.

对于(5)式中A=Gm形式的方程，考虑对模型的光滑约束和方差约束，采用正则化反演方法计算.将反演计算看作是以下的最优化问题：

(10)

式中，m_ref是参考模型，λ、γ为正则化因子，W_d是数据权重矩阵，W_m是模型约束矩阵，一般令W_d=diag(1/σ_j)，σ_i由观测数据质量计算得到的权重值，模型约束矩阵写为

(11)

其中α_x、α_z为模型约束的权重系数.W_x、W_z分别是模型m在x、z方向的差分矩阵.

在第n次迭代中，利用高斯-牛顿迭代计算模型的修正量，令，可以得到：

(12)

其中，对(12)式计算采用预条件共轭梯度法求解，构造预优矩阵M，使得，改善方程的条件数，提高共轭梯度法的收敛速度.

4 模型实验 4.1 水平圆柱体模型

在均匀介质中基于射线理论的电磁波正反演都是准确的，而当介质中存在电性变化较大的非均匀体时，电磁波除了透射，还会出现反射、折射、衍射等现象，此时射线理论将与实际情况出现偏差.为了分析射线层析成像理论和上述研究方法的有效性，在图 1均匀介质模型的基础上，设置一个中心坐标为(15 m, -40 m)，半径4 m的水平圆柱体模型(图 5a).其介电常数ε=6ε₀，磁导率μ=μ₀，电导率σ=10 mS·m^-1，发射天线依然是发射频率为30 MHz的半波偶极天线，可得到模型的吸收系数β=6.08 dB·m^-1.采用对调观测的模式，依次在2个钻孔中接收由另一钻孔发射的电磁波数据.分别利用射线理论和时域有限差分法计算观测点处的辐射场强，然后对正演数据采用传统射线成像方法和本文研究方法进行反演计算.

图 5b、图 5c是利用射线理论正演数据成像的结果，传统射线理论成像图(图 5b)和估计辐射参数的成像图(图 5c)均能准确定位高吸收异常体的位置和形态，异常体的吸收系数和理论模型几乎完全相同；图 5f中的红色曲线是初始辐射场强和方向性因子拟合结果的和，与理论值也非常吻合.可以看出，在射线正演理论的基础上，这两种方法都非常准确，并且估计辐射参数的层析成像技术还能够同时准确计算出初始辐射场强和方向性因子，由此也说明了正则化反演的有效性.

图 5d、图 5e是利用时域有限差分法模拟数据成像的结果，传统射线理论成像图(图 5d)和估计辐射参数的成像图(图 5e)能够准确显示异常体的位置，但异常体形态相对变小，传统射线理论的数值与理论模型的偏差达到4 dB·m^-1；图 5f中初始辐射场强和方向性因子的拟合结果(绿色曲线)和理论值还存在一定差异，但其成像结果更接近理论模型.以上对比可知，在非均匀介质中射线理论成像图由于射线理论近似的影响，导致成像的分辨率和结果数值降低，而估计辐射参数的层析成像技术通过增加未知量从一定程度上改善了这一问题，获得了准确的结果.

为了分析方向性因子变化时估计辐射参数层析成像的效果，将图 5a模型的发射天线长度设为之前的3倍(长度2l=5.0 m)，然后对比模拟数据的成像效果.传统射线理论成像图(图 6a)和估计辐射参数成像图(图 6b)的差异非常明显的，前者无法显示出井间的异常体，而估计辐射参数的成像结果依然能准确显示异常体的位置和形态.图 6c中反演得到的初始辐射场强和方向性因子不再和图 5f一致，而是接近长偶极天线(2l=1.5λ)的值，说明估计辐射参数的层析成像技术能适应方向性因子的变化以得到准确的成像图.

4.2 复杂层状模型

为了进一步分析估计辐射参数层析成像技术的反演效果，在图 5a的观测条件下，采用相同的发射天线，设置了一个复杂的层状模型.设置两个中心坐标分别为(10 m, -30 m)、(22 m, -60 m)，半径为4 m的无限长水平圆柱体模型(图 7a中的红色区域)，介电常数ε=6ε₀，磁导率μ=μ₀，电导率σ=10 mS·m^-1，其吸收系数为6.08 dB·m^-1；在高程-40~-15 m和-15~0 m设置了两层高吸收的介质.上层介质的介电常数ε=30ε₀，磁导率μ=μ₀，电导率σ=10 mS·m^-1，吸收系数为2.97 dB·m^-1；下层介质的介电常数ε=4ε₀，磁导率μ=μ₀，电导率σ=2 mS·m^-1，吸收系数为1.62 dB·m^-1.

图 7b、图 7c分别是传统射线理论和估计辐射参数的成像图.对比这两组结果，均显示层状的地质特征，传统射线理论的成像图比较杂乱，存在很多虚假异常，不能准确显示高吸收的水平圆柱体；图 7c显示的结果相对平滑，并能够比较准确地显示两个水平圆柱体.在复杂介质中，射线理论的近似的影响导致图 7c中也存在虚假异常，并且反演数值低于理论值，可见在种情况下，相比传统射线成像方法而言，估计辐射参数的层析成像技术依然能够获得准确的结果，说明了该方法的有效性.

5 实例应用

在广西某市地铁详勘中，为了调查地下空间岩溶的发育情况，在水平距离为22.8 m的两个钻孔间开展了电磁波测量.采用定点发射的观测方式，发射天线位于ZK6钻孔20.5~45.5 m的深度范围内，接收天线位于ZK7钻孔20~45 m的深度范围内，发射点距和接收点距均为1 m，共采集了26组数据，观测射线基本覆盖了井间的目标区域.

选择信号较强的10 MHz数据进行处理，剔除小于-140 dB的数据，然后分别采用传统射线成像方法和本文研究方法进行反演计算.传统射线成像方法计算时，设置初始辐射场强E₀=-5 dB，采用联合迭代计算50次后数据拟合的均方误差为±1.7 dB，结果如图 8a所示，吸收系数分布图和钻探所揭露的剖面地质特征基本一致，但是对ZK6钻孔在深度36~38 m附近钻遇岩溶没有明显反映，并很少有一些连续的异常区域，在成像图的左上角还出现了局部的高值异常，难以和地质资料对应.

图 8b是估计辐射参数方法的成像图，从图中能够清楚地显示ZK6孔钻遇的岩溶，并能容易划分出几个连续性的异常区域.结合编录资料分析，深度在25~28 m之间的中等吸收异常区对应于粉质黏土层，其下方高吸收的带状异常区则对应于高吸收的强风化粉质黏土，这与图 8a中显示的2个独立高吸收异常区别明显；在深度31~35 m之间还存在一个倾斜的高吸收带状异常区，推测可能是粉质黏土充填的岩溶，并和ZK6孔钻遇的岩溶存在连通的可能性；其他相对低吸收的区域则与石灰岩、白云质灰岩对应.然而由于射线理论近似、噪声等因素的影响，结果中的一些局部细节并不可靠，但从反演结果与地质编录的吻合程度上看，估计辐射参数的层析成像技术优于传统的射线理论成像方法.

6 结论

本文采用柱坐标下的时域有限差分法计算了电偶极子天线的辐射场，通过分析均匀介质，充填钻孔和层状模型中初始辐射强度和方向性因子的变化特征，提出了估计辐射参数的井间电磁波正则化反演方法，通过对理论模型和实例数据进行解释和分析得出以下认识：

(1) 天线长度与波长的比值、钻孔的充填情况、钻孔周围介质的物性均会影响天线的初始辐射场强或方向性因子，采用不合理的初始辐射场强和方向性因子计算将会导致成像结果出现偏差.

(2) 传统射线成像方法由于简化波的传播特征，导致反演中出现分辨率降低的现象，而估计辐射参数的层析成像技术通过增加未知数，一定程度上改善了射线理论近似的影响，同时还能适应初始辐射场强和方向性因子的变化，从而获得准确的成像图.

(3) 通过模型和实例验证，估计辐射参数的层析成像技术能够结合钻孔资料等约束条件，获取相对传统射线成像方法更准确、合理的结果.同时，由于受到射线理论近似的影响，得到的辐射参数与实际存在偏差，成像结果中还有一些虚假的异常，因此需要在更准确的正演基础上研究高精度的电磁波层析成像方法.