黨亞民，楊 強，王 偉，梁玉可

1. 中國測繪科學研究院，北京 100830； 2. 山東科技大學，山東 青島 266590

青藏高原是世界屋脊，與世界上的其他高原不同，青藏高原地殼具有復雜、不均一、多地體拼貼、巨厚(60～70 km)等特點[1]，是地學研究的寶庫。青藏高原地殼水平方向變形主要呈現(xiàn)東西向拉伸、南北向壓縮趨勢，自南向北運動速度逐漸衰減，高程方向表現(xiàn)為分階段隆升特征[2-3]。

青藏高原的變形與隆升，與印度板塊、歐亞板塊的匯聚碰撞過程密切相關。喜馬拉雅造山帶位于印度板塊與歐亞板塊交接處，是青藏高于最南端，長約2500 km，寬300～500 km，珠穆朗瑪峰(以下簡稱“珠峰”)是其中最高峰，也是我國和世界最高峰。喜馬拉雅造山帶是印度板塊和歐亞板塊匯聚的主要變形帶，是研究大陸構造形變機制最具代表性的區(qū)域。

青藏高原地殼隆升和變形機制目前仍存在很多爭論。文獻[4]提出“雙倍地殼模式”，隨后眾多學者先后提出多種理論和動力學假說，主要有“地殼增厚模式”[5-7]“大陸逃逸模式”[8-11]“雙向俯沖作用模式”[12]，我國科學家也先后提出“多邊大陸俯沖和地幔底辟模式”[13-15]“大陸注入模式”[16-17]等理論假說。上述多種假說均是基于不同地質或地球物理觀測現(xiàn)象提出的青藏高原隆升演化變形模式，目前仍存在爭議。

隨著空間大地測量觀測技術的發(fā)展，GNSS呈現(xiàn)出高精度、大范圍、準實時的優(yōu)點，自1988年我國首次開展GNSS地殼形變監(jiān)測以來，GNSS就越來越廣泛地應用于地殼形變監(jiān)測中，使得地殼運動與形變監(jiān)測逐步實現(xiàn)高精度、連續(xù)動態(tài)監(jiān)測，實現(xiàn)了由定性向定量監(jiān)測的飛躍。青藏高原地殼隆升與變形也逐步實現(xiàn)連續(xù)、定量監(jiān)測。

我國于1988年在滇西地震試驗場，首次將GNSS技術應用到地殼運動與形變監(jiān)測中[18]。1991年，中美合作首先在青藏高原東部及鄰區(qū)，開展了GPS觀測。通過1991—1995年的4期觀測，定量地監(jiān)測了青藏高原東南緣現(xiàn)今構造變形特征[19]。同期，中美還合作開展了跨喜馬拉雅GPS觀測。美方在尼泊爾布設28個GPS站進行了觀測，中方在西藏南部布設了5個GPS站，并聯(lián)合格爾木地震臺進行了同步觀測[20]。

我國獨立開展青藏高原GNSS測量始于1992年，國家攀登計劃項目“現(xiàn)代地殼運動和地球動力學研究”于1992和1994年兩期聯(lián)測在青藏高原獲取了拉薩-溫泉等基線[21]。1993—1995年，國家地震局在“八五”期間，在青藏高原開展了兩期GPS觀測[22]?！熬盼濉逼陂g，我國開展了“中國地殼運動觀測網(wǎng)絡”(CMONOC-I)工程建設，2000年12月通過國家驗收，正式投入運行[23-25]?！笆晃濉逼陂g，在中國地殼運動觀測網(wǎng)絡的基礎上，開展了“中國大陸構造環(huán)境監(jiān)測網(wǎng)絡”(CMONOC-Ⅱ)工程建設，該網(wǎng)絡于2011年正式投入運行，包括260個基準站和近2000個區(qū)域站[24-26]。

1992—2020年，我國與意大利、美國合作，先后于1992、1998、1999年開展了珠峰高程GNSS測量，我國于2005、2020年自主開展了珠峰高程GNSS測量觀測[27-29]。2020年珠峰高程測量中首次在峰頂使用了國產GNSS接收機，利用北斗衛(wèi)星信號進行了峰頂高程測量。

