From 8b69681d4a8e29606e3efbd3c93d36473d24e51e Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Fri, 7 Aug 2026 15:28:25 +0800 Subject: [PATCH] FEAT: unify PTAM wavenumber integration with DWM --- docs/source/Advanced/integ_converg/ptam.rst | 2 +- pygrt/C_extension/include/grt/common/const.h | 1 - pygrt/C_extension/include/grt/integral/ptam.h | 5 +- .../C_extension/src/integral/integ_process.c | 2 +- pygrt/C_extension/src/integral/ptam.c | 162 +++++++++++------- 5 files changed, 102 insertions(+), 70 deletions(-) diff --git a/docs/source/Advanced/integ_converg/ptam.rst b/docs/source/Advanced/integ_converg/ptam.rst index e3cfa7cc..1981b031 100644 --- a/docs/source/Advanced/integ_converg/ptam.rst +++ b/docs/source/Advanced/integ_converg/ptam.rst @@ -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代码如下， diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index f578392c..ac9e639c 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -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错误 diff --git a/pygrt/C_extension/include/grt/integral/ptam.h b/pygrt/C_extension/include/grt/integral/ptam.h index a8bca4ee..6688c423 100755 --- a/pygrt/C_extension/include/grt/integral/ptam.h +++ b/pygrt/C_extension/include/grt/integral/ptam.h @@ -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 震中距数组 * @@ -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); diff --git a/pygrt/C_extension/src/integral/integ_process.c b/pygrt/C_extension/src/integral/integ_process.c index ee0a0c48..9945af00 100644 --- a/pygrt/C_extension/src/integral/integ_process.c +++ b/pygrt/C_extension/src/integral/integ_process.c @@ -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){ diff --git a/pygrt/C_extension/src/integral/ptam.c b/pygrt/C_extension/src/integral/ptam.c index f869200e..59dca8b4 100755 --- a/pygrt/C_extension/src/integral/ptam.c +++ b/pygrt/C_extension/src/integral/ptam.c @@ -104,31 +104,30 @@ 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); } @@ -136,8 +135,7 @@ static void process_peak_or_trough( * 在输入被积函数的情况下，对不同震源使用峰谷平均法 * * @param[in] ir 震中距索引 - * @param[in] nr 震中距个数 - * @param[in] precoef 积分值系数 + * @param[in] waits_max 当前震中距允许的最大等待采样数 * @param[in] k 波数 * @param[in] dk 波数步长 * @param[in,out] SUM3 被积函数的幅值数组 @@ -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; jjQWV, 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]; + } } // 做缩减序列，赋值最终解 @@ -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