全文快速搜索:   高级搜索

  中国石油大学学报(自然科学版)  2017, Vol. 41 Issue (1): 60-68  DOI:10.3969/j.issn.1673-5005.2017.01.007
0

引用本文 [复制中英文]

闫伟超, 孙建孟, 崔利凯, 等. 一种估算砂岩核磁T2谱新方法[J]. 中国石油大学学报(自然科学版), 2017, 41(1): 60-68. DOI: 10.3969/j.issn.1673-5005.2017.01.007.
[复制中文]
YAN Weichao, SUN Jianmeng, CUI Likai, et al. A new method of estimating NMR T2 spectrum of sandstone[J]. Journal of China University of Petroleum (Edition of Natural Science), 2017, 41(1): 60-68. DOI: 10.3969/j.issn.1673-5005.2017.01.007.
[复制英文]

基金项目

国家自然科学基金项目(41374124, 41574122);国家科技重大专项(2011ZX05014-004-002H)

作者简介

闫伟超(1991-), 男, 博士研究生, 研究方向为数字岩石物理实验、核磁共振实验及测井应用。E-mail: yanweichaoqz@163.com

文章历史

收稿日期:2015-12-09
一种估算砂岩核磁T2谱新方法
闫伟超1 , 孙建孟1 , 崔利凯1 , 闫国亮2     
1. 中国石油大学地球科学与技术学院, 山东青岛 266580;
2. 中国石油勘探开发研究院西北分院, 甘肃兰州 730020
摘要: 根据核磁共振弛豫机制, 基于三维数字岩心技术建立一种估算砂岩核磁共振横向弛豫(T2)谱的新方法。采用过程模拟法重建砂岩三维数字岩心, 通过Delaunay剖分将砂岩微观孔隙结构划分为多个四面体孔隙单元, 再用几何方法和Monte Carlo方法相结合求出每个孔隙单元的表面积与体积之比(S/V), 进而获得砂岩模拟的核磁T2谱。经过大量模拟试验筛选参数, 减小计算误差。考虑单峰态和双峰态核磁T2谱, 模拟结果和实验结果基本吻合。该方法有助于研究粒度分布、孔隙结构等对岩石核磁共振响应的影响。
关键词: 核磁T2    数字岩心    过程模拟法    Delaunay剖分    Monte Carlo方法    
A new method of estimating NMR T2 spectrum of sandstone
YAN Weichao1 , SUN Jianmeng1 , CUI Likai1 , YAN Guoliang2     
1. School of Geosciences in China University of Petroleum, Qingdao 266580, China;
2. Research Institute of Petroleum Exploration & Development-Northwest, PetroChina, Lanzhou 730020, China
Abstract: Based on nuclear magnetic resonance(NMR)relaxation mechanism and three-dimensional digital core technology, a new method of estimating the spectrum of sandstone transverse relaxation time (T2) is established. A process-based method was used to reconstruct 3D digital rock models of sandstone. Sandstone microstructures were divided into tetrahedral pore units using Delaunay tessellation. Considering the geometry and using a Monte Carlo approach, the surface-to-volume ratio (S/V) of each pore unit was calculated, and then the T2 spectrum of sandstone could be simulated. After a large number of simulations, suitable parameters were selected in order to reduce calculation errors. Unimodal and bimodal T2 spectrum were considered; both results agreed well with laboratory tests. This new method is helpful in studies of different NMR response effects with varying grain sizes or pore structures.
Keywords: NMR T2 spectrum    digital core    process-based method    Delaunay tessellation    Monte Carlo method    

核磁共振(NMR)实验通过测量岩石孔隙流体的核磁弛豫特性提供有关岩石孔隙度及其分布、渗透率、流体性质及含量等信息[1-2], 已经成为无损检测岩石物性和含油性特征的重要手段。岩石骨架表面、孔隙中流体的物理特征以及孔隙内部微观结构决定了岩石的宏观物理性质, 因此从微观尺度进行岩石的核磁共振响应特性的模拟对于核磁共振岩石物理基础研究具有重要的价值。目前, 核磁共振T2谱的微观尺度模拟主要采用随机行走算法[3-4], 它通过模拟粒子的随机行走, 得到磁化强度衰减曲线, 并通过解谱方法得到核磁共振T2谱。这种算法与离散体素分辨率和统计粒子数量有关, 分辨率越高、统计粒子数量越多, 模拟结果越接近实验结果, 但同时会延长计算时间。笔者根据核磁共振弛豫机制, 基于三维数字岩心技术提出一种估算砂岩核磁T2谱的新模拟方法, 并通过单峰态和双峰态核磁T2谱对比进行验证。

