Near-surface formation longitudinal wave velocity inversion method based on surface wave information

Abstract

The invention provides a surface wave information-based near-surface formation longitudinal wave velocity inversion method. The surface wave information-based near-surface formation longitudinal wave velocity inversion method comprises the following steps of: 1, extracting surface wave frequency dispersion data through a phase shift method; 2, according to the frequency dispersion curve, direct inversion is carried out to obtain the near-surface shear wave speed and thickness; step 3, calculating a W-D curve and a W / D-D curve; and 4, adding the W-D data and the W / D-D data into an objective function, and carrying out longitudinal wave velocity direct inversion. According to the near-surface formation longitudinal wave velocity inversion method based on the surface wave information, direct inversion of the near-surface formation longitudinal wave velocity can be achieved, information such as the transverse wave velocity, the formation thickness, the longitudinal wave velocity and the Poisson's ratio is obtained in the inversion process, and powerful research basis and support are provided for application of the surface wave information in oil exploration and engineering exploration.

Description

Technical Field

[0001] The invention relates to the technical field of earthquake source exploration, and in particular to a near-surface stratum P-wave velocity inversion method based on surface wave information. Background Art

[0002] For mineral exploration, P-wave velocity is often more important than S-wave velocity in applications such as static correction. However, the dispersion data currently used is very insensitive to P-wave velocity, and its influence on dispersion data is almost negligible compared with S-wave velocity and layer thickness. It is often given directly during S-wave velocity structure inversion, while it is almost impossible to obtain P-wave velocity directly from dispersion data.

[0003] Although P-wave velocity can be obtained through other methods such as tomographic inversion, micro-logging, small refraction, etc., its velocity structure is often not as fine as the S-wave velocity structure obtained by surface waves. Another method is to obtain the P-wave velocity by combining the Poisson's ratio with the S-wave velocity structure through independent data. Experiments have shown that the given Poisson's ratio in the inversion process will affect the effect of the S-wave velocity structure calculation, so scholars began to study the information in the dispersion information that may be sensitive to the Poisson's ratio.

[0004] A key breakthrough is that Pelekis et al. found that the detection depth parameter is sensitive to Poisson's ratio, and Pan et al. proposed an empirical linear formula between the detection depth coefficient and Poisson's ratio. These works did not directly obtain the P-wave velocity, but found the key information hidden in the dispersion data that is sensitive to the P-wave velocity. Socco et al. made a specific analysis of the relationship between the surface wave wavelength (W) and the detection depth (D), explained the stable relationship between the WD relationship and the Poisson's ratio, proposed the concept of the apparent Poisson's ratio, and inferred the P-wave velocity from the calculated Poisson's ratio. Anjom et al. proposed a specific process for reconstructing the near-surface velocity structure based on the apparent Poisson's ratio, compared the uncertainty, and obtained the P-wave velocity model at 10-70m underground through this method.

[0005] However, the P-wave velocity model obtained by the apparent Poisson's ratio method still has certain errors, and the apparent Poisson's ratio is very insensitive to shallow layers, so it is impossible to obtain a more accurate near-surface P-wave velocity during inversion.

[0006] Because the dispersion curve extracted by surface waves is very sensitive to shear wave velocity, at present, what is obtained by surface waves is mainly shear wave information near the surface, and it is impossible to invert the longitudinal wave velocity through dispersion data. In oil exploration, the application range of longitudinal wave velocity is wider than that of transverse wave. Although the longitudinal wave velocity can be obtained through some other methods such as tomographic inversion, micro-logging, small refraction, etc., its velocity structure is often not as fine as the transverse wave velocity structure obtained by surface waves. Another method is to obtain the longitudinal wave velocity by combining the Poisson's ratio with the shear wave velocity structure through independent data. Experiments have shown that the given Poisson's ratio in the inversion process will affect the effect of the shear wave velocity structure calculation.

