Как я могу правильно восстановить сигнал, используя signal.deconvolve, примененный к сигналу, свернутому с помощью вейвлета Рикера, используя signal.convolve?

Итак, у меня есть массив, хранящийся в матрице с размерами (251, 240). Затем я создал вейвлет Рикера, который я сворачиваю с каждым столбцом (временным рядом). Кажется, это работает нормально. Следующим шагом в моем процессе будет деконволюция результата свертки с тем же вейвлетом Рикера. Я бы ожидал, что мой исходный сигнал будет реконструирован, однако это не так. Что я делаю неправильно и как правильно развернуть вейвлет Рикера?

Я прикрепляю часть своего кода ниже

# the array 'time' and and 'seismic_data' with dimensions (251,) and (251,240) respectively, where created in previous cells.
from scipy import signal
ricker2 = signal.ricker(time.size, 4) #create ricker wavelet with same dimensions as time series
seismogram = []
seismic_trace = [signal.convolve(ricker2,seismic_data[:,i], mode='same')for i in range(offsets.size - 1)] #creates array with each row being the time series convolved
seismogram = np.array(seismic_trace).T #transpose to get time series as columns 
# PlottingI 
fig,ax = plt.subplots(figsize=(8,10))
ax.imshow(seismic_data,  aspect=1400,  extent= [np.min(offsets), np.max(offsets),np.max(time),np.min(time)]) 
plt.title('Central Shot Seismic Gather \n ', fontweight="bold")
plt.xlabel('Offset [m]', fontweight="bold")
plt.ylabel('Time [s]', fontweight="bold")
plt.show()

График, показанный здесь, - это то, что я ожидал от свертки. (Я пока не могу добавлять изображения. Извините за это).

Следующий шаг (деконволюция) — это то, что не работает нормально.

deconvolved_seismogram = []
for i in range (offsets.size-1):
    rems, deconvolved_trace = signal.deconvolve(seismogram[:,i],ricker2)
    deconvolved_seismogram.append(deconvolved_trace)
deconvolved_seismogram = np.array(deconvolved_seismogram).T

fig,ax = plt.subplots(figsize=(8,10))
ax.imshow(deconvolved_seismogram,  aspect=1400,  extent= [np.min(offsets), np.max(offsets),np.max(time),np.min(time)]) 
plt.title('Deconvolved Central Shot Seismic Gather \n ', fontweight="bold")
plt.xlabel('Offset [m]', fontweight="bold")
plt.ylabel('Time [s]', fontweight="bold")
plt.show()

Массив deconvolved_seismogram имеет правильные размеры, но сигнал совсем не похож на исходный сигнал (seismic_data). Что я делаю не так? Как я могу это исправить?

Заранее спасибо!


person Fer Berumen    schedule 19.10.2020    source источник


Ответы (1)


Основная проблема заключается в том, что свертка с помощью фильтра может удалить информацию. Если частотная характеристика фильтра содержит нули (или числовые точки, близкие к нулю), то эти частотные компоненты невосстановимы, и полная деконволюция невозможна.

Другая проблема заключается в том, что scipy.signal.deconvolve пытается выполнить точную деконволюцию (путем полиномиального деления); это не шумоустойчивый метод. Деконволюция, по сути, представляет собой точечное деление в частотной области на частотную характеристику фильтра. Если отклик где-то мал, это деление усилит любой шум, включая числовую ошибку округления — я полагаю, это объясняет ваше наблюдение, что результат деконволюции с scipy.signal.deconvolve совсем не похож на исходный сигнал.

Вот график rickert2, а внизу график его частотной характеристики.

Сюжет ricker2

ricker2 очень плавный! Это определенно трудно деконволюционировать. Нижний график показывает, что частотная характеристика имеет ноль на постоянном токе и быстро падает (обратите внимание, что ось Y в дБ) выше 0,2 цикла/выборка или около того. Это означает, что некоторое среднечастотное содержимое может быть восстановлено с помощью помехоустойчивого метода деконволюции, но постоянное и более высокочастотное содержание исчезает.

Я предлагаю попробовать деконволюцию Винера. Это классический метод надежной деконволюции. В этот гитхаб-список есть реализация Python (см. функцию wiener_deconvolution).

person Pascal Getreuer    schedule 21.10.2020