董仁東 姚強(qiáng) 施賀青 王琪

1)中國地質(zhì)大學(xué)(武漢)，地球物理與空間信息學(xué)院，武漢 430074 2)中國地質(zhì)大學(xué)(武漢)，地球科學(xué)學(xué)院，武漢 430074 3)中國地震局地震研究所，武漢 430071

0 引言

2001年11月14日9時(shí)26分10秒，青海昆侖山口西發(fā)生MS8.1 地震，該地震是自三分量寬頻帶數(shù)字記錄以來，青藏高原內(nèi)部發(fā)生的最大地震。美國地質(zhì)調(diào)查局(USGS)給出的震中位置為35.946°N，90.541°E，矩震級(jí)MW為7.9，該特大地震導(dǎo)致昆侖山斷層西段地表大范圍、高強(qiáng)度破裂。昆侖山斷裂是一條近東西走向、距離長達(dá)數(shù)千千米以上的大型左旋走滑斷層，是調(diào)節(jié)青藏高原受印度板塊擠壓而造成物質(zhì)東向滑移的關(guān)鍵部位(Tapponnier et al，2001；Wang et al，2001； 王敏等，2020)。歷史上該斷裂大震頻繁，如1963年MS7.0 都蘭地震、1973年MS7.3 和1997年MS7.9 瑪尼地震(Triep et al，1997； 劉剛等，2020)分別位于該地震東西兩側(cè)，2001年昆侖山口西地震則填補(bǔ)了昆侖山斷裂帶西段的部分地震空區(qū)。

Xu等(2002)的地質(zhì)調(diào)查結(jié)果表明，此次地震破裂長度約400km，產(chǎn)生了最大達(dá)7～8m的地表破裂。喬學(xué)軍等(2002)首次公布GPS同震觀測結(jié)果，其中距斷層最近的點(diǎn)位位移超過1m。任金衛(wèi)等(2005)系統(tǒng)地分析震后GPS觀測資料，給出了該地震同震變形場及震后位移時(shí)間序列。Lasserre等(2005)通過獲取橫跨斷層兩側(cè)的4幅InSAR影像數(shù)據(jù)，更為準(zhǔn)確地圈定地震破裂及影響范圍。綜合上述觀測資料顯示，地震自西向東導(dǎo)致太陽湖、庫塞湖、昆侖山口3段集中破裂區(qū)。

Bouchon等(2003)利用近場區(qū)域?qū)掝l帶數(shù)據(jù)估算地震破裂的平均速度為3.7～3.9km/s，并探測到其最大速度達(dá) 5km/s，超過當(dāng)?shù)丶羟胁ㄋ俣?。Antolik等(2004)利用遠(yuǎn)震體波反演地震破裂時(shí)空過程，揭示此次地震可分為3個(gè)子事件，初次子事件發(fā)生在拉張型的走滑地塹上，后續(xù)兩次為走滑型斷層破裂，總破裂持續(xù)時(shí)間約120s； 當(dāng)破裂傳播到達(dá)主斷層時(shí)，平均破裂速度達(dá)到 3.6km/s。Robinson等(2006)利用遠(yuǎn)場SH波資料，顯示破裂主要集中在20km以上脆性地殼內(nèi)，指出最大破裂速度、最大應(yīng)力降、最大滑移量之間具有一致性，且昆侖山斷層在西大灘段的分叉導(dǎo)致超剪切破裂迅速終止。

不同學(xué)者曾給出一些昆侖山口西地震的破裂模型。近年來的代表性成果有屠泓為等(2016)利用GPS/InSAR反演的破裂分布模型及Wen等(2009)利用遠(yuǎn)震地震波與地表破裂資料聯(lián)合反演的破裂時(shí)空過程，但迄今尚無綜合所有大地測量及地震觀測資料的聯(lián)合反演模型。鑒于利用不同資料、不同算法構(gòu)建的模型間存在一定的差異，為了深入、細(xì)致地揭示此次地震的發(fā)震過程，本研究融合全球遠(yuǎn)震動(dòng)態(tài)波形、近場GPS和InSAR靜態(tài)位移數(shù)據(jù)，基于有限斷層位錯(cuò)理論和模擬退火算法開展聯(lián)合反演，構(gòu)建統(tǒng)一的破裂模型； 同時(shí)利用反投影方法，不依賴斷層幾何模型和地球內(nèi)部結(jié)構(gòu)模型等先驗(yàn)信息，直接確定地震破裂過程，分析破裂速度變化與有限斷層反演破裂空間分布的內(nèi)在關(guān)系。

