基于二维图像和多点统计方法的多孔介质三维重构方法

专利名称：基于二维图像和多点统计方法的多孔介质三维重构方法
技术领域：
本发明涉及一种构图方法，具体涉及一种多孔介质三维图像的重构方法。
背景技术：
以达西定律为基础的渗流力学取得了很大的进展，在众多工程领域有广泛的应 用。但是由于达西定律根本不涉及孔隙介质本身的细观空间结构，因此对细观层次上的渗 流形貌研究无能为力。为克服达西定律的局限性，研究者们提出了孔隙网络模型、毛管模型 和球形颗粒堆积模型等多孔介质孔隙模型。这些模型中的孔隙和喉道常常被简化成一些简 单的几何体。如孔隙体会被简化成球形、立方体或棱柱体；喉道会被简化成圆管、方形管、三 角形管道或星形的管道。但是由于多孔介质内部结构十分复杂，而上述模型的建立并不以 真实的多孔介质三维结构数据为基础，因此并不能很好地定量描述多孔介质内部不规则的 拓扑结构。可见获得真实的多孔介质三维结构已成为制约渗流机理研究的重要瓶颈。
随着实验技术的进步，诸如聚焦离子光束法，激光扫描显微法和X射线CT扫描等 方法被应用于获取三维孔隙图像。虽然一些实验设备的精度可以达到微米甚至是亚微米量 级，然而，一些多孔介质的纳米量级结构(例如石灰岩和硅藻土)却无法直接三维成像。相 比三维孔隙图像，多孔介质的二维薄片图像较易获取并可得到较高的分辨率。例如美国加 州大学在2004年利用实验设备获得精度达到10纳米的石灰石二维薄片图像。
利用二维薄片图像来重构三维孔隙图像的方法目前主要有三类。第一类是过程 法，它通过模拟地质成岩的过程来重构三维图像。过程法虽然可以再现多孔介质内部的长 连通性，但是该方法过于复杂。例如在模拟石灰岩中不规则形状的沉积、压实和成岩过程 时，该方法就显得过于复杂而难以实现。 第二类方法利用从二维图像中获取的孔隙度和两点自相关函数等统计信息来重 构三维孔隙图像，但是由于该类方法只是使用了二维图像中的低阶统计信息，因此难以再 现孔隙空间的长连通性，尤其是难以再现低孔隙度或具有特定孔隙几何形状(例如颗粒和 球形)的多孔介质拓扑结构。 第三类方法由0kabe等人提出，其利用多点统计方法来重构多孔介质的孔隙结 构，但是该方法在构建训练图像时只是通过旋转二维平面获得三维情况下的条件概率分布 函数(cpdf， conditional probability distribution function)，并没有获得真正的三维 训练图像。 利用多点统计方法重构图像的过程需要反映真实几何结构及其分布模式的训练 图像，它强调使用训练图像把先验模型明确而定量地引入到建模当中。先验模型包含了被 研究的真实物质中确信存在的样式，而训练图像则是该模型的定量化表达。通过再现高阶 统计量，多点统计方法能够从训练图像中捕捉复杂的(非线性)特征样式并把它们复制到 重构图像中。可以说训练图像中的概率信息决定了最终的模拟结果。而0kabe等人假定多 孔介质为各向同性，认为铅直方向的cpdf完全等同于水平方向的cpdf，因此该方法重构的 多孔介质三维结构无法反映真实情况下铅直方向的孔隙和骨架结构信息。

