ArcGis地统计插值方法

ArcGis地统计插值方法

地统计（如什么是地统计？专题介绍中所述）是一个方法集，用于估计尚未进行任何采样的位置的值以及评估这些估计的不确定性。这些函数在很多决策过程中至关重要，因为实际操作中不可能对感兴趣区域中的每个位置都进行采样。

但要特别注意，这些方法只是用于构造现实模型（即您感兴趣的现象）的一种手段。至于如何构建模型能够满足您的特定需要并能够为正确合理制定决策提供必要的信息，由您自己（实践者）决定。要构建良好的模型，很大程度上取决于您对现象的理解、采样数据的获取方式和它所表示的内容，以及希望模型提供的信息。构建模型过程的一般步骤在地统计工作流中进行了介绍。

插值方法有很多。有些方法十分灵活，可适用于采样数据的不同方面。有些方法则具有更大的限制性，要求数据满足特定条件。例如，克里金方法十分灵活，但在克里金系列方法中必须满足不同程度的条件才能使输出有效。Geostatistical Analyst 提供了以下插值方法

· 全局多项式

· 局部多项式

· 反距离权重法

· 径向基函数 (RBF) 插值法

· 扩散核

· 核平滑

· 普通克里金法

· 简单克里金法

· 泛克里金法

· 指示克里金法

· 概率克里金法

· 析取克里金法

· 高斯地统计模拟

全局多项式插值法的工作原理

全局多项式插值法可根据输入采样点拟合出一个由数学函数（多项式）定义的平滑表面。全局多项式表面会逐渐变化并捕捉数据中的粗尺度模式。

从概念上讲，全局多项式插值法类似于取出一张纸，然后将其插入凸起点（凸起到一定高度）之间的一个最贴合的位置。下图展示的是从平缓山丘采集而来的一组高程采样点（这张纸是洋红色的）。

但是平整的纸张无法精确贴合带有山谷地形的地表。不过，如果可以将纸张弯曲一下，就会更贴合。为数学公式添加一个项也可以达到类似的效果，即平面的弯曲。平面（纸张无弯曲）是一个一阶多项式（线性）。二阶多项式（二次）允许一次弯曲，三阶多项式（三次）允许两次弯曲，依此类推；在 Geostatistical Analyst 中最多允许 10 次弯曲。下图在概念上展示出一个与山谷拟合的二阶多项式。

纸张几乎无法穿过各实际测量点，这使得全局多项式插值法成为一种不精确的插值器。有些测量点位于纸张上方，而其他点则位于纸张下方。但是，如果将测量点高出纸张的距离相加，并将测量点低于纸张的距离也相加，得到的这两个和值应该相近。以洋红色表示的表面是通过最小二乘回归拟合得到的结果。该生成表面将使凸起点与纸张之间的平方差最小化。

何时使用全局多项式插值法

使用全局多项式插值法获得的是一个可表示感兴趣区域表面渐进趋势的平滑表面。

全局多项式插值法用于下列情况：

· 感兴趣区域（例如，工业区的污染情况）的表面在各位置间出现渐变时，可根据采样点拟合出一个表面。

· 检查和/或消除长期趋势或全局趋势的影响。此类情况下，采用的方法通常为趋势面分析。

在全局多项式插值法中，将利用可描述某种物理过程（例如，污染情况和风向）的低阶多项式创建渐变表面。不过，应注意的是，使用的多项式越复杂，为其赋予物理意义就越困难。此外，计算得出的表面对异常值（极高值和极低值）非常敏感，尤其是在表面的边缘处。

 

局部多项式插值法的工作原理

全局多项式插值法可以根据整个表面拟合多项式，而局部多项式插值法可以对位于指定重叠邻域内的多个多项式进行拟合。通过使用大小和形状、邻域数量和部分配置，可以对搜索邻域进行定义。或者，可以使用探索性趋势面分析滑块同步更改带宽、空间条件数（如果已启用）和搜索邻域值。还可以通过优化按钮  基于交叉验证统计信息为 LPI 进行参数优化。

一阶全局多项式可以根据数据对单平面进行拟合；二阶全局多项式可以对包含一个弯曲的表面进行拟合（表面可以表示山谷）；三阶全局多项可以对包含两个弯曲的表面进行拟合；依此类推。但是，当表面具有多种形状时（如延绵起伏的地表），单个全局多项式将无法很好地拟合。多个多项式平面能够更加准确地体现表面（参见下图）。

另一方面，局部多项式插值法仅使用既定邻域内的点对指定阶数（0 阶、1 阶、2 阶、3 阶，等等）的多项式进行拟合。邻域相互重叠，每次预测所使用值是位于邻域中心的拟合的多项式的值。

