致密砂岩储层测井数字岩石物理研究需求、进展与挑战

李华阳，李潮流，周灿灿，胡法龙，李霞，邹友龙

（中国石油勘探开发研究院，北京 100083）

0 引 言

据翟光明、邱中建院士等［1－2］的研究，中国陆上未来油气勘探的主要领域包括前陆盆地、断陷盆地富油气凹陷、坳陷盆地大面积地层岩性带和大型古隆起。从测井专业角度看，目前及未来油气勘探主体与以往常见油气藏的区别主要是储层的复杂化，即岩性、孔隙空间类型和流体分布特征多样化，给测井采集、解释评价与岩石物理研究带来很大挑战。上述4类勘探对象中，复杂储层主要包括复杂碎屑岩（特低孔隙度渗透率、致密砂岩、页岩等）、缝洞型碳酸盐岩、火山岩储层。其中，以鄂尔多斯盆地、四川盆地为代表的致密砂岩油气区勘探开发近2年取得重要突破，已成为目前最为现实的非常规油气资源［3］。该类储层也是本文论述的对象。

致密砂岩储层主要特征是孔喉细小、孔隙度渗透率条件极差、孔隙结构极复杂。针对常规砂岩储层条件研发的岩石物理分析技术和手段已难以满足致密砂岩储层实验精度和时效性的要求，特别是常规驱替实验不适用于此类储层［3］。以CT扫描和数值模拟为核心的测井数字岩石物理性质研究已逐渐成为解决相关难题的重要手段，并在渗透率、电性模拟等方面取得重要进展，但同时也遇到了很多挑战。

1 测井数字岩石物理研究历程与主要内容

利用数值模拟的方法分析岩石渗流等宏观物理性质的研究早在20世纪50年代就由Purcell等［4］提出，早期有毛细管束模型、二维电路网络模型等。20世纪70年代由于计算机及图像分析技术的出现使得孔隙网络模型的研究得到进一步发展，但直到20世纪90年代以后随着CT技术的发展才真正开始了模拟砂岩三维孔隙的数字岩石物理研究。

美国科学院院士Amos Nur教授1977年在斯坦福大学组建SRB（Stanford Rock Physics ＆ Borehole Geophysics Project）并致力于地层条件下的岩石物理性质研究，1995年左右开始数字岩石物理性质研究，目前在石油天然气行业处于领先地位，引领相关学科的发展方向，并于2007年创建In Grain公司进行相关技术服务和咨询。Schlumberger、Shell和Chevron等服务公司和油公司都较早地开展了数字岩石物理性质方面的研究工作。

中国石油大学孙建孟教授、陶果教授和姚军教授等在2000年以后开始利用格子气自动机、LBM及有限元等方法模拟电学、声学及渗透率等性质以及两相流模拟［5］。毕业于英国帝国理工大学的董虎博士于2010年在厦门注册成立了数岩科技公司，专门从事相关技术研发和服务。

综合国内外关于数字岩石物理研究的发展历程，测井数字岩石物理是依据物理学原理和方法，基于反映岩石内部孔隙特征的图像采用数字模拟手段分析储层岩石的物理性质和输运特性新兴学科，主要内容包括2部分。

（1）数字孔隙格架的构建。目前用于构建数字孔隙格架的方法主要有物理实验法和数值重建法2种。物理实验法利用CT或扫描电镜对岩心内部结构进行三维灰度成像，再通过图像二值化处理分离孔隙和骨架，以获取岩石真实三维孔隙格架。随着CT扫描技术和成像技术的进步，该方法的应用已经越来越普遍。国外还有学者提出采用环氧树脂灌注孔隙并通过聚焦扫描进行孔隙成像的实验方法，但只能处理岩心薄片［6］。

数值重建法基于具体的岩心薄片图像分析数据，利用统计模型或模拟沉积与成岩过程的方法建立数字孔隙格架，典型算法有高斯场法、模拟退火法、沉积过程重建法和马尔可夫随机重建法等［7－11］。需要指出的是，沉积过程重建方法采用不同尺寸球形、椭球形颗粒堆积并模拟沉积岩的沉积、压实、胶结等成岩过程，最终得到近似于真实砂岩的三维孔隙格架，与其他几种数值重建方法相比较，其结果更接近真实岩心特征［8］。

无论采用哪种方法，在获取三维数字孔隙格架的同时还必须有效区分孔隙和喉道，这是对孔隙网络的一种近似简化以保证模拟计算成功的关键。

