地球物理学报  2021, Vol. 64 Issue (4): 1375-1388   PDF    
引用本文
吴成梁, 王华忠, 冯波, 等. 2021. 基于CLG光学流和波场分解的逆时偏移角度道集提取方法研究. 地球物理学报, 64(4): 1375-1388, doi: 10.6038/cjg2021O0088.  
Wu C L, Wang H Z, Feng B, et al. 2021. RTM angle gathers based on the Combining Local and Global (CLG) optical flow method and wavefield decomposition method. Chinese J. Geophys. (in Chinese), 64(4): 1375-1388, doi: 10.6038/cjg2021O0088.  
基于CLG光学流和波场分解的逆时偏移角度道集提取方法研究
吴成梁, 王华忠, 冯波, 盛燊     
同济大学海洋与地球科学学院波现象与智能反演成像研究组, 上海 200092
收稿日期 2020-03-05, 2021-01-09 收修定稿
基金项目: 国家重点研发计划（2018YFA0702503，2019YFC0312004），国家自然科学基金（42074143，41774126），上海市浦江人才计划、中国石化地球物理重点实验室项目（33550006-20-ZC0699-0011，33550006-19-FW0399-0041）和南方海洋科学与工程广东省实验室（湛江）资助项目（ZJW-2019-04）资助
第一作者简介: 吴成梁, 男, 1990年生, 现就读于同济大学海洋与地球科学学院, 博士在读.主要研究方向为地震波偏移反演成像方法与速度建模方面的研究工作.E-mail: wuchengliang1990@163.com
摘要：角度域共成像点道集是衔接叠前地震数据与储层特征的重要桥梁，对地震偏移成像与储层描述具有重要意义.与克希霍夫偏移和单向波动方程偏移相比，逆时偏移是复杂地区最精确的成像方法.高效稳健地生成逆时偏移角度道集目前仍然是一个挑战.本文主要讨论如何采用光学流方法高效、高质量地提取角度道集.在逆时偏移波场外推过程中，光学流方法可以估计波场传播方向.其中Lucas-Kanade（LK）和Horn-Schunck（HS）方法是光学流方法中两种典型的方法.LK光学流方法是一种局部方法，该方法依赖于局部点的梯度值，但是容易出现奇异现象，HS光学流方法属于全局方法，波场方向估计依赖于整个波场，易受噪声影响，对异常值比较敏感，导致整体波场方向计算精度不高.本文提出采用局部和整体结合（Combining Local and Global，CLG）的光学流方法估计波场传播方向.该方法可以有效地提高波场方向的精度，并且简单高效，便于并行处理.对比HS光学流方法，CLG光学流方法几乎不增加额外的计算量.另外，为了弱化光学流方法无法处理波前重叠问题，本文利用解析波场和方向滤波对波场进行方向分解，仅需波场的空间傅里叶变换即可实现任意波场方向分解，将分解后的波场分别估计波场反向，提取成像结果.进一步地，在估计反射张角和方位角时，本文提出有效的归一化方法和改进的最小二乘除法，提高角度估计的精度和稳定性.最后，理论和实际资料例证了本文提出方法的有效性.
关键词: 逆时偏移      角度道集      光学流      波场分解      波场方向      角度估计     
RTM angle gathers based on the Combining Local and Global (CLG) optical flow method and wavefield decomposition method
WU ChengLiang, WANG HuaZhong, FENG Bo, SHENG Shen     
Wave Phenomena and Intelligent Inversion Imaging Group(WPI), School of Ocean and Earth Science, Tongji University, Shanghai 200092, China
Abstract: Angle-Domain Common-Image Gathers (ADCIGs) is an important bridge connecting the pre-stack seismic data and reservoir characteristics, which is of great significance for seismic migration and reservoir characterization. Compared with Kirchhoff migration and one-way wave-equation migration, Reverse-Time Migration (RTM) is the most accurate imaging method for complex areas. Generating the RTM angle gathers with efficiency and robustness still is a challenge. In the paper, we will mainly discuss how to efficiently generate high-quality angle gathers with the optical flow method. The optical flow method can be applied for the wavefield direction estimation, during the process of RTM wavefield extrapolation. The Lucas-Kanade (LK) technique and the Horn-Schunck (HS) algorithm are two fundamental approaches. However, the LK method relies on local wavefield gradient and suffers from the singular phenomenon. The HS method is sensitive to noise, and the accuracy of the estimated wavefield direction is not high. We propose using the Combining Local and Global (CLG) optical flow method to estimate the wavefield direction. The accuracy of wavefield direction can be improved effectively by the proposed CLG optical flow method, especially for the reflected interface. And the proposed method is simple, efficient and convenient for parallel processing. Moreover, compared with the HS optical flow method, the proposed method has not much extra calculation. Then, we apply the wavefield decomposition method with the analytical wavefield and directional filtering to weaken the inaccurate estimation of optical flow method at the wave-front overlapping situation. In the analytical wavefield, the directional decomposition can be realized by just only the spatial Fourier transform of the wavefield. Then we calculate wavefield direction and extract angle gathers based on the decomposed wavefield. Furthermore, an effective normalization method and an improved least square division method are proposed to improve the accuracy and stability of angle estimation. Finally, theoretical and practical data illustrate the effectiveness of the proposed method.
Keywords: Reverse-time migration    Angle gathers    Optical flow    Wavefield decomposition    Wavefield direction    Angle estimation    
0 引言

