基于非常快速模拟退火法的页岩岩心双重区域重构方法

邹孟飞，隋微波，王旭东，张　爽

（中国石油大学（北京）石油工程学院，北京102249）

基于非常快速模拟退火法的页岩岩心双重区域重构方法

邹孟飞，隋微波，王旭东，张爽

（中国石油大学（北京）石油工程学院，北京102249）

页岩储层微观结构在储层开发评价和渗流机理研究方面的重要性日益凸显。但是页岩储层具有低孔、特低渗透、强非均质性和各向异性等特征，应用常规物理实验方法在研究微观渗流方面存在明显的局限性。数字岩心技术已成为研究储层微观结构的重要方法之一，重构方法有其较好的准确性和可靠性。以真实页岩岩心扫描电镜二维图像为基础，利用软件对图像进行分析和数据提取，应用传统模拟退火法重构页岩岩心。针对重构中的问题总结出基于非常快速模拟退火法和双重区域的页岩岩心重构方法。通过改进算法，大幅度的提高重构速度；针对页岩孔隙分布特征，对页岩进行分区域重构，使重构结果更加接近真实页岩岩心。

数字岩心重构页岩双重区域非常快速模拟退火法模拟退火法

现今油气资源日益紧缺，非常规油气资源尤其是页岩气的开采显得尤为重要［1-4］。页岩储层的微观孔隙结构是进行储层评价的基础，也是计算页岩气产能的关键。页岩气具有低孔、特低渗透的储层特性［5］，非均质性和各向异性较强，主要发育纳米级的孔喉系统［6］，局部发育微米—毫米级孔隙。现有的室内岩心渗流实验技术对于页岩岩心的测定缺乏统一的实验方法和技术标准。近年来，数字岩心技术在多孔介质微观孔隙结构研究及表征、渗流参数计算数值实验，以及多孔介质微观渗流机理研究方面都体现了较高的研究和应用价值。与常规岩石物理实验相比，数字岩心可利用很小的钻屑，甚至无法进行渗透率测试的被破坏的岩心进行重构，对岩心无损伤，并且可重复进行多次数值实验［7］。但是目前的数字岩心技术是基于砂岩和碳酸盐岩储层建立的，不适用于页岩岩心。为此，笔者针对页岩岩心孔隙类型多样、孔隙尺度相差较大的特点，总结出基于非常快速模拟退火算法的页岩岩心双重区域重构方法，为快速有效地建立页岩数字岩心提供了一种新的思路。

1　传统数字岩心技术

传统数字岩心技术主要分为2大类：直接成像方法和基于薄片分析的图像重构方法。直接成像方法又包括激光共聚焦扫描显微镜、序列切片法和CT扫描法，这3种方法在砂岩和碳酸盐岩中已得到广泛应用。Fredrich［8］等曾用序列切片法建立数字岩心，分辨率可达到亚微米级别，但该方法并不是真正意义上扫描岩心，只能处理岩心薄片，得到的数字岩心厚度有限、规模小，对于研究页岩的微观孔隙结构意义不大。Tomutsa等［9］较早地应用序列切片法建立了砂岩的数字岩心并进行了岩心孔隙结构分析。Curtis等［10］利用聚焦离子束—扫描电镜（FIB-SEM）系统得到页岩的序列切片图片，并重构出页岩三维数字岩心。但序列切片法在岩心切割和抛光过程中耗时巨大且易破坏岩心结构，实用性较差。CT扫描法能提供微米—亚微米级别的高分辨率图像，对岩心无损伤，但其构建数字岩心的造价和实验成本较高，而且现有的Micro-CT分辨率无法满足页岩孔隙纳米级的精度要求。

基于薄片分析的图像重构方法通常只需要几张岩心切片图像即可进行岩心重构，目前已发展形成高斯场法、模拟退火法、过程法、多点统计法和蒙特卡洛—马可夫链随机法等重构方法，同样多用于重构砂岩、碳酸盐岩［11-12］。

高斯场法1974年由Joshi［13］提出并用于数字岩心的重构，受到计算机技术的限制，只建立出二维数字岩心。1984年由Quiblier［14］建立了第1个真正的三维数字岩心。1995年Ioannidis等［15］改进了Quiblier的算法，使建模速度有所提高。但该方法重构结果的连通性和各向异性都较差。

模拟退火法（SA） 1997年由Hazlett［16］首次提出的另一种随机法，并于1998年由Yeong等［17-18］进一步改进。2001年Hidajat［19］综合运用高斯场法和模拟退火法，提高了运算速度，并在当年由Talukdar等重构出轻度固结的玻璃球、北海白垩岩［20-21］。2007年Capek等重构出三氧化二铝样本［22］，是一种有效的非线性组合优化算法，在重构数字岩心时可以考虑任意多的约束条件，并且可以跳出局部极值。虽然在约束函数较多时重构时间较慢，但改进空间大。

过程法1997年由Bakke等［23］提出，并于2002年重构了Fontainebleau砂岩模型，2003年对其算法进一步改进［24］，提高了过程法重构的能力，并建立了Berea砂岩的孔隙图像。该方法主要以岩石的真实形成过程为基础，模拟岩石沉积、压实和成岩过程，建立的岩心各向异性较好。

多点统计法Strebelle等［25-28］首次提出，由Okabe等［29-30］用于重构砂岩和碳酸盐岩岩心的孔隙结构。为了反映真实多孔介质的各向异性，2009年张挺等［31］考虑垂直方向孔隙和骨架结构信息，重构出更符合真实情况的多孔介质。但该方法计算耗时较大。

蒙特卡洛—马尔科夫随机场法（MCMC） 由Wu等［32-33］提出并对其进行改进，将算法由二维扩展为三维，重构出泥土、砂岩和泥岩的三维岩心模型。但该方法建立的数字岩心约束函数少，只能建立非均质性较弱的数字岩心。

2　页岩岩心电镜实验

应用传统的模拟退火法对页岩岩心进行模拟重构时，要获取一定数量的页岩岩心的扫描电镜图像，并对图像进行处理、分析和统计，为岩心重构提供基础数据和图像；并对重构结果进行分析比对。

实验所用页岩岩心来源于四川盆地下志留统龙马溪组，其扫描电镜（SEM）图像来源于中国石油大学（北京）能源材料微结构实验室，利用Quanta200F场发射环境扫描电镜拍摄，该扫描电镜分辨率可达到1.2 nm。页岩储层发育的主要孔隙类型有微裂缝、矿物质孔、有机质纳米孔和生物化石孔，其中有机质纳米孔隙被众多学者认为对页岩气的储集及产能具有重要影响，因此将其作为研究重点。无定形有机质形态特征多是不规则的、渐变的模糊边缘，形成凝块或絮团状，其主要成因是浮游植物残骸沉积形成［34］。为了获取页岩岩心孔隙特征的基础数据，首先对其SEM图像进行孔隙定量分析。对图1采用多次对比原图像的方法获取最佳孔隙选区效果［35］（图2）。针对图1中所包含的有机质范围内的孔隙进行统计分析知，孔隙数目为574个，孔隙直径为0.006～0.255 μm，孔隙面积为0.04×10-3～50.78×10-3μm2，孔隙周长为0.006～2.533 μm，有机质覆盖总面积为8.39 μm2。

图1　页岩样品有机质孔隙Fig.1　Organic pores in the shale samples（red part）

图2　页岩样品有机质区域Fig.2　Selected organic area in the shale samples

3　改进算法重构页岩岩心

应用模拟退火法重构页岩孔隙时，虽然两点概率函数和线性路径函数曲线拟合较好，但是不同尺度的孔隙随机分布于整个重构区域，不能体现页岩真实的有机质孔隙分布特征，而且计算时间长。因而，需要一种既能重构出多尺度的孔隙，又能使孔隙满足页岩的特殊孔隙分布特征的算法。为此，根据页岩岩心孔隙特征的分析以及模拟退火法对页岩岩心的重构，进行了2个方面的改进：①根据岩心中孔隙的多尺度性和分布不均匀性，提出双重区域重构模式；②因重构算法耗时较大，为提高重构速度应用非常快速的模拟退火法（简称VFSA）对岩心进行重构，总结出基于非常快速模拟退火法的页岩岩心双重区域重构方法。