（2）岩石物理特性的模拟。在获取三维孔隙格架的基础上利用多种算法模拟岩石特征参数（见图1）。图1表明，目前的测井数字岩石物理性质研究已经可以实现包括电、声、核磁共振及渗流等主要参数的模拟，基本涵盖了常规实验测试项目，可对致密砂岩的岩石物理特征进行较全面的分析研究。

图1 基于三维孔隙格架的岩石物理性质数值模拟对象

2 近2年研究进展概述

（1）扫描测试仪器配套能够获取各种分辨率图像。以美国Xradia公司为代表针对岩心的CT扫描设备逐渐配套完善，目前已商业化应用的包括UltraXRM－L200型、Micro XCT－400型、MicroXCT－200型和VersaXRM－500型。美国GE公司也有类似CT设备，它们的最高成像分辨率分别为50 nm、0.5μm和1μm，所使用的X射线能量一般为120～160 keV，纳米CT的X射线能量更低，为8 ke V。此外，美国FEI公司最新推出了分辨率高达0.8～3.5 nm级聚焦离子束扫描电镜FIBSEM（Focused Ion Beam Scanning Electron Microscopes），通过对岩心表面原子级别的逐层切割输出几个纳米级的高清图像。

（2）发展了模拟沉积与成岩过程的数字岩心构建方法。迄今已发展了多种较为成熟的数值重建方法。图2给出了沉积过程重建方法的示意图。该方法综合考虑了形成沉积岩的颗粒尺寸分布、颗粒形状、压实与胶结及次生加大等成岩过程，能够得到与致密砂岩类似的孔隙格架。

（3）孔隙格架的数学表征方法不断完善。岩心CT扫描的图像是灰度图像，灰度级别反映了X射线衰减程度，它主要受矿物密度和X射线传播距离的影响。致密砂岩一般具有岩屑矿物种类杂、含量多的特点［3］，不同类型矿物的密度不同，因而其CT图像一般具有较宽的灰度分布范围。通过图像处理技术可较好地将孔隙和骨架进行分割，得到只包含孔隙（像素值为0）、骨架背景（像素值为1）的二值图像。

基于分割后的二值图像，目前已形成多种比较有效识别孔隙和喉道的算法，如多向扫描法、中轴线法、多面体法和最大球法。应用这些处理方法可以建立具有真实拓扑结构的孔隙格架。图3为最大球法区分孔隙和喉道的示意图。在二值图像中任意孔隙位置放置一个球体，不断增加球体半径直至遇到孔隙壁为止，在每一个孔隙像素位置均获得这样的球体，用局部最大球表示孔隙体［见图3（a）］，连接相邻孔隙体的所有球体表示喉道［见图3（b）］，从而得到三维孔隙格架信息数据［见图3（c）］。

引入欧拉数表征三维孔隙网络的连通性，

式中，N为网络中直径大于某一尺寸的所有孔隙个数；C为大于该尺寸的所有连通喉道个数；H为溶洞个数，砂岩中一般取0；V为样品体积。

图2 过程模拟法构建数字孔隙格架的方法原理示意图

图3 最大球法区分孔隙、喉道的原理示意图

分析式（1）可知，当网络中连通喉道个数C大于孔隙个数N时，λ为负值，表示连通程度高；反之，当N＞C时，网络中存在孤立孔隙，λ为正值，连通程度低。通过不断改变观测尺度得到一组λ值，绘制出如图4所示的欧拉特征图谱。图4中欧拉曲线与横坐标的交点对应观测尺寸定义为欧拉特征半径d0，单位为μm，它是表征复杂砂岩孔喉网络连通程度的临界值，尺寸小于d0的孔喉连通性好，大于d0的孔喉连通性差。欧拉特征半径是反映复杂孔隙结构碎屑岩储层渗流能力和储层品质的关键参数。

图4 利用欧拉数表示孔隙连通性的原理示意图

除欧拉数以外，描述孔隙结构的参数还包括连通孔隙体积比、孔隙尺寸分布函数、局部孔隙度分布函数和局部渗流概率函数等［7］。

（4）形成多种基于数字孔隙格架的岩石物理属性模拟方法。基于三维数字孔隙格架模拟储层渗透性、电性主要采用基尔霍夫电路节点法、格子玻尔兹曼方法、随机步行者法等，核磁共振性质的模拟采用随机步行者法，岩石弹性力学的模拟主要采用有限元法逐像素点分析骨架在应力作用下的应变分布。基尔霍夫电路节点法的基本原理如图5所示。

