基于三维电测量的复杂地质构造识别及参数估计方法*

2020-09-10毕鹏程车爱兰冯少孔

工程地质学报 2020年4期

毕鹏程车爱兰冯少孔

（上海交通大学船舶海洋与建筑工程学院，上海200240，中国）

0 引言

高密度电法勘探效率高，采集数据丰富，勘探能力强，在工程地质勘察中得到广泛应用。传统的高密度电法是二维的，其研究与应用已经非常成熟，但是二维高密度电法的基本假设是电阻率仅沿测线方向以及深度方向上变化，而沿着与测线垂直的方向电阻率固定，即假定介质是二维的，因此针对某一方向上延伸的二维地质体（如断层）的探测效果比较理想（崔永圣等，2014；郑智杰，2017；底青云等，2018）。然而在实际工程中，探测的地质体通常呈现复杂的三维空间分布，此时二维电法存在局限性（胡静云等，2014；叶懿尉，2018）。三维高密度电法是在二维高密度电法的基础上发展而来的，它采集的信息大幅增加，经过反演后可反应空间整体的电阻率分布，因此一定程度上可以解决二维电法难以直观反映异常体空间形态的问题（高卫富等，2011；黄真萍，2015）。

目前三维高密度电法反演仍然处于研究阶段，反演速度及准确性仍需不断提高，并且实际工程中数据采集耗时长，其应用受到较大的限制。国内外对于三维高密度电法的研究主要以小范围的现场试验及反演方法为主。Dahlin et al.（2002）在瑞典地下污泥的勘探中应用了三维高密度电法，指出必须使用反演模型进行数据解析才能保证勘探的精度，表明了反演是实现三维高密度电法的关键。Loke et al.（2010）采用了49个电极进行了三维高密度电法试验，数据采用平滑约束最小二乘方法，大大缩短了计算时间。Wilkinson et al.（2015）通过将电阻率结构的估计值引入正演响应模型，提高了反演结果的精度。Alile et al.（2017）将三维高密度电法应用于地下填埋垃圾的勘探，初步探测出了填埋物的位置，但未能验证其准确性。

在国内，三维高密度电法在工程上的应用，大部分用二维数据反演，通过三维成图软件把视电阻剖面连续成三维视电阻率分布立体图，而同时跨测线及同测线的真三维测量并反演的实例很少（陈康等，2015）。张振勇（2015）将三维高密度电法应用于积水采空区的探测中，通过实例验证了三维高密度电法正演反演的可行性。李颜贵等（2009）在小场地进行了96道电极三维高密度电法工程勘查试验，发现三维高密度电法能直观地确定地表浅层的管线位置和埋深。张彬等（2011）将三维高密度电法应用于煤矿采空区的勘查，采用切片技术分析了拟建煤管站下方采空区分布情况，但其勘探精度有限。三维高密度电法的应用较少，对其结果的验证难度大，目前也尚未有研究分析其与二维高密度电法的结果差异。

本文首先介绍了三维高密度电法的基本原理，在研究其反演方法的基础上，以福建省寿宁县古银矿采空区勘探项目为背景，利用二维和三维高密度电法进行勘探，利用三维激光扫描对采空区的位置及形态进行精确测量，从而检验三维高密度电法在复杂地质构造勘探中的适用性以及其结果的准确性，提高其在实际工程中的勘探能力。

1 高密度电法原理及方法

1.1 基本工作原理

1.1.1 二维高密度电法

电阻率法勘探的原理基于地下不同岩土体的电阻率差异，首先将直流电源两极接入大地激发稳定的人工地电场，通过研究电流场在地下空间的分布规律，根据电阻率异常来判别异常体目标，从而实现地质构造调查的目的（查甫生等，2010）。

高密度电法采用阵列式电极进行电阻率观测，其本质仍然是电阻率法，它能反映地下空间剖面上岩土体电阻率分布，并且同时具有测深和剖面观测的能力。在野外进行高密度电法采集时，首先将若干电极沿测线以一定间隔垂直插入大地中，主机包括直流激电法仪和多路电极转换器，其作用是控制电极供电与测量，从而自动采集数据，二维高密度电法工作示意图如图1所示。高密度电法的实际应用过程中，其效果与各采集参数密切相关，主要包括装置类型，电极间距，电极数量等。这些采集参数直接影响着反演数据的优劣以及最终结果的精确性和分辨率。

图1 二维高密度电法示意图Fig.1 Diagram of 2D high-density resistivity method

1.1.2 三维高密度电法

