CN102884448A - 在地球物理数据集中进行异常检测的窗口统计分析 - Google Patents

Abstract

使用窗口主分量分析(22)，或独立分量分析，或扩散映射分析(61)从地球物理或属性数据鉴定地质特征的方法。部分数据体或残差数据体中的细微特征可被鉴定。通过(36)消除没有被最突出的主分量(14)捕获的数据而产生残差数据体(24)。通过将数据(21)投影(35)在选择的主分量(22、61)上产生部分数据体。也可从用可变尺度数据相似性矩阵(73)产生的图形分析(77)或异常数据体(62,79)鉴定地质特征。该方法适于鉴定指示烃潜能的物理特征。

Description

在地球物理数据集中进行异常检测的窗口统计分析

相关申请的交叉引用

本申请要求于2010年5月6日提交的，题目为“在地球物理数据集中进行异常检测的窗口统计分析(Windowed Statistical Analysis forAnomaly Detection in Geophysical Datasets)”的权益。美国申请号12/775,226是2009年9月30日提交的PCT国际申请号PCT/US09/059044的部分继续，其要求于2008年11月14日提交的美国临时申请61/114,806和于2009年7月31日提交的美国临时申请61/230,478的权益，它们的公开内容以其整体通过引用并入本文。

发明领域

本发明主要地和一般地涉及地球物理勘探领域，更特别地涉及处理地球物理数据的方法。具体而言，本发明是不使用先前的训练数据而突出一个或多个地质学数据集或地球物理数据集如地震数据集中的区域的方法，所述区域代表包括潜在烃聚集(hydrocarbon accumulations)的真实世界的地质特征，并且其中期望的物理特征可以仅以细微形式出现在未处理的数据中，被更突出的异常所掩盖。

发明背景

地震数据集常常包含这样的复杂图形（pattern），其是细微的并且在多种地震数据体或属性/衍生(derivative)数据体(volume)中并以多重空间尺度出现。几十年来，地质学家和地球物理学家已经研发了一系列的技术以选取指示烃存在的许多重要的图形。但是，这些方法中的大部分涉及在一个数据体或者最多两个数据体中寻找具有预先指定特征的已知的或不严格限定的图形。这些“基于模板”或“基于模型”的方法常常错过不符合这类规定的细微的或意外的异常。在此并不进一步讨论这些方法，因为除了它们致力于相同的技术问题外，它们与本发明基本没有共同之处。

这些已知方法中的大部分涉及在一个数据体或者最多两个数据体中寻找具有预先指定特征的已知的或不严格限定的图形的人类解释者。这些“基于模板”或“基于模型”的方法常常错过不符合这类规定的细微的或意外的异常。因此期望开发能够在多重空间尺度的一个或多个地震数据体中自动突出异常区域的统计分析方法，而无需它们是什么和它们在何处的在先知识。本发明满足该需要。

发明内容

在一个实施方式中，本发明是用于在代表地下区域的一个或多个地球物理数据或数据属性的离散集(每个这样的数据集被称作“原始数据体”)中鉴定地质特征的方法，包括：(a)选择数据窗口的形状和大小；(b)对于每个原始数据体，将所述窗口移动至多个位置，并为每个窗口位置形成数据窗口向量，其分量由来自该窗口内的体素值构成；(c)进行数据窗口向量的统计分析，在多个原始数据体的情况下所述统计分析共同进行；(d)使用统计分析以鉴定数据中的离群值或异常；和(e)使用所述离群值或异常预测所述地下区域的地质特征。在本发明的一种实施方式中，统计分析技术是扩散映射（diffusion mapping），在该扩散映射中计算表示数据的基础。该基础是非线性变换的结果，其提供定义尺度概念的参数。

然后使用本发明方法鉴定的所述地质特征可以被用于预测烃聚集的存在。

附图简述

由于专利法限制，图1A-C、2和6是彩色原件的黑白复印件。彩色原件已经在对应的美国申请，美国申请号12/775,226中提交。在请求并支付必要的费用后，可从美国专利和商标局获得具有彩色附图的该申请的公开复印件。

通过参考以下详述和附图，本发明及其优势将更好地被理解，在附图中：

作为本发明方法的测试实例应用，图1A显示来自合成地震数据的3D数据体的图像(2D时间切片)；图1B显示通过本发明方法生成的、由前十六个主分量限定的原始图像的残差（residual），其占信息的90%；以及图1C以30×30窗口形式图解前十六个主分量；

图2是使用残差分析的本发明方法的一个实施方式中的基本步骤的示意图；

图3是示出使用单个窗口大小对多个数据体应用本发明的窗口PCA实施方式的基本步骤的流程图；

图4A-B示出数据体(大的矩形)和该数据在窗口中不同像素样本(较小的矩形)的2D切片图，图4A显示像素(1,1)的数据样本，图4B显示第i个像素的数据样本；和

图5A-B示出不在图4A-B的用于有效计算协方差阵的2D数据集的样本中的数据的细分。

图6是显示尺度参数的两个不同值的本发明扩散映射实施方式的示意图；和

图7是显示本发明扩散映射实施方式的基础步骤的流程图。

本发明将结合实例实施方式进行描述。其程度为以下描述具体到本发明的特定实施方式或特定应用，这些描述意欲仅是示例性的，并且不解释为是对本发明范围的限制。相反，意图覆盖所有的可被包括在本发明范围内的可选项、修改和等同物，如被所附权利要求限定的。

实例实施方式详述

