由于受到地表障碍、坏道、禁区和成本等因素的制约, 野外地震数据在空间上经常不满足采样定理, 不完备的地震数据会影响波动方程偏移^[1]、表面多次波消除^[2]和AVO分析^[3]等数据处理的效果, 因此须重建地震数据以获得完备的数据。基于稀疏变换的重建方法计算效率高, 数值效果稳定, 因此是重建方法重要的研究方向。该方法利用信号在变换域的稀疏性为约束建立反演模型, 并且采用稀疏反演方法求解。稀疏变换的性质是决定这类方法质量的首要因素。地震信号在变换域中的表达越稀疏, 则越有利于提高重建质量。在地震信号处理领域经典的稀疏变换有傅里叶变换^[4]和Radon类变换^[5-6]。最近曲波变换^[7]等多尺度变换被用于地震数据重建且展示了很好的数值效果。多尺度变换比单一尺度的变换能够展示更多的信息, 因此是信号分析的有力工具。曲波变换打破了传统傅里叶变换和小波变换只能表示单一尺度或者不能稀疏表示多尺度信号的缺点。对于地震数据有非常稀疏的表达, 在地震数据重建^[8-10]、去噪^[11-12]、多次波消除^[13-14]、反褶积^[15]、偏移^[16]和反射数据的增强^[17]等问题中都有应用并且展示了很好的数值效果。曲波变换对曲线或曲面形状信号进行多方向和各项异性的表示, 为边缘不连续的目标提供了几乎最优的表达^[18-21]。虽然曲波变换有多尺度性、多方向性以及无须对数据划分时间窗口等优点, 但是曲波变换的冗余度(冗余度指的是变换域系数的维数与原始信号系数维数之比)很大。在二维情况下的原始曲波变换的冗余度约为8, 三维数据变换的冗余度为24~32, 需要很多的内存和计算时间, 影响大规模数据计算的效率。笔者将一种低冗余三维曲波变换用于三维地震数据重建, 这种新的低冗余曲波变换^[22]比原始曲波变换的冗余度有很大的下降, 因此称之为低冗余度曲波变换, 在求解时采用一种解基于分析的1范数正则化模型的迭代软阈值算法, 对模拟和真实地震数据的试验证明方法的有效性。

曲波^[23]变换有两种快速实现方式, 分别是基于USFFT(unequally-space fast transform)的快速算法和基于wrapping的快速算法。由于基于wrapping的方法更加简单和容易执行, 本文中采用基于wrapping的曲波变换的实现过程。基于wrapping的三维曲波^[23]变换包括一个低通逼近子带和一些包含曲波系数的曲波基子带, 这些曲波系数由尺度、位置和方向索引^[22]。该变换主要由两部分组成:多尺度分解和角度分解。首先, 假设输入的三维数据为N=(N_x, N_y, N_z), 其中，N_x、N_y、N_z分别表示x、y、z三个方向上的维数, 在傅里叶域通过三维Meyer小波将其分解, 得到尺度为N, N/2, …, N/2^J的以原点为中心的立方体, 其中J为尺度的个数。第二步, 将这些立方体和角度窗口相乘, 分割成各向异性的楔形体, 这些楔形体符合抛物线尺度准则, 将楔形体包装成立方体后采用三维傅里叶逆变换得到曲波系数。三维低冗余度曲波变换的具体实现方式如下:①输入N=(N_x, N_y, N_z)的三维数据体, 定义划分的尺度个数J, 在最小尺度的每个面上划分的角度个数N_a。②多尺度分解。在傅里叶域应用三维Meyer小波变换, 得到尺度为N, N/2, …, N/2^J的立方体。③角度分解。在每一个尺度的每个方向上, 在傅里叶域将小波立方体和角度窗相乘得到楔形体; 将这些楔形体包装成中心在原点的平行六面体; 对平行六面体应用三维傅里叶逆变换得到曲波系数。

在原始曲波变换中, 冗余度主要来源于角度分解。此外, Meyer小波的应用方法产生了额外冗余性, 而低冗余曲波变换中Meyer小波没有产生额外冗余性^[22]。新的曲波变换和原始曲波变换在实现上存在许多不同之处, 但是其能够降低额外冗余度的主要原因是Meyer小波变换的应用方法不同。

使用Curvelab软件包^[20], 曲波变换的额外冗余度主要来自放射窗口的实现, 尤其是最小尺度的窗口。以一维Meyer小波变换为例, 假设ψ_j为尺度j∈{0, …, J-1}的Meyer小波, ϕ_J-1为最大尺度上的尺度函数, 用${M_j} = {{\hat \psi }_j} = {2^{ -3j2}}\hat \psi ({2^{ -j}})$和${M_J} = {{\hat \phi }_{J -1}} = {2^{ - 3\left( {J - 1} \right)2}}\hat \phi ({2^{ - 2(J - 1)}})$分别代表它们的傅里叶变换。Meyer小波${\hat \psi }$ (ξ)定义在傅里叶域, 如下所示:

其中ν从0到1逐渐变化并且满足ν((x)+ν(1-x)=1。Meyer尺度函数定义为

最小尺度的小波在傅里叶域的支撑为[-2/3, -1/6]∪[1/6, 2/3], 因此超出了Shannon带。原始曲波变换隐含的假设周期边界条件, 计算一个周期信号的小波变换等价于将信号在周期小波基上分解。最小尺度小波超出Shannon带的部分被其以|ξ|=1/2为中心沿着竖轴做镜像代替, 数据的支撑是原始数据的4/3倍, 对于三维数据来说, 产生的额外的冗余度为(4/3)³。在低冗余度曲波变换中, 尺度函数和小波函数的支撑首先被缩小为原来的3/4, 然后为了保持单位划分的一致性, 对最小尺度的小波进行修改来压制其下降的尾部, 小波在[-1/2, -1/4]∪[1/4, 1/2]之间变为一个固定值。这样对于三维Meyer小波来说, 在傅里叶域内没有增加额外的冗余度。

对于原始曲波变换, 三维数据的冗余度为24~32, 而低冗余度曲波变换的冗余度约为10^[22]。除了在最小尺度上保持方向选择性和低冗余性, 该低冗余曲波变换保持了等距性, 并且能够快速精确重建, 即该低冗余度曲波变换是一个Parseval紧框架, C^*C=I, 其中C为曲波分解算子, C^*为其伴随。C^*也是逆算子, 因此该变换具有快速的重构算法。关于低冗余度曲波变换的详细实现见参考文献[22]。

地震采集过程是地震信号的采样过程, 因此地震重建问题可以根据压缩感知理论求解^[7]。压缩感知理论证明了可以根据少量的测量数据高精度地重建原始数据, 从而降低地震勘探的成本。地震采集过程可以表示为如下数学模型:

式中, R为采样矩阵; x为完整的地震数据; ε为可加性噪声; d为采样的数据。由于数据d是不完备且不满足采样定理的, 理论上存在无穷多x满足式(3)。根据压缩感知理论, 假设地震数据在某个变换域中是稀疏的, 即s=Ψx是稀疏的或者可压缩的, 其中Ψ是一个稀疏变换算子, 令Ψ为上述的低冗余曲波变换, 由于其是紧框架, 其共轭转置也是其逆, 即Ψ^*=Ψ^-1, 则公式(3)变为

其中Ψ^*为Ψ的共轭转置。基于s的稀疏性假设, 可以通过解以下优化问题求出s:

其中σ≥0是对噪声能量的估计, 该问题称为基追踪问题^[24]。当σ=0时, 上述模型针对高信噪比数据或者不含噪声数据的重建; 当σ>0时, 上述模型可以求解含随机噪声数据的重建, 即可以实现同时重建和去噪, 因此上述模型既可以求解不含噪声数据的重建, 也可以求解含随机噪声数据的重建。当s被解出后, 对于低冗余曲波变换来说, x可以由x=Ψ^*s求出。由于问题(5)是一个约束优化问题, 对于该问题的求解方法很少^[10], 一般求解其等价的无约束优化问题

其中λ是一个正则参数, 用来平衡残差的2范数和解的1范数, 当λ选择合适时, 问题(5)和(6)的解相等。为了使问题(6)的解与数据模型相匹配, 对于不含噪声数据的重建, 只须选择λ很小, 接近于0即可; 对于含随机噪声的重建, 若σ已知, 则可以通过偏差原则, 广义偏差原则^[25]等方法求解与σ匹配的正则参数λ, 若σ已知, 则可以采用类似于广义交叉检验原则(GCV)或者L曲线准则的方法求解合适的λ。由于广义交叉检验原则(GCV)或者L曲线准则不能保证一定获得可靠的λ^[26], 因此λ常常通过人工经验选取来获得。

地震数据重建属于大规模计算问题, 因此需要快速稳定的计算方法。对于上述1范数正则化的问题(6)的求解, 不能采用经典的方法如共轭梯度法或者牛顿法^[27-28], 传统的内点法等由于用到矩阵求逆, 因此计算速度无法与经典方法相比, 须研究解1范数正则化模型的快速算法。阈值类方法^[29]是针对1范数正则化模型求解的一类重要方法, 该方法根据1范数可分的性质得到, 很多解问题(6)的方法, 如不动点法、逼近点法等都含有软阈值运算。最近研究发现采用如下的模型的数值效果比基于模型(6)的效果要好:

该模型称为基于分析的模型, 以时间空间域的数据x为变量, 问题(6)被称为基于合成的模型, 以变换域中的系数s为变量。虽然两者之间可以相互推导, 但是基于分析的模型(7)的数值效果要稍好于前者。给出一个解问题(7)的简单易行的阈值算法, 该算法以x为变量进行迭代求解。算法如下:

(1) 输入采样的数据d, 采样矩阵R, 稀疏变换Ψ, 最大迭代步数N, 控制停机的参数ω和最小阈值τ_min, 令k=0。

(2) 令初始解为x⁰=0, 定义残差为r⁰=d, 定义初始阈值τ₀=max||Ψ^*d||_∞。

(3) 计算预测残差${{\mathit{\boldsymbol{\hat r}}}^{k + 1}} = \mathit{\boldsymbol{d}} - \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{x}}^k}$, 通过软阈值运算升级变换域系数α^k+1=T_τk(Ψ^*(x^k+κ_k${{\mathit{\boldsymbol{\hat r}}}^k}$)), 其中T_τk(·)为软阈值运算, κ_k为步长, 升级迭代解x^k+1=Ψα^k+1。

