核函数 \(f-v\) 谱图
在 积分收敛性 部分介绍了程序可以导出积分过程中不同频率的核函数值， 这使得我们可以绘制核函数的 \(f-v\) 谱图 [1] ，以观察其中的频散特征 。
在本页文档的旧版本中（ 核函数 f-v 谱图（旧版） ），我们使用了格林函数计算过程中导出的核函数文件， 其中每个采样点是在波数上等距采样，需要经过插值才能转到相速度，且定义不同频率、不同相速度的采样点不方便。
为了更方便地导出核函数文件，在 0.14.0 版本起，我新增了模块 kernel 直接面向这一任务。
以薄层模型为例，
0.01 1.5 0.18 1.78 1e4 1e4
0.01 1.7 0.35 1.85 1e4 1e4
0.02 1.6 0.25 1.80 1e4 1e4
0.00 2.0 0.6 1.94 1e4 1e4
核函数文件
使用 grt 程序调用 kernel 计算核函数，
# +w 控制虚频率，从而调整谱图中高亮部分的粗细
grt kernel -Mmod1 -D0.03/0 -F0/25/0.1+w2 -C0.1/0.6/5e-3 -OKERNEL
得到的每个频率的核函数文件名开头为 C ，表示采样位置直接记录的相速度，
且 积分收敛性 介绍的核函数文件的读取方法依然适用，
只是对应波数的列会改为等距的相速度，列名用 “c” 表示。
备注
注意其中使用了虚频率，这使得 \(f-v\) 谱图中的高亮区域变宽，方便观察。
绘制核函数 \(f-v\) 谱图
以下将使用 Python 进行图件绘制。
import numpy as np
import matplotlib.pyplot as plt
from typing import Union
import pygrt
# 读取所有频率的核函数
# 不指定ktypes，默认返回全部核函数，均以2D数组的形式保存，shape=(nfreqs, nvels)
kerDct = pygrt.utils.read_kernels_freqs("KERNEL/mod1_0.03_0")
print(kerDct.keys())
# dict_keys(['_vels', '_freqs', 'EX_q', 'EX_w', 'VF_q', 'VF_w', 'HF_q', 'HF_w', 'HF_v', 'DD_q', 'DD_w', 'DS_q', 'DS_w', 'DS_v', 'SS_q', 'SS_w', 'SS_v'])
# 绘制图像
def plot_kernel(kerDct:dict, RorI:bool, out:Union[str,None]=None):
funcRorI = np.real if RorI else np.imag
ktypes = []
for key in kerDct:
if key[0] == '_':
continue
ktypes.append(key)
srctypes = ['EX', 'VF', 'HF', 'DD', 'DS', 'SS']
vels = kerDct['_vels']
freqs = kerDct['_freqs']
fig = plt.figure(figsize=(12, 3*len(srctypes)))
gs = fig.add_gridspec(len(srctypes), 3)
qwvLst = ['q', 'w', 'v']
for ik, key in enumerate(ktypes):
srctype, qwv = key.split("_")
ax = fig.add_subplot(gs[srctypes.index(srctype), qwvLst.index(qwv)])
# 对不同速度间取归一化
data = kerDct[key].copy()
data[...] = data/np.max(np.abs(data), axis=1)[:,None]
pcm = ax.pcolormesh(freqs, vels, np.abs(funcRorI(data)).T, vmin=0, vmax=1, shading='nearest', rasterized=True)
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Velocity (km/s)")
ax.set_title(key)
fig.colorbar(pcm, ax=ax)
if RorI:
fig.suptitle("Real parts of Kernels", fontsize=20, x=0.5, y=0.99)
else:
fig.suptitle("Imag parts of Kernels", fontsize=20, x=0.5, y=0.99)
if out is not None:
fig.savefig(out, bbox_inches='tight')
plot_kernel(kerDct, False, "imag.svg")
plot_kernel(kerDct, True, "real.svg")
虚部
实部
从图上可看到， 核函数 \(f-v\) 谱图的虚部峰值正是频散曲线的位置 ，尽管其中不同震源的核函数幅值不同，但 \(q_m，w_m\) 呈现的频散特征一致，这对应的是 Rayleigh波频散 ，而 \(v_m\) 对应的是 Love波频散 。
用户还可调整震源场点深度观察不同部分的能量变化，还可调整相速度范围观察泄露模 （leaky-mode）。 程序包为理论研究提供了便捷工具。