隨機振動有限斷層方法計算衡水市強地面運動參數(shù)

劉博研，申文豪，史保平

（1.中國地震局地殼應力研究所（地殼動力學重點實驗室），北京 100085；2.中國科學院大學，北京 100049）

0 引言

研究強地面運動對于減少震害，保護人民生命財產(chǎn)安全具有重要意義。強地面運動研究的科學目的是揭示地震引起地震動的機理，為強地震動的預測提供理論基礎。20世紀90年代末至今的二十幾年間，強地面運動觀測記錄累積得很快，包括1989年 美 國Loma PrietaMW7.0 地 震、1994年 美 國NorthridgeMW6.7 地震、1999年中國臺灣MW7.6地震以及2008年汶川MW7.9地震、2010年墨西哥BajaMW7.2地震等，每一次地震都積累了幾十乃至上百條強震記錄，這些記錄對于獲得強地面運動規(guī)律具有重要的科研價值。盡管這些強震記錄對強震觀測數(shù)據(jù)庫做了極大的補充，但是目前仍然沒有從根本上解決強地面運動觀測缺乏的問題。

由于中強地震觀測數(shù)據(jù)缺乏，發(fā)展強地面運動的模擬技術已成為定量化估計強地面運動的重要手段。地震學家在建立模型的探索中發(fā)現(xiàn)，模型建立過程必須把震源破裂過程、地震波傳播過程以及場地效應作為一個系統(tǒng)來考慮。準確地預測強地面運動，不僅對了解地震的震源物理和力學過程有重要的理論價值，而且在工程抗震和地震災害的定量化評估、災后救援、政府決策中有著重要的實際應用價值。強地面運動預測還是地震風險評估和未來潛在地震危險性區(qū)劃工作的重要組成部分。利用區(qū)域歷史地震活動性資料、地質(zhì)（斷層）資料、活動斷層探測成果、GPS（大地測量）數(shù)據(jù)以及其他地球物理資料，結(jié)合計算機網(wǎng)絡化GIS可視化技術，實施基于物理的地震動數(shù)值模擬技術，可為未來中強地震潛在危險性分析提供合理參數(shù)，為區(qū)域地震風險估算提供新的技術支撐平臺，彌補傳統(tǒng)方法的局限性。未來潛在中強地震活動區(qū)內(nèi)地面運動的模擬計算成果可為區(qū)域內(nèi)地震危險性分析提供合理的約束。

事實上，地震波的產(chǎn)生與傳播是極其復雜的過程，真實的地震波并沒有固定的振動周期，其振幅也是隨時間變化的［1］。在地震工程學中，常常把任意一條真實的地震波看成有限持時的隨機振動波，對于任何實測的地震波每一個頻率均有相應的傅里葉譜值［2］。很多學者在模擬地震波的過程中考慮了這種隨機性［3－4］，在地震工程中利用隨機理論模擬強震地面運動有了很大的發(fā)展。隨機方法的基本原理是把遠場地震動看成是有限帶寬、有限持時的高斯白噪聲，并把影響地震動參數(shù)的一系列因素（包括震源、路徑和場地）在頻率域里表達成簡單的乘積形式，然后與通過窗函數(shù)調(diào)整的噪聲傅里葉譜結(jié)合，最后通過傅里葉反變換得到地震動參數(shù)的時程曲線。隨機方法的震源理論模型是Brune［5］在1970年首次提出的垂直斷層面輻射剪切波的震源模型；Hanks和McGuire［6］等人在1981年提出了隨機振動理論計算峰值加速度和峰值速度的公式，為隨機性模擬地震動峰值做出了理論上的貢獻；1983年，Boore［4］在Hanks 和McGuire［6］等 人 的 工 作 基 礎上，開創(chuàng)性的將蒙特卡羅方法應用其中，給出了隨機點源模型，實現(xiàn)了地震動時程及反應譜的估計；1984年，Atkinson［7］在Hanks和McGuire［6］與Boore［4］等人的研究基礎上給出了三段幾何衰減模型，進一步修正了隨機振動方法在遠場的應用。至此，隨機性方法的框架已經(jīng)基本構建完善。若以震源類型進行分類，可將計算模型簡單分為點源模型和有限斷層模型，前者將震源簡化為一個點源，這類模型多用于對強地面運動的簡單估算及遠場強地面運動的模擬。而在近斷層范圍內(nèi)，地面運動直接同斷層破裂過程相關，斷層的走向、傾角、滑移方向等都會對地表運動的形態(tài)產(chǎn)生影響。因此，采用點源模型模擬強地面運動的時空變化特征顯然無法滿足工程抗震設計和地震災害預測的需要。對于中強地震采用有限斷層模型更為合理。

針對強地面運動的研究涉及多種震源理論和斷層模型，在本研究中，我們著重介紹用隨機振動的有限斷層模型來計算衡水市強地面運動參數(shù)以及設定地震的時程曲線。

1 方法原理

1.1 隨機振動有限斷層

Beresnev和Atkinson［8］等對隨機點源方法做了改進，提出隨機有限斷層方法。其基本思路是將斷層面分成N個大小相等的矩形子斷層，每個子斷層即為1個點源。通過設定震源模式和破裂傳播速度，就可以得到子源破裂的時間順序，然后根據(jù)子源與場地的幾何關系，計算每個子源對場地的影響，即可疊加成場地的地震動時程。該方法在我國得到了廣泛的應用，李啟成［9］比較了隨機點源和隨機有限斷層2種方法，指出2種方法的適用性和局限性，并模擬了1982年ML5.1地震和1976年MS6.9寧河地震。王國新和史家平［10］、王海云［11］利用隨機有限斷層模型分別進行了汶川地震和1994年Northridge地震的強地震動模擬，模擬結(jié)果與臺站記錄一致性較好。陶夏新和王國新［12］論證了這種隨機方法在合成的時程中破裂的方向性效應和上盤效應的有效性。在國外，隨機振動的有限斷層模型也已被廣泛應用于已發(fā)地震的強地面運動模擬，以及對未來地震或者設定地震強地面運動的預測。Atkinson和Boore［13］用改進后的隨機有限斷層模型模擬了美國東部地區(qū)的強地面運動，結(jié)合模擬數(shù)據(jù)和臺站記錄發(fā)展了美國東部的衰減關系。Chopra等［14］基于隨機振動的有限斷層模型模擬了2001年MW7.6印度Bhuj大地震強地面運動參數(shù)，并利用模擬得到PGA 計算出了理論MMI的分布圖。Zafarani等［15］用該模型評估了伊朗8個重要地震，并預測了大德黑蘭地區(qū)未來地震的強地面運動參數(shù)。

在隨機有限斷層模型中，斷層破裂面被分成N個大小相等的矩形子源，每個子源即為1個點源，每個子源滿足ω2衰減規(guī)律。通過設定震源模式和破裂傳播速度，可以得到子源破裂的時間順序，然后根據(jù)子源與場地的幾何關系，計算每個子源對場地的影響。所有子源在觀測點引起的地震動在時域中以適當?shù)难舆t時間疊加，即可獲得場地的地震動時程a（t）：

式中：NL、NW分別為斷層沿長度和寬度劃分的子斷層數(shù)；i和j分別代表沿斷層走向和傾向方向的斷層指標（Index）；Δtij為子源（下標為ij）到場地的滯后時間；aij（t）為子源（下標為ij）利用隨機點源方法得到的在觀測點的加速度。

在有限斷層中，為了保證斷層輻射的總輻射能量守恒，子斷層上劃分有嚴格的要求，Atkinson和Silva［16］指出用隨機模型合成地震動時，子斷層大小為：

