青藏高原及周邊區(qū)域巖石圈重力勢(shì)能及其產(chǎn)生的偏應(yīng)力場(chǎng)

李忠亞， 胡敏章， 王勇， 汪健2,， 申重陽(yáng)， 李輝

1 中國(guó)科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院， 大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 武漢 430077 2 中國(guó)科學(xué)院大學(xué)， 北京 100049 3 中國(guó)地震局地震研究所， 中國(guó)地震局地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室， 武漢 430071 4 引力與固體潮國(guó)家野外觀測(cè)研究站， 武漢 430071 5 湖北省地震局， 武漢 430071

0 引言

印度和歐亞板塊持續(xù)碰撞擠壓，形成了大規(guī)模的青藏高原造山帶.青藏高原是全球海拔最高的高原，平均高程4 km以上，被稱(chēng)為“世界屋脊”和“地球第三極”(高原及周邊區(qū)域地形見(jiàn)圖1).地震學(xué)和重力學(xué)研究表明，青藏高原內(nèi)部地殼較厚，高原周邊區(qū)域地殼較薄，且相同深度高原內(nèi)外密度存在差異(滕吉文，1979；滕吉文等，1983；Molnar, 1988; 沈旭章等，2013；李紅蕾等，2017；Shen et al., 2016; Tian et al., 2021).當(dāng)巖石圈質(zhì)量存在橫向差異時(shí)，重力作用會(huì)驅(qū)使大質(zhì)量塊向小質(zhì)量塊運(yùn)動(dòng)，從而影響地表地形和密度結(jié)構(gòu)橫向分布不均勻的巖石圈的構(gòu)造應(yīng)力場(chǎng)和構(gòu)造變形(Artyushkov, 1973; Fleitout and Froidevoux, 1982, 1983; Coblentz et al., 1994; Sandiford and Coblentz, 1994; Iaffaldano et al., 2006).因此，青藏高原及周邊地區(qū)巖石圈顯著的橫向密度和地表地形變化是產(chǎn)生其構(gòu)造應(yīng)力場(chǎng)及形變不可忽視的力源(Jones et al., 1996; Flesch et al., 2001).

圖1 青藏高原及周邊區(qū)域地形，黑線為二級(jí)活動(dòng)地塊邊界，紅線表示喜馬拉雅主推覆斷裂，藍(lán)色箭頭表示水平GNSS速度(Zheng et al., 2017)

單位面積巖石圈柱體重力勢(shì)能為單位面積質(zhì)量和其相對(duì)參考深度(一般取巖石圈厚度)的高度的乘積沿整個(gè)巖石圈積分(Jones et al., 1996)，其大小與巖石圈質(zhì)量分布和參考位置有關(guān).巖石圈重力勢(shì)能差即為相鄰柱體間因巖石圈地表地形和密度橫向變化而引起的重力勢(shì)能之差.巖石圈重力勢(shì)能可以根據(jù)大地水準(zhǔn)面高數(shù)據(jù)和巖石圈密度結(jié)構(gòu)數(shù)據(jù)計(jì)算.Turcotte和Schubert(1982)推導(dǎo)了重力勢(shì)能差和大地水準(zhǔn)面高之間簡(jiǎn)單解析公式，該公式在完全均衡條件下成立，且計(jì)算結(jié)果與大地水準(zhǔn)面的濾波方法關(guān)系密切(Ghosh et al., 2009).Coblentz 等(1994)給出了不同巖石圈類(lèi)型的重力勢(shì)能解析計(jì)算公式，后來(lái)隨著多個(gè)地殼模型發(fā)布，現(xiàn)在最常用的方法是根據(jù)重力勢(shì)能定義采用密度結(jié)構(gòu)數(shù)據(jù)計(jì)算.現(xiàn)有關(guān)于青藏高原及周邊區(qū)域重力勢(shì)能計(jì)算研究，大多采用Crust 2.0地殼模型和EGM96地球重力場(chǎng)模型(Flesch et al.，2001; Ghosh et al., 2006, 2009, 2013)，且很多研究是從全球尺度來(lái)討論，我們十分必要采用相對(duì)較新的Crust 1.0地殼模型和EGM2008地球重力場(chǎng)模型來(lái)研究青藏高原及周緣重力勢(shì)能分布規(guī)律.

在地表地形和巖石圈密度結(jié)構(gòu)橫向分布不均勻引起的構(gòu)造應(yīng)力場(chǎng)計(jì)算方面，薄板模型是處理地球重力場(chǎng)與構(gòu)造應(yīng)力場(chǎng)之間關(guān)系的常用模型，其基本內(nèi)涵是當(dāng)巖石圈橫向尺度遠(yuǎn)大于垂向尺度時(shí)，可將巖石圈近似為薄板.國(guó)際上許多學(xué)者是通過(guò)求解薄板模型下重力勢(shì)能和偏應(yīng)力之間平衡方程來(lái)研究重力場(chǎng)和構(gòu)造應(yīng)力場(chǎng)之間的關(guān)系(Flesch et al., 2000, 2001, 2007; Ghosh et al., 2006, 2009; Jay et al., 2017)，盡管該方法計(jì)算量大，求解復(fù)雜，然而其平衡方程求解巧妙，不需要流變信息作為先驗(yàn)信息，理論公式嚴(yán)密，這些使得該方法成為處理地球重力場(chǎng)與構(gòu)造應(yīng)力場(chǎng)問(wèn)題的常用方法.國(guó)內(nèi)學(xué)者提出了采用重力異常數(shù)據(jù)(如布格異常)直接計(jì)算不同深度構(gòu)造應(yīng)力的方法(游永雄，1994；向文和李輝，1999；郭飛霄等，2015；毛經(jīng)倫等，2019；Xu et al., 2020)，但該方法采用了過(guò)多假設(shè)近似，大大限制了其應(yīng)用范圍.

