1. 引言
中国岩溶地区分布广袤，面积达344万平方千米，约占国土面积的三分之一[1]。石灰岩在广东、广西、云南、贵州等南方省区呈广泛连续分布态势，孕育了复杂多样的溶洞系统[2]。岩溶发育的隐蔽性与空间形态的不规则性，给工程建设及运维带来严峻挑战，诸如岩溶水突涌[3]、地面塌陷[4]、溶洞垮塌[5]等灾害事故频发，造成重大人员伤亡与经济损失。因此，实现岩溶空间分布的精准探测，成为保障工程安全的关键前提。
当前的岩溶探测主要依赖地质钻探与地球物理探测两大技术体系。地质钻探虽能获取直观可靠的岩芯样本，但存在成本高昂、效率低下、信息离散等局限性，仅能反映钻孔点位的局部特征。地球物理探测技术凭借高效、经济、无损的优势，成为大范围勘察的首选，其细分技术中，井中物探虽具备高分辨率，但作业效率较低，难以满足大面积场地探测需求[6]。地面物探法中的浅层地震反射波法，通过分析人工激发地震波的反射信号，可获取地下介质波速分布，但在城市环境中，炸药震源审批严格，且沿海地区高水位与松散沉积物覆盖层导致锤击震源效果受限[7]；瞬变电磁法利用二次电磁场反演电阻率，易受城市电磁干扰，分辨率不足[8]；探地雷达法探测效率高、分辨率可达0.08~3米，但在高水位地层中信号衰减显著，且城市建筑物产生的电磁干扰严重影响数据解译。
高密度电法融合电阻率剖面与测深技术，通过采集视电阻率数据构建地下电阻率模型，对低阻异常体响应敏感，且具有抗干扰能力强、数据处理便捷、成本可控等特点，在煤矿采空区探测、滑坡隐患识别、水库渗漏排查等领域已广泛应用[9]。在岩溶探测方面，国内外学者已开展诸多探索：Fu [10]等人利用该方法在城市环境中成功定位含水溶洞，反演结果与钻探数据高度吻合；吴荣新等通过数值模拟揭示地形、充填类型及埋深对电阻率响应的影响；刘道涵结合高密度电法与钻探，实现岩溶地下水通道的电性结构成像；划分出岩溶塌陷高危地质区域；在半干旱地区识别出天坑与空洞，为灾害治理提供依据。然而，沿海地区地层结构复杂，高水位与厚层第四系覆盖层显著改变介质电性特征，系统研究该环境下岩溶的电阻率响应规律，对提升探测准确性具有重要意义。本研究以广州市花都区某园区建设项目为依托，采用高密度电法开展地下岩溶探测。通过电阻率剖面数据反演，圈定多处潜在岩溶发育区域；结合地质钻探验证，成功揭露两处大型充填溶洞，且反演所得基岩面与钻孔岩芯结果高度一致，有效验证了高密度电法在沿海城市岩溶探测中的可靠性与实用性。
2. 高密度电法原理
高密度电法作为地球物理探测的重要手段，其技术原理建立在地下介质电阻率差异的基础上。该方法通过向地下施加稳定电流场，同步分析电流在介质中的传导分布规律，从而实现对地质结构的探测。在数据的采集过程中，高密度电法沿着测线将一组电极等间距的一次性布设完成，再根据选定的点极阵列进行数据采集，获得地质视电阻率分布(见图1)。
Figure 1. High-Density electrical method detection principle schematic diagram
图1. 高密度电法探测原理示意图
高密度电阻率法中常用的四极阵列包括温纳装置、斯伦贝尔装置(又称施贝阵列)、偶极–偶极装置、微分装置及环形装置。各类装置因电极排列原理差异，具备不同的性能优势与适用局限，实际勘察中需结合具体探测目标(如探测深度、异常体分辨率需求、地质复杂度等)，选择单一或组合式电极排列方案。
在上述装置中，偶极–偶极装置的性能特性尤为鲜明：其横向分辨率突出，可精准识别地下横向延展的地质异常体(如小尺度断层、孤立溶洞等)；但短板同样显著——垂直分辨率相对薄弱，难以精细刻画地层垂向分层结构(如薄互层岩性界面、含水层顶底板)，且数据质量随电极距增大呈明显衰减趋势，这一特性直接制约了其有效探测深度(通常适用于埋深小于50 m的中浅部区域)，因此更适配中浅部横向地质异常的探测场景。
与偶极–偶极装置形成鲜明对比的是温纳装置与斯伦贝尔装置：二者均具备优异的垂直分辨率，可精准捕捉地层垂向结构的细微变化，且信噪比高、数据稳定性强，受地表电磁干扰、地形起伏等外界因素的影响较小；但二者的横向分辨率均弱于偶极–偶极装置，对横向小尺度异常体(如宽度小于10 m的小型溶洞、薄层破碎带)的识别能力存在局限。值得关注的是，在上述三种装置中，斯伦贝尔装置的探测深度优势最为显著，其有效探测范围更大(可覆盖埋深超100 m的区域)，因此更适用于深部地质目标(如基岩风化界面、深部区域性断层、深埋含水层)的勘察任务。
在数据采集阶段，技术人员需沿预设测线将多组电极按等间距(通常间距为2~10 m，具体根据探测精度需求与勘察区地质条件确定)标准一次性布设完成，随后依托高密度电阻率仪的自动化采集功能，按照预先设定的电极阵列模式(如温纳装置、偶极–偶极装置等)开展系统化数据采集作业，确保数据覆盖的完整性与连续性。
3. 正演模型的建立
本文使用COMSOL6.2 (试用版)建立正演模型，核心是围绕物理场精准匹配加参数真实化加网格适配，先明确正演目标并选择对应物理场模块(电流)，再通过自定义绘制或导入外部文件构建与真实场景尺寸一致的几何模型，为各部件分配匹配物理场的材料参数，接着设置模拟真实约束与激励的边界条件和载荷，然后根据物理场特性划分网格并做独立性分析，本文正演模型是以广州市花都区花山镇南部的L4真实高密度电法反演数据(图2)为例进行模型设计[11]。根据临近区域的勘察资料，研究区域下覆基岩为石炭系灰岩，岩溶发育，岩性较为复杂。主要岩土层包括：第四系填土层、冲洪积层(细砂、中砂、粗砂、砾砂、粉质黏土)、残坡积层(粉质黏土、角砾)、石炭系灰岩(中风化、微风化)和石炭系砂岩(全风化、强风化、中风化、微风化)等。测线L4中基岩面有明显抬升，西侧平均深度约为13 m。东侧横向110 m~135 m处存在U型基岩面异常A5，指示着潜在的岩溶发育区。
Figure 2. High-Density electrical method inversion data
图2. 高密度电法反演数据
根据反演结果图的电阻率实例以及结合前期勘察资料，将电阻率大于210 Ω·m的区域解释为中–微风化石灰岩，小于210 Ω·m的区域解释为填土、饱和砾砂、粉质黏土等第四系覆盖层。通过COMSOL6.2 (试用版)对其进行模型的构造其中设置长度150 m，宽度50 m，深30 m (图3)，该模型划分为四个区域，左下角区域设置电阻率为300 Ω，右上角红色区域设置电阻率为400 Ω，中间部分设置电阻率为200 Ω，将红色异常体区域设置为600 Ω，沿着y = 0布置59个电极，电极间距为2.5 m，通过不同的装置来对该模型进行正演模拟。
Figure 3. Model schematic diagram
图3. 模型示意图
4. 网格剖分
在有限单元法中，计算域的网格剖分是正演问题求解的关键环节。由于正演计算需将连续微分方程离散化，整个计算域的物理量分布最终需通过有限个离散点的数值表征，因此其计算速度与结果准确性直接取决于网格剖分质量及单元数量。非结构化网格作为常用剖分方式，主要通过Delaunay三角剖分与前沿推进法生成，高质量的网格剖分可有效提升正演数值模拟的精度与效率，故在实际应用中，网格生成方法的选择与网格质量优化均需重点关注。
COMSOL Multiphysics 6.2 (试用版)具备强大的网格生成能力，支持六面体与非结构化四面体网格剖分，并允许用户灵活实施局部网格加密，这使得COMSOL能够高效处理各类复杂结构与地质模型的网格剖分需求，为电阻率法正演模拟提供了可靠工具。具体而言，用户可通过软件“网格”模块直接选择四面体或六面体单元进行剖分；针对电阻率法正演计算的特点，用户常需在模型边界设置非均匀网格分布以适配计算需求，而COMSOL的网格工具可对边界网格进行单独设定。为进一步平衡计算精度与效率，还可借助扫掠网格、网格细化及梯度分布等功能，对目标计算域的最浅部、异常体位置及测点区域进行网格加密，同时将远离目标的边界区域剖分为稀疏网格，此举可显著减少非目标区域的单元数量，降低计算负担，最终提升整体计算效能。
尤其在实际地球物理勘探场景中，常需面对起伏地形、滑坡地形等复杂地貌，此时电阻率法正演模拟需对模型进行大范围剖分，而选择适配的区域离散方法成为保障计算精度、控制计算量的关键。基于此，本文采用非结构化四面体网格剖分：该类型网格的灵活性可支持在模型不同区域采用差异化大小与形状的单元，能够更好地适配复杂地质模型的剖分需求。综上，在有限单元法正演模拟中，通过COMSOL的丰富功能实现网格剖分的合理设置，对满足电阻率法正演的精度要求与效率需求具有重要意义。
通过对模型采用非结构化四面体网格剖分模型剖分后包含199,929个域单元，15,317个边界单元619个边单元(图4)。
Figure 4. Model after grid meshing
图4. 网格剖分后的模型
5. 正演结果
5.1. 对称四极装置
该装置特点是rAM = rNB，其记录点位于测量电极MN的中点，沿着测线x方向布设电极并不断移动记录点位，最后得出所测数据。对称四极装置电阻率剖面图见图5。
Figure 5. Symmetric quadrupole resistivity profile
图5. 对称四极装置电阻率剖面
5.2. 温纳装置
温纳装置是较为特殊的对称四极装置，其基本原理为：rAM = rMN = rNB = a，通过同时将供电电极A、B和测量电极M、N向右逐点移动，可以测量得到测深为a/2的一条测线数据；随后逐步增大a的大小即增大电极间距，并从初始A点开始逐点向右移动四个电极，可以测量得到增大电极间距后另外一条测线数据；重复上述测量过程，可以得到倒梯形的断面。温纳装置电阻率剖面图见图6。
Figure 6. Wenner array schematic diagram
图6. 温纳装置电阻率剖面图
5.3. 偶极–偶极装置
这种装置的特点是供电电极AB和测量电极MN均为偶极。通常取偶极长度AB = MN = a，偶极装置间隔BM = na，其中n为正整数称为偶极间隔数。偶极装置的记录点取在BM的中点，这也是整个装置的中点。偶极–偶极装置电阻率剖面图见图7。
Figure 7. Dipole-Dipole resistivity profile
图7. 偶极–偶极装置电阻率剖面
6. 结论
在岩溶地区的地质探测研究中，通过构建数学模型开展数值模拟，并与实地真实测量结果进行对比分析，发现不同高密度电阻率法装置的模拟结果在岩溶区分辨率上存在显著差异。根据对称四极装置的电阻率拟断面图可以看出异常体分布在x轴85~105 m，异常的顶部在AB/2的30处，对称四极装置对异常反应的灵敏度适中，能够有效探测到一定深度范围内的地质异常，但在深部探测时，受限于电流的扩散和衰减，其分辨率会有所下降，常用于初步探测地下大致的地质结构。温纳装置的电阻率拟断面图可以看出异常体分布在x轴90~105 m，异常体的顶部在AB/2的35处，温纳装置的探测结果与实际情况较为接近，尤其在浅部地层的分辨上具有优势，能清晰反映浅部岩溶的分布情况，在浅层地质调查中应用广泛。偶极装置的电阻率拟断面图显示异常部分在x轴80~100 m，异常体顶部在AB/2的45处，偶极装置通常会使异常形态复杂化，不过它对横向边界的识别精度更高，在确定岩溶边界时能提供更详细的信息。
尽管不同装置对异常均有不同程度的响应，但在实际应用中，都无法仅依据视电阻率拟断面图直接对模拟模型进行解释。这是因为在数值模拟过程中，虽然采用了三维正演计算等方法来模拟电流在地下介质中的传播，但实际地质环境复杂，存在多种地电干扰因素，同时数据采集过程中也不可避免地会产生误差，这些因素都会导致视电阻率拟断面图存在一定的不确定性。
本次数值模拟仅选取了高密度电阻率法中较为常用的对称四极装置、温纳装置和偶极装置这三种进行研究，具有一定的局限性。为了更全面地了解高密度电阻率法在岩溶地区的应用效果，其他装置的适用性有待进一步探究。实时调整模型参数，提高解释的准确性，从而推动高密度电阻率法在地质勘探领域的发展与应用。
NOTES
*第一作者。
#通讯作者。