基于统计和聚类分析的黄土台塬塬边裂缝危险性分级

王 卓，许 强，魏 勇，李骅锦

(成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059)

黄土因其孔隙大、垂直节理发育、结构疏松、水敏性高，故易发生滑坡地质灾害[1-2]。据实地调查显示，仅陕西省黄土高原发生的地质灾害就约有1.5万起，其密度超过6起/km2，且其中85%为黄土滑坡[3]。黄土滑坡不仅频繁发生，造成的后果也极其严重，如：1971年5月5日，陕西省卧龙寺黄土滑坡共造成28人死亡，并毁坏了宝鸡峡引渭渠[4]；1983年3月7日，甘肃省东乡县发生著名的洒勒山高速黄土滑坡，仅30 s内体积约3.1×107m3的滑坡体急剧向南滑动1 km左右，瞬间淹没了1.3 km2范围内的农田，摧毁了3个村庄及一座小型水库，共造成237人死亡、22人受伤，成为20世纪80年代最具灾难性的滑坡之一[5]。近年来虽然国家已经对黄土地区地质灾害的治理给予了高度重视，但由于黄土特殊的结构性和水敏性，致使灾难性的滑坡事件仍在发生，如：2011年9月17日，西安灞桥区白鹿塬北坡发生黄土滑坡，滑坡体在滑动 320 m的过程中冲毁了大量厂房与宿舍，共造成32人死亡、5人受伤[6]；2015年8月中旬，陕西省山阳县黄土滑坡共造成64人被埋，经济财产损失无法估量[7]。由此可见，黄土滑坡已成为地质灾害研究的重中之重。

在黄土滑坡的演化过程中，塬边裂缝成为了坡体变形最直观的反映，不仅改变了地表水的入渗方式，更在一定意义上控制着黄土边坡的变形破坏[8-10]。赵宽耀等[11]、许强等[12]研究认为黄土中的节理裂隙为地表水的主要入渗通道，并通过试验验证了当甘肃黑方台地区地下水水位占黄土厚度比例为0.4时会引起坡体失稳；Zhang等[13]通过模型试验发现黄土塬边裂缝不但为入渗提供了优势通道，还会演化为每次新破坏的黄土滑坡后缘边界；伊龙等[14]结合模型试验和数值模拟总结出黄土塬边距离15 m之内不同深度的张拉裂缝对黄土边坡稳定性有着不同的影响；Xu等[15]也发现黄土塬边裂缝的深浅控制着黄土边坡失稳的模式。

随着对黄土滑坡塬边裂缝的深入研究，实现了基于塬边裂缝变形扩展对黄土滑坡的监测预警，近年来的几起成功预警案例极大地保障了人们的生命与财产安全[16-17]，同时也突显了研究黄土台塬裂缝的重要性。目前，对黄土滑坡塬边裂缝的塬边研究主要集中于塬边裂缝类型、成因机制等的探讨[18-20]，也有对塬边裂缝变形的力学分析[21-24]，但没有涉及塬边裂缝分级的研究。因此，为了研究塬边裂缝的发育阶段和分布与黄土滑坡的内在关系，本文以陕西省泾阳南源地区为研究区，在现场勘测获取的塬边裂缝资料的基础上，统计分析了2007—2016年以来该地区裂缝属性值的变化，并结合统计分类与K-means聚类两种方法对研究区塬边裂缝进行了危险性分级，划分出了几处危险性极高的塬边裂缝集中区，最后通过新滑坡事件验证了塬边裂缝危险性分级的可靠性。本研究对陕西省泾阳南塬地区防灾减灾工作具有一定的指导意义。

1 现场调查与数据收集

1.1 研究区黄土滑坡概况

本次研究区陕西省泾阳南塬地区以前为丝绸之路经济带西安自由贸易试验区以及西安国际港务区的核心地带，现今属于国家深入实施西部大开发战略所建设的西咸新区，自1976年大规模农业灌溉以来，区内共发生42处、60余起黄土边坡滑动。其中，1984年12月2日发生于泾阳南塬河滩村的蒋刘黄土滑坡代价最为惨痛，体积约为1.13×106m3的滑坡体摧毁了86间房屋，共造成20人死亡、20人重伤，致使河滩村行政级别取消[25]；而2003年7月发生的东风黄土滑坡和2004年3月发生的寨头黄土滑坡的方量大，分别约为1.6×106m3和1.5×106m3，造成近500亩土地被埋，经济损失严重，但所幸未造成人员伤亡[26]。

1.2 塬边裂缝现状

