CN104345340B - 一种地震数据规则化中抖动稀疏采样方法 - Google Patents

Abstract

本发明提供了一种地震数据规则化中抖动稀疏采样方法，属于地震勘探数据规则化领域。本方法包括：(1)获取地震信号的稀疏表达式；(2)地震道进行分组重采样，形成满足抖动Bernoulli分布规律的欠采样地震数据；(3)应用重构算法对步骤(2)得到的欠采样地震数据进行重构，生成完整地震数据。利用本发明方法能通过较少地震道用压缩感知方法获取完整地震道重构结果，为不规则地震数据提供一种地震数据规则化方法。

Description

一种地震数据规则化中抖动稀疏采样方法

技术领域

本发明属于地震勘探数据规则化领域，具体涉及一种地震数据规则化中抖动稀疏采样方法。

背景技术

地震数据规则化是地震信号处理的一个重要步骤，也是近年地震资料处理的一个研究热点。近年来受到广泛关注的压缩感知技术已经被应用到地震数据规则化中，可以有效地规则化不规则数据体。

压缩感知(Compressed sensing，简称CS)理论与技术是2006年美国科学家基于信号稀疏性或可压缩性提出的一种介于数学和信息科学之间的全新信号采集、编解码理论。该理论突破了传统信号采集必须满足Shannon/Nyquist采样定理的限制，可以用远少于传统采样定理所需的采样点数或观测点数高质量重构出原信号或图像。

压缩感知理论包含三个关键技术：信号的稀疏表达、随机测量矩阵的设计、高效恢复算法的设计与实现。

将压缩感知理论与技术应用到地震数据规则化需要解决三个问题：寻找适合地震信号的稀疏变换、如何丢道才能使丢失的地震数据被高质量地恢复出来、寻找适合地震信号的高效率恢复算法。本专利提出的方法是应用压缩感知技术实现地震数据规则化需要解决的三个关键技术之一，即如何缺道才能使丢失数据被高质量重构出来。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种地震数据规则化中抖动稀疏采样方法，对应用该方法得到的满足抖动Bernoulli分布规律的随机稀疏采样的不完整地震道，应用压缩感知方法获取完整地震道重构结果。

本发明是通过以下技术方案实现的：

一种地震数据规则化中抖动稀疏采样方法，所述方法包括：

(1)获取地震信号的稀疏表达式；

(2)对地震道进行分组重采样，形成满足抖动Bernoulli分布规律的欠采样地震数据；

(3)应用重构算法对步骤(2)得到的欠采样地震数据进行重构，生成完整地震数据。

所述步骤(1)是这样实现的：

选取Curvelet变换作为地震信号的稀疏表达方式，具体如下：

将长度为N的原始地震信号表示为x＝[x₁,...,x_N]^T(原始单一信号离散数据)，Ψ＝[ψ₁...ψ_N]是原始地震信号的稀疏表达的N×N的一组正交基矩阵，则将原始地震信号表示成：

x＝Ψs (1-1)

式中，s＝[s₁,...,s_N]为变换系数；

用x和M×N的测量矩阵Φ计算出M<N的y：

y＝Φx＝ΦΨs＝Θs (1-2)

其中Θ＝ΦΨ是M×N矩阵；

将一个规则的2D地震炮记录用矩阵X表示为(x是单一个地震道数据，X是多个地震道数据组成炮记录)：

式中表示第i道地震信号，N_p是采样点数，N_t表示原始采集地震道数。

所述步骤(2)是这样实现的：

将用公式(1-3)表示的2D地震炮记录按照缺道比例η，每P道分成一组，应用公式(1-4)和公式(1-5)确定欠采样地震道的位置：

R是N_s×N_t的约束矩阵，R中第i行，仅第j列为1，其余为零，表示仅第j列被采样；当第j列数值全部为零时，表示第j列地震道丢失了；在这种情况下，将方程(1-2)中的测量矩阵表示成：

(1-5)式表示两个矩阵的Kronecker积，I为N_p×N_p的单位矩阵；将方程(1-3)中X重新排列成一个列矢量该列矢量即为欠采样地震数据。

所述步骤(3)是这样实现的：

将方程(1-2)重写为：

