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

冯 波 WU Ru-Shan 罗 飞 许荣伟 王华忠

(①同济大学海洋与地球科学学院，上海 200092；②同济大学海洋高等研究院，上海 200092；③Modeling and Imaging Laboratory，University of California，Santa Cruz，CA 95060，USA)

0 引言

地震波旅行时常用于反演地下速度结构。在射线层析中，由于引入了高频近似假设，旅行时只与连接炮检对的几何射线上的速度扰动有关。对于有限频带的地震波，若速度非均匀体的尺度小于菲涅耳体的宽度时，地震波散射起主导作用。由此发展了绕射层析理论[1-2]，可以较好地处理地震波的一阶绕射效应。基于Born近似， Luo等[3]给出了地震波旅行时Fréchet导数(旅行时敏感度核函数)，可以定量描述速度扰动对地震波旅行时扰动的影响。通过波动方程线性化近似(Born近似或Rytov近似)，Woodward[4]提出了“波路径”的概念，作为无限高频假设下的射线路径向有限频地震波传播的推广。借助伴随状态理论，可以实现(波形、旅行时、振幅等)Fréchet导数的高效计算[5-11]。

通常旅行时Fréchet导数由Born近似导出[12-18]，但Born近似成立条件较为苛刻。相比于Born近似， Rytov近似可以更好地描述由于速度扰动引起的前向散射波的相位扰动[19]，因而更适用于旅行时反演。基于Rytov近似也可以构造相应的旅行时敏感度核函数[20-27]。

然而，Born近似与Rytov近似均属于弱散射近似，要求速度扰动远小于背景速度[28-29]。为克服弱散射近似的制约， Feng等[30]提出了可以适用于强扰动模型的广义Rytov近似理论(Generalized Rytov approximation，GRA)。对于前向散射及小角度传播，GRA可以较为准确地预测地震波的相位扰动，因而更适用于旅行时反演。本文将GRA理论应用于传统的有限频层析，构造一种适用范围更广的旅行时敏感度核函数并用于初至旅行时层析。模型数据测试结果表明，本文提出的广义Rytov近似旅行时层析可以建立高精度的近地表速度模型，且收敛速度较快。

1 理论

1.1 广义Rytov近似有限频旅行时敏感度核函数

频率域标量声波方程可以表示为

(2+k2)u0=0

(1)

式中：k(x)=ω/v0(x)，为背景速度场v0(x)中的地震波数；u0(x，ω)为背景波场；2为Laplace算子。

记扰动后的速度场为v(x)，波场(总场)为u(x，ω)，总场与背景场的关系为

u(x,ω)=u0(x,ω)exp[ψ(x,ω)]

(2)

式中ψ(x，ω)为复相位。

在前向散射(散射角定义如图1所示)及小角度传播假设下，复相位可以用广义Rytov近似[31]表示为

图1 散射角定义：地震波入射方向与散射方向的夹角前向散射对应小散射角，背向散射对应大散射角

(3)

式(3)引起的相位延迟可以表示为

(4)

有限频带地震信号的旅行时扰动可以表示为单频谐波的相位延迟的加权叠加[23-26，32]，有

(5)

结合式(3)～式(5)，可以建立带限信号旅行时扰动与慢度扰动的线性关系

(6)

式中：xs和xr分别为源和检波点坐标；KGRA(x;xr,xs)为广义Rytov近似旅行时敏感度核函数(Generalized Rytov approximation based traveltime sensitivity kernel， 简称为GRA核函数)，可表示为

(7)

显然，显式计算式(7)需要计算频率域格林函数[25]，因而计算量较大，尤其是对于三维问题。为了避免显式计算核函数，本文采用冯波等[26]给出的方法，通过引入如下辅助函数

(8)

利用Parseval定理，将式(7)中对角频率的积分转化为对时间积分，则

KGRA(x;xr,xs)

G0(x,ω;xr)]dω}

(9)

式中：T表示地震记录时长；g0表示时间域格林函数；“*”表示褶积(与上标“*”表示复共轭不同)。

若定义旅行时伴随场为

(10)

式(9)可以重写为

KGRA(x;xr,xs)=

(11)

1.2 旅行时层析反问题求解

基于旅行时残差L2范数的误差泛函可表示为

(12)

式中：m=(m1,m2,…,mNv)T∈M，为模型空间M中的向量，本文为慢度模型(采用网格参数化方式)，其中Nv为网格点数；Δt=(Δt1,Δt2,…,ΔtNd)T∈D，为数据空间D中的向量(旅行时残差)，其中Nd为炮检对数。

利用Gauss-Newton反演算法，式(12)可以迭代求解

(13)

式中：k为迭代次数；α为步长，可以通过线性搜索方法计算；J为目标函数梯度；为线性Hessian矩阵(忽略旅行时对模型参数的二阶导数)，其中K为KGRA的矩阵形式；P为光滑算子。

根据式(12)，泛函梯度可表示为

(14)

基于广义Rytov近似，∂Δt/∂m即为式(11)给出的核函数。将式(11)代入式(14)，则泛函数梯度可表示为

J(mk)(x)=(KTΔtk)(x)

(15)

考虑到显式计算Hessian矩阵需要巨大的计算量及存储量，且求逆困难，因此，将对Hessian矩阵的求逆转化为求解如下法方程的近似解

Ha(mk)Δmk=J(mk)

(16)

本文用共轭梯度(Conjugate gradient，CG)法求解式(16)获得模型更新方向，并采用隐式矩阵向量乘方法[26]直接计算Hessian矩阵与向量的乘，因而无需显式计算及存储Hessian矩阵。具体计算公式详见附录A。

2 数值试验

2.1 高斯扰动模型

图2为一个含有高斯异常体的二维10km×5km速度模型，网格间距为10m。其背景速度v0=3000m/s，异常体的速度扰动可表示为

图2 含有高斯扰动(ε= 50%)的速度模型及观测系统红色倒三角为检波点位置

(17)

首先，用时—空域高阶有限差分求解声波方程，用真实模型(包含高斯异常体的速度模型)和背景模型(即均匀背景速度模型)进行正演(震源是主频为15Hz的Ricker子波)，得到的地震记录分别作为观测数据和模拟数据。然后，用互相关(Cross-correlation，简记为CC)函数测量观测数据和模拟数据的时差，作为真实时差的测量值(即观测时差)。随后，分别用射线(Ray)核函数(速度网格内的射线长度)、Rytov近似旅行时敏感度核函数[26]及本文提出的GRA旅行时敏感度核函数预测旅行时扰动(即核函数与模型扰动做内积)。图3为震源在地表(5.0km，0)和(0，0)处、201个检波器均匀分布在z=5.0km处(图2中红色倒三角所示)旅行时扰动对比。由图3可以看出：①射线和GRA核函数的预测时差完全重合，但与Rytov近似预测结果有较大偏差。②对于小角度前向散射(图3a中对应为：震源位于地表(5.0km，0)，检波器横坐标在(5.0km，5.0km)附近)，GRA核函数的预测时差与观测时差基本一致。证明了GRA核函数在前向小角度散射时具有较好的线性特征，克服了传统波动方程线性化近似中的速度小扰动假设的制约。这与GRA理论预测相符(GRA理论指出，对于小角度前向散射，GRA可以较为精确地预测地震波的相位扰动)。③随着散射角增大，GRA核函数的预测时差逐渐偏离观测时差，并逐渐与Rytov近似预测结果趋于一致。根据GRA理论推导可知[30]，Rytov近似为GRA的一阶近似，本质上都隐含了小角传播假设。因此，随着散射角增大，GRA的成立条件(即前向小角散射)不再成立，导致预测时差偏离观测值。对于图3b中的震源位置(左上角(0，0)处)，前向小角度散射对应检波器横坐标在10km附近，三类核函数预测结果与图3a一致。

图3 震源在地表不同位置三种旅行时敏感度核函数的预测时差与测量时差对比(a)(xs，zs)=(5.0km，0)；(b)(xs，zs)=(0，0)

上述测试表明：GRA核函数精度与散射角密切相关。对于前向小角度散射，GRA可以较为准确地预测地震波的旅行时扰动。随着散射角不断增大，GRA核函数精度逐渐降低，但仍好于Rytov近似。虽然对于大角度散射，GRA核函数的预测结果并不准确，甚至与观测时差偏离较远，但并不会导致层析反问题求解失败。原因在于：层析反问题通过最优化算法迭代求解时差目标函数，而非直接对线性方程组直接求逆。即使正问题的线性化近似不严格成立，只要反问题的凸性较好，依然可以采用最优化算法求解目标泛函。

为了验证GRA初至层析方法的有效性、尽可能消除观测系统对反演结果的影响，设计了一个理想观测系统如图4a所示，100个检波器均匀分布高斯异常体模型的四周，围绕模型四周激发100次，震源是主频为15Hz的Ricker子波(对应波长λ0=0.2km，与高斯异常体尺度参数关系为a= 5λ0)。

图4 高斯模型观测系统(a)及GRA旅行时反演结果(b)黑色正三角形表示震源位置；红色倒三角形表示检波点位置

反演采用的初始模型为均匀背景速度场(v0=3000m/s)，并用互相关方法测量观测数据和模拟数据的时移量作为初至旅行时残差。反演采用SEIS-COPE数值优化软件包(SEISCOPE optimization toolbox)中的截断牛顿算法(对于旅行时层析，截断牛顿算法退化为高斯牛顿算法)实现最优化，其中，梯度和Hessian矩阵与向量乘由附录A中给出的隐式矩阵与向量乘方法计算。在迭代反演中，采用2次CG内迭代求解式(16)，获得模型更新方向，并用线性搜索方法计算迭代步长。反演收敛准则为判断相对旅行时误差(J(mk)/J(m0))是否小于预设的阈值(σ=1.0×10-4)。经过10次迭代(图5为反演收敛曲线)，反演结果如图4b所示，与真实速度模型在视觉上几乎没有差别。为定量对比反演结果，对反演的扰动模型(反演结果减去初始模型)与真实的高斯异常体进行横向和纵向抽线对比(图6)，反演的速度扰动与真实的高斯扰动基本吻合，证明了对于完备的观测系统及简单高斯速度模型，GRA初至层析可以得到高精度的反演结果，证实了方法的有效性。

图5 GRA初至层析的收敛曲线

图6 高斯模型层析(蓝色)与真实(红色)速度扰动曲线(a)x=2.5km；(b)z=2.5km

2.2 Marmousi-2模型

应用Marmousi-2 模型验证GRA初至层析对复杂速度模型及不完备的观测系统的有效性。为模拟陆上地震采集，切除了原始速度模型中上覆的水体部分，并在模型右侧扩边3km(图7a)。x和z方向网格数目分别为2001和301，网格间距均为10m。

图7 Marmousi-2模型(a)及GRA初至层析的速度更新量(b)

设计了一个单边观测系统：震源和检波器均放置在地表，起始炮点位于x=0，炮间距为100m，共171炮；最小和最大炮检距分别为400m和5000m，道间距为20m(每炮共231道)。最大炮检距较大，因此可以用初至波旅行时反演近地表速度。采用有限差分求解波动方程正演地震信号，震源是主频为15Hz的Ricker子波。由于观测地震记录较为复杂，采用自动方法[31]拾取初至，作为观测记录的初至旅行时。

初始模型采用随深度线性递增的速度场(z=0，v=1000m/s；z=3.5km，v=3000m/s)。为降低计算量，采用20m×20m网格剖分进行速度反演，并用相同的初至拾取方法拾取模拟数据的旅行时，进而计算初至时差。反演实现框架与高斯模型一致(相对旅行时误差阈值σ=1.0×10-2)。经过5次迭代反演收敛，得到的速度更新量如图7b所示，浅层的高速异常得到了较好的恢复。其收敛曲线如图8所示，可以看出，高斯—牛顿算法收敛速度较快。

图8 GRA初至层析的收敛曲线

为定量对比反演结果，抽取z=0.1、0.2、0.4、0.6km处速度横向变化曲线，如图9所示。可以看出，在近地表(z=0.1、0.2km)，反演结果(蓝色曲线)逼近真实的高速隆起，反演精度较高(图9a和图9b)。随着深度增加，反演精度逐渐降低(图9c和图9d)，但反演结果仍然与真实速度模型的整体趋势相吻合。

图9 不同深度处初始(绿色)、真实(红色)、GRA层析(蓝色)速度横向变化曲线对比(a)z=0.1km；(b)z=0.2km；(c)z=0.4km；(d)z=0.6km

Marmousi-2模型实验结果表明，利用中、小炮检距的初至波旅行时信息，GRA初至层析可以建立高精度的近地表速度模型。

3 讨论

旅行时(差)的精确测量是所有旅行时反演类方法所共有的问题。本文的重点为讨论如何将广义Rytov近似方法与有限频层析理论的结合，是笔者过去研究工作[26]的理论推广和应用，因而没有讨论如何准确测量初至波时差。若震源子波已知且模拟数据与观测数据的初至波形基本一致，此时测量初至波时差相对简单。反之，若震源子波未知或由于地下介质复杂导致初至波形畸变严重，此时需要采用与震源无关的反演类算法消除子波波形畸变对旅行时测量的影响，如双差旅行时(double-difference traveltime)[33-34]测量等。

4 结论

本文提出了一种基于广义Rytov近似的有限频旅行时敏感度核函数，克服了传统有限频层析方法中弱速度扰动假设的制约。对于前向散射及小角度传播，广义Rytov近似旅行时敏感度核函数可以较为准确地预测地震波的旅行时扰动，从而可以提高反演精度并加速反演收敛。在计算方面，本文采用了隐式矩阵向量乘方法直接计算Hessian矩阵与向量的乘，无需显式计算及存储Hessian矩阵，进而可以提高高斯—牛顿迭代反演算法的效率。复杂的Marmousi-2模型测试结果表明，GRA初至层析可以建立高精度的近地表速度模型。

附录A 隐式矩阵向量乘

Hessian矩阵与向量积可表示为

Hap=KTKp=KTh

(A-1)

式中：p=p(x)为模型空间中的向量；h=h(xr,xs)为数据空间中的向量。因此需要计算两次矩阵向量积。

首先， ∀(xr,xs)，Kp可以表示为

[G0(x,ω;xr)q*(xr,ω;xs)]dω}dx

(A-2)

式中up(xr,ω;xs)满足如下Born近似

(A-3)

利用Parseval定理，式(A-2)可改写为

(A-4)

式中

uq(xr,t;xs)=2πRe[q*(xr,t;xs)]

=2πRe[FT-1q(xr,ω;xs)]*

(A-5)

up(xr,t;xs)=FT-1up(xr,ω;xs)

(A-6)

同时，式(A-3)在时—空域满足如下方程组

(A-7)

式中f为震源函数。

其次，根据式(15)，对于空间任意一点x，KTh可以表示为

(A-8)

式中

(A-9)

为单炮旅行时伴随波场。