九寨溝MS7.0地震強(qiáng)地震動(dòng)模擬及參數(shù)敏感性分析

李明會(huì) 景月嶺,2 張宇璇 鄭司元 宛明晨 何李浩

1 合肥工業(yè)大學(xué)土木與水利工程學(xué)院，合肥市屯溪路193號(hào)，230001 2 水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京市清華園1號(hào)，100089

2017-08-08四川省九寨溝縣發(fā)生MS7.0地震。關(guān)于九寨溝地震強(qiáng)地面運(yùn)動(dòng)數(shù)值模擬的研究取得了一些成果[1-4]，但仍有不足之處，如應(yīng)用復(fù)合震源模型進(jìn)行模擬及其震源參數(shù)敏感性分析的研究較少、隨機(jī)有限斷層法計(jì)算地震動(dòng)分量方向單一、地震峰值加速度擬合精度不高等。為此，本文應(yīng)用復(fù)合震源模型對(duì)九寨溝MS7.0地震強(qiáng)地面運(yùn)動(dòng)數(shù)值模擬展開研究，并探討不同模型參數(shù)對(duì)合成地震動(dòng)的影響規(guī)律，為九寨溝地區(qū)工程場(chǎng)地地震安全評(píng)價(jià)提供參考依據(jù)。

1 復(fù)合震源模型基本理論

復(fù)合震源模型是用于計(jì)算寬頻帶強(qiáng)地面運(yùn)動(dòng)的有限源模型[5]。該模型通過描述震源的運(yùn)動(dòng)學(xué)特征，結(jié)合地震波在層狀介質(zhì)中的傳播過程，將廣義反射和透射系數(shù)法[6]得到的格林函數(shù)與震源時(shí)間函數(shù)進(jìn)行卷積，計(jì)算得到寬頻帶地震動(dòng)。模型中，強(qiáng)震震源被認(rèn)為是大量具有恒定應(yīng)力降圓形子震的疊加。首先將斷層面劃分成一系列尺寸相同的正方形子斷層網(wǎng)格，之后用一定數(shù)量和相應(yīng)半徑大小分布服從Frankel自相似模型[7]的圓形子震填充該斷層面，并允許其部分重疊。圓形子震相當(dāng)于一系列小震，在震源破裂過程中相繼被激發(fā)，最終釋放的能量在斷層面上不斷被疊加，達(dá)到模擬強(qiáng)震的目的。設(shè)圓形子震半徑為R、數(shù)目為N，則有：

(1)

式中，分形維數(shù)D=2，控制斷層面上圓形子震數(shù)量的分布及圓形子震的上升時(shí)間，P為比例常數(shù)：

(2)

式中，M0為地震矩，Δσ為應(yīng)力降，Rmax為最大圓形子震半徑，Rmin為最小圓形子震半徑。對(duì)式(1)進(jìn)行積分，可得半徑大于R且小于Rmax圓形子震的分布數(shù)量N(R)：

(3)

(4)

(5)

(6)

式中，β為剪切波的速度。斷層面上滑動(dòng)位移量的空間分布由圓形子震滑動(dòng)位移量疊加得到，第i個(gè)圓形子震的滑動(dòng)位移量ΔUi為[1]：

(7)

式中，C為依賴斷層類型的常數(shù)，μ為介質(zhì)剪切模量。

2 地震數(shù)據(jù)與震源參數(shù)

本文研究區(qū)域(103°～105°E，31°～35°N)位于巴顏喀拉塊體，該塊體近年來地震頻發(fā)。九寨溝地震后，中國地震臺(tái)網(wǎng)中心等機(jī)構(gòu)公布了其震源位置、地震矩、震源深度、斷層節(jié)面的走向、傾角和滑動(dòng)角等震源參數(shù)(表1)，其中，根據(jù)美國地質(zhì)調(diào)查局的研究，斷層節(jié)面Ⅱ?yàn)閷?shí)際斷層面[1]，因此本文取節(jié)面Ⅱ進(jìn)行研究。

表1 九寨溝地震震源參數(shù)

基于孟令媛等[9]對(duì)余震分布情況及斷層面上滑動(dòng)位移的反演結(jié)果，建立長30 km、寬20 km的斷層模型，并將其劃分為150個(gè)2 km×2 km的正方形子斷層。斷層面上圓形子震的分布如圖1所示，其中圓形子震地震矩的總和等于總地震矩。通過式(7)計(jì)算得到斷層面上滑動(dòng)位移量的不均勻分布，如圖2所示。九寨溝地震震中距150 km區(qū)域范圍內(nèi)臺(tái)站信息如表2所示，九寨溝地區(qū)介質(zhì)速度信息(圖3)來自全球地殼模型Crust1.0(https:∥igppweb.ucsd.edu/～gabi/crust1.html)，臺(tái)站位置分布如圖4所示。對(duì)所有臺(tái)站觀測(cè)數(shù)據(jù)進(jìn)行基線校準(zhǔn)和濾波處理，帶通濾波范圍為0.1～20 Hz。

圖1 斷層面上圓形子震分布Fig.1 Distribution of circular subevents on the fault plane

圖2 斷層面上滑移量的不均勻分布Fig.2 Non-uniform distribution of slip on the fault plane

表2 臺(tái)站信息

圖3 九寨溝地區(qū)速度模型Fig.3 Velocity model of the Jiuzhaigou region

圖4 臺(tái)站分布Fig.4 Distribution of the stations

3 參數(shù)敏感性分析

在運(yùn)動(dòng)學(xué)震源模型中，斷層的破裂面積、長度、寬度、埋置深度、走向、傾角和滑移角等是描述斷層特征的參數(shù)；地震矩、應(yīng)力降、平均滑動(dòng)時(shí)間、平均破裂速度、破裂方式和初始破裂點(diǎn)位置等是描述平均破裂過程的參數(shù)。這些震源參數(shù)控制著地震動(dòng)強(qiáng)度和分布[10]，分析這些參數(shù)的影響規(guī)律可為地震動(dòng)合成精度提供保證。本文首先以無強(qiáng)震觀測(cè)記錄為前提，研究地震矩、應(yīng)力降、破裂速度、子斷層尺寸、震源深度、斷層走向、傾角和滑移角等8個(gè)震源參數(shù)在所有合理取值范圍內(nèi)對(duì)復(fù)合震源模型計(jì)算結(jié)果的影響程度；然后根據(jù)九寨溝地震觀測(cè)數(shù)據(jù)確定適合的震源參數(shù)研究范圍，并進(jìn)行敏感性規(guī)律分析，在其他參數(shù)保持不變的條件下，僅改變單個(gè)參數(shù)變量，利用評(píng)價(jià)指標(biāo)比較分析計(jì)算結(jié)果，總結(jié)震源參數(shù)對(duì)地震動(dòng)計(jì)算的影響規(guī)律。

3.1 評(píng)價(jià)參數(shù)和評(píng)價(jià)指標(biāo)

為表征震源參數(shù)對(duì)臺(tái)站NS和EW向峰值加速度(peak ground acceleration，PGA)的影響，用誤差值E表示計(jì)算的PGA相對(duì)于觀測(cè)的PGA的變化：

E=|AOBS-ASYN|/AOBS×100%

(8)

式中，AOBS為觀測(cè)的PGA，ASYN為計(jì)算的PGA。

地震動(dòng)強(qiáng)度主要用PGA、峰值速度或譜加速度等地震動(dòng)參數(shù)表達(dá)。本文以PGA為例，根據(jù)王德才等[11]的參數(shù)敏感性分析模型，通過計(jì)算模型偏差和模型標(biāo)準(zhǔn)差討論不同震源參數(shù)對(duì)PGA的影響。

