非完整滑坡编目三参数反Gamma概率分布模型

胡元鑫¹，刘新荣¹，蒋洋²，葛华³

 (1. 重庆大学 土木工程学院，重庆，400045；

2. 长安大学 公路学院，陕西 西安，710064；

3. 成都地质调查中心，四川 成都， 610081)

摘要：中国“5.12”地震震级高，震源浅，于震区内触发了大量滑坡。为了研究地震滑坡的概率分布特征，对四川映秀极震区(约109.3 km²)地震滑坡进行系统编目，该编目属非完整滑坡编目。利用该编目建立三参数反Gamma概率分布模型，其拟合参数分别为b=1.11，a=2.17×10^{-3 km2}及s=-3.925×10^{-4 km2}。基于该模型对编目区内滑坡面积均值、总面积及最大滑坡面积进行计算，并引入滑坡标度及建立指数修正函数分别对“5.12”四川汶川地震触发滑坡总数及滑坡总面积进行预测，同时建议模型建立原则。研究结果表明：模型的幂律衰减参数值与已有研究结果相近；模型对编目区内滑坡面积均值、总面积及最大滑坡面积的计算结果与实测值之间的相对误差均为10%左右；模型对“5.12”汶川地震触发滑坡总数及滑坡总面积的预测结果均与已有研究结果比较接近。

关键词：滑坡编目；Gamma函数；概率分布；“5.12”地震；四川映秀极震区

中图分类号：P642          文献标志码：A         文章编号：1672-7207(2011)10-3176-06

Three-parameter inverse-gamma probability distribution model of incomplete landslide inventory

HU Yuan-xin¹, LIU Xin-rong¹, JIANG Yang², GE Hua³

(1. College of Civil Engineering, Chongqing University, Chongqing 400045, China;

2. School of Highway, Chang’an University, Xi’an 710064, China;

3. Chengdu Center of China Geological Survey, Chengdu 610081, China)

Abstract: Because the “5.12” earthquake is of high magnitude and shallow hypocenter, it triggered a large number of landslides in seismic area. In order to research probability distribution of these earthquake-induced landslides, systematic earthquake-induced landslides inventory in Yingxiu meizoseismal area in Sichuan Province (about 109.3 km²), which is an imcomplete landslide inventory compared with those throughout seismic area, was finished. On the basis of the incomplete landslides inventory, a three-parameters inverse-gamma probability distribution model was presented, in which the three parameters are b=1.11, a=2.17×10^{-3 km2}, s=-3.925×10^{-4 km2} respectively. By means of the model, the mean, total and maximum landslide areas in inventoried area were calculated, and the total number and area of landslide triggered by the “5.12” earthquake were predicted by introducing landslide magnitude scale and exponential correcting function. Meanwhile, some suggestions for building three-parameters inverse-gamma probability distribution model were proposed. The results show that power-law decay parameters of the model are comparable to the previous research results. Computed relaive errors, with respect to the mean, total and maximum landslide areas in inventoried zone, are all about 10%. The total number and area of “5.12” earthquake induced landslide resulting from prediction of the model are both comparable to the previous research results.

Key words: landslide inventory; Gamma function; probability distribution; “5.12” earthquake; Yingxiu meizoseismal area in Sichuan Province

系统的滑坡编目可提供滑坡数量、类型、位置、规模等有效数据^[1-3]。结合地形、地质单元等环境因  子^[4-5]与地震、降雨、融雪等触发因素^[6-9]，滑坡编目能反映滑坡的分布特征^[10-11]、滑坡几何特征与频率关系^[12]及与各触发因子之间的关系^{[2, 5-6]}。由于滑坡为复杂地质体，含有大量非确定性信息^[13]，因此，对于大区域的滑坡常采用统计分析方法进行研究^[2]。滑坡概率密度函数多具有幂指数关系^{[2, 10, 12]}。在单触发因素触发滑坡编目和历史滑坡编目的统计分析方面，常用的滑坡面积-频率模型有累积或非累积频率统计、Double Pareto概率分布，Den等^[10]对此进行了详细列举和分析。Malamud等^[2]基于不同地点并分别由地震、融雪及降雨触发的三类滑坡编目数据引入Γ函数对滑坡的频率、面积、体积等进行研究，其成果获得了广泛应用^{[2, 10, 14]}。上述研究成果均体现了具较大规模滑坡概率密度分布的幂律衰减形式，可用于确定滑坡面积与规模的期望概率。但Malamud等^[2]所建立的三参数反Gamma概率分布模型是基于相对完整滑坡编目，当应用于非完整滑坡编目时有较大差别^{[2, 10]}。为快速对滑坡灾情进行评估，本文以2008年“5.12”四川映秀极震区地震滑坡编目为数据源建立三参数反Gamma概率分布模型，从滑坡面积入手探讨滑坡规模与频率的关系，并对“5.12”地震触发滑坡总体进行预测。其结果可用于震区及区域级滑坡灾害风险    评价。

1  模型数据源

滑坡编目区坐标范围为31°1′33″N—31°6′58″N和103°23′12″E—103°29′58″E，其面积约109.3 km²，属极震区。“5.12”地震于区内触发了大量滑坡(图1)。地震滑坡编目数据源于震后的高分辨率卫星影像(分辨率为2.5 m)、航片(比例尺为1:3 000)解译及现场调查结果。编目区共有1 948处地震滑坡，滑坡总面积为40.74 km²，为编目区面积的37.2%，其中最小的滑坡面积约为80 m²，最大的滑坡面积约为1.7 km²。现场调查结果表明：编目区内地震滑坡绝大部分为浅层分裂型滑坡(滑体平均厚度＜5 m)，滑体体积相差极大，最小的岩坠为数立方米，最大的岩石崩泻(莲花新沟和牛圈沟滑坡)超过1 000万m³。

2  模型的建立

2.1  滑坡概率密度计算

地震滑坡概率密度与滑坡面积之间的关系可通过计算其概率密度函数获得^[2]：

式中：N_LT为编目区内地震滑坡总数；δN_L为处于滑坡面积A_L~A_L+δA_L之间的地震滑坡数量；δA_L为地震滑坡面积增量(km²)。

图1  编目区内的地震滑坡

Fig.1  Earthquake-induced landslides in inventory area

利用区内地震滑坡编目数据，地震滑坡概率密度函数与滑坡面积之间的相关性见图2。

图2  地震滑坡概率密度p与面积A_L之间的相关性

Fig.2  Dependence of earthquake-induced landslide probability density p on landslide area A_L

2.2  三参数反Gamma概率分布模型

基于式计算的地震滑坡概率密度可按Malamud等提出的函数形式^[2]进行拟合，其拟合曲线见图2中粗实线，判定系数R²＝0.875 6。拟合模型   如下：

式中：b=1.11；a=2.17×10^{-3 km2}；s=-3.925×10^{-4 km2}。其方差分析见表1。

表1　 p(A_L;b,a,s)的方差分析

Table 1　 ANOVA for p(A_L;b,a,s)

式表明中等面积与大面积滑坡具有典型的幂律衰减形式，这与文献[2,10,12,14]中的结论一致；而小面积滑坡则具有指数翻转形式，其最大概率分布位于6.36×10^{-4 km2}处。式中参数b控制中等面积与大面积滑坡的幂律衰减幅度，-(b+1)=-β=-2.11，与文献[10]中的统计值(β=2.3)类似；参数a主要影响最大概率密度出现的滑坡面积范围；参数s则主要控制小面积滑坡的指数翻转模式。

在正常情况下，式在s≤A_L＜+∞时对A_L的积分应当为1，但本文所定义的概率分布模型仅定义 0≤A_L＜+∞。由于该式中参数s为负，故有：

根据不完全Γ函数定义，有：

令t=a/(A_L-s)，对t求微分并改变式积分限，则：

因此，式变为：

式右端为修正项。代入式中各参数值，则该修正项值为0.004，可忽略不计。因此，可认为式即为映秀极震区地震滑坡的三参数反Gamma概率分布模型。显然，对于具有大面积值的滑坡，式可简化为：

3  讨论

3.1  编目区滑坡面积的均值、总和及最大值

理论上，滑坡面积均值是其概率密度函数的一阶矩，即：

将式代入式，有：

由式及忽略修正项，则式为：

令t=a/(A_L-s)，则有：

由于-a/s=5.53，exp(-5.53)≈0，故式可近似为：