GNSS應用于青藏高原地殼形變監(jiān)測，從多層次、多方面研究了青藏高原地殼形變及其動力學機制[30-45]。文獻[36—37]利用GNSS數(shù)據(jù)，發(fā)現(xiàn)印歐板塊碰撞使青藏高原地殼縮短并增厚，但碰撞90%以上被喜馬拉雅山及其邊緣形變所吸收，從而為地殼增厚學說提供了支持。文獻[38—39]利用青藏高原的553個GNSS站點，研究了青藏高原構造運動特征，發(fā)現(xiàn)青藏高原現(xiàn)今主要構造運動屬于連續(xù)介質變形過程。文獻[40]基于GNSS觀測結果研究了青藏高原的整體運動狀態(tài)，揭示了青藏高原物質運移機理。文獻[41]利用10多年累積的GNSS觀測資料，對青藏高原現(xiàn)今地殼運動進行研究，探討了青藏高原現(xiàn)今隆升狀況及其動力學含義。文獻[42]通過GNSS應變率場，研究了青藏高原內部更加重要的地殼減薄過程。文獻[43]利用GNSS數(shù)據(jù)結合GRACE重力衛(wèi)星數(shù)據(jù)研究了青藏高原東緣地殼流速度。文獻[44]利用GPS結合GRACE監(jiān)測了喜馬拉雅地區(qū)地殼水平形變。文獻[45]利用GRACE研究了青藏高原地殼垂直形變，并與GPS結果進行了對比。

本文利用1992—2020年歷次珠峰高程測量GNSS復測點監(jiān)測數(shù)據(jù)，并結合自1999年以來，20多年的青藏高原GNSS觀測數(shù)據(jù)，對青藏高原及周邊地殼形變、地殼隆升進行了定量研究，分析了青藏高原及周邊地殼三維形變特征及其動力學機制。

1 GNSS數(shù)據(jù)處理

本文收集了1992—2020年我國歷次珠峰高程測量建立的區(qū)域地殼監(jiān)測網(wǎng)GNSS數(shù)據(jù)6個復測站GNSS數(shù)據(jù)，1999—2020年陸態(tài)網(wǎng)絡位于青藏高原及周邊30個連續(xù)觀測的CORS站、352個定期觀測的區(qū)域站，2009—2020年尼泊爾14個CORS站GNSS觀測數(shù)據(jù)。站點分布如圖1所示。本文利用GAMIT軟件對GNSS數(shù)據(jù)基線解算，解算方案主要參數(shù)設置見表1。

圖1 青藏高原及周邊GNSS站點分布

表1 基線解算設置

解算獲取各站單日解時間序列，單天解精度水平方向優(yōu)于3 mm，高程方向優(yōu)于5 mm。利用GLOBK軟件進行整體平差，獲取ITRF2014框架下GNSS站點三維線性速率。

2 GNSS監(jiān)測點垂直形變精化

為研究青藏高原地殼隆升，首先需要精確獲取地殼垂直運動速率。GNSS時序不僅反映線性的構造運動，還包含了非線性變化特征，這種非線性變化一般是由大氣、陸地水等地表質量負荷引起，屬于非構造負荷形變，特別是垂直方向時序主要特征是非構造負荷導致的周期性變化。GNSS觀測站單天解時序如圖2所示。有學者研究發(fā)現(xiàn)，GRACE能夠監(jiān)測一定空間尺度上的陸地水變化，利用GRACE扣除陸地水負荷形變，可以消除約64%的GNSS共模誤差[46-47]。本文利用GRACE時變重力場模型，分析非構造負荷對GNSS高程時序的影響，排除陸地水等非構造負荷影響，進一步精化站點垂直運動速率。

圖2 GNSS時序

