An elastic wave inversion method based on natural frequency-time domain energy spectrum

Abstract

The application relates to an elastic wave inversion method based on natural frequency division time-frequency energy spectrum, which extracts single-frequency energy spectrum information of seismic data by using S transformation and realizes natural frequency division multi-scale inversion by embedding a frequency point cycle, establishes a more accurate initial model of P-wave and S-wave velocity under the condition of low-frequency loss, and on the basis, uses a multi-scale strategy to carry out frequency division group iterative updating of elastic wave full waveform inversion, so that an accurate P-wave and S-wave velocity model is obtained, and the result accuracy and reliability of the elastic wave full waveform inversion under the condition of low-frequency loss are improved.

Description

Technical Field

[0001] This invention relates to the field of seismic exploration inversion technology, specifically to an elastic wave inversion method based on natural frequency division time-frequency domain energy spectrum. Background Technology

[0002] Velocity modeling is one of the most prominent problems in seismic exploration. The accuracy of velocity modeling affects the realism and reliability of seismic imaging and interferes with the interpretation of seismic data. Establishing a high-precision velocity model is key to solving this problem. There are three main types of velocity modeling methods: the first is travel-time tomography, which only utilizes travel-time information from seismic records and can only obtain a low-precision velocity model; the second is velocity analysis, which is affected by the subjective factors of the operator; and the third is the full waveform inversion method, which utilizes all information in the seismic record and can obtain a high-precision velocity model.

[0003] Elastic wave full waveform inversion is more nonlinear than acoustic wave full waveform inversion. The difficulty of elastic wave full waveform inversion lies in the greater dependence of the shear wave velocity model inversion on the initial model and the huge amount of computation involved in elastic wave inversion. The former affects the accuracy of the shear wave velocity model in elastic wave full waveform inversion, while the latter is directly related to whether the relevant technology can be applied to production.

[0004] Sirgue proposed a hybrid-domain full-waveform inversion method based on time-domain forward modeling and frequency-domain inversion. This method combines the advantages of both methods by incorporating Fourier transform into the time-domain forward modeling process to obtain the frequency-domain wave field. Compared to frequency-domain full-waveform inversion, this method saves memory space and is easily implemented in large-scale parallel computation, making it suitable for practical production work. Wang Yuwei et al. proposed an envelope-based multi-scale inversion method. This method utilizes envelope inversion to provide an initial model containing large-scale information for hybrid-domain elastic wave full-waveform inversion, greatly alleviating the "cycle jump" phenomenon in hybrid-domain elastic wave full-waveform inversion.

[0005] The envelope-based multi-scale inversion method first obtains the seismic record envelope and uses it for inversion to obtain a large-scale initial model. Then, using the obtained initial model, Sirgue's hybrid domain full waveform inversion method is applied to establish a high-precision velocity model. Envelope inversion is based on the full frequency band of the seismic record. If frequency-band envelope inversion is to be performed, additional low-pass filtering processing is required on the seismic record. The setting of the low-pass filter in this process will also affect the effect of envelope inversion. Therefore, although this approach reduces the dependence of elastic wave full waveform inversion on low-frequency information to a certain extent, the additional low-pass filter increases the workload and makes the frequency-division framework of multi-scale inversion unnatural. Summary of the Invention

[0006] To address the problem that frequency division inversion in the envelope-based full waveform inversion technique for elastic waves requires the introduction of an additional low-pass filter in the above-mentioned existing technical solutions, this invention provides an elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum, which reduces the generation of errors, lowers the probability of inversion getting trapped in local extrema, and improves the accuracy of the velocity model.

[0007] The technical solution adopted in this invention is: an elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum, comprising the following steps:

[0008] S1: Input the measured three-component seismic record, read the field observation system from the track head of the measured data, read the seismic wavelet, inversion parameters and the initial P-wave velocity model mvp1(0) and S-wave velocity model mvs1(0), and perform grid subdivision of the computational space based on the inversion parameters;

[0009] S2: Perform forward modeling of the elastic wave equation using the set parameters and the current model to obtain the simulated seismic records currecordx and currecordz of the x and z components under the current velocity model;

