基于LBM模拟低场核磁共振分析流体的方法及装置

摘要

本发明实施例提供一种基于LBM模拟低场核磁共振分析流体的方法及装置。该方法包括：获取样品岩石的数字岩心；依据所述数字岩心获得所述样品岩石中孔隙和岩石骨架对应的三维数据体；依据所述三维数据体建立单相多组分模型，并基于所述单相多组分模型模拟所述孔隙中流体的低场核磁共振过程，分析获得所述流体的组分浓度，根据所述流体的组分浓度获得磁化矢量的模值。本发明实施例通过单相多组分模型模拟样品岩石孔隙中流体的低场核磁共振过程，分析所述流体的组分，通过所述流体的组分浓度获得磁化矢量的模值，实现了模拟低场核磁共振分析流体分子的技术，为研究低场核磁共振分析流体分子提供了理论基础。

说明书

技术领域

本发明实施例涉及石油勘探技术领域，尤其涉及一种基于LBM模拟低场 核磁共振分析流体的方法及装置。

背景技术

核磁共振测井是一种适用于裸眼井的测井新技术，是目前唯一可以直接 测量任意岩性储集层自由流体渗流体积特性的测井方法，有明显的优越性。 核磁共振技术利用原子核的顺磁性以及原子核之间相互作用的外加磁场实现 测井。

低场核磁共振是一种非侵入性技术，常用于测量储集层岩心的孔隙度、 孔隙大小分布、孔隙中流体类型，以及流体饱和度，同时还用于测量含油岩 心的物理性质，观测分析流体分子的动态行为，为提高油气采收率、优化采 油过程提供理论依据和实践可能。

但是，现有技术中没有能够模拟低场核磁共振分析流体分子的技术，导 致研究低场核磁共振分析流体分子缺少理论基础。

发明内容

本发明实施例提供一种基于LBM模拟低场核磁共振分析流体的方法及装 置，以实现模拟低场核磁共振分析流体分子的技术，提高研究低场核磁共振 分析流体分子的理论基础。

本发明实施例的一个方面是提供一种基于LBM模拟低场核磁共振分析流 体的方法，包括：

获取样品岩石的数字岩心；

依据所述数字岩心获得所述样品岩石中孔隙和岩石骨架对应的三维数据 体；

依据所述三维数据体建立单相多组分模型，并基于所述单相多组分模型 模拟所述孔隙中流体的低场核磁共振过程，分析获得所述流体的组分浓度， 根据所述流体的组分浓度获得磁化矢量的模值。

本发明实施例的另一个方面是提供一种基于LBM模拟低场核磁共振分析 流体的装置，包括：

数字岩心获取模块，用于获取样品岩石的数字岩心；

三维数据体获取模块，用于依据所述数字岩心获得所述样品岩石中孔隙 和岩石骨架对应的三维数据体；

建模分析模块，用于依据所述三维数据体建立单相多组分模型，并基于 所述单相多组分模型模拟所述孔隙中流体的低场核磁共振过程，分析获得所 述流体的组分浓度，根据所述流体的组分浓度获得磁化矢量的模值。

本发明实施例提供的基于LBM模拟低场核磁共振分析流体的方法及装 置，通过单相多组分模型模拟样品岩石孔隙中流体的低场核磁共振过程，分 析所述流体的组分，通过所述流体的组分浓度获得磁化矢量的模值，实现了 模拟低场核磁共振分析流体分子的技术，为研究低场核磁共振分析流体分子 提供了理论基础。

附图说明

图1为本发明实施例提供的基于LBM模拟低场核磁共振分析流体的方法 流程图；

图2为本发明实施例提供的CPMG脉冲序列示意图；

图3为本发明实施例提供的基于LBM模拟低场核磁共振分析流体的装置 的结构图；

图4为本发明另一实施例提供的基于LBM模拟低场核磁共振分析流体的 装置的结构图。

具体实施方式

图1为本发明实施例提供的基于LBM模拟低场核磁共振分析流体的方法 流程图。本发明实施例针对现有技术中没有能够模拟低场核磁共振分析流体 分子的技术，提供了基于LBM模拟低场核磁共振分析流体的方法，该方法具 体步骤如下：

步骤S101、获取样品岩石的数字岩心；

本发明实施例选择一块岩石为样品，选用X光全岩心CT、微米X光CT、 纳米X光CT、聚焦离子束扫描电镜和岩心反射光谱扫描仪中的任一一种设备 获取样品岩石的数字岩心，该数字岩心具体为样品岩石切片的扫描图像，且 不同的设备其最高分辨率，以及适用的样品尺寸、适用的样品特性如表1所 示：

表1

步骤S102、依据所述数字岩心获得所述样品岩石中孔隙和岩石骨架对应 的三维数据体；

