Earthquake quality factor acquisition method, system and device and storage medium

Abstract

The invention discloses an earthquake quality factor acquisition method, system and device and a storage medium, and the method comprises the steps: S1, obtaining an earthquake observation record, a speed parameter and an initial quality factor model, and building an earthquake quality factor model; s2, establishing a trust region type least square objective function according to the wave field residual error; s3, updating the trust region type least square objective function by adopting the optimized damping coefficient; s4, updating the seismic quality factor model, and performing iterative calculation of S2-S3 until the trust region type least square objective function tends to be stable or the maximum number of iterations, thereby obtaining an updated seismic quality factor model; and S5, setting a plurality of frequency groups from low frequency to high frequency, carrying out iterative calculation from S2-S4, and calculating to obtain the seismic quality factor by adopting the seismic quality factor model under the maximum frequency group. The calculation amount can be reduced, the modeling efficiency can be improved, and the modeling precision of the seismic quality factor can be maintained.

Description

Technical Field

[0001] The invention belongs to the field of geophysical exploration and relates to a method, system, equipment and storage medium for obtaining a seismic quality factor. Background Art

[0002] When seismic waves propagate in the stratum, the viscoelasticity and other factors cause energy attenuation and loss of high-frequency information, making seismic inversion and interpretation difficult. Accurate modeling of the formation quality factor Q value can not only accurately describe the formation viscoelasticity, but also compensate for the high-frequency information lost in the propagation of seismic waves. The traditional Q value modeling method based on seismic data mainly estimates the Q value by selecting a suitable time window, extracting wavelets, and then estimating the Q value by comparing the spectral change relationship between the two wavelets, such as the spectral ratio method (Bath, 1974), the centroid frequency shift method (Quan and Harris, 1997), and the peak frequency method (Zhang, 2002). The inversion effect of such methods is limited by how to choose the best time window type and length. As an effective seismic attribute modeling technology, waveform inversion technology can solve nonlinear problems well. By introducing complex velocity to establish the frequency domain viscous acoustic wave equation to simulate seismic records with absorption attenuation effects, the waveform inversion method is used to extract the Q value, which can avoid the selection of time windows and eliminate the influence of scattering attenuation on Q value estimation. However, for nonlinear inversion problems, due to their high nonlinearity and pathological nature, traditional methods such as the steepest descent method and Newton's method are all local optimization algorithms, and their core strategies are to k , calculate the descending direction d k , and then from x k Departure along d k Direction adjustment step size, linear search x k+1 These methods are prone to causing the search results of the objective function to fall into a local minimum rather than a global optimal solution.

[0003] The trust region algorithm can be traced back to the solution of nonlinear least squares problems proposed by Levenberg in 1944 and Marquardt in 1963. It is also called damped least squares method (LM algorithm). Its strategy is to give an initial iteration point x k , establish a x k Centered on R k The search radius R is a spherical area whose range is called the trust region. Then the search radius R is continuously optimized by solving subproblems. kUntil the accuracy requirement is met. The trust region algorithm is equivalent to giving the descent direction and step size at the same time. Compared with the linear search strategy, it has stronger convergence and can significantly improve the ill-posed problem of the local optimization algorithm. Yuan (2000) and Wang and Yuan (2006) gave the convergence proof of the trust region algorithm. In the field of geophysics, ZengQi et al. used the trust region algorithm for magnetotelluric inversion of layered media; Ju-Wonoh et al. (2013) applied the damped least squares method to full waveform inversion in the frequency domain; Zhang Huan (2017) also constructed a subspace adaptive trust region algorithm for full waveform inversion.

[0004] Although the traditional damped least squares method has good convergence and robustness, solving sub-problems often requires iterative calculations through the quasi-Newton method, which consumes a lot of computing resources and reduces the overall computing efficiency. Summary of the invention

[0005] The purpose of the present invention is to overcome the shortcomings of the above-mentioned prior art and provide a method, system, device and storage medium for obtaining earthquake quality factor, which can reduce the amount of calculation, improve the modeling efficiency and maintain the modeling accuracy of the earthquake quality factor.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for obtaining a seismic quality factor includes the following steps:

[0008] S1, obtaining earthquake observation records, velocity parameters and an initial quality factor model, and establishing an earthquake quality factor model based on the observation records, velocity parameters and the initial quality factor model;

[0009] S2, forward modeling the earthquake quality factor model to obtain simulated earthquake records, subtracting the simulated earthquake records from the earthquake observation records to obtain the wave field residual, and establishing the trust region type least squares objective function based on the wave field residual;

[0010] S3, using the optimized damping coefficient to update the trust region type least squares objective function;

[0011] S4, updating the earthquake quality factor model, and performing iterative calculations of S2-S3 until the trust region type least squares objective function tends to be stable or the maximum number of iterations is reached, at which time the updated earthquake quality factor model is obtained;

[0012] S5, setting multiple frequency groups from low frequency to high frequency, performing iterative calculations of S2-S4, and using the earthquake quality factor model under the maximum frequency group to calculate the earthquake quality factor.

[0013] Preferably, the viscoacoustic wave equation of the two-dimensional frequency domain additional perfect matching layer of the seismic quality factor model is:

[0014]

[0015] Where x is the horizontal direction; z is the vertical direction; ω is the circular frequency; K(x,z) is the bulk modulus; ρ(x,z) is the density; is the complex velocity, Q(x,z) is the quality factor; v(x,z) = v x (x,z)+v z (x,z), v(x,z) is the velocity, v x (x,z),v z (x,z) is the component of v(x,z) at x and z; U(x,z,ω)=U x (x,z,ω)+U z (x,z,ω), U(x,z,ω) is the stress wave field in the frequency domain, U x (x,z,ω),U z (x,z,ω) is the component of U(x,z,ω) at x and z; S(x,z,ω) is the frequency domain source term;

[0016]

[0017] Where g x (x,z),g z (x,z) is the attenuation function of the perfectly matched layer in the x and z directions; C PML is the absorption coefficient; L is the thickness of the perfectly matched layer.

[0018] Furthermore, the viscoacoustic wave equation with an additional perfectly matched layer in the two-dimensional frequency domain is established based on the finite difference equation of the 9-point hybrid staggered grid difference format:

[0019]

[0020] Where a = 0.5461, c = 0.6248, d = 0.09381 are the differential coefficients a priori through the optimization algorithm; Δ = Δx = Δz is the spatial grid spacing; U x-1,z-1 ,U x,z-1 ,U x+1,z-1 ,…etc. are the frequency domain wave fields at the point (x,z) and its eight adjacent points.

[0021] Preferably, the trust region type least squares objective function is:

[0022]

[0023] In the formula T stands for transpose; * represents conjugate; δU represents the wavefield residual; m represents the model parameter.

[0024] Preferably, the specific process of S3 is: converting the problem of finding the minimum value of the trust region type least squares objective function into a sub-problem:

[0025]

[0026] Where B k is the approximation matrix B of the Hessian matrix k =(J k ) T J k , if B k is the trust region type least squares objective function C(m) at m k The Hessian matrix at D k The error of (d) is o||d|| 3 , R k represents the trust region radius, and the subproblem is updated by optimizing the damping coefficient. The trust region range and the objective function gradient are calculated according to the updated trust region type least squares objective function.

[0027] The above subproblems are solved using the optimality conditions of the constrained optimization problem, which are:

[0028]

[0029] In the formula represents the gradient of the trust region least squares objective function; λ k Represents the damping coefficient.

[0030] Further, in solving the subproblem, we get the descending direction d k After that, can the approximate solution of the original unconstrained optimization least squares problem be updated to m k +d k , it is also necessary to judge D by the ratio of the actual decrease in the objective function value to the estimated decrease k (d) Approximation:

[0031]

[0032] According to ρ k The size of λ is adjusted k , if ρ k Close to 1, reduce λ k , trust region radius R k Increase, otherwise λ should be increased k , trust region radius R k Reduce, so that each step iterates to solve a sub-problem once, and then dynamically adjust λ according tok Get the new descent direction d and step size.

[0033] Preferably, the specific process of S4 is: setting multiple frequency groups of different sizes respectively, performing multiple iterations in each frequency group in order of size, and using the low-frequency inversion result as the initial model of the high-frequency inversion, and when the maximum frequency is reached, outputting the inversion result to obtain the final seismic quality factor model.

[0034] A seismic quality factor acquisition system, comprising:

[0035] A model building module, used for obtaining earthquake observation records, velocity parameters and an initial quality factor model, and building an earthquake quality factor model based on the observation records, velocity parameters and the initial quality factor model;

[0036] A function building module is used to forward model the earthquake quality factor model to obtain simulated earthquake records, subtract the simulated earthquake records from the earthquake observation records to obtain the wave field residual, and establish a trust region type least squares objective function based on the wave field residual;

[0037] An updating module, used for updating the trust region type least squares objective function by using an optimized damping coefficient;

[0038] The optimization iteration module is used to update the earthquake quality factor model and perform iterative calculations of the function establishment module and the update module until the objective function value of the trust region type least squares objective function tends to be stable or the maximum number of iterations is reached, at which time the updated earthquake quality factor model is obtained;

[0039] The frequency iteration module is used to set multiple frequency groups from low frequency to high frequency, perform iterative calculations of the function establishment module, update module and optimization iteration module, and use the earthquake quality factor model under the maximum frequency group to calculate the earthquake quality factor.