在下图中，截取了样本高程数据的横截面（样带）。在第一幅图中，使用三个相邻点（红色点）对一阶多项式进行拟合，使用一条线（红色线）预测由蓝色点标识的位置所对应的未知值。在第二幅图中，通过另一个一阶段多项式对另一个位置（黄色点）进行预测。由于此位置与第一个位置之间的距离很小，因此在预测过程中使用了相同的测量点，但权重稍有不同，从而使得多项式的拟合效果存在些许差异（蓝线）。

随着这一过程的继续，将侧重于后续的预测位置，拟合局部多项式以对值进行预测。以下两幅图显示了为创建最终表面而被预测的另外两个任意点。橙色点是使用经测量的绿色采样点根据拟合的多项式（绿色线）预测而来的，而褐色点是根据浅紫色多项式预测而来的。

 

在以下两幅图中，为预测另外两个位置（蓝绿色点和绿色点）对另外两个多项式（黄色线和灰色线）进行了拟合。

 

将针对各位置重复执行上述过程。您可以看到如何为以下采样点创建表面（紫色表面线）。

与在 IDW 中选择幂 (p) 值相似，也需要在局部多项式中选择最佳参数，以将均方根预测误差 (RMSPE) 降至最小。

准确性的度量

局部多项式插值法提供了两种用于度量准确性的方法，但这两种方法对于 ArcGIS Geostatistical Analyst 中提供的其他确定性插值方法并不适用。

· 预测标准误差用于表示与各位置的预测值相关的不确定性。

· 空间条件数是一种表示预测方程的稳定或不稳定解对特定位置的适应程度的度量方式。如果条件数较大，则很小的矩阵系数变化便会导致解向量（回归系数）的较大变化。由于预测标准误差表面是在假定模型是正确的情况下创建的，因此空间条件数表面可显示数字模型稳定性的变化并提供与预测不确定性相关的其他信息。

局部多项式插值法依赖以下假设：

· 在格网上采样（即，样本的间距相等）。

· 搜索邻域内的数据值呈正态分布。

事实上，大多数数据集都不符合上述假设；在此类情况下，预测值将受到影响，但误差不会像预测标准误差那样大。为帮助确定特定区域内的结果是否可靠，LPI 提供了一个空间条件数表面。下表中显示了一些经验值，这些关键值在“条件数”表面中被渲染成黄色。

多项式的阶	关键空间条件数阈值
1	10
2	100
3	1000
大于 3	大多数情况下都不建议使用

 

低于关键空间条件数阈值的值表示解在哪个位置可靠。接近或等于关键值的值可能存在问题（应用心检验），高于关键限值的值不可靠。

通过估测当线性预测方程的系数发生微小变化时预测值面对此变化的敏感程度，可以生成空间条件数。如果空间条件数较小，则表示解可靠；如果空间条件数较大，则表示解不可靠。如果解的不稳定性出现在特别感兴趣的区域，则应对此问题格外关注。这是因为，输出数据中（包括数据的值、位置和空间排列）的较小变化会导致预测值发生较大变化。这意味着，与输入数据（例如，属性测量中的误差或进行测量的位置处的坐标不精确）相关的任何不确定性都可能对预测值产生重大影响。此外，搜索邻域的变化会修改在预测中使用的数据点的数量（如果平滑搜索邻域发生变化，还会修改权重）并可能影响位置的空间条件数。

为一阶、二阶和三阶多项式创建空间条件数表面。在假定 LPI 模型正确无误的前提下，估测预测标准误差（即，局部加权最小二乘法回归是合适的算法，并且空间条件数值小于上表中的空间条件数阈值）。

通过在 LPI 对话框中将使用空间条件数阈值设置为 True，可以选择将预测中出现较高条件数的区域和预测标准误差图排除在外。条件数仅与数据位置配置相关。换言之，无论是否将同一个数据集中的臭氧值或高程值用作 LPI 的输入，条件数表面都将保持不变。

如果数据规则分布，那么从理论角度来看，最适用于零阶、一阶和二阶多项式的分别为常数核、Epanechnikov 核和四次核。如果数据不规则分布，则应根据验证和交叉验证诊断以及空间条件数值来选择最佳核。

在含障碍的核插值法中提到的核平滑是 LPI 的一个变型。通过使用类似于岭回归的一种方法对这些结果中的局部不稳定性进行校正。取舍问题指的是预测值存在微小偏差，但在大多数实际情况中，该偏差并不足以对根据预测做出的决定造成影响。

表面中存在“孔洞”。

