CN110097637A - 一种三维地质属性模型时空插值方法及系统 - Google Patents

Abstract

本发明公开了一种三维地质属性模型时空插值方法及系统。该方法包括：依据地形数据、钻孔数据和地质剖面数据，生成三维地质结构模型；对待插值区域进行分块，得到多个三维地质分块；对各三维地质分块对应的地质采集点属性数据依次进行水平空间插值和垂直空间插值，得到三维空间属性数据；构建完整时间序列，并采用求和自回归移动平均模型对完整时间序列进行建模分析和插补，得到时间属性数据；采用协方差模型将三维空间属性数据和对应的时间属性数据进行融合，得到三维地质属性模型。本发明能够提高插值精度，保证插值后时空数据的完整性。

Description

一种三维地质属性模型时空插值方法及系统

技术领域

本发明涉及地质勘探技术领域，特别是涉及一种三维地质属性模型时空插值方法及系统。

背景技术

随着经济的高速发展，快速的城市化进程加剧了城市中土地资源的供需矛盾，为解决“城市病”带来的交通拥堵、生态环境恶化等问题，迫切需要合理开发利用地下空间资源。地下空间开发离不开对地质条件的充分了解和对地下空间数据的定量化计算。三维地质建模通过真三维空间中对地质实体几何形态和地质实体内部属性参数建模，能直观了解地质体空间分布和变化规律。根据建模对象和建模方法的不同，三维地质建模分为三维地质结构建模和三维地质属性建模。其中，三维地质属性建模通过对地质体属性参数空间位置和动态变化的模拟，能直观展示地质体内部属性参数的空间分布、变化规律及参量间的关联关系，是地下空间数据定量化分析的重要手段。由于建模需要的采样信息大多由钻孔获得，常常数量有限，且分布不均，需要合理的插值方法来进行地质参数模拟。

目前，国际地学领域常用的插值方法主要包括反距离加权法、克里金插值法和序贯高斯模拟。其中，反距离加权法主要依赖于反距离的幂值，幂参数可基于距输出点的距离来控制已知点对内插值的影响。通过定义更高的幂值，可进一步强调最近点，邻近数据将受到更大影响，表面会变得更不平滑。指定较小的幂值将对距离较远的周围点产生更大的影响，从而导致平面更加平滑。克里金插值法首先考虑空间属性在空间位置上的变异分布，确定对一个待插点有影响的距离范围，用此范围内的采样点来估计待插点的属性值。根据样品空间位置不同、样品间相关程度的不同，对每个样品点赋予不同的权，进行加权平均来估计待插点的属性值。序贯高斯模拟在综合考虑了样品的大小、形状、待插值区域相互之间的空间分布以及区域变量的空间结构信息后，为了反映变量空间波动性，沿着随机路径序贯求取各节点的条件累计分布函数，并从条件累计分布函数中提取模拟值。

常用的空间插值方法多针对二维、三维空间，未考虑地下空间数据的时空属性和地质空间的异质性。地下空间数据不仅具有三维空间属性，还具有时间属性，是研究数据动态变化规律和发展趋势的基础。时间维度的丢失将严重破坏数据的完整性。由于不同专业领域数据的采样周期各异，数据的时间属性不连续，需要进行插值计算。不同地质条件的区域之间的属性值存在较大差异，如断层两侧邻近两个点的属性可能差异巨大，导致待插点的属性估计错误，严重影响插值精度。

发明内容

基于此，有必要提供一种三维地质属性模型时空插值方法及系统，以提高插值精度，保证插值后时空数据的完整性。

为实现上述目的，本发明提供了如下方案：

一种三维地质属性模型时空插值方法，包括：

获取待插值区域的地形数据、钻孔数据、地质剖面数据和地质采集点属性数据；

依据所述地形数据、所述钻孔数据和所述地质剖面数据，生成三维地质结构模型；所述三维地质结构模型表示所述待插值区域内断层、地层和岩体三者之间的接触关系；

依据所述三维地质结构模型对所述待插值区域进行地质条件分块，得到多个三维地质分块；

采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据；所述改进的反距离加权插值算法为考虑高程修正的反距离加权插值算法；

采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据；

构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据；所述完整时间序列是由不同时刻的三维空间属性数据构成的；

采用协方差模型将所述三维空间属性数据和对应的时间属性数据进行融合，得到三维地质属性模型；所述三维地质属性模型用于确定待插值区域的地质空间分布和变化规律。

可选的，所述采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据，具体包括：

对所述地质采集点属性数据进行正态分布校验，确定各所述三维地质分块对应的地质采集点属性校验数据；所述地质采集点属性校验数据为服从正态分布的地质采集点属性数据；

采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性校验数据分别进行水平空间插值，得到各所述三维地质分块的水平插补属性值；所述水平插补属性值为在水平空间上的待插值点的估计值；所述改进的反距离加权插值算法的计算公式为

