1. 确定性空间插值

确定性插值方法将兴趣变量视为确定性变量，因此不考虑其随机性，其输出是确定性的值。插值方法主要基于数学公式或几何规则，不涉及概率或统计模型。

(1) 反距离加权法（ Inverse Distance Weighting, IDW ）

  • 原理：未知点的值由已知点的值加权平均得到，权重与距离成反比。
  • 公式：

    z^(x)=∑i=1nwizi∑i=1nwi,wi=1d(x,xi)p, \hat{z}(x) = \frac{\sum_{i=1}^n w_i z_i}{\sum_{i=1}^n w_i}, \quad w_i = \frac{1}{d(x, x_i)^p},

其中，Z(xi)Z ( x_i ) 是已知点的值，d(x,xi)d ( x, x_i ) 是未知点与已知点的距离，pp 是幂参数。

  • 特点：
    • 简单易用，计算速度快。
    • 对数据分布敏感，可能导致“牛眼效应”。

(2) 泰森多边形法（ Voronoi Diagram / Thiessen Polygons ）

  • 原理：将空间划分为多个多边形，每个多边形内的点值由最近的已知点决定。
  • 特点：
    • 适用于离散数据。
    • 插值结果呈阶梯状，不连续。

(3) 三角剖分线性插值法（ Triangulated Irregular Network, TIN ）

  • 原理：将已知点连接成三角形网络，在三角形内进行线性插值。
  • 特点：
    • 适用于不规则分布的数据。
    • 插值结果连续但不平滑。

(4) 自然邻域法（ Natural Neighbor Interpolation ）

  • 原理：基于 Voronoi 图，利用自然邻域关系进行插值。
  • 特点：
    • 适用于不规则分布的数据。
    • 插值结果平滑且连续。

2. 基于函数的插值方法

此类方法主要来自于计算机图形学，通过拟合空间曲面的函数来估计未知点的值。常见方法有：

(1) 多项式插值法（或拉格朗日插值）

  • 原理：通过构造基函数使得多项式在给定点处的函数值与指定数据匹配。此方法在曲面插值中被广泛应用，尤其在不规则网格数据拟合和快速近似场景中具有独特优势。
  • 常见方法：
    • 双线性插值：基于四个邻域点的加权平均，常用于图像处理和纹理映射，计算简单但平滑性有限。
    • 双三次插值：使用16个邻域点的加权平均，生成更平滑的曲面，广泛用于计算机图形学（如光照贴图）。
    • 张量积多项式插值：通过构造多项式基函数（如双二次、三次）进行高阶插值，适合规则网格数据。

(2) 径向基函数法（ Radial Basis Function, RBF ）

  • 原理：以数据点为中心构建基于距离的径向基函数，然后通过所有基函数的线性组合来计算未知点处的曲面值。在具体实施时，假设基函数的权重系数是未知的确定值，利用所有已知数据点求解权重系数的最优值。
  • 常见方法：
    • 高斯核：ϕ(r)=e−(ϵr)2\phi ( r ) = e^{- ( \epsilon r ) ^2}
    • 多项式核：ϕ(r)=1+(ϵr)2\phi ( r ) = \sqrt{1 + ( \epsilon r ) ^2}

(3) 样条插值法（ Spline Interpolation ）

  • 原理：使用分段多项式函数（ 如三次样条 ）拟合数据。具体来说，通过两个方向（如 uu 和 vv）的样条基函数组合来构造曲面：

    S(u,v)=∑i=0m∑j=0nPij⋅Bi,p(u)⋅Bj,q(v),\mathbf{S}(u, v) = \sum\limits_{i=0}^m \sum\limits_{j=0}^n \mathbf{P}_{ij} \cdot B_{i,p}(u) \cdot B_{j,q}(v),

    其中：
    • Pij\mathbf{P}_{ij} 是控制点网格（定义曲面形状）。
    • Bi,p(u)B_{i,p}(u) 和 Bj,q(v)B_{j,q}(v) 是 pp 次和 qq 次的 B 样条基函数。
    • 节点向量 U={u0,u1,…,um+p+1}U = \{u_0, u_1, \dots, u_{m+p+1}\} 和 V={v0,v1,…,vn+q+1}V = \{v_0, v_1, \dots, v_{n+q+1}\} 决定基函数的局部支撑域。
  • 常见样条：
    类型 公式/特点 适用场景
    均匀 B 样条曲面 节点等距分布，基函数周期性重复。 规则数据（如网格化点云）
    非均匀 B 样条（NURBS） 引入权重系数，支持更复杂的几何形状（如圆锥、球体）。 CAD建模、工业设计
    薄板样条（Thin Plate Spline） 最小化弯曲能量，全局平滑但计算成本高。 地质建模、医学图像重建
    双三次样条（Bicubic Spline） 三阶连续性（C2），基函数为分段三次多项式。 图像插值、纹理映射

