广义Rytov近似有限频初至层析

2021年2月哀油狄旅在探第56卷第1期.偏移成像•文章编号：1000-7210(2021)01-0092-08广义R y t o v近似有限频初至层析
冯波①② W U R u-S h a n③罗飞*①许荣伟①王华忠①
(①同济大学海洋与地球科学学院，上海200092;②同济大学海洋高等研究院，上海200092;
<3)Modeling and Imaging Laboratory,University of California%Santa C r u z.C A 95060,U S A)
冯波，W U R u-S h a n,罗飞，许荣伟，王华忠.广义R y t o v近似有限频初至层析.石油地球物理勘探，2021,56(1):
92-99.
摘要在传统的有限频层析理论中，采用B o r n或R y t o v近似推导地震波旅行时对模型参数的一阶F K c h e t微 商，隐含了弱散射近似假设。

对于强速度扰动，有限频层析正问题的线性化近似不再成立。

近年来发展的广义
R y t o v近似理论可以较为精确地预测散射波的相位扰动，因此更适用于有限频旅行时层析。

为突破弱散射近似
的限制，将广义R y t o v近似应用于传统的有限频层析理论，导出了基于广义R y t o v近似理论的有限频旅行时敏
感度核函数。

由数值试验可知，对于前向小角度散射，无论速度扰动强度如何，新的旅行时核函数均可较为准
确地预测地震波的旅行时扰动。

针对旅行时层析反问题，提出了一种基于隐式矩阵向量乘的高斯一牛顿反演
算法，具有准二阶收敛速度且无需显式计算和存储H e s s i a n矩阵。

复杂模型数据测试结果表明，提出的广义
R y t o v近似旅行时层析可以建立高精度的近地表速度模型，且收敛速度较快。

关键词有限频旅行时敏感度核函数广义R y t o v近似初至层析高斯一牛顿反演算法
中图分类号：P631文献标识码：A doi： 10. 13810/j. cnki. issn. 1000-7210. 2021.01.011
〇引言
地震波旅行时常用于反演地下速度结构。

在射 线层析中，由于引人了高频近似假设，旅行时只与连 接炮检对的几何射线上的速度扰动有关。

对于有限 频带的地震波，若速度非均匀体的尺度小于菲涅耳体的宽度时，地震波散射起主导作用。

由此发展了 绕射层析理论[12:;，可以较好地处理地震波的一阶绕 射效应。

基于B o m近似，L u o等[3]给出了地震波旅行时Fr6c h e t导数（旅行时敏感度核函数），可以 定量描述速度扰动对地震波旅行时扰动的影响。

通 过波动方程线性化近似（B o rn近似或R y tov近似），W〇〇dward[l1提出了“波路径”的概念，作为无限高频假设下的射线路径向有限频地震波传播的推广。

借助伴随状态理论，可以实现（波形、旅行时、振幅 等）F K ch et导数的高效计算「5-11]。

通常旅行时F r6c h e t导数由B o r n近似导出[1218]，但B o rn近似成立条件较为苛刻。

相比于 B o rn近似，R y to v近似可以更好地描述由于速度扰 动引起的前向散射波的相位扰动n9]，因而更适用于 旅行时反演。

基于R y to v近似也可以构造相应的旅行时敏感度核函数[2°27]。

然而，B o m近似与R y to v近似均属于弱散射近似，要求速度扰动远小于背景速度128293。

为克服弱散 射近似的制约，Feng等[3°:1提出了可以适用于强扰动 模型的广义R y to v近似理论（Generalized Rytov ap­pro x imation，GRA)。

对 于前向 散射及 小角度 传播，G R A可以较为准确地预测地震波的相位扰动，因而 更适用于旅行时反演。

本文将GRA理论应用于传统
*上海市四平路1239号同济大学海洋与地球科学学院，200092。

Email:丨uofeil9901217@126. c om
本文于2020年3月20日收到，最终修改稿于同年11月1日收到。

本项研究受国家重点研发计划"深海关键技术与装备”重点专项“深海油气地震勘探高精度速度建模技术”（2019Y F C0312004)、国家自然科 学基金项目“强散射介质中的波动方程线性化及速度反演”（42074143)和"针对勘探地震介质及波场特征的高波数参数估计方法研究”(41774126)、国家重点研发计划“变革性技术关键科学问题”重点专项“多元信息联合驱动的地震成像方法研究"（2018Y F A0702503)、上海 市浦江人才计划（20PJ141300)、南方海洋科学与工程广东省实验室（湛江）项目“特征波域地震成像方法研究’’（ZJW-2019-04)、中国石化地 球物理重点实验室项目“基于波动方程的高精度近地表速度建模方法研究”（33550006-19-FW0399-0041)联合资助。

第56卷第1期冯波，等：广义Rytov近似有限频初至层析93
的有限频层析，构造一种适用范围更广的旅行时敏感 度核函数并用于初至旅行时层析。

模型数据测试结 果表明，本文提出的广义Rytov近似旅行时层析可以 建立高精度的近地表速度模型，且收敛速度较快。

1理论
1.1广义Rytov近似有限频旅行时敏感度核函数
频率域标量声波方程可以表示为
(▽2+々2)“。

