数字地形模型上机实习二.docx

数字地形模型课程实习2
栅格数据
栅格数据是通过指定栅格类型的方式添加到镶嵌数据集中的。栅格类型用于与栅格格式一起标识元数据，例如地理配准、采集日期和传感器类型。栅格格式用于定义像素的存储方式，例如，行数和列数、波段数、实际像素值，以及其他栅格格式特定的参数。但是，根据栅格类型添加栅格数据时，会读取相应的元数据并将其用于定义任何需要应用的处理。
Esri Grid 格式
格网是 Esri 栅格数据的原生存储格式。通常包含以下两种类型的格网：整型和浮点型。整型格网多用于表示离散数据，浮点型格网则多用于表示连续数据。
整型格网的属性存储在它的值属性表 (VAT) 中。格网中的每个唯一值对应于表中的一条 VAT 记录。该记录存储了这个唯一值（VALUE 是表示特定类或像元分组的整数）和它所表示的格网像元数 (COUNT)。例如，如果栅格中共有 50 个代表森林的值是 1 的像元，则在 VAT 中,这些像元将显示为一条 VALUE = 1 和 COUNT = 50 的记录。

图1 Esri Grid 格式
浮点型格网没有 VAT，因为格网中的像元可以是给定范围的任意值。此格网类型中的像元不能整齐地落在各个离散类别中。像元值用于描述其所在位置的属性。例如，在使用米作为单位的表示高程的高程数据格网中，像元值 10.1662 代表其所在位置高于海平面大约 10 米。
可存储为格网值的数据值的范围如下：浮点型格网可存储值的范围为 -3.4 x 1038 至 3.4 x 1038。整型格网可存储值的范围为 -2147483648 至 2147483647（-231 至 231-1）。
对于整型格网，此信息仅适用于 VALUE 项。整型格网可将其他 INFO 项添加到值范围取决于该项定义的 VAT。
格网的坐标系与其他地理数据的坐标系相同。行和列分别与坐标系的 x 和 y 轴平行。由于格网中的每个像元都与其他像元具有相同的尺寸，因此通过行和列可轻松地确定任意像元的位置和其所覆盖的区域。这样，格网的坐标系就可由像元大小、行和列的数目，以及左上角的 x,y 坐标定义。格网也可传递其他信息，例如与格网相关联的坐标系。

不同格式的比较
格式
描述
扩展模块
读/写
ESRI Grid
支持 32 位整型和 32 位浮点型栅格格网的专有 ESRI 格式。格网可用于表示在空间上连续变化的地理现象，及用于执行流向、趋势和表面（如水文）的空间建模和空间分析。
目录
色彩映射表文件 - 扩展名 *.clr
可读写
标记图像文件格式 (TIFF)（支持 GeoTIFF 标记）
广泛应用于桌面出版领域。它可用作某些扫描仪和图形艺术包的界面。TIFF 支持黑白、灰度、伪彩色以及真彩色图像，所有这些图像都可以以压缩或者解压缩的格式存储。
支持 BigTIFF 格式。
单个文件 - 可能的文件扩展名 *.tif、*.tiff 和 *.tff
坐标文件 - 扩展名 *.tfw
默认情况下，ArcCatalog 只能识别 .tif 文件扩展名。要在不重新命名 .tiff 或 .tff 文件的情况下将这两种文件添加到 ArcMap，可将这些文件扩展名添加到 ArcCatalog，或将这些文件从 Windows 资源管理器拖入地图。
可读写
ENVI 文件头格式
使用 ENVI处理栅格数据集时，会创建一个包含软件需求信息的头文件。可针对多个栅格文件格式创建此类头文件。
头文件 - 扩展名 *.hdr
多个数据文件 - 扩展名为 *.raw、*.img、*.dat、*.bsq 等
可读写
（.dat - 通过 UI）
（.bsq、.img 和 .raw - 仅开发人员）