发明内容
本发明针对上述现有方法在重构多孔介质三维结构时所存在的缺陷，而提供一种 基于二维图像和多点统计信息的方法对多孔介质进行三维重构。利用同步辐射光技术获取 多孔介质的体数据，将该体数据的一个二维切面图像作为训练图像。提取训练图像的孔隙 和骨架点作为硬数据(条件数据)，利用多点统计信息重构该训练图像的下一层二维图像。 再将重构得到的新的二维图像作为训练图像模拟其下层图像。如此重复多次，即利用每次 重构得到的新的二维图像作为训练图像来模拟其下层图像。将每层重构图像依次叠加，最 终得到多孔介质的三维结构。 为了达到上述目的，本发明采用如下的技术方案 基于二维图像和多点统计方法的多孔介质三维重构方法，该方法包括如下步骤
(1)利用二维的多重数据模板扫描二维训练图像，建立搜索树；
(2)利用采样点提取模板遍历多孔介质的二维训练图像提取采样点，并将提取得 到的采样点作为初始条件数据分配到最近的网格节点上；(3)定义一条随机路径访问所有 待模拟节点；对随机路径上的每一个节点，利用与步骤(1)中相同的多重模板提取其条件 数据事件，然后从搜索树上获取该点的条件概率分布函；利用计算机随机模拟方法提取该 点的随机模拟值，并将该模拟值作为后续模拟新增的条件数据；继续模拟随机路径上的其 他节点，直到模拟生成一幅新的二维图像； (4)将由步骤(3)模拟得到的新的二维图像作为步骤(1)中的训练图像，重复步骤 (1)到步骤(3)生成其下一层图像； (5)循环步骤(1)到步骤(4)直至得到N幅新的二维模拟图像，并将每幅二维图像
的每个像素对应成三维空间的一个体素，而体素状态值即其对应像素点的状态值；将步骤
(1)中最初的训练图像和这N幅图像依次叠加，形成N+l层的多孔介质三维结构。 所述步骤(1)中最初的训练图像为微米精度的多孔介质二维薄片图像。 所述步骤(1)中进行扫描二维训练图像时，可采用逐步减少数据模板大小的方法
以使数据事件在训练图像中获得足够多的重复。 所述步骤(3)定义随机路径时，可以先利用数据模板遍历整个待模拟区域，统计
每个待模拟点周围条件数据的数目，优先模拟那些条件数据较多的节点。 根据上述技术方案得到的本发明在进行重构多孔介质三维结构时能够反映真实
情况下铅直方向的孔隙和骨架结构信息，并能重现真实情况下空隙复杂的长连通性特征，
够获得真实的多孔介质三维结构，突破制约渗流机理研究的瓶颈。 本发明使用具有更高精度的二维薄片图像作为训练图像，从而可以克服三维直接 成像分辨率较低的不足。 本发明提供的方法进行重构时具有较快的模拟速度，大大的提高效率。同时本发 明还可以广泛应用于地球科学、生物学、医学等相关领域中。


