A method for determining the propagation angle of a wave field based on sliding window-by-wave polarization analysis

Abstract

The present invention discloses a method for determining the wavefield propagation angle based on sliding window-by-wave polarization analysis. Aiming at the calculation requirement of wave propagation angles in three-component multi-wave wavefields, data preprocessing is performed based on VSP three-component seismic data. Then, polarization analysis is carried out on the processed hodograph. After selecting a time window, amplitude information is extracted from the corresponding original components according to the size and position of the window respectively. The data of different components in the same window are integrated to construct a hodograph, and by fitting a straight line, the polarization angle and the wave propagation direction are determined. By sliding the window, the propagation direction of the entire wavefield is simply and effectively obtained. The method flow of the present invention is concise and efficient. Using the sliding window can quickly process large-scale VSP three-component seismic data, significantly improve the data processing efficiency, and can provide high-fidelity wavefield propagation characteristic analysis to assist velocity modeling and wavefield separation, improving the reliability of the results.

Description

Technical Field

[0001] The present invention belongs to the technical field of wave field propagation angle estimation, and particularly relates to a method for determining the wave field propagation angle based on sliding window-by-wave polarization analysis. Background Art

[0002] With the improvement of computing power and the continuous increase in the accuracy requirements of seismic exploration, the quality of seismic wave velocity modeling and wave field processing has become increasingly important in modern refined exploration. In this context, efficiently and accurately calculating the propagation angle of seismic waves has become one of the crucial technologies in seismic wave velocity modeling and wave field analysis. Currently, on the premise of a known underground medium model, the most commonly used technology is the ray tracing method, which can approximate the propagation path of seismic waves by simulating high-frequency rays. Specifically, through the existing velocity model, seismic waves are regarded as high-frequency rays. After calculating their propagation paths, the propagation angles are estimated by recording the characteristics of these paths. The advantage of this method is its relatively high computational efficiency and its ability to effectively simulate the wave propagation process. However, in actual exploration work, the underground medium model is usually not completely clear. Especially in complex geological conditions, there are often certain errors in the inversion results of the velocity model, resulting in limitations in the accuracy of the calculated propagation angles. The uncertainty of underground structures limits the accuracy of ray tracing. Other existing methods based on direction vector techniques (such as the Poynting vector method) are widely used to estimate the propagation direction of seismic waves, but these methods generally have problems such as large computational amounts and complex operation processes. Facing complex areas with unknown geological conditions and the processing requirements of massive exploration data, their applicability and application efficiency are gradually difficult to meet the requirements of production applications.

[0003] Calculating the propagation angle of multi-component seismic data plays a key role in seismic exploration. Seismic exploration aims to identify and depict underground rock layer structures and mineral resources by analyzing the propagation characteristics of seismic waves in underground media. By calculating the propagation angle of multi-component seismic data, it can effectively support the analysis of wave field propagation characteristics, assist in the separation and fine processing of complex wave fields, and is also of great significance for velocity modeling.

[0004] However, the existing methods for calculating the propagation angle of seismic waves pose great challenges in terms of computational amount and velocity model accuracy. In actual exploration, the existing methods face limitations such as low computational efficiency and strong dependence on the velocity model. Therefore, in order to break through these limitations, it is necessary to design a seismic wave propagation angle calculation method that is more suitable for actual exploration needs, has higher computational efficiency, and stronger applicability. Summary of the Invention

[0005] To solve the above technical problems, the present invention provides a method for determining the wave field propagation angle based on sliding window-by-wave polarization analysis, which can not only quickly obtain the polarization directions of waves in different modes, but also calculate the propagation direction of the waves in this mode based on the polarization information.

[0006] The technical solution adopted by the present invention is as follows: A method for determining the wave field propagation angle based on sliding window-by-wave polarization analysis, and the specific steps are as follows:

[0007] S1. Generate seismic data to prepare for data reading, and then perform preprocessing;

[0008] First, construct a velocity model, perform forward modeling on the model data through the code of the finite difference method to obtain the wave field data on the x and z components, that is, the original wave field data. Then, for the original wave field data obtained by forward modeling, use median filtering to separate the wave field according to different apparent velocities to obtain different sub-wave field data.