1 資料與方法

1.1 有限斷層反演數(shù)據(jù)

本研究使用美國地震學(xué)研究聯(lián)合會(huì)(IRIS)數(shù)據(jù)中心(1)http：//www.iris.washington.edu/wilber3震中距為30°～90°的遠(yuǎn)震體波數(shù)據(jù)，該區(qū)域數(shù)據(jù)無其他波形干擾，依據(jù)信噪比和臺(tái)站空間分布最終選取35個(gè)P波臺(tái)站和33個(gè)SH波臺(tái)站(圖1)?？紤]到地震持續(xù)時(shí)間，本研究拾取P波、SH波初至前30s至初至后130s，共計(jì)160s的時(shí)程數(shù)據(jù)，通過去除儀器響應(yīng)和帶通濾波(0.0033～0.5Hz)，獲得用于模型反演的速度波形。

注： (a)區(qū)域構(gòu)造背景； 紅色五角星表示震中，黃色虛線方框?yàn)檠芯繀^(qū)域，綠色虛線方框?yàn)镮nSAR覆蓋區(qū)域，藍(lán)色和紅色箭頭分別表示同震位移觀測值和模擬值，插圖中灰色實(shí)線為研究區(qū)域的斷層分布，沙灘球?yàn)槔錾娇谖鞯卣鸷蜌v史地震空間分布及震源機(jī)制解，五角星為震后半年余震分布，最大震級(jí)為MW5.6，紅色實(shí)線為聯(lián)合反演采用的斷層幾何模型沿地表的跡線； (b)遠(yuǎn)震體波臺(tái)站分布； 藍(lán)色三角形表示P波臺(tái)站、紅色三角形表示SH波臺(tái)站、黃色三角形表示P波、SH波臺(tái)站； (c)淺藍(lán)色三角形表示所選用的歐洲區(qū)域臺(tái)站； 紫色三角形表示參考站圖1 2001年昆侖山口西地震區(qū)域構(gòu)造背景和臺(tái)站分布

在變形數(shù)據(jù)方面，本研究使用Lasserre等(2005)處理的4幅ERS-1/2衛(wèi)星SAR干涉影像數(shù)據(jù)，從中提取了4468個(gè)樣點(diǎn)(樣點(diǎn)最大視線向位移121cm)，完全覆蓋破裂區(qū)域； 另外，從任金衛(wèi)等(2005)公布的GPS位移數(shù)據(jù)中選取34個(gè)近場GPS站點(diǎn)。GPS站點(diǎn)雖然對(duì)破裂區(qū)域覆蓋遠(yuǎn)不及InSAR均勻、完整，但最靠近斷層的XZ02和BS33兩點(diǎn)位于斷層南北兩側(cè)，均有約1m的同震位移，方向相反，清楚地顯示了本次地震導(dǎo)致的左旋走滑破裂。由于近場資料時(shí)效性較弱，測量時(shí)間跨度達(dá)數(shù)年，故包含有地震快速震后變形信息。

1.2 斷層幾何形態(tài)及有限斷層反演

本研究幾何模型基于Xu等(2002)的地質(zhì)調(diào)查結(jié)果，并按萬永革等(2008)的方式將斷層自西向東分為3段。前兩段主要考慮彎曲的斷層幾何結(jié)構(gòu)及拉張型破裂，第一段位于破裂最西端的太陽湖，長40km，走向?yàn)?05°； 第二段位于布格達(dá)坂，長48km，走向?yàn)?6°； 第三段為純走滑的主破裂段，長360km，走向?yàn)?7°。各段的傾角均為81°，寬度均為32km，較萬永革等(2008)的處理略有不同。將各段按8km×4km網(wǎng)格大小分割為子斷層，估算每個(gè)子斷層的滑動(dòng)量、滑動(dòng)角、上升時(shí)間及平均破裂速度(Ji et al，2002、2003)，使用公式為

(1)

為求解各類參數(shù)，本研究采用模擬退火方法，利用前文中的遠(yuǎn)震波形和近場變形數(shù)據(jù)，在參數(shù)空間內(nèi)按一定步長搜尋最優(yōu)值，使得建模目標(biāo)函數(shù)最小。目標(biāo)函數(shù)由觀測值的誤差函數(shù)、權(quán)重比和模型參數(shù)約束函數(shù)構(gòu)成。對(duì)于靜態(tài)位移數(shù)據(jù)，采用如下誤差函數(shù)

(2)

對(duì)于遠(yuǎn)震體波數(shù)據(jù)，在小波域內(nèi)定義如下誤差函數(shù)

errwf=el+eh

(3)

其中，el，eh分別為低、高頻波形數(shù)據(jù)的誤差函數(shù)，具體形式上與靜態(tài)位移誤差函數(shù)類似(Ji et al，2002)。

為使反演結(jié)果穩(wěn)定且滿足一定的物理?xiàng)l件，反演計(jì)算設(shè)定針對(duì)模型參數(shù)的約束函數(shù)，強(qiáng)制模型參數(shù)滿足如下條件：①相鄰子斷層上滑移幅度差值應(yīng)盡量小(光滑約束Scon)；②由滑移量估算的總地震矩逼近某個(gè)特定值(震級(jí)約束Wcon)，實(shí)現(xiàn)參數(shù)解算的正則化。

綜上，動(dòng)態(tài)、靜態(tài)觀測值誤差函數(shù)及模型參數(shù)約束函數(shù)的建模目標(biāo)函數(shù)為

err=errwf+Wst*errst+γScon+Wcon

(4)

其中，Wst為動(dòng)態(tài)波形與靜態(tài)位移數(shù)據(jù)間的權(quán)重比；γ為平滑因子，是控制子斷層滑移幅度及其分布的關(guān)鍵參數(shù)，一般通過試錯(cuò)的方式來選取。本研究選取平滑因子時(shí)，重點(diǎn)參考地表破裂的地質(zhì)考察結(jié)果，使模型在地表處的預(yù)測值與地質(zhì)結(jié)果一致。

將模型的初始破裂點(diǎn)設(shè)定在USGS確定的震中位置，模型的總地震矩接近GCMT結(jié)果(5.9×1020Nm)。依據(jù)野外地質(zhì)調(diào)查結(jié)果，模型的滑動(dòng)量大小限定在0～9m范圍，搜索步長為30cm； 滑動(dòng)角范圍為-30°～30°(步長2°)； 平均破裂速度變化范圍為2.0～4.0km/s(步長 0.1km/s)； 上升時(shí)間變化范圍為2～10s； 搜索間隔為1s。本研究采用Galvé等(2002)的波速模型計(jì)算格林函數(shù)，青藏高原北部的巖石圈厚度在50～65km。相較于P波，SH波對(duì)破裂速度具有更高的敏感性，在波形誤差函數(shù)中將P波與SH波的權(quán)重比設(shè)為1︰1。GPS數(shù)據(jù)點(diǎn)數(shù)較少，與InSAR數(shù)據(jù)點(diǎn)數(shù)相差約100倍，為避免反演模型過分依賴單一資料，故將位移誤差函數(shù)中GPS與InSAR權(quán)重比設(shè)為100︰1。此外，本研究將目標(biāo)函數(shù)中Wst設(shè)為1，等權(quán)對(duì)待遠(yuǎn)場和近場觀測。

1.3 反投影方法

本研究將反投影方法應(yīng)用于歐洲區(qū)域的臺(tái)網(wǎng)，來求解該地震在P波高頻段的能量分布特征，使用德國地學(xué)中心(2)https：//www.gfz-potsdam.de/歐洲區(qū)域臺(tái)網(wǎng)的110個(gè)寬頻帶垂直分量數(shù)據(jù)(圖1)，將其重采樣至100Hz。臺(tái)站震中距在40°～75°之間，此段距離無其他震相對(duì)直達(dá)P波的干擾； 方位角覆蓋范圍為278°～346°。

