CN103487833A - 一种三维地震体矢量量化压缩体绘制方法 - Google Patents

Abstract

本发明公开了一种三维地震体矢量量化压缩体绘制方法，该方法包括数据预处理、矢量量化、解压绘制三部分。通过引入目标特征，减少目标区域的压缩损失，提高非目标区域的压缩率；引入数据分级的概念，利于数据压缩；引入层位特征，结合该特征确定目标区域，利用该信息进行体数据的压缩；引入误差系数的概念，通过该误差系数可以有效的降低目标区域的失真度，提高整体的压缩率。通过理论和仿真分析表明，本发明在提高整体数据压缩率的情况下，提高了目标区域的峰值信噪比。

Description

一种三维地震体矢量量化压缩体绘制方法

技术领域

本发明涉及海量三维标量场体绘制可视化方法，具体是一种三维地震体矢量量化压缩体绘制方法。

背景技术

体绘制技术是根据三维数据场产生二维图像的技术。相对于面绘制技术而言，体绘制技术更注重的是展示数据的体细节而不是面细节，并在医学、地质勘探、气象分析、分子模型构造等多个领域具有广泛的应用前景。

随着体绘制技术的不断完善，人们提出了多种体绘制的算法，包括：光线投射法直接体绘制，错切变换法体绘制，基于纹理预积分的体绘制技术等。而随着对体绘制交互需求的日益提高，体绘制的加速技术成为该领域研究的热点。图形硬件的可编程性的发展，给可视化领域带来了显著的变化。体绘制已经有原来的软件体绘制发展成为图形硬件体绘制，并且基于图形硬件的体绘制技术在精度，速度方面的优势使GPU成为标准的体绘制平台。Cabral等提出的基于3D纹理映射硬件的直接体绘制，因具有较高的绘制效率也成为直接体绘制的重要方法之一。LEE B等和Chen Shihao等分别提出的算法是在前人的基础之上进行了继承和发展。Kiss G等提出通过在GPU上实现数据预处理和光线投影法来实现实时的高质量的体绘制技术。

然而随着测量技术的发展，三维数据场其数据规模越来越大。导致传统的基于图形硬件体绘制加速算法难以支撑应用的需求。主要体现在两个方面：其一是显存和主存本身容量限制；其二是数据在多级存储结构：磁盘-CPU-GPU之间传输，传输带宽的限制等。解决问题的一个重要思路是如何降低体数据的规模。基于此，Nathaniel Fout等提出并实现了压缩体绘制的方法。其基本思想是将体数据压缩方法和体绘制方法结合，通过数据压缩以降低数据的规模，从而解决显存容量和主存到显存带宽的限制。Ning等首次提出将矢量量化压缩算法应用于体数据压缩中。吴仲乐等在此基础之上使用基于区域分裂的LBG算法对大规模体数据进行矢量量化压缩。上述算法仅仅通过牺牲压缩比或增加解压难度来获得较好的图像质量，在海量数据体绘制中，难以获得较好的图像质量。Schneider等提出并实现了多层次矢量量化压缩（hierarchical vector quantization），赵利平等在上述算法的基础之上提出分类分层矢量量化算法（flag based classical hierarchical vector quantization）根据数据分块后的梯度信息对分块数据进行分类。该方法避免了空体元带来的影响，从而提高了压缩效率。

上述体绘制压缩方法在不同程度上有效解决了海量数据的体绘制问题。但针对如图1所示同质区域小，数值分布比较散乱的情况下（如地震勘探数据等），不同分块的相似性较低。统一的数据压缩方法在保证压缩率的情况下，图像质量较差。此外，通常情况下，人们只关注海量数体中的特定目标，其余大部分是冗余数据。针对这种情况Chao等提出的信息熵驱动的压缩体绘制方法。Guoqing Wu等在此基础之上进一步优化，通过在进行数据压缩和数据传输时跳过非重要性数据，减小传递到内存的数据量。但在解决数值散乱分布的体数据情况下，单纯的信息熵无法标定体数据中的目标区域。

