近年来, 压缩感知(compressive sensing, CS)技术^[1]已在地震勘探领域引起极大关注并取得了诸多成功应用^[2-6]。多维数据规则化、多维数据去噪及压缩是地震数据处理中的重要步骤^[7-8], 其中最为重要的是寻找合适的稀疏变换, 这同时也是压缩感知的重要步骤。XU等^[9]于2005年提出了反泄露傅里叶变换(anti-leakage Fourier transform, ALFT)并将其应用于地震数据规则化, 2010年又将其推广到高维数据规则化^[10]。随后几年, ALFT及其各种发展形式在地震数据处理中得到了广泛应用^[11-13]。梁东辉^[14]将非均匀快速傅里叶变换引入反泄露傅里叶变换, 以提高反泄露傅里叶变换在实际地震数据处理中的运算速度。

在实际应用ALFT中, 尤其是在处理高维地震数据时, 通常进行分窗处理以减少计算量。虽然采用分窗处理后, 窗内的数据可以近似认为是少量平面波叠加, 但是, 即使在均匀网格的前提下, 较小的窗仍然会明显降低傅里叶变换域的分辨率, 同时造成频谱泄露, 不能利用傅里叶基函数稀疏表示信号。另外, 采用分窗处理虽然能降低计算量, 但仍大量采用离散傅里叶变换(discrete Fourier transform, DFT), 对于TB规模的地震数据, 计算量仍然巨大。本文针对均匀网格的ALFT, 给出一种频率域快速实现方法, 大幅减少ALFT的运算量, 提高ALFT应用于数据压缩或多维去噪的计算效率。

1 方法原理 1.1 均匀网格下傅里叶变换的频谱泄露现象

假设g[n]为一个64点离散合成信号, 采样步长为1, 采样频率为1Hz。由于实信号可通过Hilbert变换转换为解析信号, 因此这里仅采用复解析信号, 其表达式为:

$ \mathit{\boldsymbol{g}}\left[n \right] = {{\rm{e}}^{{\rm{j2}}\pi {f_0}n}} + 0.15{{\rm{e}}^{{\rm{j2}}\pi ({f_1}n + 0.25)}} $

(1)

式中:f₀和f₁为两个复正弦信号的频率。当f₀和f₁均为1/64Hz的整数倍时(例如为1/16Hz和5/32Hz), 其实部波形及利用快速傅里叶变换(fast Fourier transform, FFT)获得的傅里叶系数模如图 1a和图 1b所示, 图 1b能够真实地反映信号中的两个频率成分。当f₀和f₁都不是1/64Hz的整数倍时, 例如为9/128Hz和21/128Hz时, 其实部波形及利用FFT获得的傅里叶系数模如图 2a和图 2b所示。

虽然图 1a和图 2a均为正弦信号, 但是其振幅谱(图 1b和图 2b)差别非常大。究其原因, 图 1a所示信号的两个分量均能够用64个相互正交的傅里叶基函数中的一个表示(当f₀和f₁为1/64Hz的整数倍时), 因此可以准确稀疏地表示该信号; 而图 2a所包含的两个分量均不能用64个相互正交的傅里叶基函数中的一个表示, 必须采用这64个相互正交的傅里叶基函数来表示图 2a所示信号, 造成单频信号的能量泄露到其它频率, 即频谱泄露的假象。

1.2 均匀网格下的反泄露傅里叶变换

尽管XU等^[9-10]是针对数据规则化中的非均匀网格提出的反泄露傅里叶变换, 但是该算法对于均匀网格同样适用。

假设g[n]为待分析信号, G[f]为利用DFT计算的g[n]傅里叶系数:

$ \mathit{\boldsymbol{G}}\left[f \right] = \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {\mathit{\boldsymbol{g}}\left[n \right]{{\rm{e}}^{{\rm{ - j2}}\pi fn}}} $

(2)

这里, 为了能够更加精准、更加稀疏地表示g[n], 对f的离散化间隔可以小于1/N(可通过时域补零的方法来减小f的离散化间隔)。设$\mathit{\boldsymbol{\tilde g}}\left[n \right]$为重构信号, r[n]为残量信号, 则均匀网格下的ALFT实现步骤可以简述如下:

1) 将$\mathit{\boldsymbol{\tilde g}}\left[n \right]$初始化为0向量, r[n]初始化为g[n], 同时设定停止迭代的能量阈值ε;

2) 利用DFT计算r[n]的傅里叶系数R[f];

3) 找出R[f]中的最大模值及其对应的f_k, 并记录f_k和对应的R[f_k],

$ {f_k} = \mathop {{\rm{arg}}}\limits_f {\rm{max}}\left| {\mathit{\boldsymbol{R}}\left[f \right]} \right| $

(3)

4) 将其对应的时域信号加到重构信号$\mathit{\boldsymbol{\tilde g}}\left[n \right]$中,

$ \mathit{\boldsymbol{\tilde g}}\left[n \right] = \mathit{\boldsymbol{\tilde g}}\left[n \right] + \mathit{\boldsymbol{R}}[{f_k}]{{\rm{e}}^{{\rm{j}}2\pi {f_k}n}} $

(4)

5) 更新r[n],

$ \mathit{\boldsymbol{r}}\left[n \right] = \mathit{\boldsymbol{r}}\left[n \right] - \mathit{\boldsymbol{R}}[{f_k}]{{\rm{e}}^{{\rm{j}}2\pi {f_k}n}} $

(5)

6) 如果r[n]的能量小于设定的能量阈值ε, 则退出迭代过程, 否则返回步骤2)。

