基于直立长方体模型的山体隆升与重力变化关系数值模拟

梁晨昊 申重阳 王嘉沛

摘要：基于形变与密度变化耦合运动理论，利用时变场内重力垂直梯度的计算方法，采用直立长方体模型，根据青藏高原平均降升速率，模拟计算在艾黎地壳均衡模式下，地表形变所引起重力及其垂直梯度的变化。结果显示，在山体抬升过程中，伴随着地表的隆升，重力值亦逐渐减小，导致其减小的原因为介质体密度减小与测点的高度增加。伴随着山体最高点抬升了5 cm，在最高点处重力变化为-14 μGal，对应的重力梯度约为-2.6 E。重力垂直梯度与静态场重力梯度存在一定的差异，其原因为在重力梯度场中加入了时间效应。

关键词：耦合运动;地表形变;重力变化;重力垂直梯度

中图分类号：P315.721   文献标识码：A   文章编号：1000-0666（2021）02-0162-08

0 引言

20世纪70年代以来，重力、水准、GNSS等观测技术不断进步，观测精度不断提高，目前，GNSS观测以及水准测量精度可达到毫米级别，重力测量精度亦可达到微伽。近年来，我国的大陆重力与GNSS观测网的时空分辨率有了大幅度提升（李强等，2012;申重阳，2005）。现阶段的GNSS和重力观测所提供的数据，可以很好地将重力与形变相结合，用来探究地壳运动及地震机理，取得了许多研究成果，特别是在青藏高原地区隆升与重力变化关系的研究中（邢乐林等，2017;段虎荣等，2020），形变与重力数据的结合发挥了很好的作用。但在研究过程中，依然存在两者数据结合不足的情形。

因此，将形变与重力数据更好地结合，研究两者之间的关系，探究重力变化机理及其与形变的关系成为现代地球物理学与大地测量学的研究热点。自20世纪以来，已有不少学者在理论上对其进行了初步的探究，Walsh（1975）首先推导出形变引起重力变化的计算公式，率先从理论上进行了分析。Reilly和Hunt（1976）指出了Walsh所得结果的错误，并给出了地球表面固定不变时形变所引起的重力变化，在该研究的基础上，陈运泰等（1980）完善了形变引起的重力变化理论，给出了区域形变和物质迁移引起的重力变化效应公式。李瑞浩（1988）采用不同原理，推导出和陈运泰等（1980）一样的结果。上述研究结果只适用于准静态变形情形，尚没有真正涉及形变时变模型。为此，申重阳和李辉（2005）提出了地壳变形与密度变化的耦合运动的思想与理论，给出了地壳变形与密度变化的耦合运动产生的重力场时变公式。王嘉沛等（2015）结合形变与密度变化耦合运动理论，利用直立长方体模型和川滇地区GNSS水平运动观测结果，模拟计算了该地区水平运动对空间固定点所产生的重力效应，但尚未考虑垂直运动效应和地表固定观测点伴随地表运动的情形。同时，对于山脉地区的形变研究，邢乐林等（2017）和段虎荣等（2020）分别利用重力变化数据研究了青藏高原的地壳增厚和隆升速率，但没有考虑到在伴随着青藏高原隆升过程中，地壳内部介质的密度是否发生变化。

因此，本文基于形变与密度变化耦合理论，运用直立长方体重力异常模型，模拟计算在山体隆升运动（以青藏高原为例）下产生的重力时变效应，计算在时变场内重力垂直梯度的变化，并对形变、密度与重力之间的时变关系作进一步研究。

1 基础理论

1.1 介质体形变所引起的重力变化的一般表达式

地壳形变会引起介质体质量的重新分布，从而引起重力变化。一些学者（Walsh，1975;陈运泰等，1980;申重阳等，2005，2007）从理论上分析了形变与重力变化的关系，为后面的模拟计算提供了理论依据，即以地心为原点（0，0，0）（定点），设置相对地球的惯性直角坐标系（x，y，z）。地球内部时刻发生运动，对于任意时刻t，设地球物质的集合为Ω（t）， 地球表面用S（t）表示， S（t）同时包含了地球表面和地球内部洞穴的表面。假设，在S（t）上存在任意一点P0， 其坐标设置为（rs，t），其中rs=xsi+ysj+zsk。在地球内部存在一个任意点Q（r，t）∈Ω（t），其中r=xi+yj+zk，点的密度减和位移可表示为ρ（r，t）≠0和u（r，t），在Q点的形变速度可用位移u（r，t）对时间的一阶导数来计算u·（r，t）。将Q到P之间的向量关系表示为：R=r0-r，其数值大小表示为：R=r0-r。