精密的現(xiàn)代大地測(cè)量技術(shù)(如GNSS技術(shù))是研究地殼構(gòu)造變形的有效手段.Wang等(2001)采用中國(guó)大陸GPS觀測(cè)數(shù)據(jù)，首次給出青藏高原及周邊構(gòu)造形變特征.隨著GNSS技術(shù)的發(fā)展、相關(guān)觀測(cè)網(wǎng)的布設(shè)、測(cè)站陸續(xù)增多和觀測(cè)數(shù)據(jù)不斷累積，許多關(guān)于青藏高原及其周邊構(gòu)造形變研究成果相繼發(fā)表(Gan et al., 2007; Liang et al., 2013; Zheng et al., 2017; Jin et al., 2019; Zhang et al., 2019)，豐富的GNSS觀測(cè)資料和研究成果為構(gòu)建青藏高原及鄰區(qū)現(xiàn)今應(yīng)變場(chǎng)提供了堅(jiān)實(shí)基礎(chǔ).構(gòu)造大地測(cè)量學(xué)研究結(jié)果表明，巖石圈重力勢(shì)能差和板塊邊界作用力是青藏高原構(gòu)造形變作用力的重要來(lái)源(Flesch et al., 2000, 2001; Thatcher, 2009).但已有研究沒(méi)有定量分析重力勢(shì)能差產(chǎn)生的構(gòu)造應(yīng)力場(chǎng)和現(xiàn)今GNSS應(yīng)變場(chǎng)之間的關(guān)系，因此本文將探討重力勢(shì)能差對(duì)青藏高原及周緣構(gòu)造形變的貢獻(xiàn)，以便深入認(rèn)識(shí)該區(qū)域構(gòu)造變形方式.

本文基于Crust 1.0地殼模型數(shù)據(jù)，并在其基礎(chǔ)上通過(guò)均衡調(diào)整機(jī)制約束地幔密度，計(jì)算青藏高原及鄰區(qū)重力勢(shì)能.以研究區(qū)域平均重力勢(shì)能作為參考構(gòu)造狀態(tài)重力勢(shì)能，構(gòu)建重力勢(shì)能差模型.利用有限元方法求解重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng).接著，討論地球重力場(chǎng)和構(gòu)造應(yīng)力場(chǎng)之間關(guān)系，并直接利用EGM2008地球重力場(chǎng)模型構(gòu)建重力勢(shì)能差及其產(chǎn)生的偏應(yīng)力場(chǎng).最后，根據(jù)青藏高原及周邊GNSS速度場(chǎng)數(shù)據(jù)計(jì)算應(yīng)變場(chǎng)，分析現(xiàn)今GNSS應(yīng)變和重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)之間相關(guān)性，探討重力勢(shì)能差對(duì)青藏高原及鄰區(qū)構(gòu)造形變的貢獻(xiàn).本文的研究結(jié)果對(duì)深入理解青藏高原及周邊地區(qū)的大地構(gòu)造動(dòng)力學(xué)和構(gòu)造形變機(jī)制等具有參考意義.

1 方法

應(yīng)力和重力之間平衡方程為

(1)

式中ρ為密度，σij為應(yīng)力張量的ij分量(符號(hào)為正表示張裂)，gi為重力加速度第i分量，式中采用愛(ài)因斯坦求和約記.σij可以表示為

σij=τij+1/3(σxx+σyy+σzz)δij,

(2)

式中τij為偏應(yīng)力張量的ij分量，δij為Kronecker符號(hào).

當(dāng)巖石圈橫向尺寸遠(yuǎn)大于其垂向尺寸時(shí)，可將巖石圈抽象成薄板(Jones et al., 1996; 許才軍和尹智，2014).采用薄板模型研究巖石圈密度結(jié)構(gòu)和地表地形橫向不均勻?qū)?gòu)造應(yīng)力場(chǎng)影響時(shí)，可通過(guò)求解平衡方程研究整個(gè)薄板垂向平均偏應(yīng)力場(chǎng)(Flesch et al., 2000, 2001).為此，將式(2)代入式(1)并將各物理量從地表積分至參考深度L，取參考深度的平均值，有

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

2 數(shù)據(jù)與計(jì)算結(jié)果

2.1 數(shù)據(jù)

研究重力勢(shì)能橫向變化對(duì)巖石圈構(gòu)造應(yīng)力場(chǎng)的影響需要地表地形數(shù)據(jù)和巖石圈密度結(jié)構(gòu)數(shù)據(jù).本文采用Crust 1.0地殼模型(Laske et al., 2013)數(shù)據(jù)計(jì)算巖石圈重力勢(shì)能.Crust 1.0地殼模型分辨率為1°×1°，將全球劃分為64800個(gè)網(wǎng)格，每個(gè)網(wǎng)格將地殼分成8層，即水層、冰層、上沉積層、中沉積層、下沉積層、上結(jié)晶地殼層、中結(jié)晶地殼層和下結(jié)晶地殼層.模型提供各層的邊界數(shù)據(jù)和密度數(shù)據(jù)，同時(shí)提供莫霍面以下地幔密度，本文中將該數(shù)據(jù)作為地幔初始參考值.為了與Crust 1.0模型分辨率一致，本文地表地形數(shù)據(jù)采用該模型提供的第一層上邊界數(shù)據(jù).

2.2 未補(bǔ)償重力勢(shì)能計(jì)算結(jié)果

