0 引言

城市地下空洞是一类特殊的地下结构，在军用和民用领域都有广泛的用途，具有非常显著的实用价值。目前，地球物理探测技术被广泛用于地下空洞探测。地球物理探测技术指的是利用岩石之间的物理特性差异(比如密度、弹性等)进行探测的方法，典型的方法包括重力勘探、电法勘探、地震勘探、磁法勘探、放射性勘探、地热勘探等^[1-2]。其中，重力勘探具有能够相对快速地进行大区域物理场探测而不受外界电磁信号干扰和测区环境影响的抗干扰能力，对地下空洞、坑道及人工构造的探测具有明显优势。

现有的地下空洞探测方法多采用静态重力测量方法。Romaides等^[3]以加州范登堡空军基地的地下导弹发射设施作为研究对象，通过将实测重力异常与模型正演异常和区域布格重力异常相加进行比较的方式揭示了地下空洞的存在。王庆宾等^[4]通过假定地下空洞的位置和体积, 分析不同埋深及体积的人造空洞对地面重力及重力垂直梯度的影响，最后利用高精度相对重力仪CG-5进行了初步实验，表明重力垂直梯度对浅层地下空洞具有很好的探测能力。胡强等^[5]利用CG-5进行微重力观测实验，成功探测出同济大学校园内地下防空洞的位置。静态重力测量方法虽然能够探测出地下空洞，但应用于城市探测时，由于需要在地面上进行布线逐点式测量，并且为了保证测量精度需要在每个测点上停留足够长的时间来消除温差和震动等产生的噪声; 在未知地下异常体位置时，存在测量周期长、影响城市正常生活等诸多限制因素，从而导致探测难度大、效率低。

为了提高勘探效率，各国都十分重视移动平台重力测量的研究和开发，且随着海空领域快速、蓬勃的发展，移动平台重力测量在这些领域逐渐推广应用，展现出了较好的发展前景^[6-8]。在城市地下空洞探测时，由于重力测量范围较小、城市环境存在复杂的重力背景场源干扰，导致基于移动平台的重力探测方法在地下空洞探测鲜有应用。

因此本文针对静态重力测量方法的不足，提出基于车载重力测量平台的城市地下空洞快速探测方法。首先对城市环境中存在的楼房进行干扰校正，之后采用比值DEXP(depth from extreme points)方法^[9-10]和Tilt梯度^[11-13]方法对测得的重力异常进行深度成像和边界识别，探索车载移动平台重力测量技术在城市中开展的可能性，为实现快速探测城市地下空洞奠定基础。

1 楼房干扰校正

城市重力背景场干扰主要包括测点周围房屋建筑、车辆和人等引起的重力异常。测区周围的建筑物大多都是规则的长方体，对于连成一体的复合建筑物，我们把它看成规则长方体的组合。测点周围房屋对测点重力异常的影响可以等效为多个长方体对测点重力异常的叠加(图 1)，然后计算它们在测量点产生的重力效应，并从实测数据中去除。

图 1中实线部分为长方体。假设测点O的坐标为(0, 0, 0)，那么长方体在测点引起的重力异常^[14]为

式中：x₁、x₂、y₁、y₂、z₁、z₂表示长方体角点空间位置坐标；G为引力常数；σ为建筑物的平均密度。

在考虑建筑物的影响时，由于距离测点较远的房屋对测点的重力异常影响较小，通常只考虑距离测点250 m以内建筑物的影响^[15]。

若测区附近楼房建筑分布均匀，则可以认为每个测点受到250 m内楼房建筑产生的重力异常干扰是相同的。

进行测量时，车道为双向车道，包括车道r₁、r₂，车道中心线为x轴。测点O₁、O₂分别处于车道r₁和r₂上。对距离测点250 m正方形区域内的房屋进行建模，测点周围建筑物俯视图如图 2所示，其中a、b、s、s₁、s₂分别为楼房的长、宽、道路宽度、x方向楼房间隔、y方向楼房间隔。

通过式(1)可计算出每个长方体对测点的重力异常值，经过叠加可得由测点附近房屋引起的在测点处的重力异常值。从图 2中可知，由于建筑物的上下对称性，测点O₁、O₂在车道上受到的建筑物重力异常干扰是对称的，并且保持在一定范围。