当使用优化按钮  时，通常需要更改带宽值、空间条件数和搜索邻域值。如果空间条件数阈值标记被设置为“假”，则只会更改带宽值和邻域值。如果“邻域”类型被设置为“标准”，则会使用所有数据（最多达 1000 个点）并将最小相邻点数设置为 0。将针对等于 1 个分区的“分区”类型以及平滑邻域进行优化。此优化过程可以导致输出表面出现“孔洞”。出现“孔洞”的位置是指超过空间条件数阈值或搜索邻域过小的那些区域。为填平这些区域，可以对搜索邻域参数和空间条件数阈值进行调整。但必须注意到，之所以创建这些“孔洞”，是因为在计算预测值的过程中这些位置可能存在不稳定性。

何时使用局部多项式插值法

全局多项式插值法适用于在数据集中创建平滑表面及标识长期趋势。然而，在地球科学中，除了长期趋势之外，感兴趣的变量通常还具有短程变化。当数据集显示出短程变化时，局部多项式插值法地图可捕获这种变化。

局部多项式插值法对邻域距离很敏感；较小的搜索邻域可能会在预测表面内创建空区域。因此，可以在生成输出图层之前预览表面。

 

反距离权重插值的工作原理

反距离权重 (IDW) 插值可以明确地验证这样一种假设：彼此距离较近的事物要比彼此距离较远的事物更相似。当为任何未测量的位置预测值时，反距离权重法会采用预测位置周围的测量值。与距离预测位置较远的测量值相比，距离预测位置最近的测量值对预测值的影响更大。反距离权重法假定每个测量点都有一种局部影响，而这种影响会随着距离的增大而减小。由于这种方法为距离预测位置最近的点分配的权重较大，而权重却作为距离的函数而减小，因此称之为反距离权重。在以下示例中对分配给数据点的权重进行了说明：

“权重”窗口包含分配给各数据点的权重的列表，其中各数据点用于生成通过十字光标进行标记的位置处的预测值。

幂函数

如上所述，权重与反距离（数据点与预测位置之间）的 p 次幂成正比。因此，随着距离的增加，权重将迅速降低。权重下降的速度取决于值 p。如果 p = 0，则表示距离没有减小，因为每个权重 λi 均相同，预测值将是搜索邻域内的所有数据值的平均值。随着 p 值的增大，较远数据点的权重将迅速减小。如果 p 值极大，则仅最邻近的数据点会对预测产生影响。

Geostatistical Analyst 使用大于或等于 1 的幂值。当 p = 2 时，此方法称为反距离平方权重插值。将 P = 2 用作默认值。尽管没有理论依据证明该值优于其他值，但应通过预览输出和检验交叉验证统计信息来调查更改 p 值时产生的影响。

可通过将均方根预测误差 (RMSPE) 降至最小值来确定最佳幂值。RMSPE 是在交叉验证过程中计算出的统计数据。RMSPE 用于对预测表面的误差进行量化。Geostatistical Analyst 会对若干个不同的幂值进行评估，从而确定出可以生成最小 RMSPE 的幂值。下图说明了 Geostatistical Analyst 是如何计算出最佳幂值的。针对几个不同的幂值绘制 RMSPE（使用相同的数据集）。根据点拟合一条曲线（局部二次多项式插值法），然后从该曲线上将提供最小 RMSPE 的幂确定为最佳幂。

搜索邻域

由于彼此距离较近的事物比彼此距离较远的事物更加相似，因此，随着位置之间的距离增大，测量值与预测位置的值的关系将变得越来越不密切。为缩短计算时间，可以将几乎不会对预测产生影响的较远的数据点排除在外。因此，通过指定搜索邻域来限制测量值的数量是一种常用方法。邻域的形状限制了要在预测中使用的测量值的搜索距离和搜索位置。其他邻域参数限制了将在该形状中使用的位置。在下图中，在为没有测量值的位置（黄色点）预测值时，将使用五个测量点（相邻点）。

输入数据以及尝试创建的表面会影响邻域的形状。如果数据中不存在方向影响，则需要考虑在各个方向都均等的数据点。为此，将搜索邻域定义为圆。但是，如果数据中存在方向影响，如盛行风，则可能需要对此影响进行调整，方法是：将搜索邻域的形状更改为主轴与风向平行的椭圆。此方向影响的调整是合理的，因为众所周知，对比于与风向垂直但距离预测位置较近的位置，与预测位置风向相逆的位置在远距离位置处会更加相似。

指定完邻域形状之后，便可限制应在该形状中使用哪些数据位置。可以定义要使用的最大位置数和最小位置数，并可以将邻域划分成若干分区。如果将邻域划分成若干分区，则会将最大值和最小值限制应用到各个分区。以下显示了可以使用和显示的几个不同分区。

