CN102879823B - 一种基于快速独立分量分析的地震属性融合方法 - Google Patents

Abstract

本发明涉及独立分量分析（ICA）技术领域和地震多属性融合领域，提供了一种基于快速独立分量分析（FICA）的地震多属性融合方法。本发明方案是通过把参与融合的属性都分成相同大小和块数的属性块，从这些块里面选出一定数量的块，应用FICA的原理对选出的块做处理，求得分离矩阵和与之互逆的混合矩阵，用分离矩阵把所有的块映射到ICA域。根据融合规则分别对各属性对应的块在ICA域做融合，最后用把ICA域的融合结果映射到空间域得到融合结果。融合结果有助于分析复杂的地层信息，提高储层预测的精度。本方法可广泛应用于地震属性分析、综合解释、地震储层预测、岩性与流体识别。

Description

一种基于快速独立分量分析的地震属性融合方法

技术领域

本发明涉及独立分量分析技术领域和地震多属性融合领域，是一种采用快速独立分量分析(FICA)技术进行地震属性融合的方法，本方法可广泛应用于地震属性分析、综合解释、地震储层预测、岩性与流体识别。

背景技术

近年来，地震属性技术发展迅速，已广泛应用于地层分析、油藏特征描述以及油藏动态检测等各个领域，成为了油藏地球物理的核心部分。地震属性也从早期的振幅属性发展到目前常用的数百种，然而由于单一的地震属性往往受观测条件、测量精度、地质的复杂性、主观因素、采集和处理过程中引起的多解性等因素的影响，缺乏可信度。因此，地震多属性融合技术研究的开展有利于降低多解性，为降低油气勘探开发风险做出科学合理的指导。而且随着属性分析理论及技术本身的发展，多属性综合或融合研究已经在国内外地震属性分析及储层预测应用中兴起。

属性融合方法是基于单一属性在储层预测中多解性的基础上提出来的。地震多属性融合就是利用计算机自动实现两种或两种以上属性整合的有效手段，也是地震属性分析、综合解释必备手段，是一项高级地震属性分析和综合解释技术。多属性融合可以说是一种多属性信息融合。信息融合的目标是通过对信息的优化组合导出更多的有效信息，它的最终目的是利用多信息共同或联合的优势来提高整个系统的有效性。较早而且也是比较简单的地震多属性融合方法是加权融合法。加权融合算法由于其实现起来特别容易,所以它在实际工程上得到了广泛的应用。但由于其加权因子直接影响融合结果,导致算法性能很不稳定，一般说来都根据平时经验制定其加权因子。该方法的基本原理是给要融合的n个地震属性分别赋予一个权值a_i(0<a_i<1，i＝1,…,n),这些权值还满足然后把与对应权值相乘之后的所有属性加和，得到一个融合的属性。加权融合虽然能得到融合的属性，但是由于融合权重需要人为干预，即需要预先确定优势属性，会使融合结果比较粗糙，精度不高，从而影响储层预测工作的开展。

地震多属性融合的方法还有很多，除了上面介绍的最简单的加权融合外，还有多属性线性回归法、主成分分析(PCA)法、人工神经网络(ANN)法、小波多分辨分析法和RGBA彩色融合法等属性融合的方法。尽管上述融合方法中的部分方法应用已经较为成熟，但是某些方法在应用中依然存在一些缺陷，如提及到的加权融合。地震多属性融合技术作为目前国内外许多石油公司及科研机构的研究重点之一，其完善的理论体系以及有效的广义融合模型和算法是必不可少的。所以，不断地完善融合理论体系、拓展融合方法，建立有效的融合方法模型和研究相应的算法是促进多属性融合技术发展的关键。当然，伴随着多属性融合方法的研究，建立多属性的融合规则和融合结果的评价标准也是重点，只有这样才能有助于在地震储层预测、岩性与流体识别方面取得良好的效果。

发明内容

本发明的内容是针对技术背景中地震属性融合各方法存在的缺陷及作为一种丰富地震属性融合的方法而提出来的。研究设计了一种基于快速独立分量分析(FICA)的地震多属性融合方法，通过应用快速独立分量分析理论对多个属性进行融合，以达到在快速得到融合结果的同时，有效提高所融合属性的准确度，进而有效地提高地质储层预测精度的目的。

本发明为了实现上述目的采用以下技术方案：

一种基于快速独立分量分析的地震属性融合方法，其特征在于包括：

