CN106600537B - 一种反距离权重的异向性三维空间插值方法 - Google Patents

Abstract

本发明公开了一种反距离权重的异向性三维空间插值方法，包含以下步骤：输入三维空间采样数据；三维空间采样数据旋转变换；三维空间采样数据各向异性探索；三维空间拉伸变换；反距离权重的异向性三维空间插值；三维空间插值可视化。本发明能够在使用反距离权重插值时顾及三维地理空间现象的各向异性特征，更好的反映三维地理空间场的重建现象。

Description

一种反距离权重的异向性三维空间插值方法

技术领域

本发明涉及一种三维空间插值方法，特别是一种反距离权重的异向性三维空间插值方法。

背景技术

三维地理现象具有各向异性特征。如何利用有限的三维观测采样数据直接对具有各向异性三维地理现象进行可靠的三维空间插值不仅是三维空间分析的重要问题也是三维地理信息系统的基本功能需要。

反距离权重插值方法作为一种精确性插值方法，具有形式简单和不受维度限制等优点，可以用于三维空间空间插值。但是目前反距离权重插值方法是基于地理学第一定律的，默认插值点与多个参考点在距离相等的情况下，每个参考点对插值点的权重贡献度相同，但是实际上地理现象受多种因素的影响和制约，与插值点相同距离的参考点对插值点的权重是不同的，直接采用反距离权重插值方法插值会导致插值精度有所下降，因此需要在采用反距离权重插值时，分析三维空间数据的各向异性特征，根据各向异性特征进行三维反距离权重插值，以获取更高的插值精度和可靠的三维空间场重建。

发明内容

本发明所要解决的技术问题是提供一种反距离权重的异向性三维空间插值方法，它具有更高的插值精度。

为解决上述技术问题，本发明所采用的技术方案是：

一种反距离权重的异向性三维空间插值方法，其特征在于包含以下步骤：

步骤一：输入三维空间采样数据；

步骤二：三维空间采样数据旋转变换：

步骤三：三维空间采样数据各向异性探索；

步骤四：三维空间拉伸变换；

步骤五：反距离权重的异向性三维空间插值；

步骤六：三维空间插值可视化。

进一步地，所述步骤一中采样数据包含，

三维空间采样数据集记为：S＝{(P_i，f_i)，i＝1，2，3，…，n}，其中P_i表示第i个采样点的三维坐标(x_i，y_i，z_i)，f_i为第i个采样点的属性值，记S中所有三维坐标P_i的集合为P:

P中各坐标的分量记为：

记

S中所有属性值f_i的集合为f，表达式如下：

记μ_f为所有采样数据属性值的平均值，

记R为三维空间的旋转矩阵，L为三维空间的拉伸矩阵，

记S_R＝{(P_i′,f_i),i＝1,2,3,...,n}，其中P_i′表示第i个采样点通过旋转矩阵R变换后的三维坐标(x_i′,y_i′,z_i′)，其中属性值f_i不参与旋转和拉伸变换，

记S_R·L＝{(P_i″,f_i),i＝1,2,3,...,n}，其中P_i″表示第i个采样点通过旋转矩阵R和拉伸矩阵L变换后的三维坐标(x_i″,y_i″,z_i″)，其中f_i属性值不参与旋转和拉伸变换。

记三维空间插值点集为：S_{Interpolation}＝{(I_j，value_i)，j＝1，2，3，…，m}，其中I_j表示第j个插值点的三维坐标(u_j，v_j，w_j)，value_j为第j个插值点的属性值，I_j″表示第j个插值点通过旋转矩阵R和拉伸矩阵L变换后的三维坐标(u_j″,v_j″,w_j″)，其中value_j为插值点的属性值不参与旋转和拉伸变换。

进一步地，所述步骤二具体为，

根据平均主海森矩阵表达式

∑_fppb_j＝λ_j∑_pb_j

(1)

求取∑_fpp的广义特征值λ_j及其对应的特征向量b_j，其中∑_fpp为加权协方差矩阵，

其中n为采样点的总数，∑_p表示P的协方差矩阵，

