刘 佳，赵海军，马凤山，郭 捷，孙琪皓，段学良

(1.中国科学院地质与地球物理研究所/中国科学院页岩气与地质工程重点实验室，北京 100029；2.中国科学院地球科学研究院，北京 100029；3.中国科学院大学，北京 100049)

0 引言

青藏公路是世界上海拔最高、线路最长的公路，自建成后便加快了内陆与高原地区的物资运输速度，促进了西藏经济建设和人民生活改善。公路沿线泥石流灾害是损毁路基、路面及桥梁、隧道等构造物的重要外在动力，也是造成公路长时间中断的主要原因之一[1]。所以，利用数学方法和计算机技术，建立更高效的评价模型进行危险性分区，对泥石流的预测和防治至关重要。

国内外对泥石流危险性评价的研究自20世纪伊始。HOLLINGSWORTH等[2]采用因子叠加法和打分方法构建滑坡泥石流危险度评价的基本框架，实现了滑坡危险度的评估；MANTOVANI等[3]基于RS技术对欧洲滑坡灾害进行研究和分区；CHUNG等[4]针对调查区的滑坡灾害情况，采用GIS的专业模型扩展分析功能进行了多指标综合区划；YALCIN[5]利用层次分析法和二元统计法研究了滑坡灾害易发性；PRADHAN等[6]利用神经网络方法结合逻辑回归模型对马来西亚槟岛地区滑坡灾害进行危险性分区；YANG等[7]对比了两种集对分析方法在泥石流危险性评价中的应用；殷坤龙等[8]应用GIS调查研究后发现滑坡灾害发育具有区域性和重复性的规律；黄润秋等[9]利用GIS技术和数学模型方法进行地灾区划；LIU[10]提出了泥石流灾害等级与发生频率之间指数关系的理论解；张春山等[11]研究了区域地质灾害风险评价要素权值的计算方法；薛东剑等[12]借助遥感数据解译和地面调查以及GIS技术，运用信息量模型和逻辑回归模型完成了对青川、平武县的地质灾害危险性评价。REN等[13]建立了SPA-GRA模型并利用了变异系数法来确定矿山参数；由此看来，随着科技进步、社会发展，GIS等计算机科学技术和多种数学方法逐渐被应用于泥石流危险性评价。

尽管评价方法不断增加，但目前泥石流评价仍存在问题[14]：如评价指标赋权过程中，层次分析法主观性过强，客观变异系数法具有一定局限性、未考虑主观因素的影响等等。本文以国道G109沿线拉萨至那曲地区路段为例，通过利用层次分析法来改进变异系数法，弥补了原方法适用性差的缺陷，并利用改进的变异系数法对研究区进行了泥石流危险性分区，从而确定危险区域。

1 评价模型构建

1.1 评价指标的选取

评价指标是整个评价体系最基本的组成元素。泥石流灾害的形成与其所处的地貌形态、物质、构造、水文、气象和植被条件等多种因素密切相关[15]，在充分考虑调查区灾害特点的基础上，并着重考虑泥石流灾害形成所需要的三个基本条件[16]，选取地层岩性、地形地貌、地质构造(活动断裂)、气候气象与地震五类因子作为公路沿线泥石流危险性评价的指标因子，以此作为评价体系中的准则层，选取依据如下：

(1)区内主要分布有花岗岩、砂岩、页岩及第四系堆积物等，易风化的岩浆岩和较软弱岩层在温度变化剧烈、日照强烈的青藏高原，更易于风化、崩解堆积于泥石流沟两岸，往往能为泥石流提供更加丰富的物质基础。

(2)地形条件的影响主要体现在重力势能以及能量转化梯度方面，而重力势能的高低和能量转化率实质上受坡体地形坡度和起伏度(泥石流坡体物源区至沟口的垂直距离)的影响：坡度不仅影响斜坡内的应力分布，同时对斜坡上松散物质的堆积厚度起着决定性作用，进而决定着斜坡的稳定性[17]。相对高差大和地表切割强烈地区的岩土体剥蚀作用强烈，水流速度相对较快，具备了泥石流的赋存条件[18]。

(3)由于欧亚板块与印度板块的不断挤压，青藏高原构造活动剧烈。断裂构造越发育则岩体中次生结构面越多、岩体越破碎。构造作用主要体现在断层等宏观表现上，主要根据采样区内各类构造线的密度分布状况[19]，对所发育的泥石流沟的固体物储备方量、补给形式以及泥石流的性质有很大影响。地质构造的复杂程度将直接影响到岩层的稳定性[20]；进而为泥石流发育提供丰富的物源基础。

(4)气候也是泥石流产生的重要影响因素，它不仅影响风化作用，还会因为暴雨、连续降雨诱发泥石流滑坡等地质灾害。降雨过程中，孔隙水压力不断升高，直接影响滑裂面以下土体向泥石流转化的速度[21]。

