基于GNSS觀測(cè)研究2015年尼泊爾MW7.8地震震后地殼形變特征及其機(jī)制

姚未正 徐克科 朱緒林 邵振華

1 河南理工大學(xué)測(cè)繪與國(guó)土信息工程學(xué)院，河南省焦作市世紀(jì)路2001號(hào)，454003

相對(duì)于同震形變而言，大地震震后形變的持續(xù)時(shí)間更長(zhǎng)、影響范圍更廣。地震震后形變的研究有助于認(rèn)識(shí)地震的震后形變機(jī)制，探索地球深部巖石圈流變結(jié)構(gòu)，為判定后續(xù)地震危險(xiǎn)性提供重要依據(jù)[1]。前人研究結(jié)果表明，震后形變機(jī)制主要包括發(fā)生在地殼淺部的孔隙彈性形變、發(fā)生在同震破裂面上的震后余滑形變和發(fā)生在下地殼和上地幔的粘彈性松弛形變[2-4]。

據(jù)中國(guó)地震臺(tái)網(wǎng)測(cè)定，北京時(shí)間2015-04-25尼泊爾首都加德滿都西北約80 km發(fā)生MW7.8地震，震中位于84.70°E、28.20°N，震源深度約為20 km[5-7]，美國(guó)地質(zhì)調(diào)查局(USGS)通過(guò)對(duì)全球數(shù)字地震臺(tái)網(wǎng)地震波形數(shù)據(jù)進(jìn)行反演，確定了斷層破裂模型的滑動(dòng)幅度和滑動(dòng)方向，發(fā)現(xiàn)此次地震是一次以逆沖為主，兼少量右旋走滑運(yùn)動(dòng)的低傾角事件，最大同震位移達(dá)到3～4 m。地震發(fā)生后，各國(guó)研究機(jī)構(gòu)迅速在尼泊爾及其周邊地區(qū)新增了許多GPS站點(diǎn)，并對(duì)舊GPS站點(diǎn)進(jìn)行復(fù)測(cè)，獲得了大量的震后形變觀測(cè)資料，為研究震后形變機(jī)制提供了數(shù)據(jù)基礎(chǔ)。目前已有的關(guān)于尼泊爾地震震后形變機(jī)制的研究大多采用的是早期震后形變資料[8-14]，研究對(duì)象也主要是早期的震后余滑，雖然有部分學(xué)者在研究尼泊爾地震震后形變機(jī)制時(shí)綜合考慮了粘彈性松弛效應(yīng)的影響，但使用的觀測(cè)數(shù)據(jù)時(shí)間尺度相對(duì)較短；部分學(xué)者采用了較長(zhǎng)時(shí)間的震后觀測(cè)數(shù)據(jù)，但使用的觀測(cè)站點(diǎn)較少，無(wú)法很好地約束反演模型。因此，本文收集了尼泊爾境內(nèi)和中國(guó)藏南區(qū)域共32個(gè)GPS連續(xù)站震后3 a的觀測(cè)資料，采用對(duì)數(shù)函數(shù)模型對(duì)位移時(shí)間序列進(jìn)行擬合，分離區(qū)域構(gòu)造運(yùn)動(dòng)和同震、余震活動(dòng)的影響，獲得尼泊爾地震震后3 a的累積形變場(chǎng)，再分別使用孔隙彈性回彈模型、震后余滑模型、粘彈性松弛模型和組合機(jī)制模型研究尼泊爾地震的震后形變機(jī)理。

1 震后GPS觀測(cè)數(shù)據(jù)處理

本文收集了美國(guó)內(nèi)華達(dá)州大地測(cè)量實(shí)驗(yàn)室(http:∥geodesy.unr.edu/)處理的尼泊爾境內(nèi)28個(gè)測(cè)站的原始觀測(cè)數(shù)據(jù)及中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(http:∥www.cgps.ac.cn)處理的中國(guó)藏南區(qū)域4個(gè)測(cè)站的原始觀測(cè)數(shù)據(jù)。圖1為32個(gè)測(cè)站的位置分布，圖中的震源機(jī)制球分別為USGS測(cè)定的尼泊爾地震主震和最大余震的震中位置，黑色三角形為尼泊爾境內(nèi)測(cè)站，黑色正方形為中國(guó)境內(nèi)測(cè)站。

