基于Living Atlas数据为木里山体滑坡敏感性建模

滑坡是一种地质过程，山体滑坡时重力会导致岩石、土壤、碎屑或其混合物沿下坡移动。山体滑坡会严重破坏建筑物、道路及其他基础设施，并破坏野生野生动物栖息地。一场野火之后，由于植被的突然损失，山体滑坡的风险通常更大。
2020年3月28日19时30分，四川省凉山州木里县乔瓦镇锄头湾村与项脚蒙古族乡项脚村交界处发生森林火灾。

2020年7月1日四川省木里县项脚乡因连日持续强降雨引发大股泥石流，造成项脚乡阿牛窝子组境内山体滑坡、道路中断、农作物大面积受灾，村民房屋、集中安置点被淹没。

本文探索根据ArcGIS Living Atlas of the World 的地形、影像和降雨栅格图层生成滑坡敏感性表面，并根据此表面找到火灾后山体滑坡风险最高的区域。

查找数据

许多因素可导致山体滑坡风险增加，包括土壤组成、降雨、植被、坡度和坡向。本文将重点关注三个因素：植被密度、地形陡度和降雨量。
ArcGIS Living Atlas of the World是不断发展的精选和权威地理信息的集合。该集合中可以找到我们需要的影像数据和地形等数据。
由于笔者并未找到合适的降雨量数据，选择通过NOAA下载该区域部分测站的水文数据，通过差值形成降雨量数据。数据下载方法在ArcGIS pro通过聚合点创建时空立方体一文中有详细说明。

1.通过美国NOAA网站下载中国测站的水文数据，点击https://gis.ncdc.noaa.gov/maps/，选择日测数据

2. 定位到四川，框选木里县一定范围内的监测站，Add to Cart

3. 选择csv格式，数据日期本次选择了6-8月的降雨量数据，Continue，后续选择数据属性带有Precipitation降水量数据。

稍后会在邮箱中收到邮件，数据结果如下图（需要注意将DATE转换为yyyy-mm-dd格式）

降水量数据差值

1. 首先通过添加XY点数据，将降雨量数据添加到地图中
   
   

2. 在工具中搜索反距离权重法（IDW），将点插值成栅格表面
   

3. Z值字段选择降雨量字段PRCP，点击运行。
   
   需要注意NOAA的PRCP单位为英寸，需要乘25.4转为毫米
   

4. 得到插值结果
   
   拉伸后，降雨较多的区域以白色显示
   

搜索 Living Atlas of the World

获取Landsat

在ArcGIS pro中依次点击视图，目录窗格，门户。

选择Living Atlas，搜索landsat，选择Multispectral Landsat，并单击确定

在Mutispectral Landsat图层属性中，新建定义查询，选择 AcquisitionDate，选择灾后某一日的数据。（具体影像覆盖信息在USGS中查看会比较快，这里选择5月6日-8日间的数据）

处理模板可以根据需要选择，为了不影响分析，此次将处理模板更改为无，因此图像中的所有光谱波段均可用。

以下是几种常用的处理模板对比

Bands shortwave IR-2, shortwave IR-1, red (7, 6, 4)

bands red, green, blue (4, 3, 2)

Bands near-IR, red, green (5, 4, 3)

*Bands shortwave IR-1, near-IR, blue (6, 5, 2) *

区划过火范围

在工程数据库中新建要素类，命名为“木里328火灾”，以Landsat为底图，勾绘出过火范围（本次未精细区划）。

获取Terrain

在目录窗格中，单击门户选项卡，然后单击 Living Atlas 选项卡。

在搜索栏中，输入 terrain，然后按 Enter 键。
右键单击名为 Terrain 的项目并选择添加至当前地图。

Terrain 影像图层以默认灰度的可视化效果显示在地图上。较高的高程以浅灰色显示。

处理数据

现在已经收集了建模所需的三个栅格图层，其中landsat和terrain来自在线的ArcGIS Living Atlas of the World，因此未将任何数据下载到本地计算机。

将影像隔离到一个感兴趣的场景

Landsat 图层包含许多覆盖不同年份的世界场景，先前已经通过定义查询将图层更新为仅显示2020年5月6日至2020年5月8日之间捕获的场景，这是在木里森林火灾被完全遏制后的1个月后。你可能仅看到地图上显示的一个场景，但是那几天世界上的其他地点也捕获了数百个其他场景。接下来，需要使用空间查询来进一步过滤图层，使其仅显示覆盖木里森林火灾区域的一个场景。

• 单击功能区上的地图选项卡。在选择组中，单击按位置选择。

• 在地理处理窗格中，对于的输入要素，选择 Multispectral Landsat。

• 对于选择要素，选择之前区划的过火范围**“木里328火灾“”**
  

• 单击确定。所选场景的数量显示在地图视图的右下方。
  

• 在内容窗格中，右键 Multispectral Landsat，指向选择，然后单击锁定至选择内容。
  

• 此操作将确保仅显示覆盖木里森林火灾区域的场景。

通过将定义查询与空间查询相结合，Multispectral Landsat 影像图层已被过滤到一个满足日期和位置标准的场景。

更改 Terrain 处理模板和符号系统

与 Multispectral Landsat 图层相似，发布的 Terrain 图层也具有不同的处理模板。其中一些模板（例如 Slope_Degrees）可以将数据转换为可用于分析的图层；其他模板（例如 Slope_Degrees_Map）是处理后数据的可视表示，并且仅用于制图。

• 在内容窗格中，关闭 Multispectral Landsat 图层并打开 Terrain 图层。
• 选择 Terrain 图层后，单击功能区上的数据选项卡。
• 在处理组中，单击处理模板，然后选择 Slope_Degrees。
  
  Terrain 图层在地图上的外观随即发生变化。在内容窗格中，图例也会更新以显示 0 到 90 之间的值范围。较亮的区域具有更陡的坡度，最高可达 90 度。最暗的区域是平坦的，坡度为 0 度。
• 右键单击 Terrain 图层，然后选择符号系统。
• 在符号系统窗格中，单击配色方案菜单并选中显示名称。
• 向下滚动并选择坡度配色方案
  

现在，已经以绿色-黄色-红色配色方案描绘了坡度。下一步将统计数据更改为仅使用当前视图中的数据，以适应感兴趣区域。

• 在符号系统窗格中的统计数据中，对于统计数据，选择 DRA。
  
  

更改降雨量图层的符号系统

为了更好地查看降雨量图层值的变化，我们更改该图层的符号系统。

• 对于降雨量图层，在符号系统窗格中，将配色方案更改为降水量。
• 将统计信息更改为 DRA。
  
  地图的外观会发生变化。蓝色区域降雨最多，而红色区域最少。
  

以上通过应用定义查询、更改处理模板和修改符号系统，将所需的三个影像图层进行了可用于分析的处理。

构建栅格函数模板

创建一个栅格函数模板，以分析山体滑坡敏感性。该模板中将包含栅格函数，以对输入数据进行重分类并合并结果。该模板将生成一个输出图层，可用于标识更容易发生山体滑坡的区域。

栅格函数是一种更快、更轻量化的地理处理栅格工具替代方案。可将其动态应用于数据集的显示像素。栅格函数的输出是内存栅格图层。这样可以缩短处理时间，并免去冗余数据的创建和存储。

计算植被指数NDVI

将栅格函数合并到栅格函数模板中，就像 ModelBuilder 模型一样。

• 在影像选项卡的分析组中，单击函数编辑器。
  
  随即显示空白的函数编辑器视图。

• 在函数编辑器工具栏上，单击添加栅格变量按钮。
  
  一个绿色块随即显示在编辑器中。这是输入栅格的占位符。

• 右键单击栅格变量块，然后单击重命名。将现有文本替换为多光谱影像，然后按 Enter 键。

• 单击功能区上的影像选项卡。在分析组中，单击栅格函数按钮。
  随即显示栅格函数窗格。

• 在系统选项卡上的分析组中，找到 NDVI 函数并将其拖动到函数编辑器上。
  

• 指向多光谱影像占位符。单击OUT标注，然后将其拖动到 NDVI 函数上的栅格标注。
  这两个块随即会通过灰色箭头链接在一起。现在，将多光谱影像设置为 NDVI 函数的输入栅格。
  

• 双击 NDVI 函数以打开其属性窗口。
  NDVI 函数使用可见光的红色波长和近红外波长来计算绿色植被的密度。在多光谱 Landsat 影像图层中，这些波长是通过波段 4 和 5 进行测量的。
  关于Landsat 8的波段组合定义可以点击这里

• 在参数选项卡的可见波段 ID 中，输入 4。对于红外波段 ID，输入 5。

• 单击变量选项卡。在栅格行中，选中 IsPublic 框。
  公开栅格输入参数，以便在运行模板时可以选择不同的多光谱影像。用于具有不同数据的不同区域，而不仅限于针对此分析找到的特定 Landsat 场景。

• 单击确定关闭 NDVI 属性窗口。
  NDVI 函数将评估植被健康状况，但是要在山体滑坡敏感性分析中使用这些值，则需要将其重映射。
• 在栅格函数窗格中，搜索重映射，然后将重映射函数拖动到函数编辑器上。
• 从 NDVI 函数上的OUT标注拖动到重映射函数上的栅格标注以将两者进行连接。
  
