摘要: 本贴从零开始学习正演的数值模拟方法.
1. 偏微分基础
1.1 导数
引例: 物体从一维坐标的原点开始移动, 在 t t t 时刻, 它在坐标轴的位置由函数 s ( t ) s(t) s(t) 确定, 则速度为位置变化量与时间的比值:
v ( t ) = d s ( t ) d t = lim Δ t → 0 s ( t + Δ t ) − s ( t ) Δ t (1) v(t) = \frac{\mathrm{d} s(t)}{\mathrm{d} t} = \lim_{\Delta t \to 0} \frac{s(t + \Delta t) - s(t)}{\Delta t} \tag{1} v(t)=dtds(t)=Δt→0limΔts(t+Δt)−s(t)(1)
加速度为速度变化量与时间的比值:
a ( t ) = d v ( t ) d t = lim Δ t → 0 v ( t ) − v ( t − Δ t ) Δ t = lim Δ t → 0 s ( t + Δ t ) − 2 s ( t ) + s ( t − Δ t ) Δ t 2 (2) a(t) = \frac{\mathrm{d} v(t)}{\mathrm{d} t} = \lim_{\Delta t \to 0} \frac{v(t) - v(t - \Delta t)}{\Delta t} = \lim_{\Delta t \to 0} \frac{s(t + \Delta t) - 2 s(t) + s(t - \Delta t)}{\Delta t^2} \tag{2} a(t)=dtdv(t)=Δt→0limΔtv(t)−v(t−Δt)=Δt→0limΔt2s(t+Δt)−2s(t)+s(t−Δt)(2)
推广 1: 给定一个单变量函数
y = f ( x ) (3) y = f(x) \tag{3} y=f(x)(3)
其一阶导数记为
y ′ = d f ( x ) d x (4) y' = \frac{\mathrm{d} f(x)}{\mathrm{d} x} \tag{4} y′=dxdf(x)(4)
二阶导数记为
y ′ ′ = d 2 f ( x ) d x 2 (5) y'' = \frac{\mathrm{d}^2 f(x)}{\mathrm{d} x^2} \tag{5} y′′=dx2d2f(x)(5)
1.2 偏导
给定一个二变量函数
z = f ( x , y ) (6) z = f(x, y) \tag{6} z=f(x,y)(6)
其针对 x x x 偏导的为
∂ z ∂ x = lim Δ x → 0 f ( x + Δ x , y ) − f ( x , y ) Δ x (7) \frac{\partial z}{\partial x} = \lim_{\Delta x \to 0} \frac{f(x + \Delta x, y) - f(x, y)}{\Delta x} \tag{7} ∂x∂z=Δx→0limΔxf(x+Δx,y)−f(x,y)(7)
即 x x x 发生了变化, 而 y y y 并没变化. 二阶偏导为
∂ 2 z ∂ x 2 = lim Δ x → 0 f ( x + Δ x , y ) − 2 f ( x , y ) + f ( x − Δ x , y ) Δ x 2 (8) \frac{\partial^2 z}{\partial x^2} = \lim_{\Delta x \to 0} \frac{f(x + \Delta x, y) - 2 f(x, y) + f(x - \Delta x, y)}{\Delta x^2} \tag{8} ∂x2∂2z=Δx→0limΔx2f(x+Δx,y)−2f(x,y)+f(x−Δx,y)(8)
另外有:
∂ 2 z ∂ x ∂ y = lim Δ x → 0 , Δ y → 0 f ( x + Δ x , y + Δ y ) − f ( x , y + Δ y ) − f ( x + Δ x , y ) + f ( x , y ) Δ x Δ y (9) \frac{\partial^2 z}{\partial x \partial y} = \lim_{\Delta x \to 0, \Delta y \to 0} \frac{f(x + \Delta x, y + \Delta y) - f(x, y + \Delta y) - f(x + \Delta x, y) + f(x, y)}{\Delta x \Delta y} \tag{9} ∂x∂y∂2z=Δx→0,Δy→0limΔxΔyf(x+Δx,y+Δy)−f(x,y+Δy)−f(x+Δx,y)+f(x,y)(9)
∂ 2 z ∂ y ∂ x = ∂ 2 z ∂ x ∂ y (10) \frac{\partial^2 z}{\partial y \partial x} = \frac{\partial^2 z}{\partial x \partial y} \tag{10} ∂y∂x∂2z=∂x∂y∂2z(10)
在进行数值模拟的时候, 不可能取 Δ x → 0 \Delta x \to 0 Δx→0, 因此 (8) 式简化为
∂ 2 z ∂ x 2 ≈ f ( x + Δ x , y ) − 2 f ( x , y ) + f ( x − Δ x , y ) Δ x 2 (11) \frac{\partial^2 z}{\partial x^2} \approx \frac{f(x + \Delta x, y) - 2 f(x, y) + f(x - \Delta x, y)}{\Delta x^2} \tag{11} ∂x2∂2z≈Δx2f(x+Δx,y)−2f(x,y)+f(x−Δx,y)(11)
其中 Δ x \Delta x Δx 越小越准确, 但涉及的计算量越大, 我们只能取一个折中.
注 1: 为统一起见, 即使一元函数, 以后也常使用 ∂ \partial ∂ 而不是 d \mathrm{d} d.
1.3 泰勒级数
然而, (11) 式本质上是有问题的. 注意到 Δ x → 0 \Delta x \to 0 Δx→0 已经不成立, 且在实际地震数据中, 它可能是 5m 甚至 20m. 因此, 需要用到泰勒级数. 当函数 f ( x ) f(x) f(x) 在 x 0 x_0 x0 处存在直到 n n n 阶的导数, 则
f ( x ) = ∑ i = 0 n f ( i ) ( x 0 ) ( x − x 0 ) i i ! + o ( ( x − x 0 ) n ) (12) f(x) = \sum_{i = 0}^n \frac{f^{(i)}(x_0)(x - x_0)^i}{i!} + o((x - x_0)^n) \tag{12} f(x)=i=0∑ni!f(i)(x0)(x−x0)i+o((x−x0)n)(12)
其中 f ( i ) ( x 0 ) / i ! f^{(i)}(x_0) / i! f(i)(x0)/i! 称为泰勒展开式的系数.
直观的解释: 如果函数 f ( x ) f(x) f(x) 不是线性的, 则它的变化量不仅与斜率有关, 而且与斜率的变化率也有关. 更多内容就只有自己去找数学书啃了.
这里需要继续写 …
2. 波动方程
2.1 弦振动 (横波) 方程
参见全波形反演的深度学习方法: 第 2 章 正演, 根据牛顿第二定律
F = m a (12) F = ma \tag{12} F=ma(12)
弦振动方程为
∂ 2 u ( x , t ) ∂ t 2 = c 2 ∂ 2 u ( x , t ) ∂ x 2 + f ( x , t ) (13) \frac{\partial^2 u(x, t)}{\partial t^2} = c^2 \frac{\partial^2 u(x, t)}{\partial x^2} + f(x, t) \tag{13} ∂t2∂2u(x,t)=c2∂x2∂2u(x,t)+f(x,t)(13)
其中 c 2 = T / ρ c^2 = T / \rho c2=T/ρ, f ( x , t ) = F ( x , t ) / ρ f(x, t) = F(x, t) / \rho f(x,t)=F(x,t)/ρ, 左式的物理意义是瞬时加速度 a a a, 右式第一项的物理意义是 单位质量所受的力 F F F, c c c 的物理意义是速度.
进一步忽略重力 F ( x , t ) F(x, t) F(x,t) 的作用, 可以推出一维齐次波动方程的解:
∂ 2 u ( x , t ) ∂ x 2 = 1 c 2 ∂ 2 u ( x , t ) ∂ t 2 (14) \frac{\partial^2 u(x, t)}{\partial x^2} = \frac{1}{c^2} \frac{\partial^2 u(x, t)}{\partial t^2} \tag{14} ∂x2∂2u(x,t)=c21∂t2∂2u(x,t)(14)
2.2 声波 (纵波) 方程
声波仅有纵波. 考虑二维的情况, 它满足
1 v 2 ∂ 2 U ∂ t 2 = ∂ 2 U ∂ x 2 + ∂ 2 U ∂ z 2 (15) \frac{1}{v^2} \frac{\partial^2 U}{\partial t^2} = \frac{\partial^2 U}{\partial x^2} + \frac{\partial^2 U}{\partial z^2} \tag{15} v21∂t2∂2U=∂x2∂2U+∂z2∂2U(15)
其中 U U U 指压力.
为了便于数值模拟, 将平面进行离散化, 仅考虑某些网格交叉点, 质量、压力等仅存在于这些点 (称为质点, 不知是否专业). 这样, 我们只考察第 i i i 行第 j j j 列的质点在时间 k k k 的压力
U i , j k (16) U_{i, j}^k \tag{16} Ui,jk(16)
在最简化的情况, 将 (11) 式按照变量名改造后代入 (15) 式可得
1 v 2 U i , j k + 1 − 2 U i , j k + U i , j k − 1 Δ t 2 = U i + 1 , j k − 2 U i , j k + U i − 1 , j k Δ x 2 + U i , j + 1 k − 2 U i , j k + U i , j − 1 k Δ y 2 (17) \frac{1}{v^2} \frac{U_{i, j}^{k + 1} - 2 U_{i, j}^{k} + U_{i, j}^{k - 1}}{\Delta t^2} = \frac{U_{i + 1, j}^k - 2 U_{i, j}^{k} + U_{i - 1, j}^k}{\Delta x^2} + \frac{U_{i, j + 1}^k - 2 U_{i, j}^{k} + U_{i, j - 1}^k}{\Delta y^2} \tag{17} v21Δt2Ui,jk+1−2Ui,jk+Ui,jk−1=Δx2Ui+1,jk−2Ui,jk+Ui−1,jk+Δy2Ui,j+1k−2Ui,jk+Ui,j−1k(17)
其中 k + 1 k + 1 k+1 表示下一个时间点, i + 1 i + 1 i+1 表示下一个质点.