(4) 升级残差r^k+1=d-Rx^k+1, 减小阈值为τ_k+1, 令k=k+1。如果k≥N或者||r^k||₂≤ω, 转步骤(5), 否则转步骤(3)。

(5) 输出最终解x^final=x^k。

如果将算法中软阈值运算变为硬阈值或者其他改进的阈值运算, 上述算法也能得到稀疏解。本文中步长κ_k固定为1, 同时还有更加先进的步长选择方法如BB法、回溯法等方法, 步长的选择对算法有较大的影响。阈值的下降可以是线性的, 也可以是指数下降, 本文中采用指数下降τ_k+1=τ₀e^c(k-1)/(N-1), 其中c=ln(τ₀/τ_min), 实际上阈值的下降方式对解的改善不大, 最小阈值τ_min的大小决定着解的质量, 因为τ_min与噪声的能量相关。当噪声很小时, 选择τ_min很小, 当噪声较大时τ_min应适当调节。由于噪声的能量实际上也是很难估计的, 本文中采用尝试法选择τ_min。本文中算法与固定阈值的阈值类算法不同, 实际上是带有连续化策略的阈值法, 其中令阈值逐渐下降就是优化算法中的连续化策略, 这比固定阈值参数的算法需要更少的迭代次数和更好的数值效果, 因此本算法实际上由两层组成, 外层控制阈值的下降, 内层只用1次迭代, 已经证明本算法比固定最小阈值的方法数值表现要改进很多^[26]。

第一个例子用来验证基于低冗余度曲波变换数据重建的计算效率和数值效果。为了公平地比较两种变换, 所用的重建算法采用相同的参数。数据由有限差分法模拟得到, 假设地下只有两层速度, 合成数据在inline方向和crossline方向的检波器个数都为64, 道间距为12 m, 时间采样个数为64, 时间采样间隔为4 ms。完整采样的数据见图 1(a), 20%随机采样的结果见图 1(b)。基于低冗余度曲波变换的重建结果见图 1(c), 图 1(d)是原始数据与图 1(c)的误差。更多不同采样比例的重建结果详见表 1, 其中的信噪比的定义为${R_{{\rm{SN}}}} = 10{\rm{lo}}{{\rm{g}}_{10}}\frac{{\parallel {\mathit{\boldsymbol{x}}_{{\rm{orig}}}}\parallel _2^2}}{{\parallel {\mathit{\boldsymbol{x}}_{{\rm{orig}}}} - {\mathit{\boldsymbol{x}}_{{\rm{rest}}}}\parallel _2^2}}$, x_orig代表完整的数据, x_rest表示重建的数据。表 1还列出了基于原始曲波变换的重建时间和信噪比, “Candes”表示原始曲波变换, “Woiselle”表示低冗余度曲波变换。由表 1可知, 基于低冗余度曲波变换的计算时间约为原始曲波变换计算时间的1/4, 此外这两种变换重建的结果是相当的。

图 1

图 1 三维合成数据重建实验 Fig.1 3D synthetic data reconstruction experiment

表 1

表 1 合成数据不同采样比例的重建结果比较 Table 1 Comparison of interpolation results for synthetic data with different sampling ratios 参数 采样率50% 采样率20% 采样率10% 采样率5% Candes Woiselle Candes Woiselle Candes Woiselle Candes Woiselle 计算时间t/s 1066 176 756 190 759 176 1082 271 信噪比/dB 20.3365 20.2013 13.5222 17.1638 8.7660 10.3631 4.1165 7.6232

