Fracture-cave body signal feature enhancement and connectivity detection method

Abstract

The invention provides a signal feature enhancement and connectivity detection method based on a fracture-cavity body, and the method comprises the steps: carrying out the preprocessing of the seismic data from the original seismic data, mainly carrying out the noise suppression through an SVD denoising method based on a Hankel matrix, and carrying out the detection of the connectivity of the seismic data; and then a multi-window coherence scanning method is combined with quadric surface fitting to obtain a high-precision seismic apparent dip angle. Signal separation is carried out by means of a three-dimensional filtering window constructed through inclination angle guidance, a fracture-cavity body signal after feature enhancement is obtained, an envelope of the fracture-cavity body signal is solved, connectivity detection is carried out, and finally the overall trend description of the fracture-cavity body is obtained. Based on the technical scheme of the invention, the anti-noise capability of the data during the calculation of the seismic dip angle is ensured, and the seismic data dip angle attribute with higher precision and higher robustness and the fracture-cavity feature enhancement effect are obtained, so that the fracture-cavity trend is more accurately depicted finally.

Description

Technical Field

[0001] The present invention relates to the technical field of signal pattern recognition, and particularly to a method for enhancing the characteristics of fracture-vug body signals and detecting the connectivity to identify the morphology of underground fracture-vug bodies. Background Art

[0002] In the field of geophysical exploration, fracture-vug bodies, as typical secondary pore structures in carbonate reservoirs, are formed through multiple geological processes such as karstification, tectonic movement, and fluid migration. Such complex geological bodies developed in three-dimensional space exhibit significant multi-scale characteristics - from millimeter-scale dissolution pores to kilometer-scale cave systems, their spatial distribution is both random and regular. This dual property causes them to generate unique energy attenuation patterns and phase distortion characteristics in the seismic wave field. In the oil and gas resource evaluation system, the effective identification of fracture-vug bodies directly affects the quantitative calculation of reservoir parameters. Their advantages are as follows: compared with matrix pores, fracture-vug systems have higher porosity and permeability gradients. Especially under the action of vertical stacking effects, they can form significant oil and gas enrichment areas. However, this geological advantage is also accompanied by exploration challenges - when fracture-vug bodies show a gradual transition with the background rock layer, conventional seismic exploration methods are easily affected by surface interference, underground medium anisotropy, and noise superposition, resulting in insufficient imaging resolution. Statistics show that the recognition accuracy of fracture-vug bodies in traditional seismic data processing workflows is usually low, and there are serious problems of missed detection and misjudgment.

[0003] Fracture-vug bodies generally refer to structures or regions with fractures and holes in underground rock layers. These structures are very important targets in research such as oil and gas exploration and geological surveys because they usually have good reservoir spaces and permeability, making fracture-vug type oil and gas reservoirs a current hot topic in oil and gas resource research and development. However, due to their complex geological structures and special reservoir characteristics, accurately identifying their underground morphology poses certain challenges. It is mainly manifested in the fact that the structure has a certain degree of complexity, heterogeneity, and concealment, resulting in low efficiency and accuracy in the exploration work of fracture-vug bodies. As a main means of describing the morphology of fracture-vug bodies, fracture-vug body signals are usually based on the preprocessed signal attributes in research, and fracture-vug bodies are identified by the differences between fracture-vug bodies and background attributes.

[0004] The prior art usually identifies fracture-vug bodies based on the attributes of preprocessed signals by differentiating the attributes of fracture-vug bodies from those of the background. In recent years, with the rise of machine learning, some people use the method of building neural networks to train network models to achieve the detection and identification of fracture-vug bodies. Deep learning technology is used to predict fracture-vug bodies, and fracture-vug bodies are characterized by automatically interpreting seismic data and multi-attribute analysis. For the description and identification of special underground structures, in addition to machine learning methods, predecessors directly start from seismic attributes such as curvature and dip angle, and separate special structures such as fracture-vug bodies from the formation background through specific algorithms. The prior art also uses a two-dimensional non-normalized cross-correlation scanning method to perform dip angle scanning along two-dimensional seismic lines, abstracting the three-dimensional surface to two dimensions for dip angle estimation, with low accuracy. Summary of the Invention

[0005] In view of the problems in the above prior art, the present invention proposes a method for enhancing the signal characteristics and detecting the connectivity of fracture-vug bodies.

