地震数据受到施工条件、采集孔径、仪器故障等影响, 难免会出现缺失及不规则的情况, 对后续处理效果产生很大影响, 特别是表面多次波衰减、波动方程偏移等基于多道处理的算法^[1-4]。KABIR等^[5]利用抛物Radon变换^[6]重建近偏移距缺失数据, 取得了良好效果, 但只能处理规则采样数据, 且需要多次迭代。SACCHI等^[7]提出高分辨率Radon变换方法, 并将其成功应用于近偏移距地震道重建, 但在提高分辨率的同时增加了计算量。王维红等^[8]提出并应用加权抛物线Radon变换法进行数据规则化和重建, 在近偏移距地震数据重建中取得了良好效果。目前应用Radon变换进行数据重建的研究主要集中在提高分辨率和减小计算量两个方面, 且多针对各向同性条件下的近偏移距地震数据。海上地震数据多为远偏移距, 且深层陡倾角成像需要准确的远偏移距数据, 准确重建偏移距较远的地震数据同样十分重要^[9-10]。

相比于抛物Radon变换, 双曲Radon变换分辨率更高, 进行数据重建具有更高的精度。但目前使用双曲Radon变换重建远偏移距地震数据面临以下问题:①常规双曲Radon变换公式基于水平或中等倾角地质构造的反射同相轴在CMP道集上呈双曲线走时的假设, 而偏移距较远的反射波同相轴已经不能用简单的双曲线走时公式进行准确描述;②双曲Radon变换具有时变属性, 变换算子只能在时间域表示, 利用共轭梯度法求解计算量巨大, 这也是时间域求解Radon变换的固有问题^[11-13]。本文提出各向异性Radon变换数据重建方法, 将常规双曲变换公式改为由偏移距、慢度、非椭圆率三参数控制的各向异性变换公式^[14];用FISTA算法^[15]取代常规共轭梯度求解方法;采用稀疏约束控制各向异性Radon变换, 通过正则化参数λ平衡反演误差与模型稀疏化程度^[16]。理论模型和实际资料重建结果表明, 本文方法对中远偏移距地震数据有着很好的重建效果。

1 方法原理 1.1 各向异性Radon变换原理

任意水平层状介质模型反射波走时方程可由级数展开得到与层厚、层速度及偏移距2阶、4阶、6阶等有关的复杂函数。当模型为各向同性且偏移距较小时(最大偏移距小于反射界面深度), 可忽略4阶及更高阶项, 用简单的双曲线公式来描述纵波反射同相轴。常规双曲Radon变换定义为:

$m\left( \tau ,\text{ }p \right)={{\int }_{x}}d(t=\text{ }\sqrt{{{\tau }^{2}}+{{x}^{2}}{{p}^{2}}}~,\text{ }x)dx~$

(1)

式中:τ是截距时间;p是慢度(速度倒数);x是炮检距;d(t, x)为道集数据。

但长偏移距纵波反射同相轴时距曲线不满足双曲线规律, 需用非双曲公式进行精确描述^[17]。TSVANKIN^[18]给出了水平单层各向异性介质长偏移距动校时差公式, ALKHALIFAH等^[19]给出了该公式用动校速度和非椭圆率参数表示的形式。权衡计算精度与效率, 本文采用ALKHALIFAH等^[19]给出的非双曲时差公式定义了各向异性Radon变换, 其正变换公式如下。

$\begin{array}{*{35}{l}} m\left( \tau ,\text{ }p \right)={{\int }_{x}}d\text{ }\!\![\!\!\text{ }t= \\ \sqrt{{{\tau }^{2}}+{{x}^{2}}~{{p}^{2}}-\text{ }\frac{2\eta {{x}^{4}}~{{p}^{4}}~}{{{\tau }^{2}}+\left( 1+2\eta \right){{x}^{2}}~{{p}^{2}}}}~,\text{ }x]dx \\ \end{array}$

(2)

与常规Radon变换相比, 各向异性Radon变换在包含高阶项的同时, 增加了非椭圆率参数η, 该参数与每个采样时刻一一对应。要想获得高精度的各向异性Radon变换域数据, 需要准确的η值, 因此如何求取实际数据的η值显得尤为重要。理论模型中的η值是已知的, 对于实际地震数据, 本文采用双参数速度分析的方法来求取准确的η值^[20-21]。虽然各向异性Radon变换公式更加复杂, 但是离散求解时矩阵的维数、大小未变, 计算量与常规双曲Radon变换相当。

为避免数据重建时产生假频^[22], 得到高分辨率的采样结果, 需要合理选择参数, 下面讨论参数的选取准则。参数的采样率为:

$\Delta p={{p}_{k}}-{{p}_{k-1}}$

(3)

已知p=1/v, 则(3) 式可化为:

$\Delta p=\text{ }\frac{1}{{{v}_{k}}}\text{ }~-\text{ }\frac{1}{{{v}_{k-1}}}$

(4)

令Δv=v_k-1-v_k, 得:

$\Delta p=\text{ }\frac{\Delta v}{{{v}_{k}}{{v}_{k-1}}}$

(5)

当v_k和v_k-1比较接近时, v_k≈v_k-1, 公式(5) 变为:

$\Delta p=\text{ }\frac{\Delta v}{{{v}_{k}}^{2}}\text{ }$

(6)

设原始数据中最大有效反射速度为v_max, 要使(6) 式对所有的有效波场均成立, 则式中v_k应替换为v_max。同时, 在CSP或CMP道集上做变换时, Δv选取的范围通常为50～100 m/s, 因而得到临界值关系如下:

$\Delta p=\text{ }\frac{\Delta v}{{{v}_{max}}^{2}}$

(7)

确定了Δp的值后, 根据p=1/v给出p的初始值p₀和最大值p_max, 就可以在不增加数据重建计算量的同时保证重建的准确性^[15]。

1.2 稀疏约束反演

为了得到更准确的重建结果, 本文在常规最小二乘反演的基础上, 提出了基于稀疏约束的反演方法, 采用正则化参数λ来平衡模型的稀疏化与反演误差, 从而控制重建精度。高分辨率Radon变换通常先定义反变换, 再通过求解最小二乘约束下的目标函数求取正变换。Radon变换反变换公式如下。

$d\text{ }=\text{ }Lm$

(8)

式中: d 和 m 分别表示离散原始数据与Radon域数据的矩阵形式; L 和 L ^T分别为Radon变换算子和伴随变换算子。因Radon变换的非正交性, 建立最小平方意义下的线性化反演误差目标函数, 其范数条件的解为:

$m\text{ }={{(\text{ }L{{~}^{T}}~L\text{ })}^{-1}}~L{{~}^{T}}~d$

(9)

在时间域内进行稀疏约束可以很好地提高Radon域数据的分辨率。将Radon域数据作为稀疏条件进行约束反演, 对反演误差取l₂范数, 对稀疏约束模型取l₁范数, 即通过线性反演问题的l₁-l₂混合范数求解, 建立如下的目标函数:

$J=\text{ }\left\| d\text{ }-\text{ }Lm \right\|{{~}_{2}}^{2}+\lambda \left\| \text{ }m \right\|{{~}_{1}}$

(10)

式中:λ是正则化参数, 用于平衡反演误差与模型的稀疏化, 其值越大, Radon域内数据越稀疏。对于规则的模拟数据, λ通常取0, 而对实际数据, 合适的λ值可以使重建结果更加准确。

1.3 各项异性Radon变换求解算法

计算量大是时间域求解Radon变换的固有问题:时间域算子 L 和 L ^T是大型的稀疏矩阵, 每次迭代都需要计算大型稀疏矩阵 L 和 L ^T与向量的乘法^[16]。目前通常采用共轭梯度法求解 d - Lm ², 对超定问题给出最小平方解, 对欠定问题给出最小范数解。求解过程可理解为用共轭梯度法求解N×N阶线性方程组 A · x = b 的问题, 求解时需要构造4个向量序列 r _k, r_k, p _k, p_k, k=1, 2, …。只要提供初始向量 r _k和 r_k, 并设 p ₁= r ₁, p₁= r₁, 即可得到递推公式:

$\alpha {{~}_{k}}=\text{ }\frac{r{{~}_{k}}\cdot ~{{{\bar{r}}}_{k}}~}{p{{~}_{k}}\cdot \text{ }A\text{ }\cdot \bar{p}{{~}_{k}}}$

(11a)

$r{{~}_{k+1}}=\text{ }r{{~}_{k}}-\text{ }\alpha {{~}_{k}}~A\text{ }\cdot \text{ }p{{~}_{k}}$

(11b)

$_{k+1}=~{{{\bar{r}}}_{k}}-\text{ }\alpha {{~}_{k}}~A\text{ }\cdot \bar{p}{{~}_{k}}$

(11c)

$\beta \text{ }=\text{ }\frac{r{{~}_{k+1}}\cdot {{{\bar{r}}}_{k+1}}}{r{{~}_{k}}\cdot \bar{r}{{~}_{k}}}$

(11d)

$p{{~}_{k+1}}=\text{ }r{{~}_{k+1}}+\text{ }\beta {{~}_{k}}~p{{~}_{k}}$

(11e)

${{{\bar{p}}}_{k+1}}=\bar{r}{{~}_{k+1}}+\text{ }\beta {{~}_{k}}~{{{\bar{p}}}_{k}}$

(11f)

可见共轭梯度算法的递推步骤较多, 在实际运算中收敛速度较慢。

为了提高计算效率, 我们采用快速迭代软阈值算法(FISTA)加快反演的收敛速度, 迭代更新公式如下:

$m{{~}_{0}}={{(\text{ }L{{~}^{T}}~L\text{ })}^{z-1}}~L{{~}^{T}}~d$

(12a)

$x{{~}_{0}}=\text{ }m{{~}_{0}}$

(12b)

${{t}_{0}}=1$

(12c)

当k≥0时:

$x{{~}_{k+1}}=soft\text{ }\!\![\!\!\text{ }m{{~}_{k}}+\frac{\text{ }L{{~}^{T}}(\text{ }d\text{ }-\text{ }L\text{ }m{{~}_{k}})}{a}\text{ },\text{ }\frac{\lambda \text{ }}{2a}]\text{ }$

(13a)

${{t}_{k+1}}=\text{ }\frac{1+\text{ }\sqrt{1+4{{t}_{k}}^{2}}}{2}$

(13b)

$m{{~}_{k+1}}=\text{ }x{{~}_{k+1}}+\text{ }\frac{{{t}_{k}}-1}{{{t}_{k+1}}}\text{ }~(\text{ }x{{~}_{k+1}}-\text{ }x{{~}_{k}})$

(13c)

式中:k是当前迭代次数;soft是取软阈值算子;a是Lipschitz常数, a≥max L ^T L )]。这样经过数次迭代就可以获得混合域高分辨率的Radon变换结果。

2 模型试算及实际资料处理 2.1 理论模型

为检验本文各向异性Radon变换数据重建方法的效果, 建立了一个远偏移距数据模型, 模拟地震记录共81道, 道间距100 m, 每道500个采样点, 采样间隔4 ms, 各向异性参数η取固定值0.1, 如图 1a所示。图 1b为该模型不规则采样及缺失远偏移距数据的情况, 即对原始数据前50道间隔抽样, 且使第66至75道(共10道)远偏移距数据缺失。采用高分辨率双曲Radon变换和各向异性Radon变换分别进行数据重建, 用FISTA算法迭代100次, 其Radon域结果分别如图 1c和图 1d所示, 重建后的结果分别如图 1e和图 1f所示。可以看出, 各向异性变换使Radon域内聚焦效果更好, 重建结果更接近原始数据, 特别是远偏移距缺失数据的重建结果。

