地震反演

约 9218 字大约 31 分钟

2026-05-28

地震反演的基本概念

地震学中的正演问题是：给定震源、地下介质模型和边界条件，计算地震波走时、振幅或完整波形。反演问题则是利用观测到的地震数据，估计地下介质或震源参数。

一般可以写成

\mathbf{d}=\mathbf{g}(\mathbf{m})

其中 $\mathbf{d}$ 是数据向量， $\mathbf{m}$ 是模型参数， $\mathbf{g}$ 是正演算子。

观测数据 $\mathbf{d}$	模型参数 $\mathbf{m}$
P、S 波走时	P、S 波速度
面波相速度、群速度	剪切波速度结构
波形、振幅和相位	速度、密度、衰减
接收函数	间断面深度、波速比
剪切波分裂	各向异性参数

若数据和模型之间是线性关系，则

\mathbf{d}=\mathbf{G}\mathbf{m}

其中 $\mathbf{G}$ 为正演矩阵，又称设计矩阵或敏感矩阵。

大多数地震反演问题是非线性的。选取参考模型 $\mathbf{m}_0$ ，在其附近作一阶泰勒展开：

\mathbf{g}(\mathbf{m}_0+\delta\mathbf{m}) \approx \mathbf{g}(\mathbf{m}_0) +\mathbf{G}\delta\mathbf{m}

G_{ij} = \left. \frac{\partial g_i}{\partial m_j} \right|_{\mathbf{m}_0}

定义数据残差

\delta\mathbf{d} = \mathbf{d}^{\mathrm{obs}} -\mathbf{g}(\mathbf{m}_0)

可以得到线性化反演方程

\delta\mathbf{d} = \mathbf{G}\delta\mathbf{m}

其中 $\mathbf{G}$ 描述每一个数据对不同模型参数的敏感程度。连续介质中常将其写成敏感核：

\delta d_i = \int_V K_i(\mathbf{x})\delta m(\mathbf{x}) \mathrm dV

超定、欠定与混定问题

设数据个数为 $N$ ，模型参数个数为 $M$ 。

当 $N>M$ 时，称为超定问题，通常不存在严格满足所有数据的模型，需要寻找最佳拟合解；
当 $N<M$ 时，称为欠定问题，能够拟合数据的模型不唯一；
即使 $N>M$ ，若某些模型参数之间存在强相关，或数据对部分模型参数不敏感，问题仍可能表现为混定问题。

地震反演通常同时具有数据有限、覆盖不均、噪声存在和正演模型不准确等问题，因此反演结果一般具有非唯一性。

最小二乘反演

普通最小二乘

对于线性方程

\mathbf{d} = \mathbf{G}\mathbf{m} +\boldsymbol{\varepsilon}

定义数据残差

\mathbf{r} = \mathbf{G}\mathbf{m}-\mathbf{d}

最小二乘反演使目标函数

\Phi_d(\mathbf{m}) = \|\mathbf{G}\mathbf{m}-\mathbf{d}\|_2^2

达到最小。

令目标函数对 $\mathbf{m}$ 的梯度为零，得到正规方程

\mathbf{G}^{\mathrm T}\mathbf{G}\mathbf{m} = \mathbf{G}^{\mathrm T}\mathbf{d}

若 $\mathbf{G}^{\mathrm T}\mathbf{G}$ 可逆，则

\hat{\mathbf{m}} = \left(\mathbf{G}^{\mathrm T}\mathbf{G}\right)^{-1} \mathbf{G}^{\mathrm T}\mathbf{d}

加权最小二乘

不同数据的误差和量纲通常不同。设数据协方差矩阵为 $\mathbf{C}_d$ ，加权目标函数为

\Phi_d(\mathbf{m}) = \left(\mathbf{G}\mathbf{m}-\mathbf{d}\right)^{\mathrm T} \mathbf{C}_d^{-1} \left(\mathbf{G}\mathbf{m}-\mathbf{d}\right)

其解为

\hat{\mathbf{m}} = \left( \mathbf{G}^{\mathrm T}\mathbf{C}_d^{-1}\mathbf{G} \right)^{-1} \mathbf{G}^{\mathrm T}\mathbf{C}_d^{-1}\mathbf{d}

若各数据误差相互独立，则 $\mathbf{C}_d$ 为对角矩阵。误差越小的数据具有越大的权重。

阻尼反演

当地震数据不足或设计矩阵病态时，最小二乘解可能具有很大的模型振幅。可以加入阻尼项：

\Phi(\mathbf{m}) = \|\mathbf{W}_d(\mathbf{G}\mathbf{m}-\mathbf{d})\|_2^2 +\lambda^2 \|\mathbf{m}-\mathbf{m}_{\mathrm{ref}}\|_2^2

其中 $\mathbf{m}_{\mathrm{ref}}$ 为参考模型， $\lambda$ 为阻尼因子。对应的正规方程为

