夏 晨 趙伯明
(北京交通大學(xué)土木建筑工程學(xué)院,北京 100044)
近年來全球范圍大地震頻發(fā),比如1989年美國Lama Prieta地震、1995年日本Kobe地震、1999年臺灣集集地震、1999年土耳其Kocaeli地震、2005年印尼蘇門答臘島地震、2008年中國汶川大地震以及2011年東日本大地震。這些地震給人類帶來巨大災(zāi)難的同時,也極大豐富了近場強(qiáng)震記錄。地震學(xué)及地震工程學(xué)界普遍認(rèn)識到,近場強(qiáng)地面預(yù)測對減輕未來大地震造成的災(zāi)害和指導(dǎo)工程建設(shè)的抗震設(shè)計具有重要作用。如何對未來大地震的強(qiáng)地面運(yùn)動進(jìn)行可靠預(yù)測,是地震學(xué)和地震工程所面臨的重大挑戰(zhàn)之一,也是地震學(xué)發(fā)展的成果應(yīng)用到地震工程的一個橋梁。
經(jīng)驗格林函數(shù)法是地震動數(shù)值模擬計算中常用的方法之一,其基本思想是:大震產(chǎn)生的地面運(yùn)動是由大震斷層中一系列子源所產(chǎn)生的地面運(yùn)動疊加而成,子源的地震動由大震附近的余震記錄表示。該方法由Hartzell(1978)首次提出,其優(yōu)點(diǎn)在于余震記錄包含震源、傳播路徑和場地效應(yīng)(地震動模擬中最重要的3個要素)的影響,避免理論數(shù)值計算的復(fù)雜過程以及對地下傳播介質(zhì)信息的依賴。運(yùn)用經(jīng)驗格林函數(shù)法合成強(qiáng)地震動,關(guān)鍵在于小震記錄的選取,作為格林函數(shù)的小震應(yīng)與大震震源位置相近,震源機(jī)制解相似。從理論上考慮,小震的震級越小越接近點(diǎn)源,疊加合成時越符合大震的震源時間函數(shù)。但在實際應(yīng)用中,小震級的地震記錄是有限的,另外,震級太小時,記錄的信噪比不高,將會影響結(jié)果的可靠性。
Hartzell提出經(jīng)驗格林函數(shù)法之后,得到廣泛關(guān)注,多位學(xué)者對此方法進(jìn)行改進(jìn)并用于近場強(qiáng)地震動模擬(Kanamori,1979;Tanaka,1982;Irikura,1983,1986;Joyner et al.,1986;Heaton et al.,1989;Hutchings,1991;Kamae et al.,1992;Frankel,1995;Hough,2001;羅齊峰等,1994;金星等,2002)。在疊加小震記錄時,對于所疊加的小震記錄數(shù)量選取,有2種常用的準(zhǔn)則。第1種是所疊加的小震數(shù)為N(假設(shè)大小震地震矩之比為N);第2種是所疊加的小震數(shù)為N2/3,即等于大小震的斷層面面積之比。為后文表述方便,記第1種為A方法,第2種為B方法。在實際應(yīng)用中,這2種方法均有缺陷。A方法合成的地震動在低頻部分符合觀測記錄,但是高頻部分偏高;B方法合成的地震動在高頻部分符合觀測記錄,但是低頻部分偏低。對這些缺陷,Joyner等(1986)提出隨機(jī)疊加N4/3個小震記錄,對疊加后的結(jié)果乘以系數(shù)N-1/3,以便得到高頻和低頻均較為理想的結(jié)果。一些學(xué)者提出將小震記錄與隨機(jī)震源時間函數(shù)求卷積,再疊加N個小震記錄,不改變低頻結(jié)果,降低高頻結(jié)果。另一些學(xué)者提出疊加N2/3個記錄,再采用相應(yīng)的過濾器在不改變高頻結(jié)果的情況下提高低頻結(jié)果(Frankel,1995)。
在傳統(tǒng)經(jīng)驗格林函數(shù)方法中,未考慮大小地震拐角頻率的差異,由經(jīng)驗格林函數(shù)法合成的大震記錄,其拐角頻率等于所用的小震拐角頻率,而實際中拐角頻率隨地震震級的增大而減小。這是引起經(jīng)驗格林函數(shù)法合成結(jié)果不能在高頻和低頻同時符合觀測記錄的一個重要因素。考慮大小震拐角頻率的不同,基于ω2震源譜模型,通過理論推導(dǎo)構(gòu)建一個新的函數(shù),將動力學(xué)拐角頻率引入到經(jīng)驗格林函數(shù)法(B方法)中,改進(jìn)其合成的結(jié)果,并運(yùn)用改進(jìn)后的方法合成汶川大地震近場強(qiáng)地震動,通過與觀測記錄對比,驗證動力學(xué)拐角頻率在經(jīng)驗格林函數(shù)法中的有效性。
以Tanaka(1982)提出的一個經(jīng)驗格林函數(shù)法為B方法的代表。Tanaka(1982)基于大震與小震的定標(biāo)關(guān)系,導(dǎo)出如下關(guān)系(Kanamori et al.,1975):
式(1)中:L,W,D,τ,M0分別為大震的震源長度、寬度、斷層面平均位移、上升時間和地震矩;Le,We,De,τe,M0e分別為小震的相應(yīng)參數(shù)。大地震的斷層面由n2個大小相同的子源構(gòu)成,每個子源與小震大小相同。將小震記錄ae(t)作為子源格林函數(shù),疊加合成大地震的地震動記錄a(t),如下所示:
式(2)、(3)中,c和v分別為地震波傳播速度和破裂傳播速度;Re表示小震震源到觀測點(diǎn)的距離;R0表示大震斷層面破裂點(diǎn)到觀測點(diǎn)的距離;Rij表示第ij個子源到觀測點(diǎn)的距離;rij表示第ij個子源到破裂起始點(diǎn)的距離。
在實際計算中,式(2)的疊加過程是在頻域中完成的,其表達(dá)式如下:
式(4)中:Ae(f)為小震記錄ae(t)的傅立葉變換。將疊加之和A(f)通過傅立葉逆變換轉(zhuǎn)換到時域,即可得到大地震的地震動記錄a(t)。
考慮到大小震應(yīng)力降Δσ的差異,在文中采用Tanaka(1982)經(jīng)驗格林函數(shù)法時,式(4)改寫為如下形式:
傳統(tǒng)的經(jīng)驗格林函數(shù)方法未考慮大小地震拐角頻率的差異,合成的大地震具有與小地震相同的拐角頻率,這不符合拐角頻率隨震級增大而減小的實際規(guī)律。這也是導(dǎo)致經(jīng)驗格林函數(shù)法(B方法)模擬結(jié)果與觀測記錄在低頻出現(xiàn)差異的主要因素之一。
Motazedian等(2005)提出動力學(xué)拐角頻率,將第ij個子源的拐角頻率表示為
式(6)中,N為主斷層面上的子源總數(shù),NR為破裂傳播到第ij個子源時累計滑動的子源個數(shù),β為介質(zhì)剪切波速,Δσ為主斷層應(yīng)力降,M0為主斷層地震矩。用動力學(xué)拐角頻率f0ij替代Boore(1983)隨機(jī)點(diǎn)源模型方法中的拐角頻率f0,再將隨機(jī)點(diǎn)源模型方法計算的點(diǎn)源地震動作為子源的格林函數(shù),運(yùn)用隨機(jī)有限斷層模型程序EXSIM疊加合成主斷層地震動。由于拐角頻率在每個子源的計算中是顯式的,因此動力學(xué)拐角頻率可方便地引入到程序EXSIM中進(jìn)行應(yīng)用。
孫曉丹等(2009)提出一種與震源譜模型相關(guān)的動力學(xué)拐角頻率,改善原來的動力學(xué)拐角頻率定義中對破裂面上地震波頻率輻射的不均勻性表達(dá),考慮到拐角頻率隨破裂面積增加而有所下降的趨勢。
但是,在經(jīng)驗格林函數(shù)法中,每個子源的信息包含在小震記錄里,拐角頻率并不體現(xiàn)在計算過程中。若將動力學(xué)拐角頻率引入到經(jīng)驗格林函數(shù)法中,需要做一些處理。
假設(shè)大小地震的震源譜均符合Brune(1970)的ω2模型:
圖1 經(jīng)驗格林函數(shù)法采用的主震斷層模型(Tanaka,1982)Fig.1 The main fault plane model described by Tanaka(1982).
對于小地震,其拐角頻率f0small通過地震動記錄反演得到;對于每個子源,假設(shè)f0ij與累計滑動的子源個數(shù)NR的平方根成反比,表示如下:
式(8)中,f0main為主斷層地震的拐角頻率,其余參數(shù)與上文一致。在Boore(1983)隨機(jī)點(diǎn)源模型方法中,拐角頻率由式(9)計算獲得:
由式(9),可得
假設(shè)把已經(jīng)滑動的NR個子源組成一個斷層面,這NR個子源的總地震矩為M0ruptured,基于大震與小震斷層參數(shù)的定標(biāo)關(guān)系(Kanamori et al.,1975),可以得到如下關(guān)系式:
由式(8),(11)和(12),可以得到如下關(guān)系:
由式(13)可以看出,動力學(xué)拐角頻率f0ij與已發(fā)生破裂子源的總地震矩M0ruptured立方根成反比。將NR個子源看作一個斷層面,式(13)具有與式(9)同樣的意義,即
基于ω2模型,為引入動力學(xué)拐角頻率f0ij到經(jīng)驗格林函數(shù)法中,構(gòu)建一個函數(shù)如下:
為保證Hij(f)在高頻(高于小震拐角頻率的頻段)時接近1,即函數(shù)Hij(f)不影響經(jīng)驗格林函數(shù)法高頻計算的結(jié)果,參數(shù)Cij必須滿足以下條件:
由式(8),(10)和(15),可得
進(jìn)而,函數(shù)Hij(f)可表達(dá)為
函數(shù)Hij(f)與Frankel(1995)提出的相對滑移速度函數(shù)具有相似的形式,這是因為這2個函數(shù)都是基于ω2模型構(gòu)建的,但是,二者在計算參數(shù)時基于不同的假設(shè)。Frankel(1995)基于所有疊加的子源地震矩之和與主斷層地震矩相等這一假設(shè),計算參數(shù)C和f0ij,得出f0ij與成正比。本文對于每個子源,假設(shè)f0ij與累計滑動的子源個數(shù)NR的平方根成反比,且Hij(f)在高頻(高于小震拐角頻率的頻段)時接近1,即函數(shù)Hij(f)不影響經(jīng)驗格林函數(shù)法高頻計算的結(jié)果,由此得出動力學(xué)拐角頻率f0ij與已發(fā)生破裂子源的總地震矩M0ruptured立方根成反比,即f0ij∝M0ruptured-1/3。可以看出,在Frankel(1995)的研究中,f0ij是不變量,本文中f0ij隨破裂傳播而變化。
函數(shù)Hij(f)可以方便地與經(jīng)驗格林函數(shù)法結(jié)合,以引入動力學(xué)拐角頻率f0ij。引入函數(shù)Hij(f)后的式(5)改寫為
改進(jìn)后的經(jīng)驗格林函數(shù)法可方便地實現(xiàn),其過程如圖2所示。
為檢驗動力學(xué)拐角頻率引入經(jīng)驗格林函數(shù)法的有效性,分別運(yùn)用改進(jìn)前后的方法對汶川大地震主震斷層周圍6個臺站點(diǎn)的地震動進(jìn)行模擬計算。通過對比改進(jìn)前后合成結(jié)果與觀測記錄的加速度反應(yīng)譜,驗證動力學(xué)拐角頻率對合成結(jié)果低頻成分的提高效果。
2008年5月12日14點(diǎn)28分(北京時間),在青藏高原東緣的龍門山斷裂帶上發(fā)生特大地震,簡稱5·12汶川地震,其震級為MW7.9;起始破裂點(diǎn)(震中)位于30.986°N,103.364°E;震源深度19km;震源破裂長度達(dá)300km,寬40km。
多位學(xué)者通過遠(yuǎn)場體波波形記錄結(jié)合近場強(qiáng)地震動數(shù)據(jù),以及GPS和SAR等數(shù)據(jù),利用反演技術(shù)重建汶川地震的震源破裂過程(Ji et al.,2008;Shen et al.,2008;王衛(wèi)民等,2008;Yarai et al.,2008;Koketsu,2009)。Koketsu(2009)的反演結(jié)果顯示,震源斷層面由西南段和東北段兩部分構(gòu)成,如圖3所示。整個斷層面長300km,寬40km。兩部分的走向均接近230°,西南段的傾角為35°,東北段的傾角65°。
Kurahashi等(2010)根據(jù)近斷層6個臺站的強(qiáng)震記錄,提出震源斷層存在4個凹凸體(asperity),其中3個在西南段,1個在東北段上,地震動主要由西南段上的3個凹凸體引起。本文采用Koketsu(2009)的反演結(jié)果作為震源模型,并采用Kurahashi等(2010)建議的凹凸體參數(shù)(表1),其分布如圖4所示。選取3個余震的記錄分別作為3個凹凸體區(qū)域的格林函數(shù),余震參數(shù)見表2,其分布如圖4所示。圖4中,3個黑色圓點(diǎn)為所取3個余震的震中位置;HSL為黑水雙溜索臺站,JYC為江油重華臺站,LDD為瀘定得妥臺站,LXT為理縣桃坪臺站,MXD為茂縣疊溪臺站,SPA為松潘安宏臺站。
選取主震斷層周圍6個臺站(HSL,JYC,LDD,LXT,MXD,SPA),空間分布如圖4所示。參照上文中的震源模型,分別采用未修改的經(jīng)驗格林函數(shù)法和引入動力學(xué)拐角頻率后的經(jīng)驗格林函數(shù)法合成汶川地震近場強(qiáng)震記錄。合成結(jié)果和觀測記錄都采用四階Butterworth濾波器進(jìn)行0.5~20Hz的帶通濾波。改進(jìn)前后的合成結(jié)果與觀測記錄的加速度反應(yīng)譜如圖5所示,加速度時程如圖6所示。
圖2 改進(jìn)后的經(jīng)驗格林函數(shù)法實施過程示意圖Fig.2 The schematic diagram of improved empirical Green's function method.
圖3 汶川地震斷層面上的破裂滑動分布(引自Koketsu et al.,2009)Fig.3 Inverted slip distribution on the fault plane of Wenchuan earthquake(after Koketsu et al.,2009).
從圖5所示的加速度反應(yīng)譜對比可以看出,未修改的經(jīng)驗格林函數(shù)法合成結(jié)果(黑色虛線)在低頻部分明顯低于觀測記錄的結(jié)果(黑色粗實線),引入動力學(xué)拐角頻率后的合成結(jié)果(灰色細(xì)實線)提高了低頻部分結(jié)果,高頻部分結(jié)果基本不變化,低頻部分提高后與觀測記錄接近。尤其在HSL、LDD和SPA 3個臺站,低頻部分的結(jié)果與觀測記錄符合度較高。部分臺站模擬的結(jié)果與觀測記錄有差異,比如JYC的合成結(jié)果,引起差異最重要的原因是破裂傳播效應(yīng),而經(jīng)驗格林函數(shù)法不能包含這部分因素對合成結(jié)果的影響。本文重點(diǎn)在于驗證引入動力學(xué)拐角頻率后,經(jīng)驗格林函數(shù)法合成的低頻地震動結(jié)果是否得到了有效提高。從這方面看,圖5中6個臺站的結(jié)果顯示,改進(jìn)后的方法對于改善經(jīng)驗格林函數(shù)法的低頻結(jié)果是有效的。
表1 凹凸體參數(shù)(引自Kurahashi et al.,2010)Table 1 Parameters of the asperities(Kurahashi et al.,2010)
圖4 震源斷層凹凸體,余震震中,以及若干臺站分布圖(修改自Kurahashi et al.,2010)Fig.4 The locations of asperities,epicenters of aftershocks and some of the stations(adapted from Kurahashi et al.,2010).
表2 選取的余震相關(guān)參數(shù)Table 2 List of parameters of selected aftershocks
圖5 觀測記錄與合成結(jié)果的加速度反應(yīng)譜對比圖Fig.5 Comparison of acceleration spectra of observed and synthesized ground motions.
圖6 觀測記錄與合成結(jié)果的加速度波形對比圖(1)Fig.6 Comparison of observed and synthesized acceleration waveforms(1).
圖6 觀測記錄與合成結(jié)果的加速度波形對比圖(2)Fig.6 Comparison of observed and synthesized acceleration waveforms(2).
圖6對比了引入動力學(xué)拐角頻率之后合成的加速度波形圖與實際記錄的加速度波形圖。大部分臺站合成的加速度波形接近觀測的波形,加速度峰值比在0.5~2。LDD和LXT臺站的水平分量結(jié)果偏大。JYC臺站的結(jié)果在65~80s與實際記錄有明顯差異,這是由于在本文的模擬中,只取西南段的發(fā)震斷層作為震源,因此JYC的結(jié)果在70~80s缺少西北段斷層引起的地震動。但是西南段的結(jié)果與實際記錄較為符合。
本文從大小震拐角頻率差異為出發(fā)點(diǎn),假設(shè)主斷層上每個子源的拐角頻率與累計滑動的子源個數(shù)平方根成反比,基于ω2震源譜模型,通過理論推導(dǎo)構(gòu)建一個新的函數(shù),將動力學(xué)拐角頻率引入傳統(tǒng)經(jīng)驗格林函數(shù)法中。以Koketsu(2009)的震源過程反演結(jié)果作為震源模型,并參考Kurahashi(2010)建議的凹凸體參數(shù),分別用改進(jìn)前和引入動力學(xué)拐角頻率后的經(jīng)驗格林函數(shù)法合成5·12汶川地震近場強(qiáng)地面運(yùn)動加速度時程。合成結(jié)果分別與觀測記錄的加速度反應(yīng)譜,以及加速度波形圖進(jìn)行了對比。
(1)通過比較合成和觀測的加速度反應(yīng)譜,未修改經(jīng)驗格林函數(shù)法的大部分結(jié)果在低頻相對觀測結(jié)果偏低;而引入動力學(xué)拐角頻率的合成結(jié)果在低頻部分明顯提高,并與觀測結(jié)果較為符合,其高頻部分基本保持不變。動力學(xué)拐角頻率引入經(jīng)驗格林函數(shù)法后有效地改進(jìn)了合成的結(jié)果。
(2)汶川大地震地震動合成計算顯示,從加速度反應(yīng)譜以及加速度波形圖可看出,大部分合成結(jié)果與觀測記錄是符合的,說明經(jīng)驗格林函數(shù)法對于MW7.9地震的地震動計算同樣有效。經(jīng)驗格林函數(shù)法可對未來可能發(fā)生的大地震地震動進(jìn)行預(yù)測。