基于GIS和logistic模型的地震滑坡致死人數(shù)快速評(píng)估方法

白仙富 聶高眾 葉燎原 戴雨芡 余慶坤

1)云南師范大學(xué)地理學(xué)部， 昆明 650500 2)云南省地震局， 昆明 650224 3)中國(guó)地震局地質(zhì)研究所， 北京 100029

0 引言

在高山峽谷地區(qū)， 地震通常會(huì)引發(fā)嚴(yán)重的山體滑坡(Pain， 1972)， 且滑坡體將對(duì)下方的居民點(diǎn)造成沖擊， 使其成為地震導(dǎo)致人員傷亡的主要因素之一。美國(guó)阿拉斯加大地震導(dǎo)致130人死亡， 其中48人死于地震滑坡(Keeper， 1984)。1969—1993年20多a間， 日本M≥6.9地震造成的死亡人數(shù)中一半以上源于地震滑坡(Kobayashi， 1994)。1970年的秘魯欽博特7.6級(jí)地震引發(fā)了迄今為止國(guó)外最大的地震滑坡災(zāi)難性事件， 此次地震滑坡至少造成2萬(wàn)人死亡(Plafker， 1971)。地震滑坡在中國(guó)也十分常見， 特別是西部山區(qū)， 在該區(qū)域發(fā)生的幾乎所有的強(qiáng)烈地震都會(huì)引發(fā)大量滑坡， 并造成不同程度的人員傷亡， 簡(jiǎn)要統(tǒng)計(jì)結(jié)果如表1 所示。

表1 中國(guó)部分地震滑坡災(zāi)害導(dǎo)致的人員死亡列表Table1 Death toll from the historical earthquake-induced landslides disasters in China

從歷史統(tǒng)計(jì)來(lái)看， 在中國(guó)西部山區(qū)地震導(dǎo)致的人員傷亡中， 地震滑坡致死人數(shù)可能占總死亡人數(shù)的25%～30%。在個(gè)別震例中， 地震滑坡導(dǎo)致的死亡人員數(shù)甚至超過建筑物破壞導(dǎo)致的死亡人員數(shù)(高慶華等， 2011)； 而一些震級(jí)較小的死亡人數(shù)為1至數(shù)人的地震事件， 人員死亡可能僅由地震滑坡所致， 說明地震滑坡已經(jīng)成為中國(guó)西部山區(qū)人民面臨的主要風(fēng)險(xiǎn)隱患之一。結(jié)合《中國(guó)地震動(dòng)參數(shù)區(qū)劃圖(GB18306-2015)》推測(cè)， 未來(lái)中國(guó)西部地區(qū)仍面臨著較高的地震滑坡人員傷亡風(fēng)險(xiǎn)。要提高中國(guó)西部山區(qū)地震人員傷亡快速評(píng)估的準(zhǔn)確性， 除了繼續(xù)發(fā)展和完善建筑物破壞導(dǎo)致人員傷亡的評(píng)估技術(shù)之外， 必須盡快建立和發(fā)展針對(duì)地震地質(zhì)災(zāi)害導(dǎo)致人員傷亡的評(píng)估方法(白仙富等， 2014)?？茖W(xué)評(píng)估地震滑坡災(zāi)害可能造成的人員傷亡， 是提高地震人員傷亡評(píng)估的準(zhǔn)確性進(jìn)而完善地震災(zāi)害損失評(píng)估體系的關(guān)鍵科學(xué)問題之一， 同時(shí)也是一直困擾地震災(zāi)害風(fēng)險(xiǎn)評(píng)估和震后高精準(zhǔn)快速評(píng)估的技術(shù)難題之一。

發(fā)展地震滑坡災(zāi)害導(dǎo)致人員傷亡的評(píng)估方法主要包括2個(gè)階段： 第1階段是探討地震發(fā)生后評(píng)估地震滑坡災(zāi)害危險(xiǎn)性的方法， 即回答“什么樣的地震會(huì)導(dǎo)致滑坡災(zāi)害”以及“可能在什么地點(diǎn)發(fā)生滑坡災(zāi)害”的問題。眾多研究者對(duì)地震滑坡災(zāi)害的危險(xiǎn)性評(píng)估進(jìn)行了多年探索和研究， 發(fā)展出多種方法， 已經(jīng)從定性評(píng)估逐漸走向半定量、 定量化評(píng)估(陳曉利， 2011； Yuanetal.， 2013； 劉甲美等， 2016； 許沖等， 2019)。第2階段是研究地震后出現(xiàn)的這些滑坡災(zāi)害是否會(huì)導(dǎo)致人員傷亡， 以及可能導(dǎo)致多少人員傷亡， 即地震滑坡災(zāi)害的危害性評(píng)估問題。歷史上在進(jìn)行地震人員死亡統(tǒng)計(jì)時(shí)， 通常只關(guān)注建筑物造成的人員死亡情況； 偶有針對(duì)地震滑坡致死的考察， 但也幾乎不記載滑坡遇難者的具體位置等信息， 使得我們難以進(jìn)行以地震參數(shù)、 滑坡位置、 居民點(diǎn)位置與地震滑坡致死人數(shù)的關(guān)系為對(duì)象的研究。也正是因?yàn)闅v史記錄的不完善， 至今仍很少見到關(guān)于地震滑坡災(zāi)害危害性方面的正式研究文獻(xiàn)， 導(dǎo)致長(zhǎng)期以來(lái)在地震損失快速評(píng)估或重點(diǎn)危險(xiǎn)區(qū)損失預(yù)評(píng)估等工作中， 對(duì)地震滑坡致死人數(shù)的預(yù)測(cè)更多還只能停留在專家經(jīng)驗(yàn)估計(jì)的層面上， 無(wú)法保證其準(zhǔn)確性。對(duì)該問題的關(guān)注和研究正是本文的主要目標(biāo)。