[0007] In the Chinese patent application with application number: CN202111572999.3, a method for constructing a near-surface velocity model is involved, which belongs to the field of seismic exploration. The method includes: obtaining linear surface wave data in the three-dimensional surface wave data of the seismic gun in the target area, processing to obtain dispersion curve data, establishing a variable layer initial model of the target area, setting the variable layer initial model as the target model, calculating the predicted dispersion curve data, based on the difference between the dispersion curve data and the predicted dispersion curve data and the correction value of the velocity of each stratum at each detection point in the target model, optimizing the target model, processing the variable layer optimization model, obtaining a three-dimensional variable layer velocity model of the target area, establishing a work area elevation surface model of the target area, fusing the three-dimensional variable layer velocity model with the work area elevation surface model, and obtaining a near-surface velocity model of the target area. The use of this application helps to solve the problem of low accuracy in depicting near-surface details when constructing a velocity model using longitudinal wave seismic data.

[0008] In the Chinese patent application with application number: CN202211090913.8, a LSTM surface wave inversion method, device and medium are involved. The method includes determining the fuzzy formation parameter interval of the work area to be inverted; randomly generating different models according to the fuzzy formation parameter interval, calculating the theoretical surface wave dispersion curve of each model using the generalized reflection-transmission coefficient method, and constructing training sample data pairs; preprocessing the training sample data pairs; constructing an LSTM network based on the FHLV loss function, training the LSTM network based on the preprocessed training sample data pairs and saving the trained model; using unsupervised learning to automatically pick up the dispersion curve of the actual dispersion imaging data, and using the trained model to predict the dispersion curve after the unified dimension size to obtain the near-surface shear wave velocity model. The invention improves the processing efficiency and inversion accuracy, and is suitable for large-scale data processing.

[0009] In the Chinese patent application with application number: CN201310520027.9, a method of combined static correction of P-wave and S-wave based on inversion of surface S-wave velocity is involved. This method effectively improves the accuracy of direct inversion of S-wave velocity through high-precision wave field interpolation, combined application of P-wave and converted wave data surface wave information, and weighted superposition, and effectively solves the problem of inaccurate static correction in conventional methods by indirectly obtaining S-wave static correction, thereby effectively improving the imaging level of multi-wave data. After application tests in the processing of three-component seismic data in Luojia area of ​​Shengli Oilfield, good processing results were achieved, proving that the combined static correction method of P-wave and S-wave based on inversion of surface S-wave velocity is a feasible method technology.

[0010] The above existing technologies are greatly different from the present invention and fail to solve the technical problem we want to solve. Therefore, we have invented a new method for inverting the P-wave velocity of near-surface strata based on surface wave information. Summary of the invention

[0011] The object of the present invention is to provide a near-surface stratum P-wave velocity inversion method based on surface wave information, which directly inverts the P-wave velocity using dispersion data extracted from surface wave information.

[0012] The object of the present invention can be achieved by the following technical measures: a near-surface stratum P-wave velocity inversion method based on surface wave information, the near-surface stratum P-wave velocity inversion method based on surface wave information comprises:

[0013] Step 1, extracting surface wave dispersion data by phase shift method;

[0014] Step 2, directly invert the near-surface shear wave velocity and thickness based on the dispersion curve;

[0015] Step 3, calculate the WD curve and the W / DD curve;

[0016] Step 4: Add both WD data and W / DD data to the objective function to perform direct inversion of P-wave velocity.

[0017] The purpose of the present invention can also be achieved by the following technical measures:

[0018] The near-surface stratum P-wave velocity inversion method based on surface wave information further includes, before step 1, acquiring surface wave data from actual seismic acquisition data.

[0019] In step 1, assuming that the simulated moving source surface wave data is d(x, t), its spectrum can be obtained through Fourier transform:

[0020]

[0021] Where A(x,ω) is the amplitude spectrum, P(x,ω)=exp(-iωx / c R ) is the phase spectrum, x is the detector position, c R (ω) is the actual phase velocity.

[0022] In step 1, since the phase spectrum contains dispersion information, it is necessary to first eliminate the influence of the amplitude spectrum and obtain the normalized spectrum data:

[0023]

[0024] In the formula, Normalized spectrum data, that is, phase spectrum data, is obtained by performing multiple phase velocity values ​​c corresponding to a certain frequency. T Substituting and summing, we can get the phase velocity spectrum corresponding to this frequency:

[0025]

[0026] Where V(c T ,ω) is the phase velocity spectrum, N is the total number of detectors, n is the detector number, x n is the distance between the nth detector and the shot point;