三维高密度电法数据采集原理与二维方法相同，利用主机的多路电极转换器，用电缆链接所有的电极，利用电极进行供电与测量。但是其测线布设是把地表上的电极按网格布置，在网格点上插入电极，同时进行测线间和测线内的数据采集。由于数据量庞大，野外数据采集耗时长，因此需考虑到夜间工作的安全性以及采集效率的影响。图2为三维高密度电法工作示意图。

图2 三维高密度电法示意图Fig.2 Diagram of 3D high-density resistivity method

1.2 数据处理方法

1.2.1 二维高密度电法

高密度电法采集的数据电极测得的地下空间的视电阻率，可据此对地下的构造进行初步的分析。还需要对数据进一步处理并反演，得出地下空间剖面上的真实电阻率分布，从而对电阻率图像进一步解释。数据反演流程如图3所示。对于二维高密度电法，佐迪法以及最小二乘法是最常见的反演算法。

图3 高密度电法反演流程图Fig.3 Inversion of high-density resistivity method

1.2.2 三维高密度电法

三维电法反演流程与二维电法相同，但是由于计算过程更为复杂且数据量过于庞大，对计算机性能要求高，因此在修正模型时需要改进反演算法。目前主要的反演算法有基于Born近似的三维反演、层析成像反演、Tarantola反演以及最小二乘反演等方法（吴小平等，2000）。最小二乘迭代方法得到非常广泛、有效的应用，其中基于光滑约束最小二乘法的最优化方法提高了反演速度，同时也保证了解的稳定性，其原理基于如下方程。

式中：J为偏导数矩阵；J′为其转置矩阵；u为阻尼系数；d为模型参数修改矢量；g为残差矢量；f为各方向上的平滑滤波系数矩阵。

本文三维数据处理采用日本DIA公司开发的反演分析软件Etomo3d，该软件采用光滑约束的最小二乘法，约束因子自动选择。具体的数据处理步骤分为以下部分：

（1）预处理：首先对采集到的数据进行噪音剔除、平滑与插值等预处理，再转换成本软件可读取的数据格式。

（2）划分网格：根据现场的GPS定位信息编辑电极位置信息及地形信息，设置合理的网格大小参数，软件自动划分网格。本软件的特点在于可设置足够大的模型，减小边界效应的影响，在电极布设区域网格加密，越远离布设区域网格尺寸越大。

（3）反演分析：输入网格模型及转换好的数据，选取合适的参数（主要包括约束条件及类型、自由度、阻尼系数等）。首先计算电阻率分布初始模型理论电位，与测量电位进行比较差异，然后修正模型，通过多次迭代后逼近最优解。测量电阻值与模型参数关系为：

式中：φ为测量电阻值；p为模型参数向量；i为数据编号。代入初始模型参数，可得：

式中：p0为初期模型参数。因此观测值φobs与计算值 Фi（p）之差可表示为：

经过对数变换 θi＝ln pi及xi＝lnφi后可得观测方程式：Δx＝AΔθ，式中 Δx为对数化后的观测值与计算值之差，A为感度矩阵。为求可靠解，必须提供边界条件，本软件采用光滑约束的最小二乘法，在光滑约束条件下，观测方程式变为：

式中：ω为阻尼系数，采用AIC统计准则取最优值；Δθ为观测值与计算值之差。因此解得：Δθ＝（ATA＋ω2LTL）－1ATΔx。并根据Δθ调整模型，迭代计算求解。

2 现场试验

现场试验以福建省寿宁县溪乾村古银矿洞调查开发项目为背景。当地政府决定对溪乾村一带分布的古银矿洞采空区进行调查，以便开发成旅游景点。由于采空区洞内状况、延伸与走向等不详，加之山坡上到处都有洞口露出，直接施工存在很大安全隐患，因此需要通过勘探把握采空区古矿洞的内部状况、分布以及各洞口间的连通性。现场试验的主要目的在于检验三维高密度电法的适用性，并比较二维与三维高密度电法结果的差异，因此选取了本工程中采空区内部空间大，易于验证结果的部分来进行试验。

2.1 测区的地质资料

通过资料搜集、现场的初步踏勘及测绘，获取了采空区初步的地质资料，包括地质信息、地形图等，把握其总体的分布范围，并依此进行测区、测线的方案设计。

本项目位于华南活动带之东南沿海中生代火山断陷带中段，闽东火山断坳带之北段，中生代以来，区内岩浆活动频繁、强烈，其中火山活动更为突出，形成了分布广泛、厚度巨大的一套陆相沉积－火山岩系。根据前期现场踏勘，采空区古矿洞可进入洞长约60im，入口处洞高1.5im，宽约1im，洞内高度及宽度约4～6im，洞内状态如图4所示。洞室周边火山岩体破碎，存在较丰富节理裂隙。洞内有滴水及渗水现象，未见大面积积水。

