不同數值模擬方法之比較

        許多種數值方法陸續發展應用於三維地震波模擬，例如前述介紹的牽引力鏡像有限差分法(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”，波形相位擬和極佳。本研究研判，該較大波形差異測點的產生，可能與本研究半盲測試過程中，模型(速度，地形，座標)仍有些微不同所致。

        經由兩組模型測試指出，有限差分法與譜元素法在適當的操作下，其計算結果非常接近，唯須注意操作過程的各項參數之設計，可能影響計算結果產生差異。