基于直立長方體模型的山體隆升與重力變化關(guān)系數(shù)值模擬

梁晨昊 申重陽 王嘉沛

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

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

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

0 引言

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

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

1 基礎(chǔ)理論

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

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

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

1.3 模擬山體形變模型

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

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

Z=4.45H（14）

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

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

2 模擬計(jì)算方法

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

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

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

3 模擬計(jì)算結(jié)果

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

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

4 結(jié)論

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

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

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

參考文獻(xiàn)：

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

段虎榮，康明哲，吳紹宇，等.2020.利用 GRACE時(shí)變重力場反演青藏高原的隆升速率[J].地球物理學(xué)報(bào)，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)動(dòng)探析[J].大地測量與地球動(dòng)力學(xué)，25（3）：11-16.

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

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

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

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

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

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

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