Local scale traveltime inversion method based on mixed source data

Abstract

The application belongs to the field of seismic data processing, and is a local scale travel time inversion method based on mixed source data. The low wave number component of the underground velocity model is constructed by using the travel time information of the seismic data, and the calculation efficiency of the travel time inversion is greatly improved by combining the mixed source coding strategy. The limitation that the conventional travel time inversion method cannot use mixed source data is broken. The local scale transformation is performed on the seismic signal to construct a local scale travel time inversion objective function. Then, the mixed source coding strategy is introduced to improve the calculation efficiency of the travel time inversion, avoid the crosstalk noise between shot and shot data, increase and decrease the surface coverage times, and improve the accuracy of the seismic wave velocity modeling in the underground space. The gradient operator corresponding to the local scale travel time inversion objective function of the mixed source data is derived, and the velocity parameter model is updated and iterated by using an optimization algorithm. The reliability of the method is verified by testing the Marmousi velocity model data.

Description

Technical Field

[0001] This invention belongs to the field of earthquake data processing, specifically relating to a local-scale travel time inversion method based on mixed source data. More specifically, it is a method that uses a local-scale travel time inversion method to perform velocity modeling on mixed source data, thereby improving computational efficiency and mitigating period jumps in the full waveform inversion process. Background Technology

[0002] The full waveform inversion method was first proposed by Lailly (1983) and Tarantola (1984) in the 1980s. This method fully utilizes the kinematic and dynamic information of seismic data to construct the objective function of the inversion algorithm. Through continuous optimization and iteration, the algorithm ultimately finds the physical property parameters that minimize the objective function value. Due to the enormous computational cost of full waveform inversion, the adjoint state method is used to calculate the gradient to avoid directly calculating the Fréchet derivative, greatly accelerating the industrialization of full waveform inversion. Subsequently, Pratt et al. proposed a frequency domain full waveform inversion method, pointing out that only a few discrete seismic data points from low to high frequencies are needed to obtain high-precision inversion results (Pratt et al., 1998; Pratt, 1999), which also greatly improves the computational efficiency of the forward and reverse propagation wavefields in the full waveform inversion process. It wasn't until the 21st century, with the significant improvement in computer computing performance, that the full waveform inversion data processing workflow officially entered the industrialization stage.

[0003] Considering the strong nonlinearity of the objective function in full waveform inversion, although the full waveform inversion method has been developed for nearly 40 years, the periodic jump problem remains the biggest obstacle to practical application. When the waveforms of observed and simulated data differ in phase by more than half a period, a periodic jump phenomenon will occur during the full waveform inversion process. To address this, Luo and Schuster (1991) proposed a wave equation travel-time inversion method based on the gradient calculation using the adjoint state method in full waveform inversion, and combined the travel-time inversion method with the full waveform inversion method to obtain high-precision inversion results for inter-well seismic data. Fichtner et al. (2008) developed a time-frequency domain full waveform inversion method, constructing envelope and phase objective functions in the time and frequency domains respectively, and obtained good inversion results in actual data experiments of natural earthquake full waveform inversion. Warner and Guasch proposed an adaptive full-waveform inversion method, which uses Wiener filtering to reshape simulated data to more closely approximate the phase changes of observed data. Numerical experiments show that this method can effectively overcome period jumps (Guasch et al., 2014) and correct seismic waveform travel time differences even when seismic data lacks low frequencies, source wavelet errors occur, and the initial velocity model differs significantly from the true velocity model. This prevents the phase difference between observed and simulated data waveforms from exceeding half a period. Zhu and Fomel (2016) used an adaptive matched filter to measure the phase difference between observed and simulated data and pointed out that using waveform-corrected seismic data to construct the objective function has the advantage of smoothing the second-order functional. To avoid excessive phase differences between observed and simulated data waveforms, the most direct method is to use waveform correction techniques to make the waveform change characteristics of observed data closer to those of simulated data. Hu et al. (2018) constructed a time-frequency domain phase correction objective function and pointed out that the time-frequency domain phase-corrected observed data can be used to gradually guide the simulated data out of the local minima in the full-waveform inversion objective function. Guasch et al. (2019) tested the adaptive full-waveform inversion method on real data and pointed out that the method has good performance in practical applications. Hu Yong et al. (2020) used a local-scale travel-time inversion method based on the wave equation to transform the observed data at a local scale, fully considering the travel-time error characteristics within the local scale range, and making good use of reflected wave information to improve the accuracy of travel-time inversion. Therefore, based on previous research, this paper fully considers the travel-time information of seismic waves and uses mixed source data to construct macroscopic information of underground velocity structure, providing a better initial velocity model for full-waveform inversion. Summary of the Invention

