CN103871064B - 一种火山岩ct图像的预处理和确定分割阈值的方法 - Google Patents

Abstract

本发明公开了一种火山岩CT图像的预处理和确定分割阈值的方法，该火山岩CT图像的预处理和确定分割阈值的方法包括以下步骤：火山岩CT选样；火山岩CT2D图像的采集；2DCT图像数据的预处理；叠加图像形成三维数据体；三维重构体二值化阈值方法的改进及确立；建立火山岩CT图像的三维重构模型，并进行三维显示。本发明克服了以往分析火山岩储集空间的尺度限制的难题，拓宽了CT技术的应用范围，使得广大油田科研人员可以更方便地定量进行火山岩CT图像的三维重构，火山岩CT图像的预处理和确定分割阈值方法具有科学性和普适性。

Description

一种火山岩CT图像的预处理和确定分割阈值的方法

技术领域

本发明属于火山岩储集特征研究技术领域，尤其涉及一种火山岩CT图像的预处理和确定分割阈值的方法。

背景技术

自从伦琴1895年发现X射线以来，CT技术发展迅速，这以后的一百多年，发展了许多成像技术，如计算机断层成像技术（CT）、核磁共振成像技术（MRI），正电子发射断层成像技术（PET），单光子发射断层成像技术（SPECT），超声成像。

计算机断层成像技术（CT）是一种无损伤检测技术，能检测出岩石材料内部结构的变化，已经成为观测岩石类材料内部结构变化和裂隙演化过程重要实验手段之一。Raynaud（2008）、杨更社（1999）、葛修润（1999）等研究了岩石材料受力情况与裂纹结构、CT数的关系，并建立了相关的理论与公式，并研究了岩石内部微裂纹演化与CT数间的关系；康志勤（2009）等进行了油页岩从常温到600℃高温下的热破裂过程显微观测和分析，确定了原油的热破裂阈值温度；杨栋（2011）等研究了不同地区油页岩在100℃和600℃下裂缝变化情况，揭示了不同地区油页岩成因不同，在高温下会形成不同孔隙连通模式；郭和坤（2010）等研究了岩石样品的CT图像读取、显示控制、一维分析、二维分析等问题；周枫（2013）等将CT技术和显微镜下岩石薄片鉴定技术结合起来，通过建立鉴定的岩心样品的矿物成分和含量，与三维灰度图像之间的关系，实现研究岩心的各项异性特征。然而国内没有专用于岩石CT图像预处理和确定分割阈值的方法，只能借助医用CT的方法，而医用CT图像的预处理和确定分割阈值有其临床的独特性，不能满足岩石的CT图像的非均质性，不适用岩石样品分析，而国内关于岩石CT图像的预处理和确定灰度阈值方面研究工作还鲜见报道。

现有的分析火山岩储集空间的尺度受限制。

发明内容

本发明实施例的目的在于提供一种火山岩CT图像的预处理和确定分割阈值的方法，旨在解决现有的分析火山岩储集空间的尺度受限制的问题。

本发明实施例是这样实现的，一种火山岩CT图像的预处理和确定分割阈值的方法，该火山岩CT图像的预处理和确定分割阈值的方法包括以下步骤：

步骤一，综合考虑反映火山岩岩石的特征，选择具有代表性的岩样，从喷溢相、爆发相、浸入相、火山通道相和火山沉积相五种岩相中，选出与油气分布关系较大的岩相/岩性，做CT扫描实验，工业CT、微米CT或纳米CT；

步骤二，火山岩CT2D图像的采集，按完全投影数据方式进行图像序列采集，获得CT三维重构的一手2DCT图像资料；

步骤三，2DCT图像数据的预处理，将较差的原始图像变的清晰、信息量 丰富的可用图像，有效地去掉图像中的噪声、突出图像中边缘或感兴趣的区域，引进了邻域内像素的方差作为两种方法的衔接纽带，既消除了噪声的影响又增强了图像的边缘；

步骤四，叠加图像形成三维数据体；以预处理过的2DCT图像为支撑点，利用Amira软件建立三维数据体，为火山岩CT三维重构奠定基础；

步骤五，三维重构体二值化阈值方法的改进及确立；综合学习了解阈值分割方法、边缘分割方法和区域分割方法，进行改进提升，重新修订Ostu法，建立了火山岩CT图像三维重构的二值化阈值分割方法；