[0009] By obtaining the original wave field and the separated sub-wave fields to prepare for data reading, and then performing preprocessing, restore the weak-band signals at the interference positions generated due to the limitations of the separation technology in the separated sub-wave fields. Use the signals of the two adjacent normal channels within the upper and lower ranges where the sliding window is located to perform interpolation calculation, that is, within this range, take the result of adding the signals of the two normal channels with their respective weights of 0.5 and assign it to the band range of the abnormal channel to achieve the effect of restoring the interference signal.

[0010] S2. Based on step S1, establish an initial sliding window by using the in-phase axis with the earliest time sorting, and extract the in-phase axis by sliding the window;

[0011] S3. Based on step S2, construct an arrow diagram according to the bands on the two components intercepted by each window, and obtain the polarization angle according to the arrow diagram, that is, dynamically visualize the arrow diagram, and then perform linear fitting and calculate the included angle of the straight lines to achieve the calculation of the polarization angle wave by wave;

[0012] S4. Based on step S3, calculate the polarization angle wave by wave by using the angle of the fitted straight line, and then determine the propagation direction of the reflected wave received by the detector corresponding to each window according to the mutual relationship between the propagation direction and the polarization direction of the waves in different modes, so as to achieve the calculation of the propagation angle wave by wave, and finally complete the calculation of the propagation angle of the full wave field.

[0013] Further, the specific steps of step S2 are as follows:

[0014] An initial sliding window is established by using the in-phase axis with the earliest time sorting. The entire in-phase axis is located at the center of the window. The window contains the waveforms of s periods, where 2 ≤ s ≤ 3. After processing each in-phase axis, the window slides down by a certain number of sampling points until it contains the next in-phase axis. The calculation expression for the sliding distance is as follows:

[0015]

[0016] Among them, x(j) represents the amplitude at the j-th time sampling point, m and k respectively represent the best position of the previous window and the window moving down by k time sampling points, and n represents that there are n traces in the trace gather.

[0017] When maxenergy reaches the maximum value for the first time as the window slides down, the current window position contains the second in-phase axis. At this time, it is k time sampling points away from the original position. And so on, until the sliding window marks all the in-phase axes.

[0018] Furthermore, the specific steps of step S3 are as follows:

[0019] S31. Construct an arrow diagram based on the wavebands on the two components intercepted by each window;

[0020] Each selected time window is applied to two different components respectively to generate two vectors of equal length. The numerical values of one vector are used as the x-axis coordinates, and the corresponding numerical values of the other vector are used as the z-axis coordinates to form a set of points distributed on the x-z plane.

[0021] S32. Based on step S31, dynamically visualize the arrow diagram;

[0022] In the arrow diagram, different colors are assigned to each point according to the time sequence of the data points in the seismic trace, and a color bar is added to create an arrow diagram containing time information.

[0023] S33. Based on step S32, perform linear fitting, calculate the angle between the lines, and realize the calculation of the polarization angle of each wave;

[0024] Perform linear fitting on the points distributed in the x-z coordinate system to generate an optimal zero-crossing fitting line. The fitting algorithm expression is as follows:

[0025] [a1,a2] = PCA([x(i),z(i)]) (2)

[0026] k′ = a2 / a1 (3)

[0027]

[0028] Among them, x(i) and z(i) represent the signals in the corresponding window intervals on two components. The eigenvalues are obtained by performing the PCA (Principal Component Analysis) method on them. An initial slope is determined by the ratio of the two principal component coefficients a1 and a2. This slope is rotated in both positive and negative directions within a certain range to find a best-fit line, and the sum of the straight-line distances from this line to all points is the smallest. The arctangent calculation is performed on the slope k″ of the best-fit line to obtain the polarization angle.

[0029] Furthermore, in the step S4, the full-wavefield propagation angle calculation includes: P-wave and S-wave, specifically as follows:

[0030] (1) P-wave:

[0031] The propagation direction of the P-wave is consistent with its polarization angle. The angle of the fit line represents the incident angle of the P-wave within this time window. The final propagation angle is determined according to whether the wave is an up-going wave or a down-going wave.

[0032] (2) S-wave:

[0033] The polarization direction of the S-wave is perpendicular to the propagation direction, that is, the polarization direction minus π / 2. This line represents the propagation direction. The final propagation angle is determined according to whether the wave is an up-going wave or a down-going wave.

