李得博 李井岡 趙艷南

1 中國地震局地震研究所,武漢市洪山側(cè)路40號,430071 2 中國地震局地震大地測量重點實驗室,武漢市洪山側(cè)路40號,430071

地震危險性分析是對工程場地未來可能遭遇的強地震動進行估計,是在工程場址周圍地震活動性和地震地質(zhì)環(huán)境研究的基礎(chǔ)上估計工程場址的設(shè)計地震動參數(shù)。評定某個工程場址的地震危險性須考慮3個方面:1)確定該場址周圍未來一定時間內(nèi)可能發(fā)生破壞性地震的地點、強度和性質(zhì);2)通過強震觀測數(shù)據(jù)和歷史地震資料,了解地震波傳播途徑的區(qū)域特性,確定地震動衰減關(guān)系;3)詳細勘察場地條件,分析局部場地條件對地震動的影響[1]。目前地震危險性分析方法大體上分為確定性方法和概率方法[2]兩大類,其中概率地震危險性分析方法由章在墉等[3]率先引入國內(nèi),經(jīng)高孟潭[4]和時振梁等[5]改進。改進后的方法能夠反映我國地震活動的時空不均勻性,被廣泛應用于地震安全性評價和地震小區(qū)劃等工作中。

地震災害普查中要求按國家級30″×30″網(wǎng)格、省級6″×6″網(wǎng)格、市級3″×3″網(wǎng)格、縣級1″×1″網(wǎng)格劃分經(jīng)緯度計算地震危險性,計算控制點數(shù)量達到千萬量級。若每個計算點耗時0.5 s,則省級計算點約1 000萬個,耗時約1 300 h,因此優(yōu)化相關(guān)計算方法是一個重要的研究課題。大空間網(wǎng)格的計算結(jié)果顯示,各相鄰點之間的危險性計算結(jié)果變化較為平穩(wěn),若能夠通過大間隔網(wǎng)格控制點的危險性結(jié)果插值獲得可靠的小空間間隔控制點的地震危險性結(jié)果,將大大節(jié)省計算時間?；诖?本文先計算120″、60″、30″、15″、9″和6″經(jīng)緯度劃分網(wǎng)格的地震危險性,然后插值得到較小網(wǎng)格的地震危險性結(jié)果,最后比較直接計算的小網(wǎng)格地震危險性結(jié)果和大空間網(wǎng)格插值得到的小網(wǎng)格地震危險性結(jié)果,評估由較大空間間隔插值所獲結(jié)果的可靠性、適用性及使用要求,以討論地震災害風險普查中兼顧可靠性和效率的計算地震危險性的最優(yōu)插值網(wǎng)格。

1 概率地震危險性分析方法

概率地震危險性分析方法(PSHA)于20世紀80年代引入我國,之后逐漸形成了具有中國特色的考慮地震時空分布不均勻性的概率性地震危險性分析方法(CPSHA),GB 17741-2005《工程場地地震安全性評價》和GB 18036-2015《中國地震動參數(shù)區(qū)劃圖》(簡稱五代區(qū)劃圖)都采用CPSHA方法,其計算步驟如圖1所示,包含確定地震區(qū)帶及其地震活動性參數(shù)、選定地震動衰減關(guān)系和進行地震危險性計算等主要步驟。

圖1 概率地震危險性分析步驟Fig.1 Steps of probabilistic seismic hazard analysis

根據(jù)地震風險普查規(guī)范,區(qū)域的標準網(wǎng)格分為國家、省、市、縣4種,國家級網(wǎng)格為30″×30″(約1 km),省級網(wǎng)格為6″×6″(約200 m),市級網(wǎng)格為3″×3″(約100 m),縣級網(wǎng)格為1″×1″(約30 m)。其中30″標準網(wǎng)格文件可通過軟件平臺下載,6″標準網(wǎng)格嚴格按照標準制作,地震、地質(zhì)、氣象等所有災種統(tǒng)一使用。標準網(wǎng)格基準點(西北角)為73°E、54°N,行列號向東及向南逐漸增加;標準網(wǎng)格制作中,行列號根據(jù)給定公式進行計算,作為對計算點的標簽,使用網(wǎng)格點中心點坐標(xcenter,ycenter)計算地震危險性;標準網(wǎng)格要求經(jīng)緯度保留小數(shù)點后8位,空間位置與屬性信息字段的經(jīng)緯度位置保持一致;分別計算50 a 63%、50 a 10%、50 a 2%和100 a 1%共4個概率水準的峰值加速度(PGA),不計算反應譜,給出基巖和場地調(diào)整的PGA,其中場地條件采用全國1∶100萬宏觀場地類別數(shù)據(jù),場地調(diào)整系數(shù)參考《地震危險性圖編制規(guī)范FXPC/DZ P-01》,場地類別根據(jù)標準網(wǎng)格中心點進行判別,地震危險性分級方法按照《地震危險性圖編制規(guī)范FXPC/DZ P-01》以100 a 1%概率結(jié)果作為分級指標(不使用50 a 10%),分級排序采用降序(危險性最高為1級,最低為4級)。該策略的目的一方面是為了避免使用50 a 10%造成混淆誤解;另一方面,萬年一遇可反映計算點可能遭受的最大地震動。潛在震源區(qū)、地震活動性參數(shù)和衰減關(guān)系使用五代區(qū)劃圖資料,不考慮近幾年的潛源調(diào)整。

因本文主要關(guān)注計算效率和插值方法的有效性問題,故區(qū)域范圍統(tǒng)一采用計算點周圍150 km范圍截取潛在震源的方法確定,暫不考慮150 km范圍外遠大震的影響,PGA計算結(jié)果統(tǒng)一保留小數(shù)點后1位。

2 測試區(qū)域選取

測試區(qū)范圍為34.55°～35.88°N、114.80°～116.42°E(圖2),主要基于以下考慮:1)根據(jù)五代區(qū)劃圖,區(qū)域場地地震動分類較全,0.05～0.20g都有,便于比較不同等級分解處插值是否合理;2)地震安全性評價報告資料豐富,便于對比;3)區(qū)域大小合適,采用不同網(wǎng)格間隔進行劃分,網(wǎng)格點數(shù)跨越幾千到幾百萬量級,利于全面研究不同量級的計算效率和精度。

圖2 測試區(qū)位置及120″間隔網(wǎng)格點分布Fig.2 Location of test area and distribution of 120″ interval grid

計算區(qū)間及網(wǎng)格劃分參數(shù)為:全國網(wǎng)格劃分統(tǒng)一起點緯度為54°,起點經(jīng)度為73°,網(wǎng)格間隔δ為3″ ,分別使用120″、60″、30″、15″、9″、6″及3″進行網(wǎng)格劃分;行號起點為10 001+int[((ycenter-54)×3 600-δ/2)/δ+0.5],列號起點為10 001+int[((xcenter-73)×3 600-δ/2)/δ+0.5],式中,int表示取整數(shù)。按照規(guī)范,不同級別的網(wǎng)格行列號的起點數(shù)不同,本文不作區(qū)分,只要能夠通過行列號找到對應的點即可。

單個控制點地震危險性分析計算輸出結(jié)果保存在以“點序號-列號-行號”命名的單獨文件夾中,其中包括5個周期(0 s,0.2 s,1.0 s,2.0 s,6.0 s)的危險性結(jié)果及相應潛源的貢獻、150 km范圍內(nèi)五代區(qū)劃圖中潛源分布(圖3)、各潛源的空間分布函數(shù)和風險普查要求的4個概率水準PGA。前期不同周期的危險性結(jié)果和潛源貢獻分別存放在不同的文件中;后期為減少文件讀寫數(shù)目,盡可能提高計算效率,將危險性結(jié)果和潛源貢獻結(jié)果合并存儲在同一個文件中。實際計算表明,該措施對計算速度的影響不大,但可減少小文件的數(shù)目,提高計算機存儲空間的利用率。風險普查要求計算4個概率水準(50 a 63%、50 a 10%、50 a 2%和100 a 1%)下5個周期的譜加速度值,由于數(shù)據(jù)量巨大,且不同周期、不同場地條件下結(jié)果類似,本文僅展示基巖場地50 a超越概率10%的PGA比較結(jié)果。

