一种基于张量近似的海量数据多分辨率体绘制方法

摘要

本发明公开一种基于张量近似的海量数据多分辨率体绘制方法，首先对原始数据进行分块处理，得到若干数据块，然后对每一个数据块进行张量分解和多分辨率处理，最后对经张量分解和多分辨率处理后的每一个数据块进行重构处理，并创建二维纹理，完成地震数据的绘制。通过采用秩截断的方法有效的过滤原始数据中的噪声，并采用因子矩阵和核心张量秩截断的方式来替代对每一次秩试探所做的张量分解，准确的对每一个数据块确定秩的大小，节省数据块最佳秩选取的时间，并根据秩截断后的得到的数据，选择每个数据块的细节水平，快速有效地降低数据的整体分辨率大小，减少处理时间，从而实现多分辨率处理，并得到比传统基于信息熵的多分辨率更好的绘制效果。

说明书

技术领域

本发明属于图像处理领域，具体涉及一种体绘制技术。

背景技术

对体数据的可视化技术是一项非常常用的技术，能够被广泛应用于许多领域中，例如：医学领域、流体物理领域、气象领域、地质勘探领域等等。由于人类对于视觉信号的刺激最为敏感，因此可视化技术可以将含有大量信息的文件、图片、表格等转换为三维图像，便于研究人员直观地去观察和分析它们。

由于在现实生活中，常见的三维物体都是其表面，因此，在基于信息熵的模型中，常采用表面表示的方式去绘制一个三维体。然而很多时候，人们所关心的恰恰是物体的内部结构。例如在医学领域，人们更希望通过可视化技术观察器官内部的组织是否有病变；在地质勘探领域，人们更希望通过可视化技术观察地层以下的各种地质结构。因此此时，表面绘制就有了很大的局限性，它并不能帮助研究人员获得他们需要的信息。所以，人们需要一个能够观测物体内部结构信息的可视化技术——体绘制技术。

体绘制技术的本质，就是通过技术手段，用一个三维的标量数据生成一张二维的图片，并将其显示在电脑屏幕上。体绘制能够对物体的内部结构进行透视，使观察者能够一览物体的整体而不仅仅是表面，因此具有广泛的应用价值。

然而，随着数据探测技术的发展，人们获得的数据量成几何倍数增长，加之现在全球进入了互联网时代，互联网上的数据更是呈爆炸式增长。因此，用于体绘制的体数据量越来越大。由于计算机寻址空间的限制，体数据的大小很容易就超过了计算机的显存甚至内存大小。

于是，人们开始使用并行绘制来解决大数据量的问题。通过对体数据的不同部分进行分布式计算的方法来完成对整个体数据的绘制。但是，由于在体数据中，我们所关心的目标数据往往很小，采用并行体绘制显得得不偿失。

在这种情况下，产生了多分辨率体绘制技术。多分辨率体绘制技术是通过将数据体划分为不同的分块，每一个分块赋予不同的分辨率来实现压缩数据量和减少绘制点数。每个分块的分辨率又被称为分块的LOD(levelofdetail,细节水平)。因此，多分辨率技术的关键即为如何确定每一个分块的细节水平。

通过计算分块的信息熵(方差)来确定一个分块的细节水平是一种十分常用的方法。它通过计算分块的方差来得到分块的同质性。通常认为同质性高的分块，所包含的信息量小，因此不需要较高的分辨率。因而，方差越大的分块，同质性越低，说明其蕴含的信息量越大，因此其需要较高的分辨率。

对于某些数据，传统的依靠计算分块的数值特征并不能很好的实现多分辨率处理。例如，在地质勘探中常用的地震数据。由于地震数据具有信噪比低，变化剧烈，同质区域较少的特点，因此每个分块的方差都很高，多分辨率处理后的分块都普遍具有较高的分辨率。同时由于地震数据中存在大量的噪声和其他无用信息，因而高的方差并不能代表其具有较高的信息量。所以，单纯的多分辨率处理并不能有效的降低地震数据的数据量。另一方面，理论上的地震数据具有十分显著的微结构特征，不同的结构值之间的轮廓十分清晰，例如层位、断层。但是，由于地震数据中存在着大量的噪声和无用信息，导致实际应用中通过传感器探测到的数据在局部是混乱的，结构之间界限是模糊的，从而导致研究人员不能方便地从中分辨出他们所关心的结构。因此需要将地震数据中人们所关心的结构特征提取出来，过滤掉噪声和无用信息，然后将其可视化，提供直观的清晰的显示。此时，基于信息熵的多分辨率处理已经不能满足有效降低数据量和提取结构特征的要求了。

发明内容