1 研究思路與方法

1.1 基本思路

圖 1 全區(qū)域和分區(qū)域人口信息示意圖Fig. 1 A sketchmap showing the demographic information in the whole region and the sub-regions.

根據(jù)上述推理認(rèn)識(shí)， 本文提出一種基于GIS和logistic模型的地震滑坡致死人數(shù)快速評(píng)估方法。該方法的基本思路(圖 2)是： 各公里網(wǎng)格單元的地震滑坡致死率與其內(nèi)部的滑坡危險(xiǎn)性等級(jí)屬性和數(shù)量相關(guān)， 并據(jù)此建立公里網(wǎng)格單元的地震滑坡致死率模型； 評(píng)估區(qū)域內(nèi)的地震滑坡致死人數(shù)等于區(qū)域范圍內(nèi)所有公里網(wǎng)格單元的地震滑坡致死人數(shù)之和。

圖 2 地震滑坡致死人數(shù)評(píng)估的基本思路Fig. 2 Basic ideas of the evaluation of the death toll from earthquake-induced landslides.

本研究所提出的地震滑坡致死人數(shù)評(píng)估方法的關(guān)鍵點(diǎn)在于構(gòu)建針對(duì)公里網(wǎng)格單元尺度的地震滑坡致死率模型。由于難以準(zhǔn)確評(píng)估某個(gè)確定地點(diǎn)內(nèi)是否發(fā)生滑坡， 因此， 選取地震滑坡危險(xiǎn)性作為網(wǎng)格單元的地震滑坡致死率影響因子。地震滑坡危險(xiǎn)性作為影響因子， 采用公里網(wǎng)格單元內(nèi)不同滑坡危險(xiǎn)性屬性像元數(shù)量的線性組合表示， 通過邏輯回歸求解線性組合中的未知參數(shù)建立基于公里網(wǎng)格單元的地震滑坡致死率模型。在此基礎(chǔ)上引入公里網(wǎng)格單元內(nèi)的人口數(shù)量， 最終建立基于公里網(wǎng)格單元的地震滑坡致死數(shù)量評(píng)估模型。本研究中， 用GIS統(tǒng)計(jì)各公里網(wǎng)格單元內(nèi)不同滑坡危險(xiǎn)性屬性的像元數(shù)量并計(jì)算網(wǎng)格單元的地震滑坡致死人數(shù)及死亡率等， 通過獲取這些參數(shù)提出公里網(wǎng)格的地震滑坡致死率與其內(nèi)部各屬性值的滑坡危險(xiǎn)性像元數(shù)量之間的聯(lián)系：

Mi=f(Xil)

(1)

式中，Mi表示第i個(gè)公里網(wǎng)格單元的地震滑坡致死率，f為推斷的影響因素與Mi的關(guān)系，Xil表示第i個(gè)公里網(wǎng)格單元內(nèi)各屬性值的滑坡危險(xiǎn)性像元數(shù)量。

按照上述思路， 針對(duì)公里網(wǎng)格地域單元的地震滑坡致死人數(shù)快速評(píng)估方法一般分為3個(gè)步驟(圖 3)： 1)在GIS中將研究區(qū)劃分為1km×1km的網(wǎng)格單元， 把網(wǎng)格內(nèi)的人口數(shù)量信息賦值給網(wǎng)格作為其人口屬性值； 2)計(jì)算網(wǎng)格單元的地震滑坡致死率。通過GIS的空間分析功能從預(yù)測(cè)數(shù)據(jù)集中提取研究區(qū)的地震滑坡危險(xiǎn)性數(shù)據(jù)， 并統(tǒng)計(jì)網(wǎng)格單元內(nèi)各滑坡危險(xiǎn)性屬性值的像元數(shù)量作為網(wǎng)格單元的滑坡危險(xiǎn)性屬性， 根據(jù)滑坡危險(xiǎn)性屬性基于統(tǒng)計(jì)的logistic模型計(jì)算網(wǎng)格單元的地震滑坡致死率； 3)評(píng)估公里網(wǎng)格單元地震滑坡致死人數(shù)分布和整個(gè)研究區(qū)地震滑坡致死人數(shù)。

圖 3 地震滑坡致死人數(shù)的快速評(píng)估流程Fig. 3 Procedure of the rapid assessment of death toll from earthquake-induced landslides.

1.2 研究方法

1.2.1 公里網(wǎng)格人口屬性賦值

本研究使用中國(guó)地震局地質(zhì)研究所和中國(guó)科學(xué)院地理科學(xué)與資源研究所合作開發(fā)的中國(guó)大陸公里網(wǎng)格人口數(shù)據(jù)集， 該數(shù)據(jù)集的適用性已經(jīng)得到驗(yàn)證(陳振拓等， 2012； 丁文秀等， 2014； 安基文等， 2015)， 且數(shù)據(jù)具有良好的更新機(jī)制。

1.2.2 建立地震滑坡危險(xiǎn)性數(shù)據(jù)

