柯西波场及在角度域成像条件中的应用

王雄文， 王华忠， 钱建良

1 同济大学海洋与地球科学学院， 上海　200092 2 Department of Mathematics, Michigan State University, East Lansing, MI, 48824， U.S.A.



柯西波场及在角度域成像条件中的应用

王雄文1， 王华忠1， 钱建良2

1 同济大学海洋与地球科学学院， 上海200092 2 Department of Mathematics, Michigan State University, East Lansing, MI, 48824， U.S.A.

角度域成像道集是叠前深度偏移的重要输出结果，它是偏移速度分析、各向异性分析和AVA分析等研究工作的基础.目前存在的角度域成像道集的生成方法受计算效率或角度分辨率的影响，仍然满足不了实际生产的要求.角度域成像道集的生成方法可以大致分为直接法和间接法两大类.在直接方法中，波矢量方向计算和局部平面波分解是两个最重要的内容，它们共同决定角度域成像条件的实现效率和角度域成像道集的质量.为了完善现有的角度域成像道集生成方法，本文提出一种新的波矢量计算方法和局部平面波分解方法.本文先用波动方程任意时刻的柯西条件构造一个只含原波场负频率成分的柯西波场，然后根据柯西波场在时间波数域的振幅谱计算波场的波矢量方向.该方法仅在需要计算波矢量方向的时刻合成柯西波场，不需要增加额外的数据读写操作，是一种高效的波矢量计算方法.本文还以柯西波场为基础提出一种高效的局部平面波分解方法，保证角度域成像条件的实现效率.结合柯西波场和局部平面波分解方法，本文最后给出一种新的角度域成像方法.文中最后的数值实验证明该方法得到的角度域成像道集具有理想的角度分辨率，可以反映地下构造的角度照明情况.

柯西波场； 波矢量方向； 局部平面波分解； 角度域成像条件； 角度域成像道集

1　引言

共成像点道集(Common-Imaging Gathers, CIGs)不仅反映偏移速度、正演算法和散射模型的准确性，还反映地下介质的各向异性、AVO和AVA等特征，是偏移技术从构造成像走向波阻抗反演的一个关键环节.受波场多传播路径的影响，叠前深度偏移方法产生的CIGs按照地表记录属性(如共方位角、共偏移距等)进行排序很容易出现偏移假象，降低CIGs的质量.而利用角度域成像条件得到的角度域共成像点道集(Angle-Domain Common-Imaging Gathers, ADCIGs)不受波场多传播路径的影响，是信噪比比较高的CIGs(Xu et al., 2001).

基于ADCIGs的速度分析方法和AVA分析方法已经存在丰富的理论和应用成果(Jin and Madariaga, 1993; Symes, 1993; Beydoun et al., 1993; Liu and Bleistein, 1995; Tura et al.,1998; Dong and Ponton, 1999; Chauris et al., 2002).角度域成像条件大致可以分为以下几种方法：局部平面波分解方法(Xie and Wu, 2002; Yan and Xie, 2009; Xu et al., 2011)；Poynting矢量方法(Yoon and Marfurt, 2006;Vyaset al., 2011a, 2011b; Dickens and Winbow, 2011; Yoon et al., 2011)；扩展成像道集转化方法(Sava and Fomel, 2003, 2006; Fomel, 2004).其中局部平面波分解方法和Poynting矢量方法通过计算波矢量方向将偏移结果投影到对应的ADCIGs上，属于直接方法.扩展成像道集转化方法先用扩展成像条件得到扩展成像道集，然后再用线性Radon变换将扩展成像道集转化为ADCIGs，属于间接方法.

局部平面波分解方法在频率空间域中对波场进行平面波分解并生成角度域成像道集，包括单程波方法和双程波方法.其中单程波方法在频率空间域中沿深度方向对波场进行外推并对波场做局部平面波分解，利用得到的局部平面波进行成像并提取ADCIGs(Xie and Wu, 2002; Soubaras, 2003; Yan and Xie, 2009).双程波方法通过记录所有时刻的正演波场，然后用Fourier变换将时间空间域的波场变换为频率空间域的波场，在频率空间域中实现局部平面波分解和角度域成像条件并提取ADCIGs(Xu et al., 2011).在双程波方法中，为了实现时间域的Fourier变换，需要记录所有时刻的震源波场和检波点波场并以时间为快维对波场重新排序.该重排序操作会增加大量的数据读写操作.

Poynting矢量方法在波场的正演过程中先用Poynting矢量计算波矢量方向，然后用波矢量方向判断波场的传播方向并实现角度域成像条件，生成ADCIGs.Poynting矢量是目前常用的波矢量计算方法，利用Poynting矢量可以快速得到波矢量方向并实现角度域成像条件(Yoon and Marfurt, 2006;Vyaset al., 2011a, 2011b; Dickens and Winbow, 2011; Yoon et al., 2011).但Poynting矢量得到的是波场的能流方向.当波场的波前存在交叉情况时，Poynting矢量无法区分沿不同方向传播的波前，得到的矢量方向实际上是所有波前叠加后的能流方向.因此，Poynting矢量方法在计算复杂波场的波矢量方向时具有一定的局限性.另外，Poynting矢量方法是一种抗噪能力比较差的波矢量方向计算方法.在计算检波点波场的波矢量方向时，Poynting矢量方法会因为地震数据中的噪音而得不到准确的波矢量方向.

扩展成像道集转化方法先用扩展成像条件得到局部偏移距道集或时延相关道集，然后再用线性Radon变换将该扩展成像道集转化为ADCIGs(Sava and Fomel, 2003, 2006; Fomel, 2004).在二维情况下，时延相关道集和局部偏移距道集都可以较为简单地转化为ADCIGs.但在三维情况下，时延相关道集是一个四维数据体(空间的x、y、z坐标和时延量τ)，在转化为五维的ADCIGs(空间的x、y、z坐标，入射角φ和方位角φ)时会损失部分信息.而局部偏移距道集与ADCIGs之间的关系在三维情况下会变得很复杂，大大增加转化过程所需要的计算量，降低该方法的计算效率.另外，将扩展成像道集从局部偏移距域(或时延量)投影到角度域时还存在采样不足、采样不均匀等现象，这些现象都会降低ADCIGs的质量，影响正确的AVA特征.

在上述的几种方法中，局部平面波分解方法和Poynting矢量方法的物理意义更明确，得到的ADCIGs的物理含义更清楚.这两种方法均需要计算局部波场的传播方向.其中局部平面波分解方法得到的波矢量方向信息更丰富、更准确，但需要更多的计算量和数据读写操作；Poynting矢量方法利用Poynting矢量计算波矢量方向，计算效率高，但抗噪能力和适用范围不及局部平面波分解方法.高效、准确的波矢量方向计算方法和局部平面波分解方法是角度域成像条件两个最重要的研究内容，它们共同决定角度域成像条件的计算效率和ADCIGs的质量.任何一种不准确的波矢量方向计算方法或局部平面波分解方法均会在ADCIGs中产生偏移假象，影响ADCIGs的质量.签于上述原因，本文提出一种用波动方程的柯西条件构造柯西波场的方法，并以柯西波场为基础提出一种高效的波矢量方向计算方法和局部平面波分解方法，提出一种新的角度域成像方法.最后的数值实验证明本文提出的角度域成像方法可以得到高质量的ADCIGs，为下一步的波阻抗反演和AVA分析等研究工作提供了良好的基础.

2　理论与方法

2.1柯西波场的构造方法

考虑如下波动方程：

(1)

该方程为波动方程的纯初值问题，描述了声波介质在无外力情况下的波传播现象.由于方程(1)的解取决于方程在t0时刻的Cauchy条件f(x)和g(x)，因此利用f(x)和g(x)可以得到u(x,t)在t0时刻的传播方向.假设u(x,t)存在如下的形式表达式：

(2)

其中A和φ分别表示u(x,t)的振幅函数和相位函数.沿时间方向对该形式表达式进行求导可得：

(3)

其中At和φt分别表示A和φ对时间的导数，ω表示波场的频率.在上式推导中，本文假设u(x,t)满足高频假设条件，因此有At≪Aφt.根据波动方程的频散关系可以得到u(x,t)的频率和波数具有如下关系：

(4)

其中v(x)表示x处的声波速度，k表示波场的波数.结合(3)和(4)得到的结果，用u(x,t)和ut(x,t)可以构造如下的复值波场c(x,t)：

(5)

上面介绍的波矢量方向计算方法用波动方程的柯西条件构造一个只含原波场负频率成分的柯西波场c(x,t)，并用c(x,t)在空间域的Fourier变换结果计算波矢量方向.CWF的构造方法利用快速Fourier变换提高其计算效率，并保留了波场所有方向的波矢量信息，可以高效、准确地计算出波场的波矢量方向.

2.2局部平面波分解方法

β(k)=arctan(kz/kx),

(7)

(8)

集合W给出了波场u(x,t)包含的所有平面波传播角度的集合.

(9)

公式(9)即为本文在实现角度域成像条件时所采用的局部平面波分解方法.该方法最大的优点是计算复杂度小、计算成本低，保证实现角度域成像条件的计算效率.利用集合W和平面波分解方法(9)，可以将原波场u(x,t)分解为一系列沿不同角度传播的平面波之和：

(10)

2.3角度域成像条件

波矢量方向的计算和波场的局部平面波分解是实现角度域成像条件的两个关键环节.这两个环节即决定实现角度域成像条件的计算效率，也决定最终生成的ADCIGs的质量.利用CWF本文可以高效地得到波矢量方向、实现波场的平面波分解，总波场则可以表示为这些平面波的叠加结果.将波场的平面波表达形式(10)代入到单炮叠前偏移公式可得I(x,xs)=∫us(x,xs,t)ur(x,xs,t)dt

(11)

Rloc(x,xs,θ,t)=

Iloc(x,xs,θ)=∫Rloc(x,xs,θ,t)dt,

(13)

其中θ表示入射波与反射波之间的张角，Rloc(x,xs,θ,t)表示t时刻在张角θ上的成像结果，Iloc(x,xs,θ)表示第xs炮在局部窗内得到的ADCIGs.由于计算Rloc(x,xs,θ,t)需要对Ws(t)和Wr(t)内的所有角度进行循环，因此波场越复杂、局部平面波的个数越多，Rloc(x,xs,θ,t)的计算成本越高.

公式(13)给出了角度域成像条件的计算方法，利用该公式可以得到波场在任意局部空间窗内的ADCIGs，将所有局部空间窗内的成像结果叠加后即可得到全空间的ADCIGs.图1展示了角度域成像条件的实现流程，在该实现流程中虚框内的计算过程是角度域成像条件与常规零延迟相关成像条件的主要区别.

图1　基于炮道集的角度域成像条件Fig.1　The angle domain imaging condition of shot domain migration

3　数值实验

用本文提出的角度域成像条件分别计算水平层状模型和Marmousi模型的ADCIGs.这两个数值实验证明了本文提出的角度域成像方法可以准确地获得复杂模型的ADCIGs，是一种适用性比较强的成像方法.

图2　角度域成像条件实现流程示意图

(a、b) 某一时刻炮点端和检波点端的波场快照； (c、d) 该时刻图a、b红框内的波场以及波场对时间的一阶导数； (e、f) 由图c、d构造得到的柯西波场所对应的振幅谱； (g、h) 振幅谱e、f在角度域的能量分布函数，横坐标表示角度，单位°， 纵坐标表示能量大小； (i、j) 局部平面波分解结果； (k) 局部平面波的成像结果.

Fig.2The realization flow of angle domain imaging condition

(a, b) Snapshots of source wave field and receiver wave field； (c, d) Local wave field and its derivative with respect to time in the red rectangles of figure 2a and 2b； (e, f) Spectrum of Cauchy Wave Field constructed from Cauchy condition shown in figure 2c and 2d； (g, h) Energy distribution function of Cauchy Wave Field in angle domain, the axises mean energy and angle respectively； (i, j) Local plane wave decomposition result； (k) Imaging result of local plane waves shown in figure 2i and 2j.

3.1水平层状模型