本文所采用的GRACE數(shù)據(jù)產品為美國德克薩斯大學空間研究中心(CSR)發(fā)布的GRACE Level-2(RL06)月重力場模型數(shù)據(jù)。主要計算過程為：利用衛(wèi)星激光測距(SLR)獲得的低階重力場球諧系數(shù)，替換時變重力場模型(GSM)數(shù)據(jù)的C20項及一階項；利用Fan組合濾波P3M15處理，去除噪聲和泄露誤差等影響；利用前后月份數(shù)據(jù)內插，補足時變重力場模型中缺失月份數(shù)據(jù)；以所有月份的數(shù)據(jù)均值為基準，從月重力場模型中扣除基準；根據(jù)扣除基準后的重力場球諧系數(shù)文件反演獲得陸地水儲量變化；根據(jù)負荷形變理論和相關計算公式，利用陸地水儲量變化，計算獲得其產生的垂直形變。

圖3展示了GARCE獲取地殼垂直運動時序與GNSS高垂直方向時序，可以看出兩者在周期、振幅上基本一致，這表明GNSS年周期主要是由非構造負荷形變引起的。

圖3 GNSS垂直時序GRACE非構造負荷形變改正

GRACE對于水平方向線性運動影響很小，在研究水平方向構造運動時可以忽略不計。但是對于垂直方向影響很大，特別是對于大部分定期觀測的區(qū)域站，由于數(shù)據(jù)量較少，且每期觀測時間并不完全固定，利用GRACE對GNSS垂直方向時序進行改正，可以有效提升線性速率提取的準確度和精度。

圖4展示了利用GRACE對垂直方向時序進行改正前后，GNSS垂直方向線性速率及其精度，可以看出，改正后GNSS站點垂直方向線性速率精度明顯提升，平均提升約40%，其中對于定期觀測的區(qū)域站提升更為明顯，平均達到約49%。為從不同角度研究青藏高原地殼運動特征，本文還將ITRF2014框架下速率轉換為相對歐亞板塊速率。站點三維運動速率如圖5所示。

圖4 改正前后GNSS站點垂直運動速率和精度

圖5(a)包括了ITRF2014框架下各個站點水平方向運動速率(藍色箭頭)和相對歐亞板塊的水平方向運動速率(紅色箭頭)。ITRF框架下站點水平速率顯示出明顯的分區(qū)塊特征，自南向北運動速率逐漸遞減，東向運動速率逐漸增加。

圖5 GNSS觀測站三維速率

而在扣除歐亞板塊整體運動趨勢后，水平方向運動量級明顯減小，運動方向分化更為明顯，東部右旋更加明顯，西部則出現(xiàn)左旋趨勢。

3 基于塊體模型的地殼形變分析

3.1 塊體模型構建

GNSS站點分布較為分散，且分布不均勻。為研究區(qū)域地殼整體運動特征，需要通過建立模型，實現(xiàn)格網(wǎng)化。青藏高原地殼運動最顯著的特征之一是具有明顯的分區(qū)塊特征[48-49]，因此本文采用三維塊體模型，通過構建地殼三維運動與形變模型，研究地殼三維形變特征[50-53]

Vs=Rt·Ω

(1)

(2)

(3)

式中，vn、ve、vu表示塊體任一點(λ,φ)南北向、東西向和垂直方向速率；r為地球平均半徑；(λo,φo)為塊體幾何中心坐標；Δx、Δy、Δz表示塊體x,y,z方向的整體位移；ωx、ωy、ωz是歐拉旋轉矢量；εe、εn、γen分別為塊體相對于塊體幾何中心的東向、北向主應變和剪應變。

塊體劃分采用參考我國CPM-CGCS2000板塊模型及地質研究成果[53-55]。本文主要研究范圍主要包括拉薩塊體、羌塘塊體、巴顏喀拉塊體和柴達木塊體。

本文分別利用ITRF2014框架下速度場和相對歐亞板塊速度場構建了塊體模型，利用塊體模型擬合站點三維速度場(紅色箭頭)，與實測GNSS速度場(藍色箭頭)進行對比，評估模型精度。結果如圖5、圖6所示，精度統(tǒng)計結果見表2、表3。