重力勢(shì)能計(jì)算時(shí)，巖石圈參考深度L選為100 km(England and Molnar, 1997; Flesch et al., 2001; Ghosh et al., 2006)，在Crust 1.0地殼模型每層內(nèi)部，密度為常數(shù).莫霍面至參考深度間的地幔部分，其密度仍然采用Crust 1.0模型提供的參考數(shù)據(jù).我們根據(jù)Crust 1.0地殼模型數(shù)據(jù)，首先積分計(jì)算全球每個(gè)格網(wǎng)中心點(diǎn)的重力勢(shì)能，其余點(diǎn)重力勢(shì)能通過(guò)格網(wǎng)中心點(diǎn)樣條插值得到.直接根據(jù)Crust 1.0地殼模型數(shù)據(jù)計(jì)算得到的巖石圈重力勢(shì)能分布見(jiàn)圖2，由于該重力勢(shì)能計(jì)算過(guò)程中未對(duì)Crust 1.0地殼模型提供的地幔參考密度進(jìn)行改正，稱(chēng)該計(jì)算結(jié)果為未補(bǔ)償模型.圖2顯示，青藏高原及周邊區(qū)域重力勢(shì)能總體分布極其不均勻，高重力勢(shì)能幾乎全部集中在青藏高原內(nèi)部，青藏高原周邊區(qū)域重力勢(shì)能相對(duì)較低，并且在青藏高原邊界附近出現(xiàn)高-低重力勢(shì)能變化梯度帶.重力勢(shì)能分布特征還呈現(xiàn)與地形正相關(guān)特性，即地形高區(qū)域重力勢(shì)能大，反之則為低重力勢(shì)能分布.研究區(qū)域內(nèi)高重力勢(shì)能分布主要集中在柴達(dá)木—隴西塊體和巴顏喀拉塊體西部，最大重力勢(shì)能達(dá)到了1579.99 MPa，在拉薩塊體西南部和柴達(dá)木—隴西塊體東部出現(xiàn)次高重力勢(shì)能分布，但在柴達(dá)木盆地、羌塘塊體東部和川滇塊體南部出現(xiàn)了青藏高原內(nèi)部相對(duì)低重力勢(shì)能分布.上述分布使青藏高原內(nèi)部呈現(xiàn)了以高重力勢(shì)能分布為主、高低重力勢(shì)能相間的特征，但就青藏高原內(nèi)部總體分布而言，大致為由西向東重力勢(shì)能逐漸減少.在塔里木塊體、阿拉善塊體、鄂爾多斯塊體西部、華南塊體西北部和印度板塊等青藏高原周邊區(qū)域重力勢(shì)能較低，最低重力勢(shì)能位于印度板塊內(nèi)喜馬拉雅主斷裂西部，其數(shù)值為1420.24 MPa.

2.3 考慮均衡補(bǔ)償?shù)闹亓?shì)能計(jì)算結(jié)果

本節(jié)將討論采用均衡補(bǔ)償機(jī)制對(duì)Crust 1.0地殼模型提供的地幔初始參考密度數(shù)據(jù)進(jìn)行約束改正.常見(jiàn)均衡調(diào)整機(jī)制有兩種方式，即艾里均衡補(bǔ)償模式和普拉特均衡補(bǔ)償模式.前者調(diào)整方式是海平面以上地形部分盈余質(zhì)量(或海平面以下虧損質(zhì)量)以山根(或反山根)形式補(bǔ)償；普拉特模式則是假設(shè)補(bǔ)償面地形沒(méi)有起伏，通過(guò)改變柱體密度進(jìn)行均衡補(bǔ)償.我們?cè)趫D2未補(bǔ)償重力勢(shì)能的基礎(chǔ)上分別計(jì)算艾里均衡調(diào)整和普拉特均衡調(diào)整補(bǔ)償后的重力勢(shì)能.

圖2 青藏高原及周邊區(qū)域未補(bǔ)償重力勢(shì)能

計(jì)算考慮艾里均衡補(bǔ)償?shù)闹亓?shì)能，國(guó)際上有兩種策略：(1)Flesch等(Flesch et al., 2000, 2007; Flesch and Kreemer, 2010)和Jay等(2017)研究青藏高原、北美西部和智利中部重力勢(shì)能分布時(shí)，艾里均衡補(bǔ)償對(duì)重力勢(shì)能的影響是通過(guò)山根或反山根方式計(jì)算；(2)Ghosh等(2009)和Neves等(2014)計(jì)算艾里均衡補(bǔ)償?shù)闹亓?shì)能時(shí)，在地表地形中扣除動(dòng)力地形部分貢獻(xiàn)，該方法的不足是容易混淆均衡補(bǔ)償和巖石圈下部密度差作用(Naliboff et al., 2012).我們計(jì)算艾里均衡補(bǔ)償?shù)闹亓?shì)能模型時(shí)，采用前者，即青藏高原及周緣地形部分盈余質(zhì)量通過(guò)地殼底部插入地幔部分的山根進(jìn)行補(bǔ)償，山根的深度采用Crust 1.0地殼模型提供的地形數(shù)據(jù)、上沉積層密度、下結(jié)晶地殼密度和地幔密度等數(shù)據(jù)計(jì)算.計(jì)算重力勢(shì)能時(shí)山根的密度采用艾里均衡補(bǔ)償后的數(shù)據(jù)，其余數(shù)據(jù)仍采用Crust 1.0地殼模型結(jié)果，由此計(jì)算的重力勢(shì)能分布見(jiàn)圖3，我們稱(chēng)該結(jié)果為艾里均衡補(bǔ)償重力勢(shì)能模型.與圖2中未補(bǔ)償重力勢(shì)能模型相比，圖3重力勢(shì)能數(shù)值稍許減小，其變化區(qū)間為1413.22～1533.80 MPa，平均重力勢(shì)能為1476.04 MPa.艾里均衡補(bǔ)償?shù)幕居^點(diǎn)是海平面以上質(zhì)量被其下面低密度物質(zhì)所補(bǔ)償，因此重力勢(shì)能計(jì)算中加上艾里均衡補(bǔ)償后數(shù)值必然降低.但是，考慮艾里均衡調(diào)整后，重力勢(shì)能整體分布規(guī)律并未改變，仍然是青藏高原內(nèi)部為高重力勢(shì)能分布區(qū)，周邊區(qū)域重力勢(shì)能較低.

