0 引言

全波形反演(Full Waveform Inversion，FWI)利用波动方程模拟得到的合成地震记录与实际地震记录的残差构建目标函数求取地下介质的参数，包括纵横波速度、各向异性参数和密度等。相较于其他简化波动方程的地震反演技术，如基于射线追踪的旅行时层析反演、基于褶积正演模型的地震反演等，FWI运用波动方程能更精确地模拟波场的传播过程，特别是地下构造复杂时，FWI可以获得更高精度和分辨率的反演结果。因此，近些年FWI成为了热点研究领域^[1-3]。

为了保证复杂构造模型的正演模拟精度，通常需要使用细小空间网格对低速层采样。然而，对整个介质进行细网格剖分会导致巨大的计算量以及存储需求，同时在模型的高速区也会产生过采样现象。为解决该问题，最优的方法是采用不同的空间网格步长对模型的不同区域进行剖分，在模型的浅层低速区使用细网格剖分，在中深层高速区使用粗网格剖分。基于变网格的波场模拟技术通常通过对粗、细网格过渡带的波场插值技术和基于变差分系数的变网格波场模拟技术实现^[4-5]。

粗细网格过渡带的波场插值技术在井间管波数值模拟^[6-7]、变网格逆时偏移成像^[8]、黏弹介质波场正演模拟^[9]等方面都有应用。该方法使用规则的矩形网格对模型中的高、低速体进行剖分，粗、细网格步长比只能是整数倍，在粗、细网格过渡带，需要对细网格波场下采样赋值给粗网格，并且需要对粗网格波场值进行插值以确定细网格各点波场值，通常采用的插值方法有线性插值、三角函数多项式插值等^[10-11]。目前，基于插值法的变网格技术依然面临着两大难题：一是地震波传播到粗、细网格分界面处时存在波场虚反射问题，二是波场在传播时间较长时会出现不稳定现象^[12-13]。Adriano等^[14]分析并总结了在变网格波动方程数值模拟中粗、细网格分界面存在虚反射的问题。Kang等^[15-16]提出时间、空间双变网格弹性波波动方程正演模拟方法，通过合成数据测试验证了该方法的可靠性和可行性。在此基础上，许多学者进行了进一步的研究。虽然在抑制虚反射上有一定效果，但是当介质复杂时，依旧会存在虚反射现象^[17]。针对粗、细网格过渡带插值不稳定的现象，许多学者改进了插值方法，如采用时空双变插值^[15-16]、波数域插值^[9]、转置矩阵法^[18]，或采用Lanczos滤波^[19]、Gaussian滤波^[20]、时间域滤波^[21]等方法，提高了波场模拟的稳定性。但对于介质变化剧烈、长时间模拟，依旧存在不稳定问题。

基于变差分系数的变网格技术不需要对波场进行插值，因此可以有效避免粗、细网格分界面造成的虚反射问题^[22]。鉴于变差分系数波场模拟方法的优势，许多学者进行了深入的研究和应用，例如基于变差分系数的起伏地表弹性波波动方程波场模拟^[23]、声学孔隙模型波场模拟^[24]、倾斜界面声波模拟^[25]等。

基于变差分系数的变网格波场模拟技术对网格进行非规则四边形的划分，通过直接求取各网格节点处的差分系数实现变网格波场模拟。目前，计算差分系数的方法主要包括基于Taylor展开式和基于最小二乘法求解两种方法。本文采用基于Taylor展开式的方法求取差分系数并进行变网格波场模拟。在此基础上，本文基于Zhang等^[26]提出的变网格方法，提出了一种基于变差分系数的变网格弹性波全波形反演方法。该方法在全波形反演的正演模拟和残差反向延拓过程中，针对不同尺度和速度的目标体，分别采用不同尺寸的网格进行剖分。对近地表低速带或小尺度目标体，使用细网格剖分，以保证波场模拟和反演精度；对中深层大尺度高速目标体，使用粗网格剖分，从而提高计算效率，并有效抑制正演模拟中出现的空间频散现象。

1 方法原理 1.1 弹性波动方程有限差分格式

