中国岩溶和采空区分布较为广泛，是世界上岩溶和采空区最发育的国家之一。近年来，随着我国经济快速发展，高速公路建设投入比重不断提高，岩溶和采空区日益增多，遭受的地质灾害也越来越多。在众多高速公路线路上，岩溶和采空区塌陷、地面沉降等问题突出，严重威胁着人们的生命安全^[1-2]。据不完全统计，有详细记录的岩溶和采空区塌陷共1 500多起，形成了超过4万个塌陷坑，对工程建设造成严重影响。为此，在高速公路建设前期及运营过程中进行岩溶和采空区探测具有重要的工程价值和现实意义。

目前，岩溶和采空区的探测方法主要为工程钻探和物探^[3]。钻探方法直观，但施工时间长，且往往只是“一孔之见”，由于岩溶和采空区发育具有复杂性和多变性，少量钻孔难以把地下岩溶和采空区完全揭露出来。物探方法可以为岩溶和采空区提供有用的地下信息，在岩溶和采空区的探测中得到广泛的应用。Legchenko等^[4]用核磁共振方法对死海沿岸地区岩溶充水洞穴的体积进行了估算。Redhaounia等^[5]使用高密度电法成功地探测到地下不同深度的石灰岩空洞。Duan等^[6]利用井间地震CT成像探测技术研究深部岩层中的洞穴，探明了岩溶和采空区的数量及规模，为各地区的注浆方案提供参考。Imposa等^[7]利用地震折射法研究了填满碎石材料的掩埋洞穴。Caselle等^[8]对地质雷达在石膏采石场岩溶构造探测中的应用进行了研究。随着工程建设的需要，对勘探技术的精度和准确度的要求不断提高，而多解性却一直是地球物理探测的固有难题^[9]，单一的探测方法对地质预测的准确性不够可靠，不同方法对不同地质缺陷的探测效果也不一样，目前还没有哪种预报方法能够准确预测各种地质缺陷。各种探测手段相互结合，取长补短，相互补充，相互验证，可提高探测的可靠性^[10-11]。为此，不少学者致力于采用综合地球物理探测的方法来研究高速公路路基下的岩溶和采空区。

当前，对岩溶和采空区的研究大多是多种地球物理探测方法的简单搭配和组合，不同方法之间约束联合反演的研究比较少。工程实例研究较多，对各种类型的岩溶或采空区的正反演成像规律研究较少。地质雷达(GPR)和高密度电阻率法(ERT)是电法勘探中最为成熟的2种方法^[12-13]，二者在刻画岩溶和采空区地质特征方面具有高效性和互补性^[14]，采用地质雷达和高密度电法相结合的综合探测方法，对常见的岩溶和采空区类型进行研究，通过数值模拟，得到了高密度电法和地质雷达对不同岩溶和采空区的响应特征，总结了不同岩溶和采空区的正反演成像规律，根据地质雷达对异常体边界识别和定位效果好的特点，将地质雷达探测结果与已知地质资料等先验信息作为形态约束和不等式约束条件施加到高密度电法反演方程中，通过对比不同类型岩溶和采空区不施加约束和施加约束的高密度电法反演成像结果，突出了施加约束反演方法的优越性；将上述结果应用于山东章丘模型试验探测中，成功探明了区域内的岩溶和采空区的类型及分布，为高速公路的治理方案提供了依据和参考。

1 岩溶和采空区综合探测方法 1.1 传统高密度电阻率反演

传统的电阻率反演方法为含光滑约束的最小二乘法，其目标函数为

$\varPhi=\varphi_d+\lambda \varphi_m=(\Delta \boldsymbol{d}-\boldsymbol{A} \Delta \boldsymbol{m})^{\mathrm{T}}(\Delta \boldsymbol{d}-\boldsymbol{A} \Delta \boldsymbol{m})+\lambda(\boldsymbol{C} \Delta \boldsymbol{m})^{\mathrm{T}}(\boldsymbol{C} \Delta \boldsymbol{m}),$

(1)

式中，Δd为实际观测数据与理论观测数据的差向量；Δm为模型参数的增量向量；A为偏导数矩阵；C为光滑度矩阵；φ_d为实际观测数据与理论观察数据的方差；φ_m为相邻网格电阻率的差异；λ为拉格朗日常数，表示φ_d与φ_m的权重。

传统的电阻率反演成像效果同时受到实际观测数据和光滑约束的控制，除施加光滑约束外，并未施加其他先验信息约束；光滑约束是一种比较宽松的约束，对反演多解性压制有限，探测效果不理想，需要携带更多先验信息的约束形式引导反演计算。

