非凸Lp范数三维地震数据重建

刘 群 付丽华 张婉娟

(中国地质大学(武汉)数学与物理学院，湖北武汉 430074)

0 引言

地震数据为探测地质构造和地层分布提供了有效的地球物理信息。地震数据的质量直接影响地震处理和解释结果[1-4]。地震数据在空间上可能出现不规则采样[5-6]，造成有效信号的缺失，从而增加地震数据处理难度[7-8]。因此，对地震数据进行有效重建，使后续数据处理能够更好地刻画复杂地质结构非常必要[9]。随着地震勘探研究的不断深入,更多的数学方法被用于地震数据重建[10-12]。其中一类是基于滤波的重建方法，该类方法通过褶积插值滤波器实现[13-16]。另一类重要方法是基于变换域的地震数据重建，主要有Radon域变换[17-19]、Fourier域变换[20-22]、Curvelet域变换[23-25]和Shearlet变换[26]等，即基于数学变换理论，把数据转换到不同的域进行重建。基于波动方程的地震数据重建方法也得到广泛研究，如炮域延拓法[27]、炮检距域延拓法[28]、DMO法[29]和拓展成像条件法[30]，该类方法计算量大，运算耗时，且当地下信息未知或精度较低时会影响重建效果，所以实际应用较少。

近年来，基于降秩重建的方法引起广泛关注。Trickett等[31]和Naghizadeh等[32]提出基于Cadzow滤波的三维地震数据随机噪声衰减方法，建立了Hankel矩阵并对其进行奇异值分解，通过降秩压制线性干扰。Oropeza等[33]提出多道奇异谱分析(Multichannel Singular Spectrum Analysis,MSSA)将地震数据变换到频率域，然后对每个频率切片构造块Hankel矩阵，通过随机奇异值分解对块Hankel矩阵降秩，得到重建的地震数据。Kreimer等[34]将该框架直接运用于五维地震数据重建，通过高阶奇异值分解(Higher Order Singular Value Decomposition)进行降秩。Abma等[35]提出凸集投影迭代方法，在频率切片上进行傅里叶变换后，对幅值进行阈值处理，最后进行反傅里叶变换并插值。Jia等[36]提出正交矩阵匹配追踪Hankel重建方法(Orthogonal Matrix Pursuit HankelReconstruction,OMPHR)，利用秩1正交匹配追踪方法(Orthogonal Rank-one Matrix Pursuit)代替SVD分解进行降秩，大幅提升运算效率。Siahsar等[37]提出最优奇异值收缩法，根据随机矩阵理论产生最优低秩估计量，并利用衰减因子约束奇异值，提高低秩估计的鲁棒性[38]。Gao等[39]引入一种基于构建块Toeplitz矩阵的快速降秩方法，运用Lanczos双对角化算法代替截断奇异值分解操作，提升了计算效率。Ma[40]提出基于构建纹理块的低秩矩阵补全方法，通过三维纹理块映射将三维数据变换到二维空间进行地震数据重建。

对于稀疏信号的恢复和低秩矩阵补全，一般建立秩最小化模型并将该模型凸松弛到核范数最小化模型。相对于凸核范数松弛模型，非凸Lp范数松弛模型可以以较低的采样率恢复缺失数据，通过在迭代中设置权重约束奇异值，很好地保证了数据结构的低秩性[41-42]。Lu等[43]成功地将非凸Lp范数松弛模型应用到彩色图像恢复上，并取得较好的效果。

基于降秩理论，本文提出非凸Lp范数Hankel重建(Non-convex LpNorm Hankel Reconstruction，NLPHR)方法对三维地震数据进行重建。非凸Lp范数比核范数更接近秩函数，因此通过求解最小Lp范数问题可以得到更接近于真实值的解。此外，MSSA方法和OMPHR方法均要求给定秩的大小，通过多次重复实验选取合适的秩。本文方法不需要给定秩的大小，而是通过设置权重约束奇异值，迭代减小奇异值，保证了重建数据的低秩性。

1 基于非凸Lp范数的地震数据重建

1.1 三维地震数据低秩重建

