2015, Vol. 58
2. 中国石油大学(北京)CNPC物探重点实验室, 北京 102249;
3. 石家庄经济学院, 石家庄 050031
2. CNPC Key Laboratory of Geophysical Prospecting, China University of Petroleum, Beijing 102249, China;
3. Shijiazhuang University of Economics, Shijiazhuang 050031, China
The curvelet transform and projection onto convex sets (POCS) method are used to achieve the reconstruction of 3D digital cores. The curvelet transform is a sparse transform which has been widely used in seismic data denoising and interpolation and image denoising. The POCS method is an efficient method for seismic data interpolation and can be used in the reconstruction of 3D digital cores. This method is applied on each vertical slice and the 3D digital cores can be obtained after all the vertical slices are processed. Besides, the proposed method is superior to the spgl1 method.
With the proposed method, we achieve the 3D digital cores from the cores which were sampled one per two slices in depth for the 3D volume of sandstone obtained by X ray scanner. The reconstruction result is consistent with the original one and superior to the spgl1 method, which proves the validity and superiority of the proposed method. Then the proposed method is applied to the shale cores, of which the 2D horizontal slices are obtained using focused ion beam scanning electron microscopy (FIB-SEM). Because of instrumental limitations, the vertical resolution is almost halved compared with the horizontal resolution. With the proposed method, we can double the horizontal slice in depth, which can help improve the vertical resolution to make it consistent with the horizontal one, weakening the discontinuity of shale slices in depth caused by the instrument limitation, resulting in a clearer fracture distribution.
3D digital cores are reconstructed from the 2D shale slices based on the projection onto convex sets (POCS) method in the curvelet domain. The sand reconstruction test of 3D digital cores demonstrates that the proposed method is more suitable for rock slice reconstruction with high efficiency and accuracy compared with the popular spgl1 method. The shale reconstruction test of 3D digital cores from nano-scale shale slices obtained by FIB-SEM indicates that vertical spatial distribution of fractures is more clear and the minimum interval between adjacent vertical slices is basically consistent with the horizontal resolution after reconstruction, which can lay the foundation for the subsequent data processing and related simulation analysis of shale digital cores. 3D digital cores reconstruction numerical tests on 2D sand and shale slices demonstrate the validity of the proposed method. The methods based on direct 3D datasets and efficient sparse transform will be developed in future work.
随着勘探程度的日益加深，勘探难度逐渐增加.近几年页岩气的研究成为非常规能源的研究热点，页岩气可成为日益减少的石油能源的替代能源之一.但是页岩气的勘探与开发具有很大的挑战性，页岩气主要富集于页岩之中，而页岩的孔隙性低、渗透率低成为页岩气开发的主要瓶颈.因此研究页岩裂隙的精细空间展布，是勘探开发页岩气十分重要的环节.三维数字岩心技术可在孔隙尺度上描述岩石的微观结构，为精细研究页岩裂隙的空间展布提供技术支持.姚军等(2005)对数字岩心的现状作了总结，指出现有数字岩心重建方法的不足，包括过程法无法对复杂的沉积体系进行模拟以及随机法建立的数字岩心不具有大范围的传导性等缺点；还探讨了数字岩心技术在微观渗流研究、岩心驱替模拟实验以及评价驱油剂效果等方面的应用.常用的三维数字岩心重建方法有两种：X射线层析扫描成像(CT)和基于岩石二维切片的三维重建技术，基于二维切片的三维重建技术又分为过程法以及随机法(ren and Bakke，2002; Okabe et al.，2004;Liu et al.，2009a，b; 刘学锋等，2013).刘学锋等(2013)基于砂岩分析比较了不同方法重建三维数字岩心的效果，并对三维数字岩心的发展进行了展望，指出X射线扫描方法精度的局限性以及聚焦离子束扫描电镜纳米级高精度的优势；同时，将孔隙和骨架进行二值化处理，虽然重建出的三维数字岩心孔隙展布的分辨率较高，但是忽略了裂隙的强度，会对后续的解释产生影响.X射线方法可以得到三维数字岩心，但其精度较低，对裂隙的空间展布表征较粗糙，难以满足页岩纳米-微米级孔隙成像的要求.而高精度的聚焦离子束扫描电镜可以实现页岩纳米级孔隙的精确成像，但该设备在保证岩石纵向上成像最小间隔与横向分辨率一致方面存在困难，影响了页岩数字岩心数据处理及其相关模拟分析精度的提高.页岩气勘探与有效开发的深入迫切需要对页岩裂隙的纵向展布进行精确成像，因此，对聚焦离子束扫描电镜扫描获得二维页岩切片进行高精度三维重建，显得尤为重要.
页岩孔隙结构复杂，传统的线性数据重建方法很难满足其三维重建精度要求.Curvelet变换是一种新兴的稀疏变换，属于多尺度数学变换的范畴；相对于小波变换，其增加了方向性，更适合表征曲线奇异性，已被应用于图像去噪(Candès，2001; Starck et al.，2002)、图像重构(Starck et al.，2001)、地震数据去噪(Neelamani et al.，2008)、地震数据插值(Hennenfent and Herrmann，2008; Herrmann and Hennenfent，2008; Naghizadeh and Sacchi，2010; 曹静杰等，2012)以及多次波的去除(Herrmann et al.，2007)等方面，体现出其巨大的优势.凸集投影方法(POCS)首先由Bregman(1965，1967)提出，并被用于图像重建(Stark and Oskoui，1989; Youla and Webb，1982)；Abma和Kabir(2006)将POCS方法引入非规则地震数据插值方面，此后，若干学者在计算量和阈值选取等方面对该方法进行了改进，大幅提高了收敛速度和计算精度(Gao et al.，2010，2012; Yang et al.，2012; 张华和陈小宏，2013).
针对页岩孔隙结构的复杂性，研究将POCS方法与curvelet稀疏变换有效结合，充分发挥两者的优势，对高精度聚焦离子束扫描电镜获得的纳米级二维页岩切片进行三维图像重建，得到三维数字岩心，分析表明其裂隙的纵向空间展布更连续，可为后续页岩数字岩心数据处理及其相关模拟分析精度提高奠定基础.为了验证新方法的有效性，研究首先对X射线扫描砂岩得到的三维数字岩心进行处理试验，分析比较重建前后的三维岩心数据体的异同，验证新方法的有效性；其次利用本文方法对聚焦离子束扫描电镜扫描得到的纳米级页岩切片进行三维重建，使得裂隙的纵向空间展布更连续.砂岩以及实际页岩的数据处理均验证了本文方法的有效性.
2 方法原理 2.1 Curvelet变换Curvelet变换相对小波变换增加了角度的参数，使得其更适合于曲线奇异性的表征，适合页岩的复杂孔隙结构特征.Curvelet变换经历了ridgelet变换、第一代curvelet变换以及第二代curvelet变换，现在常用的curvelet变换为第二代curvelet变换(Candès et al.，2006)，其表达式为
由于仪器设备的限制，聚焦离子束扫描电镜获得的纳米级二维页岩切片在深度上具有不连续性，两片之间存在一定的间隔.为了更好地满足页岩数字岩心数据处理及其相关模拟分析精度要求，需要将其考虑成图片重建问题，对三维数据体的每一纵切片进行处理，计算公式为
方程(2)的求解是不适定的，考虑到 d 0在curvelet域的稀疏性，采用稀疏促进策略，构建目标泛函为
Blumensath和Davies(2008，2009)以及Loris(2010)指出，当P(x)=‖ x ‖1时，Tλk满足条件为
 | 图 1 结合POCS方法与Curvelet稀疏变换三维数字岩心重建流程图 Fig. 1 Flow chart of the 3D digital core reconstruction with the POCS method in the curvelet domain |
