非均匀介质中的高频近似解

在非均匀介质中，弹性波的传播受到材料属性随空间位置变化的影响。当波的频率很高（即波长远小于材料属性变化的特征长度）时，我们可以使用 WKB 近似 (Wentzel-Kramers-Brillouin approximation) 或射线理论 (Ray Theory) 来推导其高频近似解。这种方法将复杂的波动方程转化为求解相位的程函方程 (Eikonal Equation) 和求解振幅的输运方程 (Transport Equation)。

一维非均匀介质

考虑一根截面积恒定的一维非均匀杆。设 $x$ 为空间坐标， $t$ 为时间， $u(x, t)$ 为纵向位移。杆的材料属性随 $x$ 变化：线密度 $\rho(x)$ 和弹性模量 $E(x)$。根据牛顿第二定律和胡克定律，无源情况下的波动方程为

$$\rho(x) \frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial x} \left( E(x) \frac{\partial u}{\partial x} \right)$$

假设解具有简谐时间依赖性 $u(x, t) = U(x) e^{i \omega t}$，其中 $U(x)$ 是复振幅，$\omega$ 是角频率。将此解代入波动方程，得到

$$E(x) U'' + E'(x) U' + \omega^2 \rho(x) U = 0$$

在高频极限下，波的相位变化剧烈，而振幅变化缓慢。引入 WKB 假设：

$$U(x) = A(x) e^{i \omega \tau(x)}$$

其中 $A(x)$ 是振幅，$\tau(x)$ 是相位函数（程函）。将 $U(x)$ 的表达式代入方程，并收集不同阶数的 $\omega$ 项：

