基于Synchrosqueezing变换的地震资料时频分析和衰减估计方法

摘要

本发明公开了一种基于Synchrosqueezing变换的地震资料时频分析和衰减估计方法，首次将新的时频分析工具Synchrosqueezing变换用于地震资料时频分析，该变换通过对变换域系数的重排，获得一个更加聚集的时频表示，时频分辨率大大提高，将其用于实际地震资料分析和致密砂岩模型含气性检测，能够准确界定储层的位置，指示河道与断层等地质结构，进而有利于进一步的资料解释和井位确定；提出基于Synchrosqueezing变换的地震衰减估计方法，并给出具体实现流程，对某油田致密砂岩储层三维数据体的衰减估计结果和钻井结果有着较好的一致性，该方法可以帮助地质人员指示含气储层，确定钻井位置。

说明书

技术领域

本发明属于地球物理勘探中的信号处理领域，涉及地震资料的时频分析和 地震衰减估计，具体涉及一种基于Synchrosqueezing变换的地震资料时频分析 和衰减估计方法。

背景技术

地震信号是复杂的非平稳信号，时频分析可以用来描述地震信号随时间的 变化规律，刻画地震信号的局部特征，从而反映这些特征所对应的地质结构和 储层信息。因此，时频分析是地震资料处理与解释的重要手段。随着近代信号 处理理论的发展，涌现出大量的时频分析工具和时频分析方法，其中很多也被 应用在地震信号分析中。

地震波经过含气储层后，高频成分衰减迅速，造成地震波在该区域的局部 主频降低，这种异常可以用来指示碳氢储层的分布。然而，该异常在原始地震 资料上并不明显，却可以在频率分解的剖面上得到增强。因此，不同的时频分 析工具被用来检测这种异常。在这一过程中，时频分析工具的时间-频率分辨率 就成为该问题的关键。

Fourier变换作为经典的谱分析工具，在地震信号处理的很多方面得到应用， 比如：频谱分析，噪声压制等，它也是目前地球物理软件中的必备模块。将整 道数据进行Fourier变换以后，可以得到数据的频谱分布，但是无法反映出频率 的局域化特征，即无法刻画频率随时间的变化，而这恰恰是地震信号的一个重 要特征。因此，大量的时间-频率域联合分析工具被用在地震信号分析中，时间 域上的地震信号经过变换以后被展开到时间-频率域，得到的时频图反映了频率 成分随着时间变化的情况。

加窗傅里叶变换(WFT)用移动的窗函数截取信号，然后对每个时间窗里面的 信号做Fourier变换，得到信号的时间-频率分布。Partyka等人将加窗傅里叶变 换用于墨西哥湾某区块三维地震数据体的谱分解，确定薄层的位置和厚度，刻 画地质体的不连续性。Marfurt和Kirlin用加窗傅里叶变换分析薄层调谐效应。

加窗傅里叶变换中窗函数的选择是时间分辨率和频率分辨率的折中，一旦 窗函数选定，时频图的时间分辨率和频率分辨率也随之确定。作为典型的非平 稳信号，地震信号的频率随着时间不断变化，需要进行多尺度分析，即用较宽 的窗分析信号的缓变成分，用较窄的窗分析信号的突变成分，加窗傅里叶变换 的窗宽度固定，显然无法胜任多尺度分析的要求。为了解决这一问题，小波变 换应运而生。最早将小波变换作为时间-尺度局域化工具提出的是地球物理学家 Morlet，此后，物理学家Grossmann和Morlet合作，给出了小波变换的严格定 义。小波变换通过尺度控制窗函数的宽度，用大尺度获得宽窗函数，小尺度获 得窄窗函数，在高低频具有不同的时间-频率分辨率，实现了对信号的多分辨率 分析。由于窗宽通过尺度控制，所以信号经过小波变换得到的是时间-尺度分布， 称作小波尺度谱(scalogram)，尺度和频率没有确定的对应关系，依赖于母小波 的选择。Chakraborty和Okaya将小波变换用于地震信号分析，并与加窗傅里 叶变换作对比，展示了其优越性。Siha等人提出时频连续小波变换的概念，将 小波变换产生的时间-尺度谱转化为时间-频率谱，用于检测和碳氢储层相关的低 频阴影和增强细小地质结构的可视性。高静怀等人研究了地震资料处理中小波 函数的选取问题，提出用匹配地震子波的函数作为基本小波，对地震资料进行 去噪及分频解释的方法。高静怀等人提出三参数小波，并用于薄互层地震资料 分析。Kazemeini等人将小波变换用于德国某区块的三维数据体分析，结果能 够反映河道的沉积和油气的运移。

