Seismic geophysics relies heavily on subsurface modeling based on the numerical analysis of data collected in the field. The computational processing of a large amount of data generated in a typical seismic experiment can take an equally large amount of time before a consistent subsurface model is produced. Electromagnetic reservoir data, like CSEM (Controlled Source Electromagnetic), petrophysical techniques, such as electrical resistivity and magnetic resonance on multi-wells, and engineering optimization problems like reservoir flux simulators, well field design and oil production maximization also need a strong computational apparatus for analysis.

On the other hand, in the past decade, there has been much progress in the development of quantum computers: Machines exploiting the laws of quantum mechanics to solve hard computational problems faster than conventional computers. A concrete example of such progress is the so-called quantum supremacy, that has been recently demonstrated using specific purpose quantum computers [1–3]. The Geoscience field and related industries, such as the hydrocarbon industry, are strong candidates to benefit from those advances brought by quantum computing.

Currently, different quantum technologies and computational models are being advanced. Giant companies like IBM, Google, and Intel are developing quantum computers based on superconducting technologies [4]. Other companies are also putting considerable effort into building a fully functional quantum computer based on Josephson junctions, such as the North American Rigetti, whereas, the also American IonQ and the Austrian AQT are working on computers based on trapped ions [5]. The Canadian company D-Wave, leader in the computational model known as quantum annealing [6], is already trading quantum machines, and the also Canadian Xanadu is providing cloud access to their photonic quantum computer [7, 8]. Recent reviews comparing superconducting and trapped-ion technologies and different cloud-based platforms can be found in Refs. [9, 10], respectively.

In the field of Geoscience, recent works have used quantum annealers to hydrology inversion problems [11, 12]. In those works, it was shown that, although the size of the problem that can be solved on a third-generation D-Wave quantum annealer is considered small for modern computers, they are larger than the problems solved with similar methodology with Intel’s third and fourth generation chips. It is also important to mention that optimization techniques are widely used in seismic inversions, but usually classical algorithms get stuck in local minima. Previous works [13–15] have indicated that quantum annealing can be advantageous to solve seismic problems. However, the potential applications of quantum computing in Geoscience has so far been largely unexploited in the specialized literature.

In this work, we present a formulation of a seismic problem as a binary optimization, and one small scale subsurface seismic problem is solved using the D-Wave quantum annealer available in the Amazon Braket service. We have evaluated the performance of the quantum computer comparing the results obtained to those derived from a classical computer. Our analysis was focused on the accuracy. It was found that the accuracy achieved by the quantum computer is at least as good as that of the classical computer for the problem we have studied.

This paper is organized as follows. In Section 2 we introduce the basic idea of quantum annealing. In Sections 3, 4, we present the formulation of a seismic inversion as a binary optimization problem and the methods used, respectively. In Section 5, the results obtained in the D-Wave quantum annealer are shown. In the last section, we draw the conclusion.

In the literature, there are many different quantum computational models developed. Currently, three main models of quantum computing are being considered: The logic gate model, or circuit model, boson sampling, and the adiabatic model. The gate model is an universal quantum computation model that is performed programming a step-by-step instruction build from basic building blocks, known as quantum gates, similar to the classical circuit model [16]. This model is exploited by companies such as IBM, Intel, Rigetti, IonQ, AQT, and Google. Boson sampling computers consists in sampling from the output distribution of bosons in a linear interferometer [8, 17]. There is strong evidence that such an experiment is hard to be simulated in classical computers, but it is efficiently solved by special purpose photonic chips. The adiabatic computation model [18, 19] is a model in which the computational problem is mapped into a Hamiltonian, in such a way that the solution of the computational problem is encoded in the ground state of the quantum system represented by the Hamiltonian H _final . The computation is performed starting from the ground state of a known Hamiltonian H _initial . The Hamiltonian is slowly modified towards to the target Hamiltonian H _final . During the process, the total Hamiltonian of the system is given by H p = 1 − p H i n i t i a l + p H fi n a l , (1) where p ∈ [0, 1]. According to the adiabatic theorem, if the evolution is performed adiabatically, the quantum state of the system remains in the instantaneous ground state throughout the entire process.

The adiabatic model is equivalent to the gate model [20], i.e., any problem that can be computed in the gate model will also be computable in the adiabatic model. This statement is valid for certain types of k − local Hamiltonians [20]. For example, considering a system of n qubits, we can perform universal adiabatic computation by choosing [21]. H i n i t i a l = ∑ i δ i σ i x , (2) H fi n a l = ∑ i H i σ i z + ∑ j > i J i , j σ i z σ j z + ∑ j > i K i , j σ i x σ j x , (3) where σ i k , with k = x, y, z, is one of the Pauli matrices of the i ^th qubit, δ _i and H _i are local transverse and longitudinal fields, respectively, and J _i,j and K _i,j are coupling constants.

The adiabatic model can be viewed as a special case of the quantum annealing computing. In quantum annealers, the Hamiltonian change is not adiabatic. Therefore, quantum annealing is a heuristic type of computation. The D-Wave computer is a quantum annealer that uses Ising chains H fi n a l = ∑ i H i σ i z + ∑ j > i J i , j σ i z σ j z . (4)

The quantum annealing performed with Ising chains is unlikely to implement universal quantum computation [22]. Therefore, D-Wave annealers are more restrictive than the universal adiabatic model. Quantum annealing with Ising chains can be applied to a class of computational problems known as NP-hard problems [23]. It is believed that quantum annealing will be able to find better approximate solutions or find such approximate solutions faster then classical computers [22]. Here, it is also important to mention that the advantage of quantum annealers over classical methods is still under debate [24–26]. Recent works have proposed that, in some cases, quantum annealing is advantageous over classical computing [27–30], on the other hand, no advantage was reported in Refs. [24, 31]. The origin of the possible speedup is also under debate. Quantum tunneling is often claimed to be the key mechanism underlying possible speedups of quantum annealing. However, recent work has found numerical evidence that quantum tunneling processes can be efficiently simulated by Monte Carlo methods [32]. There is also evidence to suggest that it is unlikely to achieve exponential speedups over classical computing solely by the use of quantum tunneling [33]. The role of the temperature in the performance of quantum annealing has been also studied in [34].

To solve a problem in the D-Wave quantum annealer [6], it is necessary first to express the problem to be solved as an Ising problem or as a quadratic unconstrained binary optimization (QUBO), which is equivalent to Ising but defined on the binary values 0 and 1, whereas the Ising problem is defined on the binary values ±1/2. The QUBO problem can be written as the minimization of the quadratic function f q = q T Q q , (5) where q ∈ {0,1} ⁿ , Q is a n × n upper (lower) triangular matrix and the vector q T = ( q 1 , q 2 , … , q n ) contains n binary variables. QUBO problems are commonly used in machine learning and many important computational problems can be translated to a QUBO formulation as well [35]. Examples of problems that have been addressed with a D-Wave quantum annealer are: the classification of human cancer types [36], traffic optimization [37], transcription factor DNA binding [29], metamaterial designing [38] and Higgs boson data analysis [39].

Recently, there has been a growing interest in quantum algorithms for systems of linear equations, Ax = b, where A is a n × n matrix and b is a unit vector. Such algorithms may find applications in different research areas, including Geoscience. In the quantum gate model, the quantum version of such problem is called Quantum Linear Systems Problem (QLSP) [40, 41], and it is defined as the problem of preparing the state | ψ 〉 = ∑ i n x i | i 〉 ∑ i n | x i | 2 , (6) where x = ( x 1 , x 2 , … ) T is the solution of Ax = b.

In 2008, Harrow, Hassidim, and Lloyd (HHL) proposed a quantum algorithm for the QLSP problem [42]. Given some assumptions [43], the run time of HHL is O(k ² s ²  log(n)/ϵ) where k is the condition number of the matrix A, defined as the absolute value of the ratio between the largest and smallest eigenvalues of A, s is the sparsity of A, defined as the number of nonzero entries per row, and ϵ is the desired precision.

After the initial HHL proposal, several improvements were achieved: the condition number dependence was reduced from k ² to k log ³(k) [44], the error dependency was reduced from 1/ϵ to a polynomial function in log(1/ϵ) [45], and a sparsity-independent runtime scaling was achieved in [46]. The QLSP problem can also be solved using iterative quantum solvers in runtime O(k ²  log(k)/log(1/ϵ)) [47] and with runtimes O(k log(k)/ϵ) and O (k ²  log(k)/ϵ) using the evolution randomization method, a simple variant of adiabatic quantum computing where the parameter p in (1) varies discretely, rather than continuously [48]. The best general-purpose classical conjugate gradient algorithm to solve Ax = b has the runtime O(nks log(1/ϵ)). Here, we must emphasize a fundamental difference between classical and quantum algorithms. While the conjugate gradient returns the solution vector x, quantum algorithms return a quantum state, Eq. 6, that approximately contains all the components x _i of the solution vector x. It is possible to obtain any specific entry x _i by measuring the output state (6), but in general, it will require repeating the algorithm many times, which would kill the exponential speedup. Still, quantum algorithms can be used as subroutines in different applications [40, 43].

In quantum annealers, the problems of solving a system of linear equations and a system of polynomial equations were previously studied in [49–51] and [52], respectively. Unlike the previously mentioned quantum algorithms, quantum annealing solves Ax = b completely, i.e., it returns the vector x. To compare the performances of a quantum annealer and a classical computer, we must take into account the cost to prepare the problem, i.e., the procedure to map Ax = b into a QUBO or Ising problem, the cost to perform the annealing, and also the cost to post-process the results. The performance of quantum annealers to solve linear systems was studied in detail in Ref. [50]. It was shown that quantum annealers might be competitive if there exists a post-processing method that is polynomial in the size of the Matrix A with a degree less then 3.

We considered the propagation of sound waves in a multi-layered medium, as shown in Figure 1. Multiple sources produce sound waves that can be reflected in the interface of each layer. Assuming that the wave propagation can be modeled as narrow beams or rays, the sound trace originated in the i ^th source reaches the i ^th detector after the time interval t i = 2 ∑ j = 1 i d i j v j , (7) where d _ij and v _j are the distance traveled by the sound waves and the sound speed in the j ^th layer, respectively. If we consider the thickness of each layer as h _j and the distance between two consecutive sources (detectors) as Δ _i , we can write d _ij = h _j /cos θ _j where θ _j = arctan(Δ _j /h _j ).

The layered model described above is commonly used in seismic explorations, either offshore or onshore [53]. In seismic experiments the goal is to determine the velocities {v _j } from the time intervals {t _i }, by solving the system M s = t , (8) where t T = ( t 1 , t 2 , … , t m ) , s T = ( 1 / v 1 , 1 / v 2 … , 1 / v m ) is the slowness vector, M is a m × m lower triangular matrix with nonzero elements given by M _.,j = 2h _j /cos θ _j , and m is the number of layers.

In order to use a quantum annealer to solve the above seismic problem it is necessary to translate the problem into a QUBO formulation. To proceed, first, we rewrite the system (8) as a minimization problem with the objective function f s = ‖ M s − t ‖ 2 , (9) where s ∈ R m . Next we write the slowness vector as s = s ₀ + L(x − I), where L defines the bounding limits of s, 0 ≤ x _i < 2 ∀ x _i , I T = ( 1,1 , … ) and the vector s ₀ is an initial guess for s. The objective function f(s) is rewritten as f x = ‖ M x − b ‖ 2 , (10) where x ∈ R m and b = (t + L MI − Ms ₀ )/L. The matrix M and the vector b are parameters of the objective function while the vector x must be converted into a binary format. Here we discretized each x _i variable with the R-bit approximation x i = ∑ r = 0 R − 1 q i , r 2 − r . (11)

To formulate our QUBO problem, we construct a new binary vector q and a new real matrix A in order to form the binary system of equations Aq = b. It is straightforward to reformulate this system as a binary optimization problem [49–51], whose solution vector, q*, is given by q ∗ = arg min q ∈ 0,1 R × m ‖ A q − b ‖ 2 (12) = arg min q ∈ 0,1 R × m q T Q q + C , (13) where Q i i = ∑ k = 1 m A k i A k i − 2 A k i b k , (14) Q i j = ∑ k = 1 m 2 A k i A k j ,  for  i < j , (15) and C = ∑ k = 1 m b k 2 is an additive constant that does not change the ground state.

The precision of the solution depends on how many binary digits are used to represent the real variables of the problem, a solution with good precision would consume a large number of qubits of the quantum hardware. Here we have used a recursive approach similar to what was used in [49] to improve the precision of floating-point division. Our recursive approach is described in the Algorithm 1, using such an approach we could improve the solution of a system of linear equations with 46 real variables, using just a few qubits to represent each variable. Next, we will show in our example that using R = 3, in Eq. 11, and carrying 20 iterations is sufficient to reach a good solution.

Algorithm 1

Function to solve the system of equations Ms = t. The vector s ₀ is the initial guess for s, ϵ is the tolerance, N _max is the maximum number of iterations, and L defines the interval where we expect to find the solutions (see text).

We have performed the quantum computation with the D-Wave computer provided by the Amazon Braket service. Currently, there are two versions of quantum annealer available in Amazon Braket. The first is the D-wave 2000Q version, this computer contains 2041 working qubits. The connections among the qubits are represented by a graph called Chimera [54]. In this topology, each qubit is coupled to no more than 6 other qubits. We can call Chimera as a graph of degree 6. The second version is the D-Wave Advantage system, it is a more advanced computer with 5,436 working qubits disposed in a Pegasus graph with degree 15 [55]. A QUBO problem can also is represented by a graph, where each vertex of the graph corresponds to a binary variable q _i . When the QUBO problem is represented by a graph with degree greater then 6, for Chimera, or 15, for Pegasus, it is necessary to embed the QUBO graph onto the chip topology. The present seismic problem, for example, is represented by a full graph, where each vertex is connected to all other vertices. In this work, the embeddings were obtained using a heuristic algorithm [56] provided by D-Wave.

Errors during the computation are important issues to be considered. In quantum annealing, the computational problem is encoded into the ground state of the Ising Hamiltonian (4), the gap between the ground state and the excited states is a key property of the system. If the gap is too small, thermal excitations and non-adiabatic transitions can induce transitions, as a consequence, the computer will output an excited state, which can be viewed as a computational error [57]. In addition, the wrong implementation of the Hamiltonian (4) may result in the wrong ground state [57]. We have noticed that the ground state was achieved with high probability. Therefore, unwanted transitions due to thermal fluctuations or non-adiabatic evolution do not represent an important issue in our case. The annealing time used was 20 μs, the minimal and default value of the D-Wave machine. To post-process the output we use default D-Wave’ routines.

In our implementation, analog errors are the most important, i.e., inaccurate implementations of the parameters H _i and J _i,j described in the Eq. 4. To reduce the impact of analog errors, we have used 10 spin-reversal transforms, also known as gauge transformations [58, 59]. This type of transformation is based on the fact that the structure of the Ising problem is not affected when the following transformations are applied: H _i → g _i H _i and J _i,j → g _i g _j J _i,j , where g _i = ±1. The original and transformed problems have identical energies. However, the sample statistics are affected by the spin-reversal transform because the quantum hardware is a physical object subject to errors.

