CN115657125B - 地震数据重建与去噪的方法、装置和电子设备 - Google Patents
地震数据重建与去噪的方法、装置和电子设备 Download PDFInfo
- Publication number
- CN115657125B CN115657125B CN202211399773.2A CN202211399773A CN115657125B CN 115657125 B CN115657125 B CN 115657125B CN 202211399773 A CN202211399773 A CN 202211399773A CN 115657125 B CN115657125 B CN 115657125B
- Authority
- CN
- China
- Prior art keywords
- seismic data
- frequency
- target
- frequency component
- domain seismic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 86
- 239000011159 matrix material Substances 0.000 claims abstract description 165
- 238000012545 processing Methods 0.000 claims abstract description 85
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 57
- 230000009466 transformation Effects 0.000 claims abstract description 44
- 230000009467 reduction Effects 0.000 claims abstract description 32
- 238000013016 damping Methods 0.000 claims abstract description 18
- 238000010183 spectrum analysis Methods 0.000 claims abstract description 14
- 238000005070 sampling Methods 0.000 claims description 12
- 238000000354 decomposition reaction Methods 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 7
- 238000003860 storage Methods 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 6
- 230000007423 decrease Effects 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000008569 process Effects 0.000 description 11
- 238000010586 diagram Methods 0.000 description 10
- 230000000694 effects Effects 0.000 description 10
- 238000012217 deletion Methods 0.000 description 7
- 230000037430 deletion Effects 0.000 description 7
- 238000004891 communication Methods 0.000 description 5
- 239000010755 BS 2869 Class G Substances 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 2
- 230000000903 blocking effect Effects 0.000 description 2
- 230000006835 compression Effects 0.000 description 2
- 238000007906 compression Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 238000005520 cutting process Methods 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000802 evaporation-induced self-assembly Methods 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000011946 reduction process Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种地震数据重建与去噪的方法、装置和电子设备，涉及地震数据处理的技术领域，包括：获取时间域地震数据；对时间域地震数据进行离散傅里叶变换，得到频率域地震数据；利用阻尼多道奇异谱分析DMSSA算法对频率域地震数据进行处理，得到目标频率域地震数据；DMSSA算法执行矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的；对目标频率域地震数据进行傅里叶反变换，得到数据重建与去噪处理后的时间域地震数据。本发明方法能够自动确定矩阵降秩时所保留的目标奇异值个数，不需要人工预估，缓解了现有的地震数据重建与去噪方法存在的数据处理效率低的技术问题。
Description
技术领域
本发明涉及地震数据处理的技术领域，尤其是涉及一种地震数据重建与去噪的方法、装置和电子设备。
背景技术
地震勘探数据处理的最终目的是获得高信噪比、高保真度、高分辨率的地下结构，而提高地震数据信噪比是获得高保真和高分辨率地震图像的基础。在地震数据的采集过程中，由于环境以及坏道等因素，经常会出现地震道缺失以及随机噪声干扰的问题，降低了地震数据的信噪比。地震数据重建与去噪是获得高信噪比、完备地震数据的主要手段。
现有技术中通常基于矩阵降秩的方法进行地震数据重建和去噪，但目前的降秩方法都是基于地震线性同相轴的基本假设，时间域信号先经过离散快速傅里叶变换后转换成频率域信号，然后通过人工预估的方法对所有的频率截断相应的奇异值个数来实现降秩处理，因此存在数据处理效率低的技术问题。
发明内容
本发明的目的在于提供一种地震数据重建与去噪的方法、装置和电子设备，以缓解了现有的地震数据重建与去噪方法存在的数据处理效率低的技术问题。
第一方面，本发明提供一种地震数据重建与去噪的方法，包括：获取时间域地震数据；对所述时间域地震数据进行离散傅里叶变换，得到频率域地震数据；利用阻尼多道奇异谱分析DMSSA算法对所述频率域地震数据进行处理，得到目标频率域地震数据；其中，所述DMSSA算法执行矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的；对所述目标频率域地震数据进行傅里叶反变换，得到数据重建与去噪处理后的时间域地震数据。
在可选的实施方式中，利用阻尼多道奇异谱分析DMSSA算法对所述频率域地震数据进行处理，包括：步骤一，利用算式计算所述频率域地震数据中频率为wi的频率分量进行第n次迭代后的结果，得到更新后的频率分量；其中，Fn(wi)表示所述更新后的频率分量，n取值1至nmax，nmax表示预设迭代次数；wi∈{w1,wJ}，{w1,wJ}表示所述频率域地震数据的有效频率分量范围；an表示随循环线性减少的标量，a1＝1,anmax＝0，F0(wi)表示所述频率域地震数据中频率为wi的原始频率分量，/>表示预设采样算子，DMSSA(·)表示对当前处理的频率分量进行去噪处理；Fn-1(wi)表示所述频率域地震数据中频率为wi的频率分量进行第n-1次迭代后的结果；步骤二，判断所述更新后的频率分量是否符合预设收敛条件；若符合，则将所述更新后的频率分量作为所述目标频率域地震数据中的频率分量；若不符合，则返回执行上述步骤一至步骤二，直至达到所述预设迭代次数，并将第nmax次迭代后得到的频率分量作为所述目标频率域地震数据中的频率分量。
在可选的实施方式中，对当前处理的频率分量进行去噪处理，包括：对所述当前处理的频率分量进行Toeplitz变换，得到所述当前处理的频率分量对应的块Toeplitz矩阵；对所述块Toeplitz矩阵进行奇异值分解，得到所述当前处理的频率分量对应的酉矩阵和对角矩阵；基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数；基于所述目标奇异值个数、所述当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵；对所述降秩后的块Toeplitz矩阵进行Toeplitz反变换，得到去噪后的频率分量。
在可选的实施方式中，对所述当前处理的频率分量进行Toeplitz变换，包括：对所述当前处理的频率分量中的每一行元素进行Toeplitz变换，得到多个子Toeplitz矩阵；基于所述多个子Toeplitz矩阵构建所述当前处理的频率分量对应的块Toeplitz矩阵。
在可选的实施方式中，基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数；利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，得到所述目标频率分量对应的有效奇异值个数；其中，所述目标频率分量表示所述当前处理的频率域地震数据中所有频率分量中的任一频率分量；基于所有频率分量的对角矩阵处理结果，统计每种有效奇异值个数的数量，并按照数量大小进行降序求和处理，直至参与求和的数量达到预设百分比；将参与求和处理的数量所对应的有效奇异值个数中，最大的有效奇异值个数作为所述目标奇异值个数。
在可选的实施方式中，利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，包括：步骤1.1，基于所述目标频率分量对应的对角矩阵中的奇异值构建目标二维数组；其中，所述目标二维数组中的每个元素由奇异值所属的行数和奇异值本身构成；步骤1.2，将所述目标二维数组中的每个元素作为一类，并计算每两类之间的欧式距离；步骤1.3，合并欧式距离最小的两个类，得到新类；步骤1.4，在类的总数大于2的情况下，计算所述新类与其余各类之间的距离，并返回执行所述步骤1.3；步骤1.5，在类的总数为2的情况下，将元素数量最少的类的元素总数作为目标频率分量对应的有效奇异值个数。
在可选的实施方式中，基于所述目标奇异值个数、所述当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵，包括：截取所述当前处理的频率分量对应的对角矩阵中的前P个奇异值，得到目标对角矩阵；其中，P的取值为所述目标奇异值个数；截取所述当前处理的频率分量对应的酉矩阵中的前P行元素，得到目标酉矩阵；基于所述目标酉矩阵和所述目标对角矩阵确定降秩后的块Toeplitz矩阵。
第二方面，本发明提供一种地震数据重建与去噪的装置，包括：获取模块，用于获取时间域地震数据；第一变换模块，用于对所述时间域地震数据进行离散傅里叶变换，得到频率域地震数据；处理模块，用于利用阻尼多道奇异谱分析DMSSA算法对所述频率域地震数据进行处理，得到目标频率域地震数据；其中，所述DMSSA算法执行矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的；第二变换模块，用于对所述目标频率域地震数据进行傅里叶反变换，得到数据重建与去噪处理后的时间域地震数据。
第三方面，本发明提供一种电子设备，包括存储器、处理器，所述存储器上存储有可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现前述实施方式中任一项所述的地震数据重建与去噪的方法的步骤。
第四方面，本发明提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令被处理器执行时实现前述实施方式中任一项所述的地震数据重建与去噪的方法。
本发明提供的地震数据重建与去噪的方法，能够根据DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵自动确定算法执行矩阵降秩时所保留的目标奇异值个数，不需要人工预估，进而提升了地震数据重建与去噪的处理速度，缓解了现有的地震数据重建与去噪方法存在的数据处理效率低的技术问题。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种地震数据重建与去噪的方法的流程图；
图2a为采用主频为30Hz的雷克子波正演得到一个大小为300×60×60的模拟地震数据；
图2b为对图2a中的数据随机抽取50％地震道且增加10％随机噪声后的数据的示意图；
图3为图2a和图2b中的数据在某频率下Toeplitz矩阵的特征值分布图；
图4为图2b中数据的1-150Hz进行层次聚类法测得的奇异值的数量的示意图；
图5a为采用层次聚类法确定的奇异值P＝3的情况下图2b的重建结果；
图5b为图2a与图5a的差；
图6a为一个大小为300×60×60的三维真实地震数据的示意图；
图6b为对图6a进行随机30％地震道缺失后的地震数据示意图；
图7为本发明实施例提供的一种对图6b中的一个小块通过层次聚类法确定奇异值的结果；
图8a为基于层次聚类法确定奇异值的方法对图6a的重建结果；
图8b为图6a与图8a的差；
图9a为另一种三维真实地震数据的示意图；
图9b为对图9a进行随机30％地震道缺失后的地震数据示意图；
图10a为图9b基于层次聚类法确定奇异值的重建结果；
图10b为图10a与图9a的差；
图11是本发明实施例提供的一种地震数据重建与去噪的装置的功能模块图；
图12是本发明实施例提供的一种电子设备的示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合附图，对本发明的一些实施方式作详细说明。在不冲突的情况下，下述的实施例及实施例中的特征可以相互组合。
现有技术中通常基于矩阵降秩的方法进行地震数据重建和去噪，但目前的降秩方法都是基于地震线性同相轴的基本假设，时间域信号先经过离散快速傅里叶变换后转换成频率域信号，然后通过人工预估的方法对所有的频率截断相同的奇异值个数来实现降秩处理。对于大规模数据来说，基于降秩的方法需要将地震数据分成不同的块，然而每个块对应的奇异值个数不同，目前需要人工对每个块的数据估计合适的奇异值个数，计算效率低，无法实现工业化。
实施例一
图1为本发明实施例提供的一种地震数据重建与去噪的方法的流程图，如图1所示，该方法具体包括如下步骤：
步骤S102，获取时间域地震数据。
在对大规模真实地震数据进行去噪和重建时需要分块处理，通过滑动分块将大块地震数据分成各个小块，然后针对每个小块，再利用本发明实施例提供的方法进行重建与去噪处理。用户可以根据实际需求设置每个小块时间域地震数据的大小，本发明实施例不对其进行具体限定。本步骤获取到的小块的时间域地震数据可理解为是大小为M×N×O的三维含噪地震数据D(x,y,t)，其中，x，y，t分别表示空间x方向，y方向和时间t方向的坐标，x＝(x1,x2,…,xM)，y＝(y1,y2,…,yN)，t＝(t1,t2,…,tO)。
步骤S104，对时间域地震数据进行离散傅里叶变换，得到频率域地震数据。
为了对时间域地震数据进行去噪与重建，本发明实施例需要将地震数据从时间域变换到频率域进行处理，因此，在获取到时间域地震数据D(x,y,t)之后，对时间域地震数据D(x,y,t)进行离散傅里叶变换，得到频率域地震数据F(x,y,w)，其中，w＝(w1,w2,…,wJ)，wJ表示频率域信号的最高频率。
步骤S106，利用阻尼多道奇异谱分析DMSSA算法对频率域地震数据进行处理，得到目标频率域地震数据。
其中，DMSSA算法执行矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的。
通过上文中的描述可知，现有的矩阵降秩方法中，进行频率截断时所使用的奇异值个数是人工预估的，导致数据处理效率低，且无法确保其准确性。为了解决上述问题，本发明实施例在使用DMSSA算法对频率域地震数据进行处理时，矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的，也就是说，本发明实施例能够自动确定可靠的奇异值个数，以提高数据处理效率，以及确保高信噪比的数据重建与去噪效果。经过对频率域地震数据进行处理之后，即可得到目标频率域地震数据FDMSSA(x,y,w)，该目标频率域地震数据FDMSSA(x,y,w)即在频率域对步骤S104所得到的频率域地震数据F(x,y,w)进行数据重建与去噪处理后所得到的结果。
步骤S108，对目标频率域地震数据进行傅里叶反变换，得到数据重建与去噪处理后的时间域地震数据。
也就是说，利用算式对FDMSSA(x,y,w)进行傅里叶反变换处理，即可得到数据重建与去噪处理后的时间域地震数据/>其中，Ψ-1表示傅里叶反变换算子。
本发明实施例所提供的地震数据重建与去噪的方法，能够根据DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵自动确定算法执行矩阵降秩时所保留的目标奇异值个数，不需要人工预估，进而提升了地震数据重建与去噪的处理速度，缓解了现有的地震数据重建与去噪方法存在的数据处理效率低的技术问题。
在一个可选的实施方式中，上述步骤S106，利用阻尼多道奇异谱分析DMSSA算法对频率域地震数据进行处理，具体包括如下步骤：
步骤一，利用算式计算频率域地震数据中频率为wi的频率分量进行第n次迭代后的结果，得到更新后的频率分量。
其中，Fn(wi)表示更新后的频率分量，n取值1至nmax，nmax表示预设迭代次数；wi∈{w1,wJ}，{w1,wJ}表示频率域地震数据的有效频率分量范围；an表示随循环线性减少的标量，a1＝1,anmax＝0，F0(wi)表示频率域地震数据F(x,y,w)中频率为wi的原始频率分量F(x,y,wi)，表示预设采样算子，DMSSA(·)表示对当前处理的频率分量进行去噪处理；Fn-1(wi)表示频率域地震数据中频率为wi的频率分量F(x,y,wi)进行第n-1次迭代后的结果。
步骤二，判断更新后的频率分量是否符合预设收敛条件。
若符合，则将更新后的频率分量作为目标频率域地震数据中的频率分量。
若不符合，则返回执行上述步骤一至步骤二，直至达到预设迭代次数，并将第nmax次迭代后得到的频率分量作为目标频率域地震数据中的频率分量。
缺失的地震道通常在地震数据中以0的形式存在，在对地震数据进行重建时，地震道缺失部分可看作有效信号与随机噪声相加为0的情况，因此，三维地震数据的重建方法与去噪的方法相似，只是在重建空白地震道过程中，采用了类似加权凸集投影的方式。基于上文关于利用阻尼多道奇异谱分析DMSSA算法对频率域地震数据进行处理的具体步骤描述可知，其数据处理流程的伪代码如下：
上述伪代码中，F0(w)＝F(x,y,w)，Fn(w)的算式表示同时进行地震数据重建与去噪的计算过程，表示上述预设收敛条件，ε表示预设收敛阈值。
本发明实施例中，对于频率频率w1到wJ范围内的每个频率w，设定最大迭代次数为nmax，每次迭代执行Fn(w)的计算公式，如果迭代次数未达到nmax，但是满足预设收敛条件则返回Fn(w)为数据结果，将频率范围{w1,wJ}都执行上述运算，就得到频率域的同时重建和去噪结果，再通过傅里叶反变换即可得到时间域的重建和去噪结果。
本发明方法通过循环迭代的方式完成地震数据的重建和去噪，每次迭代中，在地震道缺失位置处由DMSSA方法进行重建，在地震道没有缺失的位置由对数据进行去噪，F0(w)为原始观测数据，Fn-1(w)为第n-1次迭代后得到的数据，an控制两者所占的比率，由1线性递减至0。该方法可以有效降低直接重建时由于原始观测数据的噪声对数据相干性造成的影响，完成地震数据的同时重建与去噪。
已知DMSSA(·)表示对当前处理的频率分量进行去噪处理，在一个可选的实施方式中，DMSSA(·)去噪处理具体包括如下步骤：
步骤S201，对当前处理的频率分量进行Toeplitz变换，得到当前处理的频率分量对应的块Toeplitz矩阵。
基于上文中循环迭代处理的流程可知，DMSSA算法需要对F(x,y,w)(也即F0(w))中的每一个频率分量分别进行迭代处理，以实现去噪与重建，上述当前处理的频率分量可以是任一个频率分量在任一次迭代中的数据。以DMSSA(·)当前处理的频率分量为F0(wi)(也即，F(x,y,wi))为例，表示为：DMSSA(F0(wi))，对F0(wi)进行处理时，首先需要对F0(wi)进行Toeplitz变换，得到其对应的块Toeplitz矩阵，块Toeplitz矩阵表示为：其中，Aq表示对F0(wi)中的第q行元素进行Toeplitz变换所得到的子Toeplitz矩阵，/>符号/>表示下取整运算。块Toeplitz矩阵的行数与列数的取值是根据本领域的常用取值算法确定的，其目的是尽量保持块Toeplitz矩阵接近方阵。
步骤S202，对块Toeplitz矩阵进行奇异值分解，得到当前处理的频率分量对应的酉矩阵和对角矩阵。
具体的，对块Toeplitz矩阵进行奇异值分解，可得到：/>其中，/>表示Lx×Kx阶酉矩阵，/>表示Ly×Ky阶酉矩阵，/>Ky＝N-Kx+1，Σ＝diag(σ1,σ2,…,σs)表示对角矩阵，且σ1＞σ2＞…＞σs。也就是说，F0(wi)对应的酉矩阵为：/>和对角矩阵为：Σ。
步骤S203，基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数。
在对块Toeplitz矩阵进行奇异值分解之后，区别于传统根据线性同相轴的个数确定对角矩阵中所保留的奇异值个数，本发明实施例利用当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数，若当前处理的频率分量为F0(wi)，那么当前处理的频率域地震数据就是F0(w)＝F(x,y,w)。
因此，本步骤需要对F(x,y,w)中的每一个频率分量执行上述步骤S201-S202的处理，得到F(x,y,w)中每一个频率分量对应的对角矩阵，然后再通过对上述所有对角矩阵中的元素值进行处理，以得到目标奇异值个数P，其中，目标奇异值的个数表示对块矩阵进行矩阵降秩处理时所保留的奇异值个数。
步骤S204，基于目标奇异值个数、当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵。
已知那么在确定了目标奇异值个数P之后，如果对F0(wi)对应的对角矩阵Σ中的奇异值进行截断处理，则保留其前P个奇异值，并将其余奇异值设为0，即可得到Σ′＝diag(σ1,σ2,…,σP,0,…,0)，虽然通过/>也可以达到降秩的效果，但是这样处理并不能完全去除随机噪声，因为当地震数据中只有有效信号(也即，没有噪声)时，Σ′中包含与地震图像(也即，D(x,y,t))中线性同相轴数量相同个数的非零奇异值。当地震数据混入随机噪声，所有的奇异值大小都会发生改变，非零奇异值个数将会增加。传统的去噪方法保留前P个奇异值，其余奇异值置为0，但保留的前P个奇异值中依然残留随机噪声引起的变化量，其去噪结果中依然残存着随机噪声。
本发明实施例提出使用阻尼多道奇异谱分析DMSSA去噪方法，将阻尼因子引入去噪方法中，通过阻尼因子减弱随机噪声带来的奇异值的增量，以有效压制残存的随机噪声。具体的，本发明实施例中，可使用算式计算降秩后的块Toeplitz矩阵，其中，/>I表示单位矩阵，D表示阻尼因子，T表示阻尼算子。D的值越小，阻尼效果越强。通过上述算式可以利用阻尼因子D对第P+1个奇异值进行放大或缩小，使用前P个奇异值与其相减，达到压制前P个奇异值中与随机噪声相关部分的目的。
在一个可选的实施方式中，在已知目标奇异值个数的情况下，为了减少矩阵降秩处理的计算量，上述步骤S204，基于目标奇异值个数、当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵，包括：
步骤S2041，截取当前处理的频率分量对应的对角矩阵中的前P个奇异值，得到目标对角矩阵。
其中，P的取值为目标奇异值个数。
步骤S2042，截取当前处理的频率分量对应的酉矩阵中的前P行元素，得到目标酉矩阵。
步骤S2043，基于目标酉矩阵和目标对角矩阵确定降秩后的块Toeplitz矩阵。
具体的，步骤S2041执行结束可得到目标对角矩阵：Σ′＝diag(σ1,σ2,…,σP,0,…,0)，已知F0(wi)对应的酉矩阵为：和/>分别将以上两个酉矩阵进行分解，可得到前P行元素组成的目标矩阵以及剩余元素组成的剩余矩阵，得到/>和/>也即，目标酉矩阵为：/>和/>由于对角矩阵中只有P个有效元素，其余为0，因此，为了避免不必要的计算(与元素0相乘)，本发明实施例将/>的算式进行优化，得到：
步骤S205，对降秩后的块Toeplitz矩阵进行Toeplitz反变换，得到去噪后的频率分量。
在得到降秩后的块Toeplitz矩阵之后，利用算式/>对/>进行Toeplitz反变换，即可得到对F(x,y,wi)去噪后的频率分量/>其中，Γ-1表示Toeplitz矩阵变换的逆。/> 的集合即为目标频率域地震数据FDMSSA(x,y,w)。
在一个可选的实施方式中，上述步骤S201，对当前处理的频率分量进行Toeplitz变换，包括：
步骤S2011，对当前处理的频率分量中的每一行元素进行Toeplitz变换，得到多个子Toeplitz矩阵。
步骤S2012，基于多个子Toeplitz矩阵构建当前处理的频率分量对应的块Toeplitz矩阵。
已知时间域地震数据D(x,y,t)的大小为M×N×O，因此，F(x,y,w)频率为wi的频率分量可表示为对F(x,y,wi)的每一行元素进行Toeplitz变换，即可得到M个子Toeplitz矩阵，其中，对F(x,y,wi)的第q行元素进行Toeplitz变换，即可得到/>其中，q＝1，2，…，M，/>子Toeplitz矩阵的行数与列数的取值是根据本领域的常用取值算法确定的，其目的是尽量保持子Toeplitz矩阵接近方阵。
在得到子Toeplitz矩阵A1至AM之后，利用所有的子Toeplitz矩阵即可构建出F(x,y,wi)对应的块Toeplitz矩阵
在一个可选的实施方式中，上述步骤S203，基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数，包括：
步骤S2031，利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，得到目标频率分量对应的有效奇异值个数。
其中，目标频率分量表示当前处理的频率域地震数据中所有频率分量中的任一频率分量。
随机噪声的增加与地震道的缺失使奇异值中前N项奇异值变小，非零项增多，但前N项奇异值依然与N项之后的奇异值有较大区别。含噪声地震数据的奇异值可以分为两类，一类是与有效信号有关的奇异值，数值较大，另一类是与噪声和道缺失有关的奇异值，数值较小。本发明实施例使用层次聚类的方法对对角矩阵中的元素进行处理，将对角矩阵中的所有元素划分为两类，由于有效奇异值个数远远小于奇异值的总数量，故选取数量较少的一类，统计其包含奇异值的数量即可得到目标频率分量对应的有效奇异值个数N。
步骤S2032，基于所有频率分量的对角矩阵处理结果，统计每种有效奇异值个数的数量，并按照数量大小进行降序求和处理，直至参与求和的数量达到预设百分比。
步骤S2033，将参与求和处理的数量所对应的有效奇异值个数中，最大的有效奇异值个数作为目标奇异值个数。
为了提升本发明实施例所提供方法的稳定性，本发明实施例首先利用步骤S2031所提供的方法计算出所有频率分量对应的有效奇异值个数，例如，得到{(w1,5),(w2,6),(w3,5),…,(wJ,6)}，也即，DMSSA算法当前处理的频率域地震数据Fn(x,y,w)中频率为w1的频率分量F(x,y,w1)对应的有效奇异值个数为5，频率为w2的频率分量F(x,y,w2)对应的有效奇异值个数为6，以此类推。
接下来，统计每一种有效奇异值个数的数量，然后按照按照数量大小进行降序求和处理，并监测参与求和的数量是否达到预设百分比。在确定达到预设百分比之后，再将参与求和处理的数量所对应的有效奇异值个数中，最大的有效奇异值个数作为目标奇异值个数。
为了便于理解，下面举例说明，假设预设百分比为80％，统计出所有频率分量对应的有效奇异值个数共有3种，分别为5，6，7，并且，经层次聚类算法处理后得到有效奇异值个数为5的频率分量共有4个，经层次聚类算法处理后得到有效奇异值个数为6的频率分量共有1个，经层次聚类算法处理后得到有效奇异值个数为7的频率分量共有3个，也即，有效奇异值个数5的数量为4，有效奇异值个数6的数量为1，有效奇异值个数7的数量为3。
接下来，按照数量的大小对有效奇异值个数进行降序求和，也即，首先按照数量大小对所有有效奇异值个数进行降序排序，已知4>3>1，因此有效奇异值个数的排序为{5,7,6}。随后按顺序进行求和处理，首先，数量最多的有效奇异值个数5的占比为：4/(4+1+3)<80％，未达到预设百分比；因此需要将有效奇异值个数为5的数量4与有效奇异值个数为7的数量3进行求和，(4+3)/(4+1+3)>80％，也即，当前参与求和的数量达到了预设百分比，因此，将参与求和处理的数量(分别为4和3)所对应的有效奇异值个数(分别为5和7)中最大的有效奇异值个数7作为目标奇异值个数。
在一个可选的实施方式中，上述步骤S2031，利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，包括：
步骤1.1，基于目标频率分量对应的对角矩阵中的奇异值构建目标二维数组。
其中，目标二维数组中的每个元素由奇异值所属的行数和奇异值本身构成。
若目标频率分量对应的对角矩阵为Σ＝diag(σ1,σ2,…,σs)，那么由该对角矩阵中的奇异值构建的目标二维数组为：((1,σ1)；(2,σ2)；…；(s,σS))，分别用(z1,z2,…,zS)进行表示，也即，zS＝(s,σS)。
步骤1.2，将目标二维数组中的每个元素作为一类，并计算每两类之间的欧式距离。
本发明实施例中，定义类Gp与类Gq之间的距离Dpq为：Dpq＝min{dij|zi∈Gp,zj∈Gq}，其中，zi为类Gp中的二维向量，zj为类Gq中的二维向量，dij表示zi和zj之间的欧式距离。
步骤1.3，合并欧式距离最小的两个类，得到新类。
步骤1.4，在类的总数大于2的情况下，计算新类与其余各类之间的距离，并返回执行步骤1.3。
步骤1.5，在类的总数为2的情况下，将元素数量最少的类的元素总数作为目标频率分量对应的有效奇异值个数。
按照上述步骤，即可计算出所有频率分量对应的有效奇异值个数。在实际处理中，有效奇异值的数量必然大于1，且在含噪数据中第一个奇异值往往远大于其余奇异值。为了保证聚类效果的稳定性，还可以选择在进行聚类处理时可忽略第一个奇异值对应的二维向量z1，直接对(z2,…,zS)进行聚类处理，其聚类后对得到的奇异值数量再加1即可得到目标频率分量对应的有效奇异值个数。
为了验证本发明实施例所提出的方法对实际数据的效果，发明人做了以下实验。图2a是采用主频为30Hz的雷克子波正演得到一个大小为300×60×60的模拟地震数据的示意图，时间采样率为2ms，该数据有3个线性同相轴。图2b为对图2a中的数据随机抽取50％地震道且增加10％随机噪声后的数据的示意图，信噪比为-3.9002dB。模拟数据使用的子波是主频为40HZ的雷克子波，信噪比(SNR)定义为：其中，d表示不含噪声数据，r表示去噪后的数据，||.||2表示L2范数。
图3为图2a和图2b中的数据在某频率下Toeplitz矩阵的特征值分布图，其中实线为图2a数据对应的奇异值，虚线为待重建地震数据(图2b)对应的奇异值。通过图3可看出，图2a地震数据P+1项之后的特征值远小于前P项特征值，且绝大多数为0。随机噪声的增加与地震道的缺失使奇异值中前N项奇异值变小，非零项增多，但前N项奇异值依然与N项之后的奇异值有较大区别。
图4为图2b中数据的1-150Hz进行层次聚类法测得的奇异值的数量的示意图，对于图2b数据来说，10-90Hz在信号的有效频率范围内，层次聚类法确定的奇异值个数稳定在3，而超出这个范围外，测得的奇异值数量不断上下波动，其原因是超出这个频率范围的可能为噪声。根据MSSA的理论，可验证在10-90Hz内测得的奇异值数量是准确的。实际数据情况会更加复杂，为了算法稳定使用有效信号频率范围求取的奇异值个数的最大值来作为最终使用的值。
首先用模拟数据验证方法的有效性。图5a为采用层次聚类法确定的奇异值P＝3的情况下图2b的重建结果，阻尼因子D取2，重建后信噪比为19.8893dB，图5b为残留的噪声，也即图2a与图5a的差。从图5a和图5b可看出，重建后线性同相轴清晰可见，噪声残留较少，说明本发明方法预测的P值准确，能重建出有较高信噪比的地震图像。
为了验证本发明方法对实际数据的处理效果，发明人还使用两个三维叠后地震数据进行验证。图6a是一个大小为300×60×60的三维真实地震数据的示意图，时间上有276个采样点，时间采样率为2ms。主测线方向与联络测线方向上各有40个采样点。图6b为对图6a进行随机30％地震道缺失后的地震数据示意图。该数据噪声能量强，信噪比低，地震道缺失严重。
在计算大规模真实地震数据去噪和重建时需要分块处理的方法，通过滑动分块将大块地震数据分成各个小块，对每个小块通过DMSSA方法进行重建。通过实验发现，空间分块对奇异值个数的计算存在一定的影响，通过多次比较发现空间分块为10×10左右时能取得稳定的效果。
图7为对图6b中的一个小块通过层次聚类法确定奇异值的结果。由图7可得，真实数据相比于模拟数据，奇异值波动更大。在有效频率范围内，奇异值测得的数量相对较为稳定，超出有效频率范围后，测得的奇异值个数出现严重波动，故取此小块在有效频范围内测得的奇异值的个数为3。
图8a为基于层次聚类法确定奇异值的方法对图6a的重建结果，阻尼因子D取2，图8b为去噪后残留的噪声，也即，图6a与图8a的差。图8a的重建与去噪效果良好，有效信号损失较少，剖面的同相轴轮廓恢复清晰，噪声清除干净，证明了基于层次聚类的自适应多道奇异谱分析算法在真实地震数据中的有效性与可行性。
图9a为另一种三维真实地震数据的示意图，大小为276×180×180，时间上有276个采样点，时间采样率为2ms。主测线方向与联络测线方向上各有180个采样点。图9b为对图9a进行随机30％地震道缺失后的地震数据示意图。该数据同样存在噪声能量强，信噪比低，地震道缺失严重的特点。
图10a为图9b基于层次聚类法确定奇异值的重建结果，阻尼因子D取2，图10b为去噪后残留的噪声，也即，图10a与图9a的差。图10a的重建与去噪效果良好，同相轴轮廓恢复清晰，噪声清除干净，有效信号损失较少，进一步证明了基于层次聚类的自适应多道奇异谱分析算法在真实地震数据中的有效性与可行性。
综上所述，本发明实施例提供的方法进行地震数据重建与去噪时，能够自动确定可靠的奇异值个数，不需要人工预估，进而提升了地震数据重建与去噪的处理速度，并且可获得高信噪比的重建效果。数值实验也证明了本发明实施例所提供的方法的有效性和可靠性，在工业化应用中具有巨大潜力。
实施例二
本发明实施例还提供了一种地震数据重建与去噪的装置，该地震数据重建与去噪的装置主要用于执行上述实施例一所提供的地震数据重建与去噪的方法，以下对本发明实施例提供的地震数据重建与去噪的装置做具体介绍。
图11是本发明实施例提供的一种地震数据重建与去噪的装置的功能模块图，如图11所示，该装置主要包括：获取模块10，第一变换模块20，处理模块30，第二变换模块40，其中：
获取模块10，用于获取时间域地震数据。
第一变换模块20，用于对时间域地震数据进行离散傅里叶变换，得到频率域地震数据。
处理模块30，用于利用阻尼多道奇异谱分析DMSSA算法对频率域地震数据进行处理，得到目标频率域地震数据；其中，DMSSA算法执行矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的。
第二变换模块40，用于对目标频率域地震数据进行傅里叶反变换，得到数据重建与去噪处理后的时间域地震数据。
本发明实施例提供的地震数据重建与去噪的装置，能够根据DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵自动确定算法执行矩阵降秩时所保留的目标奇异值个数，不需要人工预估，进而提升了地震数据重建与去噪的处理速度，缓解了现有的地震数据重建与去噪方法存在的数据处理效率低的技术问题。
可选地，处理模块30包括：
计算单元，用于执行步骤一，利用算式计算频率域地震数据中频率为wi的频率分量进行第n次迭代后的结果，得到更新后的频率分量。
其中，Fn(wi)表示更新后的频率分量，n取值1至nmax，nmax表示预设迭代次数；wi∈{w1,wJ}，{w1,wJ}表示频率域地震数据的有效频率分量范围；an表示随循环线性减少的标量，a1＝1,anmax＝0，F0(wi)表示频率域地震数据中频率为wi的原始频率分量，表示预设采样算子，DMSSA(·)表示对当前处理的频率分量进行去噪处理；Fn-1(wi)表示频率域地震数据中频率为wi的频率分量进行第n-1次迭代后的结果。
判断单元，用于执行步骤二，判断更新后的频率分量是否符合预设收敛条件。
确定单元，用于在更新后的频率分量符合预设收敛条件的情况下，将更新后的频率分量作为目标频率域地震数据中的频率分量。
返回单元，用于在更新后的频率分量不符合预设收敛条件的情况下，返回执行上述步骤一至步骤二，直至达到预设迭代次数，并将第nmax次迭代后得到的频率分量作为目标频率域地震数据中的频率分量。
可选地，计算单元包括：
变换子单元，用于对当前处理的频率分量进行Toeplitz变换，得到当前处理的频率分量对应的块Toeplitz矩阵。
分解单元，用于对块Toeplitz矩阵进行奇异值分解，得到当前处理的频率分量对应的酉矩阵和对角矩阵。
第一确定子单元，用于基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数。
第二确定子单元，用于基于目标奇异值个数、当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵。
Toeplitz反变换单元，用于对降秩后的块Toeplitz矩阵进行Toeplitz反变换，得到去噪后的频率分量。
可选地，变换子单元具体用于：
对当前处理的频率分量中的每一行元素进行Toeplitz变换，得到多个子Toeplitz矩阵。
基于多个子Toeplitz矩阵构建当前处理的频率分量对应的块Toeplitz矩阵。
可选地，第一确定子单元包括：
处理子模块，用于利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，得到目标频率分量对应的有效奇异值个数；其中，目标频率分量表示当前处理的频率域地震数据中所有频率分量中的任一频率分量。
统计子模块，用于基于所有频率分量的对角矩阵处理结果，统计每种有效奇异值个数的数量，并按照数量大小进行降序求和处理，直至参与求和的数量达到预设百分比。
确定子模块，用于将参与求和处理的数量所对应的有效奇异值个数中，最大的有效奇异值个数作为目标奇异值个数。
可选地，处理子模块具体用于执行下述步骤：
步骤1.1，基于目标频率分量对应的对角矩阵中的奇异值构建目标二维数组；其中，目标二维数组中的每个元素由奇异值所属的行数和奇异值本身构成。
步骤1.2，将目标二维数组中的每个元素作为一类，并计算每两类之间的欧式距离。
步骤1.3，合并欧式距离最小的两个类，得到新类。
步骤1.4，在类的总数大于2的情况下，计算新类与其余各类之间的距离，并返回执行步骤1.3。
步骤1.5，在类的总数为2的情况下，将元素数量最少的类的元素总数作为目标频率分量对应的有效奇异值个数。
可选地，第二确定子单元具体用于：
截取当前处理的频率分量对应的对角矩阵中的前P个奇异值，得到目标对角矩阵；其中，P的取值为目标奇异值个数。
截取当前处理的频率分量对应的酉矩阵中的前P行元素，得到目标酉矩阵。
基于目标酉矩阵和目标对角矩阵确定降秩后的块Toeplitz矩阵。
实施例三
参见图12，本发明实施例提供了一种电子设备，该电子设备包括：处理器60，存储器61，总线62和通信接口63，所述处理器60、通信接口63和存储器61通过总线62连接；处理器60用于执行存储器61中存储的可执行模块，例如计算机程序。
其中，存储器61可能包含高速随机存取存储器(RAM，Random Access Memory)，也可能还包括非易失性存储器(non-volatile memory)，例如至少一个磁盘存储器。通过至少一个通信接口63(可以是有线或者无线)实现该系统网元与至少一个其他网元之间的通信连接，可以使用互联网，广域网，本地网，城域网等。
总线62可以是ISA总线、PCI总线或EISA总线等。所述总线可以分为地址总线、数据总线、控制总线等。为便于表示，图12中仅用一个双向箭头表示，但并不表示仅有一根总线或一种类型的总线。
其中，存储器61用于存储程序，所述处理器60在接收到执行指令后，执行所述程序，前述本发明实施例任一实施例揭示的过程定义的装置所执行的方法可以应用于处理器60中，或者由处理器60实现。
处理器60可能是一种集成电路芯片，具有信号的处理能力。在实现过程中，上述方法的各步骤可以通过处理器60中的硬件的集成逻辑电路或者软件形式的指令完成。上述的处理器60可以是通用处理器，包括中央处理器(Central Processing Unit，简称CPU)、网络处理器(Network Processor，简称NP)等；还可以是数字信号处理器(Digital SignalProcessing，简称DSP)、专用集成电路(Application Specific Integrated Circuit，简称ASIC)、现成可编程门阵列(Field-Programmable Gate Array，简称FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。结合本发明实施例所公开的方法的步骤可以直接体现为硬件译码处理器执行完成，或者用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器，闪存、只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器61，处理器60读取存储器61中的信息，结合其硬件完成上述方法的步骤。
本发明实施例所提供的一种地震数据重建与去噪的方法、装置和电子设备的计算机程序产品，包括存储了处理器可执行的非易失的程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个处理器可执行的非易失的计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
此外，术语“水平”、“竖直”、“悬垂”等术语并不表示要求部件绝对水平或悬垂，而是可以稍微倾斜。如“水平”仅仅是指其方向相对“竖直”而言更加水平，并不是表示该结构一定要完全水平，而是可以稍微倾斜。
在本发明的描述中，还需要说明的是，除非另有明确的规定和限定，术语“设置”、“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
Claims (7)
1.一种地震数据重建与去噪的方法，其特征在于，包括：
获取时间域地震数据；
对所述时间域地震数据进行离散傅里叶变换，得到频率域地震数据；
利用阻尼多道奇异谱分析DMSSA算法对所述频率域地震数据进行处理，得到目标频率域地震数据；其中，所述DMSSA算法执行矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的；
对所述目标频率域地震数据进行傅里叶反变换，得到数据重建与去噪处理后的时间域地震数据；
其中，利用阻尼多道奇异谱分析DMSSA算法对所述频率域地震数据进行处理，包括：
步骤一，利用算式计算所述频率域地震数据中频率为wi的频率分量进行第n次迭代后的结果，得到更新后的频率分量；
其中，Fn(wi)表示所述更新后的频率分量，n取值1至nmax，nmax表示预设迭代次数；wi∈{w1，wJ}，{w1，wJ}表示所述频率域地震数据的有效频率分量范围；an表示随循环线性减少的标量，a1＝1，anmax＝0，F0(wi)表示所述频率域地震数据中频率为wi的原始频率分量，表示预设采样算子，DMSSA(·)表示对当前处理的频率分量进行去噪处理；Fn-1(wi)表示所述频率域地震数据中频率为wi的频率分量进行第n-1次迭代后的结果，I表示单位矩阵；
步骤二，判断所述更新后的频率分量是否符合预设收敛条件；
若符合，则将所述更新后的频率分量作为所述目标频率域地震数据中的频率分量；
若不符合，则返回执行上述步骤一至步骤二，直至达到所述预设迭代次数，并将第nmax次迭代后得到的频率分量作为所述目标频率域地震数据中的频率分量；
其中，对当前处理的频率分量进行去噪处理，包括：
对所述当前处理的频率分量进行Toeplitz变换，得到所述当前处理的频率分量对应的块Toeplitz矩阵；
对所述块Toeplitz矩阵进行奇异值分解，得到所述当前处理的频率分量对应的酉矩阵和对角矩阵；
基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数；
基于所述目标奇异值个数、所述当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵；
对所述降秩后的块Toeplitz矩阵进行Toeplitz反变换，得到去噪后的频率分量；
其中，基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数，包括：
利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，得到所述目标频率分量对应的有效奇异值个数；其中，所述目标频率分量表示所述当前处理的频率域地震数据中所有频率分量中的任一频率分量；
基于所有频率分量的对角矩阵处理结果，统计每种有效奇异值个数的数量，并按照数量大小进行降序求和处理，直至参与求和的数量达到预设百分比；
将参与求和处理的数量所对应的有效奇异值个数中，最大的有效奇异值个数作为所述目标奇异值个数。
2.根据权利要求1所述的方法，其特征在于，对所述当前处理的频率分量进行Toeplitz变换，包括：
对所述当前处理的频率分量中的每一行元素进行Toeplitz变换，得到多个子Toeplitz矩阵；
基于所述多个子Toeplitz矩阵构建所述当前处理的频率分量对应的块Toeplitz矩阵。
3.根据权利要求1所述的方法，其特征在于，利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，包括：
步骤1.1，基于所述目标频率分量对应的对角矩阵中的奇异值构建目标二维数组；其中，所述目标二维数组中的每个元素由奇异值所属的行数和奇异值本身构成；
步骤1.2，将所述目标二维数组中的每个元素作为一类，并计算每两类之间的欧式距离；
步骤1.3，合并欧式距离最小的两个类，得到新类；
步骤1.4，在类的总数大于2的情况下，计算所述新类与其余各类之间的距离，并返回执行所述步骤1.3；
步骤1.5，在类的总数为2的情况下，将元素数量最少的类的元素总数作为目标频率分量对应的有效奇异值个数。
4.根据权利要求1所述的方法，其特征在于，基于所述目标奇异值个数、所述当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵，包括：
截取所述当前处理的频率分量对应的对角矩阵中的前P个奇异值，得到目标对角矩阵；其中，P的取值为所述目标奇异值个数；
截取所述当前处理的频率分量对应的酉矩阵中的前P行元素，得到目标酉矩阵；
基于所述目标酉矩阵和所述目标对角矩阵确定降秩后的块Toeplitz矩阵。
5.一种地震数据重建与去噪的装置，其特征在于，包括：
获取模块，用于获取时间域地震数据；
第一变换模块，用于对所述时间域地震数据进行离散傅里叶变换，得到频率域地震数据；
处理模块，用于利用阻尼多道奇异谱分析DMSSA算法对所述频率域地震数据进行处理，得到目标频率域地震数据；其中，所述DMSSA算法执行矩阵降秩时所保留的目标奇异值个数是基于算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定的；
第二变换模块，用于对所述目标频率域地震数据进行傅里叶反变换，得到数据重建与去噪处理后的时间域地震数据；
其中，所述处理模块包括：
计算单元，用于执行步骤一，利用算式计算所述频率域地震数据中频率为wi的频率分量进行第n次迭代后的结果，得到更新后的频率分量；
其中，Fn(wi)表示所述更新后的频率分量，n取值1至nmax，nmax表示预设迭代次数；wi∈{w1,wJ}，{w1，wJ}表示所述频率域地震数据的有效频率分量范围；an表示随循环线性减少的标量，a1＝1，anmax＝0，F0(wi)表示所述频率域地震数据中频率为wi的原始频率分量，表示预设采样算子，DMSSA(·)表示对当前处理的频率分量进行去噪处理；Fn-1(wi)表示所述频率域地震数据中频率为wi的频率分量进行第n-1次迭代后的结果，I表示单位矩阵；
判断单元，用于执行步骤二，判断所述更新后的频率分量是否符合预设收敛条件；
确定单元，用于在所述更新后的频率分量符合所述预设收敛条件的情况下，将所述更新后的频率分量作为所述目标频率域地震数据中的频率分量；
返回单元，用于在所述更新后的频率分量不符合所述预设收敛条件的情况下，返回执行上述步骤一至步骤二，直至达到所述预设迭代次数，并将第nmax次迭代后得到的频率分量作为所述目标频率域地震数据中的频率分量；
其中，所述计算单元包括：
变换子单元，用于对所述当前处理的频率分量进行Toeplitz变换，得到所述当前处理的频率分量对应的块Toeplitz矩阵；
分解单元，用于对所述块Toeplitz矩阵进行奇异值分解，得到所述当前处理的频率分量对应的酉矩阵和对角矩阵；
第一确定子单元，用于基于DMSSA算法当前处理的频率域地震数据中所有频率分量对应的对角矩阵确定目标奇异值个数；
第二确定子单元，用于基于所述目标奇异值个数、所述当前处理的频率分量对应的酉矩阵和对角矩阵确定降秩后的块Toeplitz矩阵；
Toeplitz反变换单元，用于对所述降秩后的块Toeplitz矩阵进行Toeplitz反变换，得到去噪后的频率分量；
其中，所述第一确定子单元包括：
处理子模块，用于利用层次聚类算法对目标频率分量对应的对角矩阵进行处理，得到所述目标频率分量对应的有效奇异值个数；其中，所述目标频率分量表示所述当前处理的频率域地震数据中所有频率分量中的任一频率分量；
统计子模块，用于基于所有频率分量的对角矩阵处理结果，统计每种有效奇异值个数的数量，并按照数量大小进行降序求和处理，直至参与求和的数量达到预设百分比；
确定子模块，用于将参与求和处理的数量所对应的有效奇异值个数中，最大的有效奇异值个数作为所述目标奇异值个数。
6.一种电子设备，包括存储器、处理器，所述存储器上存储有可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现上述权利要求1至4中任一项所述的地震数据重建与去噪的方法的步骤。
7.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机指令，所述计算机指令被处理器执行时实现上述权利要求1至4中任一项所述的地震数据重建与去噪的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211399773.2A CN115657125B (zh) | 2022-11-09 | 2022-11-09 | 地震数据重建与去噪的方法、装置和电子设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211399773.2A CN115657125B (zh) | 2022-11-09 | 2022-11-09 | 地震数据重建与去噪的方法、装置和电子设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115657125A CN115657125A (zh) | 2023-01-31 |
| CN115657125B true CN115657125B (zh) | 2024-03-26 |
Family
ID=85016557
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211399773.2A Active CN115657125B (zh) | 2022-11-09 | 2022-11-09 | 地震数据重建与去噪的方法、装置和电子设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115657125B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106646612A (zh) * | 2016-12-20 | 2017-05-10 | 中国地质大学（北京） | 基于矩阵降秩的地震数据重建方法 |
| CN112764103A (zh) * | 2020-07-09 | 2021-05-07 | 五季数据科技(北京)有限公司 | 一种基于稀疏编码特征dbscan聚类地震相分析方法 |
| CN114200525A (zh) * | 2021-12-10 | 2022-03-18 | 河北地质大学 | 一种自适应的多道奇异谱分析地震数据去噪方法 |
| CN115061203A (zh) * | 2022-05-25 | 2022-09-16 | 华北科技学院 | 一种基于频域奇异值分解的矿山单通道微震信号降噪方法及应用 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3593173A1 (en) * | 2017-03-08 | 2020-01-15 | Saudi Arabian Oil Company | Automated system and methods for adaptive robust denoising of large-scale seismic data sets |
| CN109164483B (zh) * | 2018-08-29 | 2020-04-03 | 中国科学院地球化学研究所 | 多分量地震数据矢量去噪方法及多分量地震数据矢量去噪装置 |
| US20200217979A1 (en) * | 2019-01-08 | 2020-07-09 | King Fahd University Of Petroleum And Minerals | Observation-driven method based on iir wiener filter for microseismic data denoising |
-
2022
- 2022-11-09 CN CN202211399773.2A patent/CN115657125B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106646612A (zh) * | 2016-12-20 | 2017-05-10 | 中国地质大学（北京） | 基于矩阵降秩的地震数据重建方法 |
| CN112764103A (zh) * | 2020-07-09 | 2021-05-07 | 五季数据科技(北京)有限公司 | 一种基于稀疏编码特征dbscan聚类地震相分析方法 |
| CN114200525A (zh) * | 2021-12-10 | 2022-03-18 | 河北地质大学 | 一种自适应的多道奇异谱分析地震数据去噪方法 |
| CN115061203A (zh) * | 2022-05-25 | 2022-09-16 | 华北科技学院 | 一种基于频域奇异值分解的矿山单通道微震信号降噪方法及应用 |
Non-Patent Citations (5)
| Title |
|---|
| An open-source Matlab code package for improved rank-reduction 3D seismic data denoising and reconstruction;Yangkang Chen等;Computers & Geosciences;第2页 * |
| 基于收缩多道奇异谱分析的含噪地震信号重建;杨志鹏等;四川地质学报;第485-487页 * |
| 基于改进阻尼多通道奇异谱分析的地震数据去噪研究;朱跃飞等;中国石油学会 2021 年物探技术研讨会论文集;第459-461页 * |
| 基于阻尼多道奇异谱分析算法的二维不规则缺失地震道重建;李崇明等;中国地球科学联合学术年会 2018;全文 * |
| 朱跃飞等.基于改进阻尼多通道奇异谱分析的地震数据去噪研究.中国石油学会 2021 年物探技术研讨会论文集.2021,第459-461页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115657125A (zh) | 2023-01-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Li et al. | Image classification by a two-dimensional hidden Markov model | |
| CN105843919A (zh) | 一种基于多特征融合与聚类集成的移动对象轨迹聚类方法 | |
| CN110096956B (zh) | 基于eemd和排列熵二阶差分的信号去噪方法及装置 | |
| CN110909480B (zh) | 一种水轮机振动信号的去噪方法与装置 | |
| CN112364291B (zh) | 一种前置滤波极值点优化集合经验模态分解方法及装置 | |
| CN111105068B (zh) | 基于序回归学习的数值模式订正的方法 | |
| CN103473755B (zh) | 基于变化检测的sar图像稀疏去噪方法 | |
| CN113537102B (zh) | 一种微震信号的特征提取方法 | |
| CN112990139A (zh) | 基于变模态分解加权重构信号结合小波阈值去噪方法 | |
| CN112765550A (zh) | 一种基于Wi-Fi信道状态信息的目标行为分割方法 | |
| CN115657125B (zh) | 地震数据重建与去噪的方法、装置和电子设备 | |
| CN112285793B (zh) | 一种大地电磁去噪方法及系统 | |
| CN102903104B (zh) | 一种基于减法聚类的快速图像分割方法 | |
| CN117335841A (zh) | 一种用于光伏组件间载波通信信号去噪的方法及系统 | |
| CN110248325B (zh) | 一种基于信号多重消噪的蓝牙室内定位系统 | |
| CN114118177B (zh) | 一种基于奇异谱分析的位平面降噪方法、系统及存储介质 | |
| CN104320659A (zh) | 背景建模方法、装置及设备 | |
| Pizurica et al. | Image de-noising in the wavelet domain using prior spatial constraints | |
| CN114898767A (zh) | 基于U-Net的机载语音噪音分离方法、设备及介质 | |
| Chang et al. | Applying prior correlations for ensemble-based spatial localization | |
| CN107607993B (zh) | 一种确定叠加速度的方法、装置及系统 | |
| CN118157696B (zh) | 一种车载5g天线阵列信号计量方法及系统 | |
| CN111008356A (zh) | 一种基于WTSVD算法扣除背景的γ能谱集分析方法 | |
| CN117454095B (zh) | 一种桥梁动挠度数据降噪方法 | |
| CN116500553B (zh) | 雷达互相干扰抑制方法、装置、计算机设备及存储介质 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |