专利名称：具有增强计算效率的多点统计(mps)仿真的制作方法
相关专利申请的交叉引用
本申请通过在此引用下列共同未决申请而全部并入本文
共同提交的“Method for Creating Facies Probability CubesBased Upon Geologic Interpretation(基于地质解释创建相概率立方的方法)”，代理人申请案编号T-6359。
共同提交的“Method of Making a Reservoir Facies ModelUtilizing a Training Image and a Geologically Interpreted FaciesProbability Cube(利用训练图像和地质解释的相概率立方建立储集相模型的方法)”，代理人申请案编号T-6404。
背景技术：
模拟油气层的常规地质统计工作流包括模拟相，然后使用基于方差图的算法为每个相填充岩石物理性质，典型地多孔性和渗透性。因为方差图仅是空间变化性的两点测量，基于方差图的地质统计不允许模拟逼真的曲线或其他几何上复杂的相型，例如蜿蜒的沙道，而这对于连通性和流动特性评价是关键性的。
称作多点统计仿真或MPS仿真的最近模拟方法已经由Guardiano和Srivastava，Multivariate GeostatisticsBeyond BivariateMoments(多变量地质统计超出二变量力矩)Geostatistics-Troia，in Soares，A.，ed.，Geostatistics-TroiaKluwer，Dordrecht，V.1，p.133-144，(1993)提出。MPS仿真是使用概念地质模型作为3D训练图像以产生地质逼真的储集层模型的储集相模拟技术。训练图像基于储集层体系结构模拟中的测井解释和一般经验提供地下地质体的概念描述。MPS仿真从训练图像中提取多点模式并且将模式锚定到井中数据。
已经发表了关于MPS及其应用的许多其他发表物。Caers，J.和Zhang，T.，2002，Multiple-point GeostatisticsA Quantitative Vehiclefor Integrating Geologic Analogs into Multiple Reservoir Models(多点地质统计一种将地质模拟结合到多储集层模型的定量载体)，在Grammer，G.M et al.，eds.，Integration of Outcrop and ModernAnalog Data in Reservoir Models(储集层模型中露头岩和现代模拟数据的结合)AAPG研究报告。Strebelle，S.，2000，SequentialSimulation Drawing Structures from Training Images(从训练图像中绘制结构的顺序仿真)博士论文，斯坦福大学。Strebelle，S.，2002，Conditional Simulation of Complex Geological Structures UsingMultiple-Point Statistics(使用多点统计的复杂地质结构的条件仿真)数学地质学，V.34，No.1。Strebelle，S.，Payrazyan，K.和J.Caers，J.，2002，Modeling of a Deepwater Trubidite ReservoirConditional to Seismic Data Using Multiple-Point Geostatistics(使用多点地质统计局限于地震数据的深海浊积岩储集层的模拟)，在2002SPE技术年会和展览会中提出的SPE 77425，San Antonio，Sept.29-Oct.2。Strebelle，S.和Journel，A，2001，Reservoir ModelingUsing Multiple-point Statistics(使用多点统计的储集层模拟)在2001SPE技术年会和展览会中提出的SPE 71324，New Orleans，Sept.30-Oct.3。
SNESIM(单个正则方程仿真)是特别对相和储集层模拟领域技术人员众所周知的一种MPS仿真程序。特别地，SNESIM仿真在Strebelle，S.，2000，Sequential Simulation of Complex GeologicalStructures Using Multiple-Point Statistics(使用多点统计的复杂地质结构的顺序仿真)博士论文，斯坦福大学和Strebelle，S.，2002，Conditional Simulation of Complex Geological Structures UsingMultiple-Point Statistics(使用多点统计的复杂地质结构的条件仿真)数学地质学，V.34，No.1中详细描述。基本SNESIM代码也可以在网站http://pangea.stanford.edu/～strebell/research.html获得。同样包括在网站中的是提供SNESIM背后的理论并且包括各种用例研究的PowerPoint演示文稿senesimtheory.ppt。PowerPoint演示文稿senesimprogram.ppt通过底层SNESIM代码提供指导。再次，这些发表物对于在创建相和储集层模型时使用多点统计的相模拟者是众所周知的。Strebelle的这些发表物在此引用其全部内容作为参考。
经验已经显示MPS仿真程序SNESIM相当好地再现训练模式。但是，SNESIM比同样在斯坦福大学研制的可比较的基于方差图的仿真程序SISIM需要更多的CPU。SNESIM需要非常大量的内存以从3D数百万节点训练立方中提取，然后存储模式。
在Strebelle学术论文(2000，pp.40-53)中描述的MPS仿真程序SNESIM与常规指示基于方差图的程序SISIM基于相同的顺序仿真范例。该描述SNESIM的浓缩包含在该说明书的附录A中。使用SISIM，仿真网格节点沿着随机路径访问仅一次。SISIM在Deutsch，C.和Journel，A.(1998)GSLIBGeostatistical Software Library andUser’s Guide(GSLIB地质统计软件库和用户手册)，第二版，牛津大学出版社，纽约中描述。一旦被仿真，节点值变成将限制按顺序随后访问的节点的仿真的硬数据。虽然在基于方差图的算法中，在每个未采样节点执行克里格法以估计仿真值将从其中绘制的局部条件分布，但是在MPS仿真中，局部条件分布通过使用调节数据的给定数据模板扫描训练图像而估计(Guardiano和Srivastava，1992)。Strebelle学术论文(2000)的主要贡献是通过在称作搜索树的动态数据结构中提前存储所有必需的条件概率分布函数(cpdf)来减少Srivastava原始代码的cpu时间。为了方便，该说明书的附录B描述如何生成Strebelle(2002)的搜索树。更确切地说，W(u)表示以位置u为中心的数据搜索窗口，并且τn表示由定义W(u)的n个位置u+hα的n个矢量{hα，α＝1...n}构成的数据模板。在仿真之前，训练图像使用数据模板τn扫描，然后与数据模板τn相关联的所有可能数据事件的出现次数存储在搜索树中。在MPS仿真期间，局部cpdf直接从该搜索树中检索。因此，训练图像不需要为每个节点仿真而重新扫描。
搜索树方法的一个主要限制是数据模板τn不能包括太多网格节点。对于这种限制有两个原因 1.用来构造搜索树的内存量随着数据模板的大小按指数规律增长对于取K个可能值的属性，例如K相值，与数据模板τn相关联的可能数据事件的最大数目是Kn。幸运地，该最大数目很少达到。
2.从搜索树中检索cpdf所需的cpu时间随着大数据模板τn而急剧增长。在任何未采样节点u，数据模板τn中仅n′(≤n)个数据位置由调节数据(原始采样数据或先前仿真的节点)填充有信息。以由那些n′个调节数据构成的数据事件dn′为条件的概率分布的推理需要计算训练图像中dn′的出现次数c(dn′)。该次数通过将与数据模板τn相关联并且包含dn′的所有数据事件dn的出现次数求和而获得
任何这种数据事件dn的出现次数c(dn)可以直接从搜索树中读出。调节数据的数目n′越小，包含dn′的可能数据事件dn的数目越大，从搜索树中检索所有相关次数c(dn′)所需的cpu时间越多。对于取K个可能值的属性，包含dn′的可能数据事件dn的数目可以多达Kn-n′。
数据模板因内存和cpu时间的考虑而不能包括太多网格节点。然而，数据搜索窗口应当足够大以捕获训练图像的大规模结构。缓和该冲突的一种解决方法是使用多网格方法，由此仿真G个嵌套且越来越细的网格，参看Tran，T.，Improving Variogram Reproduction onDense Simulation Grids(改进密集仿真网格上的方差图再现)，计算机与地质科学，20(7)1161-1168(1994)和Strebelle学术论文(2000，p.46)。第g(1≤g≤G)网格由最后仿真网格的每第2g-1个节点构成。应用于大小8×8＝64节点的2D网格的多网格仿真顺序的实例在图2中显示，其中G＝3嵌套网格。
数据模板与待仿真网格中节点的间隔成比例地重新调节，这允许在粗网格级使用大尺寸模板捕获大规模训练结构。第g网格由下一个第(g+1)网格的每2nd个节点构成。那么在2D中，第(g+1)网格中的节点数目是先前仿真的第g网格中节点数目的22＝4倍。这意味着，在第(g+1)网格仿真的开始，该网格中节点的1/4是属于第g网格的先前仿真(带有信息)的节点。因此，为了找到至少，比方说，20个数据以调节未采样(不带有信息的)节点的仿真，应当考虑由n＝4＊20＝80个节点位置构成的搜索数据模板τn。(在3D中，第(g+1)网格中带有信息的节点的数目是第g网格中不带有信息的节点数目的23＝8倍，并且将需要由n＝8＊20＝160个节点构成的数据模板。)对于取K＝2个可能值的属性，与这种模板τn相关联的可能数据事件的最大数目是Kn＝280＝1.3×1024，并且包含由n′＝20个数据构成的调节数据事件的可能完全带有信息的事件的最大数目可能多达幸运地，所有这些非常大的数目以正在扫描的训练图像的节点总数为上限。
本发明解决在SNESIM中使用的MPS方法的计算效率方面的前述缺点。
发明内容
在下文提出两种解决方法以减小数据搜索模板的大小(1)使用多网格仿真方法中的中间子网格，以及(2)选择由先前仿真的节点优先构成的数据模板。这两种解决方法的组合允许节省大量内存和cpu时间，仍然保证大规模的训练结构被捕获并输出到仿真练习。
本发明的目的在于在使用搜索树的多点统计仿真中使用减小大小的数据模板以增强计算效率。
另一个目的在于在MPS仿真中使用嵌套的网格，其中嵌套的网格通常同样由带有信息和不带有信息的节点组成。
再一个目的在于使用主要由带有信息的节点和很少数不带有信息的节点组成的数据模板以减小数据模板大小并增强计算效率。
本发明的这些和其他目的、特征和优点将关于下面的描述、未决的权利要求和附随附图而变得更好理解，其中 图1是在根据本发明的MPS仿真中执行的步骤的总体流程图，其中减小大小的数据模板用来增强仿真的计算效率； 图2A-C说明使用3层嵌套且越来越细的网格对于大小8×8＝64个节点的仿真网格，现有技术的多网格仿真顺序，其中先前仿真的节点是黑色，当前网格中待仿真的节点是灰色，以及该步骤中没有仿真的节点是白色； 图3A-E描绘使用中间子网格对于大小8×8＝64个节点的仿真网格，根据本发明的多网格仿真顺序(先前仿真的节点是黑色并且当前子网格中待仿真的节点是灰色)； 图4A-C显示使用中间子网格在3D中的多网格仿真顺序(先前仿真的节点是黑色并且当前子网格中待仿真的节点是灰色)； 图5A-C说明与细的8×8仿真网格相关联的第二子网格(先前仿真的节点是灰色，当前子网格中待仿真的节点是白色)，以及轮廓为黑色的搜索窗口；b)由包含在搜索窗口中的所有网格节点构成的现有技术数据模板；以及c)由先前仿真的节点以及四个与待仿真其属性的网格最接近的未仿真节点构成的新的数据模板； 图6A-D显示a)训练图像；b)用于原始仿真的80数据模板；c)使用原始多网格仿真方法产生的仿真实现；以及d)使用本发明的新的多网格仿真方法产生的仿真实现； 图7说明搜索树实例，包括a)训练图像；b)保留的数据模板；以及c)使用搜索数据模板从训练图像中获得的搜索树。每个节点下面的斜体数字是在图8中显示的表格中使用的节点参考数字； 图8显示给出扫描以构造图8C的搜索树的所有训练重复的中心节点u坐标的表格；以及 图9是图1的仿真步骤800的子步骤的流程图，其中网格值从搜索树中仿真。
具体实施例方式 图1是显示在本发明中执行的并且与常规SNESIM算法中使用的那些步骤相匹配的步骤的流程图。附录A是SNESIM算法如何工作的一般理论说明，并且是Strebelle学术论文的p.40-53的浓缩。本发明中在步骤500和600中进行显著改进，如下面将更详细描述的。比较和区别将在常规SNESIM算法和本发明的增强MPS仿真之间进行。
在常规SNESIM算法和本发明中，节点的地层网格在步骤100中创建，其用来模拟所关注的储集层。训练图像在步骤200中创建，其反映模拟者对可能在储集层中存在的地层模式和异质性的概念化。为了捕捉储集层的广阔概观以及限制计算中节点的数目，节点的初始粗网格在步骤300中选择，其对应于地层网格。然后在步骤400中使用井中数据和训练图像仿真属性，也就是相值。优选地，这些属性使用在下面更详细描述的MPS仿真步骤仿真。
在确定粗网格的这些初始属性之后，在步骤500中通过添加另外的节点到网格中来精制网格。该更细的网格，或者工作网格，包括先前仿真其属性的节点。这些节点称作带有信息的节点，因为已经指定属性。带有信息的节点可能包括已知的井中数据否则仿真值。添加到属性仍然未知的工作网格的另外节点称作不带有信息的节点。步骤500是在本发明中对基本SNESIM算法进行增强的步骤之一。
然后在步骤600中节点的数据模板从节点的该精制工作网格中选择。该步骤也已经在本发明中改进优于常规SNESIM实现。
然后在步骤700中使用步骤600的数据模板扫描训练图像以创建搜索树。然后不带有信息的节点的属性在步骤800中从搜索树中顺序地确定。关于搜索树的创建以及不带有信息的节点的属性如何可以确定的细节可以在附录B中找到。附录B来自Strebelle 2000学术论文。
然后检查工作网格以查看它的细度是否与地层网格的细度匹配。如果是，那么所有节点都已经指定属性。如果不是，那么增强工作网格的细度并且通过重复步骤500至800确定另外不带有信息的节点的属性，直到工作网格在数目上与地层网格匹配，且所有属性例如相类型已经确定。
常规SNESIM MPS仿真遭受计算密集。在下文提出两种解决方法以减小数据搜索模板的大小(1)在多网格仿真方法中使用中间子网格(步骤500)，以及(2)选择由先前仿真的节点优先构成的数据模板(步骤600)。这两种解决方法的组合允许节省大量内存和cpu时间，仍然保证大规模的训练结构被捕获并输出到仿真练习。
中间子网格的仿真 为了减小搜索数据模板的大小，每个网格中的中间子网格在先前描述的原始多网格方法中使用，除了第一(最粗)网格。参看图3A-E。在2D中，考虑两个子网格。与第(g+1)网格相关联的第一子网格由前一个第g网格的节点加上位于由这些节点构成的正方形中心的节点构成，如图3中对于大小8＊8＝64个节点的网格的仿真所示。该第一子网格中节点的一半是在第g网格中先前仿真的节点。因此，为了在每个未采样节点的附近找到至少20个调节数据，搜索模板应当包含2＊20＝40个节点位置。注意，这与像常规SNESIM算法一样直接仿真第(g+1)网格时80个必需节点相比较。第二子网格是第(g+1)网格自身，其节点数目是第一子网格的节点数目的两倍，参看图3。从而，再次，该子网格的节点的一半是先前仿真的节点，并且40个节点位置的数据搜索模板足够大以在每个未采样位置找到至少20个调节数据。
在本发明中，对于数据模板，具有相对高的带有信息的节点与总节点的比例是期望的。注意，在图2A中该比例是1/4。在图3中该比例更高，也就是1/2，其高于在常规SNESIM中发现的1/4比例。对于3D网格，常规SNESIM的比例是1/8。在本发明中，比例保持1/2。
在原始SNESIM多网格仿真方法中，原始采样数据，也就是井中数据，重新定位到当前仿真网格的最接近节点；当网格的仿真完成时，它们重新指定到它们的原始位置。在本发明的新的多网格仿真方法中，相同的数据重定位过程应用于每个子网格。
在3D中，提出以减小搜索模板大小的新的多网格方法需要每个嵌套网格3个子网格，以便像2D中一样，每个子网格中的节点数目是前一个子网格中的两倍，参看图4A-C。
由先前仿真的节点位置优先构成的数据模板 数据模板的大小可以通过优先保留与先前仿真的节点相对应的节点位置来进一步减小。考虑与图5A中显示的细的8＊8网格相关联的第二子网格的节点仿真。与图5A中所示的搜索窗口相对应的数据模板将由40个网格节点构成，如果保留该搜索窗口中的所有网格节点，参看图5B。代替地，提出数据模板由下列构成 ·属于前一个子网格的已经仿真的位置。在图5C中的实例中存在它们中20个；以及 ·不属于前一个子网格，但是接近待仿真中心位置的少量节点位置，比方说4个最接近的节点。在当前子网格的仿真开始， 这些位置不带有信息，除非它们对应于原始采样数据位置。
注意新提出的模板，如图5C中显示的，仍然提供原始搜索窗口限制内空间的充分均衡的覆盖。
对于取K＝2个可能值的属性，与图4C的24数据模板相关联的可能数据事件的数目是224＝1.7*107，对比图5B的原始40数据模板的240＝1.1*1012。
在每个未采样节点，可以找到至少20个调节数据。可以包含特定调节20数据事件的24数据事件的最大数目是224-20＝24＝16，与可以包含相同调节20数据事件的40数据事件的最大数目240-20＝210＝106相比较。
内存和CPU时间节省技巧 使用提出的中间仿真子网格和由先前仿真的节点位置优先构成的数据搜索模板允许减小数据模板的大小，从而节省内存需求和cpu时间。但是，精确地估计该节省量是困难的，因为下列原因 ·新的搜索树必须为每个子网格而构造。在2D中，由新的多网格仿真方法所需的搜索树的数目是由原始方法所需的搜索树数目的大约两倍，以及3D中该数目的大约三倍。但是，在内存和cpu时间方面那些另外搜索树的一次成本与从减小数据模板的大小中期望的节省相比较期望是较小的。回想搜索树的一次构造所需的cpu时间比仿真3D数百万节点仿真网格的多个实现所需的cpu时间少得多。
·在先前的分析中，提供可能的数据事件数目的数量级。因为训练图像显示重复模式并且具有有限大小N(典型地小于3D中几百万节点)，与数据模板τn相关联的训练数据事件的数目小于N，因此比与τn相关联的可能数据事件的总数Kn小得多。但是，对于取K＝2个可能值的属性，使用由24个节点位置(代替80)构成的模板仍然导致内存需求的显著降低，如由在下文提出的用例研究说明的。
·搜索树实际上包括n级。假设τn′＝{hα，α＝1...n′}是由最接近其中心的τn的n′个数据位置构成的模板。在搜索树的n′(≤n)级存储有与τn′相关联的大小n′的所有可能数据事件dn′的训练重复的数目。如果nmax是在以u为中心的数据搜索模板中找到的最远调节数据的索引，以数据事件dn′为条件的概率估计仅需要检索包含dn′的所有可能数据事件的出现次数dnmax。这种数据事件的次数dnmax可以多达Knmax-n′，其可能比包含dn′的数据事件dn的次数Kn-n′小得多。
·伪随机路径在原始SNESIM中使用以便首先仿真由最大数目调节数据填充信息的节点，参看Strebelle(2000，p.47)。但是，除了存在非常密集的原始采样数据集之外，对于这种最佳带有信息的节点找到的调节数据的数目在2D中保持接近包含在数据搜索窗口中的节点的1/4。提出的新的子网格方法将使该分数加倍。
从搜索椭圆体中的数据模板构造 在原始SNESIM算法中，数据模板由用户作为模板节点的相对坐标在其中指定的格式化文件而提供，参看GSLIB软件中的Geo-EAS格式(Deutsch和Journel，1998，p.21)；这种文件的实例在Strebelle论文(2000，p.169)中提供。使用新的多网格方法，每个子网格一个模板将需要由用户提供。另外，用户找出然后输入由先前仿真的节点位置优先构成的模板可能非常单调乏味。由于那两个原因，在由本发明提供的新的SNESIM版本中，模板优选地使用由用户提供的下列参数自动生成 ·由3个半径和3个角度参数定义的数据搜索椭圆体，如GSLIB用户手册的图II.3中描述为radius，radius1，radius2，sang1，sang2和sang3； ·与先前仿真的节点相对应的模板位置的数目nclose。在本发明的优选实施方案中，还没有仿真的接近节点的数目在2D中以及在3D中任意地设置为4。
使用这些参数，数据搜索模板如下构造 ·使用与该搜索椭圆体相关联的各向异性距离测量计算搜索椭圆体的每个网格节点到椭圆体中心的距离。
·以距离递增的顺序排列节点，并且仅保留属于当前仿真子网格的那些节点。
·保留属于前一个仿真子网格的nclose个最接近的节点，也就是nclose个最接近的先前仿真节点，以及四个还没有仿真的最接近节点。对于第一(最粗)仿真网格，因为不存在先前仿真的节点，简单地保留该网格的nclose个最接近的节点。
多网格仿真方法中嵌套网格数目的计算 在原始SNESIM算法中，用户必须指定在多网格方法中使用的嵌套网格的数目。根据本发明，基于由用户提供的数据搜索椭圆体和模板大小，使用下面的递归子例程计算该数目 ·设置嵌套网格的数目nmult为1。
·使用与搜索椭圆体相关联的各向异性距离测量，构造数据模板τn＝{hα，α＝1...nclose}使得在任何位置u，nclose个位置u+hα与最接近u的nclose个网格节点相对应。
·如果所有nclose个位置u+hα都在搜索椭圆体内，nmult递增1，通过将其所有矢量乘以2重新调节模板τn，并且检查该重新调节的模板的所有位置是否仍然位于搜索椭圆体内。如果是，nmult再次递增1，重新调节数据模板...直到至少一个(重新调节的)模板位置在搜索椭圆体外。
实例到河流储集层的2D水平截面的仿真的应用。
为了提供当减小数据模板的大小时节省的内存和cpu时间的估计，改进的SNESIM算法的性能与原始版本关于河流储集层的水平2D截面的仿真相比较。图6A中显示的训练图像描绘期望在地下存在的蜿蜒沙道的几何形状；该训练图像的大小是250＊250＝62,500像素，并且通道比例是27.7％。
在原始SNESIM的参数文件中，在多网格仿真方法中使用的嵌套且越来越细的网格的数目nmult设置为6。图6B中显示的数据模板τn用来构造6个搜索树；该模板包括80个网格节点。使用该模板获得的较大搜索树由2,800,000个节点组成(每个搜索树节点对应于与τn的子模板相关联并且其至少一次出现可以在训练图像中找到的一个数据事件)。
由原始SNESIM生成的一个无条件仿真实现在图5C中显示。使用400MHz的Silicon Graphic的Octane机，仿真所需的总时间是60.6cpu秒，其中构造然后删除6个搜索树花费28.8秒，并且从搜索树中检索局部cpdf花费27.8秒。
对于根据本发明的改进SNESIM版本，与先前仿真的节点相对应的模板数据位置的数目nclose设置为20。考虑各向同性搜索椭圆体radius＝radius1＝50，radius2＝1，sang1＝0，sang2＝0以及sang3＝0，距离单位与最后仿真网格的节点间隔相对应。使用那些参数值，在新的多网格仿真方法中使用的嵌套网格的数目计算为6。图6D显示使用本发明的改进SNESIM版本生成的一种无条件仿真实现。训练模式的再现类似于使用原始SNESIM生成的实现。一个实现花费总共16.0秒，其中 ·11.3秒用于构造然后删除11个搜索树。回想每个子网格一个搜索树，也就是除了(第一)最粗网格之外每个仿真网格2个搜索树被构造。因此，虽然改进的SNESIM要求更多搜索树，但是cpu时间除以2.5。另外，较大的搜索树仅使用100,000个节点，对比原始SNESIM的2,800,000个节点。
·3.6秒用于从搜索树中检索所有必需的cpdf，这比原始SNESIM少八倍。
虽然在前述说明书中本发明已经关于其某些优选实施方案而描述，并且许多细节已经为了说明而陈述，但是对于本领域技术人员将显然，本发明容易更改并且在这里描述的某些其他细节可以相当多地改变，而不背离本发明的基本原理。
附录A 使用搜索树的多点统计仿真-一般理论 术语 考虑取K个可能状态{sk，k＝1，...，K}的属性S。S可以是分类变量，或者其可变性间隔由(K-1)个阈值离散化成K类的连续变量。使用下列术语 以待仿真的位置u为中心的大小n的数据事件dn由下列构成 -由n个矢量{hα，α＝1，...，n}定义的数据几何形状 -n个数据值{s(u+hα)＝α＝1，...，n 该数据事件的中心值是待估计的未知值，表示为s(u)。
数据模板τn仅包括前一个数据几何形状。τn的子模板是由τn的矢量的任何子集n′构成的模板，其中n′≤n。
数据事件dn所述与几何模板τn“相关联”。与τn相关联的cpdf(条件概率分布函数)是以与τn相关联的特定数据事件dn为条件的中心值s(u)的概率分布。
扩展正则方程 考虑取K个可能状态{sk，k＝1，...，K}的属性S。我们希望估计以最接近的n个硬数据α＝1，...，n为条件的变量S(u)的局部概率分布。
A0表示与状态sk在位置u出现相关联的二进制(指示)随机变量
类似地，假设D是与由联合考虑的n个调节数据α＝1，...，n构成的数据事件dn的出现相关联的二进制随机变量
D可以分解成与每个调节数据相关联的二进制随机变量的乘积 其中
先前的分解允许(广义)指示克里格体系或者扩展正则方程的应用。
Prob{A0＝1|D＝1}+E(A0)+λ[1-E{D}](1) 因此，由指示克里格法提供的精确解标识条件概率的定义，如由贝叶斯关系给出的。
扫描训练图像 精确解(2)要求超出传统两点方差图或协方差模型很多的(n+1)点统计。从实际(地下)数据中推理这种多点统计通常是没有希望的，因此通过在平稳性的事前决策下扫描一个或多个训练图像来借用它们的想法(输出许可) ·表达式(2)的分母α＝1，...，n可以通过计数在训练图像中找到的调节数据事件α＝1，...，n的重复次数c(dn)来推理。重复应当具有相同的几何构造和相同的数据值。
·分子Prob{S(u)＝Sk且α＝1，...，n通过计数c个先前重复中与中心值S(u)等于Sk相关联的重复次数ck(dn)而获得。然后所需的条件概率识别为训练比例ck(dn)/c(dn) SNESIM算法 为提出的仿真算法而选择的snesim名字回想它仅涉及一个单个正则方程(sne)，这导致条件概率的真正定义(2)。SNESIM已经研制以仿真分类属性，例如地质相，但是可以扩展以使用一系列嵌套仿真来仿真连续属性例如渗透性场 ·首先将连续属性值离散化成有限数目K个类。在大部分储集层应用中，因为流动仿真由与特定相(例如沙和页岩)相关联的极端渗透性的连通性控制，不需要细的离散化。四类渗透性值应当足够，例如，两个极端类对应于最低和最高十分位数，并且两个中间类每个具有边缘概率40％。
·使用snesim仿真作为结果的分类变量，然后在每个分类(值的类)中使用传统两点算法例如GSLIB顺序高斯仿真程序sgsim(Deutsch和Journel，1998，p.170)仿真原始连续变量。SNESIM算法基于顺序仿真范例，由此每个仿真值变成调节按顺序随后访问的节点值的仿真的硬数据值(Goovaerts，1997，p.376)。Guardiano和Srivastava(1993)提出在每个未采样节点重新扫描全部训练图像以推理类型(3)的局部条件概率分布。这种重复扫描可能是非常需要cpu的，特别是当考虑大的训练图像时或者当产生每个具有许多节点的大量实现时。
一种备选实现将包括提前将所有必需的cpdf(条件概率分布函数)制成表格。但是，大量可能调节数据事件排除任何蛮力制表，实际上 ·每个取K个可能值或K类值的n个数据变量的单个模板产生K数据事件；例如K＝10类和n＝10导致Kn＝1010种可能的调节数据事件！ ·在顺序仿真模式中，因为调节数据包括先前仿真的节点，并且网格节点沿着随机路径访问，调节数据事件的几何形状从一个节点到另一个节点而有所不同。因此应当保留非常大量的数据模板(几何构造)。
在该论文中提出的算法实现比Guardiano和Srivastava的实现cpu需求少得多，且没有太多内存(RAM)需求。该新的实现基于下面两个性质 ·性质1给定n个数据变量的模板τn，可以实际地从训练图像中推理的与τn相关联的cpdf数目与训练图像尺寸有关，因此通常比与τn相关联的cpdf的总数Kn小得多。
考虑n个数据变量的模板τn。与τn相关联的cpdf可以从训练图像中推理，只要相应调节数据事件的至少一次出现可以在训练图像中找到。Nn表示可以由τn扫描的侵蚀训练图像τn的大小。因为仅Nn个数据事件在Tn中扫描，可以实际地从训练图像推理的与τn相关联的cpdf的最大数目必然小于Nn，因此通常是和与τn相关联的全部cpdf的巨大数目Kn相对比的适当数目。
·性质2以与τn的子模板τn′(n′≤n)相关联的数据事件dn′为条件的概率分布可以从以与τn相关联且dn′是其子集的数据事件dn为条件的概率分布中检索。
假设dn′是与τn的子模板τn′(n′≤n)相关联的数据事件。在附录A中，我们显示dn′的重复次数c(dn′)等于与τn相关联且dn′是其子集的所有数据事件dn的重复的总和
类似地对于中心值S(u)等于Sk的dn重复次数Ck(dn′)
然后ck(dn′)和c(dn′)的知识允许使用关系(3.6)估计以dn′为条件的概率分布。
W(u)表示以位置u为中心的数据搜索邻域。考虑由n个矢量{hα，α＝1，...，n}构成的数据模板τn，其中定义n个矢量使得n个位置u+hα，α＝1，...，n对应于在W(u)中存在的全部n个网格节点。snesim算法在两个步骤中继续 1.首先，在动态数据结构(搜索树)中仅存储可以实际地从训练图像推理的与τn相关联的那些cpdf。更确切地说，在搜索树中存储实际地在训练图像上找到的数据事件和中心值的出现次数(ck(dn))，并且训练比例(3.6)可以从其中计算。附录A的部分1提供该搜索树的详细描述。因为性质1，由搜索树所需的RAM量不会太大，如果保留具有适当数目n个节点，比方说，小于100个节点的数据模板τn。该搜索树的构造需要在图像仿真之前扫描训练图像仅一次，因此它非常快，参看附录A的部分2。
2.接下来，通过沿着随机路径访问每个网格节点一次执行仿真。在待仿真的每个节点u，调节数据在W(u)中搜索，因此局部调节数据事件与τn的子模板相关联。根据性质2，局部cpdf然后可以从以与τn相关联且dn′是其子集的数据事件dn为条件的概率分布中检索，这些cpdf直接从搜索树中读取。任何局部cpdf的快速检索在附录A的部分3中详细描述。训练图像不需要在每个未采样节点重新扫描，这使得该snesim算法比Guardiano和Srivastava的实现快得多。
流程图显示snesim仿真算法的主要步骤。
1.扫描训练图像以为对应于保留的数据搜索邻域的数据模板τn＝{hα，α＝1，...，n}构造搜索树。可以定义τn使得n个位置u+hα，α＝1，...，n对应于最接近u的n个网格节点位置，但不是必要地。
n个矢量hα根据递增模数排序 |h1|≤|h2|≤...≤|hn| 注意，各向异性方差图距离可以用于该排序。
2.将原始采样数据指定到最接近的网格节点。定义访问一次且仅一次所有未采样节点的随机路径。
3.在每个未采样位置u，保留在由以u为中心的模板τn定义的搜索邻域中实际存在的调节数据。假设n′是那些调节数据的数目(n′≤n)，并且dn′是相应调节数据事件。从搜索树中检索u的中心值等于sk，k＝1，...，K的训练dn′重复的次数ck(dn′)，如附录A的部分3中所示。识别局部cpdf为对应于数据事件dn′的类型(3.6)的比例。
为了保证这些比例有意义，从而避免局部cpdf的较差推理，如果训练dn′重复的总数小于输入的最小值Cmin，丢弃最远的调节数据，将调节数据的数目减小到(n′-1)；以该较少的数据事件dn′-1为条件的概率分布再次从搜索树检索，等等...。如果数据的数目降至n′＝1，并且c(dn′)仍然低于cmin，条件概率p(u；Sk|(n))；由边缘概率pk代替。
4.从搜索树检索的cpdf中绘制节点u的仿真s值。然后该仿真值增加到s数据以用于调节所有随后节点的仿真。
5.沿着随机路径移动到下一个节点，并且重复步骤3和4。
6.循环直到所有网格节点被仿真。一副随机图像已经产生。使用不同的随机路径从步骤2开始重复整个过程以产生另外的实现。
硬数据的调节 必须满足两个条件以保证硬数据的适当调节 ·硬数据必须在它们的位置准确再现。
·随着我们接近任何硬数据位置，条件方差应当变小，缩减到块金方差。更确切地说，节点u的L个仿真值{S(l)，l＝1，...，L}的方差应当随着节点接近硬数据位置uα而减小。
将采样数据重新定位到最接近的仿真网格节点并且冻结它们的值保证第一条件。
随着节点u接近硬数据位置u1，训练概率分布{p(u；Sk|(n))，k＝1，...，K}变得接近S(u)＝sk的单原子分布，如果S(u1)＝Sk。实际上，根据训练图像的空间连续性，随着u1接近u任何数据模板的中心值s(u)将越来越经常地处于与调节数据值s(u1)相同的状态中。训练图像的小规模空间连续性通过训练比例传递到仿真实现，因此当它确实在训练图像上时，随着仿真节点接近硬数据，它的条件方差减小。
多网格仿真方法 由数据模板τn定义的数据搜索邻域不应当取得过小，否则将不能再现训练图像的大规模结构。另一方面，包括太多网格节点的大搜索模板将导致在搜索树中存储大量cpdf，增加cpu成本和内存需求。
捕获大规模结构同时考虑包括适当少量网格节点的数据模板τn的一种解决方法是多网格方法(
1991；Tran，1994)。在snesim中实现的多网格方法包括仿真G个嵌套且越来越细的网格。第g(1＜g＜G)网格由最后仿真网格(g＝1)的每第2g-1个节点构成。数据模板τn＝{hα，α＝1，...，n}与待仿真网格内节点的间隔成比例地重新调节。假设是第g网格的作为结果的数据模板α＝1，...，n。较粗仿真网格的较大搜索邻域(由τn9定义)允许捕获训练图像的大规模结构。
一个搜索树需要根据嵌套的仿真网格，可能使用反映该规模特有的异质性的不同训练图像构造。当第g网格的仿真完成时，它的仿真值冻结为数据值以用于在下一个更细的仿真网格上的调节。
注意，原始采样数据不需要位于当前仿真的网格的节点。适当调节需要将原始采样数据指定到当前仿真网格的最接近节点。当网格的仿真完成时，这些采样数据重新指定到它们的原始位置。采样数据从其中去除的网格节点作为下一个更细网格的节点随后仿真。
但是，将原始采样数据重新定位到当前仿真网格的最接近节点可能影响仿真实现的局部准确度。因此仅更细的网格应当使用搜索树仿真。对于较粗的网格，每个未采样节点的局部cpdf可以使用Guardiano和Srivastava的原始实现，也就是通过在每个节点重新扫描全部训练图像来推理，在该情况下不需要数据重定位。因为全部训练图像的这种重新扫描是需要cpu的，仅非常粗的网格，比方说由最后仿真网格的每第8个(或更多)节点(其代表网格节点总数的仅1.5％)应当不使用搜索树而仿真。
在多网格方法中，可以对于每个仿真网格使用不同的训练图像。训练图像可以具有不同于初始仿真网格大小的大小(像素数目)，但是它们的分辨率(像素大小)必须相同。搜索树的构造非常快，但是当考虑非常大的训练图像或者具有大量节点的模板时可能非常需要内存。因此我们在相应仿真网格的仿真之前为任何单个搜索树分配内存，然后一旦网格的仿真完成则解除分配该内存。
附录B 在搜索树中存储CPDF 该部分表现动态数据结构(搜索树)，从训练图像推理的条件概率分布存储在其下，然后在使用snesim执行的仿真期间检索。
B.1使用的数据结构 考虑取K个可能状态{Sk，k＝1，...，K}的属性S。W(u)表示以任何未采样位置u为中心的数据搜索邻域，并且考虑由n个矢量{hα，α＝1，...，n}构成的数据模板τn，其中定义n个矢量使得n个位置u+{hα，α＝1，...，n}对应于在W(u)中存在的全部n个网格节点。n个矢量hα根据递增的模数排序 |h1|≤|h2|≤...≤|hn| 各向异性方差图距离可以用于该排序。
在下文呈现的动态数据结构允许检索在训练图像中存在的所有条件概率分布(cpdf)，假设调节数据配置包括在τn中。
该动态数据结构(搜索树)包括连接的一组节点，其中每个节点指向下一级的多达K个节点。作为实例，图A.1c显示使用图A.1b的数据模板τ4，从图A.1a中显示的大小5＊5＝25个像素的二进制训练图像中获得的搜索树。状态值0对应于白色像素，状态值1对应于黑色像素。
每个搜索树节点对应于一个数据事件。假设τn′＝{hα，α＝1，...，n′}是由τn的前n′个矢量构成的子模板，定义子模板τn′的n′个位置u+{hα，α＝1，...，n}对应于数据搜索邻域W(u)中最接近u的n′个位置。位于搜索树的n′(∈
)级的节点对应于与τn′相关联的数据事件。特别地，树从那里开始生长的0级的单个节点称作‘根’，并且对应于模板中不存在调节数据的数据事件d0。
仅对应于至少一个重复在训练图像中找到的数据事件的节点在搜索树中存在。假设考虑关联到数据事件α＝1，...，n′的节点 ·该节点包含K个整数的阵列{ck(dn′)，k＝1，...，K}，其中ck(dn′)是中心值S(u)等于sk的dn′的训练重复次数。那么dn′的总数是 ·一组K个指针{Pk(dn′)，k＝1，...，K}连接到dn′节点。Pk(dn′)指向与数据事件dn′+1＝dn′，且s(u+hn′+1)＝sk相对应的n+1级的节点，如果dn′+1节点在搜索树中存在，也就是，如果至少一个dn′+1重复在训练图像中找到。如果该dn′+1节点缺少，Pk(dn′)是‘空’指针，这意味着它不指向任何另外的节点。
用图解法，在每个节点，搜索树分成多达K个分支。
在图A.1c的二进制搜索树中，关联到任何数据事件dn′的节点包含中心节点是白色像素(c0(dn′))或黑色像素(c1(dn′))的训练dn′的次数。两个指针连接到dn′+1节点左指针P0(dn′)指向与数据事件{dn′，且s(un′+1)＝0}相关联的n′+1级的节点(如果该节点不缺少)，而右指针P1(dn′)指向与数据事件{dn′，且s(un′+1)＝1}相关联的n′+1级的节点(如果不缺少)。
例如，在根(关联到没有调节数据的数据事件d0的节点1)，c0(d0)＝14且c1(d0)＝11，对应于在图A.1a中显示的训练图像的14个白色和11个黑色像素。P0(d0)指向关联到一点数据事件{s(u+h1)＝0}的1级的节点(节点2)，而P1(d0)指向关联到{s(u+h1)＝1}的节点(节点3)。
T表示图A.1a的训练图像，并且T表示由τn的前n′个矢量构成的子模板τn′扫描的侵蚀训练图像Tn′＝{u∈Tst.u+hα∈T，α＝1，...，n′。由子模板τ1＝{h1}扫描的侵蚀训练图像T1由T的前四行(j＝1，...，4)构成。d1＝{s(u+h1)＝0}的12次重复可以在T1中找到。中心值s(u)对于那些d1重复中5个是白色，对于7个其他d1重复是黑色，因此在节点2中c0({s(u+h1)＝0})＝5且c1({s(u+h1)＝0})＝7。所有d1重复的中心节点u的坐标在图A.2中显示的表格中给出。
类似地，{s(u+h1)＝1}的8次重复可以在侵蚀训练图像T1中找到5个具有白色中心值，3个具有黑色中心值，因此在相应搜索树节点3中c0({s(u+h1)＝1})＝5且c1({s(u+h1)＝1})＝3。
B.2预扫描以构造搜索树 为数据模板τn＝{hα，α＝1，...，n}构造搜索树需要扫描由T表示的训练图像一次；它如下继续 1.为根(对应于n′＝0的数据事件d0)分配内存，并且对于k＝1，...，K初始化次数ck(d0)＝0。
2.在(全部)训练图像T的每个网格节点u，递增与在u处观察的状态Sk相对应的次数ck(d0)。
Tn′＝{u∈Tst.u+hα∈T，α＝1，...，n′表示由τn的前n′个矢量构成的子模板τn′扫描的侵蚀训练图像 TnTn-1…T1T0＝T 假设nmax是该侵蚀训练图像的下标，使得 保留子模板τnmax的nmax个位置uαuα＝u+hα，α＝1。因为nmax个矢量hα根据递增的模数排序，nmax个位置uα根据到u的递增距离而排序 表示以u为中心的nmax数据事件。
3.从根开始，考虑与由dnmax的第一数据构成的数据事件相对应的1级的节点。如果d1节点缺少，在搜索树中分配内存以创建它。如果u的中心值是Sk，则移动到该节点并且递增次数ck(d1)。然后，考虑与由dnmax的前两个数据构成的数据事件相对应的2级的节点，依此类推。循环通过数据s(uα)，α＝1，...，nmax的序列，a＝1，…，直到到达与dnmax相对应的nmax级的节点并且递增次数ck(dnmax)。
4.移动到训练图像T的下一个网格节点u并且重复步骤2和3。
5.循环直到训练图像的所有网格节点u都已经访问。
考虑与由τn的前n′个矢量n′∈
构成的子模板τn′相关联的任何数据事件dn′。假如在扫描训练图像之后，与dn′相对应的节点在搜索树中缺少。这意味着没有dn′的重复在训练图像中找到。在这种情况下，ck(dn′)＝0，k＝1，...，K。例如，在图A.1c的搜索树中，与由4个白色像素(d4＝{s(uα)＝0，α＝1，0，...，4})构成的数据事件相对应的节点缺少，因为没有该数据事件的重复在训练图像中找到c0(d4)＝c1(d4)＝0。
不存储与缺少的节点相对应的空次数ck(dn′)对于最小化内存需求是关键的。实际上，Nn1表示侵蚀训练图像T1的大小(像素数目)。因为仅Nn1数据事件在Tn1中扫描，搜索树的n′级的节点数目必然小于Nn1例如，由图A.1b的数据模板τ4扫描的T4的侵蚀训练图像由图A.1a中显示的原始训练图像的9个中心像素构成。因此，作为结果的搜索树的最后一级的节点数目必然小于N4＝9，虽然，使用二进制变量，四数据模板可能产生总共24＝16种可能数据事件。实际上，仅6个节点存在于搜索树的4级。实践中，训练图像构造得越多，搜索树中节点数目越少。
B.3从搜索树中检索CPDF 考虑以与τn的子模板相关联的任何n′数据事件dn′为条件的概率分布的推理。数据事件dn′可以与由τn的前n′个矢量定义的子模板τn′的一个相关联例如，与τ1＝{h1}相关联的但是这不需要如此例如 ·如果dn′与由τn的前n′个矢量定义的特定子模板τn′相关联，相应dn′节点可以在搜索树中如下搜索 从根(0级)移动到与由dn′的第一数据值构成的数据事件相对应的1级的节点开始，如果d1节点在搜索树中存在。如果d1节点缺少，那么dn′节点也缺少，这意味着不存在dn′的训练重复，因此以dn′为条件的概率分布不能从训练图像推理。如果d1节点不缺少，移动到与由dn′的前2个数据值构成的数据事件相对应的2级的节点，如果d2节点存在，依此类推，直到dn′节点到达搜索树的n′级。然后以dn′为条件的概率分布可以使用次数ck(dn′)估计 ·如果数据事件dn′不与这种子模板τn′相关联，那么在搜索树中不存在相应dn′节点，dn′重复在训练图像中找到或没有找到。例如，关联到数据事件d1＝{s(u2)＝0}的节点在图A.1c的搜索树中不存在，虽然d1重复可以在图A.1a的训练图像中找到。对于这种数据事件考虑包含dn′的最小特定子模板τnmin(由τn的前nmin个矢量构成) 可以应用一般等式Prob{A}＝∑bProb{AandB＝b}，其中求和在由事件B可能采取的所有状态b上执行，以写出下面的关系
其中求和在与子模板τnmin相关联且dn′是其子集的所有数据事件dnmin上执行。根据训练数据事件
因为在求和中涉及的数据事件dnmin关联到特定子模板τnmin，次数ck(dnmin)可以直接从搜索树的nmin级的相应dnmin节点检索。
注意这种计算允许仅检索在侵蚀训练图像τnmin中找到，而在由dn′数据配置扫描的侵蚀训练图像中没有找到的那些dn′重复。但是，当考虑大训练图像时，这是微小的限制。
作为实例，假设使用图A.1c的搜索树估计以数据事件d1＝{s(u2)＝0}为条件的概率分布
次数ck(d2)可以直接从搜索树检索
因此c0(d2)＝6且c1(d2)＝1，并且以{s(u2)＝0}为条件的概率分布可以由下面估计
权利要求
1.一种仿真地层网格中属性的方法，该方法包括
(a)创建节点的地层网格；
(b)创建代表地下地质异质性的训练图像；
(c)创建与地层网格的节点相对应的节点的粗网格；
(d)利用井中数据仿真粗网格的节点的属性以获得带有信息的节点；
(e)通过添加不带有信息的节点到带有信息的节点来精制节点的网格，从而创建节点的工作网格；
(f)从工作网格中选择节点的数据模板并且使用数据模板和训练图像构建搜索树；
(g)使用搜索树仿真工作网格的不带有信息的节点的属性；以及
(h)重复步骤(e)-(g)直到地层网格中节点的属性都已被仿真；
其中在精制步骤(e)的至少一个中，带有信息的节点与节点总数的比例对于2D网格大于1/4，并且对于3D网格大于1/8。
2.根据权利要求1的方法，其中
带有信息的节点与节点总数的比例在全部精制步骤(e)中对于2D网格大于1/4，并且对于3D网格大于1/8。
3.根据权利要求1的方法，其中
在精制步骤(e)的至少一个中，带有信息的节点与节点总数的比例是1/2。
4.根据权利要求1的方法，其中
在多数精制步骤(e)中带有信息的节点与节点总数的比例是1/2。
5.根据权利要求1的方法，其中
对于精制步骤(e)的每个，带有信息的节点与节点总数的比例是1/2。
6.根据权利要求1的方法，其中
选择节点的数据模板的步骤包括创建搜索窗口并且选择搜索窗口中带有信息的节点和仅一部分不带有信息的节点以在步骤(f)的至少一个中创建节点的数据模板。
7.根据权利要求6的方法，其中
搜索窗口中不带有信息的节点的少于1/2在步骤(f)的至少一个中被选择为数据模板的一部分。
8.根据权利要求1的方法，其中
选择节点的数据模板的步骤包括创建搜索窗口并且选择搜索窗口中带有信息的节点和不带有信息的节点，其中多数所选节点是带有信息的节点。
9.根据权利要求1的方法，其中
选择节点的数据模板的步骤包括创建搜索窗口并且选择搜索窗口中全部带有信息的节点和仅一部分不带有信息的节点。
10.根据权利要求1的方法，其中
搜索窗口是椭圆体。
11.根据权利要求1的方法，其中
属性是相。
12.一种仿真地层网格中属性的方法，该方法包括
(a)创建节点的地层网格；
(b)创建代表地下地质异质性的训练图像；
(c)创建与地层网格的节点相对应的节点的粗网格；
(d)利用井中数据仿真粗网格的节点的属性以获得带有信息的节点；
(e)通过添加不带有信息的节点到带有信息的节点来精制节点的网格；
(f)从工作网格中选择节点的数据模板并且使用数据模板构建搜索树，该数据模板通过创建搜索窗口并且选择搜索窗口中带有信息的节点和不带有信息的节点来选择，其中多数所选节点是带有信息的节点；
(g)使用搜索树仿真网格的不带有信息的节点的属性；以及
(h)重复步骤(e)-(g)直到地层网格中所有节点的属性都已被仿真。
13.根据权利要求12的方法，其中
在精制步骤(e)的至少一个中，带有信息的节点与节点总数的比例对于2D网格大于1/4，并且对于3D网格大于1/8。
14.根据权利要求12的方法，其中
在搜索窗口中找到的不带有信息的节点的少于1/2被选择为数据模板的一部分。
15.根据权利要求12的方法，其中
选择节点的数据模板的步骤，多数所选节点是带有信息的节点。
16.根据权利要求12的方法，其中
选择节点的数据模板的步骤，选择搜索窗口中全部带有信息的节点和仅一部分不带有信息的节点。
17.根据权利要求12的方法，其中
搜索窗口是椭圆体。
全文摘要
一种增强的多点统计(MPS)仿真。使用一种多网格仿真方法，其从5常规MPS方法中改进以减小数据搜索模板的大小，6在仿真期间节省大量内存和cpu时间。用来减小数据搜索模板的大小的特征包括(1)在多网格仿真方法中使用中间子网格，以及(2)选择由先前仿真的节点优先构成的数据模板。这些特征的组合允许优于先前的MPS算法而节省大量内存和cpu时间，仍然保证大规模训练结构被捕获并输出到仿真练习。
文档编号G06G7/48GK101133422SQ200580028073
公开日2008年2月27日 申请日期2005年8月16日 优先权日2004年8月20日
发明者瑟巴斯汀·B·斯卓伯乐 申请人:切夫里昂美国公司