空间插值是空间分析的重要基础^{[1, 2, 3]}，空间采样的不均匀性是空间插值算法所面临的关键问题之一^{[4, 5]}。结构保持的空间插值应能兼顾样本空间分布不均与空间结构约束，自适应生成待插控制点并进行估值。常用的空间插值方法多直接基于不规则分布样点进行规则格网插值，易导致样点分布特征的均匀化，进而引起空间结构和统计特征的畸变。如基于距离的插值点估值易出现边缘效应和牛眼现象^{[4, 6]}；函数逼近法对数据噪音敏感而出现局部最小^[7]；地统计学方法基于空间自相关特性进行插值估计和误差不确定性分析，但要求数据分布满足相应的假设条件^{[8, 9, 10]}。因此可从样点空间分布入手，基于样点代表性及分布特征构建几何结构约束的空间插值算法^{[11, 12]}。其策略可归纳为：①根据样本点分布特征生成非均匀的几何空间结构；②构建几何结构约束下的插值控制点生成规则；③建立顾及几何与结构约束的插值控制格网或插值函数。

Voronoi-邻域(简称V-邻域)结构作为一种几何优化结构，可有效表征空间样点分布的结构关系，并能唯一确定不同样点间的空间影响域边界。Voronoi-Delaunay的对偶特性及空间剖分特性提供了样点空间分布范围划分的结构约束规则，可作为兼顾几何结构约束插值点自适应生成的基础^{[13, 14]}。本文基于Voronoi-Delaunay剖分的几何特性作为空间插值的结构约束限定，利用多次迭代Voronoi实现空间影响域逐层划分。基于高阶V-邻域的几何剖分特性，建立基于自然邻域的控制点自适应生成规则及结构约束的空间插值算法。基于中国气象台站日均气温数据进行插值试验，给出了不同层次插值对空间统计量的影响，并与常用插值方法进行了结果对比，讨论了该方法在空间参数特征保持等方面的优势。

设Ω是任意维度空间Rⁿ上的凸空间，给定有限点集S={p₁,p₂,…,p_n}∈Ω，当存在i≠j,p_i≠p_j 时，则对任意节点x∈S，均可定义一个子空间V_x 满足：

其中，d(p_i,p_j)为点p_i、p_j间的欧氏距离；V_x为节点x的一个Voronoi区域。

在点x的Voronoi区域V_x上增加任意新节点y∈Ω，则可生成新的Voronoi区域V_y，此时V_y与原有的点x的Voronoi区域V_x的重叠部分会形成新的Voronoi区域，称之为2层Voronoi区域V_xy^{[14, 15]}，数学表达为：

Voronoi的剖分特性表明，对于新的Voronoi区域V_xy自然邻近点坐标可由所划分的新Voronoi区域V_xy的面积与原Voronoi区域V_x的面积比直接确定。更高层次的Voronoi区域同样可由式(2)迭代获得。

Voronoi区域的空间剖分特性可为空间插值提供结构性空间约束^[16]。对于任意点x，其自然邻近点定义为与它具有共同Voronoi区域边界的节点，即对应Delaunay三角形中与之有连线的节点^[17]。自然邻近点定义了几何特性上距离某点最近的节点系，其节点数目和位置仅取决于节点的局部分布，可作为几何结构约束表征的理想控制点。以任意K层-Voronoi的自然邻近点为控制点，基于自然近邻规则^[18]的插值定义如下：

其中，f(x)为待插点x处的插值结果；i是点x的自然邻近点序号；f_i为给定的样本点值；φ_i(x)为对应节点x_i的插值基函数。自然邻点插值方法中φ_i(x)为点x的自然邻近坐标的线性插值，常依据插值前后空间几何剖分的面积比确定。而Voronoi区域的面积大小与节点分布的密度成反比关系，生成单元的疏密程度取决于节点分布的密度^{[13, 14]}。

基于空间分布特征的迭代插值过程如图 1所示。

图 1 基于空间分布特征的迭代插值过程 Fig. 1 Iterative Interpolation Process Based on Data Spatial Distribution

图选项

以基于原始点生成的K层-Voronoi的自然邻近点作为空间插值的控制点，即基于原始样点构建Delaunay三角网，各三角形的空间范围反映了已知数据的几何分布特征，并作为样点分布的结构约束初始结构。基于Delaunay和Voronoi的对偶关系，将Delaunay三角网中三角形的质心作为下一层Voronoi区域的自然邻近点。以此类推，可逐层生成K层-Voronoi自然邻近点的坐标。Voronoi空间剖分的不规则性导致插值控制点生成的不规则性。在任意K层，原始数据分布 稀疏的区域，三角形的分布个数较少，而分布较密集的地区，生成的待插点相对集中。考虑控制点与插值分辨率的关系，基于Voronoi插值的分辨率由迭代次数控制，而非传统插值算法的单元大小。迭代次数增加时，原始样点密集区插值点增加越多，而稀疏区的插值点增加相对较慢，从而保持对原有空间样本分布结构的表征。

基于前述K层-Voronoi的定义，可以将每次生成的K层-Voronoi的自然邻近点作为φ_i(x)的基础构建插值方法。由于反距离加权函数空间上更符合地理学第一定律，且在局部区域其迭代多次,插值稳定性较好，更适合于新待插点值的估计^[19]。 对K层-Voronoi生成的1，…，K层自然邻近节点逐层迭代，并进行反距离加权获得φ_i(x)。由于K层-Voronoi的生成具有样点分布自适应性、规则唯一性和局部保形性，因此基于K-Voronoi插值同样具备上述特点，从而为构造几何结构约束的空间插值提供了基础。

