地球重力场是地球的基本物理场之一,其时空变化反映了地球表面和内部物质分布变化、地球的形变及质量迁移、各圈层之间物质交换和耦合等(孙和平等,2017).目前大范围长时间尺度的全球重力场是通过卫星进行观测的,卫星观测到的地球重力场空间分辨率可达400km,时间间隔约30天,观测精度为毫伽(mGal,10-5m·s-2)量级(Tapleyetal.,2004;Visser,1999).由于卫星重力的空间分辨率是有限的,所以不能反映区域重力场的精细空间分布.地表流动重力测量、航空或车载重力测量是目前陆地上测定区域重力场精细空间分布的主要观测手段(柴华等,2015;高晟俊等,2018;胡平华等,2017;于瑞航等,2015;祝意青等,2017).由于卫星和航空或车载重力测量的观测精度在毫伽量级,使得许多较小的重力变化不能在观测中显现,但流动重力测量的观测精度可达到十几微伽(μGal,10-8m·s-2)的量级甚至更高(FuandShe,2017;毛经伦和祝意青,2018),能够观测到较小的重力场变化,因此除非在观测比较困难的地区(例如山区、沙漠等)采用航空重力观测外,目前国内外高空间分辨率和高精度的区域地表重力场观测所采用的主要观测技术手段仍为地表流动重力测量.近十几年来,随着各类地表重力测量的快速发展,形成了一种新的地表流动重力测量手段,称为混合(微)重力观测(HybridGravimetry或HybridMicro-gravimetry)(Portieretal.,2018).混合重力观测采用绝对重力仪控制长期重力变化,流动重力测量提高空间分辨率,连续重力观测提高时间分辨率,GPS/GNSS测定测点位置和高程,已成功用于地震、地热、地下结构反演等观测研究中.混合(微)重力测量的主要特点是观测区域相对较小,测点间距离较短,多种观测技术同期进行.本研究也是混合(微)重力测量实验研究,主要目的是研究城市内的流动重力观测的观测能力,验证其是否能够为城市地表下沉提供重力约束.

引起城市内的地面沉降的原因可分为两类:人类活动和自然因素.人类活动主要是指城市施工、地下水开采等;而自然因素包括降雨、地下水非人为变化、特殊地质构造(例如:喀斯特地形)等因素(杨天亮和许言,2017;袁铭等,2016).传统监测地面沉降的方法主要是水准测量(崔效文,2018;吴红让,2018)和GPS/GNSS监测(罗建忠等,2017),目前由于InSAR技术的快速发展,城市地面沉降能利用InSAR技术精确实施观测(王平豪,2017),全国的大约几十个城市(北京、上海等)都利用该技术进行了地面沉降观测(刘凯斯等,2018;秦晓琼等,2016),武汉市也有观测结果(Baietal.,2016).这三种地面形变观测技术都有各自的优势和缺点,因此目前各项观测技术在地面沉降监测中都有应用.GPS/GNSS监测属于点观测,优点是可实现实时监测,广泛应用在地面形变的动态监测中(何杰,2017),缺点是无法实现整个区域的面观测,当然多个测点可以插值成面观测,但是毕竟不是面观测.水准测量是线观测,虽然无法像GPS/GNSS一样实现实时监测,但是也可获得较短时间内的地面形变,在工程形变和长期形变监测中应用较多(崔效文,2018),缺点是无法实现面观测和实时监测.利用InSAR技术进行地面形变观测属于面观测,由于卫星重复周期的原因,所以无法实现实时监测,因为能够对整个观测区域进行面观测,目前在地面沉降观测中得到广泛应用,缺点是无法实时监测.上面提到的三种地面形变的观测技术,观测量都是几何形变量.但是在讨论地面形变产生机制方面,只有地面几何形变观测量,不能充分地解释形变产生的原因.重力观测是物理观测量,可以观测到由于地下物质迁移或地下结构形变引起的重力变化,因此在地面沉降观测中,重力观测技术是非常有必要的,是对上面三种形变观测技术的有益补充.但是目前利用重力观测技术进行地面形变观测的工作非常少(章传银等,2018),主要有以下三个方面的原因:(1)重力仪价格昂贵,维护、运输和观测的费用非常高,这是最主要的原因;(2)重力观测属于线观测,无法实现面观测,也无法实现实时监控;(3)重力虽然是物理观测量,但是也是一种综合效应,需要结合其他观测资料来共同解释观测结果.尽管如此,重力能反映地下质量迁移或结构变化引起的地面形变,是上面三种形变观测技术所不能的,因此可为解释地面形变提供几何观测量之外的物理观测约束.

InSAR的观测结果表明,在2009年至2010年约一年的时间里,由于城市施工等原因,在武汉城区有最大67.3mm的沉降(Baietal.,2016).如此大的地表沉降在城市中可能会引起较大的破坏,造成巨大的经济损失和环境灾害,因此地表沉降监测是城市不可缺少的.观测和研究城市中的地表沉降,可为防范灾害和城市建设提供重要参考.鉴于InSAR观测到的地表沉降为几何观测量,在解释引起沉降的机制方面无法对来自地下的因素提供约束,因此本研究试验利用重力观测资料来弥补InSAR在此方面的不足,为更科学合理的解释地表沉降机制提供依据.本研究基于武汉市内InSAR观测结果(Baietal.,2016),在武汉市内地表沉降较大的区域布设了流动重力观测网,进行了7期流动重力观测实验.观测方案中减小了测点间的距离来减小噪声可能对观测结果的影响,测线中绝大多数测点间的距离都小于3km.实验中采用绝对重力观测对整个观测网进行控制,采用考虑了CORS基站的GPS精确测定了测点的经纬度和高程,并采用GLDAS水模型进行陆地水储量变化的重力影响改正,最后用D-InSAR的垂直位移观测结果进行了比较验证,重力变化的解释中还采用了长江水位观测数据.

1、观测方案和方法

本实验在武汉市InSAR观测地表沉降比较大的区域布设了流动重力观测网,如图1所示.考虑到城市内测点离交通道路比较近,仪器受来往交通工具引起地面振动的影响较大,另外人为噪声(例如居民区空调的外机、建筑工地的施工噪声等)对仪器的影响也比较大,测点的选择尽量避开了这些因素的影响,并且测点间的距离尽量均匀并尽量小于3km.由图1中可知,布设的观测网着重考虑了地表沉降较大的区域和长江沿岸,在长江沿岸设点主要是为了能够观测长江水的季节性变化导致的重力变化,因为已有的研究表明重力观测受局部水质量迁移的影响非常大(Güntneretal.,2017;LampitelliandFrancis,2010).另外为了减少重力长期变化对观测结果的影响,在测网中的一个固定测点(点位编号:wh1700,图1中用★表示)进行了绝对重力测量.最终观测方案测点个数为64个(与图1中给出的观测方案仅少数几个点不一致),并且每个测点的位置和高程都利用GPS进行了精确测定.

参照“中国大陆构造环境监测网络相对重力联测技术规程(2010年修订)”中的观测方法,本实验的观测方法为:(1)采用3台CG-5重力仪同时进行观测;(2)采用串式对称观测,即A→B→C→…→C→B→A;(3)一条测线要求在1天内闭合.主要技术要求:(1)一条测线结束后,及时进行验算,验算时应加入格值、潮汐改正、气压改正、仪器高改正和零漂改正,重力段差的三台仪器互差限差优于40×10-8m·s-2.(2)平差后联测精度优于20×10-8m·s-2,点值精度平均值优于15×10-8m·s-2.

2、数据处理

整个实验观测网络的数据处理主要包括以下三个方面:(1)绝对重力控制;(2)陆地水储量变化重力影响改正;(3)流动重力观测数据处理.每一期的观测值考虑点位变动等因素对异常观测点进行了剔除,剔除依据为上文给出的“主要技术要求”.整个观测的时间段是从2016年3月10日到2017年5月21日,每期观测大约为10天.

2.1 绝对重力控制

为了对整个观测网的重力变化进行控制,实验中采用对同一测点进行绝对重力测量获得观测期间整个观测网络的重力变化.绝对重力控制测点在图1中用红色★表示,在编号为IGG014-1基墩上,绝对重力观测值的仪器高为1.0m,观测仪器为FG5-112绝对重力仪,绝对重力观测值(已去掉一个常数)见表1.选定的wh1700测点就在隔壁的观测室内,距IGG014-1基墩的直线距离约15m,实验中假设这两测点处的绝对重力变化是一样的,因此用IGG014-1基墩的绝对重力变化就可以对整个流动重力观测网进行控制.从表1中可以看出,在观测期间内,IGG014-1基墩的绝对重力变化非常小,最大为234.9×10-8m·s-2,最小为230.0×10-8m·s-2,两者相差仅4.9×10-8m·s-2,约为整网平差后点值平均精度的一半,相对来讲重力变化比较小.