有限差分正演模拟需要对地下整体空间进行网格离散化。规则网格有限差分方法将弹性介质参数、振动速度和应力设置在相同的网格点上进行正演模拟；而交错网格有限差分方法则是将它们分别放置于整网格节点、半网格点以及面元中心点进行正演模拟。交错网格有限差分方法相较于规则网格有限差分方法数值频散更小，模拟精度更高。三维弹性波动方程交错网格有限差分不同量的网格节点离散化分布如图 1所示。

使用交错网格有限差分模拟一阶速度—应力弹性波方程时，速度对时间的偏导可以通过应力在空间各个方向的偏导获得，而应力的时间偏导则可以通过速度的空间偏导求得。计算当前时间节点上的应力(速度)值，只需要前一个时间节点上的应力(速度)值以及两个时间节点之间的速度(应力)值，可以大量减少存储空间。当使用二阶时间精度差分时，一阶速度—应力弹性波方程的第一项可以写成

$ {v}_{x}\left(t+\frac{\mathrm{\Delta }t}{2}\right)={v}_{x}\left(t-\frac{\mathrm{\Delta }t}{2}\right)+\frac{\mathrm{\Delta }t}{\rho }\left(\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xy}}{\partial y}+\frac{\partial {\sigma }_{xz}}{\partial z}\right) $

(1)

式中：Δt为时间采样间隔。由Taylor级数展开得到应力空间导数的任意2M阶空间精度的差分格式($ \partial {\sigma }_{xx}/\partial x $为例)为

$ \frac{\partial {\sigma }_{xx}}{\partial x}=\sum\limits_{i=1}^{2M}{C}_{i}{\sigma }_{xx}\left[x+{\left(-1\right)}^{i+1}\mathrm{\Delta }x\left(i-\frac{1}{2}\right)\right] $

(2)

式中：Δx为空间采样间隔；$ {C}_{i} $为交错网格2M阶有限差分系数。

1.2 基于变差分系数的变网格正演

在基于变差分系数的变网格有限差分正演模拟方法中，差分系数的推导方式主要有两种：第一种基于Taylor展开式^[27-29]，第二种通过最小二乘拟合法求取最优差分系数。本文使用前一种方法。在任意方向的网格剖分如图 2所示，其中，蓝色方框表示整节点的应力值p，红色三角形表示半网格节点的质点振动速度v，i表示网格节点号，$ \mathrm{d}{x}_{i} $表示网格剖分的空间步长，$ {d}_{i}\left(i=\mathrm{1, 2}, \cdots \right) $表示空间差分步长增量，用以描述参与差分计算点(蓝色方块)和计算点(红色三角形)的间隔，可以由空间步长计算得到。

基于Taylor展开的任意波场分量对空间方向的导数为

$ \frac{\partial f}{\partial x}=\sum\limits_{i=1}^{2M}{C}_{i}f\left[x+{\left(-1\right)}^{i+1}{d}_{i}\right] $

(3)

将式(3)变换到波数域，有

$ ik=\sum\limits_{i=1}^{2M}{C}_{i}\mathrm{e}\mathrm{x}\mathrm{p}\left[{\left(-1\right)}^{i+1}ik{d}_{i}\right] $

(4)

式中k表示有效波数。

根据Pitarka^[29]的推导，有限差分系数能够展开到2M阶精度，本文基于空间四阶精度推导式(4)。对$ \mathrm{e}\mathrm{x}\mathrm{p}\left[{\left(-1\right)}^{i+1}ik{d}_{i}\right] $展开，有

$ \begin{array}{l}\mathrm{e}\mathrm{x}\mathrm{p}\left[{\left(-1\right)}^{i+1}ik{d}_{i}\right]=1+{\left(-1\right)}^{i+1}ik{d}_{i}+\frac{1}{2}{\left[{\left(-1\right)}^{i+1}ik\right]}^{2}{d}_{i}^{2}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \begin{array}{cc}& \end{array}\frac{1}{6}{\left[{\left(-1\right)}^{i+1}ik\right]}^{3}{d}_{i}^{3}+O\left({d}_{i}^{4}\right)\begin{array}{cc}& i=\mathrm{1, 2}, \mathrm{3, 4}\end{array}\end{array} $

(5)

将式(5)代入式(4)，推导出四阶差分系数满足方程组

$ \left(\begin{array}{cccc}1& 1& 1& 1\\ {d}_{1}& -{d}_{2}& {d}_{3}& -{d}_{4}\\ {d}_{1}^{2}& {d}_{2}^{2}& {d}_{3}^{2}& {d}_{4}^{2}\\ {d}_{1}^{3}& -{d}_{2}^{3}& {d}_{3}^{3}& -{d}_{4}^{3}\end{array}\right)\left(\begin{array}{c}{C}_{1}\\ {C}_{2}\\ {C}_{3}\\ {C}_{4}\end{array}\right)=\left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right) $

(6)

对于任意2M阶有限差分系数矩阵可以统一表示为

$ \boldsymbol{A}\boldsymbol{C}=\boldsymbol{b} $

(7)

式中：$ \boldsymbol{C}={\left({C}_{1}, {C}_{2}, \cdots , {C}_{2M}\right)}^{\mathrm{T}} $；$ \boldsymbol{b}={\left(\mathrm{0, 1}, 0, \cdots , 0\right)}^{\mathrm{T}} $。变换矩阵A由空间差分增量$ {d}_{i} $计算

$ {A}_{n, i}={\left[{\left(-1\right)}^{i+1}\right]}^{n-1}{d}_{i}^{n-1} $

(8)

图 2所示网格剖分的空间差分步长增量$ {d}_{i} $可以由空间步长$ \mathrm{d}{x}_{i} $计算得到，二者之间的关系为