本发明是不使用现有的培训数据，在多重空间尺度的多重地震或其它地球物理数据(例如，电磁数据)体中，检测异常图形的方法。本发明方法基于窗口统计分析，在本发明的一个实施方式中其涉及以下基本步骤：

1.提取(extracting)用户指定大小和形状的窗口内数据的统计分布。可使用标准的统计技术如主分量分析(PCA)、独立分量分析(Independent Component Analysis(ICA))、聚类分析。

2.通过(a)计算提取的分布中每个数据窗口的发生概率(或等价度量)(b)鉴定低概率数据区域作为可能的异常来提取数据中的异常区域。

提取统计分布不是本发明的必需步骤。直接或通过统计分布之外的其他技术，可从统计分析确认异常或离群值；例如，本发明的一些实施方式使用扩散映射。见图6和7。

本发明的特别便利的实施方式包括窗口主分量分析(“WPCA”)、残差分析和聚类分析的组合，其将在以下详细描述。但是，本技术领域的任何普通技术人员将容易理解其它统计分析技术可以如何被使用或适当地改变以达到相同的目标。

主分量分析(“PCA”)的有用推广是称为独立分量分析(“ICA”)的方法，当数据大大地不同于标准的多维高斯分布时，该方法是优选的。在该情况下，本发明方法相应地概括为使用窗口ICA(“WICA”)，随后使用残差分析——称为离群值检测——进行推广。在一个实施方式中，本发明在移动窗上使用PCA，随后是从主分量(“PC”)计算内积和数据残差，据信这不但在地震应用中而且在较宽的多维数据处理的领域中都是可有利应用的。这包括图像、语言和信号处理领域。

主分量分析(“PCA”)是用于数据分析的众所周知的经典技术，首先由Pearson(“On Lines and Planes of Closest Fit to Systems of Points inSpace,”Philos.Magazine v.2,pp.559-572(1901))提出，并且被Hotelling(“Analysis of a Complex of Statical Variables Into Principal Components,”Journal of Education Psychology v.24,pp.417-441(1933))进一步发展。据信首次已知的将主分量分析应用于地震数据以Karhunen-Loeve变换的形式进行，所述变换以Kari Karhunen和Michel Loeve命名(Watanabe,“Karhunen-Loeve Expansion and Factor Analysis,”Transactions of theFourth Prague Conference,J.Kozesnik,ed.,Prague,CzechoslovakiaAcademy of Science(1967))。该方法使用PCA描述一组地震道中的信息内容，形式是输入数据集为全部地震道，而不是可变大小的多维窗口。Watanabe的主要应用是分解全部地震道，并使用前几个主分量道重建最相干的能量，从而过滤掉非地质学噪音。

在地震分析中PCA最常用于将测量特征的数量减少至统计上独立的属性组(参见例如，Fournier&Derain,“A Staticla Methodology forDeriving Reservoir Properties from Seismic Data,”Geophysics v.60,pp.1437-1450(1995);和Hagen,“The Application of Principal ComponentsAnalysis to Seismic Data Sets,”Geoexploration v.20,pp.93-111(1982))。地震诠释过程常常从原始数据产生众多导数积。由于这些属性不同程度相关联，所以PCA已成为减少属性数量、同时保留大量信息的一流方法。

迄今，据信没有如此基于移动窗的统计离群值检测技术，该技术致力于在认真查看和勘测(scoping and reconnaissance)的基础上在地质学和地球物理数据中发现关注的地质特征。但是，这类技术已被应用至具体的地震数据的子集或域，用于专门的信号处理或油藏表征应用。Key和Smithson(“New Approach to S eismic Reflection Event Detectionand Velocity Determination,”Geophysics v.55,pp.1057-1069(1990))将PCA应用在叠前(pre-stack)地震数据中的2D移动窗上，并按比例缩放得到的本征值作为信号相干的度量。没有应用由主分量本身组成来检测叠前地震数据中的特征。Sheevel和Payrazyan(“Principal ComponentAnalysis Applied to 3D Seismic Data for Resevoir Property Estimation,”Society of Petroleum Engineers Annual Conference and Exhibition(1999))使用小的、1D移动垂直窗计算基于地震道（trace-based）的主分量，并将那些显得最地质学的PC输入预测远离井标定（well calibration）的油藏性能的分类算法。再一次地，该1D、单数据集方法没有尝试自动鉴定数据中的异常或离群值。Cho和Spencer(“Estimation ofPolarization and Slowness in Mixed Wavefields,”Geophysics v.57,pp.805-814(1992))以及Richwalski等人(“Practical Aspects of WavefieldSeparation of Two-Component Surface Seismic Data Based on Polarizationand Slowness Estimates,”Geophysical Prospecting v.48,pp.697-722(2000))使用频域中的2D窗口PCA以模拟预先限定数量的P-&S-波的传播。

Wu等人(“Establishing Spatial Pattern Correlations Between WaterSaturation Time-Lapse and Seismic Amplitude Time-Lapse,”PetroleumSociety’s 6^th Annual Canadian International Petroleum Conference (56^thAnnual Technical Meeting)(2005))的目标是使单次地震数据体或时间推移地震数据体与油藏模型中的流动模拟数据最佳地相关联以估计空间图形的实际饱和时间推移值(actual saturation time-lapse values)。其方法是进行点对点比较，其不是对原始数据体，而是对来自PCA分析的这些数据在第一主本征向量上的投影进行。因此，其目标是使地震数据与已知的模型相关联而不是鉴定地震数据中的异常图形。

