17  Physics-Based ML for Inversion

Key references
  • Physics-guided neural networks (PGNN) — hybrid models that compose a physics-based prediction with a neural correction, with physical-consistency penalties in the loss (Karpatne et al., 2017).
  • Implicit neural representations (INRs) — coordinate-based networks that represent continuous signals, originally popularized for 3D scene reconstruction (Mildenhall et al., 2022; Sitzmann et al., 2020).
  • Fourier feature encodings — overcoming spectral bias of MLPs by mapping inputs through Fourier features, enabling learning of high-frequency detail (Tancik et al., 2020).
  • PINNs for inversion — embedding PDE constraints in the loss to simultaneously solve forward and inverse problems (Raissi et al., 2019).
  • Physics-informed review — broad treatment of physics-informed ML including PDE-constrained optimization and hybrid approaches (Karniadakis et al., 2021).
  • Self-adaptive PINNs — learning the loss-weighting factors alongside the network parameters to balance PDE residual and data terms (McClenny & Braga-Neto, 2023).

In the previous chapter we used PINNs to solve a PDE: the network represented the solution field and was trained purely from the governing equation. Here we tackle a harder and more practically important problem: geophysical inversion. The goal is to recover an unknown physical property (e.g. subsurface conductivity, seismic velocity, or density) from measured data, while honoring a known PDE that relates the property to the observations.

Two complementary styles of physics-based ML have emerged for this task. They differ in what the network represents and how the physics enters the model.

PGNN — Physics-Guided Neural Network INR — Implicit Neural Representation
What is the network? A hybrid: classical physics model + neural correction A continuous coordinate-based MLP for the unknown field
How does physics enter? Through (i) an explicit physics prediction blended into the architecture and (ii) a physical-consistency penalty Through a PDE-based forward operator that turns predicted properties into predicted observations
Typical motivation Improve a known imperfect physics model with data Replace voxel/mesh inversion with a mesh-free continuous prior
Where it shines Calibration, model-error correction, hybrid forecasting High-resolution recovery from sparse data, mesh-free regularization

We work through both on the same 1D toy problem, so that the difference is methodological, not problem-dependent.

17.1 Shared forward problem

Throughout this chapter we use steady-state 1D diffusion as a stand-in for elliptic geophysical PDEs (electrical resistivity, Darcy flow, steady heat conduction):

\[ -\frac{d}{dx}\!\left[\kappa(x)\,\frac{du}{dx}\right] = f(x), \qquad x \in [0, 1], \qquad u(0) = u(1) = 0 \]

with a known source \(f(x) = \sin(2\pi x)\) and a true but unknown coefficient field

\[ \kappa_\star(x) = 1 + 0.5\sin(2\pi x). \]

We observe noisy values of \(u\) at a small number of sensor locations and try to recover \(\kappa(x)\).

using Lux, Random, Optimisers, Zygote, Statistics, LinearAlgebra, Printf, CairoMakie

rng = Xoshiro(42)

#-----shared discretization + forward solver------
nx     = 200
x_grid = Float32.(range(0, 1, length = nx))
dx     = x_grid[2] - x_grid[1]

κ_true(x)   = 1.0f0 + 0.5f0 * sin(2π * x)
f_source(x) = sin(2π * x)

# Differentiable FD solver: -(κ u')' = f, Dirichlet BCs.
# Returns u on the same grid for any κ-on-grid vector.
function poisson_kappa(κvec, fvec, dx)
    n  = length(κvec)
    A  = zeros(eltype(κvec), n, n)
    A[1, 1] = 1; A[n, n] = 1
    for i in 2:n-1
        κp = 0.5 * (κvec[i] + κvec[i+1])
        κm = 0.5 * (κvec[i] + κvec[i-1])
        A[i, i-1] = -κm / dx^2
        A[i, i]   = (κp + κm) / dx^2
        A[i, i+1] = -κp / dx^2
    end
    b = copy(fvec); b[1] = 0; b[n] = 0
    return A \ b
end
#-----generate synthetic observations------
κ_vec   = κ_true.(x_grid)
f_vec   = f_source.(x_grid)
u_truth = poisson_kappa(κ_vec, f_vec, dx)

n_obs   = 20
obs_idx = sort(rand(rng, 2:nx-1, n_obs))
x_obs   = x_grid[obs_idx]
u_obs   = u_truth[obs_idx] .+ 0.005f0 .* randn(rng, Float32, n_obs)

The two inversion methods below will both try to recover \(\kappa_\star\) from u_obs.

17.2 Method 1 — PGNN: a hybrid physics + neural model

A physics-guided neural network (Karpatne et al., 2017) does not throw the existing physics away. It parameterizes the unknown property as a small multiplicative correction to a trusted baseline value \(\bar\kappa\),

\[ \kappa_\phi(x) \;=\; \bar\kappa \cdot \exp\!\bigl(\mathrm{MLP}_\phi(x)\bigr), \]

and trains \(\phi\) to fit data while also penalizing violations of a physical-consistency criterion — here, drift away from the baseline and lack of smoothness.

For our inverse problem we instantiate this idea as follows:

  1. Physics baseline. Pick a trusted constant guess \(\bar\kappa = 1\).
  2. Neural property network. \(\kappa_\phi(x) = \bar\kappa \cdot \exp(\mathrm{MLP}_\phi(x))\) guarantees positivity and is centred on the baseline.
  3. Data loss through physics. Solve the PDE with \(\kappa_\phi\) to get \(u_{\mathrm{pred}}\) and compare to u_obs.
  4. Physics-consistency penalty. Penalize the deviation of \(\kappa_\phi\) from \(\bar\kappa\) (a Tikhonov-style “stay close to physics unless data says otherwise”) and the squared \(L^2\) norm of \(\kappa_\phi'(x)\) for smoothness.

The neural network is guided by the baseline physics; it is not free to override the PDE.

#-----PGNN: NN predicts log-correction to a baseline κ̄=1------
pgnn = Chain(
    Dense(1 => 32, tanh),
    Dense(32 => 32, tanh),
    Dense(32 => 1),
)
ps_p, st_p = Lux.setup(rng, pgnn)

function kappa_pgnn(xs, ps, st)
    z, _ = pgnn(reshape(xs, 1, :), ps, st)
    return vec(exp.(z))      # κ = κ̄ · exp(MLP(x)) with κ̄ = 1
end

function pgnn_loss(ps_p, st_p)
    κ_pred = kappa_pgnn(x_grid, ps_p, st_p)
    u_pred = poisson_kappa(κ_pred, f_vec, dx)

    loss_data = mean((u_pred[obs_idx] .- u_obs) .^ 2)

    # physics-consistency: stay close to baseline κ̄=1 + smoothness penalty
    loss_baseline = mean((κ_pred .- 1.0f0) .^ 2)
    loss_smooth   = mean(diff(κ_pred) .^ 2) / dx^2

    return loss_data + 1f-2 * loss_baseline + 1f-4 * loss_smooth
end
st_o_p = Optimisers.setup(Adam(5f-3), ps_p)
for epoch in 1:1500
    (loss_v,), grads = Zygote.withgradient(p -> pgnn_loss(p, st_p), ps_p)
    st_o_p, ps_p = Optimisers.update(st_o_p, ps_p, grads[1])
    (epoch == 1 || epoch % 300 == 0) && @printf "PGNN epoch %4d  loss = %.5f\n" epoch loss_v
end
κ_pgnn = kappa_pgnn(x_grid, ps_p, st_p)
u_pgnn = poisson_kappa(κ_pgnn, f_vec, dx)

The PGNN ends up close to \(\kappa_\star\) because (i) the data term pulls it towards solutions explaining u_obs through the same governing physics, and (ii) the consistency terms prevent unphysical excursions when the data is silent.

17.3 Method 2 — INR: an implicit neural representation of the property field

An INR (Mildenhall et al., 2022; Sitzmann et al., 2020) takes a different stance: there is no baseline physics blended into the architecture. Instead, the unknown property field itself is reparameterized as a coordinate-based MLP,

\[ \kappa_\theta(x) = \mathrm{softplus}\!\bigl(\mathrm{MLP}_\theta(\gamma(x))\bigr), \]

where \(\gamma\) is a fixed Fourier feature encoding (Tancik et al., 2020),

\[ \gamma(x) = \bigl[\sin(2\pi B x),\; \cos(2\pi B x)\bigr], \]

with \(B \in \mathbb{R}^{m \times 1}\) a fixed random frequency matrix. The frequency scale of \(B\) directly controls the maximum spatial frequency the network can resolve.

The physics enters not through the architecture but only through the forward operator used in the loss:

\[ \mathcal{L}(\theta) = \frac{1}{N_d}\sum_{i=1}^{N_d}\bigl(u_{\mathrm{pred}}(x_i;\kappa_\theta) - u_i^{\mathrm{obs}}\bigr)^2. \]

#-----INR: Fourier-feature MLP, no baseline physics in the architecture------
n_feat   = 32
B_matrix = randn(rng, Float32, n_feat, 1) .* 4.0f0

function fourier_encode(xs)
    proj = B_matrix * reshape(xs, 1, :)
    return vcat(sin.(2π .* proj), cos.(2π .* proj))
end

inr = Chain(
    Dense(2 * n_feat => 64, tanh),
    Dense(64 => 64, tanh),
    Dense(64 => 1),
)
ps_i, st_i = Lux.setup(rng, inr)

function kappa_inr(xs, ps, st)
    enc   = fourier_encode(xs)
    z, _  = inr(enc, ps, st)
    return vec(Lux.NNlib.softplus.(z))
end

function inr_loss(ps_i, st_i)
    κ_pred = kappa_inr(x_grid, ps_i, st_i)
    u_pred = poisson_kappa(κ_pred, f_vec, dx)
    return mean((u_pred[obs_idx] .- u_obs) .^ 2)
end
st_o_i = Optimisers.setup(Adam(1f-3), ps_i)
for epoch in 1:2000
    (loss_v,), grads = Zygote.withgradient(p -> inr_loss(p, st_i), ps_i)
    st_o_i, ps_i = Optimisers.update(st_o_i, ps_i, grads[1])
    (epoch == 1 || epoch % 400 == 0) && @printf "INR  epoch %4d  loss = %.5f\n" epoch loss_v
end
κ_inr = kappa_inr(x_grid, ps_i, st_i)
u_inr = poisson_kappa(κ_inr, f_vec, dx)

The INR converges to \(\kappa_\star\) purely through implicit regularization — the spectral bias of MLPs together with a controlled Fourier feature bandwidth — without ever specifying an explicit prior.

17.4 Side-by-side comparison

fig = Figure(size = (820, 320))

ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = "κ(x)",
           title = "Recovered diffusion coefficient")