其中，Z为水平空间上的待插值点的估计值，Z_i为第i个采样点的实测值，λ_i为采样点的权重系数，D_i为第i个采样点与待插值点的高程差，d_i为第i个采样点与待插值点的距离，m为参与计算的实测采样点个数，n为幂指数，x_i为第i个采样点的横坐标，x为待插值点的横坐标，y_i为第i个采样点的纵坐标，y为待插值点的纵坐标；

依据各所述三维地质分块对应的地质采集点属性校验数据和各所述三维地质分块的水平插补属性值，确定各所述三维地质分块对应的水平空间属性数据。

可选的，所述采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据，具体包括：

采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块的垂直插补属性值；所述垂直插补属性值为在立体空间上的待插值点的估计值；所述空间圆弧插补方法的计算公式为

其中，为立体空间上的待插值点的估计值，g为参与计算的水平空间属性数据对应的采样点的个数，w_j为第j个采样点的属性值对缺失数据的空间贡献权重，y_j为第j个采样点的属性值，水平空间属性数据由h个采样点的属性值构成；g≤h；

依据各所述三维地质分块对应的水平空间属性数据和各所述三维地质分块的垂直插补属性值，确定各所述三维地质分块对应的三维空间属性数据。

可选的，所述构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据，具体包括：

构建各所述三维空间属性数据对应的完整时间序列；所述完整时间序列为

X＝{x(t₀),x(t₁),x(t₂),…,x(t_s),…,x(t_r)}，

其中，X为采样间隔为Δt的完整时间序列，Δt＝t_s-t_s-1，x(t₀)为在t₀时刻的三维空间属性数据，x(t_s)为在t_s时刻的三维空间属性数据，x(t_r)为在t_r时刻的三维空间属性数据；

采用差分法对所述完整时间序列进行平稳化处理，得到平稳时间序列；所述平稳时间序列为

W_t＝{w₀,w₁,w₂,…,w_t,…,w_v}，

其中，W_t为平稳时间序列，w₀为在初始时刻的平稳三维空间属性数据，w_t为在t时刻的平稳三维空间属性数据；

采用训练好的求和自回归移动平均模型对所述平稳时间序列中时间插值点的属性值进行预测，得到时间插值点属性预测值；所述求和自回归移动平均模型具体为

其中，为t时刻的时间插值点属性预测值，实参数φ₁,φ₂,…,φ_p为自回归系数，非负整数p为自回归项数，实参数θ₁,θ₂,…,θ_q为移动平均系数，非负整数q为移动平均阶次，d为使所述完整时间序列成为平稳序列所做的差分阶次，{ε_t}为白噪声序列；

依据所述完整时间序列和所述时间插值点属性预测值，得到时间属性数据。

可选的，所述采用协方差模型将各所述三维地质分块对应的三维空间属性数据和对应的时间属性数据进行拟合，得到三维地质属性模型，具体为：

其中，A表示空间维度权重，B表示时间维度权重，为t＝0时刻的时间插值点属性预测值。

本发明还提供了一种三维地质属性模型时空插值系统，包括：

数据获取模块，用于获取待插值区域的地形数据、钻孔数据、地质剖面数据和地质采集点属性数据；

地质结构模型构建模块，用于依据所述地形数据、所述钻孔数据和所述地质剖面数据，生成三维地质结构模型；所述三维地质结构模型表示所述待插值区域内断层、地层和岩体三者之间的接触关系；

分块模块，用于依据所述三维地质结构模型对所述待插值区域进行地质条件分块，得到多个三维地质分块；

水平插值模块，用于采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据；所述改进的反距离加权插值算法为考虑高程修正的反距离加权插值算法；

垂直插值模块，用于采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据；

时间插值模块，用于构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据；所述完整时间序列是由不同时刻的三维空间属性数据构成的；

拟合模块，用于采用协方差模型将所述三维空间属性数据和对应的时间属性数据进行融合，得到三维地质属性模型；所述三维地质属性模型用于确定待插值区域的地质空间分布和变化规律。

可选的，所述水平插值模块，具体包括：

正态校正单元，用于对所述地质采集点属性数据进行正态分布校验，确定各所述三维地质分块对应的地质采集点属性校验数据；所述地质采集点属性校验数据为服从正态分布的地质采集点属性数据；

水平插值单元，用于采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性校验数据分别进行水平空间插值，得到各所述三维地质分块的水平插补属性值；所述水平插补属性值为在水平空间上的待插值点的估计值；所述改进的反距离加权插值算法的计算公式为

其中，Z为水平空间上的待插值点的估计值，Z_i为第i个采样点的实测值；λ_i为采样点的权重系数，D_i为第i个采样点与待插值点的高程差，d_i为第i个采样点与待插值点的距离，m为参与计算的实测采样点个数，n为幂指数，x_i为第i个采样点的横坐标，x为待插值点的横坐标；y_i为第i个采样点的纵坐标，y为待插值点的纵坐标；

第一确定单元，用于依据各所述三维地质分块对应的地质采集点属性校验数据和各所述三维地质分块的水平插补属性值，确定各所述三维地质分块对应的水平空间属性数据。

可选的，所述垂直插值模块，具体包括：

垂直插值单元，用于采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块的垂直插补属性值；所述垂直插补属性值为在立体空间上的待插值点的估计值；所述空间圆弧插补方法的计算公式为