2 基于比值DEXP成像的深度估计

成像反演技术能很好地对地质体的深度和分布信息进行估计。为了解决传统成像技术计算速度慢、趋肤效应明显的问题，提出了快速成像反演方法。主要的快速成像方法包括：极值点深度估计(DEXP)方法、小波成像方法、概率成像方法、广义线性反演方法等^[16-19]。

Fedi提出了DEXP方法进行地质体的深度估计及分布成像^[11]。DEXP成像反演技术基于简单的位场数据向上延拓得到，不需要任何滤波处理方法对测量数据进行噪声压制，一定程度上对城市地面上车辆、行人等的干扰进行了滤除，因此保证了成像结果的稳定性。但是该方法包含了与地质体形状相关的构造指数，当地下异常体类型未知时，构造指数的不准确将会给估计的深度结果带来很大误差。因此，本文采用基于传统DEXP方法进行改进的比值DEXP变换方法。

DEXP成像理论的关键概念尺度函数τ_P定义为重力异常数据f的对数或者f的P阶垂向导数f_P的对数相对于深度z的对数的导数，表达式为

在均衡场中可以简化为

式中：z₀为地质体的深度；N_P为f_P的构造指数。N_P定义为

式中, N为构造指数，与地质体几何形状相关。

比值DEXP变换利用重力异常不同阶垂向导数比值的DEXP方法进行重力数据的解释，在成像过程中独立于构造指数，使得方法完全自动化。公式为

式中：W_nR_{m, n}是关于R_{m, n}的DEXP变换；R_{m, n}为重力异常数据的不同阶垂向导数的比值；m和n为重力异常的不同阶数。

式中：0＜γ＜1；f₀通常取最大重力值。

在进行比值DEXP变换时，分别对重力数据不同阶垂向导数f_m和f_n进行向上延拓，根据式(6)对这两组延拓数据在不同高度上的比值来构建R_{m, n}，输出W_n(R_{m, n})的极大值即对应地质体的深度。当利用式(5)估计出地质体的深度后，就可以利用式(2)、(3)、(4)在z=-z₀的尺度函数值得到构造指数。

3 基于Tilt梯度的边界提取

梯度是重力场变化最剧烈的地方，重力水平梯度的极大值和垂直梯度的零值基本反映地质体边界，边界提取算法也是基于此提出的。

常见的边界识别算法有水平梯度模、Tilt梯度法、解析信号、Tilt梯度的水平导数等。本文采用简单典型的Tilt梯度法进行边界识别。

Tilt梯度是场的垂直梯度比上水平梯度模的arctan角度：

式中：g_x、g_y、g_z是场在x、y、z方向的一阶导数；T_dr为Tilt角度。

g_z可能为正值也可能为负值，在地下空洞处为负，在地下空洞外为正，在边界位置为零；g_x、g_y和的模总为正值，在物体边界位置处取最大值；T_dr在空洞处为负值，在边界位置为零，在物体外为正值。所以对于深部场源，即使它的垂直梯度和水平梯度都很小，两者的比值仍然会很大，边界的位置在零位置处取得，不受埋深深浅的影响。把Tilt看成是一个角度而不是一个比例的话，由于arctan值的属性，不管场的垂直方向一阶导数和水平方向一阶导数和的模的振幅为多大，T_dr所有的振幅都限制在(-π/2，π/2)。这个事实使T_dr的运行效果像自动增抑控制过滤器(automatic gain control, AGC)一样，均衡重力异常的输出振幅。

通过对重力异常数据进行Tilt梯度边界识别可以得到较为清晰的地下空洞边界信息。

4 实验及分析 4.1 测区概况及仪器

本次实验区域位于长春市南关区卫星广场附近，测区周围环境复杂，地下存在轻轨三号线轻轨站等多个地下空洞建筑，测线图如图 3所示。实验所用仪器是北京航天控制仪器研究所研制的捷联式重力仪SAG，该仪器的主要性能参数见表 1。

4.2 实验步骤

1) 测量系统搭建

移动平台重力测量系统由SAG重力仪和控制箱共同组成。为了减小车体行进过程中震动对测量产生的影响，测量前首先将重力仪和控制箱安装在减震板上，再用3M牌胶带将减震板固定在别克GL8型车中。

2) 仪器准备

