PyLops / pylops

PyLops – A Linear-Operator Library for Python
https://pylops.readthedocs.io
GNU Lesser General Public License v3.0
420 stars 101 forks source link

It seems to be the two-dimensional hyperbolic Radon transformation is not invertible? #560

Closed ZhangWeiGeo closed 8 months ago

ZhangWeiGeo commented 8 months ago
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 27 16:55:37 2023

@author: navie test
"""

import numpy as np
import scipy.io as scio

import pylops

import matplotlib.pyplot as plt

seismic = scio.loadmat('Seismic.mat')
seismic = seismic['Seismic'].T

dt=0.001 #(s)
dx=10 #(m)

taxis = np.arange(0, 1000)*dt
xaxis = np.arange(0, 200)*dx

npx =200
pxmax = 5e-5# s/m
px = np.linspace(-0, pxmax, npx)

Rop = pylops.signalprocessing.Radon2D(
    taxis, xaxis, px, kind="hyperbolic", interp="nearest", centeredh=True, dtype="float64"
)

Rseismic = Rop.H*seismic

# sparse Radon transform
xinv, niter, cost = pylops.optimization.sparsity.fista(
    Rop, seismic.ravel(), niter=15, eps=1e1
)
xinv = xinv.reshape(Rop.dims)

# Inverse transformation
yfilt = Rop * xinv

plt.subplot(141)
plt.imshow(seismic.T,cmap="seismic")
plt.subplot(142)
plt.imshow(Rseismic.T,cmap="seismic")
plt.subplot(143)
plt.imshow(xinv.T,cmap="seismic")
plt.subplot(144)
pclip=0.7
plt.imshow(yfilt.T,cmap="seismic",vmin=-pclip * np.abs(yfilt).max(),
    vmax=pclip * np.abs(yfilt).max(),
)

Why is this transformation so poorly reversible?

The data is Seismic.zip

Seismic.zip

mrava87 commented 8 months ago

Hi @ZhangWeiGeo, I’ll take a look into this but I guess it’s because of the parameter choice that is not suitable to the data. In other words, you cannot find a sparse representation that fully explains the data: you could relax eps in Fista or play with pxmax to ensure it spans the entire range of hyperbolas in the data :)

mrava87 commented 8 months ago

Ok, I checked the data... you will NEVER be able to find a sparse representation for this data with a hyperbolic Radon transform since the events are hyperbolas whose apex have different x-location. You will need to use something like this https://www.researchgate.net/profile/Daniel-Trad/publication/229003858_Apex_Shifted_Radon_transform/links/02e7e52323558d24a7000000/Apex-Shifted-Radon-transform.pdf (which is currently not available in PyLops)

mrava87 commented 8 months ago

@ZhangWeiGeo shall I close or you need any further clarification?

ZhangWeiGeo commented 8 months ago

@ZhangWeiGeo shall I close or you need any further clarification?

Thank for your reply. I don‘t have any further clarification.