相干体裂缝特征描述是裂缝性油气藏勘探的重要手段之一,而获取用于裂缝预测所需高分辨率数据体以及裂缝参数的求取一直是裂缝油气藏勘探中重点研究的问题。

薄层反射系数谱反演^[1-3]是近些年发展起来的一种能有效分辨薄层、提高分辨率的方法。该方法的主要特点是不需要任何先验信息和数学假设条件,利用反射系数奇偶分解理论和部分频率谱信息反演地层稀疏反射系数,进一步提高分辨低于调谐厚度的地震薄层能力。PURYEAR等^[4]、CHOPRA等^[5]对谱反演算法理论进行了详细推导,模型测试和实际资料应用结果证明该方法能稳健地分辨薄层厚度。YUAN等^[6]详细分析了谱反演的不适定性,并提出了快速的粒子群Levenberg-Marquardt混合反演算法,提高了反演的精度和效率。柴新涛等^[7]利用最小二乘QR分解(Least Squares QR,LSQR)算法进行谱反演,得到较好的模型测试效果。陈祖庆等^[8]利用地震资料的谱信息和基追踪算法进行稀疏脉冲反射系数谱反演,模型试验结果表明该方法能分辨原始数据中无法识别的薄层。

三维相干体技术是20世纪90年代后期兴起的一项十分有效的地震解释技术,该技术主要从相邻地震道相互之间的相干性出发,给出一种定量描述,能突出断层处地震波的变化情况,形成相干值异常区域,从而预测断裂展布。BAHORICH等^[9]提出在振幅上进行互相关相干分析C₁算法,但该算法计算出的断点非相干性点模糊,噪声大、能量弱。MARFURT等^[10]提出基于相似系数相干C₂算法,具有稳定性好、断层刻画精度高的特点,但横向分辨率低的劣势无法改变。GERSZTENKORN等^[11]提出基于本征结构分析的C₃算法,利用Hilbert变换计算相干体,提高其稳定性和抗噪性,但该算法用于高陡复杂地区时,如不考虑倾角则效果不如第二代算法的效果好,如考虑倾角则因计算效率低而无法应用于生产。张军华等^[12]、王玉学等^[13]、刘振峰等^[14]、齐晴等^[15]对C₃算法及其优化算法在断层解释、裂缝预测方面进行了讨论和改进,取得不错的应用效果。

本文在前人研究的基础上,利用基追踪算法是全局优化算法,能使地震信号的稀疏分解更加稳定的优势^[16],进行基于基追踪算法理论的交替方向法(ADM)谱反演算法研究,求取高精度反射系数,再利用该反射系数与含低频的宽频子波重构得到分辨率更高的三维地震数据体,在该数据体上利用改进的C₃相干算法进行高精度断层识别和裂缝预测,避免了解释人员在断层解释和组合的随意性,使断层解释、裂缝预测的精度大大提高。

1 ADM谱反演的高频拓展原理 1.1 稀疏谱反演基本原理

根据不含噪声的褶积模型,地震合成记录在频率域可表示为:

$S\left( f \right)=W\left( f \right)R\left( f \right)$

(1)

式中:S(f),W(f)和R(f)分别表示频率域地震记录、地震子波、反射系数。

稀疏谱反演利用频率域地震记录与子波的丰富信息,采用合适的反演算法消除子波影响,从记录中得到所需的反射系数。利用反射系数奇偶分解原理可得到频率域下谱反演目标函数表达式:

$\begin{align} & O=\int \left[ {{a}_{e}}Re\left[ \frac{S\left( f \right)}{W\left( f \right)} \right]-Re\left[ R\left( f \right) \right]+ \right. \\ & \left. {{a}_{o}}Im\left[ \frac{S\left( f \right)}{W\left( f \right)} \right]-Im\left[ R\left( f \right) \right] \right]df \\ \end{align}$

(2)

式中:a_e和a_o分别代表反射系数R(f)的偶、奇分量比例。

为建立更符合实际地层情况的反射系数模型,必须推导多层反射系数模型下的目标函数,由反射系数奇偶分解原理得到多层反射系数偶、奇分量,表示为:

$\left\{ \begin{align} & {{g}_{e}}\left( t \right)=\int {{r}_{e}}\left( t \right)\left\{ \delta \left[ t-\frac{T\left( t \right)}{2} \right]+\delta \left[ t+\frac{T\left( t \right)}{2} \right] \right\}dt \\ & {{g}_{o}}\left( t \right)=\int {{r}_{o}}\left( t \right)\left\{ \delta \left[ t-\frac{T\left( t \right)}{2} \right]+\delta \left[ t+\frac{T\left( t \right)}{2} \right] \right\}dt \\ \end{align} \right.$

(3)

式中:T表示反射系数时间厚度；r_e和r_o分别代表反射系数r(t)的偶、奇分量。

对(3)式进行傅里叶变换得到其频率域实部与虚部表达式:

$\left\{ \begin{align} & Re\left[ R\left( t,f \right) \right]=\int {{r}_{e}}\left( t \right)cos[\pi fT\left( t \right)]dt \\ & Im\left[ R\left( t,f \right) \right]=\int {{r}_{o}}\left( t \right)sin[\pi fT\left( t \right)]dt \\ \end{align} \right.$

(4)

利用公式(2)和公式(4)得到多层反射模型的目标函数表达式为:

$\begin{align} & O({{r}_{e}},{{r}_{o}},T,t)=\int \left\{ {{a}_{e}}\left\{ Re\left[ \frac{S\left( t,f \right)}{W\left( t,f \right)}{{e}^{-i2\pi f\Delta t}} \right] \right. \right.- \\ & \left. \int {{r}_{e}}\left( t \right)cos[\pi fT\left( t \right)]dt \right\}+{{a}_{o}}\left\{ Im\left[ \frac{S\left( t,f \right)}{W\left( t,f \right)}{{e}^{-i2\pi f\Delta t}} \right] \right.- \\ & \left. \left. \int {{r}_{o}}\left( t \right)sin[\pi fT\left( t \right)]dt \right\} \right\}df \\ \end{align}$

(5)

式中:Δt表示时移量。

为了便于求解,将公式(5)写成矩阵形式:

$O({{r}_{e}},{{r}_{o}})=sum\left| \left[ \begin{matrix} {{a}_{e}}({{b}_{e}}-{{A}_{e}}\times {{R}_{e}}) \\ {{a}_{o}}({{b}_{o}}-{{A}_{o}}\times {{R}_{o}}) \\ \end{matrix} \right] \right|$

(6)

其中,

$\begin{align} & {{A}_{e}}= \\ & \left[ \begin{matrix} cos(\pi {{T}_{1}}{{f}_{1}}) & cos(\pi {{T}_{2}}{{f}_{1}}) & \ldots & cos(\pi {{T}_{N/2}}{{f}_{1}}) \\ cos(\pi {{T}_{1}}{{f}_{2}}) & cos(\pi {{T}_{2}}{{f}_{2}}) & \ldots & cos(\pi {{T}_{N/2}}{{f}_{2}}) \\ \vdots & \vdots & {} & \vdots \\ cos(\pi {{T}_{1}}{{f}_{M}}) & cos(\pi {{T}_{2}}{{f}_{M}}) & \ldots & cos(\pi {{T}_{N/2}}{{f}_{M}}) \\ \end{matrix} \right] \\ & {{A}_{o}}= \\ & \left[ \begin{matrix} sin(\pi {{T}_{1}}{{f}_{1}})s & sin(\pi {{T}_{2}}{{f}_{1}}) & \cdots & sin(\pi {{T}_{N/2}}{{f}_{1}}) \\ sin(\pi {{T}_{1}}{{f}_{2}}) & sin(\pi {{T}_{2}}{{f}_{2}}) & \cdots & sin(\pi {{T}_{N/2}}{{f}_{2}}) \\ \vdots & \vdots & {} & \vdots \\ sin(\pi {{T}_{1}}{{f}_{M}}) & sin(\pi {{T}_{2}}{{f}_{M}}) & \cdots & sin(\pi {{T}_{N/2}}{{f}_{M}}) \\ \end{matrix} \right] \\ & {{b}_{e}}={{[\begin{matrix} {{r}_{e}}({{t}_{1}}) & {{r}_{e}}({{t}_{2}}) & \cdots & {{r}_{e}}({{t}_{N}}) \\ \end{matrix}]}^{T}} \\ & {{b}_{o}}={{[\begin{matrix} {{r}_{o}}({{t}_{1}}) & {{r}_{o}}({{t}_{2}}) & \cdots & {{r}_{o}}({{t}_{N}}) \\ \end{matrix}]}^{T}} \\ & {{R}_{e}}=\left[ \begin{matrix} Re\frac{S(t,{{f}_{1}})}{W(t,{{f}_{1}})}{{e}^{-i2\pi {{f}_{1}}\Delta t}} \\ Re\frac{S(t,{{f}_{1}})}{W(t,{{f}_{1}})}{{e}^{-i2\pi {{f}_{2}}\Delta t}} \\ \vdots \\ Re\frac{S(t,{{f}_{1}})}{W(t,{{f}_{1}})}{{e}^{-i2\pi {{f}_{M}}\Delta t}} \\ \end{matrix} \right] \\ & {{R}_{o}}=\left[ \begin{matrix} Im\frac{S(t,{{f}_{1}})}{W(t,{{f}_{1}})}{{e}^{-i2\pi {{f}_{1}}\Delta t}} \\ Im\frac{S(t,{{f}_{1}})}{W(t,{{f}_{1}})}{{e}^{-i2\pi {{f}_{2}}\Delta t}} \\ \vdots \\ Im\frac{S(t,{{f}_{1}})}{W(t,{{f}_{1}})}{{e}^{-i2\pi {{f}_{M}}\Delta t}} \\ \end{matrix} \right] \\ \end{align}$

(7)

1.2 ADM谱反演算法

求解目标函数((6)式)的算法有很多,如共轭梯度法、模拟退火法、匹配追踪法、基追踪法等,这些算法及混合算法各有优劣,其中,基追踪去噪算法由于其数值计算过程中表现出的稳健性、高效性和精确性越来越受到地球物理学者的青睐,本文采用以基追踪去噪理论为基础的ADM谱反演算法来求解目标函数最优解。

由于该目标函数((6)式)有多个全局最优解,为了降低求解问题时的多解性和不确定性,需要结合实际地震反射系数为稀疏的这一假设条件。因此,加上一个稀疏化约束,将目标函数改为范数求解形式为:

$\begin{align} & O({{r}_{e}},{{r}_{o}})=\left\| \left[ \begin{matrix} {{a}_{e}}({{b}_{e}}-{{A}_{e}}\times {{R}_{e}}) \\ {{a}_{o}}({{b}_{o}}-{{A}_{o}}\times {{R}_{o}}) \\ \end{matrix} \right] \right\|_{2}^{2}+ \\ & \lambda \|{{r}_{e}}+{{r}_{o}}{{\|}_{1}} \\ \end{align}$

(8)

式中:λ为稀疏因子;‖·‖₁表示L₁范数;‖·‖₂²表示Euclidean范数的平方。

根据基追踪去噪算法理论,可将目标函数(8)的求解问题转化为:

$\underset{x\in {{G}^{n}}}{\mathop{min}}\,\frac{1}{2}\upsilon \|Ax-b\|_{2}^{2}+\|x{{\|}_{1}}$

(9)

式中:υ为极小权重值;G为复数集合；n为向量x的维数。

如果引入变量r∈G^m,那么公式(9)等价于公式(10):

$\underset{x\in {{G}^{n}},r\in {{G}^{m}}}{\mathop{min}}\,\left\{ \frac{1}{2\upsilon }\|r\|_{2}^{2}+\|x{{\|}_{1}}:Ax+r=b \right\}$

(10)

式中:m为向量r的维数。

公式(10)对应的增广拉格朗日子问题可表示为:

$\begin{align} & \underset{x\in {{G}^{n}},r\in {{G}^{m}}}{\mathop{min}}\,\left\{ \frac{1}{2\upsilon }\|r\|_{2}^{2}+\|x{{\|}_{1}}-{{y}^{T}}\left( Ax+r-b \right)+ \right. \\ & \left. \frac{\beta }{2}\|Ax+r-b\|_{2}^{2} \right\} \\ \end{align}$

(11)

式中:y^T表示乘子的共轭转置运算;β为罚参数。

利用公式(11)进行基追踪ADM谱反演算法,实现步骤如下。

1) 令k=0对r^k,x^k,y^k赋初始值,并给定υ,β常数值,然后进行算法运算,如果满足终止准则,则完成运算,否则,进行步骤2);