图1 泾阳南塬地区塬边裂缝调查区及局部示意图Fig.1 The crack investigation area of the South Jingyang Plateau area and partial sketch(a)塬边裂缝调查区；(b)67#裂缝;(c)4#裂缝;(d)4#裂缝的局部

研究区以泾河大桥为界，泾河南岸西至太平镇庙店村，东至高庄镇大堡子村[见图1(a)]。经2016年6月现场调查，黄土台塬塬边距10 m以内共见地表裂缝68条，其中太平段64条，高庄段4条；塬边裂缝总长为4 683.36 m，平均延展长度为68.87 m，最长可达577.9 m，平均张拉宽度为26.1 cm，平均错台高差为35.2 cm。其中，太平段的塬边裂缝均位于庙店、魏村黄土滑坡的后缘，多似67#裂缝[见图1(b)]，为滑坡牵引所致，分布密集且杂乱交错，错台高差大；高庄段的塬边裂缝位于蒋刘黄土滑坡后缘及塬边水渠旁，张拉宽度小、错台不明显；但蒋刘黄土滑坡后缘的4#裂缝[图1(c)、(d)]延伸约91 m，平均张拉宽度为52 cm，平均错台为130 cm，自身稳定性较差。

与2007年许领等[20]统计的研究区黄土台塬塬边裂缝发育情况相比(见表1)，裂缝平均每年新增4条，增长速率约为326 m/a，而裂缝平均宽、高却有所降低；但从2007年和2016年两期黄土台塬塬边裂缝对比柱状图(见图2)可以较直观地看出，塬边裂缝张拉宽度与错台高差的最大值仍在持续增长，塬边裂缝处于发育扩展阶段，活动强烈，一旦形成整体贯通，势必影响台塬黄土边坡的稳定性。

表1 2007年和2016年两期研究区黄土台塬塬边裂缝的统计对比

图2 研究区2007年和2016年两期黄土台塬塬边裂缝对比柱状图Fig.2 Comparison histogram of the two-phase cracks in the loess platform edge in 2007 and 2016

2 研究方法

以往研究裂缝时常以整条裂缝的平均值作为基础，而现实中野外裂缝发育不均匀，平均值并不能准确代表裂缝的实际特征。在裂缝数据处理时，本文分别选用了条、段(以RTK数据点为基础，有拐点处分为一段)、点3种控制模式对塬边裂缝的延展长度、张拉宽度和错台高差进行了统计分析，对比发现裂缝段的属性特征最为符合塬边裂缝的实际特征，故本文以裂缝段的形式对塬边裂缝进行描述，并结合统计分类和K-means聚类两种方法对研究区塬边裂缝危险性进行了分级。

2.1 统计分类

为了研究裂缝段属性值与频数分布的相关性，本文采用频数直方图对裂缝段数据进行了处理，但由于直方图会丢失部分潜在信息，且常用组距的频数直方图往往无规律呈现，故在本研究中希望通过不断调整直方图的组距，使其呈现出常用组距无法反映的潜在信息，具体的研究思路见图3。

图3 裂缝段数据的统计分类流程图Fig.3 Classification flow chart for crack segment data

以裂缝段张拉宽度的频数分布图为例，通过对组距的不断调整，发现无论组距如何变化，裂缝段张拉宽度在大于10～12 cm后频数开始大幅下降，而其在大于54～60 cm后频数趋于平缓，且频数分布在组距为6 cm时裂缝段张拉宽度阶梯性降低现象最为明显，见图4。根据裂缝段张拉宽度的阶梯性下降现象可知：裂缝段张拉宽度的频数分布不仅代表了裂缝的整体情况，更能反映某一条裂缝段在其张拉宽度扩大的过程中自身稳定性降低的过程，即裂缝段张拉宽度在0～12 cm时稳定性最高，其在中间部分时稳定性逐渐降低，而当裂缝段张拉宽度大于某一阈值后稳定性达到最低状态，随台塬局部失稳裂缝段条数减少。

图4 不同组距条件下裂缝段张拉宽度频数分布图Fig.4 Analysis of interval adjustment of the tension width of crack segments under different group spacing conditions

运用相同的方法可以分析裂缝段延展长度和错台高差的组距调整，最终得到了组距分别为1.2 m、9 cm时裂缝段属性频数阶梯性降低的直方图，见图5。

图5 裂缝段属性频数直方图Fig.5 Attribute frequency histogram of crack segments

由图5可见，在一定范围内，裂缝段属性频数较高且平稳，而随着塬边裂缝的发育和扩展，台塬出现局部崩塌或滑坡，使发育充分的裂缝条数减少，极大值段的裂缝段属性频数降低，这实际上就是塬边裂缝的发育及贯通导致台塬失稳的过程。

通过分析裂缝段属性频数阶梯状分布的现象，发现裂缝段属性频数在100以上及20以下时均较稳定，而在100至20间逐渐降低，由此将裂缝段属性频数为100和20的两个横坐标点定为裂缝段属性频数分布的特征值点，频数小于20的属性值对应为极大值段，频数在20～100间的属性值属于过渡段，频数大于100的属性值对应为极小值段(图5中裂缝段延展长度前部也应属于极小值段)，并根据裂缝段数据的实际意义将这三个值段分为了3个危险性等级，见表2。

由于塬边裂缝的长、宽、错台3种属性相互独立且对裂缝稳定性的影响均较明显，故认为此3种属性权重相等，分级时只考虑各值段的权重问题。本文运用MATLAB工具，对裂缝段3种属性的危险性按照表3的规则进行组合，将研究区68条裂缝的773个裂缝段划分为5个危险性等级，得到极高危险性至极低危险性裂缝段条数分别为6个、35个、109个、146个、477个。

表2 裂缝段属性各值段危险性等级

表3 裂缝段属性组合危险性分级

注：A表示裂缝段延展长度；B表示裂缝段拉张宽度；C表示裂缝段错台高差；脚标1、2、3分别表示极大值段、过渡段、极小值段。

2.2 聚类分析

聚类分析作为一种最典型的非监督式机器学习方法，考虑了样本数据的物理意义，并根据样本间的相似度以统计形式将数据进行自动聚类，该聚类算法定义清晰，可用于研究对象的定量评价[27]。聚类分析在滑坡危险性评估方面运用广泛，如毛伊敏等[28]采用UM-Chameleon聚类算法，以延安宝塔区为例，建立了适用于大区域的滑坡危险性评估模型；桂蕾等[29]运用K-means聚类法对巴东县新城区滑坡灾害进行了危险性分区；黄发明等[30]结合支持向量机与聚类分析，得出了准确的万州区滑坡易发性分区图。可见，聚类分析在滑坡地质灾害的研究中适用性较强。

聚类分析的本质是根据样本距离的远近将样本数据按照一定的方法分为若干个类别，使数据达到组内距离最小、组间距离最大。而聚类效果的好坏主要依赖于两个因素：距离度量方法(distance measurement)和聚类算法(algorithm)。

在本研究中，根据裂缝数据的变量类型，选择了适用于度量数值变量距离的euclidean距离法对样本数据的相似度进行计算，即：

(1)

式中:xi、zj分别为第i个样本数据和第j次聚类中心；n为样本指标属性个数，即裂缝属性个数；xai为第i个样本数据的第a种属性的值；zai为第a种属性的第j次聚类中心值。

由于样本数据量较大，故选择目的明确而又简便快捷的K-means聚类算法，其准则函数为误差平方和(sum of the squared errors，SSE)，其计算公式为

(2)

式中：Ci为第i个簇;x为Ci中的样本点;zi为Ci的质心(聚类中心)通过不断迭代趋于稳定;SSE为所有样本的聚类误差，代表了聚类效果的好坏。

随着聚类数k的增大，样本划分更加精细，每个簇的聚合程度逐渐提高，而误差平方和SSE则逐渐变小。一般来说，存在一个最优k值使得SSE的下降曲线由骤减趋于平缓，而本文中由于k值较明确，故省略探讨，直接取k=5。K-means聚类算法的具体步骤如下：

(1) 任意选择5个对象作为初始聚类中心。

(2) 计算其他样本数据到每个初始聚类中心的距离，并按照最小距离原则按下式将数据分配至最临近的初始聚类中心：

dij=min(‖xi-zj‖) (xi∈X,zj∈Z)

(3)

式中：X、Z分别为xi、zj的集合。

(3) 重新计算每个簇里所有数据的平均值作为新的聚类中心。

(4) 计算SSE，并判断SSE是否收敛，若未收敛则重复步骤(2)～(4)，直至收敛则结束聚类。

3 研究区裂缝段危险性分级对比与验证

3.1 裂缝段危险性两种分级结果的对比

本文以SPSS统计分析软件为操作平台，迭代17次后实现了对裂缝段的K-means聚类。据K-means聚类分析结果显示：危险性极高的2条裂缝段分别为LFD57和LFD177，其中与人为统计分类结果一致的为LFD57裂缝段。分析两种分级结果中极高危险性裂缝段条数不同的原因发现：LFD57和LFD177裂缝段的延展长度分别为53.26 m和38.78 m，是裂缝段延展长度最大的两个值，在聚类时尽管进行了标准化处理，但与裂缝段张拉宽度和错台高差的数值相比仍然较大，对聚类中心的选择有一定影响；而按照人为统计分类规则对裂缝段进行危险性分级，能较全面、合理地考虑裂缝段各属性对危险性分级的影响。