角度域共成像点道集(ADCIGs)包含背景速度和地下角度反射系数的信息，对于速度建模和储层描述具有重要意义.与克希霍夫偏移、Beam偏移和单向波偏移相比，逆时偏移(Baysal et al., 1983；Whitmore，1983)能够精确地处理横向变速情况以及多波至现象，克服陡倾角限制，成为了复杂介质的首选成像方法.逆时偏移角度道集可以采用扩展成像条件方法(Sava and Fomel，2003，2006；Fomel，2004；Duveneck，2013)或局部平面波分解方法(Xie and Wu, 2002；Xu et al., 2011；Yan and Xie, 2011，Tang et al., 2013；Hu et al., 2015)实现.扩展成像条件方法首先生成地下局部偏移距道集，然后借用倾斜叠加等方法转换为地下角度道集.局部平面波分解方法需要将波场分解为局部平面波分量，应用角度域成像条件提取角度道集，计算代价较大.在逆时偏移角度道集中，算法的高效性和并行性是制约高效提取角度道集的重要因素.波矢量方向估计方法(Yoon et al., 2004；Yoon and Marfurt, 2006；Dickens and Winbow, 2011；Vyas et al., 2011a, b；Zhang and McMechan，2011a，b；王保利等，2013；Jin et al., 2014；Yoon，2017)在逆时偏移角度道集生成方法中是一类比较高效的方法.

在波矢量方向估计方法中，可以通过任意空间和任意时刻的振幅和相位梯度来计算传播方向.一旦估计出震源波场和检波点波场的方向矢量，就可以计算出反射张角和方位角，并利用互相关成像条件提取成像道集.在波矢量估计方法中，可以采用坡印廷矢量方法计算波场传播方向.由于坡印廷矢量方法计算高效，具有较高的角度分辨率，在逆时偏移角度道集生成方法中被广泛地应用.但是由于在检波点端波场比较复杂，此时计算的坡印廷矢量通常会出现不稳定情况，Yoon等(2011)和Vyas等(2011a)提出利用坡印廷矢量计算震源端波场传播方向和地下反射界面倾角来估计反射角；Zhao等(2012)提出通过偏移剖面预测反射层法向向量，利用稳定的震源端方向波场结合反射层法向向量来计算角度道集.虽然坡印廷矢量方法比较高效，但是该方法无法处理波前重叠问题(Patrikeeva and Sava, 2013).

在生成逆时偏移角度道集中，外推波场的方向也可由光学流方法计算.光学流最先由Gibson于1950年提出，被广泛应用于解决连续图像帧之间的视运动问题.目前有许多种计算光学流场的方法(Sobey and Srinivasan, 1991；Barron et al., 1994；Brox et al., 2004；Bruhn et al., 2005；Szeliski，2010).其中Lucas-Kanade方法(LK，或称局部方法)(Lucas and Kanade, 1981)和Horn-Schunck方法(HS方法，或称全局方法)(Horn and Schunck, 1981)是两类主要的方法，其他方法大都是基于上述方法的补充.Vyas等(2011a)应用LK光学流方法生成角度道集，Zhang(2014)采用HS光学流方法估计波场方向，提取角度道集.然而，LK光学流方法容易出现奇异现象，依赖于局部点的梯度值，而HS光学流方法依赖于整体的能量约束，距离较远的波场也会参与到波场方向计算中，导致计算的波场方向精度较低.

本文提出采用CLG(Combining Local and Global)光学流方法计算波场方向，提高波场方向的精度.CLG光学流方法结合局部LK方法和全局HS方法的优点，在施加全局约束的基础上，考虑局部点的加权作用，在保持波场方向估计稳定的基础上，提高估计的波场方向精度.另外，为了弱化光学流方法无法处理波前重叠问题，本文提出利用解析波场对波场进行方向分解，将分解后的行波分别成像，提取角度道集.进一步地，在提取角度道集过程中，需要将光学流场计算的波场方向转换为反射张角和方位角，本文提出一种有效的归一化方法和改进的最小二乘除法，提高角度计算的准确性，避免带限子波引起的角度计算的不稳定.最后，采用角度道集面元化和规整化来优化角度道集.

1 方法理论 1.1 波场方向估计

光学流问题的基本假设可以描述为后续时刻外推的波场值不随时间变化：

(1)

其中(x, y, z)是地下空间坐标，t代表时间，u是波场.基于小位移假设，采用一阶泰勒展开可得到光学流方程.

(2)

其中是空间导数，ut=是时间导数，vx，vy和vz分别是x，y方向和z方向的光学流场分量.将方程(2)写成矩阵形式：

(3)

其中x=(x, y, z)，v=(vx, vy, vz)，ux=(ux, uy, uz).由于在方程(3)中未知量的个数大于方程的个数，该问题是一个欠定问题，需要引入额外的约束求解该问题.

在局部LK光学流方法中，通过假设未知的光学流矢量在局部邻域Ω内是恒定不变的，构建如下加权最小二乘泛函求解上述光流问题：

(4)

其中WΩ为在局部邻域Ω内的加权系数，最小化ELK(v)得到如下线性方程组：

(5)

在每个外推时刻逐个空间点求解上述3×3的线性方程组即可得到光学流场.LK光学流方法计算量小，只受局部范围的影响，距离当前点较远区域的误差不会影响当前点的计算，误差不具有传播性.但是LK方法不稳定，当梯度消失时，线性方程组会出现奇异现象.因此，LK光学流方法通常在稀疏场中实现，然后采用不同的插值方法插密.此外，光学流的精度依赖于方程组中梯度张量的准确性.

在全局HS光学流方法中，通过引入额外的全局能量的信息约束，求解上述欠定问题，构建如下的误差泛函：

(6)

其中▽为梯度算子，C(v, u)是正则化项，λ为加权因子，在本文的测试中，λ取值为1.正则化项用于惩罚光学流场中较大的空间导数，保证在波场中相邻的点具有相似的方向.最小化误差泛函(式(6))，可转化为求解相应的Euler-Lagrange方程：

(7)

其中Δ为拉普拉斯算子.全局HS光学流方法通过引入整体的能量约束，避免出现奇异现象，即使在某些点上没有梯度值，但是由于全局能量约束，也能比较稳定地求解光学流场，因此该方法能产生较密的光学流场.然而波场的传播具有局部方向性，整体的能量约束带来的问题是，距离计算点较远区域的波场也会参与到当前点的计算中，导致计算的波场方向精度有所下降.

结合局部邻域变化特征和全局的能量约束，本文提出构建如下的CLG光学流误差泛函：

(8)

其中ω为光学流场构成的向量，‖▽ω‖2是光学流梯度场的正则化项，▽4u为波场导数构成的向量，为局部邻域内的加权系数，为局部邻域内构建的结构张量矩阵.在保持全局能量约束的情况下，通过采用对局部邻域的加权，强调波场的局部方向连续性，避免较远区域的影响，在保持波场方向估计稳健情况下，提高波场方向计算的精度.采用变分原理，公式(8)对应的Euler-Lagrange方程为(具体推导过程见附录A)：

(9)

CLG光学流方法不仅克服了梯度为零时无法估计光学流问题，而且实现了局部约束以提高抗噪性.采用CLG光学流方法估计的波场方向更加精确和稳定.对比HS光学流方法(方程(7))，CLG光学流方法增加的计算量并不大.方程(9)和方程(7)的迭代格式是一致的，比较容易实施并行化处理，方便融合在逆时偏移的波场外推过程中.

1.2 波场方向分解

采用CLG光学流方法可以计算波场方向.但是光学流方法的假设前提是在每个成像点每个成像时刻只有一个传播方向.因此，该方法不能处理波前重叠的情况.本文采用解析波场和空间傅里叶变换对外推的波场进行行波方向分解，将分解后的行波分别计算光学流场，提取波场传播方向构建角度道集.

在解析波场中，仅包含正频率信息，仅需对波场实施空间域的傅里叶变换即可实现波场分解，从而可以避免时间维的傅里叶变换.解析波场是复数波场，其实部可通过波动方程外推获得，其虚部可通过对地下波场进行Hilbert变换得到.但是由于Hilbert变换是关于时间方向的卷积，理论上需要将整个波场存盘并进行卷积，对内存和存储要求比较大.借助于波动方程的源项和波场具有线性关系(Liu et al., 2006)，把对波场的Hilbert变换转化为震源项的Hilbert变换，利用新的“震源”进行传播，从而构建出解析波场的虚部(Hu et al., 2015).以震源端为例，震源端的解析波场如下所示：

(10)

(11)

其中c(x)为偏移速度场，δ(·)克罗内克函数，fH(t)是震源子波f(t)的希尔伯特变换，usH是源端波场us的希尔伯特变换.是对应的解析波场.对解析波场实施空间域傅里叶变换，变换到波数域：

(12)

其中F算子代表傅里叶变换，ks为波数.因此，可以根据波数关系提取任意方向、任意角度宽度传播的波场.

(13)

其中k1β, k2β分别为波数(角度)边界矢量，▽k=(k2β－k1β)为对应的波数(角度)宽度矢量，▽α为阈值系数.iF代表逆傅里叶变换.k1β, k2β可根据地下波场的复杂程度而变化，波场越简单，角度宽度可取得大一些，若波场交叉比较严重，则角度宽度需要取小一些，但是对应增加的计算量较多.

虽然理论上可以提取任意角度的波场.但是在逆时偏移过程中，波场方向分解需要付出较多的计算代价.为了平衡计算代价和方向估计精度，我们可以根据地下介质的特征进行适应性的选择.若把地下介质抽象为在空间上广泛分布的缓变的层状沉积层，加上火山活动，构造运动等引起的大、小尺度的速度异常体(或波阻抗变化的异常体)(王华忠等，2015).则只需要采用上下行波分解即可满足大部分的成像需要，能够在一定程度上弱化波前交叉问题.此时行波分解可以简单的写成下式：