2) 令x=x^k,y=y^k,求解r的子问题得到r^k+1:

${{r}^{k+1}}=\frac{\upsilon \beta }{1+\upsilon \beta }\left[ \frac{{{y}^{k}}}{\beta }-(A{{x}^{k}}-b) \right]$

(12)

3) 令r=r^k+1,y=y^k则关于x的极小化问题公式(11)等价于公式(13):

$\underset{x\in {{G}^{n}}}{\mathop{min}}\,\|x{{\|}_{1}}+\frac{\beta }{2}A\left\| x+{{r}^{k+1}}-b-\frac{{{y}^{k}}}{\beta } \right\|_{2}^{2}$

(13)

那么可以通过公式(14)近似求解来完成对公式(13)的精确求解:

$\left\{ \begin{align} & \underset{x\in {{G}^{n}}}{\mathop{min}}\,\|x{{\|}_{1}}+\beta {{\left( {{g}^{k}} \right)}^{T}}\left( x-{{x}^{k}} \right)+\frac{1}{2\tau }\|x-{{x}^{k}}\|_{2}^{2} \\ & {{g}^{k}}\underline{\underline{\Delta }}{{A}^{T}}\left( A{{x}^{k}}+{{r}^{k+1}}-b-\frac{{{y}^{k}}}{\beta } \right) \\ \end{align} \right.$

(14)

式中:τ为大于零的邻近参数,g^k为x=x^k时二次项的梯度,则公式(14)通过公式(15)求解得到x^k+1:

$\begin{align} & {{x}^{k+1}}=Shrink\left( {{x}^{k}}-\tau {{g}^{k}},\frac{\tau }{\beta } \right) \\ & \underline{\underline{\Delta }}max\left( |{{x}^{k}}-\tau {{g}^{k}}|-\frac{\tau }{\beta },0 \right)\circ sign\left( {{x}^{k}}-\tau {{g}^{k}} \right) \\ \end{align}$

(15)

式中:“°”表示逐元素相乘;Shrink(·)表示一维收缩算子;sign表示符号函数。

4) 令x=x^k+1,r=r^k+1,求解y的子问题得到y^k+1:

${{y}^{k+1}}={{y}^{k}}-\gamma \beta (A{{x}^{k+1}}+{{r}^{k+1}}-b)$

(16)

式中:γ为大于零的常数。

5) 令k=k+1,重复步骤1)至步骤4)进行迭代运算。

1.3 保持低频信息的宽频子波提取

为得到保持低频信息的高分辨率资料,需对求取的反射系数与保持低频信息的宽频子波重构。

零相位子波是具有较高分辨率的子波,但采用不同主频的零相位子波,主频的变化会使其低频相应发生变化,重构后无法保持提频前的低频信息,所以必须获取能保持低频信息的宽频子波振幅谱,从而得到零相位宽频子波。

我们利用谱模拟解析法^[17]从提频前地震数据的振幅谱进行高频拓宽,把该宽频振幅谱作为零相位子波振幅谱,采用傅里叶反变换得到真正具有提频前低频信息的宽频子波。

假设地震子波的振幅谱可以表示为下面光滑的解析函数:

$w\left( f \right)={{f}^{2}}exp\sum\limits_{n=0}^{L}{{{a}_{n}}{{f}^{n}}}$

(17)

式中:f表示频率;L为多项式阶数;a_n为依赖于实际地震记录的待定常数。

对于确定的L,将计算出的a_n(n=0,1,2,…,L)代入(17)式,可以得到一条光滑的拟合曲线|w(f)|,此曲线即是我们所认为的子波振幅谱。

利用公式(17)求出提频前地震数据振幅谱作为子波振幅谱(图 1a),然后通过交互扩展该振幅谱高频部分,低频振幅保持不变,得到宽频子波振幅谱(图 1b),以该振幅谱与零相位谱重构经过傅里叶反变换后得到保持低频信息的宽频子波(图 1c)。

