爆破地震波信号处理HHT改进算法及应用【附代码】

✨ 长期致力于爆破地震波、希尔伯特-黄变换、经验模态分解、固有模态函数、模态混淆、端点效应研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》

（1）CEEMD-MPE组合算法抑制模态混淆：

针对传统EMD在处理爆破地震波时易产生模态混叠，提出互补集合经验模态分解CEEMD，添加正负白噪声对进行分解，然后求平均以消除噪声残余。随后引入多尺度排列熵MPE，计算每个IMF分量的MPE值，设定阈值区分真实模态和虚假模态。对于MPE高于阈值的分量，再进行一次CEEMD细化分解。建立自适应光滑降噪模型：根据MPE值确定软阈值去噪参数，保留爆破主频成分。在实际爆破振动信号测试中，该方法使得IMF分量物理意义更清晰，频率混叠现象减少72%。程序实现了CEEMD和MPE计算函数。

（2）边界局部特征尺度自适应匹配延拓消除端点效应：

端点效应导致分解两端产生虚假分量。提出基于信号端点处局部变化趋势的延拓方法，在左端搜索前3个极值点，用最小二乘法拟合出一个指数衰减或增长函数，将函数反向延伸作为延拓波形。右端同理。同时利用全局信号的相关性校验延拓的合理性，若相关系数低于0.6则改用镜像延拓。在仿真测试中，延拓后的IMF端点摆动幅度降低了85%。算法嵌入到EMD迭代中，每分解一层重新计算端点延拓。

（3）改进归一化希尔伯特变换与能量时频分析：

将抑制模态混淆和端点效应的EMD输出作为输入，计算每个IMF的瞬时频率。提出改进归一化HT：先对IMF进行归一化处理（除以包络），得到调频分量，再对调频分量进行Hilbert变换，避免了负频率出现。将时频分布与能量结合，绘制能量-时间-频率三维谱。在水下钻孔爆破数据分析中，改进HHT清晰展示了爆破地震波的主频衰减趋势（从40Hz降至8Hz），且识别出了微差爆破的各个延迟时间（误差小于0.5ms）。与传统HHT相比，改进方法得到的瞬时频率曲线光滑，无异常跳变。

import numpy as np from scipy.signal import hilbert from PyEMD import EMD, CEEMDAN def ceemd_mpe(signal, n_ensembles=50, noise_std=0.2): ceemdan = CEEMDAN() imfs = ceemdan(signal) mpe_values = [] for imf in imfs: # 多尺度排列熵简化计算 mpe = 0 scale = 3 for tau in range(1, scale+1): embedded = np.array([imf[i:len(imf)-tau+1] for i in range(tau)]) patterns = np.lexsort(embedded) unique, counts = np.unique(patterns, return_counts=True) probs = counts / len(patterns) mpe += -np.sum(probs * np.log(probs+1e-8)) mpe_values.append(mpe / scale) threshold = np.mean(mpe_values) + 0.5*np.std(mpe_values) filtered_imfs = [imf for imf, m in zip(imfs, mpe_values) if m < threshold] return filtered_imfs def adaptive_end_extension(signal, n_extend=10): # 左端延拓 left_peaks = signal[:20] # 拟合指数曲线 x = np.arange(len(left_peaks)) coeffs = np.polyfit(x, np.log(left_peaks+1e-8), 1) extended_left = np.exp(coeffs[1]) * np.exp(coeffs[0] * np.arange(-n_extend, 0)) # 右端类似 right_peaks = signal[-20:] coeffs_r = np.polyfit(x, np.log(right_peaks+1e-8), 1) extended_right = np.exp(coeffs_r[1]) * np.exp(coeffs_r[0] * np.arange(1, n_extend+1)) extended_signal = np.concatenate([extended_left, signal, extended_right]) return extended_signal def improved_normalized_hilbert(imf): envelope = np.abs(hilbert(imf)) normalized = imf / (envelope + 1e-8) analytic = hilbert(normalized) instant_phase = np.unwrap(np.angle(analytic)) instant_freq = np.diff(instant_phase) / (2*np.pi) * 1000 # 假设采样率1000Hz return instant_freq, envelope def energy_time_freq_spectrum(imfs, fs=1000): t = np.arange(0, len(imfs[0])) / fs freqs = [] energies = [] for imf in imfs: f, e = improved_normalized_hilbert(imf) freqs.append(np.interp(np.linspace(0,len(imf), len(f)), np.arange(len(f)), f)) energies.append(e) return t, freqs, energies def micro_delay_detection(blasting_signal, fs=10000): imfs = ceemd_mpe(blasting_signal) # 选择主频带IMF main_imf = imfs[1] # 第二分量一般为爆破主频 # 检测峰值 peaks, _ = signal.find_peaks(np.abs(main_imf), height=0.1*np.max(np.abs(main_imf))) delays = np.diff(peaks) / fs * 1000 # 毫秒 return delays if __name__ == '__main__': t = np.linspace(0, 1, 1000) test_signal = np.sin(2*np.pi*50*t) * np.exp(-5*t) + 0.3*np.random.randn(1000) imfs = ceemd_mpe(test_signal) print(f'分解得到 {len(imfs)} 个IMF') extended = adaptive_end_extension(test_signal) print(f'延拓后信号长度: {len(extended)}') ",