基于时域加权的拉普拉斯—频率域弹性波全波形反演

刘张聚 童思友 方云峰 贾君莲

(①中国海洋大学海底科学与探测技术教育部重点实验室，山东青岛 266100；②青岛海洋科学与技术试点国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266061；③浙江大学地球科学学院，浙江杭州 310027；④自然资源部第二海洋研究所海底科学重点实验室，浙江杭州310012；⑤东方地球物理公司物探技术研究中心，河北涿州 072751；⑥东方地球物理公司研究院，河北涿州 072751)

0 引言

全波形反演(Full Waveform Inversion，FWI)是利用地震记录的全波形信息，通过匹配模拟记录与实际记录反演地下介质参数的一种高分辨率建模方法。Tarantola[1-2]最先提出了时间域全波形反演方法，该方法基于广义最小二乘理论，使用理论记录与实际记录的残差来构建目标函数。此后，全波形反演得到了全面的发展。Bunks等[3]提出了时间域多尺度反演，降低了反演的非线性影响；Pratt[4]提出了在频率域进行全波形反演，不需要对地震记录进行预处理，且使用有限个频率信息即可得到高精度的反演结果；Sirgue等[5]综合时间域反演和频率域反演的优势，提出了时间—频率混合域反演方法，避免了因时间域求梯度所需要的震源波场重构以及频率域正演所需计算机内存过大等问题。

另外，还有一种研究思路是降低全波形反演的非线性，降低陷入局部极小值的可能性，从而在缺失低频或者初始模型不准确的情况下，反演依然能逼近正确结果。Choi等[28]提出的基于全局互相关全波形反演可有效降低反演的非线性，其强调的是记录之间整体的相关程度，并非是记录中点与点的匹配程度；Liu等[29]详细对比了全局互相关目标函数与L2范数目标函数之间的差异，说明了互相关目标函数的优势；Tao等[30]将全局互相关应用到了全波形层析，其反演结果同样优于基于旅行时的层析成像；李青阳等[31]将互相关与反射波波形反演相结合降低了反演的非线性。Luo等[32]建立了反褶积目标函数，与互相关目标函数类似，降低了陷入局部极小值的可能性，但相对而言对带限或非脉冲源的灵敏度降低。Warner等[33-35]建立了自适应波形反演理论，通过将维纳滤波器引入目标函数，寻找将两种记录相匹配的滤波器，再建立滤波器的目标函数，使其逐渐逼近零延迟的单位脉冲函数；Guasch等[36]将该方法推广到三维，并且在缺失低频的情况下得到了较好的反演结果。van Leeuwen等[37]提出了波形重构反演，将波动方程本身作为惩罚项引入到目标函数中，通过重构波场更新模型参数，将解空间扩展到了数据空间和模型空间，有效减小了局部极小值的影响；Li等[38]将该方法引入时间域，提出了时间域波形重构法，相比于频率域可有效减少计算成本。

针对如何在缺失低频信息的情况下建立较高精度的初始模型这一难题，本文建立了一种基于时域加权的拉普拉斯—频率域弹性波全波形反演方法(Time-Weighted Laplace-Frequency Domain Elastic Full Waveform Inversion，TWLF-EFWI)。该方法在时间—频率域弹性波全波形反演的基础上，通过引入拉普拉斯阻尼因子，降低了全波形反演对于初始模型的依赖程度，形成了时间—拉普拉斯—频率域的混合域全波形反演方法。TWLF-EFWI方法仍使用传统的L2范数作为目标函数，为消除因引入拉普拉斯阻尼因子带来的地震道能量不均衡、远炮检距地震道能量衰减严重的影响，在目标函数中引入了时域加权因子，消除了拉普拉斯阻尼因子的负面影响。Marmousi 2部分模型测试结果表明，在缺失低频信息以及初始模型精度较低的情况下，本文方法可得到含有模型宏观构造的反演结果，足以满足后续反演对初始模型的精度要求。

1 全波形反演基本理论

1.1 频率域弹性波全波形反演

全波形反演可以描述为一个基于地震全波场模拟的数据拟合过程，目的是实现正演记录和实际观测记录的完全拟合，也可以将全波形反演重新转换为线性化的最小二乘问题。

在频率域，波动方程矩阵形式可以表示为

AU=S

(1)

式中：A为阻抗矩阵；U为全部网格点的波场值所组成的列向量；S为震源项。

基于广义最小二乘的全波形反演目标函数定义为

(2)