本文利用“十一五”國(guó)家科技支撐課題“重大地震災(zāi)害及其災(zāi)害鏈綜合風(fēng)險(xiǎn)評(píng)估技術(shù)”產(chǎn)出的地震滑坡危險(xiǎn)性評(píng)估模塊繪制地震災(zāi)區(qū)滑坡危險(xiǎn)性分布預(yù)測(cè)圖。該模塊可給出不同烈度條件下的地震滑坡危險(xiǎn)性空間分布狀況， 其產(chǎn)出的數(shù)據(jù)集為90m×90m柵格數(shù)據(jù)， 根據(jù)地震滑坡概率值從小到大將滑坡危險(xiǎn)性劃分為5個(gè)等級(jí)， 危險(xiǎn)性等級(jí)從低到高的屬性值依次為1、 2、 3、 4、 5。屬性值從1到5， 表示地震滑坡危險(xiǎn)性等級(jí)從極低到極高(圖 4)， 屬性值越高， 發(fā)生滑坡的可能性越大(白仙富等， 2015)。圖 4 中， 每個(gè)柵格的屬性值即代表該空間位置的地震滑坡危險(xiǎn)性等級(jí)。以往的地震應(yīng)急和地震災(zāi)害風(fēng)險(xiǎn)預(yù)評(píng)估中的試驗(yàn)性應(yīng)用表明該套數(shù)據(jù)具有很好的合理性和實(shí)用性(Anetal.， 2015； 白仙富等， 2015； 高娜， 2015； 張如前， 2015； 魏連雨等， 2016； 張超， 2016； 邱丹丹等， 2017； 曹彥波等， 2018； 裴強(qiáng)等2018； 易樹健等， 2018； 鄧樹榮等， 2019； 楊根云等， 2019)。本文在利用該數(shù)據(jù)集生成研究區(qū)地震滑坡危險(xiǎn)性數(shù)據(jù)時(shí)， 使用現(xiàn)場(chǎng)調(diào)查后確定的烈度從數(shù)據(jù)集中提取對(duì)應(yīng)烈度區(qū)的地震滑坡危險(xiǎn)性等級(jí)數(shù)據(jù)， 然后通過GIS鑲嵌功能建立整個(gè)研究區(qū)地震滑坡危險(xiǎn)性空間分布數(shù)據(jù)(圖 5)。

圖 4 Ⅵ～Ⅺ度地震烈度設(shè)定下地震滑坡危險(xiǎn)性等級(jí)空間分布(90m柵格數(shù)據(jù))Fig. 4 Landslide susceptibility distribution assuming different seismic intensities， panels Ⅵ～Ⅺ show the landslide susceptibility according to intensity Ⅵ～Ⅺ(90m grid)， respectively.

圖 5 建立地震滑坡危險(xiǎn)性數(shù)據(jù)流程圖Fig. 5 A flow chart showing the susceptibility dataset of the earthquake-induced landslide.

1.2.3 公里網(wǎng)格地震滑坡屬性賦值

建立研究區(qū)的地震滑坡危險(xiǎn)性等級(jí)數(shù)據(jù)后， 在GIS中通過空間關(guān)聯(lián)統(tǒng)計(jì)功能模塊對(duì)每個(gè)公里網(wǎng)格單元內(nèi)的地震滑坡危險(xiǎn)性像元信息進(jìn)行統(tǒng)計(jì)， 統(tǒng)計(jì)對(duì)象為像元中心點(diǎn)完全落在公里網(wǎng)格單元內(nèi)的滑坡危險(xiǎn)性數(shù)據(jù)。 統(tǒng)計(jì)各個(gè)公里網(wǎng)格單元內(nèi)每類屬性值的滑坡危險(xiǎn)性像元數(shù)量， 并將統(tǒng)計(jì)結(jié)果作為屬性值賦值給對(duì)應(yīng)的公里網(wǎng)格作為其地震滑坡屬性(表2)。公里網(wǎng)格單元屬性示例表中， FID表示字段的機(jī)器自動(dòng)編碼， 若研究區(qū)共有N個(gè)公里網(wǎng)格單元， 編碼從0開始到N-1結(jié)束， 共N個(gè)字段編碼。P為人口字段名， 表中的PN表示第N個(gè)公里網(wǎng)格單元內(nèi)的人口數(shù)為PN人；X1表示屬性為1的地震滑坡危險(xiǎn)性像元字段名，X2等以此類推；X1，1表示第1個(gè)公里網(wǎng)格單元內(nèi)屬性值為1的滑坡危險(xiǎn)性像元個(gè)數(shù)為X1，1個(gè)，X1，2表示第1個(gè)公里網(wǎng)格單元內(nèi)屬性值為2的滑坡危險(xiǎn)性像元個(gè)數(shù)為X1，2個(gè)，X1，3等以此類推。在實(shí)際的屬性表中，X1，1—XN，5均為非負(fù)的整數(shù)型數(shù)值。

表2 公里網(wǎng)格屬性示例表Table2 A sample table of kilometer-grids’ properties

1.2.4 公里網(wǎng)格地震滑坡致死率函數(shù)

將公里網(wǎng)格單元的地震滑坡致死率與其內(nèi)部不同屬性值的滑坡危險(xiǎn)性像元數(shù)量之間的關(guān)系描述為函數(shù)f。 將滑坡危險(xiǎn)性的屬性看作地震滑坡致死率的影響因子， 根據(jù)影響因子的數(shù)據(jù)屬性， 函數(shù)f可以有不同的選擇， 例如：

Mi=f(Xi，1，Xi，2，Xi，3，Xi，4，Xi，5)

(2)

式中，Mi表示第i個(gè)公里網(wǎng)格單元的地震滑坡致死率，f表示推斷的影響因子與Mi的關(guān)系，Xi，1～Xi，5為地震滑坡致死率影響因子，Xi，1表示第i個(gè)公里網(wǎng)格單元內(nèi)地震滑坡危險(xiǎn)性屬性為1的像元數(shù)量，Xi，2表示第i個(gè)公里網(wǎng)格單元內(nèi)地震滑坡危險(xiǎn)性屬性為2的像元數(shù)量， 其余以此類推。理論上， 1個(gè)公里網(wǎng)格單元的地震滑坡致死率最低為0， 最高為1， 隨著地震滑坡危險(xiǎn)性屬性這一影響因子的變化而相應(yīng)地變化， 這種范圍增長(zhǎng)率的變化和影響因子之間的關(guān)系通常可以用logistic模型進(jìn)行描述， 則將式(2)表達(dá)為常用的logistic模型：

