fmi_inversion: Invert fluvial data set based on reference spectra catalogue

In eseis: Environmental Seismology Toolbox

Invert fluvial data set based on reference spectra catalogue

Description

The fluvial model inversion (FMI) routine uses a predefined look-up table with reference spectra to identify those spectra that fit the empirical data set best, and returns the corresponding target variables and fit errors.

Usage

fmi_inversion(reference, data, n_cores = 1)

Arguments

reference	List containing lists with precalculated model spectra.
data	eseis object or numeric matrix (spectra organised by columns), empiric spectra which are used to identify the best matching target parameters of the reference data set.
n_cores	Numeric value, number of CPU cores to use. Disabled by setting to 1. Default is 1.

Details

Note that the frequencies of the empiric and modelled data sets must match.

Value

List object containing the inversion results.

Author(s)

Michael Dietze

Examples

## NOTE THAT THE EXAMPLE IS OF BAD QUALITY BECAUSE ONLY 10 REFERENCE 
## PARAMETER SETS AND SPECTRA ARE CALCULATED, DUE TO COMPUATATION TIME
## CONSTRAINTS. SET THE VALUE TO 1000 FOR MORE MEANINGFUL RESULTS.

## create 100 example reference parameter sets
ref_pars <- fmi_parameters(n = 10,
                           h_w = c(0.02, 1.20),
                           q_s = c(0.001, 8.000) / 2650,
                           d_s = 0.01,
                           s_s = 1.35,
                           r_s = 2650,
                           w_w = 6,
                           a_w = 0.0075,
                           f_min = 5,
                           f_max = 80,
                           r_0 = 6,
                           f_0 = 1,
                           q_0 = 10,
                           v_0 = 350,
                           p_0 = 0.55,
                           e_0 = 0.09,
                           n_0_a = 0.6,
                           n_0_b = 0.8,
                           res = 100)

## create corresponding reference spectra
ref_spectra <- fmi_spectra(parameters = ref_pars)

## define water level and bedload flux time series
h <- c(0.01, 1.00, 0.84, 0.60, 0.43, 0.32, 0.24, 0.18, 0.14, 0.11)
q <- c(0.05, 5.00, 4.18, 3.01, 2.16, 1.58, 1.18, 0.89, 0.69, 0.54) / 2650
hq <- as.list(as.data.frame(rbind(h, q)))

## calculate synthetic spectrogram
psd <- do.call(cbind, lapply(hq, function(hq) {

  psd_turbulence <- eseis::model_turbulence(h_w = hq[1],
                                            d_s = 0.01,
                                            s_s = 1.35,
                                            r_s = 2650,
                                            w_w = 6,
                                            a_w = 0.0075,
                                            f = c(10, 70),
                                            r_0 = 5.5,
                                            f_0 = 1,
                                            q_0 = 18,
                                            v_0 = 450,
                                            p_0 = 0.34,
                                            e_0 = 0.0,
                                            n_0 = c(0.5, 0.8),
                                            res = 100, 
                                            eseis = FALSE)$power

  psd_bedload <- eseis::model_bedload(h_w = hq[1],
                                      q_s = hq[2],
                                      d_s = 0.01,
                                      s_s = 1.35,
                                      r_s = 2650,
                                      w_w = 6,
                                      a_w = 0.0075,
                                      f = c(10, 70),
                                      r_0 = 5.5,
                                      f_0 = 1,
                                      q_0 = 18,
                                      v_0 = 450,
                                      x_0 = 0.34,
                                      e_0 = 0.0,
                                      n_0 = 0.5,
                                      res = 100,
                                      eseis = FALSE)$power

  ## combine spectra
  psd_sum <- psd_turbulence + psd_bedload

  ## return output
  return(10 * log10(psd_sum))
}))

graphics::image(t(psd))

## invert empiric data set
X <- fmi_inversion(reference = ref_spectra, 
                   data = psd)

## plot model results
plot(X$parameters$q_s * 2650, 
     type = "l")
plot(X$parameters$h_w, 
     type = "l")

eseis documentation built on Aug. 10, 2023, 5:08 p.m.