通过平均主海森矩阵求取的三个特征值大小为λ₁>λ₂>λ₃，所对应的三个非正交向量为将最大特征值λ₁所对应的特征向量定义为第一主轴，将其单位化后记为第二大的特征值λ₂所对应的特征向量定义为第二主轴，将其单位化后记为 最小特征值λ₃所对应的特征向量定义为第三主轴，将其单位化后记为

保持第一主轴方向不变，将第二主轴进行投影正交变换使之与第一主轴正交，记为：

求取与组成平面的法向量，即为新的第三主轴通过投影正交变换之后，将第一主轴新的第二主轴新的第三主轴组成为坐标旋转矩阵

采用坐标旋转矩阵R对三维空间采样数据集S进行旋转变换，变换后的三维空间采样数据集记为S_R：

S_R＝S·R

(4)

建立三维空间采样点坐标P及其属性值f_i与旋转拉伸变换之后P_i′及其f_i之间的映射关系，即(x_i,y_i,z_i,f_i)与(x_i′,y_i′,z_i′,f_i)映射关系，即：(x_i,y_i,z_i,f_i)→(x_i′,y_i′,z_i′,f_i)，其中属性值f_i不变。

进一步地，所述的正交计算过程为：

A)求取和2个向量的叉积，

B)求得与的叉积，

进一步地，所述步骤三具体为，

以变换后的三维空间采样数据集S_R为对象，采用地统计学方法分别沿三维坐标轴计算三个方向上变异函数值，选择相应的变异函数拟合模型，一般可以选择球状模拟、指数模型、高斯模型、幂函数模型、对数函数模型等对各方向上的变异函数值进行拟合，分别求取拟合后的三个方向上的变程值，记为：a_x、a_y、a_z，其所对应的变异函数值为：γ(h_x)、γ(h_y)、γ(h_z)；

求取γ(h_x)、γ(h_y)、γ(h_z)三个数值的最小值，记为γ(h_min)；

将γ(h_min)设置为变异函数拟合模型公式的函数值，即变异函数拟合模型公式的因变量，分别求取x、y、z三个方向上变异函数拟合模型公式的自变量，记为h_x′、h_y′、h_z′。

进一步地，所述步骤四具体为，

求取h_x′、h_y′、h_z′中的最大值，记为h′_max，构建将各向异性空间转换为向各向同性空间的拉伸转换L：

采用拉伸矩阵L对三维空间采样数据集S′进行旋转变换，变换后的三维空间采样数据集记为：S_R·L，

S_R·L＝S_R·L

(5)

建立三维空间采样点坐标P及其属性值f_i与旋转拉伸变换之后P_i″及其f_i之间的映射关系，即(x_i,y_i,z_i,f_i)与(x_i″,y_i″,z_i″,f_i)的映射关系，(x_i,y_i,z_i,f_i)→(x_i″,y_i″,z_i″,f_i)，其中属性值f_i不变。

进一步地，所述步骤五具体为，

通过集合S中P的边界，构建所有三维空间数据的插值点；

对于任何一个三维空间插值点I_j，其三维坐标(u_j,v_j,w_j)，将其通过旋转矩阵R和拉伸矩阵L变换记为：I″_j，其三维坐标记为：(u_j″,v_j″,w_j″)，而插值点的属性值不参与旋转和拉伸，具体计算公式如下：

I″_j＝I_j·R·L

(6)

建立旋转和拉伸变换之前I_j和旋转和拉伸变换之后I″_j的映射关系，即(u_j,v_j,w_j,f_j)与(u_j″,v_j″,w_j″,f_j)的映射关系，即：(u_j,v_j,w_j,value_j)→(u_j″,v_j″,w_j″,value_j)，其中插值点属性值不变；

对于任一插值点(u_j,v_j,w_j)的属性值value_j，采用旋转及变换后的三维空间坐标(u_j″,v_j″,w_j″)和旋转及变换后的参考点三维空间坐标(x_i″,y_i″,z_i″)进行反距离权重插值方法计算，首先计算插值点三维空间坐标(u_j″,v_j″,w_j″)与所有三维空间采样数据三维空间坐标(x_i″,y_i″,z_i″)的欧式距离，取与(u_j″,v_j″,w_j″)距离最近的前k个三维空间采样数据的坐标及其属性值参与反距离权重计算，具体公式如下：

其中

进一步地，所述构建所有三维空间数据的插值点具体步骤为：