[0040] A computer device comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the seismic quality factor acquisition method when executing the computer program.

[0041] A computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the seismic quality factor acquisition method are implemented.

[0042] Compared with the prior art, the present invention has the following beneficial effects:

[0043] The method of the present invention is different from the local optimization algorithm used in general domestic earthquake attribute modeling, such as the gradient method or the Newton method. The initial model is first used to find the descent direction, and then the line search criterion of the step length is given. The method of the present invention introduces the trust region damped least squares method into the field of earthquake quality factor Q value modeling for the first time. The algorithm directly solves the approximate model in a bounded area by giving the descent direction and step length of the initial model at the same time, and then iterates continuously, which solves the problem of serious dependence on the initial model selection and easy to fall into local extreme values ​​in the current general domestic earthquake attribute modeling, and improves the stability and accuracy of modeling.

[0044] Furthermore, in the classical trust region damped least squares method, each optimization search radius R k A sub-problem needs to be solved every time, and this solution process will occupy a considerable amount of calculation. The method of the present invention dynamically adjusts the trust region radius R by directly changing the damping coefficient. k The size of makes it necessary to solve only one subproblem in each iteration. This strategy can reduce the amount of calculation and improve the modeling efficiency while maintaining the original modeling accuracy of the seismic quality factor Q value. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a flow chart of a method for modeling seismic quality factor based on optimizing damping coefficient λ provided in an embodiment;

[0046] Figure 2 is a speed parameter model provided by an embodiment;

[0047] Figure 3 is an initial seismic quality factor model provided by an embodiment;

[0048] Figure 4 is an observation record for seismic quality factor inversion provided by an embodiment;

[0049] Figure 5 is a schematic diagram of a trust region subproblem provided in an embodiment;

[0050] Figure 6 The seismic quality factor model provided in the embodiment is calculated by using the trust region least square method based on the optimized damping coefficient λ;

[0051] Figure 7 It is a quality factor Q value curve extracted from the earthquake quality factor calculation model at x=500m provided in the embodiment. DETAILED DESCRIPTION

[0052] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments; based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0053] It should be noted that the words "front", "rear", "left", "right", "up" and "down" used in the following description refer to directions in the drawings, and the words "inside" and "outside" refer to directions toward or away from the geometric center of a specific component, respectively.

[0054] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as those generally understood by those skilled in the art of the present invention. The terms used herein in the specification of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The term "and / or" used herein includes any and all combinations of one or more related listed items.

[0055] The seismic quality factor acquisition method of the present invention adopts the following technical scheme: establish a two-dimensional frequency domain viscoacoustic medium wave equation containing the seismic quality factor Q (absorption attenuation), and then perform seismic numerical simulation based on a 9-point hybrid staggered grid difference format and a perfectly matched layer (PML) absorption boundary condition to reduce grid dispersion and error anisotropy and improve the accuracy of forward simulation. The objective function gradient and trust domain range are estimated by the trust domain type least squares method based on the optimization damping coefficient λ update, solving the problem that the local optimization algorithm is prone to fall into the local extreme value along the linear search direction, and ensuring the modeling accuracy of the seismic quality factor Q value.

[0056] The present invention optimizes the inversion strategy to reduce the search radius R k The problem is transformed into updating the damping coefficient λ k The problem of the damped least squares method is solved so that each iteration only needs to solve the sub-problem once, thereby improving the problem of excessive computational complexity when searching for optimal solutions to sub-problems in the traditional trust region damped least squares method, thereby improving computational efficiency without affecting accuracy. In addition, the damped least squares method, as an algorithm with excellent robustness, can also improve the problem of singular and non-positive definite Hessian matrix.

[0057] like Figure 1 As shown, the method includes:

[0058] (i) Actual earthquake observation records, velocity parameters and initial quality factor model are input and the viscoacoustic wave equation and 9-point hybrid staggered grid difference format of the two-dimensional frequency domain with an additional perfectly matched layer (PML) are established.

[0059] The seismic observation records of the actual formation quality factor required for seismic quality factor modeling, velocity parameters and the smoothed initial seismic quality factor model are input into the seismic quality factor Q value modeling program, and the seismic quality factor model is established. The viscoacoustic wave equation of the two-dimensional frequency domain of the seismic quality factor model with a perfectly matched layer (PML) is:

[0060]

[0061] Where x is the horizontal direction; z is the vertical direction; ω is the circular frequency; K(x,z) is the bulk modulus; ρ(x,z) is the density; is the complex velocity, Q(x,z) is the quality factor; v(x,z) = v x (x,z)+v z (x,z), v(x,z) is the velocity, v x (x,z),v z (x,z) is the component of v(x,z) at x and z; U(x,z,ω)=U x (x,z,ω)+U z (x,z,ω), U(x,z,ω) is the stress wave field in the frequency domain, U x (x,z,ω),U z (x,z,ω) is the component of U(x,z,ω) at x and z; S(x,z,ω) is the frequency domain source term;