1 基础理论

对于核磁共振实验, 在测量岩样时其横向弛豫过程受到自由弛豫、表面弛豫、扩散弛豫3种机制的作用, 即

$ \frac{1}{{{T_2}}} = \frac{1}{{{T_{2{\rm{B}}}}}} + \frac{1}{{{T_{{\rm{2S}}}}}} + \frac{{D{{(\gamma G{T_{\rm{E}}})}^2}}}{{12}}. $ (1)

式中, T2为横向弛豫时间; T2B为流体的自由弛豫时间; T2S为表面弛豫时间; γ为核旋磁比; G为磁场梯度; TE为回波间隔; D为流体的扩散系数。

实验测量表明, T2B>3 000 ms, 自由弛豫比表面弛豫慢很多, 因此式(1)右边第一项可以忽略; 在没有磁场梯度或(GTE)值很小时, 式(1)右边第三项也可以忽略, 此时式(1)简化为

$ \frac{1}{{{T_2}}} \approx \frac{1}{{{T_{2{\rm{S}}}}}} = {\rho _2}\left( {\frac{S}{V}} \right). $ (2)

式中, ρ2为岩石的横向表面弛豫强度; S为孔隙的表面积; V为孔隙的体积。

由公式(2)可以看出, 对于饱和水的岩心, 如果已知ρ2值, 求出每个孔隙的表面积与体积比(S/V), 就能得到核磁T2谱。该方法将弛豫时间分布与孔径大小关联, 对两个条件进行假设:一是孔隙中的水扩散在快扩散域中, 二是孔隙的表面弛豫或整个岩心中顺磁离子在每个孔隙表面平均密度相差不大。由于扩散距离短, 自旋游动范围多数局限于一个孔隙, 认为每个孔隙或一组相邻孔隙的弛豫相互独立, 不考虑孔隙之间的耦合。通过构建岩心三维空间模型, 使用Delaunay剖分方法可以确定不同类型的孔隙, 几何方法和Monte Carlo方法结合则可以求出每个孔隙的S/V, 进而获得核磁T2谱, 这是估算砂岩核磁T2谱新方法的主要流程。

2 核磁T2谱估算新方法 2.1 砂岩岩心重建

采用过程模拟法构建砂岩岩心, 模拟岩石的形成过程, 可以重现真实砂岩的孔隙结构和连通特性[5-8]。该方法结合岩石的颗粒粒度分布, 通过对沉积类岩石形成过程(包括沉积、压实和成岩作用)的模拟建立数字岩心。

首先, 通过实验室方法(筛析法、沉降法、激光衍射法、薄片分析法等)或者数学模型获得砂岩的颗粒粒径累积分布, 将颗粒等效为不同直径的球体, 采用颗粒堆积算法进行随机沉积。假设颗粒在重力的作用下沉积, 判断该颗粒下落时所有可能遇到的其他颗粒, 最终将落点最高的颗粒作为落点颗粒。如果在下落过程中接触到其他已经沉积的颗粒, 则旋转一定角度继续下落, 直到形成3个颗粒稳定的接触关系, 记录此时颗粒中心点的坐标。执行下一个颗粒沉积, 重复该过程直到颗粒填满整个指定空间。

初始颗粒堆积为砂岩压实和成岩模拟奠定了基础。根据计算公式, 将每一个颗粒的垂向坐标向下移动指定的长度, 以此模拟砂岩的压实过程。计算公式为

$ z = {z_0}(1-\lambda ). $ (3)

式中, z为下移后新的垂向坐标; z0为原始垂向坐标; λ为压实因子, 且0≤λ≤1。

成岩算法模拟是在压实算法模拟结果基础上进行。以压实算法模拟结果作为输入数据, 采用各向同性增长算法模拟砂岩成岩过程, 最终得到与实验测得的孔隙度相近的三维数字岩心模型。与常规过程模拟法构建数字岩心不同, 本研究不仅须输出“0”、“1”格式的三维数字岩心数据体用于目标孔隙度控制, 而且须输出每个颗粒的空间位置及相应的颗粒半径, 以便于后续的四面体剖分。

2.2 Delaunay剖分

由公式(2)可以看出, 为了计算每个孔隙的S/V, 必须先把模型划分为多个孔隙单元。此时, 可以利用三维Delaunay剖分方法将每个孔隙单元等效为独立不重叠的四面体单元[9-12]