整个地球运动过程均满足质量守恒定律，我们将以地球内部任意点Q为中心的介质体元设置为dv（r），其密度变化为：

1.3 模拟山体形变模型

山体在隆升过程中，对其周围地区的地貌会产生很大的影响，同时也会发生形变作用，从而引起相应的重力变化。根据形变引起的重力变化的一般表达式以及质量守恒定律可知，介质体在发生形变作用的同时会伴随着密度的变化。研究山体隆升运动所产生的重力变化时应结合形变与密度变化耦合理论，应遵循质量守恒定律。

地球板块在相互运动过程中，由于板块之间的碰撞（例如印度洋板块和亚欧板块）导致板块边界挤压形成山体。假设在运动过程中满足艾黎地壳均衡模型，即可把地壳视为较轻的均质岩石柱体漂浮在较重的均质岩浆之上，处于静力平衡狀态。根据阿基米德浮力原理，山越高，增加的质量越多，陷入岩浆越深，形成山根。根据艾黎模型的均衡理论，可以知道在山体隆升过程中，山体的高度越高，其补偿深度也越深，莫霍面的深度也就越深。假设山体的海拔高度为H，地壳的原始厚度为T，山根的大小Z与山体海拔H之间的关系为：

Z=4.45H（14）

此时，山体的地壳厚度可以表示为：T+Z+H。

如图3所示，在t0时刻，地表和莫霍面处于水平状态，并未发生形变。此后，在受到板块运动等某些因素的影响时，其发生了变形，形成了山体，根据艾黎地壳均衡模型的均衡理论，在山体隆升的同时，莫霍面深度也在加深，就形成了t5时刻的状态。随后板块不断运动，山体不断隆升，莫霍面不断加深，形成了t10时刻的状态。因此艾黎地壳均衡模型是后续模拟计算山体隆升时所产生的重力变化的基础。

2 模拟计算方法

本文采用地壳形变与密度变化耦合运动理论来计算地壳形变引起的重力变化。该方法结合了地壳形变与密度两种信息，综合性地研究了地球重力场变化特征，改善了利用单一资料研究重力变化的方法，更加有利于对地球动力学进一步的研究。对于计算地壳形变引起的地

表观测点的重力变化公式，可利用直立长方体（图4）进行近似计算。

利用直立长方体计算产生的重力异常为：

3 模拟计算结果

本文通过数值计算分析青藏高原以5 mm/a（邢乐林等，2017;段虎荣等，2020）的速率隆升时所产生的重力变化，时间步长取10 a，初始时刻为t0，t5时刻表示以t0时刻为基础山体隆升5 cm之后的状态，t10时刻表示以t0时刻为基础山体隆升10 cm之后的状态，分别测量t5和t10时刻重力以及重力梯度分布。在计算时由于密度会随时间和位置变化，设置初始密度为ρ=2.67 g/cm3，在隆升过程中所产生的密度变化可利用式（2）计算求得。假定地表隆升形变范围为（50×50）km2。观测网大小定位（70×70）km2，这样可以将形变区域完全覆盖，在计算时将地表划分为（0.5×0.5）km2的網格。在计算过程中假设每个块体密度都是均匀的，其密度变化也是均匀的。

本文计算得到了形变过程中地表介质重力分布。在初始状态（t0时刻），地表未发生形变，密度均匀，其地表重力分布如图6所示。t5时刻，由于板块运动等作用，地表发生隆升形变。t5时刻对应的重力及重力变化情况如图7a所示。从图中可以看出，在山体隆升5 cm所产生的重力差值为-14 μgal。t10时刻，地表隆升形变还在持续，隆升速度与前一个时间段相同，其对应的重力及重力变化值如7b所示。同时，利用式（13）分别计算了在地表隆升过程中（t0～t5，t5～t10时段）山体的地表重力垂直梯度，如图8所示。从图8中可以看出山体的大致形状。在山体地区，重力垂直梯度约为-2.6 E，在山体周围的地区重力垂直梯度较大，约为-5 E。在t0～t5和t5～t10时段山体隆升运动过程中，山体隆升的量级相同，所以在运动过程中所产生的重力变化和重力梯度相同，同时重力变化和重力垂直梯度的空间分布情况同模拟的山体空间分布情况一致，因此导致两次运动过程的重力差值和梯度图像相同。

4 结论

本文基于形变与密度变化耦合运动理论，给出了在时变场里重力垂直梯度计算公式的表达形式，并利用直立长方体模型模拟山体在抬升过程中所产的重力变化以及重力垂直梯度的变化情况，主要得到以下结论：