授权给Bishop的美国专利5,848,379(“Method for CharacterizingSubsurface Petrophysical Properties Using Linear Shape Attributes,”(1998))公开了预测地下岩石性能和分类地震数据的方法用于表面或质地分析，不是在认真查看和勘测的基础上鉴定关注的地质特征——这是本发明致力于解决的技术问题。Bishop使用PCA进行统计分析以将地震道分解为正交波形基的线性组合——称作预先指定时间或深度间隔内的线性形状。线性形状属性(LSA)被定义为用于重建特定地震道形状的权重(或本征值)的子集。同样，Bishop没有公开同时分析多个数据体的交叠窗口，也没有公开使用统计分布来检测异常数据区域。

用于统计分析地质学和地球物理数据的其它方法已使用如人工神经网络、遗传算法和多点统计的方法，但目标不是自动检测异常图形。此外，这些方法一般具有有限的成功，因为它们的内运算（innerworking）常常是模糊的，并且它们常常需要并且高度依赖于大量的培训数据。

如前所述，PCA和ICA是常用于将高维(即，多变量或多属性)信号分离成统计上无关联的(即，独立的)分量的方法。本发明的窗口PCA和ICA将分量分析应用至通过将原始数据中的每个点(在本发明的一些实施方式中)表示为其相邻的点的集合(即，窗口)从原始数据衍生的数据集。为参照图3的流程图阐述该概念，在单个3维数据体上使用固定窗口大小进行WPCA在以下描述。同一程序或其ICA等同物可以被应用至2D数据，或同时应用至多个2D或3D数据体(参见图3的步骤31)。考虑大小N_x×N_y×N_z的3D地震数据体：

(步骤32)选择窗口形状(例如，椭圆体或立方体)和大小(例如，半径r，n_x×n_y×n_z)

3D地震数据体中的每个体素I_i，j,k，被表示为n_x×n_y×n_z维向量

其包括每个体素的窗口附近内的体素值。尽管在实施方式中使用的被描述，但是对于本发明不需要穷尽采样。可代替地使用其他采样策略，比如潜在随机采样或平铺(tiled)采样。

(步骤33)计算所有n维向量的平均矩阵和它们的协方差阵(n＝n_x×n_y×n_z)

(N=(N_x-n_x)×(N_y-n_y)×(N_z-n_z))如下：

$\stackrel{\‾}{\stackrel{\→}{I}}=\frac{1}{N}\underset{i,j,k}{\mathrm{\Σ}}{\stackrel{\→}{I}}_{i,j,k},$ $W=\frac{1}{N}\underset{i,j,k}{\mathrm{\Σ}}({\stackrel{\→}{I}}_{i,j,k}-\stackrel{\‾}{\stackrel{\→}{I}}){({\stackrel{\→}{I}}_{i,j,k}-\stackrel{\‾}{\stackrel{\→}{I}})}^T$

计算相关矩阵为

其中t和k是向量I的两个指数，因此代表三维的两个不同集的空间坐标。

(步骤34)计算

的本征值(主值){λ₁>λ₂>…>λ_n}和本征向量(主分量){v₁,v₂,…,v_n}。可选地，可以计算协方差阵的本征值；它们将与相关矩阵的本征值仅比例因子不同。这些本征向量大小为n_x×n_y×n_z，并且当从其向量形式变形回窗口形式时，表示数据中的不同的(独立的)空间图形(spatial patterns)，其顺序为从最常见的到最不常见的。相应的本征值表示每个本征向量所占的原始数据(即，方差的量)的多少。

产生一个或多个以下的部分的地震或属性数据体，然后检查其异常，这些异常可能从原始数据体来看并不明显：

(a)(步骤35)投影：可以使用每个主分量或主分量组(例如，选自聚类分析)重新产生的原始数据的一部分。这可以通过对每个主分量或主分量组的平均中心和归一化的地震数据体取内积获得。因此，向量A在向量B上的投影指proj(A)＝(A·B)B/|B|²并且是B方向的向量。

(b)(步骤36)残差：没有被第一k-1(即，最常见的)主分量捕获的原始数据体中的残留信号。在本发明优选的实施方式中，这通过将平均中心和归一化的地震数据体投影到由{v_k,v_k+1,…,v_n}生成的子空间来完成，所以

其中R是用户定义的介于0和1之间的阈值。可选地，可以颠倒地添加投影，但在大多情况下这将成为更大的计算负担。

(c)离群值：条目(b)的残差分析是在本发明的一个实施方式中测定每个体素的“异常程度”的方式。(a)和(b)的属性数据体在计算每个体素的“异常程度”的可选方式中并不需要，该“异常程度”被表示为R′(因为其与上述定义的残差R相关但不相同)，并通过下式给出：

$R_{i,j,k}^{\′}={(I_{i,j,k}-\stackrel{\‾}{I})}^T{\hat{W}}^{-1}(I_{i,j,k}-\stackrel{\‾}{I})$

使用该异常程度的量度，研究部分数据体。该量度还挑选位于前几个本征向量生成的空间内的“离群值”，但在一些情况下计算强度可以高于以上两个步骤。但是，可以注意到，在该情况下，以上步骤34可以被跳过，或者简单地替换为相关矩阵的楚列斯基分解，这确保R′的更快估算。

存在上述基本方法的变体，其使用不同的数据归一化方案。该方法可以扩展至任意数量的地震数据体。用户可以进行实验的可调节参数是(1)窗口形状，(2)窗口大小，和(3)残差投影的阈值，R。

在地震数据的2维切片上运用3×3WPCA的结果示于图1A-C中。图1A显示来自合成3D地震数据体的图像(2D时间切片)。在实际应用中，该显示典型地是彩色的，其中色彩指示地震反射振幅(例如，蓝色=正，红色=负)。图1B显示前十六个主分量占信息的90%之后原始图像的残差。在异常图形下残差具有高值，在这种情况下是有错误的。在图1B的彩色版本中，蓝色可能指示少量的残差，较暖的色彩可能突出异常错误系统，其现在在图1B的残差显示中清晰可见。在图1C中，以其30×30窗口形式，示出顶部(即前)十六个主分量14。在底端两行中可见几个主分量中的错误被捕获。

在2-维合成地震横截面上运用9×9WPCA的结果示于图2的示意性流程图中。在21，来自合成3D地震数据体的2D横截面被显示。色彩通常被用于表示地震反射振幅。太细微而不能被肉眼观察到的小的8ms的背斜（anticline）被包埋在背景水平反射率中。在22，输入图像的前四个主分量(本征向量)被显示。显示23示出原始图像在前四个本征向量上的投影，其占信息的99%。显示24示出从原始图像减去投射图像之后的残差。包埋的细微特征现在在介于地震道数字(一维测量横向位置)30-50之间约440ms的深度（双向移动时间）被展现。在彩色显示中，'热'色彩可能被用于展现包埋的细微特征的位置。

图3的流程图概述了本发明方法的实施方式，其中使用单个窗口大小将WPCA运用到多个数据体。

典范图形构造中的概括（推广）和效率

以下部分描述对本发明的窗口主分量分析的改进，其通过减少计算能确保更便利的适用性，以及通过解释主分量或独立分量及其选择性的保留或删除能够更好地使用结果。

计算效率：对于大的数据集来说，直接计算以上协方差阵的方法在内存和处理器要求方面都是计算繁重的。因此，本文公开的可选的方法利用这样的事实：PCA的单独向量是移动越过所有数据的窗口。例如，考虑具有值{I₁,I₂，...,I_N}的1-D数据集。为评估大小为K<N的窗口的协方差阵，表值（矩阵的元，entry）的平均和二阶矩(the mean and secondmoment)可以如下计算：

$E\left(X_i\right)={\stackrel{\&OverBar;}{X}}_i=\frac{1}{N-K}\underset{k=i}{\overset{N-K+i}{\mathrm{\&Sigma;}}}{I}_k$ 对于1≤i≤K

$E\left(X_i{X}_j\right)=\frac{1}{N-K}\underset{k=i}{\overset{N-K+i}{\mathrm{\&Sigma;}}}{I}_k{I}_{k+j-i}$ 对于1≤i≤j≤K

可以注意到，本方法仅包括采用数据的子向量的平均值和内积(较高维数的子矩阵)，因此避免储存和操作源自原始数据的诸多较小尺寸的窗口。因此计算方法的该改进允许具有有效阵列索引的面向对象软件(如Matlab和使用区域求和表，由Crow在“Summed-Area Tables forTexture Mapping,”Computer Graphics 18,207(1984)中描述的数据结构)以最小储存和计算工作（effort）计算协方差阵。

可选地，通过将协方差阵的计算表示为逐渐减小的区域上的一系列交叉相关操作，可以增进计算效率。为阐述该方法，考虑图4A-B所示的大小为n＝n_x*n_y的二维数据集和大小为m＝m_x*m_y的二维窗口。然后，通过首先计算每个数据样本的平均值，然后计算内积矩阵，再然后归一化该矩阵并减去平均值，可以获得相关矩阵W(t,k)。

首先，可以通过用由都等于1/(DS1中的像素数)的表值构成的数据样本(例如DS1)大小的核卷积(convolving)数据体来计算平均值。该运算的结果产生大的矩阵，但平均值是落在位于输出的左上角的大小m的窗口内的值。一般而言，该类型的运算被表示为corrW(核，数据)，并且其结果是以上读出的大小m的窗口。使用快速傅里叶变换(FFT)进行运算花费的时间与n*log(n)成比例，并且独立于采样窗口的大小。当m充分地大于log(n)时，该FFT方法比显式方法快。

其次，通过在数据集的子样本上进行一系列的corrW运算来计算内积矩阵U(t,k)。可以注意到，该矩阵的行i——表示为U(i，:)——可以被计算为U(i,:)=corrW(DSi,数据)。因此，以这种方式增加(populating)矩阵花费的时间与m*nlog(n)成比例或更好。但是，通过在数据的各个子区域进行几个corrW运算来计算U(t,k)是更有利的。特别地，可以改写

corrW(DSi,数据)＝corrW(数据,数据)-corrW(数据,DNSi)

其中corrW(数据,DNSi)表示DNSi与DNSi附近数据的交叉相关，所述附近数据在DNSi位置的m_x或m_y内。仅需要对所有行进行一次运算corrW(数据,数据)，然后corrW(数据,DNSi)需要被计算m次。优势来自于DNSi典型地远小于数据集的大小这一事实，所以corrW(数据,DNSi)是输入远小于corrW(数据，数据)的交叉相关。类似地，corrW(数据,DNSi)的计算可以被分解为对甚至更小的子区域的几个corrW运算。

对于不同的样本来说，DNSi的大部分是相同的，仅沿时间的一维采样窗口不同。例如，考虑图5A-B的图解。由A、B和C表示的图5A中的区域，一起采用形成没有被像素1采样的数据体的整个区域。即，区域可以被进一步细分以进行较少的计算。考虑由A和C生成的“垂直”区域，并比较图5B中所示的不同的采样区域DSi。通过联合几个较小的区域C1+C2+C3+C4+A1+A2生成类似的垂直区域(图5A中区域B的等同裂片是图5B中的联合B1+B2)。一般而言，仅m_x不同于这些可能的区域，每个区域相应于DSi的唯一横向位置。换言之，对于许多不同的数据样本DSi，A+C中包含的数据是相同的，所以其需要仅被操作m_x次——在该区域节约了m_y次计算。因此，corrW(数据,DNSi)的计算可以以该模式被最优化并根据以下进行计算，

corrW(数据,DNS1)=corrW(数据,A+C)+corrW(数据,B+C)-corrW(数据,C)

其中由字母表示的区域意指由该字母和数字标注的所有区域的联合；例如，方程式中的C指图5A中的区域C和图5B中的C1+C2+C3+C4，所以A+C由图5B中的A1+A2+C1+C2+C3+C4表示。由于corrW(数据,A+C)的计算需要对U(t,k)的m_y行进行仅一次，并且对于corrW(数据,B+C)相似处理，所以对于每行来说需要计算的唯一部分是corrW(数据,C)。效率增加来自于由C表示的区域典型地明显小于其它区域这一事实。以该方式进行的算法扩展至3-D数据集和窗口(并且实际上至任何维数)。

最后，通过适当地归一化矩阵U并减去平均值来获得交叉-相关矩阵W(t,k)。

W(t,k)=U(t,k)/nDS-平均值(DSt)*平均值(DSk)

其中nDS是每个数据样本中的元的数量。

掩码(Masks)的使用：对于非常大的数据集，甚至以上描述的计算效率可能并不足以用于可得的计算资源以便以适时方式产生结果。在这样的情况下，可以在预先限定的掩码上应用(a)用本征向量的内积计算或(b)主分量计算。掩码是对其进行计算的数据的空间子集。掩码可以(a)通过用户交互产生或(b)使用衍生属性自动产生。(b)的实例是使用梯度估计算法(gradient estimation algorithms)进行具有高的局部梯度的数据区域的预先选择。内积计算比主分量的计算更繁重，这促使在需要时对一个或两个计算运用掩码。

典范图形的应用

此外，计算的主分量/独立分量可以被聚类成代表通过结构（texture）、混沌或其它特征测量的类似图形的组。与残差数据体一起，原始地震数据在单独主分量或主分量组上的投影将产生多数个具有突出异常图形的衍生地震数据体。本发明方法的这些实施方式将在下文更详细地描述。

多重窗口/空间比例：此外，相较于同时直接计算它们的方式，在以分级序列针对多重嵌套窗口大小计算协方差阵中简化工作量是可能的。再次，考虑带有两个窗口大小K₁<K₂的一维实例。使用以上方法首先计算K₂的平均和二阶矩，随后K₁的同样的量可以如下计算：

$E_{K_1}\left(X_i\right)=\frac{1}{N-K_1}[E_{K_2}\left(X_i\right)*(N-K_2)+\underset{k=N-K_2+i}{\overset{N-K_1+i}{\mathrm{\&Sigma;}}}{I}_k]$ 对于1≤i≤K₁

$E_{K_1}\left(X_i{X}_j\right)=\frac{1}{N-K_1}[E_{K_2}\left(X_i{X}_j\right)*(N-K_2)+\underset{k=N-K_2+i}{\overset{N-K_1+i}{\mathrm{\&Sigma;}}}{I}_k{I}_{k+j-i}]$ 对于1≤i≤j≤K₁

注意，以上公式允许以增加的工作量计算较小窗口的量。该方法直接被扩展至较高维数的嵌套系列的窗口。

主分量和投影的利用：有许多可能的方式，其中本发明方法产生的主分量和投影可以被使用、组合和可视化。一个优选的实施包括使用残差鉴定异常，如上述。同样有效的方法是对选择的PC的子集进行原始数据的选择性投影。可以(a)由用户交互地或(b)使用PC上的计算矩阵自动地选择子集。(b)的实例可以是使用自动几何算法(automaticgeometric algorithm)选择这样的PC，所述PC具有类似于“通道”或管状结构的特征。另一实例可以是通过使用噪音检测算法或离差矩阵（dispersion matric）产生排出了“噪音”PC的投影来减少输入数据中的噪音。本领域技术人员将从该描述中联想到其它例子。

在各种窗口大小可视化投影结果的可选有用方法包括可视化(a)用户选择或自动选择的PC投影的组合，(b)各种残差阈值的残差，或(c)噪音分量。另一有用的变型包括可视化“分类数据体”，其包括用唯一确定在每个数据位置哪个PC投影具有最高值的颜色对该位置进行颜色编码。

迭代WPCA：已经发现图3中列出的工作流程产生的残差数据体在包含更多异常图形的区域中显示出较大值。因而，输入数据中的更细微图形常常被残差数据体中更明显的异常所掩盖。为增加WPCA对极细微图形的敏感度，可以使用两种可选的迭代方法：

迭代本征向量删除：该第一可选方法可以包括以下步骤：

1.进行图3流程图的前四个步骤(到本征向量和本征值产生)。

2.鉴定那些其投影重建大量的背景信号和最明显的异常的本征向量。

3.仅将数据投影在前述步骤中没有被鉴定的本征向量的子集上(在该投影图像中背景信号和最明显异常的信号应当被削弱)。

4.对前述步骤中产生的投影图形进行WPCA。

5.需要时重复步骤1-3。

迭代掩码或数据删除：该第二可选方法可以包括以下步骤：

1.进行图3的前四个步骤(到本征向量和本征值产生)。

2.通过检查各种残差数据体，鉴定输入数据中与最明显异常相对应的那些区域。

3.对数据进行WPCA，通过以下步骤排除那些鉴定的区域：

a.在WPCA分析之前将那些区域中的所有属性值设定为零，或者

b.当输入到WPCA时不包括那些区域。

4.对新的数据集进行WPCA。

5.需要时重复步骤1–3。

WPCA分类：基于投影强度，主分量可以被用于分类图像。这样的分类有助于通过便利的可视化，鉴定在所选的主分量中表示有特定图形的区域，尤其是当原始数据由多个数据体构成时。该变型可以包括以下步骤：

1.进行图3的步骤31-34(到本征向量和本征值产生)。

2.给数据中的每个点分配与本征向量相对应的数字，所述本征向量重建环绕该点的窗口中的最多信号。这构成了分类数据体，其中每个点包含介于1(即，第一个本征向量)和N＝n_x×n_y×n_z(即，最后一个本征向量)之间的数。

3.然后，通过给每个从1-N的值(或值的组)分配独特的颜色或透明度(或其组合)可视化分类结果。该方法是N-维图像的基于图形分类的一种形式。通过输出类别，仍基于投影图像中信号的大小，而非连续谱残差或投影值，该方法对细微特征的敏感度变大。

因此，本发明方法有利于从大的高维数据集如地震数据中提取特征。例如，将PCA应用至地震数据的大部分公开的方法与本发明的相似处仅在于它们对数据窗进行本征模式分解。一个实例是以上提到的Wu等的方法。在几个基本方面，他们的方法不同于本发明。首先，他们仅将小的1D垂直移动窗应用至地震数据作为对PCA的输入。仅对流程模拟数据（flow simulation data）使用3D移动窗。其次，仅第一个PC被用于重建时间推移地震数据和流程模拟数据。没有进行其它投影或数学组合如建立残差数据体。最后，没有尝试同时检查多个地震数据体，更不用说提取地震数据的固有图形(即，没有结合至预先存在的地质学模型)。

地震数据的扩散映射

从地震数据中提取地质学有意义图形的一个方法由计算在一些线性空间中数据的适当表示组成。通常地，这是主分量分析(PCA)的结果，借此数据变换成通过该方法获得的基础元的线性组合。地质学感兴趣的一些图形违背了数个假设——PCA利用：同等重要的图形可出现在不同的尺度中，它们的分布没有必要是高斯分布，并且在数据中收集它们的流形(manifold)可以不是线性的。这里描述的是解决所有的这些考虑同时保留PCA好处的方法。该方法基于所谓的R.R.Coifman等的Diffusion Map；见“Geometric diffusions as a tool for harmonic analysisand structure definition of data:Diffusion maps,”Proceedings of theNational Academy of Sciences,102(21),7426-7431(2005),其通过引用以其所有被允许的权限并入本文。与PCA一样，计算表示数据的基础的(图6中的61)。不像PCA，该基础是非线性变换的结果，其提供定义尺度概念的参数(ε)。因此，数据的非线性以可控的方式捕获。感兴趣地，尺度参数可被调整以产生与PCA的那些类似的结果，并且我们在此采用的归一化已经被A.Singer示出与独立分量分析(ICA)相关(“Spectralindependent component analysis”,Applied Computational HarmonicAnalysis,Anal.21(2006)135-144)。

进行地震数据扩散映射的步骤

在一种实施方式中，参考图7的流程图，扩散映射可用如下的基础步骤进行：

1.对于给定的2D或3D地球物理数据体(一个或多个)(71)和任何给定的采样策略，例如潜在随机采样、穷举采样或平铺采样，用样本设计{s_x,s_y，s_z}采样数据体(一个或多个)：每个样本由收集的数据点构成，如通过任意的、但尺寸固定的3-D窗口(例如，如果s_x＝s_y＝s_z，则为立方体)给出的，所述3-D窗口根据采样策略从一个位置移动至另一位置(步骤72)。

2.收集随机样本，即数据窗口向量,进入数据阵列A_m，n，其中m=1,...,M，M是每个样本中数据体素的数量(例如，对于矩形窗口，M=s_x×s_y×s_z)，n＝1,...,N，并且N是样本的数量，因此M<<N。

3.在步骤73，计算对称相似性矩阵L_i，j=exp(-||a_i-a_j||²/ε)，以便L=M×M矩阵，并且其中||...||指示选择的范数(norm)。这里a_i和a_j是来自数据阵列A_m，n和

的长度N的行向量。Epsilon(ε)是预先定义的尺度因子(比例因子)。

4.在步骤74，用

计算对角归一化矩阵D=Diag(D₁，...,D_I)。

5.在步骤75，通过归一化相似性矩阵计算扩散矩阵M=D^-1L。

6.在步骤76，计算对称扩散矩阵M_sym=D^1/2MD^-1/2。

地震扩散映射的使用

该对称归一化相似性矩阵可用于以下列方式进行数据分析：

1.图形分析：

a.在图7的步骤77，将M_sym通过eig(M_sym)分解成其本征值和本征向量，具有非零本征值的本征向量表示图形分析的尺度(ε)依赖性基础。在一些情况下，本征向量的子集完整定义感兴趣的图形。

2.异常检测(异常属性)

a.在步骤78，使用与当计算M_sym时相同的窗口，在数据体的所有可能位置收集样本

b.产生各处初始化为零的异常数据体。

c.在步骤79，对于异常数据体中每个类似位置n，设置值为

d.使用不同的尺度参数(ε)值产生至少一个异常数据体，并且观察异常中尺度依赖性差异；图6图解在不同尺度下产生的两个异常数据体(62)。

前述应用为阐述本发明的目的而涉及本发明的特定实施方式。但是，本文所述实施方式的许多改进和变化是可能的，这对于本领域技术人员来说是明显的。所有这些改进和变化意图在本发明的范围内，如在所附权利要求中限定的。本领域技术人员将认识到对于实际应用，本发明方法的至少一些步骤必须使用根据本文教导编程的计算机执行。

Claims (31)

1.用于在代表地下区域的一个或多个地球物理数据或数据属性的离散集——每个这样的数据集被称作“原始数据体”——中鉴定地质特征的方法，所述方法包括下述，至少其中之一是使用计算机执行的：

(a)选择数据窗口的形状和大小；

(b)对于每个原始数据体，将所述窗口移动至多个位置，并为每个窗口位置形成数据窗口向量，其分量由来自该窗口内的体素值构成；

(c)进行所述数据窗口向量的统计分析，在多个原始数据体的情况下所述统计分析共同进行；

(d)使用所述统计分析以鉴定所述数据中的离群值或异常；和

(e)使用所述离群值或异常预测所述地下区域的地质特征。

2.根据权利要求1所述的方法，其中所述统计分析使用由以下组成的统计分析技术组中的一个来进行：

(i)计算所有所述数据窗口向量的平均矩阵和协方差阵；

(ii)独立分量分析；

(iii)使用聚类方法聚类所述数据；

(iv)扩散映射；和

(v)另一统计分析方法。

3.根据权利要求2所述的方法，其中使用(i)进行统计分析，进一步包括使用主分量分析。

4.根据权利要求3所述的方法，其中所述协方差阵的本征值和本征向量被计算，所述本征向量是相应原始数据体的主分量集；

并且其中步骤(d)和(e)包括将原始数据体投影在选择的本征向量子集上以产生部分投影数据体，所述本征向量子集基于其相应的本征值选择，并且测定残差数据体，所述残差数据体是在所述投影数据体中没有被捕获的原始数据体的一部分；然后鉴定所述残差数据体中的异常特征，并使用它们预测所述地下区域的物理特征。

5.根据权利要求1所述的方法，其中所述数据窗口是N-维的，其中N是整数以使1≤N≤M，其中M是数据集的维数。

6.根据权利要求3所述的方法，其中所选择窗口大小和形状的所述平均矩阵和所述协方差阵使用补充窗口进行计算，其中相应于在(a)选择的窗口中的每个位置的补充窗口表示当所述窗口移动通过原始数据体时在该位置出现的数据值集。

7.根据权利要求4所述的方法，其中所选择的子集基于通过结构、混沌或其它数据或几何学属性测量的图形的内在相似性进行选择。

8.根据权利要求4所述的方法，其中所选择的所述本征向量子集通过以下确定：按照从最大到最小的顺序对本征值求和，直至最大N个本征值的和除以所有本征值的和超过预先选择的R值，其中0<R<1，然后选择与N个最大本征值相关的N个本征向量。

9.用于从代表地下区域的地球物理数据或数据属性的2D或3D离散集(“原始数据体”)鉴定地质特征的方法，包括下述，至少其中之一是使用计算机执行的：

(a)选择数据窗口的形状和大小；

(b)将所述窗口移动至所述原始数据体中多个交叠或不交叠的位置，以使每个数据体素被包括在至少一个窗口中，并为每个窗口形成数据窗口向量I，其分量由来自该窗口内的体素值构成；

(c)计算所有所述数据窗口向量的协方差阵；

(d)计算所述协方差阵的本征向量；

(e)将所述原始数据体投影在选择的所述本征向量的子集上以产生部分投影数据体；和

(f)鉴定所述部分投影数据体中的离群值或异常，并使用它们预测所述地下区域的地质特征。

10.根据权利要求9所述的方法，其中产生部分投影数据体的所述选择的本征向量的子集通过基于与其相关的本征值消除本征向量来确定。

11.根据权利要求9所述的方法，其中所述选择的本征向量的子集由用户交互地选择或基于自动化鉴定的噪音或几何特征进行选择。

12.根据权利要求9所述的方法，其中所述选择的本征向量的子集通过以下步骤来确定：设计用于确定所述原始数据体中的明显异常的标准，使用所述标准选择一个或多个明显异常，以及鉴定一个或多个本征向量，所述本征向量的相关数据分量(所述原始数据体在所述本征向量上的投影)对所选择的明显异常有贡献，或比预先设定量的背景信号占更多，然后选择保留的本征向量的一些或全部；其中步骤(f)能够发现比用于确定所述选择的本征向量的子集的所述明显异常更细微的异常。

13.根据权利要求12所述的方法，进一步包括在步骤(e)后使用所述部分投影数据体代替所述原始数据体，重复步骤(a)-(e)，产生更新的部分投影数据体，其然后被用于步骤(f)。

14.用于在代表地下区域的地球物理数据或数据属性的2D或3D离散集(“原始数据体”)中鉴定地质特征的方法，包括下述，至少其中之一是使用计算机执行的：

(a)选择数据窗口的形状和大小；

(b)将所述窗口移动至所述原始数据体中多个交叠或不交叠的位置，以使每个数据体素被包括在至少一个窗口中，并为每个窗口形成数据窗口向量I，其分量由来自该窗口内的体素值构成；

(c)计算所有所述数据窗口向量的协方差阵；

(d)计算所述协方差阵的本征值和本征向量；

(e)选择计算体素的异常程度的方法，并使用其确定由计算的比预先确定的阈值更异常的体素构成的部分数据体；和

(f)鉴定所述部分数据体中的一个或多个异常特征，并使用它们预测所述地下区域的地质特征。

