用多重罗吉斯回归和gis技术预测美国堪萨斯东北的滑坡灾害.doc

过去的几十年间，沿堪萨斯东北的堪萨斯河和密苏里河多山地区的滑坡造成了上百万美元的财产损失。为了解决这一问题，用称为多重罗吉斯回归的统计方法对艾奇逊堪萨斯州、堪萨斯州及其周边地区进行了滑坡灾害制图。包括数字化地质、斜坡和滑坡等资料，用ArcView GIS进行处理。罗吉斯回归将地理单元内数据预测变量和滑坡是否发生联系起来，用这种联系编制的图可以显示给定斜坡和地质条件的区域将来发生滑坡的概率。结果表明在研究区进行滑坡评价时，斜坡是最重要的变量。由页岩、粉砂岩和砂岩组成的地质单元最易发生滑坡。考虑了土壤类型和滑坡长宽比因素，但在最后的分析中这些因素被排除了，原因是这些变量对罗吉斯回归预测功效作用不明显。土壤类型与地质单元具有显著相关性，滑坡和斜坡形状之间没有明显联系。
一、概述
现有两种常用的滑坡预测方法。一种是基于控制边坡失稳物理机制的数学模型确定或工程方法，确定性方法的主要局限性是在大区域内难以获得特性数据。第二种方法是统计学方法，这样不用关心控制边坡失稳的机制，但是假定了过去发生的滑坡与地形的可测量特征有关。反过来，这些特征可以用来预测将来滑坡的发生。由于在研究区不易获得必要的特性数据，因而选择统计学方法来编制美国艾奇逊堪萨斯州滑坡灾害图。
1995年，在堪萨斯Overland公园的一处滑坡，毁坏了两所房屋，每所价值都超过40万美元，这一破坏突出了在堪萨斯城市地区编制滑坡灾害图的必要性。堪萨斯地质调查局目前正在调查堪萨斯州部分大城市地区的滑坡。研究区包括202.3km2大部分的农耕地和艾奇逊城区，总人口约1万。
研究区具有滑坡发生的必要条件，包括陡坡、充足的降雨、松软土层和岩石地层。尽管堪萨斯在轻微起伏不平的美国中部是相对较为平坦的地区，但在整个州都存在有陡坡的多山地区。艾奇逊地区的地形主要是由于密苏里河的切割作用形成，支流进入宾夕法尼亚基岩，由砂岩、粉砂岩、页岩和石灰岩组成。第四纪沉积物上覆于基岩上，由冰碛、黄土和冲积物组成。基岩风化壳和上覆第四纪沉积物生成富含粘粒的土壤，易于发生滑坡。
研究区包括发达区和不发达区，具备天然和人为的滑坡诱发因素。天然诱发因素包括大暴雨、潮湿的气候和河流侵蚀；人为因素包括已被改造的陡峭斜坡、不恰当的施工或填埋以及透水性差等。

为用量化方法表达滑坡发生的可能性，填图必须采用概率的概念，这是一个滑坡发生相对频率的评价指标。总的来说，当在一定的时间间隔，也许是10年或25年，特定的位置或小区域范围，包括人口或结构脆弱性评价，进行滑坡发生概率评价时要用到“风险”。判断滑坡发生概率需要有以前滑坡发生时间的信息，而这一信息通常不易获得。相反，可以通过计算过去滑坡数量，不管其发生时间来判断概率。成果图指的是“灾害图”，表示滑坡发生的概率，并不表示发生率。概率图表示概率（低，中等，高）的顺序分类，而不是明确发生的概率可以根据主观概率评价获得。
在艾奇逊研究区，如果已知给定点的斜坡和地质情况，编制的滑坡灾害图，可以表示将来滑坡发生的可能性。使用者包括规划师、工程师、开发者、土地所有者，都可以用这一图件来评价灾害程度以采取合适的行动。本文包括选择输入数据、运用GIS技术处理数据、用罗吉斯回归生成滑坡灾害图。
二、艾奇逊研究区滑坡和区域地质概况
艾奇逊位于堪萨斯州东北65km处。密苏里河西部分布在堪萨斯州，西部是平缓起伏山地，东部是多岩石的山地，沿密苏里河及其支流陡岸分布。研究区西部的冰川沉积坡度较缓，坡角为从近水平山顶和洪泛区向峭壁变化。
（一）滑坡
“滑坡”定义为大量的土壤和岩体向下运动。本研究区的滑坡类型包括土崩、土体滑动、岩崩和危崖。
研究区的滑坡按滑坡特征分为近期滑坡和老滑坡。“近期”滑坡明显的特征是有确定的边界，最近几年发生，包括活动的、暂停的和休眠的土崩和土体滑动。老滑坡具有丘状地形、平缓特征和不明确的边界。这一类型包括休眠的、残余的和古老的土崩和土体滑动。
本文用近期滑坡和老滑坡资料讨论分析滑坡灾害。滑坡主要发生在土壤和风化壳浅层。由于较深的土质和岩质滑坡也常发生，因而浅滑坡常叠加于深部滑坡上。
（二）区域地质