表 1 合成数据不同采样比例的重建结果比较 Table 1 Comparison of interpolation results for synthetic data with different sampling ratios

用一个海洋地震数据验证基于低冗余曲波变换的真实数据的重建质量。该海洋数据为512×128×128, 即时间轴上的采样个数为512, 时间采样率2 ms, 在inline和crossline方向上的采样个数为128, 相邻地震道间隔为25 m。图 2给出海洋数据的一个切片, 图 3是随机采样的数据, 采样的道数为总道数的50%。图 4是采用低冗余曲波变换和迭代软阈值法的重建结果, 其信噪比为16.5596 dB。为了更加清晰的显示重建的效果, 将其中一道原始数据和其相应的重建道进行对比, 如图 5所示, 重建数据和原始数据吻合得很好。

图 2

图 2 真实海洋数据 Fig.2 Field marine data

图 3

图 3 随机采样道数为50%的海洋数据 Fig.3 Randomly sampled marine data with 50% traces sampled

图 4

图 4 基于低冗余曲波变换的海洋数据重建结果 Fig.4 Reconstructed marine data based on low-redundancy curvelet transform

图 5

图 5 一道原始海洋数据与重建数据对比 Fig.5 Comparison of one original trace and reconstructed trace of marine data

为了验证不同采样比例下该低冗余曲波变换的重建效果, 对一个含噪声的三维海洋数据分别采样20%和10%的地震道进行重建, 该试验同时也验证了低冗余曲波变换对含噪声数据重建的可靠性。由于计算机内存的限制, 截取一小部分数据进行计算, 数据为200×200×200, 即在时间上有200个采样点, 空间两个方向上分别有200个采样, 时间采样率为4 ms, 道间距为12.5 m。图 6(a)为原始数据的一个切片显示, 图 6(b)为随机采样20%地震道的显示, 图 6(c)为图 6(b)采用低冗余曲波变换得到的重建结果。由图可见, 不仅未采样的地震道被重建出来, 而且噪声的能量得到了压制。图 6(d)为原始数据和重建结果的误差, 视觉上有效信号很少。为了测试更低采样率的数据能否得到较好的重建, 将采样率降至10%, 即随机选择10%的地震道进行重建, 如图 7(b)所示。图 7(c)为基于低冗余曲波变换的重建结果, 图 7(d)为原始数据和重建数据的差。由此可见, 图 7(c)的结果是可以接受的, 并且图 7(d)中有效信号很少。为了在频率中比较, 图 8(a)为原始数据图 7(a)的一个空间切片的频谱, 图 8(b)为重建数据(图 7(c))相同空间切片的频谱, 图 8(c)为图 7(d)相同空间切片的频谱, 可以看出基于低冗余曲波变换的重建可以滤掉高波数分量和一些高频分量。本试验说明:①低冗余曲波变换对于非常低的采样比例的数据也能够得到稳定的结果, 当采样比例为5%时, 数值效果会变差。②通过选择合适的参数, 采用的模型和算法能够实现含噪声数据的重建, 即问题(6)和(7)的模型可以解决含随机噪声的重建问题。③相比于二维变换, 三维变换需要的采样比例更少。之前的二维数据重建研究表明, 当以曲波变换为稀疏变换, 并且采用阈值类算法如凸集投影法^[26]求解时, 二维数据的采样比例在50%左右才能得到较好的结果, 一旦低于这个比例, 则数值效果会变差, 三维变换由于将两个空间维数的信息整合起来, 具有更稀疏的表达, 因此能够适合更低的采样比例。同样, 五维数据重建也适合更低的采样比例。

图 6

图 6 三维真实数据随机采样20%的重建试验 Fig.6 3D field data reconstruction with 20% traces sampled randomly

图 7

图 7 三维真实数据随机采样10%的重建 Fig.7 3D field data reconstruction with 10% traces sampled randomly

图 8

图 8 真实数据和重建数据的频谱比较 Fig.8 Spectrum comparation of original data and reconstructed data

将一种低冗余度曲波变换用于加速地震数据重建, 给出了一种快速的迭代阈值法实现地震数据重建和去噪, 提高了基于曲波变换的地震数据重建的计算效率和重建质量。数据试验表明，基于低冗余度曲波变换的重建时间约为基于原始曲波变换的1/4, 且能够达到可靠的重建结果, 可以适合更低采样比例的数据重建。基于多尺度变换的重建存在一定的局限性, 由于曲波变换是局部的、多尺度变换, 如果没有对小尺度构造充分的采样, 则基于这种多尺度变换的重建不能得到满意的效果。