在数据视图中高亮显示的点显示了在椭圆中心预测某个位置（十字光标的位置）时将使用的位置和权重。搜索邻域被限制在椭圆的内部。在如下示例中，将为南部分区内的一个点（红色）和西部分区内的两个点（红色）分配大于 10% 的权重。在南部分区内，将为一个点（褐色）分配介于 5% 和 10% 之间的权重。将为搜索邻域中的其他点分配较小的权重。

何时使用反距离权重法

使用反距离权重法计算出的表面取决于幂值 (p) 的选择和搜索邻域策略。反距离权重法是一个精确插值器，其中插值表面内的最大值和最小值（参见上图）只能出现在采样点处。输出表面对拓扑和异常值的出现十分敏感。反距离权重法假定正在被建模的现象会受到局部变化的影响，而局部变化可以通过定义合适的搜索邻域进行捕获（已建模）。由于反距离权重法不提供预测标准误差，因此证明使用此模型是否合理可能有点困难。

径向基函数工作原理

RBF 方法是一系列精确插值方法的组合；即表面必须通过每一个测得的采样值。有以下五种基函数：

· 薄板样条函数

· 张力样条函数

· 规则样条函数

· 高次曲面函数

· 反高次曲面函数

在不同的插值表面中，每种基函数都有不同的形状和结果。RBF 方法是样条函数的一个特例。

从概念上讲，RBF 类似于在最小化表面的总曲率时通过测得的样本值拟合橡皮膜。所选基函数确定如何在值之间拟合橡皮膜。下图从概念上演示了 RBF 表面如何通过一系列高程采样值进行拟合。请注意，在横截面上表面穿过数据值。

作为精确插值器，RBF 方法不同于全局和局部多项式插值器，它们都不是精确插值器（不要求表面穿过测量点）。比较 RBF 和 IDW（也是精确插值器）来看，IDW 从不预测大于最大测量值或小于最小测量值的值，如以下样本数据的样带横截面所示。

而 RBF 却可预测大于最大测量值和小于最小测量值的值，如以下横截面所示。

使用交叉验证确定最佳参数，方法与针对 IDW 和局部多项式插值法所述的方法相似。

何时使用径向基函数

RBF 用于根据大量数据点生成平滑表面。这些函数可为平缓变化的表面（如高程）生成很好的结果。

但在表面值在短距离内出现剧烈变化和/或怀疑样本值很可能有测量误差或不确定性时，这些方法不适用。

径向基函数涉及的概念

在 Geostatistical Analyst 中，每个数据位置都会形成 RBF。RBF 是一个随着距某一位置的距离变化而变化的函数。

例如，假设径向基函数只是距每个位置的距离，因此它在每个位置处都形成一个倒圆锥。如果在 y = 5 处绘制 x,z 平面的横截面，将看到每个径向基函数的一部分。现在假设要在 y = 5、x = 7 处预测值。每个径向基函数的值均可取自上图中的每个预测位置，这些值由仅依赖于距每个数据位置的距离的 Φ1、Φ2 和 Φ3 值给出。预测器通过接收加权平均值 w1Φ1 + w2Φ2 + w3Φ3 + … 而形成

现在的问题是如何确定权重？到目前为止，根本没有使用过数据值。在将预测移到含有测量值的位置时，通过要求精确预测测量值找到了权重 w1、w2 和 w3 等。这形成了含有 N 个未知数的 N 个方程，它们只有一种求解方法。因此，表面穿过数据值以做出精确预测。

此例中的径向基函数是高次曲面 RBF 的一个特例。Geostatistical Analyst 还允许使用其他 RBF，例如规则样条函数、薄板样条函数、张力样条函数和反高次曲面样条函数。有时，这些函数的差别不大，但可出于一定原因而选择某一个函数，也可尝试多个函数并使用交叉验证选择一个。每个 RBF 都有一个控制表面“平滑度”的参数。

对于除反高次曲面外的所有方法，参数值越高，地图越平滑；对于反高次曲面则正相反。

“含障碍的扩散插值”的工作原理

扩散插值指热方程的基本解，该方程描述热或粒子如何在同类介质中随着时间扩散。使用此方法得到的预测结果缓慢地在障碍周围流动。

在没有障碍的情况下，通过扩散插值得到的预测结果与通过使用高斯核的核插值法得到的预测结果大致相同。

扩散插值可以使用通过成本面定义的复制的距离度量，它是一个常用栅格函数，用于计算从一个栅格像元到另一个栅格像元的行程成本。

扩散插值为自动选择的格网（像元）生成预测结果，而 Geostatistical Analyst 中的所有其他模型使用大小可变的三角形。

如果未指定成本面，相邻像元中心之间的距离为欧氏距离。如果已经指定成本面，则使用以下公式之一定义成本面上像元之间的距离：

· 附加障碍

（相邻像元的平均成本值）x（像元中心间的距离）

· 累积障碍

