地震災害風險普查中的危險性計算與數據處理

邵 磊，胡 剛

（1. 湖南省地震局, 長沙 410004；2. 湖南省震災風險防治中心, 長沙 410001；3. 中國地震局地球物理研究所, 北京 100081）

0 引言

地震危險性分析是地震災害風險普查工作中的重要環節，可為地震災害風險評估與區劃提供危險性輸入[1]。根據地震危險性圖編制與技術規范的要求，采用《中國地震動參數區劃圖：GB18306—2015》使用的潛源劃分方案、地震動衰減關系、地震區帶劃分方案及地震活動性等參數，進行地震危險性概率計算[2]，并根據1∶100 萬宏觀場地類別進行場地調整、建立基巖和場地峰值加速度矢量數據庫[3]。此次工作任務有國家級、省級、市級以及區縣級，按照要求需建立國家級30″、省級30″或6″的標準格網；市級和區縣級因地區而異，有些要求3″、甚至1″格網。格網精度的不同往往對應不同級別的數據量，在計算與后續數據處理上可能會產生一系列問題，如計算點經緯度坐標是否為格網中心位置、地震危險性概率分析是選擇逐點計算還是通過插值實現、矢量數據shp 文件的制作、矢量數據間拓撲檢查是否滿足要求等問題。與常規的地震安全性評價中地震危險性分析工作相比，本次風險普查工作的數據計算與處理工作量巨大，往往都是十萬、百萬甚至千萬級別，且需要按照要求建立矢量數據庫，需要熟悉相關軟件和一定的編程處理數據經驗才可得到滿足要求的結果。本文沒有涉及地震危險性概率分析方法的原理，而是根據筆者在從事風險普查地震危險性計算中積累的經驗、做法以及可能出現的問題，運用Matlab、Python 代碼結合ArcGIS 等軟件平臺，主要從數據處理角度進行梳理，可省時、高效的完成該部分工作。在后續地震災害風險普查成為常態化工作中，可為繼續從事該方向的科研工作者提供一些研究思路，部分做法在地震安全性評價、地震區劃工作中也有一定參考價值。

1 標準格網的建立

根據地震災害風險普查工作的要求，一般是以某個省、市、區縣的普查行政邊界為工作目標區（以下簡稱目標區），生成所要求精度的標準格網，以各格網中心經緯度位置為計算點，進行地震危險性計算與分析。推薦使用Matlab 程序編寫較短的代碼即可實現標準格網的自動生成與計算點經緯度坐標的輸出，主要使用函數shaperead 讀取目標區邊界數據矢量特征和屬性，結合函數inpolygon 加以判斷控制點是否在目標區內，同時使用函數shapewrite 可輸出所需目標區格網矢量數據shp 文件[4]。使用程序代碼操作十分便捷、運行速度快、可重復性操作效率高。為滿足質檢要求，需按照風險普查所要求的某個省份格網范圍開始起算，程序在循環上設置滿足所要求格網大小的一定倍數即可快速判斷?？紤]到與邊界相交的格網中心點在邊界內或外的問題是否滿足數據質檢規則，此處建議按照比目標區行政邊界稍大的矩形區域或對目標區行政邊界生成一定距離的緩沖區范圍所有計算點，進行整體計算與數據處理，最后剪切到目標區范圍所需的上傳至普查系統平臺的矢量文件。在生成標準格網的同時輸出格網中心點坐標文件，坐標經緯度小數點后建議不小于10 位，便于后續處理滿足質檢不少于8 位且不四舍五入的數據要求。主要使用的matlab 代碼如下。

NDP = shaperead(‘矢量邊界shp 文件’);

BoundaryLon=[NDP.X];BoundaryLat=[NDP.Y];

inpolygon(Longitude,Latitude,BoundaryLon,Bound aryLat);

P(i).Geometry = 'Point';

P(i).X = Longitude;

P(i).Y = Latitude;

P(i).xcenter = Longitude;

P(i).ycenter = Latitude;

shapewrite(P,‘輸出shp 文件’).

2 地震危險性計算與數據處理

2.1 地震危險性概率計算

