0 引言

随着石油工业从增加产量转向价值优化，石油行业对储层预测的精度要求越来越高。由于储层的渗流能力反映储层品质，因此在评价中、深层优质储层时，不能仅考虑孔隙度参数，而要同时兼顾渗透率、流体性质等其他因素。但迄今为止，人们还基本停留在通过研究渗透率的地震响应特征评价储层^[1-2]，而很难基于地震资料反演渗透率或预测储层。因此，近年来人们将目光逐渐由渗透率转向储层流体的可动性，如流体流度属性可很好地表征流体的可动性。有人^[3-5]以渗流理论为基础首次推导了饱和流体弹性介质波动方程，并证明该方程与试井分析中常用的Frenkel-Gassmann-Biot多孔弹性模型、压力扩散模型^[6]相关，最终得到饱和流体储层低频域反射系数渐近表达式，认为低频反射系数正比于流体流度(Fluid Mobility，又称迁移率)、岩石体积密度、地震信号频率的平方根。Goloshubin等^[7-8]应用流体流动性和散射机制，推导出与地震频率相关的流体流度属性，并用该属性预测储层中流体的流动能力和渗透率。Chen等^[9]基于实际地震资料的流体流度属性识别流体，取得了较好的效果。张生强等^[10]基于高分辨率反演谱分解近似计算流体流度属性，进一步提高了流体流度属性剖面的分辨率。其中高分辨率反演谱分解使用稀疏反演算法将地震资料投影到分频子波库，再将所得结果进行重排，从而得到一个稀疏的时频谱, 其实质仍是一种分频方法。

由于在流体流度属性计算过程中，需要计算地震数据的瞬时振幅谱，不同时频分析方法获得的流体流度属性剖面的分辨率存在差异，进而影响储层预测精度。因此为提高储层预测精度，有必要选择时频分析方法以获得高分辨流体流度属性剖面。小波变换(WT)继承和发展了短时傅里叶变换局部化的思想，同时能够提供一个随频率改变的“时—频”窗口，是信号时频分析的较理想工具，但小波变换的结果是时间—尺度谱，而非时频谱。为此，Stockwell等^[11]提出了S变换，S变换不仅反映信号在频率域的变化，而且可以实现无损重构。为了使窗函数更灵活，人们又提出了广义S变换，如陈学华等^[12-14]引入了两个调节参数λ、p，可根据实际地震资料分解尺度的需要灵活调节，以使地震资料的分辨率达到所需的分析尺度。与小波变换相比，S变换、广义S变换等引入信号频率参数控制窗口尺度，其变换结果是更为直观的时间—频率谱，便于分析信号的频率特征。受测不准原理的制约，以上变换的时频谱分辨率仍然有限。Wigner-Ville分布作为一种非线性时频分析方法，摆脱了测不准原理的制约，但对于复杂信号会产生交叉项干扰^[15]。Lu等^[16-17]提出反褶积短时傅里叶变换，不仅提高了时、频分辨率，且减少了交叉项干扰。朱恒等^[18]首次利用反褶积短时傅里叶变换提取单频剖面预测薄储层。张懿疆等^[19]首次利用反褶积广义S变换分析薄储层的低频伴影，取得较好的效果。本文利用反褶积广义S变换提取流体流度属性，通过仿真试验和实例试算验证了方法的有效性和实用性。

1 方法原理 1.1 流体流度属性的定义及提取方法

Korneev等^[3]、Silin等^[4]、Goloshubin等^[8]推导了饱和流体介质在低频域地震响应的反射系数线性渐近表达式。平面纵波在干岩石与饱含流体弹性介质界面上的反射系数为

$ R = {R_0} + {R_1}(1 + {\rm{i}})\sqrt {\frac{{\kappa {\rho _{\rm{b}}}}}{\eta }\omega } $

(1)