A)求取集合S中P的各分量的最大值和最小值，分别记为：x_max，x_min，y_max，y_min，z_max，z_min；

B)分别对各分量做差计算，三个方向上的分量差值记为：x_distance，y_distance，z_distance，求取三者之间的最小值min_distance；

C)计算构建插值点的间距interval＝min_distance/number，number的数值自行设定；

D)分别从x_min、y_min、z_min开始按照interval等间距构建插值点I_i，i＝1，2，3，...，m。

进一步地，所述步骤六具体为，根据插值点的属性值value_j大小进行颜色渲染，对颜色渲染后的结果进行多个截面显示。

本发明与现有技术相比，具有以下优点和效果：本发明对原始三维空间采样数据进行平均主海森矩阵计算，求取矩阵的广义特征向量，基于最小投影原理将广义非正交特征向量转换为标准正交特征向量，将标准正交特征向量组成三维坐标旋转矩阵，对原始三维空间采样数据进行旋转变化使其变换为以标准正交特征向量的x、y、z坐标系下的三维数据，分别沿x、y、z三个方向计算变异函数值并选择相应的模型进行拟合，求取x、y、z三个方向上变程及所对应的最小变异函数值，以最小变异函数为变异函数拟合模型函数的函数值，分别求取x、y、z三个方向上的自变量，并根据自变量比值大小，构建3*3的拉伸矩阵，将各向异性空间变换为各向同性空间，采用反距离权重进行三维空间插值计算。通过这样的方法，能够在使用反距离权重插值时顾及三维地理空间现象的各向异性特征，更好的反映三维地理空间场的重建现象。

附图说明

图1是本发明的一种反距离权重的异向性三维空间插值方法的流程图。

图2是本发明的实施例采用的矿石品位数据图。

图3是本发明的x方向上变异函数的球状模型拟合曲线。

图4是本发明的y方向上变异函数的球状模型拟合曲线。

图5是本发明的z方向上变异函数的球状模型拟合曲线。

图6是本发明的最终插值结果可视化效果。

具体实施方式

下面结合附图并通过实施例对本发明作进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于以下实施例。

符号说明：

将三维空间采样数据集记为：S＝{(P_i，f_i)，i＝1，2，3，…，n}，其中P_i表示第i个采样点的三维坐标(x_i，y_i，z_i)，f_i为第i个采样点的属性值，记S中所有三维坐标P_i的集合为P，表达式如下所示：

P中各坐标的分量记为：

记

S中所有属性值f_i的集合为f，表达式如下：

记μ_f为所有采样数据属性值的平均值，

记R为三维空间的旋转矩阵，L为三维空间的拉伸矩阵，

记S_R＝{(P_i′,f_i),i＝1,2,3,...,n}，其中P_i′表示第i个采样点通过旋转矩阵R变换后的三维坐标(x_i′,y_i′,z_i′)，其中属性值f_i不参与旋转和拉伸变换；

记S_R·L＝{(P_i″,f_i),i＝1,2,3,...,n}，其中P_i″表示第i个采样点通过旋转矩阵R和拉伸矩阵L变换后的三维坐标(x_i″,y_i″,z_i″)，其中f_i属性值不参与旋转和拉伸变换。

记三维空间插值点集为：S_{Interpolation}＝{(I_j，value_i)，j＝1，2，3，…，m}，其中I_j表示第j个插值点的三维坐标(u_j，v_j，w_j)，value_j为第j个插值点的属性值，需要通过插值方法计算的属性值。I″_j表示第j个插值点通过旋转矩阵R和拉伸矩阵L变换后的三维坐标(u_j″,v_j″,w_j″)，其中value_j为插值点的属性值不参与旋转和拉伸变换。

如图1所示，本发明的一种反距离权重的异向性三维空间插值方法，包含以下步骤：

(1)根据平均主海森矩阵表达式，如下(1)式为计算表达式，求取∑_fpp的广义特征值λ_j及其对应的特征向量b_j，

∑_fppb_j＝λ_j∑_pb_j

(1)

其中∑_fpp为加权协方差矩阵，具体计算表达式如下2所示

其中n为采样点的总数，∑_p表示P的协方差矩阵，表达式如下3所示，

