VSP Reflection Wavefield Rotation and Orientation Method Based on Edge Detection and In-phase Axis Tracking

Abstract

The present invention proposes a VSP reflected wavefield rotation and orientation method based on edge detection and in-phase axis tracking. By simulating the rotation process of three-component geophones, it can achieve the true amplitude migration of waves. Using edge detection technology, the position of the in-phase axis can be intuitively obtained, making the whole analysis process more concise and easy to understand, facilitating understanding and further analysis. The present invention adopts a windowing method to apply corresponding polarization information to rotate each band of the in-phase axis. Combining with edge detection technology, the positions of other in-phase axes can be accurately located, thus significantly improving the signal-to-noise ratio of the calculated wavefield. Through the above improvements, the present invention not only has a significant improvement in operation speed, but also has a more concise process. This method performs excellently in the rotation and orientation processing of multi-component seismic data, providing solid technical support for subsequent steps such as velocity modeling.

Description

Technical Field

[0001] The present invention belongs to the field of oil and gas geophysical exploration, and particularly relates to a VSP reflected wavefield rotation orientation method based on edge detection and isochron tracing. Background Art

[0002] In the fine exploration of deep and complex exploration targets in China, with the continuous improvement of the requirements for wavefield processing accuracy and efficiency, wavefield noise suppression and wavefield separation methods have been continuously optimized and innovated. High-fidelity elastic wavefield separation has become a key technology in multi-wave and multi-component seismic exploration. Existing wavefield separation methods can be roughly divided into two categories: methods based on the kinematic characteristics of the wavefield and methods based on the dynamic characteristics of the wavefield. The first category of methods mainly separates the wavefield at different propagation apparent velocities by analyzing the wavefield in the time-space domain or the transform domain. This type of method usually requires manual setting of the wavefield apparent velocity parameter, but it is difficult to extract accurate parameters under complex wave path conditions. The second category of methods focuses on the polarization characteristics of waves and can effectively separate the wavefield by keeping the phase, amplitude, and polarization characteristics unchanged, usually achieving good results. However, this type of method is usually limited to processing the first arrival wave (i.e., the direct P wave) and cannot fully extract and analyze the propagation characteristics of subsequent complex wavefields (such as reflected P waves and S waves). Such methods often only apply to the polarization characteristics of specific waves and do not cover the entire wavefield comprehensively, resulting in limitations in the analysis of the propagation characteristics of the entire wavefield.

[0003] The vector orientation rotation of the multi-component wavefield is a key technology in seismic exploration. Seismic exploration helps us detect the underground rock formation structure and mineral resources by analyzing the propagation characteristics of seismic waves in underground media. In this process, the spatial rotation orientation technology of the reflected wavefield plays a crucial role in fine seismic imaging and reservoir characterization.

[0004] In the vertical seismic profile (VSP) three-component wavefield, the geophones record the amplitude information of the incident wave on three orthogonal components according to the incident angle. The rotation orientation of the wavefield aims to accurately reconstruct the true amplitude information of the incident wave from these three components, that is, to achieve the signal repositioning of the three-component wavefield. This process not only improves the signal-to-noise ratio (SNR) of the wavefield but also provides solid technical support for subsequent seismic imaging and wavefield data processing. In the vector rotation process of the wavefield, how to effectively trace the isochron and perform accurate three-component rotation orientation directly determines the quality of the final processing result. Traditional time-varying rotation methods usually target the first arrival wave and cannot obtain effective rotation orientation for subsequent other types of wavefields, thus limiting their application in supporting the vector decoupling separation of complex wavefields. Summary of the Invention

[0005] The present invention proposes a VSP reflection wavefield rotation and orientation method based on edge detection and event tracking. This method is data-driven, uses edge detection technology to accurately locate all target events, performs three-component polarization analysis one by one to complete wave-by-wave rotation and orientation, and finally realizes global single-type wavefield decoupling through local expansion. Specifically, by analyzing the propagation angles of shear waves and longitudinal waves of the reflected waves, the rotation directions of the single-type wavefields on each component are determined, so as to obtain a reflection wavefield with maximized energy. Compared with traditional methods, the method of the present invention has significant advantages such as a simple process and high accuracy. Through this method, not only can the accuracy of wavefield separation be improved, but also a more reliable reference basis can be provided for lithology identification and oil and gas exploration in seismic exploration.

[0006] The object of the present invention is achieved by the following technical solutions: A VSP reflection wavefield rotation and orientation method based on edge detection and event tracking includes the following steps:

[0007] Step 1, seismic data preprocessing: Use a wavefield separation method based on apparent velocity to separate the reflected P-wavefield and the reflected P-to-S wavefield on the X and Z components, and read and store all reflected wavelet fields into a matlab matrix; the specific process is as follows:

[0008] (1) Obtain seismic data by geophones. The data obtained by each geophone is called a trace of data. All geophone data are put together to form seismic trace gather data;

[0009] (2) Separate wavefields: Use the physical trend information of events to filter out the reflected P-wavefield and the reflected P-to-S wavefield;

[0010] (3) Data reading: Use the load function in matlab to load the seismic data of the two components respectively and store them in two matrices;

[0011] Step 2, calculate the rotation angles of P-waves and S-waves respectively: The calculation formula is as follows:

[0012] α p =β p

[0013]

[0014] where α p , α s represent the rotation angles of P-waves and S-waves, and β p , β s represent the propagation angles of P-waves and S-waves; for S-waves, when its propagation angle is greater than 0 degrees, subtraction operation is adopted, otherwise addition operation is adopted. All angles represent the angle with the positive direction of the horizontal axis;

[0015] Step 3: Track the event axis using edge detection; this includes the following three steps:

[0016] Step 3-1: enhance weak signals in the wavefield image;

[0017] Step 3-2, perform image edge detection: respectively image the enhanced reflected sub-wavefield using a grayscale image, perform edge detection using the Roberts operator, and obtain a scatter plot of the event axis;

[0018] Step 3-3: Conduct event tracking: Traverse each point in the scatter plot. If the upper and lower distances in the same point are within the 30-sample window, they are grouped together, i.e., divided into the same interval. The positions of the points in the same point group are averaged to obtain a new point. Then, curve fitting is performed on the new points obtained in each point group to obtain the event trajectory.

[0019] Step 4. Construct a rotating event window: Using the events obtained by event tracking, take 30 sampling points above and below the event of each trace, and together with the points on the event of the trace, form a rotating window with a window size of 61 sampling points.

[0020] Step 5: Perform vector rotation and orientation on the local band: intercept the wavefield signal as the local band through the rotation window, and rotate the local band using the rotation angle calculated in step 2. Use the rotation information of each band in the wavefield to rotate the entire wavefield.

[0021] Step 6: Perform wavefield spatial rotation and orientation on the flat layer model: For each event in the sub-wavefield, perform vector rotation based on the rotation angle calculated on the two components, and rotate the energy to the component with the largest energy before rotation;

[0022] Step 7: Perform wavefield spatial rotation and orientation on the oblique layer model: perform vector rotation on each event axis in the sub-wavefield based on the rotation angle calculated on the two components, and rotate the energy to the component with the largest energy before rotation.

[0023] The beneficial effects of the present invention are as follows: Starting from the need for vector separation and synthesis of three-component multi-wave fields, the present invention performs data preprocessing based on VSP three-component seismic data. After preprocessing, the rotation angles of P and S waves are calculated based on the propagation angle. After selecting a window, amplitude information is obtained from the corresponding original components according to the window size and position. The amplitudes of the different component bands are rotated to return to their original position. By rotating each band, the three-component vector rotation orientation is finally achieved. This invention can achieve the following excellent performance:

[0024] (1) The present invention can rotate and orient the three components of the same type of wave field to achieve the effect of restoring the true amplitude of the wave, which is beneficial for wave field analysis.

[0025] (2) The present invention can quickly track the event axis through edge detection and efficiently achieve three-component rotational orientation.

[0026] (3) Under the condition of being similar to the traditional time-varying rotation effect, this method is data-driven and does not require prior knowledge of information such as the underground medium structure. To a certain extent, it can obtain a relatively faithful propagation angle of the reflected wave in real data, which is helpful for subsequent imaging. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 is a flowchart of the VSP reflected wavefield spatial rotational orientation method of the present invention;