由于重力仪在测量过程中温度和气压的变化将影响仪器内部重力敏感单元零漂误差，为了尽可能减小零漂误差，在实际测量之前对仪器进行了4 h充电预热。为了确定重力基准点，在出发地进行静态前校和静态后校各60 min。

3) 获取数据

为了减小测量中车辆和人的影响，实验选择在车流量、人流量少的清晨4:00—6:00进行，车辆尽量保持匀速行进以减小车载移动平台运动时引起的震动影响。测量路线选择沿道路左右两侧往返进行，最终回到测量起始点。

4.3 精度评价

本次实验通过对所有重复测点求观测误差的方法进行精度评价，精度评价公式为

式中：ε_g为观测误差；x_ij为第j个重复测点的第i点观测值；l_j为第j个重复测点的观测点数；μ_j为第j个重复测点的观测数据均值；n为重复测点数。

在城市的车载移动中，仪器在不停记录测量，由于车辆等待交通指示灯时静止不动，仪器在原地采样会提供大量的重复观测数据，此时测点的经纬度坐标同时保持不变，测点可以视为重复测点。

经过统计，本次实验的重复观测点共有20个，每个重复测点上平均产生19个测量值。根据式(8)的计算方法，求得本次实验数据测点观测误差为±0.115 1 mGal，满足精度要求。

4.4 数据处理及解释

选取具有地下轻轨站处道路两侧的测线A、B作为分析对象。首先进行楼房干扰校正。地面楼房的长、宽、高平均尺寸为55、15、25 m，长方向平均间隔为10 m，宽方向平均间隔15 m，车道宽度20 m，楼房剩余密度取2.6 g/cm³。将由式(1)得到的所有长方体在各个测点产生的重力异常分别进行叠加，得到车道上不同位置测点的重力异常范围为-0.24~-0.63 mGal。图 4为测线A和B经过楼房校正后重力异常值与原始重力异常值的对比，可以看出，楼房校正值保持在一定范围内变化，楼房干扰对原始数据总体趋势影响较小。

然后，利用垂向重力梯度异常与重力异常比值的DEXP变换对地下轻轨站附近两条车道测线上校正后的重力异常数据进行处理，构建比值R_{2, 1}。垂向重力梯度异常与重力异常分别向上延拓300 m，步长10 m，γ为0.2。将两条测线单独进行分析，成像结果见图 5。由图 5可见，测线A、B分别在800和1 000 m附近存在负异常，推测为地下空洞。经计算，其深度分别为24和28 m。查相关资料^[20]得知轻轨站深度约为27 m，与反演成像深度相符。

图 6为测线A、B垂向重力梯度估计的尺度函数。测线A、B在对应深度处的尺度函数值分别是-0.59和-0.53，确定垂向重力梯度的构造指数分别为0.18和0.06，接近于岩墙或者台阶的构造指数^[21-22]。

同时利用式(7)对校正的重力异常数据进行Tilt梯度边界识别，得到如图 7所示的重力异常图。图 7显示，在700~1 300 m区域附近存在负值异常，推测有地下空洞。经调研，地下轻轨站位置与图中测线下方的长方形负值异常范围接近。而图 7中右上角存在的负值重力异常处存在穿过卫星广场的地下隧道，由于仪器精度对不同大小和深度的地下空洞探测能力不同，需要做进一步的精确验证。

5 结论与建议

1) 本文利用车载重力测量对城市地下空洞进行快速探测，经过对实验数据进行反演成像和边界识别等处理，定位出地下轻轨站所在位置，对车载移动平台重力测量在城市地下空洞探测展开的可能性和有效性进行了验证。

2) 传统上对地下异常体的探测需要进行大量布线静态测量，但是对于城市重力探测不便于实现；而本文采用车载动态重力测量可以大大减少测量和布线时间，提高效率，在不知道具体地下空洞位置时可以起到预测性探测效果。

3) 研究中虽然对于地下轻轨站探测较为准确，但是在测量过程中并不能完全去除轻轨站附近部分地下隧道的影响，需要进一步实验研究出城市车载重力探测分离距离较近两个地下空洞的方法。

4) 在城市车载重力探测时，只能沿车道进行测量，为了提高测量准确性，应该尽量保持车辆慢速行驶，并且尽可能一条道路多次测量，以增加测量数据。