一种基于多属性超体素图割的河道三维建模方法与流程
本发明属于地质体三维模型构建领域，特别涉及一种基于多属性超体素图割的河道三维建模方法。
背景技术：
地质体(河道)三维模型构建是地震数据解释中最重要的任务之一，因为大多数重要的油气藏都存在于地质体周围。地质体的三维模型不仅能够直观地表现地质体的构造形态以及在三维空间的分布状况，使解释人员能够对地质体进行定量的分析，同时也为油藏数值模拟、储量计算井位部署提供重要依据。
边缘检测的方法是用于地质体检测的一种常见方法，采用基于边缘检测的方法识别盐丘边界，但由于地震数据中存在较大的噪声，这种方法并不能达到理想的效果。现有技术中还有一种将3d边缘检测器与倾角导向结合的方法用于地质体的检测，提高了信噪比与边缘连续性，使地质体的边界形状更清晰。但基于边缘的方法对振幅的变化依赖性较强，在边缘的瞬时振幅变化不明显时，无法得到很好的检测效果。为了克服这些问题，引入了基于纹理属性的地质体识别方法。可使用一组纹理属性来预测3d立方体中每个像素属于盐丘的概率，再通过分割找到盐丘边界。现有技术针对边缘处像素变化较大这一特点，使用纹理梯度来检测盐丘的边界。但对于基于纹理的方法，窗口的大小对识别结果有很大的影响，对不同工区来说不具有普适性。有一种组合的边缘与纹理属性的混合分类方法，兼具边缘法和纹理法的优点，能较好的检测盐丘。进一步将边缘检测、几何和纹理三种属性结合起来，再通过支持向量机(svm)进行半自动的断层检测，检测结果与原始图像较匹配。
但与盐丘、溶洞等地质体不同，河道本身并没有明显、确定的边界，同时河道经常变化、交叉，连续性较差，因此以上地质体识别分割的方法并不能很好的应用在河道分割中。
很多学者从地震属性的角度来分割河道。地震属性是通过一定的算法从原始三维地震数据中提取出的数据体，地震属性从不同的角度反映了地震数据的特点，不同的地质构造在属性值上显示出不同的特征。首先用地震属性来辅助分析河道，通过对多个地震属性切片图的对比分析，综合比较划分河道边界。但在此方法中，数据的信噪比极大地影响了地震属性分辨率与质量，且部分地震属性对河道的划分并没有帮助。将rgb多属性融合技术应用到河道的识别中，使河道边缘更加清晰，分辨率提高。通过相干性和谱分解将地震属性分解为三个频段，分别用不同颜色表示，再通过rgb混频技术进行可视化，能够快速有效的识别河道。但rgb融合技术只能融合三种属性，具有一定的局限性，在需要融合更多种属性时不适用，并且大多数研究只是使用可视化技术辅助主观判断，并没有定性的分割出河道边界。
近年来，很多学者把地质问题与图像结合起来以达到准确分割目的。将改进的prc分割算法应用于二维和三维的地震图像中，能够半自动的精确分割盐丘边界。使用最优路径拾取算法在新图像中提取盐丘边界，该算法通过选择具有全局最大包络值的最佳路径来跟踪高度不连续的盐丘边界，能够快速更新得到盐边界。
基于图像的方法在河道地质体的问题上也同样适用。为了解决河道结构复杂、连续性差等问题，通过图像处理中steerablepyramid的方法来增强河道的局部线性特征。有一种置信度和曲率引导的水平集方法从三维地震数据中分割河道。在融合多种地震属性的基础上，用水平集的方法对溶洞、河道等多种地质体进行了分割建模。但基于水平集的方法比较依赖初始形状，且计算效率比较低。
综合来看，现有的对河道地质体的研究大多停留在二维平面的识别检测上，对河道三维模型构建的研究有些不足。
技术实现要素：
由于单属性数据所携带的地质信息并不完整，所以仅仅在单属性数据上对河道进行解释往往不够准确。因此需要合理的选取几种可以信息互补的属性，并采用非线性降维的方法，在降维的同时，不至于损失数据信息。在进行河道分割时，如果仅仅是基于像素级的分割，分割出的河道边缘连续性较差，不能很好的解释河道。为了解决现有技术中的问题，本发明提供了一种基于多属性超体素图割的河道三维建模方法。
一种基于多属性超体素图割的河道三维建模方法，包括以下步骤：
步骤1，采用isolle的非线性降维算法融合河道地震属性，得到河道的属性数据体；
步骤2，根据所述数据体，通过简单线性迭代聚类算法(slic)，生成地质超体素；
步骤3，通过k-means聚类方法，建立目标区域和非目标区域的高斯混合模型，构建网络图和能量函数，基于最小割准则对地质数据进行分割，得到二值化的分割结果，并通过提取等值面的方法得到河道地质体的三维模型。
进一步地，所述步骤1包括以下流程：
步骤11，根据测地距离搜寻与地震数据样本点近邻的k个样本；
根据两点之间的测地距离，在三维数据中搜寻与每个地震数据样本点i测地距离最相近的k个数据点
dg(xi,xj)＝min{lg(xi,xj)}
其中，lg为两点之间某路径的长度，de为欧氏距离，dg为两点之间的测地距离；
步骤12，构造局部最优化重建权值矩阵；
引入误差函数以衡量重构误差大小，其为
其中，xij(j＝1,2,...,k)为地震数据点i的k个近邻点，wij为xi和xij之间的权值，wij符合
对于每个地震数据点，误差为
构造局部协方差矩阵
结合所述局部协方差矩阵和
当所述局部最优化重建权值矩阵为奇异矩阵时，进行正则化处理
qi＝qi+r·i
其中，r为正则化参数，i为k×k单位矩阵；
步骤13，将所有地震数据点从高维向低维空间进行映射，通过该地震数据点的局部重建权值矩阵以及它的k个近邻点计算出该地震数据点在低维空间的值；
映射条件满足
其中，ε(y)为损失函数值，yi为xi输出向量，yij(j＝1,2,…,k)为通过测地距离寻找到的的k个近邻点，同时还满足
其中，i为m×m的单位矩阵；损失函数为
其中，m为n×n的对阵矩阵，表示为
m＝(1-w)t·(i-w)。
进一步地，所述步骤2包括以下流程：
步骤21，初始化聚类中心；
设置地质超体素的初始个数为k，将河道工区中地震数据点i的标签初始化设为-1，即labeil＝-1，点i与聚类中心j的距离初始化设为无穷大，即distij＝+∞；
步骤22，在聚类中心cj的邻域内，计算各点到cj的距离；
距离为
其中，
当dij＜disti，更新disti＝dij，labeli＝j；
步骤23，更新聚类中心
步骤24，计算残差
步骤25，更新聚类中心，令cj＝c'j，若e小于预设阈值或超过最大迭代次数，进入下一流程；若e不满足条件，则流程回到所述步骤22；
步骤26，对于分割出的地质超体素，通过在三维空间中遍历地震数据点的连通区域，建立邻接矩阵a,将属性直方图作为地质超体素特征；
对体积不大于阈值的地质超体素i，计算地质超体素i与相邻超体素j之间的巴氏距离
其中，m为地质超体素灰属性直方图的维度，将地质超体素i合并到距离t最小的相邻地质超体素中；
步骤27，将所有离散的、小于阈值的地质超体素合并完成后，更新邻接矩阵及地质超体素的属性直方图。
进一步地，所述步骤3包括以下流程：
步骤31，通过k-means聚类方法，建立目标区域和非目标区域的高斯混合模型；
标记属于河道的超体素αi＝1，其余超体素αi＝0；
通过glcm四个属性作为特征对地质超体素进行k-means聚类；
根据聚类结果初始化高斯混合模型中每一个高斯分量的参数，权重ωk，样本均值uk和协方差∑k，高斯混合模型的密度函数为
其中，
计算后验概率
计算高斯分量参数的最大似然估计
通过迭代至似然函数收敛，得到河道高斯混合模型和非河道高斯混合模型；
步骤32，构建网格图和能量函数；
采用灰度共生矩阵来计算每个超地质体素的纹理属性。glcm是图像中两个灰度值的联合概率分布，选取glcm的熵、相异性、能量三种属性作为超体素的特征值，融合了地质超体素区域信息和边缘信息的能量函数为
其中，c为网格图的分割；er为区域项，代表网格图汇总t-link边的权值，反映了地质超体素的区域信息；eb为边界项，代表了网络图中n-link边的权值，体现出分割的边界属性；通过构建能量函数，将河道的分割转换成能量函数最小化问题；
步骤33，基于最小割准则对地质数据进行分割，并更新高斯混合模型中的地震数据点和数据点标记αi，迭代结束后输出分割结果，得到河道地质体的二值化分割结果，通过提取等值面的方法得到河道地质体的三维模型。
本发明的有益效果：本发明提供了一种基于多属性超体素图割的河道三维建模方法，通过融合几种优选地震属性得到新的数据体，再通过超体素图割的方法进行二值化分割，最后通过等值面提取得到河道表面。本发明中isolle的融合方法能有效地融合多种地震属性，合理的地震属性选取有助于信息互补，更全面地刻画河道地质体，且多属性融合的方式能够保持河道地震数据中的非线性关系，得到的新属性体为下一步的分割和重建打下了良好的基础；河道地质超体素具有良好的边缘保持性，地质超体素与图割相结合能准确、快速地分割河道，通过等值面提取实现了由三维地震数据河道地质体三维模型的转化，更加直观地展示了河道地质体的空间特征；本发明提出的基于多属性融合超体素图割的方法实现了对河道地质体的精准刻画，最终得到的三维模型符合地质规律，且与地质人员得到的结果基本一致，为后续工作奠定了基础。
附图说明
图1为本发明的流程图。
图2为图1中步骤1的流程图。
图3为图1中步骤2的流程图。
图4为初始化的聚类中心图。
图5为图1中步骤3的流程图。
具体实施方式
下面结合附图对本发明的实施例做进一步的说明。
请参阅图1，本发明提供了一种基于多属性超体素图割的河道三维建模方法，通过以下步骤实现：
步骤1，采用isolle的非线性降维算法融合河道地震属性，得到河道的属性数据体。
本实施例中，在地震解释的过程中不同的属性数据从不同角度刻画地质构造，通过融合不同属性，可以对储层的地质构造有更精确的解释。由于地震属性与地质特征的关系通常是非线性的，基于线性变换的pca方法不能充分地反映这种非线性关系，降低了预测识别的精度。非线性的lle降维算法采用的是欧式距离计算点与点之间的距离，并不能反映点之间的真实结构，且对近邻个数的选择十分敏感。
本实施例中，采用isolle算法，引入测地距离的概念，将欧式距离用测地距离替换，在保持了lle算法处理高维流形数据优势的同时，又提高了在相应低维空间中数据的紧密性和局部邻域数据的线性特征。本发明将isolle方法首次应用到三维河道地震数据的多属性融合中，选取几种属性融合成为新的属性体，扩展到三维空间。
本实施例中，选取了均方根振幅属性、能量属性、纹理同质性，频谱属性，瞬时频率属性来研究某地震工区的多属性融合问题。再经过归一化处理和傅立叶变换去噪后得到不同属性的切片图。同一属性在不同时间切片上的表现有好有坏，不同属性可以弥补彼此的不足。
请参阅图2，步骤1通过以下流程实现：
步骤11，根据测地距离搜寻与地震数据样本点近邻的k个样本；
根据两点之间的测地距离，在三维数据中搜寻与每个地震数据样本点i测地距离最相近的k个数据点
dg(xi,xj)＝min{lg(xi,xj)}
其中，lg为两点之间某路径的长度，de为欧氏距离，dg为两点之间的测地距离；
步骤12，得到地震数据点的局部重建矩阵。局部重建矩阵代表了多属性之间的局部线性关系，融合后的地震数据点之间也要能够保持这种线性关系。
引入误差函数以衡量重构误差大小，其为
其中，xij(j＝1,2,...,k)为地震数据点i的k个近邻点，wij为xi和xij之间的权值，wij符合
对于每个地震数据点，误差为
构造局部协方差矩阵
结合局部协方差矩阵和
当局部最优化重建权值矩阵为奇异矩阵时，进行正则化处理
qi＝qi+r·i
其中，r为正则化参数，i为k×k单位矩阵；
步骤13，将所有地震数据点从高维向低维空间进行映射，通过该地震数据点的局部重建权值矩阵以及它的k个近邻点计算出该地震数据点在低维空间的值；
映射条件满足
其中，ε(y)为损失函数值，yi为xi输出向量，yij(j＝1,2,…,k)为通过测地距离寻找到的的k个近邻点，同时还满足以下两个条件
其中，i为m×m的单位矩阵。这里的
其中，m为n×n的对阵矩阵，表示为
m＝(1-w)t·(i-w)。
矩阵m具有稀疏、半正定的特点。若要使得损失函数达到最小值，那么应该取为的最小个非零特征值所对应的特征向量。同时,把的特征值根据由小到大的顺序进行排列。第一个特征值约等于零,因此舍弃第一个特征值。
步骤2，根据所述数据体，通过简单线性迭代聚类算法(slic)，生成地质超体素。
本实施例中，基于像素级的分割，并未充分利用像素间的局部关系，在图像离散化过程中，会带来误差。超像素利用像素之间特征的相似性将像素分组,用少量的超像素代替大量的像素来表达图片特征，与像素级分割相比，它能保持原图像的边缘特征。超体素是超像素在三维空间中的扩展。在处理三维地震数据时，地质超体素的方法克服了在二维切片使用超像素最后合成三维模型时不平滑等问题，同时减少了分块个数，降低了计算量。
本实施例中，在地震数据体中，距离越近的两个地震数据点，属于同一个地质体的概率也越大。根据这一特点，我们将相似的地震数据点基于slic算法进行聚类处理生成地质超体素。生成的地质超体素既能保持内部均匀紧凑，又能保持地质体的边缘特性，降低了后续分割的计算复杂度。在slic算法中，主要步骤就是比较各点与聚类中心cj的距离进行分类，再更新聚类中心。假设融合后为crossline*inlinc*time大小的三维数据体，则所有地震数据点的集v＝{v1,v2,...,vn},其中vi＝{gi,xi,yi,zi}，gi为地震数据点i融合后的属性值，xi，yi，zi为网格点i在空间中的坐标。
请参阅图3，步骤2通过以下流程实现：
步骤21，初始化聚类中心；
设置地质超体素的初始个数为k，则超体素初始大小为
生成的种子点可能会落在河道的边缘，为了避免这种情况发生，计算出在当前种子点3×3×3邻域中所有网格点的梯度值，将种子点的位置移动到此邻域梯度最小处。即将聚类中心移动到g(x,y,z)＝ming(x,y,z)处。
以种子点为中心的3×3×3窗口内，计算河道内各地震数据点的属性梯度。计算方式如下：
g(x,y,z)＝||i(x+1,y,z)-i(x-1,y,z)||2
+||i(x,y+1,z)-i(x,y-1,z)||2
+||i(x,y,z+1)-i(x,y,z-1)||2
式中，i(x,y,z)表示在河道工区内坐标为(x,y,z)时的属性值。g(x,y,z)表示在河道工区内坐标为(x,y,z)时基于属性值计算出的梯度值。
将河道工区中地震数据点i的标签初始化设为-1，即labeli＝-1；
点i与聚类中心j的距离初始化设为无穷大，即distij＝+∞。
步骤22，在聚类中心cj的邻域内，计算各点到cj的距离；
对所有的聚类中心cj，在2s×2s×2s的邻域内，各点到聚类中心cj的距离为
其中，
当dij＜disti，更新disti＝dij，labeli＝j；
步骤23，更新聚类中心
步骤24，计算残差
步骤25，更新聚类中心，令cj＝c'j，若e小于预设阈值或超过最大迭代次数，进入下一流程；若e不满足条件，则流程回到所述步骤22；
步骤26，对于分割出的地质超体素，通过在三维空间中遍历地震数据点的连通区域，建立邻接矩阵a,将属性直方图作为地质超体素特征；
对体积不大于阈值的地质超体素i，计算地质超体素i与相邻超体素j之间的巴氏距离
其中，m为地质超体素灰属性直方图的维度，将地质超体素i合并到距离t最小的相邻地质超体素中；
步骤27，将所有离散的、小于阈值的地质超体素合并完成后，更新邻接矩阵及地质超体素的属性直方图。
本实施例中，考虑将相似的数据点进行聚类，看做河道的最小单位来进行分割。将融合后的河道数据体进行三维超体素分割，设定不同的地质超体素个数k，可以得到不同的结果。当设定地质超体素的数量k较小时，分割出的地质超体素在河道边缘不够贴合。当k较大时，河道内外的区分模糊，并且不能有效的降低图割的计算量。因此我们将超体素的数量k设置为150，分割出的超体地体素能够很好的贴合河道的边缘，并且地质超体素的同质性也更好。当然，k也可以设为其它数值。
步骤3，通过k-means聚类方法，建立目标区域和非目标区域的高斯混合模型，构建网络图和能量函数，基于最小割准则对地质数据进行分割，得到二值化的分割结果，并通过提取等值面的方法得到河道地质体的三维模型。
本实施例中，由于地质超体素是一个小的区域，区域存在着隐藏的纹理特征。因此在传统图割算法的基础上，通过融合超体素的属性值和纹理特征代替灰度属性来代表每个超体素的特性，首次将基于超体素的图割算法应用在河道地质体的分割上。与传统图割算法相比，本发明用高斯混合模型(gmm)替代属性直方图，来精确表达概率模型。高斯混合模型是多个高斯模型的线性叠加。gmm的核心在于求得每个高斯分量的均值、方差，以及权重比。首先对gmm中每一个高斯分量的参数进行随机初始化，计算其密度函数、后验概率以及高斯分量参数的最大似然估计。最后通过不断重复迭代，直至似然函数收敛，得到河道地质体的高斯混合模型，请参阅图5，具体通过以下流程实现：
步骤31，通过k-means聚类方法，建立目标区域和非目标区域的高斯混合模型；
标记属于河道的超体素αi＝1，其余超体素αi＝0；
设定k＝5，通过glcm四个属性作为特征对地质超体素进行k-means聚类；
根据聚类结果，随机初始化高斯混合模型中每一个高斯分量的参数，权重ωk，样本均值uk和协方差∑k，高斯混合模型的密度函数为
其中，
计算后验概率
计算高斯分量参数的最大似然估计
通过迭代至似然函数收敛，得到河道高斯混合模型和非河道高斯混合模型；
步骤32，构建网格图和能量函数；
采用灰度共生矩阵来计算每个超地质体素的纹理属性。glcm是图像中两个灰度值的联合概率分布，选取glcm的熵、相异性、能量三种属性作为超体素的特征值，融合了地质超体素区域信息和边缘信息的能量函数为
其中，c为网格图的分割；er为区域项，代表网格图汇总t-link边的权值，反映了地质超体素的区域信息；eb为边界项，代表了网络图中n-link边的权值，体现出分割的边界属性；通过构建能量函数，将河道的分割转换成能量函数最小化问题；
步骤33，基于最小割准则对地质数据进行分割，并更新高斯混合模型中的地震数据点和数据点标记αi，迭代结束后输出分割结果，得到河道地质体的二值化分割结果，通过提取等值面的方法得到河道地质体的三维模型。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。
- 还没有人留言评论。精彩留言会获得点赞！