(14)

(15)

其中为震源端解析波场的傅里叶变换结果，usup和usdown分别是震源端分解后的上下行波.同样地，检波点端波场采用相同的方式实现波场方向分解.当地下介质中存在陡倾角构造时，此时可以把行波分解修改为左右行波分解，减弱波场交叉对后续CLG光学流方法估计方向的影响.因此可以根据地下介质的波场传播特征，选择合适的波场分解方向，尽量避免波场交叉的影响.

1.3 反射张角和方位角估计

光学流方法计算的波场方向矢量代表了波场的运动差异.在估计反射张角和方位角之前，需要实施归一化处理.本文在归一化过程中通过采用对光学流和波场同时约束，避免噪音干扰和弱值的影响.

(16)

其中Ps是震源端的光学流矢量场，是归一化之后的方向矢量.epsps和epsus分别是阈值系数，Aps(t)和Aus(t)分别是当前时刻光学流场和波场的最大振幅.

反射张角θ和方位角φ一般可通过下式计算：

(17)

(18)

其中nx=(1, 0, 0)和nz=(0, 0, 1)分别为单位矢量.但是由于地震波是振荡的.在计算反射张角和方位角时，涉及到除法，分母较小的位置会使得除法变得非常不稳定.为了提高角度估计的稳定性和精度，以及保护弱反射信号，我们采用最小二乘反演方法估计反射张角和方位角：

(19)

(20)

其中Ω={x, t}是关于时间和空间的局部平滑窗，在本文的测试中，仅采用了关于空间的平滑窗.

当方向矢量值较小时会导致反射张角计算错误，并且计算错误的反射张角几乎为90°.对于偏移成像，较大的入射角通常对应需要滤除的低波数噪音.因此，本文采用改进的最小二乘法来估计反射张角，把因方向矢量过小计算错误的角度值转移到更大的入射角度中，避免污染成像道集.最终计算反射张角的法方程为

(21)

1.4 角度域成像条件

经过正演和逆时外推后，可分别得到炮端解析波场和检波点解析波场，将分解后的行波分别估计波场方向，然后采用角度域互相关成像条件提取角度道集，最后叠加在一起获得最终的角度道集结果.以上下行波分解为例，采用波场分解得到的上下行波分别为炮端上行波场usup(x, t)、炮端下行波场usdown(x, t)和检波点端上行波场urup(x, t)以及检波点端下行波场urdown(x, t).对应的采用CLG光学流方法计算得到的归一化后的波场方向矢量分别为(x, t)、(x, t)和(x, t)、(x, t).我们采用如下的角度域互相关成像条件提取对应的角度道集：

(22)

其中θ为提取的地下反射张角，φ为方位角，θ1和φ1分别为(x, t)和(x, t)由公式(21)和(20)计算得到.θ2和φ2分别由(x, t)和(x, t)计算得到.θ3和φ3分别由(x, t)和(x, t)计算得到.θ4和φ4分别由(x, t)和(x, t)计算得到.

1.5 角度道集面元化和规整化

由于采集的不规则地震数据，复杂的地下传播路径，在地下角度域中成像道集是不规则的.在成像点处，不同角度面元和地下不同位置的角度道集覆盖情况是不同的.随着地下空间位置的不同，角度道集会出现欠采样和过采样问题，实际的地下角度采样是散乱不规则的.另外，由于数值计算的离散化问题，实际的角度道集不能放在正确的放置，不同的角度面元间隔会使得最终的不同角度的叠加次数不同，影响AVA分析.离散的、不规则的角度道集需要在角度域实施面元化处理.在本文中，我们采用角度插值方法对角度道集进行面元化和平滑.进一步地，还需要考虑角度照明和覆盖问题.

另外，由于偏移速度的不准确、地震子波的空间变化以及地下介质的复杂性，角道集可能会出现剩余深度时差和子波拉伸等现象.成像道集后处理技术，如动态时间规划和局部相似叠加(Wu et al., 2019)，可用于优化角度道集，提高成像道集质量.

2 数值实验

基于CLG光学流和波场分解提取角道集流程如图 1所示.考虑实际的多炮地震数据特征和当前的计算机集群的结构特点，本文采用MPI主从模式对炮数据并行，在每个单炮偏移过程中，采用OpenMp多线程并行策略.在单炮偏移过程中，首先将震源子波和该炮观测的地震数据进行希尔伯特变换，然后分别进行正演传播，得到炮点端和检波点端的解析波场.对解析波场进行波场分解，得到分解后的上下行波(在本文测试中，波场分解仅实施上下行波分解).然后采用CLG光学流方法分别计算上下行波场中的传播方向，并将得到的光学流矢量场实施归一化处理.然后采用公式(21)和(20)计算反射张角和方位角.最后采用角度域成像条件提取对应的角度道集.然后将生成的角度道集进行角度插值，避免角度面元化和数值化引起的误差.每个单炮成像道集结果存放在每个处理器附带的局部盘上.把所有炮成像道集数据体收集规约在一起，放到全局盘上，就可以产生多炮角度道集.最后，对角度道集进行必要的后处理，形成最终的结果.

