地质灾害引发土地利用响应机制研究——以四川绵竹为例

李璇琼，梁延龙

(1.四川建筑职业技术学院，四川 德阳 618000；2.西南交通大学地球科学与环境工程学院，四川 成都 610031)

中国是一个地质灾害频发的国家，50%以上的国土面积受到地震、崩塌、滑坡、泥石流等地质灾害的影响，造成耕地损毁、房屋损毁和水土流失加剧。而土地利用类型、程度的变化又会间接引发地质灾害。因此，地质灾害与土地利用关系的研究具有重要意义(谭梦，2013；李磊，2019；彭义，2014；曹慧，2007；刘惠敏等，2001)。国内外学者多数是在某一时刻就某一具体地质灾害实例研究二者之间关系，而较少进行地质灾害分布规律与土地利用时空耦合关系的挖掘(Burby R.J,2000；Smith,S.M,2001；D.H.Bell,1985；Abuk,A,2001；Dai,F.C,2001；Michael,J.C,2004)。

笔者选取灾害频发区绵竹市在汶川地震期间及2009～2018年震后恢复有代表性的时段，从宏观和微观2个尺度进行地质灾害土地利用机制分析。以遥感影像和实测地灾点数据为基础，采用数字图像处理、监督分类等遥感技术，结合地理信息系统的区域统计、缓冲区分析、叠加分析、核密度分析、地理统计等方法探究地质灾害分布特征。通过多空间多时间尺度探索地质灾害与土地利用结构、幅度、速度、利用程度、地震烈度、土地敏感性及水土流失的耦合机制，总结出一套地质灾害对土地利用响应的研究方法，并对灾区土地利用合理规划及地质灾害的防治提供决策依据。

1 研究区概况

绵竹市位于四川盆地西北部，背倚龙门山脉，地处北纬31°09′～31°42′，东经103°54′～104°20′；面积为1 245.3 km2，自西北向东南伸展，东西宽约42 km，南北长约61 km；地势西北高(属龙门山地区)，多为山地，河流纵横，切割强烈，切割深度为500～1 000 m，河床狭窄，河谷陡峻，东南低，(属成都平原的一部分)，多为平原，由西北至东南逐渐倾斜，海拔高度为504～4 405 m，相对高差3 900余米，高差悬殊；山地、平原界限分明，地貌类型多样，有“六山一水三分田”的特点。

研究区地质灾害频发，尤其是2008年5月12日的汶川大地震引发了滑坡、崩塌、泥石流次生地质灾害。调查结果表明，震后新增灾害点142处，其中滑坡78处，崩塌129处，泥石流41处。造成房屋毁埋，耕地及植被遭到破坏，环境承载力明显下降。

2 研究区地质灾害分布特征

2.1 地震烈度分布

根据中国地震局公布的汶川地震烈度，绵竹烈度包含了Ⅶ度、Ⅷ度、Ⅸ度、Ⅹ度4个等级，总体较高。根据面积统计出Ⅹ级烈度区域占市总面积41.74%；Ⅸ级烈度区域占31.41%；Ⅷ级烈度区域占26.34%；Ⅶ级烈度区域占0.51%。研究区总体地震烈度由南至北逐渐升高(图1)。将地灾点叠加地震烈度分布(图1)，显示崩塌、滑坡、泥石流等地质灾害主要分布在地震烈度等级Ⅹ级和Ⅺ级的区域。

图1 绵竹地震烈度分布图Fig.1 Seismic intensity distribution of Mianzhu

2.2 地灾点分布

研究区地灾点分布见图2，显示地灾点主要分布于研究区北部和中部地震烈度较高的区域。

2.3 地质灾害分布与地形相关性

2.3.1 地质灾害点与地形关系

将地灾点叠加在高程分布图上(图2)，显示地灾点主要分布在高程600～3 000 m的中高地形区域。同样叠加坡度分布图(图3)，可知研究区地灾点主要分布在坡度20°～60°的中高坡度区域。

图2 地灾点与高程分布关系图Fig.2 Relationship between disaster and elevation

图3 地灾点与坡度分布关系图Fig.3 Relationship between disaster and slope

2.3.2 地震烈度与地形关系

研究区地震烈度随高程、坡度升高而增加。北部地形高，地震烈度大，南部地势低，地震烈度小。

2.4 地质灾害分布与沟壑密度相关性

沟壑密度的定义为流域中干支流总长度和流域面积之比，是地形发育阶段和地表抗蚀能力的重要特征值，是描述地形切割破碎程度的一个重要指标。沟壑密度越大，地面径流和冲蚀侵蚀越强烈，沟蚀切割发展越快。因此，沟壑密度对地质灾害的发育有较为重要的影响。

