0 引言

地下介质的速度建模是地震资料成像的重要环节。全波形反演(FWI)充分利用了地震波的运动学和动力学信息，是目前速度建模精度和分辨率最高的方法。Tarantola^[1]最早提出了时间域FWI，它通过最小化计算地震数据和观测地震数据之间的残差迭代更新速度模型^[2-5]。

拖缆资料的FWI是海洋地震勘探具有潜力的速度建模方式之一。由于受仪器设备的限制，采集到的拖缆资料通常缺少可靠的低频成分。而FWI是复杂的非线性问题，其效果严重依赖于有效的低频地震数据。这导致拖缆资料的FWI发生周期跳跃现象，目标函数收敛到局部极小值，无法获得准确的速度模型。为此，学者们提出了许多改进的波形反演方法，这些方法可以概括为两类。第一类方法为构建非线性程度较低的目标函数。Shin等^[6]提出Laplace域FWI方法，将地震数据变换到Laplace域，并与观测的地震数据进行拟合计算，构建FWI的目标函数^[7-8]。地震数据在Laplace域具有更多低频成分，目标函数非线性程度更低，但此类方法严重依赖于Laplace衰减常数的选取和大炮检距地震数据。Choi等^[9]利用观测地震数据与计算地震数据的全局互相关构建FWI的目标函数，在反演中使观测地震数据和计算地震数据的相位相互匹配。与常规目标函数相比，相位匹配的非线性程度较低^[10]。Bozdağ等^[11]使用地震数据的包络建立FWI的目标函数，即使地震数据缺少低频信息，地震数据的包络仍然具有丰富的低频信息^[12-14]。

第二类方法是重构观测的低频地震数据。Warner等^[15]提出通过非线性插值方法从高频地震数据中合成低频数据，并应用于三维海底电缆资料的FWI，获得了较准确的低波数模型。但当噪声能量增大时，由插值生成的低频数据质量急剧下降。Zhang等^[16]通过提取被动源地震数据中的低频信息重构主动源地震数据的低频部分，但是被动源地震受噪声严重干扰，重构效果无法保证^[17]。陈生昌等^[18]提出对散射波场进行时间二阶积分运算重构低频地震数据，由于对波场的二阶时间积分相当于低通滤波，此方法仅适用于散射波场具有一定的低频信息的情况。Zhang等^[19]提出通过盲反褶积重构低频地震数据的方法，但是当初始地震子波不够准确时，盲反褶积计算出的反射系数序列存在较大误差，重构的地震数据不够可靠。上述方法都是通过数学运算方式恢复的“数值低频”，与物理意义上的低频数据仍有一定差异。

与人工重构得到的低频数据相比，实际采集到的低频数据更具可靠性。OBN（海底节点）/OBS（海底地震仪）资料的另一个优势是宽方位角采集，能够提供大炮检距地震数据，对高陡构造以及地下介质深部的反演具有关键作用^[20-21]。

另外，FWI的效果强烈依赖于准确的震源子波。当观测地震数据的震源子波与反演使用的震源子波不匹配时，反演效果急剧变差。然而，实际采集到的地震资料中，震源子波是未知的，在拖缆与OBN资料的联合反演中，由于两种地震资料的震源子波存在差异，因此反演结果受震源子波的影响更加严重。因为OBS资料中含有丰富的低频成分，Yang等^[22]使用OBS资料与拖缆资料进行联合波形反演，但因未考虑子波差异的影响导致反演结果不佳。为此，学者们开发了不依赖震源子波的FWI方法。Lee等^[23]提出了使用单一参考道对频率域地震记录进行归一化处理；Zhou等^[24]提出了使用单炮平均振幅进行归一化处理。在频率域进行归一化处理，能够消除目标函数中的震源子波项，达到不依赖震源子波的目的。Choi等^[25]提出在频率域将模拟数据与观测数据分别与对方的参考道相乘，使目标函数中具有相同的震源子波项，消除震源子波不匹配的影响。后来，Choi等 ^[26]将参考道方法从频率域推广到时间域，模拟数据与观测数据分别与对方的参考道褶积。此外，辛天亮等^[27]提出了基于数据相似性的频率域全波形反演方法，也达到了不依赖震源子波的效果。

本文提出一种不依赖震源子波的海洋拖缆与OBN资料联合全波形反演方法。通过含有丰富低频数据的OBN资料与拖缆资料联合反演，消除因拖缆资料低频缺失造成的周期跳跃。利用时间域的平均道方法进行联合反演不依赖于震源子波，并且与传统的参考道方法相比，平均道方法具有更好的抗噪性。用Marmousi模型合成的拖缆资料与OBN资料完成模型测试，结果验证了本文方法的可行性和有效性。

1 方法原理 1.1 目标函数

传统FWI的目标函数定义为观测地震数据与计算地震数据残差的二范数^[28]。由于实际采集的拖缆资料常常缺少可靠的低频数据，单独进行FWI时容易发生周期跳跃，因此需要与OBN资料联合反演，此时目标函数为

$ \begin{array}{l}E\left(v\right)=\frac{{\alpha }_{1}}{2}\sum\limits_{s}\sum\limits_{g}{‖{\boldsymbol{d}}_{\mathrm{c}\mathrm{a}\mathrm{l}}^{\mathrm{T}\mathrm{S}}({\boldsymbol{x}}_{g}, {\boldsymbol{x}}_{s})-{\boldsymbol{d}}_{\mathrm{o}\mathrm{b}\mathrm{s}}^{\mathrm{T}\mathrm{S}}({\boldsymbol{x}}_{g}, {\boldsymbol{x}}_{s})‖}_{2}^{2}+\\ \frac{{\alpha }_{2}}{2}\sum\limits_{s}\sum\limits_{o}{‖{\boldsymbol{d}}_{\mathrm{c}\mathrm{a}\mathrm{l}}^{\mathrm{O}\mathrm{B}\mathrm{N}}({\boldsymbol{x}}_{o}, {\boldsymbol{x}}_{s})-{\boldsymbol{d}}_{\mathrm{o}\mathrm{b}\mathrm{s}}^{\mathrm{O}\mathrm{B}\mathrm{N}}({\boldsymbol{x}}_{o}, {\boldsymbol{x}}_{s})‖}_{2}^{2}\end{array} $

(1)

式中：v为速度；s为震源序号；g为拖缆的检波器序号；x_g为检波器的空间坐标；o为OBN序号；x_o为OBN的空间坐标；d_cal与d_obs分别是计算地震数据与观测地震数据；上标TS和OBN分别表示拖缆资料与OBN资料；x_s为震源的空间位置；α₁和α₂分别为拖缆资料和OBN资料在联合反演中的权重。

由于采用了从低频到高频的多尺度反演策略，为了让OBN资料在反演前期发挥作用，弥补拖缆资料的低频缺失，权重设置如下

$ {\alpha }_{1}=\frac{k}{n} $

(2)

$ {\alpha }_{2}=\frac{n-k}{n} $

(3)

式中：k为当前迭代次数；n为反演的总迭代次数。这样设置能够使OBN资料的权重α₂在反演速度模型的低波数分量时取较大的值，对速度模型的更新影响较大。

上述联合反演是建立在震源子波已知的情况下，而实际勘探中震源子波是未知的。传统不依赖震源子波的反演方法基于参考道实现，将观测数据与计算数据分别与对方的参考道进行褶积。但是参考道方法严重依赖于参考道的选取，并且与参考道褶积的过程会导致反演抗噪性降低。因此本文用平均道代替参考道，一方面避免了参考道的选取问题，另一方面求取平均道的过程一定程度上压制了地震数据中的随机噪声。

基于平均道的全波形反演目标函数定义为