[0006] The present invention provides a complete process for enhancing and identifying the characteristics of fracture-vug body signals in exploration to determine a method for describing the morphological trend of underground fracture-vug bodies. The present invention first preprocesses the three-dimensional seismic data, mainly using SVD singular value decomposition combined with the Hankel matrix for noise suppression. Subsequently, the present invention separates the signals according to median filtering guided by dip angle to obtain the seismic signals of fracture-vug bodies with enhanced characteristics. Among them, a more robust multi-window coherence scanning algorithm is used for dip angle calculation, and a more accurate apparent dip angle is obtained by using quadratic surface fitting. In the signal separation link, the present invention also adds the Bresenham algorithm for dip angle window guidance on the basis of traditional two-dimensional median filtering, and ascends to the three-dimensional level, and calculates using the apparent dip angles in the x and y directions. Due to the influence of the strong background axis, the energy value of the separated fracture-vug body signal is extremely strong or extremely weak, resulting in segmentation of the results. The present invention performs Hilbert transform on it to obtain its energy envelope, making the energy of the fracture-vug body tend to be integral. Finally, connectivity detection is performed on it to connect the points that belong to the trend of the fracture-vug body but are not connected to the main body, and exclude the noise or error scatter points that do not belong to the trend of the fracture-vug body. Make the results more accurate.

[0007] In one embodiment, the method for enhancing the signal characteristics and detecting the connectivity of fracture-vug bodies includes: Step 1), preprocessing the three-dimensional seismic data, and performing noise suppression based on the SVD noise suppression method of the Hankel matrix;

[0008] Step 2), perform singular value decomposition on the three-dimensional Hankel matrix synthesized in Step 1) to obtain the three-dimensional seismic signal after noise suppression; Step 3), after obtaining the three-dimensional seismic signal after noise suppression in Step 2), calculate the discrete dip angle of the three-dimensional seismic signal; Step 4), perform quadratic surface fitting based on the discrete dip angle obtained in Step 3), and combine the least squares method to obtain the fine apparent dip angle; Step 5), enhance the signal characteristics of the fracture-vug body by dip angle guidance according to the fine dip angle calculated in Step 4); Step 6), perform Hilbert transform on the fracture-vug body feature enhancement result obtained in Step 5) to obtain the energy envelope; Step 7), perform connectivity detection on the fracture-vug body section envelope energy body obtained in Step 6).

[0009] In one implementation, in Step 1), select a section in the three-dimensional seismic data and denote it as X i (x, t), and transform it into the frequency domain according to the Fourier transform. The principle is as follows:

[0010] For the seismic signal section X i , after obtaining the frequency domain data through Fourier transform, extract the frequency slice at its frequency f. Then the function value corresponding to this frequency slice s f is expressed as follows: s f = s(f) = [s 1 , s 2 , …, s y T , and use s f to construct a Hankel matrix, then the matrix H can be obtained as follows: In the above formula, k = y - l + 1; for all frequencies, construct Hankel matrices in sequence according to the above steps to obtain H 1 , H 2 , H 3 , …, H x , and merge all the obtained Hankel matrices to obtain the three-dimensional Hankel matrix as follows:

[0011]

[0012] Among them,

[0013] In one implementation, in Step 2), the specific singular value decomposition operation steps are as follows: Denote the synthesized matrix F imn , that is, F i has a size of m×n, where m = p×l and n = p×k; then the m×n order matrix F i ​The singular value decomposition (SVD) of can be transformed into the product of an orthogonal matrix U, an m×n diagonal matrix Σ, and an n×n orthogonal matrix V. The expression is as follows: F i = UΣV T , where: U is composed of the eigenvector of F i F i T . V is composed of the eigenvector of, and Σ is composed of singular values. The singular values are arranged in descending order on the main diagonal of the matrix. Perform singular value decomposition on the Hankel matrix F i , and select a certain proportion of singular values for reconstruction to obtain the reconstructed Hankel matrix F i '. Inverse-transform F i ' to obtain the Hankel matrices H 1 ', H 2 ', H 3 ', …, H x ' of each frequency slice. Then, inverse-transform each transformed frequency slice Hankel matrix to obtain the frequency slices s 1 ', s 2 ', …, s y '. Finally, restore these frequency slices to X i ', and X i ' is the denoising result after the singular value decomposition of the three-dimensional Hankel matrix. The above operations are for a certain profile in the three-dimensional seismic signal. By simply traversing the operations, the three-dimensional seismic signal S'(x, y, t) after noise suppression can be obtained.

