Near-surface structure prediction method for extremely thick loess accumulation area

Abstract

The invention relates to a near-surface structure prediction method for an extremely thick loess accumulation area. The method comprises the following steps: calculating a static correction value; applying the static correction value to an earthquake acquisition original single shot record to obtain single shot first arrival time after static correction, and further obtaining a direct wave offset propagated in a weathered layer and a distribution range of a refracted wave offset propagated in a high-speed layer; the single shot first arrival time is placed in a plane coordinate system with the abscissa being the shot-geophone offset and the ordinate being the first arrival time according to the distribution of the direct wave shot-geophone offset, the reciprocal of the slope is the weathered layer velocity of the shot points, and the weathered layer velocities of all the shot points and detection points are obtained repeatedly; refracted wave velocity analysis is carried out after refracted wave velocity layering, and delay time of each shot point and each detection point is separated; and taking the weathered layer velocities of all the shot points and the weathered layer velocities of all the detection points as initial conditions, and obtaining the weathered layer thickness of each shot point and the weathered layer thickness of each detection point according to delay inversion. According to the method, the requirement of high-quality pre-stack depth migration on the precision of the near-surface velocity model can be met.

Description

Technical Field

[0001] The present invention relates to geophysical exploration methods, and in particular to near-surface structure prediction methods in seismic exploration processing technology. Specifically, the present invention relates to a calculation method for use in onshore two-dimensional and three-dimensional seismic exploration, particularly when conducting near-surface modeling in thick loess areas with stable refraction velocities, to ensure surface model accuracy and reduce construction costs. Background Art

[0002] As we all know, seismic imaging is both a key and challenging aspect of seismic exploration, and the accuracy of the velocity field is crucial for seismic imaging. Since wave-equation-based prestack depth migration became the mainstream processing technology, the importance of velocity models has become increasingly prominent. Their accuracy directly impacts the quality of seismic imaging, and the accuracy of near-surface velocity models is particularly crucial. If errors exist in the accuracy of the near-surface velocity model, these errors will be cumulatively transmitted to the underlying strata, affecting the accuracy of the underlying deep-layer velocity model and further reducing the accuracy of the overall velocity model. This situation not only severely restricts the quality of prestack depth migration seismic imaging but also potentially affects the credibility of the imaging results.

[0003] There are four most commonly used methods for near-surface modeling: model interpolation based on surface survey data, time-depth relationship curve method, first-arrival refraction inversion method, and tomographic inversion method. (1) The near-surface modeling method of interpolating micro-well interpretation results can accurately obtain the near-surface structure at the micro-well survey location, but the accuracy of the near-surface model between different micro-wells cannot be guaranteed. Moreover, for areas with extremely thick weathering layers, micro-well surveys require deep wells (basically between 100 and 350 meters, and most of them exceed 150 meters) to achieve surface survey work. Due to the huge exploration investment, it is difficult to promote and apply it on a large scale. (2) The time-depth relationship curve method is generally applicable to continuous medium areas, especially in desert areas with stable water tables, and the near-surface modeling accuracy is relatively high. This method requires the establishment of a time-depth curve relationship plate, and ensures that the plate can better statistically analyze the laws of change of different near-surface structures in the work area. In thick loess areas, time-depth relationship curves are usually established based on micro-logging survey results. This requires a sufficient number of sample data points to ensure statistical results, that is, a certain number of micro-logging surveys with different weathering layer thicknesses are required. Because the weathering layer in loess areas is extremely thick, micro-logging surveys require huge investment. There are often only a few micro-logging survey points within a work area, which makes it difficult to ensure the sample data requirements. The time-depth relationship curves constructed in this way often cannot meet the surface modeling accuracy requirements. In this case, it is not appropriate to use time-depth relationship curves to construct near-surface models in thick loess areas. (3) The first-arrival refraction inversion method is a relatively mature near-surface modeling and static correction method. In areas where the refractive interface is stable, given a reasonable weathering layer velocity / thickness, the weathering layer thickness / velocity can be obtained with high accuracy, and a high-precision unified reference static correction value can be obtained. However, in actual production, the velocity and thickness parameters are the parameters that exploration workers try their best to obtain, and it is difficult to give a reasonable or accurate weathering layer velocity / thickness. (4) Travel-time ray tomography is not restricted by surface conditions. Based on the first arrival time of surface observations, it can reconstruct a near-surface velocity model that conforms to the actual spatial variation law and obtain high-quality unified datum static corrections. It is currently the most widely used near-surface modeling and static correction method. However, due to the relatively large track spacing within the near-shot offset and the insufficient number of spatial sample points, the travel-time ray tomography inversion is ultimately difficult to accurately depict the extremely shallow layer velocity, resulting in the tomography inversion results being higher than the actual weathering layer velocity and the weathering layer thickness being larger. Moreover, the travel-time ray tomography inversion results need to interpret the near-surface velocity field, and the extracted high-speed top interface has the problem of multi-solution, which makes it difficult to ensure the accuracy of the near-surface model. In exploration areas with extremely thick weathering layers, there is currently no economical, practical, and effective near-surface modeling method.

[0004] A high-quality near-surface model is crucial for ensuring the quality of raw seismic data, improving the accuracy of near-surface models, and enhancing seismic profile imaging. In seismic exploration, particularly in areas with thick regolith, developing a method to construct a high-precision near-surface model without relying on extensive deep-well micro-logging surveys remains a pressing technical challenge. Summary of the Invention

[0005] The present invention aims to address the above-mentioned problems and proposes a near-surface modeling method for thick weathering zones. This method does not rely on too much deep-well micro-logging survey work and can obtain a high-precision near-surface model while reducing exploration investment and saving project costs, providing high-quality basic data for subsequent pre-stack depth migration processing.

[0006] The specific scheme of the present invention is: A method for predicting near-surface structure in thick loess accumulation areas is as follows: Step 1: Calculate the static correction amount to achieve first-arrival flattening of the seismic waves in the weathered layer range in the original single-shot records of early seismic acquisition when the first-arrival time is greater than 0; Step 2: Apply the static correction to the original single-shot record of seismic acquisition to obtain the first arrival time of the single shot after static correction, and then obtain the offset of the direct wave propagating in the weathering layer and the distribution range of the offset of the refracted wave propagating in the high-speed layer; Step 3: Place the first arrival time of a single shot in a plane coordinate system with the offset as the horizontal axis and the first arrival time as the vertical axis based on the distribution of the direct wave offset. Select the offset range and use the least squares method to fit the slope. The inverse of the slope is the weathering layer velocity at the shot point. Repeat this process to obtain the weathering layer velocity at all shot points, and then obtain the weathering layer velocity at all detection points. Step 4: According to the distribution range of the refracted wave offset, the refracted wave velocity is stratified and then analyzed to separate the delay time of each shot point and receiver point. The weathering layer velocity of all shot points and the weathering layer velocity of all receiver points are used as initial conditions, and the weathering layer thickness of each shot point and receiver point is obtained by inversion based on the delay time.

[0007] The specific process of step 1 is as follows: first, calculating the reference static correction value of each physical point; wherein the physical points include shot points and receiver points; performing an overall time shift on the reference static correction values ​​of all physical points so that the first arrival time in the shot record is greater than 0; The reference plane static correction value is obtained by a static correction method, which includes but is not limited to a travel time ray tomography method, a first arrival refraction inversion method, and a time-depth curve method.

[0008] Specifically, the specific process of calculating the reference surface static correction value by travel time ray tomography is as follows: (1) Where: T is the unified reference plane static correction, s; E d is the unified datum elevation, m; E w is the elevation of the high-speed top interface, m; v s is the reference correction speed, in m / s; n is the total number of regolith layers, dimensionless; h i For the i The thickness of the weathering layer, in m; v i For the i The velocity of the regolith, in m / s; τ is the wellhead time, unit is s; The specific time shift value of the overall time shift is an integer of the maximum value of the reference surface static correction.

[0009] Additionally, the weathering layer velocity of the single shot point is the arithmetic mean of the left-branch weathering layer velocity and the right-branch weathering layer velocity.

[0010] The weathering layer velocity of all the detection points is obtained by interpolation using the weathering layer velocity of all the shot points as basic data through the inverse distance weighted method.

[0011] The refraction velocity analysis is achieved by the exchange method, and then the delay time of each shot point and the detection point is separated by the Gauss-Seidel method.

[0012] The calculation formula for the delay is: (2) Where: t Delay is the delay time, unit is s; h is the thickness of the weathering layer, in m; v is the weathering layer velocity, in m / s; v R is the refraction velocity, in m / s; is the critical angle and satisfies sinθ=v / v R ; The calculation formula for the thickness of the weathered layer is: (3).

