孫 強(qiáng)，陳時(shí)軍，張正帥

（1.淄博市應(yīng)急救援指揮保障服務(wù)中心，山東 淄博 255086；2.淄博遙測地震臺(tái)網(wǎng)，山東 淄博 255086；3.山東省地震局，山東 濟(jì)南 250014）

地震災(zāi)害是眾多自然災(zāi)害中的一種，且一旦發(fā)生容易形成較長的災(zāi)害鏈，造成極其嚴(yán)重的損失，因此地震災(zāi)害風(fēng)險(xiǎn)研究是災(zāi)害風(fēng)險(xiǎn)研究的重要內(nèi)容。聯(lián)合國國際十年減災(zāi)委員會(huì)秘書處為了促進(jìn)全球尤其是發(fā)展中國家的城市減災(zāi)活動(dòng)，實(shí)施了城市地區(qū)抗御地震災(zāi)害風(fēng)險(xiǎn)評估項(xiàng)目[1]，自1997年開始，在全球選取10個(gè)大城市作為研究對象，利用風(fēng)險(xiǎn)評價(jià)工具來評估建筑物、生命線工程的震害情況。1997年美國斯坦福大學(xué)的Davidon R[2]提出利用地震災(zāi)害風(fēng)險(xiǎn)指數(shù)來評價(jià)不同城市潛在地震災(zāi)害相對嚴(yán)重程度的方法，這是歷史上第一次采用定量的方法來比較各城市之間的地震風(fēng)險(xiǎn)水平以及分析不同因素對地震風(fēng)險(xiǎn)的貢獻(xiàn)；同年，陳棋福[3]也開展了類似研究，評價(jià)了不同城市潛在地震災(zāi)害的相對嚴(yán)重程度；黃崇福[4]等用模糊集理論建立了城市地震災(zāi)害風(fēng)險(xiǎn)評價(jià)的數(shù)學(xué)模型；陳棋福[5-7]等提出了基于GDP和人口資料的地震損失預(yù)測評估方法，并對全球震災(zāi)資料進(jìn)行了分析，得到了全球地震損失分布圖，該方法可在沒有詳細(xì)建筑分類數(shù)據(jù)的地區(qū)進(jìn)行震損評估；聶高眾[8]等提出了地震災(zāi)害絕對風(fēng)險(xiǎn)計(jì)算公式，但多數(shù)情況下工作量較大，估計(jì)地區(qū)地震風(fēng)險(xiǎn)較困難；姚清林[9]等基于建筑物的破壞程度和區(qū)域破壞率，通過烈度的模糊推理構(gòu)建了風(fēng)險(xiǎn)評價(jià)因子——震損度（DA），并通過模糊信息分配技術(shù)給出了DA值的計(jì)算方法，較好地實(shí)現(xiàn)了評價(jià)因子不確定性信息的傳遞；唐麗華[10]等分別利用指標(biāo)體系法和地震災(zāi)害風(fēng)險(xiǎn)系數(shù)法對烏魯木齊市開展了地震災(zāi)害風(fēng)險(xiǎn)評估，并進(jìn)行了對比；侯林鋒[11]等運(yùn)用基于指標(biāo)體系的方法對浙江省地震災(zāi)害風(fēng)險(xiǎn)進(jìn)行了評估。自然災(zāi)害風(fēng)險(xiǎn)建模和評估方法主要包括指標(biāo)體系法、基于GIS的方法、概率法和情景建模法[12]4種，前人研究成果多采用情景建模法和指標(biāo)體系法，情景建模法適用于建筑物分類等基礎(chǔ)資料詳細(xì)、豐富的地區(qū)進(jìn)行中小尺度評估，指標(biāo)體系法則需選擇眾多細(xì)化的指標(biāo)來進(jìn)行大中尺度的風(fēng)險(xiǎn)評估。

