Nonlinear amplitude inversion using a hybrid quantum genetic algorithm and the exact zoeppritz equation

Ji-Wei Cheng ,Feng Zhng ,*,Xing-Yng Li ,b

a College of Geophysics,State Key Laboratory of Petroleum Resources and Prospecting,China University of Petroleum(Beijing),Beijing,102249,China

b British Geological Survey,Edinburgh,UK

Keywords:

ABSTRACT

1.Introduction

Prestack amplitude variation with offset(AVO)inversion is one of the commonly used techniques to estimate elastic parameters of the subsurface(Ostrander,1982).The recovery of the P-wave velocity,S-wave velocity,and density from the PP-wave can be difficult due to a number of reasons that include the nonlinear nature of the forward problem,large trade-offs between density and velocities,and band-limited and noisy characteristics of seismic data(Nicolao et al.,1993;Downton,2005;Zhang et al.,2018;Wang et al.,2019;Zhang and Li,2020).Several approximations of the reflection coefficient are proposed to simplify the relationship between the seismic response and elastic parameters assume weak contrast and limited offset/angle(e.g.,Aki and Richard,1980;Shuey,1985;Smith and Gidlow,1987;Fatti et al.,1994).Thus,to account for different types of contrasts and wide range of offsets,the nonlinear Zoeppritz equation has been used in AVO inversion(e.g.,Lu et al.,2015;Zhi et al.,2016;Pan et al.,2017).

Varied inversion methods have been proposed to solve the nonlinear inverse problem(Stoffa and Sen,1991;Wei et al.,2011;Sen and Stoffa,2013;Yin et al.,2013;Li and Mallick,2015;Mallick and Adhikari,2015;Zhang et al.,2015;Aleardi and Mazzotti,2017;Huang et al.,2017;Zu et al.,2017b;Aleardi et al.,2019;Liu et al.,2020a,2020b,2021;Ma et al.,2020;Yao et al.,2020;Yu et al.,2020;Zhou et al.,2020).Linear inversion schemes can solve nonlinear inverse problems by carrying out linearization about an initial model,and then iteratively modifying the initial model until the convergence conditions are achieved,and these linear inversion methods are computationally efficient and require a relatively smaller memory(Buland and Omre,2003;Yin et al.,2016;Zu et al.,2017a;Zhang et al.,2019).However,they are heavily dependent on the initial/starting model which may result in convergence to a local solution(Liu et al.,2020c).Nonlinear inversion methods have a higher likelihood of obtaining the global solution and hence,more suitable for nonlinear seismic inverse problems based on the exact Zoeppritz equation(Lu et al.,2015;Lehochi et al.,2015;Zhi et al.,2016;Pan et al.,2017).Although the ideal inversion method for such problems is the exhaustive search method,it is unrealistic as the optimization technique requires infinite computation time irrespective of the computing power of available resources.The Monte Carlo method is a global optimization technique that has found applications in different fields of study.An inherent weakness of this technique is its requirement of trials in the order of thousands or even millions to obtain acceptable precision(Gentle,2003).With the availability of modern powerful computing resources,the implementation of parallel computation has enabled the deployment of some classic meta-heuristic algorithms based on Monte Carlo ideas,such as simulated annealing(which originated from statistical mechanics),and genetic algorithm(an abstraction of biological evolution)(Stoffa and Sen,1991;Dorigo,1992;Mallick,1995),in the inversion of seismic data.The genetic algorithm(GA)is an inherently highly parallel search heuristic inspired by the theory of natural selection.It simulates a range of solutions and has no strict restrictions on the initial model.Since its efficiency is mainly influenced by the population size,convergence may be premature,reducing the likelihood of obtaining the global solution if a small population size is used(Laboudi and Chikhi,2012).

Quantum computing methods have been hailed as the future of computing sciences(Moradi et al.,2018).This can be attributed to their theoretically proven supercomputing speed,better stability and effectiveness.Although many studies have been conducted to examine the applicability of quantum computing in general(e.g.,Moradi et al.,2018;Liu et al.,2018),its potential in geophysics has received limited attention.Han and Kim(2000)proposed the quantum genetic algorithm(QGA)which combines quantum computing with the concept of biological evolution.Compared with the classical GA,the QGA uses a small population size to generate a larger and more diverse population,and a quantum rotating gate to update the solutions.This improvement in QGA results in improved convergence efficiency and accuracy of results in optimization problems(Lahoz-Beltra,2016).In spite of these improvements,the QGA shows two major deficiencies when solving relatively more complex problems(Han et al.,2001):Firstly,slow rate of convergence.The angle of rotation which is usually fixed and the direction of the quantum rotating gate are obtained from a look-up table in Han et al.(2001),making the algorithm less flexible and resulting in a slower convergence rate in the early stages.Secondly,high tendency of been trapped in local minima.The QGA uses the quantum rotating gate to update genes with best fit.This increases the population of the best fit genes but ultimately results in loss of diversity which increases the tendency of been trapped in local minima.

To address these shortcomings in the QGA method,in this study,we propose an improved version called the hybrid quantum genetic algorithm(HQGA).It combines a self-adaptive search strategy and the operations of quantum mutation and catastrophe,enjoying the advantages of quantum computing and the GA.This positions the HQGA as a great tool for global optimization using a small population size.It also enjoys the advantages of the results been independent of the initial model and high likelihood of obtaining the global solution.We verified its reliability and stability by conducting synthetic tests with a model based on actual logging data.Similar tests were also conducted with real field data.The results show that the HQGA is a strong candidate for global optimization with fast convergence speed and robust stability.

2.Methodology

2.1.Forward modeling

Based on the three-dimensional wave equation with boundary conditions of continuous displacement and stress at the interface of two media,Zoeppritz(1919)deduced expressions for the reflection and transmission coefficients for a seismic wave that impinges on the interface.These expressions are known as the Zoeppritz equations.We assume a solid-solid interface between two homogeneous isotropic elastic half-spaces,where P-wave velocity,S-wave velocity and density of the upper and lower half-spaces are denoted by VP1,VS1,ρ1,and VP2,VS2,ρ2,respectively.The incident angle and angle of transmission of the P-wave and S-wave are i1and i2,and j1and j2,respectively.The ray parameter is constant and given as p=sin i1/VP1=sin i2/VP2=sin j1/VS1=sin j2/VS2.Aki and Richards(1980)derived an accurate solution of the Zoeppritz equation for Rppas

where

and

As shown in Equations(1)-(3),although the Zoeppritz equations are highly nonlinear functions with respect to these properties,the formulations are explicitly valid in isotropic media,and they use few approximations,which allows the prediction of reflection coefficients to be accurate from near to far incident angles.We deploy global optimization schemes to perform the AVO inversion using the exact Zoeppritz equation.This combination can improve the chances of obtaining the global solution.

2.2.Quantum genetic algorithm

Nonlinear inversions can be performed by minimizing an L2 norm objective function:

where r is the n×N-by-1 observed data vector consisting of seismic amplitudes for different incidence angles at each time sample,and argmin returns the value of m which minimizes the function.The number of angle traces and data samples for each trace are n and N+1,respectively,while R is the nonlinear forward operator.The model parameter m=[VP,VS,ρ]Tis the 3(N+1)-by-1 vector consisting of the three model parameters of interest.

Using either GA or QGA,the population of any model parameter for each time sample for the tthgeneration is expressed as

where mtqrepresents the qthindividual chromosome of a model parameter,and v is population size.A larger v implies higher probability for obtaining the global solution but with a longer computation time.

GA represents each chromosome as a bit string where s represents the length of each chromosome and is also the number of genes;whileτtqxtakes either 0 or 1.Coding the physical properties using genes limits the search space and defines the resolution of each of the model parameters.A larger s means higher resolution but longer computation time.When bit string representations of integers are used,Grey coding is often employed.

On the other hand,a chromosome at the tthgeneration in QGA can be described as

A gene called qubit is the basic unit of information in quantum computation,which can be represented by a unit vector of a twodimensional Hilbert space(Williams,2010)and described by a superposition of the basis states:

where|0〉and|1〉denote two basic states(Han and Kim,2000)and αandβare complex numbers that satisfy:

|α|2and|β|2are called the probability amplitude of corresponding states|0〉and|1〉of qubit,respectively.Since the application of the quantum genetic algorithm is mostly performed by a classical computer,the gene qubits are expressed in the sine and cosine form:,as it appears more concise and intuitive(Silveira et al.,2017).

The quantum rotating gate is the core operator in the evolution operation of the quantum genetic algorithm.Hence,it directly affects the performance of the algorithm.Since model parameters are in superposition states,the genetic qubits in the population should be updated by quantum rotating gates to adjust the probability amplitude and constitute new individuals.The rotation strategy adopted is given by the following equation:

whereθUis the rotation angle,|φ′〉is the updated quantum superposition state,and cosθ′and sinθ′are the probability amplitudes of the quantum state after rotation.Generally,the rotation angle and rotation direction of the quantum gates are empirically determined in advance(Han and Kim,2000;Layeb and Saidouni,2007),while their adjustment strategy is updated in accordance with the adjustment table.

The main challenge of optimization algorithms is its computational complexity(Laboudi and Chikhi,2012).The genetic algorithm has a computational complexity of O(wv),where w is the number of function convergence and v is the population size,while the proposed algorithm finds a solution in O(w log v)time(Nielsen et al.,2010).We also note that unlike a classical bit in GA,a qubit is not a value of 0 or 1,but a probability to be represented with values,which reduce the correlation of individuals and the number of iterations required for the globally optimal value.These show QGA can save much computational expense especially for large population sizes,where traditional GA requires a larger population size and more computational expense to reach a similar accuracy as the new method.

We use the Rastrigin function to evaluate the GA and QGA.As shown in Fig.1,many local minima make it difficult to obtain the global minimum of the Rastrigin function.The global minimum of the function is located as indicated by the vertical line,while local minima regularly distribute around the global solution.At any local minima other than[0 0],the Rastrigin function is greater than 0.The allowable parameter range is[-5,5],the accuracy is 0.1,and simulation tests are performed 50 times.From Table 1,we observe that QGA outperforms GA for all model space dimensions and the number of convergence increase with increasing dimensions of the model space as expected.The QGA is a robust stochastic global search method capable of handling multiple minima,but its computational expense increases and inversion accuracy decreases with increasing nonlinearity of the inversion problem.There are,therefore,some potentials for QGA.

Table1 The average number of model convergence required to converge to the global minimum using the GA and QGA for different model space dimensions.

2.3.A hybrid quantum genetic algorithm

Here,we present a new hybrid quantum genetic algorithm(HQGA)which is an improved variant of the traditional QGA to solve the nonlinear AVO inversion problem.A self-adaptive rotation strategy is proposed to adjust the rotation angle according to the convergence.Quantum mutation and catastrophe are incorporated to enhance the local search ability,and to obtain the globally optimal model parameters,respectively.The Rastrigin function with dimension 2 is used to test the effect of the different proposed improvement measures on accuracy and efficiency.See Table 2 for list of the average convergence number to reach the termination condition and the optimal solution of the results associated with the different improvement measures.

Table2 The average number of model convergence required to converge to the global minimum and the optimal result after improvements.

2.3.1.Self-adaptive rotating strategy

A fixed rotation angle of the quantum rotating gate makes the QGA less flexible.It does not maximize the full potentials of quantum computing(Han and Kim,2000).For the HQGA,we propose a self-adaptive rotating strategy which takes care of this shortcoming inherent in the QGA.

The angles of qubits[αB,βB]Tand[αi,βi]Tare denoted asθbandθiin a unit circle.Then,we use the new representation:

Fig.1.2D Rastrigin function.The vertical line indicates the position of the global minimum.

When A≠0,if 0＜|θi-θb|＜π,sgn(Δθ)=-sgn(θi-θb)=-sgn(sin(θi-θb))=-sgn(A),while ifπ＜|θi-θb|＜2π,sgn(Δθ)=-sgn(θi-θb)=-sgn(sin(θi-θb))=-sgn(A).When A=0,|θi-θb|=0 orπ.At this time,the positive and negative rotating effects are the same,so that the direction of the rotation angle,θ,can be either positive or negative.

We applied a dynamic rotation method to adjust the rotation angle of the quantum rotating gate based on the evolutionary process(Li et al.,2005).When the difference in fitness between the current best individual and the individual in the current generation is large,the current individual is updated at a large angle to facilitate the speed of convergence.On the other hand,when this difference is small,the individual in the current generation is finetuned to obtain the best individual.The rotation angle is defined as:

whereθ?[0?001π,0?05π].θmaxandθminrepresent the maximum and minimum values in the interval,respectively.Fminrepresents the fitness of the best individual while FCurrepresents the fitness of the current individual,and C is a positive integer constant.C is a trade-off constant to adjust the rotation angle.When C is large,the rotation angle will start with a relatively small value and converge with slow speed,while for small C,the rotation angle become too large and convergence to the global minimum becomes more unlikely.Accordingly,extreme care should be taken in assigning a value to C which should be informed by the requirements of the experiments.Our improved method has simplified the determination of the rotation direction by dynamically adjusting the rotation angle,thereby improving the convergence rate without losing accuracy.Results listed in Table 2 show that the average convergence number required for the global minimum reduces from 301 to 228 with the application of the self-adaptive rotation strategy,and the optimal solution is closer to the global minimum.

2.3.2.Quantum mutation

Quantum mutation is an assistant operation whose aim is to avoid the loss of important information on the population and enhance the local search ability of the quantum genetic algorithm(Han and Kim,2000).It is required because the diversity maintained by the superposition state in practical applications is not sufficient to reflect the superiority of the genetic algorithm which could lead to a premature convergence aided by the effect of noise.

The quantum NOT gate(Pauli-X gate)is commonly adopted to perform quantum mutation(Suo,2020).The probability amplitude of each gene|0〉or|1〉after measurement is exchanged,and it rotates the angle of the gene by-2θeach time.The rotation angle is slightly increased and the convergence speed of the algorithm is accelerated in the early stages.This may cause a premature phenomenon in later stages.The matrix representation of the Hadamard gate is:(Djordjevic,2012),and the gene mutation is given as:

From Equation(13),we find that the angle of the qubit is rotated by-2θ.It is clear that the Hadamard gate is more suitable as a mutation operator when the convergence is in the later stages which prevents the information for the current optimal individual from disappearing.The rotation processes of Pauli-X and Hadamard gates are shown in Fig.2.

Based on the value of the threshold K,the quantum mutation strategy can be described as follows:1.For K＜FCur-Fmin,the algorithm is said to be in the primary stage,where only a small number of individuals are having optimal fitness values,and others requiring significant update.In this case or stage,the quantum NOT gate is used as the mutation operator,and the large rotation angle of quantum mutation increases the rate of convergence of the algorithm.2.For K>FCur-Fmin,the algorithm is considered to be in the middle and late stages.The high-fitness individuals in the population are in the majority and only a small number of individuals are required to adjust the rotation angle of quantum mutation,ensuring that the information for the current optimal individual in the population remain unaffected by mutation and the population is stable.It is,therefore,appropriate to use the Hadamard gate as the mutation operator at this stage.Table 2 shows that using quantum mutation leads to higher efficiency,decreases the average calculation number to 157,higher accuracy,and the optimal solution closes to zero.

2.3.3.Quantum Catastrophe

Quantum catastrophe is an essential operation that ensures that the search is not trapped in local minima,especially for algorithms having strong likelihood of convergence to local minima(Sen,2010).When solving for functions like the Zoeppritz equation with multiple peaks and high nonlinearity,the quantum genetic algorithm may be trapped in a local optimal solution.In such cases,the quantum catastrophe strategy ensures the search escapes the current local solution space and continues to the global space.The quantum catastrophe pseudocode program is called Algorithm 1.

The Rastrigin function is used to illustrate the advances of the HQGA.We set model space dimension 2,and the global solution of the function is[0,0].The average number of model convergence required to obtain the minimum value and the optimal results after improvements to the QGA are recorded in Table 2.The result shows that the optimal solution using the HQGA is approximately equal to the global solution and has highest accuracy among the three different improved variants.Although the average convergence generation of HQGA is larger than that of the case without catastrophe operation,it is much smaller compared to other improved cases.The catastrophe operation always happens late in the convergence and inevitably increases the speed of convergence,making the HQGA a promising global optimization algorithm with high efficiency.Therefore,the HQGA is a promising global optimization algorithm with high efficiency.

Fig.2.Rotation process of Pauli-X gate and Hadamard gate.

Algorithm 1.Quantum Catastrophe.

2.3.4.Procedure for nonlinear inversion based on the hybrid quantum genetic algorithm

The procedure for a nonlinear inversion using the HQGA is displayed in brief in Table 3 and described as follows:

Table3 Procedure for nonlinear inversion using HQGA.

(i)Initialization.This entails the initialization of the inversion parameters according to the number and complexity of inverted attributes.Using the application to real data implemented in this study as an example,we set the maximum generation to 3000,population size(v)to 40,chromosomes of length(s)to 10,the range of rotation angle(Δθ)to[0?001π,0?05π],initial probability amplitude of qubits(αandβ)to 1/√---2,the stopping criterion to 10-5,the probability of mutation(pm)to 0.2,the threshold of mutation gate(K)to 0.01,the condition of catastrophe to 10% of the maximum generation.Once catastrophe occurs,3/4 population and 2/3 chromosome will be initialized.With these parameters an initial population m(t)=｛mt1,…,mtq,…mtv｝is randomly generated.

(ii)Measurement.Make m2bit(t)by observing m(t)states and a number between 0 and 1 is randomly generated.If the number is greater than the probability amplitude,the bit is taken as 1,otherwise,it is 0.

(iii)Evaluation.The measured binary string is then decoded to a decimal value according to the constraints of the parameter.The group model parameters are used for forward modeling.We calculate the fitness value from the objective or fitness function,which is defined as follows:

Fig.3.Real well logs are used to generate the synthetic data.P-wave velocity,S-wave velocity and density are shown from left to right.

(iv)Select and update memory.The best solution b in the current generation is located and recorded before updating the memory by saving other individual solutions.

(v)Update quantum gate.The current population associated withtheconventionalquantumrotatinggateis updated using the proposed rotation strategyin consonance with the difference in fitness between the current solution and the best solution.

(vi)Quantum mutation.The quantum mutation is performed using the Quantum NOT or Hadamard gate depending on the performance of the algorithm.

(vii)Quantum catastrophe.The best individuals in the current population are saved.One-tenth of this population is used for re-initialization in order to preserve the diversity of the population.

Fig.4.Synthetic angle gathers with different levels of noise:(a)without noise,(b)signal-to-noise ratio=4,and(c)signal-to-noise ratio=2.

Fig.5.Comparison between the real logs(black solid)and inversion results for a noise-free angle gather using the(a)HQGA(red),(b)QGA(blue),and(c)GA(green)with 2000 generations.The upper and lower bounds of search windows are denoted by grey curves.(d)The differences between inversion results using the three methods and the true values.(e)The errors distribution of predicted P-wave velocity(left),S-wave velocity(middle),and density(right)for a noise-free angle gather with 2000 generations using the GA(top),QGA(middle)and HQGA(bottom),respectively.

3.Experiments with synthetic and field data sets

3.1.Synthetic test

The real well logs in Fig.3 are used to construct synthetic PPwave angle gathers displayed in Fig.4 based on the exact Zoeppritz equation.Two different levels of random noise(SNR=2 and 4)are added to the synthetic gathers as shown in Fig.4(b)and(c).

Fig.6.Comparison between real logs(black solid)and inversion results for noise-free angle gather using the(a)HQGA(red),(b)QGA(blue),and(c)GA(green)with 8000 generations.The upper and lower bounds of search windows are denoted by grey curves.(d)The differences between inversion results using the three methods and the true values.(e)The errors distribution of predicted P-wave velocity(left),S-wave velocity(middle),and density(right)for a noise-free angle gather with 8000 generations using the GA(top),QGA(middle)and HQGA(bottom),respectively.

In order to evaluate the convergence performance of different algorithms,we set the number of generations as 2000.The corresponding inversion results for the GA,QGA,and HQGA are shown in Fig.5(a)-(c).The errors captured as the difference between the predicted and real properties are plotted in Fig.5(d).In general,theHQGA shows smaller errors compared with the QGA and GA.A larger number of these errors are near zero as seen in the calculated distribution of errors shown in Fig.5(e).Dynamic search windows are constructed by smoothing the calibrated well logs.They are applied in all three methods to improve the computational efficiency and resolution.The relative errors(REs=)and correlation coefficients()of results are calculated and shown in Table 4.

Table4 Relative Errors(REs)and Correlation Coefficients(CCs)between the inversion results and the true models using GA,QGA and HQGA with 2000 generations.

Fig.7.Comparison between real logs(black solid)and inversion results for a noise-free angle gather using the(a)HQGA(red),(b)QGA(blue),and(c)GA(green)with 3000 generations.P-wave velocity,S-wave velocity and density are shown from left to right in each subfigure.The upper and lower bounds of search windows are denoted by grey curves.(d)The differences between inversion results using the three methods and the true values.(e)The errors distribution of predicted P-wave velocity(left),S-wave velocity(middle),and density(right)for a noise-free angle gather with 3000 generations using the GA(top),QGA(middle)and HQGA(bottom),respectively.

Fig.8.Comparison between real logs(black solid)and inversion results of a SNR=4 angle gather using the(a)HQGA(red),(b)QGA(blue),and(c)GA(green)with 3000 generations.P-wave velocity,S-wave velocity and density are shown from left to right in each subfigure.The upper and lower bounds of search windows are denoted by grey curves.(d)The differences between inversion results using the three methods and the true values.(e)The errors distribution of predicted P-wave velocity(left),S-wave velocity(middle),and density(right)for a SNR=4 angle gather with 3000 generations using the GA(top),QGA(middle)and HQGA(bottom),respectively.

Fig.9.Comparison between real logs(black solid)and inversion results of a SNR=2 angle gather using the(a)HQGA(red),(b)QGA(blue),and(c)GA(green)with 3000 generations.P-wave velocity,S-wave velocity and density are shown from left to right in each subfigure.The upper and lower bounds of search windows are denoted by grey curves.(d)The differences between inversion results using the three methods and the true values.(e)The errors distribution of predicted P-wave velocity(left),S-wave velocity(middle),and density(right)for a SNR=2 angle gather with 3000 generations using the GA(top),QGA(middle)and HQGA(bottom),respectively.

Table6 REs and CCs between the inversion results and the true models using GA with different SNR.

Table7 REs and CCs between the inversion results and the true models using QGA with different SNR.

Table8 REs and CCs between the inversion results and the true models using HQGA with different SNR.

The inversion results and absolute errors using 8000 generations are shown in Fig.6,and their REs and CCs are calculated and displayed in Table 5.Inversion using the HQGA,QGA and GA ran for a period of approximately 25 h 25 min,25 h 11 min and 30 h 35 min,respectively.Further analysis showed that although the inversion results of the QGA and GA can be improved by increasing the generation number,their associated errors remained larger than those of HQGA,which gives optimal solutions with 2000 generations(Fig.6(d)-(e)).The QGA and GA need more generations to achieve a similarly accurate result of the HQGA,thus,present a huge computation cost.

Table5 REs and CCs between inversion results and the true models using GA,QGA and HQGA with 8000 generations.

Then three methods are further tested using data contaminated with different amount of noise(no noise,SNR=4 and 2).The iteration number is set as 3000.The inversion results and errors using the GA,QGA,and HQGA are shown in Figs.7-9,respectively.The results of noise-free data using the HQGA are closely consistent with real model parameters,i.e.,the global solutions,when the number of generation is 3000.For the two other noise levels,using the three algorithms,the differences between the observed and predicted models are proportional to the magnitude of the noise.The errors are smallest for the HQGA and largest for the GA.The errors associated with the HQGA are closer to zero compared with those of the other two methods.It is also deduced from the results that P-wave can be accurately inverted for in the noise-free case and the errors of inverted parameters significantly increase with increase in the amount of noise.The REs and CCs of the results are calculated and displayed in Tables 6-8.In the case of SNR=2,the relative errors become larger while correlation coefficients become smaller.

The convergence performance graphs are used to show the average error performance over all runs.Convergence performance of the GA,QGA and HQGA are investigated by plotting a graph of fitness(the square of the data residual)varied with generation.For the case of noise-free data,although HQGA has higher convergence compared with QGA,the fitness of both QGA and HQGA can reach a small value after 2000 generations.On the other hand,the fitness of GA remains quite large even after 8000 generations(Fig.10a).The differences in the results of the GA,QGA and HQGA increase with decrease in SNR.As shown in Fig.10(b)and(c),the fitness using HQGA can reach a small value after 3500 and 5000 generations for SNR=4 and 2,respectively.Even after 8000 generations,the QGA and GA do not achieve this small fitness value.These results clearly rank the HQGA higher than the QGA and GA on the basis of convergence performance and accuracy.In summary,the HQGA has performed better in terms of accuracy,efficiency,and robustness compared with the other two algorithms,and in practice,we recommend the Bayesian inference to ensure stability of the density inversion.

3.2.Field seismic data application

Fig.10.Convergence performance of the GA(green),QGA(blue)and HQGA(red)are compared under different noise levels((a)no noise,(b)SNR=4 and(c)SNR=2).Fitness is calculated with objective function to assess how close a solution is to achieve the set aims.

Fig.11.Constant-angle seismic profile of(a)10?,(b)20?,(c)30?,(d)40?and(e)50?,and well location.(f)The angle gather near the well location at CDP 400.(g)Five wavelets extracted from the constant-angle seismic profiles.

Fig.12.Inversion results of field data using the HQGA:(a)P-wave velocity,(b)S-wave velocity,and(c)density.

We applied the new improved algorithm to real seismic dataset from the Sichuan Basin in southwest China.The constant-angle seismic profiles(10oto 50o)are used as the input data for the inversion as shown in Fig.11(a)-(e).The angle gatherat CDP 400 near well location is displayed in Fig.11(f).It consists of five traces corresponding to 10oto 50o.Five seismic wavelets are extracted from the constant-angle seismic profiles and are used in the inversion(Fig.11(g)).Similar as the synthetic test,a dynamic search window acting as an initial model is constructed using the interpreted horizons and smoothed well logscalibratedwith theseismic data.With the assistance of dynamic time windows,the global solutions can be determined with higher resolution and faster convergence speed.Seismic inversion results using HQGA with 3000 generations are shown in Fig.12(a)-(c).The results are consistent with the real logs(Fig 13(a)).REs of HQGA are 0.0104,0.0476 and 0.0466 for P-wave velocity,S-wave velocity and density,respectively,and the CCs are 0.9267,0.9199 and 0.7237,respectively.In order to validate the superiority of HQGA in field data application,the real well logs are also compared with the inversion results of QGA and GA with 3000 generations as shown in Fig.13(d)-(f)and 13(g)-(i).The REs of QGA are 0.0214,0.0740 and 0.0861 for P-wave velocity,S-wave velocity and density,respectively,and the CCs of QGA are 0.8222,0.7569 and 0.5178,respectively.The REs of GA are 0.0416,0.1040 and 0.1245 for P-wave velocity,S-wave velocity and density,respectively,and the CCs of GA are 0.8195,0.7464 and 0.5060,respectively.The inverted results using the HQGA have the highest accuracy,while those of GA have the lowest accuracy.The inverted P-wave velocity using QGA and GA have acceptable accuracy,while the corresponding inverted density have great disagreement with the real log.Finally,the inverted parameters of HQGA are used to generate synthetic angle gathers for comparison with the real seismic gathers(Fig.14(a)-(c)).The results suggest that HQGA is a suitable method for the real data application.Nonlinear inversion using HQGA based on the Zoeppritz equation can effectively characterize the variation of complex formations and responses.To improve the practical applicability of the HQGA,extreme care must be taken in setting up the objective function as the objective loss is closely related to noise.The stop condition is also very important,considering its role in determining the trade-off between total computation time and accuracy.In general,robust data processing procedures preserving the amplitude information play an important role in improving the accuracy of inversion.

Fig.13.Comparison between real logs(black solid)and inversion results using HQGA(red solid),QGA(blue solid)and GA(green solid).The upper and lower bounds of search windows are denoted by grey curves.(a),(d),(g)P-wave velocity,(b),(e),(h)S-wave velocity and(c),(f),(i)density are shown from left to right in each inversion result.

Fig.14.(a)The real seismic angle gather at CDP 400.(b)Synthetic angle gather generated using the inverted parameters.(c)Difference between the real seismic angle gather and synthetic angle gather.

4.Discussion and conclusion

In this paper,we propose the HQGA for nonlinear AVO inversion method based on the exact Zoeppritz equation.To demonstrate its improvements over other similar algorithms,we applied it to synthetic and real field data.The results have clearly established the superiority of the new inversion method which can be broadly attributed to the relatively high accuracy of the exact Zoeppritz equation and,strong global search ability and high efficiency of the HQGA.The self-adaptive strategy helps the HQGA to achieve the optimum solution efficiently by updating the search step,and the quantum mutation and quantum catastrophe further improved the local and global search capabilities.Dynamic search windows are used in all of the inversions and can be constructed by combining the interpreted horizons and smoothed well log data calibrated with seismic data.The inversion using the dynamic search windows have higher efficiency than that of the static search windows.Even though it is advised to set the range of the search window with extreme care,in general,a large search window requires more population and has a better likelihood of obtaining the global solution,while a small search window may miss the true solution.In this study,the minimum and maximum values of the dynamic search window were set to±20%,according to the measured well logs.With fewer searches,nonlinear inversion using HQGA can efficiently achieve higher accuracy than the GA and QGA.The proposed HQGA inversion method can be also useful in other geophysics problems and disciplines.Relevant studies considering anisotropy and complex wave propagation effects will be considered in the future work.

Acknowledgements

This work was supported by the National Natural Science Foundation of China(U19B6003,42122029)and the Strategic Cooperation Technology Projects of CNPC and CUPB(ZLZX 2020-03).Jiwei Cheng is partially supported by SEG/WesternGeco Scholarship,SEG Foundation/Chevron Scholarship,and SEG/Norman and Shirley Domenico Scholarship.We gratefully acknowledge the helpful comments from the editors and anonymous reviewers,which greatly improved this manuscript.