[0027] From formula (3), we can see that only when c T Determined as the actual phase velocity c R When , the phase velocity spectrum corresponding to this frequency can reach its maximum value, showing the strongest superposition energy on the dispersion image; repeating the above calculation for all frequencies can perform dispersion extraction on the moving source surface wave seismic records.

[0028] In step 2, the near-surface shear wave velocity and thickness are directly inverted based on the dispersion curve, and a certain artificially set theoretical near-surface shear wave velocity profile is sought. The dispersion data is obtained by forward calculation through the shear wave velocity, layer thickness and other information of each layer. If the curve and the dispersion data extracted from the measured data achieve the best fit, the profile is the shear wave velocity and thickness information.

[0029] In step 2, nonlinear inversion is used to perform a global search to obtain the optimal model parameters. To improve computational efficiency, the objective function definition method used by the fast inversion algorithm is adopted:

[0030]

[0031] Where N is the number of frequency points of the extracted dispersion data, m is the total number of models, and f i obs is the frequency, Vr i obs is the phase velocity corresponding to the extracted dispersion curve, Vp j , Vs j,ρ,h j is the longitudinal wave velocity, shear wave velocity, density and thickness parameters corresponding to the jth model, w i is the weight coefficient corresponding to different frequencies, Ψ j is the objective function value corresponding to the jth model; F vr This is the forward operator obtained by the Haskell-Thomson transfer matrix method.

[0032] In step 2, a genetic algorithm is used as an inversion method. The steps of using the genetic algorithm to invert the shear wave velocity and thickness of each layer in the target area include:

[0033] ①First, the model parameter variation range is given, and the combination of shear wave velocity and layer thickness is determined and encoded within this range;

[0034] ② Substitute the combined parameters into the objective function (4) and determine their fitness, that is, their ability to survive in the entire population, according to the size of each combined objective function. According to the objective function setting method, the smaller the objective function value, the higher the fitness;

[0035] ③ Determine the selection probability, crossover probability and mutation probability, and perform selection, crossover and mutation operations on the entire population to produce a new generation of population, so that the entire population can achieve survival of the fittest;

[0036] ④ Repeat the above operations for the new generation population until the optimal solution that meets the conditions is found.

[0037] In step 3, the average shear wave velocity Vsz at a certain depth is calculated based on the acquired shear wave velocity profile:

[0038]

[0039] Among them, n is the depth number, h i With V i is the thickness and shear wave velocity corresponding to the depth number.

[0040] In step 3, the fundamental dispersion curve is displayed through the relationship between wavelength and phase velocity, and the depth D and wavelength W corresponding to the same velocity in the shear wave average velocity curve and the dispersion curve are found to obtain the WD relationship curve of the model; for any point on the curve, the ratio between the ordinate W and the abscissa D is calculated to form a new W / DD curve.

[0041] In step 4, both WD data and W / DD data are added to the objective function:

[0042]

[0043] For problems of different orders of magnitude, normalization is performed when calculating the objective function; where N is the number of depth points, Wri obs With Wr i cal are the ratios of the wavelength and depth corresponding to the observed and calculated i-th depth point, respectively; max(W i obs ) and max(Wr i obs ) are the maximum values ​​of the observed wavelength and the ratio of wavelength to depth respectively; the objective function is calculated and iterated by genetic algorithm, and the optimal solution is the required longitudinal wave velocity.

[0044] In step 4, the genetic algorithm inversion of the P-wave velocity of each layer in the target area includes the following steps:

[0045] ①First, the model parameter variation range is given, and the P-wave velocity of each layer is randomly determined within this range and combined to generate initial parameters, which are binary-coded to generate the initial population;

[0046] ② Substitute the combined parameters into the objective function (6) and determine their fitness, that is, their ability to survive in the entire population, according to the size of each combined objective function. According to the objective function setting method, the smaller the objective function value, the higher the fitness;

[0047] ③ Determine the selection probability, crossover probability and mutation probability, and perform selection, crossover and mutation operations on the entire population according to the size of fitness to produce a new generation of population, so that the entire population can achieve survival of the fittest;

[0048] ④ Decode the new generation population and calculate the fitness of the new population, that is, the size of the objective function. If the objective function is less than the set value, that is, it is sufficiently converged, the combination of longitudinal wave velocity with the smallest objective function is selected as the optimal solution, or if the iteration exceeds the set maximum number of iterations, the loop ends;

