邹友龙，谢然红，郭江峰，谭茂金，胡法龙，李潮流，李长喜，周灿灿

（1.中国石油大学油气资源与探测国家重点实验室，北京102249；2.中国石油大学地球探测与信息技术北京市重点实验室，北京102249；3.中国地质大学（北京）地球物理与信息技术学院，北京100083；4.中国石油勘探开发研究院测井与遥感技术研究所，北京100083）

致密储层数字岩心重构及核磁共振响应模拟

邹友龙1，2，谢然红1，2，郭江峰1，2，谭茂金3，胡法龙4，李潮流4，李长喜4，周灿灿4

（1.中国石油大学油气资源与探测国家重点实验室，北京102249；2.中国石油大学地球探测与信息技术北京市重点实验室，北京102249；3.中国地质大学（北京）地球物理与信息技术学院，北京100083；4.中国石油勘探开发研究院测井与遥感技术研究所，北京100083）

基于物理过程法模拟沉积岩的沉积、压实和胶结过程，构建致密岩石的三维数字岩心。利用随机行走法模拟不同成岩过程岩石的核磁共振响应以及不同润湿性岩石孔隙中流体的核磁共振响应。模拟结果表明，岩石胶结成岩后孔隙半径减小导致核磁共振横向弛豫（T2）分布的峰值向短弛豫方向移动，流体视扩散系数略小于自由扩散系数；致密岩石中随润湿相流体饱和度减小，润湿相T2分布向短弛豫方向偏移，受限扩散越来越明显。

致密储层；数字岩心重构；孔隙尺度；核磁共振响应模拟；流体受限扩散

致密储层岩石物理实验困难，核磁共振（NMR）响应不同于常规储层，需要开展致密储层数字岩心的重构方法及孔隙尺度下NMR响应的数值模拟研究。数字岩心重构方法主要分为4类：实验方法［1-3］、统计方法［4-8］、几何方法［9-10］和过程法［11-17］。实验方法存在实验条件限制和分辨率的问题。统计方法构建的多孔介质间的几何特征存在明显的差异。几何方法适用于模拟松散的沉积岩。过程法模拟岩石微观结构时考虑了岩石形成的基本过程：沉积、压实和成岩过程［11-12］，包括几何过程法［11］和物理过程法［12］。几何过程法只考虑了岩石形成的最终几何形态结果，而物理过程法还将颗粒的物理受力运动考虑到岩石成岩的模拟过程中，模拟的数字岩心与实际岩石更加接近。孔隙尺度下的NMR响应模拟通常采用随机行走法［18-21］，它相对于有限元方法［22］和有限差分方法［23］具有更好的灵活性，且操作简单容易实现。笔者基于物理过程法重构致密储层数字岩心，采用随机行走法模拟岩石的NMR响应。

1 数字岩心重构

1.1 物理过程法原理

物理过程法［12］根据颗粒受力和力矩平衡方程，利用有限差分原理计算每一时刻的颗粒运动状态参数（受力、位置、速度等）。

1.1.1 颗粒受力分析

整个岩石颗粒系统的行为主要依赖于颗粒之间的相互作用力，多颗粒系统中单个颗粒均满足牛顿运动方程：

式中，Fi和Mi分别为颗粒在i（i=x，y，z）方向上受到的合力与力矩；位移加速度u¨i和角加速度θ¨i分别为位置ui与旋转角θi对时间求二次导数。

颗粒受多种力作用，包括：体力Fb，i，作用于颗粒上的外力与重力；接触力Fc，i，颗粒与颗粒和颗粒与边界之间的接触力；阻尼力Fd，i，与颗粒速度成正比的阻尼力Fdm，i和与颗粒刚度成正比的阻尼力

接触力Fc，i由法向接触力与切向接触力构成，即法向接触力Fnc，i可由下式计算：

阻尼力Fdm，i表征颗粒在黏性流体中运动时的能量耗散，阻尼力Fdk，i表征颗粒与颗粒、颗粒与边界碰撞时的能量耗散，两种阻尼力可由下式计算：

式中，α和β为对应的阻尼系数；m为颗粒质量；.ui为颗粒的位移速度；.ur，i为颗粒与颗粒接触、颗粒与边界接触时的相对位移速度；k为刚度系数，包括法向刚度kn和切向刚度kt。

