CN110555189B - 一种基于反向计算思维的空间插值方法 - Google Patents

Abstract

本发明公开了一种基于反向计算思维的空间插值方法，属于空间插值模型构建的技术领域。本发明基于空间物体自相关定律，反向思维地提出了构建采样点与未知格网点约束方程、以及格网点和格网点之间的约束方程，通过预处理共轭梯度法以及最小二乘原理对方程进行求解，得出格网点插值属性值。该方法综合考虑采样点对格网点和格网点之间的影响，并解决了传统经典插值方法中，在采样点较少时，单个采样点影响范围较大问题，极大程度上提高了空间插值精度。

Description

一种基于反向计算思维的空间插值方法

技术领域

本发明属于空间插值模型构建的技术领域，涉及采样点与邻近插值格网点的约束关系构建，以及插值格网点之间的约束关系构建，特别是提出了通过反向计算思维的空间插值方法。

背景技术

空间插值是在一定区域内，基于有限采样点数据来计算区域内任意位置对应变量值的一项地理空间分析技术。借助空间插值技术可以将某种空间变量由点扩展面，得到空间变量的空间分布，因此在地理、环境、生态等领域都有非常广泛的应用。插值精度是决定空间插值技术应用价值的重要因素，因此许多学者针对不同空间插值技术的精度特征展开讨论，同时也不断提出、优化新的空间插值技术。特别地，考虑到数据采样的成本，如何基于有限的采样点获取精度较高的插值结果，已经成为空间分析技术领域的重要课题。

根据利用采样点信息方式的不同，可将传统的插值方法分为确定性插值方法和地统计插值方法，其中确定性插值是基于一定的数学原理或局部相似性规则直接利用采样点的值计算其他位置的属性值，主要包括反距离插值法、样条函数插值法、多项式插值法等；而地统计插值则是要根据采样点先分析变量值的空间变化特征，依据半变异函数计算其他位置的属性值，主要包括克里格家族中的各种插值方法。对于以上列举的插值方法，其提出的时间较早，已经在许多领域中得到应用，且已经集成到一些地理空间分析软件中，便于用户使用。但是这些方法存在的缺点也是明显的，例如：反距离加权插值过分依赖小邻域内采样点的信息，当研究区域内采样点信息较少时插值效果一般较差；样条函数插值法过于追求插值结果的光滑度，会在采样点较少的区域形成异常的凹凸特征；多项式插值将研究变量的区域分布看作一个用多项式描述的数学函数，其物理意义不明显，且当多项式次数较高时采样点的误差会带来较大的插值误差；而对于克里格插值法，尽管其精度较高，但是插值过程复杂，需要计算若干插值参数，限制了普通用户的使用。

为解决以上各种插值方法的问题，我国学者岳天祥研究员等基于曲面论基本定理开创性地建立了高精度曲面建模方法。高精度曲面建模方法将研究变量在区域内的分布理想成一个光滑的数学曲面，则根据微分几何学中的曲面论基本定理，在曲面中任意位置的微分曲面内不同的变量之间存在一定的约束关系，即高斯方程组。高精度曲面建模方法就是根据这种约束关系，再加上已知采样点的值的约束，将空间插值问题转换为一个大型线性方程组的求解问题。高精度曲面建模方法巧妙地将研究区域内不同位置的点(包括采样点和要计算的非采样点)联系到一起，非常适合对稀疏样本下的采样点进行空间插值，且精度优于经典插值方法，近年来已经在气候、土壤、大气、生态等领域中得到成功应用。然而，高精度曲面建模方法仍然存在一些不可解决的问题：第一是该方法将变量在空间中的分布理想化为光滑的数学曲面，这在一些工业设计中的插值问题中是合适的，但是对于绝大部分地理、环境、生态等学科中的插值问题都是不适用的；第二是该方法依据曲面论基本定理，本质上是考虑了邻近点之间数学上的约束关系，而忽略了邻近点之间的自相似性，导致在对地理变量进行空间插值时过分关注结果的光滑性，结果的合理性反而受到影响；第三个问题，虽然高精度曲面建模方法经过多年发展与更新，但是该方法相对复杂，涉及到的模型参数较多，且需要有一定的数学功底方能理解，这在一定程度上限制了该方法的进一步应用。针对上述问题，在高精度曲面建模方法构建思路的基础上，必须要发展一种新的空间插值技术，新的插值技术首先应当充分顾及变量的自相似性，即位置越近的变量值越相似，从而保证插值技术在地理、环境、生态等领域的应用广度；其次，新的插值方法应当将变量的局部相似性与区域变化的趋势特征结合，基于有限的采样点获取最优的插值效果；最后，应当简单易操作，尽可能避免过多的复杂参数，从而使新技术易于被大众用户接受。

