数字地形模型上机实习二分解

数字地形模型课程实习2

栅格数据

栅格数据是通过指定栅格类型的方式添加到镶嵌数据集中的。栅格类型用于与栅格格式一起标识元数据，例如地理配准、采集日期和传感器类型。栅格格式用于定义像素的存储方式，例如，行数和列数、波段数、实际像素值，以及其他栅格格式特定的参数。但是，根据栅格类型添加栅格数据时，会读取相应的元数据并将其用于定义任何需要应用的处理。

Esri Grid 格式

格网是Esri 栅格数据的原生存储格式。通常包含以下两种类型的格网：整型和浮点型。整型格网多用于表示离散数据，浮点型格网则多用于表示连续数据。

整型格网的属性存储在它的值属性表(V A T) 中。格网中的每个唯一值对应于表中的一条V AT 记录。该记录存储了这个唯一值（V ALUE 是表示特定类或像元分组的整数）和它所表示的格网像元数(COUNT)。例如，如果栅格中共有50 个代表森林的值是 1 的像元，则在V AT 中,这些像元将显示为一条V ALUE = 1 和COUNT = 50 的记录。

图1 Esri Grid 格式

浮点型格网没有VAT，因为格网中的像元可以是给定范围的任意值。此格网类型中的像元不能整齐地落在各个离散类别中。像元值用于描述其所在位置的属性。例如，在使用米作为单位的表示高程的高程数据格网中，像元值10.1662 代表其所在位置高于海平面大约10 米。

可存储为格网值的数据值的范围如下：浮点型格网可存储值的范围为-3.4 x 1038至3.4 x 1038。整型格网可存储值的范围为-2147483648 至2147483647（-231至231-1）。

对于整型格网，此信息仅适用于V ALUE 项。整型格网可将其他INFO 项添加到值范围取决于该项定义的VAT。

格网的坐标系与其他地理数据的坐标系相同。行和列分别与坐标系的x 和y 轴平行。由于格网中的每个像元都与其他像元具有相同的尺寸，因此通过行和列可轻松地确定任意像元的位置和其所覆盖的区域。这样，格网的坐标系就可由像元大小、行和列的数目，以及左上角的x,y 坐标定义。格网也可传递其他信息，例如与格网相关联的坐标系。

不同格式的比较

不同类型的数字高程模型

数字高程模型(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)插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越多广

泛应用。

三次样条插值函数的定义：

给定区间],[b a 上的个节点b x x x a n =<<<= 10和这些点上的函数值

),,1,0()(n i y x f i i == 若)(x S 满足：

（1）),,2,1,0()

(n i y x S i i ==；

（2）在每个小区间],[b a 上至多是一个三次多项式； （3）)(),(),(x S x S x S '''在],[b a 上连续。

则称)(x S 为函数)(x f 关于节点的n x x x ,,,10 三次样条插值函数。

边界问题的提出与类型

单靠一个函数表是不能完全构造出一个三次样条插值函数。我们分析一下其条件个数，条件（2）三次样条插值函数)(x S 是一个分段三次多项式，若用)(x S i 表示它在第i 个子区间],[1i i x x -上的表达式，则)(x S i 形如

],[,)(1332210i i i i i i i x x x x a x a x a a x S -∈+++=

其中有四个待定系数)3,2,1,0(=j a ij ，子区间共有n 个，所以)(x S 共有n 4个待定系数。

由条件（3）)(),(),(x S x S x S '''在],[b a 上连续，即它们在各个子区间上的连接点110,,,-n x x x 上连续即可，共有)1(4-n 个条件，即

???

???

?==-=+''=-''-=+'=-'-=+=-),2,1,0()()1,,2,1)(0()0()1,,2,1)(0()0()

1,,2,1)(0()0(n i y x S n i x S x S n i x S x S n i x S x S i i i i i i i i 共有241)1(3-=++-n n n 个条件，未知量的个数是n 4个。

这样需要加2个条件。

这两个条件通常在插值区间],[b a 的边界点

b a ,处给出，称为边界条件。

边界条件的类型很多，常见有： （1） 给定一阶导数值n

n y x S y x S '='=')(,)(/00； （2）给定一阶导数值n

n y x S y x S ''=''=

'')(,)(//00；（特别地0)(,0)(0=''=''n x S x S 时称为自然边界条件，满足自然边界条件的次样条插值函数称为自然样条插值函数）

（3） 当)(x f y =是周期为a b -的函数时，要求)(x S 及其导数都是以a b -为周期的函数，相应的边界条件为：

)0()0(),0()0(00-''=+''-'=+'n n x S x S x S x S

注：1. 由)(x f 是周期函数

n y y =0，必有)0()0(0-=+n x S x S 。

2. 虽然可利用边界条件及方程组可求出待定系数)3,2,1,0(=j a ij ，从而得到三次样条插值函数)(x S 在各个子区间],[1i i x x -的表达式)(x S i ，但计算量大，不便于应用，我们介绍一种简便的方法（M 方法或三弯方法）： 设

),,2,1,0()(n i M x S i i ==''

因为在子区间],[1i i x x -上)()(x S x S i =是不高于三次的多项式，其二阶导数

)(x S ''必是线性函数，所以

],[,)(11

1

11

i i i i i i i i i i i x x x x x x x M x x x x M x S -----∈--+--='' (1)

设1--=i i i x x h ，则有

],[,)(11

1

i i i

i i i i i i x x x h x x M h x x M x S ---∈-+-='' (2) 连续二次积分得：

i i i i

i i i i i i B x x A h x x M h x x M x S +-+-+-=---)(6)(6)()(13

131 (3)

其中i i B A ,是积分常数，

利用插值条件 i i i i i i y x S y x S ==--)(,)(11，带入（3）中得： )(6

11-----=

i i i

i i i i M M h h y y A

2116

1

i i i i h M y B ---

=， 带入（3），整理得：

i

i i i i i i i i i i i i i

i i i h x x h M y h x

x h M y h x x M h x x M x S 1221133

1

)

6()6(6)(6)()(------+--+-+-= ).,2,1];,[(1n i x x x i i =∈- (4) 从而只要确定n M M M ,,10这1+n 个未知数，即可定出三次样条插值函数)(x S 。

根据)1,,1,0)((-=n i x S i 的连续条件，即：)0()0(-'=-'i i x S x S 等同于: )0()0(1+'=-'+i i i i x S x S , （*） 由（4）得

)(6

2)(2)()(112121-------+-+--='i i i

i i i i i i i i i i M M h h y y h x x M h x x M x S （5）

所以 i i i i i i i i i M h

M h h y y x S 3

6)0(11-+-=

-'-- （6） 将（5）式中的1+→i i ，得)(x S '在子区间],[1+i i x x 上的表达式为)(1x S i +'，

11

11113

6)0(++++++-+-=+'i i i i i i i i i M h M h h y y x S （7）

将（5）（7）代入（*）式中，整理得：

i

i i i i i i i i i i i i h y y h y y M h M h h M h 1

111111636-+++++----=+++，

两编同乘

1

6

++i i h h 得

)(6

211111111i

i i i i i i i i i i i i i i i i h y y h y y h h M h h h M M h h h -+++++-+---+=++++

（8）

令???

??

?

???-+=

-=+=+=-+++++],[],[(61111

1

11

i i i i i i i i i i i i i i i

i x x f x x f h h g h h h h h h μλμ（)1,,2,1-=n i （9） 则所得方程可简写为:

)1,,2,1(211-==+++-n i g M M M i i i i i i λμ

即???????=++=++=++-----1

1121232212121101222n n n n n n g M M M g M M M g M M M λμλμλμ (10) 在（1）边界条件下，即n

n y x S y x S '='=

')(,)(/00， 由于)(x S 在子区间],[10x x 上的导数为

)(6

2)(2)()(011101120112101M M h h y y h x x M h x x M x S ---+-+--='

由 /

00)(y x S ='得 )(62/01

1110y h y y h M M --=

+ （11） 同理由 /)(n

n y x S ='得 )(62/1

11n n

n n n n y h y y h M M +--

=+-- （12） （11）（12）与（10）合在一起得下列的关于n M M M ,,,10 的线性方程组

?????

??

?

???

?????=???????????

???

?

?????????????????----n n n n n n g g g g M M M M 110110

1

111212212 λμλμ （13）

其中???

?

???--=-=-)],[(6

)],[(6//

01010n n n n n y x x f h g y x x f h g （14）

在（2）边界条件下，即n n

n M y x S M y x S =''=''=='')(,)(0//00，在方程中实际上只包含有1-n 未知数

1

21,,-n M M M ，并且可以改写为:

???????

?

??

?

??

???--=???????????

???

?

?????????????????--------//

1122//

1112211

2

2

221

2

2

22n n n n o n n n n n y g g g y g M M M M λμμλμλμλ

(15) 在（3）边界条件下，)0()0(),0()0(00-''=+''-'=+'n n x S x S x S x S 由)0()0(0-''=+''n x S x S 得 n M M =0 （15） 由)0()0(0-'=+'n x S x S

)(6

2)(6211011

10110-----+=---+-n n n n n n n n M M h h y y h M M M h h y y h M

注意到：n n M M y y ==00,，将上式整理得：

)(6

21101111111n

n n n n n n n n h y y h y y h h M h h h M M h h h -----+=++++

记

h

h h n

n +=

1μ=

n λ

]),[],[(6

,1110111n n n

n n n x x f x x f h h g h h h --+=-=+μ

即

n n n n n g M M M =++-211μλ （16）

结合（15）（16）（10）得n M M M ,,21的线性方程组为：

??????

?

?

????????=???????????

???

?

?????????????????----n n n

n n

n n n g g g g M M M M 1211211

1

2211

2

222

μλλμλμμλ

由代数知识知：三对角方程组由于主对角占优，方程有唯一解，具体方法一

般应用三对角的“追赶法”求解。

参数设置

克里金法

克里金法假定采样点之间的距离或方向可以反映可用于说明表面变化的空间相关性。克里金法工具可将数学函数与指定数量的点或指定半径内的所有点进行拟合以确定每个位置的输出值。克里金法是一个多步过程；它包括数据的探索性统计分析、变异函数建模和创建表面，还包括研究方差表面。当您了解数据中存在空间相关距离或方向偏差后，便会认为克里金法是最适合的方法。该方法通常用在土壤科学和地质中。

由于克里金法可对周围的测量值进行加权以得出未测量位置的预测，因此它与反距离权重法类似。这两种插值器的常用公式均由数据的加权总和组成：

(17)

其中：Z(s i)= 第i个位置处的测量值, λi= 第i个位置处的测量值的未知权重, s0= 预测位置, N= 测量值数

在反距离权重法中，权重λi仅取决于预测位置的距离。但是，使用克里金方法时，权重不仅取决于测量点之间的距离、预测位置，还取决于基于测量点的整体空间排列。要在权重中使用空间排列，必须量化空间自相关。因此，在普通克里金法中，权重λi取决于测量点、预测位置的距离和预测位置周围的测量值之间空间关系的拟合模型。以下部分将讨论如何使用常用克里金法公式创建预测表面地图和预测准确性地图。

参数设置

反距离权重法

反距离权重法主要依赖于反距离的幂值，幂参数可基于距输出点的距离来控制已知点对内插值的影响。幂参数是一个正实数，默认值为2。（一般0.5到3的值可获得最合理的结果）。

通过定义更高的幂值，可进一步强调最近点。因此，邻近数据将受到更大影响，表面会变得更加详细（更不平滑）。随着幂数的增大，内插值将逐渐接近最近采样点的值。指定较小的幂值将对距离较远的周围点产生更大的影响，从而导致平面更加平滑。

由于反距离权重公式与任何实际的物理过程都不关联，因此无法确定特定幂值是否过大。作为常规准则，认为值为30的幂是超大幂，因此不建议使用。此外还要牢记一点，如果距离或幂值较大，则可能生成错误结果。

在IDW插值之前,我们可以事先获取一个离散点子集，用于计算插值的权重；

原因1：离散点距离插值点越远，其对插值点的影响力越低，甚至完全没有影响力；

原因2：离散点越少可以加快运算速度；

参数设置