地震数据采集过程中, 受野外采集地质环境、仪器故障、人为干扰、噪声以及经济条件等因素的影响, 采集系统很难记录完整的地震波场信息, 采集到的数据往往会出现缺失道和野值道。这些不规则的或者稀疏的地震数据不仅会产生假频, 也为偏移成像和多次波衰减带来了影响, 从而不能提供精确、规则的地震信号, 给数据的处理与解释工作造成了一定的困扰。因此, 地震数据重建成为地震数据处理中的关键步骤之一。

在实际地震勘探中, 很多数据重建方法已经取得了良好的效果, 大致可以将其分为3类:①基于信号处理的方法^[1], 将地震数据从时间域转换到变换域, 根据变换域内信号的特征进行重建, 如:傅里叶域、小波域、Curvelet域等, 这类方法的处理效果在很大程度上依赖于数学变换的好坏; ②基于波动方程的方法^[2], 根据地震波场在地下介质中传播的物理特征完成重建, 这类方法需要地下的速度信息, 并且计算量很大, 因此并没有广泛应用; ③基于预测滤波的数据重建方法, 这类方法主要利用分频预测的思想, 用低频信息预测高频信息, 主要有f-x域预测滤波^[3]、f-k域预测滤波^[4]以及t-x域预测滤波。

基于信号处理的方法常见的是傅里叶变换和小波变换, 这两种方法对于线性同相轴能取得很好的重建效果, 但对于弯曲同相轴的效果很难令人满意。近年来发展起来的Shearlet变换与小波变换不同, Shearlet变换有各向异性基, 近似满足抛物线比例, 且在频率域严格局部化, 因而Shearlet变换能够更好地表示地震波场特征。与Curvelet变换类似, Shearlet变换是一种多尺度稀疏变换, 并且具有更好的方向性、稀疏性以及更加简单的数学结构。目前, Shearlet变换已经在图像处理领域得到广泛应用, 但尚未广泛应用到地震数据处理领域。在地震数据处理领域中, Curvelet变换在地震数据去噪^[5-6]与稀疏重建中得到了广泛的应用。

地震数据重建问题可以转化为一个稀疏优化问题, 压缩感知理论^[7-8]的提出, 使得在大量数据缺失的情况下仍然可以很好恢复出地震数据(即使在低于香农采样定理的情况下)。其重构的关键在于数据的稀疏性, 但是地震信号往往并不是稀疏的, 如果将这些信号转换到某些变换域, 使得信号在这些变换域内的大多数系数为零或者很小, 并且绝大多数重要的信号在变换域内的系数较大, 较大的系数保留了信号绝大部分信息, 使得信号在变换域是稀疏的。因此, 稀疏变换成为压缩感知理论应用的关键, 通过某些稀疏优化解法可以对缺失数据进行重构。HERRMANN等^[9]利用Curvelet变换通过L₁范数最优化解法对地震数据进行重建, 相对常规重建方法效果有很大的提升。Curvelet变换成为压缩感知关键的组成部分, 之后这种稀疏变换迅速发展。周亚同等^[10]采用离散余弦变换(DCT)与曲波字典组合进行稀疏重建。路交通等^[11]在Curvelet变换的基础上进行了改进, 取得了较好的效果。

本文对地震数据重建问题的求解过程可以理解为压缩感知理论范畴, 压缩感知的3个组成部分分别是稀疏变换、采样方式以及求解方法。首先, 采用Shearlet变换作为数据重建的稀疏变换代替了Curvelet变换, Shearlet变换比Curvelet变换有更好的稀疏性, 同时对地震数据采取jitter欠采样方式(jitter欠采样方式很大程度上模拟了野外地震数据采集中遇到的实际问题, 避免采用随机采样引起间隔过大, 同时又保证一定的随机性, 规则采样的假频在一定程度上比jitter采样更加严重, 对地震数据的重构影响更大), 与常规的规则采样和随机采样方法相比，有效地提高了重构精度; 其次, 将地震数据重建描述为图像修复的过程, 通过Shearlet变换的L₁范数最优化求解, 采用凸集投影算法(POCS)对地震数据进行重建, 有效地提高了重建效率; 最后, 对地震数据进行动校正预处理后进行重建, 动校正预处理的Shearlet变换的重建效果较好, 与传统的小波重建和Curvelet变换重建相比具有明显的优势。

1 Shearlet变换

2005年以来, 许多学者^[12-17]将复合小波理论和几何多尺度分析联系到一起, 通过特殊形式的具有合成膨胀的仿射系统构造了一种接近最优的多维函数稀疏表示方法——Shearlet变换, 有着比曲波和轮廓波更敏感的方向性和更好的稀疏表示性质。

在二维情况下合成仿射系统的形式为:

$\begin{array}{*{35}{l}} {{\varphi }_{AS}}\left( \psi \right)=\{{{\psi }_{a,s,t}}\left( x \right)={{\left| detA \right|}^{a/2}}\cdot \\ \psi \left[ {{A}_{a}}{{S}_{s}}\left( x-t \right) \right]:a,s\in Z,t\in {{Z}^{2}}\} \\ \end{array}$

(1)

式中:A和S为2×2可逆矩阵, |detS|=1。那么剪切波系统的基本单元复合小波为φ_AS(ψ), 由复合小波构成了如下Parseval框架(紧框架条件):

$\sum\limits_{j,l,k}{|\langle f,{{\psi }_{j,l,k}}\rangle {{|}^{2}}}=\|f{{\|}^{2}}$

(2)

式中:f∈L²(R²)。在上述Shearlet系统中, 各向异性对角矩阵A_a为Shearlet尺度矩阵, S_s为Shearlet剪切矩阵, 其中, a, s分别为尺度参数和剪切参数, 对∀a＞0, s∈R, 有:

$\begin{matrix} {{A}_{a}}=\left[ \begin{matrix} a & 0 \\ 0 & \sqrt{a} \\ \end{matrix} \right] & {{S}_{s}}=\left[ \begin{matrix} 1 & s \\ 0 & 1 \\ \end{matrix} \right] \\ \end{matrix}$

(3)

在上述系统下, 令ψ∈L²(R²)满足:

$\hat{\psi }\left( \xi \right)={{{\hat{\psi }}}_{1}}({{\xi }_{1}}){{{\hat{\psi }}}_{2}}\left( \frac{{{\xi }_{2}}}{{{\xi }_{1}}} \right)$

(4)

式中:为${\hat{\psi }}$的傅里叶变换; ψ₁为连续小波, ${{{\hat{\psi }}}_{1}}$∈C^∞(${\hat{R}}$)且supp${\hat{\psi }}$₁∈[-2, -1/2]∪[1/2, 2]; ${\hat{\psi }}$₂∈C^∞(${\hat{R}}$)且supp${\hat{\psi }}$₂∈[-1,1]。其中, 为复数域。

在上述仿射系统下, 由Shearlet母函数ψ, 各向异性膨胀矩阵A_a和剪切矩阵S_s构建了剪切波系统:

$\begin{align} & {{\psi }_{ast}}\left( x \right)={{a}^{-\frac{3}{4}}}\psi \left[ {{A}^{-1}}_{a}{{S}^{-1}}_{s}\left( x-t \right) \right] \\ & \begin{matrix} a\in {{R}^{+}} & s\in R & t\in {{R}^{2}} \\ \end{matrix} \\ \end{align}$

(5)

基于数学变换的信号处理, 其结果好坏主要取决于数学变换对有效信号的逼近程度和信号在变换域的稀疏性。

首先, 图 1a, 图 1b和图 1c分别展示了Meyer小波、Curvelet变换以及Shearlet变换的频率域剖分^[18]:Meyer小波基——正则化的小波基可以很好地表示各向同性数据; Curvelet基——具有方向性的紧支撑框架, 精细尺度可以很好地表示各向异性数据; Shearlet基——具有方向性的紧支撑框架, 精细尺度可以很好地表示各向异性数据; 且在连续系统和离散系统之间的转变也更加灵活自然, 是目前多尺度领域唯一可以将连续系统和数字系统统一处理的变换方法。

其次, 相对于傅里叶变换、小波变换、Curvelet变换在它们各自的域中得到的稀疏系数, Shearlet变换系数具有更加稀疏的特性^[19], 但在地震数据处理领域还没有得到广泛的应用。

2 稀疏变换域L₁范数约束的缺失道地震数据重建 2.1 结合Shearlet变换的POCS地震数据重建

假设完整的地震数据为A, y为缺失的地震数据, M为地震数据的一个缺失算子, 与数据缺失位置对应的元素为0, 其它元素为1。那么地震数据重建的问题可以归结为L₁范数最优化问题, 即:

$x=argmi{{n}_{x}}\|x{{\|}_{1}}或\|y-MA{{\|}_{2}}\le \varepsilon $

(6)

式中:A=S^-1x, 其中, S^-1为Shearlet逆变换, x为一组Shearlet系数。即将问题转换到Shearlet域重构出地震数据A, 反演出一组L₁范数最小的Shearlet系数x, 可以利用如下约束最优化问题来表示稀疏约束情况下重建问题:

$\left\{ \begin{matrix} x=argmi{{n}_{x}}\|y-M{{S}^{-1}}x\|_{2}^{2}+\lambda \|x{{\|}_{1}} \\ A={{S}^{-1}} \\ \end{matrix} \right.$

(7)

为了解决上述约束最优化问题, GERSZTENKORN等^[20]提出迭代重加权最小平方法(IRLS), BERG等^[21]提出SPGL₁算法以及各种阈值迭代方法求解。但是上述求解方法过于繁琐, 计算量很大, 利用凸集投影理论(POCS)^[22-23]求解此类问题使求解过程变得简单。结合POCS算法的Shearlet域稀疏约束地震数据重建基本步骤如下。

1) 初始化。

设置误差参数ε, 迭代次数N以及阈值参数δ;

初始化X_m=XM, X_i=0(X为完整地震数据, X_m为缺失数据);

令λ=1, δ为Shearlet系数的最大归一化系数。

2) 阈值迭代。

计算缺失数据R=M(X_m-X_i);

进行Shearlet变换x=S(R+X_i);

作阈值, Shearlet逆变换求重建后地震数据X_i=S^T[T(x)]。

3) 更新阈值迭代。

δ_i+1=λ_iδ_i, λ_i=ε^1i-1; 迭代求解, 直到δ＜λ。其中, 阈值函数为:

${{T}_{\lambda }}\left( x \right)=\left\{ \begin{matrix} x\left| x \right|\ge \delta \\ 0\left| x \right|＜\delta \\ \end{matrix} \right.$

(8)

为了加快收敛速度, 本文选用指数衰减的阈值代替线性阈值, 其中λ_i=ε^1/(i-1)。求解上述优化问题就是反演一组Shearlet系数x, 且该组数据满足经缺失算子作用后与原始数据在二范数约束下最接近。因此, 最终的输出就是这组Shearlet系数逆变换的结果, 代表了重建后的数据。

2.2 观测系统的采样方式

采样方式是压缩感知的重要组成部分, 直接影响到数据的稀疏性。地震勘探通常在地面观测, 采样方式有限, 常见的有规则采样和随机采样。由于地质条件等因素的限制, 往往不允许采取规则采样, 且规则采样产生的假频信息很多, 不利于后续的数据处理工作; 而随机采样能减少假频, 但是很难保证采样间隔, 容易出现采样点过于聚集或者过于分散的状况。当数据不存在随机性或者缺失过大时, Shearlet稀疏重建效果与Curvelet稀疏重建效果类似, 数据重构结果都不够理想。对于规则采样和随机采样数据重建, 很多学者进行了研究^[24], 这里不再赘述。HENNENFENT等^[25]引进了一种jitter采样方式, 这种方法兼具了规则采样方式和随机采样方式的优势, 既能控制采样间隔又能保持采样点的随机性。首先将目标区域等间隔采样, 然后以这些采样点为中心随机扰动, 为了说明jitter采样方式, 假设采样因子γ为奇数(γ=1, 3, 5, …), 总采样点是γ的倍数, 采样的样点数为n=N/γ取整数, 这里, 可以给出的jitter采样点为:

$\begin{align} & \begin{matrix} y\left[ i \right]={{f}_{0}}\left[ j \right] & i=1,2,\ldots ,n \\ \end{matrix} \\ & 且j=\underbrace{\frac{1-\gamma }{2}+\gamma \cdot i}_{确定值}+\underbrace{{{\varepsilon }_{i}}}_{随机值} \\ \end{align}$

(9)

离散随机变量ε_i是独立的整数, 且在采样间隔[-(ξ-1)/2, (ξ-1)/2]内均匀分布, jitter参数ξ(0≤ξ≤γ)表示粗糙网格下的jitter大小, 当ξ=γ时为最优jitter采样^[26]。

2.3 动校正预处理下的地震数据重建

野外采集地震数据时, 当地面水平、反射界面为平面且界面内介质均匀时, 反射时距曲线为一条双曲线, 不能直接反映地下界面的起伏情况, 只有在激发点处接收的t₀时间, 才能直接反应界面的真深度。因此, 我们可以对其进行动校正, 即把炮检距不同的道上来自同一界面、同一点的反射波到达时间, 校正为共中心点处的回声时间, 从而实现同相叠加。

由于Shearlet变换是一种多尺度变换, 它的结构元素除了包括尺度参数和平移参数, 还包括剪切参数, 所以Shearlet函数具有很好的方向选择性和各向异性, 可以很好地表示各向异性地震波场。基于Shearlet变换稀疏约束重建的关键在于地震数据在Shearlet域的稀疏性, 稀疏性越好意味着我们可以利用越少的系数对地震数据达到精确重建。因此, 为了使地震数据在Shearlet域更加稀疏, 利用Shearlet变换重建的地震数据做动校正预处理, 动校正之后的数据同相轴更加平缓, 在一定程度上降低了假频影响, 更重要的是, 相比于曲率过大的同相轴, 在Shearlet域中可以用更少的系数表示地震同相轴信息, 其稀疏性将会更好, 相应地, Shearlet域中的重构误差将会减小, 重建之后再反动校正回来, 从而达到更好的重建效果。WANG等^[27]将这种预动校正处理的方法与Curvelet变换结合对地震数据进行重建, 也取得了较好的效果。

3 模型及实际数据测试 3.1 模型数据测试

为了验证Sheratlet变换在地震数据重建中的效果, 先对一个合成地震数据采用jitter欠采样后进行重建。该地震记录共361道, 道间距25m, 采样间隔4ms。首先分别利用Curvelet变换迭代法、Shearlet变换迭代法对缺失地震数据重建, 其次对缺失的地震数据进行动校正预处理, 然后再利用Curvelet变换迭代法和Shearlet变换迭代对缺失的地震数据进行重建。

图 2显示了原始地震数据以及缺失比为50%, 75%的地震数据(设定最大连续缺失道为5道)。

图 3显示了采用不同方法对缺失比为50%的地震数据(jitter欠采样)进行重建的结果。下面比较不同方法的重建效果, 并分析产生的假频以及重建后对假频的压制效果。图 3a和图 3b是缺失比为50%地震数据(图 2b)经过Curvelet重建后的地震记录及差剖面, 可见, 重建取得了一定的效果, 但是同相轴连续性不好, 差剖面中同相轴剩余能量很强。图 3c和图 3d分别为经动校正预处理过后再经Curvelet变换重建后的地震记录及差剖面, 可以明显看出, 相比未经动校正处理的地震数据重建后同相轴变得光滑连续, 并且差剖面残余能量更少。图 3e和图 3f分别为经Shearlet变换重建后地震记录及差剖面, 相对于图 3a和图 3b, 同相轴几乎完美重构, 差剖面也少见残余能量, 而经动校正预处理(图 3g, 图 3h)后, 重建效果仍然好于未经动校正预处理过的重建效果。

为了进一步验证基于Shearlet变换对大量数据缺失情况下重建的效果, 对缺失比为75%的地震数据做同样的处理, 得到的结果如图 4所示。

从图 4可以很明显看出, 基于Shearlet变换的地震数据重建效果要好于Curvelet变换的重建效果, 并且经动校正预处理后的重建效果更好，尤其在同相轴曲率大的部分。这是由于稀疏变换时利用各向异性小波对同相轴逼近, 对同一数据而言, Shearlet的基函数可以取得比Curvelet的基函数更好的逼近效果, 而当同相轴曲率很大时, 逼近效果都有所欠缺, 经动校正预处理后的地震记录同相轴相对较为平滑, 能取得较好的逼近效果。

为了展示Shearlet变换重建对假频信息的压制效果, 选取了缺失比为50%的地震数据, 并对重建后的地震数据作f-k谱分析。图 5a, 图 5b分别是原始地震数据和缺失比为50%的地震数据的f-k谱, 可见, 缺失地震数据的f-k谱出现了严重的假频。采用不同方法对缺失比为50%的地震数据进行重建, 结果如图 6所示。从图 6可以看出, 经过Curvelet方法重建后, 假频信息得到了较好的压制, 但是部分频率的信息有所缺失; 而Shearlet变换重建对假频的压制效果明显更好; 经过动校正预处理的重建效果好于未经动校正处理的重建效果。

3.2 实际数据测试

为了进一步验证Shearlet变换稀疏约束地震数据重建在实际地震数据处理中的效果, 选取某实际地震记录(图 7a)进行重建, 缺失的地震数据(图 7b)缺失比为50%, 重建效果如图 8至图 11所示。

图 8与图 9分别展示了Curvelet变换与Shearlet变换重构数据的结果与图 7a的差剖面, 可见, Shearlet变换重构效果要好于Curvelet变换重构效果。图 10与图 11分别展示了经过动校正预处理的Curvelet变换与Shearlet变换的重建效果以及与图 7a差剖面, 同样可以看出, Shearlet变换重建数据与实际数据差值最小。

4 结束语

1) 本文采用POCS算法将Shearlet变换算子引进到地震数据重建中, 验证了Shearlet变换能够在jitter欠采样中地震数据缺失的情况下, 比Curvelet变换更能有效重构出地震数据; 从理论数据与实际数据的f-k谱中能够看出有效压制了由缺失道引起的假频, 经过动校正预处理的重构效果更佳, 为地震数据波场重构提供了新的技术思路。

2) Shearlet变换较Curvelet变换的数学表达更加简洁, 因此运算时间比Curvelet变换减少一半, 但要大规模应用到野外海量地震数据处理, 还需要进一步提高运算效率。

3) 目前野外地震数据大规模采用了三维采集, 将三维Shearlet变换应用于三维地震数据重建将是未来该方法研究的方向之一。