$ \left\{\begin{array}{l}{\boldsymbol{d}}^{\mathrm{h}}={\boldsymbol{A}}^{\mathrm{h}}{\boldsymbol{d}}_{x}^{\mathrm{h}}\\ {\boldsymbol{d}}^{\mathrm{f}}={\boldsymbol{A}}^{\mathrm{f}}{\boldsymbol{d}}_{x}^{\mathrm{f}}\end{array}\right. $

(9)

式中：$ {\boldsymbol{d}}^{\mathrm{h}}={\left({d}_{1}^{\mathrm{h}}, {d}_{2}^{\mathrm{h}}, \cdots , {d}_{2M}^{\mathrm{h}}\right)}^{\mathrm{T}} $和$ {\boldsymbol{d}}^{\mathrm{f}}={\left({d}_{1}^{\mathrm{f}}, {d}_{2}^{\mathrm{f}}, \cdots , {d}_{2M}^{\mathrm{f}}\right)}^{\mathrm{T}} $分别表示整网格节点和半网格节点的空间差分步长增量；$ {\boldsymbol{d}}_{x}^{\mathrm{h}}={\left({d}_{x-M+1}^{\mathrm{h}}, \cdots , {d}_{x}^{\mathrm{h}}, \cdots {d}_{x+M-1}^{\mathrm{h}}\right)}^{\mathrm{T}} $和$ {\boldsymbol{d}}_{x}^{\mathrm{f}}={\left({d}_{x-M}^{\mathrm{f}}, \cdots , {d}_{x}^{\mathrm{f}}, \cdots , {d}_{x+M-1}^{\mathrm{f}}\right)}^{\mathrm{T}} $分别表示整网格节点和半网格节点的空间步长。当空间精度为四阶时，整网格节点的变换系数矩阵A为

$ {\boldsymbol{A}}_{4}^{\mathrm{h}}=\left(\begin{array}{ccc}0& \frac{1}{2}& 1\\ 1& \frac{1}{2}& 0\\ 0& \frac{1}{2}& 0\\ 0& \frac{1}{2}& 0\end{array}\right) $

(10)

空间高阶精度变换矩阵可以通过空间四阶精度变换矩阵获得

$ {\boldsymbol{A}}_{2M}^{\mathrm{h}}=\left[{\boldsymbol{\alpha }}_{1}, \left(\begin{array}{c}{\boldsymbol{B}}^{\mathrm{h}}\\ {\boldsymbol{A}}_{2\left(M-1\right)}^{\mathrm{h}}\end{array}\right), {\boldsymbol{\alpha }}_{2}\right] $

(11)

$ \left\{\begin{array}{l}{\boldsymbol{\alpha }}_{1}={\left(\mathrm{0, 1}, 0, \cdots , 0\right)}^{\mathrm{T}}\\ {\boldsymbol{B}}^{\mathrm{h}}=\left(\begin{array}{c}0, \cdots , 0, \frac{1}{2}, 1, \cdots , 1\\ 1, \cdots , 1, \frac{1}{2}, 0, \cdots , 0\end{array}\right)\\ {\boldsymbol{\alpha }}_{2}={\left(\mathrm{1, 0}, \cdots , 0\right)}^{\mathrm{T}}\end{array}\right. $

(12)

式中：B^h为矩阵，其维度与空间差分阶数有关，为2×(2M-3)；α₁和α₂均为2M×1的一维向量。空间四阶精度半网格节点上变换系数矩阵可以表示为

$ {\boldsymbol{A}}_{4}^{\mathrm{f}}=\left(\begin{array}{cccc}0& 0& 1& \frac{1}{2}\\ \frac{1}{2}& 1& 0& 0\\ 0& 0& \frac{1}{2}& 0\\ 0& \frac{1}{2}& 0& 0\end{array}\right) $

(13)

任意空间高阶半网格节点变换矩阵，可由空间四阶精度获得

$ {\boldsymbol{A}}_{2M}^{\mathrm{f}}=\left[{\boldsymbol{\beta }}_{1}, \left(\begin{array}{c}{\boldsymbol{B}}^{\mathrm{f}}\\ {\boldsymbol{A}}_{2\left(M-1\right)}^{\mathrm{f}}\end{array}\right), {\boldsymbol{\beta }}_{2}\right] $

(14)

$ \left\{\begin{array}{l}{\boldsymbol{\beta }}_{1}={\left(0, \frac{1}{2}, 0, \cdots , 0\right)}^{\mathrm{T}}\\ {\boldsymbol{B}}^{\mathrm{f}}=\left(\begin{array}{c}0, \cdots , \mathrm{0, 1}, \cdots , 1\\ 1, \cdots , \mathrm{1, 0}, \cdots , 0\end{array}\right)\\ {\boldsymbol{\beta }}_{2}={\left(\frac{1}{2}, \mathrm{0, 0}, \cdots , 0\right)}^{\mathrm{T}}\end{array}\right. $

(15)

式中：B^f的维度为2×(2M-2); β₁和β₂的维度均为2M×1。在获得变换矩阵和差分空间步长增量后，可通过式(6)求取差分系数。

1.3 基于变差分系数的弹性波全波形反演

全波形反演的目标函数可以表示为模拟记录和观测记录残差的L₂范数

$ \chi \left(\boldsymbol{m}\right)=\frac{1}{2}{||{\boldsymbol{D}}^{\mathrm{s}\mathrm{y}\mathrm{n}}-{\boldsymbol{D}}^{\mathrm{o}\mathrm{b}\mathrm{s}}||}_{2}^{2} $

(16)

式中：$ \boldsymbol{D}^{\rm{syn}} $为合成数据；$ \boldsymbol{D}^{\rm{obs}} $为观测数据。在反演过程中通过迭代来更新速度模型，使合成数据不断接近观测数据。模型更新公式为

$ {\boldsymbol{m}}_{k+1}={\boldsymbol{m}}_{k}+{\gamma }_{k}{\rm{ \mathsf{ δ} }}{\boldsymbol{m}}_{k} $

(17)

式中：m_k表示第k次迭代的反演结果；$ {\gamma }_{k} $表示第k次迭代的梯度更新步长；$ {\rm{ \mathsf{ δ} }}{\boldsymbol{m}}_{k} $表示第k次迭代梯度搜索方向。本文采用共轭梯度法求解$ {\rm{ \mathsf{ δ} }}{\boldsymbol{m}}_{k} $。

2 数值模拟和反演算例 2.1 基于变差分系数的变网格有限差分数值模拟

本文首先采用一维变网格波场模拟验证方法的正确性。其中，纵波速度为4000 m/s，震源为30 Hz主频的Ricker子波。当网格变化比r=1.0时，所有网格步长均为4 m；当r=2.0时，表示网格步长在0~400 m深度范围内为2 m，在400~1600 m内为4 m, 在1600~3500 m内为8 m；当r=2.5时，网格步长在0~320 m内为1.6 m, 在320~1520 m内为4 m, 在1520~3500 m内为10 m。图 3为三种网格模拟的200、400 ms时的波形，可以看出，三种网格模拟结果非常一致，说明了变网格模拟方法的正确性。

本文对变网格剖分下波场对空间方向求导的精度也进行了测试。当网格步长变化比为1.0和5.0时，正弦函数离散化差分求导与真实导数曲线几乎没有差异(图 4a、图 4b)，当网格步长变化比为10.0时，二者的差异也非常小(图 4c)。

对二维常速度模型(纵波速度为4000 m/s，横波速度为2353 m/s)在浅层进行变网格剖分(图 5a)，其中横向网格步长固定为10 m，深度0~400 m内纵向网格步长为5 m，400~2000 m为10 m。变网格模拟的v_z分量0.24 s的波场快照如图 5b所示，与细网格(两个方向网格步长均为5 m)模拟结果的差如图 7a所示。常速度模型中层变网格剖分(图 6a，在深度500~900 m纵向网格步长为5 m，其他为10 m)，其变网格模拟的v_z分量0.24 s的波场快照如图 6b所示，与细网格(两个方向网格步长均为5 m)模拟结果的差如图 7b所示。由图 5~图 7可见，两种变网格剖分模拟结果非常接近，与细网格的残差非常小，进一步验证了方法的正确性和灵活性。三种网格剖分正演的计算耗时和内存占用如表 1所示，可见，变网格正演模拟可以在保证模拟精度的同时，提高计算效率、减小内存的占用。

2.2 基于变差分系数的变网格全波形反演

将基于变差分系数正演方法应用于多尺度弹性波全波形反演。图 8为真实的二维推覆体速度模型，用120 m×120 m的平滑半径其平滑结果如图 9所示。首先对初始模型进行8 m×8 m的均匀网格剖分，分别采用主频为2.5、3.5、5.5、9.0、13.0 Hz和20 Hz的雷克子波进行多尺度弹性波全波形反演，从低频到高频逐步更新模型，将低频反演结果作为高频反演的初始模型，主频20 Hz的子波反演结果如图 10~图 11所示。然后采用空间变网格进行剖分(横向网格步长固定为8 m，深度0~320 m范围内纵向网格步长为4 m，其余为8 m)，其反演结果如图 12~图 13所示。图 14~图 15为均匀网格反演结果和变网格反演结果的浅层放大对比，可以看出，变网格反演结果中的浅层更准确，并且层界面也比均匀网格反演结果更清晰。

为了测试本文方法的抗噪性，分别对合成炮集数据添加不同信噪比的噪声，纵波速度的反演相对误差曲线如图 16所示。与无噪声的全波形反演结果相比，噪声会降低反演结果精度，但由于炮集数据的多次覆盖，反演结果精度降低有限。

3 结论

为了提高近地表低速异常体的反演精度，本文在多尺度全波形反演方法的基础上进行了基于变差分系数的变网格正演和反演研究。基于变差分系数的变网格正演模拟可以在一定条件下实现任意网格变化比的波场模拟。该方法几乎不存在虚反射问题，并且波场能够稳定传播。变差分系数的变网格正演模拟方法在达到细网格模拟精度的基础上，可以有效减少计算量和内存占用。将基于变差分系数的变网格正演模拟应用于时间域弹性波多尺度全波形反演，针对浅地表低速体进行细网格剖分，可以在保证计算效率的同时提高浅层空间采样率，更精确地刻画近地表的速度变化。

本文采用的变差分系数变网格模拟方法局限性在于无法根据速度模型自适应剖分网格。