式中：m为模型参数；Ns为总炮数；Δd=dcal-dobs，为模拟地震记录dcal与实际地震记录dobs的差；上标“*”表示矩阵的伴随形式。全波形反演就是从目标函数出发，求解模型的更新量，然后通过对模型的迭代更新，使目标函数达到最小，或者使前后两次的目标函数的相对误差达到限定值，进而得到与其对应的模型参数的最优解。

将式(2)中的目标函数对m求导，可以导出目标函数在频率域关于模型参数的梯度为

(3)

式中J为Fréchet导数矩阵，可由波动方程两边同时对参数m求导得到

(4)

式中R表示波场U到接收记录d的映射关系。则频率域梯度可表示为

(5)

上式表明目标函数的梯度由正演波场U、阻抗矩阵对模型参数的导数∂A/∂m和以地震记录残差共轭为震源的正演波场A-1RTΔd*三部分组成。

对于二维弹性波方程，在频率域表达式为

(6)

式中：Ufx和Ufz分别为频率域速度波场的水平和垂直分量；Fx和Fz为对应的频率域震源项；ω为角频率；ρ为密度；λ和μ为拉梅常数。

联合式(5)和式(6)，可得目标函数关于拉梅常数的梯度在频率域的表达式

(7)

(8)

式中vP和vS分别为纵波速度和横波速度。

1.2 拉普拉斯—频率域弹性波全波形反演

拉普拉斯—频率域的二维弹性波方程与频率域的二维弹性波方程极其相似，为

(9)

式中：ULx和ULz分别为拉普拉斯—频率域波场的水平和垂直分量；Lx和Lz为对应的拉普拉斯—频率域震源；s=σ+iω，其中σ为拉普拉斯衰减常数。

对比式(9)与式(6)可以发现，二者的唯一区别是ω与s，所以拉普拉斯—频率域也可以看作是“复频率域”。两个域的目标函数梯度也具有相同的形式，在拉普拉斯—频率域具体形式为

(10)

同样，根据式(8)，可得目标函数关于纵波速度和横波速度的梯度。

1.3 混合域弹性波全波形反演

全波形反演在时间域和频率域里的基本理论是相通的，二者并无差异。时间—频率混合域全波形反演就是在时间域模拟正、反传波场，然后转换到频率域计算梯度。两种域的结合只需使用离散傅里叶变换(DFT)即可实现，即在时间域求取波场时，在时间循环中加入DFT，当循环结束时也就得到了对应的频率域的单频波场，之后再用频率域单频波场计算梯度。混合域反演结合了时间域模拟波场的高效性和频率域计算梯度的灵活性，其流程如图1所示。

DFT和逆离散傅里叶变换(IDFT)的表达式为

(11)

式中：u(kΔt)和UF(ωn)分别为时间域波场和频率域波场；Δt为时间采样间隔；NT为时间域序列长度；NF为频率域序列长度。

图1 时间—频率域全波形反演流程

在上述混合域的基础上，使用拉普拉斯—频率域代替频率域，形成了时域—拉普拉斯—频率域的全波形反演(图2)。与图1流程不同的是，需要将DFT改为离散拉普拉斯变换(DLT)进行时域与拉普拉斯—频率域间的转换，并且在求取伴随震源时需要做一次DLT和逆离散拉普拉斯变换(IDLT)。另外本文将伴随震源在时域中进行了加权处理，后文会对此进行详细描述。DLT和IDLT的表达式为

(12)

式中：UL(sn)为拉普拉斯—频率域波场；NL为拉普拉斯—频率域序列长度。

图2 时间—拉普拉斯—频率域全波形反演流程

2 时域指数衰减波场

相对于傅里叶变换，拉普拉斯变换只是增加了一个时间指数衰减项。正是这一指数衰减项，在拉普拉斯—频率域反演降低了对低频信息的依赖。

时间域信号u(t)的拉普拉斯变换为

(13)

将拉普拉斯变换中的衰减项单独取出，在时间域分析该项对信号频率的影响。时间域信号衰减可表示为

(14)

在时间域的乘积等价于频率域的卷积，衰减后信号的离散傅里叶变换与原始信号和阻尼因子在频域的循环卷积相同，所以衰减后信号的零频、低频分量是原始信号频率分量的加权和。Shin等[20]曾描述衰减波场中包含 “mirage-like”(海市蜃楼式)低频成分；Ha等[23]详细证明了这种衰减波场中确实存在零频和低频信息，并证明了这些信息的有效性，指出这些信息可以看作是原始信号的模糊处理结果。

