2010, Vol. 53
地下介质的密度是地球物理探测的重要参数之一，基于重力异常的重力勘探是常用的获取地下介质密度参数的地球物理方法.随着研究的深入，重力探测给出的密度估计已不能满足油气勘探生产的需要，从精细程度更高的地震勘探数据中获取密度信息已成为人们努力的一个方向[1~3].
目前，从地震数据中获取密度参数主要利用AVO反演.由于小偏移距的振幅数据主要受控于地震波速度，对密度参数敏感的是中、远偏移距振幅，用AVO反演获取密度参数需要大偏移距地震数据[4~6].一般来说，超过常规地震孔径的大偏移距数据资料品质差、信噪比低、动校拉伸严重、各向异性明显，因此，只利用随偏移距变化的振幅数据的AVO密度反演难度大[7].基于波动方程的全波形反演直接利用炮集数据，减少了常规叠前地震反演需要更多数据处理而对地震信号产生的畸变，同时全波形反演能充分利用接收到的地震波场运动学（如时差）和动力学（如振幅）信息，在一定程度上缓解了常规叠前地震反演只利用振幅而使反演信息量严重不足的矛盾，可以给出比较可靠的密度参数估计[8~10]，因此，利用描述地震波在地下传播规律的波动方程叠前全波形反演受到了越来越多的关注[8, 11].
直接利用地震数据估计孔隙度是非常困难的.目前常规的利用神经网络技术建立地震属性与地层孔隙度[12~14]关系估计孔隙度的方法只是一种数学映射，很难解释其物理意义.依据阻抗、速度或密度等弹性参数估计孔隙度是一种具有物理基础的可靠方法.弹性参数直接与地层孔隙度、孔隙流体类型、饱和度和骨架矿物成分有关.已有的研究表明，相对于阻抗和速度，密度参数更适合于孔隙度估计[15].
本文比较详细地分析了孔隙度和流体饱和度变化对地层密度值的影响，确定了密度的敏感性储层物性参数和密度-孔隙度估计，讨论了与密度相关的地震物理参数：地震波速度、阻抗和振幅，在此基础上开展了基于波动方程的全波形密度反演和孔隙度估计研究，以及速度误差和噪声对反演结果的影响.理论模型研究验证了反演方法的可行性和实现过程的可靠性.我国西部某气田的实际数据叠前密度反演和孔隙度估计结果与测井基本一致，进一步证明了方法的有效性.
2 储层物性与密度为了讨论方便，这里首先将泥质归于固体骨架，用有效固体骨架密度代替骨架物质密度，这时有效介质密度就是孔隙度、饱和度、孔隙流体密度和固体骨架密度的函数，用下式表示：
|
(1a) |
|
(1b) |
其中，ρ，ρs，ρf分别为有效介质、固体骨架和孔隙流体密度，ϕ为孔隙度，Sw为含水饱和度，ρw、ρo（g）分别为水、油（气）的密度.
对于一个研究区来说，ρs和ρw、ρo（g）是确定的，剩下的是孔隙度和含水饱和度.这两个参数对密度值的影响的大小是合理解释密度参数的依据.
2.1 孔隙度与密度式（1a）表示了孔隙度与密度的关系，对该式求导可以得出密度随孔隙度的变化率：
|
(2) |
由（2）式可知，孔隙度对密度的影响程度取决于骨架固体密度和孔隙流体密度.对于我国典型碎屑岩油气藏来说，储层骨架固体密度一般不超过2.8g/cm3.流体密度在完全水饱和时最大，一般在1.05g/cm3左右，完全饱和气时最小，约为0.000187g/cm3.因此，碎屑岩储集层密度对孔隙度的变化率一般在-1.75~-2.80范围（负号表示密度随孔隙度增大而减小）.
图 1显示了完全饱和水、油和气三种情况下密度随孔隙度的变化曲线.可以看出，完全饱和气时曲线的斜率最大，其次是完全饱和油，当孔隙中完全饱和水时斜率最小，说明气藏密度-孔隙度估计准确度更高.
|
图 1 密度随孔隙度变化的曲线 Fig. 1 Density versus porosity |
饱和度对密度的影响体现在式（1b）中.类似于孔隙度，密度对饱和度的变化率表示为
|
(3) |
它与水和油（或气）的密度差值成正比，同时还受控于孔隙度，孔隙度越大，饱和度对密度的影响越大. 图 2显示了孔隙度分别为10%、20%、30%和40%四种情况下气-水系统密度随含水饱和度变化的曲线.可见，孔隙度越大，曲线变化的梯度越大.孔隙度为10%时，梯度为0.1048，表明饱和度每变化1%，密度变化量只有0.001048g/cm3，相对变化量约为0.037%；孔隙度为40%时梯度为0.4192，即饱和度每变化1%，密度变化0.004192g/cm3，相对变化量约为0.15%.对于碎屑岩油气藏，孔隙度一般不超过40%，所以1%的饱和度变化引起密度变化上限不超过0.15%.
|
图 2 密度随饱和度变化的曲线 Fig. 2 Density versus saturation |
比较密度与孔隙度、密度与饱和度两种情况，孔隙度变化对密度参数的影响比饱和度变化对密度参数的影响大的多，其相对变化量约为饱和度的4倍.
3 密度与地震参数与密度直接相关的地震参数主要有速度和阻抗等，表现形式为
|
(4) |
式中，Vp、Vs、Z分别为纵波速度、横波速度和纵波阻抗，ρ为密度，K、μ分别为体积模量和剪切模量.
通过速度和阻抗等参数，密度对地震波的振幅、相位和走时产生影响.
式（4）体现了速度和阻抗与密度的关系，那么密度对速度和阻抗影响见（5）式：
|
(5) |
式中可见，在纵、横波速度和阻抗三个参数中密度对阻抗的影响最大，其次是纵波速度，横波速度影响最小.但密度引起的纵、横波速度和阻抗的相对变化基本相同，约为
地震波的反射系数取决于界面两侧弹性参数的变化.零入射角情况，密度与反射系数（r）的关系可表示为
|
(6) |
保持上覆地层速度和密度不变（V1=5490m/s，ρ1=2.6g/cm3），界面下方地层密度从2.8g/cm3变化到2.0g/cm3，由（4）式的第一式和（6）式计算的反射系数（图 3所示）可见，反射系数相对变化量的大小与密度值有关.一般情况下，当地层含气后，密度变化不会超过0.2g/cm3，假设含气后地层密度从2.55g/cm3变化到2.40g/cm3，即密度相对变化约6%，由（5）和（6）式可得，引起的反射系数变化约为28%，速度变化3%，走时变化6%.因此，利用有效的方法，我们可以从地震信息中提取介质的密度参数.
|
图 3 零入射角反射（a）和反射系数随密度变化曲线（b） 黑线为反射系数随密度变化曲线，浅灰色为反射系数相对变化量. Fig. 3 Reflection at zero angle of incidence (a) and versus density (b) |
常规地震勘探假设地下是声学介质，地震波在地下传播可用下列方程描述：
|
(7) |
式中，“·”表示对时间求导，K为有效介质体积模量，ρ为密度，P为声波场，s为震源，f为牵引力.方程（7）中描述介质的参数有密度和体积模量，因此，利用声波方程反演可以获得这两个参数.
4.1 方法原理在常规地震勘探中，能接收到的地震波场信息非常有限，一般只有地表数据，这时直接从波动方程（7）中求解介质参数是非常困难的.这里我们基于扰动理论，利用Tarantola提出的波场回传技术[8]，由震源正演波场和剩余波场（野外观测数据与合成数据的差逆时传播获得的波场）时间导数的零延迟互相关及正演波场与剩余波场空间梯度的零延迟互相关，分别计算体积模量和密度的共轭修改量：
|
(8a) |
|
(8b) |
式中
有了体积模量和密度的共轭修改量，利用可控共轭梯度法[16]，通过迭代修改密度和体积模量的初始模型，至观测数据与合成记录达到最佳拟合.
4.2 初始模型和先验模型可控共轭梯度法是一种局部优化方法，需要初始模型，同时要求初始模型至少反映模型参数的低频变化.这里，我们首先用速度分析数据和Dix公式计算层速度，并进行时深转换，得到深度域速度模型；再由Gardner公式计算密度[17]，最后由密度和速度计算体积模量，从而建立密度和体积模量初始模型.对于地球物理反演，利用先验模型对反演结果进行约束也是非常重要的，这里，我们用声波测井和密度测井构建先验模型.
4.3 模型和实际数据反演 4.3.1 模型反演我们设计了一个三层理论模型，上、下层速度和密度相同，中间夹有厚度为40 m的低速、低密度层.该层含有一垂直断距为10m的断层，图 4中黑线显示了断层左侧模型密度值随深度的变化，红线为5次迭代后反演的密度曲线，除了一些高频抖动外，基本与黑色曲线一致，说明反演结果是真实模型的一个比较好的估计.
|
图 4 密度随深度变化曲线 Fig. 4 Densities versus depth |
图 5a显示了该模型合成的一个炮集数据，震源为胀缩源.图中除反射纵波，还可见转换波和断层处的绕射波.我们将合成地震记录作为反演的输入数据.图 5b为用初试模型（图 4蓝线所示）合成的地震数据.与图 5a比较，初始模型合成的地震记录与真实模型合成的记录有比较明显的差异.图 5c显示了由最终反演结果合成的地震数据，与图 5a比较，图 5c比图 5b有了明显的改善.图 5d为图 5c与图 5a的差，差值数据的能量比较弱，其均方根振幅与输入数据之比约为15%，表明基于反演得到的模型参数合成的记录与实际数据非常接近.
|
图 5 合成炮集数据及残差 （a）输入数据；（b）初始模型合成数据；（c）反演结果合成数据；（d）（a）和（c）数据差. Fig. 5 Synthetic shot gather and difference (a) Input data; (b) Synthetic data for initial model; (c) Synthetic data for final model; (d) Difference of data between (a) and (c). |
图 6a为经过预处理（包括静校正、地表一致性振幅补偿和去噪）后的单炮数据.频谱宽度可达5~170Hz，主频在60 Hz左右，数据中还有明显的低频和线性斜干扰.反演前，我们对此数据进行了二次处理，包括叠前去噪、室内组合、重采样和数据切除，二次处理后的数据（图 6b）作为反演的输入.图 6c显示了去噪前（深红色）和去噪后（绿色）的频谱，二次去噪有效地消除了低频干扰.
|
图 6 去噪前（a）去噪后（b）数据及频谱（c） Fig. 6 Data before (a) and after (b) denoise and spectra (c) |
图 7显示了井孔处两次迭代反演的密度数据（红线）和测井数据（黑线），可以看出两者基本一致，进一步验证了全波形密度反演的可靠性.
|
图 7 反演密度和测井密度 Fig. 7 Densities of inversion and log |
地震响应特征主要取决于速度和密度，是两者的共同贡献.利用地震数据反演时，如果给出的速度不准确，则速度误差将由密度来祢补.通常情况下获得的速度都有一定的误差，这里我们讨论这种误差对密度反演的影响.
为了简单，考虑零入射情况.将（6）式用微分形式表示：
|
(9) |
其中t为走时，用于表示深度.
（9）式两边对速度求导，令反射系数的导数为零，可得
|
(10) |
式中Δρ和ΔV分别表示密度和速度的变化量；rρ=dρ/ρ和rV=dV/V分别为密度和速度反射率，且r=（rV+rρ）/2.
由（10）式可以看出速度误差对反演密度的影响程度还与密度反射率和速度反射率有关.表 1给出了不同反射率情况下因速度误差而引起的反演密度误差.
|
表 1 速度误差引起的密度误差 Table 1 Densidty error caused by velocity error |
反射率为正情况下，速度误差引起的密度误差随反射率增加而增大，且大于速度误差，如反射率为0.01时，约是速度误差的1.01倍，反射率为0.3时，约是速度误差的1.24倍.反射率为负情况下，速度误差引起的密度误差随反射率绝对值增大而减小，且小于速度误差，如反射率为-0.15时，引起的密度误差只是速度误差的0.92倍，反射率为-0.3时，约是速度误差的0.87倍.
4.5 噪声对密度反演的影响式（8）为基于波动方程的全波形密度反演的基本公式，体积模量和密度的共轭修改量都含有波场（正演波场和剩余波场）导数（时间或空间的导数）对时间的积分和对炮的求和，相当于数据的叠加，因此对噪声具有比较强的压制效应.
我们对前面三层介质模型的合成数据分别加入不同水平的随机噪声，形成反演的输入数据.输入的初始模型同前（如图 4蓝线所示），不同噪声水平的反演结果差别不大.图 8显示了其中信噪比为2的输入地震数据（图 8a）及经过5次迭代后的反演结果（图 8b）.比较图 8b和图 4，可以看出用加入噪声的地震数据进行反演，其结果除了高频抖动幅度较大外，与无噪情况基本一致，反演结果都很好地给出了三层介质密度的分布特征.
|
图 8 信噪比为2的输入数据（a）及反演结果（b） Fig. 8 Data at ratio of S/N=2(a) and inversion (b) result |
地震反演的密度参数在油气勘探中只是一个半成品，它与地下多种因素有关，因此，合理的解释亦是充分发挥密度反演作用的重要一步.
利用公式（1）可得孔隙度估计式：
|
(11) |
式中，ρs可以由实验室测量，也可利用测井数据获得；ρ为反演的有效介质密度；ρf用测井数据和式（1b）计算.
实际工作中，一般不直接利用（11）式计算孔隙度，而将泥质含量和骨架矿物分开，按下列方法进行泥质校正：
|
(12) |
式中sh为泥质含量，可由伽玛测井获得，ρsh为泥质密度，其他符号的意义同前.
图 9显示了由反演密度（图 7）估计孔隙度和测井孔隙度交会图.交会点虽然比较分散，但基本沿45°线（红线）分布，表明全波形密度反演可以给出比较好的地层孔隙度估计.
|
图 9 孔隙度交会图 Fig. 9 Cross plot of porosity |
上述方法用于我国西部某气田密度反演和孔隙度估计.该气田储集层主要特点是：纵向上厚度大，平面上大面积分布，物性差，孔隙度小，属于典型的低孔隙度、低渗透率储层.有效砂岩单层厚度薄、横向变化大，非均质性强，为大面积分布的砂岩储集层中相对高孔、高渗透带，但最大孔隙度不超过16%.
该区地震数据为高密度采集：道间距10 m，炮间距40m，最大偏移距4000m，记录长度5s，1ms采样，双边接收.我们对采集的数据进行了静校正、地表一致性振幅补偿、叠前去噪和室内组合等预处理，并对数据进行频谱分析，确定用于正演模拟的震源子波主频和有效频宽.图 10显示了处理后用于反演的输入地震数据（10a）及其频谱（10b）.
|
图 10 炮集数据（a）及频谱（b） Fig. 10 Shot gather (a) and spectrum (b) |
图 11a显示了两次迭代反演得到的密度剖面（两口井处覆盖了密度测井曲线），紫色至绿色表示密度逐渐增大，变化范围为2.1~2.8g/cm3.图中可以看出密度分布成层状，但横向上很不均匀，这在一定程度上反映了该区辨状的水下分流河道的沉积特征.图 11b是由反演密度预测的孔隙度剖面.与密度剖面不同，孔隙度剖面上有三个明显的成层状高孔隙度带，分别位于层X6~J1z、层X4~X5和J1d的上方，除J1d上覆层外（没有测井孔隙度解释结果），其他两个高孔带与钻井结果基本一致.三个高孔带上孔隙度横向变化也充分反映了该区有效油气储集层横向分布的不均匀性.
|
图 11 全波形密度反演剖面（a）及孔隙度剖面（b） Fig. 11 Inversion density section (a) and porosity section (b) |
地震勘探是提供地下介质精细密度信息的有效工具.与速度参数相比，密度参数对地震响应特征的影响程度小，常规地震方法很难获得稳定的密度参数估计.基于波动方程的全波形地震反演，自动考虑了地震波的反射、透射、绕射、波型转换等行为，同时充分利用了接收的地震波场信息，能给出比较精确的密度参数估计.另一方面，地震反演的密度参数只是一个半成品，合理的解释是充分发挥地震方法在矿藏勘探中作用的一个关键.在油气勘探中，密度参数与岩性、孔隙度和孔隙流体有关.研究结果表明，对于一个研究地区，油气储集层和围岩的岩性基本一定，这时引起密度值变化的主要因素为孔隙度和饱和度.相对于饱和度，孔隙度变化引起的密度变化梯度大，因此利用密度估计孔隙度比估计饱和度更可靠.理论模型结果和实际数据应用效果均表明，全波形地震反演可以获得比较精确的地下介质密度参数，同时利用全波形反演的密度可以很好地估计地层孔隙度，从而为预测油气藏分布和计算油气储量提供了可靠的参数.
| [1] | Yongyi L. Density anomaly indicators and application. 76th SEG Expanded Abstracts, 2006. 299~303 |
| [2] | Li Yongyi. A study on applicability of density inversion in defining reservoirs. 75th SEG Expanded Abstracts, 2005.1646~1650 |
| [3] | 黄平, 涂涛, 刘兴刚. 川东石炭系岩石物性圈闭探索性研究. 石油地球物理勘探 , 2000, 35(3): 326–332. Huang P, Tu T, Liu X G. Research on lithologic trap in Carbonifer-oils system in east Sichuan basin. OGP (in Chinese) , 2000, 35(3): 326-332. |
| [4] | Behural J, Kabir N, Crider R, et al. Density extraction from P-wave AVO inversion: Tuscaloosa Trend example. 77th SEG Expanded Abstracts, 2007. 1466~1470 |
| [5] | Khare V. Density Inversion using Joint PP/PS data: Sensitivity to the Angle Range. 77th SEG Expanded Abstracts, 2007.965~969 |
| [6] | Tarantola A. A strategy for nonlinear elastic inversion of seismic reflection data. Geophysics , 1986, 51(10): 1893-1903. DOI:10.1190/1.1442046 |
| [7] | Roy B, Anno P, Gurch M. Imaging oil-sand reservoir heterogeneities using wide-angle prestack seismic inversion. The Leading Edge , 2008, 67(9): 1192-1201. |
| [8] | Tarantola A. Inversion of reflection data in the acoustic approximation. Geophysics , 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
| [9] | Zhou C, Schuster G. T, Hassanzadeh S, et al. Elastic wave equation traveltime and waveform inversion of crosswell data. Geophysics , 1997, 62(3): 853-868. |
| [10] | Barnes C, Charara M. full-waveform inversion results when using acoustic approximation instead of elastic medium. 78th SEG Expanded Abstracts, 2008. 1895~1899 |
| [11] | Singh S, Sears T, Roberts M, et al. Full Elastic Waveform Inversion: Future of Quantitative Seismic Imaging. 78th SEG Expanded Abstracts, 2008. 1905~1909 |
| [12] | Pramanik A G, et al. Estimation of effective porosity using geostatistics and multiattribute transforms: A case study. Geophysics , 2004, 69(2): 352-372. DOI:10.1190/1.1707054 |
| [13] | 李忠, 贺振华, 巫芙蓉, 等. 地震孔隙度反演技术在川西砂岩储层中的应用与比较. 天然气工业 , 2006, 26(3): 50–52. Li Z, He Z H, Wu F R, et al. Application of seismic porosity inversion techniques to sandstone reservoirs in western sichuang basin and its effect comparsion. Natual Gas Idustry (in Chinese) , 2006, 26(3): 50-52. |
| [14] | 陈遵德, 祝文亮, 王卫华, 等. 孔隙度预测中的地震特征优化方法及应用. 石油学报 , 1999, 20(1): 50–55. Cheng Z D, Zhu W L, Wang W H, et al. Method of seismic characteristic optimization for porosity prediction and its application. Acta Petolei Sinica (in Chinese) , 1999, 20(1): 50-55. |
| [15] | Young G R, Mrinal K S. A comparision of porosity estimates obtained using post-, partial-, and prestack seismic inversion methods: Marco Polo Field, Gulf of Mexico. 77th SEG Expanded Abstracts, 2007.1790~1794 |
| [16] | Yumei S, Wenzhi Z, Hong C. Nonlinear process control of wave-equation inversion and its application in the detection of gas. Geophysics , 2007, 72(1): R9-R18. DOI:10.1190/1.2399450 |
| [17] | Gardner G H F, Gardner L W, Gregory A R. Formation velocity and density: The diagnostic basis for stratigraphic traps. Geophysics , 1974, 39(6): 770-780. DOI:10.1190/1.1440465 |