不同类型的数字高程模型
数字高程模型(Digital Elevation Model，DEM)，是国家基础空间数据的重要组成部分，它表示地表区域上地形的三维向量的有限序列，即地表单元上高程的集合，数学表达为：z=f(x，y)。
Contour
等高线通常被存储成一个有序的坐标点序列，可以认为是一条带有高程值属性的简单多边形或多边形弧段。由于等高线模型只是表达了区域的部分高程值，往往需要一种插值方法来计算落在等高线以外的其他点的高程，又因为这些点是落在两条等高线包围的区域内，所以，通常只要使用外包的两条等高线的高程进行插值。
缺点：
1. 数字化现有等高线地图产生的DEM比直接利用航空摄影测量方法产生的DEM质量要差；
2. 数字化的等高线对于计算坡度或生成着色地形图不十分适用。

TIN to Contour
1. 等高线追踪，利用TIN点的高程内插出格网边上的等高线点，并将这些等高线点排序；
2. 等高线光滑，进一步加密等高线点并绘制光滑曲线。

图2 某地区Mo的含量
与TIN点转为等高线相似，采样点内插得出网格点，再连接等值点，光滑后得到等值面

Grid to Contour
1. 等高线追踪，将网格中的等高线点排序；
2. 等高线光滑，进一步加密等高线点并绘制光滑曲线。

Grid
规则格网法是把DEM表示成高程矩阵，此时，DEM来源于直接规则矩形格网采样点或由不规则离散数据点内插产生。它结构简单，计算机对矩阵的处理比较方便，高程矩阵已成为DEM最通用的形式。高程矩阵特别有利于各种应用。
缺点：
1. 地形简单的地区存在大量冗余数据；
2. 如不改变格网大小,则无法适用于起伏程度不同的地区；
3. 对于某些特殊计算如视线计算时，格网的轴线方向被夸大；
4. 由于栅格过于粗略，不能精确表示地形的关键特征,如山峰、洼坑、山脊等；
优点：计算机处理以栅格为基础的矩阵很方便，使高程矩阵成为最常见的DEM；
TIN to Grid
利用TIN点的高程内插出格网边上的点的高程
Contour to Grid
利用等高线内插出格网边上的点的高程

TIN
TIN（Triangulated Irregular Network)表示法利用所有采样点取得的离散数据，按照优化组合的原则，把这些离散点(各三角形的顶点)连接成相互连续的三角面(在连接时，尽可能地确保每个三角形都是锐角三角形或是三边的长度近似相等—Delaunay)。因为TIN可根据地形的复杂程度来确定采样点的密度和位置，能充分表示地形特征点和线，从而减少了地形较平坦地区的数据冗余。
表示方法：将区域划分为相邻的三角面网络，区域中任意点都将落在三角面顶点、线或三角形内。落在顶点上其高程与顶点相同；落在线上则由两个顶点线性插值得到；落在三角形内则由三个顶点插值得到。生成方法：由不规则点、矩形格网或等高线转换而得到。
TIN允许在地形复杂地区收集较多的信息，而在简单的地区收集少量信息，避免数据冗余。对于某些类型的运算比建立在数字等高线基础上的系统更有效，如坡度、坡向等的计算。
优缺点：不规则三角网(TIN)表示法克服了高程矩阵中冗余数据的问题，而且能更加有效地用于各类以DTM为基础的计算；但其结构复杂。

Grid to TIN
利用网格的高程内插出TIN点的高程
Contour to TIN
利用等高线的高程内插出TIN点的高程

插值方法原理
三次样条插值
分段线性插值的优点：计算简单、稳定性好、收敛性有保证且易在计算机上实现
缺点：它只能保证各小段曲线在连接点的连续性，却无法保证整条曲线的光滑性，这就不能满足某些工程技术的要求。
三次Hermit插值优点：有较好的光滑性，缺点：要求节点的一阶导数已知。
从２０世纪６０年代开始，首先由于航空、造船等工程设计的需要而发展起来所谓样条(Spline)插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越多广泛应用。
三次样条插值函数的定义：
给定区间上的个节点和这些点上的函数值 若满足：
（1）；
（2）在每个小区间上至多是一个三次多项式；
（3）在上连续。
则称为函数关于节点的三次样条插值函数。
边界问题的提出与类型
单靠一个函数表是不能完全构造出一个三次样条插值函数。我们分析一下其条件个数，条件（2）三次样条插值函数是一个分段三次多项式，若用表示它在第个子区间上的表达式，则形如

