✍️ 朱邓达  •  📅 2025-04-17  •  ⚡ 2025-09-22

计算动态格林函数

Python中计算动态格林函数的主函数为 compute_grn() ，C模块为 greenfn。

核心计算逻辑来自 初稿 ，具体代码可见与C API中对应同名 *.c 文件，其中计算格林函数频谱的主函数为 grn.c 里的 integ_grn_spec()。输出结果的坐标系见下图。

../../_images/coord.svg

警告

推导过程中定义的Z分量与C/Python输出的Z分量方向相反。因为从Z轴向下为正的坐标系对于 半无限水平分层均匀介质 来讲更便于公式推导；而在结果输出后，会对Z分量取相反数，因为Z轴向上和目前实采数据使用的Z轴方向一致，这更便于后续数据处理。

除非特别指明（例如描述震源机制），且一般使用 PyGRT 时也不会涉及底层程序内部， 在实际使用过程中提及的Z轴均以向上为正 ，仅在讨论代码内部细节和公式推导时对此加以区分。

示例程序

假设在 milrow 模型中，震源深度2km，接收点位于地表，震中距为5km，8km和10km，计算格林函数，要求采样点数500个点，采样间隔0.02s。

# 不同震源深度、接收点深度和震中距的格林函数会在 
# GRN/milrow_{depsrc}_{deprcv}_{dist}/ 路径下，使用SAC格式保存。 
grt greenfn -Mmilrow -D2/0 -N500/0.02 -OGRN -R5,8,10

不同震源深度、接收点深度和震中距的格林函数会在 GRN/milrow_{depsrc}_{deprcv}_{dist}/ 路径下，使用SAC格式保存。

一些基本信息（包括源点和场点的物性参数）保存在SAC头段变量中，其中 t0 和 t1 分别代表初至P波和初至S波的到时。在不设定其它参数时，程序中使用0（发震时刻）作为参考时间，故其等价于走时。

  • travt 模块可以显式地再计算初至波走时，

grt travt -Mmilrow -D2/0 -R5,8,10

走时结果会输出到终端，

------------------------------------------------
 Distance(km)     Tp(secs)         Ts(secs)     
 5                1.270            2.395          
 8                1.881            3.465          
 10               2.256            4.120          
------------------------------------------------

  • 如果你没有安装SAC软件，可以使用Python的ObsPy库读取生成的SAC数据，或者使用 sac2asc 模块临时将SAC格式文件转为如下的文本文件：

grt sac2asc GRN/milrow_2_0_10/HFZ.sac > HFZ

输出的文本文件如下，两列分别为时间点和幅值。这种输出仅保留波形信息，缺失SAC文件中的头段变量。

0.0000000e+00  -7.0514270e-06
2.0000000e-02  8.7115750e-06
3.9999999e-02  -7.6045740e-06
5.9999999e-02  9.3538965e-06
7.9999998e-02  -7.9860201e-06
9.9999994e-02  9.7270295e-06
1.2000000e-01  -8.4296462e-06
1.4000000e-01  9.9532535e-06
1.6000000e-01  -9.1152588e-06
1.7999999e-01  1.0513032e-05
...

当时窗长度 nt*dt 太小“包不住”有效信号，或时窗长度足够但时延不合适，输出的波形会发生混叠， 此时需调整相关参数。

备注

格林函数计算结果的单位：

  • 爆炸源： \(10^{-20} \, \frac{\text{cm}}{\text{dyne} \cdot \text{cm}}\)

  • 单力源： \(10^{-15} \, \frac{\text{cm}}{\text{dyne}}\)

  • 剪切源： \(10^{-20} \, \frac{\text{cm}}{\text{dyne} \cdot \text{cm}}\)

离散波数积分

格林函数频谱的计算本质转化为求以下积分:

\[P_m(\omega) = \int_0^\infty F_m(k, \omega)J_m(kr)kdk\]

其中 \(F_m(k,\omega)\) 称为核函数，它是和介质属性相关的量，与震中距无关。我们可以使用离散波数积分法 (Bouchon, 1981) 将上式积分转变为求和:

\[P_m(\omega) = \Delta k \sum_{j=0}^{\infty} F_m(k_j,\omega)J_m(k_j r)k_j\]

其中 \(\Delta k = 2\pi/L, k_j=j\Delta k\)，\(L\) 为特征长度，要求满足:

\[\begin{split}\left\{ \begin{aligned} & r < L/2 \\ & (L-r)^2 + z_s^2 > (\alpha T)^2 \end{aligned} \right.\end{split}\]

其中 \(\alpha\) 为参考P波速度， \(T\) 为所需计算的理论地震图的总时间长度。常见的保守经验值为 \(L=20r\) ，但也应依具体情况而定 。为了避开附加源以及奇点的影响，(Bouchon, 1981) 在频率上添加微小虚部，具体推导过程详见 (Bouchon, 1981) 和 (张海明, 2021)。

积分形式分类

通过在面谐矢量坐标系中建立波函数进行公式推导，最终格林函数的三分量频谱 \(W_m(\omega), Q_m(\omega), V_m(\omega)\) (分别为垂向，径向，切向)可以表达为:

\[\begin{split}\left\{ \begin{aligned} W_m(\omega) &= \int_0^\infty w_m J_m(kr)kdk \\ Q_m(\omega) &= \int_0^\infty (q_m J_m^{\prime}(kr) - v_m \frac{m}{kr} J_m(kr)) kdk \\ V_m(\omega) &= \int_0^\infty (q_m \frac{m}{kr} J_m(kr) - v_m J_m^{\prime}(kr)) kdk \end{aligned} \right.\end{split}\]

备注

初次推导该公式可能会对虚数 \(i\) 及公式中的正负号感到疑惑，但其实这里的设计是将虚数 \(i\) 和方向因子 \(e^{im\theta}\) 合并，所以在后续合成理论地震图时你会发现，\(m=0,1,2\) 阶的 \(W_m, Q_m\) 的方向因子对 \((m\theta)\) 的偏导就是 \(V_m\) 的方向因子。

公式来自 初稿 (5.6.22)式，其中阶数 \(m=0,1,2\)。核函数 \(q_m,w_m,v_m\) 根据广义反射透射系数矩阵法(GRTM)求得，当震源比场点深时，有如下公式（震源浅于场点时有类似公式，这里不再展示）。

\[\begin{split}\begin{aligned} \begin{bmatrix} q_m \\ w_m \end{bmatrix} &= \mathbf{R}_{EV}(z_R) \left(\mathbf{I} - \mathbf{R}_D^{RS}\mathbf{R}_U^{FR}\right)^{-1} \mathbf{T}_U^{RS} \left(\mathbf{I} - \mathbf{R}_D^{SL}\mathbf{R}_U^{FS}\right)^{-1} \left[ \mathbf{R}_D^{SL} \begin{pmatrix} P_m^+ \\ SV_m^+ \end{pmatrix} + \begin{pmatrix} P_m^- \\ SV_m^- \end{pmatrix} \right] \\ v_m &= R_{EV,L}(z_R) \left(I - R_{D,L}^{RS}R_{U,L}^{FR}\right)^{-1} T_{U,L}^{RS} \left(I - R_{D,L}^{SL}R_{U,L}^{FS}\right)^{-1} \left[ R_{D,L}^{SL} SH_m^+ + SH_m^- \right] \end{aligned}\end{split}\]

为了方便程序实现，根据积分形式，我们对待求积分进行如下分类，其中每一阶都分为4类( \(p=0,1,2,3\) )，除了0阶只需两类，此时 \(v_0=0\) :

\[\begin{split}\left\{ \begin{aligned} p=0 & \rightarrow - \int q_0(k, \omega) J_1(kr)kdk \\ p=2 & \rightarrow \int w_0(k, \omega) J_0(kr)kdk \end{aligned} \right.\end{split}\]

  • \(m=1,2\)

\[\begin{split}\left\{ \begin{aligned} p=0 & \rightarrow \int q_m(k, \omega) J_{m-1}(kr)kdk \\ p=1 & \rightarrow - \int (q_m(k, \omega) + v_m(k, \omega)) \frac{m}{kr} J_m(kr)kdk \\ p=2 & \rightarrow \int w_m(k, \omega) J_m(kr)kdk \\ p=3 & \rightarrow - \int v_m(k, \omega) J_{m-1}(kr)kdk \end{aligned} \right.\end{split}\]

以上每个积分都形成 \(\int_0^\infty F_m(k, \omega)J_m(kr)kdk\) 的形式，便可逐个使用离散波数积分（或Filon积分、峰谷平均法等）求解每个积分。

格林函数分类

程序会输出15个格林函数（也可以选择输出哪些震源），但并不是每个震源类型对应的每一阶每种积分类型都存在。以下为15个格林函数定义的名称，以及对应上述的阶数以及积分类型:

名称

格林函数类型

对应阶数

对应积分类型

EXZ

爆炸源Z分量

\(m=0\)

\(p=2\)

EXR

爆炸源R分量

\(m=0\)

\(p=0\)

VFZ

垂直向下力源Z分量

\(m=0\)

\(p=2\)

VFR

垂直向下力源R分量

\(m=0\)

\(p=0\)

HFZ

水平力源Z分量

\(m=1\)

\(p=2\)

HFR

水平力源R分量

\(m=1\)

\((p=0)+(p=1)\)

HFT

水平力源T分量

\(m=1\)

\(-(p=1)+(p=3)\)

DDZ

倾角45度倾滑Z分量

\(m=0\)

\(p=2\)

DDR

倾角45度倾滑R分量

\(m=0\)

\(p=0\)

DSZ

倾角90度倾滑Z分量

\(m=1\)

\(p=2\)

DSR

倾角90度倾滑R分量

\(m=1\)

\((p=0)+(p=1)\)

DST

倾角90度倾滑T分量

\(m=1\)

\(-(p=1)+(p=3)\)

SSZ

倾角90度走滑Z分量

\(m=2\)

\(p=2\)

SSR

倾角90度走滑R分量

\(m=2\)

\((p=0)+(p=1)\)

SST

倾角90度走滑T分量

\(m=2\)

\(-(p=1)+(p=3)\)