首先基于X射线扫描产生的砂岩三维数据体进行处理试验与分析，证明本文方法的有效性.先对X射线扫描产生的砂岩三维数据体在深度上进行隔片抽稀，利用新方法以及常用的spgl1方法(Van Den Berg and Friedl and er，2008)对每一纵切片进行重建，将重建后的纵切片与原始三维数据体中的纵切片进行比较，对比不同方法重建结果的优劣，分析新方法的优势以及有效性.其次利用新方法对聚焦离子束扫描电镜扫描得到纳米级的二维实际页岩切片进行三维重建，对比分析三维重建前后裂隙纵向连续性以及空间展布，证明方法的有效性.
3.1 砂岩由于页岩孔隙通常在纳米-微米级，目前X射线扫描方法的精度还无法满足精细刻画页岩裂隙空间展布的要求.但X射线扫描方法可以对砂岩进行全方位扫描，得到砂岩全三维数据体.由X射线方法扫描得到的砂岩的横、纵切片如图 2所示，可以观测到砂岩颗粒以及裂隙的空间展布，其分辨率为13.4923 μm，即采样点之间距离为13.4923 μm，精度较低.图 2中砂岩样品横切片的实际采样点数为500×500，实际尺寸为6.74615×6.74615 mm，纵切片的实际采样点数为500×100，实际尺寸为6.74615×1.34923 mm.
 | 图 2 砂岩的横、纵切片 Fig. 2 Horizontal and vertical slices of the sandstone |
为验证本文重建方法的有效性，首先将砂岩三维数据体在深度上进行隔片抽稀，抽稀后的横、纵切片如图 3所示.对抽稀后的岩石纵切片分别利用本文研究方法以及spgl1方法进行重建，最大迭代次数设置为50次，重建后砂岩的横、纵切片如图 4a以及4b所示.将砂岩原始、隔片抽稀以及重建后的岩石纵切片进行放大显示，如图 5所示，其中图 5a表示砂岩原始纵切片，图 5b表示深度上抽稀的纵切片，图 5c表示本文方法重建得到的纵切片，图 5d表示spgl1方法重建得到的纵切片.
 | 图 3 抽稀后的砂岩横纵、纵切片 Fig. 3 Horizontal and vertical slices of the sandstone sampled in depth |
 | 图 4 重建得到的砂岩横、纵切片
(a) 本文方法； (b) spgl1方法. Fig. 4 Horizontal and vertical sandstone slices after reconstruction (a) Using the proposed method； (b) Using the spgl1 method. |
 | 图 5 砂岩纵切片的放大显示 Fig. 5 Enlarged vertical slice of the sandstone |
观察完整三维数据体图 2、纵切片重建结果图 4、以及纵切片的放大显示图 5可以看出，本文方法以及spgl1方法均能实现对抽稀后的砂岩进行重建，但本文方法的精度更高；另外每次迭代中，本文方法仅做正反变换以及阈值处理，而spgl1方法需要做正反变换、计算更新步长等运算，因此，本文方法计算效率更高.抽稀后砂岩三维数据重建，验证了本文方法的有效性，且与spgl1方法的效果进行对比，体现了本文方法的优势，为后续将其应用于具有纳米-微米级复杂孔隙结构页岩三维数字岩心重建奠定了基础.
3.2 页岩三维数字岩心重建聚焦离子束扫描电镜可以实现页岩纳米级孔隙的精确成像，但由于设备的局限，页岩切片在深度上具有不连续性，即设备在保证岩石纵向上成像最小间隔与横向分辨率一致方面存在困难，限制了裂隙纵向空间展布的表征，影响了页岩数字岩心数据处理及其相关模拟分析精度的提高，因此利用二维页岩切片进行三维图像重建获得三维数字岩心，成为数字岩石物理研究的最重要步骤之一.研究基于聚焦离子束扫描电镜得到的192个页岩切片进行三维图像重建获得三维数字岩心，其中纵切片间的距离为10 nm，横切片分辨率为5.56 nm，垂直方向和水平方向像素数分别为2048与1768.页岩岩石横、纵切片如图 6所示.
 | 图 6 原始页岩切片组成三维数据体的横、纵切片 Fig. 6 Horizontal and vertical slices of the 3D dataset composed of raw shale slices |
为了更好表征页岩裂隙空间展布，提高页岩数字岩心数据处理及其相关模拟分析精度，研究利用本文方法对原始页岩切片在深度上进行加密处理，最大值迭代次数设置为50次，形成纵向上成像最小间隔与横向分辨率基本一致的三维数字岩心数据，页岩切片是原始切片数的两倍，其纵切片如图 7与图 8所示.对比分析图 7和图 8可知，重建得到的纵切片中，切片之间的距离减半，纵向连续性得到保持，裂隙展布更加清晰，由于仪器限制引起的页岩切片深度上的不连续性得到减弱.
 | 图 7 页岩深度上加密前后的纵切片