[0004] The technical problem to be solved by this invention is to overcome the difficulty of travel time inversion for mixed seismic sources by using local scale transformation. Furthermore, by utilizing the local scale travel time information of seismic data, a local scale travel time inversion method based on mixed seismic source data is proposed to invert the low wavenumber component of the subsurface velocity model, alleviate the period jump problem of full waveform inversion, and improve the computational efficiency of full waveform inversion.

[0005] This invention provides a local scale travel-time inversion method based on hybrid source data, comprising:

[0006] Step 1. Use the earthquake data processing software Geoeast to denoise the earthquake data to obtain mixed source earthquake data, and use the processed mixed source earthquake data as the data input for the local scale travel time inversion method based on the wave equation;

[0007] Step 2. Construct an initial model of seismic wave velocity parameters as the model input for forward modeling using the local scale travel-time inversion method based on hybrid source data.

[0008] Step 3. Define the observation system and set the source wavelet; perform forward modeling using the initial model of seismic wave velocity parameters to obtain the propagating wave field, where the wave equation for a two-dimensional acoustic constant-density isotropic medium is expressed as:

[0009]

[0010] Where v is the velocity of sound; u is the sound wave field; x and z represent the horizontal and vertical directions, respectively; t represents time. Based on the sound wave equation and combined with the finite difference algorithm, the propagation process of seismic waves underground is simulated.

[0011] Step 4. Acquisition of Hybrid Source Data: The seismic sources are encoded using a random source coding strategy. The coding includes: spatial location, excitation delay, etc. The wavelet matrix based on the dynamic hybrid source is represented as:

[0012] S p =D y ·T·L·f

[0013] Where D y It is a random dynamic coding function; T is a random excitation delay coding function; L is a random location coding function; f represents the source wavelet, which combines multiple sources into a super source, and the dynamic coding strategy D... y This means that the source coding method will be updated every ten velocity inversion iterations;

[0014] Step 5. Obtain the local scale waveform information of the observed and simulated data using local scale transformation, where the local scale transformation of the observed and simulated data is as follows:

[0015] U(k,t)=u(t)h(tk)

[0016] u(t)=∫U(k,t)h(tk)dk

[0017] Where u(t) represents the mixed source seismic data, U(k,t) represents the seismic waveform data after local scaling transformation, h(t) represents the Gaussian window function, and k is the shift of the Gaussian window function along the time axis. The original mixed source seismic data is obtained by using the inverse transformation formula, which requires that ||h(t)||2=1;

[0018] Step 6. Calculate the local travel time difference between the observed data and the simulated data: Use the cross-correlation algorithm to calculate the travel time difference between the two waveforms. The corresponding cross-correlation algorithm calculates the travel time difference between the local-scale observed data and the simulated data as follows:

[0019]

[0020] Where U(k,t) and D(k,t) represent the observed data and simulated data obtained by local scaling transformation, respectively, and the travel time difference Δτ between the two waveforms is obtained by taking the maximum value of the cross-correlation result of the two waveforms away from the center point position:

[0021]

[0022] Where T0 is the maximum travel time difference between the pre-set observed data and the simulated data, and Δτ represents the travel time difference. If Δτ = 0, it means that there is no travel time difference between the observed data and the simulated data, and the waveform matching is successful.

[0023] Step 7. Constructing a local-scale travel-time inversion objective function based on hybrid source data: Since travel-time information has a better linear correspondence with the subsurface medium, the travel-time information is used to construct the objective function for subsurface velocity inversion. The specific objective function is expressed as follows:

[0024]

[0025] Where W represents the weighting factor, Δτ∈[-T0,T0] represents the local scale travel time information of the observed data and the simulated data, ns represents the number of seismic sources, and nr represents the number of detectors;

[0026] Step 8. Calculate the partial derivatives of the objective function with respect to the velocity parameters:

[0027]

[0028] Using the partial derivative of the cross-correlation function (C(k,τ)) with respect to τ as a bridge, the relationship between travel time difference information and local scale wavefield is established. When τ = Δτ, the first derivative of the cross-correlation function with respect to time is 0, and the following relationship holds:

[0029]

[0030] By establishing the relationship between travel time difference and wave field using the above formula, the partial derivative of travel time difference with respect to velocity can be expressed as an implicit function as follows:

[0031]

[0032] in Finally, the partial derivative of the objective function with respect to the velocity parameter can be expressed as:

[0033]

[0034] Step 9. Define the accompanying source for the local scale travel-time inversion method of the hybrid source data, and backpropagate the accompanying source to the model space to obtain the backpropagated wavefield, which is used to calculate the gradient. The corresponding accompanying source is:

[0035]

[0036] Step 10. Calculate the update direction of the model perturbation using the L-BFGS optimization algorithm. The iterative formula is as follows:

[0037] m k+1 =m k -α k H k g k

[0038] Where m k Let α be the update amount of the model perturbation at step k. k H is the step size. k It is the inverse of the approximate Hessian matrix. This is the gradient update for model perturbation. In L-BFGS optimization updates, only a few vector pairs need to be saved for updating the Hessian matrix, and the update formula is as follows:

[0039] H k+1 =V k T H k V k +ρ k s k s k T

[0040]

[0041] s k =m k+1 -m k y k =g k+1 -g k

[0042] Where H k+1 It is based on the vector pair {s k ,y k} and H k Calculated; H k g k The product can be obtained through the gradient g. k With vector pair {s k ,y k The approximate Hessian matrix H is obtained by the inner product of a series of vectors and the sum of the vectors. k The following update formula must be met:

[0043]

[0044] Step 11. Determine if the termination condition is met. If it is, output the results of the local scale travel-time inversion method based on the hybrid source data. If the termination condition is not met, continue to use the current inversion results as the initial velocity parameter model for the next cycle until the termination condition is met.

[0045] This invention fully utilizes travel time information from seismic data to construct the low wavenumber component of the subsurface velocity model. Simultaneously, by incorporating a hybrid source coding strategy, it significantly improves the computational efficiency of travel time inversion and overcomes the limitation of conventional travel time inversion methods that cannot utilize hybrid source data. First, a local scale transformation is performed on the seismic signal to construct a local-scale travel time inversion objective function. Second, a hybrid source coding strategy is introduced to improve the computational efficiency of travel time inversion, avoid crosstalk noise between shot data, and increase or decrease the number of surface overlays to improve the accuracy of subsurface seismic wave velocity modeling. Finally, the gradient operator corresponding to the local-scale travel time inversion objective function of the hybrid source data is derived, and an optimization algorithm is used to update and iterate the velocity parameter model. Testing with Marmousi velocity model data verifies that this invention can obtain the low wavenumber component inversion results of the subsurface velocity model. Attached Figure Description

[0046] Figure 1 This is a flowchart of the local scale travel time inversion method based on hybrid source data according to the present invention;

[0047] Figure 2 It is hybrid source data after hybrid source coding;

[0048] Figure 3It is the Marmousi velocity model, (a) the actual velocity model; (b) the initial velocity model;

[0049] Figure 4 The results are conventional multi-scale full waveform inversion results based on hybrid source data: (a) low-frequency band inversion results (4-6Hz); (b) mid-frequency band inversion results (4-10Hz); (c) high-frequency band inversion results (4-15Hz).

[0050] Figure 5 The results are based on local-scale travel-time inversion of hybrid source data: (a) Local-scale travel-time inversion results (4-6Hz); (b) Low-frequency band inversion results (4-6Hz); (c) Mid-frequency band inversion results (4-10Hz); (d) High-frequency band inversion results (4-15Hz). Detailed Implementation

[0051] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0052] The steps of the local-scale travel-time inversion method based on hybrid source data are as follows:

[0053] Step 1. Use the seismic data processing software Geoeast to denoise the seismic data, and use the processed mixed source seismic data as the data input for the local scale travel-time inversion method based on the wave equation.

[0054] Step 2. Construct an initial model of seismic wave velocity parameters as the model input for forward modeling using the local scale travel-time inversion method based on hybrid source data.

[0055] Step 3. Define the observation system and set the source wavelet; perform forward modeling using the initial velocity model and store the propagating wave field. The wave equation for a two-dimensional acoustic constant-density isotropic medium can be expressed as:

[0056]

[0057] Where v is the velocity of sound; u is the sound wave field; x and z represent the horizontal and vertical directions, respectively; and t represents time. Based on the sound wave equation and combined with the finite difference algorithm, the propagation process of seismic waves underground can be simulated.

[0058] Step 4. Acquisition of Hybrid Source Data: The seismic sources are encoded using a random source coding strategy. The coding includes: spatial location and excitation delay. Therefore, the wavelet matrix based on dynamic hybrid sources can be expressed as:

[0059] S p=D y ·T·L·f

[0060] Where D y is a random dynamic coding function; T is a random excitation delay coding function; L is a random location coding function; f represents the source wavelet. Mixing multiple sources into a single super-source can significantly improve computational efficiency, and this coding strategy can effectively reduce the impact of crosstalk noise on the inversion results. Dynamic coding strategy D y This means that the source coding method is updated every ten velocity inversion iterations, which can effectively suppress the influence of crosstalk noise. In addition, by continuously changing the spatial location of the source excitation, after multiple iterations, the surface of each grid can be covered by the seismic source, which can effectively increase the source coverage density and improve the inversion accuracy.

[0061] Step 5. Obtain the local scale waveform information of the observed and simulated data using local scale transformation. The local scale transformation for the observed and simulated data is as follows:

[0062] U(k,t)=u(t)h(tk)

[0063] u(t)=∫U(k,t)h(tk)dk

[0064] Where u(t) represents the mixed-source seismic data, U(k,t) represents the seismic waveform data after local scaling, h(t) represents the Gaussian window function, and k is the shift of the Gaussian window function along the time axis. To ensure that the original mixed-source seismic data can be obtained using the inverse transformation formula after local scaling of the seismic data, ||h(t)||² = 1 must be satisfied. In the above formula, if the length of the Gaussian window is too long, it may not only cover the waveform information of the entire seismic event, but also include waveform changes from multiple seismic events, resulting in severe crosstalk noise. If a short Gaussian window is selected to perform local scaling of the seismic data, it can effectively alleviate crosstalk noise, but it may divide a seismic event into different parts, making it impossible to extract the travel time information of a complete waveform.

[0065] Step 6. Calculate the local travel time difference between the observed data and the simulated data: Use the cross-correlation algorithm to calculate the travel time difference between the two waveforms. The corresponding cross-correlation algorithm for calculating the travel time difference between the local-scale observed data and the simulated data can be expressed as:

[0066]

[0067] Where U(k,t) and D(k,t) represent the observed data and simulated data obtained from the local scaling transformation, respectively. Furthermore, the travel time difference (Δτ) between the two waveforms can be obtained by using the deviation of the maximum value of the cross-correlation result of the two waveforms from the center point position:

[0068]

[0069] Where T0 is the maximum travel time difference between the pre-set observed data and the simulated data, and Δτ represents the travel time difference. If Δτ = 0, it means that there is no travel time difference between the observed data and the simulated data, and the waveform matching is successful.

[0070] Step 7. Constructing a local-scale travel-time inversion objective function based on hybrid source data: Since travel-time information has a better linear correspondence with the subsurface medium, travel-time information is used to construct the objective function for subsurface velocity inversion. The specific objective function can be expressed as:

[0071]

[0072] Where W represents the weighting factor, Δτ∈[-T0,T0] represents the local scale travel time information of the observed data and the simulated data, ns represents the number of seismic sources, and nr represents the number of detectors;

[0073] Step 8. Calculate the partial derivatives of the objective function with respect to the model parameters:

[0074]

[0075] However, the partial derivatives of the objective function contain This information cannot be directly calculated using observed and simulated data. Therefore, the partial derivative of the cross-correlation function (C(k,τ)) with respect to τ is needed as a bridge to establish the relationship between travel time difference information and the local scale wavefield. When τ = Δτ, the first derivative of the cross-correlation function with respect to time is 0, and the following relationship holds:

[0076]

[0077] The relationship between travel time difference and wave field can be easily established using the above formula. The partial derivative of travel time difference with respect to velocity can then be expressed as an implicit function:

[0078]

[0079] in Finally, the partial derivative of the objective function with respect to the velocity parameter can be expressed as:

[0080]

[0081] Step 9. Define the accompanying source for the local scale travel-time inversion method of the hybrid source data, and backpropagate the accompanying source to the model space to obtain the backpropagated wavefield, which is used to calculate the gradient. The corresponding accompanying source is:

[0082]

[0083] Step 10. Calculate the update direction of the model perturbation using the L-BFGS optimization algorithm. The iterative formula is as follows:

[0084] m k+1 =m k -α k H k g k

[0085] Where m k Let α be the update amount of the model perturbation at step k. k H is the step size. k It is the inverse of the approximate Hessian matrix. This is the gradient update for model perturbation. In L-BFGS optimization updates, only a few vector pairs need to be saved for updating the Hessian matrix, and the update formula is as follows:

[0086] H k+1 =V k T H k V k +ρ k s k s k T

[0087]

[0088] s k =m k+1 -m k y k =g k+1 -g k

[0089] Where H k+1 It is based on the vector pair {s k ,y k} and H k Calculated; H k g k The product can be obtained through the gradient g. k With vector pair {s k ,y k It is obtained by the inner product of a series of vectors and the sum of the vectors. The inverse matrix H of the approximate Hessian matrix is... k The following update formula must be met:

[0090]

[0091] Step 11. Determine if the termination condition is met. If it is, output the results of the local scale travel-time inversion method based on hybrid source data. If the termination condition is not met, continue to use the current inversion results as the initial velocity parameter model for the next cycle until the termination condition is met.

[0092] Example:

[0093] Because the objective function of full-waveform inversion exhibits strong nonlinear characteristics, it often gets trapped in local minima during the search for the minimum using local optimization algorithms, resulting in the inability of the full-waveform inversion method to obtain high-precision inversion results for the subsurface velocity model. Furthermore, the full-waveform inversion process requires continuous finite-difference forward modeling, resulting in a huge computational burden and further limiting the industrialization of this method. To address the problem of high computational cost in full-waveform inversion, this invention utilizes a hybrid source coding strategy. However, since it is difficult to extract the first arrival travel time or first arrival waveform from the seismic records after hybrid source coding, conventional travel time inversion methods based on wave equations cannot obtain ideal inversion results. Therefore, this invention proposes a local-scale travel time inversion method based on hybrid source data, which fully utilizes the first arrival waveform information of the seismic data to reconstruct the low wavenumber components of the subsurface velocity model, alleviating the period jump problem of full-waveform inversion. The specific process is as follows: Figure 1 As shown. A hybrid seismic source was set up on the Earth's surface, comprising 13 seismic sources with uneven spatial locations. Figure 2 Each seismic source has a different excitation delay, and the source code is changed every 10 iterations of the inversion results to maximize the number of times the seismic source covers the Earth's surface. Seismic data is recorded for 4 seconds at 2-ms intervals. The seismic source uses a Ricker wavelet with a dominant frequency of 10Hz, and the L-BFGS optimization algorithm is used to calculate the descent direction and update the model perturbations iteratively.

[0094] This invention uses Marmousi model data for testing, firstly utilizing a true velocity model ( Figure 3 a) Observational data is obtained by combining forward modeling of the wave equation and regarded as seismic data used in actual production processes. According to the process described in the claim, the seismic data is first preprocessed to construct an initial velocity model (3b), define the observation system, extract the source wavelet information corresponding to the seismic data, and define the hybrid source coding method. Local waveform information of the hybrid source is obtained by using local scale decomposition. Combined with the cross-correlation algorithm, the objective function of the local scale travel time inversion method based on hybrid source data is established, and the partial derivative of the objective function with respect to the model perturbation is obtained to obtain the update amount of the model perturbation. Finally, the L-BFGS optimization algorithm is used to continuously update the model perturbation based on the initial model. This invention also compares in detail the conventional multi-scale full waveform inversion method based on hybrid source data ( Figure 4 ) and local scale travel-time inversion method based on hybrid source data ( Figure 5 ).

[0095] The model parameters are as follows:

[0096] Table 1: Test parameters of the time-frequency domain amplitude-phase joint multi-scale elastic wave full waveform inversion method based on total variation constraints

[0097]

[0098] Will Figure 4 a and Figure 5 A comparison shows that the local-scale travel-time inversion method based on hybrid source data has a more uniform velocity structure than the conventional full-waveform inversion method based on hybrid source data, especially on the left side of the Marmousi velocity model, where the velocity inversion results are significantly better than the conventional full-waveform inversion method. Therefore, the local-scale travel-time inversion results based on hybrid source data are closer to the actual model perturbation. As the frequency of seismic data gradually increases, the detailed information of the Marmousi velocity model is gradually recovered. Figure 4 c and Figure 5 A comparison of the two methods reveals that the local-scale travel-time inversion method based on hybrid source data, combined with the conventional full-waveform inversion method, yields better inversion results, especially on the left side of the Marmousi velocity model. Due to the better linear correspondence between the travel-time information of the seismic data and the subsurface velocity model, the dependence of the full-waveform inversion method on the initial velocity model is significantly reduced. Numerical experimental results demonstrate that the local-scale travel-time inversion method based on hybrid source data has certain advantages in low-wavenumber component velocity modeling.

Claims

1. A local-scale traveltime inversion method based on hybrid source data, characterized in that, include: Step 1. Use the earthquake data processing software Geoeast to denoise the earthquake data to obtain mixed source earthquake data, and use the processed mixed source earthquake data as the data input for the local scale travel time inversion method based on the wave equation; Step 2. Construct an initial model of seismic wave velocity parameters as the model input for forward modeling using the local scale travel-time inversion method based on hybrid source data; Step 3. Define the observation system and set the source wavelet; perform forward modeling using the initial model of seismic wave velocity parameters to obtain the propagating wave field, where the wave equation for a two-dimensional acoustic constant-density isotropic medium is expressed as: ， in The speed of sound waves; For sound wave field; and These represent the horizontal and vertical directions, respectively. The time is represented by the sound wave equation, and the process of seismic wave propagation underground is simulated by combining the finite difference algorithm. Step 4. Acquisition of Hybrid Source Data: The seismic sources are encoded using a random source coding strategy. The coding includes: spatial location, excitation delay, and wavelet matrix representation based on dynamic hybrid sources: ， in It is a random dynamic coding function; It is a randomly triggered delay coding function; It is a random position encoding function; Representing the source wavelet, multiple sources are mixed into a super source, using a random dynamic coding function. This means that the source coding method will be updated every ten velocity inversion iterations; Step 5. Obtain the local scale waveform information of the observed and simulated data using local scale transformation, where the local scale transformation of the observed and simulated data is as follows: ， ， in This indicates mixed-source earthquake data. This represents seismic waveform data after local scaling. Represents the Gaussian window function; Let be the shift of the Gaussian window function along the time axis. The original mixed-source seismic data is obtained using the inverse transform formula, which needs to satisfy... ; Step 6. Calculate the local travel time difference between the observed data and the simulated data: Use the cross-correlation algorithm to calculate the travel time difference between the two waveforms. The corresponding cross-correlation algorithm calculates the travel time difference between the local-scale observed data and the simulated data as follows: ， in and These represent the travel time difference between the observed data and the simulated data obtained through local scaling transformation, respectively. The value of the cross-correlation result of two waveforms is obtained by measuring the deviation of the maximum value from the center point. ， in It is the maximum travel time difference between the pre-set observed data and the simulated data. Indicates time difference, if This indicates that there is no time difference between the observed data and the simulated data, and the waveform matching is successful; Step 7. Constructing a local-scale travel-time inversion objective function based on hybrid source data: Since travel-time information has a better linear correspondence with the subsurface medium, the travel-time information is used to construct the objective function for subsurface velocity inversion. The specific objective function is expressed as follows: ， in Indicates the weighting factor. This represents the local-scale travel time information of observed and simulated data. Indicates the number of earthquake focal points. Indicates the number of detectors; Step 8. Calculate the partial derivatives of the objective function with respect to the velocity parameters: ， Using cross-correlation function ( ) right The partial derivative is used as a bridge to establish the relationship between travel time difference information and local scale wavefield. The partial derivative of the cross-correlation function with respect to time is... When the first derivative of the cross-correlation function is 0, the following relationship holds: ， By establishing the relationship between travel time difference and wave field using the above formula, the partial derivative of travel time difference with respect to velocity can be expressed as an implicit function as follows: ， in , , Finally, the partial derivative of the objective function with respect to the velocity parameter can be expressed as: ， Step 9. Define the accompanying source for the local scale travel-time inversion method of the hybrid source data, and backpropagate the accompanying source to the model space to obtain the backpropagated wavefield, which is used to calculate the gradient; the corresponding accompanying source is: ， Step 10. Calculate the update direction of the model perturbation using the L-BFGS optimization algorithm. The iterative formula is as follows: ， in Let be the update amount of the model perturbation at step k. Step size, It is the inverse of the approximate Hessian matrix. The gradient is used to update the model perturbation. In L-BFGS optimization updates, only a few vector pairs need to be saved for updating the Hessian matrix, and the update formula is as follows: ， ， ， ， ， in It is based on vector pairs and Calculated; The product can be obtained through the gradient. With vector pairs It is obtained by taking the inner product of a series of vectors and the sum of the vectors, where the inverse of the approximate Hessian matrix is ​​the matrix. The following update formula must be met: ， Step 11. Determine if the termination condition is met. If it is, output the results of the local scale travel-time inversion method based on the hybrid source data. If the termination condition is not met, continue to use the current inversion results as the initial velocity parameter model for the next cycle until the termination condition is met.

Images