[0014] In one implementation, in step 3), the specific calculation formula is as follows:

[0015]

[0016] where C is the coherence value, T is the number of samples in the time direction of the volume element for calculating the coherence value, J is the total number of channels for coherence analysis, θ x , θ y are the apparent dips of the formation in the main survey line direction x and the connecting survey line direction y respectively, u is the number of original seismic channels, u H is the seismic channel data after Hilbert transform, x j , y j are the coordinates from the analysis sample point to the j-th channel.

[0017] In one implementation, in step 4), the maximum coherence value and its corresponding apparent dips in the x and y directions are obtained through quadratic surface fitting of the following formula: C(θ x , θ y ) = aθ x 2 + bθ y2 + cθ x θ y + dθ x + eθ y + f, where the coefficients a, b, c, d, e, f can be determined by the least squares method, and θ can be accurately obtained by the method of finding the extreme value x 、θ y ; among them, θ x and θ y are obtained as follows: θ is obtained by the method of finding the extreme value through quadratic surface fitting x 、θ y , as the apparent dips of the formation at this point in the x and y directions.

[0018] In one embodiment, in step 5), first calculate the window end coordinates: where x 1 , y 1 are the horizontal and vertical coordinates of the current seismic point, l is half of the given filter window length, and θ x is the apparent dip in the x direction of the current point; then calculate d x = |x 2 - x 1 | and d y = |y 2 - y 1 |, compare their magnitudes and take the larger value as the moving direction; initially take the positive x direction as the moving direction, then the moving distance in y is d i , x is incremented by 1 in turn, and d i+1 = d i + k is calculated in turn. When d i is greater than or equal to 0.5, the y coordinate is incremented by 1, otherwise the y coordinate remains unchanged; repeat this step until the end point is reached; determine half of the filter window in the positive x-axis direction and determine half of the filter window in the negative x-axis direction; the algorithm is the same; when a one-dimensional median filter window is formed through the above steps, the window is translated up and down by (n - l) / 2 units each to form a two-dimensional filter window in the x apparent dip direction, where n is the window width; after forming the two-dimensional filter window in the x direction, according to the Bresenham direction determined by the apparent dip in the y direction at the current point, the two-dimensional filter window is translated in this direction to form the final three-dimensional filter window; finally, use the three-dimensional window to traverse all seismic data to obtain the separated seismic background signal data S bak (x, y, t), and after taking the difference, the characteristic enhancement result S cave (x, y, t) of the fracture-vug body is obtained.

[0019] In one embodiment, in step 6), according to the fracture-vug body characteristic enhancement result S obtained in step 5 cave(x, y, t), perform Hilbert transform on it, and find its energy envelope; select the fracture-vug body data S cave A certain profile X of (x, y, t) cave (x, t), and then perform Hilbert transform on each trace data of this profile; let x(t) be a seismic signal of X cave (x, t), and perform Hilbert transform on it according to the following formula: Then the analytic signal of x(t) can be obtained: And its energy envelope is the amplitude of g(t) in the above formula: Traverse all traces of profile X cave (x, t) to obtain the envelope energy body of X cave (x, t), denoted as X A (x, t).

[0020] In one implementation, in step 7), for X A (x, t), perform binarization, and give a threshold threshold; make all points greater than the threshold value be 1, and the rest of the points be 0, denoted as X A2 (x, t); P i represents the coordinates (x A , t i ) of the i-th point in X i (x, t), give a connectivity window with a window size of W ab , indicating that the window size is a×b, and another given num represents the number of points with values in the window; for each point P i , construct a connectivity window W ab centered on it in the calculated dip direction; count the number N ab of points in its window W i ; if N i ≥num, it means that this point is in the main trend area of the fracture-vug body, then connect P i with other points in the window W ab to form a connected region. The connection method is to regard the two points as the starting point and the ending point, and use the Bresenham algorithm to record the values of the connecting points between the two points as 1; if N i <num, it means that there are fewer surrounding points of P i , and it does not belong to the main trend area of the fracture-vug body. Therefore, regard P i as a discrete point, do not connect it, and change the value of P i to 0; repeat the above steps until all points in X A2 (x, t) have passed the connectivity detection, and the connectivity detection is completed. Compared with the prior art, it has at least the following beneficial effects:

[0021] 1. By performing SVD decomposition noise suppression on the Hankel matrix after dimensionality elevation of seismic data, the coherence between seismic effective data is better considered, ensuring the anti-noise ability of data when calculating seismic dip angles.

[0022] 2. The multi-window coherence scanning algorithm is adopted and combined with quadratic surface fitting of apparent dip angles to obtain more accurate and robust seismic data dip angle attributes.

[0023] 3. A three-dimensional median filtering window that simultaneously considers the apparent dip angles in the x and y directions is used to separate the background formation signal and the fracture-vug body signal, and a more sufficient signal separation result can be obtained, thereby achieving a more accurate fracture-vug body feature enhancement effect.

[0024] 4. By obtaining the envelope body of the fracture-vug body signal and performing connectivity detection, the influence of strong axes and discrete points in the formation background is excluded, making the final characterization of the fracture-vug body trend more accurate. Description of the Drawings

[0025] Figure 1 This is a comparison diagram of the SVD decomposition noise suppression effect of the Hankel matrix of seismic data calculated in the present invention after dimensionality elevation;

[0026] Figure 2 This is a schematic diagram of seismic dip angles obtained by the coherence scanning of multiple windows combined with the quadratic surface fitting method in the present invention;

[0027] Figure 3 This is a schematic diagram of fracture-vug body feature enhancement after signal separation using a three-dimensional window and subtracting from the original seismic signal in the present invention;

[0028] Figure 4 This is a schematic diagram of the connectivity detection result after obtaining the envelope of the fracture-vug body signal in the present invention;

[0029] Figure 5 This is a certain seismic profile of the three-dimensional seismic data collected;

[0030] Figure 6 This is the result of dip angle calculation, fracture-vug body signal separation, and connectivity detection after applying the overall process method of the present invention;

[0031] Figure 7 This is the method flow chart of this case. Detailed Embodiment

[0032] In the following, the present invention will be described in more detail based on embodiments and with reference to the drawings. Among them:

[0033] Figures 1 to 4 Applying the method of the present invention to simulated data. Among them Figure 1 (a) is the noise image, Figure 1In (b) is the noise suppression result, Figure 1 In (c) is the removed noise, Figure 2 In (a) is the original seismic data, Figure 2 In (b) is the apparent dip in the x direction calculated from the multi-window coherence scanning dip, Figure 2 In (c) is the dip angle in the y direction calculated from the multi-window coherence scanning dip, Figure 3 In (a) is the original formation signal, Figure 3 In (b) is the fracture-vug body separation signal, Figure 3 In (c) is the fracture-vug body signal label, Figure 4 In (a) is the fracture-vug body signal, Figure 4 In (b) is the signal envelope energy body, Figure 4 In (c) is the fracture-vug body label, Figure 4 In (d) is the connectivity detection result, Figure 5 In the left figure is the measured seismic data, and in the right figure is the enlarged view of the fracture-vug body part, Figure 6 In (a) is the measured seismic data, Figure 6 In (b) is the calculation result of the multi-window coherence scanning dip; Figure 6 In (c) is the separation result of the three-dimensional window fracture-vug body signal; Figure 6 In (d) is the connectivity detection result, Figure 7 is the method flow chart of this case.

[0034] This embodiment describes in detail a method for enhancing and identifying fracture-vug body characteristics based on seismic dip orientation:

[0035] Step 1) First, preprocess the obtained three-dimensional seismic data S(x, y, t), mainly for noise suppression. For example, a SVD noise suppression method based on the Hankel matrix is used. Select a profile in the three-dimensional seismic data and denote it as X i (x, t), and transform it into the frequency domain according to the Fourier transform, Figure 1 are the noise image and the noise suppression result; its principle is as follows:

[0036]

[0037] For the seismic signal profile X i , obtain the frequency domain data through the Fourier transform After that, extract the frequency slice at its frequency f, then the function value corresponding to this frequency slice s f is expressed as follows:

[0038] s f = s(f) = [s 1 , s 2 ,…, sy T

[0039] In one embodiment, using s f Construct a Hankel matrix, then the matrix H can be obtained as shown in the following equation:

[0040]

[0041] In the above equation, k = y - l + 1.

[0042] Through this process method, after performing Fourier transform on the frequency slices of a certain profile of the three-dimensional seismic data, the Hankel matrix H is constructed f . Then for All frequencies are used to construct the Hankel matrix in sequence according to the above steps to obtain H 1 , H 2 , H 3 , …, H x , and all the obtained Hankel matrices are combined to obtain the three-dimensional Hankel matrix as follows:

[0043]

[0044] Among them,

[0045] Step 2) Perform singular value decomposition on the three-dimensional Hankel matrix F i synthesized in step 1). In some embodiments, the specific singular value decomposition operation steps are as follows:

[0046] Denote the synthesized matrix F imn , that is, F i is of size m×n, where m = p×l and n = p×k.

[0047] Then the singular value decomposition SVD of the m×n matrix F i can be reduced to the product of an orthogonal matrix U, an m×n diagonal matrix Σ, and an n×n orthogonal matrix V, and its expression is as follows:

[0048] F i = UΣV T

[0049] Among them: U is composed of the eigenvector of F i F i T , V is composed of the eigenvector of, Σ is composed of singular values, and its singular values are arranged from large to small on the main diagonal of the matrix.

[0050] Define the square root of the matrix eigenvalue as the matrix singular value σ i ​, it is easy to know that the strength of the effective signal energy of the reconstructed 3D Hankel matrix seismic signal is related to the magnitude of the singular value. It can be seen that in the seismic signal, the larger the component corresponding to the singular value, the greater its contribution rate. Then, if we take the components corresponding to the larger singular values to reconstruct F i , and these larger singular values represent the aggregate of the principal components with certain attributes, then the retention of the main effective seismic signal can be obtained, and the removal of the components with less noise or coherence can be achieved.

[0051] In one implementation, the Hankel matrix F i is subjected to singular value decomposition, and a certain proportion of singular values are selected for reconstruction to obtain the reconstructed Hankel matrix F i '. Then, the Hankel matrices H i ' of each frequency slice are inverse-transformed from F 1 ', H 2 ', H 3 ', …, H x '. Then, the frequency slices s 1 ', s 2 ', …, s y ' are inverse-transformed from each transformed frequency slice Hankel matrix. Finally, these frequency slices are restored to X i ', and X i ' is the denoising result after the singular value decomposition of the 3D Hankel matrix.

[0052] The above operations are for a certain profile in the 3D seismic signal. By simply traversing the operations, the 3D seismic signal S'(x, y, t) after noise suppression can be obtained.

[0053] Step 3) After obtaining the 3D seismic signal after noise suppression in step 2), perform discrete dip angle calculation on it.

[0054] Figure 2 are the apparent dip angles in the x and y directions calculated from the original seismic data and the multi-window coherence scanning dip angle. In this embodiment, first, a high-precision discrete dip angle scanning calculation is performed according to the discrete dip angle scanning coherence value calculation formula obtained by extending the similarity scanning method to 3D seismic data processing by Marfurt et al. The specific formula is as follows:

[0055]

[0056] where C is the coherence value, T is the number of samples in the time direction of the coherence value calculation volume element, J is the total number of channels used for coherence analysis, θ x , θ y are the apparent dip angles of the formation in the x and y directions (i.e., the main survey line direction and the connecting survey line direction) respectively, u is the number of original seismic channels, uH is the seismic trace data after Hilbert transform, x j , y j are the coordinates from the analysis sample point to the j-th trace.

[0057] Step 4) Perform quadratic surface fitting based on the discrete dip angles obtained in Step 3), and combine with the least squares method to obtain the fine apparent dip angles.

[0058] The above discrete dip angle scanning method can only calculate discrete θ x 、θ y , which will lead to the loss of fine features and cause certain errors in the calculated apparent dip angles. When calculating the coherence value of a 2D seismic profile, when the formation dip in the window is consistent with the dip of the formation reflection wave background trend surface, the coherence value is the largest, and the coherence value decreases when deviating from this direction. It can be known that it conforms to the functional relationship. Therefore, in order to accurately estimate the dip angle of the formation, in one embodiment, the maximum coherence value and its corresponding apparent dip angles in the x and y directions can be obtained through quadratic surface fitting of the following formula.

[0059] C(θ x , θ y ) = aθ x 2 + bθ y 2 + cθ x θ y + dθ x + eθ y + f

[0060] The coefficients a, b, c, d, e, f in the formula can be determined by the least squares method, and θ x 、θ y can be accurately obtained by the method of finding the extreme value. Among them, θ x and θ y are obtained as follows:

[0061]

[0062] θ x 、θ y are obtained by the method of finding the extreme value through quadratic surface fitting, and are used as the apparent dip angles of the formation at this point in the x and y directions.

[0063] By the above method of multi-window coherence scanning combined with quadratic surface fitting, the apparent dip angles of the main survey line and the connecting side line directions of a certain profile of 3D seismic data can be obtained.

[0064] Step 5) Enhance the signal characteristics of the fracture-cavity body by dip angle guidance according to the fine dip angles calculated in Step 4).

[0065] Figure 3The original formation signal and the fracture-vug body separation signal are shown. In order to be able to utilize the apparent dips in the x and y directions calculated in step 4) above simultaneously to achieve a better separation effect, in the embodiment, it is preferred to construct a three-dimensional median filtering window, and the specific operation is as follows:

[0066] In this embodiment, the Bresenham algorithm idea is adopted so that all points are calculated as integers in the process of obtaining the inclination-guided window, and it can also ensure that the window selection approximately conforms to the inclination angle. Thus, the median filtering accuracy is improved. First, calculate the window end coordinates according to the following formula:

[0067]

[0068] where x 1 , y 1 are the horizontal and vertical coordinates of the current seismic point, l is half of the given filtering window length, and θ x is the apparent dip in the x direction of the current point.

[0069] Furthermore, calculate d x =|x 2 -x 1 | and d y =|y 2 -y 1 |, compare their magnitudes, and take the larger value as the moving direction. Assume that the positive x direction is the moving direction, then the moving distance in y is d i , x is incremented by 1 in turn, and d i+1 =d i +k is calculated in turn. When d i is greater than or equal to 0.5, the y coordinate is incremented by 1, otherwise the y coordinate remains unchanged. Repeat this step until the end point is reached. In this way, half of the filtering window in the positive x-axis direction is determined, and half of the filtering window in the negative x-axis direction also needs to be determined, and the algorithm is the same.

[0070] When a one-dimensional median filtering window is formed through the above steps, the window is translated up and down by (n - l) / 2 units each to form a two-dimensional filtering window in the x apparent dip direction, where n is the window width.

[0071] After forming the two-dimensional filtering window in the x direction, according to the Bresenham direction determined by the apparent dip in the y direction at the current point, the two-dimensional filtering window is translated in this direction to form the final three-dimensional filtering window.

[0072] Finally, use the three-dimensional window to traverse all seismic data to obtain the separated seismic background signal data S bak (x, y, t), and after taking the difference, the characteristic enhancement result S cave (x, y, t) of the fracture-vug body is obtained.

[0073] Step 6) Based on the enhanced result S of the fracture-vug body characteristics obtained in Step 5) cave (x, y, t), perform Hilbert transform on it to obtain its energy envelope.

[0074] In this embodiment, select the fracture-vug body data S cave A certain profile X of (x, y, t) cave (x, t), and then perform Hilbert transform on each trace data of this profile.

[0075] Let x(t) be a trace of seismic signal of X cave (x, t), and perform Hilbert transform on it according to the following formula:

[0076]

[0077] Then the analytic signal of x(t) can be obtained:

[0078]

[0079] And its energy envelope is the amplitude of g(t) in the above formula:

[0080]

[0081] Traverse all traces of profile X cave (x, t), and obtain the envelope energy body of X cave (x, t), denoted as X A (x, t).

[0082] Step 7) Based on the fracture-vug body profile envelope energy body obtained in 6), perform connectivity detection on it. Figure 4 For the fracture-vug body signal, signal envelope energy body, and connectivity detection result.

[0083] Specifically, perform binarization on X A (x, t), and given a threshold threshold. Note that the selection of the threshold threshold here needs to be repeatedly experimented to achieve the best. Make all points greater than the threshold value be 1, and the rest of the points be 0, denoted as X A2 (x, t).

[0084] P i Represents the coordinates (x A , t i ) of the i-th point in X i (x, t), given a connectivity window, the window size is W ab , indicating that the window size is a×b, and another given num represents the number of points with values in the window.

[0085] For each point P i, calculate the dip direction based on this point, and construct a connectivity window W centered on this direction ab . Count the number of points N ab within its window W i .

[0086] If N i ≥num, it indicates that this point is within the main trend area of the fracture-vug body. Then connect P i to other points within the window W ab to form a connected area. The connection method is to regard the two points as the starting point and the ending point, and use the Bresenham algorithm to record the point values on the line between the two points as 1.