計(jì)算考慮普拉特均衡補(bǔ)償?shù)闹亓?shì)能時(shí)，我們的計(jì)算方法是將海平面以上的盈余質(zhì)量(或海平面以下虧損質(zhì)量)通過(guò)調(diào)整地殼底部至參考深度L間的地幔密度來(lái)補(bǔ)償，然后采用均衡調(diào)整后的地幔密度計(jì)算重力勢(shì)能，計(jì)算結(jié)果見(jiàn)圖4，我們稱(chēng)該結(jié)果為普拉特均衡補(bǔ)償重力勢(shì)能模型.與未補(bǔ)償重力勢(shì)能模型相比，普拉特均衡補(bǔ)償模型重力勢(shì)能分布特征與其基本一致，但補(bǔ)償后重力勢(shì)能數(shù)值略微下降，其平均值為1481.50 MPa，變化區(qū)間為1415.00～1552.15 MPa，造成該種變化的原因是普拉特均衡調(diào)整機(jī)制通過(guò)地殼底部至參考深度間低密度補(bǔ)償海平面以上盈余質(zhì)量，因此會(huì)導(dǎo)致重力勢(shì)能數(shù)值比未補(bǔ)償模型小.比較圖3和圖4可以發(fā)現(xiàn)，兩種均衡調(diào)整補(bǔ)償后重力勢(shì)能分布特征基本一致，但是艾里均衡補(bǔ)償對(duì)于重力勢(shì)能數(shù)值影響更大.

2.4 青藏高原重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)

求解重力勢(shì)能橫向變化產(chǎn)生的偏應(yīng)力場(chǎng)，首先需要計(jì)算重力勢(shì)能差，即重力勢(shì)能與構(gòu)造參考狀態(tài)重力勢(shì)能之差(Flesch et al.,2000,2001; Ghosh et al., 2009).構(gòu)造參考狀態(tài)重力勢(shì)能可采用研究區(qū)域重力勢(shì)能模型的平均值(Coblentz et al., 1994)，將重力勢(shì)能模型數(shù)值減去構(gòu)造參考狀態(tài)重力勢(shì)能即得重力勢(shì)能差.根據(jù)前面構(gòu)建的三種重力勢(shì)能模型(圖2、圖3和圖4)分別計(jì)算，得到重力勢(shì)能差，采用FEMCalculation程序(李忠亞，2017)計(jì)算了對(duì)應(yīng)的重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)(圖5).未補(bǔ)償重力勢(shì)能差產(chǎn)生的偏應(yīng)力(圖5a)分布特征整體表現(xiàn)為地形高的青藏高原內(nèi)部為張裂應(yīng)力，青藏高原周邊的低地形區(qū)域(如塔里木盆地、四川盆地和印度板塊北部)則是壓縮應(yīng)力.青藏高原及周緣偏應(yīng)力數(shù)值主要集中在10～30 MPa區(qū)間內(nèi)，最大主張應(yīng)力方向在青藏高原西部為NNE-SSW，青藏高原中部(拉薩塊體中北部與羌塘塊體中部和巴顏喀拉塊體中部)幾乎接近N-S分布，而在青藏高原東部則較復(fù)雜，具體表現(xiàn)為東北部由高原中部接近N-S樣式逐漸過(guò)渡為NNE-SSW，川滇塊體、滇南塊體北部和滇西塊體北部接近E-W方向，在滇南塊體南部主要表現(xiàn)為NEE-SWW.在塔里木盆地，重力勢(shì)能差產(chǎn)生的最大主壓應(yīng)力方向由盆地東邊NNE-SSW轉(zhuǎn)變至盆地西邊的NNW-SSE，四川盆地最大主壓應(yīng)力方向?yàn)镹WW-SEE.印度板塊東北部至西北部，最大主壓應(yīng)力方向由近N-S 逐漸過(guò)渡至NW-SE.將艾里均衡補(bǔ)償和普拉特均衡補(bǔ)償兩種重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)(圖5b 和5c)與未補(bǔ)償重力勢(shì)能差得到的偏應(yīng)力場(chǎng)比較可知，三種偏應(yīng)力場(chǎng)分布樣式幾乎一致，但是經(jīng)均衡調(diào)整改正的偏應(yīng)力場(chǎng)應(yīng)力數(shù)值變小，艾里均衡補(bǔ)償和普拉特均衡補(bǔ)償兩種偏應(yīng)力場(chǎng)模型的最大偏應(yīng)力數(shù)值大小分別是25.96 MPa和25.30 MPa.與未補(bǔ)償重力勢(shì)能模型相比，均衡調(diào)整使重力勢(shì)能數(shù)值變小并降低重力勢(shì)能水平變化幅度，最終導(dǎo)致均衡補(bǔ)償模型求解的偏應(yīng)力比未補(bǔ)償模型結(jié)果減小.

圖3 青藏高原及周邊區(qū)域艾里均衡補(bǔ)償后的重力勢(shì)能