其中y是N_s×N_p的长度矢量，N_s是采样地震道，在方程(1-6)中应用SPGL1重构之后，依据方程(1-1)重构方程(1-3)中规则地震记录X，即用公式(1-1)计算单个地震道的数据，然后用公式(1-3)形成多地震道组成的完整地震数据。

与现有技术相比，本发明的有益效果是：

①提供一种基于压缩感知技术的地震道抖动稀疏选取方案。

②得到满足抖动Bernoulli分布规律的随机稀疏采样的不完整地震道。

③由较少地震道用压缩感知方法获取完整地震道重构结果。

④为不规则地震数据提供一种地震数据规则化方法。

附图说明

图1(a)是2D地震数据规则化中的重构丢失地震道。

图1(b)是2D地震数据规则化中的基于CS的2D地震数据规则化图解。

图2是基于Bernoulli分布的地震道抖动稀疏采样方式示意图(采集20％数据)。

图3是本发明方法的步骤框图。

图4a是效果评估中的原始炮记录。

图4b-1是效果评估中的满足抖动Bernoulli随机分布规律的采集1/2道观测系统。

图4b-2是效果评估中的采集1/2道炮记录。

图4c是效果评估中的由1/2道数据应用曲波变换和SPGL1算法重构的完整炮记录。

图4d是效果评估中的重构误差(SNR为61.52)。

具体实施方式

下面结合附图对本发明作进一步详细描述：

在地震数据采集过程中，地震数据在时间上的采样通常是规则的，而在空间上由于受到自然采集条件(如河流、道路、村庄、建筑等)的限制或为了降低勘探成本，沿检波器方向或震源方向采集的数据体不规则或不满足Shannon-Nyquist采样定理。此外，也可能原始采集数据本身就不完整和不规则，缺少地震道/炮。而大量地震信号处理算法(如多次波消除，波动方程偏移，NMO，AVO等)要求勘探区域满足密集且规则的覆盖次数，以保证地震信号的横向分辨率。为了满足这一需求，有必要在处理前期进行地震数据规则化(或插值)处理，将空间不规则的地震数据通过室内处理形成规则数据体。

压缩感知技术突破了传统的Shannon-Nyqiust采样定理的限制，用采集的少量地震数据重构完整数据，可以用于解决地震数据规则化问题。

基于压缩感知技术的地震数据规则化的质量极大地取决于地震道随机稀疏采样方式。从已采集地震道数据中二次采样哪些地震道数据或从野外地震采集中只采集哪些地震道数据才能更好地用压缩感知技术重构完整数据，就是本专利的内容，即确定参与重构运算的地震道选取方式问题。

本发明提出了新的地震道随机稀疏采样方法-抖动Bernoulli分布地震道随机稀疏采样方法。该方法是Bernoulli分布地震道随机稀疏采样方法的改进，能较好地控制稀疏采集地震道的分布范围，克服局部位置缺道过于集中的问题，提高数据重构质量。

(1)基本原理

压缩感知理论表明，如果信号在某个变换域是稀疏或可压缩的，就可以利用一个与变换基不相关的随机测量矩阵将变换系数投影到一个低维空间上，再根据这些少量观测值，通过求解优化问题实现信号的高质量恢复。

本发明只涉及一种随机测量矩阵的选取，即欠采样地震道选取方法，也就是从规则地震道中随机丢失哪些地震道或从不规则数据体中选取哪些道参与重构运算。

①压缩感知技术

令长度为N的原始离散时间信号为x＝[x₁,...,x_N]^T(列向量)。Ψ＝[ψ₁...ψ_N]是信号的稀疏表达(如Curvelet变换)的N×N的一组正交基矩阵。则信号可以表示成：

x＝Ψs (1-1)

式中s＝[s₁,...,s_N]为变换系数。如果变换系数中有K个非零值，(N-K)个零值，则称信号x是K稀疏的，并且希望K<<N。对于大多数实际信号，这种理想情况很少出现。作为替代，要求信号是可压缩的，即只有少量的大系数和大量的小系数。

考虑一个线性测量过程，用x和M×N的测量矩阵Φ计算出M<N的y：

y＝Φx＝ΦΨs＝Θs (1-2)

其中Θ＝ΦΨ是M×N矩阵。由于M<N，与按照Shannon-Nyqiust采样定理采样的数据相比，压缩感知技术可以采集更少数据。

