From e88678c17032c531e60ded349cee9ee2e8b24ad3 Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Mon, 12 Jan 2026 22:31:18 +0800 Subject: [PATCH] FEAT: support skip amplitude compensation from imaginary frequency --- docs/source/Module/greenfn.rst | 10 ++++--- pygrt/C_extension/src/dynamic/grt_greenfn.c | 30 ++++++++++++++------- pygrt/pygrn.py | 5 ++-- pygrt/pymod.py | 8 +++--- 4 files changed, 36 insertions(+), 17 deletions(-) diff --git a/docs/source/Module/greenfn.rst b/docs/source/Module/greenfn.rst index 263622d..664171b 100644 --- a/docs/source/Module/greenfn.rst +++ b/docs/source/Module/greenfn.rst @@ -15,7 +15,7 @@ greenfn **grt greenfn** |-M|\ *model* |-D|\ *depsrc/deprcv* -|-N|\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**] +|-N|\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**][**+f**] |-R|\ *file*\|\ *r1,r2,...* |-O|\ *outdir* [ |-H|\ *f1/f2* ] @@ -56,7 +56,7 @@ greenfn .. _-N: -**-N**\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**] +**-N**\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**][**+f**] 采样点数 *nt* 和采样时间间隔 *dt* (secs) ，这将决定计算的最高频。还可设置： + **+w**\ *zeta* - 虚频率系数 :math:`\zeta` [0.8]。 |bouchon1981| 提出在频率上添加一个虚频率偏移， @@ -70,7 +70,11 @@ greenfn 默认情况下，程序会跳过非常低频的几个点以避免引入误差， 详见 :doc:`/Advanced/k_integ/drift/waveform_drift` 的介绍。 - 当时窗长度 nt\*dt 太小“包不住”有效信号，或时窗长度足够但时延(|-E|)不合适，输出的波形会发生混叠， + + **+f** - 不进行由虚频率引起的振幅补偿。 + 默认情况下，程序会对时间域结果乘上 :math:`\exp(\frac{\zeta \pi}{T} t)` 进行振幅补偿， + 但当时窗长度远超有效信号的长度时，这种补偿会放大尾部噪声。 + + **注：** 当时窗长度 nt\*dt 太小“包不住”有效信号，或时窗长度足够但时延(|-E|)不合适，输出的波形会发生混叠， 此时需调整 |-N| 和 |-E| 。 .. include:: explain_-R.rst_ diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index 2ccce70..23b6185 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -72,7 +72,8 @@ typedef struct { real_t wI; ///< 虚频率 zeta*PI/r real_t *freqs; size_t upsample_n; ///< 升采样倍数 - bool keepAllFreq; + bool keepAllFreq; ///< 计算所有频率，不论频率多低 + bool skipImagComps; ///< 跳过虚频率的补偿 } N; /** 输出目录 */ struct { @@ -316,7 +317,7 @@ printf("\n" " : source depth (km).\n" " : receiver depth (km).\n" "\n" -" -N/
[+w][+n][+a] \n" +" -N/
[+w][+n][+a][+f] \n" " : number of points. (NOT requires 2^n).\n" "
: time interval (secs). \n" " +w: define the coefficient of imaginary \n" @@ -328,6 +329,8 @@ printf("\n" " and calculated frequencies stay unchanged.\n" " +a: All frequencies are calculated regardless of\n" " how low the frequency is.\n" +" +f: skip the amplitude compensation from \n" +" imaginary frequency.\n" "\n" " -R,[,...]\n" " Multiple epicentral distance (km), \n" @@ -508,7 +511,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; - // 点数,采样间隔,虚频率 -Nnt/dt[+w][+n][+a] + // 点数,采样间隔,虚频率 -Nnt/dt[+w][+n][+a][+f] case 'N': Ctrl->N.active = true; { @@ -546,6 +549,10 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'a': Ctrl->N.keepAllFreq = true; break; + + case 'f': + Ctrl->N.skipImagComps = true; + break; default: GRTBadOptionError(N, "+%s is not supported.", token); @@ -848,7 +855,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ static void write_one_to_sac( const char *srcname, const char ch, GRT_FFTW_HOLDER *fh, const real_t wI, SACTRACE *sac, const char *s_output_subdir, const char *s_prefix, - const int sgn, const cplx_t *grncplx) + const int sgn, bool skipImagComps, const cplx_t *grncplx) { snprintf(sac->hd.kcmpnm, sizeof(sac->hd.kcmpnm), "%s%s%c", s_prefix, srcname, ch); @@ -879,12 +886,17 @@ static void write_one_to_sac( // 归一化，并处理虚频 // 并转为 SAC 需要的单精度类型 real_t fac, coef; - coef = fh->df * exp(sac->hd.b * wI); - fac = exp(wI*fh->dt); + coef = fh->df; + fac = 1.0; + if (! skipImagComps) { + coef *= exp(sac->hd.b * wI); + fac = exp(wI*fh->dt); + } for(size_t i = 0; i < fh->nt; ++i){ sac->data[i] = fh->w_t[i] * coef; coef *= fac; } + // 以sac文件保存到本地 grt_write_SACTRACE(s_outpath, sac); @@ -1117,11 +1129,11 @@ int greenfn_main(int argc, char **argv) { char ch = GRT_ZRT_CODES[c]; - write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "", sgn, grn[ir][im][c]); + write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "", sgn, Ctrl->N.skipImagComps, grn[ir][im][c]); if(Ctrl->e.active){ - write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "z", sgn*(-1), grn_uiz[ir][im][c]); - write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "r", sgn, grn_uir[ir][im][c]); + write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "z", sgn*(-1), Ctrl->N.skipImagComps, grn_uiz[ir][im][c]); + write_one_to_sac(GRT_SRC_M_NAME_ABBR[im], ch, fh, Ctrl->N.wI, sac, s_output_subdir, "r", sgn , Ctrl->N.skipImagComps, grn_uir[ir][im][c]); } } diff --git a/pygrt/pygrn.py b/pygrt/pygrn.py index 7bc33e2..3dde0c9 100755 --- a/pygrt/pygrn.py +++ b/pygrt/pygrn.py @@ -109,7 +109,7 @@ def plot_response(self): return fig, (ax1, ax2) - def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0): + def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0, skipImagComps:bool=False): ''' Convert the Green's function from the frequency domain to the time domain and return it in the form of :class:`obspy.Trace` @@ -142,7 +142,8 @@ def freq2time(self, T0:float, travtP:float, travtS:float, mult:float=1.0): # 实序列的傅里叶变换 data = irfft(cmlx_grn, nt, norm='backward') * (1/dt) # *(1/dt)和连续傅里叶变换幅值保持一致 # 抵消虚频率的影响 - data *= np.exp((np.arange(0,nt)*dt + T0)*wI) + if not skipImagComps: + data *= np.exp((np.arange(0,nt)*dt + T0)*wI) # 保存sac头段变量 sac.o = 0.0 diff --git a/pygrt/pymod.py b/pygrt/pymod.py index 37398dc..8e82562 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -410,6 +410,7 @@ def compute_grn( converg_method:Union[str,None]=None, delayT0:float=0.0, delayV0:float=0.0, + skipImagComps:bool=False, calc_upar:bool=False, gf_source=['EX', 'VF', 'HF', 'DC'], statsfile:Union[str,None]=None, @@ -439,6 +440,7 @@ def compute_grn( :param safilonTol: precision of Self-Adaptive Filon's Integration Method :param filonCut: The splitting point of DWM and (SA)FIM, k*=/rmax, default is 0 :param converg_method: The method of explicit convergence, you can set "DCM", "PTAM" or "none". Default use "DCM" when abs(depsrc-deprcv) <= 1.0 km + :param skipImagComps: skip the amplitude compensation from imaginary frequency. :param calc_upar: whether calculate the spatial derivatives of displacements. :param gf_source: The source type to be calculated :param statsfile: directory path for saving the statsfile during k integral, used to debug or observe the variations of :math:`F(k,\omega)` and :math:`F(k,\omega)J_m(kr)k` @@ -513,10 +515,10 @@ def compute_grn( continue sgn = -1 if ZRTchs[c]=='Z'=='Z' else 1 - stream.append(pygrnLst[ir][im][c].freq2time(delayT, travtP, travtS, sgn )) + stream.append(pygrnLst[ir][im][c].freq2time(delayT, travtP, travtS, sgn, skipImagComps)) if(calc_upar): - stream.append(pygrnLst_uiz[ir][im][c].freq2time(delayT, travtP, travtS, sgn*(-1) )) - stream.append(pygrnLst_uir[ir][im][c].freq2time(delayT, travtP, travtS, sgn )) + stream.append(pygrnLst_uiz[ir][im][c].freq2time(delayT, travtP, travtS, sgn*(-1), skipImagComps)) + stream.append(pygrnLst_uir[ir][im][c].freq2time(delayT, travtP, travtS, sgn , skipImagComps)) # 在sac头段变量部分