利用绵竹1∶5万DEM数据，通过水文分析工具进行河网和流域提取，共提取出小流域1 345个。基于ARCGIS空间分析平台统计每个小流域内的河网长度及流域面积，二者相除计算出沟壑密度。计算结果及分布情况见图4。其显示区内沟壑密度集中分布在0.60～2.96 km/km2区间内，整体较高，说明工作区沟壑密集，地形较为破碎。将沟壑密度重分类成5个等级，0～0.6 km/km2为低沟壑密度区域；0.6～1.19 km/km2为中低沟壑密度区域；1.19～1.91 km/km2为中沟壑密度区域；1.19～2.96 km/km2为中高沟壑密度区域；2.96～5.41 km/km2为高沟壑密度区域。

图4 小流域沟壑密度分布图Fig.4 Gully density distribution of watershed

叠加研究区地质灾害分布点可知，崩塌、滑坡、泥石流分布在沟壑密度较高的区域。统计各沟壑密度等级区域面积比例，得到饼状图(图5)。表明研究区中低和中沟壑密度区域范围最大，高沟壑密度区域范围最小。统计2007、2009、2018三期各沟壑密度等级区域未利用地面积比例，得到折线图(图6)。揭示各沟壑密度等级中的未利用地比例的都为2007～2009年增大，2009～2018年减小。其中，高沟壑密度等级和中高沟壑密度等级变化的幅度最大。说明地震期间中高和高沟壑密度区域土地破坏最为严重，震后得到了有利的恢复，这与地质灾害点分布在中高和高沟壑密度等级区域相对应。利用空间分析工具提取各泥石流、滑坡和崩塌点的沟壑密度值，生成盒须图(图7)。显示泥石流灾点分布的沟壑密度最高，崩塌和滑坡次之。表明研究区沟壑密度对地质灾害的影响程度为泥石流>崩塌>滑坡。

图5 各沟壑密度等级面积统计图Fig.5 The area statistics of each gully density range

图6 各沟壑密度等级未利用地变化趋势图Fig.6 Unused land variation trend of each gully density range

图7 地灾点沟壑密度分布图Fig.7 Gully density distribution of disaster spots

3 地质灾害土地利用宏观响应

土地利用宏观响应分为地震期及震后恢复2个阶段，研究整个绵竹市土地利用结构、幅度、速度、程度、地震烈度、土地敏感性、沟壑密度及水土流失响应(李璇琼,2010,2015；李鹏杰,2012；魏巍,2012；金雄兵,2003)。

3.1 研究区土地利用信息提取

笔者数据源采用绵竹2007年5月6日 、2009年6月12日LandSat5影像及2018年6月5日LandSat8共3期遥感影像，轨道编号为129-38。遥感图像的分类方法有很多，监督分类应用较为广泛。监督分类又称训练场地法，是以建立统计识别函数为理论基础，依据典型样本训练方法进行分类的技术。可通过反复检验训练样本，来提高分类精度。

本研究结合第三次全国土地调查土地分类系统和《土地利用现状分类(GB/T21010-2007)》，根据研究区具体情况，利用监督分类的最大似然算法将研究区土地利用类型分为6类。第一类为耕地，对应2个分类体系的一级地类耕地；第二类为林地，对应2个分类体系中一级地类林地；第三类为水域，对应2个分类体系中的一级地类水域及水利设施用地；第四类为建设用地，对应《土地利用现状分类(GB/T21010-2007)》中的一级地类建设用地，包含了第三次全国土地调查土地分类系统体系中商服用地、工矿仓储用地、住宅用地、公共管理与公共服务用地、特殊用地等一级地类；第五类为道路，对应2个分类体系中的交通运输用地；第六类为未利用地，对应2个分类体系中的一级地类其他土地。分类结果见图8～图10，图中显示研究区林地、耕地资源较为丰富，河流分布较为密集。

图8 2007.05.06土地利用分类图Fig.8 Land use classification on 2007.05.06

图9 2009.06.12土地利用分类图Fig.9 Land use classification on 2009.06.12

图10 2018.06.05土地利用分类图Fig.10 Land use classification on 2018.06.05

3.2 土地利用结构响应

利用GIS的区域统计功能，计算研究区不同时期各土地利用类型面积，分析面积变化，结果见表1。显示地震期间绵竹土地利用结构发生了较大的变化。林地、耕地、建设用地大幅度减少，未利用地增加，主要是由于地震导致土地损毁，转化为未利用地。震后随着绵竹城市灾后重建、城镇化进程及经济快速发展，土地利用结果也明显变化，耕地与未利用地减少，建设用地、道路增加明显。

