汶川地震震后余滑和震后粘彈性松弛數(shù)值模擬分析

劉志軍 譚 凱 王 琪 汪 雷 張 健 趙 斌 喬學(xué)軍

1 中國地震局地震研究所，武漢市洪山側(cè)路40號，430071 2 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，武漢市魯磨路388號，430074

2008-05-12龍門山中央斷裂帶和山前斷裂上發(fā)生汶川MW7.9地震，研究汶川地震震后形變不僅對理解震后機制具有重要意義，而且對青藏高原東部地區(qū)的區(qū)域流變動力學(xué)研究具有重要價值[1]。在現(xiàn)有針對汶川地震的研究中，部分研究為早期1～2 a短期變化，部分研究為長期持續(xù)跟蹤分析[1-3]，缺乏將短期數(shù)據(jù)與長期數(shù)據(jù)進(jìn)行聯(lián)合分析的研究。

本文首先采用位置分布均勻、廣泛且時間范圍為2010～2015年的GNSS形變觀測數(shù)據(jù)，建立精細(xì)的三維粘彈性模型，基于有限元方法進(jìn)行震后粘彈性松弛模擬，獲得最佳的中下地殼粘滯系數(shù)，并研究單一粘彈性松弛機制對震后形變場的影響；然后采用2008～2009年的GNSS形變觀測數(shù)據(jù)，利用具有上述最佳中下地殼粘滯系數(shù)的三維粘彈性模型計算震后1 a內(nèi)的粘彈性松弛效應(yīng)。將觀測值減去粘彈性松弛模擬值的剩余形變量作為約束，建立粘彈性松弛和余滑的雙機制聯(lián)合模型，并與運動學(xué)反演獲得的單一震后余滑模型進(jìn)行比較，分析震后1 a內(nèi)粘彈性松弛與震后余滑2種機制對汶川地震震后形變的相對貢獻(xiàn)。

1 數(shù)據(jù)來源

本文收集的形變觀測資料包含52個GNSS觀測站數(shù)據(jù)，包括2008～2009年汶川地震震后形變監(jiān)測網(wǎng)中21個觀測站資料[4]和2008～2015年陸態(tài)網(wǎng)中31個觀測站資料[2]。目前震后形變理論主要包含震后余滑和震后粘彈性松弛，震后余滑主要在震后1～2 a內(nèi)存在，2 a后一般以震后粘彈性松弛為主。但汶川地震震后1 a內(nèi)的形變影響因素既包含震后余滑也包含震后粘彈性松弛，為了能在震后1 a內(nèi)的形變研究中將兩者區(qū)分開，本文首先進(jìn)行震后2～7 a的粘彈性松弛模擬研究，再進(jìn)行震后1 a內(nèi)的余滑研究。

2 震后粘彈性松弛

本文采用2010～2015年GNSS震后形變觀測數(shù)據(jù)，考慮到該期間的震后形變主要受粘彈性松弛機制影響，因此可進(jìn)行單一粘彈性松弛的有限元模擬計算，通過計算得到的震后形變場研究汶川地震震后2～7 a的形變特征。

2.1 模型構(gòu)建方法

首先建立汶川地震的斷層模型，可分為青川-南壩、南壩-北川、北川-映秀及彭灌區(qū)域4個部分，其中前3個部分區(qū)域均為鏟狀曲面，第4部分區(qū)域為光滑平面，且彭灌區(qū)域斷層平面與北川-映秀區(qū)域鏟狀斷層在深度約20 km處相交。斷層模型的幾何特征見表1。

表1 斷層模型的幾何特征

根據(jù)地質(zhì)資料及斷層破裂位置和幾何特征建立長和寬均為3 000 km、深度為400 km的模型塊體(圖1(a))，基于汶川斷層模型對塊體進(jìn)行切割，在青藏高原東部區(qū)域設(shè)置彈性上地殼、粘彈性下地殼和上地幔，在四川盆地區(qū)域設(shè)置彈性上地殼和粘彈性上地幔。對網(wǎng)格進(jìn)行劃分，斷層附近的網(wǎng)格比遠(yuǎn)場網(wǎng)格密集，既可滿足精度需求又可提高計算效率，其中斷層區(qū)域網(wǎng)格劃分單元最小為3 km，邊界網(wǎng)格單元最大為40 km，共包含334 651個單元和59 535個節(jié)點。保持模型上表面自由無應(yīng)力，固定底部和四周的節(jié)點，模型距離斷層邊界位置足夠遠(yuǎn)以減少固定條件引起的影響。