再利用Γ函数的定义及性质，式改写为：

将式各参数代入式，得即映秀极震区地震滑坡面积均值为0.019 km²。而编目的实际滑坡面积均值为0.021 km²(表2)，与式计算值差别不大。

由于式中N_LT仅为编目区内的地震滑坡数量，因此，编目区内的地震滑坡总面积可用下式获得：

式中：A_LT为滑坡的总面积(km²)。则图1中的地震滑坡总面积为0.019×1 948=37.01 km²，与编目区内实测滑坡总面积(40.74 km²)也非常接近。

令N_LC(≥A_L)为面积不小于A_L的滑坡累积数，据式和，有：；当N_LC=1时，有：。代入式中参数，则A_Lmax=1.90 km²，与实际编目最大滑坡面积1.70 km²较吻合，见表2。

3.2  非完全滑坡编目的外推预测

式系映秀极震区地震滑坡编目的拟合结果，尽管该式在编目区内具有较高精度，但该编目相对“5.12”地震触发滑坡的分布空间而言属非完全滑坡编目。因此，应对式的外推预测能力加以分析，其预测指标包括地震触发滑坡总量及其面积。为此，引入滑坡标度的概念^[2]，即：。式中：m_L为滑坡标度，系某事件触发滑坡数量的表征；N_LT为某事件触发滑坡总数。据式(14)，则有：

     (15)

根据式，则滑坡数与面积的频率密度函数为：

            (16)

表2  模型的计算值与实测值

Table 2  Computed value of model and measured value

结合式与式(16)，有：

         (17)

式(17)表明频率密度可通过式与N_LT相乘得到，并可与m_L产生联系。当m_L分别为1，2，3，4和5时，f(A_L)的理论计算曲线见图3。

图3表明：映秀滑坡编目的幂律衰减尾端滑坡标度m_L=4.0±0.1，即“5.12”地震触发滑坡总数N_LT为10 270± 2 322。该数据与已有研究结果的比较见表2。本文所得N_LT与文献[7，15-16]中的值较接近，但与文献[16]中的值有较大差别。文献[7，15-18]在统计面积上略有不同，但均属汶川“5.12”地震的重灾区，并且文献[7,15,16]均采用现场调查进行校核，因此，其滑坡总数具有较高可信度。文献[16-17]中的滑坡总数远高于文献[7，15-16]中及本文预测数，究其原因主要在于以下2个方面：(1) 震区地震滑坡中岩坠比例较大，由于震区地形陡峻，岩坠易在山壁形成弹跳翻滚，故留下多处破坏痕迹，这是造成遥感解译中滑坡数增大的原因之一；(2) 现场调查发现，震区地震滑坡常呈连片状分布，不同物源区形成1处滑坡或1处物源区形成多处滑坡，现场调查常将这类情形视滑坡机理计为1处滑坡复合体，而按文献[17-18]，则可能计为多处滑坡。同时，考虑到文献[7，15-18]中统计面积的差别，本文模型对“5.12”地震触发滑坡总数的预测是可信的。

图3  滑坡频率密度f与滑坡面积A_L的关系

Fig.3  Relationship between landslide frequency density and area

因图1所示滑坡编目为非完整编目，式(15)不能用于直接计算完整编目的滑坡总面积，但基于式及式(17)中的指数衰减形式，可在式(15)的基础上利用文献[17-18]中的编目与滑坡面积增加1个指数修正函数，对假定完整编目的滑坡总面积进行预测，即：

    (18)

       (19)

式中：。

综合式(18)和(19)并代入式，有：

  (20)

式中：r为非完整编目面积与待编目面积的比值，其中非完整编目面积为109.3 km²，震区平均滑坡密度按3处/km^-2计^[15]，则根据图3确定的m_L所表征的滑坡总量应分布于不小于4 200 km²的区域，因此，式(20)的定义域上限取为109.3/4 200＝0.025。由于式(18)中m_L为本文编目区内“5.12”地震触发滑坡标度，因此，该式仅适用于在本文已有非完整编目基础上计算待编目滑坡面积，即r越小，式(20)基于预测的待编目面积就越大。式(20)表示的A_LT与r的关系曲线见图4。

