地震信号降噪与事件提取（Python）

摘要：import pandas as pdimport matplotlib.pyplot as plt# Load your datasetdata = pd.read_csv('xa.s12.00.mhz.1970-01-19HR00_evid00002exa

import pandas as pdimport matplotlib.pyplot as plt# Load your datasetdata = pd.read_csv('xa.s12.00.mhz.1970-01-19HR00_evid00002example.csv')# Display the first few rows to understand its structureprint(data.head)# The dataset tested had the columns named as "time_rel(sec)" and veclocity(sec), change as required time_column = 'time_rel(sec)' value_column = 'velocity(m/s)' # Plot the original dataplt.figure(figsize=(12, 6))plt.plot(data[time_column], data[value_column], color='blue', label='Seismic Data')plt.title('Seismic Data Plot')plt.xlabel('Time (seconds)')plt.ylabel('Amplitude')plt.legendplt.gridplt.show time_abs(%Y-%m-%dT%H:%M:%S.%f) time_rel(sec) velocity(m/s)0 1970-01-19T00:00:00.665000 0.000000 -6.150000e-141 1970-01-19T00:00:00.815943 0.150943 -7.700000e-142 1970-01-19T00:00:00.966887 0.301887 -8.400000e-143 1970-01-19T00:00:01.117830 0.452830 -8.100000e-144 1970-01-19T00:00:01.268774 0.603774 -7.100000e-14

import pandas as pdimport numpy as npimport pywt# Apply Wavelet transform to the relevant column (adjust the column name as needed)coeffs = pywt.wavedec(data['velocity(m/s)'], 'db4', level=6)# Create time absolute and relative columnstime_abs = data['time_abs(%Y-%m-%dT%H:%M:%S.%f)']time_rel = (pd.to_datetime(time_abs) - pd.to_datetime(time_abs.iloc[0])).dt.total_seconds# Create a dictionary to store wavelet coefficients by level (VERYYYYYY IMPORTANT)wavelet_data = {'time_abs(%Y-%m-%dT%H:%M:%S.%f)': time_abs,'time_rel(sec)': time_rel}# Add each level's wavelet coefficients to the DataFramefor i, coeff in enumerate(coeffs):wavelet_data[f'wavelet_coefficients_level_{i}'] = np.pad(coeff, (0, len(time_abs) - len(coeff)), mode='constant', constant_values=np.nan)# Create the DataFrame with the wavelet coefficients organized by levelwavelet_df = pd.DataFrame(wavelet_data)# Save the DataFrame to a CSV filewavelet_df.to_csv('wavelet_coefficients_by_level.csv', index=False)print("Wavelet coefficients by level saved to 'wavelet_coefficients_by_level.csv'.")Wavelet coefficients by level saved to 'wavelet_coefficients_by_level.csv'.import matplotlib.pyplot as plt# Create a figure to plot the coefficientsplt.figure(figsize=(14, 10))# Plot each level of coefficientsfor i, coeff in enumerate(coeffs):plt.subplot(len(coeffs), 1, i + 1)plt.plot(coeff, label=f'Wavelet Coefficients - Level {i}')plt.title(f'Wavelet Coefficients - Level {i}')plt.xlabel('Sample Index')plt.ylabel('Coefficient Value')plt.legendplt.tight_layoutplt.show

import pandas as pdimport numpy as np# Load the wavelet coefficients CSV filewavelet_df = pd.read_csv('wavelet_coefficients_by_level.csv')# List to hold coefficient arrays for each levelcoeffs = max_levels = 7 # Assuming you want to check for levels 0 to 6 (the 'resolution' so to say of this program goes from 0-6)# Check and create arrays for each levelfor level in range(max_levels):col_name = f'wavelet_coefficients_level_{level}'if col_name in wavelet_df.columns:# Get the coefficient array for this levelcoeff_array = wavelet_df[col_name].values# Find the last non-empty value's indexlast_valid_index = np.where(~np.isnan(coeff_array))[0][-1] if np.any(~np.isnan(coeff_array)) else -1if last_valid_index >= 0:# Slice the coefficient array to keep only the valid valuestrimmed_coeff_array = coeff_array[:last_valid_index + 1]coeffs.append(trimmed_coeff_array)print(f'Shape of {col_name}: {trimmed_coeff_array.shape}') # Print the shape of each coefficient arrayelse:print(f'No valid values found in {col_name}.')else:print(f'Column {col_name} does not exist in the DataFrame.')# At this point, you have coeffs as a list containing arrays for each level#Arrays are padded to avoid dimension errorsShape of wavelet_coefficients_level_0: (8950,)Shape of wavelet_coefficients_level_1: (8950,)Shape of wavelet_coefficients_level_2: (17894,)Shape of wavelet_coefficients_level_3: (35782,)Shape of wavelet_coefficients_level_4: (71558,)Shape of wavelet_coefficients_level_5: Shape of wavelet_coefficients_level_6: import pandas as pdimport numpy as npimport pywtimport matplotlib.pyplot as pltfrom scipy.signal import find_peaksfrom scipy.ndimage import gaussian_filter1d# Load wavelet coefficientswavelet_df = pd.read_csv('wavelet_coefficients_by_level.csv')# Extract coefficients for levels 4, 5, and 6 (most common for quake levels)level_4 = wavelet_df['wavelet_coefficients_level_4'].dropna.valueslevel_5 = wavelet_df['wavelet_coefficients_level_5'].dropna.valueslevel_6 = wavelet_df['wavelet_coefficients_level_6'].dropna.values# Ensure the same length for all levels (padding)min_length = min(len(level_4), len(level_5), len(level_6))level_4 = level_4[:min_length]level_5 = level_5[:min_length]level_6 = level_6[:min_length]# Combine coefficientscombined_coeffs = level_4 + level_5 + level_6# Smooth the combined datasmoothed_data = gaussian_filter1d(combined_coeffs, sigma=2)# Create time series for plottingtime_rel = wavelet_df['time_rel(sec)'].dropna.values[:min_length]# Set a higher threshold for peak detectionmean_value = np.mean(smoothed_data)std_value = np.std(smoothed_data)######### for the start point# Increase the multiplier significantly to further reduce false positivesthreshold_value = mean_value + 4 * std_value # Try increasing to 4 or 5 (Change as required)minimum_peak_height = threshold_value # Set a minimum height for detected peaksminimum_peak_width = 10 # Set minimum width for detected peaks# Find peaks to mark the start of seismic eventspeaks, properties = find_peaks(smoothed_data, height=minimum_peak_height, distance=50, width=minimum_peak_width)# Check if any peaks are found and use the first one for markingevent_start = peaks[0] if len(peaks) > 0 else None########## for the end point# Set a threshold for end detectionthreshold = mean_value - 2.5 * std_value # Adjust this value to your data characteristics (same as above)# Start from the end of the signal and look for the end of the eventfor i in range(min_length - 1, 0, -1):# Check if the amplitude is below the threshold for stabilityif smoothed_data[i]

Smoothed wavelet coefficients saved to 'smoothed_wavelet_coefficients.csv_with_start_and_end'.The suspected start of the seismic event is at: 9187.924528 secondsThe suspected end of the seismic event is at: 9676.830189 secondsimport pandas as pdimport numpy as npimport pywtimport matplotlib.pyplot as pltfrom scipy.signal import find_peaksfrom scipy.ndimage import gaussian_filter1d# Load wavelet coefficientswavelet_df = pd.read_csv('wavelet_coefficients_by_level.csv')# Extract coefficients for levels 4, 5, and 6level_4 = wavelet_df['wavelet_coefficients_level_4'].dropna.valueslevel_5 = wavelet_df['wavelet_coefficients_level_5'].dropna.valueslevel_6 = wavelet_df['wavelet_coefficients_level_6'].dropna.values# Ensure the same length for all levels (padding)min_length = min(len(level_4), len(level_5), len(level_6))level_4 = level_4[:min_length]level_5 = level_5[:min_length]level_6 = level_6[:min_length]# Combine coefficientscombined_coeffs = level_4 + level_5 + level_6# Smooth the combined datasmoothed_data = gaussian_filter1d(combined_coeffs, sigma=2)# Create time series for plottingtime_rel = wavelet_df['time_rel(sec)'].dropna.values[:min_length]# Set a higher threshold for peak detectionmean_value = np.mean(smoothed_data)std_value = np.std(smoothed_data)# Increase the multiplier significantly to further reduce false positivesthreshold_value = mean_value + 4 * std_value # Try increasing to 4 or 5minimum_peak_height = threshold_value # Set a minimum height for detected peaksminimum_peak_width = 10 # Set minimum width for detected peaks# Find peaks to mark the start of seismic eventspeaks, properties = find_peaks(smoothed_data, height=minimum_peak_height, distance=50, width=minimum_peak_width)# Check if any peaks are found and use the first one for markingevent_start = peaks[0] if len(peaks) > 0 else None# Set a threshold for end detectionthreshold = mean_value - 2.5 * std_value # Adjust this value to your data characteristics# Start from the end of the signal and look for the end of the eventfor i in range(min_length - 1, 0, -1):# Check if the amplitude is below the threshold for stabilityif smoothed_data[i]

Modified wavelet coefficients saved to 'modified_wavelet_coefficients.csv'.
The suspected start of the seismic event is at: 9187.924528 seconds
The suspected end of the seismic event is at: 9676.830189 seconds

知乎学术咨询：

https://www.zhihu.com/consult/people/792359672131756032?isMe=1担任《Mechanical System and Signal Processing》《中国电机工程学报》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测。

分割线分割线

MATLAB环境下一种改进的瞬时频率（IF）估计方法