圖6 實測GNSS站點水平速率與模型擬合速率

表2 ITRF 2014框架下塊體模型精度

表3 相對歐亞板塊運動框架下塊體模型精度

圖6顯示，模型擬合水平方向速度場，與站點實測速率趨勢、量級基本一致。表2和表3統(tǒng)計了擬合速度場與實測速度場之間差值的均值和均方差，均值反映模型無偏性，均方差則反映模型有效性?？梢钥闯觯鱾€塊體差值N、E、U 3個方向均值接近零，說明本文構建的模型具有很好的無偏性；水平方向均方差最大約6 mm/a，最小僅有1 mm/a，顯示模型在水平方向具有良好的有效性。模型在巴彥喀拉塊體、羌塘塊體在南北向均方差較大的主要原因是塊體運動出現(xiàn)明顯的右旋趨勢，導致南北向速度場發(fā)生較大變化，反映到模型上，表現(xiàn)為站點南北向速率離散度較大。與其他學者采用歐拉矢量法[55]、有限元法[56]構建的速度場模型對比，本文構建的模型在無偏性更好，有效性方面在東西方向均方差優(yōu)于明顯優(yōu)于其他結果，南北方向均方差則基本一致。這表明，青藏高原地殼采用非連續(xù)活動塊體模型，考慮塊體內部應變構建速度場模型，更能體現(xiàn)地殼實際運動與應變特征。

圖7顯示，模型擬合垂直方向速率與站點實測速率，部分站點差別較大，對模型在垂直方向有效性有一定影響，但表2、表3顯示，模型垂直方向均值接近于0，表明模型無偏性較好，能夠反映地殼整體垂直運動趨勢。

圖7 實測GNSS站點垂直速率與模型擬合速率

3.2 基于塊體模型的地殼水平形變特征

塊體模型可以反映地殼三維運動與應變特征。圖8展示了模型構建的地殼水平形變場，包括格網(wǎng)速度場、塊體間相對運動。圖8藍色箭頭表示地殼水平速度場，紅色箭頭表示塊體間相對運動，底圖為水平運動角度。

圖8 基于塊體模型的地殼水平運動

圖8(a)展示了ITRF框架下塊體之間的相對運動。新生喜馬拉雅塊體和拉薩塊體之間西段主要為東西向相對運動，中段則具有明顯的擠壓狀態(tài)，東段又轉為東西向滑動；拉薩塊體與羌塘塊體之間同樣表現(xiàn)為明顯的擠壓，尤其是東段量級較大；羌塘塊體和巴顏喀拉塊體之間則主要表現(xiàn)為東西向相對滑動；巴顏喀拉塊體和柴達木塊體之間主要呈現(xiàn)擠壓運動。

圖8(b)顯示了去除歐亞板塊整體運動趨勢后，青藏高原地殼水平運動空間分布特征，總體趨勢基本一致，相對運動速率量級有所減小。

本文通過截取ITRF框架下東經(jīng)90°、北緯33°剖面展示不同區(qū)域地殼運動量級、方向分布(圖9)，喜馬拉雅山脈南部運動量級約為50 mm/a，方向大概在NE40°，向北量級有所減弱，拉薩板塊和羌塘板塊運動為40～45 mm/a，運動方向則有明顯變化，由NE30°逐漸旋轉為大概在NE0°，至巴顏喀拉塊體、柴達木塊體則北向速率略有減小，東向速率基本穩(wěn)定，量級在30～40 mm/a，方向在NE10～20°。

圖10展示相對歐亞板塊地殼水平運動空間變化特征，喜馬拉雅山脈南部運動量級約為40 mm/a，方向達到NE70°，向北量級逐漸減弱，拉薩板塊和羌塘板塊運動為20～30 mm/a，運動方向則有明顯變化，逐漸旋轉為大概在NE30°，往北至巴顏喀拉塊體、柴達木塊體則北向速率略有減小，東向速率基本穩(wěn)定，量級約在20 mm/a左右，方向則明顯增大。

