技术领域
[0001] 本
发明属于
地震勘探领域，特别是有关于采用匹配地震子波的物理小波，对地震二维或三维资料进行瞬时频率分析的方法。
背景技术
[0002] 地震
信号的瞬时属性是地震资料分析的重要工具。瞬时属性包括：瞬时振幅、瞬时频率、瞬时
相位、瞬时带宽等。瞬时属性可以用于分析
地层岩性变化、地下构造及反演地层中的岩性参数。分析地震信号的瞬时频率有多种途径：包括时间-频率域方法、解析信号方法等。前者应用的时-频分析方法有：短时Fourier变换、
小波变换、S变换及广义S变换、Wigner分布、匹配追踪方法等；后者主要应用复信号方法。
[0003] 由于地震信号是实函数，在解析信号方法中需计算其对应的复信号。计算地震信号对应的复信号方法有多种(B.Boashash，1992，L.Conhen，1994；Tanner，1979)，但在地震
信号处理领域最常用的方法还是Hilbert变换法。这种方法对地震数据进行单道处理，首先计算待分析
地震道信号的Hilbert变换，将其结果作为
虚部，该道地震信号作为
实部，构成复信号，然后利用该复信号计算瞬时频率(Tanner，1979)。但是这种方法存在如下
缺陷：1)不能提供多
分辨率瞬时属性，即无法对地下地质体进行最佳分辨率解释等(高静怀等，
1997)；2)对噪声敏感，难以用于低
信噪比资料(Gao Jinghuai，et al，1999；高静怀，汪文秉，朱光明等，1997)；3)由复信号计算瞬时频率时，采用二阶中心差商，当地震资料频带较宽时误差大。
[0004] 为了克服缺陷1)，小波变换的创始人之一Morlet提出了基于小波变换的地震资料多分辨率解释方法(Morlet，1989)，目的是提高地震资料解释的可靠性。众所周知，信号的小波变换结果既与待分析的信号有关，也与采用的小波函数有关，但是Morlet提出的小波变换难以匹配待分析的地震反射信号，使得该方法分析地震信号的能
力受到限制。为了克服缺陷2)，在Morlet等人工作的
基础上，高静怀等人讨论了小波变换用于地震资料多分辨率分析时小波函数的选择问题，并构造出了适用于这类问题的小波-匹配地震子波的小波，下文称这类小波为匹配地震子波的物理小波，简称物理小波(高静怀等，1996；2001)。为了能在小波域计算地震信号对应的解析信号，高静怀等建立了待分析信号的解析小波变换与其Hilbert变换之间的关系，给出了在小波域计算地震信号对应的解析信号的方法(Gao Jinghuai et al，1999；2001)。为了克服缺陷3)，高静怀等提出了基于极平坦
滤波器的瞬时频率计算方法。然而对于实际资料，如何构造恰当的物理小波、如何针对具体的地质目标选择最佳分辨率以及极平坦滤波器中参数等，没见到相关报道。
发明内容
[0005] 针对上述问题，本发明目的是提供一种基于匹配地震子波的物理小波的地震资料瞬时频率分析方法。该方法弥补了基于Hilbert变换方法对噪声敏感的问题，也弥补了基于小波变换的多分辨率解释方法难以匹配待分析的地震反射信号的问题。
[0006] 为解决上述问题，本发明采取的技术方案是：一种基于匹配地震子波的物理小波的地震瞬时频率分析方法，包括如下步骤：
[0007] 1)首先获取二维或三维的经偏移或
叠加处理后的地震资料；
[0008] 2)根据研究的对象和目的对所获取的二维或三维地震资料进行空间分区，在所划分的区域里获得
测井资料或零偏VSP资料、井旁地震记录、储层的地质构造及其它先验信息；
[0009] 3)通过测井资料或零偏VSP资料，以及井旁地震记录反演地震子波，确定匹配该子波的母物理小波，按如下步骤进行：
[0010] ①反演地震子波
[0011] 利用分区内的测井资料和井旁地震记录，反演出地震子波，或利用垂直地震剖面VSP资料得到地震子波；
[0012] ②确定匹配地震子波的母物理小波
[0013] 基于如下表达式：
[0014] g(t；α)＝Aexp[-τ(t-β)2]exp(iσt)+R(t；α) (1)[0015] 式中g(t，α)为解析小波，为书写简便，下文把g(t；α)简记为g(t)，R(t；α)为修正项，表达式为：
[0016]
[0017] α为一矢量，定义为α＝(A，σ，τ，β)，A为地震子波的幅度，σ为母小波的调制频率，τ为母小波的
能量衰减率，β为母小波的能量延迟时间，
[0018] 取(1)式的实部和第1)步所得到的地震资料构造如下目标函数：
[0019]
[0020] (2)式中，real表示取实部，取(2)式达到极小值时对应的A、σ、τ、β四个参数，代入到(1)式即得到匹配地震子波的物理小波；
[0021] 4)在物理小波域计算地震信号对应的解析信号，采用多尺度解析信号计算和在小波域有效信号能量分布空间计算两种方法：
[0022] ①多尺度解析信号计算
[0023] (1)式定义的g(t)满足：当ω＜0时 因此，g(t)为解析小波，并且[0024] g(t)∈L1(R，dt)∩L2(R，dt)
[0025] 和
[0026]
[0027] 所以任给一个地震信号s(t)∈L2(R，dt)，s(t)相对于g(t)的小波变换定义为：
[0028]
[0029] 这里t∈R，a∈R\{0}，b∈R，t，b都表示时间，函数 表示对其取复共轭；对任一尺度因子a(a＞0)，S(b，a)即为在尺度因子为a(a＞0)时，s(t)对应的解析信号，a(a＞0)决定了分辨率。
[0030] ②在小波域有效信号能量分布空间计算地震信号对应的解析信号[0031] 此方法中采用下式计算地震信号s(t)对应的解析信号：
[0032]
[0033] 这里，S(b，a)由(3)式定义，H[s(t)]表示s(t)的Hilbert变换，Ω表示s(t)中的有效信号的能量分布空间，将待分析信号s(t)对应的复信号记为 的虚部记为sI(t)，则：
[0034]
[0035] 式中Im{f(t)}表示f(t)的虚部；
[0036] 5)根据得到的解析信号，基于极平坦滤波器计算瞬时频率：
[0037] 用具有极平坦频率特性因果的FIR滤波器，并用该滤波器的延时器和微分器组成地震信号瞬时频率估计器，具体讲，用x[n]和y[n]记地震记录的实部和虚部，x[n]＝s(n)，y[n]＝sI(n)，fiN[n]是滤波器长度为N时估计的瞬时频率，具体计算公式如下：
[0038]
[0039] 式中 表示卷积，hN，α(n)表示延时器，dN，α(n)表示微分器，[0040]
[0041]
[0042] n＝1，2，...，N-1，其中N为FIR因果滤波器的长度，α表示延时量，当N为奇数时，α＝(N-1)/2，当N为偶数时，α＝N/2+r，r为一小数；
[0043] 6)根据得到的瞬时频率，进行最佳分辨率瞬时频率分析：
[0044] ①对应于多尺度解析信号计算得到的瞬时频率，对单频带瞬时频率分析的方法为：
[0045] 根据实际地震数据及先验信息，确定尺度的变化范围，在该范围内利用(3)式计算所有的多种分辨率瞬时频率，并将它们按分辨率由低到高顺序排列，利用测井或其它信息在多分辨率瞬时频率剖面上选择最佳分辨率，实现对地震记录的最佳分辨率分析；
[0046] ②对应于在小波域有效信号能量分布空间计算解析信号的方法，在有效信号能量分布空间进行瞬时频率分析，分以下步骤：
[0047] a)利用小波
阈值法或利用有效信号和噪声统计特性差别，确定有效信号能量分布空间；
[0048] b)利用上述第4)步②方法及第5)步计算瞬时频率，利用第2)步中提供的信息，判断得到的瞬时频率是否满足需要，如果满足，输出结果，否则回到a)，重复上述过程，直到得到满足要求的瞬时频率。
[0049] 该方法与Hilbert变换方法相比具有如下优点：1、具有多分辨率特性。2、适用于低信噪比资料。3、对宽频带地震资料可得到高
精度瞬时频率。
附图说明
[0051] 图2a是利用物理小波模拟出的零相位地震子波；
[0052] 图2b是利用物理小波模拟出的混合相位地震子波；
[0053] 图3是采用常规的方法计算瞬时频率时，与精确值的对比图；
[0054] 图4a是叠后地震资料；
[0055] 图4b是基于Hilbert变换的方法计算结果图；
[0056] 图4c是基于本发明的计算结果图。
具体实施方式
[0057] 下面通过附图和
实施例对本发明进一步说明。
[0058] 图1所示是本发明的总流程图，从流程图中可看出，本发明按以下步骤实现：
[0059] 1、首先获取二维或三维的经偏移或叠加处理后的地震资料，这些资料是地震勘探领域技术人员都可以方便获得的。
[0060] 2、根据研究的对象和目的对所获取的二维或三维地震资料进行空间分区，一般按储层深度及储层的空间非均质性情况进行划分。在所划分的区域里能方便获得测井资料或零偏VSP资料、井旁地震记录特别是典型地震道数据、储层的地质构造及其它先验信息。
[0061] 3、通过测井资料或零偏VSP资料，以及井旁地震记录反演地震子波，确定匹配该子波的母物理小波。这一步中本发明采取了特殊的技术手段。分两步来完成：
[0062] 1)反演地震子波
[0063] 利用分区内的测井资料和井旁地震记录，反演出地震子波，或利用垂直地震剖面VSP资料得到地震子波。这项工作在常见的地震资料处理及解释系统中均可完成，不再赘述。
[0064] 2)确定匹配地震子波的母物理小波
[0065] 本发明中，匹配地震子波的母物理小波采用如下表达式：
[0066] g(t；α)＝Aexp[-τ(t-β)2]exp(iσt)+R(t；α) (1)[0067] (1)式中g(t，α)为解析小波，R(t；α)为修正项，其表达式为：
[0068]
[0069] α为一矢量，定义为α＝(A，σ，τ，β)，A为子波的幅度，σ为母小波的调制频率，τ为母小波的能量衰减率，β为母小波的能量延迟时间，它决
定子波的形状。
[0070] (1)式表示的是匹配地震子波的母物理小波，含有四个待定参数，这四个参数确定后物理小波就确定下来。(1)式可模拟各种常见的地震子波，如图2a、2b所示，给出了用(1)式模拟出的零相位及混合相位两
种子波。
[0071] 确定匹配地震子波的母物理小波的方法是：取(1)式的实部和第1)步所得到的地震子波构造如下目标函数：
[0072]
[0073] 式中，real表示取实部，目标函数的值由(1)式中所含的四个待定参数来决定。我们可用最优化方法，通过使(2)式达到极小值，这时对应的四个参数即是待求参数。有了这四个参数，由(1)式可得我们需要的匹配地震子波的物理小波。求解(2)式是个数学问题，利用现成解题方法即可求得。下文为了书写简便，把g(t；α)简记为g(t)。
[0074] 4、在物理小波域计算地震信号对应的解析信号。有两种
算法：
[0075] 1)多尺度解析信号计算
[0076] (1)式定义的g(t)满足：当ω＜0时 (或近似为零)，因此，g(t)为解析小波。另外它还具有如下性质：
[0077] g(t)∈L1(R，dt)∩L2(R，dt)
[0078] 和
[0079]
[0080] 因此满足小波函数的容许条件，这保证了由它进行小波分解，信息不丢失。
[0081] 任给一个地震信号s(t)∈L2(R，dt)(能量有限函数空间)，s(t)相对于g(t)的小波变换定义为：
[0082]
[0083] 这里t∈R，a∈R\{0}，b∈R，t，b都表示时间。函数 表示对其取复共轭。地震信号s(t)经过(3)式就被变换到小波域(也即时频域)。
[0084] 根据解析小波的性质，对任一尺度因子a(a＞0)，S(b，a)的虚部是其实部的Hilbert变换。换句话说，S(b，a)为解析信号。
[0085] 2)在小波域有效信号能量分布空间计算地震信号对应的解析信号[0086] 我们采用下式可计算地震信号s(t)对应的解析信号：
[0087]
[0088] 这里，S(b，a)由(3)式定义，H[s(t)]表示s(t)的Hilbert变换。地震信号s(t)中一般含有噪声，Ω表示s(t)中的有效信号的能量分布空间。(4)式提供了在时间-尺度域计算H[s(t)]的新途径。
[0089] 下文将待分析信号s(t)对应的复信号记为 的虚部记为sI(t)，则：
[0090]
[0091] 式中Im{f(t)}表示f(t)的虚部。
[0092] 5、根据解析信号计算瞬时频率。
[0093] 1)常规方法：
[0094] s(t)对应的瞬时频率计算如下：
[0095]
[0096]
[0097] 其中e(t)和f(t)分别代表s(t)的瞬时幅度和瞬时频率。为了提高抗噪性能，实际应用中常采用如下的阻尼瞬时频率，这个频率是更好的瞬时频率：
[0098]
[0099] 其中 0＜ε＜1称为阻尼因子。
[0100] 2)基于极平坦滤波器(FIFM)计算瞬时频率
[0101] 常规的方法由解析信号计算瞬时频率时采用二阶中心差商来近似微分运算，因此会引起误差，如图3所示。
[0102] 我们用
鸟叫信号来测试其进算精度。鸟叫信号s(t)可表示为：
[0103] s(t)＝sin(c1t2+c2t)
[0104] 图3中的1号线是该信号准确的瞬时频率，2号线是用(7)式计算的瞬时频率，与精确值相比可见，在频率较高时(大于80Hz)时有明显的误差。
[0105] 本发明中采用具有极平坦频率特性因果的FIR滤波器，并用该滤波器的延时器和微分器组成地震信号瞬时频率估计器。具体讲，用x[n]和y[n]记复地震记录的实部和虚部，即x[n]＝s(n)，y[n]＝sI(n)，fiN[n]是滤波器长度为N时估计的瞬时频率。具体计算公式如下：
[0106]
[0107] 式中 表示卷积，fi，N[n]表示归一化瞬时频率。上式中hN，α(n)表示延时器，dN，α(n)表示微分器，其表达式如下：
[0108]
[0109]
[0110] 这里，n＝1，2，...，N-1，其中N为FIR因果滤波器的长度，α表示延时量，当N为奇数时，α＝(N-1)/2；当N为偶数时，α＝N/2+r；r为一小数，实际中N取20或21。
[0111] 6、最佳分辨率瞬时频率分析。此步对应于第4步的两种算法，分别有两种方法：
[0112] 1)对应于多尺度解析信号，进行单频带瞬时频率分析：
[0113] (3)式中定义的小波变换实际上是对原地震信号的
带通滤波，不同的尺度因子a得到的小波变换结果，主频及频宽是不同的，也即其分辨率是不同的。任给一个尺度因子，由(3)式可得到地震信号在该分辨率下的解析信号。把该信号作为输入，利用步骤5可得到该分辨率下的瞬时频率(称这样得到的瞬时频率为单频带瞬时频率)。
[0114] 根据实际地震数据及先验信息，确定尺度的变化范围，在该范围内利用(3)式计算所有的多种分辨率瞬时频率，并将它们按分辨率由低到高顺序排列，利用测井或其它信息在多分辨率瞬时频率剖面上选择最佳分辨率，实现对地震记录的最佳分辨率分析。
[0115] 2)对应于在小波域有效信号能量分布空间计算得到解析信号，进行瞬时频率分析，分以下步骤：
[0116] (a)利用小波阈值法或利用有效信号和噪声统计特性差别，确定有效信号能量分布空间；
[0117] (b)利用上述第4步2)方法及第5步2)方法计算瞬时频率，利用第2步中提供的信息，判断得到的瞬时频率是否满足需要，如果满足，输出结果，否则回到(a)，重复上述过程，直到得到满足要求的瞬时频率。图4a～c是实际资料的算例，并与常用的基于Hilbert变换的方法处理结果的对比图，其中图4a为叠后地震资料，图4b为利用Hilbert变换的方法计算结果，图4c是利用本发明计算的结果，图上标注的A及B处，是气层的顶部和底部，通过对比图4b和图4c可见，本发明明显优于常规方法。