[0013] The specific effects of the present invention are: The present invention proposes a high-precision near-surface modeling method for thick loess accumulation areas, which effectively solves the technical bottleneck problem of large investment and low model accuracy in near-surface modeling exploration in thick loess accumulation areas. After adopting the present invention, in the process of field exploration, the entire process of deep well micro-logging (the depth of micro-logging is generally between 100~350m, and the average depth is greater than 150m) survey work that is time-consuming, labor-intensive and expensive can be reduced. Indoor surface modeling can be completed conveniently, quickly and with high quality only by the first arrival time of the cannon, thereby improving the accuracy of the surface model and reducing production costs. The results of seismic numerical simulation analysis show that the error of loess thickness obtained by the method of the present invention is between -8.8~8.4m; the error of loess velocity is between -17~25m / s (-2.1%~2.9%) ( Figure 8 ), the surface model accuracy can meet the requirements of high-quality prestack depth migration for the near-surface velocity model accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 Schematic diagram of the tomographic inversion velocity field and high-speed top interface.

[0015] Figure 2 Comparison chart before and after time shift to unify the reference plane.

[0016] Figure 3 Comparison chart of the original single shot (left) and the single shot after static calibration (right).

[0017] Figure 4 This is the fitting diagram of the near-offset and first arrival time after static calibration.

[0018] Figure 5 This is a comparison chart between this method and the theoretical weathering layer velocity.

[0019] Figure 6 This is a comparison chart of the weathering layer thickness between this method and the theoretical one.

[0020] Figure 7 Comparison of the theoretical (left) and near-surface models constructed using this method (right).

[0021] Figure 8 Comparison of weathering layer thickness (left) and velocity error (right) with theoretical data.

[0022] Figure 9 This is a geological digital model diagram for earthquake forward simulation. DETAILED DESCRIPTION

[0023] Specific application cases To verify the effectiveness of this invention and conduct surface model accuracy analysis, actual near-surface structural parameters are required for comparison. However, current technologies are not yet able to accurately obtain near-surface models for field exploration areas. Therefore, this embodiment utilizes a digital model for seismic forward modeling and the synthetic seismic data generated through numerical simulation.

[0024] A seismic forward model of the Ordos Basin's thick loess plateau was constructed based on surface elevation data from a 2D seismic line, field micro-logging results, and velocity fields derived from tomographic inversion. The model's surface elevations are based on field-measured elevation data, ranging from 506.7 to 1254.3 m above sea level. The high-velocity top interface is the 2000 m / s velocity top interface from tomographic inversion of the 2D seismic data. The velocity between the surface and this interface, corrected for the weathered layer's thickness using field micro-logging results, is 514 to 1005 m / s, adopted by the model. The model's weathered layer thickness ranges from 9.1 to 331.0 m. The velocity of high-velocity layer 1 is 2500 m / s, while the velocities of the underlying two layers are 3500 m / s and 4500 m / s, respectively. The geological digital model is 20km long, 3km deep, and has a grid size of 5m×5m ( Figure 9 The excitation and receiving points were evenly distributed on the surface. The numerical simulation used an intermediate excitation and bilateral reception method with a roll-in / roll-out pattern. The observation system used a 15m track spacing, a 60m shot spacing, and an 8092.5-7.5-15-7.5-8092.5 band. Wave equation forward modeling was performed using Ricker wavelets at a 25Hz dominant frequency. The sampling interval was 2ms, and the record length was 6s. A total of 333 synthetic seismic records were simulated and collected. The first arrival time of each shot recorded in the forward simulation was accurately picked, and near-surface modeling was performed using this technical invention. The specific implementation details are as follows.

[0025] 1) Calculate the reference plane static correction values ​​for all physical points (shot points and receiver points); The purpose of calculating the datum statics is to ensure that the first arrivals within the weathered layer of the single-shot record are flattened. In this example, travel-time ray tomography is used to calculate the datum statics. The specific inversion parameters for travel-time ray tomography are as follows: a longitudinal grid size of 15m, a vertical grid size of 5m, a model bottom elevation of -505.7m, a maximum first arrival offset range of 2500m, and 10 inversion iterations. Travel-time ray tomography inversion is performed to obtain the inverted near-surface velocity field. Based on the characteristics of the work area, a velocity top interface of 2000m / s is selected as the high-speed top interface for interpreting the tomographic inversion velocity field. Figure 1 ), thus obtaining the near-surface velocity model, and calculating the unified datum static correction according to formula (1): T , and its value range is between -79.9~390.7ms.