為了更準(zhǔn)確地模擬震后形變場，本文根據(jù)文獻(xiàn)[5]中的同震滑動分布模型，采用同震形變場數(shù)據(jù)，利用Okada 彈性半空間位錯模型計算格林函數(shù)G，然后將反演的同震滑動分布模型作為計算粘彈性松弛的源斷層模型(圖1(b))。

圖1 汶川地震三維有限元模型Fig.1 Three-dimensional finite element model of the Wenchuan earthquake

2.2 彈性層厚度選取

在確定地震的破裂源后，地球表面的震后形變將取決于深度的流變結(jié)構(gòu)。分析表明，四川盆地的巖石圈結(jié)構(gòu)和力學(xué)性質(zhì)與青藏高原東部存在很大差異。參考文獻(xiàn)[6]中的建模分層方法，本文共建立4個三維有限元模型，設(shè)置青藏高原東部塊體的彈性層分別為25 km、30 km、33 km、35 km，粘性塊體介質(zhì)主要采用麥克斯韋體，粘滯系數(shù)搜索范圍為1017～1019Pa·s。根據(jù)介質(zhì)的差異性將上地殼、中下地殼和上地幔賦予不同的介質(zhì)參數(shù)(表2)，其中各層參數(shù)主要參考文獻(xiàn)[2]。利用有限元軟件Pylith[7]模擬計算因粘彈性松弛引發(fā)的汶川地震震后地表位移，通過網(wǎng)格搜索方法分別對4個模型的最佳粘滯系數(shù)進(jìn)行搜索，并計算擬合誤差RMS，擬合誤差為最小值時對應(yīng)的彈性層厚度即為本文青藏高原塊體彈性層的最優(yōu)估值(圖2)。

表2 三維粘彈模型分層介質(zhì)參數(shù)

圖2 搜索龍門山斷裂帶彈性層厚度與中下地殼粘滯系數(shù)的最優(yōu)估值Fig.2 Search for the optimal estimation of the elastic layer thickness and the viscosity coefficient of the middle and lower crust in the Longmenshan fault zone

本文利用2010～2015年的震后形變數(shù)據(jù)約束青藏高原東部的巖石圈流變結(jié)構(gòu)，認(rèn)為該時段的震后形變主要受粘彈性松弛效應(yīng)控制。計算最優(yōu)估值時所用的GNSS站點中有2個測站位于斷層附近，其余29個測站均距離斷層區(qū)域較遠(yuǎn)。研究表明[2]，增加100 km以內(nèi)的近場測站進(jìn)行反演得到的粘滯系數(shù)值與去除近場測站單獨反演獲取的結(jié)果差異較小，因此本文的測站分布不會影響反演結(jié)果。比較4個彈性層模型的擬合誤差可知(圖2)，最優(yōu)彈性層厚度為25 km，對應(yīng)的最優(yōu)粘滯系數(shù)為4×1018Pa·s，RMS值為6.01 mm。

根據(jù)采用有限元方法建立的彈性層厚度為35 km的三維模型[1]計算單一粘彈性松弛，結(jié)果表明，青藏高原東部中下地殼最優(yōu)粘滯系數(shù)為1×1018Pa·s，然后建立聯(lián)合模型計算得到青藏高原東部中下地殼對應(yīng)的最優(yōu)粘滯系數(shù)為2×1018Pa·s。溫?fù)P茂等[8]采用將同震和震后形變進(jìn)行聯(lián)合的方法，獲得中下地殼粘滯系數(shù)下限為2.0×1018Pa·s；Huang等[6]采用震后1 a的GNSS觀測數(shù)據(jù)及震后1.5 a的InSAR觀測資料，基于有限元方法構(gòu)建龍門山斷裂帶兩側(cè)粘滯系數(shù)差異的巖石圈流變結(jié)構(gòu)三維模型，得到四川盆地粘滯系數(shù)不低于1.0×1020Pa·s，而川西高原下地殼瞬態(tài)和穩(wěn)態(tài)粘滯系數(shù)分別為4.4×1017Pa·s和1.0×1018Pa·s。本文建立的4個三維有限元模型計算粘彈性松弛對應(yīng)的最佳粘滯系數(shù)為2×1018～4×1018Pa·s，與前人的研究結(jié)果基本一致。