\left( \mathbf{G}^{\mathrm T}\mathbf{W}_d^{\mathrm T}\mathbf{W}_d\mathbf{G} +\lambda^2\mathbf{I} \right)\mathbf{m} = \mathbf{G}^{\mathrm T}\mathbf{W}_d^{\mathrm T}\mathbf{W}_d\mathbf{d} +\lambda^2\mathbf{m}_{\mathrm{ref}}

阻尼因子越大，反演模型越接近参考模型，模型振幅越小，但数据拟合程度也可能下降。

平滑反演

层析成像中常要求相邻网格之间的模型变化平缓。令 $\mathbf{L}$ 为空间一阶或二阶差分算子，定义

\Phi(\mathbf{m}) = \|\mathbf{W}_d(\mathbf{G}\mathbf{m}-\mathbf{d})\|_2^2 +\lambda_s^2\|\mathbf{L}\mathbf{m}\|_2^2

得到

\left( \mathbf{G}^{\mathrm T}\mathbf{W}_d^{\mathrm T}\mathbf{W}_d\mathbf{G} +\lambda_s^2\mathbf{L}^{\mathrm T}\mathbf{L} \right)\mathbf{m} = \mathbf{G}^{\mathrm T}\mathbf{W}_d^{\mathrm T}\mathbf{W}_d\mathbf{d}

阻尼约束模型的绝对幅值，平滑约束模型的空间梯度或曲率。实际反演中常同时使用阻尼和平滑：

\begin{bmatrix} \mathbf{W}_d\mathbf{G}\\ \lambda_d\mathbf{I}\\ \lambda_s\mathbf{L} \end{bmatrix} \mathbf{m} = \begin{bmatrix} \mathbf{W}_d\mathbf{d}\\ \lambda_d\mathbf{m}_{\mathrm{ref}}\\ \mathbf{0} \end{bmatrix}

正则化参数

正则化参数控制数据拟合与模型复杂度之间的权衡。

正则化过小：模型可能拟合噪声，出现幅值过大或空间振荡；
正则化过大：模型过于平滑或过分接近初始模型，真实异常被削弱；
常利用 L 曲线、交叉验证、数据误差水平或多组参数试验选择正则化参数。

奇异值分解

设计矩阵可以分解为

\mathbf{G} = \mathbf{U}\boldsymbol{\Lambda}\mathbf{V}^{\mathrm T}

其中 $\boldsymbol{\Lambda}$ 的对角元素为奇异值。广义逆为

\mathbf{G}^{\dagger} = \mathbf{V}\boldsymbol{\Lambda}^{\dagger}\mathbf{U}^{\mathrm T}

模型解为

\hat{\mathbf{m}} = \mathbf{G}^{\dagger}\mathbf{d}

较小的奇异值对应数据难以约束的模型组合，直接求逆会放大数据噪声。截断奇异值分解通过舍弃较小奇异值提高反演稳定性，但同时损失模型分辨率。

大规模层析成像通常不显式计算矩阵逆，而采用共轭梯度、LSQR 等迭代方法求解。

非线性反演

对于非线性正演关系

\mathbf{d} = \mathbf{g}(\mathbf{m})

可以从初始模型 $\mathbf{m}_0$ 出发，重复进行正演、线性化和模型更新：

\mathbf{G}_k\delta\mathbf{m}_k = \mathbf{d}^{\mathrm{obs}} -\mathbf{g}(\mathbf{m}_k)

\mathbf{m}_{k+1} = \mathbf{m}_k +\alpha_k\delta\mathbf{m}_k

其中 $\alpha_k$ 为步长。

Gauss-Newton 方法

定义最小二乘目标函数

\Phi(\mathbf{m}) = \frac{1}{2} \|\mathbf{g}(\mathbf{m})-\mathbf{d}\|_2^2

忽略正演函数的二阶导数，可以得到 Gauss-Newton 更新方程

\left( \mathbf{G}_k^{\mathrm T}\mathbf{G}_k +\lambda^2\mathbf{I} \right) \delta\mathbf{m}_k = \mathbf{G}_k^{\mathrm T} \left[ \mathbf{d} -\mathbf{g}(\mathbf{m}_k) \right]

非线性反演结果依赖初始模型。若初始模型距离真实模型较远，线性化可能失效，反演可能收敛到局部极小值。

全局搜索方法

遗传算法、模拟退火、邻域算法和蒙特卡洛方法可以在较大的模型空间中搜索，降低陷入局部极小值的风险。它们适用于参数数量较少的一维结构或震源参数反演，但计算量随模型维数迅速增大。

贝叶斯反演

贝叶斯反演将模型视为随机变量：

p(\mathbf{m}|\mathbf{d}) \propto p(\mathbf{d}|\mathbf{m})p(\mathbf{m})