表1 绵竹土地利用分类统计结果表(m2)Tab.1 Statistical classification results of TM images(m2)

3.3 土地利用幅度响应

通过某时间段一种土地利用类型变化量与起始时间面积比来研究土地利用幅度演化，体现不同类型的土地利用情况在数量上的波动。通过运算，具体变化幅度结果见表2。

表2显示2007～2009年建设用地、林地、耕地损毁分别为8.57%、6.99%和2.16%；未利用地大幅度增长了一倍。2009～2018年建设用地增加40.93%，道路增加14.94%。

表2 土地利用类型动态幅度变化(m2)Tab.2 Change of dynamic amplitude of land use type(m2)

3.4 土地利用速度响应

3.4.1 单一土地利用类型动态度

为了量化研究时段内土地利用类型的变化情况以及土地类型的动态度，引入单一土地利用动态度的概念，在数量上反映研究区土地利用类型面积的年变化。根据相关公式计算出绵竹市各时间段单一土地利用动态度数值(表3)。

表3揭示绵竹建设用地、耕地、林地在地震期间呈现减弱的趋势，地震的破坏导致了未利用地的极大增长，建设用地损毁速度最快。震后林地、建设用地、道路、水体都呈“涨势”，耕地和未利用地呈“落势”，道路增长速度最为迅速。总体时段道路、未利用地与建设用地增加，耕地、林地、水体减少，建设用地增加趋势最快，耕地减少趋势最快。

表3 单一土地利用动态度表(年)Tab.3 Single dynamic amplitude of land use(Year)

3.4.2 综合土地利用类型动态度

综合土地利用动态度的值反映该区土地利用综合年变化率。经计算得出地震期间、震后恢复、总体时段年综合土地利用动态度分别为1.19%、0.36%、0.11%。说明地震期间年土地利用变化速度最快，地震为主要影响因素。

3.5 土地利用程度响应

土地利用程度综合反映了土地自然属性与人文因素的综合效应。根据刘纪远等提出的土地利用程度的综合分析方法(Liu J-Y,2000)，计算出研究区2007年、2009年、2018年土地利用综合程度指数分别为249.64、243.06、251.49。地震期间耕地、建设用地、林地转化为未利用的灾毁土地，土地利用程度也就处于衰退的趋势。震后随着城市建设、交通用地增加，土地利用程度总体增强。

3.6 土地利用地震烈度响应

分析统计地震期间每个烈度区域土地利用转化面积及占该区域面积比例见图11。其显示地震烈度Ⅺ级区域耕地破坏比例为3.49%，烈度最大的X级区域破坏比例为6.23%，依次增加。建设用地和林地损毁的趋势都随着地震烈度增强而升高。地震烈度的土地利用响应程度依次为建设用地>林地>耕地。

图11 各地震烈度级别土地利用转化图Fig.11 Land use transformation in each earthquake intensity range

3.7 地质灾害引发的水土流失响应

笔者利用ARCGIS与MATLAB相结合，对水土流失相关的土地利用、植被盖度、坡度、沟壑密度及降雨量等5个因子编程，实现遗传神经网络(GA-BP)算法(李璇琼，2012)，得到最终的水土流失分结果图。算法流程见图12，结果见图13～图15。

图12 GA-BP模型建立流程图Fig.12 Flow chart of GA-BP model

图13 2007.05.06水土流失分级图Fig.13 Soil and water erosion on 2007.05.06

图14 2009.06.12水土流失分级图Fig.14 Soil and water erosion on 2009.06.12

图15 2018.06.05水土流失分级图Fig.15 Soil and water erosion on 2018.06.05

叠加地质灾害分布点，统计每个水土流失等级崩塌、滑坡、泥石流灾害点的频度，结果见表4。显示2007年震前，轻度和中度区域灾害点分布较多，2009年减少，2018年经过恢复呈“涨势”。2007年震前，剧烈、极强和强度范围地灾点频度较低，2009年震后增加，2018年经过恢复频度呈“落势”。

表4 各水土流失等级地灾点分布表(个)Tab.4 Distribution of disaster spots of each soil and water erosion range (Quantity)

4 地质灾害土地利用微观响应

微观响应利用崩塌、滑坡、泥石流点缓冲区为研究区，分析地质灾害土地利用结构、类型、幅度、速度、程度、沟壑密度及水土流失响应。

4.1 地灾点缓冲区建立