将K-means聚类结果与表3人为统计分类结果进行了对比(见表4)，发现对同一条裂缝段的危险性分级而言，两种裂缝段危险性分级结果的相似度达到0.804 7，满足数学意义上的显著相关(大于0.8)，说明根据K-means聚类算法的校正，人为统计分类结果比较合理。

表4 裂缝段危险性两种分级结果的对比

3.2 裂缝段危险性分级的验证

根据上述裂缝段的统计分类结果，运用ArcMAP对裂缝段危险性5个等级赋予不同色标，得到了泾阳南塬地区裂缝段危险性分级图，见图6。其中，危险性极高的6条裂缝段编号自东向西依次为LFD57、LFD75、LFD77、LFD186、LFD144和LFD285。

图6 泾阳南塬地区裂缝段危险性分级图Fig.6 Crack risk classification map of the South Jingyang Plateau area

由图6可见，高庄段调查区仅含一条极高危险性裂缝段，编号为LFD57，与LFD56和LFD58裂缝段共同组成了如图1(c)、(d)所示的4#裂缝；其余极高危险性裂缝段均分布于太平段调查区。

2018年10月对研究区6条极高危险性裂缝段进行现场调查发现，位于高庄段蒋刘黄土滑坡后缘危险性极高的裂缝段LFD57和危险性中等的裂缝段LFD56、LFD58(即4#裂缝)已随新滑坡事件滑落，现存台塬与4#裂缝的延展方向近似平行[见图7(a)]，推测4#裂缝在降雨和灌溉条件下为台塬塬边斜坡提供了渗流的优势通道，水体从裂缝流入后引起局部孔隙水压力的迅速上升，降低了边坡整体的稳定性，最终引发滑坡[31]。此外，在新滑坡两侧发现了两条塬边新生裂缝[图7(b)、7(c)]，其中如图6(b)所示新裂缝的延展长度达到了80 m，平均张拉宽度为40 cm，平均错台高差为140 cm，深不见底。可见，泾阳南塬地区塬边裂缝仍在不断发育扩展，其与黄土滑坡耦合作用的机理亟待研究。

蒋刘村塬边的新滑坡事件中危险性极高的裂缝段LFD57滑落，较好地验证了本文裂缝段危险性分级结果的可靠性。因此，对于新生裂缝[图7(b)、图7(c)]和位于太平段调查区的5个极高危险性裂缝段集中区(见图8，其中LFD453～LFD459裂缝段均为高危险性裂缝段)，今后需要对其引起足够重视并进行重点监测。

图7 研究区LFD56、LFD57、LFD58号裂缝段和塬边 新生裂缝位置图Fig.7 Location of crack segments LFD56,LFD57, LFD58 and new cracks in the study area(a)已滑落的4#裂缝;(b)新滑坡西侧塬边新生裂缝; (c)新滑坡东侧塬边新生裂缝

图8 泾阳南塬地区太平段调查区极高危险性裂缝段集中区Fig.8 Concentration area of extremely high-risk crack segments in Taiping section investigation area in the south Jingyang Plateau area

4 结 论

本文通过对陕西省泾阳南塬塬边裂缝的调查和统计，分析了黄土滑坡后缘裂缝的危险性等级及发育分布特点，得到的结论如下：

(1) 泾阳南塬地区共发育有塬边裂缝68条，总延展长度为4 683.36 m，主要集中在0～90 m之间，平均延展长度为68.87 m；裂缝张拉宽度和错台高差均集中于0～50 cm，平均张拉宽度为26.1 cm，平均错台高差为35.2 cm。

(2) 采用统计分类和K-means聚类两种方法对研究区塬边裂缝危险性分级结果进行了对比分析，两种分级结果的相似度可达到显著相关(大于0.8)，认为统计分级结果真实、可靠，对黄土滑坡灾害及塬边裂缝的研究具有一定的参考价值。

(3) 首次以裂缝段形式对塬边裂缝进行研究，基于裂缝段属性特征的危险性分级并通过新滑坡事件的验证表明：泾阳南塬地区塬边裂缝的危险性等级与黄土边坡的失稳有一定的内在关系，并且从裂缝段危险性分级结果可知，太平段调查区的5个极高危险性裂缝段集中区以及塬边新生裂缝位置有再次发生滑坡的可能性，需要对其引起高度重视。

(4) 本研究中未考虑裂缝的深度属性值和塬边距等对裂缝危险性分级的影响，分级结果具有一定的片面性。