圖1 測(cè)站分布

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

GPS時(shí)間序列中主要包含構(gòu)造運(yùn)動(dòng)信號(hào)、非構(gòu)造運(yùn)動(dòng)信號(hào)和季節(jié)性周期信號(hào)，其中非構(gòu)造運(yùn)動(dòng)信號(hào)包含地震的同震和震后形變及其他因素引起的形變，一般可用式(1)表示：

y(ti)=y0+vti+Asin(2πti)+Bcos(2πti)+

Fln(1+(ti-teq)/τi)

(1)

式中，y(ti)為GPS站點(diǎn)在歷元ti時(shí)刻的觀測(cè)值，ti為歷元數(shù)，v為長(zhǎng)期線性速度值，A、B、C、D分別為周期信號(hào)振幅，E為由其他原因引起的階躍，H(ti-tcj)為階躍函數(shù)，F(xiàn)為震后形變振幅，τi為松弛時(shí)間常數(shù)。其中，F(xiàn)和τi存在較強(qiáng)的負(fù)相關(guān)性，不能同時(shí)估計(jì)。

為從原始坐標(biāo)時(shí)間序列中分離出地震的震后形變信號(hào)，需要對(duì)時(shí)間序列中的構(gòu)造運(yùn)動(dòng)信號(hào)、季節(jié)性周期信號(hào)和其他階躍信號(hào)進(jìn)行修正。由于周期性水平形變的幅度較小，本文在計(jì)算尼泊爾地震震后水平形變場(chǎng)時(shí)不考慮季節(jié)性周期信號(hào)的影響。而在構(gòu)造運(yùn)動(dòng)信號(hào)獲取方面，由于部分GPS測(cè)站是震后新增的，沒(méi)有震前觀測(cè)數(shù)據(jù)，可先利用最小二乘參數(shù)擬合的方法獲取其他測(cè)站的震前速度場(chǎng)，再采用克里金插值的方法獲得這些新增測(cè)站的震前速度場(chǎng)，計(jì)算結(jié)果見(jiàn)圖2(ITRF2014框架下，圖中藍(lán)色箭頭表示利用震前時(shí)間序列計(jì)算的長(zhǎng)期運(yùn)動(dòng)速率，紅色箭頭表示利用克里金插值估算的震前運(yùn)動(dòng)速率)；對(duì)于同震及余震產(chǎn)生的階躍信號(hào)的影響，可采用在擬合函數(shù)中添加初始常數(shù)的方法解決。

圖2 震前GPS速度場(chǎng)

1.2 GPS觀測(cè)的震后形變與函數(shù)擬合

利用計(jì)算的震前背景速度場(chǎng)可以對(duì)GPS原始坐標(biāo)時(shí)間序列進(jìn)行構(gòu)造形變改正，得到改正后的震后坐標(biāo)時(shí)間序列。由于各個(gè)GPS測(cè)站的觀測(cè)時(shí)段不一，為計(jì)算每個(gè)測(cè)站震后3 a的位移量，需要對(duì)改正后的震后坐標(biāo)時(shí)間序列進(jìn)行函數(shù)擬合。本文采用式(2)對(duì)數(shù)模型進(jìn)行擬合：

(2)

式中，F(xiàn)為震后形變振幅，Δt為震后逝去時(shí)間，τ為震后松弛時(shí)間，k為去除同震及余震階躍影響添加的時(shí)間序列初始常數(shù)。

外采低價(jià)資源雖然利潤(rùn)豐厚，但其前提條件是通過(guò)加油站零售出去。受制于網(wǎng)絡(luò)能力限制，每年外采2000 萬(wàn)噸資源只有一部分可以進(jìn)入零售環(huán)節(jié)，其余資源基本相當(dāng)于賺取批發(fā)差價(jià)。

由于F和τ具有較強(qiáng)的負(fù)相關(guān)性，很難同時(shí)反演確定，本文通過(guò)大量試算，最終固定τ為38 d，將非線性方程轉(zhuǎn)化為線性方程，擬合每個(gè)GPS站的觀測(cè)分量，擬合的平均水平誤差為2.6 mm，計(jì)算的震后水平位移如圖3所示?？梢钥闯?，震后形變主要集中在尼泊爾北部區(qū)域，且東西方向形變較小，南北方向形變較大，整體繼續(xù)向南運(yùn)動(dòng)，最大震后地表水平位移發(fā)生在CHLM站，約10.93 cm，為西南方向。

