diff --git a/.github/workflows/testbuild.yml b/.github/workflows/testbuild.yml index fc34a9e5..d561fea4 100644 --- a/.github/workflows/testbuild.yml +++ b/.github/workflows/testbuild.yml @@ -236,19 +236,20 @@ jobs: # don't use $(pwd) to get abspath, not valid on Windows echo "EXAMPLE_COPY_PATH=${{ env.PACK_NAME }}/test_tmp" >> $GITHUB_ENV - - name: Test 1 compare_results - working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compare_results + - name: Test compare_c_py + working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compare_c_py run: | + # 检查C程序和Python程序的计算结果是否一致 chmod +x *.sh - ./run_milrow_grt.sh - python plot_cps_pygrt.py + ./run_grt.sh + python run_pygrt.py - - name: Test 2 site_effect - working-directory: ${{ env.EXAMPLE_COPY_PATH }}/site_effect - run: | - chmod +x *.sh - ./run1.sh - python run2.py + # - name: Test 2 site_effect + # working-directory: ${{ env.EXAMPLE_COPY_PATH }}/site_effect + # run: | + # chmod +x *.sh + # ./run1.sh + # python run2.py # - name: Test 3 multi_traces # working-directory: ${{ env.EXAMPLE_COPY_PATH }}/multi_traces @@ -258,12 +259,12 @@ jobs: # ./run1.sh # continue-on-error: true # 即使失败，仍然标记为成功 - - name: Test 4 lamb_problem - working-directory: ${{ env.EXAMPLE_COPY_PATH }}/lamb_problem - run: | - chmod +x *.sh - ./run1.sh - python run2.py + # - name: Test 4 lamb_problem + # working-directory: ${{ env.EXAMPLE_COPY_PATH }}/lamb_problem + # run: | + # chmod +x *.sh + # ./run1.sh + # python run2.py # - name: Test 5 far_field # working-directory: ${{ env.EXAMPLE_COPY_PATH }}/far_field @@ -274,24 +275,24 @@ jobs: # python plot_compare_pygrt.py # continue-on-error: true # 即使失败，仍然标记为成功 - - name: Test signal convolution - working-directory: ${{ env.EXAMPLE_COPY_PATH }}/convolve_signals - run: | - chmod +x *.sh - ./run_grt.sh - python run_pygrt.py + # - name: Test signal convolution + # working-directory: ${{ env.EXAMPLE_COPY_PATH }}/convolve_signals + # run: | + # chmod +x *.sh + # ./run_grt.sh + # python run_pygrt.py - - name: Test strain stress (dynamic) - working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/dynamic - run: | - chmod +x *.sh - ./run_grt.sh + # - name: Test strain stress (dynamic) + # working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/dynamic + # run: | + # chmod +x *.sh + # ./run_grt.sh - - name: Test static displacement (static) - working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/static - run: | - chmod +x *.sh - ./run_stgrt.sh + # - name: Test static displacement (static) + # working-directory: ${{ env.EXAMPLE_COPY_PATH }}/compute_strain_stress/static + # run: | + # chmod +x *.sh + # ./run_stgrt.sh - name: Remove the test files run: | diff --git a/example/compare_c_py/README b/example/compare_c_py/README new file mode 100644 index 00000000..be13a40b --- /dev/null +++ b/example/compare_c_py/README @@ -0,0 +1 @@ +Compare results from C API and Python API. \ No newline at end of file diff --git a/example/compare_c_py/milrow b/example/compare_c_py/milrow new file mode 100644 index 00000000..a8a74b18 --- /dev/null +++ b/example/compare_c_py/milrow @@ -0,0 +1,9 @@ +0.2 3.4 1.7 2.3 9e10 9e10 +0.6 3.7 1.9 2.4 9e10 9e10 +0.5 4.2 2.1 2.4 9e10 9e10 +0.5 4.6 2.3 2.5 9e10 9e10 +0.7 4.9 2.8 2.6 9e10 9e10 +0.5 5.1 2.9 2.7 9e10 9e10 +6.0 5.9 3.3 2.7 9e10 9e10 +28. 6.9 4.0 2.8 9e10 9e10 +-0.1 8.2 4.7 3.2 9e10 9e10 diff --git a/example/compare_c_py/run_grt.sh b/example/compare_c_py/run_grt.sh new file mode 100755 index 00000000..52a1693d --- /dev/null +++ b/example/compare_c_py/run_grt.sh @@ -0,0 +1,90 @@ +#!/bin/bash + + +dist=10 +depsrc=2 +deprcv=3.3 + +nt=1024 +dt=0.01 + +modname="milrow" +out="GRN" + + +#-------------------------- Dynamic ----------------------------------------- +grt -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${dist} -e + +# convolve different signals +G=$out/${modname}_${depsrc}_${deprcv}_${dist} +P=cout +S=1e24 +az=39.2 +for N in "-N" "" ; do +grt.syn -G$G -Osyn_exp -P$P -A$az -S$S -Dt/0.2/0.2/0.4 -e $N +grt.strain syn_exp/$P +grt.stress syn_exp/$P + +fn=2 +fe=-1 +fz=4 +grt.syn -G$G -Osyn_sf -P$P -A$az -S$S -F$fn/$fe/$fz -Dt/0.1/0.3/0.6 -e $N +grt.strain syn_sf/$P +grt.stress syn_sf/$P + +stk=77 +dip=88 +rak=99 +grt.syn -G$G -Osyn_dc -P$P -A$az -S$S -M$stk/$dip/$rak -Dp/0.6 -e $N +grt.strain syn_dc/$P +grt.stress syn_dc/$P + +M11=1 +M12=-2 +M13=-5 +M22=0.5 +M23=3 +M33=1.2 +grt.syn -G$G -Osyn_mt -P$P -A$az -S$S -T$M11/$M12/$M13/$M22/$M23/$M33 -Dr/3 -e $N +grt.strain syn_mt/$P +grt.stress syn_mt/$P +done + + +#-------------------------- Static ----------------------------------------- +x1=-3 +x2=3 +nx=11 + +y1=-4 +y2=4 +ny=16 + +mkdir -p static +cd static +stgrt -M../${modname} -D${depsrc}/${deprcv} -X$x1/$x2/$nx -Y$y1/$y2/$ny -e > grn + +for N in "-N" "" ; do +stgrt.syn -S$S -e $N < grn > stsyn_exp$N +stgrt.strain < stsyn_exp$N > strain_exp$N +stgrt.stress < stsyn_exp$N > stress_exp$N + +stgrt.syn -S$S -F$fn/$fe/$fz -e $N < grn > stsyn_sf$N +stgrt.strain < stsyn_sf$N > strain_sf$N +stgrt.stress < stsyn_sf$N > stress_sf$N + +stgrt.syn -S$S -M$stk/$dip/$rak -e $N < grn > stsyn_dc$N +stgrt.strain < stsyn_dc$N > strain_dc$N +stgrt.stress < stsyn_dc$N > stress_dc$N + +stgrt.syn -S$S -T$M11/$M12/$M13/$M22/$M23/$M33 -e $N < grn > stsyn_mt$N +stgrt.strain < stsyn_mt$N > strain_mt$N +stgrt.stress < stsyn_mt$N > stress_mt$N +done + +cd - + +for p in $(ls static/*); do + echo "----------------------------- $p -------------------------" + cat $p +done \ No newline at end of file diff --git a/example/compare_c_py/run_pygrt.py b/example/compare_c_py/run_pygrt.py new file mode 100644 index 00000000..5a7a35c9 --- /dev/null +++ b/example/compare_c_py/run_pygrt.py @@ -0,0 +1,201 @@ +import numpy as np +import pygrt +from obspy import * + + +def compare3(st_py:Stream, c_prefix:str, ZNE:bool=False, dim2:bool=False): + """return average relative error""" + + if ZNE: + pattern = "[ZNE]" + else: + pattern = "[ZRT]" + + if dim2: + pattern = pattern*2 + + st_c = read(f"{c_prefix}{pattern}.sac") + + error = 0.0 + nerr = 0 + for tr_c in st_c: + tr_py = st_py.select(channel=tr_c.stats.channel)[0] + rerr = np.mean(np.abs(tr_c.data - tr_py.data) / np.abs(tr_c.data)) + if np.isnan(rerr) or np.isinf(rerr): + rerr = 0.0 + error += rerr + nerr += 1 + + return error/nerr + + +def static_compare3(resDct:dict, c_prefix:str, ZNE:bool=False, dim2:bool=False): + """return average relative error""" + midname = c_prefix.split("_")[-1][:2].upper() + + if ZNE: + chs2 = ['Z', 'N', 'E'] + else: + chs2 = ['Z', 'R', 'T'] + + chs1 = chs2.copy() + + nx = len(resDct['_xarr']) + ny = len(resDct['_yarr']) + + c_res = np.loadtxt(c_prefix) + c_resDct = {} + if not dim2: + for i, c in enumerate(chs2): + c_resDct[f"{midname}{c}"] = c_res[:,i+2].reshape((nx, ny), order='F') + else: + chs1 = [] + for i1, c1 in enumerate(chs2): + for i2, c2 in enumerate(chs2): + if i2 < i1: + continue + chs1.append(f"{c1}{c2}") + + for i, cc in enumerate(chs1): + c_resDct[cc] = c_res[:,i+2].reshape((nx, ny), order='F') + + error = 0.0 + nerr = 0 + + for k in c_resDct.keys(): + val1 = resDct[k] + val2 = c_resDct[k] + rerr = np.mean(np.abs(val1 - val2) / np.abs(val2)) + if np.isnan(rerr) or np.isinf(rerr): + rerr = 0.0 + error += rerr + nerr += 1 + + return error/nerr + + + +dist=10 +depsrc=2 +deprcv=3.3 + +nt=1024 +dt=0.01 + +modname="milrow" + +modarr = np.loadtxt(modname) + +pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) + +AVGRERR = [] + +#-------------------------- Dynamic ----------------------------------------- +# compute green functions +st_grn = pymod.compute_grn(dist, nt, dt, calc_upar=True)[0] + +S=1e24 +az=39.2 + +fn=2 +fe=-1 +fz=4 + +stk=77 +dip=88 +rak=99 + +M11=1 +M12=-2 +M13=-5 +M22=0.5 +M23=3 +M33=1.2 + +for ZNE in [False, True]: + # synthetic + + st = pygrt.utils.gen_syn_from_gf_EXP(st_grn, S, az, ZNE=ZNE, calc_upar=True) + sigs = pygrt.sigs.gen_triangle_wave(0.4, dt) + pygrt.utils.stream_convolve(st, sigs) + AVGRERR.append(compare3(st, "syn_exp/cout", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(st) + ststress = pygrt.utils.compute_stress(st) + AVGRERR.append(compare3(ststrain, "syn_exp/cout.strain.", ZNE=ZNE, dim2=True)) + AVGRERR.append(compare3(ststress, "syn_exp/cout.stress.", ZNE=ZNE, dim2=True)) + + + st = pygrt.utils.gen_syn_from_gf_SF(st_grn, S, fn, fe, fz, az, ZNE=ZNE, calc_upar=True) + sigs = pygrt.sigs.gen_trap_wave(0.1, 0.3, 0.6, dt) + pygrt.utils.stream_convolve(st, sigs) + AVGRERR.append(compare3(st, "syn_sf/cout", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(st) + ststress = pygrt.utils.compute_stress(st) + AVGRERR.append(compare3(ststrain, "syn_sf/cout.strain.", ZNE=ZNE, dim2=True)) + AVGRERR.append(compare3(ststress, "syn_sf/cout.stress.", ZNE=ZNE, dim2=True)) + + st = pygrt.utils.gen_syn_from_gf_DC(st_grn, S, stk, dip, rak, az, ZNE=ZNE, calc_upar=True) + sigs = pygrt.sigs.gen_parabola_wave(0.6, dt) + pygrt.utils.stream_convolve(st, sigs) + AVGRERR.append(compare3(st, "syn_dc/cout", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(st) + ststress = pygrt.utils.compute_stress(st) + AVGRERR.append(compare3(ststrain, "syn_dc/cout.strain.", ZNE=ZNE, dim2=True)) + AVGRERR.append(compare3(ststress, "syn_dc/cout.stress.", ZNE=ZNE, dim2=True)) + + st = pygrt.utils.gen_syn_from_gf_MT(st_grn, S, [M11,M12,M13,M22,M23,M33], az, ZNE=ZNE, calc_upar=True) + sigs = pygrt.sigs.gen_ricker_wave(3, dt) + pygrt.utils.stream_convolve(st, sigs) + AVGRERR.append(compare3(st, "syn_mt/cout", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(st) + ststress = pygrt.utils.compute_stress(st) + AVGRERR.append(compare3(ststrain, "syn_mt/cout.strain.", ZNE=ZNE, dim2=True)) + AVGRERR.append(compare3(ststress, "syn_mt/cout.stress.", ZNE=ZNE, dim2=True)) + + +#-------------------------- Static ----------------------------------------- +xarr = np.linspace(-3, 3, 11) +yarr = np.linspace(-4, 4, 16) +static_grn = pymod.compute_static_grn(xarr, yarr, calc_upar=True) + +for ZNE in [False, True]: + suffix = "-N" if ZNE else "" + static_syn = pygrt.utils.gen_syn_from_gf_EXP(static_grn, S, ZNE=ZNE, calc_upar=True) + AVGRERR.append(static_compare3(static_syn, f"static/stsyn_exp{suffix}", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + AVGRERR.append(static_compare3(ststrain, f"static/strain_exp{suffix}", ZNE=ZNE, dim2=True)) + AVGRERR.append(static_compare3(ststress, f"static/stress_exp{suffix}", ZNE=ZNE, dim2=True)) + print(static_syn, ststrain, ststress) + + static_syn = pygrt.utils.gen_syn_from_gf_SF(static_grn, S, fn, fe, fz, ZNE=ZNE, calc_upar=True) + AVGRERR.append(static_compare3(static_syn, f"static/stsyn_sf{suffix}", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + AVGRERR.append(static_compare3(ststrain, f"static/strain_sf{suffix}", ZNE=ZNE, dim2=True)) + AVGRERR.append(static_compare3(ststress, f"static/stress_sf{suffix}", ZNE=ZNE, dim2=True)) + print(static_syn, ststrain, ststress) + + static_syn = pygrt.utils.gen_syn_from_gf_DC(static_grn, S, stk, dip, rak, ZNE=ZNE, calc_upar=True) + AVGRERR.append(static_compare3(static_syn, f"static/stsyn_dc{suffix}", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + AVGRERR.append(static_compare3(ststrain, f"static/strain_dc{suffix}", ZNE=ZNE, dim2=True)) + AVGRERR.append(static_compare3(ststress, f"static/stress_dc{suffix}", ZNE=ZNE, dim2=True)) + print(static_syn, ststrain, ststress) + + static_syn = pygrt.utils.gen_syn_from_gf_MT(static_grn, S, [M11,M12,M13,M22,M23,M33], ZNE=ZNE, calc_upar=True) + AVGRERR.append(static_compare3(static_syn, f"static/stsyn_mt{suffix}", ZNE=ZNE)) + ststrain = pygrt.utils.compute_strain(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + AVGRERR.append(static_compare3(ststrain, f"static/strain_mt{suffix}", ZNE=ZNE, dim2=True)) + AVGRERR.append(static_compare3(ststress, f"static/stress_mt{suffix}", ZNE=ZNE, dim2=True)) + print(static_syn, ststrain, ststress) + + +AVGRERR = np.array(AVGRERR) +print(AVGRERR) +print(np.mean(AVGRERR), np.min(AVGRERR), np.max(AVGRERR)) +if np.mean(AVGRERR) > 0.01: + raise ValueError + diff --git a/example/convolve_signals/run_pygrt.py b/example/convolve_signals/run_pygrt.py index 1ef92f38..670990b0 100644 --- a/example/convolve_signals/run_pygrt.py +++ b/example/convolve_signals/run_pygrt.py @@ -15,7 +15,7 @@ pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) # compute green functions -st_grn = pymod.gen_gf_spectra(dist, nt, dt)[0] +st_grn = pymod.compute_grn(dist, nt, dt)[0] # synthetic S=1e24 diff --git a/pygrt/C_extension/include/common/attenuation.h b/pygrt/C_extension/include/common/attenuation.h index 9e3d270c..2d3bdfd6 100755 --- a/pygrt/C_extension/include/common/attenuation.h +++ b/pygrt/C_extension/include/common/attenuation.h @@ -22,4 +22,9 @@ * @param omega (in)复数频率\f$ \tilde{\omega} =\omega - i\zeta \f$ * @return atncoef 系数因子，作用在 \f$ k=\omega / c(\omega)\f$的计算 */ -MYCOMPLEX attenuation_law(MYREAL Qinv, MYCOMPLEX omega); \ No newline at end of file +MYCOMPLEX attenuation_law(MYREAL Qinv, MYCOMPLEX omega); + +/** + * attenuation_law函数在python中被调用的版本，长度2的数组分别表示复数的实部和虚部 + */ +void py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]); \ No newline at end of file diff --git a/pygrt/C_extension/include/common/const.h b/pygrt/C_extension/include/common/const.h index 3c325363..2709c1b6 100755 --- a/pygrt/C_extension/include/common/const.h +++ b/pygrt/C_extension/include/common/const.h @@ -120,4 +120,7 @@ typedef int MYINT; ///< 整数 #define INIT_C_IDENTITY_2x2_MATRIX {{CONE, CZERO}, {CZERO, CONE}} ///< 初始化复数单位阵 #define MIN_DEPTH_GAP_SRC_RCV 1.0 ///< 震源和台站的最小深度差（不做绝对限制，仅用于参考波数积分上限） -#define GCC_ALWAYS_INLINE __attribute__((always_inline)) ///< gcc编译器不改动内联函数 \ No newline at end of file +#define GCC_ALWAYS_INLINE __attribute__((always_inline)) ///< gcc编译器不改动内联函数 + +#define GRT_STATIC_STRING_FMT "%18s" ///< stgrt系列程序字符串输出格式 +#define GRT_STATIC_REAL_FMT "%18.8e" ///< stgrt系列程序浮点数输出格式 \ No newline at end of file diff --git a/pygrt/C_extension/src/common/attenuation.c b/pygrt/C_extension/src/common/attenuation.c index 2b8ad23d..20ff9edc 100755 --- a/pygrt/C_extension/src/common/attenuation.c +++ b/pygrt/C_extension/src/common/attenuation.c @@ -15,4 +15,12 @@ MYCOMPLEX attenuation_law(MYREAL Qinv, MYCOMPLEX omega){ return RONE + Qinv/PI * CLOG(omega/PI2) + RHALF*Qinv*I; // return RONE; +} + +void py_attenuation_law(MYREAL Qinv, MYREAL omg[2], MYREAL atte[2]){ + // 用于在python中调用attenuation_law + MYCOMPLEX omega = omg[0] + I*omg[1]; + MYCOMPLEX atte0 = attenuation_law(Qinv, omega); + atte[0] = CREAL(atte0); + atte[1] = CIMAG(atte0); } \ No newline at end of file diff --git a/pygrt/C_extension/src/static/stgrt_main.c b/pygrt/C_extension/src/static/stgrt_main.c index b57f829c..1a732e9f 100644 --- a/pygrt/C_extension/src/static/stgrt_main.c +++ b/pygrt/C_extension/src/static/stgrt_main.c @@ -490,8 +490,8 @@ int main(int argc, char **argv){ MYREAL rcv_rho = pymod->Rho[pymod->ircv]; // 输出物性参数 - fprintf(stdout, "# %15.5e %15.5e %15.5e\n", src_va, src_vb, src_rho); - fprintf(stdout, "# %15.5e %15.5e %15.5e\n", rcv_va, rcv_vb, rcv_rho); + fprintf(stdout, "# "GRT_STATIC_REAL_FMT" "GRT_STATIC_REAL_FMT" "GRT_STATIC_REAL_FMT"\n", src_va, src_vb, src_rho); + fprintf(stdout, "# "GRT_STATIC_REAL_FMT" "GRT_STATIC_REAL_FMT" "GRT_STATIC_REAL_FMT"\n", rcv_va, rcv_vb, rcv_rho); // 定义标题数组 const char *titles[15] = { @@ -506,11 +506,13 @@ int main(int argc, char **argv){ }; // 输出标题 - fprintf(stdout, "#%14s", "X(km)"); - fprintf(stdout, "%15s", "Y(km)"); - for(int i=0; i<15; ++i) fprintf(stdout, "%15s", titles[i]); + char XX[20]; + sprintf(XX, GRT_STATIC_STRING_FMT, "X(km)"); XX[0]='#'; + fprintf(stdout, "%s", XX); + fprintf(stdout, GRT_STATIC_STRING_FMT, "Y(km)"); + for(int i=0; i<15; ++i) fprintf(stdout, GRT_STATIC_STRING_FMT, titles[i]); if(calc_upar) { - for(int i=0; i<30; ++i) fprintf(stdout, "%15s", upar_titles[i]); + for(int i=0; i<30; ++i) fprintf(stdout, GRT_STATIC_STRING_FMT, upar_titles[i]); } fprintf(stdout, "\n"); @@ -518,7 +520,7 @@ int main(int argc, char **argv){ for(int iy=0; iy` 在频率上添加微小虚部， 更多测试见 :ref:`(张海明, 2021) ` :param vmin_ref: 最小参考速度，默认vmin=max(minimum velocity, 0.1)，用于定义波数积分上限，小于0则在达到积分上限后使用峰谷平均法 - （默认当震源和场点深度差<=0.5km时自动使用峰谷平均法） + （默认当震源和场点深度差<=1km时自动使用峰谷平均法） :param keps: 波数k积分收敛条件，见 :ref:`(Yao and Harkrider, 1983) ` :ref:`(初稿) `， 为负数代表不提前判断收敛，按照波数积分上限进行积分 :param ampk: 影响波数k积分上限的系数，见下方 :param iwk0: k0是否取随频率变化的线性关系，即 :math:` k_{0} = k_{0} * f/f_{max}` :param k0: 波数k积分的上限 :math:`\tilde{k_{max}}=\sqrt{(k_{0}*\pi/hs)^2 + (ampk*w/vmin_{ref})^2} ` , 波数k积分循环必须退出, hs=max(震源和台站深度差,1.0) :param Length: 定义波数k积分的间隔 `dk=2\pi / (L*rmax)`, 选取要求见 :ref:`(Bouchon, 1981) ` - :ref:`(张海明, 2021) `，默认自动选择 + :ref:`(张海明, 2021) `，默认自动选择；负数表示使用Filon积分 :param calc_upar: 是否计算位移u的空间导数 :param gf_source: 待计算的震源类型 :param statsfile: 波数k积分（包括Filon积分和峰谷平均法）的过程记录文件，常用于debug或者观察积分过程中 :math:`F(k,\omega)` 和 :math:`F(k,\omega)J_m(kr)k` 的变化 @@ -253,7 +259,7 @@ def gen_gf_spectra( :param print_runtime: 是否打印运行时间 :return: - - **gfspecList** - 列表，每个元素以字典的形式存有每个震中距的格林函数频谱(:class:`pygrt.pygrn.PyGreenFunction` ) 实例) + - **dataLst** - 列表，每个元素为 :class:`obspy.Stream` 类型 ) ''' @@ -325,7 +331,7 @@ def gen_gf_spectra( # 参考最小速度 if vmin_ref == 0.0: vmin_ref = max(self.vmin, 0.1) - if abs(depsrc - deprcv) <= 0.5: + if abs(depsrc - deprcv) <= 1.0: vmin_ref = - abs(vmin_ref) # 自动使用PTAM @@ -460,7 +466,15 @@ def gen_gf_spectra( #================================================================================= #///////////////////////////////////////////////////////////////////////////////// - + # 震源和场点层的物性，写入sac头段变量 + rcv_va = self.modarr[self.ircv, 1] + rcv_vb = self.modarr[self.ircv, 2] + rcv_rho = self.modarr[self.ircv, 3] + rcv_qainv = 1.0/self.modarr[self.ircv, 4] + rcv_qbinv = 1.0/self.modarr[self.ircv, 5] + src_va = self.modarr[self.isrc, 1] + src_vb = self.modarr[self.isrc, 2] + src_rho = self.modarr[self.isrc, 3] # 对应实际采集的地震信号，取向上为正(和理论推导使用的方向相反) dataLst = [] @@ -496,25 +510,36 @@ def gen_gf_spectra( if calc_upar: if i<2: if calc_EXP: - stream.append(EXPgrn_uiz[ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(EXPgrn_uiz[ir][i].freq2time(delayT, travtP, travtS, sgn*(-1) )) stream.append(EXPgrn_uir[ir][i].freq2time(delayT, travtP, travtS, sgn )) if calc_VF: - stream.append(VFgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(VFgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn*(-1) )) stream.append(VFgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) if calc_DC: - stream.append(DDgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(DDgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn*(-1) )) stream.append(DDgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) if calc_HF: - stream.append(HFgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(HFgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn*(-1) )) stream.append(HFgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) if calc_DC: - stream.append(DSgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(DSgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn*(-1) )) stream.append(DSgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) - stream.append(SSgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn )) + stream.append(SSgrn_uiz [ir][i].freq2time(delayT, travtP, travtS, sgn*(-1) )) stream.append(SSgrn_uir [ir][i].freq2time(delayT, travtP, travtS, sgn )) + # 在sac头段变量部分 + for tr in stream: + SAC = tr.stats.sac + SAC['user1'] = rcv_va + SAC['user2'] = rcv_vb + SAC['user3'] = rcv_rho + SAC['user4'] = rcv_qainv + SAC['user5'] = rcv_qbinv + SAC['user6'] = src_va + SAC['user7'] = src_vb + SAC['user8'] = src_rho dataLst.append(stream) @@ -522,3 +547,164 @@ def gen_gf_spectra( return dataLst + + def compute_static_grn( + self, + xarr:Union[np.ndarray,List[float],float], + yarr:Union[np.ndarray,List[float],float], + vmin_ref:float=0.0, + keps:float=-1.0, + k0:float=5.0, + Length:float=15.0, + calc_upar:bool=False, + statsfile:Union[str,None]=None): + + r""" + 调用C库计算静态格林函数，以字典的形式返回 + + :param xarr: 北向坐标数组，或单个浮点数 + :param yarr: 东向坐标数组，或单个浮点数 + :param vmin_ref: 最小参考速度（具体数值不使用），小于0则在达到积分上限后使用峰谷平均法 + （默认当震源和场点深度差<=0.5km时自动使用峰谷平均法） + :param keps: 波数k积分收敛条件，见 :ref:`(Yao and Harkrider, 1983) ` :ref:`(初稿) `， + 为负数代表不提前判断收敛，按照波数积分上限进行积分 + :param k0: 波数k积分的上限 :math:`\tilde{k_{max}}=(k_{0}*\pi/hs)^2` , 波数k积分循环必须退出, hs=max(震源和台站深度差,1.0) + :param Length: 定义波数k积分的间隔 `dk=2\pi / (L*rmax)`, 默认15；负数表示使用Filon积分 + :param calc_upar: 是否计算位移u的空间导数 + :param statsfile: 波数k积分（包括Filon积分和峰谷平均法）的过程记录文件，常用于debug或者观察积分过程中 :math:`F(k,\omega)` 和 :math:`F(k,\omega)J_m(kr)k` 的变化 + + :return: + - **dataDct* - 字典形式的格林函数 + """ + + depsrc = self.depsrc + deprcv = self.deprcv + + if isinstance(xarr, list): + xarr = np.array(xarr) + elif isinstance(xarr, float) or isinstance(xarr, int): + xarr = np.array([xarr*1.0]) + + if isinstance(yarr, list): + yarr = np.array(yarr) + elif isinstance(yarr, float) or isinstance(yarr, int): + yarr = np.array([yarr*1.0]) + + nx = len(xarr) + ny = len(yarr) + nr = nx*ny + rs = np.zeros((nr,), dtype=NPCT_REAL_TYPE) + for iy in range(ny): + for ix in range(nx): + rs[ix + iy*nx] = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) + c_rs = npct.as_ctypes(rs) + + # 参考最小速度 + if vmin_ref == 0.0: + vmin_ref = max(self.vmin, 0.1) + if abs(depsrc - deprcv) <= 1.0: + vmin_ref = - abs(vmin_ref) # 自动使用PTAM + + # 设置波数积分间隔 + if Length == 0.0: + Length = 15.0 + + # 积分状态文件 + c_statsfile = None + if statsfile is not None: + c_statsfile = c_char_p(statsfile.encode('utf-8')) + + # 初始化格林函数 + EXPgrn = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_EXPgrn = npct.as_ctypes(EXPgrn.reshape(-1)) + VFgrn = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_VFgrn = npct.as_ctypes(VFgrn.reshape(-1)) + HFgrn = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_HFgrn = npct.as_ctypes(HFgrn.reshape(-1)) + DDgrn = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_DDgrn = npct.as_ctypes(DDgrn.reshape(-1)) + DSgrn = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_DSgrn = npct.as_ctypes(DSgrn.reshape(-1)) + SSgrn = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_SSgrn = npct.as_ctypes(SSgrn.reshape(-1)) + + # 位移u的空间导数 + EXPgrn_uiz = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_EXPgrn_uiz = npct.as_ctypes(EXPgrn_uiz.reshape(-1)) + VFgrn_uiz = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_VFgrn_uiz = npct.as_ctypes(VFgrn_uiz.reshape(-1)) + HFgrn_uiz = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_HFgrn_uiz = npct.as_ctypes(HFgrn_uiz.reshape(-1)) + DDgrn_uiz = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_DDgrn_uiz = npct.as_ctypes(DDgrn_uiz.reshape(-1)) + DSgrn_uiz = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_DSgrn_uiz = npct.as_ctypes(DSgrn_uiz.reshape(-1)) + SSgrn_uiz = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_SSgrn_uiz = npct.as_ctypes(SSgrn_uiz.reshape(-1)) + + EXPgrn_uir = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_EXPgrn_uir = npct.as_ctypes(EXPgrn_uir.reshape(-1)) + VFgrn_uir = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_VFgrn_uir = npct.as_ctypes(VFgrn_uir.reshape(-1)) + HFgrn_uir = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_HFgrn_uir = npct.as_ctypes(HFgrn_uir.reshape(-1)) + DDgrn_uir = np.zeros((nr,2), dtype=NPCT_REAL_TYPE); C_DDgrn_uir = npct.as_ctypes(DDgrn_uir.reshape(-1)) + DSgrn_uir = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_DSgrn_uir = npct.as_ctypes(DSgrn_uir.reshape(-1)) + SSgrn_uir = np.zeros((nr,3), dtype=NPCT_REAL_TYPE); C_SSgrn_uir = npct.as_ctypes(SSgrn_uir.reshape(-1)) + + if not calc_upar: + C_EXPgrn_uiz = C_VFgrn_uiz = C_HFgrn_uiz = C_DDgrn_uiz = C_DSgrn_uiz = C_SSgrn_uiz = None + C_EXPgrn_uir = C_VFgrn_uir = C_HFgrn_uir = C_DDgrn_uir = C_DSgrn_uir = C_SSgrn_uir = None + + + # 运行C库函数 + #///////////////////////////////////////////////////////////////////////////////// + # 计算得到的格林函数的单位： + # 单力源 HF[ZRT],VF[ZR] 1e-15 cm/dyne + # 爆炸源 EX[ZR] 1e-20 cm/(dyne*cm) + # 剪切源 DD[ZR],DS[ZRT],SS[ZRT] 1e-20 cm/(dyne*cm) + #================================================================================= + C_integ_static_grn( + self.c_pymod1d, nr, c_rs, vmin_ref, keps, k0, Length, + C_EXPgrn, C_VFgrn, C_HFgrn, C_DDgrn, C_DSgrn, C_SSgrn, + calc_upar, + C_EXPgrn_uiz, C_VFgrn_uiz, C_HFgrn_uiz, C_DDgrn_uiz, C_DSgrn_uiz, C_SSgrn_uiz, + C_EXPgrn_uir, C_VFgrn_uir, C_HFgrn_uir, C_DDgrn_uir, C_DSgrn_uir, C_SSgrn_uir, + c_statsfile + ) + #================================================================================= + #///////////////////////////////////////////////////////////////////////////////// + + # 震源和场点层的物性 + rcv_va = self.modarr[self.ircv, 1] + rcv_vb = self.modarr[self.ircv, 2] + rcv_rho = self.modarr[self.ircv, 3] + src_va = self.modarr[self.isrc, 1] + src_vb = self.modarr[self.isrc, 2] + src_rho = self.modarr[self.isrc, 3] + + # 结果字典 + dataDct = {} + dataDct['_xarr'] = xarr.copy() + dataDct['_yarr'] = yarr.copy() + dataDct['_src_va'] = src_va + dataDct['_src_vb'] = src_vb + dataDct['_src_rho'] = src_rho + dataDct['_rcv_va'] = rcv_va + dataDct['_rcv_vb'] = rcv_vb + dataDct['_rcv_rho'] = rcv_rho + + # 整理结果，将每个格林函数以2d矩阵的形式存储，shape=(nx, ny) + for i, ch in enumerate(['Z', 'R', 'T']): + sgn = -1 if ch=='Z' else 1 + if i<2: + dataDct[f'EX{ch}'] = sgn * EXPgrn[:,i].reshape((nx, ny), order='F') + dataDct[f'VF{ch}'] = sgn * VFgrn[:,i].reshape((nx, ny), order='F') + dataDct[f'DD{ch}'] = sgn * DDgrn[:,i].reshape((nx, ny), order='F') + + dataDct[f'HF{ch}'] = sgn * HFgrn[:,i].reshape((nx, ny), order='F') + dataDct[f'DS{ch}'] = sgn * DSgrn[:,i].reshape((nx, ny), order='F') + dataDct[f'SS{ch}'] = sgn * SSgrn[:,i].reshape((nx, ny), order='F') + + if calc_upar: + if i<2: + dataDct[f'zEX{ch}'] = sgn * EXPgrn_uiz[:,i].reshape((nx, ny), order='F') * (-1) + dataDct[f'rEX{ch}'] = sgn * EXPgrn_uir[:,i].reshape((nx, ny), order='F') + dataDct[f'zVF{ch}'] = sgn * VFgrn_uiz[:,i].reshape((nx, ny), order='F') * (-1) + dataDct[f'rVF{ch}'] = sgn * VFgrn_uir[:,i].reshape((nx, ny), order='F') + dataDct[f'zDD{ch}'] = sgn * DDgrn_uiz[:,i].reshape((nx, ny), order='F') * (-1) + dataDct[f'rDD{ch}'] = sgn * DDgrn_uir[:,i].reshape((nx, ny), order='F') + + dataDct[f'zHF{ch}'] = sgn * HFgrn_uiz[:,i].reshape((nx, ny), order='F') * (-1) + dataDct[f'rHF{ch}'] = sgn * HFgrn_uir[:,i].reshape((nx, ny), order='F') + dataDct[f'zDS{ch}'] = sgn * DSgrn_uiz[:,i].reshape((nx, ny), order='F') * (-1) + dataDct[f'rDS{ch}'] = sgn * DSgrn_uir[:,i].reshape((nx, ny), order='F') + dataDct[f'zSS{ch}'] = sgn * SSgrn_uiz[:,i].reshape((nx, ny), order='F') * (-1) + dataDct[f'rSS{ch}'] = sgn * SSgrn_uir[:,i].reshape((nx, ny), order='F') + + return dataDct \ No newline at end of file diff --git a/pygrt/signals.py b/pygrt/signals.py index ed4975e9..16abef9d 100755 --- a/pygrt/signals.py +++ b/pygrt/signals.py @@ -11,6 +11,7 @@ import numpy as np import numpy.ctypeslib as npct +from ctypes import byref, cast from .c_interfaces import * diff --git a/pygrt/utils.py b/pygrt/utils.py index e09ce027..e5adbca3 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -23,10 +23,12 @@ from scipy.fft import rfft, irfft import math import os +from typing import List, Union +from copy import deepcopy from numpy.typing import ArrayLike -from .c_structures import * +from .c_interfaces import * __all__ = [ @@ -35,6 +37,9 @@ "gen_syn_from_gf_EXP", "gen_syn_from_gf_MT", + "compute_strain", + "compute_stress", + "stream_convolve", "stream_integral", "stream_diff", @@ -47,20 +52,29 @@ ] -def gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:str, M0:float, az:float, **kwargs): +#================================================================================================================= +# +# 根据辐射因子合成地震图 +# +#================================================================================================================= + +def _gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:str, M0:float, az:float, ZNE=False, **kwargs): r""" 一个发起函数，根据不同震源参数，从格林函数中合成理论地震图 :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型 :param calc_upar: 是否计算位移u的空间导数 - :param M0: 标量地震矩, 单位dyne*cm :param compute_type: 计算类型，应为以下之一: - 'COMPUTE_EXP'(爆炸源), 'COMPUTE_SF'(单力源), - 'COMPUTE_DC'(双力偶源), 'COMPUTE_MT'(矩张量源) + 'COMPUTE_EXP'(爆炸源), 'COMPUTE_SF'(单力源), + 'COMPUTE_DC'(双力偶源), 'COMPUTE_MT'(矩张量源) + :param M0: 标量地震矩, 单位dyne*cm + :param az: 方位角(度) + :param ZNE: 是否以ZNE分量输出? + """ chs = ['Z', 'R', 'T'] sacin_prefixes = ["", "z", "r", ""] # 输入通道名 - sacout_prefixes = ["", "z", "r", "t"] # 输入通道名 + sacout_prefixes = ["", "z", "r", "t"] # 输出通道名 srcName = ["EX", "VF", "HF", "DD", "DS", "SS"] allchs = [tr.stats.channel for tr in st] @@ -70,6 +84,8 @@ def gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:str, M0:float, az:fl if baz > 360: baz -= 360 + azrad = np.deg2rad(az) + calcUTypes = 4 if calc_upar else 1 stall = Stream() @@ -77,12 +93,12 @@ def gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:str, M0:float, az:fl dist = st[0].stats.sac['dist'] upar_scale:float = 1.0 for ityp in range(calcUTypes): - if ityp > 1: + if ityp > 0: upar_scale = 1e-5 if ityp == 3: upar_scale /= dist - srcCoef = set_source_coef(ityp==3, upar_scale, compute_type, M0, az, **kwargs) + srcCoef = _set_source_coef(ityp==3, upar_scale, compute_type, M0, azrad, **kwargs) inpref = sacin_prefixes[ityp] outpref = sacout_prefixes[ityp] @@ -108,12 +124,200 @@ def gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:str, M0:float, az:fl tr.data += coef*tr0.data stall.append(tr) + + if ZNE: + stall = _data_zrt2zne(stall) return stall -def set_source_coef( - par_theta:bool, coef:float, compute_type:str, M0:float, az:float, +def _gen_syn_from_static_gf(grn:dict, calc_upar:bool, compute_type:str, M0:float, ZNE=False, **kwargs): + r""" + 一个发起函数，根据不同震源参数，从静态格林函数中合成理论静态场 + + :param grn: 计算好的静态格林函数, 字典类型 + :param calc_upar: 是否计算位移u的空间导数 + :param compute_type: 计算类型，应为以下之一: + 'COMPUTE_EXP'(爆炸源), 'COMPUTE_SF'(单力源), + 'COMPUTE_DC'(双力偶源), 'COMPUTE_MT'(矩张量源) + :param M0: 标量地震矩, 单位dyne*cm + :param ZNE: 是否以ZNE分量输出? + + """ + chs = ['Z', 'R', 'T'] + sacin_prefixes = ["", "z", "r", ""] # 输入通道名 + srcName = ["EX", "VF", "HF", "DD", "DS", "SS"] + allchs = list(grn.keys()) + + s_compute_type = compute_type.split("_")[1][:2] + + calcUTypes = 4 if calc_upar else 1 + + xarr:np.ndarray = grn['_xarr'] + yarr:np.ndarray = grn['_yarr'] + + # 结果字典 + resDct = {} + + # 基本数据拷贝 + for k in grn.keys(): + if k[0] != '_': + continue + resDct[k] = deepcopy(grn[k]) + + XX = np.zeros((calcUTypes, 3, len(xarr), len(yarr)), dtype='f8') + dblsyn = (c_double*3)() + dblupar = (c_double*9)() + + for iy in range(len(yarr)): + for ix in range(len(xarr)): + # 震中距 + dist = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) + + # 方位角 + azrad = np.arctan2(yarr[iy], xarr[ix]) + + upar_scale:float = 1.0 + for ityp in range(calcUTypes): + if ityp > 0: + upar_scale = 1e-5 + if ityp == 3: + upar_scale /= dist + + srcCoef = _set_source_coef(ityp==3, upar_scale, compute_type, M0, azrad, **kwargs) + + inpref = sacin_prefixes[ityp] + + for c in range(3): + ch = chs[c] + for k in range(6): + coef = srcCoef[c, k] + if coef==0.0: + continue + + # 读入数据 + channel = f'{inpref}{srcName[k]}{ch}' + if channel not in allchs: + raise ValueError(f"Failed, channel=\"{channel}\" not exists.") + + XX[ityp, c, ix, iy] += coef*grn[channel][ix, iy] + + if ZNE: + for i in range(3): + dblsyn[i] = XX[0, i, ix, iy] + if calc_upar: + for k in range(3): + dblupar[k + i*3] = XX[i+1, k, ix, iy] + if calc_upar: + C_rot_zrt2zxy_upar(azrad, dblsyn, dblupar, dist*1e5) + else: + C_rot_zxy2zrt_vec(-azrad, dblsyn) + + for i in range(3): + XX[0, i, ix, iy] = dblsyn[i] + if calc_upar: + for k in range(3): + XX[i+1, k, ix, iy] = dblupar[k + i*3] + + + # 将XX数组分到字典中 + if ZNE: + chs = ['Z', 'N', 'E'] + + for ityp in range(calcUTypes): + c1 = '' if ityp==0 else chs[ityp-1].lower() + for c in range(3): + resDct[f"{c1}{s_compute_type}{chs[c]}"] = XX[ityp, c].copy() + + return resDct + + +def _data_zrt2zne(stall:Stream): + r""" + 将位移分量和位移空间导数分量转为ZNE坐标系 + + :param stall: 柱坐标系(zrt)下合成地震图 + + :return: + - **stream** - :class:`obspy.Stream` 类型 + """ + + chs = ['Z', 'R', 'T'] + + midname = stall[0].stats.channel[-3:-1] + + synLst:List[Trace] = [] # 顺序要求Z, R, T + uparLst:List[Trace] = [] # 顺序要求zXXZ, zXXR, zXXT, rXXZ, rXXR, rXXT, tXXZ, tXXR, tXXT + stsyn_upar = Stream() + for ch in chs: + st = stall.select(channel=f"{midname}{ch}") + if len(st) == 1: + synLst.append(st[0]) + + for ch2 in chs: + st = stall.select(channel=f"{ch.lower()}{midname}{ch2}") + if len(st) == 1: + uparLst.append(st[0]) + + if len(synLst) != 3: + raise ValueError(f"WRONG! synLst should have 3 components.") + if len(stsyn_upar) != 0 and len(stsyn_upar) != 9: + raise ValueError(f"WRONG! stsyn_upar should have 0 or 9 components.") + + + # 是否有空间导数 + doupar = (len(uparLst) == 9) + + nt = stall[0].stats.npts + azrad = np.deg2rad(stall[0].stats.sac['az']) + dist = stall[0].stats.sac['dist'] + + dblsyn = (c_double * 3)() + dbleupar = (c_double * 9)() + + # 对每一个时间点 + for n in range(nt): + # 复制数据 + for i1 in range(3): + dblsyn[i1] = synLst[i1].data[n] + if doupar: + for i2 in range(3): + dbleupar[i2 + i1*3] = uparLst[i2 + i1*3].data[n] + + if doupar: + C_rot_zrt2zxy_upar(azrad, dblsyn, dbleupar, dist*1e5) + else: + C_rot_zxy2zrt_vec(-azrad, dblsyn) + + # 将结果写入原数组 + for i1 in range(3): + synLst[i1].data[n] = dblsyn[i1] + if doupar: + for i2 in range(3): + uparLst[i2 + i1*3].data[n] = dbleupar[i2 + i1*3] + + znechs = ['Z', 'N', 'E'] + # 修改通道名 + for i1 in range(3): + ch1 = znechs[i1] + tr = synLst[i1] + tr.stats.channel = tr.stats.sac['kcmpnm'] = f'{midname}{ch1}' + if doupar: + for i2 in range(3): + ch2 = znechs[i2] + tr = uparLst[i2 + i1*3] + tr.stats.channel = tr.stats.sac['kcmpnm'] = f'{ch1.lower()}{midname}{ch2}' + + stres = Stream() + stres.extend(synLst) + if doupar: + stres.extend(uparLst) + + return stres + + +def _set_source_coef( + par_theta:bool, coef:float, compute_type:str, M0:float, azrad:float, fZ=None, fN=None, fE=None, strike=None, dip=None, rake=None, MT=None, **kwargs): @@ -123,10 +327,10 @@ def set_source_coef( :param par_theta: 是否求对theta的偏导 :param coef: 比例系数 :param compute_type: 计算类型，应为以下之一: - 'COMPUTE_EXP'(爆炸源), 'COMPUTE_SF'(单力源), - 'COMPUTE_DC'(双力偶源), 'COMPUTE_MT'(矩张量源) + 'COMPUTE_EXP'(爆炸源), 'COMPUTE_SF'(单力源), + 'COMPUTE_DC'(双力偶源), 'COMPUTE_MT'(矩张量源) :param M0: 地震矩 - :param az: 方位角(度) + :param azrad: 方位角(弧度) - 其他参数根据计算类型可选: - 单力源需要: fZ, fN, fE, @@ -134,10 +338,6 @@ def set_source_coef( - 矩张量源需要: MT=(Mxx, Mxy, Mxz, Myy, Myz, Mzz) """ - # 定义常数: 1度对应的弧度值 - DEG1 = np.pi / 180.0 - - azrad = az*DEG1 caz = np.cos(azrad) saz = np.sin(azrad) @@ -171,9 +371,9 @@ def set_source_coef( elif compute_type == 'COMPUTE_DC': # 双力偶源情况 (公式4.8.35) # 计算各种角度值(转为弧度) - stkrad = strike * DEG1 # 走向角 - diprad = dip * DEG1 # 倾角 - rakrad = rake * DEG1 # 滑动角 + stkrad = np.deg2rad(strike) # 走向角 + diprad = np.deg2rad(dip) # 倾角 + rakrad = np.deg2rad(rake) # 滑动角 therad = azrad - stkrad # 方位角与走向角差 # 计算各种三角函数值 @@ -234,73 +434,546 @@ def set_source_coef( return src_coef -def gen_syn_from_gf_DC(st:Stream, M0:float, strike:float, dip:float, rake:float, az:float, calc_upar:bool=False): +def gen_syn_from_gf_DC(st:Union[Stream,dict], M0:float, strike:float, dip:float, rake:float, az:float=-999, ZNE=False, calc_upar:bool=False): ''' 双力偶源，角度单位均为度 - :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型 + :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型，或者静态格林函数（字典类型） :param M0: 标量地震矩, 单位dyne*cm :param strike: 走向，以北顺时针为正，0<=strike<=360 :param dip: 倾角，以水平面往下为正，0<=dip<=90 :param rake: 滑动角，在断层面相对于走向方向逆时针为正，-180<=rake<=180 - :param az: 台站方位角，以北顺时针为正，0<=az<=360 + :param az: 台站方位角，以北顺时针为正，0<=az<=360（静态情况不需要） + :param ZNE: 是否以ZNE分量输出? :param calc_upar: 是否计算位移u的空间导数 :return: - - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 + - **stream** - 三分量地震图, :class:`obspy.Stream` 类型 ''' - - return gen_syn_from_gf(st, calc_upar, "COMPUTE_DC", M0, az, strike=strike, dip=dip, rake=rake) + if isinstance(st, Stream): + if az > 360 or az < -360: + raise ValueError(f"WRONG azimuth ({az})") + return _gen_syn_from_gf(st, calc_upar, "COMPUTE_DC", M0, az, ZNE, strike=strike, dip=dip, rake=rake) + elif isinstance(st, dict): + return _gen_syn_from_static_gf(st, calc_upar, "COMPUTE_DC", M0, ZNE, strike=strike, dip=dip, rake=rake) + else: + raise NotImplementedError -def gen_syn_from_gf_SF(st:Stream, S:float, fN:float, fE:float, fZ:float, az:float, calc_upar:bool=False): +def gen_syn_from_gf_SF(st:Union[Stream,dict], S:float, fN:float, fE:float, fZ:float, az:float=-999, ZNE=False, calc_upar:bool=False): ''' 单力源，力分量单位均为dyne - :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型 + :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型，或者静态格林函数（字典类型） :param S: 力的放大系数 :param fN: 北向力，向北为正 :param fE: 东向力，向东为正 :param fZ: 垂向力，向下为正 - :param az: 台站方位角，以北顺时针为正，0<=az<=360 + :param az: 台站方位角，以北顺时针为正，0<=az<=360 （静态情况不需要） + :param ZNE: 是否以ZNE分量输出? :param calc_upar: 是否计算位移u的空间导数 :return: - - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 + - **stream** - 三分量地震图, :class:`obspy.Stream` 类型 ''' - return gen_syn_from_gf(st, calc_upar, "COMPUTE_SF", S, az, fN=fN, fE=fE, fZ=fZ) + if isinstance(st, Stream): + if az > 360 or az < -360: + raise ValueError(f"WRONG azimuth ({az})") + return _gen_syn_from_gf(st, calc_upar, "COMPUTE_SF", S, az, ZNE, fN=fN, fE=fE, fZ=fZ) + elif isinstance(st, dict): + return _gen_syn_from_static_gf(st, calc_upar, "COMPUTE_SF", S, ZNE, fN=fN, fE=fE, fZ=fZ) + else: + raise NotImplementedError -def gen_syn_from_gf_EXP(st:Stream, M0:float, az:float, calc_upar:bool=False): +def gen_syn_from_gf_EXP(st:Union[Stream,dict], M0:float, az:float=-999, ZNE=False, calc_upar:bool=False): ''' 爆炸源 - :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型 + :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型，或者静态格林函数（字典类型） :param M0: 标量地震矩, 单位dyne*cm - :param az: 台站方位角，以北顺时针为正，0<=az<=360 [不用于计算] + :param az: 台站方位角，以北顺时针为正，0<=az<=360 [不用于计算] （静态情况不需要） + :param ZNE: 是否以ZNE分量输出? :param calc_upar: 是否计算位移u的空间导数 :return: - - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 + - **stream** - 三分量地震图, :class:`obspy.Stream` 类型 ''' - return gen_syn_from_gf(st, calc_upar, "COMPUTE_EXP", M0, az) - + if isinstance(st, Stream): + if az > 360 or az < -360: + raise ValueError(f"WRONG azimuth ({az})") + return _gen_syn_from_gf(st, calc_upar, "COMPUTE_EXP", M0, az, ZNE) + elif isinstance(st, dict): + return _gen_syn_from_static_gf(st, calc_upar, "COMPUTE_EXP", M0, ZNE) + else: + raise NotImplementedError + -def gen_syn_from_gf_MT(st:Stream, M0:float, MT:ArrayLike, az:float, calc_upar:bool=False): +def gen_syn_from_gf_MT(st:Union[Stream,dict], M0:float, MT:ArrayLike, az:float=-999, ZNE=False, calc_upar:bool=False): ''' 矩张量源，单位为dyne*cm - :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型 + :param st: 计算好的时域格林函数, :class:`obspy.Stream` 类型，或者静态格林函数（字典类型） :param M0: 标量地震矩 :param MT: 矩张量 (M11, M12, M13, M22, M23, M33),下标1,2,3分别代表北向，东向，垂直向上 - :param az: 台站方位角，以北顺时针为正，0<=az<=360 + :param az: 台站方位角，以北顺时针为正，0<=az<=360 （静态情况不需要） + :param ZNE: 是否以ZNE分量输出? :param calc_upar: 是否计算位移u的空间导数 :return: - - **stream** - 三分量ZRT地震图, :class:`obspy.Stream` 类型 + - **stream** - 三分量地震图, :class:`obspy.Stream` 类型 ''' - return gen_syn_from_gf(st, calc_upar, "COMPUTE_MT", M0, az, MT=MT) + if isinstance(st, Stream): + if az > 360 or az < -360: + raise ValueError(f"WRONG azimuth ({az})") + return _gen_syn_from_gf(st, calc_upar, "COMPUTE_MT", M0, az, ZNE, MT=MT) + elif isinstance(st, dict): + return _gen_syn_from_static_gf(st, calc_upar, "COMPUTE_MT", M0, ZNE, MT=MT) + else: + raise NotImplementedError + + +#================================================================================================================= +# +# 根据几何方程和本构方程合成应力和应变 +# +#================================================================================================================= + + +def _compute_strain(st_syn:Stream): + r""" + Compute dynamic strain from synthetic spatial derivatives. + + :param st_syn: synthetic spatial derivatives. + + :return: + - **stream** - dynamic strain, in :class:`obspy.Stream` class. + """ + + midname = st_syn[0].stats.channel[-3:-1] + + # 检查是否每个trace是否具有相同midname + for tr in st_syn: + if tr.stats.channel[-3:-1] != midname: + raise ValueError("WRONG INPUT! inconsistent component names.") + + zrtchs = ['Z', 'R', 'T'] + znechs = ['Z', 'N', 'E'] + chs = zrtchs + + # 判断是否有标志性的trace + if len(st_syn.select(channel=f"n{midname}N")) > 0: + chs = znechs + + dist = st_syn[0].stats.sac['dist'] + + # ---------------------------------------------------------------------------------- + # 循环6个分量 + stres = Stream() + for i1 in range(3): + c1 = chs[i1] + for i2 in range(i1, 3): + c2 = chs[i2] + + channel = f"{c1.lower()}{midname}{c2}" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + tr = st[0].copy() + + channel = f"{c2.lower()}{midname}{c1}" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + tr.data = (tr.data + st[0].data) * 0.5 + + # 特殊情况加上协变导数 + if c1=='R' and c2=='T': + channel = f"{midname}T" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + tr.data -= 0.5*st[0].data / dist * 1e-5 + + elif c1=='T' and c2=='T': + channel = f"{midname}R" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + tr.data += st[0].data / dist * 1e-5 + + # 修改通道名 + tr.stats.channel = tr.stats.sac['kcmpnm'] = f"{c1}{c2}" + + stres.append(tr) + + return stres + + +def _compute_static_strain(syn:dict): + r""" + Compute static strain from synthetic spatial derivatives. + + :param syn: synthetic spatial derivatives. + + :return: + - **res** - static strain, in dict class. + """ + + midname = "" + # 检查是否每个分量是否具有相同midname + for k in syn.keys(): + if k[0] == '_': + continue + if len(midname)==0: + midname = k[-3:-1] + + if k[-3:-1] != midname: + raise ValueError("WRONG INPUT! inconsistent component names.") + + zrtchs = ['Z', 'R', 'T'] + znechs = ['Z', 'N', 'E'] + chs = zrtchs + + # 判断是否有标志性的分量名 + if f"n{midname}N" in syn.keys(): + chs = znechs + + xarr:np.ndarray = syn['_xarr'] + yarr:np.ndarray = syn['_yarr'] + + # 结果字典 + resDct = {} + + # 基本数据拷贝 + for k in syn.keys(): + if k[0] != '_': + continue + resDct[k] = deepcopy(syn[k]) + + # 6个分量建立数组 + for i1 in range(3): + c1 = chs[i1] + for i2 in range(i1, 3): + c2 = chs[i2] + channel = f"{c1}{c2}" + resDct[channel] = np.zeros((len(xarr), len(yarr)), dtype='f8') + + + for iy in range(len(yarr)): + for ix in range(len(xarr)): + # 震中距 + dist = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) + + # ---------------------------------------------------------------------------------- + # 循环6个分量 + for i1 in range(3): + c1 = chs[i1] + for i2 in range(i1, 3): + c2 = chs[i2] + + channel = f"{c1.lower()}{midname}{c2}" + v12 = syn[channel][ix, iy] + + channel = f"{c2.lower()}{midname}{c1}" + v21 = syn[channel][ix, iy] + + val = 0.5*(v12 + v21) + + # 特殊情况加上协变导数 + if c1=='R' and c2=='T': + channel = f"{midname}T" + v0 = syn[channel][ix, iy] + val -= 0.5*v0 / dist * 1e-5 + + elif c1=='T' and c2=='T': + channel = f"{midname}R" + v0 = syn[channel][ix, iy] + val += v0 / dist * 1e-5 + + channel = f"{c1}{c2}" + resDct[channel][ix, iy] = val + + return resDct + + +def compute_strain(st:Union[Stream,dict]): + r""" + Compute dynamic/static strain from synthetic spatial derivatives. + + :param st: synthetic spatial derivatives + :class:`obspy.Stream` class for dynamic case, dict class for static case. + + :return: + - **stres** - dynamic/static strain, in :class:`obspy.Stream` class or dict class. + """ + if isinstance(st, Stream): + return _compute_strain(st) + elif isinstance(st, dict): + return _compute_static_strain(st) + else: + raise NotImplementedError + + +def _compute_stress(st_syn:Stream): + r""" + Compute dynamic stress from synthetic spatial derivatives. + + :param st_syn: synthetic spatial derivatives. + + :return: + - **stream** - dynamic stress (unit: dyne/cm^2 = 0.1 Pa), in :class:`obspy.Stream` class. + """ + + # 由于有Q值的存在，lambda和mu变成了复数，需在频域进行 + + midname = st_syn[0].stats.channel[-3:-1] + + # 检查是否每个trace是否具有相同midname + for tr in st_syn: + if tr.stats.channel[-3:-1] != midname: + raise ValueError("WRONG INPUT! inconsistent component names.") + + zrtchs = ['Z', 'R', 'T'] + znechs = ['Z', 'N', 'E'] + chs = zrtchs + rot2ZNE:bool = False + + # 判断是否有标志性的trace + if len(st_syn.select(channel=f"n{midname}N")) > 0: + chs = znechs + rot2ZNE = True + + nt = st_syn[0].stats.npts + dt = st_syn[0].stats.delta + dist = st_syn[0].stats.sac['dist'] + df = 1.0/(nt*dt) + nf = nt//2 + 1 + va = st_syn[0].stats.sac['user1'] + vb = st_syn[0].stats.sac['user2'] + rho = st_syn[0].stats.sac['user3'] + Qainv = st_syn[0].stats.sac['user4'] + Qbinv = st_syn[0].stats.sac['user5'] + + # 计算不同频率下的拉梅系数 + mus = np.zeros((nf,), dtype='c16') + lams = np.zeros((nf,), dtype='c16') + omega = (REAL*2)(0.0, 0.0) + atte = (REAL*2)(0.0, 0.0) + for i in range(nf): + freq = 0.01 if i==0 else df*i + w = 2.0*np.pi*freq + omega[0] = w + C_py_attenuation_law(Qbinv, omega, atte) + attb = atte[0] + atte[1]*1j + mus[i] = vb*vb*attb*attb*rho*1e10 + C_py_attenuation_law(Qainv, omega, atte) + atta = atte[0] + atte[1]*1j + lams[i] = va*va*atta*atta*rho*1e10 - 2.0*mus[i] + del omega, atte + + # ---------------------------------------------------------------------------------- + # 先计算体积应变u_kk = u_11 + u22 + u33 和 lamda的乘积 + lam_ukk = np.zeros((nf,), dtype='c16') + for i in range(3): + c = chs[i] + channel = f"{c.lower()}{midname}{c}" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + lam_ukk[:] += rfft(st[0].data, nt) + + # 加上协变导数 + if not rot2ZNE: + channel = f"{midname}R" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + lam_ukk[:] += rfft(st[0].data, nt) / dist * 1e-5 + + lam_ukk[:] *= lams + + # ---------------------------------------------------------------------------------- + # 循环6个分量 + stres = Stream() + for i1 in range(3): + c1 = chs[i1] + for i2 in range(i1, 3): + c2 = chs[i2] + + channel = f"{c1.lower()}{midname}{c2}" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + tr = st[0].copy() + fftarr = np.zeros((nf,), dtype='c16') + + channel = f"{c2.lower()}{midname}{c1}" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + fftarr[:] = rfft(tr.data + st[0].data, nt) * mus + + # 对于对角线分量，需加上lambda * u_kk + if c1==c2: + fftarr[:] += lam_ukk + + # 特殊情况加上协变导数 + if c1=='R' and c2=='T': + channel = f"{midname}T" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + fftarr[:] -= mus*rfft(st[0].data, nt) / dist * 1e-5 + + elif c1=='T' and c2=='T': + channel = f"{midname}R" + st = st_syn.select(channel=channel) + if len(st) == 0: + raise NameError(f"{channel} not exists.") + fftarr[:] += 2.0*mus*rfft(st[0].data, nt) / dist * 1e-5 + + # 修改通道名 + tr.stats.channel = tr.stats.sac['kcmpnm'] = f"{c1}{c2}" + + tr.data = irfft(fftarr, nt) + + stres.append(tr) + + return stres + + +def _compute_static_stress(syn:dict): + r""" + Compute static stress from synthetic spatial derivatives. + + :param syn: synthetic spatial derivatives. + + :return: + - **res** - static stress (unit: dyne/cm^2 = 0.1 Pa), in dict class. + """ + + midname = "" + # 检查是否每个分量是否具有相同midname + for k in syn.keys(): + if k[0] == '_': + continue + if len(midname)==0: + midname = k[-3:-1] + + if k[-3:-1] != midname: + raise ValueError("WRONG INPUT! inconsistent component names.") + + zrtchs = ['Z', 'R', 'T'] + znechs = ['Z', 'N', 'E'] + chs = zrtchs + rot2ZNE:bool = False + + # 判断是否有标志性的分量名 + if f"n{midname}N" in syn.keys(): + chs = znechs + rot2ZNE = True + + xarr:np.ndarray = syn['_xarr'] + yarr:np.ndarray = syn['_yarr'] + va = syn['_rcv_va'] + vb = syn['_rcv_vb'] + rho = syn['_rcv_rho'] + mu = vb*vb*rho*1e10 + lam = va*va*rho*1e10 - 2.0*mu + + # 结果字典 + resDct = {} + + # 基本数据拷贝 + for k in syn.keys(): + if k[0] != '_': + continue + resDct[k] = deepcopy(syn[k]) + + # 6个分量建立数组 + for i1 in range(3): + c1 = chs[i1] + for i2 in range(i1, 3): + c2 = chs[i2] + channel = f"{c1}{c2}" + resDct[channel] = np.zeros((len(xarr), len(yarr)), dtype='f8') + + + for iy in range(len(yarr)): + for ix in range(len(xarr)): + # 震中距 + dist = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) + + # ---------------------------------------------------------------------------------- + # 先计算体积应变u_kk = u_11 + u22 + u33 和 lamda的乘积 + lam_ukk = 0.0 + for i in range(3): + c = chs[i] + channel = f"{c.lower()}{midname}{c}" + lam_ukk += syn[channel][ix, iy] + + # 加上协变导数 + if not rot2ZNE: + channel = f"{midname}R" + lam_ukk += syn[channel][ix, iy] / dist * 1e-5 + + lam_ukk *= lam + + # ---------------------------------------------------------------------------------- + # 循环6个分量 + for i1 in range(3): + c1 = chs[i1] + for i2 in range(i1, 3): + c2 = chs[i2] + + channel = f"{c1.lower()}{midname}{c2}" + v12 = syn[channel][ix, iy] + + channel = f"{c2.lower()}{midname}{c1}" + v21 = syn[channel][ix, iy] + + val = mu*(v12 + v21) + + # 对于对角线分量，需加上lambda * u_kk + if c1==c2: + val += lam_ukk + + # 特殊情况加上协变导数 + if c1=='R' and c2=='T': + channel = f"{midname}T" + v0 = syn[channel][ix, iy] + val -= mu*v0 / dist * 1e-5 + + elif c1=='T' and c2=='T': + channel = f"{midname}R" + v0 = syn[channel][ix, iy] + val += 2.0*mu*v0 / dist * 1e-5 + + channel = f"{c1}{c2}" + resDct[channel][ix, iy] = val + + return resDct + + +def compute_stress(st:Union[Stream,dict]): + r""" + Compute dynamic/static stress from synthetic spatial derivatives. + + :param st: synthetic spatial derivatives + :class:`obspy.Stream` class for dynamic case, dict class for static case. + + :return: + - **stres** - dynamic/static stress (unit: dyne/cm^2 = 0.1 Pa), in :class:`obspy.Stream` class or dict class. + """ + if isinstance(st, Stream): + return _compute_stress(st) + elif isinstance(st, dict): + return _compute_static_stress(st) + else: + raise NotImplementedError + def __check_trace_attr_sac(tr:Trace, **kwargs): ''' @@ -313,6 +986,14 @@ def __check_trace_attr_sac(tr:Trace, **kwargs): tr.stats.sac = AttribDict(**kwargs) +#================================================================================================================= +# +# 卷积、微分、积分、保存SAC +# +#================================================================================================================= + + + # def stream_convolve(st0:Stream, timearr:np.ndarray, inplace=True): # ''' # 频域实现线性卷 @@ -479,6 +1160,11 @@ def stream_write_sac(st:Stream, path:str): +#================================================================================================================= +# +# 积分过程文件读取及绘制 +# +#================================================================================================================= def read_statsfile(statsfile:str):