式中：Δl為子斷層尺度，單位為km；M為地震震級。斷層總面積除以子斷層的面積就是斷層上子斷層的個數(shù)。由于子斷層的大小有嚴格的規(guī)定，單次觸發(fā)的子斷層地震距的總和往往小于總地震距。為了保證總地震距守恒，常常采用部分子斷層多次觸發(fā)的方式。針對有限斷層模型的這些不合理之處，Motazedian和Atkinson［17］提出了“動力學拐角頻率”的概念。他們認為，在斷層面破裂過程中拐角頻率是一個時間的函數(shù)，與破裂面積成反比，破裂從高拐角頻率開始，隨著破裂面積的擴大，拐角頻率逐漸降低，子斷層之間相互獨立。動力學拐角頻率為：

式中：fcij（t）為第ij個子斷層的動力學拐角頻率；t是第ij個子源被觸發(fā)的時刻；NR（t）是在時刻t已破裂子斷層的累計數(shù)；M0ave＝M0／N，是子斷層的平均地震矩；S是表示子斷層輻射強度的一個常數(shù)。隨著破裂的傳播，破裂的子斷層數(shù)增加，子斷層拐角頻率減小。為了使斷層總輻射能守恒，引入一個標度因子H，動力學拐角頻率對子斷層總輻射能的遞減影響到標度因子H的補償，總輻射能保持不變。這樣，在保證輻射能守恒的前提下避免了子斷層多次觸發(fā)這一非物理現(xiàn)象。至此，隨機振動的有限斷層方法的理論框架基本得到完善。

1.2 淺層橫波速度V30對地震動模擬結(jié)果的影響

在宏觀地震震害調(diào)查中發(fā)現(xiàn)，一般軟弱地基上的震害較堅硬地基上的震害重。隨著進一步研究，人們認識到淺層速度結(jié)構對場地強地面運動有很大影響，表現(xiàn)為對地震動的放大或縮小，常被看作引起場地效應的重要原因，直接影響到地震災害的嚴重程度。因此，在數(shù)值計算中淺層速度結(jié)構對強地面運動的影響必須考慮。

在隨機有限斷層模型中，若要得到淺層橫波速度對某場點地震動的影響，必須將該點淺層各個深度的速度結(jié)構加入到程序中運算。對于大規(guī)模、高精度的網(wǎng)格運算而言，由于構成場地的巖土及其所含各種結(jié)構面是在漫長的地質(zhì)歷史中形成的，并經(jīng)歷了復雜的地殼構造運動的改造，所以在一般情況下場地條件都比較復雜，顯然不可能獲得每個場點不同深度的淺層速度結(jié)構。Wald和Allen［17］根據(jù)全球地形數(shù)據(jù)給出了全球范圍內(nèi)淺層速度結(jié)構數(shù)據(jù)V30，V30不是嚴格意義的地表30 m 以下的橫波速度，而是全球平均意義下的淺層橫波速度，能夠反映地表淺層橫波速度在橫向上的變化，因此，可以用來進行強地面運動的模擬。V30數(shù)據(jù)范圍為180～760 m／s，通常認為760 m／s為一般情況下基巖的速度結(jié)構。接下來的問題就是如何處理V30對地震動參數(shù)的影響。Boore和Atkinson［19］在建立衰減關系NGA（Next Generation Attenuation）時將地震動參數(shù)Y 的構成寫作（4）式：

式中：FM、FD和FS分別代表震級影響、距離影響和場地放大效應。FS是由線性部分FLIN和非線性部分FNL兩部分組成，并且線性部分占絕大部分，非線性部分可以忽略不計。如果地震動參數(shù)Y為PGA，那么FLIN＝－0.36ln（V30／760）；如果Y是PGV，則FLIN＝－0.6ln（V30／760）。如果將基巖上的強地面運動參數(shù)設為Y，那么增加了場地放大效應后的強地面運動參數(shù)可設為Y′，可得：

因此，如果我們能夠獲取計算區(qū)域內(nèi)網(wǎng)格點上的V30分布，那么增加場地影響后的地震動參數(shù)為：

這樣，我們就能將V30的影響加入到強地面運動模擬的結(jié)果中，使得模擬結(jié)果能夠充分反映淺層速度結(jié)構的影響。