式中：i为虚数单位；R₀和R₁是实系数，为反映岩石和流体特性的函数，包括孔隙度、密度和弹性系数信息；κ为储层渗透率；η为流体黏滞系数；ρ_b为储层内流体的密度；ω为低频地震信号的角频率。

将式(1)对ω求导，得

$ \frac{{{\rm{d}}R}}{{{\rm{d}}\omega }} = \frac{{{R_1}(1 + {\rm{i}})}}{2}\sqrt {\frac{{{\rho _{\rm{b}}}}}{\omega }} \sqrt {\frac{\kappa }{\eta }} $

(2)

Goloshubin等^[7]在低频渐近分析理论的基础上提出了成像属性的概念，即

$ A(x, y) = \frac{{{\rm{d}}R\left( {{\omega _{{\rm{low}}}}} \right)}}{{{\rm{d}}\omega }} \sim \frac{{{\rm{d}}a\left( {{\omega _{{\rm{low}}}}} \right)}}{{{\rm{d}}\omega }} $

(3)

式中：R(ω_low)为低频反射系数；a(ω_low)为低频端地震振幅；~代表等价关系；x、y为储层的横向坐标。令$C = \frac{{{R_1}(1 + {\rm{i}})}}{2}\sqrt {\frac{{{\rho _{\rm{b}}}}}{\omega }} $(孔隙流体和岩石骨架力学性质以及频率的复函数)，则

$ A(x, y) \approx C{\left( {\frac{\kappa }{\eta }} \right)^{\frac{1}{2}}} = C{m^{\frac{1}{2}}} $

(4)

式中$m = \frac{\kappa }{\eta }$为流体流度，其与振幅对频率的一阶导数$\frac{{{\rm{d}}a}}{{{\rm{d}}\omega }}$的关系为

$ m = \frac{\kappa }{\eta } \sim \frac{1}{{{C^2}}}{\left( {\frac{{{\rm{d}}a}}{{{\rm{d}}\omega }}} \right)^2} $

(5)

Silin等^[4-5]指出，当含流体孔隙介质其结构属性值在合理范围内时，依赖频率的地震反射系数在低频端达到峰值，而对地震数据进行瞬时谱分解后，单频瞬时振幅或能量能准确地刻画该频率的地震反射振幅或能量，且利用具有高分辨率和良好时频局部化性能的时频分析方法可减少薄互层调谐对地震瞬时振幅计算结果的影响。因此，文中依据文献[9]、[10]的近似计算流体流度属性的方法，对地震资料的低频成分进行时频分析获得$\frac{{{\rm{d}}a}}{{{\rm{d}}\omega }}$，再将$\frac{{{\rm{d}}a}}{{{\rm{d}}\omega }}$代入式(5)，最终获得储层流体流度属性m。式中C可通过岩石物理测试得到。在了解已知井的产能或渗透率参数的情况下，通过抛物线拟合求得参数C。在缺少井资料的探区，C相当于一个比例系数，可在不依赖井的情况下计算储层相对流体流度，定性分析储层流体^[9-10]。

1.2 优势频率及基质体积模量的求取

在求取流体流度属性过程中，需要求取峰值频率，Silin等^[5]在推导纵波入射渗透地层的反射理论公式时，获得了快纵波反射共振峰值频率

$ {f_{\max }} = \frac{\kappa }{{2\pi \eta {H^2}}}\frac{{2M}}{{\gamma \beta + \gamma {K^2}}} $

(6)

式中：H为目的层厚度；$M = K + \frac{{4\mu }}{3}$, μ是干岩石剪切模量；$\gamma \beta = K\left( {{\beta _{\rm{f}}}\phi + \frac{{1 - \phi }}{{{K_{{\rm{fg}}}}}}} \right)$; $\gamma K = 1 - \frac{{(1 - \phi )K}}{{{K_{{\rm{sg}}}}}}$；而

$ {K_{{\rm{sg}}}} = \frac{{{K_{\rm{g}}}}}{{1 - \phi }} $

$ {K_{{\rm{fg}}}} = \frac{{{K_{\rm{g}}}}}{{1 - \frac{K}{{{K_{\rm{g}}}}}}} $

式中：K_sg和K_fg分别为骨架的体应变和流体压力变化造成的基质颗粒压缩，把颗粒体应变与骨架应力、流体压力联系起来；K_g为基质体积模量；ϕ为孔隙度；K为干岩石体积模量。以上各参数中，K_g最难求取。为此，本文采用贺锡雷^[20]的LRM线性拟合法求取K_g。通过对Biot-Gassmann流体替换方程进行合理改造，可获得其近似公式。考虑到难以直接由测井资料获得骨架模量或干岩石模量以及孔隙结构明显影响岩石有效弹性模量等因素，将含有孔隙结构信息的干岩石椭球包体近似公式引入流体替换方程近似式，与常规计算岩石基质模量的方法相比，LRM法简单易行且精度较高。图 1为根据南海A区岩石物理参数得到的压缩系数(体积模量的倒数)—孔隙度交会图。根据LRM方法可知，压缩系数为截距B(在负半轴，取绝对值)与斜率A的比值。表 1是据A区54块岩样计算的压缩系数及拟合系数A、B。由表可见，岩样饱和不同流体求出的压缩系数基本是一致的，在实际计算过程中取平均值即可。

1.3 反褶积广义S变换

陈学华等^[12-13]引入高斯窗函数的调节参数λ＞0、p＞0，通过调节λ、p改变高斯窗随频率的变化趋势，从而灵活地分析、处理具体信号。对于信号x(t)，其广义S变换可表示为

$ {\mathop{\rm GST}\nolimits} (\tau , f) = \int_{ - \infty }^\infty x (t)\frac{{\lambda |f{|^p}}}{{\sqrt {2\pi } }}{{\rm{e}}^{\frac{{ - {\lambda ^2}{f^{2\beta }}{{(\tau - t)}^2}}}{2}}}{{\rm{e}}^{ - {\rm{i}}2\pi ft}}{\rm{d}}t $

(7)

式中：τ为用于确定分析时窗的时间位移；f为频率。

广义S变换谱为原始信号与高斯窗各自的Wigner-Ville分布的二维褶积，即

$ \begin{array}{l} {{\mathop{\rm SPEC}\nolimits} _x}{\left| {{{{\mathop{\rm GST}\nolimits} }_x}(t, f)} \right|^2} = \int {\int_{ - \infty }^\infty {{\rm{WV}}{{\rm{D}}_x}} } (u, v) \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{WV}}{{\rm{D}}_h}(t - u, f - v){\rm{d}}u{\rm{d}}v \end{array} $

(8)

式中：WVD_x为原始信号x(t)的Wigner-Ville分布(简写为W_x )；WVD_h为高斯窗h(u)的Wigner-Ville分布(简写为W_h)。本文采用非线性的迭代复原反褶积算法——Lucy-Richardson反褶积算法^[16-17]计算广义S变换谱，其表达式为

$ W_x^ - (k + 1) = W_x^ - (k)\left( {{W_h}*\frac{{{S_x}}}{{{W_h} \otimes W_x^ - (k)}}} \right) $

(9)

式中：k+1是目前的迭代次数；$W_{x}^{-}(0)=S_{x}$为原始信号的广义S变换；*为相关算子；⊗为褶积算子。

2 仿真实验

首先利用合成信号验证广义S变换、反褶积短时傅里叶变换、反褶积广义S变换的时频分析效果。合成信号x(t)为由3个线性信号叠加而成的时变信号，即