（相邻单元的成本值之差）+（单元中心间的距离）

· 流动障碍

指示符（至相邻像元的成本值 > 自相邻像元的成本值）*（至相邻像元的成本值 - 自相邻像元的成本值）+（像元中心间的距离）

其中，Indicator(True) = 1，Indicator(False) = 0。

流动障碍可用于使用数据变化的主方向插入数据。

可将下图左侧没有成本面的“含障碍的扩散插值”结果与下图右侧使用“含障碍的核插值法”创建的预测图相比较。

在扩散插值的情况下，核的形状根据扩散方程在障碍附近变化；而在核插值的情况下，点之间的距离根据点之间的最短距离变化。

含障碍的核插值法的工作原理

ArcGIS 10 

核插值是一阶局部多项式插值法的一个变型，此方法使用一种类似于在用于估算回归系数的岭回归中使用的方法来防止在计算过程中出现不稳定性。当评估值仅存在较小偏差且比无偏差评估值更加精确时，可以将其作为首选评估值。例如，可以在 Hoerl and Kennard (1970) 中了解到岭回归方面的详细信息。

两个模型之间的另一个区别在于，核插值模型使用两点间的最短距离，因此，位于指定的不透明（绝对）障碍的侧面上的点会通过一系列直线相连。

核插值法使用以下径向对称核：指数、高斯、四次式、Epanechnikov、5 阶多项式和常数。使用交叉验证诊断根据数据对核的带宽进行估算。

当使用一阶多项式时，Epanechnikov 核通常会生成比较理想的结果。不过，根据数据，交叉验证和验证诊断可能会建议使用另一种核，Fan and Gijbels (1996)。

下面对具有绝对障碍的含障碍的核插值法预测（左侧）和不具有绝对障碍的含障碍的核插值法预测（右侧）进行比较：

在水文和气象应用中，可优先采用基于两点间最短距离的模型。

参考材料和更多阅读材料

Fan, J. and Gijbels, I. (1996).Local Polynomial Modelling and Its Applications, Chapman & Hall.London.

Hoerl, A.E. and Kennard, R.W. (1970), Ridge regression:biased estimation for nonorthogonal problems, Technometrics, 12, 55-67.

Yan, Xin.(2009) Linear regression analysis :theory and computing.Published by World Scientific Publishing Co. Pte. Ltd. 5 Toh Tuck Link, Singapore 596224.

了解普通克里金法

普通克里金法假设模型为

Z(s) = ? + ε(s),

其中，? 是一个未知常量。对于普通克里金法，我们所关心的主要问题之一就是对常量平均值的假设是否合理。有时有很充分的科学依据来拒绝该假设。不过，作为一种简单的预测方法，它具有显著的灵活性。下图所举的是处于某一空间维度中的示例：

从图上看，数据好像是从山谷或山体的线横断面中采集的高程值。而且，好像数据在左侧变化更显著，而在右侧则变得更平滑。事实上，该数据是在平均值 ? 为常量的情况下基于普通克里金法模型模拟得到的。虚线给出的是平均值，该平均值是是真值但是是未知的。因此，普通克里金法可用于似乎带有某种趋势的数据。单凭数据无法确定已观测到的模式是否是自相关（? 为常量的情况下，在误差 ε(s) 之间）或趋势（?(s) 随 s 变化）所造成的。

普通克里金法可以使用半变异函数或协方差（用于表达空间自相关的数学形式），使用变换和移除趋势，还允许测量误差。

了解简单克里金法

ArcGIS 10 

简单克里金法假设模型为

Z(s) = ? + ε(s),

其中，? 是已知常量。例如，下图中使用了和普通克里金法和泛克里金法概念所使用的相同的数据，观测数据以实心圆的形式给出。

已知常量（实线）为 ?。这点可与普通克里金法进行比较。对于简单克里金法，因为假设确切已知 ?，那么也确切已知数据位置上的 ε(s)。对于普通克里金法，如果估算了 ?，因此也会估算 ε(s)。如果已知 ε(s)，可以比估算 ε(s) 时更好地估算自相关。通常，假设确切已知平均值 ? 是不现实的。但是，有时假设基于物理的模型能够给出已知趋势却是有意义的。由此可以使用模型和观测值的差值（称为残差），并且假设残差中的趋势已知为零，可以在残差上使用简单克里金法。

简单克里金法可以使用半变异函数或协方差（用于表达自相关的数学形式）和变换，并且允许测量误差。

了解泛克里金法

ArcGIS 10 

泛克里金法假设模型为

Z(s) = ?(s) + ε(s),

其中 ?(s) 为某些确定性函数。例如，下图中使用了和普通克里金法概念所使用的相同的数据，观测数据以实心圆的形式给出。