图1InSAR观测的地表沉降(Baietal.,2016)和流动重力观测网

因为绝对重力的变化较小,数据处理中按时间采用线性插值得到每期观测值的绝对重力改正值,每期取插值时间为该期观测中间的日期,插值得到的每期改正值见表2.研究中也对绝对重力改正前后的数据做了比较,发现测点的最大差别约为2.6×10-8m·s-2,对最终各期间结果的影响较小,但是此项改正因为会影响整个观测网的重力值,所以需要考虑此项改正来控制整个重力观测网.

2.2 陆地水重力影响改正

本实验中采用GLDAS全球水模型计算测区整体的陆地水重力影响,模型时间间隔为3小时,空间分辨率为1°×1°,而流动重力测量单条测线1天闭合,因此计算中把每一天的重力影响的平均值对相应的流动重力测量的同一天的观测数据进行改正.根据流动重力测量的观测时间段,计算了wh1700测点处陆地水重力影响,时间段为2016年3月1日到2017年5月30日,计算的重力变化见图2.从图2可以看出,计算的陆地水重力影响约在(24～29)×10-8m·s-2之间,因此影响最大约5×10-8m·s-2,在一年内有明显的季节性特征.最终把计算的陆地水重力影响在每天的观测数据中扣除,然后再整网平差.关于陆地水重力影响改正需要说明的是,因为GLDAS全球水模型空间分辨率为1°×1°,并不能精确反映测区内水质量迁移的空间变化,因此是区域平均的结果.测区内有长江、汉江和一些湖泊,其季节性的水储量变化对重力观测的影响比较大,本研究在后面有关于长江沿岸测点的重力变化的分析讨论.这里加入GLDAS全球水模型计算的陆地水储量变化重力影响,是为了去掉测区整体大范围的水质量变化的重力影响,改善测区整体的重力变化,提高整网平差精度.

图2由GLDAS水模型数据计算的wh1700测点处陆地水重力影响

2.3 流动重力观测数据处理

流动重力观测数据预处理的主要过程包括:格值系数转换、固体潮改正、气压改正、仪器高改正和零漂改正.上述预处理过程的改正方法详见“中国大陆构造环境监测网络重力观测技术规程”附录3.经过上述预处理各项改正得到各个测站的相对重力预处理值,便可进行重力网平差处理,主要过程为:(1)根据上文由GLDAS水模型计算得到的陆地水储量变化的重力影响按日对原始观测进行陆地水引起的重力效应改正;(2)以测区中点位稳定的绝对重力控制点(wh1700)作为统一起算的拟稳基准,对每一期的观测进行拟稳平差计算,以获取统一起算基准下的重力变化;(3)先对七期的重力观测资料平差计算结果进行分析,初步了解各台仪器的观测精度后,确定合理的仪器先验方差后再重新平差计算,以得到最合理的解算结果;(4)按表2给出的插值后各期的绝对重力改正值,对各期平差后的结果进行绝对重力改正;(5)以第1期的观测结果为背景重力场,用第2到7期的观测结果减去第1期的观测结果,查看观测时间段内的重力变化;(6)选取重力变化明显的测点,查看其重力随时间变化的特征.其中陆地水储量变化重力影响改正是按日进行的,因此是在流动重力各条测线(每条测线当日完成全部观测)的观测值上进行改正的,因此需要各期观测每条测线的观测日期,本实验中流动重力观测各期的观测时间段见表3,由于天气和其他原因,各期中并不是每天都连续观测的,具体各期每条测线的观测日期要查证实际观测日期,最终各期的平差精度和平差后点值精度的平均值见表3.

3、结果及其分析