(a) 页岩加密前的纵切面； (b) 本文方法重建后页岩的纵切片. Fig. 7 Vertical shale slice before and after reconstruction (a) Before reconstruction; (b) After reconstruction by the proposed method |
 | 图 8 另一方向页岩深度上加密前后的纵切片 (a) 页岩加密前的纵切面； (b)本文方法重建后页岩的纵切片. Fig. 8 Same asFig.7 but in another direction |
将重建前后的岩石切片进行局部三维显示，如图 9所示，重建后裂隙的连续性更好，岩石纵向上成像最小间隔与横向分辨率基本一致，且重建过程中没有引入虚假裂隙，证明了本文方法利用页岩切片实现高精度三维数字岩心重建的有效性.
 | 图 9 页岩切片局部三维显示 (a) 重建前； (b) 重建后. Fig. 9 Local display of the 3D shale dataset (a) Before reconstruction; (b) After reconstruction. |
基于curvelet变换与凸集投影(POCS)方法利用岩石切片实现了三维数字岩心重建.对砂岩的三维数字岩心重建试验表明，与目前广泛应用的spgl1方法相比较，本文方法能更好对岩石切片进行了有效的重建，且耗时较少、精度更高.对聚焦离子束扫描电镜扫描得到的纳米级页岩切片进行三维数字岩心重建试验表明，重建后裂隙的纵向展布更加清晰，保证了数字岩心纵向上成像最小间隔与横向分辨率的一致性，为后期页岩数字岩心数据处理及其相关模拟分析精度提高奠定基础.基于砂岩岩石切片和纳米级页岩岩石二维切片的三维数字岩心重建与效果分析均验证了本文方法的有效性.发展更高效的稀疏变换及直接对三维数据体进行处理的方法将是下一步研究工作的重点.
| [1] | Abma R, Kabir N. 2006. 3D interpolation of irregular data with a POCS algorithm. Geophysics, 71(6): E91-E97, doi: 10.1190/1.2356088. |
| [2] | Blumensath T, Davies M E. 2008. Iterative thresholding for sparse approximations. Journal of Fourier Analysis and Applications, 14(5-6): 629-654, doi: 10.1007/s00041-008-9035-z. |
| [3] | Blumensath T, Davies M E. 2009. Iterative hard thresholding for compressed sensing. Applied and Computational Harmonic Analysis, 27(3): 265-274, doi: 10.1016/j.acha.2009.04.002. |
| [4] | Bregman L. 1965. The method of successive projection for finding a common point of convex sets (Theorems for determining common point of convex sets by method of successive projection). Soviet Mathematics, 6: 688-692. |
| [5] | Bregman L. 1967. The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming. USSR Computational Mathematics and Mathematical Physics, 7(3): 200-217. |
| [6] | Candès E. 2001. The curvelet transform for image denoising. Proceedings 2001 International Conference on Image Processing. Thessaloniki: IEEE. |
| [7] | Candès E, Demanet L, Donoho D, et al. 2006. Fast discrete curvelet transforms. Multiscale Modeling & Simulation, 5(3): 861-899. |
| [8] | Cao J, Wang Y F, Yang C C. 2012. Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization. Chinese J. Geophys. (in Chinese), 55(2): 596-607, doi: 10.6038/j.issn.0001-5733.2012.02.022. |
| [9] | Gao J J, Chen X H, Li J Y, et al. 2010. Irregular seismic data reconstruction based on exponential threshold model of POCS method. Applied Geophysics, 7(3): 229-238, doi: 10.1007/s11770-010-0246-5. |
| [10] | Gao J J, Stanton A, Naghizadeh M, et al. 2012. Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction. Geophysical Prospecting, 61(Suppl. 1): 138-151, doi: 10.1111/j.1365-2478.2012.01103.x. |
| [11] | Hennenfent G, Herrmann F J. 2008. Simply denoise: Wavefield reconstruction via jittered undersampling. Geophysics, 73(3): V19-V28, doi: 10.1190/1.2841038. |
| [12] | Herrmann F J, Baniger U, Verschuur D J E. 2007. Non-linear primary-multiple separation with directional curvelet frames. Geophysical Journal International, 170(2): 781-799, doi: 10.1111/j.1365-246X.2007.03360.x. |
| [13] | Herrmann F J, Hennenfent G. 2008. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International, 173(1): 233-248, doi: 10.1111/j.1365-246X.2007.03698.x. |
| [14] | Liu X F, Sun J M, Wang H T. 2009a. Reconstruction of 3-D digital cores using a hybrid method. Applied Geophysics, 6(2): 105-112, doi: 10.1007/s11770-009-0017-y. |
| [15] | Liu X F, Sun J M, Wang H T. 2009b. Numerical simulation of rock electrical properties based on digital cores. Applied Geophysics, 6(1): 1-7, doi: 10.1007/s11770-009-0001-6. |
| [16] | Liu X F, Zhang W W, Sun J M. 2013. Methods of constructing 3-D digital cores: A review. Progress in Geophys. (in Chinese), 28(6): 3066-3072, doi: 10.6038/pg20130630. |
| [17] | Loris I, Douma H, Nolet G, et al. 2010. Nonlinear regularization techniques for seismic tomography. Journal of Computational Physics, 229(3): 890-905, doi: 10.1016/j.jcp.2009.10.020. |
| [18] | Naghizadeh M, Sacchi M D. 2010. Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data. Geophysics, 75(6): WB189-WB202, doi: 10.1190/1.3509468. |
| [19] | Neelamani R, Baumstein A, Gillard D, et al. 2008. Coherent and random noise attenuation using the curvelet transform. The Leading Edge, 27(2): 240-248, doi: doi: 10.1190/1.2840373. |
| [20] | Okabe H, Blunt M J. 2004. Prediction of permeability for porous media reconstructed using multiple-point statistics. Physical Review E, 70(6): 066135, doi: 10.1103/PhysRevE.70.066135. |
| [21] | Øren P E, Bakke S. 2002. Process based reconstruction of sandstones and prediction of transport properties. Transport in Porous Media, 46(2-3): 311-343. |
| [22] | Starck J L, Donoho D L, Candes E J. 2001. Very high quality image restoration by combining wavelets and curvelets. //Proc. SPIE 4478, Wavelets: Applications in Signal and Image Processing IX, 9, doi:10.1117/12.449693 . |
| [23] | Starck J L, Candès E J, Donoho D L. 2002. The curvelet transform for image denoising. IEEE Transactions on Image Processing, 11(6): 670-684. |
| [24] | Stark H, Oskoui P. 1989. High-resolution image recovery from image-plane arrays, using convex projections. JOSA A, 6(11): 1715-1726. |
| [25] | Van Den Berg E, Friedlander M P. 2008. Probing the Pareto frontier for basis pursuit solutions. SIAM Journal on Scientific Computing, 31(2): 890-912. |
| [26] | Yang P L, Gao J H, Chen W C. 2012. Curvelet-based POCS interpolation of nonuniformly sampled seismic records. Journal of Applied Geophysics, 79: 90-99, doi: 10.1016/j.jappgeo.2011.12.004. |
| [27] | Yao J, Zhao X C, Yi Y J, et al. 2005. The current situation and prospect on digital core technology. Petroleum Geology and Recovery Efficiency (in Chinese), 12(6): 52-54. |
| [28] | Youla D C, Webb H. 1982. Image Restoration by the Method of Convex Projections: Part 1-Theory. IEEE Transactions on Medical Imaging, 1(2): 81-94. |
| [29] | Zhang H, Chen X H. 2013. Seismic data reconstruction based on jittered sampling and curvelet transform. Chinese J. Geophys. (in Chinese), 56(5): 1637-1649, doi: 10.6038/cjg20130521. |
| [30] | 曹静杰, 王彦飞, 杨长春. 2012. 地震数据压缩重构的正则化与零范数稀疏最优化方法. 地球物理学报, 55(2): 596-607, doi: 10.6038/j.issn.0001-5733.2012.02.022. |
| [31] | 刘学锋, 张伟伟, 孙建孟. 2013. 三维数字岩心建模方法综述. 地球物理学进展, 28(6): 3066-3072, doi: 10.6038/pg20130630. |
| [32] | 姚军, 赵秀才, 衣艳静等. 2005. 数字岩心技术现状及展望. 油气地质与采收率, 12(6): 52-54. |
| [33] | 张华, 陈小宏. 2013. 基于jitter采样和曲波变换的三维地震数据重建. 地球物理学报, 56(5): 1637-1649, doi: 10.6038/cjg20130521. |