通过上述迭代, 由各个单频信号相加的信号可获得稀疏表示(高维情况下通常进行分窗处理, 在分析窗内可以近似认为仅有少数几个平面波混合在一起), 并可恢复出任意其它网格点的信号。

如果一维信号包含有N个采样点, 频率域采用4N个采样点, 则迭代一次的运算复杂度为O[4N·lg(4N)]。假设共迭代M次, 则上述过程的运算复杂度大约为O[4MNlg(4N)]。

1.3 频率域均匀网格ALFT的快速实现方法

观察上述ALFT实现过程可知, 步骤2)为运算量最大的步骤。如果能将该步骤的运算量减少, 则整个算法的运算量会大幅降低。上述过程中之所以需要步骤2), 是因为每次迭代后, 信号r[n]都有所变化, 导致需要利用DFT重新计算R[f]。如果R[f]在每一次迭代之前是已知的, 则该步骤可以省去, 从而大幅减少计算量。

对公式(5)做DFT, 得:

$ \mathit{\boldsymbol{R}}\left[f \right] = \mathit{\boldsymbol{R}}\left[f \right] - \mathit{\boldsymbol{R}}[{f_k}\left] \mathit{\boldsymbol{S}} \right[f-{f_k}] $

(6)

其中,

$ \mathit{\boldsymbol{S}}[f-{f_k}] = \frac{1}{N}{{\rm{e}}^{ - j\pi \left( {f - {f_k}} \right)(N - 1)}}\frac{{{\rm{sin}}[\pi (f-{f_k})N]}}{{{\rm{sin}}[\pi (f-{f_k})]}} $

(7)

为SINC函数。因此, 我们可以在频率域快速实现均匀网格的ALFT, 其步骤简述如下:

1) 将$\mathit{\boldsymbol{\tilde g}}\left[n \right]$初始化为0向量, r[n]初始化为g[n], 同时设定停止迭代的能量阈值ε;

2) 利用DFT计算g[n]的傅里叶系数G[f], 并将R[f]初始化为G[f];

3) 寻找R[f]中模值最大的一个及其对应的f_k, 并记录f_k和对应的R[f_k],

$ {f_k} = \mathop {{\rm{arg}}}\limits_f {\rm{max}}\left| {\mathit{\boldsymbol{R}}\left[f \right]} \right| $

(8)

4) 将其对应的时域信号加到重构信号$\mathit{\boldsymbol{\tilde g}}\left[n \right]$中

$ \mathit{\boldsymbol{\tilde g}}\left[n \right] = \mathit{\boldsymbol{\tilde g}}\left[n \right] + \mathit{\boldsymbol{R}}[{f_k}]{{\rm{e}}^{{\rm{j}}2\pi {f_k}n}} $

(9)

5) 通过减去SINC函数更新R[f],

$ \mathit{\boldsymbol{R}}\left[f \right] = \mathit{\boldsymbol{R}}\left[f \right]\mathit{\boldsymbol{ - R}}[{f_k}\left] \mathit{\boldsymbol{S}} \right[f-{f_k}] $

(10)

6) 如果R[f]的能量小于设定阈值ε, 则退出迭代过程, 否则返回步骤3)。

如果一维信号包含有N个采样点, 频率域采用4N个采样点, 由于上述实现方法避免了每次迭代中使用DFT, 因此迭代一次的运算复杂度为O(4N), 迭代M次的运算复杂度大约为O(4MN)。与前述常规实现方法相比, 本文给出的实现方法运算量大幅下降。

2 应用算例

利用本文ALFT频率域高效实现方法计算了图 2a所示合成信号②的傅里叶系数模, 结果如图 3所示。与图 2b相比, 图 3明显改善了频谱泄露现象, 同时对原始信号的表示更为稀疏。采用35Hz的Ricker子波合成了一个三维数据体, 该数据体包含有8个传播方向不同的平面波, 时间方向采样间隔为0.5ms, 共有512个采样点; x方向和y方向采样间隔均为5m, 每个方向均有64个采样点。整个数据体共有2097152个采样点。在能量阈值ε为2.0×10^-5的限制下, 利用本文ALFT高效实现方法对该三维数据体进行稀疏表示, 仅需要2616个系数(为原始数据采样点数的1/800)就可以基本恢复出三维合成数据。在同样的误差限制下, 采用DFT则需要53796个系数才能基本恢复出三维合成数据, 所需系数的个数大约为本文ALFT方法的20倍。图 4和图 5对比了利用本文ALFT方法和DFT方法压缩三维合成数据并进行恢复的误差。其中图 4a与图 5a为三维合成数据的两个原始剖面, 图 4b和图 5b为利用本文ALFT方法保留2616个系数恢复的两个剖面的误差; 图 4c和图 5c为利用DFT方法保留53796个系数恢复的两个剖面的误差。可以看出, 两种方法恢复的误差基本在一个量级上, 但是采用本文ALFT方法所需要的系数仅为DFT方法的1/20。

3 结束语

本文提出了一种均匀网格反泄露傅里叶变换的频率域高效实现方法。该方法直接在频率域实现迭代搜索过程, 大大提高了ALFT的计算效率。合成数据应用算例表明, 利用本文ALFT高效实现方法可以明显减少频谱泄露现象, 提高稀疏表示性能。

尽管本文是以一维数据为例推导了ALFT的频率域高效实现方法, 但是要将其推广到高维并非难事。下一步的研究工作是针对五维数据去噪及数据压缩进行测试, 并在以后的工作中不断发现问题。

致谢: 感谢西安交通大学陈文超教授提供的建议及帮助。