2015, Vol. 58
2. 成都理工大学油气藏地质与开发工程国家重点实验室, 成都 610059
2. State Key Laboratory of Oil and Gas Reservoir Geology and Exploration, Chengdu University of Technology, Chengdu 610059, China
We use Gassmann equation to build up the relationship between prestack seismic data and petrophysical parameters. We treat the joint posterior probability of elastic and petrophysical parameters as the objective function under the Bayesian architecture, by expanding the the objective function with Taylor formula, the joint equations composed of physical and petrophysical parameters can thus be obtained, and the derivatives of the elastic parameters with respect to petrophysical parameters are obtained by differentiating. The conjugate gradient method is then used to find the optimal solutions of P-wave velocity, shear velocity, density, porosity, water saturation and clay content.
We apply this inversion method to prestack seismic data set in Sulige gas field, which is one of the largest gas fields in central Ordos Basin in China. The main gas-producing interval is the H8 sandstone reservoir at the bottom of Lower Shihezi Formation, Lower Permian, Upper Paleozoic. Drilling results indicate that this region is a composite sandstone gas reservoir and its heterogeneous reservoirs are characterized by low permeability, low pressure and low abundance, which bring great challenge to the exploitation of this field. From the comparison of the inversion result at the well location using our method and two-step method, we can see that there are slight but noticeable differences between these two methods, especially for the elastic parameters. Our joint inversion is significantly closer to the "true" elastic and reservoir model than the two-step inversion. Later we use 2-D field prestack seismic data to assess the performance of the proposed method, and selected a 160 ms time interval of interest for the reservoir description. Three wells that intercept the section are previously characterized in two-way travel time, and used to calibrate the rock physics model and generate low-frequency starting model. The major seismic events in the inversion result show more continuity along the reference structural horizon direction, and the estimated property have a good match with the corresponding well-log curves, and commensurate with the frequency content of the seismic data, as expected.
We have developed a method for inverting prestack seismic data under well-log constraints derived from Gassmann fluid substitution relations and prior petrophysical models. In our model, we use Gassmann fluid substitution relations calibrated to the well-log data for more robust modeling, which is consistent with common petrophysical knowledge. In this specific setting, the extended application to invert prestack seismic data can estimate water saturation, porosity, clay volume, P-wave velocity, shear wave velocity and density jointly. Compared with the stochastic optimization method, our method utilize gradient information of objective function, the optimal solution can be obtained more quickly in favor of large-scale production mission. Compared with Bosch's poststack deterministic optimization, prestack seismic amplitude have more information varying with the angle of incidence. And compared with Wyllie equation, Gassmann equation is better for establishing the relationship between the reservoir parameters and prestack seismic observation data, which helps a lot to understand the impact of reservoir parameters on the elastic properties.
地震勘探观测数据是地下地质体和岩性流体等信息的综合反映，通过地震观测数据反演获得储层岩石的弹性与物性参数是储层预测和油气藏描述的重要技术手段.自1984年Ostr and er提出从地震观测数据估算岩石弹性参数的AVO技术之后，地震叠前反演方法获得了突飞猛进的发展，已被广泛应用于油气勘探和开发的多个阶段.随着勘探开发难度的不断增加，单纯的弹性参数反演已不能完全满足储层预测和油气藏描述的需要，利用叠前地震资料开展储层弹性与物性参数反演已成为学术界和工业界共同关注的前沿科学问题.
储层弹性与物性参数地震反演有两类方法，即分步反演法和同步反演法(Bosch et al., 2010).分步反演法首先利用地震观测数据反演得到弹性参数，然后根据弹性参数与物性参数之间确定性或统计性的岩石物理关系，利用分类等方法估算孔隙度、饱和度和泥质含量等储层物性参数(Avseth et al., 2001；Mukerji et al., 2001；Saltzer et al., 2005； Bachrach，2006；Sengupta and Bachrach， 2007；Buland et al., 2008；Grana and Rossa， 2010).由于地下储层的复杂性和多变性，不同岩性的弹性参数常常存在一定范围的重叠，同时，弹性参数反演产生的误差可能会进一步导致储层物性参数的估算错误，这些缺点限制了分步反演法在储层预测和油气藏描述中的广泛应用.同步反演法将储层物性参数的地质先验信息、岩石物理模型和地震正演模型相结合，通过迭代应用正演模型直接从地震观测数据反演得到储层弹性与物性参数，受到人们越来越多的关注.同步反演法通常采用随机优化方法作为主要的技术手段，将模型参数和观测数据同时看作随机变量，随机地产生大量可选择的模型，并对这些模型进行计算，将其结果与实际观测数据进行比较，根据预先给定的信息来确定是否接受该模型.随机优化方法的优点是不完全依赖于初始猜测，理论上在反演过程中不会陷入局部极值.利用随机优化方法对油气区进行储层弹性与物性参数的综合推断，既可以产生大量对地震观测数据和岩石物理之间相关联的认识，又可以计算储层物性参数的边缘概率，定量评价反演结果的不确定性(Eidsvik et al., 2004；Gunning and Glinsky， 2007；Bosch et al., 2007；Loures et al., 2006；Larsen et al., 2006；Kjønsberg et al., 2010；Rimstad and Omre， 2010；Ulvmoen and Omre， 2010；Ulvmoen et al., 2010；Grana，2012).但是随机优化方法的缺点也是十分明显的，就是计算工作量大、效率低，难以完全满足大规模地球物理反演的要求.
基于确定性优化方法的同步反演法常采用数理方程或代数的方法求解反演问题，通过反复迭代，求解反演问题的最优解(Spikes et al., 2007；González et al., 2007).确定性优化方法的优点是收敛速度快，适用于大规模的生产任务.Bosch等(2004，2009)在贝叶斯理论的框架下采用数学形式简单的Wyllie方程建立储层物性参数与波阻抗的联系，通过牛顿迭代法对信噪比相对较高的叠后地震资料进行了储层物性参数反演.尽管该方法取得了较好的效果，但Bosch等(2009)指出了该方法存在的两个缺陷：一是地震叠后资料不包含地震反射振幅随入射角变化的信息，难以获得横波速度和密度等弹性参数，对储层物性参数的估计精度也难以保证；二是Wyllie时间平均方程不适用于含泥质的砂岩，在计算固结较差的砂岩孔隙时需要做压实校正，而且还 不能用于含气砂岩的孔隙度计算(Elliott and Wiley， 1975; Saleh and Castagna， 2004).
本文利用Gassmann方程建立储层物性参数与叠前地震观测数据的联系，在贝叶斯反演框架下以储层物性参数和弹性参数联合的后验概率为目标函数，通过将目标函数的梯度用泰勒公式近似展开得到储层弹性与物性参数联合的方程组，其中储层弹性参数对物性参数的梯度用差分形式来表示，最后通过共轭梯度算法同步反演得到了孔隙度、含水饱和度、泥质含量、纵波速度、横波速度和密度等储层弹性与物性参数.与其他岩石物理模型相比，Gassmann 方程能够更好地建立流体饱和岩石有效弹性参数与孔隙大小、孔隙流体之间的联系，有助于认识储层物性参数对储层弹性性质的影响；与随机优化方法相比，利用确定性优化方法求解后验概率分布，在保证方法有效性和准确性的基础上，提升了反演速度，适用于大规模的生产任务.理论试算和苏里格某工区实际资料的应用验证了该方法的有效性和可靠性. 2 地震叠前正演模型
Zoeppritz方程描述了平面波在水平分界面上发生的反射和透射，是AVO技术的理论基础.为了克服由Zoeppritz方程导出的反射系数在形式上复杂且不易进行数值计算的困难，许多学者对Zoeppritz方程进行了简化.其中，Aki和Richard(1980)在假设相邻地层介质弹性参数变化较小的情况下对Zoeppritz方程进行近似，用速度和密度的相对变化量近似地表示反射系数.在此基础上，Gidlow等(1993)用波阻抗的变化量近似地表示了反射系数，该反射系数方程表示为
在Gidlow的基础之上，Hampson等(2005)假设纵波阻抗的对数LP与横波阻抗的对数LS以及密度的对数Ld之间分别存在线性关系，该线性关系可以表示为
岩石的弹性性质是地震波在弹性介质中传播的基础，不同的岩石具有不同的弹性特征.当同一岩石中含有不同流体时，其弹性性质会发生相应的变化，这就是流体替换技术的物理基础.Gassmann方程是流体替换技术的理论支撑，常见的情况是利用岩石骨架速度预测流体饱和岩石的地震速度.本文首 先结合接触Hertz-Mindlin理论和Hashin-Strikman 上边界计算砂岩骨架的模量Kdry和μdry，然后利用Gassmann方程估算流体饱和岩石的模量Ksat和μsat，建立储层物性参数[φ，Sw，clay]T与弹性参数[LP，ΔLS，ΔLd]T之间的联系，进而计算[LP，ΔLS，ΔLd]T和纵波速度、横波速度以及密度.
胶结物沉积在颗粒接触之外的砂岩骨架的体积模量和剪切模量可以用未胶结砂岩模型来计算.假设开始未胶结砂岩的骨架的孔隙度为φ0=0.36和平均配位数C=9的等同球体颗粒紧密任意排列.在该孔隙度下，接触Hertz-Mindlin理论给出静水压力P条件下等同球体颗粒紧密任意排列的骨架有效体积模量KHM和剪切模量μ<sub>HM：
在计算骨架有效体积模量KHM和剪切模量μHM的基础上，要求砂岩骨架不同孔隙度φ的有效模量Kdry和μdry ，可以采用改进的Hashin-Strikman上边界：
低频Gassmann-Biot理论(Mavko et al., 2009)通过下面的公式预测所产生的流体饱和岩石的模量的增加：
在储层弹性与物性参数的地震叠前反演问题中，已知观测数据为叠前地震观测数据dobs，未知模型参数为储层弹性参数 x和物性参数 y，其中，x 是纵波波阻抗的对数LP、横波波阻抗的对数的相对变化量ΔLS、密度的对数相对变化量ΔLd构成的向量，y 是孔隙度φ、含水饱和度Sw、泥质含量clay构成的向量.本文把观测数据与模型参数同时看作为随机变量，反演目标函数是以观测数据为条件的关于模型参数的条件概率，反演结果是条件概率最大时对应的模型参数，可以表示为
根据贝叶斯原理，式(15)可以表示为地震似然函数与模型先验函数的乘积：
在已知储层物性参数 y 的情况下，结合岩石物理模型，储层弹性参数x与物性参数 y 的联合分布P(y，x)可以表示为
在求解式(20)时，为了方便描述，通常假设地质先验信息、岩石物理误差和地震噪声同时服从高斯分布.因此，式(20)可以进一步表示为
又因为式(21)等价于：
令式(22)中的∇S(y n+1，x n+1)=0，得到：
式中，GT是叠前正演模型对弹性参数x的梯度，GT=
是储层弹性参数对物性参数y 的梯度，
.储层弹性参数对孔隙度的梯度可以表示为F′φ，如果岩石物理方程f(y)是线性的，则F′ φ为常数.由于岩石物理方程f(y)是非线性的，其数学表达式较为复杂，一般难以直接对其求导，F′ φ可以采用差分形式来计算：
在文中，Δφ=ΔSw=Δclay=0.01.
为了将式(24)中的∇2S(y n，x n)化简为上三角 矩阵，在式(24)的左右两边分别左乘储层物性参数 y和弹性参数x联合的协方差矩阵C y，x=
的逆矩阵
，得到：
通过共轭梯度法迭代求解式(31)可以得到 Δ x，将Δx代入式(32)可以计算得到Δ y . 4 算法应用
为了检验本文提出的储层弹性与物性参数地震叠前同步反演方法的性能，我们分别用合成地震数据和实际地震资料对其进行测试. 4.1 合成地震数据测试
在合成地震数据测试中，首先利用测井曲线解释得到储层物性参数曲线(图 1d—f)；然后利用岩石物理方程将物性参数转换得到弹性参数，并在弹性参数中将加入均值为0、方差为ε1的高斯噪声(图 1a—c)；最后利用弹性参数合成角道集数据(图 2a)，并在角道集中加入均值为0、方差为ε2的高斯噪声(图 2b)，其中，地震子波为30 Hz零相位的雷 克子波，时间采样间隔为2 ms，最小炮检距为249 m，最大炮检距为3929 m.
 |
图 1 利用测井曲线得到的储层参数 (a)纵波速度;(b)横波速度;(c)密度;(d)孔隙度;(e)含水饱和度;(f)泥质含量. Fig. 1 Reserver parameters derived from well-log (a)P wave velocity;(b)S wave velocity;(c)Density;(d)Porosity;(e)Water saturation;(f)Clay content. |
 |
图 2(a)利用测井数据合成的角道集;(b)加入高斯噪声的角道集Fig. 2(a)The angle gather generated by well-log data;(b)The angle gather affected by a Gaussian noise |
图 3是利用分步法反演得到的结果：首先利用含有高斯噪声的角道集数据反演得到弹性参数，然后根据弹性参数与物性参数之间的岩石物理关系估算孔隙度、含水饱和度和泥质含量.图 4是利用本文方法反演得到的结果，其中，黑色实线表示测井曲线，红色实线表示估计值.从图中可以看出，在高斯噪声环境下，利用分步法反演得到的弹性参数与真实数据之间存在一定的误差，而且弹性参数反演产生的误差进一步导致储层物性参数的估算错误.本文提出的储层弹性和物性参数同步反演方法在目标函数中引入了储层物性参数的先验地质信息和岩石物理约束条件，使其更加符合真实地质情况，反演结果和测井曲线几乎完全吻合.图 5是地震残差随迭代次数的变化关系，可以看出，本文方法稳态性能较好，且具较快的收敛速度.
 |