如将文献[17]中的编目面积48 678 km²视为完整编目面积，则式(20)的计算值为701.5 km²，与实测值711.8 km²的相对误差仅为1.45%。若将文献[17]中的编目面积视为完整编目面积，则式(20)的计算值为684.98 km²，与实测值687.09 km²也非常接近。同时，当r→0时，即滑坡编目面积趋向无穷小时，式(20)趋近于某一定值，具有明显的物理意义，即“5.12”地震触发滑坡总面积不会随编目面积增大而无限增大。

图4  基于编目面积A_LT与待编目面积比r的滑坡总面积A_LT预测

Fig.4  Predicted total landslide area based on ratio of inventoried area to expected inventory area

3.3  三参数反Gamma概率分布模型的建立原则

以上分析表明：三参数反Gamma概率分布模型可有效应用于滑坡编目研究。但Malamud等^[2]通过位于不同区域并分别由地震、强降雨及融雪等因素触发滑坡建立的三参数反Gamma概率分布模型极易给人以错觉，即认为所有的滑坡编目都具有相同的概率分布特征，如具有相同的滑坡面积均值、相同的最大概率分布点等。事实上，由于不同地质环境、不同因素的触发效应(不同的震级、降雨强度)不同，滑坡数量、规模、分布区域等均可能不同。另外，对于分布于编目区域内的老滑坡，其触发因素、机理及分布特征也可能不同。因此，式中参数与文献[2]中的参数稍有差别，这种差异系客观原因引起，但仍处于同一量级水平。除上述地质环境、触发因素等客观原因外，这种差异还源于以下原因：(1) 编目区内分布有大量小滑坡，这些小滑坡使数据离散性大并造成其概率密度发生波动；(2) 映秀地区岷江两侧及渔子溪北侧的滑坡常呈连片状分布，常出现不同物源、不同主滑方向的滑坡有相互重叠的区域，形成滑坡复合体，这给滑坡周界辨识和编目带来一定误差。为此，本文建议利用已有滑坡编目建立三参数反Gamma概率分布模型时应遵循以下2条原则以保证滑坡编目精度：(1) 对于不同因素触发的滑坡应分别建立模型，这对研究滑坡机理与滑坡概率分布的关系有重要意义；(2) 明确区分编目区内老滑坡和新近触发滑坡并分别编目及建立模型。据上述2条原则，必须建立完善的滑坡编目数据库并进行多时相不同因素触发滑坡编目信息的更新，以获知滑坡与触发因素不同强度(如震级、降雨强度)之间的关系。

4  结论

(1) 本文提出的三参数反Gamma概率分布模型基于“5.12”四川映秀极震区非完整地震滑坡编目数据建立的，其参数分别为b=1.11，a=2.17×10^{-3 km2}及s=-3.925×10^{-4 km2}，其中模型的幂律衰减参数与已有研究成果相近，滑坡的最大概率分布位于6.36×10^{-4 km2}处。

(2) 模型得到的编目区滑坡面积均值、总面积及最大滑坡面积分别为0.019，37.01和1.9 km²，均与编目区内实测值较接近，其相对误差分别为9.5%，9.1%和11.7%，说明模型在编目区内具有较高精度。

(3) 模型基于滑坡标度与滑坡面积频率分布对“5.12”地震触发滑坡总数的外推预测数据与已有研究数据均在同一数量级内，具有较高可信度；模型基于本文非完整滑坡编目与指数修正函数对假定完整编目区内滑坡总面积的预测结果与已有研究成果非常接近。上述预测数据均说明本文所建立的非完整滑坡编目三参数反Gamma概率分布模型具有实用意义。

参考文献：

[1] Corominas J, Moya J. A review of assessing landslide frequency for hazard zoning purposes[J]. Engineering Geology, 2008, 102(3/4): 193-213.

[2] Malamud B D, Turcotte D L, Guzzetti F, et al. Landslide inventories and their statistical properties[J]. Earth Surface Processes and Landforms, 2004, 29(6): 687-711.

[3] Dai F C, Lee C F, Ngai Y Y. Landslide risk assessment and management: an overview[J]. Engineering Geology, 2002, 64(1): 65-87.

[4] Chau K T, Sze Y L, Fung M K, et al. Landslide hazard analysis for Hong Kong using landslide inventory and GIS[J]. Computers & Geosciences, 2004, 30(4): 429-443.

