第八章 波动方程数值解

波动方程数值解

本章节介绍了求解声波方程和弹性波方程数值解的一些方法。这些方法求取的解中包含有初至和多次散射等所有的同相轴，因此常被用于逆时偏移（Reverse time migration, RTM）和波形反演（Waveform inversion）。我们将概述规则网格（Kelly et al.，1976）和交错网格（Virieux，1984和1986；Levander，1988）下的有限差分方法。为了获得更高的精度，需要更高阶的格式，因此现已经开发出比有限差分法更精确的谱方法（Spectral collocation）。Trefehen（2000）写道：“如果想要在一个简单的域上高精度地求解ODE或PDE，并且如果定义问题的数据是平滑的，那么谱方法通常是最好的选择。它们通常可以达到十位数的精度，而有限差分法或有限元法则仅可以达到两位数或三位数精度。在精度较低的情况下，其计算机内存消耗与其他方案相比更少。”因此，本章还介绍了伪谱（Pseudospetral）和谱元法（Spectral element）。伪谱法（Kosloff和Baysal，1982；Kosloff等人，1984）精度最高，但当界面两侧阻抗差异较大时可能会出现明显的数值误差。这一问题主要可以通过谱元法（Komatitsch et al.，2005；Kudela et al.，2007）来解决，并且其求解仍然可保持高精度和可观的效率。

8.1 有限差分法

表8.1中展示了一阶、二阶偏导数不同格式、不同精度的有限差分（Finite difference, FD）近似，其精度的阶数可通过泰勒展开来估计。

$$\begin{align} \notag \frac{f(x+\Delta x)-f(x)}{\Delta x} &=\frac{1}{\Delta x}[f(x) + \frac{df(x)}{dx} \Delta x \\ \notag &+ \frac{1}{2} \frac{d^2f(x)}{dx^2} (\Delta x)^2 + ... - f(x)] \\ \notag &= \frac{df(x)}{dx} +1/2 \frac{d^2f(x)}{dx^2} \Delta x + ..., \\ \notag &= \frac{df(x)}{dx} + O(\Delta x) \end{align} \tag{8.1}$$

这意味着中点在$x$处的FD近似在采样间隔$\Delta x$中是一阶精确的，其中项$O(x)$表示截断误差。类似地，对二阶导数$d[df(x)/dx]/dx$的中心FD近似可以通过对括号中的项应用一阶精确的前向FD近似和对括号外项应用一阶准确的后向FD近似来导出。

8.1.1 波动方程的有限差分近似

常密度介质中的二维声波方程公式为：

$$\frac{\partial^2p(\mathbf x, t^{\prime})}{\partial x^2} + \frac{\partial^2p(\mathbf x, t^{\prime})}{\partial z^2} - \frac {1}{c(x,z)^2} \frac{\partial^2p(\mathbf x, t^{\prime})}{\partial t^{\prime}} =f(\mathbf x, t^{\prime}), \tag{8.2}$$

其中，$c(x,z)^2$为速度模型，$p(\mathbf x, t^{\prime})$为压力场，$f(\mathbf x, t^{\prime})$为非均匀震源项。连续的波动方程及其解可以离散在均匀分布的时间-空间域网格上：

$$(x,z,t^{\prime}) \leftrightarrow (i\Delta x, j\Delta z, t\Delta t), \\ p(x,z,t^{\prime}) \leftrightarrow p^t_{ij}, \\ f(x,z,t^{\prime}) \leftrightarrow f^t_{ij}, \\ c(x,z) \leftrightarrow c_{ij} \tag{8.3}$$

方便起见，假设垂向空间步长$\Delta z$等于水平空间步长$\Delta x$。可以使用一个以$x, y, t$为坐标轴的三维数据体将离散解可视化（如图8.1所示）。数据体的深度、宽度和时间长度分别为$M\Delta z=D, N\Delta x = X, L\Delta t = T$，其中$M,N,L$均为整数。
求解$(i,j,t)$的压力场，并使用二阶中心差分来近似公式8.2中的二阶偏导数可得:

$$\begin{align} \notag \frac{\partial p^2}{\partial t^{\prime}} \leftrightarrow [p^{t+1}_{ij}-2p^{t}_{ij}+p^{t-1}_{ij}]/(\Delta t)^2,\\ \notag \frac{\partial p^2}{\partial x^{\prime}} \leftrightarrow [p^{t}_{i+1j}-2p^{t}_{ij}+p^{t}_{i-1j}]/(\Delta x)^2,\\ \notag \frac{\partial p^2}{\partial z^{\prime}} \leftrightarrow [p^{t}_{ij+1}-2p^{t}_{ij}+p^{t}_{ij-1}]/(\Delta z)^2 \end{align} \tag{8.4}$$