图 3 利用分步反演法得到的结果(黑色曲线代表测井曲线，红色曲线代表估计值) (a)纵波速度;(b)横波速度;(c)密度;(d)孔隙度;(e)含水饱和度;(f)泥质含量. Fig. 3 The inversion result by two-step inversion(well-log in blue line，estimated values in red line) (a)P wave velocity;(b)S wave velocity;(c)Density;(d)Porosity;(e)Water saturation;(f)Clay content. |
 |
图 4 利用同步反演法得到的结果(黑色曲线代表测井曲线，红色曲线代表估计值) (a)纵波速度;(b)横波速度;(c)密度;(d)孔隙度;(e)含水饱和度;(f)泥质含量. Fig. 4 The inversion result by proposed method(well-log in blue line，estimated values in red line) (a)P wave velocity;(b)S wave velocity;(c)Density;(d)Porosity;(e)Water saturation;(f)Clay content. |
 |
图 5 地震残差随迭代次数的变化关系Fig. 5 Seismic residual evolution with iterations |
研究区位于苏里格气田，目的层是上古生界二叠系底部盒8—山1砂层组，属于沼泽背景下的辫状河沉积，天然气成藏条件复杂，构造对天然气聚集不起主要控制作用.储层经过强烈的成岩作用改造，有效孔隙以次生孔隙为主，具有低孔、低渗和非均质性严重的特点.主要表现为储层段和围岩段的纵横波速度、密度和波阻抗等信息存在一定范围的重叠，利用传统的储层预测方法难以有效地进行储层预测和油藏描述.由于储层物性参数可直接应用于储层流体检测和岩性识别，在解决此类地质问题时具有巨大优势，因此利用测井和地震资料，开展了基于Gassmann方程的储层弹性和物性参数地震叠前同步反演可行性研究.
图 6是该工区地震叠后记录的一条连井剖面，共有601道，道间间距为25 m，时间间隔为2 ms.该剖面经过3口井：S1井、S2井和S3井，各井所在 CMP分别为1560、1691和2007.我们对这条剖面进行了反演，图 7是S2井井旁道反演结果与测井曲线的对比图，图中黑色曲线为测井曲线，红色曲线为反演结果，可以看出反演得到的纵波速度、横波速度、密度、孔隙度、含水饱和度和泥质含量与测井曲线吻合较好.
 |