本发明为解决的上述技术问题，提出一种基于张量近似的海量数据多分辨率体绘制方法，通过采用秩截断的方法有效的过滤原始数据中的噪声，并通过自适应秩截断的方法，准确的对每一个分块确定秩的大小，根据秩截断后的数据，选择分块的细节水平，从而实现多分辨率处理，并得到比传统基于信息熵的多分辨率更好的绘制效果。

本发明采用的技术方案是：一种基于张量近似的海量数据多分辨率体绘制方法，首先对原始数据进行分块处理，得到若干数据块，然后对每一个数据块进行张量分解和多分辨率处理，最后对经张量分解和多分辨率处理后的每一个数据块进行重构处理，并创建二维纹理，完成地震数据的绘制。

进一步地，具体包括以下步骤：

S1：对原始数据进行分块，得到若干数据块；

S2：对步骤S1得到的每一个数据块进行张量分解；

S3：对步骤S2得到的每一个数据块进行多分辨率处理；

S4：对步骤S3得到的每一个数据块进行重构，并创建二维纹理，根据重构得到的数据块进行地震数据绘制。

进一步地，所述步骤S2包括以下分步骤：

S21：对每个数据块根据各自的初始秩进行张量分解，得到因子矩阵和对应核心张量，并设置初始秩截断参数秩为1；

S22：对数据块的因子矩阵和对应的核心张量根据当前秩截断参数秩进行秩截断，得到秩截断后的因子矩阵和对应的核心张量；

S23：根据步骤S22得到的因子矩阵和对应的核心张量，进行该数据块的重构，并计算重构误差；

S24：判断步骤S23得到的重构误差是否满足收敛条件，若是，则进行步骤S25，否则，进行步骤S26；

S25：输出数据块在当前秩截断参数秩得到的因子矩阵和对应的核心张量；

S26：当前秩截断参数秩进行自加一操作，重复步骤S22至S25，得到每个数据块在各自对应当前秩截断参数秩下的因子矩阵和对应的核心张量。

更进一步地，步骤S21所述的对每个数据块的秩进行初始化，该初始秩根据每个数据块的分块尺寸确定。

更进一步地，所述步骤S22计算重构误差，具体为：

$e=\frac{\left|\right|A-\tilde{A}|{|}_F}{\left|\right|A|{|}_F}$

其中，e表示重构误差，A表示原始张量，表示重构后的近似张量，||||_F表示矩阵的范书。

更进一步地，步骤S23所述的判断步骤S22得到的重构误差是满足收敛条件具体为：判断当前秩截断参数秩是否小于或等于数据块的初始化秩R；或者，判断当前秩截断参数秩的重构误差e是否小于或等于当前秩截断参数秩的归一化重构误差T_e；或者，判断当前秩截断参数秩的重构误差e是否满足下式：

$\frac{e^{\prime }-e}{e^{\prime }}\le T_p$

其中，e′表示上一秩截断参数秩的重构误差，T_p表示当前秩截断参数秩的重构误差的提升值。

进一步地，所述步骤S3具体包括以下分步骤：

S31：根据数据块的分块个数为2的k次幂，得到分辨率为k个级别；

S32：根据初始的秩截断参数秩为1以及步骤S21中的每个数据块的秩的初始化值R，得到秩的取值范围为[1，R]，共有R个取值；

S33：根据(R-1+1)/(k+1)，得到商为a，余数为b，则秩的分组为b组包括a+1个秩，剩下的k+1-b组包括a个秩，随机组合得到分组排列；

S34：将秩的取值按照从大到小进行排序，并根据步骤S33得到的分组排列进行分组。

进一步地，所述步骤S4根据下式对数据块进行重构：

$\tilde{A}=\underset{r_1}{\Sigma }\underset{r_2}{\Sigma }\underset{r_3}{\Sigma }{b}_{r_1{r}_2{r}_3}·{\overrightarrow{u}}_{r_1}^{\left(1\right)}·{\overrightarrow{u}}_{r_2}^{\left(2\right)}·{\overrightarrow{u}}_{r_3}^{\left(3\right)}$

其中，表示核心张量位于(r₁,r₂,r₃)位置的值，表示矩阵U⁽ⁿ⁾的第r_n列的列向量。

本发明的有益效果：本发明的一种基于张量近似的海量数据多分辨率体绘制方法，通过采用秩截断的方法有效的过滤原始数据中的噪声并通过自适应秩截断的方法，采用因子矩阵和核心张量秩截断的方式来替代对每一次秩试探所做的张量分解，准确的对每一个数据块确定秩的大小，节省数据块最佳秩选取的时间，并根据秩截断后的得到的数据，选择每个数据块的细节水平，快速有效地降低数据的整体分辨率大小，减少处理时间，从而实现多分辨率处理，并得到比传统基于信息熵的多分辨率更好的绘制效果。

