基于隨機有限斷層法的地震烈度計算研究*

王 俊，霍祝青，詹小艷，宋 浩，江昊琳，徐 戈

(江蘇省地震局，江蘇南京210014)

基于隨機有限斷層法的地震烈度計算研究*

王 俊，霍祝青，詹小艷，宋 浩，江昊琳，徐 戈

(江蘇省地震局，江蘇南京210014)

采用基于動力學拐角頻率的隨機有限斷層法，首先模擬計算了東臺、如皋及大豐3個臺站記錄的2006年江蘇東臺3.8級地震的地震動，并與其實際記錄進行比較研究，驗證了該方法的計算結果能夠較好地反映出江蘇地區主震剪切波的主要特征。隨后運用該方法模擬計算了1979年江蘇溧陽6.0級地震震源區內156個虛擬觀測點的地震動，并通過內插值的方式得到此次地震的加速度場分布情況;隨后依據峰值加速度統計方法和模糊評定方法，分別計算此次地震的理論烈度分布，結果顯示:兩種烈度劃分方法的結果與實際調查的烈度分布都基本一致，但模糊評定方法的劃分結果更加接近真實的烈度分布。

動力學拐角頻率;隨機有限斷層法;合成地震動;地震烈度;溧陽6.0級地震

0 引言

地震烈度是用來描述地震對地面的影響和破壞程度的參數。破壞性地震發生后，除震源的基本參數外，政府和社會公眾最迫切想了解的是極震區的震害分布情況，從而為緊急救援爭取時間。目前在我國獲取震害分布情況的主要途徑仍是震后實際的震害考察和評估 (崔建文等，2008)，最短時間也需要數天，不利于救援工作的開展。據統計，在一次地震中，90%以上的建 (構)筑物的破壞是由地震動造成的 (王海云，謝禮立，2009)。地震動是由震源破裂、波在地殼介質中的傳播以及場地反應等共同形成的復雜地面運動，在低頻段具有確定性而在高頻段則表現出強烈的隨機性的特征 (Somerville et al，1991;Beresnev，Atkinson，1997)。

對合成地震動的研究，特別是對近場強地震動的研究一直是地震學和地震工程學領域的熱點問題。目前國內外已經發展了多種模擬方法，如Aki和Richaids(1986)提出的將地震動表示為震源時間函數和格林函數卷積的方法，可用于長周期地震動的合成;高斯帶限白噪聲的隨機點源方法 (Hanks，McGuire，1981;Boore，1983) 和經驗格林函數方法 (Hartzell，1978;Irikura，1983)，可用于高頻地震動的合成。其中經驗格林函數方法是最為有效的，主要思路是將研究區域內與主震震源機制相似的、信噪比較高的小震記錄作為子源響應，從而建立起合成地震動的經驗格林函數。該方法的優點是在地震動的模擬過程中不需要再計算地震波的傳播路徑效應和場地響應，因為小震記錄中已包含了這些信息，但不足之處是無法在地震記錄稀疏或缺乏余震記錄的地區使用。

隨機有限斷層震源模型方法解決了經驗格林函數方法中對余震記錄要求較高這一問題 (Beresnev，Atkinson，1997，1998)。對于遠源和小震，采用隨機點源模型來模擬;對于近源或強震，將主震的破裂過程視為破裂面上一系列子源破裂的結果 (Beresnev，Atkinson，2002)，并用隨機方法產生滿足Brune震源模型的子源 (Brune，1970)，最后通過疊加各子源的地震動合成具有主震破裂特征的地震動。該方法已被應用于國內外許多地區的地震動參數估計中，在多次大震地震動的模擬中顯示出良好的效果 (Atkinson，Beresenev，2002;Zafeiria，Beresnev，2003;石玉成等，2005;崔建文等，2008)。但該方法基于靜力學拐角頻率，合成地震動會隨著子斷層的尺寸變化發生顯著的變化，Motazedian和Atkinson(2005)引入了動力學拐角頻率將其進一步改進，從而彌補了原方法中的缺陷。

場點的地震動與地震烈度之間是怎樣的一種關系?如何依據地震動快速地評定出震區的烈度分布，是震后緊急救援更加關注的實際問題。目前，地震烈度的計算方法主要歸納為統計回歸法和模糊評定法 (袁一凡，1998;李山有等，2002)。統計回歸法是通過統計回歸的方法，僅考慮10 Hz左右的地震動峰值加速度與地震烈度之間的函數關系。袁一凡 (1998)考慮到宏觀烈度定義的模糊性，引入模糊評定方法來建立地震動參數與地震烈度之間的關系，該方法可以將更多的地震動參數納入評定系統中，能夠更好地反映出兩者的內在關系。