(1)对待融合属性数据分块：选取大小为N×N的窗W对各大小均为M1×M2待融合属性切片平滑分块，把分出来的每一个小块存储成N²×1的向量，放到各属性块的集合{Ii}里，得到一个块的集合C，C＝{{I1},{I2},…,{In}}。分块的原则是从属性切片的第一个值开始，按从左往右，从上到下的顺序平滑分块。C这个集合的大小为行是M1-N+1，列是(M2-N+1)*n，它的每一个元素是N²×1的向量。

其中

(2)选块：根据集合C的维数和各小块的均值，从集合C里面随机的选取P块(近似1000块)属性块作为ICA算法的输入来求取ICA变换的变换核T{·}。选块的过程描述如下(在Matlab环境中)：求出集合C里面每个元素所代表的属性块的均值，把所有均值存储到矩阵c里面(c的行和列与集合C相同)，把矩阵c里面的元素以列为单位，组成一个按列接连起来的向量c1，对向量c1的元素按从小到大排序，记录下排序后元素的索引值组成的向量l。对k0＝[(M1-N+1)*(M2-N+1)*n]/1000四舍五入，以k0/n+k0*(i-1)<＝(M1-N+1)*(M2-N+1)*n作为判决条件(i的初始值取1)，取m＝l(k/n+k*(i-1))，然后取p＝m％(M1-N+1)来确定选取的块在集合C中所对应的行，q＝(m-p)/M+1来确定选取的块在集合C中随对应的列，然后选中集合C的第p行、第q列的元素所对应的属性块放到矩阵mixedX中(作为mixedX的一列)，i这个变量自加，接着再回到判决条件，以此类推，直到不满足判决条件为止，选块完成。

(3)求取ICA域的变换核：由步骤2得到一个由P块属性块组成的矩阵mixedX,其维数为N²×P，把mixedX作为快速独立分量分析算法的输入，求出mixedX的协方差矩阵CovX，CovX的维数为N²×N²，经过特征分解得到白化矩阵WhiteMatrix和去白化矩阵deWhiteMatrix，去白化矩阵与白化矩阵是互逆矩阵，在融合部分都要用到，通过求出的白化矩阵，求出白化信号Whitesig，把求得的白化信号作为快速独立分量核心算法中的x，求出ICA域的变换核T{·}即分离矩阵W和ICA域的变换核的逆T^-1{·}即混合矩阵A，去白化矩阵deWhiteMatrix的表达式如下：

$\mathrm{deWhiteMatrix}={\mathrm{ED}}^{\frac{1}{2}}$

得到白化信号Whitesig＝WhiteMatrix×mixedX。把求得的白化信号作为快速独立分量核心算法中的x，求出ICA域的变换核T{·}即分离矩阵W。ICA域的变换核 $T\{\&CenterDot;\}=W={[w_1,w_2,...,w_p,...,w_{N^2}]}^T\&times;\mathrm{WhiteMatrix},$ 其维数为N²×N²。ICA域的变换核的逆 $T^{-1}\{\&CenterDot;\}=A=\mathrm{deWhiteMatrix}\&times;[w_1,w_2,...,w_p,...,w_{N^2}],$ 其维数为N²×N².

(4)空间域映射到ICA域：把集合C＝{{I1},{I2},…,{In}}中{I1},{I2},…,{In}对应的块I1i,I2i,…Ini(i＝1,2,…,m,m为一个属性所分的块数)去均值(各块的属性值均减去各块的均值)，得到I₀1i,I₀2i,…I₀ni，均值为m1i,m2i,…mni，然对I₀1i,I₀m2i,…I₀ni分别做如下变换，SS1i＝W×I₀1i，SS2i＝W×I₀2i，…，SSni＝W×I₀ni，得到各属性第i块的独立分量。依照上述的方法，求出各属性对应块的独立分量，即把空域的块映射到ICA域，SS1i,SS2i,…,SSni均为N²×1。

(5)ICA域融合：由步骤4把各属性的所有块均映射到ICA域里面，用如下的融合规则得到第i块在ICA域的融合结果：

$\mathrm{Fi}\left(l\right)=\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\left|\mathrm{SSki}\left(l\right)\right|\mathrm{SSki}\left(l\right)}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}\left|\mathrm{SSki}\left(l\right)\right|},l=1,...,N^2---\left(2\right)$

依照公式(2)求出所有块在ICA域的融合结果；

(6)ICA域映射到空域：第i块在空域融合的结果,

$\mathrm{MFi}=A\&times;\mathrm{Fi}+\frac{1}{n}\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{msi}---\left(3\right)$

依照公式(3)得到所有块在空域的融合结果，然后把(M1-N+1)*(M2-N+1)*n块属性块的空间域融合结果(融合的属性块有(M1-N+1)*(M2-N+1)块)组合成M1*M2大小的融合属性。对得到的融合结果采用求平均值法分块去重叠处理，得到最终的融合结果。

本发明具有以下有益效果：

本发明利用快速独立分量分析对多个2维属性做融合，与加权融合或某些融合方法比起来，快速独立分量分析保持了各属性特征的高阶相互独立性，其融合过程最大可能地保证了把各属性的主要信息体现在其融合结果上，且矩形窗(即分的块大小)越小，数据处理量少，具有简捷、快速可靠，效用准确度高，处理能力强等特点，可有效提高地震储层预测的精度。

附图说明

图1FICA原理结构流程图；

图2为本发明的流程示意图(方框图)；

图3为具体实施方式的融合属性，图中3a为均方根振幅地震属性，3b为相干体切片地震属性，3c为FICA融合的属性，3d为PCA融合属性，3e为加权融合属性。

具体实施方式

本发明用到的基本原理是快速独立分量分析，它是独立分量分析(ICA)的一种快速算法。独立分量分析技术源于盲信号分离问题，它与主成分分析(PCA)，奇异值分解同属于线性变换技术。由于后两者按能量大小对数据进行分解，因此只能消除数据之间的二阶相关性。而在地震属性融合原理及应用中，通常需要提取多组特征，特征之间的相关性隐藏在高阶统计特性中。使用独立分量分析的方法能够有效约减特征维数，保持特征的高阶相互独立性，比仅消除二阶相关性的主分量分析和奇异值分解方法更为有效。然而在运用中，独立分量分析的时间复杂度较高，算法效率低，基于这一不足，快速独立分量分析算法被提出。快速独立分量分析算法是芬兰赫尔辛基工业大学计算机及信息科学实验室Hyvarinen等人提出并发展起来的。它是一种基于非高斯性最大化原理，使用固定点迭代理论寻找w^Tx(FICA算法理论因子)的非高斯性最大值，该算法采用牛顿迭代算法对观测变量x的大量采样点进行批处理，每次从观测信号中分离出一个独立分量。快速独立分量分析算法继承了独立分量分析的优点，且具有较高的算法效率。综合某些融合方法没有的优点，我们首次提出了基于快速独立分量分析的地震属性融合技术。快速独立分量分析原理的流程为：

A.输入数据：输入同一个量级的地震属性数据，这些数据组成一个行为属性个数，列为各属性的对应的属性数据的矩阵X；

B.去均值：对X的各行的数据均减去对应行数据的均值，得到一个去均值矩阵X₀；

C.白化(球化)处理：对去均值的矩阵X₀做白化处理，这样做的目的是确保各属性间是互不相关的。具体过程如下：

(1)求矩阵X₀的协方差矩阵Cov，其表达式为

CovX＝(X₀×X₀ ^T)/MN

式中MN为一个地震属性所含有的属性数据数。

(2)对矩阵X₀的协方差矩阵Cov做特征值分解，得到一个特征值从小到大的特征值矩阵D和特征向量矩阵E，其中E的各列均与矩阵D的特征值相对应，且E的各列是归一化正交的。

(3)求白化(球化)矩阵Q，其表达式如下：

$Q=D^{\frac{1}{2}}{E}^T$

(4)求得白化信号Whitesig，表达式如下：

Whitesig＝Q×X₀

白化信号是快速独立分量分析实现的关键，它确保了输入信号之间是互不相关的。

D.FICA算法的核心，求取分离矩阵W；

FICA算法的步骤如下：

(1)选择一个初始的权矢量w；

(2)令w⁺＝E[x^Tg(x^Tw)]-E[g′(x^Tw)]w；

(3)令w＝w⁺/||w⁺||；

(4)如不收敛，则返回(2)。

上面的算法步骤用到的x表示的是步骤C中求得的白化信号Whitesig。该过程求取的是一个独立分量，对于求多个独立分量的算法步骤是求一个独立分量的重复，具体如下：

(1)得到前p个向量w₁，w₂，…，w_p后，再次调用上面的算法得到w_p+1，

(2)令： $w_{p+1}=w_{p+1}-\underset{j=1}{\overset{p}{\mathrm{\&Sigma;}}}{w^T}_{p+1}{w}_j{w}_j;$

(3)令： $w_{p+1}=w_{p+1}/\sqrt{{w_{p+1}}^T{w}_{p+1}};$

上述过程求出的w₁，w₂，…，w_p，…，w_n均为归一化正交的。FICA算法步骤中，w⁺＝E[x^Tg(x^Tw)]-E[g′(x^Tw)]w中的函数g(.)是一些非二次型函数G(.)的一阶导数，g’(.)是G(.)的二阶导数。较好地选择G，可以得到稳健的估计器，常用的非二次型函数有：