本文采用聯(lián)合國對自然災(zāi)害風(fēng)險(xiǎn)的定義[13]，即地震災(zāi)害風(fēng)險(xiǎn)度可表達(dá)為地震危險(xiǎn)度和承災(zāi)體易損度的乘積。地震危險(xiǎn)度的估計(jì)結(jié)果對地震災(zāi)害風(fēng)險(xiǎn)評估的影響較大。在地震安全性評價(jià)中，通常采用美國Cornell 1968年提出的地震危險(xiǎn)性分析的概率估計(jì)法[14]計(jì)算地震危險(xiǎn)度，其工作繁雜，需要的原始資料較多，且在潛在震源區(qū)劃分、地震活動(dòng)空間分布系數(shù)確定等方面存在較大的不確定性。本文利用CF模型和Logistic回歸模型分析了研究區(qū)發(fā)生烈度6度及以上地震的概率，并從地震烈度發(fā)生概率的角度計(jì)算了研究區(qū)地震危險(xiǎn)度，方法與過程較簡明，其結(jié)果通過檢驗(yàn)證明了方法的合理性。在承災(zāi)體易損度的計(jì)算中，本文首先利用衛(wèi)星遙感數(shù)據(jù)提取研究區(qū)的土地利用類型，利用騰訊大數(shù)據(jù)的定位工具獲得研究區(qū)某時(shí)段的最大人口密度，利用研究區(qū)2019年統(tǒng)計(jì)年鑒資料獲取當(dāng)?shù)氐腉DP數(shù)據(jù)，數(shù)據(jù)獲取方式較簡便，且保證了一定的精度；然后劃分?jǐn)?shù)據(jù)等級，利用GIS手段疊加成圖；最終完成研究區(qū)地震災(zāi)害風(fēng)險(xiǎn)的評估。

1 地震災(zāi)害危險(xiǎn)度

1.1 研究方法

1.1.1 確定性系數(shù)（CF） 模型

CF模型是由Shortliffe E H[15]等提出并由Heckerman D[16]改進(jìn)的一個(gè)概率函數(shù)，假設(shè)將來發(fā)生地震災(zāi)害的條件和過去發(fā)生地震災(zāi)害的條件相同，則其計(jì)算公式為：

式中，PPa為地震災(zāi)害在因子分類數(shù)據(jù)a中發(fā)生的條件概率，在實(shí)際研究中通常用因子分類a的地震災(zāi)害面積與數(shù)據(jù)分類a的面積比表示；PPs為研究區(qū)的災(zāi)害總面積與研究區(qū)總面積的比值。

由式（1）可知，CF的變化區(qū)間為[-1,1]，其中正值表示事件發(fā)生的確定性增長，即易于發(fā)生地震災(zāi)害；負(fù)值表示事件發(fā)生的確定性降低，即不易于發(fā)生地震災(zāi)害。

1.1.2 Logistic回歸模型

Logistic回歸模型是研究二分類結(jié)果與其影響因子之間關(guān)系的一種多元統(tǒng)計(jì)分析方法，描述的是二元因變量（通常0表示地震災(zāi)害不發(fā)生，1表示地震災(zāi)害發(fā)生）和一系列自變量（x1，x2，…，xn）之間的關(guān)系。自變量可以是連續(xù)的也可以是離散的，無需滿足正態(tài)頻率分布。Logistic回歸函數(shù)的表達(dá)式為：

式中，P為地震災(zāi)害發(fā)生的概率，取值范圍為0～1；α為截距；β為回歸系數(shù)。

1.2 研究區(qū)概況

淄博市位于山東省中部（圖1）、我國東部著名的郯廬地震帶西側(cè)，境內(nèi)地質(zhì)構(gòu)造復(fù)雜，斷裂發(fā)育。淄博市域面積為5 965 km2，人口為433.96萬人，人口平均密度為728人/km2，包括張店、臨淄、淄川、博山、周村5個(gè)市轄區(qū)和桓臺(tái)、高青、沂源3個(gè)市轄縣，2018年全市地區(qū)生產(chǎn)總值為3 572.98億元。淄博市工業(yè)發(fā)達(dá)，人口眾多，人口密度大。隨著全球地震形勢的嚴(yán)峻，地震災(zāi)害風(fēng)險(xiǎn)的不確定性也在增大，為了更好地把握淄博市的震情形勢和地震災(zāi)害風(fēng)險(xiǎn)情況，對淄博市進(jìn)行地震災(zāi)害風(fēng)險(xiǎn)評估研究很有必要。