目前，江蘇省雖然已形成了布局較為合理的微震、強震動觀測網絡，但仍然沒有正式開展地震烈度速報的相關工作。1979年江蘇溧陽發生M6.0地震，造成了巨大的經濟損失和人員傷亡，震后專家學者進行了大量的科學考察和震害評估，繪制了震區真實烈度分布圖 (胡連英等，1997)。為了在將來破壞性地震發生后，實現快速劃分地震烈度，筆者以溧陽地震為例，首先采用基于動力學拐角頻率的隨機有限斷層法，合成研究區內一定數量的虛擬場點的地震動;然后根據峰值加速度統計法和模糊評定方法計算虛擬場點的理論烈度分布，并與真實烈度分布進行對比研究，以判斷基于合成地震動的地震烈度計算的有效性。

1 基于動力學拐角頻率的隨機有限斷層法

Beresnev和Atkinson(1997，1998)將地震斷層面視為由若干個大小相等的矩形子斷層所構成的，每個子斷層即為一個點源 (或子源)，如圖1所示。地震的破裂過程是從破裂起始點以一定的破裂速度 (一般為剪切波速的0.8倍)向外呈輻射狀傳播，當傳播到每個子源的中心時，該子源即被觸發，且子源的震源譜符合 Brune模型(Brune，1970)。因此，某一時刻所有子源在某一基巖觀測點上由剪切波引起的地震動，可在時域中經過相位校正后進行疊加而求得，表示為

式中，NL、NW分別為沿著斷層走向和下傾方向的子斷層數;NS為子源被觸發的次數;Δtijk為地震波從破裂起始點傳播至第ij個子源的時間延遲與該子源至觀測點的傳播時間延遲之和;aij(t)為第ij個子源在觀測點由剪切波引起的地震動。

圖1 隨機有限斷層法的幾何示意圖Fig.1 Geometric schematic diagram of stochastic finite faults method

Beresnev和Atkinson(1997，1998)在合成地震動時，使用的是基于Brune模型 (Brune，1970)的靜力學拐角頻率，其對子斷層的大小有嚴格的限制，為了保證地震矩的能量守恒，要求一個子源被多次觸發。在合成大地震的地震動時要求子斷層的尺寸在5～15 km范圍內，限制了對小震地震動的模擬，并且難以模擬出大震的長周期成分(Motazedian，Atkinson，2005)。為了彌補使用靜力學拐角頻率合成方法的缺陷，Motazedian和Atkinson(2005)引入了動力學拐角頻率，其優點是通過給各個子源賦予不同的拐角頻率，來表示破裂面上地震波頻率的非均勻性，同時達到拐角頻率隨破裂面積的增大而下降的趨勢。將第ij個子源的震源加速度譜定義改寫為

式中，C是標量常數 (Boore，1983);Hij=1/{exp(－πfk)exp(－πfRij/Qβ)/Rij}，是保證子源高頻輻射守恒的標度因子;其中，Q是介質的彈性衰減系數;β是剪切波速;Rij是體波的幾何衰減模型;foij是子源的拐角頻率，可以表示為在t時刻時子斷層破裂總數NR(t)的函數:

式中，E為整個斷層在高頻段的總輻射能;Moave=M0/N，為子斷層的平均地震矩。

2 地震動的合成

為了檢驗基于動力學拐角頻率的隨機有限斷層法，筆者首先以2006年12月26日江蘇東臺MW3.8地震為實例進行模擬。該次地震后分布在不同方位的東臺、如皋及大豐3個強震臺清晰地記錄到此次地震，震中距分別為22.4 km、36.4 km、53.3 km，其震中位置及震源機制分布如圖2所示。

圖2 溧陽M6.0地震、東臺MW3.8地震的震中位置及震源機制Fig.2 Epicenter locations and focal mechanisms of Liyang M6.0 and Dongtai MW3.8 earthquakes

孫業君 (2010)①孫業君.2010.2006年東臺地震M3.8地震破裂參數研究.利用P波初動反演得到了江蘇東臺MW3.8地震的震源機制解，并且根據地震多普勒效應計算了震源的破裂參數 (表1)。隨機有限斷層法的全局基本參數參照Motazedian和Atkinson(2005)的研究結果 (表2)。地殼初始模型參考黃耘等 (2006)根據人工地震探測的研究結果，在10～15 km深度范圍內江蘇地區的剪切波速約3.5 km/s;(2)式中標度因子里的彈性衰減系數Q值采用詹小艷等 (2011)利用S波反演得到的結果Q(f)=271.8×f0.5671。

表1 兩次地震的震源破裂參數Tab.1 Source rupture parameters of two earthquakes

表2 隨機有限斷層法的基本參數Tab.2 Basic parameters of stochastic finite faults method

取蘇中沿海1～10 Hz范圍內的平均場地響應值2.739(詹小艷等，2011)，對合成地震動進行場地效應校正，最終得到東臺、如皋及大豐3個臺站的合成地震動峰值加速度分別為2.44、1.23、1.37 cm/s2，實際記錄的地震動峰值則分別為2.62、1.28、1.85 cm/s2，兩者十分接近;在剪切波的卓越持續時間、波形包絡線的近似程度方面，東臺、如皋2個臺的結果與真實記錄也十分相似(圖3)，大豐臺的相似度則較低，分析認為這可能與其距此次震中的距離較遠和臺基地下結構差異有關;3個臺的合成地震動在5～15 Hz頻帶范圍內、阻尼比為2%的反應譜計算結果與實際記錄結果均基本一致，并且隨著震中距的增加，反應譜在低頻段的一致性更加明顯，說明基于動力學拐角頻率的合成地震方法具有更寬頻帶的優點，能夠較好地反映出主震剪切波的主要特征，可適用于江蘇地區。

1979年7月9日江蘇溧陽發生M6.0地震，極震區 (Ⅷ度)實際的長軸為北西西向，而重度破壞區 (Ⅶ度)的長軸卻是北北東向，兩者的明顯差異引起一些學者對地震斷層判斷的爭議 (孫受成等，1980)，震后眾多學者對其震源破裂參數做了詳細的研究。筆者對溧陽地震合成地震動計算所需的震源破裂參數和斷層滑動模型，采用被普遍接受的研究結果和方法來進行構建:胡連英等 (1997)對其震源機制解的研究認為，主震的斷層面沿矛東斷裂帶發生右旋走滑錯動、兼有正斷層性質的傾滑錯動;林邦慧等 (1982)利用P波頻譜分析方法推算了地震斷層的破裂長度約為19 km;王煒 (1984)由破裂長度和破裂速度進一步地推算了此次地震的地震矩，平均位錯為0.53 m。子斷層的尺度依據 Beresnev和 Atkinson(2002)提出的下式進行劃分:

可得子斷層的尺度為2×2 km，并結合McGinty和Robinson(1999)的方法建立起主震斷層面的滑動模型，如圖4所示。

圖3 東臺(a)、如皋(b)及大豐臺(c)的合成地震動和實際記錄波形及阻尼比為2%的反應譜Fig.3 Synthetic ground motions and the actual waveforms recorded by Dongtai(a)，Rugao(b)，Dafeng(c)Stations，and their response spectrum of damping ratio is 2%

圖4 溧陽M6.0地震的子斷層及滑動模型Fig.4 Sub-faults and the sliding model of Liyang M6.0 earthquake

為了得到震源區 (31.1°～31.9°N，119.0°～119.7°E)基于合成地震動的烈度分布，如圖2中矩形方框區域所示，需要獲取震源區內一定數量的虛擬觀測點的地震動記錄。于是我們采用上述方法和震源參數來模擬得到溧陽地震震源區內156個節點處 (以0.06°為間隔)的地震動數值，并通過樣條曲面插值方法進行內插值，最終得到其加速度場的分布。圖5b顯示峰值加速度大于180 cm/s2的區域長軸走向為近北北西向;而大于125 cm/s2的區域長軸走向為近北東東向，兩者差異明顯，這可能反映出極震區的合成地震動分布與震源的破裂模式有關。

3 地震烈度的劃分

溧陽M6.0地震后，一些學者進行了大量的震害考察和科學評估，以《新的中國烈度表(1980)》為烈度線劃分的依據得到了真實烈度分布 (胡連英等，1997):上沛西塘至慶豐為極震區，烈度達Ⅷ度，等震線長軸約7 km，呈北西西向，短軸為4.5 km，呈北北東向，總面積約為20 km2;Ⅶ度區主要是河口—舊縣—東岳廟區域的一個規則橢圓，長軸走向為北北東向，面積約為389 km2;Ⅵ度區主要是社渚—金壇的似橢圓區域內，面積約為2 194 km2;Ⅴ度區的總面積約為10 945 km2，長軸為北北西向，如圖5a所示。

為了檢驗基于合成地震動的地震烈度計算是否可行，是否可以運用在江蘇地區今后的破壞性地震中，本文將采用2種劃分方法來進行對比研究。一種為參考《新的中國烈度表 (1980)》中的物理指標來劃分，即根據水平向的峰值加速度值，Ⅴ度時平均為31 cm/s2、Ⅵ度時平均為63 cm/s2、Ⅶ度時平均為125 cm/s2和Ⅷ度時平均為250 cm/s2，其結果如圖5b中黑色線條所示。另一種則是采用袁一凡 (1998)提出的模糊評定方法，共分為8檔，即:小于Ⅳ、Ⅳ、Ⅴ、Ⅵ、Ⅶ、Ⅷ、Ⅸ和大于Ⅸ。參與烈度計算的地震動參數主要為:(1)峰值加速度A:垂直峰值加速度和水平峰值加速度的比值作為一個獨立的因子，當兩者之比小于0.3時，烈度肯定小于Ⅸ度，因此它可以作為否定因子;(2)地震動的卓越頻率; (3)20%的相對持時，即最初和最后到達峰值20%的兩點間時間;(4)1、2、5、8 Hz頻率點的反應譜值，其中5 Hz和1 Hz為主要值。但該方法的步驟較為復雜，且本文中經過內插值后虛擬場點的數量較多。因此，我們只計算了破壞最嚴重的Ⅵ和Ⅶ度區的分布，結果如圖5b中紅色線條所示。

由圖5可以看出，基于合成地震動的兩種烈度劃分在相同的烈度區內基本重合 (圖5b)，且與真實的烈度分布 (圖5a)基本相符。通過比較可以看出采用模糊評定方法劃分的Ⅵ度區，其長軸走向為北北東向、與真實烈度Ⅵ度區的長軸走向更為接近;Ⅶ度區的方向性不明顯，真實烈度的Ⅶ度區長軸走向和極震區 (Ⅷ度區)的長軸走向差異明顯，難以準確綜合地判斷其方向性;而基于峰值加速度評定的Ⅶ度區則是近“十”字形、長軸近北東東向;這表明模糊評定方法能更為真實地反映出各地震動參數與真實烈度之間的內在關系，更接近真實烈度的分布，體現了能綜合不同因子來進行判定的優勢。

圖5 溧陽M6.0地震的實際烈度 (a)(據胡連英等，1997)和合成地震動的理論烈度 (b)Fig.5 Actual intensity(a)(according to Hu et al，1997)and the theory intensity of synthetic ground motion based on two methods(b)of Liyang M6.0 earthquake

4 結論與討論

從圖5的計算結果來看，采用基于動力學拐角頻率的合成地震動來劃分烈度的方法，能夠較為合理地反映出震源區的真實烈度分布。相比于地震烈度速報臺網，通過合成地震動來實現震源區的烈度劃分可能需要花費更長的時間，但仍可在震后獲取到可靠的震源破裂參數后在數小時內完成，并可根據震源破裂參數不斷地對結果進行修正，這比通過震后的震害考察評估來獲得烈度分布速度快;目前我國大部分地區并不具備烈度速報的條件。因此，本文研究的快速理論烈度劃分方法，能為將來破壞性地震發生后獲取烈度分布提供參考，其現實意義是顯著的。

地震動的模擬結果取決于震源、傳播路徑及場地3個物理過程的模型建立，本文中震源的斷層滑動模型是根據McGinty和Robinson(1999)計算靜態庫侖應力觸發時所采用的平均錯動方法而建立的，模型的簡化可能會造成計算結果的一些誤差;在對合成地震動進行場地效應校正時，采用的是區域平均場地效應值，這與場點的真實場地響應有一定差異。此外，本文只計算了一種破裂方法的地震，在結果的普遍性上受到了限制，因此在今后的研究中還需要對不同震源模型、破裂模式、震級規模的地震進行研究。

崔建文，盧大偉，高東，等.2008.基于合成地震的震區烈度劃分［J］.地震研究，31(4):387－393.

胡連英，徐學思，孫壽成，等.1997.溧陽地震與茅東斷裂帶［M］.北京:地震出版社.

黃耘，李清河，孫業君，等.2006.江蘇及鄰區地殼上地幔結構研究［J］.西北地震學報，28(4):369－376.

李山有，金星，陳先，等.2002.地震動強度與地震烈度速報研究［J］.地震工程與工程震動，22(6):1－7.

林邦慧，魏富勝，劉萬琴，等.1982.溧陽—介林—五原北西地震帶強震的破裂特征［J］.地震學報，4(2):14－24.

石玉成，陳厚群，李敏，等.2005.隨機有限斷層法合成地震動的研究與應用［J］.地震工程與工程震動，25(4):18－23.