其中，为立体空间上的待插值点的估计值，g为参与计算的水平空间属性数据对应的采样点的个数，w_j为第j个采样点的属性值对缺失数据的空间贡献权重，y_j为第j个采样点的属性值，水平空间属性数据由h个采样点的属性值构成；g≤h；

第二确定单元，用于依据各所述三维地质分块对应的水平空间属性数据和各所述三维地质分块的垂直插补属性值，确定各所述三维地质分块对应的三维空间属性数据。

可选的，所述时间插值模块，具体包括：

时间序列构建单元，用于构建各所述三维空间属性数据对应的完整时间序列；所述完整时间序列为

X＝{x(t₀),x(t₁),x(t₂),…,x(t_s),…,x(t_r)}，

其中，X为采样间隔为Δt的完整时间序列，Δt＝t_s-t_s-1，x(t₀)为在t₀时刻的三维空间属性数据，x(t_s)为在t_s时刻的三维空间属性数据，x(t_r)为在t_r时刻的三维空间属性数据；

平稳化处理单元，用于采用差分法对所述完整时间序列进行平稳化处理，得到平稳时间序列；所述平稳时间序列为

W_t＝{w₀,w₁,w₂,…,w_t,…,w_v}，

其中，W_t为平稳时间序列，w₀为在初始时刻的平稳三维空间属性数据，w_t为在t时刻的平稳三维空间属性数据；

属性值预测单元，用于采用训练好的求和自回归移动平均模型对所述平稳时间序列中时间插值点的属性值进行预测，得到时间插值点属性预测值；所述求和自回归移动平均模型具体为

其中，为t时刻的时间插值点属性预测值，实参数φ₁,φ₂,…,φ_p为自回归系数，非负整数p为自回归项数，实参数θ₁,θ₂,…,θ_q为移动平均系数，非负整数q为移动平均阶次，d为使所述完整时间序列成为平稳序列所做的差分阶次，{ε_t}为白噪声序列；

第三确定单元，用于依据所述完整时间序列和所述时间插值点属性预测值，得到时间属性数据。

可选的，所述拟合模块，具体为：

其中，A表示空间维度权重，B表示时间维度权重，为t＝0时刻的时间插值点属性预测值。

与现有技术相比，本发明的有益效果是：

本发明提出了一种三维地质属性模型时空插值方法及系统，所述方法包括：依据地形数据、钻孔数据和地质剖面数据，生成三维地质结构模型；对待插值区域进行分块，得到多个三维地质分块；对各三维地质分块对应的地质采集点属性数据依次进行水平空间插值和垂直空间插值，得到三维空间属性数据；构建各三维地质分块对应的完整时间序列，并采用采用差分法对完整时间序列进行平稳化处理，得到平稳时间序列，采用求和自回归移动平均模型对时间插值点的属性值进行预测，得到时间属性数据；采用协方差模型将各三维地质分块对应的三维空间属性数据和对应的时间属性数据进行拟合，得到三维地质属性模型。本发明解决了传统的平面插值法(如反距离加权法、克里金插值法等)未考虑时间和空间要素影响的问题，提高了插值精度，保证了插值后时空数据的完整性，进而更准确的确定待插值区域的地质空间分布和变化规律。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例一种三维地质属性模型时空插值方法的流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明实施例一种三维地质属性模型时空插值方法的流程图。

参见图1，实施例的三维地质属性模型时空插值方法，包括：

步骤S1：获取待插值区域的地形数据、钻孔数据、地质剖面数据和地质采集点属性数据。

步骤S2：依据所述地形数据、所述钻孔数据和所述地质剖面数据，生成三维地质结构模型。

所述三维地质结构模型表示所述待插值区域内断层、地层和岩体三者之间的接触关系。

所述步骤S2，具体包括：

21：根据地形数据形成地形表面；根据钻孔数据和地质剖面形成各地层面、断层面和岩体表面。

22：根据断层的性质和断层之间的接触关系，建立断层接触模型；根据断层与地层面之间的接触关系，建立地层-断层接触模型；根据断层与岩体之间的接触关系，建立岩体-断层接触模型(区域若无岩体可不考虑)；根据岩体与地层之间的接触关系，建立岩体-地层接触模型(区域若无岩体可不考虑)。

23：判断以上生成的模型与地质规律有无矛盾，对于不符合地质规律的部分进行修改，生成三维地质结构模型。

步骤S3：依据所述三维地质结构模型对所述待插值区域进行地质条件分块，得到多个三维地质分块。

所述步骤S3，具体包括：

由于断层两侧的地层、岩性可能差别巨大，在插值过程中对待插点估值影响强烈，因此首先利用断层接触模型对待插值三维区域进行构造分区。在此基础上，对地层的岩性进行梳理，对岩性差异较小的地层进行合并处理；根据地层的空间分布情况，对构造分区的结果再次进行分块处理。岩体主要指岩浆岩，其岩性与围岩有明显差异，因此需要对岩体进行分区处理；根据岩体三维结构模型将岩体范围单独分区；通过构造分区、地层分区和岩体分区，待插值区域被分成若干地质条件各异的三维地质分块。

