0 引言

根据地震波在地下介质中传播速度的差异，通常可以利用速度建模来还原地下的地层构造，如速度分析、层析成像和全波形反演等.但在实际的建模过程中，无可避免地会引入“被光滑的”速度界面(下文中称这样的界面为“软界面”，常规的速度界面为“硬界面”)，例如迭代反演中每次更新后的速度模型(王连坤等, 2016; 张兵和王华忠, 2019; 张盼等, 2019)，其对下次模拟数据的影响势必有别于硬界面模型.此外，尽管多数情况下地质体与地质体之间有着截然不同的速度，但当两者之间存在过渡带时(Gómez and Ravazzoli, 2012; Han et al., 2012; 韩复兴等, 2014)，软件界反射也成为不得不考虑的问题.软界面本身并不是一个新的概念，关于其更多的描述来源于过渡层反射系数(Berryman et al., 1958; Brekhovskikh, 1980; Kuhn, 1988; Robins, 1990).区别于硬界面，由于层厚引起的传播效应，软界面反射系数(下文未经特殊说明的，默认为平面波反射系数)因此具有频变性质(Meyer, 1975; Wolf, 1987).以最简单的单界面反射为例，界面光滑前后反射系数的这种变化将是光滑处理对反射信号影响的最直接体现.

前人为了模拟反射信号做了大量工作(Wiggins and Helmberger, 1974; Chapman, 1981; Porter and Reiss, 1984; Červený, 1989; Verweij and De Hoop, 1990)，但对于介质光滑以后反射信号是否发生变化的问题却鲜有关注.一方面是因为后者涉及到连续介质(Brekhovskikh, 1980)下反射波场的理论计算，其本身就是一个亟待解决的难题；另一方面是因为许多反演方法迫切需要的是初始模型，目标函数及迭代算法的有效性(刘璐等, 2013; 王连坤等, 2016; 邢贞贞等, 2019)，而非模拟数据的十分精准.关于光滑处理对透射波场影响分析的工作已然在偏移成像领域得到体现(Gray, 2000; Jones et al., 2010; 韩复兴等, 2015)，随着勘探技术的进一步发展，光滑处理对反射波场的影响也成为不可回避的问题，特别是对于众多模拟算法(如有限差分法，伪谱法和有限元法等)在高精度计算中的校正和推广具有借鉴意义.近几年来，连续介质或过渡层反射的问题多数以散射的方式进行描述(Fuck et al., 2011; Ross and Lavery, 2012; 李辉和王华忠, 2015; Gibson, 2018; 李懿龙, 2019)，其在实际的数值模拟过程中与有限差分同样存在着一定程度的数值模拟误差(张永刚, 2003; 严红勇和刘洋, 2013)，而由于光滑处理引起的介质差异往往是极其微弱的，这就需要较为精确的理论计算.文中以平面波作为入射波场，利用过渡层反射系数推导了反射信号计算的理论公式，并通过反射系数方程控制反射系数，数值求解了软界面下的反射系数，与经典Epstein层的对比说明了算法的精度，最终的实验分析揭示了速度光滑对反射信号到时、主频与能量的不同影响.

1 方法原理 1.1 反射信号表示

设平面波只沿深度z传播(零炮检距)，波场记为φ(z, t)，其与谱密度函数Φ(z, ω)的关系由下述Fourier变换决定，

(1)

(2)

其中，i为虚数单位，ω为角频率.

将波场φ(z, t)分解为入射波场p(z, t)与反射波场r(z, t)，对应的谱密度函数记为P(z, ω)和R(z, ω).定义深度z处的反射系数(Goupillaud, 1961; Brekhovskikh, 1980)：

(3)

如果传播介质均匀，则在接收点z₀处的反射信号为

(4)

