Primary wave and multiple wave combined full-waveform inversion method

Abstract

The invention relates to the technical field of geological exploration, in particular to a primary wave and multiple combined full-waveform inversion method, which comprises the following steps of: based on an initial velocity model, performing numerical simulation by adopting a virtual surface source to obtain multiple simulation records of different orders; distinguishing according to a multiple order in the analog data, and constructing a multiple full-waveform inversion objective function of the corresponding order; an inversion gradient is calculated based on the objective function, and only the same-order multiples are matched and inversion is executed; sequentially carrying out inversion according to multiple orders from high to low, and taking a result of a previous order as an initial model of a next order; taking a multiple wave inversion result as a primary wave inversion initial model, and executing primary wave full-waveform inversion; and when a convergence condition is satisfied, outputting a final underground velocity model. According to the method, the dependence on the initial velocity model is reduced, and the stability and convergence of the inversion process and the construction precision of the underground velocity model are improved at the same time.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geological exploration, and particularly relates to a method for joint full waveform inversion of primary waves and multiple waves. BACKGROUND

[0002] As a high-resolution seismic inversion method based on wave equation, full waveform inversion can comprehensively utilize various waveform information such as amplitude, phase and travel time in seismic records to realize fine construction of the velocity model of underground medium. In theory, full waveform inversion can fully exploit all wave information contained in seismic data, and has significant advantages in seismic imaging and underground structure description, and has become one of the important technical means in the field of seismic exploration.

[0003] In practical application, when there is a large deviation between the initial velocity model and the real underground structure, full waveform inversion is prone to difficulty in waveform matching, and the inversion process is difficult to stably converge, especially in the case of multiple wave development, the complexity of different order multiple waves in the propagation path and reflection mechanism is easy to introduce crosstalk artifacts in the imaging or inversion process, thereby reducing the reliability of the inversion result. In order to avoid the above problems, the prior art usually removes multiple waves in advance and uses only primary wave data for inversion, but this approach has the problem of insufficient illumination in deep targets or complex geological structure areas, which limits the construction accuracy of the underground velocity model. SUMMARY

[0004] The present application provides a method for joint full waveform inversion of primary waves and multiple waves, by distinguishing the order of multiple waves and orderly introducing into the inversion process, on the premise of avoiding mutual crosstalk of multiple waves of different orders, progressive joint of multiple wave full waveform inversion and primary wave full waveform inversion is carried out, so as to improve the stability and convergence of the full waveform inversion process, and effectively improve the construction accuracy of the underground velocity model.

[0005] A method for joint full waveform inversion of primary waves and multiple waves, comprising the following steps: S1, in the initial stage of inversion, based on the initial velocity model, using a virtual surface source to numerically simulate the current velocity model, and obtaining simulated seismic records of multiple waves of different orders; S2, according to the order of multiple waves in the simulated seismic records, multiple wave full waveform inversion objective functions corresponding to the orders are respectively constructed; S3, based on the multiple wave full waveform inversion objective function, the inversion gradient is calculated, and the multiple wave full waveform inversion is executed under the condition of matching only multiple waves of the same order; S4, in order of multiple wave order from high to low, multiple wave full waveform inversion of each order is executed in turn, and the velocity model obtained by the previous order multiple wave full waveform inversion is taken as the initial model of the next order multiple wave full waveform inversion. S5, after the full waveform inversion of each order multiple wave is completed, taking the obtained velocity model as the initial model of the full waveform inversion of primary wave, performing the full waveform inversion of primary wave; S6, when the inversion process meets the preset convergence condition, outputting the final subsurface velocity model to obtain the joint full waveform inversion result of primary wave and multiple wave.

[0006] Optionally, the S1 comprises: S11, under the condition of the current velocity model, selecting the last order multiple wave observation record as the virtual surface source input, and loading the multiple wave observation record to the spatial position of the corresponding receiver point in time to construct the virtual source condition for the numerical simulation of multiple wave; S12, under the virtual source condition, based on the wave field propagation operator corresponding to the current velocity model, forward propagating the wave field excited by the virtual surface source to simulate the propagation process of the multiple wave field in the subsurface medium; S13, receiving the simulated multiple wave field at the receiver point position to obtain the next order multiple wave simulation seismic record under the constraint of the current velocity model, which is represented as: ; Wherein, is the receiver constraint operator, which is related to the spatial position of the receiver point, is the order multiple wave observation record (the zero order multiple wave is regarded as the primary wave), is the wave field propagation operator under the current velocity model, is the order multiple wave simulation data.

