劉志鵬， 蓋增喜

北京大學(xué)地球與空間科學(xué)學(xué)院地球物理系， 北京 100871

?

利用多臺陣壓縮傳感方法反演尼泊爾MW7.9地震破裂過程

劉志鵬， 蓋增喜*

北京大學(xué)地球與空間科學(xué)學(xué)院地球物理系， 北京 100871

2015年4月25日，尼泊爾地區(qū)發(fā)生MW7.9地震，震中位于28.1°N，84.7°E.為了詳細(xì)地研究此次破壞性極強(qiáng)的地震的破裂過程，本文利用多臺陣壓縮傳感方法，使用了阿拉斯加、歐洲和澳大利亞三個臺網(wǎng)的共計(jì)179個臺站的遠(yuǎn)場P波垂直分量的數(shù)據(jù)來反演，結(jié)果表明本次地震的破裂過程是一個清晰的南東東方向的單側(cè)破裂，破裂尺度約為105 km，整體持續(xù)時間約為58 s.在破裂初始的前15 s，能量輻射基本圍繞在震源附近，16 s后破裂開始向南東東方向以1.9 km·s-1的速度破裂.釋放能量最大的時間為第38 s，位于距震中70 km處.該位置從第29秒開始破裂，并持續(xù)釋放能量長達(dá)30 s之久.

多臺陣壓縮傳感； 尼泊爾地震； 地震破裂過程

1 引言

北京時間2015年4月25日14時11分，在尼泊爾博克拉發(fā)生MW7.9地震，震中位于28.1°N，84.7°E，震源深度15 km(由USGS(United States Geological Survey，美國地質(zhì)調(diào)查局)給出).此次地震位于人口相對密集、建筑物抗震性能較差的地區(qū)，因此造成了尼泊爾境內(nèi)的大量人員傷亡.中國西藏、印度、孟加拉國、不丹等地也均出現(xiàn)人員傷亡.此次地震發(fā)生在印度板塊和歐亞板塊的陸陸碰撞帶上并且震級較大，所以也引發(fā)了雪崩、滑坡等次生災(zāi)害，造成重大人員傷亡和財(cái)產(chǎn)損失.本次地震發(fā)生在印度板塊向北俯沖到歐亞板塊的匯聚邊界上，在其西部遠(yuǎn)處1505年發(fā)生過一次更大的地震，而緊鄰本次地震的東部，在1934年發(fā)生過一次M8.0地震，該地震附近在1255年還發(fā)生過一次更大的地震.考慮到本次地震對加德滿都乃至尼泊爾全境造成的嚴(yán)重?fù)p失，系統(tǒng)研究本次地震的破裂過程并準(zhǔn)確給出能量釋放極大的位置，對于揭示當(dāng)?shù)匕鍓K構(gòu)造活動的規(guī)律和抗震救災(zāi)都有著非常重要的意義.

地震發(fā)生后，中國地震局地球物理研究所給出了有限斷層反演的解(http:∥www.cea-igp.ac.cn/ tpxw/272110.shtml)，其結(jié)果表明主要能量釋放區(qū)域在震中與加德滿都之間的位置.不同于有限斷層反演，在當(dāng)?shù)財(cái)鄬拥膸缀涡螒B(tài)未知的情況下，反投影方法是近十年來發(fā)展起來的獲取大地震破裂過程的穩(wěn)定有效的方法(Zhang and Ge，2010; Ishii et al.，2005; Koper et al.，2012).IRIS(Incorporated Research Institutions for Seismology，美國地震學(xué)研究聯(lián)合會)網(wǎng)站給出了不同臺陣的反投影結(jié)果(http:∥ds.iris.edu/spud/backprojection/9924467)，但由于單一臺陣的方位角局限，這些結(jié)果表現(xiàn)出了較大的差異.為了提高破裂的分辨率，多重信號分類算法(Schmidt，1986; Goldstein and Archuleta，1991)也被引入了地球物理學(xué)領(lǐng)域來得到震源破裂的信息(Meng et al.，2011，2012).作為反投影方法和反演方法的結(jié)合，壓縮傳感算法(Donoho，2006; Candès et al.，2006; Malioutov et al.，2005)基于反演計(jì)算而不是直接的波形聚束，可以復(fù)現(xiàn)出稀疏分布的源而不依賴于任何的先驗(yàn)信息.這種方法已經(jīng)被廣泛使用在了核磁共振成像(Lustig et al.，2007)和計(jì)算攝影學(xué)(Duarte et al.，2008)等領(lǐng)域.在一個給定的時間，對于網(wǎng)格化的空間中，源的分布應(yīng)該較為離散，數(shù)量較少，這一點(diǎn)符合稀疏性的前提條件，并滿足震源的物理假設(shè).因此，基于某一個臺陣的大量數(shù)據(jù)，可以將壓縮傳感方法應(yīng)用在地震學(xué)領(lǐng)域來對大地震的時空演化過程進(jìn)行成像(Yao et al.，2011).但是，在傳統(tǒng)的壓縮傳感方法中，為了數(shù)據(jù)之間更好的相關(guān)性，一般只會使用由一個臺陣記錄到的P波數(shù)據(jù)來進(jìn)行反演.

考慮到一個臺陣的方位角的局限性可能會帶來一些誤差，本文將傳統(tǒng)的壓縮方法進(jìn)行改進(jìn)，使用多個不同方位角的臺陣的數(shù)據(jù)進(jìn)行反演，并將其稱為多臺陣壓縮傳感方法.本文利用該方法，使用阿拉斯加、歐洲和澳大利亞三個臺網(wǎng)共179個臺站的遠(yuǎn)場P波垂直分量的數(shù)據(jù)來反演本次MW7.9尼泊爾地震的破裂過程.

2 方法和數(shù)據(jù)

地震位錯點(diǎn)源產(chǎn)生的位移場可以表示為(Aki and Richards，2002)

(1)

對于有限斷層，考慮震源的有限破裂尺度效應(yīng)，則(1)式可修改為

u(x,t)=?Smpq(ξ,τ)*Gp,q(x,t;ξ,τ)dS(ξ),

(2)

寫成離散的形式，則有

(3)

其中下標(biāo)k為子破裂元的編號，ξk，τk，Sk分別為子破裂元k的位置坐標(biāo)、破裂發(fā)生時間和破裂面積.由于考慮的是遠(yuǎn)場的P波垂直分量，并使用歸一化的波形而非振幅的絕對值來計(jì)算每個子破裂元在頻率域釋放能量的相對強(qiáng)度，則意味著格林函數(shù)Gip,q中有關(guān)輻射花樣和幾何衰減因子的項(xiàng)可以忽略，從而格林函數(shù)Gip,q在頻率域可簡化為相移項(xiàng)eiω t,其中t為從子破裂元到地震臺站的走時.子破裂元的矩張量mp qk=μSkuk.則(3)式可重新改寫為矩陣形式

(4)

因?yàn)槭褂玫氖嵌鄠€臺陣的數(shù)據(jù)，因此對于基于單臺陣的傳統(tǒng)壓縮傳感方法的解的形式(Yaoetal.，2011)就轉(zhuǎn)化為

+λ‖x‖1].

(5)

這里，j代表的是臺陣的數(shù)量，i是某個臺陣的序號.λ是一個正值的拉格朗日乘子.wi是每一個臺陣的權(quán)重因子，它的大小取決于這個臺陣的孔徑和數(shù)據(jù)質(zhì)量.因此，這個反演問題的解是數(shù)據(jù)殘差和解的稀疏性之間的折中.由于遠(yuǎn)震P波對深度幾乎沒有分辨，所以我們只在震源的深度平面內(nèi)進(jìn)行反演(Xuetal.，2009).

為了使用多臺陣壓縮傳感方法反演此次地震的破裂過程，我們使用了來自三個臺陣(歐洲 79臺站，阿拉斯加 81臺站，澳大利亞19臺站)，共計(jì)179個臺站的遠(yuǎn)場P波垂直分量的數(shù)據(jù)(圖1，圖2).首先，所有波形都進(jìn)行0.2～5Hz的帶通濾波，然后分別根據(jù)每一個臺站的數(shù)據(jù)最大值進(jìn)行歸一化.之所以選擇這個頻段進(jìn)行濾波是因?yàn)檫^于低頻的信息空間的分辨率太低，很難獲得精細(xì)的細(xì)節(jié)信息；高頻能量由于噪聲過大可能導(dǎo)致結(jié)果的不穩(wěn)定.然后，為了減少由于臺站所在的接收區(qū)域的淺層速度結(jié)構(gòu)的三維非均勻性對于走時的影響，我們先利用PREM模型(DziewonskiandAnderson，1981)計(jì)算出了由震源(28.147°N，84.708°E，15km)到各個臺站的理論到時，并通過計(jì)算最早的7s數(shù)據(jù)的互相關(guān)函數(shù)，將同一臺陣內(nèi)的數(shù)據(jù)對齊到零時刻.隨后，由于到達(dá)不同臺陣的射線路徑不同，也會影響到時，因此需要將每個臺陣近似看作一個觀測點(diǎn)，對齊所有來自不同臺陣的數(shù)據(jù).為了達(dá)到此目的，我們分別將每個臺陣已經(jīng)對齊好的數(shù)據(jù)分別疊加起來獲得三條疊加的數(shù)據(jù)，從而提高信噪比，達(dá)到P波初動可以清晰識別的目的.最后，將疊加好的數(shù)據(jù)通過人工標(biāo)定的P波初動時間作為參考時間來對齊并進(jìn)行相應(yīng)的極性翻轉(zhuǎn)(圖2).

