Excel PSHA: 基于 Excel 的概率式地震危害度分析 程式之開發及研究 王瑞斌 香港科技大學土木及環境工程學系 摘要 : 概率式地震危害度分析以地震發生之时空關係為基礎, 以概率的觀點評估研究區域在未來某一時間內可能發生之強地動值大小與地震災害發生幾率 在學界及工程實務中, 此分析多藉由计算机程式輔助完成 以美國地質調查局研发的计算机辅助分析软体为例, 傳統 PSHA 程式以 Fortran 語言撰寫居多 本研究開發出基於 Microsoft Excel 的地震危害度计算程式 Excel PSHA 我們選擇 Excel VBA 語言開發此技術程式之主要原因在於其作為微軟視窗內置程式所伴隨的廣大用戶群及友好的使用界面, 未來用戶僅需進一步將自定義 VBA 程式移植至 Excel PSHA 中即可進行個性化修改及分析計算 本項研究另一項重要技術則為以經由 Excel PSHA 程式計算及輸出的地震危害圖 在區域性災害分析中,Excel PSHA 可依據不同風險水平, 合理評估潛在地震災害之分佈 此外, 以瞭解場地之總地震危害主要來自于附近哪些主要震源為目的,Excel PSHA 成功整合了參數解構方法, 將超越某一地動值的年超越频率拆解成各規模區間與距離區間的貢獻情形, 并于 Excel 界面實時輸出三維圖表 此項研究成果為用戶提供直觀 易於使用的地震危害度分析工具, 藉由少量震源參數亦可作為工程實務中耐震設計之重要參考 關鍵字 : 概率式地震危害度分析,Excel, 地震災害
1. 引言 地震是人類所面臨的最大天然災害之一, 但迄今為止, 科學家仍無法有效地預測地震發生的時間 位置和其大小 因而在工程應用中, 普遍採用地震危害度分析對工程場地可能遭受到的地震影響程度做出科學評價, 並以概率式地震危害度分析及確定性地震危害度分析最具代表性 發展于 1960 年代的概率式地震危害度分析, 在近數十年內迅速獲得廣泛認同并于地震災害評估中普遍應用 [1-4], 特別地, 美國核管理委員會更以此種方法作為核電站耐震設計規範 [5] 確定性地震危害度分析作為計算機輔助分析, 現有的計算工具如 SEISRISK III 及 EZ-FRISK 皆由 FORTRAN 程式編譯 相較于基於 Excel VBA 之程式, FORTRAN 程式之用戶使用難度遠甚於前者, 而 Excel 更可通過內置工具箱實現實時結果及圖像顯示 誠然, 我們不否認 Excel VBA 計算效率之劣勢, 而若不執行大規模計算, 實際應用中其完成一次運算所需時間與 FORTRAN 無甚差異 藉于此, 在地球科學領域出現較多基於 Excel VBA 的程式以實現計算機輔助分析及計算 [6-9] 特別是 Excel 作為微軟視窗環境的捆綁應用程序, 其用戶可用性具有極大優勢 本研究開發出基於 Excel 的計算程式 Excel PSHA, 以實現地震危害度分析這一地球科學中日益重要的課題 本文涵蓋計算分析方法之闡述, 繼而深入探討 Excel PSHA 程式之開發方法 此外, 本程式可靠性亦由概率式地震危害度分析基準範例驗證 [10] 2. 概率式地震危害度分析之概述 概率式地震危害度分析考慮震源位置, 震源深度, 震源規模及強地動衰減式 之不確定性, 對場地做出基於概率的地震災害評估 其表達式如下 : S M D N S NM ND vi i 1 j 1 k 1 Y y* m, d Pr M m Pr D d ( y*, N, N, N ) Pr (1) j k j k 其中地震災害以年平均超越概率 λ 為表達, 作為地動值 (y*) 震源數目(N S ) 震源規模之概率密度分佈組段數量 (N M ) 震中距之概率密度分佈組段數量(N D ) 的函數 並且, 地震年平均概率 (v) 及 Pr[M = m j ] 由震源規模閥值 (m 0 ) 最大震源規模 (m max ), 以及地震再現參數 (a 值與 b 值 ) 控制, 地震災害由上述八個參數確定 而震源位置 震源區邊界 強地動衰減式 震源規模之概率密度分佈組距及震源區劃分亦決定概率密度方程 綜上所述, 概率式地震危害度分析由 13 個參數共同確定
3. Excel PSHA 電子錶格界面之闡述 圖一為概率式地震危害度分析基準範例三個震源區之建置 [10] 圖二為 Excel PSHA 程式輸入參數之屏幕視圖 用戶可以通過點擊程式中 Write Headings 按鍵, 便捷地在電子錶格相應位置填寫輸入參數 在輸入完成后, 用戶可以點擊 Excel PSHA 在 VBA 編程環境中運行此程序 圖三為程式輸出參數及實時災害評估曲綫之屏幕視圖 相較于 FORTRAN 語言開發之軟體, 此程式直觀的用戶介面和圖表顯示具有極大優勢 4. Excel PSHA 程式開發之闡述 由於地震發生之極大隨機性 [10], 此程式開發之首要步驟即為評估震中距之概率密度分佈 : 即通過對震源區的離散化得到一系列震中距, 並以統計學方法計算出相應之震中距概率密度分佈 此步驟需要一個判斷方程來實現, 用以確定空間中一點是否在給定震源區邊界內 圖四為此判斷方程的示意圖, 通過判斷交點 ( 原始點引出的射線與震源區邊界之交點 ) 之奇偶性, 進而判斷原始點是否在震源區域內 ( 奇數在內, 偶數在外 ) 由於面震源區于幾何上具有較大不確定性特徵, 此判斷方程須合理處理邊界條件及特殊情況, 以提高其可用性及可靠性 圖五為四種複雜面震源區情況下的驗證示意圖, 可以明顯看出, 本程序之判斷方程可以合理地離散震源區域, 并準確判斷原始點是否位於震源區內 代碼 此文對線形震源之離散化算法不作贅述, 具體細節可參考 Excel PSHA 程式 5. Excel PSHA 程式之驗證 本文通過一地震災害評估基準範例對 Excel PSHA 程式之可用性及可靠性做出驗證 由 Excel PSHA 程式計算出線形震源 面震源 點震源對於場地 (0,0) 構成的地震災害風險, 超越 0.01 g 地動值的年平均概率相應為 1.953,1.000 及 0.004, 與基準範例 [10] 提供的參考值 1.923,1.106 及 0.005 相近 本研究亦通過對更多不同場地的地震災害風險評估, 進而驗證 Excel PSHA 程序之可靠性 圖六顯示了 12 個檢驗場地 (S1 至 S12), 兩兩對稱于三個震源區域, 圖七為相應災害曲綫, 結果顯示, 對稱于震源區域的任意兩個場地災害曲綫相同 此驗證結果與預期設想一致, 因而 Excel PSHA 程式亦得到驗證
6. 參數解構 參數解構以瞭解場址的總地震危害度主要來自於附近哪些主要震源規模及位置為目的, 將超越某一地動值的年超越幾率拆解成各地震規模區間與距離區間的危害度貢獻情形 Excel PSHA 亦參考了美國核管理委員會技術規範, 以 0.0001 作為參數拆解中默認超越概率水平, 而用戶亦可以在電子錶格界面 Q5 位置上自行定義此數值 我們在 Excel PSHA 程式中整合了對 PSHA 結果的參數拆解分析 圖八為參數解構輸入項及輸出項的屏幕視圖 其中輸入項分別為兩行震級區間劃分及兩列距離區間劃分 ; 輸出項為兩個矩陣 : 用於生成三維參數拆解視圖的百分比矩陣, 以及用於後期邏輯樹分析的參數拆解災害矩陣 從基準範例 PSHA 結果可以看出, 在 0.0001 的超越概率水平下, 場址 (0,0) 近百分之四十的地震災害由震中距位於 25 公里到 50 公里範圍內震級高於七級的地震貢獻 Excel PSHA 程式最高可以實現對九個震級區間 十五個距離區間的參數拆解, 同時我們認識到由於地震數據庫的局限性, 過於細化的參數拆解不具有實際意義 7. 地震危害度圖 如前文所述,Excel PSHA 程式可以成功實現對給定場址的地震災害分析 災害曲綫輸出及參數拆解 在區域性地震災害分析中, 通常使用地震危害度圖將此地區地震危害度分析之成果予以清楚呈現 以美國地質調查局发布的美国地震危害度图为例, 分别反映了此区域内 50 年内超越概率为 2% 和 10% 水平下的地震灾害 此地震灾害与年超越概率 λ 相关, 并服从泊松概率分布模型 [10], 其表达式如下 : Pr( Y y*) 1 e y* t (2) 其中 Pr(Y > y*) 給定地動值 y* 的超越概率 ;t 為時間 ( 以年計 ) 整理表達 式 (2) 可得到年平均超越幾率為 : ln(1 Pr( Y y * t y*)) (3) 0.0021 由表達式 (3) 可以計算出 50 年內相應於超越概率 10% 的年平均超越幾率為
我們由此擴展了 Excel PSHA 之功能, 以實現區域性地震災害分析 圖九為 121 個場址建置以實施地震災害分析 圖十程式輸入項及輸出項的屏幕視圖, 類似于參數解構, 矩陣化的輸出項可以使 Excel 工具箱實現的三維視圖清晰可見 誠然,Excel 軟體中捆綁的等高線製圖功能不及 Origin 等商業軟件強大 因此, Excel PSHA 亦提供了數據結果輸出 ( 屏幕視圖中 D 列至 F 列 ), 為商業軟件後處理輸出高質量等高線製圖提供數據準備 圖十一顯示了此基準範例之 PSHA 計算結果, 并藉由 Origin 後處理的地震危害度圖 與參數解構類似,Excel PSHA 提供了地震災害結果的完整數據列 ( 屏幕視 圖中 G 列至 AH 列 ), 以待後續邏輯樹分析之用 8. 宏列表與功能闡述 由於 Excel PSHA 結構浩大, 功能繁多, 在其開發過程中, 作者以編寫子程序的方式實現各個功能模塊, 例如震源區離散化, 參數解構等等 具體言之, Excel PSHA 共計編寫了八百余行代碼以實現十五個 VBA 子程序模塊 在此, 我們不對子程序作逐行逐句贅述 9. 程序性能及其兼容性 儘管 Excel PSHA 之容量受限於其開發環境, 即 Excel 電子錶格之容量, 然而實際使用中此容量亦已充足 2003 版 Excel 電子錶格之行 列容量分別爲 65,536 行及 256 列, 因此 Excel PSHA 可以實現 100 余個震源區之分析 ( 每個震源佔用兩列電子錶格空間 ), 並且每個面震源區 線震源區可離散化為 65,000 余個微分區 此程式由 Excel 2003 開發, 經調試在 Excel 2007 及 Excel 2010 中亦可順利運行 同時, 伴隨更新版本之 Excel 電子錶格容量擴充,Excel PSHA 之容量亦顯著提高 在 Excel PSHA 中我們使用子程序模塊 PSHA_compute 編寫強地動衰減 式, 簡要示意如下 : 第一行 Sub PSHA_compute().. 第四十三行 gmm = 6.74 + 0.859 * m(k) - 1.8 * Log(d(j) + 25)..
第五十三行 End Sub 可以明顯看出, 用戶只需自行調整第四十三行便可實現對強地動衰減式之修改 與之類似, 用戶稍作修改 參數解構 子程序模塊中的強地動衰減式, 便可作自定義相應參數 由於上述基準範例使用笛卡爾座標系之建置, 若計算球面上兩點距離, 須修改 距離 子程序模塊并以大圓公式代之 在經緯座標系下, 用戶只需激活程序中的大圓公式以替代笛卡爾距離公式, 便可進行後續 PSHA 計算 近年來概率式地震危害度分析通常使用邏輯樹來考慮所需的參數及模式之不確定性, 如地震規模上限 地震規模分佈 強地動衰減公式等等 由於對於每個地震危害度分析, 其邏輯樹架構不盡相同, 因而, 我們建議用戶在不用工作簿中分別運行 Excel PSHA 程序, 最後進行整合以實現邏輯樹分析 最後, 我們在 Excel PSHA 電子錶格首頁提供了程式使用說明, 以供用戶參 考之用 附錄 : 專業名詞譯法說明 概率式地震危害度分析 : 確定性地震危害度分析 : 美國地質調查局 : 地震危害图 : 参数解构 : 地震年超越频率 : 规模区间 : 距离区间 : 耐震設計 : 面震源區 : 线形震源区 : 邏輯樹分析 : Probabilistic seismic hazard analysis Deterministic seismic hazard analysis United States Geological Survey Seismic hazard map Deaggregation Annual exceedance rate Magnitude bin Distance bin Earthquake-resistant design Area source Line source Logic tree analysis
Y-Coordinate (km) 80 60 40 (-15, 75) Point Source at (0, 60) a = 2.7 b = 1.2 m 0 = 4.0 m max = 5.0 (20, 78) Area Source a = 3.5 b = 0.8 m 0 = 4.0 m max = 7.7 (100, 78) 20 0-20 Line Source a = 4.4 b = 1.0 m 0 = 4.0 m max = 7.3 (20, 18) Site at (0, 0) (100, 18) (-15, -30) -40-60 -40-20 0 20 40 60 80 100 X-Coordinate (km) 圖一 概率式地震危害度分析基準範例三個震源區之建置 圖二 Excel PSHA 程式參數輸入界面之屏幕視圖
圖三 Excel PSHA 參數輸出及實時災害評估曲綫之屏幕視圖 Outside of the Polygon Crossing Points = 2 Inside of the Polygon Crossing Points = 1 Origin Outside of the Polygon Crossing Points = 0 圖四 判斷方程示意圖, 透過交點個數奇偶性判斷原始點是否位於震源區域內
Y-Coordinate (km) Y-Coordinate Y-Coordinate 18 (a) 18 (b) 16 16 14 14 12 12 10 10 8 10 12 14 16 8 10 12 14 16 18 (c) 18 (d) 16 16 14 14 12 10 12 8 6 10 4 8 10 12 14 16 X-Coordinate 8 10 12 14 16 X-Coordinate 圖五 Excel PSHA 程序驗證示意圖 80 S3 S11 S8 S6 60 S9 S10 40 20 0 S2 S1 S12 S5 S7-20 S4-40 -60-40 -20 0 20 40 60 80 100 X-Coordinate (km) 圖六 震源區域與 12 個檢驗場地建置之示意圖
Annual Rate of Exceedance Annual Rate of Exceedance Annual Rate of Exceedance (a) Test Points for the Line Source 1x10 0 1x10-1 S1 (-20, 20.67) S2 (-45, 18.33) S3 (-15, 30) S4 (-50, 75) 1x10-2 1x10-3 1x10-4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Peak Ground Acceleration (g) (b) Test Points for the Area Source 1x10 0 1x10-1 S5 (20, 18) S6 (100, 78) S7 (60, 10) S8 (60, 86) 1x10-2 1x10-3 1x10-4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Peak Ground Acceleration (g) 1x10-2 (c) Test Points for the Point Source 1x10-3 S9 (-10, 60) S10 (10, 60) S11 (0, 80) S12 (0, 40) 1x10-4 1x10-5 0.0 0.2 0.4 Peak Ground Acceleration (g) 圖七 12 個檢驗場址之地震災害曲綫,(a) 線形震源,(b) 面震源,(c) 點震 源
圖八 參數解構輸入項及輸出項之屏幕視圖 : 紅色為輸入區域, 藍色為 輸出區域
Y-Coordinate (km) 150 125 100 75 50 Point Source 25 0 Area Source -25-50 Line Source -75-100 -100-75 -50-25 0 25 50 75 100 125 150 X-Coordinate (km) 圖九 121 個場址建置, 用以實施地震災害分析并繪製災害圖
圖十 地震災害圖模塊輸入項及輸出項屏幕視圖 : 紅色為輸入區域, 藍色為輸 出區域, 等高線地震災害圖亦實時輸出
Y Coordinates (km) Y Coordinates (km) 150 100 50 0-50 (a) 2% Eexceedance 0.10 Pr. in 50 yrs 0.50 0.30 0.50 0.60 0.20 0.70 0.40 0.10 Source Boundary Source Location 0 0.10 0.10-100 -100-50 0 50 100 150 X Coordinates (km) PGA in g 0.80 0.70 0.60 0.50 0.40 0.30 0.20 150 (b) 2% Eexceedance 0.16 Pr. in 50 yrs 0.16 1.3 100 50 1.1 0.81 0.49 1.1 0.98 0.81 0.65 0-50 0.16 0.33 0.81 0.98 0.65 Source Boundary 0.16 Source Location 0-100 -100-50 0 50 100 150 X Coordinates (km) PGA in g 0.49 0.33 圖十一基準範例之地震災害圖 :(a) 地動值於五十年內超越概率為 10%; (b) 地動值於五十年內超越概率為 2%
參考文獻 1 Cheng, C.T., Chiou, S.J., Lee, C.T., Tsai, Y.B., 2007. Study on probabilistic seismic hazard maps of Taiwan after Chi-Chi earthquake. Journal of GeoEngineering 2, 19-28. 2 Mualchin, L., 2005. Seismic hazard analysis for critical infrastructures in California. Engineering Geology 79, 177-184. 3 Sokolov, V.Y., Wenzel, F., Mohindra, R., 2009. Probabilistic seismic hazard assessment for Romania and sensitivity analysis: A case of joint consideration of intermediate-depth (Vrancea) and shallow (crustal) seismicity. Soil Dynamics and Earthquake Engineering 29, 364-381. 4 Wiemer, S., Giardini, D., Fah, D., Deichmann, N., Sellami, S., 2009. Probabilistic seismic hazard assessment of Switzerland: best estimates and uncertainties. Journal of Seismology 13 (4), 449-478. 5 U.S. Nuclear Regulatory Commission, 2007. A performance-based approach to define the site-specific earthquake ground motion. Regulatory Guide 1.208. 6 Brock, B.W., Arnold, N.S., 2000. A spreadsheet-based (Microsoft Excel) point surface energy balance model for glacier and snow melt studies. Earth Surface Processes and Landforms 25 (6), 649-658. 7 Ozcep, F., 2010. Soil Engineering: a Microsoft Excel spreadsheet program for geotechnical and geophysical analysis of soils. Computers & Geosciences 36, 1355-1361. 8 Ozkaya, S.I., Mattner, J., 1996. Drag An excel visual basic program for modeling fault drag using cubic splines and calculation of minimum dip and strike separation. Computers & Geosciences 22 (5), 517-524. 9 Wang, J.P., Huang, D.R., 2012. Rosenpoint: A Microsoft Excel-Based Program for the Rosenblueth Point Estimate Method and an Application in Slope Stability Analysis. Computers & Geosciences, DOI: 10.1016/j.cageo.2012.01.009. 10 Kramer, S.L., 1996. Geotechnical Earthquake Engineering. Prentice Hall Inc., New Jersey, 116-117.