2 強地面運動模擬及時程曲線的計算

衡水市位于河北省東南部，是重要的區(qū)域經(jīng)濟文化重鎮(zhèn)。衡水地區(qū)在大地構造上屬中朝準地臺的二級構造單元——華北斷坳區(qū)，它是我國最古老的陸塊之一，經(jīng)歷了太古代和早元古代基底形成、晚元古代至三疊紀蓋層發(fā)育和中、新生代以來地殼活化3個構造演化階段。在地殼活化階段，華北斷坳發(fā)生大規(guī)模的伸展裂陷作用。早第三紀發(fā)育了眾多的NE 向或NNE 向斷陷盆地，形成典型的盆嶺構造；晚第三紀華北斷坳大面積沉降，上第三系超覆于整個盆地；第四紀進一步發(fā)展，但構造分異比較顯著，形成了一系列?。ㄛ辏?、凸（凹）相間結(jié)構。平原地整體由東向西緩慢掀斜，內(nèi)部地勢平坦，一般海拔高100m 以下，大部分小于50m。衡水市位于華北平原拗陷區(qū)中的邢衡隆起、冀中坳陷和滄縣隆起交匯地帶，第四系厚度為350～400m。

根據(jù)衡水市城市發(fā)展和防震減災的需要，以及主要斷層的實際情況，確定本研究的目標區(qū)東西長30km，南北寬25km，面積約750km2。地理坐標為：北緯37.66°～37.81°，東經(jīng)115.52°～115.78°。為了考慮衡水市周邊斷層的影響，將工作區(qū)域擴大為北緯37.0°～38.5°，東經(jīng)114.5°～116.5°區(qū)域內(nèi)，包含（1）滄西斷裂，（2）舊城北斷裂，（3）護駕池斷裂，（4）新河斷裂，（5）前磨頭斷裂，（6）明化鎮(zhèn)斷裂，（7）衡水斷裂，（8）德州斷裂（圖1）。鑒于舊城北斷裂（2）和德州斷裂（8）只有極小一部分分布在近場區(qū)內(nèi)，我們只對其他6條斷層進行設定震級條件下的強地面運動評估。

考慮斷層活動性，我們將依據(jù)斷層長度得到各條斷層的最大震級與第5代潛在震源區(qū)劃圖進行對比。綜合評定后，得到活動斷裂潛在最大震級（表1）和隨機有限斷層模型輸入?yún)?shù)（表2）。

圖1 地震構造圖

3 結(jié)果分析與討論

本研究將研究區(qū)域劃定在（37.4°～38.0°N，115.2°～160.0°E），區(qū)域網(wǎng)格大小為0.6°×0.8°，網(wǎng)格精度為0.005°（0.5km），利用隨機振動有限斷層模型計算每個網(wǎng)格點的峰值加速度（PGA）、峰值速度（PGV）以及網(wǎng)格點到斷層的最短距離D。我們擬采用的計算流程是：

第1步：選定研究區(qū)域。根據(jù)地震地質(zhì)和地震活動性調(diào)查以及地球物理資料確定發(fā)震斷裂的位置、產(chǎn)狀、埋藏深度等性質(zhì)。根據(jù)震級大小以及震后可能造成的影響，選取合適的計算范圍。

第2步：利用隨機振動的有限斷層模型計算研究區(qū)域PGA、PGV 以及網(wǎng)格點到斷層的最短距離。根據(jù)第1步劃定的研究區(qū)域大小，我們會根據(jù)斷層的長度、寬度劃分合適的子斷層數(shù)目以及網(wǎng)格點的精度（為確保計算結(jié)果精確，設定計算網(wǎng)格點間距為0.005°），根據(jù)構造地質(zhì)結(jié)果以及歷史地震等資料確定發(fā)震斷層的應力降和幾何衰減因子等參數(shù)。

