CN106646612B - 基于矩阵降秩的地震数据重建方法 - Google Patents
基于矩阵降秩的地震数据重建方法 Download PDFInfo
- Publication number
- CN106646612B CN106646612B CN201611182694.0A CN201611182694A CN106646612B CN 106646612 B CN106646612 B CN 106646612B CN 201611182694 A CN201611182694 A CN 201611182694A CN 106646612 B CN106646612 B CN 106646612B
- Authority
- CN
- China
- Prior art keywords
- matrix
- toeplitz
- data
- quadruple
- rank
- 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.)
- Expired - Fee Related
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 426
- 238000000034 method Methods 0.000 title claims abstract description 147
- 230000008602 contraction Effects 0.000 title abstract 3
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 86
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 61
- 238000012935 Averaging Methods 0.000 claims abstract description 52
- 238000012545 processing Methods 0.000 claims abstract description 44
- 230000009467 reduction Effects 0.000 claims description 174
- 239000013598 vector Substances 0.000 claims description 33
- 230000006835 compression Effects 0.000 claims description 26
- 238000007906 compression Methods 0.000 claims description 26
- 230000009466 transformation Effects 0.000 claims description 14
- 238000001514 detection method Methods 0.000 claims description 4
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 40
- 238000013500 data storage Methods 0.000 abstract description 20
- 238000010183 spectrum analysis Methods 0.000 description 20
- 238000005516 engineering process Methods 0.000 description 14
- 230000001788 irregular Effects 0.000 description 12
- 238000010586 diagram Methods 0.000 description 10
- 238000003860 storage Methods 0.000 description 10
- 238000005070 sampling Methods 0.000 description 9
- 230000008569 process Effects 0.000 description 8
- 238000004590 computer program Methods 0.000 description 7
- 230000008901 benefit Effects 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 238000001914 filtration Methods 0.000 description 4
- 230000006870 function Effects 0.000 description 4
- 230000005012 migration Effects 0.000 description 4
- 238000013508 migration Methods 0.000 description 4
- 239000000463 material Substances 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 230000002411 adverse Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- VWDWKYIASSYTQR-UHFFFAOYSA-N sodium nitrate Chemical compound [Na+].[O-][N+]([O-])=O VWDWKYIASSYTQR-UHFFFAOYSA-N 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 239000002699 waste material Substances 0.000 description 2
- 238000012300 Sequence Analysis Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/32—Transforming one recording into another or one representation into another
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种基于矩阵降秩的地震数据重建方法，包括：获取由地震频率切片数据构建的四重块Toeplitz矩阵，并存储所述四重块Toeplitz矩阵的部分行列元素以表示所述四重块Toeplitz矩阵，所述部分行列元素包括所述四重块Toeplitz矩阵的各不同元素；利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理；利用无展开求平均算法对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均值，得到所述地震频率切片数据的重建数据；对所述地震频率切片的重建数据做傅里叶反变换，得到时间域地震重建数据，用于油气勘探。本发明的方法能够减少数据储存量和计算量且不降低重建质量。
Description
技术领域
本发明涉及油气勘探开发技术领域，尤其涉及一种基于矩阵降秩的地震数据重建方法。
背景技术
随着我国油气勘探开发程度的不断提高，地震勘探对象逐渐从构造油气藏转向岩性油气藏，地震勘探目标也变得日益复杂，这对地震数据的处理质量提出了更高要求。然而，在地震数据的采集和处理过程中，经常会遇到不规则缺失道地震数据。它的存在会对地震数据后续多道处理技术的正确运行产生不良影响，降低地震数据的处理品质。
不规则地震数据产生的原因主要来源于两个方面：一是在地震数据的采集阶段，由于存在障碍物、禁采区、海上地震拖缆羽状漂移以及采集的经济成本考虑等因素，使得地震数据在空间方向通常呈现稀疏或不规则分布；二是在地震数据预处理阶段，由于剔除废炮和废道等因素也会引起地震数据在空间方向的不规则分布。不规则地震数据不但会对地震后续处理产生噪声干扰，而且会对地震多道处理技术的正确运行产生不良影响，其中以波动方程叠前深度偏移，自由表面多次波消除，谱估计和时移地震可重复性处理最为严重。对于空间稀疏采样数据，在偏移处理中会产生空间假频，这不但影响地震偏移速度场的精确建立，而且会降低偏移成像质量，使构造模糊，断点不清晰，最终达不到油藏精确解释和描述的目的。对于不规则采样地震数据，会引起地下面元覆盖次数不均匀，在叠加成像时出现采集脚印现象。此外，在时移地震油藏监测中由不规则采样引起的前后两次监测差异也会反映在处理剖面上，但这种差异实质并不代表地下流体的真实变化。总之，以上面临的问题都需要借助地震数据插值重建技术来解决。
此外，在我国3D地震勘探施工方式已经普及，3D地震勘探必然涉及叠前5D不规则缺失道地震数据的重建问题(5D分别指时间t、炮点坐标sx和sy、检波器rx和ry，也可以指时间t、共中心点cmpx和cmpy、偏移距hx和hy，还可以指时间t、共中心点cmpx和cmpy、偏移距h和方位角θ)。
与2D和3D地震数据重建相比，5D数据重建可以利用更多维数中的已知道信息，并将这些信息作为约束条件来对缺失道进行重建，获得更加精确的波场重建结果。因此，叠前5D地震数据重建是不规则缺失道重建领域的重点研究对象。
地震数据重建方法在经历了基于函数变换重建、预测滤波重建和波动方程重建之后，目前已进入第四代基于降秩法的重建，基于矩阵降秩的地震数据重建是降秩法重建中最具代表性的一类方法。矩阵降秩法重建的基本原理为假设原始完整采样的理想无噪声地震数据可以被一个低秩矩阵表示，而随机噪声和不规则缺失采样会增加矩阵的秩。因此，地震数据重建问题可以被看做一个高秩矩阵的降秩问题。目前矩阵降秩的方法主要有Cadzow滤波法(Cadzow filtering，简称CF)(Rank-reduction-based trace interpolation,SEGTechnical Program Expanded Abstracts 2010,3829-3833)、基于随机SVD分解的多道奇异谱分析方法(Multichannel singular spectrum analysis，简称MSSA)(Simultaneousseismic data denoising and reconstruction via multichannel singular spectrumanalysis,Geophysics,2011,76(3):V25-V32)以及基于Lanczos双对角分解的多道奇异谱分析方法(A fast reduced-rank interpolation method for prestack seismicvolumes that depend on four spatial dimensions,Geophysics,2013,78(1):V21–V30)。Cadzow滤波方法与多道奇异谱分析方法在方法原理上是一致的，区别在于方法来源不同。Cadzow滤波方法来源于语音信号增强领域，多道奇异谱分析方法来源于古气候时间序列分析领域。这两种方法在进行数据重建时需要在频率域对每个频率切片数据构建Hankel矩阵或者Toeplitz矩阵，然后运用截断SVD分解算法，保留前几个较大奇异值，舍弃较小的奇异值来实现矩阵降秩。由于截断SVD分解计算速度慢，计算成本高，Oropeza和Sacchi提出运用随机SVD分解算法取代截断SVD分解来抬高降秩分解计算效率，对于2D和3D等低维数据重建而言这种方法有一定效果(Simultaneous seismic data denoising andreconstruction via multichannel singular spectrum analysis,Geophysics,2011,76(3):V25-V32；基于多道奇异谱分析的三维地震数据规则化方法,石油地球物理勘探，2014，49(5)：846-851)。但是，对于高维地震数据重建尤其是叠前5D地震数据重建而言，需要构建大型四重块Hankel矩阵或者四重块Toeplitz矩阵，随机SVD分解算法因计算量大，计算成本高无法用于该大型四重块Hankel矩阵或者四重块Toeplitz矩阵的降秩分解。Gao等人(Afast reduced-rank interpolation method for prestack seismic volumes thatdepend on four spatial dimensions,Geophysics,2013,78(1):V21–V30)提出运用Lanczos双对角分解取代随机SVD分解实现对四重块Hankel或者四重块Toeplitz矩阵的降秩，虽然运用Lanczos双对角分解方法可以提高降秩计算效率，但是在对降秩后的四重Toeplitz矩阵沿对角线求平均时仍需先生成该大型低秩四重块Toeplitz矩阵，导致数据的存储量大，重建效率仍很低。
因此，5D地震数据重建除面临计算量大的问题之外，还面临着数据存储量大的挑战。如何寻找一种计算效率高，降秩计算速度快的降秩重建方法且能减小或者压缩数据存储量也是5D重建面临的迫切问题。
发明内容
本发明提供一种基于矩阵降秩的地震数据重建方法，以提高地震数据重建速度且减小数据存储量。
本发明提供一种基于矩阵降秩的地震数据重建方法，包括：获取地震频率切片数据的四重块Toeplitz矩阵，并存储所述四重块Toeplitz矩阵的部分行列元素以表示所述四重块Toeplitz矩阵，所述部分行列元素包括所述四重块Toeplitz矩阵的各不同元素；利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理；利用无展开求平均算法对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均，得到所述地震频率切片数据的重建数据；对所述地震频率切片数据的重建数据做傅里叶反变换，得到时间域地震重建数据，用于油气勘探。
一个实施例中，获取地震频率切片数据的四重块Toeplitz矩阵，包括：对时间域地震数据做关于时间的傅里叶变换，得到频率域地震数据，并从所述频率域地震数据中提取得到所述地震频率切片数据；利用所述地震频率切片数据构建所述四重块Toeplitz矩阵，或利用所述地震频率切片数据构建四重块Hankel矩阵，并将所述四重块Hankel矩阵变换为所述四重块Toeplitz矩阵。
一个实施例中，对时间域地震数据做关于时间的傅里叶变换，得到频率域地震数据之前，还包括：通过坐标变换将地震勘探数据由炮-检域转换到共中心点-偏移距域，得到所述时间域地震数据。
一个实施例中，所述时间域地震数据为5D地震数据。
一个实施例中，所述部分行列元素为所述四重块Toeplitz矩阵的第一行元素及第一列元素或最后一行元素及最后一列元素。
一个实施例中，利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理，包括：通过将表示的四重块Toeplitz矩阵乘以一随机矩阵对表示的四重块Toeplitz矩阵进行矩阵压缩，得到压缩矩阵，所述随机矩阵的行数和列数分别为所述四重块Toeplitz矩阵的列数和秩；利用随机QR降秩分解算法对所述压缩矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理。
一个实施例中，通过将表示的四重块Toeplitz矩阵乘以一随机矩阵对表示的四重块Toeplitz矩阵进行矩阵压缩，包括：用向量形式表示所述随机矩阵，并利用基于四维傅里叶快速变换的快速乘法将表示的四重块Toeplitz矩阵和向量形式表示的随机矩阵相乘，以对表示的四重块Toeplitz矩阵进行矩阵维数压缩。
一个实施例中，利用随机QR降秩分解算法对所述压缩矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理，包括：对所述压缩矩阵实施QR分解，得到由正交矩阵和上三角矩阵的乘积表示的压缩矩阵，并存储所述正交矩阵；将所述正交矩阵的共轭转置矩阵和表示的四重块Toeplitz矩阵相乘得到子矩阵，存储所述子矩阵，并利用所述正交矩阵和所述子矩阵的乘积表示降秩后的所述四重块Toeplitz矩阵。
一个实施例中，利用无展开求平均算法对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均，包括：利用所述正交矩阵和所述子矩阵计算降秩后的所述四重块Toeplitz矩阵的对角线元素的平均。
一个实施例中，对所述地震频率切片数据的重建数据做傅里叶反变换，得到时间域地震重建数据之前，还包括：计算对角线元素求平均的结果和所述地震频率切片数据的差值并判断该差值是否在设定误差范围内；若否，利用对角线元素求平均的结果重新获取所述地震频率切片数据的四重块Toeplitz矩阵，并存储重新获取的四重块Toeplitz矩阵的部分行列元素以表示重新获取的四重块Toeplitz矩阵，重新获取的部分行列元素包括重新获取的四重块Toeplitz矩阵的各不同元素；利用随机QR降秩分解算法对重新获取的部分行列元素表示的重新获取的四重块Toeplitz矩阵进行降秩处理，以对重新获取的四重块Toeplitz矩阵进行降秩处理；利用无展开求平均算法对降秩处理后的重新获取的四重块Toeplitz矩阵的对角线元素求平均，得到重新获取的对角线元素平均；计算重新获取的对角线元素平均和所述对角线元素求平均的结果的差值并判断该差值是否在所述设定误差范围内，若是，将重新获取的对角线元素平均作为所述地震频率切片数据的重建数据，若否，利用重新获取的对角线元素平均再次重新获取所述地震频率切片数据的四重块Toeplitz矩阵，并存储再次重新获取的四重块Toeplitz矩阵的部分行列元素以表示再次重新获取的四重块Toeplitz矩阵，再次重新获取的部分行列元素包括再次重新获取的四重块Toeplitz矩阵的各不同元素，利用随机QR降秩分解算法对再次重新获取的部分行列元素表示的再次重新获取的四重块Toeplitz矩阵进行降秩处理，以对再次重新获取的四重块Toeplitz矩阵进行降秩处理，利用无展开求平均算法对降秩处理后的再次重新获取的四重块Toeplitz矩阵的对角线元素求平均，计算再次重新获取的对角线元素平均和所述重新获取的对角线元素平均的差值，依次迭代进行，直到所得差值在所述设定误差范围内。
本发明实施例的基于矩阵降秩的地震数据重建方法，通过仅存储四重块Toeplitz矩阵的部分行列元素来表示整个四重块Toeplitz矩阵能够降低地震数据重建过程中的数据存储量；利用随机QR降秩分解算法可以实现对部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以此可以减少降秩计算量，提高地震数据重建速度；利用无展开求平均算法对降秩处理后的四重块Toeplitz矩阵的对角线元素求平均值，无需计算并存储降秩处理后的所述四重块Toeplitz矩阵的具体形式，无需存储大型的低秩四重块Toeplitz矩阵，可以减少数据存储量，提高地震数据重建效率。本发明的方法能够减少数据储存量和计算量且不降低重建质量。本发明实施例的方法不仅能够用于五维地震数据的重建，还可以很好地用于N维地震数据的重建，N≥2。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1是本发明实施例的基于矩阵降秩的地震数据重建方法的流程示意图；
图2是本发明一实施例中获取地震频率切片数据的四重块Toeplitz矩阵的方法流程示意图；
图3是本发明另一实施例中获取地震频率切片数据的四重块Toeplitz矩阵的方法流程示意图；
图4是本发明一实施例中利用随机QR降秩分解算法对部分行列元素表示的四重块Toeplitz矩阵进行降秩处理的方法流程示意图；
图5是本发明一实施例中利用随机QR降秩分解算法对压缩矩阵进行降秩处理的方法流程示意图；
图6是本发明另一实施例的基于矩阵降秩的地震数据重建方法的流程示意图；
图7是本发明另一个实施例的基于矩阵降秩的地震数据重建方法的流程示意图；
图8是分别利用现有方法和本发明一实施例的基于矩阵降秩的地震数据方法进行重建所耗费计算时间的对比示意图；
图9是分别利用现有方法和本发明一实施例的基于矩阵降秩的地震数据方法进行重建的重建质量因子对比示意图；
图10是本发明一实施例中原始完整数据的cmpx面元切片示意图；
图11是本发明一实施例中50％不规则缺失道数据的cmpx面元切片示意图；
图12是本发明一实施例中重建数据的cmpx面元切片示意图；
图13是图10所示原始完整数据与图12所示重建数据的差剖面；
图14是本发明一实施例中原始实际采集的观测数据的cmpy面元切片示意图；
图15是利用现有基于Lanczos双对角分解的多道奇异谱分析方法的重建结果的cmpy面元切片示意图；
图16是利用本发明实施例的方法的重建结果的cmpy面元切片示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
在现有矩阵降秩类重建方法中最具代表性的是基于随机SVD分解的多道奇异谱分析方法，该方法可以归结为：(1)对(t,x1,x2)域3D地震数据做关于时间t的傅里叶变换到频率-空间域，得到(f,x1,x2)域数据，这里x1和x2表示空间维变量；(2)对每一个频率切片数据构建二重块Hankel矩阵；(3)对二重块Hankel矩阵利用随机SVD分解算法进行降秩，得到降秩后的二重块Hankel矩阵；(4)对降秩后的低秩二重块Hankel矩阵沿反对角线元素求平均得到重建数据；(5)对该频率切片重复步骤(2)，(3)和(4)，直到该切片满足重建误差要求；(6)对所有频率成分按步骤(2)-(5)进行重建，最后对频率f做傅里叶反变换到时间域完成重建。
发明人发现，基于随机SVD分解的多道奇异谱分析方法存在如下缺陷：
1、该技术方案的第二歩需要生成并存储二重块Hankel矩阵，对于实际数据在重建时需要构建大型块Hankel矩阵，导致存储量大幅增加，占用大量内存。然而，Hankel矩阵具有沿反对角线元素相等的特殊结构，只需存储矩阵的第1行和最后1列就可以表示整个矩阵，不需要将整个矩阵生成并且存储它的全部元素。例如，对于一个大小为10×10的Hankel矩阵，若直接生成和存储该矩阵需要存储100个元素，而利用本发明的方案只需存储19个元素。
2、在步骤(3)矩阵降秩阶段，该技术采用随机SVD分解算法，虽然该算法与标准的截断SVD分解算法相比计算量减少，但是仍需依赖SVD分解，对大型块Hankel矩阵而言，降秩计算量依然很大。例如，对一个秩为k,大小为m×n的二重块Hankel矩阵(m≤n)利用该技术进行降秩，计算量为8mn+8nk2+5mk2,而本发明提出的随机QR分解技术计算量为km2+k(m+n-1)log2(m+n-1)。若取m＝n＝10，k＝3，则该技术的运算量为3570次，而本发明的技术只需585次运算。
3、在步骤(4)沿反对角线求平均的过程中，该技术是先生成和存储降秩后的低秩二重块Hankel矩阵，这会占用大量内存。其实该二重块Hankel矩阵可以不用存储，沿反对角线进行求平均的元素可以直接从降秩分解后的子矩阵中提取。
4、该技术只能对叠前3D地震数据进行插值，无法对高维的5D数据进行插值重建处理。
基于Lanczos双对角分解的多道奇异谱分析重建方法(A fast reduced-rankinterpolation method for prestack seismic volumes that depend on four spatialdimensions,Geophysics,2013,78(1):V21–V30)，采用Lanczos双对角分解算法来取代截断SVD分解算法，提高降秩计算效率，并且采用四维快速傅里叶变换来实现四重块Toeplitz矩阵与向量的快速乘法，避免矩阵与向量直接相乘。但是，该技术在对降秩后的低秩四重块Toeplitz矩阵沿对角线元素求平均时需要先生成并存储该大型四重块Toeplitz矩阵，然后再沿对角线求平均，得到重建数据。生成并存储该大型四重块Toeplitz矩阵必然会占用大量内存，而这一过程完全可以省略，直接从降秩分解得到的子矩阵中提取出对角线元素并求平均，从而节约数据存储量。该技术的实施方案可以归结为：(1)将(t,x1,x2,x3,x4)域5D地震数据做关于时间t的傅里叶变换到频率域得(f,x1,x2,x3,x4)域地震数据，这里x1,x2,x3,x4表示空间维变量，既可以指炮点x坐标Sx，炮点y坐标Sy，检波器x坐标rx和检波器y坐标ry，也可以指共中心点x坐标CMPx，共中心点y坐标CMPy，偏移距x坐标hx和偏移距y坐标hy)；(2)对每一个频率切片数据构建四重块Toeplitz矩阵；(3)对该四重块Toeplitz矩阵利用Lanczos双对角分解进行降秩，得到降秩后的低秩四重块Toeplitz矩阵；(4)对低秩四重块Toeplitz矩阵沿对角线元素求平均，得到该地震频率切片对应的重建数据；(5)对该地震频率切片重复步骤(2)，(3)和(4)，直到该切片满足重建误差要求，完成对该频率切片的数据重建；(6)对所有的频率成分按步骤(2)-(5)进行重建，最后对频率f做傅里叶反变换到时间域，完成重建。
发明人发现，基于Lanczos双对角分解的多道奇异谱分析重建方法存在如下缺陷：
1、该技术直接对由频率切片数据构建的四重块Toeplitz矩阵实施降秩，而四重块Toeplitz矩阵元素存在冗余，完全可以采用类似于现有技术1的做法，先对该大型四重Toeplitz矩阵乘以一个随机矩阵进行矩阵压缩，然后对压缩后的矩阵进行降秩分解，进而可以减小计算量，但是该技术未采用该做法。
2、该技术步骤(3)采用Lanczos双对角分解算法进行降秩，与传统的截断SVD分解算法相比，计算效率有提高，但还可以进一步被提高。对于一个秩为k的m×n型四重块Toeplitz矩阵，Lanczos分解的降秩计算量为O(km2)阶，若采用本发明的随机QR分解算法，计算量仅为O(mk2)，这里k＜＜m。
3、在步骤(4)对降秩后的低秩四重块Toeplitz矩阵沿对角线求平均时该技术仍需先生成一个大型四重块Toeplitz矩阵，然后再沿对角线求平均，此处生成和存储该四重块Toeplitz矩阵需要占用大量存储空间。本发明采用无展开求平均技术，在求平均步骤可以不用生成和存储大型该四重块Toeplitz矩阵，而是直接从降秩分解的子矩阵中提取出对角线元素并对它们求平均。
基于上述发现，针对叠前5D重建面临的降秩分解计算量大，计算成本高以及数据存储量大的问题，本发明提出了一种基于矩阵降秩的地震数据重建方法。该方法是一种非SVD分解方法，可以有效解决矩阵降秩因严重依赖SVD分解算法导致的计算量大、计算成本高、无法应用于大规模工业实际生产数据插值重建等问题。
图1是本发明实施例的基于矩阵降秩的地震数据重建方法的流程示意图。如图1所示，本发明实施例的基于矩阵降秩的地震数据重建方法，可包括步骤：
S110：获取地震频率切片数据的四重块Toeplitz矩阵，并存储所述四重块Toeplitz矩阵的部分行列元素以表示所述四重块Toeplitz矩阵，所述部分行列元素包括所述四重块Toeplitz矩阵的各不同元素；
S120：利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理；
S130：利用无展开求平均算法对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均，得到所述地震频率切片数据的重建数据；
S140：对所述地震频率切片数据的重建数据做傅里叶反变换，得到时间域地震重建数据，用于油气勘探。
在上述步骤S110中，该地震频率切片数据可以是频率域某频率下的地震数据。该四重块Toeplitz矩阵可以由地震频率切片数据直接构建，或由地震频率切片数据构建的其他类型矩阵，例如四重块Hankel矩阵，通过变换得到。可以根据四重块Toeplitz矩阵的结构特点，仅存储四重块Toeplitz矩阵中的部分行列元素，该部分行列元素中的元素可以覆盖四重块Toeplitz矩阵中的各不同元素。利用该部分行列元素可以表示整个四重块Toeplitz矩阵。
在上述步骤S120中，利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理时，可以仅涉及上述部分行列元素，例如仅涉及四重块Toeplitz矩阵的第一行及第一列，以此可以实现对整个四重块Toeplitz矩阵进行降秩处理。
在上述步骤S130中，对降秩处理后的所述部分行列元素表示的四重块Toeplitz矩阵的对角线元素求平均，既可以实现对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均。降秩处理后的低秩四重块Toeplitz矩阵可以用两个或多个矩阵相乘的形式表示，利用无展开求平均算法可以直接计算得到降秩处理后的低秩四重块Toeplitz矩阵对角线元素的平均值，而无需先计算得到降秩处理后的低秩四重块Toeplitz矩阵的具体形式后再对具体形式中的对角线元素求平均值。若降秩处理后的低秩四重块Toeplitz矩阵的对角线元素的平均值满足重建条件，则可得到所述地震频率切片数据的重建数据。
在上述步骤S140中，地震勘探数据可以对应多个上述地震频率切片数据。每个地震频率切片数据可以利用上述步骤S110～步骤S130得到相应的重建数据。对各地震频率切片数据的重建数据做傅里叶反变换，可得到时间域地震重建数据。该时间域地震重建数据可以填补最初的地震缺失道数据，并可以压制包含在最初地震勘探数据中的随机噪声，从而使该时间域地震重建数据相比于原来含有缺失道地震数据具有较高的品质，将重建后的数据应用于油气勘探，勘探结果会更准确。
本发明实施例中，通过仅存储四重块Toeplitz矩阵的部分行列元素来表示整个四重块Toeplitz矩阵能够降低地震数据重建过程中的数据存储量；利用随机QR降秩分解算法可以实现对部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以此可以减小降秩计算量，提高地震数据重建速度；利用无展开求平均算法对降秩处理后的四重块Toeplitz矩阵的对角线元素求平均，无需计算并存储降秩处理后的低秩四重块Toeplitz矩阵的具体形式，无需存储大型的低秩四重块Toeplitz矩阵，可以减少数据存储量，提高地震数据重建效率。
一些实施例中，所述部分行列元素可为所述四重块Toeplitz矩阵的第一行元素及第一列元素或最后一行元素及最后一列元素。利用四重块Toeplitz矩阵的第一行元素及第一列元素或四重块Toeplitz矩阵的最后一行元素及最后一列元素即可表示整个四重块Toeplitz矩阵。
本实施例中，仅存储四重块Toeplitz矩阵的第一行元素及第一列元素或最后一行元素及最后一列元素，可以显著降低数据存储量。
其他一些实施例中，在上述步骤S110中，可获取由地震频率切片数据构建的四重块Hankel矩阵，可以只需存储四重块Hankel矩阵的第1行元素和最后1列元素。本实施例中，针对四重块Hankel矩阵沿反对角线元素相等的结构特征，只存储矩阵的第一行和最后一列，可以减少数据存储量。
图2是本发明一实施例中获取地震频率切片数据的四重块Toeplitz矩阵的方法流程示意图。如图2所示，在上述步骤S110中，获取地震频率切片数据的四重块Toeplitz矩阵的方法，可包括步骤：
S111：对时间域地震数据做关于时间的傅里叶变换，得到频率域地震数据，并从所述频率域地震数据中提取得到所述地震频率切片数据；
S112：利用所述地震频率切片数据构建所述四重块Toeplitz矩阵，或利用所述地震频率切片数据构建四重块Hankel矩阵，并将所述四重块Hankel矩阵变换为所述四重块Toeplitz矩阵。
在上述步骤S111中，可以通过各种不同方式从所述频率域地震数据中提取得到所述地震频率切片数据，例如可以根据傅里叶变换的长度参数、时间域地震数据沿时间维的采样点数和采样间隔等确定频率切片的提取间隔，可以按频率从小到大的顺序提取地震频率切片数据。
在上述步骤S112中，通过多种不同方法得到的所述地震频率切片数据的所述四重块Toeplitz矩阵，可以充分利用四重块Toeplitz矩阵的结构特点，可仅存储四重块Toeplitz矩阵的第一行元素及第一列元素或最后一行元素及最后一列元素，从而减少四重块Toeplitz矩阵的元素存储量。
图3是本发明另一实施例中获取地震频率切片数据的四重块Toeplitz矩阵的方法流程示意图。如图3所示，图2所示的获取地震频率切片数据的四重块Toeplitz矩阵的方法，在步骤S111之前，即对时间域地震数据做关于时间的傅里叶变换，得到频率域地震数据之前，还可包括步骤：
S113：通过坐标变换将地震勘探数据由炮-检域转换到共中心点-偏移距域，得到所述时间域地震数据。
本实施例中，通过将地震勘探数据由炮-检域转换到共中心点-偏移距域，并基于共中心点-偏移距域的时间域地震数据进行后续降秩处理，可以充分利用共中心点-偏移距域的时间域地震数据的低秩特征，进一步减少四重块Toeplitz矩阵降秩处理的计算量。
一些实施例中，所述时间域地震数据为5D地震数据。本实施例中，对5D(五维)地震数据进行重建可以利用更多维的已知地震道信息，将这些信息作为约束来进行重建，可以获得更加精确的地震波场重建结果。在其他实施例中，本发明的方法可以用于其他各种维数的地震数据的重建，例如3D(三维)地震数据。
图4是本发明一实施例中利用随机QR降秩分解算法对部分行列元素表示的四重块Toeplitz矩阵进行降秩处理的方法流程示意图。如图4所示，在上述步骤S120中，利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理的方法，可包括步骤：
S121：通过将表示的四重块Toeplitz矩阵乘以一随机矩阵对表示的四重块Toeplitz矩阵进行矩阵压缩，得到压缩矩阵，所述随机矩阵的行数和列数分别为所述四重块Toeplitz矩阵的列数和秩；
S122：利用随机QR降秩分解算法对所述压缩矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理。
在上述步骤S121中，该随机矩阵可需满足两个条件，其一是该矩阵的各列线性独立，互不相干；其二是该矩阵的行数和列数分别为所述四重块Toeplitz矩阵的列数和秩。该随机矩阵具有较小的维数，通过将上述部分行列元素表示的四重块Toeplitz矩阵乘以该随机矩阵，可以对上述部分行列元素表示的四重块Toeplitz矩阵进行矩阵压缩，得到压缩矩阵。在上述步骤S122中，与直接对四重块Toeplitz矩阵进行降秩处理相比，利用随机QR降秩分解算法对所述压缩矩阵进行降秩处理，能够减少对四重块Toeplitz矩阵降秩的冗余，从而可以进一步减少将秩处理的计算量。
一些实施例中，在上述步骤S122中，通过将表示的四重块Toeplitz矩阵乘以一随机矩阵对表示的四重块Toeplitz矩阵进行矩阵压缩的方法，具体实施方式可为：用向量形式表示所述随机矩阵，并利用基于四维傅里叶快速变换的快速乘法将表示的四重块Toeplitz矩阵和向量形式表示的随机矩阵相乘，以对表示的四重块Toeplitz矩阵进行矩阵维数压缩。
本实施例中，利用基于四维傅里叶快速变换的快速乘法计算上述部分行列元素表示的四重块Toeplitz矩阵和向量形式表示的随机矩阵的乘积，能够提高对四重块Toeplitz矩阵进行矩阵压缩的计算速度。
一个实施例中，上述实施例的基于四维傅里叶快速变换的快速乘法可以参考“Afast reduced-rank interpolation method for prestack seismic volumes thatdepend on four spatial dimensions”(Geophysics,2013,78(1)：V21–V30)中基于4D FFT的四重块Toeplitz矩阵和向量的快速乘法原理实施。
图5是本发明一实施例中利用随机QR降秩分解算法对压缩矩阵进行降秩处理的方法流程示意图。如图5所示，在上述步骤S122中，利用随机QR降秩分解算法对所述压缩矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理的方法，可包括步骤：
S1221：对所述压缩矩阵实施QR分解，得到由正交矩阵和上三角矩阵的乘积表示的压缩矩阵，并存储所述正交矩阵；
S1222：将所述正交矩阵的共轭转置矩阵和表示的四重块Toeplitz矩阵相乘得到子矩阵，存储所述子矩阵，并利用所述正交矩阵和所述子矩阵的乘积表示降秩后的所述四重块Toeplitz矩阵。
在上述步骤S1221中，对上述压缩矩阵实施QR分解，可以分解得到正交矩阵，记为Q和上三角矩阵，记为R，利用该正交矩阵Q和该上三角矩阵R的乘积形式可以表示该压缩矩阵。在上述步骤S1222中，利用该正交矩阵R、该正交矩阵的共轭转置矩阵，记为Q*，及上述部分行列元素表示的四重块Toeplitz矩阵三者的乘积可以表示降秩后的低秩四重块Toeplitz矩阵。而该正交矩阵和上述子矩阵维数较小，通过存储该正交矩阵和该子矩阵，并利用该正交矩阵和该子矩阵表示降秩后的所述四重块Toeplitz矩阵，可以避免直接计算和存储该大型低秩四重块Toeplitz矩阵，以此可以进一步减少数据存储量和计算量。
一些实施例中，可以利用基于四维傅里叶快速变换的矩阵和向量快速乘法将所述正交矩阵的共轭转置矩阵Q*和表示的四重块Toeplitz矩阵相乘得到子矩阵，以此可以进一步提高地震数据重建的速度。
一些实施例中，在上述步骤S130中，利用无展开求平均算法对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均的方法，具体实施方式可以表示为：利用所述正交矩阵和所述子矩阵计算降秩后的所述四重块Toeplitz矩阵的对角线元素的平均。
本实施例中，在不生成和不存储降秩后的大型低秩四重块Toeplitz矩阵的情况下，利用事先存储的上述正交矩阵和上述子矩阵直接计算降秩后的所述低秩四重块Toeplitz矩阵的对角线元素平均值，无需计算及存储大型的降秩后的低秩四重块Toeplitz矩阵，以此可以进一步减少数据存储量和计算量。
一个实施例中，对降秩后的所述四重块Toeplitz矩阵沿对角线元素求平均值可表示为：
其中，
在上式中，D(fi,i1,i2,i3,i4)表示降秩后沿所述四重块Toeplitz矩阵的对角线元素求取的平均值，fi表示第i个地震频率切片数据的频率，i1、i2、i3及i4表示第i个地震频率切片数据分别在共中心点x坐标cmpx、共中心点y坐标cmpy、偏移距x坐标hx及偏移距y坐标hy的数值；为根据正交矩阵Q构建的四重块Toeplitz矩阵，为正交矩阵Q的共轭转置矩阵Q*和四重块Toeplitz矩阵可记为A(4)相乘得到的子矩阵M的列向量形式，M＝Q*A(4)，r为所述四重块Toeplitz矩阵的秩，l为正整数，1≤l≤r；四重块Toeplitz矩阵大小为N4×K4，其中的元素为根据正交矩阵构建的N3×K3的三重块Toeplitz矩阵，p4＝1,2,...,L4；该三重块Toeplitz矩阵的每个元素为根据正交矩阵构建的N2×K2的二重块Toeplitz矩阵，p3＝1,2,...,L3；该二重块Toeplitz矩阵的每个元素为根据正交矩阵构建的N1×K1的Toeplitz矩阵，p2＝1,2,...,L2；该Toeplitz矩阵的每个元素是正交矩阵Q中第p1行且第l列的元素，p1＝1,2,...,L1；列向量大小为K4×1，其中的每个元素是根据子矩阵M构建的大小为K3×1的列向量，ν4＝1,2,...,K4；向量的每个元素是根据子矩阵M构建的大小为K2×1的列向量，ν3＝1,2,...,K3；向量中的每个元素是根据子矩阵M构建的大小为K1×1的列向量，ν2＝1,2,...,K2；向量中的每个元素是子矩阵M的第l行且第ν1列元素，ν1＝1,2,...,K1；Lj+Kj-1＝Nj，符号表示取整数部分，Nj为地震数据D(fi,i1,i2,i3,i4)在第j维空间的方向长度，j＝1,2,3,4。
其他一些实施例中，由地震频率切片数据构建的四重块Hankel矩阵，可以不存储该四重块Hankel矩阵，只需存储该四重块Hankel矩阵的第一行和最后一列。而沿四重块Hankel矩阵反对角线进行求平均的元素可以直接从降秩分解后的低秩四重块Hankel矩阵的子矩阵中提取，然后求平均。本实施例中，直接从降秩分解得到的正交矩阵Q和子矩阵M中提取出对角线元素并求平均，可以减少数据存储量和计算量。
图6是本发明另一实施例的基于矩阵降秩的地震数据重建方法的流程示意图。如图6所示，图1所示的基于矩阵降秩的地震数据重建方法，在步骤S140之前，即对所述地震频率切片数据的重建数据做傅里叶反变换，得到时间域地震重建数据之前，还可包括步骤：
S150：计算对角线元素求平均的结果和所述地震频率切片数据的差值并判断该差值是否在设定误差范围内；
S160：若否，利用对角线元素求平均的结果重新获取所述地震频率切片数据的四重块Toeplitz矩阵，并存储重新获取的四重块Toeplitz矩阵的部分行列元素以表示重新获取的四重块Toeplitz矩阵，重新获取的部分行列元素包括重新获取的四重块Toeplitz矩阵的各不同元素；
S170：利用随机QR降秩分解算法对重新获取的部分行列元素表示的重新获取的四重块Toeplitz矩阵进行降秩处理，以对重新获取的四重块Toeplitz矩阵进行降秩处理；
S180：利用无展开求平均算法对降秩处理后的重新获取的四重块Toeplitz矩阵的对角线元素求平均，得到重新获取的对角线元素平均；
S190：计算重新获取的对角线元素平均和所述对角线元素求平均的结果的差值并判断该差值是否在所述设定误差范围内，若是，将重新获取的对角线元素平均作为所述地震频率切片数据的重建数据，若否，利用重新获取的对角线元素平均再次重新获取所述地震频率切片数据的四重块Toeplitz矩阵，并存储再次重新获取的四重块Toeplitz矩阵的部分行列元素以表示再次重新获取的四重块Toeplitz矩阵，再次重新获取的部分行列元素包括再次重新获取的四重块Toeplitz矩阵的各不同元素，利用随机QR降秩分解算法对再次重新获取的部分行列元素表示的再次重新获取的四重块Toeplitz矩阵进行降秩处理，以对再次重新获取的四重块Toeplitz矩阵进行降秩处理，利用无展开求平均算法对降秩处理后的再次重新获取的四重块Toeplitz矩阵的对角线元素求平均，计算再次重新获取的对角线元素平均和所述重新获取的对角线元素平均的差值，依次迭代进行，直到所得差值在所述设定误差范围内。
本实施例中，当后一次计算得到的对角线元素平均值与前一次计算得到的对角线元素平均值(或地震频率切片数据)的差值不满足重建条件，即不小于或等于设定误差范围，例如10-4，则利用迭代求解方法，继续对重建的结果进行重建直到满足重建条件，以此可以获得更佳的地震数据重建结果。
图7是本发明另一个实施例的基于矩阵降秩的地震数据重建方法的流程示意图。如图7所示，本发明实施例的基于矩阵降秩的地震数据重建方法，可包括步骤：
(1)对时间域5D地震数据做关于时间t的傅里叶变换到频率域；
(2)对每一个频率切片数据构建四重块Toeplitz矩阵并只存储该四重块Toeplitz矩阵的第一行和第一列元素来表示整个矩阵，以此可以达到节约计算机内存，减小数据存储量的目的；
(3)采用随机QR快速降秩分解算法对四重块Toeplitz矩阵进行降秩；
(4)采用无展开求平均算法对降秩后的四重块Toeplitz矩阵沿对角线元素求平均，得到该地震频率切片的重建数据；
(5)对该地震频率切片按照步骤(2)、(3)和(4)进行迭代重建，直至满足重建误差条件；
(6)对所有频率成分按步骤(2)-(5)进行重建，最后做关于频率f的傅里叶反变换到时间域得到时间域重建数据，完成重建。
一个实施例中，利用随机QR分解进行矩阵降秩的地震数据重建方法，可包括步骤：
(1)在共中心点-偏移距域进行地震数据重建。首先可通过坐标变换将3D地震勘探中采集的叠前炮-检域5D数据D(t,sx,sy,rx,ry)转换到共中心点-偏移距域的5D数据D(t,cmpx,cmpy,hx,hy)。例如，坐标变换公式可表示为 hx＝sx-rx和hy＝sy-ry，然后再进行重建；
(2)对于共中心点-偏移距域的5D数据D(t,cmpx,cmpy,hx,hy)做关于时间t的傅里叶变换到频率域得到频率域5D地震数据D(f,cmpx,cmpy,hx,hy)；
(3)提取第i个频率切片数据D(fi,cmpx,cmpy,hx,hy)，设置迭代次数初始值k＝1，对该频率切片数据进行迭代重建，fi是第i个频率切片数据对应的频率值；
(4)将由频率切片数据D(fi,cmpx,cmpy,hx,hy)构建的四重块Toeplitz矩阵记为A(4)，充分挖掘和利用四重块Toeplitz矩阵沿对角线元素相等的结构特征，只存储矩阵A(4)的第一行和第一列来表示整个四重块Toeplitz矩阵A(4)，以此可以达到压缩数据存储量的目的；
(5)对第一行和第一列表示的四重块Toeplitz矩阵A(4)运用随机QR分解进行降秩，得到频率切片fi对应的重建数据。具体过程例如可为：设A(4)是一个秩为r(r<m≤n)，大小为m×n的四重块Toeplitz矩阵，Ω是一个列向量互不相干且大小为n×r的随机矩阵。将四重块Toeplitz矩阵A(4)与随机矩阵Ω相乘得到压缩后的矩阵B，B＝A(4)Ω，此时压缩后的矩阵B是一个维数较小的m×r型压缩矩阵。在计算压缩后的矩阵B时，可运用基于四维傅里叶快速变换的四重块Toeplitz矩阵和向量的快速乘法来实现矩阵A(4)和Ω的相乘，而不直接将A(4)与Ω相乘，以此可提高计算效率。对压缩后的矩阵B实施QR分解得B＝QR，Q为正交矩阵，R为上三角矩阵。将降秩后的低秩四重块Toeplitz矩阵记为则取M＝Q*A(4)，则低秩四重块Toeplitz矩阵需要说明的是此处只是一个符号标记，具体实施时不需要真正计算和存储该低秩四重块Toeplitz矩阵而只需计算和存储子矩阵Q和M。在计算子矩阵M时可运用基于四维快速傅里叶变换的四重块Toeplitz矩阵和向量的快速乘法技术来实现正交矩阵Q的共轭转置矩阵Q*和四重块Toeplitz矩阵A(4)相乘。该快速乘法可参照“A fast reduced-rank interpolation method for prestack seismic volumesthat depend on four spatial dimensions”(Geophysics,2013,78(1)：V21–V30)。对低秩四重块Toeplitz矩阵可采用无展开快速求平均算法，沿低秩四重块Toeplitz矩阵对角线元素求平均值，得到重建数据。无展开快速求平均算法例如可以表示为，
其中，变量i1,i2,i3和i4分别表示空间坐标cmpx,cmpy,hx和hy，设Nj表示地震数据D在第j维空间方向的长度，Lj和Kj满足Lj+Kj-1＝Nj，符号表示取整数部分，j＝1,2,3,4；为一个四重块Toeplitz矩阵，其首行元素由一个大小为K4×1的向量构成，其首列元素由一个大小为N4×1的向量 构成，为一个大小为K4×1的列向量 矩阵和向量中元素具体表达式为，
其中中的每个元素p4＝1,2,...,L4为一个三重块Toeplitz矩阵，该三重块Toeplitz矩阵中每个元素又可表示为p3＝1,2,...,L3。为一个二重块Toeplitz矩阵，该二重块Toeplitz矩阵中的每个元素可表示为p2＝1,2,...,L2。为一个普通Toeplitz矩阵，该Toeplitz矩阵中的每个元素可表示为p1＝1,2,...,L1。是位于正交矩阵Q中第p1行和第l列的元素。同样，每个ν4＝1,2,...,K4都是一个K3×1的列向量。向量中每个元素ν3＝1,2,...,K3又是一个K2×1的列向量，中每个元素ν2＝1,2,...,K2是一个K1×1的列向量。向量中每个元素ν1＝1,2,...,K1取自子矩阵M的第l行和第ν1列元素。计算可运用基于四维快速傅里叶变换的四重块Toeplitz矩阵和向量快速乘法技术来实现，该快速乘法的具体计算原理可参照“A fast reduced-rank interpolation method for prestack seismic volumesthat depend on four spatial dimensions”(Geophysics,2013,78(1)：V21–V30)。对D(fi,i1,i2,i3,i4)中的缺失点数据而言，做一次降秩分解重建后该点的振幅能量值并不能完全被恢复，需要不断按照步骤(3)、(4)和(5)进行反复迭代重建，直至满足重建误差条件。
一个实施例中，对频率fi对应的地震频率切片数据而言，第k次迭代重建的表达式可表示为，
D(k+1)(fi,i1,i2,i3,i4)＝αDobs(fi,i1,i2,i3,i4)+(I-αS)Dk(fi,i1,i2,i3,i4),k＝1,...,Niter.
其中，Niter为最大迭代次数，Dobs(fi,i1,i2,i3,i4)为原始输入的观测数据，α为权值因子，α∈(0,1]，S为采样算子，其元素为0和1，其中0表示缺失点，1表示已知点；Dk(fi,i1,i2,i3,i4)表示第k次迭代的重建数据，D(k+1)(fi,i1,i2,i3,i4)表示第k+1迭代的重建数据，I表示元素值全为1的四维数组。
(6)对所有频率切片数据D(fi,i1,i2,i3,i4)按步骤(2)-(5)进行重建。最后，对频率f做傅里叶反变换到时间域，得到时间域重建数据，最终完成重建。
本发明实施例的方法，因对四重块Toeplitz矩阵采用了随机QR快速分解降秩方法，无展开求平均方法，基于四维傅里叶快速变换的四重块Toeplitz矩阵和向量快速乘法计算技术以及依据四重块Toeplitz矩阵结构特点只存储矩阵的第一行和第一列来表示整个四重块Toeplitz矩阵的有效存储方法，使得矩阵降秩重建的计算效率得到明显提升。在提升计算速度的同时还能保持重建质量。本发明实施例的方法是一种非SVD降秩重建方法，可以有效解决目前矩阵降秩重建方法因严重依赖SVD分解，计算量大而无法应用于大规模工业数据插值重建处理的瓶颈问题。
一个实施例中，设计不同大小的叠前5D无噪声模型数据D(t,cmpx,cmpy,hx,hy)，分别用基于随机SVD分解的多道奇异谱分析技术(Rsvds方法)，基于Lanczos双对角分解的多道奇异谱分析技术(Lanczos方法)和本发明实施例的基于无展开随机QR分解矩阵降秩技术进行重建计算效率比较。每个数据在时间方向含有nt个采样点，例如nt＝256，采样间隔dt＝0.004s。在空间方向的大小为N×N×N×N，共包含三个不同斜率的同相轴。可取N＝6,7,…,13，最大迭代次数Niter＝10，秩r＝6，加权因子α＝1.0，迭代停止误差ε＝10-4，然后用三种方法进行重建。图8是分别利用现有方法和本发明一实施例的利用随机QR分解进行矩阵降秩的地震数据重建方法进行重建所耗费计算时间的对比示意图。如图8所示，本发明实施例的无展开随机QR分解降秩重建方法要明显快于现有的Rsvds方法和Lanczos方法两种重建方法，并且随着数据维数N的逐渐增大，三种重建方法的计算时间差异也更加明显。取N＝10时的5D模型数据D(t,cmpx,cmpy,hx,hy)，然后随机剔除10％,20％,…,80％的地震道形成不规则缺失道地震数据，分别运用三种方法进行重建并比较重建效果。定义重建质量因子Q，Drecon表示重建数据，Dtrue表示原始真实数据，F表示四维数组的斐波那契范数。图9是分别利用现有方法和本发明一实施例的利用随机QR分解进行矩阵降秩的地震数据重建方法进行重建的重建质量因子对比示意图。如图9所示，本发明实施例的无展开随机QR分解降秩重建方法和现有的基于随机QR分解的多道奇异谱分析方法具有相似的重建质量，但是前者比后者计算速度快，计算效率高。而且从图9还可以看出现有的基于Lanczos双对角分解的多道奇异谱分析方法的重建质量明显低于本发明实施例的方法。
一个实施例中，合成一个大小为256×10×10×10×10的5D无噪声模型数据，随机剔除50％的地震道形成不规则缺失道地震数据，然后用本发明实施例的方法进行重建。该数据在时间方向含有256个时间采样点，在空间方向包含10个CMP面元，10个偏移距，设置最大迭代次数Niter＝10次，秩r＝6，加权因子α＝1.0。图10是本发明一实施例中原始完整数据的cmpx面元切片示意图。如图10所示，共中心点x坐标cmpx＝1，3，5，7，9，共中心点y坐标cmpy＝2，偏移距x坐标hx＝1：10，偏移距y坐标hy＝2。图11是本发明一实施例中随机缺失50％地震道数据的cmpx面元切片示意图。图12是本发明一实施例中重建数据的cmpx面元切片示意图。图13是图10所示原始完整数据与图12所示重建数据的cmpx面元切片差剖面，如图13所示，从差剖面上可以看出缺失道数据的振幅值被有效恢复。
一个实施例中，将本发明实施例的方法应用于某油田陆上区块叠前5D实际数据插值重建处理。该数据大小为301×15×15×13×13，在时间方向含有301个采样点，在空间方向含有15个CMP面元，13个偏移距，数据整体缺失89％。设置最大迭代次数Niter＝100次，秩r＝10，加权因子α＝0.8，然后分别运用现有的基于Lanczos双对角分解的多道奇异谱分析方法和本发明实施例的基于无展开随机QR分解矩阵降秩重建方法对该共中心点-偏移距域不规则缺失道数据进行重建。图14是本发明一实施例中原始实际采集的观测数据的cmpy面元切片示意图。如图14所示，共中心点x坐标cmpx＝5，共中心点y坐标cmpy＝5，6，7，8，9，偏移距x坐标hx＝8和偏移距y坐标hy＝1:13。图15是利用现有基于Lanczos双对角分解的多道奇异谱分析方法的重建结果的cmpy面元切片示意图。图16是利用本发明实施例方法的重建结果的cmpy面元切片示意图。从图15和图16的矩形窗中可以看出，在相同的重建条件下本发明实施例的方法在重建后缺失道数据得到有效恢复并且同相轴的振幅能量也更强光滑连续。
本发明实施例的基于随机QR分解矩阵降秩的地震数据重建方法，通过仅存储四重块Toeplitz矩阵的部分行列元素来表示整个四重块Toeplitz矩阵能够降低地震数据在重建过程中的数据存储量；利用随机QR降秩分解算法可以实现对部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以此可以减小降秩计算量，提高地震数据重建效率；利用无展开求平均算法对降秩处理后的四重块Toeplitz矩阵的对角线元素求平均值，无需计算并存储降秩处理后的四重块Toeplitz矩阵的具体形式，无需存储该大型低秩四重块Toeplitz矩阵，可以有效减少数据存储量，提高地震数据重建效率。本发明的方法能够减少数据储存量和计算量且不会降低重建质量。本发明实施例的方法不仅能够用于五维地震数据的重建，还可以很好地用于N维地震数据的重建，N≥2。
在本说明书的描述中，参考术语“一个实施例”、“一个具体实施例”、“一些实施例”、“例如”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。各实施例中涉及的步骤顺序用于示意性说明本发明的实施，其中的步骤顺序不作限定，可根据需要作适当调整。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (8)
1.一种基于矩阵降秩的地震数据重建方法，其特征在于，包括：
获取地震频率切片数据的四重块Toeplitz矩阵，并存储所述四重块Toeplitz矩阵的部分行列元素以表示所述四重块Toeplitz矩阵，所述部分行列元素包括所述四重块Toeplitz矩阵的各不同元素；
利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理；
利用无展开求平均算法对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均，得到所述地震频率切片数据的重建数据；
对所述地震频率切片数据的重建数据做傅里叶反变换，得到时间域地震重建数据，用于油气勘探；
获取地震频率切片数据的四重块Toeplitz矩阵，包括：
对时间域地震数据做关于时间的傅里叶变换，得到频率域地震数据，并从所述频率域地震数据中提取得到所述地震频率切片数据；
利用所述地震频率切片数据构建所述四重块Toeplitz矩阵，或利用所述地震频率切片数据构建四重块Hankel矩阵，并将所述四重块Hankel矩阵变换为所述四重块Toeplitz矩阵；
利用随机QR降秩分解算法对所述部分行列元素表示的四重块Toeplitz矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理，包括：
通过将表示的四重块Toeplitz矩阵乘以一随机矩阵对表示的四重块Toeplitz矩阵进行矩阵压缩，得到压缩矩阵，所述随机矩阵的行数和列数分别为所述四重块Toeplitz矩阵的列数和秩；
利用随机QR降秩分解算法对所述压缩矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理。
2.如权利要求1所述的基于矩阵降秩的地震数据重建方法，其特征在于，对时间域地震数据做关于时间的傅里叶变换，得到频率域地震数据之前，还包括：
通过坐标变换将地震勘探数据由炮-检域转换到共中心点-偏移距域，得到所述时间域地震数据。
3.如权利要求1所述的基于矩阵降秩的地震数据重建方法，其特征在于，所述时间域地震数据为5D地震数据。
4.如权利要求1所述的基于矩阵降秩的地震数据重建方法，其特征在于，所述部分行列元素为所述四重块Toeplitz矩阵的第一行元素及第一列元素或最后一行元素及最后一列元素。
5.如权利要求1所述的基于矩阵降秩的地震数据重建方法，其特征在于，通过将表示的四重块Toeplitz矩阵乘以一随机矩阵对表示的四重块Toeplitz矩阵进行矩阵压缩，包括：
用向量形式表示所述随机矩阵，并利用基于四维傅里叶快速变换的快速乘法将表示的四重块Toeplitz矩阵和向量形式表示的随机矩阵相乘，以对表示的四重块Toeplitz矩阵进行矩阵维数压缩。
6.如权利要求1所述的基于矩阵降秩的地震数据重建方法，其特征在于，利用随机QR降秩分解算法对所述压缩矩阵进行降秩处理，以对所述四重块Toeplitz矩阵进行降秩处理，包括：
对所述压缩矩阵实施QR分解，得到由正交矩阵和上三角矩阵的乘积表示的压缩矩阵，并存储所述正交矩阵；
将所述正交矩阵的共轭转置矩阵和表示的四重块Toeplitz矩阵相乘得到子矩阵，存储所述子矩阵，并利用所述正交矩阵和所述子矩阵的乘积表示降秩后的所述四重块Toeplitz矩阵。
7.如权利要求6所述的基于矩阵降秩的地震数据重建方法，其特征在于，利用无展开求平均算法对降秩处理后的所述四重块Toeplitz矩阵的对角线元素求平均，包括：
利用所述正交矩阵和所述子矩阵计算降秩后的所述四重块Toeplitz矩阵的对角线元素的平均。
8.如权利要求1所述的基于矩阵降秩的地震数据重建方法，其特征在于，对所述地震频率切片数据的重建数据做傅里叶反变换，得到时间域地震重建数据之前，还包括：
计算对角线元素求平均的结果和所述地震频率切片数据的差值并判断该差值是否在设定误差范围内；
若否，利用对角线元素求平均的结果重新获取所述地震频率切片数据的四重块Toeplitz矩阵，并存储重新获取的四重块Toeplitz矩阵的部分行列元素以表示重新获取的四重块Toeplitz矩阵，重新获取的部分行列元素包括重新获取的四重块Toeplitz矩阵的各不同元素；
利用随机QR降秩分解算法对重新获取的部分行列元素表示的重新获取的四重块Toeplitz矩阵进行降秩处理，以对重新获取的四重块Toeplitz矩阵进行降秩处理；
利用无展开求平均算法对降秩处理后的重新获取的四重块Toeplitz矩阵的对角线元素求平均，得到重新获取的对角线元素平均；
计算重新获取的对角线元素平均和所述对角线元素求平均的结果的差值并判断该差值是否在所述设定误差范围内，若是，将重新获取的对角线元素平均作为所述地震频率切片数据的重建数据，若否，利用重新获取的对角线元素平均再次重新获取所述地震频率切片数据的四重块Toeplitz矩阵，并存储再次重新获取的四重块Toeplitz矩阵的部分行列元素以表示再次重新获取的四重块Toeplitz矩阵，再次重新获取的部分行列元素包括再次重新获取的四重块Toeplitz矩阵的各不同元素，利用随机QR降秩分解算法对再次重新获取的部分行列元素表示的再次重新获取的四重块Toeplitz矩阵进行降秩处理，以对再次重新获取的四重块Toeplitz矩阵进行降秩处理，利用无展开求平均算法对降秩处理后的再次重新获取的四重块Toeplitz矩阵的对角线元素求平均，计算再次重新获取的对角线元素平均和所述重新获取的对角线元素平均的差值，依次迭代进行，直到所得差值在所述设定误差范围内。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611182694.0A CN106646612B (zh) | 2016-12-20 | 2016-12-20 | 基于矩阵降秩的地震数据重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611182694.0A CN106646612B (zh) | 2016-12-20 | 2016-12-20 | 基于矩阵降秩的地震数据重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106646612A CN106646612A (zh) | 2017-05-10 |
| CN106646612B true CN106646612B (zh) | 2018-11-30 |
Family
ID=58834946
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201611182694.0A Expired - Fee Related CN106646612B (zh) | 2016-12-20 | 2016-12-20 | 基于矩阵降秩的地震数据重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106646612B (zh) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109143341A (zh) * | 2018-08-20 | 2019-01-04 | 成都理工大学 | 基于Hampel范数的降秩滤波方法 |
| CN109164483B (zh) * | 2018-08-29 | 2020-04-03 | 中国科学院地球化学研究所 | 多分量地震数据矢量去噪方法及多分量地震数据矢量去噪装置 |
| CN111337972B (zh) * | 2018-12-19 | 2022-08-30 | 中国石油天然气集团有限公司 | 用于地震数据重建的规则化域选取方法、系统 |
| CN111025385B (zh) * | 2019-11-26 | 2020-11-27 | 中国地质大学（武汉） | 一种基于低秩和稀疏约束的地震数据重建方法 |
| CN111830560B (zh) * | 2020-07-24 | 2022-06-28 | 河北工业大学 | 一种基于降秩算法的地震数据重建方法 |
| CN112861074B (zh) * | 2021-03-09 | 2022-10-04 | 东北电力大学 | 基于Hankel-DMD的电力系统机电参数提取方法 |
| CN113009560B (zh) * | 2021-03-23 | 2022-03-29 | 中国地质大学（武汉） | 一种地震数据重建方法、装置、设备及存储介质 |
| CN113094648B (zh) * | 2021-04-02 | 2022-08-09 | 算筹(深圳)信息科技有限公司 | 外积累加求解三角矩阵与矩阵内积的方法 |
| CN114509805B (zh) * | 2021-05-14 | 2023-02-28 | 中国地质大学（北京） | 矢量凸集投影多分量三维地震数据重建方法及装置 |
| CN114441174B (zh) * | 2022-02-09 | 2024-07-09 | 上海电气集团股份有限公司 | 滚动轴承复合故障的诊断方法、系统、设备及介质 |
| CN114966861B (zh) * | 2022-05-17 | 2024-03-26 | 成都理工大学 | 基于Lp伪范数和γ范数稀疏低秩约束的地震去噪方法 |
| CN115598702B (zh) * | 2022-10-25 | 2023-11-28 | 中国矿业大学(北京) | 一种地热资源热储空间构造分布的探测方法及装置 |
| CN115657125B (zh) * | 2022-11-09 | 2024-03-26 | 河北地质大学 | 地震数据重建与去噪的方法、装置和电子设备 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102854533B (zh) * | 2011-07-01 | 2016-02-03 | 中国石油化工股份有限公司 | 一种基于波场分离原理提高地震资料信噪比的去噪方法 |
| CN104459770B (zh) * | 2013-09-24 | 2017-06-16 | 中国石油化工股份有限公司 | 一种高维地震数据规则化方法 |
| CN105549078B (zh) * | 2015-12-31 | 2019-06-11 | 中国石油天然气股份有限公司 | 不规则地震数据的五维插值处理方法及装置 |
| CN105425301B (zh) * | 2016-01-08 | 2018-10-19 | 东华理工大学 | 一种频率域三维不规则地震数据重建方法 |
-
2016
- 2016-12-20 CN CN201611182694.0A patent/CN106646612B/zh not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN106646612A (zh) | 2017-05-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106646612B (zh) | 基于矩阵降秩的地震数据重建方法 | |
| CN106772583B (zh) | 一种地震绕射波分离方法和装置 | |
| Kaur et al. | Seismic data interpolation using deep learning with generative adversarial networks | |
| CN107728211B (zh) | 基于张量核范数正则化的地震信号恢复算法 | |
| Zu et al. | Hybrid-sparsity constrained dictionary learning for iterative deblending of extremely noisy simultaneous-source data | |
| CN110568486B (zh) | 基于同步稀疏低秩张量补全模型的地震信号补全方法 | |
| CN109001802B (zh) | 基于Hankel张量分解的地震信号重构方法 | |
| AU2023202187B2 (en) | Use nuos technology to acquire optimized 2d data | |
| Cheng et al. | Fast dual-domain reduced-rank algorithm for 3D deblending via randomized QR decomposition | |
| Wang et al. | Fast dictionary learning for high-dimensional seismic reconstruction | |
| CN105425301A (zh) | 一种频率域三维不规则地震数据重建方法 | |
| GB2460175A (en) | Interpolating seismic data by projection onto convex sets | |
| CN105549078B (zh) | 不规则地震数据的五维插值处理方法及装置 | |
| Chen et al. | Five-dimensional seismic data reconstruction using the optimally damped rank-reduction method | |
| CN104459770B (zh) | 一种高维地震数据规则化方法 | |
| Lu et al. | A fast SVD for multilevel block Hankel matrices with minimal memory storage | |
| Jumah et al. | Dimensionality‐reduced estimation of primaries by sparse inversion | |
| Zhang et al. | 2-D seismic data reconstruction via truncated nuclear norm regularization | |
| WO2010014118A1 (en) | Statistical decoding and imaging of multishot and single-shot seismic data | |
| Cavalcante et al. | Low‐rank seismic data reconstruction and denoising by CUR matrix decompositions | |
| Cheng et al. | A fast rank-reduction algorithm for 3D deblending via randomized QR decomposition | |
| CN114325821A (zh) | 基于3d-snacnn网络的叠前地震资料中强散射噪声压制方法及系统 | |
| Trad | Five dimensional seismic data interpolation | |
| CN109782346B (zh) | 一种基于形态成分分析的采集脚印压制方法 | |
| Gou et al. | Complex seismic wavefield interpolation based on the Bregman iteration method in the sparse transform domain |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20181130 Termination date: 20191220 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |