CN103295208A - 面向地质体数据可视化的保特征量化方法 - Google Patents

Abstract

本发明公开一种面向地质体数据可视化的保特征量化方法，引入双边滤波方法对初始地质体数据做平滑处理，分别统计初始地质体数据及平滑后地质体数据的直方图并得到这两个直方图的差异直方图，利用高斯拟合方法获取高频细节特征区间；定义局部相关性计算方法，统计初始地质体数据的局部相关性直方图，利用高斯拟合方法获取地质体数据中局部相关性较强的特征区间；除去初始地质体数据的高动态范围区间中的高频细节特征区间和局部相关特征区间，获得初始地质体数据中的奇异点特征区间。用户根据自身需求交互定义上述特征的重要性比例，利用分段量化方案，将初始地质体数据由32位高动态范围量化至8位地质体数据，实现初始地质体数据的保特征量化。

Description

面向地质体数据可视化的保特征量化方法

技术领域

本发明涉及一种面向地质体数据可视化的保特征量化方法，属于体数据可视化及图像处理技术领域。

背景技术

二维地质数据可视化技术仅能够显示有限的地质切片数据，严重制约了对地质体数据中三维构造变形、泥沙沉积和流体流动等地质现象的理解。体绘制技术是浏览、分析、解释地质体数据的有效工具，不仅能够帮助地质专家直观、真实地观察和理解地质结构与地层学特征，还能够帮助油气预测专家分析油气蕴藏规律，指导钻井设计和油气开发。

但是，实际采集的地质体数据一般存储为32位浮点型标量数据，即为高动态范围体数据。而高动态范围数据可视化主要存在如下两点不足：（1）数据规模较大，相比于8位体数据，32位体数据虽然精度有所提高，但是体数据规模也提升数倍，可视化效率较低；（2）对于高动态范围体数据的分析相对困难，例如直方图统计、地质属性计算等，不利于体数据的进一步分析和处理。因此，将初始高动态范围体数据量化为计算机容易处理的8位体数据是十分必要而且有意义的。

传统的量化方法是通过量化函数将实数值范围分成有限间隔的过程。量化函数将所有落入量化单元的高动态数据量化为单个值。将高动态的数值量化到任意给定的间隔内，这个过程是不可逆的。量化过程有可能将原始高动态范围地质体数据的目标信号和随机噪声量化到同一量化区间。量化后数据较原始数据丢失一些结构细节特征，量化往往被看作有损压缩算法。GerHardt等人对地质体数据量化处理的研究中强调了量化对于地质体数据可视化的重要性（AndréGerhardt,Anselmo Paiva,Ana Elisa Schmidt,etc.Requisites of3-D Seismic DataVolume Rendering[C].GOCAD ENSG Conference3D Modeling of NaturalObjects,1998.），指出地质体数据量化处理，可以参照计算机图形学与图像处理领域的颜色量化处理算法，如常见颜色量化算法—中值切割（Median-cut）量化、密度（Popularity）量化与均匀（Uniform）量化等算法。多数颜色量化算法为三维空间聚类的过程，三个坐标轴代表三个颜色通道，每个点保存原始图像的颜色。最著名的颜色量化算法为Hechbert提出的中值分割算法（Heckbert,P.,Colorimage quantization for frame buffer display[C].SIGGRAPH Comput.Graph.,1982.16:297-307.）。在此之前，颜色量化主要采用密度量化算法，根据数据直方图，分配量化颜色。颜色量化算法还引入抖动技术（Dithering）,通过插值产生一些在调色板无法找到的颜色，使得量化图像的颜色过度比较自然。Ketterer等人提出了空间颜色量化算法（Jens Ketterer,Jan Puzicha,Marcus Held,etc.On spatialquantization of color images[C].Computer Vision,ECCV'98,Lecture Notes inComputer Science,1998,1406:563-577.），通过结合颜色模板生成过程中的抖动技术和简单的人类视觉感知模型，能够以更少的颜色表示更加逼真的图像。调和映射是图像处理领域将高动态范围图像绘制在普通显示器的另外一类方法。Yuan等人提出高动态范围体绘制算法（Yuan,X.,et al.,HDR VolVis:HighDynamic Range Volume Visualization[J].IEEE Transactions on Visualization andComputer Graphics,2006.12:433--445.），采用非线性高分辨率传输函数将高动态范围体数据映射为高动态范围中间图像，然后采用动态调和映射算法，满足普通显示设备绘制要求。

虽然上述算法一定程度上实现了高动态范围数据的量化与可视化分析与处理，但是针对地质体数据的有效量化方法仍然是地质数据可视化与分析领域的研究热点问题。

发明内容

本发明的目的是提供一种面向地质体数据可视化的保特征量化方法。

为实现上述目的，本发明所采取的技术方案是：

本发明面向地质体数据可视化的保特征量化方法包括如下步骤：

（1）利用直方图分析方法获得初始地质体数据中的高频细节特征的区间范围、局部相关特征的区间范围和奇异点特征的区间范围；

（2）用户根据自身需求交互定义高频细节特征、局部相关特征和奇异点特征的重要性比例，利用分段量化策略，将初始地质体数据由32位高动态范围量化至8位地质体数据；

（3）利用直接体绘制方法对所述8位地质体数据进行绘制，获得地质体数据中的感兴趣特征增强的可视化结果。

进一步地，本发明所述直方图分析方法包括如下步骤：

1）利用双边滤波方法对初始地质体数据做平滑处理，分别统计初始地质体数据及平滑后地质体数据的直方图并得到这两个直方图的差异直方图，利用高斯拟合方法对所述差异直方图做高斯拟合处理，获得初始地质体数据中的高频细节特征的区间范围；

2）利用局部相关性计算方法统计初始地质体数据的局部相关直方图，使用高斯拟合方法对所述局部相关直方图做高斯拟合处理，获得初始地质体数据中的局部相关特征的区间范围；

3）除去初始地质体数据的高动态区间范围中的高频细节特征的区间范围和局部相关特征的区间范围，获得初始地质体数据中的奇异点特征的区间范围。

进一步地，本发明的步骤（2）的所述“分段量化策略”为如下方案中的任一种：

方案一：根据用户定义的高频细节特征、局部相关特征和奇异点特征的重要性比例，将所述高频细节特征、局部相关特征和奇异点特征的分段区间利用线性量化方式进行量化，将初始地质体数据由32位高动态范围量化至8位地质体数据；

方案二：以高频细节特征、局部相关特征和奇异点特征的分段区间为横坐标，以用户定义的高频细节特征、局部相关特征和奇异点特征的重要性比例为纵坐标，并以所述横纵坐标为Bézier曲线的控制点确定Bézier曲线方程，再利用Bézier曲线方程将将初始地质体数据由32位高动态范围量化至8位地质体数据；

方案三：根据用户定义的高频细节特征、局部相关特征和奇异点特征的重要性比例，将高频细节特征、局部相关特征和奇异点特征的分段区间利用色调映射技术将初始地质体数据由32位高动态范围量化至8位地质体数据并保留特征对比度。

进一步地，本发明的步骤1）中“利用局部相关性计算方法统计初始地质体数据的局部相关直方图”的方法如下：

对初始地质体数据中的各个体素以及与其相应的26邻域体素的相似性进行分析，交互定义相似性阈值，根据所述相似性阈值分别记录初始地质体数据中与各体素满足相似关系的体素的数目；

遍历初始地质体数据中的每个体素，若初始地质体数据中与当前体素满足相似关系的体素的数目大于用户交互设置的局部相关阈值，则视当前体素为局部相关性较强体素；否则，即认为当前体素为非局部相关性较强体素；统计所有局部相关性较强体素的数目，进而获得初始地质体数据的局部相关直方图。

与现有技术相比，本发明的有益效果是：引入双边滤波方法对初始地质体数据做平滑处理，分别统计初始地质体数据及平滑后地质体数据的直方图并得到这两个直方图的差异直方图；利用高斯拟合方法获取高频细节特征的区间范围；定义局部相关性计算方法，统计初始地质体数据的局部相关性直方图，利用高斯拟合方法获取地质体数据中局部相关性较强的特征区间；除去初始地质体数据的高动态区间范围中的高频细节特征的区间范围和局部相关特征的区间范围，获得初始地质体数据中的奇异点特征的区间范围。用户根据自身需求交互定义高频细节特征、局部相关特征和奇异点特征的重要性比例，利用分段线性量化方案、Bézier曲线量化方案或色调映射量化方案，将初始地质体数据由32位高动态范围量化至8位地质体数据，实现初始地质体数据的保特征量化。相比于传统的线性量化技术，本发明将初始高动态范围的32位地质体数据量化为计算机图形处理器容易处理的8位体数据，有效保留地质体数据中重要的特征信息，有利于地质解释人员对初始地质体数据做进一步的分析与处理。

附图说明

图1是本发明方法的流程示意图；

图2是获得高频细节特征的区间的示意图，其中，（a）是初始32位地质体数据的统计直方图，（b）是初始32位地质体数据经双边滤波方法平滑获得地质体数据的直方图，（c）是初始地质体数据及平滑后地质体数据的直方图并得到这两个直方图的差异直方图；

图3是获得局部相关性较强特征区间的示意图，其中，（a）是初始32位地质体数据的局部相关性较强体素统计直方图，（b）是局部相关性较强体素直方图的平滑直方图，（c）和（d）是利用高斯拟合方法获取局部相关性较强特征的范围区间，根据用户设定的高斯拟合曲线峰值阈值，分别获得了5个和3个局部相关性较强特征区间；

图4是分段线性量化方案示意图；

图5是Bézier曲线量化方案示意图；

图6是色调映射量化方案示意图；

图7是不同量化方案获得数据的可视化结果对比图，其中，（a）是初始32位地质体数据可视化结果；（b）是利用直接线性量化策略获得的8位地质体数据可视化结果；（c）是利用分段线性量化策略获得的8位地质体数据可视化结果；（d）是利用Bézier量化策略获得的8位地质体数据可视化结果。

具体实施方式

下面结合附图，对本发明的面向地质体数据可视化的保特征量化方法作进一步的说明。参见图1，本发明的具体步骤如下：

步骤（1）：对初始的高动态范围地质体数据做直方图统计。由于初始的地质体数据存储类型为32位浮点型数据，为能够有效表达初始地质体数据的分布信息，将每个体素对应的浮点型数据取整，统计于对应的整型区间分段中，图2（a）出示了初始高动态范围地质体数据的直方图。

步骤（2）：利用双边滤波方法对初始地质体数据做平滑处理，统计平滑后体数据的直方图，如图2（b）所示，进而获得初始地质体数据的直方图和平滑后地质体数据的直方图的差异直方图，即地质体数据高频细节特征的分布直方图，如图2（c）所示。

双边滤波方法是一种可以保边去噪的组合式滤波器，主要由两个函数构成，一个是由几何空间距离决定滤波器系数，而另一个是由振幅值差异决定的滤波器系数，如下：

其中，p表示当前体素的空间位置，q表示与当前体素相邻的体素的空间位置，||p-q||表示当前体素与相邻体素的几何空间距离，G_σs表示对应的几何空间距离高斯核函数，控制几何空间距离对平滑尺度的影响；I_p表示当前体素的振幅值，I_q表示与当前体素相邻的体素的振幅值，|I_p-I_q|表示当前体素与相邻体素的振幅值差异，G_σi表示对应的振幅值差异高斯核函数，控制振幅值差异对平滑尺度的影响；

表示相邻体素的取值范围，通常由用户交互指定，如体素周围的26邻域；N_p是单位化权值因子；I′_p是当前体素p经双边滤波器平滑后的振幅值。

步骤（3）：利用高斯曲线拟合方法对差异直方图进行拟合，进而以高斯曲线的均值μ和方差σ确定高频细节特征区间，即为μ±λσ，其中，λ为用户交互指定的参数，在本发明中，λ可取经验值为1.5。

图2（c）中的黑色曲线即为高斯拟合曲线。

步骤（4）：利用局部相关性分析方法，统计当前体素的邻域体素中满足用户交互定义相似性阈值的体素数目。若当前体素邻域中满足相似关系的体素数目大于用户交互定义的数目阈值threhold₂，本发明中，threhold₂可取经验值为18，则视当前体素为局部相关性较强体素。遍历初始地质体数据中所有体素，***部相关性较强体素直方图（如图3（a）所示）。