附图说明

图1为本发明提供的方法流程图。

图2为本发明提供的地震数据多分辨率绘制效果图；

其中，a图为基于信息商的多分辨率绘制效果图；b图为基于张量近似的多分辨率绘制效果图。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

如图1所示，本发明的一种基于张量近似的海量数据多分辨率体绘制方法，包括以下步骤：

S1：对原始数据进行分块，得到若干数据块；

S2：对步骤S1得到的每一个数据块进行张量分解；

S3：对步骤S2得到的每一个数据块进行多分辨率处理；

S4：对步骤S3得到的每一个数据块进行重构，并创建二维纹理，根据重构得到的数据块进行地震数据绘制。

所述步骤S1具体为：分块尺寸设置的大小，直接关系到每个数据块内的信息量大小，如果数据块尺寸设置的大小过小，那么会导致每个数据块内的信息量太少，信息的连续性太低，整体的近似效果离散化严重；如果数据块尺寸设置得过大，那么会导致整体数据块个数的减小，有可能会影响整体数据的压缩效果。通过大量的仿真试验，得出分块的边长大小设置为32或64是最为理想的。

所述步骤S2对步骤S1得到的每一个数据块进行张量分解；采用自适应秩的方法，自适应地秩截断能够在保证一定的准确度的条件下，根据每一个分块的不同特点，自适应地选择每一个分块的秩大小。从而实现对每一个分块进行不同程度的数据压缩。具体为：

S21：对每个数据块根据各自的初始秩进行张量分解，得到因子矩阵和对应核心张量，并设置初始秩截断参数秩为1；例如本申请的分块大小设置为32，则本申请初始时对数据块做秩16分解，因为数据块的分块大小只能是2的幂，例如本申请数据块分块大小为2的5次幂。本身请初始秩选择为数据块分块大小的一半，也就是2的4次幂；即将每个数据块的秩进行初始化为16，并进行张量分解，得到每个数据块的因子矩阵和核心张量，具体计算过程如下：

在确定完张量分解的秩以后，就可以对每个分块进行张量分解。对三维数据块的张量分解就是上述对n阶张量分解在n＝3时的特殊情况。一个三维数据块被分解成一个核心张量和三个因子矩阵和的TTM乘积：

$A_{I_1\times I_2\times I_3}\approx B_{R_1\times R_2\times R_3}\times {}_{1}U{{}^{\left(1\right)}}_{I_1\times R_1}\times {}_{2}U{{}^{\left(2\right)}}_{I_2\times R_2}\times {}_{3}U{{}^{\left(3\right)}}_{I_3\times R_3}.---\left(1\right)$

确定初始的秩截断参数秩大小，理论上，秩截断参数秩初始取值为1，在本申请中，对地震数据的张量近似实验，得到当R＝4时，张量近似的绘制效果开始大幅下降。因此，本申请设定初始的秩截断参数秩大小为R＝4。并以R＝4对因子矩阵与和核心张量进行秩截断，得到当前秩截断参数秩下的因子矩阵和对应的核心张量。

S22：对数据块的因子矩阵和对应的核心张量根据当前秩截断参数秩进行秩截断，得到秩截断后的因子矩阵和对应的核心张量。

S23：根据步骤S22得到的因子矩阵和对应的核心张量，进行该数据块的重构，并计算重构误差；根据当前因子矩阵和对应的核心张量，进行该数据块的重构，并计算重构误差；重构误差的计算可以根据实际的需要选择合适的标准。本申请采用的是矩阵的Frobenius范数作为重构误差的标准：

$e=\frac{\left|\right|A-\tilde{A}|{|}_F}{\left|\right|A|{|}_F}---\left(2\right)$

其中，e代表重构后的归一化误差，A代表原始张量，代表重构后的近似张量，对大小为M×N的矩阵B做||||_F运算定义为：