圖1 研究區(qū)位置與地震震中圖

1.3 地震災(zāi)害危險(xiǎn)度計(jì)算的相關(guān)因素與分級

1.3.1 淄博市歷史地震綜合烈度圖（6度及以上）

本文根據(jù)山東省地震預(yù)報(bào)研究中心提供的地震目錄，選取有歷史記載以來（公元408年）的淄博及其周邊Ms3.0級以上的地震284個(gè)，對于有明確發(fā)震斷裂的地震，采用式（3）～（6）[17]繪制地震烈度6度及以上（PGA≥0.05 g）的橢圓區(qū)域；其余地震采用式（7）[18]繪制地震烈度6度及以上（PGA≥0.05 g）的圓形區(qū)域，最終形成6度及以上的綜合烈度圖，如圖2所示。

圖2 淄博市6度及以上綜合烈度圖

1.3.2 50 km內(nèi)的最大地震

地震的發(fā)生都有一定的烈度范圍，根據(jù)較大地震的烈度影響范圍，某點(diǎn)周邊50 km內(nèi)的最大地震基本可以確定6度及以上烈度的影響范圍。淄博市東側(cè)的郯廬斷裂帶歷史上曾多次發(fā)生強(qiáng)震，因此50 km內(nèi)的最大地震這一指標(biāo)可反映研究區(qū)的地震危險(xiǎn)性。根據(jù)該指標(biāo)將淄博市范圍劃分為50 km內(nèi)最大地震為5～6級和50 km內(nèi)最大地震為6級以上兩個(gè)區(qū)域（圖3a）。

1.3.3 地震頻度

地震頻度可反映某地點(diǎn)地震活動(dòng)水平的高低。本文計(jì)算地震頻度參考王健[19]給出的公式，我國地震觀測精度e約為2～3 km，小于該距離已無意義，因此當(dāng)震中距Di＜e時(shí)，取Di=e，這樣既可體現(xiàn)地震數(shù)量，又可顯示震中位置的空間集中程度。在GIS中計(jì)算各網(wǎng)格點(diǎn)的地震頻度，再利用自然間斷點(diǎn)法將其劃分為3個(gè)等級：低頻度0.318 3～1.799 6，中等頻度1.799 7～2.568 9，高頻度2.569 0～4.571 4（圖3b）。

式中，Vj為第j個(gè)網(wǎng)格點(diǎn)的地震頻度；ln為自然對數(shù)；Di為第i個(gè)地震到網(wǎng)格點(diǎn)j的距離；rm為以網(wǎng)格點(diǎn)為中心的搜索半徑，這里取30 km；n為網(wǎng)格點(diǎn)周邊30 km范圍內(nèi)的地震個(gè)數(shù)。

1.3.4 斷層長度值

通常情況下，斷層越長其孕震能力越強(qiáng)，即該斷層的控震震級越高、地震危險(xiǎn)性越大；反之則越小。本文首先按照斷層長度對其危險(xiǎn)性進(jìn)行10分制評分（表1），然后利用式（9）[20]計(jì)算斷層長度值，最后將結(jié)果等間距分為5級，即[0,2]、(2,4]、(4,6]、(6,8]、(8,10]，如圖3c所示。

表1 斷層長度評分表

式中，Yj為第j個(gè)網(wǎng)格的斷層長度值；Lj為網(wǎng)格周邊20 km內(nèi)的斷層長度評分值；Dj為網(wǎng)格中心到斷層的距離。

1.3.5 斷層時(shí)代值