所述依据所述数字岩心获得所述样品岩石中孔隙和岩石骨架对应的三维 数据体包括：对所述数字岩心进行三维重建获得三维数据体；对所述三维数 据体进行滤波处理，依据阈值对所述三维数据体进行二值化处理区分出所述 样品岩石中的孔隙和岩石骨架对应的三维数据体。

步骤S101获取的数字岩心的格式一般为tiff图片，利用二维图片的三 维重建技术对tiff图片进行重建获得三维数据体，二维图片的三维重建技术 具体采用现有技术中的任一种。对获得的三维数据体进行滤波处理以消除噪 点，依据预设的阈值对滤波后的三维数据体进行二值化处理，将三维数据体 中大于阈值的点判断为样品岩石的孔隙、小于阈值的点判断为样品岩石的岩 石骨架，或者将三维数据体中小于阈值的点判断为样品岩石的孔隙、大于阈 值的点判断为样品岩石的岩石骨架，经过二值化处理便可区分出所述样品岩 石中的孔隙和岩石骨架。

步骤S103、依据所述三维数据体建立单相多组分模型，并基于所述单相 多组分模型模拟所述孔隙中流体的低场核磁共振过程，分析获得所述流体的 组分浓度，根据所述流体的组分浓度获得磁化矢量的模值。

利用单相多组分模型具体为格子Boltzmann方法(LatticeBoltzmann Method,简称LBM)的Shan-Doolen模型模拟所述孔隙中流体的低场核磁共振 过程，在该过程中分析流体的组分，具体对流体的组分进行标记，获得被标 记的组分和未被标记的组分，其中被标记的组分是流体中能够对CPMG脉冲序 列产生反应的分子，未被标记的组分是流体中不能对CPMG脉冲序列产生反应 的分子，通过所述流体的组分浓度获得磁化矢量的模值，具体用被标记的组 分的浓度表示磁化矢量的模值。

本发明实施例通过单相多组分模型模拟样品岩石孔隙中流体的低场核磁 共振过程，分析所述流体的组分，通过所述流体的组分浓度获得磁化矢量的 模值，实现了模拟低场核磁共振分析流体分子的技术，为研究低场核磁共振 分析流体分子提供了理论基础。

在上述实施例的基础上，所述依据所述三维数据体建立单相多组分模型， 并基于所述单相多组分模型模拟所述孔隙中流体的低场核磁共振过程，分析 获得所述流体的组分浓度，根据所述流体的组分浓度获得磁化矢量的模值， 包括：

对所述三维数据体的X轴、Y轴和Z轴进行离散划分获得多个晶格；

针对所述多个晶格中的任一晶格，将所述晶格内的流体的组分划分为第 一类组分g和第二类组分f，所述第一类组分g为所述晶格内在XY平面内对 应有磁化矢量分量的流体分子，所述第二类组分f为所述晶格内除所述第一 类组分g之外的流体分子；

依据公式(1)和(2)分别获得所述第一类组分g在平衡态的分布函数和所述第二类组分f在平衡态的分布函数f_i^eq：

$g_i^{eq}=w_i{\rho }_g[1+\frac{e_i*u_g^{eq}}{c_s^2}+\frac{{(e_i*u_g^{eq})}^2}{2c_s^4}-\frac{u_g^{eq}*u_g^{eq}}{2c_s^2}]---\left(1\right)$

$f_i^{eq}=w_i{\rho }_f[1+\frac{e_i*u_f^{eq}}{c_s^2}+\frac{{(e_i*u_f^{eq})}^2}{2c_s^4}-\frac{u_f^{eq}*u_f^{eq}}{2c_s^2}]---\left(2\right)$

其中，i∈[1,2,....,N],N-1表示与所述晶格相邻的所有晶格的数目，w_i表示 权重系数，ρ_g表示所述第一类组分g的浓度，ρ_f表示所述第二类组分f的浓 度，c_s表示晶格声速，e_i表示离散速度单位矢量，表示所述第一类组分g在 平衡态的宏观速度，表示所述第二类组分f在平衡态的宏观速度；

依据公式(3)和(4)分别获得所述第一类组分g的宏观速度u_g和所述 第二类组分f的宏观速度u_f：

$\underset{i=1}{\overset{N}{\Sigma }}{c}_i{g}_i(x,t)={\rho }_g{u}_g---\left(3\right)$

$\underset{i=1}{\overset{N}{\Sigma }}{c}_i{f}_i(x,t)={\rho }_f{u}_f---\left(4\right)$

其中，c_i为常数，g_i(x,t)表示所述第一类组分g在位置x、时刻t的分布函 数，f_i(x,t)表示所述第二类组分f在位置x、时刻t的分布函数；

依据公式(5)和(6)分别获得所述第一类组分g在位置x、时刻t的平 衡态的宏观速度和所述第二类组分f在位置x、时刻t的平衡态的宏观 速度

$u_g^{eq}(x,t)=u_g(x,t)+{\mathrm{\tau F}}_g(x,t)/{\rho }_g(x,t)---\left(5\right)$

