基于电磁波层析成像技术的岩溶探测正演模拟及应用研究

陈春飞，沈晓武，张秉政

(浙江省工程勘察设计院集团有限公司 综合物探研究所，浙江 宁波 315012)

1 引 言

电磁波层析成像技术(CT技术)是将一定频率的电磁波作为发射源,利用接收机获取有效参数的地球物理方法,已经被普遍地运用于城市地下工程探测[1-3]。利用不同介质对电磁波吸收存在差异,造成在传播路径上其能量衰减的不均匀性，电磁波CT成像技术通过对实测数据进行反演计算可获取介质的吸收系数分布,从而得到地下的精细结构和性质差异图像[4-6]。

溶洞的形成是灰岩或碳酸盐地区地下水长期溶蚀的结果[7]。溶蚀的程度不同可能会形成破碎的溶洞及含水溶洞或者是空洞。这些对工程的施工与建设具有一定的威胁，可能会造成坍塌、涌水，然后导致地表开裂、下沉等等，不仅对工程的正常建设有着重大影响，还严重威胁施工人员的安全，造成重大的安全事故[8]。

因不同充填状态存在的空洞、破碎含泥沙溶洞以及充水溶洞等，其对电磁波的吸收存在一定的差异，本文将电磁波CT成像技术运用于溶洞探测。本文首先阐述了电磁波CT成像技术的基本原理与装置形式，然后构建不同填充溶洞模型，并进行正反演计算，讨论总结了电磁波CT成像技术在溶洞模型上的响应特征及规律。最后，通过野外工程实例，验证了电磁波CT成像技术在溶洞探测工作中的有效性，为以后野外实际工作提供更多的参考。

2 基本原理

2.1 电磁波衰减理论及层析原理

地下介质的不同物性分布对电磁波的作用主要表现在对电磁波能量的吸收，这种吸收作用与地下介质的裂隙和破碎带的分布、含水程度、矿物质的含量，以及不同的岩性分布等因素有关[9]。通过两个钻孔之间电磁波扫描观测，利用层析成象反演算法，将不同岩性导致的电磁波能量上的差异分布转变成二维介质分布图像，进而推断地下的地质结构情况[10,11]。

根据电磁波传播衰减理论，电磁波在介质中传播时将损失部分能量。衰减的大小与介质的导电性、密度等物性差异都有很大关系，也跟其本身的波长(频率)有关[12,13]。对于配置偶极子天线的电磁波仪，在不考虑电磁波绕行的情况下，接收到的电场强度，即在接收天线方向上的分量Ei可近似为：

(1)

式(1)中，E0为初始场强；r为发射源与接收点之间的距离；β为介质的电磁波吸收系数；θ则为接收天线与电场方向间夹角。

电磁波CT技术，其原理与医学CT相似。因而，电磁波CT就是根据实测数据，依照电磁波在地下有耗介质中传播规律及一定的物理和数学关系进行反演计算，得出被测区域内地质介质的电磁波物理参数的分布，最后以图像的形式表现出来。在实际工作中主要是用吸收系数这一参数来表征[14,15]。将式(1)变换，可得

(2)

式(2)中,ω为电磁波角频率；ε为介电常数；μ为介质的磁导系数；σ为介质电导率。

电磁波CT主要是在两相邻钻孔分别布置接收天线和发射天线，通过两个天线的上下移动，多次接收，得到若干条射线。每条射线记录了该点的电磁波场强，它能反映沿射线方向介质电磁波吸收情况。当测量区域中某点附近有数条射线通过时，即可通过公式求出该点的物性参数[16-21]。

电磁波CT观测系统如图1所示，在两钻孔间进行电磁波的层析成像，来探测两钻孔之间地质体的异常情况。图中ZK1为发射孔，ZK2为接收孔。不同深度按一定距离分别布置发射点和接收点。对数据进行层析化处理时，根据不同的发射源和接收装置，运用射线追踪原理，把成像区域离散为k×s网格，见图2。

图1 电磁波CT观测系统示意图

图2 成像区域网格化示意图

网格化后的单元总数为N(N=k×s),射线总数为M，第j个单元格的吸收系数为βj，j=1, 2, 3…,N。aij为第i条射线经过第j个单元格的长度。Pi为第i条射线的电磁波损失的能量,i=1, 2, 3, …,M。由此，可表示成矩阵形式为：

(3)

针对上述矩阵，用不同的数学方法进行层析处理，实现图像的重建。目前使用较多的是代数重建法(ART),联合迭代重建技术(SIRT)、最小二乘正交分解法(LSQR)等。

2.2 模型正演原理

正演是电磁波CT的基础，在反演后图像重建中有着非常重要的作用。因此，系统地对典型地质体进行正演模拟，选取一定方法对此反演，对比总结出其规律，为实际生产工作中能够更为有效地运用此方法提供一定的依据。

在模型计算中，为更好地与实际生产接轨，在假设介质的电磁性质连续的情况下，对公式(1)进行一定处理，把公式两端都取对数然后乘以20，把原电场强度的单位由V/m定义为分贝(dB)。经推导与处理，表示为：

根据公式(4)可知，在已知电磁波初始场强，介质吸收系数及接收机天线与场强方向夹角的情况下，求出观测到的场强值(为了成图方便，计算出的场强值使用绝对值处理)，运用射线理论进行层析分析，得到测试剖面的电磁波吸收系数等值线图。

2.3 反演方法及特点

电磁波CT的反演主要是求解一些大型线性方程组，使用方法也较多。如求解阻尼最小二乘法、加权最小二乘法和迭代类算法。前者的算法效果较好、精度高，但占计算机内存量很大，运算时间也较长。因此，现在应用较多的还是迭代类算法，如代数重建法(ART)、最小二乘正交分解法(LSQR)及联合迭代法(SIRT)等。

代数重建法(ART)是迭代法中较经典的算法，计算时所要求的内存少，但它的迭代收敛效果较差，对初值的选择依赖性也较大。最小二乘正交分解法(LSQR)是一种共轭梯度法，计算时所需要的内存比ART少，对于病态问题时效果也不错，但当数据误差较大时往往容易发散，一般需加入阻尼因子。联合迭代法(SIRT)在电磁波层析成像中的应用较为普遍，其算法是对ART算法收敛性的改善，内存所需量虽然要比ART多，但无论线性方程组具有超定性或欠定性，都可使用联合迭代法(SIRT)，而且收敛效果不错。因此，本文选用联合迭代法(SIRT)作为本次研究的反演方法。

3 正演模型研究

模型两孔间距为25 m，测量深度为4 m到35 m，测量的点间距为1 m，本次选取单层模型研究，完整基岩的吸收系数假定0.4 dB/m，溶洞半径为3 m，初始电场强度为-5 dB,定点发射深度则为5 m，10 m，15 m，20 m，25 m，30 m，35 m。为保证成像质量，需对发射孔和接收孔进行调换(定点接收5 m，10 m，15 m,20 m，25 m，30 m，35 m)，进行同样的测量。

首先对空洞模型进行研究，因空洞里一般为空气，电磁波衰减比完整基岩还小，所以其吸收系数比基岩还小，因此模拟吸收系数假定0.2 dB/m，模型图与正演后的电场曲线如图3所示。

图3 溶洞模型与正演电场曲线

据图3(b)可看出，曲线1～7(1～7分别代表发射深度5 m，10 m，15 m，20 m，25 m，30 m，35 m时的电场曲线)在不同深度都有极值点，说明存在异常，但电场的变化量较小，其它地方变化亦较为平缓，主要是路径衰减。反演处理后结果如图4所示。

图4 反演后视吸收系数等值线

图4中显示，异常体的形态和埋深与模型较为吻合，但异常体的异常边界不够明显，不易圈定，效果不甚理想。

接着研究破碎含泥沙溶洞和含水溶洞，水的吸收系数假定1.5 dB/m，破碎含泥沙溶洞则假定0.8 dB/m。模型与电场变化曲线分别如图5与图6所示。

图5 含水溶洞模型与正演电场曲线

图6 含泥沙溶洞模型与正演电场曲线

图5与图6中电场曲线变化都较大，极值点明显，且含水溶洞曲线峰值更高。反演结果如，图7所示。

从图7可知：异常体可根据等值线的走势和密集程度基本可圈定出来，无论是埋深、大小及形态都与模型较为吻合，且后者较前者效果更佳。

图7 破碎含泥沙溶洞与含水溶洞反演后视吸收系数等值线

综合以上三种模型可知：无论溶洞里含水或破碎含泥沙，或是空洞，电磁波CT均能反映出其异常，但值得一提的是，空洞因吸收系数与完整围岩较接近，因此效果不是非常理想，不能较为准确地判断出其空间形态及物性参数，最好能结合其他地球物理方法综合判断为宜。

4 工程应用实例

4.1 工程及工区概况

工程为拟建的中国石油云南1 000万吨/年炼油项目，该工程位于云南省安宁市，厂区属滇中高原构造侵蚀溶蚀切割地貌之丘陵地带，岩溶裂隙等较为发育。为了建设工程安全，需要探明场地内的地下岩溶分布状态。

4.1.1 地质条件

场地范围地形起伏较大，场地内北部、西部及南部边缘地貌多为丘陵，冲沟发育，中部及东部比较平缓，地表植被较发育，丘陵地带为树林，平缓地带为农田。地层岩性主要分布有寒武系中谊村组白云岩、含磷白云岩、泥质白云岩、石英砂岩、硅质岩等；筇竹寺组砂岩、泥岩、炭质泥岩、炭质砂岩等；泥盆系中统海口组、上统宰割组的白云岩、灰岩、石英砂岩、硅质岩、泥岩、砂岩等；二叠系下统倒石头组的泥岩、砂岩等。“昆明隆起”缺失中、后寒武纪和奥陶纪、志留纪、石炭纪的地层沉积。上覆第四系黏性土、砂土、含砾黏土、含有机质黏土等。岩溶发育区分布在泥盆系中统海口组、上统宰割组的白云岩、灰岩地层。

4.1.2 水文地质条件