[0010] S3: Read in the measured seismic records orirecordx and orirecordz for the x and z components, and determine the frequency range for time-frequency domain inversion. , Setting the inversion frequency interval Number of inversion frequency points The single-frequency S-transform of the simulated seismic records of the x and z components is used to obtain curx and curz, and the single-frequency S-transform of the measured seismic records of the x and z components is used to obtain orix and oriz. The squares of the amplitudes of curx, curz, orix, and oriz are taken respectively to obtain their respective single-frequency energy spectra.

[0011] S4: Using the time-frequency domain energy spectrum elastic wave full waveform inversion with natural frequency division, multi-scale inversion is performed to obtain the longitudinal and transverse wave layer velocity model containing large-scale information;

[0012] S5: Using the velocity model containing large-scale information as the initial model, perform hybrid domain multi-scale elastic wave full waveform inversion to obtain the final high-precision inversion result.

[0013] In this technical solution, the seismic record is transformed to the time-frequency domain by S-transform and inversion is performed using the energy spectrum of different single frequency points. This eliminates the step of designing a low-pass filter in conventional envelope inversion, reducing workload and difficulty. At the same time, the hybrid domain multi-scale full waveform inversion is used to gradually invert from the low-frequency band to the high-frequency band, reducing the probability of the inversion getting trapped in local extrema and improving the accuracy of the velocity model.

[0014] Preferably, step S4 utilizes the time-frequency domain energy spectrum elastic wave full waveform inversion based on natural frequency division to perform multi-scale inversion and obtain a longitudinal and transverse wave layer velocity model containing large-scale information. The specific steps are as follows:

[0015] Assuming the current frequency point is f, and the maximum number of iterations at the current frequency point is Niter1, the initial longitudinal wave model mvp1(j,f) and the initial transverse wave model mvs1(j,f) for the j-th iteration at the f-th frequency point are obtained according to S3, and the single-frequency energy spectra of curx, curz, orix, and oriz are obtained respectively. , , , The associated sources adjx and adjz of the full waveform inversion are calculated as follows:

[0016]

[0017] ;

[0018] in, In order to seek the truth, Indicates the recording time. Represents the imaginary unit. Indicates a time delay;

[0019] The adjoint wave field is obtained by backpropagating adjx and adjz using the elastic wave equation. The gradient is obtained by using the adjoint wave field and the simulated wave field according to the gradient class or Newton class method. The longitudinal wave update direction directionvp1 and the transverse wave update direction directionvs1 are determined for this iteration, and the update step size a is determined.

[0020] The longitudinal wave model update formula is mvp1(j+1,f)=mvp1(j,f)+a*directionvp1;

[0021] The update formula for the transverse wave model is mvs1(j+1,f)=mvs1(j,f)+a*directionvs1;

[0022] Update the velocity model, and determine if the current iteration number j of frequency point f is equal to Niter1. If it is, proceed to the next frequency point; otherwise, proceed to the next iteration of the current frequency point. If the values ​​are equal, output the final P-wave velocity model and S-wave velocity model, then f = f + 1, j = 0, and proceed to the iteration of the next frequency point.

[0023] Preferably, step S5 includes using a velocity model containing large-scale information as an initial model to perform hybrid domain multi-scale elastic wave full waveform inversion to obtain a high-precision velocity model. The specific steps are as follows:

[0024] S51: Set the range of each frequency band [ifmin, ifmax], the number of frequency bands groupnum, the number of frequencies in the frequency band fnum, the current frequency band number k, the initial model minit containing large-scale information read in, the source wavelet wavelet, set the horizontal grid spacing dx, the vertical grid spacing dz, the time sampling interval dt, the total time sampling points Tn, the number of horizontal grid points Xn, the number of vertical grid points Zn, the maximum number of iterations Nitermax, the current number of iterations j, and assume that the P-wave model and S-wave model of the k-th frequency band in the iteration process are mvp(k,j) and mvs(k,j) respectively;

[0025] S52: Use wavelet, minit and set parameters to perform forward modeling of the elastic wave equation to obtain the simulated records CURX and CURZ of the x and z components, and record the x and z components ux and uz of the wave field value at each moment, and read in the measured records orix and oriz;

[0026] S53: Perform discrete Fourier transform on ux and uz according to the current frequency band range [ifmin ifmax] to obtain fux and fuz. Subtract the simulated record from the measured record to obtain the record residuals Residualx and Residualz.

[0027] S54: Perform discrete Fourier transform on ux and uz according to the current frequency band range [ifmin ifmax] to obtain fux and fuz, perform discrete Fourier transform on CURX and orix and calculate the difference to obtain Residualx, and perform discrete Fourier transform on CURZ and oriz and calculate the difference to obtain Residualz.

[0028] S55: Calculate the gradients gradvp and gradvs using the cross-correlation of fux, fuz, frux, and fruz;

[0029] S56: Gradient preprocessing determines the update directions directionvp and directionvs, and calculates the update step size a;

[0030] S57: Update the model: mvp(k,j+1)=mvp(k,j)+a*directionvp;

[0031] mvs(k,j+1)=mvs(k,j)+a*directionvs;

[0032] S58: Determine if j is equal to Nitermax; <1> If so, proceed to step S59; <2> Otherwise, let j = j + 1 and proceed to step S52;

[0033] S59: Determine if k is equal to groupnum; <1> If so, output the final P-wave velocity model and S-wave velocity model; <2> Otherwise, let k=k+1, j=0, and proceed to step S52.

[0034] Compared with existing technologies, the advantages are as follows: This invention proposes an initial model construction method based on time-frequency domain energy spectrum elastic wave inversion. By transforming the seismic record to the time-frequency domain through S-transform and taking the energy spectrum of the lowest usable frequency for inversion, the step of designing a low-pass filter in conventional envelope inversion is eliminated, reducing the generation of errors. At the same time, by using hybrid domain multi-scale full waveform inversion to gradually invert from the low-frequency band to the high-frequency band, the probability of inversion getting trapped in local extrema is reduced, and the accuracy of the velocity model is improved. Attached Figure Description

[0035] Figure 1 This is a flowchart of an elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum according to the present invention;

[0036] Figure 2 This is a detailed flowchart of the elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum of the present invention;

[0037] Figure 3 This is a precise velocity model diagram of an elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum according to the present invention;

[0038] Figure 4 This is the initial linear model diagram of an elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum of the present invention;

[0039] Figure 5 This is a traditional hybrid domain multi-scale full waveform inversion result diagram of an elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum according to the present invention;

[0040] Figure 6 This is a time-frequency domain energy spectrum inversion result diagram of an elastic wave inversion method based on natural frequency division time-frequency domain energy spectrum according to the present invention;

[0041] Figure 7 The energy spectrum basis of the elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum of this invention is to further perform hybrid domain multi-scale inversion effect diagram. Detailed Implementation

[0042] The accompanying drawings are for illustrative purposes only and should not be construed as limiting this patent. To better illustrate this embodiment, some components in the drawings may be omitted, enlarged, or reduced, and do not represent the actual product dimensions. It is understandable to those skilled in the art that some well-known structures and their descriptions may be omitted in the drawings. The positional relationships described in the drawings are for illustrative purposes only and should not be construed as limiting this patent.

[0043] In the accompanying drawings of the embodiments of the present invention, the same or similar reference numerals correspond to the same or similar parts. In the description of the present invention, it should be understood that if terms such as "upper," "lower," "left," "right," "long," and "short" indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, they are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, the terms used to describe positional relationships in the drawings are only for illustrative purposes and should not be construed as limiting the present patent. For those skilled in the art, the specific meaning of the above terms can be understood according to the specific circumstances.

[0044] The technical solution of the present invention will be further described in detail below through specific embodiments and in conjunction with the accompanying drawings:

[0045] Example 1

[0046] like Figures 1-2 The following is an embodiment 1 of an elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum, which includes the following steps:

[0047] S1: Input the measured three-component seismic record, read the field observation system from the track head of the measured data, read the seismic wavelet, inversion parameters and the initial P-wave velocity model mvp1(0) and S-wave velocity model mvs1(0), and perform grid subdivision of the computational space based on the inversion parameters;

[0048] S2: Perform forward modeling of the elastic wave equation using the set parameters and the current model to obtain the simulated seismic records currecordx and currecordz of the x and z components under the current velocity model;