步骤S4：采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据。

所述改进的反距离加权插值算法为考虑高程修正的反距离加权插值算法。

所述步骤S4，具体包括：

41：对所述地质采集点属性数据进行正态分布校验，确定各所述三维地质分块对应的地质采集点属性校验数据；所述地质采集点属性校验数据为服从正态分布的地质采集点属性数据。具体为：

首先，对地质采集点属性数据进行正态分布检验，可通过数理统计如峰度、偏度、QQ图法等，判断地质采集点属性数据是否服从正态分布。对于不服从正态分布的数据进行异常值处理，如剔除异常值、替代异常值等。根据钻孔信息求得区域化变量的空间变异规律-变异函数，变异函数云图是针对样品对(地质属性值/间距)的变异函数值的分布，根据分布图直观的查看样品对(地质属性值/间距)分布情况，找出异常数据的位置，在插值时剔除或替代异常值，提升预测准确性。剔除异常值是指将异常值所在属性点的值删除；替代异常值是指用其他值替代异常值，如均值、最大值、中值等，应根据插值属性的具体类型进行设定。对于异常值处理后的数据再次进行正态分布检验，若数据接近正态分布，且数据相关性高，说明该数据进行空间插值有意义。

42：对于不同的三维地质分块，应依次利用反距离加权(IDW)插值算法进行插值计算。由于不同采样点所在高程不同，待插点与采样点常常不在同一高程平面上，需要进行高程修正。

本实施例中，采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性校验数据分别进行水平空间插值，得到各所述三维地质分块的水平插补属性值。所述水平插补属性值为在水平空间上的待插值点的估计值。所述改进的反距离加权插值算法的计算公式为

其中，Z为水平空间上的待插值点的估计值；Z_i为第i个采样点的实测值；λ_i为采样点的权重系数；D_i为第i个采样点与待插值点的高程差；d_i为第i个采样点与待插值点的距离；m为参与计算的实测采样点个数；n为幂指数，控制权重系数随待插点与采样点之间距离的增加而下降的程度，n越大，邻近的采样点赋予更高的权重，n越小，权重更均匀分配给各采样点；x_i为第i个采样点的横坐标；x为待插值点的横坐标；y_i为第i个采样点的纵坐标；y为待插值点的纵坐标。

43：依据各所述三维地质分块对应的地质采集点属性校验数据和各所述三维地质分块的水平插补属性值，确定各所述三维地质分块对应的水平空间属性数据。

经过步骤41-43，将分布不均匀的地质采集点属性数据转化为在水平方向上分布均匀的数据，即水平空间属性数据。所述水平空间属性数据可能是由部分地质采集点属性校验数据和全部的水平插补属性值构成的，也可能是由全部的水平插补属性值构成。这取决于在插值过程中，地质采集点属性校验数据是否位于网格交点，如果是，则地质采集点属性校验数据属于水平空间属性数据。

步骤S5：采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据。

所述这步骤S5，具体包括：

51：采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块的垂直插补属性值。所述垂直插补属性值为在立体空间上的待插值点的估计值。

下面对空间圆弧插补方法进行介绍。空间圆弧插补方法是根据空间不在一条直线上的三点确定下来的圆弧，通过空间矢量的变换、计算来进行插补实现的一种手段。其基本流程为：

1)通过不在一条直线上的空间三点确定圆弧所在圆的圆心坐标；

2)通过圆心与三点中任意一点的距离确定所在圆的半径；

3)通过圆心与起点、终点所在矢量确定圆弧的圆心角大小；

4)通过圆心角、速度、插补周期、半径、三点坐标、圆心坐标等标量与矢量计算出当前时刻上的插补坐标。

所述空间圆弧插补方法的计算公式为

其中，为立体空间上的待插值点的估计值，g为参与计算的水平空间属性数据对应的采样点的个数，w_j为第j个采样点的属性值对缺失数据的空间贡献权重，y_j为第j个采样点的属性值，水平空间属性数据由h个采样点的属性值构成；g≤h。

本实施例中，g＝12，根据三维空间中与待插值点距离最近的12个采样点的水平空间属性数据进行插值计算，选择12个采样点的数据，是兼顾算法效率和插值效果的经验值。若参与计算的数据过多，则导致插值效率低下；若数据过少，则影响插值精度，导致插值结果不准确。

利用待插值点周围12个相关系数最大的周围空间采样数据并按照上述方法可插补计算得出每个三维地质分块在不同时刻的的垂直插补属性值。

52：依据各所述三维地质分块对应的水平空间属性数据和各所述三维地质分块的垂直插补属性值，确定各所述三维地质分块对应的三维空间属性数据。

经过步骤51-52，将垂直方向上分布不均匀的水平空间属性数据转化为在水平方向上和垂直方向上均分布均匀的数据，即三维空间属性数据。所述三维空间属性数据可能是由部分水平空间属性数据和全部的垂直插补属性值构成的，也可能是由全部的垂直插补属性值构成。这取决于在插值过程中，水平空间属性数据是否位于网格交点，如果是，则水平空间属性数据属于三维空间属性数据。