反投影方法是將震源區(qū)網(wǎng)格化，將每個(gè)網(wǎng)格節(jié)點(diǎn)看作一次可能的破裂子事件，通過疊加波形的高頻信號(hào)來求解地震破裂時(shí)空分布。該方法不依賴高精度震中位置信息或其他先驗(yàn)信息(例如經(jīng)驗(yàn)格林函數(shù)、破裂持續(xù)時(shí)間、斷層幾何等特征)，可以用來快速地估算破裂長度、破裂持續(xù)時(shí)間及破裂速度等運(yùn)動(dòng)學(xué)參數(shù)(Ishii et al，2005；Wang et al，2017；Yao et al，2019)。

本研究在震源區(qū)域布置了3600個(gè)網(wǎng)格點(diǎn)，格點(diǎn)間距5km，深度30km。采用上述的震中和速度模型計(jì)算P波到時(shí)，通過波形的互相關(guān)運(yùn)算(時(shí)窗長度10s，滑動(dòng)長度5s，相關(guān)系數(shù)閾值0.6)，將波形對(duì)齊，并使用1s的滑動(dòng)時(shí)窗疊加波形(總時(shí)窗長度200s)，頻率范圍為0.5～2Hz，最終得到該事件的時(shí)空展布特征。

2 反演結(jié)果

圖2(a)為聯(lián)合遠(yuǎn)震和近場永久位移資料反演的昆侖山口西地震同震滑動(dòng)模型。反演結(jié)果表明，盡管昆侖山口西地震具有向西發(fā)展的趨勢，但主要表現(xiàn)單側(cè)破裂。如考慮深度15km以上、幅度大于1m的破裂，其累計(jì)長度已超過400km。反演得到的地震矩大小為6.1×1020Nm，對(duì)應(yīng)的矩震級(jí)為MW7.78，略大于GCMT給出的結(jié)果。

注： (a)中白色實(shí)線為破裂時(shí)間等值線(單位：s)，矢量箭頭大小和方向分別代表塊體破裂大小和方向，走向方向比例尺為深度方向0.5倍，紅色五角星為震中； (c)中紅色虛線為模型地表破裂值，黑色柱體為地質(zhì)考察結(jié)果(Xu et al，2002)圖2 昆侖山口西地震同震破裂模型(a)、地震矩能量釋放函數(shù)(b)和與地質(zhì)調(diào)查結(jié)果對(duì)比(c)

模型清楚地顯示：在主破裂段上3個(gè)最大破裂幅度超過6m的集中破裂區(qū)，其深度不超過10km。其中，最大破裂位于庫塞湖附近，距震中240～280km，滑動(dòng)量約8m。與地表考察結(jié)果對(duì)比，在距震中40～150km的區(qū)段模型顯示的破裂幅度(小于1m)低于現(xiàn)場實(shí)測和基于光學(xué)遙感觀測的結(jié)果(Xu et al，2006；Klinger et al，2005)。

模型顯示滑移量大于2m的破裂區(qū)域一般不深于20km。值得注意的是，庫塞湖下方30km處的小塊區(qū)域，其破裂幅度達(dá)到3～4m。Mechie等(2004)推斷該區(qū)域地下20km溫度可達(dá)700℃，發(fā)生剛性破裂的可能性較小，而萬永革等(2008)和屠泓為等(2016)單獨(dú)反演GPS/InSAR資料得到的破裂模型并無此特征，因此其結(jié)果可能并不真實(shí)可靠，推測是遠(yuǎn)震波形中包含地表反射波成分導(dǎo)致的假象。

