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

1 中國地震局地震研究所， 地震大地測量重點(diǎn)實(shí)驗(yàn)室， 武漢 430071 2 中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院， 武漢 430074 3 國家基礎(chǔ)地理信息中心， 北京 100044 4 中國地震局地質(zhì)研究所， 北京 100029 5 武漢市應(yīng)急管理局， 武漢 430010

0 引言

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

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

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

1 資料

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

1.1 cGPS觀測

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

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

1.2 坐標(biāo)解算

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

圖2 GPS單日解垂向分量的時間序列圖中虛線為2004年和2015年發(fā)震時刻，綠色代表經(jīng)歷2015年尼泊爾地震的臺站.淺藍(lán)色代表震前臺站，深藍(lán)色為同一站點(diǎn)的GAMIT時序，紅色為震后布設(shè)的臺站.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 垂直坐標(biāo)時序分析

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

(1)

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

1.4 誤差模型

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

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

2 垂直運(yùn)動速率及精度評估

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

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

續(xù)表2

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

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

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

圖3 喜馬拉雅地區(qū)垂向速率剖線左圖：InSAR圖像、水準(zhǔn)線路和GPS站分布圖，紫色和藍(lán)色方框代表A、B剖線，色棒標(biāo)示InSAR和水準(zhǔn)點(diǎn)速率，黑色箭頭為相對歐亞板塊的GPS水平運(yùn)動速率(Wang and Shen，2020)；右圖：GPS(紅或紫色方框)，InSAR(棕色圓點(diǎn))和水準(zhǔn)(藍(lán)色圓點(diǎn))速率及地形合成圖.GPS速率未做任何修正.圖中近場和遠(yuǎn)場GPS采用兩種距離尺度表示，水準(zhǔn)和InSAR誤差棒代表1 mm·a-1的假設(shè)誤差，GPS誤差為真實(shí)誤差.綠色虛線為震間運(yùn)動模型的預(yù)測速率(Grandin et al.，2012)，棕色圓環(huán)代表此前的GPS結(jié)果(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 內(nèi)符合精度分析

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

圖4 4組GPS速率結(jié)果比對與此前公布的4組速率結(jié)果對比，西藏地區(qū)陸態(tài)站點(diǎn)用綠色點(diǎn)位表示.虛線陰影區(qū)代表互差在+/-1 mm·a-1以內(nèi)區(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.

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

2.3 模型對比及外符合精度

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

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

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

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

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

2.4 垂直速率修正

GPS測定的垂向變動主要包含三種運(yùn)動成分：1)板塊擠壓引起的構(gòu)造變形；2)地表水變化導(dǎo)致的升降，近年來喜馬拉雅冰川加速融化和恒河平原地下水抽取加劇(Rodell et al.，2009；Yi and Sun，2014)，將導(dǎo)致區(qū)域地殼大幅調(diào)整；3)末次冰期以來北半球巖石圈松弛變形及青藏高原冰后回彈.依據(jù)近年來GPS和水準(zhǔn)觀測，印度—?dú)W亞碰撞帶內(nèi)板塊水平移動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)，其垂向運(yùn)動受環(huán)境影響更大.

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

末次冰盛期導(dǎo)致的青藏高原隆升同樣存在不確定問題.按極端區(qū)域性冰川模型推測，青藏高原中部地殼冰后回彈的升幅可達(dá)4～7 mm·a-1(Kuhle，1998；Kaufmann，2005)，但另有研究表明青藏高原冰蓋規(guī)模不大、影響極小(Derbyshire，1991；Matsuo and Heki，2010；Liang et al.，2013)，本文也不加考慮.但北半球巖石圈松弛變形不容忽視，據(jù)現(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)，本文還考慮地心運(yùn)動對GPS速率的影響.根據(jù)全球大地測量結(jié)果(Argus，2012；Métivier et al.，2020)，末次冰期以來北半球巖石圈松弛變形以及內(nèi)部質(zhì)量遷移導(dǎo)致地球形心相對其質(zhì)心以0.5～1.0 mm·a-1速率向北極運(yùn)移，由于ITRF框架的原點(diǎn)為地球質(zhì)心，GPS實(shí)測速率相對于質(zhì)心，而構(gòu)造運(yùn)動相對于地球形心，地心的視運(yùn)動將導(dǎo)致中低緯度地區(qū)真實(shí)的垂向運(yùn)動偏大0.4 mm·a-1左右(Ding et al.，2019)，必須加以修正.本文按0.8 mm·a-1地心速率乘測站緯度正弦值計(jì)算修正值，修正值的精度大致在0.2 mm·a-1左右.綜合以上因素，現(xiàn)有GPS速率中非構(gòu)造影響在1 mm·a-1左右，按保守估計(jì)，修正值的不確定度大致在0.3～0.4 mm·a-1.