模拟过程中颗粒的动态运动过程在其达到稳定平衡态时停止，通常用以下准则来衡量颗粒系统的平衡态：所有颗粒所受合力的均值与所有颗粒接触点平均接触力的比值或者最大合力与最大接触力的比值满足先验条件（如小于某一阈值0.01）。处于平衡态时，颗粒所受合力趋近于零，颗粒运动速度也趋近于零，可认为系统达到动态平衡，从而结束模拟。

1.1.2 数值模拟方法

利用有限差分原理对公式（1）和公式（2）采用固定时间步长Δt求解，t时刻颗粒的加速度和速度可用下式近似计算：

其中

最后，颗粒的位置ui与旋转角θi可由下式计算：

对于时间步长的选择，ITASCA［24］给出了具体的选取方法，其在每个时间步长之后都重新计算新的时间步长，计算过程复杂且计算量大，为了操作简单可以选取某一足够小的常数值作为时间步长。

1.2 岩心重构过程

1.2.1 沉积过程

岩石颗粒的沉积包括两种处理方法：①相互之间不重叠的所有颗粒同时沉积；②颗粒逐个沉积，当前一颗粒沉积达到平衡态后，新的颗粒开始沉积，这种方法相对耗时，主要用于在已有颗粒堆积体上添加新的颗粒。

假设岩石颗粒半径服从截断正态分布，半径最小值为20 μm，最大值为40 μm，颗粒的密度为2.65 g/cm3，法向刚度kn和切向刚度kt均为8.0×1010N/ m。采用颗粒同时沉积的方法，首先根据颗粒半径分布在0.6 mm×0.6 mm×1.4 mm的立方体容器中生成随机分布的2000个颗粒，且颗粒与颗粒、颗粒与边界之间均不重叠，如图1（a）所示。然后，模拟在重力作用下颗粒沉积下降，如图1（b）所示。模拟过程利用式（3）～（14）计算颗粒在每一时刻的位置、速度及受力情况，其中时间步长为3.77×10-10s。颗粒在沉积的过程中受到流体阻尼力和碰撞产生的阻尼力作用，能量逐渐耗散，速度逐渐趋近于零，最终达到平衡态完成沉积过程，结果如图1（c）所示。

图1 数字岩心构建Fig.1 Digital core reconstruction

1.2.2 压实过程

随着岩石颗粒的不断堆积，岩石的上部会受到巨大的地层压力作用。为此，通过在颗粒堆积体的上部施加一道不断向下移动的平板来实现岩石的压实过程，数值计算方法与沉积过程的模拟方法相同，但岩石的上部受上覆地层压力作用，如图1（d）所示。压实过程通过岩石应变进行控制，随着压力的不断增大，平板的垂向位移不断减小。压实后的岩石孔隙度控制在30%左右较为合理。本文中取应变率为1.2×10-2s-1，当颗粒堆积体在垂向上应变为0.12时停止，完成岩石的压实过程，结果如图1（e）所示。

1.2.3 胶结过程

当流体流过岩石孔隙时，可能发生化学成岩作用。在颗粒的表面将发生岩石的溶解和矿物的沉淀，使得这些单个松散的颗粒胶结在一起形成多孔的固体岩石。采用如下颗粒表面胶结物生长公式［12］：

式中，L（r）为颗粒中心到胶结物表面的距离；R为颗粒半径；¯R为原始颗粒半径分布的平均半径；l（r）为颗粒表面与多面体（颗粒内切于该多面体）平面的距离；ξ为控制颗粒半径影响的参数；ζ为控制颗粒优势胶结生长方向的参数；κ为控制孔隙度的参数。

为了获得致密岩石，利用公式（15）取ξ=1.0，ζ =0.5，κ=1.0，对图1（e）压实后的岩石模拟其颗粒表面的石英胶结物生长作用，模拟结果如图1（f）所示。颗粒堆积体胶结之后，其孔隙度进一步减小，且颗粒的形状也发生了很大变化，颗粒表面呈不规则的多边弧形。

2 NMR响应模拟方法

2.1 随机行走法

通常情况下，孔隙介质中流体的核磁共振横向弛豫时间T2可写成如下形式：

式中，T2S为表面弛豫，主要与孔隙结构和岩石类型有关；T2B为自由弛豫，主要与流体性质有关，与孔隙结构无关；T2D为扩散弛豫，主要与磁场梯度G、回波间隔TE和流体的扩散系数D有关；γ为旋磁比。

