Inputs¶
Bayesian1DEM requires three types of inputs:
Data
Metadata
Inversion parameters
each of which is discussed below.
Data¶
Bayesian1DEM is set up to invert three types of EM data:
Magnetotelluric (MT)
Controlled-source electromagnetic (CSEM)
DC resistivity (DCR)
The data are stored in the main structure, S, as follows:
MT
MT apparent resistivity and phase are stored in the substructure S.MTdat, whose elements are
freqs — an array (1 x nFrequencies) of the frequencies at which the MT data has been processed
TEappRes –— an array (1 x nFrequencies) of the TE apparent resistivity
TEappResErr –— an array (1 x nFrequencies) of the uncertainty in the TE apparent resistivity
TEphase –— an array (1 x nFrequencies) of the TE phase
TEphaseErr –— an array (1 x nFrequencies) of the uncertainty in the TE phase
The forward modeling of MT apparent resistivity and phase is done in MT1D.m, which uses the standard Ward and Hohmann [WardHohmann87] recursive forumla.
CSEM
Two types of CSEM data is supported: surface-towed data like that collected by the Scripps Porpoise system [Sherman15], and generic dipole-dipole CSEM data. Both assume point or finite length dipole transmitters, and the forward response in both cases is computed using Dipole1D [Key09]. The difference is simply that the surface-towed system assumes a number of receivers for each transmitter location, while the generic system assumes a number of transmitters for each receiver. Since this is 1D modeling, it is the source-receiver offset that is relevant.
In either case, the inverted data are inline electric field amplitude and phase, which are stored in the substructure S.stDat or S.obDat (for “ocean-bottom”), whose elements are identical:
Freqs –— an array (nFrequencies x 1) of the frequencies at which the data has been collected/processed
Er –—— an array (nFrequencies x nReceivers) of the log10 of the amplitude of the in-line electric field
ErErr –—— an array (nFrequencies x nReceivers) of the error (uncertainty) in the log10 of the amplitude of the in-line electric field
Phase –—— an array (nFrequencies x nReceivers) of the phase of the in-line electric field
PhaseErr –—— an array (nFrequencies x nReceivers) of the phase of the error (uncertainty) in the in-line electric field
The in-line electric field amplitude at the second inversion frequency at the fourth receiver is then S.Dat.Er(2,4), and its uncertainty is S.Dat.ErErr(2,4). Thus, each row corresponds to a frequency and each column to a receiver.
DC resistivity
DC apparent resistivity data are for the Schlumberger system setup and are computed using . These data are stored in S.dcrDat:
es — an array (nElectrodeSpacings x 1) of the electrode spacings at which the data were collected
appRes — an array (nElectrodeSpacings x 1) of the DC apparent resistivity
appResErr — an array (nElectrodeSpacings x 1) of the error in the DC apparent resistivity
Metadata¶
The structure S contains an array indicating which data are being inverted — Bayesian1DEM natively supports joint inversion; simply indicate which data are included and the algorithm will jointly invert all of them (see [Blatter19]).
S.dataTypes — an array (4 x 1) of logicals indicating which data types are being inverted.
The positions refer to surface-towed CSEM, ocean bottom CSEM, MT, and DC resistivity, in that order. Thus, the array [0 0 1 0] indicates to Bayesian1DEM that only MT data are to be inverted.
Additionally, both CSEM data types require two metadata substructures—one for transmitter geometries and the other for receiver geometries.
S.stTx or S.obTx — a substructure containing all the necessary information about the CSEM transmitter geometry:
X—an array of size 1 x nTx where nTx is the number of transmitter locations; contains the x-coordinates (strike direction) of the transmitter locations
Y—an array of size 1 x nTx; contains the y-coordinates (in-line direction) of the transmitter locations
Z—an array of size 1 x nTx ; contains the z-coordinates (depth) of the transmitter locations
Dip—float; the dip of the transmitter from horizontal (degrees).
Azimuth—float; the angle of the transmitter from the strike direction (degrees). An azimuth of 90 degrees puts the transmitter in-line with the tow direction
Length—float; length (meters) of the transmitter. A length of 0 tells the forward solver to treat the transmitter as a point dipole.
WaterDepth—float; water depth (meters).
Soundings—an array of size 1 x nTx. The values are not important, only the length, which must be the same as the lengths of the X, Y, and Z arrays.
S.stRx or S.obRx — a substructure containing all the necessary information about the CSEM receiver geometry (Note: receiver rotations to the in-line/cross-line coordinate system must be done prior to inversion. This code is only setup to handle the in-line electric field component and assumes the receiver data is rotated into the in-line/cross-line coordinate system):
X—an array of size nRx x nTx where nRx is the number of surface-towed receivers and nTx is the number of transmitter locations; contains the x-coordinates (strike direction) of the receiver locations (in meters).
Y—an array of size nRx x nTx; contains the y-coordinates (in-line direction) of the receiver locations (in meters).
Z—an array of size nRx x nTx; contains the z-coordinates (depth) of the receiver locations (in meters).
Inversion parameters¶
Bayesian1DEM requires that the user set a handful of parameters prior to running an inversion. These are stored in the structures S, S1, and S2. Many of them are true or false flags used to toggle optional functionality. If any of these options remain murky or you are uncertain how to use them after reading this section, consult the examples. Each of the Step1 scripts in the Examples folder defines all the options necessary to invert those data. Copying them in the creation of your own Step1 script is a great way to start.
S
S.fixedDimension — a logical (set to “true” or “false”/1 or 0) that tells Bayesian1DEM not to sample over the number of layer interfaces or their locations. If you set this true, you must also provide the list of fixed interfaces, including S.regionBoundaryDepth, in the array S.TrueModel.z.
S.fixedWater — similar to S.fixedDimension, a logical that tells Bayesian1DEM not to sample over the water column (model region 1) resistivities. This assumes that S.fixedDimension is also set to true and that the fixed water column layers have been provided in S.TrueModel.z. It also requires that you provide a list of water column resistivities in S.TrueModel.rho.
S.debug_prior — a logical that causes Bayesian1DEM to ignore the likelihood (the forward codes are not called), such that the posterior equals the prior. This is for testing/debugging the prior.
S.transform01_ab — a logical that enables a depth-dependent prior distribution. For this feature to work, you must also provide the following:
S.maxRho — an array of maximum allowed resistivity (the upper bound of the depth-dependent prior)
S.minRho — an array of minimum allowed resistivity (the lower bound of the depth-dependent prior). Must be the same length as S.maxRho.
S.zRhoLim — an array (same length as S.maxRho and S.minRho) that specifies the depth (m) associated with each upper and lower bound resistivity. Bayesian1DEM linearly interpolates to all other depths when needed.
S.prior — an array (same length as each of the above) that specifies the value of the prior distribution at each depth: 1/(S.maxRho - S.minRho).
S.dz — bin width (m) used to generate an interpolation from the upper and lower resistivity bounds given above. If S.dz is really fine, the forward calculations will take longer to compute; if it is too coarse, the depth-dependent prior will likewise be coarsely defined.
S.zMin, S.zMax — the minimum and maximum model depths over which the depth-dependent prior is defined. Should cover the entire inversion domain.
S.logZ — “true” means Bayesian1DEM will sample the base 10 logarithm of depth; “false” means it will sample linear depth.
S.regionBoundaryDepth — the depth of the boundary between model region 1 and model region 2
S.nTemps — integer number of parallel tempering temperatures to use. More temperatures means better convergence properties but more chains in parallel and therefore more computational cost.
S.nChains — Since the set of temperatures (the “temperature ladder”) includes T = 1, nTemps - nChains + 1 is the number of chains at T = 1.
S.Tmax — the temperature ladder starts at T = 1 and ends at T = Tmax, with nTemps temperatures spaced logarithmically in between the two end members
S.B — the temperature ladder, which in the Step1 scripts in the examples auto-populates using nChains, nTemps, and Tmax.
S.numIterations — total number of iterations Bayesian1DEM will run before terminating
S.saveEvery — iteration interval between saves/writing model ensemble to disk. This can be a very time consuming step, so make it as long as you think reasonably safe. I typically keep it the same as numIterations, so that the save step occurs only once, at the end.
S.displayEvery — how often Bayesian1DEM will write progress output to the command window. I usually set this to 100 or so when first running an inversion to check that the RMS misfit is converging, then start the inversion over again with a higher value
S.beta – the fraction of the time that Bayesian1DEM will spend sampling model region 1. By implication (1-beta) is the fraction of the time Bayesian1DEM will spend sampling model region 2.
S.TrueModel — if fixedDimension is set to “true”, this holds the depths of the interfaces. In the case of synthetic inversion, this holds (in S.TrueModel.z and S.TrueModel.rho) the “true model”.
S1
S1.isotropic — must be set to 1 or “true”. Currently, Bayesian1DEM only inverts for isotropic resistivity models.
S1.birth_death_from_prior — set to 1 to propose new layer resistivity during a “birth” step (adding a new layer) from the prior distribution. Set to 0 to propose the new layer resistivity from a Gaussian whose mean is the current layer resistivity of the layer the new interface splits.
S1.kInit — the number of layers in the initial model regions at the start of the algorithm. Can’t be less than 1, and should probably be set to 1 unless you have a good reason to do otherwise.
S1.kMin, S1.kMax — the minimum and maximum allowed number of interfaces in model region 1.
S1.zMin, S1.zMax — the minimum and maximum depths at which interfaces can be placed by Bayesian1DEM
S1.rhoMin, S1.rhoMax — the minimum and maximum allowed resistivity of layers in model region 1. These define a uniform prior distribution for layer resistivity.
S1.UstepSize – the variance of the Gaussian distribution used when perturbing the resistivities of the model layers (during an update step). There should be S.nChains of these variances, increasing with temperature.
S1.BstepSize – the variance of the Gaussian distribution used when determining the resistivity of a newly created layer (during a “birth” step), unless S1.birth_death_from_prior is set to 1, in which case this array is ignored. Should be an array S.nChains elements long (one for each parallel tempering chain, as above)
S1.BstepSize – the variance of the Gaussian distribution used when determining the new location of an interface during a move step. As usual, there should be S.nChains of these in this array.
S2
The input parameters for S2 are identical to those for S1, but apply to model region 2 instead.