CN104335042B - 从多孔介质的数字表征中选择代表性单元体积的高效方法 - Google Patents

从多孔介质的数字表征中选择代表性单元体积的高效方法 Download PDF

Info

Publication number
CN104335042B
CN104335042B CN201380029060.XA CN201380029060A CN104335042B CN 104335042 B CN104335042 B CN 104335042B CN 201380029060 A CN201380029060 A CN 201380029060A CN 104335042 B CN104335042 B CN 104335042B
Authority
CN
China
Prior art keywords
volume
sub
volumes
sample
objective function
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201380029060.XA
Other languages
English (en)
Other versions
CN104335042A (zh
Inventor
朱塞佩·德·普里斯科
乔纳斯·特尔克
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ingrain Inc
Original Assignee
Ingrain Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ingrain Inc filed Critical Ingrain Inc
Publication of CN104335042A publication Critical patent/CN104335042A/zh
Application granted granted Critical
Publication of CN104335042B publication Critical patent/CN104335042B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N33/00—Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
    • G01N33/24—Earth materials
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N33/00—Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
    • G01N33/24—Earth materials
    • G01N33/241—Earth materials for hydrocarbon content
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/40—Imaging
    • G01N2223/419—Imaging computed tomograph
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/60—Specific applications or type of materials
    • G01N2223/649—Specific applications or type of materials porosity

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Physics & Mathematics (AREA)
  • Pathology (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Medicinal Chemistry (AREA)
  • Environmental & Geological Engineering (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • Food Science & Technology (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Theoretical Computer Science (AREA)
  • Radiology & Medical Imaging (AREA)
  • Pulmonology (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)
  • Image Analysis (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明涉及用以估计多孔介质的样本中的代表性单元体积即REV的方法，其中与现有方法相比，所选择的子体积是单元体积的更好近似。可以定义诸如岩石等的多孔介质的样本中的REV，其中相对于经由该多孔介质的流体流动的预期方向来选择该REV。该方法可以鉴定岩石的数字表示的良好程度以及通过达西定律对流体流动进行描述的准确程度，并且使得能够在沿不同方向应用该方法的情况下，对REV的沿不同方向的不同长度尺度进行评价并且对孔隙结构的各向异性进行评估。该方法还可以确定用以推定由于子样本的大小不足而导致孔隙率‑渗透率的趋势何时崩溃的鲁棒标准。

Description

从多孔介质的数字表征中选择代表性单元体积的高效方法
本申请根据35U.S.C.§119(e)而要求2012年7月11日提交的在先美国专利申请13/546,053和2012年3月30日提交的美国临时专利申请61/618,265的优先权，在此通过引用包含其全部内容。
背景技术
本发明涉及用于预测流体经由诸如多孔岩石等的多孔介质的流动的属性的方法和系统，尤其涉及用于从异质多孔介质的数字表示中选择最具代表性的子样本以用于预测诸如孔隙率、渗透率和/或相关特性等的属性的这些方法和系统。
可以经由x射线计算机断层图像扫描、扫描电子显微镜、共焦显微镜和其它技术来产生诸如岩石、骨、土壤和其它材料等的多孔介质的数字表示。这些数字表示用于使用计算机模拟来表征多孔介质(Knackstedt,M.A.等人，“Digital Core Laboratory:Properties of Reservoir Core Derived from 3DImages”,Society of Petroleum Engineers,2004；和Vermeulen,J.P.,“NewDevelopments in FESEM Technology”,Carl Zeiss nano-technology SystemsDivision,http://www.zeiss.com/C1256E4600307C70/EmbedTitelIntern/NewDevelopmentinFESEMTechnology/$File/New_Development_FESEM_Technology.pdf.)。
多孔介质特性的数字模拟中的重要问题是样本大小。诸如多孔岩石等的许多实际关注的样本是异质的，并且大体积的多孔介质的一般属性将需要对非常大量的样本进行数字化。诸如绝对渗透率等的许多岩石特性需要大量的计算资源来进行模拟，结果与该关注体积的代表性特性相比，样本大小通常小得多。可以由训练有素的地质学家从视觉上选择子样本，但该方法主观且充满变数。此外，基于岩石特性的数字模拟所作出的诸如投资井、井射孔方案、可回采储量估计和其它这种决策等的商业和技术决策通常涉及巨额费用。结果，需要去除表现这种多孔介质的特性时的主观性、错误和变化。
用以识别适合的子样本的一个方法是识别代表性单元体积(REV)。REV是可以进行将得到代表整体的值的特定测量的最小体积。REV以下的体积在该特定测量中呈现出变化，这使得小于REV的样本不适合进行模拟。在Bear发表的文献(Bear,J.,Dynamics of Fluids in Porous Media；General PublishingCompany Ltd.,Canada,1972,pp.19-21)中描述了利用体积孔隙率作为测量值来计算REV的方法。标记为REV的许多方法结果并不真正是“单元”。也就是说，使用中的许多方法可以得到代表较大体积的较大体积的子体积，但这些方法可能无法产生最小可能体积或最小单元体积。
Razavi等人描述了与REV共通的方法(Razavi等人，“RepresentativeElementary Volume Analysis of Sands Using X-Ray Computed Tomography,”Geotechnical Testing Journal,Vol.30,No.3,Paper ID GTJ100164,2006)。在本发明的图1中示出Razavi等人所描述的方法所用的流程图。在Razavi等人所示的方法中，选择样本的近似中心的点，然后在该中心点周围检查球形子样本体积。针对球形子样本计算样本属性。增大子样本的半径并重新计算这些属性。逐级地增大子样本体积，直到达到REV为止。该方法存在多个缺点。该方法在某些异质样本中可能无法得到合适的结果。尽管这可以得到可接受的RV，但可能无法得到REV。如上所述，对岩石样本的数字表示的计算会需要大量的计算机时间来完成，因而确定样本内的最小的REV很有价值。
美国专利6,516,080(Nur)公开了用于从代表性区域中选择REV的方法。本发明的图2示出位于样本的一个面的中心的方形区域的大小被增大，直到得到代表性区域为止。然后，选择方形代表区域的边长作为三维样本中居中的立方体的边的长度。该方法依赖于均质且各向同性的样本。这种样本无法代表诸如井岩心等的许多现实世界样本。
美国专利申请公开2011/0004447(Hurley等人)涉及如下方法，其中该方法用于使用在样本的两个或更多个深度处取回两组或更多组所发送的测量数据的至少一个测量工具来表现多孔介质的三维样本的特性。在该方法中，通过以下操作来估计孔隙率代表性单元体积(REV)：(1)从所测量或所建模的样本中随机选择大小均匀的多个非重叠块；(2)标绘各个块孔隙率vs相应的块体积；以及(3)针对给定块体积确定针对各样本所测量到的孔隙率之间的方差。孔隙率是所选择的样本内的孔隙率的平均值。在所测量到的孔隙率的方差降至所选择阈值以下的情况下，相应体积是研究中的岩石的孔隙率REV。Hurley等人的该方法并不从点开始增大体积，并且如此将覆盖将有效地减小样本大小的更多可能子体积。该方法的缺点在于，该方法被设计成使用多个子样本以使得可以获得统计相关方差，并且需要使用大的样本以实现期望收敛，其中这两个需求并非始终可能，并且可能给出原始样本整体作为RV。本研究者已认识到这是对于经过激光扫描共焦显微镜(LSCM)的样本的情况。Hurley等人的方法也可能无法识别样本内的最小可能REV。
用以表现多孔介质的微观结构的有趣且强大的方法是Hilfer所提出的局部孔隙率理论的随机分析(1992)。该方法以尺度依赖方式被公式化，并且该方法给出REV的积分长度尺度的良好估计。然而，局部孔隙率方法没有给出与多孔空间的各向异性有关的结果。Liu等人对该方法作了如下改进(2009和2010)：对以尺度依赖方式进行评价的局部各向异性分布进行评价。该改进需要应用Ketcham所提出的方法(2005)，其中在该方法中，各向异性是孔隙结构的特性的定向变化的函数。
使用达西定律(Darcy’s Law)来进行对诸如岩石等的多孔介质的属性的许多估计。达西定律是以现象学推导的、对流体经由多孔介质的流动进行近似的方程式。Henry Darcy基于他对水经由沙床的流动所进行的实验的结果将该定律公式化。达西定律本质上是动量守恒的表达。由于达西定律经常应用于经由诸如岩石样本等的多孔介质的流动，因此达西定律可用于利用具有诸如图24所示等的达西流动参数的以下方程式1来估计体积流量。
Q = - k A μ P b - P a L
其中，
Q＝试样的一个相内的每单位时间的体积流速
k＝多孔介质的绝对渗透率
A＝流动的截面面积
μ＝动态粘性
Pb,Pa＝体积的入口和出口处的压力
L＝样本的长度
形式上，为了推导达西定律以定义渗透率，例如，必须根据第一原理验证一些假设。特别地，如Whitaker,S.,Transport in Porous Media 1,1986,pp.3-25所示，根据Navier Stokes方程式(即，动量方程式)推导达西定律的方式是应用Gary的分解：
P = P ‾ + p ~
其中该分解从根本上是尺度的分解：是假定相对于(如横向或纵向维度那样、可以是样本的长度尺度的)平均积分尺度“良好表现”的平均量(在这种情况下为压力)。换句话说，假定这些平均函数充分描述它们所表示的量。例如，在与平均长度尺度相当的长度尺度内快速改变的压力信号无法表示该长度内的压力。量是压力的波动部分，其中该量表示函数的变化。假设是平均量在允许波动部分具有小的变化的小尺度上不改变。为了推导达西定律，必须向Navier-Stokes方程式应用平均运算(例如，体积平均法)与Gary的分解。然而，在这种情况下，获得了场的梯度的平均值，而该平均值是所期望的平均量的梯度(如达西定律所示)。可以很容易地证明，两个算子(梯度和平均值)在应用于在平均长度尺度内不快速改变的函数的情况下互换。特别地，如果孔隙率均匀，则可以将达西定律写成如以下方程式2那样的更一般表示法。
q ( x ) ‾ = - k ( x ) μ ▿ P ( x ) ‾
其中：
k＝位置x处的多孔介质的绝对渗透率
μ＝动态粘性
使用该方程式，可以看到平均量随位置而改变的尺度，并且不再看到波动部分。该方程式可用于模拟储层内的流动。
在利用上述任意方式选择REV的情况下，存在如下可能性：存在REV内的孔隙率的变化，从而使与达西流动有关的假定无效或容易出错。此外，压力梯度可以沿着流动方向快速地改变，这样使得无法定义与特定子样本相关联的渗透率。这特别适用于诸如在真实世界岩石形成中可以发现的样本等的高度异质样本。
本研究者认识到需要用以估计包括异质样本的多孔介质的样本中的代表性单元体积(REV)的更加高效的方法。此外，该分析必须考虑孔隙结构的定向变化，从而考虑多孔介质的各向异性、以及在考虑到流动那样的定向属性的情况下的流动的方向。
发明内容
本发明的特征是提供用以估计诸如岩石等的多孔介质的样本中的代表性单元体积(REV)的高效方法，其中与先前方法相比，所选择的子体积是单元体积的更好近似。
本发明的另一特征是提供用以定义诸如岩石的多孔介质的样本中的REV的高效方法，其中相对于经由多孔介质的流体流动的预期方向来选择REV。
本发明的另一特征是提供用以对岩石的数字表示的好(或坏)程度、以及通过达西定律对流体流动进行描述的准确程度进行定量的高效方法。
本发明的另一特征是提供一种确定用于推定由于子样本的大小不足而导致孔隙率-渗透率的趋势何时崩溃的鲁棒标准的方法。
本发明的另一特征是提供用于以尺度依赖方式来分析(包括异质的变化的定向信息的)孔隙结构的方法。
为了实现这些和其它优点，并且根据本发明的目的，如这里所体现的以及从广义上所述，本发明在一部分中涉及一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：a)获得表现孔隙空间和至少一个固相的特性的分割体积；b)针对所述分割体积整体，推导第一目标函数P1的平均属性值<P1>；c)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差σvol；d)在所述体积内定义多个子体积；e)针对各所述子体积，计算所述第一目标函数P1的属性值P相对于平均属性值<P1>的标准偏差σi；f)求出标准偏差σi与σvol良好地匹配的所有候选代表性子体积；g)从所述候选中选择并存储代表性子体积；以及h)使用所述代表性子体积来推导至少一个关注属性值。
本发明还涉及一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：a)获得表现孔隙空间和至少一个固相的特性的分割体积；b)使所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向；c)通过对与所述所定义的流动方向垂直的数字切片进行分析，来针对所述分割体积整体推导作为至少第一目标函数P1的一个或多个函数的值；d)在所述体积内定义多个子体积；e)关于所述所定义的流动方向，针对各所述子体积计算至少第一目标函数P1的一个或多个函数的值；f)求出函数识别体积值和子体积值之间的匹配的所有代表性子体积候选；g)从所述候选中选择代表性体积；h)存储所述代表性子体积；以及i)使用所述代表性子体积来进行模拟或推导至少一个关注属性值。
本发明还涉及一种用于从多孔样本的较大三维数字图像中获得代表性单元体积的高效估计的方法，包括以下步骤：a)获得表现孔隙空间和至少一个固相的特性的分割体积；b)针对所述分割体积整体，推导作为至少第一目标函数P1的至少一个函数的值；c)在所述体积内定义多个子体积，其包括：定义子体积的初始大小，利用具有所定义的初始大小的子体积填充所述体积整体，以及对其它子体积的大小进行迭代，以及利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止；d)针对各所述子体积，计算作为至少第一目标函数的至少一个函数的值；e)针对良好地匹配的所述体积和所述子体积的值，求出所有代表性子体积候选；f)从所述候选中选择并存储代表性子体积；以及g)使用所述代表性子体积来进行模拟或推导至少一个关注属性值。
本发明还涉及一种用于从多孔样本的较大三维数字图像中获得代表性单元体积的高效估计的方法，包括以下步骤：a)获得表现孔隙空间和至少一个固相的特性的分割体积；b)将所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向；c)使用针对以与所述所定义的流动方向垂直的方式所截取的样本体积的多个数字切片的分析，来针对所述分割体积整体推导第一目标函数P1的平均属性值<P1>；d)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差；e)在所述体积内定义多个子体积，其包括：定义子体积的初始大小，利用具有所定义的初始大小的子体积填充所述体积整体，以及对其它子体积的大小从大到小地进行迭代，利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止；f)关于所述所定义的流动方向，针对各所述子体积计算属性P相对于平均属性值<P1>的标准偏差σi；g)求出σi与σvol良好地匹配的所有候选代表性子体积；h)选择最小的候选并且存储所述最小的候选作为代表性单元体积；以及i)使用所述代表性单元体积来推导至少一个关注属性值。
本发明还涉及一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：1)将多孔介质的分割三维图像加载至计算机系统中，其中，所述分割三维图像包括各自分配有灰度值的体素；2)选择被定义为Z方向的流动方向；3)定义探询体积的大小，其中：i.所述探询体积是维度为Xi、Yi和Zi的原始分割三维图像的子样本，其中整个样本的维度为Xs,Ys,Zs，ii.定义探询体积的最大值imax，iii.设置各所述探询体积的以体素为单位的维度(Xi,Yi,Zi)，其中Xi、Yi和Zi是针对值为1～imax的i所定义的，以及iv.将i的初始值设置为1；4)针对所述探询体积的各切片来计算所选择属性Ps(0,0,0)～Ps(0,0,Zs)；5)计算σs(0,0,0)；6)设置大小为Xi,Yi,Zi的所述探询体积在大小为Xs,Ys,Zs的整个样本内所占据的最大坐标，其中：amax＝Xs–Xi+1，bmax＝Ys–Yi+1，以及cmax＝Zs–Zi+1；7)将当前探询体积的位置坐标设置为a＝b＝c＝0；8)针对所述当前探询体积的切片计算所选择属性Pi(a,b,c)～Pi(a,b,c+Zi)，其中，所述所选择属性包括孔隙率、表面积与体积之比、相似属性或它们的任意组合；9)计算σi(a,b,c)，其中，可选地，对计算σi的值所使用的Pi的值进行滤波，以及其中，可选地，设置Pi的平均值；10)使所述探询体积的位置在X方向上移动1个体素、即a＝a+1；11)重复步骤8)～10)并且存储Pi和σi的所有值，直到所述当前探询体积的X坐标的值a等于所述当前探询体积能够占据的最大值amax为止；12)将所述当前探询体积的X坐标设置为零即a＝0，并且使当前位置体积的Y坐标增加1个体素、即b＝b+1；13)重复步骤8)～12)并且存储Pi和σi的所有值，直到所述当前探询体积的Y坐标的值b等于所述当前探询体积能够占据的最大值bmax为止；14)将所述当前探询体积的X坐标设置为零、即a＝0，将所述当前探询体积的Y坐标设置为零、即b＝0，并且使当前位置体积的Z坐标增加1个体素、即c＝c+1；15)重复步骤8)～14)并且存储Pi和σi的所有值，直到所述当前探询体积的Z坐标的值c等于所述当前探询体积能够占据的最大值cmax为止；16)增大所述当前探询体积的大小，其包括：i.通过使指针增加至下一探询体积、即i＝i+1来选择下一组探询体积，以及ii.将当前探询大小设置为Xi,Yi,Zi；17)重复步骤6)～16)，直到选择了所有探询体积并且计算并存储了Pi和σi的所有值为止；18)选择用以匹配的一个或多个所选择属性；19)针对各探询体积计算λi；20)选择λi为最小值的探询体积，其中所选择的探询体积是代表性单元体积即REV的大小和位置；以及21)计算所述多孔介质的属性。
还提供了用于进行这些方法的计算机化系统、计算机可读介质和程序。
应当理解，以上的一般说明和以下的详细说明仅是示例性和解释性的，并且仅意图提供本发明的进一步解释。
包含在本申请中并构成本申请一部分的附图示出了本发明的一部分实施例，并和说明书一起用来解释本发明的原理。附图未必是按比例绘制的。附图中的相同附图标记指代各视图中的相同元件。
附图说明
图1是示出用于使用M-REV程序来识别REV的先前实践的流程图。
图2示出用于通过从代表性区域中选择REV来识别REV的另一先前采样方案。
图3示出用于使用孔隙率REV方法来识别REV的另一先前采样方案。
图4是示出REV的先前定义的所测量属性vs样本体积的曲线图。
图5A和5B示出根据本申请的示例的、利用具有经由小的限制性导管所连接的大导管的管的流体流动系统所建模的孔隙空间连通性中的子样本选择。
图6示出根据本申请的示例的、用以基于从统计上鉴定的子体积来估计REV的方法的流程图。
图7示出根据本申请的示例的、包括与样本和探询体积有关的术语的定义的样本体积和探询体积。
图8A和8B示出根据本申请的示例的具有明显方向流体流动特性的所建模流体流动系统中的子体积选择，其中，图8A是流体流动导管的俯视图，并且图8B是沿着图8A的线8B-8B所截取的流体流动导管的截面图。
图9示出根据本申请的示例的探询体积的数字切片。
图10是示出根据本申请的示例的用以估计REV的方法的流程图，其中该方法还包括以下步骤：使网格定向成流动；测试与多个属性的子体积适当性；以及有方法地移动经过子体积。
图11A和11B示出根据本申请的示例的更加复杂的所建模流体流动系统中的子体积选择，其中，在图11A中重新对齐笛卡尔(Cartesian)网格，并且图11B是沿图11A的线11B-11B所截取的截面图。
图12示出根据本申请的示例的表示具有基本异质的特征的天然岩石样本的分割体积。
图13示出根据本申请的示例的、表示具有与图12所示的异质结构相比变少的异质结构的天然岩石样本的分割体积。
图14是概述根据本申请的示例的一个实施例的详细流程图。
图15A和15B示出根据本申请的示例的、分别关于图11A和11B中所建模的流体流动系统的针对表面积/体积比和孔隙率的标准偏差分布，其中探询体积的大小与对应于样本整体内的周期性的单元晶格相匹配。
图16A～16E示出根据本申请的示例的、针对图11A～11B中所建模的流体流动系统中的不同的探询体积大小的流体流动系统内的孔隙率的标准偏差。
图17A～17E示出根据本申请的示例的、针对图11A～11B中所建模的流体流动系统中的不同的探询体积大小的流体流动系统内的孔隙空间的表面积与体积之比的标准偏差。
图18A～18B示出根据本申请的示例的、针对450×450×450的探询体积的图13的岩石样本的孔隙率(图18A)和表面积与体积之比(图18B)的目标函数的标准偏差。
图19A～19B示出根据本申请的示例的、针对300×300×300的探询体积的图13的岩石样本的孔隙率(图19A)和表面积与体积之比(图19B)的目标函数的标准偏差。
图20A～20B示出根据本申请的示例的、针对200×200×200的探询体积的图13的岩石样本的孔隙率(图20A)和表面积与体积之比(图20B)的目标函数的标准偏差。
图21A～21B示出根据本申请的示例的、针对450×450×450的探询体积的图12的岩石样本的孔隙率(图21A)和表面积与体积之比(图21B)的目标函数的标准偏差。
图22A～22B示出根据本申请的示例的、针对300×300×300的探询体积的图12的岩石样本的孔隙率(图22A)和表面积与体积之比(图22B)的目标函数的标准偏差。
图23A～23B示出根据本申请的示例的、针对200×200×200的探询体积的图12的岩石样本的孔隙率(图23A)和表面积与体积之比(图23B)的目标函数的标准偏差。
图24是达西流动的示意例示。
图25示出根据本申请的示例的、针对子样本维度为95×95×95(灰色三角形)、190×190×190(灰色圆形)和285×285×285(灰色十字形)的Fountainebleau岩石样本的绝对渗透率(mD)vs(作为0～1.0的小数值的)孔隙率的标绘图中的三个poro-perm趋势，其中该标绘图针对趋势线(包括中空菱形符号数据点的“UL”粗实线)包括500×500×500的原始样本大小的poro-perm值，并且黑色符号是通过以函数表面积/体积和孔隙率这两者作为目标所得出的最优选择。这两条线分别是通过针对岩石所进行的实验而得到的下限和上限。
图26示出根据本申请的示例的、针对子样本维度为300×300×300(灰色十字形)、200×200×200(灰色圆形)和100×100×100(灰色三角形)的非固结砂岩样本的绝对渗透率(mD)vs孔隙率的标绘图中的孔隙率/渗透率相关关系(poro-perm)趋势。“1_100”和“1_200”以及“2_100”和“2_200”这两组数据都是针对两个不同样本(样本1和2)。这两个样本非常相似。
图27示出根据本申请的示例的、针对与图25的样本相比孔隙率变低的Fountainebleau样本的绝对渗透率(mD)vs孔隙率的标绘图中的poro-perm趋势曲线，其中该标绘图包括针对子样本维度为190×190×190(灰色三角形)、285×285×285(灰色圆形)和380×380×380(灰色十字形)的poro-perm趋势，并且针对趋势线UL(与中空菱形符号交叉的粗灰色线)包括针对500×500×500的原始样本大小的poro-perm值，并且黑色符号是通过以函数表面积/体积和孔隙率这两者作为目标所得到的最优选择。“实验室下限(Lowr Lab)”曲线是通过针对这些岩石进行的实验所得到的下限。
图28A～28H示出根据本申请的示例的、分别针对图25和图27中所涉及的两个不同Fontainebleau岩石的、孔隙率的标准偏差(图28A)、表面积与体积之比(图28B)、孔隙率的相同分布的方差(V)(图28C)、表面积与体积之比的方差(图28D)、孔隙率的偏度(图28E)、偏度的方差(图28F)、孔隙率的峰度(图28G)和峰度的方差(图28H)vs子体积样本维度(大小)的分布的平均比值(A)。图25中所涉及的Fontainebleau岩石由图28A～28H中的黑色圆形所定义的曲线来表示，图27中所涉及的岩石由图28A～28H中的灰色圆形所定义的曲线来表示。在图28A～28H中，与插图图例以及图25和27所给出的编号相对应地对这些标绘图进行编号。
图29A～29H示出根据本申请的一个示例的、分别针对两个不同的碳酸盐岩石的孔隙率标准偏差的分布(图29A)、表面积与体积之比(图29B)、孔隙率的相同分布的方差(V)(图29C)、表面积与体积之比的方差(图29D)、孔隙率的偏度(图29E)、偏度的方差(图29F)、孔隙率的峰度(图29G)和峰度的方差(图29H)vs子体积样本维度(大小)的平均比值(A)，其中在这些标绘图中，利用灰色圆形表示一个样本并且利用黑色圆形表示另一样本。在图29A～29H中，与插图图例以及图29I和29J所给出的编号相对应地对这些标绘图进行编号。
图29I～29J示出根据本申请的示例的针对图29A～29H中所涉及的两个不同碳酸盐岩石的poro-perm趋势。图29I包括针对子样本维度为95×95×95、190×190×190和285×285×285的图29A～29H中的灰色圆形所表示的样本的poro-perm趋势，其中该poro-perm趋势针对趋势线D1(粗灰色线、中空菱形符号)包括500×500×500的原始样本维度大小的poro-perm值，并且黑色符号是通过以函数表面积/体积和孔隙率这两者作为目标所得到的最优选择。图29J包括针对子样本维度为190×190×190、285×285×285和380×380×380的图29A～29H中由黑色圆形所表示的样本的poro-perm趋势，其中该poro-perm趋势针对趋势线D2(粗灰色线、中空菱形符号)包括500×500×500的原始样本维度大小的poro-perm值，并且黑色符号是通过以函数表面积/体积和孔隙率这两者作为目标所得到的最优选择。
图30A～30H示出根据本申请的示例的、分别针对两个不同的相对均质岩石的孔隙率的标准偏差的分布(图30A)、表面积与体积之比(图30B)、孔隙率的相同分布的方差(V)(图30C)、表面积与体积之比的方差(图30D)、孔隙率的偏度(图30E)、偏度的方差(图30F)、孔隙率的峰度(图30G)和峰度的方差(图30H)vs子体积样本维度(大小)的平均比值(A)，其中在这些标绘图中，利用灰色圆形来表示一个样本并且利用黑色圆形来表示其它样本。在图30A～30H中，与插图图例中所给出的编号相对应地对这些标绘图进行编号。与先前的图相同，这些图示出针对两个不同的均质岩石样本的标绘图。
图31A～31H示出根据本申请的示例的、分别针对两个附加岩石的孔隙率的标准偏差的分布(图31A)、表面积与体积之比(图31B)、孔隙率的相同分布的方差(V)(图31C)、表面积与体积之比的方差(图31D)、孔隙率的偏度(图31E)、偏度的方差(图31F)、孔隙率的峰度(图31G)和峰度的方差(图31H)vs子体积样本维度(大小)的平均比值(A)，其中在这些标绘图中，利用灰色圆形来表示一个样本并且利用黑色圆形来表示其它样本。在图31A～31H中，与插图图例以及图31I和31J中所给出的编号相对应地对这些标绘图进行编号。用作样本的岩石是砂岩(Fontainbleau)。
图31I～31J示出根据本申请的示例的图31A～31H中所涉及的两个不同岩石的poro-perm趋势。图31I包括针对子样本维度为190×190×190、285×285×285和380×380×380的图31A～31H中利用灰色圆形标识的样本的poro-perm趋势，其中该poro-perm趋势针对趋势线D3(粗灰色线、中空棱形符号)包括原始样本维度大小为500×500×500的poro-perm值，并且黑色符号是通过以函数表面积/体积和孔隙率这两者作为目标所得到的最优选择。图31J包括针对子样本维度为95×95×95、190×190×190和285×285×285的图31A～31H中利用黑色圆形所标识的样本的poro-perm趋势，其中该poro-perm趋势针对趋势线D4(粗灰色线、中空棱形符号)包括原始样本维度大小为500×500×500的poro-perm值，并且黑色符号是通过以函数表面积/体积和孔隙率这两者作为目标所得到的最优选择。
图32A～32B、33A～33B和34A～34B分别示出根据本申请的示例的、在原始维度为550×550×550的情况下所分析的砂岩岩石样本的孔隙率的目标函数(图32A、33A、34A)和表面积与体积之比的目标函数(图32B、33B、34B)的标准偏差的分布，其中标准偏差的分布是在子样本为200×200×200的情况下所获得的，并且分割分辨率对于图32A～32B为10X、对于33A～33B为20X并且对于图34A～34B为40X。
图35A～35B和36A～36B分别示出根据本申请的示例的、在原始维度为550×550×550的情况下所分析的砂岩岩石样本的孔隙率的目标函数(图35A、36A)和表面积与体积之比的目标函数(图35B、36B)的标准偏差的分布，其中标准偏差的分布是在子样本为200×200×200的情况下所获得的，并且分割分辨率对于图35A～35B为4X并且对于图36A～36B为10X。
图37示出根据本申请的示例的、将多孔介质的三维(3D)扫描成像分析与应用于多孔介质的3D数字表示的方法集成的系统。
具体实施方式
本申请部分涉及用以估计诸如岩石等的多孔介质的样本中的代表性单元体积(REV)的高效方法，其中与先前方法所提供的子体积相比，所选择的子体积是单元体积的更好近似。
本发明还至少部分涉及用于通过使用较小的子样本来表现诸如储集岩等的多孔样本的特性的方法，其中该较小的子样本具有相同或非常相似的所选择特性以及这些所选择特性在经由样本的预期流体流动的方向上的相同或非常相似的变化。如果样本过大，则这些样本可能会损害计算机的存储器，并且可能需要过多的计算机时间来完成计算。因此，本发明部分涉及挑选将代表整个样本的子样本所用的REV、因而可以减少计算时间并且不会损坏计算机存储器的方法。REV在原始样本内具有样本大小和特定位置。REV例如可以是原始样本内的子样本的物理大小和位置、或者REV可以是原始样本的表示内的子样本的数字大小和位置。该方法产生样本内该位置处的、与较大样本的诸如孔隙率和绝对渗透率等的关注多孔介质特性最匹配的子样本。
可以针对多孔介质的样本的数字表示进行本发明的方法。可以通过首先生成样本的计算机断层成像X射线图像、然后对该数字表示进行分割以将各体素识别为颗粒或孔隙，来产生多孔介质的样本的数字表示。然后，可以通过选择施加有压力的入口面来选择样本的主流动方向以进行后续的岩心分析(RCAL)测量和特殊岩心分析(SCAL)测量。针对子样本垂直于流动方向的的各切片，计算诸如孔隙率、表面积与孔隙体积之比(这里还标记为表面积/体积，并且作为孔隙和立体空间之间的边界的长度(2d)或面积(3d)与孔隙空间的面积(2d)或体积(3d)之比来计算)、其它样本属性或它们的组合等的属性，以使得针对该子样本获得仅依赖于该流动方向坐标的属性。对于这种属性函数f(z)，可以通过如下方程式来计算相对于整个样本的平均值fV的标准偏差(数值)：
&sigma; = &Sigma; i = 1 N z ( f ( z i ) - f V ) 2 N Z f V .
如果上述方程式中的σ是小的数值，则样本的函数f相对于在大的原始域(fV)中被评价的相同函数偏差了小的量，因而函数f是沿着主流动方向的该函数的良好表示(因为函数的变化在该方向上小)。在理想情况下(即，对于无限大的介质)，σ的值变为零。将初始子样本维度选择成接近原始样本的大小。针对子样本位置i计算相对于整个样本的平均值fV的标准偏差。注意，在该过程中，穷尽地使用函数f中所包含的“信息”：对于各子样本，沿流动方向提取统计信息。在一些先前专利中，对于各子样本仅使用平均体积。然后，使子样本移动到原始样本内的所有可能位置x_i处，并且针对各位置计算标准偏差。这样得到利用f描述的所选择属性的标准偏差S_i的分布T。所有子样本之间的频率定义事件的分布。在以下说明中，将分布的方差(其标准偏差)定义为V，将平均值定义为A，并且将众数定义为M。
子样本维度例如在各侧上减小了一个体素以上或者仅在某个方向上减小，并且针对所有可能的子样本位置计算所选择属性。重复该处理，直到对所有可能的子样本大小进行了评价为止或直到满足停止标准为止。
通过使用标准偏差的分布T的众数M、或者平均值A和方差V来选择REV。T的众数或者平均值和方差是较大样本的特性的良好指标。如果分布T的众数接近较大样本的体积的标准偏差σ、并且分布T的方差小，则大量子样本具有与原始大体积相同的所选择属性的变化(例如，在所选择属性是孔隙率的情况下，此变化为异质性)，因而子体积的长度尺度大到足以表示整个原始体积。在原始大体积的所选择属性的标准偏差小、并且分布的方差V也小的情况下，可以得出以下两个结论：1)体积整体的原始大小相对于利用函数f描述的流动方向上的变化(例如，异质性)足够大(该尺度是相对于原始体积的所选择属性的积分尺度)；以及(2)流动方向上的异质性对于大多数子体积而言也为小，因而这些样本是用以表示较大体积的良好候选。如果所选择属性例如是孔隙率和表面积/体积，则预期与和原始体积相同的变化相匹配的子样本具有接近原始体积的绝对渗透率的绝对渗透率。在原始体积的所选择属性的标准偏差为零并且分布T的方差也为零的极端情况下，这意味着原始大体积是由在流动方向上以周期性方式复制子体积所形成的：在这种情况下，子体积表示利用f描述的特定量的单元体积。子样本的最佳位置即REV是标准偏差为1、并且两个或更多个所选择属性以尽可能接近的方式匹配原始样本整体的标准偏差的位置。如果子样本的所选择属性的标准偏差小于原始样本的所选择属性的标准偏差，则子样本的变化小于原始样本的变化，这意味着子样本丢失了原始样本所具有的一些异质性并且子样本从人造角度更好。如果标准偏差大于原始样本的标准偏差，则子样本与原始样本相比具有更强的异质性，并且应丢弃该子样本。结果，本发明的方法可以识别接近单元大小的最具代表性的子样本，并且还可以判断原始样本的异质性是否过大而导致由于无法应用达西定律而无法使用代表性子样本。
如以上在背景技术部分所论述的，图1～3示出在数字模拟和分析中用以识别应用于表示诸如岩石样本等的多孔材料方面的REV的先前努力。图1示出用于研究直径不断增大的同心球形子样本体积的流程图300。图2示出范围不断增大且转换成三维立方体以试图选择适当REV的同心正方形302a、302b和302c。图3示出建模的体积310，其中各种均匀大小的几组子体积(这里为312a、312b和312c)随机地(但不重叠地)配置在体积310内。这些子体积被示出为具有立方体或长方体的几何形状。体积310和子体积312a、312b和312c分别具有图3所示的维度。例如，体积310的维度为600×600×150并且子体积312a的维度为150×150×150。针对各样本，选择并计算诸如孔隙率和/或渗透率等的参数，并且计算方差或可变性。选择方差限制，并且通过符合该限制、例如相对于平均值的正或负(±)5％来确定REV的适当性。
先前尝试没有产生用于近似最小REV的高效方法，并且没有很好地处理天然岩石或其它多孔样本的异质性。此外，先前尝试没有提供与REV是否适合应用达西定律有关的指导。
图6是本发明的用于更好地处理现实世界样本的异质性的方面的处理的流程图。获得样本的3D数字图像作为分割体积110，其中诸如步骤112所示等，根据该分割体积，在体积整体内推导一个或多个属性值“P”并对这些属性值“P”求平均，以生成例如如图所示指定为<Pvol>或<P>的体积整体的平均属性值。如以下将进一步论述的，孔隙率和表面积/体积是在鉴定REV时要应用的便利的目标函数或属性，然而本发明并不局限于此。在步骤114中，还针对体积整体计算所选择属性的标准偏差σvol(“σ”)。在步骤116中定义一组子体积，并且在步骤118中使这一组子体积移动经过整体体积，其中在各子体积处针对各目标函数计算标准偏差σi。将步骤118的输出与停止标准120进行比较，并且如果不满足停止标准，则在步骤122中调整子体积的大小，并且迭代地重复定义子体积的步骤116、移动经过这些子体积并针对各目标函数计算标准偏差σi的步骤118和与停止标准进行比较的步骤120，直到满足停止标准120为止。停止标准例如可以是子体积的给定大小，其中调整子体积的大小包括连续且渐进地缩小或放大子体积。这里稍后在真实应用的例示内说明具体的适当停止标准。在满足120的停止标准的情况下，方法进入求出最小适合的REV子体积的步骤124。例如，通过将子样本的标准偏差σi与整体体积的标准偏差σvol进行一致性比较来测试适合性。在步骤126中存储REV并且使用该REV以推导在步骤128中输出的关注属性值。本发明的该方面的实践提供了更具代表性的REV的选择、即与整体样本体积的异质性以及一个或多个选择标准属性的平均属性值相匹配的REV的选择。
图4是示意性示出代表性单元体积或REV的先前一般定义的曲线图。相对于样本体积的大小来标绘测量属性。利用曲线320所追踪的测量属性的波动随着样本体积的大小而减小，直到减小至子样本中的该属性的值可被视为整体体积的代表的点为止。在该例示中，这也适用于超出样本体积大小322的区域。REV是观察到该现象的最小样本大小。
尽管图4所示的REV定义是开始所用的理想化模型，但该REV定义在样本为同质且各向同性的情况下最适合。通常不是这种情况。考虑例如在图5A和5B中进行建模的情形。在该示例中，在附图中轮廓描画为立方体体积的体积130具有贯穿该体积130的管131。该管是多孔空间，并且具有从大导管134到小导管132的限制的多个直径。针对表面积/体积(而不是针对孔隙率)包含代表性结构的单元晶格在这种情况下是管131的、包括从大导管134到小导管132的转变的区段。在这种情况下，体积130的整体样本的表面积/体积具有与单元晶格(子体积136)的表面积/体积相同的值。实际上，整体体积通过在流动方向上重复的单元晶格的整数次复制而完成。如果如图5A那样、REV的分析中所使用的探询体积136恰好是单元晶格的体积、并且目标函数是表面积/体积，则提供最佳选择。因而，标准偏差σi与沿流动方向的体积整体130的表面积/体积的标准偏差σvol严格相匹配的子体积136与单元晶格的体积相同。
如果单元体积小于单元晶格(参见图5B的子体积138)，则最终结果将是裁剪管131的一部分从而以尽可能接近的方式匹配体积整体的表面积/体积的标准偏差的体积。
现有技术的用以求出受约束的REV以随机地或同心地研究所选点的努力没有很好地倾向于发现或识别这种单元晶格。因而，本发明的一些实践中的另一特征是顺次且递增地移动经过整个样本体积的方法性子样本，其中该整个样本体积具有大小以递增方式改变的子体积。在以下针对图10的论述中介绍了该方面，并且通过使用诸如图7所示等的笛卡尔网格140以利用坐标a,b,c和大小(Xs,Ys,Zs)来设置样本体积142以及移动迭代设置大小(Xi,Yi,Zi)的探询体积144通过样本体积142，来促进该方面。
返回图5B，子体积138最多使其自己的表面积/体积的标准偏差σi和体积整体的表面积/体积的标准偏差σvol之间的差最小。注意，对于该示例，定义孔隙率的标准偏差与试样整体的孔隙率的标准偏差相匹配的子体积138没有意义。使本发明的一些实践受益的另一特征是使用多个目标函数或属性来识别在模拟和属性偏差的宽范围内有用的更具鲁棒性的REV。例如，利用相匹配或最优组合来进行约束以满足孔隙率和表面积/体积这两者将产生更有用的REV。在以下针对图10的进一步论述中将处理该情况。
图5A和5B介绍另一问题、即流动方向的重要性和相对于固定流动方向的REV的定义。在图8A和8B中更特别地例示出该情况。如图8A所示，平行导管152的阵列横跨样本体积150。图8B是沿图8A的线8B-8B所截取的图8A的截面图。子样本154是具有孔隙率的标准偏差和表面积/体积的标准偏差这两者作为目标函数的所识别的REV。包括流体输送属性的重要属性在诸如天然岩石样本等的多孔材料中是各向异性，这很常见。即，属性具有定向性。对齐笛卡尔网格以考虑到该情况便于求解真实的REV。返回图7，使流动方向与例如Z轴对齐将便于该求解。分割体积的视觉检查通过辨识裂缝或提供孔隙空间连通性的导管之间的图案可足以对齐网格、或者可建议通过孔隙不对称性来进行该对齐。可选地，可以通过根据样本或子样本初步推导估计值来确定属性的各向异性。此外，流动方向可以固定作为相对储集层中的岩心的位置的约束。所示示例具有固定流动方向，但应当理解，针对不同流动方向的组合实现本发明也可以是有用的。
图9示出如下笛卡尔网格160的应用，其中该笛卡尔网格160与流动方向(箭头162)对齐，并且穿过以与流动方向垂直的方式截取的探询体积166的数字切片164行进。探询体积166和数字切片164包括个体体素168。顺次处理以与流动方向垂直的方式截取的数字切片不仅便于计算探询体积的标准偏差σi，而且还可应用于计算诸如与图6的步骤112和114的论述有关等的、样本整体的目标函数或属性<P1>、<Pn>以及标准偏差σ1、σn的平均值。
用于进行与流动方向对齐的REV选择所用的处理是使以下的图10所采用并例示的本发明的一些实践受益的另一特征。
图10呈现从功能上示出包括以上介绍的特征的本发明的实施例的流程图170并且将更详细地进行论述。用于获得分割体积的步骤172从通过x射线计算机断层成像、聚焦离子束扫描电子显微镜、磁性共振成像、同步成像或者其它显微层析或显微放射工艺等所制作的天然岩石样本或其它多孔材料的3D灰度图像开始。用于制作根据本发明的方法可使用的图像的适合CT扫描器的示例例如包括Xradia,Inc.(美国加利福尼亚州普莱森顿市)所制造的诸如MicroXCT-200和Ultra XRM-L200CT等的3D断层成像x射线透射显微镜。灰度图像在被分割成表示孔隙空间和一个或多个固相的多个相(诸如颗粒和可能的一个或多个基体相)之前，例如可以经过滤波或其它预处理。并且，初始分割例如可以具有用以呈现最初成像的材料样本的更好表示的后处理步骤。
然而，如以上所论述的，可以最佳地理解样本的属性和行为所依据的许多模拟和推导在计算量和存储量方面极大，并且对于针对样本整体进行这些模拟和推导既不高效也不可行。因而，估计最小REV非常实用。
如图10所示，在172中所获得的分割体积面向与流动方向对齐的笛卡尔网格(步骤174)。使Z轴与通过在表示样本的分割体积中从视觉上检查孔隙空间连通性是明显的流动方向对齐很便利。
可以使用以与流动方向垂直的方式截取的连续垂直切片、针对多个目标函数或属性Pi～Pn来在步骤176中产生平均值<Pi>、<Pn>并且在步骤178中产生标准偏差σ1、σn。在这方面参考以上针对图9的论述。优选地，目标函数从计算量上不具有挑战性，并且还选择这些目标函数以提供输入的多样性，从而得到对于从计算量和存储量上更具挑战性的应用中的宽范围的模拟和属性偏差而言有用的鲁棒REV估计。孔隙率(φ)和孔隙空间的表面积与体积之比是良好候选。如这里所使用的，简单地将孔隙率(φ)计算为分配至数字切片的孔隙空间的体素数除以该切片中的总体素数，并且孔隙率(φ)提供基本目标函数。计算孔隙空间与固相和基体相之间的界面的表面积并且将该表面积除以数字切片的孔隙空间的总面积提供了第二个有用的数字函数。对于许多应用，这两个属性良好地满足合适目标函数的期望标准，但是受益于本发明的人员以及本领域技术人员应当理解可以替换和/或添加其它属性。对于各属性，对数字切片的值求平均以建立仅依赖于流动方向的目标属性的函数，并且对样本体积整体的值进行评价。另一选项是以某种方式向目标函数应用滤波器以沿固定流动方向在特定位置处修改该目标函数。例如，在入口或出口处，子样本可以期望比原生岩石的孔隙率大的孔隙率。
在步骤180中定义子样本或子体积的大小。然后，使具有所定义的大小的子体积传播经过整个体积，从而完成定义子体积的步骤。参考图4所示的REV的基本定义，这些子体积可以从非常小的大小开始并且逐步增大，或者这些子体积可以从大的大小开始并且向较小的子样本变化。在定义了子体积的连续网格的情况下，步骤182有方法地移动经过子样本，并且针对所选择的各目标函数计算并存储标准偏差σ1、σn。步骤184轮询是否满足停止标准，并且判断是从可用子样本中鉴定并选择REV还是替代为增大候选池。停止标准可以如达到预选大小那样简单、或者例如可以基于针对与前一组子样本有关系的最后一组子样本所计算出的方差(V)和平均值(A)中的一个或多个的分析。如果不满足停止标准，则通过在返回至步骤180之前在步骤186中进行迭代do循环并递增地调整子体积的大小、并且在整个样本中定义这种大小的子体积的连续网格，来增大子样本候选池。再次计算并存储标准偏差并且再次轮询停止标准184。在满足停止标准的情况下，诸如步骤186所示等，REV选择通过首先识别子体积的标准偏差满意地与样本体积的标准偏差相匹配的所有子体积而继续。在例如孔隙率和孔隙空间的表面积与体积之比等的多个目标函数的情况下，可能没有样本提供合适的匹配。在这种情况下，可能期望将两个函数组合并应用最小化过程以选择以尽可能接近的方式匹配所有功能的子样本。针对各步骤188定位来自该匹配子样本池的最小子样本，并且在步骤190中，将该最小子样本用于模拟或推导关注属性值、例如需要较大的存储器和/或计算要求的属性值。
图11A和11B示出样本体积192的更加复杂的流体流动模型。这里，通过在单元晶格194的所有方向上进行周期复制来得到分析所用的体积。参考图5A和5B。可以通过在流动的横向方向上交错的、从大196向着小198的转变来得到单元晶格。这里，重新对齐笛卡尔网格200以定向成流动方向，并且与流动方向垂直的单元晶格内的各X-Y平面在孔隙率和孔隙空间的表面积与体积之比方面具有相同值，因而这些量在流动方向上的标准偏差为零(无变化)。由于样本整体是通过单元晶格的整数次复制形成的，因此该样本整体在单元晶格的孔隙率和表面积/体积方面具有相同值。
与图5A和5B的示例相同，如果探询体积具有维度相同的单元晶格，则由于周期性，因此样本整体内的该维度的所有可能体积将在孔隙率和表面积/体积方面具有相同值，其中标准偏差为零。在图11A和11B(实际为以上所论述的图8A和8B所示的相同系统的不同视图)中进行建模的这种情况下，表面积/体积或孔隙率的标准偏差的分布将为零，其中方差也为零。例如，参见分别示出针对孔隙率和孔隙空间的表面积与体积之比这两者的相同标准偏差曲线的图15A和15B。
在探询体积的维度相对于单元晶格开始改变的情况下，标准偏差的分布示出特定维度在整个区域内不再是周期性的：方差较大的分布是与单元晶格相比探询体积较小的分布。在这种情况下，由于在流动方向上孔隙率或表面积/体积的变化将较大，因此预期到大的变化。图16A～16E示出针对不同的探询体积大小的孔隙率的标准偏差的分布。图11B所示的单元晶格的维度为80×80×40。图16A涉及20×20×10的探询体积。图16B涉及40×40×20的探询体积，并且图16C～16E分别涉及79×79×39、81×81×41和120×120×60的探询大小。图17A～17E分别涉及针对与图16A～16E相对应的探询体积的大小的表面积/体积的分布。
通过以上示例，应当清楚，在分布的众数接近零并且其方差也小的情况下，关于特定目标函数(孔隙率或表面积/体积)，探询体积是样本整体内的准周期结构。
参见图12和13，将相同的分析应用于真实岩石是进一步有用的。在图12中，与示出具有较大异质性的样本的图13相比，样本呈现出非常小的异质性。这两个岩石的样本维度为500×500×500。对于这些岩石各自，基于大小不断改变的探询体积，针对孔隙率和孔隙空间的表面积/体积比这两个不同的目标函数，推导出三个分布。在分别表示针对体积大小为450×450×450、300×300×300和200×200×200的目标函数对的图18A～18B、19A～19B和20A～20B中展示参见图12的异质性较小的岩石的分布。在呈现目标函数对的图18A～18B、19A～19B和20A～20B的曲线图中，在图18A、19A和20A中示出孔隙率，并且在图18B、19B和20B中示出孔隙空间的表面积/体积比。在图21A～21B、22A～22B和23A～23B中，针对图12所示的异质性较大的岩石应用相同大小的样本和呈现。各曲线图将样本整体的特定属性<Pvol>的标准偏差呈现为点并且呈现探询子样本的标准偏差分布σi。
因此，很明显，随着子体积的维度减小，分布的方差增大并且其众数开始在较大值的范围内移动。这意味着从统计上预期到维度较小的子样本内的目标函数的沿着流动方向的变化与原始体积相比更大。在这两种岩石中，如果子岩石的维度非常接近原始岩石的维度、或者对于所选的子岩石的维度而言目标函数在流动方向上的异质性小，则分布的众数非常接近原始岩石的标准偏差的值。该后者情况用于异质性较小的岩石(例如，图12)。在图13的岩石的情况下，可以看出，分布的方差非常大、并且众数与大小已为300×300×300的原始岩石整体的值相差很大。
图14是示出本发明的实施例的流程图。关于本发明的本实施例的该具体实施方式部分，使用以下定义。
1)流动方向垂直于X-Y面
2)Xs＝以体素为单位的样本的宽度
3)Ys＝以体素为单位的样本的高度
4)Zs＝以体素为单位的样本的深度
5)所选择属性可以是Φ、Sv等
6)i＝指向第i个探询体积的指针
7)imax＝探询体积的数量
8)Xi＝以体素为单位的探询体积i的宽度
9)Yi＝以体素为单位的探询体积i的高度
10)Zi＝以体素为单位的探询体积i的深度
11)Xmin,Ymin,Zmin＝探询体积的最小维度
12)Xmax,Ymax,Zmax＝探询体积的最大维度
13)a,b,c＝探询体积的坐标。坐标a,b,c分别是如图6所示的探询体积的左上角的X、Y和Z坐标。
14)Ps(a,b,c)＝位置a,b,c处的整个样本的切片的所选择属性
15)σs＝该组所选择属性Ps(a,b,c)～Ps(a,b,c+Zi)的标准偏差
16)Pi(a,b,c)＝位置a,b,c处的探询体积i的切片的所选择属性
17)σi＝相对于整个样本的一组所选择属性Pi(a,b,c)～Pi(a,b,c+Zi)的标准偏差
18)其中：μ＝所有σj的平均值、即分布(A)的平均值；σs是样本整体的标准偏差、或者在分布的最小值大于原始样本的值的情况下是该最小值。λ的索引i用于特定目标函数、例如孔隙率。如果存在多个目标函数，则可以考虑λi的叠加(或组合)，其中I是目标函数的索引。
本发明的该例示示例可以使用诸如包括以下步骤等的以上组合论述的许多特征，其中放置在括号内的数字是针对在图14中标识出的相关处理流程图框的参考。
1)可以将诸如储集岩等的多孔介质的分割后的三维图像加载至用于处理图像并计算岩石属性的计算机系统(10)。
i.可以使用本领域技术人员所使用的任何分割技术来对分割后的三维图像进行分割。这里可以使用美国专利8,170,799、美国专利8,155,377、美国专利8,085,974、8,081,802和美国专利8,081,796所提及的一个或多个分割技术，并且这些专利通过引用全部包含于此。分割后的三维图像可以包括各自可以分配有灰度值的体素，其中各值表示体素的相对浓度。
ii.可以利用来自计算机断层成像x射线扫描器的原始图像来产生分割后的三维图像，然后利用适当的软件程序对该分割后的三维图像进行分割以将体素分类为颗粒或孔隙等。
2)随后，使用分割后的三维图像来进行用以估计经由多孔介质的流体的流动的模拟。选择流动方向并将该流动方向定义为Z方向(11)。
3)定义探询体积的大小。在图6中示出该命名法的详情。
i.探询体积是维度为Xi、Yi和Zi的、原始分割后的三维图像的子样本。整个样本的维度为Xs、Ys、Zs(12)。
ii.定义探询体积的最大数imax。
iii.设置各探询体积的以体素为单位的维度(Xi,Yi,Zi)。针对i的值(1～imax(12))定义Xi、Yi和Zi。
iv.将i的初始值设置为1(12)。
4)针对探询体积的各切片计算所选择属性Ps(0,0,0)～Ps(0,0,Zs)(13)。在图7中，阴影区域表示5×5×5的探询体积的切片。该切片的角的坐标为(0,0,0)、(0,5,0)、(5,5,0)和(5,0,0)。在该示例中，存在Z＝0～Z＝5的5个切片。最后一个切片的角的坐标为(0,0,5)、(0,5,0)、(5,0,5)和(5,5,5)。
5)计算σs(0,0,0)(14)。
6)设置大小为Xi,Yi,Zi的探询体积在大小为Xs,Ys,Zs的整个样本内可能占据的最大坐标(15)。
i.amax＝Xs–Xi+1
ii.bmax＝Ys–Yi+1
iii.cmax＝Zs–Zi+1
7)将当前探询体积的位置坐标设置为a＝b＝c＝0(16)。
8)针对当前探询体积的切片计算所选择属性Pi(a,b,c)～Pi(a,b,c+Zi)(17)。
i.所选择属性包括孔隙率、表面积与体积之比、相似属性或它们的任何组合。
9)计算σi(a,b,c)(18)。
i.可选地，可以对用于计算σi的值的Ps的平均值进行滤波(19)。
ii.可选地，可以设置Pi的平均值(20)。
10)使探询体积的位置在X方向上移动1个体素(a＝a+1)(21)。
11)重复该方法的步骤8)～10)，存储Pi和σi的所有值，直到当前探询体积的X坐标的值a等于当前探询体积可以占据的最大值amax为止(22)。
12)将当前探询体积的X坐标设置为零(a＝0)，并且使当前位置体积的Y坐标增加1个体素(b＝b+1)(23)。
13)重复该方法的步骤8)～12)，存储Pi和σi的所有值，直到当前探询体积的Y坐标的值b等于当前探询体积可以占据的最大值bmax为止(24)。
14)将当前探询体积的X坐标设置为零(a＝0)，将当前探询体积的Y坐标设置为零(b＝0)，并且使当前位置体积的Z坐标增加1个体素(c＝c+1)(25)。
15)重复该方法的步骤8)～14)，存储Pi和σi的所有值，直到当前探询体积的Z坐标的值c等于当前探询体积可以占据的最大值cmax为止(26)。
16)增大(或减小)当前探询体积的大小(27)。
i.通过使指针增加至下一探询体积(i＝i+1)来选择下一组探询体积。
ii.将当前探询大小设置为Xi,Yi,Zi。
17)重复步骤6)～16)，直到选择了所有的探询体积并且计算并存储了Pi和σi的所有值为止(28)。
18)选择用以匹配的一个或多个所选择属性(29)。
19)针对各探询体积计算λi(30)。
20)选择具有最小λi值的探询体积。这是REV的大小和位置(31)。
21)计算多孔介质的期望属性。
i.期望属性可以包括常规岩心分析(RCAL)和特殊岩心分析(SCAL)。RCAL分析包括但不限于多个轴(x,y,z)中的孔隙率(连接、分离、合计)、油母质含量、绝对渗透率。SCAL分析包括但不限于相对渗透率(两相相对渗透率：水-油、气-油或水-气置换)、毛细管压力(主排水、吸液和二次排水循环的各饱和时的毛细管压力值)、颗粒大小分布、电气属性(形成因素、电阻率指数、a、m、n)、弹性属性(Vp.Vs、E、K、G、泊松(Poisson)比)和相似分析。
参考图24，本发明的另一特征是REV针对达西定律的应用的适合性的分析。如上所述，在方程式1中示出经常应用于经由诸如岩石样本等的多孔介质的流动时的达西定律。经常通过经由采用Gary的分解的Navier Stokes方程式(动量方程式)应用达西定律来获得渗透率。然而，如以上所论述的，该方法依赖于假定是相对于平均积分尺度“表现良好”的平均量(在这种情况下为压力)。不幸地，在与平均长度标度相当的长度标度内快速改变的压力信号无法表示这些应用提供可靠结果的该长度内的压力。在根据上述内容选择REV的情况下，仍存在如下可能性：可能存在子样本内的孔隙率的变化，从而得出达西定律无效或可能出错的假设。此外，压力梯度可以沿着流动方向快速改变，这导致无法定义与特定子样本相关联的渗透率。这特别适用于诸如在现实世界岩石形成中所发现的样本等的高度异质性的样本。
因而，本发明的另一特征是如下的高效方法，其中该高效方法用以量化岩石的数字表征的好(差)程度以及通过达西定律对流体流动的描述的精确程度、即以鲁棒且高效的方式预测孔隙率/渗透率相关关系(“poro-perm”)趋势由于数字子样本变得过小而崩溃。图25涉及以渗透率vs孔隙率的记录的交绘图示出poro-perm趋势的该问题。
这里，原生岩石是有据可查的Fontainebleau岩石的样本，并且原始数字样本的维度为500×500×500。根据数字样本整体所推导出的poro-perm值是大的中空菱形，并且恰好在针对Fontainebleau岩石所示出的“实验室上限(Upperlab)”实验趋势(粗灰色线“UL”)上。“LL”是下限。这证明了原始大小足够大而具有正确的poro-perm关系，从而确认该原始大小是RV(代表性体积)。得知是否可以追踪到将初始岩石整体子分割成较小样本的趋势poro-perm也可以是有用的。如果较小的样本大到足以被视为RV，则这些样本将正确地符合利用中空菱形所示出的poro-perm趋势。目前的问题是一个子样本变得小于REV、换句话说子样本不再是代表性体积的点。灰色十字形符号和灰色圆形符号分别是根据维度为285×285×285和190×190×190的子样本所推导出的poro-perm趋势。通过这些维度，满足“实验室上限”实验趋势。然而，如维度为95×95×95的趋势(灰度三角形)所示，该趋势针对～100x100x100的维度崩溃。将针对分别由黑色圆形符号和黑色十字形符号所示的、大小为190×190×190和285×285×285的子样本分别示出的最佳poro-perm值进行比较。
图26是示出来自作为非固结砂岩的另一种岩石的示例，其示出针对子样本维度分别为300×300×300、200×200×200和100×100×100的由灰色十字形、灰色圆形和灰色三角形所示的数据点的研究表明该岩石的单元体积显然等于或小于100×100×100。
图27示出来自与图25的样本相比孔隙率变低的Fontainebleau样本的结果。原始样本大小为500×500×500的poro-perm值(中空菱形)在“实验室上限”实验曲线(上方粗灰色线“UL”)的上方，以使得大小是代表性大小。“LL”是下限。在poro-perm相关关系在大小约为300×300×300处崩溃的情况下，在这方面参考针对样本大小为285×285×285的灰色圆形所示的值，并且分别与针对样本大小为190×190×190和380×380×380的灰色三角形和灰色十字形所示的趋势进行比较。
因而，图25～27的孔隙率-渗透率交叉绘图示出非常不同的行为。显然，能够精确且高效地预测该行为很有用。
图28A～28H分别示出针对图25和27所涉及的两个Fontainebleau岩石的、孔隙率的标准偏差的分布(图28A)、表面积/体积比的标准偏差的分布(图28B)、孔隙率的相同分布的方差(V)(图28C)和表面积/体积比的方差(图28D)vs子样本维度(大小)的平均值(A)。还对偏度和峰度进行评价，其中在图28E～28H示出这些结果。在图25、26和27中，黑色菱形符号、黑色圆形符号和黑色十字形符号是通过以函数表面积/体积和孔隙率这两者作为目标的工具所得到的最佳选择。在图28A～28H中，由黑色圆形定义的线用于图25的孔隙率较大的样本，并且由灰色圆形定义的线用于图27的孔隙率较小的样本。根据数据的这些趋势，以下变得显而易见：分布的这两个平均值都随着子岩石维度的不断增大而不断减小。对于特定大小，平均值随子岩石的维度的变化率变小。对于190的大小(对于高孔隙率岩石、黑色线)和380的大小(对于低孔隙率岩石、灰色线)，正发生该情况。该分布的方差也不断减小，并且针对维度相同的子岩石达到非常小的值。高阶的动量给出了在子样本的大小足够大的情况下的分布的对称性的指示：该分布变为高斯状。
回忆分布的平均值和标准偏差(或方差)的含义，其中：与在分布为高斯的情况下的分布的众数相同的平均值给出分布关于“零”的“位置”，并且方差解释其关于平均值的分布。图11A～11B和15A～15B例如表明以零为中心且方差为零的分布意味着所选择的子体积的完美周期性。因而，平均值和方差是岩石整体内的子体积的目标函数(孔隙率或表面积/体积)的“周期性”的度量。
在子岩石的维度从原始大小起缓慢地不断减小(在极限情况下，针对每次减小1个体素)的情况下，该分布发生不同的情况：首先，方差开始增大；其次，在大小减小了特定阈值以上的情况下，分布的平均值开始改变，并且分布关于在原生岩石整体中评价的目标函数的值变为不对称(平均值增大，这意味着目标函数在流动方向上的较大的变化)。基本上，分布正向原始位置的右侧移动(参见这里针对两个不同岩石所参考的先前标绘图)。
显而易见，对于给出平均值的偏移的子岩石的维度，相关关系poro-perm崩溃。由于在平均值大并且方差也大的情况下、挑选孔隙率的变化和关于该变化的原始值的表面积/体积大的子样本的可能性高，因此这是有意义的。注意，在该平均值具有与原始样本整体的平均值相同的值的情况下，挑选孔隙率和表面积/体积的变化相同的子样本的可能性高，但这并不暗示值相同的孔隙率或表面积/体积(因而值相同的渗透率)。换句话说，仍可能存在poro-perm趋势。为了进一步表明这些特征，以下提供针对碳酸盐和砂岩的其它情况。
作为进一步证明，例如，图29A～29H示出针对两个不同碳酸盐样本的本发明的方法的应用。对于这两个碳酸盐，认为结果与先前岩石的情况相似。也就是说，在针对特定大小、孔隙率和标准偏差的平均值停止而减小(或以缓慢速率减小)的情况下，方差也小并且可以预期到良好的趋势poro-perm。在图29A～29H中，由黑色圆形定义的线用于孔隙率较大的碳酸盐样本，并且由灰色圆形定义的线用于孔隙率较小的碳酸盐样本。这两个碳酸盐样本在孔隙率、特别是渗透率方面不同。孔隙率较低的碳酸盐样本的渗透率低于100mD，并且孔隙率较高的碳酸盐样本的渗透率低于几百mD。在图29I中，灰色三角形符号、灰色圆形符号和灰色十字形符号涉及分别与图29A～29H中的灰色圆形所标识的样本的poro-perm趋势有关的、针对子样本维度95×95×95、190×190×190和285×285×285的poro-perm趋势。在图29J中，灰色三角形符号、灰色圆形符号和灰色十字形符号涉及分别与图29A～29H中的黑色圆形所标识的样本的poro-perm趋势有关的、针对子样本维度190×190×190、285×285×285和380×380×380所提供的poro-perm趋势。注意，在图29I和29J的各标绘图中，黑色三角形符号、黑色圆形符号和黑色十字形符号是通过以函数表面积/体积和孔隙率作为目标的工具所做出的最优选择。
图30A～30H示出这里为如下情况的本发明的其它应用：岩石是相对同质的，并且可以提供从子样本大小100×100×100开始的良好poro-perm趋势。对于该研究，使用两个不同的相对同质的岩石。在图30A～30H中，由黑色圆形定义的线用于孔隙率较大的样本，并且由灰色圆形定义的线是孔隙率较小的样本。在图26中示出图30A～30H的poro-perm。
图31A～31H还提供附加示例，这里为研究在维度为100×100×100和约200×200×200处poro-perm崩溃的两个岩石。在图31A～31H中，由黑色圆形定义的线用于孔隙率较大的样本，并且由灰色圆形定义的线是孔隙率较小的样本。这些样本是砂岩。在图31I中，(与外侧包围标绘图相对应的)灰色三角形符号、(与内侧包围标绘图相对应的)灰色圆形符号和灰色十字形符号涉及与图31A～31H中利用灰色圆形标识的样本的poro-perm趋势有关的、分别针对子样本维度190×190×190、285×285×285和380×380×380的poro-perm趋势。在图31J中，(与最大包围标绘图相对应的)灰色三角形符号、(与中间包围标绘图相对应的)灰色圆形符号和(与最小包围标绘图相对应的)灰色十字形符号涉及与图31A～31H中利用黑色圆形标识的样本的poro-perm趋势有关的、分别针对子样本维度95×95×95、190×190×190和285×285×285的poro-perm趋势。在图31I和31J的各标绘图中，黑色三角形符号、黑色圆形符号和黑色十字形符号是利用以函数表面积/体积和孔隙率为目标的工具所做出的最优选择。
用以使用本发明的另一方式是估计何种分辨率和视场最适合岩石。实际上，子样本的维度可以是固定的(例如，高达400×400×400)，并且发生改变的是扫描的分辨率和视场。通常，可以利用扫描器固定所使用的点数并且可以在不同大小的体积中分配这些点。这样给出了针对岩石的扫描的不同分辨率：对于较小的视场，与大视场相比，分辨率将变大。一个问题是推定何种视场(因而分辨率)适合岩石。可以使用目标函数的标准偏差的分布来解决该问题，其中子样本的维度对于所有视场而言均是固定的。例如，在图32A～23B、33A～33B和34A～34B中，利用原始维度550×550×550来分析砂岩。利用子样本200×200×200来获得标准偏差的分布(右侧为孔隙率、左侧为表面积/体积)。从图32A～23B至33A～33B至34A～34B发生改变的是分割的分辨率，其中分割后的样本的分辨率分别为10X、20X和40X，这些分辨率可能意味着一个体素分别为2、1和0.5微米。根据图33A～33B明显的是，20X的分辨率的分布的方差非常大，并且其平均值比如图32A～32B所示的10X的分辨率的平均值大得多。为了利用20X的分辨率进行工作，在分割中应使用较大的视场。对于该示例而言，10X的分辨率是可接受的。
在图35A～35B和36A～36B所示的下一示例中，分析另一砂岩，并且利用与之前相同的子样本200×200×200来获得标准偏差的分布(右侧为孔隙率、左侧为表面积/体积)。在这种情况下，相对于先前示例的情况，图35A～35B所示的分辨率4X与图36A～36B所示的10X的分辨率相比更加适当，其中10X的分辨率的平均值和方差过大。
参考图37，示出可被配置为进行本方法的系统100。如该示例所示，利用扫描器102生成从源101所获得的多孔介质样本的三维(3D)图像。可以将该扫描器的3D图像输出103传送至具有用于执行3D图像分析以及所指示的数据和模拟分析的程序指令的计算机104，以生成可以发送至诸如显示器、打印机、数据存储介质或它们的组合的一个或多个装置105的样本建模输出/结果。3D图像分析和CFD计算和模拟建模所使用的计算机程序能够作为程序产品存储在与被配置为运行该程序的至少一个处理器104A(例如，CPU)相关联的至少一个计算机可用存储介质104B(例如，硬盘、闪速存储器装置、致密盘、磁带/磁盘或其它介质)上，或者可以存储在计算机处理器可访问的外部计算机可用存储介质(未示出)上。计算机104可以包括至少一个存储器单元104C，其中该存储器单元104C用于存储程序、输入数据和输出数据和其它程序结果或者这些的组合。对于输出显示，装置105例如可以是显示监视器、CRT或其它视觉显示部件(未示出)。计算机104可以包括可作为单个个人计算机或计算机的网络来实现的一个或多个系统计算机。然而，本领域技术人员应当理解，可以在各种计算机系统结构中实践这里所述的各种技术的实现，其中，所述计算机系统结构包括超文本传输协议(HTTP)服务器、手持式装置、多处理器系统、基于微处理器或可编程消费电子产品、网络PC、小型计算机和大型计算机等。系统100的包括扫描器102、计算机104以及输出显示器、打印机和/或数据存储装置/介质105的各单元可以经由硬线、射频通信、电信、因特网连接或其它通信方式彼此连接以进行通信(例如，数据传输等)。
本发明按任意顺序和/或以任意组合包括以下的方面/实施例/特征。
1.本发明涉及一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)针对所述分割体积整体，推导第一目标函数P1的平均属性值<P1>；
c)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差σvol；
d)在所述体积内定义多个子体积；
e)针对各所述子体积，计算所述第一目标函数P1的属性值P相对于平均属性值<P1>的标准偏差σi；
f)求出标准偏差σi与σvol良好地匹配的所有候选代表性子体积；
g)从所述候选中选择并存储代表性子体积；以及
h)使用所述代表性子体积来推导至少一个关注属性值。
2.根据任何前述或以下实施例/特征/方面所述的方法，其中，在所述体积内定义多个子体积还包括：
定义子体积的初始大小；
利用具有所定义的初始大小的子体积填充所述体积整体；以及
对其它子体积的大小进行迭代，并且利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止。
3.根据任何前述或以下实施例/特征/方面所述的方法，其中，对所述大小进行迭代是以小的增量从大到小而进行的。
4.根据任何前述或以下实施例/特征/方面所述的方法，其中，选择并存储代表性体积还包括：求出最小的代表性数字体积。
5.根据任何前述或以下实施例/特征/方面所述的方法，其中，所述停止标准包括所述子体积的给定大小。
6.根据任何前述或以下实施例/特征/方面所述的方法，其中，还包括以下步骤：
将所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向，
其中，
针对所述分割体积整体推导第一目标函数P1的平均属性值<P1>包括对以与所述所定义的流动方向垂直的方式所截取的所述样本体积的多个数字切片进行分析；以及
针对各所述子体积计算所述第一目标函数P1的属性值P相对于平均属性值<P1>的标准偏差σi是相对于所述定义的流动方向而进行的。
7.根据任何前述或以下实施例/特征/方面所述的方法，其中，还包括以下步骤：
针对所述分割体积整体，推导第二目标函数P2的平均属性值<P2>；
针对所述分割体积整体，计算相对于平均属性值<P2>的标准偏差σvol；
在所述体积内定义多个子体积；
针对各所述子体积，计算第二目标函数P2的属性值P相对平均属性值<P2>的标准偏差σi；以及
针对第一目标函数P1和第二目标函数P2的组合，求出标准偏差σi与σvol良好地匹配的所有代表性子体积。
8.根据任何前述或以下实施例/特征/方面所述的方法，其中，所述第一目标函数P1是孔隙率，并且所述第二目标函数P2是所述孔隙空间的表面积与体积之比。
9.根据任何前述或以下实施例/特征/方面所述的方法，其中，还包括以下步骤：在选择之前鉴定候选子体积，包括确定所述候选子体积用来根据达西定律推导流体输送属性的适合性，所述步骤包括：
构建目标函数的标准偏差的分布；
对目标函数的标准偏差的分布的平均值或者可选地任何其它一阶特性以及所述分布的方差、峰度或偏度进行评价；
对一阶以上的动量相对于所述子体积的维度的趋势进行评价；以及
在一阶动量相对于该一阶动量针对在较大子体积上构建的分布的值的变化至少为0.1的情况下、以及/或者在针对方差而言高阶动量高于特定阈值0.1的情况下，停止使所述子体积的维度减小。
10.本发明涉及一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)使所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向；
c)通过对与所述所定义的流动方向垂直的数字切片进行分析，来针对所述分割体积整体推导作为至少第一目标函数P1的一个或多个函数的值；
d)在所述体积内定义多个子体积；
e)关于所述所定义的流动方向，针对各所述子体积计算至少第一目标函数P1的一个或多个函数的值；
f)求出函数识别体积值和子体积值之间的匹配的所有代表性子体积候选；
g)从所述候选中选择代表性体积；
h)存储所述代表性子体积；以及
i)使用所述代表性子体积来进行模拟或推导至少一个关注属性值。
11.本发明涉及一种用于从多孔样本的较大三维数字图像中获得代表性单元体积的高效估计的方法，所述方法包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)针对所述分割体积整体，推导作为至少第一目标函数P1的至少一个函数的值；
c)在所述体积内定义多个子体积，其包括：
定义子体积的初始大小，
利用具有所定义的初始大小的子体积填充所述体积整体，以及
对其它子体积的大小进行迭代，以及利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止；
d)针对各所述子体积，计算作为至少第一目标函数的至少一个函数的值；
e)针对良好地匹配的所述体积和所述子体积的值，求出所有代表性子体积候选；
f)从所述候选中选择并存储代表性子体积；以及
g)使用所述代表性子体积来进行模拟或推导至少一个关注属性值。
12.根据任何前述或以下实施例/特征/方面所述的方法，其中，还包括以下步骤：在选择之前鉴定候选子体积，包括确定所述候选子体积用来根据达西定律推导流体输送属性的适合性，所述步骤包括：
构建目标函数的标准偏差的分布；
对目标函数的标准偏差的分布的平均值或者可选地任何其它一阶特性以及所述分布的方差、峰度或偏度进行评价；
对一阶以上的动量相对于所述子体积的维度的趋势进行评价；以及
在一阶动量相对于所述一阶动量针对在较大子体积上构建的分布的值的变化至少为0.1的情况下、以及/或者在针对方差而言高阶动量高于特定阈值0.1的情况下，停止使所述子体积的维度减小。
13.本发明涉及一种用于从多孔样本的较大三维数字图像中获得代表性单元体积的高效估计的方法，包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)将所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向；
c)使用针对以与所述所定义的流动方向垂直的方式所截取的样本体积的多个数字切片的分析，来针对所述分割体积整体推导第一目标函数P1的平均属性值<P1>；
d)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差；
e)在所述体积内定义多个子体积，其包括：
定义子体积的初始大小，
利用具有所定义的初始大小的子体积填充所述体积整体，以及
对其它子体积的大小从大到小地进行迭代，利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止；
f)关于所述所定义的流动方向，针对各所述子体积计算属性P相对于平均属性值<P1>的标准偏差σi；
g)求出σi与σvol良好地匹配的所有候选代表性子体积；
h)选择最小的候选并且存储所述最小的候选作为代表性单元体积；以及
i)使用所述代表性单元体积来推导至少一个关注属性值。
14.根据任何前述或以下实施例/特征/方面所述的方法，其中，还包括以下步骤：
针对所述分割体积整体，推导第二目标函数P2的平均属性值<P2>；
针对所述分割体积整体，计算相对于平均属性值<P2>的标准偏差σvol；
在所述体积内定义多个子体积；
针对各所述子体积，计算第二目标函数P2相对于平均属性值<P2>的标准偏差σi；以及
针对第一目标函数P1和第二目标函数P2的组合，求出σi与σvol良好地匹配的所有代表性子体积。
15.根据任何前述或以下实施例/特征/方面所述的方法，其中，第一目标函数P是孔隙率，并且第二目标函数P2是所述孔隙空间的表面积与体积之比。
16.根据任何前述或以下实施例/特征/方面所述的方法，其中，还包括以下步骤：在选择之前鉴定候选子体积，包括确定所述候选子体积用来根据达西定律推导流体输送属性的适合性，所述步骤包括：
构建目标函数的标准偏差的分布；
对目标函数的标准偏差的分布的平均值或者可选地任何其它一阶特性以及所述分布的方差、峰度或偏度进行评价；
对一阶以上的动量相对于所述子体积的维度的趋势进行评价；以及
在所述一阶动量相对于所述一阶动量针对在较大子体积上构建的分布的值的变化至少为0.1(或者至少为0.5、1、2、5或任何值)的情况下、以及/或者在针对方差而言高阶动量高于特定阈值0.1(或其它值)的情况下，停止使所述子体积的维度减小。
17.一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：
1)将多孔介质的分割三维图像加载至计算机系统中，
其中，所述分割三维图像包括各自分配有灰度值的体素；
2)选择被定义为Z方向的流动方向；
3)定义探询体积的大小，其中：
i.所述探询体积是维度为Xi、Yi和Zi的原始分割三维图像的子样本，其中整个样本的维度为Xs,Ys,Zs，
ii.定义探询体积的最大值imax，
iii.设置各所述探询体积的以体素为单位的维度(Xi,Yi,Zi)，其中Xi、Yi和Zi是针对值为1～imax的i所定义的，以及
iv.将i的初始值设置为1；
4)针对所述探询体积的各切片来计算所选择属性Ps(0,0,0)～Ps(0,0,Zs)；
5)计算σs(0,0,0)；
6)设置大小为Xi,Yi,Zi的所述探询体积在大小为Xs,Ys,Zs的整个样本内所占据的最大坐标，其中：
i.amax＝Xs–Xi+1，
ii.bmax＝Ys–Yi+1，以及
iii.cmax＝Zs–Zi+1；
7)将当前探询体积的位置坐标设置为a＝b＝c＝0；
8)针对所述当前探询体积的切片计算所选择属性Pi(a,b,c)～Pi(a,b,c+Zi)，
其中，所述所选择属性包括孔隙率、表面积与体积之比、相似属性或它们的任意组合；
9)计算σi(a,b,c)，
i.其中，可选地，对计算σi的值所使用的Pi的值进行滤波，以及
ii.其中，可选地，设置Pi的平均值；
10)使所述探询体积的位置在X方向上移动1个体素、即a＝a+1；
11)重复步骤8)～10)并且存储Pi和σi的所有值，直到所述当前探询体积的X坐标的值a等于所述当前探询体积能够占据的最大值amax为止；
12)将所述当前探询体积的X坐标设置为零即a＝0，并且使当前位置体积的Y坐标增加1个体素、即b＝b+1；
13)重复步骤8)～12)并且存储Pi和σi的所有值，直到所述当前探询体积的Y坐标的值b等于所述当前探询体积能够占据的最大值bmax为止；
14)将所述当前探询体积的X坐标设置为零、即a＝0，将所述当前探询体积的Y坐标设置为零、即b＝0，并且使当前位置体积的Z坐标增加1个体素、即c＝c+1；
15)重复步骤8)～14)并且存储Pi和σi的所有值，直到所述当前探询体积的Z坐标的值c等于所述当前探询体积能够占据的最大值cmax为止；
16)增大所述当前探询体积的大小，其包括：
i.通过使指针增加至下一探询体积、即i＝i+1来选择下一组探询体积，以及
ii.将当前探询大小设置为Xi,Yi,Zi；
17)重复步骤6)～16)，直到选择了所有探询体积并且计算并存储了Pi和σi的所有值为止；
18)选择用以匹配的一个或多个所选择属性；
19)针对各探询体积计算λi；
20)选择λi为最小值的探询体积，其中所选择的探询体积是代表性单元体积即REV的大小和位置；以及
21)计算所述多孔介质的属性。
18.根据任何前述或以下实施例/特征/方面所述的方法，其中，所述分割三维图像是如下产生的：通过利用计算机断层成像x射线扫描器扫描所述样本获得所述样本的图像，并且利用软件程序对所述图像进行分割，以将体素分类为颗粒、孔隙和可选地其它相。
19.根据任何前述或以下实施例/特征/方面所述的方法，其中，所述属性包括常规岩心分析属性即RCAL属性、特殊岩心分析属性即SCAL属性、或这两者。
20.根据任何前述或以下实施例/特征/方面所述的方法，其中，所述RCAL分析属性是沿多个轴的孔隙率、油母质含量、多个轴的绝对渗透率，并且所述SCAL属性是相对渗透率、毛细管压力、颗粒大小分布、电气属性、弹性属性和它们的任意组合。
21.一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的系统，包括：
a)扫描器，其能够产生所述多孔介质的三维数字图像；
b)计算机，其包括至少一个处理器，其中所述至少一个处理器用于执行能够获得表现孔隙空间和至少一个固相的特性的分割体积的计算机程序，
c)与b)相同或不同的计算机，其包括至少一个处理器，其中所述至少一个处理器用于执行能够进行计算的计算机程序，其中所述计算包括：i)针对所述分割体积整体，推导第一目标函数P1的平均属性值<P1>；ii)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差σvol；iii)在所述体积内定义多个子体积；iv)针对各所述子体积，计算第一目标函数P1的属性值P相对于平均属性值<P1>的标准偏差σi；v)求出标准偏差σi与σvol良好地匹配的所有候选代表性子体积；vi)从所述候选中选择并存储代表性子体积；以及vii)使用所述代表性子体积来推导至少一个关注属性值；以及
d)用以显示、打印或存储所述计算的结果的至少一个装置。
22.一种计算机可读介质(例如，非瞬态)上的计算机程序产品，其中，在所述计算机程序产品在计算机化装置中的处理器上执行的情况下，提供用于进行根据前述方法和系统的所示步骤中的一个或多个步骤或所有步骤的计算的方法。
本发明可以包括如句子和/或段落中所记载的以上和/或以下的这些各种特征或实施例的任何组合。这里所公开的特征的任何组合均被视为本发明的一部分并且并不意图局限于可组合的特征。
通过上述说明和所附权利要求书将明白其它特征、方面和优点。此外，在没有背离要求保护的发明的范围的情况下，并非所有特征、方面和优点均需存在于本发明的各实施例中，并且可以单独或与其它特征、方面或优点相组合地出现在各实施例中。

Claims (21)

1.一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)针对所述分割体积整体，推导第一目标函数P1的平均属性值<P1>；
c)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差σvol；
d)在所述体积内定义多个子体积；
e)针对各所述子体积，计算所述第一目标函数P1的属性值P相对于平均属性值<P1>的标准偏差σi；
f)求出标准偏差σi与σvol良好地匹配的所有候选代表性子体积；
g)从所述候选中选择并存储代表性子体积；以及
h)使用所述代表性子体积来推导至少一个关注属性值。
2.根据权利要求1所述的方法，其中，在所述体积内定义多个子体积还包括：
定义子体积的初始大小；
利用具有所定义的初始大小的子体积填充所述体积整体；以及
对其它子体积的大小进行迭代，并且利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止。
3.根据权利要求2所述的方法，其中，对所述大小进行迭代是以小的增量从大到小而进行的。
4.根据权利要求3所述的方法，其中，选择并存储代表性子体积还包括：求出最小的代表性数字体积。
5.根据权利要求4所述的方法，其中，所述停止标准包括所述子体积的给定大小。
6.根据权利要求2所述的方法，其中，还包括以下步骤：
将所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向，
其中，
针对所述分割体积整体推导第一目标函数P1的平均属性值<P1>包括：对以与所述所定义的流动方向垂直的方式所截取的样本体积的多个数字切片进行分析；以及
针对各所述子体积计算所述第一目标函数P1的属性值P相对于平均属性值<P1>的标准偏差σi是相对于所述所定义的流动方向而进行的。
7.根据权利要求6所述的方法，其中，还包括以下步骤：
针对所述分割体积整体，推导第二目标函数P2的平均属性值<P2>；
针对所述分割体积整体，计算相对于平均属性值<P2>的标准偏差σvol；
在所述体积内定义多个子体积；
针对各所述子体积，计算第二目标函数P2的属性值P相对于平均属性值<P2>的标准偏差σi；以及
针对第一目标函数P1和第二目标函数P2的组合，求出标准偏差σi与σvol良好地匹配的所有代表性子体积。
8.根据权利要求7所述的方法，其中，所述第一目标函数P1是孔隙率，并且所述第二目标函数P2是所述孔隙空间的表面积与体积之比。
9.根据权利要求8所述的方法，其中，还包括以下步骤：在选择之前鉴定候选子体积，包括确定所述候选子体积用来根据达西定律推导流体输送属性的适合性，该步骤包括：
构建目标函数的标准偏差的分布；
对目标函数的标准偏差的分布的平均值或者可选地任何其它一阶特性以及所述分布的方差、峰度或偏度进行评价；
对一阶以上的动量相对于所述子体积的维度的趋势进行评价；以及
在一阶动量相对于该一阶动量针对在较大子体积上构建的分布的值的变化至少为0.1的情况下、以及/或者在针对方差而言高阶动量高于特定阈值0.1的情况下，停止使所述子体积的维度减小。
10.一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)使所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向；
c)通过对与所述所定义的流动方向垂直的数字切片进行分析，来针对所述分割体积整体推导作为至少第一目标函数P1的一个或多个函数的值；
d)在所述体积内定义多个子体积；
e)关于所述所定义的流动方向，针对各所述子体积计算至少第一目标函数P1的一个或多个函数的值；
f)求出函数识别体积值和子体积值之间的匹配的所有代表性子体积候选；
g)从所述候选中选择代表性子体积；
h)存储所述代表性子体积；以及
i)使用所述代表性子体积来进行模拟或推导至少一个关注属性值。
11.一种用于从多孔样本的较大三维数字图像中获得代表性单元体积的高效估计的方法，包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)针对所述分割体积整体，推导作为至少第一目标函数P1的至少一个函数的值；
c)在所述体积内定义多个子体积，其包括：
定义子体积的初始大小，
利用具有所定义的初始大小的子体积填充所述体积整体，以及
对其它子体积的大小进行迭代，以及利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止；
d)针对各所述子体积，计算作为至少第一目标函数的至少一个函数的值；
e)针对良好地匹配的所述体积和所述子体积的值，求出所有代表性子体积候选；
f)从所述候选中选择并存储代表性子体积；以及
g)使用所述代表性子体积来进行模拟或推导至少一个关注属性值。
12.根据权利要求11所述的方法，其中，还包括以下步骤：在选择之前鉴定候选子体积，包括确定所述候选子体积用来根据达西定律推导流体输送属性的适合性，该步骤包括：
构建目标函数的标准偏差的分布；
对目标函数的标准偏差的分布的平均值或者可选地任何其它一阶特性以及所述分布的方差、峰度或偏度进行评价；
对一阶以上的动量相对于所述子体积的维度的趋势进行评价；以及
在一阶动量相对于所述一阶动量针对在较大子体积上构建的分布的值的变化至少为0.1的情况下、以及/或者在针对方差而言高阶动量高于特定阈值0.1的情况下，停止使所述子体积的维度减小。
13.一种用于从多孔样本的较大三维数字图像中获得代表性单元体积的高效估计的方法，包括以下步骤：
a)获得表现孔隙空间和至少一个固相的特性的分割体积；
b)将所述分割体积的笛卡尔网格的所选择轴定向为所定义的流动方向；
c)使用针对以与所述所定义的流动方向垂直的方式所截取的样本体积的多个数字切片的分析，来针对所述分割体积整体推导第一目标函数P1的平均属性值<P1>；
d)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差σvol；
e)在所述体积内定义多个子体积，其包括：
定义子体积的初始大小，
利用具有所定义的初始大小的子体积填充所述体积整体，以及
对其它子体积的大小从大到小地进行迭代，利用具有这样大小的子体积填充所述体积整体，并且重复该步骤，直到满足停止标准为止；
f)关于所述所定义的流动方向，针对各所述子体积计算属性P相对于平均属性值<P1>的标准偏差σi；
g)求出σi与σvol良好地匹配的所有候选代表性子体积；
h)选择最小的候选并且存储所述最小的候选作为代表性单元体积；以及
i)使用所述代表性单元体积来推导至少一个关注属性值。
14.根据权利要求13所述的方法，其中，还包括以下步骤：
针对所述分割体积整体，推导第二目标函数P2的平均属性值<P2>；
针对所述分割体积整体，计算相对于平均属性值<P2>的标准偏差σvol；
在所述体积内定义多个子体积；
针对各所述子体积，计算第二目标函数P2相对于平均属性值<P2>的标准偏差σi；以及
针对第一目标函数P1和第二目标函数P2的组合，求出σi与σvol良好地匹配的所有代表性子体积。
15.根据权利要求14所述的方法，其中，第一目标函数P1是孔隙率，并且第二目标函数P2是所述孔隙空间的表面积与体积之比。
16.根据权利要求15所述的方法，其中，还包括以下步骤：在选择之前鉴定候选子体积，包括确定所述候选子体积用来根据达西定律推导流体输送属性的适合性，该步骤包括：
构建目标函数的标准偏差的分布；
对目标函数的标准偏差的分布的平均值或者可选地任何其它一阶特性以及所述分布的方差、峰度或偏度进行评价；
对一阶以上的动量相对于所述子体积的维度的趋势进行评价；以及
在一阶动量相对于所述一阶动量针对在较大子体积上构建的分布的值的变化至少为0.1的情况下、以及/或者在针对方差而言高阶动量高于特定阈值0.1的情况下，停止使所述子体积的维度减小。
17.一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的方法，包括以下步骤：
1)将多孔介质的分割三维图像加载至计算机系统中，
其中，所述分割三维图像包括各自分配有灰度值的体素；
2)选择被定义为Z方向的流动方向；
3)定义探询体积的大小，其中：
i.所述探询体积是维度为Xi、Yi和Zi的原始分割三维图像的子样本，其中整个样本的维度为Xs,Ys,Zs，
ii.定义探询体积的最大值imax，
iii.设置各所述探询体积的以体素为单位的维度(Xi,Yi,Zi)，其中Xi、Yi和Zi是针对值为1～imax的i所定义的，以及
iv.将i的初始值设置为1；
4)针对所述探询体积的各切片来计算所选择属性Ps(0,0,0)～Ps(0,0,Zs)；
5)计算σs(0,0,0)，其中σs是一组所选择属性Ps(a,b,c)～Ps(a,b,c+Zi)的标准偏差，a、b、c表示当前探询体积的位置坐标；
6)设置大小为Xi,Yi,Zi的所述探询体积在大小为Xs,Ys,Zs的整个样本内所占据的最大坐标，其中：
i.amax＝Xs–Xi+1，
ii.bmax＝Ys–Yi+1，以及
iii.cmax＝Zs–Zi+1；
7)将当前探询体积的位置坐标设置为a＝b＝c＝0；
8)针对所述当前探询体积的切片计算所选择属性Pi(a,b,c)～Pi(a,b,c+Zi)，
i.其中，所述所选择属性包括孔隙率、表面积与体积之比、相似属性或它们的任意组合；
9)计算σi(a,b,c)，其中σi是相对于整个样本的一组所选择属性Pi(a,b,c)～Pi(a,b,c+Zi)的标准偏差，
i.其中，可选地，对计算σi的值所使用的Pi的值进行滤波，以及
ii.其中，可选地，设置Pi的平均值；
10)使所述探询体积的位置在X方向上移动1个体素、即a＝a+1；
11)重复步骤8)～10)并且存储Pi和σi的所有值，直到所述当前探询体积的X坐标的值a等于所述当前探询体积能够占据的最大值amax为止；
12)将所述当前探询体积的X坐标设置为零即a＝0，并且使当前位置体积的Y坐标增加1个体素、即b＝b+1；
13)重复步骤8)～12)并且存储Pi和σi的所有值，直到所述当前探询体积的Y坐标的值b等于所述当前探询体积能够占据的最大值bmax为止；
14)将所述当前探询体积的X坐标设置为零、即a＝0，将所述当前探询体积的Y坐标设置为零、即b＝0，并且使当前位置体积的Z坐标增加1个体素、即c＝c+1；
15)重复步骤8)～14)并且存储Pi和σi的所有值，直到所述当前探询体积的Z坐标的值c等于所述当前探询体积能够占据的最大值cmax为止；
16)增大所述当前探询体积的大小，其包括：
i.通过使指针增加至下一探询体积、即i＝i+1来选择下一组探询体积，以及
ii.将当前探询大小设置为Xi,Yi,Zi；
17)重复步骤6)～16)，直到选择了所有探询体积并且计算并存储了Pi和σi的所有值为止；
18)选择用以匹配的一个或多个所选择属性；
19)针对各探询体积计算λi，其中，μ是所有σj的平均值，λ的索引i用于特定目标函数；
20)选择λi为最小值的探询体积，其中所选择的探询体积是代表性单元体积即REV的大小和位置；以及
21)计算所述多孔介质的属性。
18.根据权利要求17所述的方法，其中，所述分割三维图像是如下产生的：通过利用计算机断层成像x射线扫描器扫描所述样本获得所述样本的图像，并且利用软件程序对所述图像进行分割，以将体素分类为颗粒、孔隙和可选地其它相。
19.根据权利要求17所述的方法，其中，所述属性包括常规岩心分析属性即RCAL属性、特殊岩心分析属性即SCAL属性或这两者。
20.根据权利要求19所述的方法，其中，所述RCAL属性是沿多个轴的孔隙率、油母质含量、绝对渗透率，并且所述SCAL属性是相对渗透率、毛细管压力、颗粒大小分布、电气属性、弹性属性和它们的任意组合。
21.一种用于识别与多孔介质的样本相对应的子样本代表性数字体积的系统，包括：
a)扫描器，其能够产生所述多孔介质的三维数字图像；
b)计算机，其包括至少一个处理器，其中所述至少一个处理器用于执行能够获得表现孔隙空间和至少一个固相的特性的分割体积的计算机程序，
c)与b)相同或不同的计算机，其包括至少一个处理器，其中所述至少一个处理器用于执行能够进行计算的计算机程序，其中所述计算包括：i)针对所述分割体积整体，推导第一目标函数P1的平均属性值<P1>；ii)针对所述分割体积整体，计算相对于平均属性值<P1>的标准偏差σvol；iii)在所述体积内定义多个子体积；iv)针对各所述子体积，计算第一目标函数P1的属性值P相对于平均属性值<P1>的标准偏差σi；v)求出标准偏差σi与σvol良好地匹配的所有候选代表性子体积；vi)从候选中选择并存储代表性子体积；以及vii)使用所述代表性子体积来推导至少一个关注属性值；以及
d)用以显示、打印或存储所述计算的结果的至少一个装置。
CN201380029060.XA 2012-03-30 2013-02-04 从多孔介质的数字表征中选择代表性单元体积的高效方法 Expired - Fee Related CN104335042B (zh)