[0034] Advantages of the present invention: The method of the present invention aims at the calculation requirement of the wave propagation angle in the three-component multi-wavefield. Based on the VSP three-component seismic data, data preprocessing is carried out, and then polarization analysis is performed on the processed hodograph. After selecting the time window, the amplitude information is extracted from the corresponding original components respectively according to the size and position of the window. The data of different components in the same window are integrated to construct a hodograph, and the polarization angle and the wave propagation direction are determined through a fit line, and the window is slid, so as to simply and effectively obtain the propagation direction of the entire wavefield. The method flow of the present invention is simple and efficient. Using a sliding window can quickly process large-scale VSP three-component seismic data, significantly improve the data processing efficiency, and can provide high-fidelity wavefield propagation characteristic analysis, assist in velocity modeling and wavefield separation, and improve the reliability of the results.

[0035] The method of the present invention is completely data-driven, fully utilizes the characteristics of multi-component seismic data, extracts wave information of different components, performs polarization analysis, and infers the propagation angle from the polarization information. Algorithms such as weak signal enhancement and angle fitting are also adopted to ensure that the finally calculated propagation angle is more accurate. Compared with the existing methods, the method of the present invention has a significant improvement in the operation speed, the process is more concise, and the accuracy is higher. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 It is a flowchart of a method for determining the wavefield propagation angle based on sliding window-by-wave polarization analysis of the present invention.

[0037] Figure 2 Schematic diagram of velocity model A in the embodiments of the present invention.

[0038] Figure 3 Schematic diagram of velocity model B in the embodiments of the present invention.

[0039] Figure 4 Schematic diagram of the waveform before signal enhancement in the embodiments of the present invention.

[0040] Figure 5 Schematic diagram of the waveform after signal enhancement in the embodiments of the present invention.

[0041] Figure 6 Schematic diagram of the original wave field obtained by forward modeling of velocity model A in the embodiments of the present invention.

[0042] Figure 7 Schematic diagram of different sub-wave fields separated from velocity model A in the embodiments of the present invention.

[0043] Figure 8 Schematic diagram of the original wave field obtained by forward modeling of velocity model B in the embodiments of the present invention.

[0044] Figure 9 Schematic diagram of different sub-wave fields separated from velocity model B in the embodiments of the present invention.

[0045] Figure 10 Schematic diagram of the initial position of the sliding window in the embodiments of the present invention.

[0046] Figure 11 Schematic diagram of the result of the first sliding of the sliding window in the embodiments of the present invention.

[0047] Figure 12 Schematic diagram of a group of corresponding wave bands in the embodiments of the present invention.

[0048] Figure 13 Schematic diagram of the time series vector end graph in the embodiments of the present invention.

[0049] Figure 14 Schematic diagram of the time series vector end graph fitted with a straight line in the embodiments of the present invention.

[0050] Figure 15 Schematic diagram for calculating the propagation angle of the reflection wave isophase axis of velocity model A in the embodiments of the present invention.

[0051] Figure 16 Schematic diagram for calculating the propagation angle of the transmitted wave isophase axis of velocity model A in the embodiments of the present invention.

[0052] Figure 17 Schematic diagram for calculating the propagation angle of the reflection wave isophase axis of velocity model B in the embodiments of the present invention.

[0053] Figure 18 It is a schematic diagram for calculating the propagation angle of the in-phase axis of the transmitted wave of velocity model B in the embodiments of the present invention. Specific embodiments

[0054] The method of the present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0055] As Figure 1 shown, it is a flowchart of a method for determining the wave field propagation angle based on sliding window-by-wave polarization analysis of the present invention, and the specific steps are as follows:

[0056] S1. Generate seismic data to prepare for data reading, and then perform preprocessing;

[0057] First, construct a velocity model, perform forward modeling on the model data through the code of the finite difference method, and obtain the wave field data on the x and z components, that is, the original wave field data. Then, for the original wave field data obtained by forward modeling, use median filtering to separate the wave field according to different apparent velocities to obtain different sub-wave field data (subsequent processing is all performed on the sub-wave field data).

[0058] In this embodiment, the geophone interval is set to 20 m, the depth of the first geophone is 400 m, the sampling time interval is 2 ms, and the frequency of the source wavelet is 40 hz. After the parameter settings are completed, finite difference forward modeling is performed. This embodiment uses two velocity models for simulation verification. The specific velocity models are as Figure 2 and Figure 3 shown, that is, velocity model A and velocity model B, Figure 3 and the velocity model B is more complex.