NMR的弛豫信号强度M（t）随时间t的变化可由下式计算：

式中，MS（t）为t时刻表面弛豫信号强度；MD（t）为t时刻扩散弛豫信号强度。

随机行走法模拟岩石孔隙中流体的表面弛豫信号强度MS与扩散弛豫信号强度MD的具体步骤如下：

（1）将固定数目的质子随机分布于岩石孔隙中。

（2）计算质子与最近固体表面的距离d，当d（＜3ε，ε为传统方法的扩散半径）较小时，采用传统方法，即扩散半径r=ε；当d（≥3ε）较大时，采用第一旅行时方法［25］，即扩散半径r=d。

（3）计算时间间隔Δt=r2/6D和质子下一时刻的位置［x（t+Δt），y（t+Δt），z（t+Δt）］：

其中，cos φ从［-1，1］中随机选取；θ从［0，2π］中随机选取；［x（t），y（t），z（t）］为当前时刻质子的位置。

（4）判断质子是否与固体表面碰撞，若发生碰撞，有以下两种处理方法：

①质子以δ概率殒灭，若质子未陨灭，则质子发生反弹［21］；

②质子弛豫强度以（1-δ）或exp（-δ）衰减并发生反弹［19］；

概率δ的计算公式为δ=（2ρr）/（3D）.

其中，ρ为表面弛豫率；r为扩散半径。

第②种处理方法相对第①种处理方法运算速度更慢，但模拟结果相对光滑。当质子数足够多时，两种方法模拟结果之间的误差可以忽略，因此通常都采用第①种处理方法。

（5）判断质子是否走出岩石，若质子走出岩石，则质子下一时刻的位置为在岩石的相反面随机选取一个孔隙位置。

（6）计算质子的相位偏移φ，

（7）当t=nTE时，记录质子的扩散弛豫强度（相位余弦）和表面弛豫强度得到质子信号强度总和。重复步骤（2）～（6），直至采样时间大于设置的阈值时终止。

2.2 方法验证

假设在均匀场中，球形孔隙饱含水，孔隙半径为5 μm，表面弛豫率为30 μm/s，水的扩散系数为2.1 μm2/ms，此时满足快扩散条件，可知其表面弛豫时间理论值为55.56 ms。图2（a）为模拟的表面弛豫衰减信号（蓝色实线）与理论值（红色点线）的对比图，从图中可以看到，两者完全重合，验证了随机行走法模拟表面弛豫的有效性。

在0.3 T/m的梯度场下，模拟自由状态下水的扩散弛豫并与理论值对比，其中水的扩散系数为2.1 μm2/ms，回波间隔为1.8 ms，自由扩散条件下水的扩散弛豫时间理论值为273.82 ms。从图2（b）可以看到，模拟值（蓝色实线）与解析值（红色点线）吻合，验证了随机行走法模拟扩散弛豫的有效性。

图2 随机行走法模拟值与理论值对比Fig.2 Simulated values using random-walk method versus theoretical values

3 NMR响应模拟结果

3.1 不同成岩过程岩石的NMR响应

为了解不同成岩过程中岩石的NMR响应，分别从各成岩过程模拟的岩石中取出一个300 μm× 300 μm×300 μm的多孔介质进行NMR响应模拟。图3（a）、（b）和（c）分别为从沉积、压实和胶结过程的岩石中提取的多孔介质，其孔隙度分别为40.86%、33.58%和9.14%。图3（d）、（e）和（f）分别为沉积、压实和胶结过程多孔介质的切片图。从图中可以看到，随着成岩过程的进行，岩石颗粒由最初沉积过程的点接触，到压实和胶结过程变成了线接触。沉积和压实过程岩石孔隙结构差别不是特别明显，而岩石胶结后颗粒呈不规则多面体，孔隙结构复杂，接近实际地层岩石。