首先用一个水平层状模型测试本文提出的角度域成像方法.图3展示了一个水平层状速度模型，该模型的观测系统如下：炮点间隔100 m，检波点间隔20 m；炮点范围从地表x=0 km到x=10 km；每一炮的接收排列均从地表x=0 km到x=10 km.图4为xs=6 km处的单炮数据得到的不同地表位置的ADCIGs，其中横坐标的快维表示ADCIGs对应的反射角，慢维表示ADCIGs对应的地表位置.在图4中，位于炮点正下方的ADCIGs的能量集中在0°的位置，表示该炮对炮点正下方的照明能量主要集中在入射角为0°的位置.随着地表偏移距的增加，该炮对地下层位照明的入射角越来越大.图4所示的ADCIGs准确地反映了该炮对地下层位的角度照明情况，为后续AVA分析和角度照明分析等研究工作提供了重要基础.

图5展示了用所有炮数据生成的ADCIGs.在图5中，浅层界面的角度照明比较均衡，除了边界处的ADCIGs受观测系统的影响外，其他地方的ADCIGs基本达到全角度照明.在生成的ADCIGs中，小角度的照明能量对应着地震波的反射现象，是成像中的有效能量；大角度的照明能量对应着地震波的全反射现象，是成像中的噪音部分，这部分能量会降低最后成像结果的分辨率.受地表观测孔径和速度模型的影响，模型深层界面的照明角度主要集中在小反射角范围内，数据对界面的照明情况基本表现为：界面的深度越深，有效照明角度的范围越小，照明角度的采样间隔越小.

图3　水平层状速度模型Fig.3　The layer velocity model

图4　单炮数据的ADCIGs，炮点位置位于xs=6 km处Fig.4　The ADCIGs of single shot gather, the position of source is xs=6 km

图5　所有炮数据得到的ADCIGsFig.5　The ADCIGs generated by all shot gathers

图6　不同入射角范围的ADCIGs的叠加结果(a) 入射角为0°～10°； (b) 入射角为10°～20°； (c) 入射角为20°～30°； (d) 入射角为30°～40°； (e) 入射角为40°～50°； (f) 入射角为0°～50°.Fig.6　The stack results of ADCIGs with different incident angle (IA)(a) From 0 to 10 degree； (b) From 10 to 20 degree； (c) From 20 to 30 degree； (d) From 30 to 40 degree； (e) From 40 to 50 degree；(f) From 0 to 50 degree.

图7　不同入射角范围的ADCIGs的叠加结果(a) 入射角为0°～10°； (b) 入射角为10°～20°； (c) 入射角为20°～30°； (d) 入射角为30°～40°.Fig.7　The stack results of ADCIGs with different incident angle (IA)(a) From 0 to 10 degree； (b) From 10 to 20 degree； (c) From 20 to 30 degree； (d) From 30 to 40 degree.

图8　入射角为0°～45°的ADCIGs的叠加结果Fig.8　The stack result of ADCIGs with IA from 0 to 45 degree

图9　不同地表位置的ADCIGs(最大入射角为60°)Fig.9　The ADCIGs of different place (the maximum incident angle is 60 degree)

图6展示的是不同入射角范围的ADCIGs的叠加结果.图6所示的偏移结果与地下构造的实际角度照明情况一致(如图6(a—e)中第一个反射界面的照明情况所示)，进一步证明本文提出的角度域成像方法具有较高的角度分辨率，可以准确地得到地下界面的角度照明情况，生成高质量的ADCIGs.成像中的有效信号(入射波与反射波的相关结果)主要集中在ADCIGs的中、小入射角范围内，ADCIGs的大入射角的能量多为炮点的反射波、折射波与检波点的反射波相关后的结果，是成像中的低频噪音部分.因此，用ADCIGs的中、小入射角范围内的数据进行叠加即可得到最终的偏移结果.图6f展示了入射角范围为0～50°的ADCIGs的叠加结果，该图即为最终的偏移结果.

3.2Marmousi模型

