張 健 趙 斌 王東振 王海斌 劉志軍

1 中國地震局地震大地測量重點實驗室，武漢市洪山側(cè)路40號，430071

巖石對地球內(nèi)部力學(xué)作用的變形響應(yīng)取決于巖石的流變屬性，而巖石圈及軟流圈的粘滯系數(shù)是理解板塊以及塊體邊界帶演化過程的重要參數(shù)[1]。地震可以看作是突然發(fā)生的自然巖石力學(xué)實驗，其引起的應(yīng)力變化導(dǎo)致下地殼及上地幔粘彈性弛豫，從而使地面發(fā)生長期形變[2-3]。目前，采用大地測量技術(shù)獲取的震后形變已被廣泛用于探測地球巖石圈的流變結(jié)構(gòu)[4-5]。

2015-04-25青藏高原南緣的主喜馬拉雅逆沖斷裂帶(main Himalayan thrust, MHT)發(fā)生MW7.8地震。該地震還伴隨7級以上強余震，其中最大余震發(fā)生在主震震中東南約150 km處，矩震級為MW7.3[6-7]。布設(shè)在尼泊爾及藏南的GPS臺站在震后記錄到顯著形變，為探測青藏高原南緣下地殼和上地幔流變特性提供了條件，也為進一步認識喜馬拉雅逆沖斷層的地震演化和青藏高原隆升機制提供了重要參考資料。

尼泊爾MW7.8地震發(fā)生后，大量學(xué)者利用大地測量資料對該區(qū)域的流變結(jié)構(gòu)進行研究[3,8-9]，但使用的GPS觀測數(shù)據(jù)局限于尼泊爾附近，缺少距離震中200 km以外的青藏高原南部地區(qū)的數(shù)據(jù)；且假設(shè)的印度板塊和青藏高原板塊垂直邊界過于簡單，因此模擬的下地殼及上地幔粘滯系數(shù)具有較大的不確定性。另外，尼泊爾-喜馬拉雅地區(qū)的大地電磁、地震波層析成像、地?zé)峤Y(jié)構(gòu)、有效巖石圈厚度等相關(guān)研究顯示，喜馬拉雅地區(qū)兩側(cè)的結(jié)構(gòu)具有顯著的橫向不均勻性[10]，而這種特性會對傳統(tǒng)解析或半解析方法的結(jié)果造成很大影響[8]。有限元模擬方法在解決該類問題上具有較大優(yōu)勢。本文基于更多藏南遠場GPS站點震后5 a的觀測數(shù)據(jù)，采用有限元方法模擬尼泊爾地震的震后粘彈性松弛效應(yīng)，并通過垂直邊界和俯沖帶結(jié)構(gòu)邊界來約束區(qū)域下地殼和上地幔的流變特性，該研究對評估喜馬拉雅地震帶的地震災(zāi)害、認識地震機理和巖石圈流變特性具有重大意義。

1 GPS數(shù)據(jù)與處理

1.1 數(shù)據(jù)處理

本文共收集51個GPS站點的觀測數(shù)據(jù)，其中包含近期新采集的距離破裂區(qū)150～250 km的13個位于藏南地區(qū)的連續(xù)站數(shù)據(jù)[11]和7個新增的流動站數(shù)據(jù)。采用GAMIT/GLOBK 10.7版本對所有的原始RINEX數(shù)據(jù)進行統(tǒng)一處理。首先利用IGS發(fā)布的精密星歷估計測站坐標(biāo)、衛(wèi)星軌道、大氣延遲等參數(shù)，得到無基準(zhǔn)的自由網(wǎng)解；然后采用卡爾曼濾波方法計算區(qū)域單日松弛解與全球單日松弛解的組合解，并通過50個IGS追蹤站將組合解歸算至ITRF2014參考框架。由此得到的坐標(biāo)時間序列中仍含有因長期構(gòu)造運動、季節(jié)性的周年、半周年周期變化、同震或儀器變更等引起的階躍信號[3]，表達式為：

X(t)=A+Bt+C1sin(2πt)+C2cos(2πt)+

(1)

式中，A為震后形變；B為線性運動速率；C1、C2、C3、C4為描述季節(jié)性的周年、半周年周期信號的系數(shù)；Di為地震同震或儀器變更(如替換天線)等原因引起的階躍信號；H(t)為階躍函數(shù)。在提取震后形變時應(yīng)將上述信號從原始時間序列中去除。