综上，矢量量化压缩体绘制是解决海量三维标量场体绘制可视化的一种主要方法。其基本思想是将海量体数据进行矢量量化压缩至显存可以容纳的空间，以提高绘制效率，达到交互体绘制的目的。现有的方法在解决同质区域小，数值分布比较散乱的地震数据体绘制时，压缩率较低。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种地震数据的压缩提绘制方法，减少目标区域的压缩损失，提高非目标区域的压缩率，同时在提高整体数据压缩率的情况下，提高目标区域的峰值信噪比。

为解决上述技术问题，本发明所采用的技术方案是：一种三维地震体矢量量化压缩体绘制方法，该方法为：

1)将原始地震数据分割成大小为N*N*N大小的数据块；

2)计算每一个数据块的梯度值，根据给定的阈值（该阈值的绝对值小于100），将梯度值低于所述阈值的数据块标记为零块，并计算所述零块的均值；将梯度值不低于所述阈值的数据块标记为非零块；

3)对每一个数据块中相邻的（N/2）*（N/2）*（N/2）的数据求平均值，获得大小为（N/2）*（N/2）*（N/2）的小数据块；计算每一个小数据块的均值；

4)对所述数据块与所述小数据块进行差分处理，得到最高层数据，即大小为N*N*N的块差分数据；对所述小数据块与所述小数据块的均值进行差分处理，得到次高层数据，即大小为（N/2）*（N/2）*（N/2）的差分矢量；将所述最高层数据经矢量量化后，得到码书一；将所述次高层数据经矢量量化后，得到码书二；

5)将上述步骤4）处理后得到的N*N*N的块差分数据或（N/2）*(N/2)*(N/2)的差分矢量作为输入矢量I，将所述输入矢量I的集合作为一个量化单元V₁，计算量化单元V₁的失真度D₁，进入步骤6）；V₁的形心矢量Y₁就是块差分数据对应的码书一的一个码字或差分矢量对应的码书二的一个码字；

6)根据上述量化单元V₁构造双向链表，所述双向链表中每一个分组对应步骤5）的一个码字，遍历所述双向链表，将所述双向链表分裂为若干个分组；

7)确定所述双向链表每一个分组的形心，所述形心即为其对应的分组的码字；根据码字的个数，将所有码字重新排列为一个体数据codebook；

8)为所述双向链表每一个分组中的输入矢量赋给代表所述输入矢量的码字在码书一或码书二中的索引，计算所述分组对应码字在codebook中的索引，并建立索引体数据，其中码书一对应索引数据体一，码书二对应索引数据体二，两个索引数据体中均包含输入矢量对应的码字在码书一或码书二中的索引信息，索引数据体二中还包括了每一个输入矢量的对应的数据块或小数据块的类别信息，即零块或非零块；

9)确定所述最高层数据或次高层数据中每一个数据块或小数据块在所述码书一或码书二中的索引值信息，根据索引数据体二中的数据块或小数据块的类别信息分别处理非零块和零块：对于零块，直接用步骤2）中计算的均值代替该数据块的标量值；对于非零块，根据索引信息从码书一中查找最高层数据对应的码字，再从码书二中查找次高层数据对应的码字，最后将上述两个码字与该非零块对应的均值相加，即为该非零块的标量值；所述标量值的取值范围是0～255；根据标量值的范围，对每一个标量值值赋给设定颜色信息，然后用这组颜色信息生成地震数据体色标纹理；

10)根据上述步骤9）计算的标量值，查找原始地震数据的色表纹理，获取对应的颜色信息，完成绘制。

所述步骤5）中，量化单元V₁的失真度D₁的计算公式如下：

$D_1={\mathrm{\Σ}}_{I\∈V_1}{\mathrm{\λ}}_j\mathrm{\δ}(I,Y_1),$

其中，Y₁是量化单元V₁的形心，即Y1为V₁中所有输入矢量的平均值；δ(I,Y₁)表示针对输入矢量I，在形心为Y₁的情况下的均方误差；λ_j是为每一个数据块或小数据块设定的误差系数，该误差系数是根据该数据块或小数据块是否为重要分块来设定的（根据用户对重要分块的需求程度设定，如果用户认为重要分开非常重要，且只想要重要分块，则可以将误差系数设置为无求大）；表示对量化单元V₁内的所有输入矢量I的均方误差求和。

重要分块的确定过程为：根据给定的层位的空间位置特征信息，计算每一个数据块的数据范围是否在给定的层位上，若数据块的数据范围在给定的层位上，则该数据块为重要分块。

