diff --git a/docs/source/Tutorial/dynamic/gfunc.rst b/docs/source/Tutorial/dynamic/gfunc.rst index e7b07087..7adc37ce 100644 --- a/docs/source/Tutorial/dynamic/gfunc.rst +++ b/docs/source/Tutorial/dynamic/gfunc.rst @@ -36,10 +36,22 @@ Python中计算动态格林函数的主函数为 :func:`compute_grn() HFZ +# END grt.b2a +head -n 10 HFZ > HFZ_head +echo "..." >> HFZ_head + # BEGIN SYN EXP # 合成结果在 syn_exp/ 目录下，以SAC格式保存。 grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -Osyn_exp diff --git a/docs/source/Tutorial/dynamic/run_upar/milrow b/docs/source/Tutorial/dynamic/run_upar/milrow new file mode 100644 index 00000000..a99f4912 --- /dev/null +++ b/docs/source/Tutorial/dynamic/run_upar/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 +1. 6.9 4.0 2.8 9e10 9e10 +0.0 8.2 4.7 3.2 9e10 9e10 \ No newline at end of file diff --git a/docs/source/Tutorial/dynamic/run_upar/run.py b/docs/source/Tutorial/dynamic/run_upar/run.py new file mode 100644 index 00000000..a648f2fd --- /dev/null +++ b/docs/source/Tutorial/dynamic/run_upar/run.py @@ -0,0 +1,173 @@ +# BEGIN GRN +import numpy as np +import pygrt + +modarr = np.loadtxt("milrow") + +pymod = pygrt.PyModel1D(modarr, depsrc=2.0, deprcv=0.0) + +# 传入calc_upar=True计算空间导数 +stgrn = pymod.compute_grn(distarr=[10], nt=500, dt=0.02, calc_upar=True)[0] +print(stgrn.__str__(extended=True)) +# 45 Trace(s) in Stream: +# .SYN..EXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..VFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DDZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..HFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..SSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zEXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rEXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zVFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rVFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDDZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDDZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zHFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rHFZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zSSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rSSZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..EXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..VFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DDR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..HFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..SSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zEXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rEXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zVFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rVFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDDR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDDR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zHFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rHFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zSSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rSSR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..HFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..SST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zHFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rHFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zSST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rSST | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END GRN + + + +# BEGIN SYN DC +# 传入calc_upar=True计算空间导数 +stsyn = pygrt.utils.gen_syn_from_gf_DC(stgrn, M0=1e24, strike=33, dip=50, rake=120, az=30, calc_upar=True) +print(stsyn) +# 12 Trace(s) in Stream: +# .SYN..DCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..rDCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..tDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..tDCR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..tDCT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END SYN DC + + +# BEGIN ZNE +# 传入ZNE=True可返回ZNE分量 +stsyn = pygrt.utils.gen_syn_from_gf_DC(stgrn, M0=1e24, strike=33, dip=50, rake=120, az=30, calc_upar=True, ZNE=True) +print(stsyn) +# 12 Trace(s) in Stream: +# .SYN..DCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..DCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..zDCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..nDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..nDCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..nDCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..eDCZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..eDCN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..eDCE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END ZNE + +from obspy import * +import matplotlib.pyplot as plt + +def plot6(st6:Stream, title:str, out:str|None=None): + nt = st6[0].stats.npts + dt = st6[0].stats.delta + t = np.arange(nt)*dt + + MAX, MIN = 0, 9e30 + for tr in st6: + d = np.abs(tr.data) + if MAX < np.max(d): + MAX = np.max(d) + if MIN > np.min(d): + MIN = np.min(d) + + travtP = stsyn[0].stats.sac['t0'] + travtS = stsyn[0].stats.sac['t1'] + + fig, axs = plt.subplots(6, 1, figsize=(10, 7), gridspec_kw=dict(hspace=0.0), sharex=True) + for i in range(6): + tr = st6[i] + ax = axs[i] + + ax.plot(t, tr.data, c='k', lw=0.5, label=tr.stats.channel) + ax.legend(loc='upper left') + + ylims = ax.get_ylim() + if np.max(np.abs(np.array(ylims)))/MAX < 1e-5: + ylims = [-1, 1] + ax.set_ylim(ylims) + + # 绘制到时 + ax.vlines(travtP, *ylims, colors='b') + ax.text(travtP, ylims[1], "P", ha='left', va='top', color='b') + ax.vlines(travtS, *ylims, colors='r') + ax.text(travtS, ylims[1], "S", ha='left', va='top', color='r') + + ax.set_xlim([t[0], t[-1]]) + ax.set_ylim(np.array(ylims)*1.2) + + axs[0].set_title(title) + + if out is not None: + fig.savefig(out, dpi=100) + + + +# BEGIN STRAIN +st_strain = pygrt.utils.compute_strain(stsyn) +print(st_strain) +# 6 Trace(s) in Stream: +# .SYN..ZZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..ZN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..ZE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..NN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..NE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..EE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END STRAIN +plot6(st_strain, "Strain", "strain.png") + +# BEGIN STRESS +st_stress = pygrt.utils.compute_stress(stsyn) +print(st_stress) +# 6 Trace(s) in Stream: +# .SYN..ZZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..ZN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..ZE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..NN | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..NE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# .SYN..EE | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:09.980000Z | 50.0 Hz, 500 samples +# END STRESS +plot6(st_stress, "Stress", "stress.png") diff --git a/docs/source/Tutorial/dynamic/run_upar/run.sh b/docs/source/Tutorial/dynamic/run_upar/run.sh new file mode 100755 index 00000000..8e86afd4 --- /dev/null +++ b/docs/source/Tutorial/dynamic/run_upar/run.sh @@ -0,0 +1,28 @@ +#!/bin/bash + + +# BEGIN GRN +# -e 表示计算空间导数 +grt -Mmilrow -D2/0 -N500/0.02 -OGRN -R10 -e +# END GRN + +# BEGIN SYN DC +# -e 表示计算空间导数 +grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50/120 -Osyn_dc -e +# END SYN DC + +# BEGIN ZNE +# -N 表示输出ZNE分量 +grt.syn -GGRN/milrow_2_0_10 -S1e24 -A30 -M33/50/120 -Osyn_dc_zne -e -N +# END ZNE + + +# BEGIN STRAIN +# 指定文件夹以及中间名，会在原文件夹内输出SAC格式的应变 +grt.strain syn_dc_zne/out +# END STRAIN + +# BEGIN STRESS +# 指定文件夹以及中间名，会在原文件夹内输出SAC格式的应变 +grt.stress syn_dc_zne/out +# END STRESS \ No newline at end of file diff --git a/docs/source/Tutorial/dynamic/run_upar/strain.png b/docs/source/Tutorial/dynamic/run_upar/strain.png new file mode 100644 index 00000000..8cc34bd2 Binary files /dev/null and b/docs/source/Tutorial/dynamic/run_upar/strain.png differ diff --git a/docs/source/Tutorial/dynamic/run_upar/stress.png b/docs/source/Tutorial/dynamic/run_upar/stress.png new file mode 100644 index 00000000..fd779dd5 Binary files /dev/null and b/docs/source/Tutorial/dynamic/run_upar/stress.png differ diff --git a/docs/source/Tutorial/dynamic/strain_stress.rst b/docs/source/Tutorial/dynamic/strain_stress.rst index fc6815ec..7d5a811d 100644 --- a/docs/source/Tutorial/dynamic/strain_stress.rst +++ b/docs/source/Tutorial/dynamic/strain_stress.rst @@ -1,3 +1,5 @@ +.. _strain_stress_rst: + 计算动态应变和应力 =================== @@ -6,4 +8,183 @@ ----------------------------------------------------------- -TODO + +计算应变张量与应力张量需要计算位移的空间导数 :math:`\partial z,\partial r,\partial \theta` [#]_ 。 从前文关于格林函数的计算 （ :ref:`gfunc_rst` ） 和位移合成的公式 （ :ref:`syn_rst` ） 出发，位移空间导数分两个步骤计算，再根据几何方程和本构方程合成应变和应力。 + +.. [#] :math:`\partial z,\partial r,\partial \theta` 为 :math:`\dfrac{\partial}{\partial z},\dfrac{\partial}{\partial r},\dfrac{\partial}{\partial \theta}` 的简写。 + +计算格林函数阶段，求出 :math:`\partial z,\partial r` +----------------------------------------------------------------------------------------------- + +假设在 :file:`milrow` 模型中，震源深度2km，接收点位于地表，震中距为10km，计算格林函数及其空间导数，要求采样点数500个点，采样间隔0.02s。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run_upar/run.sh + :language: bash + :start-after: BEGIN GRN + :end-before: END GRN + + 在 :rst:dir:`GRN/milrow_{depsrc}_{deprcv}_{dist}/` 路径下，文件名开头有"r"和"z"就代表 :math:`\partial r` 和 :math:`\partial z`。 + + .. tab:: Python + + .. literalinclude:: run_upar/run.py + :language: python + :start-after: BEGIN GRN + :end-before: END GRN + + 通道名开头有"r"和"z"就代表 :math:`\partial r` 和 :math:`\partial z`。 + + +.. note:: + + 计算格林函数阶段求出的空间导数的单位： + + + 爆炸源： :math:`10^{-25} \, /\text{dyne} \cdot \text{cm}` + + 单力源： :math:`10^{-20} \, /\text{dyne}` + + 剪切源： :math:`10^{-25} \, /\text{dyne} \cdot \text{cm}` + + 矩张量源： :math:`10^{-25} \, /\text{dyne} \cdot \text{cm}` + + + +合成位移阶段，求出 :math:`\partial \theta / r` 并合成位移空间导数 +----------------------------------------------------------------------------- + +假设震源为剪切源，断层走向33°，倾角50°，滑动角120°，标量矩 1e24 dyne·cm，方位角30°。待求项中有 :math:`r` 是为了保证量纲统一。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run_upar/run.sh + :language: bash + :start-after: BEGIN SYN DC + :end-before: END SYN DC + + 在 :rst:dir:`syn_dc/` 路径下，文件名开头有"r"，"z"，"t"分别代表 :math:`\partial r`， :math:`\partial z`，:math:`\partial \theta / r`。 + + .. tab:: Python + + .. literalinclude:: run_upar/run.py + :language: python + :start-after: BEGIN SYN DC + :end-before: END SYN DC + + 通道名开头有"r"，"z"，"t"分别代表 :math:`\partial r`， :math:`\partial z`，:math:`\partial \theta / r`。 + + +以上计算的空间导数在柱坐标系下，即 :math:`\dfrac{\partial u(z,r,\theta)}{\partial(z,r,\theta)}`。如果需要 :math:`\dfrac{\partial u(z,x,y)}{\partial(z,x,y)}` （ :math:`x` 表示北向， :math:`y` 表示东向），传入适当参数即可，则输出结果开头的标识符从"z/r/t"变为"z/x/y"。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run_upar/run.sh + :language: bash + :start-after: BEGIN ZNE + :end-before: END ZNE + + 在 :rst:dir:`syn_dc_zne/` 路径下，文件名开头有"z"，"n"，"e"分别代表 :math:`\partial z`， :math:`\partial x`，:math:`\partial y`。 + + .. tab:: Python + + .. literalinclude:: run_upar/run.py + :language: python + :start-after: BEGIN ZNE + :end-before: END ZNE + + 通道名开头有"z"，"n"，"e"分别代表 :math:`\partial z`， :math:`\partial x`，:math:`\partial y`。 + + .. warning:: + + 不可以使用ObsPy库中提供的Stream.rotate()函数进行位移空间导数的旋转。 + + +.. note:: + + 合成的位移空间导数单位为1。 + + +合成应变和应力 +----------------------- +应变和应力均为二阶对称张量，故将输出6个独立分量。程序选择应变应力的坐标系的依据是——根据文件名/通道名判断是否存在 :math:`\dfrac{\partial u_x}{\partial x}`，如果有则使用ZNE坐标系，否则使用ZRT坐标系。 **所以建议保存结果的文件夹中只使用同种坐标系**，就像上面分为 :rst:dir:`syn_dc/` 和 :rst:dir:`syn_dc_zne/` 两个文件夹保存。 + +应变张量 +~~~~~~~~~~~~~~ +根据几何方程 [#]_ + +.. math:: + + e_{ji} = \dfrac{1}{2} \left( u_{i,j} + u_{j,i} \right) = \dfrac{1}{2} \left( \dfrac{\partial u_i}{\partial x_j} + \dfrac{\partial u_j}{\partial x_i} \right) + + +.. [#] 这只适用于ZNE坐标系，对于ZRT坐标系需考虑协变导数。程序中已考虑，这里不再做公式介绍。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run_upar/run.sh + :language: bash + :start-after: BEGIN STRAIN + :end-before: END STRAIN + + 在 :rst:dir:`syn_dc_zne/` 原路径下，生成 :file:`*.strain.??.sac`，文件名中包括分量名，如ZZ、ZN等。 + + .. tab:: Python + + .. literalinclude:: run_upar/run.py + :language: python + :start-after: BEGIN STRAIN + :end-before: END STRAIN + + 返回的 |Stream| 通道名即为分量名，如ZZ、ZN等。 + +.. image:: run_upar/strain.png + :align: center + + +应力张量 +~~~~~~~~~~~~~~ +根据各向同性介质的本构方程 [#]_ + +.. math:: + + \sigma_{ji} = \lambda \delta_{ij} e_{kk} + 2 \mu e_{ij} = \lambda \delta_{ij} u_{kk} + \mu \left( u_{i,j} + u_{j,i} \right) + + +.. [#] 这只适用于ZNE坐标系，对于ZRT坐标系需考虑协变导数。程序中已考虑，这里不再做公式介绍。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run_upar/run.sh + :language: bash + :start-after: BEGIN STRESS + :end-before: END STRESS + + 在 :rst:dir:`syn_dc_zne/` 原路径下，生成 :file:`*.stress.??.sac`，文件名中包括分量名，如ZZ、ZN等。 + + .. tab:: Python + + .. literalinclude:: run_upar/run.py + :language: python + :start-after: BEGIN STRESS + :end-before: END STRESS + + 返回的 |Stream| 通道名即为分量名，如ZZ、ZN等。 + +.. image:: run_upar/stress.png + :align: center + + +.. note:: + + 计算得到的应力的单位为 :math:`\frac{\text{dyne}}{\text{cm}^2} = 0.1 \text{Pa}`。 + + +由于场点位于地表（自由表面），过Z平面的应力均为0（由于浮点数计算误差，呈极小非零数），结果和理论保持一致。 \ No newline at end of file diff --git a/docs/source/Tutorial/dynamic/syn.rst b/docs/source/Tutorial/dynamic/syn.rst index 1ba0eb84..0eea0ea5 100644 --- a/docs/source/Tutorial/dynamic/syn.rst +++ b/docs/source/Tutorial/dynamic/syn.rst @@ -1,3 +1,5 @@ +.. _syn_rst: + 合成动态位移 ================= @@ -6,9 +8,31 @@ ----------------------------------------------------------- +Python中合成动态位移的主函数为 :func:`gen_syn_from_gf_*() ` （\*表示对不同震源） ，C程序为 :command:`grt.syn`。 使用上节计算的格林函数，合成动态位移（理论地震图）。方便起见，这里统一使用milrow模型，震源深度2km，场点位于地表，震中距10km的格林函数，方位角30°。 +在已知三分量格林函数 :math:`W_m(t), Q_m(t), V_m(t)` 后，合成三分量位移 :math:`u_z(t), u_r(t), u_\theta (t)` 的公式为 + +.. math:: + + \left\{ + \begin{aligned} + u_z(t) &= D(t) * \left[ \sum_{m=0}^{m=2} A_m W_m(t) \right] \\ + u_r(t) &= D(t) * \left[ \sum_{m=0}^{m=2} A_m Q_m(t) \right] \\ + u_\theta (t) &= D(t) * \left[ \sum_{m=1}^{m=2} A_{m+3} V_m(t) \right] + \end{aligned} + \right. + +其中 :math:`D(t)` 为震源时间函数，:math:`*` 表示卷积，:math:`A_m` 为与方位角和震源机制相关的方向因子，其中 :math:`u_z, u_r` 的方向因子相同，而 :math:`u_\theta` 的方向因子满足 + +.. math:: + + A_{m+3} = \frac{d A_m}{d (m\theta)}, m=1,2 + +其中 :math:`m` 为阶数，:math:`\theta` 为方位角。 + + .. note:: 合成位移的结果单位为 :math:`\text{cm}`。 diff --git a/docs/source/Tutorial/static/run_upar/milrow b/docs/source/Tutorial/static/run_upar/milrow new file mode 100644 index 00000000..a99f4912 --- /dev/null +++ b/docs/source/Tutorial/static/run_upar/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 +1. 6.9 4.0 2.8 9e10 9e10 +0.0 8.2 4.7 3.2 9e10 9e10 \ No newline at end of file diff --git a/docs/source/Tutorial/static/run_upar/run.py b/docs/source/Tutorial/static/run_upar/run.py new file mode 100644 index 00000000..996bf0fd --- /dev/null +++ b/docs/source/Tutorial/static/run_upar/run.py @@ -0,0 +1,62 @@ +import numpy as np +import pygrt + +modarr = np.loadtxt("milrow") + +pymod = pygrt.PyModel1D(modarr, depsrc=2.0, deprcv=0.0) + +xarr = np.linspace(-3, 3, 41) +yarr = np.linspace(-2.5, 2.5, 33) +# 传入calc_upar=True可计算空间导数 +static_grn = pymod.compute_static_grn(xarr, yarr, calc_upar=True) + +# 传入calc_upar=True可计算空间导数 +# 传入ZNE=True返回ZNE分量 +static_syn = pygrt.utils.gen_syn_from_gf_DC(static_grn, M0=1e24, strike=33, dip=50, rake=120, ZNE=True, calc_upar=True) + +# 计算应变 +static_strain = pygrt.utils.compute_strain(static_syn) + +# 计算应力 +static_stress = pygrt.utils.compute_stress(static_syn) + + +# 简单绘制结果 +import matplotlib.pyplot as plt + +def plot6(data:dict, title:str, out:str|None=None): + chs = ['ZZ', 'ZN', 'ZE', 'NN', 'NE', 'EE'] + xarr = data['_xarr'] + yarr = data['_yarr'] + fig, axs = plt.subplots(2, 3, figsize=(10, 6)) + axs = axs.ravel() + + MAX = 0 + for i in range(6): + ch = chs[i] + m = np.max(np.abs(data[ch])) + if m > MAX: + MAX = m + + for i in range(6): + ax = axs[i] + ch = chs[i] + vmin = vmax = None + if np.max(np.abs(data[ch]))/MAX < 1e-5: + vmin = -1 + vmax = 1 + + pcm = ax.pcolormesh(yarr, xarr, data[ch], shading='nearest', vmin=vmin, vmax=vmax) + ax.set_title(ch) + cbar = fig.colorbar(pcm, ax=ax) + cbar.formatter.set_powerlimits((0, 0)) + cbar.update_normal(pcm) + + fig.suptitle(title) + + if out is not None: + fig.tight_layout() + fig.savefig(out, dpi=100) + +plot6(static_strain, "Strain", 'static_strain.png') +plot6(static_stress, "Stress", 'static_stress.png') \ No newline at end of file diff --git a/docs/source/Tutorial/static/run_upar/run.sh b/docs/source/Tutorial/static/run_upar/run.sh new file mode 100755 index 00000000..09a15ff1 --- /dev/null +++ b/docs/source/Tutorial/static/run_upar/run.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +depsrc=2 +deprcv=0 + +x1=-3 +x2=3 +nx=41 + +y1=-2.5 +y2=2.5 +ny=33 +# -e 表示计算空间导数 +stgrt -Mmilrow -D${depsrc}/${deprcv} -X$x1/$x2/$nx -Y$y1/$y2/$ny -e > grn + +# -N 表示输出ZNE分量 +stgrt.syn -S1e24 -M33/50/120 -e -N < grn > syn_dc_zne + +# 计算应变 +stgrt.strain < syn_dc_zne > strain + +# 计算应力 +stgrt.stress < syn_dc_zne > stress \ No newline at end of file diff --git a/docs/source/Tutorial/static/run_upar/static_strain.png b/docs/source/Tutorial/static/run_upar/static_strain.png new file mode 100644 index 00000000..ef3802b1 Binary files /dev/null and b/docs/source/Tutorial/static/run_upar/static_strain.png differ diff --git a/docs/source/Tutorial/static/run_upar/static_stress.png b/docs/source/Tutorial/static/run_upar/static_stress.png new file mode 100644 index 00000000..c819468b Binary files /dev/null and b/docs/source/Tutorial/static/run_upar/static_stress.png differ diff --git a/docs/source/Tutorial/static/static_strain_stress.rst b/docs/source/Tutorial/static/static_strain_stress.rst index 42e34eac..aba1e09d 100644 --- a/docs/source/Tutorial/static/static_strain_stress.rst +++ b/docs/source/Tutorial/static/static_strain_stress.rst @@ -6,4 +6,31 @@ ----------------------------------------------------------- -TODO +除了使用不同的函数名/程序名，输出文件不同之外，流程基本和 :ref:`strain_stress_rst` 类似。这里直接给出脚本。运行后注意观察和 :ref:`static_syn_rst` 的输出有了什么变化。 + +.. tabs:: + + .. tab:: C + + .. literalinclude:: run_upar/run.sh + :language: bash + + + .. tab:: Python + + .. literalinclude:: run_upar/run.py + :language: python + +------------- + +.. image:: run_upar/static_strain.png + :align: center + + +------------- + +.. image:: run_upar/static_stress.png + :align: center + + +由于场点位于地表（自由表面），过Z平面的应力均为0（由于浮点数计算误差，呈极小非零数），结果和理论保持一致。 \ No newline at end of file diff --git a/docs/source/Tutorial/static/static_syn.rst b/docs/source/Tutorial/static/static_syn.rst index a2c85325..3345dbe5 100644 --- a/docs/source/Tutorial/static/static_syn.rst +++ b/docs/source/Tutorial/static/static_syn.rst @@ -1,3 +1,5 @@ +.. _static_syn_rst: + 合成静态位移 ================= diff --git a/docs/source/index.rst b/docs/source/index.rst index 426ade16..00fe3a0d 100755 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -58,6 +58,8 @@ 以下包括 **Python** 和 **C** 两部分的API介绍。如果只关心程序的使用，可以只查看 **Python** 部分；若对具体实现过程感兴趣，尤其是在学习 **广义反射透射系数矩阵、离散波数法、线性Filon积分法、峰谷平均法** 等内容感到迷惑难解时，可以在运行程序的同时结合阅读相关 **C代码**，我在C代码中给出了 **明确的中文注释以及相关公式索引** ，希望能帮助你在学习相关内容时解答一些疑惑。 +本文档页面源码位于 Github主页的 `docs/ `_，可查看示例代码和画图脚本。 + .. toctree::