用本文的角度域成像方法计算Marmousi模型的ADCIGs.为了改善地下构造的角度照明情况，这里设计了如下观测系统：炮点间隔100 m，检波点间隔25 m；炮点范围从地表x=3.6 km到x=9.0 km；每一炮的接收排列均从地表x=3.0 km到x=9.0 km.根据前面分析得到的结论可知，ADCIGs的中、小入射角范围内的能量是成像中的有效信号成分，大入射角的能量是成像中的低频噪音成分.因此，在成像过程中本文只生成中、小入射角范围的ADCIGs.图7(a—d)分别展示了入射角范围为0～10°、10～20°、20～30°和30～40°的ADCIGs的叠加结果.图8是所有入射角(0～45°)的ADCIGs的叠加结果，图9是不同地表位置的ADCIGs，ADCIGs的入射角范围为0～60°.图7、图8和图9证明了本文提出的角度域成像条件可以准确地获得地下构造的角度照明情况，生成高质量的ADCIGs.

4　结论与讨论

波矢量方向计算和局部平面波分解是角度域成像条件的两个核心内容.本文提出一种利用波动方程的柯西条件构造柯西波场(CWF)的方法，并以CWF为基础实现波矢量方向计算和局部平面波分解，从而实现角度域成像条件并生成ADCIGs.CWF是根据波动方程的柯西条件构造的只含原波场负频率成分的复值波场，因此可以根据CWF在波数域的振幅谱计算波矢量方向.另外，本文还提出一种以CWF为基础的局部平面波快速分解方法.该方法根据CWF在角度域的能量分布情况设计相应的角度滤波器，并用该滤波器实现局部平面波的快速分解，保证实现角度域成像条件的计算效率.根据测不准原理，局部平面波的角度分辨率主要受局部空间窗的孔径和波前的曲率影响，而局部平面波的角度分辨率将决定生成的ADCIGs的角度分辨率.在本文提出的角度域成像方法中，CWF是最核心的内容，它是本文计算波矢量方向和实现局部平面波分解的基础.本文提出的CWF构造方法以快速Fourier变换为基础，具有较高的计算效率，保证了该方法的实用性.本文相关的数值实验证明本文提出的角度域成像条件可以生成高质量的ADCIGs.

Beydoun W, Hanitzsch C, Jin S. 1993.Why migrate before AVO? A simple example.∥55th EAEG Meeting.EAGE.

Chauris H, Noble M S, Lambaré G, et al. 2002.Migration velocity analysis from locally coherent events in 2-D laterally heterogeneous media, Part I: Theoretical aspects.Geophysics, 67(4): 1202-1212.

Dickens T A, Winbow G A. 2011.RTM angle gathers using Poynting vectors.∥81st SEG Annual International Meeting. SEG,3109-3113.

Dong W J, Ponton M. 1999.AVO inversion and interpretation via localized 3D migrations.∥6th International Congress of the Brazilian Geophysical Society. SEG, 816-819.

Fomel S. 2004. Theory of 3-D angle gathers in wave-equation imaging.∥ 74th SEG Annual International Meeting. SEG,1053-1056.

Jin S, Madariaga R. 1993. Background velocity inversion with a genetic algorithm.GeophysicalResearchLetters, 20(2): 93-96. Liu Z Y, Bleistein N. 1995. Migration velocity analysis:Theory and an iterative algorithm.Geophysics, 60(1): 142-153.

Sava P C, Fomel S. 2003.Angle-domain common-image gathers by wavefield continuation methods.Geophysics, 68(3): 1065-1074.

Sava P, Fomel S. 2006.Time-shift imaging condition in seismic migration.Geophysics, 71(6): S209-S217.

Soubaras R. 2003. Angle gathers for shot-record migration by local harmonic decomposition.∥SEG Technical Program Expanded Abstracts. SEG, 889-892.

Symes W W. 1993.A differential semblance criterion for inversion of multioffset seismic reflection data.JournalofGeophysicalResearch:SolidEarth, 98(B2): 2061-2073.

Tura A, Hanitzsch C, Calandra H. 1998.3-D AVO migration/inversion of field data.TheLeadingEdge, 17(11): 1578-1578.

Vyas M, Du X, Mobley E, et al. 2011a. Methods for computing angle gathers using RTM.∥73rd EAGE Conference & Exhibition. EAGE.