假设岩石饱含水，其表面弛豫率为30 μm/s，水的自由扩散系数为2.1 μm2/ms，体弛豫为3.0 s。在0.3 T/m的梯度场中分别设置回波间隔为0.45、0.9、1.8、3.6和7.2 ms，对不同成岩过程岩石（图3（a）、（b）和（c））各自模拟得到5组回波串。对模拟的回波串联合反演后得到图4所示的不同成岩过程岩石孔隙流体扩散系数D与横向弛豫时间T2的D-T2分布，图中白色实线代表水的自由扩散系数。图5所示为将不同成岩过程岩石的D-T2分布分别投影至T2轴和D轴。从图4和图5可以看出，沉积和压实过程岩石的D-T2分布基本一致，T2分布的峰值约在190 ms。胶结过程岩石的T2分布峰值则约在65 ms（图5（a）），且视扩散系数相对自由扩散系数略微偏小（图5（b）），这是由于胶结成岩后岩石的孔隙半径变小，导致表面弛豫加快，且随孔隙半径减小，流体扩散受限，使得视扩散系数略小于自由扩散系数造成的。

3.2 两相流的NMR响应

为了观察两相流的NMR响应，考虑水湿和油湿两种情况，分别模拟润湿相饱和度不同的岩石NMR响应。对图3（c）所示的胶结后的岩石分别构造润湿相饱和度分别为25%、50%和75%的数字岩心，如图6所示，其中黑色为骨架，红色为润湿相流体，蓝色为非润湿相流体。

设水的表面弛豫率为30 μm/s，自由扩散系数为2.1 μm2/ms，体弛豫为3.0 s；中等黏度油的表面弛豫率为10 μm/s，自由扩散系数为0.11 μm2/ms，体弛豫为0.2 s。在0.3 T/m的梯度场中，分别设置回波间隔为0.45、0.9、1.8、3.6和7.2 ms，分别对润湿相饱和度不同的岩石模拟得到5组回波串并反演得到其D-T2分布。

图3 不同成岩过程岩石的多孔介质与切片图Fig.3 Porous mediums and slices of rock of different diagenetic processes

图4 不同成岩过程岩石的D-T2分布Fig.4 D-T2distributions of rock of different diagenetic processes

图5 不同成岩过程岩石的D-T2分布分别在T2和D轴的投影分布Fig.5 D-T2distributions of rock of different diagenetic processes projected on T2and D axis

图6 润湿相饱和度不同的数字岩心Fig.6 Digital cores with different saturations of wetting phase

图7 和图9分别为水润湿和油润湿情况下，不同饱和度岩石的D-T2分布，图中白色实线代表水的自由扩散系数，红色实线表示油的自由扩散系数。图8所示为水润湿情况下，不同含水饱和度岩石的D-T2分布分别投影至T2轴和D轴。从图8可以看到，随着含水饱和度减小，水的T2分布向短弛豫时间偏移（图8（a）），视扩散系数减小（图8（b）），这是由于水占据的孔隙空间的表面积与体积的比值（S/V）越来越大，流体受限扩散越来越明显；然而非润湿相油的T2和视扩散系数分布的峰值位置基本没有变化。图10所示为油润湿情况下，不同含油饱和度岩石的D-T2分布分别投影至T2轴和D轴。从图10可以看出，随着含油饱和度减小，油的T2分布向短弛豫时间偏移，油的视扩散系数减小，这是由于油占据的孔隙空间的S/V增大，油的受限扩散明显。从图7～10可以发现，随润湿相流体饱和度减小，润湿相受限扩散增强。

图7 水润湿情况下不同含水饱和度岩石的D-T2分布Fig.7 D-T2distributions of rocks with different water saturations at water-wet condition

图8 水润湿情况下不同含水饱和度岩石的D-T2分布分别在T2和D轴的投影Fig.8 D-T2distributions of rocks with different water saturations projected on T2and D axis at water-wet condition

图9 油润湿情况下不同含油饱和度岩石的D-T2分布Fig.9 D-T2distributions of rocks with different oil saturations at oil-wet condition

图10 油润湿情况下不同含油饱和度岩石的D-T2分布分别在T2和D轴的投影Fig.10 D-T2distributions of rocks with different oil saturations projected on T2and D axis at oil-wet condition

4 结 论

（1）物理过程法能够重构致密岩石的三维数字岩心，岩石成岩过程中孔隙度和孔隙半径逐渐减小，最终获得的数字岩心的孔隙表面呈不规则的几何形状，与实际岩心接近。

（2）岩石胶结成岩后孔隙半径减小，导致NMR弛豫信号衰减速率加快，T2分布的峰值向短弛豫方向移动，流体视扩散系数相对自由扩散系数略微偏小。

（3）致密岩石中润湿相流体随饱和度减小，润湿相的T2分布向短弛豫方向偏移，流体视扩散系数减小，受限扩散越来越明显。

［1］ SERIGHT R S，LIANG J，LINDQUIST W B，et al.Use of x-ray computed microtomography to understand why gels reduce relative permeability to water more than that to oil［J］.Journal of Petroleum Science and Engineering，2003，39（3/4）：217-230.

［2］ SOK R M，VARSLOT T，GHOUS A，et al.Pore scale characterization of carbonates at multiple scales：integration of micro-CT，BSEM，FIBSEM［J］.Petrophysics，2010，51（6）：379.

［3］ TIWARI P，DEO M，LIN C L，et al.Characterization of oil shale pore structure before and after pyrolysis by using X-ray micro CT［J］.Fuel，2013，107：547-554.

［4］ MANWART C，TORQUATO S，HILFER R.Stochastic reconstruction of sandstones［J］.Physical Review E，2000，62（1）：893.

［5］ OKABE H，BLUNT M J.Pore space reconstruction using multiple-point statistics［J］.Journal of Petroleum Science and Engineering，2005，46（1）：121-137.

［6］ HAJIZADEH A，SAFEKORDI A，FARHADPOUR F A. A multiple-point statistics algorithm for 3D pore space reconstruction from 2D images［J］.Advances in Water Resources，2011，34（10）：1256-1267.

［7］ SHEIDAEI A，BANIASSADI M，BANU M，et al.3-D microstructure reconstruction of polymer nano-composite using FIB-SEM and statistical correlation function［J］. Composites Science and Technology，2013，80：47-54.

［8］ 王晨晨，姚军，杨永飞，等.碳酸盐岩双孔隙数字岩心结构特征分析［J］.中国石油大学学报（自然科学版），2013，37（2）：71-74.WANG Chenchen，YAO Jun，YANG Yongfei，et al. Structure characteristics analysis of carbonate dual pore digital rock［J］.Journal of China University of Petroleum（Edition of Natural Science），2013，37（2）：71-74.

［9］ JIA X，WILLIAMS R A.A packing algorithm for particles of arbitrary shapes［J］.Powder Technology，2001，120：175-186.

［10］ JIA X，GAN M，WILLIAMS R A，et al.Validation of a digital packing algorithm in predicting powder packing densities［J］.Powder Technology，2007，174（1）：10-13.

［11］ BAKKE S，ØREN P E.3-D pore-scale modelling of sandstones and flow simulations in the pore networks［R］.SPE 35479，1997.

［12］ JIN G D，PATZEK T W，SILIN D B.Physics-based reconstruction of sedimentary rocks［R］.SPE 83587，2003.

［13］ JIN G D，LAN C，KOPAL M，SHAO W，et al.Predicting formation permeability from rock models constructed using log and/or geological data：SPWLA 52nd Annual Logging Symposium，14-18 May，Colorado Springs，Colorado，2011［C］.Texas：Society of Petrophysicists and Well-Log Analysts，c2011.

［14］ 刘洋.过程法构建数字岩心技术［D］.青岛：中国石油大学，2007. LIU Yang.Construction of digital core by process-based simulation method［D］.Qingdao：China University of Petroleum，2007.

［15］ 刘学锋.基于数字岩心的岩石声电特性微观数值模拟研究［D］.青岛：中国石油大学，2010. LIU Xuefeng.Numerical simulation of elastic and electrical properties of rock based on digital cores［D］. Qingdao：China University of Petroleum，2010.

［16］ 盛金昌，叶辉明，周治荣，等.三维数字岩心的重构及孔隙网络模型的孔隙度分析［J］.水电能源科学，2012，30（10）：65-68. SHENG Jinchang，YE Huiming，ZHOU Zhirong，et al. Reconstruction of 3D digital core and porosity analysis of pore network model［J］.Water Resources and Power，2012，30（10）：65-68.

［17］ 孙建孟，赵建鹏，闫伟超，等.应用核磁T2谱与数字岩心技术计算粒度分布方法［J］.中国石油大学学报（自然科学版），2013，37（3）：57-62. SUN Jianmeng，ZHAO Jianpeng，YAN Weichao，et al. Calculation of grain size distribution using NMR T2spectrum and digital rock technology［J］.Journal of China University of Petroleum（Edition of Natural Science），2013，37（3）：57-62.

