梁晨昊 申重陽 王嘉沛

摘要：基于形變與密度變化耦合運(yùn)動理論，利用時變場內(nèi)重力垂直梯度的計算方法，采用直立長方體模型，根據(jù)青藏高原平均降升速率，模擬計算在艾黎地殼均衡模式下，地表形變所引起重力及其垂直梯度的變化。結(jié)果顯示，在山體抬升過程中，伴隨著地表的隆升，重力值亦逐漸減小，導(dǎo)致其減小的原因為介質(zhì)體密度減小與測點的高度增加。伴隨著山體最高點抬升了5 cm，在最高點處重力變化為-14 μGal，對應(yīng)的重力梯度約為-2.6 E。重力垂直梯度與靜態(tài)場重力梯度存在一定的差異，其原因為在重力梯度場中加入了時間效應(yīng)。

關(guān)鍵詞：耦合運(yùn)動;地表形變;重力變化;重力垂直梯度

中圖分類號：P315.721?? 文獻(xiàn)標(biāo)識碼：A?? 文章編號：1000-0666（2021）02-0162-08

0 引言

20世紀(jì)70年代以來，重力、水準(zhǔn)、GNSS等觀測技術(shù)不斷進(jìn)步，觀測精度不斷提高，目前，GNSS觀測以及水準(zhǔn)測量精度可達(dá)到毫米級別，重力測量精度亦可達(dá)到微伽。近年來，我國的大陸重力與GNSS觀測網(wǎng)的時空分辨率有了大幅度提升（李強(qiáng)等，2012;申重陽，2005）?，F(xiàn)階段的GNSS和重力觀測所提供的數(shù)據(jù)，可以很好地將重力與形變相結(jié)合，用來探究地殼運(yùn)動及地震機(jī)理，取得了許多研究成果，特別是在青藏高原地區(qū)隆升與重力變化關(guān)系的研究中（邢樂林等，2017;段虎榮等，2020），形變與重力數(shù)據(jù)的結(jié)合發(fā)揮了很好的作用。但在研究過程中，依然存在兩者數(shù)據(jù)結(jié)合不足的情形。

因此，將形變與重力數(shù)據(jù)更好地結(jié)合，研究兩者之間的關(guān)系，探究重力變化機(jī)理及其與形變的關(guān)系成為現(xiàn)代地球物理學(xué)與大地測量學(xué)的研究熱點。自20世紀(jì)以來，已有不少學(xué)者在理論上對其進(jìn)行了初步的探究，Walsh（1975）首先推導(dǎo)出形變引起重力變化的計算公式，率先從理論上進(jìn)行了分析。Reilly和Hunt（1976）指出了Walsh所得結(jié)果的錯誤，并給出了地球表面固定不變時形變所引起的重力變化，在該研究的基礎(chǔ)上，陳運(yùn)泰等（1980）完善了形變引起的重力變化理論，給出了區(qū)域形變和物質(zhì)遷移引起的重力變化效應(yīng)公式。李瑞浩（1988）采用不同原理，推導(dǎo)出和陳運(yùn)泰等（1980）一樣的結(jié)果。上述研究結(jié)果只適用于準(zhǔn)靜態(tài)變形情形，尚沒有真正涉及形變時變模型。為此，申重陽和李輝（2005）提出了地殼變形與密度變化的耦合運(yùn)動的思想與理論，給出了地殼變形與密度變化的耦合運(yùn)動產(chǎn)生的重力場時變公式。王嘉沛等（2015）結(jié)合形變與密度變化耦合運(yùn)動理論，利用直立長方體模型和川滇地區(qū)GNSS水平運(yùn)動觀測結(jié)果，模擬計算了該地區(qū)水平運(yùn)動對空間固定點所產(chǎn)生的重力效應(yīng)，但尚未考慮垂直運(yùn)動效應(yīng)和地表固定觀測點伴隨地表運(yùn)動的情形。同時，對于山脈地區(qū)的形變研究，邢樂林等（2017）和段虎榮等（2020）分別利用重力變化數(shù)據(jù)研究了青藏高原的地殼增厚和隆升速率，但沒有考慮到在伴隨著青藏高原隆升過程中，地殼內(nèi)部介質(zhì)的密度是否發(fā)生變化。

因此，本文基于形變與密度變化耦合理論，運(yùn)用直立長方體重力異常模型，模擬計算在山體隆升運(yùn)動（以青藏高原為例）下產(chǎn)生的重力時變效應(yīng)，計算在時變場內(nèi)重力垂直梯度的變化，并對形變、密度與重力之間的時變關(guān)系作進(jìn)一步研究。

1 基礎(chǔ)理論

1.1 介質(zhì)體形變所引起的重力變化的一般表達(dá)式

地殼形變會引起介質(zhì)體質(zhì)量的重新分布，從而引起重力變化。一些學(xué)者（Walsh，1975;陳運(yùn)泰等，1980;申重陽等，2005，2007）從理論上分析了形變與重力變化的關(guān)系，為后面的模擬計算提供了理論依據(jù)，即以地心為原點（0，0，0）（定點），設(shè)置相對地球的慣性直角坐標(biāo)系（x，y，z）。地球內(nèi)部時刻發(fā)生運(yùn)動，對于任意時刻t，設(shè)地球物質(zhì)的集合為Ω（t）， 地球表面用S（t）表示， S（t）同時包含了地球表面和地球內(nèi)部洞穴的表面。假設(shè)，在S（t）上存在任意一點P0， 其坐標(biāo)設(shè)置為（rs，t），其中rs=xsi+ysj+zsk。在地球內(nèi)部存在一個任意點Q（r，t）∈Ω（t），其中r=xi+yj+zk，點的密度減和位移可表示為ρ（r，t）≠0和u（r，t），在Q點的形變速度可用位移u（r，t）對時間的一階導(dǎo)數(shù)來計算u·（r，t）。將Q到P之間的向量關(guān)系表示為：R=r0-r，其數(shù)值大小表示為：R=r0-r。

整個地球運(yùn)動過程均滿足質(zhì)量守恒定律，我們將以地球內(nèi)部任意點Q為中心的介質(zhì)體元設(shè)置為dv（r），其密度變化為：

1.3 模擬山體形變模型

山體在隆升過程中，對其周圍地區(qū)的地貌會產(chǎn)生很大的影響，同時也會發(fā)生形變作用，從而引起相應(yīng)的重力變化。根據(jù)形變引起的重力變化的一般表達(dá)式以及質(zhì)量守恒定律可知，介質(zhì)體在發(fā)生形變作用的同時會伴隨著密度的變化。研究山體隆升運(yùn)動所產(chǎn)生的重力變化時應(yīng)結(jié)合形變與密度變化耦合理論，應(yīng)遵循質(zhì)量守恒定律。

地球板塊在相互運(yùn)動過程中，由于板塊之間的碰撞（例如印度洋板塊和亞歐板塊）導(dǎo)致板塊邊界擠壓形成山體。假設(shè)在運(yùn)動過程中滿足艾黎地殼均衡模型，即可把地殼視為較輕的均質(zhì)巖石柱體漂浮在較重的均質(zhì)巖漿之上，處于靜力平衡狀態(tài)。根據(jù)阿基米德浮力原理，山越高，增加的質(zhì)量越多，陷入巖漿越深，形成山根。根據(jù)艾黎模型的均衡理論，可以知道在山體隆升過程中，山體的高度越高，其補(bǔ)償深度也越深，莫霍面的深度也就越深。假設(shè)山體的海拔高度為H，地殼的原始厚度為T，山根的大小Z與山體海拔H之間的關(guān)系為：

Z=4.45H（14）

此時，山體的地殼厚度可以表示為：T+Z+H。

如圖3所示，在t0時刻，地表和莫霍面處于水平狀態(tài)，并未發(fā)生形變。此后，在受到板塊運(yùn)動等某些因素的影響時，其發(fā)生了變形，形成了山體，根據(jù)艾黎地殼均衡模型的均衡理論，在山體隆升的同時，莫霍面深度也在加深，就形成了t5時刻的狀態(tài)。隨后板塊不斷運(yùn)動，山體不斷隆升，莫霍面不斷加深，形成了t10時刻的狀態(tài)。因此艾黎地殼均衡模型是后續(xù)模擬計算山體隆升時所產(chǎn)生的重力變化的基礎(chǔ)。

2 模擬計算方法

本文采用地殼形變與密度變化耦合運(yùn)動理論來計算地殼形變引起的重力變化。該方法結(jié)合了地殼形變與密度兩種信息，綜合性地研究了地球重力場變化特征，改善了利用單一資料研究重力變化的方法，更加有利于對地球動力學(xué)進(jìn)一步的研究。對于計算地殼形變引起的地

表觀測點的重力變化公式，可利用直立長方體（圖4）進(jìn)行近似計算。

利用直立長方體計算產(chǎn)生的重力異常為：

3 模擬計算結(jié)果

本文通過數(shù)值計算分析青藏高原以5 mm/a（邢樂林等，2017;段虎榮等，2020）的速率隆升時所產(chǎn)生的重力變化，時間步長取10 a，初始時刻為t0，t5時刻表示以t0時刻為基礎(chǔ)山體隆升5 cm之后的狀態(tài)，t10時刻表示以t0時刻為基礎(chǔ)山體隆升10 cm之后的狀態(tài)，分別測量t5和t10時刻重力以及重力梯度分布。在計算時由于密度會隨時間和位置變化，設(shè)置初始密度為ρ=2.67 g/cm3，在隆升過程中所產(chǎn)生的密度變化可利用式（2）計算求得。假定地表隆升形變范圍為（50×50）km2。觀測網(wǎng)大小定位（70×70）km2，這樣可以將形變區(qū)域完全覆蓋，在計算時將地表劃分為（0.5×0.5）km2的網(wǎng)格。在計算過程中假設(shè)每個塊體密度都是均勻的，其密度變化也是均勻的。