[0059] Among them, the vertical well is located at 3500 m, and the black five-pointed star represents the source.

[0060] Prepare for data reading by obtaining the original wave field and the separated sub-wave field, and then perform preprocessing. The limitation of the separation method results in local energy discontinuity of the in-phase axis of the sub-wave field. As Figure 4 shown, then restore the weak band signals at the interference positions generated by the limitation of the separation technology in the separated sub-wave field. Since different waves interfere with each other during wave field separation, the wave field separation method based on apparent velocity (such as median filtering, etc.) cannot effectively distinguish the attribution of the energy of the two waves, and the wave signal intensity decays with depth, which will cause several amplitude abnormal bands to appear in some in-phase axes. At this time, considering that the energy loss of seismic waves propagating in the underground medium is a gradual process. Figure 4Only the energy of the last few channels is very weak. To avoid affecting subsequent polarization calculations, interpolation is used to recover this weak signal. That is, interpolation calculations are performed using the signals of the two adjacent normal channels within the upper and lower ranges of the sliding window. Specifically, the result of adding the signals of the two normal channels with respective weights of 0.5 is assigned to the band range of the abnormal channel within this range, achieving the effect of recovering the interference signal. The change in the energy of the wave propagating in the underground medium is continuous, so it is reasonable to perform interpolation on the last few channels, obtaining a waveform diagram of the enhanced signal as shown in Figure 5 . The enhanced data is stored in a new matrix in Matlab, which only participates in the calculation of the polarization angle.

[0061] In this embodiment, according to the Figure 2 shown velocity model A, the forward modeled wavefield is as shown in Figure 6 , and the separated different sub-wavefields are as shown in Figure 7 .

[0062] Among them, Figure 6 (a) is the wavefield data in the z-component, Figure 6 (b) is the wavefield data in the x-component; Figure 7 (a) is the reflected P-to-S wave in the z-component, Figure 7 (b) is the reflected P-to-S wave in the x-component, Figure 7 (c) is the reflected P wave in the z-component, Figure 7 (d) is the reflected P wave in the x-component, Figure 7 (e) is the transmitted P-to-S wave in the z-component, Figure 7 (f) is the transmitted P-to-S wave in the x-component, Figure 7 (g) is the transmitted P wave in the z-component, Figure 7 (h) is the transmitted P wave in the x-component.

[0063] In this embodiment, according to the Figure 3 shown velocity model B, the forward modeled wavefield is as shown in Figure 8 , and the separated different sub-wavefields are as shown in Figure 9 .

[0064] Among them, Figure 8 (a) is the wavefield data in the z-component, Figure 8 (b) is the wavefield data in the x-component; Figure 9 (a) is the reflected P-to-S wave in the z-component, Figure 9 (b) is the reflected P-to-S wave in the x-component, Figure 9 (c) is the reflected P wave in the z-component, Figure 9 (d) is the reflected P wave in the x-component, Figure 9 (e) is the transmitted P-to-S wave in the z-component, Figure 9 (f) is the transmitted P-to-S wave in the x-component, Figure 9 (g) is the transmitted P wave in the z-component, Figure 9(h) is the transmitted P-wave in the x-component.

[0065] S2. Based on step S1, establish an initial sliding window by using the in-phase axis with the earliest time sorting, and extract the in-phase axis by sliding the window.

[0066] S3. Based on step S2, construct an arrowhead diagram according to the wavebands on the two components intercepted by each window, and obtain the polarization angle according to the arrowhead diagram, that is, dynamically visualize the arrowhead diagram, then perform linear fitting and calculate the included angle of the straight line to realize the calculation of the polarization angle wave by wave.

[0067] S4. Based on step S3, calculate the polarization angle wave by wave using the angle of the fitted straight line, and then determine the propagation direction of the reflected wave received by the geophone corresponding to each window according to the mutual relationship between the propagation direction and the polarization direction of waves in different modes, realize the calculation of the propagation angle wave by wave, and finally complete the calculation of the propagation angle of the full wavefield.

[0068] In this embodiment, the specific steps of step S2 are as follows:

[0069] Establish an initial sliding window by using the in-phase axis with the earliest time sorting. The entire in-phase axis is located at the center of the window. The window contains the waveforms of s periods, and 2 ≤ s ≤ 3. After processing each in-phase axis, the window slides down a certain number of sampling points until it contains the next in-phase axis. The calculation expression of the sliding distance is as follows:

[0070]

[0071] Where x(j) represents the amplitude at the j-th time sampling point, m and k respectively represent the best position of the previous window and the window moving down k time sampling points, and n represents that there are n channels in the trace gather.

[0072] When maxenergy reaches the maximum value for the first time as the window slides down, the current window position contains the second in-phase axis. At this time, the distance from the original position is k time sampling points. And so on until the sliding window marks all the in-phase axes.

[0073] In this embodiment, for the separated sub-wavefield, according to the obtained in-phase axis, 30 time sampling points are taken above and below respectively to obtain the initial position of the sliding window as Figure 10 shown, where Figure 10 (a) is the transmitted P-wave in the z-component, Figure 10 (b) is the transmitted P-wave in the x-component. The two red dashed lines above and below in the figure represent the upper and lower edges of the window.

[0074] After calibrating the initial window, the window slides upward. Based on the process of step S2 above, until the window slides to the range of the next in-phase axis, that is, the first sliding result of the sliding window, asFigure 11 as shown, where Figure 11 (a) is the transmitted P-wave in the z-component, Figure 11 (b) is the transmitted P-wave in the x-component. Compared with Figure 10 , Figure 11 It shows that the sliding window moves upward adaptively to mark the range of the second group of event axes. The details of the actual application are as follows: The initial window slides step by step. Since the window slowly moves out of the range of the event axis, the total energy in the window shows a downward trend at this time. When the lower edge of the window reaches the upper interval of the next event axis, the energy in the window should stop decreasing and instead show an upward trend. The window continues to slide. When the energy in the window reaches the maximum value, it means that the event axis here has been found. It should be noted that due to the fact that in VSP multi-component seismic data, there is often an event axis with very weak energy characteristics in the first component and very strong energy in the second component. At this time, if the sliding window is used to calibrate the range of the event axis for only one component of the wave field, there is a risk of missing the calibration interval. It is necessary to perform the operation on both components simultaneously, and then complement the results with each other to finally obtain a series of complete event axis interval ranges.

[0075] In this embodiment, step S3 is specifically as follows:

[0076] S31. Construct an arrow diagram based on the wave bands on the two components intercepted by each window;

[0077] Each selected time window is applied to two different components respectively to generate two vectors of equal length. The numerical values of one vector are used as the x-axis coordinates, and the corresponding numerical values of the other vector are used as the z-axis coordinates to form a set of points distributed on the x-z plane.

[0078] S32. Based on step S31, dynamically visualize the arrow diagram;

[0079] In the arrow diagram, different colors are assigned to each point according to the time sequence of the data points in the seismic trace, and a color bar is added to create an arrow diagram containing time information.

[0080] In this embodiment, for a certain trace in the area delimited by the window, the band information intercepted in the window is used to extract the corresponding bands on other components as well, as Figure 12 shown, indicating the extraction of Figure 2 the band information of the two components of the 81st trace in the window interval of the reflected P-wave field under the velocity model A. The abscissa represents the number of the time sampling point in the window, and the ordinate represents the amplitude. The bands of the two different components are mapped into a two-dimensional amplitude coordinate system, and the colors of the points are marked in chronological order and connected as Figure 13 shown.

[0081] S33. Based on step S32, perform linear fitting, calculate the angle between the lines, and achieve the calculation of the polarization angle of each wave.

[0082] Perform linear fitting on the points distributed in the x-z coordinate system to generate an optimal zero-crossing fitting line. The expression of the fitting algorithm is as follows:

[0083] [a1, a2] = PCA([x(i), z(i)]) (2)

[0084] k′ = a2 / a1 (3)

[0085]

[0086] Among them, x(i) and z(i) represent the signals in the corresponding window intervals on the two components. Perform the PCA principal component analysis method on them to obtain the eigenvalues, determine an initial slope through the ratio of the two principal component coefficients a1 and a2, rotate this slope direction forward and backward within a certain range, and find an optimal fitting line. The sum of the straight-line distances between this line and all points is the smallest. Perform the arctangent calculation on the slope k″ of the optimal fitting line to obtain the polarization angle. This process aims to determine the polarization direction of the seismic wave in the coordinate system composed of two components.

[0087] On the basis of Figure 13 , in this embodiment, all discrete points are linearly fitted according to the method described in step S33. As Figure 14 shown, convert the slope of the fitted line into an angle and save it. By calculating the polarization angle of each trace, the polarization information of the entire event can be obtained. By sliding the window, the polarization information of the full wavefield can be further extracted.