(5)地震烈度越高(地震加速度越大)则岩土体越易遭受破坏，为泥石流发育提供丰富物质条件；地震通过为泥石流提供固体物质来源、水源、动力条件和激发条件等影响泥石流的形成和发展[22]。

1.2 评价指标体系的建立

根据层次分析法的相关理论，评价指标体系一般包括目标层、准则层和指标层。目标层是指泥石流危险性的评价；准则层即评价危险性的指标因子，选取原则如上所述；指标层是对准则层的细化，综合考虑调查区内外部因素，将准则层细分为六个指标：岩石易风化程度、平均坡度、起伏度、断层线密度(单位面积内断裂线的长度)、地震加速度、汛期(6～9月)月均降水量。对三个评价分析层次分别用字母A、B、C表示，最终建立危险性评价指标体系(图1)。

1.3 改进变异系数法的提出

泥石流危险性评价因子权重的确定有主观和客观两种方法。对于主观方法来说：传统层次分析法在计算权重系数时，相对权的最小平方法复杂一些，但其所需的一致性检验步骤保证了判断矩阵的合理性和一致性。使用层次分析法主要是为了定性定量地研究所赋权重问题，具有试探性，所构造判断矩阵具有很大主观性，若是专家经验不足或者对所构造矩阵意见不一就会导致准确度较低甚至错误的结论。对于客观变异系数法来说：此方法具有一定局限性，仅考虑数据的变异离散程度而确定权重系数，这点对于不可量化的评价指标因子不适用。所以本文对变异系数法灵活创新运用，利用层次分析法来修正可量化评价的指标因子所占权重，从而得到综合考虑主观、客观两方面因素的评价模型。

1.3.1层次分析法

本文将采用两种层次分析法，即传统层次分析法和权的最小平方法，要使上述两种方法所得的权重值之间的差异程度与其相对应的分配系数间的差异程度相一致,引入距离函数的概念[23]，得到相应比例系数从而取得最终层次分析法所得权重值。

(1)传统层次分析法

此方法大致可以分为四个步骤(详细计算过程及一致性检验方法参见文献[24])：在这一过程中根据专业人员建议，考虑西藏地区特殊的地质环境条件，构造判断矩阵见表1和表2。

表1 指标判断矩阵(1)

表2 指标判断矩阵(2)

(2)权的最小平方法

权的最小平方法是由A.T.W.Chu等于1979年提出，于1988年介绍到我国[25]。权的最小平方法的基本原理是在一定的约束条件下，建立目标函数，对其一致性判断矩阵A求得最小值，即一致性判断矩阵排序向量的最优解[26]。由权的最小平方法构成的最优化问题，用向量形式可表示为[27]：

(1)

②约束条件：

(2)

通过求解可得最优解为：

(3)

式中：W*——一致性判断矩阵排序权向量。

(4)

e=(1,1,1,…,1)T

(5)

式中：aij——判断矩阵中的数值；

wi或wj——权重系数。

1.3.2对变异系数法进行改进

变异系数又称为标准差率，它是利用求取一组数据的标准差与平均值之比，来作为评价体系指标权重的客观赋权方法。该方法可以直接利用ArcGIS软件中关于指标因子的栅格数据文件的平均值和标准差，通过简单计算就可以得到各权重系数，具体计算公式如下：

(6)

(7)

式中：S——某组指标因子数据的标准差；

R——平均值；

CVi——某组指标因子的变异系数；

Mi——指标因子所占的客观权重。

此处，首先选取传统层次分析法和权的最小平方法进行组合赋权，各自的权重系数依据距离函数所得，目的是减少误差、增大可信度。

W12=αW1+βW2

(8)

α=0.511，β=0.489

(9)

式中：W12——聚合赋权法所确定的权重；

W1——传统层次分析法所确定的权重；

W2——权的最小平方法所得结果。

其次，利用组合赋权后的权重系数，对变异系数法进行改进。对于不可量化的评价指标，权重系数取层次分析法所得结果，这一过程可以避免变异系数法不能衡量不可量化指标权重的缺点；对于可量化评价指标，本着总权重和为1的原则，权重系数为由层次分析法所得的总权重系数与变异系数法所得的各个评价指标权重系数的乘积：

(10)

式中：Vi——评价指标改良后最终权重系数；

W12i——利用上述组合赋权法确定的权重系数；

Mi——变异系数法确定的指标权重。

1.4 综合评价模型

确定各个评价指标因子的权重，仅仅是完成了该指标因子在致灾评价中作用高低的评价，而要想获得最终的综合评价结果以及建立综合风险性分区，需要利用各个指标因子的权重与分级分数的乘积确定。本文采用加权综合评分模型进行综合评价，公式如下：

(11)

式中：Q——泥石流危险性系数；

xi——第i个指标因子在某一评价单元的分级赋值分数；

vi——该指标因子的权重系数。

2 研究区域危险性分析

2.1 研究区概况

研究区域位于西藏自治区拉萨市至那曲市的国道G109(毗邻青藏铁路)附近10 km范围内，地处青藏高原东南部，北邻唐古拉山脉，南靠喜马拉雅山脉，地势总体北高南低，西高东低；构造活动剧烈，多高山河谷地貌；区域内包括西藏拉萨市和那曲市的众多重要地区：拉萨、当雄、那曲、安多和唐古拉山口(其为西藏自治区和青海省的分界)。气候属于高原气候，平均气温低，昼夜温差大。据已有调查资料显示：那曲地区已发生地质灾害110多处，造成几千亿元经济损失。

2.2 计算过程

2.2.1数据来源

本文研究数据主要涉及泥石流评价指标因子的数据。其中，地层岩性和断层线密度数据主要来自于1∶50 万地质图资料；地形坡度和起伏度的DEM数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)，借鉴前人泥石流危险性评价经验[28-29]采用30 m精度(考虑到数据获取的难易程度和最大利用效果)；汛期月均降雨量数据来源于全国降雨统计资料；地震加速度数据来自于《中国地震动峰值加速度区划图》(GB 18306—2015)。将这些数据导入ArcGIS软件并划分网格，共有11 917 个面积为1 km×1 km的网格，形成栅格数据文件，并将其分为五级。

2.2.2权重系数的计算

根据前面章节内容提到的评价指标权重的确定方法，代入数据整理计算后可得结果。传统层次分析法(W1)和权的最小平方法(W2)、组合后的层次分析法(W12)、变异系数法(M)以及改进变异系数法(V)的计算结果见表3。

表3 各种方法所得权重

图2 评价指标分级图Fig.2 Classifying picture of evaluation index

2.2.3确定评价指标分级分数

泥石流的形成离不开特定的自然地理、地质和气候条件：利用ArcGIS软件，根据上述所选评价指标因子，对西藏地区拉萨至那曲国道G109附近10 km范围(前人研究主张距离铁路5 km以上的泥石流均不会对青藏铁路交通干线产生破坏性影响[30]，为了综合评价更广泛的区域灾害特点，本文决定扩大为10 km范围)内的区域进行指标因子等级划分(表4)，结合“《中国公路综合自然区划》体系框架研究”和文献[31]的研究成果，将六类评价指标分别划分为五个等级(图2)；为便于更直观地了解等级划分结果，对其中部分指标数据采取适当取整处理，分别赋值1、3、5、7、9；与大部分研究中都选择15作为分级数值相比，此方法使得分级分数间拉开差距，更具明显性(g为地震加速度，取值10 m/s2)。

2.3 结果验证分析

利用ArcGIS数据处理功能，得到研究区域泥石流灾害各个指标因子的分级图，在此基础上，利用软件的加权叠加功能，形成泥石流综合评价危险性分区图(图3)。借助其数据统计功能，对改进变异系数法所建模型进行分析：泥石流危险度中高区域(中等、较高、高)面积约2.3×103km2，占调查区总面积的20%，主要分布在：国道G109拉萨至羊八井路段约1 000 km2范围内；羊八井—当雄路段道路右侧部分区域；当雄县东北方向约30 km处的道路左侧拐角处约170 km2范围；以及安多县周围部分区域。泥石流危险度低区域(较低、低)约占总调查面积的80%：大部分分布在那曲县城北部。将实际灾害数量分布特点(图3)和危险性分区图对照，结合灾害点数量在各个等级危险区的分布情况(图4)，综合来看评价可靠度较高。根据现场调查结果及四川省广汉地质工程勘察院所撰写的地灾调查与区划报告资料得出，23处(约占泥石流总量的79%)泥石流位于评价模型中高危险区，其余6处(约占泥石流总量的21%)处于低、较低危险区。根据现场调查情况分析，图3中高危险区往往处于高山峡谷段，地势较陡，而低危险区往往地势平坦，这就验证了改进变异系数方法具有较好的准确度和合理性。

表4 泥石流评价指标等级划分

图3 危险性分布图Fig.3 The distribution map of hazards

图4 灾害数量分布图Fig.4 The quantity distribution map of disaster

3 结论

本文针对西藏拉萨至那曲地区的国道G109附近10 km范围内的泥石流危险性评价问题，利用改进变异系数法借助ArcGIS的数据统计和加权叠加等功能，构建了高效合理的泥石流危险性评价模型，并进行泥石流危险性分区及评价，所得主要结论包括：

(1)提出了一种改进变异系数法，该方法克服了层次分析法的主观性及客观变异系数法的局限性，并将其应用于泥石流危险性综合评价模型，以此确定的中高危险区与泥石流高发区重合度较大，此方法可靠性较高。

(2)研究区泥石流中高危险区域面积约2.3×103km2(仅占调查区总面积的20%)，泥石流发生数量占总量的79%，国道G109东侧比西侧危险性更高，拉萨至当雄路段尤其明显；泥石流低危险区约9.6×103km2，占总调查面积的80%，但灾害数量仅占21%。

(3)变异系数法在地质灾害危险度研究中还不是很广泛，利用改进变异系数法，可克服原方法只能计算可量化指标的缺点，具有更广泛的应用空间。