[0049] S3: Read in the measured seismic records orirecordx and orirecordz for the x and z components, and determine the frequency range for time-frequency domain inversion. , Setting the inversion frequency interval Number of inversion frequency points The single-frequency S-transform of the simulated records of the x and z components is used to obtain curx and curz, and the single-frequency S-transform of the measured seismic records of the x and z components is used to obtain orix and oriz. The squares of the amplitudes of curx, curz, orix, and oriz are taken respectively to obtain their respective single-frequency energy spectra.

[0050] S4: Using the time-frequency domain energy spectrum elastic wave full waveform inversion with natural frequency division, multi-scale inversion is performed to obtain the longitudinal and transverse wave layer velocity model containing large-scale information;

[0051] S5: Using the velocity model containing large-scale information as the initial model, perform hybrid domain multi-scale elastic wave full waveform inversion to obtain the final high-precision inversion result.

[0052] The beneficial effects of this embodiment are as follows: In this technical solution, the seismic record is transformed to the time-frequency domain by S-transform and the energy spectrum of different single frequency points is used for inversion. This eliminates the step of designing a low-pass filter in conventional envelope inversion, reducing the workload and the difficulty of the work. At the same time, by using hybrid domain multi-scale full waveform inversion to gradually invert from the low frequency band to the high frequency band, the probability of the inversion getting trapped in local extrema is reduced, and the accuracy of the velocity model is improved.

[0053] Example 2

[0054] Example 2 of an elastic wave inversion method based on natural frequency division time-frequency domain energy spectrum further defines the steps in Example 1 based on Example 1.

[0055] Specifically, step S4 utilizes the time-frequency domain energy spectrum elastic wave full waveform inversion based on natural frequency division to perform multi-scale inversion and obtain the P-wave and S-wave layer velocity model containing large-scale information. The specific steps are as follows:

[0056] Assuming the current frequency point is f, and the maximum number of iterations at the current frequency point is Niter1, the initial longitudinal wave model mvp1(j,f) and the initial transverse wave model mvs1(j,f) for the j-th iteration at the f-th frequency point are obtained according to S3, and the single-frequency energy spectra of curx, curz, orix, and oriz are obtained respectively. , , , The associated sources adjx and adjz of the full waveform inversion are calculated as follows:

[0057]

[0058] ;

[0059] The adjoint wave field is obtained by backpropagating adjx and adjz using the elastic wave equation. The gradient is then calculated using the adjoint wave field and the simulated wave field according to gradient-based or Newton-based methods. This determines the P-wave update direction (directionvp1) and the S-wave update direction (directionvs1) for this iteration, and finally, the update step size. a ;

[0060] Longitudinal wave model update formula: mvp1(j+1,f)=mvp1(j,f)+ a *directionvp1;

[0061] The shear wave model update formula is mvs1(j+1,f)=mvs1(j,f)+ a*directionvs1;

[0062] Update the velocity model, and determine if the current iteration number j of frequency point f is equal to Niter1. If it is, proceed to the next frequency point; otherwise, proceed to the next iteration of the current frequency point. If the values ​​are equal, output the final P-wave velocity model and S-wave velocity model, then f = f + 1, j = 0, and proceed to the iteration of the next frequency point.

[0063] Specifically, step S5 includes using the velocity model containing large-scale information as the initial model to perform hybrid domain multi-scale elastic wave full waveform inversion to obtain a high-precision velocity model. The specific steps are as follows:

[0064] S51: Set the range of each frequency band [ifmin ifmax], the number of frequency bands groupnum, the number of frequencies in the frequency band fnum, the current frequency band number k, the initial model minit containing large-scale information read in, the source wavelet wavelet, set the horizontal grid spacing dx, the vertical grid spacing dz, the time sampling interval dt, the total time sampling points Tn, the number of horizontal grid points Xn, the number of vertical grid points Zn, the maximum number of iterations Nitermax, the current number of iterations j, and assume that the P-wave model and S-wave model of the k-th frequency band in the model iteration process are mvp(k,j) and mvs(k,j) respectively;

