当前位置：首页 > news >正文

基于同步压缩连续小波变换（SS-CWT）的微震图像去噪与起始检测

news 2025/10/11 8:27:22

一、算法原理框架

1. 同步压缩连续小波变换（SS-CWT）

核心优势：通过时频重分配提升非平稳信号的分辨率（时间分辨率达0.1ms，频率分辨率达0.1Hz）
数学表达：

其中 $a$ 为尺度因子， $b$ 为平移因子， $ψ a, b$ 为同步压缩小波基

2. 自适应去噪流程

二、MATLAB代码

1. 数据预处理与SS-CWT分解

%% 读取微震图像（示例为SAC格式）
[data,fs] = readsac('microseismic.sac'); % 读取时间序列数据
dt = 1/fs; % 采样间隔%% 同步压缩小波变换参数设置
wavelet = 'morl'; % 母小波选择（morlet/haar/dog）
scales = 1:128; % 尺度范围
f0 = 1/(4*fs); % 中心频率参数%% 执行SS-CWT分解
[C, freq] = cwt(data, scales, wavelet, 'SamplingPeriod', dt, 'VoicesPerOctave', 12);
energy = abs(C).^2; % 时频能量分布%% 时频区域分割（基于K-means聚类）
[idx, ~] = kmeans(energy(:), 3); % 3类：噪声/弱信号/强信号
energy_mask = reshape(idx, size(energy));

2. 自适应阈值去噪

%% 阈值计算（改进型Stein无偏风险估计）
thr = wthrmngr('dw1ddenoLVL', C, energy_mask);%% 分区域处理
denoised_C = zeros(size(C));
for i = 1:numel(scales)for j = 1:size(C,2)if energy_mask(j,i) == 1 % 噪声区域denoised_C(j,i) = wthresh(C(j,i),'h',thr*0.8);else % 信号区域denoised_C(j,i) = wthresh(C(j,i),'s',thr*0.5);endend
end%% 逆变换重构
denoised_data = icwt(denoised_C, scales, wavelet, 'SamplingPeriod', dt);

3. 微震起始时间检测

%% 包络特征提取
envelope = hilbert(denoised_data);
env_amp = abs(envelope);%% 改进能量比函数（ER2）
ER2 = movmean(env_amp.^2, [5,5]) ./ movstd(data, [5,5]).^2;%% 初至时间检测（AIC准则）
[~, aic_peaks] = findpeaks(-ER2, 'MinPeakHeight', 0.1*max(ER2));
onset_time = aic_peaks(1)/fs; % 首个显著峰值对应时间

三、工程应用优化

1. 实时处理加速方案

% GPU并行加速
gpu_data = gpuArray(data);
[C_gpu, ~] = cwt(gpu_data, scales, wavelet, 'SamplingPeriod', dt);% 分块处理策略
block_size = 1024;
num_blocks = ceil(length(data)/block_size);
denoised_data = zeros(size(data));
parfor i = 1:num_blocksstart_idx = (i-1)*block_size +1;end_idx = min(i*block_size, length(data));block = data(start_idx:end_idx);denoised_block = process_block(block); % 自定义处理函数denoised_data(start_idx:end_idx) = denoised_block;
end

2. 多尺度融合策略

% 多分辨率融合算法
low_freq = denoise_band(data, 0.1, 10); % 低频段去噪
mid_freq = denoise_band(data, 10, 100); % 中频段去噪
high_freq = denoise_band(data, 100, 500); % 高频段去噪
denoised_data = reconstruct_band(low_freq, mid_freq, high_freq);