二阶多项式为趋势（长虚线）为 ?(s)。如果从原始数据减去二阶多项式将得到误差 ε(s)，并假设该误差是随机的。所有 ε(s) 的平均值为 0。从概念上讲，自相关现在通过随机误差 ε(s) 来建模。当然，也可以拟合线性趋势、三阶多项式或任何数目的其他函数。上图看上去就像基础统计课程的多项式回归。实际上其为泛克里金法。使用作为解释变量的空间坐标进行回归。但是，没有假设误差 ε(s) 是独立的，而是将它们建模为自相关的。对泛克里金法的建议和普通克里金法的建议相同：只根据数据是无法确定的，还需要采用适当的分解。

泛克里金法可以使用半变异函数或协方差（用于表达自相关的数学形式）和变换，并且允许测量误差。

了解指示克里金法

ArcGIS 10 

指示克里金法假设模型为

I(s) = ? + ε(s),

其中，? 是一个未知常量，I(s) 是一个二进制变量。二进制数据的创建可利用连续数据的阈值实现，或者观测数据可以为 0 或 1。例如，假设存在一个由某点是否为森林栖息地或非森林栖息地的相关信息组成的样本，则其中二进制变量用来指示这两种类别。使用二进制变量时，指示克里金法的处理过程与普通克里金法相同。

在下图中，已使用了解阈值中介绍的阈值将数据转换为二进制值。

用空心方块给出了观测的二进制数据。虚线表示所有指示变量的未知平均值，即 ?。可以将这一点与普通克里金法进行比较。在使用普通克里金法时，假设 ε(s) 是自相关的。请注意，因为指示变量为 0 或 1，所以插值结果将位于 0 和 1 之间，而且基于指示克里金法的预测可解释为变量是 1 的概率或属于 1 所指示的类别的概率。如果创建指示变量时使用了阈值，则生成的插值地图会显示超出（或低于）阈值的概率。

通过选择多个阈值可以为同一数据集创建多个指示变量。在本例中，一个阈值创建主要指示变量，而另一个指示变量则用作协同克里金法中的次要变量。

指示克里金法可使用半变异函数或协方差，它们都是用于表达自相关的数学形式。

了解概率克里金法

ArcGIS 10 

概率克里金法假设模型为

I(s) = I(Z(s) > ct) = ?1 + ε1(s)  Z(s) = ?2 + ε2(s),

其中 ?1 和 ?2 为未知常量，I(s) 是通过使用阈值指示 I(Z(s) > ct) 创建的二进制变量。请注意，现在有两种类型的随机误差：ε1(s) 和 ε2(s)，因此它们各自存在自相关，并且它们之间存在互相关。概率克里金法要实现指示克里金法相同的功能很吃力，而使用协同克里金法进行尝试则可更好地实现。

例如，在下图中普通克里金法、泛克里金法、简单克里金法和指示克里金法概念使用相同的数据，请注意标注为 Z(u=9) 的基准的指示变量为 I(u) = 0，标注为 Z(s=10) 的基准的指示变量为 I(s) = 1。

如果要预测它们中间的位于 x 坐标 9.5 处的值，单独使用指示克里金法将给出接近 0.5 的预测值。但是，可以看出 Z(s) 刚好高于阈值，而 Z(u) 却远低于阈值。因此，有理由相信位置 9.5 处的指示预测值应该小于 0.5。概率克里金法尝试利用原始数据中除二进制变量之外的其他信息。但是，这也存在一些代价。必须要进行更多的估算，包括估算每个变量的自相关和互相关。然而，每次估算未知的自相关参数时，都会引入更多的不确定性，因此概率克里金法可能不值得付出额外努力。

概率克里金法可以使用半变异函数或协方差（用于表达自相关的数学形式）、交叉协方差（用于表达互相关的数学形式）和变换，但是不允许测量误差。

了解析取克里金法

ArcGIS 10 

析取克里金法假设的模型为

f(Z(s)) = ?1 + ε(s),

其中，?1 是一个未知常量，f(Z(s)) 是 Z(s) 的一个任意函数。请注意，您可以写成 f(Z(s)) = I(Z(s) > ct)，这样，指示克里金法就成为析取克里金法的一种特殊情况。在 Geostatistical Analyst 中使用析取克里金法时，您既可预测值本身，也可预测指示器。

在 Geostatistical Analyst 中，提供的 g(Z(s0)) 函数其实就是 Z(s0) 本身和 I(Z(s0) > ct)。一般来说，相比普通克里金法，析取克里金法可以做更多事情。尽管回报更丰厚，但成本也更高。析取克里金法要求接受二元正态分布假设和 fi(Z(si)) 函数的近似值；但是很难对假设进行验证，而且从数学和计算角度来看，解决方案都很复杂。

