Statistical Learning for Scientific Imaging and Physical Measurements: Geosciences

Contents

• Introduction and overview
• 0. Review: linear algebra, signals, probability and optimisation
• 1. Measurement systems: images and logs
• 2. Compression
• 3. Super-resolution
• 4. Quality measures and ranking
• 5. Data alignment
• 6. Contrast enhancement and inpainting
• 7. Classification and segmentation
• 8. Inversion
• 9. From prototype to deployment: the industrial challenges
• 10. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: students with a good M1 in mathematics, computer science or physics who want to learn how statistical learning is adapted to images and signals produced by physical sensors, with subsurface exploration and energy as the running case.

Empirical first: on synthetic well logs built to mimic real ones, a random forest reaches 97 % accuracy under random cross-validation and 49 % on wells it has never seen (code 7.3). Neither the model nor the data are at fault; the evaluation is. This course is about such traps, and about the physics that explains them.

Geophysical tools never photograph the rock. They send a seismic, electromagnetic, ultrasonic or mechanical field into the ground, record a response, and leave the rest to inference. The course follows that chain from acquisition to interpretation: describe the measurement, compress it, sharpen it, judge its quality, align it, enhance and repair it, classify and segment it, and finally invert it.

At each step it puts a classical method, a statistical-learning method and a deep-learning method side by side, and it asks what changes when the model leaves the laboratory. Closely related MVA courses are Deep Learning, Image Denoising, Sub-Pixel Image Processing and Machine Learning for Time Series; in semester 2, Inverse Problems and Imaging and Remote Sensing Data take several threads further.

What you will be able to do

• Write the forward model of a measurement system (well log, borehole image, seismic trace, ultrasonic echo) and read its resolution limits.
• Compress images and logs with transform coding, PCA and learned codecs, and predict the distortion from the rate.
• Formulate super-resolution as an ill-posed estimation problem and explain why mean-squared-error networks blur.
• Choose and compute quality measures (PSNR, SSIM, noise estimates) and ranking metrics, and fit scores from pairwise expert preferences.
• Align logs and images with cross-correlation, phase correlation, dynamic time warping and Procrustes registration.
• Enhance contrast with histogram methods and fill gaps by variational and learned inpainting.
• Train and validate facies classifiers and segmentation networks without leaking information across wells.
• Solve linear and sparse inverse problems (Tikhonov, ISTA) and combine them with learning and physical constraints.
• Diagnose distribution shift, quantify uncertainty with conformal prediction, and fit a model to compute and memory budgets.
• Plan the path from a research prototype to a tool that interpreters can trust.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com presents modern statistical-learning and deep-learning methods for scientific imaging and for data from physical sensors. Its problems (image reconstruction, super-resolution, data alignment, contrast enhancement, segmentation, inverse problems) are shared by medical imaging, space, industrial vision, robotics and the geosciences.

Through real applications from subsurface exploration and the energy sector, it shows how to adapt AI methods to complex data acquired by advanced measurement systems. These data are images or signals and rest on several physics: electromagnetic, ultrasonic, mechanical and seismic.

The description stresses the challenges met by AI engineers in industry: limited or imperfectly annotated data, distribution shifts between environments, compute and memory constraints, and the need for robustness, interpretability and integration with physical models. It also covers the practices that turn a research prototype into a solution deployable at scale, presented as a concrete introduction to the job of AI engineer.

The content covers the whole processing chain, from data acquisition to interpretation, and combines theoretical foundations with industrial case studies. Several data sets and applications are announced in geology, corrosion and geo-energy.

The Université Paris-Saclay course detail (“Introduction à l’apprentissage statistique pour les géosciences”) describes an introduction to statistical learning for the geosciences with applications in energy. Geophysical measurements are acquired with high-technology sensors in extreme conditions; they describe the subsurface to support critical decisions, and their several physics give interpreters complementary information.

It announces theoretical and practical content on traditional methods, conventional statistical-learning methods and deep-learning methods, and lists the module in the I&I track. Its learning objectives repeat the six themes of master-mva.com, which the table below matches to this page.

Beyond the official syllabus, this page adds the review chapter (ch. 0), formal guarantees that the description does not name (rate–distortion bound, conformal coverage, quantisation error), and a self-study chapter on research directions (ch. 10).

Where this course fits in the programme

Official prerequisites (both pages): knowledge of linear algebra, probability and statistics, optimisation and Python programming. The course is a Modelling course of the I&I track. Master-mva.com announces a cap of 20 students and teaching in French or English with materials in English.

The two official sources disagree on the semester: master-mva.com lists the course in semester 1, the Université Paris-Saclay table in semester 2. This page follows the catalogue (semester 1); check the timetable of the year.

• Refreshers — Probability and Statistics (Gaussian models, estimation), Linear Algebra, Analysis and Optimization (SVD, least squares, gradient methods) and Scientific Python, NumPy and PyTorch (every lab).
• Related courses of the same semester — Deep Learning for networks and training, Introduction to Statistical Learning for validation, Convex Optimization for ISTA and duality, Image Denoising and Sub-Pixel Image Processing for image models, Machine Learning for Time Series for alignment and sequences, Deep Learning and Signal Processing for industrial signals, and Fundamentals of Reproducible Research for deployment practice.
• Courses that build on it or share its tools — Inverse Problems and Imaging, Learning Physical Phenomena (surrogates and uncertainty), Remote Sensing Data, Deep Learning for Image Restoration and Synthesis, Generative Models for Images, Detection Theory and Industrial Applications and Responsible Machine Learning.

Course outline

• Review: linear algebra, signals, probability and optimisation — the SVD, convolution and sampling, MAP estimation, gradient descent and convolutional layers.
• Measurement systems: images and logs — forward models of well logs, borehole images, seismic traces and ultrasonic echoes.
• Compression — transform coding, PCA, quantisation and the rate–distortion limit, then learned codecs.
• Super-resolution — why the problem is ill-posed, what a trained network estimates, and how deep networks are built for it.
• Quality measures and ranking — full-reference and no-reference quality, ranking metrics and learning to rank.
• Data alignment — shifts, warps and rotations between logs, images and seismic.
• Contrast enhancement and inpainting — histogram methods and the filling of gaps.
• Classification and segmentation — facies from logs, features in images, and honest validation.
• Inversion — from measurements back to rock properties, with regularisation, sparsity, physics and learning.
• From prototype to deployment: the industrial challenges — imperfect labels, distribution shift, uncertainty, budgets and deployment.
• Research directions and open problems — self-study on current research.

Session plan

30 h in class = 10 sessions of 3 h. The volume is the Université Paris-Saclay programme table (12 h of lectures, 18 h of labs); master-mva.com announces the same “12 h of lectures + 18 h of labs”. The split into sessions is a proposal: each session pairs about 1 h 12 min of lecture with about 1 h 48 min of lab.

Level 1 = foundations, 2 = core methods, 3 = advanced material; the level never decreases from one session to the next.

Self-study, outside the 30 h (about 3 h in total): ch. 10 Research directions and open problems (≈ 3 h, with Ex. 15 and 16).

How to study this module

• Read ch. 0 before the first session, even quickly: the SVD, the convolution theorem and the MAP formula come back in almost every chapter.
• After each session, run the chapter’s code block, change one parameter, and predict the result before looking.
• Do the exercises listed in the session plan in the same week; the ★ ones check the session, the ★★ ones prepare the exam.
• Keep one notebook per chapter with a common data generator (synthetic logs and borehole images): the labs reuse it.
• The environment is Python 3 with NumPy, SciPy, scikit-learn and PyTorch on CPU; every code block runs in less than about 30 s on a laptop.
• Start the project (see Projects and assessment) at the kick-off of session 10, once classification (session 8) and inversion (session 9) have been seen.

0. Review: linear algebra, signals, probability and optimisation

Every method of this course is a variation on one situation: an unknown property, a measurement, and an operator that links them. This review gathers the M1 tools that the later chapters use without stopping: the singular value decomposition, convolution and sampling, the Bayesian reading of regularisation, gradient descent, and the convolutional layer.

Every tool samples, in depth or in time. The sampling interval fixes the highest frequency that can be recovered; super-resolution (ch. 3) is precisely about what lies beyond it.

Two-dimensional layers work the same way on images. Encoder–decoder networks with skip connections (U-Net, ch. 7.2) combine a large receptive field with precise localisation.

Code 0.9 — checking the convolution theorem, aliasing and least squares.

import numpy as np

rng = np.random.default_rng(0)
n = 64
x, h = rng.standard_normal(n), rng.standard_normal(n)

# circular convolution, directly and through the DFT (Proposition 0.4)
direct = np.array([sum(x[m] * h[(k - m) % n] for m in range(n)) for k in range(n)])
via_fft = np.real(np.fft.ifft(np.fft.fft(x) * np.fft.fft(h)))
print("max |direct - FFT| =", np.abs(direct - via_fft).max())

# aliasing (Theorem 0.5): 70 Hz sampled at 100 Hz looks like 30 Hz
t = np.arange(50) / 100.0
gap = np.abs(np.cos(2 * np.pi * 70 * t) - np.cos(2 * np.pi * 30 * t)).max()
print("max |cos(70 Hz) - cos(30 Hz)| on the grid =", gap)

# least squares through the pseudo-inverse (Definition 0.1)
A = rng.standard_normal((20, 5))
y = rng.standard_normal(20)
x_ls = np.linalg.pinv(A) @ y
print("normal equations residual:", np.abs(A.T @ (A @ x_ls - y)).max())

All three printed numbers are at the level of rounding errors (about 10^{-14}): the two convolutions agree, the 70 Hz and 30 Hz cosines coincide on a 100 Hz grid, and the pseudo-inverse solves the normal equations.

1. Measurement systems: images and logs

A geoscientist never observes porosity, a fluid or a fault directly. A tool emits a physical field, records a response, and the property is inferred. This chapter models that chain, so that every later task (compressing, enhancing, aligning, classifying, inverting) is stated on the right object.

1.1 Measurement as a forward model

Computing y from x is the forward problem, a simulation that is usually well-posed. Recovering x from y is the inverse problem of ch. 8. Most “images” of this course are not photographs: each pixel is a physical measurement with a unit, a resolution and an error.

1.2 Well logs

1.3 Borehole images

An imaging tool records many logs around the borehole wall at once. Resistivity tools press several pads against the wall and cover part of the circumference; rotating ultrasonic transducers cover all of it. Unrolling the cylinder gives an image I(z, \varphi), with depth z downwards and azimuth \varphi \in [0, 2\pi) across.

1.4 Seismic data: the convolutional model

The scale gap is the central fact of subsurface data. Seismic resolves tens of metres over kilometres; logs resolve decimetres along one well; borehole images resolve centimetres around it. Tying these scales together is the job of alignment (ch. 5) and inversion (ch. 8).

1.5 Ultrasonic pulse–echo and casing corrosion

