基于并行矩阵分解的低秩五维地震数据重构方法及系统

1.本发明属于信号处理技术领域，具体涉及一种基于并行矩阵分解的低秩五维地震数据重构方法及系统。


背景技术：

2.由于地震数据采集过程中某些障碍物，如不规则的地表条件和设备的限制，地震数据往往沿空间坐标不规则地采样。地震记录中缺失的迹线会对后续处理、成像和涉及多通道反卷积、速度分析、全波形反演、同步源分离和故障检测的处理、成像和相互伪装产生负面影响。因此，地震数据重建是一个持续存在的重要问题，引起了学术界和工业界的广泛关注。
3.对于地震数据的重建问题，通常假设完整数据是低秩的，并且在缺少迹线和存在噪声中会导致秩的增大。
4.近些年，基于张量网络的张量分解理论已经获得了越来越多的关注，在处理高阶张量方面显示出非凡的前景，特别是在三阶以上的张量。其中，tensor traintt分解最具代表性，具有高数据压缩能力和计算效率，在许多任务中得到了深入研究。通常使用svd，交替最小二乘法，以及梯度下降算法计算tt分解的潜在核心张量，以利用其低秩结构。然而，尽管存在许多优化算法，但由于计算成本高，显式tt分解对于大规模应用仍然难以解决。本发明通过最小化tt秩，以隐式方式利用低秩结构，对不同的张量展开执行并行低秩矩阵分解，较传统方法大大降低了计算成本。
5.现有技术：
6.数据驱动技术，包括字典学习方法，涉及k奇异值分解，和数据驱动的紧帧方法，可以以数据驱动的方式直接从输入数据中学习特征，以数据驱动的方式直接从输入数据中学习特征，具有很强的适应性。
7.现有技术的缺点：
8.1高维情况下，数据驱动技术的计算复杂性高，参数选择对训练结果过于灵敏，使其不能稳定完成计算任务。
9.2现有计算机对于数据驱动的方法的处理能力仅限于三维数据，不足以利用五维地震数据中固有的全面高维关系。


技术实现要素：