3.1双重区域重构模式

根据SEM图像观察分析知，页岩岩心孔隙分布极不均匀，多集中发育于有机质和部分粘土矿物中，这将影响重构出的页岩岩心的性质，使得重构孔隙随机的分布于整个重构区域，为此提出了双重区域重构模式。该模式将页岩岩心分为2大区域：①孔隙发育程度较好的孔隙发育区，包括有机质和部分孔隙较发育的粘土矿物区域；②孔隙发育程度较低的基质区。

3.2非常快速模拟退火法

非常快速模拟退火法是一种依赖于温度的似Cauchy分布产生新的扰动模型和快速的降温方法。该方法可以提高岩心重构的速度与规模，已广泛应用于地震谱反演方法［24］、地球物理反演的不确定性估计［36］、土地数据同化方法问题［37］等领域。实验研究选择单点概率函数、两点概率函数和线性路径函数3个输入数据的约束函数，并生成随机岩心系统。根据约束函数计算出系统的能量值关系式为

式中：E为系统能量；i为相的序号；αi和βi为对应不同自变量的函数的权重值；Sr为新系统两点概率函数；ri为第i相中某一点；S为原系统两点概率函数；Lr为新系统线性路径函数；L为原系统线性路径函数。

3.2.1约束函数

单点概率函数单点概率函数是描述岩心孔隙空间的最基本特征。实验所建的岩心只考虑岩石和骨架2相系统，因而相函数可简化为

式中：Z（r）为相函数；r为系统中某一点。

岩心系统的单点概率函数可通过统计相函数的平均值得到。

两点概率函数岩心系统中第j相的两点概率函数的定义式为

两点概率函数表征了在岩心系统中随机选取的2点同时分布于同一相中的概率，但其不能反映相的连通性信息。

线性路径函数线性路径函数是描述多孔介质的重要函数，其表征了同一相的连通性。线性路径函数关系式为

其中

式中：

式中：T（k）为迭代过程中每个迭代步对应的T值；k为迭代次数；T0为初始温度；c为给定常数；N为待反演参数的个数；α为温度衰减率，通常为0.7～1.0。为第i相的线性路径函数；rx为连接2点的线段上的任意点；Vj为系统的第j相。

3.2.2扰动方式

传统的模拟退火法中对当前模型进行随机扰动得到新模型，与传统模拟退火法中用到的高斯分布不同，VFSA的扰动方式依赖于温度的似Cauchy分布法，其关系式为

其中

式中：mi′为扰动后的模型中第i个变量；mi是当前模型的第i个变量；yi为变化系数；Bi和Ai为mi取值范围；T为当前系统温度；u为［0，1］内均匀分布的随机数。

应用似Cauchy分布法在高温情况下整个模型范围内进行搜索，低温情况下在当前模型附近进行搜索，易于迅速跳出局部极值，从而加快收敛速度。

3.2.3退火计划

系统温度关系式为

其中

3.3页岩岩心重构方法步骤

电镜图像处理首先对拍摄的页岩岩心的SEM图像进行筛选，对合适的图像进行降噪、阈值分割后得到二值化图像。通过设置合适的阈值，将图像中的孔隙发育区和区域内的孔隙分别提取出来以便分别进行重构。

系统初步优化为确定双重区域的位置，首先选择SA优化系统，设置初始温度，在不断降温的过程中，系统得到初步优化。但单纯使用SA算法重构时间较长，尤其是降到较低温度后，降温的过程异常缓慢甚至难以降温，不能快速重构页岩岩心。

VFSA算法快速建模在重构页岩岩心中的双重区域时，当温度降到某个临界值，重构区域具备了一定的形状，继续在整个重构区域内随机交换待重构页岩岩心中的点，会增加不必要的计算。因此，设置一个较低的温度阈值，应用非常快速模拟退火法优化待重构页岩岩心，这需要根据不同的重构图像尺度和所需的重构精度进行调整，该方法在扰动机制中发挥优势。在低温情况下，似Cauchy分布法会在当前模型附近进行搜索，这样就使整个重构区域的随机点互换转化成了有目的在重构区块附近搜索，迅速的确定区域位置，以最快的速度优化系统，当温度降低到足够低或相邻温度下系统能量变化非常小时，系统优化过程结束，最终获得双重区域的位置。