(3)

式中，Mi表示第i個(gè)公里網(wǎng)格單元的地震滑坡致死率， 取值范圍為0～1，k表示關(guān)系式函數(shù)的系數(shù)，a0是常數(shù)項(xiàng)，a1、a2、a3、a4、a5表示影響因子(滑坡危險(xiǎn)性屬性)的參數(shù)，Xi，1～Xi，5分別為第i個(gè)公里網(wǎng)格單元內(nèi)地震滑坡危險(xiǎn)性屬性值為1～5的像元數(shù)量。模型中的未知參數(shù)k、a0、a1、a2、a3、a4、a5可通過多次觀察樣本數(shù)據(jù)獲得。求解時(shí)， 把樣本數(shù)據(jù)中每條數(shù)據(jù)作為一次觀察數(shù)據(jù)， 一共觀察N次。假設(shè)每條樣本數(shù)據(jù)都滿足式(3)， 則：

(4)

式中，i表示第i次觀察， 共觀察N次。觀察的目標(biāo)就是利用數(shù)據(jù)和這些函數(shù)關(guān)系估計(jì)未知參數(shù)k、a0、a1、a2、a3、a4、a5。在統(tǒng)計(jì)學(xué)上， 把致死人數(shù)di看作隨機(jī)變量， 且具有如下分布規(guī)律：

(5)

N次觀察的似然函數(shù)為

(6)

為便于計(jì)算， 通常取對(duì)數(shù)， 因此有對(duì)數(shù)似然函數(shù)：

(7)

將式(4)代入式(7)可得：

(8)

該函數(shù)中，N、Pi、Di、Xi，1、Xi，2、Xi，3、Xi，4、Xi，5是已知的， 即為我們觀察到的樣本數(shù)據(jù)， 故該函數(shù)為未知參數(shù)k、a0、a1、a2、a3、a4、a5的函數(shù)。未知參數(shù)k、a0、a1、a2、a3、a4、a5的最大似然估計(jì)應(yīng)使得對(duì)數(shù)l(k，a0、a1、a2、a3、a4、a5)似然函數(shù)取得最大值。函數(shù)l(a0、a1、a2、a3、a4、a5)表達(dá)式的最后一項(xiàng)和k、a0、a1、a2、a3、a4、a5無(wú)關(guān)， 求解優(yōu)化問題時(shí)可以省略， 即只需求解如下函數(shù)的最大值點(diǎn)：

(9)

該函數(shù)優(yōu)化問題的解需采用數(shù)值模擬方法求得， 即通過樣本數(shù)據(jù)進(jìn)行模擬求取。

1.2.5 地震滑坡致死人數(shù)評(píng)估

整個(gè)區(qū)域的地震滑坡致死人數(shù)等于區(qū)域內(nèi)各公里網(wǎng)格單元致死人數(shù)的總和。根據(jù)1.1節(jié)所述的基本思想， 若將區(qū)域R劃分為N個(gè)1km×1km的網(wǎng)格單元， 總死亡人數(shù)D可以通過式(10)進(jìn)行評(píng)估：

(10)

式中，D表示區(qū)域R內(nèi)地震滑坡致死人數(shù)，N表示區(qū)域R內(nèi)共劃分為N個(gè)公里網(wǎng)格單元，Di表示區(qū)域內(nèi)第i個(gè)公里網(wǎng)格單元的地震滑坡致死人數(shù)。評(píng)估時(shí)通常需要對(duì)計(jì)算得到的直接結(jié)果D進(jìn)行取整處理。

2 震例案例

2.1 研究區(qū)域

本文選取2014年魯?shù)镸S6.5地震災(zāi)區(qū)、 2012年彝良MS5.7、MS5.6地震災(zāi)區(qū)和2008年汶川MS8.0地震災(zāi)區(qū)3個(gè)區(qū)域(圖 6)為研究區(qū)， 檢驗(yàn)所提出的基于GIS和logistic模型的地震滑坡致死人數(shù)快速評(píng)估方法。3個(gè)研究區(qū)是川滇乃至整個(gè)西南高原山地的典型區(qū)域， 對(duì)中國(guó)其他山區(qū)也有一定的代表性。這3個(gè)研究區(qū)中， 使用地震死亡人員樣本數(shù)據(jù)最詳實(shí)的魯?shù)闉?zāi)區(qū)的相關(guān)數(shù)據(jù)構(gòu)建地震滑坡致死人數(shù)快速評(píng)估方法并評(píng)價(jià)其有效性。彝良和汶川研究區(qū)數(shù)據(jù)則用來(lái)對(duì)構(gòu)建方法在更廣研究范圍內(nèi)開展地震滑坡致死人數(shù)評(píng)估的外延適用性進(jìn)行評(píng)價(jià)。

圖 6 案例研究區(qū)域Fig. 6 The case study areas.

2.2 建立公里網(wǎng)格人口數(shù)據(jù)

在汶川研究區(qū)， 用2007年的人口數(shù)據(jù)對(duì)研究區(qū)的公里網(wǎng)格單元進(jìn)行人口屬性賦值， 對(duì)彝良研究區(qū)和魯?shù)檠芯繀^(qū)使用2011年的人口數(shù)據(jù)對(duì)公里網(wǎng)格單元進(jìn)行人口屬性賦值， 得到研究區(qū)具有人口屬性值的公里網(wǎng)格空間數(shù)據(jù)(圖 7)。

圖 7 研究區(qū)公里網(wǎng)格人口數(shù)量的空間分布Fig. 7 The kilometer-grid population distribution in the study areas.

