Variable discretization

There are a few cases when we can simplify the subsurface structure by expressing it using only a few layers. In those cases, it may be desirable to vary the cell size (layer thickness) along with the model parameters. This allows us to move the layer interfaces up and down while parameterizing the model space in a different way wherein there is more flexibility.

Let's denote the model parameters, eg., conductivity, by m, and the layer thickness by h. Therefore, in a N-layer case, we will have

$$\begin{gathered} m=[m_1,m_2,m_3,...,m_N] \\ h=[h_1,h_2,h_3,...,h_{N_1}] \end{gathered}$$

such that

$$\begin{gathered} m_i\in {\mathcal{D}}_{m_i}\text{ ; where }{\mathcal{D}}_{m_i}=\mathit{\text{a priori}}\text{ distribution for }{m}_i\text{ ; and} \\ {h}_i\in {\mathcal{D}}_{h_i}\text{ ; where }{\mathcal{D}}_{h_i}=\mathit{\text{a priori}}\text{ distribution for }{h}_i \end{gathered}$$

In the example that follows, we demonstrate MCMC inversion for a 6-layered earth, including the half-space being imaged using Rayleigh waves. THe prior distribution assumes all layers have uncorrelated shear wave velocities bounded between $3.5$ and $5km/s$, defined using a uniform distribution, and layered thickness values vary between 15 and 25 km for all the layers above the half-space.

In the following example, we demonstrate MCMC inversion for a 3-layered earth, including the half-space being imaged using DC resistivity method. The prior distribution assumes all layers have uncorrelated resistivities bounded between ${10}^{-1}$ and ${10}^5$, defined using a uniform distribution.

!!! tip Important

The most important thing to be noted here is the specification of the prior distribution, done via:

```julia
modelD = RWModelDistribution(
    Product(
        [Uniform(-1.0, 5.0) for i in eachindex(z)]
    ),
    Product(
        [Uniform(h_lb[i], h_ub[i]) for i in eachindex(h)]
    )
    ...
);
```

in the example.

Copy-Pasteable code

Let's create a synthetic dataset first, with 1% error floors:

vs = [4.1, 4.4, 4.8]
vp = [7.0, 7.5, 7.8]
ρ = [3.0, 3.2, 3.3]
h = [20.0, 20.0] .* 1e3
m_test = RWModel(vs, h, ρ, vp)

T = exp10.(0:0.1:3)

r_obs = forward(m_test, T);

err_c = 0.01 * r_obs.c;
err_resp = SurfaceWaveResponse(err_c)

Now, let's define the a priori distribution with the variable grid points. Note that we only invert for one physical property, i.e., shear wave velocity.

z = collect(0e3:20e3:40e3)
h = diff(z)
h_ub = h .+ 2.5e3
h_lb = h .- 2.5e3

# variable discretization
modelD = RWModelDistribution(Product([Uniform(4.0, 5.0) for i in eachindex(z)]),
    Product([Uniform(h_lb[i], h_ub[i]) for i in eachindex(h)]),
    fill(3.3, length(z)), fill(7.5, length(z)))

then define the likelihood

respD = SurfaceWaveResponseDistribution(normal_dist)

Put everything together for MCMC

n_samples = 10_000
mcache = mcmc_cache(modelD, respD, n_samples, MH())

rw_chain = stochastic_inverse(r_obs, err_resp, T, mcache; progress=true)

Chains MCMC chain (10000×6×1 Array{Float64, 3}):

Iterations        = 1:1:10000
Number of chains  = 1
Samples per chain = 10000
Wall duration     = 12.73 seconds
Compute duration  = 12.73 seconds
parameters        = m0[:m][1], m0[:m][2], m0[:m][3], m0[:h][1], m0[:h][2]
internals         = lp

Summary Statistics
  parameters         mean         std       mcse   ess_bulk   ess_tail      rh ⋯
      Symbol      Float64     Float64    Float64    Float64    Float64   Float ⋯

   m0[:m][1]       4.0527      0.0159     0.0030    24.4100    21.2059    2.63 ⋯
   m0[:m][2]       4.3601      0.0450     0.0070    26.4309    25.3130    2.56 ⋯
   m0[:m][3]       4.8271      0.0196     0.0017    33.8827    47.7080    1.92 ⋯
   m0[:h][1]   18737.6564   1172.3405   224.2910    29.5801    29.9702    2.24 ⋯
   m0[:h][2]   18737.4695   1189.1129   256.0957    21.8875    20.5146    2.16 ⋯
                                                               2 columns omitted

Quantiles
  parameters         2.5%        25.0%        50.0%        75.0%        97.5% 
      Symbol      Float64      Float64      Float64      Float64      Float64 

   m0[:m][1]       4.0388       4.0388       4.0506       4.0506       4.0886
   m0[:m][2]       4.3358       4.3372       4.3372       4.3827       4.4938
   m0[:m][3]       4.8069       4.8237       4.8320       4.8320       4.8389
   m0[:h][1]   17822.1176   18291.8485   18291.8485   18534.2299   22049.3088
   m0[:h][2]   17585.0796   17585.0796   17617.7733   19985.9975   20459.2005

The obtained rw_chain contains the a posteriori distributions that can be saved using JLD2.jl.

using JLD2
JLD2.@save "file_path.jld2" rw_chain

Code for this figure

fig = Figure()
ax = Axis(fig[1, 1])
hm = get_kde_image!(ax, rw_chain, modelD; kde_transformation_fn=log10,
    grid=(m=collect(4.0:0.01:5.0), z=collect(1:1e3:60e3)),
    trans_utils=(m=no_tf, h=no_tf),
    colormap=:binary, colorrange=(-4, 0.0))
Colorbar(fig[1, 2], hm; label="log pdf")

mean_kws = (; color=:seagreen3, linewidth=2)
std_kws = (; color=:red, linewidth=1.5)
# get_mean_std_image!(ax, rw_chain, modelD; confidence_interval=0.9, trans_utils=(m=no_tf, h=no_tf), mean_kwargs=mean_kws,
#     std_plus_kwargs=std_kws, std_minus_kwargs=std_kws, z_points=collect(1:1e3:60e3))
ylims!(ax, [6e4, 0])

plot_model!(ax, m_test; color=:black, linestyle=:dash, linewidth=2, label="true")
Legend(fig[2, :], ax; orientation=:horizontal)

fig # hide

The list of models can then be obtained from chains using, which can then be used to check the fits and perform other diagnostics:

model_list = get_model_list(rw_chain, modelD)

Code for this figure

fig = Figure()
ax1 = Axis(fig[1, 1])

resp_post = forward(model_list[1], T);
for i in 1:(length(model_list) > 1000 ? 1000 : length(model_list))
    forward!(resp_post, model_list[i], T)
    plot_response!([ax1], T, resp_post; alpha=0.4, color=:gray)
end

plot_response!([ax1], T, r_obs; errs=err_resp, plt_type=:errors, whiskerwidth=10)
plot_response!([ax1], T, r_obs; plt_type=:scatter, label="true")