Vyas M, Nichols D, Mobley E. 2011b. Efficient RTM angle gathers using source directions.∥81st SEG Annual International Meeting. SEG, 3104-3108.

Xie X B, Wu R S. 2002. Extracting angle domain information from migrated wavefield.∥72nd SEG Annual International Meeting.SEG, 1360-1363.

Xu S, Chauris H, Lambaré G, et al. 2001.Common-angle migration: A strategy for imaging complex media.Geophysics, 66(6): 1877-1894.

Xu S, Zhang Y, Tang B. 2011.3D angle gathers from reverse time migration.Geophysics, 76(2): S77-S92.

Yan R, Xie X B. 2009. A new angle-domain imaging condition for prestack reverse-time migration.∥79th SEG Annual International Meeting. SEG, 2784-2788.Yoon K, Marfurt K J. 2006.Reverse-time migration using the Poyntingvector.ExplorationGeophysics, 37(1): 102-107.

Yoon K, Guo M H, Cai J, et al. 2011. 3D RTM angle gathers from source wave propagation direction and dip of reflector.∥81st SEG Annual International Meeting. SEG, 3136-3140.

(本文编辑胡素芳)

Cauchy wave field and its application in angle-domain imaging condition

WANG Xiong-Wen1， WANG Hua-Zhong1， QIAN Jian-Liang2

1SchoolofOceanandEarthScience,TongjiUniversity,Shanghai200092，China2DepartmentofMathematics,MichiganStateUniversity,EastLansing,MI, 48824，U.S.A.

Angle Domain Common-Imaging Gathers (ADCIGs) is an important output of pre-stack depth migration. It is the basis of migration velocity analysis, anisotropy analysis and AVA analysis. However, there are still a lot of problems that reduce the efficiency of the existing ADCIGs generating methods and hinder the application of ADCIGs.Methods used to generate ADCIGs are mainly divided into two types: direct and indirect methods. For the direct methods, wave vector computing method and local plane wave decomposing method are the most important for realizing angle domain imaging condition.They decide the efficiency of imaging condition and the quality of ADCIGs. In order to perfect the exiting methods, we propose a new wave vector computing method and an efficient local plane wave decomposing method in this paper. The new method proposed in this paper utilizes the Cauchy condition of wave equation at any given time to generate a complex-valued wave field called Cauchy Wave Filed (CWF), which contains the negative frequencies of original wave fields only. Since CWF contains negative frequencies only, the energy distribution of CWF in wave number domain reveals the propagating direction of original wave fields. The only added cost of new method is the generation of CWF when imaging condition is applied. Therefore, it avoids additional In Out cost and is an efficient method to computing the wave vector. In addition, an efficient local plane wave decomposing method is also proposed in this paper to reduce the computation cost of ADCIGs and make the new method applicable. With CWF and local plane wave, we give an efficient angle-domain imaging condition. The numerical examples given in this paper prove that the ADCIGs obtained by the new proposed method has an ideal angle resolution, and have the ability to reveal the angle illumination of subsurface.KeywordsCauchy Wave Field (CWF); Polarization of wave vector; Local plane wave decomposition; Angle-domain imaging condition; Angle-domain common-imaging gathers

10.6038/cjg20161024.

国家“973”重点基础研究发展计划项目(2011 CB201002)，国家自然科学基金(41374117)以及国家重大专项(2011ZX05005-005-008HZ, 2011ZX05006-002, 2011ZX05023)联合资助.

王雄文，男，1983年出生，同济大学海洋与地球科学学院博士毕业生. 主要从事地震信号处理、叠前深度偏移和速度分析等研究. E-mail: xiongwenwang@gmail.com

10.6038/cjg20161024

P631

2015-07-03，2016-01-28收修定稿

王雄文， 王华忠， 钱建良. 2016. 柯西波场及在角度域成像条件中的应用. 地球物理学报，59(10):3798-3809,

Wang X W， Wang H Z， Qian J L. 2016. Cauchy wave field and its application in angle-domain imaging condition.ChineseJ.Geophys. (in Chinese),59(10):3798-3809,doi:10.6038/cjg20161024.