所述步骤6）中，所述双向链表分裂的过程如下：

1)选择所述双向链表中失真度最大的分组ζ_j；

2)计算所述分组ζ_j的自相关协方差矩阵M；

其中，X_i表示分组ζ_j中的第i个矢量，Y_j表示分组ζ_j的形心；ζ_j表示第j个分组；

3)计算上述自相关协方差矩阵M的最大特征向量和该最大特征向量对应的最大特征值λ_max；

4)将所述分组ζ_j分割为左分组ζ_left和右分组ζ_right，分割公式如下分割公式如下：

ζ_left={i∈ζ_j,<(Y_j-X_i),λ_max><0}；ζ_right={i∈ζ_j,<(Y_j-X_i),λ_max>≥0}

5)计算上述左分组ζ_left和右分组ζ_right的形心和失真度；

6)将上述左分组ζ_left和右分组ζ_right加入到所述双向链表中；

7)将步骤6）得到的双向链表作为下一次分裂的对象，重复步骤1）～6），直到最后

得到的双向链表中的分组个数满足用户设定的码书一或码书二码字个数的要求。

与现有技术相比，本发明所具有的有益效果为：本发明的方法在解决同质区域小，数据分布散乱的三维地震数据体绘制上，压缩比高，失真度明显降低，在大量压缩数据的同时，保证了目标区域的绘制精度，减少目标区域的压缩损失，提高非目标区域的压缩率，同时在提高整体数据压缩率的情况下，提高了目标区域的峰值信噪比；本发明引入了层位特征，结合该特征确定目标区域，利用该信息进行体数据的压缩，可以有效的减少压缩后的数据量；在体数据的预处理中，引入数据分级的概念，利于数据压缩；在矢量量化压缩算法中，引入误差系数的概念，通过该误差系数可以有效的降低目标区域的失真度，提高整体的压缩率；本发明的方法在解决海量数据体绘制上有明显的优势，该方法同样适用于同质区域小，数据分布散乱的其余海量数据压缩体绘制。

附图说明

图1为地震数据效果图；

图2为本发明的方法流程图；

图3为地震数据显示图；

图4为层位特征的地震数据示意图；

图5为本发明数据预处理流程图；

图6为码书大小为25478，层位约束范围分别为12,20,30，误差计算系数逻辑无穷大情况下本发明方法与传统的压缩体绘制方法效果图；图6(a)层位范围12、码书个数为25478的无约束效果图；图6(b)层位约束为12、码书个数为25478的有约束效果图；图6(c)层位范围12、码书个数为25478的无约束改变色标效果图；图6(d)层位范围12、码书个数为25478的有约束改变色效果图；图6(d)层位范围20、码书个数为25478的无约束改变色效果图；图6(e)层位范围20、码书个数为25478的有约束改变色效果图；图6(f)层位范围30、码书个数为25478的无约束改变色效果图；图6(g)层位范围30、码书个数为25478的有约束改变色效果图。

具体实施方式

本方法的总体流程如图2所示，主要包括：数据预处理、矢量量化、解压绘制。

本发明结合地震数据特征中的层位特征（如图3和图4所示），约束原始数据待压缩数据量，提出分级的分类分层矢量量化算法。该方法可以将待压缩数据集中在目标区域，剔除了更多的数据冗余，在相同码书大小的情况下，显著地降低了压缩的失真度。

数据预处理流程如图5对数据进行如下处理：

1）原始地震数据分块处理：本发明采用大小为4*4*4的分块方式将原始地震数据分割成数据量为64个点的效果（该分块大小可根据实际数据大小进行调整）。

2）确定分块数据的重要级别：本发明根据用户给定的层位特征信息，通过计算每一个分块数据范围是否在该层位上，确定该分块是否为重要分块。

3）原始地震数据分类处理：计算步骤1）中每个分块的梯度值，根据用户给定的阈值，将低于该阈值的分块标记为零块，并计算其均值。其他块标记为非零块。

4）分块数据下采样处理：本发明通过对每一个数据块中相邻的2*2*2的数据求平均的方式对数据块进行下采样处理，获得大小为2*2*2的小块。

5）小数据块块求均值处理：计算每一个小数据块的均值。