步骤S6：构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据。

所述完整时间序列是由不同时刻的三维空间属性数据构成的。

所述步骤S6，具体包括：

61：构建各所述三维地质分块对应的完整时间序列；所述完整时间序列是一种将具有统计指标的数值按时间顺序排列形成的数列；所述完整时间序列为

X＝{x(t₀),x(t₁),x(t₂),…,x(t_s),…,x(t_r)}，

其中，X为采样间隔为Δt的完整时间序列，Δt＝t_s-t_s-1，x(t₀)为在t₀时刻的三维空间属性数据，x(t_s)为在t_s时刻的三维空间属性数据，x(t_r)为在t_r时刻的三维空间属性数据。

通过编制和统计分析原始地质属性数据的采样特征和规律，将某一变量的属性值，按照一系列时刻(t₀，t₁，t₂，…，t_r)排序，建立时间序列。

62：分析完整时间序列平稳性，采用差分法对完整时间序列(非平稳时间序列)数据进行平稳化处理，得到平稳时间序列。所述平稳时间序列为

W_t＝{w₀,w₁,w₂,…,w_t,…,w_v}，

其中，W_t为平稳时间序列，w₀为在初始时刻的平稳三维空间属性数据，w_t为在t时刻的平稳三维空间属性数据。

所述平稳时间序列W_t中包括X序列中t_r个时间点对应的三维空间属性数据和插值后的时间插值点(缺失时刻)对应的时间属性值。

步骤62，具体包括：

对于地质数据的时间要素来说，建立完整时间序列需考虑四个因素：

(1)长期趋势：是指在一个相当长的时间内表现为一种近似直线的持续向上或向下或平稳的趋势，比如沉降变形、地下水水位、地层结构、含水量、孔隙比。

(2)季节变动：是指收季节变动影响所形成的一种长度和幅度固定的短周期波动。比如：降雨量、地下水径流量、降水入渗补给量、地表水入渗补给量等。

(3)循环变动：指一种较长时间的上下起伏周期性波动。比如：地下水质量中微生物指标、放射性指标、毒理性指标。

(4)不规则变动：指受偶然因素影响所形成的不规则波动，在时间序列中无法预计。比如地质灾害中地震灾害、泥石流灾害、坍塌灾害等。

基于上述四大因素，时间插值点的属性值的预测步骤主要分为四步：

(1)绘制时间序列图观察趋势；

(2)分析序列平稳性并进行平稳化。对于地质数据来说，先对数据准确性和完整性进行分析，在对(1)建立的时间序列进行平稳性识别，在实际的地质应用中大多数时间序列都是不平稳的，故需要做平稳化处理，采用差分法，利用序列中前后相邻的两期数据之差对时间序列进行插补，形成一个完整的平稳时间序列。

一次差分表达式：

二阶差分法表达式：

其中，w_t为在t时刻的三维空间属性数据，w_t-1为在t-1时刻的三维空间属性数据，为差分算子。

63：时间序列建模分析。采用训练好的求和自回归移动平均模型对所述平稳时间序列中时间插值点的时间属性值进行预测，得到时间插值点属性预测值；所述求和自回归移动平均模型具体为

其中，为t时刻的时间插值点属性预测值，实参数φ₁,φ₂,…,φ_p为自回归系数，非负整数p为自回归项数，实参数θ₁,θ₂,…,θ_q为移动平均系数，非负整数q为移动平均阶次，d为使所述完整时间序列成为平稳序列所做的差分阶次，{ε_t}为白噪声序列。

64：依据所述完整时间序列和所述时间插值点属性预测值，得到时间属性数据。所述时间属性数据由t_r个时间点对应的三维空间属性数据和时间插值点(缺失时刻)对应的时间插值点属性预测值构成。

步骤S7：采用协方差模型将所述三维空间属性数据和对应的时间属性数据进行融合，得到三维地质属性模型。

所述步骤S7中，协方差分析模型是称带有协变量的方差分析，是将线性回归与方差分析综合运用的一种统计方法。时间和空间纬度的估算值采用线性加权方法对其进行数据融合，步骤S7的具体计算公式为：

其中，A表示空间维度权重，B表示时间维度权重，为t＝0时刻的时间插值点属性预测值。

本实施例通过上述步骤S1-S7，实现了对地质模型时空属性的融合，即得到三维地质属性模型，所述三维地质属性模型在后续的实际应用中，可确定待插值区域的地质空间分布和变化规律。

本实施例的三维地质属性模型时空插值方法，解决了传统的平面插值法(如反距离加权法、克里金插值法等)未考虑时间和空间要素影响的技术问题，先根据三维地质结构模型中断层、地层、岩体的空间位置和接触关系，进行地质条件分区；利用异常值处理让时空属性数据接近标准正态分布，利用高程修正的反距离加权法依次对不同的地质条件分区进行水平方向上的插值计算，利用空间圆弧插补方法进行垂直方向上的插值计算；基于时间序列分析，利用差分法对非平稳时间序列进行平稳化处理，利用求和自回归移动平均模型进行时间序列建模分析、评估和预测；利用协方差模型实现时空属性的融合，完成地质模型的属性赋值。解决了传统算法精度低的技术问题，将地质属性模型融入更多要素，使其表现更贴近真实地质结构，解决了目前地质领域插值后时间属性丢失的技术问题。该方法能够基于地质分布规律，提高插值精度，保证插值后时空数据的完整性。