Applications Claiming Priority (5)

Application Number Priority Date Filing Date Title
US201261618265P 2012-03-30 2012-03-30
US61/618,265 2012-03-30
US13/546,053 US20130262028A1 (en) 2012-03-30 2012-07-11 Efficient Method For Selecting Representative Elementary Volume In Digital Representations Of Porous Media
US13/546,053 2012-07-11
PCT/US2013/024593 WO2013147995A1 (en) 2012-03-30 2013-02-04 An efficient method for selecting representative elementary volume in digital representations of porous media

Publications (2)

Publication Number Publication Date
CN104335042A CN104335042A (zh) 2015-02-04
CN104335042B true CN104335042B (zh) 2016-08-24

Family

ID=49236161

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201380029060.XA Expired - Fee Related CN104335042B (zh) 2012-03-30 2013-02-04 从多孔介质的数字表征中选择代表性单元体积的高效方法

Country Status (10)

Country Link
US (1) US20130262028A1 (zh)
EP (1) EP2831578A1 (zh)
CN (1) CN104335042B (zh)
AU (1) AU2013240570B2 (zh)
BR (1) BR112014024357B1 (zh)
CA (1) CA2868872C (zh)
CO (1) CO7081149A2 (zh)
MX (1) MX336642B (zh)
RU (1) RU2586397C2 (zh)
WO (1) WO2013147995A1 (zh)