（1）在山体抬升运动过程中，由于形变与密度变化耦合，形变即发生密度变化。因此，利用质量守恒定律所求得的山体密度在海拔高的地区较小。在山体抬升过程中，密度值逐渐减小，观测点与介质的距离逐渐增加，导致重力值逐渐减小。第一个时间段内（t0～t5），地表抬升5 cm，其重力梯度约为-2.6 E。在第二个时间段内（t5～t10），地表隆升5 cm，其重力梯度约为-2.6 E。在整体运动过程中地表隆升10 cm，重力垂直梯度约为-2.6 E，重力梯度分布与山体形态分布相同，重力垂直梯度与静态场重力梯度存在一定的差异，其原因为在重力梯度场里考虑了时间效应。

（2）结合形变与密度耦合定律，模拟计算了在青藏高原隆升时，在密度变化的情况下，所产生的重力变化和重力垂直梯度的变化情况。计算结果显示，青藏高原在两次隆升过程中，重力梯度没有发生变化，重力变化同时受到地表隆升和地下介质变动的影响。在研究青藏高原隆升时所产生的重力变化时，应当考虑密度变化和形变效应的双重影响。

参考文献：

陈运泰，顾浩鼎，卢造勋.1980.1975年海城地震与1976年唐山地震前后的重力变化[J].地震学报，2（1）：21-30.

段虎荣，康明哲，吴绍宇，等.2020.利用 GRACE时变重力场反演青藏高原的隆升速率[J].地球物理学报，63（12）：4345-4360.

李强，游新兆，杨少敏，等.2012.中国大陆构造变形高精度大密度GPS监测——现今速度场[J].中国科学：地球科学，42（5）：629-632.

李瑞浩.1988.重力学引论[M].北京：地震出版社.

申重阳.2005.地壳形变与密度变化耦合运动探析[J].大地测量与地球动力学，25（3）：11-16.

申重阳，李辉.2007.研究现今地壳运动和强震机理的一种方法[J].地球物理学进展，22（1）：49-56.

王嘉沛.2015.地壳变形与密度变化耦合运动引起的重力变化效应研究[D].北京：中国地震局地震研究所.

王嘉沛，申重阳，玄松柏.2015.全球地壳模型CRUST1.0在青藏高原东南部的重力检核[J].大地测量与地球动力学，35（4）：621-626.

邢乐林，王林海，胡敏章，等.2017.时变重力测量确定青藏高原地壳隆升与增厚速率[J].武汉大学学报（信息科学版），42（5）：569-574.

Reilly W I，Hunt T M.1976.Comment on An analysis of local changes in gravity due to deformation by J B Walsh[J].Pure & Applied Geophysics，114（6）：1131-1133.

Walsh J B.1975.An analysis of local changes in gravity due to deformation[J].Pure & Applied Geophysics，113（1）：97-106.

Simulation of the Relationship between the Surface Deformationand the Gravity Variation Based on Vertical Cuboid

LIANG Chenhao，SHEN Chongyang，WANG Jiapei

（Key Laboratory of Earthquake Geodesy，Institute of Seismology，China Earthquake Administration，Wuhan 430071，China）

Abstract

Based on the coupling theory of deformation and density variation，the vertical gradient of gravity in time-varying field is calculated.By the help of the vertical cuboid and average uplift rate of the Tibetan plateau，the process of the mountain movement is simulated，and the variation of gravity and its vertical gradient caused by surface deformation is simulated in accordance with the Airy equilibrium model.The results show that in the process of mountain uplift，the gravity value gradually decreases.It is speculated that the gravity decrease is caused by the decrease of medium density and the increase of the height of the measuring point.The mountaintop where the maximum gravity value appears rises for 5 cm，as a result，the maximum gravity changes for-14 μGal，and the corresponding gravity gradient is about -2.6 E.Due to the time effect in the gravity gradient field，there is a certain difference between the vertical gravity gradient and the static gravity gradient.

Keywords：coupled motion;surface deformation;gravity variation;vertical gravity gradient

收稿日期：2020-12-25.

基金項目：国家自然科学基金项目（41674018）资助.

第一作者简介：梁晨昊（1996-），硕士，主要从事形变重力解释方面研究.E-mail：lch158768@163.com

通讯作者简介：申重阳（1963-），研究员，主要从事形变重力解释、地震前兆机理、地球物理场及其反演、工程地震等方面的研究.E-mail：scy907@163.com.