[0049] ⑤If the objective function is greater than the set value and the number of iterations is insufficient, repeat the above operations until the optimal solution that meets the conditions is found.

[0050] The method for inverting the P-wave velocity of near-surface strata based on surface wave information in the present invention solves the problem that the P-wave velocity model obtained by the apparent Poisson's ratio method still has certain errors. On the basis of this method, combined with the Monte Carlo inversion method, a direct P-wave velocity inversion method combining the WD relationship with the W / DD relationship is proposed, which can directly invert the P-wave velocity through surface wave information. The innovation of the present invention lies in:

[0051] (1) The longitudinal wave velocity can be obtained only through the surface wave dispersion information, without the need for other information such as the first arrival wave. The surface waves that need to be removed can be effectively turned into treasures, and the near-surface can be finely modeled through surface waves, providing a reference for static correction, etc.;

[0052] (2) By changing the objective function, the P-wave velocities of all target layers from the surface to the deep layer can be accurately inverted simultaneously;

[0053] (3) In the process of calculating the longitudinal wave velocity, the shear wave velocity and thickness are inverted at the same time, and other information such as Poisson's ratio can be obtained.

[0054] Compared with the existing technology, the near-surface P-wave velocity inversion method based on surface wave information has the following advantages:

[0055] Through the calculation of this method, direct inversion of the P-wave velocity of near-surface strata can be achieved. During the inversion process, information such as S-wave velocity, stratum thickness, P-wave velocity, Poisson's ratio, etc. can be obtained simultaneously, providing a strong research foundation and support for the application of surface wave information in petroleum exploration and engineering exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 It is a flow chart of a specific embodiment of the method for inverting the near-surface stratum compressional wave velocity based on surface wave information of the present invention;

[0057] Figure 2 It is a flow chart of realizing the genetic algorithm of the present invention;

[0058] Figure 3 It is a schematic diagram of model parameters, average shear wave velocity and dispersion curve in a specific embodiment of the present invention;

[0059] Figure 4 In a specific embodiment of the present invention Figure 3 WD relationship curves of the three models;

[0060] Figure 5 In a specific embodiment of the present invention Figure 3 Schematic diagram of the WD curve fitting of the three models;

[0061] Figure 6 In a specific embodiment of the present invention Figure 3 Schematic diagram of the W / DD curve fitting of the three models;

[0062] Figure 7 In a specific embodiment of the present invention Figure 3 Schematic diagram of the P-wave velocity inversion results of three models;

[0063] Figure 8 A schematic diagram of actual surface roll raw data in a specific embodiment of the present invention;

[0064] Fig. 9 It is a schematic diagram of dispersion data extracted by a phase shift method in a specific embodiment of the present invention;

[0065] Fig.10 It is a schematic diagram of the dispersion curve fitting when performing inversion by genetic algorithm in a specific embodiment of the present invention;

[0066] Fig.11 is a schematic diagram of the inversion results of shear wave velocity and thickness in a specific embodiment of the present invention;

[0067] Fig.12 A schematic diagram of WD data calculation in a specific embodiment of the present invention;

[0068] Fig.13 It is a schematic diagram of the WD data corresponding to the P-wave velocity structure obtained by performing P-wave velocity inversion by genetic algorithm in a specific embodiment of the present invention and the fitting of the WD data with the actual WD data;

[0069] Fig.14 It is a schematic diagram of the fitting of W / DD data corresponding to the P-wave velocity structure obtained by performing P-wave velocity inversion by genetic algorithm and actual W / DD data in a specific embodiment of the present invention;

[0070] Fig.15 It is a schematic diagram of the longitudinal wave velocity inversion result in a specific embodiment of the present invention. DETAILED DESCRIPTION

[0071] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

[0072] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprising" and / or "including" are used in this specification, it indicates the presence of features, steps, operations and / or combinations thereof.