其中 $p(\mathbf{m})$ 为先验分布， $p(\mathbf{d}|\mathbf{m})$ 为似然函数， $p(\mathbf{m}|\mathbf{d})$ 为后验分布。贝叶斯方法能够同时描述最可能模型和模型不确定性，但高维问题的后验采样计算量很大。

模型分辨率与不确定性

数据拟合良好并不意味着模型唯一或可靠。反演结果需要从射线覆盖、模型分辨率、数据误差和参数权衡等方面进行评价。

模型分辨率矩阵

若线性反演解为

\hat{\mathbf{m}} = \mathbf{G}^{\dagger}\mathbf{d}

且数据由真实模型产生

\mathbf{d} = \mathbf{G}\mathbf{m}_{\mathrm{true}}

则

\hat{\mathbf{m}} = \mathbf{R}_m\mathbf{m}_{\mathrm{true}}

\mathbf{R}_m = \mathbf{G}^{\dagger}\mathbf{G}

其中 $\mathbf{R}_m$ 为模型分辨率矩阵。理想情况下 $\mathbf{R}_m=\mathbf{I}$ 。对角元素反映对应模型参数的恢复程度，非对角元素反映不同位置或参数之间的混叠。

后验协方差

在线性高斯问题中，若模型先验协方差为 $\mathbf{C}_m$ ，则后验协方差可以写为

\mathbf{C}_{m|\mathbf{d}} = \left( \mathbf{G}^{\mathrm T}\mathbf{C}_d^{-1}\mathbf{G} +\mathbf{C}_m^{-1} \right)^{-1}

其对角元素给出模型参数方差，非对角元素反映参数间相关性。

检测板测试

层析成像中常利用检测板测试评估空间恢复能力：

构造正负相间的速度异常模型；
使用实际震源、台站和射线路径计算合成数据；
向合成数据加入适当噪声；
采用与实际数据完全相同的参数化和反演设置；
比较恢复模型与输入模型。

检测板测试只能说明给定异常尺度和反演设置下的恢复能力，不能证明实际反演模型一定正确。异常形状、网格尺度、平滑强度和噪声水平都会影响测试结果。

体波走时层析成像

走时正演问题

在射线理论下，地震波沿射线路径 $L$ 的走时为

t = \int_L s(\mathbf{x})\mathrm dl = \int_L\frac{1}{v(\mathbf{x})}\mathrm dl

其中 $s=1/v$ 为慢度。

在参考模型 $s_0(\mathbf{x})$ 附近，忽略射线路径变化的一阶影响，走时扰动为

\delta t = \int_{L_0}\delta s(\mathbf{x})\mathrm dl

将地下介质离散为多个网格，若第 $i$ 条射线在第 $j$ 个网格内的长度为 $l_{ij}$ ，则

\delta t_i = \sum_jl_{ij}\delta s_j

矩阵形式为

\delta\mathbf{t} = \mathbf{G}\delta\mathbf{s}

其中

G_{ij}=l_{ij}

走时层析成像的基本任务是利用大量具有不同传播路径的走时数据，反演地下介质的慢度或速度分布。

射线追踪

实际地球介质非均匀，射线路径依赖速度模型。常用射线追踪方法包括：

射线射击法；
射线弯曲法；
动力射线追踪；
快速行进法；
高斯束法。

程函方程为

|\nabla T(\mathbf{x})| = s(\mathbf{x}) = \frac{1}{v(\mathbf{x})}

其中 $T(\mathbf{x})$ 为从震源到空间位置 $\mathbf{x}$ 的首波走时。快速行进法可以直接求解整个空间的首波走时场，再由走时梯度追踪射线路径。

由于射线路径随速度模型改变，体波层析通常需要迭代：

在当前速度模型中追踪射线并计算理论走时；
计算观测走时残差；
反演速度扰动；
更新模型并重新追踪射线；
直到走时残差和模型变化达到要求。

局地地震层析成像

局地地震层析成像使用研究区域内部的地震和台站，主要反演地壳和上地幔浅部结构。地震发震时刻和位置未知，走时方程应同时考虑速度结构和震源参数：

\delta t_{ij} = \delta t_i^0 +\frac{\partial T_{ij}}{\partial x_i}\delta x_i +\frac{\partial T_{ij}}{\partial y_i}\delta y_i +\frac{\partial T_{ij}}{\partial z_i}\delta z_i +\int_{L_{ij}}\delta s(\mathbf{x})\mathrm dl

因此震源位置误差和速度结构误差存在权衡。局地层析通常联合反演震源参数和三维速度结构。

P 波和 S 波走时可以分别约束 $V_P$ 和 $V_S$ ，也可以进一步计算

\frac{V_P}{V_S}

高 $V_P/V_S$ 可能与流体、部分熔融或基性组分有关，但解释时需要结合绝对波速和其他地质资料。

双差层析成像