3 討論

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

3.1 印度次大陸

沿A剖線印度板塊相對歐亞板塊的水平速率從北到南為36～38 mm·a-1(Wang and Shen，2020)，如扣除非構(gòu)造因素，現(xiàn)有GPS表明印度板塊前緣部分已下沉2～3 mm·a-1，使其以3°～4°傾角向下插入歐亞板塊.根據(jù)尼泊爾地震的同震破裂產(chǎn)狀，位于低喜馬拉雅底部的印度板塊以7°傾角向北運(yùn)移(Wang and Fialko，2015)，遠(yuǎn)震接收函數(shù)成像顯示大約在高喜馬拉雅到雅魯藏布之間，印度板塊的傾角陡變到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結(jié)果還表明，德干高原北部為約500 km寬的前緣隆起區(qū)(從BHUP到JBPR)，德干高原南部(IISC、HYDE)幾乎沒有構(gòu)造性垂向運(yùn)動，GPS觀測與造山理論的預(yù)測一致(Molnar and Lyon-Caen，1988)，展示了印度板塊在青藏高原的重壓下，撓曲由淺入深的幾何結(jié)構(gòu)(圖5).

圖5 青藏高原南緣垂向運(yùn)動速度場GPS 實(shí)測速率經(jīng)地心(坐標(biāo)框架原點(diǎn)視)運(yùn)動、冰川均衡調(diào)整模擬和地表水變化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觀測到的印度次大陸前緣部分沉降可視為永久性、不可恢復(fù)變形，但GPS觀測的喜馬拉雅變形大部分是彈性的、與板塊邊界錯動及大震周期有關(guān)的變形.作為板塊邊界的喜馬拉雅主沖斷層最南緣的100 km段在兩次大震間完全閉鎖，低喜馬拉雅地區(qū)的震間變形緩慢，高喜馬拉雅地區(qū)快速隆升、變形累積(Bilham et al.，1997).這種隆升不是單調(diào)持續(xù)的，特大地震將導(dǎo)致大部分隆升回落，例如2015年尼泊爾地震時海拔3000 m以下的地帶全面隆升，而海拔3000 m以上地區(qū)整體下降(Elliott et al.，2016)，凸顯出地震在彈性變形從高海拔向低海拔區(qū)域的遷移過程中的關(guān)鍵作用.

GPS反映的彈性變形狀態(tài)與瞬時地震觀測的結(jié)果完全相反，加德滿都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年尼泊爾地震不能一次性將大震間積累的全部變形從高喜馬拉雅轉(zhuǎn)移至斯瓦里克，只有1934年比哈爾MW8.2地震乃至更大地震才有可能錯斷整個主沖斷裂，至最南端主前緣斷裂(Sapkota et al.，2013)，將震間絕大部分變形遷移至前緣地帶，轉(zhuǎn)換為永久性地形抬升.但高喜馬拉雅地帶GPS隆升速率變化與地形變化高度一致，表明這種轉(zhuǎn)換并不徹底，仍有部分彈性變形永久存留在高海拔地帶，維持其陡變地形(Meade，2010).

3.3 青藏高原

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

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

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

4 結(jié)論

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

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