圖4 青藏高原及周邊區(qū)域普拉特均衡補(bǔ)償后的重力勢(shì)能

圖5 青藏高原及周邊區(qū)域重力勢(shì)能差產(chǎn)生的偏應(yīng)力

3 討論

3.1 地球重力場(chǎng)和構(gòu)造應(yīng)力場(chǎng)之間關(guān)系

地球重力是由地球?qū)ξ矬w的引力和地球自轉(zhuǎn)產(chǎn)生的離心力的合力，空間域中每一點(diǎn)均有惟一的重力與之對(duì)應(yīng)，即地球重力場(chǎng).地球重力場(chǎng)反映了地球系統(tǒng)中物質(zhì)的分布和運(yùn)動(dòng)狀態(tài)，包括地球內(nèi)部質(zhì)量、密度的分布和變化和地球物質(zhì)空間分布、運(yùn)動(dòng)和變化.構(gòu)造應(yīng)力場(chǎng)表示一定范圍內(nèi)某一瞬時(shí)的應(yīng)力狀態(tài)(徐開(kāi)禮和朱志澄，1989)，其時(shí)空變化是導(dǎo)致地殼形變、破裂、褶皺等地質(zhì)現(xiàn)象的直接動(dòng)因(姚瑞等，2017).地球重力場(chǎng)和巖石圈構(gòu)造應(yīng)力場(chǎng)之間聯(lián)系是：巖石圈密度結(jié)構(gòu)和地表地形橫向不均勻，在地球重力場(chǎng)作用下可以產(chǎn)生橫向作用力形成構(gòu)造應(yīng)力，這種橫向作用力在薄板模型下可以用重力勢(shì)能差定量描述(Turcotte and Schubert, 1982; Jones et al., 1996; Flesch et al., 2000, 2001; Ghosh et al., 2009).國(guó)內(nèi)在20世紀(jì)90年代，游永雄(1994)給出了由布格重力異常計(jì)算構(gòu)造應(yīng)力場(chǎng)的方法.此后，向文和李輝(1999)和Xu等(2020)對(duì)該方法進(jìn)行了發(fā)展，但核心思路依然是采用重力異常(如均衡異常)計(jì)算構(gòu)造應(yīng)力；郭飛霄等(2015)和毛經(jīng)倫等(2019)采用該方法計(jì)算了川西地區(qū)水平構(gòu)造應(yīng)力場(chǎng).然而，該方法數(shù)理模型不夠嚴(yán)謹(jǐn)，具體表現(xiàn)為：用矢量來(lái)表示不同深度水平構(gòu)造應(yīng)力；將重力水平矢量方向來(lái)近似代替構(gòu)造應(yīng)力最大主方向，無(wú)法表示壓縮應(yīng)力和張裂應(yīng)力；采用布格重力或均衡重力異常無(wú)法考慮地表地形橫向差異對(duì)構(gòu)造應(yīng)力場(chǎng)的貢獻(xiàn).許才軍和尹智(2014)對(duì)國(guó)內(nèi)外利用重力數(shù)據(jù)計(jì)算構(gòu)造應(yīng)力場(chǎng)的研究進(jìn)展進(jìn)行了綜述，建議構(gòu)建彈性、黏彈性和彈塑性的利用重力數(shù)據(jù)計(jì)算地殼構(gòu)造應(yīng)力應(yīng)變場(chǎng)的解析模型.本文第1節(jié)中闡述了地球重力場(chǎng)和構(gòu)造應(yīng)力場(chǎng)之間理論關(guān)系，該方法從重力與應(yīng)力平衡方程出發(fā)，基于薄板模型，建立了偏應(yīng)力和重力勢(shì)能之間的平衡方程，并在偏應(yīng)力第二不變量最小約束下求解平衡方程，得到偏應(yīng)力場(chǎng)，該方法是研究地球重力場(chǎng)和構(gòu)造應(yīng)力場(chǎng)之間關(guān)系的嚴(yán)密方法.

重力勢(shì)能差可以采用大地水準(zhǔn)面高估算，從而可由重力場(chǎng)模型通過(guò)求解第1節(jié)的平衡方程得到偏應(yīng)力場(chǎng).下面將利用EGM2008地球重力場(chǎng)模型計(jì)算青藏高原及周邊區(qū)域重力勢(shì)能差及其產(chǎn)生的偏應(yīng)力場(chǎng)，并將計(jì)算結(jié)果與第3節(jié)中Crutst 1.0 地殼模型結(jié)果進(jìn)行比較.大地水準(zhǔn)面高ΔN和重力勢(shì)能差ΔU之間關(guān)系(Sandiford and Coblentz, 1994)為

(12)

式中，g為重力加速度，G為萬(wàn)有引力常數(shù).為與第1節(jié)中定義的重力勢(shì)能一致，通過(guò)式(12)計(jì)算的重力勢(shì)能差需要除以巖石圈參考深度100 km.由于Crust 1.0地殼模型分辨率為1°×1°，因此我們將EGM2008地球重力場(chǎng)模型的最大階次數(shù)截取至200.重力場(chǎng)模型低階次項(xiàng)包含巖石圈下面地球物質(zhì)對(duì)大地水準(zhǔn)面的影響，需要扣除其貢獻(xiàn)，但兩者之間沒(méi)有一一具體對(duì)應(yīng)關(guān)系.Ghosh等(2009)計(jì)算全球重力勢(shì)能差時(shí)減去7階次以下系數(shù)，Neves等(2014)研究伊比利亞地區(qū)重力勢(shì)能差時(shí)是去除12階次以下系數(shù).本文將待去除的階次系數(shù)設(shè)置在區(qū)間[5,15]，然后逐一試算.以高地形重力勢(shì)能差為正、低地形為負(fù)作為篩選條件，最終確定去除10階次以下系數(shù).采用上述策略，由大地水準(zhǔn)面高計(jì)算的重力勢(shì)能差及其產(chǎn)生的偏應(yīng)力場(chǎng)見(jiàn)圖6.青藏高原南部為高重力勢(shì)能差分布區(qū)域，最大重力勢(shì)能差達(dá)到了45.367 MPa，塔里木盆地、柴達(dá)木盆地、四川盆地和印度板塊北部均為低重力勢(shì)能差區(qū)域.青藏高原內(nèi)部正重力勢(shì)能差產(chǎn)生張裂應(yīng)力，負(fù)重力勢(shì)能差則是壓縮應(yīng)力，應(yīng)力大小主要分布在1～20 MPa.

