图1是海洋地震测量的示意图，其中地震能量由震源1发出并被海 面6之下深度h处的地震接收器2检测到。由震源1发出的能量被海底3 和海底3下的反射体4反射，然后被接收器检测到。在图1中这一地震 能量路径标为5。由入射到接收器上的反射地震能量能导出关于地球内部 地质结构的信息。
图1中所示地震接收器2是一个海上拖缆，它是在海洋地震测量中 经常使用的一种地震接收器。一个海上拖缆含有多个传感器S1、S2、…SN 沿其长度分布，如压强传感器和/或质量速度传感器，海上拖缆的长度可 以是几百米，从而能同时在多个点测量被反射的地震能量。一个海上拖 缆可由一个或多个浮子8悬挂，从而使海上拖缆的所有接收器在平静海 中的同一深度。
除了图1中所示所希望的地震能量路径5之外，还将发生其他地震 能量路径，这是地震能量从海面6反射或散射的结果。这些路径被称作“虚 反射”。例如，图1中的参考数字7显示一个虚反射，其中被反射体4反 射的地震能量不是直接入射到接收器2上，而是在到达接收器前在海面6 受到进一步反射。下行海面虚反射是不希望的地震数据污染源，因为它 们模糊对来自地球内部的希望的上行反射的解释。如众所周知的那样， 虚反射将在地震能量频谱中产生一个或多个“凹槽”，这些凹槽所在频率 依赖于接收器在海面下的深度。
从海面反射的虚信号相对于直达信号有一个延时。造成这一延时的 有两个分量：第一，在海面反射时有一个180°相位变化，第二，存在一 个对应于附加路径长度(对于垂直方向发射的能量，这是2Z，或者说接 收器深度的二倍)的时间延迟。实际的垂直远场信号是直达信号和虚信 号之和。直达信号和虚信号将发生干涉，这造成远场信号振幅的变化。 对于某些频率，这种干涉是破坏性的，在频谱中造成零振幅或“凹槽”。 这发生在其接收器深度是四分之一波长的偶数倍的频率上：
fnotch＝(nc/2z)，n＝0，1，2，3…
(这里C是水中的声速，n是整数，给出谐振阶数，Z是接收器在水面 下的深度)。
建设性干涉刚好发生在相邻凹槽频率的中间位置。这导致在这些频 率上有振幅极大值，由下式给出：
fpeak＝(2n+1)c/4z，n＝0，1，2，3…
直达信号和虚信号之间的干涉效应可被认为是对直达信号应用了一个频 率域“虚滤波器”。这一虚滤波器有如下形式：
g(f)＝1+|r2|-2|r|cos(2zf/c)
(这里r是在海面的反射系数)。该虚滤波器的振幅谱类似于一个被全波 整流的正弦波，其零点在虚凹槽频率，而振幅峰值2.0(6dB)在那些峰 值频率上。
图2显示一个典型的虚滤波器振幅作为频率的函数。该图显示 Z＝12m、c＝1500m/s而且在海面的反射系数为-1.0时的虚滤波器。可以看 到，在凹槽频率0Hz、62.5Hz、125Hz…，振幅降为0，而在峰值频率 31.25Hz、93.75Hz…有振幅极大值。
地震消虚(seismic deghosting)是一个古老的和长期存在的问题。 传统上，消除虚反射的主要目标一直是加宽通过凹槽频率的数据的带宽。 对地震数据消除虚反射的许多先有途径所使用的假定是完全平坦的海 面、在海面下已知深度的海上拖缆以及在海面处压强消失。然而，在实 际的地震数据采集中，海面经常是崎岖不平的。在消虚计算中使用一个 显式平均拖缆深度的确定性算法在强浪条件下将不能可靠地工作。即使 是统计方法也不能在强浪条件下完全纠正虚反射，这是由于该问题的时 变性质，如E.Kragh和R.Laws在“强浪和统计反褶积”(EAGE第62届 年会(2001))一文中说明的那样。
强浪对地震数据的影响已是近年来相当多研究的主题。特别是 R.Laws和E.Kragh已在“时间延迟地震和强浪小波”(勘探地球物理学会 第70届国际会议扩展摘要第1603-1606页(2000))一文中已说明，对强 浪条件下采集的数据的处理过程中所产生的误差对于时移地震测量和对 于为地层反演而可靠采集可重复数据的影响是显著的。
如在英国专利申请0015810.5中指出的那样，对强浪海面虚反射的纠 正需要波场被完全地分离成其上行和下行分量。通常这种分离需要记录 压强和质点运动速度(或等效的垂直压强梯度)二者(例如，见 L.Amundsen在Geophysics(地球物理学)第58卷第1335-1348(1993) 中给出的公式)。
如J.O.A.Robertsson和E.Kragh在“使用单个海上拖缆和压强梯度近 似在强浪条件下消除虚反射”(提交给Geophysics(2001))一文中建议的 并由Rsten等在“使用垂直质点速度场近似在强浪条件下消除虚反射”(提 交EAGE第63届年会)一文中进一步研究的那样，另一种可能性是利用 压强在海面消失这一事实从在接收器采集的压强导出质点速度垂直分量 的估计。沿垂直位于海上拖缆上方一个等值线的零压强可被看作是在双 拖缆配置中的第二个拖缆。
Robertsson和Kragh(上文)以及UK专利申请0015810.5建议使 用来自地震波传播的会成有限差分模型的一种技术，也称作 Lax-Wendroff校正，来导出强浪表面附近拖缆处垂直压强梯度的方程式。 这些方程式利用沿拖缆的压强梯度和压强的时间导数表示垂直压强梯 度，而沿拖缆的压强梯度和压强的时间导数二者可容易地从拖缆采集的 数据中得到。然后，可使用压强和对垂直压强梯度的近似来对地震数据 消除虚反射。与地震消虚的其他建议相反，Robertsson和Kragh的方法 是局部的：一个短空间滤波器(通常为长度三)被用于对单个接收器获 取的压强数据消除虚反射。