昆侖山口西地震的破裂持續(xù)時(shí)間約130s。前20s發(fā)生雙側(cè)破裂，破裂傳播速度約為 2km/s。在第一段上表現(xiàn)為明顯的走滑特征，在第二段拉張型斷層上，釋放能量較小。20s后破裂向東傳播到主斷層，主破裂區(qū)內(nèi)的平均破裂速度達(dá)到 3.5km/s。主破裂區(qū)上的擴(kuò)展可分為3個(gè)階段：①20～55s，破裂的集中程度和地震釋放能量均較弱，平均破裂速度約 3km/s；②55～80s，矩能量釋放速率達(dá)到峰值，此過程破裂向東傳播約125km，平均速度約 5km/s，大幅度超過局部剪切波傳播速度(3.2～4.0km/s)，表現(xiàn)為明顯的超剪切破裂，該階段釋放了整個(gè)地震約70%的能量，破裂速度與破裂尺度具有一致性； ③80s后，破裂速度和能量釋放速率迅速下降，100s時(shí)矩能量釋放速率達(dá)到一較小峰值，約120s時(shí)，整個(gè)地震破裂過程基本結(jié)束。該地震發(fā)生后6個(gè)月內(nèi)發(fā)生了5次強(qiáng)余震，參考CMT的定位結(jié)果，余震主要集中在庫塞湖東段，與反演模型中的集中破裂區(qū)域位置具有互補(bǔ)關(guān)系，證實(shí)了反演模型的可靠性。

聯(lián)合模型能較好地?cái)M合遠(yuǎn)震波形及近場變形資料。對(duì)于遠(yuǎn)震波形資料而言(圖3)，除個(gè)別位于走向方向上的臺(tái)站對(duì)破裂速度不敏感，整體擬合效果較好，波形擬合的均方根誤差為0.6。對(duì)于GPS數(shù)據(jù)，距離斷層較近的XZ02和BS33兩點(diǎn)擬合效果較好，遠(yuǎn)場擬合度下降，模擬得到的同震位移均小于觀測值。InSAR數(shù)據(jù)整體擬合效果較好(圖4)，LOS向位移平均殘差約0.05m，昆侖山口段殘差稍大(0.10～0.15m)。

注： 黑色實(shí)線為觀測值； 紅色實(shí)線為理論值； 左側(cè)數(shù)字分別表示方位角和震中距； 字母表示體波類型和臺(tái)站名； 右側(cè)數(shù)字表示時(shí)程內(nèi)觀測值峰值圖3 遠(yuǎn)震體波擬合情況

圖4 InSAR數(shù)據(jù)擬合情況

本研究同樣嘗試僅依賴大地測量資料得到反演，得到的破裂滑動(dòng)分布特征與聯(lián)合反演模型類似，但反演得到的地震矩大小為7.9×1020Nm，大于聯(lián)合反演結(jié)果6.1×1020Nm。GPS和InSAR數(shù)據(jù)殘差有一定減小。

反投影結(jié)果如圖5所示。相比于有限斷層反演結(jié)果，反投影從破裂速度變化的角度詮釋了昆侖山口西地震過程的階段性特征。其中，反投影得出最大能量釋放在距震中250km處附近，與有限斷層反演顯示的最大破裂區(qū)完全一致。其區(qū)別僅在超剪切破裂階段，反投影顯示破裂速度可達(dá)到 6km/s，大于有限斷層反演得到的 5km/s的破裂速度，且矩能量釋放更為集中，破裂80s后能量釋放速率迅速下降，但破裂速度未明顯減小。

注： 紅色五角星為震中位置； 灰色圓圈為震源區(qū)域網(wǎng)格； 彩色圓圈為反投影得到的地震能量信息； 圓圈半徑大小代表能量相對(duì)大小； 顏色代表色標(biāo)尺所示的破裂時(shí)間； 左上角為破裂過程在能量-時(shí)間上的投影； 右上角為地震破裂過程在距離-時(shí)間上的投影； 圓圈位置代表能量點(diǎn)在時(shí)間和距離上的分布； 紅色背景線為速度參考線圖5 反投影地震破裂過程

3 討論與結(jié)論