通过平均主海森矩阵求取的三个特征值大小为λ₁>λ₂>λ₃，所对应的三个非正交向量为将最大特征值λ₁所对应的特征向量定义为第一主轴，将其单位化后记为第二大的特征值λ₂所对应的特征向量定义为第二主轴，将其单位化后记为最小特征值λ₃所对应的特征向量定义为第三主轴，将其单位化后记为

(2)保持第一主轴方向不变，将第二主轴进行投影正交变换使之与第一主轴正交，记为：

具体步骤如下：

1.求取和2个向量的叉积，

2.求得与的叉积，

这样是与原始第二轴偏移量最小且与第一主轴正交的向量。

(3)求取与组成平面的法向量，即为新的第三主轴，记变换后新的第三主轴为通过投影正交变换之后，将第一主轴新的第二主轴新的第三主轴组成为坐标旋转矩阵

(4)采用坐标旋转矩阵R对三维空间采样数据集S进行旋转变换，变换后的三维空间采样数据集记为S_R，具体计算公式如下：

S_R＝S·R

(4)

建立三维空间采样点坐标P及其属性值f_i与旋转拉伸变换之后P_i′及其f_i之间的映射关系，即为(x_i,y_i,z_i,f_i)与(x_i′,y_i′,z_i′,f_i)的映射关系：(x_i,y_i,z_i,f_i)→(x_i′,y_i′,z_i′,f_i)，其中属性值f_i不变。

(5)以变换后的三维空间采样数据集S_R为对象，采用地统计学方法分别沿三维坐标轴计算三个方向上变异函数值，选择相应的变异函数拟合模型，一般可以选择球状模拟、指数模型、高斯模型、幂函数模型、对数函数模型等，对各方向上的变异函数值进行拟合，分别求取拟合后的三个方向上的变程值，记为：a_x、a_y、a_z，其所对应的变异函数值为：γ(h_x)、γ(h_y)、γ(h_z)；

(6)求取γ(h_x)、γ(h_y)、γ(h_z)三个数值的最小值，记为γ(h_min)；

(7)将γ(h_min)设置为所步骤(5)变异函数拟合模型公式的函数值，即变异函数拟合模型公式的因变量，分别求取x、y、z三个方向上变异函数拟合模型公式的自变量，记为h_x′、h_y′、h_z′；

(8)求取h_x′、h_y′、h_z′中的最大值，记为h′_max，构建将各向异性空间转换为各向同性空间的拉伸转换矩阵L，表达式如下：

(9)采用拉伸矩阵L对三维空间采样数据集S′进行旋转变换，变换后的三维空间采样数据集记为：S_R·L，具体计算公式为：

S_R·L＝S_R·L

(5)

建立三维空间采样点坐标P及其属性值f_i与旋转拉伸变换之后P_i″及其f_i之间的映射关系，即(x_i,y_i,z_i,f_i)与(x_i″,y_i″,z_i″,f_i)的映射关系，(x_i,y_i,z_i,f_i)→(x_i″,y_i″,z_i″,f_i)，其中属性值f_i不变。

(10)通过集合S中P的边界，构建所有三维空间数据的插值点，具体步骤如下：

1：求取集合S中P的各分量的最大值和最小值，分别记为：x_max，x_min，y_max，y_min，z_max，z_min；

2：分别对各分量做差值计算，三个方向上的分量差值记为：x_distance，y_distance，z_distance，求取三者之间的最小值min_distance；

3：计算构建插值点的间距interval＝min_distance/number，number的值可以自行设定，一般不超过250；

4：分别从x_min、y_min、z_min开始按照interval等间距构建插值点I_i，i＝1,2,3，...，m。

(11)对于任何一个三维空间插值点I_j，其三维坐标(u_j,v_j,w_j)，将其通过旋转矩阵R和拉伸矩阵L变换记为：I″_j，其三维坐标记为：(u_j″,v_j″,w_j″)，而插值点的属性值不参与旋转和拉伸，具体计算公式如下：

I″_j＝I_j·R·L

(6)

建立旋转和拉伸变换之前I_j和旋转和拉伸变换之后I″_j的映射关系，即(u_j,v_j,w_j,f_j)与(u_j″,v_j″,w_j″,f_j)的映射关系，即：(u_j,v_j,w_j,value_j)→(u_j″,v_j″,w_j″,value_j)，其中插值点属性值不变。