[0028] Figure 2 is a flat layer velocity model diagram;

[0029] Figure 3 is an original wavefield diagram;

[0030] Figure 4 is a reflected P-to-S wave wavefield diagram; [[ID=2,5]]

[0031] Figure 5 is a reflected P wave wavefield diagram;

[0032] Figure 6 is a waveform diagram before weak signal enhancement;

[0033] Figure 7 is a waveform diagram after weak signal enhancement;

[0034] Figure 8 is a schematic diagram of the SVI technology;

[0035] Figure 9 is a schematic diagram of the event axis tracking process;

[0036] Figure 10 is a time-varying rotation local waveform diagram;

[0037] Figure 11 is a sub-wavefield diagram before rotation;

[0038] Figure 12 is a vector rotation effect diagram;

[0039] Figure 13 is a corresponding event axis schematic diagram;

[0040] Figure 14 is a time-varying rotation result diagram of the reflected P wave wavefield of the flat layer model; <^

[0041] Figure 15 is a vector rotation result diagram of the reflected P-to-S wave wavefield of the flat layer model;

[0042] Figure 16It is a diagram of the dipping layer velocity model;

[0043] Figure 17 It is a diagram of the reflection wavelet field;

[0044] Figure 18 It is a diagram of the original wave field of the dipping layer model;

[0045] Figure 19 It is a diagram of the rotation result of the reflection P - wave to S - wave field vector of the dipping layer model;

[0046] Figure 20 It is a diagram of the time - varying rotation result of the reflection P - wave field of the dipping layer model. Detailed implementation manners

[0047] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0048] As Figure 1 shown, a VSP reflection wave field rotation and orientation method based on edge - detection and event tracking of the present invention includes the following steps:

[0049] Step 1, seismic data pre - processing: Use a wave field separation method based on apparent velocity to separate the reflection P - wave field and the reflection P - to - S - wave field on the X and Z components, and read and store all the reflection wavelet fields into a matlab matrix; the specific process is as follows:

[0050] (1) Forward modeling parameter setting: Set the geophone interval to 20m, the depth of the first geophone to 400m, the sampling time interval to 2ms, and the frequency of the source wavelet to 40hz; obtain seismic data by getting the source wavelet through the geophones. The data obtained by each geophone is called a trace of data, and all the geophone data are put together to form a seismic trace gather data.

[0051] (2) Wave field separation: The present invention analyzes different components under the same type of wavelet field of the VSP three - component wave field. After obtaining the seismic trace gather data by wave field forward modeling, use the physical trend information of the events to filter out the reflection P - wave field and the reflection P - to - S - wave field. This step mainly performs median filtering on the original wave fields on two components, and separates the reflection wavelet field using the event trend. If it is necessary to process the wave fields on three components, randomly select the wave fields on two components and process them according to the method of the present invention, and then process the processed wave fields and the wave fields on the third component according to the method of the present invention. The specific velocity model of this embodiment is as Figure 2 shown, where the vertical well is located at 3500m, the black five - pointed star represents the source, and the original wave field is as Figure 3As shown in the figure, in the figure, (a) is the original wave field of the Z component, and (b) is the original wave field of the X component. The blue and black lines represent the isophase axes representing the trends of the selected reflected P waves and reflected P-to-S converted waves. The wave field diagrams of the reflected P-to-S converted waves and the reflected P waves obtained after final filtering are respectively as Figure 4 and Figure 5 shown in the figure, where (a) is the Z component and (b) is the X component.

[0052] (3) Data reading: In Matlab, use the load function to load the seismic data of the two components respectively and store them in two matrices.

[0053] (4) Window selection: With the isophase axis as the axis, a total of 60 time sampling points are taken up and down as the window size. This window will be used for local calculations of vector rotation.

[0054] Step 2: Calculate the rotation angles of P waves and S waves respectively: Calculate the rotation angles of waves of different modes from the known propagation angles. For example, for P waves, its propagation angle is the same as the rotation angle; for S waves, its propagation angle differs from the rotation angle by 90°.

