Evolution of Inversion Algorithms for Rayleigh Wave Dispersion: From Haskell-Thomson to Neighborhood Algorithms

Surface wave inversion is a fundamental procedure in geophysics and civil engineering used to characterize the elastic properties of the subsurface. The process involves measuring the dispersion characteristics of Rayleigh and Love waves—where different frequencies travel at different velocities—and applying mathematical models to reconstruct the vertical profile of shear-wave velocities (Vs), density, and layer thickness. This discipline has evolved through several distinct phases, moving from the foundational matrix calculations of the mid-20th century to the sophisticated stochastic optimization algorithms used in modern site characterization.

The efficacy of these algorithms directly impacts the accuracy of subsurface imaging, affecting everything from earthquake hazard assessment to the non-destructive testing of bridge decks and pavements. As computational power has increased, the focus of the scientific community has shifted from finding a single ‘best-fit’ solution toward understanding the non-uniqueness and uncertainty inherent in seismic data. This progression reflects a broader trend in computational physics: the move from deterministic, linearized approximations to global search techniques that explore a vast parameter space to identify high-probability geological models.

Timeline

• 1950:William Thomson publishes a matrix method for calculating the transmission of elastic waves through a stratified solid medium, establishing the ground for modern dispersion modeling.
• 1953:Norman A. Haskell refines Thomson’s work, introducing a systematic matrix formulation that becomes known as the Haskell-Thomson method, which allows for the calculation of dispersion in an arbitrary number of layers.
• 1964:The introduction of the ‘delta matrix’ formulation by Harkrider addresses numerical instabilities found in the original Haskell-Thomson approach at high frequencies.
• 1970s:The development of the first linearized inversion programs based on the work of Dorman and Ewing, allowing for the iterative adjustment of a starting model to fit observed data.
• 1990s:Researchers begin applying global optimization techniques, including Genetic Algorithms (GA) and Simulated Annealing (SA), to circumvent the local minima problems of linearized inversion.
• 1999:Malcolm Sambridge introduces the Neighborhood Algorithm (NA), a derivative-free global search method that uses Voronoi cells to sample the parameter space efficiently.
• 2000s–Present:Integration of Bayesian frameworks and multi-modal inversion (combining Rayleigh and Love waves) becomes the standard for high-resolution geotechnical site investigation.

Background

The study of surface wave dispersion is rooted in the physical reality that seismic waves of different wavelengths sample different depths of the Earth’s crust. In a homogeneous medium, surface waves are non-dispersive; however, in a vertically heterogeneous medium—such as a layered soil profile or a road structure—the velocity depends on the frequency. Rayleigh waves, which involve both longitudinal and transverse motion in a vertical plane, are particularly sensitive to the shear-wave velocity of the material, making them ideal for determining the stiffness of geological layers.

The mathematical challenge lies in the ‘inverse problem.’ While the ‘forward problem’ (calculating the dispersion curve from a known soil profile) is relatively straightforward using matrix methods, the inverse problem (determining the soil profile from an observed dispersion curve) is non-linear and ill-posed. This means that multiple different soil configurations could potentially produce the same observed dispersion curve. Consequently, the history of inversion algorithms is essentially a quest for methods that can reliably handle this ambiguity while remaining computationally feasible.

The Haskell-Thomson Matrix Era

The modern era of surface wave analysis began with the landmark papers by Thomson (1950) and Haskell (1953). Before their work, solving the wave equation for more than two or three layers was mathematically cumbersome and prone to error. The Haskell-Thomson method introduced a transfer matrix approach where each layer of a stratigraphic column is represented by a matrix. By multiplying these matrices together, the boundary conditions at the surface and the deep half-space could be satisfied for a system of any number of layers.

Despite its brilliance, the original formulation suffered from ‘precision loss’ when dealing with high frequencies or thick layers. As frequencies increased, the terms within the matrices grew exponentially, leading to numerical instability in computers of that era. This necessitated the development of the ‘reduced delta matrix’ method and other variants like the Schwab-Knopoff method in the 1960s and 70s. These refinements ensured that the forward calculation was strong enough to be used as the engine for automated inversion routines.

