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

劉志軍 譚 凱 王 琪 汪 雷 張 健 趙 斌 喬學軍

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

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

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

1 數據來源

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

2 震后粘彈性松弛

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

2.1 模型構建方法

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

表1 斷層模型的幾何特征

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

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

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

2.2 彈性層厚度選取

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

表2 三維粘彈模型分層介質參數

圖2 搜索龍門山斷裂帶彈性層厚度與中下地殼粘滯系數的最優估值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年的震后形變數據約束青藏高原東部的巖石圈流變結構，認為該時段的震后形變主要受粘彈性松弛效應控制。計算最優估值時所用的GNSS站點中有2個測站位于斷層附近，其余29個測站均距離斷層區域較遠。研究表明[2]，增加100 km以內的近場測站進行反演得到的粘滯系數值與去除近場測站單獨反演獲取的結果差異較小，因此本文的測站分布不會影響反演結果。比較4個彈性層模型的擬合誤差可知(圖2)，最優彈性層厚度為25 km，對應的最優粘滯系數為4×1018Pa·s，RMS值為6.01 mm。

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

2.3 震后粘彈性松弛位移場

基于青藏高原東部區域最優彈性層厚度和最優中下地殼粘滯系數的比較，本文選擇彈性層厚度為25 km的三維粘彈性模型，青藏高原東部區域中下地殼、上地幔和四川盆地粘彈介質均采用麥克斯韋體，青藏高原東部中下地殼和上地幔粘滯系數分別設置為4×1018Pa·s和1×1020Pa·s，四川盆地上地幔粘滯系數設置為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)可以看出，多數測站分布于龍門山斷裂帶上盤中遠場、鮮水河斷裂帶附近及岷江斷裂以東區域，少數測站分布于四川盆地區域，龍門山斷裂帶上盤區域站點的震后形變擬合值最大約為3.2 cm，最小約為2 mm，與GNSS站點震后形變觀測值最大約為3.9 cm、最小約為2 mm一致。結合圖3(b)可知，龍門山斷裂帶遠場區域的擬合效果比中場區域好，且距離斷裂區域越近，擬合效果越差，考慮到距離斷裂區域越近其站點位移量越大，從而導致擬合誤差也偏大，也可能由該區域深部結構的流變性質分布不同所致。鮮水河斷裂帶附近站點的擬合剩余形變場偏大，可能與鮮水河斷裂帶附近復雜的構造環境有關，也可能受該區域其他地震及震后影響。站點H043和H044位于龍門山斷裂帶下盤區域，考慮到站點位置處于相對穩固的四川盆地塊體，站點震后形變量約為4 mm，形變量偏小且擬合效果不明顯，表明龍門山斷裂帶上下盤存在明顯的不對稱性，反映青藏高原東部區域和四川盆地區域地下深部的流變結構性質存在很大差異。在岷江斷裂附近及其以東區域，站點形變擬合值比觀測值小，可能受該區域介質粘滯系數偏大的影響；從站點形變擬合值的方向來看，存在沿斷裂帶向北逆時針旋轉的運動趨勢，符合汶川地震右旋走滑的特征。

3 震后余滑

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

3.1 單一余滑機制

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

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

(1)

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

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

3.2 粘彈性松弛與余滑聯合機制

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

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

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

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

4 討 論

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

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

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

因此，本文采用2010～2015年中遠場區域站點數據研究震后粘彈性松弛效應，采用2008～2009年數據研究震后余滑和震后粘彈性松弛的短期作用特征。

4.2 粘彈性松弛結果分析

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

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

4.3 單一余滑與多機制聯合結果分析

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

5 結 語

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

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

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

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

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