研究区的地质组成为冰碛物、黄土和覆盖于宾夕法尼亚纪基岩上的第四纪冲积物。宾夕法尼亚纪岩石地层由滨海砂岩、粉砂岩、页岩和灰岩组成，重复序列或旋回反应了宾夕法尼亚时代的水位变化。正常旋回中的层序始于大陆砂岩和粉砂岩，通过浅海页岩和灰岩海进层序延续后进入海退层序。
研究区较老的地层单元是Lawrence地层，其组成为8m多厚的绿灰页岩，红粉砂质页岩和煤。劳伦斯地层沿艾奇逊南部山区分布，认为是堪萨斯东北一个最易发生滑坡的因素，许多近期和老滑坡都发生在Lawrence地层。
Oread灰岩覆盖在Lawrence地层上，由16m厚的灰岩和灰黑页岩组成。在研究区有四种灰岩，最大厚度是6m。在研究区Oread灰岩组成最陡的斜坡。Oread灰岩中2m厚的底部灰岩上是4m厚的灰海相页岩。海相页岩风化成富含粘粒的土壤；在一些高速公路切割面上已经观察到小的土崩和土体滑动。
Kanwaka页岩上覆于Oread灰岩之上，由18m厚的灰页岩、粉砂岩和砂岩组成，在地层中间为1m厚的灰岩床。在Kunwaka页岩上观察到了滑坡，但是主要是发生在1m厚的灰岩中。Lecompton灰岩上覆于Kanwaka页岩，由12m厚的灰岩和灰页岩组成，但在研究区出露的不太好。Tecumseh页岩上覆于Lecompton灰岩上，由约20m厚的灰页岩，粉砂岩和砂岩组成。在出露在高速公路切割面和陡坡的Tecumseh页岩中观测到了滑坡。Deer Greek 石灰岩上覆在Tecumseh页岩上，由11m厚的石灰岩和灰页岩组成。Calhoun 页岩上覆于Deer Creek灰岩上，由5m厚的粉砂岩和页岩组成。
将地质单元分为五类。根据碎屑含量（粘粒、粉砂和砂）将基岩分为两类，“页岩”类包括Lawrence地层，Kanwaka页岩，Tecumseh页岩和Calhoun页岩；“石灰岩”类包括Oread石灰岩,Lecompton石灰岩，Deer Creek 石灰岩和Topeka石灰岩。其余的三类是冰碛物，黄土和冲积物。
三、用GIS编制滑坡灾害图
本项研究用地理信息系统ArcView3.2编译和处理数据，并编制最终滑坡图。数据包括斜坡、坡形、地质和土壤。土壤类型和滑坡发生之间具有明显联系，然而土壤类型与地质单元之间显著相关，当将地质状况包括在内时土壤对判断滑坡发生的可能性作用不大。因为土壤是冗余的，在分析中将其排除在外。