图 1 本文提出的基于CLG光学流和波场分解的逆时偏移角度道集方法流程图 Fig. 1 The flow chart for extracting angle gathers in RTM based on the proposed CLG optical flow method with the wavefield decomposition
2.1 二维合成数据测试

首先采用层状模型说明本文提出的波场方向估计方法的有效性.其中水平层状速度模型如图 2a所示，采用有限差分正演模拟方法生成的震源端波场和检波点端波场快照如图 2b和2c所示.在该模型中，反射层深度为1500 m，分析可知只有反射界面附近的波前信息才对反射角道集有效.因此，该区域的波场方向必须被正确的估计.采用不同方法计算的反射界面附近的局部放大的检波点端波场估计的方向如图 3所示.其中图 3a是坡印廷矢量方法的结果，采用HS和CLG光学流方法估计的波场方向分别如图 3b和3c所示.在光学流方法中，考虑到LK方法可能出现奇异现象，因此在本文中，仅对比HS光学流方法和与本文提出的CLG光学流方法.从图 3a中可以看出，在大多数区域，采用坡印廷矢量估计的波场方向都是不正确的，波场方向不能被有效的识别，特别是图 3a中的红圈所示，该部分的波矢量方向杂乱无章，无法获得正确的波场角度值.坡印廷矢量方法计算的波方向矢量中，有效的波矢量较少，存在较多幅值较小的波方向矢量(图 3中，箭头之后线段的长短代表了波矢量的幅值大小)，该部分波矢量无法用于波场方向计算，如图 3a中绿圈所示.而HS光学流方法的计算结果则优于坡印廷矢量方法，其中HS光学流方法计算的波矢量幅值比坡印廷矢量方法更加平衡，波矢量方向更加准确.这是因为坡印廷矢量的结果只是对光学流场的初始估计(Zhang，2014).但是在图 3b绿圈所示中，仍然存在较多幅值不均的波方向矢量，导致波场方向估计的精度降低.另外，在图 3b中红圈中，杂乱无章的波方向矢量仍然存在，在该区域，计算的波场方向是错误的.对比HS光学流方法和CLG光学流方法，CLG光学流方法能更准确地估计波场方向，特别在图 3c中红圈位置，差异更为明显.CLG光学流方法计算的波矢量幅值比HS光学流方法更加平衡(如图 3c绿圈所示)，CLG光学流方法可以同时改善估计的方向和幅度.地震波场的方向在一定的范围具有相似性(子波的带限作用)，CLG光学流方法估计的波场方向更具有连续性，更符合地下介质波场传播的特征.其中0~60°角道集叠加的成像结果如图 4所示.由于坡印廷矢量方法波场方向不准确，导致在对应的成像剖面中出现了不该出现的大角度的低频噪音.而在HS光流法方法中，低频噪音被有效的去除.但是在大偏移距位置(如图 4中白圈所示)，成像结果不聚焦，弥散在成像剖面中.这是由于在该区域，波场的方向计算不准确导致(如图 3中所示).而在本文提出的方法中，成像效果较好，优于HS方法，从而说明了方向估计的有效性.

图 2 水平层状速度模型及相应的波场快照 (a) 层状速度模型；(b) 1.0 s时刻的震源端的波场快照；(c) 1.0 s时刻的检波点端的波场快照. Fig. 2 The layered velocity model and corresponding wave field snapshot (a) The layered velocity model; (b) The snapshot of the source wavefield at 1.0 s; (c) The snapshot of the receiver wavefield at 1.0 s.
图 3 采用不用方法计算的检波点端的波场方向 (a) 坡印廷矢量方法；(b) HS光学流方法；(c) 本文提出的CLG光学流方法. Fig. 3 The direction of the receiver wavefield estimated by different methods (a) Poynting vector method; (b) The HS optical flow method; (c) The proposed CLG optical flow method.
图 4 0~60°角度道集叠加结果 (a) 坡印廷矢量方法；(b) HS光学流方法；(c) 本文提出的CLG光学流方法. Fig. 4 The imaging results stacking by the 0~60°angle gathers (a) Poynting vector method; (b) The HS optical flow method; (c) The proposed CLG optical flow method.

接下来测试波场方向分解的有效性.图 5为采用上下行波波场分解方法得到的检波点端波场快照.图 6a为在上行波场(图 5a)中应用CLG光学流方法估计的局部波场方向(如图 5a中白框所示).图 6b为在下行波场(图 5b)中采用CLG光学流方法估计的局部波场方向(如图 5b中白框所示).图 6c为在未采用波场分解的检波点端波场中应用CLG光学流方法估计的波场方向(如图 2c中白框所示).可以看到，在图 6c中，波形交叉部分(图 6c中红圈所示)，仍然存在一些计算错误的波矢量方向.而在波场分解之后，该区域的波场方向得到正确的估计(图 6b中红圈所示)，有效提高了波矢量计算的精度，从而说明了波场分解的有效性.另外，图 6a中所示的波场方向与炮点端波场方向是相反的，互相关成像结果为大角度的低频噪音，在偏移成像中是需要去除的.因此，采用上下行波分解，可以有效地减少这些噪音的影响.