We have applied the above formulation to solve a small scale underwater seismic inversion problem. Artificial data were generated by simulating sound traces in ocean. We have used the sound speed profile of the Philippine Sea, available to public [60], as shown in Figure 2. From the simulation, we obtained the travel times between the sources and detectors. The seismic inversion model was constructed using 46 layers while the real variables were digitalized with R = 3 bits. This model yields 138 binary variables, for the present type of problem, it is the maximum number of variables that we can embed in the working graph of the D-Wave Advantage system available in Amazon.

To perform the seismic inversion we have considered that all layers in the model match the position where the sound wave is reflected, as shown in Figure 2. The solution is presented in Figure 3, as can be seen, good results are obtained with 20 iteractions. We have also performed the inversion with a classical computer running the forward substitution algorithm to invert the lower triangular matrix in Eq. 8. When classical and quantum inversions are compared, we found that the relative error between them was ≈ 1 0 − 4 , as shown in Figure 4. This result shows that the quantum computer using the recursive approach to solve a system of linear equations has enough control to find solutions with good precision.

The results could also be compared to the benchmark provided by the condition number of the numerical problem at hand. When performing a numerical inversion procedure on the lower triangular matrix M to recover the solution s in Eq. 8, a straightforward estimation of the lower bound of the relative error ϵ _s in s arising from the relative error ϵ _t in the righthand side vector t due to the numerical conditioning of M is given by [61]. ϵ s = κ ∞ ϵ t , (16) where κ _∞ is the condition number given by κ ∞ ≥ max j h j / cos θ j min j h j / cos θ j . (17)

Therefore, the condition number relates the error associated with t, ϵ _s , and the error of the solution s, ϵ _s . In this particular application of seismic inversion, where distances are generally in the scale of 10³ m, it is reasonable to assume that the figures are accurate up to the order of ≈ 10 m , with the next order of magnitude (1m) giving the scale of the error ϵ _t . Given that the error ϵ _t can be estimated as ≈ 1 0 − 3 and the condition number is κ _∞ ≈ 1, then we can estimate ϵ _s ≈ 10⁻³ from (16). This shows that the error of the quantum approach comes from the inversion problem by itself and is of the same order as the classical approach.

Quantum computing represents a fundamentally different paradigm, an entirely different way to perform calculations. The Geoscience field and related industries are strong candidates to benefit from it. However, the performance of Geoscience inversion problems on current available quantum computers has so far been largely unexploited.

In this paper, we have solved a seismic inversions in a D-Wave quantum annealer. The seismic problem was written as a system of linear equations and then translated into a QUBO formulation. The results presented here indicate that the current available quantum annealers can solve a seismic inversion at a relatively small size with nearly the same accuracy as a classical computer. The proof-of-principle computations performed here show some promise for the use of quantum annealing in Geosciences.

The practical use of quantum annealing in Geosciences will require the ability to solve large problems. To address this issue, decomposer tools have been proposed to divide a large problem into small subproblems which can be solved individually by the quantum hardware [62–64]. Using such approach, it is possible to solve a large-sized problem using just a limited number of available qubits. Another interesting approach is the reverse annealing. Within this method, one starts from known local solutions which can be obtained in a classical computer. The annealing is performed backward from the known classical state to a state of quantum superposition, then proceeding forward it is possible to reach a new classical state that is a better solution than the initial one. Recently, it has been shown that it is possible to refine local solutions with recursive applications of reverse annealing [65–68]. We believe that the development of hybrid quantum-classical methods, such as mentioned above, will be essential to solve complex seismic problems on quantum annealers in the near future.

Finally, we should mention that the problem solved here is well-conditioned, however, often Geoscience problems are ill-conditioned. An interesting question is whether quantum computers can solve ill-conditioned problems efficiently. In Ref. [69], it was theoretically proposed that preconditioning methods can expand the number of linear systems problems that can achieve exponential speedup over classical computers. In future works, we plan to study the performance of the quantum annealers to solve ill-conditioned problems. Another attractive prospect for future work is the implementation of Geoscience problems in gate model computers, based either on superconducting or trapped ions technologies.