$ E=\sum\limits_{i=1}^{{N}_{\mathrm{s}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}{‖{\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}‖}_{2}^{2} $

(4)

式中：$ \mathrm{*} $表示褶积；N_s为炮点数；N_r为检波器数；A^obs表示观测数据的平均道；A^cal表示计算数据的平均道。平均道定义为

$ {\boldsymbol{A}}_{i}=\frac{1}{{N}_{\mathrm{r}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}{\boldsymbol{d}}_{i, j} $

(5)

平均道方法不依赖震源子波的原理与单一参考道类似。将式(4)中的地震数据写成震源子波与格林函数褶积的形式，即

$ E=\sum\limits_{i=1}^{{N}_{\mathrm{s}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}‖{\boldsymbol{s}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{G}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}\left(\frac{1}{{N}_{\mathrm{r}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}{\boldsymbol{s}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{G}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right)- \\{{\boldsymbol{s}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{G}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}\left(\frac{1}{{N}_{\mathrm{r}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}{\boldsymbol{s}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{G}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\right)‖}_{2}^{2} $

(6)

式中：s_i为震源子波；$ \boldsymbol{G} $_{i, j}为格林函数。

由于褶积运算满足分配律和交换律，式(6)可变形为

$ E=\sum\limits_{i=1}^{{N}_{\mathrm{s}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}‖{\boldsymbol{s}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{s}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{G}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}\left(\frac{1}{{N}_{\mathrm{r}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}{\boldsymbol{G}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right)- \\ {{\boldsymbol{s}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{s}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{G}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}\left(\frac{1}{{N}_{\mathrm{r}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}{\boldsymbol{G}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\right)‖}_{2}^{2} $

(7)

式(7)中做差的两项具有相同的新震源子波s_i^obs*s_i^cal，因此消除了反演中震源子波的不匹配影响，达到不依赖震源子波的效果。由此可以建立基于平均道的不依赖震源子波的拖缆和OBN资料联合FWI目标函数

$ \begin{aligned} E(v)= & \frac{\alpha_1}{2} \sum_s \sum_g \| \boldsymbol{d}_{\mathrm{cal}}^{\mathrm{TS}}\left(\boldsymbol{x}_g, \boldsymbol{x}_s\right) * \boldsymbol{A}_{\mathrm{obs}}^{\mathrm{TS}}\left(\boldsymbol{x}_s\right)- \\ & \boldsymbol{d}_{\mathrm{obs}}^{\mathrm{TS}}\left(\boldsymbol{x}_g, \boldsymbol{x}_s\right) * \boldsymbol{A}_{\mathrm{cal}}^{\mathrm{TS}}\left(\boldsymbol{x}_s\right) \|_2^2+ \\ & \frac{\alpha_2}{2} \sum_s \sum_o \| \boldsymbol{d}_{\mathrm{cal}}^{\mathrm{OBN}}\left(\boldsymbol{x}_o, \boldsymbol{x}_s\right) * \boldsymbol{A}_{\mathrm{obs}}^{\mathrm{OBN}}\left(\boldsymbol{x}_s\right)- \\ & \boldsymbol{d}_{\mathrm{obs}}^{\mathrm{OBN}}\left(\boldsymbol{x}_o, \boldsymbol{x}_s\right) * \boldsymbol{A}_{\mathrm{cal}}^{\mathrm{OBN}}\left(\boldsymbol{x}_s\right) \|_2^2 \end{aligned} $

(8)

1.2 正演模拟

本文使用常密度声波方程进行正演模拟

$ \frac{1}{{v}^{2}\left(\boldsymbol{x}\right)}\frac{{\partial }^{2}p(\boldsymbol{x}, t\left|{\boldsymbol{x}}_{s}\right.)}{\partial {t}^{2}}-{\nabla }^{2}p(\boldsymbol{x}, t\left|{\boldsymbol{x}}_{s}\right.)={f}_{\mathrm{s}}(\boldsymbol{x}, t\left|{\boldsymbol{x}}_{s}\right.) $

(9)

式中：x为空间位置；t为时间；p为波场；f_s为地震子波。使用与褶积重构低频地震数据相同的地震子波。正演时使用规则网格有限差分正演，时间上为2阶精度，空间上为12阶精度。使用完全匹配层方法以消除边界反射。

1.3 局部优化算法

最小化目标函数时，使用共轭梯度法^[29]，即通过求取目标函数的共轭梯度计算目标函数的下降方向。目标函数的共轭梯度由下式求得

$ \boldsymbol{d}{}_{k}(\boldsymbol{x})=\left\{\begin{array}{l}-\nabla {E}_{0}\left(\boldsymbol{x}\right) \;\;\;\;\;\;\;k=0\\ -\nabla {E}_{k}\left(\boldsymbol{x}\right)+{\beta }_{k}{\boldsymbol{d}}_{k-1}\left(\boldsymbol{x}\right)k > 0\end{array}\right. $

(10)

式中参数$ {\beta }_{k} $的作用是使$ {\boldsymbol{d}}_{k}\left(\boldsymbol{x}\right) $和$ {\boldsymbol{d}}_{k-1}\left(\boldsymbol{x}\right) $正交。$ {\beta }_{k} $由下式求得^[30]

$ {\beta }_{k}=\frac{〈\nabla {E}_{k}, \nabla {E}_{k}〉}{〈{\boldsymbol{d}}_{k-1}, \nabla {E}_{k}-\nabla {E}_{k-1}〉} $

(11)

1.4 梯度计算

联合反演目标函数的梯度为拖缆与OBN资料各自的梯度加权后求和

$ \nabla E\left(\boldsymbol{x}\right)={\alpha }_{1}{\boldsymbol{G}}^{\mathrm{T}\mathrm{S}}\left(\boldsymbol{x}\right)+{\alpha }_{2}{\boldsymbol{G}}^{\mathrm{O}\mathrm{B}\mathrm{N}}\left(\boldsymbol{x}\right) $

(12)

式中$ \boldsymbol{G} $^TS、$ \boldsymbol{G} $^OBN分别为由拖缆资料和OBN资料计算出的梯度。单一地震资料的梯度计算在频率域进行，这样可以更好地实现多尺度方法，降低反演的非线性。频率域中，通过将相同频率的正传波场与反传残差波场相乘计算目标函数的梯度，即

$ \nabla E\left(x\right)=\mathrm{R}\mathrm{e}\left[\sum\limits_{s}\sum\limits_{\omega }{\omega }^{2}{p}_{\mathrm{b}}^{\mathrm{*}}(\boldsymbol{x}, \omega \left|{\boldsymbol{x}}_{s}\right.){p}_{\mathrm{f}}(\boldsymbol{x}, \omega \left|{\boldsymbol{x}}_{s}\right.)\right] $

(13)

式中：Re表示取复数的实部；$ \omega $为角频率；p_b和p_f分别是正传波场和反传波场；上标$ \mathrm{*} $表示复共轭。频率域波场可以由时间域波场通过离散傅里叶变换(DFT)得到。

正传波场的震源在炮点位置，震源子波是重构低频地震数据使用的子波。反传波场的震源在检波点位置，反传震源由伴随状态法^[31]推导(推导过程见附录A)得到

$ {\boldsymbol{r}}_{i, j}^{\text{'}}={\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\otimes \left({\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\right) $

(14)

$ {\boldsymbol{r}}_{i, j}^{″}=-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\otimes \left({\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\right) $

(15)

式中：$ \otimes $表示互相关；$ {\boldsymbol{r}}_{i, j}^{\text{'}} $和$ {\boldsymbol{r}}_{i, j}^{″} $分别是反传震源的第一项和第二项。反演中，对拖缆与OBN资料分别进行反传波场计算，单独计算两种资料的梯度。

计算出共轭梯度之后，对速度模型$ {v}_{k}\left(\boldsymbol{x}\right) $进行更新

$ {v}_{k+1}\left(\boldsymbol{x}\right)={v}_{k}\left(\boldsymbol{x}\right)+{\alpha }_{k}{\boldsymbol{d}}_{k}\left(\boldsymbol{x}\right) $

(16)

式中α_k是速度模型第$ k $次迭代的步长，可以通过抛物线插值法或线性搜索法^[32-33]求得。

2 模型测试

为了验证本文方法的可行性与有效性，用Marmousi模型（图 1）合成数据进行测试。模型尺寸为5.0 km×1.5 km，模型的速度范围为1.5~4.5 km/s。使用的初始速度模型与真实速度模型相差较大，速度随深度呈线性变化，速度范围为1.50~4.03 km/s。在模型上表面共激发50炮，炮间距为100 m。拖缆资料由500个水听器以10 m间距在水面同时接收，OBN资料由10个OBN以500 m间距在海底同时接收。进行正演模拟时，合成拖缆资料的震源子波主频为15 Hz，合成OBN资料的震源子波主频为8 Hz，因为实际勘探中OBN资料主频比拖缆资料低，为了模拟低频缺失的拖缆资料，对合成的拖缆资料进行5 Hz的高通滤波。

图 2为模拟的拖缆资料以及OBN资料。从频谱可知，合成的OBN地震资料含有丰富的5 Hz以下的低频数据。

2.1 拖缆与OBN资料联合FWI

对滤除低频的拖缆资料进行FWI，反演使用的8个频率分别为2、3、4.5、6.3、8.82、12.348、17.2872和24.2021 Hz，每个频率迭代10次，得到图 3a所示结果。对不同间距的OBN数据进行FWI，反演使用的8个频率分别为1、1.4、1.96、2.744、3.8416、5.3782、7.5295和10.5414 Hz，每个频率迭代10次。OBN资料间距分别为500、1000、1500 m时，反演结果如图 3b~图 3d所示。由图可见，使用缺失低频部分的拖缆资料单独进行FWI，反演结果中出现了局部极小值现象，无法得到准确的反演结果。单独使用OBN资料进行反演时，由于其覆盖次数不足，无法反演速度模型高分辨率部分，并且随着OBN间距增大到1500 m，反演出现严重错误。不同间距的OBN资料与拖缆资料进行联合反演时使用的频率和迭代次数与单一资料反演相同。OBN资料间距分别为500、1000、1500 m时，联合反演结果如图 4所示。由图可见，联合反演结果的效果得到显著提高。但当OBN间距过大时，OBN资料对反演的改善作用不明显。

本文FWI采用基于MPI的炮并行运算策略，模型测试中，单一资料FWI约耗时2.1 h，联合FWI约耗时3.5 h。

图 5对比了理论模型、初始速度模型和几种反演方法的反演结果于水平方向2.5 km处速度随深度变化的曲线。从速度—深度曲线对比也能够看出：拖缆与OBN资料联合FWI结果比两种资料单独反演结果更接近真实模型。随着OBN间距增大，联合FWI效果逐渐变差。

另外，为了验证联合FWI中OBN资料大炮检距的作用，使用尺寸为15 km×4 km的Marmousi模型进行了测试。拖缆资料采用单边接收，缆线长度为6 km，OBN以500 m间距布设在模型的整个海底位置。拖缆资料道间距为12.5 km，炮点设置在水平方向6~15 km位置的海面，炮间距125 m，拖缆与OBN同时采集。图 6为理论模型与初始模型。分别使用拖缆资料、OBN资料以及拖缆与OBN联合三种方式进行反演，得到图 7所示结果。由图 7a的拖缆资料FWI结果可以看出，在炮检距受限制的情况下，模型深部以及水平方向7~12 km处的高陡构造反演效果较差。图 7b的OBN资料FWI结果中，由于OBN布设稀疏，导致反演模型的分辨率较低。当两种资料进行联合FWI，得到图 7c所示结果，较好地反演了模型深部和高陡构造，相比图 7b结果分辨率更高。

2.2 不依赖震源子波的联合FWI

上述测试是建立在合成资料与反演使用相同震源子波的基础上，为了测试基于平均道的不依赖震源子波的联合反演方法的可行性与有效性，在反演中使用与合成地震数据不同的子波。采用三组子波进行正演，分别合成拖缆与OBN数据。第一组为将不同振幅（图 8a上)和相位（图 8a下)的Ricker子波叠加得到的新子波；第二组为高斯函数的一阶导数（图 8b上)及其相反数（图 8b下)；第三组为复杂子波（图 8c上)及其相反数（图 8c下)。对每组子波合成的地震数据分别用式（1）目标函数的传统联合FWI方法和本文提出的不依赖震源子波联合FWI方法进行对比实验，反演中正演模拟拖缆数据和OBN数据的震源子波使用主频分别为15 Hz和8 Hz的Ricker子波。反演结果如图 9所示。

由图 9可知，在联合反演中引入平均道方法，即使反演使用的子波与地震数据子波不一致，仍然能够得到正确的反演结果，说明反演可以不依赖震源子波。基于三组子波合成地震数据的对比实验证明了本文所提方法能够解决常规联合FWI中子波不匹配的问题。

2.3 含噪数据测试

为了测试基于平均道的联合FWI方法在含噪数据中的优势，在合成的拖缆数据中加入不同信噪比(S/N)的高斯噪声。然后在合成地震数据子波与反演子波不一致时，分别使用单一参考道方法、平均道方法进行全波形反演，对两种反演方法的结果进行对比。单一参考道方法的反演实验中，选取最小炮检距的地震道作为参考道。图 10为加入不同比例噪声后的拖缆地震数据。图 11为基于单一参考道方法和基于平均道方法的反演结果。由图可知，当地震数据含噪时，使用单一参考道方法进行反演，受噪声影响严重。使用平均道方法时，由于求平均道的过程一定程度上压制了随机噪声，因此反演结果优于单一参考道方法。结果对比证明，在各噪声条件下，基于平均道的全波形反演方法对于随机噪声均具有更好的抗噪性。当噪声过强时，两种方法效果均受影响严重，平均道方法效果仍优于单一参考道方法。

3 结束语

全波形反演是一个强非线性问题，局部优化算法容易陷入局部极小值。当观测地震数据缺少低频，初始速度模型不够准确时，往往得到错误的反演结果。现有的拓频方法大多是通过数学变换获取地震数据的数值低频，缺少物理意义上的联系。使用本身含有丰富低频数据的少量OBN资料与拖缆资料进行联合全波形反演，能够很好地解决拖缆地震资料的低频缺失，压制周期跳跃现象。本文提出了基于平均道方法的海洋拖缆与OBN资料联合波形反演方法。Marmousi模型测试表明，利用本文方法引入OBN资料解决了拖缆资料低频缺失问题，利用平均道褶积方法，使反演不依赖震源子波，压制了子波不匹配引起的反演错误，达到了不依赖震源子波的效果。另外，测试还表明与传统的基于单一参考道反演相比，基于平均道的反演方法在面对含噪数据时具有更好的抗噪性。

因无同时采集的拖缆和OBN资料，未进行实测资料测试。

附录A 不依赖震源子波的目标函数梯度计算

基于平均道的不依赖震源子波的目标函数为

$ E=\sum\limits_{i=1}^{{N}_{\mathrm{s}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}{‖{\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}‖}_{2}^{2} $

(A-1)

将上式对模型求偏导可得

$ \frac{\partial E}{\partial \boldsymbol{m}}=\sum\limits_{i=1}^{{N}_{\mathrm{s}}}\sum\limits_{j=1}^{{N}_{\mathrm{r}}}\left[\left(\frac{\partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right)\cdot {\boldsymbol{r}}_{i, j}-\left({\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}\frac{\partial {\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\right)\cdot {\boldsymbol{r}}_{i, j}\right] $

(A-2)

式中：m为模型参数；$ {\boldsymbol{r}}_{i, j}={\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}} $。将上式右端第一项写成积分形式

$ \begin{array}{l}\left(\frac{\partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right)\cdot {\boldsymbol{r}}_{i, j}\\ \;\;\;\;\;\;\;\;={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\left[\frac{\partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\left(t-\tau \right){\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\left(\tau \right){\boldsymbol{r}}_{i, j}\left(t\right)\right]\mathrm{d}\tau \mathrm{d}t\end{array} $

(A-3)

令$ t-\tau =\xi $，则有$ \mathrm{d}\tau =-\mathrm{d}\xi $，上式可以写成

$ \left(\frac{\partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right)\cdot {\boldsymbol{r}}_{i, j}=-{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{\partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\left(\xi \right){\boldsymbol{r}}_{i, j}^{\text{'}}\left(\xi \right)\mathrm{d}\xi $

(A-4)

式中：$ {\boldsymbol{r}}_{i, j}^{\text{'}}\left(\xi \right) $是地震数据残差$ {\boldsymbol{r}}_{i, j} $与观测数据参考道进行互相关的结果；$ \partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}/\partial m $可以根据虚拟源理论由正传格林函数$ {\boldsymbol{G}}_{x, j}^{\mathrm{c}\mathrm{a}\mathrm{l}} $与虚拟源$ {\boldsymbol{V}}_{i, x}^{\mathrm{c}\mathrm{a}\mathrm{l}} $的卷积获得，因此有

$ \begin{array}{l}\left(\frac{\partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right)\cdot {\boldsymbol{r}}_{i, j}\\ =-{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\left[{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\boldsymbol{V}}_{i, x}^{\mathrm{c}\mathrm{a}\mathrm{l}}\left(\xi -\tau \right){\boldsymbol{G}}_{x, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\left(\tau \right)\mathrm{d}\tau \right]{\boldsymbol{r}}_{i, j}^{\text{'}}\left(\xi \right)\mathrm{d}\xi \end{array} $

(A-5)

再令$ \xi -\tau =t $，则有$ \mathrm{d}\tau =-\mathrm{d}t $，上式变为

$ \begin{array}{l}\left(\frac{\partial {\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right)\cdot {\boldsymbol{r}}_{i, j}\\ ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\boldsymbol{V}}_{i, \boldsymbol{x}}^{\mathrm{c}\mathrm{a}\mathrm{l}}\left(t\right)\left[{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\boldsymbol{G}}_{\boldsymbol{x}, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\left(\xi -t\right){\boldsymbol{r}}_{i, j}^{\text{'}}\left(\xi \right)\mathrm{d}\xi \right]\mathrm{d}t\end{array} $

(A-6)

根据互易定理，格林函数$ {\boldsymbol{G}}_{\boldsymbol{x}, j}^{\mathrm{c}\mathrm{a}\mathrm{l}} $满足$ {\boldsymbol{G}}_{\boldsymbol{x}, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}={\boldsymbol{G}}_{j, \boldsymbol{x}}^{\mathrm{c}\mathrm{a}\mathrm{l}} $，因此可以将$ {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\boldsymbol{G}}_{\boldsymbol{x}, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\left(\xi -t\right){\boldsymbol{r}}_{i, j}^{\text{'}}\left(\xi \right)\mathrm{d}\xi $视为残差波场在检波器位置以时间变量t反传的过程。类似地，对式(A-2)$ \mathrm{等} $式右端第二项进行类似的推导，可以得到

$ \begin{array}{l}-\left({\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}\frac{\partial {\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}}{\partial \boldsymbol{m}}\right)\cdot {\boldsymbol{r}}_{i, j}\\ ={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\boldsymbol{V}}_{i, \boldsymbol{x}}^{\mathrm{c}\mathrm{a}\mathrm{l}}\left(t\right)\left[{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\boldsymbol{G}}_{\boldsymbol{x}, k}^{\mathrm{c}\mathrm{a}\mathrm{l}}\left(\xi -t\right){\boldsymbol{r}}_{i, j}^{″}\left(\xi \right)\mathrm{d}\xi \right]\mathrm{d}t\end{array} $

(A-7)

式中$ {\boldsymbol{r}}_{i, j}^{″}\left(\xi \right)=-{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\left[{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\left(t-\xi \right)\cdot {\boldsymbol{r}}_{i, j}\left(t\right)\right]\mathrm{d}t $。

综上所述，求取梯度所需要的反传震源可以根据$ {\boldsymbol{r}}_{i, j}^{\text{'}} $和$ {\boldsymbol{r}}_{i, j}^{″} $得到

$ {\boldsymbol{r}}_{i, j}^{\text{'}}={\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}\otimes \left({\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\right) $

(A-8)

$ {\boldsymbol{r}}_{i, j}^{″}=-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\otimes \left({\boldsymbol{d}}_{i, j}^{\mathrm{c}\mathrm{a}\mathrm{l}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\boldsymbol{d}}_{i, j}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{*}{\boldsymbol{A}}_{i}^{\mathrm{c}\mathrm{a}\mathrm{l}}\right) $

(A-9)