2.3 震后粘彈性松弛位移場

基于青藏高原東部區(qū)域最優(yōu)彈性層厚度和最優(yōu)中下地殼粘滯系數(shù)的比較，本文選擇彈性層厚度為25 km的三維粘彈性模型，青藏高原東部區(qū)域中下地殼、上地幔和四川盆地粘彈介質(zhì)均采用麥克斯韋體，青藏高原東部中下地殼和上地幔粘滯系數(shù)分別設(shè)置為4×1018Pa·s和1×1020Pa·s，四川盆地上地幔粘滯系數(shù)設(shè)置為1×1020Pa·s。將2010～2015年的GNSS震后形變觀測值作為模擬值，計算震后2～7 a粘彈性松弛影響的震后形變位移(圖3)。

圖3 2010～2015年震后形變模擬值和觀測值對比及剩余形變場分布Fig.3 Comparison of the observed and simulated values of post-seismic deformation and distribution of the residual deformation filed during 2010-2015

從圖3(a)可以看出，多數(shù)測站分布于龍門山斷裂帶上盤中遠(yuǎn)場、鮮水河斷裂帶附近及岷江斷裂以東區(qū)域，少數(shù)測站分布于四川盆地區(qū)域，龍門山斷裂帶上盤區(qū)域站點的震后形變擬合值最大約為3.2 cm，最小約為2 mm，與GNSS站點震后形變觀測值最大約為3.9 cm、最小約為2 mm一致。結(jié)合圖3(b)可知，龍門山斷裂帶遠(yuǎn)場區(qū)域的擬合效果比中場區(qū)域好，且距離斷裂區(qū)域越近，擬合效果越差，考慮到距離斷裂區(qū)域越近其站點位移量越大，從而導(dǎo)致擬合誤差也偏大，也可能由該區(qū)域深部結(jié)構(gòu)的流變性質(zhì)分布不同所致。鮮水河斷裂帶附近站點的擬合剩余形變場偏大，可能與鮮水河斷裂帶附近復(fù)雜的構(gòu)造環(huán)境有關(guān)，也可能受該區(qū)域其他地震及震后影響。站點H043和H044位于龍門山斷裂帶下盤區(qū)域，考慮到站點位置處于相對穩(wěn)固的四川盆地塊體，站點震后形變量約為4 mm，形變量偏小且擬合效果不明顯，表明龍門山斷裂帶上下盤存在明顯的不對稱性，反映青藏高原東部區(qū)域和四川盆地區(qū)域地下深部的流變結(jié)構(gòu)性質(zhì)存在很大差異。在岷江斷裂附近及其以東區(qū)域，站點形變擬合值比觀測值小，可能受該區(qū)域介質(zhì)粘滯系數(shù)偏大的影響；從站點形變擬合值的方向來看，存在沿斷裂帶向北逆時針旋轉(zhuǎn)的運動趨勢，符合汶川地震右旋走滑的特征。

3 震后余滑

研究表明[3,6,9]，余滑模型能夠很好地擬合近場觀測數(shù)據(jù)，但在遠(yuǎn)場區(qū)域擬合效果不佳；粘彈性松弛模型則相反，能夠較好地擬合中遠(yuǎn)場觀測數(shù)據(jù)，卻無法擬合近場數(shù)據(jù)。因此對于汶川地震，單一機制無法很好地解釋觀測到的震后形變場，本文將聯(lián)合震后余滑與粘彈性松弛，建立合理的雙機制聯(lián)合模型進(jìn)行反演。

3.1 單一余滑機制

震后余滑是地震后同震應(yīng)力變化引起無震滑移的過程，可通過對GNSS測量的累積震后位移進(jìn)行反演來分析破裂斷層表面的余滑分布。假設(shè)可通過彈性半空間位錯模型來描述余滑，反演的目標(biāo)函數(shù)可表示為：

F(s,β)=‖W(Gs-d)‖+β‖?2s‖

(1)