對于在尼泊爾地震之前已經(jīng)開始觀測的GPS站，本文直接估算其長期運動速率。部分GPS站由于損壞或在震后才建立，無震前觀測數(shù)據(jù)，無法直接估計其長期運動速率，因此采用內(nèi)插方法進行計算并從震后形變中去除[12]。此外，采用Zhao等[3]的方法改正2015-05-12 MW7.3余震引起的階躍項以及季節(jié)性的周年、半周年周期變化。經(jīng)過上述改正，獲得每個GPS站點的震后位移時間序列。

1.2 震后形變

由于不同測站觀測資料時長不同，為得到GPS站點震后5 a的累積位移，對改正后的時間序列采用對數(shù)函數(shù)進行擬合：

(2)

式中，u(t)為t時刻累積觀測位移；a1為常數(shù)；a2為震后形變幅度；τ為震后松弛時間。圖1為本文所采用的GPS站點累積5 a的震后位移，其中藍色箭頭代表水平位移，顏色變化代表垂直位移。從圖中可以看出，尼泊爾地震引起地面向南-西南方向運動，震后形變主要集中在尼泊爾北部和中尼邊界區(qū)域，并向北衰減，靠近主前緣斷裂帶的GPS站點的位移較小。CHLM站的水平位移最大，約為121 mm，NAST站的垂直位移最大，為83 mm。

圖1 尼泊爾地震構(gòu)造背景及震后5 a累積位移

2 模型與方法

2.1 有限元模型

參考Elliott等[13]提出的“雙斷坡”斷層模型，以斷層面中心點在地面的投影為笛卡爾坐標(biāo)系原點，建立三維粘彈性有限元模型(圖2)。模型南北方向長度為2 000 km，東西方向長度為1 500 km，厚度為800 km。為更好地計算斷層破裂附近的應(yīng)力場和位移，首先將斷層面劃分為邊長約5 km的三角形，然后設(shè)置格網(wǎng)尺寸隨距斷層面的距離以1∶1的比例因子逐漸增大進行整體劃分，最終共劃分445 210個四面體單元。模型共包含95 130個節(jié)點，最小格網(wǎng)尺寸為5 km，最大格網(wǎng)尺寸為158 km。

圖2 有限元模型及格網(wǎng)劃分

保持模型側(cè)邊界垂直及下邊界水平，上邊界可進行調(diào)整，采用文獻[14]中P波和S波速度結(jié)構(gòu)，并設(shè)置上地殼、下地殼和上地幔的密度分別為2 700 kg/m3、3 300 kg/m3和3 800 kg/m3。

2.2 巖石圈流變結(jié)構(gòu)配置

本文采用2種不同的印度板塊與青藏高原板塊邊界(BL)來構(gòu)建模型：一種為垂直邊界模型(模型1，圖3(a))，另一種為印度彈性俯沖板片模型(模型2，圖3(b))。對于模型1，假設(shè)印度板塊彈性地殼厚度為50 km，其下是粘滯系數(shù)為1×1020Pa·s 的Maxwell粘彈性上地幔；而在青藏高原南部，根據(jù)該地區(qū)地震震源深度分布，取彈性上地殼厚度為30 km，下地殼由30 km延伸至青藏高原莫霍面平均深度70 km[15]，其下為Maxwell粘彈性上地幔，粘滯性系數(shù)固定為1×1020Pa·s。對于模型2，根據(jù)Nábělek等[16]的研究成果，設(shè)定印度俯沖板片在MFT以南厚度為50 km，而在MFT以北，隨著向青藏高原下地殼俯沖，其厚度逐漸減小為25 km，印度板塊和青藏高原上地幔的邊界由俯沖板片下邊界擴展劃分得到，其他參數(shù)與模型1保持一致。通過上述方法可構(gòu)建考慮喜馬拉雅區(qū)域尼泊爾中部的橫向不均勻性的巖石圈流變結(jié)構(gòu)模型。

圖3 喜馬拉雅區(qū)域尼泊爾中部的橫向不均勻流變結(jié)構(gòu)

2.3 模擬方法

采用Pylith有限元分析軟件模擬尼泊爾地震震后粘彈性松弛效應(yīng)[17]，從而約束青藏高原南緣的流變參數(shù)。將Pylith中GenMaxwell體轉(zhuǎn)換為雙粘性Burgers體來表示青藏高原下地殼的流變屬性[18]，Burgers體可同時捕捉到震后早期形變的快速衰減和長期線性變化。本文定義穩(wěn)態(tài)Maxwell粘滯系數(shù)ηm為瞬態(tài)Kelvin粘滯系數(shù)ηk的10倍[3]?？紤]到主震和MW7.3最大余震的影響，為更可靠地模擬震后形變，本文基于“雙斷坡”斷層模型重新反演主震和MW7.3余震的同震滑動。通過最鄰近內(nèi)插方法將同震滑動施加到斷層面上以驅(qū)動粘彈性松弛的發(fā)生。巖石對粘彈性松弛效應(yīng)的響應(yīng)由區(qū)域流變結(jié)構(gòu)決定，通過評估模擬值和觀測值之間的加權(quán)均方根(WRMS)來搜索最優(yōu)的流變參數(shù)，計算公式為：

(3)

式中，N為GPS點數(shù)量;T為觀測時間; obs和cal分別為觀測值和模擬值;σ為GPS觀測值的不確定度。

以往研究認為，尼泊爾地震的震后形變受多種機制影響(震后余滑、粘彈性松弛)，而粘彈性松弛主要控制遠場位移及整體震后形變的長期效應(yīng)[3,9]。因此，為更準(zhǔn)確地探測青藏高原下地殼的粘滯系數(shù)，本文采用2種方案進行約束并相互驗證：1)假設(shè)余滑在震后2 a內(nèi)完全衰減[9]，采用所有GPS站點震后2～5 a的觀測數(shù)據(jù)進行約束；2)采用距同震破裂區(qū)200 km以外的9個遠場GPS點(XZZB, XZAR、CANA、QIRE、YANI、J030、JZLS、J038、J338)震后5 a的觀測數(shù)據(jù)進行擬合約束，由于遠場點J337的觀測位移比周圍其他點偏大，從而使粘滯系數(shù)減小，因此計算時不采用該點數(shù)據(jù)。

3 結(jié)果分析

3.1 模型1結(jié)果

采用模型1模擬粘彈性松弛并進行最優(yōu)化分析，該過程中使用2個參數(shù)：青藏高原南緣下地殼的穩(wěn)態(tài)Maxwell粘滯系數(shù)ηm(變化范圍為1×1017～1×1019Pa·s)以及MFT到印度板塊與青藏高原板塊垂直邊界的距離D(變化范圍為120～180 km)。為提高計算效率，首先設(shè)置較大的參數(shù)間隔進行搜索，然后在誤差較小的范圍內(nèi)搜索最優(yōu)的巖石圈流變參數(shù)。圖4(a)為采用所有GPS站點震后2～5 a的觀測數(shù)據(jù)約束的最優(yōu)參數(shù)D和ηm，圖4(b)為采用遠場GPS站點震后5 a的數(shù)據(jù)約束的結(jié)果，2種方法約束的參數(shù)結(jié)果相等：D為130 km，ηm和ηk分別為3×1018Pa·s和3×1017Pa·s。圖4(c)為觀測的GPS站點震后5 a累積位移(藍色箭頭)和模擬的地面節(jié)點累積位移(綠色箭頭)，圖中顏色變化代表模擬的垂直位移。從圖中可以看出，模擬結(jié)果整體顯示為南-西南的水平地面運動模式，與觀測數(shù)據(jù)基本一致，并且可以很好地解釋遠場點的觀測位移，但中近場模擬的水平位移嚴(yán)重偏小，這可能是因為未考慮震后余滑的影響。

圖4 垂直邊界模型的最優(yōu)網(wǎng)格搜索結(jié)果和震后5 a累積位移

3.2 模型2結(jié)果

對于模型2，由于印度板塊與青藏高原板塊的邊界由印度俯沖板片的下邊界來確定，因此只約束青藏高原下地殼穩(wěn)態(tài)粘滯系數(shù)ηm。首先完成ηm變化范圍在1×1017～1×1019Pa·s之間的粘彈性松弛模擬，然后采用與§3.1相同的方法進行最優(yōu)化網(wǎng)格搜索，2種方案的搜索結(jié)果如圖5(a)和5(b)所示。圖5(a)為采用所有GPS站點震后2～5 a的觀測數(shù)據(jù)約束的結(jié)果，圖5(b)為采用遠場GPS站點震后5 a的數(shù)據(jù)約束的結(jié)果，紅色圓點代表最優(yōu)參數(shù)，2種方法約束的青藏高原下地殼穩(wěn)態(tài)粘滯系數(shù)均為1×1018Pa·s，對應(yīng)的瞬態(tài)粘滯系數(shù)為1×1017Pa·s，比模型1約束的結(jié)果更小。圖5(c)為震后5 a的累積位移，藍色箭頭代表GPS觀測位移，綠色箭頭為模擬的地面格網(wǎng)節(jié)點位移，顏色變化代表模擬的垂直位移。從圖中可以看出，其與觀測數(shù)據(jù)和模型1具有非常接近的水平運動模式，可合理地解釋遠場GPS站點的水平觀測位移，但中近場的模擬值遠小于觀測值。

圖5 印度俯沖板片模型的最優(yōu)網(wǎng)格搜索結(jié)果和震后5 a累積位移