本发明还提供了一种三维地质属性模型时空插值系统，包括：

数据获取模块，用于获取待插值区域的地形数据、钻孔数据、地质剖面数据和地质采集点属性数据。

地质结构模型构建模块，用于依据所述地形数据、所述钻孔数据和所述地质剖面数据，生成三维地质结构模型；所述三维地质结构模型表示所述待插值区域内断层、地层和岩体三者之间的接触关系。

分块模块，用于依据所述三维地质结构模型对所述待插值区域进行地质条件分块，得到多个三维地质分块。

水平插值模块，用于采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据；所述改进的反距离加权插值算法为考虑高程修正的反距离加权插值算法。

垂直插值模块，用于采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据。

时间插值模块，用于构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据；所述完整时间序列是由不同时刻的三维空间属性数据构成的。

拟合模块，用于采用协方差模型将所述三维空间属性数据和对应的时间属性数据进行融合，得到三维地质属性模型；所述三维地质属性模型用于确定待插值区域的地质空间分布和变化规律。

作为一种可选的实施方式，所述水平插值模块，具体包括：

正态校正单元，用于对所述地质采集点属性数据进行正态分布校验，确定各所述三维地质分块对应的地质采集点属性校验数据；所述地质采集点属性校验数据为服从正态分布的地质采集点属性数据。

水平插值单元，用于采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性校验数据分别进行水平空间插值，得到各所述三维地质分块的水平插补属性值；所述水平插补属性值为在水平空间上的待插值点的估计值；所述改进的反距离加权插值算法的计算公式为

其中，Z为水平空间上的待插值点的估计值，Z_i为第i个采样点的实测值；λ_i为采样点的权重系数，D_i为第i个采样点与待插值点的高程差，d_i为第i个采样点与待插值点的距离，m为参与计算的实测采样点个数，n为幂指数，x_i为第i个采样点的横坐标，x为待插值点的横坐标；y_i为第i个采样点的纵坐标，y为待插值点的纵坐标。

第一确定单元，用于依据各所述三维地质分块对应的地质采集点属性校验数据和各所述三维地质分块的水平插补属性值，确定各所述三维地质分块对应的水平空间属性数据。

作为一种可选的实施方式，所述垂直插值模块，具体包括：

垂直插值单元，用于采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块的垂直插补属性值；所述垂直插补属性值为在立体空间上的待插值点的估计值；所述空间圆弧插补方法的计算公式为

其中，为立体空间上的待插值点的估计值，g为参与计算的水平空间属性数据对应的采样点的个数，w_j为第j个采样点的属性值对缺失数据的空间贡献权重，y_j为第j个采样点的属性值，水平空间属性数据由h个采样点的属性值构成；g≤h。

第二确定单元，用于依据各所述三维地质分块对应的水平空间属性数据和各所述三维地质分块的垂直插补属性值，确定各所述三维地质分块对应的三维空间属性数据。

作为一种可选的实施方式，所述时间插值模块，具体包括：

时间序列构建单元，用于构建各所述三维空间属性数据对应的完整时间序列；所述完整时间序列为

X＝{x(t₀),x(t₁),x(t₂),…,x(t_s),…,x(t_r)}，

其中，X为采样间隔为Δt的完整时间序列，Δt＝t_s-t_s-1，x(t₀)为在t₀时刻的三维空间属性数据，x(t_s)为在t_s时刻的三维空间属性数据，x(t_r)为在t_r时刻的三维空间属性数据。

平稳化处理单元，用于采用差分法对所述完整时间序列进行平稳化处理，得到平稳时间序列；所述平稳时间序列为

W_t＝{w₀,w₁,w₂,…,w_t,…,w_v}，

其中，W_t为平稳时间序列，w₀为在初始时刻的平稳三维空间属性数据，w_t为在t时刻的平稳三维空间属性数据。

属性值预测单元，用于采用训练好的求和自回归移动平均模型对所述平稳时间序列中时间插值点的属性值进行预测，得到时间插值点属性预测值；所述求和自回归移动平均模型具体为

其中，为t时刻的时间插值点属性预测值，实参数φ₁,φ₂,…,φ_p为自回归系数，非负整数p为自回归项数，实参数θ₁,θ₂,…,θ_q为移动平均系数，非负整数q为移动平均阶次，d为使所述完整时间序列成为平稳序列所做的差分阶次，{ε_t}为白噪声序列。

第三确定单元，用于依据所述完整时间序列和所述时间插值点属性预测值，得到时间属性数据。

作为一种可选的实施方式，所述拟合模块，具体为：

其中，A表示空间维度权重，B表示时间维度权重，为t＝0时刻的时间插值点属性预测值。