［18］ TOUMELIN E，TORRES-VERDIN C，CHEN S H.Modeling of multiple echo-time NMR measurements for complex pore geometries and multiphase saturations［J］.SPE Reservoir Evaluation&Engineering，2003，6（4）：234-243.

［19］ TOUMELIN E，TORRES-VERDIN C，SUN B Q，et al. Random-walk technique for simulating NMR measurements and 2D NMR maps of porous media with relaxing and permeable boundaries［J］.Journal of Magnetic Resonance，2007，188：83-96.

［20］ JIN G D，TORRES-VERDIN C，TOUMELIN E.Comparison of NMR simulations of porous media derived from analytical and voxelized representations［J］.Journal of Magnetic Resonance，2009，200：313-320.

［21］ TALABI O，ALSAYARI S，IGLAUER S，et al.Porescale simulation of NMR response［J］.Journal of Petroleum Science and Engineering，2009，67：168-178.

［22］ HAGSLATT H，JONSSON B，NYDEN M，et al.Prediction of pulsed field gradient NMR echo-decays for molecules diffusing in various restrictive geometries，simulations of diffusion propagators based on a finite element method［J］.Journal of Magnetic Resonance，2003，161（2）：138-147.

［23］ ZIENTARA G Z，FREED J H.Spin-echoes for diffusion in bounded，heterogeneous media：numerical study［J］. J Chem Phys，1980，72：1285-1292.

［24］ ITASCA.PFC3D：theory and background［M］.1st ed. Minnesota：Itasca Consulting Group Inc，1999.

［25］ ZHENG L H，CHIEW Y C.Computer simulation of diffusion-controlled reactions in dispersions of spherical sinks［J］.The Journal of Chemical Physics，1989，90（1）：321-327.

（编辑 修荣荣）

Reconstruction of digital core of tight reservoir and simulation of NMR response

ZOU Youlong1，2，XIE Ranhong1，2，GUO Jiangfeng1，2，TAN Maojin3，HU Falong4，LI Chaoliu4，LI Changxi4，ZHOU Cancan4
（1.State Key Laboratory of Petroleum Resources and Prospecting in China University of Petroleum，Beijing 102249，China；2.Key Laboratory of Earth Prospecting and Information Technology in China University of Petroleum，Beijing 102249，China；3.School of Geophysics and Information Technology in China University of Geosciences，Beijing 100083，China；4.Department of Well Logging and Remote Sensing Technology，PetroChina Research Institute of Petroleum Exploration&Development，Beijing 100083，China）

A 3D digital core model of tight rock was constructed by simulating the sedimentation，compaction and cementation processes of sedimentary rocks.The nuclear magnetic resonance（NMR）responses of different diagenetic processes of the sedimentary rock，and the NMR responses of fluid in tight rock under different wetting conditions were simulated using a random walk method.The simulation results show that the decreasing of pore radius caused by cementation process leads to the peaks of NMR transverse relaxation time（T2）distribution shifting to short relaxation time，and results in slightly smaller apparent diffusion coefficient of fluids than its free diffusion coefficient.In tight rock，when the saturation of wetting phase decreases，the peaks of its T2distributions shift to short relaxation time and the restricted diffusion becomes more obvious.

tight reservoir；digital core reconstruction；pore-scale；NMR response simulation；restricted diffusion of fluid

P 313.1

A

邹友龙，谢然红，郭江峰，等.致密储层数字岩心重构及核磁共振响应模拟［J］.中国石油大学学报（自然科学版），2015，39（6）：63-71.

ZOU Youlong，XIE Ranhong，GUO Jiangfeng，et al.Reconstruction of digital core of tight reservoir and simulation of NMR response［J］.Journal of China University of Petroleum（Edition of Natural Science），2015，39（6）：63-71.

1673-5005（2015）06-0063-09

10.3969/j.issn.1673-5005.2015.06.008

2014-12-18

国家自然科学基金委员会-中国石油天然气集团公司石油化工联合基金（U1262114）；高等学校博士学科点专项科研基金（20130007110012）

邹友龙（1988-），男，博士研究生，研究方向为岩石物理、核磁共振测井方法及应用。E-mail：zoyolo_ok@126.com。