CN106054247A - 基于转换波地震数据求取高精度反射系数的方法 - Google Patents

Abstract

本发明提供了一种基于转换波地震数据求取高精度反射系数的方法，所述方法包括：采集转换波地震数据，进行预处理和叠前偏移后得到转换波偏移剖面；在偏移剖面上建立稀疏约束和波阻抗约束的多薄层目标函数；对多薄层目标函数利用频域映射与L1范数联合优化算法求取高精度反射系数。本发明的方法能够基于转换波地震数据构建更加合理的多薄层目标函数并获得高精度反射系数，提高了目标函数求解精度和抗噪性，而且所获得的高精度反射系数能够为后续高分辨率剖面重构奠定坚实基础。

Description

基于转换波地震数据求取高精度反射系数的方法

技术领域

本发明属于转换波地震资料处理技术领域，具体来讲，涉及一种基于转换波地震数据求取高精度反射系数的方法。

背景技术

目前，相干体裂缝特征描述是探寻裂缝性油气藏的一项重要手段，而获取裂缝预测所需高分辨率数据体、裂缝参数的求取一直是裂缝油气藏勘探中重点研究的问题，尤其转换波低信噪比和分辨率的特点，需要提高其分辨率，便于后续解释的薄层分辨。但是，最小相位和白噪假设为前提的常规高分辨率处理方法，在主频带范围限制下，提高转换波分辨率的能力有限，无法达到后续裂缝预测分辨多薄层最小厚度的要求，所以需要采取新的高分辨率处理思路，即：转换波高精度反射系数求取与高频子波高分辨率重构的新思路，该思路不受任何假设和频带限制，能完成常规手段无法分辨多薄层最小厚度的要求，获取裂缝预测满意的高品质剖面。

通常，反射系数的求取需要构建多薄层目标函数，然而，由于常规方法构建的目标函数存在诸多不合理之处，故而，在求取过程中的反射系数存在极值和位置、大小上的误差，这种误差导致后续高分辨率剖面的错像。

此外，目标函数的求解方法有很多，比如模拟退火法、共轭梯度法、匹配追踪法等这些算法，虽然都能求取目标函数的反射系数，但精度和稳定性较差，且对新构建的目标函数适应性较差。

综上，亟需一种能够基于转换波地震数据求取高精度反射系数的方法。

发明内容

本发明的目的在于解决现有技术存在的上述不足中的至少一项。例如，本发明的目的之一在于改善现有转换波地震数据领域中为求取反射系数所构建的多薄层目标函数的合理性，并求解出高精度反射系数。

为了实现上述目的，本发明提供一种基于转换波地震数据求取高精度反射系数的方法。所述方法包括：采集转换波地震数据，进行预处理和叠前偏移后得到转换波偏移剖面，在转换波偏移剖面上取n道地震记录x(x,t)，x＝1,2,…,n；在偏移剖面上建立稀疏约束和波阻抗约束的多薄层目标函数；对多薄层目标函数利用频域映射与L1范数联合优化算法求取高精度反射系数；

其中，

所述目标函数为下式(1)：

$O(r_e,r_o)=\left|\right|\left[\begin{array}{c}{a}_e(b_e-A_e\×r_e)\\ a_o(b_o-A_o\×r_o)\end{array}\right]|{|}_2^2+\mathrm{\λ}||r_e+r_o|{|}_1+\∂\left|\right|\frac{{(Cr-{\mathrm{\ξ}}_p)}^T(Cr-{\mathrm{\ξ}}_p)}{2}|{|}_1,$

其中，r表示反射系数，a_e代表反射系数的偶分量比例，a_o代表反射系数的奇分量比例，r_e为反射系数的偶分量，r_o为反射系数的奇分量，A为变换矩阵，b为地震数据和子波相关矩阵，λ为稀疏因子，为阻抗因子，C为积分算子矩阵，ξ_p为纵波波阻抗矩阵。

与现有技术相比，本发明的有益效果包括：能够基于转换波地震数据构建更加合理的多薄层目标函数，并能有效提高目标函数求解精度和抗噪性，能够获得高精度反射系数，进而能够为后续高分辨率剖面重构奠定坚实基础。

具体实施方式

在下文中，将结合示例性实施例来详细说明本发明的基于转换波地震数据求取高精度反射系数的方法。

在本发明的一个示例性实施例中，基于转换波地震数据求取高精度反射系数的方法可通过以下步骤来实现：

(1)对采集到的转换波地震数据，经过常规预处理和叠前偏移后得到的转换波偏移剖面，取n道地震记录x(x,t)，x＝1,2,…,n；其中，x表示道数，t表示时间。

(2)在偏移剖面上建立由稀疏约束和波阻抗约束双重约束的多薄层目标函数；

步骤(2)可按下列方式完成：

(a)对于不含噪声的褶积模型，地震合成记录在频率域中可表示为：

S(f)＝W(f)·R(f) (1)

其中，S(f)表示频率域下地震记录(即，上面的x(x，t))，W(f)表示频率域下地震子波，R(f)表示频率域下反射系数；f为频率，时间t经过傅里叶变换后变换到频率f中。

利用频率域下地震记录与地震子波的丰富信息，利用合适的常规反演算法(如，共轭梯度法、匹配追踪法、模拟退火等)消除子波(即，地震子波)影响(如，子波主频、旁瓣干扰等会影响反演求取反射系数精度的内容)，从记录中得到所需的反射系数(这里，为了得到反射系数R(f))。利用反射系数奇偶分解原理，可得到频率域下目标函数表达式：

$O=\∫\[a_e\{\mathrm{Re}\[\frac{S\left(f\right)}{W\left(f\right)}\]-\mathrm{Re}\[R\left(f\right)\]\}+a_o\{\mathrm{Im}\[\frac{S\left(f\right)}{W\left(f\right)}\]-\mathrm{Im}\[R\left(f\right)\]\}\]df---\left(2\right)$

其中，a_e代表反射系数R(f)的偶分量比例，a_o代表反射系数R(f)的奇分量比例；Re表示R(f)的实部，Im表示R(f)的虚部。

(b)为建立更符合实际地层情况的反射系数模型，对步骤(a)中的多薄层反射系数模型下的目标函数进行推导或改进。

由奇偶反射系数分解原理，得多薄层反射系数奇偶分量可表示为：

$\left\{\begin{array}{c}{g}_e\left(t\right)=\{r_e\left(t\right)\[\mathrm{\δ}(t-\frac{T\left(t\right)}{2})+\mathrm{\δ}(t+\frac{T\left(t\right)}{2})\]dt\\ g_o\left(t\right)=\{r_o\left(t\right)\[\mathrm{\δ}(t-\frac{T\left(t\right)}{2})+\mathrm{\δ}(t+\frac{T\left(t\right)}{2})\]dt\end{array}\right.---\left(3\right)$

r_e(t)代表反射系数r(t)的偶分量，r_o(t)代表反射系数r(t)的奇分量，t表示时间，T表示时间厚度；r(t)为时间域反射系数，r(t)经过傅里叶变换后到上面的频率域反射系数R(f)；g(t)当前分析位置反射系数，g_e(t)代表g(t)的偶分量，g_o(t)代表g(t)的奇分量。

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

$\left\{\begin{array}{c}\mathrm{Re}\[R(t,f)\]=\∫r_e\left(x\right)cos\left(\mathrm{\π}fT\right(t))dt\\ \mathrm{Im}\[R(t,f)\]=\∫r_o\left(t\right)\sin \left(\mathrm{\π}fT\right(t))dt\end{array}\right.---\left(4\right)$

其中，π为圆周率。

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

$\begin{array}{c}O(r_e,r_o,T,t)=\∫\[a_e\{\mathrm{Re}\[\frac{S(t,f)}{W(t,f)}{e}^{-i2\mathrm{\π}f\mathrm{\Δ}t}\]-\∫r_e\left(t\right)\cos \left(\mathrm{\π}fT\right(t))dt\}\\ +a_0\{\mathrm{Im}\[\frac{S(t,f)}{W(t,f)}{e}^{-i2\mathrm{\π}f\mathrm{\Δ}t}\]-\∫r_o\left(t\right)\sin \left(\mathrm{\π}fT\right(t))dt\}\]df\end{array}---\left(5\right)$

其中，r_e代表反射系数r(t)的偶分量，r_o代表反射系数r(t)的奇分量，Δt表示时移量；i无物理意义，是指数与正弦函数变化标识符。

为了便于求解，将式(5)可写成矩阵形式：

$O(r_e,r_o)=sum\left(\left[\begin{array}{c}{a}_e(b_e-A_e\×r_e)\\ a_o(b_o-A_o\×r_o)\end{array}\right]\right)---\left(6\right)$

其中，

其中，t为时间，t₁、t₂、……、t_N分别表示第1、2、……、N个时间，f为频率，f₁、f₂、……、f_M分别表示第1、2、……、M个分析频率范围内的频率，T₁、T₂……、T_N/2分别表示第1、2、……、N/2个时间厚度，Δt表示时移量

(c)由于该目标函数(6)(即，式(6))有多个全局最优解，为了降低求解问题时的多解性和不确定性，需要加上稀疏化约束，将目标函数改为范数求解形式为：

$O(r_e,r_o)=\left|\right|\left[\begin{array}{c}{a}_e(b_e-A_e\×r_e)\\ a_o(b_o-A_o\×r_o)\end{array}\right]|{|}_2^2+\mathrm{\λ}||r_e+r_o|{|}_1---\left(7\right)$

其中，λ为稀疏因子，λ一般在0-1范围内取值，||||₁表示L1范数，表示欧几里得(Euclidean)范数的平方。

(d)由于波阻抗可看作是反射系数对时间的积分，因此发明人把它引入到目标函数中作为反射系数r的先验约束。

由于转换波相对于纵波来说其信噪比和分辨率较低,为了提高反射系数精度,本发明使用分辨率更高的纵波波阻抗作为转换波的约束条件。

若I(t)为纵波波阻抗值，则纵波相对波阻抗ε_t为：

${\mathrm{\ϵ}}_t=\frac{1}{2}ln\frac{I\left(t\right)}{I\left(t_0\right)}={\∫}_{t_0}^{t}r\left(\mathrm{\η}\right)d\mathrm{\η},t=t_0,...,t_{M-1}---\left(8\right)$

其中，I(t₀)表示纵波初始波阻抗值，r为反射系数，η表示在时间t₀与时间t中取值符合，t₀表示初始时间，t_M-1表示结束时间。

把上式(8)写成矩阵形式：

ξ_p＝Cr (9)

其中C为积分算子矩阵，其离散形式可表示为：

则由最小平方定义波阻抗约束得到新的目标函数为：

$O(r_e,r_o)=\left|\right|\left[\begin{array}{c}{a}_e(b_e-A_e\×r_e)\\ a_o(b_o-A_o\×r_o)\end{array}\right]|{|}_2^2+\mathrm{\λ}||r_e+r_o|{|}_1+\∂\left|\right|\frac{{(Cr-{\mathrm{\ξ}}_p)}^T(Cr-{\mathrm{\ξ}}_p)}{2}|{|}_1---\left(11\right)$

其中，r表示反射系数，a_e和a_o分别代表反射系数r(t)的偶奇分量比例，r_e是反射系数的偶分量，r_o是反射系数的奇分量；A为变换矩阵，b是地震数据和子波相关矩阵，λ是稀疏因子，是阻抗因子(其大于零即可)，C为积分算子矩阵，ξ_p是纵波波阻抗矩阵)

(3)对于步骤(2)中的式(11)，利用频域映射与L1范数联合优化算法求取高精度反射系数。

本步骤(3)可按下列方式完成：

(a)对式(11)进行求解中，由于波阻抗约束项是由测井信号得到的已知先验约束条件，所以式(11)的求解类似于求矩阵方程Ax＝b在含有噪声情况下的稀疏解。

可将目标函数(11)的求解问题转化为下式：

$\underset{x\∈G^n}{\min }\frac{1}{2\mathrm{\υ}}\left|\right|Ax-b|{|}_2^2+|\left|x\right|{|}_1---\left(12\right)$

其中，υ为极小权重值，G为复数集合；本式(12)中，A表示A_e时b表示b_e，X表示r_e；A表示A_o时b表示b_o，X表示r_o。

如果引入变量q∈G^m，那么式(12)等价于式(13)：

$\underset{x\∈G^n,r\∈G^m}{\min }\left\{\frac{1}{2\mathrm{\υ}}\right|\left|q\right|{|}_2^2+\left|\right|x|{|}_1:Ax+q=b\}---\left(13\right)$

式(13)对应的增广拉格朗日子问题可表示为：

$\underset{x\∈G^n,r\∈G^m}{\min }\left\{\frac{1}{2\mathrm{\υ}}\right|\left|q\right|{|}_2^2+\left|\right|x|{|}_1-y^T(Ax+q-b)+\frac{\mathrm{\β}}{2}||Ax+q-b|{|}_2^2---\left(14\right)$

其中，y^T表示拉格朗日乘子的共轭转置运算，β为罚参数。

利用式(14)进行频域映射与L1范数联合优化算法实现步骤如下：

1)令k＝0对q^k，x^k，y^k赋初始值，并给定υ，β常数值(这里，υ，β为人为给定，可为大于零的常值)，然后进行以下步骤2)至5)的算法运算，如果满足终止准则(例如，终止准则可以为：最优解值或者迭代次数)，则算法完成运算；否则，进行步骤2)；

2)令x＝x^k，y＝y^k，求解r的子问题得到q^k+1：

$q^{k+1}=\frac{\mathrm{\υ}\mathrm{\β}}{1+\mathrm{\υ}\mathrm{\β}}(\frac{y^k}{\mathrm{\β}}-({\mathrm{Ax}}^k-b\left)\right)---\left(15\right)$

3)令q＝q^k+1，y＝y^k则关于x的极小化问题式(14)等价于式(16)：

$\underset{x\∈G^n}{\min }\left|\right|x|{|}_1+\frac{\mathrm{\β}}{2}||Ax+q^{k+1}-b-\frac{y^k}{\mathrm{\β}}|{|}_2^2---\left(16\right)$

那么可以通过式(17)近似求解来完成对式(16)的精确求解：

$\left\{\begin{array}{c}\underset{x\∈G^n}{\min }\left|\right|x|{|}_1+\mathrm{\β}\left({\left(g^k\right)}^T\right(x-x^k)+\frac{1}{2\mathrm{\τ}}||x-x^k|{|}_2^2)\\ g^k\stackrel{\mathrm{\Δ}}{=}{A}^T({\mathrm{Ax}}^k+q^{k+1}-b-\frac{y^k}{\mathrm{\β}})\end{array}\right.---\left(17\right)$

其中，τ为大于零的邻近参数(这里，τ为人工给定，其大于零即可)，g^k为x＝x^k时二次项的梯度，那么式(17)通过下式(18)求解得到x^k+1：

其中，ο表示逐元素相乘，Shrink(,)表示一维收缩算子，sign表示符号函数。

4)令x＝x^k+1，q＝q^k+1，求解y的子问题得到y^k+1：

y^k+1＝y^k-γβ(Ax^k+1`+q^k+1-b) (19)

其中，γ为大于零的常数(这里，γ人工给定，其大于零即可)。

5)令k＝k+1，重复步骤1)～步骤4)进行迭代运算。

本发明的方法能够基于转换波地震数据构建更加合理的多薄层目标函数并获得高精度反射系数，提高了目标函数求解精度和抗噪性，而且所获得的高精度反射系数能够为后续高分辨率剖面重构奠定坚实基础。

尽管上面已经结合示例性实施例描述了本发明，但是本领域普通技术人员应该清楚，在不脱离权利要求的精神和范围的情况下，可以对上述实施例进行各种修改。

Claims (4)

1.一种基于转换波地震数据求取高精度反射系数的方法，所述方法包括：采集转换波地震数据，进行预处理和叠前偏移后得到转换波偏移剖面，在转换波偏移剖面上取n道地震记录x(x,t)，x＝1,2,…,n；

在偏移剖面上建立稀疏约束和波阻抗约束的多薄层目标函数；

对多薄层目标函数利用频域映射与L1范数联合优化算法求取高精度反射系数；

其中，

所述目标函数为下式(1)：

$O\left(r_e,r_o\right)=\left|\right|\left[\begin{array}{c}{a}_e(b_e-A_e\×r_e)\\ a_o(b_o-A_o\×r_o)\end{array}\right]|{|}_2^2+\mathrm{\λ}||r_e+r_o|{|}_1+\∂\left|\right|\frac{{(Cr-{\mathrm{\ξ}}_p)}^T(Cr-{\mathrm{\ξ}}_p)}{2}|{|}_1,$

其中，r表示反射系数，a_e代表反射系数的偶分量比例，a_o代表反射系数的奇分量比例，r_e为反射系数的偶分量，r_o为反射系数的奇分量，A为变换矩阵，b为地震数据和子波相关矩阵，λ为稀疏因子，为阻抗因子，C为积分算子矩阵，ξ_p为纵波波阻抗矩阵。

2.根据权利要求1所述的基于转换波地震数据求取高精度反射系数的方法，其中，

$\left\{\begin{array}{c}{r}_e={\[r_e\left(t_1\right),r_e\left(t_2\right),\mathrm{.....},r_e\left(t_N\right)\]}^T\\ r_o={\[r_o\left(t_1\right),r_o\left(t_2\right),\mathrm{.....},r_o\left(t_N\right)\]}^T\end{array}\right.,$

$\left\{\begin{array}{c}{A}_e=\left[\begin{array}{cccc}\cos \left({\mathrm{\πT}}_1{f}_1\right)& \cos \left({\mathrm{\πT}}_2{f}_1\right)& ...& \cos \left({\mathrm{\πT}}_{N/2}{f}_1\right)\\ \cos \left({\mathrm{\πT}}_1{f}_2\right)& \cos \left({\mathrm{\πT}}_2{f}_2\right)& ...& \cos \left({\mathrm{\πT}}_{N/2}{f}_2\right)\\ ...& ...& ...& ...\\ \cos \left({\mathrm{\πT}}_1{f}_M\right)& \cos \left({\mathrm{\πT}}_2{f}_M\right)& ...& \cos \left({\mathrm{\πT}}_{N/2}{f}_M\right)\end{array}\right]\\ A_o=\left[\begin{array}{cccc}\sin \left({\mathrm{\πT}}_1{f}_1\right)& \sin \left({\mathrm{\πT}}_2{f}_1\right)& ...& \sin \left({\mathrm{\πT}}_{N/2}{f}_1\right)\\ \sin \left({\mathrm{\πT}}_1{f}_2\right)& \sin \left({\mathrm{\πT}}_2{f}_2\right)& ...& \sin \left({\mathrm{\πT}}_{N/2}{f}_2\right)\\ ...& ...& ...& ...\\ \sin \left({\mathrm{\πT}}_1{f}_M\right)& \sin \left({\mathrm{\πT}}_2{f}_M\right)& ...& \sin \left({\mathrm{\πT}}_{N/2}{f}_M\right)\end{array}\right]\end{array}\right.,$

$\left\{\begin{array}{c}{b}_e=\left[\begin{array}{c}\mathrm{Re}\[\frac{S\left(t,f_1\right)}{W\left(t,f_1\right)}{e}^{-i2{\mathrm{\πf}}_1\mathrm{\Δ}t}\]\\ \mathrm{Re}\[\frac{S\left(t,f_2\right)}{W\left(t,f_2\right)}{e}^{-i2{\mathrm{\πf}}_2\mathrm{\Δ}t}\]\\ ...\\ \mathrm{Re}\[\frac{S\left(t,f_M\right)}{W\left(t,f_M\right)}{e}^{-i2{\mathrm{\πf}}_M\mathrm{\Δ}t}\]\end{array}\right]\\ b_o=\left[\begin{array}{c}\mathrm{Im}\[\frac{S\left(t,f_1\right)}{W\left(t,f_1\right)}{e}^{-i2{\mathrm{\πf}}_1\mathrm{\Δ}t}\]\\ \mathrm{Im}\[\frac{S\left(t,f_2\right)}{W\left(t,f_2\right)}{e}^{-i2{\mathrm{\πf}}_2\mathrm{\Δ}t}\]\\ ...\\ \mathrm{Im}\[\frac{S\left(t,f_M\right)}{W\left(t,f_M\right)}{e}^{-i2{\mathrm{\πf}}_M\mathrm{\Δ}t}\]\end{array}\right]\end{array}\right.,$

其中，t为时间，t₁、t₂、……、t_N分别表示第1、2、……、N个时间，f为频率，f₁、f₂、……、f_M分别表示第1、2、……、M个分析频率范围内的频率，T₁、T₂……、T_N/2分别表示第1、2、……、N/2个时间厚度，Δt表示时移量，Re表示频率域下反射系数R(f)的实部，Im表示频率域下反射系数R(f)的虚部，S表示频率域下地震记录x(x,t)，W表示频率域下地震子波；

其中，t₀表示初始时间，t_M-1表示结束时间；

C的离散形式表示为

3.如权利要求1所述的基于转换波地震数据求取高精度反射系数的方法，其中，所述求取高精度反射系数的步骤包括：

A、将多薄层目标函数式(1)的求解问题转化为式(2)：

其中，υ为极小权重值，G为复数集合；

B、引入变量q∈G^m，则式(2)等价于式(3)：

$\underset{x\∈G^n,r\∈G^m}{\min }\left\{\frac{1}{2\mathrm{\υ}}\right|\left|q\right|{|}_2^2+\left|\right|x|{|}_1:Ax+q=b\};$

C、式(3)对应的增广拉格朗日子问题表示为式(4)：

其中，y^T表示乘子的共轭转置运算，β为罚参数；

D、利用式(4)进行频域映射与L1范数联合优化算法求取高精度反射系数。

4.如权利要求3所述的基于转换波地震数据求取高精度反射系数的方法，其中，所述步骤D包括步骤：

1)令k＝0对r^k，x^k，y^k赋初始值，并给定υ，β常数值，然后进行下面步骤2)至5)的算法运算，若满足终止准则，则完成运算，否则，进行步骤2)；

2)令x＝x^k，y＝y^k，求解r的子问题得到是(5)：

$r^{k+1}=\frac{\mathrm{\υ}\mathrm{\β}}{1+\mathrm{\υ}\mathrm{\β}}(\frac{y^k}{\mathrm{\β}}-({\mathrm{Ax}}^k-b\left)\right);$

3)令r＝r^k+1，y＝y^k则关于x的极小化问题式(4)等价于式(6)：

$\underset{x\∈G^n}{min}\left|\right|x|{|}_1+\frac{\mathrm{\β}}{2}||Ax+r^{k+1}-b-\frac{y^k}{\mathrm{\β}}|{|}_2^2,$

那么可以通过式(7)近似求解来完成对式(6)的精确求解，式(7)为：

$\left\{\begin{array}{c}\underset{x\∈G^n}{\min }\left|\right|x|{|}_1+\mathrm{\β}\left({\left(g^k\right)}^T\right(x-x^k)+\frac{1}{2\mathrm{\τ}}||x-x^k|{|}_2^2)\\ g^k\stackrel{\mathrm{\Δ}}{=}{A}^T({\mathrm{Ax}}^k+r^{k+1}-b-\frac{y^k}{\mathrm{\β}})\end{array}\right.,$

其中，τ为大于零的邻近参数，g^k为x＝x^k时二次项的梯度，那么式(7)通过式(8)求解得到x^k+1，式(8)为：

其中，ο表示逐元素相乘，Shrink(,)表示一维收缩算子；

4)令x＝x^k+1，r＝r^k+1，求解y的子问题得到式(9)：

y^k+1＝y^k-γβ(Ax^k+1`+r^k+1-b)，

其中，γ为大于零的常数。

5)令k＝k+1，重复步骤1)～步骤4)进行迭代运算。