[0007] Optionally, the S2 comprises: S21, distinguishing the obtained multiple wave simulation seismic record according to the order of multiple wave to determine the multiple wave observation record corresponding to the order; S22, for each order multiple wave, respectively constructing the corresponding multiple wave full waveform inversion objective function for measuring the difference between the simulation data of the order multiple wave and the observation multiple wave data; S23, in the inversion process, only the simulation record and the corresponding observation record of the same order multiple wave are matched.

[0008] Optionally, the multiple wave full waveform inversion objective function is represented as: ; Wherein, is the full waveform inversion objective function corresponding to the order multiple wave, represents the two norm.

[0009] Optionally, the S3 comprises: S31, calculating an inversion gradient of the multiple wave full waveform inversion objective function with respect to the current velocity model; S32, in the inversion gradient calculation process, only introducing the simulated record and the observed record corresponding to the first order multiple wave as the matching object; S33, updating the current velocity model based on the inversion gradient, completing the full waveform inversion iteration for the first order multiple wave.

[0010] Optionally, the inversion gradient is expressed as: ; Wherein, is the gradient of the first order multiple wave full waveform inversion objective function with respect to the current velocity model, is the partial derivative of the first order multiple wave full waveform inversion objective function with respect to the current velocity model, is the data information residual of the first order multiple wave, is the sensitivity term of the first order multiple wave simulation record to the current velocity model, represents the cross-correlation calculation, represents the wave field reverse propagation operator.

[0011] Optionally, the S4 comprises: S41, taking the initial subsurface velocity model as the starting model of the multiple wave full waveform inversion, determining the highest order multiple wave participating in the multiple wave full waveform inversion; S42, taking the objective function of the current to be inverted order as the constraint condition, calling the corresponding full waveform inversion update operator to update the current velocity model, obtaining the velocity model after the full waveform inversion of the order multiple wave; S43, after completing the full waveform inversion of the current order multiple wave, taking the obtained velocity model as the initial model of the next adjacent low order multiple wave full waveform inversion, sequentially executing the multiple wave full waveform inversion of each order according to the order of the multiple wave from high to low; S44, after completing all the multiple wave orders participating in the inversion, outputting the subsurface velocity model after the full order multiple wave full waveform inversion, expressed as: ; Wherein, is the subsurface velocity model after the full order multiple wave full waveform inversion, is the initial subsurface velocity model, , ,​ They are respectively the 1st, 2nd, The objective function corresponding to the first-order multiple wave is the full waveform inversion update operator with constraints.

[0012] Optionally, S5 includes: S51. After completing the full waveform inversion of each order of multiple waves, the obtained underground velocity model is used as the initial model for the full waveform inversion of the first wave, and the source wavelet is determined as the forward modeling input. S52, the source wavelet is subjected to wavefield extension, and a detector constraint is applied to extract the wavefield response, resulting in a primary wave simulated seismic record, expressed as: ; in, For the source wavelet, This is a simulated earthquake record for a single wave. S53, Match the simulated first-wave seismic record with the observed first-wave data to construct the objective function for the full waveform inversion of the first wave; S54. Calculate the inversion gradient based on the objective function of the first wave full waveform inversion, and iteratively update the underground velocity model to complete the first wave full waveform inversion.

[0013] Optionally, the convergence condition of the multi-wave full waveform inversion objective function is that when the change in the multi-wave full waveform inversion objective function is less than the multi-wave threshold in two adjacent iterations, the inversion process is determined to have reached the convergence condition, or when the number of inversion iterations of the multi-wave full waveform inversion objective function reaches the maximum number of multi-wave iterations, the inversion process is terminated. The convergence condition of the primary wave full waveform inversion objective function is as follows: when the change of the primary wave full waveform inversion objective function is less than the primary wave threshold in two adjacent iterations, the inversion process is determined to have reached the convergence condition; or when the number of inversion iterations of the primary wave full waveform inversion objective function reaches the maximum number of primary wave iterations, the inversion process is terminated.

[0014] The beneficial effects of this invention are: This invention distinguishes multiples by order during the full waveform inversion process and constructs corresponding objective functions for full waveform inversion of multiples of the same order. It only matches the simulated records and observation records of multiples of the same order, thereby effectively avoiding crosstalk between multiples of different orders. This allows the complex propagation path information contained in the multiples to be introduced into the inversion process in an orderly manner, improving the stability of full waveform inversion of multiples and the physical consistency of the inversion results.