(4) Hermite 插值

  • 原理：与上述方法仅关注节点处的函数值不同，Hermite 插值同时利用节点处的 “函数值” 和 “一阶导数信息” 来构造平滑的多项式曲面。它不仅要求曲面经过给定的点，还要求曲面在这些点处的一阶导数（斜率）与给定值匹配，从而实现更高的平滑性和物理意义上的自然过渡。
  • 常见方法：
    • 双 Hermite 插值：这通过 16 个控制点（包含位置和切向量信息）定义一个双三次多项式曲面。
    • 局部 Hermite 插值：在细分曲面或非均匀有理 B 样条（NURBS）中，Hermite 条件可用于局部调整曲面形状。
    • 基于物理的 Hermite 插值：通过引入力或能量约束（如弯曲能量最小化）优化曲面形状，同时满足 Hermite 边界条件。

(5) 特点

  • 计算效率高
  • 插值结果平滑且连续。
  • 适用于规则或不规则分布的数据。

3. 物理驱动的插值

  • 原理：基于物理模型驱动的曲面插值方法通过结合物理规律（如力学平衡、能量最小化、流体动力学等）实现插值，确保生成的曲面不仅满足数学平滑性，还符合实际物理行为。这类方法在工程仿真、医学建模、计算机图形学等领域广泛应用。
  • 主要方法：
    • 有限元方法：将曲面视为弹性薄壳或膜，通过能量（如弯曲能、拉伸能）最小化约束曲面形状，并利用有限元网格离散化求解。
    • 有限差分方法：将插值问题转化为偏微分方程（PDE）的边值问题，并在离散网格上求解该方程。此方法适用于规则网格数据，并假设场分布满足某种物理规律（如平滑性、能量最小化）。
    • 质点弹簧系统：将曲面建模为质点网格，质点间通过弹簧连接，受物理力（弹性力、阻尼力、外力）驱动达到平衡状态。
    • 流体力学：将曲面视为流体界面，通过求解 Navier-Stokes 方程模拟流动，利用速度场驱动曲面变形。
    • 能量最小化曲面：通过最小化特定能量泛函（ 如 Dirichlet 能量、Willmore 能量 ）生成平滑曲面。
    • 薄板样条：通过最小化曲面的弯曲能量实现全局平滑插值，也是一种 RBF 方法（因此无需网格），但强制要求曲面的二阶导数为零以保证刚性。
    • 基于物理的样条：在传统样条（如 B 样条、NURBS）中引入物理约束（如弹性力、惯性力），动态调整控制点。
  • 小结：
方法 适用场景 计算成本 物理一致性
有限元法 高精度工程仿真、复杂边界 高 高
薄板样条 散乱数据插值、全局平滑曲面 中 中
质点弹簧系统 实时动画、交互式建模 低 低
流体力学驱动 流体界面追踪、动态拓扑变化 极高 高
能量最小化曲面 数学平滑曲面生成、去噪 中-高 高
物理样条 参数化设计与物理模拟结合 中 中-高

4. 地统计插值方法

地统计插值方法基于统计学原理，考虑了数据的空间自相关性。地统计方法的插值结果不仅是确定值，而且包含不确定性量化结果。