式中，s為估計的滑動分量；W為根據(jù)觀測不確定性構(gòu)造的權(quán)重矩陣；d為大地測量數(shù)據(jù)矢量；G為Green函數(shù)矩陣；?2為用于平滑滑動模型的Laplacian算子的有限差分，且β為調(diào)整平滑度的因子。通過使用有界可變最小二乘(BVLS)算法獲得最佳模型解，β最佳值應(yīng)在模型粗糙度和數(shù)據(jù)失配之間進(jìn)行權(quán)衡[10]。

本文單一余滑模型采用運動學(xué)反演方法，認(rèn)為汶川地震震后1 a內(nèi)的形變均由余滑作用引起，故將汶川地震震后1 a 的GNSS形變觀測數(shù)據(jù)作為輸入值。利用反演函數(shù)計算震后1 a內(nèi)破裂斷層面的余滑分布(圖4(d))，得到地面站點1 a內(nèi)的震后形變值，進(jìn)而分析汶川地震在單一余滑機制下震后1 a內(nèi)的形變特征(圖4(a))。

3.2 粘彈性松弛與余滑聯(lián)合機制

在研究粘彈性松弛時，利用汶川地震震后2～7 a的形變數(shù)據(jù)測得青藏高原東部彈性層最優(yōu)厚度為25 km，龍門山斷裂帶上盤中下地殼最佳粘滯系數(shù)為4.0×1018Pa·s，此時擬合誤差最小，RMS值為6.01 mm。這些結(jié)果與前人的研究結(jié)論相近[1,11]，故可認(rèn)為得到的中下地殼粘滯系數(shù)具有可靠性。假設(shè)汶川地震震后7 a內(nèi)龍門山斷裂帶中下地殼的粘滯系數(shù)不變且為4.0×1018Pa·s，選擇彈性層厚度為25 km的三維有限元模型對震后1 a內(nèi)的粘彈性松弛進(jìn)行正演計算，再將震后1 a的觀測數(shù)據(jù)減去粘彈性松弛的正演結(jié)果，得到剩余形變場。在只考慮粘彈性松弛與余滑2種機制影響的條件下，去除粘彈性松弛形變后的剩余形變場可認(rèn)為均由余滑作用引起，將剩余形變場作為輸入值，采用上述反演程序得到破裂斷層面的余滑分布(圖4(c))及由余滑影響的地面震后形變結(jié)果。將正演粘彈性松弛的地表形變結(jié)果與剩余形變場反演得到的余滑影響下的地表形變結(jié)果相加，得到2種機制的聯(lián)合結(jié)果，然后與震后1 a內(nèi)GNSS觀測結(jié)果進(jìn)行對比，研究雙機制下汶川地震震后1 a內(nèi)的形變場特征(圖4(a))。

(a)黑色箭頭為震后形變觀測值，紅色箭頭為粘彈性松弛與余滑聯(lián)合機制下的震后形變擬合值，藍(lán)色箭頭為單一余滑機制下的震后形變擬合值；(b) 紅色箭頭為粘彈性松弛與余滑聯(lián)合機制下的剩余形變量，藍(lán)色箭頭為單一余滑機制下的剩余形變量；(c)為通過去除粘彈性松弛后的剩余形變場反演的斷層面上余滑分布結(jié)果；(d)為通過單一余滑機制下的震后形變觀測值反演的破裂斷層面上余滑分布結(jié)果圖4 兩種方法反演的震后形變模擬值與觀測值比較Fig.4 Comparison of post-seismic deformation simulated values and observed values of two methods