步骤六，建立火山岩CT图像的三维重构模型，并进行三维显示；基于步骤三获得预处理后的火山岩2DCT图像，进一步结合步骤四、步骤五，实现火山岩CT图像二值化的三维重构，在利用Amira软件进行体显示。

进一步，在步骤三中，图像预处理将均值滤波和中值滤波两种方法进行结合，重新考虑了两组的自身相似性和两组间的相似性的关系，并建立相应公式。

进一步，在步骤三中，CT图像预处理是对CT图像进一步处理与分析的先行步骤，设M是3×3的正方形邻域，当正方形邻域内的方差大于等于ε_临界值（ε _临界值=18）时，用中值滤波进行图像预处理，当正方形邻域内的方差小于ε_临界值（ε _临界值=18）时，用均值滤波进行预处理。

进一步，均值滤波的预处理方法为：对于给定的图像f（x,y）中的每一个点（m,n），取邻域A，设A含有N个像素，取平均值作为预处理后所得图像像素点（m,n）处的灰度，设M是3×3的正方形邻域，点（m,n）位于M中心，具体公式如下：

中值滤波的预处理方法为：给定图像f（x,y）中的每一个点（m,n）的邻域A，设A含有N个像素{a₁,a₂,┄,a_N}，将其按像素大小排序，若N是奇数时， 则位于中间的那个像素值就是修改后图像g(x,y)在f（m,n）处的像素值，若N是偶数时，则取中间两个像素的平均值就是修改后图像g(x,y)在f（m,n）处的像素值，具体表达式如下：

进一步，建立火山岩CT图像预处理的方法，具体公式如下

建立相关函数，当Z值最大，确定火山岩CT三维重构体二值化的最优阈值。

进一步，灰度图像的具处理公式如下：

进一步，在步骤五中阈值分割法基于Otsu法进行改进提升，方法为：

给定一幅火山岩CT图像灰度级的取值范围为G={0，1，2┄，L-1，}，大小为M×N的图像，假设灰度级为i的像素点数为m_i，各灰度级的分布概率为：

$P_i=\frac{m_i}{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{m}_i},i=0,\·\·\·\·\·\·,L-1$

设定一个阈值T，将图像灰度级分为C₀和C₁两组，u₀和ω₀为C₀组的灰度平均值和产生的概率；u₁和ω₁为C₁组的灰度平均值和产生的概率，两组间的方差表达式如下式：

$\mathrm{\δ}\left(t\right)=\sqrt{{\mathrm{\ω}}_0*{(u_2-u)}^2+{\mathrm{\ω}}_1*{(u_1-u)}^2}$

其中u是整幅图像的灰度平均值，如下式：

${\mathrm{\ω}}_0{u}_0+{\mathrm{\ω}}_1{u}_1={\mathrm{\ω}}_0\underset{i=0}{\overset{T}{\mathrm{\Σ}}}\frac{{\mathrm{iP}}_i}{{\mathrm{\ω}}_0}+{\mathrm{\ω}}_1\underset{i=T+1}{\overset{L-1}{\mathrm{\Σ}}}\frac{{\mathrm{iP}}_i}{{\mathrm{\ω}}_1}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{iP}}_i=u$

进而将两组间的方差和每一组内的方差都考虑进来，进行改进提升，建立了火山岩CT图像的阈值分割方法，得到公式，如下：

$Z=\frac{\mathrm{\δ}\left(t\right)}{{\mathrm{\δ}}_{u_0}\left(t\right)+{\mathrm{\δ}}_{u_1}\left(t\right)}=\frac{\sqrt{{\mathrm{\ω}}_0*{(u_2-u)}^2+{\mathrm{\ω}}_1*{(u_1-u)}^2}}{\sqrt{\underset{0}{\overset{T}{\mathrm{\Σ}}}{({\mathrm{\ω}}_0-P_i)}^2}+\sqrt{\underset{T+1}{\overset{L}{\mathrm{\Σ}}}{({\mathrm{\ω}}_1-P_i)}^2}}$

当Z值取最大时，此时，两组间的相似性最小，方差值最大，而其中的每一组组内方差值最小，找到了最优化的阈值T，此时就可取得最优阈值T：

T＝argmaxσ²(t)

基于上述确定的阈值，利用Amira软件进行火山岩CT图像三维重构体的二值化。