发明内容

本发明为了弥补现有插值方法，在采样点较少时，插值误差较大的问题，面向空间采样点数据，提供了一种基于反向计算思维的空间插值方法，该方法可以显著提高建模精度。

为实现上述发明目的，本发明采用的技术方案如下：

一种基于反向计算思维的空间插值方法，具体步骤如下：

S1：确定采样点的坐标属性和插值属性，获得采样点的三元坐标；根据采样点和空间插值分辨率确定空间插值区域，并将插值区域划分为若干网格单元；

S2：基于反距离加权插值法，构建采样点附近待插值格网点的未知插值属性值与采样点的已知插值属性值的约束方程：

其中p_i为第i个采样点插值属性值，s_i表示第i个采样点附近待插值格网点插值属性值，W_i表示第i个采样点附近待插值格网点所占插值权重；i＝1，…，n，n为采样点个数，j＝1，…，m，m表示采样点附近待插值格网点的个数；

计算各个待插值格网点所占插值权重；

S3：基于反距离加权插值法，构建当前待插值格网点s_o的未知插值属性值与其附近插值格网点s_r的未知插值属性值的约束方程：

其中为，表示T_o为第o个待插值格网点附近插值格网点所占插值权重；o＝1，…g，g为待插值格网点个数，m′为插值格网点s_r的个数；

计算各个待插值格网点附近插值格网点所占插值权重；

S4：依据最小二乘法则，将第二步和第三步的约束方程联合成的大型线性方程组转化成最小值问题：

其中，s表示待插值格网点插值属性值的集合，k为迭代次数，λ为权重因子，W为采样点附近各待插值点所占插值权重的集合，T为待插值格网点附近插值格网点所占插值权重的集合，p为采样点插值属性值的集合；

S5：利用预处理共轭梯度法求解所述大型线性方程组，求得待插值格网点的未知插值属性值；

S6：将待插值格网点的插值属性值输出为TIFF格式的DEM文件。

进一步地，所述步骤S1中，遍历所有采样点，寻找采样点数据的最大与最小经纬坐标，确定采样点数据的外包矩形，将采样点外包矩形拓宽二分之一分辨率大小作为空间插值区域范围；然后将空间插值区域范围依据分辨率划分为m行和n列，得到m*n个网格单元，作为待插值格网点的计算单元。

进一步地，所述步骤S2中，计算各个待插值格网点所占插值权重的具体过程为：

S21：根据插值区域起始位置坐标、采样点的三元坐标以及插值分辨率，计算采样点所在的位置，即采样点的行号和列号；

S22：根据采样点的行号和列号计算采样点周围的待插值格网点的行号和列号；

S23：根据待插值格网点的行号和列号，结合插值区域起始点坐标及插值分辨率，计算待插值格网点的坐标；

S24：分别计算每个待插值格网点与采样点之间的距离；

S25：根据待插值格网点与采样点之间的距离，计算每个待插值格网点所占插值权重。

进一步地，所述步骤S3中，计算各个待插值格网点附近插值格网点所占插值权重的具体过程为：

S31：根据划分的网格单元，逐行构建每个待插值格网点的约束方程；根据当前处理的待插值格网点的行号和列号，判断待插值格网点周围采样点是否满足推出待插值格网点的插值属性计算需求，如果满足则转入步骤S4，否则执行步骤S32；

S32：基于待插值格网点的位置，确定构建约束方程的周围格网点m′个数，取值为3、5和8；计算周围格网点到待插值格网点距离；

S33：根据待插值格网点与周围格网点之间的距离，计算每个周围点所占插值权重。

本发明不同于传统经典插值方法根据距离利用已知采样点对未知格网点进行插值，而是反向思维地提出了构建采样点与未知格网点约束方程、以及格网点和格网点之间的约束方程，通过预处理共轭梯度法以及最小二乘原理对方程进行求解，得出格网点插值属性值。本发明方法综合考虑采样点对格网点和格网点之间的影响，并解决了传统经典插值方法中，在采样点较少时，单个采样点影响范围较大问题，极大程度上提高了空间插值精度。本发明具有以下优势：

