Kirchhoff型高保真反偏移理论

孙建国

(1.吉林大学地球探测科学与技术学院,吉林长春130026;2.吉林大学国土资源部应用地球物理综合解释理论开放实验室-波动理论与成像技术实验室,吉林长春130026)

0 引言

反射地震资料反偏移成像的基本目的是将深度偏移像场(深度偏移的输出结果)再反映射到时间域中去[1-2]。因此,在实际功效上,反偏移成像相当于是一种以深度偏移像场作为输入数据的数值模拟。然而,利用数值模拟只能得到与真实物理过程相对应的数据体(例如共炮点或共场源数据体),而利用反偏移成像既可以构造出与真实物理过程相对应的数据体,也可以构造出不对应任何真实物理过程的数据体(例如共炮检距或其他类型的数据体)。再者,利用数值模拟可以计算出波动过程中的各种现象(例如传播、透射、反射以及散射等)及其随时间的演化历史,而利用反偏移成像只能构造出反射地震观测数据而不能再现形成这些数据的物理过程。

如果在反偏移成像计算中所使用的速度模型和观测系统与在偏移成像中所使用的速度模型和观测系统完全相同,则对于相同的输出波场类型,反偏移成像将转化为偏移成像的逆运算[1-3]。反之,如果在反偏移成像过程中所使用的速度模型、观测装置以及输出波场类型与在偏移过程中所使用的有所不同,则反偏移成像的输出将形成一个与偏移输入完全不同的数据体。这一特点奠定了利用反偏移成像解决若干反射地震资料处理问题的基础。事实上,在均匀介质中,基于微分方程的DMO理论就是利用偏移与反偏移串联才得以建立的[4],而在一般非均匀介质条件下向零炮检距偏移的 Kirchhoff型理论更是利用了反偏移成像的直接结果[5]。

对于实际应用来讲,反偏移成像研究的重要意义主要体现在两个方面。首先,反偏移研究可以为速度分析、波型转换、观测系统变换、观测数据规格(则)化以及数值模拟提供非常规的理论基础和计算技术[6-8]。其次,通过与深度偏移的反复串联应用,反偏移理论可以为提出和发展新的反射地震成像方法提供出发点[1-2]。

反偏移的基本思想与偏移类似,可以通过不同的方式实现。事实上,反偏移既可以表示成为沿等时面的加权叠加积分[1-2,8],又可以表示成为频率-波数域中的积分[9]。如果反偏移的基本思想是利用沿等时面的加权叠加实现的,则称其为等时面加权叠加反偏移或 Kirchhoff型反偏移。如果在等时面加权叠加中所利用的权函数是真振幅权函数,则称其为 Kirchhoff型真振幅反偏移。在文献中,关于Kirchhoff型真振幅反偏移的基本理论发表于20世纪90年代。与 Kirchhoff型偏移类似,在建立Kirchhoff型反偏移的过程中利用了一种与反射面无关的策略,即假设地下的每一个点都是反射点。另外,为了分析反偏移的输出结果假设偏移像场的子波具有很短的延续。如果偏移像场的子波是δ脉冲,则真振幅反偏移结果与数值模拟结果完全一致。然而,如果偏移像场的子波延续是有限的,则反偏移的输出结果将与数值模拟有所不同。事实上,根据对 Kirchhoff型真振幅反偏移的输出场所做的稳相分析[10],反偏移场由两个因子的乘积组成,其中一个是反偏移信号(震源子波与反射系数和几何扩散因子的乘积),另一个是振幅畸变因子。在数学上,真振幅反偏移信号是利用了真振幅加权函数的直接结果,而振幅畸变因子的出现完全是由与反射面无关的实现过程所引起的。这意味着即使在不使用任何加权函数的条件下,振幅畸变现象也是存在的。因此,靠反偏移本身无法消除振幅畸变效应。如果一定要消除振幅畸变效应所带来的振幅失真现象,必须在反偏移之后再对反偏移振幅进行校正。只有在经过了这种校正之后,反偏移的输出结果才与数值模拟的输出结果完全一致。