$\left|\right|B|{|}_F=\sqrt{\underset{i=1}{\overset{M}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{b}_{ij}^2}---\left(3\right)$

S24：判断步骤S23得到的重构误差是否满足收敛条件，若是，则进行步骤S25，否则，进行步骤S26；对于收敛条件来说，首先，做秩截断的参数秩大小显然不应该超过初始张量分解时设定的初始秩大小。即，最终确定的秩截断参数秩大小应该小于等于16。其次，按当前秩截断参数秩大小重构后的误差大小，应该小于所设定的当前秩的归一化重构误差。最后，如果误差大小仍大于当前秩截断参数秩的归一化重构误差，但是每增加一次秩截断参数秩所带来的准确度的提升十分的小，那么继续增加秩截断参数秩对于减少误差来说意义很小，只会地增加数据量。因此，在收敛条件中，还应该判断当前秩截断参数秩的重构误差与上一秩截断参数秩的重构误差相比，准确度的提升是否显著。如果提升不显著，则停止继续增加秩截断参数秩。因此，本申请将收敛条件归纳如下：

R＝16(4)

$e=\frac{\left|\right|A-\tilde{A}|{|}_F}{\left|\right|A|{|}_F}\le T_e---\left(5\right)$

$\frac{e^{\prime }-e}{e^{\prime }}\le T_p---\left(6\right)$

只要满足式(4)，(5)和(6)中任意其一，即满足收敛条件。其中，e′为上一秩截断参数秩的重构误差，T_e和T_p分别为当前秩截断参数秩的归一化重构误差和当前秩截断参数秩的重构误差的提升。根据实际需要，用户可以设置合适的T_e和T_p值。

S25：输出数据块在当前参数秩得到的因子矩阵和对应的核心张量；输出根据步骤S23得到的当前秩截断参数秩对因子矩阵和核心张量做秩截断，得到的数据块对应的当前秩截断参数秩下的因子矩阵以及对应的核心张量。

S25：秩截断的秩进行自加一操作，并根据得到的秩截断参数秩对因子矩阵与和核心张量进行秩截断，得到当前秩截断参数秩下的因子矩阵和对应的核心张量，并转至步骤S22。

S26：当前秩截断参数秩进行自加一操作，重复步骤S22至S25，得到每个数据块在各自对应当前秩截断参数下的因子矩阵和对应的核心张量。

所述步骤S3对步骤S2得到的每一个数据块进行多分辨率处理，主要采用基于分块秩大小的细节水平选择算法，能够有效地降低整体的分辨率大小，并且在细节水平选择过程中能够直接确定其分块的分辨率，减少处理时间。具体包括以下分步骤：

S31：根据数据块的分块尺寸为2的k次幂，得到分辨率为k+1个级别；例如本申请中数据块的分块尺寸为2的5次幂，得到分辨率为0—5，共6个级别。

S32：根据初始的秩截断参数秩为1以及步骤S21中的每个数据块的秩的初始化值R，得到秩的取值范围为[1，R]，共有R个取值；例如，本申请中的秩截断的初始值为4，每个数据块第一次张量分解时的秩为16，则秩的取值范围为[4，16].

S33：根据(R-1+1)/(k+1)，得到商为a，余数为b，则秩的分组为b组包括a+1个秩，剩下的k+1-b组包括a个秩，随机组合得到分组排列；根据(16-4+1)/(5+1)＝13/6，商2，余1，则秩的分组为5组包括两个秩，剩下的一组包括3个秩。随机取组合为前五组各包括两个秩，最后一组包括3个秩。

S34：将秩的取值按照从大到小进行排序，并根据步骤S33得到的分组排列进行分组。根据本申请的取值得到的秩分组与分辨率级别对应关系如表1所示。

表1基于秩大小的细节水平选择

秩大小 4或5 6或7 8或9 10或11 12或13 14,15或16 细节水平 第0级 第1级 第2级 第3级 第4级 第5级

所述步骤S4：对步骤S3得到的每一个数据块进行重构，并创建二维纹理，根据重构得到的数据块进行地震数据绘制。由于式(1)的并行计算度不高，在实时重构的绘制中，不利于GPU的并行加速。因此，我们可以将其改写为并行度较高的等价形式，如式(7)所示：

$\tilde{A}=\underset{r_1}{\Sigma }\underset{r_2}{\Sigma }\underset{r_3}{\Sigma }{b}_{r_1{r}_2{r}_3}·{\overrightarrow{u}}_{r_1}^{\left(1\right)}·{\overrightarrow{u}}_{r_2}^{\left(2\right)}·{\overrightarrow{u}}_{r_3}^{\left(3\right)}---\left(7\right)$

其中，表示核心张量B的位于(r₁,r₂,r₃)位置的值，代表矩阵U⁽ⁿ⁾的第r_n列列向量。具体的体绘制技术为常规技术，故，在此不做过多解释。

通过验证，本发明的方法的效果与基于信息熵的未进行张量近似的多分辨率体绘制效果进行了对比。如图2所示，上面两图为基于信息熵多分辨率和基于张量近似多分辨率体绘制的整体绘制效果，下面两图分别为矩形框范围的细节放大。显然，从右图中可以清晰地看到地震数据中的每一个小层位的结构，而在左图中则十分模糊。从中可以看出，相较于基于信息熵多分辨率体绘制，基于张量近似的多分辨率体绘制的绘制图像对于地震数据的结构特征表现得更明显。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。