方智偉， 鄒蓉， 李志才， 王敏， 譚凱， 楊少敏， 王琪

1 中國地震局地震研究所， 地震大地測量重點實驗室， 武漢 430071 2 中國地質大學地球物理與空間信息學院， 武漢 430074 3 國家基礎地理信息中心， 北京 100044 4 中國地震局地質研究所， 北京 100029 5 武漢市應急管理局， 武漢 430010

0 引言

青藏高原隆升是新生代大陸活動構造最引人注目的地質事件(Molnar and Tapponnier，1975)，有關此次陸內造山的起始、過程、機制及力源一直是大陸動力學研究的熱點問題(Tapponnier et al.，2001).青藏高原隆升研究涉及對其歷史高程的推演及不同階段隆升速率的厘定，目前前者主要借助高原內第三紀地層湖相沉積物的氧同位素分析(Rowley and Currie，2006)，后者包括從中新世以來出露變質巖的低溫年代記錄推算中下地殼折返率(Herman et al.，2010)，或利用山前第四紀沖積盆地測定抬升速率(Lavé and Avouac，2000)，以及利用大地測量(精密水準、GPS、InSAR)測定現(xiàn)今的地表垂直運動速率等(張青松等，1991； Jackson and Bilham，1994；Xu et al.，2000；張全德等，2001；姜衛(wèi)平等，2008；Grandin et al.，2012)，觀測手段雖多，但迄今獲取的可靠資料卻嚴重不足，限制了青藏高原構造演化研究的發(fā)展.

與地質學手段相比，用大地測量監(jiān)測垂向變形具有量化清晰、物理意義明確、受客觀因素制約較少的特點，對判斷青藏高原隆升狀態(tài)、動力作用具有指示意義.近期的連續(xù)GPS觀測揭示了青藏部分地區(qū)的三維形變場，初顯用大地測量方法解決大陸構造演化中動力學問題的良好發(fā)展勢頭(Sun et al.，2009；Fu and Freymueller，2012；Ader et al.，2012；Liang et al.，2013).盡管如此，現(xiàn)有觀測結果仍有較大的不一致，例如早期跨青藏高原精密水準網復測資料顯示其最大隆升位于藏南一帶，隆升幅度達到10 mm·a-1(張青松等，1991)，這與尼泊爾境內水準資料推算的結果并不一致(Jackson and Bilham，1994).早期利用4年的GPS資料推測的高原中部唐古拉山一帶隆升速率高達16 mm·a-1(Xu et al., 2000)，而用中國大陸構造環(huán)境監(jiān)測網絡(陸態(tài)網絡，下同)GPS區(qū)域站10年資料解算的高原中部的構造隆升不過2 mm·a-1(Liang et al.，2013)，差異十分明顯.

相對于精密水準和InSAR而言，GPS觀測具有基準統(tǒng)一、精度均勻、觀測簡便，受地形、地貌限制較小等優(yōu)勢，但GPS站點的高程(24 h)觀測精度不高于10 mm，而大陸內部地殼升降速率低于10 mm·a-1，垂向位移監(jiān)測的信噪比低，定期、流動的GPS監(jiān)測往往難以勝任.青藏地區(qū)連續(xù)GPS站點的增加，使得短時間內高精度測定垂直形變場逐漸成為可能(Sun et al.，2009；Ader et al.，2012；Hao et al.，2016).本文基于青藏高原南緣GPS連續(xù)站1995—2020年間觀測數據，系統(tǒng)解算測站坐標時間序列(時序，下同)，求取站點垂直運動速率，合成具有信噪比、點密度較高的區(qū)域垂向位移場，為分析區(qū)域構造變形模式提供觀測約束.

1 資料

本文研究匯集了印度次大陸、喜馬拉雅及青藏南部地區(qū)固定觀測98個GPS站點資料(圖1)，這些站點大多設立在穩(wěn)定巖石上，采用扼徑圈天線，觀測持續(xù)、質量穩(wěn)定，可在相對較短的(5年)時間內，以優(yōu)于1 mm·a-1的精度監(jiān)測垂向位移(Lau et al.，2020).而區(qū)域內不定期觀測的流動站的數據質量參差不齊，且觀測10年左右的陸態(tài)網絡區(qū)域站，其垂向位移觀測精度一般在1～3 mm·a-1(Liang et al.，2013)，本文不予采用.

1.1 cGPS觀測

