高精度重力在城市地下空间探测研究中的应用

徐燕君，陈 明，伍卓鹤，张轩国

(1.广州市政府投资项目研究评审中心，广东 广州 510000；2.广东省地质物探工程勘察院，广东 广州 510000)

0 引 言

城市环境特殊，各种人为干扰如交通震动、电磁波干扰及地面建筑物覆盖、城市地面硬化等对传统的地球物理探测带来了困难。高精度重力测量以密度差异为探测基础，无需在地面布设电缆、打孔布设电极或检波器，无损高效，绿色环保，适合在城市进行高密度、高精度数据采集，为城市地下空间探测和工程地质勘察提供了重要支撑。

1∶1 000或更大比例尺的高密度、高精度重力探测是一种高新探测方法技术，由于其涉及的高精度重力仪及其测量技术、高精度三维坐标获取及地形改正技术、后期数据整理及异常提取技术等难度大，技术标准不明确，且均源于中大比例尺重力测量规程(张燕等，2015),不适合城市高密度、高精度重力测量，可执行性较差，技术盲点较多，实际工程应用很少。

通过城市地下空间探测技术研究，形成一套可操作的高密度高精度重力探测技术方法流程、数据处理、成果显示和仪器应用技术，以期对今后深化研究有所帮助。

1 城市高密度高精度重力探测方法

1.1 仪器应用

近年来引进的新一代高精度重力仪(精度大于±0.001×10-5m/s2)(肖凡等，2011)不仅适用于传统领域，还适用于城市地下空间探测，可识别微伽级异常(±0.050×10-5m/s2以内)，用于寻找地面浅表特殊地质体如岩溶洞、地裂缝、地下地道、地下空洞、基岩面起伏等。新型井中重力仪可划分地层，确定地层密度，圈定井旁溶洞或低密度地质体等(曾华霖，1999)。

1.2 高精度重力测地技术

1.2.1 测地技术 作为高密度、高精度重力探测组成部分，测地精度要求非常高，传统的经纬仪和水准仪虽能获得高精度数据，但实施难度大，普及率差。随着卫星定位及省级COSR系统的应用，可利用RTK快捷获取高精度三维坐标数据，为高精度重力探测提供了基础。由于城市空间探测需研究0～200 m以浅的地层结构，特殊目标体不大，其引起的重力异常在±50微伽级以内，中大比例尺重力测量规范中的测量技术适应性差，要想获取准确数据，需对测地工作方法技术进行适当调整。

测地工作采用RTK测量即可确保精度而满足要求，但城市施工时，城市街道两侧的高楼会产生干扰。实验表明，当RTK测量位置距高楼30 m以内时，高程相对误差在±0.04 m范围内，超过高精度重力测量的要求。因此，距高楼30 m以内时采用光学水准测量可满足要求，可使用RTK快捷获得高精度三维坐标数据，或布置的重力测点距高楼30 m以外(于红娟等，2017)。

1.2.2 重力仪底面与地面距离 重力仪读数时用三脚架固定仪器并保持水平，仪器底面至地面在不同测量点高差不同，对高密度、高精度重力测量的误差影响较大，忽视该误差可能导致假异常。中大比例尺重力测量规范中没有涉及重力仪及其固定三脚架测量时与地面高差的问题，如常规重力测量的精度均>0.015×10-5m/s2，重力测地一般以地面或桩头作为重力测量的起算点，对重力观测精度的影响是可以接受的。但高密度、高精度重力测量时若存在65 mm的高程差，引起的重力场值误差最大可达±20.1×10-8m/s2(表1)，会影响后续成果解释。

表1 不同高程差引起的重力场值统计结果Table 1 Statistics of gravity field values caused by different height differences

为消除不同测点重力仪底面与地面的高程差，通过试验在重力仪测量前丈量地面至仪器底面的高差并取准到0.1 cm，在重力异常计算中消除这部分距离引起的场值，可确保每个测点以地面同一高程作起点计算场值，提高重力测量精度。消除三脚架测量误差的影响，对提高高密度、高精度重力测量精度具重要意义，以往忽视该误差导致高精度重力测量出现较多假异常，影响了高精度使用效果。

1.3 重力仪读数影响因素

1.3.1 操作员与重力仪的距离 操作员属活动质量体，读数时与重力仪距离太近会影响读数准确性，为了解其影响程度，选择一地基稳定的场地进行试验。结果表明：操作员与仪器距离1 m时对重力仪读数的影响较大，误差可达0.006×10-5m/s2；距离扩至2 m后，操作员对重力仪的影响不大。因此，在高精度重力测量时，操作员应距仪器1 m以上，并确保不要走动。