[0088] In this embodiment, in step S4, the calculation of the propagation angle of the full wavefield includes: P-wave and S-wave, specifically as follows:

[0089] (1) P-wave:

[0090] Since the propagation direction of the P-wave is consistent with its polarization angle, the angle of the fitting line represents the incident angle of the P-wave within this time window. Determine the final propagation angle according to whether the wave is an up-going wave or a down-going wave.

[0091] (2) S-wave:

[0092] Since the polarization direction of the S-wave is perpendicular to the propagation direction, that is, subtract π / 2 from the polarization direction. This line represents the propagation direction. Determine the final propagation angle according to whether the wave is an up-going wave or a down-going wave.

[0093] In this embodiment, the corresponding propagation angle is calculated based on the acquired polarization angle information to distinguish the different relationships between the polarization directions and propagation directions of P and S waves. For the S wave, it is necessary to rotate 90° based on its polarization angle. On the basis of the velocity model A shown in Figure 2 , the propagation angle of the reflection wave isochrone is calculated. As shown in Figure 15 , Figure 15 (a)(b)(c)(d) represent calculating the propagation angle for each trace of the isochrone within the calibrated window range, and the final result is as shown in Figure 15 (e). Calculate the propagation angle of the transmitted wave isochrone. As shown in Figure 16 , Figure 16 (a)(b)(c)(d) represent calculating the propagation angle for each trace of the isochrone within the calibrated window range, and the final result is as shown in Figure 16 (e). On the basis of the velocity model B shown in Figure 3 , the same calculation process is carried out, and the results are as shown in Figure 17 and Figure 18 , Figure 17 (a)(b)(c)(d) represent calculating the propagation angle for each trace of the isochrone within the calibrated window range, and the final result is as shown in Figure 17 (e). Calculate the propagation angle of the transmitted wave isochrone. As shown in Figure 18 , Figure 18 (a)(b)(c)(d) represent calculating the propagation angle for each trace of the isochrone within the calibrated window range, and the final result is as shown in Figure 18 (e). It can be seen from the figure that when the continuity of the isochrone is relatively good, the propagation angles calculated for each trace are continuous, which also conforms to the propagation law of seismic waves underground, verifying that the method of the present invention can quickly solve the propagation direction of waves through simple calculations in the face of complex geological conditions. Processing the above series of calibrated sliding windows according to the idea here, the propagation angles of the full wave field are finally obtained.

[0094] In summary, compared with the existing model-based propagation angle calculation methods, the method of the present invention has many advantages such as being applicable to unknown velocity model data, being intuitive and easy to understand, having simple calculations, and high robustness. The hodograph analysis can intuitively present the vibration modes of multi-component seismic data in a rectangular coordinate system. By adopting the strategy of a sliding window, the calculation window slides up and down in time or space, making the whole calculation process clearer and easier to understand. Compared with the existing Poynting vector method, the method steps of the present invention are more concise and the amount of calculation is also greatly reduced. The hodograph calculation of the sliding window only depends on data driving, and the propagation angle can be efficiently obtained without cumbersome steps. The present invention calculates the propagation direction of the entire event by using a sliding window, which can take into account the continuity of wave propagation underground. Even if there are serious waveform anomalies in a few traces, the propagation angles of the front and rear windows can be reasonably estimated to ensure the accuracy of the calculation. By combining the sliding window and hodograph analysis, the method of the present invention can batch calculate the propagation angles of the separated sub-wave fields, realizing fast and efficient calculation and analysis, and providing strong technical support for subsequent steps such as velocity modeling.

[0095] Those of ordinary skill in the art will realize that the embodiments described herein are for helping the reader understand the principles of the present invention and should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those of ordinary skill in the art can make various other specific deformations and combinations without departing from the essence of the present invention according to these technical revelations disclosed in the present invention, and these deformations and combinations are still within the protection scope of the present invention.

Claims