[0026] 2) Static correction of all physical points in time shift In order to ensure that the first arrival time in the single shot record is greater than 0 after the unified reference static correction obtained in step 1) is applied to the single shot record, it is necessary to perform an overall time shift on the static correction of all physical points. The specific time shift value of the static correction is determined by the unified reference static correction. T The maximum value of the static correction is determined by: Generally, the integer of the maximum value of the unified reference plane static correction is used. T The maximum value is 390.7ms, so the static correction value is selected as 400ms, and the value range of the static correction value after time shift is between 320.1~790.7ms ( Figure 2 ).

[0027] 3) Static correction applied to single shot records Apply the time-shifted static correction obtained in step 2) to the original single shot record, and obtain the first arrival time of the single shot after static correction ( Figure 3 The offset distribution ranges of the direct wave propagating in the weathering layer and the refracted wave propagating in the high-velocity layer were determined by using single-shot records (first arrival time of single shots after static calibration). The offset distribution ranges of the direct wave are between -622.5 and 487.5 m, while the offset distribution ranges of the refracted wave are less than -637.5 m and greater than 502.5 m, respectively.

[0028] 4) Calculate the weathering layer velocity at the shot point The single shot first arrival time obtained in step 3) is placed in a plane coordinate system with the offset as the horizontal axis and the first arrival time as the vertical axis according to the distribution of the direct wave offset. The offset range is selected and the slope is obtained by least square fitting. The inverse of the slope is the weathering layer velocity ( Figure 4 A single shot record has a left branch (small direction) and a right branch (large direction). The arithmetic mean of the left and right branch weathering layer velocities is the weathering layer velocity of the shot point. Repeat the operation for each shot point to obtain the weathering layer velocity of all shot points; This example uses a shot point at a horizontal distance of 16747.5m as an example to illustrate the calculation process. The different shot offsets and corresponding first arrival time data at this location are shown in Appendix 1: The least squares method was used to perform linear fitting on the data with offset ranges of -622.5~0m and 0~487.5m, and the independent variable was obtained as offset. x The dependent variable is the first arrival time t The calculation formula is: Left branch: t 左 = -1.0221 x+ 1068.7, with a coefficient of determination of 0.9992; Right branch: t 右 = 1.3149 x + 1042.7, with a coefficient of determination of 0.9998.

[0029] The reciprocal of the slope is the direct wave velocity, so the direct wave velocities are 997.9 m / s (left branch) and 760.5 m / s (right branch), and the arithmetic mean of the two, 879.2 m / s, is the weathering layer velocity at that location. The theoretical weathering layer velocity at the shot point is 877.4 m / s. The weathering layer velocity error obtained by this method is: 879.2-877.4 = 1.8 (m / s), and the error percentage is 1.8÷877.4 = 0.2%. It can be seen that this method has high accuracy.

[0030] 5) Interpolate the weathering layer velocity of all detection points The weathering layer velocity of all the shot points is used as the basic data, and the weathering layer velocity of all the detection points is obtained by interpolation through the inverse distance weighted method. The weathering layer velocity of the detection points obtained by this method ranges from 509 to 1007 m / s, which is almost completely consistent with the theoretical velocity of 514 to 1005 m / s ( Figure 5 ), the overall error of the weathering layer velocity is very small, ranging from -17 to 25 m / s (-2.1% to 2.9%).

[0031] 6) Obtaining the thickness of the weathering layer by refraction inversion The offset range of the refraction velocity distribution is determined in the first-arrival data set. This embodiment focuses on the distribution within the range of 500-3000 m. After accurate refraction velocity stratification, the refraction velocity analysis is performed using the interchange method, and the Gauss-Seidel method is then used to separate the delay time of each shot point and receiver point. Using the weathering layer velocity obtained in step 4) and the weathering layer velocity of the receiver point obtained in step 5) as the initial conditions, the weathering layer thickness calculation formula (3) is derived based on the delay time calculation formula (2), thereby inverting the weathering layer thickness of each shot point and receiver point.

[0032] (2) (3) Taking the shot point at a horizontal distance of 16747.5m as an example, after the refraction velocity analysis using the interchange method, the refraction velocity at this point is v R The delay time of the shot point is 2046 m / s. The Gauss-Seidel method is used to separate the delay time of the shot point. t Delay is 0.204074s, and the weathering layer velocity at this location is vis 879.2 m / s. Substituting these parameters into equation (3), the calculated weathering layer thickness at this shot point is 198.7 m. The theoretical weathering layer thickness at this shot point is 195.3 m. The error in the weathering layer thickness calculated using this method is: 198.7 - 195.3 = 3.4 m, and the error percentage is 3.4 ÷ 198.7 = 1.7%. This shows that this method has high accuracy.