(12)对于任一插值点(u_j,v_j,w_j)的属性值value_j，采用旋转及变换后的三维空间坐标(u_j″,v_j″,w_j″)和旋转及变换后的参考点三维空间坐标(x_i″,y_i″,z_i″)进行反距离权重插值方法计算，首先计算插值点三维空间坐标(u_j″,v_j″,w_j″)与所有三维空间采样数据三维空间坐标(x_i″,y_i″,z_i″)的欧式距离，取与(u_j″,v_j″,w_j″)距离最近的前k个三维空间采样数据的坐标及其属性值参与反距离权重计算，具体公式如下：

其中

(13)根据插值点的属性值value_j大小进行颜色渲染，对颜色渲染后的结果进行多个截面显示。

实施例的三维空间采样数据选择4125个矿石品位数据，根据说明书发明内容的符号，则S＝{(P_i，f_i)，i＝1，2，3，…，4125}，铁品位的范围为(0，100)，如图2所示。

(1)矿石品位数据的导入

采用Matlab软件导入4125个矿石品位数据，第i个矿石品位数据的三维坐标信息(x_i，y_i，z_i)和属性信息f_i，数据导入后，以前5个矿石品位为例说明，数据具体如下：

(35389，19335.5，-53.83，0.4310)

(35389，19335.5，-57.24，0.3188)

(35389，19335.5，-71.37，0.5875)

(35389，19335.5，-76.26，0.5855)

(35389，19335.5，-79.54，0.4894)

每行数据前三列代表某个矿石品位数据的三维空间数据坐标，是需要后续旋转和拉伸变换的，最后一列表示属性值f_i，本实例为矿石品位，不参与旋转和拉伸变换。

(2)三维空间采样数据的旋转变换

步骤一：根据平均主海森矩阵表达式公式(1)可计算出4125个矿石品位数据的三个广义正交特征向量和特征值，对特征值大小进行排序，最大特征值所对应的特征向量记为：第二大特征值所对应的特征向量记为：最小特征值所对应的特征向量记为

步骤二：分别对进行单位化，单位化后分别记为： 具体计算公式为：计算后各值为：

步骤三：将单位化后的进行投影正交变换，使其三个向量两两正交，首先求取和2个向量的叉积，记为然后求取与的叉积，记为：最后求取与组成平面的法向量，记为：通过投影正交变换之后，为相互正交的向量，将其组成为坐标旋转矩阵具体计算后R值为：

步骤四：采用坐标旋转矩阵R对4125个矿石品位数据S进行旋转变换，变换后的4125个矿石品位数据记为S_R，具体计算公式为：S_R＝S·R，前5个数据旋转之后的三维空间坐标及其属性值具体信息如下：

(-27088，3897，29521，0.4310)

(-27090，3897，29519，0.3188)

(-27101，3897，29510，0.5875)

(-27104，3897，29506，0.5855)

(-27107，3897，29504，0.4894)

(3)三维空间采样数据的拉伸变换

步骤一：以变换后4125个矿石品位数据S_R为对象，采用地统计学方法分别沿三维坐标轴计算三个方向上变异函数值，选择相应的变异函数拟合模型，以球状模型为例，计算公式为：

式中c₀表示块金，a表示变程，c₀+c₁表示基台值，h表示计算变异函数的分组距离，本实例设置为20，分别设置各个方向的角度分组为45°，γ(h)为相距h的点对变异函数值，选取球状模型对各方向上的变异函数值进行拟合，分别求取拟合后的三个方向上的变程值，分别记为：a_x、a_y、a_z，具体值为：a_x＝543.87，a_y＝441.46，a_z＝524.31，其所对应的变异函数值为：γ(h_x)＝0.0357、γ(h_y)＝0.0303、γ(h_z)＝0.053；

步骤二：求取γ(h_x)、γ(h_y)、γ(h_z)三个数值的最小值，记为γ(h_min)＝0.0303；

步骤三：将γ(h_min)＝0.0303设置为步骤一中变异函数拟合模型公式的函数值，即变异函数拟合模型公式的因变量，分别求取x、y、z三个方向上变异函数拟合模型公式的自变量，记为h_x′、h_y′、h_z′，其值分别为：308.35、404.46、218.79