根据前人研究经验结果，结合研究区实际情况，通过实验最终选取1 500 m作为缓冲半径(杨德生，2011)。以采样点为中心，建立半径为1 500 m的缓冲区，再分析缓冲区内地质灾害土地利用响应(图16、图17、图18)。

图16 崩塌点缓冲区图Fig.16 Buffer of collapse points

图17 滑坡点缓冲区图Fig.17 Buffer of landslide points

4.2 土地利用结构响应

对崩塌、滑坡、泥石流地灾点缓冲区内的土地利用结构进行统计分析，结果见表5。说明崩塌、滑坡、泥石流等地质灾害对林地的影响最大，影响程度为林地>耕地>建设用地。

表5 地灾缓冲区土地利用类型结构表Tab.5 Land use type structure of disaster buffer range

4.3 土地利用幅度响应

对崩塌、滑坡、泥石流地灾点缓冲区内的土地利用结构进行统计分析，结果见表6。表中显示地震期间崩塌、滑坡、泥石流缓冲区未利用地分别增加了1 803.48%、2 315.46%和3 716.82%，高出整个绵竹未利用地比例(表2)的18～37倍。说明对地质灾害区域土地利用影响非常大。耕地、林地受滑坡影响最强烈，建设用地损毁比例在泥石流缓冲区最大。

表6 土地利用类型幅度变化表(年)Tab.6 Change of land use type amplitude(Year)

4.4 土地利用速度响应

土地利用速度响应主要通过年综合土地利用动态度分析。具体见公式1。

(1)

公式(1)中，LUi为监测时间段初始状态下第i类土地利用类型面积；△LUi-j为研究时段内第i类土地利用类型转为非i类土地利用类型面积和非第i类土地利用类型转为i类土地利用类型面积之和；T为研究时段长度。

利用公式计算出地灾点缓冲区地震期间、震后恢复及总体时段的年综合土地利用动态度，结果见表7。显示地灾点缓冲区动态度整个绵竹的25～27倍主要由于地震影响引发次生灾害。响应速度为泥石流>滑坡>崩塌。说明泥石流引发灾害区域土地利用类型速度变化最快。

表7 地灾点缓冲区综合土地利用动态度表(年)Tab.7 Dynamic attitude of comprehensive land use in disaster range(Year)

4.5 土地利用程度响应

土地利用程度综合指数模型是将土地利用的综合指标进行数学整合，形成一个连续分布的综合指数，其数值的大小综合反映了某一地区土地利用的程度。计算公式为:

(2)

公式2中，Li为某研究区域土地利用程度综合指数;Ai为研究区域第i级土地利用程度分级指数；Ci为研究区域内第i级土地利用程度分级面积百分比；n为土地利用程度分级数。据公式2计算出地质灾害点缓冲区的土地利用程度综合指数见表8。显示灾害点缓冲区土地利用程度均低于研究区综合水平。地震期间地灾点缓冲区土地利用程度减弱，灾害恢复阶段增强。滑坡区地震期间土地利用程度衰退幅度最大，泥石流区震后增加幅度最强。

表8 地灾点缓冲区土地利用综合指数表Tab.8 Land use composite index of disaster point buffer

4.6 地质灾害引发土地敏感性响应

为了从土地利用的角度探讨人类活动与地质灾害的相关性，引入了敏感系数的概念，以敏感性定量分析不同土地利用类型地质灾害的敏感性(何易平,2005;陈和平,2002;吴常润,2020;仲佳鑫,2018)。采用灾点法分析，具体见公式3。

SCi=In(DensClassi/DensMap)

(3)

公式3中,SCi为第I类土地利用类型的敏感性；DensClassi为某类灾害在类土地利用类型中分布的面积比率或出现的频率；DensMap为所有该类灾害总的面积比率或频率。对每一类缓冲区范围内的敏感性系数进行计算，结果见表9。

表9 各土地利用类型地质灾害敏感系数表Tab.9 Sensitivity coefficient of geological hazard of each land use type range

计算结果表明：地震期间耕地对滑坡敏感性增强；林地对崩塌、滑坡、泥石流敏感性减弱；建设用地对滑坡敏感性增强，对泥石流敏感性减弱。震后恢复阶段耕地、未利用地对各地质灾害的敏感性减弱；林地对各地质灾害敏感性增强；建设用地对滑坡敏感性减弱，对泥石流敏感性增加。耕地和林地对地质灾害敏感性最强。

4.7 地质灾害引发水土流失响应