(1) 经典克里金法（ Kriging ）

  • 原理：基于数据的空间自相关性进行插值，其中空间自相关性由变异函数（或协方差函数）建模。经典克里金认为变异函数（或协方差函数）的超参数是一个未知的确切值，可以在已知数据基础上，利用最小二乘法求解其最优值。
  • 常见方法：
    • 普通克里金（ Ordinary Kriging ）：假设均值未知但恒定。
    • 简单克里金（ Simple Kriging ）：假设均值已知且恒定。
    • 泛克里金（ Universal Kriging ）：假设均值符合趋势面分析。
    • 协同克里金（ Co-Kriging ）：引入协变量进行插值。
  • 主要步骤：
    • ① 模型假设： 经典克里金法假设观测数据满足以下条件：
      • 平稳性假设：数据的统计性质（如均值、方差、协方差）在空间上是平稳的。
      • 正态性假设：误差项 ϵi\epsilon_i 服从均值为0、方差为 σ2\sigma^2 的正态分布。
      • 无偏性假设：插值结果是无偏估计。
    • ② 克里金模型构建： 假设观测数据 Z={zi,xi}i=1n\mathbf{Z} = \{z_i, \mathbf{x}_i\}_{i=1}^n 满足以下模型：

      zi=f(xi)+ϵi,z_i = f(\mathbf{x}_i) + \epsilon_i,

      其中：
      • f(x)f(\mathbf{x}) 是目标函数（如地下水位、污染物浓度）。
      • ϵi\epsilon_i 是独立同分布的误差项，通常假设服从正态分布 N(0,σ2)\mathcal{N}(0, \sigma^2)。
    • ③ 协方差函数建模：协方差函数（又称半方差函数）描述了空间点之间的相关性强度与距离的关系，表示为 C(xi,xj)=Var[f(xi)−f(xj)]C(\mathbf{x}_i, \mathbf{x}_j) = \text{Var}[f(\mathbf{x}_i) - f(\mathbf{x}_j)]，常用形式包括：
      • 指数模型：C(d)=Ae−λd+σ2C(d) = A e^{-\lambda d} + \sigma^2，
      • 球形模型：C(d)=A(1−γd)2C(d) = A(1 - \gamma d)^2，
      • 高斯模型：C(d)=Aexp⁡(−d22α2)C(d) = A \exp\left(-\frac{d^2}{2\alpha^2}\right)，
        其中 dd 是两点间距离，A,λ,γ,αA, \lambda, \gamma, \alpha 等为待估计的超参数。
    • ④ 超参数的点估计（矩量法）：
      • 计算所有点对的半方差 sij=12E[(Zi−Zj)2]s_{ij} = \frac{1}{2}\text{E}[(Z_i - Z_j)^2]，得到拟合协方差函数所需的样本。
      • 通过最小二乘法拟合协方差函数的超参数，如 A,λ,γ,αA, \lambda, \gamma, \alpha 等。
    • ⑤ 权重计算：
      目标是找到一组权重 λi\lambda_i，使得预测值 z^∗\hat{z}^* 能够最小化误差：

      z^∗=∑i=1nλizi,\hat{z}^* = \sum\limits_{i=1}^n \lambda_i z_i,

      约束条件为：

      ∑i=1nλi=1,∑i=1nλixi=x∗,∑i=1nλiyi=y∗.\sum\limits_{i=1}^n \lambda_i = 1, \quad \sum\limits_{i=1}^n \lambda_i x_i = x^*, \quad \sum\limits_{i=1}^n \lambda_i y_i = y^*.

      通过拉格朗日乘子法求解线性方程组，得到最优权重 λi\lambda_i。
    • ⑥ 预测与不确定性量化
      • 预测值：

        z^∗=∑i=1nλizi.\hat{z}^* = \sum_{i=1}^n \lambda_i z_i.

      • 预测方差：

        Var(z^∗)=σ2(1−∑i=1nλi),\text{Var}(\hat{z}^*) = \sigma^2 \left(1 - \sum_{i=1}^n \lambda_i\right),

        反映插值的不确定性。若 Var(z^∗)<σ2\text{Var}(\hat{z}^*) < \sigma^2，说明插值有效；否则数据可能存在偏差或模型不适用。