Transition to Linearized Inversion

By the 1970s, the focus shifted toward automating the search for the best-fit model. The first generation of automated tools utilized linearized inversion. These methods require an initial ‘guess’ or starting model. The algorithm then calculates the derivatives (the Jacobian matrix) of the phase velocity with respect to the model parameters (velocity, thickness, density) and uses these derivatives to iteratively nudge the model toward a better fit.

While computationally fast, linearized methods possess a significant flaw: they are highly dependent on the initial model. If the starting guess is too far from the actual geological reality, the algorithm often converges to a ‘local minimum’—a mathematical solution that looks good locally but is globally incorrect. In the context of engineering, this could lead to an overestimation or underestimation of soil stiffness, with potentially dangerous implications for infrastructure design.

The Rise of Global Optimization in the 1990s

To overcome the limitations of linearized methods, the 1990s saw the introduction of global search algorithms. These methods do not rely on derivatives and do not require a precise starting model. Instead, they explore the entire range of physically plausible parameters using stochastic (randomized) techniques. Two primary methods dominated this era: Genetic Algorithms (GA) and Simulated Annealing (SA).

• Genetic Algorithms:Inspired by biological evolution, GA maintains a ‘population’ of models. The best-performing models are ‘bred’ together, using crossover and mutation operators to create a new generation of models that increasingly fit the observed data.
• Simulated Annealing:Based on the cooling process of metals, SA allows the algorithm to occasionally accept ‘worse’ models early in the search to escape local minima, gradually narrowing its focus as the ‘temperature’ of the simulation drops.

These global methods provided a more objective way to handle complex layered systems where the depth to bedrock or the presence of velocity inversions (a soft layer trapped between two hard layers) made linearized methods unreliable.

The Neighborhood Algorithm and Modern Benchmarking

In 1999, Malcolm Sambridge introduced the Neighborhood Algorithm (NA), which represented a significant leap in the efficiency of global searches. Unlike GA or SA, which often require thousands of forward model evaluations to reach convergence, NA uses the geometry of the previously sampled points to guide the search. It partitions the parameter space into Voronoi cells—regions where every point is closer to its generating node than to any other node.

The Neighborhood Algorithm prioritizes sampling in cells where the data fit is highest, but it maintains a level of exploration in less-fit regions to ensure global coverage. This dual focus allows it to map out the ‘uniqueness’ of the solution more clearly than its predecessors.

Comparative performance reviews using synthetic datasets—where the ‘true’ subsurface model is known—have consistently shown that the Neighborhood Algorithm is superior in identifying multiple high-probability solutions. In a typical benchmark, NA can identify the correct velocity profile in fewer iterations than a Genetic Algorithm while providing a more strong estimate of the uncertainty associated with each layer thickness.

Comparative Methodology Review

Contemporary Applications and Surface Wave Hub Focus

Today, the practical application of these algorithms at the Surface Wave Hub involves the integration of passive and active source data. Passive methods, such as the Analysis of Microtremors (MAM), use the background vibrations of the Earth and human activity. Active methods, like the Spectral Analysis of Surface Waves (SASW) or Multichannel Analysis of Surface Waves (MASW), use controlled sources like sledgehammers or weight drops. The inversion of this combined data requires algorithms that can handle broad frequency ranges and varying signal-to-noise ratios.

Modern research is also pushing into the area of 2D and 3D inversion. While the Haskell-Thomson foundation assumed perfectly horizontal layers, current inversion algorithms are being adapted to handle lateral variations in soil properties. This is critical for detecting localized anomalies such as buried utilities, voids, or the edge of a landslide slip plane. By refining the inversion of Rayleigh and Love wave dispersion, researchers continue to provide clearer, more reliable windows into the hidden structures beneath the surface.