Early arrival wave equation travel time inversion method for optimizing error function

Abstract

The invention provides an early-arrival wave equation travel time inversion method for optimizing an error function. The early-arrival wave equation travel time inversion method for optimizing the error function comprises the following steps: step 1, inputting field observation early-arrival wave seismic data, seismic wavelets and an initial speed model; step 2, carrying out numerical solution on the wave equation to obtain a seismic source wave field and simulated seismic data; 3, calculating the travel time difference between the observed seismic data and the simulated seismic data, and constructing an accompanying seismic source; 4, the wave equation is subjected to numerical solution, and an accompanying wave field is obtained; step 5, calculating a gradient and obtaining a speed update quantity; step 6, constructing an optimized error function according to the travel time difference and the seismic data; and step 7, when the inversion result converges, outputting the inversion result. According to the early arrival wave equation travel time inversion method for optimizing the error function, the defect that wave equation inversion is prone to falling into a local minimum value is overcome, and inversion convergence and speed inversion precision are improved.

Description

Technical Field

[0001] This invention relates to the field of seismic exploration technology, and in particular to a method for travel-time inversion of early arrival wave equations with optimized error functions. Background Technology

[0002] Early-arrival wave equation travel-time inversion, based on wave theory, establishes an objective function by matching the travel-time information of simulated early-arrival wave seismic data with that of field-observed early-arrival wave seismic data. Iteratively finding the global minimum of the objective function using local optimization methods allows for the high-precision recovery of seafloor geological velocity parameters, representing one of the cutting-edge directions in seismic exploration. However, wave equation inversion exhibits strong nonlinearity. When the initial velocity is not accurate enough, seismic data matching is prone to cycle skipping, and the error function is susceptible to getting trapped in local minima, leading to inversion failure. Therefore, prioritizing the matching of seismic data travel-time information and optimizing the error function based on the differences in travel-time between the matches not only reduces the likelihood of the error function getting trapped in local minima but also improves the convergence speed and accuracy of the inversion.

[0003] The full waveform inversion method has extremely strong nonlinearity. The conventional early arrival wave equation inversion matches the waveform information of the early arrival wave seismic data observed in the field and the simulated early arrival wave seismic data. When the initial velocity is not accurate enough, the seismic data matching is prone to cycle skipping, and the error function is prone to getting trapped in local minima, leading to inversion failure.

[0004] Chinese patent application CN109085643A discloses a step-by-step joint inversion method for early arrival waves, comprising: inputting observation data and an initial velocity model; defining a travel-time tomography objective functional expression and obtaining the gradient expression of the objective functional with respect to the model using the adjoint state method; updating the model background field using the obtained travel-time tomography gradient; using the obtained accurate model background field as the initial model for full waveform inversion; establishing the objective functional for full waveform inversion using the obtained velocity model and obtaining the gradient expression of the objective functional with respect to the model using the adjoint state method; updating the model background field using the obtained full waveform inversion gradient; and performing full waveform inversion using the obtained accurate model background field to obtain high-precision velocity field modeling results. However, this application does not disclose the specific form of the error function. Although a step-by-step inversion strategy is adopted, its convergence and stability may still face challenges under certain complex geological conditions.

[0005] Chinese patent application CN103913768A discloses a method and apparatus for modeling shallow subsurface layers based on seismic wave data. The method includes: preprocessing raw seismic wave data to obtain a first model; determining the propagation time and offset of an early-arriving wave based on the first model and the depth of interest; extracting observational data about the early-arriving wave from the raw data based on the propagation time and offset; generating calculated data about the early-arriving wave based on the first model and the determined propagation time ts and offset; performing data domain fitting on the observational and calculated data to obtain a gradient, thereby further obtaining a velocity update model; superimposing the first model and the velocity update model to obtain a second model; repeating the above steps until the obtained velocity update model reaches the global optimum. By using only partial early-arriving wave information, the nonlinearity of the objective function is reduced; and by employing small offset information, error accumulation is reduced. Therefore, this invention is more likely to converge to the global optimum. However, this application focuses on using seismic wave data to model shallow subsurface layers, which is highly dependent on the accuracy of the initial model. If the initial model differs significantly from the actual geological structure, it may affect the convergence of the iterative process and the accuracy of the final model.

