2017, Vol. 32
2. 页岩油气富集机理与有效开发国家重点实验室, 北京 100083
3. 中国石化页岩油气勘探开发重点实验室, 北京 100083
4. 中国石化石油勘探开发研究院, 北京 100083
2. State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development, Beijing 100083, China
3. SinoPEC Key Laboratory of Shale Oil/Gas Exploration and Production Technology, Beijing 100083, China
4. SinoPEC Petroleum Exploration and Production Research Institute, Beijing 100083, China
地震横波速度在叠前地震属性分析、叠前地震反演、测井岩性解释及多波资料解释等过程中是重要的参数.由于横波测井成本较高，解释难度大等因素，在很多的实际生产资料中都缺少横波测井信息.
许多学者通过建立岩石物理测量、试验和分析建立多种经验关系及理论模型来进行预测横波速度，解决实际工区中常规测井资料缺失横波速度的情况.Castagna等 (1985)提出了“泥岩线”公式——含水碎屑硅酸盐中纵、横波速度之间的经验关系；Han等 (1986)通过不同压力下的岩样测量并考虑了砂岩中孔隙度和粘土含量的影响，给出了纵、横波速度之间的线性回归方程；Greenberg和Castagna (1992)基于Biot-Gassmann理论，假设了纵、横波速度的经验关系和岩石固体组分的线性混合定律给出了一种半物理方法估算横波速度；Xu和White (1996)基于Kuster-Toksöz理论，考虑了砂岩孔隙和泥岩孔隙几何形态，提出了适用于砂泥岩储层的等效介质模型——Xu-White模型；Keys和Xu (2002)引入了干岩石近似对Xu-White模型进行了改进；Lee (2006)通过改进Pride等提出的方法，通过计算压实砂岩干岩石的固结参数估算横波速度；Xu和Payne (2009)利用Xu-White模型对碳酸盐岩估算了孔隙形状，成功把该模型引入到横波速度的预测中，同时对孔隙形状与Gassmann流体替代的问题进行了讨论；Ruiz和Dvorkin (2010)研究了包裹体模型在非碎屑岩中估算弹性性质的应用，并且对横波速度的估算进行了数值实验.上述反演过程使用基于遍历搜索的传统方法，用时长、精度较低，而一些研究将智能优化算法引入反演过程.Li等 (2013)基于DEM解析模型的研究，提出了基于非线性模拟退火法岩石的等效孔隙纵横比反演及应用，并做出了可行性分析.同年刘欣欣等 (2013)等提出了基于自适应遗传算法反演岩石矿物组分的体积、剪切模量从而计算出碳酸盐岩的岩性参数的方法.刘倩等 (2016)提出了一种拓展的岩石物理逆建模型方法，对描述储层特征的参数进行了定量预测.同年罗水亮等将Pride模型与Lee模型进行变形，然后借鉴了地震反演的思路，将目标函数最优化问题转换成为线性矩阵方程迭代求解问题，对横波速度进行了精细预测.
粒子群算法 (Particle Swarm Optimization, PSO) 是由Kennedy和Eberhart于1995年根据鸟群的社会行为提出的一种群集智能优化算法.该算法因其本身具有实现简单，参数设置少和收敛速度快等特点，在很多领域都得到了应用，也能应用于岩石物理孔隙结构反演及横波速度的预测中.本文就将会采用这种智能优化算法进行岩石孔隙结构反演与横波速度的预测.
上述所提出的经验公式与岩石物理理论模型均是针对于常规碎屑岩、碳酸盐岩地层所提出的，对于页低孔低渗岩储层并不适用，而针对页岩储层的岩石物理建模与横波速度预测研究较少.本文首先通过建立页岩岩石物理模型，计算出岩石纵、横波速度与密度、孔隙度和矿物组分等参数之间的关系，其次利用传统遍历搜索方法、粒子群算法两种方法计算出一个最佳的孔隙纵横比使得理论纵波速度与实际纵波速度的误差最小，利用孔隙纵横比作为预测横波速度的约束条件进行横波速度的预测，最后利用已知的横波速度作为对比条件，进行两种方法的误差分析，同时对比两种方法的准确性.
1 方法原理 1.1 自相容近似岩石物理模型等效介质自相容近似模型 (Self-Consistent Approximation, SCA) 是由Berryman (1980)年提出来的.其中心思想就是用暂时未知的等效介质作为背景介质，其中模型中包含物之间的弹性互动会受到固体颗粒形状和孔隙形状的影响，通过调节背景介质的弹性参数，令多相弹性介质参数与背景介质的弹性参数相匹配，这时所得到的背景介质的弹性模量即为多相弹性介质的等效弹性模量.
这种方法给出的n相矿物和孔隙空间的自相容弹性模量计算公式为
|
(1) |
|
(2) |
其中，每个j代表一种矿物相或孔隙空间，对应每个相的体积分数xj、体积 (Kj) 和剪切 (μj) 模量.
(1) 式和 (2) 式中P和Q的上标*j指的是此几何因数是针对具有自相容等效模量KSC*和μSC*的背景介质中包含物材料j.纵横比α≤1的椭球体包含物的P和Q值为
|
(3) |
|
(4) |
上式中，张量Tijkl将均匀远场应变场与椭球包含物内的应变联系起来.Berryman给出了计算P和Q所需的相关的标量为
|
(5) |
|
(6) |
其中：
|
(7) |
上式中的A, B, R, θ和f分别为
|
(8) |
在P和Q的公式中，所有的Km和μm分别用KSC*和替代μSC*，通过同时迭代的方法来求解耦合的公式可以得到SCA模型下岩石的等效体积模量、剪切模量和密度.
1.2 基于岩石物理模型的横波速度预测原理前人研究表明，随着孔隙纵横比的增大，孔隙的刚度也会随之增加，当岩石的总刚度增加会导致地震波速度的增大.这说明了，孔隙纵横比会影响岩石的弹性模量，从而影响岩石的纵、横波速度，所以在横波速度的预测中，我们需要利用岩石的孔隙纵横比作为约束条件，通过已知的纵波速度预测未知的横波速度.
设定岩石物理模型为低孔低渗情况下的SCA模型，模型计算了矿物组分和孔隙空间的等效弹性模量，其中矿物组分由实际测井资料给出，孔隙空间混合流体的密度ρfl和体积模量Kfl用Reuss平均理论计算，公式为
|
(9) |
式中，xi、Ki和ρi分别是各孔隙流体的体积分数、体积模量和密度.
由SCA计算的等效弹性模量计算饱和岩石的纵、横波速度，公式为
|
(10) |
K、μ和ρ分别为饱和岩石的体积模量、剪切模量和密度.
(1)-(2) 式和 (9)-(10) 式给出了纵、横波速度与孔隙纵横比之间的定量关系.确定了纵、横波速度与孔隙纵横比之间的定量关系，将目标函数设定为min (Vp观测值-Vp计算值)，寻找符合该目标函数的孔隙纵横比的值，即为反演的过程.
对于该问题传统反演方法基于遍历搜索，步骤如下：
1) 给出孔隙纵横比的初始数值.设置孔隙纵横比的取值范围为0.001到1，综合考虑计算精度和计算效率将计算增量设为0.005.
2) 输入该深度矿物组分、孔隙度、密度等测井资料，基于岩石物理模型计算各个孔隙纵横比取值时页岩的弹性参数进而计算出纵波速度.
3) 输入该深度纵波速度测井数据，寻找min (Vp观测值-Vp计算值) 时的孔隙纵横比值，由该值计算得出预测的横波速度.
4) 依次对每个测井深度进行上述操作.
1.3 基于粒子群算法的横波速度预测粒子群算法是由Kennedy和Eberhart (1995)所提出的一种新的智能优化算法.在“Particle Swarm Optimization”的文章中解释了这种新的优化算法的思想是基于鸟类觅食行为的研究所得出的.粒子群算法的计算方式是先在可行解空间中随机初始化一群粒子，这群粒子中的每一个都是优化问题的一个可行解，再由目标函数提供出一个确定的适应度值，这个适应度值为1/(Vp观测值-Vp计算值), 同时还有一个速度来解决其方向和距离.基础粒子会跟随最优粒子运动，通过逐代搜索的方式得到最优解.在每一代中，粒子通过跟踪两个“极值”来更新自己，一个为粒子本身目前找到的最优解Pbest，另一个为全种群目前找到的最优解Gbest.
找到这两个最优值后, 粒子根据如下公式来更新自己的速度和位置，公式为
|
(11) |
|
(12) |
上式中，下标i表示第i个粒子；下标m表示粒子的第m维；ω为惯性权值；常数c1和c2为学习因子；r1和r2是取值范围为[0，1]之间的随机数；vimk和vimk+1分别表示第k次和第k+1次迭代时第i个粒子的速度矢量的第m维分量；和分别表示第k次和第k+1次迭代时第个粒子的位置矢量的第m维分量；pim是第i个粒子的个体最佳位置Pbest的第m维分量；pgm是种群最好位置Gbest的第m维分量.
同样，对于目标函数min (Vp观测值-Vp计算值)，寻找符合该目标函数的孔隙纵横比值.常规方法基于遍历搜索，计算速度缓慢，应用上述智能优化算法可使计算速度更快，结果更为准确.图 1给出了采用粒子群算法后根据纵波速度反演孔隙纵横比预测横波速度的流程.其中岩石物理模型为低孔低渗情况下的SCA模型，模型计算了矿物组分和孔隙空间的等效弹性模量.
|
图 1 采用粒子群算法的反演孔隙纵横比预测横波速度流程图 Figure 1 The inversion of porosity by using PSO algorithms flow chart of aspect ratio to predict shear wave velocity |
为了检验本文所采用的横波速度预测方法的效果，选用了所研究区域带有实测横波速度资料的页岩储层井位，但该井位对其进行了实际计算，同时对预测结果和实测数据进行了误差分析.研究区域所选用的井位，其目的层段位于四川盆地下志留统龙马溪组的底，岩性主要为页岩、含粉砂质页岩和碳质页岩.富有机质页岩厚度介于100~150 m，底部40多米厚优质页岩分布稳定，有机碳含量高，现场测试含气量高.
图 2中给出了该井位地层测井曲线，(从左到右依次为) 伽马射线 (GR)、纵波速度 (VP)、横波速度 (Vs)、密度 (ρ) 和孔隙度 (Φ).从图中可以看到，纵、横波速度曲线的相关性很强，与伽马射线测井曲线呈负相关变化规律.图中标出了龙马溪组页岩地层，位于含气储层位于页岩层底部.与周围地层相比，页岩层具有GR 值较大，密度值较小的特征.
|
图 2 研究区域井位测井曲线 Figure 2 Study area of well logging curves |
图 3给出了测井深度范围内地层各矿物组分，该研究区域内的页岩矿物主要有石英和黄铁矿，含少量伊利石，储层内的孔隙流体主要是水和气.岩石物理模型中用到的矿物和流体的弹性性质在表 1中给出.
|
图 3 研究区域井位地层矿物组分 Figure 3 Study area well formation mineral composition |
|
表 1 矿物和孔隙流体经验弹性参数 Table 1 Mineral and pore fluid experience elastic parameters |
通过上述的实际测井资料与矿物和孔隙立体经验弹性参数，进行实际计算.图 4为该研究区域井位孔隙纵横比和横波速度计算结果.这里使用的横波速度预测采用的是传统方法.从图中可以看出，页岩部分的孔隙纵横比数值比较稳定，但围岩部分的孔隙纵横比变化范围较大.从中可以得出，页岩的孔隙结构与围岩相比更为稳定.在页岩层范围内，比较预测的横波速度曲线 (红色) 与实测的横波速度曲线 (黑色) 吻合较好，两者的变化趋势一致.估算的横波速度与实测的横波速度数值误差较小.
|
图 4 采用传统算法的研究区域井位孔隙纵横比和横波速度计算结果 (红色为传统遍历搜索预测横波速度曲线，黑色为实测的横波速度曲线) Figure 4 Using traditional algorithm study area location pore aspect ratio and shear wave speed calculation results (Red for the traditional ergodic search to predict shear velocity curve, black for the measured shear velocity curve) |
图 5给出了研究区域井位的孔隙纵横比频数分布直方图，它的形态接近于正态分布.
|
图 5 孔隙纵横比频数分布直方图 Figure 5 Pore aspect ratio frequency distribution histogram |
为了对比采用了粒子群算法的反演孔隙纵横比预测横波速度与传统方法的精度，仍然选择上述研究区域井位的实际测井资料，并对其实际计算.图 6为采用了粒子群算法的反演孔隙纵横比预测横波速度的计算结果，从图中可以看出，研究区域井位页岩部分的孔隙纵横比值较为稳定，而围岩部分的纵横比变化范围较大，相比于传统方法，孔隙纵横比在图像上面更加稳定.在页岩范围内，预测的横波速度曲线 (红色) 与实测的横波速度曲线 (黑色) 吻合程度较好，预测的横波速度与实测的横波速度数值误差较小.
|
图 6 采用PSO算法的研究区域井位孔隙纵横比和横波速度计算结果 (红色为粒子群算法预测横波速度曲线，黑色为实测的横波速度曲线) Figure 6 Using PSO algorithm study area location pore aspect ratio and shear wave speed calculation results (Red for predicting shear velocity curve of PSO algorithm, black for the measured shear velocity curve) |
为了便于对比两种方法的横波速度预测结果，将两种方法预测的横波速度结果与实际资料的横波速度放置一起，针对我们所研究的页岩部分进行截取，截取范围为深度2320~2420 m.
如图 7所示，对比图中黑色是实测横波速度曲线，蓝色是传统遍历搜索方法预测横波速度曲线，红色是粒子群算法预测横波速度曲线.
|
图 7 采用遍历搜索与粒子群算法的横波速度预测结果曲线对比图 (黑色是实测横波速度曲线，蓝色是传统遍历搜索方法预测横波速度曲线，红色是粒子群算法预测横波速度曲线) Figure 7 Using traversal search and shear velocity of PSO algorithm to predict results figure curve (Black is the measured shear velocity curve, blue is a traditional traversal search method to predict shear velocity curve, red is PSO algorithm to predict shear velocity curve) |
通过对比采用粒子群算法与传统方法对研究区域井位孔隙纵横比和横波速度计算的结果，使用传统横波速度预测方法的预测与实际差值均值与方差为:0.1452、0.0229，误差均方根为:0.2098，而使用基于粒子群算法的横波速度预测方法的预测与实际差值均值与方差为0.1352、0.0219，误差均方根为0.2068.对比图 4、6，可以看出，采用了粒子群算法计算出的孔隙纵横比，无论是在围岩部分还是在页岩部分都要比采用传统方法得出的结果更加的稳定.在预测与实际差值均值与方差、误差均方根上面采用粒子群算法的计算结果比采用传统算法的计算结果，准确度有所提高，更接近于实际横波速度的真实值.
3 结论 3.1本文提出了基于粒子群算法的页岩孔隙纵横比和横波速度预测方法.利用岩石物理建模得出孔隙纵横比与横、纵波速度，岩石孔隙度与矿物组分等参数之间的定量关系，利用孔隙纵横比作为横波速度预测的约束条件，并采用传统算法与基于粒子群算法进行了预测横波速度.
3.2通过对实际井位数据的反演结果证明了基于粒子群算法估算横波速度的有效性，同时对比了采用传统算法估算横波速度的结果，使用基于粒子群算法的横波速度预测方法的预测与实际差值均值与方差、误差均方根都优于采用传统遍历搜索方法计算的预测结果，证明了基于粒子群算法估算横波速度的准确性.对孔隙纵横比与横波速度预测值得误差分析表明，计算得到的横波速度与实测值的误差较小，能很好的拟合出横波速度的变化曲线.
3.3由于模型建立基于各向同性假设，存在理想化近似，不可能完全符合实际岩石的多样性和复杂性，加之围岩部分由于孔隙度的影响，反演结果不佳.但加入了粒子群算法后，相比传统算法有所提高.也再次证明了加入粒子群算法后，对于整体的实际井位数据反演结果的准确性有所提高.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [] | Berryman J G. 1980. Long-wavelength propagation in composite elastic media Ⅰ. Spherical inclusions[J]. The Journal of the Acoustical Society of America, 68(6): 1809–1831. DOI:10.1121/1.385171 |
| [] | Castagna J P, Batzle M L, Eastwood R L. 1985. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks[J]. Geophysics, 50(4): 571–581. DOI:10.1190/1.1441933 |
| [] | Greenberg M L, Castagna J P. 1992. Shear-wave velocity estimation in porous rocks:Theoretical formulation, preliminary verification and applications[J]. Geophysical Prospecting, 40(2): 195–209. DOI:10.1111/gpr.1992.40.issue-2 |
| [] | Han D H, Nur A, Morgan D. 1986. Effects of porosity and clay content on wave velocities in sandstones[J]. Geophysics, 51(11): 2093–2107. DOI:10.1190/1.1442062 |
| [] | Kennedy J, Eberhart R C. 1995. Particle swarm optimization[C].//Proceedings of the Institute of Electrical and Electronics Engineers Conference on Neural Networks, Ⅳ. Piscataway, NJ, 1943-1948. |
| [] | Keys R G, Xu S Y. 2002. An approximation for the Xu-White velocity model[J]. Geophysics, 67(5): 1406–1414. DOI:10.1190/1.1512786 |
| [] | Lee M W. 2006. A simple method of predicting S-wave velocity[J]. Geophysics, 71(6): F161–F164. DOI:10.1190/1.2357833 |
| [] | Li H B, Zhang J J, Yao F C. 2013. Inversion of effective pore aspect ratios for porous rocks and its applications[J]. Chinese Journal of Geophysics, 56(1): 43–51. DOI:10.1002/cjg2.20004 |
| [] | Liu Q, Yin X Y, Li C. 2016. Quantitative prediction of reservoir characteristics based on inverse modeling theory[J]. Chinese Journal of Geophysics, 59(9): 3491–3502. DOI:10.6038/cjg20160931 |
| [] | Liu X X, Yin X Y, Zhang F. 2013. S-wave velocity estimation method in carbonate reservoir[J]. Journal of China University of Petroleum, 37(1): 42–49. |
| [] | Luo S L, Yang P J, Hu G M, et al. 2016. S-wave velocity prediction based on the modified P-L model and matrix equation iteration[J]. Chinese Journal of Geophysics, 59(5): 1839–1848. DOI:10.6038/cjg20160527 |
| [] | Mavko G, Mukerji T, Dvorkin J. 2009. The Rock Physics Handbook[M]. 2nd ed. Cambridge: Cambridge University Press. |
| [] | Pride S R, Berryman J G, Harris J M. 2004. Seismic attenuation due to wave-induced flow[J]. Journal of Geophysical Research, 109(B1): B01201. |
| [] | Reuss A. 1929. Berechnung der flieβgrenze von mischkristallen auf grund der plastizit?tsbedingung für einkristalle[J]. Zeitschrift für Angewandte Mathematik und Mechanik, 9(1): 49–58. DOI:10.1002/(ISSN)1521-4001 |
| [] | Ruiz F, Dvorkin J. 2010. Predicting elasticity in nonclastic rocks with a differential effective medium model[J]. Geophysics, 75(1): E41–E53. DOI:10.1190/1.3267854 |
| [] | Xu S Y, Payne M A. 2009. Modeling elastic properties in carbonate rocks[J]. The Leading Edge, 28(1): 66–74. DOI:10.1190/1.3064148 |
| [] | Xu S Y, White R E. 1996. A physical model for shear-wave velocity prediction[J]. Geophysical Prospecting, 44(4): 687–717. DOI:10.1111/gpr.1996.44.issue-4 |
| [] | 刘倩, 印兴耀, 李超. 2016. 基于逆建模理论的储层特征定量预测方法[J]. 地球物理学报, 59(9): 3491–3502. DOI:10.6038/cjg20160931 |
| [] | 刘欣欣, 印兴耀, 张峰. 2013. 一种碳酸盐岩储层横波速度估算方法[J]. 中国石油大学学报 (自然科学版), 37(1): 42–49. |
| [] | 罗水亮, 杨培杰, 胡光明, 等. 2016. 基于变形P-L模型的矩阵方程迭代精细横波预测[J]. 地球物理学报, 59(5): 1839–1848. DOI:10.6038/cjg20160527 |