9.5. Grid的生成¶
9.5.1. 网格化插值计算¶
将离散 DEM 数据经插值计算转换为格网 DEM 数据的过程称为 DEM 数据网格化,即生成Grid。 也就是说,需要用离散的观测点值去估算未知的格网点的值(图9-14)。 插值的算法有许多种,但不论用哪种方法,如果希望通过插值加密数据点,使一个粗糙的初始地形模型变成一个新的精确模型, 那肯定是不现实的。 相反,原始模型经过插值计算得到的新模型总是可能带来一定的变形,因为原始数据点是实测的,因而是精确的, 而插值点可能会有一定的误差。 此外,应注意到地形起伏变化的不均衡性,有的地区平坦和缓,有的地区起伏骤变,用统一网格较难适应这种变化的差异。 但是设计和选用适宜的插值算法仍然是必要的。 使原始数据中包括的地理特征能够无明显损失地传递到内插计算的数字高程模型中去。

图 9.13 格网化过程¶
数字化 格网化
自五十年代末 DEM 方法提出后,已从基本理论到内插方法进行了大量研究,相继出现了多种插值算法,积累了丰富的经验。 概括他们的较为统一的认识,将有助于设计和选择插值算法,建立起精度较高的 DEM 数据,这就是:
(1)对于地形模型,因无一定数学规律可循,其精度主要决定于原始数据的获取(密度和分布)。
(2)一般情况下, DEM 精度与原始数据密度基本上呈线性关系。
(3)在原始数据一定密度的条件下,尤其密度较高条件下,在同一样点网上展铺不同类型的数学面, 产生的 DEM 精度无明显差别。
(4)由于地形的千变万化──平原、丘陵、山地、黄土、喀斯特等等,各种算法只有一定的适应性。
(5)由于地形在地性线(山脊、山谷、断崖等)上发生转折,在分块拟合曲面时,跨越地性线的数据点不宜参加曲面方程的拟合。
9.5.2. 网格尺寸的确定¶
这是 DEM 数据网格化的一个重要问题,它首先关系到派生数据的密度,并直接影响地形模型的精度。 网格过细不仅不会提高 DEM 精度,反而产生冗余的“游离”数据,即相邻格网点值差别微小且不包含有效的地形特征信息。 网格过大则不仅丢失了地形特征信息且会造成地形扭曲。 并且,若网格过大,派生的网格点数明显低于样点采集数,那意味着样点的采集工作就是失败的,采集了许多实际无用的样点。
因此,一般情况下,样点的密度基本上决定了网格点密度。 网格点数宜大于或接近样点数。 在采集优选点情况下,可考虑:
n〈 N〈 2n
(公式9-6)
式中N为网格点数,n为采样点数。
此外,网格尺寸的选定还应分析地形形态特征,如图9-15所示,原等高线明显地反映了微地形特征,若是采用L边长的网格, 则这些特征将消失,再现的等高线必定是一条平直的曲线。 为此,可选取原图上一定数量的有代表性的等高线,分析其曲率变化情况,而后决定网格的适宜尺寸。 图9-16a表示一条n点组成的等高线,图9-16b显示其中三点构成一段曲线。 当Di≥K时(K根据 DEM 精度要求确定),计算Si值,i=1,2,···,m,m≤n。 计算:
- ::
Smin=min{S1,S2,···,Sm }
(公式9-7)
Smin值可以作为确定网格尺寸的参考值。
但为了顾及整个制图区域,宜在Si中取若干个小值S小值,并且以其平均值小值为参考值,可能比Smin更适宜。
此时可考虑;
l≤\ |image2|\ 小值
(公式9-8)
式(9-6)、(9-7)、(9-8)仅可作为网格尺寸的参考值使用。 网格尺寸主要决定于制图目标、 DEM 数据使用方法和精度要求等因素,并应通过试验,选择适宜的尺寸。

图 9.14 网格尺寸的影响¶

图 9.15 图9-16¶
9.5.3. 空间插值方法¶
1. 移动平均法
移动平均法认为任一点上场的趋势分量可以从该点一定邻域内其它各点的值及其分布特点平均求得,参加平均的邻域称做窗口。 窗口的形状可以是方形或圆形。 圆形比较合理,但方形更方便计算机取数。 求平均时可以用算术平均值、众数或其它加权平均数。 选用大小不同的窗口,可以实现数据的分解,大窗口使区域趋势成分比重增大,小窗口则可突出一些局部异常。 逐格移动窗口逐点逐行计算直到覆盖全区,就得到了网格化的数据点图。
移动平均在求趋势分量时只涉及一定范围,因此在接图的边界上,只要适当的扩边,把相邻图幅影响范围内的数据接上一起处理, 即可以方便地实现拼接。 这类方法已成为各国区域化探数据处理的标准方法。 我国区域化探中常用的是8×8km的窗口。
移动平均保持了对一般趋势的反映,而且很容易填补一些小的数据空缺,使图面完整,但移动平均有一定的平滑效应和边缘效应。
当原始取样点分布较稀且不规则时,可以采用定点数而不定范围的取数方法,即搜索邻近的点直到预定的数目为止。 搜索方法可以是四方搜索或八方搜索等。 此时由于距离可能相差较大,因此常同时采用距离倒数或距离平方倒数加权的办法,以便压低远处的点的影响。
2. 距离平方倒数加权法
距离平方倒数加权法主要原理是某点(或某待估网格)的估计值与周围已知点值的距离平方倒数成一定关系, 以空间位置的加权平均来计算。
(公式9-9)
设平面上分布一系列离散点,已知其坐标和高程为Xi、Yi、Zi、(i=1,2,…,n),
(X,Y)为任一格网点,根据周围离散点的高程,通过距离加权插值求
点高程,
周围点与
点因分布位置的差异,对
影响不同,我们把这种影响称为权函数
。
权函数主要与距离有关,有时也与方向有关。
若是在
点周围四个方向上均匀取点,那么可不考虑方向因素,这时
(公式9-10)
实践证明, 是较优的选择,
为离散点至
点的距离:
(公式9-11)
这样,在每一格网点周围搜索若干靠近离散点,用以逐一内插格网点高程,建立一个网格 DEM 。 这种算法的前提是离散点均匀分布,且每一离散点具有同等的意义。
3. 趋势面拟合技术
地球表面起伏变化、千姿百态。
多年来为了准确合理地表现这一形态,人们曾就地形曲面数学模型以及相应的数值逼近方法作过不少努力。
多项式回归分析是描述大范围空间渐变特征的最简单方法。
多项式回归的基本思想是用多项式表示线(数据是一维时)或面(数据是二维时)按最小二乘法原理对数据点进行的拟合,
一维的拟合称为线拟合技术,但地理信息系统研究的对象在空间上和时间上都有复杂的分布特征,
在空间上的分布常是不规则的曲面,数据往往是二维的,而且以更为复杂的方式变化,一维拟合技术不能反映区域性趋势变化,
因此必须利用趋势面分析技术。
其基本思想是:用函数所代表的面来逼近(或拟合)现象特征的趋势变化。
拟合时假定数据点的空间坐标、
为独立变量,
而表征特征值的
坐标为因变量。
当数据处在一维空间时,一元回归函数为:
(一元线性回归) (公式9-12)
(一元非线性回归)
(公式9-13)
当数据处在二维空间时,二元回归函数为:
(一次多项式回归) (公式9-14)
(二次多项式回归) (公式9-15)
上述公式中,、
、
,
、
、
、
、
、
为多项式系数。
当n个采样点上观测值
和估计值
的离差平方和为最小时,即
(公式9-16)
则认为回归方程与被拟合的线或面达到了最佳配准(图9-17,图9-18)。 且可计算出多项式系数。

图 9.16 一元线性回归分析¶

图 9.17 一元非线性回归分析¶
回归函数的次数并非越高越好,在实际工作中,一般只用到二次,超过三次的复项多项式往往会导致解的奇异, 高次趋势面不仅计算复杂(如六次趋势面方程系数达28个),而且次数高的多项式在观测点逼近方面效果虽好, 但在内插、外推的效果上则常常降低分离趋势的作用,使整体趋势分离,降低趋势规律的反映。 趋势面是一种平滑函数,很难正好通过原始数据点,这就是说在多重回归中的残差属正态分布的独立误差, 而且趋势面拟合产生的偏差几乎都具有一定程度的空间非相关性。
整体趋势面拟合除应用整体空间的独立点内插外,另一个最有成效的应用之一是揭示区域中不同于总趋势的最大偏离部分。 因此,在利用某种局部内插方法以前,可以利用整体趋势面拟合技术从数据中去掉一些宏观特征(例如最小二乘配置法)。
4.样条函数
最小二乘曲面拟合假设了所有样品值被观测到的概率相等,而没有考虑样品间的相对位置。
当观测点数比较大时,需要用高阶多项式去拟合,这不但使计算复杂化,并且高阶多项式还可能在观测点之间产生振荡。
因此,多采用分块拟合的办法,用低阶多项式进行局部拟合。
样条函数拟合即是常用的方法,它将数据平面分成若干单元,在每一单元上用低阶多项式,
通常为三次多项式(三次样条函数)构造一个局剖曲面,对单元内的数据点进行最佳拟合,
并使由局部曲面组成的整个表面连续。
Akima在1978年提出了用双五次多项式和连续的一阶偏导数进行光滑曲面拟合和内插的方法,
称为Akima样条插值法。
将平面分割为三角形格网,各三角形以三个数据点在(
)平面上的投影点为顶点,
三角形内的某点(
)的值用下列公式内插得出:
=
(公式9-17)
其中,i,j=1,2,…,5,qij表示多项式系数矩阵元素。
根据三个顶点的场值、一阶偏导数和二阶偏导数值,可得到18个不相关的条件, 三角形三条边两侧的一阶偏导数相等给出另外三个边界条件,这样可求出方程的21个系数。
5. 克立金法
克立金法最初是由南非金矿地质学家克立金(D.G.Krige)根据南非金矿的具体情况提出的计算矿产储量的方法; 按照样品与待估块段的相对空间位置和相关程度来计算块段品位及储量,并使估计误差为最小。 然后,法国学者马特隆(G.Matheron)对克立金法进行了详细的研究,使之公式化和合理化。
克立金法基本原理是根据相邻变量的值(如若干样品元素含量值), 利用变差函数所揭示的区域化变量的内在联系来估计空间变量数值的方法。
地质变量是区域化变量,具有空间结构性,即在空间点x和x+h处的变量值z(x)和z(x+h)具有自相关性。 这种相关依赖于两点间的向量h和矿化特征。 区域化变量的空间特征由变差函数来描述。
变差函数为区域变量z(x)的增量平方的数学期望,即区域化变量增量的方差。 变差函数即是距离h的函数,又是方向a的函数,通式可写成:
=
=E
![]()
(公式9-18)
地质统计学中,多用半变差函数表示,即将上式左边的因子2移至右端作分母。
变差函数一般以变差曲线表示(图9-19)。
由图可见,随着 的增大,r(h)趋于稳定值,这时的h称为变程,记为a,
它表示了变量从空间相关状态到不相关状态的转折点,而r
(a)称为基台值(c+c0)。
变程揭示了变量空间自相关性的影响范围,基台值反映了变量变异的强弱。
当h趋于零时,r(h)的极限值即曲线在纵坐标的截距为块金常数(c0),
它反映了变量随机性的大小(图9-19)。

图 9.18 变程(a),基台值(c+ c0),块金常数(c0)在变差曲线的位置¶
对不规则分布的数据,
按一定的方向
和距离 将数据点组合成角组和距离组(图9-20)。
和
分别为角度和距离的允许误差限。
在
范围内的数据都被认为是沿着
方向。
在每个方向上,以一维规则分布的计算公式为基础,用求
的算术平均值方法计算实验变差函数
,作出实验曲线,并用理论曲线拟合,求出a、c、c0等特征参数值。
(公式9-19)
其中, 为信息值,
为对应于某个
的数据对的个数,
为滞后距。
例如,如图9-21所示,沿x方向分布9个数据点,则对应于不同的 可求出不同的
值。
如对于
,
,
2.4。
时,
,
2.7。
图9-21一维规则数据实验半变差函数计算数据构形
图9-20 二维不规则数据
在实际工作中, 当怀疑有异向性的特定方向β存在时,可以确定一个角度允许误差限dβ,以准确地估计这一方向上的变差函数。
克立金估值要求满足无偏和最优两个条件。
任一估计域的主量 可以通过该域影响范围内几个有效信息值
的线性组合得到,即
|image95|
(公式9-20)
是与
有关的加权系数,
用来表示各个信息值
对估值
的贡献。
对于任一给定的块段和数据信息
存在一组加权系数
,
要求满足:
(1) 无偏估计
所有估计域的实际值与预测值之间的偏差ε平均为零,即估计误差的期望值为零
(公式9-21)
(2)最优估计
估计域的估计值与真实值之间的单个偏差应尽可能地小,用估计方差表示
(公式9-22)
可以有效地用来度量估计值的精度。
根据无偏条件和估计方差最小,得到下列方程组
(公式9-23)
也可用协方差函数
表示
(公式9-24)
要使 在无偏条件下最小,则求条件极值。
采用标准拉格朗日乘数法求解得到克立金方程组
(公式9-25)
克立金方差
(公式9-26)
用矩阵表示为
(公式9-27)
其中
,
,
(公式9-28)
为对称阵,
。
当数据点较多时,用克立金法估值,需解高维克立金方程组。 计算中采用超级块段思想减少克立金方程组的数目和维数。
9.5.4. 几种典型数据网格化插值方法选择¶
遥感数据是按影像方式记录的栅格数据,内插放大或重采样时,常用矩形网格内插法, 如最邻近点法、双线性插值法或立方卷积法。
地球物理数据,特别是位场数据,是典型的空间连续型数据。 一般多用样条函数插值方法,使生成的曲面具有连续的二阶导数和最小的平方曲率。 三次样条插值比较适合于高频成分较多的场,对于台阶异常,Akimn样条插值可能显得更为合理, 也可以用最小二乘曲面拟合法和距离反比加权法。
一些资料的测线间距大于探测目标的埋深,形成欠采样资料。 当测线与场源地质体不垂直时,常规插值常常形成虚假孤立异常,这时可用方向增强插值的方法弥补采样的缺陷。 在对水系沉积物或分散流数据加权插值时,可考虑沿水系的方向给予较大的权。 对测线与目标地质体走向不垂直时的数据,可选取长矩形插值窗口, 矩形的长边平行于地质体走向并同时加大走向方向各点的权重。
化探异常数据具有较强的随机性和采样点稀疏不规则的特点, 因此网格化估值方法常用滑动平均法、距离平方倒数法和克立金法。