步骤四：求取h_x′、h_y′、h_z′中的最大值，记为h′_max＝404.46，根据h_x′、h_y′、h_z′、h′_max构建将各向异性空间转换为各向同性空间的拉伸转换矩阵L，具体值为：

步骤五：采用拉伸矩阵L对三维空间采样数据集S′进行旋转变换，变换后的三维空间采样数据集记为：S_R·L，具体计算为公式(5)，可使三维空间数据从各向异性变换为各向同性，对前5个数据的三维空间坐标进行旋转和拉伸变换后的具体信息如下：

(-20651，3897，15969，0.4310)

(-20653，3897，15968，0.3188)

(-20661，3897，15963，0.5875)

(-20663，3897，15961，0.5855)

(-20666，3897，15960，0.4894)

步骤六：建立旋转和拉伸后的S_R·L三维空间采用数据坐标及其属性值与原始三维采用数据坐标及其属性值的映射关系，其属性值不变化，以前5个数据为例具体如下：

(35389，19335.5，-53.83，0.4310)→(-20653，3897，15968，0.4310)

(35389，19335.5，-57.24，0.3188)→(-20651，3897，15969，0.3188)

(35389，19335.5，-71.37，0.5875)→(-20661，3897，15963，0.5875)

(35389，19335.5，-76.26，0.5855)→(-20663，3897，15961，0.5855)

(35389，19335.5，-79.54，0.4894)→(-20666，3897，15960，0.4894)

(4)通过集合S中P的边界，构建所有三维空间数据的插值点，具体步骤如下：

1：求取集合S中P的各分量的最大值和最小值，各分量的最大最小值为：x_max＝36316.45，x_min＝35119.78，y_max＝19822.76，y_min＝19052.14，z_max＝-0.23，z_min＝-665.89；

2：分别对各分量做差计算，三个方向上的分量差值为：x_distance＝1196.67，y_distance＝770.62，z_distance＝665.66；

3：计算构建插值点的间距x_interval＝x_distance/number，y_interval＝y_distance/number,z_interval＝z_distance/number,number的值可以自行设定，一般不超过250，本实例采用250；

4：从x_min、y_min、z_min开始，按照x_interval＝4.79、y_interval＝3.08、z_interval＝2.66等间距构建插值点I_j，j＝1,2,3，...m。

(5)三维空间插值点坐标旋转与拉伸变换

对于任何一个三维空间插值点I_j，其三维坐标(u_j,v_j,w_j)，将其通过旋转矩阵R和拉伸矩阵L变换记为：I″_j，其三维坐标记为：(u_j″,v_j″,w_j″)，而插值点的属性值不参与旋转和拉伸，具体计算公式如下：

I″_j＝I_j·R·L

建立旋转和拉伸变换之前I_j和旋转和拉伸变换之后I″_j的映射关系，即(u_j,v_j,w_j,value_j)与(u_j″,v_j″,w_j″,value_j)的映射关系，即：(u_j,v_j,w_j,value_j)→(u_j″,v_j″,w_j″,value_j)，其中插值点属性值不变，取j＝10000插值点为例，映射关系如下：

(35306.46，19699.46，-665.89，value₁₀₀₀₀)→(-21026.66，426.50，15773.75，value₁₀₀₀₀)

(6)反距离权重的异向性三维空间插值方法

对于任一插值点(u_j,v_j,w_j)的属性值value_j，采用旋转及变换后的三维空间坐标(u_j″,v_j″,w_j″)和旋转及变换后的参考点三维空间坐标(x_i″,y_i″,z_i″)进行反距离权重插值方法计算，首先计算插值点三维空间坐标(u_j″,v_j″,w_j″)与所有三维空间采样数据三维空间坐标(x_i″,y_i″,z_i″)的欧式距离，取与(u_j″,v_j″,w_j″)距离最近的前k个三维空间采样数据的坐标及其属性值参与反距离权重计算，具体公式如下：

其中

j＝10000插值点的属性值采用公式(7)计算为0.64。

