基于压缩感知的L1范数谱投影梯度算法地震数据重建
兰天维;韩立国;张良
【摘 要】随着油气勘探的发展,采集的数据规模与复杂度越来越大,对这些数据进行重建的精度与效率影响到后续地震资料的处理效果.常用于地震数据重建的压缩感知理论与重建算法各有精度与效率的优势,因此对于大规模、复杂地震数据,综合考虑重建精度与计算时间,提出了一种基于压缩感知理论和L1范数谱投影梯度算法(SPGL1)的地震数据重建方法.首先根据地震数据的缺失情况选择采样矩阵,然后在contourlet域中采用L1范数谱投影梯度算法重建缺失的稀疏系数,最后进行contourlet反变换实现地震数据的重建.合成地震数据实验结果表明,基于压缩感知和L1范数谱投影梯度算法重建的地震数据精度较好,计算效率高.通过实际地震资料处理,对比了相同稀疏变换基情况下常用的贪婪算法中的正交匹配追踪(OMP)、梯度投影稀疏重建算法(GPSR)及L1范数谱投影梯度算法(SPGL1)的应用效果,发现基于压缩感知的L1范数谱投影梯度算法鲁棒性较好,受噪声影响小,重建精度高,并且兼顾了计算效率的需求.
【期刊名称】《石油物探》
【年(卷),期】2019(058)002
【总页数】11页(P219-228,244)
【关键词】压缩感知;测量矩阵;contourlet变换;地震数据重建;贪婪算法;绕射波;SPGL1
【作 者】兰天维;韩立国;张良
【作者单位】吉林大学地球探测科学与技术学院,吉林长春130026;吉林大学地球探测科学与技术学院,吉林长春130026;吉林大学地球探测科学与技术学院,吉林长春130026
【正文语种】中 文
【中图分类】P631
随着我国油气勘探程度的加深,勘探区域的范围与地质构造的复杂程度愈发加大,采集到的地震数据的规模与复杂度不断增加。复杂、大范围的勘探区域与地形、障碍物等环境因素以及仪器等经济因素一起,加剧了地震数据的不规则性或稀疏分布,这将不同程度地影响地震数据的后续处理效果,因此需要对不完整采样下的缺失地震道进行插值处理,即地震数据重建。目前地震数据重建的方法主要分为3类:基于滤波器的策略、基于波动方程算子的方法以及基于
某种变换的重建技术。基于滤波器的重建技术利用褶积插值滤波器将不规则数据当做规则数据处理,常用的方法是预测误差滤波法[1-2],此类方法会造成较大误差,其结果也具有不确定性。基于波动方程的插值方法主要通过正演与反演算子来迭代求解一个反问题,如NMO、反DMO、方位角时差校正AMO等[3-4],此类方法需要地下结构的先验信息,计算量大,对采样率要求较高。基于某种变换的方法通常是对地震数据进行某种变换,然后在变换域对地震数据进行重建。主要有傅里叶变换、小波变换[5]、Randon变换[6]、curvelet变换[7]、Seislet变换、contourlet变换[8]等方法。
压缩感知理论[9-10]是一种以低于Nyquist采样率采样的理论,其本质是当信号具有稀疏性或者在某个变换域可以稀疏表示时,利用一个与稀疏变换基不相关的观测矩阵,将该信号从高维映射到低维,然后利用稀疏促进算法求解最优化问题,实现数据的高概率重构。地震信号可以通过某种变换进行稀疏表示满足了压缩感知理论应用的条件。基于压缩感知的地震数据重建常用的变换方法主要有离散余弦变换(DCT)、傅里叶变换、小波变换、curvelet变换和学习型超完备字典等。DCT与傅里叶变换是全局变换,无法有效识别地震信号局部时频特征;短时傅里叶变换时频局部化能力有限,无法满足不同频率与时间的信号分辨要求;小波变换虽然弥补了短时傅里叶变换的缺点,但是无法很好地描述二维信号中的线性特征;curvelet变换具有方向识
别能力,可近似最优地表示二维信号中的曲线特征,然而其冗余度较高,计算效率低;学习型超完备冗余字典的变换能够很好地稀疏表示信号,但迭代次数多,训练时间长。相较于傅里叶变换与小波变换,contourlet变换以其方向性、局部性与各向异性,可较好地处理信号的线性、曲线特征,善于处理信号的奇异性,捕捉地震记录中由直达波、反射波、绕射波等波场构成的轮廓信息,从而能够较清晰地描述信号中因复杂构造产生的各种波场。与学习型超完备字典相比,contourlet变换运算时间短,算法简单,更加符合实际应用中的计算效率要求。
在基于压缩感知的地震数据重建方法中,重建算法是重要的一环。近年来基于压缩感知的重建算法主要有3种,分别是贪婪算法、组合算法和凸优化算法。其中贪婪算法通过迭代寻一组与观测值匹配的最稀疏原子,实现信号重建,包括匹配追踪(MP)[11]和正交匹配追踪(OMP)[12-13]等。组合算法要求原始信号能够支持快速分组测试重建,包括组测试和数据流草图[14-15],该算法需要对观测矩阵进行进一步设计。凸优化算法是将非凸问题转化为凸问题求解,寻信号的逼近,该算法的一个极为突出的优点是,在目标函数为严格凸函数时,用此算法求解得到的局部极大值(极小值)就是全局最大值(最小值),并且所求结果只有一个。凸优化算法包括基追踪[16]、迭代阈值[17]、内点法和梯度投影法[18-19]。其中基追踪算法为凸优化算法中的基本算法,由于计算量大,基本上用于一维信号处理;迭代阈值法运算简单,但收敛速度慢;内点法一
般用于中小规模问题;投影梯度算法[20-22]应用简单并且适合大规模与复数域问题,其中SPGL1收敛速度较普通的梯度投影法快[19]。
综合以上3种重建方法,考虑各稀疏变换与重建算法的精度及计算效率,本文提出了一种基于压缩感知的地震数据重建方法,该方法以contourlet变换为稀疏变换基,以L1范数谱投影梯度法为重建算法。该方法应用于模型数据和实际地震数据的处理结果验证了其在大规模、复杂地震数据情况下的有效性。
1 方法理论
基于压缩感知的L1范数谱投影梯度算法地震数据重建方法由三部分组成:压缩感知理论、contourlet变换与SPGL1算法,其中压缩感知理论为主体,contourlet变换为压缩感知的稀疏变换阶段,SPGL1为压缩感知的重建阶段。
1.1 压缩感知理论
压缩感知理论为地震数据的压缩采样提供了可能。在该理论框架下,地震数据的采样表达式为:
(1)
式中:f∈Rn为原始地震数据;y∈Rm(m≪n)为采样后地震数据;φ∈Rm×n为采样矩阵。地震数据的重建即将y恢复为f。若f在某个变换域中可稀疏表示,即:
(2)
式中:φ为稀疏变换算子;x为f的稀疏表示系数。将(1)式和(2)式结合,那么压缩感知理论框架下的数据重建可以表述成:
(3)
式中:φH为稀疏逆变换算子。以感知矩阵A代替φφH,则(3)式可写为:
(4)
该问题的求解方法为:
(5)
此方法也称之为L0优化问题,但是由于离散性与L0范数的非凸性,它是一个NP难题。为了解决此问题,人们用L1替代L0[23]:
(6)
在压缩感知中,为了保证(4)式的可解性,对矩阵A进行了限定,即矩阵A应满足有限等距性(restricted lsometry property,RIP)[23]。矩阵A是否满足RIP,由φ与φH是否相干决定,两者不相干时,A大概率满足RIP。在地震勘探中,常常把缺失的地震道看做0,未缺失的地道看做1[24],由此组成的对角矩阵被作为采样矩阵φ时,φ与φH不相干,满足了矩阵A的限定条件。因此,本文将此对角矩阵作为采样矩阵。
1.2 contourlet变换
随着地震勘探的深入,采集到的信号中包括越来越多的反射波、绕射波等轮廓信息,但是傅里叶变换与小波变换无法稀疏表示轮廓线,难以适应包含较多轮廓信息的复杂地震信号,因此,DO[25]提出了contourlet变换,这是一种二维可采集信号内在几何结构的变换,具有灵活的方向性和各向异性。它用类似于轮廓段的基结构来逼近几何结构,基的支撑区间是长宽比随尺度变化的“长条形”结构,可以对分段函数实现最优近似。
在压缩感知理论下,稀疏表示地震信号的contourlet变换步骤如下:
1) 使用拉普拉斯金字塔滤波器(LP)对地震信号进行多尺度分解,捕获不同尺度上的边缘奇异点。每一级LP分解首先对上一级低频信息进行下采样产生本级低频信息,然后对此低频信息进行上采样并与上一级低频信息相减得到高频信息,最后对低频信息进行LP分解迭代即实现了多尺度分解。
2) 对由LP分解产生的高频信息使用方向滤波器组(DFB)进行方向分解,并连接同方向的奇异点,得到contourlet变换系数。如DO[25]提出的contourlet变换中的DFB,将双通道梅花滤波器组与剪切算子结合,采用简化的二叉树分解方法,可获得较好的方向信息。
图1为contourlet变换流程图[26]。
图1 contourlet变换流程
1.3 SPGL1算法
针对凸优化问题一般具有规模大、目标函数非光滑的特点,本文采用L1范数谱投影梯度算法(SPGL1)来解决凸问题,得到重建的信号。
(6)式又被称为基追踪问题(basis pursuit,BP),其中A为m×n的矩阵,y为m列的列向量。当数据含有噪声或数据不完整时,基追踪问题变为基追踪去噪(basis pursuit denoise,BPDN):
(7)
式中:正参数σ为对噪声水平的估计,σ=0时为BP问题的解。
BPDN仅为L1范数正则化最小二乘问题中的一种,它还包括Lasso问题,其表达式如下:
(8)
式中:τ为阈值。
BPDN常用在惩罚最小二乘问题上,其形式是:
正则化正交匹配追踪
(9)
其中,λ是一个与Lasso问题中的约束拉格朗日乘数和BPDN问题中的约束乘数的倒数有关的参数。参数σ,λ,τ的适当选择,使得BPσ问题、LSτ问题、QPλ问题的解在某些情况下一致。在σ已知的情况下,λ可用同伦算法求取[19]。
因此,SPGL1算法的思想是将BPDN问题转化为一系列的Lasso子问题,通过谱投影梯度法(SPG)求解Lasso问题,以得到BPDN问题的解。
令xτ为LSτ问题的最优解,对于任意τ≥0,Lasso问题的最优值为:
(10)
Pareto曲线定义了残差r的L2范数与解x的L1范数之间的平衡,BPDN与Lasso是同一条曲线的两个不同特征,可以用参数τ参数化Pareto曲线。对于定义域内所有的点,ψ都是连续可微的,这样可以使用牛顿法求解(11)式:
(11)
(11)式定义了一系列正则化参数τk→τσ,可表示为:
(12)
式中:ψ′(τk)为ψ(τk)的导数。这样Lasso问题的相关解xτσ收敛到xσ,参数τσ使得BPDN与Lasso得到相同解。
SPGL1算法流程如下:
1) 初始化参数,步长α0∈[αmin,αmax],充分下降常数γ∈(0,1),初始迭代x0←Pτ[x],初始残差r0←y-Ax0,初始梯度方向g0←-ATr0,l←0,整数线性搜索步长M≥1。
2) 计算对偶间隙δl←‖rl‖2-(yTrl-τ‖gl‖∞)/‖rl‖2,收敛则退出,转步骤6)。
3) 起始步长α←αl,预备搜索迭代并更新相应残差,若则退出,否则α←α/2,重复步骤3)。
4) 更新迭代若ΔxTΔg≤0,则更新步长αl+1←αmax,否则αl+1←min{αmax,max[αmin,(ΔxTΔx)/(ΔxTΔg)]},l←l+1。
5) 返回步骤2)。
6) 输出xτ←xl,rτ←rl。
其中,右上角标T表示矩阵的转置;为x的近似解;为下的残差;Δx为xl+1与xl之差,Δg同理;l,j,k均为下标索引值;←为赋值号,表示将箭头右边符号表示的数值赋予左边符号。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。