利用ArcGIS 軟件添加格網中心點坐標或程序代碼直接輸出的計算點坐標數據，使用《中國地震動參數區劃圖：GB18306—2015》的潛源等數據資料[5]，進行4 個超越概率（50年超越概率63%、10%、2%和年超越概率10-4）的基巖地震動峰值加速度的計算。由于計算點數量較大，需考慮計算效率與計算機存儲問題，可從幾個方面考慮：①由于較遠的潛源對基巖峰值加速度的貢獻很小或可忽略不計，故不要使用過多潛源參與計算，一般選擇目標區周邊300～400 km 范圍內的潛源即可滿足要求；②盡量使用較高配置的計算機，尤其是硬盤與內存空間充足，能滿足大量數據的一次性輸出與存儲需求；③若計算點過多，超過百萬級別，為了保持計算與輸出數據的穩定性，可考慮分批計算，最后對計算結果數據進行合并。計算完成后，使用簡短的Matlab 或Python 程序代碼提取出不同超越概率水平的基巖峰值加速度等數據，建議輸出為逗號分隔符的csv 格式，便于后續生成矢量數據文件。主要使用的matlab 代碼如下（以50年63%的峰值加速度為例）。

fidin = fopen('反應譜數據文件','r'); i=0;

while ～feof(fidin)

i=i+1;line1 = fgetl(fidin1);

PGA50a63=str2num(line1(12:20));

fprintf(fidout,'%10d,%8.3f ',[i PGA50a63]);

end

2.2 場地調整與地震危險性分級

根據地震危險性圖編制技術規范（FXPC/DZ P-01），地震危險性概率分析計算得到基巖峰值加速度，對應Ⅰ1類場地峰值加速度，使用全國1∶100 萬宏觀場地類別，按照表1 所給值采取分段線性插值確定調整系數，進行調整轉換為相應場地的峰值加速度值[3]。

表1 場地地震動峰值加速度調整系數Table 1 Adjustment coefficient of peak acceleration of site ground motion

根據年超越概率10-4的地震動峰值加速度（ax），將場地地震危險性分為4 級，其中1 級的危險程度最高，4 級的危險程度最低。分級原則為：1 級為ax≥760 gal、2 級為380 gal≤ax＜760 gal、3 級為190 gal≤ax＜380 gal、4 級為ax＜190 gal。值得注意的是，此次地震災害風險普查全國30″格網的地震危險性結果，分級原則在此基礎上有所調整。

地震動場地調整大致有以下3 種處理方法。

1）若完全基于ArcGIS 軟件平臺實現操作，則首先將基巖地震危險性計算結果轉換為矢量數據shp 點文件[6]，或直接存儲在gdb 數據庫中會提高運行速度。新建場地類別與各個超越概率場地峰值加速度、地震危險性等級等字段，使用空間連接工具賦值場地類別，編寫Python 代碼實現根據場地類別和線性插值獲取調整系數的函數，運用字段計算器對各個概率場地加速度值進行計算，再根據年超越概率10-4的地震動峰值加速度對地震危險性分級。當計算數據量達百萬級別時，在矢量數據shp 文件中處理字段有時會不易操作且耗時較長，后續對加速度值進行歸檔用于風險評估與區劃的輸入還需進行一系列處理。

2）完全使用Matlab 程序來處理，讀取宏觀各個場地類別數據的邊界，使用函數inpolygon 逐點判斷場地類別，編寫子函數線性插值獲得場地調整系數對峰值加速度進行調整、并判斷地震危險性等級。同時為了滿足后續繪圖分析以及地震風險評估所需地震危險性輸入，可根據場地峰值加速度按照一定間隔賦值等值線值和進行烈度歸檔，輸出所需數據文件，利用函數shapewrite 可生成矢量數據shp 文件。

3）由于不同場地類別都是由多個面對象組成，若在程序代碼架構方面不是特別完善，特別是處理沿海島嶼較多的地區，循環中嵌套循環，在判別場地類別存在一定的困難，且耗時較長。此時使用ArcGIS 軟件自帶的空間連接工具賦值該字段屬性效率更高，再提取場地類別并與峰值加速度數據進行合并處理，后續使用程序處理數據并生成矢量數據shp 文件。

2.3 矢量文件

按照地震災害風險普查工作的要求，地震危險性分析部分需要生成3 個矢量數據shp 文件即地震動PGA 圖（矢量點）、地震危險性等級圖（矢量點）、地震危險性等級圖（矢量面），在全國綜合風險普查系統中通過質檢才可上傳成功，且各shp 文件對字段與順序都有要求。

若前述數據操作均基于ArcGIS 軟件，則可直接選擇相應的字段數據輸出，低版本（如ArcGIS 10.2）在選取字段時不支持調整順序，可在其他軟件（如Access）中調整表字段順序，或者直接使用ArcGIS10.6 及以上的較高版本。若使用Matlab 程序代碼，利用函數shapewrite 則可快速生成所要求的數據；帶有后綴為shx、shp、dbf 三個文件是矢量數據必不可少的，加上后綴為prj 的投影文件即可，投影文件實際上是和其他數據同名的文件，程序操作時將所需投影文件內容寫入即可（如CGCS2000）[6]，故在程序中只需要輸出這4 個文件；輸出數據字段名可自行設置，建議為不超過8 位字符的英文名且符合shp 數據字段命名規則。需要注意的是矢量數據shp 文件存儲數據的上限為2 GB，超過則無法保存，需要存儲在gdb 數據庫中操作。為了順利生成矢量文件，需預先構建所需字段的結構體，否則當數據量超過十萬級別時生成shp 數據較為困難，所使用的matlab 代碼示例如下（其中Geometry、X、Y 為shp 文件必須，longitude、latitude、risk 為輸出字段，可根據需要增減）。

STR_P='struct(''Geometry'',values,''X'',values,''Y'',v alues,''longitude'',values,''latitude'',values,''risk'',values)';

values = cell(Grid_Points,1); P = eval(STR_P).

3 批量繪圖與生成文檔

3.1 批量繪圖

根據地震危險性分析結果，需要基于ArcGIS軟件平臺編制峰值加速度和地震危險性等級等圖件且制圖需滿足地震災害普查相關標準的要求。若要求繪制某省各個區縣相應的圖件，則可按照省級行政邊界整體進行地震危險性分析與計算處理，在ArcMap 中按照各區縣范圍實現自動出圖。具體操作步驟如下：加載地震危險性分析結果、行政邊界、水系、公路、鐵路等矢量圖層，啟動數據驅動并進行設置，以區縣矢量數據為索引圖層；勾選區縣名稱字段，在數據框中選擇剪裁至當前數據驅動頁面范圍，此時可以根據制圖需求排除不需剪裁的圖層，切換到布局視圖；設置動態文本標題，插入邊框、風險普查logo、指北針、比例尺等，調整到預期出圖效果，保存該地圖文檔mxd 文件；最后選擇pdf 格式批量導出地圖。4 個不同超越概率水平的基巖和場地地震危險性結果以及地震危險性等級單獨保存為不同的地圖文檔mxd 文件，逐個批量導出圖件。

為提高批量出圖效率需要注意以下幾點：①不同的地圖文檔mxd 中需設置統一的圖例；②水系、道路等數據量較大的圖層，剪裁至比目標區稍大范圍即可；③制圖所用矢量圖層均導入gdb 數據庫中會提高效率；④區縣數量較多時，在ArcMap 中批量導出較高分辨率的圖件，耗時較長且可能無法導出，尤其是低版本的ArcGIS，此時可在ArcGIS Pro 中導入mxd 文件再批量導出地圖；⑤若對ArcGIS Pro 較為熟悉，可按照上述類似步驟直接繪圖并批量導出，較為高效；⑥設置地圖導出至一個pdf 文件后，可批量轉換為圖片，輸出順序與數據驅動索引圖層字段內容一致，此時可用簡潔的代碼新建以各區縣名稱命名的文件夾，將對應的地圖圖件重命名并導入。

3.2 批量生成word 文檔

對于地震災害風險普查工作而言，同一地區不同區縣的地震危險性分析報告，內容上相對來說較為固定，主要包括任務概述、地震危險性分析原理、潛源劃分、地震活動性參數的確定、衰減關系、地震危險性計算與場地調整等幾個方面?？梢栽O置一個通用的word 文檔模板，利用Matlab 或Python等程序，構建代碼啟動Word 調用功能（Word.Application），批量修改報告名稱、替換區縣名稱、增加數據統計表格等內容。同時調用上述各區縣文件夾中的地震危險性圖件，設置剪裁圖件空白部分的參數、調整大小等處理后，插入相應的地震危險性word 報告中去，并添加各圖件標題。運行程序之前，需在Microsoft Word 中將當前文件夾設置為受信任位置，可避免反復出現信任提醒而導致中斷循環的操作。主要使用的matlab 代碼如下。

copyfile('word 模板文件名','所需word 文件名');

try

Word = actxGetRunningServer('Word.Application');

catch

Word = actxserver('Word.Application');

end

Word.Visible = 1;

Selection = Word.Selection;

Selection.InlineShapes.AddPicture('圖片路徑+名稱','True','True').

4 示例分析

以湖南省地震危險性計算與數據處理工作為例，需要生成上傳至風險普查系統中的3 個地震危險性矢量數據shp 文件，以及湖南省14 個地市州與122 個區縣地震危險性圖件與報告，下面從計算過程與處理平臺、效率等方面進行簡要分析。采用30″標準格網（約28 萬個計算點），計算出50年超越概率63%、10%、2%以及年超越概率10-4基巖峰值加速度并進行場地調整；批量繪制4 個超越概率的場地地震危險性圖和地震危險性等級圖（每個市或區縣均5 張圖件、共計680 張圖），批量生成每個市或區縣的地震危險性分析word 報告；所使用的計算機配置為Windows10 操作系統、內存容量32 GB、硬盤容量2 TB、inteli7 處理器。