本发明的第一方面提供一种方法，用于由水柱内放置的接收器所采 集的地震数据确定质点速度的垂直分量，该方法包含：使用一个算子由 接收器采集的压强确定质点速度的垂直分量，对于在水柱表面下直至达 到至少是为采集地震数据所用地震能量最小波长0.4倍的深度所采集的 地震数据，该算子是准确的。在一个优选实施例中，对于在水柱表面下 基本上是任何深度所采集的地震数据，该算子都是准确的。
本发明提供了估计垂直质点速度的一种改进的方法。估计质点速度 垂直分量的先前方法只适用于水柱表面下接收器深度的一个受限制范 围。如下文讨论的那样，先有技术的方法只当接收器在水柱表面下的深 度不大于λ/3.5时才是准确的，这里λ是采集过程中所用地震能量的最小 波长。这严重地限制了对能使用该方法的接收器阵列的选择。本发明提 供了估计质点速度垂直分量的方法，它能用于在比以往都大的深度采集 的数据。
该算子可以是空间的和/或依赖波数的算子。
在优选实施例中，该方法包含：使用下列方程或由此导出的方程由 接收器采集的压强确定质点速度的垂直分量：
$v_2(\omega ,x,y,z)=-\frac{i}{\mathrm{\rho \omega z}}\underset{m}{\Sigma }{F}_m(k,z)D_2^m(x,y)*p(\omega ,x,y,z)---\left(2\right)$
这里p是采集的压强波场，Vz是质点速度的垂直分量，ω是角频率，k＝ (ω/α)是波速，α是P波速度，ρ是密度，D2是差分算子， $\mathrm{Fm}=\mathrm{Fm}(k,z)=k^{-2m}\underset{n=m}{\overset{\infty }{\Sigma }}\mathrm{Anm}{\left(\mathrm{kz}\right)}^{2n}$以及*是2D(二维)空间褶积运算符。
方程式(2)是一个严格表达式，对任何接收器深度都是准确的。使 用方程式(2)使得能由在任何接收器深度采集的地震数据得到垂直质点 速度。
为减少应用式(2)所需计算，求解由式(2)导出的一个方程是可 能的，其中求和被截断，而不是严格地求解式(2)。解方程式(2)的截 断形式减少了所需计算，但其代价是准确性。式(2)的被截断版本通常 不是对任何接收器深度都准确，而是只对小于某一极大值的接收器深度 才准确。
在另一个优选实施例中，该方法包含使用下列方程或由此导出的方 程由接收器采集的压强确定质点速度的垂直分量：
    vz(ω，x，z)≈Fk*p(ω，x，z)
${\tilde{F}}_k=\underset{b_m,a_n}{\min }\underset{k=-(K-1)}{\overset{K-1}{\Sigma }}{W}_k{|\frac{{\Theta }_{\mathrm{km}}{b}_m}{{\Phi }_{\mathrm{kn}}{a}_n}-F_k|}^2$
这里an和bm是后向和前向滤波器系数，Wk是一组正加权，*是1-D(一 维)空间褶积运算符，Fk＝F(kΔkx)代表所希望的离散水平波数响应。
这一滤波器是下文附录中给出的式(A34)和(A36)的一个近似。 方程式(A34)和(A36)是严格的，但对强浪表面的情况无效，这一近 似允许它们应用于强浪表面的情况。
在一个优选实施例中使用如下滤波器：
${\tilde{F}}_k=\underset{g_m}{\min }\underset{k=0}{\overset{\kappa -1}{\Sigma }}{W}_k{|{\Gamma }_{\mathrm{km}}{g}_m-F_k|}^2$
这里
Γkm＝cos(kΔkxmΔx)，bm＝0.5gm for gm＝1，2，…，M/2 and bo＝go
本发明提供的对垂直质点速度的估计可用于对数据消除虚反射。另 一种作法是，由本发明提供的对质点垂直速度的估计可用于其他地震数 据处理应用中，如估计震源特性(Signature)。
本发明的第二方面提供一种方法，用于处理由水柱内放置的海上拖 缆上的接收器所采集的地震数据，该方法包含由沿着拖缆方向的压强导 数确定沿垂直于拖缆的水平方向的压强导数的步骤。这可由下式完成：
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{1}{x_0}\frac{\partial }{\partial x}p(x_0,y_0)$
这里拖缆沿x轴延伸，震源位于(O，y0)。
在该式中，接收器在x-y平面中的坐标是(x0，y0)。这使得即使对单 个拖缆也能估计出在交叉线方向的压强导数。相反，先有技术的方法至 少两个横向分离的拖缆来确定沿交叉线方向的导数。
本发明的第三方面提供一个装置，用于由水柱内放置的接收器所采 集的地震数据确定质点速度的垂直分量，该装置包含使用一个算子由接 收器采集的压强确定质点速度垂直分量的手段，对于在水柱表面下直至 达到至少是为采集地震数据所用地震能量最小波长0.4倍的深度所采集 的地震数据，该算子是准确的。
本发明的第四方面提供一种装置，用于处理由水柱内放置的海上拖 缆上的接收器所采集的地震数据，该装置包含由沿着拖缆方向的压强导 数确定沿垂直于拖缆的水平方向的压强导数的手段。该装置可根据下式 确定沿垂直方向的导数：
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{1}{x_0}\frac{\partial }{\partial x}p(x_0,y_0)$
这里拖缆沿x轴延伸。
在一个优选实施例中，该装置包含一个可编程数据处理器。
本发明的第五方面提供一个存储介质，其中含有上文定义的装置中 的数据处理器所使用的程序。
本发明的又一方面提供一种确定质点速度垂直分量的方法，包含使 用下列方程或由此导出的方程由接收器采集的压强确定质点速度的垂直 分量：
$v_z(\omega ,x,y,z)=-\frac{i}{\mathrm{\rho \omega z}}\underset{m}{\Sigma }{F}_m(k,z)D_2^m(x,y)*p(\omega ,x,y,z)$
本发明的优选特性在独立权利要求中定义。

现在将参考附图以举例方式描述本发明的优选实施例，附图中：
图1是拖曳海上地震测量的示意剖面图；
图2显示虚反射对图1所示测量中由震源发射的震能量的振幅谱的 影响；
图3是本发明方法的方块流程图；
图4是图1的地震测量安排的示意平面图；
图5(a)至5(d)显示对于深度6m的海上拖缆本发明的消除虚反 射方法的结果；
图6(a)至6(d)显示对于深度10m的海上拖缆本发明的消除虚 反射方法的结果；以及
图7是根据本发明的装置的示意方块图。

从水柱内放置的地震接收器采集的压强数据中去掉虚反射的影响等 效于把所采集的压强分离成它的上行和下行分量。所采集的压强的消除 虚反射后的分量是上行分量。已经提出了各种方法把所采集的压强数据 分离成上行和下行分量，包括下式：
$\tilde{P}=\frac{1}{2}[P-\frac{\mathrm{\rho \omega }}{k_z}{V}_z]$
在这一方程中， 代表所采集的压强消除虚反射后的(上行)分量的空间 付立叶变换，P代表所采集的压强p的空间付立叶变换，ρ代表水柱的密 度，ω代表角频率，kz是垂直波数，Vz是质点速度垂直分量vz的空间付 立叶变换。
Vz、P和 一般是水平波数kx和ky、角频率以及接收器深度z的函 数。上述方程式(1)对于所有波数都是有效的。
对水柱内放置的地震接收器所采集的压强数据消除虚反射的一种先 有技术方法是使用方程式(1)确定消除虚反射后的(上行)压强分量。 这一方法直到现在仍存在的缺点是尚未得到vz的足够准确的表达式。尽 管已提议计算vz的各种近似公式，但通常这些公式只可应用于有限范围 的接收器深度或地震能量频率。例如，vz的一个先有技术近似是在附录 的式(A1)和(A2)中给出的近似，这一先有技术的近似只当接收器深 度小于约λ/3.5时才是准确的，这里λ是地震能量的最小波长。这一近似 还要求一个拖缆上相邻接收器之间的间距约3m或更小。再有，这一先有 方法基本上局限于地震能量频率低于第一虚凹槽频率的情况，因此，例 如对于拖缆深度6m的情况，方程式(A1)只在达到约70Hz之前是准确 的。
本发明提供垂直质点速度Vz的更准确表达式，从而能使用式(1) 对压强数据更准确地消除虚反射。
根据本发明，所要求的质点速度的垂直分量由下式确定，或由从中 导出的方程式确定：