2.3 反演地震滑坡危險(xiǎn)性數(shù)據(jù)

按照1.2.2節(jié)所述， 建立研究區(qū)的地震滑坡危險(xiǎn)性空間分布數(shù)據(jù)(圖 8)， 并按照1.2.3節(jié)所述方法， 對(duì)研究區(qū)內(nèi)的公里網(wǎng)格單元進(jìn)行滑坡危險(xiǎn)性賦值， 最終建立具有人口、 滑坡屬性的公里網(wǎng)格單元數(shù)據(jù)， 數(shù)據(jù)屬性如表2 所示。

圖 8 研究區(qū)地震滑坡危險(xiǎn)性空間分布反演Fig. 8 Inversion of the distributions of the landslide susceptibility in the study areas(assessed results).

2.4 魯?shù)檠芯繀^(qū)模型的建立及其評(píng)價(jià)

2.4.1 地震滑坡致死率影響因子分析

2014年8月3日魯?shù)镸S6.5地震滑坡直接造成134人死亡和116人失蹤， 我們將失蹤人員等同于死亡人員， 把250人作為該次地震滑坡致死總?cè)藬?shù)。根據(jù)地震滑坡致死和失蹤人員的具體空間信息， 計(jì)算出研究區(qū)每個(gè)公里網(wǎng)格單元的地震滑坡致死人數(shù)和致死率。魯?shù)檠芯繀^(qū)共包括10395個(gè)公里網(wǎng)格單元， 故有10395條地震滑坡致死率信息， 限于篇幅， 此處僅列出部分結(jié)果(表3)。列出的數(shù)據(jù)中不包括地震滑坡致死人數(shù)為0的網(wǎng)格單元信息， 但在整個(gè)研究區(qū)， 絕大多數(shù)網(wǎng)格單元內(nèi)實(shí)際地震滑坡致死人數(shù)均為0， 且這些滑坡致死人數(shù)為0的網(wǎng)格單元中大部分為無(wú)人區(qū)。有滑坡致死的網(wǎng)格單元主要分布在烈度為Ⅸ度的區(qū)域。屬性表中，P表示公里網(wǎng)格單元的人口數(shù)字段名；D表示公里網(wǎng)格單元的地震滑坡致死人員數(shù)量字段名；M表示公里網(wǎng)格單元的地震滑坡人員致死率字段名， 由各網(wǎng)格單元對(duì)應(yīng)的D/P得到；X1～X5分別表示公里網(wǎng)格單元內(nèi)地震滑坡危險(xiǎn)性屬性為1～5的字段名， 表中X1對(duì)應(yīng)的數(shù)字是該公里網(wǎng)格單元內(nèi)屬性為1的地震滑坡危險(xiǎn)性像元數(shù)量，X2對(duì)應(yīng)的數(shù)字是屬性為2的地震滑坡危險(xiǎn)性像元數(shù)量， 以此類推。

表3 魯?shù)檠芯繀^(qū)公里網(wǎng)格單元屬性表(部分)Table3 Property of the kilometer-grids in Ludian study area

為了判明公里網(wǎng)格單元地震滑坡致死率與網(wǎng)格內(nèi)每類地震滑坡危險(xiǎn)性像元數(shù)量之間有無(wú)關(guān)聯(lián)， 即是否獨(dú)立， 采用列聯(lián)表(contingency table)對(duì)魯?shù)檠芯繀^(qū)10395個(gè)公里網(wǎng)格單元的地震滑坡致死率與地震滑坡危險(xiǎn)性屬性的相合性進(jìn)行分析與推斷。我們把滑坡死亡率分為“>0”和“=0” 2類， 分別考慮不同地震滑坡危險(xiǎn)性屬性值與死亡率分類之間的相合關(guān)系。通過建立滑坡危險(xiǎn)性屬性像元數(shù)量與死亡率的列聯(lián)表(表4)， 對(duì)X1～X5與M分別進(jìn)行列聯(lián)分析。分析方式為： 原假設(shè)H0，X1～X5與M相互獨(dú)立； 備擇假設(shè)H1，X1～X5與M正相合， 檢驗(yàn)數(shù)據(jù)支持哪個(gè)命題。X5與M的列聯(lián)分析檢驗(yàn)p值為1.773×10-93， 相合系數(shù)為0.2359， 當(dāng)網(wǎng)格單元的X5取值=0時(shí)， 傾向于M的取值=0；X5的取值>0時(shí)， 傾向于M的取值>0。以上結(jié)果表明X5與M相互獨(dú)立的原假設(shè)不成立， 兩者存在正相合關(guān)系。X4與M的列聯(lián)分析檢驗(yàn)p值為4.67×10-184， 相合系數(shù)為0.3718， 當(dāng)網(wǎng)格單元的X4取值=0時(shí)， 傾向于M的取值=0；X4的取值>0時(shí)， 傾向于M的取值>0。以上結(jié)果表明X4與M相互獨(dú)立的原假設(shè)不成立， 兩者存在正相合關(guān)系。X3與M的列聯(lián)分析檢驗(yàn)p值為8.854×10-34， 相合系數(shù)為0.1389， 當(dāng)網(wǎng)格單元的X3取值=0時(shí)， 傾向于M的取值=0；X3的取值>0時(shí)， 傾向于M的取值>0。以上結(jié)果表明X3與M相互獨(dú)立的原假設(shè)不成立， 兩者存在正相合關(guān)系。X2與M的列聯(lián)分析檢驗(yàn)p值為7.139×10-83， 相合系數(shù)為-0.03289， 當(dāng)網(wǎng)格單元的X2取值=0時(shí)， 傾向于M的取值>0；X2的取值>0時(shí)， 傾向于M的取值=0。以上結(jié)果表明X2與M相互獨(dú)立的原假設(shè)不成立， 兩者存在負(fù)相合關(guān)系。X1與M的列聯(lián)分析檢驗(yàn)p值為9.26×10-106， 相合系數(shù)為-0.03371， 當(dāng)網(wǎng)格單元的X1取值=0時(shí)， 傾向于M的取值>0；X1的取值>0時(shí)， 傾向于M的取值=0。以上結(jié)果表明X1與M相互獨(dú)立的原假設(shè)不成立， 兩者存在負(fù)相合關(guān)系。

