井间电磁波CT不同初始模型成像效果对比

2022-12-29郑军

工程地球物理学报 2022年6期

郑军

(湖北煤炭地质物探测量队，湖北武汉 430200)

1 引言

井间电磁波 CT(Electromagnetic Computed Tomography)是20世纪80年代发展起来的地球物理探测技术[1]，以其抗干扰能力强，不受地形限制等优势，近年来在工程地质勘察中得到了广泛的应用[2]，大量的工程实例表明，电磁波 CT 探测技术在工程地质勘察中具有较好的应用效果和应用前景[3]。

图像重建是井间电磁波 CT技术中最关键的步骤，其理论基础为Radon 变换和 Radon 逆变换[4,5]，实质是借助于介质间的收发装置采集数据，按照数学和物理关系反演介质内部物理量的分布，然后借助各种数据处理方法，最终以图像形式呈现。图像重建实现的过程可以归结为求解一个大型的、稀疏的、常常是病态的线性方程组[6]，比较常用的方法是迭代算法，包括：代数重建法(Algebraic Reconstruction Technique, ART)，联合迭代法(Simultaneous Iterative Reconstruction Technique, SIRT)，正交分解最小二乘算法(Least Square QR Decomposition,LSQR)等[7]，以上这些方法均需要首先从初始模型开始，因此初始模型关系到反演的效率甚至成像的正确性。

通过收集资料来看，有关电磁波CT的理论研究多集中于各种反演算法的研究，对初始模型的研究很少，陈国金等[8]对井间地震层析成像中自动生成初始速度模型的方法进行了研究，提出了一种自动求取层析成像初始速度模型的方法；段心标等[9]对井间地震层析成像初始速度模型研究，提出了一种生成模型网格节点初始速度方法，均为地震层析成像中初始模型的相关研究[10]。参照地震层析成像中初始模型的给定方法，本文尝试以溶洞模型为例分别利用平均值法、水平测线值法以及BPT反投影法给出电磁波CT层析成像初始模型，通过层析成像将基于三种初始模型的反演结果进行对比，之后结合工程应用实测数据，验证实际成像结果与模拟成果的一致性，为今后层析成像的结果解释和判断提供借鉴和依据[11]。

2 电磁波CT成像原理

井间电磁波 CT 成像的基本思想，是在利用高频电磁波穿透两孔间介质得到的投影数据，通过反演方法来精确地界定地下目标体的吸收系数和电导率等属性，最终得到两个钻孔之间的高分辨率地质剖面图像，由此来判断地下目标体的属性[12]。直射线理论是目前电磁波CT成像中常用、实用的技术，本文将电磁波在地下介质的传播路径看作是直射线。

根据电磁波理论，电磁波传输过程中，耗散介质的初始场强E0(单位：V/m)与接收场强E(单位：V/m)关系为[13]：

(1)

其中，α为介质对电磁波吸收系数；r为收发点间距，单位为m；f(θ)为收发天线方向因子；θ为接收点处天线与电场方向夹角,单位为°;e为自然常数。

由上述理论，假设R是n条电磁波射线在m个网格单元中的射线元组成的m×n维矩阵；X是电磁波吸收系数组成的矩阵X=(α1，α2，α3，α4，…，αn)T；D为接收场强的观测值D=(d1,d2,d3,d4,…,dn)T；则所有观测点的发射和接收之间得到的线性方程组可以写成：RX=D；求解该线性方程组，就能够得到井间介质的吸收系数，进而得到吸收系数的等值线衰减图[16]。求解上述线型方程通常采用的方法为：首先给定吸收系数向量的初始值，求取系数矩阵，然后将吸收系数向量视为未知向量，其次用误差很小的迭代反演算法求解该线性方程组的解，即新的吸收系数向量，求出新的吸收系数向量对应的场强，根据此时的场强与已知的场强进行比较，修正吸收系数向量，并将修正后的吸收系数向量作为新的初始值；重复以上的过程，直至场强达到足够小，满足要求为止[17]。

3 模拟分析

3.1 原始模型建立

利用 Matlab 软件分别建立了包含溶洞的均匀介质和层状介质模型，模型尺寸大小均设置为 20 m×30 m，两个模型中均设置溶洞一处，整个模型剖面以1.0 m为单元网格进行剖分，共生成20×30=600个网格。图2(a)为均匀介质模型，整个模型剖面电磁吸收系数相同，取吸收系数为0.1 Np/m，相当于地下介质为完整的灰岩；溶洞设置为2 m×2 m正方形，中心与模型中心重合，取吸收系数为0.3 Np/m。图2(b)为层状介质模型，模型给出了一个二层介质模型，由上往下0～10 m，取吸收系数为0.3 Np/m，相当于第四系黏土层；10～30 m取吸收系数为0.1 Np/m，相当于地下为完整的灰岩；溶洞设置在模型中心位置，设置为2 m×2 m正方形，取吸收系数为0.3 Np/m。