根據(jù)圖4(a)所示，將所用的52個地表站點劃分為3個區(qū)域，距離斷層破裂區(qū)域最遠(yuǎn)的站點擬合效果最差，距離較近的站點擬合效果較好，位于汶川地震破裂面附近的站點擬合結(jié)果與震后形變觀測值基本一致，2種擬合方法均表現(xiàn)出由近到遠(yuǎn)擬合誤差越來越大的特點。H027、H028、H029等遠(yuǎn)場站點的擬合效果較差，這可能是因為斷層面滑脫層的水平寬度與延伸深度較小，使龍門山斷裂帶上盤區(qū)域遠(yuǎn)場位移受到較大影響[1]；而在鮮水河斷裂帶區(qū)域的站點，如H053、H054、H056、H062等站點的擬合值普遍偏小，這可能是鮮水河斷裂帶區(qū)域地殼結(jié)構(gòu)復(fù)雜、介質(zhì)屬性與附近區(qū)域存在較大差異所致。從單一余滑機制來看，余滑存在于斷層面上，其影響范圍一般覆蓋近場及中場區(qū)域，對于遠(yuǎn)場站點的形變影響較小，因此出現(xiàn)藍(lán)色擬合值由近到遠(yuǎn)逐漸減小的情況。從粘彈性松弛與余滑聯(lián)合機制來看，雖然也存在擬合值漸變的情況，但從中遠(yuǎn)場擬合效果及圖4(b)可以看出，聯(lián)合機制的結(jié)果比單一余滑機制下的擬合值偏大，且更接近震后形變觀測值，表明多機制聯(lián)合比單一余滑機制能更好地解釋汶川地震震后1 a內(nèi)的形變特征。

由圖4(c)和4(d)可以看出，2種擬合方法的余滑結(jié)果基本一致，且主要分布在同震過程中未滑動或滑動量較小的區(qū)域，包括青川東北部斷坡區(qū)域、南壩-北川段斷坡區(qū)域、虹口-映秀深部滑脫層區(qū)域及映秀西南部斷坡區(qū)域等，其中單一余滑機制反演的余滑模型最大滑動量約為0.62 m，雙機制聯(lián)合反演的余滑模型最大滑動量約為0.54 m。從圖4還可以看出，斷層淺部和深部均存在余滑分布，說明斷層的震后余滑不僅發(fā)生在斷層淺部，同時也可能發(fā)生在深部。青川東北部斷坡區(qū)域可能是由斷層破裂后東北部區(qū)域向北走滑形成；北川斷坡區(qū)域可能是由同震破裂后存在的余震作用引起；虹口-映秀深部滑脫層區(qū)域可能是由于地震發(fā)生后其逆沖趨勢依然存在，導(dǎo)致深部繼續(xù)出現(xiàn)余滑運動。上述分析表明，余滑形變繼承了同震形變的特征，從西南部逆沖過渡到中部逆沖兼走滑，東北部為純走滑運動。

4 討 論

4.1 震后粘彈性松弛與震后余滑時間劃分分析

根據(jù)文獻(xiàn)[2]計算震后粘彈性松弛效應(yīng)的方法，本文采用震后2～7 a的形變數(shù)據(jù)研究該時間段的震后粘彈性松弛效應(yīng)，同時采用雙機制聯(lián)合方法對震后2 a內(nèi)的粘彈性松弛和余滑效應(yīng)的時間序列進(jìn)行模擬分析。由圖5可知，H035站點距離破裂面較近，震后粘彈性松弛形變極小，大部分為震后余滑成分；SEEG站點距離破裂面較遠(yuǎn)，震后粘彈性松弛所占分量較大；WARI站點距離破裂面最遠(yuǎn)，震后粘彈性松弛成分甚至比震后余滑多。從增長趨勢可以看出，震后1 a內(nèi)余滑比重較大且增長明顯，震后1～2 a內(nèi)余滑增長量逐漸減少并趨于平緩；而震后粘彈性松弛一直呈遞增趨勢，尤其在中遠(yuǎn)場區(qū)域增長趨勢明顯。各站點粘彈性松弛殘差可近似認(rèn)為由余滑引起，根據(jù)震后形變時間序列估計，在震后2～7 a內(nèi)中場區(qū)域(以SEEG站點為例)余滑效應(yīng)占比約36%，遠(yuǎn)場區(qū)域(以TAGO站點為例)余滑效應(yīng)占比約25%。綜合圖3和5可知，2010～2015年中遠(yuǎn)場區(qū)域震后形變以粘彈性松弛效應(yīng)為主，因此利用震后2～7 a的形變數(shù)據(jù)研究震后粘彈性松弛效應(yīng)具有合理性。由于在研究震后粘彈性松弛效應(yīng)時采用的測站均分布于中遠(yuǎn)場區(qū)域，故對近場區(qū)域不作考慮。

圖5 震后2 a內(nèi)形變時間序列觀測值和模擬值Fig.5 Observed and simulated values of deformation time series after the earthquake within 2 years

因此，本文采用2010～2015年中遠(yuǎn)場區(qū)域站點數(shù)據(jù)研究震后粘彈性松弛效應(yīng)，采用2008～2009年數(shù)據(jù)研究震后余滑和震后粘彈性松弛的短期作用特征。

4.2 粘彈性松弛結(jié)果分析

由圖2可知，彈性層厚度為25 km時擬合誤差最小，為6.01 mm。根據(jù)crust2.0地殼模型可知，北川-青川地區(qū)上地殼彈性層厚度為17 km，映秀地區(qū)為22 km。夏婷婷[12]通過建模得出龍門山斷裂帶上地殼厚度為22 km，且自北向南厚度逐漸變薄，說明龍門山斷裂帶上盤彈性層厚度并非均勻分布，北川-青川地區(qū)較淺，映秀地區(qū)較深。考慮到改變龍門山斷裂帶上盤彈性層厚度會影響中下地殼最優(yōu)粘滯系數(shù)估計和震后粘彈性松弛效應(yīng)，從而影響震后形變特征，建立非均勻厚度的彈性層具有重要意義。

利用有限元方法對單一粘彈性松弛進(jìn)行震后2～7 a的正演模擬，結(jié)果表明，最大形變量為32 mm，最小為2 mm，與震后形變觀測值較為接近。擬合結(jié)果與觀測值存在誤差，可能是受計算誤差的影響，也可能是由三維建模中橫向、縱向的分層粗糙導(dǎo)致，還可能與站點位置區(qū)域的地殼構(gòu)造影響有關(guān)。下一步可根據(jù)汶川地震的層析成像結(jié)果建立更精細(xì)的三維模型，在斷層兩側(cè)施加不均勻分布的介質(zhì)，進(jìn)而完善模擬效果。

4.3 單一余滑與多機制聯(lián)合結(jié)果分析

以往研究中，直接由地表震后位移通過線性反演得到震后余滑是解釋觀測結(jié)果最方便、最有效的方法，但該方法缺乏物理基礎(chǔ)，需要設(shè)置大量的自由參數(shù)。此外，由于忽略了粘彈性松弛的貢獻(xiàn)，獲得的余滑在大小和深度上往往偏大[1]。通過對比單一余滑機制與雙機制聯(lián)合的模擬結(jié)果及剩余形變量(圖4(a)、4(b))可知，斷層面附近測站的模擬結(jié)果均與觀測值高度擬合，說明在斷層面附近的震后形變粘彈性松弛作用較小，余滑起主導(dǎo)作用；對于中遠(yuǎn)場區(qū)域，雙機制聯(lián)合的模擬結(jié)果均比單一余滑機制的模擬結(jié)果更接近震后形變觀測值，說明在中遠(yuǎn)場區(qū)域粘彈性松弛起到一定作用。從整體上看，多機制聯(lián)合的擬合效果比單一余滑機制的擬合效果好，單一的形變機制不能很好地解釋汶川地震震后形變場，多機制聯(lián)合的方法能更好地說明汶川地震震后的形變特征。

5 結(jié) 語

本文采用位置分布均勻、廣泛且時間范圍分別為2010～2015年和2008～2009年的31個GNSS形變觀測資料及21個GNSS形變觀測資料進(jìn)行震后形變分析研究，得到以下結(jié)論：

1)龍門山斷裂帶兩側(cè)的站點震后形變具有明顯的不對稱性，下盤形變值小于上盤，表明斷層兩側(cè)上下盤的深部流變性質(zhì)存在巨大差異。

2)震后2～7 a粘彈性松弛模擬結(jié)果表明，龍門山斷裂帶上盤最佳彈性層厚度為25 km，中下地殼最佳粘滯系數(shù)為4.0×1018Pa·s，此時擬合誤差RMS為6.01 mm。

3)震后2～7 a粘彈性松弛模型在中近場區(qū)域擬合誤差偏大，在遠(yuǎn)場區(qū)域擬合效果較好。

4)震后1 a內(nèi)的形變影響中余滑占主導(dǎo)作用；雙機制聯(lián)合方法得到的震后形變結(jié)果比單一余滑機制更接近實際觀測值，粘彈性松弛和余滑的聯(lián)合機制能更有效地解釋震后形變特征。