整理可得波动方程离散时间递推格式：

$$\begin{align} \notag for &\space t = 2:nt \\ \notag &p^{t+1}_{ij} = 2p^{t}_{ij}-p^{t-1}_ij \\ \notag &+a[p^t_{i+1j}-2p^t_{ij}+p^t_{i-1j}] \\ \notag &+a[p^t_{ij+1}-2p^t_{ij}+p^t_{ij-1}]-(\Delta x)^2af^t_{ij}\\ \notag end \end{align} \tag{8.5}$$

其中，$a=(c_{ij}\Delta t/\Delta x)^2$。这里假定初始条件为$p(\mathbf x, t^{\prime})=0$，所以源波场是由body-force项$f^{t}_{ij}$所激发产生的。
从$t$变换到坐标$−t$，即反转时间，波动方程8.2的形式将保持不变。这意味着时间反转后，波动方程的响应将与原始方程相同，只是时间变量的符号发生了变化。在这种情况下，图8.2中的前向光锥将变成后向光锥。重要的一点是，通过在时间上向后求解而不是在时间上向前求解，可以用有限差分方程生成后向光锥。这种后向锥体可用于从后来的数据重建早期的波场，这些数据被称为后向外推（backward-extrapolated）波场（第11章）。

8.1.2 稳定性和精度分析

经验发现，如果每个最小波长至少有10个点，并且传播距离不超过十几个波长，则2-2（即时间和空间中二阶精度）差分格式的精度（Kelly et al.，1976）在均匀介质中是可接受的。在非均质介质中，每个波长通常使用15-20个点。更高阶精度，比如2-4差分格式需要大约五个点/波长（Levander，1988），但在非均匀介质中每波长需要10至15点。如果网格点间距太粗，则倾斜界面显示为阶梯，每个阶梯的边缘都会成为强散射点。如果波传播超过30个波长，经验表明需要2–8或更高的差分精度。事实上，Stork（2013）声称，在传播50个波长后，波的时间误差不应超过四分之一周期，否则波形反演将会出问题。

CFL（Courant Friedrichs-Lewy）稳定性条件中，满足波动方程和$t=0$初始条件的有限差分解的数值依赖域DOD（Domain of dependence, DOD）必须比解析DOD（Mitchell和Griffiths，1980）有更大（如图8.3所示）来确定。否则，FD解将忽略影响$(x_0, t_0)$处解的一些初始条件。这种忽略将导致不稳定的有限差分解。为了避免这种ignorance，数值传播速度（$\Delta x / \Delta t$）必须要比实际传播速度$c$更快。这一条件等价于一维的CFL稳定性准则：

$$\frac{1}{c} > \frac{\Delta x}{\Delta t}, \tag{8.6}$$

在更高维度的空间中，具有如下形式：

$$\frac{1}{\sqrt N} > \frac{c\Delta x}{\Delta t}, \tag{8.7}$$

其中，$N=1,2,3$表示空间维度，$\Delta x$需要满足精度条件，$\Delta t$需要满足条件稳定性条件。

有限差分格式不稳定指的是误差随着时间递推而逐渐增大。例如，假设一阶双曲波动方程的解是归一化指数$p^t_j=e^{ikx_j}$由指数$t$表示的特定时间步长。将该指数代入一维波动方程的FD递推格式中，可得到其迭代解$p^{t+p}_j=\lambda ^p e^{ikx_j}$,其中$\lambda$是放大系数，取决于时间和空间采样间隔。如果放大系数|λ|＞1，差分解不稳定。

