处理地震数据的高分辨率拉冬变换 【技术领域】
本专利申请要求于2002年5月24日提交的、名称为“处理数据的高分辨率拉冬变换”(High Resolution Radon Transform for ProcessingData)的美国临时专利申请系列号60/383209的优先权权益。
背景技术
本发明涉及用于地震数据处理的高分辨率拉冬变换的计算。
拉冬变换，与傅立叶以及其他几种变换一起，是地球物理学家可用来模拟和分析地震信号的工具的一部分。在地球物理学和其他应用领域中，通过克服由于数据的采样和噪声内容造成的限制，改进的拉冬变换已被获得。
用于地球物理学的最早的拉东变换仅仅是在大的间隔上采样的连续函数得到的结果的离散型式。稍晚引入的离散拉冬变换算法(Beylkin，1987，Hampson，1986)提供离散采样、孔径限制数据的完全可逆变换。现在被称为常规的这些变换，对分离信号和噪声是有效的，假如输入数据和它们地线积分(即拉冬变换域中的数据)被采样而不混淆。实际上，被混淆的输入数据的处理是一个重要问题，而对现行实践的改进已通过两种方法来寻求。
第一种方法包括得到常规拉冬变换上的采样要求(Schonewille和Duijndam，2001，Hugonnet和Canadas，1995，Marfurt，1996)以及通过数据采集、数据内插和预处理来满足这些要求(Manin和Spitz，1995)。
在第二种方法中，关于被混淆的数据部分的先验信息被包含在拉冬变换的计算中。在其最简单形式中，先验信息作为拉冬变换的最小二乘解中的对角线正则化项而引入。当沿着对角线的所有值都相等时，所得到的解一般不解决由于被混淆的输入数据而引起的多义性。Nichol(1992a，1992b)、Herrmann等人(2000)，以及Hugonnet(2001)已证明适当的对角线正则化产生准确而有效的拉冬变换，即使输入数据被混淆。
在贝叶斯估计的框架内，对角线正则化可以解释为把高斯先验分布指定给变换参数，这些参数也被假设为统计独立的(Tarantola，1987，Ulrych等人，2001)。其他算法通过避开高斯分布的假设，而需要更少的先验信息，但趋向于计算上更复杂(Sacchi和Ulrych，1995，Harlan等人，1984，Thorson和Claerbout，1985)。
关于高分辨率拉冬变换计算中的对角线正则化的现有技术包括由Nichols(1992a)，Herrmann等人(2000)，以及Hugonnet等人(2001)描述的方法。在这些情况中，为了计算效率，变换通常应用于具有频率相关对角线正则化的空间频率域。但是，正则化项并不作为所计算的拉冬变换的函数而更新。
Nichols(1992a)从沿着数据的瞬时频率轴的外观测量得到正则化权。非混淆地震能量的外观通常随频率缓慢变化。相反地，外观中的强而短规模的变化经常用来识别混淆能量。为了得到权，Nichols(1992a)公开在包含比正处理的频率低和高的频率的频率间隔上外观的平滑。
Herrmann(2000)和Hugonnet(2001)从最低频率处的常规拉冬变换开始，递归于频率而得到权。正处理的频率的权从较低频率处的高分辨率拉冬变换的结果建立。
图2，3和4中的例子说明用不同现有技术正则化方案计算的拉冬变换。图2a-2d对应于用标量乘以单位矩阵来正则化的常规拉冬变换。在图3a-3d中，正则化权是沿着对角线的可变项，并且相应的解具有较高分辨率。但是，在图4中，我们看到应用于被噪声污染的数据的同一正则化方案产生次最优结果。
图4e中的权不检测低和高频率处的抛物线事件。权可靠地指示仅在15Hz以上的数据中的时差(moveouts)，而在如图3e中所示无噪声情况下，时差从大约5Hz被检测。图4e中的权在变换域的边缘还具有虚假的高值，可能是由于从常规拉冬变换开始变换产生的人为信号(即等权)。这些高值在数据域中产生强的人为信号。
【发明内容】
一种用于得到高分辨率拉冬变换计算的正则化权的新方法被描述。
在本发明的一种实施方案中，一种处理地震数据的方法被提出。高分辨率拉冬变换被定义在地震数据上使用。高分辨率拉冬变换使用地震数据的外观测量来正则化。地震数据使用高分辨率拉东变换来处理，以增强地震数据的所希望特征。所处理的地震数据的有形表示被给出。
地震的外观测量可以包括沿着地震数据的一个维的外观测量。地震数据可以分成二维阵列，其包括时间或深度的一个维，以及空间表面位置或角度的另一个维。地震数据也可以分成多维阵列，其包含包括时间或深度的一个维，以及选自空间表面位置和角度的其他维。
地震数据的处理可以包括使用地震数据的外观测量，在地震数据上执行高分辨率拉冬变换。被拉冬变换的地震数据中第一区域和第二区域可以分离。在分离的、被拉冬变换的地震数据上的高分辨率拉冬变换的逆可以执行。
第一区域可以是信号区域，而第二区域可以是噪声区域。地震数据的维可以包括地震数据的多个维。地震数据的维也可以包括频域。
在另一种实施方案中，使用高分辨率拉冬变换以增强地震数据的所希望特征的地震数据处理还可以包括分段地逼近复杂时差轨道(moveouttrajectory)。地震数据可以分成与分段一致的局部数据窗口。高分辨率拉冬变换可以在每个局部数据窗口上执行，以增强与分段一致的地震数据的所希望特征。在每个局部数据窗口上执行高分辨率拉冬变换可以包括步骤：以周期边界条件计算变换，以及通过使用模型分量之间的已知的时移，施加与零值边界条件一致的信号噪声分离。
使用地震数据的外观测量来正则化高分辨率拉冬变换可以包括计算正则化矩阵，这包括将时差轨道所确定的相移施加到地震数据上。它也可以包括沿着时差轨道来归一化(normalize)叠加能量，以获得外观测量。它也可以包括在地震数据的第二维上削尖外观测量或平滑外观测量。地震数据的第二维可以包括选自空间位置和角度的维。
在另一种实施方案中，本发明可以在计算机系统上实现，其中计算机包括存储器和处理器，以及驻留于计算机存储器中的可执行软件。软件可用处理器运行，以定义供地震数据上使用的高分辨率拉冬变换。高分辨率拉冬变换使用地震数据的外观测量来正则化。地震数据使用高分辨率拉冬变换来处理，以增强地震数据的所希望特征。
软件也可以规定在地震数据上使用地震数据的外观测量的高分辨率拉冬变换的性能。第一区域和第二区域可以在被拉冬变换的地震数据中分离。在分离的、被拉冬变换的地震数据上的高分辨率拉冬变换逆可以使用。
软件也可以包括分段地逼近复杂的时差轨道。地震数据可以分成与分段一致的局部数据窗口。高分辨率拉冬变换可以在每个局部数据窗口上执行，以增强与该多个分段一致的地震数据的所希望特征。
所提出的方法可以保持与高分辨率拉冬变换关联的分辨率和去混淆品质，以及与Nichol和Herrmann方法关联的大部分效率，但是可以比Herrmann方法更稳定，并且提供比Nichol方法更好的分辨率。
鉴于下面的描述，本发明的其他特征和优点将更加明白。
【附图说明】
图1：(a)在最大偏移2000m包含具有0和400ms时差的抛物线事件的综合测试数据(b)加上随机噪声的相同数据集。
图2：施加到图1a的无噪声数据的常规拉冬多次衰减。(A)估计的主量(B)估计的多量(C)剩余量：输入数据同施加到输入数据的正和逆拉东变换的序列之差；以及(D)拉冬域。
图3：使用递归法，对图1a的无噪声数据施加高分辨率拉冬多次衰减，以得到权。图(a)到(d)与图2相同；最后一个图(e)：用于高分辨率拉冬变换的权。权显示为频率(垂直轴)和时差参数(水平轴)的函数。
图4与图3一样，但是对应图1b中所示的噪声综合数据。
图5与图3一样，但是使用外观得到的权。
图6：放大图4b(左边，权的递归更新)和5b(右边，基于外观的权)。
图7：与使用一个全局数据窗口和抛物线拉冬变换来处理双曲线事件相关联的聚焦不良和变换人为信号(a)剩余量。(b)在抛物线拉冬变换域中的数据。(c)用于高分辨率拉冬变换的权。
图8：说明从局部时间偏移窗口到抛物线拉冬变换域的映射的图。
图9：在三个分离的偏移窗口内，处理与图7中相同的输入数据的结果。(a)输入数据同施加到输入数据的正和逆拉冬变换的序列之差。(b，c，d)三个输入数据窗口的每一个的拉冬变换域中的数据。(e，f，g)用于三个输入数据窗口的每一个的高分辨率拉冬变换的权。
图10：现场数据实例，用常规拉冬变换来处理。(a)输入数据。(b)估计的主量。(c)估计的多量。(d)输入数据减去估计的主量。(e)正和逆变换的剩余量：输入数据减去估计的主量和多量之和。(f)抛物线拉冬变换域中的数据。
图11：现场数据实例，与图10中一样，其中权递归地计算。
图12：现场数据实例，与图10中一样，具有基于外观的权。
图13：权的比较：(13a，左边)基于外观的方法，以及(13b，右边)递归方法。
图14：拉冬变换方法的流程图。
图15：(a)准备输入到2D拉冬变换的数据区域，具有在每个区域内的局部数据窗口的示意指示。(b)拉冬变换的地震数据事件(曲线形状)和时差轨道(直线)的示意表示。(c)拉冬变换域的示意表示。在沿着直线求和之后，在变换域中的垂直线将被填入，并且大部分能量将贡献自与地震事件相切的时差轨道。
图16：处理地震数据的方法的流程图。
【具体实施方式】
本发明引入一种用于计算被正则化的高分辨率拉冬变换的新方法。
参考图16，地震数据使用行业标准方法来记录，并根据公共表面(例如，CMP集合)，或地面下位置(例如，在迁移过程中形成的角度集合)来分组1601。图15a显示通过用一个维作为时间或深度，并且另一个维作为空间位置(偏移)或角度而创建的2D阵列1502。时间或深度维通常密集地采样，从而不会混淆，然而由于沿着数据的空间位置(偏移)或角度维的不充足采样，混淆有可能发生。在其他实施方案中，不同类型的分组可能创建供处理的三维或更高维阵列。
为了变换数据，拉冬变换的时差轨道使用数据中事件的时差并考虑计算效率来选择1602。时间轨道的常见选择是线性的和抛物线的，如图15b和15c中所示。其他时间轨道，例如双曲线或更复杂的轨道也可以选择。
在拉冬变换过程中，数据在良好采样的维(时间或深度)上变换到频域1605。频率ω处拉冬变换的输入数据向量标记为d(x，ω)，简写为d。数据按通常表示为源和接收器之间的距离(偏移)的空间位置x来索引。数据向量d可以表示图15a中所示来自整个集合或来自局部(时间偏移)数据窗口的数据1501。与局部数据窗口中的处理相关的特定问题在下面名称为“在局部数据窗口中分离信号和噪声的准则”的章节中讨论。
频率ω处的拉冬变换的结果是向量m(q，ω)，或简写为m，也称作模型向量。m的元素按定义拉冬变换的所选时差轨道的参数q来索引。在抛物线拉冬变换的情况下，参数q解释为零偏移处的抛物线的曲率。通常，输入到变换的数据是2D区域。但是，在本发明的一些实施方案中使用的拉冬变换方法的描述可适用于更高维的数据，在这些情况下，曲率参数不再是标量，而是向量。
拉冬变换是将数据和模型向量相联系的线性变换。模型和数据向量之间的关系可以通过方程1(Beylkin，1987)用矩阵形式来表示：
d＝Am 方程(1)
其中，A是大小为(M，N)的矩阵，M是数据向量d中元素(复标量)的个数，而N是模型向量m中元素的个数。
拉冬变换可以定义为方程(1)中给出的线性方程系统的解。为了保证线性方程(1)的系统具有唯一且稳定的解，通常定义与方程(1)相关的目标函数J(m)，它是正则化项和数据拟合项的和：
J(m)＝mHW-1m+(d-Am)H(d-Am) 方程(2)其中，上标H表示向量或矩阵的共轭转置，而W-1是正则化矩阵。
W-1的最常见选择是沿着对角线具有恒定元素的对角线矩阵，例如W-1＝αI，其中I是N×N单位矩阵。正则化的这一选择与常规拉冬变换相关联，并被称作“阻尼最小二乘法”。它提供稳定的解，但是需要非混淆的大孔径输入数据。
在本发明的一种实施方案中，正则化权通过沿着偏移轴计算外观测量来得到1603：
S(p，ω)＝|AH(p，ω)d(ω)|α/(|d(ω)|α+ε(ω)) 方程(3)其中，AH是A的共轭转置矩阵；数据向量d(ω)可以包含来自整个集合或来自集合的局部时间偏移窗口的数据；α是用来削尖外观的参数；而ε是稳定化和/或归一化因子。矩阵A的共轭转置根据时差轨道来定义，并对数据施加适当的相移。
然后，外观测量可以使用标准加权平均方法(Marple，1987)，在频率范围上平滑，由此对角线权矩阵W计算为：
W(p，ω0)＝∑kαkS(p，ωk)/∑kαk 方程(4)其中，系数αk取决于正用于平滑的频率数，而不是取决于ω0，ω1等的实际频率值。系数αk的常见选择是恒定的、三角的和高斯型的权。外观测量在其上平滑的频率的范围优选地是宽频带，并且包括高频。本领域技术人员已知的其他类型的平滑函数可以用来得到基于外观的对角线权矩阵。
另外，类似于地震数据处理中的标准操作，外观测量可以通过在少量相邻CMP位置上空间平均化，或者通过滤去低于阈值的外观值来过滤，以提高稳定性。
一旦地震数据已使用高分辨率拉冬变换来变换，被变换数据的信号和噪声区域可以分离(也称作净噪)1605。分离这些区域的方法在本领域是已知的。最简单的方法是手动选择信号和噪声区域之间的划分几何(净噪函数)，并且将几何一侧的数据移走。如果地震数据分离成局部窗口，那么净噪可以通过在标题为“自局部数据窗口的变换的信号噪声分离(净噪)的有效实现”的章节中描述的方法来实现。
一旦噪声已被抑止，地震数据使用逆高分辨率拉冬变换来再变换1607。然后，地震数据可以表示为有形的表示，例如计算机打印输出或显示，或者在向用户显示或提交之前进一步处理1608。
在本发明的不同实施方案中，计算机或计算机系统中的软件可以用来实施本发明中给出的和本领域中已知的高分辨率拉冬变换。
抛物线拉冬变换在NMO修正之后，依赖于具有抛物线时差和恒定幅度的数据中的事件。但是，实际上，数据中的时差轨道是复杂的，不一定是抛物线的。在本发明的另一种实施方案中，复杂的时差轨道可以通过局部抛物线段的组合来很好地逼近。下面的讨论涉及抛物线，但是适用于通常用于拉冬变换的任何段。
局部数据窗口通过将数据的每个维分成覆盖该维的整个范围的段来定义。抛物线拉冬变换可以应用于每个局部窗口中的数据。当在每个局部窗口中应用抛物线拉冬变换时，抛物线时差和恒定幅度假设比在该维的整个范围上能更好地满足。不像常规拉冬变换，高分辨率抛物线拉冬变换能够把局部窗口中的信号和噪声数据分离。使用在局部窗口上很好地逼近复杂轨道的抛物线段，信号和噪声分离可以根据复杂时差轨道来有效地表达。通过模拟复杂的时差轨道，选择数据中所希望特征的准则可以更加有效。
图7说明具有不是很好聚焦给在整个偏移范围上计算的抛物线拉冬变换的剩余双曲线时差的地震事件。当将偏移范围分裂成三个部分并且在每个部分上分别计算抛物线拉冬变换时，每个变换被很好地聚焦，如图9中所示。
但是注意到，图9中偏移范围相关变换的每个将双曲线聚焦于模型空间中稍微不同的位置。在每个偏移窗口中的最佳拟合抛物线具有稍微不同的参数(曲率和零偏移时间)。这一观测意味用来分离变换域中的事件的静噪函数也作为用于计算变换的数据的偏移范围的函数而变化。
为了管理静噪参数随偏移的变化，我们引入下面的方法。首先，以包含零偏移的窗口的静噪函数q(τ，x＝0)开始。在给定零偏移时间τ，由q0表示的相应的静噪位置q0＝q(τ，x＝0)。参数(τ，q0)定义抛物线t＝τ+q0x2/xmax2，以及具有与抛物线相同的零偏移时间和零偏移曲率的双曲线：
θ2＝τ2+2q0τx2/xmax2 方程(4)
那么，偏移窗口之间的静噪关系为：
qref＝q0*τ/θ(τ，xref) 方程(5)其中，qref是与偏移xref关联并对应于偏移时间τ的数据窗口的期望静噪位置。时差参数qref通过在位置xref处匹配抛物线(τ，qref)和双曲线(τ，q0)的斜率来获得。
一般地说，上述过程定义与全局信号噪声分离准则一致的并且根据复杂的(例如双曲线)静噪函数来表示的局部信号噪声分离准则。关于现有技术的新特征是使用全局信号噪声分离准则的时差轨道(例如，双曲线)，它并不是用于拉冬变换计算的时差轨道(例如，抛物线)中的一部分。这一新特征为信号噪声分离增添适应性和改进的准确性，同时了保持计算效率。实际上，改进的信号噪声分离可以改善地震数据中所希望特征的分离。
在本发明的另一种实施方案中，数据域(x-t)中和拉冬变换域(q，τ)中的窗口之间的最佳时移可以被计算，如图8中所示。
考虑由时间间隔[t1，t2]和偏移间隔[x1，x2]定义的时间偏移窗口，并且假设xref是窗口的适当参考偏移(典型地中心偏移)。对于给定的q值，跨越窗口内的参考偏移的抛物线的零偏移时间处于间隔[τ1，τ2]中，其中τ1＝t1-Δτ(q)，并且τ2＝t2-Δτ(q)，其中Δτ(q)＝qxref2/xmax2。可以看到，τ中的最佳窗口比t中的要早时移Δτ(q)(对于正的q)。注意到，时移在q中是线性的，因此t-x中的矩形窗口在τ-q中由平行四边形来表示。当q和xref很大，时移Δτ可以明显地比时域中的窗口长度要长。
因为频域中的处理意味着数据和模型两者都是周期性的，具有周期t2-t1，模型总是有效地定义，因此我们仅从该模型中选择适当的(平行四边形)窗口，已知时移Δτ(q)。实际上，这里所描述的过程允许使用周期边界条件来有效地计算变换，当实施信号噪声分离时，就好像边界条件是在输入数据间隔[t1，t2]外的零值(优选描述)。
为了说明本发明的实施方案，图1b的合成(加上随机噪声)已被显示。对于具有附加噪声的数据的情况下，结果显示在图5a到5e中。基于外观的权在70Hz到90Hz的频率范围中得到，并假设在整个频率上是恒定的。由于更准确的权，主量和多量的分离现在被改善了。在变换域的边缘处没有高值。
图6是图4b和5b中的结果的放大比例显示。它清楚地说明当使用基于外观的权(图6b)而不是递归更新的权(图6a)时，主量和多量改善的分离。
该实例的现场数据是CMP集合(图10a)，具有100m轨迹间距和大约6km的最大偏移。常规拉冬变换的结果(图10b到10f)说明遇到的常见问题：变换域中差的聚焦，变换域边缘处显著的能量，变换中和数据域中混叠的人为信号。
本发明的结果在图12b到12f中显示。相应的基于外观的权在图13a中显示。正如所料，与常规拉冬变换的结果(图10b)相比，混淆和变换人为信号减少，并且估计的主事件(图12b)被改善。
使用权的递归更新的高分辨率拉冬变换的结果在图11a到11f中显示，相应的权在图13b中显示。这些结果并不像用基于外观的权所获得结果一样好。这在图11f(变换域)和11b(估计主量)中尤其明显，其中在变换域边缘处的一些混淆能量和人为信号现在再次可见，并可能与由递归方法(图13b)估计的次最佳权有关。
本发明的许多实施方案已经描述。然而，应当明白，可以不背离本发明的本质和范围而做各种修改。因此，其他实施方案处于下面的权利要求的范围内。