基于大地測量資料的模型反演不能反映震源時(shí)間過程信息(Lasserre et al，2005； 萬永革等，2008)，而完全利用遠(yuǎn)震體波反演模型的空間分辨率較低(Antolik et al，2004)，難以揭示破裂的細(xì)節(jié)特征。大地測量數(shù)據(jù)對(duì)斷層淺層部位(深度小于15km)的破裂敏感，對(duì)深部破裂的分辨率逐步降低，但遠(yuǎn)震波形對(duì)深淺破裂的分辨率基本相同。因此，融合遠(yuǎn)震波形及近場大地測量數(shù)據(jù)建模不僅將震源時(shí)空過程統(tǒng)一表述，同時(shí)也提高了深部破裂分布分辨率。由于近場變形觀測數(shù)據(jù)的約束，聯(lián)合反演模型更少受地球結(jié)構(gòu)復(fù)雜性和地震波傳播的影響，同時(shí)，遠(yuǎn)震波形的約束還可將大地測量觀測中可能存在的震后形變信息分離出來，獲得準(zhǔn)確的地震矩能量。萬永革等(2008)根據(jù)GPS和InSAR數(shù)據(jù)反演得到的總地震矩為9.3×1020Nm，本研究反演得到的地震矩為6.1×1020Nm，更接近于GCMT的結(jié)果，也能解釋部分GPS臺(tái)站及InSAR觀測值大于模擬值。

Wen等(2009)以地質(zhì)調(diào)查結(jié)果作為約束，僅用遠(yuǎn)震體波獲得斷層滑動(dòng)分布，與本研究結(jié)果大體相似，但也有所差別。本文聯(lián)合反演結(jié)果顯示在主破裂段西端地表(91°E～91.5°E)未發(fā)生明顯破裂，而Wen等(2009)強(qiáng)制該部位破裂大小與地表觀測值一致。此外，萬永革等(2008)反演顯示該區(qū)域地表破裂為1～2m，相比于Lasserre等(2005)的結(jié)果(滑動(dòng)分量3～4m)，萬永革等(2008)的結(jié)果由于GPS數(shù)據(jù)的加入，避免了反演結(jié)果的隨意性。Xu等(2006)地質(zhì)調(diào)查結(jié)果表明在該區(qū)域地層具有不確定性，地表觀測到斷層位移可能并非實(shí)際破裂，很大程度上可能是深處破裂產(chǎn)生的強(qiáng)地震運(yùn)動(dòng)導(dǎo)致的地表裂縫，地表測量結(jié)果夸大了實(shí)際破裂值。

有限斷層反演結(jié)果表明，沿?cái)鄬幼呦?55～280km之間的破裂速度為 5.0km/s，且滑移量大于5m的破裂集中在該區(qū)域。280km后破裂速度和能量釋放速率迅速減小，根據(jù)地質(zhì)調(diào)查結(jié)果，斷層在此處分叉且破裂沿南部的昆侖山口斷層繼續(xù)傳播，余震主要集中在該區(qū)域，相較于主震，余震震級(jí)相對(duì)較小，結(jié)合歷史地震分布，昆侖山口東段仍具有較強(qiáng)的地震危險(xiǎn)性。反投影結(jié)果的最大破裂速度約 6km/s，破裂穿過該分叉口能量釋放速率減小，但仍以超過 5km/s的速度繼續(xù)傳播到距震中320km處，與Wang等(2016)利用Hi-net得到的結(jié)果基本一致。

Vallée等(2008)通過分析破裂速度轉(zhuǎn)換帶的高頻地震輻射能量，結(jié)合地質(zhì)調(diào)查資料，表明在320km處，斷層走向存在7.8°的改變，阻礙了超剪切破裂的傳播。同樣，在距震中150km處，斷層走向存在5.7°的改變，使得破裂達(dá)到超剪切速度，表明斷層的幾何結(jié)構(gòu)變化對(duì)超剪切破裂具有較大控制作用。由于反投影方法相對(duì)于有限斷層反演方法使用頻率更高的信號(hào)，對(duì)速度變化具有更高的敏感性，且在有限斷層反演中，破裂速度與滑移量大小具有折中關(guān)系(Fan et al，2014)，故有限斷層反演模型未能識(shí)別出由于走向細(xì)微改變而導(dǎo)致的破裂速度快速變化。