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)分别利用重力变化数据研究了青藏高原的地壳增厚和隆升速率,但没有考虑到在伴随着青藏高原隆升过程中,地壳内部介质的密度是否发生变化。

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

本文采用地壳形变与密度变化耦合运动理论来计算地壳形变引起的重力变化。该方法结合了地壳形变与密度两种信息,综合性地研究了地球重力场变化特征,改善了利用单一资料研究重力变化的方法,更加有利于对地球动力学进一步的研究。对于计算地壳形变引起的地表观测点的重力变化公式,可利用直立长方体(图4)进行近似计算。

图4 直立长方体模型
Fig.4 Vertical cuboid model

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

式中:(x,y,z)为长方体外任意点的观测坐标;(ε,η,τ)是直立长方体内任意一点的坐标,r为长方体内任意点到长方体外任意点的距离; r=[(ε-x)²+(η-y)²+(τ-z)²]^1/2,ρ为直立长方体密度。利用直立长方体运动来计算重力变化的过程中,首先假设其运动的参数,然后根据质量守恒定律,计算运动过程中所产生的密度变化,以及在形变前后重力变化:

Δg=g₁(x₁,y₁,z₁,t₁)-g₀(x₀,y₀,z₀,t₀) (16)

本文利用直立长方体将山体进行剖分,用n个长方体近似模拟山体。如图5所示,将研究区域分成m×n个小块体单元块体,假设P[KG0](x,y,z)点在t时刻第i行、第j列的小单元块体产生的重力变化为g_ij(x,y,z,t),则P[KG0](x,y,z)在t时刻受到研究区域的大块体的重力变化为:

图5 地下分块模型
Fig.5 Underground block model

本文通过数值计算分析青藏高原以5 mm/a(邢乐林等,2017; 段虎荣等,2020)的速率隆升时所产生的重力变化,时间步长取10 a,初始时刻为t₀,t₅时刻表示以t₀时刻为基础山体隆升5 cm之后的状态,t₁₀时刻表示以t₀时刻为基础山体隆升10 cm之后的状态,分别测量t₅和t₁₀时刻重力以及重力梯度分布。在计算时由于密度会随时间和位置变化,设置初始密度为ρ=2.67 g/cm³,在隆升过程中所产生的密度变化可利用式(2)计算求得。

假定地表隆升形变范围为(50×50)km²。观测网大小定位(70×70)km²,这样可以将形变区域完全覆盖,在计算时将地表划分为(0.5×0.5)km²的网格。在计算过程中假设每个块体密度都是均匀的,其密度变化也是均匀的。

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

图6 t0时刻地表重力等值线图
Fig.6 Gravity contours at t0

图7 t5(a)、t10(b)时刻重力等值线与重力差值分布图
Fig.7 Distribution of gravity and gravity difference at t5(a),t10(b)

图8 t0～t5(a)、t5～t10(b)地表重力垂直梯度等值线图
Fig.8 Vertical gradient cotour map of surface gravity at t0～t5(a)and t5～t10(b)

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

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

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