析取克里金法可使用半变异函数或协方差（用于表达自相关的数学公式）以及变换，但是它不允许出现测量误差。

地统计模拟的重要概念

ArcGIS 10 

模拟概念

模拟在广义上是指使用模型复制现实的过程。在地统计中，模拟是随机函数（表面）的实现，其与生成该模拟的样本数据拥有相同的地统计要素（使用均值、方差和半变异函数来度量）。更具体地说，高斯地统计模拟 (GGS) 适用于连续数据，并假设数据或数据的变换具有正态（高斯）分布。GGS 所依托的主要假设是数据是静态的 - 均值、方差和空间结构（半变异函数）在数据空间域上不发生改变。GGS 的另一个主要假设是建模的随机函数为多元高斯随机函数。

同克里金法相比，GGS 具有优势。由于克里金法是基于数据的局部平均值的，因此，其可生成平滑的输出。另一方面，GGS 生成的局部变异性的制图表达比较好，因为 GGS 将克里金法中丢失的局部变异性重新添加到了其生成的表面中。对于由 GGS 实现添加到特定位置的预测值中的变异性，其平均值为零，这样，很多 GGS 实现的平均值会趋向于克里金预测。下图对此概念进行了说明。各种实现以一组堆叠输出图层的形式表示出来，并且特定坐标位置的值服从高斯分布，其平均值等于该位置的克里金估计值，而扩散程度则由该位置上的克里金法方差给出。

提取值到表工具可以用来为上图中的图形生成数据，在对 GGS 生成的输出进行后处理时该工具也很有用。

对 GGS 的使用在地统计实际操作中日益呈现出一种趋势，它不是追求获得每个未采样位置的最佳无偏预测结果（正如克里金法所体现的），而是强调对决策分析和风险分析的不确定性的特证描述，这样更适合于呈现数据中的全局趋势 (Deutsch and Journel 1998, Goovaerts 1997)。模拟还会克服克里金估计值中的条件偏差带来的问题（高值区域预测值通常偏低，而低值区域预测值通常偏高）。

对于所研究属性的空间分布，地统计模拟可为其生成多个具有同等可能性的制图表达。可基于这些制图表达来测量未采样位置的不确定性，这些未采样位置在空间上被一起选取，而不是逐个被选取（如同通过克里金法方差进行测量一样）。此外，克里金法方差通常独立于数据值，且通常不能用作估计精度的测量值。另一方面，可以通过使用多个模拟实现（该实现用呈正态分布的输入数据通过简单克里金模型进行构建，即，数据呈正态分布或已使用常态得分变换或其他类型的变换对数据进行了变换）为未采样位置的估计值构建分布来测量估计精度。对于使用估计数据值的风险评估和决策分析而言，这些不确定性的分布很关键。

GGS 假设数据呈正态分布，但在实际中，很少会出现这种情况。对数据执行常态得分变换，使得数据符合标准正态分布（均值 = 0，方差 = 1）。然后对此正态分布数据进行模拟，并对结果做反向变换，以便以原始单位获得模拟输出。对正态分布数据使用简单克里金法时，该克里金法所提供的克里金估计值和方差可完全定义研究区域中每个位置的条件分布。这样，您可以在只知道每个位置的这两个参数的情况下绘制随机函数（未知采样表面）的模拟实现，这也是 GGS 基于简单克里金模型和正态分布数据的原因。

“高斯地统计模拟”工具支持两种类型的模拟：

· 条件模拟遵循数据值（除非克里金模型中包含测量误差）。由于模拟会在格网像元中心生成值，因此，如果此值与采样点的位置不完全对应，则采样位置的测量值与模拟值可能会不同。条件模拟也将以平均方式（即，在很多实现上平均）复制数据的均值、方差和半变异函数。模拟表面看起来很像克里金预测地图，但其将显示更多的空间变异性。

· 非条件模拟不遵循数据值，但会以平均方式复制数据的均值、方差和半变异函数。模拟表面所显示的空间结构类似于克里金地图，但输入数据中存在高值或低值的地方不一定会出现高值和低值。

模拟示例

示例 I

