Python微震波频散相速分析

🎯要点

1. 在二维均匀介质均匀源中合成互相关函数以便构建波层析成像。
2. 闭环系统中微积分计算情景：完美弹性体震波、随机外力对模式的能量分配。
3. 开环系统中微积分计算情景：无数震源激发波方程、闭合曲线上的随机源、不相关平面波事件。
4. 整理地震波场频谱数据，选择局部平稳状态数据，对数据归一化去噪平滑。
5. 通过频率-时间分析等方法了解频散曲线。

📜地震波信号解析算法及深度学习

• Python地震波逆问题解构算法复杂信号分析

• Python空间地表联动贝叶斯地震风险计算模型

• Python噪声敏感度和沉积区震动波

• 互相关导图
  

🍪语言内容分比

🍇Python瑞利波

瑞利波是一种沿固体表面传播的表面声波。它们可以通过多种方式在材料中产生，例如通过局部冲击或压电传导，并且经常用于无损检测以检测缺陷。瑞利波是地震在地球上产生的地震波的一部分。当它们在层中传播时，它们被称为兰姆波、瑞利-兰姆波或广义瑞利波。

在由拉梅参数 $\lambda$ 和 $\mu$ 描述的各向同性线性弹性材料中，瑞利波的速度由以下方程的解给出：
$${\zeta }^3-8{\zeta }^2+8\zeta (3-2\eta )-16(1-\eta )=0$$
其中 $\zeta ={\omega }^2/k^2{\beta }^2,\eta ={\beta }^2/{\alpha }^2,\rho {\alpha }^2=\lambda +2\mu$, 且 $\rho {\beta }^2=\mu$ 。由于此方程没有固有尺度，因此产生瑞利波的边界值问题无弥散性。一个有趣的特殊情况是泊松固体，其 $\lambda =\mu$，因为这会产生一个与频率无关的相速度，等于 $\omega /k=\beta \sqrt{0.8453}$。对于具有正泊松比$(\nu >0.3)$的线弹性材料，瑞利波速可近似为$c_R=c_S\frac{0.862+1.14\nu }{1+\nu }$，其中$c_S$是剪切波速度。

由于材料性质的变化，弹性常数通常会随深度而变化。这意味着瑞利波的速度实际上取决于波长（因此也取决于频率），这种现象称为频散。受频散影响的波具有不同的波列形状。如上所述，理想、均匀和平坦的弹性固体上的瑞利波没有频散。但是，如果固体或结构的密度或声速随深度而变化，瑞利波就会变得频散。一个例子是地球表面的瑞利波：频率较高的波比频率较低的波传播得更慢。这是因为较低频率的瑞利波具有相对较长的波长。长波长波的位移比短波长波更深地穿透地球。由于地球中的波速随着深度的增加而增加，较长波长（低频）波可以比较短波长（高频）波传播得更快。因此，瑞利波在远处地震记录站记录的地震图上经常显得分散。还可以观察薄膜或多层结构中的瑞利波弥散。

瑞利波椭相移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

👉参阅、更新：计算思维 | 亚图跨际