一种利用相邻地层互相关约束的地层面模型修正方法

1.本发明涉及地统计学技术领域，尤其涉及一种利用相邻地层互相关约束的地层面模型修正方法。


背景技术：

2.多数城市建设工程都位于第四系沉积物覆盖区，为了科学合理地利用城市地下空间，有必要研究城市地下沉积层的几何形态和空间分布等信息，因而通过建立地层结构模型来定量描述地质结构特征。但是在地质调查过程中，受成本和技术上的限制，有时部分区域的测量数据会存在采样稀疏、样本分布不均等问题。基于这样的样本信息构建地层模型时，难免会出现有的地层样本较多，地层模型精度高，而有的地层样本稀疏，建模精度较低的情况。
3.由于第四系沉积层构造相对简单，很多研究和应用采用了显式建模方法对地层的结构和几何形态进行建模。显式建模方法往往直接对地层面建模，包括以下几种主要思路：1、基于钻孔和剖面数据拟合地层结构参数曲面；2、通过插值或者网格剖分的方法构建地层面的规则格网或不规则三角网；3、利用grid、tin或者两者的混合模拟各地质层分界面的多层dem。这些显式建模方法在地层建模时，往往先对各地层面单独建模，然后通过层与层之间的层序和接触关系，利用地层的拓扑关系将多层面模型组织在一起构成完整的地质体。事实上，无论是三角剖分还是参数剖分，独立创建的曲面之间的关系通常与地层叠置原理不一致。特别是，独立创建的表面往往相互交叉或在域内留下间隙，这违反了基本地貌有效性条件。这种每个地层单独建模的处理往往只利用了本地层的数据(如钻孔或剖面上对应地层面的接触位置)或借助其他类型的数据作为辅助约束，在采样较多的区域可以得到接近实际的模型，然而在采样信息较少的区域，得到的地层模型可能存在较大的不确定性。在实际生产中，样本在深层和浅层地层上的数量和空间分布并不一定对等，一般来说，地表和浅层的地层的地质调查数据和钻孔采样较多，出于成本因素考虑，深度越大的地层采样数据越少；在有些工程应用中，在部分区域也会出现深层地层有采样，而浅层位置未采样的情况(如水平井)。对于这些情况，受样本数量的限制，各地层的建模精度会有较大差异，精度较低的地层模型会降低模型的整体质量。
4.对于已经建好的模型，可以采用模型修正的方法使地质模型更接近真实情况。建模人员可以根据自身地质经验，通过增加地层面插值控制点的方式修改地质剖面，或增加包含地层位置信息的虚拟钻孔，利用修改后的剖面和虚拟钻孔约束，对模型中不符合要求的部分进行调整修正。此类基于建模初始数据的模型误差修正方法往往需要在建模数据补充新信息(如新增钻孔、虚拟钻孔、辅助剖面、航磁数据等)，再借助自动建模的方法实现模型的更新修正。此外还有基于建模中间结果的模型误差修正方法，需要用户直接对地质模型的中间结果(如建模参数)进行人工修改。
5.现有技术中，针对显式模型的修正或者更新方法是利用新增数据(新增钻孔、剖面)或其他类型的数据(坡度、产状、重力、磁场、地震数据等)作为地质建模的约束项，或将
专家知识(层厚、断层、褶皱的几何参数)转化为虚拟钻孔或剖面等数据约束，利用这些额外的信息约束重新计算模型，提升模型精度。然而，在很多时候，额外的数据难以获取，建模人员面对稀疏的样本数据，对研究区域的个人理解也存在很大的模糊性，难以通过人工交互编辑或新增经验数据的方法修正模型。综上，如何在不依赖新增数据的情况下进行显式模型的修正是亟待解决的问题。


技术实现要素：