1.3.2 场地晃动对重力仪读数的影响 重力仪的读数精度用仪器读数标准差来衡量：标准差越大，代表重力仪测量场地晃动越大，读数越不可靠；标准差越小，代表重力仪测量场地越稳定，读数越可靠。

城市地区施工时由于人多、车多，场地晃动大，重力测量时较难获得高精度可靠数据。根据城市地区场地试验结果，读数标准差<0.2时，读数可靠；标准差>0.2时，场地存在晃动干扰因素，需保持场地稳定方可读取可靠数据。

1.3.3 多种干扰条件下的重力仪读数时长 为研究复杂环境中仪器读数情况，在某汽车站附近场地进行25、30、35、40、50、60 s不同采样时长的读数试验。结果表明，读数时间越短跳动越大，读数时间越长采样率越高、读数越稳、精度越高，综合考虑，采准数据最佳读数时长为60 s。

1.3.4 工作区基点闭合时间 由于重力存在零点位移现象，工作时间越长零点位移改正越多，会对测点读数产生影响。对分段闭合与总闭合进行分析发现，在1个工作日内简单缩短闭合时间并不能保证提高读数精度，而应增加工作日内闭合次数。综合分析零点位移特征，选择符合仪器性能特征的短时间段进行零点位移改正，减少仪器零点位移带来的误差，方可提高重力观测精度。单段闭合时间应不长于3倍设计精度在仪器动态曲线上的零点位移用时。

1.4 近区地改与中间层改正

1.4.1 近区地改 地改模型有人工目估地改、仪器测量地改和网格高程节点地改(孙文珂等，2017)，如在某场地利用目估与测距仪相结合进行了人工地改，通过实际GPS测地高程数据DEM(工作区外扩20 m)生成节点1 m×1 m网格进行网格节点地改，结果二者差别明显，人工结合仪器地改有9.7%的测点超设计误差，1 m×1 m、2 m×2 m网格高程节点地改未出现超差。

因此，当探测浅层地质目标且工作场地不大时，应测定大比例尺地形图，以工作区外扩20 m作为地改区，采用1 m×1 m或2 m×2 m的网格高程点用计算机进行近区地改；当探测更深地质构造时，外扩200 m测定大比例尺地形图，按2 m×2 m或5 m×5 m网格节点作为20～200 m地改区。

一般情况下，测量区不大的中远区地改值对各测点的影响值基本一样，可作为背景场值，对结果影响较小。地改半径具体要求按设计精度参考大比例尺重力勘查规范进行。

1.4.2 中间层改正密度取值方法 重力异常值计算需进行中间层改正，采用不同中间层密度的计算结果不同。一般采用研究区中间层密度平均值，即研究层的平均密度。实际工作中，一般在现场采集密度标本求取平均值：若存在水泥路面则按水泥混凝土的密度参与计算平均密度；若各物性层层厚差别大，宜用密度与层厚计算物性层加权平均密度。

另一种方法是通过实测重力剖面，利用不同密度值计算重力场分析求取。某剖面的地表为素填土及第四系残积层(图1)，素填土密度为1.75 t/m3，第四系残积层密度为1.84 t/m3，平均密度为1.80 t/m3，有2 m高差的陡坎，通过系列中间层取值计算对比发现，当密度为1.80 t/m3时所得重力异常值与高程无关，为最佳改正结果(刘代芹等，2012)。

图1 某重力剖面不同密度中间层改正计算重力异常值与高程对比图(1.60～2.67 t/m3为计算各剖面重力异常值中间层密度)Fig. 1 Comparison of gravity anomaly values and elevations calculated by correction of middle layer with different density values in a gravity profile(1.60～2.67 t/m3 is the middle layer density range in calculating gravity anomaly values)

1.5 重力异常值与起算高程

为考察不同起算高程对重力异常值的影响，对某剖面进行0、5、10、15、20、25、30、35.16 m 8个起算高程的重力异常值计算，发现重力异常值与重力起算高程无关(图2)，所以相同测点重力场值不同，其相对变化异常的大小不受起算面的影响。

图2 不同起算高程重力异常剖面图Fig.2 Profile of gravity anomalies at different starting elevations

1.6 数据处理

1.6.1 重力异常分离 城市空间探测需查找浅而小的特殊地质体以解决工程地质问题。例如，基岩面以下的岩溶溶洞引起的重力异常幅值通常在几十微伽级，为分离岩溶溶洞的局部异常与背景异常，最好参照通过已有钻孔控制剖面进行正演获得的模型，或采用先验经验了解场地岩溶溶洞规模再进行场分离。某剖面的重力异常呈左高右低、变化较大的特征(图3)，主要反映该处基岩面为斜坡面，采用多项式进行场分离效果明显。剩余重力异常提取需结合实际情况，采用多项式、趋势面或不同频谱等非线性方法进行剩余重力异常分离。

