用cGPS研究青藏高原南緣現(xiàn)今垂向變動(dòng)

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

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

0 引言

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

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

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

1 資料

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

1.1 cGPS觀測(cè)

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

圖1 cGPS站點(diǎn)分布及觀測(cè)時(shí)長(zhǎng)Fig.1 Distribution and observation time span of cGPS stations

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

續(xù)表1

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

1.2 坐標(biāo)解算

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

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

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

(1)

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

1.4 誤差模型

GIPSY給出的單日大地高程的解算精度(名義誤差)約為3～4 mm，主要反映了觀測(cè)當(dāng)日載波相位中的隨機(jī)噪聲.其他因素，如星歷、電離層、對(duì)流層折射、測(cè)站周邊溫度變化、多路徑等因素引起的各種偏差不在其中，這些偏差對(duì)于單日坐標(biāo)而言屬于系統(tǒng)性的，相鄰數(shù)日或數(shù)周內(nèi)系統(tǒng)偏差變化不大，對(duì)站點(diǎn)坐標(biāo)的擾動(dòng)一致，在時(shí)序中表現(xiàn)為有色噪聲(不同時(shí)間觀測(cè)值的統(tǒng)計(jì)性質(zhì)相關(guān)).如簡(jiǎn)單將實(shí)測(cè)數(shù)據(jù)誤差視為白噪聲，雖不影響對(duì)各個(gè)信號(hào)的提取，但大大高估其解算精度，故本文時(shí)序分析中的誤差模型采用白噪聲+有色噪聲的組合方式(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í)際狀況.每個(gè)測(cè)站具體的有色噪聲模型主要依據(jù)其擬合殘差的功率譜指數(shù)而定，相關(guān)計(jì)算采用Hector軟件，同時(shí)給出站點(diǎn)垂向速率、周期項(xiàng)參數(shù)及其不確定性(表2).

2 垂直運(yùn)動(dòng)速率及精度評(píng)估

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

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

續(xù)表2

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

GPS測(cè)站的垂直運(yùn)動(dòng)特征概括為：恒河平原及喜馬拉雅山前(海拔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)相對(duì)印度次大陸最大隆升接近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ù)顯示高喜馬拉雅與低喜馬拉雅相對(duì)前陸盆地的抬升分別約為4±1 mm·a-1和2.0±0.3 mm·a-1(Herman et al.，2010).尼泊爾西部的InSAR觀測(cè)更凸顯這種升降變化(圖3)，高喜馬拉雅相對(duì)于恒河平原最大隆升達(dá)到8 mm·a-1(Grandin et al.，2012).盡管本文直接解算的尼泊爾cGPS垂向速率場(chǎng)及印度西北地區(qū)cGPS觀測(cè)結(jié)果(Yadav et al.，2019)的空間分辨率較低，但結(jié)合InSAR、精密水準(zhǔn)，現(xiàn)有大地測(cè)量綜合結(jié)果完整反映了大陸俯沖帶前陸盆地沉降、山前快速抬升、造山帶后緣隆升延展、幅度急劇降低這一典型特征.

圖3 喜馬拉雅地區(qū)垂向速率剖線左圖：InSAR圖像、水準(zhǔn)線路和GPS站分布圖，紫色和藍(lán)色方框代表A、B剖線，色棒標(biāo)示InSAR和水準(zhǔn)點(diǎn)速率，黑色箭頭為相對(duì)歐亞板塊的GPS水平運(yùn)動(dòng)速率(Wang and Shen，2020)；右圖：GPS(紅或紫色方框)，InSAR(棕色圓點(diǎn))和水準(zhǔn)(藍(lán)色圓點(diǎn))速率及地形合成圖.GPS速率未做任何修正.圖中近場(chǎng)和遠(yuǎn)場(chǎng)GPS采用兩種距離尺度表示，水準(zhǔn)和InSAR誤差棒代表1 mm·a-1的假設(shè)誤差，GPS誤差為真實(shí)誤差.綠色虛線為震間運(yùn)動(dòng)模型的預(yù)測(cè)速率(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é)果逐一對(duì)比(圖4).加州理工最早公布測(cè)站速率誤差一般為1～3 mm·a-1(Ader et al.，2012)，比本文大3倍.加州理工除顧及有色噪聲外，還考慮了時(shí)序不連續(xù)可能引入的階躍誤差，因此總體上精度估計(jì)偏保守.其他三組與本文采用完全一致的誤差評(píng)估算法(Fu and Freymueller，2012；Liang et al.，2013；Pan et al.，2018)，圖4顯示本文速率精度總體上最優(yōu)，這得益于某些測(cè)站已有了更長(zhǎng)時(shí)序，當(dāng)時(shí)間跨度大于5年時(shí)(XZAR、XZNM、SNDL)，PPP給出的垂向速率精度在0.5～0.8 mm·a-1，大于10年時(shí)(KKN4、KAWA、CHLM)，誤差進(jìn)一步降低到0.3 mm·a-1，而超過(guò)20年觀測(cè)(LHAS、HYDE、IISC)，誤差甚至可低至0.1 mm·a-1.與此對(duì)照，陸態(tài)網(wǎng)絡(luò)觀測(cè)時(shí)長(zhǎng)12年的區(qū)域站精度一般低于1 mm·a-1(Liang et al.，2013)，如此估計(jì)25年觀測(cè)歷史的流動(dòng)站有望將誤差降低到0.5 mm·a-1以內(nèi).

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

2.3 模型對(duì)比及外符合精度

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