笔者的目的就是要提出一种实现上述反偏移振幅校正的方法,以使得反偏移成像的结果与数值模拟结果完全一致。考虑到在经过了振幅校正之后的反偏移场应该只由反偏移信号组成这一事实,将所提出的带有振幅畸变校正的 Kirchhoff型反偏移称为高保真反偏移。与常规的真振幅反偏移相比,在高保真反偏移中考虑了反偏移实现策略对反偏移结果的影响。在理论上,高保真反偏移应包括3部分:①常规真振幅反偏移;②振幅畸变因子提取;③振幅畸变因子校正。由于振幅畸变因子与稳相点和反射点处的真振幅加权函数和深度差函数的二阶导数矩阵有关,所以为了提取振幅畸变因子必须首先确定稳相点和反射点的位置。鉴于 Kirchhoff型反偏移算子与 Kirchhoff型偏移算子在数学形式上的类似性,将采用在 Kirchhoff型偏移中提取稳相点的方法,即多重加权叠加法[11-12]。一旦确定了稳相点的位置,就可以根据反射点是一个特殊的稳相点这一事实确定出反射点的位置,进而计算出振幅畸变因子的具体数值。另外,还可以根据在反射点上振幅畸变因子恒等于1这一事实,利用不同时间点上的反偏移场与反射时间点上的反偏移场之间比值的方法提取振幅畸变因子。一旦得到了振幅畸变因子,就可以通过简单的除法实现振幅畸变校正。

根据上述讨论,在具体实施算法中高保真反偏移应包括3个步骤:①多重加权叠加反偏移;②稳相点和振幅畸变因子求取;③振幅畸变校正。为此将首先对等时面加权叠加及其输出场进行必要回顾,然后再具体讨论如何利用多重加权叠加反偏移计算振幅畸变校正因子以及如何进行振幅畸变校正。

1 等时面加权叠加及其输出场

与绕射加权叠加偏移类似,在等时面加权叠加反偏移中假设速度模型的不均匀程度符合几何射线理论的应用条件,即在一个波长的范围之内速度呈平缓变化。此外,假设地下介质为由一系列光滑曲面分开的不均匀层状介质组成;在每层内,速度为坐标的连续函数,同时具有横向和垂向不均匀性[13]。

为了能够建立等时面加权叠加的数学公式,引入2个坐标系。

(1)用于确定地下空间中任意一点位置的全局坐标系

(2)用于确定反偏移输出数据体中任意一点位置的全局坐标系

在坐标系(r,z)中,r=(r1,r2)=(x,y)代表水平坐标矢量,z代表垂直深度,其正方向指向地下(图1)。为了方便起见,将由坐标系(r,z)所表示的地下空间称为模型空间(深度域)。类似地将由坐标系(ξ,t)所描述的空间称为数据空间(时间-空间域)。其中,参数坐标矢量ξ=(ξ1,ξ2)用于确定任意一个地震道的位置,而t代表时间,向上为正(图2)。根据定义,参数矢量ξ与在模型空间中所给定的炮点坐标rs和接收点坐标 rg有一定的函数关系[12,14]。另外,时间变量t所给出的是双程传播时间,即

式中:ts(r,z)为由炮点 rs到点(r,z)的传播时间;tg(r,z)为由接收点 rg到点(r,z)的传播时间。对于一对给定的炮点和接收点,在模型空间中的每一个点上都有一个双程传播时间。由具有相同双程传播时间的地下点所组成的曲面称为等时面。在Kirchhoff型反偏移理论中,假设等时面是光滑曲面,其各阶导数为空间坐标(r,z)的连续函数[13]。一般来讲,等时面可以用隐式方程来描述

图1 模型空间Fig.1 Model Space

图2 数据空间Fig.2 Data Space

式中:下标I为等时面;zI为等时面上任意一点的深度。考虑到 rs和 rg均为ξ的函数,即

可以将上述等时面的隐式方程写为

从这个方程中解出垂直坐标zI,可以得到等时面的显式方程

在建立 Kirchhoff型反偏移理论的过程中并不需要知道 zI=zI(r;ξ,t)的具体数学形式,而在Kirchhoff型反偏移的具体实现过程中,等时面zI= zI(r;ξ,t)是利用射线追踪计算出来的。

类似地将目标反射面表示为

式中:下标R代表反射。为了方便起见,将反射面zR上的任意一点称为 PR。根据反射定律,不同的点 PR将对应于不同的炮点和接收点对(rs,rg)= (rs(ξ),rg(ξ))。因此,PR与参数矢量ξ有一一对应关系。这意味着,可以把反射面 zR(rR)写成ξ的函数,即

