一种基于地质曲面局部复杂度的自适应混合插值方法与流程

本发明属于地质曲面重构技术领域，具体涉及一种基于地质曲面局部复杂度的自适应混合插值方法的设计。

背景技术：

在工程地质领域中，从野外获取地表地形、地层界限、断层、地下水位和风化层厚度分布等各种原始地质数据大多数并不是沿规则网格分布的，而是离散不连续的数据，并且规模比较庞大，不能满足三维建模的要求，因此需要通过选择合适的空间插值技术对地质信息进行处理，重构地质空间界面。

基于多项式曲面拟合的方法是解决复杂地质曲面重构的最有效的方法之一。用多项式拟合方法重构曲面的困难在于如何确定多项式的阶，如果我们用低阶多项式拟合，在陡峭的地区将有一个较大的误差，但如果我们使用高阶多项式拟合，重构曲面将产生伪地质结构。

混合插值的目的是结合低阶和高阶多项式插值方法的优点，通过分析地质的复杂性，可以平衡高阶和低阶的权重。现有的方法通常选取一个全局的权重值来确定高阶和低阶的比值。事实上，实际的地质曲面可能会在某些地区有更多的高阶成分，而在其他地区有更多的低阶成分。因此，全局的权重值不能适应实际的地质表面重构的要求。

W.Z.SHI采用了一种新型的权重值来结合双线性和双立方插值方法的混合插值方法。不同的地质曲面的复杂度是不同的，但这种方法可以通过控制在每个曲面上的权重值满足不同类型的曲面。它比只采用线性插值法或非线性插值法更精确，即使这个权重值依赖于被内插的地形的复杂性。

然而在同一个地质曲面的不同区域的复杂度也是不同的，在一个曲面上使用一个权重值依旧会造成较大的误差，因此需要在目前研究现状的基础上提出一种自适应检测地质曲面的局部复杂度来匹配动态的权重值的方法。

技术实现要素：

本发明的目的是为了解决现有技术中基于权重值的混合多项式曲面拟合方法不能很好应用于局部复杂度变化太大的地质曲面的问题，提出了一种基于地质曲面局部复杂度的自适应混合插值方法。

本发明的技术方案为：一种基于地质曲面局部复杂度的自适应混合插值方法，包括以下步骤：

S1、构建双线性插值函数作为混合插值的低阶成分；

S2、构建双立方插值函数作为混合插值的高阶成分；

S3、根据双线性插值函数和双立方插值函数构建混合插值函数；

S4、在每一个插值区域中求取混合插值函数的高低阶比例值ρ。

进一步地，步骤S1中双线性插值函数表示为：

f₁(x,y)＝a₁+a₂x+a₃y+a₄xy (1)。

进一步地，步骤S2中双立方插值函数表示为：

f₂(x,y)＝b₁+b₂x+b₃y+b₄x²+b₅xy+b₆y²+b₇x³+b₈x²y+b₉xy²+b₁₀y³ (2)。

进一步地，步骤S3具体为：

构建混合插值函数：

I＝ρA+(1-ρ)B (3)

式中A＝f₁(x,y)，表示混合插值的低阶成分；B＝f₂(x,y)，表示混合插值的高阶成分；ρ表示混合插值函数的高低阶比例值，取值范围为0≤ρ≤1；令I＝f^ρ(x,y)，则有：

其中

进一步地，步骤S4具体为：

定义RMSE为：

其中f和f^ρ分别表示原始数据集和混合插值后的结果，f_i,j和分别是在空间坐标(i,j)上的数据值和插值结果；

将公式(3)变换为：

I＝ρ(A-B)+B (6)

也即是：

其中a_i,i＝1,2,...,n；b_j,j＝1,2,...,m；x,y都是已知的常数；然后在每一个插值区域得到了一个唯一的高低阶比例值ρ，以使得每一个构造曲面的RMSE值是最小的，具体过程为：

已知待插值利用n个已知点进行插值，则通过最小二乘法拟合曲面方程，整个插值的总误差为：

因为A和B中所有系数为常数，则令常数(A-B)为A_i，B为B_i，对系数ρ求偏导值：

当时，Q值最小，则插值区域的RMSE值最小，因此令公式(9)等于0得到一个一元一次方程，给定A_i和B_i值后可以求得ρ值并使得RMSE值最小。

本发明的有益效果是：本发明在进行地质曲面构造插值时，增加一步对于多项式方程高低阶比例控制系数ρ的计算，通过局部动态自适应的方式，针对每一个插值点计算该点的ρ值，使得在同一地质曲面上的各个局部区域能够自适应更加精确的比例值来匹配当前区域的地质复杂度。对于地质复杂度变化较大的地质曲面，能够很好的减少插值带来的误差，从而较大幅度的提高插值精度。

附图说明

图1为本发明提供的一种基于地质曲面局部复杂度的自适应混合插值方法流程图。

图2为本发明实施例的设定全局高低阶比例值ρ＝0.3构建的插值曲面示意图。

图3为本发明实施例的设定全局高低阶比例值ρ＝0.3的曲面等势线分布图。

图4为本发明实施例的设定全局高低阶比例值ρ＝0.7构建的插值曲面示意图。

图5为本发明实施例的设定全局高低阶比例值ρ＝0.7的曲面等势线分布图。