令平面域(R2)上M个离散点的集合为P={P1, P2, …, PM}, 有多种方法能够实现点集P的三角剖分(剖分的所有边界区域都是三角形), 其中俄国科学家Delaunay证明:必定存在且仅存在一种三角剖分算法, 使得所有三角形的最小内角之和最大。这种剖分方法也被称作Delaunay剖分, 它能够尽可能地避免病态三角形的出现。同样, 在三维空间, 令Q={Q1, Q2, …, QN}为空间域(R3)上N个离散点的集合, 可以实现点集Q的四面体剖分, 将空间离散为独立的四面体单元。

构建空间Delaunay四面体剖分, 首先要确定空间离散点数据的三维坐标。经过过程模拟法重建砂岩岩心后, 可以得到每个球形颗粒的球心坐标, 将球心坐标集合视为离散点集合。将离散点集合预处理后, 从第一个点P1找到离它最近的点P2, 再找P3使得P1P2P3满足平面Delaunay准则(每个三角形的外接圆弧内不包含其他点), 建立第一个Delaunay三角形。从一个三角面出发, 寻找一个点, 再使组成的四面体满足三维Delaunay准则(每个四面体的外接球均不包含除此四面体的顶点以外的网格点)。为了保证四面体剖分的完整性和正确性, 先找到第一个四面体, 然后由它的各个面为底面建立新的四面体, 这样逐层建网, 直到将空间所有离散点进行四面体剖分。每个四面体单元的4个顶点为球形颗粒的球心, 颗粒视为骨架, 颗粒不重叠的空间视为孔隙, 如图 1

图 1 四面体单元 Fig.1 Tetrahedral unit
2.3 每个孔隙的S/V计算

随机行走法模拟核磁T2谱过程中“随机粒子”碰撞岩石骨架表面, 因此计算结果受到骨架表面积的影响, 而且随机行走法是对离散化后的岩心进行模拟, 骨架表面积与离散程度直接相关, 需要高度离散化才能近似统计颗粒骨架表面积, 但这样使数字岩心数据体非常庞大。本文中结合几何方法和Monte Carlo方法求解每个孔隙的S/V, 克服了离散化后统计孔隙表面积和体积的缺陷。

经过Delaunay剖分, 得到一系列的四面体单元, 每个四面体单元包含一个孔隙。已得到每个颗粒球体的球心坐标和半径, 单纯用几何方法可以快速计算出大部分孔隙的表面积和体积。

基本思想:用四面体单元体积减去4个球扇形(球形颗粒与四面体相交部分)体积, 可以得到孔隙体积。4个球面三角形(球形颗粒被四面体截下的表面)表面积之和为孔隙表面积。如果球扇形有重叠, 去掉重叠部分体积和表面积。

主要公式:假设四面体4个顶点坐标分别为(x1, y1, z1)、(x2, y2, z2)、(x3, y3, z3)、(x4, y4, z4), 半径分别为r1r2r3r4

$ {V_{\rm{t}}} = \frac{1}{6}\det \left| {\begin{array}{*{20}{c}} 1&1&1&{{x_1}}\\ {{y_1}}&{{y_2}}&{{y_3}}&{{y_4}}\\ {{z_1}}&{{z_2}}&{{z_3}}&{{z_4}} \end{array}} \right|, $ (4)
$ {S_i} = ({\omega _i}-{\rm{ \mathit{ π} }}){r_i}^2, $ (5)
$ {V_i} = \frac{1}{3}{S_i}{r_i}. $ (6)

式中, Vt为四面体体积; Si为第i个球面三角形的表面积; ωi为第i个球面三角形内角和(相当于3个四面体两面所夹二面角之和); Vi为第i个球扇形的体积。

$ \theta = {\rm{ \mathit{ π} }}-\arccos \frac{{{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over n} }_1} \cdot {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over n} }_2}}}{{\left| {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over n} }_1}} \right| \cdot \left| {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over n} }_2}} \right|}}. $ (7)

式中, θ为四面体两面所夹二面角; ${{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over n} }_1}}$${{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over n} }_2}}$分别为两面的法向量。

重叠区域的二维投影如图 2所示。

