基于压缩感知的非规则地震勘探观测系统设计与数据重建

周 松,吕 尧,吕公河,舒国旭,石太昆,霍守东

(1.中石化石油工程地球物理有限公司,北京100020;2.中国科学院地质与地球物理研究所,北京100029)

基于压缩感知的非规则地震勘探观测系统设计与数据重建

周 松1,吕 尧2,吕公河1,舒国旭2,石太昆2,霍守东2

(1.中石化石油工程地球物理有限公司,北京100020;2.中国科学院地质与地球物理研究所,北京100029)

提出了基于压缩感知的非规则观测系统优化设计和地震数据重建方法。基于贪心序贯策略,利用逐个选取采样点位置的方法来降低感知矩阵的最大互相关值,从而完成非规则观测系统的优化;利用L0正则化和L1正则化混合迭代的方法求解欠定采样矩阵,重建地震记录,获得良好的效果。进行了TFT工区实际观测系统的设计,鉴于接收电缆长度的限制,添加最大道间距和最小道间距作为观测系统设计过程中的约束,得到了符合实际勘探要求的非规则观测系统的优化设计。该方法为三维高密度地震勘探的大规模开展提供了有效的方案。

压缩感知;稀疏;非规则;观测系统设计;数据重建

随着我国油气资源勘探向中西部发展,地球物理的研究也逐渐转向地表构造复杂的区域。在这些区域进行勘探,不论是资料采集,还是后期的资料处理和解释,都对地球物理人员提出了更高的要求。

对于地震勘探来说,在复杂地区施工,往往会遇到峭壁、河流、沟壑、村庄和工业区等不易甚至无法布置检波器的情况,从而给数据的采集带来了一定的困难,甚至导致无法采集的情况发生。如果因为这些原因导致原始数据的大量缺失,就会严重影响整个工区的勘探质量和成像效果,从而给资料解释带来困难。

经典信号采集理论一般都会遵循Shannon-Nyquist采样定理,即:若要使采集到的数据能不失真地保持原信号中的信息,采样频率必须是原信号频带宽度的两倍以上。地球物理的野外资料采集也都遵循该定理。为了描述目标体的精细特征,往往需要采集更精确的信号,因此,需要通过在空间上加密采样点数和在时间上减小采样间隔来实现。这样的采集方式大大增加了采集的数据量,从而导致采集成本的急剧增高。

基于信号稀疏性的采样理论[1-2]的压缩感知(Compressed Sensing,CS)在图像压缩、无线通信、模式识别和医疗成像等研究领域受到了格外的关注。2006年,CANDES等[3]提出基于压缩感知的信号恢复技术,大大推动了压缩感知在地球物理勘探方面的应用。

压缩感知技术的发展和应用基于大部分信号都可以进行稀疏表示。在稀疏域对信号进行重建,可以获得原始的真实信号。压缩感知技术突破了Shannon-Nyquist采样定理的局限,大大节省了采集数据所占的空间,成为提高地震勘探效率的重要方法。

目前,压缩感知理论在地球物理勘探领域得到广泛应用[4-15]。基于压缩感知理论,HERRMANN等[4]提出了一种多震源的采集框架;MOLDOVEANU[5]探索了海上数据采集的随机观测方法;MOSHER等[6-7]提出并完善了一种在观测系统设计中基于约束条件的更优选择激发点和接收点位置的非均匀采样方法;陈生昌等[8]提出了一种针对地球物理数据的高效采集方法,并通过归一化各列数据的感知矩阵的互相关系数最小作为约束条件来设计最优的采集方案。

基于压缩感知的地震采集技术,获得的数据需要进行重建[9-10]。HENNENFENT等[11]提出了利用抖动欠采样进行波场重建的方法,并基于曲波变换进行了地震数据的重建;WU等[12]提出了Dreamlet域地震成像和数据重建的方法;LI等[13]提出了基于插值压缩感知的数据重建方法(interpolated Compressive Sensing,interpolated CS),设计了非均匀优化采样观测系统(Non-Uniform Optimal Sampling,NUOS),并通过地震数据重建,分别在海上和陆上取得了良好的试验效果;王华忠等[14]详细讨论了压缩感知技术在地震勘探应用中的问题和难点,从不同方面展示了压缩感知的有效应用。本文通过降低感知矩阵的最大互相关值,从而进行非规则观测系统的优化设计;并利用L0正则化和L1正则化混合迭代的方法进行欠定采样矩阵的求解,完成了地震数据的重建。

1 基于压缩感知的非均匀观测系统设计方法

压缩感知理论的基本假设是目标信号具有稀疏性或者可压缩性,也就是说,目标信号或者其在某个变换域中只有有限个成分不等于0(对应稀疏性),或者只有有限个成分远大于0(对应可压缩性)。若目标信号有K个成分不等于0,则称该信号为K稀疏的。

假设具有稀疏性或者可压缩性的目标信号为x,稀疏变换采用傅里叶变换F,则正交稀疏基为傅里叶基函数或分窗傅里叶基函数,于是有:

(1)

式中:s为信号x在傅里叶域的稀疏表示;(·)H表示共轭转置。

这样采集得到的数据可以看作是采样函数或采样矩阵与目标信号相乘的结果。记采样矩阵为Φ,当进行满采样时,有Φ=I,其中,I表示单位矩阵。当不能满足满采样条件(如采样数不足)时,Φ为I抽出的若干列组成的矩阵,只有采样位置对应的列向量被保留,则有采样数据:

(2)

式中:Ψ=ΦFH记为感知矩阵。

(3)

JAMAIL-RAD等[16]指出,感知矩阵列向量间的最大互相关值即为非规则采样归一化频谱的最大非零频率振幅。在频率域,采样数据y的频谱是感知矩阵Ψ的频谱和目标信号的稀疏表达s的频谱之间卷积作用的结果。μ就是因为非规则采样造成傅里叶基(稀疏基)的正交性被破坏所引起的最大频谱泄漏,压制μ使感知矩阵Ψ的频谱近似Delta函数,否则,感知矩阵Ψ的频谱出现多个峰,与稀疏表达s卷积作用后,将在频率域(稀疏域)产生“假频”噪声。CANDES等[17]指出,μ越小,则信号在非规则采样后能够重建的概率就越高。由于Ψ=ΦFH,其中F是固定的,因此,在采样数目不足时,可以通过改变采样矩阵Φ,即改变非规则采样点的分布来降低最大互相关值,从而优化采样矩阵,提高目标信号重建的概率,即:

(4)

JAMAIL-RAD等[16]给出了μ达到下界时的一种确定性采样方式,然而这种采样方式只在一些特定的条件下才能成立。注意到该优化问题是非凸的,具体优化时只要能够找到局部最优解就足够满足需求[18-19]。这里我们采用贪心序贯策略来优化采样点的位置,进而构建采样矩阵。具体策略如下:

1) 首先确定目标信号的网格点数目和间距,将所有目标信号的网格点作为候选采样点;

2) 遍历所有候选采样点,计算将其加入采样矩阵后的μ值,取使μ值最小的候选点作为新的采样点,更新采样矩阵;

3) 重复步骤2)的操作,直到采样点数达到规定数目;

4) 在初步选定规定数目的采样点位置后,利用时间抖动(Jitter)方法[18-19],以(4)式为目标,进行采样微调,最终确定采样矩阵。

以一维简单情况为例,假设在751个网格点上,选取375个点(近50%)作为采样点。当规则采样时,采样点分布及其归一化后的傅里叶谱如图1所示。可以看出,欠采样的规则采样函数会产生很强的频谱泄漏,造成假频混叠。同样的网格和采样点数,当随机采样时,采样点分布及其归一化后的傅里叶谱如图2所示。可以看到,随机产生的非规则采样降低了采样矩阵与傅里叶矩阵的相关性,没有假频产生,但是在整个频率范围都出现了很低的频谱泄漏,最大幅值为0.7949×10-1。保持网格和采样点数目不变,采用贪心序贯策略获得优化的非规则采样,采样点分布及其归一化后的傅里叶谱如图3所示。图3 中采样点非均匀分布,压制了假频,降低了整个频率范围内的频谱泄漏,频谱泄漏的最大幅值0.5382×10-1,比随机采样的频谱泄漏小。

由压缩感知原理可知,更小的频谱泄漏表明由非规则采集而造成的漫布在整个频带内的杂乱噪声更弱,这样采样频谱就更加逼近规则满采样时的频谱,更有利于目标信号的重构,优化后的非规则采样比随机采样所得到的数据能够完美重建的概率更高。当采样点数较少或者采样率较低时,这种差异更加明显。而当采样点较多时,采样条件往往会受其它因素如地形、实际施工条件限制等影响,此时可在采样点优化设计过程中加入这些约束,从而得到满足施工要求的采集方案,这是随机采样很难实现的。另外,对于某些二次施工工区,采样设计时还需要考虑其它指导目标如地下照明、数值模拟恢复信号的信噪比等,此时可在优化目标(4)式中加入其它约束项,然后优化完成满足要求的非规则采样点设计,这也是随机采样很难实现的。

图1 规则采样点分布(a)及其频谱(b)

图2 随机采样点分布(a)及其频谱(b)

图3 非规则优化采样点分布(a)及其频谱(b)

(5)

式中:Fx和Fy分别为x方向和y方向上的傅里叶变换矩阵。

引入向量化操作[16],将二维目标信号转化为一维信号。向量化算子:

(6)

表示将任意矩阵A各列叠加排列成一个很长的列向量a。向量化算子具有如下性质:

(7)

(8)

(9)

于是,可以通过最小化最大互相关值来确定采样算子G,即:

(10)

采用贪心序贯策略解决该问题时,需要增加候选网格的维度。考虑到非规则采样在滚动排列上的复杂性,实际施工的试验区采用了非滚动全检波点接收的方案,降低了炮点和检波点的联系。这里,将炮点和检波点的分布进行独立设计,减少了计算量。同时,假设没有地下构造的先验信息,仅凭目标信号的采样网格精度和采样数目来约束设计方案。

在实际观测系统设计优化时,需要根据实际情况添加一定的约束条件。特别是在设计检波点位置时,考虑到设备的实际情况,需要添加约束条件。第一是检波线距非规则,检波点非规则地分布在这些检波线上,而不是所有检波点都完全分散分布;第二是相邻检波器之间的最大距离受实际施工时电缆线长度的限制。考虑到这些因素,在设计检波点位置时,先优化检波线的位置,以(4)式为优化目标;然后固定检波线的位置,将M个采样点位置平均分布在Nx×Ny的网格上,然后利用Jitter采样[17]的方式,让采样点在初始位置附近抖动取值,以(10)式为目标进行迭代优化。每次迭代只抖动其中一个采样点,而其它点保持不变,取使μ最小的位置为新的采样点,更新该点位置,然后选择下一个点进行抖动,进行下一次迭代,直到抖动任何一个采样点均不能使μ变小为止[18]。为加速迭代收敛,可以先让各采样点在平均分布位置附近随机抖动,然后以此为初始条件进行迭代优化,同样可以达到满足要求的优化设计结果。

对炮点进行设计时,考虑最大炮间距和最小炮间距,在炮间距限制范围之内进行优化设计,使得非规则炮点的位置在目标网格上完全打散分布。此时可以根据序贯思想,利用贪心算法,逐个添加炮点位置以使得当前的μ最小。每步迭代先固定之前已有的所有采样点,然后在剩余未被选中的网格点中,选出使μ最小的网格点作为采样点。要加速运算,可以控制变量,每次只在一个方向取得最小值,然后换方向接着搜索,直到收敛。在所有M个采样点全部确定之后,再以(10)式为目标,进行全局微调,最终得到炮点的优化非规则采样设计。

2 TFT工区的实际应用

我们利用TFT工区的一块区域进行实际工区的非规则观测系统设计,所选工区范围为14km×8km,检波点布设到整个工区范围,炮点区域为10km×2km。规则检波点设计为32条接收线,线距240m,每条线456道,道间距30m。非规则设计是在原始规则检波点设计的基础上,保持大致的覆盖范围,但是只取总采样点数的75%,即线数32条不变,平均每条线取342道,平均道间距40m,共10944道。重建的目标网格为42条接收线,线距180m,每条线911道,道间距15m。优化后的检波点分布如图4所示,每个蓝色的十字型代表一个检波点的位置,最小点距为15m,最大点距为75m。同时,线距也进行了非规则优化设计,最小线距为150m,最大线距为705m。这里由于检波线数较少,为保证重建质量,降低频谱泄漏,采样点没有落在实际的整数网格上。检波点位置归一化后的频谱如图5所示,在垂直于检波线方向上由于检波线条数较少,频谱泄漏较强,最大值为0.0930,但仍然要比规则采样导致的假频小很多。

实际工区的非规则炮点设计为将1800个炮点分散在10km×2km的范围内,目标重建网格为炮点距15m,炮排距90m。优化后的炮点分布如图6所示,每个红色的方框代表一个炮点位置,其归一化后的频谱如图7所示,可以看到,频谱泄漏很小,最大值为0.0444。

图4 优化设计的非规则检波点分布

图5 优化设计的非规则检波点频谱

最终得到的该工区优化的非规则观测系统如图8 所示。

图6 优化设计的非规则炮点分布

图7 优化设计的非规则炮点频谱

图8 优化设计的非规则观测系统(红色表示炮点位置,蓝色表示检波点位置)

3 地震记录的模拟与重建

对于地球物理问题来说,基于压缩感知的数学重建模型可表示为:

(11)

式中:x∈RN表示原始地震数据;Φ表示采样矩阵;y∈RM表示观测地震数据。这里,M

(12)

式中:Ψ表示感知矩阵,Ψ=ΦCH。

对此,可以通过求解反问题:

(13)

来获取稀疏系数s,进而通过x=CHs求取原始数据x。但是,由于采样矩阵是一个欠定矩阵,可知上述问题是一个典型的NP-hard问题。对此,需要对上述问题增加约束条件,本文利用L0正则化和L1正则化混合迭代的方法求解上述问题。首先通过软阈值算法求解L1范数约束问题:

(14)

通过对(14)式的求解,可在进行下一步算法前获得一个好的迭代初值,并估算信号的稀疏水平。然后,通过硬阈值算法求解L0范数约束问题:

(15)

反复迭代,直至获得最终结果。

为了验证该算法重建地震记录的有效性,采用主频为30Hz的Ricker子波合成500道的单炮地震记录,道间距10m,每道3001个采样点,采样间隔为2ms,如图9a所示。图9b为采样50%地震道后的不完整地震记录,图9c为使用本文算法恢复的地震记录,图9d为本文算法恢复的地震记录与原始模型记录的误差。

图9 地震记录处理结果a 原始地震记录; b 随机缺失的地震记录; c 重建结果; d 重建结果与原始记录之差

从图9c的重建地震记录与原始记录(图9a)的对比以及两者的误差(图9d)可以看出,两者同相轴清晰,连续性相同,振幅一致性准确,重建误差小。因此,该算法具有良好的数据恢复效果。

图10是该地震记录的频率-波数谱。从原始地震记录和重建地震记录的频率-波数谱对比可以看出,重建数据的频率-波数谱与原始记录的频率-波数谱非常接近,两者的误差很小(图10d),进一步表明了该算法良好的重建效果。以上实验均表明,利用基于压缩感知理论的重建方法,可以较好地恢复缺失的地震数据。

图10 地震记录处理结果频率-波数谱a 原始地震记录频率-波数谱; b 随机缺失的地震记录频率-波数谱; c 重建结果的频率-波数谱; d 重建结果与原始地震记录频率-波数谱之差

4 结论与建议

基于压缩感知理论,本文进行了地震勘探非规则观测系统的优化设计和地震数据的模拟重建。通过降低感知矩阵的最大互相关值的方法优化设计采样点的位置,并应用于实际工区,该过程通过添加部分约束条件来使得最终的设计在实际可行的情况下达到最优解。利用L0正则化和L1正则化混合迭代的方法进行地震数据的重建,取得了良好的效果。

随着勘探精度的提高,野外地震资料的采集都要求能够达到“两宽一高”的标准,这就给野外采集带来了相当大的工作量。基于压缩感知的非规则观测系统的地震数据采集和信号重建,可以达到“相同野外采集工作量的条件下获得更高分辨率的地震资料”或者“获得相同分辨率的地震资料需要较少的野外工作量”的目的。

因此,利用压缩感知技术进行地震数据的采集和处理,可以大大降低野外资料采集的工作量,对目前“两宽一高”要求的地震勘探具有重要意义。

致谢:感谢中石化西北石油分公司在正演模型建立方面给予的大力支持,感谢中石化休斯敦研发中心的罗明秋博士在目标模型的建立和观测系统的设计方面给出的一些好的建议！

[1] DONOHO D L.Compressed sensing[J].IEEE Transactions on Information Theory,2006,52(4):1289-1306

[2] CANDES E J,ROMBERG J,TAO T.Robust uncertainty principles:exact signal reconstruction from highly incomplete frequency information[J].IEEE Transactions on Information Theory,2006,52(2):489-509

[3] CANDES E J,ROMBERG J.Quantitative robust uncertainty principles and optimally sparse decompositions[J].Foundations of Computational Mathematics,2006,6(2):227-254

[4] HERRMANN F J,ERLANGGA Y A,LIN T T Y.Compressive simultaneous full-waveform simulation[J].Geophysics,2009,74(4):A35-A40