[0065] S52: Use wavelet, minit and set parameters to perform forward modeling of the elastic wave equation to obtain the simulated records CURX and CURZ of the x and z components, and record the x and z components ux and uz of the wave field value at each moment, and read in the measured records orix and oriz;

[0066] S53: Perform discrete Fourier transform on ux and uz according to the current frequency band range [ifmin ifmax] to obtain fux and fuz. Subtract the simulated record from the measured record to obtain the record residuals Residualx and Residualz.

[0067] S54: Perform discrete Fourier transform on ux and uz according to the current frequency band range [ifmin ifmax] to obtain fux and fuz, perform discrete Fourier transform on CURX and orix and calculate the difference to obtain Residualx, and perform discrete Fourier transform on CURZ and oriz and calculate the difference to obtain Residualz.

[0068] S55: Calculate the gradients gradvp and gradvs using the cross-correlation of fux, fuz, frux, and fruz;

[0069] S56: Gradient preprocessing determines the update directions directionvp and directionvs, and calculates the update step size a;

[0070] S57: Update the model: mvp(k,j+1)=mvp(k,j)+a*directionvp;

[0071] mvs(k,j+1)=mvs(k,j)+a*directionvs;

[0072] S58: Determine if j is equal to Nitermax;

[0073] If yes, proceed to step S59; if no, let j = j + 1 and proceed to step S52.

[0074] S59: Determine if k is equal to groupnum;

[0075] If yes, output the final P-wave velocity model and S-wave velocity model; if no, let k=k+1, j=0, and proceed to step S52.

[0076] The beneficial effects of this embodiment are as follows: This invention proposes an initial model construction method based on time-frequency domain energy spectrum elastic wave inversion. By transforming the seismic record to the time-frequency domain through S-transform and taking the energy spectrum of the lowest usable frequency for inversion, the step of designing a low-pass filter in conventional envelope inversion is eliminated, reducing the generation of errors. At the same time, by using hybrid domain multi-scale full waveform inversion to gradually invert from the low-frequency band to the high-frequency band, the probability of inversion getting trapped in local extrema is reduced, and the accuracy of the velocity model is improved.

[0077] Example 3

[0078] Example 3 of an elastic wave inversion method based on natural frequency division time-frequency domain energy spectrum, as follows: Figures 3-7 As shown, based on Example 1 or Example 2, the advantages of this method are illustrated by example.

[0079] Specifically, such as Figure 3 As shown, the accurate model used employs Ricker wavelet excitation at a dominant frequency of 8 Hz, with both the transverse and longitudinal grid step sizes being 10 m, a time sampling interval of 1 ms, a total sampling time of 3 s, a shot spacing of 100 m, and full-array reception. The shot and receiver points are located on the Earth's surface, resulting in seismic records. Low-frequency information below 2.5 Hz is filtered out from the seismic records to obtain the measured records lacking low frequencies. Note that the ratio of P-wave velocity to S-wave velocity is fixed at 1.6.

[0080] like Figure 4As shown, the 1D linear initial model used is illustrated. It can be seen that the model does not contain any construction information, which is a significant challenge for the full waveform inversion of elastic waves.

[0081] like Figure 5 As shown, this demonstrates the direct use of Figure 4 The results of the hybrid domain multi-scale elastic wave full waveform inversion are shown. The L_BFGS algorithm was used to optimize and accelerate the inversion process. It can be seen that thanks to the use of the multi-scale method, the shallow structure of the model was restored to a certain extent. However, due to the lack of low frequencies and the poor initial model, the conventional hybrid domain multi-scale elastic wave full waveform inversion is severely affected by the periodic jump phenomenon. Therefore, the inversion effect is poor in the deep part of the model and near the fracture surface indicated by the arrow. Since the transverse wave velocity model inversion is more dependent on low-frequency information and the initial model, its inversion effect is even worse when low frequencies are missing and the initial model is poor.

[0082] like Figure 6 As shown, a velocity model containing large-scale structures obtained from time-frequency domain energy spectrum inversion is presented. During the time-frequency domain energy spectrum inversion experiment, the L_BFGS algorithm was used for optimization and acceleration, demonstrating a significant improvement compared to... Figure 4 The initial model provided by the present invention based on the multi-scale elastic wave full waveform inversion of the time-frequency domain energy spectrum is more accurate. It recovers the large-scale structure of the model in the case of low-frequency missing, which provides a good initial model for subsequent hybrid domain multi-scale inversion. Furthermore, the present invention does not introduce an additional low-pass filter in the inversion process, which reduces the workload and the difficulty of the work.

[0083] like Figure 7 As shown, the results of further hybrid-domain multi-scale elastic wave full-waveform inversion based on the natural frequency division time-frequency domain energy spectrum inversion are presented. It can be seen that the fracture surface, high-velocity layer, etc., in the model are effectively reconstructed. Compared to... Figure 5 , Figure 7 The inversion accuracy and resolution are higher within the same region, and it can be seen that the present invention significantly improves the inversion effect on deeper parts of the model. This is mainly due to the better initial model provided by the natural frequency division time-domain multi-scale energy spectrum inversion used in the present invention, which avoids the period jump phenomenon to a certain extent. The present invention improves the inversion effect of both P-wave velocity model and S-wave velocity model with missing low frequencies, and the improvement in S-wave inversion effect is particularly significant.

[0084] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art can make other variations or modifications based on the above description. It is neither necessary nor possible to exhaustively describe all embodiments here. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the claims of the present invention.

Claims

1. An elastic wave inversion method based on the natural frequency division time-frequency domain energy spectrum, characterized in that, Includes the following steps: S1: Input the measured three-component seismic record, read the field observation system from the track head of the measured data, read the seismic wavelet, inversion parameters and the initial P-wave velocity model mvp1(0) and S-wave velocity model mvs1(0), and perform grid subdivision of the computational space based on the inversion parameters; S2: Perform forward modeling of the elastic wave equation using the set parameters and the current model to obtain the simulated seismic records currecordx and currecordz of the x and z components under the current velocity model; S3: Read in the measured seismic records orirecordx and orirecordz for the x and z components, and determine the frequency range for time-frequency domain inversion. , Setting the inversion frequency interval Number of inversion frequency points The single-frequency S-transform of the simulated seismic records of the x and z components is used to obtain curx and curz, and the single-frequency S-transform of the measured seismic records of the x and z components is used to obtain orix and oriz. The squares of the amplitudes of curx, curz, orix, and oriz are taken respectively to obtain their respective single-frequency energy spectra. S4: Using the time-frequency domain energy spectrum elastic wave full waveform inversion with natural frequency division, multi-scale inversion is performed to obtain the longitudinal and transverse wave layer velocity model containing large-scale information; S5: Using the velocity model containing large-scale information as the initial model, perform hybrid domain multi-scale elastic wave full waveform inversion to obtain the final high-precision inversion result; In S4, the specific implementation method is as follows: Assuming the current frequency point is f, and the maximum number of iterations at the current frequency point is Niter1, the initial longitudinal wave model mvp1(j,f) and the initial transverse wave model mvs1(j,f) for the j-th iteration at the f-th frequency point are obtained according to S3, and the single-frequency energy spectra of curx, curz, orix, and oriz are obtained respectively. , , , The associated sources adjx and adjz of the full waveform inversion are calculated as follows: ； in, In order to seek the truth, Indicates the recording time. Represents the imaginary unit. Indicates a time delay; The adjoint wave field is obtained by backpropagating adjx and adjz using the elastic wave equation. The gradient is then calculated using the adjoint wave field and the simulated wave field according to gradient-based or Newton-based methods. This determines the P-wave update direction (directionvp1) and the S-wave update direction (directionvs1) for this iteration, and finally, the update step size. a ; Longitudinal wave model update formula: mvp1(j+1,f)=mvp1(j,f)+ a *directionvp1; The shear wave model update formula is mvs1(j+1,f)=mvs1(j,f)+ a *directionvs1; Update the velocity model, and determine if the current iteration number j of frequency point f is equal to Niter1. If it is, proceed to the next frequency point; otherwise, proceed to the next iteration of the current frequency point. If the values ​​are equal, output the final P-wave velocity model and S-wave velocity model; if they are not equal, then f = f + 1, j = 0, and proceed to the iteration of the next frequency point.

Images