图 2 重叠区域二维投影 Fig.2 2D projection of overlapping area
$ {V_{\rm{o}}} = \frac{{\rm{ \mathit{ π} }}}{{12d}}{\left( {R + r-d} \right)^2}({d^2} + 2dr-3{r^2} + 2dR + 6rR-3{R^2}), $ (8)
$ {h_1} = R-\frac{{{R^2}-{r^2} + {d^2}}}{{2d}}, $ (9)
$ {h_2} = r-\frac{{{r^2}-{R^2} + {d^2}}}{{2d}}, $ (10)
$ {S_{{\rm{o}}1}} = 2{\rm{ \mathit{ π} }}R{h_1}, $ (11)
$ {S_{{\rm{o}}2}} = 2{\rm{ \mathit{ π} }}r{h_2}. $ (12)

式中, Vo为重叠区域的体积; So1为重叠区域的第一部分表面积; So2为重叠区域的第二部分表面积; d为两球心坐标距离; Rr分别为重叠的两球半径; h1h2为重叠球冠的高。

该几何方法可以快速求出孔隙的S/V, 但有一小部分结果出现了负值。出现这种现象的原因是球体半径过大, 球扇形有部分超出了四面体单元范围。对于此类问题, 本文中采用Monte Carlo方法近似求得孔隙的表面积和体积。

Monte Carlo方法亦称随机模拟(random simulation)方法[13], 首先建立一个概率模型或随机过程, 使它的参数等于问题的解; 然后通过对模型或过程的观察或抽样试验计算所求参数的统计特征, 最终给出所求解得近似值。

根据四面体坐标构建立方体, 在立方体中产生N个随机数。判断随机数是否在孔隙中, 即根据四面体顶点坐标和有向体积概念判断随机数是否在四面体内部, 且在球体外部。统计孔隙中随机数个数K, 孔隙体积为

$ {V_{\rm{p}}} \approx \frac{K}{N}{V_1}. $ (13)

式中, Vp为孔隙体积; Vl为立方体体积。

同样, 统计孔隙表面薄层的随机数个数, 可得孔隙的表面积。

Monte Carlo方法不仅可以近似求出孔隙的体积和表面积, 而且不必单独考虑球体是否重叠。随着产生的随机数N越大, 模拟结果更接近准确值, 但计算量也会增大, 因此为了提高计算精度同时减少计算量, 它仅用于几何方法产生负值的情况。用几何方法和Monte Carlo方法相结合, 最终可以求出每个孔隙单元的S/V

2.4 获得T2

根据公式(2), 现已求出每个孔隙的S/V, 若已知ρ2值, 便可得到每个孔隙对应的T2值。“表面弛豫强度”ρ2的值与岩石颗粒表面及胶结物的性质有关, 同时也受到温度的影响[14]。有多种方法可以获得该参数, 包括孔隙图像数值分析、脉冲场梯度核磁共振、压汞毛管压力曲线等, 本文中采用压汞毛管压力曲线方法求解ρ2[15]

每个孔隙对应一个T2值, 同时可以算出每个孔隙占整块岩心的体积分数。以T2弛豫时间为横坐标, 区间孔隙度为纵坐标, 便可得到模拟核磁T2谱。类似常规核磁T2谱反演处理, 上述利用每个孔隙单元的S/V得到的核磁T2谱也需要进行布点, 本文中采用对数布点方式, 并经过三次样条插值函数平滑得到最终的核磁T2[16]

对于Delaunay剖分得到的独立的孔隙, 会有极个别四面体由于较长较窄, 导致S/V偏小, T2值过大, 本文中将其按充填胶结物处理。这样会使模拟得到的岩心孔隙度略小于真实孔隙度, 因此在过程模拟法建立砂岩岩心时可以取孔隙度略大于真实孔隙度, 以便最终模拟得到的孔隙度接近真实值。与随机行走算法不同, 新方法不必处理磁化强度衰减曲线数据, 而是直接得到核磁T2谱。这种模拟方法不仅可以划分不同类型的孔隙, 也可以确定每个孔隙的表面积和体积。

3 方法验证与模拟分析 3.1 Monte Carlo方法准确性验证

在本文砂岩微观结构核磁T2谱模拟方法中, 几何方法求孔隙的表面积和体积完全是解析解, 但Monte Carlo方法是统计模拟方法, 求解结果会有误差, 因此须统计相对误差验证方法的准确性并筛选参数。

构建边长为2r μm的立方体, 完全包含了边长为r μm的球体, 以球心坐标为原点。利用Monte Carlo方法在球体表面内侧取长度为10-a μm的厚度, 与球体表面组成球环。在立方体中生成随机数, 作为随机点的坐标(x, y, z), 一共生成N个随机点。对于每个随机点的三维坐标值, 如果满足

$ \left| {\sqrt {{x^2} + {y^2} + {z^2}}-r} \right| < {10^{-a}}, $ (14)