對(duì)比圖5和圖6可知，整體上EGM2008地球重力場(chǎng)模型計(jì)算的重力勢(shì)能差與Crust 1.0地殼模型結(jié)果較相似，即青藏高原內(nèi)部除柴達(dá)木盆地外均為正重力勢(shì)能差，青藏高原周邊大部分區(qū)域重力勢(shì)能差小于青藏高原內(nèi)部.但是，兩者具體細(xì)節(jié)分布規(guī)律存在較大差異，如EGM2008 地球重力場(chǎng)模型結(jié)果最大重力勢(shì)能差主要分布在青藏高原南端，并向北遞減，而Crust 1.0模型在青藏高原西北部和東北部均存在高重力勢(shì)能差分布.兩種模型結(jié)果重力勢(shì)能差分布差異導(dǎo)致偏應(yīng)力分布存在相應(yīng)的區(qū)別.造成兩者細(xì)節(jié)差異較大的主要原因是，式(12)僅在滿足完全均衡條件下才成立(Turcotte and Schubert, 1982; Ghosh et al., 2009)，實(shí)測(cè)重力結(jié)果顯示青藏高原并非處于完全均衡狀態(tài)(張永謙等，2010；陳石等，2011；付廣裕等，2015).此外，重力場(chǎng)模型系數(shù)的選取對(duì)結(jié)果有一定影響，但難以在球諧系數(shù)中完全分離出巖石圈物質(zhì).綜上，式(12)表示的方法使用時(shí)前提條件較苛刻且對(duì)重力場(chǎng)系數(shù)選取存在難度，因此該方法只能作為重力勢(shì)能差的一種近似估算方法.但該方法可直接得到重力勢(shì)能差，在實(shí)際計(jì)算中可作為一種獨(dú)立方法，用于與其他結(jié)果比較.

圖6 利用EGM2008地球重力場(chǎng)模型計(jì)算的青藏高原及鄰區(qū)重力勢(shì)能差(a)及其產(chǎn)生的偏應(yīng)力(b)

3.2 重力勢(shì)能差產(chǎn)生的偏應(yīng)力與GNSS現(xiàn)今構(gòu)造形變

現(xiàn)代空間大地測(cè)量觀測(cè)技術(shù)(如GNSS技術(shù))能夠精確測(cè)量板塊之間相對(duì)運(yùn)動(dòng)，其觀測(cè)結(jié)果可揭示板塊之間和板塊內(nèi)部的構(gòu)造變形方式.我們收集Zheng等(2017)和Kreemer等(2014)發(fā)表的青藏高原及其周邊共計(jì)2576個(gè)測(cè)站GNSS速度場(chǎng)資料(圖1)，利用Shen等(1996)提出的應(yīng)變計(jì)算方法，采用SSPX程序(Cardozo and Allmendinger, 2009)計(jì)算了青藏高原及周緣現(xiàn)今應(yīng)變場(chǎng).計(jì)算結(jié)果(圖7)表明，青藏高原及其周邊區(qū)域GNSS現(xiàn)今構(gòu)造應(yīng)變場(chǎng)和該地區(qū)長(zhǎng)期地質(zhì)構(gòu)造背景相符合.在喜馬拉雅和龍門(mén)山等地區(qū)以壓縮應(yīng)變率為優(yōu)勢(shì)分布，與這些造山帶的地殼壓縮變形一致；青藏高原內(nèi)部羌塘塊體中南部和拉薩塊體中北部存在近東西向張裂應(yīng)變，與這些區(qū)域發(fā)育的正斷層分布一致；青藏高原周邊的華南塊體和塔里木塊體作為穩(wěn)定的構(gòu)造單元，其應(yīng)變明顯小于青藏高原內(nèi)部.

圖7 青藏高原及周邊區(qū)域GNSS主應(yīng)變率

青藏高原及周邊區(qū)域GNSS應(yīng)變場(chǎng)較好地反映了該地區(qū)現(xiàn)今構(gòu)造形變.為了分析重力勢(shì)能作為一種構(gòu)造運(yùn)動(dòng)驅(qū)動(dòng)力對(duì)于巖石圈構(gòu)造形變的貢獻(xiàn)，我們將重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)和GNSS應(yīng)變場(chǎng)進(jìn)行比較.偏應(yīng)力場(chǎng)與應(yīng)變場(chǎng)之間的相關(guān)程度可以采用Flesch等(2007)定義的相關(guān)系數(shù)C來(lái)描述，相關(guān)系數(shù)公式為

(13)