圖3 測試區(qū)內(nèi)某個計算點周邊150 km范圍內(nèi)五代區(qū)劃圖潛源分布Fig.3 Distribution of potential sources within 150 km of a calculation point in the test area in five generations zoning map

3 提高計算效率的方法及相關(guān)結(jié)果

將測試區(qū)劃分出120″、60″、30″、15″、9″、6″和3″間隔的網(wǎng)格,計算所有網(wǎng)格點各種概率水準下的PGA。如表1所示,不同間隔劃分的總網(wǎng)格點數(shù)隨著劃分間隔的減小呈指數(shù)增長,相應的計算耗時也大幅增長。首先分別計算不同空間間隔的各網(wǎng)格點的地震危險性,然后根據(jù)較大空間間隔的計算結(jié)果,按照空間距離平方倒數(shù)加權(quán)的插值方式得到所有較小空間間隔點的PGA,最后將插值所得PGA與直接計算的網(wǎng)格點PGA(準確值)進行比較,計算誤差最大值、最小值及標準差(表2,單位Gal)。

表1 不同空間間隔的計算量和耗時Tab.1 The computation amount and time of different space intervals

表2 不同大網(wǎng)格插值結(jié)果的誤差最大值、最小值、標準差Tab.2 The maximum, minimum and standard deviation error of interpolation results with different larger grids

表1給出不同處理器的計算耗時,其中,CPU為i7-11700的臺式機計算了全部數(shù)據(jù),其他幾種處理器只處理了15″以上間隔的數(shù)據(jù),更小間隔的數(shù)據(jù)是根據(jù)單點速度乘以計算量估算出來的。

由于每個點地震危險性計算相互獨立,不存在數(shù)據(jù)交叉,可以完全并行計算而不沖突。程序采用多線程,總體來看,每個程序的線程數(shù)對速度影響并不大。單個程序運行時,CPU為i7-9850H的移動工作站計算速度最慢,計算效率約5 000點/h,CPU為i7-11700的臺式機計算速度最快,計算效率約10 000點/h。為進一步提高計算速度,測試同一臺計算機采用2、4、8個相同程序同時計算不同區(qū)域的方法。結(jié)果表明,計算速度總和有所提高,但無論采用多少個應用程序,單臺計算機總最高計算速度約為18 000點/h,這個速度上限可能是由計算機的CPU處理能力決定的。

從表2和部分大網(wǎng)格插值得到的小網(wǎng)格結(jié)果的誤差分布(圖4～7)可以看出:1)利用120″網(wǎng)格插值結(jié)果的誤差最大值都達到2.0 Gal,這可能是由大網(wǎng)格相鄰2點危險性計算值的差值決定的;60″網(wǎng)格插值結(jié)果的最大誤差為1.9 Gal;30″網(wǎng)格插值結(jié)果的最大誤差為0.7 Gal;15″ 網(wǎng)格插值結(jié)果的最大誤差為0.4 Gal;9″ 網(wǎng)格插值結(jié)果的最大誤差為0.3 Gal;6″ 網(wǎng)格插值結(jié)果的最大誤差為0.2 Gal。2)所有大網(wǎng)格插值結(jié)果的誤差的標準差都較小,說明絕大多數(shù)點的誤差較小,插值結(jié)果誤差分布也說明了這一點。3)對于國家級30″ 網(wǎng)格,采用60″ 網(wǎng)格的計算值進行插值,誤差小于1 Gal;對于省級6″ 網(wǎng)格,采用30″、15″ 和9″ 網(wǎng)格的計算值進行插值,誤差小于1 Gal;對于市級3″ 網(wǎng)格,采用30″、15″、9″ 和6″ 網(wǎng)格的計算值進行插值,誤差小于1 Gal;同理,對于縣級1″網(wǎng)格,有理由推測采用30″、15″、9″、6″和3″ 網(wǎng)格的計算值進行插值,誤差也小于1 Gal。