则说明此点位于球环内部。

Monte Carlo方法的模拟过程引入3个变量参数:球体边长r、球环厚度指数a和随机点数N, 采用控制变量的方法研究不同参数对模拟实验结果的影响, 如表 1所示。试验编号1-7为固定rN, 研究不同球环厚度指数a的试验统计相对误差平均值变化; 试验编号8-13为固定ra, 研究不同随机点数N的试验统计相对误差平均值变化; 试验编号14-20为固定Na, 研究不同球体边长r的试验统计相对误差平均值变化。

表 1 Monte Carlo方法相对误差平均值变化 Table 1 Variation of mean value of relative error by Carlo Monte method

通过多次变N、r、a三参数的模拟试验, 发现Monte Carlo方法在不同参数条件下得到结果的相对误差具有一定规律性。随着模拟随机点N增大, 相对误差越小, 因此在计算条件允许时尽可能增大N值。对于砂岩粒度中值r, 当0 < r < 100 μm时, a取-1相对误差较小; 当100 < r < 1 000 μm时, a取0相对误差较小。

3.2 单峰态核磁T2谱模拟

由于不同类型砂岩孔隙结构有较大的差异, 在核磁T2谱上一般反映为单峰态和双峰态出现[17]。为了验证本文中核磁T2谱模拟方法对于多种类型砂岩的适用性, 分别讨论了单峰态核磁T2谱和双峰态核磁T2谱两种类型。

对于单峰态核磁T2谱, 本文中分别采用人造砂岩和实际砂岩样品, 结合常规物性分析、粒度分析、低磁场核磁共振实验、高压压汞实验等资料进行核磁T2谱新方法验证。人造砂岩由74 μm单粒径的石英砂混合环氧树脂制备, 实际砂岩样品选取某地区两块不同核磁T2谱单峰形态的砂岩样品, 分别标记为A1、A2、A3。表 2为3块岩心的基本物性参数。

表 2 单峰态核磁T2谱模拟三块砂岩岩心的基本物性参数 Table 2 Basic petrophysical parameters of three sandstones of NMR unimodal T2 spetrum simulation

本文中采用过程模拟法建立3组模拟岩心, 数字岩心A1沉积颗粒粒径均为75 μm, 数字岩心A2、A3沉积颗粒粒径根据切割的部分实际岩心粒度分析结果获得, 如图 3所示。

图 3 单峰态核磁T2谱模拟岩心粒度分析结果 Fig.3 Grain size analysis results of NMR unimodal T2 spetrum simulation

首先将粒度频率分布转为粒度累积分布, 为了保证构建的数字岩心稳定, 截掉累积分布曲线两端5%的粒径, 对其余的粒径随机选取。构建模型孔隙度分别为31.35%、14.18%、8.20%。图 4为用ImageJ软件生成的三维数字岩心图像, 黄色代表砂岩颗粒, 透明空间为孔隙。

图 4 重建砂岩岩心图 Fig.4 Reconstruction of sandstone core pictures

由于缺少相同深度的岩心, 因此须进行高压压汞测试, 如图 5所示。将压汞曲线与核磁共振测量T2谱曲线重叠, 求得岩心A1、A2和A3的表面弛豫率ρ2值分别为12.45、10.34、21.35 μm/s。经过Delaunay剖分将三维数字岩心模型每个孔隙单元等效为独立的不重叠的四面体单元, 结合几何方法与Monte Carlo方法求解每个孔隙的S/V, 调整表面弛豫率使模拟结果峰值与实际测量结果峰值点对应的弛豫时间重合, 最终得到模拟T2谱, 如图 6所示。

图 5 单峰态核磁T2谱模拟岩心压汞分析结果 Fig.5 Mercury penetration analysis results of NMR uninodal T2 spetrum simulation
图 6 单峰态核磁T2谱数值模拟与实验测量结果对比 Fig.6 Comparison of NMR unimodal T2 spectrum simulation results and measurement results

图 6数值模拟核磁T2谱与实验测量结果对比, 可以发现利用新方法同样能够很好地模拟单峰态核磁共振T2谱, 无论T2谱形态还是峰值对应的T2弛豫时间, 模拟结果与实验结果基本吻合。说明新方法在微观尺度进行砂岩核磁共振响应特性的模拟有很好的效果。同时, 该方法也能够输出每个孔隙的空间几何特征数据(位置、表面积和体积), 便于对不同孔隙进行研究。

3.3 双峰态核磁T2谱模拟