图3 某场地10线剖面多项式提取剩余重力异常图Fig.3 Map showing polynomial extraction of residual gravity anomaly from Line 10 in a site

面积性工作区剩余重力异常的提取不只考虑数据本身的特征，还要考虑背景场、异常体规模及深度等多因素。例如，在某灰岩地区(图4)探测隧道，通过趋势分析在隧道处提取的剩余重力负异常与隧洞十分吻合。图4b、c分别为三阶趋势分析的剩余重力异常图及区域重力异常图。

图4 趋势面法提取剩余重力异常示意图(a)某汽车客运站地下人行隧道重力异常图；(b)剩余重力异常图；(c)区域重力异常图Fig. 4 Schematic diagram of residual gravity anomaly extracted by trend surface method(a)Map showing the gravity anomaly of an underground pedestrian tunnel,near a bus terminal;(b)Map showing the residual gravity anomaly;(c)Map showing the regional gravity anomaly

1.6.2 层状密度界面反演 (1)线性回归法确定单密度界面。这是密度分界面近似计算的一种方法，属于定量解释的统计方法。

h=a+b·Δg

(1)

式(1)中：h为界面深度，Δg为界面起伏引起的重力异常，a、b为2个待求系数。基于该二元公式，工作区最少应有2个已知深度值，待求系数方有解，若已知深度较多，可用最小二乘法原理计算a、b值。

拟反演的界面为工作区主要密度界面，上下密度差较大，横向密度均匀性一般，上下两层有可能从实测数据中用数据处理的方法分离出所需重力异常。常用于高精度探测第四系与基岩分界面的计算。

对于寻找基岩起伏还可以采用单界面深度计算的Parker法。

(2)多界面多参数迭代拟合法。RGIS重磁电数据处理软件提供了三维多参数拟合的计算模块，在分离出目标重力异常后，根据合适的界面密度差、界面平均深度、迭代次数、滤波因子几个参数，求取最佳界面深度。基岩面重力异常可通过异常分离求取，当地表第四系或所求基岩面有钻孔控制，可采用正演计算基岩面以上地层的重力异常，然后用实测重力异常与其相减，所得目标重力异常会更可靠。

(3)人工交互反演界面。可通过对目标地质体进行重力异常提取后再进行二维、三维人工交互反演获得，特别是二维人工交互反演。国产RGIS及MAGS等软件均比较成熟，RGIS软件的三维人工交互反演可操作性也相当强。

不同的工作条件使用的界面反演方法不同，在现实条件许可时，具有详细物性条件及钻孔控制条件的人工交互反演界面为最优反演界面方法。

1.6.3 重力剖面反演法 为准确提取局部目标体异常信息，可利用重力异常的倾斜角、θ图边界识别、归一化标准差边界、物性反演法、归一化水平总梯度法、匹配滤波法、人机交互反演法等综合处理方法，确定局部目标体埋深及规模(伍卓鹤等，2015)。

根据剖面试验成果，物性反演法、归一化水平总梯度法能圈定目标体。

(1)自动物性剖面反演成像模型。利用3个不同埋深的低密度体模型正演出重力异常曲线，对应的3处低密度体重力异常呈明显重力低值(图5)。根据重力异常剖面曲线进行自动物性反演，断面(图5b)显示低密度区及重力低处与低密度体位置对应，说明可利用物性反演功能求取目标体，如灰岩中溶蚀溶洞、地裂缝、地下地道、地下空洞等。

图5 重力剖面自动物性反演剖面模型(a)重力异常曲线;(b)重力异常物性反演断面图Fig. 5 Profile model of automatic physical property inversion of the gravity profile(a)Map showing the gravity anomaly curve;(b)Inversion section of physical properties of gravity anomalies

(2)重力归一化水平总梯度法处理剖面成像模型。利用3个不同埋深的低密度体模型正演出重力异常曲线，对应的3处低密度体重力异常呈明显重力低值(图6)。利用该方法进行反演，其归一化水平总梯度显示局部高值对应于目标体顶界埋深及其投影范围(孙文珂等，2001)，可有效确定低密度体投影范围及其顶界埋深。若以3个不同模型体顶界面作为地层界面，也可用于不同密度界面的划分。

图6 归一化水平总梯度法圈定低密度体模型图(a)重力异常曲线;(b)重力异常水平总梯度模型断面图Fig. 6 Low density volume model delineation by normalized horizontal total gradient method(a)Map showing the gravity anomaly curve;(b)Sectional view of horizontal total gradient model of gravity anomalies

2 城市探测应用实例

2.1 某建筑场地高精度重力探测

在某建筑场地开展高精度重力测量工作，测网线距为10 m，点距为3 m，布格重力总精度达±0.014×10-5m/s2。

