DOC: “波形在漂移”的成因以及解决方法

docs/source/Advanced/index.rst

@@ -1,6 +1,6 @@
 :author: 朱邓达
 :date: 2025-04-23
-:updated_date: 2025-05-05
+:updated_date: 2025-12-05
 
 进阶教程
 =============
@@ -15,3 +15,4 @@
     filon/self_adaptive_filon
     integ_converg/integ_converg
     kernel/kernel
+    waveform_drift/waveform_drift

docs/source/Advanced/waveform_drift/run/milrow

@@ -0,0 +1,9 @@
+0.2  3.4  1.7  2.3
+0.6  3.7  1.9  2.4
+0.5  4.2  2.1  2.4
+0.5  4.6  2.3  2.5
+0.7  4.9  2.8  2.6
+0.5  5.1  2.9  2.7
+6.0  5.9  3.3  2.7
+28   6.9  4.0  2.8
+0.0  8.2  4.7  3.2

docs/source/Advanced/waveform_drift/run/run.py

@@ -0,0 +1,139 @@
+# BEGIN 1
+import numpy as np
+import pygrt
+
+modarr = np.loadtxt("milrow")
+
+pymod = pygrt.PyModel1D(modarr, depsrc=10, deprcv=0.0)
+
+nt = 500
+dt = 10
+
+st_grn = pymod.compute_grn(5000, nt=nt, dt=dt, keepAllFreq=True, statsfile="pygrtstats")[0]
+# END 1
+
+# 仅绘制一个分量做示例
+CHANNEL = "SSR"
+
+# =================================================================
+# 绘制波形
+import matplotlib.pyplot as plt
+from obspy import *
+plt.rcParams.update({
+    "font.sans-serif": "Times New Roman",
+    "mathtext.fontset": "cm"
+})
+
+t = np.arange(nt)*dt
+freqs = np.arange(nt//2 + 1)/(nt*dt)
+
+def plot_waves(tr:Trace):
+    fig, ax = plt.subplots(figsize=(8, 1.5))
+    ax.plot(t, tr.data, 'k-', lw=0.5)
+    ax.grid()
+    ax.text(0.02, 0.9, CHANNEL, transform=ax.transAxes, fontsize=10, ha='left', va='top', bbox=dict(lw=1, fc='w'))
+    ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+    ax.set_xlim(0, nt*dt)
+    ax.set_xlabel("Time (s)", fontsize=12)
+    return fig, ax
+
+fig, ax = plot_waves(st_grn.select(channel=CHANNEL)[0])
+fig.savefig("grn.png", dpi=300, bbox_inches='tight')
+# =================================================================
+
+# =================================================================
+# 回退虚频率补偿
+st_grn2 = st_grn.copy()
+wI = st_grn2[0].stats.sac['user0']
+for tr in st_grn2:
+    tr.data[:] *= np.exp(-t*wI)
+
+fig, ax = plot_waves(st_grn2.select(channel=CHANNEL)[0])
+fig.savefig("grn2.png", dpi=300, bbox_inches='tight')
+# =================================================================
+
+
+# =================================================================
+# 绘制频谱
+from scipy.fft import rfft, irfft
+
+def plot_freqs(tr:Trace):
+    fig, ax = plt.subplots(figsize=(8, 1.5))
+    data = tr.data.copy()
+
+    freqs0 = freqs.copy()
+    freqs0[0] *= 0.1
+    ax.plot(freqs0[:-1], np.abs(rfft(data))[:-1], 'k-', lw=0.8)
+    ax.grid()
+    ax.text(0.02, 0.9, CHANNEL, transform=ax.transAxes, fontsize=10, ha='left', va='top', bbox=dict(lw=1, fc='w'))
+    ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+
+    ax.set_xscale('log')
+    ax.set_yscale('log')
+    ax.set_xlim(1e-4, 0.5/dt)
+    ax.set_xlabel("Frequency (Hz)", fontsize=12)
+
+    return fig, ax
+
+fig, ax = plot_freqs(st_grn.select(channel=CHANNEL)[0])
+fig.savefig("grn_freqs.png", dpi=300, bbox_inches='tight')
+fig, ax = plot_freqs(st_grn2.select(channel=CHANNEL)[0])
+fig.savefig("grn2_freqs.png", dpi=300, bbox_inches='tight')
+# =================================================================
+
+
+# =================================================================
+# 读入核函数
+import glob
+paths = glob.glob("pygrtstats/K_000[0-5]_*")
+paths.sort()
+print(paths)
+
+fig, axs = plt.subplots(len(paths), 1, figsize=(8, len(paths)*1), gridspec_kw=dict(hspace=0), sharex=True)
+for i in range(len(paths)):
+    ax = axs[i]
+    path = paths[i]
+    freq = float(path.split("_")[-1])
+    D = pygrt.utils.read_statsfile(path)
+    ax.plot(D['k'], np.real(D['SS_q']), 'k-', lw=1, label='Real')
+    ax.plot(D['k'], np.imag(D['SS_q']), 'k--', lw=1, label='Imag')
+    ax.text(0.98, 0.9, f"{freq:.1e} Hz", transform=ax.transAxes, fontsize=10, ha='right', va='top', bbox=dict(lw=1, fc='w'))
+    ax.grid()
+    ax.set_xlim(xmax=D['k'][-1])
+    if i == 0:
+        ax.legend(loc='lower center', ncols=2, bbox_to_anchor=(0.5, 1.05))
+
+fig.savefig("kernels.png", dpi=300, bbox_inches='tight')
+# =================================================================
+
+
+
+# =================================================================
+# 跳过频段，重新计算
+st_grn3 = pymod.compute_grn(5000, nt=nt, dt=dt)[0]
+
+srctypes = ['EX', 'VF', 'HF', 'DD', 'DS', 'SS']
+chlst = ['Z', 'R', 'T']
+def plot_all_waves(st_grn:Stream):
+
+    fig = plt.figure(figsize=(14, 2*len(srctypes)))
+    gs = fig.add_gridspec(len(srctypes), 3, hspace=0.3, wspace=0.1)
+    for ik, tr in enumerate(st_grn):
+        channel = tr.stats.channel
+        srctype = channel[:2]
+        ch = channel[-1]
+
+        ax = fig.add_subplot(gs[srctypes.index(srctype), chlst.index(ch)])
+
+        data = tr.data.copy()
+
+        ax.plot(t, data, 'k-', lw=0.6)
+        ax.grid()
+        ax.text(0.1, 0.9, channel, transform=ax.transAxes, ha='center', va='top', bbox=dict(lw=1, fc='w'))
+        ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+
+    return fig
+
+fig = plot_all_waves(st_grn3.copy())
+fig.savefig("grn3.png", dpi=100, bbox_inches='tight')
+# =================================================================

docs/source/Advanced/waveform_drift/run/run.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+
+rm GRN* -rf
+
+# BEGIN 1
+grt greenfn -Mmilrow -D10/0 -N500/10+a -R5000 -OGRN -S
+# END 1
+

docs/source/Advanced/waveform_drift/waveform_drift.rst

@@ -0,0 +1,117 @@
+:author: 朱邓达
+:date: 2025-12-05
+
+如何解决“波形在漂移”的问题？
+================================
+
+问题描述
+--------------------
+
+如下所示，以 MILROW 模型为例，
+计算震源深度为 10 km，台站在地表，震中距为 5000 km 的格林函数，
+采样点数 500，采样间隔 10 秒。仅使用传统的离散波数积分，不使用其他技巧。
+
+脚本下载：:download:`run.py <run/run.py>`
+
+MILROW 模型定义如下，与文档中其它部分一致，
+
+.. literalinclude:: run/milrow
+   :language: text
+
+.. tabs::  
+
+   .. group-tab:: C 
+
+      使用如下命令计算格林函数，
+
+      .. literalinclude:: run/run.sh
+         :language: bash
+         :start-after: BEGIN 1
+         :end-before: END 1
+
+   .. group-tab:: Python
+
+      使用如下 Python 脚本计算格林函数，
+
+      .. literalinclude:: run/run.py
+         :language: python
+         :start-after: BEGIN 1
+         :end-before: END 1
+
+绘制 SSR 波形如下，发现整个波形发生“漂移”。事实上如果绘制所有格林函数会发现它们都发生了不同程度的漂移。
+
+.. figure:: run/grn.png
+   :align: center
+   :width: 100%
+
+问题成因
+-----------------
+
+从形态来看，漂移的形态呈 e 指数，很明显是在时域补偿虚频率产生的效果。
+
+回退虚频率的补偿
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+我们可以将以上波形再乘以 :math:`e^{-\omega_I t}` 以回退虚频率的补偿，观察以下波形形态，
+发现整体发生了偏移，即零频结果不准。
+
+.. figure:: run/grn2.png
+   :align: center
+   :width: 100%
+
+观察频谱
+~~~~~~~~~~~~~~~~
+我们再观察频谱， **发现在极低的频段呈现不正常的高值，而虚频率的补偿会进一步放大这一现象，**
+这表明不止是零频的一个单独问题，而是一系列低频点发生类似问题。
+
++ **虚频率补偿前**
+
+.. figure:: run/grn2_freqs.png
+   :align: center
+   :width: 100%
+
++ **虚频率补偿后**
+
+.. figure:: run/grn_freqs.png
+   :align: center
+   :width: 100%
+
+观察核函数
+~~~~~~~~~~~~~~~~~~~~
+低频段计算发生偏差，显然问题出在核函数。上述在计算格林函数时已经指定了相应的参数以保存核函数。
+如下，我们可以绘制几个低频段的核函数，发现核函数曲线存在很多毛刺，这就导致了积分结果发生偏差。
+
+.. figure:: run/kernels.png
+   :align: center
+   :width: 100%
+
+虚频率的定义
+~~~~~~~~~~~~~~~~~~~~~
+低频的核函数曲线出现毛刺，这是计算机在做数值计算时发生溢出的表现。
+但这个现象在震中距不算太大的情况下并没有出现。事实上，震中距的大小是间接原因，
+根本原因是时窗长度。
+
+根据 Bouchon (1981) 的设计，虚频率定义为 :math:`\omega_I = \zeta \dfrac{\pi}{T}` ，
+其中 :math:`T` 为时窗长度。虚频率的加入其中一个目的就是使核函数的极点略微偏离实轴，
+从而可以进行常规的波数积分。但当震中距增大，模型速度减小，计算时就需要更宽的时窗来容纳所有有效信号以避免混叠，
+代价是虚频率减小，核函数的极点逐渐靠近实轴，导致核函数计算的不稳定性风险增加。
+
+核函数有多个极点，且随频率变化，而从实际计算来看，虚频率的减小对极低频的核函数计算影响很大。
+这是因为在核函数的计算中，无时延的广义 R/T 系数矩阵仅和相速度有关，当引入虚频率后，对应的复数形式的相速度定义为
+:math:`\tilde{c}=\dfrac{\omega - j \omega_I}{k}` ，当实频率 :math:`\omega` 和虚频率 :math:`\omega_I` 
+都很小时， :math:`\tilde{c} \rightarrow 0` ，这就导致核函数计算不稳定，也就导致低频结果有偏差，整体波形在漂。
+
+解决方法
+-----------------
+既然这是方法本身的局限性导致这些低频的几个点计算不出精确结果，不如直接放弃它们。
+
+在以上计算格林函数时，可以看到为了复现上述问题，我添加了一些参数 （ ``-N+a`` 和 ``keepAllFreq`` ），
+这是因为从 0.13.0 版本起， **程序默认会根据复频率的大小跳过低频的计算。**
+具体而言，程序不计算 :math:`|\omega - \omega_I| < 0.1` 的频谱，并会在终端打印对应的Warning。
+这本质就是自动进行了高通滤波。
+如果觉得阈值不合适或需要进行数值实验，可加上上述参数以计算所有频率，然后再手动设置频段。
+
+以下是跳过这些频段后的波形（计算时不使用 ``-N+a`` 和 ``keepAllFreq`` ）。
+
+.. figure:: run/grn3.png
+   :align: center
+   :width: 100%

docs/source/Module/greenfn.rst

@@ -1,6 +1,6 @@
 :author: 朱邓达
 :date: 2025-09-22
-:updated_date: 2025-11-24
+:updated_date: 2025-12-05
 
 .. include:: common_OPTs.rst_
 
@@ -16,7 +16,7 @@ greenfn
 **grt greenfn** 
 |-M|\ *model*
 |-D|\ *depsrc/deprcv*
-|-N|\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*]
+|-N|\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**]
 |-R|\ *file*\|\ *r1,r2,...*
 |-O|\ *outdir*
 [ |-H|\ *f1/f2* ]