表2 為湖南省地震危險性分析工作的具體操作步驟、處理平臺與處理效率，其中地震危險性分析計算使用SEC R2019 軟件，其他數據處理與分析步驟主要使用Matlab、Python、ArcGIS 等軟件?？梢钥闯?，利用程序編寫代碼處理數據效率很高、靈活性強且易于操作，生成矢量數據shp 文件也比在ArcGIS 軟件中加載數據再導出更加高效；利用ArcGIS 數據驅動批量出圖對這種可重復性工作值得推薦；最后利用程序代碼實現批量生成14 個市和122 個區縣的word 報告，并將共計680 張圖件批量插入相應的word 中相應位置并增加圖名，極其高效地完成了工作。

表2 湖南省地震危險性分析與數據處理一覽表Table 2 List of seismic risk analysis and data processing in Hunan province

5 討論與結論

本文沒有涉及地震危險性概率分析方法原理及相關內容，但在實際工作中需要深入分析，比如潛源參數、地震帶參數取值[7]、衰減關系適用性、不同超越概率峰值加速度間比例關系等。在大量數據計算和處理之前，可先在目標區內選取一定數量的控制點進行地震危險性計算，并對不同超越概率水平的峰值加速度與潛源貢獻等結果進行分析，綜合判定計算結果的合理性。若需要計算反應譜，可根據數據計算量情況對不同周期點單獨或合并計算。部分地區還需要特殊處理，如湖南省南部跨越長江中游地震帶與華南沿海地震帶，分別對應中強地震區與東部活躍區的衰減關系[8]，可分別計算并取峰值加速度的較大值做為最終結果。

地震危險性等級圖（矢量面）可根據計算結果，使用ArcGIS 軟件的克里金插值和等值線工具結合實現，也結合其他專業軟件（如Surfer）提供的不同算法生成等值線。對于危險性等級不復雜的情況，使用ArcGIS 軟件的剪切面工具手工操作，可自動滿足拓撲檢查的要求，也可結合ArcGIS 二次開發生成矢量面，此問題需要進一步研究解決。

需要具備一定的編程基礎才能解決可能出現的各種問題并高效地完成工作。Matlab 和Python程序在數據處理領域便于使用，且ArcGIS 軟件內嵌Python 程序易于二次開發。如對于湖南省地震危險性分析工作，需要在全省數據基礎上得到14 個地市州與122 個區縣的地震動PGA 和地震危險性等級的矢量數據shp 文件，可以使用python 代碼調用剪裁工具arcpy.Clip_analysis ('輸入shp','剪裁面shp','輸出shp')循環實現[9]。另外，大量繪圖工作均基于ArcGIS 軟件平臺，當危險性矢量點數據量很大時操作不便，可將點轉為柵格進行處理。

地震災害風險普查工作中的危險性概率計算數據量往往很大，為了避免存儲空間不足、內存溢出等問題，使用高配置的計算機是十分必要的，在計算與后續數據處理中均需考慮是否分批進行。筆者在前述配置計算機上做過測試，單次成功計算并處理約500 萬個計算點的數據，且利用Matlab 程序代碼生成此級別的矢量點數據shp 文件只需幾分鐘，更多的數據只要內存不溢出均可處理。需要注意的是單個矢量數據shp 文件的大小上限為2 GB[10]，超出則需要在gdb 數據庫中存儲和處理。另外，在較粗格網地震危險性計算結果基礎上，使用克里金插值等方式獲取較細格網數據，可大大減少地震危險性計算的時間，但需要對結果的合理性進行論證。

文中所述數據處理、批量出圖與生成word 文檔等方法，在地震安全性評價工作也可參考。如線狀工程的沿線區劃工作，可參考標準格網建立所使用的方法，外加線狀工程一定距離緩沖區加以約束，可快速生成區劃范圍所需間隔的計算點。在區域性地震安全性評價工作中，大量的地震動時程與規準譜標定圖件，亦可利用程序構建簡潔的代碼調用Word.Application、設置圖片大小等參數批量快速的導入word 報告。

文中所述來源于筆者實際工作中的探索，并寫了相應的Matlab 和Python 程序代碼，由于水平有限可能存在不足，在后續工作中與感興趣科研工作者交流分享。

致謝 非常感謝地震災害風險普查專題培訓工作中專家們的悉心指導及工作中同事們的幫助，也感謝審稿專家提出的寶貴建議。