[0055] When the propagation angle is known, calculate the rotation angles of P and S waves. The calculation formula is as follows:

[0056] α p =β

[0057]

[0058] where α p and α s represent the rotation angles of P waves and S waves, and β p and β s represent the propagation angles of P waves and S waves; for S waves, when its propagation angle is greater than 0 degrees, subtraction operation is adopted, otherwise addition operation is adopted. All angles represent the angles with the positive direction of the horizontal axis.

[0059] Step 3: Use edge detection to track the isophase axis; including the following three steps:

[0060] Step 3-1: Enhance the weak signals in the wave field image; after separating the target sub-wave field, it is found that due to the limitations of wave field separation based on apparent velocity itself, the interference effects of different waves cannot be avoided, and the wave signal intensity decays with depth, which will cause incomplete or energy-approaching-zero bands in the isophase axis, and then lead to local energy discontinuity of the isophase axis in the sub-wave field. Since edge detection is an image-based method, if the isophase axis in the wave field image is discontinuous, it will affect the effect of isophase axis tracking. As Figure 6As shown, where there is an obvious sudden weakening of energy in the last few tails. To address this issue, weak signal enhancement is used, and the enhanced signal is as shown in Figure 7 As shown. At this time, the SVI super-virtual refraction interference method is adopted, and the principle is shown in Figure 8 As shown. The main idea is as follows: When the seismic wave at the source position k propagates to the geophones p and q, if the signal at the q position is interfered, resulting in a very small received signal energy, at this time, p′ can be regarded as a virtual source. The paths traveled by the waves reaching p and q before reaching p′ can be approximately regarded as the same. At this time, only the cross-correlation function of the signals received by the two geophones needs to be calculated, and this function is convolved with the virtual source of p′. Finally, the signal received by the q geophone after enhancement can be obtained. The steps here are to ensure that the edge detection can obtain a relatively complete result and supplement the weak signal, just to make the grayscale image complete and do not participate in the subsequent rotation and orientation calculation.

[0061] Step 3-2: Perform image edge detection. Respectively, the enhanced reflected wavelet field is imaged using a grayscale image, as shown in Figure 9 (a); Since there is almost no noise interference in the simulated data, the Roberts operator is used for edge detection to obtain a scatter plot of the event axis, and the result is shown in Figure 9 (b); The approximate calculation formula of the Roberts operator is as follows:

[0062] G(x,y) = |f x (x,y)| + |f y (x,y)|

[0063] f x (x,y), f y (x,y) respectively represent the partial derivatives of the function with respect to x and y. The specific calculation method is as follows:

[0064] f x (x,y) = f(x,y) - f(x - 1,y - 1)

[0065] f y (x,y) = f(x - 1,y) - f(x,y - 1)

[0066] The above two formulas are the operators for image discretization.

[0067] Figure 9 The red dots in (b) represent the edges of the detected figure. It can be seen from the figure that the upper and lower edges of each event axis are outlined by red dots.

[0068] Step 3-3: Perform event axis tracking. Through steps 3-1 and 3-2, a scatter plot of the event axis can be obtained, as shown in Figure 9(c). Since there may be large offsets in the scatter points, it is necessary to eliminate the points with an offset exceeding 15 sampling points between two traces. Traverse the red points in each trace of the scatter plot. If the vertical distance within the same trace is within a window of 30 sampling points, they are classified into the same point group, that is, divided into the same interval. Take the average of the point positions in the same point group to obtain a new point, which is the green point in Figure (b). As can be seen from the figure, the green points can outline the event itself. Then, perform curve fitting on the new points obtained from each point group to obtain the event trajectory, that is Figure 9 the 4 events in (d). Through this step, the trend of the event can be expressed more accurately, preparing for subsequent rotation. Store the event information in the vector ap. Each position in the ap vector represents the corresponding trace number in the wavefield, and the value corresponds to the number of the sampling point where the target event is located on that trace. Although the SVI technology overcomes the discontinuity of the wavefield events, there are still cases where the edge detection cannot detect continuous weak signals. For example Figure 9 there is actually a slightly weaker event under the first event in (d), but it is not detected. At this time, the strategy adopted is to perform another round of amplitude enhancement on the signal reconstructed by SVI to make its energy close to that of other strong-energy events, and then perform edge detection. The result is as Figure 9 shown in (e). It should be noted that when using the edge detection method, it is necessary to appropriately consider whether to perform SVI and amplitude enhancement according to the situation of the wavefield data, and further improve the accuracy of the edge detection results by flexible adjustment.