首先，分析当前体素与其相邻体素的相关性c(p,q)，如下：

$c(p,q)=\left\{\begin{array}{cc}1& |I_p-I_q|\≤{\mathrm{threhold}}_1\\ 0& \mathrm{otherwise}\end{array}\right.$

其中，p表示当前体素的空间位置，q表示与当前体素p相邻的体素。threhold₁表示振幅值差异阈值，即相似性阈值，本发明中，threhold₁可取经验值为0.5。如果邻域体素q与当前体素p的振幅值差异小于用户交互定义的相似性阈值，则c(p,q)=1，否则c(p,q)=0。如果c(x)=1，则表示当前体素p与邻域体素q相关。如果c(x)=0，则表示当前体素p与邻域体素q无关。

进而，分析当前体素的局部相关性r(p)，如下：

其中，

表示在一定邻域范围内，与当前体素p满足相关性的体素数目之和；如果则表示当前体素为局部相关性较强的体素，否则表示当前体素为非局部相关性较强的体素，阈值threhold₂可由用户根据需求交互设置，本发明中，threhold₂可取经验值为18。

遍历初始地质体数据中的所有体素，分析每个体素的局部相关性，进而***部相关性较强体素直方图LRS(s)如下：

$\mathrm{LRS}\left(s\right)=\underset{f\left(p\right)=s}{\mathrm{\Σ}}r\left(p\right)/\left|V\right|$

其中，s表示初始高动态范围地质体数据对应的整数区间值，p表示初始地质体数据中体素的空间位置，f(p)表示初始地质体数据中体素的空间位置p对应的振幅值，|V|表示初始地质体数据的体素数目。遍历初始地质体数据中每一体素，若体素满足局部相关性，则累加于其对应振幅值的整数区间中，最终完成局部相关性较强体素直方图统计。

步骤（5）：利用高斯平滑技术对局部相关性较强体素直方图做平滑处理，如图3（b）所示。同样利用高斯曲线拟合方法对平滑后的局部相关性较强体素直方图做拟合，获得对应的局部相关特征区间为μ±λσ，从全局最大值开始拟合，迭代拟合剩余区间，直到峰值小于交互指定的阈值为止。因此，获得了不同的、连续的局部相关特征区间，如图3（c）和图3（d）所示。

步骤（6）：除去初始地质体数据的高动态区间范围中的高频细节特征的区间范围和局部相关特征的区间范围，获得初始地质体数据中的奇异点特征的区间范围。

步骤（7）：根据上述步骤确定的初始高动态范围中高频细节特征、局部相关性较强特征及奇异点特征区间，用户根据需求交互指定不同特征的重要性比例，即对应的目标量化区间，进而利用分段量化策略将初始的高动态范围地质体数据量化至计算机容易处理和分析的8位地质体数据，具体策略可选择如下三种方案中的任一种：

（a）分段线性量化方案：将不同特征对应的初始高动态特征区间及用户交互指定的目标量化区间的特征点进行连线，对不同区间的数据按照对应的线段做线性量化，可以获得按照用户定义的重要性比例的、重要特征得以保留的8位范围体数据，如图4所示。

（b）Bézier曲线量化方案：将不同特征对应的初始高动态特征区间及用户交互指定的目标量化区间的特征点视为Bézier曲线的控制点，生成Bézier曲线，进而根据Bézier曲线方程将初始高动态范围地质体数据量化至8位体数据，由于控制点的位置是根据交互定义的重要性比例而确定的，因此最终的量化结果保留了感兴趣特征信息，且具有较强的连续性，如图5所示。

（c）色调映射量化方案：为了进一步保留特征的对比度，将不同特征对应的初始高动态特征区间及用户交互指定的目标量化区间，利用色调映射技术进行量化，进一步保留了有效特征信息，如图6所示；

$C_{\mathrm{new}}=(L_{\max }-L_{\min })*\frac{\log (C+\mathrm{\γ}*(C_{\max }-C_{\min }))-\log (H_{\min }+\mathrm{\γ}*(C_{\max }-C_{\min }))}{\log (H_{\max }+\mathrm{\γ}*(C_{\max }-C_{\min }))-\log (H_{\min }+(C_{\max }-C_{\min }))}+L_{\min }$

其中，H_max、H_min分别表示初始高动态范围特征区间的最大值与最小值，而L_max、L_min则分别表示8位量化目标区间的最大值、最小值，γ是由用户交互指定的色调映射系数，C代表高动态范围内某一振幅值，而C_new则代表了目标区间中与其对应的量化结果。

步骤（8）：对保特征量化获得的8位地质体数据进行可视化，采用经典的直接体绘制技术对其进行绘制，并利用经典的层位提取技术验证所述保特征量化方法获取的8位地质数据的有效性，如图7所示。

使用本发明方法后，不仅可以将规模庞大的初始高动态范围的地质体数据量化至计算机容易处理的8位体数据，而且根据用户交互定义的特征重要性比例，确定分段量化策略，以有效保留初始地质体数据的重要特征信息。相比于传统的直接线性量化技术，本发明不仅获得了规模较小、容易分析与处理的8位体数据，而且体数据的内部特征及对比度得到了有效保留，具有较强的实用价值。

Claims (4)

1.一种面向地质体数据可视化的保特征量化方法，其特征在于，包括如下步骤：

（1）利用直方图分析方法获得初始地质体数据中的高频细节特征的区间范围、局部相关特征的区间范围和奇异点特征的区间范围；

（2）用户根据自身需求交互定义高频细节特征、局部相关特征和奇异点特征的重要性比例，利用分段量化策略，将初始地质体数据由32位高动态范围量化至8位地质体数据；

（3）利用直接体绘制方法对所述8位地质体数据进行绘制，获得地质体数据中的感兴趣特征增强的可视化结果。

2.根据权利要求1所述的保特征量化方法，其特征在于，所述直方图分析方法包括如下步骤：

1）利用双边滤波方法对初始地质体数据做平滑处理，分别统计初始地质体数据及平滑后地质体数据的直方图并得到这两个直方图的差异直方图，利用高斯拟合方法对所述差异直方图做高斯拟合处理，获得初始地质体数据中的高频细节特征的区间范围；

2）利用局部相关性计算方法统计初始地质体数据的局部相关直方图，使用高斯拟合方法对所述局部相关直方图做高斯拟合处理，获得初始地质体数据中的局部相关特征的区间范围；

3）除去初始地质体数据的高动态区间范围中的高频细节特征的区间范围和局部相关特征的区间范围，获得初始地质体数据中的奇异点特征的区间范围。

3.根据权利要求1或2所述的保特征量化方法，其特征在于，步骤（2）的所述“分段量化策略”为如下方案中的任一种：

方案一：根据用户定义的高频细节特征、局部相关特征和奇异点特征的重要性比例，将所述高频细节特征、局部相关特征和奇异点特征的分段区间利用线性量化方式进行量化，将初始地质体数据由32位高动态范围量化至8位地质体数据；

方案二：以高频细节特征、局部相关特征和奇异点特征的分段区间为横坐标，以用户定义的高频细节特征、局部相关特征和奇异点特征的重要性比例为纵坐标，并以所述横纵坐标为Bézier 曲线的控制点确定Bézier 曲线方程，再利用Bézier 曲线方程将将初始地质体数据由32位高动态范围量化至8位地质体数据；

方案三：根据用户定义的高频细节特征、局部相关特征和奇异点特征的重要性比例，将高频细节特征、局部相关特征和奇异点特征的分段区间利用色调映射技术将初始地质体数据由32位高动态范围量化至8位地质体数据并保留特征对比度。

4.根据权利要求2所述的保特征量化方法，其特征在于，步骤1）中“利用局部相关性计算方法统计初始地质体数据的局部相关直方图”的方法如下：

对初始地质体数据中的各个体素以及与其相应的26邻域体素的相似性进行分析，交互定义相似性阈值，根据所述相似性阈值分别记录初始地质体数据中与各体素满足相似关系的体素的数目；

遍历初始地质体数据中的每个体素，若初始地质体数据中与当前体素满足相似关系的体素的数目大于用户交互设置的局部相关阈值，则视当前体素为局部相关性较强体素；否则，即认为当前体素为非局部相关性较强体素；统计所有局部相关性较强体素的数目，进而获得初始地质体数据的局部相关直方图。

Images