图4 洞内状态Fig.4 Condition of the cave

2.2 高密度电法数据采集

图5 洞口位置Fig.5 Position of the entrances

采用GPS定位确定测区所有露头洞口的位置信息，得出其在地形图上的分布，如图5所示，标注出所有露头洞口的位置共3个，进而确定勘探区域的范围。另外在电法测线布设时，利用GPS对测线首尾及地形突变点定位，获取其大地坐标及高程，可将电法结果准确地与地形图信息结合。

采空区往往存在顶板塌陷导致内部通路堵塞等情况，因此可能存在很多难以进入或未知的地下空洞。由于洞室沿矿脉分布，并富含水，可期望利用洞室与周围岩体的电阻率差异推测出洞室分布范围。因此决定采用高密度电法进行勘探，根据前期初步勘测中高密度电法勘探的结果，表层电阻率在100Ω·m以下，山坡的地表电阻率随含水程度在数十～数百欧姆米范围内变化。火山碎屑岩的电阻率随风化程度而变化，中－强风化火山碎屑岩电阻率约在1000Ω·m以下，弱风化－未风化火山碎屑岩的电阻率约为数千欧姆米。

高密度电法数据采集系统由主机、多路电极转换器（电极连接盒）、电极系3部分组成，主机负责发送控制命令、电流和电压测量、接收信号等部分；多路电极转换器根据主机发送的命令进行电极转换；电缆由多芯电线组成，主要作用是信号传输。电极系由数量不等的接地电极组成，负责传送和接收电流信号（图6）。

图6 电法仪器Fig.6 Instrument of high-density resistivity method

由于洞室尺寸较小，大部分在2im内，洞室属于高阻异常体，具体分布深度未知，因此电法勘探采用二极装置（Pole-Pole法）进行采集，二极装置探测深度大，适应于高阻异常体探测。远电位电极布设在测线2ikm之外。三维电法采集时，同时测量4条测线（4条相邻的二维测线），即包括这些测线的测线内数据及测线间数据，因此每次采集布设240个电极，三维采集示意图如图7所示。为使二维及三维电法结果具有可对比性，在数据采集时两者的观测方式均采用二极装置，在测线内电极间距均设置为2im，供电电压为360iV，供电时间为100ims，断电时间为50ims。不同之处在于，由于三维观测同时测量4根测线，其测线间距为5im。

图7 三维高密度电法采集示意图Fig.7 Diagram of 3D high-density resistivity method

在采空区上方共布设8条测线，测线间距5im，测线布设如图8所示，定义水平面上沿测线为x方向，垂直测线为y方向，竖直为z方向。三维高密度电法的测量区域共进行两组，第1组包括Line1～Line4，第2组包括Line5～Line8。根据当地的地形地质特点，采集工作采取了以下措施保证数据的质量：在电极难以插入的岩石处，就近取土埋入电极，撒淡盐水使电极与岩体耦合；由于本项目的勘探区域主要在山坡，坡度较陡，电极布设时根据地形变化在突变点上采用GPS定位记录地形信息，在反演时输入。高密度电法测试中采用360iV电压供电以加大电流密度，数据采集过程中加强监控，一旦发现仪器出现非正常反应或者数据出现非正常变化现象，立即停止观测，直到排除故障后方可重新开始观测。通过采取以上若干措施，保证了物探现场采集资料的质量。通过处理数据和解译结果了解隧址区地层分布情况并布置前沿勘探测线进行现场踏勘，详细记录地形、地貌、地质变化，特别是断层露头、基岩出露等关键部位，为勘探结果的解释提供辅助和验证依据。

图8 高密度电法测线布设图Fig.8 Layout of survey lines

2.3 数据处理

2.3.1 二维高密度电法

目前国内外高密度电法数据二维反演方法研究基本成熟，本次二维高密度电法数据反演采用美国GEOTOMO公司的商用反演分析软件Res2idinv进行处理，这种算法最大的优点是可以调节阻尼系数和平滑滤波器以适应不同类型的资料，是国内外最主流的商用处理程序。图9展示了测线Line1和Line2的剖面结果。为便于将二维和三维电法结果进行比较，将8条测线的二维结果进行了伪三维处理，即将各剖面依次展布于空间上，利用插值法计算剖面间数据，从而得出伪三维效果图（图10）。

图9 二维高密度电法结果Fig.9 Result of 2D high-density resistivity method

2.3.2 三维高密度电法