6）差分处理：分别进行数据块块与小数据块，小数据块与均值之间的差分处理。差分的结果矢量集作为后续矢量量化的输入数据。其中前者的差分结果作为最高层数据，经过矢量量化处理后的结果对应为码书一，后者的差分结果作为次高层数据，相应的经过矢量量化后的结果为码书二。

对于梯度值非零的块，需要对次高层和最高层的数据进行矢量量化。矢量量化的过程可以分为三个阶段：码书设计、编码和解码，其中码书设计会对压缩性能产生十分重要的影响。

本发明矢量量化的具体步骤：

1）将数据预处理阶段获得的每一个4*4*4的块差分数据作为一个输入矢量I，构造一个量化单元V₁，该量化单元包含整个层次的块差分数据，即输入矢量I的集合。Y₁是V₁的形心。V₁的失真度可以通过下式计算获得。

$D_1={\mathrm{\&Sigma;}}_{x\&Element;V_1}{\mathrm{\&lambda;}}_j\mathrm{\&delta;}(X,Y_1)$

其中，D₁为V₁的失真度，δ(X,Y₁)表示针对输入矢量X，在形心为Y₁的情况下起均方误差。λ_j是根据数据预处理阶段确定的每一个块的级别用户设定的误差系数。D₁是通过对每一个输入矢量x与形心的均方误差的和来计算得到的。

2）根据步骤一中获得的分组（V₁）构造双向链表，链表中每一个分组对应矢量量化后的一个码字，然后遍历该链表，进行分裂，具体分裂过程如下

a)选择双向链表中失真度D_j最大的分组ζ_j；

b)计算该分组的自相关协方差矩阵；

c)计算M的最大特征向量和相对应的最大特征值λ_max；

d)ζ_left={i∈ζ_j,<(Y_j-X_i),λ_max><0}；ζ_right={i∈ζ_j,<(Y_j-X_i),λ_max>≥0}根据上述公式将最大分组ζ_jj进行分割为左和右分组；

e)计算新分组的形心和失真度；

f)将新分组加入到分组链表中；

g)将步骤6）得到的双向链表作为下一次分裂的对象，重复步骤1）～6），直到最后得到的双向链表中的分组个数满足用户设定的码书个数要求；码书个数可以设为无穷大，也可能为无穷小。

3）对步骤二中生成的双向链表进行处理：

a)确定每一个分组的形心，为该分组的码字；

b)根据码字的个数，将所有的码字重新排列为一个体数据codebook；

c)分别对每一个分组中的输入矢计算该分组对应码字在codebook中的索引，并建立索引体数据，存储每一个输入矢量对应的在码书中的索引；

4）对次高层2*2*2的差分矢量重复步骤2和3的操作；

最终获得的数据如图5所示分别为最高层索引+均值，次高层索引+均值，最高层码书，次高层码书。

本发明的解压算法如下：对每一个块分别确定该块在码书一和码书二中的索引值信息，然后根据码书二索引中的类别信息分别处理非零值块和零值块。最后根据原始数据空间中数据点在其所在块中的偏移量，计算该点在码字中的偏移量，进而确定数据点的标量值。

结合GPU纹理映射的具体解压算法如下：

Step1：根据纹理坐标从索引纹理中取出码书一和码书二纹理索引；

Step2：从码书二索引中取出块的类别信息，判断如果当前块的是零值块，则直接取均值作为整个块的数据。否则转到Step3；

Step3：根据两份索引信息取出对应次高层和最高层的对应块的索引；

Step4：计算空间点在块内偏移量，并结合块的索引，计算当前点对应到码书一和码书二中的值。

Step5：根据取出的标量值查找色表，获取当前点的颜色和透明信息。

本发明的应用测试均是在Windows7，2.6GHz双核CPU，4G内存，英伟达GeForce210下进行的。编程环境为VS2008，编程语言为C++，OpenGL。压缩后图像的恢复质量使用峰值信噪比来衡量，压缩比使用原始数据与压缩后的数据的商来衡量，另外针对本发明方法的特殊性，为了分析本发明方法性能的优劣，本发明在仿真结果中加入了索引空间的计算，重要块比重等信息。仿真使用的地震数据的范围为xline方向0～400，inline方向为0～255，time方向为0～1751，标量数据类型为char型，故整个仿真数据大小为400×250×1751×8bit。矢量量化过程中码书的大小根据目标特征的实际情况而定，但在仿真中至少保证对每组效果对比图其码书大小是确定的。