• 双击重映射函数。
  在参数选项卡上，确保将重映射定义类型设置为列表，并根据下表更新列表：

最小值	最大值	输出
0	50	5
50	100	4
100	133	3
133	166	2
166	202	1

该操作会将 NDVI 值重分类或重映射为五个新的指数值。由于植被稀疏会导致更高的山体滑坡风险，因此最低的 NDVI 值（几乎没有绿色植被）将接收最高重映射值 5。

分类的值可以根据区域特点进行修改，推荐使用自然间断法把中断值求算出来，然后通过表导入到定义类型中，精度可能会更好一些。

• 单击常规选项卡。对于名称，输入 NDVI Remap。
• 单击确定关闭重映射属性窗口。
• 在函数编辑器工具栏上，单击保存按钮。
• 在保存窗口的名称中，键入山体滑坡敏感性分析。
• 对于描述，可以对该模型进行说明。
• 单击确定。
  模板的名称将在函数编辑器选项卡上更新。
  

计算坡度指数值

与植被指数类似，将坡度图层重映射为五个指数值。

• 将栅格变量添加到编辑器画布。

• 将栅格变量重命名为坡度。

• 从栅格函数窗格中，添加重映射函数并将坡度连接到该函数。
  

• 双击重映射函数以打开重映射属性窗口。

• 在参数选项卡上，针对重映射定义类型，选择列表，并根据下表更新列表：

最小值	最大值	输出
0	10	1
10	25	2
25	40	3
40	55	4
55	91	5

由于较陡的坡度更容易发生山体滑坡，因此将较高的坡度值映射到较高的指数值。

• 单击变量选项卡。对于栅格参数，选中 IsPublic 复选框。
• 单击常规选项卡。修改名称为Slope Remap。
• 单击确定。

计算降雨量指数值

与坡度指数值计算一样，将降雨量图层重映射为指数值。

• 将栅格变量添加到函数编辑器中，并将其重命名为降雨量。
• 添加重映射栅格函数，并将降雨量连接到该函数。
• 双击重映射以打开其属性窗口。
• 在参数选项卡上，针对重映射定义类型，选择列表，并根据下表更新列表：

最小值	最大值	输出
0	3	1
3	6	2
6	9	3
9	15	4
15	22	5

高日均降水量值将重映射为高指数值。

• 单击变量选项卡。对于栅格参数，选中 IsPublic 复选框。
• 单击常规选项卡。修改名称为Rainfall Remap。
• 单击确定。

合并指数值以生成敏感性表面

使用加权总和函数将其合并。此函数将对每个输入栅格进行加权，然后叠加所有输入并计算其值的总和。

• 在栅格函数窗格中，搜索加权总和函数并将其拖动到函数编辑器上。
  

• 将 NDVI Remap、Slope Remap 和 Rainfall Remap 连接到加权总和函数。

• 在函数编辑器工具栏上，单击自动布局按钮。
  
  接下来，需要为每个输入分配权重，以指示哪些变量应被视为分析中最重要的变量。

• 双击加权总和栅格函数。
  在烧毁区域，植被缺乏程度将替代降雨或坡度成为评估山体滑坡风险的更为重要的因素。针对此分析，赋予 NDVI 输入较高的权重，以指示植被指数对结果的影响是坡度或降雨输入影响的两倍。

• 在参数选项卡上的NDVI Remap.OutputRaster表行中，将权重值更改为 2。

• 对于像元大小类型，选择最小值。
  这将确保输出栅格采用最精细的输入图层的分辨率而非最粗糙的分辨率。
  

• 单击确定，然后保存栅格函数模板。

栅格函数模板的结果将是一个内存中的图层，其值的范围可能在 4 到 20 之间，较高的值表示山体滑坡敏感性更强。

指定栅格输出的位大小

由于输出值均为 4 到 20 之间的整数，结果可以用 8 位无符号数据类型表示。

-在栅格函数窗格中，搜索转为整型栅格函数并将其拖动到函数编辑器上。

• 加权总和函数拖动到转为整型函数以将两者进行连接。

• 双击转为整型栅格函数以打开其属性窗口。

• 单击常规选项卡。对于输出像素类型，选择 8 位无符号。
  

• 单击确定。

将结果裁剪至感兴趣区域

将结果裁剪至木里森林火灾边界范围。

• 在栅格函数窗格中，搜索裁剪栅格函数并将其拖动到函数编辑器上。
  裁剪函数具有两种可能的输入连接：一种用于输入栅格，一种用于裁剪几何。本次需将转为整型的输出连接到裁剪的栅格输入。

• 从转为整型函数拖动到裁剪函数上的栅格标注以将两者进行连接。
  