從1997年以來，美國科羅拉多大學、加州理工學院與尼泊爾有關機構合作，布設基本覆蓋尼泊爾全境的GPS連續(xù)臺網(表1)，該網觀測情況復雜，部分站點已不復存在或不再觀測，有些則經歷了2004年蘇門答臘地震和2015年尼泊爾地震，維持至今.加州理工學院分析了該網早期站點的三維運動狀況(Ader et al.，2012)，國內也有相關研究(Liang et al.，2013；Pan et al.，2018).本文系統(tǒng)處理該網從早期建站到今年春季數據，包括美方在加德滿都等地為監(jiān)測尼泊爾地震震后形變增補的新站(Zhao et al.，2017).

圖1 cGPS站點分布及觀測時長Fig.1 Distribution and observation time span of cGPS stations

表1 本文所用cGPS站觀測情況Table 1 Observations of continuous GPS stations used in this paper

續(xù)表1

本文將周邊國際GNSS服務組織(International GNSS Services，IGS)站點觀測數據納入處理，包括位于印度北方邦的LCK1-LCK4、印度中部的HYDE和IISC、不丹王國的TIMP、RBIT及位于拉薩的LHAZ(與陸態(tài)網絡站LHAS并址)，這些測站中最長持續(xù)觀測了25年.陸態(tài)網絡12個連續(xù)觀測基準站分布于青藏高原南部地區(qū)，此前一些研究多集中這些站點水平運動狀況(Liang et al.，2013；Wang and Shen，2020)，與此前研究比，本文增加了最近幾年觀測資料.為增加西藏境內位移場的空間分辨率，本文還收集中國氣象局于藏南布設的3個準連續(xù)站觀測資料(表1)，這些站點觀測時長約5年，尚未用于垂向位移的研究(Wang and Shen，2020).此外，本文引用了中國科學院青藏高原研究所于藏南布設的4個準連續(xù)站2007—2012年觀測的處理成果(Liang et al.，2013)，以及印度國家地球物理研究所在喜馬拉雅西段22個連續(xù)GPS站2009—2015年觀測的處理成果(Yadav，2017；Yadav et al.，2019).

1.2 坐標解算

本文大部分GPS測站的原始數據處理基于噴氣推進實驗室(Jet Propulsion Laboratory，JPL)研制的GIPSY/OASIS Ⅱ-5.0軟件.GIPSY處理普遍采用精密單點定位(Precise Point Positioning, PPP)模式(Zumberge et al.，1997)，獲得ITRF2014框架下各站點的單日(UTC：0～24 h)時段坐標(Altamimi et al.，2016).本文具體處理策略包括：載波相位和偽距數據來自RINEX格式原始觀測文件(采樣率30 s)，L1和L2雙頻載波相位經過基于精碼偽距輔助的周跳修復后，轉換為無電離層折射L3組合相位(其截止高度角為10°)，并用IGS絕對天線相位中心模型修正測站天線中心相位偏差.對L3相位觀測值的處理采用卡爾曼濾波/平滑算法，同步估算測站坐標、衛(wèi)星單差的組合相位模糊度，以及作為隨機游走變量的測站天頂方向大氣延遲濕分量.受軟件版本限制，沒有借用JPL提供的衛(wèi)星初始相位改正文件以輔助固定衛(wèi)星單差組合相位模糊度.衛(wèi)星軌道固定為JPL發(fā)布的無基準事后精密星歷，同時采用JPL提供的衛(wèi)星鐘差文件.解算中各類物理模型分別是固體潮和極潮改正的IERS2010模型，以地球質心為基準的TPXO7.0(CM)海潮模型及對流層折射改正的全球經驗大氣+斜向延遲映射函數(GMF/GPT)組合模型.最后利用JPL發(fā)布的單日坐標轉換參數，將無基準的單日坐標統(tǒng)一轉換至ITRF2014參考框架下，形成各個站點可供變形分析的時序(圖2).本文氣象局測站(也包括喜馬拉雅西段的cGPS站)采用GAMIT軟件處理原始數據，基于雙差分定位模式解算出包含這些測站及周邊測站的單日區(qū)域網解，然后用GLOBK軟件將區(qū)域網與全球IGS網聯(lián)網平差，求取這些測站在ITRF2008下坐標時序，相關信息詳見有關文獻(Zhao et al.，2017；Wang and Shen，2020).

