✍️ 朱邓达 • 🗓️ 2025-04-19 (创建时间)

峰谷平均法（Peak-Trough Averaging Method, PTAM）

具体原理很简洁易懂，从以下图像中你也能了解大概。详见 (Zhang et al., 2003; 张海明, 2021) 。

具体流程为，程序中在进行完离散波数积分后，继续增加k（同一频率下所有震中距复用DWM的统一dk），使用PTAM寻找足够数量的波峰和波谷（内部设定为36个），再对这些波峰波谷取缩减序列 \(M_i \leftarrow 0.5\times(M_i + M_{i+1})\) ，得到估计的积分收敛值。取缩减序列的C代码如下，

for(int n=n1; n>1; --n){
    for(int i=0; i<n-1; ++i){
        arr[i] = 0.5*(arr[i] + arr[i+1]);
    }
}

以下通过显式指定相关参数来使用峰谷平均法进行收敛。

# 使用 -Cp 指定使用 PTAM 进行收敛
grt greenfn -Mmilrow -D0/0 -N500/0.02 -OGRN -R5,8,10 -Cp -K+k2+f+s1.2 -S50,100
# 绘制图像部分见Python

输出的核函数文件会在 GRN_grtstats/milrow_{depsrc}_{deprcv}/ 路径下。

import numpy as np
import pygrt 
from pygrt.cli import format_float

depsrc = 0.0
deprcv = 0.0
pymod = pygrt.PyModel1D(grn="GRN", modelpath="milrow")

dists = [5,8,10]
# 设置 converg_method='PTAM' 进行收敛
pymod.greenfn(
    depsrc=depsrc, deprcv=deprcv,
    dists=dists, nt=500, dt=0.02, converg_method='PTAM', k0=2, ampk=1.2, use_kmax_ref=True,
    statsidxs=[50,100],
)

输出的核函数文件会在 GRN_grtstats/milrow_{depsrc}_{deprcv}/ 路径下。

在 K_{iw}_{freq} 文件同级目录下，程序把 PTAM过程中的核函数以及积分峰谷位置分为两个文件 保存在 PTAM_{ir}_{dist}/ 目录下（ {ir} 为震中距索引， {dist} 为震中距），其中 PTAM_{ir}_{dist}/K_{iw}_{freq} 为核函数文件（格式不变）， PTAM_{ir}_{dist}/PTAM_{iw}_{freq} 为峰谷文件，其中记录积分值的峰谷。

备注

ker2asc 模块也支持将 PTAM_{ir}_{dist}/PTAM_{iw}_{freq} 文件转为文本格式，

grt ker2asc GRN_grtstats/milrow_0_0/PTAM_0002_1.00000e+01/PTAM_0050_5.00000e+00 > ptam_stats

输出的文件如下，

#       sum_EX_0_k                          sum_EX_0        sum_EX_2_k                          sum_EX_2        sum_VF_0_k                          sum_VF_0        sum_VF_2_k                          sum_VF_2        sum_HF_0_k                          sum_HF_0        sum_HF_1_k                          sum_HF_1        sum_HF_2_k                          sum_HF_2        sum_HF_3_k                          sum_HF_3        sum_DD_0_k                          sum_DD_0        sum_DD_2_k                          sum_DD_2        sum_DS_0_k                          sum_DS_0        sum_DS_1_k                          sum_DS_1        sum_DS_2_k                          sum_DS_2        sum_DS_3_k                          sum_DS_3        sum_SS_0_k                          sum_SS_0        sum_SS_1_k                          sum_SS_1        sum_SS_2_k                          sum_SS_2        sum_SS_3_k                          sum_SS_3
69735334e+01    4.21892682e+02 -6.56930483e+03    6.68163835e+01    1.42862473e+04  7.52643521e+02    6.69735406e+01   -9.16220458e+01  1.37280443e+03    6.68163900e+01   -1.91333772e+03 -3.80404893e+01    6.68163896e+01   -6.09406253e+02 -3.76789375e+02    6.69735461e+01   -1.00020779e+01  1.00663153e+01    6.69735406e+01    9.16220458e+01 -1.37280443e+03    6.68163898e+01    1.45529681e+03  8.56571876e+02    6.69735334e+01   -1.68757073e+03  2.62772193e+04    6.68163789e+01   -5.20668710e+04 -3.09182918e+03    6.67274280e+01    0.00000000e+00  0.00000000e+00    6.67274280e+01    0.00000000e+00  0.00000000e+00    6.67274280e+01    0.00000000e+00  0.00000000e+00    6.67274280e+01    0.00000000e+00  0.00000000e+00    6.69735334e+01    8.43785365e+02 -1.31386097e+04    6.68160944e+01   -4.16930671e+02 -2.42645907e+02    6.68160889e+01    2.69000079e+04  1.25777733e+03    6.69735336e+01   -1.86041359e+04  2.77878622e+04
72875193e+01    2.75012932e+03 -6.60766958e+03    6.71307183e+01    1.19585872e+04  7.91065210e+02    6.72875227e+01   -1.10460065e+02  1.37314075e+03    6.71307214e+01   -1.84179266e+03 -3.92556369e+01    6.71307212e+01   -5.39972622e+02 -3.77933678e+02    6.72875254e+01   -9.94871418e+00  1.00654080e+01    6.72875227e+01    1.10460065e+02 -1.37314075e+03    6.71307213e+01    1.56061705e+03  8.54817236e+02    6.72875193e+01   -1.10005173e+04  2.64306783e+04    6.71307162e+01   -5.29295871e+04 -3.08273181e+03    6.68112038e+01    0.00000000e+00  0.00000000e+00    6.68112038e+01    0.00000000e+00  0.00000000e+00    6.68112038e+01    0.00000000e+00  0.00000000e+00    6.68112038e+01    0.00000000e+00  0.00000000e+00    6.72875193e+01    5.50025863e+03 -1.32153392e+04    6.71304204e+01   -4.24108032e+02 -2.42523839e+02    6.71304179e+01    2.56358048e+04  1.28035939e+03    6.72875194e+01   -1.15413441e+04  2.76702193e+04
76018525e+01    4.17097502e+02 -6.56923730e+03    6.74447019e+01    1.42910140e+04  7.52577613e+02    6.76018596e+01   -9.16861532e+01  1.37280592e+03    6.74447083e+01   -1.91312973e+03 -3.80447448e+01    6.74447079e+01   -6.09224679e+02 -3.76792712e+02    6.76018651e+01   -1.00016748e+01  1.00663080e+01    6.76018596e+01    9.16861532e+01 -1.37280592e+03    6.74447081e+01    1.45558320e+03  8.56566407e+02    6.76018525e+01   -1.66839001e+03  2.62769492e+04    6.74446974e+01   -5.20621078e+04 -3.09194685e+03    6.68949796e+01    0.00000000e+00  0.00000000e+00    6.68949796e+01    0.00000000e+00  0.00000000e+00    6.68949796e+01    0.00000000e+00  0.00000000e+00    6.68949796e+01    0.00000000e+00  0.00000000e+00    6.76018525e+01    8.34195004e+02 -1.31384746e+04    6.74444155e+01   -4.16951636e+02 -2.42645481e+02    6.74444100e+01    2.69015978e+04  1.25777262e+03    6.76018526e+01   -1.86179482e+04  2.77880456e+04
79158384e+01    2.75491682e+03 -6.60773711e+03    6.77590367e+01    1.19538276e+04  7.91131169e+02    6.79158418e+01   -1.10396695e+02  1.37313928e+03    6.77590397e+01   -1.84199859e+03 -3.92514325e+01    6.77590396e+01   -5.40152606e+02 -3.77930375e+02    6.79158444e+01   -9.94911201e+00  1.00654153e+01    6.79158418e+01    1.10396695e+02 -1.37313928e+03    6.77590396e+01    1.56033329e+03  8.54822647e+02    6.79158384e+01   -1.10196673e+04  2.64309485e+04    6.77590346e+01   -5.29343133e+04 -3.08261539e+03    6.69787554e+01    0.00000000e+00  0.00000000e+00    6.69787554e+01    0.00000000e+00  0.00000000e+00    6.69787554e+01    0.00000000e+00  0.00000000e+00    6.69787554e+01    0.00000000e+00  0.00000000e+00    6.79158384e+01    5.50983365e+03 -1.32154742e+04    6.77587415e+01   -4.24087275e+02 -2.42524261e+02    6.77587390e+01    2.56342071e+04  1.28036455e+03    6.79158384e+01   -1.15275475e+04  2.76700355e+04
82301716e+01    4.12314067e+02 -6.56916970e+03    6.80730202e+01    1.42957704e+04  7.52511553e+02    6.82301786e+01   -9.17488452e+01  1.37280737e+03    6.80730266e+01   -1.91292566e+03 -3.80489026e+01    6.80730263e+01   -6.09046112e+02 -3.76795984e+02    6.82301840e+01   -1.00012819e+01  1.00663008e+01    6.82301786e+01    9.17488451e+01 -1.37280738e+03    6.80730264e+01    1.45586461e+03  8.56561050e+02    6.82301716e+01   -1.64925627e+03  2.62766788e+04    6.80730159e+01   -5.20574140e+04 -3.09206214e+03    6.70625312e+01    0.00000000e+00  0.00000000e+00    6.70625312e+01    0.00000000e+00  0.00000000e+00    6.70625312e+01    0.00000000e+00  0.00000000e+00    6.70625312e+01    0.00000000e+00  0.00000000e+00    6.82301716e+01    8.24628135e+02 -1.31383394e+04    6.80727365e+01   -4.16972207e+02 -2.42645064e+02    6.80727311e+01    2.69032041e+04  1.25776701e+03    6.82301717e+01   -1.86317393e+04  2.77882298e+04
85441574e+01    2.75969241e+03 -6.60780473e+03    6.83873551e+01    1.19490784e+04  7.91197272e+02    6.85441608e+01   -1.10334713e+02  1.37313784e+03    6.83873581e+01   -1.84220064e+03 -3.92473239e+01    6.83873579e+01   -5.40329627e+02 -3.77927135e+02    6.85441633e+01   -9.94949974e+00  1.00654224e+01    6.85441608e+01    1.10334713e+02 -1.37313785e+03    6.83873580e+01    1.56005442e+03  8.54827947e+02    6.85441574e+01   -1.10387696e+04  2.64312189e+04    6.83873531e+01   -5.29389712e+04 -3.08250130e+03    6.71463070e+01    0.00000000e+00  0.00000000e+00    6.71463070e+01    0.00000000e+00  0.00000000e+00    6.71463070e+01    0.00000000e+00  0.00000000e+00    6.71463070e+01    0.00000000e+00  0.00000000e+00    6.85441574e+01    5.51938481e+03 -1.32156095e+04    6.83870626e+01   -4.24066907e+02 -2.42524673e+02    6.83870601e+01    2.56325936e+04  1.28037059e+03    6.85441575e+01   -1.15137729e+04  2.76698511e+04
88584906e+01    4.07542471e+02 -6.56910204e+03    6.87013386e+01    1.43005163e+04  7.52445358e+02    6.88584975e+01   -9.18101735e+01  1.37280879e+03    6.87013449e+01   -1.91272540e+03 -3.80529665e+01    6.87013446e+01   -6.08870468e+02 -3.76799195e+02    6.88585029e+01   -1.00008989e+01  1.00662939e+01    6.88584975e+01    9.18101734e+01 -1.37280879e+03    6.87013447e+01    1.45614120e+03  8.56555801e+02    6.88584906e+01   -1.63016988e+03  2.62764082e+04    6.87013345e+01   -5.20527874e+04 -3.09217515e+03    6.72300828e+01    0.00000000e+00  0.00000000e+00    6.72300828e+01    0.00000000e+00  0.00000000e+00    6.72300828e+01    0.00000000e+00  0.00000000e+00    6.72300828e+01    0.00000000e+00  0.00000000e+00    6.88584906e+01    8.15084942e+02 -1.31382041e+04    6.87010575e+01   -4.16992396e+02 -2.42644656e+02    6.87010521e+01    2.69048258e+04  1.25776055e+03    6.88584907e+01   -1.86455085e+04  2.77884146e+04
91724765e+01    2.76445614e+03 -6.60787239e+03    6.90156734e+01    1.19443397e+04  7.91263502e+02    6.91724798e+01   -1.10274069e+02  1.37313645e+03    6.90156764e+01   -1.84239895e+03 -3.92433074e+01    6.90156763e+01   -5.40503767e+02 -3.77923956e+02    6.91724823e+01   -9.94987775e+00  1.00654292e+01    6.91724798e+01    1.10274069e+02 -1.37313645e+03    6.90156763e+01    1.55978032e+03  8.54833142e+02    6.91724765e+01   -1.10578246e+04  2.64314896e+04    6.90156715e+01   -5.29435632e+04 -3.08238945e+03    6.73138586e+01    0.00000000e+00  0.00000000e+00    6.73138586e+01    0.00000000e+00  0.00000000e+00    6.73138586e+01    0.00000000e+00  0.00000000e+00    6.73138586e+01    0.00000000e+00  0.00000000e+00    6.91724765e+01    5.52891228e+03 -1.32157448e+04    6.90153836e+01   -4.24046914e+02 -2.42525076e+02    6.90153812e+01    2.56309651e+04  1.28037746e+03    6.91724765e+01   -1.15000206e+04  2.76696661e+04
94868096e+01    4.02782821e+02 -6.56903432e+03    6.93296570e+01    1.43052514e+04  7.52379047e+02    6.94868165e+01   -9.18701877e+01  1.37281017e+03    6.93296633e+01   -1.91252883e+03 -3.80569401e+01    6.93296629e+01   -6.08697659e+02 -3.76802345e+02    6.94868218e+01   -1.00005255e+01  1.00662871e+01    6.94868165e+01    9.18701877e+01 -1.37281017e+03    6.93296630e+01    1.45641310e+03  8.56550655e+02    6.94868096e+01   -1.61113128e+03  2.62761373e+04    6.93296530e+01   -5.20482256e+04 -3.09228597e+03    6.73976344e+01    0.00000000e+00  0.00000000e+00    6.73976344e+01    0.00000000e+00  0.00000000e+00    6.73976344e+01    0.00000000e+00  0.00000000e+00    6.73976344e+01    0.00000000e+00  0.00000000e+00    6.94868096e+01    8.05565642e+02 -1.31380686e+04    6.93293784e+01   -4.17012216e+02 -2.42644257e+02    6.93293730e+01    2.69064618e+04  1.25775328e+03    6.94868098e+01   -1.86592549e+04  2.77886000e+04
...

