基于雷克子波和反射系数的地震记录生成算法

import numpy as np import matplotlib.pyplot as plt def ricker_wavelet(frequency, length, dt): num_samples = int(length / dt) # 计算样本数量 t = np.linspace(-length / 2, length / 2, num_samples) wavelet = (1 - 2 * (np.pi ** 2) * (frequency ** 2) * (t ** 2)) * np.exp(-(np.pi ** 2) * (frequency ** 2) * (t ** 2)) time_axis = np.arange(0, num_samples) * dt return wavelet, time_axis def convolution(reflection_coefficients, wavelet): num_samples = len(reflection_coefficients) num_wavelet_samples = len(wavelet) output_length = num_samples + num_wavelet_samples - 1 convolution_result = np.convolve(reflection_coefficients, wavelet, mode='full') output_time_axis = np.arange(0, output_length) return convolution_result, output_time_axis # 反射系数 reflection_coefficients = [0.3, -0.5, 0.7, -0.2, 0.1] # 雷克子波参数 frequency = 25 # 频率 length = 0.128 # 波形长度（秒） dt = 0.001 # 采样间隔（秒） # 生成雷克子波 wavelet, time_axis = ricker_wavelet(frequency, length, dt) # 进行褶积操作 convolution_result, output_time_axis = convolution(reflection_coefficients, wavelet) # 绘制原始波形和褶积结果 plt.figure(figsize=(10, 6)) plt.subplot(2, 1, 1) plt.plot(time_axis, wavelet) plt.title('Ricker Wavelet') plt.xlabel('Time (seconds)') plt.ylabel('Amplitude') plt.subplot(2, 1, 2) plt.plot(output_time_axis, convolution_result) plt.title('Convolution Result') plt.xlabel('Time (samples)') plt.ylabel('Amplitude') plt.tight_layout() plt.show() #