圖2 GPS單日解垂向分量的時間序列圖中虛線為2004年和2015年發(fā)震時刻，綠色代表經歷2015年尼泊爾地震的臺站.淺藍色代表震前臺站，深藍色為同一站點的GAMIT時序，紅色為震后布設的臺站.Fig.2 Time series of the vertical component of GPS daily solutionsThe dotted lines in the figure show the time of earthquakes in 2004 and 2015, and the green represent stations that cross the 2015 Nepal earthquake. The light blue represents the stations deployed before the earthquake, the dark blue represents the time series of the same stations obtained by GAMIT, and the red represents the stations deployed after the earthquake.

1.3 垂直坐標時序分析

本文假定GPS垂向時序中包含三種可分離的信號成分：1)與構造運動或與晚更新世冰后回彈有關的線性位移信號；2)與測站周邊環(huán)境有關(如大氣、水文負荷的四季輪替)，可簡化為具有一整年、半年周期變動特征的三角函數信號；3)大地震引起的同震位移和持續(xù)震后位移，前者可用單邊函數表示，后者以指數或對數函數表示.綜合以上，本文垂向時序展開為以下方程式：

(1)

式中：h(t)為測站大地高觀測值，t為觀測歷元，t0為首次觀測時刻，ci為季節(jié)性波動幅度，φi為季節(jié)峰值相位，H(t)為單邊函數，dj為同震位移，tj為發(fā)震時刻，ej為震后變形幅度，τj為震后變形的半衰期.本文以式(1)擬合實測數據，解算時依據每個數據的解算精度定權，建立相應的誤差模型，以加權最小二乘法估算各運動參數，準確反映測站垂直速率的信噪比.本文只關注2004年蘇門答臘MW9.2、2015年尼泊爾MW7.8兩次大震對時序的擾動，忽略其他地震的影響.另外，本文季節(jié)性波動項只考慮周年和半年周期項.

1.4 誤差模型

GIPSY給出的單日大地高程的解算精度(名義誤差)約為3～4 mm，主要反映了觀測當日載波相位中的隨機噪聲.其他因素，如星歷、電離層、對流層折射、測站周邊溫度變化、多路徑等因素引起的各種偏差不在其中，這些偏差對于單日坐標而言屬于系統(tǒng)性的，相鄰數日或數周內系統(tǒng)偏差變化不大，對站點坐標的擾動一致，在時序中表現(xiàn)為有色噪聲(不同時間觀測值的統(tǒng)計性質相關).如簡單將實測數據誤差視為白噪聲，雖不影響對各個信號的提取，但大大高估其解算精度，故本文時序分析中的誤差模型采用白噪聲+有色噪聲的組合方式(Bos et al.，2013)，并結合赤池信息準則和貝葉斯準則選擇其中最合適的模型.

借鑒此前尼泊爾地區(qū)GPS數據分析現(xiàn)狀(Ader et al.，2012)，本文的有色噪聲模型為廣義高斯-馬爾科夫模型，白噪聲、冪率噪聲、閃爍噪聲和隨機游走噪聲模型都是其特例(Bos et al.，2013)，因此理論上廣義高斯-馬爾科夫模型能更好反映有色噪聲的實際狀況.每個測站具體的有色噪聲模型主要依據其擬合殘差的功率譜指數而定，相關計算采用Hector軟件，同時給出站點垂向速率、周期項參數及其不確定性(表2).

2 垂直運動速率及精度評估

表2和圖2所展示各類參數及信息可簡單歸納如下：1)經GIPSY/OASIS Ⅱ-5.0處理的各測站時序的標準差(RMS)在6～10 mm之間，時序擬合殘差限于[-20 mm, +20 mm]條帶內.內華達大學基于高版本的GIPSY X-1.0處理的同為ITRF2014框架下的模糊度固定解時序(http：∥www.geodesy.unr.edu)的標準差為3～6 mm，比本文模糊度浮點解略好，但對比本文與內華達大學計算的IISC、HYDE、LHAZ、CHLM的長期速率, 互差小于0.2 mm·a-1；2)經GAMIT處理的ITRF2008框架下時序標準差在3 mm左右，雖然雙差精密定位較單點精密定位能更好消除一些共模誤差，能更好分辨出微弱變形信號，但在提取長期速率上，GIPSY與GAMIT時序的一致性保持在0.05 mm·a-1以內；3)區(qū)域內升降速率約束在[-2 mm·a-1, +6 mm·a-1]區(qū)間(相對ITRF框架原點——地球質心)，估算精度多半優(yōu)于0.5 mm·a-1；4)各測站的季節(jié)性波動信號與地震位移信號得到較好的分離，殘差時序沒有明顯的扭曲和跳變；5)除少數幾個測站外，各站的周年項波動幅度一般在6～10 mm，估值誤差在2 mm以內，周年項峰值位于5～6月的印度季風期，喜馬拉雅山前測站波動(最大15 mm)明顯大于高原內部測站；6)半年項波動幅度不超過周年項的一半，但估值誤差類似；7)尼泊爾地震導致震區(qū)內測站(最大1.3 m)永久性垂向位移，震后變形也導致部分測站短期升降10～60 mm，此外蘇門答臘地震導致尼泊爾及印度次大陸測站隆升大致10～15 mm.