记录了不同震源、不同积分类型的峰谷位置。

故为了绘制完整的波数积分+PTAM过程，涉及三个文件。不过Python函数已经做了优化，由于程序输出的目录结构固定，文件命名格式固定，只需传入 PTAM_{ir}_{dist}/PTAM_{iw}_{freq} 文件，Python函数会自动读取对应的三个文件。

ir = 2
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(
    f"GRN_grtstats/milrow_{format_float(depsrc)}_{format_float(deprcv)}/PTAM_{ir:04d}_*/PTAM_0050_*"
)

srctype="SS"
ptype="0"
fig, ax = pygrt.utils.plot_statsdata_ptam(statsdata1, statsdata2, ptamdata, dist=dist, srctype=srctype, ptype=ptype, RorI=2)
fig.savefig(f"{srctype}_{ptype}_{depsrc}_ptam_RI.svg", bbox_inches='tight')

红十字为选取的波峰波谷，再取缩减序列即可得到估计的积分收敛值。

静态解的积分收敛性

以上部分是以动态解为例，静态解从积分类型、收敛特征、文件格式、绘图完全类似，只是不再有频率索引值。

备注

核函数文件中记录的值非理论核函数值。对于静态解，还需乘 \(\left(\dfrac{1}{4\pi\mu}\right)\)。