对于复杂电路网络，基尔霍夫电流定律（KCL）指出，对于任意一个闭合回路或闭合区域，其中的各段电压和电流满足一定的约束关系［见图5（a）］，是电路分析的基本理论依据。根据该原理，在三维孔隙格架中［见图5（b）］，假设两端施加一定电压ΔV，应用欧姆定律对于任一孔隙节点可以计算出其电流强度I［见图5（c）］，相邻孔隙之间满足KCL定律。对所有孔隙节点都得到一个方程，求解超大方程组就可以计算电阻率。

图5 基尔霍夫电路节点法模拟电阻率的原理示意图

上述模拟过程中如果假设所有孔隙均饱含水（电阻率为Rw），则模拟计算的是岩石饱含水的电阻率R0。如果进一步假设孔隙由大到小依次饱含油，则依次计算出不同含水饱和度对应的电阻率Rt，从而模拟得到岩石的I－Sw关系（见图6）。图6中，圆圈代表某一低孔隙度低渗透率砂岩样品的I－Sw模拟结果，红、蓝、绿色3条直线分别代表饱和度指数n分别取1、2、3的理论结果，表明该样品的ISw关系曲线呈非线性变化。但为了简化，其n值也可以近似取2.0。

格子玻尔兹曼方法是采用人工微观模型，将流体抽象为大量的微观粒子并按照某种简单规则在离散的格子上碰撞和迁移。通过对粒子的运动进行统计分析可以模拟流体运动的宏观特性。基于这一原理，该方法常被用于模拟岩石的电性和渗透性，图7为一种格子玻尔兹曼方法模拟计算常用的D3Q19计算模型［9］，它假设粒子在每个格点上可能存在19个迁移方向，实际模拟运算时须对每一个方向进行分析计算。

根据以上介绍可以看出，基尔霍夫方法是通过简化为孔隙、喉道的三维格架建立方程组计算电阻率，每一节点上电流的方向是确定的，具体取决于与孔隙连接的喉道个数，即使在复杂孔隙结构的致密砂岩中孔喉比一般都比较小。格子玻尔兹曼方法在每一节点都要考虑多达19个方向，因此，基尔霍夫方法的模拟计算速度要明显占优势，但二者的计算精度基本相当。

图7 格子玻尔兹曼模拟常用的D3Q19计算模型示意图

随机步行者方法常用来模拟岩石的核磁共振特性，基本原理是在孔隙格架中随机放置若干个质子，这些质子以布朗运动与孔隙壁碰撞并湮没或反弹。通过设置一定的时间步长计算质子逐渐消亡导致的宏观磁矩不断减小的过程就可以模拟核磁共振信号。该方法原理简单，只要CT分辨率足够，其模拟的结果与实测岩心的核磁共振信号基本相当。

3 测井数字岩石物理研究面临的挑战

（1）分辨率与储层非均质性之间存在矛盾严重影响模拟精度。致密砂岩储层主要发育纳米级孔喉［12］，为提高分辨率以观测更细微孔隙就必须减小样品尺寸，纳米级扫描要求样品尺寸小于0.1 cm，而对具有强烈微观非均质性的砂岩而言，小尺寸样品的模拟结果与储层宏观岩石物理性质相差较大。另一方面，大尺寸样品的CT扫描只能揭示控制渗流的主要喉道，却无法刻画对电性影响显著的微细孔喉，该矛盾导致的结果可以用图8解释。图8是4块直径2.54 cm的砂岩样品采用基尔霍夫电路节点法模拟的渗透率、地层因素与实验测量结果对比。可以看出，在14.27μm的CT分辨率下（对应样品直径2.54 cm）依据所能够观察到的大孔喉模拟的渗透率与实验值在数量级上基本接近［见图8（a）］，说明大孔喉是控制储层渗流能力的关键；但地层因素F的模拟结果与实验值相差达2～3个数量级［见图8（b）］，这主要是由于其中大量纳米级微细孔喉未被CT探测到、模拟计算时其导电作用被忽略所致。

（2）图像分割与构建数字孔隙格架过程的多解性强。图9（a）、（b）分别是1块常规中高孔隙度渗透率砂岩样品和1块致密砂岩样品2 000张CT图片灰度分布对比。表明致密砂岩的图像灰度分布表现为强烈的非均质性，图像处理过程中很难选取准确的灰度阈值以实现二值化分割。图9（c）显示典型的致密砂岩中发育极不规则的孔隙和喉道，很难用三角形、球形等简单的几何特征描述，因此在致密砂岩中，基于图像准确构建三维数字孔隙格架的过程存在很强的多解性，其结果也将影响模拟精度。

图8 4块样品基尔霍夫电路节点法模拟渗透率、地层因素的精度分析

图9 典型致密砂岩CT图片灰度分布与复杂几何形状的孔隙特征