$$\omega^2 \left[\rho A - E (\tau')^2 A\right] + i \omega \left[E (2 \tau'A'+\tau''A) + E'\tau'A \right] + \left[E A'' + E' A' \right] = 0$$

---

令 $\omega^2$ 的系数为零：

$$\rho(x)-E(x) (\tau'(x))^2 = 0 \implies \tau'(x) =\pm \sqrt{\frac{\rho(x)}{E(x)}}= \pm \frac{1}{c(x)}$$

其中 $c(x) = \sqrt{E(x)/\rho(x)}$ 是局部波速。这个方程被称为程函方程，描述了波的相位如何随空间变化。

$$\tau(x) = \pm \int_{x_0}^x \frac{1}{c(\xi)} d\xi$$

---

令 $\omega$ 的系数为零：

$$2E \tau' A' + E \tau'' A + E' \tau' A = 0$$

注意到 $(E \tau')' = E' \tau' + E \tau''$，因此上式可以重写为：

$$2E\tau' A' + (E \tau')' A = 0 \implies \frac{d}{dx} (E \tau' A^2) = 0$$

这个方程被称为输运方程，描述了波的振幅如何随空间变化。因此，$E \tau' A^2 = C$。结合 $\tau' = \pm 1/c$，得到振幅的表达式：

$$A(x) = \frac{C}{\sqrt[4]{E(x)\rho(x)}}=\frac{C}{\sqrt{Z(x)}}$$

其中 $Z(x) = \sqrt{E(x)\rho(x)}$ 是局部波阻抗。因此，高频近似解的形式为

$$u(x,t) \approx \frac{1}{\sqrt[4]{E(x)\rho(x)}} \left[f\left(t - \int_{x_0}^x \frac{d\xi}{c(\xi)}\right)+g\left(t - \int_{x_0}^x \frac{d\xi}{c(\xi)}\right)\right]$$

三维非均匀各向同性介质

在三维非均匀各向同性介质中，材料的密度 $\rho(x)$ 以及拉梅常数 $\lambda(x), \mu(x)$ 均随空间变化，纳维尔方程为

$$\rho \ddot{\mathbf{u}} = (\lambda + \mu) \nabla (\nabla \cdot \mathbf{u}) + \mu \nabla^2 \mathbf{u} + \nabla \lambda(\nabla \cdot \mathbf{u}) + \nabla \mu \cdot [\nabla \mathbf{u} + (\nabla \mu)^T ]$$

引入三维射线拟设：

$$\mathbf{u}(\mathbf{r}, t) = \mathbf{A}(\mathbf{r}) e^{i \omega (\tau(\mathbf{r})-t)}$$

将其代入纳维尔方程，并收集不同阶数的 $\omega$ 项。

---

令 $\omega^2$ 的系数为零，得到 Christoffel 方程：

$$(\lambda + \mu)(\mathbf{A}\cdot\nabla \tau) \nabla \tau + \mu |\nabla \tau|^2 \mathbf{A} - \rho \mathbf{A} = 0$$

将 $\mathbf{A}$ 分解为平行于 $\nabla \tau$ 和垂直于 $\nabla \tau$ 的方向：

$$\mathbf{A} = \mathbf{A}_{\parallel} + \mathbf{A}_{\perp}$$

代入上面方程，经过一些具有技巧性的代数运算后，得到

$$\left((\lambda + 2\mu) |\nabla \tau|^2 - \rho^2\right)\mathbf{A}_{\parallel} + \left(\lambda |\nabla \tau|^2 - \rho^2\right)\mathbf{A}_{\perp} = 0$$

为了使 $\mathbf{A}$ 有非零解，有两种情况，分别对应于 P 波和 S 波：

$$\mathbf{A}_{\parallel} \ne 0, \mathbf{A}_{\perp} = 0 \implies (\lambda + 2\mu) |\nabla \tau|^2 - \rho^2 = 0 \implies |\nabla \tau| = \frac{1}{\alpha} \\ \mathbf{A}_{\parallel} = 0, \mathbf{A}_{\perp} \ne 0 \implies \lambda |\nabla \tau|^2 - \rho^2 = 0 \implies |\nabla \tau| = \frac{1}{\beta}$$

---

令 $\omega$ 的系数为零，得到

$$(\lambda + \mu)[\nabla(\mathbf{A}\cdot \nabla \tau)+(\nabla\cdot\mathbf{A})\nabla\tau] + \mu [2(\nabla\tau\cdot \nabla)\mathbf{A} + (\nabla^2\tau)\mathbf{A}] + \\ \nabla \lambda (\mathbf{A}\cdot\nabla\tau) + (\nabla\mu\cdot\mathbf{A})\nabla\tau + (\nabla\mu\cdot\nabla\tau)\mathbf{A} = 0$$

下面分别以 P 波和 S 波的情况，求解运输方程。

在上面等式两边点乘 $\mathbf{A}$，经过系列恒等变换，得到

$$\nabla\cdot[(\lambda+\mu)(\mathbf{A}\cdot\nabla\tau)\mathbf{A}+\mu |\mathbf{A}|^2 \nabla\tau] = 0$$

对于 P 波，设其标量振幅为 $A_P$，定义单位波向矢量 $\mathbf{n}=\alpha \nabla \tau$，由于 $\mathbf{A}$ 与 $\nabla \tau$ 平行，则复振幅矢量可以表示为

$$\mathbf{A} = A_P \alpha \nabla \tau$$

将其代入上式，化简得到

$$\nabla \cdot (\rho\alpha^2A_P^2\nabla\tau) = 0$$

注

矢量 $\rho\alpha^2A_P^2\nabla\tau$ 表示 P 波的能流密度，此散度方程表明高频近似下沿射线束的能量通量在稳态下是散度为零（守恒）的。

弹性波射线管

应用高斯散度定理于一个由一束相邻射线包围而成的射线管。由于管壁的法向量垂直于 $\nabla \tau$，能量只能从管的两端流出和流入。设两端的截面积为 $J(\mathbf{r}_0)$ 和 $J(\mathbf{r})$，能量守恒要求

$$\int_{J(\mathbf{r}_0)} \rho\alpha^2A_P^2\nabla\tau \cdot d\mathbf{S} = \int_{J(\mathbf{r})} \rho\alpha^2A_P^2\nabla\tau \cdot d\mathbf{S}$$

由于 $\nabla \tau$ 的方向与射线束的传播方向一致，$\nabla\tau\cdot\mathbf{n} = |\nabla\tau|=1/\alpha$。因此上式可以重写为

$$(\rho\alpha A_P^2J)|_{\mathbf{r}_0} = (\rho\alpha A_P^2J)|_{\mathbf{r}}$$

因此，P 波的振幅的解析解为

$$A_P(\mathbf{r}) = A_P(\mathbf{r}_0) \sqrt{\frac{\rho(\mathbf{r}_0)\alpha(\mathbf{r}_0) J(\mathbf{r}_0)}{\rho(\mathbf{r})\alpha(\mathbf{r}) J(\mathbf{r})}}$$

可见振幅的衰减或放大由两个因素控制：

• 三维特有的几何发散效应（射线管截面积 $J(\mathbf{r})$ 的改变）。
• 材料非均匀性带来的波阻抗效应（$\sqrt{\rho\alpha}$ 的改变）。

对于 S 波，设其标量振幅为 $A_S$，定义单位波向矢量 $\mathbf{n}=\beta \nabla \tau$，由于 $\mathbf{A}$ 与 $\nabla \tau$ 垂直，有

$$\mathbf{A} \cdot \nabla \tau = 0$$

代入之前的散度表达式，得到

$$\nabla \cdot (\rho\beta^2A_S^2\nabla\tau) = 0$$

射线方程

我们知道，在高频近似下，$\nabla \tau$ 描述了波相位的空间变化率。程函方程只给出了梯度的模，没有给出方向。我们可以从程函方程出发，推导出描述射线路径的微分方程，称为射线方程。

已知

$$|\nabla \tau|^2 = \frac{1}{c^2(\mathbf{r})}=n^2(\mathbf{r})$$

其中 $c(\mathbf{r})$ 是介质波速，$n(\mathbf{r}) = 1/c(\mathbf{r})$ 是介质慢度（折射率）。

定义射线的路径函数为 $\mathbf{x}(s)$，其中 $s$ 是沿射线的弧长参数。根据几何光学，围绕射线单位切矢量可以得到等式：

$$\frac{d\mathbf{x}}{ds}=\frac{\nabla \tau}{|\nabla \tau|} = \frac{\nabla \tau}{n(\mathbf{r})}$$

考察 $\nabla \tau$ 沿着射线路径 $s$ 的变化率，根据方向导数的定义，有

$$\frac{d}{ds} \nabla \tau = \left( \frac{d\mathbf{x}}{ds} \cdot \nabla \right) \nabla \tau$$

将 $\frac{d\mathbf{x}}{ds} = \frac{\nabla \tau}{n}$ 代入，得到

$$\frac{d}{ds}(\nabla \tau) = \frac{1}{n} (\nabla \tau \cdot \nabla) \nabla \tau$$

利用矢量恒等式

$$(\nabla \tau \cdot \nabla) \nabla \tau = \frac{1}{2} \nabla (|\nabla \tau|^2) = \frac{1}{2} \nabla (n^2) = n \nabla n$$

代入上式，得到

$$\frac{d}{ds} \left( n \frac{d\mathbf{x}}{ds} \right) = \nabla n$$

也即

$$\frac{d}{ds} \left( \frac{1}{c} \frac{d\mathbf{x}}{ds} \right) = \nabla \left( \frac{1}{c} \right)$$

这就是射线方程，描述了波在非均匀介质中的传播路径。它表明，射线的曲率由介质慢度的空间梯度决定。当介质均匀时，射线方程退化，射线为直线。

斯涅耳定律

我们假设介质是分层介质，慢度场可以表示为 $n(\mathbf{r}) = n(z)$。

$$\nabla n = \left( 0, 0, \frac{dn}{dz} \right)$$

射线轨迹为 $\mathbf{x}(s) = [x(s), y(s), z(s)]^T$，将其代入射线方程，得到

$$\dfrac{d}{ds} \left( n \dfrac{dx}{ds} \right) = 0$$

$$\dfrac{d}{ds} \left( n \dfrac{dy}{ds} \right) = 0$$

$$\dfrac{d}{ds} \left( n \dfrac{dz}{ds} \right) = \dfrac{dn}{dz}$$

于是我们得到两个常量：

$$n \frac{dx}{ds} = C_x, \quad n \frac{dy}{ds} = C_y$$

定义 $\theta$ 为射线轨迹与 $z$ 轴（即介质变化方向的法线）之间的夹角，则

$$n\sin \theta = n \sqrt{\left(\frac{dx}{ds}\right)^2 + \left(\frac{dy}{ds}\right)^2} = const$$

对于分层介质中的任意两点（或两个界面），上述等式依然成立，即：

$$n_1 \sin \theta_1 = n_2 \sin \theta_2$$

这就是斯涅耳定律。

射线曲率

曲率 $\kappa$ 定义为射线单位切矢量 $\mathbf{u} = d\mathbf{x}/ds$ 随弧长 $s$ 的变化率。

$$\kappa = \left| \frac{d\mathbf{u}}{ds} \right| = \left| \frac{d^2\mathbf{x}}{ds^2} \right|$$

从射线方程出发：

$$\frac{d}{ds} \left( n \frac{d\mathbf{x}}{ds} \right) = \nabla n$$

展开左侧：

$$n \frac{d^2\mathbf{x}}{ds^2} + \frac{dn}{ds} \frac{d\mathbf{x}}{ds} = \nabla n$$

改写为：

$$n \frac{d^2\mathbf{x}}{ds^2} = \nabla n - \frac{dn}{ds} \mathbf{u}$$

再将

$$\frac{dn}{ds} = \nabla n \cdot \frac{d\mathbf{x}}{ds} = \nabla n \cdot \mathbf{u}$$

代入上式：

$$n \frac{d^2\mathbf{x}}{ds^2} = \nabla n - (\nabla n \cdot \mathbf{u}) \mathbf{u}$$

因此：

$$\kappa = \left| \frac{d^2\mathbf{x}}{ds^2} \right| = \frac{1}{n} | \nabla n - (\nabla n \cdot \mathbf{u}) \mathbf{u} | = \frac{|\nabla n|}{n} \sin \phi = \frac{|\nabla c|}{c} \sin \phi$$

其中 $\phi$ 是慢度梯度 $\nabla n$ 与射线切方向 $\mathbf{u}$ 之间的夹角。

水平层状介质中的走时计算

设介质由 $N$ 个水平均匀层组成，第 $k$ 层参数为：速度 $v_k$（层内为常数）、厚度 $h_k$，速度满足 $v_1 < v_2 < \cdots < v_N$。

由斯涅尔定律可知

$$p=\frac{\sin i_1}{v_1}=\frac{\sin i_2}{v_2}=\cdots = \frac{\sin i_N}{v_N}=const$$

$p$ 称为射线参数，全路径保持不变。各层入射角由 $p$ 确定：

$$\sin i_k = p v_k, \quad \cos i_k = \sqrt{1 - p^2 v_k^2}$$

直达波

直达波沿地表直线传播：

$$T_{direct} = \frac{X}{v_1}$$

反射波

由几何关系不难得到

$$T_{reflected} = \frac{2}{v_1} \sqrt{h_1^2 + \left(\frac{X}{2}\right)^2}$$

两边平方得到

$$T_{reflected}^2 = T_0^2 + \frac{X^2}{v_1^2}$$

其中 $T_0 = 2h_1/v_1$ 是零距反射走时。

提示

上面讨论的仅是单层情形，对于多层情形，反射波的走时方程会更复杂，对于小偏移距，反射波走时可以由 NMO 近似表示：

$$T_{reflected}^2 \approx T_0^2 + \frac{X^2}{v_{rms}^2}$$

其中 $v_{rms}$ 是层间速度的均方根，定义为

$$v_{rms}^2 = \frac{\sum_{k=1}^N v_k^2 \Delta t_k}{\sum_{k=1}^N \Delta t_k}$$

$\Delta t_k=2h_k/v_k$ 是第 $k$ 层的单层双程走时。

折射波

当第一层入射角增大到某个临界值 $\sin i_c = v_1/v_2,p=1/v_2$ 时，第二层入射角 $i_2$ 达到 $90^\circ$，此时射线沿第一层和第二层的界面传播，形成临界折射波。其走时方程为

$$T_{refracted} = \frac{X}{v_2} + \frac{2h_1\cos i_c}{v_1}$$

提示

对于多层情形，设头波在第 $N$ 层上方界面传播，$p=1/v_N$。各层的入射角 $i_k=\arcsin(p v_k)$，走时为

$$T_{N} = \frac{X}{v_N} + 2\sum_{k=1}^{N-1} \frac{h_k \cos i_k}{v_k}$$

转折点

当射线向深处传播时，若速度随深度增大，则 $\sin i$ 逐渐增大。当满足：

$$pv(z_t) = 1 \implies i(z_t) = 90^\circ$$

时，射线水平传播，随后折返地表。深度 $z_t$ 称为转折点。

走时参数方程

考虑地表震源，记 $X$ 表示震中距，$T$ 表示总走时。对于穿过 $N$ 层的射线：

$$X(p)=2\sum_{k=1}^N h_k \tan i_k = 2\sum_{k=1}^N \frac{h_k p v_k}{\sqrt{1-p^2 v_k^2}}$$

$$T(p)=2\sum_{k=1}^N \frac{h_k}{v_k \cos i_k} = 2\sum_{k=1}^N \frac{h_k}{v_k \sqrt{1-p^2 v_k^2}}$$

定义截距时间 $\tau(p) = T(p) - p X(p)$。

$$\tau(p) = 2\sum_{k=1}^N \frac{h_k\sqrt{1-p^2 v_k^2}}{v_k}$$

将离散模型推广到连续模型，我们可以写出积分形式的走时参数方程：

$$X(p) = 2\int_0^{z_t} \frac{p v(z)}{\sqrt{1-p^2 v^2(z)}} dz$$

$$T(p) = 2\int_0^{z_t} \frac{1}{v(z) \sqrt{1-p^2 v^2(z)}} dz$$

$$\tau(p) = 2\int_0^{z_t} \frac{\sqrt{1-p^2 v^2(z)}}{v(z)} dz$$

可以推导出下面关系：

$$\frac{dT}{dX} = p$$

$$\frac{d\tau}{dp} = -X$$

这说明在 $T-X$ 曲线上任取一点 $A = (X_0, T_0)$，该点处切线的斜率即为对应的射线参数 $p$，切线在 $T$ 轴上的截距即为 $\tau_0 = T_0 − pX_0$；在 $\tau-p$ 曲线上任取一点 $B = (p_0, \tau_0)$，该点处切线的斜率为 $-X_0$。

走时曲线的几何形态

线性介质

高速带

当射线进入一个速度比上层更快的高速夹层时，走时曲线会打折并出现三重线。此时走时就不再是震中距的单值函数，向前和向后的过渡点为走时曲线的尖点 (cusp)，称为焦散点 （caustic）。

低速带

当射线进入一个速度比上层更慢的低速夹层时，走时曲线会产生间隙，即在一定距离范围内没有射线到达，这个区域被称为影区 (shadow zone)。

球对称介质中的走时计算

球对称介质中的斯涅耳定律

在球坐标系下，由于模型球对称，介质速度 $c=c(r)$，慢度的梯度也只有径向分量：

$$\nabla n = \frac{dn}{dr}\mathbf{\hat{e}}_r$$

射线方程变为

$$\frac{d}{ds} \left( n \frac{d\mathbf{x}}{ds} \right) = \frac{dn}{dr}\mathbf{\hat{e}}_r$$

两边同时叉乘射线路径矢量 $\mathbf{x}$，得到

$$\mathbf{x} \times \frac{d}{ds} \left( n \frac{d\mathbf{x}}{ds} \right) = 0$$

经变换后得到

$$\frac{d}{ds} \left( \mathbf{x} \times n\frac{d\mathbf{x}}{ds} \right) = 0$$

即

$$\mathbf{x} \times n\frac{d\mathbf{x}}{ds} = const$$

设 $\mathbf{r}$ 与 $d\mathbf{x}/ds$ 之间的夹角为 $i$，因此有

$$r n \sin i =\frac{r \sin i}{c}= p_{sph} = const$$

走时参数方程

由前面推导的几何关系

$$\frac{d\theta}{dr} = \frac{1}{r} \tan i = \frac{1}{r} \frac{\sin i}{\sqrt{1 - \sin^2 i}}$$

代入 $p = \frac{r \sin i}{v}$，即 $\sin i = \frac{pv}{r}$：

$$\frac{d\theta}{dr} = \frac{pv}{r^2 \sqrt{1 - (\frac{pv}{r})^2}}$$

对整条射线积分（从地表 $R$ 到转折点 $r_t$ 再回到地表），总转角 $\Delta$ 为：

$$\Delta(p) = 2 \int_{r_t}^{R} \frac{p v(r)}{r \sqrt{r^2 - p^2 v^2(r)}} dr$$

同理，利用 $dt = \frac{ds}{v} = \frac{dr}{v \cos i}$：

$$\frac{dt}{dr} = \frac{1}{v \sqrt{1 - \sin^2 i}} = \frac{1}{v \sqrt{1 - (\frac{pv}{r})^2}} = \frac{r}{v \sqrt{r^2 - p^2 v^2}}$$

积分得到总走时：

$$T(p) = 2 \int_{r_t}^{R} \frac{r}{v(r) \sqrt{r^2 - p^2 v^2(r)}} dr$$

在地震学中，为了简化分析，常引入截距 $\tau = T - p\Delta$。将上述 $T$ 和 $\Delta$ 的方程代入：

$$\tau(p) = 2 \int_{r_t}^{R} \left( \frac{r}{v \sqrt{r^2 - p^2 v^2}} - \frac{p^2 v}{r \sqrt{r^2 - p^2 v^2}} \right) dr$$

合并同类项后得到极其简洁的形式：

$$\tau(p) = 2 \int_{r_t}^{R} \frac{\sqrt{r^2 - p^2 v^2(r)}}{r v(r)} dr$$

地球展平变换

地球平面展开 (Earth Flattening Transformation, EFT) 通过一组保角映射，将弯曲的球面射线路径“拉直”为平面层状介质中的射线。展平变换的条件是：

• 在同样的震中距条件下射线的走时不变
• 变换后的射线几何保角

根据保角条件，有

$$\frac{p_{sph}c(r)}{r} = p c(z)$$

$$p_{sph}=\frac{rpc(z)}{c(r)}=\frac{rpn(z)}{n(z)}$$

将其带入球对称介质的走时参数方程得到

$$T(p_{sph}) = 2 \int_{r_t}^{R} \frac{n(z)n(r)}{\sqrt{n^2(z)-p^2}} dr$$

$$\Delta(p_{sph}) = 2p \int_{r_t}^{R} \frac{1}{\sqrt{n^2(z)-p^2}} \frac{dr}{r}$$

设存在变换 $r=r(z)$，那么上式变为

$$T(p_{sph}) = -2 \int_{0}^{z_p} \frac{n(z)n(r)}{\sqrt{n^2(z)-p^2}} \frac{dr}{dz} dz$$

$$\Delta(p_{sph}) = -2p \int_{r_t}^{R} \frac{1}{r\sqrt{n^2(z)-p^2}} \frac{dr}{dz} dz$$

根据走时不变性有等式

$$R\Delta(p_{sph}) = X(p)$$

$$T(p_{sph}) = T(p)$$

可以由此推导出

$$z=-R\ln\frac{r}{R}$$

$$n(z) = \frac{r}{R} n(r)$$

射线命名规则

地震震相命名遵循 IASPE 标准，采用字母拼接方式描述射线从震源到接收器的完整传播路径。

基本字母含义

符号	物理含义	备注
P	纵波	地壳/地幔中传播
S	横波	地壳/地幔中传播，不能通过流体核
K	外核中的 P 波	源自德语 Kern
I	内核中的 P 波
J	内核中的 S 波
L/R/G	面波/长周期波	非体波，通常单独分类

边界作用修饰符

符号	位置/作用	示例
p, s	震源处向上出射的初始分支	pP 表示先向上传到地表反射再向下
m	Moho 面地壳侧反射	PmP
c	核幔边界地幔侧反射	PcP
i	内外核边界外核侧反射	PKiKP
n	沿速度间断面的首波/折射波	Pn
g	地壳内直达波	Pg

常见震相分类与典型路径

地壳震相

• Pg, Sg：地壳直达波，速度结构敏感
• Pn, Sn：地幔顶盖首波，常用于 Moho 深度与上地幔顶速度反演
• PmP, SmS：Moho 面反射波

深震震相

• P, S
• PP, SS, PS, SP：地表反射/转换震相
• sP, pS：地表反射/转换震相

全球体波震相

版权所有

版权归属：Staaaaaaaaar

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

编辑此页

贡献者: Staaaaaaaaar

上一页弹性动力学基础

下一页射线的振幅和相位