4　实例分析

非常快速模拟退火法重构获得的页岩岩心与矿场实际页岩岩心性质相似，且重构速度明显提高，该方法更具适用性和可行性。研究区的页岩岩心样本孔隙发育区多集中于有机质中，粘土矿物主要为黄铁矿，且黄铁矿中孔隙发育程度不高，对其具有代表性的扫描电镜图像（图3）进行重构与分析。首先对SEM图像进行处理分析得出，放大倍数为10 000，尺寸为200×200像素点，像素点大小为0.128×0.128 μm，有机质含量为18.56%，孔隙度为1.82%，孔隙数为1 170个，孔隙直径为0.025～10.582 μm，孔隙周长为0.025～1.343 μm。重构运算中，SA初始温度设为10，SA系统降温比为0.95，中间阈值设为10-8，VFSA温度衰减率为0.97，反演参数个数的倒数为0.7，VFSA终止温度为10-20。

重构结果表明，深灰色部分为有机质，黑色部分为孔隙，白色部分为黄铁矿，其中有机质和黄铁矿分区域，与岩心图像有较好的相似性。通过实例数据计算两点概率函数和线性路径函数，并绘制曲线（图3）。从图3可以看出，矿场岩心和重构岩心函数曲线完全重合。

图3　VFSA重构页岩图像函数拟合图像Fig.3　Function fitting results of the reconstructed shale images by VFSA

5　结论

基于非常快速模拟退火法的页岩岩心双重区域重构方法有效地区分出孔隙发育区与基质区，耦合孔隙特征定量分析获得表征参数，应用单点概率函数、两点概率函数和线性路径函数作为约束条件，与传统模拟退火法相比，所获得的重构结果更符合页岩岩心孔隙特征，并缩短了重构时间。目前由于受到建模输入数据的限制，获得数字岩心没有具体包含页岩的天然裂缝，这一问题将在下一阶段的研究工作中加以解决。

致谢：感谢导师隋微波副教授对我在学术上孜孜不倦的教诲以及在生活上的关心，感谢实验室同学的帮助。

参考文献：

［1］ 葛忠伟，樊莉.页岩气研究中应注意的问题［J］.油气地质与采收率，2013，20（6）：19-22. Ge Zhongwei，Fan Li.Some notable problems about shale gas in the scientific research［J］.Petroleum Geology and Recovery Efficiency，2013，20（6）：19-22.

［2］ 苗文培，姜汉桥，葛洪魁，等.页岩气储层蠕变特性及其对页岩气开发的影响［J］.油气地质与采收率，2014，21（4）：97-100. Miao Wenpei，Jiang Hanqiao，Ge Hongkui，et al.Gas shale creep and its influence on the shale gas development［J］.Petroleum Geology and Recovery Efficiency，2014，21（4）：97-100.

［3］ 康毅力，陈强，游利军，等.页岩气藏水相圈闭损害实验研究及控制对策——以四川盆地东部龙马溪组露头页岩为例［J］.油气地质与采收率，2014，21（6）：87-91. Kang Yili，Chen Qiang，You Lijun，et al.Laboratory investigation of water phase trapping damage in shale gas reservoir-a case of Longmaxi shale in the eastern Sichuan Basin［J］.Petroleum Geology and Recovery Efficiency，2014，21（6）：87-91.

当下加拿大天然气行业需要解决基础设施建设、降低成本和跨地区沟通协调等问题。中国石油参与加拿大LNG Canada项目需注意风险防范。

［4］ 李海，白云山，王保忠，等.湘鄂西地区下古生界页岩气保存条件［J］.油气地质与采收率，2014，21（6）：22-25. Li Hai，Bai Yunshan，Wang Baozhong，et al.Preservation conditions research on shale gas in the lower Paleozoic of western Hunan and Hubei area［J］.Petroleum Geology and Recovery Efficiency，2014，21（6）：22-25.