圖3 尼泊爾地震震后3 a水平形變場(chǎng)

2 孔隙彈性回彈模型

一般來(lái)說(shuō)，地球介質(zhì)中都會(huì)有小到分子量級(jí)的孔隙，而這些孔隙的存在可以使某些分子順利通過(guò)，這種存在大量密集微小孔隙的固體物質(zhì)被稱(chēng)為孔隙介質(zhì)或多孔介質(zhì)。地震發(fā)生前，地殼介質(zhì)處于未排空水狀態(tài)；地震發(fā)生時(shí)，累積的地震應(yīng)變能量迅速爆發(fā)，斷層發(fā)生破裂，導(dǎo)致周?chē)目紫端畨毫ι?，巖體中孔隙水的平衡狀態(tài)被破壞；震后短時(shí)間內(nèi)孔隙水由高壓力地區(qū)不斷流向低壓力地區(qū)，孔隙水壓再次恢復(fù)到另一種平衡狀態(tài)[15]。在這個(gè)過(guò)程中，孔隙水的流動(dòng)使巖石間的應(yīng)力發(fā)生變化，進(jìn)而導(dǎo)致地表形變或發(fā)生余震。

考慮到孔隙流體的反彈時(shí)間尺度非常依賴(lài)于巖石的可滲透率，通常會(huì)持續(xù)數(shù)月到數(shù)年，且其時(shí)間衰減是非線性的，使得孔隙回彈形變的計(jì)算過(guò)程非常復(fù)雜。因此，通常利用同震破裂模型分別計(jì)算地震在未排空水狀態(tài)和完全排空水狀態(tài)下產(chǎn)生的地表形變，兩者的差值即可近似為孔隙彈性回彈造成的地表形變，具體計(jì)算公式為：

ui(x)=ui(tq,v)-ui(tq,vu)

(3)

式中，ui(x)為x點(diǎn)在i方向上的位移，ui(tq,v)和ui(tq,vu)分別為未排空水狀態(tài)和完全排空水狀態(tài)下計(jì)算得到的地表形變。假定未排空水狀態(tài)下的泊松比為0.27，完全排空水狀態(tài)下的泊松比為0.25，采用Hayes等[16]公布的同震破裂模型和地殼分層模型(圖4)分別計(jì)算2種狀態(tài)下的地表形變，并對(duì)其進(jìn)行差分，得到孔隙彈性回彈引起的地表位移場(chǎng)，計(jì)算結(jié)果見(jiàn)圖5(圖中紅色箭頭表示GPS觀測(cè)的震后3 a位移，藍(lán)色箭頭表示孔隙彈性回彈引起的地表理論位移)。

圖4 斷層滑動(dòng)模型與地殼分層模型

圖5 孔隙彈性回彈引起的地表形變

從圖5可以看出，孔隙彈性回彈引起的地表理論位移主要集中在同震破裂周?chē)?，最大形變量約為2.8 mm(27.72°N，85.46°E)，向四周呈發(fā)散狀，且離震中區(qū)域越遠(yuǎn)，形變量越小。此外，通過(guò)對(duì)比孔隙彈性回彈引起的地表理論位移與GPS觀測(cè)的震后形變值可以發(fā)現(xiàn)，在多數(shù)站點(diǎn)二者沒(méi)有一致性，盡管在有些站點(diǎn)二者的運(yùn)動(dòng)方向一致，但孔隙彈性回彈產(chǎn)生的形變量遠(yuǎn)小于GPS觀測(cè)值。因此可以認(rèn)為，使用該模型不能解釋尼泊爾地震的震后形變，且由于孔隙彈性回彈引起的地表位移貢獻(xiàn)量較小，在后續(xù)研究中不予考慮。

3 震后余滑模型

設(shè)置反演的目標(biāo)函數(shù)為[10]：

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

(4)

式中，W為觀測(cè)值的權(quán)矩陣，G為格林函數(shù)，d為震后位移觀測(cè)值，s為滑動(dòng)量，β為平滑因子，用于調(diào)節(jié)模型粗糙度和數(shù)據(jù)擬合誤差之間的關(guān)系，?為有限差分拉普拉斯算子。