@@ -56,7 +56,7 @@ greenfn
 
 .. _-N:
 
-**-N**\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*]
+**-N**\ *nt/dt*\ [**+w**\ *zeta*][**+n**\ *fac*][**+a**]
     采样点数 *nt* 和采样时间间隔 *dt* (secs) ，这将决定计算的最高频。还可设置：
 
     + **+w**\ *zeta* - 虚频率系数 :math:`\zeta` [0.8]。Bouchon (1981) 提出在频率上添加一个虚频率偏移， 
@@ -66,6 +66,10 @@ greenfn
     + **+n**\ *fac* - 频率域插值倍数[1]。即在做逆傅里叶变换时，在频域上最高频后补零，
       相当于 *nt* ← *nt* \* *fac* , *dt* ← *dt* / *fac* ，使计算的波形更平滑。
 
+    + **+a** - 计算所有频点，不论频率多低。除非做数值实验，否则不建议使用该选项。
+      默认情况下，程序会跳过非常低频的几个点以避免引入误差，
+      详见 :doc:`/Advanced/waveform_drift/waveform_drift` 的介绍。
+
     当时窗长度 nt\*dt 太小“包不住”有效信号，或时窗长度足够但时延(|-E|)不合适，输出的波形会发生混叠，
     此时需调整 |-N| 和 |-E| 。
 

docs/source/run_all.sh

@@ -69,4 +69,10 @@ if [[ $1 == '8' || $1 == '' ]]; then
     python lamb1_plot_time.py
     python lamb1_plot_freq_time.py
     cd -
+fi
+
+if [[ $1 == '9' || $1 == '' ]]; then
+    cd Advanced/waveform_drift/run
+    python run.py
+    cd -
 fi

pygrt/C_extension/src/dynamic/grt_greenfn.c

@@ -281,7 +281,7 @@ printf("\n"
 "Usage:\n"
 "----------------------------------------------------------------\n"
 "    grt greenfn -M<model> -D<depsrc>/<deprcv> \n"
-"        -N<nt>/<dt>[+w<zeta>][+n<fac>] \n"
+"        -N<nt>/<dt>[+w<zeta>][+n<fac>][+a] \n"
 "        -R<r1>,<r2>[,...]     -O<outdir>     [-H<f1>/<f2>] \n"
 "        [-L<length>]        [-E<t0>[/<v0>]] \n" 
 "        [-K[+k<k0>][+s<ampk>][+e<keps>][+v<vmin>]]\n"