[0073] The surface wave information-based near-surface stratum P-wave velocity inversion method of the present invention extracts dispersion information and inverts to obtain shear wave velocity and stratum thickness information through surface wave information in seismic records. According to the sensitivity of Poisson's ratio to the wavelength (W)-depth (D) relationship, the WD curve corresponding to the dispersion curve is calculated, and the corresponding W / DD curve is calculated at the same time. The WD information and the W / DD information are combined to form an objective function. Through the Monte Carlo full-space search algorithm, in a nonlinear inversion method, the P-wave velocity profile corresponding to the minimum value of the objective function is obtained, which is the required P-wave function information.

[0074] In the above-mentioned direct inversion method of near-surface P-wave velocity, the required initial information is the surface wave information in the seismic record.

[0075] like Figure 1 As shown, Figure 1 The flowchart of the method for inverting the P-wave velocity of near-surface strata based on surface wave information of the present invention is as follows. The method for inverting the P-wave velocity of near-surface strata based on surface wave information comprises:

[0076] Step 1, obtaining seismic surface wave data;

[0077] Step 2, extracting dispersion curve using phase shift method;

[0078] Step 3, inversion to obtain shear wave velocity and thickness;

[0079] Step 4, calculating the wavelength W-detection depth D curve and the W / DD curve;

[0080] Step 5, construct the WD data and W / DD data objective function;

[0081] Step 6: Perform P-wave velocity inversion.

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

[0083] Example 1

[0084] In this embodiment 1, the near-surface P-wave velocity inversion method based on surface wave information of the present invention includes the following contents:

[0085] 1) Extract surface wave dispersion data by phase shift method;

[0086] Assuming that the simulated moving source surface wave data is d(x, t), its spectrum can be obtained through Fourier transform:

[0087]

[0088] Where A(x,ω) is the amplitude spectrum, P(x,ω)=exp(-iωx / c R ) is the phase spectrum, x is the detector position, c R (ω) is the actual phase velocity. Since the phase spectrum contains dispersion information, it is necessary to first eliminate the influence of the amplitude spectrum and obtain the normalized spectrum data:

[0089]

[0090] By performing multiple phase velocity values ​​c corresponding to a certain frequency T Substituting and summing, we can get the phase velocity spectrum corresponding to this frequency:

[0091]

[0092] Where N is the total number of detectors, and n is the detector number. From (3), we can see that only when c T Determined as the actual phase velocity c R When the phase velocity spectrum corresponding to this frequency can reach the maximum value, the strongest superposition energy is displayed on the dispersion image. Repeat the above calculation for all frequencies to extract the dispersion of the moving source surface wave seismic record.

[0093] 2) Directly invert the near-surface shear wave velocity and thickness based on the dispersion curve, and seek a certain artificially set theoretical near-surface shear wave velocity profile. The dispersion data of the shear wave velocity, layer thickness and other information of the profile obtained by forward calculation are best fitted with the actual dispersion data extracted by field reception. This profile is the shear wave velocity and thickness information.

[0094] Nonlinear inversion is used to perform a global search to obtain the best model parameters. In the later stage of inversion, the wavelength-depth curve (WD curve) needs to be calculated to achieve the inversion of the longitudinal wave velocity. If the Poisson's ratio is a fixed empirical value given in advance to fit the dispersion curve, this will affect the accurate extraction of the WD relationship, so the search range needs to be expanded. In the inversion, the Poisson's ratio is set as the target parameter and added to the calculation process. The search range is set to be wide (0.2–0.49), but this will result in a significant increase in the number of models. In order to improve the calculation efficiency, the objective function definition method used by the fast inversion algorithm is adopted:

[0095]

[0096] Where N is the number of frequency points of the extracted dispersion data, m is the total number of models, and f i obs is the frequency, Vr i obs is the phase velocity corresponding to the extracted dispersion curve, Vp j , Vs j ,ρ,h j is the longitudinal wave velocity, shear wave velocity, density and thickness parameters corresponding to the jth model, w i is the weight coefficient corresponding to different frequencies, Ψ j is the objective function value corresponding to the jth model. vr It is the forward operator obtained by the Haskell-Thomson transfer matrix method. The genetic algorithm is used as the inversion method. The genetic algorithm inversion of the target geological structure stratum shear wave velocity and thickness can be divided into the following steps:

[0097] ①First, the model parameter variation range is given, and the combination of shear wave velocity and layer thickness is determined and encoded within this range;

