2. 同济大学海洋与地球科学学院波现象与智能反演成像研究组, 上海 200092;
3. 同济大学海洋高等研究院, 上海 200092
2. Wave Phenomena and Intelligent Inversion Imaging Group (WPI), School of Ocean and Earth Science, Tongji University, Shanghai 200092, China;
3. Institute for Advanced Study, Tongji University, Shanghai 200092 China
地震波走时多用于中、大尺度的速度结构特征反演。对于有限频带的地震信号, 地震波走时(绝对走时)受震源子波的初相位、走时测量准则(如起跳时刻、包络极值时刻, 最大峰值时刻)等多种因素的影响。地震波走时的测量误差会降低反演结果的精度及可靠性。为尽可能降低绝对走时测量误差对反演的影响, BRUNE等[1]提出采用“相对测量值”代替“绝对测量值”以减少测量误差。基于此方法测量的地震走时亦称为双差走时(double-difference traveltime), 可以用于提高震源定位精度[2-4], 反演高精度速度模型[5]。
传统的双差走时层析反演方法多基于射线理论, 存在焦散及阴影区等问题[6], 且反演精度不高。对于小尺度异常体(速度非均匀体的尺度小于菲涅尔体的宽度), 有限频理论[7-14]可以更好地处理地震波的一阶绕射效应[15-16]。YUAN等[17]将双差走时测量方法引入伴随层析(adjoint tomography), 得到了更高精度的反演结果, 但双差走时敏感度核函数仍采用Born近似计算方法计算得到[11]。相较于Born近似, Rytov近似可以更好地描述由于速度扰动引起的前向散射波相位扰动[18], 因而更适用于走时反演。本文利用Rytov近似构造走时敏感度核函数[19-26], 即将双差走时测量方法与Rytov走时敏感度核函数相结合, 提出了一种新的波动方程初至层析反演方法。该方法在保留Rytov近似优点的同时, 可以降低甚至消除走时测量误差对反演结果的影响。
1 理论 1.1 双差相位延迟敏感度核函数频率域标量声波方程可以表示为:
| $ \left(\nabla^{2}+k^{2}\right) u_{0}=0 $ | (1) |
式中: k(x)为背景速度场v0(x)中的波数, k(x)=ω/v0(x); u0(x; ω)为背景场; ▽2为Laplace算子。
令扰动后的速度场为v(x), 总场为u(x; ω), 总场与背景场的关系表示为:
| $ \frac{u(\boldsymbol{x} ; \omega)}{u_{0}(\boldsymbol{x} ; \omega)}=\mathrm{e}^{\psi(\boldsymbol{x} ; \omega)} $ | (2) |
式中: ψ(x; ω)为复相位。
利用Rytov近似, 将单频谐波的相位延迟表示为:
| $ \begin{array}{c} \Delta t^{\mathrm{p}}\left(x_{\mathrm{r}}, \omega ; x_{\mathrm{s}}\right)=-\operatorname{Im}\left[\frac{\psi\left(x_{\mathrm{r}}, \omega ; x_{\mathrm{s}}\right)}{\omega}\right] \\ =\left\langle K^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{\mathrm{r}}, x_{\mathrm{s}}\right), \Delta m(\boldsymbol{x})\right\rangle_{\mathrm{M}} \end{array} $ | (3) |
式中: Im表示取复数的虚部; Δm(x)为慢度平方扰动(属于模型空间M), Δm(x)=v(x)-2-v0(x)-2; xs, xr分别代表震源和检波器横坐标; Kp(x, ω; xr, xs)为相位延迟敏感度核函数[25]。Kp(x, ω; xr, xs)可表示为:
| $ \begin{array}{l} K^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{\mathrm{r}}, x_{\mathrm{s}}\right)= \\ \quad \operatorname{Im}\left[\frac{-\omega G_{0}\left(\boldsymbol{x}, \omega ; x_{\mathrm{r}}\right) u_{0}\left(\boldsymbol{x}, \omega ; x_{\mathrm{s}}\right)}{u_{0}\left(\boldsymbol{x}_{\mathrm{r}}, \omega ; x_{\mathrm{s}}\right)}\right] \end{array} $ | (4) |
(3) 式所预测的相位延迟假定了u和u0由相同的震源产生。若震源子波未知, 复相位中包含子波剩余相位差。在全波形反演中, 通常用褶积[27]或反褶积[28-31]类方法消除震源子波对反演的影响。本文采用反褶积方法, 引入如下反褶积波场:
| $ v_{i, j}\left(\omega ; x_{\mathrm{s}}\right)=\frac{u^{w}\left(x_{i}, \omega ; x_{\mathrm{s}}\right)}{u^{w}\left(x_{j}, \omega ; x_{\mathrm{s}}\right)} $ | (5) |
式中: xi和xj为同一炮集中任意不重合的两个检波器横坐标; uw(xi, ω; xs)和uw(xj, ω; xs)分别为初至波形加窗前、后的地震信号(频率域)。uw(xi, ω; xs)可表示为:
| $ \begin{array}{c} u^{w}\left(x_{i}, \omega ; x_{\mathrm{s}}\right)=\int_{0}^{\mathrm{T}} h(t) u\left(x_{i}, t ; x_{\mathrm{s}}\right) \mathrm{e}^{i \omega t} \mathrm{~d} t \\ \equiv h(\omega) * u\left(x_{i}, \omega ; x_{\mathrm{s}}\right) \end{array} $ | (6) |
式中: h(t)为初至波形窗函数, 用于提取初至震相; *为卷积运算符号; h(ω)为窗函数的频率响应。
令uw(xi, ω; xs)=f(ω)A(xi, ω; xs)eiφ(xi, ω; xs), uw(xj, ω; xs)=f(ω)A(xj, ω; xs)eiφ(xj, ω; xs), 其中, f(ω)为震源子波频谱, A(ω), φ(ω)分别为格林函数的振幅谱和相位谱, 初至波形的相位差Δφi, j(ω; xs)可以表示为:
| $ \begin{array}{c} \Delta \varphi_{i, j}\left(\omega ; x_{\mathrm{s}}\right)=\varphi\left(x_{i}, \omega ; x_{\mathrm{s}}\right)-\varphi\left(x_{j}, \omega ; x_{\mathrm{s}}\right) \\ \equiv \operatorname{Im}\left[\ln v_{i, j}\left(\omega ; x_{\mathrm{s}}\right)\right] \end{array} $ | (7) |
显然, 子波频率对地震信号相位的影响被抵消了。采用双差测量方法[5], 得到的观测数据和模拟数据的双差相位延迟可以表示为:
| $ \begin{array}{l} \Delta \Delta t_{i, j}^{\mathrm{p}}\left(\omega ; x_{\mathrm{s}}\right)=\frac{1}{\omega}\left[\Delta \varphi_{i, j}^{\mathrm{obs}}\left(\omega ; x_{\mathrm{s}}\right)-\Delta \varphi_{i, j}^{\text {cal }}\left(\omega ; x_{\mathrm{s}}\right)\right] \\ =\frac{1}{\omega} \operatorname{Im}\left[\ln \frac{u^{w}\left(x_{i}, \omega ; x_{\mathrm{s}}\right)}{u_{0}^{w}\left(x_{i}, \omega ; x_{\mathrm{s}}\right)}-\ln \frac{u^{w}\left(x_{j}, \omega ; x_{\mathrm{s}}\right)}{u_{0}^{w}\left(x_{j}, \omega ; x_{\mathrm{s}}\right)}\right] \end{array} $ | (8) |
式中: Δφi, jobs(ω; xs)和Δφi, jcal(ω; xs)分别为观测数据和模拟数据的初至波相位差。
综合考虑(2)式至(7)式, 初至波形的相位延迟和模型扰动的线性关系可表示为:
| $ \begin{array}{c} -\frac{1}{\omega} \operatorname{Im} \ln \frac{u^{w}\left(x_{\mathrm{r}}, \omega ; x_{\mathrm{s}}\right)}{u_{0}^{w}\left(x_{\mathrm{r}}, \omega ; x_{\mathrm{s}}\right)}=\left\langle K^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{\mathrm{r}}, x_{\mathrm{s}}\right),\right. \\ \Delta m(\boldsymbol{x})\rangle_{\mathrm{M}} \end{array} $ | (9) |
将(9)式代入(8)式, 双差相位延迟与模型扰动的线性关系可表示为:
| $ \Delta \Delta t_{i, j}^{\mathrm{p}}\left(\omega ; x_{\mathrm{s}}\right)=\int_{V} \Delta m(\boldsymbol{x}) K_{\mathrm{dd}}^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{i}, x_{j}, x_{\mathrm{s}}\right) \mathrm{d} x $ | (10) |
式中: Kddp(x, ω; xi, xj, xs)为双差相位延迟敏感度核函数。Kddp(x, ω; xi, xj, xs)可表示为:
| $ \begin{array}{c} K_{\mathrm{dd}}^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{i}, x_{j}, x_{\mathrm{s}}\right)=K^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{i}, x_{\mathrm{s}}\right)- \\ K^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{j}, x_{\mathrm{s}}\right) \end{array} $ | (11) |
有限频带地震信号的走时扰动可以表示为单频谐波的相位延迟的加权叠加[22], 即:
| $ \Delta t=\int W(\omega) \Delta t_{p}(\omega) \mathrm{d} \omega $ | (12) |
式中: W(ω)为归一化加权函数, W(ω)=ω2P(ω)/∫ω2P(ω)dω; P(ω)为震源子波能谱, P(ω)=f(ω)f*(ω)。
同理, 双差走时扰动与双差相位延迟的关系可表示为:
| $ \Delta \Delta t_{i, j}\left(x_{\mathrm{s}}\right)=\int W(\omega) \Delta \Delta t_{i, j}^{\mathrm{p}}\left(\omega ; x_{\mathrm{s}}\right) \mathrm{d} \omega $ | (13) |
将(10)式、(11)式与(13)式相结合, 可以建立双差走时扰动与模型扰动的线性关系:
| $ \Delta \Delta t_{i, j}\left(x_{\mathrm{s}}\right)=\int_{V} \Delta m(\boldsymbol{x}) K_{\mathrm{dd}}\left(\boldsymbol{x} ; x_{i}, x_{j}, x_{\mathrm{s}}\right) \mathrm{d} x $ | (14) |
式中: Kdd(x; xi, xj, xs)为(带限信号或有限频)双差走时敏感度核函数。Kdd(x; xi, xj, xs)的表达式如下:
| $ \begin{array}{c} K_{\mathrm{dd}}\left(\boldsymbol{x} ; x_{i}, x_{j}, x_{\mathrm{s}}\right)=\int W(\omega)\left[K^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{i}, x_{\mathrm{s}}\right)-\right. \\ \left.K^{\mathrm{p}}\left(\boldsymbol{x}, \omega ; x_{j}, x_{\mathrm{s}}\right)\right] \mathrm{d} \omega=K\left(\boldsymbol{x} ; x_{i}, x\right)- \\ K\left(\boldsymbol{x} ; x_{j}, x_{\mathrm{s}}\right) \end{array} $ | (15) |
式中: K(x; xr, xs)为基于Rytov近似的有限频走时敏感度核函数[24-25]。
冯波等[25]给出了Rytov近似走时敏感度核函数K(x; xr, xs)在时间-空间域的显式计算公式为:
| $ K\left(\boldsymbol{x} ; x_{\mathrm{r}}, x_{\mathrm{s}}\right)=\int_{0}^{T} \partial_{t}^{2} u_{0}\left(x, t ; x_{\mathrm{s}}\right) \tilde{\lambda}_{0}\left(\boldsymbol{x}, T-t ; x_{\mathrm{r}}\right) \mathrm{d} t $ | (16) |
式中:
基于走时残差L2范数的误差泛函可以表示为:
| $ \min \limits_{m \in M} J(\boldsymbol{m})=\frac{1}{2} \sum\limits_{x_{\mathrm{s}}} \sum\limits_{j \neq i}\left[\Delta \Delta t_{i, j}\left(x_{\mathrm{s}}\right)\right]^{2} $ | (17) |
假定每个震源有Nr道地震记录, 随机采用第i个地震道作为参考道, 其它地震道与参考道计算得到的双差走时为ΔΔti, j, 满足1≤j≤Nr且j≠i。
采用Gauss-Newton反演算法, 将模型参数更新过程表示为:
| $ \boldsymbol{m}^{k+1}=\boldsymbol{m}^{k}-\alpha^{k} P\left[\boldsymbol{H}^{-1}\left(m^{k}\right) g\left(m^{k}\right)\right] $ | (18) |
式中: αk为步长, 可以通过线性搜索方法计算得到; g(mk)为目标函数梯度, g(mk)=KddTΔΔt; H为线性Hessian矩阵, H=KddTKdd; P为光滑算子(本文采用多尺度高斯平滑算子)。
由(15)式可知, 泛函梯度g(mk)可以用核函数与走时残差表示为:
| $ \begin{array}{c} g\left(\boldsymbol{m}^{k}\right)=\boldsymbol{K}_{\mathrm{dd}}^{\mathrm{T}} \Delta \Delta t^{k}=\sum\limits_{x_{s}} \sum\limits_{j \neq i} K_{\mathrm{dd}}\left(\boldsymbol{x} ; x_{i}, x_{j}, x_{\mathrm{s}}\right) \cdot \\ \Delta \Delta t_{i, j}^{k}\left(x_{\mathrm{s}}\right)=\sum\limits_{x_{\mathrm{s}}}\left[K_{\mathrm{dd}}\left(\boldsymbol{x} ; x_{i}, x_{\mathrm{s}}\right) \sum\limits_{j \neq i} \Delta \Delta t_{i, j}^{k}\left(x_{\mathrm{s}}\right)-\right. \\ \left.\sum\limits_{j \neq i} K_{\mathrm{dd}}\left(\boldsymbol{x} ; x_{j}, x_{\mathrm{s}}\right) \Delta \Delta t_{i, j}^{k}\left(x_{\mathrm{s}}\right)\right] \end{array} $ | (19) |
若定义双差走时伴随震源f†dd(x, t; xs)为:
| $ \begin{array}{c} f_{\mathrm{dd}}^{\dagger}\left(\boldsymbol{x}, t ; x_{\mathrm{s}}\right)=\delta\left(\boldsymbol{x}-x_{i}\right) \tilde{f}^{\dagger}\left(x_{i}, t ; x_{\mathrm{s}}\right) \sum\limits_{j \neq i} \Delta \Delta t_{i, j} \cdot \\ \left(x_{\mathrm{s}}\right)-\sum\limits_{j \neq i} \delta\left(\boldsymbol{x}-x_{j}\right) \tilde{f}^{\dagger}\left(x_{j}, t ; x_{\mathrm{s}}\right) \Delta \Delta t_{i, j}\left(x_{\mathrm{s}}\right) \end{array} $ | (20) |
(19) 式可以表示为:
| $ g\left(\boldsymbol{m}^{k}\right)=\sum\limits_{x_{\mathrm{s}}} \int_{0}^{T} \partial_{t}^{2} u_{0}\left(\boldsymbol{x}, t ; x_{\mathrm{s}}\right) \lambda_{\mathrm{dd}}\left(\boldsymbol{x}, T-t ; x_{\mathrm{s}}\right) \mathrm{d} t $ | (21) |
式中: λdd(x, T－t; xs)为双差走时伴随震源逆时传播产生的伴随波场。
为避免直接Hessian矩阵求逆, 可以通过求解如下方程的近似解获得模型更新方向:
| $ \boldsymbol{H}\left(\boldsymbol{m}^{k}\right) \Delta \boldsymbol{m}^{k}=g\left(\boldsymbol{m}^{k}\right) $ | (22) |
本文采用共轭梯度(conjugate gradient, CG)法求解方程, 采用冯波等[25]提出的隐式矩阵向量乘法得到高效的Hessian向量积, 无需显式计算及存储Hessian矩阵。具体计算公式见附录A。
2 数值试验 2.1 高斯扰动模型测试我们设计了一个含有高斯异常体的速度模型v(x, z)=v0+δv(x, z), 其中背景模型v0为匀速模型, 高斯异常体δv描述如下:
| $ \left\{\begin{array}{l} \delta v(x, z)=\varepsilon v_{0} \exp \left(-\frac{r^{2}}{2 a^{2}}\right) \\ r=\sqrt{\left(x-x_{0}\right)^{2}+\left(z-z_{0}\right)^{2}} \end{array}\right. $ | (23) |
式中: ε为速度扰动百分比; v0=2500m/s为均匀背景速度; a为高斯异常体的尺度参数; (x0, z0)为高斯异常体的中心坐标, (x0, z0)=(2500m, 2500m)。
采用10m×10m的网格离散化含有高斯异常体的速度模型, x和z方向网格数目均为501。高斯异常体的尺度参数a=500m, 速度扰动百分比ε=10%(图 1a)。为消除采集孔径对反演结果的影响, 本文设计了一个四周激发-接收的观测系统, 在模型四周放置100个均匀分布的震源, 每边25个震源。每炮由100个检波器接收, 检波器均匀分布在模型四周, 并与震源重合。正演采用10Hz主频的Ricker子波(对应波长λ0=250m, 与异常体尺度参数满足a=2λ0)。
|
图 1 含有高斯异常体的速度模型(a)、采用双差走时(b)与绝对走时(c)层析反演方法得到的速度扰动 |
为了验证本文双差走时反演方法的有效性, 我们首先对每个震源子波引入随机延迟, 并采用有限差分方法模拟观测地震记录。本文采用SEISCOPE数值优化软件包(SEISCOPE optimization toolbox[31])中的截断牛顿算法(对于走时层析, 截断牛顿算法退化为高斯牛顿算法)求解双差走时目标函数。根据附录A中的隐式矩阵向量积计算公式计算梯度与Hessian向量积。初始模型采用匀速背景(v0=2500m/s), 震源子波采用0.5s延迟的Ricker子波(10Hz主频)。目标函数的终止准则为检测目标函数的相对误差(σ=J(mk)/J(m0))小于预先给定的门槛值σ=1.0×10-4。步长采用线性搜索估计, 每个高斯牛顿方向采用2次CG内迭代求解(16)式。经过10次高斯牛顿方向迭代求解, 得到的速度扰动(图 1b)最大值为249.7m/s, 与真实速度扰动的最大值250.0m/s基本一致。
接着我们根据相同的反演参数, 利用绝对走时[25]层析反演方法进行反演, 将正演观测数据得到的随机延迟Ricker子波作为震源, 得到的速度扰动如图 1c所示。可以看出, 采用双差走时与绝对走时层析反演方法得到的速度扰动均与真实速度模型吻合较好。因此对于完备的观测系统, 双差走时层析反演与常规反演结果基本一致, 证明了双差走时层析反演方法的有效性。
2.2 Overthrust模型测试为进一步测试本文方法在近地表速度建模中的效果, 我们用SEG/EAGE Overthrust速度模型[32]正演观测地震数据。原始速度模型在x和z方向网格数目分别为801和187, 网格间距为25m。为更好地反演速度模型左侧的推覆构造, 我们将速度模型左右侧分别拓展, 得到一个横向长度为25000m的SEG/EAGE Overthrust近地表速度模型(图 2a)。我们设计了一个陆上单边偏移距观测系统, 在地表放置91个均匀分布的震源, 炮间距为200m, 每炮由117道检波器接收, 道间距为25m(最小偏移距100m, 最大偏移距3000m)。观测地震记录由有限差分算法计算得到, 地震子波主频为8Hz且存在随机延迟的Ricker子波。
|
图 2 横向长度为25000m的SEG/EAGE Overthrust近地表速度模型(a)、采用双差走时层析反演方法得到的速度模型(b)及速度扰动(c) |
由于观测地震记录较为复杂, 本文采用自动初至拾取算法得到观测数据和模拟数据的初至, 并计算其双差走时(参考道为最小偏移距地震道)。初始模型采用线性递增的速度模型, 当z=0时, v=2400m/s, 当z=3000m时, v=4000m/s, 震源子波采用8Hz主频的Ricker子波(子波延迟时间为0.2s), 采用截断牛顿算法计算模型更新方向, 每个方向1次CG内迭代求解(16)式, 并用线性搜索方法估计迭代步长。在反演过程中, 采用了多偏移距反演策略(从最长偏移距(3000m)逐渐减少到最小偏移距(500m)), 最终采用双差走时层析反演方法得到的速度模型如图 2b所示。可以看出, 大尺度的近地表速度结构特征恢复得较好。图 2c为利用该反演方法得到的速度扰动, 主要表现为近地表中-大尺度的速度更新, 最大有效反演深度约为150m。
图 3为不同深度处速度横向抽线对比结果, 显然双差走时层析反演结果在浅层的分辨率较高, 甚至能分辨中、小尺度的速度异常。随着深度增加, 其反演精度逐渐降低(有效反演深度由最大偏移距及速度结构共同决定)。
|
图 3 不同深度处速度横向抽线对比结果 a深度25m; b深度100m |
对于实际地震资料反演而言, 反演效果的优劣受到多种因素综合影响, 如观测系统的完备性、时差测量的精度、反演策略、数据加权矩阵的设计、模型正则化方法及先验信息等。本文主要考虑如何消除测量方式对地震波时差计算的影响。对于复杂的实际地震数据, 本文利用ΔΔti, j=Δti, jobs-Δti, jcal公式近似计算双差走时, 其中Δti, jobs和Δti, jcal分别为观测数据与模拟数据基于相同走时测量准则得到的相对时差。
采用中国某地二维陆上实际地震资料(共1638炮, 双边接收, 最大偏移距7132m;炮间隔60m, 道间距10m)进行测试, 该地震资料采用初至自动拾取方法得到(拾取起震时刻)。初始模型采用线性递增速度模型(起伏地表之上用空气速度填充), 震源子波为主频为10Hz的Ricker子波, 基于相同的初至拾取标准得到模拟数据的初至, 并采用多偏移距反演策略(最大偏移距分别为3000m, 2000m, 1000m, 500m, 每个尺度迭代8次), 得到如图 4所示的最终反演结果, 可以看出最大有效反演深度约为700m。从反演结果可以看出, 山间低速带得到了较好的刻画(近地表速度低至800m/s, 低速带厚度约为200m), 这与野外实际地表露头结果相符, 速度模型左侧的平原地带, 反演结果表现出明显的成层性。
|
图 4 二维陆上实际地震资料反演结果 |
为减小未知地震子波波形(或延迟时)及走时测量方法对(绝对)初至时间检测带来的误差, 本文对观测数据和模拟数据采用相同的初至拾取方法并计算双差走时, 以降低甚至消除走时测量误差对反演结果的影响。虽然根据公式推导, 严格的双差走时应该采用加权相位延迟计算得到, 但在实际应用中直接采用(自动或人工)拾取的初至计算双差走时得到的结果仍然是稳健的。
相较于传统波动方程走时反演方法中引入的一阶Born近似(要求速度异常体的尺度和扰动强度都足够小), 本文结合双差走时测量及Rytov近似构建相位延迟敏感度核函数, 可以更好地预测由于(大尺度)速度扰动引起的前向散射波的相位扰动, 因此降低了对初始模型精度的要求。
在反问题数值求解方面, 本文由于引入了基于隐式矩阵向量积的高斯-牛顿迭代算法, 故仅需波动方程Born模拟及逆时偏移算法即可高效计算高斯-牛顿搜索方向, 无需显式计算和存储Hessian矩阵, 因此适用于求解大规模计算问题。
附录A 隐式矩阵向量乘法Hessian向量积表示为: Hp=KddTKddp=KddTh, (其中, h=Kddp)因此需要分别计算如下两类矩阵向量积。
1) 矩阵向量积Kddp(p为模型空间中的向量p=p(x))。
根据(15)式, Kdd中任意一行与模型空间中的向量p的内积可以表示为:
| $ \left(K_{\mathrm{dd}} \boldsymbol{p}\right)\left(x_{i}, x_{j}, x_{\mathrm{s}}\right)=\left\langle K\left(\boldsymbol{x} ; x_{i}, x_{\mathrm{s}}\right)-K\left(\boldsymbol{x} ; x_{j}, x_{\mathrm{s}}\right), \boldsymbol{p}\right\rangle_{M}=\left\langle K\left(\boldsymbol{x} ; x_{i}, x_{\mathrm{s}}\right), \boldsymbol{p}\right\rangle_{M}-\left\langle K\left(\boldsymbol{x} ; x_{j}, x_{\mathrm{s}}\right), \boldsymbol{p}\right\rangle_{M} $ | (A1) |
冯波等[25]证明, 走时敏感度核函数与模型空间向量的内积可以转化为波场空间中的内积:
| $ (K \boldsymbol{p})\left(x_{\mathrm{r}}, x_{\mathrm{s}}\right)=\left\langle u_{q}\left(x_{\mathrm{r}}, t ; x_{\mathrm{s}}\right), u_{p}\left(x_{\mathrm{r}}, t ; x_{\mathrm{s}}\right)\right\rangle_{T} $ | (A2) |
将(A2)代入(A1), 有:
| $ \left(K_{\mathrm{dd}} \boldsymbol{p}\right)\left(x_{i}, x_{j}, x_{\mathrm{s}}\right)=\left\langle u_{q}\left(x_{i}, t ; x_{\mathrm{s}}\right), u_{p}\left(x_{i}, t ; x_{\mathrm{s}}\right)\right\rangle_{T}-\left\langle u_{q}\left(x_{j}, t ; x_{\mathrm{s}}\right), u_{p}\left(x_{j}, t ; x_{\mathrm{s}}\right)\right\rangle_{T} $ | (A3) |
2) 矩阵向量积KddTh(h为数据空间中的向量h=h(xi, xj, xs))。
根据(21)式, KddTh可以表示为:
| $ \boldsymbol{K}_{\mathrm{dd}}^{\mathrm{T}} \boldsymbol{h}=\sum\limits_{x_{\mathrm{s}}} \int_{0}^{T} \partial_{t}^{2} u_{0}\left(\boldsymbol{x}, t ; x_{\mathrm{s}}\right) \lambda_{\mathrm{dd}}\left(\boldsymbol{x}, T-t ; x_{\mathrm{s}}\right) \mathrm{d} t $ | (A4) |
其中λdd(x, T－t; xs)由如下双差走时伴随震源产生:
| $ f_{\mathrm{dd}}^{\dagger}\left(\boldsymbol{x}, t ; x_{\mathrm{s}}\right)=\delta\left(\boldsymbol{x}-x_{i}\right) \tilde{f}^{\dagger}\left(x_{i}, t ; x_{\mathrm{s}}\right) \sum\limits_{j \neq i} h\left(x_{i}, x_{j}, x_{\mathrm{s}}\right)-\sum\limits_{j \neq i} \delta\left(\boldsymbol{x}-x_{j}\right) \tilde{f}^{\dagger}\left(x_{j}, t ; x_{\mathrm{s}}\right) h\left(x_{i}, x_{j}, x_{\mathrm{s}}\right) $ | (A5) |
顺序计算(A3)及(A4)式, 无需显式计算和存储敏感度核函数及Hessian矩阵, 可以直接获得Hessian向量积。
| [1] |
BRUNE J, DORMAN J. Seismic waves and earth structure in the Canadian shield[J]. Bulletin of the Seismological Society of America, 1963, 53(1): 167-210. |
| [2] |
WALDHAUSER F, ELLSWORTH W L. A double-difference earthquake location algorithm: Method and application to the Northern Hayward Fault[J]. Bulletin of the Seismological Society of America, 2000, 90(6): 1353-1368. DOI:10.1785/0120000006 |
| [3] |
ZHAO L F, XIE X B, WANG W, et al. The 9 September 2016 North Korean underground nuclear test[J]. Bulletin of the Seismological Society of America, 2017, 107(6): 1-8. |
| [4] |
谢小碧, 赵连锋. 朝鲜地下核试验的地震学观测[J]. 地球物理学报, 2018, 61(3): 889-904. XIE X B, ZHAO L F. The seismic characterization of North Korea underground nuclear tests[J]. Chinese Journal of Geophysics, 2018, 61(3): 889-904. |
| [5] |
ZHANG H, THURBER C H. Double-difference tomography: The method and its application to the Hayward Fault[J]. Bulletin of the Seismological Society of America, 2003, 93(5): 1875-1889. DOI:10.1785/0120020190 |
| [6] |
罗飞, 王华忠, 冯波, 等. 透射波旅行时Beam层析成像方法[J]. 石油物探, 2019, 58(3): 356-370. LUO F, WANG H Z, FENG B, et al. Beam tomography based on transmission traveltime[J]. Geophysical Prospecting for Petroleum, 2019, 58(3): 356-370. DOI:10.3969/j.issn.1000-1441.2019.03.005 |
| [7] |
LUO Y, SCHUSTER G T. Wave-equation travel time inversion[J]. Geophysics, 1991, 56(5): 645-653. DOI:10.1190/1.1443081 |
| [8] |
MARQUERING H, DAHLEN F A, NOLET G. Three-dimensional sensitivity kernels for finite-frequency traveltimes: The banana-doughnut paradox[J]. Geophysical Journal International, 1999, 137(3): 805-815. DOI:10.1046/j.1365-246x.1999.00837.x |
| [9] |
DAHLEN F A, HUNG S H, NOLET G. Fréchet kernels for finite-frequency traveltimes-Ⅰ.Theory[J]. Geophysical Journal International, 2000, 141(1): 157-175. DOI:10.1046/j.1365-246X.2000.00070.x |
| [10] |
HUNG S H, DAHLEN F A, NOLET G. Fréchet kernels for finite-frequency traveltime-Ⅱ.Examples[J]. Geophysical Journal International, 2000, 141(1): 175-203. DOI:10.1046/j.1365-246X.2000.00072.x |
| [11] |
TROMP J, TAPE C, LIU Q. Seismic tomography, adjoint methods, time reversal, and banana-doughnut kernels[J]. Geophysical Journal International, 2005, 160(1): 195-216. |
| [12] |
LIU Q, TROMP J. Finite-frequency kernels based upon adjoint methods[J]. Bulletin of the Seismological Society of America, 2006, 96(6): 2383-2397. DOI:10.1785/0120060041 |
| [13] |
LIU Q, TROMP J. Finite-frequency sensitivity kernels for global seismic wave propagation based upon adjoint methods[J]. Geophysical Journal International, 2008, 174(1): 265-286. DOI:10.1111/j.1365-246X.2008.03798.x |
| [14] |
TAPE C, LIU Q, TROMP J. Finite frequency tomography using adjoint methods—methodology and examples using membrane surface waves[J]. Geophysical Journal International, 2007, 168(3): 1105-1129. DOI:10.1111/j.1365-246X.2006.03191.x |
| [15] |
SLANEY M, KAK A C, LARSEN L. Limitations of imaging with first-order diffraction tomography[J]. IEEE Transactions On Microwave Theory And Techniques, 1984, 32(8): 860-873. DOI:10.1109/TMTT.1984.1132783 |
| [16] |
WU R S, TOKSÖZ M N. Diffraction tomography and multisource holography applied to seismic imaging[J]. Geophysics, 1987, 52(1): 11-25. DOI:10.1190/1.1442237 |
| [17] |
YUAN Y O, SIMONS F J, TROMP J. Double-difference adjoint seismic tomography[J]. Geophysical Journal International, 2016, 206(3): 1599-1618. DOI:10.1093/gji/ggw233 |
| [18] |
WU R S. Wave propagation, scattering and imaging using dual-domain one-way and one-return propagators[J]. Pure and Applied Geophysics, 2003, 160(3/4): 509-539. DOI:10.1007/PL00012548 |
| [19] |
SPETZLER J, SNIEDER R. The effects of small-scale heterogeneity on the arrival time of waves[J]. Geophysical Journal International, 2001, 145(3): 786-796. DOI:10.1046/j.1365-246x.2001.01438.x |
| [20] |
SPETZLER J, SNIEDER R. The Fresnel volume and transmitted waves[J]. Geophysics, 2004, 69(3): 653-663. DOI:10.1190/1.1759451 |
| [21] |
JOCKER J, SPETZLER J, SMEULDERS D, et al. Validation of first-order diffraction theory for the traveltimes and amplitudes of propagating waves[J]. Geophysics, 2006, 71(6): T167-T177. DOI:10.1190/1.2358412 |
| [22] |
XIE X B, YANG H. The finite-frequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis[J]. Geophysics, 2008, 73(6): S241-S249. DOI:10.1190/1.2993536 |
| [23] |
刘玉柱, 董良国, 李培明, 等. 初至波菲涅尔体地震层析成像[J]. 地球物理学报, 2009, 52(9): 2310-2320. LIU Y Z, DONG L G, LI P M, et al. Fresnel volume tomography based on the first arrival of the seismic wave[J]. Chinese Journal of Geophysics, 2009, 52(9): 2310-2320. DOI:10.3969/j.issn.0001-5733.2009.09.015 |
| [24] |
XU W, XIE X B, GENG J. Validity of the Rytov approximation in the form of finite-frequency sensitivity kernels[J]. Pure and Applied Geophysics, 2015, 172(6): 1409-1427. DOI:10.1007/s00024-014-0991-8 |
| [25] |
冯波, 罗飞, 王华忠. 一阶Rytov近似有限频旅行时层析[J]. 地球物理学报, 2019, 62(6): 2217-2226. FENG B, LUO F, WANG H Z. Wave equation traveltime tomography using Rytov approximation[J]. Chinese Journal of Geophysics, 2019, 62(6): 2217-2226. |
| [26] |
FENG B, XU W J, LUO F, et al. Rytov-approximation-based wave-equation traveltime tomography[J]. Geophysics, 2020, 85(3): R289-R297. DOI:10.1190/geo2019-0210.1 |
| [27] |
CHOI Y, ALKHALIFAH T. Source-independent time-domain waveform inversion using convolved wavefields: Application to the encoded multisource waveform inversion[J]. Geophysics, 2011, 76(5): R125-R134. DOI:10.1190/geo2010-0210.1 |
| [28] |
CHOI Y, SHIN C, MIN D J, et al. Efficient calculation of the steepest descent direction for source-independent seismic waveform inversion: an amplitude approach[J]. Journal of Computational Physics, 2005, 208: 455-468. DOI:10.1016/j.jcp.2004.09.019 |
| [29] |
CHOI Y, MIN D J. Source-independent elastic waveform inversion using a logarithmic wavefield[J]. Journal of Applied Geophysics, 2012, 76: 13-22. DOI:10.1016/j.jappgeo.2011.10.013 |
| [30] |
SNIEDER R, LOMAX A. Wavefield smoothing and the effect of rough velocity perturbations on arrival times and amplitudes[J]. Geophysical Journal International, 1996, 125(3): 796-812. DOI:10.1111/j.1365-246X.1996.tb06024.x |
| [31] |
MÉTIVIER L, BROSSIER R. The SEISCOPE optimization toolbox: A large-scale nonlinear optimization library based on reverse communication[J]. Geophysics, 2016, 81(2): F1-F15. DOI:10.1190/geo2015-0031.1 |
| [32] |
AMINZADEH F, JEAN B, KUNZ T. 3-D Salt and overthrust models[J]. Society of Exploration Geophysicists 3-D Modeling Series, 1997, 1(1): 247-256. |