1.2 携带先验信息约束的高密度电阻率反演

综合地球物理方法是目前解决多解性问题的主要研究方向，地震反射法、地质雷达等勘探地球物理方法在识别和定位异常体界面方面效果较好，可以获得丰富的地质信息，对于改善电阻率的多解性具有重要意义。目前的综合探测大多是对多种探测结果进行对比分析，以“相互验证、相互补充”的思路进行综合评价，而对将其他勘探方法所获得的先验信息施加到反演过程中的研究比较少，未能充分发挥已知探测信息的约束作用。采用形态约束的方法，将地质雷达测得的异常体的位置、形态等先验信息施加到常规电阻率反演的目标函数中，构成携带形态约束的反演目标方程^[15]：

$\varPhi=\varphi_d+\lambda \varphi_m+\eta \varphi_s=(\Delta \boldsymbol{d}-\boldsymbol{A} \Delta \boldsymbol{m})^{\mathrm{T}}(\Delta \boldsymbol{d}-\boldsymbol{A} \Delta \boldsymbol{m})+\lambda(\boldsymbol{C} \Delta \boldsymbol{m})^{\mathrm{T}}(\boldsymbol{C} \Delta \boldsymbol{m})+\eta(\boldsymbol{F} \Delta \boldsymbol{m})^{\mathrm{T}}(\boldsymbol{F} \Delta \boldsymbol{m}),$

(2)

式中, η为形态约束的权重。

为了进一步降低反演问题的多解性，可以利用携带模型电阻率信息的不等式约束^[16]，在反演方程中加入已知信息

$\rho_{\min _i} \leqslant m_i \leqslant \rho_{\max _i}, \quad i=1, 2, \cdots, m,$

(3)

式中，m_i为第i个网格的电阻率；ρ_mini和ρ_maxi分别为第i网格的电阻率的下限和上限。需要指出的是，探测区域电阻率的变化范围可以通过钻探、地质分析和地质常识等手段获得。

从光滑约束、形态约束和不等式约束3个方面，得到电阻率反演的目标函数：

$\begin{gathered} \varPhi=(\Delta \boldsymbol{d}-\boldsymbol{A} \Delta \boldsymbol{m})^{\mathrm{T}}(\Delta \boldsymbol{d}-\boldsymbol{A} \Delta \boldsymbol{m})+\lambda(\boldsymbol{C} \Delta \boldsymbol{m})^{\mathrm{T}}(\boldsymbol{C} \Delta \boldsymbol{m})+\eta(\boldsymbol{F} \Delta \boldsymbol{m})^{\mathrm{T}}(\boldsymbol{F} \Delta \boldsymbol{m}), \rho_{\min _i} \leqslant m_i \leqslant \rho_{\max _i}, \\ i=1, 2, \cdots, m, \end{gathered}$

(4)

式中, η为形态约束的权重。

通过在电阻率反演方程中施加形态约束和不等式约束，提高了反演的可靠性，改善了异常区域的成像效果，抑制了多解性，降低了地质解释的难度。

高速公路路基下岩溶和采空区综合探测过程如图 1所示，利用地质雷达能准确地刻画出岩溶和采空区边界的特征，以地质雷达探测结果与已知地质资料等先验信息作为形态约束和不等式约束条件，施加到高密度电法反演方程中，通过综合分析地质和物探结果，探明了高速公路路基下岩溶和采空区的分布情况。

2 数值模拟

为了研究高密度电法和地质雷达对岩溶和采空区的响应特征，验证综合探测约束反演方法的有效性，采用合成观测数据进行数值模拟。利用高密度电阻率法，将60个电极按2 m的极距布置在被探测目标表面，以温纳装置的形式进行探测。地质雷达采用一体式天线，其天线的中心频率为100 MHz。

2.1 不同类型岩溶和采空区地质原型

采用5种岩溶和采空区地质模型(充填(水)、无充填(水)、倾斜充填(水)、倾斜无充填(水)、半充填)，研究高密度电法和地质雷达对于不同类型岩溶和采空区探测的一般规律，地质原型如图 2所示，其物性参数如表 1和表 2所示。

2.2 不同类型岩溶和采空区的高密度电法反演成像规律

利用传统的高密度电法反演方法，对合成观测数据进行反演，结果如图 3所示。

2.2.1 充填

反演剖面如图 3(a)所示，背景电阻率较高，在40~80 m范围内有1个椭圆形低电阻区(20 Ω·m)，由于高密度电阻率法受体积效应影响，低阻异常区域的范围比实际地质原型要大，不能准确地识别充填岩溶和采空区的发育位置、边界及规模。

