How can I properly reconstruct a signal using signal.deconvolve applied to a signal convolved with a ricker wavelet using signal.convolve?

Question

So, I have an array stored in a matrix with dimensions (251, 240). Next I created a ricker wavelet which I convolve with each column (time series). This seems to work fine. The next step in my process would be to deconvolve the outcome of the convolution with the same ricker wavelet. I would expect to get my original signal reconstructed, however this is not the case. What am I doing wrong and how can I deconvolve the ricker wavelet properly?

I'm attaching some of my code below

# 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()

The plot displayed here is what I would expect from the convolution. (I can't add images yet. Sorry for that).

The next step (deconvolution) is what doesn't seem to be working fine.

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()

The deconvolved_seismogram array has the right dimensions, but the signal is not at all similar to the original signal (seismic_data). What am I doing wrong? How can I correct this?

Thanks in advance!

Answer 1

The underlying problem is that convolution by a filter may remove information. If the filter's frequency response contains zeros (or points numerically near zero), then those frequency components are unrecoverable and a complete deconvolution is impossible.

The other problem is that scipy.signal.deconvolve attempts to deconvolve exactly (by polynomial division); it is not a noise robust method. Deconvolution is essentially pointwise division in the frequency domain by the filter's frequency response. If the response is small anywhere, this division will amplify any noise, including numerical round-off error—I believe this explains your observation that the deconvolved result with scipy.signal.deconvolve is "not at all similar to the original signal".

Here is a plot of rickert2, and on the bottom, a plot of its frequency response.

ricker2 is extremely smooth! This is definitely difficult to deconvolve. The bottom plot shows that the frequency response has a zero at DC and rapidly falls off (note the y-axis is in dB) above 0.2 cycles/sample or so. This means some mid-range frequency content could be recovered by a noise-robust deconvolution method, but DC and higher frequency content are gone.

I suggest to try Wiener deconvolution. This is a classic method for robust deconvolution. There is a Python implementation in this github gist (see the "wiener_deconvolution" function).