[0006] Chinese patent application CN105549080A discloses a method for inverting undulating surface waveforms based on an auxiliary coordinate system, belonging to the field of petroleum geophysical exploration technology. This invention first employs an early-arrival wave waveform inversion method under an auxiliary coordinate system to invert the curved grid region to obtain updated near-surface velocities. Then, it uses a full-waveform inversion method under the auxiliary coordinate system to update the global velocity field, overcoming the influence of inaccurate near-surface velocities on deep velocity field inversion. Finally, it employs a time-domain multi-scale inversion method to invert velocities from low to high frequencies, overcoming the dependence of waveform inversion methods on the initial velocity model. However, the introduction of an auxiliary coordinate system and curved grid partitioning in this application may increase the complexity of the algorithm and potentially require higher computational resources and time.

[0007] The existing technologies described above are significantly different from this invention. A search reveals no literature in the XY category, indicating that this invention possesses novelty. Since existing technologies lack solutions to the technical problem we aim to address, we have invented a novel method for travel time inversion of the early arrival wave equation with an optimized error function. Summary of the Invention

[0008] The purpose of this invention is to provide an early-arrival wave equation travel-time inversion method that solves the defect of wave equation inversion easily falling into local minima, and improves the convergence and velocity inversion accuracy of the inversion by optimizing the error function.

[0009] The objective of this invention can be achieved through the following technical measures: a method for travel time inversion of the early-arrival wave equation with optimized error function, the method comprising:

[0010] Step 1: Input early arrival wave seismic data, seismic wavelet, and initial velocity model from field observations;

[0011] Step 2: Numerically solve the wave equation to obtain the source wavefield and simulated earthquake data;

[0012] Step 3: Calculate the travel time difference between observed seismic data and simulated seismic data and construct the associated seismic source;

[0013] Step 4: Numerically solve the wave equation to obtain the adjoint wave field;

[0014] Step 5: Calculate the gradient and obtain the velocity update amount;

[0015] Step 6: Construct an optimized error function based on the travel time difference and seismic data;

[0016] Step 7: When the inversion results converge, output the inversion results.

[0017] The objective of this invention can also be achieved through the following technical measures:

[0018] In step 1, input the early arrival wave seismic data d from field observations. obs (x,t), seismic wavelet f(t), and initial velocity parameter model m0; the initial velocity parameter model m0 is known, and its velocity value gradually increases from shallow to deep layers.

[0019] In step 2, based on the early arrival wave seismic data observed in the field, d obs (x,t) is used to extract the seismic wavelet f(t). Using the seismic wavelet f(t) as the source, the time-domain wave equation is numerically solved to obtain the source wavefield u(x,t) at all time points. Data at the receiving point is taken from the forward modeling simulation wavefield u(x,t), i.e., the simulated early arrival wave seismic data d. cal (x,t).

[0020] In step 2, the time-domain wave equation is expressed as follows:

[0021]

[0022] In the formula, u(x,t) represents the forward modeling wave field, x represents the spatial position vector, m represents the velocity parameter model, and the initial value is m0; It represents the sum of spatial derivatives.

[0023] In step 3, the travel time difference between the early arrival wave seismic data observed in the field and the simulated early arrival wave seismic data is calculated and denoted as Δτ(x,t). Based on this, the source of the accompanying wave field is constructed, i.e., the accompanying source.

[0024] In step 3, the accompanying earthquake source η(x,t) is represented as follows:

[0025]