本次三维高密度电法采用日本DIA公司开发的反演软件Etomo3id进行处理，此程序的优势在于采用基于光滑约束的最小二乘法进行反演，提高了解的稳定性，同时提高了计算效率。图11展示了三维电法反演后的空间结果。

图10 插值处理后的二维电法结果图Fig.10 Result of 2D high-density resistivity method after interpolation

图11 三维维高密度电法结果Fig.11 Result of 3D high-density resistivity method

图12 三维激光扫描仪主机Fig.12 3D scanner

3 勘探结果分析及验证

3.1 三维激光扫描

三维激光扫描的主要作用是可进入洞室轮廓进行精确的测量及三维建模，并结合GPS定位得出洞室的实际位置，据此结果验证电法结果的准确性，同时还可为之后本项目的设计和施工提供依据。设备采用瑞士徕卡三维脉冲式激光扫描仪ScanStation C10（图12），其具有高达50i000点／秒的扫描速度，可以进行360°×270°全视场角扫描，获取目标顶部，水平方向以及垂直方向的数据。测量距离精度为2imm，角度精度为12″。测量时首先在洞内初步勘探地形，根据需要确定扫描站数、中转基站的位置和控制标靶（用来匹配每站扫描的点）的个数和位置。然后安置三维激光扫描仪，调整好方向和水平。扫描仪预热后，设置好扫描参数（分辨率、曝光度和全景扫描等），扫描仪自动进行扫描。

采空区内部洞室空间较大，可进入洞长63im，由于洞室较长，需设置多个中转站，总共设置7站。

三维激光扫描数据处理是通过cyclone提供的坐标匹配功能，将各站的点云数据拼合成一个完整的测量目标点云模型。删除干扰的云点，通过分布框选点云，最终完成了对洞室的建模，并可对模型进行加工，生成剖面图，计算高程等，根据应用需求以AutoCAD和GIS的格式输出数据。原始点云包括了部分误差、错误和无关信息。这些误差信息的删减也是扫描数据处理的一个关键。影响三维激光扫描仪精度的要素较多，主要包括环境因素、激光信号的信噪比、激光信号的反射率、激光脉冲接受器的灵敏度、仪器和被测点间的距离、仪器和被测目标间的角度等。这些因素直接导致了误差的产生。通常情况下，三维激光扫描仪数据的模型精度显著高于单点的精度，可通过设立球靶、球状目标的坐标来改正点云的坐标，提高扫描仪模型的精度。图13为三维激光扫描数据处理后的结果图

图13 三维激光扫描结果图Fig.13 Result of 3D laser scanning

3.2 结果对比分析

为更直观比较二维与三维电法结果，对图10及图11的空间结果取临界值，仅考虑高于临界值的高阻区域，如图14所示，高阻区域可能解释为高阻岩体或洞室。可以发现，两者都显示在探测区域左侧存在明显的高阻区，根据现场踏勘结果认为，此处为高阻岩体。在探测区域中部，三维高密度电法探测出连续的高阻体，此处与可进入洞室的位置基本吻合；二维高密度电法对浅层的洞室的反映不甚明显，难以判断洞室的位置，而在深部存在大片的高阻区，可推断为高阻岩体。

图14 二维与三维高密度电法结果对比Fig.14 Comparison of the results between 2D and 3D high-density resistivity method

图15 结果对比（Line 6）Fig.15 Comparison of the results（Line 6）

经插值处理后的二维结果（伪三维）与真三维结果存在较大差别，主要体现在二维结果反映出深层的高阻岩体，而三维结果主要反映出浅层的高阻洞室区域。为进一步分析其原因，对三维结果进行切片处理，与同一位置的二维结果的切片进行对比。取测线Line6进行分析，两者的结果如图15所示。二维高密度电法在此剖面中的浅层探测出3处高阻区，三维电法结果探测出两处高阻，根据踏勘结果，均可推断其为洞室。但二维电法结果探测出多处假象，同时由于二维电法对深层的高阻岩体敏感性更高，经多次迭代后，浅层的高阻洞室区域被弱化。而三维高密度电法由于考虑了空间上的连续性，突出了浅层的高阻异常体，从而排除了深层的高阻干扰信息，降低了多解性。

三维激光扫描采用基站作为定位点，可得出洞室的所有轮廓点距离基站的距离，再结合GPS定位能准确获取地下洞室的空间位置信息，进行切片处理后可获取在任意剖面上实际的洞室位置、大小及形状，此结果精确可靠，可作为基准。将三维激光扫描确定的洞室位置与三维高密度电法推测的洞室位置信息叠加，进行对比，如图16所示，可以发现，三维高密度电法推测的洞室位置与实际洞室位置基本吻合。