圖9、圖10中，北緯33°剖面圖顯示在90°～98°之間羌塘塊體有一個明顯速率下滑的過程，結合圖8可以看出，該區(qū)域正是地殼沿玉樹-鮮水河斷裂逐漸右旋轉向的過程，2010年玉樹Ms7.1級地震就發(fā)生在這個區(qū)域。

圖9 ITRF2014框架下地殼水平方向運動量級與角度

圖10 相對歐亞板塊框架下地殼水平方向運動量級與角度

不同框架下GNSS運動速率存在一定差異，而基于速度場的地殼應變張量則與框架選取無關，可以從不同方面反映地殼形變特征。圖11展示了面膨脹率、最大剪應變率和最大最小主應變。

圖11(a)面膨脹率顯示，在印度板塊與歐亞板塊交界的喜馬拉雅造山區(qū)域，存在明顯的擠壓變形，特別是在2015年尼泊爾地震所在地區(qū)到加德滿都周邊，擠壓變形尤其明顯。而在青藏高原中部的拉薩、羌塘、巴顏喀拉中部則面膨脹率不明顯，平均約4.0±1.0nonstrain/a。

圖11(b)顯示，在剪應變率和最大最小主應變方面，加德滿都周邊形成一個明顯的應變集中區(qū)域。在拉薩塊體，主應變表現(xiàn)為南北向壓縮與東西向拉伸，特別是在加德滿都-日喀則方向是最大主應變方向，該方向面膨脹率為負值，表明該方向是拉薩塊體主要擠壓變形區(qū)域。

圖11 地殼應變率場

3.3 基于塊體模型的地殼垂直形變特征

圖12展示了模型構建的地殼垂直形變場。不同框架下，垂直形變基本一致，相對歐亞板塊框架下，模型垂直方向形變速率略低于ITRF2014框架下，而分布一致，包括塊體之間垂直方向相對運動也基本一致。因此，在研究垂直形變時，基本可以不考慮框架影響。

圖12顯示，青藏高原總體隆升趨勢明顯。位于印度板塊部分，垂直方向速率為負值，向北垂直方向逐步變?yōu)檎?，地殼逐漸呈隆升特征，拉薩、羌塘塊體整體隆升速率為2～3 mm/a，巴顏喀拉塊體隆升不明顯，平均速率低于1 mm/a，柴達木塊體同樣隆升不明顯，而祁連塊體則出現(xiàn)明顯上升的趨勢。

圖12 基于塊體模型地殼垂直形變

為研究青藏高原地殼整體形變特征，本文通過計算垂直速率與水平速率比值，研究了地殼運動傾斜度(圖13)。傾斜度可以在一定程度上反映地殼在運動過程中出現(xiàn)的彎曲程度。由圖13可以看出，在青藏高原中部拉薩塊體和羌塘塊體，特別是在兩個塊體的西部和中部，傾斜度明顯高于周邊地區(qū)，表明在該區(qū)域地殼運動中出現(xiàn)一定程度向上彎曲的變形特征。

圖13 地殼運動傾斜度

本文截取了90°E、95°E，31°N、34°N等剖面，展示不同剖面地殼垂直運動速率分布，以及傾斜度等(圖14)。

圖14 不同經(jīng)緯度剖面地殼運動垂直速率與傾斜度

從各個經(jīng)度剖面可以看出，青藏高原在中部都有一個明顯的隆起部分，可能的解釋是青藏高原受南北壓縮出現(xiàn)中部隆升大于兩端的向上彎曲現(xiàn)象。從緯度剖面可以看出，青藏高原西部主要表現(xiàn)為整體隆升，而在東部則向東地殼隆升逐漸下降。

4 結論與討論

目前，包括CORS站和定期觀測站在內的GNSS監(jiān)測網(wǎng)已覆蓋了青藏及其周邊地區(qū)，為精確定量研究青藏高原地殼形變、隆升及其動力學機制提供了豐富、可靠的數(shù)據(jù)資料。