对于双峰态核磁T2谱, 本文中结合砂岩常规物性分析、粒度分析、低磁场核磁共振实验、高压压汞实验等资料进行核磁T2谱新方法验证。实验选取某地区3块具有不同核磁T2谱双峰形态的砂岩样品, 分别标记为B1、B2、B3, 由于粒度分析与高压压汞测试对岩心具有破坏性, 因此这两类实验选取相同深度不同位置的岩心得到样品相近的粒度分布和孔隙结构特征。表 3为砂岩岩心的基本物性参数, 3块岩心均为低渗透砂岩。

表 3 双峰态核磁T2谱模拟砂岩岩心的基本物性参数 Table 3 Basic petrophysical parameters of NMR biomodal T2 spectrum simulation

图 7图 8分别为3块岩心的低磁场核磁共振分析结果和粒度分析结果, 从图中可以看出, 核磁T2谱与颗粒粒径分布均呈现双峰态, 认为颗粒分布不均匀影响了岩石的孔隙结构, 导致了核磁T2谱双峰态的出现。邓克俊等[18-19]也认为岩石颗粒粒度分布与核磁T2谱形态有明显的相关关系。同时, 岩石越致密, 大颗粒粒径和长T2弛豫时间频率越低, 小颗粒粒径和短T2弛豫时间频率越高。

图 7 双峰态核磁T2谱模拟岩心核磁共振分析结果 Fig.7 NMR analysis results of biomodal T2 spectrum simulation
图 8 双峰态核磁T2谱模拟岩心粒度分析结果 Fig.8 Grain size analysis results of NMR biomodal T2 spectrum simulation

图 9为3块岩心的高压压汞结果, 根据进汞压力与孔径的关系, 除了相对较低的压力点出现峰值代表较大孔隙外, 在高压端汞饱和度频数均出现抬升, 说明小尺寸孔隙也占有一定的比例, 这与核磁T2谱双峰态反映的孔隙结构关系基本一致, 压汞高压端汞饱和度含量较小是由于低渗透岩心中部分小孔隙难以充注导致。

图 9 双峰态核磁T2谱模拟岩心压汞分析结果 Fig.9 Mercury penetration analysis results of biomodal T2 spectrum simulation

本文中采用过程模拟法分别构建以上3块砂岩的三维数字岩心, 通过多次模拟实验, 发现沉积体整体边长设置在3 000 μm以上, 沉积颗粒数在15 000个以上, 控制目标孔隙度所需的分辨率值小于3 μm, 才能够得到相对较好的实验结果, 否则会导致核磁T2谱双峰不明显、剖分后计算总孔隙度偏小等异常。经过Delaunay剖分以及结合几何方法与Monte Carlo方法求解每个孔隙的S/V, 只要确定岩心ρ2值即可得到模拟核磁T2谱。对于岩心ρ2值的选取, 根据Kleinberg的方法, 结合核磁T2谱与压汞曲线, 得到3块砂岩的ρ2值分别为2.0、1.4、4.1 μm/s。但在核磁T2谱新方法模拟中, 由于构建的砂岩岩心骨架颗粒均为球形, 无法反映实际孔隙粗糙表面对孔隙表面积的影响, 导致核磁T2谱的整体偏移, 因此该参数的选取须进行调试。图 10为经过调试后新方法ρ2值分别选9.5、10.0、12.5 μm/s得到的3块砂岩模拟核磁T2谱以及随机行走法模拟结果。

图 10 双峰态核磁T2谱模拟结果 Fig.10 NMR biomodal T2 spectrum simulation results

从图中可以看出, 新方法可以模拟双峰态的核磁T2谱, 但相对较长的弛豫时间谱峰幅度低于实际值, 这主要是由于构建砂岩岩心时沉积的颗粒空间位置与实际不符导致, 而随机行走法对于双峰态的核磁T2谱表征并不明显。对于单峰态砂岩, 由于孔隙结构单一, 沉积颗粒的空间位置对结果影响小; 而对于双峰态砂岩, 尽管颗粒粒度分布与实际一致, 但沉积颗粒空间位置不同造成的孔隙结构非均质性程度差异比较明显(如沉积颗粒按照粒径从小到大顺序进行沉积会导致相对较长的弛豫时间谱峰幅度远高于实际值)。除非采用X射线CT扫描等技术获得不同粒径的颗粒对应的空间位置, 常规的岩石物理分析手段很难得到不同颗粒的真实三维空间分布。采用孔隙网络模型计算的岩心渗透率进行约束[20], 对比实验渗透率值筛选不同颗粒三维空间分布砂岩岩心模型, 使模型的物理属性更加准确, 进而得到相对较好的核磁T2谱模拟结果, 但由于多解性的存在, 在缺乏实验核磁T2谱的条件下如何筛选与实际岩心一致的模型依然困难。