对于空间上相近的两个地震 $i$ 和 $j$ ，它们到同一台站 $k$ 的射线路径在震源区外大体相同。走时差为

\delta t_{k}^{ij} = \left(t_k^i-t_k^j\right)^{\mathrm{obs}} -\left(t_k^i-t_k^j\right)^{\mathrm{cal}}

路径共同部分的模型误差在差分中部分抵消，因此双差数据对地震相对位置更敏感。TomoDD 等方法利用绝对走时和差分走时，同时改善地震相对定位和速度结构。

远震层析成像

远震到达区域台阵时，同一事件到相邻台站的射线在震源端和大部分下地幔中十分接近，而在台阵下方逐渐分开。通常使用相对走时残差：

\delta t_{ij} = t_{ij} -\overline{t_i}

其中 $\overline{t_i}$ 为同一事件在台阵内的平均走时或平均残差。这样可以减弱震源时间、震源区结构和远端传播路径的共同影响，主要反映台阵下方的速度异常。

远震层析适合研究岩石圈和上地幔结构，但对台阵以外和浅层结构的约束较弱，垂向分辨率通常也受有限入射角范围限制。

Pn 波层析成像

Pn 波是在莫霍面附近折射并沿上地幔顶部传播的震相，常使用约 $2^\circ$ 到 $15^\circ$ 的区域地震记录。Pn 走时主要约束：

上地幔顶部 P 波速度；
莫霍面起伏；
Pn 波方位各向异性；
地壳平均慢度和台站延迟。

Pn 层析对上地幔顶部具有较好的水平覆盖，但其结果依赖对地壳厚度和地壳传播时间的处理。

有限频层析成像

射线理论建立在无限高频近似上，将数据敏感性集中在一条无限细射线上。实际地震波具有有限频带，走时对射线周围一定体积内的结构都可能敏感：

\delta t = \iiint_V K(\mathbf{x}) \frac{\delta v(\mathbf{x})}{v(\mathbf{x})} \mathrm dV

其中 $K(\mathbf{x})$ 为有限频走时敏感核。

对通过互相关测量的走时，敏感核通常呈“香蕉—甜甜圈”形态：

敏感区域沿几何射线周围分布；
在几何射线附近可能存在低敏感区域；
敏感区域宽度与第一菲涅耳区尺度相近；
周期越长，敏感范围越宽，空间分辨率越低。

有限频理论能够描述波前愈合等射线理论不能表示的现象，但需要更多正演计算，并依赖参考模型和 Born 近似。

面波层析成像

面波沿地表传播，具有明显的频散特征。不同周期的面波对不同深度范围敏感，因此可以利用相速度和群速度随周期的变化反演剪切波速度结构。

相速度与群速度

相速度定义为

c(\omega) = \frac{\omega}{k}

群速度定义为

U(\omega) = \frac{\mathrm d\omega}{\mathrm dk}

相速度描述单一相位的传播速度，群速度描述波包或能量包络的传播速度。在频散介质中

U = c+k\frac{\mathrm dc}{\mathrm dk} = c-\lambda\frac{\mathrm dc}{\mathrm d\lambda}

群速度受传播路径上不同频率成分的干涉影响，相速度则需要确定相位周数或进行相位解缠绕。

群速度测量

峰谷法

对于具有清晰频散特征的面波记录，可以根据振幅谱中的相邻峰谷和震源距估计群速度。该方法直观，但容易受到噪声、多阶模式和震源谱的影响。

移动时间窗分析

对地震记录进行短时傅里叶变换：

X(\tau,\omega) = \int_{-\infty}^{+\infty} x(t)w(t-\tau)e^{-i\omega t}\mathrm dt

可以得到能量随时间和频率的分布。对每个频率寻找能量最大时刻，再根据传播距离计算群速度。

多重滤波分析

多重滤波分析使用一系列中心频率不同的窄带高斯滤波器处理地震记录，计算各窄带信号的包络：

A(t,\omega_0) = \left| x_{\omega_0}(t) +i\mathcal{H}[x_{\omega_0}(t)] \right|

其中 $\mathcal{H}$ 表示 Hilbert 变换。包络最大值对应频率 $\omega_0$ 的群到时，从而得到群速度频散曲线。频时分析常简称为 FTAN。

相速度测量

设震源到台站距离为 $\Delta$ ，观测相位为

\phi(\omega) = k(\omega)\Delta +\phi_0(\omega)

则

c(\omega) = \frac{\omega\Delta} {\phi(\omega)-\phi_0(\omega)+2\pi n}

其中 $\phi_0$ 是震源和仪器等产生的初始相位， $n$ 为相位周数。单台相速度测量依赖震源相位和正确的相位解缠绕。

双台法

若两个台站大致位于同一大圆路径上，并记录同一地震，则两台之间的相位差为

\Delta\phi(\omega) = k(\omega)\Delta x

从而

c(\omega) = \frac{\omega\Delta x} {\Delta\phi(\omega)+2\pi n}

双台法可以消去共同的震源相位和部分传播路径效应，但要求事件和两个台站近似位于同一大圆路径上。

台阵法

对于线性或二维台阵，可以利用频率—波数变换估计波数矢量：

P(\omega,\mathbf{k}) = \left| \sum_j U_j(\omega) e^{-i\mathbf{k}\cdot\mathbf{x}_j} \right|^2

谱峰对应主要传播方向和波数，由此得到

c(\omega) = \frac{\omega}{|\mathbf{k}|}

$f-k$ 变换适合规则台阵。F-J 变换利用频率—距离域信息，能够提取基阶和高阶面波频散。

二维相速度和群速度成像

对固定周期 $T$ ，第 $i$ 条传播路径的面波走时为

t_i(T) = \int_{L_i} \frac{1}{c(\mathbf{x},T)} \mathrm dl

在参考相速度 $c_0(T)$ 附近，定义慢度扰动

\delta s(\mathbf{x},T) = \frac{1}{c(\mathbf{x},T)} -\frac{1}{c_0(T)}

则走时残差为

\delta t_i(T) = \int_{L_i}\delta s(\mathbf{x},T)\mathrm dl

对不同周期分别反演，可以得到一系列二维相速度或群速度分布。短周期主要反映浅部结构，长周期主要反映深部结构。

天然地震面波层析常使用单台法、双台法或双平面波方法。双平面波方法用两个平面波近似区域内复杂入射波场，可以同时考虑入射方向、振幅和相位变化。

深度敏感核

面波频散对一维介质参数的扰动可以写为

\frac{\delta c(T)}{c(T)} = \int_0^\infty \left[ K_\alpha(z,T) \frac{\delta\alpha(z)}{\alpha(z)} +K_\beta(z,T) \frac{\delta\beta(z)}{\beta(z)} +K_\rho(z,T) \frac{\delta\rho(z)}{\rho(z)} \right]\mathrm dz

其中 $K_\alpha$ 、 $K_\beta$ 和 $K_\rho$ 分别为相速度对 P 波速度、S 波速度和密度的深度敏感核。

面波对 S 波速度最敏感。瑞利波同时对 P、SV 波速度和密度敏感，Love 波主要对 SH 波速度敏感。不同周期和不同阶次的面波具有不同的深度敏感范围。

面波三维成像通常分为两步：

对每个周期反演二维相速度或群速度分布；
在每个水平位置利用局地频散曲线反演一维速度结构。

将所有位置的一维结果组合，可以得到三维剪切波速度模型。

面波频散曲线反演

设观测频散数据为 $\mathbf{d}^{\mathrm{obs}}$ ，一维分层模型为

\mathbf{m} = (V_{S1},V_{S2},\cdots,h_1,h_2,\cdots)^{\mathrm T}

在参考模型附近，

\delta\mathbf{d} = \mathbf{G}\delta\mathbf{m}

速度与界面深度之间通常存在权衡，有限周期范围也限制了可恢复的深度。实际反演常固定或弱约束 $V_P$ 和密度，重点求解 $V_S$ 。

背景噪声面波层析成像

对于两个台站 $A$ 和 $B$ ，连续噪声记录的互相关函数为

C_{AB}(\tau) = \int u_A(t)u_B(t+\tau)\mathrm dt

当噪声源在时间和方位上足够均匀时，长期叠加的互相关函数可以近似恢复两台站之间的经验格林函数：

C_{AB}(\tau) \propto G_{AB}(\tau)+G_{AB}(-\tau)

从互相关函数中提取瑞利波或 Love 波的相速度和群速度，即可进行层析成像。

背景噪声层析的特点包括：

不依赖地震事件，路径主要由台站对决定；
容易获得区域内较均匀的路径覆盖；
常用于地壳和浅部结构研究；
结果会受到噪声源分布、季节变化、台站响应和预处理方法影响。

天然地震面波包含较长周期信息，适合研究地壳和上地幔；背景噪声通常在短周期更稳定，适合研究浅部和地壳结构。二者可以联合使用。

高阶面波

实际记录可能同时包含基阶和多个高阶模式。高阶模式通常对更深结构敏感，可以提高深部速度结构的分辨率。

但不同模式之间可能发生干涉，模式识别错误会导致频散曲线和反演结果偏差。台阵变换、F-J 变换和波数域方法有助于分离不同阶次。

有限频效应

面波也具有有限频敏感区域。二维射线层析将敏感性集中在大圆路径上，而有限频面波敏感核描述相位或振幅对路径周围结构的敏感性。

有限频效应在长周期、短路径或强横向非均匀区域更明显。若数据分辨率不足，使用复杂有限频理论不一定能得到更可靠的结果。

瑞利波椭圆率

瑞利波质点在垂直—径向平面内作椭圆运动。椭圆率常写为