為了獲得此次地震震源區(qū)域不同位置能量釋放的時間分布，我們使用了一個長度為10 s的滑動窗，滑動的步長為2 s.第一個時間窗被設(shè)定為只包含P波初動的波形，這是為了更好地在結(jié)果中識別出破裂發(fā)生的初期的性質(zhì).在反演計(jì)算的過程中，每一個時間窗的數(shù)據(jù)都轉(zhuǎn)換到頻率域，并計(jì)算在0.25～1 Hz頻段范圍內(nèi)的解.之所以選擇這個頻段是由于數(shù)據(jù)的頻譜決定的(圖2d).如果頻段選擇的過低，則無法看到破裂的細(xì)節(jié)；如果選擇的過高，則會由于高頻的能量較少信噪比較低而使得解不夠精確.在反演過程結(jié)束后，對于每一個選定的時間窗，我們得到了不同頻率下能量釋放的空間分布，并將這些分布疊加起來，得到在0.25～1 Hz頻帶范圍內(nèi)的能量分布.最終，將每一個時間窗的解進(jìn)行二維的空間高斯平滑，得到了震源區(qū)域能量的釋放位置隨時間的變化過程，即此次地震的時空演化過程.

圖1 尼泊爾地震的震中位置以及所使用的臺陣數(shù)據(jù)的分布黑色五角星代表震中的位置，藍(lán)色圓形代表臺站.Fig.1 Map of the Nepal earthquake location and distribution of stations from three arrays Black star is the epicenter. Blue circles are the stations.

圖2 本文計(jì)算中所使用的數(shù)據(jù)的信息(a)歐洲,(b)阿拉斯加,(c)澳大利亞三個臺陣上方的彩色圖案代表地震波形數(shù)據(jù)及各臺陣以P波初值對齊后疊加的地震圖，箭頭標(biāo)示P波初動；(d)圖為臺網(wǎng)中代表性臺站的數(shù)據(jù)的頻譜的范例.Fig.2 The information about used data(a) Europe, (b) Alaska, (c) Australia. In the figure of each array, the top colored image is the aligned waveforms. The lower image is the stacked data, while the arrow indicates the first pulse of P wave. (d) shows the frequency spectrum of three records from three respectively as an example.

由于多臺陣壓縮傳感方法是基于反演的算法，所以數(shù)據(jù)中的噪聲和一些反射波的震相(例如pP)可能使得反演的結(jié)果中有一些假象，因此就非常有必要在通過反演直接求出的所有解中挑選出真實(shí)的解.為了達(dá)到判定解是否真實(shí)可信的目的，設(shè)定了一個判據(jù).首先，定義：

(6)

這里，對于第q個時間窗，Pm q是按由大到小排列的第m個區(qū)域極大值，Mq是相對殘差，定義為：

(7)

本文取m=1，也就是每個時間窗所求得的解只取極大值作為待選.通過公式(6)，可以計(jì)算出α隨時間窗變化的圖像(圖3).然后，再定義一個鑒別因子(DF).這意味著，所有α值小于DF的時間窗的解將會被認(rèn)為是假象而舍去.在本文研究中，DF=0.03，這與所使用的數(shù)據(jù)的信噪比(SNR=30)的倒數(shù)是幾乎一致的.但是由于在震源開始的破裂階段能量釋放較小，因此使用此判據(jù)會丟失掉解中初期的信息，所以此判據(jù)對于剛開始破裂的前7s并不適用.也就是說，雖然前7s的解的α小于DF，但是依然保留作為最后的解.圖4和圖5分別給出的是篩選前和篩選后的結(jié)果.可以看到，經(jīng)過篩選，一些能量相對較小，殘差相對較大的解被刪去，這有助于我們更精確地估計(jì)震源破裂持續(xù)時間和斷層破裂的尺度，獲取破裂特征信息.

圖3 藍(lán)線代表α隨時間變化的曲線，即不同時間窗中的解的α的變化.紅線代表本次選定的鑒別因子的值,這意味著除前7 s外紅線下方的解是被舍棄的.Fig.3 α versus rupture time. The red line is the value of the chosen DF which means that the results below the red line are removed except the results of first 7 s.

