爆破地震波信号处理HHT改进算法及应用【附代码】
✨ 长期致力于爆破地震波、希尔伯特-黄变换、经验模态分解、固有模态函数、模态混淆、端点效应研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》1CEEMD-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_ensembles50, noise_std0.2): ceemdan CEEMDAN() imfs ceemdan(signal) mpe_values [] for imf in imfs: # 多尺度排列熵简化计算 mpe 0 scale 3 for tau in range(1, scale1): embedded np.array([imf[i:len(imf)-tau1] for i in range(tau)]) patterns np.lexsort(embedded) unique, counts np.unique(patterns, return_countsTrue) probs counts / len(patterns) mpe -np.sum(probs * np.log(probs1e-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_extend10): # 左端延拓 left_peaks signal[:20] # 拟合指数曲线 x np.arange(len(left_peaks)) coeffs np.polyfit(x, np.log(left_peaks1e-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_peaks1e-8), 1) extended_right np.exp(coeffs_r[1]) * np.exp(coeffs_r[0] * np.arange(1, n_extend1)) 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, fs1000): 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, fs10000): imfs ceemd_mpe(blasting_signal) # 选择主频带IMF main_imf imfs[1] # 第二分量一般为爆破主频 # 检测峰值 peaks, _ signal.find_peaks(np.abs(main_imf), height0.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)}) ,