[0069] Step 4: Construct an event rotation window: For the same type of sub-wavefields of two components, use the events obtained by event tracking. Take 30 sampling points above and below the event in each trace, together with the points on the event of that trace, to form a rotation window with a window size of 61 sampling points, preparing for subsequent vector rotation.

[0070] Step 5: Perform vector rotation orientation on the local band (the signal intercepted in the rotation window): Intercept the wavefield signal through the rotation window as the local band, and rotate the local band using the rotation angles calculated in Step 2 respectively. The result is as Figure 10 shown. It can be observed that by simulating the rotation of the geophone, the horizontal direction of the geophone is exactly parallel to the polarization direction, so that the wave is maximally detected in the x′ direction, and almost no energy is detected in the z′ direction because it is perpendicular to the polarization direction. The calculation formula is:

[0071] x′ = xcosα + zsinα

[0072] z′ = -xsinα + zcosα

[0073] Among them, x and z are the positions of a point on the X and Z axes in the original coordinate system. When the coordinate system is rotated counterclockwise by α degrees to obtain a new coordinate system, x′ and z′ are the corresponding positions in the new coordinate system. The so-called vector rotation is exactly the same, except that it is simulated as rotating the geophone, so as to obtain a new wavefield record.

[0074] Through local experiments, it can be generalized to the entire wavefield. Using the rotation information of each waveband in the wavefield, the entire wavefield can be rotated to obtain a better three-component amplitude migration effect. The specific operation is to perform vector rotation on the signals in the rotation window of each trace on the isochrone according to the set isochrone rotation window, so as to obtain the rotation result of the wave of the entire isochrone. The core idea is to decompose the vector rotation of the reflected wavefield into the rotation of signal segments as shown in Figure 10 shown. Specifically, each trace on the isochrone intercepts the signal according to the rotation window delimited in the previous step. The amplitudes of each point on the x and z components respectively correspond to the x and y variables in the formula in this step, and the rotation angle within this interval obtained in step 2 corresponds to the α variable (each point within the interval uses this angle value). Each point within the interval is calculated. Finally, the amplitudes of each point will be redistributed on the x and z components, and the amplitude values are x′ and z′. After each point is calculated, the signal segment within the interval is rotated. The result is as shown in Figure 11 shown, where (a) and (b) represent the original reflected P-wave wavefields on two components. After rotation and merging, the reflected P-wave is almost transferred to the component with the maximum energy, that is, as shown in Figure 12 (a) and (b). One additional thing to note is that when performing rotation orientation, first compare the energy magnitudes of the two components. The wavefield of the weak energy component can be rotated to the strong energy component first. This orientation strategy is better because when there is a certain error in the given propagation angle, the remaining energy generated by the rotation of the weak energy component is generally smaller.

[0075] Step 6: Perform wavefield spatial rotation orientation on the flat layer model: Obtain the reflected wavefield through forward modeling and median filtering; perform vector rotation on each isochrone in the sub-wavefield based on calculating the rotation angles on two components, and rotate the energy to the component with the maximum energy before rotation; for the flat layer model, as shown in Figure 2 shown, the reflected wavefield obtained through forward modeling and median filtering is as shown in Figure 4 、 Figure 5 shown, and the corresponding isochrones are as shown in Figure 13 shown. Perform vector rotation on each isochrone in the sub-wavefield based on calculating the rotation angles on two components, and rotate the energy to the component with the maximum energy before rotation. The effect is as shown in Figure 14 、 Figure 15 shown, Figure 4 、Figure 5 The wave fields representing two components before rotation, in Figure 14 (a) and (b) in Figure 15 show the effects on the two rotated components. Similarly, it is not difficult to find that on the wave field generated by the flat layer model, after rotation, the energy is well transferred to one of the components, and the signal has a high fidelity.