（3）缺乏能够模拟黏土矿物影响的方法模型。黏土矿物一般以分散状、结构状或层状3种形式分布在砂岩储层中，不同的黏土矿物类型和分布形式对储层渗透率、电阻率的影响非常显著［13－14］，特别是类似于渤海湾盆地歧口凹陷中深层沙河街组、岩性较细、泥质含量和压实程度高、成岩作用强的致密砂岩储层［15］，这种影响就更加突出。但综合目前国内外的研究情况，考虑储层润湿特征影响多相流的模型较多，但尚欠缺考虑黏土矿物影响岩石物理属性的数字孔隙格架构建及相关模拟方法的研究。

（4）现有模拟方法的时效性亟需优化完善。受数据量大、计算过程复杂等因素制约，现有模拟方法耗时长，采用当前市场上最好的高性能微机工作站模拟运算1次一般需要近20 h，时效性差是制约现有方法推广应用的关键，亟需进行算法的优化完善，开发并行算法程序。

4 研究方向展望

（1）趋势云分析技术。针对分辨率与储层非均质性之间的矛盾，对一块常规尺寸致密砂岩样品先进行粗扫，根据粗扫图像选择孔隙发育程度或连通程度不同的多个部位切割出小尺寸样品再进行精细扫描，分别基于粗扫图像和细扫图像模拟储层的岩石物理参数，最后结合不同分辨率的模拟结果进行趋势分析，判断渗透率、电阻率等属性的变化规律和主要控制因素。

（2）多分辨率三维孔隙网络整合。与趋势云分析技术相对的是多分辨率孔隙网络整合，基本思路是对大尺寸岩心按照相同的思路进行逐级精细扫描，对每个级别的扫描图像分别统计其欧拉数及孔隙尺寸分布函数，最后在不改变拓扑结构特征的前提下将高分辨率图像整合到粗扫图像中，最终获得能够反映大尺寸岩心的高分辨率图像，据此进行模拟分析。

5 结论与认识

（1）为满足致密砂岩等非常规储层岩石物理研究需要，急需深化和完善测井数字岩石物理模拟与分析技术。

（2）以致密砂岩为代表的复杂储层岩石物理研究不能沿用常规思路，将高精度实验模拟与数值模拟相结合是一种必然趋势。

（3）当前的测井数字岩石物理模拟技术尚不完备，不能完全代替实验模拟，亟需进行理论模型和计算方法的优化完善。

［1］ 翟光明，何文渊.从区域构造背景看我国油气勘探方向 ［J］.中国石油勘探，2005，2（10）：1－8.

［2］ 邱中建，康竹林，何文渊.从近期发现的油气新领域展望中国油气勘探发展前景 ［J］.石油学报，2002，4（23）：1－6.

［3］ 赵政璋，杜金虎，等.致密油气 ［M］.北京：石油工业出版社，2012.

［4］ Purcell W R.Capillary Pressures——Their Measurement Using Mercury and the Calculation of Permeability［C］∥ Trans AIME，1949，186：39－48.

［5］ 陶果，岳文正，谢然红，等.岩石物理的理论模拟和数值实验新方法 ［J］.地球物理学进展，2005，1（20）：4－11.

［6］ Fredrich J T，Menendez B，Wong T F.Imaging the Pore Structure of Geomaterials ［J］.Science，1995，268：276－279.

［7］ 赵秀才.数字岩心及孔隙网络模型重构方法研究［D］.青岛：中国石油大学（华东），2009.

［8］ 姚军，赵秀才，衣艳静，等.数字岩心技术现状及展望 ［J］.油气地质与采收率，2005，12（6）：52－54.

［9］ 孙建孟，姜黎明，刘学锋，等.数字岩心技术测井应用与展望 ［J］.测井技术，2012，36（1）：1－7.

［10］王晨晨，姚军，杨永飞，等.基于CT扫描法构建数字岩心的分辨率选取研究 ［J］.科学技术与工程，2013，13（4）：1050－1052.

［11］Jin Guodong，Carlos Torres－Verdin，Chun Lan.Porelevel Study of Grain－shape Effects on Petrophysical Properties of Porous Media［C］∥SPWLA，2009，Paper TTT.

［12］邹才能.非常规油气地质 ［M］.北京：地质出版社，2011.

［13］雍世和，张超谟.测井数据处理与综合解释 ［M］.东营：石油大学出版社，1996.

［14］欧阳健，毛志强，修立军，等.测井低对比度油层成因机理与评价方法 ［M］.北京：石油工业出版社，2009.

［15］中国石油勘探与生产分公司.低阻油气藏测井评价技术及应用 ［M］.北京：石油工业出版社，2009.