✍️ 朱邓达 • 📅 2025-04-19

积分收敛性

通过输出核函数文件，观察源点和场点深度接近时，积分收敛性的变化。 具体积分表达式以及分类详见 计算动态格林函数。

核函数文件

该文件记录了波数积分过程中不同波数对应的核函数值。以下用法展示了如何输出该文件。

# -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
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
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
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
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
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
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
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
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
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时，积分收敛速度明显变慢，积分值振荡，这要求增加波数积分上限，但这必然降低计算效率。

你可以尝试 当震源和场点位于同一深度，收敛振荡现象也很明显。要注意的是这和是否在地表无关，即使震源和场点都在地下，深度接近时积分收敛也会变慢。