[0062]

[0063] Where g x (x,z),g z (x,z) is the attenuation function of the perfectly matched layer (PML) in the x and z directions; C PML is the absorption coefficient; L is the thickness of the perfectly matched layer.

[0064] The finite difference equation based on the 9-point hybrid staggered grid difference format is established for formula (1):

[0065]

[0066] Where a = 0.5461, c = 0.6248, d = 0.09381 are the differential coefficients a priori through the optimization algorithm; Δ = Δx = Δz is the spatial grid spacing; U x-1,z-1 ,U x,z-1 ,U x+1,z-1 ,…etc. are the frequency domain wave fields at the point (x,z) and its eight adjacent points.

[0067] (ii) Obtain the wave field residual and establish the objective function.

[0068] The seismic quality factor model is forward modeled using formula (2) to obtain simulated seismic records. The wave field residual is obtained by subtracting the simulated seismic records from the seismic observation records, and the L2 norm is used to establish the objective function in the least squares sense without constraints.

[0069]

[0070] In the formula T stands for transpose; * represents conjugate; δU represents the wavefield residual; m represents the model parameter.

[0071] (iii) Estimate the trust region range and objective function gradient by using the trust region least squares method based on optimizing the damping coefficient λ.

[0072] Using the second-order approximation of C(m) to characterize C(m) in m k Features of the place:

[0073]

[0074] Where B k is the approximation matrix B of the Hessian matrix k =(J k ) T J k , if B k It happens that the objective function C(m) is k The Hessian matrix at , then D k The error of (d) is o||d|| 3 .

[0075] Due to D k (d) describes the function C(m) at m k The local properties at +d are to make D k (d) can accurately represent the objective function, so it is necessary to limit the size of d, so the following provisions are made:

[0076] Ω k ={m k +d|||d||≤R k} (5)

[0077] Where Ω k represents the trust region; R k represents the trust region radius.

[0078] In this domain D k (d) can well approximate C(m k+d), so the problem of finding the minimum value of the function C(m) is transformed into the following sub-problems:

[0079]

[0080] To solve the above sub-problems using the optimality conditions of the constrained optimization problem, the following conditions need to be met:

[0081]

[0082] In the formula represents the gradient of the objective function; k Represents the damping coefficient.

[0083] From formula 7, we can deduce that when B k +λ k When I is reversible, we can get d(λ k )=-(B k +λ k I) -1 g k , that is, the optimal solution of d can be expressed as a family of vectors with λ as a parameter, so we only need to find a suitable λ so that the second and third equations in formula 7 are valid; through analysis, we can find that ||d k ||=R k The solution exists and is unique, so we only need to use Newton's method to iteratively solve the optimal solution λ, and then according to d(λ k )=-(B k +λ k I) -1 g k The iteration direction can be found.

[0084] In solving the subproblem, we get the descending direction d k After that, can the approximate solution of the original unconstrained optimization least squares problem be updated to m k +d k , it is also necessary to judge D by the ratio of the actual decrease in the objective function value to the estimated decrease k (d) Approximation:

[0085]

[0086] In the classical trust region damped least squares method, we can use ρ k The value of dynamically adjusts the trust region radius R k , if ρ k Close to 1, indicating that D k (d) It approximates C(m k ), then the trust region radius R is enlarged k , if ρ k Very small or even negative, indicating that D k(d) It is not possible to describe C(m k ), at this time, reduce the trust region radius R k This mechanism can be used to dynamically adjust R k , so that the approximate model D k (d) can well describe C(m k ) characteristics. In the classical trust region damped least squares method, each optimization search radius R k Each time a sub-problem needs to be solved, this solution process will take up a considerable amount of computation.

[0087] From the previous inference, we can know that d(λ k )=-(B k +λ k I) -1 g k , and ||d(λ)|| 2 is a strictly decreasing function of λ, that is, the modulus length of the descending direction d decreases with the increase of λ, which means that the trust region radius R k is implicitly determined by λ, and directly changing the size of λ is equivalent to adjusting R k Therefore, the present invention uses this relationship to establish a Q value modeling optimization strategy based on optimizing the damping coefficient. k The size of λ is adjusted k , in contrast to the classical trust region method, if ρ k Close to 1, λ should be reduced k , then the trust region radius R k Increase, otherwise λ should be increased k , then the trust region radius R k This makes it possible to solve the subproblem only once in each iteration, and then dynamically adjust λ according to the optimization strategy of the present invention. k The new descent direction d and step size can be obtained, thereby improving the computational efficiency while maintaining the modeling accuracy.

[0088] (iv) Update the earthquake quality factor model under the current frequency group.

[0089] Repeat steps (iii) to (iv) and perform iterative calculations until the trust region type least squares objective function becomes stable or reaches the maximum number of iterations, at which time the seismic quality factor model under the current frequency group is obtained.

[0090] (V) Frequency group updating: using the updated earthquake quality factor model as the initial model to obtain the earthquake quality factor.

[0091] Specifically, multiple frequency groups of different sizes from low frequency to high frequency are set respectively, and the updated seismic quality factor model at low frequency is used as the initial model of the next frequency group. Steps (ii) to (iv) are repeated until the maximum frequency is reached, and the inversion results are output to obtain the final seismic quality factor model.

[0092] Specifically, multiple frequency groups of different sizes are set, and multiple iterations are performed in each frequency group in order of size. The low-frequency inversion results are used as the initial model for high-frequency inversion. When the iterative relative error is less than 0.01% or the maximum frequency is reached, the inversion results are output to obtain the final seismic quality factor model.

[0093] The process of obtaining the earthquake quality factor of the present invention is demonstrated below based on actual cases.

[0094] Step (i), input observation records, velocity parameters and initial seismic quality factor model.

[0095] Step (ii), obtain the wave field residual and establish the objective function.

[0096] Step (iii), establish a trust region subproblem, and calculate the objective function gradient and the trust region range by using the trust region least squares method based on the optimized damping coefficient λ.

[0097] Step (iv), updating the earthquake quality factor model.

[0098] Step (five), using the updated earthquake quality factor model to obtain the earthquake quality factor.

[0099] In the specific implementation, step (i) inputs observation records, velocity parameters and an initial seismic quality factor model, including:

[0100] like Figure 2 The figure shows the velocity parameter model for inversion input. The model grid size is 508×345, and the velocity range is from 2000 m·s -1 To 3500m·s -1 .

[0101] like Figure 3 Shown is the initial seismic quality factor model of the inversion input. The model grid number is also 508×345, and the quality factor ranges from 52.4 to 90.2.

[0102] like Figure 4The figure shows the seismic records obtained in the time domain by the finite difference real earthquake quality factor model. The specific parameters of the differential format are: the source uses the Ricker wavelet with a main frequency of 40Hz, the boundary uses the perfectly matched layer (PML) with 50 absorption layers, the calculation accuracy is 2nd order in time and 8th order in space, the grid interval is 5m, the time interval is 0.1ms, and the forward modeling time is 1.5s. Since the inversion is performed in frequency, the observation records in the time domain need to be transformed into the frequency domain through fast Fourier transform during the specific implementation.

[0103] In the specific implementation, step (ii) obtains the wave field residual and establishes the objective function, including: subtracting the simulated record of the frequency domain forward modeling of the seismic quality factor model from the observed record to obtain the wave field residual; establishing the objective function, the L1 norm, L2 norm, Huber norm and mixed norm can be used. Since the L2 norm is relatively simple in form and its derivative is continuous, the L2 norm is used in this case to establish the least squares objective function under no constraints.

[0104] In the specific implementation, step (iii) establishes a trust region subproblem, and calculates the trust region range and the objective function gradient by using a trust region damped least squares method based on an optimized damping coefficient λ, including: k )The problem of minimization is transformed into solving the following trust region subproblem:

[0105]

[0106] The conditions to be met are:

[0107]

[0108] Where D k (d) represents the sub-problem, B k is the approximation matrix of the Hessian matrix, d = mm k is the descending direction, λ is the damping coefficient, Represents the gradient direction of the objective function. Solving the above subproblems yields d(λ k )=-(B k +λ k I) -1 g k .

[0109] like Figure 5 The figure is a schematic diagram of the trust region subproblem provided by an embodiment of the present invention. The solid line in the figure represents the contour line of C(m), and the dotted line represents D k (d) The contour lines, represents the descent direction obtained by solving the subproblem, Represents the Newton direction. It can be seen that the trust region algorithm limits the size of d, so that the iteration effect can be improved when the Newton direction is poor.

[0110] Solve the subproblem to get the descent direction d k Then, the decline rate ρ is calculated by the ratio of the actual decline of the function value to the estimated decline k , and then according to ρ k Size judgment D k (d) approximation, changing the damping coefficient λ k , dynamically adjust the trust region radius R k .

[0111]

[0112] In this case, we set the accuracy parameters to 0≤ρ1<ρ2<1,β1<1<β2, where ρ1=0.25,ρ2=0.75,β1=0.25,β2=1.5. k >ρ2, then reduce the damping coefficient λ k+1 =β1λ k (Expand the radius of the trust region); if ρ1<ρ k <ρ2, the damping coefficient remains unchanged (the trust region radius remains unchanged); otherwise, let m k+1 =m k , expand the damping coefficient λ k+1 =β2λ k (Reduce the trust region radius) and expand the damping coefficient λ k+1 Substitute it into the subproblem (Equation (2)) and calculate the new descent direction d k and the drop rate ρ k , repeat the above process.