根据前期钻探揭露，场地地下水类型主要为含有机质黏土、黏性土、全风化基岩等层中的上层滞水。上层滞水主要补给方式为大气降水及附近水体径流补给，排泄方式以大气蒸发及向邻区径流排泄为主；基岩裂隙水不发育。钻探期间，场地内上覆第四系的滞水、潜水水位及基岩裂隙水水位变化幅度较大。上层滞水埋深为1.3～5.0 m，分布不匀，大部分地段已枯竭；基岩裂隙水埋深在5.0～60.0 m，主要受季节性降水及场地起伏、构造裂隙影响，据初勘钻孔及物探验证钻孔揭露的溶洞资料，溶洞中大部分未见地下水，故场地地下水没有统一稳定的潜水水位。

4.2 地球物理特征

电磁波跨孔CT层析成像是根据电磁波在介质中的传播特性来进行探测，获得的层析结果是介质吸收系数β值的二维分布图像，据此推断地下的异常、构造等。吸收系数与地下介质的磁导率、电阻率、介电常数、发射电磁波频率等均有关，当层析介质不均匀或有缝洞破碎等异常时，其电阻率、介电常数、磁导率等均会发生变化。使用介电常数与电阻率测试仪(HF657EVF)对工区钻孔岩芯进行测试，获得场地的介质物性参数如表1所示。

表1 拟建场地介质物性参数

4.3 数据采集及解释

4.3.1 工作方法

工作方式采用跨钻孔采集，电磁波CT工作原理如图8所示。采用定点发射方式，点距为0.5 m，定收发点的间距定为2 m。选用的频率根据收发孔间距的不同，选用4 MHz至16 MHz，一般孔距越大，选用频率越小。

图8 电磁波CT工作原理示意图

4.3.2 数据处理解释

根据前期对工区已知钻孔异常点的探测及验证并结合以往工作经验，总结、归纳出本工区物探异常特征及解释原则如下：

1)第四系土层的电磁波吸收系数一般大于0.7 dB/m；

2)破碎含泥沙的白云岩、灰岩的电磁波吸收系数一般在0.4～0.7 dB/m之间；

3)含水溶洞一般吸收系数在0.7 dB/m以上；

4)空洞吸收系小于围岩，一般在0.4 dB/m以下；

5)根据视吸收系数剖面等值线的走势、疏密程度、封闭情况进行综合圈定。

因探测剖面太多，本文只选取其中典型剖面来解释与研究。选取26#剖面，测孔分别为30和28两孔，孔间距为20 m，测量深度为10～30 m。实际测量电场曲线如图9所示。

图9 26#剖面实际观测电场曲线

图9中曲线1～7分别为定发点在深度12 m，14 m，16 m，18 m，20 m，22 m，24 m的电场曲线图。如图所示，深度19～21 m处有极值点，曲线上下分界效果较为明显，上层曲线电场值随深度变化大，说明吸收较强，衰减大。部分曲线在24～26 m也有较大变化，说明也存在高吸收介质。

电场数据经过反演得到其视吸收系数等值线图及其地质钻孔解译剖面图如图10和图11所示。从图10可看出：剖面0～20.0 m点位，约-17.0～-19.0 m深度以上视吸收系数在0.75 dB/m以上，推测为第四系覆盖层；剖面4.0～9.0 m点位，约-23.0～-26.0 m深度，在16.0～20.0 m点位，约-18.5～-26.0 m深度视吸收系数在0.4 dB/m以下，结合等值线延伸趋势及钻孔情况，分析推测为不含水空洞；剖面其他区域视吸收系数在0.4～0.6 dB/m之间，推测为破碎基岩。以上对比钻探资料都较为吻合。

图10 26#剖面视吸收系数等值线

图11 26#剖面地质钻孔解译剖面

通过以上实际剖面资料可得出：资料中观测电场曲线对地层的分界效果良好，且能定性判断是否存在异常体；通过对剖面视吸收系数等值线图与钻孔资料对比可知，异常体的空间形态及位置均能较好地反映出来，说明电磁波CT技术在地层分界、特别是不同充填溶洞探测等方面均有良好的效果。

5 结 论

通过工程数据解释分析，得出以下结论：

1)无论溶洞里含水或者是空洞、破碎含泥沙等，电磁波CT方法均能反映出其异常，但因空洞与完整围岩的吸收系数较接近，因此效果不是非常理想，不能较为准确地判断出其空间形态及物性参数，最好能结合其他地球物理方法综合判断为好。

2)只要地层的吸收系数存在一定的差异，电磁波层CT即可较清晰、有效地反映出地层的分界情况。可根据观测电场曲线判断地层分界位置，地层的分界位置一般在曲线的极值点附近，而且可根据极值点上下曲线的斜率来定性判断地层吸收系数大小关系，从而进一步推断出其岩性；而根据反演所得电磁波视吸收系数等值线图则能更为直观地反映出其分层情况，分界线附近往往等值线密集，且走向一致。但若要比较精确地找出分层位置，结合观测电场曲线图与吸收系数等值线图综合分析往往效果更好。

致 谢

感谢中国地质大学(武汉)许顺芳副教授、范志雄高级工程师及刘嵘博士的帮助！