Edit on GitHub

Seismic Preconditioners

JUDI provides a selected number of preconditioners known to be beneficial to FWI and RTM. We welcome additional preconditioners from the community. Additionnaly, any JOLI operator can be used as a preconditiner in conbination with JUDI operator thanks to the fundamental interface between JUDI and JOLI.

Model domain preconditioners

Model space preconditioners acts on model size arrays such as the velocity model or the FWI/RTM gradient. These preconditioners are indepenedent of the number of sources and therefore should not be indexed.

Water column muting (top mute)

Create a linear operator for a 2D model topmute, i.e. for muting the water column:

JUDI.TopMute — Type
TopMute{T, N, Nw}

Mute top of the model in N dimensions

Constructor

judiTopmute(model; taperwidht=10)
judiTopmute(n, wb, taperwidth)   # Legacy
source

Usage:

# Forward
m_mute = Mr*vec(m)

# Adjoint
m_mute = Mr'*vec(m)

As Mr is self adjoint, Mr is equal to Mr'.

legacy:

The legacy constructor judiTopmute(n, wb, taperwidth) is still available to construct a muting operator with user specified muting depth.

Model depth scaling

Create a 2D model depth scaling. This preconditionenr is the most simple form of inverse Hessain approximation compensating for illumination in the subsurface. We also describe below a more accurate diagonal approximation of the Hessian with the illlumination operator. Additionnaly, as a simple diagonal approximation, this operator is invertible and can be inverted with the standard julia inv function.

JUDI.DepthScaling — Type
DepthScaling{T, N, K}

Depth scaling operator in N dimensions scaling by depth^K.

Constructor

judiDepthScaling(model::AbstractModel; K=.5)
source

Usage:

# Forward
m_mute = Mr*vec(m)

# Adjoint
m_mute = Mr'*vec(m)

Illumination

The illumination computed the energy of the wavefield along time for each grid point. This provides a first order diagonal approximation of the Hessian of FWI/LSRTM helping the ocnvergence and quality of an update.

JUDI.judiIllumination — Type
judiIllumination(model; mode="u", k=1, recompute=true)

Arguments

  • model: JUDI Model structure

  • mode: Type of ilumination, choicees of ("u", "v", "uv")

  • k: Power of the illumination, real number

  • recompute: Flag whether to recompute the illumination at each new propagation (Defaults to true)

    judiIllumination(F; mode="u", k=1, recompute=true)

Arguments

  • F: JUDI propagator
  • mode: Type of ilumination, choicees of ("u", "v", "uv")
  • k: Power of the illumination, real positive number
  • recompute: Flag whether to recompute the illumination at each new propagation (Defaults to true)

Diagonal approximation of the FWI Hessian as the energy of the wavefield. The diagonal contains the sum over time of the wavefield chosen as mode.

Options for the mode are "u" for the forward wavefield illumination, "v" for the adjoint wavefield illumination, and "uv" for the pointwise product of the forward and adjoint wavefields illuminations. Additionally, the parameter "k" provides control on the scaling of the daiagonal raising it to the power k.

Example

I = judiIllumination(model)

Construct the diagonal operator such that I*x = x ./ |||u||_2^2

source

Usage:

# Forward
m_mute = I*vec(m)

# Adjoint
m_mute = I'*vec(m)

Data preconditioners

These preconditioners are design to act on the shot records (data). These preconditioners are indexable by source number so that working with a subset of shot is trivial to implement. Additionally, all DataPreconditionner are compatible with out-of-core JUDI objects such as judiVector{SeisCon} so that the preconditioner is only applied to single shot data at propagation time.

Data topmute

Create a data topmute for the data based on the source and receiver geometry (i.e based on the offsets between each souurce-receiver pair). THis operator allows two modes, :reflection for the standard "top-mute" direct wave muting and :turning for its opposite muting the reflection to compute gradients purely based on the turning waves. The muting operators uses a cosine taper at the mask limit to allow for a smooth transition.

JUDI.DataMute — Type
struct DataMute{T, mode} <: DataPreconditioner{T, T}
    srcGeom::Geometry
    recGeom::Geometry
    vp::Vector{T}
    t0::Vector{T}
    taperwidth::Vector{Int64}
end

Data mute linear operator a where {T, N}sociated with source srcGeom and receiver recGeom geometries used to compute the distance to the source for each trace in the data. Data mute preconditionner. Supports two modes (:reflection, :turning) that mute either the turning waves (standard direct wave mute) or mutes the reflections. A cosine tapr is applied with width taperwidth to avoid abrupt change and infinite frequency jumps in the data.

Constructors

judiDataMute(srcGeom, recGeom; vp=1500, t0=.1, mode=:reflection, taperwidth=floor(Int, 2/t0))

Construct the data mute operator from the source srcGeom and receiver recGeom geometries.

judiDataMute(q, d; vp=1500, t0=.1, mode=:reflection, taperwidth=floor(Int, 2/t0))

Construct the data mute operator from the judivector source q and judivector data d.

Parameters

The following optional paramet where {T, N}rs control the muting operator

  • vp: P wave velocity of the direct wave (usually water velocity). Can be a constant or a Vector with one value per source position. Devfaults to 1500m/s
  • t0: Time shift in seconds (usually width of the wavelet). Defaults to $.1 sec$
  • mode: :reflection to keep the reflections and mute above the direct wave (i.e for RTM) :turning to keep the turning waves and mute below the direct wave (i.e for FWI)
  • taperwidth: Width of the cosine taper in number of samples. Defaults to 2 / t0
source

Band pass filter

While not purely a preconditioner, because this operator acts on the data and is traditionally used fro frequency continuation in FWI, we implemented this operator as a source indexable linear operator as well. Additionally, the filtering function is available as a standalone julia function for general usage

JUDI.FrequencyFilter — Type
struct FrequencyFilter
    recGeom

Bandpass filter linear operator. Filters the input judiVector or Vector

Constructor

judiFilter(geometry, fmin, fmax) 
judiFilter(judiVector, fmin, fmax)
source
JUDI.filter_data — Function
filter(Din, dt_in; fmin=0, fmax=25)

Performs a causal filtering [fmin, fmax] on the input data bases on its sampling rate dt. Automatically perfroms a lowpass if fmin=0 (default)

source

Data time derivative/intergraton

A TimeDifferential{K} is a linear operator that implements a time derivative (K>0) or time integration (K<0) of order K for any real K including fractional values.

JUDI.TimeDifferential — Type
struct TimeDifferential
    recGeom

Differential operator of order K to be applied along the time dimension. Applies the ilter w^k where k is the order. For example, the tinme derivative is TimeDifferential{1} and the time integration is TimeDifferential{-1}

Constructor

judiTimeIntegration(recGeom, order)
judiTimeIntegration(judiVector, order)

judiTimeDerivative(recGeom, order)
judiTimeDerivative(judiVector, order)
source

Inversion wrappers

For large scale and practical cases, the inversions wrappers fwi_objective and lsrtm_objective are used to minimize the number of PDE solves. Those wrapper support the use of preconditioner as well for better results.

Usage:

For fwi, you can use the data_precon keyword argument to be applied to the residual (the preconditioner is applied to both the field and synthetic data to ensure better misfit):

fwi_objective(model, q, dobs; data_precon=precon)

where precon can be:

Similarly, for LSRTM, you can use the model_precon keyword argument to be applied to the perturbation dm and the data_precon keyword argument to be applied to the residual:

lsrtm_objective(model, q, dobs, dm; model_precon=dPrec, data_precon=dmPrec)

where dPrec and dmPrec can be: