基于匹配地震子波的物理小波的地震瞬时频率分析方法

摘要

本发明涉及一种基于匹配地震子波的物理小波的地震瞬时频率分析方法，包括如下步骤：1)获取二维或三维的经偏移或叠加处理后的地震资料；2)根据研究对象对所获取的地震资料进行空间分区，区域里获得测井资料或零偏VSP资料、井旁地震记录、储层的地质构造及其它先验信息；3)通过测井资料或零偏VSP资料，以及井旁地震记录反演地震子波，确定匹配该子波的母物理小波；4)在物理小波域计算地震信号对应的解析信号；5)根据得到的解析信号，基于极平坦滤波器计算瞬时频率；6)根据得到的瞬时频率，进行最佳分辨率瞬时频率分析。本发明具有多分辨率特性，适用于低信噪比资料，对宽频带地震资料可得到高精度瞬时频率。

主权项

1.一种基于匹配地震子波的物理小波的地震瞬时频率分析方法，包括如下步骤：1)首先获取二维或三维的经偏移或叠加处理后的地震资料；2)根据研究的对象和目的对所获取的二维或三维地震资料进行空间分区，在所划分的区域里获得测井资料或零偏VSP资料、井旁地震记录、储层的地质构造及其它先验信息；3)通过测井资料或零偏VSP资料，以及井旁地震记录反演地震子波，确定匹配该子波的母物理小波，按如下步骤进行：①反演地震子波利用分区内的测井资料和井旁地震记录，反演出地震子波，或利用零偏VSP资料得到地震子波；②确定匹配地震子波的母物理小波基于如下表达式：g(t；α)＝Aexp[-τ(t-β)²]exp(iσt)+R(t；α)    (1)式中g(t；α)为解析小波，为书写简便，下文把g(t；α)简记为g(t)，R(t；α)为修正项，表达式为：$R(t;\alpha )=-A\sqrt{2}\exp [-{\sigma }^2/\left(8\tau \right)]\exp [-2\tau {(t-\beta )}^2]\exp \left(\mathrm{i\sigma t}\right),$α为一矢量，定义为α＝(A，σ，τ，β)，A为地震子波的幅度，σ为母小波的调制频率，τ为母小波的能量衰减率，β为母小波的能量延迟时间，取(1)式的实部和第1)步所得到的地震资料构造如下目标函数：$\Phi \left(\tilde{\alpha }\right)=\underset{\alpha }{\min }{\int (w\left(t\right)-\mathrm{real}[g(t;\alpha \left]\right)}^2\mathrm{dt}---\left(2\right)$(2)式中，real表示取实部，取(2)式达到极小值时对应的A、σ、τ、β四个参数，代入到(1)式即得到匹配地震子波的物理小波；4)在物理小波域计算地震信号对应的解析信号，采用多尺度解析信号计算和在小波域有效信号能量分布空间计算两种方法：①多尺度解析信号计算(1)式定义的g(t)满足：当ω＜0时因此，g(t)为解析小波，并且g(t)∈L¹(R，dt)∩L²(R，dt)和$\hat{g}\left(\omega \right)\in (R/\{0\},\frac{\mathrm{d\omega }}{\left|\omega \right|}\cap L^2(R/\{0\},\frac{\mathrm{d\omega }}{\left|\omega \right|}),$所以任给一个地震信号s(t)∈L²(R，dt)，s(t)相对于g(t)的小波变换定义为：$S(b,a)=\frac{1}{\sqrt{\left|a\right|}}{\int }_{-\infty }^{\infty }s\left(t\right)\overline{g\left(\frac{t-b}{a}\right)}\mathrm{dt}---\left(3\right)$这里t∈R，a∈R/{0}，b∈R，t，b都表示时间，函数表示对其取复共轭；对任一尺度因子a且a＞0，S(b，a)即为在尺度因子为a且a＞0时s(t)对应的解析信号；②在小波域有效信号能量分布空间计算地震信号对应的解析信号此方法中采用下式计算地震信号s(t)对应的解析信号：$\frac{1}{C_g}{\int }_{\Omega }S(t,a)\frac{\mathrm{da}}{a}=s\left(t\right)+\mathrm{iH}\left[s\left(t\right)\right]---\left(4\right)$这里，S(b，a)由(3)式定义，H[s(t)]表示s(t)的Hilbert变换，Ω表示s(t)中的有效信号的能量分布空间，将待分析信号s(t)对应的复信号记为的虚部记为s₁(t)，则：$s_I\left(t\right)=H\left[s\left(t\right)\right]=\mathrm{Im}\{1/C_g\underset{v}{\int }S(t,a)\mathrm{da}/a\}---\left(5\right)$式中Im{f(t)}表示f(t)的虚部；5)根据得到的解析信号，基于极平坦滤波器计算瞬时频率：用具有极平坦频率特性因果的FIR滤波器，并用该滤波器的延时器和微分器组成地震信号瞬时频率估计器，具体讲，用x[n]和y[n]记地震记录的实部和虚部，x[n]＝s(n)，y[n]＝s_I(n)，f_iN[n]是滤波器长度为N时估计的瞬时频率，具体计算公式如下：$f_{i,N}\left[n\right]=\frac{1}{2\pi }\frac{\left\{h_{N,\alpha }\right[n]\otimes x[n\left]\right\}\left\{d_{N,\alpha }\right[n]\otimes y[n\left]\right\}-\left\{d_{N,\alpha }\right[n]\otimes x[n\left]\right\}\left\{h_{N,\alpha }\right[n]\otimes y[n\left]\right\}}{{\left\{h_{N,\alpha }\right[n]\otimes x[n\left]\right\}}^2+{\left\{h_{N,\alpha }\right[n]\otimes y[n\left]\right\}}^2}---\left(9\right)$式中表示卷积，h_N，α(n)表示延时器，d_N，α(n)表示微分器，$h_{N,\alpha }\left[n\right]=\exp \left[j{w}_0(n-\alpha )\right]\underset{k\ne n}{\overset{N-1}{\underset{k=0}{\Pi }}}\frac{\alpha -k}{n-k},---\left(10\right)$$d_{N,\alpha }\left[n\right]=\exp \left[{\mathrm{jw}}_0(n-\alpha )\right]\frac{{\mathrm{jw}}_0\underset{k\ne n}{\overset{N-1}{\underset{k=0}{\Pi }}}(\alpha -k)-\underset{l\ne n}{\overset{N-1}{\underset{l=0}{\Sigma }}}\underset{k\ne l,n}{\overset{N-1}{\underset{k=0}{\Pi }}}(\alpha -k)}{\underset{k\ne n}{\overset{N-1}{\underset{k=0}{\Pi }}}(n-k)}---\left(11\right)$n＝1，2，...，N-1，其中N为FIR因果滤波器的长度，α表示延时量，当N为奇数时，α＝(N-1)/2，当N为偶数时，α＝N/2+r；r为一小数，N取20或21；6)根据得到的瞬时频率，进行最佳分辨率瞬时频率分析：①对应于多尺度解析信号计算得到的瞬时频率，对单频带瞬时频率分析的方法为：根据实际地震数据及先验信息，确定尺度的变化范围，在该范围内利用(3)式计算所有的多种分辨率瞬时频率，并将它们按分辨率由低到高顺序排列，利用测井或其它信息在多分辨率瞬时频率剖面上选择最佳分辨率，实现对地震记录的最佳分辨率分析；②对应于在小波域有效信号能量分布空间计算解析信号的方法，在有效信号能量分布空间进行瞬时频率分析，分以下步骤：a)利用小波阈值法或利用有效信号和噪声统计特性差别，确定有效信号能量分布空间；b)利用上述第4)步②方法及第5)步计算瞬时频率，利用第2)步中提供的信息，判断得到的瞬时频率是否满足需要，如果满足，输出结果，否则回到a)，重复上述过程，直到得到满足要求的瞬时频率。