2017, Vol. 60
随着油气勘探开发任务难度的加大，研究有效解决实际问题的地震地质理论技术变得越来越重要.不管是新油田还是老油田，增储上产是油田企业发展的根本问题.大型的整装储层几乎很少，大部分是小而零散的储层；还有一类就是非常规储层.面临这些储层勘探预测问题，从现有地震资料中，挖掘弱小的地震地质有效信息，具有重要的研究意义.
提高地震资料分辨率是地震勘探的永恒话题，对应的理论技术浩如烟海.利用理论技术提取地震弱信息已取得了一些研究成果.当然，地下地质体情况的复杂性决定了地震资料处理的复杂性，虽然各种理论技术取得了一定的效果，但是离有效揭示弱小地震反射信息提取要求还有很大的差距.利用地震资料研究揭示复杂而又相对较小地质体地震信息提取识问题，具有很大的研究潜力.
压缩感知理论(Donoho, 2006a, b；Candès and Wakin, 2008)是近年来引进的新的信号处理理论，在数据压缩传输领域得到了大量应用.在勘探地震领域，不管国外还是国内也有一些学者进行了大量的研究，压缩感知理论信号采集突破了传统采样的不足(Candès，2006)，利用变空间描述信号，直接采集包含数据全部信息的压缩数据，将信号的采样变成信息的采样，需要使用原始信号时，可以通过一个优化问题从压缩数据中恢复得到(Wadayama，2010).在满足信号可以压缩和系统与观测系统不相关两大条件下，从低分辨率观测中恢复出高分辨信号.
本文结合反射地震信号的特征分析(Han and Zhang, 2011)，在简要讨论压缩感知理论基础上，分析了多道联合稀疏建模思想，从稀疏基、稀疏度及稀疏重建法进行了较深入研究，提出了新的方法，推导出了新方法的计算公式，建立形成多道联合稀疏重建地震信息的新方法.
2 压缩感知理论简介 2.1 压缩感知的数学模型将N维实信号x∈RN×1在某组正交基{ψi}i=1N(ψi为N维列向量)下进行展开，即:
|
(1) |
其中展开系数θi≤x, ψi≥ψiTx，写成矩阵形式可以得到：
|
(2) |
这里ψ=[ψ1, ψ2, …, ψN]∈RN×N为正交基矩阵(满足ψψT=ψTψ=I)，展开系数向量θ=[θ1, θ2, …, θN]T.假设系数向量θ是K-稀疏的，即其中非零系数的个数K≪N，那么采用另一个与正交基ψ不相关的观测矩阵ϕ:M×N(M≪N)(这里ϕ的每一行可以看作是一个传感器，它与系数相乘，获取了信号的部分信息)，对信号x执行一个压缩观测:
|
(3) |
就可以得到M个线性观测y∈RM，这些线性投影中包含了重构信号x的足够信息(Baraniuk et al., 2010).
从y中恢复x是一个解线性方程组的问题，但从方程(3) 的角度看这是一个超定方程(Donoho, 2006a, b)，存在无穷多个解，将式(2) 代入式(3) CS信息算子ACS=ϕψ，可以得到：
|
(4) |
这样使得信号重构成为可能，通过求解一个非线性优化问题就能从观测y、观测矩阵ϕ和矩阵ψ中可以很好的重建信号x.
2.2 压缩感知的条件压缩感知信号重建(Dai et al., 2011)是要满足一定条件：首先，信号在ψ下具有稀疏性或可压缩性(Pennec and Mallat, 2005)，即信号需要在变换空间下的展开系数足够的稀疏；其次，系统与观测系统不相关.在这两个条件都同时满足时，就可以通过求解如下问题:
|
(5) |
获得一个唯一确定的解，即稀疏系数向量θ，就可以得到信号x=ψθ.在求取稀疏系数时θ有
|
(6) |
这就将式(5) 的优化问题变成了一个凸优化问题，可以方便地转化为线性规划问题求解.
2.3 压缩感知流程压缩感知理论流程(Zhou，2010)主要包括三步：第一步，信号的稀疏化表示；第二步，观测矩阵的设计；第三步，信号的重构.信号的稀疏化是压缩感知理论的一个重要前提，并且直接影响着信号感知的效率.由统计理论和组合优化理论可知:在满足重构条件时，通过选择合适的观测方式和重建算法，仅需要K+1次观测就可将N维空间的K-稀疏信号精确地重建.观测矩阵的设计是压缩感知的关键，观测矩阵设计中的两个关键内容就是观测波形和采样方式，设计的主要原则是: (1) 观测波形在理论上的最优性能，即ACS要具有良好的性质; (2) 观测波形的普适性; (3) 实用性，包括快速计算、低存储量、硬件易实现等.信号的重构是压缩感知的核心，压缩感知信号的重构可以实现从低分辨观测中恢复出高分辨信号，其付出的代价就是在信号重建时的软件代价.人们在信号重构中提取出了l0范数下的凸化压缩感知框架(Davies and Gribonval, 2009；Gholami，2014；Wadayama，2010)，lp(0<p<1) 范数下的凸化压缩感知框架，l1范数下的凸化压缩感知框架等重构的方法，并在实际应用中取得了一定的效果.
3 多道地震资料联合稀疏重建方法 3.1 联合稀疏建模利用压缩感知理论进行地震资料信息重建过程中，地震记录的稀疏表示是影响重构质量的关键因素.目前在压缩感知地震资料重构中只利用单一的基函数表示地震记录.但是地震记录同时包含地层的岩性、厚度、结构模式及边缘信息，是各种因素的综合响应，并且各种信息成分所占的比重不同.从压缩感知原理出发，根据地质体在纵向时间方向和横向空间方向地震响应变化特征，构造联合稀疏地震记录响应信息表示模型.
把一个信号在某一基上稀疏的概念扩展到一组信号的联合稀疏上，构成了联合稀疏模型.联合稀疏模型中，每个原信号包括两部分稀疏的公共部分(区域)和特有部分(局部).其中，所有的信号中都包含稀疏的公共部分，每一个原信号的特有部分都各不相同，并且稀疏的公共部分和特有部分都可以在某一个基上稀疏表示.
|
(7) |
其中Zc=Ψθc，Zj=Ψθj
|
(8) |
对原始信号而言，信号Zc就是稀疏的公共部分，并且信号Zc在稀疏基Ψ上的稀疏度为Kc；信号Zj就是稀疏的局部部分，并且信号Zj在稀疏基Ψ上的稀疏度为Kj.
3.2 多道地震资料稀疏重建到目前为止，关于多道稀疏重建问题还没有形成较系统和有效的方法，多道稀疏大多数是跟分块稀疏基相关的，而多道稀疏信号可恢复性能与信号结构本身的关系，还有待进一步研究.与一般稀疏信号不同，分块稀疏结构涉及的参数更多，与可恢复性能的关系也更加复杂.算法方面，目前求解多道稀疏重构问题的算法还不多，多数限于理论研究，真正更多更有效的算法还有待进一步提出.
本文结合反射地震信号，特别是叠后反射地震信号地质体和地震响应的特点，尝试提出了以下方法，该方法大致思想是：对任意点采用滑动椭圆窗口，椭圆的长轴为时间方向，短轴为空间横向方向.纵向采用联合稀疏重建方法，横向采用常规稀疏重建方法，把两者的重建结果进行加权相加.利用压缩感知理论，实现多道联合稀疏重建地震信号.
3.2.1 纵向地震信号的稀疏重建方法地震记录在纵向时间方向是反射稀疏和地震子波的褶积，包含了地震记录信息的公共(区域)部分和特有(局部)部分，这两部分可以联合表示为如下形式：
|
(9) |
Xt, l表示纵向联合稀疏重建结果，Xtc, l纵向联合稀疏区域(公共部分)重建结果，Xtj, l表示纵向联合稀疏局部(特有部分)重建结果.
3.2.2 横向多道地震信号的稀疏重建方法地震记录在横向空间方向反映了地层的横向信息，也包含了地震记录信息的公共(区域)部分和特有(局部)部分，这两部分可以联合表示为
|
(10) |
Xh, l表示横向联合稀疏重建结果，Xhc, l表示横向联合稀疏区域(公共部分)重建结果，Xhj, l表示横向联合稀疏局部(特有部分)重建结果.
3.2.3 空间横向数据的取值方式地面一个采样点的地震信息是地下菲涅尔带的地震响应信息，横向信息体现了地质体横向地震响应的加权叠加，地震记录的横向信息，体现了地震体的横向变化.因此横向地震数据的取值，应沿着地震同相轴变化方向进行取值.
3.2.4 联合重建地震信息具有多方向特性，为了充分挖掘地震信息特别是弱地震信息，综合纵横方向的信息进行联合重建，对任意点纵横向稀疏重建的信号进行联合，有：
|
(11) |
Xsth表示纵横向联合稀疏重建结果，Xsthc表示纵横向联合稀疏区域重建结果，Xsthj表示纵横向联合稀疏局部重建结果，c1, c2为加权系数.
4 多道联合稀疏重建的区域和局部部分的稀疏表示 4.1 小波分解稀疏化为了获得地震记录的区域和局部部分，采用小波分解的方法，把地震信号进行小波分解得到低频部分和高频部分，即区域部分和特有部分.分别对这两部分进行稀疏重建，得到每一道地震记录的区域和局部稀疏重建结果.
4.2 高斯基函数稀疏化将地震信号进行小波分解，就是对地震信号进行了小波域的稀疏分解.地震信号在小波域的分解，不可能完全稀疏表示.为了对信号最大程度地稀疏分解，许多人研究了字典稀疏分解.理论上信号字典稀疏分解能够最大程度地表示信号的全信息.但是，实际实现起来有一定的难度.存在很多实际中难以解决的问题，如信号的原子字典库如何建立，建立的字典库原子基函数集，到底多少合适，何种类型的基函数集能够表示信号的有效信息等等.本文结合叠后地震资料的特点，地震资料一般经过零相位化处理.高斯基函数就有一定代表性，因此选取高斯族函数作为稀疏基，其中[0, 1]组成的基函数就是一类高斯基函数.
对地震信号一次小波稀疏表示基础上，再利用高斯基函数做进一步稀疏表示.纵向时间分析采用渐变高斯基函数，横向空间方向采用(0, 1) 基函数，这类基函数对于突出断层的信息效果更好.对于渐变的地质体边缘，空间方向利用高斯基函数效果比较理想.
5 稀疏度估计压缩感知稀疏度选取与确定是决定信号重建质量的又一关键问题.在讨论稀疏度确定的理论方法后，结合弱地震信号的特征分析，研究了确定稀疏度的方法.
5.1 纵向稀疏度估计方法信号的稀疏或近似稀疏是根据确定的基分析讨论.信号可以表示为：
|
(12) |
|
(13) |
这里
|
因此有
|
(14) |
并且
|
根据匹配追踪算法的迭代条件
|
(15) |
|
(16) |
这样当fk*≤ρ, ρ为一常数，ρ<1且接近1，K被称作稀疏度.
信号的稀疏度是能够表示信号的非零系数的个数.稀疏重建时非零系数是按从大到小排序的.稀疏度大小就是截取的降序排列系数的多少.对于地震衰减信号，稀疏度小意味取非零系数个数少，重建结果主要是振幅较大的部分信号，扔掉了后面相对振幅较小的地震信号，而这些振幅较小的信号，恰恰又是我们需要提取的信号，因此由式fk*≤fk－1*判断选取稀疏度不合适.其实这也是压缩感知存在的问题，直接把强弱排序的信号进行进行稀疏重建，往往把局部的弱信号压制了.换言之，压缩感知矩阵运算，大的数据会淹没相对较小的数据，重建结果只能是保留相对较大的数据，压制了相对较小的数据.为了把相对较小的地震信号，通过压缩感知重建提取出来，最常用的方法是对数据进行一系列处理，如规格化处理.但是这些处理提高弱信号水平的同时，又放大了噪声信号，有时甚至弱信号的水平低于噪声水平.因此，直接利用上述方法进行稀疏度的判断选取是不合适的.
对于含噪声地震信号，稀疏度应受噪声背景影响.特别是弱的地震信号往往和噪声的能量水平相当，当稀疏度选取不合适时信号重建结果就达不到预期的效果.因此，为了提取弱信号，在局部部分稀疏度选取时，应该考虑横向相关约束条件.首先在满足fk*≤ρ条件下，求得不同稀疏度下的重建信号的自相关函数，对自相关函数进行横向互相关函数计算，自相互函数的互相关系数最大时的稀疏度就是要求取的稀疏度.
5.2 横向稀疏重建时稀疏度的选取对于叠后地震资料，由于地震信号是纵波垂直分量信息，加上地震资料处理的改变.使得横向地震信号的变化特征不同于纵向.稀疏度太小，地震记录横向连续性很差，人为造成了间断，稀疏度太大，又淹没了横向地层边界信息，如断层边界.因此，横向稀疏度的选取合适与否对于准确提取横向边界信息极为重要.断层边界和地层边界的横向变换特征不同，断层边界变化特征同相轴错段，不连续点组成面是直截面；地层边界点组成的面是渐变的面.要同时凸显这两种边界信息，选取能同时揭示这些边界信息的稀疏度是很困难的，顾此会失彼.鉴于以上分析讨论，本文采用经验试算方法确定横向稀疏度，首先对某个地区的实际地震资料，从井资料或其他已知地震资料获取地层和断层的边界信息，通过试算的办法，选取能够揭示断层和地层横向边界信息的稀疏度.
6 实际资料测试结果分析本文在讨论压缩感知基本理论基础上，针对叠后地震资料，研究了有效而实用的地震弱信息压缩感知处理方法，研究提出的地震资料多道联合压缩感知处理新方法，形成了一些新的理念和新的设计思想.为了验证研究新设计思想的实际应用效果，对实际资料进行了试算和验证分析.
图 1a为叠后地震剖面，图 1b到图 1d为小波分解稀疏重建结果，从重建结果看到，重建后地震资料弱的有效信息得到较好的揭示.同时，在重建过程中利用本研究的纵向稀疏度确定方法，把区域和局部分解结果进行了有效弱信息提取的同时，又较好地压制了一些无用的噪声成分.
|
图 1 小波分解联合稀疏重建结果 (a)原始叠后地震剖面; (b)区域部分重建结果; (c)局部部分重建结果; (d)联合重建结果. Fig. 1 Results of wavelet decomposition combined with sparse reconstruction (a) Original post-stack seismic profile; (b) Regional partial reconstruction results; (c) Local partial reconstruction results; (d) Combined reconstruction results. |
图 2为不同稀疏度情况的结果比较：图 2a为稀疏度取较小值时结果，取值较小时，凸显了信号水平或能量大的部分，压制或丢掉了信号水平低或能量弱的部分.图 2b为稀疏度取较大值时的结果，稀疏度取值较大时，当稀疏度等于原信号长度时，完全恢复了原信号，但是没有提取重建弱有效信息部分.
|
图 2 不同稀疏度结果 (a)小稀疏度; (b)大稀疏度. Fig. 2 The results of different sparseness (a) Small sparseness; (b) Large sparseness. |
图 3为横向不同稀疏度情况的结果比较：图 3a为稀疏度取较大值时结果，稀疏度取值较大时，虽然增大了地层横向连续性，但是淹没了地层的横向边界信息；图 3b为当稀疏度取较小值时，增强了地层横向信息特别是断层信息的揭示能力；图 3c为当稀疏度取较更小值时，地层横向信息过分增大，造成了横向假的边界信息.分析比较结果看到，稀疏度取值大小决定了信号的重建质量，特别是对弱信号的提取，选取合适稀疏度尤为重要.
|
图 3 横向稀疏度结果 (a)大稀疏度; (b)较小稀疏度K=7; (c)小稀疏度K=5. Fig. 3 The results of transverse sparseness (a) Large sparseness; (b) Small sparseness K=7; (c) Smaller sparseness K=5. |
稀疏基的类型和特性，是决定压缩感知稀疏重建地震信号质量的又一关键因素.本文对此进行了讨论和试算分析，稀疏的目的是把信号投影到众多稀疏基上，然后进行全信息重建.稀疏基的数量、类型基特征，又决定了投影信息的最优表示程度.限于研究内容，本文结合地震信号的特点，只对高斯基和(0, 1) 组成的稀疏基进行试算结果的对比分析，平滑稀疏基重建地震渐变信号效果较好，如砂体边缘信息重建，突变稀疏基重建突变地震信号能力较强，如断层信号重建，见图 4.
|
图 4 不同稀疏基结果 (a)高斯基; (b)正交基. Fig. 4 The results of different sparse base functions (a) Gaussian function; (b) Orthonormal basis function. |
图 5为单道和多道联合稀疏重建结果，图 5c为多道联合稀疏重建结果和单道联合稀疏重建结果差值结果.分析图 5结果看出，多道联合压缩感知稀疏重建的效果明显优于单道结果.
|
图 5 单道与多道联合稀疏重建比较结果 (a)单道稀疏重建结果; (b)多道稀疏重建结果; (c)差值结果. Fig. 5 Comparison results of single channel and multi-channel combined sparse reconstruction (a) Longitudinal sparse reconstruction; (b) Joint sparse reconstruction; (c) Difference result. |
地震数据稀疏重建，纵向采用联合稀疏重建方法，横向采用常规稀疏重建方法，把两者的重建结果进行加权相加.重建结果看到，重建后地震资料弱的有效信息得到较好的揭示.同时，在重建过程中利用本研究的纵向稀疏度确定方法，把区域和局部分解结果进行了有效弱信息提取的同时，一些无用的噪声成分得以较好地压制.
将地震信号进行小波分解，就是对地震信号进行了小波域的稀疏分解.对地震信号一次小波稀疏表示基础上，再利用高斯基函数进行进一步稀疏表示.纵向空间方向采用高斯基函数，横向空间方向采用(0, 1) 基函数，这类基函数对于突出断层的信息效果更好.对于渐变的地质体边缘，空间方向利用高斯基函数效果比较理想.
弱地震信号往往和噪声的能量水平相当，为了提取弱信号，在局部部分稀疏度选取时，求得不同稀疏度下的重建信号的自相关函数，对自相关函数进行横向互相关函数计算，自相互函数的互相关系数最大时的稀疏度就是要求取的稀疏度.
平滑稀疏基重建地震渐变信号效果较好，突变稀疏基重建突变地震信号能力较强，多道联合压缩感知稀疏重建的效果明显优与单道结果.研究方法通过实际资料验证方法的合理正确性，具有显著的应用效果.
| Baraniuk R G, Cevher V, Duarte M F, et al. 2010. Model-based compressive sensing. IEEE Transactions on Information Theory, 56(4): 1982-2001. DOI:10.1109/TIT.2010.2040894 | |
| Candès E J, Romberg J, Tao T. 2006. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory, 52(2): 489-509. DOI:10.1109/TIT.2005.862083 | |
| Candès E J, Tao T. 2006. Near-optimal signal recovery from random projections: Universal encoding strategies. IEEE Transactions on Information Theory, 52(12): 5406-5425. DOI:10.1109/TIT.2006.885507 | |
| Candès E J, Wakin M B. 2008. An introduction to compressive sampling. IEEE Signal Processing Magazine, 25(2): 21-30. DOI:10.1109/MSP.2007.914731 | |
| Dai Q H, Fu C J, Ji X Y. 2011. Research on compressed sensing. Chinese Journal of Computers (in Chinese), 34(3): 3425-3434. | |
| Davies M E, Gribonval R. 2009. Restricted isometry constants where lp sparse recovery can fail for 0 < < p≤ 1. IEEE Transactions on Information Theory, 55(5): 2203-2214. DOI:10.1109/TIT.2009.2016030 | |
| Donoho D L. 2006a. Compressed sensing. IEEE Transactions on Information Theory, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582 | |
| Donoho D L. 2006b. For most large underdetermined systems of equations, the minimal-norm near-solution approximates the sparsest near-solution. Communications on Pure and Applied Mathematics, 59(7): 907-934. DOI:10.1002/(ISSN)1097-0312 | |
| Gholami A. 2014. Non-convex compressed sensing with frequency mask for seismic data reconstruction and denoising. Geophysical Prospecting, 62(6): 1389-1405. DOI:10.1111/gpr.2014.62.issue-6 | |
| Han W G, Zhang J H. 2011. Theoretical study on characteristic of weak signal and its identification. Oil Geophysical Prospecting (in Chinese), 46(2): 232-236. | |
| Le Pennec E, Mallat S. 2005. Sparse geometrical image representation with bandelets. IEEE Transactions on Image Processing, 14: 423-438. DOI:10.1109/TIP.2005.843753 | |
| Tsaig Y, Donoho D. 2006. Extensions of compressed sensing. Signal Processing, 86(3): 549-571. DOI:10.1016/j.sigpro.2005.05.029 | |
| Wadayama T. 2010. Sufficient condition for perfect recovery in compressed sensing. The Brain & Neural Networks, 17(2): 63-69. | |
| Yin H P, Li J X, Chai Y, et al. 2014. A survey on distributed compressed sensing: theory and applications. Frontiers of Computer Science, 8(6): 893-904. DOI:10.1007/s11704-014-3461-7 | |
| Zhou C M. 2010. Research on signal reconstruction algorithms based on compressed sensing[Master's thesis] (in Chinese). Beijing: Beijing Jiaotong University. | |
| 戴琼海, 付长军, 季向阳. 2011. 压缩感知研究. 计算机学报, 34(3): 3425–3434. | |
| 韩文功, 张军华. 2011. 弱反射地震信号特征及识别方法理论研究. 石油地球物理勘探, 46(2): 232–236. | |
| 周灿梅. 2010. 基于压缩感知的信号重建算法研究[硕士论文]. 北京: 北京交通大学. http://cdmd.cnki.com.cn/Article/CDMD-10004-2010119669.htm | |