三维波动方程与施图姆–刘维尔（S-L）理论

在三维空间中，弹性波方程（Navier 方程）是一个矢量方程。施图姆–刘维尔（S-L）理论的应用，通常是通过分离变量法，将复杂的矢量场分解为一系列互相正交的标量场或矢量基函数的过程。

在地震学中，这最典型的应用是处理球对称地球模型（如 PREM 模型）的简正振型（Normal Modes）或分层介质中的波传播。

1. 弹性波方程的算子化
   各向同性均匀介质中的三维弹性波方程为：

ρ∂2u∂t2=(λ+μ)∇(∇⋅u)+μ∇2u\rho \frac{\partial^2 \mathbf{u}}{\partial t^2} = (\lambda + \mu) \nabla (\nabla \cdot \mathbf{u}) + \mu \nabla^2 \mathbf{u}ρ∂t2∂2u​=(λ+μ)∇(∇⋅u)+μ∇2u

其中u(x,t)\mathbf{u}(\mathbf{x}, t)u(x,t)是位移矢量。

为了应用 S-L 理论，我们考虑定态解（简谐波），令u(x,t)=U(x)e−iωt\mathbf{u}(\mathbf{x}, t) = \mathbf{U}(\mathbf{x}) e^{-i\omega t}u(x,t)=U(x)e−iωt，得到空间部分的特征方程：

−ω2ρU=(λ+μ)∇(∇⋅U)+μ∇2U- \omega^2 \rho \mathbf{U} = (\lambda + \mu) \nabla (\nabla \cdot \mathbf{U}) + \mu \nabla^2 \mathbf{U}−ω2ρU=(λ+μ)∇(∇⋅U)+μ∇2U

这可以写成线性算子形式：LU=ω2U\mathcal{L} \mathbf{U} = \omega^2 \mathbf{U}LU=ω2U。这里的ω2\omega^2ω2就对应 S-L 问题中的本征值。

2. 分离变量与标量化
   对于三维问题，我们通常利用**亥姆霍兹分解（Helmholtz Decomposition）**将位移矢量U\mathbf{U}U分解为标量位ϕ\phiϕ（纵波 P）和矢量位Ψ\boldsymbol{\Psi}Ψ（横波 S）：

U=∇ϕ+∇×Ψ,∇⋅Ψ=0\mathbf{U} = \nabla \phi + \nabla \times \boldsymbol{\Psi}, \quad \nabla \cdot \boldsymbol{\Psi} = 0U=∇ϕ+∇×Ψ,∇⋅Ψ=0

以球坐标系(r,θ,ϕ)(r, \theta, \phi)(r,θ,ϕ)为例，对于球对称介质，位移场可以分解为：

类扭转振荡（Toroidal modes）：只有横波分量，对应SHSHSH波。

类球振荡（Spheroidal modes）：包含纵波和垂直横波分量，对应P−SVP-SVP−SV波。

3. 径向方程：标准的 S-L 问题
   以类扭转振荡为例，其位移场可以表示为：

UT=W(r)[r^×∇Ylm(θ,ϕ)]\mathbf{U}_T = W(r) [ \hat{\mathbf{r}} \times \nabla Y_{lm}(\theta, \phi) ]UT​=W(r)[r^×∇Ylm​(θ,ϕ)]

其中YlmY_{lm}Ylm​是球谐函数。将此形式代入波动方程，角度部分被球谐函数抵消，剩下的径向部分W(r)W(r)W(r)满足如下方程：

ddr[μr4dWdr]+[ω2ρr4−(l−1)(l+2)μr2]W=0\frac{d}{dr} \left[ \mu r^4 \frac{dW}{dr} \right] + \left[ \omega^2 \rho r^4 - (l-1)(l+2)\mu r^2 \right] W = 0drd​[μr4drdW​]+[ω2ρr4−(l−1)(l+2)μr2]W=0

如果我们对比标准的 S-L 形式：ddr[p(r)dWdr]+[q(r)+λρ(r)]W=0\frac{d}{dr} [p(r) \frac{dW}{dr}] + [q(r) + \lambda \rho(r)] W = 0drd​[p(r)drdW​]+[q(r)+λρ(r)]W=0，可以清晰地看到：

本征值：λ=ω2\lambda = \omega^2λ=ω2（频率的平方）。

系数函数：p(r)=μ(r)r4p(r) = \mu(r) r^4p(r)=μ(r)r4。

权重函数：ρSL(r)=ρ(r)r4\rho_{SL}(r) = \rho(r) r^4ρSL​(r)=ρ(r)r4。

势函数项：q(r)=−(l−1)(l+2)μ(r)r2q(r) = -(l-1)(l+2)\mu(r) r^2q(r)=−(l−1)(l+2)μ(r)r2。

4. 边界条件与奇异性
   在地球物理应用中，S-L 问题必须配合边界条件才能定解：

自由表面条件：在地球表面r=Rr=Rr=R，应力为零，即Tr=μ(dWdr−Wr)=0T_r = \mu (\frac{dW}{dr} - \frac{W}{r}) = 0Tr​=μ(drdW​−rW​)=0。

内部连续性：在介质分界面处（如核幔边界），位移和应力必须连续。

地心奇点：在r=0r=0r=0处，p(0)=0p(0)=0p(0)=0，这正是你图片中提到的奇异边界点。根据 S-L 理论，我们要求解在r=0r=0r=0处保持有限。

5. 为什么这样做？（应用价值）
   正交性关系：

S-L 理论保证了不同频率的振型Un\mathbf{U}_nUn​和Um\mathbf{U}_mUm​在全地球体积VVV上关于质量密度ρ\rhoρ正交：

∭Vρ(x)Un∗⋅UmdV=δnm\iiint_V \rho(\mathbf{x}) \mathbf{U}_n^* \cdot \mathbf{U}_m dV = \delta_{nm}∭V​ρ(x)Un∗​⋅Um​dV=δnm​

这意味着我们可以把任何地震位移记录分解为这些固有振型的叠加。

合成地震图（Synthetic Seismograms）：

在 CPS 等软件中，计算理论地震图的核心逻辑就是：先通过求解 S-L 方程组得到本征值（色散曲线）和本征函数（振型深度分布），然后进行模态叠加（Mode Summation）。

反演地球内部结构：

如果我们观测到了地球的固有频率ωobs\omega_{obs}ωobs​，我们可以根据 S-L 方程建立ω\omegaω对参数μ(r),ρ(r)\mu(r), \rho(r)μ(r),ρ(r)的敏感核，从而反推地球深部的物理性质。