本章節闡述山腳斷層所彙整之震源參數,並以此作為依據,以地震學、地體動力學與工程地震學角度,考慮對北臺灣較具衝擊之情境,擬定單一特徵震源模型。另對於速度、地形構造之建構及場址放大效應函數進行說明。最終針對該單
一特徵震源模型進行情境模擬地震動,模擬完整情境地震歷時,計算地震動參數,包含(1) 峰值地動加速度(PGA)、(2) 峰值地動速度(PGV)、(3) 峰值地動位移(PGD)、(4) 0.3 秒譜加速度(SA0.3)與(5) 1.0 秒譜加速度(SA1.0)等五項參數作為
產出,以利後續災損推估小組推估北臺灣震災情境。綜整上述說明,模擬分析程序可詳見下圖。
本研究基於臺灣地震模型(Taiwan Earthquake Model, TEM)組織於 2015 發表之全台 38 條孕震構造,以及中央地質調查所 2012 年公布之臺灣活動斷層分布等資訊,以北臺灣之重要控制斷層-山腳斷層作為大規模地震情境模擬之假想對象。山腳斷層陸域部分,TEM 版本與地調所版本資訊無較大歧異,線型分布類似,山腳斷層為呈北北東走向之正斷層;此外,TEM 組織以高解析度地形資料審視構造線位置,海域延伸段採用國家地震工程中心 SHACC level 3 計畫海域地形資料檢視海陸域構造線連續性與延伸位置。綜合上述,本研究計畫彙整 TEM 版本與地調所版本山腳斷層資訊與線型作為情境模擬依據,並考慮山腳斷層於海域之延伸段(圖 2-2)
本研究之低頻地震波模擬方法,採用有限差分法(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) 地表地形,進行震波模擬。該方法採用 MacCormack 格式(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,低頻模擬方法無法反應 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)。
考量前述地栓滑移量分布之隱沒帶震源情境特徵震源模型,並根據「本研究設定點位」所述採行之 24,755 網格點位,以三維有限差分法與隨機式分析方法分別模擬各點位之低頻與高頻地震波,並以剪力波波速輔以互相關(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 子斷層與地表面 24,755 虛擬測站,其運算量共223,191,080 次點震源模擬,以 96 計算核心進行計算,低頻模擬約莫耗時 23 小時,高頻模耗時約 5.4 小時,混合法與場址修正耗時約 64 小時。後續若對單一斷層進行多次情境模擬,仍須審慎評估計算量與計算時間,合理使用並擴增計算資源使計算效益最大化。
考量前述地栓滑移量分布之隱沒帶震源情境特徵震源模型,並根據「本研究設定點位」所述採行之 24,755 網格點位,以三維有限差分法與隨機式分析方法分別模擬各點位之低頻與高頻地震波,並以剪力波波速輔以互相關(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 子斷層與地表面 24,755 虛擬測站,其運算量共223,191,080 次點震源模擬,以 96 計算核心進行計算,低頻模擬約莫耗時 23 小時,高頻模耗時約 5.4 小時,混合法與場址修正耗時約 64 小時。後續若對單一斷層進行多次情境模擬,仍須審慎評估計算量與計算時間,合理使用並擴增計算資源使計算效益最大化。
許多種數值方法陸續發展應用於三維地震波模擬,例如前述介紹的牽引力鏡像有限差分法(traction-image finite-difference method, Zhang et al.,2012),及其他常見方法如,交錯式網格有限差分法(staggered-grid finite-difference method, Open-source Seismic Wave Propagation Code, OpenSWPC, Maeda et al.,2017) 和譜元素法 (spectral-element method, Komatitsch and Vilotte, 1998 ;Komatitsch and Tromp, 1999)等方法。
關於不同數值方法之間的量化比較工作,第三屆的地表地質效應研討會(International symposium on the effects of surface geology, ESG 2006)已有團隊(Chaljub et al., 2010)進行詳盡的測試。該團隊比較四種不同數值方法,包括任意高階導數非連續伽遼金法(arbitrary high-order derivative discontinuous Galerkin
method, ADER-DGM, Käser et al., 2006) ,有限差分法(finite-difference method, FDM, Kristek et al., 2009),譜元素法(spectral-element method, SEM, Chaljub, 2009),以及第二種譜元素法(spectral-element method, SEM, Stupazzini, 2009 and Stupazzini et al., 2009)。結論認為,雖然現今數值方法仍在發展階段,但是在適當地使用下,四種方法皆可以正確處理複雜的三維速度模型異質性,獲得可靠的地動預估結果。此外,也建議採用兩種以上相對準確的方法進行模擬,例如有限差分法以及譜元素法,如此可以獲得更加準確可靠的結果。
本研究使用牽引力鏡像有限差分法(以下簡稱為 FDM1)、交錯式網格有限差分法(以下簡稱為 FDM2)以及譜元素法(以下簡稱為 SEM),進而比較三種數值方法採用相同參數時所模擬之合成波形,透過量化合成波形差異,理解不同的數值方法及不同的操作過程所可能出現的差異。比較過程中共測試兩組模型,第一組模型為速度均勻不具側向變化,且地形平坦之簡化模型,模型剪力波速為 2.0 km/s,壓縮波速為 3.5 km/s,密度為 2.34 g/cm3;採用點震源,震源走向 90°,傾角 90°,滑移角 180°,震源深度 8.0 公里位於模型原點(0,0);模型尺度範圍及測站配置下圖所示。第一組模型採用三種數值方法進行比較,三種方法之合成波形比較,以 LINE_001 測站陣列為例。 根據震源機制,理論上 X 與Z 分量振幅為 0,僅 Y 分量有值存在;三種方法之合成波皆呈現非常相近的結果,包含振幅與波形相位。
而為了量化合成波形之間的差異,本研究導入 Kristekova et al. (2009)所提出的時間頻率域擬合度 (Time-frequency misfit) 評估方法,以 TF_MISFIT_GOF_CRITERIA 套件,對合成波形之間進行差異的量化計算。Kristekova et al. (2009)對波形之間的差異在頻率域與時間域中進行量化,並給出量化分數介於 0~10 之間,當擬合度高於 8.5 時屬於”excellent”的擬合度,擬合度高於 6.5 時屬於”good”的擬合度,擬合度高於 4.5 時屬於”fair”的擬合度,擬合度低於 4.5 時屬於”poor”的擬合度。本研究皆以 FDM1 之合成波形做為參考波形,分別對 FDM2及 SEM 合成波形進行時間頻率域擬合度計算。FDM2 與 FDM1 之相位擬合度及振幅擬合度,其中相位擬合度只有在測站接近震源時,擬合度略為下降,或位於震源上方(Station ID=41)時之奇異點大幅下降,而其他測站之波形相位擬合度皆非常高;而振幅擬合度則呈現 Y 分量整體略為下降,然而平均擬合度仍在 8.5 分以上,X 與 Z 分量則皆趨近滿分。SEM 與 FDM1 之相位擬合度及振幅擬合度,其中相位擬合度只有在震源上方(Station ID=41)時之奇異點大幅下降,其餘測站皆趨近滿分擬合度;而振幅擬合度也呈現 Y 分量整體略為下降,平均擬合度高於 9 分,X 與 Z 分量則皆趨近滿分。
根據測試過程得知,FDM2 與 FDM1 之合成波形差異來源主要包含兩點:第一點為非交錯式網格有限差分法與交錯式網格有限差分法之數值方法本身對於座標系統之定義所帶來的誤差,即 FDM2 之實際測站位置與 FDM1 差距為半個網格點,以本測試來說大約差 12.5 公尺;第二點為震源時間函數定義不同,FDM1 採用 GAUSSIAN 函數,而 FDM2 採用 HEERMAN 函數做為震源時間函數。而 SEM 與FDM1 之合成波形差異比例極低,此差異來源屬不同的數值方法間的正常計算誤差;惟仍須注意 SEM 與 FDM 模型網格化時可能產生之座標誤差。
測試之第二組模型為考慮側向速度變化及地形起伏之模型,速度模型採用(Kuo-Chen et al., 2012)之三維速度分布數值,而地形起伏採用 ETOPO1 (Amante and Eakins, 2009)做為初始地形模型。震源則採用前期報告之山腳斷層南段破裂有限斷層模型。在第二組測試中僅以 FDM1 與 SEM 進行比較。選取 6 個測站與斷層的相對位置,而經選取 6 個測站之 FDM1 與 SEM 合成波形比較如下圖。合成波形比較顯示,在多數測站中,FDM1 與 SEM 之合成波形皆十分相似(圖 a、b、c、d);僅在近震源處(圖 e、f)有極少數測站在振幅上出現較明顯的差異。本研究在空間上均勻選取 1492 個測站進行時間頻率域擬合度計算,以量化評估 FDM1 與 SEM 之合成波形差異,振幅擬合度之分布,整體平均之擬合度達”excellent”,僅少數測點屬”good”。相位擬合度之分布,整體平均之擬合度均達”excellent”,波形相位擬和極佳。本研究研判,該較大波形差異測點的產生,可能與本研究半盲測試過程中,模型(速度,地形,座標)仍有些微不同所致。
經由兩組模型測試指出,有限差分法與譜元素法在適當的操作下,其計算結果非常接近,唯須注意操作過程的各項參數之設計,可能影響計算結果產生差異。