其中有四个待定系数，子区间共有个，所以共有个待定系数。
由条件（3）在上连续，即它们在各个子区间上的连接点上连续即可，共有个条件，即

共有个条件，未知量的个数是个。
这样需要加2个条件。
这两个条件通常在插值区间的边界点处给出，称为边界条件。
边界条件的类型很多，常见有：
（1） 给定一阶导数值；
（2）给定一阶导数值；（特别地时称为自然边界条件，满足自然边界条件的次样条插值函数称为自然样条插值函数）
（3） 当是周期为的函数时，要求及其导数都是以为周期的函数，相应的边界条件为：

注：1. 由是周期函数，必有。
2. 虽然可利用边界条件及方程组可求出待定系数，从而得到三次样条插值函数在各个子区间的表达式，但计算量大，不便于应用，我们介绍一种简便的方法（方法或三弯方法）：

设
因为在子区间上是不高于三次的多项式，其二阶导数必是线性函数，所以
(1)
设，则有
(2)
连续二次积分得：
(3)
其中是积分常数，
利用插值条件 ，带入（3）中得：

，
带入（3），整理得：

(4)
从而只要确定这个未知数，即可定出三次样条插值函数。
根据的连续条件，即：
等同于: , （*）
由（4）得
（5）
所以 （6）
将（5）式中的，得在子区间上的表达式为，
（7）
将（5）（7）代入（*）式中，整理得：
，
两编同乘得
（8）
令（ （9）
则所得方程可简写为:

即 (10)
在（1）边界条件下，即，
由于在子区间上的导数为

由 得 （11）
同理由 得 （12）
（11）（12）与（10）合在一起得下列的关于的线性方程组
（13）
其中 （14）
在（2）边界条件下，即，在方程中实际上只包含有未知数，并且可以改写为: (15)
在（3）边界条件下，
由得 （15）
由

注意到：，将上式整理得：

记
即 （16）
结合（15）（16）（10）得的线性方程组为：

由代数知识知：三对角方程组由于主对角占优，方程有唯一解，具体方法一般应用三对角的“追赶法”求解。
参数设置
Parameter
Explanation
Data Type
in_point_features
The input point features containing the z-values to be interpolated into a surface raster.
Feature Layer
z_field
The field that holds a height or magnitude value for each point.
This can be a numeric field or the Shape field if the input point features contain z-values.
Field
cell_size
(Optional)
The cell size at which the output raster will be created.
This will be the value in the environment if it is explicitly set; otherwise, it is the shorter of the width or the height of the extent of the input point features, in the input spatial reference, divided by 250.
Analysis Cell Size
spline_type
(Optional)
The type of spline to be used.
· REGULARIZED — Yields a smooth surface and smooth first derivatives.
· TENSION — Tunes the stiffness of the interpolant according to the character of the modeled phenomenon.
String
weight
(Optional)
Parameter influencing the character of the surface interpolation.
When the REGULARIZED option is used, it defines the weight of the third derivatives of the surface in the curvature minimization expression. If the TENSION option is used, it defines the weight of tension.
The default weight is 0.1.
Double
number_points
(Optional)
The number of points per region used for local approximation.
The default is 12.
Long

克里金法
克里金法假定采样点之间的距离或方向可以反映可用于说明表面变化的空间相关性。克里金法工 具可将数学函数与指定数量的点或指定半径内的所有点进行拟合以确定每个位置的输出值。克里金法是一个多步过程；它包括数据的探索性统计分析、变异函数建模 和创建表面，还包括研究方差表面。当您了解数据中存在空间相关距离或方向偏差后，便会认为克里金法是最适合的方法。该方法通常用在土壤科学和地质中。
由于克里金法可对周围的测量值进行加权以得出未测量位置的预测，因此它与反距离权重法类似。这两种插值器的常用公式均由数据的加权总和组成：
(17)

其中：Z(si) = 第 i 个位置处的测量值, λi = 第 i 个位置处的测量值的未知权重, s0 = 预测位置, N = 测量值数
在反距离权重法中，权重 λi 仅取决于预测位置的距离。但是，使用克里金方法时，权重不仅取决于测量点之间的距离、预测位置，还取决于基于测量点的整体空间排列。要在权重中使用空间排列，必须量化空间自相关。因此，在普通克里金法中，权重 λi 取决于测量点、预测位置的距离和预测位置周围的测量值之间空间关系的拟合模型。以下部分将讨论如何使用常用克里金法公式创建预测表面地图和预测准确性地图。
参数设置
Parameter
Explanation
Data Type
in_point_features
The input point features containing the z-values to be interpolated into a surface raster.
Feature Layer
z_field
The field that holds a height or magnitude value for each point.
This can be a numeric field or the Shape field if the input point features contain z-values.
Field
semiVariogram_props
kriging_model
The KrigingModel class defines which kriging model is to be used.
There are two types of kriging classes. The KrigingModelOrdinary method has five types of semivariograms available. The KrigingModelUniversal method has two types of semivariograms available.
· KrigingModelOrdinary ({semivariogramType}, {lagSize}, {majorRange}, {partialSill}, {nugget})
§ semivariogramType—The semivariogram model to be used. The available models include the following:
o SPHERICAL—Spherical semivariogram model. This is the default.
o CIRCULAR—Circular semivariogram model.
o EXPONENTIAL—Exponential semivariogram model.
o GAUSSIAN—Gaussian (or normal distribution) semivariogram model.
o LINEAR—Linear semivariogram model with a sill.
· KrigingModelUniversal ({semivariogramType}, {lagSize}, {majorRange}, {partialSill}, {nugget})
§ semivariogramType—The semivariogram model to be used. The available models include the following:
o LINEARDRIFT—Universal Kriging with linear drift.
o QUADRATICDRIFT—Universal Kriging with quadratic drift.
· After the {semivariogramType}, the other parameters are common between Ordinary and Universal kriging.
§ lagSize—The default is the output raster cell size.
§ majorRange—Represents a distance beyond which there is little or no correlation.
§ partialSill—The difference between the nugget and the sill.
§ nugget—Represents the error and variation at spatial scales too fine to detect. The nugget effect is seen as a discontinuity at the origin.
KrigingModel
cell_size
(Optional)
The cell size at which the output raster will be created.
This will be the value in the environment if it is explicitly set; otherwise, it is the shorter of the width or the height of the extent of the input point features, in the input spatial reference, divided by 250.
Analysis Cell Size
search_radius
(Optional)
The Radius class defines which of the input points will be used to interpolate the value for each cell in the output raster.
There are two types of radius classes: RadiusVariable and RadiusFixed. A Variable search radius is used to find a specified number of input sample points for the interpolation. The Fixed type uses a specified fixed distance within which all input points will be used for the interpolation. The Variable type is the default.
· RadiusVariable ({numberofPoints}, {maxDistance})
§ {numberofPoints}—An integer value specifying the number of nearest input sample points to be used to perform interpolation. The default is 12 points.
§ {maxDistance}—Specifies the distance, in map units, by which to limit the search for the nearest input sample points. The default value is the length of the extent's diagonal.
· RadiusFixed ({distance}, {minNumberofPoints})
§ {distance}—Specifies the distance as a radius within which input sample points will be used to perform the interpolation.
The value of the radius is expressed in map units. The default radius is five times the cell size of the output raster.
§ {minNumberofPoints}—An integer defining the minimum number of points to be used for interpolation. The default value is 0.
If the required number of points is not found within the specified distance, the search distance will be increased until the specified minimum number of points is found.
When the search radius needs to be increased it is done so until the {minNumberofPoints} fall within that radius, or the extent of the radius crosses the lower (southern) and/or upper (northern) extent of the output raster. NoData is assigned to all locations that do not satisfy the above condition.
Radius
out_variance_prediction_raster
(Optional)
Optional output raster where each cell contains the predicted semi-variance values for that location.
Raster Dataset