(7)根据插值点的属性值value_j大小进行颜色渲染，对颜色渲染后的结果进行多个截面显示，由于建立了插值点旋转和拉伸前后的映射关系，绘制图形时只需采用插值点三维空间坐标(u_j,v_j,w_j)及其对应的属性值value_j进行渲染即可。对插值数据进行可视化，在Matlab中采用slice函数进行横截面三维显示，其结果见图6所示。

对于反距离权重的异向性三维空间插值方法精度的验证采用逐点交叉验证，可与普通反距离权重的三维空间插值方法进行对比，采用平均误差和均方根误差进行对比分析，以4125个矿石品位数据为例，本发明专利的平均误差：6.63，均方根误差为，9.99而普通反距离权重的三维空间插值方法的平均误差为：6.89均方根误差为：10.76，说明本发明专利具有较好的插值精度。

本说明书中所描述的以上内容仅仅是对本发明所作的举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种修改或补充或采用类似的方式替代，只要不偏离本发明说明书的内容或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。

Claims (5)

1.一种反距离权重的异向性三维空间插值方法，其特征在于包含以下步骤：

步骤一：输入三维空间采样数据；

采样数据包含，

三维空间采样数据集记为：S＝{(P_i，f_i)，i＝1，2，3，…，n}，其中P_i表示第i个采样点的三维坐标(x_i，y_i，z_i)，f_i为第i个采样点的属性值，记S中所有三维坐标P_i的集合为P:

P中各坐标的分量记为：

记

S中所有属性值f_i的集合为f，表达式如下：

记μ_f为所有采样数据属性值的平均值，

记R为三维空间的旋转矩阵，L为三维空间的拉伸矩阵，

记S_R＝{(P′_i,f_i),i＝1,2,3,...,n}，其中P′_i表示第i个采样点通过旋转矩阵R变换后的三维坐标(x′_i,y′_i,z′_i)，其中属性值f_i不参与旋转和拉伸变换，

记S_R·L＝{(P″_i,f_i),i＝1,2,3,...,n}，其中P″_i表示第i个采样点通过旋转矩阵R和拉伸矩阵L变换后的三维坐标(x″_i,y″_i,z″_i)，其中f_i属性值不参与旋转和拉伸变换；

记S_{Interpolation}＝{(I_j，value_i)，j＝1，2，3，…，m}为三维空间插值点集，其中I_j表示第j个插值点的三维坐标(u_j，v_j，w_j)，value_j为第j个插值点的属性值，I″_j表示第j个插值点通过旋转矩阵R和拉伸矩阵L变换后的三维坐标(u″_j,v″_j,w″_j)，其中value_j为插值点的属性值，不参与旋转和拉伸变换；

步骤二：三维空间采样数据旋转变换；

步骤三：三维空间采样数据各向异性探索；

以变换后的三维空间采样数据集S_R为对象，采用地统计学方法分别沿三维坐标轴计算三个方向上变异函数值，选择相应的变异函数拟合模型，可选择的模型有球状模型、指数模型、高斯模型、幂函数模型和对数函数模型，对各方向上的变异函数值进行拟合，分别求取拟合后的三个方向上的变程值，记为：a_x、a_y、a_z，其所对应的变异函数值为：γ(h_x)、γ(h_y)、γ(h_z)；

求取γ(h_x)、γ(h_y)、γ(h_z)三个数值的最小值，记为γ(h_min)；

将γ(h_min)设置为变异函数拟合模型公式的函数值，即变异函数拟合模型公式的因变量，分别求取x、y、z三个方向上变异函数拟合模型公式的自变量，记为h′_x、h′_y、h′_z；

步骤四：三维空间拉伸变换；

求取h′_x、h′_y、h′_z中的最大值，记为h′_max，构建将各向异性空间转换为各向同性空间的拉伸转换L：

采用拉伸矩阵L对三维空间采样数据集S′进行旋转变换，变换后的三维空间采样数据集记为：S_R·L，

S_R·L＝S_R·L (5)

建立三维空间采样点坐标P及其属性值f_i与旋转拉伸变换之后P″_i及其f_i之间的映射关系，即(x_i,y_i,z_i,f_i)与(x″_i,y″_i,z″_i,f″_i)的映射关系，(x_i,y_i,z_i,f_i)→(x″_i,y″_i,z″_i,f″_i)，其中属性值f_i不变；