2.2.2 无充填

反演剖面如图 3(b)所示，背景电阻率相对较低，在43~75 m范围内有1个椭圆形高电阻区(1 100 Ω·m)，高阻异常区域的范围要比实际地质原型大，且有2处明显的低阻假异常，难以准确识别无充填岩溶和采空区的发育位置、边界及规模。

2.2.3 倾斜充填

反演剖面如图 3(c)所示，背景电阻率相对较高，在40~80 m范围内存在1个不规则的低电阻区(35 Ω·m)，异常体的范围、形状与实际地质原型相差较大，无法判断低阻异常体的真实形态与走向，难以辨识倾斜充填岩溶和采空区的发育位置和形态规模。

2.2.4 倾斜无充填

反演剖面如图 3(d)所示，背景电阻率相对较低，在25~70 m范围内有2个高阻异常区(800~1 200 Ω·m)，与实际地质原型存在明显的差异，右侧高阻异常范围与实际地质原型有明显的差异，左侧高阻异常为假异常，不能准确判别实际地质原型中倾斜无充填岩溶和采空区的形态和发育位置。

2.2.5 半充填

反演剖面如图 3(e)所示，在20、60、100 m位置分别存在1个高阻异常区(950~1 150 Ω·m)，在60 m位置处有1个低阻异常区(30~80 Ω·m)，与实际地质原型相比，存在2处明显的高阻异常(20 m和100 m位置处)。60 m处的高电阻区和低电阻区的位置、形态与地质原型相差较大，发育范围比地质原型大，不能较好地定位和识别半充填岩溶和采空区。

2.3 不同类型岩溶和采空区的地质雷达正演成像规律

对本算例中的不同类型岩溶和采空区地质原型进行地质雷达正演模拟，得到地质雷达正演模拟结果，如图 4所示。

2.3.1 充填

从图 4(a)可以看出，电磁波向地下传播过程中遇到物性分界面，且物性界面向下延伸，表明在地下存在异常体，异常体的上部反射界面位置与实际地质原型异常体的位置一致，但不能确定异常体的下部边界，由此推断，充填岩溶和采空区的雷达信号受到强烈干扰，大幅衰减，观察不到充填岩溶和采空区的下部边界。

2.3.2 无充填

从图 4(b)可以看出，电磁波向地下传播过程中遇到物性分界面，且物性界面向下延伸，表明在地下存在异常体，是异常体的上层反射界面和下层反射界面与真实地质原型异常体的位置一致，可较好地反映无充填岩溶和采空区的形态、位置和边界。

2.3.3 倾斜充填

从图 4(c)可以看出，电磁波向地下传播过程中遇到物性分界面，表明在地下存在异常体，异常体的位置及界面均与地质原型较为一致，能较好地定位和识别异常体，但无法判断异常体的充填情况。

2.3.4 倾斜无充填

从图 4(d)可以看出，物性分界面的存在表明地下存在异常体，是异常体的位置和界面均与地质原型相吻合，容易对异常体进行准确地定位和识别，但异常体的充填情况难以识别。

2.3.5 半充填

如图 4(e)所示，电磁波向地下传播过程中遇到物性分界面，且物性界面向下延伸，反映了在地下存在异常体，异常体处有3层物性分界面，与地质原型半充填岩溶和采空区的上下边界比较一致，能准确识别异常体的位置与形态，异常区域分为两层，推断上层无充填对电磁波影响小，下层含水溶洞或采空区因为厚度薄、衰减程度小，在底部可观察到物性分界面。

2.4 不同类型岩溶和采空区的已知信息约束的高密度电法反演成像规律

通过地质雷达探测得到的岩溶和采空区的位置、形态等信息以及已知的地质资料，对于高密度电法而言是重要的先验信息，可以将其作为先验约束条件施加到高密度电法反演方程中去，反演结果如图 5所示。

2.4.1 充填

如图 5(a)所示，背景电阻率较高，在55~65 m范围内存在方形低电阻区(10~50 Ω·m)，位置和形态与实际地质原型相吻合，反映出充填岩溶和采空区的位置和规模。

2.4.2 无充填

如图 5(b)所示，背景电阻率较低，在56~66 m范围内存在方形高电阻区(950~1 200 Ω·m)，位置和形态与实际地质原型相吻合，较好地反映了无充填岩溶和采空区的发育特征。

2.4.3 倾斜充填

