汶川地震震后GNSS形變分析

余建勝，趙 斌，譚 凱，王東振

中國(guó)地震局地震研究所地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430071

印度板塊以每年約5 cm的速度向北漂移，持續(xù)碰撞并擠壓中國(guó)大陸所在的歐亞板塊，致使青藏高原內(nèi)部物質(zhì)向東南方向側(cè)流[1-7]，側(cè)流過(guò)程中遇到其東邊穩(wěn)固的四川盆地阻擋，從而在龍門山斷裂帶上積累了大量應(yīng)力。一旦斷層面累積應(yīng)力超過(guò)斷層可承受最大閾值，會(huì)觸發(fā)地震事件，釋放出積蓄的巨大能量，引發(fā)各種自然災(zāi)害。2008年5月12日，位于青藏高原東北緣的龍門山斷裂帶上，以汶川縣映秀鎮(zhèn)為震中發(fā)生了一次具有右旋走滑兼逆沖運(yùn)動(dòng)性質(zhì)的Mw7.9級(jí)特大地震事件[1]。

大震發(fā)生后數(shù)年至數(shù)十年甚至更長(zhǎng)時(shí)間尺度內(nèi)，在震區(qū)及周邊會(huì)觀測(cè)到隨時(shí)間呈指數(shù)或?qū)?shù)衰減的地表位移[8-12]，即震后位移。引起震后位移效應(yīng)的機(jī)制主要有孔隙回彈、余滑和黏彈性松弛。已有文獻(xiàn)[13—20]研究發(fā)現(xiàn)，由地震導(dǎo)致流體流動(dòng)引發(fā)的震后孔隙回彈通常在含有大量地下水的區(qū)域?qū)Υ怪毕蛐巫冇休^大影響；震后余滑通常只局限在震后幾個(gè)月到1年左右的短時(shí)間尺度內(nèi)，且主要作用于同震破裂下傾的近場(chǎng)區(qū)域；震后黏彈性松弛對(duì)中、遠(yuǎn)場(chǎng)區(qū)域有持續(xù)影響，其作用時(shí)間可以延續(xù)幾年、幾十年甚至更久。

震后形變因包含有地球介質(zhì)的流變信息，蘊(yùn)含著地殼運(yùn)動(dòng)和深部構(gòu)造信息，因此是研究地球介質(zhì)流變性質(zhì)及巖石圈動(dòng)力學(xué)關(guān)系的一種重要手段[18-21]。以震后位移為地表約束，結(jié)合不同的震后形變機(jī)制，可以對(duì)斷層的動(dòng)力學(xué)性質(zhì)、地殼深部的流變結(jié)構(gòu)等地球動(dòng)力學(xué)問(wèn)題有進(jìn)一步認(rèn)知。由于地球內(nèi)部的未知性，近年來(lái)，隨著大地測(cè)量觀測(cè)技術(shù)不斷發(fā)展進(jìn)步，有關(guān)科研人員開(kāi)始利用GNSS觀測(cè)資料研究地球內(nèi)部介質(zhì)的差異性。發(fā)生在青藏高原東邊界龍門山斷裂帶的汶川Mw7.9級(jí)地震為研究青藏高原東緣的巖石圈流變結(jié)構(gòu)及龍門山斷裂帶的摩擦性質(zhì)提供了難得的機(jī)會(huì)。文獻(xiàn)[16]研究汶川震后8個(gè)月47個(gè)測(cè)站震后GNSS形變序列。文獻(xiàn)[14]計(jì)算汶川震后1年16個(gè)測(cè)站三維位移，并以此進(jìn)行反演約束。文獻(xiàn)[21]利用GPS和InSAR數(shù)據(jù)聯(lián)合反演同震滑動(dòng)分布和黏彈性松弛，基于模型反演給出汶川震后2年時(shí)間尺度的InSAR形變場(chǎng)分布。文獻(xiàn)[20]采用震后1.5年的InSAR資料和震后一年少量GNSS數(shù)據(jù)對(duì)汶川地震震后形變的一級(jí)特征進(jìn)行了分析，發(fā)現(xiàn)川西高原和四川盆地的黏滯系數(shù)存在顯著差異，相差兩個(gè)數(shù)量級(jí)。盡管汶川地震已經(jīng)過(guò)去十年了，但現(xiàn)有成果主要集中在震后早期1～2年短時(shí)間尺度內(nèi)，缺少對(duì)汶川震后形變長(zhǎng)時(shí)間尺度持續(xù)跟蹤分析，并且所使用的GNSS數(shù)據(jù)主要集中在近場(chǎng)，缺乏中遠(yuǎn)場(chǎng)觀測(cè)資料，對(duì)黏滯系數(shù)的約束存在不確定性。本文根據(jù)2008年汶川地震震后龍門山斷裂帶上、下盤區(qū)域收集到的109個(gè)測(cè)站觀測(cè)資料，通過(guò)高精度GNSS數(shù)據(jù)處理并對(duì)震前、震后時(shí)間序列進(jìn)行分析，獲取各測(cè)站震后位移。根據(jù)位移場(chǎng)分布，進(jìn)一步研究龍門山斷裂上、下盤青藏高原東部和四川盆地深部介質(zhì)差異性，推算汶川地震在未來(lái)幾十年內(nèi)對(duì)周邊區(qū)域地殼形變的影響。