在式(2)中，函数Fm(k，z)是不依赖于水平空间坐标的函数，由下式给 出：
$F_m=F_m(k,z)=k^{-2m}\underset{n=m}{\overset{\infty }{\Sigma }}{A}_{\mathrm{nm}}{\left(\mathrm{kz}\right)}^{2n}---\left(3\right)$
这里
$A_{\mathrm{nm}}=\frac{{(-1)}^n{2}^{2n}{B}_{2n}}{\left(2n\right)!}\left(\begin{array}{c}n\\ m\end{array}\right)---\left(4\right)$
在式(4)中，Anm是标量，Bn是第n个Bernoulli(伯努利)数，n！ 是n的阶乘。
在前面的式(2)中，D2是差分算子，*代表二维空间褶积运算符。 分量P、vz一般是ω、x、y和z的函数。
差分一般导致对高波数的大放大。因此，在实际情况中，优选地对 D2使用一个差分算子的限定带宽版本。
式(2)中给出的vz表达式对于水柱中所有接收器深度都有效。使用 式(2)确定质点速度的垂直分量使能对水柱中任何深度采集的压强数据 准确确定质点速度的垂直分量。这又使得能对水柱中任何深度采集的压 强数据准确地进行那些需要知道质点速度垂直分量的进一步处理步骤， 如消除虚反射。
本发明不限于使用严格的方程式(2)来得到质点速度的垂直分量。 本发明还考虑可使用由式(2)导出的方程式，如式(2)的近似解，来 得到质点速度的垂直分量。
当m为有限值且从零延伸到M时，式(2)产生一个特例。在这种 情况中，并假定无限带宽，则式(2)可重写为：

函数Fm在上文的式(3)中给出，为无限级数。函数Fm也可以分析上 利用三角函数给出，F0、F1和F2的举例下文的式(A8)、(A9)和(A10) 中给出。
将会看到，在式(6)中的x和y变量是与z不耦合的。这表明式(6) 不仅对水体的平表面有效，而且对水体的崎岖表面有效。所以，本发明 的这个优选实施例提供了进一步的好处，因为它使得能对水柱中的任何 接收器深度对强浪条件下采集的压强数据进行准确的消除虚反射。
在M＝2的情况中，式(6)减化为：
$v_z\approx -\frac{i}{\mathrm{\rho \omega z}}[F_0+F_1(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})+F_2{(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})}^2]p---\left(7\right)$
可由两个或更多个接收器采集的压强来估计式(6)和(7)中的空间导 数。例如，可由同一时刻沿拖缆的不同接收器记录的压强来确定沿x轴 (假定它与拖缆轴重合)的压强导数。
确定空间导数的优选方法将依赖于空间采样间隔，即一个拖缆上相 邻传感器之间的距离。在低采样间隔时，即传感器密集放置因而提供密 集采样的压强数据时，空间导数可由经典的中心有限差分近似来代表。 例如，在式(7)的情况中，可以使用三点有限差分算子来估计二阶空间 导数，该算子由一个接收器处以及在每一侧的相邻接收器采集的压强值 估计在那一点沿x方向的空间导数(这有时称作长度为3的迭代滤波器)。
对于稀疏的接收器放置，水平导数最优选地实现为数值优化差分算 子，或使用付立叶技术，在付立叶技术中，对压强进行付立叶变换以变 换到波数域，乘以波数平方的负值(对于二阶导数)，再对结果进行逆付 立叶变换，从而得到导数。
式(6)和(7)是从式(2)导出的并且是式(2)的特例，其中变 量m被选为有限值。当式(2)中的变量m和n二者都选为有限值时则 发生又一些特例，其中m从0变到M，n从m变到N(m)。在这种情 况中，式(2)简化为：

在式(8)中，函数Fm有与式(3)中给出的相同的一般表达式，但对n 的求和现在是从m进行到N(m)，如下文中式(A20)所示。
如果我们选择M＝1，N(0)＝2和N(1)＝1，式(8)简化为下式：