[0026] In the formula, and Let represent the first and second time derivatives of the early arrival wave seismic data observed in the field after time shifting, respectively. The early arrival wave seismic data observed in the field after time shifting is consistent with the simulated early arrival wave seismic data in time. Therefore, the early arrival wave seismic data observed in the field can be replaced by the simulated early arrival wave seismic data. Equation (2) can be further expressed as:

[0027]

[0028] In the formula, and Let represent the first and second time derivatives of the simulated early arrival wave seismic data, respectively. The main contributing terms are: and For a specific number, it represents a secondary contribution.

[0029] In step 4, the accompanying source is used as the source to replace the seismic wavelet in step 2, and then numerical forward modeling is performed again to obtain the accompanying wave field; that is, the accompanying source of equation (3) is used as the source to replace the seismic wavelet f(t) in equation (1), and then the time domain wave equation (2) is numerically solved again to obtain the accompanying wave field q(x,t).

[0030] In step 5, the gradient operator is calculated based on the source wave field and the associated wave field, and an appropriate step size is obtained to apply to the gradient operator to calculate the velocity parameter model update amount for this iteration.

[0031] In step 5, the velocity update amount for each iteration is:

[0032] The gradient operator is calculated based on the source wavefield u(x,t) and the adjoint wavefield q(x,t).

[0033]

[0034] In the formula, s and g represent the source point and receiver point of the early arrival wave seismic data generated in the field observation, respectively, and t is the time when the early arrival wave arrives at the receiver point; first, a trial step size is given, and then the appropriate step size is obtained by the one-dimensional line search method under the Armijo condition. The gradient operator is applied to calculate the velocity parameter model update amount of this iteration, and the updated velocity parameter model m is calculated using equation (5):

[0035] m=m0+λ·▽E(m)(11)

[0036] In the formula, m0 represents the initial velocity parameter model, and λ represents the step size. This indicates the amount of time the velocity parameter model has been updated in this iteration.

[0037] In step 6, an optimized error function, Error, is constructed based on the main contributions of the accompanying epicenter.

[0038]

[0039] Where Δτ(x,t) represents the travel time difference between the early arrival wave seismic data observed in the field and the simulated early arrival wave seismic data. denoted by the first time derivative of the simulated early arrival wave seismic data, s and g represent the source point and receiver point that generated the early arrival wave seismic data from field observations, respectively, and t is the time when the early arrival wave arrives at the receiver point.

[0040] In step 7, if the inversion result does not converge, return to step 2 to continue iterative updates.

[0041] The objective of this invention can also be achieved through the following technical measures: an early arrival wave equation travel-time inversion system with optimized error function, wherein the early arrival wave equation travel-time inversion system with optimized error function uses the early arrival wave equation travel-time inversion method with optimized error function to invert the seismic wave propagation velocity model.

[0042] The objective of this invention can also be achieved through the following technical measures: a computer-readable storage medium comprising a stored computer program; wherein, when the computer program is running, it controls the device where the computer-readable storage medium is located to execute an early arrival wave equation travel time inversion method with optimized error function.

[0043] The early arrival wave equation travel time inversion method of the present invention with optimized error function matches the travel time information of early arrival wave seismic data observed in the field and simulated early arrival wave seismic data, and reconstructs the associated source and error function of the full waveform inversion according to the travel time difference. This solves the defect of wave equation inversion easily falling into local minima and improves the convergence and velocity inversion accuracy. Attached Figure Description

[0044] Figure 1 This is a schematic diagram of a method for travel time inversion of early arrival wave equations with optimized error functions, as disclosed in an embodiment of the present invention.

[0045] Figure 2 This is a schematic diagram of the Marmous i2 velocity parameter model in a specific embodiment of the present invention;

[0046] Figure 3 This is a schematic diagram of the pseudo-well curve of the velocity parameter model in a specific embodiment of the present invention;

[0047] Figure 4 This is a schematic diagram of early arrival wave seismic data in a specific embodiment of the present invention;