$u_f^{eq}(x,t)=u_f(x,t)+{\mathrm{\tau F}}_f(x,t)/{\rho }_f(x,t)---\left(6\right)$

其中，u_g(x,t)表示所述第一类组分g在位置x、时刻t的宏观速度，ρ_g(x,t) 表示所述第一类组分g在位置x、时刻t的浓度， F_g(x,t)＝-ψ_g(x)∑G[ψ_g(x+e_i)-ψ_f(x+e_i)]e_i，ψ_g(x)＝exp(-1/ρ_g(x,t))， ψ_g(x+e_i)＝exp(-1/ρ_g(x+e_i,t))，τ和G为常数，u_f(x,t)表示所述第二类组分f在 位置x、时刻t的宏观速度，ρ_f(x,t)表示所述第二类组分f在位置x、时刻t的 浓度，F_f(x,t)＝-ψ_f(x)∑G[ψ_f(x+e_i)-ψ_g(x+e_i)]e_i，ψ_f(x)＝exp(-1/ρ_f(x,t))， ψ_f(x+e_i)＝exp(-1/ρ_f(x+e_i,t))；

依据公式(7)和(8)分别获得所述第一类组分g在位置x、时刻t的平 衡态的分布函数和所述第二类组分f在位置x、时刻t的平衡态的分布 函数

$g_i^{eq}(x,t)=w_i{\rho }_g(x,t)[1+\frac{e_i*u_g^{eq}(x,t)}{c_s^2}+\frac{{(e_i*u_g^{eq}(x,t\left)\right)}^2}{2c_s^4}-\frac{u_g^{eq}(x,t)*u_g^{eq}(x,t)}{2c_s^2}]---\left(7\right)$

$f_i^{eq}(x,t)=w_i{\rho }_f(x,t)[1+\frac{e_i*u_f^{eq}(x,t)}{c_s^2}+\frac{{(e_i*u_f^{eq}(x,t\left)\right)}^2}{2c_s^4}-\frac{u_f^{eq}(x,t)*u_f^{eq}(x,t)}{2c_s^2}]---\left(8\right)$

依据公式(9)和(10)分别获得所述第一类组分g在位置(x+c△t)、时刻 (t+△t)的分布函数g_i(x+c△t,t+△t)和所述第二类组分f在位置(x+c△t)、时刻 (t+△t)的分布函数f_i(x+c△t,t+△t)：

$r(x+c\Delta t,t+\Delta t)=-\frac{1}{\tau }[g_i^{eq}-g_i^{eq}(x,t)]+g_i^{eq}---\left(9\right)$

$f_i(x+c\Delta t,t+\Delta t)=-\frac{1}{\tau }[f_i^{eq}-f_i^{eq}(x,t)]+f_i^{eq}---\left(10\right)$

其中，$c^2=c_s^2/3;$

依据公式(11)和(12)分别获得所述第一类组分g的浓度ρ_g和所述第 二类组分f的浓度ρ_f：

$\underset{i=1}{\overset{N}{\Sigma }}{g}_i(x+c\Delta t,t+\Delta t)={\rho }_g---\left(11\right)$

$\underset{i=1}{\overset{N}{\Sigma }}{f}_i(x+c\Delta t,t+\Delta t)={\rho }_f---\left(12\right)$

获得所述磁化矢量的模值，所述磁化矢量的模值为所述第一类组分g的 浓度ρ_g。

在本发明实施例中，对时间进行离散化处理以使流体的时间步长为△t， 同时对三维数据体的X轴、Y轴和Z轴进行离散划分获得多个晶格，该晶格 优选为立方体，即孔隙中的流体具体以晶格的形式进行划分，并且将流体的 组分划分为第一类组分g和第二类组分f，所述第一类组分g为所述晶格内在 XY平面内对应有磁化矢量分量的流体分子即能够对CPMG脉冲序列产生反应 的被标记的组分，所述第二类组分f为所述晶格内除所述第一类组分g之外的 流体分子即不能对CPMG脉冲序列产生反应的未被标记的组分，在初始状态 下，合理假设第一类组分g的浓度初始值为1-10^-10，第二类组分f的浓度初始 值为10^-10，且第一类组分g的浓度初始值和第二类组分f的浓度初始值不限于 此，所述第一类组分g在平衡态的宏观速度为0，第二类组分f在平衡态的 宏观速度为0，通过公式(1)计算获得第一类组分g在平衡态的分布函数 此处ρ_g＝1-10^-10，通过公式(2)计算获得第二类组分f在平衡态 的分布函数f_i^eq＝w_iρ_f，此处ρ_f＝10^-10，其中，i∈[1,2,....,N],N-1表示与所述晶 格相邻的所有晶格的数目，w_i表示权重系数，ρ_g表示所述第一类组分g的浓 度，ρ_f表示所述第二类组分f的浓度，c_s表示晶格声速，e_i表示离散速度单 位矢量，且$e_i=\left(\begin{array}{ccccccccccccccccccc}0& 1& -1& 0& 0& 0& 0& 1& -1& -1& 1& 1& -1& -1& 1& 0& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0& 1& -1& -1& 1\\ 0& 0& 0& 0& 0& 1& -1& 0& 0& 0& 0& 1& 1& -1& -1& 1& 1& -1& -1\end{array}\right),$表示所述 第一类组分g在平衡态的宏观速度，表示所述第二类组分f在平衡态的宏 观速度，第一类组分g的浓度初始值为ρ_g＝1-10^-10，第二类组分f的浓度的初 始值为ρ_f＝10^-10。