图 6 叠后地震记录连井剖面Fig. 6 The seismic section taking from Sulige |
 |
图 7 S2井反演结果与测井曲线的对比图(黑色表示测井曲线，红色表示反演结果) (a)纵波速度;(b)横波速度;(c)密度;(d)孔隙度;(e)含水饱和度;(f)泥质含量. Fig. 7 The inversion results of well S2 compared with well log(well log in blue and estimated data in red) (a)P wave velocity;(b)S wave velocity;(c)Density;(d)Porosity;(e)Water saturation;(f)Clay content. |
图 8和图 9分别是图 5所示的连井剖面经过反演得到的储层弹性参数与物性参数剖面.从图中可以看出，含气砂岩段纵波速度和密度略有降低，纵波 速度变化范围为4200～4400 m·s-1，横波速度变化范围为2200~2400 m·s-1，密度变化范围为2.45~2.55 g·cm-3.由于储层段和围岩段的纵横波速度、密度和波阻抗等信息存在一定范围的重叠，利用本文的反演得到的纵横波速度和密度只能识别部分储层，而反演生成的孔隙度、含水饱和度和泥质含量意义却更加明确，对第1950 ms、1970 ms和2000 ms等3套砂体的反映较好.需要注意的是，反演得到的储层物性参数虽然对气层有一定的表征能力，但是也存在一些假象，利用测井曲线建立先验信息模型来约束反演过程，在一定程度上可以使反演结果更符合实际情况.
 |