圖4 對結(jié)果進(jìn)行篩選前的地震破裂過程黑色五角星代表本次地震的震中位置，彩色圓圈代表能量釋放位置，顏色代表時間，圓圈的大小代表釋放能量的相對大小.Fig.4 The raw rupture process of 2015 Nepal earthquake before filtering Black star is the epicenter of the 2015 earthquake. Colored circles denote the location of power radiation. Color means rupture time. The radius of circle is proportional to the relative value of power radiation.

圖5 篩選后的2015尼泊爾地震的最終破裂過程黑色五角星代表本次地震的震中位置，彩色圓圈代表能量釋放位置，顏色代表時間，圓圈的大小代表釋放能量的相對大小.灰色圓圈代表震后5天內(nèi)的余震分布.紅色五角星代表1833年和1934年的震中位置.虛框的橢圓為2015和1934兩次地震的估計(jì)斷層破裂范圍.Fig.5 The final rupture process of 2015 Nepal earthquake after filtering Black star is the epicenter of the 2015 earthquake. Colored circles denote the location of power radiation. Color means rupture time. The radius of circle is proportional to the relative value of power radiation. Grey circle denotes aftershocks within 5 days after the main shock. Red stars are the location of 1833 and 1934 Nepal earthquake. Dashed ellipses are the estimated rupture zones of 2015 and 1934 earthquake.

3 結(jié)果與討論

經(jīng)過篩選，最終得到了這次地震的破裂過程(圖5).結(jié)果顯示，本次地震從震中開始，向南東東方向破裂約105km.其中，破裂分為兩個階段.第一個階段是地震發(fā)生后的前15s，破裂所釋放的能量基本圍繞在震中附近10km以內(nèi)，并沒有向遠(yuǎn)處傳播.第二個階段從第16s開始，破裂從震中東部約40km處開始向南東東方向傳播，總的破裂持續(xù)時間為59s.值得注意的是，破裂到斷層最東邊的時間為51s，此時破裂尚未結(jié)束.在破裂的最后階段，破裂能量集中的位置又回到了此次地震能量釋放最強(qiáng)的地方(28°N，85.4°E).這意味著這個位置從29s開始破裂，直到59s才結(jié)束，持續(xù)時間長達(dá)30s，持續(xù)釋放了較長時間的能量.并且可以判斷此處應(yīng)該是受災(zāi)較為嚴(yán)重，釋放能量較大的區(qū)域.

從上述解中提取出破裂前沿的位置，并沿?cái)鄬幼呦蛲队埃梢赃M(jìn)一步得到在破裂離開震中附近到達(dá)震中東部40km處之后，破裂向南東東方向以約1.9km·s-1的速度傳播(圖6)，到達(dá)距離震中105km處終止.

圖6 破裂前沿與震源之間的距離隨破裂前沿傳播的時間的變化黑線給出的是線性擬合出的第二階段中的破裂速度，為1.9 km·s-1.Fig.6 The distance from the rupture front to the epicenter versus rupture time The black line is obtained by linear regression to fit time-distance relations of rupture process. Its slope is the mean rupture velocity with a value of 1.9 km·s-1.

對比震級極其相近且也是發(fā)生在陸地板塊間的2008年汶川地震，可以發(fā)現(xiàn)兩次具有巨大破壞力的地震有著明顯的幾個相同點(diǎn)：(1)破裂都為單側(cè)破裂:汶川地震的破裂方向是從震源開始向東北方向沿龍門山斷裂帶發(fā)展(張勇等，2008；ZhangandGe，2010)，尼泊爾地震是向東南東方向沿印度板塊與歐亞大陸板塊的俯沖帶走向的反方向上發(fā)展，均沒有向另一方向破裂的能量釋放.(2)兩次地震都存在明顯的階段性: 汶川地震在初始階段的0～20s，看不出明顯的破裂方向(Zhangetal.，2014)，尼泊爾地震前15s破裂也都圍繞在震源附近.兩次地震都是在初始的小部分能量在震源附近區(qū)域釋放后的一段時間內(nèi)才表現(xiàn)出明顯的破裂方向性和主要能量釋放特征.(3)能量釋放最大的位置均不在震中: 汶川地震能量釋放的最大位置在震源東北方向約50km處(ZhangandGe，2010)，尼泊爾地震在距震中70km處.然而，這兩次地震也有很明顯的不同點(diǎn)：(1)破裂速度和尺度：汶川地震的破裂速度約為3km·s-1，破裂尺度約為300km(張勇等，2008；Xuetal.， 2009；ZhangandGe，2010)；而尼泊爾地震的第二階段有明顯方向性的破裂速度為1.9km·s-1，長度為105km.經(jīng)過對比可以發(fā)現(xiàn)，尼泊爾地震無論是破裂尺度還是破裂速度均要小于汶川地震.這表明尼泊爾地震的能量釋放更加集中.(2)地質(zhì)資料和地震形成的地表破裂軌跡表明，汶川地震的破裂并不是發(fā)生在單一的斷層上，而是發(fā)生在龍門山構(gòu)造帶的映秀—北川斷裂和灌縣—江油斷裂，中間的能量釋放在空間上有一個跳躍(王衛(wèi)民等，2008).從大的地質(zhì)背景來看，汶川地震屬于歐亞大陸板塊內(nèi)部的地震.而尼泊爾地震并沒有發(fā)現(xiàn)跨斷層觸發(fā)的明顯證據(jù)，屬于板間地震.

