基于广义Morse标架的地震瞬时属性提取方法

基于广义Morse标架的地震瞬时属性提取方法
【技术领域】
[0001] 本发明属于地球物理勘探中的信号处理领域，设及地震资料的瞬时属性提取方 法，尤其设及一种基于广义Morse标架的地震瞬时属性提取方法。
【背景技术】
[0002] 地震属性是指由地震数据经过数学变换导出的有关地震波的几何形态、运动学特 征、动力学特征和统计学特征，该些特征能够从视觉上反映储层的形态及其含油气性。20世 纪70年代，地震属性分析开始被引入地震解释中。20世纪90年代W来，由于储层描述和= 维数据体解释的需要，地震属性分析技术迅速发展，地震属性在地层构造解释、储层岩性和 物性特征描述、油气藏预测与动态监视方面得到了广泛应用。目前还没有一个公认的地震 属性分类，按照物理意义，地震属性可W分为时间、振幅、频率、相干、衰减等几大类，按照属 性拾取方法，每类又可W分为瞬时属性、单道时窗属性、多道时窗属性。其中瞬时属性是在 地震波到达的位置上拾取的属性，包括瞬时振幅、瞬时相位、瞬时频率、瞬时带宽等等。近年 来，地震属性的数量突增，然而，瞬时属性仍然是地震数据地质解释的支柱，也是目前商业 处理软件必备的技术模块。
[0003] 在信号处理领域，W瞬时振幅、瞬时相位、瞬时频率为代表的信号瞬时属性的概念 由来已久，在G油or和Vilie的工作之后，大量学者在相关领域做了研究，Boashash对该 些工作做了综述。Taner等人于1979年引入复地震道分析，提出通过复地震道求取瞬时属 性的方法，并给出了该些属性的物理意义W及在地震解释中的应用。瞬时振幅与相邻层的 岩性变化及油气聚集有关。瞬时相位反映界面的不连续性、断层、不整合面和层序边界等。 瞬时频率的变化可W有效刻画地层的厚度和岩性变化，指示油气的分布等。Robedson和 Nogami将瞬时属性用于孔隙砂岩薄层厚度估计，化opra和Marfud利用瞬时属性进行不连 续性、断层和横向不连续性检测。Liu和Mar^d利用瞬时属性检测和刻画曲流河的分布并 确定其厚度。Zeng利用瞬时频率异常来指示薄层。Gao等人利用瞬时频率进行地震资料Q 值估计。
[0004] 常用的提取地震瞬时属性的方法可W分为W下S类：
[0005] (1化1化6的变换方法。对地震道做Hilbed变换，转化为复地震道，其中实部为 原地震道数据，虚部为其Hilbed变换。得到复地震道之后，就可W在每个采样点计算振 幅、相位和频率等属性，即瞬时属性。在Taner的工作之后，HUbed变换法被广泛用于计 算地震瞬时属性，至今大部分商业软件仍采用该方法。很多学者对该方法进行了发展和改 进，Barnes于1996年提出二维复地震道分析的概念。Luo等人于2003年提出广义化化ert 变换并给出其在地球物理方面的应用。Barnes于2007年提出加权瞬时频率的概念。Lu和 化ang将HUbed变换方法进行推广，提出基于自适应滤波器计算瞬时频率的方法。
[0006] (2)基于时频分析的方法。时频分析方法利用信号的时频分布来求取瞬时属性。 Boashash等人提出基于时频分析的自适应瞬时频率估计方法。Stankovic等人利用自适应 窗时频分布计算瞬时频率。高静怀等人提出了在相空间计算地震瞬时属性的方法。 等人提出基于窄带谱分析的地震资料瞬时属性提取方法。Stee曲s和化ijkoningen提出 基于二次时频分布的地震层序分析W及属性提取方法。化ang等人提出基于经验模式分解 (EMD)计算瞬时频率的方法。Han等人将完全总体经验模式分解（CEEMD)方法用于提取地 震资料的瞬时频率。
[0007] 做基于反演的方法。Fomel等人提出利用反演的方法来获得局部属性，相比于瞬 时频率，局部频率物理意义更加明确，且应用效果明显。Liu等人提出基于反演的方法，计算 地震资料的瞬时频率。
[0008] 但是，常用的基于HUbed变换估计瞬时参数的方法对噪声很敏感，且由于滤波 器的截断效应，使得计算出的瞬时属性精度低。