本实施例的三维地质属性模型时空插值系统，能够基于地质分布规律，提高插值精度，保证插值后时空数据的完整性。

对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种三维地质属性模型时空插值方法，其特征在于，包括：

获取待插值区域的地形数据、钻孔数据、地质剖面数据和地质采集点属性数据；

依据所述地形数据、所述钻孔数据和所述地质剖面数据，生成三维地质结构模型；所述三维地质结构模型表示所述待插值区域内断层、地层和岩体三者之间的接触关系；

依据所述三维地质结构模型对所述待插值区域进行地质条件分块，得到多个三维地质分块；

采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据；所述改进的反距离加权插值算法为考虑高程修正的反距离加权插值算法；

采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据；

构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据；所述完整时间序列是由不同时刻的三维空间属性数据构成的；

采用协方差模型将所述三维空间属性数据和对应的时间属性数据进行融合，得到三维地质属性模型；所述三维地质属性模型用于确定待插值区域的地质空间分布和变化规律。

2.根据权利要求1所述的一种三维地质属性模型时空插值方法，其特征在于，所述采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据，具体包括：

对所述地质采集点属性数据进行正态分布校验，确定各所述三维地质分块对应的地质采集点属性校验数据；所述地质采集点属性校验数据为服从正态分布的地质采集点属性数据；

采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性校验数据分别进行水平空间插值，得到各所述三维地质分块的水平插补属性值；所述水平插补属性值为在水平空间上的待插值点的估计值；所述改进的反距离加权插值算法的计算公式为

其中，Z为水平空间上的待插值点的估计值，Z_i为第i个采样点的实测值，λ_i为采样点的权重系数，D_i为第i个采样点与待插值点的高程差，d_i为第i个采样点与待插值点的距离，m为参与计算的实测采样点个数，n为幂指数，x_i为第i个采样点的横坐标，x为待插值点的横坐标，y_i为第i个采样点的纵坐标，y为待插值点的纵坐标；

依据各所述三维地质分块对应的地质采集点属性校验数据和各所述三维地质分块的水平插补属性值，确定各所述三维地质分块对应的水平空间属性数据。

3.根据权利要求1所述的一种三维地质属性模型时空插值方法，其特征在于，所述采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据，具体包括：

采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块的垂直插补属性值；所述垂直插补属性值为在立体空间上的待插值点的估计值；所述空间圆弧插补方法的计算公式为

其中，为立体空间上的待插值点的估计值，g为参与计算的水平空间属性数据对应的采样点的个数，w_j为第j个采样点的属性值对缺失数据的空间贡献权重，y_j为第j个采样点的属性值，水平空间属性数据由h个采样点的属性值构成；g≤h；

依据各所述三维地质分块对应的水平空间属性数据和各所述三维地质分块的垂直插补属性值，确定各所述三维地质分块对应的三维空间属性数据。

4.根据权利要求3所述的一种三维地质属性模型时空插值方法，其特征在于，所述构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据，具体包括：

构建各所述三维空间属性数据对应的完整时间序列；所述完整时间序列为

X＝{x(t₀),x(t₁),x(t₂),…,x(t_s),…,x(t_r)}，

其中，X为采样间隔为Δt的完整时间序列，Δt＝t_s-t_s-1，x(t₀)为在t₀时刻的三维空间属性数据，x(t_s)为在t_s时刻的三维空间属性数据，x(t_r)为在t_r时刻的三维空间属性数据；

采用差分法对所述完整时间序列进行平稳化处理，得到平稳时间序列；所述平稳时间序列为

W_t＝{w₀,w₁,w₂,…,w_t,…,w_v}，

其中，W_t为平稳时间序列，w₀为在初始时刻的平稳三维空间属性数据，w_t为在t时刻的平稳三维空间属性数据；

采用训练好的求和自回归移动平均模型对所述平稳时间序列中时间插值点的属性值进行预测，得到时间插值点属性预测值；所述求和自回归移动平均模型具体为

其中，为t时刻的时间插值点属性预测值，实参数φ₁,φ₂,…,φ_p为自回归系数，非负整数p为自回归项数，实参数θ₁,θ₂,…,θ_q为移动平均系数，非负整数q为移动平均阶次，d为使所述完整时间序列成为平稳序列所做的差分阶次，{ε_t}为白噪声序列；

依据所述完整时间序列和所述时间插值点属性预测值，得到时间属性数据。

5.根据权利要求4所述的一种三维地质属性模型时空插值方法，其特征在于，所述采用协方差模型将各所述三维地质分块对应的三维空间属性数据和对应的时间属性数据进行拟合，得到三维地质属性模型，具体为：

其中，A表示空间维度权重，B表示时间维度权重，为t＝0时刻的时间插值点属性预测值。

6.一种三维地质属性模型时空插值系统，其特征在于，包括：

数据获取模块，用于获取待插值区域的地形数据、钻孔数据、地质剖面数据和地质采集点属性数据；

