國立中央大學地震災害鏈風險評估及管理研究中心
大規模地震情境模擬對外資料
2018-隱沒帶大規模地震情境模擬
本研究計畫之工作項目，包含震源參數擬定、速度與構造模型建置、特徵震源模型建置、地震波模擬分析及地表地震動歷時與反應譜計算分析等工作。其中，參數與方法部分，模擬分析所需參數蒐集及擬定部分，陳述如何針對琉球隱沒帶地震蒐集相關文獻用以擬定震源模型參數；地形與速度構造模型建置部分，說明供後續情境地動模擬使用且覆蓋分析區域之地震波速度構造模型；特徵震源模型建置部分，考慮琉球隱沒帶地震之特徵規模及都會區可能之地動影響模式，設定可能震源破裂模式，詳加擬定震源參數，並製表呈現；地震波模擬分析章節闡述使用方法。另，後段則針對地表地震動歷時與反應譜計算分析，提出涵蓋研究區域內之模擬成果。
參數與方法
本處闡述考量琉球隱沒帶情境地震而彙整之震源參數，並以此作為依據，以地震學、地體動力學與工程地震學角度，考慮對北臺灣較具衝擊之情境， 擬定特徵震源模型。另對於速度、地形構造之建構進行考量。最終針對該單一特徵震源模型進行情境模擬地震動，模擬完整情境地震歷時，計算地震動參數， 包含(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移 (PGD)、(4) 0.3 秒譜加速度(SA0.3)與(5) 1.0 秒譜加速度(SA1.0)等五項參數作為產出，以利後續災損推估小組推估北臺灣震災情境。綜整上述說明，模擬分析程序詳見下圖。
震源參數擬定
琉球海溝具引致大規模隱沒帶地震之潛勢，歷史上於 1920 年，此區域卽發生規模 8 之地震事件，Theunissen et al. (2010)提出該區域潛勢斷層 ISZ:板塊介面孕震區域(interpolate seismogenic zone)模式為最可能發生之 1920 年地震與琉球隱沒帶地震之破裂型態。此外，Hsu et al. (2012)使用 2005 年至 2010 年間之 連續全球定位系统（Global Positioning System, GPS）資料，推估琉球隱沒帶斷層幾何及滑移虧損(slip-deficit)，該區域具有發生震矩規模 7.5 至 8.7 的災害地震潛勢。本研究綜合考量上述論文研究，並根據國家地震工程中心 SHACC level 3 計畫綜整近年臺灣海陸域斷層調查成果，作為擬定琉球隱沒帶情境震源參數之依據，詳見本報告文末「附錄：會議紀錄」。
地形與速度構造模型建置
地震波傳模擬所需速度構造模型建置，係以近期臺灣地區速度構造研究(Kuo-Chen et al., 2012)之三維速度分布數值，作為本研究模擬計算之模型建置基礎資料。該速度分布數據係藉由中央氣象局強震站及短週期地震站、中央氣象局及中央研究院寬頻地震站、大量布設臨時站所記錄到之大量地震波走時資訊，經由逆推方法計算建構而成。參考三維速度分布數值所強調之構造不均質(heterogeneity)與側向變化(lateral variation)，提升了僅考量無方向性變化的一維速度模型可能影響地震模擬評估之不足。評估本先導研究情境模擬工作所需要之涵蓋空間範圍後，整合三維速度分布數值所建構之北臺灣速度構造模型如下圖所示。
此外，因數值方法特性，本研究所採行之地震動模擬方法於低頻模擬時可考慮地形起伏造成的震波散射效應，例如震波傳遞至山峰及山脊處時所造成震幅放大，山谷處造成震幅衰減(Lee et al., 2009)。因地形效應額外產生之表面波與散射波，預期會影響測站處記錄之波形。本研究計畫以解析度為一弧秒(約1.85 km)之 ETOPO1 (Amante and Eakins, 2009)做為初始地形模型，考量前述北臺灣速度構造模型後，經穩定性測試後，計算出震波模擬適合之模型解析，將地形模型重採樣為 300 m 取一格點，作為數值計算網格水平方向格距，建立數值網格模型，模型空間範圍考量納入涵蓋琉球隱沒帶情境震源與北臺灣目標場域之範圍。
特徵震源模型建置
考量琉球隱沒帶線型有所轉折，完整詳實地震動模擬須考量具轉折之破裂面、細緻化震源模型。本研究簡化琉球隱沒帶線型，合理考量符合線型之斷層長度與寬度，進而假定特徵震源模型。
特徵震源模型之參數擬定原則，則採用日本地震調查研究推進本部之 Recipe (Irikura and Miyake, 2011)方法，取其特徵震源擬定之標準化流程作為依據，針對琉球隱沒帶情境地震破裂，探討長週期地震動對北臺灣影響性，假定斷層面上具單一地栓(asperity)為主要能量釋放位置，將地栓中心置於沿斷層走向之中線，地栓上緣與斷層上緣重合，並考慮震源(即破裂起始點)位於地栓上緣或下緣之四個不同位置，產製四組震源模型，意即評估不同斷層破裂方向性效應(directivity effect)對北臺灣之衝擊：(1) 模型 01(平面/側向)：震源位於地栓上緣最外海處，考量斷層由海床往深處破裂；(2) 模型 02(平面/側向)：震源位於地栓下緣最外海處，考量斷層由深處往海床破裂；(3) 模型 03(平面/側向)：震源位於地栓下緣中間位置，考量斷層由深處往海床破裂；(4) 模型 04(平面/側向)：震源位於地栓下緣最近海處，考量斷層由深處往海床破裂。四組琉球隱沒帶情境地震震源模型以平面與側向示意，詳見下圖；各組模型之震源位置(即破裂起始點)、斷層面之巨觀與微觀震源參數則分列於下表。
虛擬測站點位
NCDR 點位
本研究配合災損推估小組之評估網格，配合國家災害防救科技中心 (National Science and Technology Center for Disaster Reduction, NCDR)格點解析度為500公尺，分布範圍為東經118.125度至122.005度，北緯21.895度至26.385 度，範圍涵蓋臺灣本島與離島範圍，共 132,712 網格點。下圖為 NCDR 網格示意圖，黑色格點即為網格點位。
本研究設定點位
對應災損推估小組後續評估工作，空間範圍與解析上，本研究基於 NCDR 500 公尺網格點擬定虛擬測站，進行地震動模擬計算。所採行點位分布範圍為東經 121.000 度至 122.005 度，北緯 24.8 度至 25.3 度，共 13,216 點位，空間點位選取範圍系考慮琉球隱沒帶情境震源所在區域、北北基桃生活圈、北臺灣關鍵基礎設施進而擬定之範圍。如下圖所示，分布於北臺灣陸地區域之黑色格點為 NCDR 500 公尺網格點位，視為虛擬測站；紅框區域為情境地震動模擬範圍。
地震波模擬分析
低頻模擬：三維有限差分法
本研究之低頻地震波模擬方法，採用有限差分法 (finite-difference method)，於三維地形及地下速度構造模型設定中，透過擬定假想地震之情境(特徵震源模型)進行模擬，評估低頻地震波對特定空間點位、工址或測站產生之強地動影響。
震波傳遞之模擬方法係基於物理機制下，考慮地震源與震波傳遞路徑中各種耦合效應，計算特定點位的震動歷時反應，其中震源項包含地震震源機制、斷層幾何與破裂型態，而路徑項則包含不同尺度的地下構造模型(包含波速、密度與衰減因子等)。以有限差分法進行震波模擬，主要係由求解質點運動方程與 組構關係隨時間之變化，即為模擬質點受力後於空間中的運動行為與震波隨時空變化之傳遞過程。有限差分法為基於網格運算之數值方法，即為求網格系統中某格點位置之差分運算，運算過程僅需鄰近格點即可計算，因此可透過訊息傳輸介面(Message Passing Interface, Gropp et al., 1996)將數值方法改寫，並在叢集式電腦(PC cluster)進行平行化計算 (parallel computation)，進而提高效率或進行更為高頻的模擬工作。
本研究採用 Zhang et al. (2012)之牽引力鏡像有限差分法(traction-image finite-difference method)可考慮(1) 三維地下構造與(2) 地表地形，進行震波模擬。該方法採用 MacCormak 格式(MacCormack, 1969; Gottlieb & Turkel, 1976)與 optimized Dispersion Relation Preserving (DRP/opt, e.g. Tam & Webb, 1993; Hixon, 1997)能有效降低數值頻散(numerical dispersion)效應。同時，該演算法採取完美匹配吸收層(perfectly matched layer, PML)做為吸收邊界條件(absorbing boundary condition, e.g. Berenger, 1994; Marcinkovich & Olsen, 2003)，能有效降低人為設定模型邊界產生之反射波。
為考慮地表地形或地下構造反射面效應，Zhang et al. (2012)透過坐標轉換，將地表地形與地下構造界面以曲線性坐標系統(curvilinear coordinate system)描述，進而使物理空間中各曲線性坐標位置
透過 Jacobian 變換至電腦計算所使用的直角坐標系統 
。
高頻模擬：隨機式方法
針對需求所擬定強地動模擬，須計算 PGA、SA0.3 等高頻地震動參數，模擬頻段須達 10 Hz，2.6.1 節之低頻模擬方法無法反應 1-10 Hz 頻段之高頻地震動時間歷時，故本研究計畫選用 EXSIM (Atkinson and Assatourians, 2015)隨機式方法進行高頻地震動模擬，補足人工合成地震歷時之高頻響應。
EXSIM 係根據使用者給定之地震矩，以經驗公式給定斷層長寬形成一斷層面，而後離散化為一系列之子震源(subsource)，其中的每一個子震源皆遵循點震源隨機模型(point-source stochastic model)，依給定的地震歷時(duration)與剪力波輻射型態以隨機高斯噪訊(random Gaussian noise)模擬強地面加速度地動歷時。EXSIM 亦引入 Brune model，於傅氏譜(Fourier spectrum)考慮地震矩、應力參數、經驗衰減參數於頻率域之衰減效應，經快速傅立葉變換(fast Fourier transform, FFT)至時間域後，根據破裂速度疊加各子震源，並考慮場址效應修正函數，形成時間序列。本研究亦參考 Lee et al. (2015)根據場址 VS30 條件分類推估地震歷時之經驗公式，整合於高頻模擬計算。
EXSIM 為相當具代表性之隨機式強地動模擬方法，除南加洲地震中心(Southern California Earthquake Center, SCEC)已對 EXSIM 實行相關檢核(e.g.Atkinson and Assatourians, 2015)外，EXSIM 模擬結果亦已與隱沒帶地震事件檢核(e.g. Atkinson and Macias, 2009)。
本研究以 EXSIM 模組進行隨機法強地動模擬時，給定琉球隱沒帶特徵震源模型之滑移量分布作為輸入震源模型，並針對模擬範圍，考慮適用臺灣區域之隨機式模型(Sokolov et al., 2006; Sokolov et al., 2009; Huang et al., 2017)，詳細模擬參數列於下表。
混合法合成寬頻地震歷時
考量前述地栓滑移量分布之隱沒帶震源情境特徵震源模型，並根據採行之 13,216 網格點位，以三維有限差分法與隨機式分析方法分別模擬各點位之低頻與高頻地震波，並以剪力波波速輔以互相關(cross-correlation)方法修正各點位之震波到時(arrival-time)後疊加低頻與高頻地震波，以工程常用之 VS=760 m/s 基盤面做為地震動計算基準，以混合法合成寬頻地震歷時至該基盤面，後續便以該寬頻震動歷時進行 5 項地動參數計算，包含(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)、 (4) 0.3 秒譜加速度(SA0.3)與(5) 1.0 秒譜加速度(SA1.0)。
本研究以叢集式電腦進行平行運算，考慮單一情境，琉球隱沒帶情境地震之特徵震源震源模型共 9,016 子斷層與地表面 13,216 需擬測站，其運算量共119,155,456 次點震源模擬，以 120 計算核心進行計算，低頻模擬約莫耗時 18小時，高頻模耗時約 3 小時，混合法與場址修正耗時約 36 小時。後續若對單一斷層進行多次情境模擬，仍須審慎評估計算量與計算時間，合理使用並擴增計算資源使計算效益最大化。
地震動歷時與反應譜計算分析
依前述地震波模擬分析方法，考慮琉球隱沒帶情境地震之特徵震源模型，分別進行低頻與高頻模擬，進而以混合法合成於工程基盤之寬頻加速度地震歷時，爾後分別積分一次得速度歷時、積分二次得位移歷時，各地震歷時之最大震幅處即分別為(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)等地動數值。另一方面，選取加速度地震歷時作為輸入地震歷時，考慮 5% 阻尼計算該地震歷時之反應譜(response spectrum)，爾後分別取週期 0.3 秒與週期 1.0 秒之譜加速度作為反應譜地動數值。各組模型引致北臺灣地震動分述如下：
- 模型 01：考慮震源位於地栓上緣最外海處，斷層由海床往深處破裂，因斷層破裂的方向指向臺灣本島與北臺灣方向，破裂的方向性效應最為強烈，此起情境為四組震源模型中，引致北臺灣地震動最甚之情境，PGA、PGV、PGD、SA0.3 與 SA1.0 分布均較其他模型為高，PGA 分布顯示最大值為170.77 cm/s2，對應震度為 5 級，而北臺灣區域整體震度大約為 4 至 5 級；PGV 分布與 PGD 分布之最大值則分別為 59.97 cm/s 與 126.81 cm，量值較大反應長週期地震動的效應；SA0.3 分布與SA1.0 分布之最大值則分別為 0.33 g 與 0.31 g。
- 模型 02：考慮震源位於地栓下緣最外海處，斷層由深處往海床破裂，此情境事件破裂方向指向臺灣本島與南臺灣方向，地震動量值分布普遍次之於模型01，PGA 分布顯示最大值為 111.34 cm/s2，對應震度為 5 級，而北臺灣區域整體震度大約為 4 級；PGV 分布與 PGD 分布之最大值則分別為18.42 cm/s 與 30.53 cm；SA0.3 分布與 SA1.0分布之最大值則分別為 0.32 g 與 0.17 g。
- 模型 03：考慮震源位於地栓下緣中間位置，斷層由深處往海床破裂，此情境事件破裂方向為由斷層中央深部往海床，地震動量值普遍小於模型 01與模型 02，PGA 分布顯示最大值為 99.89 cm/s2，對應震度為五級，而北臺灣區域整體震度大約為 4 級；PGV 分布與 PGD 分布之最大值則分別為 12.27 cm/s 與 34.26 cm，；SA0.3 分布與 SA1.0 分布之最大值則分別為 0.29 g 與 0.15 g。
- 模型 04：考慮震源位於地栓下緣最近海處，斷層由深處往海床破裂，此情境事件破裂方向，地震動量值普遍小於其他三組模型，PGA 分布顯示最大值為 103.60 cm/s2，對應震度為 5 級，而北臺灣區域整體震度大約為 4 級；PGV 分布與 PGD 分布之最大值則分別為 12.01 cm/s 與44.02 cm；SA0.3 分布與 SA1.0 分布之最大值則分別為 0.25 g 與 0.14 g。
成果總結
過去，琉球隱沒帶地震事件曾引致地震災害，如 1920 年琉球隱沒帶地震事件、1986 花蓮外海地震事件等，接造成傷亡與建物毀損。本研究計畫擇定琉球隱沒帶作為假想目標，考慮單一地栓位置與不同震源起始破裂位置、震波傳遞等因素，擬定單一震源破裂情境，針對北臺灣地區以混合法進行地震動模擬，完成情境模擬地動參數，包含峰值地動加速度、峰值地動速度、峰值地動位移、0.3 與 1.0 秒譜加速度等，供後續災損推估小組進行評估。研究成果顯示，該震源模式對北臺灣區域造成長週期地震動，可能對高樓、長跨距橋梁、高架道路與軌道有較大影響。
然而，本研究僅考慮單一地栓震源破裂情境，倘若琉球隱沒帶破裂，不全然為本研究所假想之震源破裂型態以及推估規模；此外，場址效應對於臺北市區之地動放大結果影響甚巨。故本研究建議如下：
(1) 針對單一地栓震源破裂面，須考慮不同地栓位置之滑移量分布模式、震源破裂位置等多種破裂情境，對各情境產出之地震動參數進行統計分析，有效量化地震動中值與極值，估算地震動範圍。
(2) 詳細化大臺北地區之場址放大因子，空間高解析之準確場址放大對於提升地動空間解析有相當助益。
參考文獻
Amante, C. and Eakins, B. W. (2009). ETOPO1 1 arc-minute global relief model:
procedures, data sources and analysis, NOAA Technical Memorandum, NESDIS NGDC-24, 19 pp.
Atkinson, G. M and Assatourians, K. (2015). Implementation and validation of
EXSIM (a stochastic finite-fault ground-motion simulation algorithm) on the SCEC broadband platform, Seismol. Res.. Lett., 86(1), 48–60, DOI:10.1785/0220140097.
Atkinson, G. M. and Macias, M. (2009). Predicted ground motions for great
interface earthquakes in the Cascadia subduction zone, Bull. Seism. Soc. Am.,99(3), 1552–1578.
Berenger, J. A., (1994). Perfectly matched layer for the absorption of
electromagnetic waves, J. Comput. Geophys., 114, 185–200.
Gottlieb, D. and Turkel, E. (1976). Dissipative 2-4 methods for time-dependent
problems, Math. Comput., 30, 703–723.
Gropp, W., Lusk, E., Doss, N., and Skjellum, A. (1996). A high-performance, portable
implementation of the MPI message passing interface standard, Parallel Comput., 22(6), 789–828.
Hixon, R. (1997). On increasing the accuracy of MacCormack schemes for aero
acoustic applications, AIAA Paper, doi: 10.2514/6.1997–1586.
Hsu, Y. J., Ando, M., Yu, S. B., and Simons, M. (2012). The potential for a great
earthquake along the southernmost Ryukyu subduction zone, Geophys. Res.Lett., 39, L14302, doi:10.1029/2012GL052764.
Huang J. Y., Wen, K. L., Lin, C. M., Kuo, C. H., Chen, C. T., and Chang, S. C. (2017). Site
correction of a high-frequency strong-ground-motion simulation based on an empirical transfer function, J. Asian Earth Sci., 138, 399–415.
Irikura, K. and Miyake, H. (2011). Recipe for predicting strong ground motion from crustal earthquake scenarios, Pure Appl. Geophys., 168, 85–104, doi:10.1007/s00024-010-0150-9. Kuo-Chen, H., Wu, F. T., and Roecher, S. W. (2012). Three-dimensional P velocity
structures of the lithosphere beneath Taiwan from the analysis of Taiger and related seismic data sets, J. Geophys. Res., B6.
Lee, S. J., Komatitsch, D., Huang, B. S., and Tromp, J. (2009). Effects of topography
on seismic-wave propagation: an example from northern Taiwan, Bull. Seism. Soc. Am., 99, 314–325.
Lee, Y. T., Ma, K. F., Wang, Y. J., and Wen, K. L. (2015). An empirical equation of
effective shaking duration for moderate to large earthquakes, Nat. Hazards,75(2), 1779–1793.
MacCormack, R. W., (1969). The effect of viscosity in hypervelocity impact
cratering, AIAA Paper, 69–354.
Marcinkovich, C. and Olsen, K. (2003). On the implementation of perfectly matched
layers in a three-dimensional fourth-order velocity-stress finite difference scheme, J. Geophy. Res., 108, doi: 10.1029/2002JB002235. issn: 0148-0227.
Sokolov, V., Loh, C. H., and Jean, W. Y. (2006). Strong ground motion source scaling
and attenuation models for earthquakes located in different source zones in Taiwan, 4th International Conference on Earthquake Engineering, , Paper No.003.
Sokolov, V., Wen, K. L., Miksat, J., Wenzel, F., and Chen, C. T. (2009). Analysis of
Taipei basin response for earthquakes of various depths and locations using empirical data, Terr. Atmos. Ocean. Sci., 20, 687–702, doi:10.3319/TAO.2008.10.15.01(T).
Tam, C. and Webb, J. C. (1993). Dispersion-Relation-Preserving finite difference
schemes for computational acoustics, J. Comput. Phys., 107, 262–281.
Theunissen, T., Font, Y., Lallemand, S., and Liang, W. T. (2010). The largest
instrumentally recorded earthquake in Taiwan: Revised location and magnitude, and tectonic significance of the 1920 event, Geophys. J. Int., 183,1119–1133, doi:10.1111/j.1365-246X.2010.04813.x.
Zhang, W., Zhang, Z. G., and Chen, X. F. (2012). Three-dimensional elastic wave
numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids, Geophys. J. Int.,190, 358–378.
行政院災害防救辦公室(民 105)。105 年災害防救白皮書，行政院，191 頁，臺灣臺北。
陳亮全、劉楨業、陳宏宇、柯孝勳、張國鎮、許健智、林正洪、馬國鳳與劉兆漢(中央研究院大規模地震災害防治策略研議小組) (民 104)。大規模地震災害防治策略建議書，中央研究院，44 頁，臺灣臺北。
國家地震工程研究中心，臺灣地震危害高階模型建置，http://sshac.ncree.org.tw/
附錄
會議記錄