15.根据权利要求14所述的方法，其中通过x，y，z指数i,j,k表示的体素的所述异常程度R′由下式计算：

$R_{i,j,k}^{\&prime;}={(I_{i,j,k}-\stackrel{\&OverBar;}{I})}^T{\hat{W}}^{-1}(I_{i,j,k}-\stackrel{\&OverBar;}{I})$

其中I_i，j，k是来自(b)的数据窗口向量的分量，其包括体素i,j,k；

$\stackrel{\&OverBar;}{\stackrel{\&RightArrow;}{I}}=\frac{1}{N}\underset{i,j,k}{\mathrm{\&Sigma;}}{\stackrel{\&RightArrow;}{I}}_{i,j,k},$ $W=\frac{1}{N}\underset{i,j,k}{\mathrm{\&Sigma;}}({\stackrel{\&RightArrow;}{I}}_{i,j,k}-\stackrel{\&OverBar;}{\stackrel{\&RightArrow;}{I}}){({\stackrel{\&RightArrow;}{I}}_{i,j,k}-\stackrel{\&OverBar;}{\stackrel{\&RightArrow;}{I}})}^T;$

其中离散的原始数据体由N_x×N_y×N_z体素构成，所选择的窗口形状和大小为n_x×n_y×n_z体素，并且N=(N_x-n_x)×(N_y-n_y)×(N_z-n_z)。

16.根据权利要求14所述的方法，其中所述异常程度通过将所述原始数据体投影在选择的本征向量的子集上以产生部分投影数据体，以及确定残差数据体来确定，所述本征向量的子集基于其相应的本征值进行选择，所述残差数据体是在投影数据体中没有被捕获的所述原始数据体的一部分，所述残差数据体是(f)中用于预测所述地下区域的物理特征的部分数据体。

17.根据权利要求14所述的方法，其中所述异常程度通过将所述原始数据体投影在选择的本征向量的子集上以产生用在(f)中的所述部分数据体来确定。

18.根据权利要求1所述的方法，其进一步包括使用所预测的所述地下区域的地质特征推测石油潜能或其欠缺。

19.用于从地下区域开采烃的方法，包括：

(a)获得所述地下区域的地球物理勘测数据；

(b)至少部分地基于使用权利要求1所述的方法对地球物理勘测数据鉴定的所述区域的物理特征获得所述地下区域石油潜能的预测，权利要求1通过引用而被并入此处；

(c)响应于石油潜能的肯定性预测，在所述地下区域中钻井并开采烃。

20.根据权利要求9所述的方法，其中通过在所述数据体的逐渐变小的区域上计算一系列交叉相关运算来进行所述协方差阵的计算。

21.根据权利要求2所述的方法，其中使用(i)进行统计分析，并通过在每个窗口的逐渐变小的区域上计算一系列交叉相关运算来进行所述协方差阵的计算。

22.根据权利要求2所述的方法，其中所述统计分析技术是(iv)，包括使用非线性变换用基础集表示所述数据窗口向量，其中所述非线性变换包括定义尺度概念的参数。

23.根据权利要求22所述方法，其中有N个数据窗口位置并因此有N个数据窗口向量，每个具有M个分量，即在每个数据窗口位置处具有M个数据体素；并且其中所述扩散映射技术包括：

收集所述N个数据窗口向量

进入数据阵列A_m，n，其中m=1,...,M并且n=1,...,N；

计算M×M相似性矩阵L，其中L_i.j是a_i和a_j之间差异的量度，其中a_i和a_j是来自所述数据阵列A_m，n和

的长度N的行向量；并且其中L_i.j包括使用者选择的尺度参数；

从所述相似性矩阵形成对角矩阵D，其中

通过归一化所述相似性矩阵计算扩散矩阵：M=D^-1L；

计算对称扩散矩阵M_sym=D^1/2MD^-1/2；和

使用所述对称扩散矩阵用于所述原始数据体的图形分析或异常检测，并将一种或多种地质特征与一个或多个图形或异常相关联。

24.根据权利要求23所述的方法，其中

L_i，j＝exp(-||a_i-a_j||²/ε)

其中ε是使用者选择的尺度参数，并且||...||指示选择的范数。

25.根据权利要求24所述的方法，进一步包括为至少一个另外的ε选择进行所述统计分析。

26.根据权利要求22所述的方法，其中M<<N。

27.根据权利要求23所述的方法，其中使用所述统计分析以鉴定所述数据中的离群值或异常包括计算所述对称扩散矩阵的本征向量和本征值，并在图形分析中使用它们。

28.根据权利要求23所述的方法，其中使用所述统计分析以鉴定所述数据中的离群值或异常包括：

使用所选择的数据窗口形状和大小，从所述原始数据体收集所有可能的样本；

将所述数据样本投影在所述对称扩散矩形的逆，从而产生尺度依赖性异常属性数据体；和

鉴定所述尺度依赖性异常属性数据体中的离群值或异常。

29.根据权利要求1所述的方法，其中所述多个位置通过选自潜在随机的、穷举的和平铺的采样策略确定。

30.根据权利要求1所述的方法，其中：

所述数据窗口可移动至交叠位置；

在所述原始数据体中的每个数据体素包括在至少一个窗口中；并且

数据值的分布从所述统计分析计算并且用于鉴定所述数据中的所述离群值或异常。

31.根据权利要求30所述的方法，其中所述数据中离群值或异常的鉴定包括(i)计算所述数据值分布中每个数据窗口的发生概率或等价度量；和(ii)鉴定低概率数据区域作为可能的离群值或异常。

Images