青藏高原及周邊區域巖石圈重力勢能及其產生的偏應力場

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

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

0 引言

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

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

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

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

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

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

1 方法

應力和重力之間平衡方程為

(1)

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

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

(2)

式中τij為偏應力張量的ij分量，δij為Kronecker符號.

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

2 數據與計算結果

2.1 數據

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

2.2 未補償重力勢能計算結果

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

2.3 考慮均衡補償的重力勢能計算結果

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

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

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

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

2.4 青藏高原重力勢能差產生的偏應力場

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

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

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

圖5 青藏高原及周邊區域重力勢能差產生的偏應力

3 討論

3.1 地球重力場和構造應力場之間關系

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

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

(12)

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

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

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

3.2 重力勢能差產生的偏應力與GNSS現今構造形變

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

圖7 青藏高原及周邊區域GNSS主應變率

青藏高原及周邊區域GNSS應變場較好地反映了該地區現今構造形變.為了分析重力勢能作為一種構造運動驅動力對于巖石圈構造形變的貢獻，我們將重力勢能差產生的偏應力場和GNSS應變場進行比較.偏應力場與應變場之間的相關程度可以采用Flesch等(2007)定義的相關系數C來描述，相關系數公式為

(13)

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

圖8 青藏高原及周緣重力勢能差產生的偏應力和GNSS應變之間相關系數

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

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

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

4 結論

本文基于Crust 1.0地殼模型計算了青藏高原及周邊區域重力勢能，并根據重力勢能差求解了其產生的偏應力場.在此基礎上討論了地球重力場和構造應力場之間關系和重力勢能差對青藏高原及周緣構造形變的影響，得出的主要結論如下：

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

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

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

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

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