测量矩阵Φ需要从长度M<N的y中和长度为N的x中重构出来。值得注意的是，线性方程组(1-2)是欠定的。K稀疏或可压缩信号稳定解的一个充分条件是Θ满足限制等距性质(RIP)或Φ满足不相干性。从理论上，当M≥cK log(N/K)(其中c是一个小常数)时，通过将Φ设为随机矩阵，RIP和不相干性都可以非常高概率地得到满足。

规则的稀疏采样(欠采样)会导致假频，而随机的稀疏采样会将假频转化为随机噪声。对于地震数据规则化而言，地震道的稀疏采样需要满足一定的随机规律。

②基于压缩感知的2D地震数据规则化

设一个规则的2D地震炮记录用矩阵X表示：

式中表示第i道地震信号，N_p是采样点数，N_t表示完整采集地震道数。目前仅采样N_s道，随机丢失(N_t-N_s)道。2D地震数据规则化就是应用CS技术重构那些丢失的地震道，如图1(a)所示。设R是N_s×N_t的约束矩阵：

R中第i行，仅第j列为1，其余为零，表示仅第j列被采样。当第j列数值全部为零时，表示第j列地震道丢失了。在这种情况下，在方程(1-2)中测量矩阵可以表示成：

表示两个矩阵的Kronecker积，I为N_p×N_p的单位矩阵。将方程(1-3)中X重新排列成一个列矢量并将方程(1-2)重写为：

其中y是N_s×N_p的长度矢量，N_s是采样地震道，如图1(b)所示。在方程(1-6)中应用SPGL1重构之后，就能依据方程(1-1)重构方程(1-3)中规则地震记录X。

③抖动Bernoulli随机稀疏采样方案及应用效果

带Bernoulli分布的抖动稀疏采样方式是对Bernoulli分布稀疏采样方式的改进，利用Bernoulli分布对稀疏采样地震道位置产生微小扰动，控制实际采集地震道的分布范围，避免缺失地震道位置过于集中，从而提高缺失数据重构质量。设P是一个稀疏采样因子，它表示将需要采集的完整地震道N_t每P道分成一组，每组有1道被采样，被采样地震道的随机数为1，抖动稀疏采样道数N_s为随机数1的累加。再假设N_t是P的倍数，以确保采样地震道数N_s＝N_t/P是整数，此外稀疏采样率η＝N_s/N_t(用分数或百分数表示)，如图2所示。图中显示满足抖动Bernoulli分布规律的一种地震道稀疏采样方式，其中20个地震道分成4组，每组由5个地震道组成，在每组中采集1道，总计采集总道数的1/5，即只采集4个地震道(采集地震道用实心倒三角表示)，再用采集的这4个地震道数据重构完整20道地震记录。

图4a至图4d为对满足抖动Burnoulli分布规律采集的不完整地震道应用压缩感知技术重构效果分析。图4a为2D原始数值模拟炮记录(200道)，图4b-1图为满足抖动Bernoulli随机分布规律的采集1/2地震道观测系统(黑点为采集地震道，空白为未采集地震道)，图4b-2是相应的按照抖动Bernoulli随机分布规律稀疏采集50％地震道的炮记录(表示采集100道，缺失100道)，图4c为应用Curvelet变换和SPGL1重构算法获取的完整单炮记录(200道)。图4d为重构误差，即原始炮记录(图4a)与重构炮记录(图4c)的残差，重构误差质量用SNR表示，SNR越大重构效果越好，本实例SNR为61.52，重构误差很小，成像精度高。

(2)技术实现步骤

如图3所示，本发明的步骤如下：

①地震信号稀疏表达

针对2D完整地震数据(炮或道集记录)，选取具有多尺度、各向异性特征的Curvelet变换作为地震信号稀疏表达方式(公式1-1)。

②用抖动Bernoulli分布随机稀疏采样方案对地震道进行一定比例分组重采样

设计满足抖动Bernoulli分布随机地震道稀疏采样方案：