计算中深度z₀处入射波场的谱密度函数P(ω; z₀)通常是已知的，问题的关键是反射系数的计算.均匀介质中，根据定义(3)，深度z处反射系数与深度z′处反射系数具有如下关系(z < z′，并且取波场传播因子为e^-ikΔ，Δ>0)：

(5)

其中，k₀=ω/v₀，表示均匀介质中的波数.

考虑图 1中的硬界面与软界面速度变化曲线(黑色实线对应硬界面，红色虚线对应软界面)，其中硬界面深度为z₂，设其上层半无穷均匀介质速度与密度分别为v₀与ρ₀，下层半无穷均匀介质速度与密度分别为v₁与ρ₁，则有反射系数：

(6)

根据公式(4)和(5)，硬界面下的反射信号被表示为

(7)

为了计算软界面下的反射信号，引入了参考深度的概念(图 1中黑色虚线)，其位于软界面过渡层外上部半无穷均匀介质中，具有深度z₁，如果该深度的反射系数V(ω; z₁)可算，则软界面下的反射信号被表示为

(8)

为了具体分析速度光滑对反射信号到时与能量的影响(与硬界面相比)，分别定义相对到时与相对能量，

(9)

(10)

公式(9)中的t_peak^s(r^s)和t_peak^h(r^h)对应表示软界面反射信号r^s(t; z₀)与硬界面反射信号r^h(t; z₀)峰值的到时，公式(10)中R^s(ω; z₀)和R^h(ω; z₀)对应表示软界面反射信号与硬界面反射信号的谱密度函数.

下一部分将就参考深度处反射系数的计算问题展开讨论，实际为了计算稳定，参考深度应尽可能地接近软界面过渡层的上限深度.

1.2 反射系数计算

根据反射系数的定义，软界面过渡层外上部半无穷均匀介质中，反射系数具有递推公式(5)；而软界面过渡层外下部半无穷均匀介质中，由于没有反射波场的存在，反射系数为零.因此软界面下反射系数计算问题(零炮检距)可以表述成如下的反射系数方程(推导见附录)及相应的边界条件：

(11)

(12)

其中，波数k=ω/v，ρ为密度，z_bot为软界面的下限深度.

设反射系数具有如下形式的解(Brekhovskikh, 1980)：

(13)

其中，u(z)和w(z)为两个任意的函数.

将方程(13)代入到方程(11)中得到，

(14)

方程(14)及边界条件(12)成立的一个充分条件可以写成方程组

(15)

及相应的边界条件，

(16)

不妨取

(17)

为了简化方程组(15)的表示，引入替换，

(18)

其中，n=v₀/v表示折射率.从而可将方程组(15)改写为

(19)

其中，k₀=ω/v₀，表示上半无穷均匀介质中的波数.在密度恒定的情况下，η₁和η₂具有表达式：

(20)

(21)

接着对方程组(19)两端同时进行积分变换，

(22)

并考虑边界条件(17)，有

(23)

得到的w(z)和u(z)继续代入方程组(19)的右端，如此反复，第k+1次迭代结果与第k次迭代结果之间的关系被写成下述求和表达式：

(24)

迭代格式(24)的收敛性很容易通过其级数解收敛的性质来证明(Brekhovskikh, 1980).这样，参考深度z₀处第k+1次迭代得到的反射系数即为

(25)

为了使计算结果稳定，迭代格式(24)的计算范围[z, z_bot]不宜过大，所以计算中应尽可能地让参考深度z₀接近软界面的上限深度.

2 模型测试

为了说明迭代算法的精度，在本小节测试经典的Epstein过渡层(Epstein, 1930)，其速度变化满足下式：

(26)

其中，m与N为两个特定参数.

在距离过渡层足够远的入射端区域，反射系数(复数)的模大小表示为(Brekhovskikh, 1980)

(27)

其中，sinh表示双曲正弦函数；S≡2k₀/m, 表示过渡层的相对厚度.