表4 X與M列聯(lián)表Table4 Contingency table about X and M

從列聯(lián)表分析結(jié)果看， 每組的檢驗(yàn)p值都遠(yuǎn)<0.001， 表明地震滑坡致死率與影響因子相互獨(dú)立的原假設(shè)不成立， 即公里網(wǎng)格單元的地震滑坡致死率與地震滑坡危險(xiǎn)性屬性相關(guān)。從相合系數(shù)看，X5、X4、X3與M有顯著的正相合關(guān)系， 網(wǎng)格單元內(nèi)屬性值為5、 4、 3的地震滑坡危險(xiǎn)性像元增加時(shí)， 公里網(wǎng)格單元的地震滑坡致死率傾向于>0。X1、X2則與M的增長(zhǎng)有顯著的負(fù)相合關(guān)系， 網(wǎng)格單元內(nèi)屬性值為1和2的地震滑坡危險(xiǎn)性像元增加時(shí)， 公里網(wǎng)格單元的地震滑坡致死率傾向于=0。

2.4.2 地震滑坡致死人數(shù)評(píng)估模型

從魯?shù)檠芯繀^(qū)的10395條樣本數(shù)據(jù)中選取表3 內(nèi)的44條數(shù)據(jù)外再隨機(jī)選擇3001條沒有地震滑坡人員死亡的數(shù)據(jù)作為訓(xùn)練集， 進(jìn)行基于公里網(wǎng)格單元的地震滑坡致死人數(shù)logistic建模， 并用其余的7350條數(shù)據(jù)作為測(cè)試集對(duì)所建立模型進(jìn)行合理性檢驗(yàn)。從比例來(lái)看， 訓(xùn)練集數(shù)據(jù)約占樣本總數(shù)的30%， 測(cè)試集數(shù)據(jù)約占樣本總數(shù)的70%。我們通過R軟件的Newton方法對(duì)式(9)中的最大值點(diǎn)和未知參數(shù)k、a0、a1、a2、a3、a4、a5進(jìn)行求解。將軟件輸出結(jié)果代入式(3)， 得到地震滑坡致死率與地震滑坡危險(xiǎn)性屬性的logistic回歸關(guān)系：

(11)

式(11)即為針對(duì)公里網(wǎng)格單元的地震滑坡致死率計(jì)算模型。從回歸的公里網(wǎng)格單元地震滑坡logistic致死率模型看，X5和X1的系數(shù)絕對(duì)值較大， 表明地震滑坡危險(xiǎn)性屬性最高和最低的像元數(shù)量變化對(duì)公里網(wǎng)格單元地震滑坡致死率影響較大； 其中，X5的系數(shù)最大， 表明公里網(wǎng)格單元內(nèi)滑坡危險(xiǎn)性屬性為5的像元數(shù)量增加對(duì)網(wǎng)格單元地震滑坡致死率增加的貢獻(xiàn)最大，X1的系數(shù)最小， 表明公里網(wǎng)格單元內(nèi)滑坡危險(xiǎn)性屬性為1的像元數(shù)量增加對(duì)網(wǎng)格單元地震滑坡致死率趨于0的貢獻(xiàn)最大。死亡率模型中，X5、X4、X3的系數(shù)分別為0.04077、 0.03130、 0.01365， 說明公里網(wǎng)格單元內(nèi)屬性值為5、 4、 3的滑坡危險(xiǎn)性像元數(shù)量增加則網(wǎng)格單元的地震滑坡致死率增大， 屬性值越大的像元數(shù)量對(duì)網(wǎng)格單元致死率增大的貢獻(xiàn)越大；X2、X1系數(shù)分別為-0.01666、 -0.09652， 說明公里網(wǎng)格單元內(nèi)屬性值為2和1的滑坡危險(xiǎn)性像元數(shù)量增加則網(wǎng)格單元的地震滑坡致死率減小， 屬性值越小的像元數(shù)量對(duì)網(wǎng)格單元致死率減小的貢獻(xiàn)越大。為了評(píng)估該模型的有效性， 我們對(duì)訓(xùn)練集和測(cè)試集分別進(jìn)行了F測(cè)試值計(jì)算。F值越大， 則模型擬合度越好， 模型解釋數(shù)據(jù)的能力越強(qiáng)， 理論上， 只有當(dāng)F≥0.5時(shí)模型才通過檢驗(yàn)。在R軟件中， 通過>1-pchisq(hp.logit$deviance， hp.logit$df.residual)語(yǔ)句對(duì)訓(xùn)練集和測(cè)試集的F值進(jìn)行計(jì)算， 輸出結(jié)果均為[1]1， 即F=1。這表明， 地震滑坡致死率是地震滑坡人員致死性的良好指標(biāo)， 同時(shí)也進(jìn)一步表明， 所建立的模型可用于公里網(wǎng)格單元的地震滑坡致死率計(jì)算。

結(jié)合式(10)和式(11)可推導(dǎo)出基于公里網(wǎng)格單元的地震滑坡致死人數(shù)評(píng)估模型， 簡(jiǎn)化后的模型表達(dá)式為

(12)