反距离权重法
反距离权重法主要依赖于反距离的幂值，幂参数可基于距输出点的距离来控制已知点对内插值的影响。幂参数是一个正实数，默认值为2。（一般0.5到3的值可获得最合理的结果）。
通过定义更高的幂值，可进一步强调最近点。因此，邻近数据将受到更大影响，表面会变得更加详细（更不平滑）。随着幂数的增大，内插值将逐渐接近最近采样点的值。指定较小的幂值将对距离较远的周围点产生更大的影响，从而导致平面更加平滑。
由于反距离权重公式与任何实际的物理过程都不关联，因此无法确定特定幂值是否过大。作为常规准则，认为值为30的幂是超大幂，因此不建议使用。此外还要牢记一点，如果距离或幂值较大，则可能生成错误结果。

在IDW插值之前,我们可以事先获取一个离散点子集，用于计算插值的权重；
原因1：离散点距离插值点越远，其对插值点的影响力越低，甚至完全没有影响力；
原因2：离散点越少可以加快运算速度；

参数设置

Parameter
Explanation
Data Type
in_point_features
The input point features containing the z-values to be interpolated into a surface raster.
Feature Layer
z_field
The field that holds a height or magnitude value for each point.
This can be a numeric field or the Shape field if the input point features contain z-values.
Field
cell_size
(Optional)
The cell size at which the output raster will be created.
This will be the value in the environment if it is explicitly set; otherwise, it is the shorter of the width or the height of the extent of the input point features, in the input spatial reference, divided by 250.
Analysis Cell Size
power
(Optional)
The exponent of distance.
Controls the significance of surrounding points on the interpolated value. A higher power results in less influence from distant points. It can be any real number greater than 0, but the most reasonable results will be obtained using values from 0.5 to 3. The default is 2.
Double
search_radius
(Optional)
The Radius class defines which of the input points will be used to interpolate the value for each cell in the output raster.
There are two types of radius classes: RadiusVariable and RadiusFixed. A Variable search radius is used to find a specified number of input sample points for the interpolation. The Fixed type uses a specified fixed distance within which all input points will be used for the interpolation. The Variable type is the default.
RadiusVariable ({numberofPoints}, {maxDistance})
§ {numberofPoints}—An integer value specifying the number of nearest input sample points to be used to perform interpolation. The default is 12 points.
§ {maxDistance}—Specifies the distance, in map units, by which to limit the search for the nearest input sample points. The default value is the length of the extent's diagonal.
· RadiusFixed ({distance}, {minNumberofPoints})
§ {distance}—Specifies the distance as a radius within which input sample points will be used to perform the interpolation.
The value of the radius is expressed in map units. The default radius is five times the cell size of the output raster.
§ {minNumberofPoints}—An integer defining the minimum number of points to be used for interpolation. The default value is 0.
If the required number of points is not found within the specified distance, the search distance will be increased until the specified minimum number of points is found.
When the search radius needs to be increased it is done so until the {minNumberofPoints} fall within that radius, or the extent of the radius crosses the lower (southern) and/or upper (northern) extent of the output raster. NoData is assigned to all locations that do not satisfy the above condition.
Radius
in_barrier_polyline_features
(Optional)
Polyline features to be used as a break or limit in searching for the input sample points.
Feature Layer

结果

(a)

(b)

(c)

(d)
图3
(a)envi格式的dem (b)grid (c)TIN (d)等高线

(a)

(b)

(c)

(d)
图3
(a)TIN点 (b)克里金法插值 (c)反距离权重法插值 (d)三次样条插值