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

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

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

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

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

1 震后GPS觀測數據處理

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

圖1 測站分布

1.1 數據處理

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

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

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

(1)

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

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

圖2 震前GPS速度場

1.2 GPS觀測的震后形變與函數擬合

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

(2)

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

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

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

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

2 孔隙彈性回彈模型

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

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

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

(3)

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

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

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

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

3 震后余滑模型

設置反演的目標函數為[10]：

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

(4)

式中，W為觀測值的權矩陣，G為格林函數，d為震后位移觀測值，s為滑動量，β為平滑因子，用于調節模型粗糙度和數據擬合誤差之間的關系，?為有限差分拉普拉斯算子。

3.1 地殼分層模型

表1 尼泊爾地區地殼結構模型

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

3.2 斷層模型

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

表2 預設斷層模型參數

圖6 有限斷層模型

3.3 震后形變模擬

基于§3.1地殼分層模型計算的格林函數和斷層模型參數，采用SDM反演程序對震后的斷層余滑分布進行反演，反演過程中需要不斷調整平滑因子，以平衡模型粗糙度和數據擬合誤差之間的關系。圖7為余滑模型的模型粗糙度和數據擬合誤差的折中曲線(圖中圓點旁邊的數字為對應的平滑因子)，可以看出，當β=0.1時可較好地平衡模型粗糙度和數據擬合誤差之間的關系，此時模型值與觀測值的擬合誤差為2.71 mm。從圖8可以看出，震后余滑模型可較好地解釋GPS觀測的震后地表水平位移。

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

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

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

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

4 粘彈性松弛模型

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

4.1 地殼分層模型

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

4.2 斷層模型

使用Hayes等[16]的同震斷層滑動分布模型(圖4)，斷層幾何參數為：走向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]的同震破裂模型計算粘彈性松弛引起的震后地表形變，計算結果見圖10(圖中，藍色箭頭代表粘彈性松弛引起的地表形變理論值，紅色五角星表示同震震中和余震震中，黃色方框表示同震破裂模型)。從圖10可以看出，粘彈性松弛引起的地表形變主要分布在同震破裂區域，最大位移量約為7.9 cm(27.95°N，85.08°E)，整體表現為向北的運動趨勢，與同震破裂面引起的地表形變變化趨勢相反，但在同震破裂面的南部遠場區域向北運動，在同震破裂面的北部遠場區域向南運動，與同震引起的形變運動方向基本一致。

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

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

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

5 組合機制模型

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

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

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

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

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

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

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

6 結 語

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

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

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

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

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