上述加窗傅里叶变换、小波变换均为线性变换，其时间-频率联合分辨率受 到不确定性原理的限制。Wigner-Ville分布具有最高的时间-频率联合分辨率， 也被广泛用在地震信号分析中。Li和Zheng将Wigner-Ville分布用于地震资料 的谱分解，刻画碳酸盐储层。

近几年，很多学者针对地震信号分析的特点，将已有的时频分析工具进行 推广和发展，或者将新的相空间分析工具用于地震信号。2013年，Lu和Liu将 一个二维反褶积算子作用于加窗傅里叶变换时频谱，得到反褶积加窗傅里叶变 换时频谱，提高了时频分辨率，基于该方法的谱分解被用来刻画塔里木盆地某 区块的洞穴型的碳酸盐储层。2013年，Han和van der Baan利用经验模式分 解(EMD)进行地震信号的时频分析，并用于加拿大某盆地的白垩纪曲流河检测和 细小地质结构的刻画。

但是，传统的时频分析方法，如短时傅里叶变换和小波变换等用于地震信 号时频分析时具有局限性，主要体现在由于时频原子的影响，导致信号在时频 平面上的能量扩散和时频分辨率的下降。

发明内容

本发明目的在于克服现有技术的不足，提供了一种时频分辨率高的基于 Synchrosqueezing变换的地震资料时频分析和衰减估计方法。

为达到上述目的，本发明采用以下技术方案：

基于Synchrosqueezing变换的地震资料时频分析方法：

Synchrosqueezing变换具体包括以下步骤：

步骤a：连续小波变换

信号的小波变换时域和频域分别表示为：

$W_f(a,b)\left(\begin{array}{c}=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{+\infty }f\left(t\right){\psi }^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}\\ =\frac{\sqrt{a}}{2\pi }{\int }_{-\infty }^{+\infty }F\left(\omega \right){\Psi }^{*}\left(\mathrm{a\omega }\right)e^{\mathrm{ib\omega }}\mathrm{d\omega },\end{array}\right)---\left(7\right)$

其中ψ(t)为基本小波，a和b分别为尺度因子和平移因子，F(ω)和Ψ(ω)为 f(t)和ψ(t)的Fourier变换；假设小波函数几乎没有负频率分量，即当ω＜0时， Ψ(ω)≈0；

步骤b：FM解调频率的计算

${\omega }_f(a,b)=\left(\begin{array}{cc}\frac{-i{\partial }_b{W}_f(a,b)}{W_f(a,b)}& \left|W_f(a,b)\right|>0,\\ \infty & \left|W_f(a,b)\right|=0.\end{array}\right)---\left(8\right);$

步骤c：时间-尺度域到时间-频率域的映射

(1)连续形式

$T_f(\omega ,b)={\int }_{\{a:a>0,{\omega }_f(a,b)=\omega ,W_f(a,b)\ne 0\}}{W}_f(a,b)a^{-3/2}\mathrm{da}---\left(9\right)$

通过公式(7)，在时间-尺度域所有和频率ω对应的小波系数进行组合，在时 间-频率域重新将能量“挤压”到频率ω所在的位置；

(2)离散形式

进行数值计算时，需要对公式(7)中的尺度a和频率ω进行离散。离散化后的 尺度记为{a_k}，其中a_k＞a_k-1，尺度间隔为a_k-a_k-1＝(Δa)_k；对频率进行划分，记为 {ω_l}，其中ω_l＞ω_l-1，频率间隔为ω_l-ω_l-1＝Δω；Synchrosqueezing变换的离散形 式表示为：

$T_f({\omega }_l,b)=\underset{a_k:|{\omega }_f(a_k,b)-{\omega }_l|\le \mathrm{\Delta \omega }/2}{\Sigma }{W}_f(a_k,b){a_k}^{-3/2}{\left(\mathrm{\Delta a}\right)}_k,---\left(10\right);$

地震资料时频分析方法具体包括以下步骤：

步骤1：在三维数据体地震剖面中选取典型道，对该道数据进行 Synchrosqueezing变换，找出异常区域对应频率；

步骤2：对整个地震剖面进行Synchrosqueezing变换，提取异常区域对应 频率切片；