=0(1)式中：H jc)=〇>/x;g(jc)，为背景速度场中的地
震波数；《〇(•*■，〇；)为背景波场；V2为Laplace算子。

记扰动后的速度场为w(x)，波场（总场）为 M(X，O J)，总场与背景场的关系为
，a»)=M。

（J C，£U)eXp[0(J C，<W)] (2)式中0(J C，W)为复相位。

在前向散射(散射角定义如图1所示)及小角度传 播假设下，复相位可以用广义Rytov近似[31]表示为
</P R A(x,w)
=J a^C o c )X
2⑴、。

⑷(3)
u〇 (x,co)
式中：，〇；;幻为w C x)中从X到〇:的格林函数；
A.s( jt) =s(x) —5。

（x)为慢度扰动，其中 5。

（j c)=队""1(X)、.S(X)二"^1 (A T)分别为背景慢度和扰动后的慢度。

图1散射角定义：地震波人射方向与散射方向的夹角前向散射对应小散射角，背向散射对应大散射角
式（3)引起的相位延迟可以表示为
A,p(⑴)=—M，)] (4)
C0
有限频带地震信号的旅行时扰动可以表示为单频谐 波的相位延迟的加权叠加[23#’32]，有
A?= J W(〇))A；p(<w)dc« (5)式中W(W)=o>2S(⑴)/J'a/S(a〇da>为归一化加权函数，其中S(w) = /0)厂(〇>)为震源子波能谱，其中 上标”表示复共扼。

结合式（3)〜式（5)，可以建立带限信号旅行时扰动与慢度扰动的线性关系
A f(»-^*)=j A.v(J：)K l,RA (jc；jrr,)djc(6)式中:和a分别为源和检波点坐标;(jc;a，x s)为 广义R ytov近似旅行时敏感度核函数（Generalized Rytov approximation based traveltime sensitivity kernel，简称为G R A核函数），可表示为
KnRA U;u«)=|w(⑴）丄 X
J C i)
Im[-2co2s A x)(7) _ U〇(x r ,C O；X,)」
显然，显式计算式（7)需要计算频率域格林函数[25]，因而计算量较大，尤其是对于三维问题。