6.有鉴于此，有必要提供一种利用相邻地层互相关约束的地层面模型修正方法，用以解决如何在不依赖新增数据的情况下进行显式模型的修正的问题。
7.本发明提供一种利用相邻地层互相关约束的地层面模型修正方法，包括：
8.获取观测样本数据和待修正地层界面的已有模型；
9.根据所述观测样本分别统计所述待修正地层界面和辅助地层界面的高程边际特性函数，其中，所述辅助地层界面为所述待修正地层界面的相邻地层界面；
10.估计所述待修正地层界面的估计概率密度函数，作为对应的高程先验概率密度；
11.根据所述辅助地层界面在待修正位置上的数据精度确定辅助地层数，并根据所述辅助地层数、所述待修正地层界面的高程边际特性函数、所述辅助地层界面的高程边际特性函数以及相邻地层界面的高程变量的copula函数，确定联合分布函数；
12.根据所述待修正地层界面的高程边际特性函数和所述联合分布函数，确定所述待修正地层界面至少一个坐标点对应的似然函数；
13.根据所述至少一个坐标点对应的似然函数，对所述待修正地层界面的高程先验概率密度进行贝叶斯更新并进行归一化处理，确定所述至少一个坐标点的高程后验概率密度；
14.根据所述至少一个坐标点的高程后验概率密度，确定对应的模型修正值，对所述待修正地层界面的已有模型中的高程值进行更新。
15.进一步地，所述根据所述观测样本分别统计待修正地层界面和辅助地层界面的高程边际特性函数包括：
16.分别将所述待修正地层界面和所述辅助地层界面的高程值范围划分为多个等间距区间；
17.分别统计所述待修正地层界面和所述辅助地层界面的观测样本中的高程值落入不同的等间距区间的频率数据；
18.根据所述频率数据，分别确定所述辅助地层界面关于高程值的第一边际概率密度函数和所述待修正地层界面关于高程值的第二边际概率密度函数；
19.根据所述频率数据，分别确定所述辅助地层界面的第一累积概率密度函数和所述待修正地层界面的第二累积概率密度函数；
20.将所述第一累积概率密度函数变换为[0，1]区间内均匀分布下的分位数，确定所述辅助地层界面对应的第一分位数，将所述第二累积概率密度函数变换为[0，1]区间内均匀分布下的分位数，确定所述待修正地层界面对应的第二分位数。
[0021]
进一步地，所述估计所述待修正地层界面的估计概率密度函数，作为对应的高程先验概率密度包括：
[0022]
获取所述待修正地层界面的地层接触样本数据；
[0023]
通过克里金方法，确定所述至少一个坐标点对应的克里金方差；
[0024]
根据所述待修正地层界面的已有模型，确定所述至少一个坐标点对应的期望高程；
[0025]
将所述克里金方差作为正态分布的方差、所述期望高程作为正态分布的期望，建立正态分布函数，其中，所述正态分布函数作为所述待修正地层界面的高程值在所述至少一个坐标点下的估计概率密度函数；
[0026]
将所述待修正地层界面的估计概率密度函数作为对应的所述高程先验概率密度。
[0027]
进一步地，所述辅助地层界面包括位于所述待修正地层界面上侧的第一辅助地层界面、位于所述待修正地层界面下侧的第二辅助地层界面，所述根据所述辅助地层界面在待修正位置上的数据精度确定辅助地层数包括：
[0028]
分别获取所述第一辅助地层界面的已有模型和所述第二辅助地层界面的已有模型；
[0029]
通过模型不确定性评价方法，分别确定所述第一辅助地层界面的已有模型中的第一高程误差、所述第二辅助地层界面的已有模型中的第二高程误差；
[0030]
根据所述第一高程误差和所述第二高程误差，分别判断所述第一辅助地层界面的已有模型和所述第二辅助地层界面的已有模型是否满足数据精度条件；
[0031]
若其中之一满足，则将满足所述数据精度条件的所述第一辅助地层界面或所述第二辅助地层界面作为所述辅助地层界面，所述辅助地层数为一层；
[0032]
若都满足，则所述第一辅助地层界面和所述第二辅助地层界面都作为所述辅助地层界面，所述辅助地层数为两层。
[0033]
进一步地，所述根据所述辅助地层数、所述待修正地层界面的高程边际特性函数、所述辅助地层界面的高程边际特性函数以及相邻地层界面的高程变量的copula函数，确定联合分布函数包括：
[0034]
若所述辅助地层数为一层，则根据所述第一边际概率密度函数、所述第二边际概率密度函数、以及所述相邻地层界面的高程变量的copula函数，确定所述联合分布函数，通过以下公式得到所述联合分布函数：
[0035]
f
a,b
(z
a
,z
b
)＝f
a
(z
a
)
·
f
b
(z
b
)
·
c
a,b
(u
a
,u
b
)
[0036]
其中，f
a,b
(z
a
,z
b
)表示所述辅助地层数为一层时的所述联合分布函数，f
a
(z
a
)表示所述第一边际概率密度函数，f
b
(z
b
)表示所述第二边际概率密度函数，c
a,b
(u
a
,u
b
)表示所述相邻地层界面的高程变量的copula函数，u
a
为所述第一分位数，u
b
为所述第二分位数。
[0037]
进一步地，所述根据所述待修正地层界面的高程边际特性函数和所述联合分布函数，确定所述待修正地层界面至少一个坐标点对应的似然函数包括：
[0038]
若所述辅助地层数为一层，则在所述第一辅助地层界面或所述第二辅助地层界面的已有模型中，确定所述至少一个坐标点对应的模型高程值；
[0039]
根据所述模型高程值、所述第二边际概率密度函数和所述联合分布函数，通过以下公式确定所述似然函数：
[0040]
[0041]
其中，l(z
b
|z
a
＝v
a
)、f
a|b
(z
a
＝v
a
|z
b
)表示所述待修正地层界面在所述至少一个坐标点处对应的所述似然函数，z
a
表示所述辅助地层界面对应的高程变量，z
b
表示所述待修正地层界面对应的高程变量，v
a
表示所述模型高程值，f
a，b
(z
a
＝v
a
，z
b
)表示所述辅助地层数为一层时的所述联合分布函数，f
b
(z
b
)表示所述第二边际概率密度函数。
[0042]
进一步地，所述第一边际概率密度函数包括所述第一辅助地层界面的边际概率密度函数、所述第二辅助地层界面的边际概率密度函数，所述根据所述辅助地层数、所述待修正地层界面的高程边际特性函数、所述辅助地层界面的高程边际特性函数以及相邻地层界面的高程变量的copula函数，确定联合分布函数包括：
[0043]
若所述辅助地层数为两层，则根据第二边际概率密度函数、所述第一辅助地层界面的边际概率密度函数、所述第二辅助地层界面的边际概率密度函数、以及所述相邻地层界面的高程变量的copula函数，确定所述联合分布函数，通过以下公式得到所述联合分布函数：
[0044]
f
a,b,c
(z
a
,z
b
,z
c
)＝f
a
(z
a
)
·
f
b
(z
b
)
·
f
c
(z
c
)
·
c
a,b
(u
a
,u
b
)
·
c
b,c
(u
b
,u
c
)
·
c
a,c|b
(u
a|b
,u
c|b
)
[0045]
其中，f
a，b，c
(z
a
，z
b
，z
c
)表示所述辅助地层数为两层时的所述联合分布函数，f
b
(z
b
)表示所述第二边际概率密度函数，f
a
(z
a
)表示所述第一辅助地层界面的边际概率密度函数，f
c
(z
c
)表示所述第二辅助地层界面的边际概率密度函数，c
a，b
(u
a
，u
b
)、c
b，c
(u
b
，u
c
)以及c
a，c|b
(u
a|b
，u
c|b
)表示所述相邻地层界面的高程变量的copula函数，u
a
为所述第一辅助地层界面对应的所述第一分位数，u
b
为所述第二辅助地层界面对应的所述第一分位数，u
c
为所述第二分位数。
[0046]
进一步地，所述根据所述待修正地层界面的高程边际特性函数和所述联合分布函数，确定所述待修正地层界面至少一个坐标点对应的似然函数包括：
[0047]
若所述辅助地层数为两层，则在所述第一辅助地层界面和所述第二辅助地层界面的已有模型中，确定所述至少一个坐标点对应的模型高程值；
[0048]
根据所述模型高程值、所述第二边际概率密度函数和所述联合分布函数，通过以下公式确定所述似然函数：
[0049][0050]
其中，l(z
b
|z
a
＝v
a
，z
c
＝v
c
)、f
a，c|b
(z
a
＝v
a
，z
c
＝v
c
|z
b
)表示所述待修正地层界面在所述至少一个坐标点处对应的所述似然函数，z
a
表示所述第一辅助地层界面对应的高程变量，z
b
表示所述待修正地层界面对应的高程变量，z
c
表示所述第二辅助地层界面对应的高程变量，v
a
表示所述第一辅助地层界面对应的所述模型高程值，v
c
表示所述第二辅助地层界面对应的所述模型高程值，f
a，b，c
(z
a
，z
b
，z
c
)表示所述辅助地层数为两层时的所述联合分布函数，f
b
(z
b
)表示所述第二边际概率密度函数。
[0051]
进一步地，所述根据所述至少一个坐标点对应的似然函数，对所述待修正地层界面的高程先验概率密度进行贝叶斯更新并进行归一化处理，确定所述至少一个坐标点的高程后验概率密度包括：
[0052]
在所述至少一个坐标点处，根据所述似然函数和所述高程先验概率密度的乘积，对所述高程先验概率密度进行贝叶斯更新：
[0053]
对更新后的高程先验概率密度进行归一化处理。
[0054]
进一步地，所述根据所述至少一个坐标点的高程后验概率密度，确定对应的模型修正值，对所述待修正地层界面的已有模型中的高程值进行更新包括：
[0055]
根据所述至少一个坐标点的高程后验概率密度，确定所述高程后验概率密度对应的条件期望值；
[0056]
将所述条件期望值作为所述至少一个坐标点在所述待修正地层界面的已有模型中对应的模型修正值；
[0057]
根据所述模型修正值，对所述至少一个坐标点在所述待修正地层界面的已有模型中的高程值进行更新。
[0058]
与现有技术相比，本发明的有益效果包括：首先，进行观测样本数据和待修正地层界面的已有模型的有效获取；然后，利用观测样本数据的统计特性确定高程边际特性函数，以此表示待修正地层界面和辅助地层界面的分布特性；进而，对待修正地层界面的已有模型进行概率密度的估计，以此有效作为待修正地层界面的先验分布，便于后续求后验分布；进而，利用copula方法构建相邻多个地层界面高程变量的相关结构模型和联合分布函数，以copula函数对地层相关结构进行建模，构建多地层界面高程联合分布，准确描述多个相邻地层的复杂相关结构，不需要新增其他类型辅助数据，直接利用地层间的相关结构约束修正提高已有地层面模型的建模精度；接着，计算待修正地层界面每一个坐标点的似然函数，便于后续对修正模型的先验分布进行更新；然后，利用似然函数对已有模型中的每一个坐标点进行先验分布的更新，确定每一个坐标点对应的后验分布，并进行归一化，保证得到准确的高程后验概率密度；最后，根据每一个坐标点的高程后验概率密度求得对应的模型修正值，并对其对应的高程值进行更新，以此对已有模型进行准确快速的修正。综上，本发明不需要额外采集辅助数据，也不需要人工交互编辑，利用相邻地层互相关约束的地层面模型修正方法及系统，将相邻地层的互相关结构作为一种约束条件引入模型修正，通过copula方法描述相邻地层的复杂相关结构，通过相关性进行模型修正，提高已有模型的建模精度，大大提高了算法的实用性和快速性。
附图说明
[0059]
图1为本发明提供的利用相邻地层互相关约束的地层面模型修正方法的流程示意图；
[0060]
图2为本发明提供的确定高程边际特性函数的流程示意图；
[0061]
图3为本发明提供的高程先验概率密度的流程示意图；
[0062]
图4为本发明提供的确定辅助地层数的流程示意图；
[0063]
图5为本发明提供的当辅助地层数为一层时确定联合分布函数的流程示意图；
[0064]
图6为本发明提供的当辅助地层数为一层时确定似然函数的流程示意图；
[0065]
图7为本发明提供的当辅助地层数为两层时确定联合分布函数的流程示意图；
[0066]
图8为本发明提供的当辅助地层数为两层时确定似然函数的流程示意图；
[0067]
图9为本发明提供的更新高程后验概率密度的流程示意图；
[0068]
图10为本发明提供的更新已有模型的流程示意图；
[0069]
图11为本发明实施例提供的地质钻孔数据示意图；
[0070]
图12为本发明提供的修正示意图一；
[0071]
图13为本发明提供的修正示意图二。
具体实施方式
[0072]
下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理，并非用于限定本发明的范围。
[0073]
实施例1
[0074]
本发明实施例提供了一种利用相邻地层互相关约束的地层面模型修正方法，结合图1来看，图1为本发明提供的利用相邻地层互相关约束的地层面模型修正方法的流程示意图，上述利用相邻地层互相关约束的地层面模型修正方法包括步骤s1至步骤s7，其中：
[0075]
在步骤s1中，获取观测样本数据和待修正地层界面的已有模型；
[0076]
在步骤s2中，根据观测样本分别统计待修正地层界面和辅助地层界面的高程边际特性函数，其中，辅助地层界面为待修正地层界面的相邻地层界面；
[0077]
在步骤s3中，估计待修正地层界面的估计概率密度函数，作为对应的高程先验概率密度；
[0078]
在步骤s4中，根据辅助地层界面在待修正位置上的数据精度确定辅助地层数，并根据辅助地层数、待修正地层界面的高程边际特性函数、辅助地层界面的高程边际特性函数以及相邻地层界面的高程变量的copula函数，确定联合分布函数；
[0079]
在步骤s5中，根据待修正地层界面的高程边际特性函数和联合分布函数，确定待修正地层界面至少一个坐标点对应的似然函数；
[0080]
在步骤s6中，根据至少一个坐标点对应的似然函数，对待修正地层界面的高程先验概率密度进行贝叶斯更新并进行归一化处理，确定至少一个坐标点的高程后验概率密度；
[0081]
在步骤s7中，根据至少一个坐标点的高程后验概率密度，确定对应的模型修正值，对待修正地层界面的已有模型中的高程值进行更新。
[0082]
在本发明实施例中，首先，进行观测样本数据和待修正地层界面的已有模型的有效获取；然后，利用观测样本数据的统计特性确定高程边际特性函数，以此表示待修正地层界面和辅助地层界面的分布特性；进而，对待修正地层界面的已有模型进行概率密度的估计，以此有效作为待修正地层界面的先验分布，便于后续求后验分布；进而，利用copula方法构建相邻多个地层界面高程变量的相关结构模型和联合分布函数，以copula函数对地层相关结构进行建模，构建多地层界面高程联合分布，准确描述多个相邻地层的复杂相关结构，不需要新增其他类型辅助数据，直接利用地层间的相关结构约束修正提高已有地层面模型的建模精度；接着，计算待修正地层界面每一个坐标点的似然函数，便于后续对修正模型的先验分布进行更新；然后，利用似然函数对已有模型中的每一个坐标点进行先验分布的更新，确定每一个坐标点对应的后验分布，并进行归一化，保证得到准确的高程后验概率密度；最后，根据每一个坐标点的高程后验概率密度求得对应的模型修正值，并对其对应的高程值进行更新，以此对已有模型进行准确快速的修正。
[0083]
优选地，结合图2来看，图2为本发明提供的确定高程边际特性函数的流程示意图，上述步骤s2包括步骤s21至步骤s25，其中：
[0084]
在步骤s21中，分别将待修正地层界面和辅助地层界面的高程值范围划分为多个等间距区间；
[0085]
在步骤s22中，分别统计待修正地层界面和辅助地层界面的观测样本中的高程值落入不同的等间距区间的频率数据；
[0086]
在步骤s23中，根据频率数据，分别确定辅助地层界面关于高程值的第一边际概率密度函数和待修正地层界面关于高程值的第二边际概率密度函数；
[0087]
在步骤s24中，根据频率数据，分别确定辅助地层界面的第一累积概率密度函数和待修正地层界面的第二累积概率密度函数；
[0088]
在步骤s25中，将第一累积概率密度函数变换为[0，1]区间内均匀分布下的分位数，确定辅助地层界面对应的第一分位数，将第二累积概率密度函数变换为[0，1]区间内均匀分布下的分位数，确定待修正地层界面对应的第二分位数。
[0089]
作为具体实施例，本实施例通过统计观测样本数据，利用频率数据，对待修正的地层界面和用于提供修正约束的辅助地层界面的高程边际分布进行有效的估计，以此求得相关的高程边际特性函数，便于后续确定联合分布函数。
[0090]
在本发明一个具体的实施例中，确定高程边际特性函数的具体流程如下：
[0091]
第一步：定义已有模型中(x,y)位置处，地层界面i(i为地层界面的编号，用于区分变量和函数对应的地层面)上的接触点的高程函数为z
i
(x,y)；可以理解的是，若第i个地层界面为待修正地层界面，则第i
‑
1个地层界面或第i+1个地层界面为辅助地层界面；
[0092]
第二步：统计所有观测位置上地层界面i的观测样本数据，将界面高程值z
i
的范围划分为多个等间距区间，计算样本在各区间出现的频率，根据频率数据拟合地层界面i的高程变量z
i
的概率密度函数f
i
(z
i
)(包括第一边际概率密度函数和第二边际概率密度函数)和累积概率密度函数f
i
(z
i
)(包括第一累积概率密度函数和第二累积概率密度函数)。其中，地层界面高程概率密度函数和累积概率密度函数可以使用经验概率密度函数和经验累积概率密度函数，也可采用但不限于核密度估计方法拟合频率数据得到；
[0093]
第三步：将地质界面i的高程样本数据值z
i
根据累积概率密度函数f
i
(z
i
)变换为[0,1]均匀分布下的分位数u
i
＝f
i
(z
i
)(第一分位数和第二分位数)。
[0094]
优选地，结合图3来看，图3为本发明提供的高程先验概率密度的流程示意图，上述步骤s3还包括步骤s31至步骤s35，其中：
[0095]
在步骤s31中，获取待修正地层界面的地层接触样本数据；
[0096]
在步骤s32中，通过克里金方法，确定至少一个坐标点对应的克里金方差；
[0097]
在步骤s33中，根据待修正地层界面的已有模型，确定至少一个坐标点对应的期望高程；
[0098]
在步骤s34中，将克里金方差作为正态分布的方差、期望高程作为正态分布的期望，建立正态分布函数，其中，正态分布函数作为待修正地层界面的高程值在至少一个坐标点下的估计概率密度函数；
[0099]
在步骤s35中，将待修正地层界面的估计概率密度函数作为对应的高程先验概率密度。
[0100]
作为具体实施例，本实施例利用地质模型不确定性评价方法，估计地质模型中待修正地层界面高程的概率分布，作为待修正模型的先验分布，采用了克里金方差、期望高程
作为对应的正态分布函数的参数，以此构建了准确有效的高程先验概率密度。
[0101]
在本发明一个具体的实施例中，确定高程边际特性函数的具体流程如下：
[0102]
第一步：利用不确定性评价方法分析待修正地层界面的已有模型的不确定性，给出每个(x,y)位置处，待修正地层界面的高程z的概率密度函数f(z)；
[0103]
第二步：利用待评价地层界面的地层接触样本数据，通过克里金方法(如简单克里金)计算每个(x,y)位置上的克里金方差将(x,y)位置处的已有模型的高程值z
m
(x,y)作为地层界面位置的最佳猜测。以z
m
(x,y)为期望，以方差建立正态分布函数作为待评价地层界面高程z在(x,y)位置上的概率分布，此方法得到的地层界面高程的概率密度函数为：
[0104][0105]
需要说明的是，除不确定性评价方法外，地层界面高程的概率密度函数还可以通过其他地质统计学方法或模型不确定性评价方法获得。
[0106]
第三步：将待修正地层界面的高程的概率密度函数f(z)作为待修正模型的高程先验概率密度。
[0107]
优选地，结合图4来看，图4为本发明提供的确定辅助地层数的流程示意图，上述步骤s4包括步骤s41至步骤s45，其中：
[0108]
在步骤s41中，分别获取第一辅助地层界面的已有模型和第二辅助地层界面的已有模型；
[0109]
在步骤s42中，通过模型不确定性评价方法，分别确定第一辅助地层界面的已有模型中的第一高程误差、第二辅助地层界面的已有模型中的第二高程误差；
[0110]
在步骤s43中，根据第一高程误差和第二高程误差，分别判断第一辅助地层界面的已有模型和第二辅助地层界面的已有模型是否满足数据精度条件；
[0111]
在步骤s44中，若其中之一满足，则将满足数据精度条件的第一辅助地层界面或第二辅助地层界面作为辅助地层界面，辅助地层数为一层；
[0112]
在步骤s45中，若都满足，则第一辅助地层界面和第二辅助地层界面都作为辅助地层界面，辅助地层数为两层。
[0113]
作为具体实施例，本发明实施例判断待修正地层界面的相邻界面的数据准确性，其上侧和下侧的相邻地层界面的数据精度是否较高，可用作模型修正的辅助数据。当只有一侧相邻地层界面高程精度较高时，可用的辅助地层数n＝1，当两侧相邻地层界面的高程精度均较高时，可用的辅助地层数n＝2，以此保证辅助地层界面的选取准确性。需要说明的是，本发明提供的模型不确定性评价方法以待评价位置上的已有地质模型取值作为地层面位置的最佳猜测，以该位置处计算的克里金估值方差表达地层面位置分布的随机性，主要步骤如下：
[0114]
第一步，根据待评价位置周边邻域点，利用克里金方法计算该位置上的克里金估值方差，以该方差作为该位置上地层面位置分布的方差
[0115]
第二步，以待评价位置上的已有地质模型取值作为地层面位置分布的期望em；
[0116]
第三步，构建以已有地质模型取值em为期望，以克里金估值方差为方差的高斯分布函数描述待评价位置处地层面模型的不确定性。
[0117]
优选地，结合图5来看，图5为本发明提供的当辅助地层数为一层时确定联合分布函数的流程示意图，上述步骤s4还包括步骤s46至步骤s47，其中：
[0118]
在步骤s46中，若辅助地层数为一层，则获取第一边际概率密度函数、第二边际概率密度函数、以及相邻地层界面的高程变量的copula函数；
[0119]
在步骤s47中，根据第一边际概率密度函数、第二边际概率密度函数、以及相邻地层界面的高程变量的copula函数的乘积，确定联合分布函数。
[0120]
其中，通过以下公式得到联合分布函数：
[0121]
f
a,b
(z
a
,z
b
)＝f
a
(z
a
)
·
f
b
(z
b
)
·
c
a,b
(u
a
,u
b
)
[0122]
其中，f
a,b
(z
a
,z
b
)表示辅助地层数为一层时的联合分布函数，f
a
(z
a
)表示第一边际概率密度函数，f
b
(z
b
)表示第二边际概率密度函数，c
a,b
(u
a
,u
b
)表示相邻地层界面的高程变量的copula函数，u
a
为第一分位数，u
b
为第二分位数。
[0123]
作为具体实施例，本发明实施例在若辅助地层数n＝1的情况下，选取数据精度良好的一侧相邻地层界面作为辅助地层界面，根据待修正地层界面和辅助地层界面的和双变量copula概率密度函数，以此有效确定两者之间的结构关联性。
[0124]
在本发明一个具体的实施例中，在辅助地层数为一层的情况下，联合分布函数的确定流程如下：
[0125]
第一步：若辅助地层数n＝1，令待修正地层界面为b，相邻的辅助地层界面为a；
[0126]
第二步：根据a、b界面的高程样本分位数u
a
,u
b
，估计与a、b界面高程样本拟合度最佳的双变量copula函数c
a,b
(u
a
,u
b
)函数族和copula参数；其中，参数估计方法可先采用但不限于极大似然估计法；
[0127]
第三步：根据a、b界面的边际概率密度函数f
a
(z
a
),f
b
(z
b
)和双变量copula概率密度函数c
a,b
(u
a
,u
b
)，计算a、b界面的联合概率密度函数：
[0128]
f
a,b
(z
a
,z
b
)＝f
a
(z
a
)
·
f
b
(z
b
)
·
c
a,b
(u
a
,u
b
)
[0129]
优选地，结合图6来看，图6为本发明提供的当辅助地层数为一层时确定似然函数的流程示意图，上述步骤s5包括步骤s51至步骤s52，其中：
[0130]
在步骤s51中，若辅助地层数为一层，则在第一辅助地层界面或第二辅助地层界面的已有模型中，确定至少一个坐标点对应的模型高程值；
[0131]
在步骤s52中，根据模型高程值、第二边际概率密度函数和联合分布函数，通过以下公式确定似然函数：
[0132][0133]
其中，l(z
b
|z
a
＝v
a
)、f
a|b
(z
a
＝v
a
|z
b
)表示待修正地层界面在至少一个坐标点处对应的似然函数，z
a
表示辅助地层界面对应的高程变量，z
b
表示待修正地层界面对应的高程变量，v
a
表示模型高程值，f
a，b
(z
a
＝v
a
，z
b
)表示辅助地层数为一层时的联合分布函数，f
b
(z
b
)表示第二边际概率密度函数。
[0134]
作为具体实施例，本发明实施例在若辅助地层数n＝1的情况下，根据待修正地层界面的高程边际分布和多地层界面高程的联合分布，计算每一个坐标点(x,y)位置上，对应的似然函数，以便后续对每一坐标点对应的高程值进行更新。
[0135]
在本发明一个具体的实施例中，在辅助地层数为一层的情况下，似然函数的确定流程如下：
[0136]
若辅助地层数n＝1，根据坐标点s
j
＝(x
j
,y
j
)处的辅助地层a的数据值z
a
＝v
a
，待修正地层界面b的边际概率密度函数f
b
(z
b
)和a、b界面的联合概率密度函数f
a,b
(z
a
,z
b
)，计算待修正地层界面b在s0处的似然函数如下：
[0137][0138]
优选地，结合图7来看，图7为本发明提供的当辅助地层数为两层时确定联合分布函数的流程示意图，上述步骤s4还包括步骤s48至步骤s49，其中：
[0139]
在步骤s48中，若辅助地层数为两层，则获取第二边际概率密度函数、第一辅助地层界面的边际概率密度函数、第二辅助地层界面的边际概率密度函数、以及相邻地层界面的高程变量的copula函数；
[0140]
在步骤s49中，根据第二边际概率密度函数、第一辅助地层界面的边际概率密度函数、第二辅助地层界面的边际概率密度函数、以及相邻地层界面的高程变量的copula函数的乘积，确定联合分布函数。
[0141]
其中，通过以下公式得到联合分布函数：
[0142]
f
a,b,c
(z
a
,z
b
,z
c
)＝f
a
(z
a
)
·
f
b
(z
b
)
·
f
c
(z
c
)
·
c
a,b
(u
a
,u
b
)
·
c
b,c
(u
b
,u
c
)
·
c
a,c|b
(u
a|b
,u
c|b
)
[0143]
其中，f
a，b，c
(z
a
，z
b
，z
c
)表示辅助地层数为两层时的联合分布函数，f
b
(z
b
)表示第二边际概率密度函数，f
a
(z
a
)表示第一辅助地层界面的边际概率密度函数，f
c
(z
c
)表示第二辅助地层界面的边际概率密度函数，c
a，b
(u
a
，u
b
)、c
b，c
(u
b
，u
c
)以及c
a，c|b
(u
a|b
，u
c|b
)表示相邻地层界面的高程变量的copula函数，u
a
为第一辅助地层界面对应的第一分位数，v
b
为第二辅助地层界面对应的第一分位数，u
c
为第二分位数。
[0144]
作为具体实施例，本发明实施例在若辅助地层数n＝2的情况下，选取数据精度良好的两侧相邻地层界面作为辅助地层界面，根据待修正地层界面和辅助地层界面的和双变量copula概率密度函数，以此有效确定三者之间的结构关联性。
[0145]
在本发明一个具体的实施例中，在辅助地层数为两层的情况下，联合分布函数的确定流程如下：
[0146]
第一步：若辅助地层数n＝2，令待修正地层界面为b，修正地层上下两侧的辅助地层界面为a、c。根据a、b、c界面的高程样本分位数u
a
,u
b
,u
c
，分别估计a、b界面高程的双变量copula函数c
a,b
(u
a
,u
b
)、b、c界面高程的双变量copula函数c
b,c
(u
b
,u
c
)和a、c界面高程的条件copula函数c
a,c|b
(u
a|b
,u
c|b
)的最佳拟合函数族和对应的copula参数。其中，参数估计方法可先采用但不限于极大似然估计法。
[0147]
第二步：根据a、b、c界面的边际概率密度函数f
a
(z
a
),f
b
(z
b
),f
c
(z
c
)和双变量copula概率密度函数c
a,b
(u
a
,u
b
),c
b,c
(u
b
,u
c
),c
a,c|b
(u
a|b
,u
c|b
)，利用pair
‑
copula方法计算a、b、c界面的联合概率密度函数：
[0148]
f
a,b,c
(z
a
,z
b
,z
c
)＝f
a
(z
a
)
·
f
b
(z
b
)
·
f
c
(z
c
)
·
c
a,b
(u
a
,u
b
)
·
c
b,c
(u
b
,u
c
)
·
c
a,c|b
(u
a|b
,u
c|b
)
[0149]
优选地，结合图8来看，图8为本发明提供的当辅助地层数为两层时确定似然函数的流程示意图，上述步骤s5包括步骤s53至步骤s54，其中：
[0150]
在步骤s53中，若辅助地层数为两层，则在第一辅助地层界面和第二辅助地层界面的已有模型中，确定至少一个坐标点对应的模型高程值；
[0151]
在步骤s54中，根据模型高程值、第二边际概率密度函数和联合分布函数，通过以下公式确定似然函数：
[0152][0153]
其中，l(z
b
|z
a
＝v
a
，z
c
＝v
c
)、f
a，c|b
(z
a
＝v
a
，z
c
＝v
c
|z
b
)表示待修正地层界面在至少一个坐标点处对应的似然函数，z
a
表示第一辅助地层界面对应的高程变量，z
b
表示待修正地层界面对应的高程变量，z
c
表示第二辅助地层界面对应的高程变量，v
a
表示第一辅助地层界面对应的模型高程值，v
c
表示第二辅助地层界面对应的模型高程值，f
a，b，c
(z
a
，z
b
，z
c
)表示辅助地层数为两层时的联合分布函数，f
b
(z
b
)表示第二边际概率密度函数。
[0154]
作为具体实施例，本发明实施例在若辅助地层数n＝2的情况下，根据待修正地层界面的高程边际分布和多地层界面高程的联合分布，计算每一个坐标点(x,y)位置上，对应的似然函数，以便后续对每一坐标点对应的高程值进行更新。
[0155]
在本发明一个具体的实施例中，在辅助地层数为两层的情况下，似然函数的确定流程如下：
[0156]
若辅助地层数n＝2，根据坐标点s
j
＝(x
j
,y
j
)处的辅助地层a、c的数据值z
a
＝v
a
,z
c
＝v
c
，待修正地层界面b的边际概率密度函数f
b
(z
b
)和a、b界面的联合概率密度函数f
a,b,c
(z
a
,z
b
,z
c
)，计算待修正地层界面b在坐标点s
j
＝(x
j
,y
j
)处的似然函数：
[0157][0158]
优选地，结合图9来看，图9为本发明提供的更新高程后验概率密度的流程示意图，上述步骤s6包括步骤s61至步骤s62，其中：
[0159]
在步骤s61中，在至少一个坐标点处，根据似然函数和高程先验概率密度的乘积，对高程先验概率密度进行贝叶斯更新；
[0160]
在步骤s62中，对更新后的高程先验概率密度进行归一化处理。
[0161]
作为具体实施例，本发明实施例对待修正模型的先验分布进行贝叶斯更新，并做归一化处理，得到待修正地层界面高程的后验分布，保证后验分布的准确性。
[0162]
在本发明一个具体的实施例中，通过下式对坐标点s
j
＝(x
j
,y
j
)处的待修正界面b的高程先验概率密度f(z
b
)进行贝叶斯更新：
[0163]
f(z
b
|z
a
＝v
a
,k)
∝
f(z
b
)
·
l(z
b
|z
a
＝v
a
,k).
[0164]
式中，
∝
为正比于的符号，l(z
b
|z
a
＝v
a
,k)为待修正地层界面b的高程在相邻地层数据约束下的似然函数，若辅助地层数n＝1，l(z
b
|z
a
＝v
a
,k)＝l(z
b
|z
a
＝v
a
)；若辅助地层数n＝2，l(z
b
|z
a
＝v
a
,k)＝l(z
b
|z
a
＝v
a
,z
c
＝v
c
)。f(z
b
|z
a
＝v
a
,k)为归一化前的待修正地层界面高程后验概率密度函数，其中：
[0165]
若辅助地层数n＝1，f(z
b
|z
a
＝v
a
,k)＝f(z
b
|z
a
＝v
a
)
[0166]
若辅助地层数n＝2，f(z
b
|z
a
＝v
a
,k)＝f(z
b
|z
a
＝v
a
,z
c
＝v
c
)。
[0167]
对f(z
b
|z
a
＝v
a
,k)做归一化处理，令计算尺度因子α，得到归一化的b地层界面的后验概率分布：
[0168]
f
b|a,k
(z
b
|z
a
＝v
a
,k)＝α
·
f(z
b
|z
a
＝v
a
,k)
[0169]
优选地，结合图10来看，图10为本发明提供的更新已有模型的流程示意图，上述步骤s7包括步骤s71至步骤s73，其中：
[0170]
在步骤s71中，根据至少一个坐标点的高程后验概率密度，确定高程后验概率密度对应的条件期望值；
[0171]
在步骤s72中，将条件期望值作为至少一个坐标点在待修正地层界面的已有模型中对应的模型修正值；
[0172]
在步骤s73中，根据模型修正值，对至少一个坐标点在待修正地层界面的已有模型中的高程值进行更新。
[0173]
作为具体实施例，本发明实施例将条件期望值作为待修正地层界面b在s
j
＝(x
j
,y
j
)位置上的模型修正值，以此完成整个已有模型的有效更新，提供精准度。
[0174]
在本发明一个具体的实施例中，根据辅助地层界面的数据计算后验分布的条件期望值或者中位数作为模型修正值。对于对称的分布f
b|a,k
(z
b
|z
a
＝v
a
,k)，通常计算分布的条件期望值作为模型修正值。计算后验分布f
b|a,k
(z
b
|z
a
＝v
a
,k)的条件期望值具体公式为：
[0175][0176]
将条件期望值e
b|a,k
作为待修正地层界面b在s
j
＝(x
j
,y
j
)位置上的模型修正值。对于非对称的分布f
b|a,k
(z
b
|z
a
＝v
a
,k)，也可选取后验分布的中位数作为模型修正值；对每一个待修正的(x,y)位置，重复以上，完成整个已有模型的修正。
[0177]
在本发明一个具体的实施例中，结合图11来看，图11为本发明实施例提供的地质钻孔数据示意图，图11中以上海城市地质钻孔数据为例，采用自上而下的三个地层qh3,qh2,qh1的下界面作为实验对象，利用相邻上部地层的地层界面qh3和相邻下部地层的地层界面qh1的钻孔采样数据去修正精度较低的界面qh2的模型，为了在实验中便于标注，以1、2、3号地层界面分别指代qh3,qh2,qh1的下界面。将1，2地层和2，3地层设置为两组相邻地层，分别进行了两个实验：
[0178]
第一个实验：结合图12来看，图12为本发明提供的修正示意图一，利用上侧的1号地层接触样本数据修正下方的2号地层界面模型；
[0179]
求得1号地层与2号地层数据拟合copula函数：survival tawn copula(δ＝1,ρ＝0.53,θ＝3.57,kendallτ＝0.43)。
[0180][0181]
θ∈[1,+∞),δ∈[0,1],ρ∈[0,1].
[0182]
下面以图11中编号为29钻孔位置为例展示利用1号地层界面数据修正2号地层界面模型。1号地层界面样本高程值v1＝
‑
4.27m，2号地层界面样本高程值v2＝
‑
20.27m，克里金
方法得到的2号地层界面模型高程值v
k
＝
‑
17.37m，克里金估计方差编号为29位置上2号地层界面模型先验分布f
d29
(z2):n(
‑
17.37,9.80)；图12展示了编号为29钻孔位置模型修正的结果：2号地层界面编号为29位置的先验分布f
d29
(z2)和后验分布f
d29
(z2|z1＝v1)。将修正前的模型值v
k
＝
‑
17.37m、修正后的模型值e(z2|z1＝v1)＝
‑
19.61m与样本真实高程值(灰线)v2＝
‑
20.27m对比可知，修正后的模型更接近真实值。
[0183]
第二个实验：结合图13来看，图13为本发明提供的修正示意图二，利用1号和3号地层接触样本数据修正2号地层界面模型；
[0184]
1号地层与2号地层数据拟合copula函数：survival tawn copula(δ＝0.53,ρ＝1,θ＝3.99,kendallτ＝0.44)：
[0185][0186]
θ∈[1,+∞),δ∈[0,1],ρ∈[0,1].
[0187]
2号地层与3号地层数据拟合copula函数：tawn copula(δ＝1,ρ＝0.69,θ＝5.44,kendallτ＝0.59)。
[0188][0189]
1号地层与3号地层数据拟合条件copula函数：frank copula(θ＝1.68,kendallτ＝0.18)：
[0190][0191]
式中，δ,ρ,θ均为copula函数参数。
[0192]
下面以图11中编号为44钻孔位置为例展示利用1号地层界面数据修正2号地层界面模型。编号为44位置上1号地层界面样本高程值v1＝
‑
4.10m，3号地层界面样本高程值v3＝
‑
30.40m，2号地层界面样本高程值v2＝
‑
21.10m，克里金方法得到的2号地层界面模型高程值v
k
＝
‑
16.56m，克里金估计方差编号为44位置上2号地层界面模型先验分布f
d44
(z2):n(
‑
16.56,9.99)。图13展示了编号为44钻孔位置模型修正的结果：2号地层界面编号为44位置的先验分布f
d44
(z2)和后验分布f
d44
(z2|z1＝v1,z3＝v3)。将修正前的模型值v
k
＝
‑
16.56m、修正后的模型值e(z2|z1＝v1,z3＝v3)＝
‑
19.93m与样本真实高程值v2＝
‑
21.10m对比可知，修正后的模型更接近真实值。
[0193]
实施例2
[0194]
本发明实施例提供了一种利用相邻地层互相关约束的地层面模型修正装置，包括处理器以及存储器，存储器上存储有计算机程序，计算机程序被处理器执行时，实现如上所述的利用相邻地层互相关约束的地层面模型修正方法。
[0195]
实施例3
[0196]
本发明实施例提供了一种计算机可读存储介质，其上存储有计算机程序，计算机该程序被处理器执行时，实现如上所述的利用相邻地层互相关约束的地层面模型修正方法。
[0197]
本发明公开了一种利用相邻地层互相关约束的地层面模型修正方法，首先，进行
观测样本数据和待修正地层界面的已有模型的有效获取；然后，利用观测样本数据的统计特性确定高程边际特性函数，以此表示待修正地层界面和辅助地层界面的分布特性；进而，对待修正地层界面的已有模型进行概率密度的估计，以此有效作为待修正地层界面的先验分布，便于后续求后验分布；进而，利用copula方法构建相邻多个地层界面高程变量的相关结构模型和联合分布函数，以copula函数对地层相关结构进行建模，构建多地层界面高程联合分布，准确描述多个相邻地层的复杂相关结构，不需要新增其他类型辅助数据，直接利用地层间的相关结构约束修正提高已有地层面模型的建模精度；接着，计算待修正地层界面每一个坐标点的似然函数，便于后续对修正模型的先验分布进行更新；然后，利用似然函数对已有模型中的每一个坐标点进行先验分布的更新，确定每一个坐标点对应的后验分布，并进行归一化，保证得到准确的高程后验概率密度；最后，根据每一个坐标点的高程后验概率密度求得对应的模型修正值，并对其对应的高程值进行更新，以此对已有模型进行准确快速的修正。
[0198]
本发明技术方案，不需要额外采集辅助数据，也不需要人工交互编辑，利用相邻地层互相关约束的地层面模型修正方法及系统，将相邻地层的互相关结构作为一种约束条件引入模型修正，通过copula方法描述相邻地层的复杂相关结构，通过相关性进行模型修正，提高已有模型的建模精度，大大提高了算法的实用性和快速性。
[0199]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。