表2 cGPS時序分析結果Table 2 Time series analysis results of continuous GPS stations

續(xù)表2

2.1 垂直速率的區(qū)域特征

GPS測站的垂直運動特征概括為：恒河平原及喜馬拉雅山前(海拔500 m以下)總體下降不大于3 mm·a-1，DNGD沉降幅度最大，速率為3.1±0.4 mm·a-1；低喜馬拉雅地區(qū)(海拔500～1500 m)處于升降轉換區(qū)，高喜馬拉雅地區(qū)(海拔1500～4500 m)相對印度次大陸最大隆升接近2～8 mm·a-1，在30 km距離內，抬升速率從1 mm·a-1(NAGA)急劇增加到6 mm·a-1(GUMB)，與地形陡變完全一致. 藏南地區(qū)(海拔4500 m以上)的GPS站在大約50 km距離內抬升速率從4～6 mm·a-1(GHER、CHLM)減少到1～2 mm·a-1(LMJM、SYBC)，雅魯藏布縫合線以北(平均海拔4500 m)GPS站逐漸從上升約1 mm·a-1轉變至下降1 mm·a-1，XZNM(尼瑪)沉降最大為1.5±0.4 mm·a-1.

喜馬拉雅1975—1990年間兩期精密水準數據顯示高喜馬拉雅與低喜馬拉雅相對前陸盆地的抬升分別約為4±1 mm·a-1和2.0±0.3 mm·a-1(Herman et al.，2010).尼泊爾西部的InSAR觀測更凸顯這種升降變化(圖3)，高喜馬拉雅相對于恒河平原最大隆升達到8 mm·a-1(Grandin et al.，2012).盡管本文直接解算的尼泊爾cGPS垂向速率場及印度西北地區(qū)cGPS觀測結果(Yadav et al.，2019)的空間分辨率較低，但結合InSAR、精密水準，現(xiàn)有大地測量綜合結果完整反映了大陸俯沖帶前陸盆地沉降、山前快速抬升、造山帶后緣隆升延展、幅度急劇降低這一典型特征.

圖3 喜馬拉雅地區(qū)垂向速率剖線左圖：InSAR圖像、水準線路和GPS站分布圖，紫色和藍色方框代表A、B剖線，色棒標示InSAR和水準點速率，黑色箭頭為相對歐亞板塊的GPS水平運動速率(Wang and Shen，2020)；右圖：GPS(紅或紫色方框)，InSAR(棕色圓點)和水準(藍色圓點)速率及地形合成圖.GPS速率未做任何修正.圖中近場和遠場GPS采用兩種距離尺度表示，水準和InSAR誤差棒代表1 mm·a-1的假設誤差，GPS誤差為真實誤差.綠色虛線為震間運動模型的預測速率(Grandin et al.，2012)，棕色圓環(huán)代表此前的GPS結果(Fu and Freymueller，2012).Fig.3 Section lines of vertical velocities in HimalayaLeft: InSAR image, leveling line and distribution of GPS stations. Blue and purple boxes represent section lines A and B. Color bars indicate InSAR and leveling benchmark velocity. Black arrows represent GPS horizontal motion rates relative to the Eurasian plate (Wang and Shen, 2020). Right: Composite image of GPS (red or purple boxes), InSAR (brown dots) and leveling (blue dots) rate with terrain. No correction has been made to the GPS rate. In the figure, near-field GPS and far-field GPS are represented by two distance scales. Leveling and InSAR error bars represent hypothetical errors with 1 mm·a-1, while the GPS error is real error. The green dotted line represents the predicted velocity of the interseismic motion model (Grandin et al., 2012), and the brown rings represent the previous GPS results (Fu and Freymueller, 2012).

2.2 內符合精度分析