用模型偏差B表征計(jì)算值和觀測(cè)值之間的差異，偏差越大，計(jì)算值越偏離觀測(cè)數(shù)據(jù)：

(9)

式中，M為總臺(tái)站數(shù)。

用模型標(biāo)準(zhǔn)差S描述數(shù)據(jù)離散程度，標(biāo)準(zhǔn)差越大，數(shù)據(jù)波動(dòng)性越強(qiáng)：

(10)

3.2 參數(shù)敏感性分析結(jié)果

首先選取地震矩、應(yīng)力降、破裂速度、子斷層尺寸、震源深度、斷層走向、傾角和滑移角等8個(gè)震源參數(shù)的最大合理取值范圍，計(jì)算得到不同震源參數(shù)數(shù)值下臺(tái)站NS和EW向的PGA，然后以模型偏差B和標(biāo)準(zhǔn)差S為衡量標(biāo)準(zhǔn)，分析8個(gè)震源參數(shù)對(duì)復(fù)合震源模型計(jì)算結(jié)果的影響程度。

地震矩根據(jù)Hanks等[12]提出的地震矩與矩震級(jí)的關(guān)系，在強(qiáng)震6～8級(jí)范圍內(nèi)取值計(jì)算，間隔為1×1018Nm；應(yīng)力降參考王振宇等[13]的取值范圍10～200 bar，間隔為5 bar；斷層的平均破裂速度約為剪切波速的0.6～0.9倍[14]，本文選取此區(qū)間，并設(shè)定區(qū)間間隔為0.05；子斷層尺寸取值分別為1 km×1 km、2 km×2 km、6 km×6 km、8 km×8 km和10 km×10 km；震源深度范圍選定為0～30 km，取值間隔為1 km；斷層走向、傾角和滑移角依據(jù)理論范圍分別取0°～360°、0°～90°和-180°～180°，取值間隔分別為10°、5°和10°。

8個(gè)震源參數(shù)取值范圍和梯度各不相同，本文通過式(9)和式(10)計(jì)算出相應(yīng)PGA的偏差和標(biāo)準(zhǔn)差，然后將每個(gè)震源參數(shù)下PGA偏差、標(biāo)準(zhǔn)差的最大值及最小值作差，消除AOBS實(shí)測(cè)數(shù)據(jù)的影響，得到每個(gè)震源參數(shù)下PGA偏差和標(biāo)準(zhǔn)差的變化幅度，結(jié)果如圖5所示。由圖5可知，地震矩、應(yīng)力降、破裂速度、子斷層尺寸和震源深度對(duì)PGA偏差和標(biāo)準(zhǔn)差的影響明顯，其中地震矩和應(yīng)力降影響更突出；斷層走向、傾角和滑移角對(duì)PGA偏差和標(biāo)準(zhǔn)差的影響極不明顯。根據(jù)地震學(xué)理論，地震矩直接決定地震釋放的能量大小，地震矩越大，地震波攜帶的能量越多，使得整體加速度計(jì)算值增加；應(yīng)力降表示斷層上某點(diǎn)破裂前后的應(yīng)力強(qiáng)度變化，應(yīng)力降越大，斷層上釋放的應(yīng)力越多；破裂速度決定斷層面圓形子震的發(fā)震時(shí)間，間接影響整個(gè)斷層面的破裂過程；子斷層尺寸的劃分會(huì)對(duì)斷層面上滑移量分布產(chǎn)生很大的影響；震級(jí)相同的地震，震源越淺，造成的破壞越嚴(yán)重。

圖5 不同震源參數(shù)對(duì)PGA偏差和標(biāo)準(zhǔn)差的影響幅度Fig.5 The influence amplitude of different sourceparameters on PGA deviation and standard deviation