将完整地震道根据需要按照一定缺道比例η(如η＝25％，表示只采集完整地震数据的1/4地震道，缺失3/4地震道)每P道(P为稀疏采样因子)分成一组(如η＝25％，表示每4个地震道分成1组，在每组中采集1个地震道)，与①结合，应用公式1-4和公式1-5确定欠采样地震道位置，最后形成图4b-1的欠采样方案。所述缺道比例η要根据实际数据缺道情况和实际处理目的来确定，缺道少重构效果好，比如一般缺道比例采用25％，33％，50％都可以重构，但50％效果更好。有时为了节省成本，可能需要只采集25％地震道数据，再由25％数据重构100％数据。

③形成欠采样地震数据

形成满足抖动Bernoulli分布规律的欠采样地震数据。

④应用重构算法产生重构完整数据

对③形成的欠采样数据，利用公式1-6应用SPGL1重构S之后，再依据公式(1-1)重构公式(1-3)中规则地震记录X。

⑤重构效果评估

将基于压缩感知技术重构的完整数据④与最初完全采样数据进行对比，说明方法的有效性(图3中从第一个方框到最后一个方框的箭头就表示利用最初数据来对比)。

本发明专利针对地震数据规则化问题，设计一种满足抖动Bernoulli分布地震道随机稀疏采集方案。该方案可以控制稀疏采集地震道的分布范围，以便获得最佳重构效果。应用该采集方案对已采集地震数据进行二次采样或野外直接采集少量地震道数据，再经过室内压缩感知技术处理获取完整地震数据，实现地震数据规则化。为不能得到规则地震数据地区提供一种获得完整数据途径，可增加叠加次数，改善地震成像质量。只采集部分数据也可节约采集成本。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (2)

1.一种地震数据规则化中抖动稀疏采样方法，其特征在于：所述方法包括：

(1)获取地震信号的稀疏表达式；

(2)地震道进行分组重采样，形成满足抖动Bernoulli分布规律的欠采样地震数据；

(3)应用重构算法对步骤(2)得到的欠采样地震数据进行重构，生成完整地震数据，

所述步骤(1)是这样实现的：

选取Curvelet变换作为地震信号的稀疏表达方式，具体如下：

将长度为N的原始地震信号表示为x＝[x₁,...,x_N]^T，Ψ＝[ψ₁ ... ψ_N]是原始地震信号的稀疏表达的N×N的一组正交基矩阵，则将原始地震信号表示成：

x＝Ψs (1-1)

式中，s＝[s₁,...,s_N]为变换系数；

用x和M×N的测量矩阵Φ计算出M<N的y：

y＝Φx＝ΦΨs＝Θs (1-2)

其中,y表示对原始地震信号x稀疏后的信号,Θ＝ΦΨ是M×N矩阵；

将一个规则的2D地震炮记录用矩阵X表示为：

$X=\&lsqb;x_1,x_2,\mathrm{...},x_{N_t}\&rsqb;---(1-3)$

式中x_i表示第i道地震信号，N_p是采样点数，表示采样点数N_p的约束矩阵,N_t表示原始采集地震道数，

所述步骤(2)是这样实现的：

将用公式(1-3)表示的2D地震炮记录按照缺道比例η，每P道分成一组，应用公式(1-4)和公式(1-5)确定欠采样地震道的位置：

$R=\left[\begin{array}{ccc}1& & \\ & 1& \\ & & 1\end{array}\right]---(1-4)$

R是N_s×N_t的约束矩阵，N_s是目前采样地震道数，R中第i行，仅第j列为1，其余为零，表示仅第j列被采样；当第j列数值全部为零时，表示第j列地震道丢失了；在这种情况下，将方程(1-2)中的测量矩阵表示成：

$\mathrm{\&Phi;}=R\&CircleTimes;I---(1-5)$

(1-5)式表示两个矩阵的Kronecker积，I为N_p×N_p的单位矩阵；将方程(1-3)中X重新排列成一个列矢量该列矢量即为欠采样地震数据。

2.根据权利要求1所述的地震数据规则化中抖动稀疏采样方法，其特征在于：所述步骤(3)是这样实现的：

将方程(1-2)重写为：

$y=\mathrm{\&Phi;}x=(R\&CircleTimes;I)\mathrm{\&Psi;}s---(1-6)$

其中y是N_s×N_p的长度矢量，N_s是目前采样地震道数，在方程(1-6)中应用SPGL1重构之后，用公式(1-1)计算单个地震道的数据，然后用公式(1-3)形成多地震道组成的完整地震数据。