E和T分別表示應(yīng)變和偏應(yīng)力第二不變量，ε為GNSS應(yīng)變張量.相關(guān)系數(shù)C數(shù)值變化范圍為[-1,1].相關(guān)系數(shù)等于1，表示偏應(yīng)力場(chǎng)和應(yīng)變場(chǎng)完全匹配，即偏應(yīng)力場(chǎng)與應(yīng)變場(chǎng)的壓縮張裂樣式、偏應(yīng)力張量與應(yīng)變張量相對(duì)數(shù)值大小和主軸方向一致；若相關(guān)系數(shù)等于-1，則表明偏應(yīng)力場(chǎng)與應(yīng)變場(chǎng)負(fù)相關(guān)；若相關(guān)系數(shù)為0，則表示偏應(yīng)力場(chǎng)與應(yīng)變場(chǎng)不相關(guān).我們分別計(jì)算了未補(bǔ)償、艾里均衡和普拉特均衡三種重力勢(shì)能差模型產(chǎn)生的偏應(yīng)力場(chǎng)(圖5(a、b、c))與GNSS應(yīng)變場(chǎng)之間相關(guān)系數(shù)(圖8).三種相關(guān)系數(shù)模型結(jié)果沒(méi)有顯著差別，這主要是因?yàn)榫庹{(diào)整改正均未根本改變重力勢(shì)能分布規(guī)律，進(jìn)而相應(yīng)的偏應(yīng)力場(chǎng)和相關(guān)系數(shù)均會(huì)十分接近.在青藏高原內(nèi)部，柴達(dá)木盆地和川滇塊體中南部，重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)和GNSS應(yīng)變場(chǎng)呈高相關(guān)性，高原內(nèi)部其余區(qū)域兩者相關(guān)性則較弱，甚至在柴達(dá)木—隴西塊體東部和拉薩塊體南部還呈現(xiàn)負(fù)相關(guān)特性.青藏高原周緣，如塔里木塊體西北部、阿拉善塊體南部、鄂爾多斯塊體西部、華南塊體西北部和喜馬拉雅主斷裂以南區(qū)域的中部和西部，重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)和GNSS應(yīng)變場(chǎng)展現(xiàn)較強(qiáng)相關(guān)性.偏應(yīng)力場(chǎng)和應(yīng)變場(chǎng)之間的相關(guān)性主要與兩者的壓縮張裂樣式和主軸方向有關(guān)，對(duì)比重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)和GNSS應(yīng)變場(chǎng)，在相關(guān)性強(qiáng)的區(qū)域，兩者的壓縮張裂樣式和主軸方向均符合較好.相關(guān)性弱甚至負(fù)相關(guān)區(qū)域，如柴達(dá)木-隴西塊體東部和拉薩塊體南部，偏應(yīng)力分布樣式以張裂為主，而應(yīng)變場(chǎng)最大主方向樣式則是壓縮.

圖8 青藏高原及周緣重力勢(shì)能差產(chǎn)生的偏應(yīng)力和GNSS應(yīng)變之間相關(guān)系數(shù)

重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)和GNSS應(yīng)變場(chǎng)之間相關(guān)性強(qiáng)，表明重力勢(shì)能差對(duì)構(gòu)造形變貢獻(xiàn)大，若兩者相關(guān)性弱則表示構(gòu)造變形受其他構(gòu)造驅(qū)動(dòng)力源(如地幔對(duì)流產(chǎn)生的牽引力和邊界載荷作用等因素)所控制.大量研究結(jié)果(張健和石耀霖, 2002; England and Molnar, 1997; Flesch et al., 2001; Ghosh et al., 2006; Ghosh and Holt, 2012; Schmalholz et al., 2014, 2019)指出，重力勢(shì)能差是青藏高原構(gòu)造變形必不可少的構(gòu)造力源.我們研究結(jié)果進(jìn)一步表明，在重力勢(shì)能數(shù)值較大的區(qū)域，如青藏高原中西部和東北部，重力勢(shì)能差對(duì)構(gòu)造形變貢獻(xiàn)微乎其微，只有在重力勢(shì)能橫向梯度變化劇烈的區(qū)域，如柴達(dá)木盆地、川滇塊體中南部和青藏高原周緣，重力勢(shì)能差對(duì)構(gòu)造形變的影響才是不可或缺.

青藏高原周緣的塔里木塊體西北部，阿拉善塊體南部，鄂爾多斯塊體西部、華南塊體西北部和青藏高原內(nèi)部川滇塊體中南部，重力勢(shì)能差產(chǎn)生的偏應(yīng)力和GNSS應(yīng)變相關(guān)性強(qiáng)，板塊運(yùn)動(dòng)驅(qū)動(dòng)力在遠(yuǎn)離歐亞板塊邊界區(qū)域以重力勢(shì)能的形式影響其構(gòu)造形變.上述青藏高原周緣區(qū)域構(gòu)造形變較小，而高原內(nèi)部川滇塊體中南部卻較大，這與兩者巖石圈流變性質(zhì)有關(guān)，前者黏度大，而后者黏度小(石耀霖和曹建玲，2008; Flesch et al., 2001; Jay et al., 2017)，故導(dǎo)致形變差異.喜馬拉雅主斷裂南端，盡管重力勢(shì)能差產(chǎn)生的偏應(yīng)力和GNSS應(yīng)變相關(guān)性強(qiáng)，但是該區(qū)域重力勢(shì)能差產(chǎn)生的壓縮應(yīng)力是北部青藏高原內(nèi)部高重力勢(shì)能壓縮南部低重力勢(shì)能的結(jié)果，而GNSS壓縮應(yīng)變則是南部印度板塊向北部高原俯沖的結(jié)果，兩種機(jī)制作用力方向完全相反.