3.4 不同粒度分布对核磁T2谱影响

为了单因素研究粒度分布对核磁T2谱的影响, 需要尽可能地保持其他影响因素一致。在一定的沉积环境下, 天然砂岩的粒度分布满足确定的数学模型, 如威布尔分布、对数正态分布等[21]。本文中采用单峰威布尔分布构建砂岩粒度分布, 通过调整威布尔分布参数, 选择3组粒度分布数据, 使其中值粒径分别为200、300、400 μm。粒度分布数据作为颗粒堆积的初始信息输入过程模拟法程序中, 构建3种不同粒径的三维数字岩心, 孔隙度均为20%, 假设模拟中岩心ρ2值均为4.0×10-3cm/s, 岩心完全含水。经过四面体剖分与S/V的计算, 得到3种不同中值粒径、不同粒度分布的三维数字岩心核磁T2谱, 如图 11所示。

图 11 不同粒度分布核磁T2 Fig.11 NMR T2 spectrums of different grain size distributions

从图中可以看出, 对于不同粒度分布的岩心, 随着粒度中值增大, 核磁T2谱形态没有明显的变化, 但峰值对应的T2弛豫时间有增大的趋势。由于完全含水条件下核磁T2谱基本反映了岩石孔隙结构, 因此该结果也说明了当孔隙度相同时, 岩心粒度不同会导致孔隙结构有较大差异, 而且随着粒度中值增大, 孔喉中值半径也增大。

4 结论

(1) 根据核磁共振弛豫机制, 基于三维数字岩心技术建立一种新的核磁共振响应特性模拟方法。利用过程模拟法重建微观砂岩岩心, 通过Delaunay剖分将砂岩划分为多个四面体孔隙单元, 结合几何方法和Monte Carlo方法求出每个孔隙单元的S/V, 进而获得砂岩模拟的核磁T2谱。该方法易于实现, 计算结果较为准确, 与随机行走方法相比, 对于双峰态核磁T2谱表征更明显。

(2) Monte Carlo方法作为统计模拟方法, 求解结果会出现误差。在模拟时基于实际计算条件尽可能增大模拟随机点数, 并根据颗粒粒度中值选择合适的模拟参数。

(3) 对单峰态和双峰态的核磁T2谱模拟分析, 认为单峰态的核磁T2谱模拟时限制条件较少, 计算结果更加准确; 而双峰态核磁T2谱模拟时, 虽然容易得到双峰态的核磁T2谱结果, 但由于孔隙结构的复杂性, 构建砂岩模型颗粒空间分布无法准确地确定, 导致峰值幅度与实际值有一定的差异。

(4) 尽管本研究数值模拟方法不能完全代表真实砂岩岩心, 但对岩石物理基础研究具有重要的价值, 有助于研究粒度分布、孔隙结构等对岩石核磁共振响应的影响。当确定孔隙空间流体分布后, 可以进一步模拟不同含水饱和度条件下核磁共振响应, 对岩石的核磁共振响应特征进行更深入的研究。