基于8個(gè)震源參數(shù)的影響規(guī)律，結(jié)合表1不同研究機(jī)構(gòu)和眾多學(xué)者[2-4,8]提供的震源參數(shù)，確定九寨溝地震震源參數(shù)研究范圍，進(jìn)一步研究8個(gè)震源參數(shù)對(duì)九寨溝11個(gè)臺(tái)站NS和EW向PGA的影響大小和規(guī)律，取值間隔如圖6所示。

從圖6看出，地震矩、應(yīng)力降、破裂速度、子斷層尺寸和震源深度對(duì)臺(tái)站PGA的影響大；斷層走向、傾角和滑移角對(duì)臺(tái)站PGA的影響小。說明在應(yīng)用復(fù)合震源模型計(jì)算九寨溝地震時(shí)，即使8個(gè)震源參數(shù)取值范圍有所限制，沒有取到理論范圍內(nèi)的所有值，但其對(duì)計(jì)算結(jié)果的影響規(guī)律基本與圖5一致。

圖6 不同震源參數(shù)對(duì)臺(tái)站PGA的影響Fig.6 Influence of different source parameters on station PGA

圖7為地震矩、應(yīng)力降、破裂速度、子斷層尺寸和震源深度等5個(gè)震源參數(shù)對(duì)九寨溝地震地震動(dòng)模擬的影響規(guī)律，其中實(shí)心點(diǎn)表示本文取值。由表3和圖7可知，受地震矩影響，PGA的偏差變化范圍超過1，PGA的偏差和標(biāo)準(zhǔn)差隨地震矩的增加而減小，其中PGA的偏差下降趨勢(shì)近似為線性變化，地震矩為7×1018Nm時(shí)，PGA的偏差基本為0，同時(shí)PGA的標(biāo)準(zhǔn)差也較小。應(yīng)力降越大，臺(tái)站整體PGA越大，PGA的偏差越小，應(yīng)力降為40 bar時(shí)，PGA的偏差接近于0，這與俞瑞芳等[15]應(yīng)用隨機(jī)有限斷層法得到PGA隨著應(yīng)力降增大而逐漸增加的結(jié)論類似。破裂速度越大，PGA的偏差越小，當(dāng)破裂速度等于剪切速度0.8倍(2.96 km/s)時(shí)，偏差接近于0，此時(shí)標(biāo)準(zhǔn)差也最小，與傅磊等[16]應(yīng)用隨機(jī)有限斷層法和李宗超等[4]應(yīng)用經(jīng)驗(yàn)格林函數(shù)法在地震模擬過程中選取的破裂速度與剪切速度的比例關(guān)系一致，PGA的偏差變化范圍為0.9。子斷層尺寸越小，PGA的偏差越小，最佳尺寸為2 km×2 km，此時(shí)標(biāo)準(zhǔn)差也為最小值，與Dang等[1]利用隨機(jī)有限斷層法的研究結(jié)論相同。震源深度在20 km處PGA的偏差幾乎為0，同時(shí)標(biāo)準(zhǔn)差最小，隨著震源深度的增加，PGA的偏差增長趨勢(shì)變緩；而且如圖6所示，震源深度為25 km時(shí)，大部分臺(tái)站的峰值加速度很小，與劉啟方等[10]運(yùn)用震源位錯(cuò)模型得到的研究結(jié)果相似。

表3 不同震源參數(shù)PGA偏差和標(biāo)準(zhǔn)差統(tǒng)計(jì)結(jié)果

圖7 不同震源參數(shù)對(duì)PGA偏差和標(biāo)準(zhǔn)差的影響Fig.7 Influence of different source parameters on PGA deviation and standard deviation

斷層走向、傾角和滑移角中，PGA的偏差最大變化幅度為0.08，標(biāo)準(zhǔn)差最大變化幅度僅為0.1，對(duì)計(jì)算結(jié)果影響不明顯，與王德才等[11]采用混合方法[17]和隨機(jī)方法[18]研究的結(jié)果一致。

3.3 地震輸入?yún)?shù)組合