(1)将变量的局部相似性特征与在区域内的趋势变化特征通过构建待插值点相互之间的约束关系很好地结合起来，在保证插值结果合理性的前提下，能够在稀疏采样下获得较高的插值精度，使得本方法特别适用于地理、环境、生态等领域内不同研究变量的高分辨率插值问题。

(2)原理简单易懂，便于软件实现与操作。本发明提出的新型插值方法不涉及复杂的模型参数，简单的原理让本领域技术人员普遍都能掌握。

(3)可扩展性好。本发明提出的新型插值方法实际给出了一种构建新型插值方法的范式，可以根据研究变量的分布特征构建不同的约束方程，从而实现最优化插值。

附图说明

图1为本发明方法的流程示意图。

图2为本发明实施步骤中空间插值名称解释示意图。

图3为本发明实施步骤中插值范围确定示意图。

图4为本发明实施步骤中采样点、插值格网点约束示意图。

图5为本发明实施步骤中插值格网点之间约束示意图。

图6为本发明实施例中训练点、验证点示意图。

图7为本发明实施例中基于不同插值方法的格网DEM，(a)本发明插值方法构建的DEM，(b)反距离加权插值法构建的DEM，(c)克里金插值法构建的DEM，(d)样条插值法构建的DEM。

图8为本发明实施例中基于不同插值方法的格网DEM山体阴影，(a)本发明插值方法构建的DEM山体阴影，(b)反距离加权插值法构建的DEM山体阴影，(c)克里金插值法构建的DEM山体阴影，(d)样条插值法构建的DEM山体阴影。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细说明。

首先对本实施例中涉及的名称进行如下解释：

采样点：在研究区域内采集的包含具有真实地物地理位置和高程属性的点，称为采样点。如图2点S0所示。

待插值点：通过已知采样点插值属性信息，计算未知点插值属性的点，称为待插值点。本发明将栅格格网的四个顶点作为待插值点，如图2点I1、I2、I3、I4所示。

行号：是栅格数据的像元存储索引，每个像元都有相应行、列号，可通过像元行号，找到该像元在栅格数据中的对应行，如图2所示。

列号：是栅格数据的像元存储索引，每个像元都有相应行、列号，可通过像元列号，找到该像元在栅格数据中的对应列，如图2所示。

起始点：栅格数据读取的起点，一般为空间插值范围的左上角坐标点。如图2点O所示。

本实施例的空间插值方法的流程图见图1，具体步骤如下：

步骤一、初始化：

1.a根据采样点和插值分辨率确定空间插值范围。根据插值分辨率和采样点的最小经纬度坐标(x_min，y_min)(如图3采样点S1所示)和最大经纬度坐标(x_max，y_max)(如图3采样点S2所示)，确定空间插值范围的最小经纬度坐标(X_min，Y_min)(如图3点O1所示)和最大经纬度坐标(X_max，Y_max)(如图3点O2所示)，计算公式如下：

X_min＝x_min-l_cell/2

X_max＝x_max+l_cell/2

Y_min＝y_min-l_cell/2

Y_max＝y_max+l_cell/2

1.b通过范围和分辨率计算插值区域网格的行列数，将空间插值范围划分文若干格网单元，按照如下公式计算：

行数＝Int((x_max-x_min)/l_cell)

列数＝Int((y_max-y_min)/l_cell)

其中x_min、x_max为插值范围的起始和终止经度坐标，y_min、y_max为插值范围的起始和终止纬度坐标。l_cell为分辨率长度。

1.c确定采样点坐标属性与插值属性，用三元坐标表示(x，y，h)。

步骤二、构建采样点与插值格网点之间的约束方程：

2.a判断未读采样点是否为空，如果是，转到步骤三，否则读取采样点的三元坐标(x，y，h)，并执行2.b；

2.b根据插值区域起始位置坐标、读取的采样点的坐标以及插值分辨率，按照如下公式计算采样点所在的位置，即行号和列号；

其中x₀、y₀为插值区域起始点坐标，x_i、y_i为采样点坐标。l_cell为分辨率长度。

2.c根据2.b计算出的采样点(如图4点S0所示)的行号和列号，计算采样点周围的4个待插值格网点(如图4点I1、I2、I3、I4所示)的行号和列号。

2.d根据2.c计算得到的四个待插值格网点的行、列号，结合插值区域起始点坐标及插值分辨率，计算四个待插值格网点的坐标；

x＝jl_cell+x₀

y＝y₀-il_cell

其中i为行号，j为列号，x₀、y₀为插值区域起始点坐标，l_cell为分辨率长度。

2.e分别计算待插值格网点与采样点之间的距离，计算距离公式为：

其中xi、yi为采样点坐标，x_j、y_j为插值点坐标。

2.f列出用四个待插值格网点计算采样点插值属性值的线性表达式：

其中p为采样点插值属性值，s_j为插值点插值属性值，W_j为各插值点所占插值权重，m为参与约束方程构建的插值格网点个数。

2.g计算上步骤中的各个权重，按照如下式子计算：

其中d_j为插值点到采样点距离的平方。

2.h组成线性方程组，结束步骤二。

步骤三、构建插值格网点之间的约束方程：

3.a依据步骤1.b划分的格网，逐行构建每个插值格网点的约束方程。根据当前处理的待插值格网点的行、列号，判断插值格网点周围采样点是否满足推出插值格网点的插值属性计算需求，是则转入步骤四，否则执行步骤3.b；

3.b基于确定的待插值格网点，确定其周围、距其最近的待插值格网点，根据位置不同，寻找的待插值格网点数量不等，可能为8个如图5点S3所示；可能为5个如图5点S2所示；也可能为3个如图5点S1所示。

3.c计算格网点到待插值格网点的距离，格网点到待插值格网点的距离一般有两种，一种为距离等于格网分辨率，另一种为距离等于根号下2倍的格网分辨率。

3.d根据距离，计算权重，权重计算和步骤2.g相同。

3.e列出表达式

其中s₀为插值点插值属性值，s_j为相邻插值点插值属性值，T_j为各相邻插值点所占插值权重。m为参与约束方程构建的待插值点附近插值格网点个数。

3.f转到步骤3.a，继续判断；

3.g组成线性方程组，并结束步骤三：

步骤四、依据最小二乘法则将第二步和第三步构建的方程联合成大型线性方程组：

4.a将步骤二和步骤三得到的线性方程组转化为如下给定条件下的最小值问题，

4.b基于最小二乘原理，可将4.a中的最小值问题变为如下：

可被缩写为

其中

为大型对称正定稀疏矩阵，

为一维向量。

步骤五、调用预处理共轭梯度法求解步骤四得到的大型线性方程组。

5.a生成预处理矩阵，具体方法是λ取默认为1，预处理矩阵M取以上方程组系数矩阵A(A为4.b中

)中对角线元素。

5.b初始化，x₀，k＝0，r₀＝b-Ax₀，计算第一步原始残差r₀，设x₀为n阶零矩阵，b(b为4.b中

)，A为以上方程组系数矩阵。

5.c求解Mz_k＝r_k，并令k＝k+1，如果k＝1，转5.d，否则转5.e。

5.d计算预处理残差z₀＝M^-1r₀，初始化共轭向量p₁＝z₀。

5.e计算第k步预处理残差

搜索方向标量

搜索方向p_k＝z_k-1+β_kp_k-1。

5.f计算第k步标量

第k步解向量x_k＝x_k-1+a_kp_k；第k步原始残差r_k＝r_k-1-a_kAp_k，并判断r_k是否为0，如果不为0则转到5.d，否则执行5.g。

5.g令x＝x_k，完成求解，输出解向量。

步骤六、输出结果，具体过程如下：

6.a.依次输出各文件头信息：包括列数、行数、西南角格网单元纵坐标、西南角格网单元纵坐标、格网间距、无效数据区域值。

6.b.将求解向量中的结果按照行列号依次输出，当输出元素个数等于结果文件列数时，另起一行，并重新计数。

6.c.求解向量结果全部输出，以TIFF格式输出为DEM文件。

已有的经典插值方法在对采样点的依赖、插值精度和方法操作的简易等方面，或多或少都存在一定的缺陷。为验证本发明方法的实用性，本实施例选取了地形起伏明显区域(能体现出插值精度差异)作为研究区，该研究区为山地地区，平均海拔1208米；相对高差104米；地形起伏明显。研究区内实验高程点各数10851个。本实施例通过与经典插值方法对比，检验本发明插值方法的精度与特有优势。

步骤1：分离训练点和验证点，训练点与验证点分离比为99：1，如图6所示。

步骤2：训练点插值，利用本发明方法与经典插值方法对训练点进行插值，得到格网DEM，如图7所示。