【发明内容】

[0009] 本发明目的在于克服现有技术的不足，提供了一种基于广义Morse标架的地震瞬 时属性提取方法，具有良好的抗噪性能和准确性，得到的瞬时频率剖面能够更加清楚地反 映瞬时频率的变化。
[0010] 为达到上述目的，本发明采用W下技术方案：
[0011] 一种基于广义Morse标架的地震瞬时属性提取方法，包括W下步骤：
[0012] 步骤1 ;采集地震数据体；
[0013] 步骤2;计算各道含噪地震数据对应的广义Morse标架系数；
[0014] 实地震道S(t)的小波变换定义为：
[0015]
(22)
[0016] 式中基本小波(t)选用广义Morse小波；
[0017] 广义Morse小波频域表达式为；
[0018]
(23)
[001引其中U(W)为单位阶跃函数，0和丫为小波的参数，且0 >0, 丫 >0,ap，Y为 归一化常量，曰e, 丫三2(e丫/e)日/丫；
[0020] 离散化后的小波族可W表示为：
[0021]
(24)
[00过其中a。尺度因子a的离散化步长且a。> 1，t。为平移因子t的离散化步长；
[0023] 和公式（1)相对应的离散化后的小波变换表示为：
[0024] Cm,n= <s，IDm,n〉，侦）
[00巧]其中Cm，n为标架系数；
[0026] 步骤3;对于各道的广义Morse标架系数迭代得到有效信号对应的系数；
[0027] 定义K为从沪(护)映射到王2 (吸)的算子，将标架系数C= 映射为Z2 (吸） 上的信号，即
（26)
[002引算子K的伴随算子r为从I，(吸)映射到户(妒)的算子，将Z2 (吸)上的信号投影到 标架系数上：
[0029] K*s= <s，Vm，n〉， (27)
[0030]则s=KK*s, (28)
[0031] 算子K为合成算子，r为分析算子和小波标架{> ' m,。}相对应；
[0032] 将含噪信号表示为y=s+n=Kx+n, (29)
[0033] 其中y表示含噪信号，s为不含噪的有效信号，n为高斯白噪声，K为公式（2。中 的合成算子，X表示变换域的系数，通过求解下列优化问题得到有效信号对应的系数宝；
[0034]
(30)
[003引其中e和待分析信号的噪声水平有关，用Lagrange乘子A将该问题转化为W下 的无约束问题：
[0036]
(31)
[0037] 其中A被称为正则化参数，采用下列公式所示的迭代方法对系数X进行更新
[0038]x(w)=T,[x(k)+K*(y-Kx(k))]，k= 1，2,…，N，（32)
[0039] 其中为阔值函数，定义为
[0040]
(33)；
[00川步骤4;由公式似计算解析信号；
[004引
（34)
[0043] 其中Wt(t，a)是有效信号对应的系数，h(t)是s(t)的化化ed变换，c(t)是s(t) 对应的解析信号；
[0044] 步骤5;利用C(t)计算瞬时振幅、瞬时相位和瞬时频率：
[004引其中Re[c(t)]和Im[c(t)]分别表示c(t)的实部和虚部。
[0049] 作为本发明进一步优选方案，所述步骤3公式（26)每一次迭代过程中，采用指数 阔值下降策略，如公式（28)所示：
[0050]
(38)
[00川其中
[0052]
(39)
[005引公式（29)中Am。,和Ami。分别为正则化参数的最小值和最大值，Pm。济Pmi。为最 大和最小百分比。
[0054] 作为本发明进一步优选方案，所述步骤3中为了及时中止迭代过程，定义一个动 态停止准则
[0055]
(40)< br>[0056] 其中tolerance为迭代前给定的容许值，在该准则下，当继续迭代不会取得更好 结果时，迭代过程就自动中止。
[0057] 作为本发明进一步优选方案，所述步骤3中为了提高收敛速度，采用快速迭代萎 缩阔值算法见公式（31);
[0058]
(41)
[005引其中x?= 0,tW= 1，tW满足W下递推公式：
[0060]
(42)。
[0061] 本发明基于广义Morse标架的地震瞬时属性提取方法，采用广义Morse标架，将有 效信号能量分布空间的确定转化成一个优化问题，利用迭代萎缩阔值算法和快速迭代萎缩 阔值算法求解。在确定了有效信号能量分布空间的基础上，利用小波变换与HUbed变换 的关系，提出了含噪信号瞬时属性分析的方法。本发明的方法具有良好的抗噪性能和准确 性，得到的瞬时频率剖面能够更加清楚地反映瞬时频率的变化，指示异常区域。
[0062] 进一步，为了提高算法的收敛速度，本发明采用指数阔值策略和动态停止准则，使 得通过较少的迭代次数，就可W求得有效信号对应的系数。
【附图说明】
[0063] 图1 50化的化cker子波的时间尺度谱图；
[0064] (a)不含噪ricker子波的小波变换谱图
[0065] 化）含噪ricker子波的小波变换谱图
[0066] (C)迭代后的系数构成的谱图
[0067] 图2测试信号和不同方法计算的瞬时频率；
[0068] 图3含噪测试信号和不同方法计算的瞬时频率；
[0069] 图4不同方法计算出来的含噪测试信号的瞬时频率；
[0070] (a)ricke;r子波瞬时频率
[0071] (b)合成地震记录瞬时频率
[0072] (C)实际数据瞬时频率
[0073] 图5为S个测试信号的瞬时频率信噪比；
[0074] (a)ricker子波的FSNR
[00巧]化）合成地震记录的FSNR[00M] (C)实际地震道的FSNR
[0077]图6为实际资料算例；
[007引 （a)含噪地震剖面
[0079] 化）采用化bed得到的法瞬时频率剖面
[0080] (C)采用本发明得到的含噪地震剖面
[0081] 图7为某S维地震资料的瞬时频率沿层切片
[0082] (a)利用商业软件计算的瞬时频率沿层切片
[0083] 化）基于本发明方法计算的瞬时频率沿层切片
[0084] 图8某S维地震资料的瞬时带宽沿层切片
[0085] (a)利用商业软件计算的瞬时带宽沿层切片
[0086] 化）基于本发明方法计算的瞬时带宽沿层切片
[0087] 表1离散化后的小波族构成的标架的上下界和舒适度
[0088] 表2通过迭代求取有效信号对应系数的流程（伪代码）
[0089] 表3不同方法的瞬时频率信噪比FSNR。
【具体实施方式】
[0090] W下通过具体实施例和附图对本发明方案做具体说明：
[0091] 小波变换法计算瞬时属性
[0092] 实地震道S(t)的小波变换定义为
[0093]
(43)
[0094]当基本小波Mt)为解析小波时，高静怀等证明了W下定理：
[009引如果Mt)为解析小波，其实部in(t)为偶函数且其傅里叶变换We(w)满足如 下条件；0 < S￡中* (汾)/份d? <GO，则对任意的信号_^的e巧吸），有下列等式成立
[0096]
(44)
[0097] 其中h(t)是s(t)的HUbed变换，c(t)是s(t)对应的解析信号。该定理建立 了小波变换和化化6的变换之间的关系。因此，本发明利用c(t)来计算瞬时振幅、瞬时相 位和瞬时频率：
[0098]
(45)
[0101] 其中Re[c(t)]和Im[c(t)]分别表示c(t)的实部和虚部。
[0102] 小波函数的选取
[0103] 广义Morse小波是一类双参数解析小波族，最早由Daubechies等人研究时频局域 化算子时导出，是时频联合局域化问题的解。其频域表达式为
[0104]
(48)
[010引其中U(W)为单位阶跃函数，0和丫为小波的参数，且0 > 0, 丫 > 0。ap，Y为归一化常量：
[0106] a日Y三 2(e丫/0)日/丫.（49)
[0107] 本发明选用广义Morse小波计算信号的瞬时属性，主要原因有：
[010引 （1)广义Morse小波是严格解析的，而常用的Morlet小波只是在调制频率足够大 的时候近似解析，在时间局域化要求比较高的时候，Morlet小波会出现负频率泄露。更重 要的是，上面定理中要求基本小波为解析小波，小波偏离解析性会造成该定理不再成立，影 响所求瞬时属性的准确性。而广义Morse小波严格解析，完全满足定理要求；
[0109] (2)广义Morse小波具有更高的自由度，可W通过调整参数，表现出不同的形态和 特点，W匹配不同的地震信号；
[0110] (3)离散后的广义Morse小波更容易构成紧标架，使得小波变换的运算量大大降 低，也保证了后面求解优化问题时算法的收敛性。
[0111] 小波函数的离散化方法
[0112] 小波原子由基本小波进行伸缩和平移得到，即
[0113]
(50)
[0114] 在进行数值计算时，需要对尺度因子a和平移因子t进行离散化，本发明对尺度因 子进行指数化离散，记为{端}，其中a。为离散化步长且a。> 1，平移因子t的离散化步长为 t。。该样，离散化后的小波族可W表示为
[0115]
巧1)
[0116] 在本发明的离散化方法中，t。和待分析信号的采样间隔相同，因此离散化的小波 族满足平移不变性，可W看做由字典?{沁（x) = 口。平移产生。对于具备平移不变性 的字典，有W下定理：
[0117] 如果存在B>A> 0,使得对于we吸，有
[0118]
(52)
[0119] 则该平移不变字典构成一个标架，A和B分别为标架的上界和下界。该定理将字 典构成标架的条件转化为4m(X)的傅里叶变换〇m(?)应满足的条件。
[0120] 对于小波函数，基于本发明的离散化方法，有fl，,,,（W)="，(。:心)，则上述定理的条 件等价为
[0121]
.巧 3)
[0122] 本发明采用舒适度5刻画标架紧的程度；
[0123]
[0124] 舒适度越小，说明标架的下界A和上界B越接近，标架越接近于紧标架。当5 = 0时，上下界相等，标架为紧标架。
[0125]
(55)
[0126] 然后计算在不同的尺度离散化步长a。下，离散化后的Morlet小波和广义Morse 小波（GMW)构成的标架的上下界和舒适度，如表1所示。本发明选取广义Morse小波的参 数为0 = 1，丫 =3,为了确保Morlet小波的近似解析性，参数取0 =6。从表中可W看 出，基于本发明的离散化方式，Morlet小波族和广义Morse小波族均可W构成标架，但是广 义Morse小波的舒适度更小，更容易构成紧标架，特别是在采样步长a。比较大的时候，广义 Morse小波族已经可W被看做紧标架，该可W大大降低运算量。
[0127] 表1离散化后的小波族构成的标架的上下界和舒适度
[012 引
[0129]
[0130] 对于可W构成紧标架的小波族，和公式相对应的离散化后的小波变换和信号重构 可W表示为：
[0133] 其中Cm,。为标架系数，S为离散化的信号构成的向量，（WmJ为小波紧标架， <s，ilv。〉表示二者的内积，A是标架的界（上界与下界相 等）。离散化的小波族构成紧标架， 计算小波变换的系数和从小波系数重构信号变得很方便，同样也保证了后面求解优化问题 时算法的收敛性。
[0134] 小波标架的算子表示
[0135]令
[0136]Vm,n= 巧8)
[0137] 则公式（14)和（15)可W写成
[014引本发明定义K为从映射到王2 (化；)的算子，将标架系数C= (C,,w),,wcH映射为 王2 (化）上的信号，即 [0143]
(62)
[0144] 算子K的伴随算子r为从￡2 (吸)映射到f2 (妒)的算子，将王2 (吸)上的信号投影到 标架系数上：
[0145]K*s= <s,Vm，n〉，化3)
[0146] 则公式（19)可W写成
[0147] S=KK*s,化4)
[0148] 算子K为合成算子，r为分析算子，和小波标架{>' 相对应。
[0149] 有效信号能量分布空间的确定
[0150] 当选用合适的小波函数对含噪信号进行小波变换，将其投影到时间尺度域的时 候，有效信号的能量将会分布在较小的子空间V，噪声的能量会扩散到比较大的子空间 V'，甚至整个时间-尺度域。换句话说，有效信号和少数的系数相对应，而噪声几乎分布在 全部的系数上。如果确定了有效信号对应的子空间V，得到有效信号对应的系数，将其它系 数置零，则噪声就会在变换域被压制，信噪比得到提高。
[0151] 为了确定有效信号的能量分布，希望有效信号投影到尽可能少的小波标架原子 上，也就是说希望信号在变换域得到尽可能稀疏的表示，该样就能压制更多的噪声，提高信 噪比。具体来说，一个含噪信号可W表示为
[0152]y=s+n=Kx+n, (65)
[015引其中y表示含噪信号，s为不含噪的有效信号，n为高斯白噪声，K为公式（2。中 的合成算子，X表示变换域的系数。希望得到一个最优的系数X，使得（1)S--KS，S尽可 能稀疏；(2) ||s-4的值小。稀疏性和1。范数最小化问题相对应，是非凸问题，求解难度很 大，常常在一定条件下转化为li范数最小化问题。因此，可W通过求解下列优化问题得到 X:
[0154]
[0155] 其中e和待分析信号的噪声水平有关。根据Elad等人提出的方法，可W用 Lagrange乘子A将该问题转化为W下的无约束问题：
[0156]
[0157] 其中A被称为正则化参数，该样一来，可W直接采用迭代萎缩阔值算法（1ST)来 求解该问题。Daubechies等人指出，采用下列公式所示的迭代方法对系数X进行更新，当迭 代次数足够大的时候，迭代的结果收敛于公式（24)所示优化问题的解。
[0巧引 x(w)=T,[x(k)+K*(y-Kx(k))]，k= 1，2,…，N，化8)
[0159] 其中为阔值函数，定义为
[0160]
[0161] 需要指出的是，由于迭代萎缩阔值算法（1ST)中采用的是软阔值（如公式（27)所 示），此时阔值和正则化参数相等，用A表示。
[0162] 在求解优化问题过程中，需要采取阔值下降策略，即；随着迭代次数的增加，入的 值不断减小。常用的阔值下降策略有线性下降、指数下降等，Gao等人研究发现，采用指数 阔值下降方法能够大大提高算法收敛的速度。
[0163] 在每一次迭代过程中，采用指数阔值下降策略，如公式（28)所示：
[0164]
[0165] 其中
[0166]
[0167] 公式（29)中Amax和Amin分别为正则化参数的最小值和最大值，Pma济Pmin为最 大和最小百分比。
[016引上面的阔值策略只有在迭代次数达到最大值N的时候才会停止，为了在能够及时 中止迭代过程，本发明定义了一个动态停止准则
[0169]
[0170] 其中tolerance为迭代前给定的容许值。在该准则下，当继续迭代不会取得更好 的结果时，迭代过程就会自动中止，W减少不必要的迭代次数。
[0171] 迭代萎缩阔值算法（1ST)的收敛速度很慢，仅为0(l/k)，即具有一阶收敛速度。 Beck和Teboulle提出了快速迭代萎缩阔值算法（FIST)，见公式（31)，该方法使得算法收敛 速度提高为0(l/k2)，即具有二阶收敛速度。
[0172]
[017引其中x?= 0,tW= 1，tW满足W下递推公式：
[0174]
[01巧]表2给出了求取有效信号对应系数的流程（伪代码），通过迭代，本发明得到了优 化问题的解，即与噪声被压制后的信号相对应的系数X。图l(a)-(b)分别画出了不含噪和 信噪比为5地的50化的化cker子波的变换域谱图，可W看到噪声对有效信号的分布造成 的影响。图1(c)画出了经过迭代W后得到的系数构成的谱图，图中噪声被有效压制，准确 刻画了有效信号的能量分布。
[0176] 表2通过迭代求取有效信号对应系数的流程（伪代码）
[0177]
[0178] 图150化的化cker子波的时间尺度谱图，可W看出在含噪信号的小波变换谱图 中，有效信号能量分布受到噪声的影响，经过迭代W后，噪声得到压制，谱图清晰地反映出 有效信号的能量分布，有了系数S之后，本发明可W得到和有效信号对应的解析信号，进一 步计算瞬时属性。
[0179] 本发明的物质基础是地震数据体，采用的逐道处理办法。具体步骤为：
[0180] 步骤1 ;根据公式（17)计算各道地震数据对应的广义Morse标架系数；
[0181] 步骤2;运用表2的流程进行迭代，得到有效信号对应的系数；
[0182] 步骤3 ;由公式似计算解析信号；
[018引步骤4;根据公式（3)-巧）计算瞬时属性。
[0184] 效果分析
[0185] 合成信号算例
[0186] 本发明用=个测试信号来检验不同方法计算瞬时频率的效果。图2(a)-(c)是分 别是50化的化cker子波，用50化的化cker子波与反射系数序列卷积合成的地震记录 和实际地震道数据。首先，本发明用传统的化化6的方法计算=个信号的瞬时频率，如图 2(d)-(f)所示，然后，采用本发明提出的方法，选用0=1，丫 =3的广义Morse小波标架 来计算瞬时频率，如图2(g)-(i)所示，可见两种方法得到的结果几乎相同。
[0187] 图2为测试信号和不同方法计算的瞬时频率，可见在不含噪情况下，两种方法均 能得到准确的瞬时频率计算结果。
[018引然而，当信号被噪声干扰时，HUbed变换法计算的瞬时频率准确性大大降低。本 发明给S个测试信号加上噪声，其信噪比均为5地（图3 (a)-(c))，然后重新用两种方法计 算瞬时属性，如图3(d)-(f)和图3(g)-(i)所不。
[0189] 图3含噪测试信号和不同方法计算的瞬时频率，可见在HUbed变换法的计算结 果中，有效信号的瞬时频率完全被噪声掩盖，而采用本发明方法的计算结果，即使在强噪声 条件下，仍能够清晰地显示出有效信号的瞬时频率
[0190] 在本发明提出的方法中，通过求解优化问题，得到和有效信号对应的系数，然后计 算解析信号和瞬时频率。可见，即使在噪声比较严重的情况下，该方法仍能够准确得到测试 信号的瞬时频率。而传统的HiAed变换方法得到的结果，有效信号的瞬时频率完全被噪 声影响，无法辨识。
[0191] 下面，将本发明提出的方法和基于HUbed变换方法计算的含噪测试信号的瞬时 频率和不含噪信号的瞬时频率画在一起，并且在本发明提出的方法中选用不同的小波标架 进行计算，结果如图4所示。
[0192]图4不同方法计算出来的含噪测试信号的瞬时频率，其中藍线为含噪测试信号的 瞬时频率，红线为不含噪信号的瞬时频率。分别采用HUbed方法和本发明提出的方法计 算，在本发明提出的方法中，分别选用广义Morse标架（GMW)和Morlet小波标架。可W看 出本发明提出的方法具有更好的抗噪性能和精度，同时，采用广义Morse标架的结果 要优 于Morlet标架
[0193] 图4 (a)-(c)画出了用HUbed变换法计算的瞬时频率，图4(d)-(f)画出了用本 发明提出的方法，选用广义Morse标架计算的瞬时频率，图4(g)-(i)为采用本发明提出的 方法，选用Morlet标架计算的瞬时频率。
[0194] 为了定量衡量不同方法计算出来的瞬时频率的抗噪性能，本发明定义了瞬时频率 信噪比：
[0195]
[0196] 其中f。是不含噪信号的瞬时频率，f。为含噪信号的瞬时频率。
[0197] 各个方法的瞬时频率信噪比FSNR列在表3。
[0198] 表3不同方法的瞬时频率信噪比FSNR
[0199]
[0201] 可W看出，由于噪声的影响，基于HUbed变换计算的瞬时频率准确性大大降低， 甚至无法辨识出有效信号的瞬时频率。而本发明提出的方法计算的瞬时频率具有良好的抗 噪性能和准确度。同时可W看出，采用广义Morse标架的精度要优于常用的Morlet小波离 散化后构成的标架。
[0202] 下面，在本发明提出的方法求取瞬时频率的过程中，分别采用迭代萎缩阔值算法 (1ST)和快速迭代萎缩阔值算法（FIST)求解优化问题，在每一步迭代W后，均计算一次瞬 时频率，求出瞬时频率信噪比FSNR，如图5所示。
[0203] 结果表明，快速迭代萎缩算法（FIST)的收敛速度要优于迭代萎缩阔值算法 (1ST)，可W在较少的迭代次数之后得到优化问题的解。在实际应用中，由于计算速度的要 求，过多的迭代次数不太现实，常常需要设置一个最大迭代次数。在本次实验中，本发明设 置最大迭代次数N= 24,可W看到该方法在10次W内即可得到收敛解。由于动态停止准则 的作用，使得实际用到的迭代次数很少，该可W大大节省计算时间，便于进行大数据体的属 性分析。
[0204] 图5S个测试信号的瞬时频率信噪比，图中红线为采用迭代萎缩阔值（1ST)算法， 藍线采用快速迭代阔值萎缩（FIST)算法，可W看出后者可W在较少的迭代次数之后得到 较高的瞬时频率信噪比
[0205] 实际地震资料算例
[0206] 本发明选取了某油田叠后=维数据体的一个剖面，分别用HUbed变换法和本发 明提出的方法计算瞬时频率，如图6所示。可W看到本发明的方法具有良好的抗噪性能和 准确性，得到的瞬时频率剖面能够更加清楚地反映瞬时频率的变化，指示异常区域。
[0207] 图6实际资料算例。分别采用HUbed变换法和本发明提出的方法计算瞬时频率 剖面，后者具有良好的抗噪性能，可W更加清晰地反映出瞬时频率的变化，指示异常区域 [020引下面，本发明将基于广义Morse标架提取瞬时属性的方法用于某油田致密砂岩储 层的S维地震资料处理，分别利用商业软件和基于广义Morse标架提取瞬时属性的方法计 算该=维地震资料的瞬时频率和瞬时带宽，然后沿目标层位提取岩层切片。利用商业软件 提取的瞬时频率切片如图7 (a)所示，由于常规的商业软件基于传统的化化ed变换方法来 计算瞬时频率，对噪声极为敏感，结果受噪声干扰严重。基于广义Morse标架提取的瞬时频 率沿层切片如图7(b)所示。通过对比可W看到，本发明提出的方法具有很高的抗噪能力， 对目标层砂体的空间分布刻画更加明显。
[0209] 分别采用常规商业软件和本发明方法计算的瞬时带宽沿层切片如图8所示。同样 可W看到前者受噪声影响较大，给地震解释带来很多不便，而后者受噪声影响较小，可W清 晰地刻画地质结构。
[0210] (a)利用商业软件计算的瞬时频率沿层切片
[0211] 化）基于本发明方法计算的瞬时频率沿层切片
[0212] 图7某S维地震资料的瞬时频率沿层切片
[0213] (a)利用商业软件计算的瞬时带宽沿层切片
[0214] 化）基于本发明方法计算的瞬时带宽沿层切片
[0215] 图8某S维地震资料的瞬时带宽沿层切片。
【主权项】
1.基于广义Morse标架的地震瞬时属性提取方法，其特征在于包括以下步骤： 步骤1 :采集地震数据体； 步骤2 :计算各道含噪地震数据对应的广义Morse标架系数； 实地震道s (t)的小波变换定义为：式中基本小波φ (t)选用广义Morse小波； 广义Morse小波频域表达式为： β,y (W) ^ U { , (2) 其中υ(ω)为单位阶跃函数，β和γ为小波的参数，且β > 〇, γ > 〇, α p，γ为归一 化常量，s2(ef/々广 离散化后的小波族可以表示为： Ψ,ηΛΧ) = <'ψ\_<{χ-η^)\ (3) 其中％尺度因子a的离散化步长且a1，t ^为平移因子t的离散化步长； 和公式（1)相对应的离散化后的小波变换表示为： Cm,n= <S, Φ m；n>， (4) 其中Qn为标架系数； 步骤3 :对于各道的广义Morse标架系数迭代得到有效信号对应的系数； 定义K为从(N2)映射到Z2 (R)的算子，将标架系数C = 央射为Z2丨^)上的信 号，算子K的伴随算子矿为从I2(M)映射到f2(N2；)的算子，将上的信号投影到标架 系数上： K*s = <s, Φ rm；n>, (6) 则 s = KK% (7) 算子K为合成算子，矿为分析算子和小波标架{ Φ ' 相对应； 将含噪信号表示为y = s+n = Kx+n, (8) 其中y表示含噪信号，s为不含噪的有效信号，η为高斯白噪声，K为公式（22)中的合 成算子，X表示变换域的系数，通过求解下列优化问题得到有效信号对应的系数i : x = argmm||x||1 s.t. ||y-Kx||2 <ε ^ (9) 其中ε和待分析信号的噪声水平有关，用Lagrange乘子λ将该问题转化为以下的无 约束问题： t = argmjn臺|y-Kx||;+l|4，（ 10) 其中λ被称为正则化参数，采用下列公式所示的迭代方法对系数X进行更新 x(k+1)= I\[X(k)+K*(y-KX(k))]，k = 1，2,…，N， (11) 其中Τλ为阈值函数，定义为步骤4 :由公式（2)计算解析信号；其中Mt, a)是有效信号对应的系数，h⑴是s⑴的Hilbert变换，c⑴是s⑴对 应的解析信号； 步骤5 :利用c (t)计算瞬时振幅、瞬时相位和瞬时频率：其中Re[c(t)]和Im[c(t)]分别表示c(t)的实部和虚部。2. 根据权利要求1所述的基于广义Morse标架的地震瞬时属性提取方法，其特征在于： 所述步骤3公式（26)每一次迭代过程中，采用指数阈值下降策略，如公式（28)所示：公式（29)中λ _和λ min分别为正则化参数的最小值和最大值，pmin为最大和 最小百分比。3. 根据权利要求1所述的基于广义Morse标架的地震瞬时属性提取方法，其特征在于： 所述步骤3中为了及时中止迭代过程，定义一个动态停止准则其中tolerance为迭代前给定的容许值，在该准则下，当继续迭代不会取得更好结果 时，迭代过程就自动中止。4. 根据权利要求1所述的基于广义Morse标架的地震瞬时属性提取方法，其特征在于： 所述步骤3中为了提高收敛速度，采用快速迭代萎缩阈值算法见公式（31);其中χω= 〇, t ω= L t ω满足以下递推公式：
【专利摘要】本发明公开了一种基于广义Morse标架的地震瞬时属性提取方法，采用广义Morse标架，将有效信号能量分布空间的确定转化成一个优化问题，利用迭代萎缩阈值算法和快速迭代萎缩阈值算法求解。在确定了有效信号能量分布空间的基础上，利用小波变换与Hilbert变换的关系，提出了含噪信号瞬时属性分析的方法。本发明的方法具有良好的抗噪性能和准确性，得到的瞬时频率剖面能够更加清楚地反映瞬时频率的变化，指示异常区域。
【IPC分类】G01V1/28
【公开号】CN104880731
【申请号】CN201510141682
【发明人】高静怀, 王平
【申请人】西安交通大学
【公开日】2015年9月2日
【申请日】2015年3月27日