另一方面，可以看到在破裂的第一階段和第二階段空間位置的中間存在一片沒有明顯能量分布的空區(qū).同時這一空區(qū)的位置也和本次地震后五天內(nèi)余震的分布的空區(qū)相一致(圖5).這表明該部分地區(qū)無論是震時還是在震后在高頻部分均沒有較大幅度的能量釋放，兩者分布一致，證明了利用多臺陣壓縮傳感算法反演得到的能量釋放位置是可信的.而且，根據(jù)我們的結(jié)果，可以看到能量釋放最集中的區(qū)域就在尼泊爾首都加德滿都附近，而加德滿都也是本次受災(zāi)最為嚴(yán)重的地區(qū).由于造成房屋毀壞倒塌，破壞力較強(qiáng)的地震波為相對高頻的波動，相比利用有限斷層反演的低頻的結(jié)果，高頻能量輻射位置會與破壞較為嚴(yán)重的位置更為一致.因此，利用多臺壓縮傳感方法來反演大地震的破裂過程可以快速得到能量釋放最為集中、破壞程度較為嚴(yán)重的地區(qū)的地理位置，這給震后第一時間的快速救援提供了很大的幫助，同時體現(xiàn)出了在相對較高頻率進(jìn)行破裂過程反演在這一方面具有一定的優(yōu)越性.同時，我們將本次地震在0.25～1Hz頻段釋放的能量累加之后，可得到高頻總能量釋放的空間分布(圖7).多臺陣壓縮傳感方法給出的結(jié)果與利用INSAR數(shù)據(jù)給出的破裂范圍基本相同(http:∥e-decider.org/content/nepal-earthquake-april-2015)，這證明我們有關(guān)破裂尺度的估計(jì)是可信的.以上對比證明，通過多臺陣壓縮傳感方法得到的高頻能量釋放區(qū)域的解與其他類型數(shù)據(jù)得到的結(jié)果具有一致性.

圖7 本次地震在0.25～1 Hz能量釋放的空間分布紅色五角星代表震中位置，虛框矩形為INSAR給出的斷層破裂范圍.Fig.7 The spatial distribution of power radiation in frequency band of 0.25～1 Hz The red star denotes the epicenter of 2015 Nepal earthquake. Dashed rectangle is the range of rupture observed by INSAR.

本次研究所用到的多臺陣壓縮傳感方法使用了來自三個臺陣的數(shù)據(jù)，目的是減輕由于臺陣方位角和孔徑的局限性而導(dǎo)致的對于破裂細(xì)節(jié)的分辨率下降的問題.我們通過計(jì)算1Hz時的臺陣響應(yīng)函數(shù)(ARF)(RostandThomas，2002；Xuetal.，2009)來判斷這樣的臺站分布是否可以有效提升分辨率.由圖8可以發(fā)現(xiàn)，在選取同樣的置信區(qū)間的時候(ARF>0.7)，使用多臺陣分布可以使解的范圍更加精確，等值包絡(luò)線范圍最小，而選取其他單一臺陣的時候的精度較低，包絡(luò)范圍較大.同時，我們發(fā)現(xiàn)，若只使用單一臺陣，例如阿拉斯加或歐洲臺陣的數(shù)據(jù)，橢圓的傾斜方向即長軸方向與破裂方向較為一致，因此可能對于本次地震的破裂細(xì)節(jié)的反演精度較低，在破裂方向上分辨率較差.這表明臺陣中臺站的分布以及臺陣的方位角都可能影響到最終解的精確度.值得注意的是，在使用多臺陣數(shù)據(jù)計(jì)算時，雖然解的精度有較為明顯的提升，但是臺陣響應(yīng)函數(shù)的結(jié)果顯示其分布并非是一個圓形，而在南北方向可能出現(xiàn)一些誤差，但是考慮到此次地震的破裂方向?yàn)槟蠔|東，因此對于整體破裂信息的獲取影響不大.