Families Citing this family (31)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9047513B2 (en) 2012-08-10 2015-06-02 Ingrain, Inc. Method for improving the accuracy of rock property values derived from digital images
US20140052420A1 (en) * 2012-08-20 2014-02-20 Ingrain Inc. Digital Rock Analysis Systems and Methods that Estimate a Maturity Level
US9070049B2 (en) * 2013-03-15 2015-06-30 Bp Corporation North America Inc. Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties
US20140270393A1 (en) * 2013-03-15 2014-09-18 Bp Corporation North America Inc. Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties
US20160305237A1 (en) * 2013-12-05 2016-10-20 Schlumberger Technology Corporation Method and system of showing heterogeneity of a porous sample
AU2015241029B2 (en) * 2014-03-31 2018-02-01 Ingrain, Inc. Representative elementary volume determination via clustering-based statistics
CN107111895A (zh) * 2015-01-30 2017-08-29 惠普发展公司有限责任合伙企业 可变密度建模
US10181216B2 (en) 2015-01-30 2019-01-15 Hewlett-Packard Development Company, L.P. Generating slice data from a voxel representation
US9892321B2 (en) * 2015-06-27 2018-02-13 New England Research, Inc. Using maximal inscribed spheres for image-based rock property estimation
US9841387B2 (en) * 2015-07-22 2017-12-12 Test Research, Inc. Inspection method and device
US10394974B2 (en) * 2015-10-01 2019-08-27 Arizona Board Of Regents On Behalf Of Arizona State University Geometry based method for simulating fluid flow through heterogeneous porous media
WO2017055647A1 (es) * 2015-10-02 2017-04-06 Repsol, S.A. Método para proporcionar un modelo numérico de una muestra de roca
CN105651964B (zh) * 2015-12-29 2017-11-03 河海大学 一种确定裂隙岩体表征单元体积的方法
US10431419B2 (en) 2016-07-19 2019-10-01 Battelle Memorial Institute Sparse sampling methods and probe systems for analytical instruments
CN107290358B (zh) * 2017-06-20 2019-11-08 大连理工大学 一种应用ct测量多孔介质内co2-盐水界面面积变化的方法
US10256072B2 (en) 2017-08-01 2019-04-09 Battelle Memorial Institute Optimized sub-sampling in an electron microscope
CN107680131B (zh) * 2017-09-08 2020-09-11 北京理工大学 一种快速确定多孔介质表征单元体积尺寸的方法
CN107784169B (zh) * 2017-10-11 2020-10-13 上海电力学院 基于变差函数及孔隙度的多孔介质表征单元体确定方法
US11609174B2 (en) * 2017-11-06 2023-03-21 Khalifa University of Science and Technology Method and system for determining permeability of a porous medium
CN110029989B (zh) * 2018-01-11 2021-11-02 中国石油化工股份有限公司 一种非常规油气采出程度计算方法及系统
CN108320307B (zh) * 2018-01-19 2020-09-08 中国石油天然气股份有限公司 一种确定储层岩石样品的有效单元体积的方法及装置
CN111742207B (zh) * 2018-04-23 2023-10-17 日本碍子株式会社 确定有效或无效流路的方法和装置
US10969323B2 (en) 2018-05-30 2021-04-06 Saudi Arabian Oil Company Systems and methods for special core analysis sample selection and assessment
US11009497B2 (en) 2018-06-22 2021-05-18 Bp Corporation North America Inc. Systems and methods for estimating mechanical properties of rocks using grain contact models
CN109191423B (zh) * 2018-07-18 2019-07-02 中国矿业大学 一种基于机器图像智能学习的多孔介质渗透率预测方法
CN109145407B (zh) * 2018-08-01 2020-08-04 浙江大学 一种基于切片截面的隐式曲面多孔实体结构性能分析方法
US11454111B2 (en) 2020-01-30 2022-09-27 Landmark Graphics Corporation Determination of representative elemental length based on subsurface formation data
CN111597763B (zh) * 2020-04-09 2023-06-27 东华理工大学 十字型多尺度流孔隙介质全频段弹性波频散衰减分析方法
US11703440B2 (en) 2021-09-24 2023-07-18 General Electric Company Porosity of a part
US12049818B2 (en) 2022-01-14 2024-07-30 Halliburton Ener y Services, Inc. Upscaling of formation petrophysical characteristics to a whole core scale
CN116704008B (zh) * 2023-08-01 2023-10-17 城云科技（中国）有限公司 基于画面面积计算判断物体的方法、装置及其应用

