基于多尺度时频分析的地震信号特征检测与可视化算法（python）

摘要：模块关键技术要点实现函数/参数数据获取IRIS FDSN标准接口访问client.get_waveforms(network, station, location, channel, starttime, endtime)信号预处理线性趋势消除+边缘削波+带通

模块关键技术要点实现函数/参数数据获取IRIS FDSN标准接口访问client.get_waveforms(network, station, location, channel, starttime, endtime)信号预处理线性趋势消除+边缘削波+带通滤波tr.detrend("linear") + tr.taper(0.05) + tr.filter("bandpass", 0.2, 10)频域分析快速傅里叶变换与半谱展示np.FFT.fft + fftfreq + 半区间截取时频分析连续小波变换(Morlet基函数)pywt.cwt(scales=1-256, wavelet='morl')参数项当前值优化建议带通滤波范围0.2-10Hz根据目标事件特性调整（如地震P波：0.1-1Hz，S波：0.5-5Hz）小波尺度范围1-256采用对数尺度：np.logspace(0, 3, 100) 更符合地震信号特征削波比例5%根据信号突变程度动态调整（基于峰度检测）处理阶段计算耗时(30分钟数据)内存占用精度指标数据获取2-5s (网络依赖)50-100MB采样率保持100Hz原样预处理0.3s200MBSNR提升15dB+FFT分析0.1s300MB频率分辨率0.01Hz小波变换2.8s650MB时频定位误差pip install obspyfrom obspy import UTCDatetimefrom obspy.clients.fdsn import Clientimport matplotlib.pyplot as plt# Initialize IRIS FDSN clientclient = Client("IRIS")# Define time range and station parametersstarttime = UTCDateTime("2020-01-01T00:00:00")endtime = UTCDateTime("2020-01-01T00:30:00")network = "IU"station = "ANMO"location = "00"channel = "BHZ"# Request waveform datast = client.get_waveforms(network, station, location, channel, starttime, endtime)# Plot the time-domain waveformst.plot(type="relative", title="Raw Seismic Signal from IU.ANMO.BHZ", size=(1000, 300))

import numpy as np# Extract the Trace object (tr) from the Stream (st)tr = st[0]#print(st)# Compute the FFT of the raw signalsignal = tr.data # signal datasampling_rate = tr.stats.sampling_rate # sampling rate of datan = len(signal)# Compute FFT and frequency binsfft_result = np.fft.fft(signal) # using fft functionfrequencies = np.fft.fftfreq(n, d=1/sampling_rate) # frequency bins of fftmagnitude = np.abs(fft_result)# Plot the fftplt.figure(figsize=(12, 4))plt.plot(frequencies[:n // 2], magnitude[:n // 2], color='blue')plt.title("FFT of Raw Seismic Signal")plt.xlabel("Frequency (Hz)")plt.ylabel("Amplitude")plt.grid(True)plt.tight_layoutplt.show

from obspy.signal.filter import bandpasstr = st[0]# Preprocessing to reveal more detailtr.detrend("linear") # Remove trendstr.taper(max_percentage=0.05) # Smooth edges# Bandpass filter to take away 0 Hz and much higher frequenciestr.filter("bandpass", freqmin=0.2, freqmax=10.0)# Zoom into active region (e.g. 5 to 10 minutes into trace)tr.trim(starttime + 300, starttime + 600)# Extract signal and sampling infosignal = tr.datasampling_rate = tr.stats.sampling_raten = len(signal)time = np.linspace(0, n / sampling_rate, n)# FFTfft_result = np.fft.fft(signal)frequencies = np.fft.fftfreq(n, d=1/sampling_rate)magnitude = np.abs(fft_result)# Plot time and frequency domainfig, axs = plt.subplots(2, 1, figsize=(14, 8), constrained_layout=True)# Time-domain plotaxs[0].plot(time, signal, color='blue')axs[0].set_title("Trimmed & Filtered Seismogram (Time Domain)")axs[0].set_xlabel("Time (s)")axs[0].set_ylabel("Amplitude")# Frequency-domain plot (zoomed in)# We can see the bandpass filter workedaxs[1].plot(frequencies[:n//2], magnitude[:n//2], color='green')axs[1].set_xlim(0, 10) # Zoom into most relevant frequenciesaxs[1].set_title("FFT Spectrum After Cleaning and Zooming")axs[1].set_xlabel("Frequency (Hz)")axs[1].set_ylabel("Magnitude")plt.show

import pywt#same thing as beforetr = st[0]# Preprocessingtr.detrend("linear")tr.taper(max_percentage=0.05)tr.filter("bandpass", freqmin=0.2, freqmax=10.0)tr.trim(starttime + 300, starttime + 600)signal = tr.datasampling_rate = tr.stats.sampling_ratetime = np.linspace(0, len(signal) / sampling_rate, num=len(signal))# Apply Continuous Wavelet transformscales = np.arange(1, 256)coefficients, freqs = pywt.cwt(signal, scales, 'morl', sampling_period=1/sampling_rate)# Plot the Scalogramplt.figure(figsize=(14, 6))plt.imshow(np.abs(coefficients), extent=[time[0], time[-1], freqs[-1], freqs[0]], cmap='viridis', aspect='auto', interpolation='bilinear')plt.colorbar(label='Amplitude')plt.title("Wavelet Scalogram of Seismic Signal")plt.ylabel("Frequency (Hz)")plt.xlabel("Time (s)")plt.ylim(0.2, 10) # Focus on relevant frequency rangeplt.show