以地灾点缓冲区为分析范围，生成核密度图(图19、图20、图21)。利用空间分析工具提取每个地灾点的核密度和水土流失值，进行回归分析，探索地灾点核密度和水土流失的相关性(图22、图23、图24)。揭示水土流失程度与地灾点核密度成线性正相关，地质灾害分布越密集，水土流失程度越严重。

图19 崩塌点核密度图Fig.19 Kernel density of Collapse points

图20 滑坡点核密度图Fig.20 Kernel density of landslide points

图21 泥石流点核密度图Fig.21 Kernel density of Debris flow points

图22 崩塌点核密度与水土流失相关性图(2007.05.06)Fig.22 Correlation between collapse kernel density and soil erosion (2007.05.06)

图23 泥石流点核密度与水土流失相关性Fig.23 Correlation between Debris flow kernel density and soil erosion(2009.06.12)

图24 滑坡点核密度与水土流失相关性(2018.06.05)Fig.24 Correlation between Debris flow kernel density and soil erosion(2018.06.05)

5 土地利用优化措施建议

通过以上计算、分析及研究，借鉴国内外先进的土地规划、利用理念，结合研究区实际，提出以下土地合理规划利用建议。

(1)因在25度坡以上的地区耕作，使得斜坡上部土层结构松散，土体粘聚力显著降低，降雨渗入土体饱和后容易形成水土流失，故应改为梯田方式进行耕作，并在耕地周围加强防护林建设。

(2)因研究区耕地损毁严重，震后城市化进程也占用一部分耕地。建议将研究区中部地势相对平缓的林地进行耕作活动，以增加耕地面积。

(3)研究区泥石流对建设用地影响较大，所以需加强排水设施建设。灾后恢复时可将靠近泥石流灾点的建筑搬迁到南部地势较为平坦的区域。在陡坡坡面设置护坡，既能涵养水源，也可减弱泥石流发生。

6 结论

土地是人类赖以生存的宝贵资源，研究地质灾害的土地利用响应机制对灾后土地合理利用具有重要意义。笔者以遥感和GIS手段相结合，对绵竹市2007～2018年间地质灾害引发的土地利用响应机制进行多尺度定量分析，得到以下结论。

(1)地质灾害分布。研究区崩塌、滑坡、泥石流等地质灾害主要分布在地震烈度等级高、高程为600～3 000 m的中高地形及坡度为20°～60°的中高坡度区域。地震期间中高和高沟壑密度区域土地破坏最为严重，震后得到了有利地恢复。沟壑密度对地质灾害的影响程度:泥石流>崩塌>滑坡。

(2)宏观响应分析。汶川地震导致林地、耕地、建设用地大幅度减少，未利用地增加。随着震后重建及城镇化进程，耕地与未利用地多转化为建设用地和道路。地震期间、震后恢复、总体时段年综合土地利用动态度分别为1.19%、0.36%、0.11%。3期土地利用综合程度指数分别为249.64、243.06、251.49。地震期间耕地、建设用地、林地转化为未利用的灾毁土地，土地利用程度也就处于衰退的趋势。震后随着城市建设、交通用地的增加，土地利用程度总体增强。地震烈度的土地利用响应程度依次为：建设用地>林地>耕地。崩塌、滑坡影响程度：林地>耕地>建设用地；泥石流影响程度：林地>耕地>建设用地>道路。地质灾害导致地震期间水土流失加剧，震后逐渐恢复。2007年震前轻度和中度流失区域灾害点分布频度高；2009年震后减少；2018年经过恢复呈“涨势”。震前剧烈、极强和强度水土流失地灾点频度较低，2009年震后增加；2018年经过恢复频度呈“落势”。

(3)微观响应分析。微观缓冲区崩塌、滑坡、泥石流等引发土地利用响应程度：林地>耕地>建设用地；年综合动态度：地震期间>震后恢复阶段>整体时段；响应速度：泥石流>滑坡>崩塌。耕地、林地受滑坡影响最强烈，建筑损毁受泥石流影响最大。地震期间滑坡地区土地利用响应速度最快，泥石流区响应程度最强。地质灾害点缓冲区范围的土地利用结构、幅度、速度、程度响应均高于整个绵竹市区。地震期间耕地对滑坡敏感性增强；林地对崩塌、滑坡、泥石流敏感性减弱；建设用地对滑坡敏感性增强，泥石流敏感性减弱。耕地和林地对地质灾害敏感性最强。地灾点核密度与水土流失呈线性正相关。

(4)对绵竹市灾后土地合理利用及地质灾害防治提出合理建议。例如,坡地改梯田、建立农田防护林、加强排水设施建设、损毁耕地异地置换等。