diff --git a/docs/source/Tutorial/index.rst b/docs/source/Tutorial/index.rst index bd23bb88..5544c029 100644 --- a/docs/source/Tutorial/index.rst +++ b/docs/source/Tutorial/index.rst @@ -83,7 +83,7 @@ dynamic/strain_stress static/static_strain_stress integ_converg/integ_converg - + kernel/kernel diff --git a/docs/source/Tutorial/integ_converg/integ_converg.rst b/docs/source/Tutorial/integ_converg/integ_converg.rst index 3d5c6884..84a88f7a 100644 --- a/docs/source/Tutorial/integ_converg/integ_converg.rst +++ b/docs/source/Tutorial/integ_converg/integ_converg.rst @@ -1,3 +1,5 @@ +.. _integ_converg_rst: + 积分收敛性 =================== @@ -199,13 +201,18 @@ C和Python导出的核函数文件是一致的，底层调用的是相同的函 :start-after: BEGIN SGRN :end-before: END SGRN +---------------------------------- + -+ 只使用离散波数积分 ++ **只使用离散波数积分** .. image:: run/DC_20_0.1_static.png :align: center -+ 使用峰谷平均法 +---------------------------------- + + ++ **使用峰谷平均法** .. image:: run/DC_20_0.1_ptam_static.png :align: center diff --git a/docs/source/Tutorial/kernel/kernel.rst b/docs/source/Tutorial/kernel/kernel.rst new file mode 100644 index 00000000..de9184a2 --- /dev/null +++ b/docs/source/Tutorial/kernel/kernel.rst @@ -0,0 +1,82 @@ +核函数频谱图 +============= + +:Author: Zhu Dengda +:Email: zhudengda@mail.iggcas.ac.cn + + +----------------------------------------------------------- + +在 :ref:`integ_converg_rst` 部分介绍了程序可以导出积分过程中不同频率的核函数值，这使得我们可以绘制核函数的频谱图 [#]_ ，以观察其中的频散特征 。 + +以薄层模型为例， + +.. literalinclude:: run/mod1 + :language: text + +.. [#] 感谢席超强博士 `@xichaoqiang `_ 的建议。 + + +生成核函数文件 +---------------------- + +.. note:: + + 与之前的示例相比，这里添加了控制波数积分的参数，以显式控制核函数的采样间隔和范围。但这并非是计算格林函数的最佳参数，故此时格林函数计算结果很可能是不收敛的。 + + +.. tabs:: + + .. group-tab:: C + + .. literalinclude:: run/run.sh + :language: bash + :start-after: BEGIN GRN + :end-before: END GRN + + + .. group-tab:: Python + + .. literalinclude:: run/run.py + :language: python + :start-after: BEGIN GRN + :end-before: END GRN + + +读取，采样，获得核函数频谱 +---------------------------------- +Python端提供了 :py:func:`pygrt.utils.read_kernels_freqs` 函数来完成所有频率的核函数读取以及从波数 ``k`` 插值到特定速度点 ``v`` 的工作 ( :math:`k=\omega/v` )，最后获得保存有核函数频谱的字典。 + +.. literalinclude:: run/run.py + :language: python + :start-after: BEGIN read + :end-before: END read + +绘制核函数频谱图 +------------------------ +以下将使用Python进行图件绘制。 + +.. literalinclude:: run/run.py + :language: python + :start-after: BEGIN plot + :end-before: END plot + +---------------------------------- + ++ **虚部** + +.. image:: run/imag.png + :align: center + +---------------------------------- + ++ **实部** + +.. image:: run/real.png + :align: center + + +---------------------------------- + + +从图上可看到， **核函数频谱的虚部峰值正是频散曲线的位置** ，其中 :math:`q_m，w_m` 呈现的频散特征一致，只是不同震源的核函数幅值不同，这对应的是 **Rayleigh波频散** ，而 :math:`v_m` 对应的是 **Love波频散**。 \ No newline at end of file diff --git a/docs/source/Tutorial/kernel/run/mod1 b/docs/source/Tutorial/kernel/run/mod1 new file mode 100644 index 00000000..917e6089 --- /dev/null +++ b/docs/source/Tutorial/kernel/run/mod1 @@ -0,0 +1,4 @@ +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 \ No newline at end of file diff --git a/docs/source/Tutorial/kernel/run/run.py b/docs/source/Tutorial/kernel/run/run.py new file mode 100644 index 00000000..eb40688e --- /dev/null +++ b/docs/source/Tutorial/kernel/run/run.py @@ -0,0 +1,82 @@ +# ----------------------------------------------------------------- +# BEGIN GRN +import numpy as np +import matplotlib.pyplot as plt +from typing import Union +import pygrt + +modarr = np.loadtxt("mod1") + +pymod = pygrt.PyModel1D(modarr, depsrc=0.01, deprcv=0.0) + +# 不指定statsidx，默认输出全部频率点的积分过程文件 +# vmin_ref 显式给定参考速度（用于定义波数积分上限），避免使用PTAM +# Length 给定波数积分间隔dk +_ = pymod.compute_grn(distarr=[1], nt=500, dt=0.02, vmin_ref=0.1, Length=20, statsfile="pygrtstats") +# END GRN +# ----------------------------------------------------------------- + +# ----------------------------------------------------------------- +# BEGIN read +# 指定待采样的速度数组 +vels = np.arange(0.1, 0.6, 0.001) + +# 读取所有频率的核函数，并插值到vels +# 不指定ktypes，默认返回全部核函数，均以2D数组的形式保存，shape=(nfreqs, nvels) +kerDct = pygrt.utils.read_kernels_freqs("pygrtstats", vels) +print(kerDct.keys()) +# dict_keys(['_vels', '_freqs', 'EXP_q0', 'EXP_w0', 'VF_q0', 'VF_w0', 'HF_q1', 'HF_w1', 'HF_v1', 'DC_q0', 'DC_w0', 'DC_q1', 'DC_w1', 'DC_v1', 'DC_q2', 'DC_w2', 'DC_v2']) +# END read +# ----------------------------------------------------------------- + + +# ----------------------------------------------------------------- +# BEGIN plot +# 绘制图像 +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 = ['EXP0', 'VF0', 'HF1', 'DC0', 'DC1', 'DC2'] + + 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 = key.split("_")[0]+key.split("_")[1][1:] + qwv = key.split("_")[1][0] + + 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') + 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.tight_layout() + fig.savefig(out, dpi=100) + + +plot_kernel(kerDct, False, "imag.png") +plot_kernel(kerDct, True, "real.png") +# END plot +# ----------------------------------------------------------------- diff --git a/docs/source/Tutorial/kernel/run/run.sh b/docs/source/Tutorial/kernel/run/run.sh new file mode 100755 index 00000000..b9fc2817 --- /dev/null +++ b/docs/source/Tutorial/kernel/run/run.sh @@ -0,0 +1,10 @@ +#!/bin/bash + +# ----------------------------------------------------------------- +# BEGIN GRN +# -S-1 表示输出所有频率点的核函数 +# -V0.1 显式给定参考速度（用于定义波数积分上限），避免使用PTAM +# -L20 定义波数积分间隔dk +grt -Mmod1 -D0/0 -N500/0.02/0.8 -OGRN -R1 -V0.1 -S-1 -L20 +# END GRN +# ----------------------------------------------------------------- diff --git a/docs/source/Tutorial/run_all.sh b/docs/source/Tutorial/run_all.sh index 5ae355d8..e9b30fb6 100755 --- a/docs/source/Tutorial/run_all.sh +++ b/docs/source/Tutorial/run_all.sh @@ -42,4 +42,12 @@ if [[ $1 == '5' || $1 == '' ]]; then ./run.sh python run.py cd - +fi + +if [[ $1 == '6' || $1 == '' ]]; then + cd kernel/run + chmod +x *.sh + ./run.sh + python run.py + cd - fi \ No newline at end of file diff --git a/pygrt/utils.py b/pygrt/utils.py index b3003daa..3404acab 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -23,6 +23,7 @@ from scipy.signal import oaconvolve from scipy.fft import rfft, irfft from scipy.special import jv +from scipy.interpolate import interpn import math import os import glob @@ -48,6 +49,7 @@ "stream_diff", "stream_write_sac", + "read_kernels_freqs", "read_statsfile", "read_statsfile_ptam", "plot_statsdata", @@ -1250,6 +1252,75 @@ def read_statsfile(statsfile:str): return data +def read_kernels_freqs(statsdir:str, vels:np.ndarray, ktypes:Union[List[str],None]=None): + r""" + 读取statsdir目录下所有频率（除了零频）的积分过程文件（K_开头），再将核函数线性插值到指定的速度数组vels + + :param statsdir: 存储积分过程文件的目录 + :param vels: 待插值的速度数组(km/s)，必须正序 + :param ktype: 指定返回一系列的核函数名称，如EXP_q0，DC_w2等，默认返回全部 + + :return: + - **kerDct** - 字典格式的核函数插值结果 + """ + + if not np.all(np.diff(vels) > 0): + raise ValueError("vels must be in ascending order.") + + KLst = np.array(glob.glob(os.path.join(statsdir, "K_*"))) + freqs = np.array([float(s.split("_")[-1]) for s in KLst]) + # 根据freqs排序 + _idx = np.argsort(freqs) + freqs[:] = freqs[_idx] + KLst[:] = KLst[_idx] + del _idx + + # 去除零频 + if freqs[0] == 0.0: + freqs = freqs[1:] + KLst = KLst[1:] + + kerDct = {} + kerDct['_vels'] = vels.copy() + kerDct['_freqs'] = freqs.copy() + + for i in range(len(freqs)): + Kpath = KLst[i] + freq = freqs[i] + w = 2*np.pi*freq + + data = read_statsfile(Kpath) + v = w/data['k'] + + # 检查v范围 + v1 = np.min(v) + v2 = np.max(v) + if v1 > vels[0] or v2 < vels[-1]: + raise ValueError(f"In freq={freq:.5e}, minV={v1:.5e}, maxV={v2:.5e}, insufficient wavenumber samples" + " to interpolate on vels.") + + for key in data.dtype.names: + if key == 'k': + continue + if (ktypes is not None) and (key not in ktypes): + continue + + if key not in kerDct.keys(): + kerDct[key] = [] + + # 如果越界会报错 + F = interpn((v,), data[key], vels) + kerDct[key].append(F) + + # 将每个核函数结果拼成2D数组 + for key in kerDct.keys(): + if key[0] == '_': + continue + kerDct[key] = np.vstack(kerDct[key]) + + return kerDct + + def read_statsfile_ptam(statsfile:str): ''' 读取单个频率下峰谷平均法的记录文件