[0098] ② Substitute the combined parameters into the objective function (4) and determine their fitness, that is, their ability to survive in the entire population, according to the size of each combined objective function. According to the objective function setting method, the smaller the objective function value, the higher the fitness;

[0099] ③ Determine the selection probability, crossover probability and mutation probability, and perform selection, crossover and mutation operations on the entire population to produce a new generation of population, so that the entire population can achieve survival of the fittest;

[0100] ④ Repeat the above operations for the new generation population until the optimal solution that meets the conditions is found.

[0101] 3) Calculate the WD curve and W / DD curve. Based on the acquired shear wave velocity profile, first calculate the average shear wave velocity (Vsz) at a certain depth:

[0102]

[0103] Among them, n is the depth number, h i With V i is the thickness and shear wave velocity corresponding to the depth number. The fundamental dispersion curve is displayed by the relationship between wavelength and phase velocity. The WD relationship curve of the model can be obtained by finding the depth (D) and wavelength (W) corresponding to the same velocity in the shear wave average velocity curve and the dispersion curve. For any point on the curve, the ratio between the ordinate W and the abscissa D is calculated to form a new W / DD curve.

[0104] 4) Direct inversion of P-wave velocity, adding both WD data and W / DD data into the objective function:

[0105]

[0106] For problems of different orders of magnitude, normalization is performed when calculating the objective function. Where N is the number of depth points, Wr i obs With Wr i cal are the ratios of the wavelength and depth corresponding to the observed and calculated i-th depth point, respectively. max(W i obs ) and max(Wr i obs ) are the maximum values ​​of the observed wavelength and the ratio of wavelength to depth. The objective function is calculated and iterated by genetic algorithm, and the optimal solution is the required longitudinal wave velocity.

[0107] like Figure 2 As shown in the figure, the genetic algorithm inversion of the P-wave velocity of the target geological structure can be divided into the following steps:

[0108] ①First, the model parameter variation range is given, and the P-wave velocity of each layer is randomly determined within this range and combined to generate initial parameters, which are binary-coded to generate the initial population;

[0109] ② Substitute the combined parameters into the objective function (6) and determine their fitness, that is, their ability to survive in the entire population, according to the size of each combined objective function. According to the objective function setting method, the smaller the objective function value, the higher the fitness;

[0110] ③ Determine the selection probability, crossover probability and mutation probability, and perform selection, crossover and mutation operations on the entire population according to the size of fitness to produce a new generation of population, so that the entire population can achieve survival of the fittest;

[0111] ④ Decode the new generation population and calculate the fitness of the new population, that is, the size of the objective function. If the objective function is less than the set value, that is, it is sufficiently converged, the combination of longitudinal wave velocity with the smallest objective function is selected as the optimal solution, or if the iteration exceeds the set maximum number of iterations, the loop ends;

[0112] ⑤If the objective function is greater than the set value and the number of iterations is insufficient, repeat the above operations until the optimal solution that meets the conditions is found.

[0113] Example 2

[0114] The velocity increasing model (model 1) and the non-velocity increasing model (model 2) were designed. There are six overlying strata above the half space of each model with the same layer thickness. The shear wave velocity model is shown in Figure 3 , WD calculation results are as follows Figure 4 , the density of each layer is 2000Kg / m 3 , Poisson's ratio is 0.25, and the Poisson's ratios υ of each layer are set to 0.45, 0.3, 0.3, 0.25, 0.25, and 0.25, respectively. The iteration is set to 15 generations, the population size of each generation is 90, and the P-wave velocity search range is the actual P-wave velocity ±500m / s. The inversion results are as follows: Figure 5-7 shown.

[0115] According to the inversion results, when both curves fit well, the new objective function can be used to effectively realize the direct inversion of P-wave velocity, and the inversion result error is small. The P-wave velocity inversion accuracy is defined as:

[0116]

[0117] Where N is the total number of layers, is the P-wave velocity of the jth layer obtained by inversion, is the true P-wave velocity of the jth layer. According to formula (7), the accuracy of the inversion result can be calculated and shown in Table 1.

[0118] According to the above P-wave velocity inversion model test accuracy analysis, the new inversion algorithm proposed in this study can improve the P-wave velocity inversion accuracy to more than 90%.

[0119] Table 1 Comparison of the accuracy of inversion results of different objective functions

