地震波在地下传播时会受到介质的吸收而发生能量衰减。该衰减可分为两类:一类是由于介质非均匀性引起的衰减, 包括散射、层状地层引起的透射等; 另一类是由地层吸收引起的衰减特性, 即地层本征衰减^[1]。地层的吸收特征是一种非常重要的油气指示标志, 通常用品质因子Q来表征^[2]。

利用地震资料估计Q值的方法主要有两大类。第一类是时间域的方法, 如时间域解析信号法、上升时间法、子波模拟法和振幅衰减法等; 第二类是频率域的方法, 如谱比法、匹配技术和频率域谱模拟法等。TONN^[3]利用VSP数据比较了7种Q值估计方法, 发现振幅衰减法可靠性最差, 谱比法的效果相对较好。对于时间域的方法, 最主要的问题是地震资料振幅信息不保真及信噪比较低, 因而Q值估算精度不高; 频率域方法的问题是利用傅里叶谱估计一个时窗的地震记录频率谱时, 对窗函数和时窗的长度要求较高。尽管两种方法都各有利弊, 但实际应用中谱比法仍然是目前最常用的Q值估计方法之一^[4]。魏文等^[5]利用S变换将地震道由时域变为时频域, 再由时频域能量分布及中心频率计算出对应时间的中心频率, 进而得到Q值, 但该方法的S变换小波基函数固定, 难以适应实际地震信号; WANG^{[4, 6]}利用时频函数将时频谱转变为一维数据, 最后进行补偿分析得到Q值, 此方法虽然稳定但繁琐; 王宗俊^[7]改进了质心频移法, 推导了子波参数、质点频移与地层Q值之间的关系式, 进而得到地层Q值; 袁焕等^[8]等提出先利用VSP上、下行波的旅行时间联合反演得到地层层速度, 再利用VSP上、下波形反演Q值; 张繁昌等^[9]对地震信号进行子波分解得到匹配子波, 将匹配子波中心频率、中心时间之积与对数振幅进行投影, 利用拟合散点斜率得到Q值, 该方法使得自适应分解得到的子波能够应用谱比法解决薄层和反射波干涉的问题; 范明霏等^[10]研究发现S变换窗函数具有一定局限性, 提出基于广义S变换模极大值的薄储层刻画方法, 提高了广义S变换薄层结构剖面刻画的可行性与可信度; HAO等^[11]推导了衰减地震波的广义S变换表达式, 据此得到新的谱比法Q值反演公式, 消除了窗函数的频率响应, 提高了Q值反演精度; 金子奇等^[12]改进了衰减旅行时层析方法, 提高了衰减旅行时的计算精度, 解决了上覆地层的影响, 使得Q值求取较为准确。

为了进一步提高Q值求取的精度, 本文在传统广义S域求取Q值的基础上, 提出了基于Lucy-Richardson算法-广义S域(以下简称LR-GST)的Q值估算方法。文中先介绍了该方法的基本原理, 接着利用合成复杂信号对比了LR-GST算法与广义S域算法(以下简称GST)应用效果, 分析了两种算法的频率分辨率及稳定性, 再利用粘弹性地层模型进一步试算, 将LR-GST算法提取的Q值与GST算法提取的Q值与理论值进行对比, 分析了两种方法的精度差异, 最后将LR-GST算法与GST算法应用于实际资料的Q值估算, 将两者得到的结果与实际井资料进行对比, 验证本文方法的精度与可信度。

1 方法原理

陈学华等^[13]提出的两参数广义S变换时域表达式为:

$ \begin{array}{l} {S_{{\rm{GST}}}}\left( {\tau ,f} \right) = \int_{ - \infty }^{ + \infty } {x\left( t \right)\frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}\exp \left[ { - \frac{{{\lambda ^2}{f^{2p}}{{\left( {\tau - t} \right)}^2}}}{2}} \right]} \\ \;\;\;\;\;\;\;\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}ft} \right){\rm{d}}t \end{array} $

(1)

式中:f为频率; λ和p是高斯窗函数的尺度调节因子; τ为高斯窗函数在时间轴上滑动时中点对应的时间; x(t)表示地震信号。

同时给出了频率域广义S变换的表达式:

$ \begin{array}{*{20}{c}} {{S_{{\rm{GST}}}}\left( {\tau ,f} \right) = \int_{ - \infty }^{ + \infty } {\left[ {X\left( {f + {f_a}} \right)\exp \left( { - \frac{{2{{\rm{ \mathsf{ π} }}^2}f_a^2}}{{{\lambda ^2}{f^{2p}}}}} \right)} \right]} }\\ {\exp \left( {{\rm{i}}2{\rm{ \mathsf{ π} }}{f_a}\tau } \right){\rm{d}}{f_a}} \end{array} $