假设震源深度0.05km，场点位于地表，场点仅定义一个点(2,2)做示例，这里直接给出脚本。

# -S 表示输出核函数文件
grt static greenfn -Mmilrow -D0.05/0 -X2/2/1 -Y2/2/1 -Cp -K+k3+f -S -Ostgrn.nc

# grt.ker2asc 也可以读取静态解输出的核函数文件，格式一致
grt ker2asc stgrtstats/milrow_0.05_0/K > static_stats

# 绘制图像部分见Python

import numpy as np
import pygrt
from pygrt.cli import format_float

depsrc = 0.05
deprcv = 0.0
pymod = pygrt.PyModel1D(stgrn="stgrn.nc", modelpath="milrow")

norths = [2.0, 2.0, 1.0]
easts = [2.0, 2.0, 1.0]
pymod.static_greenfn(depsrc=depsrc, deprcv=deprcv, norths=norths, easts=easts, converg_method='PTAM', stats=True, k0=3, use_kmax_ref=True)

ir = 0
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(
    f"stgrtstats/milrow_{format_float(depsrc)}_{format_float(deprcv)}/PTAM_{ir:04d}_*/PTAM"
)

srctype="SS"
ptype="0"
# 只使用离散波数积分的积分变化
fig, ax = pygrt.utils.plot_statsdata(statsdata1, dist=dist, srctype=srctype, ptype=ptype, RorI=True)
fig.tight_layout()
fig.savefig(f"{srctype}_{ptype}_{depsrc}_static.svg", bbox_inches='tight')

# 使用了峰谷平均法的积分变化
fig, ax = pygrt.utils.plot_statsdata_ptam(statsdata1, statsdata2, ptamdata, dist=dist, srctype=srctype, ptype=ptype, RorI=True)
fig.tight_layout()
fig.savefig(f"{srctype}_{ptype}_{depsrc}_ptam_static.svg", bbox_inches='tight')

只使用离散波数积分

使用峰谷平均法