greenfn
- 简介:
使用广义反射透射系数矩阵法 + 离散波数法计算动态全波解格林函数
语法
grt greenfn -Mmodel -Ddepsrc/deprcv -Nnt/dt[+wzeta][+nfac][+a] -Rfile|r1,r2,… -Ooutdir [ -Hf1/f2 ] [ -Largs ] [ -K[+kk0][+sampk][+ekeps][+vvmin] ] [ -Et0[/v0] ] [ -Pnthreads ] [ -Ge|v|h|s ] [ -S[i1,i2,…] ] [ -e ] [ -s ] [ -h ]
描述
greenfn 模块计算每个震中距 r 的格林函数波形保存路径为
{outdir}/{model}_{depsrc}_{deprcv}_{r}/{stype}.sac ，
其中支持的格林函数类型 stype 详见 格林函数分类 。
执行 greenfn 模块的原始命令会附加式地保存到 {outdir}/command 文件。
不同震源的格林函数单位为：
爆炸源： \(10^{-20} \, \frac{\text{cm}}{\text{dyne} \cdot \text{cm}}\)
单力源： \(10^{-15} \, \frac{\text{cm}}{\text{dyne}}\)
剪切源： \(10^{-20} \, \frac{\text{cm}}{\text{dyne} \cdot \text{cm}}\)
必选选项
- -Mmodel
半空间层状模型文件的路径。模型格式如下:
Thickness(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) [Qp] [Qs]
后两列的 Qp, Qs 可省略。
如果首列首行的值为正数，则首列表示 每层厚度 (km) ，最后一层为半空间，其厚度设置为 0 即可（不会被使用）。
如果首列首行的值为 0，则首列表示 每层层顶深度 (km) 。
如果某层设置 Vs = 0，则该层为液体层。
- -Ddepsrc/deprcv
震源深度 depsrc (km) 和台站深度 deprcv (km)。 如果是在使用波数积分法求解格林函数，则当二者深度差小于 1 km 时，自动使用快速收敛算法。
- -Nnt/dt[+wzeta][+nfac][+a]
采样点数 nt 和采样时间间隔 dt (secs) ，这将决定计算的最高频。还可设置：
+wzeta - 虚频率系数 \(\zeta\) [0.8]。Bouchon (1981) 提出在频率上添加一个虚频率偏移， \(\omega \leftarrow \omega - i \zeta \dfrac{\pi}{T}\) ， 其中 \(T\) 为时窗长度 nt*dt，以使波数积分适当偏移实轴上的极点。
+nfac - 频率域插值倍数[1]。即在做逆傅里叶变换时，在频域上最高频后补零， 相当于 nt ← nt * fac , dt ← dt / fac ，使计算的波形更平滑。
+a - 计算所有频点，不论频率多低。除非做数值实验，否则不建议使用该选项。 默认情况下，程序会跳过非常低频的几个点以避免引入误差， 详见 如何解决“波形在漂移”的问题？ 的介绍。
当时窗长度 nt*dt 太小“包不住”有效信号，或时窗长度足够但时延(-E)不合适，输出的波形会发生混叠， 此时需调整 -N 和 -E 。
- -Rfile|r1,r2,…
震中距 (km)。可以传入逗号分隔的多个震中距，或者仅含有一列震中距的文件 file 。
- -Ooutdir
输出目录名，不存在会自动新建。若存在同名保存路径则直接覆盖已有结果。
可选选项
- -Hf1/f2
待计算的频率范围（闭区间）(Hz) [默认计算由 -D 指定的全部范围]。 f1 和 f2 分别为最小最大频率，若取 -1 则分别对应低通和高通。
- -Largs
控制波数积分步长 \(\Delta k =\frac{2\pi}{\textit{<length>}\cdot R_{\text{max}}}\) 。 其中 \(R_{\text{max}}\) 为多个震中距的最大值。 具体设置有以下几种情况：
若不使用 -L ，则积分步长中的 length 参数自动确定，其准则详见 Bouchon (1981) 或代码源文件 。
-Llength ，手动设置 length ，例如 -L20 。
-Lllength/Flength/Fcut ，手动设置三个参数， 例如 -Ll20/10/10 。 此时原本波数积分区间 \([0, k_{\text{max}}]\) 被划分为两段， \([0, k^*]\) 和 \([k^*, k_{\text{max}}]\) ， 其中 \(k^*=\frac{\textit{<Fcut>}}{R_{\text{max}}}\) 。 前段仍使用离散波数法求解积分，步长由 length 控制， 后段则使用 固定间隔的Filon积分法 , 对应步长由 Flength 控制 （仍为上述 \(\Delta k\) 公式）。
-Lalength/Ftol/Fcut ，手动设置三个参数，例如 -La20/1e-3/0 。 与上个情况类似，但后段使用 自适应Filon积分法 。 这在大震中距情况下能显著减少计算时间。
备注
固定间隔和自适应 Filon 积分不计算近场项，即 计算动态格林函数 中介绍的 \(p=1\) 对应的积分项。
- -K[+kk0][+sampk][+ekeps][+vvmin]
控制波数积分上限
\[k_{\text{max}} = \sqrt{ k_0 \cdot \dfrac{\pi}{\Delta h} + \textit{<ampk>} \cdot \left(\dfrac{\omega}{v_{\text{min}}}\right)^2}\]+kk0 - 控制零频的积分上限 [5.0]，其中深度差 \(\Delta h = \max(|z_s - z_r|, 1.0)\) 。
+sampk - 放大倍数 [1.15] 。
+ekeps - 用于判断提前结束波数积分的收敛精度[0.0, 默认不使用]，详见 Yao and Harkrider (1983) 和 控制波数积分 。
+vvmin - 参考最小速度，默认 \(\max{(\min\limits_{i} (\alpha_i \cup \beta_i), 0.1)}\) 。
只要设置了 vmin ，不论正负， 是否启用快速收敛算法不再受 \(\Delta h\) 控制 （见 -D ）， 相当于改为由 +vvmin 手动控制。
当 vmin 为负数时，使用快速收敛算法（这个用法更多是用于 debug 测试）。
- -Et0[/v0]
增加时延以调整波形起始时刻。输出的波形默认起始时刻为 0 时刻，即发震时刻。 对于大震中距的记录，P 波初至走时较大，此时可以适当调整波形起始时刻， 使 -N 定义的时窗长度能框住整段波形， 而不必使用 -N 定义从发震时刻开始的超宽时窗，增加不必要的计算量。支持以下用法：
t0 - 起始时刻调整为 \(t_0\) (secs) 。
t0/v0 - 每个震中距 \(r\) 的记录起始时刻不同，调整为 \(t_0 + \frac{r}{v_0}\) ， v0 为参考速度 (km/s) 。
- -Pnthreads
多线程数，默认使用全部线程。
- -Ge|v|h|s
控制输出的震源类型 [evhs]，不会减少计算量。
e - 爆炸源
v - 垂直力源
h - 水平力源
s - 剪切源
- -S[i1,i2,…]
指定频率索引值，输出对应频率下波数积分过程中的核函数文件，多个索引值用逗号分隔，文件保存目录为
{outdir}_grtstats。若使用 -S 则保存所有频点。 关于文件格式及其读取详见 积分收敛性 。
- -e
在计算格林函数的同时，也计算其空间导数 \(\dfrac{\partial (u_z, u_r, u_\theta)}{\partial (z, r)}\) 。 偏导对应在文件名/变量名开头添加了
z和r。 关于 \(\theta\) 的偏导与方向因子有关，这将由 syn 和 static_syn 模块计算。不同震源的格林函数空间导数单位为:
爆炸源： \(10^{-25} \, /\text{dyne} \cdot \text{cm}\)
单力源： \(10^{-20} \, /\text{dyne}\)
剪切源： \(10^{-25} \, /\text{dyne} \cdot \text{cm}\)
- -s
静默输出，不在终端打印参数和进度条。
- -h
打印帮助文档。
示例
详见教程：