[0113] This strategy makes each optimization search radius R k Each time, it is only necessary to solve the subproblem once, and then it is only necessary to dynamically adjust λ according to the optimization strategy of the present invention. k The new descending direction d can be obtained k and step length.

[0114] like Figure 6 The figure shows the seismic quality factor modeling result obtained by using the trust region type least square method based on the optimized damping coefficient λ provided by the embodiment of the present invention, and setting the convergence condition to an iteration relative error of less than 0.01% and a maximum number of iterations of 80. It can be seen that the update direction and step size do not fall into the local extreme value due to the constraint of the trust region, accurately restore the structural information of the model, and the modeling effect is excellent.

[0115] like Figure 7The diagram shows the trust region least square method based on the optimized damping coefficient λ and the classical damped least square method provided by the embodiment of the present invention, and extracts the Q value curve at x=500m. The diagram shows that the two different methods are basically consistent in the final modeling effect. At the same time, it can be seen that the quality factor modeling method based on the optimized damping coefficient λ of the present invention has smaller overall fluctuations and better convergence effect.

[0116] The hardware parameters used in this example are: 11th Gen Intel(R)Core(TM)i7-11700K CPU@3.60GHz (8 cores and 16 threads), RAM32GB. The forward modeling parameters are set the same, and the termination condition is set to the iteration relative error less than 0.01%, and the maximum number of iterations is 80. The statistical modeling calculation time of the trust region least squares method based on the optimized damping coefficient λ, the traditional damped least squares method, and the steepest descent method are 691s, 978s, and 1123s respectively. Compared with the traditional damped least squares method, the modeling method of this case has a 29% improvement in computational efficiency, and compared with the steepest descent method, it has a 38% improvement in computational efficiency.

[0117] In specific implementation, step (iv) updates the earthquake quality factor model under the current frequency group, including:

[0118] In this case, the maximum number of iterations for each frequency group is set to 20 times, and the iteration relative error is 0.01%. Steps (iii) to (iv) are repeated. When the trust region type least squares objective function meets the maximum number of iterations or convergence conditions, the seismic quality factor model for the current frequency group is obtained.

[0119] In the specific implementation, step (5) frequency group updating, using the updated earthquake quality factor model as the initial model, and obtaining the earthquake quality factor includes:

[0120] In this case, the frequency groups are set to 5-10 Hz, 10-15 Hz, 15-20 Hz, and 20-25 Hz, and the low-frequency inversion results are used as the initial model of the high-frequency inversion in turn, and the two-dimensional frequency domain viscoacoustic medium forward modeling is performed to obtain new observation records, and steps (ii) to (iv) are repeated until the maximum frequency is reached, and the inversion results are output to obtain the seismic quality factor model. The trust region type least squares method based on the optimization of the damping coefficient λ of the present invention can achieve excellent modeling results while reducing the amount of calculation.

[0121] The following are device embodiments of the present invention, which can be used to perform method embodiments of the present invention. For details not disclosed in the device embodiments, please refer to the method embodiments of the present invention.

[0122] In yet another embodiment of the present invention, a seismic quality factor acquisition system is provided. The seismic quality factor acquisition system can be used to implement the above-mentioned seismic quality factor acquisition method. Specifically, the seismic quality factor acquisition system includes a model building module, a function building module, an update module, an optimization iteration module and a frequency iteration module.

[0123] Among them, the model building module is used to obtain seismic observation records, velocity parameters and initial quality factor models, and build a seismic quality factor model based on the observation records, velocity parameters and initial quality factor models.

[0124] The objective function establishment module is used to obtain simulated seismic records by forward modeling the seismic quality factor model, to obtain wave field residuals by subtracting simulated seismic records from seismic observation records, and to establish a trust region type least squares objective function based on the wave field residuals.

[0125] Update module, used to update the trust region type least squares objective function using the optimized damping coefficient.

[0126] The optimization iteration module is used to update the earthquake quality factor model and perform iterative calculations of the function establishment module and the update module until the objective function value of the trust region type least squares objective function tends to be stable or the maximum number of iterations is reached, at which time the updated earthquake quality factor model is obtained.

[0127] The frequency iteration module is used to set multiple frequency groups from low frequency to high frequency, perform iterative calculations of the function establishment module, update module and optimization iteration module, and use the earthquake quality factor model under the maximum frequency group to calculate the earthquake quality factor.