3.1 地殼分層模型

表1 尼泊爾地區(qū)地殼結(jié)構(gòu)模型

根據(jù)表1的地殼分層模型，采用SDM程序附帶的格林函數(shù)計(jì)算程序，基于分層半無(wú)限空間模型計(jì)算尼泊爾地區(qū)的位移格林函數(shù)。

3.2 斷層模型

由于震后余滑通常發(fā)生在同震破裂面或其延伸面區(qū)域，因此在構(gòu)建震后余滑反演模型時(shí)可以以同震破裂模型為基礎(chǔ)，并在長(zhǎng)度和寬度上進(jìn)行適當(dāng)延伸，以保證其能覆蓋整個(gè)震后余滑的空間范圍。本文以Hayes等[16]的同震破裂模型為基礎(chǔ)，取長(zhǎng)度為260 km、寬度為200 km，并將其沿?cái)鄬幼呦蚝蛢A向離散成26×20個(gè)10 km×10 km的子斷層，具體模型參數(shù)見(jiàn)表2。斷層模型的實(shí)際地理位置如圖6所示。

表2 預(yù)設(shè)斷層模型參數(shù)

圖6 有限斷層模型

3.3 震后形變模擬

基于§3.1地殼分層模型計(jì)算的格林函數(shù)和斷層模型參數(shù)，采用SDM反演程序?qū)φ鸷蟮臄鄬佑嗷植歼M(jìn)行反演，反演過(guò)程中需要不斷調(diào)整平滑因子，以平衡模型粗糙度和數(shù)據(jù)擬合誤差之間的關(guān)系。圖7為余滑模型的模型粗糙度和數(shù)據(jù)擬合誤差的折中曲線(圖中圓點(diǎn)旁邊的數(shù)字為對(duì)應(yīng)的平滑因子)，可以看出，當(dāng)β=0.1時(shí)可較好地平衡模型粗糙度和數(shù)據(jù)擬合誤差之間的關(guān)系，此時(shí)模型值與觀測(cè)值的擬合誤差為2.71 mm。從圖8可以看出，震后余滑模型可較好地解釋GPS觀測(cè)的震后地表水平位移。

圖7 余滑模型的模型粗糙度與數(shù)據(jù)擬合誤差折中曲線

圖8 GPS觀測(cè)值與震后余滑模型模擬值對(duì)比

對(duì)應(yīng)的震后余滑模型見(jiàn)圖9，可以看出，反演的平均滑動(dòng)量約為8.6 cm，最大滑動(dòng)量為34.39 cm，位于地下約26.6 km處,震后余滑釋放的地震矩為1.09×1020Nm，對(duì)應(yīng)的矩震級(jí)為MW7.33。對(duì)比同震破裂模型可以看出，震后余滑在淺部同震未破裂區(qū)域的分布極其有限，主要分布在20～30 km的同震破裂面下傾延伸部分，且分布范圍廣泛。Mencin等[18]采用應(yīng)力驅(qū)動(dòng)余滑模型反演的震后余滑發(fā)生在同震破裂面的正下方，而出現(xiàn)這種模型斷層底端滑動(dòng)的現(xiàn)象可能是由于忽略了粘彈性松弛效應(yīng)導(dǎo)致的。

圖9 斷層滑動(dòng)三維及二維梯度

4 粘彈性松弛模型

下地殼和上地幔的物質(zhì)并不是完全的彈性體，通常視為粘性流體。地震發(fā)生后，受到同震和震后應(yīng)力變化的影響，具有粘彈性質(zhì)的中下地殼和上地幔應(yīng)力狀態(tài)隨時(shí)間緩慢變化，使得上地殼的彈性層發(fā)生應(yīng)力擾動(dòng)，進(jìn)而引發(fā)大范圍的地表形變，這種震后形變稱(chēng)為震后粘彈性松弛。由于粘彈性松弛效應(yīng)引起的震后形變的空間尺度和時(shí)間尺度都較大，因此在進(jìn)行震后余滑反演時(shí)需要考慮其影響，否則會(huì)高估斷層深部的震后余滑[19]。

4.1 地殼分層模型