斷層的活動(dòng)時(shí)代不同，其控震能力也有所差異。根據(jù)王華林[21]等的研究，全新世斷裂的控震能力在6.8級以上，晚更新世在6.1級左右，早中更新世在5.3級左右，由此可見，活動(dòng)時(shí)代越早控震能力越弱。首先據(jù)此對斷層的活動(dòng)時(shí)代進(jìn)行10分制評分，年代越久遠(yuǎn)賦值越小，具體賦值情況如表2所示；然后利用式（10）[20]計(jì)算研究區(qū)每個(gè)網(wǎng)格的斷層時(shí)代值；最后將結(jié)果等間距分為5級，[0,2]、(2,4]、(4,6]、(6,8]、(8,10]，如圖3d所示。

表2 斷層活動(dòng)時(shí)代評分表

式中，Zk為第k個(gè)網(wǎng)格的斷層時(shí)代值；Tk為網(wǎng)格周邊20 km內(nèi)的斷層活動(dòng)時(shí)代評分值；Dk為網(wǎng)格中心到斷層的距離。

1.3.6 場地分類

不同的場地類型在同等震級地震的作用下會(huì)造成不同的災(zāi)害嚴(yán)重程度，通常軟土場地的震害比基巖嚴(yán)重。本文利用Wald D J[22]等提出的地形坡度法來劃分場地類型。地形坡度法首先利用地形坡度與地表以下30 m內(nèi)的平均剪切波速（(30)）的關(guān)系來推導(dǎo)區(qū)域性(30)的分布，再按照相應(yīng)的規(guī)范劃分場地類別。本文首先利用地理空間數(shù)據(jù)云（www.gscloud.cn）提供的GDEMV2 30 m分辨率的數(shù)字高程數(shù)據(jù)提取坡度值，再根據(jù)美國的NEHRP規(guī)范（表3）來劃分場地類別（圖3f）。

圖3 地震災(zāi)害影響因素分級圖

表3 (30)和坡度值分類表

表3 (30)和坡度值分類表

