波动方程的格林函数解数学推导

1.在无限均匀介质中含震源项的标量波动方程

假设在无限均匀介质中，含震源项的标量声波方程可以表示为下式：

[frac{part{varphi}(mathbf{x},t)}{part{t^2}}=c^2nabla^2{varphi(mathbf{x},t)}+f(mathbf{x},t) tag{1.1} ]

初值条件为：

[varphi(mathbf{r},t)=0 space (tleq0) tag{1.1a} ]

由于假设其在无限大自由空间，其边界条件满足索末菲辐射条件(Sommerfeld Radiation Condition)。其中

其中声波波速(c) 是一个常数；(f(mathbf{x},t)) 是震源项。我们下面讨论两种震源函数情况的格林函数解：（1）点脉冲震源项(Impulsive sources); （2）任意形式的震源(Arbitrary sources)。

首先讨论点脉冲震源，在这种情况下，其震源项(f(mathbf{x},t)) 可以写成如下形式：

[f(mathbf{x},t)=delta(mathbf{x}-mathbf{zeta})delta(t-t_0)tag{1.2} ]

其中：(delta(t-t_0))是时间域脉冲函数，其中(t_0)是震源项发震时刻； (delta(mathbf{x}-zeta))是空间脉冲函数，其中(mathbf{zeta}=[zeta_1,zeta_2,zeta_3]^T) 是震源的空间位矢，在笛卡尔坐标的表达如下：

[delta{(mathbf{x}-mathbf{zeta})}=delta{(x_1-zeta_1)}delta(x_2-zeta_2)delta(x_3-zeta_3) tag{1.3} ]

在此情况下，其方程（1.1）中(varphi(mathbf{x},t)) 解可以表达成格林函数（Green Function），其方程可以写为形式：

[frac{part^2{G}}{part{t^2}}=c^2nabla^2G+delta(mathbf{x}-zeta)delta(t-t_0) tag{1.4} ]

其中(G=G(mathbf{x},t;zeta,t_0)) 为格林函数解，由于因果性，(tleq t_0) 时，其格林函数(G(mathbf{x},t;zeta,t_0)=0)

下面我们利用Fourier变换求解三维标量方程在均匀无限介质内含点脉冲震源项的格林函数解，其推导如下：

利用三重Fourier变换对式（1.1）方程中的(mathbf{x}) 进行傅里叶变换到波数域：

[hat{G}(mathbf{k},t;zeta,t_0)=mathcal{F}{G(mathbf{x},t;zeta,t_0)} tag{1.5} ]

其中：(mathcal{F}) 为傅里叶算子；(mathbf{k}=[k_1,k_2,k_3]^T) 为波数域向量。代入方程（1.5）结合Fourier变换的微分性质及位移性质可以得到：

[begin{equation} frac{part^2{hat{G}}}{part{t^2}}+c^2k^2hat{G}=e^{imathbf{k}cdotmathbf{zeta}}delta{(t-t_0)} tag{1.6} end{equation} ]

其中：(k) 为波数向量(mathbf{k})的模量，即(k=|mathbf{k}|=sqrt{k_1^2+k_2^2+k_3^2}) 。根据时间因果性我们得到(hat{G})也具有因果性，即

[hat{G}(mathbf{k},t;mathbf{zeta}_0,t_0)=0; space (t<t_0) tag{1.7} ]

由于Draclet‘s函数的性质 ，在(t>t_0) 时，(delta{(t-t_0)}=0) 。由此（1.6）的方程右边项为0，可以得到如下二阶常微分方程：

[frac{part^2{hat{G}}}{part{t^2}}+c^2k^2hat{G}=0; space(t>t_0) tag{1.8} ]

对于(mathbf{k}) 和 (mathbf{x}_0)而言，在式（1.8）中二阶微分方程中解(hat{G}) 可以写成如下形式：

[hat{G}(mathbf{k},t;zeta,t_0)=Acos{ck(t-t_0)}+Bsin{ck(t-t_0)}space (t>t_0) tag{1.9} ]

其中：(A)和(B) 是关于(mathbf{x}_0)，(t_0), (mathbf{k}) 的函数。为了确定(A)和(B) ，必须考虑其(t=t_0)的初值条件，由于(hat{G})在(t=t_0) 处连续：

[lim_{trightarrow t_{0}^{+}}G(mathbf{k},t;zeta,t_0)=lim_{trightarrow t_{0}^{-}}G(mathbf{k},t;zeta,t_0)=0 tag{1.10} ]

结合式（1.9）和（1.7）可得：

[G(mathbf{k},t_0;zeta,t_0)=A=0 tag{1.11} ]

由此可以确定参数(A=0). 下面考虑方程解第二待定系数(B), 为了获取(B) ,对式（1.6）在(t_0) 初的一个小领域([t_0-d{t},t_0+d{t}]) 简记为([t_0^{-},t_0^{+}])， 对时间变量(t)进行积分以将方程右边的(delta{(t-t_0)}) 转化为1：

[int_{t_0^{-}}^{t_0^{+}}frac{part^2{hat{G}}}{part{t^2}}dt+c^2k^2int_{t_0^{-}}^{t_0^{+}}hat{G}dt=e^{imathbf{k}cdotmathbf{zeta}} tag{1.12} ]

其将二阶导数积分后化为一阶导数，其形式化为：

[left.frac{part{hat{G}}}{part{t}}right|_{t_0^{-}}^{t_0^{+}}++c^2k^2int_{t_0^{-}}^{t_0^{+}}hat{G}dt=e^{imathbf{k}cdotmathbf{zeta}} tag{1.13} ]

考虑到(hat{G}) 在(t=t_0) 处的连续性，考虑式（1.7）的(t<t_0) 时(hat{G}=0) 可知式(1.13)的右边第二项当(t_0^{+},t_0^{-}rightarrow t_0) 时为0，即：

[lim_{t_0^{+},t_0^{-}rightarrow0}c^2k^2int_{t_0^{-}}^{t_0^{+}}hat{G}d{t}=0 tag{1.14} ]

同时，由此可以写为如下形式：

[lim_{t_0^{+},t_0^{-}rightarrow{t_0}}[frac{part{hat{G}(t_0^{+})}}{part{t}}-frac{parthat{G}(t_0^{-})}{part{t}}]=e^{imathbf{k}cdotmathbf{zeta}} tag{1.15} ]

由于在(t_0^{-}<t_0) 时，(hat{G}=0) ，由此可知，(frac{part{hat{G}(t_0^{-})}}{part{t}}=0)，由于(t^{+}_0>t_0)，(frac{parthat{G}(t_0^{+})}{part{t}}) 可参考式（1.9）推出：

[lim_{t_0^{+}rightarrow t_0}[ckBcos(t-t_0)]=e^{imathbf{k}cdotmathbf{zeta}} tag{1.16} ]

可以推出参数(B)

[B=frac{e^{imathbf{k}cdotmathbf{zeta}}}{ck} tag{1.17} ]

将式（1.11）和（1.17）中求出的系数(A,B)代入式（1.9）可以推出 (hat{G}(mathbf{k},t;zeta,t_0))的表达式：

[hat{G}(mathbf{k},t;zeta,t_0)=frac{e^{imathbf{k}cdotmathbf{zeta}}}{ck}sin(ck(t-t_0)) tag{1.18} ]

在此，采用三重Fourier逆变换将(hat{G}(mathbf{k},t;zeta,t_0)) 的波数域转换回空间域(G(mathbf{x},t;zeta,t_0)) :

[G(mathbf{x},t;zeta,t_0)=frac{1}{(2pi)^3}iiintfrac{sin(ck(t-t_0))}{ck}e^{-imathbf{k}cdot(mathbf{x}-mathbf{zeta})}d{k_x}d{k_y}d{k_z} tag{1.19} ]

将(e^{imathbf{k}cdot(mathbf{x}-mathbf{zeta})}) 进行展开,首先定义位置向量(mathbf{r}=mathbf{x}-zeta)，其模量为(r=|mathbf{r}|)，其(mathbf{k})与(mathbf{r})夹角为(theta)：

[e^{-imathbf{k}cdot(mathbf{x}-mathbf{zeta})}=e^{-ikrcos(theta)} tag{1.20} ]

为了求取上式的积分值，将积分空间k由直角坐标系转换为球面坐标系描述，为计算方便，选取(mathbf{r})作为坐标，由此球面坐标积分微元(dV)可以写为：

[dV=dk_xdk_ydk_z=ksintheta dkd{theta}d{phi} tag{1.21} ]

其中(phi) 的范围为([0,2pi])，(theta) 的范围为([0,pi])和(k) 的范围为([0,infty))。由此式（1.19）z可以写为如下的球坐标形式的三重积分形式：

[G(mathbf{x},t;zeta,t_0)=frac{1}{(2pi)^3}int_{0}^{2pi}int_0^{infty}int_0^{pi}kfrac{sin ck(t-t_0)}{c}e^{-ikrcos(theta)}sin{theta}d{theta}d{k}d{phi} tag{1.22} ]

首先，考虑最内侧关于(theta) 的积分：

[begin{aligned} int_{0}^{pi}ke^{-ikrcos(theta)}sin{theta}d{theta}&=int_0^{pi}ke^{-ikrcos{(theta)}}d{(-cos{theta})}\ &=left.frac{e^{-ikcos{theta}}}{ir}right|^{pi}_{0}\ &=frac{2}{r}sin{kr} end{aligned} tag{1.23} ]

