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

1 河南理工大學測繪與國土信息工程學院，河南省焦作市世紀路2001號，454003

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

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

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

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

圖1 測站分布

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

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

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

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

(1)

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

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

圖2 震前GPS速度場

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

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

(2)

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

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

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

圖3 尼泊爾地震震后3 a水平形變場

2 孔隙彈性回彈模型

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

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

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

(3)

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

圖4 斷層滑動模型與地殼分層模型

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

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

3 震后余滑模型

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

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

(4)

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

3.1 地殼分層模型

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

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

3.2 斷層模型

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

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

圖6 有限斷層模型

3.3 震后形變模擬

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

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

圖8 GPS觀測值與震后余滑模型模擬值對比

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

圖9 斷層滑動三維及二維梯度

4 粘彈性松弛模型

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

4.1 地殼分層模型

地震造成的震后粘彈性形變與上地殼彈性層厚度和粘彈性層粘滯系數(shù)的大小有很大的關(guān)系，當彈性層厚度越薄，也就是粘彈性層越接近地表或者粘滯系數(shù)越小時，震后產(chǎn)生的地表形變就越大，反之亦然。根據(jù)徐晶等[20]和萬永革等[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]的同震斷層滑動分布模型(圖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]的同震破裂模型計算粘彈性松弛引起的震后地表形變，計算結(jié)果見圖10(圖中，藍色箭頭代表粘彈性松弛引起的地表形變理論值，紅色五角星表示同震震中和余震震中，黃色方框表示同震破裂模型)。從圖10可以看出，粘彈性松弛引起的地表形變主要分布在同震破裂區(qū)域，最大位移量約為7.9 cm(27.95°N，85.08°E)，整體表現(xiàn)為向北的運動趨勢，與同震破裂面引起的地表形變變化趨勢相反，但在同震破裂面的南部遠場區(qū)域向北運動，在同震破裂面的北部遠場區(qū)域向南運動，與同震引起的形變運動方向基本一致。

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

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

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

5 組合機制模型

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

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

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

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

圖13 扣除粘彈性的GPS觀測值與模擬值對比

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

圖14 去除粘彈性后斷層滑動三維及二維梯度

6 結(jié) 語

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

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

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

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

4)組合機制模型反演的震后余滑主要發(fā)生在同震破裂面下傾部分，而同震未破裂斷層的淺部未檢測到滑動信息，說明斷層淺部未發(fā)生破裂，未來的地震危險性較高。