Executing Pylops with SEGY data containing a variable number of traces per shot

[guaacoelho]

Hello,

Peterson Nogueira and I use Pylops for geophysical applications. For this, we primarily use the AcousticWave2D operator.

Currently, we are working with the goal of running these applications using data stored in SEGY files. However, one factor is adding complexity to this process: data from SEGY files usually has a variable number of receivers for each modeled seismogram/shot. The AcousticWave2D has well-defined input (dims) and output (dimsd) data dimensions, where it is expected that the number of receivers be the same for each shot. Thus, returning the forward modeling data becomes a problem, as the resulting data matrix would not be regular.

I would like to know if the tool already has any feature that could help us overcome this issue.

[mrava87]

Hi @guaacoelho, thanks for raising this issue.

You are correct about the fact that the AcousticWave2D currently supports only a fixed number of sources and receivers. In theory we could easily allow for them to change internally (a user would provide a list of arrays for src_x, src_z, rec_x, rec_z where each array corresponds to a specific shot) but I am not sure how sustainable is this.

A possibly better option which you could test already is to do something like this:

# Variable geometry
nsrc = 3
nrec = 200
src_x = np.linspace(0, 500, num=nsrc)
src_z = [100., 200., 300.]

rec_x_rel = np.linspace(0, 3000, num=nrec)
rec_z = np.array(src_z)[:, None] + np.random.normal(0, .5, (nsrc, nrec))

Aops = []
for isrc in range(nsrc): 
    Aop = AcousticWave2D(shape, origin, spacing, v0 * 1e3,
                         src_x[isrc:isrc+1], src_z[isrc:isrc+1], 
                         rec_x_rel + src_x[isrc], rec_z[isrc], 
                         0., tn, 'Ricker',
                         space_order=6, nbl=100, f0=15,
                         dtype="float32", name="A")
    Aops.append(Aop)
Aops = VStack(Aops)

I think Devito is smart internally and does only point to numpy objects (eg velocity model), it does not copy them... so each operator Aop will be rather lightweight and can store its own geometry 😄