此前尼泊爾境內站點及藏南部分GPS站垂直速率曾有文獻記載(Fu and Freymueller，2012；Ader et al.，2012；Liang et al.，2013；Pan et al.，2018)，本文將處理結果與此前結果逐一對比(圖4).加州理工最早公布測站速率誤差一般為1～3 mm·a-1(Ader et al.，2012)，比本文大3倍.加州理工除顧及有色噪聲外，還考慮了時序不連續(xù)可能引入的階躍誤差，因此總體上精度估計偏保守.其他三組與本文采用完全一致的誤差評估算法(Fu and Freymueller，2012；Liang et al.，2013；Pan et al.，2018)，圖4顯示本文速率精度總體上最優(yōu)，這得益于某些測站已有了更長時序，當時間跨度大于5年時(XZAR、XZNM、SNDL)，PPP給出的垂向速率精度在0.5～0.8 mm·a-1，大于10年時(KKN4、KAWA、CHLM)，誤差進一步降低到0.3 mm·a-1，而超過20年觀測(LHAS、HYDE、IISC)，誤差甚至可低至0.1 mm·a-1.與此對照，陸態(tài)網絡觀測時長12年的區(qū)域站精度一般低于1 mm·a-1(Liang et al.，2013)，如此估計25年觀測歷史的流動站有望將誤差降低到0.5 mm·a-1以內.

圖4 4組GPS速率結果比對與此前公布的4組速率結果對比，西藏地區(qū)陸態(tài)站點用綠色點位表示.虛線陰影區(qū)代表互差在+/-1 mm·a-1以內區(qū)域.Fig.4 Comparison of GPS velocity results from four groupsComparison with the previously published four groups of velocity results, the CMONOC sites in Tibet are represented by green dots. The dotted shaded areas represent the areas with errors within +/-1 mm·a-1.

需要解釋的是：本文中少數站點盡管具有更長時序，但有些站點的速率精度仍低于此前結果(Pan et al.，2018).由于利用了尼泊爾震后數據，為抵消同震和震后變形影響，本文引入了附加參數，相應降低了長期速率的估計精度.此外，GHER、LMJG震后位移量較小，從時序中分離震后、震間信號需要借助一定先驗假定，比如，震后位移與同震位移方向一致，震后位移幅度不超過同震位移等，因此速率計算與先驗信息有關，有一定的不確定性.即便如此，本文與此前結果沒有系統(tǒng)性偏差，例如與最近一次的結果(Pan et al.，2018)比，速率互差在1 mm·a-1內，僅PYUT、DRCL、SNDL三站互差超過2 mm·a-1.

2.3 模型對比及外符合精度

喜馬拉雅變形模型以往主要基于GPS水平位移、精密水準、InSAR三種資料(Bilham et al.，1997；Grandin et al.，2012)，GPS垂直位移數據因精度較差，對建模貢獻甚小(Fu and Freymueller，2012).為客觀評估本文結果，將GPS測站分東西二組(圖3)，沿N17°E方向構建兩條跨喜馬拉雅、寬度近200 km的帶狀速率剖線(剖線A、B)，剖線A內GPS速率與C波段ENVISAT衛(wèi)星InSAR視線方向速率對比(Grandin et al.，2012)，剖線B則是與尼泊爾一等精密水準復測結果對比(Jackson and Bilham，1994).

跨喜馬拉雅的InSAR干涉圖像覆蓋83°～84°經度帶，從低喜馬拉雅延伸到雅魯藏布，南北長250 km(圖3)，是處理2003—2010年間29幅降軌原始圖像，將不同時段、獨立干涉圖疊加而成(Grandin et al.，2012).變形起算點位于圖像最北端，變形信號從北到南累積，南北兩端像素點速率差接近10 mm·a-1，InSAR條帶內各像素點變形精度估計為0.7～3.1 mm·a-1(Grandin et al.，2012)，變形分辨率為每公里0.2 mm·a-1.盡管理論上InSAR圖像代表地表變形在衛(wèi)星視線方向(LOS)投影，但實際上更多地反映了垂向變動的分布，就A剖線而言，LOS 速率與垂向速率的偏差不超過0.5 mm·a-1.本文將InSAR最北端像素點的速率設為JRGR速率值1.2 mm·a-1，推算整條剖線其他像素點在ITRF框架上的速率.

尼泊爾境內的精密水準線路南起低喜馬拉雅斯瓦里克山前(Jackson and Bilham，1994)，向北經加德滿都終于高喜馬拉雅山前的西藏樟木口岸，總長350 km，沿線高程起伏不超過1500 m, 該線路1977—1991年間觀測，觀測誤差為1.1 mm·km-1，按誤差傳播率估計，水準線路南、北兩端的速率差相對精度為1.6～2.8 mm·a-1(Jackson and Bilham，1994).珠峰北坡一帶的水準線路與尼泊爾獨立，起算點位于定日(陳俊勇等，2001).本文將尼泊爾水準起始點速率等同于SIMP的速率值(-1.0±0.5 mm·a-1)，而將珠峰北水準線路的起算點設為最靠近XZZF(1.6±0.5 mm·a-1)的一個測點，相應推算各自路線上其他水準點的ITRF框架速率值.

對比GPS、精密水準、InSAR資料可見(圖3)，三者在1 mm·a-1的精度范圍內保持了一致，剖線B內只有DAMA、RMJT、SYBC偏差較大(1.5～2.0 mm·a-1)，如果將垂向變形的模型值(Grandin et al., 2012)與GPS實測值對比，則DAMA、TPLJ、SYBC偏差大于1.5 mm·a-1.偏差較大站點都遠離水準線路，可能疊加了局部變形效應，與水準點揭示的變形特征有所不同.

剖線A中InSAR與GPS的一致性更佳，只有KLDN與InSAR偏差大于4 mm·a-1.原因可能是山前一帶InSAR觀測誤差較大或KLDN疊加了前沿斷裂變形的影響.與模型值(Grandin et al., 2012)對比，只有KLDN、YARE、XZZB有近2 mm·a-1的偏差，其中YARE、XZZB與所在區(qū)域的InSAR觀測值都系統(tǒng)性偏離了模型值，且GPS、InSAR的變化趨勢也相同，這種異常變形可能與2006—2008年仲巴序列地震有關.整體而言，本文比早期結果(Fu and Freymueller，2012)更接近模型值.

2.4 垂直速率修正

GPS測定的垂向變動主要包含三種運動成分：1)板塊擠壓引起的構造變形；2)地表水變化導致的升降，近年來喜馬拉雅冰川加速融化和恒河平原地下水抽取加劇(Rodell et al.，2009；Yi and Sun，2014)，將導致區(qū)域地殼大幅調整；3)末次冰期以來北半球巖石圈松弛變形及青藏高原冰后回彈.依據近年來GPS和水準觀測，印度—歐亞碰撞帶內板塊水平移動30～40 mm·a-1(Wang et al.，2001)，升降不超過6 mm·a-1(Jackson and Bilham，1994；王慶良等，2008；Ader et al.，2012；Liang et al.，2013)，其垂向運動受環(huán)境影響更大.

地表水變化引起的均衡調整一般根據時變重力場推算的區(qū)域地表載荷變化及彈性地球變形響應來估算(Fu and Freymueller，2012).本文利用最新的GRACE觀測結果，假定時變重力場完全是地表水變化決定，忽略其他因素的貢獻(Rao and Sun，2021)，估算喜馬拉雅及藏南地區(qū)均衡隆升大致在0.2～0.4 mm·a-1范圍，比此前估算0.4～0.8 mm·a-1低一半左右(Pan et al.，2018).兩組結果皆有0.05 mm·a-1的標稱精度，大大低于兩組間的系統(tǒng)性偏差，說明用不同濾波算法平滑GRACE資料，給出的修正值存在至少0.2～0.3 mm·a-1的不確定性.

末次冰盛期導致的青藏高原隆升同樣存在不確定問題.按極端區(qū)域性冰川模型推測，青藏高原中部地殼冰后回彈的升幅可達4～7 mm·a-1(Kuhle，1998；Kaufmann，2005)，但另有研究表明青藏高原冰蓋規(guī)模不大、影響極小(Derbyshire，1991；Matsuo and Heki，2010；Liang et al.，2013)，本文也不加考慮.但北半球巖石圈松弛變形不容忽視，據現(xiàn)有大地測量約束下的全球冰川模型，亞洲大陸的冰后回彈為0.3～0.5 mm·a-1(Peltier，2004；Peltier et al.，2015)，為此本文統(tǒng)一用0.4 mm·a-1作為修正值.需要說明的是，現(xiàn)有修正值的誤差難以估算，本文取修正值的一半.