1. A method for determining the wavefield propagation angle based on sliding window-by-wave polarization analysis, the specific steps are as follows: S1. Generate seismic data to prepare for data reading and then perform preprocessing; First, construct a velocity model, perform forward modeling on the model data through the code of the finite difference method to obtain the wavefield data on the x and z components, that is, the original wavefield data; then, for the original wavefield data obtained by forward modeling, use median filtering to separate the wavefield according to different apparent velocities to obtain different sub-wavefield data; By obtaining the original wavefield and the separated sub-wavefield to prepare for data reading and then performing preprocessing, restore the weak-band signals at the interference positions generated by the limitations of the separation technology in the separated sub-wavefield. Use the signals of the two adjacent normal traces within the upper and lower ranges of the sliding window to perform interpolation calculation, that is, within this range, take the result of adding the signals of the two normal traces with their respective weights of 0.5 and assign it to the band range of the abnormal trace to achieve the effect of restoring the interference signal; S2. Based on step S1, establish an initial sliding window by using the in-phase axis with the earliest time sorting, and extract the in-phase axis by sliding the window; S3. Based on step S2, construct an arrow diagram according to the bands on the two components intercepted by each window, and obtain the polarization angle according to the arrow diagram, that is, dynamically visualize the arrow diagram, then perform linear fitting, and calculate the included angle of the straight line to achieve the calculation of the polarization angle wave by wave; S4. Based on step S3, calculate the polarization angle wave by wave using the angle of the fitted straight line, and then determine the propagation direction of the reflected wave received by the geophone corresponding to each window according to the mutual relationship between the propagation direction and the polarization direction of waves in different modes, to achieve the calculation of the propagation angle wave by wave, and finally complete the calculation of the full-wavefield propagation angle.

2. The method for determining the wave field propagation angle based on sliding window-by-wave polarization analysis according to claim 1, wherein The specific steps of step S2 are as follows: Establish an initial sliding window by using the in-phase axis with the earliest time sorting. The entire in-phase axis is located at the center of the window. The window contains the waveforms of s periods, and 2 ≤ s ≤ 3. After processing each in-phase axis, the window slides down a certain number of sampling points until it contains the next in-phase axis. The calculation expression for the sliding distance is as follows: Among them, x(j) represents the amplitude at the jth time sampling point, m and k respectively represent the best position of the previous window and the window moving down k time sampling points, and n represents that there are n traces in the trace gather; When maxenergy reaches the maximum value for the first time as the window slides down, the current window position contains the second in-phase axis. At this time, the distance from the original position is k time sampling points, and so on until the sliding window marks all the in-phase axes.

3. The method for determining the wave field propagation angle based on sliding window-by-wave polarization analysis according to claim 1, characterized in that The specific steps of step S3 are as follows: S31. Construct an arrow diagram according to the bands on the two components intercepted by each window; Apply each selected time window to two different components respectively to generate two vectors of equal length; use the values of one vector as the x-axis coordinates and the corresponding values of the other vector as the z-axis coordinates to form a set of points distributed on the x-z plane; S32. Based on step S31, dynamically visualize the arrow diagram; In the hodograph, different colors are assigned to each data point according to the time sequence in the seismic trace, and a color bar is added to create a hodograph containing time information; S33. Based on step S32, perform linear fitting, calculate the angle between the lines, and implement the calculation of the polarization angle of each wave; Perform linear fitting on the points distributed in the x-z coordinate system to generate an optimal zero-crossing fitting line; the fitting algorithm expression is as follows: [a1, a2] = PCA([x(i), z(i)]) (2) k′ = a2 / a1 (3) Where x(i) and z(i) represent the signals in the corresponding window intervals on the two components. The principal component analysis (PCA) method is used to obtain the eigenvalues for them. An initial slope is determined by the ratio of the two principal component coefficients a1 and a2. Rotate this slope direction forward and backward within a certain range to find an optimal fitting line, for which the sum of the linear distances to all points is the smallest; the arctangent calculation is performed on the slope k″ of the optimal fitting line to obtain the polarization angle.

4. A method for determining the wave field propagation angle based on sliding window-by-wave polarization analysis according to claim 1, characterized in that In step S4, the calculation of the full-wavefield propagation angle includes: P waves and S waves, specifically as follows: (1) P waves: The propagation direction of P waves is consistent with their polarization angle. The angle of the fitting line represents the incident angle of P waves within this time window. The final propagation angle is determined according to whether the wave is an up-going wave or a down-going wave; (2) S waves: The polarization direction of S waves is perpendicular to the propagation direction, that is, the polarization direction minus π / 2. This line represents the propagation direction. The final propagation angle is determined according to whether the wave is an up-going wave or a down-going wave.

Images