做了这么多年GIS项目地质灾害风险区划图这个需求我接过不少但每次看到有人拿着DEM和一堆灾害点数据直接就开始出图我还是会心里一紧。这东西看着简单实际上牵扯到坐标系、栅格运算的精度、重分类的边界划定、密度分析的参数选择一个环节不谨慎出来的图就是“看起来很专业实际上根本不能用”。这篇文章我就把手头的完整流程梳理一遍。我会用ArcGIS把整条链路走通从DEM和灾害点数据的预处理开始到地形因子提取、灾害点密度分析、加权叠加最后到风险区划图的制图输出一共5个核心步骤。不管你是做自然资源相关项目的从业者还是高校里正在写相关论文的研究生只要跟着这套流程走至少能保证你的成果在技术逻辑上是站得住脚的。1. 整体设计思路风险区划图到底在“划”什么1.1 核心逻辑历史灾害密度与地形因子的叠加在做任何操作之前得先把原理捋清楚。地质灾害风险区划图本质上回答的是两个问题哪里已经发生过灾害以及哪里具备发生灾害的条件。“哪里发生过灾害”靠的是灾害点数据——通过核密度分析把离散的点状信息转化为连续的面状分布表达历史灾害的空间集聚程度。“哪里具备条件”靠的是DEM——通过提取坡度、坡向、高程等地形因子判断哪些区域处于容易发生崩滑流的地形环境中。把这两类信息按一定权重叠加起来得到的就是一张反映风险程度高低分区的图件。这里有一条重要的经验要提前说市面上很多教程会加入降雨、地层岩性、断裂带距离等更多因子但对于一个以DEM和灾害点为核心的快速评估场景地形因子加上灾害密度已经足够支撑初期决策。因子不是越多越好因子之间的相关性和数据的精度才是决定结果可信度的关键。我见过不少项目把七八个因子硬塞进模型最后主控因子被稀释结果反而不如三四个核心因子来得清晰。1.2 技术流派对比信息量模型、层次分析法与双指标叠加在动手之前还需要了解当前主流的技术方案这会影响你在ArcGIS里的具体操作路径。信息量模型基于贝叶斯原理统计每个因子分级中灾害点的占比与区域面积占比的比值然后取对数作为信息量值。这个方法对数据要求相对宽松也是国内地质灾害调查中最常用的方法之一因为它直接用灾害点去“学习”因子的贡献度。层次分析法通过专家打分构造判断矩阵计算每个因子的权重。适合数据不足但专家经验丰富的场景但主观性比较强评审时容易被质疑。双指标叠加先分别计算“易发程度”和“灾害密度”再做二维交叉分类。这个方法可视化效果好但分区边界往往比较琐碎。我这套5步流程结合了信息量模型的核心思想和实操便利性地形因子用重分类后的等级分值来表征灾害点数据用核密度来表征然后通过加权叠加完成区划。纯信息量模型需要大量统计计算而ArcGIS里的Reclassify和Weighted Overlay工具在表达“条件组合”时已经足够灵活。如果最后你发现结果分区与已知灾害点分布明显不符再考虑引入信息量模型修正权重也不迟。2. 数据准备与坐标系统一2.1 DEM数据获取与预处理要点DEM是整个流程的地形基础。常见的免费数据源包括ASTER GDEM30米分辨率、SRTM30米或90米、ALOS AW3D3030米以及部分省份提供的5米或10米高精度DEM。选择依据只有一个比例尺需求。做县级尺度的风险区划30米分辨率完全够用做乡镇级或者具体隐患点的精细评估就需要想办法拿到更高分辨率的DEM。拿到DEM之后第一步要做的不是急着提取坡度而是检查数据范围。很多从网上下载的DEM是分幅存储的覆盖范围远大于研究区。这时候需要用到【Spatial Analyst Tools - Extraction - Extract by Mask】用研究区边界矢量数据去裁剪DEM。注意掩膜数据本身必须和DEM坐标系一致否则会得到一张空图或者位置偏移的图。这是我认为整条流程中最容易出错但最容易被忽视的细节。2.2 坐标系检查定义投影和投影变换的区别凡是搞GIS的人几乎都在这上面翻过车。ArcGIS里有两个看似相近的操作Define Projection和Project Raster/Project。Define Projection是“声明”数据当前的坐标系不改变任何坐标数值。比如你下载的DEM没有投影信息ArcGIS不知道它是什么坐标系导致它在地图里位置不对这时候用Define Projection告诉系统“这是WGS84经纬度”或者“这是CGCS2000 3度带”。Project才是真正的坐标变换。比如DEM是WGS84经纬度单位度但你的灾害点是CGCS2000高斯投影单位米两者叠加时单位不同、坐标系不同就必须用Project Raster把DEM转换到CGCS2000投影坐标系下。一个实操建议在做栅格分析之前先把所有矢量数据灾害点、研究区边界和DEM都统一到同一个投影坐标系下。我的习惯是统一到CGCS2000 / 3-degree Gauss-Kruger CM 120E或者根据研究区经度选择对应的中央经线带。原因是后续计算坡度时的单位必须是米如果保持经纬度ArcGIS会用一个经验公式来估算平面距离在小范围内影响不算太大但凡是涉及面积、长度统计的地方误差就会被放大。提示统一坐标系之后务必用【Identify】工具在DEM和矢量图层的同一位置点一下对比坐标值是否一致。这一步一分钟都花不到却能避免后面所有分析结果全部偏移的灾难。2.3 灾害点数据清洗与属性检查灾害点数据一般来自历史地质灾害调查数据库常见格式是shapefile点文件属性字段包括灾害类型滑坡、崩塌、泥石流、发生时间、规模等级等。拿到数据后要做三件事坐标检查用ArcGIS显示点的空间位置和已有影像或地形做比对排除明显落在河流中心、山体外部的异常点。字段清洗删除无效值比如经纬度为0的记录保留灾害类型字段用于后续的权重调整或分类统计分析。同一位置去重有些灾害点会在不同年份被重复记录需要在属性表中做删除重复项处理。如果不处理核密度分析时这些重复点会被重复计数导致局部密度虚高最终影响风险分区的可信度。3. 地形因子提取与重分类3.1 坡度、坡向、高程三类因子的提取方法地形因子在ArcGIS里的提取非常成熟核心工具都在【Spatial Analyst Tools - Surface】下。坡度Slope使用【Slope】工具输入DEM输出栅格默认单位为度。坡度提取时需要注意Z因子Z factor设置。如果DEM是经纬度坐标系Z factor必须设置为111320即1度约等于111320米如果不设置算出来的坡度会整个偏小。这个坑我在早年间踩过当时算出来的坡度普遍差了一半还多。坡向Aspect使用【Aspect】工具输出以正北为0度、顺时针递增的坡向值。坡向本身是环形数据直接用原始值不适合参与权重叠加所以通常需要将坡向重分类为平面、北、东北、东、东南、南、西南、西、西北等类别。考虑到地质灾害尤其是滑坡多发生在阴坡或半阴坡水分条件更好、岩体风化更强烈坡向分级时阴坡的等级分往往更高。高程Elevation直接使用DEM原始值通过重分类划分为低、中低、中高、高等多个等级。高程分级区间需要结合研究区的实际海拔范围来确定常见做法是使用自然间断点法Jenks或等间距法。3.2 重分类的等级赋值原则重分类是整个流程中主观性最强、也最需要结合经验的一步。思路是把每一种因子按照“不利于稳定”的程度划分为4级或5级并赋予统一量纲的等级值比如1、2、3、4、55表示最易发。以坡度为例子0-15度为低风险区赋115-25度为较低风险赋225-35度为中等风险赋335-45度为较高风险赋4超过45度为高风险赋5。这个划分不是绝对的如果研究区整体坡度偏大比如整个区域都以30度以上的陡坡为主这时候就要把分级区间整体上调否则结果会集中在少数几个等级上失去区分度。重分类工具路径是【Spatial Analyst Tools - Reclass - Reclassify】。重分类完成后建议把每种因子的栅格逐一用【Extract Multi Values to Points】提取到灾害点上做一个简单的频率统计看看灾害点是不是真的集中分布在高等级区域。如果统计结果灾害点反而集中于低等级区域说明你的分级逻辑有问题需要回到上一步调整区间。4. 灾害点核密度分析与参数选择4.1 核密度分析的原理简述核密度分析Kernel Density可以理解为对每个待计算的栅格单元考察其邻域范围内的点状要素距离越近权重越大最终形成一个连续的光滑密度面。对于地质灾害而言一个地区灾害点越密集说明该地区的灾害活动性越强未来再度发生灾害的概率也相对较高。在ArcGIS中的操作路径是【Spatial Analyst Tools - Density - Kernel Density】输入灾害点要素选择密度字段。这里有一个细节如果灾害点属性中有规模等级字段可以把Population字段设置为规模等级值这样大型灾害点的贡献会被放大结果更符合实际风险含义。如果没有合适的字段直接使用None即可。4.2 搜索半径的选择默认值不一定合理核密度分析里最关键的参数是搜索半径Search Radius也就是统计邻域的范围。ArcGIS会默认计算一个适合当前数据范围的值但是这个默认值是基于点分布的平均中心距离推导出来的未必符合地质灾害的物理含义。我常用的做法是根据研究区面积和灾害点数量反推。假设研究区面积为S平方公里灾害点数量为N先计算平均每个灾害点对应的面积S/N然后反算搜索半径搜索半径 (S / (N × π)) 的平方根 × 一个经验调整系数系数通常取1.5到2之间。为什么要放大因为地质灾害不是点对点的独立事件而是受区域地质条件控制的群体现象一个灾害点的存在往往意味着周边地带也具有相似的条件所以搜索半径适当放大后密度面会更平滑也更符合地质规律。如果实测结果密度图过于破碎很多孤立的红色小斑块就调大搜索半径如果密度图模糊成一整片区分不出高低密度区就调小搜索半径。这个参数没有一个“放之四海而皆准”的值需要根据结果反复试。4.3 输出像元大小与面积归一化核密度分析还有一个容易被忽略的参数输出像元大小Output Cell Size。它的默认值通常取等于或略大于输入栅格的分辨率。如果DEM是30米建议把核密度输出像元也设置为30米这样后续叠加时各栅格分辨率一致不需要重采样也不会引入额外的插值误差。核密度分析完成后得到的是一个密度值栅格值的量纲是“每单位面积的点的个数”。因为这个值与点的绝对数量有关不同研究区之间缺乏可比性。如果只是做单期区划图直接用原始密度值做重分类并无不可如果打算做多期对比或者多个区域横向比较就建议做一个面积归一化密度值乘以研究区面积再除以点总数把结果映射到0-1区间。5. 风险等级加权叠加与制图5.1 权重赋予与栅格计算器公式到了这一步我们已经有了4类基础栅格坡度重分类栅格、坡向重分类栅格、高程重分类栅格、灾害点核密度重分类栅格。现在的问题是如何把它们合并成一张风险区划图。加权叠加推荐使用【Spatial Analyst Tools - Overlay - Weighted Overlay】或者直接用【Raster Calculator】写公式。Weighted Overlay的好处是工具内部会统一处理权重和等级表的映射关系不容易出错。Raster Calculator更灵活但需要手动保证各栅格量纲一致。这里给一组常用的经验权重可根据研究区实际调整坡度0.35灾害点核密度0.35高程0.15坡向0.15为什么坡度和灾害密度占大头坡度是地形因子中最直接控制斜坡稳定性的因素而灾害点密度是历史经验的基础两者合计占到70%的权重。高程虽然和灾害有相关性但在大多数区域它的作用是间接的高程变化往往伴随着坡度、植被、降水的联动变化。坡向的影响在湿润地区更明显但在干旱地区权重可以适当降低。如果使用栅格计算器公式可以写为风险指数 坡度重分类 * 0.35 核密度重分类 * 0.35 高程重分类 * 0.15 坡向重分类 * 0.15前提是各重分类栅格的等级值范围一致比如都是1到5。如果用的是Number类型这个表达式会直接输出一个浮点栅格。5.2 输出栅格的NoData处理加权叠加经常会遇到一个问题某个因子栅格在研究区边缘存在一块NoData区域导致计算结果在这些区域也为NoData。这是制图输出前必须处理掉的。处理方法有三种在【Extract by Mask】阶段确保所有因子栅格都使用同一个研究区边界裁剪边界外整体设为NoData边界内不应该有NoData。对于DEM本身存在的空洞比如部分卫星数据的云遮挡、水体区域没有高程值使用【Spatial Analyst Tools - Hydrology - Fill】工具或者【Con】条件语句做插值填补。如果NoData区域在边界内部零星分布可以用【Spatial Analyst Tools - Neighborhood - Focal Statistics】计算一个3×3邻域的均值然后用CASE语句赋值到NoData像元上。实际操作中我见过最多的情况是从不同来源下载的DEM和矢量边界坐标系底子就不一致导致掩膜裁剪后出现了一条宽窄不一的NoData“白边”。这本质上还是没把坐标系统一工作做扎实。建议回到第2章重新检查。5.3 风险等级划分与专题图制图得到风险指数连续栅格后下一步是把连续值划分为低风险、中风险、高风险、极高风险等几个等级。划分方法推荐使用【Spatial Analyst Tools - Reclass - Reclassify】里的自然间断点法Jenks。自然间断点法的优势是让每个等级内部的差异尽量小、等级之间的差异尽量大划分结果更符合数据本身的分布结构。如果没有特殊要求通常划分为4级低风险区中风险区高风险区极高风险区等级个数取决于出图比例尺和用途。如果是县级挂图4级足够如果是乡镇级排查可以考虑扩展到5级增加一个“较低风险区”。制图阶段有几个点值得单独说色彩方案风险等级建议使用渐变序列从绿色渐变为黄色再到红色这是国内地质灾害风险图的约定俗成配色不要在颜色上搞创新评审时的第一印象很重要。图例图例中必须注明数据来源、分辨率、坐标系、投影方式、制图时间等信息。灾害点叠加把原始灾害点按灾害类型用不同符号叠加上去可以直观验证分区结果与历史灾害分布是否吻合。如果大量灾害点落在低风险区说明权重配置或分级逻辑需要调整。出图分辨率如果用于印刷建议输出分辨率不低于300dpi如果用于网页发布栅格化输出时注意选择合适的像元大小避免出现过于明显的锯齿。6. 常见问题与排查技巧实录6.1 常见问题速查表问题现象可能原因排查方法坡度计算结果整体偏小DEM是经纬度坐标系未设置Z factor检查DEM坐标系设置Z factor为111320DEM和灾害点叠加后位置偏移坐标系不一致在ArcGIS中查看两个图层属性统一投影坐标系核密度图破碎成很多孤立斑块搜索半径太小调大搜索半径或用公式计算推荐值加权叠加结果边缘出现NoData各因子栅格范围不一致或DEM存在空洞统一掩膜范围填补DEM空洞灾害点全部落在低风险区重分类等级赋值逻辑有问题调整分级区间重新分配等级值出图后红色区域边界过于细碎输出栅格像元太小或分类数太多平滑处理栅格或减少分类等级重分类时等级值范围不同不同因子栅格的字段取值范围不一致在Reclassify时统一设定输出等级值范围6.2 排查思路与经验补充在实战中最浪费时间的问题往往不是技术难点而是数据源本身的底子问题。我有一次做某地级市的风险区划图DEM是网上下载的拼接数据下载时没仔细看元数据结果其中两块DEM的坐标系分别是WGS84和CGCS2000拼接后算出来的坡向在高程落差明显的区域出现了一条明显的“缝合线”。排查了将近半天才发现问题不在参数设置而在源头数据。从那之后我养成了几个习惯分享给大家动手前先建一个“数据体检”清单坐标系、分辨率、范围、NoData情况、字段完整性每一项都检查一遍全部通过再进入正式分析流程。每一步中间结果都要存成独立文件不要在一个流程里反复用临时图层ArcGIS的临时数据有时候不会立即释放容易导致后续步骤读到旧数据。最终成果出图前做一次“抽屉验证”把灾害点按时间分为两批比如前80%和后20%用前一批的历史灾害点做密度分析生成风险区划图再统计后一批灾害点落在高风险区的比例。如果比例偏低说明模型泛化能力不足需要回到权重配置调整。6.3 权重调整的定量分析方法如果需要对权重做定量调整而不是拍脑袋决定推荐用【GeoDetector】地理探测器或者简单的Spearman相关性检验来辅助。在ArcGIS中用【Extract Multi Values to Points】提取各因子值到灾害点再提取同等数量的非灾害点研究区内随机生成两套样本放在一起在统计软件里做相关性分析。相关性高说明该因子对灾害分布的贡献大权重应当提高。这个方法虽然会增加不少工作量但在写论文或者做正式的风险调查项目时它能极大提升成果的科学性。哪怕评审专家提出“为什么坡度占0.35”你也可以依据相关性分析结果给出一个定量的回答而不是笼统地说“经验值”。如果你的时间允许我非常建议在这个环节多投入一点精力。7. 最后分享一点个人体会做了这么多年的风险区划图我最大的感受是ArcGIS的操作确实是人人都能学会的但它不是一个“按钮就能出结果”的软件。每一个工具的参数背后都对应着对地质规律的理解和对数据质量的把控。在这套5步流程中最关键也最考验功力的是重分类的等级划分和权重的配置。同样是坡度因子在黄土高原和在高山峡谷地区分级的区间可能完全不同同样是核密度搜索半径在一个灾害点几百个的密集区和另一个灾害点只有几十个的稀疏区参数的取值逻辑也不一样。千万不要拿着别人的参数直接套用GIS不是复印机你研究区的数据特征才是唯一的决策依据。我的另一个心得是不管技术多熟练都要把过程记录下来。ArcGIS的Model Builder可以把整套流程模型化下次换一个研究区只需要替换输入数据就能一键跑完。即使你刚开始接触这套工具也建议在每一步操作后都截图记录参数等到出问题回溯时这些记录会帮你节省大量的排查时间。希望这篇文章能帮你少走一些弯路。等你用完这套流程做出你的第一张地质灾害风险区划图再回头来看这些参数和选择你会有比我写这篇文章时更深一层的体会。