本发明提供的火山岩CT图像的预处理和确定分割阈值的方法，采用适合火山岩CT图像预处理方法和确定三维重构体二值化的阈值分割方法，并建立相应的数学表达式，克服了以往分析火山岩储集空间的尺度限制的难题，拓宽了CT技术的应用范围，使得广大油田科研人员可以更方便地定量进行火山岩CT图像的三维重构，火山岩CT图像的预处理和确定分割阈值方法具有科学性和普适性。本发明实现了火山岩CT图像的三维重构，给出了火山岩CT图像的预处理和确定阈值分割的方法，提出了火山岩CT图像三维重构（主要是火山岩CT图像的预处理和确定阈值分割方法）的具体流程，为油田的实际火山岩储集空间评价工作提供了技术服务支持。本发明具有操作可行，计算简单方便的特点；实现了火山岩空隙结构的三维可视化，与火山岩物性分析数据相吻合，实现了火山岩CT图像的三维重构，提取切割出重建后的2D剖面，进行火山岩空隙直径统计，统计火山岩CT岩样的空隙直径大小分布规律，为今后的油气 勘探指明方向。本发明有效地将中值滤波和均值滤波结合起来，有效地去除噪声，提高了岩石CT图像的清晰度，为火山岩CT图像预处理提供了理论基础。

附图说明

图1是本发明实施例提供的火山岩CT图像的预处理和确定分割阈值的方法流程图；

图2是本发明实施例提供的松辽盆地徐家围子断陷火山岩CT样品选取井位示意图；

图3是本发明实施例提供的松辽盆地徐家围子断陷徐深5井火山岩选取剖面位置图；

图4是本发明实施例提供的松辽盆地徐家围子断陷徐深12井火山岩选取剖面位置图；

图5是本发明实施例提供的松辽盆地徐家围子断陷徐深8井火山岩选取剖面位置图；

图6是本发明实施例提供的松辽盆地徐家围子断陷汪深101井火山岩选取剖面位置图；

图7是本发明实施例提供的松辽盆地徐深8井火山岩CT图像预处理前后对比分析图

图8是本发明实施例提供的灰度值重新划分图；

图9是本发明实施例提供的松辽盆地徐家围子断陷徐深8井火山岩二值化的2D剖面图；

图10是本发明实施例提供的松辽盆地徐家围子断陷徐深12井流纹岩工业 CT三维重建结果展示图；

图11是本发明实施例提供的松辽盆地徐家围子断陷汪深101井安山岩工业CT三维重建结果展示图；

图12是本发明实施例提供的松辽盆地徐家围子断陷徐深8井火山岩角砾岩工业CT三维重建结果展示图；

图13是本发明实施例提供的松辽盆地徐家围子断陷徐深5井凝灰岩工业CT三维重建结果展示图；

图14是本发明实施例提供的松辽盆地徐家围子断陷徐深12井流纹岩工业CT二维结果展示图；

图15是本发明实施例提供的松辽盆地徐家围子断陷徐深8井火山岩角砾岩工业CT二维结果展示图；

图16是本发明实施例提供的松辽盆地徐家围子断陷汪深101井安山岩工业CT二维结果展示图；

图17是本发明实施例提供的松辽盆地徐家围子断陷徐深5井凝灰岩工业CT二维结果展示图；

图18是本发明实施例提供的松辽盆地徐家围子断陷徐深8井火山岩角砾岩三维重构结果与铸体薄片对比分析图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图及具体实施例对本发明的应用原理作进一步描述。

如图1所示，本发明实施例的火山岩CT图像的预处理和确定分割阈值的方法包括以下步骤：

S101：火山岩CT选样：综合考虑能反映各种火山岩岩石特征，选择具有代表性的岩样，从喷溢相、爆发相、浸入相、火山通道相和火山沉积相五种岩相中，选出与油气分布关系较大的岩相/岩性，做CT扫描实验，工业CT、微米CT或纳米CT；

S102：火山岩CT2D图像的采集，主要针对三个方向进行，按完全投影数据方式进行图像序列采集，获得CT三维重构的一手2DCT图像资料；

S103：2DCT图像数据的预处理，将较差的原始图像变的清晰、信息量丰富的可用图像，有效地去掉图像中的噪声、突出图像中边缘或感兴趣的区域，引进了邻域内像素的方差作为两种方法的衔接纽带，既消除了噪声的影响又增强了图像的边缘；

S104：叠加图像形成三维数据体；以:预处理过的2DCT图像为支撑点，利用Amira软件建立三维数据体，为火山岩CT三维重构奠定基础；