An ultrasonic transducer emits a short pulse and records the echo. The travel time t to a reflector at distance d in a fluid of sound speed c gives d = ct/2; rotating the transducer maps the borehole radius around the wall, which reveals ovalisation and breakouts.

With c \approx 5900 m/s in steel and e = 10 mm, f_1 \approx 295 kHz. Corrosion thins the casing, so the resonance frequency rises; measured around and along the well, it becomes a thickness image e(z, \varphi) on which metal loss is segmented (ch. 7.2).

Electromagnetic inspection rests on another physics: eddy currents induced in the metal decay over the skin depth \delta = \sqrt{\frac{2}{\omega \mu \sigma}}, which depends on the frequency \omega, the magnetic permeability \mu and the conductivity \sigma. Lower frequencies probe deeper.

1.6 What makes these data hard

• Physical units and dynamic range. Resistivity spans orders of magnitude, so it is modelled in \log scale; each curve has its own noise and artefacts.
• Several scales. Each tool sees its own resolution and depth of investigation (Example 1.9).
• Labels are interpretations. Facies, faults or fractures are drawn by experts who may disagree, and labelled wells are few (ch. 9.1).
• Every well is a new environment. Tool generation, borehole size, drilling mud and geology change the input distribution (ch. 7.1, ch. 9.2).
• Volumes and budgets. 3D surveys are large and downhole links are slow (ch. 2, ch. 9.4).
• Costly decisions. Predictions feed drilling and production decisions, so uncertainty and interpretability matter (ch. 9.3).

Code 1.11 — a synthetic seismogram, tuning, and a sinusoid fit.

import numpy as np

# Ricker wavelet, 30 Hz peak frequency, 1 ms sampling
dt, f = 0.001, 30.0
t = np.arange(-0.1, 0.1 + dt, dt)
w = (1 - 2 * (np.pi * f * t) ** 2) * np.exp(-(np.pi * f * t) ** 2)

# a bed of impedance Z2 inside a background Z1, velocity 3000 m/s
v, Z1, Z2 = 3000.0, 6.0e6, 7.5e6
R = (Z2 - Z1) / (Z2 + Z1)                  # Definition 1.7
print(f"R = {R:.4f}")
for h in [5, 10, 15, 20, 25, 30, 40, 60]:  # thickness in metres
    r = np.zeros(1000)
    k = int(round(2 * h / v / dt))         # two-way time in samples
    r[400], r[400 + k] = R, -R             # top and base reflections
    s = np.convolve(r, w, mode="same")     # convolutional model
    print(f"h = {h:2d} m   peak amplitude = {np.abs(s).max():.4f}")

# Theorem 1.5: fit a sinusoid to noisy picks on a borehole image
rng = np.random.default_rng(1)
radius = 0.108                              # 8.5 in borehole, metres
dip, azi, z0 = np.deg2rad(30), np.deg2rad(120), 1500.0
phi = np.sort(rng.uniform(0, 2 * np.pi, 40))
z = z0 + radius * np.tan(dip) * np.cos(phi - azi) + 0.003 * rng.standard_normal(40)
X = np.column_stack([np.ones_like(phi), np.cos(phi), np.sin(phi)])
a, b, c = np.linalg.lstsq(X, z, rcond=None)[0]
print("dip  =", np.rad2deg(np.arctan(np.hypot(b, c) / radius)).round(2), "deg")
print("azimuth =", (np.rad2deg(np.arctan2(c, b)) % 360).round(1), "deg")

With 3 mm of picking noise on 40 points, the fit returns a dip of 30.06° and an azimuth of 120.0°, against 30° and 120° in truth.

2. Compression

Subsurface data are large, and the channels that carry them can be tiny. A 3D seismic survey can reach terabytes; while drilling, data computed downhole travel to the surface by mud-pulse telemetry, at rates of the order of a few bits to a few tens of bits per second. This chapter asks how many bits a signal really needs for a given error.

2.1 Rate, distortion and quality

PSNR is only a starting point: ch. 4 discusses measures that track structure and the downstream task.

2.2 Transform coding

The classical recipe has three steps: transform the signal into coefficients that concentrate its energy, quantise the coefficients, and code the quantised values without loss (entropy coding).

PCA must transmit its basis, which is costly. The discrete cosine transform (DCT) needs none, and for strongly correlated stationary signals it is close to the PCA basis (Ahmed, Natarajan and Rao, 1974); this is why JPEG uses it. Wavelets (Mallat, 2009) suit piecewise-smooth signals, such as blocky logs, because only coefficients near the jumps are large.

2.3 Quantisation and the rate–distortion limit

2.4 Learned compression

• Linear autoencoders. A network x \mapsto W_2 W_1 x trained with squared error learns the PCA subspace: all its local minima are global and span the top principal directions (Baldi and Hornik, 1989). Theorem 2.4 is therefore the baseline of every learned codec.
• Nonlinear transform coding. An encoder f_\theta, a quantiser Q, an entropy model p_\psi of the quantised code, and a decoder g_\phi are trained together on a rate–distortion loss (Ballé, Laparra and Simoncelli, 2017).
• The training trick. Rounding has zero gradient, so during training Q is replaced by additive uniform noise on [-1/2, 1/2], which has the same error statistics as Proposition 2.5.

The loss reads \mathcal{L} = \mathbb{E}\big[-\log_2 p_\psi(Q(f_\theta(x)))\big] + \lambda\, \mathbb{E}\, \|x - g_\phi(Q(f_\theta(x)))\|^2 , where the first term is the code length after entropy coding and \lambda selects a point on the rate–distortion curve.

Industrial constraints shape the design. The encoder runs downhole on small, hot, low-power hardware, while the decoder runs at the surface and can be large (ch. 9.4). And a codec must not invent geology: a compressed image is judged on the features an interpreter picks from it (ch. 4.1), not only on its PSNR.

Code 2.8 — transform coding of a synthetic borehole image.

import numpy as np
from scipy.fft import dctn, idctn

rng = np.random.default_rng(0)
H, W = 256, 256
z, j = np.mgrid[0:H, 0:W].astype(float)                    # depth, azimuth index
phi = 2 * np.pi * j / W
# synthetic unrolled borehole image: dipping layers + two fractures + noise
img = np.sin(2 * np.pi * (z + 12 * np.cos(phi - 1.0)) / 23)
for z0, amp in [(80, 30), (170, 18)]:
    img += 2.0 * np.exp(-((z - z0 - amp * np.cos(phi - 2.5)) / 1.5) ** 2)
img += 0.1 * rng.standard_normal((H, W))
img = (img - img.min()) / (img.max() - img.min())          # range [0, 1]

def psnr(x, y, peak=1.0):                                  # Definition 2.1
    return 10 * np.log10(peak ** 2 / np.mean((x - y) ** 2))

c = dctn(img, norm="ortho")                                # orthonormal 2D DCT
for frac in [0.01, 0.05, 0.10]:
    thr = np.quantile(np.abs(c), 1 - frac)                 # keep the largest (Thm 2.3)
    rec = idctn(np.where(np.abs(c) >= thr, c, 0), norm="ortho")
    print(f"DCT, keep {frac:4.0%}: PSNR = {psnr(img, rec):5.2f} dB")

# PCA on 8x8 patches (Theorem 2.4): k of 64 coefficients per patch
P = img.reshape(32, 8, 32, 8).transpose(0, 2, 1, 3).reshape(-1, 64)
mu = P.mean(0)
U, s, Vt = np.linalg.svd(P - mu, full_matrices=False)
for k in [1, 3, 6]:
    rec = (P - mu) @ Vt[:k].T @ Vt[:k] + mu
    rec = rec.reshape(32, 32, 8, 8).transpose(0, 2, 1, 3).reshape(H, W)
    print(f"PCA, k = {k} ({k / 64:4.1%}): PSNR = {psnr(img, rec):5.2f} dB")

# uniform quantisation (Proposition 2.5): measured vs predicted MSE
for bits in [4, 6, 8]:
    step = 1.0 / 2 ** bits
    q = (np.floor(img / step) + 0.5) * step
    print(f"{bits} bits: MSE = {np.mean((img - q) ** 2):.2e}, step^2/12 = {step ** 2 / 12:.2e}")

3. Super-resolution

Logs are averaged over decimetres, seismic over tens of metres, and some imaging tools are coarser than the features interpreters need. Super-resolution estimates a finer signal from a coarser one. This chapter shows that the problem has infinitely many exact solutions, and that what a trained network returns depends on its loss.

3.1 The degradation model

Super-resolution is therefore a choice inside an affine space of dimension at least N(1 - 1/s), and only prior knowledge can make it. In Fourier terms (Theorem 0.5), frequencies above the coarse Nyquist limit are folded onto lower ones and cannot be told apart from them without a prior.

3.2 Classical estimators

• Interpolation (nearest, linear, cubic, sinc) assumes the signal is smooth or band-limited; it never restores what the blur removed.
• Regularised inversion chooses the smoothest consistent signal, with D a finite-difference matrix: \hat x = \arg\min_x \|y - Ax\|^2 + \lambda \|Dx\|^2 = (A^\top A + \lambda D^\top D)^{-1} A^\top y . The matrix is invertible when \ker A \cap \ker D = \{0\}, since the quadratic form then vanishes only at zero.
• Example-based regression learns the map from coarse to fine patches on pairs where both are known.

3.3 What a trained estimator converges to

Blau and Michaeli (2018) proved a perception–distortion trade-off: at a given level of realism, measured as a divergence between the distribution of outputs and that of real signals, the lowest achievable distortion increases as realism is pushed further. Adversarial networks such as SRGAN (Ledig et al., 2017) gain realism by giving up PSNR.

In geoscience, realism without fidelity is dangerous: a sharp but invented fracture is worse than a blurred real one. Two safeguards follow. First, check data consistency: the residual \|A\hat x - y\| should match the noise level. Second, prefer methods that report uncertainty, such as posterior sampling (ch. 10).

3.4 Deep super-resolution networks

• SRCNN (Dong et al., 2016) upsamples by bicubic interpolation, then applies three convolutional layers read as patch extraction, nonlinear mapping and reconstruction.
• Residual learning (VDSR, Kim, Lee and Lee, 2016) predicts the difference between the fine signal and the interpolated input, which is mostly high-frequency and easier to learn with a deep network.
• Sub-pixel convolution (Shi et al., 2016) computes in the coarse grid and rearranges s^2 channels into an image s times larger, which is cheaper.

The training pairs are the weak point. Real high-resolution geology is rarely observed, so pairs are usually synthesised by applying a model of the tool response to high-resolution data. The network then learns that model, not the real tool, and the gap between them is a distribution shift (ch. 9.2). When two tools of different resolution log the same interval, their depth-aligned data (ch. 5) give real pairs.

Code 3.7 — super-resolution of a synthetic log (about 10 s on a laptop CPU).

import numpy as np, torch
from scipy.sparse import diags, csr_matrix
from scipy.sparse.linalg import spsolve

rng = np.random.default_rng(0)
torch.manual_seed(0)
torch.set_num_threads(1)
N, s, width, sigma = 512, 4, 8, 0.05      # fine samples, factor, tool width, noise

def blocky_log():                          # piecewise-constant "true" log
    edges = np.cumsum(rng.geometric(1 / 12, size=N))
    vals = rng.standard_normal(len(edges) + 1)
    return vals[np.searchsorted(edges, np.arange(N), side="right")]

A = np.zeros((N // s, N))                  # decimation o box blur (Def. 3.1)
for i in range(N // s):
    c = i * s + s // 2
    A[i, max(0, c - width // 2):c + width // 2] = 1 / width
def measure(x):
    return A @ x + sigma * rng.standard_normal(N // s)

train = [blocky_log() for _ in range(32)]
x_true = blocky_log()
y = measure(x_true)
rmse = lambda u: np.sqrt(np.mean((u - x_true) ** 2))

# (a) linear interpolation of the coarse samples
centres = np.arange(N // s) * s + s // 2
print("linear interpolation:", round(rmse(np.interp(np.arange(N), centres, y)), 3))

# (b) Tikhonov with a first-difference penalty (closed form, see ch. 8)
As = csr_matrix(A)
D = diags([-1.0, 1.0], [0, 1], shape=(N - 1, N))
x_tik = spsolve((As.T @ As + 0.02 * D.T @ D).tocsc(), As.T @ y)
print("Tikhonov            :", round(rmse(x_tik), 3))

# (c) a small 1D CNN trained on synthetic pairs (x, A x + noise)
net = torch.nn.Sequential(
    torch.nn.Upsample(scale_factor=s, mode="linear"),
    torch.nn.Conv1d(1, 16, 9, padding=4), torch.nn.ReLU(),
    torch.nn.Conv1d(16, 16, 9, padding=4), torch.nn.ReLU(),
    torch.nn.Conv1d(16, 1, 9, padding=4))
X = torch.tensor(np.stack(train), dtype=torch.float32)[:, None]
opt = torch.optim.Adam(net.parameters(), lr=3e-3)
for it in range(400):
    Y = torch.tensor(np.stack([measure(x) for x in train]), dtype=torch.float32)[:, None]
    loss = torch.mean(torch.abs(net(Y) - X))          # L1 loss (Prop. 3.5)
    opt.zero_grad(); loss.backward(); opt.step()
with torch.no_grad():
    x_cnn = net(torch.tensor(y, dtype=torch.float32)[None, None])[0, 0].numpy()
print("1D CNN (L1 loss)    :", round(rmse(x_cnn), 3))

4. Quality measures and ranking

Every later decision rests on comparisons: is this reconstruction better than that one, is this interval of log reliable, which of ten models should be deployed? This chapter defines measures of quality, with and without a reference, and the tools that turn comparisons into rankings.

4.1 Full-reference measures

MSE and PSNR (Definition 2.1) compare a reconstruction with a reference pixel by pixel. They ignore where the error is and what it looks like, which motivated the structural similarity index.

The best quality measure is the one tied to the decision. For a super-resolved image, it may be the error on dips picked from it (Theorem 1.5); for an enhanced log, the facies accuracy downstream (ch. 7); for an inversion, the error on impedance at a blind well (ch. 8). Optimising a proxy too hard can degrade the task it stands for.

4.2 Quality without a reference

In the field there is no ground truth, yet each log and image must be checked before use. Quality control relies on indicators computed from the data alone.

For a log, the fine-scale coefficients can simply be first differences divided by \sqrt 2. Other no-reference indicators include the coverage fraction of an image, the fraction of saturated values, a sharpness score such as gradient energy, the agreement between redundant sensors, and physical plausibility (a density outside the range of minerals flags a bad interval).

Experts can also grade quality directly. Absolute grades are noisy and drift between experts, whereas a preference between two candidates is easier to give and more consistent. This turns quality assessment into a ranking problem.

4.3 Ranking

RankNet (Burges et al., 2005) replaces free scores by a network s_i = f_\theta(x_i) computed from the image or log x_i, trained with the same pairwise loss. The network then generalises: it scores new data with the preferences learned from experts.

Code 4.8 — SSIM, Kendall’s tau and Bradley–Terry scores.

import numpy as np
from scipy.ndimage import uniform_filter, gaussian_filter
from scipy.stats import kendalltau
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(0)
j = np.arange(128)
x = np.sin(2 * np.pi * (j[:, None] + 8 * np.cos(2 * np.pi * j[None, :] / 128)) / 16)
x = (x + 1) / 2                                    # reference image in [0, 1]

def ssim(a, b, L=1.0, win=7):                      # Definition 4.1, local windows
    C1, C2 = (0.01 * L) ** 2, (0.03 * L) ** 2
    mu_a, mu_b = uniform_filter(a, win), uniform_filter(b, win)
    va = uniform_filter(a * a, win) - mu_a ** 2
    vb = uniform_filter(b * b, win) - mu_b ** 2
    cov = uniform_filter(a * b, win) - mu_a * mu_b
    m = ((2 * mu_a * mu_b + C1) * (2 * cov + C2)) / ((mu_a**2 + mu_b**2 + C1) * (va + vb + C2))
    return m.mean()

psnr = lambda a, b: 10 * np.log10(1 / np.mean((a - b) ** 2))
target_mse = 0.01                                  # three degradations, same MSE
shift = x + np.sqrt(target_mse)                    # brightness offset
noise = rng.standard_normal(x.shape)
noisy = x + noise * np.sqrt(target_mse / np.mean(noise ** 2))
blur = gaussian_filter(x, 1.0)
blur = x + (blur - x) * np.sqrt(target_mse / np.mean((blur - x) ** 2))
for name, im in [("offset", shift), ("noise", noisy), ("blur", blur)]:
    print(f"{name:6s}  PSNR = {psnr(x, im):5.2f} dB   SSIM = {ssim(x, im):.3f}")

# Kendall tau between a metric ranking and an expert ranking of 6 methods
metric_scores = [31.2, 29.8, 30.5, 27.9, 28.4, 32.0]
expert_rank = [2, 4, 1, 6, 5, 3]                   # 1 = best
tau, _ = kendalltau(metric_scores, [-r for r in expert_rank])
print(f"Kendall tau(metric, expert) = {tau:.3f}")

# Bradley-Terry scores from pairwise expert votes (Theorem 4.7)
true_s = np.array([0.0, 0.8, 1.5, -0.5])
pairs, wins = [], []
for _ in range(400):
    i, k = rng.choice(4, 2, replace=False)
    d = np.zeros(4); d[i], d[k] = 1, -1
    pairs.append(d)
    wins.append(rng.random() < 1 / (1 + np.exp(-(true_s[i] - true_s[k]))))
bt = LogisticRegression(fit_intercept=False, C=10.0).fit(np.array(pairs), wins)
s_hat = bt.coef_[0] - bt.coef_[0][0]               # fix the gauge s_0 = 0
print("Bradley-Terry scores:", s_hat.round(2), " true:", true_s)

With 400 simulated votes on four items, the fitted scores are (0, 0.78, 1.56, −0.37) against (0, 0.8, 1.5, −0.5) in truth: logistic regression on difference vectors is exactly Bradley–Terry maximum likelihood, with a mild \ell_2 penalty.

5. Data alignment

Two logging runs in the same well rarely agree in depth: the cable stretches, the tool sticks and slips, and depth references differ. Cores carry their own depths, and seismic lives in time rather than depth. Alignment comes before any learning on several sources, because a misaligned label is a wrong label.

5.1 Constant shifts: cross-correlation and phase correlation

Phase correlation ignores gains and slow trends, extends to 2D translations of images, and handles rotation and scaling after a log-polar resampling. On a borehole image, a rotation of the tool around its axis is a circular shift in azimuth, so it is estimated this way.

5.2 Non-linear alignment: dynamic time warping

When the shift varies with depth, a single lag is not enough.

The constraints carry geology. Monotone steps preserve the stratigraphic order (layers cannot swap), the band encodes the largest plausible depth error, and slope limits forbid absurd stretching. Soft-DTW (Cuturi and Blondel, 2017) replaces the minimum by a smooth minimum, which makes the alignment cost differentiable and usable as a training loss.

5.3 Rigid registration

Rigid registration aligns repeated passes of a tool, matches 3D surveys, and registers images from different tools once a few common features are identified. Learned registration replaces hand-picked features by networks: spatial transformer networks (Jaderberg et al., 2015) predict a transformation inside a network, and VoxelMorph (Balakrishnan et al., 2019) learns deformable registration fields from pairs of images.

5.4 The well-to-seismic tie

The tie aligns a well, measured in depth, with seismic, measured in time. It chains the tools of ch. 1 and of this chapter.

1. Compute the impedance log Z = \rho v from the density and sonic logs.
2. Convert depth to two-way time with the sonic velocities, t(z) = 2\int_0^z dz'/v(z'), corrected by check-shot measurements when they exist.
3. Compute the reflectivity (Definition 1.7) and convolve it with a wavelet to get a synthetic trace.
4. Match the synthetic and the real trace: a bulk shift by cross-correlation (5.1), then a gentle stretch-and-squeeze with a narrow-band DTW (5.2), while estimating the wavelet.

The tie produces the time–depth relation that transfers well labels (facies, impedance) onto the seismic, and the low-frequency model that inversion needs (ch. 8.4).

Code 5.8 — constant shift versus dynamic time warping.

import numpy as np
from scipy.ndimage import gaussian_filter1d

rng = np.random.default_rng(0)
n = 400
gr = gaussian_filter1d(rng.standard_normal(n + 80), 3)       # "true" gamma-ray log
ref = gr[40:40 + n] + 0.05 * rng.standard_normal(n)          # reference run

# second run: depth shift of 5 samples plus a slow stretch (cable effects)
u = np.arange(n)
warp = u + 5 + 15 * np.sin(np.pi * u / n)                     # true depth map
run2 = np.interp(warp + 40, np.arange(n + 80), gr) + 0.05 * rng.standard_normal(n)

# (1) constant shift by cross-correlation (Proposition 5.1)
lags = np.arange(-30, 31)
cc = [np.dot(ref[40:-40], run2[40 + l:n - 40 + l]) for l in lags]
best = lags[int(np.argmax(cc))]
print("cross-correlation lag:", best, "samples")

# (2) dynamic time warping with a Sakoe-Chiba band (Theorem 5.4)
def dtw(a, b, band=25):
    na, nb = len(a), len(b)
    C = np.full((na + 1, nb + 1), np.inf)
    C[0, 0] = 0.0
    for i in range(1, na + 1):
        for j in range(max(1, i - band), min(nb, i + band) + 1):
            d = (a[i - 1] - b[j - 1]) ** 2
            C[i, j] = d + min(C[i - 1, j - 1], C[i - 1, j], C[i, j - 1])
    i, j, path = na, nb, []                                  # backtracking
    while i > 0 and j > 0:
        path.append((i - 1, j - 1))
        k = np.argmin([C[i - 1, j - 1], C[i - 1, j], C[i, j - 1]])
        i, j = (i - 1, j - 1) if k == 0 else ((i - 1, j) if k == 1 else (i, j - 1))
    return C[na, nb], path[::-1]

cost, path = dtw(ref, run2)
# ref sample i matches run2 sample j  <=>  i = warp[j]; compare in the interior
est = {}
for i, j in path:
    est.setdefault(j, []).append(i)
js = np.array([j for j in est if 60 < j < n - 60])
err_dtw = np.mean([abs(np.mean(est[j]) - warp[j]) for j in js])
err_cc = np.mean(np.abs((js - best) - warp[js]))   # run2[j] ~ ref[j - lag]
print(f"mean depth error: constant shift {err_cc:.2f}, DTW {err_dtw:.2f} samples")

The printed lag is −16 in the code’s convention (sample j of the second run matches sample j + 16 of the reference).

6. Contrast enhancement and inpainting

Raw images are often hard to read: a few extreme values crush the contrast, and parts of the wall are not measured at all. This chapter treats both problems: remapping values so that structure becomes visible, and filling gaps with values consistent with their surroundings.

6.1 Contrast enhancement

Borehole images are distributed in two renderings. Static normalisation uses one mapping for the whole well, which keeps values comparable between depths and shows large-scale changes. Dynamic normalisation recomputes the mapping in a sliding depth window, which reveals fractures and fine bedding everywhere but destroys absolute levels. Interpreters need both, for different questions.

Flat, noisy regions have a peaked histogram, so plain equalisation turns their noise into strong texture. CLAHE (Zuiderveld, 1994) clips the histogram, redistributes the excess, computes the map in tiles, and blends neighbouring tiles bilinearly. Learned enhancement networks can mimic the renderings interpreters prefer, at the risk of learning the taste of a single expert (ch. 9.1).

6.2 Inpainting

Harmonic filling is smooth, so it cannot continue thin sinusoids across wide gaps: it averages them away. Other models keep more structure:

• Total variation inpainting (Shen and Chan, 2002) minimises \int_D |\nabla u| and preserves edges, although it joins them by straight segments.
• Exemplar-based methods (Criminisi, Pérez and Toyama, 2004) copy patches from the known region along the isophotes.
• Multiple-point statistics fill borehole images with patterns learned from the measured part (Hurley and Zhang, 2011).
• Geometric models fit the sinusoids of Theorem 1.5 on the visible pieces and extend them across the gaps: physics as a prior.

Deep inpainting has two families. Partial convolutions (Liu et al., 2018) renormalise each convolution over the valid pixels of its window and update the mask from layer to layer; they are trained on large image collections with random masks. Deep Image Prior (Ulyanov, Vedaldi and Lempitsky, 2018) needs no training data. It fits a randomly initialised network g_\theta(z) to the known pixels only, \min_\theta \|M \odot (g_\theta(z) - f)\|^2 , where M is the indicator of known pixels, and stops early. The architecture itself acts as the prior.

To measure inpainting error honestly, hide known pixels with an artificial mask shaped like the real gaps and compare. No other measurement is possible, since the true values in the gaps are never observed.

Code 6.7 — histogram equalisation and harmonic inpainting.

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve

rng = np.random.default_rng(0)
H, W = 96, 128
z, j = np.mgrid[0:H, 0:W].astype(float)
phi = 2 * np.pi * j / W
img = np.sin(2 * np.pi * (z + 10 * np.cos(phi - 1.0)) / 20)       # dipping beds
img = np.exp(1.5 * img) + 0.05 * rng.standard_normal((H, W))     # skewed values

# histogram equalisation (Theorem 6.2): map each value through the empirical CDF
ranks = np.argsort(np.argsort(img.ravel()))
eq = ((ranks + 0.5) / img.size).reshape(H, W)
print("histogram of equalised image (5 bins):", np.histogram(eq, 5, (0, 1))[0])

# pad tool: 4 pads, each covering 75% of its quarter of the circumference
mask = (j % 32) < 24                                              # True = measured
known = np.where(mask, eq, 0.0)

# harmonic inpainting (Theorem 6.5): discrete Laplace equation in the gaps,
# periodic in azimuth, known pixels as Dirichlet data
idx = -np.ones((H, W), dtype=int)
idx[~mask] = np.arange((~mask).sum())
M = lil_matrix(((~mask).sum(),) * 2)
b = np.zeros((~mask).sum())
for (r, c) in zip(*np.nonzero(~mask)):
    k = idx[r, c]
    nbrs = [(r, (c - 1) % W), (r, (c + 1) % W)] + [(rr, c) for rr in (r - 1, r + 1) if 0 <= rr < H]
    M[k, k] = len(nbrs)
    for rr, cc in nbrs:
        if mask[rr, cc]:
            b[k] += known[rr, cc]
        else:
            M[k, idx[rr, cc]] -= 1
filled = known.copy()
filled[~mask] = spsolve(M.tocsr(), b)
rmse = np.sqrt(np.mean((filled[~mask] - eq[~mask]) ** 2))
print(f"gap fraction {1 - mask.mean():.0%}, harmonic RMSE in the gaps = {rmse:.3f}")
print(f"RMSE of a constant fill with the mean = {np.sqrt(np.mean((eq[~mask] - eq[mask].mean()) ** 2)):.3f}")

7. Classification and segmentation

Interpretation assigns classes: a facies to each depth of a well, a label to each pixel of an image. This chapter covers both tasks, with one theme throughout: the validation must imitate the deployment, or the scores will lie.

7.1 Facies classification from logs

A facies is a body of rock with distinctive characteristics (lithology, texture, depositional environment). The task is to predict a facies y_k \in \{1, \dots, K\} at each depth from the logs x(z_k), using a few wells where geologists have labelled facies from cores and logs.

• Features: the logs, their depth derivatives, and windows of neighbouring samples that give context.
• Models, in order of complexity: logistic regression, random forests, gradient boosting such as XGBoost (Chen and Guestrin, 2016), and 1D convolutional or recurrent networks along depth.
• Public benchmark: the tutorial and contest of Hall (2016) released labelled wells for this task.

Random K-fold cross-validation draws test samples from wells that also supply training samples, often from the next depth. Layers are thicker than the sampling interval and the tool smooths the logs (Proposition 1.3), so a test sample usually has a near copy in the training set. The estimate then targets accuracy at known wells, not R_{\mathrm{new}}. Group (well-wise) cross-validation holds out whole wells and estimates R_{\mathrm{new}}.

Code 7.3 — random versus well-wise cross-validation (about 10 s).

import numpy as np
from scipy.ndimage import gaussian_filter1d
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold, GroupKFold, cross_val_score

rng = np.random.default_rng(0)
means = np.array([[0.0, 0.0, 0.0, 0.0],          # 4 facies x 4 logs
                  [1.0, 0.6, -0.4, 0.2],
                  [1.6, -0.3, 0.6, -0.5],
                  [0.4, 1.2, 0.9, 0.8]])
T = np.full((4, 4), 0.01) + np.eye(4) * 0.96      # sticky Markov chain in depth

X, y, wells = [], [], []
for w in range(8):
    f = [rng.integers(4)]
    for _ in range(599):
        f.append(rng.choice(4, p=T[f[-1]] / T[f[-1]].sum()))
    f = np.array(f)
    offset = 0.7 * rng.standard_normal(4)          # well-specific calibration/trend
    logs = means[f] + offset + 0.5 * rng.standard_normal((600, 4))
    logs = gaussian_filter1d(logs, 2, axis=0)      # tool vertical response
    X.append(logs); y.append(f); wells.append(np.full(600, w))
X, y, wells = np.vstack(X), np.concatenate(y), np.concatenate(wells)

clf = RandomForestClassifier(n_estimators=100, random_state=0, n_jobs=1)
rand = cross_val_score(clf, X, y, cv=KFold(4, shuffle=True, random_state=0))
lowo = cross_val_score(clf, X, y, cv=GroupKFold(4), groups=wells)
print(f"random 4-fold accuracy      : {rand.mean():.3f}")
print(f"well-wise 4-fold accuracy   : {lowo.mean():.3f}")

# per-well standardisation, a simple remedy for well-to-well shift (ch. 9.2)
Xn = X.copy()
for w in range(8):
    Xn[wells == w] = (X[wells == w] - X[wells == w].mean(0)) / X[wells == w].std(0)
lowo_n = cross_val_score(clf, Xn, y, cv=GroupKFold(4), groups=wells)
print(f"well-wise, per-well scaling : {lowo_n.mean():.3f}")

Facies follow each other in a geological order: Walther’s law states that facies stacked vertically without a break were deposited side by side. A hidden Markov model with a transition matrix estimated on labelled wells smooths the sample-wise predictions along depth with the Viterbi algorithm. A useful geology-aware score is the accuracy that also accepts the neighbouring facies.

7.2 Semantic segmentation of images

Applications follow the physics of ch. 1:

• Seismic facies: Alaudah et al. (2019) published a benchmark with labelled facies on the public F3 block of the Dutch North Sea.
• Faults: Wu et al. (2019) trained a 3D network on synthetic seismic volumes only (FaultSeg3D) and applied it to field data, an answer to the lack of labels (ch. 9.1).
• Borehole images: fractures, bedding and breakouts. A classical baseline is the Hough transform in the parameter space (z_0, r\tan\theta, \varphi_d) of Theorem 1.5.
• Corrosion: metal-loss regions on casing thickness images (ch. 1.5).

Data augmentation must respect the physics. A circular shift in azimuth (a rotated tool) or a small change of gain is legitimate. Flipping depth reverses the stratigraphic order and is legitimate only if the task is symmetric.

Code 7.8 — a tiny segmentation network with a Dice loss (about 15 s).

import numpy as np, torch
import torch.nn as nn

rng = np.random.default_rng(0)
torch.manual_seed(0)
torch.set_num_threads(1)
H, W = 32, 64
z, j = np.mgrid[0:H, 0:W].astype(float)
phi = 2 * np.pi * j / W

def sample():                       # beds + one conductive fracture (sinusoid)
    beds = np.sin(2 * np.pi * (z + rng.uniform(0, 6) * np.cos(phi - rng.uniform(0, 6))) / rng.uniform(6, 12))
    z0, a, p = rng.uniform(8, 24), rng.uniform(2, 8), rng.uniform(0, 2 * np.pi)
    frac = np.abs(z - z0 - a * np.cos(phi - p)) < 1.0      # Theorem 1.5
    img = 0.6 * beds + 1.0 * frac + 0.4 * rng.standard_normal((H, W))
    return img, frac

def batch(n):
    xs, ys = zip(*[sample() for _ in range(n)])
    return (torch.tensor(np.array(xs), dtype=torch.float32)[:, None],
            torch.tensor(np.array(ys), dtype=torch.float32)[:, None])

def soft_dice_loss(logits, target, eps=1.0):           # Proposition 7.6
    p = torch.sigmoid(logits)
    inter = (p * target).sum((1, 2, 3))
    return 1 - ((2 * inter + eps) / (p.sum((1, 2, 3)) + target.sum((1, 2, 3)) + eps)).mean()

net = nn.Sequential(                                   # tiny fully convolutional net
    nn.Conv2d(1, 16, 5, padding=2, padding_mode="circular"), nn.ReLU(),
    nn.Conv2d(16, 16, 5, padding=2, padding_mode="circular"), nn.ReLU(),
    nn.Conv2d(16, 1, 5, padding=2, padding_mode="circular"))
opt = torch.optim.Adam(net.parameters(), lr=3e-3)
for it in range(150):
    x, t = batch(16)
    out = net(x)
    loss = nn.functional.binary_cross_entropy_with_logits(out, t) + soft_dice_loss(out, t)
    opt.zero_grad(); loss.backward(); opt.step()

x, t = batch(64)                                       # held-out synthetic images
with torch.no_grad():
    pred = (net(x) > 0).float()
inter, union = (pred * t).sum(), ((pred + t) > 0).sum()
iou = (inter / union).item()
print(f"fracture pixels: {t.mean().item():.1%} of the image")
print(f"test IoU = {iou:.3f}, Dice = {2 * iou / (1 + iou):.3f}")

Circular padding in azimuth matters here: the unrolled image is a cylinder, so its left and right edges are neighbours. Note also that the padding wraps in depth, a harmless simplification for these small synthetic patches.

8. Inversion

Inversion goes from measurements back to the property that produced them: reflectivity and impedance from seismic traces, resistivity of the formation from tool readings. Its difficulty is not computing a solution but choosing one that the data cannot contradict and that noise does not destroy. This chapter builds the classical tools and shows where learning and physics enter.

8.1 Ill-posed problems

In finite dimension, the third condition fails gradually rather than absolutely: the ratio \sigma_1 / \sigma_n (the condition number) measures how much noise is amplified. For deconvolution by a band-limited wavelet, the smallest singular values are numerically zero; in code 8.7 they range from 8.3 down to the level of rounding errors, below 10^{-16}. Dividing by them, as A^+ does, amplifies noise without limit.

8.2 Tikhonov regularisation

Small \lambda keeps the bias low and lets the variance explode through the small \sigma_i; large \lambda does the opposite. Three classical rules choose \lambda:

• the discrepancy principle picks \lambda so that \|Ax_\lambda - y\| \approx \sigma\sqrt{m}, the expected noise norm;
• the L-curve plots \log\|x_\lambda\| against \log\|Ax_\lambda - y\| and takes the corner;
• generalised cross-validation minimises a leave-one-out proxy that needs no knowledge of \sigma (Golub, Heath and Wahba, 1979).

The Bayesian reading of Proposition 0.6 also gives uncertainty. With the Gaussian prior and noise of that proposition, the posterior covariance is \Big(\frac{A^\top A}{\sigma^2} + \frac{I}{\tau^2}\Big)^{-1}, which is large exactly along the directions v_i where the data are weak.

8.3 Sparsity: sparse-spike deconvolution

Reflectivity is sparse: a few interfaces in a long trace. Smooth regularisers spread each spike over many samples, so seismic deconvolution uses the \ell_1 norm instead: \min_r \; F(r) = \frac{1}{2}\|y - Wr\|_2^2 + \lambda \|r\|_1 , where W is the convolution matrix of the wavelet.

8.4 From reflectivity to impedance

The recursion needs the starting value Z_0. More deeply, the wavelet has no energy at 0 Hz, so the low frequencies of \ln Z (the slow trends) are not in the seismic at all. Practical impedance inversion adds them from a low-frequency model interpolated between tied wells (ch. 5.4). Data and physics are combined here by construction, not by choice.

Code 8.7 — Tikhonov, ISTA and impedance recovery.

import numpy as np
from scipy.linalg import toeplitz

rng = np.random.default_rng(0)
n, dt, f = 300, 0.002, 25.0
t = np.arange(-40, 41) * dt
w = (1 - 2 * (np.pi * f * t) ** 2) * np.exp(-(np.pi * f * t) ** 2)   # Ricker
col = np.zeros(n); col[:41] = w[40:]
row = np.zeros(n); row[:41] = w[40::-1]
Wm = toeplitz(col, row)                         # convolution matrix, y = W r + noise

# sparse reflectivity from a blocky log-impedance (Proposition 8.5: R = tanh(d lnZ / 2))
lnZ = np.log(6e6) + np.cumsum(np.where(rng.random(n) < 0.04, 0.25 * rng.standard_normal(n), 0))
r = np.tanh(np.diff(lnZ, prepend=lnZ[0]) / 2)
y = Wm @ r + 0.01 * rng.standard_normal(n)

# Tikhonov through the SVD filter factors (Theorem 8.2)
U, s, Vt = np.linalg.svd(Wm)
def tikhonov(lam):
    return Vt.T @ (s / (s ** 2 + lam) * (U.T @ y))
print(f"||r|| = {np.linalg.norm(r):.3f}  (error of the zero estimate)")
for lam in [1e-3, 1e-2, 0.3, 10.0]:
    print(f"Tikhonov, lambda = {lam:g}: error {np.linalg.norm(tikhonov(lam) - r):.3f}")

# ISTA for 1/2 ||y - W r||^2 + lam ||r||_1 (Theorem 8.4)
L = s[0] ** 2
lam, x = 0.03, np.zeros(n)
for k in range(2000):
    g = x + Wm.T @ (y - Wm @ x) / L
    x = np.sign(g) * np.maximum(np.abs(g) - lam / L, 0)          # soft threshold
print(f"ISTA, lambda = {lam}: error {np.linalg.norm(x - r):.3f}, "
      f"{(np.abs(r) > 1e-3).sum()} true spikes, {(np.abs(x) > 1e-3).sum()} non-zeros")

# back to impedance: ln Z_k = ln Z_0 + 2 sum artanh(r_j), needs Z_0 from a well
lnZ_hat = lnZ[0] + 2 * np.cumsum(np.arctanh(np.clip(x, -0.99, 0.99)))
print(f"mean relative impedance error: {np.abs(np.exp(lnZ_hat - lnZ) - 1).mean():.3f}")

8.5 Nonlinear and learned inversion

Many forward operators are nonlinear: induction and laterolog resistivity tools, or wave propagation in full-waveform inversion. Writing J for the Jacobian of \mathcal{M} at the current estimate, the Levenberg–Marquardt step is \delta = (J^\top J + \mu I)^{-1} J^\top \big(y - \mathcal{M}(x)\big), \qquad x \leftarrow x + \delta . It is a Tikhonov problem (Theorem 8.2) for the linearised model, with \mu adapted at each step: decreased after a successful step, increased otherwise. Each iteration needs forward simulations and a Jacobian, usually by the adjoint method, which is why inversion is expensive.

Learning enters in four ways, each with its own failure mode:

1. Direct regression. A network f_\theta : y \mapsto x is trained on simulated pairs (x_i, \mathcal{M}(x_i) + \varepsilon_i). Inference is instantaneous, but the network knows only its training prior, and nothing guarantees that its output explains the data.
2. Unrolled optimisation. LISTA (Gregor and LeCun, 2010) unrolls K ISTA iterations into K layers with learned matrices and thresholds, and reaches a given accuracy with far fewer iterations than ISTA.
3. Learned priors in a variational scheme. Plug-and-play methods (Venkatakrishnan, Bouman and Wohlberg, 2013) replace the proximal step of the regulariser by a trained denoiser, while the data term keeps the physics.
4. Physics-consistent training. A data-consistency loss \|\mathcal{M}(f_\theta(y)) - y\|^2 trains on field data without labels and penalises outputs that contradict the measurement.

The survey of Arridge et al. (2019) gives the mathematical landscape. Whatever the method, two checks remain mandatory: the data residual must sit at the noise level, and predictions must match a blind well that played no part in training or tying.

9. From prototype to deployment: the industrial challenges

A method that works in a notebook on one data set can still fail on the next well, on the next tool, or on the hardware where it must run. This chapter treats the obstacles named by the course description: imperfect labels, distribution shift, robustness and interpretability, compute and memory budgets, and the path to a deployed tool.

9.1 Limited and imperfect labels

The good news is limited. The noisy probabilities are compressed towards \frac12, so they are miscalibrated. Noise that depends on the class or on the sample, which is typical of interpreters who disagree near boundaries, does change the classifier; it needs explicit corrections when flip rates are known (Natarajan et al., 2013).

Practical answers to scarce or imperfect labels:

• collect labels from several interpreters and train on consensus or soft labels;
• use metrics that tolerate small boundary shifts or neighbouring facies (ch. 7.1);
• pre-train on large unlabelled archives with self-supervised objectives, then fine-tune on the few labels;
• generate labelled synthetic data with forward models (ch. 1, ch. 7.2), then treat the synthetic-to-field gap as a distribution shift (9.2).

9.2 Distribution shift

Other answers to shift are standardising each well or tool (code 7.3), domain-adversarial training that makes features indistinguishable across domains (Ganin et al., 2016), and monitoring input statistics in production (9.5).

9.3 Uncertainty, robustness and interpretability

The guarantee rests entirely on exchangeability. In code 9.7, calibrated on the source, the 90 % intervals cover 90.2 % of new source samples but only 76.9 % of the shifted target. Weighted conformal prediction (Tibshirani et al., 2019) restores coverage under covariate shift with the weights of Proposition 9.3.

• Ensembles. Several networks trained from different initialisations (Lakshminarayanan, Pritzel and Blundell, 2017) disagree where data are scarce; their spread is a practical estimate of model uncertainty.
• Robustness. Test against realistic perturbations (noise levels, gaps, gains, a different tool) before worrying about adversarial ones.
• Interpretability. Permutation importance measures the drop in score when one log is shuffled; check that the model relies on inputs that make physical sense. When a simple model is almost as good, it is often the better product.

9.4 Compute and memory constraints

Storing weights in 8-bit integers instead of 32-bit floats divides memory by 4 and enables integer arithmetic (Jacob et al., 2018); the step is set per tensor or per channel from the range of the weights. Distillation trains a small student network on the soft outputs of a large teacher (Hinton, Vinyals and Dean, 2015). Pruning removes weights or channels of small importance.

Activations, not weights, often dominate memory. A seismic cube of 1000^3 samples in 32-bit floats occupies 4 GB, and one 32-channel feature map of the same size would need 128 GB. Large volumes are therefore processed in overlapping tiles. The overlap must exceed the receptive field (Definition 0.8), and predictions are blended with weights that vanish at tile borders to avoid seams.

9.5 From prototype to a deployable solution

• Frame the decision. Agree with the users on what the output will change and on the acceptable error, in their units.
• Own the data pipeline. Units, depth references, versions of tools and of processing; every training set is a versioned artefact.
• Validate like you deploy. Blind wells or blind fields (ch. 7.1), a classical baseline from each chapter, and task-based metrics (ch. 4.1).
• Ship uncertainty and quality flags. Conformal intervals or ensemble spread (9.3), and the quality-control indicators of ch. 4.2.
• Keep a human in the loop. Interpreters validate or correct, and their corrections become new labels.
• Monitor and retrain. Track input statistics against the training distribution (9.2), and retrain when they drift.
• Document. A model card states intended use, data, metrics and known limits (Mitchell et al., 2019); reproducibility practices are the subject of Fundamentals of Reproducible Research.

Code 9.7 — importance weighting and conformal coverage under shift.

import numpy as np
from sklearn.linear_model import LinearRegression, LogisticRegression

rng = np.random.default_rng(0)
def draw(n, mean, sd):
    x = rng.normal(mean, sd, n)
    return x[:, None], np.sin(x) + 0.1 * rng.standard_normal(n)

Xs, ys = draw(2000, 0.0, 1.0)            # source: training field
Xt, yt = draw(2000, 1.5, 0.6)            # target: new field (covariate shift)

# domain classifier -> density ratio w(x) = p_t(x) / p_s(x)  (Proposition 9.3)
feats = lambda X: np.hstack([X, X ** 2])
dom = LogisticRegression().fit(feats(np.vstack([Xs, Xt])), np.r_[np.zeros(2000), np.ones(2000)])
p = dom.predict_proba(feats(Xs))[:, 1]
w = p / (1 - p)                            # equal sample sizes
print(f"effective sample size: {w.sum() ** 2 / (w ** 2).sum():.0f} of 2000")

for name, sw in [("unweighted", None), ("importance-weighted", w)]:
    m = LinearRegression().fit(Xs, ys, sample_weight=sw)     # misspecified model
    print(f"{name:20s} target MSE = {np.mean((m.predict(Xt) - yt) ** 2):.3f}")

# split conformal prediction (Theorem 9.5), alpha = 0.1
m = LinearRegression().fit(Xs[:1000], ys[:1000])
scores = np.abs(ys[1000:] - m.predict(Xs[1000:]))            # calibration residuals
q = np.sort(scores)[int(np.ceil(0.9 * 1001)) - 1]
Xs2, ys2 = draw(2000, 0.0, 1.0)
cover = lambda X, y: np.mean(np.abs(y - m.predict(X)) <= q)
print(f"coverage, same distribution: {cover(Xs2, ys2):.3f}")
print(f"coverage, shifted field    : {cover(Xt, yt):.3f}")

10. Research directions and open problems

This chapter is for self-study. It points to active research lines that extend the course, with surveys and founding papers to start from.

Surveys. Bergen et al. (2019) review machine learning for data-driven discovery in solid Earth geoscience, and Reichstein et al. (2019) argue for combining deep learning with process understanding in Earth system science. Yu and Ma (2021) survey deep learning for geophysics, and Mousavi and Beroza (2022) do the same for seismology. Karpatne et al. (2017) set out the programme of theory-guided data science, close to the “integration with physical models” of this course.

Physics-informed networks and neural operators. Physics-informed neural networks (Raissi, Perdikaris and Karniadakis, 2019) add the residual of a differential equation to the training loss. Neural operators learn the solution map of a family of equations: the Fourier neural operator (Li et al., 2021) and DeepONet (Lu et al., 2021). As fast surrogates of forward models, they can make nonlinear inversion (ch. 8.5) and uncertainty quantification affordable. Their accuracy outside the training family remains an open question.

Generative priors and posterior sampling. Score-based and diffusion models trained on images can serve as priors for inverse problems, and their samples approximate the posterior instead of its mean (Song et al., 2022; Chung et al., 2023). This answers the blur of Theorem 3.3 and the hallucination risk of ch. 3.3 at once, provided the samples are calibrated. Checking that calibration on field data, where the truth is unknown, is open.

Full-waveform inversion. Full-waveform inversion fits entire seismic waveforms by solving the wave equation, and it suffers from local minima when the starting model is too far from the truth (Virieux and Operto, 2009). Adler, Araya-Polo and Poggio (2021) review how deep learning accelerates or replaces parts of seismic inversion workflows.

Data-driven inverse problems. Arridge et al. (2019) give a mathematical overview of learned regularisers, unrolled schemes and their theory. Guarantees of stability and data consistency for learned reconstructions are still partial.

Open problems close to this course:

• generalisation to a new basin or tool with an uncertainty that remains valid under shift;
• detecting hallucinated structures in super-resolved or inpainted images without ground truth;
• learning from interpreter disagreement rather than averaging it away;
• benchmarks with blind wells and fields, so that published scores measure R_{\mathrm{new}} (Definition 7.1);
• hybrid physics–learning models whose outputs provably fit the data at the noise level;
• self-supervised pre-training on large archives of logs, images and seismic, and its evaluation across fields.

Exercises

1. ★ A thin coal bed on a density log. A density tool with a 40 cm box response (Proposition 1.3) logs a 10 cm coal bed of density 1.40 g/cm³ inside a sandstone of density 2.45 g/cm³. What is the lowest density the log can show? How thick must a bed be for the log to show at least 90 % of the true contrast?

   Solution: \Delta = -1.05 g/cm³, and in absolute value the apparent contrast is at most \frac{10}{40} \times 1.05 = 0.26 g/cm³, so the log reads about 2.19 g/cm³ at best.

   With a box response the apparent contrast is \frac{h}{w}\Delta for h < w and \Delta for h \ge w; 90 % therefore requires h \ge 0.9\, w = 36 cm.

2. ★ Dip from a sinusoid. On the image of an 8.5-inch borehole (r = 10.8 cm), a fracture reaches its shallowest point at azimuth 200° and depth 1012.40 m, and its deepest point at azimuth 20° and depth 1012.90 m. Give its dip, dip azimuth and the parameter z_0 of Theorem 1.5.

   Solution: the amplitude is half the span, r\tan\theta = 0.25 m, so \tan\theta = 2.315 and \theta = 66.6°. The dip azimuth is the azimuth of the deepest point, \varphi_d = 20°, and z_0 = 1012.65 m is the mid-depth.

3. ★ PSNR of a uniform quantiser. An image with values in [0, 1] is quantised with b bits, that is with step \Delta = 2^{-b}. Using Proposition 2.5, show that \mathrm{PSNR} \approx 6.02\, b + 10.79 \text{ dB}. Check the formula against the 8-bit line of code 2.8.

   Solution: with \mathrm{MSE} = \Delta^2/12 and peak 1, \mathrm{PSNR} = 10\log_{10}(12 \cdot 4^b) = 10\log_{10}12 + 20\, b \log_{10}2 = 10.79 + 6.02\, b . For b = 8 this gives 58.96 dB, and the measured MSE 1.27 \times 10^{-6} gives 58.96 dB.

4. ★ Reflectivity and back. Three layers have impedances 6.0, 7.5 and 5.0 (in 10^6 SI units). Compute the two reflection coefficients, then recover the impedance ratios with Proposition 8.5, and compare with the small-contrast rule \Delta \ln Z \approx 2R.

   Solution: R_1 = 1.5/13.5 = 0.111 and R_2 = -2.5/12.5 = -0.2. Then 2\operatorname{artanh}(0.111) = 0.2231 = \ln(7.5/6) and 2\operatorname{artanh}(-0.2) = -0.4055 = \ln(5/7.5), exactly. The approximation gives 0.222 and −0.400, errors of 0.4 % and 1.3 %.

5. ★ Dynamic time warping by hand. Align a = (1, 3, 2) with b = (1, 1, 3, 2) with the squared difference (Theorem 5.4). Give the table, the optimal cost and an optimal path. Why is the cost zero although the sequences have different lengths?

   Solution: the table is \begin{pmatrix} 0 & 0 & 4 & 5 \\ 4 & 4 & 0 & 1 \\ 5 & 5 & 1 & 0 \end{pmatrix}, so the cost is 0, with path (1,1), (1,2), (2,3), (3,4). The sequence b is a with its first sample repeated, a pure stretch, which warping absorbs at no cost.

6. ★ Ranking metrics. (a) Two rankings of five items are (A, B, C, D, E) and (B, A, C, E, D); compute Kendall’s tau. (b) A retrieval system returns four intervals with relevances (2, 3, 0, 1); compute DCG@4 and NDCG@4.

   Solution: (a) two discordant pairs, (A, B) and (D, E), out of 10, so \tau = (8 - 2)/10 = 0.6. (b) \mathrm{DCG} = 3 + 7/\log_2 3 + 0 + 1/\log_2 5 = 7.85; the ideal order (3, 2, 1, 0) gives 9.39, so NDCG = 0.835.

7. ★★ PCA with one component. For zero-mean x with covariance C, show that \max_{\|q\| = 1} q^\top C q = \lambda_1 and that \mathbb{E}\|x - qq^\top x\|^2 = \operatorname{tr} C - q^\top C q . Apply to the covariance with rows (4, 2) and (2, 1).

   Solution: in the eigenbasis, q^\top C q = \sum_i \lambda_i c_i^2 with \sum_i c_i^2 = 1, at most \lambda_1; the error formula is step 1 of the proof of Theorem 2.4 with k = 1. Here the eigenvalues are 5 and 0, with first direction (2, 1)/\sqrt 5: the data lie on a line and one component codes them without error.

8. ★★ Three estimators for one ambiguous sample. Given the coarse data, a fine sample equals 0 with probability 0.7 and 1 with probability 0.3. Compute the MMSE, median and MAP estimates, and the expected squared and absolute errors of the first two.

   Solution: the mean is 0.3; the median and the mode are 0. For c = 0.3: squared error 0.21, absolute error 0.42. For c = 0: squared error 0.30, absolute error 0.30. Each estimate wins on its own loss (Theorem 3.3, Proposition 3.5), and the MMSE value 0.3 never occurs in reality.

9. ★★ Harmonic inpainting on a line. Values u_0 = a and u_{m+1} = b are known, and u_1, \dots, u_m minimise \sum_{i=0}^{m} (u_{i+1} - u_i)^2. Show that the minimiser is the linear interpolation between a and b, and deduce the maximum principle in this case. What does this predict for a thin layer crossing a wide gap?

   Solution: the optimality condition u_i = \frac12(u_{i-1} + u_{i+1}) says the second difference vanishes, so u_i = a + (b - a)\frac{i}{m+1}, which lies between a and b. A thin layer seen on one side only is not continued: the fill joins the two edges linearly, which motivates the geometric priors of ch. 6.2.

10. ★★ Dice, IoU and empty masks. For |A| = 30, |B| = 40 and |A \cap B| = 20, compute IoU and Dice and check Proposition 7.6. What happens when both masks are empty, and why does the soft Dice loss of code 7.8 add \varepsilon?

    Solution: IoU = 20/50 = 0.4 and Dice = 40/70 = 0.571 = 0.8/1.4. With two empty masks both ratios are 0/0; a convention (score 1) is needed. The \varepsilon makes the loss defined and equal to 0 in that case, and smooths the gradient on images with few positive pixels.

11. ★★ The best Tikhonov parameter in one dimension. Let y = s\, x^\star + \varepsilon with \varepsilon \sim \mathcal{N}(0, \sigma^2), and x_\lambda = s y / (s^2 + \lambda). Using Theorem 8.2, find the \lambda that minimises the mean squared error, and compare with Proposition 0.6.

    Solution: with f = s^2/(s^2 + \lambda), the error is (1 - f)^2 x^{\star 2} + f^2 \sigma^2 / s^2, minimal at f = s^2 x^{\star 2} / (s^2 x^{\star 2} + \sigma^2). Hence \lambda^\star = \frac{\sigma^2}{x^{\star 2}}, the MAP weight of Proposition 0.6 with prior variance \tau^2 = x^{\star 2}.

12. ★★ Well-wise validation, computed. Starting from code 7.3: (a) vary the standard deviation of the well offsets from 0 to 1.2 and plot random and well-wise accuracies; (b) make well 0 contain 80 % of facies 0 and measure the effect of per-well standardisation on that well; (c) propose a normalisation that survives (b).

    Hints and expected result: (a) the two curves coincide near 0 and separate as the offset grows. (b) Standardisation subtracts a mean dominated by facies 0, so that well’s accuracy drops. (c) Normalise on a reference interval or on a facies that every well contains (for example a regional shale), as interpreters do.

13. ★★ How much data do importance weights leave? Let the source be \mathcal{N}(0, 1) and the target \mathcal{N}(\mu, 1). Compute the weight w and its first two moments under the source, and deduce the limit of the effective sample size divided by n. Evaluate for \mu = 1.5 and n = 2000.

    Solution: w(x) = e^{\mu x - \mu^2/2}, \mathbb{E}_s[w] = 1 and \mathbb{E}_s[w^2] = e^{\mu^2}. So \frac{n_{\mathrm{eff}}}{n} \to e^{-\mu^2}, about 211 samples out of 2000 for \mu = 1.5, the same order as the 347 of code 9.7 (whose target is narrower).

14. ★★ Conformal coverage is also bounded above. With continuous scores (no ties), show that \mathbb{P}(|Y - f(X)| \le \hat q) \le 1 - \alpha + \frac{1}{n + 1}. Is the coverage of 0.902 printed by code 9.7 surprising?

    Solution: the event is exactly “rank \le k”, of probability k/(n+1), and k \le (1 - \alpha)(n + 1) + 1. For n = 1000 and \alpha = 0.1 the marginal coverage lies in [0.9, 0.901].

    With one fixed calibration set, the realised coverage fluctuates with a standard deviation near \sqrt{\alpha(1-\alpha)/n} \approx 0.009, so 0.902 is expected.

15. ★★★ Losses for super-resolution. From code 3.7, train the CNN with the squared loss and with the absolute loss. Compare the RMSE, the sharpness at bed boundaries (mean absolute jump over the true boundaries) and the data residual \|A\hat x - y\| relative to the noise level. Then add a small adversarial term and discuss the perception–distortion trade-off.

    Hints: expect similar RMSE, sharper boundaries with the absolute loss (Proposition 3.5), and sharper but less faithful outputs with the adversarial term. A data residual far above \sigma\sqrt{m} reveals invented structure.

16. ★★★ Physics-guided gap filling. Combine Theorem 1.5 with ch. 6.2: on the image of code 6.7 with added fractures, detect sinusoids on the measured pads (least squares on picked points, or a Hough transform), extend them across the gaps, and inpaint only the residual harmonically. Evaluate with hold-out masks.

    Hints: the target is an RMSE in the gaps below the 0.139 of harmonic filling, and a dip error below 2° for fractures that cross a gap. Discuss the failure cases: several crossing fractures, non-planar features, and breakouts.

Projects and assessment

The official assessment is “Évaluation continue non intégrale” (partial continuous assessment) on the Université Paris-Saclay course detail, with no weights displayed; master-mva.com displays no assessment modality. The proposals below are not official. They are sized for a 30 h module with 18 h of labs: solo and group projects in the MVA tradition of reproducing a paper, and a written exam that shows what a final test on these notes could look like.

Solo projects

1. [Solo · ≈ 20 h] Paper study — learned ISTA for seismic deconvolution. Does an unrolled network reproduce on seismic traces the speed-up reported by Gregor and LeCun (2010)?
   • Work: implement ISTA, FISTA and LISTA in PyTorch; train LISTA on synthetic traces generated as in code 8.7, with Ricker frequencies drawn between 20 and 35 Hz.
   • Success criterion: with 10 layers, LISTA comes within 10 % of the test reflectivity error that ISTA reaches after 500 iterations.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters: ch. 0, 1.4, 8.
   • Evaluation: the rubric below; comparisons must be made at equal compute.
2. [Solo · ≈ 20 h] Facies classification with honest validation. How much of a published facies score survives a well-wise protocol?
   • Data: the public wells of Hall (2016), or the public data of the FORCE 2020 lithology competition (wells of the Norwegian continental shelf).
   • Work: logistic regression, random forest and gradient boosting; random versus well-wise cross-validation; per-well normalisation; HMM smoothing along depth; prior correction (Proposition 7.4).
   • Success criterion: a documented gap between the two protocols, and a model that beats logistic regression in macro F1 under the well-wise protocol.
   • Deliverables: code, report, talk. Chapters: ch. 4.3, 7.1, 9.2. Evaluation: rubric, with the protocol criterion weighted as stated.
3. [Solo · ≈ 20 h] Paper study — Deep Image Prior for pad gaps. Can the untrained prior of Ulyanov, Vedaldi and Lempitsky (2018) fill borehole-image gaps better than harmonic filling?
   • Data: the generator of code 6.7, extended with sinusoidal fractures (code 7.8).
   • Work: harmonic filling, sinusoid extension (Exercise 16) and Deep Image Prior, evaluated with hold-out masks.
   • Success criterion: a gap RMSE below the harmonic baseline, and a dip error below 2° for fractures picked across the filled gaps.
   • Deliverables: code, report, talk. Chapters: ch. 1.3, 6, 4.1. Evaluation: rubric.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Seismic facies segmentation on the F3 benchmark. Reproduce and extend the baseline of Alaudah et al. (2019) on the public F3 block of the Dutch North Sea.
   • Work: a U-Net on inline and crossline sections with the benchmark’s split; augmentations that respect the physics; deep ensembles and calibration of pixel probabilities.
   • Split: member A builds the data pipeline and augmentations, member B the model and training, member C the evaluation, uncertainty and error analysis.
   • Success criterion: class-wise IoU on the official test split, compared with the paper’s baseline, plus one ablation.
   • Deliverables: shared repository, 6–8 page report, 15-minute talk. Chapters: ch. 7.2, 4.1, 9.3. A free GPU tier is enough.
2. [Group of 2–3 · ≈ 30 h each] Paper study — faults learned from synthetic data. Reproduce in 2D the idea of Wu et al. (2019): train on synthetic sections only, test on field data.
   • Work: a generator of folded reflectivity models with fault displacements, convolved with Ricker wavelets (ch. 1.4); a U-Net with cross-entropy and soft Dice; tests on held-out synthetic sections and on public F3 sections.
   • Split: member A writes the generator, member B the network, member C the study of the synthetic-to-field gap (noise models, normalisation, ch. 9.2).
   • Success criterion: Dice ≥ 0.8 on held-out synthetic sections, and a documented list of failure modes on field sections.
   • Deliverables: repository, report, talk. Chapters: ch. 1.4, 7.2, 9.1–9.2.
3. [Group of 2–3 · ≈ 30 h each] Compressing borehole images for a telemetry budget. Which codec preserves interpretation best at 0.1 to 1 bit per pixel?
   • Data: synthetic unrolled borehole images from the generators of codes 2.8 and 7.8, a few thousand 256 × 256 images whose sinusoid parameters are known, so that dip errors are exact.
   • Work: block DCT with uniform quantisation and entropy coding; PCA on patches; a small learned codec trained on the rate–distortion loss of Ballé, Laparra and Simoncelli (2017).
   • Split: member A builds the classical codecs and the entropy coder, member B the learned codec, member C the task metric (dip error of automatically picked sinusoids) and the ranking of methods (ch. 4).
   • Success criterion: rate–distortion curves in PSNR, SSIM and dip error for the three codecs, and the lowest rate keeping the dip error below 2°.
   • Deliverables: repository, report, talk. Chapters: ch. 2, 4, 9.4.

Written exam (proposal)

Two hours, one handwritten A4 sheet allowed, a basic calculator, four exercises for 20 points.

1. Measurement and geometry (4 points). (a) State the convolutional model of a seismic trace and explain tuning in two sentences. (b) On an 8.5-inch borehole image (r = 10.8 cm), a sinusoid spans 37.4 cm from top to bottom, with its deepest point at azimuth 135°. Give the dip and the dip azimuth. Expected: (a) Definition 1.7, then Proposition 1.8 and the quarter-wavelength interference. (b) Amplitude 18.7 cm, \tan\theta = 1.73, so \theta = 60°, dip azimuth 135°.
2. Tikhonov by hand (6 points). Let A = \operatorname{diag}(2, 0.1), x^\star = (1, 1) and y = (2.1, 0.3).
   • Part (a): Compute A^{-1}y and x_\lambda for \lambda = 0.01.
   • Part (b): With noise standard deviation \sigma = 0.2, compute the expected squared error of both estimators with Theorem 8.2.
   • Part (c): Which \lambda is best for the second component alone?
   • Expected: (a) (1.05, 3.0) and (1.047, 1.5). (b) 4.01 against 1.26. (c) \lambda = \sigma^2 / x_2^{\star 2} = 0.04 (Exercise 11), giving 0.80 on that component.
3. A short proof (5 points). Prove Theorem 6.2. Deduce a transform that gives an image the histogram of a target distribution function G. Expected: the proof of Theorem 6.2 with the generalised inverse; then G^{-} \circ F, because G^{-}(V) \sim G for V uniform (inverse transform sampling).
4. Validation under shift (5 points). A facies classifier scores 96 % in random 5-fold cross-validation on six wells and 60 % on a seventh. (a) Explain the gap. (b) Propose a protocol that predicts the 60 %. (c) Propose two remedies and how to test them. (d) How would you attach a 90 % interval to a regression output, and when does the guarantee fail? Expected: dependence within wells and well-to-well shift (Definition 7.1); well-wise cross-validation; per-well or reference-interval normalisation and importance weighting, tested well-wise; split conformal prediction (Theorem 9.5), which fails without exchangeability.

Grading rubric (proposal)

References

Official references. Neither master-mva.com nor the Université Paris-Saclay course detail lists references for this course. The works below are standard sources for its topics, grouped by chapter.

Measurement systems and geophysics (ch. 1, 5, 8)

• O. Serra, Fundamentals of Well-Log Interpretation, Vol. 1: The Acquisition of Logging Data, Developments in Petroleum Science 15A, Elsevier, 1984.
• D. V. Ellis, J. M. Singer, Well Logging for Earth Scientists, 2nd ed., Springer, 2007.
• R. E. Sheriff, L. P. Geldart, Exploration Seismology, 2nd ed., Cambridge University Press, 1995.
• Ö. Yilmaz, Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data, Society of Exploration Geophysicists, 2001.
• G. E. Archie, “The electrical resistivity log as an aid in determining some reservoir characteristics”, Transactions of the AIME 146, 1942.
• M. B. Widess, “How thin is a thin bed?”, Geophysics 38(6), 1973.
• M. R. J. Wyllie, A. R. Gregory, L. W. Gardner, “Elastic wave velocities in heterogeneous and porous media”, Geophysics 21(1), 1956.

Signals, compression and images (ch. 0, 2, 4, 6)

• S. Mallat, A Wavelet Tour of Signal Processing: The Sparse Way, 3rd ed., Academic Press, 2009.
• T. M. Cover, J. A. Thomas, Elements of Information Theory, 2nd ed., Wiley, 2006.
• R. C. Gonzalez, R. E. Woods, Digital Image Processing, 4th ed., Pearson, 2018.
• N. Ahmed, T. Natarajan, K. R. Rao, “Discrete cosine transform”, IEEE Transactions on Computers C-23(1), 1974.
• P. Baldi, K. Hornik, “Neural networks and principal component analysis: learning from examples without local minima”, Neural Networks 2(1), 1989.
• J. Ballé, V. Laparra, E. P. Simoncelli, “End-to-end optimized image compression”, ICLR, 2017.
• Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity”, IEEE Transactions on Image Processing 13(4), 2004.
• D. L. Donoho, I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage”, Biometrika 81(3), 1994.
• K. Zuiderveld, “Contrast limited adaptive histogram equalization”, in P. Heckbert (ed.), Graphics Gems IV, Academic Press, 1994.
• J. Shen, T. F. Chan, “Mathematical models for local nontexture inpaintings”, SIAM Journal on Applied Mathematics 62(3), 2002.
• A. Criminisi, P. Pérez, K. Toyama, “Region filling and object removal by exemplar-based image inpainting”, IEEE Transactions on Image Processing 13(9), 2004.
• N. F. Hurley, T. Zhang, “Method to generate full-bore images using borehole images and multipoint statistics”, SPE Reservoir Evaluation & Engineering 14(2), 2011.
• G. Liu, F. A. Reda, K. J. Shih, T.-C. Wang, A. Tao, B. Catanzaro, “Image inpainting for irregular holes using partial convolutions”, ECCV, 2018.
• D. Ulyanov, A. Vedaldi, V. Lempitsky, “Deep image prior”, CVPR, 2018.

Super-resolution (ch. 3)

• C. Dong, C. C. Loy, K. He, X. Tang, “Image super-resolution using deep convolutional networks”, IEEE Transactions on Pattern Analysis and Machine Intelligence 38(2), 2016.
• J. Kim, J. K. Lee, K. M. Lee, “Accurate image super-resolution using very deep convolutional networks”, CVPR, 2016.
• W. Shi, J. Caballero, F. Huszár, J. Totz, A. P. Aitken, R. Bishop, D. Rueckert, Z. Wang, “Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network”, CVPR, 2016.
• C. Ledig et al., “Photo-realistic single image super-resolution using a generative adversarial network”, CVPR, 2017.
• Y. Blau, T. Michaeli, “The perception-distortion tradeoff”, CVPR, 2018.

Ranking and alignment (ch. 4, 5)

• R. A. Bradley, M. E. Terry, “Rank analysis of incomplete block designs: I. The method of paired comparisons”, Biometrika 39(3/4), 1952.
• L. R. Ford Jr., “Solution of a ranking problem from binary comparisons”, The American Mathematical Monthly 64(8), 1957.
• C. Burges, T. Shaked, E. Renshaw, A. Lazier, M. Deeds, N. Hamilton, G. Hullender, “Learning to rank using gradient descent”, ICML, 2005.
• K. Järvelin, J. Kekäläinen, “Cumulated gain-based evaluation of IR techniques”, ACM Transactions on Information Systems 20(4), 2002.
• H. Sakoe, S. Chiba, “Dynamic programming algorithm optimization for spoken word recognition”, IEEE Transactions on Acoustics, Speech, and Signal Processing 26(1), 1978.
• M. Cuturi, M. Blondel, “Soft-DTW: a differentiable loss function for time-series”, ICML, 2017.
• W. Kabsch, “A solution for the best rotation to relate two sets of vectors”, Acta Crystallographica A 32(5), 1976.
• M. Jaderberg, K. Simonyan, A. Zisserman, K. Kavukcuoglu, “Spatial transformer networks”, NeurIPS, 2015.
• G. Balakrishnan, A. Zhao, M. R. Sabuncu, J. Guttag, A. V. Dalca, “VoxelMorph: a learning framework for deformable medical image registration”, IEEE Transactions on Medical Imaging 38(8), 2019.

Classification and segmentation (ch. 7)

• B. Hall, “Facies classification using machine learning”, The Leading Edge 35(10), 2016.
• Y. Alaudah, P. Michałowicz, M. Alfarraj, G. AlRegib, “A machine-learning benchmark for facies classification”, Interpretation 7(3), 2019.
• X. Wu, L. Liang, Y. Shi, S. Fomel, “FaultSeg3D: using synthetic data sets to train an end-to-end convolutional neural network for 3D seismic fault segmentation”, Geophysics 84(3), 2019.
• J. Long, E. Shelhamer, T. Darrell, “Fully convolutional networks for semantic segmentation”, CVPR, 2015.
• O. Ronneberger, P. Fischer, T. Brox, “U-Net: convolutional networks for biomedical image segmentation”, MICCAI, 2015.
• T. Chen, C. Guestrin, “XGBoost: a scalable tree boosting system”, KDD, 2016.
• M. Saerens, P. Latinne, C. Decaestecker, “Adjusting the outputs of a classifier to new a priori probabilities: a simple procedure”, Neural Computation 14(1), 2002.
• T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning, 2nd ed., Springer, 2009.
• I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016.

Inverse problems (ch. 8)

• P. C. Hansen, Discrete Inverse Problems: Insight and Algorithms, SIAM, 2010.
• A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, 2005.
• G. H. Golub, M. Heath, G. Wahba, “Generalized cross-validation as a method for choosing a good ridge parameter”, Technometrics 21(2), 1979.
• I. Daubechies, M. Defrise, C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint”, Communications on Pure and Applied Mathematics 57(11), 2004.
• A. Beck, M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems”, SIAM Journal on Imaging Sciences 2(1), 2009.
• K. Gregor, Y. LeCun, “Learning fast approximations of sparse coding”, ICML, 2010.
• S. V. Venkatakrishnan, C. A. Bouman, B. Wohlberg, “Plug-and-play priors for model based reconstruction”, IEEE GlobalSIP, 2013.
• S. Arridge, P. Maass, O. Öktem, C.-B. Schönlieb, “Solving inverse problems using data-driven models”, Acta Numerica 28, 2019.

Industrial challenges (ch. 9)

• H. Shimodaira, “Improving predictive inference under covariate shift by weighting the log-likelihood function”, Journal of Statistical Planning and Inference 90(2), 2000.
• J. Quiñonero-Candela, M. Sugiyama, A. Schwaighofer, N. D. Lawrence (eds.), Dataset Shift in Machine Learning, MIT Press, 2009.
• Y. Ganin, E. Ustinova, H. Ajakan, P. Germain, H. Larochelle, F. Laviolette, M. Marchand, V. Lempitsky, “Domain-adversarial training of neural networks”, Journal of Machine Learning Research 17, 2016.
• N. Natarajan, I. S. Dhillon, P. Ravikumar, A. Tewari, “Learning with noisy labels”, NeurIPS, 2013.
• V. Vovk, A. Gammerman, G. Shafer, Algorithmic Learning in a Random World, Springer, 2005.
• A. N. Angelopoulos, S. Bates, “Conformal prediction: a gentle introduction”, Foundations and Trends in Machine Learning 16(4), 2023.
• R. J. Tibshirani, R. Foygel Barber, E. J. Candès, A. Ramdas, “Conformal prediction under covariate shift”, NeurIPS, 2019.
• B. Lakshminarayanan, A. Pritzel, C. Blundell, “Simple and scalable predictive uncertainty estimation using deep ensembles”, NeurIPS, 2017.
• B. Jacob et al., “Quantization and training of neural networks for efficient integer-arithmetic-only inference”, CVPR, 2018.
• G. Hinton, O. Vinyals, J. Dean, “Distilling the knowledge in a neural network”, arXiv:1503.02531, 2015.
• M. Mitchell et al., “Model cards for model reporting”, ACM Conference on Fairness, Accountability, and Transparency (FAT*), 2019.

Research directions (ch. 10)

• K. J. Bergen, P. A. Johnson, M. V. de Hoop, G. C. Beroza, “Machine learning for data-driven discovery in solid Earth geoscience”, Science 363(6433), 2019.
• M. Reichstein, G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, Prabhat, “Deep learning and process understanding for data-driven Earth system science”, Nature 566, 2019.
• A. Karpatne et al., “Theory-guided data science: a new paradigm for scientific discovery from data”, IEEE Transactions on Knowledge and Data Engineering 29(10), 2017.
• S. Yu, J. Ma, “Deep learning for geophysics: current and future trends”, Reviews of Geophysics 59(3), 2021.
• S. M. Mousavi, G. C. Beroza, “Deep-learning seismology”, Science 377(6607), 2022.
• M. Raissi, P. Perdikaris, G. E. Karniadakis, “Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”, Journal of Computational Physics 378, 2019.
• Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar, “Fourier neural operator for parametric partial differential equations”, ICLR, 2021.
• L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis, “Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators”, Nature Machine Intelligence 3, 2021.
• Y. Song, L. Shen, L. Xing, S. Ermon, “Solving inverse problems in medical imaging with score-based generative models”, ICLR, 2022.
• H. Chung, J. Kim, M. T. McCann, M. L. Klasky, J. C. Ye, “Diffusion posterior sampling for general noisy inverse problems”, ICLR, 2023.
• J. Virieux, S. Operto, “An overview of full-waveform inversion in exploration geophysics”, Geophysics 74(6), 2009.
• A. Adler, M. Araya-Polo, T. Poggio, “Deep learning for seismic inverse problems: toward the acceleration of geophysical analysis workflows”, IEEE Signal Processing Magazine 38(2), 2021.