8.1.3 基于Matlab的二维声波方程有限差分求解

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% (NX,NZ,NT) - input- (Horizontal,Vertical) gridpt dimens. of vel
% model & # Time Steps
% FR - input- Peak frequency of Ricker wavelet
% BVEL - input- NXxNZ matrix of background velocity model
% (dx,dt) - input- (space, time) sample intervals
% (xs,zs) - input- (x,z) coordinates of line source
% RICKER(NT) - input- NT vector of source time histories
% (p2,p1,p0) - calc.- (future,present,past) NXxNZ matrices of
% modeled pressure field
% (p0,p1) -output- Old and present pressure panels at time NT.
% REALDATA(NX,NT) -output- CSG seismograms at z=2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c=4.0;FRE=20;
NX=300;NZ=NX; dx=c/FRE/20;dt=.5*dx/c;
xs=round(NX/2.3); zs=round(NX/2);NT=600;
t=[0:1:NT-1]*dt-0.95/FRE;RICKER=zeros(length(t));
RICKER= (1-t .*t * FRE^2 *pi^2 ) .*exp(- t.^2 * pi^2 * FRE^2 ) ;
plot([0:NT-1]*dt,RICKER);
title(’Ricker Wavelet’);xlabel(’Time (s)’)
BVEL=ones(NX,NZ)*c;
BVEL(NX-round(NX/2):NX,:)= BVEL(NX-round(NX/2):NX,:)*1.2;
REALDATA=zeros(NX,NT);
p0=zeros(NX,NZ);p1=p0;p2=p0;
cns=(dt/dx*BVEL).^2;
NX=200;NZ=NX;
%for it=NT:-1 %Reverse Time Looping
for it=1:1:NT %Forward Time Looping
    p2 = 2*p1 - p0 + cns.*del2(p1);
    p2(xs,zs) = p2(xs,zs) + RICKER(it);
    REALDATA(:,it) = p2(xs,:)’;
    p0=p1;p1=p2;
    if round(it/20)*20==it;p00=p0/max(abs(p0(:))+.001);
        imagesc([1:NX]*dx,[1:NX]*dx,(p00+BVEL)); colorbar;
        pause(.1);
    end
end;
p1=p0;p0=p2;
title(’Snapshot of Acoustic Waves’)
xlabel(’X (km)’)
ylabel(’Z (km)’)

8.2 伪谱法求解波动方程

利用有限差分法求解波动方程有一个问题，就是当波的传播距离较长时，数值频散会使解产生误差。为了减少这些误差，需要使用更高阶的有限差分格式，但这会导致差分格式更长且消耗更多的计算时间。一种有效地将离散导数算子的范围增加到计算网格宽度的方法是伪谱方法（Kreiss & Oliger，1972；Kosloff & Baysal，1982；Koslopf et al.，1984），其在傅里叶域中计算空间导数。理论上，每个波长只需要两个点即可在平滑的介质中对波动方程进行精确求解。
对于如下形式的一维波动方程：

$$\frac{\partial p^2}{\partial t} = \rho c^2 \frac{\partial}{\partial x} (\frac {1}{\rho} \frac{\partial p}{\partial x}) \tag{8.8}$$

其中，$\rho$为密度。
伪谱法求解的流程为：

1. 计算空间偏导数$\frac{1}{\rho} \frac{\partial p}{\partial x}$：对压力场做快速傅里叶变换（FFT）得到$P$，然后乘以$ik_x$。
2. 对$ik_xP$做逆FFT，其结果乘以$1/ \rho$得到$\frac{1}{\rho} \partial p / \partial x$。然后再对结果进行上述操作得到空间偏导数$\partial(1/\rho \partial p/\partial x)/\partial x$。
3. 通过二阶精度有限差分近似波动方程中的二阶时间偏导数。下一时刻的波长可以通过当前和前一时刻波场值计算得到。
4. 计算时间加1，重复步骤1-4直到设置的波场持续时间。

空间导数做FFT后表示波数乘以$ik_x$，其在空间上等价于与计算域一样宽的离散导数算子。它需要周期边界条件，但吸收边界条件（见附录8A和Liu & Sen，2010年；2012年）也可以。

8.2.1 基于Matlab的二维声波方程伪谱法求解

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pseudospectral solution of 1D acoustic wave equation by Kai Lu
% method=1 pseudospectral method
% method=2 finite difference method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; hold off;
method=1;
c=4.0;FRE=20;nx=2000;
dx=c/FRE/10; %set dx
dt=.3*dx/c; %set dt
nt=2000;
sx=nx/2;
thick=300; % sponge thickness
alpha=pi/5/thick;
factor=9; % sponge factor
% generate wavelet
t=(-200:1:nt-201)*dt-0.95/FRE;
ricker=(1-t.*t*FRE^2*pi^2).*exp(-t.*t*pi^2*FRE^2/2);
plot((0:nt-1)*dt,ricker);title(’Ricker Wavelet’);xlabel(’Time (s)’);
pause(1);
% initiate wavefield
p2=zeros(1,nx); p1=p2;p0=p2;
% calculate wavenumber
k=zeros(1,nx);
for i=1:nx
    k(i)=2*pi*(i-1)/nx/dx;