\eta(T) = \frac{|U_H(T)|}{|U_V(T)|}

或使用水平—垂直谱比表示。瑞利波椭圆率对浅层结构和近地表速度梯度较敏感，与相速度联合反演可以改善浅层 $V_S$ 结构的约束。

全波形反演

走时层析只使用波形中的到时信息，全波形反演同时利用地震记录的走时、振幅、相位和波形变化。

目标函数

最常用的波形差目标函数为

\chi(\mathbf{m}) = \frac{1}{2} \sum_r \int \left[ u_r^{\mathrm{syn}}(t,\mathbf{m}) -u_r^{\mathrm{obs}}(t) \right]^2 \mathrm dt

其中 $r$ 表示震源—接收点组合或记录分量。

反演过程为

\mathbf{m}_{k+1} = \mathbf{m}_k -\alpha_k\mathbf{p}_k

其中 $\mathbf{p}_k$ 是由目标函数梯度和优化算法确定的搜索方向。

伴随状态法

若直接计算每个模型参数对每条波形的偏导数，计算量极大。伴随状态法利用一次正演和一次伴随传播计算一个地震事件对应的目标函数梯度。

基本流程为：

在当前模型中进行三维地震波正演，得到合成波形；
计算合成波形与观测波形之间的残差；
将时间反转后的残差作为伴随源，从接收点反向传播；
将正向波场和伴随波场在时间和空间上相关，得到敏感核或梯度；
根据梯度更新模型；
重复迭代直到目标函数收敛。

正演可以采用有限差分法、有限元法、谱元法等数值方法。

周期跳跃

若观测波形与合成波形的相位差超过约半个周期，波形差目标函数可能将错误的波峰相互匹配，称为周期跳跃。

减弱周期跳跃的常用方法包括：

使用较准确的初始模型；
从长周期到短周期逐步反演；
先反演走时或包络，再反演完整波形；
使用互相关走时、相位或最优输运等目标函数；
分阶段加入不同震相和数据分量。

全波形反演理论分辨率高，但对震源参数、仪器响应、衰减、各向异性和三维数值模拟误差敏感。

接收函数

基本原理

远震 P 波近垂直入射到台站下方，在速度间断面上产生 P 到 S 转换波。接收函数通过去除震源时间函数和远端传播路径的共同影响，突出接收点下方的局地结构响应。

三分量记录通常旋转到垂直、径向和切向分量。对于 P 波接收函数，

Z(\omega) = S(\omega)I(\omega)E_Z(\omega)

R(\omega) = S(\omega)I(\omega)E_R(\omega)

其中 $S(\omega)$ 为震源和共同传播路径效应， $I(\omega)$ 为仪器响应， $E_Z$ 和 $E_R$ 为台站下方结构在垂直和径向分量上的响应。

利用径向分量对垂直分量反褶积：

H_R(\omega) = \frac{R(\omega)}{Z(\omega)} \approx \frac{E_R(\omega)}{E_Z(\omega)}

得到径向接收函数。切向接收函数可用于识别倾斜界面和各向异性。

水位法反褶积

直接频域相除在 $|Z(\omega)|$ 很小时会放大噪声。水位法将接收函数写为

H_R(\omega) = \frac{R(\omega)Z^*(\omega)} {\max\left[ Z(\omega)Z^*(\omega), c\max|Z(\omega)|^2 \right]} G(\omega)

其中 $c$ 为水位系数， $G(\omega)$ 为高斯低通滤波器。

水位系数过小容易放大谱谷处的噪声，过大则会损失波形细节。高斯参数控制接收函数的有效频带和时间分辨率。

时间域迭代反褶积

时间域迭代反褶积用一系列脉冲近似接收函数，使垂直分量与接收函数的褶积逐步逼近径向分量：

R(t) \approx Z(t)*H_R(t)

基本过程为：

计算垂直分量与径向分量残差的互相关；
找到最大互相关及其时间延迟；
向接收函数中加入相应时刻和振幅的脉冲；
将该脉冲与垂直分量褶积并更新残差；
重复迭代，直到残差下降到阈值。

迭代反褶积能够限制脉冲个数，通常比直接谱相除稳定，但结果仍依赖滤波、时间窗和停止条件。

P-S 转换震相

水平层状介质中，主要接收函数震相包括：

$Ps$ ：间断面上的 P 到 S 转换波；
$PpPs$ ：自由表面一次反射后产生的转换波；
$PpSs+PsPs$ ：包含两次 S 波传播的多次波。

设地壳厚度为 $H$ ，P 波射线参数为 $p$ ，则

t_{Ps} = H \left[ \sqrt{\frac{1}{V_S^2}-p^2} -\sqrt{\frac{1}{V_P^2}-p^2} \right]

t_{PpPs} = H \left[ \sqrt{\frac{1}{V_S^2}-p^2} +\sqrt{\frac{1}{V_P^2}-p^2} \right]

t_{PpSs+PsPs} = 2H \sqrt{\frac{1}{V_S^2}-p^2}

转换波到时主要约束界面深度，振幅和极性还与界面两侧速度对比、频率和入射角有关。

接收函数波形反演

给定一维层状速度模型，可以通过反射率法等方法计算理论接收函数。接收函数反演通过拟合观测与理论波形求解速度结构。

线性化后

\delta\mathbf{r} = \mathbf{G}\delta\mathbf{m}

其中 $\mathbf{r}$ 为接收函数采样值，模型可以包括各层 $V_S$ 、 $V_P/V_S$ 和层厚。

接收函数反演存在明显非唯一性：

界面深度和层内速度之间存在权衡；
速度突变与速度梯度可能产生相似波形；
不同薄层结构可能在有限频带内无法区分；
高频接收函数分辨率高但更容易受噪声和散射影响。

接收函数常与面波频散联合反演。接收函数主要约束速度间断面，面波频散主要约束绝对剪切波速度，二者具有较好的互补性。

H-κ 扫描

H-κ 扫描利用 $Ps$ 和多次波的理论到时，同时估计地壳厚度 $H$ 和波速比

\kappa = \frac{V_P}{V_S}

对每一组候选 $H$ 和 $\kappa$ 计算理论到时，并构造叠加函数

S(H,\kappa) = w_1r(t_{Ps}) +w_2r(t_{PpPs}) -w_3r(t_{PpSs+PsPs})

其中 $r(t)$ 为径向接收函数， $w_1,w_2,w_3$ 为权重。第三项前的负号来自该多次波通常与 $Ps$ 极性相反。

叠加函数最大值对应最优 $H$ 和 $\kappa$ 。该方法假设地壳近似水平层状，并需要给定平均 $V_P$ 。沉积层、倾斜莫霍面、各向异性和复杂壳内界面都可能造成偏差。

CCP 叠加

共转换点叠加将不同事件的接收函数沿理论射线路径进行时间—深度转换，再把转换点落入同一空间网格的振幅叠加。

基本步骤为：

选取速度模型并计算 P 波和转换 S 波射线路径；
将接收函数时间轴转换为深度；
确定每个时间采样对应的转换点位置；
将不同事件和台站的接收函数投影到空间网格；
对同一网格内的振幅进行叠加。

CCP 叠加能够成像莫霍面、岩石圈底界面和 410、660 km 间断面。它依赖速度模型，并受到射线弯曲、三维结构和有限频效应影响。

更严格的成像方法包括 Kirchhoff 偏移、波动方程偏移和 Born 逆散射成像。

沉积层影响

低速沉积层会产生强烈的自由表面多次波，使接收函数出现长周期振荡，并掩盖深部转换震相。常用处理方法包括：

在反演中显式加入沉积层；
利用谐振频率和多次波校正；
与面波或 H/V 数据联合约束浅层结构；
谨慎选择滤波频带和高斯参数。

倾斜界面和各向异性

水平各向同性介质中，P-S 转换能量主要出现在径向分量，切向分量应较弱。明显的切向接收函数可能来自倾斜界面或各向异性。

倾斜界面的典型特征为：

径向和切向转换波到时基本一致；
相反后方位角的切向分量极性反转；
转换波到时和振幅随方位呈周期变化。

各向异性的典型特征为：

径向和切向波形之间可能存在时间差；
切向分量可近似表现为径向分量的时间导数；
相反后方位角的切向分量不一定发生极性反转。

实际资料中二者可能同时存在，需要利用不同方位和不同射线参数的接收函数联合判断。

接收函数的应用

接收函数主要用于研究：

地壳厚度和莫霍面起伏；
地壳平均 $V_P/V_S$ ；
壳内低速层和沉积盆地；
岩石圈—软流圈边界和地幔岩石圈间断面；
410 km 和 660 km 间断面；
倾斜界面与地壳、上地幔各向异性。

衰减层析成像

实际地球介质具有非弹性，地震波振幅除几何扩散外还会随传播距离衰减。对于频率 $f$ 的简谐波，可近似写为

A(f,r) = A_0(f)G(r) \exp\left[ -\frac{\pi fr}{Qv} \right]

其中 $G(r)$ 为几何扩散项， $Q$ 为品质因数。

沿非均匀介质传播时

\ln\frac{A}{A_0G} = -\pi f \int_L \frac{1}{Q(\mathbf{x})v(\mathbf{x})} \mathrm dl

因此可以利用振幅、谱衰减、尾波或波形包络反演 $Q_P$ 和 $Q_S$ 的空间分布。

低 $Q$ 表示衰减强，可能与高温、部分熔融、流体或强散射有关；高 $Q$ 通常对应冷而稳定的岩石圈。振幅同时受到震源辐射、场地效应、散射、聚焦和离焦影响，因此衰减反演通常比走时反演更容易受到系统误差影响。