(2)

式中:X(f_a)为地震信号x(t)的傅里叶变换; f_a为频移量。

赵伟等^[14]将震源子波的频域表达式表示成高斯函数的形式, 即S₀(f)=exp[-(2πf-2πf_d)²/m], 其中, f_d为震源主频; m控制震源频带宽度, m值越大频带越宽, m值趋于无穷大时, 震源为脉冲震源; f_d和m一般通过对目的层处的子波进行分析得到。根据Futterman理论^[15]和频率域相移原理, 震源子波传播t^*时间后接收的地震子波振幅谱可表示为:

$ \begin{array}{*{20}{c}} {X\left( f \right) = P\exp \left[ {\frac{{ - {{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m}} \right] \cdot }\\ {\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}f{t^ * }} \right)\exp \left( {\frac{{ - {\rm{ \mathsf{ π} }}f{t^ * }}}{Q}} \right)} \end{array} $

(3)

式中:Q表示t^*时间内地震子波所经过地层的Q值; P表示由于几何扩散、透射、反射等引起的与频率无关的能量损失。

将公式(3)代入公式(2), 得:

$ \begin{array}{l} {S_{{\rm{GST}}}}\left( {\tau ,f} \right) = \int_{ - \infty }^{ + \infty } {\left\{ {P \cdot } \right.} \\ \;\;\;\;\;\;\exp \left\{ {\frac{{ - {{\left[ {2{\rm{ \mathsf{ π} }}\left( {f + {f_a}} \right) - 2{\rm{ \mathsf{ π} }}{f_{\rm{d}}}} \right]}^2}}}{m}} \right\} \cdot \\ \;\;\;\;\;\;\exp \left[ { - {\rm{i}}2{\rm{ \mathsf{ π} }}\left( {f + {f_a}} \right){t^*}} \right]\exp \left[ {\frac{{ - {\rm{ \mathsf{ π} }}\left( {f + {f_a}} \right){t^*}}}{Q}} \right] \cdot \\ \;\;\;\;\;\;\left. {\exp \left( { - \frac{{2{{\rm{ \mathsf{ π} }}^2}f_a^2}}{{{\lambda ^2}{f^{2p}}}}} \right)} \right\}\exp \left( {{\rm{i}}2{\rm{ \mathsf{ π} }}{f_a}\tau } \right){\rm{d}}{f_a} \end{array} $

(4)

经过推导和整理得t^*时刻广义S振幅谱^[11]为(推导见附录A):

$ \begin{array}{l} \left| {{S_{{\rm{GST}}}}\left( {{t^*},f} \right)} \right| = P\sqrt {\frac{1}{{\frac{{4{\rm{ \mathsf{ π} }}}}{m} + \frac{{2{\rm{ \mathsf{ π} }}}}{{{\lambda ^2}{f^{2p}}}}}}} \cdot \\ \;\;\;\;\;\;\exp \left\{ { - \left[ {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} + \frac{{{\rm{ \mathsf{ π} }}f{t^*}}}{Q}} \right] + } \right.\\ \;\;\;\;\;\;\left. {\frac{{{{\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{{{t^ * }}}{{2Q}}} \right)}^2}}}{{\frac{4}{m} + \frac{2}{{{\lambda ^2}{f^{2p}}}}}}} \right\} \end{array} $

(5)

高斯窗G(t)的Wigner-Ville分布:

$ \begin{array}{l} {\rm{WV}}{{\rm{D}}_x}\left( {t,f} \right) = \int_{ - \infty }^{ + \infty } {G\left( {t + \frac{\tau }{2}} \right){G^*}\left( {t - \frac{\tau }{2}} \right){e^{ - {\rm{i}}2{\rm{ \mathsf{ π} }}f\tau }}{\rm{d}}\tau } \\ \;\;\;\;\;\;\; = \int_{ - \infty }^{ + \infty } {\frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}\exp \left[ { - \frac{{{\lambda ^2}{f^{2p}}{{\left( {t + \frac{\tau }{2}} \right)}^2}}}{2}} \right]} \cdot \\ \;\;\;\;\;\;\;\frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}\exp \left[ { - \frac{{{\lambda ^2}{f^{2p}}{{\left( {t - \frac{\tau }{2}} \right)}^2}}}{2}} \right]\exp \left( {{\rm{i}}2{\rm{ \mathsf{ π} }}f\tau } \right){\rm{d}}\tau \end{array} $

