基于QSSP的震源破裂過程反演方法

王博文 張 燕 申重陽 談洪波

1 中國地震局地震大地測量重點實驗室，武漢市洪山側路40號，430071

研究地震震源破裂過程可以認識地震的發(fā)震機理，了解震源性質。震源過程反演是在時間與空間上對地震斷層進行有限離散劃分，在空間上將一個震源斷層分為多個子斷層，在時間上將一個地震事件分為多個子事件，組合全部子斷層塊的震源時間函數(shù)可獲得地震事件的時空滑動破裂分布。可以采用先將問題進行預處理再迭代計算的反演方法，該方法計算速度快，但需要選取合適的初始值[1]。此外還有全局搜索方法，如蒙特卡洛隨機搜索[2]、模擬退火法[3-4]。熱浴算法[5]是模擬退火法的改進，在待解參數(shù)數(shù)量較大時效率更高[6]。目前較常見的地震震源破裂過程研究方法是采用遠震速度數(shù)據(jù)的波形擬合方法或結合位移數(shù)據(jù)采用聯(lián)合反演方法。但這些方法的反演過程始終存在多解性，從同一地震的各個反演結果來看，不確定性較大。

理論上，聯(lián)合反演時使用的數(shù)據(jù)越多，對結果的約束力就越強。重力高頻觀測數(shù)據(jù)能夠很好地反映地震波信息[7]，在一定程度上提高反演結果的可靠性，但目前缺乏利用重力數(shù)據(jù)反演震源破裂過程的相關研究。QSSP軟件[8]是基于Fortran語言開發(fā)的地震學正演軟件，可使用基于球對稱、考慮地球自重力影響的地球模型來計算完整的合成地震圖，并通過將震源分解為多個源的組合模擬計算地表任意多個位置的速度、位移、應力、應變、重力、重力變化等參數(shù)的同震動態(tài)變化。本文基于該功能，將其開發(fā)為反演震源破裂過程的程序，并對2013-04-20蘆山7.0級地震的破裂滑動過程進行反演，以驗證該程序的有效性。

1 反演方法

對于初始值問題，在有初步震源機制解結果的情況下，對斷層模型與地震矩等初始值進行一定限制，根據(jù)矩張量、定位結果可以對破裂起始點、最大滑動量的大致區(qū)域等進行直接約束。對于大地震，單一的震源機制解不能符合實際的震源破裂過程，需要將地震分為多個子事件分別進行反演并疊加。

(1)

(2)

式中，Call為各臺站的理論計算波形，可通過調用QSSP計算獲取，Obsl為各臺站觀測波形，T為最大臺站數(shù)。計算第j個斷層模型參數(shù)的接收概率：

(3)

隨機取ε∈[0,1]，滿足Hj,(k-1)<ε

(4)

式中，β為降溫速度，Tmin為最低溫度。

熱浴算法每次迭代需要計算Q=ΠMn次，雖然相比模擬退火法計算速度更快，但依然需要大量的時間，故僅將其用于全局搜索階段。在該階段，收斂標準可相對放低，使模型參數(shù)能避開局部最優(yōu)解進入全局最優(yōu)解的參數(shù)區(qū)間即可，如圖1紅線中間部分。

圖1 全局最優(yōu)解區(qū)間Fig.1 Interval of global optimal solution

為保證斷層模型滑動的時空連續(xù)性，同時防止計算出現(xiàn)奇異或病態(tài)矩陣，需要先對待解參數(shù)進行光滑約束。待解參數(shù)包括滑動量、滑動破裂方向、滑動起始時間、滑動持續(xù)時間等。計算理論波形數(shù)據(jù)基本原理如下：

G*C=W

(5)

式中，G為格林函數(shù)系數(shù)矩陣，C為震源參數(shù)滑動量(地震矩)、斷層塊走向、傾向、位置與深度、滑動角、起始破裂時間、破裂持續(xù)時間，W為理論計算波形數(shù)據(jù)，該部分包含于QSSP計算模塊中。格林函數(shù)系數(shù)任一元素可由式(6)計算獲得：

φ[θ-(α-1)Δt1-t2]dθ

(6)

式中，guv(T,t0+xΔt;uv,θ)為在uv位置的子斷層塊地震矩引起的臺站波形數(shù)據(jù)格林函數(shù)，t0為波形數(shù)據(jù)起點，Δt為離散時間間隔，α為對應地震矩密度張量時間基本函數(shù)數(shù)量，Δt1為地震矩基本時間函數(shù)間隔，t2為破裂到達uv子斷層塊的時間。

(7)

聯(lián)立得：

(8)

熱浴算法雖然能保證全局性，但其溫度越低，收斂速度越慢，導致在全局最優(yōu)解區(qū)間內(nèi)收斂時間長。故以熱浴算法搜索結果為初始值，在此基礎上對子斷層模型參數(shù)進行區(qū)間限定。采用擬牛頓法進行快速收斂，使波形趨勢擬合程度在得到一定保障的前提下對振幅進行擬合。牛頓法是用泰勒展開保留一階項和二階項來近似非線性函數(shù)F，利用近似函數(shù)的最小值推測非線性函數(shù)的極小值，其展開式為：

F(xi+Δxi)≈F(xi)+J(xi)TΔxi+

(9)

(10)

得增量方程：

Δxi=-H(xi)-1J(xi)

(11)

沿此方向一維搜索最優(yōu)步長β，使F(xi+βΔxi)取最小。

xi+1-xi=H-1[?F(xi+1)-?F(xi)]

(12)

在F為非二次函數(shù)時，仿照式(12)構建滿足式(13)條件的近似矩陣B來替代矩陣H：

xi+1-xi=Bi+1[?F(xi+1)-?F(xi)]

(13)

根據(jù)式(13)即可按式(14)擬牛頓條件計算增量：

Δxi=Bi+1·Δji

(14)

(15)

由謝爾曼-莫里森公式可知，任意非奇異方陣A和n維向量r、s，有條件1+sTA-1r≠0，得：

(16)

采用謝爾曼-莫里森公式迭代求得矩陣B的逆的迭代計算式：

(Bi+1)-1=

(17)

式中，n為x的維數(shù)，將矩陣初始值設置為單位矩陣，再根據(jù)式(17)迭代逼近海塞矩陣。循環(huán)迭代至最小殘差項小于設定精度時停止并輸出對應待解模型參數(shù)。反演過程及擬牛頓算法流程如圖2所示。

圖2 算法流程Fig.2 Flow chart of algorithm

為了檢驗2種算法結合使用的可靠性，采用19×10棋盤模型進行正反演驗證。正演模型參數(shù)隨機給出，設定斷層滑動為45°逆沖，滑動量不超過4 m，斷層位置為100°E、40°N，斷層傾角45°，走向135°，每個子斷層塊大小為3.5 km×3.5 km，斷層起始深度為5 km。主要地震矩釋放一次，集中于5～10 s，破裂總持續(xù)時間15 s，破裂傳播速度3 km/s。反演第1階段中，用熱浴算法進行全局搜索并進行少量迭代收斂，結果如圖3(a)所示。斷層滑動模型及最終對子斷層塊滑動量、滑動角、破裂起始與持續(xù)時間的反演結果分別如圖3(b)、3(c)所示。

對比圖3(a)和3(c)可知，2個理論滑動破裂模型在滑動量、主要滑動區(qū)域及滑動角方面仍存在差距。從圖3(b)和3(c)可以看出，經(jīng)第2階段擬牛頓法在極值附近進行快速收斂后，反演滑動破裂模型與理論值主要破裂區(qū)各子斷層塊間滑動量、滑動角偏差很小，其中滑動量最大偏差不超過18 cm、滑動角偏差不超過5°。由圖3(d)可知，破裂主要集中在4～10 s，持續(xù)時間16 s。綜上，本文方法穩(wěn)定可靠，且擬牛頓法可大幅提升后期收斂速度，較單一熱浴算法搜索和收斂更迅速，因此可用于地震震源破裂過程反演。下文將使用實際震例進行驗證。

圖3 棋盤滑動分布及滑動量時空分布Fig.3 The slip distribution on checkerboard and spatio-temporal distribution of slip

2 斷層模型

本文采用與劉成利等[9]相同的斷層模型，斷層走向214°，傾角38°，每個子斷層大小為3.5 km×3.5 km，斷層面沿走向長66.5 km，沿傾向長35 km，斷層起始深度為6.073 km。滑動量不超過3 m，滑動角0～180°，相鄰斷層破裂起始時間不超過2 s，斷層模型投影如圖4(a)所示(圖中藍色矩形塊為斷層在地表的投影位置，黃色五角星為震中)。可以看出，斷層模型走向與龍門山斷層在該段投影走向一致，且龍門山斷層投影位置基本位于子斷層投影面的中央線處。

圖4 滑動分布投影及活動斷層與臺站分布Fig.4 The projection of the slip distribution and active fault and distribution of the stations

3 數(shù)據(jù)與方法

2013-04-20蘆山7.0級地震前，我國已有18個臺站使用秒采樣PET重力儀。經(jīng)過固體潮校正、零漂校正、跳點校正等處理后有13個臺站數(shù)據(jù)可以使用。首先對加速度數(shù)據(jù)進行積分并去除常數(shù)項，使波形更易于分辨和擬合。根據(jù)范文華等[7]對重力同震數(shù)據(jù)與地震儀速度同震數(shù)據(jù)的對比研究結果，本文對觀測數(shù)據(jù)與格林函數(shù)均進行0.01～0.3 Hz濾波。針對觀測數(shù)據(jù)，將QSSP內(nèi)置濾波器部分源代碼重新編譯形成單獨的濾波軟件進行濾波，選取臺站儀器接收到地震波后70 s的數(shù)據(jù)。臺站分布如圖4(b)所示。

地球模型包括大氣、海洋、地幔、液態(tài)外核、固態(tài)內(nèi)核等多層結構，地殼部分采用Crust2.0模型參數(shù)，其余深部模型使用Kamamori & Anderson地球模型參數(shù)，格林函數(shù)庫由QSSP直接計算得到，反演時從格林函數(shù)庫中自動讀取。

4 蘆山地震震源破裂過程反演

圖5為各臺站波形擬合情況。可以看出，理論計算波形與觀測數(shù)據(jù)波形擬合較好。

圖5 擬合波形Fig.5 Fitted waveform

圖6(b)為本文反演得到的斷層滑動破裂模型，箭頭方向表示滑動方向，箭頭長度表示滑移量，圖中深色子斷層滑動量極小，在實際情況下大概率未發(fā)生任何滑動破裂?？梢钥闯觯瑪鄬踊瑒臃较蛑饕憩F(xiàn)為逆沖，最大滑動量位于第6層，沿傾向距離17.5～21 km、沿走向-1.5～1.5 km，滑移量達157 cm。滑移主要發(fā)生在最大滑移子斷層塊周圍約7～9 km區(qū)域，方向主要為向上發(fā)展。圖6(c)為斷層滑動量時空分布?？梢钥闯?，蘆山7.0級地震的破裂持續(xù)時間為30 s，地震釋放的標量地震矩約1.155×1019Nm，相當于MW6.64地震，起震深度約16.5 km，與最終標定的17.5 km相比略淺。結合圖6(c)和圖7(b)可以看出，破裂滑動集中發(fā)生在10 s、21 s、29 s附近，地震矩釋放主要在前12 s，20～22 s和28～29 s發(fā)生的破裂地震矩釋放和滑動量都相對很小。

本文反演的矩震級為MW6.64，與劉成利等[9]的反演結果MW6.6接近。對于主要破裂滑動區(qū)域，劉成利等[9]反演結果(圖6(a))的最大滑移子斷層塊大約位于15 km，滑動量為141 cm，本文反演結果的最大滑移子斷層塊位于16.8 km，滑動量為157 cm，滑移主要區(qū)域均位于最大滑移子斷層塊周圍7～10 km，本文結果整體相對略深。劉成利等[9]的反演結果中，第1次地震矩釋放峰值中心位于3～7 s，本文反演結果的峰值位于6～11 s，發(fā)生主要破裂滑動的時間相對較晚。

圖6 蘆山MS7.0地震破裂模型Fig.6 The rupture model of Lushan MS7.0 earthquake

圖7 蘆山MS7.0地震地震矩率Fig.7 Seismic moment rate of Lushan MS7.0 earthquake

對比趙翠萍等[10]的反演結果可知，本文反演結果最大滑移區(qū)域中心位置接近，區(qū)域覆蓋相對??；其反演結果最大滑動量為1.8 m，與本文結果相比差異較大；其破裂傳播速度約為3 km/s，與本文結果2.89 km/s較為接近；其矩震級約為MW6.8，比本文的MW6.64大，這可能是主要滑移區(qū)域面積相對較大導致的。趙翠萍等[10]的3次地震矩釋放峰值分別集中于13 s、27 s與39 s附近，本文發(fā)生主要破裂滑動的時間與之相比均較早，破裂的整體持續(xù)時間也相對較短。鄭緒君等[11]的反演結果中，地震矩只釋放1次，整個破裂過程約為10 s，釋放峰值位于5～7 s附近，本文反演的主要地震矩釋放時間與之較為接近；其震中深度為16 km，主要滑動區(qū)域位于震中周圍約6 km，破裂滑動范圍更為集中，本文結果與之相比起震深度較為接近但范圍更大；其斷層面最大滑動量約為2.2 m，該結果相對其他結果及本文結果均較大。

各反演結果間均存在一定差異。本研究與劉成利等[9]的研究雖采用相同模型，但所得結果存在些許不同，這可能與使用QSSP及采用數(shù)據(jù)類型不同有關。重力數(shù)據(jù)與測震儀數(shù)據(jù)在頻率域的P波頻譜具有一定的相似度但仍存在差異，QSSP以分層地球模型計算理論波形且考慮海洋的影響，與半無限空間地球模型計算不同，這些因素都會使計算結果出現(xiàn)變化。

5 結 語

基于QSSP軟件，本文利用熱浴算法、擬牛頓算法和平滑約束反演震源破裂過程。從棋盤模型驗證及實際震例反演結果來看，本文反演方法收斂效率高、穩(wěn)定可靠。