與以往研究不同(Liang et al.，2013；Hao et al.，2016；Pan et al.，2018)，本文還考慮地心運動對GPS速率的影響.根據全球大地測量結果(Argus，2012；Métivier et al.，2020)，末次冰期以來北半球巖石圈松弛變形以及內部質量遷移導致地球形心相對其質心以0.5～1.0 mm·a-1速率向北極運移，由于ITRF框架的原點為地球質心，GPS實測速率相對于質心，而構造運動相對于地球形心，地心的視運動將導致中低緯度地區(qū)真實的垂向運動偏大0.4 mm·a-1左右(Ding et al.，2019)，必須加以修正.本文按0.8 mm·a-1地心速率乘測站緯度正弦值計算修正值，修正值的精度大致在0.2 mm·a-1左右.綜合以上因素，現(xiàn)有GPS速率中非構造影響在1 mm·a-1左右，按保守估計，修正值的不確定度大致在0.3～0.4 mm·a-1.

3 討論

本文研究涉及印度—歐亞板塊碰撞三個不同的構造區(qū)(Molnar and Tapponnier，1975；Tapponnier et al.，2001)，分別為：剛性的印度板塊，整體性向北擠壓；位于板塊邊界的喜馬拉雅，造山楔強烈擠壓變形，向南仰沖、推覆在印度板塊之上；青藏高原內部，處于擠壓、拉張狀態(tài)的拉薩、羌塘地塊.印度板塊的地殼厚度不到40 km，藏南(喜馬拉雅北坡至雅魯藏布)地區(qū)厚度接近75 km，增厚近一倍,高原中部地殼厚度減薄到65 km(Nábělek et al.，2009)，地殼從增厚到減薄轉換對應了GPS垂向運動從下沉到隆升、再到下沉的空間展布.

3.1 印度次大陸

沿A剖線印度板塊相對歐亞板塊的水平速率從北到南為36～38 mm·a-1(Wang and Shen，2020)，如扣除非構造因素，現(xiàn)有GPS表明印度板塊前緣部分已下沉2～3 mm·a-1，使其以3°～4°傾角向下插入歐亞板塊.根據尼泊爾地震的同震破裂產狀，位于低喜馬拉雅底部的印度板塊以7°傾角向北運移(Wang and Fialko，2015)，遠震接收函數成像顯示大約在高喜馬拉雅到雅魯藏布之間，印度板塊的傾角陡變到11°～12°(Nábělek et al.，2009)，板塊邊界面的傾角每100 km就減小3°～4°(板塊曲率(0.5～0.6)×10-6m)，為海洋俯沖帶傾角的1/10～1/3(Bletery et al.，2016).考慮到恒河平原新生代沉積層最大深度為5～6 km(Decelles et al.，1998)，大陸俯沖的前陸伸展盆地至少100 km寬，從NEPJ延伸到BHUP.現(xiàn)有GPS結果還表明，德干高原北部為約500 km寬的前緣隆起區(qū)(從BHUP到JBPR)，德干高原南部(IISC、HYDE)幾乎沒有構造性垂向運動，GPS觀測與造山理論的預測一致(Molnar and Lyon-Caen，1988)，展示了印度板塊在青藏高原的重壓下，撓曲由淺入深的幾何結構(圖5).

圖5 青藏高原南緣垂向運動速度場GPS 實測速率經地心(坐標框架原點視)運動、冰川均衡調整模擬和地表水變化GRACE觀測三種改正，改正后速率見表1.Fig.5 GPS-derived vertical velocity field in southern Tibetan PlateauThe GPS-derived rates have been corrected by geocenter motion (frame origin apparent motion), glacial isostatic adjustment simulation and GRACE-derived groundwater variation. The corrected rates are shown in Table 1.

3.2 喜馬拉雅

如果GPS觀測到的印度次大陸前緣部分沉降可視為永久性、不可恢復變形，但GPS觀測的喜馬拉雅變形大部分是彈性的、與板塊邊界錯動及大震周期有關的變形.作為板塊邊界的喜馬拉雅主沖斷層最南緣的100 km段在兩次大震間完全閉鎖，低喜馬拉雅地區(qū)的震間變形緩慢，高喜馬拉雅地區(qū)快速隆升、變形累積(Bilham et al.，1997).這種隆升不是單調持續(xù)的，特大地震將導致大部分隆升回落，例如2015年尼泊爾地震時海拔3000 m以下的地帶全面隆升，而海拔3000 m以上地區(qū)整體下降(Elliott et al.，2016)，凸顯出地震在彈性變形從高海拔向低海拔區(qū)域的遷移過程中的關鍵作用.