(6)

化简公式(6)得(推导见附录B):

$ {\rm{WV}}{D_x}\left( {t,f} \right) = \frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}\exp \left( { - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} - {\lambda ^2}{f^{2p}}{t^2}} \right) $

(7)

在此处引入Lucy-Richardson算法, 该算法源于贝叶斯的条件概率定理, 从最大似然估计的角度出发, 假设图像服从Possion分布^[16-17], 通过迭代逼近的方法, 使其得到较为精确的值。具体表达式^[18-19]为:

$ {W_x}\left( {k + 1} \right) = {W_x}\left( k \right)\left[ {W_h^*\frac{{{S_x}}}{{{W_h} \otimes {W_x}\left( k \right)}}} \right] $

(8)

式中:k+1表示迭代次数; S_x为x(t)的广义S振幅谱, 令W_x(0)=S_x; W_h为高斯窗的魏格纳分布, 上标*表示相关算子; ⊗表示褶积算子; W_x(k+1)为基于LR-GST算法经过k+1次迭代后的振幅谱。整理公式(8)得:

$ {W_x}\left( {t,f} \right) = {S_x}\left( {W_h^*\frac{{{S_x}}}{{{W_h} \otimes {S_x}}}} \right) $

(9)

将公式(5)、公式(7)代入公式(9)得(推导见附录C):

$ \begin{array}{l} {W_x}\left( {t,f} \right) = \frac{{\sqrt m {\lambda ^2}{f^{2p}}}}{{\sqrt {4{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}Q + 2{\rm{ \mathsf{ π} }}\lambda mQ} }} \cdot \\ \;\;\;\;\exp \left\{ { - \frac{{2{\rm{ \mathsf{ π} }}f}}{Q}t + \frac{{4{\rm{ \mathsf{ π} }}{\lambda ^2}{f^{2p}}\left( {f - {f_d}} \right)}}{{\left( {2{\lambda ^2}{f^{2p}} + m} \right)Q}}t + } \right.\\ \;\;\;\;\left[ {\frac{{m{\lambda ^2}{f^{2p}}}}{{2\left( {4{{\rm{ \mathsf{ π} }}^2}{f^{2p}} + 2m} \right){Q^2}}}{t^2}} \right] + \frac{{8{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}{{\left( {f - {f_d}} \right)}^2}}}{{2m{\lambda ^2}{f^{2p}} + {m^2}}} - \\ \;\;\;\;\left. {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m}} \right\} \end{array} $

(10)

将不同时刻t₁和t₂到达的反射子波的振幅谱相比并取对数, 有:

$ \begin{array}{*{20}{c}} {\ln \frac{{\left| {{W_x}\left( {{t_2},f} \right)} \right|}}{{\left| {{W_x}\left( {{t_1},f} \right)} \right|}} = - 2{\rm{ \mathsf{ π} }}f\frac{{\Delta t}}{Q} + \frac{{4{\rm{ \mathsf{ π} }}{\lambda ^2}{f^{2p}}\left( {f - {f_d}} \right)}}{{2{\lambda ^2}{f^{2p}} + m}} \cdot }\\ {\frac{{\Delta t}}{Q} + \frac{{m{\lambda ^2}{f^{2p}}}}{{8{{\rm{ \mathsf{ π} }}^2}{f^{2p}} + 4m}}\frac{{t_2^2 - t_1^2}}{{{Q^2}}} + \ln \frac{{{P_2}}}{{{P_1}}}} \end{array} $

(11)

式中:Δt=t₂-t₁; Q为t₁和t₂时间内的品质因子; P₁和P₂分别为两个反射子波与频率无关的能量损失。对地震道进行LR-GST后, 取振幅谱局部极大值^[20]进行对数谱比求取Q值时, 谱比值与频率不再符合线性关系。将公式(11)中的Q和Q²项提出并进行变量替换, 得到简化的Q值反演公式:

$ {S^*}\left( f \right) = \gamma \frac{1}{Q} + \hat P $

(12)

其中, $ \gamma=-2 {\rm{ \mathsf{ π} }} f \Delta t+\frac{4 {\rm{ \mathsf{ π} }} \lambda^{2} f^{2 p}\left(f-f_{d}\right)}{2 \lambda^{2} f^{2 p}+m} \Delta t ; S^{*}(f)=$$\ln \frac{\left|W_{x}\left(t_{2}, f\right)\right|}{\left|W_{x}\left(t_{1}, f\right)\right|} ; \hat{P}=\ln \frac{P_{2}}{P_{1}}$。