公式(3)中第一类组分g的宏观速度u_g的初始值为0，公式(4)中第二 类组分f的宏观速度u_f的初始值为0，则公式(5)中u_g(x,t)初始值为0， ρ_g(x,t)＝ρ_g＝1-10^-10；公式(6)中u_f(x,t)初始值为0， ρ_f(x,t)＝ρ_f＝10^-10；其中， F_g(x,t)＝-ψ_g(x)∑G[ψ_g(x+e_i)-ψ_f(x+e_i)]e_i，ψ_g(x)＝exp(-1/ρ_g(x,t))， ψ_g(x+e_i)＝exp(-1/ρ_g(x+e_i,t))，τ和G为常数， F_f(x,t)＝-ψ_f(x)∑G[ψ_f(x+e_i)-ψ_g(x+e_i)]e_i，ψ_f(x)＝exp(-1/ρ_f(x,t))， ψ_f(x+e_i)＝exp(-1/ρ_f(x+e_i,t))；通过公式(5)和(6)分别获得所述第一类组 分g在位置x、时刻t的平衡态的宏观速度和所述第二类组分f在位置 x、时刻t的平衡态的宏观速度

将公式(5)计算获得的和ρ_g(x,t)＝ρ_g＝1-10^-10带入公式(7)获得 第一类组分g在位置x、时刻t的平衡态的分布函数将公式(6)计算 获得的和ρ_f(x,t)＝ρ_f＝10^-10带入公式(8)获得第二类组分f在位置x、 时刻t的平衡态的分布函数f_i^eq(x,t)；其中，c_s表示晶格声速，w_i表示权重系数， e_i表示离散速度单位矢量。

将公式(1)获得的和公式(7)获得的带入公式(9)获得第 一类组分g在位置(x+c△t)、时刻(t+△t)的分布函数g_i(x+c△t,t+△t)，将公式(2) 获得的f_i^eq和公式(8)获得的f_i^eq(x,t)带入公式(10)获得第二类组分f在位 置(x+c△t)、时刻(t+△t)的分布函数f_i(x+c△t,t+△t)；其中，

将公式(9)获得的g_i(x+c△t,t+△t)带入公式(11)获得第一类组分g的浓 度ρ_g，将公式(10)获得的f_i(x+c△t,t+△t)带入公式(12)获得第二类组分f 的浓度ρ_f。并且第一类组分g的浓度ρ_g即为磁化矢量的模值。

本发明实施例提供了单相多组分模型即Shan-Doolen模型计算磁化矢量 的模值的方法。

图2为本发明实施例提供的CPMG脉冲序列示意图。在上述实施例的基础 上，所述依据公式(11)和(12)分别获得所述第一类组分g的浓度ρ_g和所 述第二类组分f的浓度ρ_f之后，还包括：

依据公式(13)和(14)分别计算所述流体在表面弛豫作用下，所述第 一类组分g的浓度ρ_g随时间的变化率s_g和所述第二类组分f的浓度ρ_f随时间 的变化率s_f，并依据公式(15)和(16)分别获得表面弛豫作用后的所述第 一类组分g的浓度ρ′_g和所述第二类组分f的浓度ρ′_f：

$\frac{\partial {\rho }_g}{\partial t}{|}_{surface}=-s_g{\rho }_g---\left(13\right)$

$\frac{\partial {\rho }_f}{\partial t}{|}_{surface}=-s_f{\rho }_f---\left(14\right)$

ρ′_g＝ρ_g-s_gρ_g(15)

ρ′_f＝ρ_f-s_fρ_f(16)

依据公式(17)和(18)分别计算所述流体在体弛豫作用下，所述第一 类组分g的浓度ρ′_g随时间的变化率b_g和所述第二类组分f的浓度ρ′_f随时间的 变化率b_f，并依据公式(19)和(20)分别获得体弛豫作用后的所述第一类 组分g的浓度ρ″_g和所述第二类组分f的浓度ρ″_f：

$\frac{\partial {\rho }_g}{\partial t}{|}_{body}=-b_g{\rho }_g---\left(17\right)$

$\frac{\partial {\rho }_f}{\partial t}{|}_{body}=-b_f{\rho }_f---\left(18\right)$

ρ″_g＝ρ′_g-b_gρ′_g(19)

ρ″_f＝ρ′_f-b_fρ′_f(20)

依据公式(21)和(22)分别计算所述流体的扩散系数D_g和D_f，并依据 公式(23)和(24)分别获得流体扩散后的所述第一类组分g的浓度ρ″′_g和所 述第二类组分f的浓度ρ″′_f：

$\frac{\partial {\rho }_g}{\partial t}=D_g{▿}^2{\rho }_g---\left(21\right)$

$\frac{\partial {\rho }_f}{\partial t}=D_f{▿}^2{\rho }_f---\left(22\right)$

ρ″′_g＝ρ″_g-D_gρ″_g(23)

ρ″′_f＝ρ″_f-D_fρ″_f(24)

将ρ″′_g作为公式(1)中的ρ_g、ρ″′_f作为公式(2)中的ρ_f，将作为公 式(1)中的将作为公式(2)中的循环执行公式(1)-(24)。

在本发明实施例中，第一类组分g的浓度与第二类组分f的浓度之和为 1，在流体中加如图2所示的CPMG脉冲序列，其中，CPMG脉冲序列的90度 脉冲是指以氢原子的拉莫尔进动频率将磁化矢量从z轴正方向搬转到xy平面 内所需要的脉冲强度和时间；CPMG脉冲序列180度脉冲是指以氢原子的拉莫 尔进动频率将磁化矢量从xy平面内轴对称反转所需要的脉冲强度和时间。在 流体中加CPMG脉冲序列的90度脉冲之前，设置流体中第一类组分g的浓度 初始值为1-10^-10，第二类组分f的浓度初始值为10^-10，在流体中加CPMG脉冲 序列的90度脉冲之后，设置流体中第一类组分g的浓度初始值为10^-10，第二 类组分f的浓度初始值为1-10^-10。

流体在表面弛豫作用下，第一类组分g的浓度ρ_g会随时间发生变化，变 化率为s_g，第二类组分f的浓度ρ_f会随时间发生变化，变化率为s_f，将公式 (11)获得的ρ_g带入公式(13)中可获得变化率s_g，将ρ_g和s_g带入公式(15) 获得表面弛豫作用后，第一类组分g的浓度ρ′_g；同理，将公式(12)获得的ρ_f 带入公式(14)中可获得变化率s_f，将ρ_f和s_f带入公式(16)获得表面弛豫 作用后，第二类组分f的浓度ρ′_f。

流体在体弛豫作用下，第一类组分g的浓度ρ_g会随时间发生变化，变化 率为b_g，第二类组分f的浓度ρ_f会随时间发生变化，变化率为b_f，将公式(11) 获得的ρ_g带入公式(17)获得b_g，将公式(12)获得的ρ_f带入公式(18)获 得b_f，将公式(15)获得的ρ′_g和公式(17)获得的b_g带入公式(19)获得弛 豫作用后的所述第一类组分g的浓度ρ″_g；同理，获得所述第二类组分f的浓 度ρ″_f。

通过公式(21)和(22)分别计算所述流体的扩散系数D_g和D_f，并依据 公式(23)和(24)分别获得流体扩散后的所述第一类组分g的浓度ρ″′_g和所 述第二类组分f的浓度ρ″′_f。

将流体扩散后的所述第一类组分g的浓度ρ″′_g作为公式(1)中的ρ_g，将公 式(5)获得的作为公式(1)中的将流体扩散后的所述第二类组 分f的浓度ρ″′_f作为公式(2)中的ρ_f，将公式(6)获得的作为公式(2) 中的循环执行公式(1)-(24)，在循环过程中第一类组分g的浓度初 始值和第二类组分f的浓度初始值可以随时改变。

所述根据所述流体的组分浓度获得磁化矢量的模值之后，还包括：在所 述流体上施加压力梯度，以使所述流体的流速发生变化；根据所述流体的组 分浓度获得不同流体流速下所述磁化矢量的模值。

在本发明实施例中，流体上施加压力梯度为F_ex(x,t)，则流体的流速将发 生变化，上述公式(5)变形为$u_g^{eq}(x,t)=u_g(x,t)+\tau \left(F_g\right(x,t)+F_{ex}(x,t))/{\rho }_g(x,t),$上述 公式(6)变形为$u_f^{eq}(x,t)=u_f(x,t)+\tau \left(F_f\right(x,t)+F_{ex}(x,t))/{\rho }_f(x,t),$其余公式和执行步 骤不变，通过计算第一类组分g的浓度ρ_g获得磁化矢量的模值。

所述依据所述数字岩心获得所述样品岩石中孔隙和岩石骨架对应的三维 数据体之后，还包括：利用氦孔隙度测量仪测量所述样品岩石的测量孔隙度； 依据所述三维数据体计算所述样品岩石的理论孔隙度；若所述理论孔隙度与 所述测量孔隙度不相等，则调节所述阈值，再次执行所述依据阈值对所述三 维数据体进行二值化处理区分出所述样品岩石中的孔隙和岩石骨架对应的三 维数据体的步骤，直至所述理论孔隙度与所述测量孔隙度相等。