图 5 采用上下行波分解方法之后的检波点端的波场快照(1.0 s时刻) (a) 上行波波场快照；(b) 下行波波场快照. Fig. 5 The receiver wavefield with the wavefield decomposition method at 1.0 s (a) The upgoing wavefield; (b) The downgoing wavefield.
图 6 采用CLG光学流方法计算的检波点端的波场方向 (a) 采用上行波波场计算的波场方向；(b) 采用下行波波场计算的波场方向；(c) 未采用波场分解方法计算的波场方向. Fig. 6 The direction of the receiver wavefield estimated by the proposed CLG optical flow method (a) With the upgoing wavefield; (b) With the downgoing wavefield; (c) Without the wavefield decomposition.

接下来测试本文提出的角度估计的有效性.图 7a为常规除结果.本文提出的最小二乘反演除的结果如图 7b所示(图 7a和图 7b均采用了常规的HS光学流方法计算的波矢量结果来估计反射张角).由于带限子波的震荡作用，常规方法计算的角度间断、不连续.采用反演的方法估计的角度值更加稳健精确.子波的带限作用引起的角度误差得到有效的消除.

图 7 层状模型中估计的反射张角结果 (a) 常规的方法；(b) 本文提出的最小二乘除方法. Fig. 7 The estimated reflection angle in the layered model by (a) conventional method and (b) the proposed least square division method

然后采用Sigsbee 2A理论模型说明本文提出方法的有效性.图 8a为采用常规的HS光学流方法生成的角度道集，采用本文提出的CLG光学流方法+波场上下行波分解生成的角度道集如图 8b所示.可以看出，常规的HS光学流方法的角度道集中含有较多噪声，模糊现象严重，角度道集弥散在比较广的角度中，特别是在盐丘的左侧和上部，成像道集质量比较差.然而在采用本文提出方法的成像道集中，角度道集比较连续、聚焦.对应地0~60°角度道集叠加结果如图 9所示.可以明显地看出，采用CLG光学流方法+波场上下行波分解之后的叠加成像中，假象更少，反射轴更加清楚.