由上文分析結(jié)果可知，震源參數(shù)整體對(duì)PGA偏差影響幅度大，對(duì)標(biāo)準(zhǔn)差影響幅度相對(duì)較小，因此以PGA的偏差為主要衡量標(biāo)準(zhǔn)，標(biāo)準(zhǔn)差為次要標(biāo)準(zhǔn)，得到一組與觀測(cè)記錄擬合效果較好的震源參數(shù)。從表1各機(jī)構(gòu)提供的數(shù)據(jù)中選取地震矩、破裂速度、子斷層尺寸、震源深度、斷層走向、傾角和滑移角；應(yīng)力降參考Sun等[3]的結(jié)果，品質(zhì)因子Q(f)和高頻衰減kappa值參考孟令媛等[9]的結(jié)果，具體數(shù)值如表4所示。

表4 九寨溝地震復(fù)合震源模型輸入的震源物理參數(shù)

4 地震動(dòng)模擬結(jié)果分析

根據(jù)式(8)計(jì)算11個(gè)臺(tái)站NS和EW向PGA的相對(duì)誤差，并分別與Dang等[1]和李宗超等[4]的計(jì)算結(jié)果進(jìn)行對(duì)比(圖8)?？梢钥闯?，11個(gè)臺(tái)站NS和EW向PGA相對(duì)誤差均小于50%。62SHW站NS和EW向PGA相對(duì)誤差最小，均小于3%；51HSL站NS向、51HSD站NS和EW向PGA誤差較大，且計(jì)算值均大于觀測(cè)值，主要原因可能是臺(tái)站震中距較大，導(dǎo)致傳播過程中散射的復(fù)雜性增加，同時(shí)，臺(tái)站局部山體地形引起的場(chǎng)地效應(yīng)也會(huì)增加計(jì)算誤差。11個(gè)臺(tái)站NS和EW向PGA相對(duì)誤差均在20%以內(nèi)，說明所有臺(tái)站PGA與實(shí)際數(shù)據(jù)相比誤差小。與Dang等[1]和李宗超等[4]的計(jì)算結(jié)果相比，除李宗超等[4]51JZY站EW向、62ZHQ站NS向和62MXT站NS向外，本文各臺(tái)站PGA相對(duì)誤差均較優(yōu)。為方便直觀地進(jìn)行地震動(dòng)持時(shí)和波形對(duì)比，所有時(shí)程均截取45 s時(shí)間窗(圖9)，臺(tái)站加速度記錄數(shù)據(jù)中最小PGA為2.75 cm/s2，為避免地震動(dòng)持時(shí)出現(xiàn)零持時(shí)這種不合理狀態(tài)，選取Bolt持時(shí)的閾值為0.5倍的PGA 作為計(jì)算相對(duì)Bolt持時(shí)的閾值[19]。從圖9看出，個(gè)別臺(tái)站誤差較大，但大多數(shù)臺(tái)站地震動(dòng)持時(shí)基本保持一致。如圖10所示，合成反應(yīng)譜整體變化趨勢(shì)與觀測(cè)反應(yīng)譜一致，說明11個(gè)臺(tái)站的峰值加速度、地震動(dòng)持時(shí)和反應(yīng)譜整體擬合效果較好，與觀測(cè)數(shù)據(jù)較為符合。

圖8 不同臺(tái)站NS和EW向PGA的相對(duì)誤差對(duì)比Fig.8 Comparison of relative errors of PGA in NS and EW directions at different stations

圖9 不同臺(tái)站NS和EW向加速度時(shí)程Fig.9 Time history of acceleration in NS and EW directions at different stations

圖10 不同臺(tái)站合成和觀測(cè)加速度反應(yīng)譜Fig.10 Synthetic and observed acceleration response spectra at different stations