在步骤S102之后，需要对预设的阈值进行校正，具体校正方法为：利用 氦孔隙度测量仪测量所述样品岩石的测量孔隙度；依据所述三维数据体计算 所述样品岩石的理论孔隙度；比较样品岩石的测量孔隙度和理论孔隙度是否 相等，若所述理论孔隙度与所述测量孔隙度不相等，则调节所述阈值，继续 依据所述三维数据体计算所述样品岩石的理论孔隙度，直至所述理论孔隙度 与所述测量孔隙度相等，此时的阈值为校准的阈值。

本发明实施例通过循环执行单相多组分模型模拟孔隙中流体的低场核磁 共振过程，获得随时间变化的磁化矢量的模值，以及通过流体上施加压力梯 度，分析不同流体流速下磁化矢量的模值，提高了模拟低场核磁共振分析流 体分子的能力，进一步提高了低场核磁共振分析流体分子的理论基础；另外， 通过对二值化处理所用的阈值进行校正，提高了孔隙度的测量精度。

图3为本发明实施例提供的基于LBM模拟低场核磁共振分析流体的装置 的结构图。本发明实施例提供的基于LBM模拟低场核磁共振分析流体的装置 可以执行基于LBM模拟低场核磁共振分析流体的方法实施例提供的处理流 程，如图3所示，基于LBM模拟低场核磁共振分析流体的装置30包括数字岩 心获取模块31、三维数据体获取模块32和建模分析模块33，其中，数字岩 心获取模块31用于获取样品岩石的数字岩心；三维数据体获取模块32用于 依据所述数字岩心获得所述样品岩石中孔隙和岩石骨架对应的三维数据体； 建模分析模块33用于依据所述三维数据体建立单相多组分模型，并基于所述 单相多组分模型模拟所述孔隙中流体的低场核磁共振过程，分析获得所述流 体的组分浓度，根据所述流体的组分浓度获得磁化矢量的模值。

本发明实施例通过单相多组分模型模拟样品岩石孔隙中流体的低场核磁 共振过程，分析所述流体的组分，通过所述流体的组分浓度获得磁化矢量的 模值，实现了模拟低场核磁共振分析流体分子的技术，为研究低场核磁共振 分析流体分子提供了理论基础。

图4为本发明另一实施例提供的基于LBM模拟低场核磁共振分析流体的 装置的结构图。在上述实施例的基础上，建模分析模块33具体用于对所述三 维数据体的X轴、Y轴和Z轴进行离散划分获得多个晶格；

针对所述多个晶格中的任一晶格，将所述晶格内的流体的组分划分为第 一类组分g和第二类组分f，所述第一类组分g为所述晶格内在XY平面内对 应有磁化矢量分量的流体分子，所述第二类组分f为所述晶格内除所述第一 类组分g之外的流体分子；

依据公式(1)和(2)分别获得所述第一类组分g在平衡态的分布函数和所述第二类组分f在平衡态的分布函数f_i^eq：

$g_i^{eq}=w_i{\rho }_g[1+\frac{e_i*u_g^{eq}}{c_s^2}+\frac{{(e_i*u_g^{eq})}^2}{2c_s^4}-\frac{u_g^{eq}*u_g^{eq}}{2c_s^2}]---\left(1\right)$

$f_i^{eq}=w_i{\rho }_f[1+\frac{e_i*u_f^{eq}}{c_s^2}+\frac{{(e_i*u_f^{eq})}^2}{2c_s^4}-\frac{u_f^{eq}*u_f^{eq}}{2c_s^2}]---\left(2\right)$

其中，i∈[1,2,....,N],N-1表示与所述晶格相邻的所有晶格的数目，w_i表示 权重系数，ρ_g表示所述第一类组分g的浓度，ρ_f表示所述第二类组分f的浓 度，c_s表示晶格声速，e_i表示离散速度单位矢量，表示所述第一类组分g在 平衡态的宏观速度，表示所述第二类组分f在平衡态的宏观速度；

依据公式(3)和(4)分别获得所述第一类组分g的宏观速度u_g和所述 第二类组分f的宏观速度u_f：

$\underset{i=1}{\overset{N}{\Sigma }}{c}_i{g}_i(x,t)={\rho }_g{u}_g---\left(3\right)$

$\underset{i=1}{\overset{N}{\Sigma }}{c}_i{f}_i(x,t)={\rho }_f{u}_f---\left(4\right)$

其中，c_i为常数，g_i(x,t)表示所述第一类组分g在位置x、时刻t的分布函 数，f_i(x,t)表示所述第二类组分f在位置x、时刻t的分布函数；