圖4 不同大網(wǎng)格插值得到15″網(wǎng)格PGA的誤差分布Fig.4 The error distribution of 15″ grid PGA interpolated from different larger grids

圖5 不同大網(wǎng)格插值得到9″網(wǎng)格PGA的誤差分布Fig.5 The error distribution of 9″ grid PGA interpolated from different larger grids

圖6 不同大網(wǎng)格插值得到6″網(wǎng)格PGA的誤差分布Fig.6 The error distribution of 6″ grid PGA interpolated from different larger grids

圖7 不同大網(wǎng)格插值得到3″網(wǎng)格PGA的誤差分布Fig.7 The error distribution of 3″ grid PGA interpolated from different larger grids

為進一步研究插值誤差產(chǎn)生的原因,將誤差較大的點投影到地圖上。圖8給出由60″網(wǎng)格地震危險性計算值插值得到3″網(wǎng)格地震危險性結(jié)果時誤差大于0.5 Gal和大于1.0 Gal點的分布,其他大網(wǎng)格計算值插值得到較小網(wǎng)格結(jié)果的情況與之類似。圖8(a)顯示,誤差較大的點主要位于分區(qū)邊界處、較復雜處及計算區(qū)域邊界處;圖8(b)顯示,邊界處誤差較大,這主要是由于待插值的點位于計算值網(wǎng)格的外側(cè),而插值約束點非常少,此時采用空間距離反距加權(quán)方法不能反映變化的趨勢性?？梢酝ㄟ^擴大計算區(qū)范圍,將所有需要插值的點都作為內(nèi)點進行插值,或?qū)τ谟嬎銋^(qū)域外邊界點采用趨勢外推插值的方法進行插值。由此可見,目前地震風險普查實踐中一般要求計算點取行政邊界外圍擴大10 km的標準格網(wǎng),以減少邊界效應對插值結(jié)果影響的做法是合適的。

圖8 由60″網(wǎng)格插值得到3″網(wǎng)格的PGA誤差分布Fig.8 The PGA error distribution of 3″ grid interpolated from 60″ grid

4 結(jié) 語

本文基于五代區(qū)劃圖的潛在震源區(qū)劃分方案,采用全國地震災害風險普查統(tǒng)一的危險性計算方法,分別獲得了120″、60″、30″、15″、9″ 和6″ 經(jīng)緯度網(wǎng)格插值得到較小網(wǎng)格的地震危險性結(jié)果,并將其與直接計算的小網(wǎng)格地震危險性結(jié)果進行比較。結(jié)果表明,30″ 網(wǎng)格計算值插值得到各種小網(wǎng)格的地震危險性結(jié)果誤差都小于1 Gal,是兼顧精度和速度比較合理的網(wǎng)格間隔。

由于誤差較大的點主要位于地震動區(qū)劃圖分區(qū)邊界處、較復雜處及計算區(qū)域邊界處,通過擴大計算區(qū)范圍,將所有需要插值的點都作為內(nèi)點進行插值,可以解決計算區(qū)域邊界處誤差較大的問題。當然,在不需要等待很長時間的情況下,還是選擇盡可能小的網(wǎng)格進行計算,特別是潛在震源區(qū)比較復雜的區(qū)域。

在實際地震災害風險普查的危險性計算中,根據(jù)精度要求和時間因素進行一些前期計算實驗,綜合評估精度和效率,再確定計算方案是很有必要的。值得注意的是,本文主要是研究地震災害危險性普查中地震危險性計算的效率問題,對地震危險性計算中的地震潛在震源、地震動衰減關(guān)系分析等都進行了簡化,因此本文結(jié)果不能直接應用于測試區(qū)的地震安全性評價。