图6为本发明实施例的采用基于局部复杂度自适应的混合插值方法构建的插值曲面示意图。

图7为本发明实施例的采用基于局部复杂度自适应的混合插值方法的曲面等势线分布图。

具体实施方式

下面结合附图对本发明的实施例作进一步的说明。

一个多项式拟合公式的一般描述如下：一组样本数据点，一个该曲面可以近似估计插值点的值的函数。插值曲面通常表达如下：

其中，a_p,q表示多项式的系数，可以由上述给定的样本数据点求得。

与单一的采用线性或非线性插值相比，混合多项式插值方式可以更好的针对地质曲面的特征调整高阶与低阶比例值ρ。

基于此，本发明提供了一种基于地质曲面局部复杂度的自适应混合插值方法，如图1所示，包括以下步骤：

S1、构建双线性插值函数作为混合插值的低阶成分。

一个地质曲面包括两个基本成分，一个是海拔变化频率较低的区域，另一个是海拔变化频率较高的区域。因此，将线性和非线性插值函数结合起来重建曲面是合理的，双线性方法是最流行的线性插值方法，其函数表达式如下：

f₁(x,y)＝a₁+a₂x+a₃y+a₄xy (1)

S2、构建双立方插值函数作为混合插值的高阶成分。

双立方插值是一种非线性插值方法，比大多数方法更有效，其函数表达式如下：

f₂(x,y)＝b₁+b₂x+b₃y+b₄x²+b₅xy+b₆y²+b₇x³+b₈x²y+b₉xy²+b₁₀y³ (2)

S3、根据双线性插值函数和双立方插值函数构建混合插值函数。

由于在现实世界中，地质曲面包含低频率和高频率，显然使用混合插值方法比双线性或双立方的其中一个更好。混合插值函数表示如下：

I＝ρA+(1-ρ)B (3)

式中A＝f₁(x,y)，表示混合插值的低阶成分；B＝f₂(x,y)，表示混合插值的高阶成分；ρ表示混合插值函数的高低阶比例值，取值范围为0≤ρ≤1，当ρ值取0或1时，混合插值方法将成为双立方或双线性插值法。

令I＝f^ρ(x,y)，则有：

其中

S4、在每一个插值区域中求取混合插值函数的高低阶比例值ρ。

现有技术中，高低阶比例值ρ是一个需要事先给定的恒定的比率。但是即使在同一地质面上，不同的区域也存在着不同的复杂性。为了得到更精确的插值结果，我们应该在每个插补过程给高低阶比例值ρ赋予不同的值。

本发明实施例中，采用最小二乘法计算高低阶比例值ρ。如果给ρ指定一个值，我们可以得到混合插值系数和地质曲面的估计值，但这一结果不够准确。如果我们能找到一种自适应率保证在每个插值区域的均方根误差(RMSE)值最小，结果将会更准确。均方根误差(RMSE)的定义为：

其中f和f^ρ分别表示原始数据集和混合插值后的结果，f_i,j和分别是在空间坐标(i,j)上的数据值和插值结果。

因此我们把ρ作为系数来求解，并采用最小二乘法再求解系数ρ，将公式(3)变换为：

I＝ρ(A-B)+B (6)

也即是：

其中a_i,i＝1,2,...,n；b_j,j＝1,2,...,m；x,y都是已知的常数。现在,我们需要在每一个插值区域得到了一个唯一的高低阶比例值ρ，以使得每一个构造曲面的RMSE值是最小的，具体过程为：

已知待插值利用n个已知点进行插值，则通过最小二乘法拟合曲面方程。整个插值的总误差为：

因为A和B中所有系数为常数，则令常数(A-B)为A_i，B为B_i，对系数ρ求偏导值：

当时，Q值最小，则插值区域的RMSE值最小，因此令公式(9)等于0得到一个一元一次方程，给定A_i和B_i值后可以求得ρ值并使得RMSE值最小。

通过不同的方法重构地质曲面，可以看出各种方法的优缺点。如图2、图4、图6所示分别为设定全局高低阶比例值ρ＝0.3，ρ＝0.7以及采用本发明提供的基于局部复杂度自适应的混合插值方法所得到的地质曲面图。图中灰色的曲面是我们构建的插值曲面，黑色的曲线被称为种子线是由已知的地质曲面的数据集构造。浅色圈代表的伪地质构造地形，深色圈反映与给定的数据集的一致程度。为区分不同的插值方法，我们同时形成了重建地质曲面的等势线一同进行比较。图3、图5、图7分别为图2、图4、图6的等势线分布图。

由图2和图3可知，当设定全局高低阶比例值ρ＝0.3时，意味着它有着更多的双立方成分，所以它比双线性插值法更准确但是也构建了一些伪地质结构。由图4和图5可知，当设定全局高低阶比例值ρ＝0.7时，意味着它具有更多的双线性成分，插值结果比ρ＝0.3时更光滑，可以从图5的等势线中可以看出；但由于高频成分较少，插值精度不太高。由图6和图7可知，当采用采用本发明提供的基于局部复杂度自适应的混合插值方法时，比全局固定ρ值的插值结果更准确。同时，我们可以从图7看到，该方法的等势线十分有序，很好的结合了双线性和双立方方法的优点。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。