孫受成，徐學思，胡連英.1980.溧陽地震［J］.地震，1(1):32－36.

王海云，謝禮立.2009.近斷層強地震動場預測［J］.地球物理學報，52(3):703－711.

王煒.1984.溧陽地震前震中分布集中度C的變化［J］.地震科學研究，2(2):89－94.

袁一凡.1998.由地震動三要素確定地震動強度(烈度)的研究［R］.北京:國家地震局工程力學研究所.

詹小艷，許紅梅，王俊，等.2011.江蘇地區的介質非彈性衰減和場地響應研究［J］.地震地磁觀測與研究，32(6)，20－27.

Aki K，Richaids P G.1986.定量地震學［M］.北京:地震出版社.

Atkinson G M，Beresenev I A.2002.Ground motion at Memphis and St.Louis from M7.5－8.0 earthquakes in the New Madrid seismic zones［J］.BSSA，92(3):1 015 － 1 024.

Beresnev I A，Atkinson G M.1997.Modeling finite-fault radiation from the ωn spectrum［J］.BSSA，87(1):67 －84.

Beresnev I A，Atkinson G M.1998.FINSIM －a FORTRAN program for simulating stochastic acceleration time histories from finite faults［J］.Seism Res Let，69(1):27 － 32.

Beresnev I A，Atkinson G M.2002.Source parameters of earthquakes in eastern and western North America based on finite-fault modeling［J］.BSSA，92(2):695 －710.

Boore D M.1983.Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra［J］.BSSA，73(6A):1 865－1 894.

Brune J.1970.Tectonic Stress and spectra of seismic shear waues from earthquake［J］.JGR，75(26):4997 ～ 5009.

Hanks T C，McGuire P K.1981.The character of high frequency strong ground motion［J］.BSSA，71(6):2 071 －2 095.

Hartzell S H.1978.Earthquake aftershocks as Green's functions［J］.Geophys Res Lett，5(1):1 － 4.

Irikura K.1983.Semi-empirical estimation of strong ground motions during large earthquakes［J］.Bull Disaster Prevention Res Inst(Kyoto Univ)，33(298):63 －104.

McGinty P，Robinson R.1999.Silp distribution of Lake Tennyson earthquake，New Zealand，as inferred from static stress changes and off fault aftershocks［J］.Geophysical Research Letters，26(13):1 961－1 964.

Motazedian D，Atkinson G M.2005.Stochastic finite-fault modeling based on a dynamic corner frequency［J］.BSSA，95(3):995 － 1 010.

Somerville P，Sen M，Cohee B.1991.Simulations of strong ground motions recorded during the 1985 Michocan，Mexico and Valparaiso，Chile，earthquake［J］.BSSA，81(1)，1 －27.

Zafeiria R，Beresnev I A.2003.Stochastic finite-fault modeling of ground motions from 1999 Chi-Chi，Taiwan，earthquake:Application to Rock and Soil sties with Implications for Nonlinear Site Response［J］.BSSA，93(4):1 691 －1 702.

Computational Research on the Seismic Intensity Based on Stochastic Finite Faults Method

WANG Jun，HUO Zhu-qing，ZHAN Xiao-yan，SONG Hao，JIANG Hao-ling，XU Ge
(Earthquake Administration of Jiangsu Province，Nanjing 210014，Jiangsu，China)

Firstly，using stochastic finite faults method based on the dynamic corner frequency，we simulated the synthetic ground motion of Dongtai M3.8 earthquake，Jiangsu in 2006 recorded by Dongtai，Rugao and Dafeng stations and compared synthetic ground motion with the actual ground motion recordings.We find that the synthetic ground motion using improved stochastic finite faults method can better reflects the main features of shear wave of main shock in Jiangsu area.Secondly，we apply this method to simulate the synthetic ground motion of 156 virtual observation points in source region and got the distribution of synthetic acceleration field of Liyang M6.0 earthquake in 1976.Thirdly，we respectively calculated the theoretical distribution of seismic intensity by peak acceleration statistical method and fuzzy evaluation method.The results showed that the intensity classification results by two methods essentially consistent with that of the actual investigation.However，intensity classification result by the fuzzy evaluation method is more consistent with distribution of actual intensity.

dynamic corner frequency;stochastic finite faults method;synthetic ground motion;seismic intensity;Liyang M6.0 earthquake

P315.9

A

1000－0666(2012)03－0374－07

2011－08－18.

國家自然科學基金項目 (41074036)和江蘇省地震局青年基金 (10420)聯合資助.