尽管在一些情况斜坡面的方向可能会影响斜坡的稳定性，本研究区的坡向和滑坡之间没有统计学上的重要联系，因此将变量坡向也排除在外。
根据美国地质调查局的数字高程模型（DEMs）获得斜坡数据。DEMs在值之间有30m的间隔，而且数据质量很好。根据DEMs采用ArcView GIS的斜坡程序计算斜坡，计算了每一单元和八个相邻单元间的最大斜坡。采用10m栅格间距进行概率分析，故每个30m的单元分成10m栅格的9个相等单元。
在本项研究中，采用了一幅改编的地质图，其中附加了一些野外数据，包括早期一些地质图中的数据。将改编的地质图建造成数字多边形输入到ArcView GIS中，将多边形转化为10m栅格间距的光栅格式。将每个栅格单元分配一个代码，代表单元中心的地质单位。初始栅格化的地质图代表了研究区内的岩层，分配为5种如前所述的地质单元。
开始准备好的滑坡清单图作为数字矢量文件也用ArcView GIS转化为10m的格栅文件。如果没有滑坡，每个单元都分配为“0”；若有滑坡，将单元分配为“1”；如果单元不在研究区内分配为“无数据”代码。将地质、斜坡和滑坡格栅文件在逻辑上相比，保证三个格栅覆盖在一个共同区域，并转化成一个包括每一单元中心UTM坐标的ASCⅡ格式文件，之后将ASCⅡ文件输入到JMP v.4统计软件进行概率分析。
四、概率方法和多重罗吉斯回归逼近
大多数滑坡预测统计研究采用了关联表分析、多元回归、判别分析或罗吉斯回归。所有这些方法与测量滑坡发生的一个或多个预测变量（斜坡和地质状况等）相关。最简单的统计预测方法采用关联表分析，在两个输出状态（滑坡和无滑坡）通过预测变量的离散分类进行交叉造表。表中单元的计数比例可以作为估计预测变量的滑坡（或无滑坡）发生条件概率。
为采用多元回归，必须将滑坡表示成连续变量。由于回归估计不受限制，结果可能是负值或是大于100%，因此不能直接解译为概率。
采用差别分析，将观测结果分为互斥的两组：它们是（1）发生滑坡、（2）不发生滑坡的一个位置或区域。每个观测结果由独立变量的矢量组成。滑坡组和不发生滑坡组的区别由设计每个观测将两组的重心连接成一直线的系数表示；如果得分距离滑坡重心比离不发生滑坡的重心更近，观测结果就划分为滑坡。使人
遗憾的是，没有其它方法将判别得分作为滑坡发生的概率。实际上，判别分析仅是简单地将多元问题转化为单变量关联表分析。
罗吉斯回归，在一个区域内如10m单元滑坡的状态由变量Y表示。任意代码如“1”可以用来表示滑坡，如果没有滑坡代码，值为“0”，但是这些值没有实际意义。不是预测任意代码，而是将观测结果划分为某一类而与另一类相对来判断概率似乎更有用。由于有无滑坡分类是互斥的也是全部的，如果已知一个概率，另一个也就知道了；也就是说，P（Y＝1）＝[1－P（Y＝0）]，其中P（Y＝1）是圆括号内值为真时的概率。
Y＝1的概率可以通过回归模型化，采用等式如P（Y＝1）＝α＋βX。然而此时存在一个问题，尽管观测概率P（Y＝1）必须处于0和1之间，预测概率不受任何限制,它们可能小于0或大于1。为解决这一问题，用单数来表示概率。单数定义为某事发生概率与不发生概率之间的比值，或单数（Y＝1）＝P（Y＝1）/［1－P（Y＝1）］。单数是比值，没有固定最大值，但是有最小值0。如果通过回归将单数模型化，判断不可能最大值的概率就消除了，但是预测奇数可能是小于0的不可能值仍然存在。
然而，如果将奇数通过取自然对数，可以获得一个从负无穷大到正无穷大的变量，该值为分对数，定义为（Y）＝{lnP（Y＝1）/［1－P（Y＝1）］}，当奇数从0变为1时该值变为负无穷大，当奇数超过1时变为正无穷大。由于分对数可以假定任何值，通过回归判断的值可能超过最大或最小值的问题就避免了。
必须强调，概率、奇数和分对数是同一事物的三种表达形式，而且相互之间可以很容易转化。然而采用分对数，没有在预测模型中不能使用回归的限制因素。
在罗吉斯回归中，必须用最大似然法估计回归参数。为此有必要了解观测数据的概率分布，特别是给定预测变量集合下独立变量的条件分布，用参数表达。如果在特定点滑坡发生率i是1，滑坡未发生率是0，那么Yi服从二项式分布，均值为πi，自由度为πi(1－πi)。为将这一分布转化为α和β系数的回归函数，也必须假定πi和预测变量之间的联系服从罗吉斯分布函数（完全等于通过线性函数将均值Y与预测值相联系的假设）。实际上，罗吉斯分布函数对于滑坡发生的条件概率πi是一个很好的模型。
通过为α和β找到估计值将结果条件分布或似然函数最大化进行罗吉斯回归，观测值为y1,y2,……yn。这样可以观测到最可能值的系数。作为一般回归不
能直接判断系数，相反开始要选择一个试验方法并稍加修改，来看一下似然是否增加了。重复进行这一程序直到从这一步到下一步的似然函数增加可以忽略不计。
似然函数对数的最大化等于似然函数本身的最大化，具有一定的优点。如果对数似然乘以－2，得出的统计值称为负对数似然，近似于分布。仅包含常数项的罗吉斯回归模型负对数似然可以与包括所有预测变量参数估计量的负对数似然进行比较；两个模型之间的差别可以用于测试回归系数的显著性。如果区别不大，预测变量和发生滑坡的条件关系是不显著的。也就是说，考虑预测变量的数量不能提高判断滑坡发生的可能性。这一方法的变化可以用于测试单个预测变量并选择那些最有效的变量。测试用一种与在一般回归中用于测试系数的变量表相似的分析格式，除了测试统计服从分布而不是F分布以外。
一般的罗吉斯回归用将预测变量进行分类、连续变量和二者相结合的方法进行滑坡敏感性编图。目前的研究，采用多元罗吉斯回归得出灾害图，罗吉斯回归基于两个预测因子：斜坡（连续变量）和地质单元（分类变量）。
对于一个单元（滑坡和未发生滑坡）的两个可能响应级，最大似然回归模型是： 

可以改写为

其中Yi是单元i的状态，Xi是单元i预测变量的矢量，b是估计系数的矢量。等式右边是分对数变换，也就是说奇数的对数。
发生滑坡的概率预测值在0和1之间。由于在该区预测滑坡的主要变量是地形斜坡和地质状况，多元罗吉斯回归等式可以写为：

其中地质k是属性变量，如果有地质单元，k等于1，否则为0。注意分类互斥，也就是说在单元内只能有一种地质单位。这表明回归的总形式由常量加上，对于每一地质单元，回归曲线通过系数进行偏移量计算校正，与单元格中的地质单位有关。由于在无地质单位的单元内所有属性变量都为0，乘积。在单元内地质单位为1属性变量为1，因此。
艾奇逊研究区包括2022861个单元，变量表示发生或不发生滑坡，单元中的滑坡值和地质单位代码。通过与负对数似然函数相比检验拟合的罗吉斯回归的显著性，对只有常数的可供选择模式进行完全回归，结果是显著相关。
表1给出了模型合适的参数，用值检验个体项的显著性。常数项和斜坡系数都显著相关，无穷小概率值为0。地质状况作为滑坡的预测因子进行检验，也是相关显著，尽管“冲积层”类型的系数不太显著但也不是0。地质单位的系数约束为0，因此最后一类（页岩）必须等于其它系数总和的负值。由于预先确定了最后一个系数，没有定义其标准偏差。相似地，除了页岩对于单个类别，其值作为地质状况值和总值的差。
表1　斜坡和五种地质单位滑坡概率的多元罗吉斯回归系数
变量
系数
标准偏差
值
显著性
交叉
8.8117
0.0332
70470
＞0.9999
斜坡
-0.1472
0.0009
24561
＞0.9999
地质状况

2257.04
＞0.9999
冲积层
0.0815
0.0632
1.66
0.8030
黄土
0.5043
0.0476
112.43
0.9998
冰碛
1.0818
0.0428
639.16
0.9998
石灰岩
-0.6261
0.0278
507.01
0.9998
页岩
-1.0415*
**
996.78*
0.9998
*由差值确定
**不能判断
必须牢记回归模式假设存在严重干扰，显著性水平受到影响。这些干扰最严重的影响是观测是独立的。由于地表上物体到处都是连续的，地形单元必须与相邻单元密切相连。应当根据适当的地球统计进行罗吉斯分析，但是还没有这样的分析工具。
根据这一点，可以计算出不同斜坡和地质单元的合适分对数回归模式，生成一系列曲线来预测研究区滑坡发生的概率。曲线形式相同，形态取决于常数项和斜坡项，但是对给定的地质单位有所不同（见图1）。按照惯例，回归的不确定性由适合回归线周围的置信区间表现出来，置信区间的间隔解译为包含特定百分比（通常为95%）的真实回归。置信区间是估计回归系数和特定概率指标标准偏
差的函数。由于值太小，在图1中没有表现出罗吉斯回归曲线周围的置信区间，最大的不确定性是冲积层，60º斜坡95%置信间隔的概率仅为45.4～51.7%的。

五、研究结果
将每一格栅单元的地质和斜坡资料输入到罗吉斯方程，生成格栅，并采用ArcView GIS建造滑坡灾害图。给出某一单元的给定斜坡和地质状况，最终的滑坡灾害栅格包括该滑坡的每一单元的概率。
滑坡图和滑坡清单图之间的可视化比较，表明在绘制的滑坡区和滑坡概率大的地区之间具有良好的一致性。另外，滑坡概率大的地区位于已滑坡区以外。这样，罗吉斯回归能够划定可能发生滑坡灾害的区域。
将滑坡灾害图的概率范围采用指数比例尺分为4个区间，最大的区间包括滑坡发生概率大于10%的区域，这些区域包括沿密苏里河的陡岸和多山地区。最低的间隔指滑坡发生概率小于0.1%的地区，这些区域是沿山脊和谷底相对较小的区域。多山分区不是典型的完整研究区；低区间范围在研究区西部较大。数据检查表明所有的区间都包括滑坡，但是在最低的区间滑坡相对较少。可能是较低的区间包括小的易于发生滑坡的区域，例如接近沿溪流的切割岸的台地，30m的DEM对于确定如此小面积的陡坡太粗糙了。图中的4个区间不可能适用于所有目标，其他的使用者可根据特定用途，定义可供选择的合适滑坡灾害等级。

DEM的栅格间距和研究区的最大当地地形起伏控制可分辨的最大坡角。测量最大可分辨坡角可以通过求当地地形起伏除栅格间距的反正切获得。这样，最大可分辨斜坡是68º。采用ArcView GIS根据DEM检验斜坡程序表明最大滑坡角是66º。这样ArcView GIS斜坡程序确定最大值附近的坡角。真实的10mDEM可能具有最大可分辨坡角82º。然而，由于可以将山坡分为更小的单元，用ArcView GIS测量的最大坡角可能大于66º，但不会高于82º。
不言而喻，罗吉斯回归方法的应用假定滑坡发生概率坡度的增加而增加。这一假定对于小坡度和中等坡度是正确的。然而，当坡角接近垂直将有争议，因为观测的是陡坡上薄层的浅层滑坡和土壤，当坡度超过特定角度时滑坡发生的概率可能会开始降低。非常陡的斜坡上可能没什么土壤，只有风化基岩。随坡角接近于垂直，岩崩和倒塌的可能性增加。浅层滑坡区也可能是未被确定的深度滑坡。这些观测结果有力地支持了假定的方法。
表1中采用系数值概率等式进行的简单敏感性分析证实坡度是确定滑坡发生概率的最重要因子。如果地质状况是连续的，坡角在0～90º，对所有的地质单位，估计概率在0～100%。如果坡度是连续的，而地质状况不同，坡角是60º时估计概率的最大差别是49%。这样，如果坡角不同，概率范围是地质状况不同的两倍。
如表1所示基岩单位（页岩和石灰岩）最易发生滑坡，页岩具有最高的敏感性（最大负系数）。将负系数转化为图1中的曲线，不难看出左边曲线在同样坡角情况下滑坡发生的概率较高。这些结果支持了观测结果，如Lawrence岩层更易发生滑坡。似乎有些令人吃惊，石灰岩比非固结体对滑坡活动更敏感。然而，石灰岩存在于陡坡区，而且石灰岩岩层包括页岩夹层，该夹层组成了整个岩层厚度的50%。
有两种情况可能说明冰碛是对滑坡最不敏感的地质单元的事实。冰碛具有很大的可变物理属性，因为它是物质和粒度的不均匀混合物。此外，冰碛物出露于坡度相对较低的研究区西部。黄土比冰碛物对滑坡活动更敏感，但是黄土和冰碛物之间有95%的置信区间重叠，表明这些单元对于滑坡的敏感性在统计学上是相同的。
表1中冲积层系数不等同于0，说明冲积层单元具有平均的敏感性，与单独采用斜坡值区别不大。然而，冲积层系数的标准偏差最大，反映了与格栅斜坡数据相联系的事实，冲积层位于坡度最小的区域。

地下水和降水是两个附加因素，在本文的分析中没有进行考虑。有学者认为在亚粘土冲积物斜坡有4个因素控制迅速的水位波动：降雨强度、降雨持续时间、水力传导性和前期降水情况。不像在临界降雨量边坡失稳概率增加的泥石流，对于土崩，前期降水数量会影响孔隙压力。这样，同样的暴雨，当土壤潮湿时会造成滑坡，但当土壤干燥时不会引起滑坡。不幸地是，不可能将降雨强度、降雨持续时间和前期降水情况结合起来，而且可能会导致无穷大的多元罗吉斯回归模型。
总之多元罗吉斯回归是绘制滑坡灾害图的有效统计方法。在艾奇逊研究区，坡度和地质状况是评价滑坡发生可能性的最好预测变量。土壤类型与地质状况密切相关，在地质单元缺乏有效的分辨率的区域内可以由预测变量代替。尽管在其它地区坡向是导致滑坡发生的重要因素，但在本研究区不是主要预测因子。