CN111929733A - 基于slice采样的地震信号正则化处理方法 - Google Patents
基于slice采样的地震信号正则化处理方法 Download PDFInfo
- Publication number
- CN111929733A CN111929733A CN202010847903.9A CN202010847903A CN111929733A CN 111929733 A CN111929733 A CN 111929733A CN 202010847903 A CN202010847903 A CN 202010847903A CN 111929733 A CN111929733 A CN 111929733A
- Authority
- CN
- China
- Prior art keywords
- data
- tensor
- sampling
- slice
- 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.)
- Pending
Links
- 238000005070 sampling Methods 0.000 title claims abstract description 35
- 238000003672 processing method Methods 0.000 title claims abstract description 13
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 81
- 238000000034 method Methods 0.000 claims abstract description 49
- 239000000126 substance Substances 0.000 claims description 7
- 230000008602 contraction Effects 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims description 2
- 238000012545 processing Methods 0.000 abstract description 28
- 238000012217 deletion Methods 0.000 abstract description 14
- 230000037430 deletion Effects 0.000 abstract description 14
- 238000004458 analytical method Methods 0.000 abstract description 10
- 238000004088 simulation Methods 0.000 abstract description 3
- 238000011084 recovery Methods 0.000 description 20
- 239000011159 matrix material Substances 0.000 description 13
- 230000008569 process Effects 0.000 description 11
- 230000009466 transformation Effects 0.000 description 11
- 238000000354 decomposition reaction Methods 0.000 description 10
- 230000006870 function Effects 0.000 description 9
- 230000000694 effects Effects 0.000 description 8
- 238000010586 diagram Methods 0.000 description 6
- 239000007789 gas Substances 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 5
- 239000000047 product Substances 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 230000009467 reduction Effects 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- VPNMENGBDCZKOE-LNYNQXPHSA-N [(1R,5S)-8-methyl-8-azabicyclo[3.2.1]octan-3-yl] 3-hydroxy-2-phenylpropanoate [(1S,2S,4R,5R)-9-methyl-3-oxa-9-azatricyclo[3.3.1.02,4]nonan-7-yl] (2S)-3-hydroxy-2-phenylpropanoate (1R,3S)-1,2,2-trimethylcyclopentane-1,3-dicarboxylic acid Chemical compound CC1(C)[C@H](CC[C@@]1(C)C(O)=O)C(O)=O.CC1(C)[C@H](CC[C@@]1(C)C(O)=O)C(O)=O.CN1[C@H]2CC[C@@H]1CC(C2)OC(=O)C(CO)c1ccccc1.CN1[C@H]2CC(C[C@@H]1[C@H]1O[C@@H]21)OC(=O)[C@H](CO)c1ccccc1 VPNMENGBDCZKOE-LNYNQXPHSA-N 0.000 description 3
- 238000013016 damping Methods 0.000 description 3
- 230000004069 differentiation Effects 0.000 description 3
- 230000003993 interaction Effects 0.000 description 3
- 239000000203 mixture Substances 0.000 description 3
- 239000004065 semiconductor Substances 0.000 description 3
- 238000013473 artificial intelligence Methods 0.000 description 2
- 239000003086 colorant Substances 0.000 description 2
- 238000007405 data analysis Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000003909 pattern recognition Methods 0.000 description 2
- 238000010183 spectrum analysis Methods 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000005183 dynamical system Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000001939 inductive effect Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 229920000642 polymer Polymers 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000012916 structural analysis Methods 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 238000011179 visual inspection Methods 0.000 description 1
Images
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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
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
本发明公开一种基于slice采样的地震信号正则化处理方法，应用于地震勘探领域，针对VSP数据中的面缺失问题，本发明提出低秩稀疏张量补全(Low‑tubal‑rank and Sparse Tensor Completion，简称LRSTC)算法来完成信号正则化处理；利用增强拉格朗日方法来求解LRSTC模型；通过在合成地震数据和实际地震数据上的实验仿真及分析，验证了本发明的LRSTC的卓越性能。
Description
技术领域
本发明属于地震勘探领域，特别涉及一种地震数据分析技术。
背景技术
进入21世纪以来，科技发展迅速，油气勘探技术也愈发成熟。经过多年的开发，在一些石油工业区勘探到的石油储量已经明显减少。为了满足对能源的迫切需求，必须加大对石油、天然气和页岩气的勘探力度，进一步提高地震勘探技术，寻找更加复杂和隐蔽性更高的断裂带、地下岩层构造区域，这就对地震勘探提出了更高的要求。
在对能源的强烈需求下，利用更加先进的地震勘测方法探测油气田已经变的越来越重要。地震勘测工作主要包括以下三个主要环节：地震数据采集、地震数据处理和地震数据解释。地震数据采集是地震勘测工作的基础，地震数据处理是地震勘测工作的重点环节，地震数据解释是地震勘测工作的具体应用，三个环节相辅相成，缺一不可。地震数据采集就是通过专业数字地震仪器按照预先设计好的数据采集方法在野外采集地震波数据。具体操作为利用可控震源车激发产生的地震波，信号接收器接收到从地下返回的地震波信号。通过接收器接收到的地震波信号，反演出地下的地质构造。地震数据处理是利用电子计算机对所采集到的地震波数据进行各种数字处理，例如正则化处理、反褶积处理等。地震数据处理是以提高采集到的地震波数据质量为目的，对地下构造和地质体进行成像，以便进行后续地震数据解释。地震数据解释是对地震数据处理后的各种数据体进行构造分析和岩性解释，以发现构造、地层和岩性圈闭各类油气藏，增加油气储量、提高油气产量。地震数据采集、地震数据处理和地震数据解释三个环节是相互紧密结合在一起的。只有结合勘探区域的地质特点和具体勘探的地质任务，选择合适的地震数据处理方法，才能得到高质量的处理结果。
理想的地震勘测应该涵盖所有可能的发射器和接收器组合。但由于经济和物理条件上的限制，这些情况很少出现。特别是，在采集数据之前已经计划好发射器和接收器的位置，但是真实的采样位置通常会受到地表障碍物(城市、村庄、桥梁等建筑物)、地形条件(高山、湖泊、沼泽等自然环境)的限制，偏离原始目标位置，从而导致不太理想的中点和方位角。此外，由于地震采集几何结构的固有性质，中点偏移方位角域也是不完全采样。这些影响在土地调查中更为明显。对于依赖于网格数据的方法，欠采样的程度取决于对数据进行分箱的方式。还应注意，即使在完全采集和给定几何形状的情况下，在执行其他过程的领域中，发射器和接收器数量稀少、放置不均匀等情况也会导致收集到的地震数据仍可能欠采样。这种空间欠采样会对诸如迁移、反演和定量解释研究(例如幅度与角度分析)之类的后续处理步骤产生不利影响。因此，在地震数据处理中，对欠采样得到的地震信号进行正则化处理就显得尤为重要。图1为实际工区地震勘测的示意图，可以明显看出检波器的布置并不规则。地震信号正则化处理就是将缺失的地震数据补充完整并且压制异常噪声，以增强地震数据成像的质量。图2为地震数据正则化处理的示意图。
多维地震信号可以借助张量来表征，由于地震数据时间和空间上的强相关性，可以对采集到的地震张量施加秩约束进行正则化处理。常见的正则化处理有张量核范数正则化处理。依据核范数正则化对地震数据进行处理，大体又可分为两种：一种是基于奇异值分解的多通道奇异谱分析方法(Multichannel Singular Spectrum Analysis，简称MSSA)或Cadzow重建方法。另一种是不采用奇异值分解计算的地震信号重构方法。
核范数正则化处理最近几年才在地震数据中广泛使用。核范数正则化处理的前提是地震数据是低秩结构，缺少数据和噪声污染会增加数据的秩。矩阵秩降低方法[Trickett.S,Burroughs.L,Milton.A,et al.Rank-reduction-based traceinterpolation[C].SEG Technical Program Expanded Abstracts 2010.2010.]介绍了一种在恒定频率切片上使用基于SVD的截断矩阵秩降低算法。MSSA算法同时用于处理地震数据重构和降噪问题，该算法将原始空间地震数据映射到块Hankel矩阵中，然后使用加权凸集投影类(Projection On Convex Sets，简称POCS)算法和随机奇异值分解(RandomSingular Value Decomposition，简称RSVD)以快速最小化目标块Hankel矩阵的秩。MSSA算法在截断奇异值分解(Truncated Singular Value Decomposition，简称T-SVD)过程中无法完全抑制掉噪声，提出通过调节引入的阻尼因子参数来控制的阻尼算子，以进一步衰减常规MSSA框架中的残留噪声，因此阻尼多通道奇异谱分析(Damped MultichannelSingular Spectrum Analysis，简称DMSSA)[Zhang.D,Chen.Y,Huang.W,et al.Multi-stepdamped multichannel singular spectrum analysis for simultaneousreconstruction and denoising of 3D seismic data[J].Journal of Geophysics andEngineering,2016,13(5):704-720.]算法应运而生，同时还提出了多步阻尼MSSA进一步提高性能。混合秩-稀疏约束(Hybrid Rank-Sparsity Constraint，简称HRSC)[Zhang.D,Zhou.Y,Chen.H,et al.Hybrid rank-sparsity constraint model for simultaneousreconstruction and denoising of 3D seismic data[J].Geophysics,2017,82(5):V351-V367.]是在MSSA的基础上进行扩展的方法，是一种数据驱动的技术，在低秩约束的基础上加入稀疏约束，在保证恢复精度的同时也提升了信噪比。Kreimer和Sacchi首次引入高阶奇异值分解(Higher-Order Singular Value Decomposition，简称HOSVD)来重构5D多线性张量或数组。叠前数据是可压缩的，因为它们具有较低的信息复杂性，数据的空间相关性和时间相关性比较强烈，表示为tSVD域中的低秩张量核范数(Tensor Nuclear Norm，简称TNN)[Kreimer.N,Stanton.A,Sacchi.M.D.Tensor completion based on nuclear normminimization for 5D seismic data reconstruction[J].Geophysics,2013,78(6):V273-V284.]。这种可压缩性源自t-product的卷积结构，这是由于基于tSVD的表示在描述数据之间的移位和缩放操作方面特别有效。因此，TNN算法可以在有限的采样下使用较少的复杂度惩罚算法可靠地恢复数据。截断核范数正则化(Truncated Nuclear NormRegularization，简称TNNR)算法[Zhang.D,Hu.Y,Ye.J,et al.Matrix completion bytruncated nuclear norm regularization[C].2012IEEE Conference on ComputerVision and Pattern Recognition.IEEE,2012:2192-2199.]提出利用张量截断核范数，让数据更好匹配低秩模型以完成数据正则化。权值核范数最小化(Weight Nuclear NormMinimization，简称WNNM)算法[Gu.S,Zhang.L,Zuo.W,et al.Weighted nuclearnormminimization with application to image denoising[C].Proceedings of theIEEE conference on computer vision and pattern recognition.2014:2862-2869.]提出带有权值的核范数，对每个分解得到的奇异值分配不同的权重，进一步完善和扩展了TNN算法。
2012年，Kumar等[Kumar.R,Aravkin.A,Herrmann F.Fast methods for rankminimization with applications to seismic-data interpolation[C].SEG expandedabstracts:82nd Annual international meeting.2012:1-5.]提出了一种无SVD的方法，该方法使用能够避免每次迭代中的T-SVD，这消耗了大多数降秩算法的主要计算成本。Gao等[Gao.J,Sacchi.M.D,Chen X.A fast reduced-rank interpolation method forprestack seismic volumes that depend on four spatial dimensions[J].Geophysics,2013,78(1):V21-V30.]提出了快速叠前地震数据内插法，通过将5D地震数据嵌入到块Toeplitz/Hankel矩阵，然后使用Lanczos双角化方法可快速降低秩。通过对比，视频数据和地震数据具有很高的相似性，在时间和空间方向具有很强的冗余性，这也对地震信号重构提供了思考。Tubal-Alt-Min算法[Liu.X.Y,Aeron.S,Aggarwal V,et al.Low-tubal-rank tensor completion using alternating minimization[C].Modeling andSimulation for Defense Systems and Applications XI.International Society forOptics and Photonics,2016,9848:984809.]将求解的low-tubal-rank张量分解为两个小因子张量的乘积，使用最小二乘最小化方法交替计算出两个因子张量，最后重构出求解张量。
发明内容
为解决上述技术问题，本发明提出一种基于slice采样的地震信号正则化处理方法，采用低秩稀疏张量补全(Low-tubal-rank and Sparse Tensor Completion，简称LRSTC)算法来完成信号正则化处理；具体的：将基于Hankel变换的张量补全引入VSP数据补全问题；通过Hankel变换构建的桥梁，将slice采样问题转换为tubal采样下的三阶张量补全问题，简化建模难度。
本发明采用的技术方案为：一种基于slice采样的地震信号正则化处理方法，包括：
S1、VSP数据经Hankel结构变换后，重构问题模型建模为：
其中，
表示字典，
表示对应的系数张量，TΩ表示经过Hankel变换的采样算子，T表示Hankel张量，
R表示实数，nt表示时间上的维数，no表示偏移距上的维数，nx表示xline上的维数；
S2、对
上施加稀疏约束，得到具有低秩和稀疏性的扩展张量补全模型；
S3、采用ADMM方法对步骤S2得到的模型进行求解。
进一步地，步骤S2所述扩展张量补全模型表达式为：
其中，
表示加权核范数，‖‖1表示1范数，
表示张量的Frobenius范数，ε表示缺失数据，
表示采样算子。
更进一步地，步骤S3还包括：将扩展张量补全模型转化为增强拉格朗日形式，增强拉格朗日式子表示如下：
其中，μ1、μ2是惩罚因子，
是拉格朗日因子，<·，·>表示两个张量的内积。
进一步地，所述扩展张量补全模型的增强拉格朗日形式通过迭代ALM-ADM方法求解，其解为：
参数
和ε通过最小化增广拉格朗日式子，利用其导数为零求解得到。
更进一步地，
的求解表达式为：
通过迭代收缩阈算法求解得到
进一步地，
的求解表达式为：
通过张量的奇异值收缩算法求解得到
更进一步地，
的更新式子为：
其中，
是
通过奇异值分解的奇异值向量，
上标
表示转置。
本发明的有益效果：本发明将基于Hankel变换的张量补全引入VSP数据正则化处理问题中；在相邻的完整工区中训练出完备的字典张量，将其作为先验条件引入到目标函数中，并通过Hankel变换将slice缺失问题转换为tubal缺失问题；然后利用增强拉格朗日方法来求解LRSTC模型；最后使用合成数据和实际数据进行实验，可视化处理后的数据，证明了数据重建方面LRSTC算法性能。
附图说明
图1为实际工区地震勘测示意图；
图2为地震数据正则化处理二维示意图；
图3为本发明的方法流程图；
图4为3D空间的随机前向切片采样示意图；
图5为六种算法恢复结果示意图；
图6为三种算法迭代收敛图；
图7为三种算法在不同采样率下恢复效果图；
图8为理论VSP数据第8个截面恢复结果示意图；
图9为实际VSP数据第12个截面恢复结果示意图。
具体实施方式
为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。
本发明的目的是提出低秩稀疏张量补全(Low-tubal-rank and Sparse TensorCompletion，简称LRSTC)算法来完成信号正则化处理。具体来讲，将基于Hankel变换的张量补全引入VSP数据补全问题。通过Hankel变换构建的桥梁，将slice采样问题转换为tubal采样下的三阶张量补全问题，简化建模难度。
VSP数据分析是对井中地震数据进行分析。地表上障碍物的存在和经济条件的限制会导致地震数据记录的丢失，因此VSP信号正则化处理就显得尤为重要。VSP信号正则化就是从采集到的VSP数据集中恢复丢失的地震记录。VSP信号正则化处理和地表地震信号正则化处理了主要的不同是采样方式的不同，VSP缺失表现为3D张量空间上的时间-偏移距切片缺失。
故本发明着眼于VSP数据中的面缺失问题，提出低秩稀疏张量补全(Low-tubal-rank and Sparse Tensor Completion，简称LRSTC)算法来完成信号正则化处理。利用增强拉格朗日方法来求解LRSTC模型。通过在合成地震数据和实际地震数据上的实验仿真及分析，LRSTC的卓越性能都得到了证明。
如图3所示，本发明的实现过程包括：
1、3D VSP数据重建的稀疏-低秩张量补全模型
本节介绍基于3D VSP数据重建的稀疏-低秩(Sparse and Low-tubal-rank，简称SLR)张量补全模型。与地震数据补全不同，我们的模型将VSP数据缺失视为图4所示的3D空间的随机前向切片采样。
张量空间下，完整的VSP数据可以表示为
nt，no和nx分别表示时间、偏移距和xline上的维数。在xline模式下，用一个前向切片集对VSP数据进行采样。
表示采样算子，在观测集
上的操作可以用下式表达出：
公式(1)表示如果k是观测集Ω里的数据，则进行采样操作，否则不采样。
在噪声条件下，观测到的VSP数据
可以由下式表示：
由于VSP数据中几乎不存在噪声，在这里我们并不对噪声数据进行处理，N表示噪声，
由于观测到的数据数量明显少于
中的元素数量，因此反问题(2)的病态严重，如果不对
施加额外约束，则无法直接解决。考虑到一种特殊情况，其中一些连续切片中的所有元素都在VSP张量数据中丢失。在这种情况下，张量核范数正则化方法通常无法恢复丢失的元素[Yokota.T.Missing Slice Recovery for Tensors Using a Low-rank Modelin Embedded Space[J].Proceedings/CVPR,IEEE Computer Society Conference onComputer Vision and Pattern Recognition.IEEE Computer Society Conference onComputer Vision and Pattern Recognition,2018.]。然而这种情况在VSP数据采集中经常发生。因此，张量核范数最小化无法提供足够的正则化约束来恢复方程式(3)中的不完整张量。
我们考虑一种情况，其中VSP数据中一些连续切片中的所有元素都丢失。在这种情况下，张量核范数正则化方法通常无法恢复丢失的元素[Yokota.T.Missing SliceRecovery for Tensors Using a Low-rank Model in Embedded Space[J].Proceedings/CVPR,IEEE Computer Society Conference on Computer Vision and PatternRecognition.IEEE Computer Society Conference on Computer Vision and PatternRecognition,2018.]。由于地面障碍物的存在，这种情况在VSP数据中非常常见。因此，核规范最小化无法提供足够的正则化以恢复方程式(3)中的不完全张量。
为了解决上述问题，本发明引入了2倍Hankel矩阵，以得到Hankel张量。形式如下：
其中
表示Hankel张量。经过Hankel变换后，切片缺失问题变成了更容易解决的tubal缺失问题，如下所示：
问题(3)就会进一步变成如下问题：
2、目标函数构建
意识到解决问题(6)的难点在于对隐张量的正则化约束不足，我们建议使用合适的先验条件进一步对张量的列进行正则化，以约束多维信号中的每个前向切片的值。许多先验可用于增强经典张量补全模型，例如总变化量，分段平滑度和某些基础(字典)下的稀疏度。但是，这些先验不能很好地解决我们的问题。Hankel结构张量在数据恢复求解过程中扮演重要角色，同样在解决VSP数据重构问题，本发明将Hankel结构引入其中。经过与slice缺失方向相交的方向上的Hankel变换，slice缺失问题转换为tubal缺失问题。经过可行性分析，发现效果并不好，这是因为转换后的tubal缺失问题为连续整块tubal缺失，并不是简单的随机tubal缺失。在处理整块数据缺失问题上，STC框架并不能很好地解决问题。
本发明选择已经学习过的高质量的张量字典。由于局部区域数据在2D VSP数据中是完整的，可以通过学习这个完整区域的字典，用这个高效字典去对邻近区域的不完整数据进行重构。字典的学习方式可以通过如下式子来表示。
本发明利用局部的区分性信息以恢复全局中丢失的信息。VSP数据
经过Hankel结构变换后，将重构问题建模为：
其中，
表示字典，
表示对应的系数张量，TΩ表示经过Hankel变换的采样算子，‖‖1表示1范数。
的每一列都具有稀疏表示，可以用
来表示。通过对
上施加稀疏约束，得到以下具有低秩和稀疏性的扩展张量补全模型。
关于稀疏表示[Yang.J,Yang.X,Ye X,et al.Reconstruction of Structurally-Incomplete Matrices With Reweighted Low-Rank and Sparsity Priors[J].IEEETransactions on Image Processing,2017,26(3):1158-1172.]和矩阵重构[Wang,Benfeng.An Efficient POCS Interpolation Method in the Frequency-Space Domain[J].IEEE Geoscience and Remote Sensing Letters,2016,13(9):1384-1387.]的最新工作表明，重新加权的先验分别显著地提高了稀疏性和低秩性。我们所提出的SLR模型是受等式约束的最小化问题。本发明将张量补全问题作为张量恢复问题的特例，即将丢失的行列设置为零，SLR模型可以重新表示为：
其中，ε表示缺失数据。
3、目标函数求解
求解问题(9)主要用到的方法仍然是ADMM方法。通过ADMM方法的思想，利用增强拉格朗日思路将约束项放入到目标函数当中。问题(9)的增强拉格朗日式子表示如下：
其中μ1和μ2是惩罚因子，
和
是拉格朗日因子，<·，·>表示两个张量的内积，
表示张量的Frobenius范数。
问题(10)通过迭代ALM-ADM方法求解，其解为：
其中，ρ1＞1和ρ2＞1用来确保参数μ1和μ2是正增长序列。我们基于上述方案开发了一种迭代交替算法来求解，该算法在算法1中进行了概述。
方程(11)其中参数
和ε可以通过最小化增广拉格朗日式子，利用其导数为零求解得到，其具体形式如下：
1)更新
该问题等效于优化求解以下问题：
上述最小化问题可以通过迭代收缩阈算法来求解[Yang.J,Yang.X,Ye X,etal.Reconstruction of Structurally-Incomplete Matrices With Reweighted Low-Rank and Sparsity Priors[J].IEEE Transactions on Image Processing,2017,26(3):1158-1172.]、[Jiang.F,Liu.X.Y,Lu.H,et al.Efficient multi-dimensional tensorsparse coding using t-linear combination[C].Thirty-Second AAAI Conference onArtificial Intelligence.2018.]。为了方便表述，我们用
上式就可以表示如下：
通过引入近似变量
问题(4-13)又可以写成如下方式：
其中Lf+1是莱布尼茨常数，
是在张量空间内
的梯度。因此方程(12)可以表示为：
为了解决上述问题，我们首先计算出
表示为
的梯度，表达如下：
接下来我们来讨论如何确定上式当中莱布尼茨常数的问题，对于
我们可以得到如下式子：
我们算法当中
的莱布尼茨常数可以表达为
为了简便，我们定义
的l1范数最小可以通过在
上的软阈值操作得到。具体表达如下：
其中soft(·，·）表示软阈值操作，其数学表达式为
其中近似操作符
按照以下策略更新。
一旦算法达到收敛条件，就得到稳定的
算法2概括了这部分过程。
2)更新
使用同样的方式求解
其问题等价于求解如下式子：
其中
问题(20)的最小化可以通过张量的奇异值收缩算法(Tensor Singular Value Thresholding，简称TSVT)求解，具体更新式子为：
其中
是
通过奇异值分解的奇异值向量，有如下关系：
算法3概括了这部分过程。
3)更新ε：
我们需要更新观测域
内和数据缺失域Ω中的数值。因此ε更新分为两个部分(即
与Ω两个部分)，ε的一阶最优条件为：
因而，我们可以得到：
经过以上的详细步骤，即可将问题(9)求解出来。在算法开始迭代时，需要初始化参数
可以通过方程(24)求解得到。
初始值
是通过将
的所有元素值都设置成
所有元素的平均值得到的。与[Yang.J,Yang.X,Ye X,et al.Reconstruction of Structurally-Incomplete MatricesWith Reweighted Low-Rank and Sparsity Priors[J].IEEE Transactions on ImageProcessing,2017,26(3):1158-1172.]、[Jiang.F,Liu.X.Y,Lu.H,et al.Efficientmulti-dimensional tensor sparse coding using t-linear combination[C].Thirty-Second AAAI Conference on Artificial Intelligence.2018.]相同，该问题可以通过迭代收缩阈值算法基于张量积(ISTA-T)求解得到。
通过将其应用于缺失3D VSP地震数据重构中，验证所提出方法的可行性和优越性。考虑两种不同类型的3D地震数据：1.理论模型数据；2.实际地震数据。在第本节中，我们首先描述了整个过程的实验设置。然后，将提出的算法在理论模型数据和实际数据上进行验证，并与常用的重构算法进行比较。为了展现模型的优越性，我们采用如下衡量指标来验证。
4、本发明方法的验证：
为了进行比较，我们选取了五种比较算法，以揭示LRSTC算法的优越性。将几种算法总结如下：
·KTSVD[Zhang.Z,Aeron.S.Denoising and Completion of 3D Data viaMultidimensional Dictionary Learning[J].2015.]：将用于1D数据的K-SVD算法扩展到用于处理2D和3D数据的KTSVD算法
·TNN[Ely.G,Aeron.S,Hao.N,et al.5D seismic data completion anddenoising using a novel class of tensor decompositions[J].Geophysics,2015,80(4):V83-V95.]：通过使用tSVD分解，对提出的tubal秩实行最小化算法对不完整数据进行恢复和降噪。
·DMSSA[Zhang.D,Chen.Y,Huang.W,et al.Multi-step damped multichannelsingular spectrum analysis for simultaneous reconstruction and denoising of3D seismic data[J].Journal of Geophysics and Engineering,2016,13(5):704-720.]、[Li.X.P,Huang.L,So.H.C,et al.A Survey on Matrix Completion:Perspectiveof Signal Processing[J].arXiv preprint arXiv:1901.10885,2019.]、[Rabanser.S,Shchur.O,Günnemann.S.Introduction to tensor decompositions and theirapplications in machine learning[J].arXiv preprint arXiv:1711.10781,2017.]：使用阻尼DMSSA算法可以更好地将数据分解为信号子空间和噪声子空间，以进行随机噪声衰减。
·PMF[Gao.J,Cheng.J,Sacchi M D.Five-Dimensional SeismicReconstruction Using Parallel Square Matrix Factorization[J].IEEETransactions on Geoscience and Remote Sensing,2016,PP(99):1-12.]：通过对张量进行不同的展开方式展开，用矩阵去表征张量。利用低秩矩阵分解的方式来恢复缺少数据的张量。
·MDT[Yokota.T.Missing Slice Recovery for Tensors Using a Low-rankModel in Embedded Space[J].Proceedings/CVPR,IEEE Computer Society Conferenceon Computer Vision and Pattern Recognition.IEEE Computer Society Conferenceon Computer Vision and Pattern Recognition,2018.]：通过从动力学系统的研究中引入“延迟嵌入”的概念，可以通过考虑嵌入式空间中的低秩模型来恢复丢失的切片。
为了检测所提出算法的有效性，将LRSTC算法应用于VSP理论数据中。对于测试数据，数据集为Nx×No×Nt的3D VSP数据，具体大小为Nx＝31，No＝30和Nt＝179。在具体的实验中，我们不考虑噪声问题。
图5展示了在连续三个界面缺失的情况下，各个算法的恢复结果示意图。从图5可以很明显的看出，TNN、KTSVD和PMF算法没有重构出数据，而MDT、DMSSA和LRSTC三种算法成功重建出数据。这是因为后续这三种算法在与slice缺失方向相交的方向上做Hankel变换，将slice缺失问题转换为tubal缺失问题，打乱了数据原有的缺失形式。同时也说明一个问题，对slice缺失的地震数据施加核范约束，并不能很好地重构出数据。下面我们将对这三种算法做一个更加细致的分析。
为了完整起见，用理论数据来演示LRSTC算法的收敛性能。为了展现模型的卓越性，本章采用两个指标来衡量恢复结果的质量。给定原始的完全采样数据
和从不完整数据中重构出的数据
两个指标定义如下所示：
相对平方误差(RSE)定义为：
信噪比(SNR)定义为：
在此，我们绘制出基于50％slice采样率的3D VSP理论数据集上不同算法的RSE收敛图。在图6显示了每种算法的RSE误差性能与迭代次数的关系。为了突出显示出现较小数字的区域，图6是半对数图，在迭代数轴上使用线性刻度，在RSE误差轴上使用对数刻度。如图6所示，实验表明，LRSTC算法提供了最快的绝对收敛速度和最低的最终RSE值。LRSTC算法在15次迭代左右就达到收敛值，SNR最终达到9.7757，相比较与其他两种算法，SNR比DMSSA算法提升2dB左右，比MDT算法提升8dB左右。
之后将这三种方法应用于不完整的3D VSP数据进行正则化处理，该数据是通过以不同的slice采样率{10,20,30,40,50,60,70,80,90％}对完整的VSP数据进行随机slice采样而生成的。图7总结了DMSSA、MDT和LRSTC算法相对于不同采样率的恢复性能，LRSTC算法优于其他两种方法。对于采样率从10％到40％的采样数据，LRSTC算法改进效果并不是很明显。但在50％到80％，LSRTC算法的处理结果要明显好于其他两种算法，证明了该算法的能力。这是因为在缺失率高的情况下，经过Hankel变换后，slice缺失问题转换为连续整块tubal缺失问题，字典的先验作用体现出来。丢失90％的3D VSP数据，数据缺失太多影响LRSTC算法中ADMM的求解，三种算法恢复效果都不理想。
图8展示了在无噪情况下对理论VSP地震数据的正则化效果比较图。图8(b)、8(c)和8(d)分别显示了DMSSA、MDT和LRSTC算法的正则化结果示意图。在图8(b)、8(c)和8(d)中的红色矩形显示的区域中观察到LRSTC算法的优点。为了更直观地进行视觉检查，使用地震工具箱SeisLab在图8(f)、8(g)和8(h)中显示重建数据地面真实数据之间波形重合程度的比较。对于整个面缺失的数据，MDT算法恢复出的结果都会与地面真实情况存在明显差异。DMSSA算法虽然恢复出整个面，但是与LRSTC算法的结果相比还是有较大差距。从图8(f)、8(g)和8(h)可以看出，LRSTC算法的优势可以解释为重建数据地面真实数据之间波形的重合。从图8(f)，8(g)和8(h)可以清楚地看到，尽管DMSSA和MDT算法都重建了丢失的数据，但我们的LRSTC算法获得了更好的结果，尤其是在显著提高SNR方面。
以下采用实际工区实验，用于说明本发明方法的技术发效果：
为了在实践中展现所提出方法的优越性能，将LRSTC算法应用于真实的3D VSP数据。对于测试数据，将数据集裁剪为大小为Nx×No×Nt的三个数据张量，此处Nx＝77，No＝51和Nt＝1500。受到实验设备的限制，我们在此数据上截取一部分数据来进行验证。对于每个测试数据，我们都采用截断奇异值分解方法对原始数据进行预处理以抑制原始数据中的噪声，从而收集地面真实数据以适用于我们的实验。图9(a)显示了原始数据。之后，从真实数据中随机删除50％前向切片，以模拟丢失的数据。图9(e)中显示了原始数据的局部区域图。
以图9为例，图9(b)、9(c)、9(d)分别显示对实际VSP数据的重建结果比较图。图9(b)、9(c)、9(d)分别显示了DMSSA、MDT和LRSTC算法的重构结果。从图9(b)、9(c)、9(d)中的红色矩形突出显示的区域中观察到LRSTC算法的优点。同样使用地震工具箱SeisLab中s_compare函数来对比重建数据和地面真实数据之间的差距。对于slice缺失的问题，MDT算法重构效果最差，与实际VSP数据存在很大差异。DMSSA算法恢复效果较理想，但也与真实数据存在差异。从图9(h)中可以看出，LRSTC算法的优势可以解释为重建数据地面真实数据之间的重合程度。
图8与图9中横坐标的Trace number表示记录道编号，纵坐标的Time表示时间。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。
Claims (7)
1.一种基于slice采样的地震信号正则化处理方法，其特征在于，包括：
S1、VSP数据经Hankel结构变换后，重构问题模型建模为：
其中，
表示字典，
表示对应的系数张量，TΩ表示经过Hankel变换的采样算子，T表示Hankel张量，
R表示实数，nt表示时间上的维数，no表示偏移距上的维数，nx表示xline上的维数；
S2、对
上施加稀疏约束，得到具有低秩和稀疏性的扩展张量补全模型；
S3、采用ADMM方法对步骤S2得到的模型进行求解。
2.根据权利要求1所述的一种基于slice采样的地震信号正则化处理方法，其特征在于，步骤S2所述扩展张量补全模型表达式为：
其中，
表示加权核范数，‖‖1表示1范数，
表示张量的Frobenius范数，ε表示缺失数据，
表示采样算子。
3.根据权利要求1所述的一种基于slice采样的地震信号正则化处理方法，其特征在于，步骤S3还包括：将扩展张量补全模型转化为增强拉格朗日形式，增强拉格朗日式子表示如下：
其中，μ1、μ2是惩罚因子，
是拉格朗日因子，<·，·>表示两个张量的内积。
4.根据权利要求1所述的一种基于slice采样的地震信号正则化处理方法，其特征在于，所述扩展张量补全模型的增强拉格朗日形式通过迭代ALM-ADM方法求解，其解为：
参数
和ε通过最小化增广拉格朗日式子，利用其导数为零求解得到。
5.根据权利要求4所述的一种基于slice采样的地震信号正则化处理方法，其特征在于，
的求解表达式为：
通过迭代收缩阈算法求解得到
6.根据权利要求4所述的一种基于slice采样的地震信号正则化处理方法，其特征在于，
的求解表达式为：
通过张量的奇异值收缩算法求解得到
7.根据权利要求6所述的一种基于slice采样的地震信号正则化处理方法，其特征在于，
的更新式子为：
其中，
是
通过奇异值分解的奇异值向量，
上标
表示转置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010847903.9A CN111929733A (zh) | 2020-08-21 | 2020-08-21 | 基于slice采样的地震信号正则化处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010847903.9A CN111929733A (zh) | 2020-08-21 | 2020-08-21 | 基于slice采样的地震信号正则化处理方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN111929733A true CN111929733A (zh) | 2020-11-13 |
Family
ID=73304966
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010847903.9A Pending CN111929733A (zh) | 2020-08-21 | 2020-08-21 | 基于slice采样的地震信号正则化处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111929733A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112988760A (zh) * | 2021-04-27 | 2021-06-18 | 北京航空航天大学 | 一种基于张量分解的交通时空大数据缺失补全方法 |
| CN113655534A (zh) * | 2021-07-14 | 2021-11-16 | 中国地质大学（武汉） | 基于多线性奇异值张量分解核磁共振fid信号噪声抑制方法 |
| CN114662045A (zh) * | 2022-03-24 | 2022-06-24 | 电子科技大学 | 基于框架集的p阶张量深度学习的多维地震数据去噪方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2781936A2 (en) * | 2013-03-22 | 2014-09-24 | CGG Services SA | System and method for interpolating seismic data |
| CN107561576A (zh) * | 2017-08-31 | 2018-01-09 | 电子科技大学 | 基于字典学习正则化稀疏表示的地震信号恢复方法 |
| CN107728211A (zh) * | 2017-08-31 | 2018-02-23 | 电子科技大学 | 基于张量核范数正则化的地震信号恢复算法 |
| CN109655890A (zh) * | 2017-10-11 | 2019-04-19 | 中国石油化工股份有限公司 | 一种深度域浅中深层联合层析反演速度建模方法及系统 |
| CN110568486A (zh) * | 2019-09-17 | 2019-12-13 | 电子科技大学 | 基于同步稀疏低秩张量补全模型的地震信号补全方法 |
| CN110855773A (zh) * | 2019-11-09 | 2020-02-28 | 北京工业大学 | 一种Web中面向服务环境中的基于张量的信任评估方法 |
-
2020
- 2020-08-21 CN CN202010847903.9A patent/CN111929733A/zh active Pending
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2781936A2 (en) * | 2013-03-22 | 2014-09-24 | CGG Services SA | System and method for interpolating seismic data |
| CN107561576A (zh) * | 2017-08-31 | 2018-01-09 | 电子科技大学 | 基于字典学习正则化稀疏表示的地震信号恢复方法 |
| CN107728211A (zh) * | 2017-08-31 | 2018-02-23 | 电子科技大学 | 基于张量核范数正则化的地震信号恢复算法 |
| CN109655890A (zh) * | 2017-10-11 | 2019-04-19 | 中国石油化工股份有限公司 | 一种深度域浅中深层联合层析反演速度建模方法及系统 |
| CN110568486A (zh) * | 2019-09-17 | 2019-12-13 | 电子科技大学 | 基于同步稀疏低秩张量补全模型的地震信号补全方法 |
| CN110855773A (zh) * | 2019-11-09 | 2020-02-28 | 北京工业大学 | 一种Web中面向服务环境中的基于张量的信任评估方法 |
Non-Patent Citations (1)
| Title |
|---|
| 张仓仓: "多维地震信号正则化处理方法研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112988760A (zh) * | 2021-04-27 | 2021-06-18 | 北京航空航天大学 | 一种基于张量分解的交通时空大数据缺失补全方法 |
| CN113655534A (zh) * | 2021-07-14 | 2021-11-16 | 中国地质大学（武汉） | 基于多线性奇异值张量分解核磁共振fid信号噪声抑制方法 |
| CN114662045A (zh) * | 2022-03-24 | 2022-06-24 | 电子科技大学 | 基于框架集的p阶张量深度学习的多维地震数据去噪方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110568486B (zh) | 基于同步稀疏低秩张量补全模型的地震信号补全方法 | |
| CN107728211B (zh) | 基于张量核范数正则化的地震信号恢复算法 | |
| Jafarpour et al. | History matching with an ensemble Kalman filter and discrete cosine parameterization | |
| CN111929733A (zh) | 基于slice采样的地震信号正则化处理方法 | |
| Sang et al. | DCNNs-based denoising with a novel data generation for multidimensional geological structures learning | |
| Zhang et al. | Adjoint-driven deep-learning seismic full-waveform inversion | |
| Zhang et al. | U-net generative adversarial network for subsurface facies modeling | |
| Li et al. | Consecutively missing seismic data interpolation based on coordinate attention unet | |
| US10436924B2 (en) | Denoising seismic data | |
| CN107561576A (zh) | 基于字典学习正则化稀疏表示的地震信号恢复方法 | |
| Zhu et al. | Sparse-promoting full-waveform inversion based on online orthonormal dictionary learning | |
| CN105277985A (zh) | 一种基于图像处理的ovt域地震数据规则化方法 | |
| Huang et al. | Self-supervised deep learning to reconstruct seismic data with consecutively missing traces | |
| Wang et al. | Seismic velocity inversion transformer | |
| Han et al. | Hyperspectral unmixing via nonconvex sparse and low-rank constraint | |
| Yang et al. | FWIGAN: Full‐Waveform Inversion via a Physics‐Informed Generative Adversarial Network | |
| Zhang et al. | 2-D seismic data reconstruction via truncated nuclear norm regularization | |
| CN114966860A (zh) | 一种基于卷积神经网络的地震数据去噪方法 | |
| Gao et al. | Incorporating structural constraint into the machine learning high-resolution seismic reconstruction | |
| Yu et al. | Data-driven geophysics: From dictionary learning to deep learning | |
| Gao et al. | Deep learning vertical resolution enhancement considering features of seismic data | |
| Liu et al. | Seismic data interpolation using generalised velocity‐dependent seislet transform | |
| Villarreal et al. | Seismic source reconstruction in an orthogonal geometry based on local and non-local information in the time slice domain | |
| Lin et al. | SeisGAN: Improving seismic image resolution and reducing random noise using a generative adversarial network | |
| CN106291676A (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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20201113 |