地震造成的震后粘彈性形變與上地殼彈性層厚度和粘彈性層粘滯系數(shù)的大小有很大的關(guān)系，當(dāng)彈性層厚度越薄，也就是粘彈性層越接近地表或者粘滯系數(shù)越小時(shí)，震后產(chǎn)生的地表形變就越大，反之亦然。根據(jù)徐晶等[20]和萬(wàn)永革等[21]的研究，青藏高原中下地殼的粘滯系數(shù)在1019～1020Pa·s之間，上地幔粘滯系數(shù)約為1.0×1020Pa·s，由此可確定尼泊爾地區(qū)的地殼分層模型：0～20 km為上地殼彈性層；20～51 km為中下地殼；51 km以下為上地幔及以下區(qū)域。使用Maxwell流體模擬中下地殼和上地幔的流變結(jié)構(gòu)，其中，中下地殼的粘滯系數(shù)取1.0×1019Pa·s，上地幔的粘滯系數(shù)取1.0×1020Pa·s。

4.2 斷層模型

使用Hayes等[16]的同震斷層滑動(dòng)分布模型(圖4)，斷層幾何參數(shù)為：走向293°，傾角7°，頂部埋深4.3 km，沿走向和傾向劃分為等間隔的23×24塊子斷層，分割成約8.4 km×7 km的子斷層塊，覆蓋范圍約193.2 km×168 km。

4.3 震后形變模擬

使用Wang等[22]的PSGRN/PSCMP程序，基于§4.1的地殼分層模型和Hyaes等[6]的同震破裂模型計(jì)算粘彈性松弛引起的震后地表形變，計(jì)算結(jié)果見(jiàn)圖10(圖中，藍(lán)色箭頭代表粘彈性松弛引起的地表形變理論值，紅色五角星表示同震震中和余震震中，黃色方框表示同震破裂模型)。從圖10可以看出，粘彈性松弛引起的地表形變主要分布在同震破裂區(qū)域，最大位移量約為7.9 cm(27.95°N，85.08°E)，整體表現(xiàn)為向北的運(yùn)動(dòng)趨勢(shì)，與同震破裂面引起的地表形變變化趨勢(shì)相反，但在同震破裂面的南部遠(yuǎn)場(chǎng)區(qū)域向北運(yùn)動(dòng)，在同震破裂面的北部遠(yuǎn)場(chǎng)區(qū)域向南運(yùn)動(dòng)，與同震引起的形變運(yùn)動(dòng)方向基本一致。

圖10 粘彈性松弛引起的地表形變

圖11為震后3 a的GPS觀測(cè)值與粘彈性松弛模擬值的對(duì)比?？梢钥闯觯鸷笮巫冎?，粘彈性松弛模型并不能解釋近場(chǎng)的地表形變，其遠(yuǎn)場(chǎng)形變的運(yùn)動(dòng)方向與GPS觀測(cè)值一致，但形變量小于GPS觀測(cè)值，這可能是由尼泊爾地區(qū)和青藏高原地區(qū)地球模型的橫向不均勻性導(dǎo)致的。

圖11 震后3 a累積形變值與粘彈性松弛模擬值比較

5 組合機(jī)制模型

本文分別使用孔隙彈性回彈模型、粘彈性松弛模型和震后余滑模型來(lái)解釋GPS觀測(cè)到的震后形變。結(jié)果發(fā)現(xiàn)，單獨(dú)的震后形變機(jī)理并不能合理地解釋尼泊爾地震的GPS觀測(cè)資料，其中孔隙彈性回彈引起的地表形變作用范圍和量級(jí)均較小，無(wú)法解釋GPS觀測(cè)的震后形變；粘彈性松弛模型則可以解釋遠(yuǎn)場(chǎng)的GPS觀測(cè)值，但同樣無(wú)法解釋近場(chǎng)的地表形變；震后余滑模型的擬合程度最高，但反演的震后余滑主要分布在同震破裂的下傾延伸部分，最大余滑深度接近地下30 km，且分布范圍廣泛。

因此，為同時(shí)兼顧近場(chǎng)、遠(yuǎn)場(chǎng)、擬合誤差和物理意義，本文采用震后余滑+粘彈性松弛組合機(jī)制模型進(jìn)行研究。首先從圖3的GPS觀測(cè)值中扣除圖11中粘彈性松弛引起的理論地表形變，然后以修正過(guò)的GPS觀測(cè)值為約束條件，結(jié)合§3的地殼分層模型和斷層滑動(dòng)模型，利用SDM反演程序重新反演震后余滑分布模型。

圖12為去除粘彈性的余滑模型的模型粗糙度與數(shù)據(jù)擬合誤差的折中曲線。可以看出，β=0.1為最佳光滑因子，此時(shí)的觀測(cè)值與模擬值的擬合誤差為2.72 mm。從圖13可以看出，組合機(jī)制模型同樣可以很好地解釋震后GPS觀測(cè)值。

圖12 去除粘彈性后余滑模型的模型粗糙度與數(shù)據(jù)擬合誤差的折中曲線

圖13 扣除粘彈性的GPS觀測(cè)值與模擬值對(duì)比

圖14為經(jīng)粘彈性松弛修正后的震后余滑模型?？梢钥闯?，修正后的震后余滑主要集中在17～26 km的同震破裂下傾部分，反演的最大滑動(dòng)量為38.67 cm，深度約為21.17 km，整個(gè)斷層面的平均滑動(dòng)量約為8.5 cm，震后余滑釋放的地震矩也下降到1.08×1020Nm，對(duì)應(yīng)的矩震級(jí)為MW7.32。修正后的震后余滑模型在保證數(shù)據(jù)擬合度的基礎(chǔ)上，平均滑動(dòng)量和地震矩較修正前都有所降低，且在空間分布上也有所升高和集中，這與Mencin等[18]采用應(yīng)力驅(qū)動(dòng)模型反演的結(jié)果更接近。

圖14 去除粘彈性后斷層滑動(dòng)三維及二維梯度

6 結(jié) 語(yǔ)

本文收集2015年尼泊爾MW7.8地震震區(qū)附近的GPS連續(xù)站觀測(cè)資料，獲取尼泊爾地震震后3 a的水平形變場(chǎng)，并以此為約束，分別采用孔隙彈性回彈模型、震后余滑模型、粘彈性松弛模型和組合機(jī)制模型研究尼泊爾地震的震后形變機(jī)制。結(jié)果表明：

1)GPS觀測(cè)到的尼泊爾地震震后形變主要集中在尼泊爾北部區(qū)域，且東西方向形變較小，南北方向形變較大，整體向南運(yùn)動(dòng)。地表最大震后位移發(fā)生在CHLM站，約為10.93 cm，向西南方向。

2)分別采用孔隙彈性回彈模型、震后余滑模型和粘彈性松弛模型研究尼泊爾地震的震后形變機(jī)理，模擬結(jié)果顯示，孔隙彈性回彈模型模擬計(jì)算的理論地表形變空間分布范圍較小，且遠(yuǎn)小于GPS觀測(cè)的震后形變；SDM程序反演的震后余滑模型雖然可以較好地?cái)M合GPS觀測(cè)值，但反演的震后余滑分布范圍廣泛，且主要集中在同震破裂面的下傾延伸部分；PSGRN/PSCMP程序計(jì)算的粘彈性松弛引起的理論地表形變無(wú)法解釋近場(chǎng)形變，只有在遠(yuǎn)場(chǎng)區(qū)域其運(yùn)動(dòng)方向才與GPS觀測(cè)值一致，這可能是尼泊爾地區(qū)和青藏高原地區(qū)地球模型的橫向不均勻性導(dǎo)致的。

3)利用震后余滑和粘彈性松弛的組合機(jī)制模型研究尼泊爾地震的震后形變機(jī)制后發(fā)現(xiàn)，組合機(jī)制模型在保證數(shù)據(jù)擬合度的基礎(chǔ)上，其平均滑動(dòng)量和地震矩較修正前都有所降低，且在空間分布上也有所升高和集中，與應(yīng)力驅(qū)動(dòng)模型的反演結(jié)果更為接近。

4)組合機(jī)制模型反演的震后余滑主要發(fā)生在同震破裂面下傾部分，而同震未破裂斷層的淺部未檢測(cè)到滑動(dòng)信息，說(shuō)明斷層淺部未發(fā)生破裂，未來(lái)的地震危險(xiǎn)性較高。