本文計算得到了形變過程中地表介質(zhì)重力分布。在初始狀態(tài)（t0時刻），地表未發(fā)生形變，密度均勻，其地表重力分布如圖6所示。t5時刻，由于板塊運(yùn)動等作用，地表發(fā)生隆升形變。t5時刻對應(yīng)的重力及重力變化情況如圖7a所示。從圖中可以看出，在山體隆升5 cm所產(chǎn)生的重力差值為-14 μgal。t10時刻，地表隆升形變還在持續(xù)，隆升速度與前一個時間段相同，其對應(yīng)的重力及重力變化值如7b所示。同時，利用式（13）分別計算了在地表隆升過程中（t0～t5，t5～t10時段）山體的地表重力垂直梯度，如圖8所示。從圖8中可以看出山體的大致形狀。在山體地區(qū)，重力垂直梯度約為-2.6 E，在山體周圍的地區(qū)重力垂直梯度較大，約為-5 E。在t0～t5和t5～t10時段山體隆升運(yùn)動過程中，山體隆升的量級相同，所以在運(yùn)動過程中所產(chǎn)生的重力變化和重力梯度相同，同時重力變化和重力垂直梯度的空間分布情況同模擬的山體空間分布情況一致，因此導(dǎo)致兩次運(yùn)動過程的重力差值和梯度圖像相同。

4 結(jié)論

本文基于形變與密度變化耦合運(yùn)動理論，給出了在時變場里重力垂直梯度計算公式的表達(dá)形式，并利用直立長方體模型模擬山體在抬升過程中所產(chǎn)的重力變化以及重力垂直梯度的變化情況，主要得到以下結(jié)論：

（1）在山體抬升運(yùn)動過程中，由于形變與密度變化耦合，形變即發(fā)生密度變化。因此，利用質(zhì)量守恒定律所求得的山體密度在海拔高的地區(qū)較小。在山體抬升過程中，密度值逐漸減小，觀測點與介質(zhì)的距離逐漸增加，導(dǎo)致重力值逐漸減小。第一個時間段內(nèi)（t0～t5），地表抬升5 cm，其重力梯度約為-2.6 E。在第二個時間段內(nèi)（t5～t10），地表隆升5 cm，其重力梯度約為-2.6 E。在整體運(yùn)動過程中地表隆升10 cm，重力垂直梯度約為-2.6 E，重力梯度分布與山體形態(tài)分布相同，重力垂直梯度與靜態(tài)場重力梯度存在一定的差異，其原因為在重力梯度場里考慮了時間效應(yīng)。

（2）結(jié)合形變與密度耦合定律，模擬計算了在青藏高原隆升時，在密度變化的情況下，所產(chǎn)生的重力變化和重力垂直梯度的變化情況。計算結(jié)果顯示，青藏高原在兩次隆升過程中，重力梯度沒有發(fā)生變化，重力變化同時受到地表隆升和地下介質(zhì)變動的影響。在研究青藏高原隆升時所產(chǎn)生的重力變化時，應(yīng)當(dāng)考慮密度變化和形變效應(yīng)的雙重影響。

參考文獻(xiàn)：

陳運(yùn)泰，顧浩鼎，盧造勛.1980.1975年海城地震與1976年唐山地震前后的重力變化[J].地震學(xué)報，2（1）：21-30.

段虎榮，康明哲，吳紹宇，等.2020.利用 GRACE時變重力場反演青藏高原的隆升速率[J].地球物理學(xué)報，63（12）：4345-4360.

李強(qiáng)，游新兆，楊少敏，等.2012.中國大陸構(gòu)造變形高精度大密度GPS監(jiān)測——現(xiàn)今速度場[J].中國科學(xué)：地球科學(xué)，42（5）：629-632.

李瑞浩.1988.重力學(xué)引論[M].北京：地震出版社.

申重陽.2005.地殼形變與密度變化耦合運(yùn)動探析[J].大地測量與地球動力學(xué)，25（3）：11-16.

申重陽，李輝.2007.研究現(xiàn)今地殼運(yùn)動和強(qiáng)震機(jī)理的一種方法[J].地球物理學(xué)進(jìn)展，22（1）：49-56.

王嘉沛.2015.地殼變形與密度變化耦合運(yùn)動引起的重力變化效應(yīng)研究[D].北京：中國地震局地震研究所.

王嘉沛，申重陽，玄松柏.2015.全球地殼模型CRUST1.0在青藏高原東南部的重力檢核[J].大地測量與地球動力學(xué)，35（4）：621-626.

邢樂林，王林海，胡敏章，等.2017.時變重力測量確定青藏高原地殼隆升與增厚速率[J].武漢大學(xué)學(xué)報（信息科學(xué)版），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.

基金項目：國家自然科學(xué)基金項目（41674018）資助.

第一作者簡介：梁晨昊（1996-），碩士，主要從事形變重力解釋方面研究.E-mail：lch158768@163.com

通訊作者簡介：申重陽（1963-），研究員，主要從事形變重力解釋、地震前兆機(jī)理、地球物理場及其反演、工程地震等方面的研究.E-mail：scy907@163.com.