［5］ 徐祖新，郭少斌.页岩气成藏机理研究［J］.内蒙古石油化工，2011，（6）：122-124. Xu Zuxin，Guo Shaobin.Research on accumulation mechanism of shale gas［J］.Inner Mongolia Petrochemical Industry，2011，（6）：122-124.

［6］ 李治平，李智锋.页岩气纳米级孔隙渗流动态特征［J］.天然气工业，2012，32（4）：50-53. Li Zhiping，Li Zhifeng.Dynamic characteristics of shale gas flow in nanoscale pores［J］.Nature Gas Industry，2012，32（4）：50-53.

［7］ Jack Dvorkin，Ayako Kameda，Amos Nur.Real time monitoring of permeability，elastic moduli and strength in sands and shales using digital rock physics［C］.SPE 82246，2003.

［8］ Fredrich J T，Menendez B，Wong T F.Imaging the pore structure of geomaterials［J］.Science，1995，268：276-279.

［9］ Tomutsa L，Silin D，Radmilovic V.Analysis of chalk petrophysical properties by means of submicron-scale pore imaging and modeling［J］.SPE Reservoir Evaluation&Engineering，2007，10（3）：285-293.

［10］Curtis M E，Ambrose R J，Sondergeld C H，et al.Structural characterization of gas shales on the micro-and nano-scales［C］.SPE 137693，2010.

［11］Pal-Eric Oren，Stig Bakke.Process based reconstruction of sandstone and prediction of transport properties［J］.Transport in Porous Media，2002，46（2/3）：311-343.

［12］黄保家，黄合庭，吴国瑄，等.北部湾盆地始新统湖相富有机质页岩特征及成因机制［J］.石油学报，2012，33（1）：25-31. Huang Baojia，Huang Heting，Wu Guoxuan，et al.Geochemical characteristics and formation mechanism of Eocene lacustrine organic-rich shales in the Beibuwan Basin［J］.Acta Petrolei Sinica，2012，33（1）：25-31.

［13］Joshi M.A class of stochastic models for porous media［D］.Kansas：University of Kansas，1974.

［14］Quiblier J A.A new three-dimensional modeling technique for studying porous media［J］.Journal of Colloid and Interface Science，1984，98（1）：84-102.

［15］Ioannidis M，Kwiecien M，Chatzis I.Computer generation and application of 3-D model porous media from pore-level geostatistics to the estimation of formation factor［C］.SPE 30201，1995：185-194.

［16］Hazlett R D.Statistical characterization and stochastic modeling of pore networks in relation to fluid flow［J］.Mathematical Geology，1997，29（6）：801-822.

［17］Yeong C L Y，Torquato S.Reconstructing random media［J］.Physical Review E，1998，57（1）：495-505.

［18］Yeong C L Y，Torquato S.Reconstructing random mediaⅡThreedimensional media from two-dimensional cuts［J］.Physical Review E，1998，58（1）：224-233.

［19］Hidajat I，Rastogi A，Singh M，et al.Transport properties of porous media from thin sections［C］.SPE 69623，2002.

［20］Talukdar M S，Torsaeter O，Ioannidis M A.Stochastic reconstruction of particulate media from two-dimensional images［J］.Journal of Colloid and Interface Science，2002，248（2）：419-428.

［21］Talukdar M S，Torsaeter O，Howard J J.Stochastic reconstruction of chalk sample containing vuggy porosity using a conditional simulated annealing technique［J］.Transport in Porous Media，2004，57（1）：1-15.

［22］ČapekP，Hejtmánek V，Brabec L，et al.Stochastic reconstruction of particulate media using simulated annealing：improving pore connectivity［J］.Transport in Porous Media，2009，76（2）：179-198.

［23］Bakke S，Oren P E.3D pore-scale modelling of sandstones and flow simulations in the pore networks［J］.SPE Journal，1997，2：136-149.

［24］Фren P E，Bakke S.Reconstruction of Berea and stone and porescale modelling of wettability effects［J］.Journal of Petroleum Science and Engineering，2003，39（3）：177-199.

［25］Strebelle S.Sequential simulation drawing structures from training images［D］.Palo Alto：Stanford University，2000.