10.本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于并行矩阵分解的低秩五维地震数据重构方法及系统，用于解决缺失迹线的地震数据重建的技术问题。
11.本发明采用以下技术方案：
12.基于并行矩阵分解的低秩五维地震数据重构方法，包括以下步骤：
13.s1、输入观测数据张量采样张量优化项与正则项的权重参数μ，以及将观
测数据张量沿各方向按tt方法分解得到的矩阵的秩r1，
…
，rn，初始化参数得到各方向初始分解矩阵和初始低秩重建数据张量初始化迭代变量k＝0，作为初始迭代步骤；
14.s2、对于第k次迭代过程，由真实数据的低秩约束，展开步骤s1得到的低秩重建数据张量利用sketch随机化采样与块坐标下降法，计算每个方向的矩阵使用及对应分解矩阵计算当前迭代步骤的分解矩阵
15.s3、将步骤s2得到的分解矩阵用于子优化问题低秩重建数据张量的求解，令k＝k 1，进入下一个迭代步骤，重复步骤s2和步骤s3，将迭代得到的作为真实数据的低秩重建结果，实现低秩五维地震数据重构。
16.具体的，步骤s2中，由真实数据的低秩约束，对真实数据的矩阵展开进行低秩分解，并引入观测数据张量对真实数据的约束项，最小化求解真实数据的逼近，具体如下：
[0017][0018]
其中，z
[n]
为对真实数据在第n个方向的tt分解结果，xn与yn为z
[n]
的低秩分解结果，f代表frobenius范数运算，n为张量的维数。
[0019]
进一步的，通过观测数据张量求解真实数据的低秩重建表示如下：
[0020][0021]
更进一步的，低秩重建的约束表示如下：
[0022][0023]
进一步的，步骤s2中，使用迭代法交替优化各tt分解矩阵的低秩分解矩阵xn，yn，得到xn，yn的优化结果如下：
[0024][0025][0026]
其中，为真实数据在第n个方向的tt分解矩阵。
[0027]
更进一步的，采用sketch算法加速低秩分解矩阵xn的计算，改进的优化子问题如下：
[0028][0029]
其中，和为经sketch算法采样yn和得到的结果。
[0030]
具体的，步骤s2中，分解矩阵分别为：
[0031]
[0032][0033]
其中，为对真实数据在第n个方向的tt分解结果，为z
[n]
的低秩分解结果的左矩阵。
[0034]
具体的，步骤s3中，使用迭代计算真实数据的低秩重建数据张量如下：
[0035][0036]
其中，为全1张量，α为重构张量与观测数据张量的比例系数，为各方向的分解矩阵和的重构张量。
[0037]
进一步的，步骤s3中，低秩重建数据张量的子优化问题：
[0038][0039]
其中，n为张量的维数，和为分解矩阵，z
[n]
为对真实数据在第n个方向的tt分解结果，f代表frobenius范数运算。
[0040]
第二方面，本发明实施例提供了一种基于并行矩阵分解的低秩五维地震数据重构系统，包括：
[0041]
数据模块，输入观测数据张量采样张量优化项与正则项的权重参数μ，以及将观测数据张量沿各方向按tt方法分解得到的矩阵的秩r1，
…
，rn，初始化参数得到各方向初始分解矩阵和初始低秩重建数据张量初始化迭代变量k＝0，作为初始迭代步骤；
[0042]
计算模块，对于第k次迭代过程，由真实数据的低秩约束，展开数据模块得到的低秩重建数据张量利用sketch随机化采样与块坐标下降法，计算每个方向的矩阵使用及对应分解矩阵计算当前迭代步骤的分解矩阵
[0043][0044]
重构模块，将计算模块得到的分解矩阵用于子优化问题低秩重建数据张量的求解，令k＝k 1，进入下一个迭代步骤，重复执行，将迭代得到的作为真实数据的低秩重建结果，实现低秩五维地震数据重构。
[0045]
与现有技术相比，本发明至少具有以下有益效果：
[0046]
本发明基于并行矩阵分解的低秩五维地震数据重构方法，首先利用观测地震数据张量初始化重建张量，利用重建地震数据具有较低的秩的特点，将地震张量利用tt方法分解成一系列矩阵，相较于传统分解方法，充分利用了观测张量的高维信息；将这一系列矩阵利用sketch方法低秩分解，利用低秩分解结果重建地震数据张量，相较于直接进行矩阵低秩分解，使用sketch方法可以显著加快运算速度；迭代求解直到符合预定的停止准则，在缺失观测地震数据上得到良好的重建结果。
[0047]
进一步的，对真实数据的矩阵展开进行低秩分解。使用tensor traintt秩最小化pmf-tt方法处理数据。为了恢复真实数据中的缺失值，需要使得的真实数据秩尽量小，同
时要求恢复出的数据的采样结果与观测数据张量尽量接近。
[0048]
进一步的，通过观测数据张量求解真实数据的低秩重建，通过引入抽样算子将对真实数据的抽样过程转化为了张量的哈达玛积运算，方便了优化问题的提出与计算过程的形式化表示。
[0049]
进一步的，为了恢复真实数据中的缺失值，需要使得的秩尽量小。低秩重建的优化问题表示，设计的优化问题将约束条件转化为优化的正则项，方便进一步的求解。
[0050]
进一步的，分别优化xn，yn，得到真实数据的优化结果，通过设置分解矩阵xn，yn的优化问题，使得优化出的xn，yn与真实数据分解的对应分量尽量接近，重建的张量结果也与真实数据尽量接近。
[0051]
进一步的，采用sketch算法可以有效减小程序的运行时间。
[0052]
进一步的，通过计算各迭代步骤的分解矩阵可以计算出该迭代步骤的低秩重建结果张量，更新低秩真实数据张量重建结果
[0053]
进一步的，利用当前步骤迭代得到的和求解真实数据构建的优化子问题，通过求解该优化子问题，得到对于真实数据的最优估计。
[0054]
进一步的，通过设计对的迭代步骤，随着迭代步骤的增加，对的重构结果也更加接近真实数据，直到符合预定的停止准则。
[0055]
可以理解的是，上述第二方面的有益效果可以参见上述第一方面中的相关描述，在此不再赘述。
[0056]
综上所述，本发明采用均衡的张量矩阵化的方法，具有优秀的重建性能，同时显著降低计算成本。
[0057]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0058]
图1为pmf-tt分解的示意图；
[0059]
图2为sketch抽样的示意图；
[0060]
图3为地震数据图，其中，(a)为合成地震数据，(b)为90％损失地震数据；
[0061]
图4为90％损失地震数据重建结果，其中，(a)为90％损失，无噪地震数据使用pmf方法的重建结果，(b)为90％损失，无噪地震数据使用pmf-sketch方法的重建结果，(c)为90％损失，无噪地震数据使用pmf-tt方法的重建结果，(d)为90％损失，无噪地震数据使用pmf-tt-sketch方法的重建结果，(e)，(f)，(g)，(h)分别为(a)，(b)，(c)，(d)与真实数据的重建误差；
[0062]
图5为带噪地震数据，其中，(a)为合成地震数据，(b)为90％损失带噪地震数据；
[0063]
图6为90％损失带噪地震数据重建结果，其中，(a)为90％损失，含1db噪声地震数据使用pmf方法的重建结果，(b)为90％损失，含1db噪声地震数据使用pmf-tt方法的重建结果，(c)为90％损失，含1db噪声地震数据使用pmf-tt方法的重建结果，(d)为90％损失，含1db噪声地震数据使用pmf-tt-sketch方法的重建结果，(e)，(f)，(g)，(h)分别为(a)，(b)，(c)，(d)与真实数据的重建误差；
[0064]
图7为本发明的流程框图。
具体实施方式
[0065]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0066]
在本发明的描述中，需要理解的是，术语“包括”和“包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
[0067]
还应当理解，在本发明说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本发明。如在本发明说明书和所附权利要求书中所使用的那样，除非上下文清楚地指明其它情况，否则单数形式的“一”、“一个”及“该”意在包括复数形式。
[0068]
还应当进一步理解，在本发明说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合，例如，a和/或b，可以表示：单独存在a，同时存在a和b，单独存在b这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。
[0069]
应当理解，尽管在本发明实施例中可能采用术语第一、第二、第三等来描述预设范围等，但这些预设范围不应限于这些术语。这些术语仅用来将预设范围彼此区分开。例如，在不脱离本发明实施例范围的情况下，第一预设范围也可以被称为第二预设范围，类似地，第二预设范围也可以被称为第一预设范围。
[0070]
取决于语境，如在此所使用的词语“如果”可以被解释成为“在
……
时”或“当
……
时”或“响应于确定”或“响应于检测”。类似地，取决于语境，短语“如果确定”或“如果检测(陈述的条件或事件)”可以被解释成为“当确定时”或“响应于确定”或“当检测(陈述的条件或事件)时”或“响应于检测(陈述的条件或事件)”。
[0071]
在附图中示出了根据本发明公开实施例的各种结构示意图。这些图并非是按比例绘制的，其中为了清楚表达的目的，放大了某些细节，并且可能省略了某些细节。图中所示出的各种区域、层的形状及它们之间的相对大小、位置关系仅是示例性的，实际中可能由于制造公差或技术限制而有所偏差，并且本领域技术人员根据实际所需可以另外设计具有不同形状、大小、相对位置的区域/层。
[0072]
本发明提供了一种基于并行矩阵分解的低秩五维地震数据重构方法，构建了通过观测信号计算重构的完整信号的完整方法；通过迭代求解，将上一步的结果进行张量分解，并进行并行矩阵分解获得低秩矩阵，利用低秩矩阵重建张量，不断迭代得到最终的低秩张量。除此之外，使用sketch随机采样方法进一步加快运算速度，采用均衡的张量矩阵化的方法，具有优秀的重建性能，同时显著降低计算成本。
[0073]
请参阅图7，本发明一种基于并行矩阵分解的低秩五维地震数据重构方法，包括以下步骤：
[0074]
s1、给定已知的数据张量采样张量优化项与正则项的权重参数μ，以及沿个方向分解按tt方法分解出的矩阵的秩r1，
…
，rn；初始化参数y0，作为算法的初始值；
[0075]
tensor traintt秩最小化pmf-tt方法处理的数据是高维张量数据。对于5维地震数据，如不考虑频率维度，对剩下的一系列4维张量进行处理。假设完整的真实数据为观测数据张量为张量的维数为n，认为观测数据张量是真实数据的抽样结果，表示为：
[0076][0077]
其中，表示哈达玛积按元素乘积，是抽样算子，观测到的位置取1，其他位置取0。
[0078]
通过求解的低秩重建的求解问题表示为：
[0079][0080]
为了恢复中的缺失值，需要使得的秩尽量小。上式的约束写为：
[0081][0082]
其中，是的低秩约束的泛化表示。
[0083]
上式的一种有效解法是对的矩阵展开进行低秩分解，将上式转化为：
[0084][0085]
其中，是在n方向上的tt分解，在n方向上的tt分解，rn是预先设定的xn与yn的秩。
[0086]
s2、展开当前步的低秩重建数据张量对每个方向的矩阵利用随机化采样与块坐标下降法计算的分解矩阵
[0087]
通过分别优化xn，yn，求解子问题来得到的优化结果：
[0088][0089][0090]
对于每一个xn和ynn＝1，
…
，n-1，通过下式迭代求解：
[0091][0092][0093]
采用sketch随机采样算法极大的减小运算时间，对于如果sn≤tn，此时求解xn是一个超定问题，由对z
[n]
的列进行抽样求得，反之亦然。抽样尺寸选择ssn＝10rnlog
10
(rn)，对yn和z
[n]
的对应列分别随机抽取ssn个向量，构成和构建对于x的优化子问题：
[0094][0095]
得到xn及yn的递推公式：
[0096][0097][0098][0099]
s3、将步骤s2得到的分解矩阵及经张量重建求得下一步的重复步骤s2和步骤s3，直到符合预定的停止准则。
[0100]
利用当前步骤迭代得到的x和y求解首先写出子优化问题：
[0101][0102]
求得真实数据的迭代为：
[0103][0104]
其中，
[0105]
在每次迭代过程中判断迭代的终止条件，当终止条件满足时停止迭代，此时获得的即视作观测数据张量的低秩重构结果。
[0106]
本发明再一个实施例中，提供一种基于并行矩阵分解的低秩五维地震数据重构系统，该系统能够用于实现上述基于并行矩阵分解的低秩五维地震数据重构方法，具体的，该基于并行矩阵分解的低秩五维地震数据重构系统包括数据模块、计算模块以及重构模块。
[0107]
其中，数据模块，输入观测数据张量采样张量优化项与正则项的权重参数μ，以及将观测数据张量沿各方向按tt方法分解得到的矩阵的秩r1，
…
，rn，初始化参数得到各方向初始分解矩阵和初始低秩重建数据张量初始化迭代变量k＝0，作为初始迭代步骤；
[0108]
计算模块，对于第k次迭代过程，由真实数据的低秩约束，展开数据模块得到的低秩重建数据张量利用sketch随机化采样与块坐标下降法，计算每个方向的矩阵使用及对应分解矩阵计算当前迭代步骤的分解矩阵
[0109][0110]
重构模块，将计算模块得到的分解矩阵用于子优化问题低秩重建数据张量的求解，令k＝k 1，进入下一个迭代步骤，重复执行，将迭代得到的作为真实数据的低秩重建结果，实现低秩五维地震数据重构。
[0111]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例
中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0112]
请参阅图1，为张量pmf-tt分解的过程；对于张量将其分解为一系列矩阵x
[n]
，类似于利用张量tucker秩来进行张量重建，也可将其推广到张量的tt分解中，利用tt分解矩阵的秩进行张量重建。张量的最小tt分解的形式：
[0113][0114]
其中，是一系列3维张量，核的大小分别为1
×
i1×
r1,1
×
i2×
r2,
…
,1
×in
×rn
，《
·
》是张量核的多线性积。
[0115]
将{r1，r2，
…
，r
n-1
}记作rank
tt
(x)，称为张量x的tt秩；将张量x分解成如下形式：
[0116][0117]
其中，是的第in个切片。
[0118]
通过最小化tt分解得到tt秩，但是由于直接计算tt分解需要大量的计算，因此有必要使用一种有效率的方法计算tt秩。
[0119]
x的tt分解和tt核具有如下关系：
[0120][0121]
其中，是部分乘积。
[0122]
分解结果和tt秩具有如下关系：
[0123][0124]
x
[n]
是rn的下界，当rank(x
[n]
)＝rn时即为最小tt分解，使用rank(x
[n]
)表示tt秩，解决tt秩最小化问题。
[0125]
请参阅图2，为sketch抽样示意图。在迭代计算低秩张量的过程中，需要计算各模式矩阵的低秩分解。研究表明，随机化方法对于加速这一过程具有显著的效果。我们提出的sketch抽样方法通过对给定矩阵zn的行或列进行抽样来加速计算。对于当sn≤tn时是一个超定问题，因此对zn的列进行抽样，反之亦然。通过调整采样分布可能达到更好的效果，但为了加速计算，本发明采用均匀分布。
[0126]
请参阅图3，为合成的五维地震数据的大小为256
×
20
×
20
×
10
×
10，其中共中点维度有20条迹线，偏置维度有10条迹线，并且每条迹线在时间上各有256个取值。对数据进行随机抽样，分别移除其中的10％，30％，50％，70％，90％的迹线，构成一系列不完整数据。对1-70hz范围的数据进行重建，并重组成五维张量，变换至时间-空间域中。
[0127]
为了量化各方法的重建优劣，提出重建质量指标q：
[0128][0129]
其中和分别代表真实和重建结果在时间-空间域的结果。
[0130]
对于无噪声的情况，设置停止准则为最大迭代次数为k
iter
为300或相对误差小于10-4
，相对误差可由计算；重新插入参数α调节观测结果对重建结果的影响程度。由于无噪声，设置其为1；最后设置优化矩阵的秩，简单起见设置r＝r1＝
…
＝r
max
＝r
n-1
，并扫描秩从1秩12递增的重建结果。实验表明大部分情况下r＝4得到最优结果。由实验结果还能得到经过sketch抽样后，可以在保持高q值的情况下扩展r的取值范围。
[0131]
从实验结果发现pmf-tt方法在损失70％和90％的情况下结果差于普通的pmf方法，因为4阶张量不足以表现tt方法的优势；实验结果同样表明sketch方法可以获得更高的q值，因此，sketch方法也具有提升无噪信号重建精度的优势。
[0132]
从运行时间来看，pmf-tt方法的运行时间小于pmf方法，但是q值有所下降。也可以看出sketch方法可以提高运行效率。
[0133]
为了探索sketch方法具有高q值和更短计算时间的原因，在20hz的频率数据块下对两方法重复50次。可以发现sketch方法通过降低迭代轮数一定程度上降低了运行时间，并且通过加速每轮迭代的运行时间使得总运行时间下降。另外可以发现，pmf-tt方法高运行时间的原因主要在于大量的迭代次数。
[0134]
请参阅图4，为了更清晰的解释结果，图4显示了90％损失信号100次的重建结果。可以看出在使用sketch算法前pmf-tt算法的误差图上可以看到轻微的瑕疵现象，在使用sketch算法后瑕疵被明显弱化，证实了之前的结论。另外也证实了pmf-tt在地震数据重建中的有效性。
[0135]
请参阅图5，对于加入1db噪声的数据，其他参数与无噪数据基本相同。相对误差改用计算。pmf-tt比pmf得到更差的重建结果，但由于避免了具有昂贵时间消耗的排列过程，从而大幅减小了运行时间。
[0136]
请参阅图6，显示了sketch算法的优越性，不仅增强了重建效果，还减小了信号的泄露。
[0137]
综上所述，本发明一种基于并行矩阵分解的低秩五维地震数据重构方法及系统，能够快速准确的重建地震数据，并减小参数选择对结果的影响，提升信号重建精度，减小程序运行时间的优点，具有良好的重建效果，减小了计算资源消耗，具有广阔的工业应用前景。
[0138]
本领域内的技术人员应明白，本技术的实施例可提供为方法、系统、或计算机程序产品。因此，本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质包括但不限于磁盘存储器、cd-rom、光学存储器等上实施的计算机程序产品的形式。
[0139]
本技术是参照根据本技术实施例的方法、设备系统、和计算机程序产品的流程图
和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0140]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0141]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0142]
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。