Skip to content

FEAT: unify PTAM wavenumber integration with DWM #284

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Dengda98 merged 1 commit into from
Aug 7, 2026
Merged

FEAT: unify PTAM wavenumber integration with DWM #284

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/source/Advanced/integ_converg/ptam.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@

具体原理很简洁易懂，从以下图像中你也能了解大概。详见 (|zhang2003|; |zhang2021|) 。

具体流程为，程序中在进行完离散波数积分后，继续增加k（不同震中距使用不同dk），
具体流程为，程序中在进行完离散波数积分后，继续增加k（同一频率下所有震中距复用DWM的统一dk），
使用PTAM寻找足够数量的波峰和波谷（内部设定为36个），
再对这些波峰波谷取缩减序列 :math:`M_i \leftarrow 0.5\times(M_i + M_{i+1})` ，得到估计的积分收敛值。
取缩减序列的C代码如下，
Expand Down
1 change: 0 additions & 1 deletion pygrt/C_extension/include/grt/common/const.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,6 @@ typedef double complex cplx_t;

#define GRT_PTAM_PT_MAX 36 ///< 36， 最后统计波峰波谷的目标数量
#define GRT_PTAM_WINDOW_SIZE 3 ///< 3， 使用连续点数判断是否为波峰或波谷
#define GRT_PTAM_WAITS_MAX 9 ///< 9, 判断波峰或波谷的最大等待次数，不能太小

#define GRT_INVERSE_SUCCESS 0 ///< 求逆或除法没有遇到除0错误
#define GRT_INVERSE_FAILURE -1 ///< 求逆或除法遇到除0错误
Expand Down
5 changes: 2 additions & 3 deletions pygrt/C_extension/include/grt/integral/ptam.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,7 @@
*
* @param[in,out] mstat `MODEL1D_STATE` 结构体指针
* @param[in] k0 先前的积分已经进行到了波数k0
* @param[in] predk 先前的积分使用的积分间隔dk，因为峰谷平均法使用的
* 积分间隔会和之前的不一致，这里传入该系数以做预先调整
* @param[in] dk 先前积分使用的统一波数间隔，PTAM继续复用该间隔
* @param[in] nr 震中距数量
* @param[in] rs 震中距数组
*
Expand All @@ -41,6 +40,6 @@
*
*/
void grt_PTA_method(
MODEL1D_STATE *mstat, real_t k0, real_t predk,
MODEL1D_STATE *mstat, real_t k0, real_t dk,
size_t nr, real_t *rs, K_INTEG *K, FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc);

2 changes: 1 addition & 1 deletion pygrt/C_extension/src/integral/integ_process.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ void grt_KPROC_init_fstats(
GRT_SAFE_ASPRINTF(&fname, "%s/K%s", statsstr, suffix);
Kproc->fstats = fopen(fname, "wb");

// PTAM的积分中间结果, 每个震中距两个文件，因为PTAM对不同震中距使用不同的dk
// PTAM的积分中间结果，每个震中距两个文件，保持现有统计文件布局
// 在文件名后加后缀，区分不同震中距
char *ptam_dirname = NULL;
if(Kproc->cvgmet == K_INTEG_CONVERG_PTAM){
Expand Down
162 changes: 98 additions & 64 deletions pygrt/C_extension/src/integral/ptam.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,40 +104,38 @@ static int _cplx_peak_or_trough(
* @param[in] im 不同震源不同阶数的索引
* @param[in] v 积分形式索引
* @param[in] k 波数
* @param[in] dk 波数步长
* @param[in] dk 波数步长
* @param[in] waits_max 当前震中距允许的最大等待采样数
* @param[in] J3 存储的积分采样幅值数组
* @param[in,out] Kpt 积分值峰谷的波数数组
* @param[in,out] Fpt 用于存储波峰/波谷点的幅值数组
* @param[in,out] Ipt 用于存储波峰/波谷点的个数数组
* @param[in,out] Gpt 用于存储等待迭次数的数组
* @param[in,out] iendk0 一个布尔指针，用于指示是否满足结束条件
*/
static void process_peak_or_trough(
size_t ir, int im, int v, real_t k, real_t dk,
size_t ir, int im, int v, real_t k, real_t dk, size_t waits_max,
cplxIntegGrid (*J3)[GRT_PTAM_WINDOW_SIZE], realIntegGrid (*Kpt)[GRT_PTAM_PT_MAX],
cplxIntegGrid (*Fpt)[GRT_PTAM_PT_MAX], sizeIntegGrid (*Ipt), sizeIntegGrid (*Gpt), bool *iendk0)
cplxIntegGrid (*Fpt)[GRT_PTAM_PT_MAX], sizeIntegGrid (*Ipt), sizeIntegGrid (*Gpt))
{
cplx_t tmp0;
if (Gpt[ir][im][v] >= GRT_PTAM_WINDOW_SIZE-1 && Ipt[ir][im][v] < GRT_PTAM_PT_MAX) {
if (_cplx_peak_or_trough(im, v, J3[ir], k, dk, &Kpt[ir][Ipt[ir][im][v]][im][v], &tmp0) != 0) {
Fpt[ir][Ipt[ir][im][v]++][im][v] = tmp0;
Gpt[ir][im][v] = 0;
} else if (Gpt[ir][im][v] >= GRT_PTAM_WAITS_MAX) { // 不再等待，直接取中点作为波峰波谷
} else if (Gpt[ir][im][v] >= waits_max) { // 不再等待，直接取中点作为波峰波谷
Kpt[ir][Ipt[ir][im][v]][im][v] = k - dk;
Fpt[ir][Ipt[ir][im][v]++][im][v] = J3[ir][1][im][v];
Gpt[ir][im][v] = 0;
}
}
*iendk0 = *iendk0 && (Ipt[ir][im][v] == GRT_PTAM_PT_MAX);
}


/**
* 在输入被积函数的情况下，对不同震源使用峰谷平均法
*
* @param[in] ir 震中距索引
* @param[in] nr 震中距个数
* @param[in] precoef 积分值系数
* @param[in] waits_max 当前震中距允许的最大等待采样数
* @param[in] k 波数
* @param[in] dk 波数步长
* @param[in,out] SUM3 被积函数的幅值数组
Expand All @@ -148,32 +146,27 @@ static void process_peak_or_trough(
* @param[in,out] Ipt 用于存储波峰/波谷点的个数数组
* @param[in,out] Gpt 用于存储等待迭次数的数组
*
* @param[in,out] iendk0 是否收集足够峰谷
*
*/
static void ptam_once(
const size_t ir, const size_t nr, const real_t precoef, real_t k, real_t dk,
cplxIntegGrid SUM3[nr][GRT_PTAM_WINDOW_SIZE],
cplxIntegGrid sumJ[nr],
realIntegGrid Kpt[nr][GRT_PTAM_PT_MAX],
cplxIntegGrid Fpt[nr][GRT_PTAM_PT_MAX],
sizeIntegGrid Ipt[nr],
sizeIntegGrid Gpt[nr],
bool *iendk0)
static bool ptam_once(
const size_t ir, const size_t waits_max, real_t k, real_t dk,
cplxIntegGrid (*SUM3)[GRT_PTAM_WINDOW_SIZE],
cplxIntegGrid *sumJ,
realIntegGrid (*Kpt)[GRT_PTAM_PT_MAX],
cplxIntegGrid (*Fpt)[GRT_PTAM_PT_MAX],
sizeIntegGrid *Ipt,
sizeIntegGrid *Gpt)
{
*iendk0 = true;

GRT_LOOP_IntegGrid(im, v){
int modr = GRT_SRC_M_ORDERS[im];
if(modr == 0 && v!=0 && v!= 2) continue;

// 赋更新量
// SUM3转为求和结果
sumJ[ir][im][v] += SUM3[ir][GRT_PTAM_WINDOW_SIZE-1][im][v] * precoef;
sumJ[ir][im][v] += SUM3[ir][GRT_PTAM_WINDOW_SIZE-1][im][v];
SUM3[ir][GRT_PTAM_WINDOW_SIZE-1][im][v] = sumJ[ir][im][v];

// 3点以上，判断波峰波谷
process_peak_or_trough(ir, im, v, k, dk, SUM3, Kpt, Fpt, Ipt, Gpt, iendk0);
process_peak_or_trough(ir, im, v, k, dk, waits_max, SUM3, Kpt, Fpt, Ipt, Gpt);

// 左移动点,
for(int jj=0; jj<GRT_PTAM_WINDOW_SIZE-1; ++jj){
Expand All @@ -183,6 +176,26 @@ static void ptam_once(
// 点数+1
Gpt[ir][im][v]++;
}

bool iendk0 = true;
GRT_LOOP_IntegGrid(im, v){
int modr = GRT_SRC_M_ORDERS[im];
if(modr == 0 && v!=0 && v!= 2) continue;
iendk0 = iendk0 && (Ipt[ir][im][v] == GRT_PTAM_PT_MAX);
}
return iendk0;
}


static bool ptam_is_complete(size_t ir, const sizeIntegGrid *Ipt)
{
bool iendk0 = true;
GRT_LOOP_IntegGrid(im, v){
int modr = GRT_SRC_M_ORDERS[im];
if(modr == 0 && v!=0 && v!= 2) continue;
iendk0 = iendk0 && (Ipt[ir][im][v] == GRT_PTAM_PT_MAX);
}
return iendk0;
}


Expand Down Expand Up @@ -212,14 +225,24 @@ static void _cplx_shrink(size_t n1, size_t ir, int im, int v, cplxIntegGrid (*F


void grt_PTA_method(
MODEL1D_STATE *mstat, real_t k0, real_t predk,
MODEL1D_STATE *mstat, real_t k0, real_t dk,
size_t nr, real_t *rs, K_INTEG *K, FILE *ptam_fstatsnr[nr][2], GRT_KernelFunc kerfunc)
{
// 需要兼容对正常收敛而不具有规律波峰波谷的序列
// 有时序列收敛比较好，不表现为规律的波峰波谷，
// 此时设置最大等待次数，超过直接设置为中间值

real_t k=0.0;
real_t rmin=0.0;
size_t nvalid=0;

for(size_t ir = 0; ir < nr; ++ir){
if(GRT_IS_ZERO(rs[ir])) continue;
if(nvalid == 0 || rs[ir] < rmin) rmin = rs[ir];
nvalid++;
}

if(nvalid == 0 || dk <= 0.0) return;

// 使用宏函数，方便定义
#define __CALLOC_ARRAY(VAR, TYP, __ARR) \
Expand Down Expand Up @@ -259,56 +282,66 @@ void grt_PTA_method(
#undef __ARR
#undef __CALLOC_ARRAY

// 对于PTAM，不同震中距使用不同dk
size_t *waits_max = (size_t *)calloc(nr, sizeof(*waits_max));
bool *iendkrs = (bool *)calloc(nr, sizeof(*iendkrs));
for(size_t ir = 0; ir < nr; ++ir){

// 零震中距：PTAM 步长含 1/r；震源-接收点重合时格林函数奇异
if(GRT_IS_ZERO(rs[ir])) continue;
waits_max[ir] = (size_t)ceil(PI/(rs[ir]*dk)) + 5; // + 5 以防万一
if(waits_max[ir] < GRT_PTAM_WINDOW_SIZE-1){
waits_max[ir] = GRT_PTAM_WINDOW_SIZE-1;
}
}

// 由最小非零震中距保证所有距离都有足够的峰谷
real_t kmax = k0 + (GRT_PTAM_PT_MAX+5)*PI/rmin;
bool iendk = false;
k = k0;
while(!iendk && k <= kmax){
k += dk;

real_t dk = PI/((GRT_PTAM_WAITS_MAX-1)*rs[ir]);
real_t precoef = dk/predk; // 提前乘dk系数，以抵消格林函数主函数计算时最后乘dk
// 根据波峰波谷的目标也给出一个kmax，+5以防万一
real_t kmax = k0 + (GRT_PTAM_PT_MAX+5)*PI/rs[ir];

bool iendk0=false;
// 计算一次核函数，供所有尚未完成的震中距复用
kerfunc(mstat, k, K->QWV, K->calc_upar, K->QWVz);
if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN;

k = k0;
while(true){
if(k > kmax) break;
k += dk;
iendk = true;
for(size_t ir = 0; ir < nr; ++ir){
if(GRT_IS_ZERO(rs[ir]) || iendkrs[ir]) continue;

// 计算核函数 F(k, w)
kerfunc(mstat, k, K->QWV, K->calc_upar, K->QWVz);
if(mstat->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN;
// 记录仍在处理的震中距对应的核函数
if(ptam_fstatsnr != NULL){
grt_write_stats(ptam_fstatsnr[ir][0], k, (K->calc_upar)? K->QWVz : K->QWV);
}

// 记录核函数
if(ptam_fstatsnr != NULL) grt_write_stats(ptam_fstatsnr[ir][0], k, (K->calc_upar)? K->QWVz : K->QWV);
bool iendk0 = ptam_is_complete(ir, Ipt);
if(!iendk0){
grt_int_Pk(k, rs[ir], K->QWV, false, SUM3[ir][GRT_PTAM_WINDOW_SIZE-1]);
iendk0 = ptam_once(
ir, waits_max[ir], k, dk, SUM3, K->sumJ, Kpt, Fpt, Ipt, Gpt);
}

// 计算被积函数一项 F(k,w)Jm(kr)k
grt_int_Pk(k, rs[ir], K->QWV, false, SUM3[ir][GRT_PTAM_WINDOW_SIZE-1]); // [GRT_PTAM_WINDOW_SIZE-1]表示把新点值放在最后
// 判断和记录波峰波谷
ptam_once(ir, nr, precoef, k, dk, SUM3, K->sumJ, Kpt, Fpt, Ipt, Gpt, &iendk0);

// -------------------------- 位移空间导数 ------------------------------------
if(K->calc_upar){
// ------------------------------- ui_z -----------------------------------
// 计算被积函数一项 F(k,w)Jm(kr)k
grt_int_Pk(k, rs[ir], K->QWVz, false, SUM3_uiz[ir][GRT_PTAM_WINDOW_SIZE-1]); // [GRT_PTAM_WINDOW_SIZE-1]表示把新点值放在最后
// 判断和记录波峰波谷
ptam_once(ir, nr, precoef, k, dk, SUM3_uiz, K->sumJz, Kpt_uiz, Fpt_uiz, Ipt_uiz, Gpt_uiz, &iendk0);

// ------------------------------- ui_r -----------------------------------
// 计算被积函数一项 F(k,w)Jm(kr)k
grt_int_Pk(k, rs[ir], K->QWV, true, SUM3_uir[ir][GRT_PTAM_WINDOW_SIZE-1]); // [GRT_PTAM_WINDOW_SIZE-1]表示把新点值放在最后
// 判断和记录波峰波谷
ptam_once(ir, nr, precoef, k, dk, SUM3_uir, K->sumJr, Kpt_uir, Fpt_uir, Ipt_uir, Gpt_uir, &iendk0);

} // END if calc_upar

bool iendk_uiz = ptam_is_complete(ir, Ipt_uiz);
if(!iendk_uiz){
grt_int_Pk(k, rs[ir], K->QWVz, false, SUM3_uiz[ir][GRT_PTAM_WINDOW_SIZE-1]);
iendk_uiz = ptam_once(
ir, waits_max[ir], k, dk, SUM3_uiz, K->sumJz,
Kpt_uiz, Fpt_uiz, Ipt_uiz, Gpt_uiz);
}

bool iendk_uir = ptam_is_complete(ir, Ipt_uir);
if(!iendk_uir){
grt_int_Pk(k, rs[ir], K->QWV, true, SUM3_uir[ir][GRT_PTAM_WINDOW_SIZE-1]);
iendk_uir = ptam_once(
ir, waits_max[ir], k, dk, SUM3_uir, K->sumJr,
Kpt_uir, Fpt_uir, Ipt_uir, Gpt_uir);
}

iendk0 = iendk0 && iendk_uiz && iendk_uir;
}

if(iendk0) break;
}// end k loop
iendkrs[ir] = iendk0;
iendk = iendk && iendkrs[ir];
}
}

// 做缩减序列，赋值最终解
Expand Down Expand Up @@ -339,6 +372,7 @@ void grt_PTA_method(
X(Kpt) X(Fpt) X(Ipt) X(Gpt) \
X(Kpt_uiz) X(Fpt_uiz) X(Ipt_uiz) X(Gpt_uiz) \
X(Kpt_uir) X(Fpt_uir) X(Ipt_uir) X(Gpt_uir) \
X(waits_max) X(iendkrs) \

#define X(A) GRT_SAFE_FREE_PTR(A);
__FREE_ALL_ARRAY
Expand Down
Loading