Family Cites Families (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4903207A (en) * 1986-05-15 1990-02-20 Restech, Inc. Method for determining reservoir bulk volume of hydrocarbons from reservoir porosity and distance to oil-water contact level
US5430291A (en) * 1992-05-01 1995-07-04 Texaco Inc. X-ray CT measurement of fracture widths and fracture porosity in reservoir core material
US6516080B1 (en) 2000-04-05 2003-02-04 The Board Of Trustees Of The Leland Stanford Junior University Numerical method of estimating physical properties of three-dimensional porous media
GB2414072B (en) * 2004-05-12 2006-07-26 Schlumberger Holdings Classification method for sedimentary rocks
US7516055B2 (en) * 2004-08-20 2009-04-07 Chevron U.S.A. Inc Multiple-point statistics (MPS) simulation with enhanced computational efficiency
US7630517B2 (en) * 2005-07-13 2009-12-08 Schlumberger Technology Corporation Computer-based generation and validation of training images for multipoint geostatistical analysis
US8510242B2 (en) * 2007-08-31 2013-08-13 Saudi Arabian Oil Company Artificial neural network models for determining relative permeability of hydrocarbon reservoirs
US8068579B1 (en) * 2008-04-09 2011-11-29 Xradia, Inc. Process for examining mineral samples with X-ray microscope and projection systems
CN101802649B (zh) * 2008-04-10 2013-01-23 普拉德研究及开发股份有限公司 利用井眼图像、数字岩石样品以及多点统计算法生成数值假岩心的方法
US8081796B2 (en) 2008-11-24 2011-12-20 Ingrain, Inc. Method for determining properties of fractured rock formations using computer tomograpic images thereof
US8155377B2 (en) 2008-11-24 2012-04-10 Ingrain, Inc. Method for determining rock physics relationships using computer tomographic images thereof
US8170799B2 (en) 2008-11-24 2012-05-01 Ingrain, Inc. Method for determining in-situ relationships between physical properties of a porous medium from a sample thereof
US8085974B2 (en) 2008-11-24 2011-12-27 Ingrain, Inc. Method for determining elastic-wave attenuation of rock formations using computer tomograpic images thereof
US8081802B2 (en) 2008-11-29 2011-12-20 Ingrain, Inc. Method for determining permeability of rock formation using computer tomograpic images thereof
US20110004447A1 (en) 2009-07-01 2011-01-06 Schlumberger Technology Corporation Method to build 3D digital models of porous media using transmitted laser scanning confocal mircoscopy and multi-point statistics
US8805616B2 (en) * 2010-12-21 2014-08-12 Schlumberger Technology Corporation Method to characterize underground formation

Also Published As

Publication number Publication date
EP2831578A1 (en) 2015-02-04
WO2013147995A1 (en) 2013-10-03
US20130262028A1 (en) 2013-10-03
MX336642B (es) 2016-01-26
MX2014011709A (es) 2014-12-08
RU2014143802A (ru) 2016-05-27
AU2013240570B2 (en) 2015-05-14
BR112014024357B1 (pt) 2021-06-22
CA2868872A1 (en) 2013-10-03
CO7081149A2 (es) 2014-10-10
RU2586397C2 (ru) 2016-06-10
CA2868872C (en) 2017-05-16
AU2013240570A1 (en) 2014-10-16
CN104335042A (zh) 2015-02-04

Similar Documents

Publication Publication Date Title
CN104335042B (zh) 从多孔介质的数字表征中选择代表性单元体积的高效方法
Silin et al. Robust determination of the pore space morphology in sedimentary rocks
Vogel et al. Quantification of soil structure based on Minkowski functions
US9046509B2 (en) Method and system for estimating rock properties from rock samples using digital rock physics imaging
WO2017041281A1 (en) Porous media anaylysis system and method
Talukdar et al. Stochastic reconstruction of chalk from 2D images
US20210182597A1 (en) Process parameter prediction using multivariant structural regression
Thiele et al. Insights into the mechanics of en-échelon sigmoidal vein formation using ultra-high resolution photogrammetry and computed tomography
Alqahtani et al. Super-resolved segmentation of X-ray images of carbonate rocks using deep learning
WO2014062947A2 (en) Method for modeling a reservoir using 3d multiple-point simulations with 2d training images
US20140052420A1 (en) Digital Rock Analysis Systems and Methods that Estimate a Maturity Level
Wu et al. Two-phase flow in heterogeneous porous media: A multiscale digital model approach
Hajizadeh et al. An algorithm for 3D pore space reconstruction from a 2D image using sequential simulation and gradual deformation with the probability perturbation sampler
Chen et al. Reconstruction of multiphase microstructure based on statistical descriptors
Chen et al. Intelligent adaptive sampling guided by Gaussian process inference
Li et al. Digital Rock Reconstruction Using Wasserstein GANs with Gradient Penalty
Qajar et al. Chemically induced evolution of morphological and connectivity characteristics of pore space of complex carbonate rock via digital core analysis
Yang et al. A new voxel upscaling method based on digital rock
Biswal et al. Towards precise prediction of transport properties from synthetic computer tomography of reconstructed porous media
Han et al. A random algorithm for 3D modeling of solid particles considering elongation, flatness, sphericity, and convexity
Franc et al. Image-based effective medium approximation for fast permeability evaluation of porous media core samples
Liu et al. Reconstruction of 3D multi-mineral shale digital rock from a 2D image based on multi-point statistics
Hu et al. Correlating recovery efficiency to pore throat characteristics using digital rock analysis
Silva et al. Well-testing based turbidite lobes modeling using the ensemble smoother with multiple data assimilation
Bashtani Random Network Modeling of Tight Formations

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20160824

Termination date: 20190204