对于城市中的地表沉降,沉降速率是关键的指标,也能从该指标中预测未来的沉降量及其危害程度.由于本次实验观测不知道观测区域的背景重力场,又因为关键指标沉降速率只与相对重力变化量和时间间隔有关,因此本实验中将第1期实验的观测值作为背景重力场,计算获得其后6次观测相对于第1期的差值作为变化量,得到观测区域的重力时空变化,结果如图3所示.图3a—f分别是第2期到第7期观测减去第1期后重力变化量的等值线图,从表2中可知各期的时间间隔大约为2个月,因此图3a—f可认为是在第1期观测后每2个月的累积重力变化,图中的日期采用表2中约为各期观测中间的日期.

从图3a可以看出,在约2个月的时间内,整个观测区域(除汉江与长江的交接口处)的重力变化非常小,重力变化量都约在±10×10-8m·s-2之内,因为实验中流动重力观测的点值精度约为10×10-8m·s-2,因此该结果说明测区内2个月内的地表重力变化用流动重力观测是无法确认能被观测到的.至于汉江与长江的交接口处的重力变化,后面的分析表明可排除地表下沉引起的重力变化,因此该重力变化可能与此处长江水和(或者)与地下水的水储量变化有关,重力增大最大超过了30×10-8m·s-2,说明此处长江水和(或)地下水储量在这2个月内是增大的,但是增大的量很难精确模拟,因为需要地形、流速、降雨量、蒸发量、渗透水量等多种观测数据和复杂的建模.从图3b可以看出,在约4个月的时间内,除了长江和汉江交汇处的重力变化外,其他陆地上的重力变化也明显的显示出来,这些重力增加明显的地方应该和地表下沉有关,因此这些区域应该在下沉观测中重点关注.从图3c—f,重力变化明显的地方越来越多的显现出来,刻画出更多与沉降有关的区域,并且某些局部区域呈现重力逐渐增加的特点,说明这些区域在这段时间内是持续下沉的,虽然观测结果受到陆地水(降雨、蒸发、地下水、植被水、江河湖等水)变化的影响较大,但是整个区域重力逐渐增加的趋势是存在的,说明测区内整体地表是持续下沉的.并且图3d、e和f整个测区重力变化的区域基本一致,只是后面的图3f重力变化的空间分布更加细化,说明在观测时间段内,利用流动重力观测手段,可以在城市内观测到稳定重力空间变化特征,并且观测时间越长,空间特征越精细.从图3a和f可知,在观测时间段内重力的最大变化可达40×10-8m·s-2,并且局部区域的重力是逐渐增加的,说明用地表混合重力观测技术能够高分辨率观测到区域重力变化.

图3第2到7期重力结果相对于第1期的重力变化

图4D-InSAR观测的垂直位移

为了进一步说明重力观测结果能对城市地表沉降提供约束,实验采用D-InSAR观测的约同一时间段同一测区的地表下沉来验证.D-InSAR数据来源于Sentinel-1A对地卫星观测数据(欧阳伦曦等,2017),数据类型为Sentinel-1SLCIW,观测时间段为2016年3月12日到2017年5月18日,这个时间段和图3f的时间段大致相同.D-InSAR观测结果为视线向位移,根据研究区卫星中心入射角约为43.9218°,将视线向位移转化为垂直位移,由D-InSAR获得的垂直位移如图4所示.从图4可看出,地表下沉的区域和图3f重力观测中重力增加的区域大部分相一致,不过因为D-InSAR是面观测,空间细节观测更为细致,下沉量最大超过6cm的地方很多,主要在城市施工密集区域;还有一些地方特别是湖边和江边有地表上升,应该是由于填湖导致的.

为了便于和重力的结果进行比较,把图4给出的D-InSAR垂向位移作等值线图,并和图3f的空间范围一样,然后放一起比较,见图5.图5a中有6个重力增加的区域,其中中间4个都对应于图5b中地表下沉的区域,说明这4个区域的重力增加主要是由地表下沉引起的;图5a中另外2个重力增加的区域(位于图中右上角和左下角),在图5b中没有找到对应的地表下沉区域,说明这两个区域的重力增加并不是主要由地表下沉引起的.既然排除了地表下沉为主要因素,本研究后面的结果表明左下角的重力变化与长江水的季节性变化具有很好的相关性,说明该重力变化可能是长江水变化导致的,但是该结论并不确定,也有可能是地下水或其他因素造成的.至于右上角应该是与地下岩溶物质有关,因为此处地下有较大的岩溶物质存在(范士凯,2006).对比结果说明流动重力观测能圈定由于地表下沉引起的区域,因此用重力观测技术可以观测城市地表下沉,并且能为城市地表下沉提供重力约束.