图 8 Sigsbee 2A模型每隔100个CDP点提取的角度道集(角度范围0~60°，角度间隔1°) (a) 常规的HS光学流方法；(b) 本文提出的CLG光学流+波场分解方法. Fig. 8 The ADCIGs of Sigsbee 2A model. The ADCIGs is selected with a 100-CDP interval and the angle range is 0~60 degrees with a 1 degree interval (a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.
图 9 Sigsbee 2A模型0~60°角度道集叠加结果 (a) 常规的HS光学流方法；(b) 本文提出的CLG光学流+波场分解方法. Fig. 9 The stacked image of Sigsbee 2A model by ADCIGs. (The angle range is 0~60 degrees) (a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.

在Sigsbee 2A模型上分析不同角度道集生成方法的单炮计算时间，如图 10所示.可以看出，坡印廷矢量方法计算效率最高，常规的HS光学流方法次之.而本文提出的CLG光学流方法+波场分解方法需要的计算时间则是最长.这是由于在本文提出的方法中，需要采用解析波场和空间傅里叶变换对外推的波场进行行波方向分解.波场方向分解依赖于解析波场，解析波场需要额外一倍的外推时间.本文的空间傅里叶变换是对整个波场进行分解，不同于局部的波场分解方法(有关局部波场分解方法和常规的波矢量方法的效率对比可参考文献(吴成梁等，2018))，需要的计算量并不大.另外，CLG光学流的计算，最小二乘角度估计，也都需要一定的计算量.对比常规的HS光学流方法和本文提出的方法，本文提出的CLG光学流+波场分解方法的计算效率是常规HS光学流方法的3倍左右.对比成像效果的提升和当前计算机的快速发展，本文认为这些计算代价是值得的.

图 10 采用不同方法计算的单炮时间对比 Fig. 10 Comparison of single shot time calculated by different methods
2.2 二维实际资料测试

接下来，采用某地区陆上实际资料说明本文提出方法的有效性.其中图 11a为常规的HS光学流方法生成的角度道集，图 11b为本文提出的方法生成的角度道集.对比可以看出，采用本文提出方法生成的角度道集假象比较少，成像道集比较聚焦.另外，从叠加成像结果(图 12)可以看出，常规的光学流方法结果中，存在较多的偏移假象，叠加剖面比较模糊，反射波同相轴不连续.在本文提出的方法中，成像结果连续性较好，剖面噪声和偏移假象较少，成像质量结果明显好于常规的光学流方法.

图 11 陆上某地区2D实际资料，每隔100个CDP点提取的角度道集(角度范围0~60°，角度间隔1°) (a) 常规的HS光学流方法；(b) 本文提出的CLG光学流+波场分解方法. Fig. 11 The ADCIGs of real data. The ADCIGs is selected with a 100-CDP interval and the angle range is 0~60 degrees with a 1 degree interval (a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.
图 12 陆上某地区2D实际资料，0~60°角度道集叠加结果 (a) 常规的HS光学流方法；(b) 本文提出的CLG光学流+波场分解方法. Fig. 12 The stacked image of real data by ADCIGs. (The angle range is 0~60 degrees) (a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.
3 讨论与结论

利用波动方程的时间演化进行波矢量方向计算在逆时偏移角度道集中是一个核心内容.局部波场分解方法能够处理复杂波场，传播方向计算准确，但是效率低下，无法在大规模的3D逆时偏移角度道集中实用化.坡印廷矢量方法可以用于计算波场传播方向，该方法计算简单，计算效率高，适合并行处理，但是坡印廷矢量通常会出现不稳定情况，在波场复杂情况时，波场方向估计不准确.

光学流方法可以提高波场方向估计的精度.光学流方程本身是欠定的，需要引入额外的约束计算光学流场.通过假设在局部邻域Ω内光学流矢量是不变的，LK光学流方法构建局部的线性方程组计算方向矢量.然而，LK光学流方法依赖于局部点的梯度值，容易出现奇异现象.全局HS光学流方法通过引入全局的能量约束，能够比较稳定地求解光学流场，但是在该方法中，距离计算点较远区域的波场也会参与到当前点的计算中，导致计算的波场方向精度有所下降.本文采用CLG光学流方法计算波场方向，在施加全局约束的基础上，考虑局部点的加权作用，在保持波场方向估计稳定的基础上，提高波场方向计算精度.该方法简单高效，便于并行处理，其迭代格式与HS光学流方法相同，几乎没有太多的额外计算量，非常适合大规模地震数据处理.

为了克服光学流方法无法处理波前重叠问题，本文利用解析波场对波场进行方向分解来减弱这种影响.通过构建解析波场，仅需波场的空间傅里叶变换即可实现任意波场方向分解，有效地避免了逆时偏移外推过程存储波场问题，减少了内存和硬盘使用.理论上可以通过波数域滤波提取任意方向的波场，但是考虑到计算量问题以及地下介质的特征主要以层状沉积层加上大、小尺度的速度异常体，本文的数值实验仅考虑了上下行波分解情况.针对不同的地质构造特征，可以任选几个方向角度进行分解，波场分解的角度需要根据地质特征确定.方向分解个数越多，波前交叉分的更仔细，但是需要的计算量越大.波场分解的个数需要在计算量和效果中达到平衡.

进一步地，本文提出采用有效的归一化方法和最小二乘除方法，提高反射角的估计精度，有效避免带限波场引起的震荡作用.并对角度道集进行了面元化和规整化处理，以获得高质量的角度道集.理论和实际资料证明了本文提出方法的有效性.

附录A：CLG光学流控制方程推导

公式(8)中所示的CLG光学流误差泛函中的可以写为

(A1)

进一步的，可表示为：

(A2)

重新整理，误差泛函ECLG(ω)可表示为

(A3)

采用变分法，上式误差泛函对应的Euler-Lagrange方程的解为

(A4)

最后，得到CLG光学流控制方程为

(A5)

其中Δ为拉普拉斯算子.

致谢  感谢中石油勘探开发研究院及西北分院、中海油研究院和湛江分公司、中石化物探技术研究院和胜利油田分公司对波现象与智能反演成像研究组(WPI)研究工作的资助与支持.感谢审稿人提供的宝贵建议.
References
Barron J L, Fleet D J, Beauchemin S S. 1994. Performance of optical flow techniques. International Journal of Computer Vision, 12(1): 43-77. DOI:10.1007/BF01420984
Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration. Geophysics, 48(11): 1514-1524. DOI:10.1190/1.1441434
Brox T, Bruhn A, Papenberg N, et al. 2004. High accuracy optical flow estimation based on a theory for warping. //Proceedings of the 8th European Conference on Computer Vision. Berlin, Heidelberg: Springer, 25-36.
Bruhn A, Weickert J, Schnörr C. 2005. Lucas/Kanade meets Horn/Schunck: Combining local and global optic flow methods. International Journal of Computer Vision, 61(3): 211-231.
Dickens T A, Winbow G A. 2011. RTM angle gathers using Poynting vectors. //81st Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3109-3113.
Duveneck E. 2013. A pragmatic approach for computing full-volume RTM reflection angle/azimuth gathers. //83rd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
Fomel S. 2004. Theory of 3-D angle gathers in wave-equation imaging. //74th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1053-1056.
Horn B K P, Schunck B G. 1981. Determining optical flow. Artificial intelligence, 1981, 17.1-3: 185-203.
Hu J T, Wang H Z, Wang X W. 2015. Angle gathers from reverse time migration using analytic wavefield propagation and decomposition in the time domain. Geophysics, 81(1): S1-S9.
Jin H, McMechan G A, Guan H M. 2014. Comparison of methods for extracting ADCIGs from RTM. Geophysics, 79(3): S89-S103. DOI:10.1190/geo2013-0336.1
Liu F Q, Zhang G Q, Morton S A, et al. 2011. An effective imaging condition for reverse-time migration using wavefield decomposition. Geophysics, 76(1): S29-S39. DOI:10.1190/1.3533914
Liu F Q, Hanson D W, Whitmore N D, et al. 2006. Toward a unified analysis for source plane-wave migration. Geophysics, 71(4): S129-S139. DOI:10.1190/1.2213933
Lucas B D, Kanade T. 1981. An iterative image registration technique with an application to stereo vision. //Proceedings of the 7th International Joint Conference on Artificial Intelligence. San Francisco, CA, United States: Morgan Kaufmann Publishers Inc., 674-679.
Patrikeeva N, Sava P. 2013. Comparison of angle decomposition methods for wave-equation migration. //83rd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3773-3778.
Sava P, Fomel S. 2003. Angle-domain common-image gathers by wavefield continuation methods. Geophysics, 68(3): 1065-1074. DOI:10.1190/1.1581078
Sava P, Fomel S. 2006. Time-shift imaging condition in seismic migration. Geophysics, 71(6): S209-S217. DOI:10.1190/1.2338824
Sobey P, Srinivasan M V. 1991. Measurement of optical flow by a generalized gradient scheme. Journal of the Optical Society of America A, 8(9): 1488-1498. DOI:10.1364/JOSAA.8.001488
Szeliski R. 2010. Computer Vision: Algorithms and Applications. New York: Springer-Verlag New York Inc. .
Tang B, Xu S, Zhang Y. 2013. 3D angle gathers with plane-wave reverse-time migration. Geophysics, 78(2): S117-S123. DOI:10.1190/geo2012-0313.1
Vyas M, Du X, Mobley E, et al. 2011a. Methods for computing angle gathers using RTM. //81st Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, F020.
Vyas M, Nichols D, Mobley E. 2011b. Efficient RTM angle gathers using source directions. //81st Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3104-3108.
Wang B L, Gao J H, Chen W C, et al. 2013. Extracting efficiently angle gathers using Poynting vector during reverse time migration. Chinese Journal of Geophysics (in Chinese), 56(1): 262-268. DOI:10.6038/cjg20130127
Wang H Z, Feng B, Wang X W, et al. 2015. Analysis of seismic inversion imaging and its technical core issues. Geophysical Prospecting for Petroleum (in Chinese), 54(2): 115-125, 141.
Whitmore N D. 1983. Iterative depth migration by backward time propagation. //53th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 382-385.
Wu C L, Zhou Y, Hu J T, et al. 2018. Efficient generation of reverse time migration angle gathers. Geophysical Prospecting for Petroleum (in Chinese), 57(3): 404-418.
Wu C L, Wang T Z, Wang H Z, et al. 2019. Improving the quality of common-image gathers using DTW and local similarity. //89th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 4336-4340.
Xie X B, Wu R S. 2002. Extracting angle domain information from migrated wavefields. //72nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1360-1363.
Xu S, Zhang Y, Tang B. 2011. 3D angle gathers from reverse time migration. Geophysics, 76(2): S77-S92. DOI:10.1190/1.3536527
Yan R, Xie X B. 2011. Angle gather extraction for acoustic and isotropic elastic RTM. //81st Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3141-3146.
Yoon K, Marfurt K J, Starr E W. 2004. Challenges in reverse-time migration. //74th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1057-1060.
Yoon K, Marfurt K J. 2006. Reverse-time migration using the Poynting vector. Exploration Geophysics, 37(1): 102-107. DOI:10.1071/EG06102
Yoon K, Guo M H, Cai J, et al. 2011. 3D RTM angle gathers from source wave propagation direction and dip of reflector. //81st Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3136-3140.
Yoon K. 2017. Reverse time migration angle gathers using Poynting vector and pseudospectral method. //87th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 4630-4634.
Zhang Q. 2014. RTM angle gathers and specular filter (SF) RTM using optical flow. //84th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3816-3820.
Zhang Q S, McMechan G A. 2011a. Common-image gathers in the incident phase-angle domain from reverse time migration in 2D elastic VTI media. Geophysics, 76(6): S197-S206. DOI:10.1190/geo2011-0015.1
Zhang Q S, McMechan G A. 2011b. Direct vector-field method to obtain angle-domain common-image gathers from isotropic acoustic and elastic reverse time migration. Geophysics, 76(5): WB135-WB149. DOI:10.1190/geo2010-0314.1
Zhao L, Feng B, Wang H Z. 2012. Efficient RTM angle gathers using source propagation direction and reflector normal direction. //82nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1-5.
王保利, 高静怀, 陈文超, 等. 2013. 逆时偏移中用Poynting矢量高效地提取角道集. 地球物理学报, 56(1): 262-268. DOI:10.6038/cjg20130127
王华忠, 冯波, 王雄文, 等. 2015. 地震波反演成像方法与技术核心问题分析. 石油物探, 54(2): 115-125, 141. DOI:10.3969/j.issn.1000-1441.2015.02.001
吴成梁, 周阳, 胡江涛, 等. 2018. 高效逆时偏移角度道集生成方法研究. 石油物探, 57(3): 404-418. DOI:10.3969/j.issn.1000-1441.2018.03.010