• 双击裁剪栅格函数以打开其属性窗口。

• 单击变量选项卡。对于 ClippingGeometry 和Extent参数，选中 IsPublic 复选框。
  

• 单击确定。

• 保存栅格函数模板并关闭函数编辑器。

• 保存工程。

执行山体滑坡敏感性分析

运行栅格函数模板

在运行函数模板前，先将处理范围设置到仅处理输入数据范围。避免无谓的计算。

• 在功能区上，单击分析选项卡。在地理处理组中，单击环境。
  

• 对于处理范围，选择木里328火灾，范围坐标将自动更新以匹配图层范围。之后在此工程中运行任何地理处理工具或栅格函数时，将应用此环境设置。
  并行处理因子根据自身机器性能选择，可以提高处理效率。
  

• 在栅格函数窗格中，单击之前保存的山体滑坡敏感性分析模型。

• 四个选项将显示在参数选项卡上。这些就是在构建栅格函数模板时设置为 IsPublic 的参数。依次选择对应的数据。
  

• 单击创建新图层按钮。
  栅格函数模板完成后，新的内存栅格图层将添加到地图。该图层包含感兴趣区域的山体滑坡敏感性值。该结果是通过动态处理在线数据创建的内存图层。
  

导出栅格

• 在内容窗格中，右键单击 Landslide Susceptibility Analysis，指向数据，然后选择导出栅格。
  导出栅格窗格随即出现。
• 像元大小默认设置为 0.25 米。对于如此大的区域，此分辨率很小。这里设置10，将其增大，以使栅格可以更快地导出。

修改栅格符合系统

导出的新图层随即会自动添加到地图。
图层具有黑色背景。可以通过符号系统，掩膜，勾选显示背景值去除黑色背景
图层以灰度绘制。白色地区发生山体滑坡的风险最高。

为了使其在视觉上更清晰，重新修改栅格符号系统。

• 在内容窗格中，右键单击 LandslideSusceptibility.tif，然后选择符号系统。
• 对于配色方案，选择 地狱红。
  当前，将以最浅的颜色显示山体滑坡敏感性最高的区域。
• 重新打开配色方案菜单，然后单击格式化配色方案。
• 在反向配色方案中，单击反向配色方案按钮，然后单击确定。

当前，将以最深的颜色显示山体滑坡敏感性最高的区域。

使用分析结果查找存在风险的位置

通过地理处理工具，可以将道路要素按等距离间隔分配敏感性值。由于没有该区域道路信息，因此选择创建随机点的方式来评估整个区域范围内的敏感性位置。

将山体滑坡敏感性值分配至随机点

• 在数据管理工具中，选择采样工具栏，选择创建随机点。本次创建2000个随机样点。
  

• 在地理处理窗格中，搜索并打开多值提取至点工具。
  该工具将为点要素添加新属性。此属性的值将等于每个点位置的山体滑坡敏感性分析输出栅格中的像元值。

• 对于输入点要素，选择Q2000。

• 对于输入栅格，选择Landslide Analysis.tif

• 对于输出字段名称，输入 Susceptibility_Score。
  

• 单击运行。每个点将根据基础栅格分配敏感性得分。
  

具有最高山体滑坡风险的位置

使用定义查询和分级符号来隔离山体滑坡风险最高的道路点。

• 在内容窗格中右键单击Q2000并选择属性。
• 单击定义查询选项卡，然后单击新建定义查询。
• 使用下拉菜单构建子句 Where Susceptibility_Score 大于 12。
  
  现在，该地图仅显示山体滑坡风险最高的道路点。
• 右键单击Q2000 图层，然后选择符号系统。
  分级符号将为具有较大敏感性得分的点绘制较大的符号。
  
• 切换为热力图显示效果如下
  
  缩放到局部。
  

森林火灾的负面影响之一是破坏了健康的绿色植被。当大雨落在已被砍掉植被的陡坡上时，即存在很高的山体滑坡风险。基于植被密度、坡度和降水量，对木里森林火灾地区的山体滑坡敏感性分析迅速确定了具有较高山体滑坡风险的潜在区域。后续可根据地质灾害风险相关模型丰富参数，研究结果可用于主动计划最关键地区的山体滑坡减灾工作。

鸣谢

本文参考了Mark Gilbert和Aileen Buckley的相关文章，使用的所有数据均来自 ArcGIS Living Atlas of the World 、ArcGIS Online和NOAA。

• 多光谱 Landsat 卫星：Esri、USGS、AWS、NASA
• 地形：Airbus、USGS、NGA、NASA、CGIAR、NLS、OS、NMA、Geodatastyrelsen、GSA、GSI 以及 GIS 用户社区
• 降雨量：NOAA

更多精彩内容请关注下方公众号