［26］Strebelle S.Conditional simulation of complex geological structures using multiple-point statistics［J］.Mathematical Geology，2002，34（1）：1-22.

［27］Strebelle S，Journel A G.Reservoir modeling using multiple-point statistics［C］.SPE 71324，2001.

［28］Strebelle S，Payrazyan K，Caers J.Modeling of a deepwater turbidite reservoir conditional to seismic data using principal component analysis and multiple-point geostatistics［C］.SPE 85962，2003.

［29］Okabe H，Blunt M J.Prediction of permeability for porous media reconstructed using multiple-point statistics［J］.Physical Review E，2004，70（6）：066135（1-10）.

［30］Okabe H，Blunt M J.Pore space reconstruction using multiplepoint statistics［J］.Journal of Petroleum Science and Engineering，2005，46（1）：121-137.

［31］张挺，李道伦，卢德唐，等.基于多点地质统计法的多孔介质重构研究［J］.中国科学G辑：物理学力学天文学，2009，39（9）：1 348-1 360. Zhang Ting，Li Daolun，Lu Detang，et al.Research on reconstruction method of porous media using multiple-point geostatistics ［J］.Science in China：Series G Physics，Mechanics&Astronomy，2009，39（9）：1 348-1 360.

［32］Wu Kejian，Naoise Nunan，John W Crawford，et al.An efficient markov chain model for the simulation of heterogeneous soil structure［J］.Soil Science Society of American，2004，68（2）：346-351.

［33］Wu Kejian，Van Dijke，Marinus I J，et al.3D stochastic modelling of heterogeneous porouos media-applications to reservoir rocks ［J］.Transport in Porous Media，2006，65（3）：443-467.

［34］王晨晨，姚军，杨永飞，等.碳酸盐岩双孔隙数字岩心结构特征分析［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.

［35］张季如，祝杰，黄丽，等.土壤微观结构定量分析的IPP图像技术研究［J］.武汉理工大学学报，2008，30（4）：80-83. Zhang Jiru，Zhu Jie，Huang Li，et al.IPP image technique used for quantitative analysis of soil microstructure［J］.Journal of Wuhan University of Technology，2008，30（4）：80-83.

［36］Sen M K，Stoffa P L.Bayesian inference，Gibbs'sampler and uncertainty estimation in geophysical inversion［J］.Geophysical Prospecting，1996，44（2）：313-350.

［37］Xin Li.A very fast simulated re-annealing（VFSA）approach for land data assimilation［J］.Computer&Geosciences，2004，30（3）：239-248.

编辑王星

Reconstruction of shale with dual-region based on very fast simulated annealing algorithm

Zou Mengfei，Sui Weibo，Wang Xudong，Zhang Shuang
（Faculty of Petroleum Engineering，China University of Petroleum（Beijing），Beijing City，102249，China）

The importance of shale reservoir microstructure on field appraisal and flow mechanism has become increasingly prominent.But microscopic flow of shale reservoir cannot be studied effectively by conventional laboratory experiments and approaches for the unique characteristics of extra-low permeability，low porosity，serious anisotropy and heterogeneity.Digital core technology has become an important method in studying reservoir microstructure.Reconstruction methods have the accuracy and reliability.Shale core has been reconstructed by using traditional simulated annealing method based on the real shale SEM 2D images and data extracted by software.In order to solve the problem in reconstruction，an approach based on dual-region strategy and very fast simulated annealing（VFSA）algorithm is proposed in this paper.Using the new algorithm makes a significant improvement of computing efficiency.The shale is reconstructed in dual-region due to the pore distribution characteristics of shale.The results show that real pore distribution features are characterized better by this new approach.

digital core reconstruction；shale；dual region；very fast simulated annealing；simulated annealing method

TE125.1

A文章编号：1009-9603（2015）05-0117-06

2015-07-27。

邹孟飞（1989—），女，黑龙江绥化人，在读硕士研究生，从事渗流机理研究。联系电话：18001229567，E-mail：zoumengfeihf@126. com。

国家自然科学基金青年基金项目“页岩气藏流体运移机理与渗流参数表征方法研究”（51104116）。