由于式（1.22）中被积表达式无(phi) 参与，可视为常数，其积分结果为(2pi)，式（1.22）可以写为如下形式：

[begin{aligned} G(mathbf{x},t;zeta,t_0)&=frac{1}{(2pi)^2cr}int_{0}^{infty}2sin(ck(t-t_0))sin(kr)d{k}\ &=frac{1}{(2pi)^2cr}int_{0}^{infty}cos[k(r-c(t-t_0))]-cos[k(r+c(t-t_0))]dk end{aligned} tag{1.24} ]

在式（1.24）中利用 (int_{0}^{infty}cos{kx}dx=pidelta{(x)})的性质，可以将（1.24）写成如下形式：

[G(mathbf{x},t;zeta,t_0)=frac{1}{(4pi)cr}{delta(r-c(t-t_0))-delta(r+c(t-t_0))} tag{1.25} ]

由于地震波实际传播过程，(t>t_0) 和空间任意一点到震源位置的距离(r>0)，根据脉冲函数(delta)的定义可知，(delta{(r+c(t-t_0))}=0)：

[G(mathbf{x},t;boldsymbol{zeta},t_0)=frac{1}{4pi cr}delta(r-c(t-t_0)) tag{1.26} ]

利用脉冲函数(delta)的性质

[delta(x-zeta)=delta(zeta-x) tag{1.27} ]

[delta{(cx)}=frac{1}{c}delta{(x)} tag{1.28} ]

可以将式（1.26）的形式写为如下形式：

[G(mathbf{x},t;boldsymbol{zeta},t_0)=frac{1}{4pi c^2r}delta(t-t_0-frac{r}{c}) tag{1.29} ]

由此，推导出波动方程在点脉冲震源的基本格林函数，其物理意义如下：位于(mathbf{zeta}) 时刻(t^{'}) 的单位点脉冲源，在位置(mathbf{x}) 时刻 (t=t_0+frac{r}{c}) 产生一个强度为(frac{1}{4pi c^2r})的脉冲响应。

（2） 利用格林函数求解任意形式源(f(mathbf{x},t))的波场

​ 通过上述获得的格林函数基本解，利用波场的线性叠加原理的来求解任意源(f(mathbf{x},t)) 产生的波场(varphi(mathbf{x},t))。其核心思想：将源分布函数(f(mathbf{x},t))看成无数个位于不同位置(mathbf{zeta}) 和不同时刻 (t_0) 的强度为(f(zeta,t_0)) 叠加而成。根据线性方程的叠加原理，总波场(u(mathbf{x},t)) 就是所有这些点源产生波的叠加（积分），在此假设(boldsymbol{zeta}in V_{zeta})和(t_0in[0,+infty)) 是连续的，其数学的表达形式如下：

[varphi(mathbf{x},t)=iiint_{V_{zeta}}int_{0}^{+infty}G(mathbf{x},t;boldsymbol{zeta},t_0)f(boldsymbol{zeta},t_0)dt_0d{V_zeta} tag{1.30} ]

将式（1.29）代入（1.30）可以得到：

[varphi(mathbf{x},t)=iiint_{V_{zeta}}int_{0}^{+infty}frac{1}{4pi c^2r}delta(t-t_0-frac{r}{c})f(boldsymbol{zeta},t_0)dt_0d{V_zeta} tag{1.31} ]

其中：(r=|boldsymbol{zeta}-mathbf{x}|)，首先考虑对(t_0) 进行积分。考虑脉冲函数的性质：

[int_{-infty}^{+infty}delta(s-t)f(s)ds=int_{-infty}^{+infty}delta(t-s)f(s)ds=f(t) tag{1.32} ]

上式可以化为如下形式得到一般震源分布的格林函数的一般形式解：

[varphi(mathbf{x},t)=frac{1}{4pi c^2}iiint_{V_{zeta}}frac{f(boldsymbol{zeta},t-frac{r}{c})}{|mathbf{x}-boldsymbol{zeta}|} d{V_{zeta}} tag{1.33} ]

式（1.33）中为由任意源(f(mathbf{x},t))在无限大自由空间中传播产生推迟势解(retarded solution)。其物理意义诠释：

1.推迟时间：解(varphi(mathbf{x},t))在位置(mathbf{x}) 和时间(t) 的值，取决于源函数分布距离及波速(c) 的值，即(t_{ret}=t-frac{|mathbf{x}-boldsymbol{zeta}|}{c})。 这体现了因果律：源的影响需要时间(frac{r}{c}) 才能以速度(c)的传播到观测点；

2.几何扩散：其解的形式是源函数在推迟时间的值(f(boldsymbol{zeta},t_{ret}))除以源点到观测点的距离(r=|boldsymbol{zeta}-mathbf{x}|) , 再对整个源区进行体积分。其中(frac{1}{r}) 项反映了波从点源向外传播时振幅的几何衰减，即球面扩散。