GPS反映的彈性變形狀態(tài)與瞬時地震觀測的結果完全相反，加德滿都KKN4沉降速率為0.2 mm·a-1，尼泊爾地震中卻抬升1.3 m，中尼邊界地帶高海拔的CHLM隆升速率2.6 mm·a-1，地震時下降0.6 m，彈性變形集中在地形陡變帶，永久性變形集中在低海拔的山前地帶，喜馬拉雅地區(qū)現(xiàn)今變形表現(xiàn)為雙重性.喜馬拉雅的大震周期平均500～1000年(Feldl and Bilham，2006)，2015年尼泊爾地震不能一次性將大震間積累的全部變形從高喜馬拉雅轉移至斯瓦里克，只有1934年比哈爾MW8.2地震乃至更大地震才有可能錯斷整個主沖斷裂，至最南端主前緣斷裂(Sapkota et al.，2013)，將震間絕大部分變形遷移至前緣地帶，轉換為永久性地形抬升.但高喜馬拉雅地帶GPS隆升速率變化與地形變化高度一致，表明這種轉換并不徹底，仍有部分彈性變形永久存留在高海拔地帶，維持其陡變地形(Meade，2010).

3.3 青藏高原

最早的跨青藏高原的精密水準顯示，藏南(喜馬拉雅北坡至雅魯藏布)地區(qū)近10 mm·a-1隆升(張青松等，1991)，這種大幅度現(xiàn)代隆升得到早期流動GPS觀測的佐證(Jackson and Bilham，1994).但本文結果，該地區(qū)GPS站僅有0～2 mm·a-1的上升幅度，顯然與10～16 mm·a-1極快速隆升不符(Xu et al.，2000).實際上，初期(1960)跨青藏高原水準觀測誤差偏大，后兩期(1980、1990)復測精度較好，其顯示拉薩以南地區(qū)隆升約為2 mm·a-1(張全德等，2001)，與現(xiàn)有藏南一帶的GPS速率十分接近.

藏南地區(qū)隆升固然來自喜馬拉雅主沖斷裂的活動，但仍有一部分與印度板塊上地殼撕裂、增生到歐亞板塊底部有關，GPS觀測還難以區(qū)分這兩種機制的具體貢獻.此外沿雅魯藏布縫合線XZAR、XZZB、YARE表現(xiàn)為較大幅度隆升(>2 mm·a-1)，難以完全依據板塊邊界的構造活動來解釋，也許是雅魯藏布縫合線附近的新構造活動所致(例如2006—2008年仲巴地震)，現(xiàn)有測站密度還遠不足以證實.

雅魯藏布縫合線以北的GPS速率扣除非構造因素后，基本表現(xiàn)為下沉1～3 mm·a-1，雅魯藏布縫合線南北兩側由南升轉為北降.現(xiàn)有結果一方面與青藏高原末次冰期巨厚冰蓋消融導致的隆升分布特征不同，可以排除區(qū)域性大冰蓋的可能性；另一方面，即使忽略北極冰蓋消融的影響，雅魯藏布縫合線以北地區(qū)也有最大2 mm·a-1下降，與拉薩地塊(雅魯藏布縫合線至班公—怒江縫合線)的南北向縮短幅度相當(Wang and Shen，2020)，依據印度板塊插入青藏高原下方的傾角在雅魯藏布縫合帶逐漸變小不足以解釋這種幅度的下降，可能是拉薩、羌塘地塊東西向拉張變形相應導致地殼減薄，其幅度大大超過擠壓，導致高原內部塌陷，相對于西部站點，88°E以東區(qū)域GPS測站普遍下沉.

4 結論

青藏高原及喜馬拉雅地區(qū)連續(xù)GPS觀測以優(yōu)于1 mm·a-1的精度約束區(qū)域垂直位移場.GPS與精密水準和InSAR觀測對比驗證，比較完整地揭示了印度—歐亞碰撞帶垂直形變特征：從印度次大陸到青藏高原中部，地表經歷沉降-隆升-沉降，對應了地殼從增厚到減薄的轉換過程.除雅魯藏布縫合線以南的區(qū)域繼續(xù)保持隆升狀態(tài)，并向外擴展外，鑒于青藏高原內部GPS測站普遍下降1～2 mm·a-1，青藏高原的中南部可能已不再隆升.

致謝本文原始GPS數據來自中國大陸構造環(huán)境監(jiān)測網絡，中國氣象局、加州理工學院，美國鮑靈-格林州立大學富宇寧提供了GRACE修正值.審稿人與編輯對本文提出了寶貴的修改意見和建議.在此一并表示感謝.