步骤3：对整个三维数据体进行Synchrosqueezing变换，得到频率数据体， 然后提取一个沿层切片，供地质人员进行地震资料解释。

基于Synchrosqueezing变换进行地震衰减估计的方法，包括以下步骤：

步骤1：确定目标层范围，对目标层附近的三维地震数据体进行频谱分析， 选取合适的高频f_H和低频f_L；

步骤2：利用Synchrosqueezing变换得到高频的单频数据体T(x,y,t,f_H)和 低频的单频数据体T(x,y,t,f_L)；

步骤3：在目标层上方的层位H_A(x,y)附近，将高频和低频的幅度差异预先 消除；

步骤4：估计目标层附近的衰减。

进一步，所述步骤1中，高频f_H的幅度和低频f_L的幅度大致相同。

进一步，所述步骤3中，利用公式(14)计算修正因子α(x,y)并做平滑；

$\alpha (x,y)=\frac{T(x,y,H_A(x,y),f_L)}{T(x,y,H_A(x,y),f_H)}.---\left(11\right).$

进一步，所述步骤4中，通过公式(15)估计目标层附近的衰减AS(x,y,t)；

AS(x,y,t)＝T(x,y,t,f_L)-α(x,y)T(x,y,t,f_H)(12)。

本发明首次将新的时频分析工具Synchrosqueezing变换用于地震资料时 频分析，该变换通过对变换域系数的重排，获得一个更加聚集的时频表示，时 频分辨率大大提高。将其用于实际地震资料分析和致密砂岩模型含气性检测， 能够准确界定储层的位置，指示河道与断层等地质结构，进而有利于进一步的 资料解释和井位确定。

本发明还提出基于Synchrosqueezing变换的地震衰减估计方法，并给出具 体实现流程，对某油田致密砂岩储层三维数据体的衰减估计结果和钻井结果有 着较好的一致性，该方法可以帮助地质人员指示含气储层，确定钻井位置。

附图说明

图1为余弦信号的小波变换与Synchrosqueezing变换图；

(a)余弦信号；(b)小波变换；(c)Synchrosqueezing变换；

图2为Synchrosqueezing变换示意图；

图3为测试信号和采用不同时频分析方法的结果；

(a)测试信号；(b)加窗傅里叶变换；(c)小波变换；(d)Synchrosqueezing 变换；

图4为地震剖面；

图5为地震道的时频分析；

(a)地震道数据；(b)加窗傅里叶变换；(c)小波变换；(d)Synchrosqueezing 变换；

图6为图2中地震剖面的30Hz频率切片示意图；

(a)小波变换频率切片；(b)Synchrosqueezing变换频率切片；

图7为30Hz数据体的沿层切片示意图；

(a)加窗傅里叶变换沿层切片；(b)Synchrosqueezing变换沿层切片；

图8为某致密砂岩储层三维数据体的地震衰减估计；

(a)工区底图；(b)目的层位；(c)基于Synchrosqueezing变换衰减估计结果 沿目的层的切片。

具体实施方式

以下通过具体实施例和附图对本发明方案做具体说明：

本发明所用工具Synchrosqueezing变换具体包括以下步骤：

步骤a：连续小波变换

信号的小波变换时域和频域分别表示为

$W_f(a,b)\left(\begin{array}{c}=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{+\infty }f\left(t\right){\psi }^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}\\ =\frac{\sqrt{a}}{2\pi }{\int }_{-\infty }^{+\infty }F\left(\omega \right){\Psi }^{*}\left(\mathrm{a\omega }\right)e^{\mathrm{ib\omega }}\mathrm{d\omega },\end{array}\right)---\left(13\right)$

其中ψ(t)为基本小波，a和b分别为尺度因子和平移因子，F(ω)和Ψ(ω)为 f(t)和ψ(t)的Fourier变换。

这里，假设小波函数几乎没有负频率分量，即当ω＜0时，Ψ(ω)≈0，小波的 能量集中在正频率ω＝ω₀附近。如图1(a)所示，以一个单频的余弦信号 h(t)＝Acos(Ωt)，其中A＝1，Ω＝100π为例，其小波变换的结果为

$W_h(a,b)=\frac{A}{2}\sqrt{a}{\Psi }^{*}\left(\mathrm{a\Omega }\right)e^{\mathrm{ib\Omega }}---\left(14\right)$

经过小波变换，余弦信号h(t)的能量在时间-尺度平面围绕尺度a＝ω₀/Ω对应 的直线进行扩散，分布在以该直线为中心的尺度带上，如图1(b)所示。可以看 出，h(t)本来是一个单频信号，由于受到小波函数的影响，在小波变换域能量发 生扩散，得到一个“模糊化”的分布。

本发明试图移除小波函数的影响，从小波变换的结果中将频率Ω“解调” 出来。换言之，尺度带上所分布的能量都是由一个单频成分引起，本发明通过 Synchrosqueezing变换，将尺度带上的能量重新“挤压”到该单频成分所对应 的尺度a＝ω₀/Ω所在的直线上。

步骤b：FM解调频率的计算

为了找出小波系数对应的频率成分，对小波变换的结果进行如下操作：

${\omega }_f(a,b)=\left(\begin{array}{cc}\frac{-i{\partial }_b{W}_f(a,b)}{W_f(a,b)}& \left|W_f(a,b)\right|>0,\\ \infty & \left|W_f(a,b)\right|=0.\end{array}\right)---\left(15\right)$

以余弦信号h(t)为例，将公式(2)代入公式(3)，可以得到当|W_h(a,b)|≠0时，

ω_h(a,b)＝Ω,(16)

也就是说，经过公式(3)的操作，可以从小波变换的结果中重新得到频率Ω。 经过该步骤，可以求出各个小波系数所对应的频率成分。

步骤c：时间-尺度域到时间-频率域的映射

得到各个小波系数对应的频率成分以后，在时间-尺度域进行“挤压”操作， 进行能量重排，将对应相同频率成分的系数组合起来，将时间-尺度平面映射到 时间-频率平面，即(a,b)→(ω(a,b),b)。

对于信号记f_a(t)为该信号对应的解析信号。当基本小波ψ(t)为 解析小波时，

$\left(\begin{array}{c}{\int }_0^{\infty }{W}_f(a,b)a^{-3/2}\mathrm{da}=\frac{1}{2\pi }{\int }_0^{\infty }{\int }_0^{\infty }f\left(\xi \right){\Psi }^{*}\left(\mathrm{a\xi }\right)e^{\mathrm{ib\xi }}{a}^{-1}\mathrm{d\xi da}\\ =\frac{1}{2\pi }{\int }_0^{\infty }F\left(\xi \right)e^{\mathrm{ib\xi }}\mathrm{d\xi }{\int }_0^{\infty }{\Psi }^{*}\left(\mathrm{a\xi }\right)a^{-1}\mathrm{da}\\ =\frac{1}{2}{\int }_0^{\infty }{\Psi }^{*}\left(\xi \right){\xi }^{-1}\mathrm{d\xi }·\frac{1}{2\pi }{\int }_0^{\infty }2F\left(\xi \right)e^{\mathrm{ib\xi }}\mathrm{d\xi }\\ =\frac{1}{2}{\int }_0^{\infty }{\Psi }^{*}\left(\xi \right){\xi }^{-1}\mathrm{d\xi }·f_a\left(b\right),\end{array}\right)---\left(17\right)$

令$C_{\psi }=\frac{1}{2}{\int }_0^{\infty }{\Psi }^{*}\left(\xi \right){\xi }^{-1}\mathrm{d\xi },$可以得到

$f\left(b\right)=\mathrm{Re}\left[C_{\psi }^{-1}{\int }_0^{\infty }{W}_f(a,b)a^{-3/2}\mathrm{da}\right].---\left(18\right)$

公式(6)中Re(·)表示实部。由此，可以将Synchrosqueezing变换定义为以 下几种形式：

(1)连续形式

$T_f(\omega ,b)={\int }_{\{a:a>0,{\omega }_f(a,b)=\omega ,W_f(a,b)\ne 0\}}{W}_f(a,b)a^{-3/2}\mathrm{da}---\left(19\right)$

通过公式(7)，在时间-尺度域，所有和频率ω对应的小波系数进行组合，在 时间-频率域重新将能量“挤压”到频率ω所在的位置。

(2)分布形式

$T_f(\omega ,b)={\int }_{\{a:a>0,W_f(a,b)\ne 0\}}{W}_f(a,b)\delta ({\omega }_f(a,b)-\omega )a^{-3/2}\mathrm{da},---\left(20\right)$

公式(8)采用单位冲激函数来表示Synchrosqueezing变换，Daubechies等 人从分布的角度来解释δ(ω_f(a,b)-ω)，该形式可被称为分布形式。

