Radon变换是地震数据处理中常用的数学算法之一, 广泛用于地震数据速度分析、多次波压制、波场分离和插值等领域。该变换可分为线性、抛物线和双曲线变换3类, 其中线性Radon变换也称τ-p变换。Radon变换由CLAERBOUT^[1]为首的斯坦福地球物理团队引入到地球物理领域, 并初步使用双曲线Radon变换进行速度分析和反演; 当时的算法是在时间域进行的, 需要求取大型矩阵的逆, 计算量大; 而HAMPSON^[2]和BEYLKIN等^[3]重新定义了抛物线Radon变换, 给出了最小二乘解, 并成功应用于地震数据多次波压制中。抛物线Radon变换在应用时, 假定动校正后的共中心点道集中, 一次波是被校正平, 由于校正所用的速度偏大, 多次波同相轴校正不足, 存在一定剩余时差, 即多次波同相轴是向时间的正方向弯曲。数学推导表明, 动校正后多次波时距曲线的形状可以近似用数学上的抛物曲线公式表示。因此, 可以通过抛物线Radon变换, 将CMP道集中拉平的一次波和弯曲的多次波变换到Radon域。Radon域中的横坐标为剩余时差大小, 一次波和多次波剩余时差存在的差异, 使得二者在Radon域能够更加有效地分离。为更好地保持有效波的能量, 多次波压制的一般做法是在Radon域切除一次波, 并反变换回时空域得到多次波模型, 最后将该多次波模型从原数据中减去。

与地震数据处理中常用的其它变换, 如傅里叶变换和小波变换等不同, Radon变换算子并不是正交的, 是不可逆的。对同一个数据做正、反变换后, 反变换后的数据会有能量的损失。最常用的解决手段是反演, 即在一定限制条件下, 将原始数据和变换回来的数据之间的差异最小化。最常用的一种最小化策略就是HAMPSON^[2]和BEYLKIN等^[3]给出的最小二乘范数解, 即让原始数据和变换回来的数据在最小二乘意义下差最小。这种算法在线性和抛物线Radon变换中最为常用。

除了非正交问题外, Radon变换还受采集空间和时间有限、离散采样等因素的影响, 存在假频和分辨率低等问题。变换域能量无法有效聚焦, 使得多次波和一次波的分离困难, 导致在压制结果中存在多次波的残余, 或者损伤一次波。因此, 提高Radon变换的分辨率一直是该方法研究的热点。提高Radon域分辨率的算法, 大致可以分为3类。第1类是在时间域提高分辨率, 该类算法可以很好地衰减变换域的假象, 得到清晰的高分辨率结果, 但由于需要引入时间域的大型稀疏算子, 计算效率低, 因此很少采用, 在此不作详细介绍。第2类是在频率域迭代提高分辨率的算法, 该算法将Radon变换视为一个稀疏反演问题, 将上一次计算得到的模型结果, 作为下一次迭代计算的加权算子, 以达到提高变换域分辨率的目的。但频率域提高分辨率的算法, 对同一个地震道的所有时间施加了相同的约束, 因此加权项对所有时间是耦合在一起的, 会使得能量较大的同相轴出现假象^[4]。频率域的算法有: SACCHI等^[5]提出的基于柯西范数概率密度函数Bayes理论的高分辨率算法、SACCHI等^[6]提出的基于迭代共轭梯度逐次更新阻尼参数的高分辨率算法、HERRMANN等^[7]提出的基于低频约束的去假频高分辨率Radon变换算法和CHEN等^[8]提出的能够保护地震数据弱信号、去除假频的基于地震数据特定频率解约束的非迭代高分辨率算法等。第3类提高Radon变换分辨率的方法是频率时间域的混合算法, 该算法一般是对最小二乘算法得到的时间域模型进行收敛处理, 而其中的正向和反向Radon变换均在频域进行。该算法综合了频率和时间域两种算法的优势, 既提高了分辨率, 又保持了数据的波形, 同时还没有明显增加计算时间。如LU^[9]提出的基于加速稀疏时不变的频率-时间混合域迭代收缩高分辨率Radon变换算法, 在没有明显增加计算时间的同时, 有效提高了分辨率; 马继涛等^[10]对二维情况下不同迭代阈值收缩高分辨率Radon变换算法进行了对比分析; 薛亚茹等^[11]对加权迭代软阈值算法的高分辨率Radon变换算法进行了探讨; 罗腾腾等^[12]对混合域高分辨率Radon变换在绕射波分离和成像中的应用进行了探讨。