[0128] In another embodiment of the present invention, a terminal device is provided, the terminal device includes a processor and a memory, the memory is used to store a computer program, the computer program includes program instructions, and the processor is used to execute the program instructions stored in the computer storage medium. The processor may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc. It is the computing core and control core of the terminal, which is suitable for implementing one or more instructions, and is specifically suitable for loading and executing one or more instructions to implement the corresponding method flow or corresponding function; the processor described in the embodiment of the present invention can be used for the operation of the seismic quality factor acquisition method, including: S1, obtaining seismic observation records, velocity parameters and initial quality factor models, and establishing a seismic quality factor model based on the observation records, velocity parameters and initial quality factor models; S2, The simulated seismic record is obtained by forward modeling, and the wave field residual is obtained by subtracting the simulated seismic record from the seismic observation record. The trust region type least squares objective function is established according to the wave field residual; S3, the trust region type least squares objective function is updated by optimizing the damping coefficient; S4, the seismic quality factor model is updated, and the iterative calculation of S2-S3 is performed until the objective function value of the trust region type least squares objective function tends to be stable or the maximum number of iterations is reached, at which time the updated seismic quality factor model is obtained; S5, multiple frequency groups are set from low frequency to high frequency, and the iterative calculation of S2-S4 is performed. The seismic quality factor is calculated using the seismic quality factor model under the maximum frequency group.

[0129] In another embodiment, the present invention further provides a computer-readable storage medium (Memory), which is a memory device in a terminal device for storing programs and data. It is understandable that the computer-readable storage medium here can include both the built-in storage medium in the terminal device and the extended storage medium supported by the terminal device. The computer-readable storage medium provides a storage space, which stores the operating system of the terminal. In addition, one or more instructions suitable for being loaded and executed by the processor are also stored in the storage space, and these instructions can be one or more computer programs (including program codes). It should be noted that the computer-readable storage medium here can be a high-speed RAM memory or a non-volatile memory, such as at least one disk memory.

[0130] One or more instructions stored in a computer-readable storage medium can be loaded and executed by a processor to implement the corresponding steps of the method for obtaining the seismic quality factor in the above embodiment; one or more instructions in the computer-readable storage medium are loaded by the processor and the following steps are executed: S1, obtaining seismic observation records, velocity parameters and an initial quality factor model, and establishing a seismic quality factor model based on the observation records, velocity parameters and the initial quality factor model; S2, forward modeling the seismic quality factor model to obtain simulated seismic records, subtracting the simulated seismic records from the seismic observation records to obtain wave field residuals, and establishing a trust region type least squares objective function according to the wave field residuals; S3, updating the trust region type least squares objective function by optimizing the damping coefficient; S4, updating the seismic quality factor model, performing iterative calculations of S2-S3 until the objective function value of the trust region type least squares objective function tends to be stable or reaches the maximum number of iterations, at which time the updated seismic quality factor model is obtained; S5, setting multiple frequency groups from low frequency to high frequency, performing iterative calculations of S2-S4, and using the seismic quality factor model under the maximum frequency group to calculate the seismic quality factor.

[0131] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application may adopt the form of a computer program product implemented in one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that include computer-usable program code.

[0132] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0133] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0134] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0135] It should be noted that, in this article, relational terms such as first and second, etc. are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device.

[0136] It should be understood that the above description is for illustration and not for limitation. Many embodiments and many applications beyond the examples provided will be apparent to those skilled in the art upon reading the above description. Therefore, the scope of the present teachings should not be determined with reference to the above description, but rather with reference to the foregoing claims and the full scope of equivalents to which such claims are entitled. For the purpose of comprehensiveness, all articles and references, including disclosures of patent applications and publications, are incorporated herein by reference. The omission of any aspect of the subject matter disclosed herein in the foregoing claims is not intended to be a waiver of such subject matter, nor should it be considered that the applicant has not considered such subject matter to be part of the disclosed inventive subject matter.

Claims

1. A method for obtaining an earthquake quality factor, characterized in that: The process includes: S1, obtaining earthquake observation records, velocity parameters and an initial quality factor model, and establishing an earthquake quality factor model based on the observation records, velocity parameters and the initial quality factor model; S2, forward modeling the earthquake quality factor model to obtain simulated earthquake records, subtracting the simulated earthquake records from the earthquake observation records to obtain the wave field residual, and establishing the trust region type least squares objective function based on the wave field residual; S3, using the optimized damping coefficient to update the trust region type least squares objective function; S4, updating the earthquake quality factor model, and performing iterative calculations of S2-S3 until the trust region type least squares objective function tends to be stable or the maximum number of iterations is reached, at which time the updated earthquake quality factor model is obtained; S5, setting multiple frequency groups from low frequency to high frequency, performing iterative calculations of S2-S4, and using the earthquake quality factor model under the maximum frequency group to calculate the earthquake quality factor.