对比图6（a）和图6（b）可见，本发明方法在目标区域（图6（a）和图6（b）中红色框区域）的体绘制效果较传统的压缩体绘制方法有明显的改善，但在非目标区域内体绘制效果较差。这是由于该方法在传统算法的基础之上引入了目标特征的概念，通过深入分析三维标量场数据来获取目标区域的特征，进而结合该特征对数据进行矢量量化压缩。这种的方法的优势在于可以有效的使码书向目标区域集中，在保证一定压缩比的前提下，使得目标区域绘制更精细。图6（e）和图6（f）、图6（g）和图6（h）通过同步改变传递函数得到的两组基于目标特征和传统矢量量化压缩体绘制算法的效果对比图，这两组效果图更加突出显示了本发明方法在改变目标区域绘制效果的优势。

表1数据分布情况（空间大小（MB））

本发明中，重要块是指当前码书中重要块数据所占的空间，该数据大小与误差计算系数相关，误差计算系数反应了对重要块精度的要求，表1误差计算系数是逻辑无穷大，结合重要块比重可以看出，在当前误差系数下，重要块几乎占据了码书的全部空间。其中少量非重要块是因为处理上需要而引入的。在一定码书大小下，可以使重要信息无失真。

压缩比数据：表1中压缩比数据分为两部分，第一部分是只以码书数据作为压缩后数据来进行的压缩比计算。在这种情况下，结合图2可以看出，基于特征约束的压缩算法只需要很少量的空间就可以达到很不错的效果。其中码书数据只需要原始数据的1/95.4544。压缩数据的第二部分是包含码书和索引两部分来进行计算。由表1可知，该压缩比骤降，其原因是因为索引数据占据了大量的空间。而索引数据与数据块数成正比，数据块数越大，索引数据所占空间越大。对同样大小的原始数据其分块数取决于块的大小，块空间越大，块数越少。但如果块空间过大，会使恢复的图像出现明显的方块状，图像失真明显。根据本论文仿真情况分析，4×4×4的块在本仿真环境下基本合适。所以通过适当增加各维的维度，如将块大小划分为16×16×16，在当前块数的情况下，可以处理原始数据变为原始数据的64倍（大约10G），此时并没有增加索引数据存储空间，索引数据对压缩比的影响明显降低。所以本发明方法在处理大数据时具有很好的效果。

峰值信噪比（PSNR）：在计算特征约束时，本发明仍然在对数据整体峰值信噪比进行衡量，这种衡量方式只可以显示数据整体的失真度，但在本发明提出的基于特征约束的矢量量化压缩算法前提下，该标准并不能体现出算法的实际恢复效果。所以本发明在此基础之上增加了重要数据峰值信噪（Important Data Peak Signal To Noise Ratio IDPSNR）。依据体绘制的初衷，即展示目标区域的体细节，IDPSNR能更好的体现数据失真情况，IDPSNR越大，说明目标区域信息失真越小，展示给人们的体细节越细腻，体绘制效果越好。如表所示在约束范围为12，码书个数为25478的情况下，重要块个数小于码书个数，此时重要块无失真，因此IDPSNR为无穷大。

Claims (5)

1.一种三维地震体矢量量化压缩体绘制方法，其特征在于，该方法为：

1)将原始地震数据分割成大小为N*N*N大小的数据块；

2)计算每一个数据块的梯度值，根据给定的阈值，将梯度值低于所述阈值的数据块标记为零块，并计算所述零块的均值；将梯度值不低于所述阈值的数据块标记为非零块；所述阈值的绝对值小于100；

3)对每一个数据块中相邻的（N/2）*（N/2）*（N/2）的数据求平均值，获得大小为（N/2）*（N/2）*（N/2）的小数据块；计算每一个小数据块的均值；

4)对所述数据块与所述小数据块进行差分处理，得到最高层数据，即大小为N*N*N的块差分数据；对所述小数据块与所述小数据块的均值进行差分处理，得到次高层数据，即大小为（N/2）*（N/2）*（N/2）的差分矢量；将所述最高层数据经矢量量化后，得到码书一；将所述次高层数据经矢量量化后，得到码书二；