[0048] Figure 5 This is a schematic diagram of the error function in a specific embodiment of the present invention;

[0049] Figure 6 This is a schematic diagram of the error convergence curve in a specific embodiment of the present invention;

[0050] Figure 7 This is a schematic diagram of the velocity parameter model results obtained in the 35th iteration of the early arrival wave equation inversion method with optimized error function in a specific embodiment of the present invention.

[0051] Figure 8 This is a schematic diagram comparing the pseudo-well curves of the velocity parameter model in a specific embodiment of the present invention;

[0052] Figure 9 This is a schematic diagram of an underground velocity parameter model for a specific work area according to a specific embodiment of the present invention;

[0053] Figure 10 This is a schematic diagram of early arrival wave seismic data for a certain work area in a specific embodiment of the present invention;

[0054] Figure 11 This is a schematic diagram of the error function in a specific embodiment of the present invention;

[0055] Figure 12 This is a schematic diagram of the error convergence curve in a specific embodiment of the present invention;

[0056] Figure 13 This is a schematic diagram of the velocity parameter model results obtained in the 17th iteration of the early arrival wave equation inversion method with optimized error function in a specific embodiment of the present invention.

[0057] Figure 14 This is a schematic diagram comparing the pseudo-well curves of the velocity parameter model in a specific embodiment of the present invention;

[0058] Figure 15 This is a schematic diagram of early arrival wave seismic data corresponding to the inversion velocity of a certain work area in a specific embodiment of the present invention. Detailed Implementation

[0059] like Figure 1 As shown, Figure 1 This is a flowchart of the early-arrival wave equation travel-time inversion method based on the optimized error function of the present invention. The early-arrival wave equation travel-time inversion method based on the optimized error function includes:

[0060] Step 1. Input the early arrival wave seismic data, seismic wavelet, and initial velocity parameter model from field observations.

[0061] Step 2. Using the seismic wavelet as the source, the wave field is simulated using the wave equation forward modeling to obtain early arrival wave seismic data.

[0062] Step 3. Calculate the travel time difference between the early arrival wave seismic data observed in the field and the simulated early arrival wave seismic data, and construct the source of the accompanying wave field.

[0063] Step 4. Replace the seismic wavelet in Step 2 with the accompanying source as the source, and then perform numerical forward modeling again to obtain the accompanying wave field.

[0064] Step 5. Calculate the gradient operator based on the source wave field and the associated wave field, and determine an appropriate step size to apply to the gradient operator to calculate the velocity parameter model update amount for this iteration.

[0065] Step 6. Construct an optimized error function as a judgment condition to determine whether the inversion result has converged, and then choose to return to step 2 to continue iterative updates or terminate the iteration and output the inversion result.

[0066] Conventional early-arrival wave equation inversion matches waveform information from both field-observed and simulated early-arrival seismic data. When the initial velocity is inaccurate, it is prone to getting trapped in local minima, leading to inversion failure. This invention prioritizes matching the travel time information of the seismic data and optimizes the error function based on the differences in travel time. This not only reduces the likelihood of the error function getting trapped in local minima but also improves the convergence speed and accuracy of the inversion.

[0067] The following are several specific embodiments of the application of the present invention.

[0068] Example 1

[0069] In a specific embodiment 1 of the present invention, the early-arrival wave equation travel-time inversion method for optimizing the error function includes the following steps:

[0070] Step 1, Input Data:

[0071] Input early arrival wave seismic data from field observations d obs (x,t), seismic wavelet f(t), and initial velocity parameter model m0; the initial velocity parameter model m0 is known, and its velocity value gradually increases from shallow to deep layers. In this embodiment, the Marmousi2 velocity parameter model is selected, such as... Figure 2 The initial velocity parameter model is shown in (b). The actual velocity parameter model is unknown, as shown in (b). Figure 2 As shown in (a) the true velocity parameter model, since the initial velocity model deviates significantly from the true velocity model, such as Figure 3As shown, the solid black line represents... Figure 2 The pseudo-well curve of the true velocity at a distance of 2.5km, with the black dashed line representing... Figure 2 The pseudo-well curve of the true velocity at a distance of 2.5 km in b shows that traditional wave equation inversion methods cannot directly obtain the true velocity parameter model. Field observations of early-arrival wave seismic data are as follows: Figure 4 As shown in (a).

[0072] Step 2, numerical calculation of the source wave field:

[0073] Based on early arrival wave earthquake data observed in the field d obs (x,t) is used to extract the seismic wavelet f(t). Using the seismic wavelet f(t) as the source, the time-domain wave equation is numerically solved to obtain the source wavefield u(x,t) at all time points. Data at the receiving point is taken from the forward modeling simulation wavefield u(x,t), i.e., the simulated early arrival wave seismic data d. cal (x,t), such as Figure 4 As shown in the simulated early arrival wave seismic data in (b), the time-domain wave equation is expressed as follows:

[0074]

[0075] In the formula, u(x,t) represents the forward modeling wave field, x represents the spatial position vector, m represents the velocity parameter model, and the initial value is m0; This represents the sum of spatial derivatives;

[0076] Step 3, constructing the accompanying seismic source:

[0077] The travel time difference between field-observed early-arrival wave seismic data and simulated early-arrival wave seismic data is calculated and denoted as Δτ(x,t). Based on this, the source of the accompanying wave field, i.e., the accompanying source, is constructed, and its representation is as follows:

[0078]

[0079] In the formula, Let represent the first and second time derivatives of the early arrival wave seismic data observed in the field after time shifting, respectively. The early arrival wave seismic data observed in the field after time shifting is consistent with the simulated early arrival wave seismic data in time. Therefore, the early arrival wave seismic data observed in the field can be replaced by the simulated early arrival wave seismic data, and equation (2) can be further expressed as:

[0080]

[0081] In the formula, Let represent the first and second time derivatives of the simulated early arrival wave seismic data, respectively. The main contributing terms are: and For a specific number, it represents a secondary contribution.

[0082] Step 4: Numerical calculation of the adjoint wave field:

[0083] The accompanying source of equation (3) is used as the source to replace the seismic wavelet f(t) in equation (1), and then the time domain wave equation (2) is solved numerically again to obtain the accompanying wave field q(x,t);

[0084] Step 5: Obtain the velocity update for each iteration:

[0085] The gradient operator is calculated based on the source wavefield u(x,t) and the adjoint wavefield q(x,t).

[0086]

[0087] In the formula, s and g represent the source point and receiver point of the early arrival wave seismic data generated in the field observation, respectively. First, a trial step size is given, and then the appropriate step size is obtained by using the one-dimensional line search method with Armijo conditions. The gradient operator is applied to calculate the velocity parameter model update amount for this iteration, ensuring that the velocity parameter model update amount for each iteration is between 20m / s and 100m / s. The updated velocity parameter model m is calculated using equation (5):

[0088]

[0089] In the formula, m0 represents the initial velocity parameter model, and λ represents the step size. This indicates the amount of time the velocity parameter model has been updated in this iteration;

[0090] Step Six: Construct an optimized error function and determine the convergence result:

[0091] An optimized error function is constructed based on the main contributions of the accompanying epicenter:

[0092]

[0093] The optimized error function not only considers the main contribution of the accompanying source but also avoids the problem of waveform errors easily getting caught in cycle skipping. Error functions for different iteration numbers are as follows: Figure 5 As shown, the part to the left of the white dashed line represents the error function of the first iteration, and the part to the right of the white dashed line represents the error function of the tenth iteration. As the inversion iterations proceed, the value of the error function gradually decreases, indicating successful convergence. Figure 6 As shown.

[0094] Returning to step two, the updated velocity parameter model m is substituted into equation (1) for numerical simulation. The simulated early arrival wave seismic data is recalculated. Equation (6) is used as the judgment condition to determine whether the inversion has converged. Then, the option is to continue iterative updates or terminate the iteration and output the inversion results.

[0095] Inversion results are shown Figure 7 The inversion of the early-arrival wave equation with optimized error function can effectively update the background velocity parameter model in the shallow and intermediate layers. A comparison of the pseudo-well curves from the inversion results is shown below. Figure 8 (a) Pseudo-well curve at a distance of 1.5 km, (b) Pseudo-well curve at a distance of 2.5 km, where the black solid line represents the true velocity, the black dashed line represents the initial velocity, and the black asterisk line represents the inversion velocity. It can be seen that the inversion velocity after the error converges smoothly is closer to the true velocity. When the initial velocity deviates far from the true velocity, the early arrival wave equation inversion with optimized error function can accurately recover the background velocity parameter model of the middle and shallow layers.

[0096] Example 2

[0097] In a specific embodiment 2 of the present invention, the early-arrival wave equation travel-time inversion method for optimizing the error function includes the following steps:

[0098] Step 1, Input Data:

[0099] Seismic data collected on land in a work area in eastern China was selected, and early arrival wave seismic data from field observations were input. obs (x,t), seismic wavelet f(t), and initial velocity parameter model m0; the initial velocity parameter model m0 is known, such as Figure 9 As shown. Field observation data for early arrival wave seismic waves are as follows: Figure 10 (a) Early arrival wave seismic data observed in the field.

[0100] Step 2, numerical calculation of the source wave field:

[0101] Based on early arrival wave earthquake data observed in the field d obs (x,t) is used to extract the seismic wavelet f(t). Using the seismic wavelet f(t) as the source, the time-domain wave equation is numerically solved to obtain the source wavefield u(x,t) at all time points. Data at the receiving point is taken from the forward modeling simulation wavefield u(x,t), i.e., the simulated early arrival wave seismic data d. cal (x,t), such as Figure 10 As shown in the simulated early arrival wave seismic data in (b), the time-domain wave equation is expressed as follows:

[0102]

[0103] In the formula, u(x,t) represents the forward modeling wave field, x represents the spatial position vector, m represents the velocity parameter model, and the initial value is m0; This represents the sum of spatial derivatives. Due to inaccurate initial velocities, there are differences between field-observed early arrival wave data and simulated early arrival wave data, such as... Figure 10(c) A comparison of the waveforms of early arrival wave seismic data observed in the field at a distance of 0.1 km with those of simulated early arrival wave seismic data is shown.

[0104] Step 3, constructing the accompanying seismic source:

[0105] The travel time difference between field-observed early-arrival wave seismic data and simulated early-arrival wave seismic data is calculated and denoted as Δτ(x,t). Based on this, the source of the accompanying wave field, i.e., the accompanying source, is constructed, and its representation is as follows:

[0106]

[0107] In the formula, Let represent the first and second time derivatives of the early arrival wave seismic data observed in the field after time shifting, respectively. The early arrival wave seismic data observed in the field after time shifting is consistent with the simulated early arrival wave seismic data in time; therefore, the early arrival wave seismic data observed in the field can be replaced by the simulated early arrival wave seismic data. Equation (20) can be further expressed as:

[0108]

[0109] In the formula, These represent the first and second time derivatives of the simulated early arrival wave seismic data, respectively. The main contributing terms are: and For a specific number, it represents a secondary contribution.

[0110] Step 4: Numerical calculation of the adjoint wave field:

[0111] The accompanying source of equation (3) is used as the source to replace the seismic wavelet f(t) in equation (19), and then the time domain wave equation (20) is solved numerically again to obtain the accompanying wave field q(x,t);

[0112] Step 5: Obtain the velocity update amount for each iteration:

[0113] The gradient operator is calculated based on the source wavefield u(x,t) and the adjoint wavefield q(x,t).