地质结构模型构建模块，用于依据所述地形数据、所述钻孔数据和所述地质剖面数据，生成三维地质结构模型；所述三维地质结构模型表示所述待插值区域内断层、地层和岩体三者之间的接触关系；

分块模块，用于依据所述三维地质结构模型对所述待插值区域进行地质条件分块，得到多个三维地质分块；

水平插值模块，用于采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性数据分别进行水平空间插值，得到各所述三维地质分块对应的水平空间属性数据；所述改进的反距离加权插值算法为考虑高程修正的反距离加权插值算法；

垂直插值模块，用于采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块对应的三维空间属性数据；

时间插值模块，用于构建各所述三维空间属性数据对应的完整时间序列，并采用求和自回归移动平均模型对所述完整时间序列进行建模分析和插补，得到时间属性数据；所述完整时间序列是由不同时刻的三维空间属性数据构成的；

拟合模块，用于采用协方差模型将所述三维空间属性数据和对应的时间属性数据进行融合，得到三维地质属性模型；所述三维地质属性模型用于确定待插值区域的地质空间分布和变化规律。

7.根据权利要求6所述的一种三维地质属性模型时空插值系统，其特征在于，所述水平插值模块，具体包括：

正态校正单元，用于对所述地质采集点属性数据进行正态分布校验，确定各所述三维地质分块对应的地质采集点属性校验数据；所述地质采集点属性校验数据为服从正态分布的地质采集点属性数据；

水平插值单元，用于采用改进的反距离加权插值算法对各所述三维地质分块对应的地质采集点属性校验数据分别进行水平空间插值，得到各所述三维地质分块的水平插补属性值；所述水平插补属性值为在水平空间上的待插值点的估计值；所述改进的反距离加权插值算法的计算公式为

其中，Z为水平空间上的待插值点的估计值，Z_i为第i个采样点的实测值；λ_i为采样点的权重系数，D_i为第i个采样点与待插值点的高程差，d_i为第i个采样点与待插值点的距离，m为参与计算的实测采样点个数，n为幂指数，x_i为第i个采样点的横坐标，x为待插值点的横坐标；y_i为第i个采样点的纵坐标，y为待插值点的纵坐标；

第一确定单元，用于依据各所述三维地质分块对应的地质采集点属性校验数据和各所述三维地质分块的水平插补属性值，确定各所述三维地质分块对应的水平空间属性数据。

8.根据权利要求6所述的一种三维地质属性模型时空插值系统，其特征在于，所述垂直插值模块，具体包括：

垂直插值单元，用于采用空间圆弧插补方法对各所述三维地质分块对应的水平空间属性数据进行垂直空间插值，得到各所述三维地质分块的垂直插补属性值；所述垂直插补属性值为在立体空间上的待插值点的估计值；所述空间圆弧插补方法的计算公式为

其中，为立体空间上的待插值点的估计值，g为参与计算的水平空间属性数据对应的采样点的个数，w_j为第j个采样点的属性值对缺失数据的空间贡献权重，y_j为第j个采样点的属性值，水平空间属性数据由h个采样点的属性值构成；g≤h；

第二确定单元，用于依据各所述三维地质分块对应的水平空间属性数据和各所述三维地质分块的垂直插补属性值，确定各所述三维地质分块对应的三维空间属性数据。

9.根据权利要求8所述的一种三维地质属性模型时空插值系统，其特征在于，所述时间插值模块，具体包括：

时间序列构建单元，用于构建各所述三维空间属性数据对应的完整时间序列；所述完整时间序列为

X＝{x(t₀),x(t₁),x(t₂),…,x(t_s),…,x(t_r)}，

其中，X为采样间隔为Δt的完整时间序列，Δt＝t_s-t_s-1，x(t₀)为在t₀时刻的三维空间属性数据，x(t_s)为在t_s时刻的三维空间属性数据，x(t_r)为在t_r时刻的三维空间属性数据；

平稳化处理单元，用于采用差分法对所述完整时间序列进行平稳化处理，得到平稳时间序列；所述平稳时间序列为

W_t＝{w₀,w₁,w₂,…,w_t,…,w_v}，

其中，W_t为平稳时间序列，w₀为在初始时刻的平稳三维空间属性数据，w_t为在t时刻的平稳三维空间属性数据；

属性值预测单元，用于采用训练好的求和自回归移动平均模型对所述平稳时间序列中时间插值点的属性值进行预测，得到时间插值点属性预测值；所述求和自回归移动平均模型具体为

其中，为t时刻的时间插值点属性预测值，实参数φ₁,φ₂,…,φ_p为自回归系数，非负整数p为自回归项数，实参数θ₁,θ₂,…,θ_q为移动平均系数，非负整数q为移动平均阶次，d为使所述完整时间序列成为平稳序列所做的差分阶次，{ε_t}为白噪声序列；

第三确定单元，用于依据所述完整时间序列和所述时间插值点属性预测值，得到时间属性数据。

10.根据权利要求9所述的一种三维地质属性模型时空插值系统，其特征在于，所述拟合模块，具体为：

其中，A表示空间维度权重，B表示时间维度权重，为t＝0时刻的时间插值点属性预测值。