[0015] This invention performs full waveform inversion of multiple waves sequentially from high to low order, and uses the velocity model obtained from the previous order multiple wave inversion as the initial model for the next order multiple wave inversion. This gradually guides the velocity model to approximate the actual underground structure, thereby fully utilizing the advantages of higher-order multiple waves, which have a wider illumination range and higher sensitivity to deep and complex structural areas. This effectively improves the problem of insufficient illumination of primary waves in deep targets and complex geological conditions, and enhances the construction accuracy of the underground velocity model.

[0016] This invention uses the velocity model obtained after completing the full waveform inversion of each order of multiple waves as the initial model for the full waveform inversion of the first wave. This allows the first wave inversion to be carried out under more reasonable initial model conditions, reducing the dependence of the full waveform inversion on the accuracy of the initial velocity model, reducing the risk of getting trapped in local minima, and improving the convergence and stability of the inversion process. As a result, a more reliable subsurface velocity model is obtained, providing a high-quality basic model for subsequent seismic imaging and geological interpretation. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only for this invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 This is a schematic diagram of the inversion method according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the actual parameter model of an embodiment of the present invention; Figure 3 This is a schematic diagram of the initial parameter model of an embodiment of the present invention; Figure 4 This is a schematic diagram of the parameter model for the traditional first-wave full waveform inversion result according to an embodiment of the present invention; Figure 5 This is a schematic diagram of the parameter model for the second-order multiple full waveform inversion result according to an embodiment of the present invention; Figure 6 This is a schematic diagram of the parameter model for the second-order and first-order multiple full waveform inversion results according to an embodiment of the present invention; Figure 7 This is a schematic diagram of the parameter model of the joint full waveform inversion result of the primary wave and multiple waves according to an embodiment of the present invention; Figure 8 This is a schematic diagram showing the single-track speed comparison results of an embodiment of the present invention. Detailed Implementation

[0019] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. Those skilled in the art may employ other alternative methods to implement some well-known technologies; moreover, the accompanying drawings are only for more specific description of the embodiments and are not intended to specifically limit the present invention.

[0020] like Figures 1-8 As shown, a method for joint full waveform inversion of primary and multiple waves includes the following steps: S1, In the initial stage of inversion, based on the initial velocity model, a virtual surface source is used to perform numerical simulation on the current velocity model to obtain simulated seismic records of different orders of multiple waves. S2, according to the order of multiple waves in the simulated seismic record, construct the objective function for the full waveform inversion of the corresponding order of multiple waves; S3, calculate the inversion gradient based on the objective function of the multiple full waveform inversion, and perform the multiple full waveform inversion under the condition that only multiples of the same order are matched; S4. Perform the full waveform inversion of each order of multiple waves in descending order, and use the velocity model obtained from the full waveform inversion of the previous order of multiple waves as the initial model for the full waveform inversion of the next order of multiple waves. S5. After completing the full waveform inversion of each order of multiple waves, the obtained velocity model is used as the initial model for the full waveform inversion of the first wave, and the full waveform inversion of the first wave is performed. S6. When the inversion process meets the preset convergence condition, the final underground velocity model is output, and the combined full waveform inversion results of primary and secondary waves are obtained.

[0021] S1 includes: S11. Under the current velocity model conditions, the previous order multiple wave observation record is selected as the virtual surface source input, and the multiple wave observation record is sequentially loaded to the spatial position of the corresponding receiver point to construct the virtual source conditions for multiple wave numerical simulation. S12, under the condition of virtual seismic source, based on the wave field propagation operator corresponding to the current velocity model, the wave field excited by the virtual surface source is positively extended to simulate the propagation process of multiple wave fields in the underground medium. S13, at the receiver location, the simulated multiple wave field is received to obtain the next-order multiple simulated seismic record under the constraints of the current velocity model, represented as: ; in, This is a detector constraint operator, which is related to the spatial location of the detector point. for Observational records of first-order multiples (zero-order multiples are considered as first-order multiples). For the wave field propagation operator under the current velocity model, for Multiple-wave simulation data.

[0022] S2 includes: S21, the obtained multi-wave simulated seismic records are distinguished according to the multi-wave order, and the observation records of the corresponding order of multi-wave are determined. S22, For each order of multiples, a corresponding objective function for the full waveform inversion of multiples is constructed to measure the difference between the simulated data and the observed multiples data of that order. S23, during the inversion process, only the simulated records of the same order multiple waves are matched with the corresponding observation records.