[5] MOLDOVEANU N.Random sampling:a new strategy for marine acquisition[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:51-55

[6] MOSHER C C,KAPLAN S T,JANISZEWSKI F D.Non-uniform optimal sampling for seismic survey design[J].Expanded Abstracts of 74thEAGE Annual Conference,2012:333-336

[7] MOSHER C C,LI C,MORLEY L C,et al.Non-uniform optimal sampling for simultaneous source survey design[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:105-109

[8] 陈生昌,陈国新,王汉闯.稀疏性约束的地球物理数据高效采集方法初步研究[J].石油物探,2015,54(1):24-35 CHEN S C,CHEN G X,WANG H C.The preliminary study on high efficient acquisition of geophysical data with sparsity constraints[J].Geophysical Prospecting for Petroleum.2015,54(1):24-35

[9] 冯飞,王征,刘成明,等.基于Shearlet变换稀疏约束地震数据重建[J].石油物探,2016,55(5):682-691 FENG F,WANG Z,LIU C M,et al.Seismic data reconstruction based on sparse constraint in the Shearlet domain[J].Geophysical Prospecting for Petroleum,2016,55(5):682-691

[10] 崔永福,郭念民,吴国忱,等.不规则观测系统数据规则化及在相干噪声压制中的应用[J].石油物探,2016,55(4):524-532 CUI Y F,GUO N M,WU G C,et al.Regularization of irregular geometry seismic data and its application in the coherent noise suppression[J].Geophysical Prospecting for Petroleum,2016,55(4):524-532

[11] HENNENFENT G,HERRMANN F J.Curvelet reconstruction with sparsity-promoting inversion:successes and challenges[J].Expanded Abstracts of 69thEAGE Annual Conference 2007:231-233

[12] WU R S,WU B Y,GENG Y.Imaging in compressed domain using dreamlets[J].Beijing 2009 International Geophysical Conference and Exposition,2009:212

[13] LI C B,MOSHER C C,KAPLAN S T.Interpolated compressive sensing for seismic data reconstruction[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-6

[14] 王华忠,冯波,王雄文,等.压缩感知及其在地震勘探中的应用[J].石油物探,2016,55(4):467-474 WANG H Z,FENG B,WANG X W,et al.Compressed sensing and its application in seismic exploration[J].Geophysical Prospecting for Petroleum,2016,55(4):467-474

[15] 陈祖庆,王静波.基于压缩感知的稀疏脉冲反射系数反演方法研究[J].石油物探,2015,54(4):459-466 CHEN Z Q,WANG J B.A spectral inversion method of sparse-spike reflection coefficients based on compressed sensing[J].Geophysical Prospecting for Petroleum,2015,54(4):459-466

[16] JAMAIL-RAD H,KUVSHINOV B,TANG Z,et al.Deterministically subsampled acquisition geometries for optimal reconstruction[J].Expanded Abstracts of 78thEAGE Annual Conference,2016:360-365

[17] CANDES E J,WAKIN M B.An introduction to compressive sampling:a sensing/sampling paradigm that goes against the common knowledge in data acquisition[J].IEEE Signal Processing Magazine,2008,25(2):21-30

[18] HENNENFENT G,HERRMANN F J.Simply denoise:wavefield reconstruction via jittered undersampling[J].Geophysics,2008,73(3):V19-V28

[19] LI C,KAPLAN S T,MOSHER C C,et al.Compressive sensing:US2014/063443[P].2015-05-07

(编辑：朱文杰)

Irregularseismicgeometrydesignanddatareconstructionbasedoncompressivesensing

ZHOU Song1,LV Yao2,LV Gonghe1,SHU Guoxu2,SHI Taikun2,HUO Shoudong2

(1.Sinopec Petroleum Engineering Geophysics Co.,Ltd.,Beijing 100020,China;2.Institute of Geology and Geophysics,Chinese Academy of Sciences,Beijing 100029,China)

This paper proposes a method of irregular seismic geometry optimal design and data reconstruction based on compressive sensing.To optimize the irregular geometry design in seismic data acquisition,the method selects sampling points one by one according to a greedy sequential strategy,which could decrease the maximum cross-correlation value of the sensing matrix.To reconstruct seismic data effectively,the method solves the underdetermined sampling matrix through anL0andL1regularization mixed iterative algorithm.When applied to the irregular seismic geometry design in the TFT exploration area,the method uses the maximum and the minimum group intervals as constraints to overcome the limitation of the receiver cable length.The results show that the proposed method provides an effective way for the development of high-density 3D seismic exploration.

compressive sensing,sparse,irregular,seismic geometry design,seismic data reconstruction

P631

：A

1000-1441(2017)05-061709DOI:10.3969/j.issn.1000-1441.2017.05.001

周松,吕尧,吕公河,等.基于压缩感知的非规则地震勘探观测系统设计与数据重建[J].石油物探,2017,56(5):625

ZHOU Song,LV Yao,LV Gonghe,et al.Irregular seismic geometry design and data reconstruction based on compressive sensing

[J].Geophysical Prospecting for Petroleum,2017,56(5):625

周松(1963—),男,高级工程师,从事地球物理数据采集、非常规勘探以及物探技术管理等工作。

中国石油化工股份有限公司科技攻关项目“基于压缩感知的地震勘探采集技术研究”及中石化石油工程技术服务有限公司科技项目“压缩感知技术在地震勘探中的应用研究”(SG16-52K)联合资助。

This research is financially supported by the project of the Ministry of Science and Technology of SINOPEC and the project of Sinopec Oilfield Service Corporation (Grant No.SG16-52K).

2017-01-24;改回日期：2017-06-15。