1 GNSS數(shù)據(jù)處理與時(shí)間序列分析

1.1 站點(diǎn)選取與數(shù)據(jù)處理

龍門山斷裂帶是南北地震帶的重要組成部分，位于巴顏喀拉塊體和四川盆地交界處，有歷史地震記載以來(lái)，地震多發(fā)生于斷裂帶中南段及其北邊岷江斷裂間[1]。本文使用109個(gè)GNSS測(cè)站觀測(cè)資料，時(shí)間跨度為1999—2015年，GNSS測(cè)站觀測(cè)資料統(tǒng)計(jì)見(jiàn)表1。龍門山斷裂帶區(qū)域構(gòu)造背景與GNSS測(cè)站分布情況如圖1所示，其中，連續(xù)觀測(cè)站點(diǎn)共計(jì)28個(gè)，流動(dòng)觀測(cè)站點(diǎn)共計(jì)81個(gè)。采用高精度數(shù)據(jù)處理軟件GAMIT/GLOBK10.4對(duì)收集的109個(gè)GNSS站點(diǎn)進(jìn)行高精密數(shù)據(jù)處理，獲取單日松弛解，并引進(jìn)IGS站作為公共測(cè)站聯(lián)合解算，通過(guò)坐標(biāo)轉(zhuǎn)換將解算結(jié)果歸算得到ITRF08參考框架下的站點(diǎn)坐標(biāo)[23]。

表1 GNSS觀測(cè)資料統(tǒng)計(jì)

1.2 時(shí)間序列分析

GNSS坐標(biāo)時(shí)間序列包含站點(diǎn)長(zhǎng)期構(gòu)造運(yùn)動(dòng)、周期性季節(jié)變化、已知和未知階躍及震后形變等。長(zhǎng)期觀測(cè)的GNSS單站、單分量坐標(biāo)時(shí)間序列可用式(1)表示

(1)

式中，x(t)為歷元t時(shí)刻測(cè)站某一分量坐標(biāo)位置；C1為坐標(biāo)初始位置；C2為測(cè)站長(zhǎng)期線性速率；C3和θ為周年變化的振幅與相位；C4和φ為半周年變化的振幅與相位；Di是ti時(shí)刻因地震、儀器變更等原因產(chǎn)生的階躍；H為階梯函數(shù)；f為震后形變。

圖1 龍門山斷裂帶區(qū)域構(gòu)造背景及GNSS站點(diǎn)分布Fig.1 Local tectonic setting and the distribution of GNSS stations located in Longmen shan fault zone

震后GNSS觀測(cè)到的坐標(biāo)時(shí)間序列包含震后形變f，從GNSS坐標(biāo)時(shí)間序列中獲取震后形變，需扣除測(cè)站長(zhǎng)期構(gòu)造運(yùn)動(dòng)、測(cè)站受周邊地震與儀器更換等產(chǎn)生的階躍、周期性季節(jié)性變化等影響，其中季節(jié)性變化主要影響垂直方向，本文研究水平方向坐標(biāo)時(shí)間序列，故季節(jié)性變化忽略不計(jì)。在汶川地震震后形變研究過(guò)程中，根據(jù)1997—2007年中國(guó)大陸連續(xù)跟蹤站和區(qū)域站的平均運(yùn)動(dòng)速率作為本文地殼運(yùn)動(dòng)長(zhǎng)期運(yùn)動(dòng)速率構(gòu)造背景值[23]，采用小區(qū)域均勻應(yīng)變率方法對(duì)震后新建GNSS觀測(cè)站點(diǎn)進(jìn)行內(nèi)插處理，得到全部測(cè)站在震前的長(zhǎng)期構(gòu)造運(yùn)動(dòng)速率值。