[0033] The operation was repeated for each shot point. Finally, the weathering layer thickness obtained by this method ranged from 2.9 to 328.0 m, which almost completely coincided with the theoretical velocity of 9.1 to 331.0 ( Figure 6 ), the error of the weathering layer thickness is very small, between -8.8 and 8.4 m.

[0034] The parameters obtained in steps 5) and 6) are the near-surface structure ( Figure 7 )parameter.

[0035] The above six steps complete the entire process of near-surface modeling of the thick loess accumulation area.

Claims

1. A method for predicting near-surface structure in thick loess accumulation areas, characterized in that: Here’s how: Step 1: Calculate the static correction amount to achieve first-arrival flattening of the seismic waves in the weathered layer range in the original single-shot records of early seismic acquisition when the first-arrival time is greater than 0; Step 2: Apply the static correction to the original single-shot record of seismic acquisition to obtain the first arrival time of the single shot after static correction, and then obtain the offset of the direct wave propagating in the weathering layer and the distribution range of the offset of the refracted wave propagating in the high-speed layer; Step 3: Place the first arrival time of a single shot in a plane coordinate system with the offset as the horizontal axis and the first arrival time as the vertical axis based on the distribution of the direct wave offset. Select the offset range and use the least squares method to fit the slope. The inverse of the slope is the weathering layer velocity at the shot point. Repeat this process to obtain the weathering layer velocity at all shot points, and then obtain the weathering layer velocity at all detection points. Step 4: According to the distribution range of the refracted wave offset, the refracted wave velocity is stratified and then analyzed to separate the delay time of each shot point and receiver point. The weathering layer velocity of all shot points and the weathering layer velocity of all receiver points are used as initial conditions, and the weathering layer thickness of each shot point and receiver point is obtained by inversion based on the delay time.

2. The method for predicting near-surface structure of a thick loess accumulation area according to claim 1, characterized in that: The specific process of step 1 is as follows: first, the reference plane static correction value of each physical point is calculated; wherein the physical points include shot points and receiver points; and the reference plane static correction values ​​of all physical points are time-shifted as a whole so that the first arrival time in the shot record is greater than 0.

3. The method for predicting near-surface structure of thick loess accumulation areas according to claim 2, characterized in that: The reference plane static correction amount is obtained by a static correction method, which includes but is not limited to a travel time ray tomography method, a first break refraction inversion method, and a time-depth curve method.

4. The method for predicting near-surface structure of a thick loess accumulation area according to claim 3, characterized in that: The specific process of calculating the reference surface static correction by travel-time ray tomography is as follows: Where: T is the unified reference plane static correction, s; E d is the unified datum elevation, m; E w is the elevation of the high-speed top interface, m; v s is the reference plane correction speed, in m / s; n is the total number of regolith layers, dimensionless; h i For the i The thickness of the weathering layer, in m; v i For the i The velocity of the regolith, in m / s; τ is the wellhead time, unit is s.

5. The method for predicting near-surface structure of thick loess accumulation areas according to claim 2, characterized in that: The specific time shift value of the overall time shift is an integer of the maximum value of the reference surface static correction value.

6. The method for predicting near-surface structure of a thick loess accumulation area according to claim 1, characterized in that: The weathering layer velocity of the single shot point is the arithmetic mean of the left branch weathering layer velocity and the right branch weathering layer velocity.

7. The method for predicting near-surface structure of a thick loess accumulation area according to claim 1, characterized in that: The weathering layer velocity of all the detection points is obtained by interpolation using the weathering layer velocity of all the shot points as basic data through the inverse distance weighted method.

8. The method for predicting near-surface structure of a thick loess accumulation area according to claim 1, characterized in that: The refraction velocity analysis is achieved by the exchange method, and then the delay time of each shot point and the detection point is separated by the Gauss-Seidel method.

9. The method for predicting near-surface structure of thick loess accumulation areas according to claim 1, characterized in that: The calculation formula for the delay is: Where: t Delay is the delay time, unit is s; h is the thickness of the weathering layer, in m; v is the weathering layer velocity, in m / s; v R is the refraction velocity, in m / s; is the critical angle and satisfies sinθ=v / v R .

10. The method for predicting near-surface structure of thick loess accumulation areas according to claim 9, characterized in that: The calculation formula for the thickness of the weathered layer is: 。

Images