如图 5(c)所示，背景电阻率较高，在60~80 m范围内存在倾斜长条状低电阻区(10~30 Ω·m)，位置、走向和形态与实际地质原型较一致，较好地反映了倾斜充填岩溶和采空区的走向和位置；但图中有2处假异常，推测右侧数据不在倒梯形内部，敏感性较差，反演可能会不准确。

2.4.4 倾斜无充填

如图 5(d)所示，背景电阻率较低，在60~80 m范围内出现倾斜长条状高电阻区(900~1 200 Ω·m)，其位置、走向和形态与实际地质原型一致，较好地反映了倾斜无充填岩溶和采空区的发育特征。

2.4.5 半充填

如图 5(e)所示，在54~65 m范围内存在1个高电阻区(1 000~1 150 Ω·m)和1个低电阻区(10~20 Ω·m)，位置和形态与实际地质原型一致，较好地反映了半充填岩溶和采空区的发育特征。

总体而言，高密度电法能够探测到地下岩溶和采空区，但由于体积效应，岩溶和采空区的范围比实际地质原型更大，容易产生过多假异常，因此，利用地质雷达能准确刻画岩溶和采空区边界的优点，将其正演结果作为已知先验信息施加到高密度电法反演方程中，得到了较好的结果，提高了地质解释的准确性。

3 模型试验 3.1 试验场区概况

岩溶和采空区模型试验的地点位于山东省济南章丘，场地环境如图 6所示。选择一处地下介质相对均匀、无明显异常、长50 m×宽8 m的矩形区域作为岩溶和采空区构造埋设和探测试验区。选取试验模铸空腔结构模拟无充填岩溶和采空区，结构外壳厚度为5 cm，采用水泥、粉煤灰、碎石和水混合而成，经调整配比试验，确定含水结构的配比(重量)为水泥：粉煤灰：碎石：水=200：1 000：300：160。在选定试验区埋放3块不同尺寸试验模铸空腔结构，置于地下设计位置，试验模铸空腔结构尺寸以及埋深如表 3所示。随后进行粘土回填和地表标记等工作，岩溶和采空区分布模型示意图如图 7所示。设计粘土电阻率为250 Ω·m，而试验模铸空腔结构电阻率大于500 Ω·m。

3.2 地质雷达试验探测

在试验区沿测线进行地质雷达探测，采用100 MHz发射天线，地质雷达探测结果如图 8所示。

可以看出，试验区地下存在3个同相轴强度变高的区域，推测可能存在异常体，并与实际填埋位置进行对比，发现异常区域与实际填埋位置一致，可以大致推测异常体的位置，但不能准确判断是否含水。

3.3 高密度电法反演

在地质雷达布设测线上开展高密度电法探测，采用温纳装置形式，共有30根电极，电极间距为2 m，观测数据总数为132个。

传统高密度电法探测数据的反演成像结果如图 9(a)所示。将地质雷达探测的异常体的位置和形态信息作为已知先验信息施加到高密度电法反演约束中去，反演结果如图 9(b)所示。

通过对比分析，发现携带已知信息约束的高密度反演成像结果与实际情况相吻合。图 9(a)异常体的位置与实际相差较大，且多余结构和假异常较多，给岩溶和采空区的识别和解释带来了很大的干扰和困难。图 9(b)准确识别试验区3个高阻异常体的形态和位置，无多余结构和假异常，结合高密度电法对岩溶和采空区的响应特征及反演成像规律，推断出3个异常体为无充填岩溶和采空区。

综上所述，高密度电法与地质雷达相结合的综合探测方法在探测岩溶和采空区方面具有优越性，将地质雷达探测的异常体的形态和位置作为已知信息约束施加到高密度电法反演中，得到的电阻率结构中多余构造和假异常大大减少，背景更加干净，改善了反演多解性。该方法利用地质雷达定位精度高、界面识别效果好等优点，有效提高了异常体的定位精度和界面的分辨率。

4 结论

1) 对岩溶和采空区的不同类型进行了系统的数值模拟，分析了高密度电法和地质雷达探测浅部岩溶和采空区的成像规律，提出了利用地质雷达获得的先验信息的高密度电阻率约束反演方法。

2) 在山东章丘模型试验场地进行了探测，成功地探明了浅层岩溶和采空区的分布情况，表明了高密度电阻率约束成像法进行浅部岩溶和采空区的探测是行之有效的，为设计岩溶和采空区的治理方案提供了依据。

3) 文中提出的综合探测方法适合浅层岩溶和采空区的探测，为了更准确识别深部岩溶和采空区，建议采用钻孔雷达和钻孔电阻率相结合的综合物探方法，以获得更好的工程探测效果。