🎯要点
- 在二维均匀介质均匀源中合成互相关函数以便构建波层析成像。
- 闭环系统中微积分计算情景：完美弹性体震波、随机外力对模式的能量分配。
- 开环系统中微积分计算情景：无数震源激发波方程、闭合曲线上的随机源、不相关平面波事件。
- 整理地震波场频谱数据，选择局部平稳状态数据，对数据归一化去噪平滑。
- 通过频率-时间分析等方法了解频散曲线。
📜地震波信号解析算法及深度学习
-
Python地震波逆问题解构算法复杂信号分析
-
Python空间地表联动贝叶斯地震风险计算模型
-
Python噪声敏感度和沉积区震动波
-
互相关导图
🍪语言内容分比
🍇Python瑞利波
瑞利波是一种沿固体表面传播的表面声波。它们可以通过多种方式在材料中产生，例如通过局部冲击或压电传导，并且经常用于无损检测以检测缺陷。瑞利波是地震在地球上产生的地震波的一部分。当它们在层中传播时，它们被称为兰姆波、瑞利-兰姆波或广义瑞利波。
在由拉梅参数
λ
\lambda
λ 和
μ
\mu
μ 描述的各向同性线性弹性材料中，瑞利波的速度由以下方程的解给出：
ζ
3
−
8
ζ
2
+
8
ζ
(
3
−
2
η
)
−
16
(
1
−
η
)
=
0
\zeta^3-8 \zeta^2+8 \zeta(3-2 \eta)-16(1-\eta)=0
ζ3−8ζ2+8ζ(3−2η)−16(1−η)=0
其中
ζ
=
ω
2
/
k
2
β
2
,
η
=
β
2
/
α
2
,
ρ
α
2
=
λ
+
2
μ
\zeta=\omega^2 / k^2 \beta^2, \eta=\beta^2 / \alpha^2, \rho \alpha^2=\lambda+2 \mu
ζ=ω2/k2β2,η=β2/α2,ρα2=λ+2μ, 且
ρ
β
2
=
μ
\rho \beta^2=\mu
ρβ2=μ 。由于此方程没有固有尺度，因此产生瑞利波的边界值问题无弥散性。一个有趣的特殊情况是泊松固体，其
λ
=
μ
\lambda=\mu
λ=μ，因为这会产生一个与频率无关的相速度，等于
ω
/
k
=
β
0.8453
\omega / k=\beta \sqrt{0.8453}
ω/k=β0.8453。对于具有正泊松比
(
ν
>
0.3
)
(\nu>0.3)
(ν>0.3)的线弹性材料，瑞利波速可近似为
c
R
=
c
S
0.862
+
1.14
ν
1
+
ν
c_R=c_S \frac{0.862+1.14 \nu}{1+\nu}
cR=cS1+ν0.862+1.14ν，其中
c
S
c_S
cS是剪切波速度。
由于材料性质的变化，弹性常数通常会随深度而变化。这意味着瑞利波的速度实际上取决于波长（因此也取决于频率），这种现象称为频散。受频散影响的波具有不同的波列形状。如上所述，理想、均匀和平坦的弹性固体上的瑞利波没有频散。但是，如果固体或结构的密度或声速随深度而变化，瑞利波就会变得频散。一个例子是地球表面的瑞利波：频率较高的波比频率较低的波传播得更慢。这是因为较低频率的瑞利波具有相对较长的波长。长波长波的位移比短波长波更深地穿透地球。由于地球中的波速随着深度的增加而增加，较长波长（低频）波可以比较短波长（高频）波传播得更快。因此，瑞利波在远处地震记录站记录的地震图上经常显得分散。还可以观察薄膜或多层结构中的瑞利波弥散。
瑞利波椭相移Python片段
import pickle
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FuncFormatter, MultipleLocator
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from obspy import Stream, UTCDateTime
from obspy.signal.rotate import rotate2zne
from obspy.clients.filesystem.sds import Client as sdsClient
from scipy.signal import correlation_lags, correlate
from scipy.signal import hilbert
def corr_phaseshift(zaux, zaur, period, srate):
corrs_ph = correlate(zaux, zaur)
lags = correlation_lags(len(zaux), len(zaur))
corrs_ph /= np.max(corrs_ph)
maxshift = round(period*srate/2)
minshift = -maxshift
boolsh = np.logical_and(lags>=minshift, lags<=maxshift)
newcorr = corrs_ph[boolsh]
newlag = lags[boolsh]
argmaxc = np.argmax(newcorr)
best_lag = newlag[argmaxc]
phase_shift = best_lag*360/(period*srate)
return phase_shift, best_lag
def calc_vars(radial, vertical, period, ref_time):
odict = {'ph_shift': [], 'char_funs': [],
'times_cf': [], 'hvs': [], 'ccs': []}
phs = np.arange(0, 182, 2)
lags = np.unique((period*vertical.stats.sampling_rate*phs/360).astype(int))
for lag in lags:
ph = lag*360/(period*vertical.stats.sampling_rate)
wwz = vertical.copy()
wwr = radial.copy()
wwz.data = wwz.data[lag:]
if lag!=0:
wwr.data = wwr.data[:-lag]
time_arr = wwr.times() + (wwr.stats.starttime - ref_time)
wlen = int(period)
step = int(round(period/5, 0))
if step<1:
step = 1
if wlen<1:
wlen = 1
t0win = []
corrs = []
aux_zr = Stream(traces=[wwz, wwr])
for window in aux_zr.slide(window_length=wlen, step=step):
auxz = window.select(component='Z')[0].data
auxr = window.select(component='R')[0].data
norm = (np.sum(auxz**2)*np.sum(auxr**2))**0.5
fcc = obcorrelate(auxz, auxr, shift=0, demean=True, normalize=None)/norm
shift, valuecc = xcorr_max(fcc)
t0win.append(window[0].stats.starttime-ref_time)
corrs.append(valuecc)
return(odict)
def grid_data(inx, inzdata):
common_x = inx[0]
length = common_x.size
in_x = inx.copy()
in_zdata = inzdata.copy()
for it, etime in enumerate(in_x):
thisl = len(etime)
in_x[it] = np.concatenate([etime, np.nan*np.ones(length-thisl)])
in_zdata[it] = np.concatenate([in_zdata[it], np.nan*np.ones(length-thisl)])
return(in_x, in_zdata)
def shift_unc(xvals, yvals, zgrid, per, evname, thres=0.8, save=False):
percs = [16, 50, 84]
rd = 0.01
xt, yph = np.meshgrid(xvals, yvals)
newxt, newyph = xt.flatten(), yph.flatten()
points = np.vstack((newxt, newyph)).T
fig, ax = plt.subplots()
CS = ax.contour(xt, yph, zgrid, levels=[0.0, 0.4, 0.6, 0.8],
colors='gray', alpha=0.7)
maxval = np.nanmax(zgrid)
thres_unc = 0.95*maxval
CSun = ax.contour(xt, yph, zgrid, levels=[thres_unc],
colors='black', alpha=0.8)
p_un = CSun.collections[0].get_paths()[-1]
bool_sel = p_un.contains_points(points, radius=rd)
ph_min, ph_med, ph_max = np.percentile(newyph[bool_sel], percs)
xt_min, xt_med, xt_max = np.percentile(newxt[bool_sel], percs)
if maxval>=thres:
CSthres = ax.contour(xt, yph, zgrid, levels=[thres],
colors='red', alpha=0.8, linestyles='dashed')
for c, contour in enumerate(CSthres.collections[0].get_paths()):
if any(contour.contains_points(p_un.vertices)) or\
any(p_un.contains_points(contour.vertices)):
break
if thres_unc<thres:
new_bool = contour.contains_points(points, radius=rd)
ph_min, ph_med, ph_max = np.percentile(newyph[new_bool], percs)
xt_min, xt_med, xt_max = np.percentile(newxt[new_bool], percs)
inph = np.logical_and(newyph<=ph_max, newyph>=ph_min)
incon = contour.contains_points(points, radius=rd)
selpoints = points[np.logical_and(inph, incon)]
minTmax = [ np.nanpercentile(selpoints[selpoints.T[1]==auxph].T[0], [0, 100])
for auxph in np.unique(selpoints.T[1]) ]
_, med_minT, _ = np.percentile(np.asarray(minTmax).T[0], percs)
_, med_maxT, _ = np.percentile(np.asarray(minTmax).T[1], percs)
else:
print(f'Max value not exceeding the threshold [{maxval}/{thres}]')
med_minT = np.nan
med_maxT = np.nan
unc_dict = dict(min_ph=ph_min, med_ph=ph_med, max_ph=ph_max,
minT_thres=med_minT, maxT_thres=med_maxT,
maxchf_mint=xt_min, maxchf_medt=xt_med, maxchf_maxt=xt_max)
if save:
ax.clabel(CS, CS.levels, inline=True,
fmt=FuncFormatter(lambda y,_: '{:g}'.format(y)), fontsize=10)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Phase shift [$^o$]')
dict_unc = dict(color='red', alpha=0.6, ls='--')
for phs in [ph_min, ph_med, ph_max]:
ax.axhline(phs, **dict_unc)
for xts in [med_minT, med_maxT, xt_min, xt_med, xt_max]:
ax.axvline(xts, **dict_unc)
ax.yaxis.set_major_locator(MultipleLocator(30))
ax.yaxis.set_minor_locator(MultipleLocator(10))
ax.grid(which='major', ls='--', alpha=0.5, color='gray')
ax.text(0.01, 0.98, 'T$_0$ = ' + f'{per} s \n{evname}', ha='left',
va='top', fontweight='bold', fontsize=10, transform=ax.transAxes)
ax.text(0.01, ph_med, '$\phi^*$ = ' + f'{ph_med:.1f}' + '$^o$', ha='right',
va='bottom', fontweight='bold', color='red',
fontsize=10, transform=ax.get_yaxis_transform())
fig.tight_layout()
fig.savefig(f'{evname}_{per}s_chf.png', dpi=200)
plt.close(fig)
return(unc_dict)
def shift_waves(auxz, auxr, blag, time_ref):
if blag>=0:
newz = auxz.data[blag:]
if blag==0:
newr = auxr.data
arr_time = auxr.times() + (auxr.stats.starttime - time_ref)
else:
newr = auxr.data[:-blag]
arr_time = auxr.times()[:-blag] + (auxr.stats.starttime - time_ref)
else:
newz = auxz.data[:blag]
newr = auxr.data[-blag:]
arr_time = auxz.times()[:blag] + (auxz.stats.starttime - time_ref)
unsh_zz = auxz.copy()
auxz.data = newz
auxr.data = newr
return unsh_zz, arr_time