5)将上述步骤4）处理后得到的N*N*N的块差分数据或（N/2）*(N/2)*(N/2)的差分矢量作为输入矢量I，将所述输入矢量I的集合作为一个量化单元V₁，计算量化单元V₁的失真度D₁，进入步骤6）；V₁的形心矢量Y₁就是块差分数据对应的码书一的一个码字或差分矢量对应的码书二的一个码字；

6)根据上述量化单元V₁构造双向链表，所述双向链表中每一个分组对应步骤5）的一个码字，遍历所述双向链表，将所述双向链表分裂为若干个分组；

7)确定所述双向链表每一个分组的形心，所述形心即为其对应的分组的码字；根据码字的个数，将所有码字重新排列为一个体数据codebook；

8)为所述双向链表每一个分组中的输入矢量赋给代表所述输入矢量的码字在码书一或码书二中的索引，计算所述分组对应码字在codebook中的索引，并建立索引体数据，其中码书一对应索引数据体一，码书二对应索引数据体二，两个索引数据体中均包含输入矢量对应的码字在码书一或码书二中的索引信息，索引数据体二中还包括了每一个输入矢量的对应的数据块或小数据块的类别信息，即零块或非零块；

9)确定所述最高层数据或次高层数据中每一个数据块或小数据块在所述码书一或码书二中的索引值信息，根据索引数据体二中的数据块或小数据块的类别信息分别处理非零块和零块：对于零块，直接用步骤2）中计算的均值代替该数据块的标量值；对于非零块，根据索引信息从码书一中查找最高层数据对应的码字，再从码书二中查找次高层数据对应的码字，最后将上述两个码字与该非零块对应的均值相加，即为该非零块的标量值；所述标量值的取值范围是0～255；根据标量值的范围，对每一个标量值值赋给设定颜色信息，然后用这组颜色信息生成地震数据体色标纹理；

10)根据上述步骤9）计算的标量值，查找原始地震数据的色表纹理，获取对应的颜色信息，完成绘制。

2.根据权利要求1所述的三维地震体矢量量化压缩体绘制方法，其特征在于，所述步骤1）中，所述数据块的大小为4*4*4。

3.根据权利要求1所述的三维地震体矢量量化压缩体绘制方法，其特征在于，所述步骤2）中，所述给定的阈值的取值范围为0～20。

4.根据权利要求1～3之一所述的三维地震体矢量量化压缩体绘制方法，其特征在于，所述步骤5）中，量化单元V₁的失真度D₁的计算公式如下：

$D_1={\mathrm{\&Sigma;}}_{I\&Element;V_1}{\mathrm{\&lambda;}}_j\mathrm{\&delta;}(I,Y_1),$

其中，Y₁是量化单元V₁的形心；δ(I,Y₁)表示输入矢量I在形心为Y₁的情况下的均方误差；λ_j是为每一个数据块或小数据块设定的误差系数；表示对量化单元V₁内的所有输入矢量I的均方误差求和。

5.根据权利要求1～3之一所述的三维地震体矢量量化压缩体绘制方法，其特征在于，所述步骤6）中，所述双向链表分裂的过程如下：

1)选择所述双向链表中失真度最大的分组ζ_j；

2)计算所述分组ζ_j的自相关协方差矩阵M；

其中，X_i表示分组ζ_j中的第i个矢量，Y_j表示分组ζ_j的形心；

3)计算上述自相关协方差矩阵M的最大特征向量和该最大特征向量对应的最大特征值λ_max；

4)将所述分组j分割为左分组ζ_left和右分组ζ_right，分割公式如下：

ζ_left={i∈ζ_j,<(Y_j-X_i),λ_max><0}；ζ_right={i∈ζ_j,<(Y_j-X_i),λ_max>≥0}，

其中，<>表示计算向量内积的操作；

5)计算上述左分组ζ_left和右分组ζ_right的形心和失真度；

6)将上述左分组ζ_left和右分组ζ_right加入到所述双向链表中；

7)将步骤6）得到的双向链表作为下一次分裂的对象，重复步骤1）～6），直到最后得到的双向链表中的分组个数满足用户设定的码书一或码书二码字个数的要求。

Images