这里标量系数有如下值：A00＝1，A10＝A11＝-1/3。
将会看到，式(9)与下文附录中给出的先有技术方程式(A1)有相 同的一般结构。然而，式(9)中的标量系数A10和A11不同于式(A1) 中的标量系数。这有显著影响，即式(9)对于比式(A1)更大的拖缆深 度有效。
(上文的式(2)对于所有拖缆深度有效，但对式(2)的近似解的 有效深度范围将因不同的近似解而不同。)式(9)在深达约λ/2.5时仍有 效，这里λ是最小波长，而式(A1)要求拖缆深度不大于约λ/3.5。
图3是说明本发明的一个方法的方块流程图。
开始时，在步骤9，地震数据在放在水柱内的多个接收器处被采集。 通常，步骤9将包含使用图1所通用型地震测量安排采集地震数据，其 中多个地震接收器S1…SN分成在一个海上拖缆2上，该拖缆放置在水柱 表面下深度h处。这些数据可由单一拖缆采集，或者它们可由横向彼此 分开(即沿y轴分开)的两个或更多个拖缆采集。
本发明还可用于处理预先存在的地震数据。所以，采集地震数据的 步骤9可被从存储器中检索预先存在的地震数据的步骤10取代。
在步骤11，对一个接收器(例如第j个接收器)计算质点速度的垂 直分量Vz。使用式(2)或从式(2)导出的任何方程式，如方程式(6)、 (7)、(8)和(9)，进行步骤11。如前文概括的那样，在这些方程式中 出现的压强导数值是使用由两个或更多个接收器采集的压强值估计出来 的。
一旦对第j个接收器确定了质点速度的垂直分量，则可在步骤12对 第j个接收器采集的压强数据进行消除虚反射。消除虚反射是使用式(1) 进行的，分量Vz是对步骤11的输出进行空间付立叶变换得出的。
在步骤13，步骤12的结果被输出，例如由操作员显示，或存储起来 供其后检索。
图1的拖缆2上的接收器Si通过在离散时间值对入射地震波场采样 来采集数据，于是将记录波场在时刻t0、t0+ts、t0+2ts…的值，这里ts是 采样时间。步骤11至13已在一个特定采样时间采集的数据上进行。所 以，在步骤14，确定是否应对该接收器在另一个采样时间采集的数据重 复消除虚反射过程。如果确定为“是”，则在步骤15选择新的采样时间， 并对以此新的时间采集的数据重复步骤11-13。这一过程重复至步骤14 得到的确定为“否”时为止。
在步骤16，确定是否应对拖缆上的另一个接收器采集的数据重复步 骤11、12和13。如果步骤16产生的确定为“是”，则在步骤17选择一个 新的接收器，并对新的接收器重复步骤11至14，必要时还有步骤15。 这一循环重复至步骤16得到的确定为“否”时为止。
如果希望对拖缆上每个接收器采集的数据系统地消除虚反射，则标 号j在开始时可设为1，步骤16可构成为确定标号j是否等于该拖缆上 的接收器总数J，而步骤17可构成为对j增量1。
当在步骤16得到的确定为“否”时，步骤18确定是否应对另一个拖 缆采集的数据重复该过程。如果这一步骤产生的确定为“是”，则在步骤 19选择一个新的拖缆，并对新选择的拖缆上的接收器重复步骤11至17。
当本发明方法应用于只使用单个拖缆采集的数据时，步骤18可被略 去，而步骤19是不适用的。
如前文指出的那样，本发明可应用于只使用单个拖缆采集的地震数 据，即应用于全在相同y坐标值处的接收器采集的地震数据。在这种情 况中，在初始时看来将难于得到沿y方向(垂直于拖缆长度)的空间导 数，这些是允许求解例如式(6)、(7)、(8)或(9)所需要的。所以， 本发明的又一方面是提供一种方法，用于由沿拖缆的压强导数估计必须 的沿y方向的导数。一个拖缆通常含有大量压强传感器以规则间隔沿拖 缆的长度分布，由这些传感器采集的压强能容易地估计沿拖缆的压强空 间导数。本发明使得有可能由单个拖缆采集的数据，由沿拖缆的压强空 间导数，估计出沿垂直于拖缆的方向的空间导数。
由单个拖缆采集的数据估计沿垂直于拖缆的方向的空间导数的方法 可以结合根据上述方法确定质点速度垂直分量的过程一起进行，但不限 于结合确定质点速度垂直分量的过程一起使用。
在一个优选实施例中，可以使用如下议程式估计沿垂直于拖缆方向 的空间导数：
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{1}{x_0}\frac{\partial }{\partial x}p(x_0,y_0)---\left(10\right)$
在下文的附录中详细解释方程式(10)的导出。
图5(a)至5(d)和6(a)至6(d)说明由本发明提供的结果。 这些结果是通过模拟强浪表面合成地震数据并使用本发明的方法对合成 地震数据消除虚反射而得到的。合成地震数据是使用由R.Law和 E.Kragh(上文)以及由E.Kragh和R.Laws(上文)发展的方法计算出 来的。概括地说，产生强浪表面，并计算出这些表面下的脉冲响应。这 些脉冲每三个构成一组，总共计算出40组，每组是由4m显著浪高海面 的不同实现计算出的。数据是对深度为6.0m的拖缆(图5(a)至5(d)) 及深度为10.0m的拖缆(图6(a)至6(d))合成的。
然后使用式(7)对合成地震数据消除虚反射以计算垂直质点速度。 在使用式(7)的计算中，参数设计如下：M＝1，N(0)＝∞，N(1)＝2， 于是给出F0＝kzcot(kz)， $F_1=-\frac{z^2}{3}(1+\frac{2}{15}{k}^2{z}^2),$F2＝0。
最大频率为50Hz的零相位Ricker小波与这120个脉冲褶积以估计 vz分量。在估计vz分量时，公知的长度3中心差分算子被用于得到沿x 方向的二阶水平空间导数。(该计算是使用2D(二维)估计进行的，因 此在交叉线方向(y方向)的导数变为零并可忽略。)
一旦使用式(7)估计出vz分量，使用上述式(1)估计消除虚反射 后的压强。在确定消除虚反射的压强之前，模拟的压强和估计的vz分量 被以因子4在时间上降低采样率，从而使时间上的采样距离从0.25ms增 至1.0ms。
图5(a)显示原始合成地震数据，总体被指为A，使用Robertsson 和Kragh(上文)的Lax-Wendroff技术对合成数据消除虚反射的结果被 指为B。
图5(b)显示原始合成地震数据，总体被指为A，使用式(7)对合 成数据消除虚反射的结果被指为C。图5(c)显示振幅谱扩展中的95％ 置信水平，而图5(d)显示相位谱扩展中的95％置信水平。这些结果是 对指为A的原始合成数据、指为B的用先有技术的Lax-Wendroff技术 对合成数据消除虚反射的结果以及指为C的使用式(7)对合成数据消除 虚反射的结果显示的。尽管没有从本发明得到的消除虚反射后数据中完 全去掉强浪扰动，但校正显示出在第一虚凹槽频率(约125Hz)以下的 明显改善。剩余误差在振幅中为0.5-1.0dB量级，在相位中约为5°，这些 误差显著低于当前对实际地震数据消除虚反射中得到的误差。
图6(a)至6(d)总体上分别对应于图5(a)至5(d)，只是它们 涉及10.0m拖缆深度。对于处在这一深度的拖缆。第一虚凹槽频率约为 75Hz。
截止频率为10-15-90-100Hz的一个零相位带通滤波器应用于图2中 所示小波，而截止频率为10-15-65-70的一个零相位带通滤波器应用于图 6(a)至6(d)所示小波。这对上文提到的也是下文中结合式(A14) 和(A15)的有限频带差分有类似的作用。这里，该滤波器在预处理步骤 直接应用于数据，这有类似的总体效果。
由本发明的方法得到的消除虚反射后的数据可受到进一步处理。例 如，消除虚反射后的数据可用于估计地震震源1的特征。例如，这可根 据L.Amundsen等在Geophysics第60卷第212-222页(1995)中提议的 方法来完成。
或者，消除虚反射后的数据可被进一步处理以去掉或至少是衰减数 据中由路径产生的事件，这些事件涉及在水柱自由表面的反射和/或散射。 这些事件通常称作与自由表面关联的多重波，因为产生它们的路径多次 穿过水柱。在一个实施例中。去掉与自由表面关联的多重波所需主要步 骤如下：
首先，拖缆采集的入射压强波场中的直接到达和与直接到达关联的 虚反射到达被隔离和从压强波场中去掉。这可以以任何适当的方式完成， 例如，通过“静噪(mute)”，直接到达及其相关联的虚反射到达被简单 地与其他到达分离。假定震源到接收器的炮检距不是太大而且水深不是 太浅，则地震数据中的直接到达及其相关联的虚反射到达通常能很好地 与其他到达分开。
然后，根据上文描述的方法估计质点速度的垂直分量。然后使用质 点速度的垂直分量将剩余的入射压强波场(即已经去掉直接到达及其关 联的虚反射到达之后剩余的压强波场)分开成它的上行和下行分量，例 如使用上文中的方程式(1)，以确定压强波场的上行分量。
一旦找出压强波场剩余部分的下行分量，它被加到原来的直接到达 及其关联的虚反射到达中。这给出总的下行波场。(在一个典型的测量 中，震源将比接收器更接近海面，因此直接到达及其关联的虚反射只含 有下行能量。)然后，可使用L.Amundsen在“在无需震源小波的情况下 消除与自由表面关联的多重波”(Geophysics，第66卷，第327-341页 (2001))一文中提议的“U除以D”方法。去掉由于自由表面多重反射造 成的影响。在这一方法中，“U除以D”是指去掉直接到达及其关联的虚 反射到达之后剩余波场的上行分量除以波场的整个下行分量。这一除法 实质上是反褶积步骤。
上文描述的本发明的实施例可在水柱内放置的一个或多个压强传感 器采集的地震数据上进行，例如在一个海上拖缆采集的压强数据上进行。 现在将描述本发明的另一些实施例，其中的地震数据包括压强数据和质 点速度二者。这类数据可由具有压强传感器又有质点速度传感器的地震 接收器采集，例如由海上拖缆采集。具有压强传感器又有质点速度传感 器的海上拖缆能直接测量流体压强P和质点速度V。
量p和v不是彼此独立的，而是由下列各式关联：
$\stackrel{·}{p}=-K▿.v---\left(11\right)$
$v=-\frac{1}{\rho }▿p---\left(12\right)$
式(11)中， 代表p的时间导数，K代表水柱的体积模量。
在本发明的一个实施例中。假定拖缆深度已知。例如，如果向拖缆 提供深度传感器，便可知道深度。如果拖缆深度已知，则可能将质点速 度传感器的输出与流体压强传感器的输出进行比较。例如，根据式(2) 或式(2)的导出方程，由压强的测量值和拖缆深度计算出质点速度，再 将计算出的质点速度值与质点速度直接测量值进行比较，从而完成这一 比较。当然，如果质点速度和流体压强测量准确的话，这些值应该是相 同的。这一实施例可被认为是相对于质点速度传感器来校准压强传感器， 并允许监视这些传感器的准确度。
在由本发明的方法处理所采集的地震数据的同时，可进行监视传感 器准确度的这一过程。或者，它可以单独进行，例如在地震测量开始时 检验传感器。这允许由测量用于校准传感器的p和v二者提供额外的自 由度。
在这一实施例中，例如通过对从浮子悬挂拖缆的结构安排的了解， 能估计拖缆的深度。然而，为改善准确度，最好是直接测量接收器的深 度。
本发明还可用于确定接收器的深度。为此，使用上文描述的方法确 定质点速度的垂直分量。可由质点速度垂直分量的测量值和质点速度垂 直分量的估计值估计出接收器深度。
在本发明的另一实施例名，使用下列方程式或其导出方程式，由接 收器采集的压强确定质点速度的垂直分量：
vz(ω，x，z)≈Fk*p(ω，x，z)，                       (13)
这里
${\tilde{F}}_k=\underset{b_m,{\alpha }_n}{\min }\underset{k=-(K-1)}{\overset{K-1}{\Sigma }}{W}_k{|\frac{{\Theta }_{\mathrm{km}}{b}_m}{{\Phi }_{\mathrm{kn}}{a}_n}-F_k|}^2---\left(14\right)$
在式(14)中，an和bm代表后向和前向滤波器系数，Wk是一组正 权重，Fk＝F(kΔkx)代表所希望的离散水平波数响应。方程式(14)的 导出在下文的附录中解释。
如前文指出的那样，方程式(14)是下文的附录中给出的式(A34) 和(A36)的近似。对任何深度都存在式(14)的解。式(14)可应用于 在强浪表面情况下得到的数据，而对于式(A34)和(A36)是不可能的 (尽管它们对于平坦海面是严格的方程式)。如附录中解释的那样，方程 式(14)可被实现以给出如下的零相位非递归滤波器：
${\tilde{F}}_k=\underset{g_m}{\min }\underset{k=0}{\overset{\kappa -1}{\Sigma }}{W}_k{|{\Gamma }_{\mathrm{km}}{g}_m-F_k|}^2---\left(15\right)$
其中，Γkm＝cos(kΔkxmΔx)，bm＝0.5gm for gm＝1，2，…，M/2 and bo＝go。
使用式(13)和(14)或使用式(13)和(15)得到的垂直质点速 度可用于这里描述的任何应用。
式(13)、(14)及(15)是1-D(一维)方程式，但也可使用相应的 2-D(二维)方程式。对于2D的情况，Fk是Kx和Ky二者的函数。
上文描述的本发明实施例与对接收器记录信号消除虚反射相关联。 然而，本发明也能用于震源一侧消除虚反射。在原理上，这能通过在震 源处直接测量流体压强来完成，但在实践中，可能利用互易原理以免得 需要这样做。
互易原理是波传播的基本原理，表述为：信号不受源与接收器位置 和特性互换的影响。例如，如果以震源在A点和接收器阵列在B点的测 量安排给出在接收器阵列处的某个信号，则使用在A点的单个接收器和 在B点的震源阵列将会导致同一个信号，如果该震源阵列对应于该接收 器阵列的话。(这里“对应”是指震源阵列中含有的震源个数与接收器阵 列含有的接收器数相同，而且在震源阵列中各震源安排的彼此相对位置 与接收器阵列中接收器的安排相同。
互易原理的结果是，能通过使用接收器阵列测量对震源阵列发出的 信号进行消除虚反射，不论这些震源同时起爆还是相继起爆。所记录的 信号类似于震源和接收器位置交换的互易安排中记录的信号。所以，上 文概括的方法也能应用于希望消除虚反射的震源阵列。测量由多个震源 组成的阵列在接收器阵列处产生的信号，并由这一测量的信号导出消除 虚反射滤波器。借助互易原理，这一滤波器能用于对震源阵列发出的信 号消除虚反射。
在本申请中描述的方法是快速的，使用本发明的方法处理完整波形 道是可能的。然而，如果希望的话，本发明可应用于一个波形道的选定 部分。
图7是能实现根据本发明的方法和装置20的示意方块图。
装置20包含一个可编程数据处理器21，带有程序存储器22，例如 只读存储器(ROM)形式的，其中存储的程序用于控制数据处理器21 以本发明的方法处理地震数据。该装置进一步包含非易失读/写存储器 23，用于存储例如在断电时必须保留的任何数据。由随机存取存储器 RAM24向数据处理器提供“工作”存储器或“暂存”存储器。提供了输入设 备25，例如用于接收用户命令和数据。提供了一个或多个输出设备26， 用于例如显示与处理过程的进展和处理结果有关的信息。输出设备可以 是例如打印机、可视显示单元或输出存储器。
可通过输入设备25或可选地通过机器可读数据存储器27提供地震 模型、地震测量结构安排的参数以及先有AVO信息。处理结果可经由输 出设备26输出或可以存储。
操作该系统和实现这里描述的方向所用的程序存储在程序存储器22 中，它可作为半导体存储器被嵌入，例如公知的ROM型存储器。然而， 该程序可以很好地存储在任何其他适当的存储介质中，如磁数据载体22a (如“软盘”)或CD-ROM。
附录
Robertsson和Kragh(上文)已表明质点速度的垂直分量vz可由压 强场的原始测量p估计出来。在典型的地震频带，所给出的一种可能近 似是：
$v_z(\omega ,x,y,z)=-\frac{i}{\mathrm{\rho \omega z}}[A_{00}^{\prime }+A_{10}^{\prime }{\left(\mathrm{kz}\right)}^2+A_{11}^{\prime }{z}^2(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})]p(\omega ,x,y,z),---\left(A1\right)$
标量
A′00＝1；A′10＝-1/2；A′11＝-1/2，                    (A2)
在式(A1)中，ω是(角)频率，K＝ω/α是波数，α是P波速度，ρ 是密度，而z＝z(x，y)是拖缆深度。根据Robertsson和Kragh，近似要 求拖缆被拖曳不深于约λ/3.5，这里λ是最小波长，还要求接收器空间采 样间隔约3m或更小。该方法基本上限定在低于第一虚凹槽的频率。例如， 对于深度6m的拖缆，对于高达约70Hz，式(A1)是准确的。
对于密集采样的压强测量，式(A1)的一个吸引人的特性是用于估 计垂直质点速度的滤波器是局部的和短的。利用三点中心差分算子近似 二阶水平空间导数，得到的空间滤波器有带五个系数的交叉结构。处理 单元拖缆压强数据，如由新的Q-海上拖缆在2.5D地球模型下提供的数 据，该滤波器只有三个系数。这样，为估计特定接收位置处的垂直质点 速度，只需要在这一接收器及两个相邻接收器处的压强场。假定(近水 平)拖缆在水面下的深度沿拖缆所有点处均为已知，则式(A1)给出在 强浪环境中估计垂直质点速度的方法。设P和Vz分别代表p和Vz的空 间付立叶变换。于是式(A1)中给出的估计能与原始压强测量结合以对 数据消除虚反射(例如，见L.Amundsen在Geophysics第56卷第 1027-1039页(1993)或UK专利申请9906456.0)。
$\tilde{P}(\omega ,k_x,k_y,z)=\frac{1}{2}[P(\omega ,k_x,k_y,z)-\frac{\mathrm{\rho \omega }}{k_z}{V}_z(\omega ,k_z,k_y,z)],---\left(A3\right)$
这里， 代表消除虚反射后的(上行)压强， $k_z=\sqrt{k^2-{\kappa }^2}$是 垂直波数，κz＝kx 2+ky 2，而kx和ky是水平波数。式(A3)对所有波数都 有效。另一种作法是，对式(A3)的近似如A.Osen等在“实现对海底地 震数据去多重(demultiple)和波场分裂的最佳空间滤波器”(提交给 Geophysics(2002))一文中建议的那些。
在本申请中，通过引入最优化滤波器，扩展了式(A1)中给出的垂 直质点速度场近似的准确度。L.Amundsen等在Geophysics第60卷第 212-222页(1995)中导出质点速度垂直分量和压强场之间的一个严格的 频率—波数域关系。对于单个拖缆数据，该关系是：
$V_z(\omega ,k_z,k_y,z)=-\frac{k_2}{\mathrm{\rho \omega }}\left(\frac{\exp (-i{k}_{z}z)+\exp \left(i{k}_{z}z\right)}{\exp (-i{k}_{z}z)-\exp \left({\mathrm{ik}}_{z}z\right)}\right)·P(\omega ,k_x,k_y,z)---\left(A4\right)$
当接收器电缆水平而且空气/水表面平坦而其压强消失时，式(A4) 对水柱中所有拖缆深度z都是有效的。实质上，该方程式表明，为从压 强场找出质点速度的垂直分量，首先该压强是要在接收器消除虚反射的， 然后要应用质点速度垂直分量的接收器虚反射算子。(考虑自由表面的 反射系数是-1，压强的虚反射算子是1-exp(2ikzz)，因为压强是上行和 下行波之和，而对垂直质点速度，虚反射算子是1+exp(2ikzz)，因为垂 直质点速度与上行和下行波之差成比例。)注意到当震源深度小于接收器 深度时，入射压强波场(直达压强波场及其虚反射)不含有接收器虚反 射，而只含有震源虚反射。由于式(A4)依赖于接收器消除虚反射，所 以它不能正确处置入射波场。在拖缆数据的接收器消除虚反射中这是不 重要的。然而，当震源深度大于接收器深度时，入射压强波场含有接收 器虚反射作为压强场的被反射部分。在这种情况中，式(A4)对全压强 波场有效。
式(A4)可被重写为
$V_z(\omega ,k_x,k_y,z)=-\frac{i}{\mathrm{\rho \omega z}}\underset{m=0}{\overset{\infty }{\Sigma }}{(-1)}^m{F}_m(k,z){\kappa }^{2m}·P(\omega ,k_x,k_y,z),---\left(A5\right)$
其中独立于水平波数(κ)的函数
$F_m=F_m(k,z)=k^{-2m}\underset{n=m}{\overset{\infty }{\Sigma }}{A}_{\mathrm{nm}}{\left(\mathrm{kz}\right)}^{2n}---\left(A6\right)$
这里
$A_{\mathrm{nm}}=\frac{{(-1)}^n{2}^{2n}{B}_{2n}}{\left(2n\right)}\left(\begin{array}{c}n\\ m\end{array}\right)---\left(A7\right)$
是标量，Bn是第n个伯努里数，n！是n的阶乘。所有函数Fm能以紧微 形式给出。例如，
$F_0=\mathrm{kz}\cot \left(\mathrm{kz}\right)=-\mathrm{ikz}\left(\frac{\exp (-\mathrm{ikz})+\exp \left(\mathrm{ikz}\right)}{\exp (-\mathrm{ikz})-\exp \left(\mathrm{ikz}\right)}\right)---\left(A8\right)$
$F_1=\frac{z}{2k}\frac{\partial }{\partial \left(\mathrm{kz}\right)}{F}_0,---\left(A9\right)$
$F_2=\frac{z}{{4k}^3}(\frac{\partial }{\partial \left(\mathrm{kz}\right)}\mathrm{kz}-2)\frac{\partial }{\partial \left(\mathrm{kz}\right)}{F}_0·---\left(A10\right)$
函数F0保证垂直传播的波(kx＝ky＝0)被正确地消除虚反射。
对水平波数取式(A5)的逆付立叶变换，给出 $v_z(\omega ,x,y,z)=-\frac{i}{\mathrm{\rho \omega z}}\underset{m=0}{\overset{\infty }{\Sigma }}{F}_m(k,z)D_2^m(x,y)*p(\omega ,x,y,z),---\left(A11\right)$
这里D2是κ2的逆付立叶变换，*是2D(二维)空间褶积算子。对于有限 带宽，
D2(x，y)＝D2(x)+D2(y)，                               (A12)
$D_2\left(x\right)=\frac{{\partial }^2}{{\partial x}^2};D_2\left(y\right)=\frac{{\partial }^2}{{\partial y}^2},---\left(13\right)$ 是二阶空间微分算子。然而，微分涉及对高波数的大放大。因此，在实 际场合，优选使用该微分的有限带宽版本。设d2代表D2的有限带宽版本。 两个可能的有限带宽微分算子是：
$d_2\left(x\right)=\frac{1}{\mathrm{\Delta x}}[\delta (x+\mathrm{\Delta x})-2\delta \left(x\right)+\delta (x-\mathrm{\Delta x})],---\left(14\right)$
这里ΔX是空间采样间隔，或
$d_2(x,\omega )=\frac{1}{\pi }\underset{0}{\overset{K}{\int }}{\mathrm{dk}}_x\cos \left(k_{x}x\right)W(\omega ,k_x)(-{k_x}^2),---\left(A15\right)$
这里W是适当限定微分算子带宽的加权函数。注意到能对每个频率分量 改变W。在式(A15)中，K是最大波数。两个可能的选择是K＝ω/C， 这里C是水的速度，或K＝π/ΔX，为Nyquist(奈奎斯特)波数。
注意到式(A11)对水柱中所有接收器深度Z有效。再有，注意到该 模式是迭代的：
$D_2^m*p=D_2^m*(D_2^{m-1}*p);D_2^0=1.---\left(A16\right)$
这一递归性质允许式(A11)以有效的数值计算实现。对于密集空间采样， 微分D2(x)和D2(y)可由长度三滤波器近似[见式(A14)]。一个适 当设计的、有限带宽的二阶微分器d2(x)对于得到稳定的迭代模式必不 可少的。
特例I：有限m
当m为有限时，例如M，则产生式(A11)的一个特例。于是，对 于有限带宽，
$v_z(\omega ,x,y,z)\approx -\frac{i}{\mathrm{\rho \omega z}}\underset{m=0}{\overset{M}{\Sigma }}{F}_m(k,z){(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})}^{m}p(\omega ,x,y,z),---\left(A17\right)$
这里Fm是在式(A6)中作为无限级数给出。然而，函数Fm也能利用 三角函数给出分析表达式。实例是式(A8)、(A9)和(A10)中给出的 F0、F1和F2。注意到用于导出式(4)的平静海平面假定现在可以被放松 了。在式(A17)中，x和y变量与z解耦，意味着它对强浪表面也是有 效的。
式(A17)能以不同的方式以数值计算来实现。优选的实现方式取决 于压强场的空间采样间隔。对于密集采样的数据，通常将利用经点的中 心有限差分近似来代表空间导数。如果M小，则滤波器是局部的和紧致 的。作为一例，选则M＝2，对无限带宽给出：
$v_z(\omega ,x,y,z)\approx -\frac{i}{\mathrm{\rho \omega z}}[F_0+F_1(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})+F_2{(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})}^2]p(w,x,y,z).---\left(A18\right)$
利用三点有限差分算子近似二阶空间导数，则迭代滤波器的长度为 三。由于函数F0、F1和F2对所有接收器深度都有效，故式(A18)比式 (A1)更具一般性，这是对强浪花消除虚反射的优选作法。
对于较稀疏的采样间隔，水平导数最方便地实现为数值最优化差分 算子，或利用付立叶技术。
特例II：有限m，n
对式(11)中的对m求和和式(A6)中的对n求和从m和n的无 限值分别限定为M和N(m)，对无限带宽给出近似：
$v_z(\omega ,x,y,z)\approx -\frac{i}{\mathrm{\rho \omega z}}\underset{m=0}{\overset{M}{\Sigma }}{F}_m(k,z){(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})}^{m}p(\omega ,x,y,z),---\left(A19\right)$
其中
$F_m\approx k^{-2m}\underset{n=m}{\overset{N\left(m\right)}{\Sigma }}{A}_{\mathrm{nm}}{\left(\mathrm{kz}\right)}^{2N}---\left(A20\right)$
选择M＝1，N(0)＝2和N(1)＝1，给出
$F_0\approx 1-\frac{{\left(\mathrm{kz}\right)}^2}{3},---\left(A21\right)$
$F_1\approx -\frac{z^2}{3}---\left(A22\right)$
以及
$v_z(\omega ,x,y,z)\approx -\frac{i}{\mathrm{\rho \omega z}}[{\Lambda }_{00}+A_{10}{\left(\mathrm{kz}\right)}^2+A_{11}{z}^2(\frac{{\partial }^2}{{\partial x}^2}+\frac{{\partial }^2}{{\partial y}^2})]p(\omega ,x,y,z),---\left(A23\right)$
其中标量
A00＝1；A10＝-1/3；A11＝-1/3                            (A24)
式(A23)有与Robertsson和Kragh(上文)提议的式(A1)相同的结 构。然而，标量A10和A11不同于式(A1)中的A’10和A’11。在许多方面， 式(A23)的限制类似于Robertsson和Kragh(上文)的公式的那些限 制。然而，式(A11)、(A17)和(A19)是一般性的，而且一个显著差 别是它们对所有拖缆深度都有效，而由Robertsson和Kragh(上文)导 出的理论要求拖缆被拖曳不深于约λ/3.5，这里λ是最小波长。
如果选择M＝1，N(0)＝∞和N(1)＝2，这造成F0＝kzcot(kz)和
$F_1=-\frac{z^2}{3}(1+\frac{2}{15}{k}^2{z}^2)·$
单拖缆数据的数据处理
上文的式(A11)、(A17)和(A19)包括沿y方向和沿交叉线方向 的导数(假定该拖缆沿X方向延伸)。在一个只包括单个拖缆的地震测量 安排中，不清楚如何可以由拖缆上的接收器采集的地震数据来确定上述 方程式中出现的沿y方向的交叉线导数。然而，将会注意到，只发生偶 数阶(二阶、四阶、六阶等)空间导数。所以，找出二阶空间导数的表 达式便足够了。
考虑图4中显示的单拖缆几何结构，它是图1的地震测量安排的平 面示意图。震源1位于x-y平面中的坐标(o，yo)处。拖缆沿X轴延伸， 并有y坐标y0。接收器之一Si位于坐标(xo，yo)，并假定震源1和接收器 Si位于一个声层中的同一Z平面。将假定模型只随深度变化，不随x和 y变化。通过对二阶导数使用有限差分近似，可能对压强P写出下式：
$\frac{{\partial }^2}{{\partial y}^2}p\left(x_{0,}{y}_0\right)=\frac{1}{{\mathrm{\Delta y}}^2}(p(x_0,y_0+\mathrm{\Delta y})-2p(x_0,y_0)+p(x_0,y_0-\mathrm{\Delta y}))+O\left({\mathrm{\Delta y}}^2\right),---\left(A25\right)$
这里O(Δy2)代表有限差分近似中的主导误差项。然而，因为径向对 称，在接收器两侧记录的压强相同，于是式(A25)变为
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{2}{{\mathrm{\Delta y}}^2}(p(x_0,y_0+\mathrm{\Delta y})-p(x_0,y_0))+O\left({\mathrm{\Delta y}}^2\right).---\left(A26\right)$
径向对称还造成
$p(x_0,y_0+\mathrm{\Delta y})=p(x_0\sqrt{1+{(\mathrm{\Delta y}/x_0)}^2},y_0),---\left(A27\right)$
于是式(A26)变为
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{2}{{\mathrm{\Delta y}}^2}(p(x_0\sqrt{1+{(\mathrm{\Delta y}/x_0)}^2},y_0)-p(x_0,y_0))+O\left({\mathrm{\Delta y}}^2\right)---\left(A28\right)$
通过沿X方向使用一阶单侧有限差分近似，有可能把式(A28)写成：
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{{2x}_0(\sqrt{1+{(\mathrm{\Delta y}/x_0)}^2}-1)}{{\Delta y}^2}(\frac{\partial }{\partial x}p(x_0,y_0)+O(\sqrt{1+{(\mathrm{\Delta y}/x_0)}^2}-1))+O\left({\mathrm{\Delta y}}^2\right)---\left(A29\right)$
取平方根的泰勒级数展开：
$\sqrt{1+{(\mathrm{\Delta y}/x_0)}^2}=1+\frac{1}{2}{(\mathrm{\Delta y}/x_0)}^2+O\left({\mathrm{\Delta y}}^4\right),---\left(A30\right)$
式(A29)可重写成：
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{{2x}_0(\frac{1}{2}{(\mathrm{\Delta y}/x_0)}^2+O\left({\mathrm{\Delta y}}^4\right))}{{\mathrm{\Delta y}}^2}(\frac{\partial }{\partial x}p(x_0,y_0)+O\left({\mathrm{\Delta y}}^2\right))+O\left({\mathrm{\Delta y}}^2\right)---\left(A31\right)$
对式(A31)进一步简化，产生
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{1}{x_0}\frac{\partial }{\partial x}p(x_0,y_0)+O\left({\mathrm{\Delta y}}^2\right)·---\left(A32\right)$
最后，令Δy→0，于是：
$\frac{{\partial }^2}{{\partial y}^2}p(x_0,y_0)=\frac{1}{x_0}\frac{\partial }{\partial x}p(x_0,y_0)·---\left(A33\right)$
对于水平分层介质上的压强数据，方程式(A33)中不存在奇异性问题， 因为 趋于零快于奇异性1/xo趋于无穷。
式(A33)与对圆柱对称有效的数据(在水平分层介质上的数据)是 一致的。所以，另一个途径是以 乘数据(这里t是时间)，并使用本申 请以及UK专利申请0015810.5号中导出的理论的2D(二维)版本。然 而，式(A33)有潜力给出更准确的结果。
应该指出，对于不是简单平面分层的介质，上述途径可能不可靠， 因为可能由于p相对于x不对称而发生问题。
复频率
在实现滤波器时，我们使用复频率以避免函数Fm中的奇异性，如 式(8)、(9)和(10)中分别对M＝0，1，2给出的那样。使用复频率由例 如S.Mallick和L.N.Frazer在Geophysics第52卷第1355-1364(1987) 中描述和由L.Amundsen和B.Ursin在Geophysics第56卷第1027-1039 (1991)中描述。
关于使用数值最优化空间滤波器Vz估计
1.引言
在这一部分，我们引入数值最优化空间滤波器，用于估计垂直质点 速度场。Amundsen等(1995，上文)导出质点速度垂直分量和压强场之 间的严格的频率—波数域关系。对于单拖缆数据，在2D(二维)情况中 该关系是：
$v_z(\omega ,k_x,z)=-\frac{k_z}{\mathrm{\rho \omega }}\left[\frac{\exp (-i{k}_{z}z)+\exp \left(i{k}_{z}z\right)}{\exp (-i{k}_{z}z)-\exp \left(i{k}_{z}z\right)}\right]·p(\omega ,k_x,z)---\left(A34\right)$
式(A34)可重写成：
vz(ω，kx，z)＝F(ω，kx，z)·p(ω，kx，z)                (A35)
这里Vz估计滤波器由下式给出：
$F(\omega ,k_x,z)=-\frac{k_z}{\mathrm{\rho \omega }}\left[\frac{\exp (-{\mathrm{ik}}_{z}z)+\exp \left({\mathrm{ik}}_{z}z\right)}{\exp (-{\mathrm{ik}}_{z}z)-\exp \left(i{k}_{z}z\right)}\right]---\left(A36\right)$
在频率—空间域中，式(A34)以符号写成：
$V_z=\tilde{F}*p---\left(A37\right)$
这里*代表空间褶积。
2.最优化问题的公式化
一个空间数字滤波器的离散水平波数响应能写成：
${\tilde{F}}_k=\tilde{F}\left(\mathrm{k\Delta }{k}_x\right)=\frac{\underset{m=-M/2}{\overset{M/2}{\Sigma }}{b}_m\exp (-\mathrm{ik\Delta }{k}_x\mathrm{m\Delta x})}{\underset{n=0}{\overset{N}{\Sigma }}{a}_n\exp (-\mathrm{ik\Delta }{k}_x\mathrm{n\Delta x})}=\frac{{\Theta }_{\mathrm{km}}{b}_m}{{\Phi }_{\mathrm{kn}}{a}_n}---\left(A38\right)$
这里an，bm是以N、M为各自滤波器阶数的后向和前向滤波器系数。 不失掉一般性，M和N假定为偶数且a0＝1。再有，Δx是空间采样间隔， Δkx＝π/[Δx(K-1)]是离散水平波数之间的距离。如果Fk＝F(kΔkx) 代表所希望的离散水平波数响应，那么空间滤波器的设计能被表述为一 个一般的加权非线性最小二乘最优化问题：
${\tilde{F}}_k=\underset{b_m,a_n}{\min }\underset{k=-(K-1)}{\overset{K-1}{\Sigma }}{W}_k{|\frac{{\Theta }_{\mathrm{km}}{b}_m}{{\Phi }_{\mathrm{kn}}{a}_n}-F_k|}^2---\left(A39\right)$
这里Wk是一组正加权。接下来，最优化问题只对传播波公式化。对于 传播波，理想的分解滤波器是纯实数的(即零相位)。令k表示与接近临 界波数的一个离散水平波数相对应的标号。
3.零相位FIR滤波器
对于零相位非递归或有限冲激响应(FIR)滤波器，前向滤波器系数 是对称的，而反向滤波器阶为零。这使得式(A39)中给出的非线性最小 二乘最优化问题简化成如下线性最小二乘最优化问题：
$\underset{g_m}{\min }\underset{k=0}{\overset{\kappa -1}{\Sigma }}{W}_k{|{\Gamma }_{\mathrm{km}}{g}_m-F_k|}^2---\left(A40\right)$
这里Γkm＝cos(ΔKxmΔx)。gm与bm关联，对于gm＝1，2，…，M/2，bm＝0.5gm， 而b0＝go。
4.二次规划
式(A40)中显示的线性最小二乘最优化问题是无约束的，唯一的和 最佳的解由伪逆变换给出。然而，可能重要的是最佳FIR滤波器的振幅 响应在一些特定水平滤数处有零误差，或者说它们的振幅响应在某一水 平波数范围小于一个给定函数。线性最小二乘最优化问题与线性等式和/ 或不等式约束相结合，称作“二次规划(QP)”，能以许多不同的算法有效 地解决。
例如，在零水平波数能包括一个单个等式约束，以保证由最佳FIR 滤波器正确地估计传播波。此外，在耗损水平波数范围中包括若干线性 不等式约束，以形成该空间滤波器的理想分解滤波器振幅响应。这防止 由数值上最优化的空间滤波器导出耗损波。