稀疏反演求解基追踪降噪问题的地震谱分解方法
周岩;韩立国;于江龙;孙慧秋;张盼
【摘 要】Spectral decomposition was described as a linear inversion problem.As this problem is underdeter-mined,it needs sparse inversion algorithm to solve the problem of Basic Pursuit De-noising (BPDN).The authors use Spectral Projected Gradient --L1 (SPGL1 )to increase the resolution of Inversion Spectral Decomposition (ISD),and further study its potential advantages.The results show that ISD based on BPDN (ISD--BPDN)has higher resolution in time-frequency representation,which can accurately distinguish the formation,and can detect precisely the presence of hydrocarbons.%谱分解描述为一线性反演问题，由于该问题欠定性，需使用稀疏反演算法，然后将该问题转化为基追踪降噪问题（BPDN），引入谱投影梯度（SPGL1）稀疏反演算法提高反演谱分解（ISD）分辨率，并进一步研究其潜在优势。结果表明：基于基追踪降噪问题的反演谱分解方法（ISD －BPDN）在烃类检测中有更高的时频分辨率，分层更准确，精细地指示了烃类的存在。
【期刊名称】《世界地质》
【年(卷),期】2016(035)002
【总页数】9页(P517-525)
【关键词】谱分解;谱投影梯度;稀疏反演算法;基追踪降噪
【作 者】周岩;韩立国;于江龙;孙慧秋;张盼
【作者单位】吉林大学 地球探测科学与技术学院，长春 130026;吉林大学 地球探测科学与技术学院，长春 130026; 国土资源部应用地球物理重点实验室，长春 130026;吉林大学 地球探测科学与技术学院，长春 130026;吉林大学 地球探测科学与技术学院，长春 130026;吉林大学 地球探测科学与技术学院，长春 130026
【正文语种】中 文
【中图分类】P631.443
传统的信号分析建立在傅里叶变换的基础上，由于傅里叶分析使用的是一种全局的变换，无法表征信号的时频局部特性和瞬时频率特征。地震反射信号是典型的非平稳信号，需要分析
的是信号的局部 （或瞬时）时频特征。因此在傅里叶分析的理论基础上，发展了一系列分析和处理非平稳信号的时频分析技术，也称为谱分解技术。在地震解释中谱分解被广泛用于确定薄层储层特征［1］、地层可视化［2］、与烃类有关的低频阴影［3］以及储层的精细预测［4］等领域。谱分解分辨率能力直接影响地震解释的水平，提高其分辨率对含油气储层研究有重要意义。
短时傅里叶变换（STFT）是一种发展较早的时频分析技术［5］，但该方法受所选时窗的限制，时窗过长会影响高频组分的时间分辨率，而过短又会影响低频组分的频率分辨率。常用的一种窗函数为高斯窗函数，此种短时傅里叶变换称为Gabor变换。后来又陆续发展了其他一些时频分析方法，如连续小波变换 （CWT）［6］、S变换［7］以及广义S变换［8］、Wigner-Ville分布（WVD）［9］和匹配追踪（MP）算法［10，11］等。但受信号截断效应和时间-频率分辨率互相制约的影响，这些方法的时间和频率分辨率受到限制。
为了同时提高谱分解的时间分辨率和频率分辨率，一种替代的思想是把地震信号谱分解看作一个反问题来求解。由Oleg提出［12］将地震信号看作是一系列单一频率的子波与一系列由脉冲组成的伪反射序列 （区别于真实地下介质反射系数，故将其称为伪反射系数）对应褶积
然后叠加得到。这样，从地震记录分解出一系列的单一频率子波后得到的伪反射序列就可以当作是时频谱［12］，这个过程可以转化为求解线性反问题Ax＝b的过程。解x跨度时间域和空间域，即Nt和Nf，Nf为频率采样率，Nt为时间采样率，而地震记录只有Nt，所以该反问题是欠定的。David用Ricker子波和L2-L1范数，通过迭代加权最小二乘［13］和快速软阈值［14］方法对拟合信号实现了这种反问题解谱分解的思想［15］。但是其计算效率和时频分辨率可以通过选择最优的反演算法来进一步提高。本文将该反问题转化为基追踪降噪问题 （BPDN），引入谱投影梯度（SPGL1）稀疏反演算法求解该问题，达到了很好的效果。
1.1 反演谱分解数学模型
经典的地震褶积模型表达式为：
式（1）中：s()t为地震记录；w()t为地震子波；r()t为地下介质的反射系数；n()t为随机噪声，*为褶积运算。
经典的地震褶积模型在地震勘探史上发挥了重要的作用，但却无法表示地震信号的频率随时间的变化规律。为此将经典地震褶积模型式 （1）进一步完善发展成非平稳地震褶积模型 （模型假设子波为时不变的）：
正则化最小二乘问题
式（2）中：wi()t为主频为下角标i对应频率的子波；ri()t为与wi()t对应的频率依赖的反射系数；N为参与计算的反射系数或子波向量的数量。
式 （2）的物理意义是指，地震记录s()t是由一系列不同频率子波wi()t和与之对应的频率依赖的地层伪反射系数ri()t褶积后相加，再加上随机噪声n()t组成。根据线性代数理论可将 （2）式写成矩阵与向量的乘积形式：
或
式（3）和（4）中：s＝b为地震记录；Wi为子波wi()t对应的褶积矩阵（选取Ricker子波作为母小波进行分析）；A＝（W1W2… WN）为子波褶积矩阵库；x＝（r1r2… rN））T为频率依赖的地层伪反射系数矩阵。
式 （4）中，当已知子波褶积矩阵库A和地震记录b时，便可以通过求解线性反演问题得到频率依赖的伪反射系数矩阵x。矩阵（r1r2… rN）的行和列分别代表时间和频率，所以将x转化为（r1r2… rN）的形式，这个数据是以时间和频率为坐标的矩阵形式，即可认为是反演得到的时频谱［12，15］。通常在地球物理反演中，由式 （4）定义的反演谱分解问题是一个欠
定性问题，为了降低其不确定性和多解性，得到稀疏的时频谱，就需要对时频谱x进行稀疏约束，进而将问题转化为基追踪降噪问题进行求解。
1.2 基追踪降噪问题
在最近几年，压缩感知（Compressive Sensing，简称CS）已经成为信号处理和优化领域的研究热点。在压缩感知理论中很多学者证明了在一般条件下，极小化L1范数得到的解也是其对应欠定性线性方程组可能的稀疏解。假设存在一个未知信号x，一个测量向量和一个测量矩阵，且A为满秩矩阵，同时满足下式关系：
因为测量向量b中的元素个数远小于未知信号x中的元素个数，所以从测量矩阵A和测量向量b中重构未知信号x就构成了一个欠定的线性反问题。基本压缩感知理论表明：如果未知信号x足够稀疏，并且测量矩阵A满足一定条件时，那么从测量向量b可以准确或完全重构出未知信号x。
为了从欠定性线性方程组 （5）中重构稀疏信号x，一种方法是在方程组 （5）所有解中寻最稀疏的解，即求解极小化L0范数问题：
式（6）中‖x‖0表示向量x中的非零元素个数。事实上，优化问题 （6）可高概率地从一个非常有限数量的随机测量数据b中完全地重构出稀疏信号x。L0问题是一个非凸问题，一般很难求解。一种可行的方法是用L1范数替代优化问题 （6）目标函数中的L0范数，即在压缩感知中被称为基追踪（Basis pursuit，简称为BP）问题［16］：
极小化式（7）中，L1范数在促进解的稀疏性中起着关键作用。事实上，在某些合理条件下［17］，优化问题 （6）和 （7）存在公共解。在很多实际应用中，当b含有噪声或x不是完全稀疏的，而仅仅是可压缩的，这时需要对式 （7）中的等式约束进行一定的松弛。在这种情况下，对式 （7）松弛方式为基追踪降噪（Basis pursuit de-noising，简称BPDN）问题［16］：
式 （8）记为BPσ，其中正参数σ是对数据中的噪声水平估计，在σ＝0时转化为式 （7）中的基追踪（BP）问题。
凸优化问题BPσ可看作是L1范数正则化最小二乘问题的表述。事实上，BPDN通常是应用在求解惩罚最小二乘问题上：
这是由Chen等［16］提出的该问题的表述，记为QPλ。
下面这个方程式：
具有明确的L1范数约束，通常称之为Lasso问题（LSτ）。参数λ与LSτ约束中的Lagrange乘子有关，也与BPσ约束中乘子的转置有关。恰当的选择参数σ，λ以及τ可使得BPσ，QPλ和LSτ的解一致，而且这些问题在某种程度上是等价的。
谱投影梯度（Spectral Projected-gradient）具有高效地解决一系列 LSτ问题的能力［18，19］。对于QPλ，它是一个被标量参数化的问题；两者 （LSτ，QPλ）重要的差别在于，LSτ的对偶解可以解决如何更新τ使得LSτ的下一个解更加接近于BPσ的解。
令xτ表示LSτ的最优解。单一参数方程：
给出了对每个τ≥0时LSτ的最优解。
对于求解凸优化问题BPσ，采取的方法是将其重塑为一个对单一变量非线性方程根的问题：
式 （12）定义了一系列正则化参数τk→τσ，此处xτσ为BPσ的一个解。换言之，参数τσ使得L
Sτ和BPσ得到相同的解。谱投影梯度算法的每次迭代中会产生一个对该变量的估计，由此变量估计所定义的凸优化问题的解可得到导数信息，可用于基于牛顿方法的根算法。
由式 （11）定义的方程φ，对每个正则化参数τ都可得到约束问题LSτ的最优解。经典的Pareto曲线定义了残差r的L2范数和解x的L1范数之间的最优权重。BPσ和LSτ是这同一条曲线的两个不同的特征。可以使用函数φ中的参数τ，将Pa-reto曲线参数化。可知对定义域内所有点，φ以及Pareto曲线是连续可微的，这一结果允许我们使用牛顿算法来寻非线性方程 （12）的根。
对于式 （12），牛顿根算法中的每次迭代，都需要一些τ时φ和φ′的估计，从而得到每次迭代的LSτ极小值。之前提到的基于牛顿迭代，算法可产生一系列正则化参数τk→τσ：
这样LSτ的相关解xτk收敛到xσ。
牛顿根算法的每一次迭代需要φ()τ的 （近似）估计，然后再进行对LSτ的极小化。使用谱投影梯度算法（SPGL1）可达到此目的。
2.1 与传统谱分解对比
图1是由雷克子波构成的模拟数据，包含5个不同频率和相位的子波。子波中心时间从上到下依次为50 ms、150 ms、270 ms、380 ms和420 ms；子波频率依次为65 Hz、45 Hz、20 Hz、30 Hz和30 Hz；相位 （不含延时引起的线性相位差）依次为0°、-90°、180°、0°和0°。

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