图2 均匀介质及二层介质原始模型Fig.2 Original model of homogeneous medium and two-layer medium

3.2 初始模型设置

电磁波CT数据成像网格初始模型的设置常用的方法包括：平均值法、水平测线值法以及BPT(Back Projection Technique,BPT)反投影法。

平均值法是将模型剖面看做是一个均匀的背景场，给出一个固定的吸收系数值作为背景场值，该数值并非随意给定，而是取模型中所有网格的算术平均值；通过计算，得到上述两个原始模型的平均值分别为0.101 Np/m和0.201 Np/m，利用平均值法所得到初始模型如图3(a)和图3(d)所示。

水平测线值法是将水平射线经过的每一个网格即每一行网格的吸收系数看做是相等的，其数值大小为该行所有网格的平均值，于是模型剖面在每个不同的深度均有一个相应的吸收系数背景值，初始模型则为层状模型。图3(b)、图3(e)为两个原始模型的水平测线值初始模型。

BPT反投影法是电磁波CT反演基本算法之一，将每一个网格都看作是独立的，网格内的介质看成是均匀的，把每条射线的贡献平均分配给射线途经的每个网格,然后对每个网格求通过其内所有射线贡献的加权平均值。BPT反投影法不需要迭代,其速度只相当于其它迭代算法一次迭代所用的时间，但通过BPT 重建图像边界模糊，异常反映不清晰，通常将其重建的图像作为其他反演方法的初始值。图3(c)、图3(f)为两个原始模型的BPT反投影法初始模型。

图3 初始模型设置Fig.3 Initial model settings

3.3 反演成像

使用常用的联合迭代法(SIRT)对前述均匀介质(图2a)和二层介质(图2b)原始模型进行反演成像，反演的初始模型分别采用平均值法、水平测线值法、BPT反投影法生成的初始吸收系数模型，为了便于对比不同的反演初始模型对成像效果的影响，其它反演参数设置均保持一致。

反演成像过程中，先进行不同初始模型条件下反演误差的统计分析，统一进行300次迭代，图4(a)为原始均匀介质模型分别采用平均值法、水平测线值法、BPT反投影法生成的初始吸收系数模型进行反演的迭代残差曲线，由于模型较为简单，三种初始模型反演残差整体差别不大，初期迭代残差最高的为平均值法，迭代残差0.078，其次为BPT反投影法，迭代残差0.067，最小的为水平测线值法，迭代残差0.06；经过多次迭代最终迭代残差均趋于一致，平均值法反演在迭代次数达到90次之后、水平测线值法反演在迭代次数达到10次之后，以及BPT反投影法反演在迭代次数达到60次之后，其迭代残差均趋于稳定的0.056，90次迭代之后三者误差曲线基本重合。

图4(b)为原始二层介质模型采用相应初始模型反演结果的误差曲线，相对于均匀介质模型，三种初始模型反演误差整体略大，初期迭代残差最高的仍为平均值法，迭代残差2.106，其次为BPT反投影法，迭代残差0.505，最小的仍然为水平测线值法，迭代残差0.06；经过多次迭代最终迭代残差均仍然趋于一致，平均值法反演在迭代次数达到200次之后、水平测线值法反演在迭代次数达到30次之后，以及BPT反投影法反演在迭代次数达到150次之后，其迭代残差均趋于稳定的0.043，200次迭代之后三者误差曲线基本重合。

可以看出，采用不同初始模型在多次迭代后均能够得到收敛的反演结果，不同之处在于其初期迭代残差与收敛速度存在明显的区别；平均值法初期迭代残差最高，BPT反投影法次之，水平测线值法初期迭代残差最低，在设定相同的收敛误差阈值情况下，平均值法迭代次数最多，BPT反投影法次之，水平测线值法迭代次数最少，即收敛速度水平测线值法>BPT反投影法>平均值法。

图4 不同初始模型迭代误差曲线Fig.4 Iteration error curves of different initial models

3.4 成像效果对比

成像效果分析分两种情况，一是根据上述反演误差曲线，取三种模型反演误差曲线重合时的迭代次，作为最终成像的迭代参数，将三种模型完成相同迭代次数的反演结果进行对比；二是限定反演误差的阈值，将每种模型反演后达到设定迭代残差阈值的结果作为成像结果进行对比。

