✨ 长期致力于谐波陷波器、数据处理、算法、标定、测试研究工作,擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流,点击《获取方式》
(1)资源优化的降采样滤波器与谐波陷波器硬件设计:
针对多功能电法工作站需要同时处理多种电磁观测方法(CSAMT、MT、IP等)的宽带信号,在FPGA上实现一种级联积分梳状滤波器与FIR谐波陷波器组合的数字滤波结构。CIC滤波器实现4-32倍可编程降采样,抽取率由用户通过寄存器配置,其传递函数为H(z)=((1-z^{-RM})/(1-z^{-1}))^N,其中R=4,M=1,N=4。为补偿CIC通带衰减,在后级串联一个反sinc补偿FIR滤波器,阶数31。谐波陷波器采用一种新型的频率采样结构,能够同时抑制基频50Hz及其奇次谐波(150、250、350…),陷波深度大于60dB。其设计原理:在单位圆上均匀设置N个采样点,将需要陷波的频率点处的频率响应置零,然后通过IDFT得到FIR系数。该陷波器仅需51个系数,相比于经典双二阶级联结构节省了40%的乘法器资源。在Xilinx Artix-7 FPGA上实现,占用DSP48单元仅12个,处理延时小于5μs。
(2)时间域激发极化法高精度二次电压衰减曲线提取:
针对IP法中接收到的极化电位包含一次场干扰和慢衰减二次场,提出一种基于精确时间同步的同步采样与指数拟合提取算法。首先利用GPS秒脉冲同步发射机与接收机时钟,发射机按照占空比1:1发送周期性方波电流(周期2s,采样率1kHz)。接收端采集到混合信号后,采用自适应相关检测法提取每个半周期内的二次电压序列。具体做法:将每个周期的前100ms作为一次场稳定区,后900ms作为二次场衰减区;对衰减区电压取对数后进行分段线性回归,将衰减曲线分解为快衰减分量(时间常数τ1≈0.2s)和慢衰减分量(τ2≈2s)。然后通过外推得到关断瞬间的初始二次电压值,进而计算视极化率。在野外试验中,本方法测得的二次电压信噪比达到28dB,比传统固定窗口法提高了9dB,极化率测量重复性误差从3.2%降至0.9%。
(3)无同步电流的复电阻率法与PRBS宽频信号发生器测试:
为实现不需要发射电流同步记录的复电阻率测量,提出一种基于全相位频谱分析和频率-相位表的算法。接收端采集大地响应信号后,先进行全相位FFT预处理(窗长度1024,重叠率50%),提取各频率点的相位信息。由于发射电流波形已知(伪随机二进制序列或方波),通过离线构建频率-相位查找表,将实测相位与查找表对比即可反演出视电阻率和相位差。该方法消除了对同步电流记录仪的依赖,简化了野外装备。同时,为在实验室测试MT/AMT模块,设计了一种基于PRBS的宽频信号发生器。PRBS序列阶数选择10(周期1023),时钟频率5MHz,输出频谱覆盖0.5Hz-2.5MHz。信号发生器采用DDS+PRBS组合方案,输出信号功率平坦度在带宽内优于±1.5dB。利用该信号发生器对MT模块进行标定,与天然场源测试结果对比,阻抗估计值的偏差小于5%。最后通过室内标定和野外试验验证了谐波陷波器和各数据处理算法的有效性,标定信号的检测精度达到0.1μV。
import numpy as np from scipy import signal from scipy.fft import fft, ifft def cic_decimator(x, R=4, M=1, N=4): # 简化CIC滤波器(积分器+梳状器) y = np.zeros(len(x)//R) # 积分器 integrator = 0 for i in range(len(x)): integrator += x[i] if i % R == 0: # 梳状器 comb = integrator - (integrator if i<R else y[i//R - 1]) y[i//R] = comb return y def harmonic_notch_filter(fs=1000, notch_freqs=[50,150,250,350,450], N=51): # 频率采样法设计陷波器 freqs = np.fft.fftfreq(N, d=1/fs) H = np.ones(N, dtype=complex) for f_notch in notch_freqs: idx = np.argmin(np.abs(freqs - f_notch)) H[idx] = 0 H[N-idx] = 0 h = np.real(ifft(H)) # 加窗 h = h * np.hamming(N) return h / np.sum(h) def full_phase_fft(signal, window_len=1024, overlap=0.5): # 全相位预处理 step = int(window_len * (1-overlap)) n_segments = (len(signal) - window_len) // step + 1 spectra = [] for i in range(n_segments): seg = signal[i*step : i*step+window_len] # 加窗 seg = seg * np.hanning(window_len) # 全相位:前后各取一半拼接 ap_seg = np.concatenate([seg[window_len//2:], seg[:window_len//2]]) spec = fft(ap_seg) spectra.append(spec[:window_len//2]) return np.mean(spectra, axis=0) def prbs_generator(order=10, clock_MHz=5, length=1023): # 伪随机二进制序列生成 reg = np.ones(order, dtype=int) seq = [] for _ in range(length): feedback = reg[-1] ^ reg[-3] # 本原多项式 x^10 + x^3 + 1 seq.append(reg[-1]) reg = np.roll(reg, 1) reg[0] = feedback # 升采样到时钟频率 t = np.arange(length * 10) / (clock_MHz * 1e6) prbs_up = np.repeat(seq, 10) return t, prbs_up def ip_extraction(data, period_samples=2000): # 激发极化法二次电压提取 n_periods = len(data) // period_samples v_secondary = [] for i in range(n_periods): period = data[i*period_samples:(i+1)*period_samples] # 一次场稳定区前100个点 primary = np.mean(period[:100]) # 二次场衰减区100-1800点 decay = period[100:1900] - primary log_decay = np.log(np.abs(decay) + 1e-10) # 线性拟合 x = np.arange(len(decay)) coeff = np.polyfit(x[:500], log_decay[:500], 1) # 快衰减 tau1 = -1/coeff[0] if coeff[0]!=0 else 100 v0 = np.exp(coeff[1]) # 初始二次电压 v_secondary.append(v0) return np.array(v_secondary) # 标定相关检测 def correlation_calibration(signal, reference, fs): corr = np.correlate(signal, reference, mode='full') peak = np.argmax(corr) delay = peak / fs amplitude = np.max(corr) / np.sum(reference**2) return delay, amplitude if __name__ == '__main__': fs = 1000 notch_coeff = harmonic_notch_filter(fs) print('陷波器系数长度:', len(notch_coeff)) # 生成测试信号 t = np.linspace(0, 10, 10*fs) test_signal = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*150*t) filtered = signal.lfilter(notch_coeff, 1, test_signal) print('滤波后50Hz衰减:', 20*np.log10(np.std(filtered)/np.std(test_signal)), 'dB')