[0023] The objective function for multi-wave full waveform inversion is expressed as: ; in, For the first The objective function for full waveform inversion corresponding to the first-order multiple wave. This represents the L2 norm.

[0024] S3 includes: S31, Calculate the inversion gradient of the objective function for the multiple wave full waveform inversion with respect to the current velocity model, used to characterize the effect of the current velocity model on the first wave. Sensitivity to mismatch in multiple wavelet data; S32, in the process of inverting gradient calculation, only the first... The simulated and observed records corresponding to the first multiple wave were used as matching objects; S33, update the current velocity model based on the inverted gradient, completing the update for the first... Full waveform inversion iteration of multiple waves.

[0025] The inversion gradient is expressed as: ; in, For the first The objective function for the full waveform inversion of the first-order multiple wave with respect to the current velocity model gradient, For the first The partial derivatives of the objective function for full waveform inversion of the first-order multiple wave with respect to the current velocity model. For the first residual information of multiple wave data, For the first The first-order multiple-wave simulation records the sensitivity terms to the current velocity model. Indicates cross-correlation calculation, This represents the wave field backpropagation operator.

[0026] S4 includes: S41, using the initial underground velocity model as the starting model for the full waveform inversion of multiples, the highest-order multiples participating in the full waveform inversion of multiples are determined. S42, with the objective function of the current order to be inverted as the constraint, call the corresponding full waveform inversion update operator to update the current velocity model and obtain the velocity model after the full waveform inversion of the current order; S43, After completing the full waveform inversion of the current order multiple, the obtained velocity model is used as the initial model for the next adjacent low-order multiple full waveform inversion, and the full waveform inversion of each order of multiple is executed in order from high to low. S44, after completing the inversion of all wave orders, outputs the subsurface velocity model that completes the full waveform inversion of all wave orders, expressed as: ; in, To complete the subsurface velocity model for full waveform inversion of all orders of multiple waves, For the initial underground velocity model, , , They are respectively the 1st, 2nd, The objective function corresponding to the first-order multiple wave is the full waveform inversion update operator with constraints.

[0027] S5 includes: S51. After completing the full waveform inversion of each order of multiple waves, the obtained underground velocity model is used as the initial model for the full waveform inversion of the first wave, and the source wavelet is determined as the forward modeling input. S52, wavefield extension is performed on the source wavelet, and a detector constraint is applied to extract the wavefield response, resulting in a primary wave simulated seismic record, expressed as: ; in, For the source wavelet, This is a simulated earthquake record for a single wave. S53, matching the simulated first-wave seismic record with the observed first-wave data to construct the objective function for the full waveform inversion of the first wave; S54 calculates the inversion gradient based on the objective function of the first wave full waveform inversion and iteratively updates the underground velocity model to complete the first wave full waveform inversion.

[0028] The convergence condition of the objective function of the multi-wave full waveform inversion is that when the change of the objective function of the multi-wave full waveform inversion in two adjacent iterations is less than the multi-wave threshold, the inversion process is determined to have reached the convergence condition, or when the number of inversion iterations of the objective function of the multi-wave full waveform inversion reaches the maximum number of iterations of the multi-wave, the inversion process is terminated. The convergence condition of the objective function for the first-wave full waveform inversion is that the change in the objective function for the first-wave full waveform inversion is less than the first-wave threshold in two adjacent iterations, at which point the inversion process is considered to have reached the convergence condition, or the inversion process is terminated when the number of inversion iterations of the objective function for the first-wave full waveform inversion reaches the maximum number of iterations for the first wave.

[0029] The specific implementation method is as follows: To verify the superiority of the method of the present invention over the traditional first-wave full-wave inversion and its applicability under complex structural conditions, the present invention was applied to, for example... Figure 2 In the complex underground structural model shown, a mobile observation system with incomplete primary wave illumination was adopted. 117 seismic sources were set up for excitation and 201 receiver points were set up for reception. No seismic sources were deployed within a lateral range of 2 to 4 kilometers, thus artificially limiting the illumination range of the primary reflected waves and forming a distinct area of ​​weak primary wave illumination.

[0030] Under the aforementioned geological conditions and observation system, the embodiment first uses the traditional first-wave full waveform inversion method to perform... Figure 3 The initial velocity model shown is subjected to 30 iterations of inversion to obtain... Figure 4 The results of the full waveform inversion for the first wave are shown. Due to the limitations of the first wave illumination, the overall error of the inversion results is relatively large, with local velocity update errors and unclear thin-layer and fault interfaces. Especially in the weak illumination region of 3 to 4 kilometers, the velocity model can hardly be updated effectively.

[0031] To address the aforementioned issues, this invention incorporates the full waveform inversion of multiples into the traditional full waveform inversion process for first-order multiples by constructing objective functions corresponding to different orders of multiples. In this embodiment, starting with the full waveform inversion of second-order multiples, a first-order multiple is used as a virtual surface source, and simulated data of second-order multiples are obtained based on wavefield propagation operators. The objective function is defined in the least squares sense, with the criterion being minimizing the L2 error between the observed and simulated data of second-order multiples. To ensure consistency in the number of model updates during the inversion process, each inversion stage undergoes 10 iterations, for a total of 30 updates. The full waveform inversion results of second-order multiples are shown below. Figure 5 As shown, the inversion effect is significantly improved in the 3 to 4 km range.

[0032] The inversion results of the second-order multiple full waveform are used as the initial model for the inversion of the first-order multiple full waveform. An objective function corresponding to the first-order multiple is constructed and iterated, yielding the following inversion results: Figure 6As shown, the accuracy of velocity modeling is further improved. Finally, the full waveform inversion results of the multiple waves are used as the initial model for the full waveform inversion of the primary wave for iteration, obtaining the joint full waveform inversion results of the primary and multiple waves (e.g., Figure 7 (As shown). Compared with the results obtained by the traditional full-waveform inversion method using only the first wave, the joint full-waveform inversion method proposed in this invention can effectively compensate for the insufficient illumination of the first wave and provide richer reflection information constraints for the inversion process. Even with a relatively small number of iterations, the obtained velocity model is significantly closer to the true velocity model, reducing the dependence of the inversion results on the initial model to a certain extent.

[0033] To further demonstrate the advantages of the present invention, Figure 8 This paper presents a comparison of single-channel velocities obtained through traditional primary wave full waveform inversion and joint full waveform inversion. The velocity versus depth curves at lateral positions of 1 km, 2 km, and 3 km show that, under different primary wave illumination conditions, the joint inversion method generally provides more accurate velocity updates, and the results show significantly higher agreement with the actual velocity model.

[0034] This invention encompasses any substitutions, modifications, equivalent methods, and solutions made within the spirit and scope of this invention. To provide the public with a thorough understanding of this invention, specific details are described in detail in the following preferred embodiments; however, those skilled in the art will fully understand the invention even without these details. Furthermore, to avoid unnecessary misunderstanding of the essence of this invention, well-known methods, processes, procedures, components, and circuits are not described in detail.

[0035] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for joint full-waveform inversion of primary and multiple waves, characterized in that, Includes the following steps: S1, In the initial stage of inversion, based on the initial velocity model, a virtual surface source is used to perform numerical simulation on the current velocity model to obtain simulated seismic records of different orders of multiple waves. S2, according to the order of multiple waves in the simulated seismic record, construct the objective function for the full waveform inversion of the corresponding order of multiple waves; S3, calculate the inversion gradient based on the objective function of the multiple full waveform inversion, and perform the multiple full waveform inversion under the condition that only multiples of the same order are matched; S4. Perform the full waveform inversion of each order of multiple waves in descending order, and use the velocity model obtained from the full waveform inversion of the previous order of multiple waves as the initial model for the full waveform inversion of the next order of multiple waves. S5. After completing the full waveform inversion of each order of multiple waves, the obtained velocity model is used as the initial model for the full waveform inversion of the first wave, and the full waveform inversion of the first wave is performed. S6. When the inversion process meets the preset convergence condition, the final underground velocity model is output, and the combined full waveform inversion results of the primary wave and multiple waves are obtained.

2. The method for joint full waveform inversion of primary and multiple waves according to claim 1, characterized in that, S1 includes: S11. Under the current velocity model conditions, the previous order multiple wave observation record is selected as the virtual surface source input, and the multiple wave observation record is sequentially loaded to the spatial position of the corresponding receiver point to construct the virtual source conditions for multiple wave numerical simulation. S12, under the virtual source conditions, based on the wave field propagation operator corresponding to the current velocity model, the wave field excited by the virtual surface source is forward extended to simulate the propagation process of multiple wave fields in the underground medium. S13, at the receiver location, the simulated multiple wave field is received to obtain the next-order multiple simulated seismic record under the constraints of the current velocity model, represented as: ； in, For the detector constraint operator, for Multiple wave observation records, For the wave field propagation operator under the current velocity model, for Multiple-wave simulation data.

3. The method for joint full waveform inversion of primary and multiple waves according to claim 2, characterized in that, S2 includes: S21, the obtained multi-wave simulated seismic records are distinguished according to the multi-wave order, and the observation records of the corresponding order of multi-wave are determined. S22, For each order of multiples, a corresponding objective function for the full waveform inversion of multiples is constructed to measure the difference between the simulated data and the observed multiples data of that order. S23, during the inversion process, only the simulated records of the same order multiple waves are matched with the corresponding observation records.

4. The method for joint full waveform inversion of primary and multiple waves according to claim 3, characterized in that, The objective function for the multiple wave full waveform inversion is expressed as: ； in, For the first The objective function for full waveform inversion corresponding to the first-order multiple wave. This represents the L2 norm.

5. The method for joint full waveform inversion of primary and multiple waves according to claim 4, characterized in that, S3 includes: S31, Calculate the inversion gradient of the objective function for the multiple full waveform inversion with respect to the current velocity model; S32, in the process of inverting gradient calculation, only the first... The simulated and observed records corresponding to the first multiple wave were used as matching objects; S33, Update the current velocity model based on the inverted gradient, completing the update for the first... Full waveform inversion iteration of multiple waves.

6. The method for joint full waveform inversion of primary and multiple waves according to claim 5, characterized in that, The inversion gradient is expressed as: ； in, For the first The objective function for the full waveform inversion of the first-order multiple wave with respect to the current velocity model gradient, For the first The partial derivatives of the objective function for full waveform inversion of the first-order multiple wave with respect to the current velocity model. For the first residual information of multiple wave data, For the first The first-order multiple-wave simulation records the sensitivity terms to the current velocity model. Indicates cross-correlation calculation, This represents the wave field backpropagation operator.

7. The method for joint full waveform inversion of primary and multiple waves according to claim 6, characterized in that, S4 includes: S41, using the initial underground velocity model as the starting model for the full waveform inversion of multiples, the highest-order multiples participating in the full waveform inversion of multiples are determined. S42, with the objective function of the current order to be inverted as the constraint, call the corresponding full waveform inversion update operator to update the current velocity model and obtain the velocity model after the full waveform inversion of the current order; S43, After completing the full waveform inversion of the current order multiple, the obtained velocity model is used as the initial model for the next adjacent low-order multiple full waveform inversion, and the full waveform inversion of each order of multiple is executed in order from high to low. S44, after completing the inversion of all wave orders, outputs the subsurface velocity model that completes the full waveform inversion of all wave orders, expressed as: ； in, To complete the subsurface velocity model for full waveform inversion of all orders of multiple waves, For the initial underground velocity model, , , They are respectively the 1st, 2nd, The objective function corresponding to the first-order multiple wave is the full waveform inversion update operator with constraints.

8. The method for joint full waveform inversion of primary and multiple waves according to claim 7, characterized in that, S5 includes: S51. After completing the full waveform inversion of each order of multiple waves, the obtained underground velocity model is used as the initial model for the full waveform inversion of the first wave, and the source wavelet is determined as the forward modeling input. S52, the source wavelet is subjected to wavefield extension, and a detector constraint is applied to extract the wavefield response, resulting in a primary wave simulated seismic record, expressed as: ； in, For the source wavelet, This is a simulated earthquake record for a single wave. S53, Match the simulated first-wave seismic record with the observed first-wave data to construct the objective function for the full waveform inversion of the first wave; S54. Calculate the inversion gradient based on the objective function of the first wave full waveform inversion, and iteratively update the underground velocity model to complete the first wave full waveform inversion.

9. The method for joint full waveform inversion of primary and multiple waves according to claim 8, characterized in that, The convergence condition of the multi-wave full waveform inversion objective function is as follows: when the change of the multi-wave full waveform inversion objective function in two adjacent iterations is less than the multi-wave threshold, the inversion process is determined to have reached the convergence condition; or when the number of inversion iterations of the multi-wave full waveform inversion objective function reaches the maximum number of multi-wave iterations, the inversion process is terminated. The convergence condition of the primary wave full waveform inversion objective function is as follows: when the change of the primary wave full waveform inversion objective function is less than the primary wave threshold in two adjacent iterations, the inversion process is determined to have reached the convergence condition; or when the number of inversion iterations of the primary wave full waveform inversion objective function reaches the maximum number of primary wave iterations, the inversion process is terminated.

Images