如图7所示，一阶趋势向右降低(图7a)不反映基岩面的起伏情况，只说明右侧基岩密度低于左侧，即右侧岩溶较为发育。通过一阶趋势分析分离出的剩余重力异常在-145、-110 m处呈现较大的低值异常，根据2.5D反演的岩溶溶洞空间位置(图7b)，并与20—20′线工程地质剖面(图7c)对比，发现钻孔ZK80与重力剖面在-145 m推断的溶洞十分吻合，只是形态略有差别；钻孔ZK89与推断的溶洞也有很好的对应性。重力测量发现在-115 m处存在较大重力低异常，根据2.5D反演结果认为，该异常为右侧深部基岩存在未被揭露的岩溶溶洞的反映。

图7 某工地20—20′线工程地质剖面综合解释图(a)重力异常曲线图；(b)重力2.5D反演断面图；(c)工程地质剖面图1-第四纪填土，密度为1.70 t/m3；2-第四纪残积层黏土，密度为1.90 t/m3；3-石炭纪灰岩及密度值，密度分别为2.70、2.67 t/m3；4-工程验证溶洞(充土)，密度为1.80 t/m3；5-重力推测溶洞(充土)，密度为1.80 t/m3；6-实测布格重力异常曲线/拟合重力异常曲线；7-区域重力异常曲线/剩余重力异常曲线；8-地形线；9-地层界线；10-钻孔位置、编号及标高Fig. 7 Comprehensive interpretation of Ling 20-20′engineering geological section of a construction site(a)Map showing the gravity anomaly curves;(b)Gravity 2.5D inversion section;(c)Engineering geological section

为确定某建筑场地基岩起伏状况，将每条测线的重力2.5D正反演基岩面结果与重力异常值进行最小二乘法拟合。重力反演基岩面与后期钻孔揭露的基岩面相近(表2)，重力反演与85个验证孔揭露的基岩面之差<1 m的占42.4%，在1～2 m之间的占21.2%，在2～3 m之间的占16.5%，在3～4 m之间的占10.6%，验证差>4 m的占9.4%，与钻孔揭露情况大体一致(表2)。因此，可以通过少数控制钻孔，利用重力资料进行分层或确定基岩面起伏情况，为后续工程提供依据。

表2 某建筑场地重力反演与验证孔基岩面对比Table 2 Comparison of gravity inversion and bedrock face of verified holes of a building site

2.2 自动物性及归一化总梯度模型反演综合应用

在某灰岩场地中寻找溶洞，该场地位于市中心，受交通和工业电磁场干扰较大，因此必需开展高精度重力测量工作。共施工L1—L6测线6条，每条长180 m，点距3 m。图8为L1测线成果，根据剩余重力异常曲线呈现的低值异常(图8a)，结合物性反演出的低密度区(图8b)和归一化总梯度法反映的偏正异常区(图8c)圈定溶洞特征的原则，推断出4处岩溶，同理在工作区推断岩溶溶洞29处。根据场地重力测量成果设计钻孔5个，经验证，4个钻孔与推断吻合，准确率达80%。

图8 广州某灰岩场地L1测线综合剖面图(a)重力异常曲线图；(b)重力剖面密度差断面图；(c)重力异常归一化总梯度模型断面图Fig. 8 Comprehensive profile of survey line No. L1 of a limestone site in Guangzhou(a)Map showing the gravity anomaly curve;(b)Density difference section of a gravity profile;(c)Sectional view of normalized total gradient model of gravity anomalies

3 结 论

为城市地区复杂环境中探测地质特征及浅层小目标体提供了一套高密度、高精度重力探测方法技术，取得下列成果。

(1)利用CORS系统的网络RTK技术可快速获取高精度三维坐标，首次提出消除重力仪底盘以下三脚架高度变化的影响，为城市地区高精度重力测地提供了一套技术方法。

(2)可从仪器读数时长、与操作员距离、基点闭合时间等角度提高了重力仪读数精度。

(3)在城市地区小面积工作区利用RTK实测高程数据DEM(工区外扩20 m)生成1 m×1 m网格节点地改，提高了近区地改精度。

(4)中间层改正密度采用工作场地内中间层密度平均值或通过实测重力剖面计算中间层改正值，起算高程对重力异常计算无影响。

(5)为突显城市地区浅层小目标体，可利用已有钻孔控制剖面进行正演获得模型进行场分离，或采用多项式、趋势面、不同频谱等非线性方法进行场值分离。

(6)单密度界面反演可利用线性回归法确定，多密度界面采用多参数迭代拟合法或人工交互反演法，人工交互反演法、自动物性反演及重力归一化水平总梯度法反演成像能有效提取浅而小的目标体。