由前节迭代残差情况可以看出，原始均匀介质模型情况下，不同初始模型在反演迭代90次之后误差曲线基本重合且稳定，原始二层介质模型情况下，误差曲线则在迭代200次之后基本重合且稳定。图5(a)、图5(b)、图5(c)为均匀介质模型情况下，三种初始模型迭代90次成像结果。通过比对可以看出，反演误差趋于稳定的情况下，三者取相同迭代次数的成像结果整体形态一致，在纵向上均可以准确地反应出异常体的上下标高，但在横向和斜向上又都有不同程度的拉伸，与真实模型差别较大，相比之下水平测线值法在横向上拉伸程度略大于平均值法和BPT反投影法。图5(d)、图5(e)、图5(f)为二层介质模型情况下，三种初始模型迭代200次的成像结果。通过对比，三种方法在对上下两层不同介质界限的反应上准确、一致，对异常体也有较好的反应，水平方向上拉伸变形也较均匀介质模型小；从成像细节来看，平均值法成像异常规模及强度略小于水平测线值法及BPT反投影法，水平测线值法成像异常在横向上拉伸略大于平均值法和BPT反投影法。另外对比发现，成像结果吸收系数在数值上比真实模型偏低，但三种初始模型的成像结果又基本统一，因此数值强度上的偏差可能主要由反演算法造成，本文不进行讨论。

图5 不同初始模型相同迭代次数成像结果Fig.5 Imaging results of different initial models with the same iteration times

同理根据前节中迭代残差曲线，原始均匀介质模型情况下，三种初始模型误差曲线均趋于稳定的0.056。原始二层介质模型情况下，三种初始模型误差曲线均趋于稳定的0.044，因此分别设置原始均匀介质模型反演误差阈值为0.056，原始二层介质模型反演误差阈值为0.044。图6(a)、图6(b)、图6(c)为均匀介质模型情况下，迭代残差达到0.056时三种初始模型成像结果。不考虑迭代次数，三者反演误差达到同一阈值，纵向上均可以准确反应异常体的上下标高，但在横向和斜向上亦都有不同程度的拉伸，与真实模型相比，BPT反投影法变形最小，平均值法略大，水平测线值法变形最大。图6(d)、图6(e)、图6(f)为二层介质模型情况下，迭代残差达到0.044时三种初始模型成像结果。通过对比，三种方法在对上下两层不同介质界限的反应上准确、一致，在对异常体的反应上，平均值法与BPT反投影法基本一致，水平测线值法成像的异常体强度明显变低。

综合比较图5与图6，平均值法和水平测线值法在迭代残差稳定的情况下各自成像结果之间仍有较大的差别，高迭代次数的结果要优于低迭代次数的结果，而BPT反投影法基本保持一致，因此反演过程中平均值法和水平测线值法以高迭代次数作为约束条件的成像结果可能优于以迭代残差阈值作为约束条件的成像结果，但过高的迭代次数往往会产生结果的局部畸变，BPT反投影法不论采用迭代次数或者残差阈值做为约束条件均可取得较好的效果。

图6 不同初始模型相同迭代残差成像结果Fig.6 Imaging results of different initial models with the same iteration error

3.5 综合分析

综合以上两种情况，在反演迭代残差收敛的情况下，取相同的迭代次数，三种初始模型成像之间差别不大，细节上BPT反投影法略好于其它两种；由于三种初始模型收敛速度不同，平均值法收敛最慢，当平均值法迭代残差稳定时水平测线值值法及BPT反投影法迭代残差早已收敛，因此水平测线值法及BPT反投影法成像时如何确定合适的迭代次数是比较困难的。实际应用中，往往取固定的迭代残差阈值或误差变化率值作为反演终止的条件，通过上述试验可知以相同的迭代残差阈值作为约束条件时，三种初始模型成像之间差别较大，从成像效果来看BPT反投影法明显优于其它两种。

4 应用实例

通过以上模拟分析，对比了三种初始模型对电磁波CT反演成像效果的影响，接下来选择某岩溶勘查项目中一对电磁波CT钻孔来说明实际应用中的效果。

图7 实测电磁波CT反演断面Fig.7 CT inversion section of measured electromagnetic wave

该项目场区内地层主要有素填土、红黏土和灰岩，素填土厚度0.4～3.8 m，平均值2.3 m。红黏土层厚度为4.5～21.5 m，平均厚度13.8 m，地下水丰富，灰岩上部区域裂隙比较发育且裂隙之间彼此连通，裂隙中含有丰富的裂隙水，下部灰岩中有溶洞存在。

使用的设备为湖南奥成科技有限公司研发的HX-JDT-02B型井下无线电波透视仪[18,19]，采用扫频4～8 MHz 进行测量，选用4 MHz 频率采集数据进行处理。通过多次迭代，发现实测数据反演迭代残差约在0.500附近收敛，因此设置反演迭代残差阈值为0.500，成像结果如图7所示；图中浅蓝色竖线为地质勘察钻孔(一组三条竖线表示以该位置为中心，0.5 m为半径的范围内施工三个钻孔)，红色闭合曲线为地质勘察钻孔揭露溶洞的位置；从图中可以看出水平测线值法和BPT反投影法对异常的反应较平均值法更加清晰，水平测线值法异常在横向上存在拉伸变形现象，剖面下部两侧两个异常，在图7(b)中则反应为一个水平条带异常。总体来看，BPT反投影法成像效果优于水平测线值法，水平测线值法优于平均值法，与模拟分析结果一致。