利用多臺(tái)陣壓縮傳感方法反演尼泊爾MW7.9地震破裂過(guò)程

劉志鵬， 蓋增喜

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

?

利用多臺(tái)陣壓縮傳感方法反演尼泊爾MW7.9地震破裂過(guò)程

劉志鵬， 蓋增喜*

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

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

多臺(tái)陣壓縮傳感； 尼泊爾地震； 地震破裂過(guò)程

1 引言

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

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

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

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

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

(1)

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

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

(2)

寫(xiě)成離散的形式，則有

(3)

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

(4)

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

+λ‖x‖1].

(5)

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

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

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

圖1 尼泊爾地震的震中位置以及所使用的臺(tái)陣數(shù)據(jù)的分布黑色五角星代表震中的位置，藍(lán)色圓形代表臺(tái)站.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)澳大利亞三個(gè)臺(tái)陣上方的彩色圖案代表地震波形數(shù)據(jù)及各臺(tái)陣以P波初值對(duì)齊后疊加的地震圖，箭頭標(biāo)示P波初動(dòng)；(d)圖為臺(tái)網(wǎng)中代表性臺(tái)站的數(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.

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

(6)

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

(7)

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

圖3 藍(lán)線代表α隨時(shí)間變化的曲線，即不同時(shí)間窗中的解的α的變化.紅線代表本次選定的鑒別因子的值,這意味著除前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 對(duì)結(jié)果進(jìn)行篩選前的地震破裂過(guò)程黑色五角星代表本次地震的震中位置，彩色圓圈代表能量釋放位置，顏色代表時(shí)間，圓圈的大小代表釋放能量的相對(duì)大小.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尼泊爾地震的最終破裂過(guò)程黑色五角星代表本次地震的震中位置，彩色圓圈代表能量釋放位置，顏色代表時(shí)間，圓圈的大小代表釋放能量的相對(duì)大小.灰色圓圈代表震后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)過(guò)篩選，最終得到了這次地震的破裂過(guò)程(圖5).結(jié)果顯示，本次地震從震中開(kāi)始，向南東東方向破裂約105km.其中，破裂分為兩個(gè)階段.第一個(gè)階段是地震發(fā)生后的前15s，破裂所釋放的能量基本圍繞在震中附近10km以內(nèi)，并沒(méi)有向遠(yuǎn)處傳播.第二個(gè)階段從第16s開(kāi)始，破裂從震中東部約40km處開(kāi)始向南東東方向傳播，總的破裂持續(xù)時(shí)間為59s.值得注意的是，破裂到斷層最東邊的時(shí)間為51s，此時(shí)破裂尚未結(jié)束.在破裂的最后階段，破裂能量集中的位置又回到了此次地震能量釋放最強(qiáng)的地方(28°N，85.4°E).這意味著這個(gè)位置從29s開(kāi)始破裂，直到59s才結(jié)束，持續(xù)時(shí)間長(zhǎng)達(dá)30s，持續(xù)釋放了較長(zhǎng)時(shí)間的能量.并且可以判斷此處應(yīng)該是受災(zāi)較為嚴(yán)重，釋放能量較大的區(qū)域.

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

圖6 破裂前沿與震源之間的距離隨破裂前沿傳播的時(shí)間的變化黑線給出的是線性擬合出的第二階段中的破裂速度，為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.

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

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

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

圖8 各個(gè)臺(tái)陣在1 Hz的臺(tái)陣響應(yīng)函數(shù)(ARF)大于0.7的包絡(luò)線紫色五角星為震中位置.箭頭代表指向臺(tái)陣的方向.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).通過(guò)GPS數(shù)據(jù)計(jì)算可以得到在MHT上，從地表到沿傾角向下的約100km是鎖定的區(qū)域，近20年的觀測(cè)也表明這一地區(qū)幾乎沒(méi)有慢滑移事件，也就是說(shuō)累積的能量都需要通過(guò)地震來(lái)釋放(Aderetal.，2012).作為兩個(gè)板塊的邊界，歷史上這個(gè)地方也曾發(fā)生過(guò)一些大地震.在本次地震的發(fā)生區(qū)域的方圓250km內(nèi)曾經(jīng)在過(guò)去的一個(gè)世紀(jì)發(fā)生過(guò)4次6級(jí)以上的大地震.最大的一次是震級(jí)大于M8.0，發(fā)生在1934年(Hough and Roger，2008；Sapkota et al.，2013；Bollinger et al.，2014)，斷層破裂范圍至少?gòu)?5°52′E至87°17′E，最小破裂長(zhǎng)度為150 km(圖5).從我們的結(jié)果可以看出，2015地震破裂以及余震活動(dòng)的分部剛好延伸至1934地震的破裂帶的西緣，兩個(gè)破裂帶緊鄰，所以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積累的地震矩缺損絕大部分能量還沒(méi)有釋放，因此整個(gè)尼泊爾地區(qū)在未來(lái)一段時(shí)間仍不排除有繼續(xù)發(fā)生更大地震的可能.

4 結(jié)論

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

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

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

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

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級(jí)地震震源過(guò)程. 地球物理學(xué)報(bào), 51(5): 1403-1410, doi: 10.3321/j.issn:0001-5733.2008.05.013.

張勇, 馮萬(wàn)鵬, 許力生等. 2008. 年汶川大地震的時(shí)空破裂過(guò)程. 中國(guó)科學(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.

國(guó)家自然科學(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. 利用多臺(tái)陣壓縮傳感方法反演尼泊爾MW7.9地震破裂過(guò)程.地球物理學(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.