本文基于GNSS長期監(jiān)測數(shù)據(jù)，利用塊體模型從多個方面分析了青藏高原地殼三維形變特征，探討了地殼形變與隆升機制。

顧及塊體內部應變的活動塊體模型，模型擬合結果與實際GNSS觀測結果差值均值接近于0，具有優(yōu)良的無偏性，其有效性方法水平和垂直方向均方差均優(yōu)于1 cm/a，可以在大區(qū)域內精確反映地殼實際運動特征。

GNSS水平速度場模型，特別是相對歐亞板塊速度場表明，青藏高原受到印度板塊南東方向擠壓、塔里木塊體北東方向擠壓，從而出現(xiàn)西部邊緣西向旋轉擠出、中部和東部東向旋轉的整體趨勢。

利用GRACE精化GNSS垂直速度場，可以更精確地獲取地殼垂直方向構造運動趨勢，其中喜馬拉雅山南北兩側臨近區(qū)域，GNSS垂直速率改正后略有降低，而在羌塘塊體、巴彥喀拉塊體、柴達木塊體等區(qū)域，改正后垂直速率有所上升，改正后垂直速率精度整體提升約40%，大部分GNSS站點垂直方向長期運動趨勢均為隆升。

基于塊體模型的地殼三維形變場顯示，青藏高原南部拉薩塊體，受到南北向擠壓，南北向擠壓收縮是塊體主要應變特征，特別是中段具有明顯的應變集中現(xiàn)象，2015年尼泊爾Ms8.1級地震就發(fā)生在該區(qū)域。同時南北向擠壓也導致地殼隆升相對較大。

位于青藏高原中部的羌塘塊體，受到拉薩塊體和塔里木塊體擠壓，其水平方向自西向東沿黎嘉斷裂右旋運動，東部則受到川滇塊體阻擋。塊體西部存在一定南北向擠壓應變，中部和東部則沒有明顯的南北向擠壓應變，地殼表面積略有膨脹；高程方向，羌塘塊體地殼隆升總體大于周邊塊體，平均約2.0 mm/a，其垂直速率與水平速率比率大于周邊地區(qū)，表現(xiàn)為青藏高原地殼運動過程中，在中部地區(qū)出現(xiàn)向上突出隆起的垂直形變特征。

巴顏喀拉塊體整體東向旋轉，其東部龍門山地區(qū)擠壓作用強烈，為強烈擠壓應變區(qū)；西部邊界以伸展應變?yōu)橹?，北邊界主要表現(xiàn)為東向滑移。水平運動整體具有順時針旋轉特征。地殼隆升，平均約為0.8 mm/a。

柴達木塊體受塔里木塊體東向擠壓，整體以東向運動為主，同時又收到祁連塊體北向擠壓，導致塊體整體略有南北向略有收縮，地殼隆升速率較巴彥喀拉塊體略大，約1 mm/a，但低于祁連塊體的約1.8 mm/a。

研究表明，青藏及周邊地殼形變具有明顯的分區(qū)特征。南北向收縮主要表現(xiàn)在拉薩塊體，東西向伸展則主要表現(xiàn)在巴顏喀拉塊體，中部羌塘塊體沒有明顯南北向水平擠壓應變，有長期東西向伸展趨勢。垂直方向總體隆升，特別是在中部羌塘塊體，地殼隆升速率最大、水平面上面積有擴大趨勢，反映出青藏高原地殼有增厚的可能，這在一定程度上支持了地殼增厚學說。

致謝：感謝自然資源部第一大地測量隊提供的珠峰高程測量區(qū)域監(jiān)測網(wǎng)GNSS數(shù)據(jù)，感謝“中國地球運動觀測網(wǎng)絡”“中國大陸構造環(huán)境監(jiān)測網(wǎng)絡”提供了青藏高原GNSS觀測數(shù)據(jù)，感謝UNAVCO提供了尼泊爾CORS站觀測數(shù)據(jù)，感謝GAMIT/GLOBK軟件和GMT軟件團隊。