lines!(ax1, x_grid, κ_true.(x_grid), color = :black,     linewidth = 2,  label = "true κ(x)")
lines!(ax1, x_grid, κ_pgnn,          color = :firebrick, linewidth = 2,
       linestyle = :dash, label = "PGNN")
lines!(ax1, x_grid, κ_inr,           color = :steelblue, linewidth = 2,
       linestyle = :dashdot, label = "INR")
axislegend(ax1, position = :rb, labelsize = 9)

ax2 = Axis(fig[1, 2], xlabel = "x", ylabel = "u(x)",
           title = "Reconstructed state field")
lines!(ax2, x_grid, u_truth, color = :black,     linewidth = 2,  label = "true u(x)")
lines!(ax2, x_grid, u_pgnn,  color = :firebrick, linewidth = 2,
       linestyle = :dash, label = "PGNN")
lines!(ax2, x_grid, u_inr,   color = :steelblue, linewidth = 2,
       linestyle = :dashdot, label = "INR")
scatter!(ax2, x_obs, u_obs, color = :black, markersize = 6, label = "obs.")
axislegend(ax2, position = :rt, labelsize = 9)
fig

17.5 How they relate to classical inversion

Classical Tikhonov PGNN INR
Unknown Mesh vector \(\mathbf{m} \in \mathbb{R}^M\) NN weights + physics baseline NN weights only
Regularization Explicit penalty on \(\|\nabla m\|^2\), TV, etc. Explicit physics-consistency penalty Implicit (architecture + spectral bias)
Resolution Fixed by mesh Continuous, tied to a base model Continuous, controlled by Fourier scale
Gradient Adjoint-state method Automatic differentiation through physics Automatic differentiation through physics
Best when Mature physics, no learning needed Trustworthy but imperfect physics Rich data, weak priors, mesh-free preferred

17.6 When to choose which

  • Use a PGNN when there is a partially trusted physics model you do not want to throw away. Calibration of geothermal reservoirs against legacy simulators, bias-correction of operational weather models, and seismic processing pipelines with empirical attribute corrections all fit this pattern.
  • Use an INR when the dominant prior is “the field should be smooth (and continuous) in space” rather than “the field should be close to model \(\bar m\).” 3D gravity/magnetic inversion, electrical resistivity tomography, and full-waveform velocity recovery all benefit from a mesh-free continuous prior with a tunable resolution.
  • Combine both when you have both a baseline physics model and spatial smoothness as priors. The PGNN’s baseline-fidelity term and the INR’s Fourier feature regularization are not in conflict; you can stack them.

17.7 Geoscience milestones

  • Original PGNN — Karpatne et al. (2017) paired a process-based lake-temperature model with a neural correction and physics-consistency penalties; this is the canonical reference for the hybrid physics + neural approach.
  • Seismic full-waveform inversion with PINN-style constraints — Rasht-Behesht et al. (2022) is a representative reference for using the wave-equation residual to regularize a neural representation of the velocity model.
  • Geothermal reservoir characterization — Sun et al. (2023) demonstrates physics-informed INR-style recovery of permeability fields from temperature and pressure observations.
Looking ahead

This chapter combined a single PDE with a single set of observations. Real geophysical inversions stack multiple data types (e.g. gravity + seismic) over complex 3D geometries. Both PGNN and INR frameworks extend naturally: the PGNN by adding more physics terms and baselines, the INR by moving to higher-dimensional coordinate inputs and richer Fourier feature spectra.