19.4. 使用 SAGA 地理信息系统进行插值#
空间插值包括从仅在少数位置可用的值样本重建域上的空间数据的值。 存在许多复杂性不同的方法。 在这个例子中，我们将使用一种简单、方便且常用的方法。 不切实际的是，当我们执行空间插值时，我们将点的形状文件转换为栅格。 我们在这里不会讨论各种方法的优点和缺点。 事实上，内插本身可能是一个地质统计学主题。
我们提请读者注意，从水文地质的角度来看，这些插值可能不相关。 这项练习不应被视为水文地质学中的练习，而应被视为空间插值的开始。
以下以 塞纳河流域含水层水质措施的插值作为实例进行说明。
19.4.1. 数据#
在本练习中，我们对塞纳河盆地地下水的质量进行了一些测量，如 图 19.36 。 塞纳河位于法国北部（北纬48°，东段2°）。
图 19.36 法国北部塞纳河盆地的位置#
这些数据存储在名为 "groundwater_oct_2006.shp" 的形状文件中， 由 法国塞纳河流域管理局 （塞纳河水务局）提供。 在此文件中，我们可以找到与碳酸氢盐（HCO 3，单位为mgC/L）、pH值和温度（°C）相关的数据。 这些变量是在2006年10月测量的。 现在我们将了解如何使用插值将这些离散数据转换为连续数据。 实际上，我们会将点的形状文件转换为网格，其中每个像素都有一个质量值。 在这种情况下，当我们想要了解无气象点的水质时，内插可能很有用。
19.4.2. 数据准备#
在本练习中，我们将仅关注碳酸氢盐数据插值 。 我们想要衡量插值的质量，以选择最佳的插值方法。 我们将239个测量样本分成两个相等的部分。 第一部分将用于内插，第二部分将用于评估内插的质量。
为了随机选择50%的点，我们使用 处理收件箱 的菜单 Q地理算法-载体选择工具 中的 随机选择 模块。
在新窗口中，如 图 19.37 。
我们设置了要从中选择点“groundwater_oct_2006.shp”的层，对于方法，
我们选择了“所选要素的百分比”，该百分比为50%。
图 19.37 选择50%的测量点#
选择点后，我们通过右键单击该层并从菜单中选择“保存为.”来保存此选择。 不要忘记勾选“仅保存所选功能”框。 在“文件名”行中，我们将新层命名为“groundwater_ oct_2006_validation.shp”， 并将其存储在我们的工作目录中。新的点层出现在QGIS中。
我们现在选择这些点的对称组。我们将使用这些点来执行插值。 为此，我们通过右键单击该层并选择菜单“开放属性表”来打开初始层的属性表（之前的选择应该仍然有效）。 与所选点对应的线呈蓝色。在窗口顶部，我们选择菜单“反转选择”，然后获得前一个菜单的对称选择。 至于第一部分，我们保存这个新选择， 并将其命名为“groundwater_oct_2006_inclusion.shp”（不要忘记勾选“仅保存选定的要素”）。 我们现在有两个互补的数据集。
19.4.3. 插值的实现#
现在，我们使用 SAGA 地理信息系统模块对碳酸氢盐测量值（字段 “HCO3” ）进行插值。
有多种插值方法可用。我们将比较三种不同的方法：B样条逼近、逆距离加权内插（IDW）和薄板样条（全局）。
让我们从第一种方法开始。在菜单 SAGA 地理信息系统中，我们选择模块 B样条逼近 。
我们首先选择要插值的点层（“groundwater_oct_2006_interpolation.shp”）。
然后我们选择要插值的字段（“HCO 3”）并将分辨率更改为500。
研究范围越广，该值应该越高。在“Grid”行中，
显示插值的结果“HCO3_interpol_B-spline.tif”，如 图 19.38 。
图 19.38 使用B样条方法设置插值#
图4.39所示的结果出现在窗口中，使用“B样条逼近”方法， 如 图 19.39 （盆地的边界为橙色）。 我们可以改变它的风格以更好地展示它。
图 19.39 对塞纳河盆地中的碳酸氢盐进行插值#
有关该图形的彩色版本，请访问见www.iste.co.uk/baghdadi/qgis1.zip
我们现在执行相同的插值，但使用“逆距离加权插值”方法。相应的模块以相同的方式命名。 在新窗口中，我们选择要内插的层“groundwater_oct_2006_interpolation.shp”和要内插的字段（“HCO 3”）。 我们将该方法保留为默认“到顶点的反向距离”，其参数默认设置为2和1。 我们将“搜索范围”设置为“无搜索半径（全球）”，因为我们的研究区域足够宽。 我们可以为其他参数保留默认值，如 图 19.40 。
图 19.40 使用“IDW”方法设置插值#
我们将这个新的插值命名为“HCO3_interpol_idw.tif”。 结果应显示在QGIS中，使用“距离倒数加权”法对塞纳河流域的碳酸氢盐进行插值，如 图 19.41 。
图 19.41 距离倒数加权法#
有关该图的彩色版本，请参阅www.iste.co.uk/baghdadi/qgis1.zip
为了完成任务，我们使用“薄板样条（全局）”方法进行相同的工作。 我们使用同名模块。在模块的窗口中，我们总是设置要插值的层和要插值的场。 我们将结果命名为“HCO3_interpol_thin_plate.shp”，如 图 19.42 。
图 19.42 薄板进行插值设置样条法#
我们得到的结果，使用“薄板样条（全局）”方法对塞纳河盆地中的碳酸氢盐进行插值，如 图 19.43 。
图 19.43 碳酸氢盐进行插值#
有关该图的彩色版本，请参阅www.iste.co.uk/baghdadi/qgis1.zip
19.4.4. 内插评估#
我们现在有三种不同的插值，我们想保留赋能。 为了评估这些插值，我们将使用第二个数据集， 即我们为验证设计的数据集“groundwater_oct_2006_validation.shp”。 目的是计算每个验证点的真实测量值和内插值之间的差异。
我们使用名为“将栅格值添加到点”的 SAGA 地理信息系统模块将相应的插值与验证点关联起来。
我们在 SAGA 地理信息菜单“Vector <-> raster”中的“Processing Toolbox”中找到该模块。
在新窗口中，我们选择验证点“groundwater_oct_2006_ validation.shp”，
我们在行“Grid”中设置插值的栅格（它们的顺序在点的最终形状文件中相同），我们选择方法“最近邻居”，
以便仅保留每个点的相应像素，
并将结果命名为“inclusion_validation_temp.shp”，如 图 19.44 。
图 19.44 内插数据与实测数据比较#
在新层的属性表中，我们可以看到三个新字段。这些字段包含每个方法的插值。 这些字段的顺序与我们设置“将栅值格添加到点”模块时选择网格的顺序相同。
我们可以重命名这些字段，以使接下来的步骤更容易。
为此，我们使用 Processing Toolbox 的“重构字段”模块。
要重命名字段，只需双击“名称”列中的名称即可。
我们将这个新层保存为名称“inclusion_validation.shp”。
现在，我们在“inclusion_validation.shp”层中创建三个新字段。 在这些字段中，我们计算每个插值与测量值之间的差。 为此，我们打开层的属性表，并通过单击属性窗口顶部的“打开字段计算器”图标打开字段计算器。 我们勾选“创建新字段”框并将其命名为“d-thin”。 在这个领域中，我们计算每个点的测量值与通过“薄板样条”方法获得的插值之间的差。 在“输出字段类型”行中，选择“小数（真实）”，因为该字段将包含真实数字。 我们将“精度”设置为3（我们可以选择其他值）。 在“表达面板”中，我们写下表达“HCO 3”-“thin”。 我们可以使用“search”栏下的面板输入公式，如 图 19.45 。
图 19.45 现场计算器计算测量值和插值之间的差异#
单击“确定”后，新字段将出现在属性表中。 我们按照相同的步骤计算剩余的两个差异。 将这些字段命名为“d_idw”和“d_bspline”。 计算完成后，我们必须保存该层。 我们现在有三个新字段，其中包含测量值和插值之间的差异。
我们简单地通过计算每个差异的平均误差来评估插值的质量。 我们使用“Processing Toolbox”中的“数字字段的基本统计数据”。 我们首先设置要分析的层“inspinning_ validation.shp”。 然后我们设定我们想要评估的领域。 我们可以从“d-thin”字段开始。 为了完成任务，我们设置包含统计信息“stats_thin”的输出文件， 如 图 19.46 。
图 19.46 计算每个观测值和内插值之间差异字段的统计量#
我们对剩下的两个差异也做同样的事情。 我们分别将输出文件命名为“stats_idw”和“stats_bspline”。 结果存储在html文档中。 我们可以使用“Firefox 4”等网络浏览器打开此文档。 我们分别获得“薄板样条”、“idw”和“B样条”方法的以下平均误差-0.66、-0.72和-13.72。 最后，我们可以选择第一种方法作为我们的情况下的最佳方法。 我们必须注意，结果取决于我们在本练习开始时进行的随机选择（4.4.2.2部分）。 由于薄板样条法和“idw”法的误差非常接近，因此可以根据开始使用的随机选择进行倒置。
最后，我们必须坚持最终结果高度取决于所使用的插值方法。 用户应该正确定义问题，以便选择最相关的方法及其最相关的参数。 最困难的是方法的定义，而不是其实施。