Seismic traveltime inversion with quantum annealing (2025)

Quantum annealing

Quantum computing is rapidly emerging as a pivotal area of scientific and technological advancement, attracting considerable investment and interest due to its profound potential^23,24,25. Unlike classical computers that use bits, which exist only in states of 0 or 1, quantum computers employ quantum bits, or qubits. Qubits possess unique properties such as superposition, entanglement, and interference^26,27,28, enabling them to perform certain complex computations beyond the capabilities of classical computers^29,30. Qubits can be constructed from various physical systems such as photons, trapped atoms, nuclear magnetic resonance, quantum dots, dopants in solids, and superconductors³¹. Previous research^{32,33,34,35,36} has provided evidence that quantum computing possibly surpasses classical computers in terms of processing speed and efficiency for certain problems.

The quantum annealing process facilitates the finding of the global minimum of a cost function efficiently. This process can be described using the real-time Schrödinger equation³⁷:

$$\begin{aligned} i \hbar \frac{d}{dt} |\Psi (t)\rangle = H(t)|\Psi (t)\rangle \end{aligned}$$

(1)

where $| \cdot \rangle$ is the ket of the Dirac notation³⁸, $i$ is the imaginary unit, $t$ is time, $\hbar$ is the reduced Planck’s constant, $\Psi (t)$ is the wave function, $|\Psi (t)\rangle$ is the quantum state vector, $H$ is the Hamiltonian representing the total energy of the quantum system^39,40. If $\hbar$ is set as 1, the Eq. 1 becomes:

$$\begin{aligned} i \frac{d}{dt} |\Psi (t)\rangle = H(t)|\Psi (t)\rangle \end{aligned}$$

(2)

The Hamiltonian in quantum annealing can be composed of two components^41,42:

$$\begin{aligned} H(t) = A(t) H_0 + B(t) H_1 \end{aligned}$$

(3)

where $H_0$ is the initial Hamiltonian, representing a system with an initial ground state. $H_1$ is the final Hamiltonian, whose ground state encodes the solution to the optimization problem. $A(t)$ and $B(t)$ are time-dependent coefficients. $A(t)$ and $B(t)$ are set in a range of 0 to 1 so that $A(t_0) \gg B(t_0)$ at the initial time $t_0$ and $B(t_1) \gg A(t_1)$ at the final time $t_0$. During the process, $A(t)$ monotonically decreases, while $B(t)$ monotonically increases. At the start of the annealing process, $H(t) \approx H_0$. At the end of the annealing process, $H(t) \approx H_1$. Thus, the system transitions from the ground state of $H_0$ to the ground state of $H_1$. If $H(t)$ changes sufficiently slowly, the state evolves adiabatically⁴³.

The problems are then often mapped onto a Quadratic Unconstrained Binary Optimization (QUBO) or Ising model⁴⁴:

$$\begin{aligned} \text {QUBO:} \quad \min _{x_j = 0,1} \left( \sum _{j \le k} x_j Q_{jk} x_k + C_1 \right) \end{aligned}$$

(4)

$$\begin{aligned} \text {Ising:} \quad \min _{s_j = \pm 1} \left( \sum _j h_j s_j + \sum _{j < k} J_{jk} s_j s_k + C_2 \right) \end{aligned}$$

(5)

where $j, k$ are indices, ranging over all qubits. In the QUBO model (Eq. 4), $Q$ is the QUBO matrix with values $Q_{jk} \in \mathbb {R}$. The binary variable vector is $\textbf{x}$ with $x_j \in \{0, 1\}$. In the Ising model (Eq. 5), the problem is defined by the biases $h_j \in \mathbb {R}$ and the couplers $J_{jk} \in \mathbb {R}$, and the binary variable vector is $\textbf{s}$ with $s_j \in \{-1, 1\}$. $C_1$ and $C_2$ are constants which do not affect the solution of the optimization problem. The Ising model and the QUBO model are mathematically equivalent, allowing them to be translated into each other. This equivalence provides a flexible approach to problem-solving, enabling the conversion of problems between these models based on the requirements and available tools. There are also tools, such as ToQUBO.jl, designed to convert standard problems into the QUBO format⁴⁵. In this paper, we utilize the quantum annealer from D-Wave Advantage Systems¹⁴ to employ direct quantum processing unit (QPU) methods for seismic travel inversion.

Seismic data acquisition

We construct a synthetic velocity model representing carbon storage applications, as shown in Fig. 4a. This model spans a depth range from 1000 to 1300 m and extends 100 m horizontally. The size of the grid cell for this model is 10 x 10 m. Within this model, the carbon storage structure is depicted as a wedge, starting from 1100 m and reaching a maximum depth of 1200 m. The velocity model is constructed with varying velocities to reflect real-world geological conditions. The average velocity within the carbon storage area ranges from 3180 to 3220 m/s, which is about 11% lower than the surrounding background velocity, which ranges from 3530 to 3640 m/s. Furthermore, the velocities increase with depth.

The carbon storage velocity model and ray coverage patterns. Red dots are sources, blue dots are receivers, and white lines represent the ray paths. (a) The synthetic velocity model with a wedge-shaped low-velocity carbon storage formation. (b) Ray coverage from sources and receivers placed in a uniform grid within two boreholes. (c) Ray coverage from sources and receivers placed in a non-uniform pattern, enhancing coverage and constraints for the quantum annealing process.

Full size image

The uniform placement (Fig. 4b) is commonly used in seismic data acquisition for simplicity in boreholes. However, this approach results in significantly lower ray coverage in the shallow and deep sections compared to the middle section. To address these limitations, in our survey, 20 pairs of sources and receivers are non-uniformly deployed within two boreholes (Fig. 4c). The non-uniform deployment is designed to introduce more constraints for the quantum annealing process, thereby improving the overall accuracy of the seismic inversion. The sources and receivers of non-uniform placement are distributed according to a quadratic polynomial distribution.

Transforming ray equations to QUBO

The set of ray equations can be represented as¹³:

$$\begin{aligned} \textbf{D} \textbf{s} = \textbf{T}, \end{aligned}$$

(6)

where $\textbf{D}$ is the matrix of distance increments $d_j$, $\textbf{s}$ is the slowness vector, and $\textbf{T}$ is the travel time vector. The size of $\textbf{D}$ can be very large, therefore solving for $\textbf{s}$ through matrix operations on $\textbf{D}$ can be computationally intensive. This challenge is exacerbated by the relatively sparse and random distribution of elements within $\textbf{D}$. Consequently, alternative methods are required to solve these problems efficiently while maintaining accuracy. Quantum annealers can provide quantum metaheuristic algorithms to address this issue. Eq. 6 can be solved by minimizing the objective function:

$$\begin{aligned} f(\textbf{s}) = \left\| \textbf{D} \textbf{s} - \textbf{T} \right\| ^2_2. \end{aligned}$$

(7)

The objective function $f(\textbf{s})$ computes the difference between the observed travel times and those predicted by the model given a slowness vector $\textbf{s}$. Minimizing this function ensures that the model’s predictions align as closely as possible with the observed data, thus achieving an optimal fit. Quantum annealers offer a direct approach to solving binary objective functions⁴⁶:

$$\begin{aligned} f(\textbf{q}) = \left\| \textbf{A}^d \textbf{q} - \textbf{b} \right\| ^2_2. \end{aligned}$$

(8)

In this formulation, $\textbf{q}$ is a binary vector, $\textbf{A}^d$ is a real-valued matrix, and $\textbf{b}$ is a real-valued vector. Because quantum computers are designed to solve QUBO problems, transforming real-valued variables to binary values is essential. However, the number of binary variables $n_{\text {binary}}$ increases with the number of bits R used for fixed-point approximation, and it is related to the number of real-valued variables $n_{\text {real}}$ as $n_{\text {binary}} = n_{\text {real}} \times R$. Higher values of R yield greater accuracy in representing the floating-point numbers, but the current limitations of quantum computer hardware restrict the number of qubits available. To address this issue without excessively increasing the number of binary variables, the initial velocity guess, variable boundaries¹⁵ and recursive methods⁴⁷ are employed. The initial guess and the boundaries are used to rescale the range of the slowness vector $\textbf{s}$ to a new vector $\textbf{x}$ such that $x_i \in [0, 2)$, facilitating easier binary representation. Recursive methods are then applied to enhance the precision of floating-point divisions. These methods iteratively refine the estimate of $\textbf{s}$, reducing the error at each step. The initial objective function Eq. 7 can be reformulated as:

$$\begin{aligned} f(\textbf{x}) = \left\| \textbf{D} \textbf{x} - \textbf{b} \right\| ^2_2, \end{aligned}$$

(9)

where $\textbf{b} = \left( \textbf{T} + L \textbf{I} - \textbf{D} \textbf{s}_0 \right) / L$, L is the variable boundaries, $\textbf{s}_0$ is the initial guess of the slowness vector, and $\textbf{I}$ is the identity vector. The slowness vector $\textbf{s}$ is then in the range of $[\textbf{s}_0 - L \textbf{I}; \textbf{s}_0 + L \textbf{I}]$. This range ensures that the solution space is adequately covered. To express this as a binary objective function, $x_i$ is represented in binary form using the R bit fixed-point approximation:

$$\begin{aligned} x_i = \sum _{r=0}^{R-1} 2^{-r} q_r, \end{aligned}$$

(10)

where $q_r \in \{0, 1\}$ is the value of the r-th bit. This transformation is essential for harnessing the computational power of quantum annealers, which are inherently designed to solve binary optimization problems. The new matrix $\textbf{A}^d$ in Eq. 8 is derived from $\textbf{D}$ and R such that $\textbf{D} \textbf{x} = \textbf{A}^d \textbf{q}$. The QUBO matrix $Q_{ij}$ in Eq. 4 is then constructed from the given matrix $\textbf{D}$ and the calculated vector $\textbf{b}$^46,48:

$$\begin{aligned} Q_{jj}&= \sum _i A_{ij} (A_{ij} - 2b_j) , \end{aligned}$$

(11)

$$\begin{aligned} Q_{jk}&= 2 \sum _i A_{ij} A_{ik} . \end{aligned}$$

(12)

The QUBO matrix is directly input into the Quantum Processing Unit (QPU). The system utilizes DWaveSampler() to employ a D-Wave system as the sampler. Subsequently, EmbeddingComposite() manages the mapping between the problem and the D-Wave system’s numerically indexed qubits, a process known as minor-embedding⁴⁹.

In this study, we perform traveltime inversion using $R=3$ for 10 iterations with quantum annealing. The total number of real-valued variables of the problem is 300. Due to quantum hardware limitations, we break down the model into 30 layers with 10 variables each. This division reduces the complexity of each sub-problem, making it manageable for the quantum processor and allowing for better control of the boundary $L$. Since the problem from Eq. 6 is depth-independent, we simplify the process by adjusting the system’s coordinates at each iteration. The approach ensures that each layer is treated independently, reducing the overall complexity of the inversion. By implementing these techniques, we can efficiently solve large-scale traveltime inversion problems using quantum annealers.

Tikhonov regularization linear least squares inversion

For ill-conditioned problems, small changes in $\textbf{D}$ or $\textbf{T}$ can lead to significant variations in the results⁵⁰. To mitigate the effects of noise in the data, we employ Tikhonov regularization methods as a benchmark for the comparision. The new objective function (Eq. 7) can be expressed as a general regularized form⁵¹:

$$\begin{aligned} f_\lambda (\textbf{s}) = \Vert \textbf{D} \textbf{s} - \textbf{T} \Vert _2^2 + \lambda g(\textbf{s}), \end{aligned}$$

(13)

where $\lambda$ is the regularization parameter controlling the trade-off between the data fidelity term $\Vert \textbf{D} \textbf{s} - \textbf{T} \Vert _2^2$ and the regularization term $g(\textbf{s})$. The Tikhonov regularization is flexible and allows different types of regularization functions. The standard Tikhonov regularization with $L_2$-norm is in the form:

$$\begin{aligned} f_\lambda (\textbf{s}) = \Vert \textbf{D} \textbf{s} - \textbf{T} \Vert _2^2 + \lambda \Vert \textbf{s}\Vert _2^2 \end{aligned}$$

where $\Vert \textbf{s}\Vert _2^2 = \textbf{s}^T \textbf{s}$ penalizes large values in the solution. Another form is the first-order Tikhonov regularization with a smoothness regularization:

$$\begin{aligned} f_\lambda (\textbf{s}) = \Vert \textbf{D} \textbf{s} - \textbf{T} \Vert _2^2 + \lambda \Vert D_1 \textbf{s} \Vert _2^2 \end{aligned}$$

where $D_1$ is the first-order difference operator which enforces smooth variation in $\textbf{s}$ by penalizing large first derivatives. Similarly, second-order Tikhonov regularization penalizes the curvature of the solution:

$$\begin{aligned} f_\lambda (\textbf{s}) = \Vert \textbf{D} \textbf{s} - \textbf{T} \Vert _2^2 + \lambda \Vert D_2 \textbf{s} \Vert _2^2 \end{aligned}$$

where $D_2$ is the second-order difference operator which enforces smooth curvature by penalizing large second derivatives. In general, $g(\textbf{s})$ can be $g(\textbf{s}) = \Vert \textbf{s}\Vert _2^2$ for standard $L_2$-norm regularization, $g(\textbf{s}) = \Vert D_1 \textbf{s} \Vert _2^2$ for first-order smoothness, or $g(\textbf{s}) = \Vert D_2 \textbf{s} \Vert _2^2$ for second-order smoothness. The choice of $g(\textbf{s})$ depends on prior knowledge and the desired properties of the solution. Here, we use a custom Tikhonov regularization $g(\textbf{s}) = \Vert \textbf{s} - \textbf{s}_0 \Vert _2^2$, where $\textbf{s}_0$ is the initial guess for the slowness and is chosen as the input for the quantum annealing process. The objective function is now expressed as:

$$\begin{aligned} f_\lambda (\textbf{s}) = \Vert \textbf{D} \textbf{s} - \textbf{T} \Vert _2^2 + \lambda \Vert \textbf{s} - \textbf{s}_0 \Vert _2^2. \end{aligned}$$

(14)

The solution to this regularized problem is given by:

$$\begin{aligned} \textbf{s} = \left( \textbf{D}^T \textbf{D} + \lambda \textbf{I} \right) ^{-1} \left( \textbf{D}^T \textbf{T} + \lambda \textbf{s}_0 \right) . \end{aligned}$$

(15)

Cost analysis

Quantum annealing directly solves $f(\textbf{q}) = \left\| \textbf{A}^d \textbf{q} - \textbf{b} \right\| ^2_2$ (Eq. 8) by providing the solution binary vector $\textbf{q}$. The quantum annealing process involves three sources of computational cost: preparing the binary problem, executing the annealing on the quantum hardware, and post-processing the results. The cost of preparing and post-processing is $\mathcal {O}(mn^2 + mnc + n^2c^2)$. For a single iteration of the loop, the cost of executing the annealing process on the quantum hardware is $\mathcal {O}(\mathcal {P}(cn))$. Here, $m$ and $n$ represent the rows and columns of the matrix $\textbf{D}$, respectively. $\mathcal {P}(cn)$ is a polynomial term. The parameter $c = |\Theta | + 1$, where $\Theta$ used for representing the variables $x_j$ as a fixed-point approximation in terms of power of 2 (Eq. 10). The general form of the set $\Theta$ is defined as⁴⁸:

$$\begin{aligned} \Theta = \{ 2^l: l \in [o, p] \wedge l, o, p \in \mathbb {Z} \}, \end{aligned}$$

(16)

where $l$ is contiguous integer values within the interval $[o, p]$, and $o$ and $p$ are the lower and upper bounds of the interval. In our scheme, instead of using a large range for $o$ and $p$, we refine the precision of $x_i$ iteratively over multiple loops to achieve higher accuracy. For a 3-bit fixed-point approximation with $x_i \in [0, 2)$, $o = 0$ and $p = -2$. For the $K$ iterations, the total cost for process is:

$$\begin{aligned} \mathcal {O}(mn^2 + mnc + n^2c^2 + K\cdot \mathcal {P}(cn)), \end{aligned}$$

(17)

For Tikhonov regularization methods with the direct solver (e.g., LU decomposition), the computational cost primarily depends on solving the modified linear system with the solution presented in Eq.15. The cost for finding $\textbf{s}$ with a given $\lambda$ is $\mathcal {O}(mn^2 + n^3)$. Selecting the optimal $\lambda$ is crucial for achieving a balance between data fidelity and regularization. One commonly used technique is the L-curve method, which plots the norm of the regularization term $\Vert R(\textbf{s})\Vert$ against the residual norm $\Vert \textbf{D} \textbf{s} - \textbf{T}\Vert$ on a log-log scale^52,53. For $N_\lambda$ values of $\lambda$, the cost of solving the Tikhonov-regularized system for each $\lambda$ is $\mathcal {O}(N_\lambda \cdot (mn^2 + n^3))$. In addition, the computation of the residual norm $\Vert \textbf{D} \textbf{s} - \textbf{T}\Vert$ for each $\lambda$ involves matrix-vector multiplications, contributing a cost of $\mathcal {O}(N_\lambda \cdot mn)$. If we want to find the optimal $\lambda$ automatically, the curvature of the L-curve must be computed. This requires additional operations such as numerical differentiation and curvature estimation, which are negligible compared to the cost of solving the regularized systems. Thus, the total cost for the L-curve method, including the automatic selection of $\lambda$, is:

$$\begin{aligned} \mathcal {O}(N_\lambda \cdot (mn^2 + n^3)). \end{aligned}$$

(18)

From Eqs. 17 and 18, achieving a computational cost of $\mathcal {O}(\mathcal {P}(cn)) < \mathcal {O}(n^3)$ would enable quantum annealing to significantly accelerate problem-solving. Research efforts, such as those utilizing multi-qubit correction techniques⁵⁴, aim to realize this improvement. These approaches can potentially achieve a substantial reduction in computational complexity. This advancement would facilitate rapid convergence to the global minimum and unlock considerable performance gains.