3.3 不同模型結(jié)果比較

通過上述分析可知，2種邊界模型均可產(chǎn)生南-西南水平方向運動模式，與觀測的地面位移具有很好的一致性，但模型2約束的青藏高原下地殼粘滯系數(shù)略小于模型1，原因有2個方面：一是D和ηm之間的折中關(guān)系會導(dǎo)致印度板塊與青藏高原板塊的垂直邊界更偏南，使得青藏高原下地殼向南延伸，范圍變大；二是印度彈性俯沖板片向青藏高原下地殼推擠會降低粘彈性層的厚度。關(guān)于垂直形變(圖4(c)和5(c)中顏色變化)，2種模型的模擬結(jié)果均顯示了中尼邊界附近及以北地區(qū)的隆升現(xiàn)象和雅魯藏布江縫合帶周圍的沉降現(xiàn)象，但大小具有很大差異，模型1的垂直位移峰值約為80 mm，而模型2約為40 mm，模型2的垂直位移與GPS觀測的垂直位移更為一致。而對于中尼邊界以南的區(qū)域，粘彈性松弛模擬的垂直位移遠小于觀測值，這可能與其他機制有關(guān)(如震后余滑)。

模型1和模型2采用所有GPS站點震后2～5 a的水平位移擬合計算的WRMS分別為3.42 mm和6.78 mm，而采用遠場點震后5 a的水平位移擬合計算的WRMS分別為2.13 mm和4.08 mm。相比之下，模型1的WRMS較模型2更小，雖然可以說明模型1能更好地解釋觀測的水平位移，但這主要是由于模型1中考慮了參數(shù)D和ηm之間的權(quán)衡問題。而事實上，由于印度板塊俯沖碰撞歐亞板塊的復(fù)雜性，僅以參數(shù)D來簡化印度板塊與青藏高原板塊的邊界位置并不合理。目前采用大地電磁、地震波層析成像及地質(zhì)學(xué)等方法揭示的青藏高原南緣主前緣逆沖斷裂帶(MFT)、主中央逆沖斷裂帶(MCT)和主喜馬拉雅逆沖斷裂帶(MHT)可為初步構(gòu)建區(qū)域流變結(jié)構(gòu)提供重大參考與約束。由于GPS垂直形變觀測誤差較大，影響因素較多，且模擬結(jié)果也難以完全解釋垂直觀測值，因此本文未計算其WRMS值，只能整體分析不同結(jié)果之間的差異。

如果只考慮粘彈性松弛機制的影響，則印度俯沖板片模型能更可靠地解釋震后觀測形變，尤其是垂直位移，且該模型也更接近其他方法(如接收函數(shù)法)推斷的喜馬拉雅-西藏碰撞帶的地下巖石圈結(jié)構(gòu)。總體來看，本文約束的青藏高原下地殼的穩(wěn)態(tài)粘滯系數(shù)ηm約為(1～3)×1018Pa·s，與前人的研究結(jié)果接近[8-9]；但Zhao等[3]采用該地震震后1 a的形變數(shù)據(jù)約束得到的青藏高原南緣的粘滯系數(shù)更大，為8×1018Pa·s，這是由于其考慮了更低粘滯系數(shù)的青藏高原上地幔，而本文僅固定其粘滯系數(shù)為1×1020Pa·s。

4 結(jié) 語

本文基于尼泊爾及青藏高原南緣的GPS站點觀測數(shù)據(jù)，采用有限元方法模擬分析2015年尼泊爾MW7.8地震震后5 a的粘彈性松弛效應(yīng)，探測青藏高原南緣的巖石圈流變結(jié)構(gòu)。采用Burgers體代表青藏高原下地殼的流變屬性，并用2種不同的印度板塊與青藏高原板塊的邊界結(jié)構(gòu)(垂直邊界結(jié)構(gòu)和印度彈性俯沖板片結(jié)構(gòu))來構(gòu)建模型。結(jié)果表明，2種結(jié)構(gòu)均能很好地反映遠場震后水平觀測位移，并且約束的穩(wěn)態(tài)粘滯系數(shù)較為接近，量級為1018，比瞬態(tài)Kelvin粘滯系數(shù)小1個數(shù)量級；但垂直形變差異較大，垂直邊界結(jié)構(gòu)的垂直位移嚴(yán)重偏離觀測值，而采用俯沖板片結(jié)構(gòu)模擬的垂直形變可以解釋中尼邊界及以北地區(qū)的震后隆升現(xiàn)象。后續(xù)研究可結(jié)合震后余滑機制進一步探測該區(qū)域的巖石圈流變特性。