依据公式(12), 将对数谱比S^*(f)与γ进行线性拟合, 拟合直线斜率值的倒数即为Q值。

2 方法验证 2.1 合成复杂地震信号模型验证

为对比LR-GST的时频分辨率与时频聚集性, 本文采用一个合成的复杂信号进行验证。该信号包含了低、中、高三个频段, 其中301~510ms包含了20Hz和100Hz两个频率分量的信息, 详细频率信息如表 1所示。该信号的表达式为:

$ x\left( t \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} \cos \left( {40{\rm{ \mathsf{ π} }}t/1024} \right)\\ \cos \left( {80{\rm{ \mathsf{ π} }}t/1024} \right)\\ \cos \left( {200{\rm{ \mathsf{ π} }}t/1024} \right)\\ \cos \left( {200{\rm{ \mathsf{ π} }}t/1024} \right) + \cos \left( {40{\rm{ \mathsf{ π} }}t/1024} \right)\\ \cos \left( {200{\rm{ \mathsf{ π} }}t/1024} \right)\\ \cos \left( {160{\rm{ \mathsf{ π} }}t/1024} \right)\\ \cos \left( {100{\rm{ \mathsf{ π} }}t/1024} \right) \end{array}&\begin{array}{l} 0 \le t \le 160\\ 161 \le t \le 300\\ 301 \le t \le 359\\ 360 \le t \le 450\\ 451 \le t \le 510\\ 511 \le t \le 720\\ 721 \le t \le 1024 \end{array} \end{array}} \right. $

(13)

图 1为合成信号, 图 2和图 3分别为对合成信号进行GST计算和LR-GST计算得到的时频谱。比较图 2、图 3可看出:①图 2箭头处端点效应明显, 出现了部分假频; ②GST时频谱上两种频率成分交界处时频分辨率较差。相对而言, LR-GST时频谱的整体分辨率较高, 不同频率成分交界处边界清晰, 时频聚集性相对优于GST算法。因此, LR-GST算法对于非平稳信号中不同频率分量的区分能力更强。

2.2 粘弹性地层模型验证

图 4a为建立的层状地层模型, 共分5层, 其中第3层是目的层(见图 4红色虚线框), 该层具有较小的速度与Q值; 每层的厚度分别为100、100、80、70和162m。图 4b和图 4c分别为基于该模型利用褶积加衰减、褶积不加衰减算法得到的叠后地震剖面。图 4d为图 4b和图 4c中第10道地震信号对比, 其中实线为褶积加衰减算法获得的信号, 即考虑了非弹性衰减; 虚线为褶积不加衰减算法获得的信号, 即未考虑非弹性衰减。

图 5a和图 5b分别是采用LR-GST算法和GST算法得到的等效Q值剖面。可以看出, 两图的目的层位置都存在异常, 为体现效果, 我们提取了单道数据进行深入分析。图 6a为提取的第10道地震数据, 图 6b为该道数据的LR-GST计算的时频图, 图中清晰地显示了反射层位置的频谱。图 6c中黑色实线为基于LR-GST算法估算的地层等效Q值, 红色实线为基于GST算法获得的地层等效Q值, 绿色为理论等效Q值。分析该图发现:①两种方法都能准确地反映反射界面位置; ②基于LR-GST算法获得的地层等效Q值与理论值更接近, 精度更高。图 6d是S^*(f)-γ线性拟合曲线, 整体上看, 数据点的线性关系比较稳定, 拟合效果较好, 拟合得到的直线斜率即为1/Q。

3 应用实例

为了验证本文提出的基于LR-GST算法Q值估计方法的实用性, 我们选取了实际地震数据进行试验, 并与基于GST算法的Q值估计方法进行了对比。

研究区目的层为二叠系中油组, 整体为由北西向南东抬高的构造背景, 区域上位于盐边构造带, 主要受北东向雁列状断裂带及盐体的塑性活动控制, 在挤压应力背景下形成北东延伸的大型圈闭群, 在此基础上形成多个独立的低幅度圈闭、岩性圈闭及复合圈闭。目的层为砂体, 埋深大于4000m, 厚度大约为9m, 构造幅度一般为5~25m, 其中A井位于砂体的中间, 为高产油气井。图 7为该地区的一条过A井叠后地震剖面, 其中黑色测井曲线为声波速度, 图中圆圈中箭头处为目的层位置(对应低速)。对于该剖面, 我们分别利用基于LR-GST算法(图 8)和基于GST算法(图 9)求取了Q值。从图 8和图 9可以看出, 基于GST提取的低Q值区域边界不明显, 而基于LR-GST算法得到的低Q值异常区边界清晰, 与储层(黑圈位置)对应较好。图 10a为提取的井旁道的地震数据, 图 10b为对应的LR-GST时频谱, 时频谱中时间分辨率和频率分辨率均较高。图 10c为目的层处Q值的切片, 从图中能看出A井处于明显的异常低值区, 表明地震波在该位置衰减明显, 可能与含油有关, 与实际情况相吻合。