2 改进型特征值法相干体裂缝预测

断裂、裂缝对碳酸盐岩储层有明显的改善作用,便于解释人员利用预测裂缝手段进行裂缝型储层预测。

相干数据体技术主要利用相邻地震道的相似性原理,描述地层及岩性的横向不均匀性。图 2中,地震信号上半部连续性好的同相轴对应于连续性较好的水平地层、倾斜地层等地质体,下半部连续性差的同相轴对应于连续性较差的断层、裂缝等地质体。采用相干数据体算法,对地震数据体的不连续性进行系统分析,使解释人员能有效识别构造和断层的分布情况,避免解释的随意性。

为了尽可能压制噪声干扰,提高相干体计算数据的分辨率,GERSZTENKORN等^[11]提出了第三代特征值相干体计算法。

首先利用MARFURT等^[10]在第二代多道相似性相干体计算中定义的椭圆和矩形分析窗口(图 3)中,以坐标为(x,y)的分析点为中心,时窗中心t=nΔt,则时间窗口内含2k+1个元素的分析道与窗口内的J道数据排列成矩阵:

$U=\left[ \begin{matrix} {{u}_{1,n-K}} & {{u}_{1,n-K+1}} & \cdots & {{u}_{1,n+K}} \\ {{u}_{2,n-K}} & {{u}_{2,n-K+1}} & \cdots & {{u}_{2,n+K}} \\ \vdots & \vdots & {} & \vdots \\ {{u}_{J,n-K}} & {{u}_{J,n-K+1}} & \cdots & {{u}_{J,n+K}} \\ \end{matrix} \right]$

(18)

式中:u_j,m表示第j道在时间t=mΔt－pΔx_j－qy_j沿着视倾角(p,q)插值后的振幅值。

那么矩阵U的第m列是分析体内每个地震道在时间t=mΔt－pΔx_j－qy_j(j=1,2,…,J)的振幅值,则沿着一对视倾角(p,q),时窗中心在t=mΔt的分析窗口内的2k+1个采样点的数据协方差矩阵为:

$\begin{align} & C\left( p,q \right)=U\cdot {{U}^{T}}= \\ & \left[ \begin{matrix} \sum\limits_{m=n-K}^{n+k}{{{u}_{1,m}}{{u}_{1,m}}} & \sum\limits_{m=n-K}^{n+k}{{{u}_{1,m}}{{u}_{2,m}}} & \cdots & \sum\limits_{m=n-K}^{n+k}{{{u}_{1,m}}{{u}_{J,m}}} \\ \sum\limits_{m=n-K}^{n+k}{{{u}_{1,m}}{{u}_{2,m}}} & \sum\limits_{m=n-K}^{n+k}{{{u}_{2,m}}{{u}_{2,m}}} & \cdots & \sum\limits_{m=n-K}^{n+k}{{{u}_{2,m}}{{u}_{J,m}}} \\ \vdots & \vdots & {} & \vdots \\ \sum\limits_{m=n-K}^{n+k}{{{u}_{1,m}}{{u}_{J,m}}} & \sum\limits_{m=n-K}^{n+k}{{{u}_{2,m}}{{u}_{J,m}}} & \cdots & \sum\limits_{m=n-K}^{n+k}{{{u}_{J,m}}{{u}_{J,m}}} \\ \end{matrix} \right] \\ \end{align}$

(19)

设λ_j为分析窗口内地震数据协方差矩阵公式(19)的特征值,且有λ₁≥λ₂≥…≥λ_J,则在以分析点为中心的分析窗口内,沿视倾角(p,q)的J道C₃特征值相干算法公式可表示为:

${{C}_{3}}\left( p,q \right)=\frac{{{\lambda }_{1}}}{\sum\limits_{j=1}^{J}{{{\lambda }_{j}}}}$

(20)

则在以分析点为中心的分析窗口内,沿视倾角(p,q=0)的J道C₃特征值相干算法公式表示为:

${{C}_{3}}={{C}_{3}}\left( p=0,q=0 \right)$

(21)

则沿着所有倾角的C₃特征值相干计算算法公式可表示为:

${{C}_{3}}=\underset{p,q}{\mathop{max}}\,{{C}_{3}}\left( p,q \right)$

(22)

第三代特征值相干算法C₃在子空间中进行计算,比前两代相干算法能更好地压制噪声,但计算量大、效率慢，影响了实际应用时的经济效益。

为此,我们对该算法进行改进,利用基于相似性算法的高效率性估算各点的视倾角和方位角,在大于10倍的分析窗口内计算视倾角和方位角的平均值(p,q);再在以分析点为中心的分析窗口内,沿着区域平滑后的视倾角(p,q)计算地震数据协方差矩阵的特征值及相干值:

${{C}_{3}}={{C}_{3}}\left( p,q \right)$

(23)

3 数值实验及实际资料应用 3.1 断层模型谱反演试验

利用一个断层模型(图 4a)进行ADM谱反演实验,图 4b是对图 4a断层模型正演得到30Hz主频的合成记录,图 4c为本文方法反演的反射系数,图 4d为主频为50Hz重构记录。分析图 4可知,ADM谱反演算法能求出极性和位置准确性较高的反射系数,而且通过重构能得到断层分辨更清晰的合成记录。

3.2 Marmousi模型谱反演试验

利用SEG提供的Marmousi模型进行谱反演反射系数实验,以验证本文方法对纵、横向均存在变化的地震数据反演的适应性。

图 5a为Marmousi模型反射系数,图 5b是35Hz主频的合成记录,图 5c为本文方法反演的反射系数。分析图 5可知,ADM谱反演算法能精确求出含断裂带、深部盐丘两侧高速体等复杂数据体的反射系数,证明了该方法的正确性和对复杂模型较好的适应性。

3.3 实际资料应用

将基于ADM谱反演重构后的高分辨率实际地震数据用于裂缝预测中,以验证该方法在储层预测中的有效性。对图 6a叠前时间域偏移剖面目的层进行ADM谱反演实验,得到反射系数剖面(图 6b),该反射系数与测井系数(蓝色曲线)匹配较好,各层位置和极性求取较佳。用提取的宽频子波进行高分辨率重构,得到高分辨率剖面(图 6c)，其分辨率明显比图 6a的高,且浅、中、深层的同相轴(如白框处)连续性更好、断点更清晰、断层更明了,更有利于后续解释人员进行裂缝预测。分析高分辨率处理前、后的振幅谱(图 7)可知,提频后频带得到明显拓宽,而且由于使用保持低频的宽频子波重构,提频后振幅谱的低频信息得到完整保留,为后续反演等解释提供了丰富的频带信息。

为了分析高分辨率处理前、后的效果,分别提取了20ms的C₃相干切片(图 8),分析图 8可知,利用ADM谱反演重构的高分辨率剖面比处理前剖面能获取分辨率更高的相干切片效果,且切片上裂缝走势更清晰、裂缝刻画细节更精细。

利用过井成像测井资料(图 9)进一步验证高分辨率处理前、后裂缝发育所在断层的真实性。从图 9中可以看出,在2278～2283m段有1个裂缝发育段,在2279m处裂缝发育,通过合成记录精细标定裂缝发育带相当于须三底向下20ms,对最小负曲率属性数据体进行须三底向下20ms切片(图 10a),与过井相干剖面(图 10b)对比后可知,过井点处发育黑色异常,也即存在裂缝发育带,这也印证了提频后剖面(图 6c)箭头所指位置存在明显断层,证明本文方法得到的高分辨率剖面为后续裂缝预测提供了精度更高的依据。

4 结论

本文研究了ADM谱反演高分辨率技术在裂缝预测中应用情况。模型试验和实际资料应用结果表明,利用ADM谱反演算法求取精度高、极性位置正确、能量强弱分明的地层反射系数；与谱模拟解析法提取的宽频子波褶积重构,得到分辨率高、同相轴连续性好、断层清晰、低频信息丰富的剖面；在此剖面上能得到精度更高、细节刻画更清晰的C₃相干切片，便于断层和裂缝的有效识别,避免了解释人员在断层解释和组合时的随意性,大大提高了断层解释、裂缝预测的精度。