本文研究2010—2015年期間震后形變，該時(shí)段內(nèi)于2013年4月20日在龍門山前緣構(gòu)造帶南段雅安地區(qū)發(fā)生Mw6.6級(jí)蘆山地震。通過(guò)計(jì)算發(fā)現(xiàn)，本次所選取站點(diǎn)中有9個(gè)測(cè)站受蘆山地震同震影響明顯。圖2(a)列舉出了其中2個(gè)連續(xù)站點(diǎn)SCXJ和SCMX扣除長(zhǎng)期線性速率后的坐標(biāo)時(shí)間序列，圖2(b)列出了其中10個(gè)區(qū)域站點(diǎn)扣除長(zhǎng)期線性速率后的坐標(biāo)時(shí)間序列，圖2中可以看出所示測(cè)站具有明顯的震后形變趨勢(shì)。

2 震后形變計(jì)算

描述震后形變過(guò)程的函數(shù)模型主要有3類，分別為冪函數(shù)、指數(shù)函數(shù)和對(duì)數(shù)函數(shù)，3種函數(shù)模型各有優(yōu)勢(shì)。分析余震活動(dòng)時(shí)常用冪指數(shù)函數(shù)；在多數(shù)巖體力學(xué)試驗(yàn)中，通常用指數(shù)函數(shù)模型來(lái)描述黏彈性松弛特性；對(duì)于理想的黏彈性Maxwell體，可用對(duì)數(shù)函數(shù)模型來(lái)描述其震后松弛變化過(guò)程。假設(shè)川西地區(qū)地殼結(jié)為Maxwell體，本文震后形變衰減過(guò)程可采用如下對(duì)數(shù)函數(shù)描述[20]

(2)

式中，t為震后某一觀測(cè)時(shí)刻；tEQ為地震發(fā)生時(shí)刻；D為測(cè)站震后形變幅度值；τ為震后松弛時(shí)間；f(t,τ)表示扣除長(zhǎng)期線性速率、已知或未知階躍等影響后反映震后形變的坐標(biāo)時(shí)間序列。不同測(cè)站震后衰減呈現(xiàn)不一致性，本文綜合考慮，選取時(shí)間序列質(zhì)量比較好且數(shù)據(jù)覆蓋時(shí)段相對(duì)較長(zhǎng)的站點(diǎn)來(lái)獲取可靠的震后松弛時(shí)間τ。通過(guò)計(jì)算發(fā)現(xiàn)，測(cè)站取各自最優(yōu)的τ值，與取最佳的公共τ=34 d時(shí)，計(jì)算的2010—2015年震后累計(jì)形變值間差異性很小，大部分在0.5 mm以內(nèi)，最大差異不超過(guò)1 mm。根據(jù)擬合殘差最小原則求解出的最優(yōu)τ值為34天，部分測(cè)站時(shí)間序列擬合如圖2所示，計(jì)算得到2010—2015年間各測(cè)站累計(jì)震后水平位移，如圖3所示。