end;
for i=nx/2+2:nx
    k(i)=k(i)/(i-1)*(i-nx-1);
end;
% extrapolate wavefield

for it=1:nt;

    if method==1;
        % pseudo-spectral method
        P1=fft(p1); DP1=-k.^2.*P1; dp1=real(ifft(DP1));
        p2=2*p1-p0+c^2*dt^2*dp1;
    elseif method==2
        % finite difference method
        p2(2:nx-1)=2*p1(2:nx-1)-p0(2:nx-1)+(c*dt/dx)^2*(p1(3:nx)-2*p1(2:nx-1)+p1(1:nx-2));
    end;
    
    % add source
    if (it==1);
        p2(sx)=p2(sx)+ricker(it);
    else
        p2(sx)=p2(sx)+ricker(it)-ricker(it-1);
    end;

    % sponge ABC
    for i=1:thick
        p2(i)=p2(i)*(cos(alpha*(thick+1-i))^0.5+factor)/(factor+1);
    end;
    for i=nx-thick+1:nx;
        p2(i)=p2(i)*(cos(alpha*(i-nx+thick))^0.5+factor)/(factor+1);
    end;

    % free surface ABC
    % p2(nx)=0;
    % p2(1)=0;
    % update wavefield
    p0=p1;p1=p2;
    plot((0:nx-1-2*thick),real(p2(thick+1:nx-thick)));
    axis([0 nx-1-2*thick -2.5 2.5]);
    title(’Pressure snapshot’);
    xlabel(’Distance’);
    pause(0.01);
end;

伪光谱法的主要缺点：当波跨越强变化界面时会产生Gibbs-like误差。这一问题限制了伪谱法在地球物理中的应用，因为其需要模型是平滑的，但谱精度的原理在谱元法中起着重要作用。

8.2.2 稳定性和精度分析

CFL稳定性条件是必要条件而不是充分条件。因此，我们需要用von Neumann来进一步进行稳定性分析。首先假设一维波动方程的解为平面波$p_j=e^{i(kj\Delta x-n\omega\Delta t)}$，其中$j$和$n$表示空间和时间索引。利用FFT可得到二阶偏导数为：

$$\frac{\partial ^2 p_j}{\partial x^2} = -k^2 e^{i(kj\Delta x-n\omega\Delta t)} \tag{8.9}$$

然后将平面波解带入二阶时间偏导数公式8.4中，可得

$$\frac{\partial ^2 p_j}{\partial t^2}=-\frac{4}{\Delta t^2} {\rm sin}^2 (\omega \Delta t/2) e^{i(kj\Delta x-n\omega\Delta t)} \tag{8.10}$$

将以上公式带入波动方程可得离散数值频散关系为$k=\frac{2}{c\Delta t} {\rm sin} \frac{\omega \Delta t}{2}$，其反函数为：

$$\omega = \frac{2}{\Delta t}{\rm sin}^{-1} (\frac{kc\Delta t}{2}) \tag{8.11}$$

当$\Delta t <<1$且如下公式成立时，公式8.11很好地近似了实际的频散关系$\omega=kc$

$$\frac{k_{max}c\Delta t}{2} <1 \tag{8.12}$$

其中，$k_{max}=\pi / \Delta x$是unaliased平面波的最大波长。这种情况下，频散非常小且对于均匀介质是稳定的。但如果速度场不是平滑的，而是有明显的阻抗变化，反射波场中会出现明显的误差（Trefethen, 2000）。

8.3 谱元法求解波动方程

在2.5D和3D地球弹性模型中，常利用谱元法（Spectral element）求解波动方程（Komatitsch et al.，2005）。与标准有限差分法相比，其可以使用不规则网格以适应不规则自由起伏地形和地下界面，并且不会产生假的阶梯散射。与有限差分使用空间差分相类似，如果谱元法积分的阶数足够高，则频散很小。 由于质量矩阵是对角阵，因此谱元法比具有集中质量近似（lumped-mass approximation）的标准有限元方法（Marfurt，1984）更精确。对于具有不规则界面的模型，由于谱元的边界可以与界面几何相适应，因此谱元法建模也比伪谱更准确。