[0087] If N i <num, it indicates that there are fewer surrounding points of P i , and it does not belong to the main trend area of the fracture-vug body. Therefore, regard P i as a discrete point, do not connect it, and change the value of P i to 0.

[0088] Repeat the above steps until all points in X A2 (x,t) have passed the connectivity detection, and the connectivity detection is completed.

[0089] Figures 5 to 6 This is the application result of a specific example of 3D seismic data of two large fracture-vug bodies with a burial depth of 3000 - 4000 meters in a certain area in the west.

[0090] Among them, Figure 5 the measured seismic data and the fracture-vug body location map. Figure 6 are the measured seismic data, the calculation results of multi-window coherence scanning dip angles, the separation results of 3D window fracture-vug body signals, and the connectivity detection results. It can be seen that the technical solution disclosed in the present invention realizes the connectivity detection of fracture-vug body seismic data and finally obtains the overall trend characterization of the fracture-vug body. Although the present invention is described with reference to specific embodiments in this article, it should be understood that these embodiments are only examples of the principles and applications of the present invention. Therefore, it should be understood that many modifications can be made to the exemplary embodiments, and other arrangements can be designed, as long as they do not deviate from the spirit and scope of the present invention defined by the appended claims. It should be understood that different dependent claims and the features described in this article can be combined in a manner different from that described in the original claims. It can also be understood that the features described in connection with a single embodiment can be used in other described embodiments.

Claims

1. A method for fracture-cavity signal feature enhancement and connectivity detection, characterized in that: include: Step 1), Preprocess the 3D seismic data and suppress the noise using the SVD noise suppression method based on the Hankel matrix; Step 2), performing singular value decomposition on the three-dimensional Hankel matrix synthesized in step 1); Step 3), after obtaining the three-dimensional seismic signal after noise suppression in step 2), performing discrete dip calculation on the three-dimensional seismic signal; Step 4), performing quadratic surface fitting according to the discrete inclination angle obtained in step 3), and combining the least square method to obtain a precise apparent inclination angle; Step 5), according to the fine dip angle calculated in step 4), the dip-guided fracture-cavity body signal feature enhancement is performed; Step 6), performing Hilbert transform according to the fracture-cavity feature enhancement result obtained in step 5), and obtaining the energy envelope; Step 7), enveloping the energy body according to the cross-section of the fracture body obtained in step 6), and performing connectivity detection.

2. The method for fracture-cavity signal feature enhancement and connectivity detection according to claim 1, characterized in that: In step 1), a section in the 3D seismic data is selected and recorded as X i (x, t), transform it into the frequency domain according to the Fourier transform, the principle is as follows: For seismic signal profile X i , frequency domain data is obtained through Fourier transform After that, the frequency slice is extracted at its frequency f, then the frequency slice s f The corresponding function value representation is as follows: s f =s(f)=[s1,s2,…,s y ] T , using s f Constructing the Hankel matrix, we get the matrix H, as shown below: In the above formula, k=y-l+1; All frequencies are Fourier transformed to construct the Hankel matrix and obtain H1, H2, H3, ..., H x , merge all the obtained Hankel matrices to obtain the three-dimensional Hankel matrix as follows: in, 3. The method for fracture-cavity signal feature enhancement and connectivity detection according to claim 1, characterized in that: In step 2), the specific singular value decomposition operation steps are as follows: The synthesized matrix F imn , that is, F i The size is m×n, where m=p×l, n=p×k; then the m×n matrix F i The singular value decomposition SVD of can be transformed into the product of the orthogonal matrix U, the m×n order diagonal matrix Σ and the n×n order orthogonal matrix V, and its expression is as follows: i =UΣV T , where: U is composed of F i F i T The eigenvalue vectors of V are composed of eigenvalue vectors, Σ is composed of singular values, and the singular values ​​are arranged from large to small on the main diagonal of the matrix; for the Hankel matrix F i Perform singular value decomposition and select a certain proportion of singular values ​​for reconstruction to obtain the reconstructed Hankel matrix F i ', F i 'Inverse transform the Hankel matrix H1', ​​H2', H3', ..., H of each frequency slice x ', and then inversely transform the Hankel matrix of each transformed frequency slice to obtain frequency slices s1', s2', ..., s y ', and finally restore these frequency slices to X i ', X i ' is the denoising result after the three-dimensional Hankel matrix singular value decomposition; after traversing the singular value decomposition operation, the three-dimensional seismic signal S'(x, y, t) after noise suppression is obtained.