如图 2所示绘制了3条不同参数m的Epstein过渡层速度剖面，图 3中的3条实线则表示各自反射系数模随频率变化的真实曲线，对应的3条虚线为40次迭代计算结果，可以看到，两者几乎完全吻合，只在相对较高的频率时迭代结果出现不稳定(图中未展示相应的计算结果)，鉴于真实解随着频率增加迅速衰减为零，利用这一特性，高频迭代解很容易得到修正，这里不再赘述.

3 实验分析

为了分析速度光滑对反射信号到时、主频与能量的影响，本小节以单界面为例进行计算对比.如图 4所示，单界面距源点(或接收点)的深度为400 m，其上层均匀介质速度为1500 m·s^-1，下层均匀介质速度为2000 m·s^-1，对速度剖面按照深度2 m的网格间距采样，并分别进行不同次数的3点光滑.假设震源子波为Ricker子波，其时间域表达及谱密度函数分别为(俞寿朋，1996)

(28)

与

(29)

其中，f=ω/(2π)，表示频率；f_M为Ricker子波主频，计算中取为20 Hz.

表 1给出了不同光滑次数下反射信号主频、相对到时与相对能量的计算结果，可以看到反射信号的到时几乎不受光滑次数的影响，而信号主频与能量随着光滑次数的增加而减少.图 5与图 6分别绘制了反射信号能量与主频随光滑次数的变化曲线，其中信号能量在低光滑次数的衰减速率明显大于高光滑次数，而信号主频随光滑次数的增加近似呈线性递减.为了具体展示光滑算子对速度剖面、反射信号及信号频谱的影响，分别在图 4、图 7及图 8中给出了4次、10次、20次和30次光滑后的结果.

4 结论与讨论

与硬界面不同，软界面反射系数是频变的，通常具有低通滤波的性质，并且两者在平面波反射信号上存在差异，其主要表现在光滑算子对信号能量的衰减和对信号主频的降低上，通过计算可以看出，只需较少的光滑次数，信号能量就被极大地衰减，而随着光滑次数的进一步增加，能量变化趋势又逐渐变缓.尽管光滑算子对信号能量和主频的影响是明显的，但信号到时却几乎不受影响，这也保证了实际数据处理中许多光滑处理的可行性.

文中只讨论了零炮检距下速度光滑对平面波反射信号的影响特点，对于非零炮检距，界面光滑前后，相同位置相同出射角度的平面波必然会被不同的观测点所接收，而相同观测点接收到的平面波信号又会有着不同的出射角度，因此在反映光滑处理的影响中意义不大.此外值得注意的是，在不考虑介质吸收的情况下，根据能量守恒原则，反射能量的减弱将导致透射能量的增强，这对于利用透射波场信息进行成像的一些方法有何借鉴意义将值得我们进一步的思考与探索.

附录  反射系数方程

非均匀介质中的声波方程可以通过连续性方程和Euler方程得到(Brekhovskikh, 1980)：

(A1)

其中，φ为波场(压力场，不区分压力和压强的概念)，c为质点振动速度矢量，ρ和v分别为介质的密度与速度，▽·和▽分别为散度与梯度算子.

只考虑沿深度z方向传播的波场，振动速度矢量c表示为沿深度的振动速度c，方程组(A1)被改写为

(A2)

对方程组(A2)两端同时进行公式(1)所述的Fourier变换得到

(A3)

其中，Φ(z, ω)表示波场φ(z, ω)的谱密度函数，其可被分解为入射谱P(z, ω)和反射谱R(z, ω)，对应的质点振动速度C(z, ω)分别为

(A4)

和

(A5)

将入射谱P(z, ω)和反射谱R(z, ω)，以及方程(A4)和(A5)代入到方程组(A3)，有，

(A6)

(A7)

根据反射系数的定义(3)，对方程(A6)和(A7)进行变形消项得到，

(A8)

其中，

方程(A8)即为连续介质的反射系数方程(零炮检距)，其具有Riccati方程的形式.