以下结合附图和具体实施方式
来进一步说明本发明。
图1为本发明的流程图。
2为数据模板扫描训练图像的过程示意图。 3为二维3重网格结构示意图。 4为与图3相对应的三重数据模板示意图。 5a为砂岩体数据的外表面图。 5b为砂岩体数据在三个方向的正交剖面图。 5c为砂岩体数据隐藏骨架后所得到的孔隙结构图。 6为图5所示图像的孔隙度曲线图。 7a为利用多点统计方法生成二维重构图像
7b为利用多点统计方法生成二维重构图像 7c为利用多点统计方法生成二维重构图像 8a为利用多点统计方法重构三维多孔介质 8b为利用多点统计方法重构三维多孔介质 8c为利用多点统计方法重构三维多孔介质 9a为多孔介质三维重构图像和体数据在X 9b为多孔介质三维重构图像和体数据在Y 9c为多孔介质三维重构图像和体数据在Z
时的初始训练图像。 时的训练图像的采样点图. 时的重构图像。 图像的外表面图。 图像的正交剖面图。 图像的孔隙结构图。 方向的变差函数曲线图。 方向的变差函数曲线图。 方向的变差函数曲线图。
具体实施例方式
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结 合具体图示，进一步阐述本发明。 本发明提出一种新的基于多点统计信息的方法对多孔介质进行三维重构的方法。 该方法先利用同步辐射光技术获取多孔介质的体数据，将该体数据的一个二维切面图像作 为训练图像；再提取训练图像的孔隙和骨架点作为硬数据(条件数据)，利用多点统计信息 重构该训练图像的下一层二维图像；接着再将重构得到的新的二维图像作为训练图像模拟 其下层图像；如此重复多次，即利用每次重构得到的新的二维图像作为训练图像来模拟其 下层图像；将每层重构图像依次叠加，最终得到多孔介质的三维结构。真实情况下的体数据 提供了与重构结果进行比较的参照数据。将该体数据与重构的多孔介质三维结构进行比较 的结果发现，重构的多孔介质的孔隙结构特征与真实情况吻合。同时该方法可以使用具有 更高精度的二维薄片图像作为训练图像，从而可以克服三维直接成像分辨率较低的不足。
多点统计方法是地质统计学的一个新的分支。地质统计学最初主要用于解决矿 床普查勘探、矿山设计到矿山开采整个过程中的各种储量计算和误差估计问题，后来它在 石油工业中得到迅速的发展，主要应用于储层表征与建模。多点统计方法是相对于传统的 两点地质统计方法而言的，它们的主要区别在于条件概率确定方法的不同。在多点统计方 法中，使用训练图像代替变差函数体现地质体的空间结构。多点统计方法通过在训练图像 中寻找与待模拟点周围条件数据分布完全相同的事件的个数来确定概率分布，因而可以反 映多点之间的相关性和变异性，克服基于变差函数(Variogram)的两点地质统计方法的不 足。 基于上述原理，本发明的具体实施如过程如下(如图1):
(1)利用二维的多重数据模板扫描二维训练图像，建立搜索树。
该步骤中利用数据模板获得多点统计信息，Strebelle等对数据模板和数据事件 进行了定义。数据模板是由n个向量仏。；a = 1，2， ...， n}构成的几何形态。设模板中 心位置是u，模板其他位置u。 =u+ha(a =l，2，...，n)。例如图2 (a)就是一个5 X 5像素 的二维模板，它由中心点u和其他24个向量所确定。假定一种属性S可取m个状态值{sk; k= l，2，...，m}。由数据模板中n个向量u。位置的n个属性值所组成的"数据事件"d(u) 可以定义为 =1， 2，…，w} (1) 其中i(ua)表示在u。位置的状态值。图2(b)是一幅15X15像素的训练图像。
图2(c)表示利用图2(a)所示的数据模板扫描训练图像时获得的一个数据事件。 利用数据模板扫描训练图像是为了统计一个数据事件d(u)出现的概率，即数据
事件中的n个数据点i (Ul) ， i (u2) ， . . . ， i (un)分别处于某个状态值、时该数据事件出现的
概率Prob{<iO)}=Prob{=、 ； a=l， 2， ...， "} (2) 在应用任一给定的数据模板对训练图像扫描的过程中，当训练图像中的一个数据 事件与数据模板的数据事件d(u)相同时，称为一个重复。在平稳假设的前提下，数据事件 d(u)在有效的训练图像中的重复数c(d(u))与有效的训练图像的大小Nn的比值，相当于该 数据事件d(u)出现的概率Pr0b{/(Wa)=\; "=1， 2，…，(3) 对于任一待模拟点u，需要确定在给定n个条件数据值i (ua)的情况下，属性i (u) 取m个状态值中任一个状态值的条件概率。根据贝叶斯条件概率公式，该概率可表达为
Prob(&|减")}= -d 二 、——"，:"、-1 (4) 上式中，分母为条件数据事件出现的概率；分子为条件数据事件和待模拟点U取 某个状态值的情况同时出现的概率，相当于在已有的C(d(U))个重复中i(U) = Sk的重复
的个数Ck(d(U))与有效的训练图像的大小Nn的比值，记为Ck(d(U))风。因此i(U)取某个 状态值的条件概率分布函数可表示成 Pr=& I /(O = ^ ;" = l,L ,"} (5) 基于上述条件概率，利用Monte Carlo方法可以提取该点的状态值。对于多孔介 质而言，待模拟点i(u)只可能取两种状态值孔隙或者骨架。利用模拟点周围最近的n个 孔隙与骨架值和扫描训练图像获得的cpdf，可以随机推断该点是孔隙还是骨架。由于采用 了概率估计的方法，故模拟结果具有随机性。这些结果是对训练图像先验模型的一种反映， 可以揭示属性空间中状态值的各种可能分布。 最初的多点统计方法在每次模拟一个节点时都要重新扫描一遍训练图像，以获得 对应节点的cpdf，严重影响了重构图像的速度。本发明利用一种叫"搜索树"的数据结构来 加速重构过程，只要对训练图像进行一次遍历搜索就可以生成该搜索树。所有可能的概率 模式被存贮在搜索树中，生成模拟图像时可以直接从搜索树上获得该模式的条件概率，从而加快了模拟速度。 本步骤在利用数据模板获得多点统计信息时，由于在训练图像中想要包含所有可 能的数据事件是不可能的，而且由于受到计算机内存和cpu速度的限制，只能选择合适的 数据模板大小。由于数据模板的大小有限，在重构的图像上其实只反映了训练图像上有限 尺寸下的结构特征。本发明采用了Tom Tran提出的多重模板的思想，利用网格逐渐密集化 的多个数据模板来替代一个大而密集的模板对训练图像进行扫描。具体方法是先使用稀 疏的粗网格数据模板扫描训练图像，得到粗网格下的多点统计信息，然后可以模拟得到粗 网格下的结果图像；将粗网格下的内容作为条件数据复制到细网格上，然后使用细网格模 板扫描训练图像，得到细网格下的多点统计信息，最后模拟得到细网格下的结果图像。具体 过程如图3至4所示 图3表示一个包含13X13个节点的二维3重网格结构，其中已模拟点用黑色表 示，待模拟点用灰色表示，白色点暂时忽略不作考虑。图4(a)、图4(b)、图4(c)分别是与图 3(a)、图3(b)、图3(c)各重网格相对应的三重数据模板。图4中的灰色点表示扫描图像时 待模拟点在数据模板中的位置，白色点不作考虑。可以看出，图4(a)与图3(a)中的水平和 垂直方向的相邻待模拟点之间的距离均为3个像素，故可以利用图4(a)所示的第一重数据 模板扫描图3(a)中的待模拟点，生成的模拟结果作为条件数据复制到第二重网格中(见图 3(b)中的黑色节点)。然后利用图4(b)所示的第二重数据模板扫描图3(b)，生成的模拟结 果和第一次模拟的结果均作为条件数据复制到第三重网格中(见图3(c)中的黑色节点)。 最后利用图4(c)所示的最细网格扫描图3(c)，可以模拟生成最终的重构图像。
在完成第l步后，进入到第2步利用边长为奇数的正方形的采样点提取模板提取 采样点，即，将采样点提取模板在多孔介质二维训练图像上遍历，如果在提取模板范围内全 部为孔隙或者骨架，则规定该模板中心为孔隙或骨架点。然后将提取得到的采样点作为初 始条件数据分配到最近的网格节点上，即，比较各采样点与其周围的网格点之间的位置，采 样点与哪个网格点最近，就分配到哪个网格点上。该步骤的具体过程如下
由于多孔介质图像由多个孔隙和骨架的连续区域所组成。可以近似认为，如果在 一层图像中有某个足够大的孔隙或者骨架区域(设该区域中心点为Center)，那么其下层 图像与Center相同位置处的状态值也应该为孔隙或者骨架。定义一个边长为奇数的正方 形的采样点提取模板T，利用T在二维训练图像中依次遍历，当在模板T中的节点状态值全 为孔隙或骨架时，就提取该模板中心节点为采样点，同时保持该采样点的状态值不变。不过 在提取采样点时应注意确保采样点中孔隙点的比例接近该训练图像的孔隙度，因为当采样 点中各状态值的比例接近其在训练图像中相同状态值的比例时，模拟效果最好。
步骤(3):定义一条随机路径访问所有待模拟节点；对随机路径上的每一个节点， 利用与步骤(1)中相同的多重模板提取其条件数据事件，然后从搜索树上获取该点的条件 概率分布函数；利用计算机随机模拟方法提取该点的随机模拟值，并将该模拟值作为后续 模拟新增的条件数据；继续模拟随机路径上的其他节点，直到模拟生成一幅新的二维图像。 该步骤中在定义随机路径时，可以先利用模板遍历整个待模拟区域，统计每个待模拟点周 围条件数据的数目，优先模拟那些条件数据较多的节点。 步骤(4):将由步骤(3)模拟得到的新的二维图像作为步骤(1)中的训练图像，重 复步骤(1)到步骤(3)生成其下一层图像。
7
步骤(5):循环步骤(1)到步骤(4)直至得到N幅新的二维模拟图像，并将每幅二 维图像的每个像素对应成三维空间的一个体素，而体素状态值即其对应像素点的状态值； 将步骤(1)中最初的训练图像和这N幅图像依次叠加，形成N+1层的多孔介质三维结构。
本发明实施时如果在数据模板中心点u周围的条件数据越多，那么数据事件包含
的有效节点就会越多，从而就难以在训练图像中找到较多的与该数据事件相同的重复。数 据事件重复的数目较少意味着该数据事件可能具有特殊性，这会导致模拟结果引入训练图
像的特殊模式，而不是其包含的一般模式。因此可以采用逐步减小数据模板大小的方法使 得数据事件可以在训练图像中找到足够多的重复。具体方法是如果重复数c(d(u))小于 某个设定的下限rmin，那么就先去除离u最远的节点，此时的条件数据数目就变为n-l。在 搜索树中寻找对应于这n-l个节点情况下的条件概率，如果此时的c (d(u))仍然小于rmin， 就继续去除现在数据模板中距u最远的节点，然后在搜索树中寻找对应这n-2个条件数 据的条件概率。如此重复下去，如果数据模板中的条件数据数目下降到n二 l，并且此时 c(d(u))仍然小于rmin，那么就用各状态值sk的边缘概率来作为u的条件概率。
根据上述技术方案得到的本发明具体实验如下(如图5-8所示)
在本实施例中，实验所用多孔介质体数据是在北京同步辐射装置4W1A束线形貌 站采集的。实验使用的X射线能量是24KeV，波长为0.lA。探测器使用的是X射线CCD成像 系统，其像素分辨率为10. 9 ii m，灰度分辨率为8bits，曝光时间100ms。实验样品为直径3mm 的圆柱形砂岩，平均孔隙度在17%左右。 本实施例取80X80X80体素的砂岩体数据进行实验，图5 (a)为砂岩体数据的外 表面图。图5(b)为该体数据在三个方向的正交剖面图(X二40，Y二40，Z二40)，图5(c) 为该体数据隐藏骨架后所得到的孔隙结构图。将体数据在水平方向的各层二维切面图像从 上向下编号(1到80)，附图6为这些图像的孔隙度曲线图，平均孔隙度为0. 1778。
在本实施例中，对体数据进行水平方向的截取，可获得其二维的水平切面图像。从 体数据中选取一幅接近其平均孔隙度的水平切面图像(80X80像素)作为初始训练图像 (该图像在体数据中的层号为35，孔隙度为0. 1787)，如图7(a)所示。采样点提取模板T设 置为3X3，提取图7(a)中0. 5%的像素点作为采样点，并且确保其中孔隙点的比例接近训 练图像的孔隙度。实施例中孔隙点占采样点的比例是O. 1563。图7(b)是训练图像的采样 点图，图7(c)是利用多点统计方法重构的结果图像。可以看出图7(c)具有与图7(a)相似 的孔隙结构，再现了图7(a)中孔隙不规则的特征模式。 在本实施例中，不断地将新产生的二维重构图像作为训练图像，并提取每幅训练 图像中的采样点作为条件数据，再利用多点统计方法模拟其下层图像。根据上述方法，共得 到79幅二维重构图像。将初始训练图像与这些重构图像依次叠加，得到80X80X80体素 的多孔介质三维重构图像。该重构图像的平均孔隙度为0.1789。图8(a)为该三维重构图 像的外表面图，图8(b)为其在三个方向的正交剖面图(X = 40， Y = 40， Z = 40)，图8(c) 为附图8(a)隐藏骨架后所得到的孔隙结构图。 通过比较重构的三维多孔介质图像(见图8a_图8c)与真实情况下的多孔介质体 数据(见图5a-图5c)，可以看出重构图像在铅直和水平方向具有与真实情况相似的不规则 的孔隙和骨架结构，并且重现了真实情况下孔隙复杂的长连通性特征。
变差函数能够反映地质变量在某个方向上空间结构变化的相关性。如果两幅图像
8中的某个状态值在同一个方向上具有相似的变差函数曲线，那么可以说明这两幅图像中的该状态值在此方向上具有相似的结构特征。可以获得多点统计方法重构的三维多孔介质图像和体数据在X、Y、Z方向的孔隙变差函数曲线，如图9(a)、图9(b)、图9(c)所示，多点统计方法重构图像和体数据在X、 Y和Z方向的孔隙变差函数呈现出相似的变化趋势，说明该重构图像与真实情况下砂岩的孔隙结构非常接近。 上述利用本发明提供的方法重构砂岩样品的三维拓扑结构，得到的结果图像再现了与真实情况相似的多孔介质的结构特征。而最初的训练图像是二维平面图像，与直接获取三维图像的方法相比，二维图像较易获取并拥有更高的分辨率，因此对于一些只能获取二维平面图像的地质结构体而言，利用本发明提供的方法都可以对其进行拓扑结构重构。
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
9
权利要求
基于二维图像和多点统计方法的多孔介质三维重构方法，其特征在于，所述方法包括如下步骤(1)利用二维的多重数据模板扫描二维训练图像，建立搜索树；(2)利用采样点提取模板遍历多孔介质的二维训练图像提取采样点，并将提取得到的采样点作为初始条件数据分配到最近的网格节点上；(3)定义一条随机路径访问所有待模拟节点；对随机路径上的每一个节点，利用与步骤(1)中相同的多重模板提取其条件数据事件，然后从搜索树上获取该点的条件概率分布函；利用计算机随机模拟方法提取该点的随机模拟值，并将该模拟值作为后续模拟新增的条件数据；继续模拟随机路径上的其他节点，直到模拟生成一幅新的二维图像；(4)将由步骤(3)模拟得到的新的二维图像作为步骤(1)中的训练图像，重复步骤(1)到步骤(3)生成其下一层图像；(5)循环步骤(1)到步骤(4)直至得到N幅新的二维模拟图像，并将每幅二维图像的每个像素对应成三维空间的一个体素，而体素状态值即其对应像素点的状态值；将步骤(1)中最初的训练图像和这N幅图像依次叠加，形成N+1层的多孔介质三维结构。
2. 根据权利要求1所述的基于二维图像和多点统计方法的多孔介质三维重构方法，其特征在于，所述步骤(1)中最初的训练图像为微米精度的多孔介质二维薄片图像。
3. 根据权利要求1所述的基于二维图像和多点统计方法的多孔介质三维重构方法，其特征在于，所述步骤(1)中进行扫描二维训练图像时，可采用逐步减少数据模板大小的方法以使数据事件在训练图像中获得足够多的重复。
4. 根据权利要求1所述的基于二维图像和多点统计方法的多孔介质三维重构方法，其特征在于，所述步骤(3)定义随机路径时，可以先利用数据模板遍历整个待模拟区域，统计每个待模拟点周围条件数据的数目，优先模拟那些条件数据较多的节点。
全文摘要
本发明公开了基于二维图像和多点统计方法的多孔介质三维重构方法，该方法将微米精度的多孔介质二维薄片图像作为初始的训练图像，利用多点统计方法重构其下层图像。将每次获得的下层图像作为新的训练图像，利用多点统计方法重构该训练图像的下层图像；二维重构图像的每个像素对应成三维空间的一个体素，依次叠加每层的二维重构图像得到最终的三维多孔介质图像。本发明能够获得真实的多孔介质三维结构，突破制约渗流机理研究的瓶颈，可广泛应用与能源和工程领域。
文档编号G06T11/00GK101706966SQ20091019843
公开日2010年5月12日 申请日期2009年11月6日 优先权日2009年11月6日
发明者张挺, 杜奕 申请人:上海第二工业大学