考虑三维地震数据体D(t,x,y)∈RNt×Ny×Nx，其中Nt为时间采样点数，Ny和Nx分别为x方向和y方向地震道数目。对每道数据进行傅里叶变换，得到频率域数据D(ω,x,y)。取每一个频率切片Dω∈RNy×Nx构成矩阵

(1)

首先对Dω的每行构造Hankel矩阵，以第j(1≤j≤Ny)行为例，格式如下

(2)

式中：Lx=floor(Nx/2)+1，其中floor(·)为向下取整函数；Kx=Nx-Lx+1。再对每一行的hj(1≤j≤Ny)进行Hankel矩阵化，建立块Hankel矩阵

(3)

(4)

式中算子PH∶RNt×Ny×Nx→R(Lx×Ly)(Kx×Ky)表示将三维数据映射为二维块Hankel矩阵。理论证明，同相轴数量无缺失的或无噪的地震数据是低秩的[33]，数据的缺失会增加矩阵的秩[40,44]。

1.2 NLPHR方法

低秩矩阵补全问题的一般形式为

(5)

(6)

对于地震数据重建问题，考虑建立非凸Lp范数松弛模型

(7)

假设1：h(x)∶R+→R+在[0,+∞)是凹函数且单调递增，可以为非光滑函数；

假设2：d(x)∶Rm×n→R+是一个光滑函数，即其梯度为Lipschitz常数Lf，满足

(8)

定义1：h(x)∶Rn→R是凹函数，如果对任意的b∈Rn，v满足

h(b)≤h(a)+〈v,b-a〉

(9)