$G_1\left(u\right)=\frac{1}{a_1}\log \cos \left(a_{1}u\right)$

$G_2\left(u\right)=-\frac{1}{a_2}\exp (-a_2{u}^2/2)$

$G_3\left(u\right)=-\frac{1}{4}{u}^4$

其中，1≤a₁，a₂≤2是适当的常量。

本发明的具体解决方案是通过把每一个参与融合的属性分成相同大小和块数的属性块，并从这些属性块里面随机选出一定数量的块，然后应用快速独立分量分析的原理对选出的属性块做处理，求得一个能把所有属性块变成相互独立状态的分离矩阵和一个与分离矩阵互为逆矩阵的混合矩阵。再运用求得的分离矩阵把所有的属性块映射到独立分量分析域，用相应的融合规则分别对各属性相应位置的属性块在ICA域做融合，最后用与分离矩阵互为逆矩阵的混合矩阵把独立分量分析域的融合结果映射到空间域，得到融合结果。本融合方法运用独立分析原理中的快速独立分量分析理论对选出的块做特征分析，不仅在效率上有所提高，还能突显出各属性的特征，使得得到融合结果能很好的突显主要信息。通过融合结果，能有效的分析复杂的地层信息，提高了储层预测的精度。本发明即以此实现其发明目的。

本实施方式以均方根振幅地震属性、相干体切片地震属性的融合为例，其实施方式的具体步骤如下：

步骤A.分别对大小均为201*301的均方根振幅属性、相干体切片属性平滑分块，矩形窗的大小取5*5，并把分出来的每一个小块存储成5²×1(即25*1)的向量，得到一个块的集合C，C这个集合的大小为行是197，列是297*2；

步骤B.从集合C中选出近似1000块的属性块作为FICA算法的输入，根据选块的原则，选出1000块属性块，得到矩阵mixedX的行是25，列是1000。

步骤C.把矩阵mixedX作为FICA算法的输入，求先求矩阵mixedX的协方差矩阵Cov，得到一个25*25的协方差矩阵Cov，之后对这一协方差矩阵Cov做特征值分解，得到一个25*25的特征值从小到大的特征值矩阵D和25*25的特征向量矩阵(对角矩阵)E，根据求出的矩阵D和矩阵E求出白化矩阵WhiteMatrix和去白化矩阵deWhiteMatrix，二者的都是25*25大小的矩阵。通过求出的白化矩阵，求出白化信号Whitesig，把白化信号Whitesig代替FICA核心算法中的x，求出ICA域的变换核T{·}即分离矩阵W和ICA域的变换核的逆T^-1{·}即混合矩阵A，分离矩阵W和混合矩阵A的维数是25*25。

步骤D.从均方根振幅属性块{I1}和相干体切片属性块{I2}中选出对应的第一小块I11,I21，对它们去均值得到I₀11,I₀21，均值为m11,mm21，然对I₀11,I₀m21分别做如下变换，SS11＝W×I₀11，SS21＝W×I₀21，分别得到均方根振幅属性块和相干体切片属性的第1块的独立分量。依照上述的方法，求这两个属性对应块的独立分量，即把空域的块映射到ICA域里面，SS11,SS21均为5²×1。

步骤E.应用ICA域的融合公式：

$F1\left(l\right)=\underset{k=1}{\overset{2}{\mathrm{\&Sigma;}}}\frac{\left|\mathrm{SSk}1\left(l\right)\right|\mathrm{SSk}1\left(l\right)}{\underset{k=1}{\overset{2}{\mathrm{\&Sigma;}}}\left|\mathrm{SSk}1\left(l\right)\right|},l=1,2,...,5^2$

求得均方根振幅属性块和相干体切片属性第1块的融合结果，从而求出均方根振幅属性块和相干体切片属性的其他对应属性块的ICA域的融合结果。

步骤F.根据步骤E求出的ICA域的融合结果，均方根振幅属性块和相干体切片属性第1块在空域融合的结果为依次求出均方根振幅属性块和相干体切片属性的其他对应属性块的空域融合的结果。把197*297*2块属性块的空间域融合结果组合成201*301大小的融合属性。对得到的融合结果采用求平均值法分块去重叠处理，得到最终的融合属性。

Claims (5)

1.一种基于快速独立分量分析的地震属性融合方法，其特征在于包括：

(1)对待融合属性数据分块：选取大小为N×N的窗W对各大小均为M1×M2待融合属性切片平滑分块，把分出来的每一个小块存储成N²×1的向量，放到各属性块的集合{Ii}里，得到一个块的集合C，C＝{{I1},{I2},…,{In}}；

