Pre-stack modeling operator for multiple wavelets

[dprjoao]

🛠 Multi-Wavelet Support in Pre-Stack Inversion Operator (explicit=False)

📖 Description

The convolution operator used by pylops for pre-stack inversion (explicit=False) is implemented using the Conv1D class. Currently, pylops supports only a single wavelet for all angles, meaning the inversion assumes the same wavelet for each angle.

I attempted to modify the code to allow a different wavelet for each angle (mid-angle in my case, as I’m using partial stacks). However, my current results are only consistent when using dense operators (explicit=True), whereas I would expect them to also work in the sparse case.

---

🔍 Observations

• ✅ The dense implementation (explicit=True) behaves as expected.
• ❌ The sparse implementation (explicit=False) does not produce comparable results.
• ⚠️ The issue seems related to how the Conv1D operator is built when handling multiple wavelets.

The code using Conv1D for multiple wavelets is only working for spatdims=1 at the moment, the idea is to extend for 2D or 3D cases when i adress this issue.

---

💡 Code Snippet

✅ Dense Implementation (Working)

# Handling multiple wavelets
if isinstance(wav, list):
    print("Dense multi-wavelet operator")
    n = len(wav)
    C = np.zeros((n * nt0, n * nt0))

    wsub_list = []
    for w in wav:
        wsub = np.asarray(convmtx(w, nt0, len(w) // 2)[:nt0])
        wsub_list.append(wsub)

    C = get_block_diag(theta)(*wsub_list)

else:
    # Create wavelet operator
    C = ncp.asarray(convmtx(wav, nt0, len(wav) // 2)[:nt0])
    C = [C] * ntheta
    C = get_block_diag(theta)(*C)

# Combine operators
M = ncp.dot(C, ncp.dot(G, D))
Preop = MatrixMult(M, otherdims=spatdims, dtype=dtype)

❌ Sparse implementation (Not-working)

    else:
        if isinstance(wav, list):
            print("Dense multi-wavelet operator")
            # Criar uma lista de operadores de convolução
            Cop_list = []

            for w in wav:
                Cop = Convolve1D(
                    (nt0),
                    h=w,
                    offset=len(w) // 2,  # Define o centro da wavelet
                    axis=0,
                    dtype='float32',
                    method='fft'
                )
                Cop_list.append(Cop.tosparse())  # Adiciona à lista de operadores'
            

            Cop = BlockDiag(Cop_list)

        else:
            print("Sparse single wavelet operator")
            # Create wavelet operator
            Cop = Convolve1D(
                dims,
                h=wav,
                offset=len(wav) // 2,
                axis=0,
                dtype=dtype,
            )

[mrava87]

@dprjoao thanks for attempting this. I think I know what is the problem, dims=(nt0, ntheta, ...) so when you do Cop = BlockDiag(Cop_list) you apply the different wavelets to different parts of the first axis, so time - this is how BlockDiag works as it assumes the input is a flattened array.

Other than re-write the entire operator to have ntheta in the first dimension (which is not an option as it will break backward compatibility), I think the easiest and most elegant solution for this is to add a Transpose operator before the Cop, so something along these lines (I did not test the code):

else:
        if isinstance(wav, list):
            print("Dense multi-wavelet operator")
            # Criar uma lista de operadores de convolução
            Cop_list = []

            for w in wav:
                Cop = Convolve1D(
                    (nt0),
                    h=w,
                    offset=len(w) // 2,  # Define o centro da wavelet
                    axis=0,
                    dtype='float32',
                    method='fft'
                )
                Cop_list.append(Cop.tosparse())  # Adiciona à lista de operadores'
            
            Top = Transpose(dims, axes=(1, 0, *(range(2, len(spatdims) + 2)))
            Cop = BlockDiag(Cop_list)
            ...
            Preop = Top.H * Cop * Top * AVOop * Dop

the Top.H is needed so the output shapes are consistent with those of the single wavelet case.

Let me know how it goes 😄

[dprjoao]

@mrava87 i first tested this solution and seems to work, i will look deeper and i will let you know for sure.

Thanks a lot!