在世界上的许多城市和地区，空气质量都是令人关注的重要健康指标之一。在美国，众所周知，洛杉矶的空气质量不是很好，分布密集的监控网络每半天就对臭氧、微粒物质和其他污染物等数据进行一次收集。基于此空气质量数据，可获得每种污染物的浓度以及污染物每年超过州空气质量标准和联邦空气质量标准的天数 (http://www.arb.ca.gov/aqd/aqdpage.htm)。由于这两个测量值均支持对在某个特定区域内生活进行感染风险的局部评估，因此，每年超过临界阈值的天数可用来建立显示超过阈值概率的内插地图。

在本示例中，对 2005 年加利福尼亚州每个监测站臭氧超过阈值的天数做了调查，并通过拟合该数据创建了一个半变异函数。并使用条件模拟生成了多个实现。每个实现都是一个地图，用于表示 2005 年污染物超过阈值的天数。然后对这些实现进行后处理，以估计污染物每年超过州阈值的天数多于 10 天、20 天、30 天、40 天、50 天、60 天和 70 天的概率（所有监测站记录的超过阈值的最大天数为 80 天）。下面的动画显示了生成的南海岸空气盆地地区（其中包括洛杉矶和内陆城市）的臭氧地图。海岸附近的空气质量明显好于内陆地区，主要是因为在这一地区，风向主要是由西向东吹。

这类地图可用于确定污染减轻策略的优先级，通过解答诸如“我可以忍受多少污染？”、“生活在某一特定区域我需要忍受多少污染？”等问题，来研究健康与环境质量之间的关系并帮助人们确定适宜居住的地点。

示例 II

在很多应用中，都使用与空间相关的变量作为模型的输入（例如，石油工程中的流动模拟）。在此类情况中，模型结果的不确定性是通过以下过程生成大量模拟来进行评估的：

· 1. 为变量模拟大量具有同等可能性的实现。

· 2. 使用模拟变量作为输入来运行模型（通常称为传输函数）。

· 3. 汇总模型运行以评估模型输出的变异性。

 

上述过程的一个实际示例是：为在新墨西哥州东南部成立一个废品隔离试验工场 (WIPP) 作为超铀废物的存储设施而进行的研究 (http://www.wipp.energy.gov/)。

科学家曾对位于地表以下 2000 多英尺的盐矿床进行了评估，以便将其用作废料的潜在存储设施。然而，矿床刚好位于蓄水层之上，因此，担心地下水可能会传输站点泄露的废弃物。要证明 WIPP 是安全的，科学家不得不使美国环境保护局确信：蓄水层中地下水的流速非常缓慢，完全不可能对 周围环境造成污染。

导水系数值决定了蓄水层中的水流流速，并针对拟建的 WIPP 站点附近蓄水层获得了多个此类值。使用以数字方式求解的水文方程为地下水流建模，该方程需要导水系数值，该值在常规格网上进行预测。如果使用了导水系数的克里金估计值，则导水系数值将基于邻近导水系数值的（加权）平均值，而已建模的地下水的流动时间将只会基于这些平均值。由于克里金法将生成平滑地图，所以插值表面会缺少导水系数值极高或极低的区域。要正确地对风险进行分析，科学家必须考虑可能出现的最坏情况，因此需要生成流动时间值的整个概率分布。通过此分布，科学家将能够使用地下水流动时间分布的较低尾值（对应极高流速），而不是平均流动时间，来评估 WIPP 的适宜性。曾使用条件模拟来生成流动时间值的概率分布。

废品通过地下水进行传输的概率只是评估 WIPP 适宜性时考虑的众多危及人类健康情形中的一种。复杂风险分析在评估 WIPP 是否适宜进行核废料处理以及使公众和政府监管部门确信其适宜性方面起了很大作用。在长达 20 多年的时间里，在进行了大量的科学研究、公众意见收集以及进行了大量监管工作之后，WIPP 最终于 1999 年 3 月 26 日开始运作。

应该生成多少实现？

模拟研究的结果不应取决于所生成实现的数量。确定生成多少实现的其中一种方法是：在一小部分数据属性域中对比不同实现数的统计数据（使用子集以节省时间）。随着实现数量的增加，统计数据将趋向于一个固定值。下面的示例中检查的统计数据是第一个分位数和第三个分位数，它们是为美国斯威康星州的一小部分（子集）模拟高程表面（在海平面以上，以英尺为单位）而计算的值。

上方的图显示的是前 100 个实现的高程波动。下方的图显示的是 1000 个实现的结果。

 

在本例中，值在大约 20 个模拟之后稳定下来。在很多情况下，至少需要运行 100 个实现才能确定超出阈值的均值和概率所需的足够信息。如果使用数量更多的实现，则可为汇总统计数据和模型输出变量提供更高程度的确定性，但所需计算时间也更长。

要了解有关如何在 ArcGIS 中实现高斯地统计模拟的详细信息，请参阅帮助部分高斯地统计模拟的工作原理。

参考书目

Deutsch, C.V., and A. G. Journel. 1998. GSLIB Geostatistical Software Library and User's Guide.2nd Ed. Oxford University Press, New York, pages 119–122.

Goovaerts, P. 1997. Geostatistics for Natural Resource Evaluation.Oxford University Press, New York, pages 369–376.

转自http://chenxehua.blog.163.com/blog/static/1074066702013921103837586/