为展示阻尼因子对低频成分的作用和效果，分别用地震子波及其振幅谱、地震记录及其相位谱加以说明。图3对比了子波经不同阻尼因子衰减后的结果。原始子波(图3a)是雷克子波经过切比雪夫滤波，去除了0～3Hz低频信息。阻尼因子为2、8和20的衰减子波如图3b～图3d所示。原始子波的最高波峰在中间，两侧分别有两个较小的波峰(图3a)；当阻尼因子为20时，子波的第一个波峰已经超过了中间的波峰值，第三个波峰几乎消失(图3d)。子波衰减是随着时间变化的，到达时间越早振幅增益作用越强，到达时间越晚振幅抑制作用越强，而且阻尼因子越大，作用强度越大。

图3 不同阻尼因子衰减后的子波归一化对比(a)原始； (b)σ=2； (c)σ=8； (d)σ=20

图4为图3四个子波的振幅谱，可以看出，随着阻尼因子的增大，子波中不存在的0～3Hz低频信息慢慢出现，当阻尼因子为20时，频谱中0～3Hz频率成分已十分明显。所以，对信号进行阻尼衰减可以重构出低频信息。

图5a为单炮地震记录，图5b为滤除3Hz以下频率成分后的记录，图6为两个记录的相位谱。相比于原始记录，滤波后记录的相位谱(图6b)中，0～3Hz的部分已经扭曲畸变。然后将两个记录分别进行阻尼因子为2、8和20的衰减，其相位谱如图7所示。滤波记录经过衰减后的相位谱(图7右)，0～3Hz部分不再扭曲，与原始记录衰减后的相位谱(图7左)明显相似。可以认为，经过衰减后的低频成分包含有地下介质的有效信息，可以使用该低频成分进行反演。

图4 不同阻尼因子衰减后的子波归一化振幅谱对比(a)原始； (b)σ=2； (c)σ=8； (d)σ=20

图5 原始地震记录(a)及其滤除0～3Hz低频成分后的记录(b)

图6 原始地震记录(a)及其滤波后记录(b)的归一化相位谱

图7 原始地震记录(左)及其滤波后记录(右)经不同阻尼因子衰减后的归一化相位谱(a)σ=2； (b)σ=8； (c)σ=20

3 基于L2范数的时域加权目标函数

当地震记录中引入时间阻尼因子后，远炮检距地震道的能量会被衰减得很小。Shin等[20]为了消除阻尼因子产生的这一负面影响，采用了对数目标函数，即

(15)

虽然对数目标函数可以消除指数衰减的负面影响，但对数目标函数本身存在固有的缺陷，即对小功率谱波场以及波场的振幅变化十分敏感，反演稳定性较差。为此，Jun等[39]提出了加权的L2目标函数，史才旺[27]提出了基于L1范数的对数目标函数。基于前人的研究，本文使用了基于时间域加权的L2目标函数。

本文的目标函数首先还是以传统的L2差值目标函数为核心，但是与前人不同的是，本方法在时间域进行加权处理，因为相较于频率域来说，时间域对于记录的操作更加方便和直观。时间域目标函数为

(16)

(17)

不论是用对数目标函数还是加权系数目标函数，其目的都是消除阻尼因子的负面影响，增强远炮检距地震道的能量，使其在反演中发挥更大的作用。本文应用地震道差值模的倒数作为加权系数，即

(18)

4 模型测试

使用Marmousi 2部分模型测试本文的反演方法。为节约反演运算时间，将Marmousi 2纵波速度模型进行了截取和抽稀，得到了网格节点数为401×160、网格尺寸为15m×15m的模型作为真实纵波速度模型(图8a)；纵横波速度比设为1.7，获得其横波速度模型(图8b)；密度设为常数2000kg/m3。为提高反演效率，默认陆地浅部水平层的速度已知(vP=1500m/s)。

正演采用主频为6Hz的雷克子波作为激发震源，记录长度为6s，采样间隔为1.5ms。炮点和检波点均匀布置在地面，炮点范围为0～6km，共41炮；检波点范围为0.15～5.85km，每炮381道。图9为模拟的第20炮两分量记录。

图8 Marmousi 2部分模型(a)纵波速度； (b)横波速度

使用了如图10所示的一维线性梯度模型作为初始模型，模型中不包含任何构造信息。采用多尺度反演策略，参照激发子波的主频，选取了4个频率组(表1)，前两组的频率低于3Hz，后两组的频率高于3Hz，每组分别迭代40次。

表1 反演频率分组参数

图9 第20炮观测记录(a)水平分量； (b)垂直分量

全频带数据常规时间—频率域反演最终结果如图11所示。与图8对比可以看出，反演结果十分接近真实模型，两侧层状介质、中部的异常体、深部的断裂和背斜构造都得到了精确重构。为了更为直观、定量地反映反演效果，抽取x=3.9km处(中部隆起处)的纵、横波速度曲线(图12)，可以发现，反演结果十分接近真实值。

滤除模拟记录3Hz以下的成分，以图10的一维线性梯度模型作为初始模型进行常规时域—频率域反演。因记录中的可用频率成分全部在3Hz以上，所以选取表1中的后两组频率进行反演。

图10 初始速度模型(a)纵波速度； (b)横波速度

图11 全频带常规时间—频率域反演结果(a)纵波速度； (b)横波速度

缺失低频成分的常规时间—频率域反演结果如图13所示。因没有精确的低波数更新，所以很难避免周波跳跃现象，以至于陷入局部极值。与真实模型(图8)对比可见，模型的大尺度构造几乎没有得到重构，只得到了一些小尺度的细节更新，并且产生了一些构造假象；断层、断块和模型深部的不整合面都不可识别。对比图13与图11可见，低频信息对常规反演的重要性。在初始模型精度不高的情况下，就需要使用记录中的低频信息重构出模型的大尺度构造，进而得到最终的高精度反演结果。若记录中缺失低频信息，反演就会极易陷入局部极小，导致反演失败。

应用图10的初始模型对缺失0～3Hz低频的模拟记录进行时间—拉普拉斯—频率域全波形反演。因低频已经通过拉普拉斯变换重构，所以可以使用0～3Hz的频率成分，经测试后选择了5个尺度，每个尺度进行10次迭代，具体参数如表2所示。

图14为时间—拉普拉斯—频率域反演的结果，可以看到，模型中大尺度的构造已经得到恢复，其形状与真实模型的轮廓十分相近。然后将该反演结果作为初始模型，进行常规时间—频率域反演，反演尺度参数与表1中的后两组相同。由图15的反演结果可以看出，与真实模型十分接近，各个层位都得到了精确的反演，断层、断块等构造也十分清晰。与图11相比，两种结果相差无几，说明该方法可以在缺失低频信息的情况下，达到含低频信息时的反演精度。

图12 x=3.9km处全频带常规时间—频率域反演速度曲线与真实及初始速度曲线的对比

为了测试本文方法对密度反演的适用性，将常密度修改为空间变化(图16)，与纵、横波速度同时参与反演。其他正、反演参数与以上相同。

图13 缺失低频成分时的常规时间—频率域反演结果(a)纵波速度； (b)横波速度

图14 缺失低频成分时的时间—拉普拉斯—频率域反演结果(a)纵波速度； (b)横波速度

对缺失0～3Hz低频成分的记录进行时间—拉普拉斯—频率域全波形反演，结果如图17所示，可以看出，纵、横波速度和密度模型的大尺度构造同样

表2 复频率域反演频率分组及参数

图15 以时间—拉普拉斯—频率域反演结果为初始模型的常规时间—频率域反演结果

图16 密度模型

得到了比较好的恢复。将该反演结果作为初始模型，进行时间—频率域反演，结果如图18所示。

从图18可以看出，纵、横波速度的反演结果同样十分接近真实模型。根据Ben-Hadj-Ali等[40]的模型误差公式，计算图15的纵、横波速度相对误差分别为0.2437、0.2355，图18的纵、横波速度相对误差分别为0.2497、0.2380，说明密度是否参与反演对纵、横波速度的反演结果几无影响，同时也说明本文方法同样适用于密度反演。对于密度项，由于与其他参数的耦合作用，其最终反演结果比速度反演略差，但是可以看出本文方法的密度反演结果仍有较高的可信度，可以体现出真实密度模型的构造形态与数值大小。上述测试结果表明，在缺失低频信息的情况下，本文方法仍然可以重构宏观构造模型，然后以此反演结果作为初始模型，进行常规时间—频率域反演，就可避免陷入局部极小，得到精确的反演结果。

图17 缺失低频成分时的时间—拉普拉斯—频率域纵波速度(a)、横波速度(b)和密度(c)反演结果

图18 以时间—拉普拉斯—频率域反演结果为初始模型的纵波速度(a)、横波速度(b)和密度(c)常规时间—频率域反演结果

5 结束语

本文建立了一种基于时域加权的拉普拉斯—频率域的弹性波全波形反演方法。该方法在结合时域和频率域各自优势的基础上，通过引入拉普拉斯阻尼因子，降低了全波形反演对于低频信息的依赖性。改进的目标函数，通过时域加权的形式消除了拉普拉斯阻尼因子对地震记录能量的衰减影响，构建形式方便、灵活。另外在变密度情况下，该方法仍能得到较高精度的纵、横波速度反演结果，密度反演结果由于受到其他参数的串扰，使得其精度略低于速度反演结果，因此如何改善该方法的密度反演精度应是继续研究的方向。