用式(12)對(duì)魯?shù)檠芯繀^(qū)訓(xùn)練集的地震滑坡致死人數(shù)進(jìn)行計(jì)算， 訓(xùn)練集中實(shí)際因地震滑坡致死的人數(shù)為135人， 模型檢測(cè)結(jié)果為129人， 模型計(jì)算的地震滑坡致死人數(shù)的精度為95.55%； 對(duì)魯?shù)檠芯繀^(qū)測(cè)試集的地震滑坡致死人數(shù)進(jìn)行計(jì)算， 測(cè)試集中實(shí)際因地震滑坡致死的人數(shù)為115人， 模型計(jì)算結(jié)果為104人， 模型計(jì)算的地震滑坡致死人數(shù)的精度為90.43%。對(duì)整個(gè)研究區(qū)而言， 魯?shù)榈卣鹬幸蚧露劳?含失蹤)的災(zāi)民共有250人， 模型檢測(cè)結(jié)果為233人， 誤差率為6.80%。

我們通過計(jì)算魯?shù)檠芯繀^(qū)各公里網(wǎng)格單元實(shí)際的地震滑坡致死人數(shù)與模型評(píng)估結(jié)果的Kappa系數(shù)評(píng)估所建立模型的統(tǒng)計(jì)學(xué)意義。計(jì)算網(wǎng)格單元的Kappa值時(shí)， 將每個(gè)網(wǎng)格的總?cè)藬?shù)拆分為2部分， 即地震滑坡致死人數(shù)和其他人數(shù)。魯?shù)榈卣馂?zāi)區(qū)的Kappa值為0.912， 表明模型計(jì)算結(jié)果與實(shí)際幾乎完全一致， 地震滑坡致死人數(shù)計(jì)算模型具有較好的統(tǒng)計(jì)學(xué)意義。從評(píng)估死亡人數(shù)空間分布(圖9a)和實(shí)際地震滑坡致死人數(shù)的空間分布(圖6a)來(lái)看， 兩者也具有較高的空間分布一致性。

圖 9 研究區(qū)公里網(wǎng)格地震滑坡致死人數(shù)空間分布反演Fig. 9 Inversion of earthquake-induced landslides casualties using GIS and the logistic method.

2.5 模型在汶川和彝良研究區(qū)的應(yīng)用

以汶川和彝良地震災(zāi)區(qū)為研究區(qū)域， 對(duì)該方法的外延性及其在更大區(qū)域范圍的適用性進(jìn)行檢驗(yàn)。將2.4.2節(jié)中建立的地震滑坡致死人數(shù)評(píng)估模型應(yīng)用到彝良和汶川研究區(qū)時(shí)， 對(duì)評(píng)估的方法流程、 模型的參數(shù)等均不做任何改動(dòng)。汶川地震滑坡致死人數(shù)的模型評(píng)估結(jié)果為18732人， 與汶川地震時(shí)因滑坡死亡和失蹤約20000人的統(tǒng)計(jì)結(jié)果(Chenetal.， 2011)相比約少千余人， 誤差率約為6.5%。此外， 汶川地震公里網(wǎng)格單元的地震滑坡致死人數(shù)的空間分布(圖9b)與滑坡分布比較一致。由于缺乏汶川地震研究區(qū)具體人員死亡地點(diǎn)信息， 故不再計(jì)算地震滑坡致死率的F值及各公里網(wǎng)格單元地震滑坡致死人數(shù)的Kappa檢驗(yàn)值。彝良研究區(qū)地震滑坡致死率的F檢驗(yàn)值為0.893， 雖然比魯?shù)檠芯繀^(qū)的檢驗(yàn)值低， 但其結(jié)果仍能滿足統(tǒng)計(jì)學(xué)要求， 模型能夠通過數(shù)學(xué)檢驗(yàn)。彝良研究區(qū)實(shí)際因地震滑坡致死59人(白仙富等， 2014)， 用本文的地震滑坡致死人數(shù)評(píng)估模型計(jì)算結(jié)果為死亡48人， 誤差率為18.64%。通過對(duì)比可知， 各公里網(wǎng)格單元死亡人數(shù)的評(píng)估結(jié)果(圖9a)與實(shí)際地震滑坡致死人數(shù)分布(圖6a)仍比較一致， 彝良研究區(qū)各公里網(wǎng)格單元地震滑坡致死人數(shù)的Kappa檢驗(yàn)值為0.889。彝良地震中因滑坡而死亡的人員有半數(shù)以上為外地戶口在礦區(qū)的務(wù)工者(白仙富等， 2014)， 而這些外地戶口的人員數(shù)量只有部分被統(tǒng)計(jì)到對(duì)應(yīng)的公里網(wǎng)格單元中， 這是彝良地震滑坡致死人數(shù)相比另外2個(gè)研究區(qū)的評(píng)估誤差率較大的原因之一。

從基于GIS和logistic模型的地震滑坡致死人數(shù)評(píng)估方法在汶川和彝良研究區(qū)的測(cè)試結(jié)果看， 本文建立的地震滑坡致死人數(shù)快速評(píng)估模型在允許一定誤差的前提下， 可以用于估計(jì)其他區(qū)域因地震滑坡災(zāi)害導(dǎo)致的人員死亡數(shù)量， 具有較好的外延適用性。

模型估計(jì)的3個(gè)研究區(qū)的地震滑坡致死人數(shù)都比實(shí)際地震滑坡致死人數(shù)偏低， 造成這一現(xiàn)象的原因與建模樣本數(shù)據(jù)和選用的方法都有一定關(guān)系。用于建模的魯?shù)榈卣馂?zāi)區(qū)樣本數(shù)據(jù)中， 因地震滑坡有致死人數(shù)的網(wǎng)格單元數(shù)量遠(yuǎn)遠(yuǎn)小于致死人數(shù)為0的網(wǎng)格單元數(shù)量。建模時(shí)， 在算法上不僅要考慮的有人員死亡的樣本的準(zhǔn)確性， 還必須考慮沒有人員死亡的樣本的正確率。從另一個(gè)角度看， 模型對(duì)沒有人員死亡樣本的預(yù)測(cè)的正確率遠(yuǎn)高于對(duì)有死亡樣本預(yù)測(cè)的正確率。如果減少?zèng)]有地震滑坡致死人數(shù)的樣本數(shù)據(jù)進(jìn)行建模， 則可能導(dǎo)致死亡人數(shù)的評(píng)估結(jié)果超過實(shí)際致死數(shù)量， 從而帶來(lái)另一種誤差， 這種誤差是所基于的統(tǒng)計(jì)方法對(duì)于數(shù)據(jù)驅(qū)動(dòng)固有的不足帶來(lái)的。在調(diào)試魯?shù)榈卣馂?zāi)區(qū)模型的過程中發(fā)現(xiàn)， 選取無(wú)死亡人數(shù)區(qū)域和有死亡人數(shù)區(qū)域的數(shù)據(jù)分別為3001條和44條樣本進(jìn)行建模是總體誤差最小的一種搭配。從模型評(píng)估結(jié)果的具體數(shù)據(jù)看， 對(duì)3個(gè)研究區(qū)的評(píng)估結(jié)果在數(shù)量級(jí)上都沒有偏差， 與實(shí)際地震滑坡致死人數(shù)也比較接近， 能滿足當(dāng)前應(yīng)急評(píng)估等工作的大部分需求。

3 討論與結(jié)論

顯然， 公里網(wǎng)格單元的地震滑坡致死率與90m地震滑坡危險(xiǎn)性數(shù)據(jù)屬性之間的相關(guān)性是本文建立地震滑坡災(zāi)害致死人數(shù)判斷模型的關(guān)鍵。如前文所述， 本文模型的成立需要2個(gè)條件： 1)地震滑坡致死人數(shù)評(píng)估對(duì)象的尺度是1km×1km的網(wǎng)格單元； 2)所使用的地震滑坡危險(xiǎn)性數(shù)據(jù)是90m×90m分辨率的分等級(jí)數(shù)據(jù)。后者又包括2個(gè)條件： 1)地震滑坡危險(xiǎn)性數(shù)據(jù)是以烈度為地震動(dòng)影響參數(shù)的產(chǎn)品數(shù)據(jù)； 2)地震滑坡危險(xiǎn)性等級(jí)從低到高分為5個(gè)等級(jí)， 相應(yīng)的像元的屬性值為1～5。當(dāng)?shù)卣鸹轮滤廊藬?shù)評(píng)估的基本單元尺度發(fā)生變化時(shí)， 本文提出的這種方法的效果可能會(huì)出現(xiàn)變化； 同樣， 如果使用的地震滑坡危險(xiǎn)性數(shù)據(jù)不是危險(xiǎn)性分級(jí)數(shù)據(jù)或不是90m分辨率的數(shù)據(jù)， 本文所建立模型的有效性也同樣會(huì)發(fā)生變化。

在地震滑坡致死人數(shù)評(píng)估方法建立中， 我們考慮了針對(duì)公里網(wǎng)格單元的地震滑坡致死率與地震滑坡危險(xiǎn)性屬性之間的關(guān)系。這種關(guān)系可用不同模型描述， 本文使用了logistic模型建立地震滑坡致死率與地震滑坡危險(xiǎn)性屬性之間的關(guān)系， 不同類型的模型將帶來(lái)不同的評(píng)估效果， 如何開發(fā)不同的模型并篩選出更有效的模型是一個(gè)亟待解決的問題。

通過對(duì)魯?shù)椤?汶川和彝良3個(gè)研究區(qū)進(jìn)行的方法構(gòu)建與應(yīng)用的初步嘗試表明， 本文提出的針對(duì)公里網(wǎng)格單元的地震滑坡致死人數(shù)快速評(píng)估方法是有效的， 且在其他多山區(qū)域的地震損失評(píng)估中也可能具有一定的應(yīng)用前景。同時(shí)， 該方法本質(zhì)上是一種統(tǒng)計(jì)方法， 當(dāng)樣本數(shù)量越多或研究區(qū)范圍越廣時(shí)， 該方法評(píng)估結(jié)果的誤差可能越?。?相反， 如果研究區(qū)的地震震級(jí)越小或范圍越狹窄， 由于人員在空間分布上固有的動(dòng)態(tài)變化等偶然性因素的影響， 用該方法評(píng)估的地震滑坡致死人數(shù)的相對(duì)誤差可能存在更大的離散性。

本文提出的地震滑坡致死人數(shù)評(píng)估方法是對(duì)現(xiàn)有地震地質(zhì)災(zāi)害致死人數(shù)評(píng)估方法存在缺失現(xiàn)狀的一種探索性研究， 嘗試在地震后快速地定量化評(píng)估地震滑坡導(dǎo)致的死亡人數(shù)。在歷史地震滑坡致死數(shù)據(jù)的豐富性、 基礎(chǔ)數(shù)據(jù)的現(xiàn)勢(shì)性、 模型方法的多樣性等方面， 都存在一定的不足之處， 需要在今后進(jìn)一步發(fā)展和完善。

致謝審稿人為本文提出了寶貴的意見和建議； 云南師范大學(xué)史正濤教授為本工作提供了支持； 昭通市防震減災(zāi)局賀莉晶老師在野外人員死亡地點(diǎn)調(diào)查工作中予以協(xié)助， 李春博士為本文工作提供了幫助。在此一并表示感謝!