另外也将重力等值线图直接覆盖在图4给出的D-InSAR垂直位移上进行了比较,如图6所示.从图6可以看出,测区中下沉的区域(底图中红色表示的区域)与重力增加的区域大部分是一致的.但在右上角和左下角重力增加的两个区域中,D-InSAR并没有观测到明显的下沉,具体原因上面已经解释过,此处不再赘述.而有些下沉比较明显的较小区域,重力也没有明显的增加,主要原因还是重力观测的空间分辨率(约3km)比D-InSAR的空间分辨率(2.3m)要大太多,较小区域的变化是无法在重力观测中反映出来的.

图5第7期相对于第1期重力变化(a)和D-InSAR观测的垂直位移等值线(b)的比较

图6第7期相对于第1期重力变化覆盖

图7长江汉口站水位变化(2016-01-01—2017-12-31)

为了验证重力观测受长江水变化的影响,研究利用长江水位汉口站的观测数据与流动重力观测的重力场变化进行了对比分析.长江水位观测站汉口站的位置见图6中红色▲所示,观测到汉口站水位变化如图7所示.提取的重力变化的点,一个是离水位观测站最为接近的wh1800测点(见图6),另一个是长江岸边离水位观测站较远的测点wh0700(见图6),并且要求这两个点在重力变化较为明显的区域,观测值对比结果分别如图8a和b所示.因为重力观测是各期平差的结果,因此需要对水位按各期的观测时间段取平均值,另外重力变化是相对于第1期的变化量,因此水位也需要减去第1期的平均值,水位变化的结果见图8a和b中红色曲线所示.

从图8a可以看出,测点wh1800的重力变化与汉口站水位变化趋势一致,说明测点wh1800的重力变化与长江水位变化紧密相关,由于图4中D-InSAR观测到此处并没有明显的地表下沉,因此可认为在图3f中观测区域左下角的重力变化最有可能的原因是长江水的储量长期变化引起的,当然也可能是地下水储量长期变化引起的,具体的原因需要实际水文观测资料来验证.从图8b可以看出,测点wh0700的重力变化与汉口站水位变化并不同,汉口水位站和测点wh0700的直线距离约为6km,这说明无论是长江水还是地下水,其变化是有一定范围的,而且引起的重力变化也是有一定范围的,对距离相对较远的点可认为没有影响,这与已有的研究结果相符().这个范围可以用图3f中观测区域左下角的重力变化区域给出一个粗略估计,尽管这个重力变化可能不仅仅是由于长江水变化引起的,可以把这个区域变化为10×10-8m·s-2的等值线圈定为该区域水质量变化对重力影响范围的边界线.由InSAR的结果可知,地表在这个区域没有明显的下沉,反而整体略微有些上升,进一步说明多种观测手段的重要性,互相补充来解释沉降机制.这一点在图3f中右上角重力变化明显的区域更有说服力,因此在该区域InSAR观测到的沉降并不明显,但是重力却有明显的增加,说明地下物质质量在增加,最有可能是该区域地下水在这段时间内是增加的.但是在将来,这些增加的地下物质有可能再次流走,这个时候就要特别注意可能有地质灾害产生,因此这个区域是在城市建设中特别注意的地方,尽量不要有大量的地表建筑,例如居民区或商圈等.

图8长江汉口站水位变化(红色)与测点wh1800(a)和wh0700(b)重力变化(黑色)的比较

4、结论

本研究在武汉市内地表下沉较为明显的部分城区进行了流动重力测量,采用了绝对重力观测控制整个观测区域的长期重力变化,采用GLDAS陆地水模型数据计算并去掉了测区整体的陆地水重力影响,观测结果和D-InSAR观测的垂直地表形变进行了比较,并且把长江沿岸测点的重力变化与长江水位变化也进行了比较,讨论了引起地表下沉和重力变化的原因,可为城市建设提供科学依据.