(2)选块：根据集合C的维数和各小块的均值，从集合C里面随机的选取P块属性块作为ICA算法的输入来求取ICA变换的变换核T{·}；

(3)求取ICA域的变换核：由步骤(2)得到一个由P块属性块组成的矩阵mixedX,其维数为N²×P，把mixedX作为快速独立分量分析算法的输入，先对mixedX去均值，然后求出mixedX的协方差矩阵CovX，CovX的维数为N²×N²，经过特征分解得到白化矩阵WhiteMatrix和去白化矩阵deWhiteMatrix，去白化矩阵与白化矩阵是互逆矩阵，通过求出的白化矩阵，求出白化信号Whitesig，把求得的白化信号作为快速独立分量核心算法中的x，求出ICA域的变换核T{·}即分离矩阵W和ICA域的变换核的逆T^-1{·}即混合矩阵A；

(4)空间域映射到ICA域：求出各属性对应块的独立分量，即把空间域的块映射到ICA域；

(5)ICA域融合：由步骤(4)把各属性的所有块均映射到ICA域里面，用如下的融合规则得到第i块在ICA域的融合结果：

$\mathrm{Fi}\left(l\right)=\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\left|\mathrm{SSki}\left(l\right)\right|\mathrm{SSki}\left(l\right)}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}\left|\mathrm{SSki}\left(l\right)\right|},l=\mathrm{1,2},...,N^2---\left(2\right)$

依照公式(2)求出所有块在ICA域的融合结果；

(6)ICA域映射到空间域：第i块在空间域融合的结果,

$\mathrm{MFi}=A\&times;\mathrm{Fi}+\frac{1}{n}\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{msi}---\left(3\right)$

依照公式(3)得到所有块在空间域的融合结果，然后把(M1-N+1)*(M2-N+1)*n块属性块的空间域融合结果组合成M1*M2大小的融合属性，对得到的融合结果采用求平均值法分块去重叠处理，得到最终的融合结果。

2.根据权利要求1所述的一种基于快速独立分量分析的地震属性融合方法，其特征在于：所述步骤(1)中分块的原则是从属性切片的第一个值开始，按从左往右，从上到下的顺序平滑分块，C这个集合的大小为行是M1-N+1，列是(M2-N+1)*n，它的每一个元素是N²×1的向量；其中窗W

3.根据权利要求1所述的一种基于快速独立分量分析的地震属性融合方法，其特征在于：对步骤(3)中协方差矩阵CovX做特征值分解，得到一个特征值从小到大的特征值矩阵D和特征向量矩阵E，根据求出的矩阵D和矩阵E求出白化矩阵WhiteMatrix和去白化矩阵deWhiteMatrix，去白化矩阵deWhiteMatrix的表达式如下：

$\mathrm{deWhiteMatrix}={\mathrm{ED}}^{\frac{1}{2}}$

式中D为特征值从小到大的特征值矩阵，E为特征向量矩阵，得到白化信号Whitesig＝WhiteMatrix×mixedX，把求得的白化信号作为快速独立分量核心算法中的x，求出ICA域的变换核T{·}即分离矩阵W和ICA域的变换核的逆T^-1{·}即混合矩阵A。

4.根据权利要求1所述的一种基于快速独立分量分析的地震属性融合方法，其特征在于：所述ICA域的变换核T{·}＝W＝[w₁,w₂,…,w_p,…,w_N2]^T×WhiteMatri x，其维数为N²×N²，ICA域的变换核的逆 $T^{-1}\{\&CenterDot;\}=A=\mathrm{deWhiteMat\; rix}\&times;[w_1,w_2,...,w_p,...,w_{N^2}],$ 其维数为N²×N²，w₁,w₂,…,w_p,…的含义是快速独立分量分析算法经过迭代求出的矩阵W的第p个列向量。

5.根据权利要求4所述的一种基于快速独立分量分析的地震属性融合方法，其特征在于：所述步骤(4)中求出各属性对应块的独立分量方法为，把集合C＝{{I1},{I2},…,{In}}中{I1},{I2},…,{In}对应的块I1i,I2i,…Ini，i＝1,2,…,m,m为一个属性所分的块数，去均值，各块的属性值均减去各块的均值，得到I₀1i,I₀2i,…I₀ni，均值为m1i,m2i,…mni，然后对I₀1i,I₀2i,…I₀ni分别做如下变换，SS1i＝W×I₀1i，SS2i＝W×I₀2i，…，SSni＝W×I₀ni，得到各属性第i块的独立分量，SS1i,SS2i,…,SSni均为N²×1。