各向异性反演

各向异性介质中的地震波速度随传播方向和偏振方向变化。地壳各向异性常与裂隙、层理和形状优选排列有关，上地幔各向异性常与橄榄石晶格优选定向和地幔流动有关。

SKS 分裂

近垂直传播的 SKS 波在各向异性介质中分裂为快、慢两个近似正交的偏振分量。常用参数为：

快波偏振方向 $\phi$ ；
快慢波时间延迟 $\delta t$ 。

通过旋转水平分量并校正时间延迟，使校正后的横向能量最小，可以估计 $\phi$ 和 $\delta t$ 。

单层各向异性中， $\phi$ 反映水平投影的快轴方向， $\delta t$ 近似与各向异性强度和层厚的乘积有关。多层各向异性会产生随事件后方位角变化的表观分裂参数。

面波方位各向异性

弱方位各向异性介质中的面波相速度可以近似展开为

c(T,\theta) = c_0(T) +A_2(T)\cos2\theta +B_2(T)\sin2\theta +A_4(T)\cos4\theta +B_4(T)\sin4\theta

其中 $\theta$ 为传播方位角。二阶项通常与水平快轴方向有关，四阶项反映更复杂的各向异性结构。

不同周期的面波对不同深度敏感，因此可以反演各向异性随深度的变化。可靠的方位各向异性成像要求良好的方位覆盖，否则各向异性可能与各向同性速度异常或波场效应混叠。

联合反演

不同类型的地震数据对模型参数具有互补敏感性。例如：

接收函数对速度间断面敏感，但对绝对速度约束较弱；
面波频散对绝对 $V_S$ 敏感，但速度与界面深度存在权衡；
体波走时对沿射线路径的速度积分敏感；
重力资料对密度结构敏感；
地震波形同时受速度、衰减、密度和震源参数影响。

联合反演可以构造统一目标函数：

\Phi(\mathbf{m}) = \sum_k w_k \left\| \mathbf{W}_k \left[ \mathbf{g}_k(\mathbf{m}) -\mathbf{d}_k \right] \right\|_2^2 +\Phi_m(\mathbf{m})

其中 $k$ 表示不同数据类型， $w_k$ 为数据类型之间的相对权重， $\Phi_m$ 为模型正则化项。

联合反演并不是简单地把不同数据放在一起。需要处理：

不同数据的量纲和误差水平；
不同数据量造成的权重失衡；
不同数据对应的空间尺度和敏感深度；
不同物理参数之间的岩石物理关系；
一种数据的系统误差向其他模型参数传播。

地震反演的一般流程

数据选择与处理

根据研究目标选择走时、频散、波形、振幅或接收函数等数据；
检查仪器响应、时钟误差、信噪比和异常记录；
统一滤波频带、采样率和测量标准；
估计数据误差并设置权重；
检查震源、台站和传播路径的空间覆盖。

模型参数化

模型可以采用：

规则或不规则网格；
分层模型；
节点插值；
块状参数化；
球谐函数；
样条或其他基函数。

参数化尺度不应明显小于数据能够分辨的尺度。网格过细不会自动提高真实分辨率，反而会增加欠定性和模型参数相关性。

初始模型与正演

初始模型应能够解释数据的主要特征，并尽量避免使非线性反演陷入局部极小值。正演方法必须与数据类型和目标频带相适应：

走时数据可以使用射线追踪或程函方程；
有限频数据需要计算敏感核；
面波频散需要求解本征值和本征函数；
全波形反演需要进行数值波场模拟；
接收函数需要计算层状或三维介质响应。

目标函数与正则化

目标函数应同时包含数据拟合和模型约束。需要测试：

阻尼和平滑强度；
各数据类型权重；
模型边界条件；
参考模型；
异常值处理和稳健损失函数。

迭代与收敛

每次迭代需要检查：

数据残差是否下降；
模型更新量是否减小；
模型是否出现不合理振荡或边界异常；
不同数据类型是否得到同时改善；
结果是否依赖少数事件或台站。

目标函数不再明显下降不一定表示获得真实模型，也可能表示反演进入局部极小值或正则化限制了进一步更新。

结果评价

最终结果应同时给出：

数据拟合程度和残差分布；
射线或敏感核覆盖；
检测板或其他合成恢复测试；
不同正则化参数下的结果；
初始模型依赖性；
模型不确定性或可信区间；
与其他地球物理、地质和岩石物理资料的比较。

反演图像中的每一个异常都同时受到真实结构、数据覆盖、正则化和模型参数化的影响。因此，对反演结果的解释必须以分辨率和不确定性分析为基础。

版权所有

版权归属：Staaaaaaaaar

许可证：署名-非商业性-相同方式共享 4.0 国际 (CC-BY-NC-SA-4.0)