专利名称：用于地震应用的参量fk技术的制作方法
背景技术：
阵列记录用于测量波场，并提取关于波场传播通过的介质的信息。在地震学中，阵列分析已用于地震散射分布估计，表面波扩散测量，与火山行为有关的地震波场测量，地震断裂延伸估计，以及其他分析。
通常根据是否按频域或时域进行，可将阵列分析技术分为两类。特别是，频域方法能够分辨经多重路径或使用多重源的情形中同时到达阵列的多重信号。
一种最简单的频率-波数方法是波束成形或传统方法。可将所估计的功率谱显示作为阵列响应与真实功率谱的卷积。此处，阵列响应完全由传感器的数量和空间分布来决定。由于阵列响应对参考信号的波数和旁瓣而言具有较宽的波束宽度，阵列响应使传统方法所估计的功率谱变得模糊。
已提出了最小偏差无失真(minimum variance distortionless)方法，用来产生对单个识别波而言比传统方法所产生分辨率明显更高的分辨率。不过，当存在多个波时，最小偏差无失真方法的分辨能力与传统方法等同。为分辨多个信号，提出了信号子空间方法，称为多重信号特征(MUSIC，multiple signal characterization)(如参见Johnson和Dudgeon的著作，于Prentice-Hall公司出版的Signal ProcessingSeries，ISBN 0-13-048513-6)。MUSIC的关键性能在于，将傅里叶系数的空间交叉协方差矩阵(spatial cross-covariance matrix)中的信息分成两个向量空间—信号子空间和噪声子空间。由于不将真实谱与阵列响应进行卷积，MUSIC方法允许以高分辨率估计功率谱。
在地震分析技术中，需要用于产生高分辨率频率-波数谱和更高纬度扩展频率-波数-波数谱的方法。本发明的方法能够满足这种需求。
发明内容
在地震数据分析中，非常需要使用少量传感器来产生高分辨率频率-波数(fk)谱和更高纬度扩展频率-波数-波数(fkk)谱，这样的技术具有广泛用途。这些应用包括本地速率测量和扩散波检测。当可获得较少样品时，与标准fk变换相比，参量方法提供较高的波数分辨率。诸如多重信号分类(MUSIC)方法的参量方法可用于在每个频率处的波数-波数或波数估计，以产生算法fkk-MUSIC和fk-MUSIC以及用于地震应用中的多种其他参量算法。另外，可将参量方法扩展到fkkk或稀疏/不规则阵列。经证明，当估计本地速率涉及到自有限数量接收器的数据的任何时候，本发明的技术都作为有用的工具。
这样，本发明提供了包括瞬时傅里叶变换和基于参量算法的空间估计的地震分析方法。在本发明的一个实施例中，参量算法为MUSIC。在本发明实施例的一个方面，将参量算法应用于等距离空间采样的数据。在本发明的另一方面，将参量算法应用于稀疏空间采样的数据。在本发明的该实施例中，将数据看作MN×L的矩阵馈送到参量算法中，其中，MN为二维阵列传感器的数量，L为每传感器采样数。此外，在本发明的该实施例中，将傅里叶变换应用到每个迹线(trace)，传递到参量算法的参量由傅里叶变换级，采样频率，在x方向上接收器之间的距离，在y方向上接收器之间的距离，number_of_sources，或f_span的组中选出。
本发明的其他实施例应用其他参量方法，包括，但不限于，最小偏差算法，本征向量算法，最大似然算法，Pisarenco算法，ARMA(自回归移动平均算法)，AR或MA算法以及最大熵算法。
以上给出了本发明的概括性描述，参照在该说明书中所描述的以及在附图中所给出的本发明的实施例，可给出关于本发明更为具体的描述。不过应注意，说明书和附图仅说明本发明的某些实施例，因此，不应将其视为对本发明范围的限制。本发明可包括同样等效的实施例。
图1表示本发明方法的一个实施例的流程图。
图2表示kx-ky的快速傅里叶变换基波束成形估计不能区分频率在9Hz处的两个事件。
图3表示参量方法MUSIC能够区分图2中所观测的两个事件。近似速率与在合成数据中两个事件的速率相一致。阵列由2×2的元素组成。
图4表示两个迹线的fk图。
图5显示出图4中所测同样两个迹线的fk-MUSIC估计。
具体实施例方式
现将参照本发明的示例性实施例详细描述本发明。尽管结合这些实施例会描述本发明，不过应理解，所述实施例并不意在将本发明完全和具体限定在仅描述的这些实施例。相反，本发明意在覆盖在所附权利要求所限定本发明精神和范围内中所包括的替代，修改，和等效为生成诸如fkk-MUSIC和fk-MUSIC之类的算法以及适用于地震应用中的其他算法，本发明对波数-波数或波数估计应用了参量方法，如MUSIC算法，最小偏差算法，本征向量算法，最大似然算法，Pisarenco算法，ARMA，AR或MA算法以及最大熵算法。如此处所使用的，术语“参量算法”指描述信号的数学模型，且该模型使用所获得数据确定模型的参量。类似于功率谱密度估计(PSD)，以及二次统计如自相关(ACF)的传统非参量谱估计方法假设，在观测窗口之外的数据为零。这种假设限制了谱估计的分辨率。另一方面，参量方法根据模型对观测窗口之外的数据提出假设，依据该方法，对未观测到的数据进行多种形式的描述。首先，对未观测到的数据的描述是根据关于观测数据的现有信息和用于产生观测数据的过程，其次，所估计的参量随假设模型而定(如参见S.Kay的Modem SpectralEstimation，于Prentice-Hall公司出版的Signal Processing Series，ISBN 0-13-598582-X)。
图1表示本发明方法的一个一般实施方式的流程图(100)。首先，在步骤110，从2维阵列的传感器获得地震迹线，在步骤120，对迹线应用傅里叶变换。其次，对频率F选择频率范围(106)，并产生空间谱相关矩阵(108)，然后，通过本征向量的方法分解出噪声和信号子空间(110)。在接下来的步骤，在波数域对参量方法应用以波束成形(112)，并生成对于具体频率F的二维参量谱(114)。此时，可选择另一频率F(116)，对此新频率F应用步骤106，108，110，112和114。或者，从114生成的参量谱产生谱立方(spectrum cube)，最后执行速率分析(120)。
通常，高分辨率fk谱或fkk谱需要较宽的时间和空间窗口(记录时期)以及合适的采样。在地震数据收集中，时间采样和记录时期通常不成问题，且采样频率通常比关心的频率高得多。然而空间采样通常受到限制。此外，有时需要减少在fkk或fk变换中因数据的空间不稳定特性造成的空间扩展以及对局部估计分析的必要。本发明对参量方法进行修改，以产生具有空间孔径明显小的高分辨率算法fkkk-MUSIC和fk-MUSIC。
首先，应考虑二维(2D)散布(spread)。二维散布可用于小kx-ky估计。一维(1D)散布或线性阵列可视为二维散布的特殊情形。此处，假定了均匀的二维阵列具有对于特殊情形N＝1的M×N个元素；实际上为包括M个元素的一维阵列。不过，可将元素扩展到稀疏或不规则二维和一维阵列。当将参量技术应用于有限孔径以及接收器非常少的散布时，与传统方法如快速傅里叶变化(FFT)基波束成形技术相比，产生极佳的kx-ky估计。使用fkk-MUSIC或一维fk-MUSIC方法的一个优点在于，这些具体图所提供事件的高分辨率。
本发明用于波数-波数(kk)估计的参量技术是基于由到达方向(DOA，Direction Of Arrival)技术所使用的类似应用(参见Nikias和Petropulu的著作，于Prentice-Hall公司出版的Signal ProcessingSeries，ISBN 0-13-678210-8)。fk-MUSIC技术与DOA技术之间的不同在于，在DOA技术中，以步进通过可能的到达角度来执行波束成形，而在本发明中，波束成形步进通过波数k空间。而且，在DOA技术中，假设在介质中的速率近似为常数，从而，空间相关矩阵使用时间数据。不过，在本发明中，为分辨每个频率的波数，在波数波束成形之前，建立空间谱矩阵。这允许产生fkk-MUSIC和fk-MUSIC估计—或任何其他fkk-或fk-算法，其中，使用参量方法对k进行估计。
瞬时时刻二维阵列的空间测量所获得的二维傅里叶变换将空间采样的数据转换成波数-波数(kk)域，各轴由相对应的空间采样来确定。如果在时间上还持续进行测量，三维傅里叶变换将数据转换成fkk域。
波数轴跨越负正波数。相对于通过阵列中点的法线(normalline)，正和负波数表示到达方向相对。由空间采样时期，即在每个方向上相邻传感器之间的距离，来确定波数轴的奈奎斯特空间。对于在x和y方向展开的均匀矩形二维阵列，波数轴表示由处在这些方向上的传感器所测量的表观(apparent)波数。
表观波数轴覆盖[-1/(2dx) 1/(2dx)]和[-1/(2dy) 1/(2dy)]的范围，其中，dx和dy为在x和y方向上相邻传感器之间的距离。如果传感器不等间距，则dx和dy为在这些方向上所发现的任何两个相邻传感器之间的最短距离。
传统波束成形技术的分辨率，主瓣，由孔径的尺寸而决定。在本发明中，使用了近代谱估计技术来提高波数估计的分辨率。分辨率的这种提高等效于有效增大孔径尺寸。
相关矩阵包括两个主要类空间-时间和空间-谱。在空间-时间相关矩阵中的元素Rij定义属于传感器i的时间样品与传感器j的时间样品之间的相关性。空间-谱相关矩阵通常通过计算相关矩阵前首先进行每个传感器时间采样序列的傅里叶变换来构建。于是，在该矩阵中的每个元素定义元素i和j的交叉谱。因而将传感器之间的相关延迟转化成在傅里叶域中与传感器有关的相移。
Y＝[Y11Y21...Yi1...YM1...Y12Y22...Yi2...YM2...Y1NY2N...YiN...YMN]T
Yij向量为属于二维阵列中ij传感器的时间样品的快速傅里叶变换(FFT)。从而，有R(ω)＝Y(ω)Y’(ω)。Y’为Y的共轭转置。R(ω)为MN×MN的矩阵，它是角频ω的函数。以e’Re形式应用操纵向量e，其中，e’为e的共轭转置，e向量由e＝[exp(-jk.r0)...exp(-jk.rMN-1)]T定义，推导出在k处的空间傅里叶变换操作符，二次方程形式e’Re。r定义在xy平面中传感器的位置。由于P(k，ω)＝e’Re＝e’YY’e＝|e’Y|2，该量称为二次项。P(k，ω)为在k和ω处的功率估计。对于k＝(kx，ky)，操纵向量e由下式定义e(kx，ky)＝[exp(-j(kxx0+kyy0)......exp(-j(kxxMN-1+kyyMN-1))]T由于地震源也为宽频带的，必须对每个频率成分的波数进行估计。在本发明中，选择空间-谱相关矩阵代替空间-时间相关矩阵，以便开发出当估计波数kx和ky时在具有现代谱估计方法的结果的空间中取代傅里叶变换的技术。
由于相关矩阵具有厄密共轭对称性(R’＝R)，且通常为正定矩阵(对于x≠0，x’Rx＞0)，相关矩阵的本征向量相对该矩阵正交，且本征值通常为正vi’Rvj＝λiδij且vi’vj＝δijR可由其本征系统定义R＝∑ij＝1...MNλivi’vj且R-1＝∑ij＝1...MNλi-1vi’vj从而，对于二维阵列，以及MN个接收器，可将传统波束成形器的受操纵响应写为以下形式e’Re＝∑i＝1...MNλi|e’vi|2在本征分析中，按照其幅值的顺序将本征值排序。然后，将最低的本征值归因于噪声。通过将排好序的值绘图并找出与噪声下限相对应的阈值，从而实现将本征值归入与信号有关或与噪声有关本征值组中的分类。在执行本征分析之前执行常规的傅里叶操作(空间)可有助于确定该阈值。然后将相应的本征向量按其本征值排列。最大的本征值还受噪声等级的影响。不过，与最大本征值相对应的本征向量仅通过信号传播的方向来确定。从而，将由所有本征向量所定义的空间划分为噪声和信号子空间。由于与信号相对应的本征值的幅度还受噪声等级的影响，通常将信号子空间称为信号加噪声子空间。在这里，将信号加噪声子空间称为信号子空间，这是为避免任何模糊，特别是与信号幅度相比更关心波数的估计时尤为如此。
根据相关矩阵R的本征值和本征向量，可将最小偏差谱写成如下形式PMV(k)＝[∑i＝1...Nsλi-1|e’(k)vi|2+∑i＝Ns+1...MNλi-1|e’(k)vi|2]-1Ns为信号子空间中不同波数的数量。在上述表示中，第二个求和包括噪声子空间。当使e(k)朝向正确的传播波数时，第一个求和值增加而第二个求和值下降。由于所关心的是括号内表达式的互易性，对第一个表达式调零提高了与传播波数有关的那些峰值。下面说明在该表达式中的第二个“求和”的行为。当e(k)与正确的波数调准时，第二个求和变得微不足道，从而其互易性变得显著。使用DOA术语，其原因是信号传播的任何方向与噪声子空间理想正交。从而本征向量方法提供PMV(k)＝[∑i＝Ns+1...MNλi-1|e’(k)vi|2]-1通过白化噪声子空间以及通过用1替换属于噪声子空间的所有本征值，创建MUSIC参量方法(再次参看Johnson和Dudgeon的著作，于Prentice-Hall公司出版的 Signal Processing Series，ISBN0-13-048513-6)。由该方法背后的假设可知，谱显示出对所有信号相等的幅度以及相等的白噪声下限；从而，丢失有关幅度的任何信息。
PMUSIC(k)＝[∑i＝Ns+1...MN|e’(k)vi|2]-1因此，本发明提供了用于地震应用的可替代三维傅里叶变换的算法。它以基于参量技术的估计取代空间傅里叶变换的同时，保持了时间的傅里叶变换。结果为使用类似fkk-MUSIC和特殊情形一维阵列技术fk-MUSIC的参量方法的fkk估计。
将数据作为MN×L的矩阵馈送到算法，其中，MN为接收器(迹线)的数量，L个传感器的采样数量。在该算法中，首先在每个迹线上实施傅里叶变换操作，分辨频率特征。将快速傅里叶变换(“NFFT”或FFT的数量)参量传递到算法中。示例为8192或512等(即，2n，n为整数)。
传递到算法的其他参量为“Fs”(采样频率)，dx(在x方向上接收器之间的距离)，dy(在y方向上接收器之间的距离)，“number_of_source”和“f_span”。在其当前形式中的算法允许等间距空间采样；不过，易于将该算法扩展到应用以稀疏空间采样以及多种几何方式。
当傅里叶变化操作分辨出时间频率后，必须分辨出波数。作为参量传递到算法的分辨率f_span决定df或在时间频率域中对于波数估计所选出的分辨率。在df范围内估计波数谱。如果选择NFFT为8192，可将f_span选择为例如10，将奈奎斯特空间划分为4096/10个部分。这说明奈奎斯特空间分辨率df＝10/4096。如果选择NFFT为512，则对f_span参量选择值为2更为合适，从而获得奈奎斯特空间的df＝2/256。在该[f f+df]频率范围内分辨波数。如果df过大，在此频率窗口内估计波数的工作会导致不精确的估计。为生成集中所有波束成形技术的空间-频率相关矩阵，使用来自每个迹线的[f f+df]。
参量“number_of_sources”说明算法预期有如何多个不同的波数。该参量不会比迹线数量更大。如果选择“number_of_sources”为1，表示认为仅1个本征值属于信号子空间。因此，修改相关矩阵，在通过不同的波束成形技术使用修改的相关矩阵。
该算法在不同阶段可生成fkk-MUSIC，fkk-最小偏差，fkk-最大似然和fkk-本征值估计。方法的选择例如依据信噪比，所用传感器数量，噪声类型及应用。除传统波束成形方法以外，当有过多的接收器时，输出可导致多个伪峰值。对于多个接收器，fkk-操作或相关传统波束成形的方法可能最佳。
示例fkk-参量(fkk-MUSIC)的应用fkk方法在每个频率处产生kx-ky分布。出于显示目的，给出在单个频率处的fkk结果截面。
此处给出的结果表示，诸如MUSIC之类的参量方法如何能够区分传统波束成形技术(快速傅里叶变化基)所不能区分的两个独立的到达。所用合成数据的两个事件以200m/s和500m/s到达。使用2×2的阵列。对于到达所用的小波为Ricker小波。将接收器之间的距离dx和dy设置为12.5米。在图2中表示波束成形技术的结果。从该图中显然不能区分这两个事件。另一方面，图3表示对该合成数据应用MUSIC参量技术的结果。该技术容易地区分出这两个事件。看到阵列仅包括2×2的接收器，该结果非常令人兴奋。
参量方法(例如，fkk-MUSIC，fk-MUSIC)例如在本地速率分析以及波数估计中能够找到广泛的应用。此外，可将该方法扩展到fkkk-或稀疏/不规则阵列。对于一维阵列的特殊情形，应用fk-MUSIC。图4表示两个迹线的fk图。图5显示出同样两个迹线的fk-MUSIC。通过比较，显然看出分辨率的不同。
尽管参照特殊实施例描述了本发明，本领域技术人员应该理解，在不偏离本发明精神和范围的前提下，可对本发明进行多种变化以及可存在等效替代。另外，为适应具体情形，材料，或过程以达到本发明目的，精神和范围，可进行多种修改。本发明的范围意在覆盖所有这样的修改。
在此引述的所有参考资料有助于理解本发明，在此引作参考。
权利要求
1.一种地震分析的方法，包括获取地震数据，对地震数据应用时间傅里叶变换以产生变换的数据，基于参量算法对变换的数据应用空间波数估计，以及从参量算法生成波数谱。
2.根据权利要求1的方法，其中，第二应用步骤包括选择在频率F周围的频率范围。
3.根据权利要求1的方法，其中，第二应用步骤包括产生空间-谱相关矩阵。
4.根据权利要求3的方法，其中，第二应用步骤还包括将空间-谱相关矩阵分解成噪声和信号子空间。
5.根据权利要求4的方法，其中，使用本征值分解方法执行分解步骤。
6.根据权利要求5的方法，其中，对本征值具有幅度阈值。
7.根据权利要求1的方法，其中，参量算法为MUSIC算法。
8.根据权利要求1的方法，其中，参量算法为最小偏差算法，本征向量算法，最大似然算法，Pisarenco算法，ARMA，AR或MA算法，或最大熵算法。
9.根据权利要求1的方法，其中，第二应用步骤包括选择在频率F周围的频率范围，产生空间-谱相关矩阵，将空间-谱相关矩阵分解成噪声和信号子空间，并且其中，参量算法为MUSIC算法，最小偏差算法，本征向量算法，最大似然算法，Pisarenco算法，ARMA，AR或MA算法，或最大熵算法。
10.根据权利要求1的方法，其中，对来自等间距空间采样的数据应用参量算法。
11.根据权利要求1的方法，其中，对来自稀疏空间采样的数据应用参量算法。
12.根据权利要求11的方法，其中，传递到参量算法的参量从由傅里叶变换级，采样频率，在x方向上接收器之间的距离，在y方向上接收器之间的距离，number_of_sources，或f_span构成的组中选出。
13.根据权利要求12的方法，其中，作为参数传递到参量算法的f_span确定分辨波数谱的时间频率域。
14.根据权利要求12的方法，其中，将number_of_sources传递到参量算法。
15.根据权利要求14的方法，其中，number_of_sources不大于信号数。
16.根据权利要求1的方法，其中，从一维传感器阵列获得地震数据。
17.根据权利要求16的方法，其中，将数据看作M×L的矩阵馈送到参量算法中，其中，M为一维阵列传感器中信号的数量，L为每传感器采样数。
18.根据权利要求16的方法，其中，所生成的波数谱对于每个频率范围为一维fk谱。
19.根据权利要求18的方法，包括生成fkk谱三维量的附加步骤。
20.根据权利要求1的方法，其中，从二维传感器阵列获得地震数据。
21.根据权利要求19的方法，其中，将数据看作MN×L的矩阵馈送到参量算法中，其中，MN为二维阵列传感器中信号的数量，L为每传感器采样数。
22.根据权利要求19的方法，其中，所生成的波数谱对于每个频率范围为二维fkk谱。
23.根据权利要求21的方法，包括生成fkk谱三维量的附加步骤。
24.根据权利要求1的方法，其中，从三维传感器阵列获得地震数据。。
25.根据权利要求24的方法，其中，将数据看作MNO×L的矩阵馈送到参量算法中，其中，MNO为三维阵列传感器中信号的数量，L为每传感器采样数。
26.根据权利要求24的方法，其中，所生成的波数谱对于每个频率范围为三维fkkk谱。
27.根据权利要求26的方法，包括生成fkkk谱四维图像的附加步骤。
28.根据权利要求1的方法，包括执行速率分析的附加步骤。
全文摘要
本发明提供了在地震数据处理中用于解决在从有限数量接收器的数据获取期间所遇到问题的方法。本发明应用参量方法提供较高波数分辨率。
文档编号G01V1/36GK1643402SQ03804106
公开日2005年7月20日 申请日期2003年3月19日 优先权日2002年3月27日
发明者凯毕兹·伊兰普尔 申请人:维斯特恩格科有限责任公司