根據(jù)圖3所示，在2010—2015年期間震后形變量較大的站點(diǎn)有H045、J417、H050等，其水平方向累計(jì)形變達(dá)到5～7 cm，主要分布在龍門山斷裂帶中段的上盤近場(chǎng)區(qū)域。龍門山斷裂上盤，鮮水河斷裂以東及岷江斷裂以西區(qū)域中、遠(yuǎn)場(chǎng)震后形變普遍較大，達(dá)到2～5 cm，并隨距離的增加呈衰減趨勢(shì)。龍門山斷裂上盤、岷江斷裂以東，除青川附近靠近地震破裂帶的少數(shù)站點(diǎn)形變稍大，如H010、H032等，其余站點(diǎn)的震后形變都普遍較小，基本在2 cm范圍內(nèi)，且有沿?cái)嗔褞虮狈较虺誓鏁r(shí)針旋轉(zhuǎn)的運(yùn)動(dòng)趨勢(shì)，與同震形變運(yùn)動(dòng)趨勢(shì)比較一致。龍門山斷裂帶南段與鮮水河斷裂、安寧河斷裂交界的三岔口區(qū)域，大部分站點(diǎn)震后形變量不足1 cm，且沒(méi)有固定的運(yùn)動(dòng)趨勢(shì)，說(shuō)明該區(qū)域內(nèi)震后形變比較微弱。相對(duì)上盤，龍門山斷裂下盤計(jì)算的測(cè)站震后形變則非常微小，除個(gè)別測(cè)站觀測(cè)到了比較明顯的震后形變，其余測(cè)站震后形變基本都在5 mm以內(nèi)，這與下盤測(cè)站位于相對(duì)穩(wěn)定的四川盆地上比較一致。對(duì)圖3所示結(jié)果進(jìn)行分析發(fā)現(xiàn)，以龍門山斷裂帶為界，斷裂帶上盤震后形變普遍比下盤大，表明龍門山斷裂帶上、下盤具有很強(qiáng)的不對(duì)稱性，反映青藏高原東部和四川盆地深部流變性質(zhì)的顯著差異性。

圖2 震后變化趨勢(shì)Fig.2 Change trend after earthquake

3 地殼結(jié)構(gòu)反演

大地震發(fā)生后引發(fā)的震后黏彈性松弛位移，是由地震破裂產(chǎn)生的同震應(yīng)力作用于下地殼/上地幔的軟流物質(zhì)上，造成深部物質(zhì)的流變過(guò)程，從而引起彈性層持續(xù)、緩慢的地表位移[16]。黏彈性是地球介質(zhì)的一種重要性質(zhì)，黏滯系數(shù)是描述介質(zhì)黏彈性質(zhì)的重要參數(shù)，介質(zhì)的黏彈性質(zhì)對(duì)地球動(dòng)力學(xué)過(guò)程及其演化具有重要作用[17]。黏彈性質(zhì)是隨時(shí)間的變化而不斷變化的，某一時(shí)間內(nèi)計(jì)算出的黏滯系數(shù)只能反映特定時(shí)間段內(nèi)該區(qū)域的流變結(jié)構(gòu)。

圖3 汶川震后2010—2015年累計(jì)震后水平形變Fig.3 Cumulated horizontal postseismic deformation after Wenchuan earthquake during 2010—2015

本文擬研究汶川地震震后2～7年，即2010—2015年間震后形變，認(rèn)為影響震后形變的主要機(jī)制為震后黏彈性松弛。已有相關(guān)研究表明，汶川地震震后松弛主要發(fā)生在上盤地下30～60 km深處[18]。根據(jù)全球地殼模型crust1.0，并參考前人研究成果[14-23]，建立龍門山上盤地殼分層結(jié)構(gòu)模型(表2)，h和H分別為中上地殼、下地殼/上地幔的厚度。根據(jù)介質(zhì)的差異性將彈性的中上地殼劃分為多層，賦予不同的介質(zhì)參數(shù)。彈性層厚度h的搜索范圍為30～60 km，黏滯系數(shù)搜索范圍為1017～1021pa·s。采用Wang等的PSGRN/PSCMP程序[24]模擬計(jì)算汶川震后，龍門山上盤區(qū)域因黏彈性應(yīng)力松弛所引發(fā)的地表位移，通過(guò)對(duì)參數(shù)區(qū)間(h，η)進(jìn)行搜索計(jì)算，并計(jì)算擬合誤差RMS，擬合誤差取最小值所對(duì)應(yīng)的彈性層厚度和黏滯系數(shù)即為本文龍門山斷裂上盤區(qū)域最優(yōu)估值。

表2 龍門山斷裂上盤地殼分層結(jié)構(gòu)模型

注：Vp、Vs分別為介質(zhì)中的P波和S波波速。