[0076] Step 7: Perform wave field spatial rotation orientation on the inclined layer model: For each isochrone in the wavelet field, vector rotation is performed based on the rotation angles calculated on two components, and the energy is rotated to the component with the maximum energy before rotation. For the inclined layer model, as Figure 16 shown, with the same source at a horizontal distance of 5000 m and the well at a horizontal distance of 3500 m, geophones are placed starting from a depth of 400 m in the well at an interval of 20 m, and a total of 161 geophones are placed. The original wave field after forward modeling is as Figure 18 shown. The two differently colored lines represent the isochrones of different types of waves. The reflected wave field obtained after median filtering is as Figure 17 shown. For each isochrone in the wavelet field, vector rotation is performed based on the rotation angles calculated on two components, and the energy is rotated to the component with the maximum energy before rotation. The effect is as Figure 19 shown, Figure 17 (a), (b), (c), and (d) in Figure 19 represent the wave fields of the two components before rotation, and (a) and (b) in Figure 20 show the effects on the two rotated components. Similarly, it is not difficult to find that on the wave field generated by the inclined layer model, after rotation, the energy is well transferred to one of the components, and the signal has a high fidelity.

[0077] Those of ordinary skill in the art will realize that the embodiments described herein are to assist the reader in understanding the principles of the present invention and should be understood that the scope of protection 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 that do not depart from the essence of the present invention based on these technical revelations disclosed in the present invention, and these deformations and combinations are still within the scope of protection of the present invention.

Claims

1. A VSP reflected wavefield rotation and orientation method based on edge detection and in-phase axis tracking, characterized in that It includes the following steps: Step 1, seismic data preprocessing: Use the wavefield separation method based on apparent velocity to separate the reflected P-wave field and the reflected P-to-S wave field on the X and Z components, and read and store all the reflected wavelet fields into the matlab matrix; the specific process is as follows: (1) Obtain the seismic data by getting the source wavelet with geophones. The data obtained by each geophone is called one trace of data, and all the geophone data are put together to form the seismic trace gather data; (2) Separate the wavefield: Use the physical trend information of the event to filter out the reflected P-wave field and the reflected P-to-S wave field; (3) Data reading: Use the load function in matlab to load the seismic data of the two components respectively and store them in two matrices; Step 2, calculate the rotation angles of P-waves and S-waves respectively: The calculation formula is as follows: α p ＝β p where α p and α s represent the rotation angles of P-wave and S-wave, and β p and β s represent the propagation angles of P-wave and S-wave; for S-wave, when its propagation angle is greater than 0 degree, subtraction operation is adopted, otherwise addition operation is adopted, and all the angles represent the included angles with the positive direction of the horizontal axis; Step 3, use edge detection to trace the event; It includes the following three steps: Step 3-1, enhance the weak signals in the wavefield image; Step 3-2, conduct image edge detection: Image the enhanced reflected wavelet fields respectively using grayscale imaging, and use the Roberts operator for edge detection to obtain the scatter plot of the event; Step 3-3, conduct event tracing: Traverse the points in each trace of the scatter plot. If the vertical distance within the same trace is within the window of 30 sampling points, then list them as the same point group, that is, divide them into the same interval, take the mean value of the point positions of the same point group to get a new point, and then perform curve fitting on the new points obtained from each point group to get the event trajectory; Step 4, construct the event rotation window: Use the event obtained by event tracing. Take 30 sampling points above and below the event in each trace, together with the points of the trace on the event, to form a rotation window with a window size of 61 sampling points; Step 5, perform vector rotation orientation on the local band: Intercept the wavefield signal through the rotation window as the local band, and rotate the local band respectively using the rotation angles calculated in Step 2; Use the rotation information of each band in the wavefield to rotate the entire wavefield; Step 6, perform wavefield spatial rotation orientation on the flat layer model: Perform vector rotation on each event in the wavelet field based on calculating the rotation angles on the two components, and rotate the energy to the component with the maximum energy before rotation; Step 7, perform wavefield spatial rotation orientation on the inclined layer model: Perform vector rotation on each event in the wavelet field based on calculating the rotation angles on the two components, and rotate the energy to the component with the maximum energy before rotation.

Images