4. The method for fracture-cavity signal feature enhancement and connectivity detection according to claim 1, characterized in that: In step 3), the specific calculation formula is as follows: Where C is the coherence value, T is the number of sample points in the time direction of the voxel for calculating the coherence value, J is the total number of channels used for coherence analysis, θ x ,θ y are the apparent dip angles of the strata in the x and y directions, u is the number of original seismic traces, and u H is the seismic trace data after Hilbert transformation, x j ,y j is the coordinate from the analysis sample point to the jth track.

5. The method for fracture-cavity signal feature enhancement and connectivity detection according to claim 1, characterized in that: In step 4), the maximum coherence value and its corresponding x- and y-direction apparent inclination angles are obtained by fitting the quadratic surface as follows: x ,θ y )=aθ x 2 +bθ y 2 +cθ x θ y +dθ x +eθ y +f, where the coefficients a, b, c, d, e, and f can be determined by the least squares method, and θ can be accurately obtained by the method of finding the extreme value. x ,θ y , where θ x and θ y The following formula is obtained: The method of finding the extreme value by fitting the quadratic surface is used to obtain θ x ,θ y , as the apparent dip angle of the formation at that point in the x and y directions.

6. The method for fracture-cavity signal feature enhancement and connectivity detection according to claim 1, characterized in that: In step 5), first calculate the window end point coordinates: Among them, x1, y1 are the horizontal and vertical coordinates of the current earthquake point, l is half of the given filter window length, θ x is the x-direction apparent inclination angle of the current point; and then calculate d x =|x2-x1| and d y =|y2-y1|, compare their sizes and take the larger value as the moving direction; initially take the positive direction of x as the moving direction, then the moving distance on y is d i , x is incremented by 1, and d is calculated sequentially i+1 =d i +k, when d i When the value is greater than or equal to 0.5, the y coordinate is increased by 1, otherwise the y coordinate remains unchanged; repeat this step until the end point is reached; determine half of the filter window in the positive direction of the x-axis, and determine half of the filter window in the negative direction of the x-axis; after forming a one-dimensional median filter window, translate the window up and down by (nl) / 2 units to form a two-dimensional filter window in the direction of the x-direction, where n is the window width; after forming a two-dimensional filter window in the x-direction, translate the two-dimensional filter window in the direction of the Bresenham direction determined by the y-direction oblique angle of the current point to form a final three-dimensional filter window; finally, use the three-dimensional window to traverse all seismic data to obtain the separated seismic background signal data S bak (x, y, t), and the feature enhancement result S of the fracture body is obtained after subtraction cave (x,y,t).

7. The method for fracture-cavity signal feature enhancement and connectivity detection according to claim 1, characterized in that: In step 6), the fracture feature enhancement result S obtained in step 5) is cave (x, y, t), perform Hilbert transform on it and find its energy envelope; Select fracture volume data S cave A certain section X of (x,y,t) cave (x, t), and then perform Hilbert transform on each data of this section; let x(t) be X cave A seismic signal (x, t) is Hilbert transformed according to the following formula: Then the analytical signal of x(t) can be obtained: And its energy envelope is the amplitude of g(t) in the above formula: Traverse Section X cave All paths of (x,t) get X cave The envelope energy body of (x, t) is denoted by X A (x,t).

8. The method for fracture-cavity signal feature enhancement and connectivity detection according to claim 1, characterized in that: In step 7), for X A (x, t) is binarized by giving a threshold value; all points greater than the threshold value are set to 1, and the rest are set to 0, denoted as X A2 (x, t); P i represents the coordinates (x A of the i-th point in X i , t i ). Given a connectivity window with a window size of W ab , which means the window size is a×b, and another given num, representing the number of points with values in the window; for each point P i , according to the dip direction calculated from this point, a connectivity window W ab is constructed centered on it in this direction; count the number of points N ab in its window W i ; if N i ≥num, it means that this point is in the main trend area of the fracture-vug body, then connect P i with other points in the window W ab to form a connected region. The connection method is to regard the two points as the starting point and the ending point, and use the Bresenham algorithm to set the values of the connecting points between the two points to 1; if N i <num, it means that there are fewer surrounding points of P i , and it does not belong to the main trend area of the fracture-vug body. Therefore, regard P i as a discrete point without connection, and change the value of P i to 0; repeat the binarization step until all points in X A2 (x, t) have passed the connectivity detection, and the connectivity detection is completed.

Images