现在引入等时面加权叠加的数学表示。令 Q和P分别代表数据空间和模型空间中的任意一点,其坐标分别为(ξ1Q,ξ2Q,tQ)(图 2)和(r1P,r2P,zP) (图1)。同时,令 ∑IQ代表由 Q点的参数坐标(ξ1Q,ξ2Q)和双程走时tQ定义的等时面

根据这个定义,在任意等时面上的点为

而在等时面ΣIQ上的点为 PIQ=(r,zIQ)。

式中:dr=dr1dr2;Ξ为反偏移孔径;Wd(r,z)为真振幅加权函数;O(r,z)为深度偏移像场。为了便于分析研究,假设O(r,z)是由 Kirchhoff型真振幅偏移产生,即[15]

式中:CR为反射点 PR处的平面波反射系数;s(z)为映射到了深度域中的震源子波;Λm为子波拉伸因子。

式(1)说明,Kirchhoff型反偏移在数学上所代表的是一个多点对一点的变换,即通过加权叠加的方式将位于等时面上所有点处的偏移像场O(r,zI)变换成数据空间中一个点(Q点)上的反偏移像场。为了求出反偏移场,现将其扩展成为z的函数[10],并使其在 z=0时与相等,即。可以证明,满足这个条件的可以表示成为下列积分(文献[10]中公式(3))

为了方便研究起见,定义 l(r,tQ)=zIQ-zR为深度差函数[10]。

式中:kz为垂直波数;rst为稳相点的水平坐标; L(rst,tQ)为深度差函数 l(r,tQ)在稳相点 rst处的二阶导数矩阵;Wd(rst,zIQ)为在点(rst,zIQ(rst;ξ;tQ))处的真振幅加权函数值,因子 B(rst,kz)和χL分别定义为(文献[10]中式(25)、(26))

式中:函数sgn[L(rst,tQ)]定义为矩阵L(rst,tQ)的正负符号差,即矩阵的正特征值个数减去负特征值的个数。

式中:L(rR,tQR)和 Wd(rR,zIQR)分别为深度差函数和真振幅加权函数在反射点(rR,zR)上的二阶导数矩阵值和真振幅加权函数值。

式(6)是在渐近分析意义下 Kirchhoff型反偏移的输出场。与数值模拟的结果相比,反偏移的输出结果多了一个因子 Ed。因此,高保真反偏移的目的就是要校正消除因子 Ed对真振幅反偏移结果的影响。换句话说,高保真反偏移的目的是为了得到反偏移信号。根据式(7),振幅畸变因子 Ed是2个因子的乘积,一个是在反射点(rR,zIQR)处的因子|det(L(rR,tQR))|1/2/|Wd(rR,zIQR)|,另一个是在稳相点(rst,zIQ)处的因子|det(L(rst,tQ))|1/2/ |Wd(rst,zIQ)|。由于反射点本身就是一个稳相点,所以为了实现高保真反偏移必须首先要确定稳相点的空间位置。一旦求出稳相点的位置,就可以根据在反偏移实现过程中所计算出的等时面和真振幅加权函数计算出振幅畸变因子 Ed。

2 稳相点和反射点坐标

为了求出稳相点坐标,采用 Bleistein[16]或Tygel等[11]提出的多重加权思想。在历史上,多重加权叠加法是Bleistein为了反演反射角等几何或物理量所提出的一种在偏移的同时进行反演的方法,其基本思想是把所要反演的几何或物理量作为一个绕射叠加偏移(Kirchhoff型偏移)中的加权函数,并将如此加权之后的偏移结果与不加权时的偏移结果进行对比。随后,Tygel等将Bleistein的基本思想用于求解反射点的水平坐标。然而,Tygel等的多重加权叠加偏移是在时间域中进行的,其有关计算公式对于物理可实现信号不再成立[12]。2003年,笔者对利用绕射加权叠加偏移求取反射点的水平坐标问题进行了进一步讨论,指出最好在频率域中进行反射点计算[12]。

正如在引言中所提到的,Kirchhoff型反偏移式(1)与 Kirchhoff型偏移公式具有类似的数学形式,因此可以利用多重加权叠加的基本思想求取反偏移中的稳相点。参照 Tygel等[11]关于多重加权绕射叠加的工作,令

如果用这个加权函数矢量中的各分量代替式(1)中的真振幅加权函数,则在经过等时面叠加后得到3个类似于式(3)所给出的反偏移场,即

式中

如果在给定的坐标系中等时面函数的水平坐标恒大于零,还可以利用反偏移场的模值计算稳相点的水平坐标。

式(13)、(14)所给出的是稳相点的水平坐标,而相应的深度坐标对应于深度差函数 l(r,tQ)的极值点深度。由于在数据空间中的一个时间子波是一族等时面共同作用的结果[9],所以可以在目标反射面附近求出一组稳相点,而其中具有最小深度的稳相点即为反射点。

3 振幅畸变因子及其校正

一旦得到稳相点和反射点的坐标,就可以利用式(7)计算振幅畸变因子。然而,式(7)不仅含有在反偏移过程中需要计算的真振幅加权函数,也含有反偏移过程中不需要计算的距离差函数的二阶导数矩阵。因此,直接利用式(7)需要额外计算深度差函数及其二阶导数矩阵。如果在计算 rst的同时记录下其所对应的双程走时,则可以避免这些额外的计算。事实上,对于一个给定的目标反射界面,反射走时tR等于等时面族中具有最小双程走时的等时面值,而与最小双程走时相对应的稳相点水平坐标即为反射点的水平坐标。因此,可以把稳相点的水平坐标和振幅畸变因子表示为反射走时的函数。

对于一个给定的偏移信号层,在经过上述计算后可以得到一个以双程传播时间t为自变量的稳相点序列,即 rst(tR),rst(tR+Δt),rst(tR+2Δt),…,rst(tR+Ts)。其中,Ts为反偏移子波s～(t)的延续时间(子波长度),Δt为2个相邻等时面之间的时间差。显然,rst(tR)对应着偏移信号层的顶界面,而 rst(tR+Ts)对应着偏移信号层的底界面。因此,与rst(tR)相对应的振幅畸变因子为1[10]。根据这一事实,对于一个给定反偏移输出道上位于QR点附近的点列Qn,n= 1,2,3,…有

式中:虽然在时间域内有 E(tQ=tR+nΔt)=(Q)/(QR),但是,对于物理可实现的信号来讲(QR)≡ 0,所以利用反偏移场提取振幅畸变因子的计算只能在r-kz域中进行。

在利用上述两种方法中的任意一个求出E(tR+ nΔt)后,可以直接在时间域内对反偏移信号的振幅进行校正,即

4 高保真反偏移的实现步骤

根据上述讨论,高保真反偏移的具体实施流程如下。

(1)对偏移场做关于深度坐标的Fourier变换。

(2)计算源自于炮点和接收点的射线族并根据射线走时计算等时面zIQ。

(3)对变换后的数据利用相移算子exp(ikzzIQ)进行相移(滤波)处理。

(4)对滤波后的反偏移场利用权函数矢量(r1, r2,1,Wd)T进行矢量加权叠加。

(5)根据式(13)、(14)求出与各等时面(线)相对应的稳相点。

(6)根据式(7)或式(17)求出的振幅畸变因子;如果利用式(7),需要先计算出深度差函数的二阶导数矩阵。

(7)根据式(18)对反偏移场进行校正。

5 结语

作为一种现代反射地震成像方法,反偏移主要被用于反射地震资料的解释阶段。通过其与偏移的串联应用,可以利用反偏移解决再偏移(remigration)、观测系统变换以及速度建模等一系列问题。如果偏移像场是通过人工构建而不是通过偏移计算得到的,还可以利用反偏移实现数值模拟。在理论上,反偏移应该与数值模拟具有同等的效果。具体地说,在不考虑透射损失的条件下,在几何射线理论的框架中,Kirchhoff型反偏移成像的结果应该是目标反射面的反射系数与震源子波的乘积再除以相应的几何扩散因子。然而,由于 Kirchhoff型反偏移所采取的与反射面无关的策略,其输出结果与数值模拟结果相差一个振幅畸变因子。为了使 Kirchhoff型反偏移成像输出结果与数值模拟的输出结果相一致,笔者提出了一种提取和消除振幅畸变因子的方法,称为高保真反偏移。与常规的真振幅反偏移相比,高保真反偏移是一种与反射面有关的反射地震成像方法,在一次反偏移计算中利用了由4个加权函数所组成的加权函数矢量而并非一个单一的加权函数。此外,高保真反偏移由3个处理过程组成,即稳相点求取、振幅畸变因子求取和振幅畸变因子消除。根据这3个处理过程,可以把高保真反偏移算法分成3个组成部分。在第一部分中,通过采用加权函数矢量(r1,r2,1,Wd)T的方式得到4个反偏移场。其中,与前2个加权函数相对应的是含有稳相点空间坐标信息的反偏移场,与第3个加权函数(单位加权函数)所对应的是不加权时的反偏移场,而与最后一个加权函数所对应的是常规的真振幅反偏移场。在第二部分中,利用这4个反偏移场相互之间的比值求取稳像点和反射点的水平坐标以及相应的振幅畸变因子。最后,在算法的第三部分中,利用反射点上的反偏移场和非反射点上的反偏移场进行振幅畸变校正(振幅畸变因子消除)。

虽然笔者给出了高保真反偏移的基本理论和算法步骤,但是并没有讨论具体的实现技术和在数值实现过程中可能出现的问题。例如,在具体的数值实现过程中,式(13)～(16)都会出现分母为零的情况。因此必须采用一定的正则化处理以保证数值计算的稳定性,至于采用什么样的正则化方法是一个需要进一步研究的问题。另外,笔者提出了两种提取振幅畸变因子的途径,一种是直接利用式(7),另外一种是利用式(17)。如果利用式(7),需要确定稳相点的深度坐标并计算距离差函数的二阶导数矩阵。如果利用式(17),则既不需要确定稳相点的深度坐标,也不需要计算深度差函数的二阶导数矩阵,但是要处理有可能出现的分母为零问题。虽然在理论上看第二种途径要优于第一种,但是在实践中是否如此还需要进行深入的讨论。

[1] Hubral P,Schleicher J,Tygel M.A Unified Approach to 3-D Seismic Reflection Imaging;Part I,Basic Concepts[J].Geophysics,1996,61(3):742-758.

[2] Tygel M,Schleicher J,Hubral P.A Unified Approach to 3-D Seismic Reflection Imaging;Part II,Theory[J].Geophysics, 1996,61(3):759-775.

[3] Santos L T,Schleicher J,Tygel M,et al.Seismic Modeling by Demigration[J].Geophysics,2000,65(4):1281-1289.

[4] Deregowski S M,Rocca F.Geometrical Optics and Wave Theory ofConstantOffsetSections in Layered Media[J]. Geophysical Prospecting,1981,29(3):374-406.

[5] Tygel M,Schleicher J,Hubral P,et al.2.5-D True-amplitude Kirchhoff Migration to Zero Offset in Laterally Inhomogeneous Media[J].Geophysics,1998,63(2):557-573.

[6] Whitcombe D N.Fast Model Building Using Demigration and Single-step Ray Migration[J].Geophysics,1994,59(3):439-449.

[7] Ferber R G.Migration to Multiple Offset and Velocity Analysis[J].Geophysical Prospecting,1994,42(2):99-112.

[8] 孙建国.论三维等时线叠加反偏移中的有关问题[J].吉林大学学报:地球科学版,2002,32(3):273-278.

[9] 孙建国.均匀介质中的F-K反偏移:基本概念、基本公式及其在非均匀介质中的应用[J].吉林大学学报:地球科学版, 2008,38(1):135-143.

[10] Sun J G.The Stationary Phase Analysis of the Kirchhoff-type Demigrated Field[J].Applied Geophysics,2010,7(1):18-30.

[11] Tygel M,Schleicher J,Hubral P,et al.Multiple Weights in Diffraction Stack Migration[J].Geophysics,1993,58(12): 1820-1830.

[12] Sun J G.True-amplitude Weight Functions in 3D Limited-aperture Migration Revisited[J].Geophysics,2004,69(4):1025-1036.

[13] 孙建国.Kirchhoff型真振幅偏移与反偏移[J].勘探地球物理进展,2002,25(6):1-5.

[14] Schleicher J,Tygel M,Hubral P.3-D True-amplitude Finiteoffset Migration[J].Geophysics,1993,58(8):1112-1126.

[15] Sun J G.Limited-aperture Migration[J].Geophysics,2000,65 (2):584-595.

[16] Bleistein N.On the Imaging of Reflectors in the Earth[J]. Geophysics,1987,52(7):931-942.