(2) 经验贝叶斯克里金（ Empirical Bayesian Kriging, EBK ）

  • 原理： 经验贝叶斯克里金方法认为变异函数（或协方差函数）的参数不是一个确切的值，而是一个概率分布，因此在已知数据基础上通过最大边缘似然法对参数进行概率推断。经验贝叶斯克里金通过融合贝叶斯统计与空间插值，显著提升了传统克里金方法的鲁棒性和可解释性。它在需要同时估计参数和量化不确定性的场景中具有独特优势，但需权衡计算成本与模型复杂度。对于小样本或高噪声数据，EBK 是一种强有力的工具。
  • 主要步骤：
    • ① 克里金模型构建： 假设观测数据 Z={zi,xi}i=1n\mathbf{Z} = \{z_i, \mathbf{x}_i\}_{i=1}^n 满足以下模型：

      zi=f(xi)+ϵi,z_i = f(\mathbf{x}_i) + \epsilon_i,

      其中：
      • f(x)f(\mathbf{x}) 是目标函数（如地下水位、污染物浓度）。
      • ϵi\epsilon_i 是独立同分布的误差项，通常假设服从正态分布 N(0,σ2)\mathcal{N}(0, \sigma^2)。
    • ② 协方差函数建模： 定义空间协方差函数 C(xi,xj)=Var[f(xi)−f(xj)]C(\mathbf{x}_i, \mathbf{x}_j) = \text{Var}[f(\mathbf{x}_i) - f(\mathbf{x}_j)]，常用形式包括：
      • 指数模型：C(d)=Ae−λd+σ2C(d) = A e^{-\lambda d} + \sigma^2，
      • 球形模型：C(d)=A(1−γd)2C(d) = A(1 - \gamma d)^2，
      • 高斯模型：C(d)=Aexp⁡(−d22α2)C(d) = A \exp\left(-\frac{d^2}{2\alpha^2}\right)，
        其中 dd 是两点间距离，A,λ,γ,αA, \lambda, \gamma, \alpha 为待估计参数。
    • ③ 超参数的贝叶斯估计： 通过最大化边缘似然函数（Marginal Likelihood）估计协方差函数的超参数和噪声方差：

      p(Z∣A,λ,σ2)=∏i=1n12πσ2exp⁡(−(zi−z^i)22σ2), p(\mathbf{Z} | A, \lambda, \sigma^2) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(z_i - \hat{z}_i)^2}{2\sigma^2}\right),

      其中 z^i\hat{z}_i 是克里金预测值。通过共轭梯度法等数值优化算法求解参数的后验分布。
    • ④ 预测与不确定性量化： 对任意未观测点 x∗\mathbf{x}^*，其预测值为：

      z^∗=E[f(x∗)∣Z]=∑i=1nλi(zi−zˉ)+zˉ,\hat{z}^* = \mathbb{E}[f(\mathbf{x}^*) | \mathbf{Z}] = \sum\limits_{i=1}^n \lambda_i(z_i - \bar{z}) + \bar{z},

      其中 λi\lambda_i 是拉格朗日乘子，zˉ\bar{z} 是观测均值。不确定性通过预测方差量化：

      Var(z∗∣Z)=σ2(1−∑i=1nλi).\text{Var}(z^* | \mathbf{Z}) = \sigma^2 \left(1 - \sum\limits_{i=1}^n \lambda_i\right).

5. 机器学习插值方法

近年来，机器学习方法逐渐应用于空间插值。

(1) 高斯过程回归（ Gaussian Process Regression, GPR ）

  • 原理：高斯过程利用核函数（可视为空间场中协方差函数的一种泛化）描述数据的相关性，通过贝叶斯推断获得核函数超参数的后验分布，进而为测试点预测兴趣变量的概率分布。高斯过程模型是对克里金方法在多维空间中的泛化，或反之，克里金法是高斯过程在二维或三维空间中的实例。两者尽管在数学上是等价的，但在超参数的假设和推断方法上存在着本质区别（参见克里金和高斯过程的对比）。
  • 特点：
    • 提供插值结果的不确定性估计。
    • 计算复杂度较高。

(2) 神经网络插值

  • 原理：使用神经网络（ 如多层感知机、卷积神经网络 ）学习数据的空间分布。
  • 特点：
    • 适用于高维数据。
    • 需要大量训练数据。

6. 总结

方法类型 代表方法 特点
确定性插值 IDW、泰森多边形、TIN、自然邻域法 简单易用，计算速度快
基于函数的插值 多项式插值、径向基函数、样条插值 插值结果平滑且连续
模型驱动的插值 有限元插值、质点弹簧系统 能够满足物理约束
地统计插值 克里金、经验贝叶斯克里金 考虑空间自相关性，提供不确定性估计
机器学习插值 高斯过程回归、神经网络 适用于高维数据，需要大量训练数据
混合插值 克里金与 IDW 结合、机器学习与地统计结合 结合多种方法，提高插值精度

选择插值方法时，需根据数据特性、应用场景和计算资源进行权衡。