關(guān)于青藏高原構(gòu)造形變的動(dòng)力學(xué)機(jī)制存在不同觀點(diǎn).England和Molnar(1997)認(rèn)為青藏高原形變遵循重力勢(shì)能差和應(yīng)力之間平衡方程，F(xiàn)lesch等(2001, 2007)聯(lián)合重力勢(shì)能差和板塊邊界作用力兩種力源解釋歐亞碰撞帶的動(dòng)力學(xué)過(guò)程，Ghosh等(Ghosh et al., 2008；Ghosh and Holt, 2012)則認(rèn)為重力勢(shì)能差和地幔對(duì)流聯(lián)合作用控制青藏高原動(dòng)力學(xué)過(guò)程.本文上述關(guān)于重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)和現(xiàn)今GNSS應(yīng)變場(chǎng)之間相關(guān)性分析表明，重力勢(shì)能差在青藏高原及鄰區(qū)不同區(qū)域?qū)π巫冇绊懖煌?在喜馬拉雅地區(qū)，印度板塊可以克服重力勢(shì)能差的影響持續(xù)俯沖，這表明印度板塊驅(qū)動(dòng)力的影響比重力勢(shì)能差更大.因此，青藏高原構(gòu)造形變的動(dòng)力來(lái)源至少包含重力勢(shì)能差和印度板塊驅(qū)動(dòng)力.

4 結(jié)論

本文基于Crust 1.0地殼模型計(jì)算了青藏高原及周邊區(qū)域重力勢(shì)能，并根據(jù)重力勢(shì)能差求解了其產(chǎn)生的偏應(yīng)力場(chǎng).在此基礎(chǔ)上討論了地球重力場(chǎng)和構(gòu)造應(yīng)力場(chǎng)之間關(guān)系和重力勢(shì)能差對(duì)青藏高原及周緣構(gòu)造形變的影響，得出的主要結(jié)論如下：

(1)本文以Crust 1.0地殼模型作為青藏高原及周邊區(qū)域重力勢(shì)能計(jì)算的基本參考模型，在此基礎(chǔ)上分別采用艾里均衡和普拉特均衡兩種方式對(duì)地幔密度進(jìn)行約束，由此構(gòu)建青藏高原及周邊區(qū)域重力勢(shì)能分布模型.未補(bǔ)償、艾里均衡補(bǔ)償與普拉特均衡補(bǔ)償三種方式計(jì)算的重力勢(shì)能變化范圍分別為1420.24～1579.99 MPa，1413.22～1533.80 MPa和1415.00～1552.15 MPa.重力勢(shì)能與地形呈正相關(guān)關(guān)系，青藏高原內(nèi)部重力勢(shì)能普遍較高，高原四周則為低重力勢(shì)能分布.高原西北和西南區(qū)域?yàn)楦咧亓?shì)能聚集區(qū)，印度板塊北部、塔里木盆地、四川盆地和柴達(dá)木盆地是主要低重力勢(shì)能區(qū)域.高原邊界重力勢(shì)能分布表現(xiàn)為正負(fù)急劇變化梯度帶特征.

(2)以研究區(qū)域平均重力勢(shì)能為構(gòu)造參考狀態(tài)重力勢(shì)能，計(jì)算青藏高原及周邊區(qū)域重力勢(shì)能差分布，然后求解其產(chǎn)生的偏應(yīng)力場(chǎng).地形高的青藏高原內(nèi)部重力勢(shì)能差為正，產(chǎn)生偏應(yīng)力以張裂樣式為優(yōu)勢(shì)分布；地形低的青藏高原周邊重力勢(shì)能差為負(fù)，偏應(yīng)力以壓縮為主要分布特征.青藏高原及周邊區(qū)域重力勢(shì)能差產(chǎn)生的偏應(yīng)力大小主要分布在10～30 MPa.

(3)地表地形和巖石圈密度結(jié)構(gòu)存在橫向差異時(shí)，在地球重力場(chǎng)作用下可以產(chǎn)生橫向構(gòu)造應(yīng)力.求解重力勢(shì)能與偏應(yīng)力之間平衡方程是研究地球重力場(chǎng)和構(gòu)造應(yīng)力場(chǎng)之間關(guān)系的有效方法.大地水準(zhǔn)面高和重力勢(shì)能差之間具有簡(jiǎn)單關(guān)系，但其僅在均衡區(qū)域才成立，該關(guān)系可作為研究重力勢(shì)能差的獨(dú)立近似方法.

(4)重力勢(shì)能橫向變化劇烈區(qū)域，重力勢(shì)能差對(duì)于構(gòu)造形變貢獻(xiàn)突出.青藏高原內(nèi)部的柴達(dá)木盆地和川滇塊體中南部，以及青藏高原周緣的塔里木塊體西北部、阿拉善塊體南部、鄂爾多斯塊體西部、華南塊體西北部和喜馬拉雅主斷裂以南區(qū)域的中部和西部，這些地區(qū)重力勢(shì)能差產(chǎn)生的偏應(yīng)力場(chǎng)和GNSS應(yīng)變場(chǎng)相關(guān)性強(qiáng)，重力勢(shì)能差對(duì)構(gòu)造形變的影響大.在遠(yuǎn)離歐亞板塊邊界的青藏高原北部和東部邊緣，板塊運(yùn)動(dòng)驅(qū)動(dòng)力以重力勢(shì)能的形式影響構(gòu)造形變，而在喜馬拉雅地區(qū)，重力勢(shì)能差與印度板塊俯沖兩種機(jī)制作用方向相反.

致謝感謝兩位匿名審稿專(zhuān)家提出的寶貴意見(jiàn)和Ghosh Attreyee博士、Stamps D. Sarah博士和張倩文同學(xué)的有益討論.