积分收敛性
通过输出核函数文件，观察源点和场点深度接近时，积分收敛性的变化。 具体积分表达式以及分类详见 计算动态格林函数。
核函数文件
该文件记录了波数积分过程中不同波数对应的核函数值。以下用法展示了如何输出该文件。
# -S<i1>,<i2>,...表示输出对应频率点索引值的核函数文件
grt greenfn -Mmilrow -D2/0 -N500/0.02 -OGRN -R5,8,10 -S50,100
输出的核函数文件会在 GRN_grtstats/milrow_{depsrc}_{deprcv}/ 路径下。
import numpy as np
import pygrt
modarr = np.loadtxt("milrow")
depsrc = 2.0
deprcv = 0.0
pymod = pygrt.PyModel1D(modarr, depsrc=depsrc, deprcv=deprcv)
# 通过statsfile参数自定义核函数文件的输出目录, statsidxs指定想输出的频率索引值
distarr = [5,8,10]
stgrnLst = pymod.compute_grn(
distarr=distarr, nt=500, dt=0.02,
statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100]
)
输出的核函数文件会在自定义路径下。
C和Python导出的核函数文件是一致的，底层调用的是相同的函数。文件名称格式为 K_{iw}_{freq}，其中 {iw} 表示频率索引值， {freq} 表示对应频率(Hz)。文件为自定义的二进制文件， 强烈建议使用Python进行读取及后续处理。这里还是给出两种读取方法。
ker2asc 模块可将单个核函数文件转为文本格式。
grt ker2asc GRN_grtstats/milrow_2_0/K_0050_5.00000e+00 > stats
输出的文件如下，
# k 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
4.18879020e-02 -2.02856187e-01 4.96041014e-01 -1.03425028e+02 -3.18035750e+01 -8.85165377e-02 -9.55284925e-02 4.93158286e+00 -1.21165651e+01 -2.26663577e+01 7.81815163e+00 9.45337168e-02 1.07242144e-01 2.26670306e+01 -7.81854269e+00 -3.25405088e+00 1.97453994e+00 -2.06838177e+02 -6.35936736e+01 -4.60765067e+01 -2.33936178e+02 -1.40769267e+00 8.13116327e-01 4.60781795e+01 2.33944595e+02 -9.49446169e-01 3.27485970e-01 3.95981907e-03 4.49214844e-03 9.49474359e-01 -3.27502350e-01
8.37758041e-02 -4.06881235e-01 9.91284229e-01 -1.03401670e+02 -3.18844819e+01 -1.77149204e-01 -1.91056283e-01 4.93534660e+00 -1.21176945e+01 -2.26687358e+01 7.80719981e+00 1.89079281e-01 2.14601619e-01 2.26714306e+01 -7.80876257e+00 -6.51103717e+00 3.94473178e+00 -2.06755819e+02 -6.37150286e+01 -4.59950051e+01 -2.33848344e+02 -2.81808839e+00 1.62610452e+00 4.60017376e+01 2.33882017e+02 -1.89909157e+00 6.54054442e-01 1.58402688e-02 1.79784232e-02 1.89931733e+00 -6.54185363e-01
1.25663706e-01 -6.13224641e-01 1.48487870e+00 -1.03361572e+02 -3.20176693e+01 -2.66015636e-01 -2.86601209e-01 4.94218337e+00 -1.21191468e+01 -2.26726441e+01 7.78893119e+00 2.83640484e-01 3.22190563e-01 2.26787183e+01 -7.79244160e+00 -9.77383475e+00 5.90611528e+00 -2.06616215e+02 -6.39138757e+01 -4.58603329e+01 -2.33701846e+02 -4.23407113e+00 2.43848237e+00 4.58756320e+01 2.33777632e+02 -2.84912849e+00 9.78785961e-01 3.56433145e-02 4.04876603e-02 2.84989179e+00 -9.79227091e-01
1.67551608e-01 -8.22988640e-01 1.97587080e+00 -1.03303260e+02 -3.22004756e+01 -3.55233324e-01 -3.82217751e-01 4.95285811e+00 -1.21201958e+01 -2.26779991e+01 7.76332549e+00 3.78206693e-01 4.30108797e-01 2.26888245e+01 -7.76955068e+00 -1.30451829e+01 7.85401462e+00 -2.06416413e+02 -6.41847550e+01 -4.56742513e+01 -2.33496556e+02 -5.65881724e+00 3.24903205e+00 4.57018091e+01 2.33631373e+02 -3.79973522e+00 1.30075767e+00 6.33691396e-02 7.20654206e-02 3.80154903e+00 -1.30180071e+00
2.09439510e-01 -1.03718423e+00 2.46316400e+00 -1.03225264e+02 -3.24288807e+01 -4.44911834e-01 -4.78013515e-01 4.96826413e+00 -1.21196804e+01 -2.26846821e+01 7.73036044e+00 4.72750387e-01 5.38433294e-01 2.27016507e+01 -7.74005327e+00 -1.63275746e+01 9.78345672e+00 -2.06153491e+02 -6.45194539e+01 -4.54392559e+01 -2.33232371e+02 -7.09576230e+00 4.05538089e+00 4.54829870e+01 2.33443271e+02 -4.75106871e+00 1.61904290e+00 9.90126096e-02 1.12769205e-01 4.75462261e+00 -1.62107297e+00
2.51327412e-01 -1.25666623e+00 2.94549729e+00 -1.03126785e+02 -3.26973550e+01 -5.35136602e-01 -5.74161488e-01 4.98921342e+00 -1.21159044e+01 -2.26925395e+01 7.69001804e+00 5.67231767e-01 6.47203866e-01 2.27170680e+01 -7.70390980e+00 -1.96231042e+01 1.16891316e+01 -2.05825887e+02 -6.49067297e+01 -4.51586056e+01 -2.32909294e+02 -8.54834304e+00 4.85357087e+00 4.52226353e+01 2.33213586e+02 -5.70325724e+00 1.93271234e+00 1.42560892e-01 1.62660073e-01 5.70942190e+00 -1.93620371e+00
2.93215314e-01 -1.48204805e+00 3.42145622e+00 -1.03008611e+02 -3.29988923e+01 -6.25943558e-01 -6.70903620e-01 5.01613464e+00 -1.21065803e+01 -2.27013835e+01 7.64229221e+00 6.61614510e-01 7.56406733e-01 2.27349160e+01 -7.66108337e+00 -2.29332759e+01 1.35654238e+01 -2.05435235e+02 -6.53324145e+01 -4.48363508e+01 -2.32527564e+02 -1.00194405e+01 5.63769904e+00 4.49249277e+01 2.32942956e+02 -6.65639329e+00 2.24083711e+00 1.93995506e-01 2.21790038e-01 6.66622554e+00 -2.24634697e+00
3.35103216e-01 -1.71360622e+00 3.88953365e+00 -1.02874248e+02 -3.33254195e+01 -7.17286702e-01 -7.68539699e-01 5.04868185e+00 -1.20888822e+01 -2.27109950e+01 7.58719728e+00 7.55896964e-01 8.65960920e-01 2.27550037e+01 -7.61154701e+00 -2.62587948e+01 1.54065499e+01 -2.04988585e+02 -6.57802801e+01 -4.44773486e+01 -2.32087839e+02 -1.15104865e+01 6.39978222e+00 4.45946555e+01 2.32632587e+02 -7.61052746e+00 2.54249421e+00 2.53303504e-01 2.90186290e-01 7.62527493e+00 -2.55065389e+00
3.76991118e-01 -1.95119908e+00 4.34825572e+00 -1.02731109e+02 -3.36688023e+01 -8.09005399e-01 -8.67395216e-01 5.08529186e+00 -1.20596866e+01 -2.27211300e+01 7.52477630e+00 8.50161369e-01 9.75716749e-01 2.27771106e+01 -7.55529200e+00 -2.95993908e+01 1.72068329e+01 -2.04500708e+02 -6.62340949e+01 -4.40872367e+01 -2.31591438e+02 -1.30202840e+01 7.13006622e+00 4.42370857e+01 2.32284442e+02 -8.56566422e+00 2.83677383e+00 3.20503286e-01 3.67836549e-01 8.58676841e+00 -2.84827798e+00
...
后续你可以选择习惯的方式读取和处理。
# 可使用通配符简化输入，因为对应索引值下只会有一个文件
# 返回的是自定义类型的numpy数组
statsdata = pygrt.utils.read_statsfile(f"pygrtstats_{depsrc}_{deprcv}/K_0050_*")
print(statsdata.dtype)
# [('k', '<f8'), ('EX_q', '<c16'), ('EX_w', '<c16'), ('VF_q', '<c16'), ('VF_w', '<c16'), ('HF_q', '<c16'), ('HF_w', '<c16'), ('HF_v', '<c16'), ('DD_q', '<c16'), ('DD_w', '<c16'), ('DS_q', '<c16'), ('DS_w', '<c16'), ('DS_v', '<c16'), ('SS_q', '<c16'), ('SS_w', '<c16'), ('SS_v', '<c16')]
其中除了波数 k 外，每条结果的命名格式均为 {srcType}_{q/w/v}，与 计算动态格林函数 部分介绍的积分公式中的核函数 \(q_m, w_m, v_m\) 保持一致。
备注
核函数文件中记录的值非理论核函数值。对于动态解，还需乘 \(\left(-\dfrac{1}{4\pi\rho\omega^2}\right)\)。
可视化
以下将使用Python进行图件绘制。 在Python函数中指定震源类型、阶数、积分类型，可自动绘制核函数、被积函数和积分值随波数的变化，其中积分类型对应 计算动态格林函数 部分介绍的4种类型。
ir = 2
dist=distarr[ir]
srctype="SS"
ptype="0"
fig, ax = pygrt.utils.plot_statsdata(statsdata, dist=dist, srctype=srctype, ptype=ptype)
fig.savefig(f"{srctype}_{ptype}.svg", bbox_inches='tight')
可以通过指定 RorI=False 参数来指定绘制虚部，传入 RorI=2 表示实虚部都绘制。
ir = 2
dist=distarr[ir]
srctype="SS"
ptype="0"
fig, ax = pygrt.utils.plot_statsdata(statsdata, dist=dist, srctype=srctype, ptype=ptype, RorI=2)
fig.savefig(f"{srctype}_{ptype}_RI.svg", bbox_inches='tight')
从图中可以看到，当波数增加到一定范围以上，积分收敛良好，无振荡现象。
当震源和场点深度接近/相等时，积分收敛速度变慢
假设其它参数不变，我们调整 震源深度为0.1km，再计算格林函数，读取核函数文件，绘制图像。
# 使用 -Cn 禁用任何收敛算法
grt greenfn -Mmilrow -D0/0 -N500/0.02 -OGRN -R5,8,10 -Cn -S50,100
grt ker2asc GRN_grtstats/milrow_0_0/K_0050_5.00000e+00 > stats
# 绘制图像部分见Python
depsrc = 0.0
deprcv = 0.0
pymod = pygrt.PyModel1D(modarr, depsrc=depsrc, deprcv=deprcv)
stgrnLst = pymod.compute_grn(
distarr=distarr, nt=500, dt=0.02, converg_method='none',
statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100]
)
statsdata = pygrt.utils.read_statsfile(f"pygrtstats_{depsrc}_{deprcv}/K_0050_*")
dist=10
srctype="SS"
ptype="0"
fig, ax = pygrt.utils.plot_statsdata(statsdata, dist=dist, srctype=srctype, ptype=ptype, RorI=2)
fig.savefig(f"{srctype}_{ptype}_{depsrc}_RI.svg", bbox_inches='tight')
从图中可以清晰地看到，相比震源深度2km时，积分收敛速度明显变慢，积分值振荡，这要求增加波数积分上限，但这必然降低计算效率。
你可以尝试 当震源和场点位于同一深度，收敛振荡现象也很明显。要注意的是这和是否在地表无关，即使震源和场点都在地下，深度接近时积分收敛也会变慢。