2. The method for obtaining earthquake quality factor according to claim 1, characterized in that: The viscoacoustic wave equation of the two-dimensional frequency domain with an additional perfectly matched layer in the seismic quality factor model is: Where x is the horizontal direction; z is the vertical direction; ω is the circular frequency; K(x,z) is the bulk modulus; ρ(x,z) is the density; is the complex velocity, Q(x,z) is the quality factor; v(x,z) = v x (x,z)+v z (x,z), v(x,z) is the velocity, v x (x,z),v z (x,z) is the component of v(x,z) at x and z; U(x,z,ω)=U x (x,z,ω)+U z (x,z,ω), U(x,z,ω) is the stress wave field in the frequency domain, U x (x,z,ω),U z (x,z,ω) is the component of U(x,z,ω) at x and z; S(x,z,ω) is the frequency domain source term; Where g x (x,z),g z (x,z) is the attenuation function of the perfectly matched layer in the x and z directions; C PML is the absorption coefficient; L is the thickness of the perfectly matched layer.

3. The method for obtaining earthquake quality factor according to claim 2, characterized in that: The establishment of the viscoacoustic wave equation in the two-dimensional frequency domain with an additional perfectly matched layer is based on a finite difference equation with a 9-point hybrid staggered grid difference format: Where a = 0.5461, c = 0.6248, d = 0.09381 are the differential coefficients a priori through the optimization algorithm; Δ = Δx = Δz is the spatial grid spacing; U x-1,z-1 ,U x,z-1 ,U x+1,z-1 ,…etc. are the frequency domain wave fields at the point (x,z) and its eight adjacent points.

4. The method for obtaining earthquake quality factor according to claim 1, characterized in that: The trust region least squares objective function is: Where T stands for transpose; * represents conjugate; δU represents the wavefield residual; m represents the model parameter.

5. The method for obtaining earthquake quality factor according to claim 1, characterized in that: The specific process of S3 is: transform the problem of finding the minimum value of the trust region type least squares objective function into a sub-problem: Where B k is the approximation matrix B of the Hessian matrix k =(J k ) T J k , if B k is the trust region type least squares objective function C(m) at m k The Hessian matrix at D k The error of (d) is o||d|| 3 , R k represents the trust region radius, and the subproblem is updated by optimizing the damping coefficient. The trust region range and the objective function gradient are calculated according to the updated trust region type least squares objective function. The above subproblems are solved using the optimality conditions of the constrained optimization problem, which are: In the formula represents the gradient of the trust region least squares objective function; λ k Represents the damping coefficient.

6. The method for obtaining earthquake quality factor according to claim 5, characterized in that: In solving the subproblem, we get the descent direction d k After that, can the approximate solution of the original unconstrained optimization least squares problem be updated to m k +d k , it is also necessary to judge D by the ratio of the actual decrease in the objective function value to the estimated decrease k (d) Approximation: According to k The size of λ is adjusted k , if ρ k Close to 1, reduce λ k , trust region radius R k Increase, otherwise λ should be increased k , trust region radius R k Reduce, so that each step iterates to solve a sub-problem once, and then dynamically adjust λ according to k Get the new descent direction d and step size.

7. The method for obtaining earthquake quality factor according to claim 1, characterized in that: The specific process of S4 is as follows: set multiple frequency groups of different sizes, perform multiple iterations in each frequency group in turn according to the size, and use the low-frequency inversion results as the initial model of high-frequency inversion. When the maximum frequency is reached, output the inversion results to obtain the final seismic quality factor model.

8. A seismic quality factor acquisition system, characterized in that: include: A model building module, used for obtaining earthquake observation records, velocity parameters and an initial quality factor model, and building an earthquake quality factor model based on the observation records, velocity parameters and the initial quality factor model; A function building module is used to forward model the earthquake quality factor model to obtain simulated earthquake records, subtract the simulated earthquake records from the earthquake observation records to obtain the wave field residual, and establish a trust region type least squares objective function based on the wave field residual; An updating module, used for updating the trust region type least squares objective function by using an optimized damping coefficient; The optimization iteration module is used to update the earthquake quality factor model and perform iterative calculations of the function establishment module and the update module until the objective function value of the trust region type least squares objective function tends to be stable or the maximum number of iterations is reached, at which time the updated earthquake quality factor model is obtained; The frequency iteration module is used to set multiple frequency groups from low frequency to high frequency, perform iterative calculations of the function establishment module, the update module and the optimization iteration module, and use the earthquake quality factor model under the maximum frequency group to calculate the earthquake quality factor.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of the method for obtaining the seismic quality factor according to any one of claims 1 to 7 are implemented.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the steps of the method for obtaining seismic quality factor according to any one of claims 1 to 7 are implemented.

Images