三维地震数据, 尤其宽方位数据, 形成的三维小面元道集, 受到地层各向异性等因素的影响, 同相轴存在抖动现象, 二维Radon变换无法解决该道集的多次波压制问题。需要考虑三维地震波场的传播特征, 利用三维Radon理论对不同方位的数据进行精确表征, 能更好地对其进行处理。三维Radon变换考虑了波场三维传播的特点, 利用纵测线和横测线两个方向的曲率参数, 在Radon域对动校正后的三维小面元道集进行表征, 可以实现三维数据多次波的有效压制。很多学者针对三维地震波场传播的特点, 对三维Radon算法展开过较为深入地研究。如DONATI等^[13]利用三维τ-p变换对地震数据进行插值重建处理; HUGONNET等^[14-15]给出了三维抛物线Radon变换算法并应用在宽方位地震数据处理中; 如ZHANG等^[16]给出了一种时间和频率混合域加速三维稀疏时不变τ-p变换算法, 即基于迭代阈值收缩的高分辨率τ-p变换算法, 并应用在三维叠前道集的插值重建中; CAO等^[17]给出了一种基于匹配追踪的高分辨率三维τ-p变换算法, 有效解决了空间假频问题, 提高了分辨率, 取得了较好的插值效果; SUN等^[18]给出了一种用于线性噪声压制的三维圆锥线性Radon变换算法, 用于压制三维道集中的面波。针对三维Radon变换过程中的保幅问题, 唐欢欢等^[19-20]给出了一种高阶高分辨率抛物线保幅Radon变换算法并应用于三维地震数据重建; 薛亚茹等^[21]给出了一种高阶3D变换地震数据重建算法, 有效解决了地震数据重建过程中的保幅问题; MA等^[22]给出了一种高阶高分辨率三维Radon变换算法, 考虑了三维Radon变换的保幅性, 并基于低频约束的方法提高地震数据的分辨率。

针对常规三维Radon变换算法分辨率低的问题, 很多学者提出了提高分辨率的算法, 包括时间域、频率域和混合域的算法, 但没有对这些算法进行系统的对比分析。我们基于Radon变换多次波压制, 对常用的频率域和时间频率混合域提高分辨率的算法进行了系统的对比分析。首先, 详细介绍了基于最小二乘的三维Radon变换算法; 然后, 分别给出了基于迭代重加权的高分辨率算法, 利用低频解约束高频运算的高分辨率算法和基于迭代阈值收缩的高分辨率算法; 最后, 利用模拟数据对每一种算法多次波压制的有效性进行验证和对比分析, 并分析各个算法的抗噪性能。

1 三维抛物线Radon变换原理 1.1 基于最小二乘的三维抛物线Radon变换

三维抛物线Radon变换沿某个抛物面对地震数据求和, 得到三维Radon域数据。设时空域三维地震数据为d(t, x, y), 变换后三维Radon域数据为m(τ, q_x, q_y), 其中, x和y分别为纵测线和联络测线两个方向的偏移距; q_x和q_y分别为沿着纵测线和联络测线方向的两个曲率参数; t和τ分别为时空域和Radon域的时间。根据Radon变换的定义, 沿抛物面路径对Radon域数据进行叠加求和, 可得到时空域地震数据, 则有:

$ \begin{gathered} \boldsymbol{d}(t, x, y)=\iint \boldsymbol{m}\left(\tau=t-q_{x} x^{2}-q_{y} y^{2}, q_{x},\right. \\ \left.q_{y}\right) \mathrm{d} q_{x} \mathrm{~d} q_{y} \end{gathered} $

(1)

公式(1)可离散为:

$ \begin{gathered} \boldsymbol{d}(t, x, y)=\sum\limits_{j=1}^{n_{q x}} \sum\limits_{k=1}^{n_{q y}} \boldsymbol{m}\left(\tau=t-q_{x j} x^{2}-\right. \\ \left.q_{y k} y^{2}, q_{x j}, q_{y k}\right) \end{gathered} $

(2)

式中: n_qx和n_qy为曲率qx和qy的个数。对公式(2)两侧关于时间做傅里叶变换, 将其转变到频率域, 可得:

$ \boldsymbol{D}(\omega, x, y)=\sum\limits_{j=1}^{n_{q x}} \sum\limits_{k=1}^{n_{q y}} \mathrm{e}^{-\mathrm{i} \omega q_{x j} x^{2}} \boldsymbol{M}\left(\omega, q_{x j}, q_{y k}\right) \mathrm{e}^{-\mathrm{i} \omega q y k y^{2}} $

(3)

式中: D为频率域的地震数据; M为频率域的Radon数据。对某个特定的x_xj和y_yk, 可写出D的矩阵计算式, 如下:

$ \begin{aligned} &\boldsymbol{D}\left(\omega, x_{x j}, y_{k y}\right)=\left[\begin{array}{llll} \mathrm{e}^{-\mathrm{i} \omega q_{x_1} x_{j x}^{2}} & \mathrm{e}^{-\mathrm{i} \omega q_{x_{2}} x_{j x}^{2}} & \cdots & \mathrm{e}^{-\mathrm{i} \omega q_{n_{q x}} x_{j x}^{2}} \end{array}\right] \cdot\\ &\ \ \ \ \ \ \ \ \left[\begin{array}{cccc} \boldsymbol{M}\left(\omega, q_{x_{1}}, q_{y_{1}}\right) & \boldsymbol{M}\left(\omega, q_{x_{1}}, q_{y_{2}}\right) & \cdots & \boldsymbol{M}\left(\omega, q_{x_{1}}, q_{n_{q y}}\right) \\ \boldsymbol{M}\left(\omega, q_{x_{2}}, q_{y_{1}}\right) & \boldsymbol{M}\left(\omega, q_{x_{2}}, q_{y_{2}}\right) & \cdots & \boldsymbol{M}\left(\omega, q_{x_{2}}, q_{n_{q y}}\right) \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{M}\left(\omega, q_{n_{q x}}, q_{y_{1}}\right) & \boldsymbol{M}\left(\omega, q_{n_{q x}}, q_{y_{2}}\right) & \cdots & \boldsymbol{M}\left(\omega, q_{n_{q x}}, q_{n_{q y}}\right) \end{array}\right]\left[\begin{array}{c} \mathrm{e}^{-\mathrm{i} \omega q_{y_{1}} y_{k y}^{2}} \\ \mathrm{e}^{-\mathrm{i} \omega q_{y_{2}} y_{k y}^{2}} \\ \vdots \\ \mathrm{e}^{-\mathrm{i} \omega q_{n_{q y}} y_{k y}^{2}} \end{array}\right] \end{aligned} $

(4)

公式(4)可以扩充至所有的数据D, 即:

$ \begin{array}{c} \left[\begin{array}{cccc} \boldsymbol{D}\left(\omega, x_{1}, y_{1}\right) & \boldsymbol{D}\left(\omega, x_{1}, y_{2}\right) & \cdots & \boldsymbol{D}\left(\omega, x_{1}, y_{n y}\right) \\ \boldsymbol{D}\left(\omega, x_{2}, y_{1}\right) & \boldsymbol{D}\left(\omega, x_{2}, y_{2}\right) & \cdots & \boldsymbol{D}\left(\omega, x_{2}, y_{n y}\right) \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{D}\left(\omega, x_{n x}, y_{1}\right) & \boldsymbol{D}\left(\omega, x_{n x}, y_{2}\right) & \cdots & \boldsymbol{D}\left(\omega, x_{n x}, y_{n y}\right) \end{array}\right]=\left[\begin{array}{cccc} \mathrm{e}^{-\mathrm{i} \omega q_{x_{1}} x_{1}^{2}} & \mathrm{e}^{-\mathrm{i} \omega q_{x_{2}} x_{1}^{2}} & \cdots & \mathrm{e}^{-\mathrm{i} \omega q_{n _{qx}} x_{1}^{2}} \\ \mathrm{e}^{-\mathrm{i} \omega q_{x_{1}} x_{2}^{2}} & \mathrm{e}^{-\mathrm{i} \omega q_{x_{2}} x_{2}^{2}} & \cdots & \mathrm{e}^{-\mathrm{i} \omega q_{n _{q x}} x_{2}^{2}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathrm{e}^{-\mathrm{i} \omega q_{x_{1}} x_{n x}^{2}} & \mathrm{e}^{-\mathrm{i} \omega q_{x_{2}} x_{n x}^{2}} & \cdots & \mathrm{e}^{-\mathrm{i} \omega q_{n_{q x}} x_{n x}^{2}} \end{array}\right]\cdot\\ \left[\begin{array}{cccc} \boldsymbol{M}\left(\omega, q_{x_{1}}, q_{y_{1}}\right) & \boldsymbol{M}\left(\omega, q_{x_{1}}, q_{y_{2}}\right) & \cdots & \boldsymbol{M}\left(\omega, q_{x_{1}}, q_{n_{q y}}\right) \\ \boldsymbol{M}\left(\omega, q_{x_{2}}, q_{y_{1}}\right) & \boldsymbol{M}\left(\omega, q_{x_{2}}, q_{y_{2}}\right) & \cdots & \boldsymbol{M}\left(\omega, q_{x_{2}}, q_{n_{q y}}\right) \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{M}\left(\omega, q_{n_{q x}}, q_{y_{1}}\right) & \boldsymbol{M}\left(\omega, q_{n_{q x}}, q_{y_{2}}\right) & \cdots & \boldsymbol{M}\left(\omega, q_{n_{q x}}, q_{n_{q y}}\right) \end{array}\right]\left[\begin{array}{cccc} \mathrm{e}^{-\mathrm{i} \omega q_{y_{1}} y_{1}^{2}} & \mathrm{e}^{-\mathrm{i} \omega q_{y_{1}} y_{2}^{2}} & \cdots & \mathrm{e}^{-\mathrm{i} \omega q 1 y_{n y}^{2}} \\ \mathrm{e}^{-\mathrm{i} \omega q_{y_{2}} y_{1}^{2}} & \mathrm{e}^{-\mathrm{i} \omega q_{y_{2}} y_{2}^{2}} & \cdots & \mathrm{e}^{-\mathrm{i} \omega q 2 y_{n y}^{2}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathrm{e}^{-\mathrm{i} \omega q_{n_{q y}} y_{1}^{2}} & \mathrm{e}^{-\mathrm{i} \omega q_{n_{q y}} y_{2}^{2}} & \cdots & \mathrm{e}^{-\mathrm{i} \omega q_{n_{q y}} y_{n y}^{2}} \end{array}\right] \end{array} $

(5)

令M两侧与q_xx²相关的矩阵为L_x, 与q_yy²相关的矩阵为L_y, 则式(5)可以写为如下的算子形式:

$ \boldsymbol{D}=\boldsymbol{L}_{x} \boldsymbol{M} \boldsymbol{L}_{y} $

(6)

式中: L_x为与q_x、偏移距x有关的算子, 大小为(n_x, n_qx); L_y为与q_y、偏移距y有关的算子, 大小为(n_qy, n_y)。在n_x≠n_qx, n_y≠n_qy情况下, L_x与L_y均不是方阵, 可采用与二维Radon变换类似的方式求取L_x和L_y矩阵的最小二乘逆。在实际数据处理时, 为保证处理精度, 总是使变换所用曲线的数目大于地震道的道数, 即n_qx>n_x, n_qy>n_y, 这就使得矩阵L_x为欠定矩阵, L_y为超定矩阵。同时, 考虑到公式(6)各矩阵的大小, 对L_x和L_y采用不同的方式求解最小二乘逆, 如下式:

$ \widetilde{\boldsymbol{M}}=\boldsymbol{L}_{x}^{\mathrm{H}}\left(\boldsymbol{L}_{x} \boldsymbol{L}_{x}^{\mathrm{H}}+\mu \boldsymbol{I}\right)^{-1} \boldsymbol{D}\left(\boldsymbol{L}_{y}^{\mathrm{H}} \boldsymbol{L}_{y}+\mu \boldsymbol{I}\right)^{-1} \boldsymbol{L}_{y}^{\mathrm{H}} $

(7)

式中: I是单位对角矩阵, μ是为提高矩阵求逆运算稳定性施加的阻尼因子, 取值一般在0.01~1.00之间; H代表矩阵的共轭转置。

1.2 基于迭代重加权的三维高分辨率Radon变换

公式(7)中, 阻尼因子μ对所有频率都是相同的, 使得矩阵求逆更加稳定, 但同时也降低了Radon变换结果的分辨率。为提高分辨率, SACCHI等^[6]提出了基于迭代重加权的高分辨率Radon变换算法, 该算法利用上一次迭代计算的模型域结果, 在信号区域给μ赋一个较小的值, 在非信号区赋一个较大的值, 从而将Radon域能量强的点加强, 能量弱的点减弱, 使Radon域能量聚焦, 达到提高Radon域数据分辨率的效果, 该变换的表达式如下。

$ \widetilde{\boldsymbol{M}}=\boldsymbol{L}_{x}^{\mathrm{H}}\left(\boldsymbol{L}_{x} \boldsymbol{L}_{x}^{\mathrm{H}}+\mu \boldsymbol{C}_{m x}\right)^{-1} \boldsymbol{D}\left(\boldsymbol{L}_{y}^{\mathrm{H}} \boldsymbol{L}_{y}+\mu \boldsymbol{C}_{m y}\right)^{-1} \boldsymbol{L}_{y}^{\mathrm{H}} $

(8)

式中: C_mx和C_my是和上一次迭代运算结果相关的加权矩阵。每次迭代都会生成新的加权矩阵C_mx和C_my, 以得到高分辨率的三维Radon域结果。以C_mx为例, 加权矩阵可通过式(9)计算得到:

$ \boldsymbol{C}_{m x}=\operatorname{diag}\left(\frac{1}{\boldsymbol{M}^{2}+\varepsilon^{2}}\right) $

(9)

式中: ε为使得求倒数稳定所加的稳定因子, diag(·)是将向量变为对角矩阵的算子。如果变换域M的值较小, 则计算出的约束项使得下一次迭代求解出的M的值更小, 亦即算法会使得变换中的假象噪声变小; 反之算法会使得变换中的信号增强, 达到提高分辨率的目的。

1.3 基于低频约束的三维高分辨率Radon变换

基于迭代重加权的高分辨率算法在处理稀疏采样数据时, 由于间距大, 空间采样会产生空间假频, 导致变换域出现假象, 降低分辨率。为解决这一问题, HERRMANN^[7]提出了一种不需要迭代的高分辨率Radon变换算法, 利用数据低频部分的解对数据高频部分的运算进行约束, 即用多次的重加权约束Radon域数据, 使其变得稀疏, 每一个加权矩阵都是由上一个频率的计算结果生成。地震数据的低频部分波长长, 变换时不会出现空间假频, 因此低频约束的方法可以阻止假频的出现; 此外, 方法无迭代处理过程, 可以提高计算效率; 方法计算过程中, 除阻尼参数外, 没有其他可调整的参数, 运算相对简单, 可一定程度上减轻处理人员对方法的测试过程。

基于低频约束的三维高分辨率Radon变换可表示为:

$ \begin{gathered} \widetilde{\boldsymbol{M}}_{n}=\boldsymbol{W}_{n} \boldsymbol{L}_{x}^{\mathrm{H}}\left(\boldsymbol{L}_{x} \boldsymbol{W}_{n} \boldsymbol{L}_{x}^{\mathrm{H}}+\mu \boldsymbol{I}\right)^{-1} \cdot \\ \boldsymbol{D}_{n}\left(\boldsymbol{L}_{y}^{\mathrm{H}} V_{n} \boldsymbol{L}_{y}+\mu \boldsymbol{I}\right)^{-1} \boldsymbol{L}_{y}^{\mathrm{H}} V_{n} \end{gathered} $

(10)

式中: M_n和D_n分别为第n个频率的三维Radon域和地震数据; W_n和V_n分别是纵测线和联络测线方向算子求逆所对应的对角加权矩阵。该加权矩阵是由上一次(n－1)频率的结果计算得到。以L_x算子的加权矩阵W_n为例, 可由下式计算得到:

$ \boldsymbol{W}_{i i}\left(\omega_{n}\right)=\left\|\boldsymbol{M}_{i}\left(\omega_{n-1}\right)\right\| $

(11)

式中: M_i(ω_n-1)为在上一次计算频率运算得到的结果, 此处W对角加权矩阵i的范围为1~n_qx; 若为联络测线方向的对角加权矩阵V, 则i的范围为1~n_qy。在没有多次波模型先验信息的情况下, 该方法利用数据的低频计算结果对高频计算进行约束, 低频数据在运算中不产生假频, 在采样稀疏的情况下仍能取得很好的高分辨率效果。

1.4 基于迭代阈值收缩的三维高分辨率Radon变换

迭代阈值收缩算法是在时间-频率混合域利用反变换后数据和原数据的残差更新模型, 进行提高分辨率的处理。每次迭代在时间域进行阈值收缩处理, 增强变换域能量团中心的能量, 减弱边缘的发散能量, 可以使得模型域能量团更加集中, 提高分辨率。

与二维算法类似, 三维算法是提前计算各频率L_x和L_y矩阵的最小二乘逆, 并存储起来在计算中直接调用, 以牺牲小的存储内存为代价提高计算效率; 利用Radon域模型反变换得到数据和原始数据的残差, 通过逐次迭代收缩阈值更新Radon域数据, 达到提高分辨率的效果。我们将LU^[9]给出的二维算法扩展到三维, 设算子L_x和L_y的最小二乘逆分别为L_{x_inv}和L_{y_inv}, 即令:

$ \left\{\begin{array}{l} \boldsymbol{L}_{x_{-} i n v}=\boldsymbol{L}_{x}^{\mathrm{H}}\left(\boldsymbol{L}_{x} \boldsymbol{L}_{x}^{H}+\mu \boldsymbol{I}\right)^{-1} \\ \boldsymbol{L}_{y_{-} i n v}=\left(\boldsymbol{L}_{y}^{\mathrm{H}} \boldsymbol{L}_{y}+\mu \boldsymbol{I}\right)^{-1} \boldsymbol{L}_{y}^{\mathrm{H}} \end{array}\right. $

(12)

则可通过公式(13)计算更新变换域模型:

$ \begin{gathered} m_{k}=T_{\alpha}\left\{m_{k-1}+\beta \boldsymbol{F}^{-1}\left[{\boldsymbol{ L} } _ { x _ { - } i n v } \cdot \left(\boldsymbol{F}[\boldsymbol{d}]-\right.\right.\right. \\ \left.\left.\left.\boldsymbol{L}_{x} \boldsymbol{F}\left[m_{k-1}\right] \boldsymbol{L}_{y}\right) \cdot \boldsymbol{L}_{y_{-} i n v}\right]\right\} \end{gathered} $

(13)

式中: L_{x_inv}和L_{y_inv}在每次迭代中直接调用; m_k是k次迭代后时间域的Radon模型数据; T_α代表对括号内数据的阈值处理; β为步长; d为时空域的地震数据; F和F^-1是傅里叶正变换和反变换算子。由公式(13)可以看出, 迭代阈值收缩算法是将地震数据和Radon域模型反变换后得到的地震数据之差, 以一定比例叠加到Radon域, 更新时空域的模型。之后基于公式(14)的收缩函数, 对Radon域模型数据进行收缩和提高分辨率的运算。

$ T_{\alpha}\{\boldsymbol{m}, k\}_{i j}=\left(\left|m_{i j}\right|-\alpha \frac{K-k}{K} \widetilde{m}_{i j}\right)+\operatorname{sgn}\left(m_{i j}\right) $

(14)

式中: α为介于0和1之间的系数, K为总迭代次数, sgn为取符号的算子, |m_ij|为对模型数据取绝对值后的结果, $ {{\tilde m}_{ij}}$是对Radon域模型数据进行均值滤波后的结果。

2 模型数据试算对比分析

利用模拟数据对方法的有效性和精确性进行验证, 分别展示三维常规最小二乘Radon变换、迭代重加权的高分辨率Radon变换、基于低频约束的高分辨率Radon变换以及基于迭代阈值收缩的高分辨率Radon变换结果。为方便陈述, 以LS(least-squares)代表最小二乘变换算法; 以IRHR(iterated reweight high-resolution)代表迭代重加权高分辨率算法; 以LOWHR(lower frequencies constrained high resolution)代表低频约束高分辨率变换算法; 以SRTIS(sparse radon transform iterative shrinkage)代表迭代阈值收缩高分辨率算法, 对变换域和多次波压制结果进行对比分析。

首先, 模拟生成三维地震数据CMP道集。假定纵测线偏移距范围为-4800~4800m, 每个道集数据有50道; 联络测线方向偏移距范围为-2400~2400m, 测线间距为200m。数据中有4个一次波和4个具有不同剩余时差的多次波, 一次波的剩余时差为0, 多次波的剩余时差分别为600ms, 400ms, 200ms和40ms, 最下方的多次波与一次波之间的时差非常小, 与一次波所在点横向只有0.04s, 即2.3个样点的距离。多次波能量有强有弱, 其振幅由上到下依次为-1.0, 0.5, 0.1, -0.9。三维小面元CMP道集如图 1所示。利用本文算法对此三维CMP道集进行多次波压制测试。

为更好地分析对比多次波压制的效果, 还抽取了原始数据和一次波的第1, 7, 13, 19和25个小面元道集(图 2)。可以看出, 不同位置的小道集中, 多次波和一次波具有不同的时差差异特点, 因而利用相同参数对所有小道集进行处理不可靠。图 3展示的是与图 1对应的按标量偏移距大小排列的三维CMP道集, 在该道集中可以看出, 一次波和多次波受方位各向异性的影响, 同相轴有抖动。基于二维Radon变换无法直接对其进行处理, 需要考虑基于三维抛物线Radon理论的多次波压制方法。

2.1 多次波压制的有效性

分别利用4种算法对该模拟数据进行了处理。变换域应为三维模型(τ, q_x, q_y), 但由于三维模型并不能直观显示各方法的聚焦效果, 将变换域中每个q_y面板的数据进行了累加求和, 并基于该结果对各方法分辨率进行了对比。

图 4为变换域的结果, 每个箭头处都有一次波或多次波聚焦的能量。可以看出, 图 4a的LS算法可以将数据变换到Radon域, 振幅为0.1的弱多次波在变换域也清晰可见, 但该算法假象严重; 一般而言, 能量越强, 同相轴假象越严重。因此在0.1s和1.0s处的多次波, 都有较严重的剪刀状发散假象。切除时会损伤一次波能量, 或导致多次波压制不彻底。而其它3种高分辨率算法都较好地减弱或压制了剪刀状发散假象, 变换后的分辨率都得到了很大提升。在弱多次波的识别方面, LOWHR和SRTIS两种算法都能看到较为清晰的聚焦点, 但IRHR算法在多次迭代后, 将弱多次波误认为假象噪声, 导致该多次波的能量减弱(图 4b), 这会给变换域多次波的识别带来一定干扰。

在图 4中的三维Radon变换结果中, 设置切除参数为0.02, 即将变换域q_x小于0.02的能量切除掉。数据下方1.0s处, 一次波位于横向第2个样点, 多次波位于横向第5个样点, 切除点位于第3个样点, 切除位置位于一次波和多次波中间位置, 且与二者相距非常近, 因此切除对方法分辨率的要求非常高。切除一次波能量后, 保留多次波并将其变换回时空域, 得到图 5所示的多次波模型, 之后将得到的多次波模型从原数据中减掉, 以最大程度保留一次波的能量, 得到图 6和图 7所示的多次波压制后的小道集和标量偏移距道集, 其中图 6为与图 2所对应的小道集压制后的结果, 图 7为与图 3相对应的按标量偏移距大小排列的多次波压制后的结果。在图 5所示的多次波模型中, LS算法有很严重的一次波能量泄漏问题, 而IRHR结果由于分辨率得到了提高, 一次波能量有所减弱; 图 5a和图 5b箭头所指处均为所泄漏的一次波能量。在LOWHR和SRTIS算法结果中未见一次波能量同相轴, 说明这两种高分辨率算法将一次波和多次波分离得更为彻底。在图 6和图 7的多次波压制结果中也可以看出, LS和IRHR算法压制结果中有多次波能量的残余(箭头所指), 同时两种算法由于估计多次波模型中有一次波能量的泄漏, 得到的一次波能量也有损伤, 即LS和IRHR算法既存在多次波压制不彻底的问题, 也存在一次波能量的损伤问题。SRTIS算法多次波压制的结果中(图 6d), 箭头所指处存在多次波能量的残余, 说明该方法虽然在变换域能够很好的分离开一次波和多次波, 但变换回的多次波与数据中的多次波振幅不一致, 导致压制结果有残余, 图 7d箭头所指也可以看到多次波的残余; 与其它方法相比, LOWHR算法无论在多次波和一次波的分离, 还是多次波压制中均取得了较好的多次波压制效果。

表 1列出了采用几种方法处理模拟数据所用的时间。从表 1中可以看出, 各方法运行时间由大到小顺序依次为: LOWHR>IRHR>SRTIS>LS。LS算法耗时最短, 但其效果差; LOWHR算法耗时最长, 但其效果最好。相对而言, SRITS和LOWHR算法均可以取得较好的压制效果, 但SRTIS算法运行时间更短; SRTIS算法为取得较好的结果, 需要调试多个参数, 而LOWHR算法与LS算法类似, 只需调整阻尼参数, 计算较为简单直接。

2.2 多次波压制方法的抗噪性

对模拟数据加入一定量的高斯噪声, 测试了各算法的抗噪性。结果表明, 噪声对各算法的变换分辨率及多次波压制结果均存在一定的影响, 但调整阻尼参数后均可取得较好的结果。本文给出了信噪比为20的情况下, 变换域及多次波压制的效果。图 8为LS、IRHR、LOWHR和SRTIS算法抗噪性的测试结果。其中, LS、LOWHR和SRTIS算法均调大了阻尼参数, IRHR算法保持原阻尼参数不变。可以看出, 各种高分辨率算法在没有噪声的情况下, 均有很好的多次波压制效果, 但弱能量多次波在SRTIS算法变换域更加清晰。说明SRTIS算法充分利用了频率域和时间域两种算法的优势, 既可以得到高分辨率的变换结果, 又保持了数据中弱同相轴的能量, 有利于数据中多次波的识别和压制。

3 结论与建议

系统介绍了三维Radon变换的理论, 给出了求解该问题的最小二乘算法, 基于迭代重加权、低频约束和迭代阈值收缩的高分辨率算法, 并利用模拟数据验证了各方法的有效性、抗噪性和对小时差多次波的压制能力。通过对比和分析研究, 可以得出如下结论。

1) 受采集空间和时间有限、离散采样的影响, 三维Radon变换的最小二乘解存在假象, 分辨率低。

2) 频率域高分辨率算法在变换域时间方向存在噪声干扰, 但可提高横向分辨率; 时间频率混合域高分辨率算法可以在双方向提高变换分辨率。

3) 在有噪声干扰的情况下, 通过调整变换参数, 各算法均可以取得较好的多次波压制效果。

4) 在小时差多次波的压制方面, 迭代收缩阈值高分辨率算法在变换域具有较好的收敛性, 可以有效识别多次波, 但分辨率的提升会影响所估计多次波模型的保幅性; 相对而言, 低频约束高分辨率算法在小时差多次波压制方面效果最佳。

随着反演精度的提高, 对地震数据叠前道集的保幅性要求越来越高, 如何在多次波压制过程中保持有效波的振幅, 需要进一步的研究。另外, 与二维变换算法相比, 三维变换算法计算效率较低, 尽管迭代阈值收缩算法在计算过程中以保存矩阵逆的方式减少了运算时间, 然而三维算法的计算量仍然较大, 不利于大规模生产应用。因此, 需要进一步研究如何提高三维算法的计算效率。