图 8 储层弹性参数反演剖面 (a)纵波速度;(b)横波速度;(c)密度. Fig. 8 The inversion section of reservoir elastic parameters (a)P-wave velocity;(b)Shear velocity;(c)Density. |
 |
图 9 储层物性参数反演剖面 (a)孔隙度;(b)含水饱和度;(c)泥质含量. Fig. 9 The inversion section of reservoir petrophysical parameters (a)Porosity;(b)Water saturation;(c)Clay content. |
本文提出了一种储层弹性与物性参数地震叠前 同步反演的确定性优化方法，该方法利用Gassmann 方程建立储层物性参数与弹性参数的关系，结合地 震叠前正演，进而建立储层物性参数与叠前地震观测数据的联系，最后通过确定性优化方法直接从地震叠前数据同步反演得到储层弹性与物性参数.与储层弹性与物性参数同步反演随机优化方法相比，本文的方法利用目标函数的梯度信息，能够更快地获得最优解，有利于大规模生产的需要；与Bosch提出的储层弹性与物性参数同步反演确定性优化方法相比，本文的方法利用地震叠前资料，能够通过地震振幅随入射角的变化，获得更多有关弹性和物 性参数的信息；同时，本文的方法可以利用Gassmann 方程建立储层物性参数与叠前地震观测数据的联系，优于Wyllie方程，能够更好地建立流体饱和岩石有效弹性参数与孔隙度、流体饱和度、泥质含量之间的联系，有助于认识储层物性参数对弹性性质的影响.
致 谢 感谢评审专家对本论文提供的建设性的建议.| [1] | Aki K, Richards P G. 1980. Quantitative Seismology: Theory and Methods. San Francisco: Company. |
| [2] | Avseth P, Mukerji T, Jørstad A, et al. 2001. Seismic reservoir mapping from 3-D AVO in a North Sea turbidite system. Geophysics, 66(4): 1157-1176. |
| [3] | Bachrach R. 2006. Joint estimation of porosity and saturation using stochastic rock-physics modeling. Geophysics, 71(5): O53-O63. |
| [4] | Bosch M. 2004. The optimization approach to lithological tomography: Combining seismic data and petrophysics for porosity prediction. Geophysics, 69(5): 1272-1282. |
| [5] | Bosch M, Cara L, Rodrigues J, et al. 2007. A Monte Carlo approach to the joint estimation of reservoir and elastic parameters from seismic amplitudes. Geophysics, 72(6): O29-O39. |
| [6] | Bosch M, Carvajal C, Rodrigues J, et al. 2009. Petrophysical seismic inversion conditioned to well-log data: Methods and application to a gas reservoir. Geophysics, 74(2): O1-O15. |
| [7] | Bosch M, Mukerji T, Gonzalez E F. 2010. Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: A review. Geophysics, 75(5): 75A165-75A176. |
| [8] | Buland A, Kolbjørnsen O, Hauge R, et al. 2008. Bayesian lithology and fluid prediction from seismic prestack data. Geophysics, 73(3): C13-C21. |
| [9] | Eidsvik J, Avseth P, Omre H, et al. 2004. Stochastic reservoir characterization using prestack seismic data. Geophysics, 69(4): 978-993. |
| [10] | Elliott S E, Wiley B F. 1975. Compressional velocities of partially saturated, unconsolidated sands. Geophysics, 40(6): 949-954. |
| [11] | Gidlow P M, Smith G C, Vail P. 1993. Hydrocarbon detection using fluid factor traces: A case study.// 3rd SAGA Biennial Conference and Exhibition. Expanded Abstracts. |
| [12] | González E F, Mukerji T, Mavko G. 2007. Seismic inversion combining rock physics and multiple-point geostatistics. Geophysics, 73(1): R11-R21. |
| [13] | Grana D, Rossa E D. 2010. Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion. Geophysics, 75(3): O21-O37. |
| [14] | Grana D, Mukerji T, Dvorkin J, et al. 2012. Stochastic inversion of facies from seismic data based on sequential simulations and probability perturbation method. Geophysics, 77(4): M53-M72. |
| [15] | Gunning J, Glinsky M E. 2007. Detection of reservoir quality using Bayesian seismic inversion. Geophysics, 72(3): R37-R49. |
| [16] | Hampson D P, Russell B H, Bankhead B. 2005. Simultaneous inversion of pre-stack seismic data.//75th Annual International Metting, Society of Exploration Geophysicists, Expanded Abstracts. |
| [17] | Kjønsberg H, Hauge R, Kolbjørnsen O, et al. 2010. Bayesian Monte Carlo method for seismic predrill prospect assessment. Geophysics, 75(2): O9-O19. |
| [18] | Larsen A L, Ulvmoen M, Omre H, et al. 2006. Bayesian lithology/fluid prediction and simulation on the basis of a Markov-chain prior model. Geophysics, 71(5): R69-R78. |
| [19] | Loures L G L, Moraes F S. 2006. Porosity inference and classification of siliciclastic rocks from multiple data sets. Geophysics, 71(5): O65-O76. |
| [20] | Mavko G, Mukerji T, Dvorkin J. 2009. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media. Cambridge: Cambridge University Press. |
| [21] | Mukerji T, Jørstad A, Avseth P, et al. 2001. Mapping lithofacies and pore-fluid probabilities in a North Sea reservoir: Seismic inversions and statistical rock physics. Geophysics, 66(4): 988-1001. |
| [22] | Rimstad K, Omre H. 2010. Impact of rock-physics depth trends and Markov random fields on hierarchical Bayesian lithology/fluid prediction. Geophysics, 75(4): R93-R108. |
| [23] | Saleh A A, Castagna J P. 2004. Revisiting the Wyllie time average equation in the case of near-spherical pores. Geophysics, 69(1): 45-55. |
| [24] | Saltzer R, Finn C, Burtz O. 2005. Predicting V shale and porosity using cascaded seismic and rock physics inversion. The Leading Edge, 24(7): 732-736. |
| [25] | Sengupta M, Bachrach R. 2007. Uncertainty in seismic-based pay volume estimation: Analysis using rock physics and Bayesian statistics. The Leading Edge, 26(2): 184-189. |
| [26] | Spikes K, Mukerji T, Dvorkin J, et al. 2007. Probabilistic seismic inversion based on rock-physics models. Geophysics, 72(5): R87-R97. |
| [27] | Ulvmoen M, Omre H. 2010. Improved resolution in Bayesian lithology/fluid inversion from prestack seismic data and well observations: Part 1—Methodology. Geophysics, 75(2): R21-R35. |
| [28] | Ulvmoen M, Omre H, Buland A. 2010. Improved resolution in Bayesian lithology/fluid inversion from prestack seismic data and well observations: Part 2—Real case study. Geophysics, 75(2): B73-B82. |