(3)近似形式

$S_{f,ϵ}^{\delta }(\omega ,b)={\int }_{\{a:a>0,W_f(a,b)>ϵ\}}{W}_f(a,b)\frac{1}{\delta }h\left(\frac{\omega -{\omega }_f(a,b)}{\delta }\right)a^{-3/2}\mathrm{da},---\left(21\right)$

公式(9)中，h(t)为局部支撑的无限光滑函数，即且∫h(t)dt＝1。当δ→0 时，该表示方法用单位脉冲的极限来近似代替 单位冲激函数。注意到当W_f(a,b)很小时，运用公式(3)计算小波系数对应的频率 成分时会出现数值不稳定的情况，因此只在|W_f(a,b)|＞ε时进行计算。本发明将公 式(9)称为具有门限ε和精度δ的Synchrosqueezing变换。

(4)离散形式

进行数值计算时，需要对公式(7)中的尺度a和频率ω进行离散。离散化后的 尺度记为{a_k}，其中a_k＞a_k-1，尺度间隔为a_k-a_k-1＝(Δa)_k。对频率进行划分，记为 {ω_l}，其中ω_l＞ω_l-1，频率间隔为ω_l-ω_l-1＝Δω。Synchrosqueezing变换的离散形 式可以表示为

$T_f({\omega }_l,b)=\underset{a_k:|{\omega }_f(a_k,b)-{\omega }_l|\le \mathrm{\Delta \omega }/2}{\Sigma }{W}_f(a_k,b){a_k}^{-3/2}{\left(\mathrm{\Delta a}\right)}_k,---\left(22\right)$

可见，通过公式(10)，所对应频率成分在ω_l附近(图2中虚线之间的部分) 的小波系数被重新组合，能量在时间-频率平面被“挤压”到ω_l所在位置。

如图1为余弦信号的小波变换与Synchrosqueezing变换，可以看出后者具 有更好的时频聚集性，更加准确地刻画了50Hz余弦信号的时频分布。

如图2为Synchrosqueezing变换示意图，对应频率在虚线之间的小波系 数被重新组合。

图1(c)画出了余弦信号h(t)的Synchrosqueezing变换结果，可见其能量在 时频平面重新聚集到ω＝Ω所在的直线上。

基于Synchrosqueezing变换的地震衰减估计方法

地震子波穿过各向同性衰减介质时，其振幅谱为

S(f)＝S₀(f)e^-αz(23)

其中S₀(f)为地震子波的振幅谱，z为传播距离，α为衰减系数，S(f)为地 震子波传播距离z以后的振幅谱。Aki和Richards假定品质因子Q不随频率变化， 给出了衰减系数和品质因子之间的关系：

$\alpha =\frac{\mathrm{\pi f}}{\mathrm{QV}},---\left(24\right)$

式中V为相速度，f为频率。将(12)代入(11)得

$S\left(f\right)=S_0\left(f\right)e^{-\frac{\mathrm{\pi f}}{\mathrm{QV}}z}.---\left(25\right)$

可见，地震子波中的高频分量衰减比低频分量要快，当地震波经过Q值较 低的碳氢储层时，衰减更加明显。高低频分量不同的衰减特性可以用来指示碳 氢储层，鉴于Synchrosqueezing变换的良好性质，本发明将其用来刻画地震波 经过致密砂岩含气储层的衰减，进行碳氢储层的指示。

研究高低频分量的衰减特性差异，需要选择合适的高频f_H和低频f_L，首先 对目的层附近的地震数据做频谱分析，然后选择f_H和f_L，使其满足以下条件：

(1)振幅谱上高频f_H的幅度和低频f_L的幅度应大致相同；

(2)高频f_H和低频f_L之间的频率间隔应该足够大，以保证高频分量和低频 分量有足够的衰减差异；

(3)高频f_H和低频f_L都应该在地震子波的频带范围内。

记目标层位为H_T(x,y)，利用Synchrosqueezing变换得到的高频分量为 T(x,y,t,f_H)，低频分量为T(x,y,t,f_L)，下面就可以采用高频分量和低频分量的差 来刻画层位H_T(x,y)附近的衰减。由于地下结构和层位介质的复杂性，地震波在 到达目标层之前已经发生了衰减，产生高低频分量的差异，为了刻画目标层对 附近的衰减，需要在目标层位上方将高频和低频的幅度差异预先消除，因此， 记目标层上方的层位为H_A(x,y)，定义如下修正因子：

$\alpha (x,y)=\frac{T(x,y,H_A(x,y),f_L)}{T(x,y,H_A(x,y),f_H)}.---\left(26\right)$

由于地下介质的复杂性，α(x,y)变化剧烈，容易使得计算结果不稳定，因此 在实际应用中需要对α(x,y)进行平滑操作。

由Synchrosqueezing变换得到的目标层位H_T(x,y)附近的衰减定义为

AS(x,y,t)＝T(x,y,t,f_L)-α(x,y)T(x,y,t,f_H)   (27)。

本发明的物质基础是地震数据体，采用的逐道处理办法。

本发明基于Synchrosqueezing变换进行地震资料时频分析具体步骤为：

步骤1：在三维数据体地震剖面中选取典型道，对该道数据进行 Synchrosqueezing变换，找出异常区域对应频率；

步骤2：对整个地震剖面进行Synchrosqueezing变换的结果中，提取频率 切片；

步骤3：对整个三维数据体进行Synchrosqueezing变换，得到频率数据体， 然后提取一个沿层切片。

本发明所提基于Synchrosqueezing变换进行地震衰减估计的方法总结如 下：

步骤1：确定目标层范围，对目标层附近的三维地震数据体进行频谱分析， 选取合适的高频f_H和低频f_L；

步骤2：利用Synchrosqueezing变换得到高频的单频数据体为T(x,y,t,f_H)和 低频的单频数据体为T(x,y,t,f_L)；

步骤3：在目标层上方的层位H_A(x,y)附近，利用公式(14)计算修正因子 α(x,y)并做平滑；

步骤4：通过公式(16)估计目标层附近的衰减AS(x,y,t)。

效果分析

一、数值仿真结果

首先，通过对比Synchrosqueezing变换与加窗傅里叶变换和小波变换对合 成信号进行时频分析的效果。合成信号s(t)如图3(a)所示，它由以下两个分量组 成：

s₁(t)＝sin(3(140πt+30sin(3πt))),(28)

s₂(t)＝sin(3(80πt+20sin(3πt))).(29)

本发明分别对s(t)采用加窗傅里叶变换、小波变换和Synchrosqueezing变 换，如图3(b)-(d)所示。其中加窗傅里叶变换采用128点的Hamming窗，小波 变换采用σ＝6的Morlet小波。尽管三种变换的结果都能够区分开两个分量，但 是加窗傅里叶变换和小波变换的分辨率明显低于Synchrosqueezing变换。由于 窗函数或者小波函数的影响，信号的能量在变换域发生扩散，使得两个分量之 间发生交叠。而Synchrosqueezing变换的结果可以很好地区分两个信号分量， 准确描述频率随时间的变化规律。

如图3，测试信号和采用不同时频分析方法的结果，和传统的加窗傅里叶变 换与小波变换相比，Synchrosqueezing变换可以得到更加精细的时频分布，清 晰地区分开两个分量，准确地刻画出每个分量的频率随时间的变化过程

二、实际地震资料

下面，将不同的时频分析工具用于实际地震数据。图4展示了一个三维数 据体的地震剖面。该剖面共有400道，每道500个采样点，采样间隔为2ms， 其中在第60道到第120道之间，1.22s左右，以及第255道到第305道之间， 1.25s左右有典型的河道特征，分别用椭圆标注。

如图4所示，地震剖面，椭圆指示了河道的位置。

通过抽取第90道进行时频分析，该道经过左边的河道。该道数据及其加窗 傅里叶变换、小波变换和Synchrosqueezing变换如图5所示。从中可以看出由 于地层的吸收作用带来的频率衰减趋势。其中位于1.22s左右的30Hz的异常区 域对应着该河道的位置。由于加窗傅里叶变换和小波变换的能量扩散，其变换 域谱图较为模糊，隐藏了不同信号分量的频率变化特征。而Synchrosqueezing 变换产生了一个更为稀疏的时频分布，有着更高的时频分辨率，揭示了信号频 率的局部变化特征，指示出和储层对应的异常区域。

如图5为地震道的时频分析。Synchrosqueezing变换的结果具有更加稀疏 的分布和更高的时频分辨率，反映出信号不同分量的频率的局部变化特征，清 晰地指示出1.25s左右和河道相关的30Hz的异常区域。

接下来，对整个地震剖面进行时频分析。从小波变换和Synchrosqueezing 变换的结果中提取30Hz的频率切片，如图6所示。Synchrosqueezing变换的 频率切片不像小波变换的频率切片一样平滑，有一些类似于“噪声”的特征， 然而，正是由于它有着更加稀疏的分布，有着较高的时频分辨率，可以更加清 晰地刻画地质体的特征。图中用椭圆指示了河道的位置，可见Synchrosqueezing 变换的结果很清晰地反映出河道的位置和边界。而小波变换由于小波函数的影 响，在时频图上的能量扩散造成不同频率的分量发生交叉、混叠，表现为单频 切片的分辨率下降，图中河道所在位置的能量和两边连成一片，无法分辨河道 的边界。

如图6为图2中地震剖面的30Hz频率切片。Synchrosqueezing变化的结 果具有更高的时频分辨率，准确地指示了河道的位置和边界。

最后，对整个三维数据体进行时频分析。分别采用加窗傅里叶变换和 Synchrosqueezing变换，得到30Hz的频率数据体，然后提取一个沿层切片， 如图7所示。可以看出河道(红色箭头指示部分)和断层(绿色箭头指示部分) 的特征在前者图上变得十分模糊，而由于很高的时频分辨率，Synchrosqueezing 变换的结果能够清晰地反映出河道和细小断层的特征。

如图7为30Hz数据体的沿层切片，其中红色箭头和绿色箭头分别指示了河 道和断层的位置，Synchrosqueezing变换的结果更加清晰地反映了断层和河道 的特征。

将本发明提出的基于Synchrosqueezing变换的地震衰减估计方法用于某 致密砂岩三维实际资料的衰减估计。该致密砂岩储层主要表现为孔隙度小于 10％的低渗透率的砂体储层特征，有效砂体分布比较分散，规模小，连续性差， 波阻抗差异小，给储层预测带来困难。该工区底图如图8(a)所示，目标层位如(b) 所示。该区域分布六口井，自北向南分别用WELL1到WELL6表示，其中东北 方向的WELL3和西南方向的WELL4、WELL5、WELL6为Ⅰ类和Ⅱ类产气井， Ⅰ类为高产井，Ⅱ类次之，在图中用蓝色标注，西北方向的WELL1和WELL2 为Ⅲ类井，即为干井或者产气量极低，在图中用黄色标注。

通过将基于Synchrosqueezing变换的地震衰减估计方法用于处理该区域 的地震资料。经过频谱分析后，分别选用10Hz作为低频，40Hz作为高频，利 用Synchrosqueezing变换得到衰减结果以后，沿(b)所示的目的层进行切片， 如(c)所示。可见WELL3、WELL4、WELL5和WELL6附近表现出强衰减特性， 这与其为Ⅰ类和Ⅱ类产气井是完全符合的，在WELL1和WELL2附近基本上观 察不到明显的衰减，这和它们是Ⅲ类井是吻合的。总之，基于Synchrosqueezing 变换的衰减估计结果和钻井结果有着比较好的一致性，在一定条件下可以作为 碳氢储层的指示因子，帮助地质人员进行井位确定和储层含气量估计。

如图8为某致密砂岩储层三维数据体的地震衰减估计图。

本发明效果总结：

(1)本发明提出基于Synchrosqueezing变换进行地震信号时频分析的方 法，测试信号的例子表明，和传统的时频分析工具，如加窗傅里叶变换、小波 变换等相比，Synchrosqueezing变换通过对变换域系数的重排，可以得到更加 聚集的时频分布，准确地刻画出每个分量的频率随时间的变化过程。

(2)在时频图上，Synchrosqueezing变换产生了一个更为稀疏的时频分布， 有着更高的时频分辨率，揭示了信号频率的局部变化特征，指示出和储层对应 的异常区域；

(3)在地震剖面的单频切片上，Synchrosqueezing变换的结果可以更加清 晰地反映出河道的位置和边界；

(4)在数据体的沿层切片上，Synchrosqueezing变换的结果更加清晰地反 映了断层和河道的特征；

(5)本发明提出基于Synchrosqueezing变换的地震衰减估计方法，给出其 实现流程，并将该方法用于某油田致密砂岩储层三维数据体的衰减估计，估计 结果和钻井结果有较好的一致性，可以作为含气储层的直接指示因子，用来帮 助地质人员进行储层含气量估计和井位确定。