圖8 各個臺陣在1 Hz的臺陣響應(yīng)函數(shù)(ARF)大于0.7的包絡(luò)線紫色五角星為震中位置.箭頭代表指向臺陣的方向.Fig.8 The envelopes of array response functions (ARF) larger than 0.7 of three arrays, and their sum, at 1 Hz The purple star denotes the epicenter. The black arrow represents the direction toward array.

本次地震發(fā)生的位置是在從南向北俯沖的印度板塊和在北邊上覆的歐亞板塊的俯沖帶上，這一區(qū)所處的斷裂帶被稱為主喜馬拉雅逆沖斷裂(MHT).通過GPS數(shù)據(jù)計(jì)算可以得到在MHT上，從地表到沿傾角向下的約100km是鎖定的區(qū)域，近20年的觀測也表明這一地區(qū)幾乎沒有慢滑移事件，也就是說累積的能量都需要通過地震來釋放(Aderetal.，2012).作為兩個板塊的邊界，歷史上這個地方也曾發(fā)生過一些大地震.在本次地震的發(fā)生區(qū)域的方圓250km內(nèi)曾經(jīng)在過去的一個世紀(jì)發(fā)生過4次6級以上的大地震.最大的一次是震級大于M8.0，發(fā)生在1934年(Hough and Roger，2008；Sapkota et al.，2013；Bollinger et al.，2014)，斷層破裂范圍至少從85°52′E至87°17′E，最小破裂長度為150 km(圖5).從我們的結(jié)果可以看出，2015地震破裂以及余震活動的分部剛好延伸至1934地震的破裂帶的西緣，兩個破裂帶緊鄰，所以2015地震可以被看作是1934地震所處斷層的另一區(qū)域的能量釋放.在更早的19世紀(jì)，1833年也曾有一次MW7.6的地震(Ambraseys and Douglas，2004)發(fā)生在與2015年地震相同的發(fā)震斷層(圖5).為了研究這兩次地震之間的關(guān)系，我們發(fā)現(xiàn)由GPS數(shù)據(jù)得到的此區(qū)域的板塊匯聚速率為17.8±0.5 mm·a-1，但是由于MHT淺部的鎖定，其滑移速率極慢(<5 mm·a-1)，所以這一區(qū)域的滑移缺損累積速率至少為13 mm·a-1(Bilham and Ambraseys，2005)；MHT在尼泊爾地區(qū)的震間的地震矩缺損累積速率為6.6±0.4×1019Nm·a-1(Aderetal.，2012).那么根據(jù)USGS的有限斷層反演結(jié)果(http:∥earthquake.usgs.gov/earthquakes/eventpage/us20002926scientific_finitefault)得到本次地震的最大滑移量約為3.3m，平均滑移量約為2m，標(biāo)量地震矩約為0.8×1021N·m.可以推斷出：本次地震的平均滑移量與自1833年積累至今的180年的缺損所積攢的滑移缺損總量幾乎相等；但本次地震的標(biāo)量地震矩僅為地震矩缺損累積總量的約十分之一，可以認(rèn)為地震矩缺損已經(jīng)釋放了一小部分.因此從該斷層破裂區(qū)域的小尺度范圍看，我們認(rèn)為本次地震可以被看作是1833年地震后，局部能量積累的再次釋放.但是由于尼泊爾境內(nèi)的MHT積累的地震矩缺損絕大部分能量還沒有釋放，因此整個尼泊爾地區(qū)在未來一段時間仍不排除有繼續(xù)發(fā)生更大地震的可能.

4 結(jié)論

本文利用臺陣壓縮傳感方法給出了2015尼泊爾MW7.9地震的破裂過程和破裂細(xì)節(jié). 結(jié)果表明此次破裂是一個典型的單側(cè)破裂，破裂方向?yàn)槟蠔|東, 破裂尺度約為105 km, 破裂傳播速度為1.9 km·s-1，破裂累積能量最大值位于尼泊爾首都加德滿都附近.都是造成加德滿都震害嚴(yán)重的原因.

通過與震級相近的2008汶川地震的細(xì)節(jié)破裂特征進(jìn)行比較，我們發(fā)現(xiàn)兩次地震都是單側(cè)破裂、存在明顯的分段特征且最大能量釋放位置并非震中；但是，與汶川地震相比尼泊爾地震的破裂尺度更短，說明其能量釋放更加集中.

根據(jù)歷史上曾經(jīng)發(fā)生在尼泊爾境內(nèi)的主喜馬拉雅逆沖斷裂(MHT)上的地震以及板塊的相對運(yùn)動速率，推測本次2015尼泊爾地震可能是繼1833年MW7.6地震后，在該斷層上又一次小尺度的180年間積累的能量釋放的過程，不排除整個尼泊爾地區(qū)發(fā)生更大地震的可能性.