为了避免显式计算核函数，本文采用冯波等[26]给出的方法，通过引入如下辅助函数
q(x r ,cv；x s):
W(o j)
(8)
jc〇u〇(x T,C O；x s).
利用Parseval定理，将式（7)中对角频率的积分转化为对时间积分，则
K(；R A(jc；jcr,j c s)
=Re|J(/*(x T,w；j c s)[ —2c u2s〇 (x)u〇(x9cu；xs)X G〇(x,c〇；x T)]dwj
27rR e|J25〇(jc)^fw〇(jc,/；jc s)X
[go(x，—，|冬，〇) * 9(冬山丨s)] * 士}
2sQ{x)d2t u〇^r〇 (a t,— ^|x r,0)* R e[2丌(x r，，;a t s ) ] | d，(9)式中：T表示地震记录时长；A表示时间域格林函数；“”表示褶积（与上标“”表示复共轭不同）。

若定义旅行时伴随场为
u r(x,T— t；x r)
=g〇(x,—t\x T,0)^R e[27T^*(ATr,/；Jfs)]
=(JC，了一，I x r，T)* R e[2丌g* (x r “;o:s)]
(10)式（9)可以重写为
K GRA(j c；x r,Jcs) =
'T
2s〇(x)^u〇(x ,t；x s)u T(x ,T —t；x r)d t
(11)
94石油地球物理勘探2021 年式中旅行时伴随场由旅行时伴随震源/(A，^A：S)=
Re[27i^ (A，T—Gxs)]逆时传播产生。

1.2旅行时层析反问题求解
基于旅行时残差U范数的误差泛函可表示为
1V d1
m inj(m) = — V]A/f =1(12)
m6M 乙i~,L
式中:m=(w丨，/w2，)T6M,为模型空间M中的 向量，本文为慢度模型(采用网格参数化方式），其中iV v 为网格点数;A f=(M,M，…，A A v d)Te i>.为数据空间 D中的向量(旅行时残差），其中Nd为炮检对数。

利用Gauss-N ew ton反演算法，式（12)可以迭 代求解
m M=m*)V7 (/«*)] (13)式中4为迭代次数；a为步长，可以通过线性搜索方法计算；V/为目标函数梯度；/f,=^= K1*：为线性H essian矩阵（忽略旅行时对模型参数的二阶导数），其中*：为K°R A的矩阵形式；P为光滑 算子。

根据式（12)，泛函梯度可表示为
V；(W l t)= l^= (^)，Afi=K，A^ (14)基于广义R y to v近似，即为式（11)给出的 核函数。

将式（11)代人式（14)，则泛函数梯度可表示为
V J C w'X jc)=(K T A^ )(x)
=2 2k gra O;h)厶"（H)
x,x,
rT
=X] 2sk(x)〇l u〇(x,t；x, )X
a(x，T—t；x^) A t(15)式中 a(x，T—f ;jcs ) = 2反r(jc，T—〇jcr)A/* (Jcr，jcs)为
X
单炮旅行时伴随波场。

考虑到显式计算H essian矩阵需要巨大的计算 量及存储量，且求逆困难，因此，将对H essian矩阵 的求逆转化为求解如下法方程的近似解
H A m k)Amt =VJCw*) (16)
本文用共扼梯度（Conjugate gradient，CG)法求 解式（16)获得模型更新方向，并采用隐式矩阵向量乘方法[26]直接计算H essian矩阵与向量的乘，因而 无需显式计算及存储H essian矩阵。

具体计算公式 详见附录A。

2数值试验
2.1高斯扰动模型
图2为一个含有高斯异常体的二维10km X 5k m速度模型，网格间距为10m。

其背景速度％= 3000m/s，异常体的速度扰动可表示为
△u(，,z) = £ii0exp(—^-7j(17)
式中：e为速度扰动百分比；a=l k m为髙斯异常体的尺度参数；r=s/(x—X〇)2+(z—Z〇)2 ,{x〇,Z〇)= (5. 0km,2.5 km)为高斯异常体的中心坐标。

o
4
图2含有高斯扰动（e=50%)的速度模型及观测系统
红色倒三角为检波点位置
首先，用时一空域高阶有限差分求解声波方程，用真实模型（包含高斯异常体的速度模型）和背景模 型（即均匀背景速度模型）进行正演（震源是主频为15H z的R icker子波），得到的地震记录分别作为观 测数据和模拟数据。

然后，用互相关（Cross-corre­lation，简记为 CC)函数测量观测数据和模拟数据的 时差，作为真实时差的测量值（即观测时差）。

随后，分别用射线（Ray)核函数（速度网格内的射线长度）、R ytov近似旅行时敏感度核函数及本文提出的G R A旅行时敏感度核函数预测旅行时扰动(即核函数与模型扰动做内积）。

图3为震源在地表 (5.01«11，0)和（0,0)处、201个检波器均匀分布在z=5. 0km处（图2中红色倒三角所示）旅行时扰动 对比。

由图3可以看出：①射线和G R A核函数的 预测时差完全重合，但与R y to v近似预测结果有较大偏差。

②对于小角度前向散射（图3a中对应为：震源位于地表（5.0km，0 )，检波器横坐标在(5. 0km，5.0km)附近），G R A核函数的预测时差与 观测时差基本一致。

证明了 G R A核函数在前向小 角度散射时具有较好的线性特征，克服了传统波
动
第56卷第1期冯波，等：广义Rytov近似有限频初至层析95
方程线性化近似中的速度小扰动假设的制约。

这与 G R A理论预测相符（G R A理论指出，对于小角度前 向散射，G R A可以较为精确地预测地震波的相位扰 动）。

③随着散射角增大，G R A核函数的预测时差逐渐偏离观测时差，并逐渐与R yto v近似预测结果趋于一致。

根据G R A理论推导可知[3«，Ryt〇v近 似为G R A的一阶近似，本质上都隐含了小角传播假设。

因此，随着散射角增大，G R A的成立条件（即前向小角散射）不再成立，导致预测时差偏离观测值。

对于图3b中的震源位置（左上角（0,0)处），前 向小角度散射对应检波器横坐标在l〇k m附近，三 类核函数预测结果与图3a—致。

(a)
图3震源在地表不同位置三种旅行时敏感
度核函数的预测时差与测量时差对比
(a)(j：s,zs)=(5. 0k m,0) ；(b)(j：s=(0,0)
上述测试表明：G R A核函数精度与散射角密切 相关。

对于前向小角度散射，G R A可以较为准确地 预测地震波的旅行时扰动。

随着散射角不断增大，G R A核函数精度逐渐降低，但仍好于R y to v近似。

虽然对于大角度散射，G R A核函数的预测结果并不 准确，甚至与观测时差偏离较远，但并不会导致层析 反问题求解失败。

原因在于：层析反问题通过最优化算法迭代求解时差目标函数，而非直接对线性方程组直接求逆。

即使正问题的线性化近似不严格成 立，只要反问题的凸性较好，依然可以采用最优化算 法求解目标泛函。

为了验证G R A初至层析方法的有效性、尽可 能消除观测系统对反演结果的影响，设计了 一个理想观测系统如图4a所示，100个检波器均匀分布高斯异常体模型的四周，围绕模型四周激发1〇〇次，震 源是主频为15H z的R ick er子波（对应波长A。

= 0.2km，与高斯异常体尺度参数关系为a = 5A(,)。

jc/km
2 4
(b)
图4高斯模型观测系统(a)及G R A旅行时反演结果（b) 黑色正三角形表示震源位置；红色倒三角形表示检波点位置
反演采用的初始模型为均匀背景速度场（％= 3000m/S)，并用互相关方法测量观测数据和模拟数据的时移量作为初至旅行时残差。

反演采用SEIS~ COPE 数值优化软件包（SEISCOPE optimization toolbox)中的截断牛顿算法（对于旅行时层析，截断 牛顿算法退化为高斯牛顿算法）实现最优化，其中，梯度和H essian矩阵与向量乘由附录A中给出的隐式矩阵与向量乘方法计算。

在迭代反演中，采用 2次C G内迭代求解式（16)，获得模型更新方向，并 用线性搜索方法计算迭代步长。

反演收敛准则为判断相对旅行时误差是否小于预设的
阈值0X10 4)。

经过10次迭代（图5为反演 收敛曲线），反演结果如图4b所示，与真实速度模型
96石油地球物理勘探
2021 年
(a )
jc/km
0 1 2
3 4 5
1500
1000
500
jc/km
(a )
x/km
A v /(m -s ')
I
400
300 200
(b)
图7 Marmousi -2模型（a )及G R A 初至层析的速度更新量（b )
迭代次数
图5 G R A 初至层析的收敛曲线
2. 2 Marmousi-2 模型
应用Marmousi -2模型验证G R A 初至层析 对复杂速度模型及不完备的观测系统的有效性。

为 模拟陆上地震釆集，切除了原始速度模型中上覆的 水体部分，并在模型右侧扩边3k m (图7a )。

x 和
z 方向网格数目分别为2001和301，网格间距均
为 10m o
图6高斯模型层析（蓝色）与真实（红色）速度扰动曲线
(a)a- = 2. 5k m ； ( b )z = 2. 5k m
设计了一个单边观测系统：震源和检波器均放 置在地表，起始炮点位于x = 0,炮间距为100m ，共 171炮；最小和最大炮检距分别为400m 和5000m ， 道间距为20m (每炮共231道）。

最大炮检距较大， 因此可以用初至波旅行时反演近地表速度。

采用有 限差分求解波动方程正演地震信号，震源是主频为 15H z 的Ricker 子波。

由于观测地震记录较为复杂，
在视觉上几乎没有差别。

为定量对比反演结果，对 反演的扰动模型（反演结果减去初始模型）与真实的 高斯异常体进行横向和纵向抽线对比（图6)，反演 的速度扰动与真实的高斯扰动基本吻合，证明了对 于完备的观测系统及简单高斯速度模型，G R A 初至 层析可以得到高精度的反演结果，证实了方法的有 效性。

Av/(m s )500 1000
1500
(
-J />
v
第56卷第1期冯波.等：广义Kytov 近似有限频初至层析97
逼近真实的高速隆起，反演精度较高（图9a 和图9b )。

随着深度增加，反演精度逐渐降低（图9c 和图9d ),但 反演结果仍然与真实速度模型的整体趋势相吻合。

M a rm o u si-2模型实验结果表明，利用中、小炮
检距的初至波旅行时信息.G R A 初至层析可以建立 高精度的近地表速度模型。

1.0 —
0.2
3 4迭代次数
采用自动方法[3]:拾取初至，作为观测记录的初至旅 行时。

初始模型采用随深度线性递增的速度场U =〇, ■〇= 1000m /s ; z = 3. 5km ，D = 3000m /s )。

为降低计 算量，采用20m X 20m 网格剖分进行速度反演，并用 相同的初至拾取方法拾取模拟数据的旅行时，进而计 算初至时差。

反演实现框架与高斯模型一致（相对旅 行时误差阈值〃=1.0父10—2)。

经过5次迭代反演收 敛，得到的速度更新量如图7b 所示，浅层的高速异常 得到了较好的恢复。

其收敛曲线如图8所示，可以看 出，高斯一牛顿算法收敛速度较快。

为定量对比反演结果，抽取z =0. 1、0.2、0.4、 0.6km 处速度横向变化曲线，如图9所示。

可以看 出，
3
旅行时（差）的精确测量是所有旅行时反演类方 法所共有的问题。

本文的重点为讨论如何将广义
R yto v 近似方法与有限频层析理论的结合，是笔者
过去研究工作[26]的理论推广和应用，因而没有讨论 如何准确测量初至波时差。

若震源子波已知且模拟 数据与观测数据的初至波形基本一致，此时测量初 至波时差相对简单。

反之，若震源子波未知或由于 地下介质复杂导致初至波形畸变严重，此时需要采
4结论
本文提出了一种基于广义R y to v 近似的有限
频旅行时敏感度核函数.克服了传统有限频层析方 法中弱速度扰动假设的制约。

对于前向散射及小角 度传播.广义R ytov 近似旅行时敏感度核函数可以 较为准确地预测地震波的旅行时扰动，从而可以提
&表U =0. 1、0. 2km ),反演结果（蓝色曲线)
>300
1500
>300 *
x /k m
(a )
j 1900
$
2700
E
2300
1900
1500
图8 G R A 初至层析的收敛曲线
j ：/k m
(b )
4 8
12
jc/k m
(d)
图9不同深度处初始(绿色）、真实（红色）、G R A 层析（蓝色）速度横向变化曲线对比
(a)z = 0. l k m ；(b)z=0. 2k m ；(c)z = 0. 4k m ；(d)z=0. 6k m
16
用与震源无关的反演类算法消除子波波形畸变对旅
行时测量的影响，如双差旅行时（double-difference
traveltim e)[33 34 测量等。

6
4
0.0.
(S .J /A
(->£)/
>
98石油地球物理勘探.2021 年
高反演精度并加速反演收敛。

在计算方面，本文采 用了隐式矩阵向量乘方法直接计算H essian矩阵与 向量的乘，无需显式计算及存储H essian矩阵，进而 可以提高高斯_牛顿迭代反演算法的效率。

复杂的 Marmousi-2模型测试结果表明，G R A初至层析可以建立高精度的近地表速度模型。

附录A隐式矩阵向量乘
H essian矩阵与向量积可表示为
H p =K JKp =KTh(A-l)式中：p=p(x)为模型空间中的向量;=为数据空间中的向量。

因此需要计算两次矩阵向量积。

首先，V(J^，jcs)，K：p可以表示为
(Kp ) (x r ,xs) —/)(x)Re| —2<ifs〇(x)ti<)(x,cu；x f)X
\_G〇{x,w,x r)q'( jc, , c o；jc s )] dco | dor
=Re Jg* (x r, co； x s)m p (x r ,cu；x5 )dt〇
(A-2)式中《[)(_»^，£0;_^)满足如下13〇1'11近似
m p (x r,w；jc s ) =— 2c〇zs〇 (x)p | , X
G〇(J C，a>;jC r)m〇(X，⑴；j c s)dx(A-3)利用Parseval定理，式（A-2)可改写为
-厂了-(Kp)(xr,x s)= 27rRe q* (xr)W p(jc r)d/
-J〇_ *T
=uq(x T ,t；x s)u p(x T ,t；x,)d t(A-4) Jo
式中
M q(x r“；x s)=2丌Re[g* (x r，，；x s)]
=2丌Re [F T^1 g(x r，w;x s) ] *(A-5) m p(A：r)=FT-1m p(jcr,w；j c s)(A-6)同时，式（A-3)在时一空域满足如下方程组
| [劣772 (JC) — V 2 ]M〇(X山X s) = /W d(X—JCS)
\\_d2t m{x) —V2~\uv{x-it；xs) =2s〇(x)p(x)d^u〇(x,t；xs)
(A-7)式中/为震源函数。

其次，根据式（15)，对于空间任意一点jc，/C T/i 可以表示为
(K J p)(x)=^2K G R A(jc；jcr,x s)/i (jcr,jc s )
X s X r
=S [2s〇(x)d^u〇(x U；x s)u h(x,T — t；x,)d t Xs J〇
(A-8)式中
uh(x,t；x s)={g〇(x,—t\x r,0)^
X T
[/ (x r, ^ ；jcs) /z (jcr, jcs) ]} (A-9)
为单炮旅行时伴随波场。

参考文献
[1 ] Slaney M,K a k A C, Larsen L. Limitations of imaging
with first-order diffraction t o m o g r a p h y [J] .I E E E
Transactions O n M i c r o w a v e T h e o r y and Techniques,
1984,32(8)：860-873.
[2]W u R S, T o k s o z M N.Diffraction t o m o g r a p h y and
multisource holography applied to seismic imaging
[J] .Geophysics, 1987,52(1) ：11-25.
[3]L u o Y,Schuster G T.Wave-equation travel time in-
version[J] .Geophysics, 1991,56(5) ：645-653.
[4] W o o d w a r d M J. Wave-equation tomography [J] .Geo-
physics, 1992,57(1)：15-26.
[5]Lailly P. T h e seismic inverse problem as a sequence
of before stack migrations [C] .E x p a n d e d Abstracts
of Conference o n Inverse Scattering, T h e o r y and
Application, 1983,206-220.
[6]Liu Q,T r o m p J.Finite-frequency kernels based u p o n
adjoint m e t h o d s[J] .Bulletin of the Seismological S o­
ciety of America, 2006,96(6) :2383-2397.
[7]Liu Q, T r o m p J.Finite-frequency sensitivity kernels
for global seismic w a v e propagation based u p o n ad­
joint m e t h o d s[J] .Geophysical Journal International,
2008,174(1)：265-286.
[8] Tape C, Liu Q,T r o m p J. Finite frequency tomography u-
sing adjoint methods —methodology and examples using
m e m b r a n e surface waves [J] .Geophysical Journal Inter­
national,2007, 168(3)：1105-1129.
[9] Tarantola A. Inversion of seismic reflection data in the a-
coustic approximation [J] .Geophysics, 1984, 49 (8 )：
1259-1266.
[10] T r o m p J,T a p e C, Liu Q.Seismic tomography, adjoint
methods, time reversal, and banana-doughnut kernels
[J] .Geophysical Journal International, 2005, 160(1)：
195-216.
[11]李海山，杨午阳，雍学善.三维一阶速度一应力声波
方程全波形反演[J].石油地球物理勘探，2018,53(4):
730-736.
LI Haishan, Y A N G W u y a n g, Y O N G Xueshan. Three-
dimensional full w a v e f o r m inversion based on the
first-order velocity-stress acoustic w a v e equation [J].
Oil Geophysical Prospecting, 2018,53(4) ：730-736.
[12] Marquering H,D a h l e n F A,Nolet G. Three-dimensional
sensitivity kernels for finite-frequency traveltimes: T h e
banana-doughnut paradox [J] .Geophysical Journal Inter­
national, 1999, 137(3)：805-815.
[13] Dahlen F A, H u n g S H,Nolet G. Frechet kernels for
finite-frequency traveltimes，I: T h e o r y [J] •G e o p h y­
sical Journal International, 2000,141(1)：157-174. [14] H u n g S H,D a hlen F A,Nolet G. Frechet kernels for
finite-frequency traveltime, II：E x a m p l e s [J] .G e o -
第56卷第1期冯波，等：广义Rytov近似有限频初至层析99
phyvsical Journal International, 2000,141(1) ：175-203. [15] Zhao L, Jordan T H, C h a p m a n C H.Thr e e-d i m e n­
sional Frechet differential kernels for seismic delay
times [J]. Geophysical Journal International, 2000,
141(3)：558-576.
[16]崔超，黄建平，王自颖，等.基于双差的波动方程反射
波旅行时反演方法[J].石油地球物理勘探，2017,
52(4)：734-742.
C U I C h ao, H U A N G Jianping, W A N G Ziying, et al.
Double difference wave-equation reflection traveltime
inversion [J] .Oil Geophysical Prospecting, 2017，
52(4):734-742.
[17]付继有，李振春，杨国权，等.声介质波动方程反射旅
行时反演方法[J].石油地球物理勘探，2015,50(6):
1134-1140.
F U Jiyou.LI Z h e n c h u n,Y A N
G G u o q u a n.e t al. W a v e
equation reflection travel-time inversion in acoustic
media [J] .Oil Geophysical Prospecting, 2015,50( 6):
1134-1140.
[18]谢春，刘玉柱，董良国，等.基于声波方程的有限频伴
随状态法初至波旅行时层析[J].石油地球物理勘
探，2015,50(2) :274-282.
X I E C h u n, L I U Y u z h u,D O N G Liangguo, et al. First
arrival t o m o g r a p h y with finite frequency adjoint-state
m e t h o d based o n acoustic w a v e equation [J] .Oil G e o­
physical Prospecting, 2015,50(2)：274-282.
[19] W u R S. W a v e propagation, scattering and imaging u-
vsing dual-domain o n e-w a y a nd one-return propagators
[J]. P ure a nd Applied Geophysics, 2003, 160(3-4):
509-539.
[20] Spetzler J，Snieder R. T h e effects of small-scale he­
terogeneity o n the arrival time of w a v e s [J] .G e o­
physical Journal International,2001,145(3) ：786-796. [21] Spetzler J, Snieder R. T h e Fresnel v o l u m e and trans­
mitted w a v e s[J] .Geophysics, 2004,69(3) ：653-663. [22] Jocker J，Spetzler J，S m e u lders D, et al. Validation of
first-order diffraction theory for the traveltimes and
amplitudes of propagating w a v e s [J] .Geophysics,
2006,71(6) ：T167-T177.
[23] Xie X B, Y a n g H.T h e finite-frequency sensitivity kernel
for migration residual moveout and its applications in
migration velocity analysis [J] .Geophysics, 2008, 73(6)：
S241-S249.
[24]刘玉柱，董良国，李培明，等.初至波菲涅尔体地震层
析成像[J].地球物理学报，2009,52(9):2310-2320.
L I U Y u z h u,D O N G Liangguo, LI Peiming, et al.
Fresnel v o l u m e t o m o g r a p h y based o n the first arrival
of the seismic w a v e[J] .Chinese Journal of G e o p h y­
sics, 2009,52(9)：2310-2320.
[25] X u W,X i e X B, G e n g J. Validity of the Rytov approxi­
mation in the form of finite-frequency sensitivity kernels
[J]. Pure and Applied Geophysics, 2015, 172 ( 6 )：
1409-1427.
[26]冯波，罗飞，王华忠.一阶R y t o v近似有限频旅行时
层析[J].地球物理学报，2019,62(6) :2217-2226.
F E N
G Bo, L U O Fei, W A N G
H u a z h o n g.W a v e equa­
tion traveltime t o m o g r a p h y using R y t o v a p p r o x i m a­
t i o n^.Chinese Journal of Geophysics, 2019,62(6)：
2217-2226.
[27]刘定进，胡光辉，蔡杰雄，等.高斯束层析与全波形反
演联合速度建模[J].石油地球物理勘探，2019,
54(5)：1046-1056.
L I U Dingjin, H U Guanghui, C A I Jiexiong, et al. A
joint velocity m o d e l building m e t h o d with Gaussian
b e a m t o m o g r a p h y a nd full w a v e f o r m inversion for
land seismic data [J] .Oil Geophysical Prospecting,
2019,54(5)：1046-1056.
[28] B e y d o u n W B,Tarantola A.First Bo r n a n d R y t o v a p­
proximations ： Model i n g a nd inversion conditions in a
canonical e x a m p l e[J] .T h e Journal of the Acoustical
Society of America, 1988,83(3), 1045-1055.
[29]罗飞，王华忠，冯波，等.透射波旅行时B e a m层析成
像方法[J].石油物探，2019,58(3):356-370.
L U O Fei, W A N G H u a z h o n g, F E N G Bo, et al. B e a m
tomography based on transmission traveltime [J] .G e o­
physical Prospecting for Petroleum, 2019, 58 (3 )：
356-370.
[30] F e n g B, W u R S, W a n g H.A generalized R y t o v ap­
proximation for small scattering-angle w a v e propaga­
tion and strong perturbation media [J] .Geophysical
Journal International, 2019,219(2)：968-974.
[31] Snieder R,L o m a x A.Wavefield smoothing and the effect
of rough velocity perturbations on arrival times and a m­
plitudes [J] .Geophysical Journal International, 1996,
125(3)：796-812.
[32] L u o F, W a n g H,e t al. Automatic traveltime picking
using structure-enhanced high-dimensional time-fre­
quency transform [C] .E x t e n d e d Abstracts of 81st
E A G E Conference a nd Exhibition, 2019,1-5.
[33] Y u a n Y (), Simons F J,T r o m p J. Double-difference ad­
joint seismic tomography [J] .Geophysical Journal Inter­
national, 2016,206(3) ：1599-1618.
[34] Z h ang H.T h u r b e r C H.Double-difference to m o g r a p h y：
T h e method and i t v S application to the H a y w a r d Fault,
California [J] .Bulletin of the Seismological Society of
America,2003,93(5) ：1875-1889.
(本文编辑：宜明理）
作者简介
冯波博士，副研究员，1984年生；
2007年获同济大学地球物理学专业学
士学位；2015年获同济大学固体地球物
理学专业博士学位；现就职于同济大学
海洋高等研究院，主要从事地震波反
演、成像方面的理论与应用研究。

I V Oil Geophysical Prospecting2021
technique of overlapping group sparsity, and utilize the non-convex Lp-pseud〇-norm for preserving weak signals. O ur new model can fully exploit lo­cal similarity of signals instead of unrelated indi­vidual data. T o solve the m ulti-constrained p ro b­lem, we adopt the alternating direction m ethod of m ultipliers to divide the whole problem into four sub-problem s. T h e iteratively re-weighted altern a­ting direction L i-norm and majorization minimiza­tion algorithm are added into the algorithm to im- prove the efficiency and accuracy. Applied to syn­thetic and field data, the improved m ethod not o n­ly reduced staircase effects and attenuated random noises, but also effectively preserved weak signals. Keywords:overlapping group sparsity, higher-order total variation, non-convex Lp-pseud〇-norm, sta ir­case effects, random noise suppression
1. School of E arth and Space Sciences, Peking U- n iv ersity, Beijing 100871, China
2. Research Institute of Petroleum E xploration and
D evelopm ent, PetroC hina, Beijing 100083, China
A method of combining multi-channel signals to sup­press the strong reflection through matching pursuit. YANG Zipeng1,SONG W eiqi1, LIU Jun 2,CHEN Ju-nan2 ,L IU qun2 ,and WU D i1. Oil Geophysical Pros­pecting, 2021,56(1):77-85.
It is usually difficult to identify useful signals in and around stro n g reflections. Relying on the sparse representation ability of the m atching p u r­suit algorithm, we propose a m ethod for su p p ress­ing strong reflections th ro u g h matching pursuit (T C-M M P) under the constraint of m ultiple traces. T o determ ine the condition for term inating iteration in the process of signal decomposition, we introduce a residual ratio threshold into the m a­tching pursuit algorithm. Compared with fixed itera­tions or fixed residual threshold, this algorithm can improve the efficiency of signal decomposition. It is a m atching pursuit algorithm under the con­straint of a term ination condition. F irst, seismic signals w ith strong reflections are decomposed by T C-M M P, and m atching wavelets with the stro n­gest energy are searched and regarded as the best m atched strong reflections. T h en the weighted en­ergy coefficient is obtained based on the strong re­flections extracted from multiple traces. It can in­crease the energy proportion of effective signals in strong reflections. Finally, the optimal param eter for suppressing strong reflections is obtained through te sts, m aking the useful signals in strong reflections retained while the suppressed signals near the strong reflections highlighted. Applica­tions to a theoretical model with a high-resistance layer and raw seismic data with well tops show that the proposed m ethod is more reliable and ef­fective than conventional m ethods. It can provide technical support to accurate reservoir description. Keywords：m ulti-channel e x tra c t, suppress the strong reflection, energy w eighted, T C-M M P al­gorithm, useful signal
1. School of Geosciences, China University of Petrole­um (E^ast China), Qingdao, Shandong 266580, China
2. Research In stitu te of E xploration and Develop­m ent, N o rth w est Oilfield Branch Co. , SIN O PEC, U ru m q i, Xinjiang 830011, China
Deblending massive OBN data acquired by efficient and blended shooting method. CHEN Yingpeng1, ZHANG Hongjun1, LIU Yong1, ZHAO M in1, SONG Jiawen2, and QI Qunli1. Oil Geophysical Prospecting,2021,56( 1): 86-91.
Compared w ith land blended data with simple noises, OBN blended data have many types of noi­ses due to the special m ethod of OBN acquisition. T w o prim ary causes for the m ultiple types of noi­ses are heavy overlapped sh o ts on the logical coor­dinate system and severe geom etry deformation.
A fter analyzing the types of noises in OBN data, we studied iterative dynamic mapping, pre-pro­cessing of logical shot coordinates and sparse inver­sion deblending technology. T he sparse inversion deblending m ethod is based on F K K. It globally maps the logical positions of all overlapped shots in a dynamic and iterative way to recongnize noises, and then it accurately and quickly deblends the massive seismic data severely blended. T he appli­cation in Block A show s that the deblended OBN data are of high fidelity, indicating that the m ethod can provide optimal results while improving de­blending efficientcy.
Keywords ： efficient and blended acquisition of OBN data, noises, iterative dynamic m apping, pre-pro­cessing of logical shot coordinates, sparse inver­sion deblending
1. Overseas Business D epartm ent, G R I, BGP I n c.,
C N PC, Zhuozhou, Hebei 072750,China
2. Research &Development C enter, BGP I n c.,
C N PC, Zhuozhou, Hebei 072750, China
Wave-equation traveltime tomography using the ge­neralized Rytov approximation. FENG Bo1,2, WU Ru-Shan , LUO Fei1, XU Rongwei1, and WANG Huazhong1 . Oil Geophysical Prospecting, 2021, 56(1): 92-99.
T h e conventional finite-frequency tom ography
Vol. 56 No. 1Abstracts V
is usually derived from the Born or Rytov approxi­m ation which implies w eak-scattering assum ption. Therefore, the linearized forward problem in the finite-frequency theory is not satisfied for strong velocity p ertu rb atio n s. In the case of forward scat­tering and small-angle propagation, the generalized Rytov approxim ation (G R A)m ethod recently de­veloped can achieve im proving phase accuracy of forw ard-scattered wavefield, m aking it more suit­able for traveltim e tom ography. In this paper, we combine the conventional finite-frequency theory with G R A and propose a G RA-based traveltime sensitivity kernel, which w orks well regardless of the m agnitude of velocity perturbations. Numerical examples show that the traveltim e perturbation of forw ard-scattered waves can be correctly handled by the G R A-based traveltim e sensitivity kernel. T hen we propose an implicit m atrix-vector product strategy which can calculate the H essian m atrix- vector product w ithout explicitly forming the H es­sian m atrix, m aking it more attractive for 3-D problem s. W e solve the traveltim e inverse problem with the G au ss-N ew to n m ethod, where the H e s­sian m atrix-vector product is obtained by the pro­posed implicit m atrix-vector product m ethod. C on­sequently, the G auss-N ew ton m ethod can be rea­lized in a m atrix-free fashion, reducing the com pu­ter memory and disk occupancy significantly. N u­merical tests dem onstrated that the prop〇v S ed GRA- based traveltim e tom ography can estim ate the near- surface velocity model with high resolution and at a fast convergence rate.
Keywords ：finite-frequency traveltim e sensitivity kernel, generalized Rytov approxim ation, first-ar­rival traveltim e tom o g rap h y, G auss-N ew ton algo­rithm
1. Wave Phenom ena Intelligent Inversion Imaging
G r o u p(W P I)，School of Ocean and E arth Science, Tongji U niversity, Shanghai 200092 China
2. Institute for Advanced Study, Tongji U niversi­ty, Shanghai 200092 China
3. Modeling and Im aging L aboratory, University of California, Santa Cruz, California 95060, U SA
Low-amplitude structural imaging based on seismic numerical simulation and physical modeling. TIAN Yancan12, SHI Wenwu% WANG Guoqing1,2, XU Zhonghua1,2, and JIANG Chunling1,2. Oil Geophys­ical Prospecting, 2021, 56(1): 100-108.
E xploration cases show that low-amplitude structures can form small but high-yield reservoirs with favorable conditions of source-reservoir-cap and oil-gas m igration. H ow ever, traditional me­thods of velocity modeling cannot eliminate lateral variation of velocity usually caused by the varying proportion of sand in shallow form ations and p ro­duce a velocity model with low precision, giving rise to poor imaging precision of low-am plitude stru ctu res in M atouying area in Jidong oilfield with big drilling deviation. In this paper, the velocity modeling technology for low-amplitude structural imaging based on seismic numerical simulation and physical modeling is proposed. A simplified model that highlights the major characteristic in this area is designed for seismic numerical simulation and physical m odeling. Seismic simulation experim ents show that correct imaging results can be obtained w ith the real velocity. A t the same time, real ve­locity and stratigraphic form are used as the prior information for velocity modeling of low-am plitude stru ctu res to validate the technology. T h e results show that the tom ographic inversion along th e for­m ation can create a precise initial velocity model, and the high-density structure-constrained grid to­mographic imaging technology can build a highly precise model that indicates the trend in velocity. T herefore, the joint application of tw o modeling technologies leads to more precise imaging of low- am plitude stru c tu re s.
Keywords：
physical simulation, numerical sim ula­tion, low am plitude stru ctu re, velocity modeling, grid tom ography
1. Key Laboratory of Reservoir Description, CNPC, Lanzhou, Gansu 730020，China
2. Research In stitu te of Petroleum E xploration Development -N o rth w est, PetroC hina, Lanzhou,
G ansu 730020，China
3. E xploration and Development Research In stitu­te, Jidong Oilfield C om pany, PetroC hina, Tang- shan, Hebei 063000, China
Research on modeling and imaging of deep low-am­plitude structures in the northwestern Sichuan Ba­sin. WANG Yanxiang1, SU Qin1, LE Xingfu1, ZHANG Junduo1, LIU Wei1, and YUAN Huan1. Oil Geophysical Prospecting, 2021, 56 ( 1): 109-117.
T o disclose the details of stru ctu ral trap s in the Longm enshan piedmont belt in the n o rth w e s t­ern Sichuan Basin, we carried out imaging and modeling of deep low-amplitude stru c tu re s, mesh- tom ography-based velocity m odeling, and physical sim ulations of stru c tu re s. L ongm enshan piedmont belt presents severe surface fluctuations and seri­ous problem s of static correction. T h e n adaptive first-arrival tom ographic inversion w ith m icro-log­。