为进一步比较两种方法的重建误差, 图 2a和图 2b分别给出了高分辨率双曲Radon变换和各向异性Radon变换重建数据与原始数据之差。图 2c, 图 2d, 图 2e, 图 2f分别给出了原始数据、不规则缺失数据、高分辨率双曲Radon变换重建数据和各向异性Radon变换重建数据的f-k谱, 对比可见, 不规则缺失数据出现严重假频, 两种Radon变换重建方法都能明显消除这些假频。为准确描述重建精度, 定义重建误差如下:

${E_\Delta } = \left| {\frac{{{\rm{ }}Lm{\rm{ }} - {\rm{ }}d}}{d}} \right| \times 100\% $

(14)

高分辨率双曲Radon变换重建误差E_Δ=9.84%, 各向异性Radon变换重建误差E_Δ=3.71%, 各向异性Radon变换重建误差明显更小。

在求解各项异性Radon变换时, 我们分别采用共轭梯度算法和FISTA算法进行了对比测试, 运行环境为Intel(R)Core(TM)i5-2300 CPU@2.80 GHz, 8 GB DDR3 1 333 MHz内存。FISTA算法在迭代100次的情况下重建误差为3.71%, 耗时159.8 s;共轭梯度法至少迭代150次才能使误差减小到4%以下, 耗时为377.5 s。可见FISTA算法不但计算效率高, 而且收敛速度快。

2.2 实际数据处理分析

采用我国某地区实际地震资料CMP道集, 验证本文各向异性Radon变换数据重建方法的实用性。该数据共126道, 道间距24 m, 每道626个采样点, 采样间隔4 ms, 记录长度2.5 s, 其中有效反射波的速度范围在1 200~1 600 m/s。进行各向异性Radon变换数据重建之前, 首先对道集进行双参数速度扫描, 得到每个地震同相轴所对应的η参数, 即扫描得到的η参数与同相轴个数相同;然后将η参数按时间方向插值, 得到每个离散时间采样点对应的η值, 如图 3所示。

图 4a是原始实际数据, 图 4b是对实际数据前50道间隔抽样并缺失第101到110道远偏移距数据的结果, 图 4c为其各向异性Radon变换的结果, 图 4d为各向异性Radon变换方法数据重建结果。可以看出实际数据重建效果较好, 说明本文方法具有较强的实用性。分别对原始实际数据及其不规则缺失地震道数据和重建后数据做f-k谱分析, 结果如图 5a, 图 5b, 图 5c所示, 可见不规则缺失地震道数据出现了明显假频, 利用本文方法重建数据后有效地消除了这些假频。

为进一步提高各向异性Radon变换重建效果, 在实际数据处理中加入一定程度的稀疏约束, 并取实际数据与重建数据的第101至110道进行对比。λ的取值范围一般为0.000 1~0.010 0, 在该区间用二分法求解λ, 仅需求解几次即可得到最佳λ值0.002。图 6a为实际数据第101至110道, 图 6b为常规Radon变换的重建结果, 图 6c为λ=0时的各向异性Radon变换重建结果, 图 6d为λ=0.002时的各向异性Radon变换重建结果。图 6b, 图 6c, 图 6d重建误差分别为11.84%, 6.52%, 5.04%, 说明各向异性Radon变换重建效果更好, 且λ=0.002时重建结果更接近实际数据。

3 结束语

本文提出的各向异性Radon变换重建方法可以更好地恢复远偏移距地震数据, 弥补了常规Radon变换方法在远偏移距数据重建方面的不足。对于不同地区的实际数据, 各向异性Radon变换前需进行双参数速度扫描, 获得相应正确的各向异性参数η, 确保Radon域内数据的高分辨率。采用稀疏约束反演的方法可以减小数据重建误差, 得到更为准确的重建结果。在时间域Radon变换的求解过程中, 采用FISTA算法取代常规的共轭梯度算法可以加快收敛速度, 提高运算效率。