文中采用2010—2015年間累計(jì)震后形變，約束川西高原的巖石圈流變一級(jí)結(jié)構(gòu)，認(rèn)為該時(shí)段的震后形變主要受震后黏彈性松弛效應(yīng)控制。此外，為了避免靠近斷層的GPS測(cè)站受震后余滑影響，還剔除距同震破裂區(qū)域約100 km以內(nèi)的近場(chǎng)測(cè)站，最終選取龍門山斷裂上盤中、遠(yuǎn)場(chǎng)40個(gè)測(cè)站進(jìn)行反演約束并估計(jì)最優(yōu)彈性層厚度和黏滯系數(shù)。震后形變的大小和分布特征除了受巖石圈流變結(jié)構(gòu)影響外，還與同震破裂模型息息相關(guān)。為此，筆者考查了3個(gè)不同同震破裂模型對(duì)模擬結(jié)果的影響。模型一為文獻(xiàn)[7]公布的斷層模型參數(shù)、模型二為文獻(xiàn)[25]公布的破裂模型、模型三為USGS公布的單斷層模型(http:∥earthquake.usgs.gov/eqcenter/eqintheneus/2008/us2008ryan/timite_fault.php,2008.)。其中模型一、二是基于GNSS、InSAR等大地測(cè)量資料約束反演得到的斷層滑動(dòng)模型，模型三是基于地震波數(shù)據(jù)模擬的斷層破裂及靜態(tài)滑動(dòng)分布。分別采用3種破裂模型進(jìn)行模擬計(jì)算，得到3種不同同震模型對(duì)應(yīng)的最優(yōu)彈性層厚度和黏滯系數(shù)估值(見(jiàn)表3)，以及震后形變觀測(cè)值與模擬值對(duì)比情況，如圖4所示。

表3 3種同震破裂模型估計(jì)的彈性層厚度和黏滯系數(shù)

圖4顯示3種不同破裂模型計(jì)算的震后形變擬合值和模擬值均具有較好的一致性，除局部個(gè)別站點(diǎn)外，整體差異性不大。從殘差分布圖可以看出，龍門山斷裂上盤、鮮水河斷裂以東及岷江斷裂以西區(qū)域測(cè)站模擬結(jié)果一致性較好，其中模型三在該區(qū)域內(nèi)的擬合一致性最佳。岷江斷裂及其周邊區(qū)域地質(zhì)構(gòu)造環(huán)境極其復(fù)雜，斷裂以東區(qū)域測(cè)站除H024、H007等個(gè)別站點(diǎn)模擬值相對(duì)觀測(cè)值偏小，其余測(cè)站均偏大，推測(cè)可能是岷江斷裂東、西兩側(cè)地下介質(zhì)流變性質(zhì)的差異性導(dǎo)致，岷江以東區(qū)域黏滯系數(shù)比岷江以西可能要稍大；H007、H004等測(cè)站擬合結(jié)果較差，可能是受地殼分層結(jié)構(gòu)的影響，也可能是受區(qū)域站觀測(cè)精度影響。鮮水河斷裂周邊測(cè)站距離汶川震中較遠(yuǎn)，不同破裂模型引起的震后形變變化不大；測(cè)站震后位移模擬值普遍小于觀測(cè)值，由于觀測(cè)值已經(jīng)扣除2013年蘆山Mw6.6地震同震形變影響，因而可能是該區(qū)域深部流變性質(zhì)差異，以及鮮水河斷裂長(zhǎng)期構(gòu)造運(yùn)動(dòng)導(dǎo)致。不同模型遠(yuǎn)場(chǎng)震后形變基本無(wú)差異，模型的差異主要影響上盤近場(chǎng)及岷江斷裂以東區(qū)域，模型一考慮滑脫層的存在，故J412、J413、H037等測(cè)站黏彈性模擬值被低估；另更加準(zhǔn)確的滑動(dòng)模型在岷江斷裂附近模擬效果改善有限，除H034擬合較好外，其余站點(diǎn)殘差依然較大。盡管本文反演前已剔除同震破裂區(qū)域100 km以內(nèi)的測(cè)站，但通過(guò)計(jì)算發(fā)現(xiàn)，模型一、二中少數(shù)靠近破裂帶的測(cè)站其黏彈性模擬值與觀測(cè)值仍然不一致，模擬值被明顯低估，尤其是J413、J412；距同震破裂區(qū)域較近測(cè)站的殘差近似認(rèn)為是由余滑引起，估算模型一中測(cè)站最大余滑量占比約30%，模型二最大余滑量占比約25%。綜合圖4得出，本文所研究的汶川地震破裂區(qū)100 km以外區(qū)域在2010—2015年的震后地殼形變以黏彈性松弛效應(yīng)為主，震后余滑只在同震破裂有限區(qū)域內(nèi)對(duì)少數(shù)測(cè)站有較明顯影響。

4 分析與討論

4.1 流變結(jié)構(gòu)

圖3揭示川西高原和四川盆地震后形變呈現(xiàn)強(qiáng)烈的非對(duì)稱性，反映青藏高原東部和四川盆地深部流變性質(zhì)的顯著差異，震后形變主要集中在龍門山斷裂上盤區(qū)域，下盤四川盆地幾乎沒(méi)有觀測(cè)到較明顯的震后形變。四川盆地內(nèi)的GNSS測(cè)站震后形變遠(yuǎn)小于川西高原觀測(cè)值，說(shuō)明川西高原的黏彈性系數(shù)相對(duì)較高，這與文獻(xiàn)[20]分析震后1年GPS觀測(cè)值得到的認(rèn)識(shí)是一致的。本文選取的3種同震破裂模型計(jì)算出的震后黏彈性松弛效應(yīng)可以較好地解釋2010—2015年期間龍門山斷裂上盤中遠(yuǎn)場(chǎng)震后形變，但不能很好地解釋該時(shí)間段內(nèi)靠近岷江斷裂附近的測(cè)站，如SCSP、H030等。岷江斷裂附近測(cè)站震后形變模擬值分布與文獻(xiàn)[20]的比較相似，測(cè)站中出現(xiàn)的觀測(cè)值與模擬值方位角偏差可能是受青藏高原東部橫向黏度不均勻影響。同時(shí)，龍門山斷裂以北區(qū)域發(fā)育有NWW向的東昆侖斷裂、近南北向的岷江斷裂、NNW向虎牙斷裂等，其中東昆侖斷裂東端向東南方向呈馬尾狀散開(kāi)，該區(qū)域構(gòu)造比較復(fù)雜，3種不同模型計(jì)算的擬合值與模擬值對(duì)比趨勢(shì)比較一致，可能與該區(qū)域自身長(zhǎng)期構(gòu)造運(yùn)動(dòng)相關(guān)。

文獻(xiàn)[28]采用三維有限元方法模擬了汶川震后初期14天的黏彈性松弛效應(yīng)，給出成都平原與川西高原地區(qū)中下地殼黏滯系數(shù)分別為9×1018pa·s和4×1017pa·s，得出龍門山斷裂兩側(cè)黏滯系數(shù)差異至少有1個(gè)量級(jí)。文獻(xiàn)[21]采用簡(jiǎn)單的二元分層、橫向均勻的地球模型，通過(guò)InSAR數(shù)據(jù)約束，采用同震和震后形變聯(lián)合模型反演，初步獲得中下地殼黏滯系數(shù)下限為為2.0×1018pa·s。該黏滯系數(shù)結(jié)果與文獻(xiàn)[28]川西高原的結(jié)果比較一致，但與成都平原地區(qū)相差一個(gè)量級(jí)。文獻(xiàn)[20]采用震后1年的GPS觀測(cè)數(shù)據(jù)以及震后一年半的InSAR觀測(cè)資料系統(tǒng)研究了汶川地震震后形變機(jī)理，通過(guò)三維有限元構(gòu)建龍門山斷裂兩側(cè)黏滯系數(shù)差異的巖石圈流變結(jié)構(gòu)，大地測(cè)量資料約束的四川盆地的黏滯系數(shù)不低于1020pa·s，而川西高原下地殼瞬態(tài)和穩(wěn)態(tài)黏滯系數(shù)分別為4.4×1017pa·s和1.0×1018pa·s。文獻(xiàn)[14]采用兩層的地殼結(jié)構(gòu)分層模型，以震后1年16個(gè)測(cè)站的GNSS形變?yōu)榧s束，得到彈性層最佳黏滯系數(shù)為1.8×1019pa·s。本文選取的3種不同破裂模型計(jì)算給出的龍門山斷裂上盤川西高原的黏滯系數(shù)量級(jí)均為1019pa·s，與文獻(xiàn)[28]和文獻(xiàn)[14]的結(jié)果比較吻合。