相对于三种经典插值方法，通过插值得到的DEM山体阴影对比(如图8)，可看出本发明插值方法在光滑性方面具有一定优势。

步骤3：精度检验，验证点提取高程值，验证点提取插值得到的格网DEM高程值与原始高程值取差的绝对值。通过与经典插值方法比较(表1)，可看出本插值方法在插值精度方面，具有较高的精度。

表1精度误差统计

结果DEM	本发明方法	反距离加权	克里金插值	样条法
					验证点数量	109	109	109	109
误差最小值	0.005981	0	0	0
					误差最大值	3.187012	3.562012	4.69873	4.753906
误差总和	68.705566	83.145264	76.683105	78.28479
					误差均值	0.630326	0.762801	0.703515	0.718209
误差标准差	0.582766	0.78877	0.661998	0.658124

Claims (4)

1.一种基于反向计算思维的空间插值方法，其特征在于，具体步骤如下：

S1：确定采样点的坐标属性和插值属性，获得采样点的三元坐标；根据采样点和空间插值分辨率确定空间插值区域，并将插值区域划分为若干网格单元；

S2：基于反距离加权插值法，构建采样点附近待插值格网点的未知插值属性值与采样点的已知插值属性值的约束方程：

其中p_i为第i个采样点插值属性值，s_i表示第i个采样点附近待插值格网点插值属性值，W_i表示第i个采样点附近待插值格网点所占插值权重；i＝1，···，n，n为采样点个数，j＝1，···，m，m表示采样点附近待插值格网点的个数；

计算各个待插值格网点所占插值权重；

S3：基于反距离加权插值法，构建当前待插值格网点s_o的未知插值属性值与其附近插值格网点s_r的未知插值属性值的约束方程：

其中为，表示T_o为第o个待插值格网点附近插值格网点所占插值权重；o＝1，···g，g为待插值格网点个数，m′为插值格网点s_r的个数；

计算各个待插值格网点附近插值格网点所占插值权重；

S4：依据最小二乘法则，将第二步和第三步的约束方程联合成的大型线性方程组转化成最小值问题：

其中，s表示待插值格网点插值属性值的集合，k为迭代次数，λ为权重因子，W为采样点附近各待插值点所占插值权重的集合，T为待插值格网点附近插值格网点所占插值权重的集合，p为采样点插值属性值的集合；

S5：利用预处理共轭梯度法求解所述大型线性方程组，求得待插值格网点的未知插值属性值；

S6：将待插值格网点的插值属性值输出为TIFF格式的DEM文件。

2.根据权利要求1所述的一种基于反向计算思维的空间插值方法，其特征在于，所述步骤S1中，遍历所有采样点，寻找采样点数据的最大与最小经纬坐标，确定采样点数据的外包矩形，将采样点外包矩形拓宽二分之一分辨率大小作为空间插值区域范围；然后将空间插值区域范围依据分辨率划分为m行和n列，得到m*n个网格单元，作为待插值格网点的计算单元。

3.根据权利要求1所述的一种基于反向计算思维的空间插值方法，其特征在于，所述步骤S2中，计算各个待插值格网点所占插值权重的具体过程为：

S21：根据插值区域起始位置坐标、采样点的三元坐标以及插值分辨率，计算采样点所在的位置，即采样点的行号和列号；

S22：根据采样点的行号和列号计算采样点周围的待插值格网点的行号和列号；

S23：根据待插值格网点的行号和列号，结合插值区域起始点坐标及插值分辨率，计算待插值格网点的坐标；

S24：分别计算每个待插值格网点与采样点之间的距离；

S25：根据待插值格网点与采样点之间的距离，计算每个待插值格网点所占插值权重。

4.根据权利要求1所述的一种基于反向计算思维的空间插值方法，其特征在于，所述步骤S3中，计算各个待插值格网点附近插值格网点所占插值权重的具体过程为：

S31：根据划分的网格单元，逐行构建每个待插值格网点的约束方程；根据当前处理的待插值格网点的行号和列号，判断待插值格网点周围采样点是否满足推出待插值格网点的插值属性计算需求，如果满足则转入步骤S4，否则执行步骤S32；

S32：基于待插值格网点的位置，确定构建约束方程的周围格网点m′个数，取值为3、5和8；计算周围格网点到待插值格网点距离；

S33：根据待插值格网点与周围格网点之间的距离，计算每个周围点所占插值权重。

Images