李忠亞， 胡敏章， 王勇， 汪健2,， 申重陽， 李輝

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

0 引言

印度和歐亞板塊持續(xù)碰撞擠壓，形成了大規(guī)模的青藏高原造山帶.青藏高原是全球海拔最高的高原，平均高程4 km以上，被稱為“世界屋脊”和“地球第三極”(高原及周邊區(qū)域地形見圖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ì)量存在橫向差異時，重力作用會驅(qū)使大質(zhì)量塊向小質(zhì)量塊運動，從而影響地表地形和密度結(jié)構(gòu)橫向分布不均勻的巖石圈的構(gòu)造應(yī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)力場及形變不可忽視的力源(Jones et al., 1996; Flesch et al., 2001).

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

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

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

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

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

1 方法

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

(1)

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

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

(2)

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

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

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

2.1 數(shù)據(jù)

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

2.2 未補償重力勢能計算結(jié)果

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

2.3 考慮均衡補償?shù)闹亓菽苡嬎憬Y(jié)果

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

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

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

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

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

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

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

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

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

3 討論

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

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

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

(12)

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

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

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

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

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

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

青藏高原及周邊區(qū)域GNSS應(yīng)變場較好地反映了該地區(qū)現(xiàn)今構(gòu)造形變.為了分析重力勢能作為一種構(gòu)造運動驅(qū)動力對于巖石圈構(gòu)造形變的貢獻，我們將重力勢能差產(chǎn)生的偏應(yīng)力場和GNSS應(yīng)變場進行比較.偏應(yīng)力場與應(yīng)變場之間的相關(guān)程度可以采用Flesch等(2007)定義的相關(guān)系數(shù)C來描述，相關(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)力場和應(yīng)變場完全匹配，即偏應(yīng)力場與應(yīng)變場的壓縮張裂樣式、偏應(yīng)力張量與應(yīng)變張量相對數(shù)值大小和主軸方向一致；若相關(guān)系數(shù)等于-1，則表明偏應(yīng)力場與應(yīng)變場負相關(guān)；若相關(guān)系數(shù)為0，則表示偏應(yīng)力場與應(yīng)變場不相關(guān).我們分別計算了未補償、艾里均衡和普拉特均衡三種重力勢能差模型產(chǎn)生的偏應(yīng)力場(圖5(a、b、c))與GNSS應(yīng)變場之間相關(guān)系數(shù)(圖8).三種相關(guān)系數(shù)模型結(jié)果沒有顯著差別，這主要是因為均衡調(diào)整改正均未根本改變重力勢能分布規(guī)律，進而相應(yīng)的偏應(yīng)力場和相關(guān)系數(shù)均會十分接近.在青藏高原內(nèi)部，柴達木盆地和川滇塊體中南部，重力勢能差產(chǎn)生的偏應(yīng)力場和GNSS應(yīng)變場呈高相關(guān)性，高原內(nèi)部其余區(qū)域兩者相關(guān)性則較弱，甚至在柴達木—隴西塊體東部和拉薩塊體南部還呈現(xiàn)負相關(guān)特性.青藏高原周緣，如塔里木塊體西北部、阿拉善塊體南部、鄂爾多斯塊體西部、華南塊體西北部和喜馬拉雅主斷裂以南區(qū)域的中部和西部，重力勢能差產(chǎn)生的偏應(yīng)力場和GNSS應(yīng)變場展現(xiàn)較強相關(guān)性.偏應(yīng)力場和應(yīng)變場之間的相關(guān)性主要與兩者的壓縮張裂樣式和主軸方向有關(guān)，對比重力勢能差產(chǎn)生的偏應(yīng)力場和GNSS應(yīng)變場，在相關(guān)性強的區(qū)域，兩者的壓縮張裂樣式和主軸方向均符合較好.相關(guān)性弱甚至負相關(guān)區(qū)域，如柴達木-隴西塊體東部和拉薩塊體南部，偏應(yīng)力分布樣式以張裂為主，而應(yīng)變場最大主方向樣式則是壓縮.

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

重力勢能差產(chǎn)生的偏應(yīng)力場和GNSS應(yīng)變場之間相關(guān)性強，表明重力勢能差對構(gòu)造形變貢獻大，若兩者相關(guān)性弱則表示構(gòu)造變形受其他構(gòu)造驅(qū)動力源(如地幔對流產(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)指出，重力勢能差是青藏高原構(gòu)造變形必不可少的構(gòu)造力源.我們研究結(jié)果進一步表明，在重力勢能數(shù)值較大的區(qū)域，如青藏高原中西部和東北部，重力勢能差對構(gòu)造形變貢獻微乎其微，只有在重力勢能橫向梯度變化劇烈的區(qū)域，如柴達木盆地、川滇塊體中南部和青藏高原周緣，重力勢能差對構(gòu)造形變的影響才是不可或缺.

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

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

4 結(jié)論

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

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

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

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

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

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