本文篩選上盤距汶川地震破裂區(qū)100 km以外的40個(gè)測(cè)站進(jìn)行反演約束，為了驗(yàn)證測(cè)站選取對(duì)計(jì)算本文黏彈系數(shù)結(jié)果的影響，同時(shí)增加了同震破裂區(qū)域10個(gè)近場(chǎng)測(cè)站的反演結(jié)果。以模型二為例，當(dāng)RMS=16.8 mm時(shí)，h=31 km，η=1.6×1019pa·s，彈性層厚度和黏滯系數(shù)取最優(yōu)估值。通過(guò)計(jì)算發(fā)現(xiàn)，增加100 km以內(nèi)的近場(chǎng)測(cè)站一起反演計(jì)算得到的黏滯系數(shù)值和扣除近場(chǎng)測(cè)站單獨(dú)反演獲取的結(jié)果變化不大。通過(guò)對(duì)格網(wǎng)搜索值結(jié)果進(jìn)行分析發(fā)現(xiàn)，黏滯系數(shù)比彈性層厚度對(duì)反演結(jié)果要更加敏感，本文研究給出川西高原地區(qū)下地殼/上地幔黏彈性層黏滯系數(shù)下限為1019pa·s。根據(jù)增加10個(gè)近場(chǎng)測(cè)站反演約束計(jì)算的結(jié)果，對(duì)比圖4分析發(fā)現(xiàn)，龍門山上盤余滑區(qū)域主要分布在汶川地震破裂帶100 km以內(nèi)，且有向東北方向的余滑趨勢(shì)，受影響的測(cè)站有H010、H033、J417、SCMX、H045、H050，主要分布在四川茂縣、青川、映秀等區(qū)域。

4.2 震后黏彈性松弛效應(yīng)預(yù)測(cè)

選取上述模型二反演獲取的最優(yōu)估值，基于震后黏彈性松弛模型，分別計(jì)算2015—2018年3年尺度、2018—2028年10年尺度、2028—2038年10年尺度、2038—2058年20年尺度水平向震后形變，如圖5所示。從圖中可以看出中場(chǎng)震后形變影響比較明顯，遠(yuǎn)場(chǎng)震后形變量級(jí)相對(duì)較小。不同時(shí)間尺度下測(cè)站累計(jì)最大/最小震后形變及年平均速率見(jiàn)表4。結(jié)果顯示，隨時(shí)間推移，震后形變年平均速率呈逐年減小趨勢(shì)，直至趨于穩(wěn)定。盡管由震后黏彈性松弛效應(yīng)引起的震后位移年變化量較小，但長(zhǎng)時(shí)間尺度累計(jì)下的震后形變絕對(duì)量依然較大，不可忽視。通過(guò)計(jì)算震后2018—2058年40年時(shí)間尺度累計(jì)位移，發(fā)現(xiàn)J413、H046、H037等測(cè)站震后位移最大能達(dá)到19 cm，遠(yuǎn)場(chǎng)測(cè)站H027、H028、H041等累計(jì)震后位移有5～8 cm。

圖4 3種不同滑動(dòng)模型計(jì)算的2010—2015年震后形變擬合值和模擬值對(duì)比及殘差分布Fig.4 Comparison of observed postseismic deformation to modeled displacements by using different coseismic slip models during 2010—2015.Right panels show residuals between observed and modeled postseismic displacements and coseismic slip models

圖5 模型計(jì)算給出的地表震后形變隨時(shí)間的變化Fig.5 Postseismic deformation changed over time according the model注：藍(lán)色：GNSS測(cè)站，品紅色：格網(wǎng)點(diǎn)。

5 結(jié) 論

本文通過(guò)對(duì)收集的1999—2015年間GNSS觀測(cè)資料進(jìn)行高精度數(shù)據(jù)處理和時(shí)間序列分析，并以2010—2015年震后形變?yōu)榧s束，基于震后黏彈性模型松弛進(jìn)行反演，得出如下結(jié)論：

