从核磁共振数据获得质子密度分布的制作方法
专利名称：从核磁共振数据获得质子密度分布的制作方法
技术领域：
本发明涉及从核磁共振数据获得质子密度分布。
背景技术：
用来操控存在于多孔介质中的分子的自旋的核磁共振系统是熟知的。这些系统通 常至少执行以下的功能来确定有关多孔介质和/或其中所包含的流体的信息通过静磁场 使得自旋极化、通过射频(“RF”)脉冲操控自旋、和通过从多孔介质发出的RF信号接收自 旋的响应。然后处理RF信号，以确定有关多孔介质的成分和/或在多孔介质内所包含或 “包围”的一种或多种流体的信息。
发明内容
本发明一方面涉及一种获得质子密度分布的计算机实施的方法。在一个实施例 中，该方法包括获取来自多孔介质的核磁共振数据；经由全局优化算法反演核磁共振数据， 以确定在多孔介质内的质子密度分布；以及输出所确定的质子密度分布。本发明的另一方面涉及一种获得质子密度分布的计算机实施的方法。在一个实施 例中，该方法包括获取来自多孔介质的核磁共振数据；从核磁共振数据确定多孔介质的质 子密度分布，其中质子密度分布包括由多个预定区域组成的一个或多个频谱，一个或多个 预定区域具有在其中的分布峰值，以及其中确定质子密度分布包括通过使得单个基本分量 适配于在其中具有分布峰值的每个预定区域中的分布而按照非线性基函数对质子密度分 布进行参数化；以及输出所确定的质子密度分布。本发明的另一个方面涉及一种获得有关质子密度分布的信息的方法。在一个实施 例中，该方法包括获取来自介质的核磁共振数据；定义一个实现所获取的核磁共振数据、并 依赖于质子密度分布的m个参数的函数，以使得当接近对于质子密度分布的m个参数的解 时函数被最小化；实施全局优化算法，以确定对于质子密度分布的m个参数的解；以及输出 对于质子密度分布的m个参数的解。本发明的这些和其它目的、特征和特性，以及操作步骤的方法和结构的相关单元 的功能和部件的组合和制造经济性，在参照附图考虑以下的说明和所附权利要求后，将变 得更加清楚，所有的这些形成本技术说明书的一部分，其中相同的附图标记指明在各个图 上的相应的部分。然而，应当明显看到，附图仅仅是用来显示和描述，而不是用来规定本发 明的范围。正如在技术说明书和权利要求中使用的那样，单数形式“一”、“一个”和“该”均 可包括多个单元，除非上下文明确地另外规定。
图1示出了按照本发明的一个实施例的、实施核磁共振以确定有关多孔介质的信 息的方法。图2提供了按照本发明的一个实施例的核磁共振数据的示例图。
图3提供了按照本发明的一个实施例的质子密度分布的示例图。图4示出了具有两个局部最小值和一个绝对最小值的三维表面。图5示出了具有与噪声有关的另一个局部最小值的图4所示表面。图6示出了按照本发明的一个实施例的实施全局优化算法的方法的一个实施例。图7示出了按照本发明的一个实施例的经由全局优化算法反演核磁共振数据，以 确定质子密度的方法。图8示出了质子密度分布。
具体实施例方式核磁共振(“NMR”)技术已广泛应用于测量包含流体的多孔介质的特性(例如地 质构造的岩石物理学特性、活体组织的各种特性、化学化合物的结构特性等等)。某些岩石 物理学特性的例子包括孔隙尺寸、表面积与体积的比率、构造渗透率和毛细管压力。在确定 这些特性时，诸如纵向弛豫时间(“T/’)、横向弛豫时间(“V，)和扩散系数(“D”)那样 的参数常常是感兴趣的。弛豫时间是与核自旋在激励后恢复到它们的平衡位置相关联的时 间。纵向弛豫时间是同自旋与外部静磁场的对准有关的。横向弛豫时间T2是标识在取 向为与主磁场成一个角度的自旋之间发生的相位相干性丢失的时间常数。这种丢失部分是 由于自旋之间的相互作用造成的。扩散系数D是孔隙流体的扩散系数。图1示出了实施NMR以确定有关多孔介质和/或包含在其中的流体的信息的方法 10。应当意识到，虽然方法10的某些方面在这里是在实施NMR以确定有关地质介质的信息 的上下文内描述的，但这并不意图限制本公开内容的范围。本公开内容的范围包括实施NMR 以检验除地质介质以外的介质(例如，活体组织等等)。在一个实施例中，方法10是在一个 或多个处理器上执行的计算机实施方法。在操作步骤12，多孔隙介质和包含在其中的流体由被称为“脉冲序列”的RF脉冲 串扰动。在一个实施例中，该扰动是按照Carr-Purcell-Meiboom-GilKCPMG)脉冲序列完 成的。在某些实施例中，施加其它序列来扰动多孔介质和包含在其中的流体。常规的CPMG脉冲序列包括90度脉冲，后面跟随一系列180度脉冲，其遵循以下的 表达式⑴图其中Ji/2是相对于旋转帧中的基准沿+和-x轴施加的90度脉冲；^代表第k 个回波串的回波间隔TEk的一半；Ji代表在旋转帧中沿y轴施加的180度脉冲；acq代表回 波的获取；nk是在第k个回波串中的回波数。 伴随常规的CPMG脉冲序列一起的是外部磁场梯度Gk和在180度RF脉冲之间施加 的可选的脉冲梯度。脉冲场梯度的宽度用Sk表示，以及在连续的梯度脉冲之间的间隔用 Ak表示。同时，核自旋的系统也受到由外部磁场造成的内部场梯度以及在细颗粒与孔隙流 体之间的敏感度反差作用。在下面，符号g代表核自旋的系统受到的总的磁场梯度。
在操作步骤14，获得来自多孔介质内的RF信号。这些信号通常被称为NMR数据、 “回波数据”或NMR记录数据。正如下面进一步讨论的，这些信号表示响应于操作步骤12的 脉冲序列在多孔介质中产生的一个或多个回波串。关于由在操作步骤14获得的信号传递
的回波串的信息使得能够测量多孔介质和/或包含在其中的流体的一个或多个特性(例如 以上枚举的特性)。在某些实施例中，操作步骤12和14是通过使用已知的NMR记录工具和 /或NMR频谱仪而被执行的。 在操作步骤14获得的对应于脉冲序列的一个回波串的NMR数据M(ti)具有使得 数据M(ti)与参数弛豫时间、T2弛豫时间和扩散系数D相关联的以下的形式
「00241 其中、代表回波之间的时间，、代表回转磁比率，WT代表在回波串之间的等待时 间，a (PARAMETERS)代表给定的参数组(例如，PARAMETERS = T” T2和/或D)的幅值，以 及^代表在NMR数据M(t》中存在的噪声。图2提供在操作步骤14得到的实际的数据 M(t》作为时间的函数的示例图。回到图1，在某些情形下，两个或更多个参数1\，1~2和/或D可以从公式(2)被消 去。例如，当磁场被认为是均勻时，磁场的梯度g等于0，以及被称为扩散核的第三个指数项 广、^^2变为接近于1，这样，公式变为与D无关的。同样地，当等待时间WT足够大时，被 称为极化因子的第二个指数项(丨-广“777”变为接近于1，这样，公式变为与无关。或者， 通过使用有效比率r = I/T2可以使得极化因子变为与无关。在这种情况下，k-和值也 消失。在公式⑵变为与1\和0无关的情况下，在M(ti)与1~2之间的关系可被表示为 其中M(ti)代表所测量的第i个回波幅值(i = 1，. . .，!!)，、代表一组等间隔的n 个衰变时间，以及^是第i个回波的噪声。T”代表一组在对数标尺上等间隔的m个预先 选择的T2弛豫时间，以及a (T2J)是本模型要解决的与弛豫时间T”相关联的T2幅值分布。此后，本公开内容主要描述在其中以NMR数据表示的回波串在计算上与和D无 关(即，M(tD主要取决于T2)的实例中在操作步骤14中获得的NMR数据的处理。这些实 例通常被称为单维NMR解。然而，这并不意图作为限制。应当意识到，这里对于单维NMR解 描述的技术和方法可以很容易地扩展到多于一维的NMR解(例如，在M(t》取决于和/ 或D，以及1~2的实例中)。在操作步骤16，从在操作步骤14得到的NMR数据获得质子密度分布。所获得的质 子密度分布是M(ti)所取决(例如，在单维NMR解中的T2)的参数(例如1\，1~2和/或D)的 函数。例如，图3提供从在图2中表示的NMR数据获得的质子密度分布a (T2J)的示例图。
回过来参照图1，在一个实施例中，从NMR数据获得质子密度分布包括执行反演以 根据在公式(3)中表示的关系来确定分布。例如，对于不带有噪声项、由公式(3)表示 的、对于分布a (T2J)的形式上的数学解通过拉普拉斯逆变换给出。然而，噪声项^使得 这种技术成为无法实施的，因为该反演是非唯一的。另一个方法是把公式(3)重新表达为 最小平方最小化
它是对于使x 2最小的最佳离散分布a (T2J)进行求解。在操作步骤17，输出在操作步骤16得到的质子密度分布。这可包括所述分布的电 子存储、所述分布的电子显示(例如，图形显示)、所述分布的电子传输、或者用于输出通过 计算机实施方法导出的信息的任何其它方法。通常，在操作步骤16被用来确定对于这个最小平方最小化的解的方法包括实施 奇异值分解(“SVD”)。为了给出这种常规技术的概略说明，我们把公式(3)和(4)重写 为
其中 h = M (、), aj = a (T2J)，并且 以矩阵形式表示，公式(5)和(6)变为(7)b = K a+ ￡ ；并且(8) x 2(a) = | K a~b|2 ；其中矩阵元素对应于公式(5)和(6)的分量。为了求解方程(8)，加上条件=0，
这是为了使x2最小而必须满足的条件。这得到(9)Kt (K a-b) = 0 ；其被简化为(10) (KT K) a = KT b ；其被进一步简化为(11) K-a = b;其中上标“T”代表转置。公式(11)是线性方程组，它可以通过使用代数方法，诸 如Gauss-Jordan消元法而进行求解。然而，由于噪声项￡ i，矩阵K是奇异的或接近于奇异。 这主要是因为矩阵K是超定的(即，K是nXm矩阵，其中n>m)，它造成在由噪声给出的不 确定性内的模糊的可能的解a，产生实际上的不确定系统。矩阵K的SVD形式是(12) K = U [diag (Wi) ] VT其中U是[mXn]列正交矩阵，V是[nXn]正交矩阵，以及W是具有正的或接近于 零的元素(奇异值)的[nXn]对角矩阵。与数据的全局噪声电平有关的截止条件被用来 使满足( wfflaxX 6 cutoff的所有对角值为零，其典型值5 cutoff = 10_5。通过使用公式(11) 中K的SVD形式，对于a的解可被表达为 其中向量U⑴，i = 1，...，m表示U的列(每个是长度为n的向量)，以及向量V⑴， i = 1，...，m表示V的列(每个是长度为m的向量)。公式(13)规定适配的参数a是V的 列的线性组合，具有通过形成U的列与加权的数据向量b的点积而得到的系数。尽管截止条件是指平滑掉所有的噪声效应，但实际上它不足以产生平滑的幅值分 布，其中可能出现某些尖峰。在公式(13)上加上被称为“规则化项”的附加项，它用来通过要由搜索方法确定的、匹配于在从操作步骤14得到的数据中的数据噪声电平的某些因子 使得所述分布的模、斜率或曲率最小化。在一个实施例中，作为对于SVD和/或用于对所得到的NMR数据进行反演的其它 常规技术的替代，实施全局优化算法。更具体地，函数x2(a)被看作为m维表面，并且全局 优化算法被实施以确定该表面的绝对最小值。所确定的该表面的绝对最小值的m空间向量 坐标对应于关于M(ti)所取决的参数(例如，T2)的质子密度分布。确定x2(a)的m维表面的绝对最小值并非是不重要的，这部分是因为该表 面很可能包括多个局部最小值。为此，使用局部梯度信息来检测表面最小值-诸如 Levenberg-Marquardt或共轭梯度-的常规优化算法是不现实的(因为它们通常产生局部 最小值而不是绝对最小值)。例如，图4示出了具有局部最小值和绝对最小值的三维表面。 从图4应当意识到，不提供全局优化的优化算法可能变为在一个局部最小值上“固定位”。在一个实施例中，为了减小回到局部最小值的机会，在操作步骤16实施的优化算 法是全局优化算法。正如这里使用的，术语“全局优化算法”是指可被执行以使得函数或函 数组优化到某个预定标准的任何算法或算法组。通常，全局优化算法需要一系列迭代，使得 算法能够收敛到最佳。例如，全局优化算法可包括随机优化(例如，模拟退火、并行回火等 等)、探试法或亚探试法优化(例如，遗传算法、渐近策略等等)和/或其它优化中的一项或 多项。通过实施全局优化算法从NMR数据得到质子密度分布而实现的反演仅仅是x2(a) 函数的绝对最小值的估计。此外，在操作步骤14得到的数据中存在的噪声还可能使得确 定x2(a)最小值进一步复杂化。例如，图5示出了带有与噪声相关联的其它局部最小值的 图4所显示的表面。正如应当从图5意识到的，当接近在不存在噪声时的x2(a)的最小值 (例如图4中的x2(a)的最小值)时，由噪声引起的、包围这个最小值的局部最小值造成对 于x2(a)最小值的等价的(或基本上等价的)解的多重性。这样，经由全局优化算法达到 的单个解不能提供x2(a)函数的绝对最小值，因为全局优化算法会返还与噪声相关联的一 个最小值，而不是实际的最小值。图6示出了实施全局优化算法来考虑(至少在某种程度 上)由噪声引起的等价解的多重性的方法18。在一个实施例中，方法18保证所得到的结果 将提供足够的精度。在某些实例中，方法18可以在方法10的操作步骤16(在图1中显示 并在以上描述)时被实施。然而，这并不意图限于此。方法18包括操作步骤20，在该操作步骤，通过使用全局优化算法得到对于x 2 (a) 的绝对最小值的解。在操作步骤22，更新得到对于x2(a)的解的次数的计数值(例如，加 上1，以反映操作步骤20的解)。在操作步骤24，判断对于x2(a)的绝对最小值的解进行 附加估值是否是必要的。在一个实施例中，这种判断包括将在操作步骤22更新的计数值与 一个预定的阈值进行比较。如果需要附加估值，则方法18回到操作步骤20。如果不需要附 加估值，则方法18进行到操作步骤26。在操作步骤26，对在操作步骤20确定的解的估值 进行合计，以提供对于x2(a)函数的绝对最小值的解的最终估值。例如，所述合计可包括 对所述估值进行平均。通过对解的多个估值进行合计，解的最终估值的精度考虑了(至少 在某种程度上)由噪声引起的等价解的多重性。应当意识到，在方法18中对于x2(a)的 最小值的多个独立的估值的合计在估值被线性参数化的实例中可以有更大的好处，见公式 (4)，因为x2(a)的最小化的非线性参数化会减小噪声对解的影响，见公式18。然而，这并
8不意图限于此，并且方法18可以结合解的非线性参数化一起使用，而不背离本公开内容的 范围。图7示出了按照本发明的一个实施例的经由全局优化算法反演NMR数据以确定质 子密度的方法28 (例如，在方法10的操作步骤16，在方法18的操作步骤20等等)。具体 地，方法28包括实施模拟退火算法。应当意识到，方法28中的模拟退火算法的说明不应当 看作为限制。正如以上阐述的，可以预期，在根据NMR数据确定质子密度时可以实施各种各 样的全局优化算法。通常，模拟退火算法是基于对于物理系统的模拟组成m空间向量a的分量的一 组变量是形成配置空间的假想物理系统的自由度。函数E(a) = x2(a)被认为是系统的能 量，问题被归结为找出系统的最小能量或基础状态配置的问题。众所周知，如果假想物理系 统被加热到非常高的温度T，然后慢慢地冷却到绝对零度(称为退火的过程)，则系统将发 现自己处在基础状态(即，在配置空间中具有最小量的维持能量的配置)。冷却速率必须 足够慢以避免陷入某个亚稳定状态(即，E(a)的局部最小值)。在温度T时，处在具有能量 E(a)的配置中的概率由Boltzmann-Gibbs因子给出(14)exp(_E(a)/T)。Boltzmann-Gibbs因子表明高的能量配置在高的温度T下可以以有限的概率出 现。如果温度降低，则那些高的能量状态变为不太可能的，并且当T — 0时，仅仅接近于E (a) 的最小值的状态才呈现非空概率。因此，在由对于x2(a)的配置空间形成的假想物理系统 内，当a的配置被随机地(或伪随机地)操控并且温度被适当地降低以使得T — 0时，降 低的温度约束a的允许的配置，以使得a的配置慢慢地收敛到在配置空间中最低能量状态 (其中x2(a)被最小化)时的a的配置。在一个实施例中，方法28包括操作步骤30，在该操作步骤，设置当前温度T的初始 值。当前温度T的初始值被设置为足够高以使得方法28能够探查x2(a)函数的相当大量 的配置空间。在操作步骤32，获得对于a的参数(即，分量)的当前的数值组。例如，如果a由 m个参数组成，则在操作步骤32获得对于参数的m个数值。在一个实施例中，所述m个数值 被随机地确定。正如在本公开内容内使用的，术语“随机”包括纯随机和伪随机确定。在操作步骤34，确定对于a的参数的当前的数值组的x 2 (a)函数的数值。在模拟 退火的字典编辑中，这个数值是通过对于a的参数的m个数值的当前的组所描述的配置的 “能量” E (a)。在操作步骤36，获得对于a的m个参数的建议的数值组。在一个实施例中，操作 步骤36包括实施Monte Carlo方法来确定建议的数值组。这并不意图限于此，因为在操作 步骤36确定建议的数值组时可以实施其它算法，诸如混合Monte Carlo方法、Hamiltonian 动态法和/或其它算法。更具体地，在一个实施例中，建议的数值组是通过按照所建议的改 变来改变对于a的m个参数的当前的数值组而被确定的。由于a的参数的数值总是正的， 可以实施适当的随机动态方法来获得建议的数值组。例如，可以实施类似于以下公式的对 数正态随机挪动(15) log(a' ) = log(a) + o G ；其中a’是具有对于a的参数的建议的数值组的m空间向量，G是其所有的分量除
9了一个分量之外均等于零的m空间向量(类似于a)，这个不等于零的分量是按照正态分布 N(0，1)(中心为零且标准差为1的高斯分布)而分布的随机数。应当意识到，尽管这个实施 方案实际上仅将单个数值从当前的数值组变化到建议的数值组，但也可以使用其它的实施 方案来从当前的数值组确定建议的数值组，其中改变一个以上的数值。在操作步骤38，确定对于a的参数的建议的数值组的x 2的数值(即，x2(a’))。 在模拟退火的字典编辑内，在操作步骤38中确定的数值是通过对于a的参数的建议的数值 组所描述的配置的“能量”，它可以被表达为E(a’ )。在操作步骤40，确定接受对于a的m个参数的建议的数值组作为对于a的m个参 数的当前数值组的概率。这个建议的数值组以概率Pa。。(a’ |a)被接受。在一个实施例中， 概率Pa。。(a，a)按照Metropolis接受算法被如下地确定 其中T再次代表当前的温度。在其它的实施例中，可以在操作步骤40实施不同于 Metropolis接受算法的接受算法来确定接受建议的数值组的概率。在操作步骤42，按照在操作步骤40确定的概率，接受或拒绝对于a的m个参数的 建议的数值组。如果接受所建议的数值组，则建议的数值组变为对于a的m个参数的当前 数值组。如果拒绝所建议的数值组，则对于a的m个参数的当前数值组保持不变。在操作步骤44，确定对于当前的温度的新的数值。对于当前的温度的新的数值可 以按照退火计划表被确定，该退火计划表描述当前的温度为已经执行的操作步骤36-42的 迭代次数的函数。这样的退火计划表的例子可被表示为如下(17)T(k) = T(0)e_Ak ；其中k是已经执行的操作步骤36-42的迭代次数。应当意识到，在某些实例中，当 前的温度可能不能对于每次迭代都在操作步骤44被更新。例如，为了节省处理资源和/或 增强算法的速度，在图7中操作步骤44的说明也给出其中当前的温度仅仅在每几次迭代才 被更新的实例。在操作步骤46，判断该算法是否应当执行另一次迭代。在一个实施例中，这种判断 是基于当前的温度是否允许能量的很大的减小。如果确定不应当进行另一次迭代，则对于 a的m个参数的当前的数值组被确定为代表质子密度分布。图8示出了另一种质子密度分布。在以上对方法28的说明中，所确定的m个参数 是线性的(即，幅值a (T2)为1~2的函数)。然而，这并不意图限于此。通过实施全局优化 算法以确定质子密度分布而使得能够进行的增强方案之一是可以按照非线性基函数对所 述解进行参数化。实施非线性基函数来对全局优化算法的解进行参数化(例如，图7中所 示的模拟退火技术)可以减小对于描述所述解所需要的参数的数目。这又意味着用非线性 基函数来对解进行参数化将减小对于达到该解所需要的计算量。参数的数目的减小在多维 解中甚至变为更显著的，因为必须被确定的参数的数目(和/或通过非线性参数化而变为 可能的必须被确定的参数的数目的减小)在添加每个附加维度时被指数地加倍。其它增强 方案也伴随有非线性参数化，其中的某些将在下面进一步讨论。此后，将讨论按照高斯基函 数对解进行非线性参数化，作为非线性参数化的一个例子。这并不意图限于此，因为也可以实施各种其它非线性基函数。例如，非线性参数化可以按照Gamma基函数、B-仿样基函数 或实验确定的基函数。如在图8中可以看到的，质子密度分布通常类似于具有一组连续峰值分布的频 谱。这些峰值分布通常对应于被包围在从其获得NMR数据的多孔介质中的不同类型的流 体。例如，在图8中所示的单维分布中的第一峰值45通常与泥土包围的流体相关联，第二 峰值50通常与毛细管包围的流体相关联，第三峰值52通常与“自由的”流体相关联。因为 峰值48、50和52的通常的形状，每个峰值48、50和52可以由单个高斯分布描述，其参数是 中心、宽度和幅值。在一个实施例中，按照非线性基函数对解进行参数化基本上包括与以上对于通过 执行方法28进行线性参数的解所描述的相同的操作步骤组。回过来参照图7，在操作步骤 32，获得对于参数向量的参数的数值组(例如，按照随机或伪随机确定法)。按照一组队 个高斯基函数分量对解进行参数化使得1~2优化问题能够被表示为(与以上的公式(4)相 比) 其中ak代表基本分量的幅值，bk代表中心，ck代表宽度，以及k代表不同系列的从 1到队的基本分量。从这个表示式可以看到，进行优化的参数可被组合成表示为x= {a” b^Ci，...，aNg，bNg，cNg}的单个m空间向量(其中m = 3 * Ng)。在操作步骤34，从公式(18) 求解对于x的m个参数的当前的数值组的x2(x)。在操作步骤36，确定对于x的m个参数的建议的数值组。由于x的每个分量(“x/’) 只代表几种可能的参数类型中的一种类型，比起以上对于线性参数的操作步骤36的说明， 这个操作步骤可能更多地被涉及到。具体地，由于x的每个分量可以是特定的参数类型(例 如，a, b，c)，对于每个参数类型可以建立可允许的范围，以便把优化以保持在系统的实际上 可能的限制范围内。例如，这可包括对于每种参数类型的最小值、最大值和范围的技术指标 (例如，adel = amax-amin, bdel = bmax-bmin, cdel = cmax-cmin)。从这些规定的约束条件，可以确定 一组m空间约束的向量xmin, x-，和xde1。在一个实施例中，操作步骤36再次包括实施受到约束的Lorentzian随机挪动形 式的Monte Carlo方法，以确定建议的数值组(下面用x’代表)。这个随机挪动可被表示 为(19) x' = x+o L其中o是缩放因子，它可以作为已经完成的操作步骤36-42的迭代次数的函 数而变化，以及L是m空间向量，其中所有的分量除了一个k以外均等于零，k是按照 Lorentzian分布而分布的随机数。必须检查新建议的数值处在对于该参数的允许的范围，否则生成新的建议。在 Lorentzian分布中的尾部大于高斯分布的尾部，所以达到有效的长范围的采样，这可以促 进模拟退火方法28的更佳性能，诸如通常被称为“快速模拟退火”的改进方案。在操作步骤38，从在操作步骤36得到的建议的数值组确定x2(x’)。在操作步骤40，确定接受建议的数值组的概率。例如，可以实施以上对于a和a’描述的Metropolis接 受概率(例如
在操作步骤42，按照在操作步骤40确定的
概率接受或拒绝建议的数值组。操作步骤44和46基本上如以上讨论的那样进行。回到图8，典型地，与泥土包围的流体相关联的第一峰值48位于由质子密度分布 形成的频谱的预定区域(例如，T2 e
)内。同样地，第二峰值50通常位于预定区域 (例如，T2 G [3,33])内，以及第三峰值52通常位于预定区域(例如，T2 G [33,300])内。 另外，有时存在第四峰值(图8中未示出)，它与第二种类型的自由流体相关联。第四峰值 通常也位于在图8中所示的频谱的预定区域(例如，T2e [300，3000]，以ys计)内。应 当意识到，其它区域划分也是可能的，并且这些数值仅仅被提供来用于说明的目的。在一个实施例中，以上讨论的优化方法(其中按照非线性参数组对所述解进行参 数化)，所述解被参数化，以使得单个基本分量(在优化期间)适配于在每个预定区域内的 质子密度分布。这进一步加强了优化，这样，在多孔介质内每种“类型”的流体的影响与单 个基本分量相关联，这使得即使在它穿过典型地与另一种流体类型相关联的频谱的区域的 情形下也能够跟踪对于给定的流体类型的质子密度分布的贡献。例如，这被显示于图8的 区域54。正如可以看到的，与第二峰值50和第三峰值52相关联的各个基本分量重叠。在 更惯用的解中(例如，实施线性参数的解)，在包含第二峰值50和第三峰值52的区域之间 的边界被用来分开对于与第二峰值50和第三峰值52相关联的流体类型的质子密度分布的 贡献。然而，通过把单个非线性基本分量适配于第二峰值50和第三峰值52的每个峰值而 能够进行的这些贡献的离散化使能穿过边界跟踪形成第二峰值50和第三峰值52的流体类 型的贡献。虽然本发明是根据当前认为的最实际和优选的实施例，为了说明的目的而详细地 描述的，但应当看到，这样的细节仅仅是为了该目的，本发明不限于所公开的实施例，而是 相反，本发明意图覆盖属于所附权利要求的主旨和范围内的修改方案和等价的设计。例如， 应当理解，本发明预期，在可能的范围内，任何实施例的一个或多个特性可以与任何另外的 实施例的一个或多个特性相组合。
权利要求
一种获得质子密度分布的计算机实施的方法，该方法包括获取来自多孔介质的核磁共振数据；经由全局优化算法反演核磁共振数据，以确定在多孔介质内的质子密度分布；以及输出所确定的质子密度分布。
2.权利要求1的方法，其中按照非线性基函数对质子密度分布进行参数化。
3.权利要求2的方法，其中非线性基函数包括高斯基函数、Gamma基函数、B-仿样基函 数或实验确定的基函数中的一个或多个基函数。
4.权利要求2的方法，其中质子密度分布包括由多个预定区域组成的一个或多个频 谱，一个或多个预定区域具有在其中的分布峰值，以及其中经由全局优化算法反演核磁共 振数据以确定质子密度分布包括使得单个基本分量适配于在其中具有分布峰值的每个预 定区域中的分布。
5.权利要求4的方法，其中各个预定区域被确定为对应于在多孔介质内的不同的流体 类型，以使得给定预定区域对应于多孔介质内的相应的流体类型。
6.权利要求1的方法，其中所述全局优化算法包括模拟退火、遗传算法、渐近策略或并 行回火中的一项或多项。
7.权利要求6的方法，其中通过模拟退火算法实施的采样技术包括MonteCarlo方法、 混合Monte Carlo方法、Hamiltonian动态法或随机挪动方法中的一项或多项。
8.权利要求1的方法，其中质子密度分布是η维分布，其中η大于2。
9.权利要求1的方法，其中经由全局优化算法反演核磁共振数据以确定质子密度分布 包括实施全局优化算法两次或多次，以确定对于质子密度分布的多个解；以及对于质子密度分布的多个解求平均。
10.权利要求9的方法，其中质子密度分布是幅值的线性表示。
11.一种获得质子密度分布的计算机实施的方法，该方法包括获取来自多孔介质的核磁共振数据；由核磁共振数据确定多孔介质的质子密度分布，其中质子密度分布包括由多个预定区 域组成的一个或多个频谱，一个或多个预定区域具有在其中的分布峰值，以及其中确定质 子密度分布包括通过使得单个基本分量适配于在其中具有分布峰值的每个预定区域中的 分布而按照非线性基函数对质子密度分布进行参数化；以及输出所确定的质子密度分布。
12.权利要求11的方法，其中非线性基函数包括高斯基函数、Gamma基函数、B-仿样基 函数或实验确定的基函数中的一个或多个基函数。
13.权利要求11的方法，其中各个预定区域被确定为对应于在多孔介质内的不同的流 体类型，以使得给定的预定区域对应于多孔介质内的相应的流体类型。
14.权利要求13的方法，其中所述预定区域包括对应于泥土包围的流体的预定区域、 对应于毛细管包围的流体的预定区域或者对应于一种或多种类型的自由的流体的区域的 一个或多个区域。
15.权利要求11的方法，其中由核磁共振数据确定介质的质子密度分布包括实施全局 优化算法。
16.权利要求15的方法，其中所述全局优化算法包括模拟退火、遗传算法、渐近策略或 并行回火中的一项或多项。
17.权利要求16的方法，其中通过模拟退火算法实施的采样技术包括MonteCarlo方 法、混合Monte Carlo方法、Hamiltonian动态法或随机挪动方法中的一项或多项。
18.权利要求11的方法，其中质子密度分布是n维分布，其中n大于2。
19.一种获得有关质子密度分布的信息的方法，该方法包括获取来自介质的核磁共振数据；定义一个实现所获取的核磁共振数据并与质子密度分布的m个参数有关的函数，以使 得当接近对于质子密度分布的m个参数的解时该函数被最小化；实施全局优化算法，以确定对于质子密度分布的m个参数的解；以及输出对于质子密度分布的m个参数的解。
20.权利要求19的方法，其中实施全局优化算法以确定对于质子密度分布的m个参数 的解包括(a)设置用于算法的当前温度的初始值；(b)随机地确定对于m个参数的当前数值组；(c)确定对于m个参数的当前数值组的函数值；(d)通过随机地调节在对于m个参数的当前数值组中的一个或多个数值而确定对于m 个参数的建议的数值组；(e)确定对于m个参数的建议的数值组的函数值；(f)根据对于m个参数的当前数值组的函数值、对于建议的参数组的函数值和当前的 温度，计算接受对于m个参数的建议的数值组作为对于m个参数的当前数值组的概率；(g)根据在(f)步骤计算的概率，接受或拒绝对于m个参数的建议的数值组作为对于m 个参数的当前数值组；以及(h)确定用于该算法的对于当前的温度的新的数值，使得当前的温度按该算法的迭代 的函数降低。
21.权利要求20的方法，其中(a)-(h)中的每个步骤以所述的次序执行，并且该方法还 包括在步骤(h)之后返回到步骤(d)。
22.权利要求19的方法，其中m个参数是非线性的。
23.权利要求19的方法，其中实施全局优化算法以确定对于m个参数的解包括实施全局优化算法两次或多次，以确定对于m个参数的多个解；以及对于m个参数的多个解求平均。
24.权利要求23的方法，其中m个参数是线性的。
25.权利要求19的方法，其中质子密度分布是n维分布，其中n大于2。
全文摘要
一种计算机实施的方法，使得能够获得质子密度分布。在一个实施例中，该方法包括获取来自多孔介质的核磁共振数据；经由全局优化算法反演核磁共振数据，以确定在多孔介质内的质子密度分布；以及输出所确定的质子密度分布。
文档编号G01V3/32GK101896834SQ200880120342
公开日2010年11月24日 申请日期2008年12月4日 优先权日2007年12月12日
发明者R·萨拉扎尔-蒂奥, 孙伯勤 申请人:雪佛龙美国公司
技术领域：
本发明涉及从核磁共振数据获得质子密度分布。
背景技术：
用来操控存在于多孔介质中的分子的自旋的核磁共振系统是熟知的。这些系统通 常至少执行以下的功能来确定有关多孔介质和/或其中所包含的流体的信息通过静磁场 使得自旋极化、通过射频(“RF”)脉冲操控自旋、和通过从多孔介质发出的RF信号接收自 旋的响应。然后处理RF信号，以确定有关多孔介质的成分和/或在多孔介质内所包含或 “包围”的一种或多种流体的信息。
发明内容
本发明一方面涉及一种获得质子密度分布的计算机实施的方法。在一个实施例 中，该方法包括获取来自多孔介质的核磁共振数据；经由全局优化算法反演核磁共振数据， 以确定在多孔介质内的质子密度分布；以及输出所确定的质子密度分布。本发明的另一方面涉及一种获得质子密度分布的计算机实施的方法。在一个实施 例中，该方法包括获取来自多孔介质的核磁共振数据；从核磁共振数据确定多孔介质的质 子密度分布，其中质子密度分布包括由多个预定区域组成的一个或多个频谱，一个或多个 预定区域具有在其中的分布峰值，以及其中确定质子密度分布包括通过使得单个基本分量 适配于在其中具有分布峰值的每个预定区域中的分布而按照非线性基函数对质子密度分 布进行参数化；以及输出所确定的质子密度分布。本发明的另一个方面涉及一种获得有关质子密度分布的信息的方法。在一个实施 例中，该方法包括获取来自介质的核磁共振数据；定义一个实现所获取的核磁共振数据、并 依赖于质子密度分布的m个参数的函数，以使得当接近对于质子密度分布的m个参数的解 时函数被最小化；实施全局优化算法，以确定对于质子密度分布的m个参数的解；以及输出 对于质子密度分布的m个参数的解。本发明的这些和其它目的、特征和特性，以及操作步骤的方法和结构的相关单元 的功能和部件的组合和制造经济性，在参照附图考虑以下的说明和所附权利要求后，将变 得更加清楚，所有的这些形成本技术说明书的一部分，其中相同的附图标记指明在各个图 上的相应的部分。然而，应当明显看到，附图仅仅是用来显示和描述，而不是用来规定本发 明的范围。正如在技术说明书和权利要求中使用的那样，单数形式“一”、“一个”和“该”均 可包括多个单元，除非上下文明确地另外规定。
图1示出了按照本发明的一个实施例的、实施核磁共振以确定有关多孔介质的信 息的方法。图2提供了按照本发明的一个实施例的核磁共振数据的示例图。
图3提供了按照本发明的一个实施例的质子密度分布的示例图。图4示出了具有两个局部最小值和一个绝对最小值的三维表面。图5示出了具有与噪声有关的另一个局部最小值的图4所示表面。图6示出了按照本发明的一个实施例的实施全局优化算法的方法的一个实施例。图7示出了按照本发明的一个实施例的经由全局优化算法反演核磁共振数据，以 确定质子密度的方法。图8示出了质子密度分布。
具体实施例方式核磁共振(“NMR”)技术已广泛应用于测量包含流体的多孔介质的特性(例如地 质构造的岩石物理学特性、活体组织的各种特性、化学化合物的结构特性等等)。某些岩石 物理学特性的例子包括孔隙尺寸、表面积与体积的比率、构造渗透率和毛细管压力。在确定 这些特性时，诸如纵向弛豫时间(“T/’)、横向弛豫时间(“V，)和扩散系数(“D”)那样 的参数常常是感兴趣的。弛豫时间是与核自旋在激励后恢复到它们的平衡位置相关联的时 间。纵向弛豫时间是同自旋与外部静磁场的对准有关的。横向弛豫时间T2是标识在取 向为与主磁场成一个角度的自旋之间发生的相位相干性丢失的时间常数。这种丢失部分是 由于自旋之间的相互作用造成的。扩散系数D是孔隙流体的扩散系数。图1示出了实施NMR以确定有关多孔介质和/或包含在其中的流体的信息的方法 10。应当意识到，虽然方法10的某些方面在这里是在实施NMR以确定有关地质介质的信息 的上下文内描述的，但这并不意图限制本公开内容的范围。本公开内容的范围包括实施NMR 以检验除地质介质以外的介质(例如，活体组织等等)。在一个实施例中，方法10是在一个 或多个处理器上执行的计算机实施方法。在操作步骤12，多孔隙介质和包含在其中的流体由被称为“脉冲序列”的RF脉冲 串扰动。在一个实施例中，该扰动是按照Carr-Purcell-Meiboom-GilKCPMG)脉冲序列完 成的。在某些实施例中，施加其它序列来扰动多孔介质和包含在其中的流体。常规的CPMG脉冲序列包括90度脉冲，后面跟随一系列180度脉冲，其遵循以下的 表达式⑴图其中Ji/2是相对于旋转帧中的基准沿+和-x轴施加的90度脉冲；^代表第k 个回波串的回波间隔TEk的一半；Ji代表在旋转帧中沿y轴施加的180度脉冲；acq代表回 波的获取；nk是在第k个回波串中的回波数。 伴随常规的CPMG脉冲序列一起的是外部磁场梯度Gk和在180度RF脉冲之间施加 的可选的脉冲梯度。脉冲场梯度的宽度用Sk表示，以及在连续的梯度脉冲之间的间隔用 Ak表示。同时，核自旋的系统也受到由外部磁场造成的内部场梯度以及在细颗粒与孔隙流 体之间的敏感度反差作用。在下面，符号g代表核自旋的系统受到的总的磁场梯度。
在操作步骤14，获得来自多孔介质内的RF信号。这些信号通常被称为NMR数据、 “回波数据”或NMR记录数据。正如下面进一步讨论的，这些信号表示响应于操作步骤12的 脉冲序列在多孔介质中产生的一个或多个回波串。关于由在操作步骤14获得的信号传递
的回波串的信息使得能够测量多孔介质和/或包含在其中的流体的一个或多个特性(例如 以上枚举的特性)。在某些实施例中，操作步骤12和14是通过使用已知的NMR记录工具和 /或NMR频谱仪而被执行的。 在操作步骤14获得的对应于脉冲序列的一个回波串的NMR数据M(ti)具有使得 数据M(ti)与参数弛豫时间、T2弛豫时间和扩散系数D相关联的以下的形式
「00241 其中、代表回波之间的时间，、代表回转磁比率，WT代表在回波串之间的等待时 间，a (PARAMETERS)代表给定的参数组(例如，PARAMETERS = T” T2和/或D)的幅值，以 及^代表在NMR数据M(t》中存在的噪声。图2提供在操作步骤14得到的实际的数据 M(t》作为时间的函数的示例图。回到图1，在某些情形下，两个或更多个参数1\，1~2和/或D可以从公式(2)被消 去。例如，当磁场被认为是均勻时，磁场的梯度g等于0，以及被称为扩散核的第三个指数项 广、^^2变为接近于1，这样，公式变为与D无关的。同样地，当等待时间WT足够大时，被 称为极化因子的第二个指数项(丨-广“777”变为接近于1，这样，公式变为与无关。或者， 通过使用有效比率r = I/T2可以使得极化因子变为与无关。在这种情况下，k-和值也 消失。在公式⑵变为与1\和0无关的情况下，在M(ti)与1~2之间的关系可被表示为 其中M(ti)代表所测量的第i个回波幅值(i = 1，. . .，!!)，、代表一组等间隔的n 个衰变时间，以及^是第i个回波的噪声。T”代表一组在对数标尺上等间隔的m个预先 选择的T2弛豫时间，以及a (T2J)是本模型要解决的与弛豫时间T”相关联的T2幅值分布。此后，本公开内容主要描述在其中以NMR数据表示的回波串在计算上与和D无 关(即，M(tD主要取决于T2)的实例中在操作步骤14中获得的NMR数据的处理。这些实 例通常被称为单维NMR解。然而，这并不意图作为限制。应当意识到，这里对于单维NMR解 描述的技术和方法可以很容易地扩展到多于一维的NMR解(例如，在M(t》取决于和/ 或D，以及1~2的实例中)。在操作步骤16，从在操作步骤14得到的NMR数据获得质子密度分布。所获得的质 子密度分布是M(ti)所取决(例如，在单维NMR解中的T2)的参数(例如1\，1~2和/或D)的 函数。例如，图3提供从在图2中表示的NMR数据获得的质子密度分布a (T2J)的示例图。
回过来参照图1，在一个实施例中，从NMR数据获得质子密度分布包括执行反演以 根据在公式(3)中表示的关系来确定分布。例如，对于不带有噪声项、由公式(3)表示 的、对于分布a (T2J)的形式上的数学解通过拉普拉斯逆变换给出。然而，噪声项^使得 这种技术成为无法实施的，因为该反演是非唯一的。另一个方法是把公式(3)重新表达为 最小平方最小化
它是对于使x 2最小的最佳离散分布a (T2J)进行求解。在操作步骤17，输出在操作步骤16得到的质子密度分布。这可包括所述分布的电 子存储、所述分布的电子显示(例如，图形显示)、所述分布的电子传输、或者用于输出通过 计算机实施方法导出的信息的任何其它方法。通常，在操作步骤16被用来确定对于这个最小平方最小化的解的方法包括实施 奇异值分解(“SVD”)。为了给出这种常规技术的概略说明，我们把公式(3)和(4)重写 为
其中 h = M (、), aj = a (T2J)，并且 以矩阵形式表示，公式(5)和(6)变为(7)b = K a+ ￡ ；并且(8) x 2(a) = | K a~b|2 ；其中矩阵元素对应于公式(5)和(6)的分量。为了求解方程(8)，加上条件=0，
这是为了使x2最小而必须满足的条件。这得到(9)Kt (K a-b) = 0 ；其被简化为(10) (KT K) a = KT b ；其被进一步简化为(11) K-a = b;其中上标“T”代表转置。公式(11)是线性方程组，它可以通过使用代数方法，诸 如Gauss-Jordan消元法而进行求解。然而，由于噪声项￡ i，矩阵K是奇异的或接近于奇异。 这主要是因为矩阵K是超定的(即，K是nXm矩阵，其中n>m)，它造成在由噪声给出的不 确定性内的模糊的可能的解a，产生实际上的不确定系统。矩阵K的SVD形式是(12) K = U [diag (Wi) ] VT其中U是[mXn]列正交矩阵，V是[nXn]正交矩阵，以及W是具有正的或接近于 零的元素(奇异值)的[nXn]对角矩阵。与数据的全局噪声电平有关的截止条件被用来 使满足( wfflaxX 6 cutoff的所有对角值为零，其典型值5 cutoff = 10_5。通过使用公式(11) 中K的SVD形式，对于a的解可被表达为 其中向量U⑴，i = 1，...，m表示U的列(每个是长度为n的向量)，以及向量V⑴， i = 1，...，m表示V的列(每个是长度为m的向量)。公式(13)规定适配的参数a是V的 列的线性组合，具有通过形成U的列与加权的数据向量b的点积而得到的系数。尽管截止条件是指平滑掉所有的噪声效应，但实际上它不足以产生平滑的幅值分 布，其中可能出现某些尖峰。在公式(13)上加上被称为“规则化项”的附加项，它用来通过要由搜索方法确定的、匹配于在从操作步骤14得到的数据中的数据噪声电平的某些因子 使得所述分布的模、斜率或曲率最小化。在一个实施例中，作为对于SVD和/或用于对所得到的NMR数据进行反演的其它 常规技术的替代，实施全局优化算法。更具体地，函数x2(a)被看作为m维表面，并且全局 优化算法被实施以确定该表面的绝对最小值。所确定的该表面的绝对最小值的m空间向量 坐标对应于关于M(ti)所取决的参数(例如，T2)的质子密度分布。确定x2(a)的m维表面的绝对最小值并非是不重要的，这部分是因为该表 面很可能包括多个局部最小值。为此，使用局部梯度信息来检测表面最小值-诸如 Levenberg-Marquardt或共轭梯度-的常规优化算法是不现实的(因为它们通常产生局部 最小值而不是绝对最小值)。例如，图4示出了具有局部最小值和绝对最小值的三维表面。 从图4应当意识到，不提供全局优化的优化算法可能变为在一个局部最小值上“固定位”。在一个实施例中，为了减小回到局部最小值的机会，在操作步骤16实施的优化算 法是全局优化算法。正如这里使用的，术语“全局优化算法”是指可被执行以使得函数或函 数组优化到某个预定标准的任何算法或算法组。通常，全局优化算法需要一系列迭代，使得 算法能够收敛到最佳。例如，全局优化算法可包括随机优化(例如，模拟退火、并行回火等 等)、探试法或亚探试法优化(例如，遗传算法、渐近策略等等)和/或其它优化中的一项或 多项。通过实施全局优化算法从NMR数据得到质子密度分布而实现的反演仅仅是x2(a) 函数的绝对最小值的估计。此外，在操作步骤14得到的数据中存在的噪声还可能使得确 定x2(a)最小值进一步复杂化。例如，图5示出了带有与噪声相关联的其它局部最小值的 图4所显示的表面。正如应当从图5意识到的，当接近在不存在噪声时的x2(a)的最小值 (例如图4中的x2(a)的最小值)时，由噪声引起的、包围这个最小值的局部最小值造成对 于x2(a)最小值的等价的(或基本上等价的)解的多重性。这样，经由全局优化算法达到 的单个解不能提供x2(a)函数的绝对最小值，因为全局优化算法会返还与噪声相关联的一 个最小值，而不是实际的最小值。图6示出了实施全局优化算法来考虑(至少在某种程度 上)由噪声引起的等价解的多重性的方法18。在一个实施例中，方法18保证所得到的结果 将提供足够的精度。在某些实例中，方法18可以在方法10的操作步骤16(在图1中显示 并在以上描述)时被实施。然而，这并不意图限于此。方法18包括操作步骤20，在该操作步骤，通过使用全局优化算法得到对于x 2 (a) 的绝对最小值的解。在操作步骤22，更新得到对于x2(a)的解的次数的计数值(例如，加 上1，以反映操作步骤20的解)。在操作步骤24，判断对于x2(a)的绝对最小值的解进行 附加估值是否是必要的。在一个实施例中，这种判断包括将在操作步骤22更新的计数值与 一个预定的阈值进行比较。如果需要附加估值，则方法18回到操作步骤20。如果不需要附 加估值，则方法18进行到操作步骤26。在操作步骤26，对在操作步骤20确定的解的估值 进行合计，以提供对于x2(a)函数的绝对最小值的解的最终估值。例如，所述合计可包括 对所述估值进行平均。通过对解的多个估值进行合计，解的最终估值的精度考虑了(至少 在某种程度上)由噪声引起的等价解的多重性。应当意识到，在方法18中对于x2(a)的 最小值的多个独立的估值的合计在估值被线性参数化的实例中可以有更大的好处，见公式 (4)，因为x2(a)的最小化的非线性参数化会减小噪声对解的影响，见公式18。然而，这并
8不意图限于此，并且方法18可以结合解的非线性参数化一起使用，而不背离本公开内容的 范围。图7示出了按照本发明的一个实施例的经由全局优化算法反演NMR数据以确定质 子密度的方法28 (例如，在方法10的操作步骤16，在方法18的操作步骤20等等)。具体 地，方法28包括实施模拟退火算法。应当意识到，方法28中的模拟退火算法的说明不应当 看作为限制。正如以上阐述的，可以预期，在根据NMR数据确定质子密度时可以实施各种各 样的全局优化算法。通常，模拟退火算法是基于对于物理系统的模拟组成m空间向量a的分量的一 组变量是形成配置空间的假想物理系统的自由度。函数E(a) = x2(a)被认为是系统的能 量，问题被归结为找出系统的最小能量或基础状态配置的问题。众所周知，如果假想物理系 统被加热到非常高的温度T，然后慢慢地冷却到绝对零度(称为退火的过程)，则系统将发 现自己处在基础状态(即，在配置空间中具有最小量的维持能量的配置)。冷却速率必须 足够慢以避免陷入某个亚稳定状态(即，E(a)的局部最小值)。在温度T时，处在具有能量 E(a)的配置中的概率由Boltzmann-Gibbs因子给出(14)exp(_E(a)/T)。Boltzmann-Gibbs因子表明高的能量配置在高的温度T下可以以有限的概率出 现。如果温度降低，则那些高的能量状态变为不太可能的，并且当T — 0时，仅仅接近于E (a) 的最小值的状态才呈现非空概率。因此，在由对于x2(a)的配置空间形成的假想物理系统 内，当a的配置被随机地(或伪随机地)操控并且温度被适当地降低以使得T — 0时，降 低的温度约束a的允许的配置，以使得a的配置慢慢地收敛到在配置空间中最低能量状态 (其中x2(a)被最小化)时的a的配置。在一个实施例中，方法28包括操作步骤30，在该操作步骤，设置当前温度T的初始 值。当前温度T的初始值被设置为足够高以使得方法28能够探查x2(a)函数的相当大量 的配置空间。在操作步骤32，获得对于a的参数(即，分量)的当前的数值组。例如，如果a由 m个参数组成，则在操作步骤32获得对于参数的m个数值。在一个实施例中，所述m个数值 被随机地确定。正如在本公开内容内使用的，术语“随机”包括纯随机和伪随机确定。在操作步骤34，确定对于a的参数的当前的数值组的x 2 (a)函数的数值。在模拟 退火的字典编辑中，这个数值是通过对于a的参数的m个数值的当前的组所描述的配置的 “能量” E (a)。在操作步骤36，获得对于a的m个参数的建议的数值组。在一个实施例中，操作 步骤36包括实施Monte Carlo方法来确定建议的数值组。这并不意图限于此，因为在操作 步骤36确定建议的数值组时可以实施其它算法，诸如混合Monte Carlo方法、Hamiltonian 动态法和/或其它算法。更具体地，在一个实施例中，建议的数值组是通过按照所建议的改 变来改变对于a的m个参数的当前的数值组而被确定的。由于a的参数的数值总是正的， 可以实施适当的随机动态方法来获得建议的数值组。例如，可以实施类似于以下公式的对 数正态随机挪动(15) log(a' ) = log(a) + o G ；其中a’是具有对于a的参数的建议的数值组的m空间向量，G是其所有的分量除
9了一个分量之外均等于零的m空间向量(类似于a)，这个不等于零的分量是按照正态分布 N(0，1)(中心为零且标准差为1的高斯分布)而分布的随机数。应当意识到，尽管这个实施 方案实际上仅将单个数值从当前的数值组变化到建议的数值组，但也可以使用其它的实施 方案来从当前的数值组确定建议的数值组，其中改变一个以上的数值。在操作步骤38，确定对于a的参数的建议的数值组的x 2的数值(即，x2(a’))。 在模拟退火的字典编辑内，在操作步骤38中确定的数值是通过对于a的参数的建议的数值 组所描述的配置的“能量”，它可以被表达为E(a’ )。在操作步骤40，确定接受对于a的m个参数的建议的数值组作为对于a的m个参 数的当前数值组的概率。这个建议的数值组以概率Pa。。(a’ |a)被接受。在一个实施例中， 概率Pa。。(a，a)按照Metropolis接受算法被如下地确定 其中T再次代表当前的温度。在其它的实施例中，可以在操作步骤40实施不同于 Metropolis接受算法的接受算法来确定接受建议的数值组的概率。在操作步骤42，按照在操作步骤40确定的概率，接受或拒绝对于a的m个参数的 建议的数值组。如果接受所建议的数值组，则建议的数值组变为对于a的m个参数的当前 数值组。如果拒绝所建议的数值组，则对于a的m个参数的当前数值组保持不变。在操作步骤44，确定对于当前的温度的新的数值。对于当前的温度的新的数值可 以按照退火计划表被确定，该退火计划表描述当前的温度为已经执行的操作步骤36-42的 迭代次数的函数。这样的退火计划表的例子可被表示为如下(17)T(k) = T(0)e_Ak ；其中k是已经执行的操作步骤36-42的迭代次数。应当意识到，在某些实例中，当 前的温度可能不能对于每次迭代都在操作步骤44被更新。例如，为了节省处理资源和/或 增强算法的速度，在图7中操作步骤44的说明也给出其中当前的温度仅仅在每几次迭代才 被更新的实例。在操作步骤46，判断该算法是否应当执行另一次迭代。在一个实施例中，这种判断 是基于当前的温度是否允许能量的很大的减小。如果确定不应当进行另一次迭代，则对于 a的m个参数的当前的数值组被确定为代表质子密度分布。图8示出了另一种质子密度分布。在以上对方法28的说明中，所确定的m个参数 是线性的(即，幅值a (T2)为1~2的函数)。然而，这并不意图限于此。通过实施全局优化 算法以确定质子密度分布而使得能够进行的增强方案之一是可以按照非线性基函数对所 述解进行参数化。实施非线性基函数来对全局优化算法的解进行参数化(例如，图7中所 示的模拟退火技术)可以减小对于描述所述解所需要的参数的数目。这又意味着用非线性 基函数来对解进行参数化将减小对于达到该解所需要的计算量。参数的数目的减小在多维 解中甚至变为更显著的，因为必须被确定的参数的数目(和/或通过非线性参数化而变为 可能的必须被确定的参数的数目的减小)在添加每个附加维度时被指数地加倍。其它增强 方案也伴随有非线性参数化，其中的某些将在下面进一步讨论。此后，将讨论按照高斯基函 数对解进行非线性参数化，作为非线性参数化的一个例子。这并不意图限于此，因为也可以实施各种其它非线性基函数。例如，非线性参数化可以按照Gamma基函数、B-仿样基函数 或实验确定的基函数。如在图8中可以看到的，质子密度分布通常类似于具有一组连续峰值分布的频 谱。这些峰值分布通常对应于被包围在从其获得NMR数据的多孔介质中的不同类型的流 体。例如，在图8中所示的单维分布中的第一峰值45通常与泥土包围的流体相关联，第二 峰值50通常与毛细管包围的流体相关联，第三峰值52通常与“自由的”流体相关联。因为 峰值48、50和52的通常的形状，每个峰值48、50和52可以由单个高斯分布描述，其参数是 中心、宽度和幅值。在一个实施例中，按照非线性基函数对解进行参数化基本上包括与以上对于通过 执行方法28进行线性参数的解所描述的相同的操作步骤组。回过来参照图7，在操作步骤 32，获得对于参数向量的参数的数值组(例如，按照随机或伪随机确定法)。按照一组队 个高斯基函数分量对解进行参数化使得1~2优化问题能够被表示为(与以上的公式(4)相 比) 其中ak代表基本分量的幅值，bk代表中心，ck代表宽度，以及k代表不同系列的从 1到队的基本分量。从这个表示式可以看到，进行优化的参数可被组合成表示为x= {a” b^Ci，...，aNg，bNg，cNg}的单个m空间向量(其中m = 3 * Ng)。在操作步骤34，从公式(18) 求解对于x的m个参数的当前的数值组的x2(x)。在操作步骤36，确定对于x的m个参数的建议的数值组。由于x的每个分量(“x/’) 只代表几种可能的参数类型中的一种类型，比起以上对于线性参数的操作步骤36的说明， 这个操作步骤可能更多地被涉及到。具体地，由于x的每个分量可以是特定的参数类型(例 如，a, b，c)，对于每个参数类型可以建立可允许的范围，以便把优化以保持在系统的实际上 可能的限制范围内。例如，这可包括对于每种参数类型的最小值、最大值和范围的技术指标 (例如，adel = amax-amin, bdel = bmax-bmin, cdel = cmax-cmin)。从这些规定的约束条件，可以确定 一组m空间约束的向量xmin, x-，和xde1。在一个实施例中，操作步骤36再次包括实施受到约束的Lorentzian随机挪动形 式的Monte Carlo方法，以确定建议的数值组(下面用x’代表)。这个随机挪动可被表示 为(19) x' = x+o L其中o是缩放因子，它可以作为已经完成的操作步骤36-42的迭代次数的函 数而变化，以及L是m空间向量，其中所有的分量除了一个k以外均等于零，k是按照 Lorentzian分布而分布的随机数。必须检查新建议的数值处在对于该参数的允许的范围，否则生成新的建议。在 Lorentzian分布中的尾部大于高斯分布的尾部，所以达到有效的长范围的采样，这可以促 进模拟退火方法28的更佳性能，诸如通常被称为“快速模拟退火”的改进方案。在操作步骤38，从在操作步骤36得到的建议的数值组确定x2(x’)。在操作步骤40，确定接受建议的数值组的概率。例如，可以实施以上对于a和a’描述的Metropolis接 受概率(例如
在操作步骤42，按照在操作步骤40确定的
概率接受或拒绝建议的数值组。操作步骤44和46基本上如以上讨论的那样进行。回到图8，典型地，与泥土包围的流体相关联的第一峰值48位于由质子密度分布 形成的频谱的预定区域(例如，T2 e
)内。同样地，第二峰值50通常位于预定区域 (例如，T2 G [3,33])内，以及第三峰值52通常位于预定区域(例如，T2 G [33,300])内。 另外，有时存在第四峰值(图8中未示出)，它与第二种类型的自由流体相关联。第四峰值 通常也位于在图8中所示的频谱的预定区域(例如，T2e [300，3000]，以ys计)内。应 当意识到，其它区域划分也是可能的，并且这些数值仅仅被提供来用于说明的目的。在一个实施例中，以上讨论的优化方法(其中按照非线性参数组对所述解进行参 数化)，所述解被参数化，以使得单个基本分量(在优化期间)适配于在每个预定区域内的 质子密度分布。这进一步加强了优化，这样，在多孔介质内每种“类型”的流体的影响与单 个基本分量相关联，这使得即使在它穿过典型地与另一种流体类型相关联的频谱的区域的 情形下也能够跟踪对于给定的流体类型的质子密度分布的贡献。例如，这被显示于图8的 区域54。正如可以看到的，与第二峰值50和第三峰值52相关联的各个基本分量重叠。在 更惯用的解中(例如，实施线性参数的解)，在包含第二峰值50和第三峰值52的区域之间 的边界被用来分开对于与第二峰值50和第三峰值52相关联的流体类型的质子密度分布的 贡献。然而，通过把单个非线性基本分量适配于第二峰值50和第三峰值52的每个峰值而 能够进行的这些贡献的离散化使能穿过边界跟踪形成第二峰值50和第三峰值52的流体类 型的贡献。虽然本发明是根据当前认为的最实际和优选的实施例，为了说明的目的而详细地 描述的，但应当看到，这样的细节仅仅是为了该目的，本发明不限于所公开的实施例，而是 相反，本发明意图覆盖属于所附权利要求的主旨和范围内的修改方案和等价的设计。例如， 应当理解，本发明预期，在可能的范围内，任何实施例的一个或多个特性可以与任何另外的 实施例的一个或多个特性相组合。
权利要求
一种获得质子密度分布的计算机实施的方法，该方法包括获取来自多孔介质的核磁共振数据；经由全局优化算法反演核磁共振数据，以确定在多孔介质内的质子密度分布；以及输出所确定的质子密度分布。
2.权利要求1的方法，其中按照非线性基函数对质子密度分布进行参数化。
3.权利要求2的方法，其中非线性基函数包括高斯基函数、Gamma基函数、B-仿样基函 数或实验确定的基函数中的一个或多个基函数。
4.权利要求2的方法，其中质子密度分布包括由多个预定区域组成的一个或多个频 谱，一个或多个预定区域具有在其中的分布峰值，以及其中经由全局优化算法反演核磁共 振数据以确定质子密度分布包括使得单个基本分量适配于在其中具有分布峰值的每个预 定区域中的分布。
5.权利要求4的方法，其中各个预定区域被确定为对应于在多孔介质内的不同的流体 类型，以使得给定预定区域对应于多孔介质内的相应的流体类型。
6.权利要求1的方法，其中所述全局优化算法包括模拟退火、遗传算法、渐近策略或并 行回火中的一项或多项。
7.权利要求6的方法，其中通过模拟退火算法实施的采样技术包括MonteCarlo方法、 混合Monte Carlo方法、Hamiltonian动态法或随机挪动方法中的一项或多项。
8.权利要求1的方法，其中质子密度分布是η维分布，其中η大于2。
9.权利要求1的方法，其中经由全局优化算法反演核磁共振数据以确定质子密度分布 包括实施全局优化算法两次或多次，以确定对于质子密度分布的多个解；以及对于质子密度分布的多个解求平均。
10.权利要求9的方法，其中质子密度分布是幅值的线性表示。
11.一种获得质子密度分布的计算机实施的方法，该方法包括获取来自多孔介质的核磁共振数据；由核磁共振数据确定多孔介质的质子密度分布，其中质子密度分布包括由多个预定区 域组成的一个或多个频谱，一个或多个预定区域具有在其中的分布峰值，以及其中确定质 子密度分布包括通过使得单个基本分量适配于在其中具有分布峰值的每个预定区域中的 分布而按照非线性基函数对质子密度分布进行参数化；以及输出所确定的质子密度分布。
12.权利要求11的方法，其中非线性基函数包括高斯基函数、Gamma基函数、B-仿样基 函数或实验确定的基函数中的一个或多个基函数。
13.权利要求11的方法，其中各个预定区域被确定为对应于在多孔介质内的不同的流 体类型，以使得给定的预定区域对应于多孔介质内的相应的流体类型。
14.权利要求13的方法，其中所述预定区域包括对应于泥土包围的流体的预定区域、 对应于毛细管包围的流体的预定区域或者对应于一种或多种类型的自由的流体的区域的 一个或多个区域。
15.权利要求11的方法，其中由核磁共振数据确定介质的质子密度分布包括实施全局 优化算法。
16.权利要求15的方法，其中所述全局优化算法包括模拟退火、遗传算法、渐近策略或 并行回火中的一项或多项。
17.权利要求16的方法，其中通过模拟退火算法实施的采样技术包括MonteCarlo方 法、混合Monte Carlo方法、Hamiltonian动态法或随机挪动方法中的一项或多项。
18.权利要求11的方法，其中质子密度分布是n维分布，其中n大于2。
19.一种获得有关质子密度分布的信息的方法，该方法包括获取来自介质的核磁共振数据；定义一个实现所获取的核磁共振数据并与质子密度分布的m个参数有关的函数，以使 得当接近对于质子密度分布的m个参数的解时该函数被最小化；实施全局优化算法，以确定对于质子密度分布的m个参数的解；以及输出对于质子密度分布的m个参数的解。
20.权利要求19的方法，其中实施全局优化算法以确定对于质子密度分布的m个参数 的解包括(a)设置用于算法的当前温度的初始值；(b)随机地确定对于m个参数的当前数值组；(c)确定对于m个参数的当前数值组的函数值；(d)通过随机地调节在对于m个参数的当前数值组中的一个或多个数值而确定对于m 个参数的建议的数值组；(e)确定对于m个参数的建议的数值组的函数值；(f)根据对于m个参数的当前数值组的函数值、对于建议的参数组的函数值和当前的 温度，计算接受对于m个参数的建议的数值组作为对于m个参数的当前数值组的概率；(g)根据在(f)步骤计算的概率，接受或拒绝对于m个参数的建议的数值组作为对于m 个参数的当前数值组；以及(h)确定用于该算法的对于当前的温度的新的数值，使得当前的温度按该算法的迭代 的函数降低。
21.权利要求20的方法，其中(a)-(h)中的每个步骤以所述的次序执行，并且该方法还 包括在步骤(h)之后返回到步骤(d)。
22.权利要求19的方法，其中m个参数是非线性的。
23.权利要求19的方法，其中实施全局优化算法以确定对于m个参数的解包括实施全局优化算法两次或多次，以确定对于m个参数的多个解；以及对于m个参数的多个解求平均。
24.权利要求23的方法，其中m个参数是线性的。
25.权利要求19的方法，其中质子密度分布是n维分布，其中n大于2。
全文摘要
一种计算机实施的方法，使得能够获得质子密度分布。在一个实施例中，该方法包括获取来自多孔介质的核磁共振数据；经由全局优化算法反演核磁共振数据，以确定在多孔介质内的质子密度分布；以及输出所确定的质子密度分布。
文档编号G01V3/32GK101896834SQ200880120342
公开日2010年11月24日 申请日期2008年12月4日 优先权日2007年12月12日
发明者R·萨拉扎尔-蒂奥, 孙伯勤 申请人:雪佛龙美国公司
最新回复(0)