第3步：加入淺層速度結(jié)構的影響。在第2步計算得到初始地震動參數(shù)后，充分考慮淺層速度結(jié)構對地震動的影響。

表1 結(jié)合斷裂活動性估計潛在最大震級

表2 各設定發(fā)震斷層使用的隨機有限斷層模型輸入?yún)?shù)1 出自河北省活斷層探測與地震危險性評價項目衡水市目標區(qū)斷層危害性評價報告

計算出區(qū)域峰值加速度以及峰值速度分布后，我們計算不同等級的加速度、速度以及位移的時程曲線。通過輸入模型參數(shù)，按照上述步驟得到6個衡水設定地震的峰值加速度以及衡水市區(qū)的時程曲線（圖2～圖7）。

圖2 滄西斷裂設定6.0級發(fā)震斷層產(chǎn)生的峰值加速度及衡水市區(qū)時程曲線圖

圖3 護駕池斷裂設定6.0級發(fā)震斷層產(chǎn)生的峰值加速度及衡水市區(qū)時程曲線圖

圖4 新河斷裂設定7.5級發(fā)震斷層產(chǎn)生的峰值加速度及衡水市區(qū)時程曲線圖

圖5 前磨頭斷裂設定6.0級發(fā)震斷層產(chǎn)生的峰值加速度及衡水市區(qū)時程曲線圖

圖6 明化鎮(zhèn)斷裂設定6.0級發(fā)震斷層產(chǎn)生的峰值加速度及衡水市區(qū)時程曲線圖

圖7 衡水斷裂設定6.5級發(fā)震斷層產(chǎn)生的峰值加速度及衡水市區(qū)時程曲線圖

由于本研究設定的斷層埋深都比較淺，地震輻射能量集中于地表淺層。衡水市區(qū)與幾條斷層的距離都比較近，因此地震輻射能量衰減較少，在衡水市的強地震動效應也很明顯。由6條設定斷層強地震動計算結(jié)果看出，對衡水市區(qū)危害相對最小的是新河斷裂和護駕池斷裂，峰值加速度值在50～90cm／s2左右，峰值速度在10cm／s以下，這一值對應到GB／T 17742－2008《中國地震烈度表》為Ⅵ～Ⅶ度；對衡水主城區(qū)危害最大的是衡水斷裂，峰值加速度值在200～360cm／s2左右，峰值速度在30～45cm／s之間，這一值對應的烈度值為Ⅷ～Ⅸ度。

本研究采用的強地震動預測方法（以下簡稱活斷層方法）的原則是基于已經(jīng)查明的斷層，確定其可能發(fā)生的最大震級的地震前提下，對該斷層相應的設定震級，通過確定性的方法（確定性模型）進行地震動預測，并考慮其一旦發(fā)生破裂時產(chǎn)生的地震動空間分布，是在設定斷層的設定震級發(fā)生概率100%前提下得到的模擬數(shù)值，是針對性、確定性的結(jié)果。這種確定性的方法考慮了地震的震源特征、傳播途徑特性以及場地的放大效應，包含了強地震動模擬應當考慮的各種關聯(lián)特性。

模擬得到設定震級地震強地面運動的結(jié)果對于未來劃定該地區(qū)地震危險區(qū)具有很好的參考意義。由于該地區(qū)缺乏足夠的強地面運動記錄，無法利用實際記錄到的地震動資料直接得到地震動衰減關系，而模擬得到的強震數(shù)據(jù)可以補充衡水地區(qū)強震數(shù)據(jù)的不足，利用這些強震數(shù)據(jù)可以建立衡水地區(qū)的衰減關系曲線，這對于未來地震動預測、工程場地地震安全性評價和地震區(qū)劃工作將具有重要意義。由于隨機振動有限斷層模型沒有考慮地震波在介質(zhì)中的傳播過程，不會像復合震源模型那樣耗費大量時間在計算格林函數(shù)上，具備時效性強的特點，而這一點正是地震烈度速報所要求的。因此，未來該模型可以應用于地震烈度速報技術上，這對于震后應急響應將具有重要意義。綜合以上結(jié)果，本研究系統(tǒng)化應用有限斷層隨機振動模型，對包括衡水地區(qū)在內(nèi)的華北地區(qū)進行了地震風險評估及潛在震源的地震危險性分析，得到的基巖峰值加速度、峰值速度分布圖可為進一步劃定衡水地區(qū)地震危險區(qū)以及強地面運動特征分析提供參考依據(jù)。

