从大型激光雷达点集合创建栅格 DEM 和 DSM
摘要
栅格（或格网化）高程模型是最常用的 GIS 数据类型之一。它们可用于多种分析方式并且可以很容易地共享。凭借激光雷达，您能够制作出两种不同风格的高质量高程模型：第一次回波和地面。第一次回波表面包括树冠和建筑物，通常被称为数字表面模型 (DSM)。地面（或裸露地表）仅包含地形，通常被称为数字高程模型 (DEM)。
左侧的图像显示了使用山体阴影表示的第一次回波表面（或 DSM），右侧的图像显示的是裸露地表模型（或 DEM）。
 |  |
制定计划
在用激光雷达创建栅格之前，应对一些基本因素进行评估：
- 激光雷达 coverage 的范围
- 激光雷达点的数量和点密度
- 所需的输出栅格分辨率
- 输出栅格的范围
- 输出栅格的格式
考虑这些因素将有助于决定是生成一个栅格还是一组栅格。作为此过程的一部分，您需要算出一个栅格中要具有多少行和列。而这取决于在数据的分析、显示和可能的共享或分发方面所希望实现的效果。但是，仅生成一个数据集来进行分析可能不可行，因为可能会与数据集大小方面的实际限制相冲突。因此，另一个需要考虑的问题是您拥有的激光雷达数据量。如果试图以一个数据集来处理 100 亿个激光雷达点，尽管有可能，但使用起来会非常不便。在此情形下，最好是将如此大量的激光雷达数据制成多个栅格，同时也应考虑将激光雷达处理过程拆分开来。这不仅可使单个数据集的大小合理，同时还可缩短这些数据集的处理时间。处理过程的时间越长，便越有可能遇到某些意外（例如断电）。
如果已决定要分割数据，则接下来的问题是操作如何进行。是基于常规的格网系统、行政边界分割，还是基于预期的应用分割？由于激光雷达点集通常都有多种用途，所以按照常规的格网系统或行政边界（例如县边界）分割最有意义。这样便可在具体项目中根据所需镶嵌不同的部分。如果所需的用途与某种应用关系非常密切（例如水文），则可采用该应用的分割逻辑。例如，如果是作水文用途，则流域边界是很好的备选方案。
ArcGIS 可支持多种栅格格式，所以您可选择要写出的格式。最好基于产品所需用途来做出决定。如果要与普通公众共享，则可考虑使用 TIFF 或 JPEG 格式来分发。对于使用 ArcGIS 平台的分析，则可考虑使用基于文件的地理数据库格式。
用激光雷达点来制作栅格的第一步是将这些点加载到地理数据库中。要将激光雷达点加载到多点要素类中，应使用 LAS 转多点或 ASCII 3D 文件转要素类地理处理工具，具体则取决于激光雷达数据的源数据格式。如果要用激光雷达点来构建 terrain 数据集，则应将多点要素类放到要素数据集中。尽管可选择使用 LAS 或 ASCII 格式的文件，但二进制文件格式的 LAS 的可接受度更高。LAS 文件包含的信息更多，而且由于采用二进制，所以可以被导入程序更有效的读取。有关将激光雷达源测量值导入地理数据库的详细信息，请参阅数据导入和加载工具。
使用点转栅格地理处理工具
如果激光雷达数据是唯一的数据源，则可使用点转栅格地理处理工具来生成栅格高程模型。使用点转栅格工具所生成结果的质量并不是最高的。但因为激光雷达点通常都非常密集，所以在很多应用中使用点转栅格地理处理工具所生成的精度足以满足需求。同时此工具还兼具方便和快速的特性，从而使它成为最称手的工具。
如果要制作一个裸露地表栅格表面（或 DEM），可以仅使用地面激光雷达点来生成栅格。将该工具的值字段参数设为 Shape，以使用多点折点的 z 值。而且，将单元分配类型设置为 MIN 或 MEAN。MIN 将使输出的高度偏向局部低点；而 MEAN 则是更适合常规用途的选项。要制作第一次回波表面（或 DSM），则可使用该工具的 MAX 选项来处理第一次回波激光雷达点，因为所希望的是使输出偏向于局部高点。
点转栅格地理处理工具可利用激光雷达点集来生成格网化高程模型。
尽管点转栅格是利用激光雷达点制作栅格的最简便且最快速的方式，但它有一个严重的缺点。它最终可能会产生许多 NoData 单元，因为只有当单元中有一个或多个点时，它才会返回具体值。因为地面上被植被和建筑物遮挡的位置会使数据产生空隙，所以如果仅使用地面点来制作 DEM，此问题会更加严重。要减小 NoData 对数据单元的影响，可增加输出单元大小相对于平均点间距的比例。以下情况下，还可通过在 Python 窗口中使用下列示例，在执行点转栅格后减少 NoData 像元数：
- 已启用 Spatial Analyst 扩展模块
- 上一次执行的焦点统计表达式的输出名称为“outfocalmNN”（outfocalm03、outfocalm04、outfocalm05 等）
- 输出（在本例中，点转栅格的输出称为 point2ras）为内容列表中的一个图层
- 有效路径为：C:\data\
from arcpy.sa import *
outfocalm01 = FocalStatistics("point2ras", NbrRectangle(3, 3, "CELL")
outfocalm02 = FocalStatistics(outfocalm01, NbrRectangle(3, 3, "CELL")
#Repeat using output (temporary raster object) in the next processes until
all nodata is gone
out = Con(IsNull("point2ras"), outfocalmNN, "point2ras")
#Save the result to disk
out.save("C:\data\myfinalDEM")
可多次运行 Python 示例来填充较大的 NoData 区域，但是建议不要超过两次。对于较大的空区域，最好还是接受它的存在，这是使用此方法的必然后果。
 |  |  |
左图是原始的点，中间的是点转栅格的输出（NoData 单元为白色），右侧的是经过后处理的栅格，其中的 NoData 已经填充。
使用 terrain 数据集来创建栅格 DEM
如果具有与激光雷达点配套的摄影测量隔断线，或需要使结果的质量比点转栅格工具所能实现的更高，则应使用 terrain 数据集。有关 terrain 数据集的概述，请参阅什么是 terrain 数据集？
 |
|
左侧的表面中没有沿河岸的隔断线。右侧的则进行了隔断线强化。隔断线对于在高程模型中保证与水相关要素的清晰界定非常重要。
隔断线的用途是捕捉表面中的线性不连续的地形。最常见的类型是路面边缘、湖岸线、小河流的单线河道和大河流的双线河道。有时，为对表面进行定义和塑型，还可通过采集隔断线的方式避免显示不连续的地形。这些应用的示例包括：形式类似等值线的线和平缓山脊的顶峰。
隔断线经常在裸露地表模型中使用，但不适合用于第一次回波表面中，因为它们可能会与地上点冲突。例如，对于捕捉路面边缘的各隔断线，其 x、y 可能与突出于路面的树冠对应点的 x、y 相同，但在 z 上不同。正因如此，应考虑排除第一次回波表面的隔断线，或至少考虑排除您已知道可能会产生冲突的隔断线。
对用于 terrain 数据集（请参阅下表）的隔断线进行组织的最有效方式是，按照表面要素类型 (SFType) 将其分配到不同的要素类中。表面要素类型不但控制着在模型中要素强化的方式，还控制着（在栅格化过程中使用的）自然邻域插值器在遇到这些要素时解释表面的方式。在遇到“硬”要素时将出现坡度明显中断的情况，但在遇到“软”要素时则不会。
|
测量类型 |
要素类类型 |
SFType |
|---|---|---|
|
激光雷达点 |
3D 多点要素类 |
离散多点 |
|
路面边缘、单线和双线河道、陡峭的山脊线 |
3D 线要素类 |
硬断线 |
|
湖泊、水库 |
将 z 值存储为属性的 2D 面要素类 |
硬断线或硬替换 |
|
被侵蚀的/平缓的山脊线、形式如等值线的线 |
3D 线要素类 |
软断线 |
|
研究区域边界 |
无 z 值的 2D 面要素类 |
软裁剪 |
出于 terrain 性能的考虑，最好将所有硬断线都放置在一个要素类中。但有时这一要求无法实现，例如，需要分别保留道路和水体要素。请记住，定义 terrain 时，使用的要素类越少越好。
“替换 SFType”用于将面中的一切都设为恒定高度。它最常应用的对象是湖泊，因为湖泊中可能会不慎加入其他数据（例如激光雷达点），而如果这些数据的高度与湖岸线的高度不完全相同，则可能使水体无法保持水平。与常规硬断线和软断线相比，“替换 SFType”确实会带来更大的处理成本，所以最好避免在 terrain 数据集中使用它。理想情况下，水体中不应存在激光雷达样本（可考虑将此作为一项约定添加到与数据提供商的合同中），但是如果存在，则可在 terrain 数据集构建后使用删除 Terrain 点地理处理工具来对其进行处理。或者，可在构建 terrain 前使用擦除点地理处理工具来消除所有不恰当的点。
如果要使用 terrain 数据集同时生成裸露地表和第一次回波表面，则可将激光雷达点载入到两个不同的多点要素类，其中一个要素类为地面点，另一个为地上点。定义裸露地表 terrain 时仅引用地面点。第一次回波 terrain 数据集则不仅可引用裸露地表 terrain 所用的地面点要素类，同时还可引用地上点。这意味着两个不同的 terrain 数据集可引用同一个要素类。
从 ArcGIS 9.3 开始，便可通过两种点细化过滤器中的一种来对 terrain 数据集进行构建金字塔操作：z 容差和窗口大小。如要生成 DEM，则使用两种金字塔类型均可。如果要对全分辨率点集进行栅格化，则可使用窗口大小过滤器来构建 terrain，因为这样速度会快很多。如果想要使用已细化的数据进行分析（如果激光雷达数据的采样超出了您的需求，此情况是很合理的），则可使用 z 容差过滤器。尽管更加耗时，但是这仍是最适合的方式，因为这样会为已细化的制图表达估计垂直精度。如要生成 DSM，则应使用窗口大小过滤器，并应用 MAX 选项。
使用 Terrain 转栅格工具来制作栅格化高程模型。可选择插值、输出单元大小和使用 terrain 数据集的哪种金字塔等级。
Terrain 转栅格工具可利用 terrain 数据集来生成格网化高程模型。
对于插值，自然邻域法选项是最佳选项。尽管速度不及线性插值，但其生成的结果通常更加美观精确。相对于激光雷达点采样密度设置输出单元大小。使用远小于平均点间距的单元大小并不会对精确性有任何帮助。另外，还应确保在适当的情况下为子集提取设置分析范围（如同使用环境进行设置）。因为可以对齐栅格输出，所以捕捉栅格也十分有用。
下面介绍了在 ArcGIS 中利用激光雷达点数据生成栅格 DEM 表面的主要步骤。
首先，应在您应用程序的目录 窗口或 ArcCatalog 中创建 terrain 数据集，然后使用地理处理工具来将此 terrain 数据集转换为栅格 DEM。
1. 创建 terrain 数据集。
-
确定源数据和使用该数据构成 terrain 数据集的方式。
有关表示 terrain 源数据的更多信息，请参阅呈现要素类中的 terrain 源数据和 terrain 数据集中支持的源数据类型。
注：在构建 terrain 数据集时所有源数据必须具有相同的空间参考。
- 在 ArcCatalog 或目录 窗口中创建文件地理数据库。右键单击要构建 terrain 的文件夹，指向新建，然后在快捷菜单中单击文件地理数据库。
-
创建要素数据集。右键单击文件地理数据库，指向新建，然后在快捷菜单中单击要素数据集。
有关正确生成要素数据集的更多信息，请参阅创建要素数据集。
- 将源测量值导入要素类。这些要素类必须在步骤 3 中创建的要素数据集内生成。有关如何为 terrain 导入源数据的详细信息，请参阅导入 terrain 数据集源测量值。
-
在 ArcCatalog 中或目录 窗口中使用新建 Terrain 向导构建 terrain 数据集。
要访问新建 Terrain 向导，右键单击要素数据集以显示菜单，指向新建，然后单击 Terrain。
有关使用新建 Terrain 向导的详细信息，请参阅使用 Terrain 向导构建 terrain 数据集。
2. 使用“Terrain 转栅格”地理处理工具。
- 在 3D Analyst 工具中，右键单击 Terrain 转栅格地理处理工具将其打开。
- 单击输入 Terrain 浏览按钮来添加 terrain 数据集。
- 单击输出栅格浏览按钮来指定要创建的栅格数据集的位置。
- 将可选输出数据类型参数设置为 32 位浮点型或 32 位整型。默认设置为浮点型。
-
将插值方法设置为线性或自然邻域。
这些方法是在三角化 terrain 表面上应用的基于 TIN 的插值方法。线性选项可找到包围每个单元中心的三角形，并会应用三角形的节点的加权平均值来执行插值计算。自然邻域选项使用单元中心的 Voronoi 邻域。
在插入 terrain 表面时，可考虑自然邻域方法。虽然自然邻域插值的处理时间更长，但其生成的表面比使用线性插值生成的表面平滑得多。而且它更不容易受到三角测量微小变化的影响。
-
将采样距离设置为观测值或单元大小，这控制着栅格的水平分辨率。在选定所需方法后，指定该选项旁的值。
观测值方法会基于此数据表示的设定值和在栅格表面最长的边缘上需要的单元数量来计算单元大小。可以使用单元大小选项直接设置单元大小。
- 设置要使用的分辨率。分辨率参数表示用于转换的 terrain 数据集的金字塔等级。要以全分辨率输出栅格数据集，需将此参数设置为 0。金字塔等级使用 z 容差或窗口大小确定，这表示 terrain 数据集相对于全分辨率数据的近似分辨率。
- 请考虑使用环境设置来明确地控制要生成的 DEM 的范围。要提取 terrain 的子集，可在该地理处理工具的底部处单击环境按钮。单击常规设置选项卡并定义输出 DEM 的范围。
不论使用点转栅格还是 Terrain 转栅格地理处理工具，您都可将数亿甚至是数十亿的激光雷达点处理成高分辨率的格网化 DEM 和 DSM。然后，这些表面模型便可供 ArcGIS 中大量的栅格工具进行分析。
它们还非常适合于制作地图（请参阅下图），而且由于它们的数据结构简单，所以也便于共享。
