本发明涉及三维地质建模及地质统计模拟领域，尤其涉及一种基于局部搜索策略的多点统计三维地质模型自动重构方法。
背景技术：
三维地质建模及可视化技术是地球空间信息科学领域的核心技术之一，也是三维地理、地质信息系统建设的关键和基础技术。地质结构的三维可视化表征可以提供对地下地质现象及构造更加真实、直观的描述，已经在地下水、石油和固体矿产资源的定量评价中发挥了重要作用。储层的空间异质性和连通性对地下流体运聚起到了重要的控制作用，这对量化和预测地下资源的形成和分布是至关重要的。然而，对于一个实际的三维应用来说，这些属性是很难描述和模型化的，因为我们可以获得的已知数据是非常稀疏的。传统两点地质统计方法已经被用来重构地下沉积结构和相分布的三维模型。然而，两点地质统计学方法难以捕获高阶统计信息，因此各向异性特征和连通性模式不能被很好地再现。基于对象的方法也由于无法参数化复杂地质结构和现象而难以胜任非均质岩相的三维重构。多点统计(multiplepointstatistics,简称mps)方法是三维地质模型自动构建与模拟领域的一个重要分支，其基本理论的发展以及各种算法的提出为各向异性的地下复杂地质结构的三维重建提供了技术保障，并且已经在储层模拟、水文地质建模、多孔介质重构等多个地学领域中取得了很好的应用。
多点统计学方法旨在通过直接从训练图像中提取空间模式而非利用变差函数来描述空间异构几何特征。训练图像是从实际现象中抽象得到的概念模型，在基于多点统计的随机模拟中起着至关重要的作用。多点地质统计的理论自从提出以来，已经发展出了多种基于该理论的随机模拟和建模方法。在这些方法中，使用固定的搜索模板扫描训练图像，并且将多点统计模式存储在树或链表式数据结构中。然后根据当前数据事件计算条件概率密度函数(conditionalprobabilitydensityfunction，简称cpdf)，依次获得当前位置的模拟结果。但是，一个值得注意的矛盾就在于数据模板的大小和大尺度结构的重现能力上。数据模板越大，数据事件的尺寸也越大，因此在训练图像上捕获该数据事件的重复次数将越少，这将导致难以推断出相应的条件概率密度函数。反之，当数据模板太小时，训练图像的大尺度结构将不会被再现。另外，搜索模板太大会增加计算成本和内存消耗。多重网格概念减轻了上述限制，但是由于数据模板固定，上述问题并没有被完全消除。作为一种近似表达，模式距离可用来表示训练图像中的模式和数据事件之间的匹配程度。对于基于概率密度函数的多点统计方法，模式距离的使用会在一定程度上增加匹配成功的几率、减少模式存储的数量，但是一些复杂的细节特征也会由于近似表达而被过滤掉。
合适、可靠的训练图像是基于多点统计的随机模拟方法的基础。对于实际的三维应用，三维训练图像的获取十分困难，而且其可靠性难以评估。为了克服这个困难，一些方法尝试直接使用低维数据(例如钻孔、剖面，露头、遥感影像等)来重建地下复杂地质结构的三维模型，而不是使用完整的三维训练图像。这类方法主要包括两种。其中一种广泛应用的方法directsampling(简称ds)直接以不完整的已知数据集为参考，从中获取多点统计模式，依次完成对整个三维建模空间的重构。这种方法的主要优点是在数据事件和已知数据之间的模式匹配中考虑了各个方向上的空间异质性。然而，精细的大范围三维模型将包含数千万个模拟网格。要想从稀疏的已知数据中找到相匹配的多点模式，每一个网格都将需要上百万次的尝试才能获得最终结果。因此，对于精细的三维应用，这种对全局数据进行搜索匹配的方法仍然受到计算效率的困扰。另一种方法则使用二维地质剖面作为训练图像，通过执行一系列的二维模拟来填充整个三维建模空间。这种方法虽然效率较高，但是其模拟结果的空间连续性难以得到保证。
从另一个角度来看，借助稀疏、分散的一维钻孔等抽样数据、适当的模拟方法以及专家知识经验来获得二维地质资料(如地质平面图、剖面图等)是地质工作人员的普遍工作流程。那么，在低维数据的基础上，使用可靠的模拟方法来重建三维地质结构模型是符合逻辑的，也将是可行的。但在实际建模实践中，由于地质现象和过程的空间异质性，直接使用真实地质剖面、遥感图像等作为训练图像进行复杂地质结构的三维重构仍然受到其非平稳性特征的限制。
概率融合是多元统计分析的基础。一个地质现象往往受到来自不同方面或者不同角度的因素的影响，这就需要利用概率融合方法获得一个最终的联合概率分布。根据其数学性质，概率融合方法可分为基于加法的方法和基于乘法的方法。linearpooling公式是一种广泛使用的基于概率加法的方法，其主要特点是灵活和简单。乘法方法包括bordley公式和log-linearpooling公式。已经证实，对于只有两个属性的情况，这两个公式是等效的，但后者对于非二元事件具有更好的性能。线性概率融合方法就像逻辑运算“or”一样，用来联合几个独立的概率分布以获得一个更大、更稳定的联合概率分布。而基于乘法的方法更像是逻辑运算“and”，通常使用这种方法来融合具有显著相关性的概率分布。
综上所述，由于受到缺少完整、可靠的三维训练图像的限制，传统的基于多点统计的随机模拟方法很难在实际三维地质结构的模拟与重构中得到很好的应用。因此，利用已有的低维数据来重构地下地质结构的三维模型是一条可行且符合地质数据处理流程的路径。已有的对所有已知数据进行全局搜索匹配的多点统计重构方法也受到了其计算效率的限制。而通过使用多点统计技术模拟二维切片来逐步填充三维空间域的方法难以保证全空间不同方向上重构结果的空间连续性。鉴于上述分析的多点统计方法在复杂地质结构三维重构中的问题和挑战，有必要提出一种基于低维数据的高效的多点统计三维地质结构模型的自动重构方法。
技术实现要素：
有鉴于此，本发明的实施例提供了一种高效且能够刻画空间异质性特征的基于局部搜索策略的多点统计三维地质模型自动重构方法，具体包括以下过程：
(1)加载数据，并将已知数据分配到模拟网格中；
(2)根据已知的各个二维横纵地质剖面在模拟网格中的位置，获取其在x、y、z方向上的索引，并计算被横纵剖面分割成的所有子区域的先验概率p0；
(3)确定多重网格g，依次对每一重网格进行处理，直至所有多重网格都被处理；获取当前网格g，如果此网格有效，则转步骤(4)，否则转步骤(13)；
(4)为当前网格g定义一个随机模拟路径，使得覆盖当前网格g上的所有未知结点；
(5)按顺序从随机模拟路径中获得当前模拟位置x，如果当前结点存在，转步骤(6)，否则转步骤(12)；
(6)根据当前结点在模拟网格中的位置，确定x、y、z三个方向上对应的搜索邻域nx、ny、nz和对应的数据事件devx、devy、devz；
(7)获取围绕着当前模拟位置x所在子区域的所有n个子截面的索引{x0,x1}、{y0,y1}、{z0,z1}，且1≤n≤6；
(8)以步骤(7)得到的子截面为训练图像，依次用数据事件devx、devy、devz对相应方向上的子截面进行扫描，获取对应的n个条件概率分布函数{cpdf1,...,cpdfn}；
(9)从所有子截面中获得了对应于当前模拟位置x的n个条件概率密度函数{cpdf1,...,cpdfn}，然后从步骤(2)得到的先验概率p0中获取当前模拟位置x对应的子区域的先验概率p0；
(10)根据概率融合策略，融合获得的n个条件概率密度函数{cpdf1,...,cpdfn}和先验概率p0，从而获得联合概率分布；
(11)从最终的联合概率分布中随机提取一个对应属性，赋值给当前模拟位置x，转步骤(5)继续其它未知位置的属性模拟；
(12)当前网格g上所有未知结点已完成模拟，转步骤(3)进行下一重网格的模拟；
(13)各个多重网格g上的未知位置全部模拟完成，保存模拟结果，结束本次重构；
整个模拟网格上的所有结点的属性值已经全部获得，从而得到了充满整个建模空间的三维地质结构模型。
进一步，所述步骤(8)中，对于其中任意一个子截面的扫描过程，执行以下步骤(8-1)至(8-8)：
(8-1)根据子截面在整个模拟网格中的索引，获取当前被扫描子截面的所有结点；
(8-2)设置一个随机扫描路径path，并初始化已匹配模式数目的计数器sum＝0；
(8-3)依次从随机扫描路径中获取子截面中的对应位置pathi，如果i＜path.size()×f成立，则转步骤(8-4)，否则转步骤(9)；其中f是给定的被扫描子截面的比例，即允许子截面的一部分被扫描而非全部；
(8-4)根据pathi所在训练图像中的位置获取对应的邻域n'和数据事件dev'；
(8-5)使用以下公式计算来自模拟网格和来自训练图像的两个数据事件之间的模式距离d{dev,dev'}：
其中k表示数据事件中包含结点的个数，z(xj)和z(yj)分别表示来自于模拟网格和训练图像中数据事件对应中心结点的属性值；
(8-6)如果d{dev,dev'}≤t成立，则转步骤(8-7)，否则转步骤(8-8)；其中t为当前网格g对应的模式距离的阈值，如果当前模式距离小于该阈值，则被接收；
(8-7)根据训练图像中数据事件中心结点对应的属性值，更新当前子截面对应的条件概率密度函数cpdf，并且使sum++；
(8-8)如果sum>maximumn成立，则转步骤(9)，否则转步骤(8-3)继续对当前子截面进行扫描；其中maximumn表示从当前训练图像中获取多点模式的最大数目。
进一步，所述步骤(8)中，在对子截面进行扫描时，多重网格的搜索邻域半径r和模式距离的阈值t允许随着网格级别的增加而线性递减；第一重网格上的邻域半径r和阈值t由用户给定，即r＝r0，t＝t0；最后一重网格上将其固定为r＝1，t＝0；对于中间任意一重网格，其邻域半径r和阈值t可表示为：
其中m表示多重网格的总数，τ(1≤τ≤m)表示当前网格的索引。
进一步，所述步骤(8-8)中，maximumn的一个推荐取值范围是[50,100]，当然也要根据模型网格的尺度以及地质结构的复杂程度做适度调整；从而通过控制匹配成功的模式的数目来控制条件概率密度函数cpdf的稳定性，避免不必要的计算耗时。
进一步，所述步骤(10)中，概率融合策略包含以下两步：
(10-1)使用基于概率加法的linearpooling公式融合来自同一方向上互相平行的两个子截面的条件概率密度函数：
其中a表示概率融合的属性类别，pg(a)表示融合后的该属性的概率，pi(a)表示来自任一子截面的属性概率，wi为对应权重；一个子区域平行方向最多有两个截面，也只有当存在两个截面时才需要进行概率融合；此时m＝2，权重w1、w2的值由以下公式求取：
其中d1、d2分别表示当前模拟位置距两个子截面的距离；
(10-2)使用基于概率乘法的log-linearpooling公式融合经步骤(10-1)处理后的三个正交方向的概率分布以及先验概率：
其中a表示概率融合的属性类别，n为参与概率融合的概率分布的数目，1≤n≤3，即最多为步骤(10-1)处理后得到的x、y、z三个正交方向的概率分布，pg(a)表示融合后的该属性的概率，pi(a)表示经步骤(10-1)处理后任一子方向的属性概率，p0(a)为该属性的先验概率，wi为各个方向上概率的权重；包括先验概率的权重w0在内，其和必须为1，即
与现有技术相比，本发明具有以下有益效果：
(1)本发明以低维数据(如二维剖面、一维的钻孔等取样数据)为基础，以横纵剖面代替完整的三维训练图像，实现了对复杂三维地质结构模型的自动重构，为从低维地质数据到三维地质模型的重构提供了新思路，同时也为基于多点统计的三维地质模拟和建模提供了新技术、新方法。
(2)本发明提供的局部搜索策略，将多点统计信息的获取固定到距离当前模拟位置最近的局部子截面，这将避免距离太远的多点模式被捕获，因此保证了其模拟结果的可靠性；而且局部子截面的搜索策略降低了对训练图像的搜索次数，从而提升了模拟效率。
(3)本发明中所采用多重网格以及搜索邻域的参数调整策略，允许搜索邻域的半径和模式距离的阈值随着网格级别的增加而递减；在最初的多重网格上，其邻域半径较大，模式距离的阈值也较大，这保证了大尺度结构模式的获取；而随着邻域半径和阈值的减小，模式匹配将趋向于完全匹配，这将保证细小的结构特征的重构；该策略在一定程度上减轻了传统基于多点统计的随机模拟方法中搜索模板固定且尺寸不宜太大的问题。
(4)本发明中的概率融合策略和权重选择方法，使得最终的联合概率分布包含了来自各个方向上子截面的概率分布；其中平行子截面的基于加法的概率融合拓展了概率分布的容量，使得融合后的联合概率更加稳定；正交方向上的基于乘法的概率融合保证了不同方向上多点模式的最终刻画，使得其模拟结果在不同方向上都表现出了很好的空间连续性。
(5)本发明中对最大匹配模式数目的限制，既保证了稳定的条件概率密度函数的获取，又避免了无效的多点统计搜索与模式匹配，这也在一定程度上提升了三维模型重构的效率。
(6)本发明可在各种三维地质信息系统、地理信息系统、地质建模与模拟系统等软件中推广应用。
附图说明
图1是本发明基于局部搜索策略的多点统计三维地质模型自动重构方法的总体流程。
图2是本发明实施例中局部搜索策略的示意图。
图3是本发明实施例中相邻子截面上搜索邻域、数据事件及模式匹配过程的示意图。
图4是本发明实施例中多重网格与搜索邻域之间的关系示意图。
图5是本发明实施例中所采用的三维参考模型、交错剖面以及两种不同方法的两个重构结果：ds方法和本发明方法。
图6是本发明实施例中建模结果的统计特性对比，其中每一种方法都包含20个重构结果：图6(a)为两种方法的重构结果在x、y、z三个方向上的空间连通性曲线；图6(b)为两种方法的重构结果的整体变差函数曲线；图6(c)为随着网格尺度的增加，两种方法计算效率的对比曲线。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步的描述。
本发明的实施例提供了基于局部搜索策略的多点统计三维地质模型自动重构方法，直接利用二维的横纵地质剖面重建三维地质结构模型，并给出了一种基于多点统计的局部搜索策略；该策略根据已知剖面的空间交错关系，将三维建模空间划分成若干个子区域，每一个子区域都被n(1≤n≤6)个子截面围绕；对于分布于每一个子区域内的未知属性位置，其多点统计信息直接从包围该子区域的n个子截面中获取：首先根据其位置获取周围的n个子截面，然后根据设定的邻域参数和已知数据点定制x、y、z三个方向上的数据事件，用对应的数据事件扫描对应方向上的子截面，获得n个条件概率密度函数；然后使用合理的概率融合策略获得联合概率密度分布，最终得到该位置的属性值。
请参照图1，该方法包括以下过程：
(1)加载数据，将已知数据点(离散后的横纵地质剖面、遥感影像或者钻孔取样数据)分配到模拟网格中；
(2)根据已知的各个二维横纵地质剖面在模拟网格中的位置，获取各个剖面在模拟网格中x、y、z方向上的索引，并计算被横纵剖面分割成的所有子区域的先验概率p0；
先验概率p0是已知数据点的属性的概率分布，由于本发明中的局部搜索策略只是求搜索围绕当前待模拟位置的若干个子截面，所以每一个子区域的先验概率也将从围绕着该子区域的这些子截面中获取。
(3)确定多重网格g，依次对每一重网格进行处理，直至所有多重网格都被处理；获取当前网格g，如果此网格有效，则转步骤(4)，否则转步骤(13)；
(4)为当前网格g定义一个随机模拟路径，使得覆盖当前网格g上的所有未知结点；
(5)按顺序从随机模拟路径中获得当前模拟位置x，如果当前结点存在，转步骤(6)，否则转步骤(12)；
(6)根据当前结点在模拟网格中的位置，确定x、y、z三个方向上对应的搜索邻域nx、ny、nz和对应的数据事件devx、devy、devz；
(7)获取围绕着当前模拟位置x所在子区域的所有n(1≤n≤6)个子截面的索引{x0,x1}、{y0,y1}、{z0,z1}；
如图2所示，一般而言每一个子区域由6个子截面所包围，但有时子截面会小于6；同时也允许某一个或者两个方向上没有剖面，但是至少一个方向上必须有；
(8)以这些子截面为训练图像，依次用数据事件devx、devy、devz对相应方向上的子截面进行扫描，获取对应的n个条件概率分布函数{cpdf1,...,cpdfn}；
如图3所示为一个子截面上搜索邻域以及数据事件的示意图，而且对于中心结点而言，其搜索窗口可以布满整个子截面，而数据事件内其它结点允许其从相邻子截面中抓取，这将保证在剖面分割附近模拟结果的空间连续性；如图4所示为一个3重网格与搜索邻域之间的关系示意图。
对于其中任意一个子截面的扫描过程，执行以下步骤(8-1)至(8-8)：
(8-1)根据子截面在模拟网格中的索引，获取当前被扫描子截面的所有结点；
(8-2)设置一个随机扫描路径path，并初始化已匹配模式数目的计数器sum＝0；
(8-3)依次从扫描路径中获取子截面中的对应位置pathi，如果i＜path.size()×f成立，则转步骤(8-4)，否则转步骤(9)；其中f是给定的被扫描截面的比例，即允许子截面的一部分被扫描而非全部；
(8-4)根据pathi所在训练图像中的位置获取对应的邻域n'和数据事件dev'；
(8-5)使用以下公式计算来自模拟网格和来自训练图像的两个数据事件之间的模式距离d{dev,dev'}：
其中k表示数据事件中包含结点的个数，z(xj)和z(yj)分别表示来自于模拟网格和训练图像中数据事件对应中心结点的属性值；
(8-6)如果d{dev,dev'}≤t成立(t为当前网格g对应的模式距离的阈值)，则转步骤(8-7)，否则转步骤(8-8)；
多重网格的搜索邻域半径r和模式距离的阈值t允许随着网格级别的增加而线性递减；第一重网格上的邻域半径r和阈值t由用户给定，即r＝r0，t＝t0；，在本发明给出的实施实例(图5所示)中，给定初始邻域半径r0＝30，初始阈值t0＝0.1，最后一重网格上将其固定为r＝1，t＝0；中间任意一重网格，其邻域半径r和阈值t可表示为：
其中m表示多重网格的总数，τ(1≤τ≤m)表示当前网格的索引。
(8-7)根据训练图像中数据事件中心结点对应的属性值，更新当前子截面对应的条件概率密度函数cpdf，并且使sum++；
(8-8)如果sum>maximumn成立，则转步骤(9)，否则转步骤(8-3)继续对当前子截面进行扫描；其中maximumn表示从当前子截面中获取多点模式的最大数目，图5中设定maximumn＝80；maximumn通过控制匹配成功的模式的数目来控制条件概率密度函数cpdf的稳定性；如果maximumn太小，将导致cpdf容量太小，各概率分布不稳定，从而增加模拟结果的不确定性；随着maximumn的增大，cpdf将趋于稳定，所以太大的maximumn将会带来毫无意义的搜索与匹配，增加了计算耗时；一个推荐取值范围是[50,100]，当然也要根据模型网格的尺度以及地质结构的复杂程度来做适度调整；
(9)从p0中获取当前模拟位置x对应的子区域的先验概率p0；
(10)根据概率融合策略，融合获得的n个条件概率密度函数{cpdf1,...,cpdfn}和先验概率p0，从而获得联合概率分布；
其概率融合策略包括以下两步：
(10-1)使用基于概率加法的linearpooling公式融合来自同一方向上互相平行的两个子截面的条件概率密度函数：
其中a表示概率融合的属性类别，pg(a)表示融合后的该属性的概率，pi(a)表示来自任一子截面的属性概率，wi为对应权重；在本发明中，对于一个子区域平行方向最多有两个截面，也只有当存在两个截面时才需要进行概率融合；此时m＝2，权重w1、w2的值由以下公式求取：
其中d1、d2分别表示当前模拟位置距两个子截面的距离；
(10-2)使用基于概率乘法的log-linearpooling公式融合经第一步处理后的三个正交方向的概率分布以及先验概率：
其中a表示概率融合的属性类别，n为参与概率融合的概率分布的数目，pg(a)表示融合后的该属性的概率，pi(a)表示经第一步处理后任一子方向的属性概率，p0(a)为该属性的先验概率，wi为各个方向上概率的权重；wi的选择由对应子截面的复杂程度(重要程度)来决定，如果某一向上的剖面结构更加复杂，能够提供更多更丰富的多点模式，则设置相对较大的对应权重，但必须保证
(11)从最终的联合概率分布中随机提取一个对应属性，赋值给当前模拟位置x；转步骤(5)继续其它未知位置的属性模拟；
(12)当前网格g上所有未知结点已完成模拟，转步骤(3)进行下一重网格的模拟；
(13)各个多重网格上的未知位置全部模拟完成，保存模拟结果，结束本次重构。
至此，整个模拟网格上的所有结点的属性值已经全部获得，从而得到了充满整个建模空间的三维地质结构模型。
如图5所示为本发明为了说明其方法优越性所实施的一个实验案例，如图分别为一个三维参考模型、所使用的交错剖面(x、y、z每个方向上各两个剖面)以及两种不同方法(ds方法和本发明方法)的两个对应重构结果。
为了突出建模结果的统计特征的对比分析，在此实验案例中每个方法都执行了20次重构，其20个输出结果的统计特征如图6所示。空间连通性是衡量三维地质结构模型重构质量的重要标准，图6(a)为两种方法的重构结果在x、y、z三个方向上的空间连通性曲线，从图中可以看住本发明方法的输出模型在各个方向上的连通性更加接近给定的三维参考模型。另一个地质结构各向异性的衡量标准是属性的变异性，图6(b)为两种方法的重构结果的整体变差函数曲线，其结果同样表明本发明方法在对地质结构的异质性重构能力上具有更强的优势。图6(c)为随着网格尺度的增加，两种方法计算效率的对比曲线，可以看出本发明方法随着模拟网格数目的增加，其计算效率的加速比也逐步增大，故该方法有能力实施大范围精细地质结构模型的自动重构。
本发明以低维数据(如二维剖面、一维的钻孔等取样数据)为基础，以横纵剖面代替完整的三维训练图像，实现了对复杂三维地质结构模型的自动重构，为从低维地质数据到三维地质模型的重构提供了新思路，同时也为基于多点统计的三维地质模拟和建模提供了新技术、新方法；本发明提供的局部搜索策略，将多点统计信息的获取固定到距离当前模拟位置最近的局部子截面，这将避免距离太远的多点模式被捕获，因此保证了其模拟结果的可靠性；而且局部子截面的搜索策略降低了对训练图像的搜索次数，从而提升了模拟效率；本发明中所采用多重网格以及搜索邻域的参数调整策略，允许搜索邻域的半径和模式距离的阈值随着网格级别的增加而递减；在最初的多重网格上，其邻域半径较大，模式距离的阈值也较大，这保证了大尺度结构模式的获取；而随着邻域半径和阈值的减小，模式匹配将趋向于完全匹配，这将保证细小的结构特征的重构；该策略在一定程度上减轻了传统基于多点统计的随机模拟方法中搜索模板固定且尺寸不宜太大的问题；本发明中的概率融合策略和权重选择方法，使得最终的联合概率分布包含了来自各个方向上子截面的概率分布；其中平行子截面的基于加法的概率融合拓展了概率分布的容量，使得融合后的联合概率更加稳定；正交方向上的基于乘法的概率融合保证了不同方向上多点模式的最终刻画，使得其模拟结果在不同方向上都表现出了很好的空间连续性；本发明中对最大匹配模式数目的限制，既保证了稳定的条件概率密度函数的获取，又避免了无效的多点统计搜索与模式匹配，这也在一定程度上提升了三维模型重构的效率；本发明可在各种三维地质信息系统、地理信息系统、地质建模与模拟系统等软件中推广应用。
在不冲突的情况下，本文中上述实施例及实施例中的特征可以相互结合。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。