致謝：感謝中國地震局地殼應力研究所馬保起研究員和吳榮輝研究員對本工作的大力支持和熱心指導。

［1］ 陳運泰，顧浩鼎.震源理論基礎［M］.北京：中國地震局地球物理研究所，2007：111－115.

［2］ 周錫元，吳育才.工程抗震的新發(fā)展［M］.北京：清華大學出版社，2006：18.

［3］ Boore D M.Simulation of ground motion using the stochastic method［J］.Pure and applied Geophysics，2003，160：635－676.

［4］ Boore D M.Stochastic simulation of high frequency ground motion based on seismological models of the radiated spectra［J］.Bull seism Soc Am，1983，73（6A）：1865－1894.

［5］ Brune J N.Tectonic stress and the spectra of seismic shear waves from earthquakes［J］.J Geophys Res，1970，75（26）：4997－5009.

［6］ Hanks T C，McGuire R K.The character of high frequency strong ground motion［J］.Bull Seism Soc Am，1981，71：2071－2095.

［7］ Atkinson G M.Attenuation of strong ground motion in Canada from a random vibration approach［J］.Bull Seism Soc Am，1984，74：2629－2653.

［8］ Beresnev I A，Atkinson G M.Stochastic finite－fault Modeling of ground motions from the 1994Northridge California Earthquake Validation on rock sites［J］.Bull Seism Soc Am，1997，88（6）：1394－1401.

［9］ 李啟成.經(jīng)驗格林函數(shù)方法模擬地震動研究［D］.哈爾濱：中國地震局工程力學研究所，2010.

［10］ 王國新，史家平.隨機有限斷層法在汶川強地震動模擬中的應用［J］.自然科學進展，2009，19（6）：664－669.

［11］ 王海云.近場強地震動預測的有限斷層震源模型［D］.哈爾濱：中國地震局工程力學研究所，2004.

［12］ 陶夏新，王國新.近場強地震動模擬中對破裂的方向性效應和上盤效應的表達［J］.地震學報，2003，25（2）：191－198.

［13］ Atkinson G M，Boore D M.Ground－motion relations for eastern North America［J］.Bull Seism Soc Am，1995，85（1）：17－30.

［14］ Chopra S，Kumar D，Rastogi B K.Estimation of Strong Ground Motions for 2001Bhuj（MW7.6），India Earthquake［J］.Pure Appl Geophys，2010，167：1317－1330.

［15］ Zafarani H，Noorzad A，Ansari A，et al.Stochastic modeling of Iranian earthquakes and estimation of ground motion for future earthquakes in Greater Tehran［J］.Soil Dynamics and Earthquake Engineering，2009：29：722－741.

［16］ Atkinson G M ，Silva W.Stochastic Modeling of California Ground Motions［J］.Bull Seism Soc Am，2000，90：255－274.

［17］ Motazedian D，Atkinson G M.Dynamic corner frequency：a new concept in stochastic finite fault modeling［C／OL］.Victoria，British Bolumbia：Seismological Society of America Meeting.［2008－12－25］.http：／／http－server.carleton.ca／～dariush／research／SSA2002.pdf

［18］ Wald D J，Allen T I.Topographic slope as a proxy for seismic site conditions and amplification［J］.Bull Seism Soc Am，2007，97（5）：1379－1395.

［19］ Boore D M，Atkinson G M.Boore－Atkinson NGA ground motion relations for the geometric mean horizontal component of peak and spectral ground motion parameters［R］.PEER Report，Pacific earthquake engineering research center，2007.