[0120]

[0121] Example 3

[0122] In a specific embodiment 3 of the present invention, Figure 8-10 The implementation process of this method is as follows:

[0123] Figure 8 is the actual surface wave raw data;

[0124] Fig. 9 is the dispersion data extracted by the phase shift method;

[0125] Fig.10 This is the dispersion curve fitting when the genetic algorithm is used for inversion. It can be seen that the dispersion curve fitting is good. The dispersion curve corresponding to the shear wave velocity and thickness structure obtained by inversion is very consistent with the actual extracted dispersion curve, achieving the inversion effect.

[0126] Fig.11 For the inversion results of shear wave velocity and thickness, a larger shear wave velocity and thickness search range was set to obtain a more accurate near-surface shear wave velocity structure;

[0127] Fig.12 Calculated for WD data;

[0128] Fig.13 , 14 The fitting of the WD data and W / DD data corresponding to the P-wave velocity structure obtained by the genetic algorithm for P-wave velocity inversion and the actual WD data and W / DD data shows that the fitting effect is good and a more accurate P-wave velocity structure is obtained.

[0129] Fig.15 is the inversion result of P-wave velocity.

[0130] Compared with the existing speed pickup technology, the present invention has three main advantages:

[0131] (1) The WD information in the traditional method was fully utilized to invert the deep layer;

[0132] (2) The W / DD information is added to improve the inversion accuracy of the surface layer;

[0133] (3) By changing the objective function, the P-wave velocities of all target layers from the surface to the deep layer can be accurately inverted simultaneously.

[0134] Finally, it should be noted that the above is only 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 aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein with equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

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

Claims

1. A near-surface stratum compressional wave velocity inversion method based on surface wave information, characterized in that: The near-surface stratum P-wave velocity inversion method based on surface wave information includes: Step 1, extracting surface wave dispersion data by phase shift method; Step 2, directly invert the near-surface shear wave velocity and thickness based on the dispersion curve; Step 3, calculate the WD curve and the W / DD curve; Step 4: Add both WD data and W / DD data to the objective function to perform direct inversion of P-wave velocity.

2. The method for inversion of near-surface stratum P-wave velocity based on surface wave information according to claim 1, characterized in that: The near-surface stratum P-wave velocity inversion method based on surface wave information further includes, before step 1, acquiring moving source surface wave data from actual seismic acquisition data.

3. The method for inversion of near-surface stratum compressional wave velocity based on surface wave information according to claim 1, characterized in that: In step 1, assuming that the surface wave data is d(x, t), its spectrum can be obtained through Fourier transform: Where x is the distance between the detector and the shot point, A(x,ω) is the amplitude spectrum, and P(x,ω)=exp(-iωx / c R ) is the phase spectrum, x is the detector position, c R (ω) is the actual phase velocity.

4. The method for inversion of near-surface stratum P-wave velocity based on surface wave information according to claim 3, characterized in that: In step 1, since the phase spectrum contains dispersion information, it is necessary to first eliminate the influence of the amplitude spectrum and obtain the normalized spectrum data: In the formula, Normalized spectrum data, that is, phase spectrum data, is obtained by performing multiple phase velocity values ​​c corresponding to a certain frequency. T Substituting and summing, we can get the phase velocity spectrum corresponding to this frequency: Where V(c T ,ω) is the phase velocity spectrum, N is the total number of detectors, n is the detector number, x n is the distance between the nth detector and the shot point; From formula (3), we can see that only when c T Determined as the actual phase velocity c R When , the phase velocity spectrum corresponding to this frequency can reach its maximum value, showing the strongest superposition energy on the dispersion image; repeating the above calculation for all frequencies can perform dispersion extraction on the moving source surface wave seismic records.

5. The method for inversion of near-surface stratum compressional wave velocity based on surface wave information according to claim 1, characterized in that: In step 2, the near-surface shear wave velocity and thickness are directly inverted based on the dispersion curve, and a certain artificially set theoretical near-surface shear wave velocity profile is sought. The dispersion data is obtained by forward calculation based on the shear wave velocity and thickness of each layer. If the curve and the dispersion data extracted from the measured data achieve the best fit, the profile is the shear wave velocity and thickness information.

6. The method for inversion of near-surface stratum P-wave velocity based on surface wave information according to claim 5, characterized in that: In step 2, nonlinear inversion is used to perform a global search to obtain the optimal model parameters. To improve computational efficiency, the objective function definition method used by the fast inversion algorithm is adopted: Where N is the number of frequency points of the extracted dispersion data, m is the total number of models, and f i obs is the frequency, is the phase velocity corresponding to the extracted dispersion curve, Vp j , Vs j ,ρ,h j is the longitudinal wave velocity, shear wave velocity, density and thickness parameters corresponding to the jth model, w i is the weight coefficient corresponding to different frequencies, Ψ j is the objective function value corresponding to the jth model; F vr This is the forward operator obtained by the Haskell-Thomson transfer matrix method.

7. The method for inversion of near-surface stratum compressional wave velocity based on surface wave information according to claim 6, characterized in that: In step 2, a genetic algorithm is used as an inversion method. The steps of using the genetic algorithm to invert the shear wave velocity and thickness of each layer in the target area include: ①First, the model parameter variation range is given, and the combination of shear wave velocity and layer thickness is determined and encoded within this range; ② Substitute the combined parameters into the objective function (4) and determine their fitness, that is, their ability to survive in the entire population, according to the size of each combined objective function. According to the objective function setting method, the smaller the objective function value, the higher the fitness; ③ Determine the selection probability, crossover probability and mutation probability, and perform selection, crossover and mutation operations on the entire population to produce a new generation of population, so that the entire population can achieve survival of the fittest; ④ Repeat the above operations for the new generation population until the optimal solution that meets the conditions is found.

8. The method for inversion of near-surface stratum P-wave velocity based on surface wave information according to claim 1, characterized in that: In step 3, the average shear wave velocity Vsz at a certain depth is calculated based on the acquired shear wave velocity profile: Among them, n is the depth number, h i With V i is the thickness and shear wave velocity corresponding to the depth number.

9. The method for inversion of near-surface stratum P-wave velocity based on surface wave information according to claim 8, characterized in that: In step 3, the fundamental dispersion curve is displayed through the relationship between wavelength and phase velocity, and the depth D and wavelength W corresponding to the same velocity in the shear wave average velocity curve and the dispersion curve are found to obtain the WD relationship curve of the model; for any point on the curve, the ratio between the ordinate W and the abscissa D is calculated to form a new W / DD curve.

10. The method for inversion of near-surface stratum P-wave velocity based on surface wave information according to claim 1, characterized in that: In step 4, both WD data and W / DD data are added to the objective function: For problems of different orders of magnitude, normalization is performed when calculating the objective function; where N is the number of depth points, Wr i obs With Wr i cal are the ratios of the wavelength and depth corresponding to the observed and calculated i-th depth point, respectively; max(W i obs ) and max(Wr i obs ) are the maximum values ​​of the observed wavelength and the ratio of wavelength to depth respectively; the objective function is calculated and iterated by genetic algorithm, and the optimal solution is the required longitudinal wave velocity.

11. The method for inversion of near-surface stratum P-wave velocity based on surface wave information according to claim 10, characterized in that: In step 4, the genetic algorithm inversion of the P-wave velocity of each layer in the target area includes the following steps: ①First, the model parameter variation range is given, and the P-wave velocity of each layer is randomly determined within this range and combined to generate initial parameters, which are binary-coded to generate the initial population; ② Substitute the combined parameters into the objective function (6) and determine their fitness, that is, their ability to survive in the entire population, according to the size of each combined objective function. According to the objective function setting method, the smaller the objective function value, the higher the fitness; ③ Determine the selection probability, crossover probability and mutation probability, and perform selection, crossover and mutation operations on the entire population according to the size of fitness to produce a new generation of population, so that the entire population can achieve survival of the fittest; ④ Decode the new generation population and calculate the fitness of the new population, that is, the size of the objective function. If the objective function is less than the set value, that is, it is sufficiently converged, the longitudinal wave velocity combination with the smallest objective function is selected as the optimal solution, or the iteration exceeds the set maximum number of iterations, then the loop ends; ⑤If the objective function is greater than the set value and the number of iterations is insufficient, repeat the above operations until the optimal solution that meets the conditions is found.

Images