$ x(t) = \left\{ \begin{array}{l} \cos \left[ {2\pi \cdot \left( {\frac{{20}}{{512}}} \right)} \right]t\;\;\, 0 \le t \le \frac{{512}}{3}\\ \cos \left[ {2\pi \cdot \left( {\frac{{60}}{{512}}} \right)} \right]t\;\;\, \frac{{512}}{3} < t \le \frac{{1024}}{3}\\ \cos \left[ {2\pi \cdot \left( {\frac{{100}}{{512}}} \right)} \right]t\;\;\, \frac{{1024}}{3} < t \le 512 \end{array} \right. $

(10)

图 2为合成信号及利用不同分析方法得到的时频谱。由图可见：利用广义S变换得到的时频谱中频率发生突变的位置(红色箭头所示)不清晰(图 2b)；对信号做反褶积短时傅里叶变换得到的时频分布的边界相对清楚(图 2c)，但仍不明显；对信号做反褶积广义S变换得到的时频分布(图 2d)较好地区分了不同频率成分的信号，显示出较高的时频分辨率与时频聚集性。

利用黏滞—弥散波动方程进行数值模拟，该方程是一种考虑流体黏滞性和弥散性的声波方程，可以很好地体现含流体储层对不同频率地震波的衰减特性^[21-23]。图 3为透镜体储层模型及由不同方法得到的流体流度属性剖面。由图可见：基于反褶积短时傅里叶变换得到的储层流体流度剖面中目的层的高流体流度属性区域与透镜体位置不完全对应(图 3c)；基于反褶积广义S变换得到的储层流体流度剖面中高流体流度属性区域与透镜体的位置较为一致，尤其对透镜体储层两侧尖灭点刻画较好，验证了方法的高分辨率特性(图 3d)。

3 实例分析

为了进一步验证文中方法的实效性，从南海A区的地震数据中提取了流体流度属性。该区位于南海北部陆坡前缘，处于长源富砂珠江组下方，是深水勘探的有利地区，已发现多个油气田和具有潜在商业价值的含油气构造。工区内储层纵波速度约为3500~4200m/s，地震资料主频为25~30Hz，砂体平均厚度约为30m。地温梯度、埋深、成岩作用类型及强度共同控制砂岩储集性能的优劣。研究区优质储层埋深变化较大，主要分布范围为1000~3000m，地温梯度G分布范围为3~7℃/100m。在低地温梯度条件下(G＜4℃/100m)，储层原生、次生孔隙发育，渗透率较大，即使埋深达3000m时孔隙度也可达到15%，渗透率可达数十毫达西；在中、高地温梯度(G≥4℃/100m)条件下，储层次生孔隙发育、原生孔隙次之，形成中孔、低渗储层。工区横向地温梯度及纵向埋深差异造成强烈的储层物性非均质性，致使储层物性评价面临巨大挑战。

图 4为A—B—C—D—E—F井连井地震剖面，图 5、图 6分别为利用广义S变换和反褶积广义S变换获得的流体流度属性剖面。由图可见：由于流体在砂岩中的流动性较高，对应于高流度属性值(图 5、图 6)；基于广义S变换的流体流度属性剖面虽然在储层位置呈流体流度属性高值，但高值区域的储层厚度大于井中解释的储层厚度(椭圆区域所示)，即剖面的分辨率不高(图 5)；基于反褶积广义S变换的流体流度属性剖面高值区与储层解释结果吻合很好，具有较高的分辨率(图 6)。

4 结束语

(1) 与广义S变换相比，反褶积广义S变换的时频分辨率与频率聚集性进一步提高，适合处理具有时变特征的地震数据。

(2) Goloshubin流体流度属性与储层渗透率和流体黏滞系数具相关性，在理论上表征了储层流体的活动能力。实例分析展示了由地震数据得到的流体流度属性预测储层的可行性。

(3) 精确获取地震频率信息是获得高精度流体流度属性的关键，不同时频分析方法获得的流体流度属性剖面的分辨率存在一定差异，本文通过引入反褶积广义S变换，提高了流体流度属性剖面的分辨率。

尚需指出，流体流度属性目前在流体性质识别方面上存在一定缺陷，因此需结合其他方法(如叠前反演或AVO分析)评价储层。