場地類型E D C B-Vs（30）/（m/s）＜180[180，360）[360，760]＞760坡度的正切值（tan）＜1.0E-4[1.0E-4，0.018）[0.018，0.138]＞0.138

1.4 地震災(zāi)害危險(xiǎn)度評價(jià)

本文以各分類級別的CF值為自變量，以是否發(fā)生6度及以上烈度地震災(zāi)害為因變量（1表示發(fā)生災(zāi)害，0表示未發(fā)生災(zāi)害），采用統(tǒng)計(jì)分析軟件進(jìn)行二項(xiàng)邏輯回歸分析；再將得到的各因素回歸系數(shù)代入回歸方程，計(jì)算各網(wǎng)格點(diǎn)發(fā)生6度及以上烈度地震災(zāi)害的概率值。本文研究的基本單元格大小為100 m×100 m，將淄博市剖分為600 781個(gè)網(wǎng)格單元，計(jì)算過程中各圖層按照100 m×100 m的柵格進(jìn)行重采樣。具體計(jì)算步驟為：

1）利用CF模型計(jì)算6個(gè)影響因素各分類級別的CF值。根據(jù)淄博市6度及以上烈度地震災(zāi)害的分布信息，利用CF模型可計(jì)算得到各分類級別在因素中的相對權(quán)重。該權(quán)重不僅可以比較同一評價(jià)因素下各分類級別的相對重要性，而且可以比較不同評價(jià)因素的分類級別。各影響因素分類級別的CF值如表4所示。

表4 各影響因素分類級別的CF值

2）利用Logistic回歸模型計(jì)算各因素的權(quán)重值。本文以各影響因素分類級別的CF值為自變量，以是否發(fā)生6度及以上烈度地震災(zāi)害為因變量（1為發(fā)生災(zāi)害，0為未發(fā)生災(zāi)害），隨機(jī)選取發(fā)生6度及以上烈度地震災(zāi)害的80%的單元格和未發(fā)生6度及以上烈度地震災(zāi)害的80%的單元格作為訓(xùn)練樣本，共計(jì)480 624個(gè)單元格，利用統(tǒng)計(jì)分析軟件進(jìn)行二項(xiàng)邏輯回歸分析。邏輯回歸分析結(jié)果如表5所示。在地震危險(xiǎn)性評價(jià)過程中，B表示各因素權(quán)重大小，每個(gè)變量在方程中的重要性通過比較Wals值（或Sig值）來判斷，Wals值大者（或Sig值小者）顯著性較高，在方程中就更重要，Sig值小于0.05才有統(tǒng)計(jì)意義。由表5可知，各因素的Sig值均小于0.05，說明所選因素都是有效的。

表5 邏輯回歸分析結(jié)果

3）地震災(zāi)害危險(xiǎn)度區(qū)劃。將計(jì)算得到的各因素回歸系數(shù)代入式（2），得到邏輯回歸方程：

式中，P為發(fā)生6度及以上烈度地震災(zāi)害的概率，取值范圍為[0,1]；X1為50 km內(nèi)最大地震指標(biāo)各分類級別的CF值；X2為地震頻度指標(biāo)各分類級別的CF值；X3為斷層長度值指標(biāo)各分類級別的CF值；X4為斷層時(shí)代值指標(biāo)各分類級別的CF值；X5為場地分類指標(biāo)各分類級別的CF值。

根據(jù)式（11）可計(jì)算得到淄博市發(fā)生6度及以上烈度地震災(zāi)害的概率，再按照等間距法將概率區(qū)劃分為極低易發(fā)區(qū)（0～0.25）、低易發(fā)區(qū)（0.25～0.5）、中易發(fā)區(qū)（0.5～0.75）和高易發(fā)區(qū)（0.75～1）4類，如圖4所示。

圖4 淄博市地震危險(xiǎn)度區(qū)劃圖

4）利用ROC曲線對地震危險(xiǎn)度區(qū)劃結(jié)果進(jìn)行驗(yàn)證。本文繪制ROC曲線所用數(shù)據(jù)為去除80%訓(xùn)練數(shù)據(jù)后剩余的20%數(shù)據(jù)，包括淄博市已發(fā)生6度及以上烈度地震的20%單元格和未發(fā)生6度及以上烈度地震的20%單元格，共計(jì)120 157個(gè)單元格。Logistic回歸模型的ROC曲線如圖5所示，可以看出，AUC值為0.804，說明采用CF模型和Logistic回歸模型能較客觀準(zhǔn)確地評價(jià)淄博市的地震危險(xiǎn)度。

圖5 Logistic回歸模型的ROC曲線

2 承災(zāi)體易損度

本文采用土地利用類型、人口熱力值（代表相對人口密度）和各區(qū)縣GDP值表示承災(zāi)體易損度。

2.1 土地利用類型

研究表明，地震災(zāi)害損失與土地利用類型關(guān)系密切[23]。人口傷亡主要是由于房屋倒塌造成的，因此居住用地、工礦用地、公共設(shè)施用地的人口易損性較高，而荒地、水域的易損性較低。本文利用Landsat8的衛(wèi)星影像提取了淄博市的土地利用類型，包括城鎮(zhèn)居住用地、工礦用地、耕地、水域和山地（圖6a）；并對土地利用類型的地震易損性進(jìn)行了評分，具體結(jié)果如表6所示。

表6 土地利用類型地震易損性評分表

2.2 人口熱力值

本文以人口熱力值代表相對人口密度，通過爬取騰訊位置大數(shù)據(jù)獲取了淄博市0.05°×0.05°網(wǎng)格的人口熱力值；利用GIS的核密度分析工具，設(shè)置輸出像元大小為100 m、搜索半徑為5 000 m，得到淄博市人口熱力圖；利用GIS的自然間斷點(diǎn)法將人口熱力值劃分為10檔，并由小到大依次賦于1～10的評分值（圖6b）。

2.3 各區(qū)縣GDP值

本文從2019年淄博統(tǒng)計(jì)年鑒[24]中獲取了淄博市各區(qū)縣2018年度的GDP值（圖6c）；并對各區(qū)縣GDP值進(jìn)行10分制評分，1 000億及以上的GDP值設(shè)置為10分，不足1 000億的GDP值按比例設(shè)置為0～9分。

圖6 承災(zāi)體易損度分級圖

3 淄博市地震災(zāi)害風(fēng)險(xiǎn)區(qū)劃

根據(jù)聯(lián)合國對自然災(zāi)害風(fēng)險(xiǎn)的定義，淄博市地震災(zāi)害風(fēng)險(xiǎn)值的計(jì)算公式為：

淄博市地震災(zāi)害風(fēng)險(xiǎn)值=淄博市地震災(zāi)害危險(xiǎn)度×土地利用分類值×人口熱力值×淄博各區(qū)縣GDP值

式中，淄博市地震災(zāi)害危險(xiǎn)度的取值范圍為[0,1]，其余3項(xiàng)的取值范圍均為1～10。

本文采用GIS的自然間斷點(diǎn)法將淄博市地震災(zāi)害風(fēng)險(xiǎn)區(qū)劃分為4個(gè)等級：極低風(fēng)險(xiǎn)區(qū)（0.203 664～43.115 580）、低風(fēng)險(xiǎn)區(qū)（43.115 581～124.237 066）、中等風(fēng)險(xiǎn)區(qū)（124.237 067～321.078 038）和高風(fēng)險(xiǎn)區(qū)（321.078 039～825.554 976）。其中，極低風(fēng)險(xiǎn)區(qū)面積占淄博市總面積的85.18%，低風(fēng)險(xiǎn)區(qū)面積占比為9.91%，中等風(fēng)險(xiǎn)區(qū)面積占比為3.55%，高風(fēng)險(xiǎn)區(qū)面積占比為1.36%，如圖7所示。

圖7 淄博市地震災(zāi)害風(fēng)險(xiǎn)區(qū)劃圖

4 結(jié)語

本文根據(jù)地震災(zāi)害風(fēng)險(xiǎn)的定義，首先利用CF模型和Logistic回歸模型計(jì)算了地震災(zāi)害危險(xiǎn)度，占用篇幅較大、闡述較多，突出了地震危險(xiǎn)度計(jì)算在地震災(zāi)害風(fēng)險(xiǎn)計(jì)算中的重要性，其中CF模型較好地解決了不同類別指標(biāo)各級別因子在同一標(biāo)準(zhǔn)下的定量比較問題，且避免了類別因子處理時(shí)的主觀干擾，Logistic回歸模型則根據(jù)各因素與歷史地震災(zāi)害數(shù)據(jù)之間的關(guān)系確定了各因素的權(quán)重[25]；然后確定了土地利用類型、人口密度和區(qū)縣GDP數(shù)據(jù)等承災(zāi)體易損度指標(biāo)，這部分?jǐn)?shù)據(jù)來源于衛(wèi)星影像數(shù)據(jù)、騰訊大數(shù)據(jù)和淄博統(tǒng)計(jì)年鑒，數(shù)據(jù)的精確度和分辨率有了很大提升；最后利用GIS的空間疊加功能得到了淄博市地震災(zāi)害風(fēng)險(xiǎn)區(qū)劃圖，其中高風(fēng)險(xiǎn)區(qū)面積占比為1.36%，中等風(fēng)險(xiǎn)區(qū)面積占比為3.55%，低風(fēng)險(xiǎn)區(qū)面積占比為9.91%，極低風(fēng)險(xiǎn)區(qū)面積占比為85.18%，該結(jié)果在一個(gè)少震地區(qū)是比較合理的，基本反映了淄博市的地震災(zāi)害風(fēng)險(xiǎn)狀況。本文得到的風(fēng)險(xiǎn)區(qū)劃圖可隨年度數(shù)據(jù)的更新而更新，且本文方法可應(yīng)用于其他地區(qū)，包括地震多發(fā)地區(qū)，為當(dāng)?shù)胤勒饻p災(zāi)工作提供參考。