Heterogeneous velocity modeling method based on seismic data Hessian optimization
By using the quasi-Newtonian method to approximate the solution of the higher-order derivative terms of the Hessian matrix in seismic data inversion, the problem of slow convergence of gradient-class methods and large calculations in the existing technology is solved, and high-precision and fast speed inversion is achieved.
Patent Information
- Application Number
- CN202311088918.1
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2023-08-28
- Publication Date
- 2025-02-28
AI Technical Summary
In the prior art, the gradient-class method is easy to implement but the convergence speed is slow and the inversion accuracy is not high; the Newtonian method is fast convergence speed and the inversion accuracy is high, but the calculation amount is large and the memory consumption is serious, which limits its application.
Based on the precise solution of the first derivative terms of the Hessian matrix, the quasi-Newtonian method is used to approximately solve the higher derivative terms of the Hessian matrix, reducing the calculation time and memory consumption.
The inversion accuracy and iterative convergence speed are improved, the calculation time and memory consumption are reduced, and high-precision and fast Newtonian earthquake velocity inversion is achieved.
Smart Images
Figure CN119535566A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of geophysical exploration velocity modeling, and more specifically to a heterogeneous velocity modeling method based on Hessian optimization of seismic data. Background Art
[0002] Nonlinear inversion of seismic data is a velocity modeling technology that can obtain underground structures and medium properties with high accuracy. This technology is achieved by fitting actual seismic records and simulated seismic records. First, an initial velocity is given, the error is calculated, and the velocity value is iteratively corrected until it converges to minimize the optimization objective function. The corresponding velocity at this time is the optimal velocity.
[0003] In the prior art, seismic waveform inversion is generally implemented based on gradient methods or Newton methods. The gradient method is easy to implement, but has a slow convergence speed and low inversion accuracy; the Newton method has a fast convergence speed and high inversion accuracy, but the huge amount of calculation and computer memory consumption limit the application of the Newton method. Summary of the invention
[0004] In order to overcome the defects existing in the above-mentioned prior art, the present invention discloses a heterogeneous velocity modeling method based on Hessian optimization of seismic data. The purpose of the present invention is to solve the problem that the gradient method in the prior art is easy to implement, but has a slow convergence speed and low inversion accuracy; the Newton method has a fast convergence speed and high inversion accuracy, but the huge amount of calculation and computer memory consumption limit the application of the Newton method. The present invention uses the quasi-Newton method to approximate the solution of the higher-order derivative terms of the Hessian matrix on the basis of accurately solving the first-order derivative terms of the Hessian matrix. The Hessian matrix reflects an implicit functional relationship between seismic wave field data and propagation medium attribute parameters. Compared with the gradient method, it can greatly improve the inversion accuracy and iterative convergence speed; compared with the conventional Newton method, it reduces the calculation time and computer memory consumption.
[0005] In order to achieve the above objectives, the technical solution adopted by the present invention is: A heterogeneous velocity modeling method based on Hessian optimization of seismic data includes the following steps: 1. Input observation data and initial model S1. Obtaining seismic observation data; establishing an initial velocity model for inversion and obtaining seismic simulation data based on the initial velocity model; In the present invention, the initial velocity model in the numerical simulation is independent of the seismic observation data, the two are in parallel and independent of each other, and the initial velocity model is obtained by linear smoothing of the real model.
[0006] In the above steps, observed seismic data are obtained, and an initial velocity model for inversion is established, which is used as the original input information of the inversion process.
[0007] Preferably, in the step S1, a real growth fault model and a staggered grid finite difference method are used to solve the acoustic wave equation to obtain the seismic observation data.
[0008] Preferably, in the step S1, the initial velocity model is obtained through a growth fault model using linear smoothing, and the acoustic wave equation is solved using the initial velocity model and a staggered grid finite difference method to obtain the seismic simulation data.
[0009] Preferably, in the step S1, when the finite difference method is used to solve the acoustic wave equation to obtain the seismic observation data and the seismic simulation data, the parameters are as follows: the calculation accuracy is 2nd order in time and 8th order in space, the non-split perfectly matched layer absorption model boundary is used, the source uses the Ricker wavelet with a main frequency of 40 Hz, the time interval is 0.5 ms, the lateral space interval is 5 m, the longitudinal space interval is 3 m, and the number of model grid points is , the recording time is 1.5s.
[0010] In the above steps, the numerical simulation in seismic nonlinear inversion is implemented using acoustic wave equation staggered grids and non-split perfectly matched layer absorbing boundaries (NPML), which has high accuracy and good absorption effect.
[0011] 2. Obtaining frequency division data based on convolution method S2, a seismic wavelet of a required main frequency is given, and the seismic wavelet is convolved with the seismic observation data and the seismic simulation data to obtain seismic frequency division data; In the above steps, given the seismic wavelet of the required main frequency, the seismic wavelet is convolved with the observed data and the simulated data to obtain the frequency division data within the seismic wavelet frequency band, as follows: in, represents the frequency division data after convolution processing, represents earthquake observation data or simulation data, represents convolution, Represents seismic wavelet.
[0012] In the above steps, different from the general method of using filters to achieve multi-scale, convolution is used to obtain frequency division data, which is simple to implement and has high computational efficiency.
[0013] 3. Obtaining seismic wave field residuals and establishing objective functions S3, using the seismic wave field residuals of the seismic simulation data and the observation data after frequency division to establish the objective function; Preferably, in the step S3, the seismic wave field residuals of the seismic simulation data and the observation data after frequency division are used, and an objective function is established based on the L2 norm, and the objective function is: in, is the objective function, is the velocity parameter in model space, To simulate earthquake data, To observe earthquake data.
[0014] 4. Calculate the high-order derivatives in the Hessian matrix based on the quasi-Newton method and obtain the iterative update direction; S4, using the seismic wave field residual and based on the quasi-Newton method to optimize the high-order derivative terms in the Hessian matrix, and then using the optimized high-order derivative terms to obtain the iterative update direction of the velocity model; Preferably, in the step S4, the iterative update direction is: in, Indicates the velocity model update direction, represents the residual of seismic data, represents the first-order Fréchet derivative of the simulated seismic wave field, represents the high-order derivative term, and k represents the number of iterations.
[0015] Preferably, in the iterative update direction of step S4, the first-order Fréchet derivative is: in, represents the change in the seismic wave field, represents the change in velocity model, Indicates time t The seismic wave field value at Indicates that at time 0, The seismic wave function generated by the earthquake source is located at time t The earthquake response, represents the velocity in model space, Represents convolution.
[0016] Preferably, in the iterative update direction of step S4, the high-order derivative term is: in, is the k+1th Hessian matrix high-order derivative, is the high-order term of the k-th Hessian matrix, k is the number of iterations, T is the symbol of the matrix transpose, , represents the velocity model difference between adjacent iterations, , where J represents the Jacobian matrix and v is the velocity in the model space. k =0, M0 is the unit matrix.
[0017] The specific inference of the iterative update direction of the above velocity model is as follows: In the present invention, the update direction of the quasi-Newton optimization method is defined as: in, Indicates the velocity model update direction, represents the residual of seismic data, represents the Fréchet derivative of the simulated seismic wave field, represents the product of the wavefield residual and the second-order Fréchet derivative.
[0018] The Hessian matrix in the quasi-Newton inversion method is defined as: Calculating the Hessian matrix consumes a lot of computing time and machine memory, and is difficult to implement directly. The main reason is that the amount of calculation and memory consumption of estimating the product of the seismic wave field residual and the second-order Fréchet derivative are very large. The present invention first accurately calculates the first-order Fréchet derivative of the Hessian matrix, and then uses the quasi-Newton method to approximate the product of the seismic wave field residual and the second-order Fréchet derivative. Therefore, memory consumption and computing time can be reduced.
[0019] First, the Green's function method of the acoustic wave equation is used to obtain the Fréchet derivative of the wave field with respect to the model parameters. The Green's function of the two-dimensional acoustic wave equation should satisfy the following equation: in, represents the Green's function, which is specifically expressed as , which is a function of spatial position and time; represents an impulse function.
[0020] Secondly, the Fréchet derivatives of the seismic wave field and velocity model parameters are calculated using the following formula: in, represents the change in the seismic wave field, represents the change in velocity model, Indicates time t The seismic wave field value at Indicates that at time 0, The seismic wave function generated by the earthquake source is located at time t The earthquake response, represents the velocity in model space, Represents convolution.
[0021] Subsequently, the quasi-Newton method is used to approximate the high-order derivatives of the Hessian matrix. In the k+1th iteration, the matrix Approximate the higher-order derivatives of the Hessian matrix and multiply by The following approximate equation can be obtained: The above equation can be approximated using the first-order Taylor expansion: make ,make , where J represents the Jacobian matrix.
[0022] The following quasi-Newton conditions need to be satisfied: Finally, using quasi-Newton conditions and related theories, we get the approximate calculation matrix The update equation is: Among them, the initial is a pre-given positive definite matrix.
[0023] After using the quasi-Newton method to approximate the high-order terms of the Hessian matrix, the update direction equation of the present invention can be obtained: In the above steps, unlike the general direct calculation of the Hessian matrix, the quasi-Newton method is used to approximate the high-order derivatives in the Hessian matrix, which can accurately solve the first-order Fréchet derivatives; the main characteristics of the Hessian matrix are retained, so the velocity inversion accuracy and convergence speed can be improved.
[0024] 5. Update step size and update model based on Barzilar-Borwein method and non-monotonicity criterion S5, using the objective function and based on the Barzilar-Borwein method and the non-monotonicity criterion to obtain the iteration step of the velocity model, and then using the iteration step to update the velocity model; Preferably, in the S5 step, in the first iteration, a constant is first given as the trial step length. If the trial step length satisfies the non-monotonicity criterion, it is used as the update step length for this time, otherwise the trial step length is changed until the non-monotonicity criterion is met; when it is not the first iteration, the step length required by the linear method is used as the trial step length. If the trial step length satisfies the non-monotonicity criterion, it is accepted, otherwise the trial step length is changed until the non-monotonicity criterion is met.
[0025] Preferably, in step S5, the non-monotonic criterion is: in, is the objective function value, For speed, is the step length, For direction, is a constant, is a real number sequence, Satisfy the recursive relationship ,in , is a constant, Satisfy the recursive relationship ,in ,parameter .
[0026] The above objective function is the objective function in step S3.
[0027] Preferably, in step S5, the Barzilar-Borwein method is used to obtain the trial step length, which is as follows: in, is the step size, k is the number of iterations, T is the sign of the matrix transpose, represents the difference between adjacent iteration speed models, represents the difference in update directions between adjacent iterations, For speed, direction; when the number of iterations is an odd number, BB1 in the Barzilar-Borwein method is used to calculate the trial step size; when the number of iterations is an even number, BB2 in the Barzilar-Borwein method is used to calculate the trial step size.
[0028] Preferably, in step S5, the updated velocity model is: in, For speed, is the step length, is the direction, and k is the number of iterations.
[0029] In the above steps, combining the Barzilar-Borwein method with the non-monotone linear search criterion can effectively reduce the amount of linear search calculation and speed up the inversion.
[0030] 6. Loop Iteration S6, judging whether the updated velocity model satisfies the termination condition or reaches the maximum number of iterations, if so, terminating the iterative inversion, otherwise entering step S3 and entering the next loop iteration.
[0031] Preferably, the step S6 comprises the following steps: S61, determining whether the updated velocity model meets the termination condition or reaches the maximum number of iterations, if yes, proceeding to step S62, if no, returning to step S3 to enter the next loop iteration; S62, determine whether the maximum frequency is reached, if so, output the result, if not, update the frequency band, use the low-frequency inversion result as the high-frequency inversion initial model, and return to step S2 for loop iteration.
[0032] In the above steps, the maximum frequency represents the main frequency of the observed seismic data.
[0033] Beneficial effects of the present invention: 1. Based on the accurate solution of the first-order derivative of the Hessian matrix, the present invention uses the quasi-Newton method to approximately solve the higher-order derivative of the Hessian matrix. The Hessian matrix reflects an implicit functional relationship between the seismic wave field data and the propagation medium attribute parameters. Compared with the gradient method, it can greatly improve the inversion accuracy and iterative convergence speed; compared with the conventional Newton method, it reduces the calculation time and computer memory consumption.
[0034] 2. Calculating the Hessian matrix consumes a lot of computing time and machine memory, and is difficult to implement directly. The main reason is that the amount of calculation and memory consumption of estimating the product of the seismic wave field residual and the second-order Fréchet derivative is very large. The present invention first accurately calculates the first-order Fréchet derivative of the Hessian matrix, and then uses the quasi-Newton method to approximately calculate the product of the seismic wave field residual and the second-order Fréchet derivative. While ensuring the calculation accuracy of the Hessian matrix, the present invention reduces the amount of calculation and computer memory consumption of solving the high-order derivative terms in the Hessian matrix. It is a high-precision, fast Newton-type seismic velocity inversion optimization method, which accelerates the convergence speed during the calculation process, and can improve the velocity inversion accuracy, and can accurately describe the underground structure and medium properties.
[0035] 3. The present invention combines the Barzilar-Borwein method with the non-monotonic linear search criterion, which can effectively reduce the amount of linear search calculation and speed up the inversion speed.
[0036] 4. The present invention is different from the general method of realizing multi-scale by using filters. Convolution is used to obtain frequency division data, which is simple to realize and has high computational efficiency.
[0037] 5. The numerical simulation in the seismic nonlinear inversion of the present invention is realized by using the acoustic wave equation staggered grid and non-split perfectly matched layer absorbing boundary (NPML), which has high accuracy and good absorption effect. BRIEF DESCRIPTION OF THE DRAWINGS
[0038] Figure 1 This is the velocity inversion algorithm flow of the example; Figure 2 is the true velocity model of the instance; Figure 3 is the seismic single shot record of the real velocity model of the instance; Figure 4 is the initial velocity model of the instance; Figure 5 is the seismic single shot record of the initial velocity model of the example; Figure 6 is the approximate Hessian matrix structure of the instance; Figure 7 is the step update search process of the instance; Figure 8 is the final inverted velocity model of the example. DETAILED DESCRIPTION
[0039] The concept, specific structure and technical effects of the present invention will be clearly and completely described below in conjunction with the embodiments and drawings to fully understand the purpose, characteristics and effects of the present invention.
[0040] Example 1 A heterogeneous velocity modeling method based on Hessian optimization of seismic data, such as Figure 1 As shown, the following steps are included: S1. Obtaining seismic observation data; establishing an initial velocity model for inversion and obtaining seismic simulation data based on the initial velocity model; S2, a seismic wavelet of a required main frequency is given, and the seismic wavelet is convolved with the seismic observation data and the seismic simulation data to obtain seismic frequency division data; S3, using the seismic wave field residuals of the seismic simulation data and the observation data after frequency division to establish the objective function; S4, using the seismic wave field residual and based on the quasi-Newton method to optimize the high-order derivative terms in the Hessian matrix, and then using the optimized high-order derivative terms to obtain the iterative update direction of the velocity model; S5, using the objective function and based on the Barzilar-Borwein method and the non-monotonicity criterion to obtain the iteration step of the velocity model, and then using the iteration step to update the velocity model; S6, judging whether the updated velocity model satisfies the termination condition or reaches the maximum number of iterations, if so, terminating the iterative inversion, otherwise entering step S3 and entering the next loop iteration.
[0041] Example 2 This embodiment further elaborates on step S1 based on embodiment 1. In step S1, the input seismic observation data and initial velocity model will include: using a real growth fault model (such as Figure 2 ) and staggered grid finite difference method to solve the acoustic wave equation and obtain seismic observation data (such as Figure 3 ); linear smoothing is used to obtain the initial velocity model through the growth fault model (as shown in Figure 4 When solving the acoustic wave equation with the finite difference method, the following parameters are used: the calculation accuracy is 2nd order in time and 8th order in space, the non-split perfectly matched layer (NPML) is used to absorb the model boundary, the source uses the Ricker wavelet with a main frequency of 40 Hz, the time interval is 0.5 ms, the horizontal space interval is 5 m, the vertical space interval is 3 m, and the number of model grid points is , the recording time is 1.5s.
[0042] In step S1, the initial velocity model and staggered grid finite difference method are used to solve the acoustic wave equation to obtain seismic simulation data (such as Figure 5 as shown).
[0043] Example 3 This embodiment further elaborates on step S2 based on embodiment 2. In step S2, seismic frequency data is obtained based on the convolution method, and the observation data and the simulation data are convolved with the seismic wavelet respectively, so as to obtain seismic frequency data. In this example, Ricker wavelets with main frequencies of 2 Hz, 5 Hz, 10 Hz, 20 Hz and 40 Hz are respectively used as convolution wavelets, and then the low-frequency inversion result is used as the initial velocity model of high-frequency inversion.
[0044] Example 4 This embodiment further elaborates on step S3 based on embodiment 3. In step S3, obtaining the seismic wave field residual and establishing the objective function will include: subtracting the seismic observation record after frequency selection and the simulation data to obtain the seismic wave field residual; and then using the wave field residual to establish the objective function based on the L2 norm. If multi-shot data is used as an example, the following formula can be used for cumulative summation.
[0045] in, is the objective function, is the velocity parameter in model space, To simulate earthquake data, To observe earthquake data.
[0046] Example 5 This embodiment further elaborates on step S4 based on embodiment 4. In step S4, the high-order derivative items in the Hessian matrix are calculated based on the quasi-Newton method, and the update direction of the cyclic iteration is obtained, which includes: using the Green function method to obtain the Fréchet derivative of the seismic wave field; using the quasi-Newton method to approximately calculate the high-order derivative items in the Hessian matrix; and obtaining the iterative update direction.
[0047] The Hessian matrix can be decomposed into a combination of first-order Fréchet derivatives and second-order Fréchet derivatives. This example accurately calculates the first-order Fréchet derivative of the Hessian matrix and approximates the higher-order derivatives of the Hessian matrix using the quasi-Newton method.
[0048] Where H represents the Hessian matrix.
[0049] The first-order Fréchet derivative of the seismic wave field with respect to the model parameters can be accurately calculated using the following formula: in, represents the change in the seismic wave field, represents the change in velocity model, Indicates time t The seismic wave field value at Indicates that at time 0, The seismic wave function generated by the earthquake source is located at time t The earthquake response, represents the velocity in model space, Represents convolution.
[0050] The high-order derivatives in the Hessian matrix can be approximated using the following formula: in, is the k+1th Hessian matrix high-order derivative, is the high-order term of the k-th Hessian matrix, k is the number of iterations, T is the symbol of the matrix transpose, , represents the velocity model difference between adjacent iterations, , where J represents the Jacobian matrix and v is the velocity in the model space.k =0, M0 is the unit matrix.
[0051] The initial M0 is a pre-given positive definite matrix. The present invention uses the unit matrix as the initial positive definite matrix of the quasi-Newton method. The approximate Hessian matrix (e.g. Figure 6 as shown).
[0052] The update direction in the example of the present invention is defined as follows: Example 6 This embodiment further elaborates on step S5 on the basis of embodiment 5. In step S5, the update step size and the update model are obtained based on the Barzilar-Borwein method and the non-monotonicity criterion, which will include: in the first iteration, a constant is first given as the trial step size. If the non-monotonicity criterion is met, it is used as the update step size for this time, otherwise the trial step size is changed until the non-monotonicity criterion is met; when it is not the first iteration, the step size obtained by the linear method is used as the trial step size, and it is accepted if the non-monotonicity criterion is met, otherwise the trial step size is changed until the non-monotonicity criterion is met. In particular, when the number of iterations is an odd number, the BB1 formula in the Barzilar-Borwein method is used for calculation, and when it is an even number, the BB2 formula in the Barzilar-Borwein method is used for calculation. The specific calculation process is as follows: Figure 7 shown.
[0053] The iteration step size for calculating the inversion using the non-monotonicity criterion and the Barzilar-Borwein method is: in, is the objective function value, For speed, is the step length, For direction, is a constant, is a real number sequence, Satisfy the recursive relationship ,in , is a constant, Satisfy the recursive relationship ,in ,parameter .
[0054] If the trial step size meets the non-monotonicity criterion, the trial step size is used as the iteration step size; otherwise, the trial step size is changed and recalculated until the non-monotonicity criterion is met.
[0055] The Barzilar-Borwein method calculates the trial step size using the following formula: in, is the step size, k is the number of iterations, T is the sign of the matrix transpose, represents the difference between adjacent iteration speed models, represents the difference in update directions between adjacent iterations, For speed, direction; when the number of iterations is an odd number, BB1 in the Barzilar-Borwein method is used to calculate the trial step size; when the number of iterations is an even number, BB2 in the Barzilar-Borwein method is used to calculate the trial step size.
[0056] Finally, the velocity model is updated. When the error of the inversion result meets the termination condition or reaches the maximum number of iterations, the inversion result is output. Figure 8 The method of the present invention can obtain good inversion results, which are basically consistent with the real model.
[0057] in, For speed, is the step length, is the direction, and k is the number of iterations.
[0058] The above is a specific description of the implementation mode of the present invention, but the present invention is not limited to the described embodiments. Those skilled in the art may make various equivalent modifications or substitutions without violating the spirit of the present invention, and these equivalents or substitutions are all included in the scope defined by the claims of the present invention.
Claims
1. A heterogeneous velocity modeling method based on Hessian optimization of seismic data, characterized in that: The following steps are involved: S1. Obtaining seismic observation data; establishing an initial velocity model for inversion and obtaining seismic simulation data based on the initial velocity model; S2, a seismic wavelet of a required main frequency is given, and the seismic wavelet is convolved with the seismic observation data and the seismic simulation data to obtain seismic frequency division data; S3, using the seismic wave field residuals of the seismic simulation data and the observation data after frequency division to establish the objective function; S4, using the seismic wave field residual and based on the quasi-Newton method to optimize the high-order derivative terms in the Hessian matrix, and then using the optimized high-order derivative terms to obtain the iterative update direction of the velocity model; S5, using the objective function and based on the Barzilar-Borwein method and the non-monotonicity criterion to obtain the iteration step of the velocity model, and then using the iteration step to update the velocity model; S6, judging whether the updated velocity model satisfies the termination condition or reaches the maximum number of iterations, if so, terminating the iterative inversion, otherwise entering step S3 and entering the next loop iteration.
2. The heterogeneous velocity modeling method according to claim 1, characterized in that: In the step S1, the acoustic wave equation is solved using a real growth fault model and a staggered grid finite difference method to obtain the seismic observation data; the initial velocity model is obtained through the growth fault model using linear smoothing, and the acoustic wave equation is solved using the initial velocity model and a staggered grid finite difference method to obtain the seismic simulation data.
3. The heterogeneous velocity modeling method according to claim 1, characterized in that: In the step S3, the seismic wave field residuals of the seismic simulation data and the observation data after frequency division are used, and an objective function is established based on the L2 norm. The objective function is: ， in, is the objective function, is the velocity parameter in model space, To simulate earthquake data, To observe earthquake data.
4. The heterogeneous velocity modeling method according to claim 1, characterized in that: In the step S4, the iterative update direction is: ， in, Indicates the velocity model update direction, represents the residual of seismic data, represents the first-order Fréchet derivative of the simulated seismic wave field, represents the high-order derivative term, and k represents the number of iterations.
5. The heterogeneous velocity modeling method according to claim 4, characterized in that: In the iterative update direction of step S4, the first-order Fréchet derivative is: ， in, represents the change in the seismic wave field, represents the change in velocity model, Indicates time t The seismic wave field value at Indicates that at time 0, The seismic wave function generated by the earthquake source is located at time t The earthquake response, represents the velocity in model space, Represents convolution.
6. The heterogeneous velocity modeling method according to claim 5, characterized in that: In the iterative update direction of step S4, the high-order derivative term is: ， in, is the k+1th Hessian matrix high-order derivative, is the high-order term of the k-th Hessian matrix, k is the number of iterations, T is the symbol of the matrix transpose, , represents the velocity model difference between adjacent iterations, , where J represents the Jacobian matrix and v is the velocity in the model space. k =0, M0 is the unit matrix.
7. The heterogeneous velocity modeling method according to claim 1, characterized in that: In the step S5, during the first iteration, a constant is first given as a trial step size. If the trial step size satisfies the non-monotonicity criterion, it is used as the update step size for this time. Otherwise, the trial step size is changed until the non-monotonicity criterion is met. When it is not the first iteration, the step size obtained by the linear method is used as the trial step size. If the trial step size satisfies the non-monotonicity criterion, it is accepted. Otherwise, the trial step size is changed until the non-monotonicity criterion is met.
8. The heterogeneous velocity modeling method according to claim 7, characterized in that: In step S5, the non-monotonic criterion is: ， in, is the objective function value, For speed, is the step length, For direction, is a constant, is a real number sequence, Satisfy the recursive relationship ,in , is a constant, Satisfy the recursive relationship ,in ,parameter .
9. The heterogeneous velocity modeling method according to claim 8, characterized in that: In the step S5, the Barzilar-Borwein method is used to obtain the trial step length, which is as follows: ， in, is the step size, k is the number of iterations, T is the sign of the matrix transpose, represents the difference between adjacent iteration speed models, represents the difference in update directions between adjacent iterations, For speed, direction; when the number of iterations is an odd number, the trial step length is calculated using BB1 in the Barzilar-Borwein method; when the number of iterations is an even number, the trial step length is calculated using BB2 in the Barzilar-Borwein method; In step S5, the updated velocity model is: ， in, For speed, is the step length, is the direction, and k is the number of iterations.
10. The heterogeneous velocity modeling method according to claim 1, characterized in that: The step S6 comprises the following steps: S61, determining whether the updated velocity model meets the termination condition or reaches the maximum number of iterations, if yes, proceeding to step S62, if no, returning to step S3 to enter the next loop iteration; S62, determine whether the maximum frequency is reached, if so, output the result, if not, update the frequency band, use the low-frequency inversion result as the high-frequency inversion initial model, and return to step S2 for loop iteration.