[5] 王秋良, 薛宏交, 杨强, 等. 北川地区地震灾害的地形因子分析[J]. 大地测量与地球动力学, 2008, 28(6): 68-72.
WANG Qiu-liang, XUE Hong-jiao, YANG Qiang, et al. Analysis on terrain factors of earthquake disaster in Beichuan area[J]. Journal of Geodesy and Geodynamics, 2008, 28(6): 68-72.

[6] Keefer D K. Statistical analysis of an earthquake-induced landslide distribution: The 1989 Loma Prieta, California event[J]. Engineering Geology, 2000, 58(3/4): 231-249.

[7] 黄润秋, 李为乐. “5.12”汶川大地震触发地质灾害的发育分布规律研究[J]. 岩石力学与工程学报, 2008, 27(12): 2585-2592.
HUANG Run-qiu, LI Wei-le. Research on development and distribution rules of geohazards induced by Wenchuan earthquake on 12th May, 2008[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(12): 2585-2592.

[8] Lee L M, Gofar N, Rahardjo H. A simple model for preliminary evaluation of rainfall-induced slope instability[J]. Engineering Geology, 2009, 108(3/4): 272-285.

[9] Salciarini D, Godt J W, Savage W Z, et al. Modeling landslide recurrence in Seattle, Washington, USA[J]. Engineering Geology, 2008, 102(3/4): 227-237.

[10] van Den E M, Poesen J, Govers G, et al. Characteristics of the size distribution of recent and historical landslides in a populated hilly region[J]. Earth and Planetary Science Letters, 2007, 256(3/4): 588-603.

[11] Jaiswal P, Van Westen C J. Estimating temporal probability for landslide initiation along transportation routes based on rainfall thresholds[J]. Geomorphology, 2009, 112(1/2): 96-105.

[12] Stark C P, Hovius N. The characterization of landslide size distributions[J]. Geophysical Research Letters, 2001, 28(6): 1091.

[13] Refice A, Capolongo D. Probabilistic modeling of uncertainties in earthquake-induced landslide hazard assessment[J]. Computers & Geosciences, 2002, 28(6): 735-749.

[14] Guzzetti F, Reichenbach P, Cardinali M, et al. Probabilistic landslide hazard assessment at the basin scale[J]. Geomorphology, 2005, 72(1/4): 272-299.

[15] 黄润秋, 李为乐. 汶川地震触发崩塌滑坡数量及其密度特征分析[J]. 地质灾害与环境保护, 2009, 20(3): 1-7.
HUANG Run-qiu, LI Wei-le. Analysis on the number and density of landslides triggered by the 2008 Wenchuan earthquanke, China[J]. Journal of Geological Hazards and Environment Preservation, 2009, 20(3): 1-7.

[16] 殷跃平. 汶川八级地震地质灾害研究[J]. 工程地质学报, 2008, 16(4): 433-444.
YIN Yue-ping. Researches on the geo-hazards triggered by Wenchuan Earthquake, Sichuan[J]. Journal of Engineering Geology, 2008, 16(4): 433-444.

[17] 许冲, 戴福初, 姚鑫. 汶川地震诱发滑坡灾害的数量与面积[J]. 科技导报, 2009, 27(11): 79-81.
XU Chong, DAI Fu-chu, YAO Xin. Incidence number and affected area of wenchuan earthquake-induced landslides[J]. Science & Technology Review, 2009, 27(11): 79-81.

[18] 许冲, 戴福初, 陈剑, 等. 汶川Ms8.0地震重灾区次生地质灾害遥感精细解译[J]. 遥感学报, 2009, 13(4): 754-762.
XU Chong, DAI Fu-chu, CHEN Jian, et al. Identification and analysis of secondary geological hazards triggered by a magnitude 8.0 Wenchuan Earthquake[J]. Journal of Remote Sensing, 2009, 13(4): 754-762.

(编辑 陈灿华)

收稿日期：2010-10-05；修回日期：2010-12-28

基金项目：中国地质调查局项目(1212010914011)

通信作者：刘新荣(1969-)，男，浙江开化人，教授，从事工程地质研究；电话：023-65120727；E-mail：liuxrong@126.com