跨喜馬拉雅的InSAR干涉圖像覆蓋83°～84°經(jīng)度帶，從低喜馬拉雅延伸到雅魯藏布，南北長(zhǎng)250 km(圖3)，是處理2003—2010年間29幅降軌原始圖像，將不同時(shí)段、獨(dú)立干涉圖疊加而成(Grandin et al.，2012).變形起算點(diǎn)位于圖像最北端，變形信號(hào)從北到南累積，南北兩端像素點(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í)際上更多地反映了垂向變動(dòng)的分布，就A剖線而言，LOS 速率與垂向速率的偏差不超過(guò)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)加德滿都終于高喜馬拉雅山前的西藏樟木口岸，總長(zhǎng)350 km，沿線高程起伏不超過(guò)1500 m, 該線路1977—1991年間觀測(cè)，觀測(cè)誤差為1.1 mm·km-1，按誤差傳播率估計(jì)，水準(zhǔn)線路南、北兩端的速率差相對(duì)精度為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)的一個(gè)測(cè)點(diǎn)，相應(yīng)推算各自路線上其他水準(zhǔn)點(diǎn)的ITRF框架速率值.

對(duì)比GPS、精密水準(zhǔn)、InSAR資料可見(jiàn)(圖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í)測(cè)值對(duì)比，則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觀測(cè)誤差較大或KLDN疊加了前沿?cái)嗔炎冃蔚挠绊?與模型值(Grandin et al., 2012)對(duì)比，只有KLDN、YARE、XZZB有近2 mm·a-1的偏差，其中YARE、XZZB與所在區(qū)域的InSAR觀測(cè)值都系統(tǒng)性偏離了模型值，且GPS、InSAR的變化趨勢(shì)也相同，這種異常變形可能與2006—2008年仲巴序列地震有關(guān).整體而言，本文比早期結(jié)果(Fu and Freymueller，2012)更接近模型值.

2.4 垂直速率修正

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

地表水變化引起的均衡調(diào)整一般根據(jù)時(shí)變重力場(chǎng)推算的區(qū)域地表載荷變化及彈性地球變形響應(yīng)來(lái)估算(Fu and Freymueller，2012).本文利用最新的GRACE觀測(cè)結(jié)果，假定時(shí)變重力場(chǎng)完全是地表水變化決定，忽略其他因素的貢獻(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)性偏差，說(shuō)明用不同濾波算法平滑GRACE資料，給出的修正值存在至少0.2～0.3 mm·a-1的不確定性.

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

與以往研究不同(Liang et al.，2013；Hao et al.，2016；Pan et al.，2018)，本文還考慮地心運(yùn)動(dòng)對(duì)GPS速率的影響.根據(jù)全球大地測(cè)量結(jié)果(Argus，2012；Métivier et al.，2020)，末次冰期以來(lái)北半球巖石圈松弛變形以及內(nèi)部質(zhì)量遷移導(dǎo)致地球形心相對(duì)其質(zhì)心以0.5～1.0 mm·a-1速率向北極運(yùn)移，由于ITRF框架的原點(diǎn)為地球質(zhì)心，GPS實(shí)測(cè)速率相對(duì)于質(zhì)心，而構(gòu)造運(yùn)動(dòng)相對(duì)于地球形心，地心的視運(yùn)動(dòng)將導(dǎo)致中低緯度地區(qū)真實(shí)的垂向運(yùn)動(dòng)偏大0.4 mm·a-1左右(Ding et al.，2019)，必須加以修正.本文按0.8 mm·a-1地心速率乘測(cè)站緯度正弦值計(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è)不同的構(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)換對(duì)應(yīng)了GPS垂向運(yùn)動(dòng)從下沉到隆升、再到下沉的空間展布.

3.1 印度次大陸

沿A剖線印度板塊相對(duì)歐亞板塊的水平速率從北到南為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)幾乎沒(méi)有構(gòu)造性垂向運(yùn)動(dòng)，GPS觀測(cè)與造山理論的預(yù)測(cè)一致(Molnar and Lyon-Caen，1988)，展示了印度板塊在青藏高原的重壓下，撓曲由淺入深的幾何結(jié)構(gòu)(圖5).

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

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

3.3 青藏高原

最早的跨青藏高原的精密水準(zhǔn)顯示，藏南(喜馬拉雅北坡至雅魯藏布)地區(qū)近10 mm·a-1隆升(張青松等，1991)，這種大幅度現(xiàn)代隆升得到早期流動(dòng)GPS觀測(cè)的佐證(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)觀測(cè)誤差偏大，后兩期(1980、1990)復(fù)測(cè)精度較好，其顯示拉薩以南地區(qū)隆升約為2 mm·a-1(張全德等，2001)，與現(xiàn)有藏南一帶的GPS速率十分接近.

藏南地區(qū)隆升固然來(lái)自喜馬拉雅主沖斷裂的活動(dòng)，但仍有一部分與印度板塊上地殼撕裂、增生到歐亞板塊底部有關(guān)，GPS觀測(cè)還難以區(qū)分這兩種機(jī)制的具體貢獻(xiàn).此外沿雅魯藏布縫合線XZAR、XZZB、YARE表現(xiàn)為較大幅度隆升(>2 mm·a-1)，難以完全依據(jù)板塊邊界的構(gòu)造活動(dòng)來(lái)解釋，也許是雅魯藏布縫合線附近的新構(gòu)造活動(dòng)所致(例如2006—2008年仲巴地震)，現(xiàn)有測(cè)站密度還遠(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)致地殼減薄，其幅度大大超過(guò)擠壓，導(dǎo)致高原內(nèi)部塌陷，相對(duì)于西部站點(diǎn)，88°E以東區(qū)域GPS測(cè)站普遍下沉.

4 結(jié)論

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

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