步骤五：反距离权重的异向性三维空间插值；

通过集合S中P的边界，构建所有三维空间数据的插值点；

对于任何一个三维空间插值点I_j，其三维坐标(u_j,v_j,w_j)，将其通过旋转矩阵R和拉伸矩阵L变换记为：I″_j，其三维坐标记为：(u″_j,v″_j,w″_j)，而插值点的属性值不参与旋转和拉伸，具体计算公式如下：

I″_j＝I_j·R·L (6)

建立旋转和拉伸变换之前I_j和旋转和拉伸变换之后I″_j的映射关系，即(u_j,v_j,w_j,f_j)与(u″_j,v″_j,w″_j,f″_j)的映射关系，即：(u_j,v_j,w_j,value_j)→(u″_j,v″_j,w″_j,value_j)，其中插值点属性值values_j不参与旋转和变换；

对于任一插值点(u_j,v_j,w_j)的属性值value_j，采用旋转及变换后的三维空间坐标(u″_j,v″_j,w″_j)和旋转及变换后的参考点三维空间坐标(x″_i,y″_i,z″_i)进行反距离权重插值方法计算，首先计算插值点三维空间坐标(u″_j,v″_j,w″_j)与所有三维空间采样数据三维空间坐标(x″_i,y″_i,z″_i)的欧式距离，取与(u″_j,v″_j,w″_j)距离最近的前k个三维空间采样数据的坐标及其属性值参与反距离权重计算，具体公式如下：

其中

步骤六：三维空间插值可视化。

2.按照权利要求1所述的一种反距离权重的异向性三维空间插值方法，其特征在于：所述步骤二具体为，

根据平均主海森矩阵表达式

∑_fppb_j＝λ_j∑_pb_j (1)

求取∑_fpp的广义特征值λ_j及其对应的特征向量b_j，其中∑_fpp为加权协方差矩阵，

其中n为采样点的总数，∑_p表示P的协方差矩阵，

通过平均主海森矩阵求取的三个特征值大小为λ₁>λ₂>λ₃，所对应的三个非正交向量为将最大特征值λ₁所对应的特征向量定义为第一主轴，将其单位化后记为第二大的特征值λ₂所对应的特征向量定义为第二主轴，将其单位化后记为 最小特征值λ₃所对应的特征向量定义为第三主轴，将其单位化后记为

保持第一主轴方向不变，将第二主轴进行投影正交变换使之与第一主轴正交，记为：

求取与组成平面的法向量，即为新的第三主轴 通过投影正交变换之后，将第一主轴新的第二主轴新的第三主轴组成为坐标旋转矩阵

采用坐标旋转矩阵R对三维空间采样数据集S进行旋转变换，变换后的三维空间采样数据集记为S_R：

S_R＝S·R (4)

建立三维空间采样点坐标P_i及其属性值f_i与旋转变换之后P_i′及其f_i之间的映射关系，即为(x_i,y_i,z_i,f_i)与(x′_i,y′_i,z′_i,f_i)的映射关系：(x_i,y_i,z_i,f_i)→(x′_i,y′_i,z′_i,f_i)，其中属性值f_i不变。

3.按照权利要求2所述的一种反距离权重的异向性三维空间插值方法，其特征在于：所述的正交计算过程为：

A)求取和2个向量的叉积，

B)求得与的叉积，

4.按照权利要求1所述的一种反距离权重的异向性三维空间插值方法，其特征在于：所述构建所有三维空间数据的插值点具体步骤为：

A)求取集合S中P的各分量的最大值和最小值，分别记为：x_max，x_min，y_max，y_min，z_max，z_min；

B)分别对各分量做差计算，三个方向上的分量差值记为：x_distance，y_distance，z_distance，求取三者之间的最小值min_distance；

C)计算构建插值点的间距interval＝min_distance/number，number的值可自行设定；

D)分别从x_min、y_min、z_min开始按照interval等间距构建插值点I_i，i＝1，2，3，...，m。

5.按照权利要求1所述的一种反距离权重的异向性三维空间插值方法，其特征在于：所述步骤六具体为，根据插值点的属性值大小进行颜色渲染，对颜色渲染后的结果进行多个截面显示。