由非凸Lp范数h(x)=xp(0

(10)

根据假设1及定义1(式(9))，有

(11)

其中

(12)

由h(x)的增量性及超梯度的反单调特性，有

(13)

因为d(x)是一个Lf光滑函数，可将d(X)在X(k)处线性变换为

(14)

用式(11)右边项代替式(7)的第一项，用式(14)右边项代替式(7)第二项，得到式(7)解的更新值

(15)

求解式(15)等价于计算加权核范数的逼近算子，根据式(15)的非凸性及式(13),可得到式(15)的闭式解。根据

定理1[45]：对任意λ>0,Y∈Rm×n,0≤w1≤w2≤…≤wq,q=min(m,n)，Y=UΣVT是Y的奇异值分解，Sλw(Σ)=Diag[max(Σii-λwi,0)]。对于问题

(16)

其全局最优解可通过加权奇异值阈值(WSVT)算法得到

X*=USλw(Σ)VT

(17)

通过求解式(15)更新X(k+1)后，根据

(18)

基于NLPHR的三维地震数据低秩重建流程(图1)如下。

(1)首先将三维数据D(t,x,y)从时间域变换至频率域，得到频率域数据D(ω,x,y)，然后对每个频率切片D(iΔω,x,y)分别构建相应的块Hankel矩阵H。

(4)迭代停止后，从降秩得到的近似矩阵H*反对角线平均值获取频率域数据D*(ω,x,y)，将频率域数据变换回时间域， 得到时间域数据D*(t,x,y)。

图1 基于NLPHR的三维地震数据重建流程

2 实验

通过对三维模拟地震数据和实际地震数据分别利用NLPHR法与MSSA方法、OMPHR方法进行对比。衡量重建质量的指标为信噪比

(19)

2.1 模拟数据实验

测试选用的模拟数据包含三条线性同相轴，x方向和y方向各有32道，道间距为1m。共有64个时间采样点，采样间隔为2ms。图2a为原始模拟三维数据，其大小为64×32×32。每个块Hankel矩阵大小为256×289。图2b为随机缺失40%道集的数据，图2c为NLPHR重建结果。

为了更清楚地说明文中所提出方法的重建效果，分别取图2中三维数据的二维剖面(y=19m)进行显示(图3)。图3直观地展示了NLPHR法重建前、后地震数据的吻合程度。图4为OMPHR与MSSA方法重建效果对比，两种方法中秩的大小均选取为3。NLPHR方法重建数据的信噪比为40.0dB，耗时3.8s； OMPHR和MSSA方法重建数据的信噪比分别为34.1dB和29.2dB，耗时分别为1.5s和94.5s。图5为分别随机缺失10%、20%、30%、40%、50%道集情况下NLPHR法、OMPHR法以及MSSA方法重建数据的信噪比曲线，可以看出，NLPHR方法重建信噪比明显高于OMPHR法和MSSA法。

图2 基于NLPHR法的三维模拟数据重建效果

图3 基于NLPHR法显示的重建剖面(y=19m)

2.2 实际数据实验

实际三维叠后地震资料(图6a)参数如下：x方向和y方向各有64道，道距为50m，时间采样间隔为2ms，采样点数为256。随机缺失40%道集数据见图6b，采用NLPHR法重建的地震数据如图6c所示。图7是从该三维数据体(图6)中抽取的二维剖面(y=1000m)。从图6和图7的重构结果可知，基于NLPHR方法的三维数据重建方法能有效恢复缺失道集，较好地实现(含缺失道)地震数据的重建。图8为随机缺失40%道集情况下OMPHR法及MSSA法重构结果。OMPHR法和MSSA法的秩均取为20。由于数据体较大，实验中实行分块并行重建。将大小为256×64×64的数据均分成四块，每块大小为256×32×32，对四块同时进行重建。NLPHR法重建数据的信噪比为11.9dB，耗时81.7s；OMPHR法重建数据的信噪比为9.3dB，耗时66.6s；MSSA重建数据的信噪比为7.7dB，耗时203.2s。可以看出，NLPHR法较好地恢复了缺失道集，而OMPHR法和MSSA法均未能完全恢复缺失道集。图9为不同缺失率下三种方法重建数据的信噪比对比。可以看出，NLPHR法具有较高的输出信噪比，比OMPHR法和MSSA法更具有稳健性。表1为三种数据重建方法在不同缺失率下三维叠后地震资料上重建数据的信噪比及计算耗时对比。从表1可以看出，NLPHR法重建数据的信噪比高于OMPHR法和MSSA法，耗时比MSSA方法短。图10和图11分别为缺失40%数据下该三维叠后数据体在参数p和η不同取值时重建数据信噪比及耗时对比。由图可见，参数p在(0.5,1)范围内重建数据的信噪比较高；p=0.6时重建信噪比最高，且重建耗时相对较短。由图11可以看出，η在(0.1,0.8)范围内，重建信噪比随η呈现增大趋势；在(0.8,0.9)范围内呈下降趋势。图12为数据缺失40%时三维叠后数据在三种重建方法下的SNR收敛曲线及收敛时间。由图可知，三种方法均较快收敛，其中NLPHR法收敛最快，重建数据的信噪比高于OMPHR法和MSSA法，迭代时间比MSSA法短。

图4 重建剖面(y=19m)结果对比

图5 三维模拟数据在不同缺失率下不同方法重建数据信噪比对比

图6 基于NLPHR法的三维叠后地震数据重建结果

图7 基于NLPHR法的重建结果二维剖面(y=1000m)显示

图8 叠后数据不同方法重建结果的二维剖面(y=1000m)对比

图9 三维叠后数据在不同缺失率下、不同方法重建数据信噪比对比

表1 三维叠后数据在不同缺失率下三种方法重建数据信噪比及耗时对比

图10 随机缺失40%道集的叠后数据在参数p取不同值下重建数据信噪比(左)及计算耗时时间(右)曲线

图11 随机缺失40%道集的叠后数据在参数η取不同值下重建数据信噪比(左)及重建时间(右)曲线

图12 随机缺失40%道集的叠后数据重建数据的SNR (左)及收敛时间(右)曲线

3 结束语

MSSA法通过随机奇异值分解直接降秩，得到的解往往不是最优解，并且在低采样率下重建效果欠佳。针对此问题，文章提出非凸Lp范数Hankel重建方法(NLPHR)，对三维地震数据进行重建。该方法将地震数据每个频率切片转换成块Hankel矩阵，再运用NLPHR法对块Hankel矩阵进行降秩处理，从而实现地震数据低秩重建。三维模拟地震数据和实际地震数据实验证明，与MSSA法和OMPHR法相比，本文方法数据重建信噪比最高，耗时较少，效果最优。