部分臺(tái)站地震動(dòng)三要素計(jì)算結(jié)果精度偏低的可能原因是：1)沒有充分考慮地形因素對(duì)地震動(dòng)的影響，部分臺(tái)站如51HSL、51HSS和51HSD等位于山體上，山體地形的高度及底部延展距離對(duì)斷層動(dòng)力學(xué)破裂過程影響較大，進(jìn)而影響到相應(yīng)的地面地震動(dòng)分布[20]，山體的坡角、地震動(dòng)入射角度會(huì)對(duì)山體地震動(dòng)產(chǎn)生抑制作用，導(dǎo)致臺(tái)站PGA變小[21]；2)本文在計(jì)算地震動(dòng)時(shí)采用各向同性散射模型，但實(shí)際地震波在地殼介質(zhì)中的傳播是各向異性的，而且隨著傳播路徑的增加，散射的影響也會(huì)逐漸加強(qiáng)；3)有些震源參數(shù)的選取進(jìn)行了簡化，在地震動(dòng)模擬過程中假定某些參數(shù)(如應(yīng)力降和破裂速度)為常數(shù)，忽略了參數(shù)與時(shí)間和空間的相關(guān)性；4)多數(shù)臺(tái)站低頻段反應(yīng)譜模擬結(jié)果高于觀測(cè)值，可能是由于本文模型在破裂過程中沒有考慮到圓形子震傾角的可變性，或是在合成格林函數(shù)中未納入小尺度地殼非均勻性的散射效應(yīng)[22]。

PGA是表征地震能量的重要指標(biāo)之一，能夠?yàn)榻ㄖ锟拐鹪O(shè)防提供合理的參考依據(jù)。綜合現(xiàn)有11個(gè)臺(tái)站PGA的觀測(cè)值和計(jì)算值可以得到九寨溝地震PGA 空間分布特征，如圖11所示。由圖11可知，隨著震中距增大，NS和EW向PGA從震源中心呈橢圓形向四周衰減擴(kuò)散，長軸方向大致呈SW走向，這可能與斷層破裂的方向性有關(guān)，且2個(gè)方向PGA最大值均出現(xiàn)在51JZB站；在51MXD和51PWM站之間有一個(gè)PGA值未連通的區(qū)域，可能是由于該區(qū)域缺少合適的臺(tái)站記錄所致。總體而言，PGA計(jì)算值和觀測(cè)值的空間分布特征較為一致。

圖11 九寨溝地震PGA空間分布特征Fig.11 PGA spatial distribution characteristics of Jiuzhaigou earthquake

5 結(jié) 語

1)在最大合理變化范圍內(nèi)，相對(duì)于斷層走向、傾角和滑移角，地震矩、應(yīng)力降、破裂速度、子斷層尺寸和震源深度對(duì)地震動(dòng)計(jì)算結(jié)果影響較大。

2)震源參數(shù)敏感性分析規(guī)律為九寨溝地震地震動(dòng)模擬的參數(shù)選取提供了依據(jù)，進(jìn)而得到一組合理的震源參數(shù)組合，同時(shí)能夠獲得九寨溝地區(qū)內(nèi)任意一點(diǎn)PGA的變化區(qū)間。這一結(jié)論可以進(jìn)一步推廣到其他缺乏強(qiáng)震觀測(cè)記錄的場(chǎng)地，計(jì)算出該區(qū)域強(qiáng)度特征的上限值，為建筑物抗震設(shè)防提供依據(jù)。

3)復(fù)合震源模型強(qiáng)地震動(dòng)模擬了特定圓形子震在斷層上的破裂過程、斷層面上滑動(dòng)位移的不均勻空間分布以及震源時(shí)間函數(shù)，計(jì)算結(jié)果可以再現(xiàn)九寨溝地震的地震動(dòng)特征，模擬區(qū)域臺(tái)站的地震動(dòng)PGA、地震動(dòng)持時(shí)、反應(yīng)譜和PGA分布特征與觀測(cè)值擬合較好，反應(yīng)譜在高頻段擬合精度更高。