重力观测的结果表明,每期整网平差的点值平均精度均小于10×10-8m·s-2,优于“重力联测技术规程”的精度要求(15×10-8m·s-2),证明采用重力观测技术能在城市内获得高精度的区域重力变化.重力观测结果与D-InSAR观测的垂直地表形变的比较表明,重力持续增加的区域与D-InSAR观测的下沉的区域大部分是一致的,说明利用流动重力能够观测到城市地表下沉.另一方面重力增加的区域不仅仅是由下沉引起的,可能是地下物质或长江水的变化引起的,因此利用多种观测技术进行城市沉降观测是有必要的.重力观测结果虽然没有D-InSAR的空间分辨率高,但是圈定的区域是一致的,并且重力还可以测到地下物质变化,能为解释地表沉降机制提供观测依据.重力观测结果与长江汉口站水位变化的比较表明,虽然不能肯定该区域长期重力增加是由于长江水在该处的储量增加引起的,但是长江水位的变化与重力观测有明显的相关性,并且重力增加有一定的范围.因此测区内其他重力变化明显的区域由于离长江距离较远,其观测结果主要还是地表形变引起的,这是重力观测结果与D-InSAR观测的地表下沉相一致的主要原因.

参考文献：

[1]柴华,王虎彪,武凛等.2015.陆地惯性重力测量静态试验研究.大地测量与地球动力学,2015,35(6):953-956,968.

[2]崔效文.2018.地面沉降监测中水准测量的应用.资源信息与工程,33(4):125-126.

[3]范士凯.2006.武汉(湖北)地区岩溶地面塌陷.资源环境与工程,20(S1):608-616.

[4]高晟俊,郝卫峰,李斐等.2018.极地航空重力测量及其应用进展.极地研究,30(1):97-113.

[5]何杰.2017.GNSS定位技术在地质灾害监测中的应用——以彭州市吴家岩不稳定斜坡为例[硕士论文].成都:西南交通大学.

[6]胡平华,赵明,黄鹤等.2017.航空/海洋重力测量仪器发展综述.导航定位与授时,4(4):10-19.

[7]刘凯斯,宫辉力,陈蓓蓓.2018.基于InSAR数据的北京地铁6号线地面沉降监测分析.地球信息科学学报,20(1):128-137.

[8]罗建忠,刘金柱,孔友谊等.2017.GPS在天津市地面沉降监测中的应用.海洋测绘,37(3):60-62.

[9]毛经伦,祝意青.2018.地面重力观测数据在地震预测中的应用研究与进展.地球科学进展,33(3):236-247.

[10]欧阳伦曦,李新情,惠凤鸣等.2017.哨兵卫星Sentinel-1A数据特性及应用潜力分析.极地研究,29(2):286-295.

[11]秦晓琼,廖明生,杨梦诗等.2016.应用高分辨率PS-InSAR技术监测上海动迁房歪斜形变.测绘通报,(6):18-21.

[12]孙和平,徐建桥,崔小明.2017.重力场的地球动力学与内部结构应用研究进展.测绘学报,46(10):1290-1299.

[13]王平豪.2017.PS-InSAR在沉降监测中的应用研究——以深圳福田区为例[硕士论文].北京:中国地质大学(北京).

[14]吴红让.2018.精密水准在城市区域地表沉降测量浅析.内江科技,(6):73,85.

[15]杨天亮,许言.2017.国际地面沉降与城市安全研究动态——第一届国际城市地质学术研讨会综述.上海国土资源,38(2):1-3.

[16]于瑞航,蔡劭琨,吴美平等.2015.基于SINS/GNSS的捷联式车载重力测量研究.物探与化探,39(S1):67-71.

[17]袁铭,白俊武,秦永宽.2016.国内外地面沉降研究综述.苏州科技学院学报(自然科学版),33(1):1-5.

[18]章传银,王伟,甘卫军等.2018.利用CORS站网监测三峡地区环境负荷引起的地壳形变与重力场时空变化.武汉大学学报(信息科学版),43(9):1287-1294.

[19]祝意青,梁伟锋,赵云峰等.2017.2017年四川九寨沟MS7.0地震前区域重力场变化.地球物理学报,60(10):4124-4131.

陈晓东,李航,邓明莉,江利明,孙和平,李德伟,张为民,郝洪涛,桑鹏,白林.用重力测量技术观测城市地表下沉的实验研究[J].地球物理学报,2020,63(08):2882-2892.

基金:国家自然科学基金(41974023,41674083,41874094,41604070);中国地震局地震研究所所长基金(IS201726162)联合资助.