参考文献
[1]
DUNN K J, BERGMAN D J, LATORRACA G A. Nuclear magnetic resonance: petrophysical and logging applications[M]. Pergamon: Elsevier Science, 2002.
[2]
COATES G R, XIAO L, PRAMMER M G. NMR logging: principles and applications[M]. Texas: Gulf Professional Publishing, 1999.
[3]
孙建孟, 姜黎明, 刘学锋, 等. 数字岩心技术测井应用与展望[J]. 测井技术, 2012, 36(1): 1-7.
SUN Jianmeng, JIANG Liming, LIU Xuefeng, et al. Log application and prospect of digital core technology[J]. Well Logging Technology, 2012, 36(1): 1-7.
[4]
TALABI O, ALSAYARI S, IGLAUER S, et al. Pore-scale simulation of NMR response[J]. Journal of Petroleum Science & Engineering, 2009, 67(3): 168-178.
[5]
∅RENP E, BAKKES. Reconstruction of Berea sandstone and pore-scale modelling of wettability effects[J]. Journal of Petroleum Science and Engineering, 2003, 39(3/4): 177-199.
[6]
∅RENP E, BAKKES. Process based reconstruction of sandstones and predictions of transport properties[J]. Transport in Porous Media, 2002, 46(2/3): 311-343.
[7]
BAKKE S, ∅REN P E. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks[J]. SPE Journal, 1997, 2(2): 136-149. DOI:10.2118/35479-PA
[8]
姚军, 赵秀才, 衣艳静, 等. 数字岩心技术现状及展望[J]. 油气地质与采收率, 2005, 12: 52-54.
YAO Jun, ZHAO Xiucai, YI Yanjing. The current situation and prospect on digital core technology[J]. Petroleum Geology and Recovery Efficiency, 2005, 12: 52-54.
[9]
王建华, 徐强勋, 张锐. 任意形状三维物体的Delaunay网络生成算法[J]. 岩石力学与工程学报, 2003, 22(51): 717-722.
WANG Jianhua, XU Qiangxun, ZHANG Rui. Delaunay algorithm and related procedure to generate the tetrahedron mesh for an object with arbitrary boundary[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(51): 717-722.
[10]
余杰, 吕品, 郑昌文. Delaunay三角网构建方法比较研究[J]. 中国图像图形学报, 2010, 15(8): 1158-1166.
YU Jie, LV Pin, ZHENG Changwen. A comparative research on methods of Delaunay triangulation[J]. Journal of Image and Graphics, 2010, 15(8): 1158-1166.
[11]
FANG T P, PIEGL L. Delaunay triangulation in three dimensions[J]. Computer Graphics and Applications, IEEE, 1995, 15(5): 62-69. DOI:10.1109/38.403829
[12]
周培德. 计算几何-算法设计与分析[M]. 北京: 清华大学出版社, 2006.
[13]
曲双石, 王会娟. MonteCarlo方法及其应用[J]. 统计教育, 2009(1): 45-55.
QU Shuangshi, WANG Huijuan. Monte Carlo method and its application[J]. Statistical Thinktank, 2009(1): 45-55.
[14]
谢然红, 肖立志, 傅少庆. 饱和水岩石核磁共振表面弛豫温度特性[J]. 中国石油大学学报(自然科学版), 2008, 32(2): 44-46.
XIE Ranhong, XIAO Lizhi, FU Shaoqing. Temperature effect of NMR surface relaxation in water saturated rocks[J]. Journal of China University of Petroleum (Edition of Natural Science), 2008, 32(2): 44-46.
[15]
KLEINBERG R L. Utility of NMR T2, distributions, connection with capillary pressure, clay effect, and determination of the surface relaxivity parameter ρ2[J]. Magnetic Resonance Imaging, 1996, 14(7/8): 761-767.
[16]
苏俊磊, 孙建孟, 张守伟. 核磁共振弛豫信号的多指数反演及应用[J]. 石油天然气学报, 2011, 32(6): 87-91.
SU Junlei, SUN Jianmeng, ZHANG Shouwei. Application of multi-exponential inversion of NMR relaxation signal[J]. Journal of Oil and Gas Technology, 2011, 32(6): 87-91.
[17]
何雨丹, 毛志强, 肖立志, 等. 利用核磁共振T2分布构造毛管压力曲线的新方法[J]. 吉林大学学报(地球科学版), 2005, 35(2): 177-181.
HE Yudan, MAO Zhiqiang, XIAO Lizhi, et al. A new method to obtain capillary pressure curve using NMR T2 distribution[J]. Journal of Jilin University (Earth Science Edition), 2005, 35(2): 177-181.
[18]
邓克俊, 谢然红. 核磁共振测井理论及应用[M]. 东营: 中国石油大学出版社, 2010.
[19]
孙建孟, 赵建鹏, 闫伟超, 等. 应用核磁T2谱与数字岩心技术计算粒度分布方法[J]. 中国石油大学学报(自然科学版), 2013, 37(3): 57-62.
SUN Jianmeng, ZHAO Jianpeng, YAN Weichao, et al. Calculation of grain size distribution using NMR T2 spectrum and digital rock technology[J]. Journal of China University of Petroleum (Edition of Natural Science), 2013, 37(3): 57-62.
[20]
VALVATNE P H, BLUNT M J. Predictive pore-scale modeling of two-phase flow in mixed wet media[J]. Water Resources Research, 2004, 40(7): 1-21.
[21]
KONDOLF G M, AAHIKARI A. Weibull vs. lognormal distributions for fluvial gravels[J]. Journal of Sedimentary Research, 2000, 70(3): 456-460. DOI:10.1306/2DC4091E-0E47-11D7-8643000102C1865D