依据公式(5)和(6)分别获得所述第一类组分g在位置x、时刻t的平 衡态的宏观速度和所述第二类组分f在位置x、时刻t的平衡态的宏观 速度

$u_g^{eq}(x,t)=u_g(x,t)+{\mathrm{\tau F}}_g(x,t)/{\rho }_g(x,t)---\left(5\right)$

$u_f^{eq}(x,t)=u_f(x,t)+{\mathrm{\tau F}}_f(x,t)/{\rho }_f(x,t)---\left(6\right)$

其中，u_g(x,t)表示所述第一类组分g在位置x、时刻t的宏观速度，ρ_g(x,t) 表示所述第一类组分g在位置x、时刻t的浓度， F_g(x,t)＝-ψ_g(x)∑G[ψ_g(x+e_i)-ψ_f(x+e_i)]e_i，ψ_g(x)＝exp(-1/ρ_g(x,t))， ψ_g(x+e_i)＝exp(-1/ρ_g(x+e_i,t))，τ和G为常数，u_f(x,t)表示所述第二类组分f在 位置x、时刻t的宏观速度，ρ_f(x,t)表示所述第二类组分f在位置x、时刻t的 浓度，F_f(x,t)＝-ψ_f(x)∑G[ψ_f(x+e_i)-ψ_g(x+e_i)]e_i，ψ_f(x)＝exp(-1/ρ_f(x,t))， ψ_f(x+e_i)＝exp(-1/ρ_f(x+e_i,t))；

依据公式(7)和(8)分别获得所述第一类组分g在位置x、时刻t的平 衡态的分布函数和所述第二类组分f在位置x、时刻t的平衡态的分布 函数f_i^eq(x,t)：

$g_i^{eq}(x,t)=w_i{\rho }_g(x,t)[1+\frac{e_i*u_g^{eq}(x,t)}{c_s^2}+\frac{{(e_i*u_g^{eq}(x,t\left)\right)}^2}{2c_s^4}-\frac{u_g^{eq}(x,t)*u_g^{eq}(x,t)}{2c_s^2}]---\left(7\right)$

$f_i^{eq}(x,t)=w_i{\rho }_f(x,t)[1+\frac{e_i*u_f^{eq}(x,t)}{c_s^2}+\frac{{(e_i*u_f^{eq}(x,t\left)\right)}^2}{2c_s^4}-\frac{u_f^{eq}(x,t)*u_f^{eq}(x,t)}{2c_s^2}]---\left(8\right)$

依据公式(9)和(10)分别获得所述第一类组分g在位置(x+c△t)、时刻 (t+△t)的分布函数g_i(x+c△t,t+△t)和所述第二类组分f在位置(x+c△t)、时刻 (t+△t)的分布函数f_i(x+c△t,t+△t)：

$g_i(x+c\Delta t,t+\Delta t)=-\frac{1}{\tau }[g_i^{eq}-g_i^{eq}(x,t)]+g_i^{eq}---\left(9\right)$

$f_i(x+c\Delta t,t+\Delta t)=-\frac{1}{\tau }[f_i^{eq}-f_i^{eq}(x,t)]+f_i^{eq}---\left(10\right)$

其中，$c^2=c_s^2/3;$

依据公式(11)和(12)分别获得所述第一类组分g的浓度ρ_g和所述第 二类组分f的浓度ρ_f：

$\underset{i=1}{\overset{N}{\Sigma }}{g}_i(x+c\Delta t,t+\Delta t)={\rho }_g---\left(11\right)$

$\underset{i=1}{\overset{N}{\Sigma }}f(x+c\Delta t,t+\Delta t)={\rho }_f---\left(12\right)$

获得所述磁化矢量的模值，所述磁化矢量的模值为所述第一类组分g的 浓度ρ_g。

建模分析模块33还用于依据公式(13)和(14)分别计算所述流体在表 面弛豫作用下，所述第一类组分g的浓度ρ_g随时间的变化率s_g和所述第二类 组分f的浓度ρ_f随时间的变化率s_f，并依据公式(15)和(16)分别获得表 面弛豫作用后的所述第一类组分g的浓度ρ′_g和所述第二类组分f的浓度ρ′_f：

$\frac{\partial {\rho }_g}{\partial t}{|}_{surface}=-s_g{\rho }_g---\left(13\right)$

$\frac{\partial {\rho }_f}{\partial t}{|}_{surface}=-s_f{\rho }_f---\left(14\right)$

ρ′_g＝ρ_g-s_gρ_g(15)

ρ′_f＝ρ_f-s_fρ_f(16)

依据公式(17)和(18)分别计算所述流体在体弛豫作用下，所述第一 类组分g的浓度ρ_g随时间的变化率b_g和所述第二类组分f的浓度ρ_f随时间的 变化率b_f，并依据公式(19)和(20)分别获得体弛豫作用后的所述第一类 组分g的浓度ρ″_g和所述第二类组分f的浓度ρ″_f：

