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

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

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

具体流程为，程序中在进行完离散波数积分后，继续增加k（不同震中距使用不同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 -S50,100
grt ker2asc GRN_grtstats/milrow_0_0/K_0050_5.00000e+00 > stats
# 绘制图像部分见Python

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

import numpy as np
import pygrt 

modarr = np.loadtxt("milrow")

depsrc = 0.0
deprcv = 0.0
pymod = pygrt.PyModel1D(modarr, depsrc=depsrc, deprcv=deprcv)

distarr = [5,8,10]
# 设置 converg_method='PTAM' 进行收敛
stgrnLst = pymod.compute_grn(
    distarr=distarr, nt=500, dt=0.02, converg_method='PTAM',
    statsfile=f"pygrtstats_{depsrc}_{deprcv}", statsidxs=[50,100]
)

输出的核函数文件会在自定义路径下。

在 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
67622465e+01    6.56540447e+02 -6.56685940e+03    2.66053579e+01    1.41196277e+04  7.46104302e+02    2.67622568e+01   -6.79904618e+01  1.37174989e+03    2.66053633e+01   -1.96993245e+03 -3.58939919e+01    2.66053604e+01   -6.45663838e+02 -3.75692774e+02    2.67622542e+01   -1.01275137e+01  1.00693959e+01    2.67622568e+01    6.79904618e+01 -1.37174989e+03    2.66053605e+01    1.39523921e+03  8.58335768e+02    2.67622465e+01   -2.62616179e+03  2.62674376e+04    2.66053768e+01   -5.32380923e+04 -3.03626722e+03    2.65098060e+01    0.00000000e+00  0.00000000e+00    2.65098060e+01    0.00000000e+00  0.00000000e+00    2.65098060e+01    0.00000000e+00  0.00000000e+00    2.65098060e+01    0.00000000e+00  0.00000000e+00    2.67622465e+01    1.31308089e+03 -1.31337188e+04    2.66045993e+01   -4.12169613e+02 -2.42779393e+02    2.66046025e+01    2.71794665e+04  1.23488996e+03    2.67622465e+01   -1.80244975e+04  2.77976722e+04
70764072e+01    2.51518050e+03 -6.60974264e+03    2.69195163e+01    1.21600718e+04  7.96128284e+02    2.70764171e+01   -1.33187553e+02  1.37413124e+03    2.69195216e+01   -1.78972293e+03 -4.12001618e+01    2.69195189e+01   -5.06555790e+02 -3.78932149e+02    2.70764150e+01   -9.82636763e+00  1.00624197e+01    2.70764171e+01    1.33187553e+02 -1.37413124e+03    2.69195190e+01    1.61614488e+03  8.53196237e+02    2.70764072e+01   -1.00607220e+04  2.64389706e+04    2.69195353e+01   -5.18251329e+04 -3.13355235e+03    2.65883458e+01    0.00000000e+00  0.00000000e+00    2.65883458e+01    0.00000000e+00  0.00000000e+00    2.65883458e+01    0.00000000e+00  0.00000000e+00    2.65883458e+01    0.00000000e+00  0.00000000e+00    2.70764072e+01    5.03036100e+03 -1.32194853e+04    2.69187667e+01   -4.28530947e+02 -2.42399315e+02    2.69187695e+01    2.54018812e+04  1.30068019e+03    2.70764073e+01   -1.21219357e+04  2.76612826e+04
73905679e+01    6.63922165e+02 -6.56776606e+03    2.72336747e+01    1.41052157e+04  7.47593107e+02    2.73905774e+01   -7.01923475e+01  1.37188772e+03    2.72336799e+01   -1.96554741e+03 -3.61457771e+01    2.72336774e+01   -6.43462639e+02 -3.75801285e+02    2.73905757e+01   -1.01192231e+01  1.00691645e+01    2.73905774e+01    7.01923475e+01 -1.37188772e+03    2.72336775e+01    1.39878126e+03  8.58185575e+02    2.73905679e+01   -2.65568866e+03  2.62710642e+04    2.72336938e+01   -5.31428681e+04 -3.04282370e+03    2.66668856e+01    0.00000000e+00  0.00000000e+00    2.66668856e+01    0.00000000e+00  0.00000000e+00    2.66668856e+01    0.00000000e+00  0.00000000e+00    2.66668856e+01    0.00000000e+00  0.00000000e+00    2.73905679e+01    1.32784433e+03 -1.31355321e+04    2.72329339e+01   -4.12437856e+02 -2.42771055e+02    2.72329362e+01    2.71381020e+04  1.23806921e+03    2.73905680e+01   -1.79994311e+04  2.77953276e+04
77047285e+01    2.50887510e+03 -6.60893377e+03    2.75478332e+01    1.21727432e+04  7.94803024e+02    2.77047378e+01   -1.31157183e+02  1.37400760e+03    2.75478383e+01   -1.79379310e+03 -4.09736351e+01    2.75478359e+01   -5.08623107e+02 -3.78833808e+02    2.77047363e+01   -9.83427476e+00  1.00626416e+01    2.77047378e+01    1.31157183e+02 -1.37400760e+03    2.75478361e+01    1.61279151e+03  8.53335164e+02    2.77047285e+01   -1.00355004e+04  2.64357351e+04    2.75478524e+01   -5.19131778e+04 -3.12765352e+03    2.67454255e+01    0.00000000e+00  0.00000000e+00    2.67454255e+01    0.00000000e+00  0.00000000e+00    2.67454255e+01    0.00000000e+00  0.00000000e+00    2.67454255e+01    0.00000000e+00  0.00000000e+00    2.77047285e+01    5.01775019e+03 -1.32178675e+04    2.75471010e+01   -4.28273667e+02 -2.42407434e+02    2.75471029e+01    2.54396905e+04  1.29782936e+03    2.77047287e+01   -1.21441967e+04  2.76634249e+04
80188892e+01    6.69270525e+02 -6.56849017e+03    2.78619918e+01    1.40940840e+04  7.48777392e+02    2.80188982e+01   -7.20703623e+01  1.37199912e+03    2.78619967e+01   -1.96175764e+03 -3.63505303e+01    2.78619945e+01   -6.41515349e+02 -3.75890835e+02    2.80188970e+01   -1.01116777e+01  1.00689522e+01    2.80188982e+01    7.20703623e+01 -1.37199912e+03    2.78619946e+01    1.40196221e+03  8.58056686e+02    2.80188892e+01   -2.67708210e+03  2.62739607e+04    2.78620111e+01   -5.30611952e+04 -3.04815521e+03    2.68239653e+01    0.00000000e+00  0.00000000e+00    2.68239653e+01    0.00000000e+00  0.00000000e+00    2.68239653e+01    0.00000000e+00  0.00000000e+00    2.68239653e+01    0.00000000e+00  0.00000000e+00    2.80188892e+01    1.33854105e+03 -1.31369803e+04    2.78612679e+01   -4.12684655e+02 -2.42763185e+02    2.78612694e+01    2.71034445e+04  1.24063681e+03    2.80188894e+01   -1.79797190e+04  2.77933658e+04
83330498e+01    2.50438129e+03 -6.60828353e+03    2.81761503e+01    1.21825071e+04  7.93740980e+02    2.83330585e+01   -1.29415098e+02  1.37390681e+03    2.81761552e+01   -1.79733203e+03 -4.07877656e+01    2.81761531e+01   -5.10462284e+02 -3.78751904e+02    2.83330576e+01   -9.84147934e+00  1.00628443e+01    2.83330585e+01    1.29415098e+02 -1.37390681e+03    2.81761532e+01    1.60976858e+03  8.53455072e+02    2.83330498e+01   -1.00175251e+04  2.64331341e+04    2.81761698e+01   -5.19891692e+04 -3.12281433e+03    2.69025051e+01    0.00000000e+00  0.00000000e+00    2.69025051e+01    0.00000000e+00  0.00000000e+00    2.69025051e+01    0.00000000e+00  0.00000000e+00    2.69025051e+01    0.00000000e+00  0.00000000e+00    2.83330498e+01    5.00876257e+03 -1.32165671e+04    2.81754346e+01   -4.28036856e+02 -2.42415036e+02    2.81754358e+01    2.54715411e+04  1.29550748e+03    2.83330500e+01   -1.21615863e+04  2.76652250e+04
86472104e+01    6.72998358e+02 -6.56907568e+03    2.84903089e+01    1.40855403e+04  7.49732857e+02    2.86472189e+01   -7.36906788e+01  1.37209065e+03    2.84903137e+01   -1.95844405e+03 -3.65199298e+01    2.84903117e+01   -6.39773974e+02 -3.75966048e+02    2.86472181e+01   -1.01047941e+01  1.00687588e+01    2.86472189e+01    7.36906788e+01 -1.37209066e+03    2.84903119e+01    1.40483988e+03  8.57944843e+02    2.86472104e+01   -2.69199343e+03  2.62763027e+04    2.84903286e+01   -5.29902910e+04 -3.05256478e+03    2.69810449e+01    0.00000000e+00  0.00000000e+00    2.69810449e+01    0.00000000e+00  0.00000000e+00    2.69810449e+01    0.00000000e+00  0.00000000e+00    2.69810449e+01    0.00000000e+00  0.00000000e+00    2.86472104e+01    1.34599672e+03 -1.31381514e+04    2.84896011e+01   -4.12911972e+02 -2.42755857e+02    2.84896021e+01    2.70741042e+04  1.24274434e+03    2.86472106e+01   -1.79644510e+04  2.77917112e+04
89613710e+01    2.50134247e+03 -6.60775503e+03    2.88044676e+01    1.21899582e+04  7.92878912e+02    2.89613792e+01   -1.27904277e+02  1.37382339e+03    2.88044722e+01   -1.80044237e+03 -4.06328061e+01    2.88044703e+01   -5.12114829e+02 -3.78682575e+02    2.89613787e+01   -9.84806065e+00  1.00630286e+01    2.89613792e+01    1.27904277e+02 -1.37382339e+03    2.88044705e+01    1.60702479e+03  8.53559648e+02    2.89613710e+01   -1.00053699e+04  2.64310201e+04    2.88044875e+01   -5.20555001e+04 -3.11878160e+03    2.70595847e+01    0.00000000e+00  0.00000000e+00    2.70595847e+01    0.00000000e+00  0.00000000e+00    2.70595847e+01    0.00000000e+00  0.00000000e+00    2.70595847e+01    0.00000000e+00  0.00000000e+00    2.89613710e+01    5.00268495e+03 -1.32155101e+04    2.88037675e+01   -4.27818551e+02 -2.42422087e+02    2.88037682e+01    2.54986276e+04  1.29358791e+03    2.89613713e+01   -1.21749115e+04  2.76667484e+04
92755316e+01    6.75415340e+02 -6.56955375e+03    2.91186262e+01    1.40790709e+04  7.50512696e+02    2.92755395e+01   -7.51026933e+01  1.37216693e+03    2.91186308e+01   -1.95551771e+03 -3.66621676e+01    2.91186290e+01   -6.38202411e+02 -3.76030176e+02    2.92755392e+01   -1.00984974e+01  1.00685831e+01    2.92755395e+01    7.51026933e+01 -1.37216693e+03    2.91186292e+01    1.40745998e+03  8.57846838e+02    2.92755316e+01   -2.70166136e+03  2.62782150e+04    2.91186465e+01   -5.29280866e+04 -3.05626541e+03    2.71381245e+01    0.00000000e+00  0.00000000e+00    2.71381245e+01    0.00000000e+00  0.00000000e+00    2.71381245e+01    0.00000000e+00  0.00000000e+00    2.71381245e+01    0.00000000e+00  0.00000000e+00    2.92755316e+01    1.35083068e+03 -1.31391075e+04    2.91179337e+01   -4.13121732e+02 -2.42749081e+02    2.91179342e+01    2.70490484e+04  1.24449832e+03    2.92755318e+01   -1.79529089e+04  2.77903065e+04
...

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

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

ir = 2
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(f"pygrtstats_{depsrc}_{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.1km，场点位于地表，场点仅定义一个点(2,2)做示例，这里直接给出脚本。

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

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

# 绘制图像部分见Python

import numpy as np
import pygrt 

modarr = np.loadtxt("milrow")

depsrc = 0.1
deprcv = 0.0
pymod = pygrt.PyModel1D(modarr, depsrc=depsrc, deprcv=deprcv)

xarr = np.array([2.0])
yarr = np.array([2.0])
static_grn = pymod.compute_static_grn(xarr, yarr, converg_method='PTAM', statsfile=f"static_pygrtstats_{depsrc}_{deprcv}")

ir = 0
statsdata1, statsdata2, ptamdata, dist = pygrt.utils.read_statsfile_ptam(f"static_pygrtstats_{depsrc}_{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')

只使用离散波数积分

使用峰谷平均法