致謝 本文地震波數(shù)據(jù)來自IRIS 數(shù)據(jù)中心，部分圖形的繪制使用的是GMT繪圖軟件. 兩位匿名審稿專家提出了非常中肯的修改意見，作者對他們致以誠摯的感謝！

Ader T, Avouac J P, Zeng J L, et al. 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: Implications for seismic hazard.JournalofGeophysicalResearch:SolidEarth(1978—2012), 117(B4): B04403.

Aki K, Richards P G. 2002. Quantitative Seismology (Vol.1): Theory and Methods. Sansalito, CA: University Science Books.Ambraseys N N, Douglas J. 2004. Magnitude calibration of north Indian earthquakes.GeophysicalJournalInternational, 159(1): 165-206.

Bilham R, Ambraseys N. 2005. Apparent Himalayan slip deficit from the summation of seismic moments for Himalayan earthquakes, 1500-2000.CurrentScience, 88(10): 1658-1663. Bollinger L, Sapkota S N, Tapponnier P, et al. 2014. Estimating the return times of great Himalayan earthquakes in eastern Nepal: Evidence from the Patu and Bardibas strands of the Main Frontal Thrust.JournalofGeophysicalResearch:SolidEarth, 119(9): 7123-7163.

Candès E J, Romberg J, Tao T. 2006. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information.IEEETransactionsonInformationTheory, 52(2): 489-509.

Donoho D L. 2006. Compressed sensing.IEEETransactionsonInformationTheory, 52(4): 1289-1306.

Duarte M F, Davenport M A, Takhar D, et al. 2008. Single-pixel imaging via compressive sampling.IEEESignalProcessingMagazine, 25(2): 83-91.

Dziewonski A M, Anderson D L. 1981. Preliminary reference Earth model.PhysicsoftheEarthandPlanetaryInteriors, 25(4): 297-356.

Goldstein P, Archuleta R J. 1991. Deterministic frequency-wavenumber methods and direct measurements of rupture propagation during earthquakes using a dense array: Theory and methods.JournalofGeophysicalResearch:SolidEarth(1978—2012), 96(B4): 6173-6185.

Hough S E, Roger B. 2008. Site response of the Ganges basin inferred from re-evaluated macroseismic observations from the 1897 Shillong, 1905 Kangra, and 1934 Nepal earthquakes.JournalofEarthSystemScience, 117(2S): 773-782.

Ishii M, Shearer P M, Houston H, et al. 2005. Extent, duration and speed of the 2004 Sumatra—Andaman earthquake imaged by the Hi-Net array.Nature, 435(7044): 933-936.

Koper K D, Hutko A R, Lay T, et al. 2012. Imaging short-period seismic radiation from the 27 February 2010 Chile (Mw8.8) earthquake by back-projection of P, PP, and PKIKP waves.JournalofGeophysicalResearch:SolidEarth(1978—2012), 117(B2): B02308.

Lustig M, Donoho D, Pauly J M. 2007. Sparse MRI: The application of compressed sensing for rapid MR imaging.MagneticResonanceinMedicine, 58(6): 1182-1195.

Malioutov D, ?etin M, Willsky A S. 2005. A sparse signal reconstruction perspective for source localization with sensor arrays.IEEETransactionsonSignalProcessing, 53(8): 3010-3022.

Meng L S, Inbal A, Ampuero J P. 2011. A window into the complexity of the dynamic rupture of the 2011Mw9 Tohoku-Oki earthquake.GeophysicalResearchLetters, 38(7): L00G07.Meng L, Ampuero J P, Stock J, et al. 2012. Earthquake in a maze: Compressional rupture branching during the 2012Mw8.6 Sumatra earthquake.Science, 337(6095): 724-726. Rost S, Thomas C. 2002. Array seismology: methods and applications.ReviewsofGeophysics, 40(3): 2-1-2-27.Sapkota S N, Bollinger L, Klinger Y, et al. 2013. Primary surface ruptures of the great Himalayan earthquakes in 1934 and 1255.NatureGeoscience, 6(1): 71-76.

Schmidt R O. 1986. Multiple emitter location and signal parameter estimation.IEEETransactionsonAntennasandPropagation, 34(3): 276-280.

Wang W M, Zhao L F, Li J, et al. 2008. Rupture process of theMS8.0 Wenchuan earthquake of Sichuan, China.ChineseJ.Geophys. (in Chinese), 51(5): 1403-1410, doi: 10.3321/j.issn:0001-5733.2008.05.013.

Xu Y, Koper K D, Sufri O, et al. 2009. Rupture imaging of theMw7.9 12 May 2008 Wenchuan earthquake from back projection of teleseismic P waves.Geochemistry,Geophysics,Geosystems, 10

(4): Q04006.Yao H J, Gerstoft P, Shearer P M, et al. 2011. Compressive sensing of the Tohoku-OkiMw9.0 earthquake: Frequency-dependent rupture modes.GeophysicalResearchLetters, 38(20): L20310.

Zhang H, Ge Z X. 2010. Tracking the rupture of the 2008 Wenchuan earthquake by using the relative back-projection method.BulletinoftheSeismologicalSocietyofAmerica, 100(5B): 2551-2560.

Zhang Y, Feng W P, Xu L S, et al. 2009. Spatiotemporal rupture process of the 2008 Wenchuan great earthquake.ScienceinChina(SeriesD:EarthSciences), 52(2): 145-154.

Zhang Y, Wang R J, Zschau J, et al. 2014. Automatic imaging of earthquake rupture processes by iterative deconvolution and stacking of high-rate GPS and strong motion seismograms.JournalofGeophysicalResearch:SolidEarth, 119(7): 5633-5650.

附中文參考文獻(xiàn)

王衛(wèi)民, 趙連鋒, 李娟等. 2008. 四川汶川8.0級地震震源過程. 地球物理學(xué)報(bào), 51(5): 1403-1410, doi: 10.3321/j.issn:0001-5733.2008.05.013.

張勇, 馮萬鵬, 許力生等. 2008. 年汶川大地震的時空破裂過程. 中國科學(xué)(D輯), 38(10): 1186-1194.

(本文編輯 何燕)

Rupturing process of theMW7.9 Nepal earthquake inverted by the multi-array compressive sensing method

LIU Zhi-Peng， GE Zeng-Xi*

DepartmentofGeophysics,SchoolofEarthandSpaceSciences,PekingUniversity,Beijing100871，China

AnMW7.9 earthquake occurred in Nepal on 25 April, 2015. Its hypocenter is located at 28.1°N, 84.7°E. In order to investigate the rupturing process of this temblor in detail, we use the multi-array compressive sensing (MCS) method to invert vertical component data of teleseismic P waves recorded by 179 stations of Alaska, Europe and Australia arrays.Traditionally, the compressive sensing (CS) method, based on data from unique array, can recover sparsely distributed source locations in the frequency domain. While the MCS could improve the resolution because of the larger azimuth coverage range than only using data from single array. To apply this method, the first thing is to align the data from each array using cross-correlation respectively. Secondly we stack the aligned data from each array to make the first P wave pulse clear enough to be recognized so that all the data could be aligned. Finally we invert the spatiotemporal rupture process with a sliding window based on the multi-array compressive sensing method. Combining with the discriminant factor (DF), we reduce artifacts and obtain a more credible result.The result indicates that there is a unilateral rupture with orientation of SEE, length of 105 km, and duration time of 58 s. During the first 15 s of failure, the power radiation is mainly distributed around the epicenter. Since 16 s, the rupture propagated toward SEE with a rate of 1.9 km·s-1. The largest energy is released at 38 s, 70 km from the hypocenter. The break of this position starts from 29 s and lasts for 30 s.Additionally, we compare the rupture details of the Nepal earthquake and 2008 Wenchuan earthquake because of their same magnitude. Moreover, MCS can offer zones where power radiated intensely, which could provide helpful information for emergency rescue after earthquakes. And we discuss the improvement of resolution because of the larger azimuth coverage range. In the end, based on slip deficit and rate of accumulation of moment deficit, we conclude that this 2015 earthquake could be a repeat of the 1833M7.7 Nepal temblor. And the possibility of the occurrence of a larger earthquake in this region cannot be excluded.

Multi-array compressive sensing; Nepal earthquake; Rupturing process

10.6038/cjg20150605.

國家自然科學(xué)基金項(xiàng)目(41374045)資助.

劉志鵬，男，在讀博士生，主要從事震源反演方面的研究. E-mail: liuzhipeng@pku.edu.cn

*通訊作者 蓋增喜，男，副教授，1996年畢業(yè)于北京大學(xué)，2005年獲北京大學(xué)博士學(xué)位.主要從事地震學(xué)方面的教學(xué)和科研工作. E-mail:zge@pku.edu.cn

10.6038/cjg20150605

P315

2015-05-17，2015-06-09收修定稿

劉志鵬，蓋增喜. 2015. 利用多臺陣壓縮傳感方法反演尼泊爾MW7.9地震破裂過程.地球物理學(xué)報(bào)，58(6):1891-1899,

Liu Z P, Ge Z X. 2015. Rupturing process of theMW7.9 Nepal earthquake inverted by the multi-array compressive sensing method.ChineseJ.Geophys. (in Chinese),58(6):1891-1899,doi:10.6038/cjg20150605.