S105：三维重构体二值化阈值方法的改进及确立；综合学习了解阈值分割方法、边缘分割方法和区域分割方法，进行改进提升，重新修订Ostu法，建立了适合火山岩CT图像三维重构的二值化阈值分割方法；

S106：建立火山岩CT图像的三维重构模型，并进行三维显示。

在步骤S103中，图像预处理将均值滤波和中值滤波两种方法进行结合，克服了噪声、较差的原始图像的影响，将火山岩CT2D原始灰度图像处理成清晰、视觉效果好的可用图像，确定火山岩CT图像二值化阈值分割法是基于阈值分割方法、边缘分割方法和区域分割方法的研究，进行改进提升，建立一种适合 火山岩CT图像二值化的阈值分割方法，该方法主要是依据Ostu法（阈值分割方法中的一种），但重新考虑了两组的自身相似性和两组间的相似性的关系，并建立相应公式，当Z值最大时，确定了二值化的最优阈值，为火山岩CT三维重构中半自动识别空隙（孔隙和裂缝）和基质奠定坚实的基础，上述几步确保火山岩CT图像的预处理和确定分割阈值能顺利进行，保证了火山岩CT图像三维重构结果合理、可靠；

CT图像三维重建的流程中，预处理处于三维重建的首要环节，是对CT图像进一步处理与分析的先行步骤，因此预处理效果的好坏直接影响到CT图像三维重建的效果，要对CT图像预处理和分析，必须对降质失真的图像进行必要的预处理，达到改善图像质量的目的，因此一种高效的预处理算法将大大提高后处理步骤的效率和质量，对火山岩图像三维重建的实现具有重要的现实意义；

将两种滤波方法结合，汲取两种方法的优点，既消除了噪声的影响又实现了增强图像质量的目的，总体改进思想：设M是3×3的正方形邻域，当正方形邻域内的方差大于等于ε_临界值（ε_临界值=18）时，用中值滤波进行图像预处理，当正方形邻域内的方差小于ε_临界值（ε_临界值=18）时，用均值滤波进行预处理；

均值滤波的预处理原理如下：对于给定的图像f（x,y）中的每一个点（m,n），取其邻域A，设A含有N个像素，取其平均值作为预处理后所得图像像素点（m,n）处的灰度，设M是3×3的正方形邻域，点（m,n）位于M中心，具体公式如下：

中值滤波的预处理原理如下：给定图像f（x,y）中的每一个点（m,n）的邻域A，设A含有N个像素{a₁,a₂,┄,a_N}，将其按像素大小排序，若N是奇数时，则位于中间的那个像素值就是修改后图像g(x,y)在f（m,n）处的像素值，若N是偶数时，则取中间两个像素的平均值就是修改后图像g(x,y)在f（m,n）处的像素值，具体表达式（2）如下：

基于上述均值滤波和中值滤波两种图像预处理算法进行改进提升，主要体现为将这两种滤波方法结合起来（这两种方法的结合是基于选定邻域M内所有像素的方差实现的），建立适合火山岩CT图像预处理的方法，具体公式（3）如下

改进提升后的算法主要表现在既有计算速度快的特点，又能滤去图像中噪声影响和增强图像的边缘和感兴趣区域，具有很好的处理效果，为图像分割奠定了坚实的基础；

确定火山岩CT图像三维重构体二值化阈值分割法是基于阈值分割方法、边缘分割方法和区域分割方法的研究，进行改进提升，建立一种确定火山岩CT图像三维重构体二值化的阈值分割方法，该方法主要是依据Ostu法在两组的自身相似性和两组间的相似性的关系进行改进提升，建立相关函数，当Z值最大，确定火山岩CT三维重构体二值化的最优阈值，为火山岩CT三维重构中半自动识别空隙（孔隙和裂缝）和基质奠定坚实的基础；

图像分割技术就是将目标区域与背景区域的分离，精确的分割才能实现最大程度的还原，才能更有效地研究火山岩储层内部的储集空间分布特征（空隙大小及其连通性），可以更逼真地了解地下火山岩储集空间分布情况，然而，关于岩石CT图像预处理和确定分割阈值方面研究还鲜见报道，这成为制约岩石CT图像技术发展的一大难题，因而，本次火山岩三维重构探索了图像阈值分割方法，并进行改进及建立了适合火山岩CT图像三维重构体二值化的阈值分割方法

阈值分割方法的原理主要利用目标物体与背景在灰度上的差异，将图像分 割为不同灰度范围的背景区域和目标区域，本次阈值分割法基于Otsu法进行改进提升，原理如下：

给定一幅火山岩CT图像灰度级的取值范围为G={0，1，2┄，L-1，}，大小为M×N的图像，假设灰度级为i的像素点数为m_i，各灰度级的分布概率（式4）为：

$P_i=\frac{m_i}{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{m}_i},i=0,\·\·\·\·\·\·,L-1---\left(4\right)$

设定一个阈值T，将图像灰度级分为C₀和C₁两组，u₀和ω₀为C₀组的灰度平均值和产生的概率；u₁和ω₁为C₁组的灰度平均值和产生的概率，两组间的方差表达式如式（5）：

$\mathrm{\δ}\left(t\right)=\sqrt{{\mathrm{\ω}}_0*{(u_2-u)}^2+{\mathrm{\ω}}_1*{(u_1-u)}^2}---\left(5\right)$

其中u是整幅图像的灰度平均值，如式（6）

${\mathrm{\ω}}_0{u}_0+{\mathrm{\ω}}_1{u}_1={\mathrm{\ω}}_0\underset{i=0}{\overset{T}{\mathrm{\Σ}}}\frac{{\mathrm{iP}}_i}{{\mathrm{\ω}}_0}+{\mathrm{\ω}}_1\underset{i=T+1}{\overset{L-1}{\mathrm{\Σ}}}\frac{{\mathrm{iP}}_i}{{\mathrm{\ω}}_1}=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{iP}}_i=u---\left(6\right)$

进而将两组间的方差和其每一组内的方差都考虑进来，进行改进提升，建立了适合火山岩CT图像的阈值分割方法，得到公式（7），如下：

$Z=\frac{\mathrm{\δ}\left(t\right)}{{\mathrm{\δ}}_{u_0}\left(t\right)+{\mathrm{\δ}}_{u_1}\left(t\right)}=\frac{\sqrt{{\mathrm{\ω}}_0*{(u_2-u)}^2+{\mathrm{\ω}}_1*{(u_1-u)}^2}}{\sqrt{\underset{0}{\overset{T}{\mathrm{\Σ}}}{({\mathrm{\ω}}_0-P_i)}^2}+\sqrt{\underset{T+1}{\overset{L}{\mathrm{\Σ}}}{({\mathrm{\ω}}_1-P_i)}^2}}---\left(7\right)$

当Z值取最大时，此时，两组间的相似性最小，方差值最大，而其中的每一组组内方差值最小，找到了最优化的阈值T，此时就可取得最优阈值T：

T＝argmaxσ²(t) （8）

基于上述确定的阈值，利用法国TGS公司的Amira软件进行火山岩CT图像三维重构体的二值化。

结合本发明的具体实施方式对本发明做进一步的说明：

实施例1：以松辽盆地盆地徐家围子断陷深层火山岩营城组为例，使用火 山岩CT图像的预处理方法和三维重构体二值化的阈值分割方法，实现火山岩CT图像三维重构体的二值化，具体步骤为：

（1）火山岩CT选样：综合考虑徐家围子断陷内能反映各种火山岩岩石特征，选择具有代表性的岩样（从喷溢相、爆发相、浸入相、火山通道相和火山沉积相五种岩相中，选出与油气分布关系较大的岩相/岩性），在本次研究区内具有代表性的岩样有火山角砾岩、流纹岩、安山岩和凝灰岩，因而选取上述四种岩石样品做工业CT扫描实验，如表1；井位分布位置如附图2，从图中可以看出，徐深5井、徐深8井和徐深12井分布在徐中地区，汪深101井分布在安达地区，附图3至附图6详细地阐明了每个岩样选自该口井位置的详细信息，选取四块火山岩样品（徐深8、汪深101、徐深5、徐深12分别为喷溢相夹爆发相、溢流相、火山通道相夹爆发相、侵入相）进行CT实验，四块样品的详细信息如下：

徐深5井样品所在地层为白垩系营城组营一段，样品深度为3671.65m，岩相为喷溢相，亚相为中部亚相和上部亚相，长度为3cm左右，颜色灰白，为凝灰岩（附图3）；

徐深12井样品深度为3667.21m，岩相为侵入相，长度为2cm左右，颜色发灰，为流纹岩（附图4），

徐深8井样品所在地层为白垩系营城组营一段，样品深度为3714.65m，岩相为喷溢相夹爆发相，亚相为上部亚相夹热碎屑流亚相，长度为2cm左右，颜色灰白，表面多孔，为火山角砾岩（附图5）；

汪深101井样品所在地层为白垩系营城组营三段，样品深度为3243.82m，岩相为溢流相，亚相为上部亚相，长度为2cm左右，颜色发暗，表面较为平整，为安山岩（附图6）；

表1不同井不同岩性工业CT岩石样品信息统计表

样品序号	井名	CT类型	深度	岩性
					1	徐深8井	工业CT	3714.65m	火山角砾岩

2	汪深101井	工业CT	3243.82m	安山岩
					3	徐深5井	工业CT	3671.65m	凝灰岩
4	徐深11井	工业CT	3667.21m	流纹岩

（2)火山岩CT2D图像的采集，主要针对三个方向（XY、XZ和YZ）进行，按完全投影数据方式进行图像序列采集，获得CT三维重构的一手2DCT图像资料；

（3)2DCT图像数据的预处理；图像的预处理主要体现在图像增强，图像增强就是将较差的原始图像变的清晰、信息量丰富的可用图像，有效地去掉图像中的噪声、突出图像中边缘或感兴趣的区域，以徐深8井的一张2DCT图像为例，进行阐释说明（图7），因而，将均值滤波和中值滤波两种方法结合起来，引进了邻域内像素的方差作为两种方法的衔接纽带，既消除了噪声的影响又增强了图像的边缘；如附图8和附图9中，可以看出，经过火山岩CT图像预处理的图像，噪声的影响基本全部消失，而且增强了图像的边缘清晰效果；

（4)叠加图像形成三维数据体；以C)步骤中预处理过的2DCT图像为支撑点，利用Amira软件建立三维数据体，为火山岩CT三维重构奠定基础；

（5）三维重构体二值化阈值方法的改进及确立；综合学习了解阈值分割方法、边缘分割方法和区域分割方法，进行改进提升，重新修订Ostu法（阈值分割法中的一种），建立了适合火山岩CT图像三维重构的二值化阈值分割方法；

（6)建立火山岩CT图像的三维重构模型，并进行三维显示(voltex或volren)，如附图10至附图17

从附图10至附图13中可以看出，徐深12井流纹岩、徐深8井火山角砾岩孔隙发育，汪深101井安山岩孔隙较发育，徐深5井凝灰岩孔隙不发育，与勘探实际相吻合；

从附图14至附图17实现了不同火山岩类型孔隙、裂缝分布对比，从图中可以看出，流纹岩空隙（孔隙、喉道）较角砾岩空隙发育，安山岩空隙较为发育，凝灰岩空隙不发育；

此外，进行了火山岩角砾岩的孔隙结构大小的统计，发现徐深8井火山角砾岩工业CT孔隙最大值800um，平均值110um，最小值20um，同时，重构出来的孔隙空间和颗粒结构，与铸体薄片的形态几乎完全相同（附图18），佐证了三维重构的准确性。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种火山岩CT图像的预处理和确定分割阈值的方法，其特征在于，该火山岩CT图像的预处理和确定分割阈值的方法包括以下步骤：

步骤一，综合考虑反映火山岩岩石的特征，选择具有代表性的岩样，从喷溢相、爆发相、浸入相、火山通道相和火山沉积相五种岩相中，选出与油气分布有关系的岩相/岩性，做CT扫描实验，工业CT、微米CT或纳米CT；

步骤二，火山岩2D CT图像的采集，按完全投影数据方式进行图像序列采集，获得CT三维重构的一手2D CT图像资料；

步骤三，2D CT图像数据的预处理，将不清晰、信息量不丰富的原始图像变成清晰、信息量丰富的可用图像，有效地去掉图像中的噪声、突出图像中边缘或感兴趣的区域，引进邻域内像素的方差作为均值滤波和中值滤波两种方法的衔接纽带，既能消除噪声的影响又能增强图像的边缘；

步骤四，叠加图像形成三维数据体；以预处理过的2D CT图像为支撑点，利用Amira软件建立三维数据体，为火山岩CT三维重构奠定基础；

步骤五，三维重构体二值化阈值方法的改进及确立；综合学习了解阈值分割方法，进行改进提升，重新修订Otsu法，建立火山岩CT图像三维重构的二值化阈值分割方法；阈值分割法基于Otsu法进行改进提升，方法为：给定一幅火山岩CT图像灰度级的取值范围为G＝{0，1，2┄，L-1}，大小为M×N的图像，假设灰度级为i的像素点数为m_i，各灰度级的分布概率为：

$P_i=\frac{m_i}{\underset{i=0}{\overset{L-1}{\Σ}}{m}_i},i=0,\mathrm{......},L-1$

设定一个阈值T，将图像灰度级分为C₀和C₁两组，u₀和ω₀为C₀组的灰度平均值和产生的概率；u₁和ω₁为C₁组的灰度平均值和产生的概率，两组间的方差表达式如下式：

$\mathrm{\δ}\left(t\right)=\sqrt{{\mathrm{\ω}}_0*{(u_2-u)}^2+{\mathrm{\ω}}_1*{(u_1-u)}^2}$

其中u是整幅图像的灰度平均值，如下式：

${\mathrm{\ω}}_0{u}_0+{\mathrm{\ω}}_1{u}_1={\mathrm{\ω}}_0\underset{i=0}{\overset{T}{\Σ}}\frac{{\mathrm{iP}}_i}{{\mathrm{\ω}}_0}+{\mathrm{\ω}}_1\underset{i=T+1}{\overset{L-1}{\Σ}}\frac{{\mathrm{iP}}_i}{{\mathrm{\ω}}_1}=\underset{i=0}{\overset{L-1}{\Σ}}{\mathrm{iP}}_i=u$

进而将两组间的方差和每一组内的方差都考虑进来，进行改进提升，建立适合火山岩CT图像的阈值分割方法，得到公式，如下：

$Z=\frac{\mathrm{\δ}\left(t\right)}{{\mathrm{\δ}}_{u_0}\left(t\right)+{\mathrm{\δ}}_{u_1}\left(t\right)}=\frac{\sqrt{{\mathrm{\ω}}_0*{(u_2-u)}^2+{\mathrm{\ω}}_1*{(u_1-u)}^2}}{\sqrt{\underset{0}{\overset{T}{\Σ}}{({\mathrm{\ω}}_0-P_i)}^2}+\sqrt{\underset{T+1}{\overset{L}{\Σ}}{({\mathrm{\ω}}_1-P_i)}^2}}$

当Z值取最大时，此时，两组间的相似性最小，方差值最大，而其中的每一组组内方差值最小，找到最优化的阈值T，此时就可取得最优阈值T：

T＝argmaxδ²(t)

基于上述确定的阈值，利用Amira软件进行火山岩CT图像三维重构体的二值化；

步骤六，建立火山岩CT图像的三维重构模型，并进行三维显示；基于步骤三获得预处理后的火山岩2D CT图像，进一步结合步骤四、步骤五，实现火山岩CT图像二值化的三维重构，再利用Amira软件进行体显示。

2.如权利要求1所述的火山岩CT图像的预处理和确定分割阈值的方法，其特征在于，在步骤三中，CT图像预处理是对CT图像进一步处理与分析的先行步骤，设M是3×3的正方形邻域，当正方形邻域内的方差大于等于ε_临界值时，用中值滤波进行图像预处理，当正方形邻域内的方差小于ε_临界值时，用均值滤波进行预处理，ε_临界值＝18。

3.如权利要求1所述的火山岩CT图像的预处理和确定分割阈值的方法，其特征在于，均值滤波的预处理方法为：对于给定的图像f(x,y)中的每一个点(m,n)，取邻域A，设A含有N个像素，取平均值作为预处理后所得图像像素点(m,n)处的灰度，设M是3×3的正方形邻域，点(m,n)位于M中心，具体公式如下：

中值滤波的预处理方法为：给定图像f(x,y)中的每一个点(m,n)的邻域A，设A含有N个像素{a1,a2,┄,aN}，将按像素大小排序，若N是奇数时，则位于中间的那个像素值就是修改后图像g(x,y)在(m,n)处的像素值，若N是偶数时，则取中间两个像素的平均值就是修改后图像g(x,y)在(m,n)处的像素值，具体表达式如下：

4.如权利要求1所述的火山岩CT图像的预处理和确定分割阈值的方法，其特征在于，

建立火山岩CT图像预处理的方法，具体公式如下

建立度量灰度差异Z值的相关函数，当Z值最大，确定火山岩CT三维重构体二值化的最优阈值。