4 结论

本文在基于GST算法的Q值估计基础上提出了基于LR-GST算法的Q值估计。与基于GST的Q值估计相比, 基于LR-GST算法的Q值估计方法能得到分辨率较高的时频谱, 使Q值估算的精度和稳定性更好。在实际地震资料处理中, 基于LR-GST算法的Q值估计结果与钻井结果较为吻合, 证明了该方法的稳定性和实用性。需要注意的是, 该方法对参数λ和p的准确性要求较苛刻, 合理选择参数λ和p能够调节LR-GST算法的分辨率, 使Q值的估算结果更准确。

附录A:广义S振幅谱的推导过程

$ \begin{array}{*{20}{c}} {{S_{{\rm{GST}}}}\left( {\tau ,f} \right) = \int_{ - \infty }^{ + \infty } {\left\{ {P \cdot \exp \left\{ {\frac{{ - {{\left[ {2{\rm{ \mathsf{ π} }}\left( {f + {f_a}} \right) - 2{\rm{ \mathsf{ π} }}{f_d}} \right]}^2}}}{m}} \right\} \cdot } \right.} }\\ {\left. {\exp \left[ { - {\rm{i}}2{\rm{ \mathsf{ π} }}\left( {f + {f_a}} \right){t^*}} \right]\exp \left[ {\frac{{ - {\rm{ \mathsf{ π} }}\left( {f + {f_a}} \right){t^*}}}{Q}} \right]\exp \left( { - \frac{{2{{\rm{ \mathsf{ π} }}^2}f_a^2}}{{{\lambda ^2}{f^{2p}}}}} \right)} \right\} \cdot \exp \left( {{\rm{i}}2{\rm{ \mathsf{ π} }}{f_a}\tau } \right){\rm{d}}{f_a}} \end{array} $

(A1)

化简公式(A1)得:

$ \begin{array}{*{20}{c}} {{S_{{\rm{GST}}}}\left( {\tau ,f} \right) = P \cdot \sqrt {\frac{1}{{\frac{{4{\rm{ \mathsf{ π} }}}}{m} + \frac{{2{\rm{ \mathsf{ π} }}}}{{{\lambda ^2}{f^{2p}}}}}}} }\times\\ {\exp \left\{ { - \left[ {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} + \frac{{{\rm{ \mathsf{ π} }}f{t^*}}}{Q}} \right] + \frac{{{{\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{{{t^*}}}{{2Q}}} \right)}^2} - \left( {{t^*} - \tau } \right)}}{{\frac{4}{m} + \frac{2}{{{\lambda ^2}{f^{2p}}}}}}} \right\} \cdot }\\ {\exp \left\{ {2i\left[ {\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{{{t^*}}}{{2Q}}} \right)\left( {{t^*} - \tau } \right) - {\rm{ \mathsf{ π} }}f{t^ * }} \right]} \right\}} \end{array} $

(A2)

则地震波GST的振幅谱为:

$ \begin{array}{*{20}{c}} {\left| {{S_{{\rm{GST}}}}\left( {{t^*},f} \right)} \right| = P \cdot \sqrt {\frac{1}{{\frac{{4{\rm{ \mathsf{ π} }}}}{m} + \frac{{2{\rm{ \mathsf{ π} }}}}{{{\lambda ^2}{f^{2p}}}}}}} }\cdot\\ {\exp \left\{ { - \left[ {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} + \frac{{{\rm{ \mathsf{ π} }}f{t^*}}}{Q}} \right] + \frac{{{{\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{{{t^*}}}{{2Q}}} \right)}^2} - \left( {{t^*} - \tau } \right)}}{{\frac{4}{m} + \frac{2}{{{\lambda ^2}{f^{2p}}}}}}} \right\}} \end{array} $

(A3)

令$\frac{\partial\left|S_{\mathrm{GST}}(\tau, f)\right|}{\partial \tau}\left|S_{\mathrm{GST}}(\tau, f)\right| \frac{2\left(t^{*}-\tau\right)}{\frac{4}{m}+\frac{2}{\lambda^{2} f^{2 p}}}=0 $得τ=t^*。

将τ=t^*代入公式(A3)化简得:

$ \left| {{S_{{\rm{GST}}}}\left( {{t^*},f} \right)} \right| = P \cdot \sqrt {\frac{1}{{\frac{{4{\rm{ \mathsf{ π} }}}}{m} + \frac{{2{\rm{ \mathsf{ π} }}}}{{{\lambda ^2}{f^{2p}}}}}}} \exp \left\{ { - \left[ {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} + \frac{{{\rm{ \mathsf{ π} }}f{t^*}}}{Q}} \right] + \frac{{{{\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{{{t^*}}}{{2Q}}} \right)}^2}}}{{\frac{4}{m} + \frac{2}{{{\lambda ^2}{f^{2p}}}}}}} \right\} $

附录B:高斯窗G(t)的Wigner-Ville分布

$ \begin{array}{*{20}{c}} {{\rm{WV}}{{\rm{D}}_x}\left( {t,f} \right) = \int_{ - \infty }^{ + \infty } {G\left( {t + \frac{\tau }{2}} \right){G^*}\left( {t - \frac{\tau }{2}} \right){e^{ - i2{\rm{ \mathsf{ π} }}f\tau }}{\rm{d}}\tau } }\\ { = \int_{ - \infty }^{ + \infty } {\frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}} \exp \left[ { - \frac{{{\lambda ^2}{f^{2p}}{{\left( {t + \frac{\tau }{2}} \right)}^2}}}{2}} \right]\frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}\exp \left[ { - \frac{{{\lambda ^2}{f^{2p}}{{\left( {t - \frac{\tau }{2}} \right)}^2}}}{2}} \right]\exp \left( { - {\rm{i}}2{\rm{ \mathsf{ π} }}f\tau } \right){\rm{d}}\tau }\\ { = \int_{ - \infty }^{ + \infty } {\frac{{{\lambda ^2}{f^{2p}}}}{{2{\rm{ \mathsf{ π} }}}}\exp \left[ {\left( { - \frac{{{\lambda ^2}{f^{2p}}}}{2}} \right) \cdot \left( {2{t^2} + \frac{{{\tau ^2}}}{2} - {\rm{i}}2{\rm{ \mathsf{ π} }}f\tau } \right)} \right]{\rm{d}}\tau } }\\ { = \frac{{{\lambda ^2}{f^{2p}}}}{{2{\rm{ \mathsf{ π} }}}} \times \exp \left( { - {\lambda ^2}{f^{2p}}{t^2}} \right) \times \int_{ - \infty }^{ + \infty } {\exp } \left( { - \frac{{{\lambda ^2}{f^{2p}}{\tau ^2}}}{4} - {\rm{i}}2{\rm{ \mathsf{ π} }}f\tau } \right){\rm{d}}\tau }\\ { = \frac{{{\lambda ^2}{f^{2p}}}}{{2{\rm{ \mathsf{ π} }}}} \times \exp \left( { - {\lambda ^2}{f^{2p}}{t^2}} \right) \times \int_{ - \infty }^{ + \infty } {\exp } \left[ { - \frac{1}{4}{\lambda ^2}{f^{2p}}\left( {{\tau ^2} + \frac{{8{\rm{ \mathsf{ π} }}if\tau }}{{{\lambda ^2}{f^{2p}}}}} \right)} \right]{\rm{d}}\tau }\\ { = \frac{{{\lambda ^2}{f^{2p}}}}{{2{\rm{ \mathsf{ π} }}}} \times \exp \left( { - {\lambda ^2}{f^{2p}}{t^2}} \right) \times \int_{ - \infty }^{ + \infty } {\exp } \left\{ { - \frac{1}{4}{\lambda ^2}{f^{2p}}\left[ {{{\left( {\tau + \frac{{4{\rm{ \mathsf{ π} }}if\tau }}{{{\lambda ^2}{f^{2p}}}}} \right)}^2}} \right] - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}}} \right\}{\rm{d}}\tau } \end{array} $

(B1)

利用$\int_{-\infty}^{+\infty} \mathrm{e}^{-a x^{2}} d x=\sqrt{\frac{{\rm{ \mathsf{ π} }}}{a}} $, 则公式(B1)可写成:

$ {\rm{WV}}{{D}_x}\left( {t,f} \right) = \frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {\rm{ \mathsf{ π} }} }} \times \exp \left( { - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} - {\lambda ^2}{f^{2p}}{t^2}} \right) $

附录C:LR-GST算法推导过程

先计算W_h⊗S_x, 暂时不计算P分离, 将公式(5)和公式(7)代入得:

$ \begin{array}{*{20}{c}} {{W_h} \otimes {S_x} = \int_{ - \infty }^{ + \infty } {{S_x}} \left( \tau \right){W_h}\left( {t - \tau } \right){\rm{d}}\tau = }\\ \int_{ - \infty }^{ + \infty }{\sqrt {\frac{1}{{\frac{{4{\rm{ \mathsf{ π} }}}}{m} + \frac{{2{\rm{ \mathsf{ π} }}}}{{{\lambda ^2}{f^{2p}}}}}}}\exp \left\{ { - \left[ {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} + \frac{{{\rm{ \mathsf{ π} }}f\tau }}{Q}} \right] + \frac{{{{\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{\tau }{{2Q}}} \right)}^2}}}{{\frac{4}{m} + \frac{2}{{{\lambda ^2}{f^{2p}}}}}}} \right\} \cdot }\\ {\frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {\rm{ \mathsf{ π} }} }}\exp \left( { - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} - {\lambda ^2}{f^{2p}}{{\left( {t - \tau } \right)}^2}} \right){\rm{d}}\tau }\\ { = \frac{{\sqrt m {\lambda ^2}{f^{2p}}}}{{\sqrt {4{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}} + 2{\rm{ \mathsf{ π} }}\lambda m} }} \cdot \exp \left[ { - \frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}}} \right] \cdot }\\ {\int_{ - \infty }^{ + \infty } {\exp \left[ { - \frac{{{\rm{ \mathsf{ π} }}f\tau }}{Q} + \frac{{{{\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{\tau }{{2Q}}} \right)}^2}}}{{\frac{4}{m} + \frac{2}{{{\lambda ^2}{f^{2p}}}}}} - {\lambda ^2}{f^{2p}}{{\left( {t - \tau } \right)}^2}} \right]{\rm{d}}\tau } }\\ { = \frac{{\sqrt m {\lambda ^2}{f^{2p}}}}{{\sqrt {4{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}} + 2{\rm{ \mathsf{ π} }}\lambda m} }} \cdot \exp \left[ { - \frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} + \frac{{8{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}{{\left( {f - {f_d}} \right)}^2}}}{{2m{\lambda ^2}{f^{2p}} + m}}} \right] \cdot }\\ {\int_{ - \infty }^{ + \infty } {\exp } \left\{ {\left[ {\frac{{m{\lambda ^2}{f^{2p}}}}{{\left( {16{\lambda ^2}{f^{2p}} + 8m} \right)Q}} - {\lambda ^2}{f^{2p}}} \right]{\tau ^2} + \left[ {\frac{{2{\rm{ \mathsf{ π} }}{\lambda ^2}{f^{2p}}\left( {f - {f_d}} \right)}}{{2{\lambda ^2}{f^{2p}} + m}} - \frac{{{\rm{ \mathsf{ π} }}f}}{Q} - 2{\lambda ^2}{f^{2p}}t} \right]\tau - {\lambda ^2}{f^{2p}}t} \right\}{\rm{d}}\tau } \end{array} $

(C1)

令$ A=\frac{m \lambda^{2} f^{2 p}}{\left(16 \lambda^{2} f^{2 p}+8 m\right) Q}-\lambda^{2} f^{2 p}$, $B=\frac{2 {\rm{ \mathsf{ π} }} \lambda^{2} f^{2 p}\left(f-f_{d}\right)}{2 \lambda^{2} f^{2 p}+m}-\frac{{\rm{ \mathsf{ π} }} f}{Q}-2 \lambda^{2} f^{2 p} t, C=\lambda^{2} f^{2 p} t $, 则上式可化为:

$ \begin{array}{*{20}{c}} {{W_h} \otimes {S_x} = \frac{{\sqrt m {\lambda ^2}{f^{2p}}}}{{\sqrt {4{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}} + 2{\rm{ \mathsf{ π} }}\lambda m} }} \cdot \exp \left[ { - \frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} + \frac{{8{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}{{\left( {f - {f_d}} \right)}^2}}}{{2m{\lambda ^2}{f^{2p}} + m}}} \right] \cdot }\\ {\int_{ - \infty }^{ + \infty } {\left\{ {A\left[ {{{\left( {\tau + \frac{B}{{2A}}} \right)}^2} - \frac{{{B^2}}}{{4{A^2}}}} \right] - C} \right\}} {\rm{d}}\tau } \end{array} $

(C2)

化简公式(C2)为:

$ \begin{array}{*{20}{c}} {{W_h} \otimes {S_x} = \sqrt {\frac{{mQe{\lambda ^2}{f^{2p}}}}{{4{\rm{ \mathsf{ π} }}{\lambda ^2}{f^{2p}} + 2\lambda m}}} \cdot }\\ {\exp \left\{ { - \frac{{{{\left[ {2{\rm{ \mathsf{ π} }}\left( {f - {f_d}} \right)} \right]}^2}}}{m} - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} + \frac{{8{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}{{\left( {f - {f_d}} \right)}^2}}}{{2m{\lambda ^2}{f^{2p}} + m}} - {\lambda ^2}{f^{2p}}{t^2}} \right\}} \end{array} $

(C3)

将公式(C3)代入到公式(9)中:

$ \begin{array}{*{20}{c}} {{W_x}\left( t \right) = {S_x}\left( {{W_h}*\frac{{{S_x}}}{{{W_h} \otimes {S_x}}}} \right) = \sqrt {\frac{1}{{\frac{{4{\rm{ \mathsf{ π} }}}}{m} + \frac{{2{\rm{ \mathsf{ π} }}}}{{{\lambda ^2}}}{f^{2p}}}}} \exp \left\{ { - \left[ {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} + \frac{{{\rm{ \mathsf{ π} }} ft}}{Q}} \right] + \left.\frac{\left(\frac{4 {\rm{ \mathsf{ π} }} f}{m}-\frac{4 {\rm{ \mathsf{ π} }} f_{d}}{m}+\frac{t}{2 Q}\right)^{2}}{\frac{4}{m}+\frac{2}{\lambda^{2} f^{2 p}}}\right\}} \right.}\cdot\\ {\left\{ {\left[ {\frac{{\lambda {{\left| f \right|}^p}}}{{\sqrt {\rm{ \mathsf{ π} }} }} \times \exp \left( { - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} - {\lambda ^2}{f^{2p}}{t^2}} \right)} \right] \times } \right.}\\ {\sqrt {\frac{1}{{\frac{{{\rm{4 \mathsf{ π} }}}}{2} + \frac{{2{\rm{ \mathsf{ π} }}}}{{{\lambda ^2}{f^{2p}}}}}}} \exp \left\{ { - \left[ {\frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m} + \frac{{{\rm{ \mathsf{ π} }}ft}}{Q}} \right] + \frac{{{{\left( {\frac{{4{\rm{ \mathsf{ π} }}f}}{m} - \frac{{4{\rm{ \mathsf{ π} }}{f_d}}}{m} + \frac{t}{{2Q}}} \right)}^2}}}{{\frac{4}{m} + \frac{2}{{{\lambda ^2}{f^{2p}}}}}}} \right\}/\sqrt {\frac{{mQe{\lambda ^2}{f^{2p}}}}{{4\pi {\lambda ^2}{f^{2p}} + 2\lambda m}}} }\cdot\\ \left.{{\rm{exp}} { - \frac{{{{\left[ {2{\rm{ \mathsf{ π} }}\left( {f - {f_d}} \right)} \right]}^2}}}{m} - \frac{{4{{\rm{ \mathsf{ π} }}^2}{f^2}}}{{{\lambda ^2}{f^{2p}}}} + \frac{{8{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}{{\left( {f - {f_d}} \right)}^2}}}{{2m{\lambda ^2}{f^{2p}} + m}} - {\lambda ^2}{f^{2p}}{t^2}} }\right\} \end{array} $

(C4)

上式省去高阶化简为:

$ \begin{array}{*{20}{c}} {{W_x}\left( {t,f} \right) = \frac{{\sqrt m {\lambda ^2}{f^{2p}}}}{{\sqrt {4{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}Q + 2{\rm{ \mathsf{ π} }}\lambda mQ} }} \cdot \exp \left\{ { - \frac{{2{\rm{ \mathsf{ π} }}f}}{Q}t + \frac{{4{\rm{ \mathsf{ π} }}{\lambda ^2}{f^{2p}}\left( {f - {f_d}} \right)}}{{\left( {2{\lambda ^2}{f^{2p}} + m} \right)Q}}t + } \right.}\\ {\left. {\left[ {\frac{{m{\lambda ^2}{f^{2p}}}}{{2\left( {4{{\rm{ \mathsf{ π} }}^2}{f^{2p}} + 2m} \right){Q^2}}}{t^2}} \right] + \frac{{8{{\rm{ \mathsf{ π} }}^2}{\lambda ^2}{f^{2p}}{{\left( {f - {f_d}} \right)}^2}}}{{2m{\lambda ^2}{f^{2p}} + {m^2}}} - \frac{{{{\left( {2{\rm{ \mathsf{ π} }}f - 2{\rm{ \mathsf{ π} }}{f_d}} \right)}^2}}}{m}} \right\}} \end{array} $