(1) 通過(guò)時(shí)間序列分析，獲取最佳震后松弛時(shí)間為34 d，各測(cè)站具有明顯的震后形變趨勢(shì)(圖2)。給出研究區(qū)域內(nèi)109個(gè)測(cè)站在震后2010—2015年間累計(jì)震后形變(圖3)，其中水平向震后形變較大測(cè)站位于斷裂帶中部的上盤近、中場(chǎng)區(qū)域，形變量最大達(dá)7 cm。

(2) 不同區(qū)域測(cè)站震后形變的差異反映地下深部流變性質(zhì)的差異，分析發(fā)現(xiàn)龍門山斷裂上、下盤震后形變不在同一量級(jí)(圖3)，具有明顯的強(qiáng)不對(duì)稱性，表明在下地殼/上地幔這一層，下盤黏滯系數(shù)要明顯高于上盤。

(3) 基于震后黏彈性松弛模型，以龍門山斷裂上盤40個(gè)測(cè)站2010—2015年累計(jì)GNSS震后形變?yōu)榧s束，根據(jù)3種不同破裂模型，分別反演上盤中上地殼彈性層厚度和下地殼/上地幔黏滯系數(shù)最優(yōu)估值，并計(jì)算3種不同模型下2010—2015年震后形變模擬值和觀測(cè)值對(duì)比情況及殘差分布(圖4)。

(4) 3種不同模型計(jì)算的2010—2015年震后形變擬合值和模擬值均具有較好的一致性，對(duì)不同破裂模型計(jì)算結(jié)果進(jìn)行比較分析發(fā)現(xiàn)，不同模型對(duì)近場(chǎng)震后形變有一定影響，對(duì)遠(yuǎn)場(chǎng)形變幾乎無(wú)影響，除個(gè)別測(cè)站有明顯變化外，整體差異性不大。龍門山斷裂上盤岷江斷裂以西、鮮水河斷裂以東區(qū)域，測(cè)站擬合情況基本一致。岷江斷裂以東區(qū)域模擬值整體有偏大趨勢(shì)，尤其在岷江斷裂附近，表明岷江斷裂東西兩側(cè)地下流變性質(zhì)可能有一定差異。模型一考慮滑脫層的存在，J412、J413、H037等測(cè)站黏彈性模擬值被低估，可能含有少量震后余滑效應(yīng)，其中最大余滑量級(jí)約30%；另更加準(zhǔn)確的滑動(dòng)模型，在岷江斷裂附近模擬效果改善有限，除H034擬合較好外，其余站點(diǎn)殘差依然較大。

(5) 選取上述模型二反演的最優(yōu)估值進(jìn)行震后黏彈性松弛效應(yīng)進(jìn)行分析預(yù)測(cè)(圖5)，其中h=33 km，η=1.6×1019pa·s。汶川震后形變隨時(shí)間呈遞減趨勢(shì)，盡管年平均速率較小，但隨時(shí)間推移累計(jì)震后形變量依然較大，不容忽視，僅2018—2058年40年尺度累計(jì)震后位移最大能達(dá)到19 cm。

龍門山斷裂帶覆蓋地域較廣，汶川地震破裂情況比較復(fù)雜，盡管本文使用的GNSS觀測(cè)資料在數(shù)量和時(shí)間尺度上有限，且受地殼分層結(jié)構(gòu)局限，但根據(jù)反演得到的中上地殼彈性層厚度和上地殼/下地幔黏彈性層黏滯系數(shù)與已有研究結(jié)果基本一致，表明本文計(jì)算、模擬具有一定可靠性。通過(guò)計(jì)算分析發(fā)現(xiàn)，汶川地震在未來(lái)幾十年內(nèi)主要對(duì)龍門山斷裂帶上盤區(qū)域地殼形變影響較大，震后余滑局限在同震破裂100 km范圍內(nèi)分布。對(duì)汶川地震震后形變的持續(xù)深入研究，可以為研究龍門山斷裂帶運(yùn)動(dòng)學(xué)和區(qū)域地球動(dòng)力學(xué)過(guò)程及其演化提供重要依據(jù)。

謹(jǐn)以此文獻(xiàn)給汶川地震10周年紀(jì)念！

致謝：感謝中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)(CMONOC)提供的GNSS觀測(cè)數(shù)據(jù)。