2.3 结构自适应的插值算法

基于K层-Voronoi剖分构建基于V-邻域迭代的结构保持空间插值算法(见图 2)。首先对原始样本点进行第一次Delaunay剖分，以落入边界的三角形质心为待插点，对应的三角形顶点为控制点进行反距离加权插值。将估值后的质心追加至原样本点，然后对追加后的新样本点进行第二次Delaunay剖分，同样对边界内的三角形质心进行估值，再将新的待插点追加到该次插值样本点中继续下层空间剖分，超出迭代次数约束后算法终止，输出并显示插值结果。

图 2 算法流程 Fig. 2 Flow Chart of Algorithm

图选项

日平均气温的空间插值要求在站网分布约束下较好地再现气温分布的空间结构^{[6, 7, 20]},为此选用WMO(World Meteorological Organization)提供的中国主要城市气象站2013年6月1日平均气温数据进行算法评估。该数据在中国大陆共377个站位，空间分布不均，东部数据明显较西部数据分布密集，尤其是青藏高原地区气象站点甚少。通过 计算原始数据及各层插值结果的基本统计量、空间方差和空间协方差，以验证算法的结构保持特性和稳定性,进而将其与样条函数法(Spline)、自然邻点法(natural neighbor interpolation,NNI)、克里金法(Kriging)和反距离权重法(inverse distance weighted interpolation,IDW)进行对比，分析不同方法插值结果的统计参数特征。 3.2 算法自身空间保型性和稳定性验证

基于原始测站数据进行多层V-邻域结构自适应插值，各次迭代的结果见图 3。三次迭代后,插值点已具有较高的空间分布密度,各层插值结果均保持了原始数据东密西疏的分布特征结构，较好地再现了气温变化的整体构型。受中国复杂地形条件的影响，温度总体上呈阶梯性变化。原始数据和不同迭代次数插值数据的常规统计量和空间统计量分别见表 1和图 4。

图 3 不同迭代次数的自适应插值结果 Fig. 3 Adaptive Interpolation Results in Different Iterations

图选项

图 4 不同插值方法的空间统计量对比 Fig. 4 Comparison of Spatial Statistical Index with Different Interpolation Methods

图选项

不同迭代次数的插值结果的最小值和最大值与原始数据一致,而均值和标准差随着插值层数的增加略有下降，这与迭代插值后数据分布的趋中性有关。变异系数上，各单次迭代插值具有一致的变异系数，且与原始数据的变异系数相差小于0.02。三次迭代无论在方差结构还是协方差结构上,均接近于原始数据，可很好地保持原有数据的空间分布特征。插值数据随着迭代次数的增加趋于平滑，导致其方差和协方差略有降低。总体上，本文方法的插值结果与原始数据的空间统计量对应性好，保持了原始数据的统计特征。

将3次迭代插值与Spline、NNI、Kriging和IDW方法进行算法对比(见表 1)。为方便对比，将所有方法的插值结果统一至10 m分辨率的格网，由于本文算法的插值结果具有不均匀性，采用不同格网内部所有插值点的均值作为格点值(见图 5)。

图 5 不同方法空间型的对比 Fig. 5 Comparison of Spatial Mode with Different Methods

图选项

从图 5可以看出,Spline、Kriging和IDW方法在样点稀疏分布区域存在明显的条带分层，结构过渡连续性差。NNI和本文方法均是利用V-邻域结构的插值方法，在结构上也具有更高的相似性，但NNI仅利用V-邻域结构的顶点作为控制点进行一次插值，不能支持根据分布特征及插值分辨率需求的自适应控制点生成^[13]。

本文算法通过限定待估点的影响范围，使权重分配符合原始数据的空间结构特征，稀疏区域的结构更精细，高、低值区域过渡更连续。统计指标对比也显示本文方法可很好地保持原始数据的空间分布信息，而其他方法的统计参量变化总体上更为显著。三次随机5%原始样点的交叉验证也显示本文算法具有较好的精度，可在样本数据结构特征保持和插值精度间取得较好平衡。 4 结 语

现有的插值方法多在理想数学假设的基础上构建，对样点数据空间分布结构的约束考虑不足，难以兼顾全局与局部精细结构的高精度再现，并很难保证对原始数据空间统计量的保持。本文基于K层-Voronoi结构建立了新的控制点生成方法，并利用Voronoi的局部保形特性和反距离加权实现了基于V-邻域结构约束的空间插值方法。本文算法的特性如下：①可以反映样本点分布结构与密度的层次自然近邻；②通过迭代插值获得不同分辨率插值数据的精细结构；③具有较好的稳定性、结构保形性和空间约束集成特性。

由于本文方法最终插值的分辨率与采样点分布相关，不同样点的分布可能导致插值结果的点位分布存在一定差异。该插值方法的层次在理论上可无限外推，因此可以根据样点分布特征及插值应用需求确定合适的插值层数；V-邻域结构作为常用的几何结构，具有高维的可拓展性，进而可支撑三维及高维数据插值；其在数据结构上的通用性和易集成性也利于V-邻域结构约束的空间插值方法与空间统计模型集成。