$\frac{\partial {\rho }_g}{\partial t}{|}_{body}=-b_g{\rho }_g---\left(17\right)$

$\frac{\partial {\rho }_f}{\partial t}{|}_{body}=-b_f{\rho }_f---\left(18\right)$

ρ″_g＝ρ′_g-b_gρ′_g(19)

ρ″_f＝ρ′_f-b_fρ′_f(20)

依据公式(21)和(22)分别计算所述流体的扩散系数D_g和D_f，并依据 公式(23)和(24)分别获得流体扩散后的所述第一类组分g的浓度ρ″′_g和所 述第二类组分f的浓度ρ″′_f：

$\frac{\partial {\rho }_g}{\partial t}=D_g{▿}^2{\rho }_g---\left(21\right)$

$\frac{\partial {\rho }_f}{\partial t}=D_f{▿}^2{\rho }_f---\left(22\right)$

ρ″′_g＝ρ″_g-D_gρ″_g(23)

ρ″′_f＝ρ″_f-D_fρ″_f(24)

将ρ″′_g作为公式(1)中的ρ_g、ρ″′_f作为公式(2)中的ρ_f，将作为公 式(1)中的将作为公式(2)中的循环执行公式(1)-(24)。

建模分析模块33还用于在所述流体上施加压力梯度，以使所述流体的流 速发生变化；根据所述流体的组分浓度获得不同流体流速下所述磁化矢量的 模值。

三维数据体获取模块32具体用于对所述数字岩心进行三维重建获得三 维数据体；对所述三维数据体进行滤波处理，依据阈值对所述三维数据体进 行二值化处理区分出所述样品岩石中的孔隙和岩石骨架对应的三维数据体。

基于LBM模拟低场核磁共振分析流体的装置30还包括校正模块34，校 正模块34用于利用氦孔隙度测量仪测量所述样品岩石的测量孔隙度；依据所 述三维数据体计算所述样品岩石的理论孔隙度；若所述理论孔隙度与所述测 量孔隙度不相等，则调节所述阈值，再次执行所述依据阈值对所述三维数据 体进行二值化处理区分出所述样品岩石中的孔隙和岩石骨架对应的三维数据 体的步骤，直至所述理论孔隙度与所述测量孔隙度相等。

本发明实施例提供的基于LBM模拟低场核磁共振分析流体的装置可以具 体用于执行上述图1所提供的方法实施例，具体功能此处不再赘述。

本发明实施例通过循环执行单相多组分模型模拟孔隙中流体的低场核磁 共振过程，获得随时间变化的磁化矢量的模值，以及通过流体上施加压力梯 度，分析不同流体流速下磁化矢量的模值，提高了模拟低场核磁共振分析流 体分子的能力，进一步提高了低场核磁共振分析流体分子的理论基础；另外， 通过对二值化处理所用的阈值进行校正，提高了孔隙度的测量精度。

综上所述，本发明实施例通过单相多组分模型模拟样品岩石孔隙中流体 的低场核磁共振过程，分析所述流体的组分，通过所述流体的组分浓度获得 磁化矢量的模值，实现了模拟低场核磁共振分析流体分子的技术，为研究低 场核磁共振分析流体分子提供了理论基础；通过循环执行单相多组分模型模 拟孔隙中流体的低场核磁共振过程，获得随时间变化的磁化矢量的模值，以 及通过流体上施加压力梯度，分析不同流体流速下磁化矢量的模值，提高了 模拟低场核磁共振分析流体分子的能力，进一步提高了低场核磁共振分析流 体分子的理论基础；另外，通过对二值化处理所用的阈值进行校正，提高了 孔隙度的测量精度。

在本发明所提供的几个实施例中，应该理解到，所揭露的装置和方法， 可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的， 例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外 的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或 一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或 直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连 接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作 为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方， 或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或 者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中， 也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单 元中。上述集成的单元既可以采用硬件的形式实现，也可以采用硬件加软件 功能单元的形式实现。

上述以软件功能单元的形式实现的集成的单元，可以存储在一个计算机 可读取存储介质中。上述软件功能单元存储在一个存储介质中，包括若干指 令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等) 或处理器(processor)执行本发明各个实施例所述方法的部分步骤。而前述 的存储介质包括：U盘、移动硬盘、只读存储器(Read-OnlyMemory，ROM)、 随机存取存储器(RandomAccessMemory，RAM)、磁碟或者光盘等各种可以 存储程序代码的介质。

本领域技术人员可以清楚地了解到，为描述的方便和简洁，仅以上述各 功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分 配由不同的功能模块完成，即将装置的内部结构划分成不同的功能模块，以 完成以上描述的全部或者部分功能。上述描述的装置的具体工作过程，可以 参考前述方法实施例中的对应过程，在此不再赘述。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对 其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通 技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改， 或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并 不使相应技术方案的本质脱离本发明各实施例技术方案的范围。