[0114]

[0115] In the formula, s and g represent the source point and receiver point of the early arrival wave seismic data generated in the field observation, respectively. First, a trial step size is given, and then the appropriate step size is obtained by using the one-dimensional line search method with Armijo conditions. The gradient operator is applied to calculate the velocity parameter model update amount for this iteration, ensuring that the velocity parameter model update amount for each iteration is between 20m / s and 100m / s. The updated velocity parameter model m is calculated using equation (5):

[0116]

[0117] In the formula, m0 represents the initial velocity parameter model, and λ represents the step size. This indicates the amount of time the velocity parameter model has been updated in this iteration;

[0118] Step Six: Construct an optimized error function and determine the convergence result:

[0119] An optimized error function is constructed based on the main contributions of the accompanying epicenter:

[0120]

[0121] The optimized error function not only considers the main contribution of the accompanying source but also avoids the problem of waveform errors easily getting caught in cycle skipping. Error functions for different iteration numbers are as follows: Figure 11 As shown, the part to the left of the white dashed line represents the error function of the first iteration, and the part to the right of the white dashed line represents the error function of the 17th iteration. As the inversion iterations proceed, the value of the error function gradually decreases, indicating successful convergence. Figure 12 As shown.

[0122] Returning to step two, the updated velocity parameter model m is substituted into equation (1) for numerical simulation. The simulated early arrival wave seismic data is recalculated. Equation (6) is used as the judgment condition to determine whether the inversion has converged. Then, the option is to continue iterative updates or terminate the iteration and output the inversion results.

[0123] Inversion results are shown Figure 13 The inversion of the early-arrival wave equation with optimized error function can effectively update the background velocity parameter model in the shallow and intermediate layers. A comparison of the pseudo-well curves from the inversion results is shown below. Figure 14 (a) Pseudo-well curve at a distance of 7.5 km, (b) Pseudo-well curve at a distance of 15 km, where the black solid line represents the true velocity and the black dashed line represents the inversion velocity. It can be seen that the background velocity in the shallow and middle layers is updated after the error is successfully converged.

[0124] Returning to step two, substitute the final inverted velocity parameter model m into equation (1) for numerical simulation, and recalculate the simulated early arrival wave seismic data, such as... Figure 15 As shown in Figure a. The waveforms of the early arrival wave observed in the field at a distance of 0.1 km are compared with those of the simulated early arrival wave, as shown in Figure a. Figure 15 As shown in (b), compared to Figure 10 c. The better agreement between the two indicates that the inversion of the early-arrival wave equation with optimized error function can accurately recover the background velocity parameter model of the middle and shallow layers.

[0125] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0126] Except for the technical features described in the specification, all other technologies are known to those skilled in the art.

Claims

1. A method for travel time inversion of the early-arrival wave equation with optimized error function, characterized in that, The early-arrival wave equation travel-time inversion method for optimizing the error function includes: Step 1: Input early arrival wave seismic data, seismic wavelet, and initial velocity model from field observations; Step 2: Numerically solve the wave equation to obtain the source wavefield and simulated earthquake data; Step 3: Calculate the travel time difference between observed seismic data and simulated seismic data and construct the associated seismic source; Step 4: Numerically solve the wave equation to obtain the adjoint wave field; Step 5: Calculate the gradient and obtain the velocity update amount; Step 6: Construct an optimized error function based on the travel time difference and seismic data; Step 7: When the inversion results converge, output the inversion results.

2. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 1, characterized in that, In step 1, input the early arrival wave seismic data d from field observations. obs (x,t), seismic wavelet f(t), and initial velocity parameter model m0; the initial velocity parameter model m0 is known, and its velocity value gradually increases from shallow to deep layers.

3. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 1, characterized in that, In step 2, based on the early arrival wave seismic data observed in the field, d obs (x,t) is used to extract the seismic wavelet f(t). Using the seismic wavelet f(t) as the source, the time-domain wave equation is numerically solved to obtain the source wavefield u(x,t) at all time points. Data at the receiving point is taken from the forward modeling simulation wavefield u(x,t), i.e., the simulated early arrival wave seismic data d. cal (x,t).

4. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 3, characterized in that, In step 2, the time-domain wave equation is expressed as follows: In the formula, u(x,t) represents the forward modeling wave field, x represents the spatial position vector, m represents the velocity parameter model, and the initial value is m0; It represents the sum of spatial derivatives.

5. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 1, characterized in that, In step 3, the travel time difference between the early arrival wave seismic data observed in the field and the simulated early arrival wave seismic data is calculated and denoted as Δτ(x,t). Based on this, the source of the accompanying wave field is constructed, i.e., the accompanying source.

6. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 5, characterized in that, In step 3, the accompanying earthquake source η(x,t) is represented as follows: In the formula, and Let represent the first and second time derivatives of the early arrival wave seismic data observed in the field after time shifting, respectively. The early arrival wave seismic data observed in the field after time shifting is consistent with the simulated early arrival wave seismic data in time. Therefore, the early arrival wave seismic data observed in the field can be replaced by the simulated early arrival wave seismic data. Equation (2) can be further expressed as: In the formula, and These represent the first and second time derivatives of the simulated early arrival wave seismic data, respectively. The main contributing terms are: and For a specific number, it represents a secondary contribution.

7. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 6, characterized in that, In step 4, the accompanying source is used as the source to replace the seismic wavelet in step 2, and then numerical forward modeling is performed again to obtain the accompanying wave field; that is, the accompanying source of equation (3) is used as the source to replace the seismic wavelet f(t) in equation (1), and then the time domain wave equation (2) is numerically solved again to obtain the accompanying wave field q(x,t).

8. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 1, characterized in that, In step 5, the gradient operator is calculated based on the source wave field and the associated wave field, and an appropriate step size is obtained to apply to the gradient operator to calculate the velocity parameter model update amount for this iteration.

9. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 8, characterized in that, In step 5, the velocity update amount for each iteration is: The gradient operator is calculated based on the source wavefield u(x,t) and the adjoint wavefield q(x,t). In the formula, s and g represent the source point and receiver point of the early arrival wave seismic data generated in the field observation, respectively; first, a trial step size is given, and then the appropriate step size is obtained by the one-dimensional line search method with Armijo conditions. The gradient operator is applied to calculate the velocity parameter model update amount of this iteration, and the updated velocity parameter model m is calculated using equation (5): In the formula, m0 represents the initial velocity parameter model, and λ represents the step size. This indicates the amount of time the velocity parameter model has been updated in this iteration.

10. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 1, characterized in that, In step 6, an optimized error function is constructed based on the terms of the main contributions from the accompanying earthquake source: Where Δτ(x,t) represents the travel time difference between the early arrival wave seismic data observed in the field and the simulated early arrival wave seismic data. denoted as the first time derivative of the simulated early arrival wave seismic data, s and g represent the source point and receiver point that generated the early arrival wave seismic data from field observations, respectively, and t is the time it takes for the early arrival wave to reach the receiver point.

11. The method for travel time inversion of the early-arrival wave equation with optimized error function according to claim 1, characterized in that, In step 7, if the inversion result does not converge, return to step 2 to continue iterative updates.

12. An early arrival wave equation travel-time inversion system with optimized error function, characterized in that the early arrival wave equation travel-time inversion system with optimized error function uses the early arrival wave equation travel-time inversion method with optimized error function as described in any one of claims 1-11 to perform seismic wave propagation velocity model inversion.

13. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored computer program; wherein, when the computer program is executed, it controls the device on which the computer-readable storage medium is located to perform a travel time inversion method for early arrival wave equations with optimized error functions as described in any one of claims 1-11.

Images