CN112698398B - Space characterization method for deep fracture system - Google Patents

Abstract

The invention discloses a space depiction method for a deep fracture system. According to the invention, through the application of various advanced geophysical means such as seismic forward model research, seismic data fracture enhancement processing, multi-scale fracture identification, structure tensor attribute, amplitude curvature attribute, compressed sensing inversion and the like, the broken solution identification mode is effectively determined, the space spread of sliding fracture is clearly depicted, the development characteristics of different types of broken solution reservoirs are finely depicted, and a sufficient basis is provided for the space depiction and favorable target prediction of the broken solution reservoirs under the control of deep fracture.

Description

Space characterization method for deep fracture system

Technical Field

The invention relates to the field of petroleum exploration, in particular to a space depiction method of a deep fracture system.

Background

In recent years, the 'broken solution' oil and gas reservoirs of the Tarim basin related to the deep and large sliding fracture zones are continuously discovered, a plurality of high-yield oil and gas wells are emerging, and the great potential of the unconventional oil and gas reservoirs is shown. The broken solution hydrocarbon reservoir of the Tarim basin is a series of special hydrocarbon reservoirs distributed along the sliding fracture zone formed by late hydrocarbon filling after a large sliding fracture zone moving for a long time (from a chilla system to a volunteer system) is transformed into a seam-hole reservoir by deep acid solution corrosion. The oil and gas target layer Oregano system is generally buried in depth exceeding 7000m, and has the characteristics of deep burial, strong fracture control, various reservoir types and complex control factors.

At present, researchers develop a series of research works aiming at the identification and description of the broken solution oil and gas reservoir, and mainly develop the research from aspects of the spatial distribution characteristic characterization of a sliding fracture zone, the broken solution earthquake identification mode research, the broken solution geophysical identification attack and the like. Because the buried depth of the broken solution oil and gas reservoir is large and the seismic reflection characteristics are weakened, the conventional technical means are difficult to fully mine the seismic data information, the spatial characterization effect of a fracture zone is affected, and the method is mainly characterized in that the ① broken solution seismic identification mode is unclear, the fracture identification mode, and the crack, hole and cave type different reservoir identification modes are difficult to determine; ② The fault identification precision is low, the polynary is strong, the faults with different scales cannot be distinguished, and the fault source is difficult to be clear; ③ The physical identification precision of the disconnected solution reservoirs is low, different types of disconnected solution reservoirs cannot be effectively distinguished, the spatial distribution characteristics and connectivity of the disconnected solution are difficult to determine, and quantitative description is difficult.

Aiming at the technical difficulties in the prior researches, the invention effectively determines the broken solution identification mode through the application of various advanced geophysical means such as seismic forward model research, seismic data fracture enhancement processing, multi-scale fracture identification, structure tensor attribute, amplitude curvature attribute, compressed sensing inversion and the like, clearly describes the spatial distribution of sliding fracture, finely describes the development characteristics of different types of broken solution reservoirs, and provides sufficient basis for the spatial description and favorable target prediction of the broken solution reservoirs under the control of deep fracture.

Disclosure of Invention

The invention aims to provide a space depiction method of a deep fracture system, which is used for clearly depicting the space spread of sliding fracture, finely describing the development characteristics of different types of broken solution reservoirs and providing a sufficient basis for the space depiction of the broken solution reservoirs under the control of the deep fracture and the prediction of favorable targets.

In order to achieve the above purpose, the present invention adopts the following technical scheme:

the invention provides a space depiction method of a deep fracture system, which comprises the following steps of:

s1, acquiring well drilling data, well logging data and seismic data of a target area;

s2, performing model forward modeling aiming at a sliding fracture-cut-off control reservoir development mode;

S3, screening out sensitive attributes in fracture and reservoir identification through attribute simulation of a model forward result;

S4, performing fault enhancement processing based on an anisotropic diffusion theory on the seismic data obtained in the S1;

s5, detecting the seismic data subjected to fault enhancement in the S4 based on amplitude gradient vector messiness so as to identify medium-small scale fracture;

s6, carrying out large-scale fracture identification based on machine learning on the seismic data subjected to fault enhancement processing in the S4;

s7, carrying out fracture zone boundary identification based on the structure tensor on the seismic data subjected to fault enhancement processing in the S4;

S8, carrying out hole-crack reservoir identification based on amplitude curvature on the seismic data subjected to fault enhancement in the S4;

s9, performing compressed sensing inversion based on the seismic data subjected to fault enhancement processing in the S4 so as to identify a cave type reservoir;

and S10, fusing the results of the steps S5-S9 to finish the three-dimensional depiction of the deep fracture system.

According to the invention, through the application of various advanced geophysical means such as seismic forward model research, seismic data fracture enhancement processing, multi-scale fracture identification, structure tensor attribute, amplitude curvature attribute, compressed sensing inversion and the like, the broken solution identification mode is effectively determined, the space spread of sliding fracture is clearly depicted, the development characteristics of different types of broken solution reservoirs are finely depicted, and a sufficient basis is provided for the space depiction and favorable target prediction of the broken solution reservoirs under the control of deep fracture.

Each step is described in detail below:

S2, performing model forward modeling aiming at the sliding fracture-cut-off control reservoir development mode.

Based on the space depiction method of the deep fracture system, preferably, the S2 specifically comprises the following steps:

Using drilling and logging data, designing seismic forward models with different sizes (the breaking distance is 10m-70m, the breaking band width is 30m-80m, the crack width is 3m-4m, the cavity width is 10m-60 m) and different filling media (the stratum filling speed is 6000m/s-6500m/s, the main trunk breaking filling speed is 5000m/s-5500m/s, the associated breaking and crack filling speed is 5500m/s, and the cavity filling speed is 4000 m/s) on the basis of fine interpretation of the seismic data, performing forward modeling by adopting a wave equation forward modeling method by referencing actual field three-dimensional observation parameters by forward modeling, and obtaining forward seismic data after offset processing; and performing attribute simulation calculation on the forward-modeling seismic data, establishing a quantitative interpretation quantity plate, and correcting the size of an effective reservoir, thereby achieving the purpose of semi-quantitatively describing the walk-slip fracture cavity reservoir.

After the seismic forward model is subjected to migration processing, forward seismic data with higher quality are obtained, attribute simulation calculation is performed on the forward seismic data, response characteristics of seismic attributes to caves, cracks and fractures can be analyzed, seismic sensitivity attributes are determined, and theoretical basis is provided for subsequent comprehensive identification of broken solution seismic attributes.

Based on the space depiction method of the deep fracture system, preferably, the S3 specifically comprises the following steps:

performing pseudo three-dimensional processing on the forward-modeling obtained seismic data to form a pseudo three-dimensional seismic data body, and performing attribute calculation on the pseudo three-dimensional seismic data body; sensitive properties in fracture and reservoir identification are screened out through simulation and selection of various properties of the seismic forward model results.

Based on the deep fracture system spatial characterization method of the invention, preferably, the sensitive attribute in S3 comprises messiness detection based on an amplitude gradient vector, amplitude curvature and structure tensor.

The specific characteristics of the different attributes are as follows: ① The messy detection result has better identification capability on fracture and crack development bands; ② The amplitude curvature range is larger than the actual cave range, the cave and the peripheral solution holes are actually contained, the reflection of the hole type reservoir is realized, and the transverse resolution capability is high; ③ The structure tensor attribute has better depicting ability for the fracture zone boundary.

S4, performing fault enhancement processing based on an anisotropic diffusion theory on the seismic data obtained in the S1.

Background noise often exists in the pure wave superposition data or the offset data after seismic processing, interference is formed on fault identification, and especially in a low signal-to-noise ratio area, targeted filtering is needed to enhance the definition of break points and sections. The main idea of fault enhancement processing based on the anisotropic diffusion theory is to analyze and filter principal components along the construction direction of seismic data acquisition, effectively improve the signal-to-noise ratio of the seismic data, and achieve high-definition imaging of fracture surfaces on the basis of the constraint of discontinuous detection results such as coherent bodies.

Based on the space depiction method of the deep fracture system, preferably, the S4 specifically comprises the following steps:

s4-1, construction azimuth analysis: determining the azimuth of the reflection phase axis through the earthquake inclination angle and the azimuth angle;

S4-2, breakpoint detection: calculating a coherent body and determining the breakpoint position;

S4-3, unfolding principal component analysis and smoothing (Kuwahara) filtering of a reserved boundary under the constraint of the construction and breakpoint detection results, and completing fault enhancement.

S5, detecting the seismic data subjected to fault enhancement in S4 based on amplitude gradient vector messiness so as to identify medium-small scale fracture.

In the past walk-slip fracture identification process, coherent and AFE attributes are generally adopted for carrying out fracture auxiliary depiction, because the angle of a walk-slip fracture section of a Tarim basin is large and seismic imaging is weak, the attributes are difficult to obtain a good fracture identification effect, especially the fracture longitudinal continuity identification difficulty is higher, and the longitudinal development extension degree of the walk-slip fracture is an important basis for judging the fracture source, so that a better fracture identification technology is urgently needed for carrying out higher-precision identification on the walk-slip fracture.

Aiming at the defects of conventional coherence and other technologies in fracture identification, the fracture identification technology based on amplitude gradient vector messiness detection assumes that a fracture surface is a surface in a local area, and searches the messiness of seismic amplitude gradient vectors in a three-dimensional space through each azimuth angle and each inclination angle to find out the surface with the strongest messiness as the fracture position.

Based on the space depiction method of the deep fracture system, preferably, the S5 specifically comprises the following steps:

S5-1, determining a search azimuth and an inclination angle range according to geological characteristics and data characteristics;

s5-2, calculating an amplitude construction tensor matrix in each search direction, and obtaining disorder degree (disorder) by solving characteristic values of the tensor matrix;

s5-3, determining the azimuth and the inclination corresponding to the maximum mess degree as the azimuth and the inclination of the section, and taking the mess degree value as the possibility of fracture.

The method starts from the seismic amplitude data body, directly searches the distribution rule of faults in the three-dimensional space, is simple and efficient, has strong interpretation of faults in horizontal slices or vertical sections, and is a good three-dimensional fault automatic tracking scheme.

The longitudinal continuity of the messy detection result is obviously enhanced, the characterization effect of the sliding fracture is greatly improved, the fault source can be accurately judged, the details of the messy detection result are rich, the method is finer than the characterization of the fracture by a coherent body, and the method has a better identification effect on the small-scale fault, the micro fracture and the crack development.

S6, performing large-scale fracture identification based on machine learning on the seismic data subjected to fault enhancement processing in the S4.

Based on the deep fracture system space depiction method, preferably, the machine learning in S6 adopts a fracture detection algorithm based on CNN image segmentation, and the fracture detection algorithm based on CNN image segmentation realizes fracture identification based on Unet CNN network.

With the rapid development of machine learning technology in recent years, convolutional Neural Network (CNN) technology is widely used in seismic reservoir prediction. Theoretically, the technology is very suitable for automatic fracture tracking. The fracture detection algorithm based on CNN image segmentation converts the image classification problem into the image segmentation problem, realizes high-precision identification of fracture by using a Unet-based CNN network, and has the following main characteristics and advantages: ① The training data can be picked up interactively in a three-dimensional space, and various possible conditions capable of simulating fracture development can be generated through a random model, so that full-automatic supervised neural network learning based on big data is realized; ② The most advanced Unet CNN network is used for solving the image segmentation problem; ③ The use of a GPU solves a number of operational problems.

S7, carrying out fracture zone boundary identification based on the structure tensor on the seismic data subjected to fault enhancement processing in the S4.

A large number of well drilling results show that the broken solution oil reservoir basically develops in a sliding fracture zone, the activity of oil gas outside the fracture zone is obviously reduced, the reservoir forming capability is poor, and even if a good reservoir development exists, the oil gas discovery is difficult to obtain. Therefore, the exploration and development of the disconnected solution oil reservoir are basically carried out around the sliding fracture zone, and the characterization of the fracture zone boundary is particularly important. The early forward attribute simulation shows that the structure tensor has obvious advantages for fracture zone boundary characterization, so that the structure tensor is adopted to identify the fracture zone boundary, and the inside of the fracture zone boundary is the collection of fault, crack, hole, cave and other storage bodies.

The method firstly utilizes the seismic data to calculate the energy gradient vector, then constructs the structure tensor T, finally carries out the eigenvalue decomposition of the structure tensor to obtain three eigenvalues lambda ₁、λ₂、λ₃ and one eigenvector v, and utilizes the obtained three eigenvalues to carry out the geological anomaly identification. The structure tensor can shield the influence of stratum, so that the fracture zone boundary with strong transverse heterogeneity can be effectively identified, and a basis is provided for identification of the fracture zone boundary.

Based on the space depiction method of the deep fracture system, preferably, the S7 specifically comprises the following steps:

S7-1, gradient calculation:

Calculating a gradient vector (directional derivative) of energy of each point of the seismic data subjected to fault enhancement processing in the step S4:

x, y and z are respectively the line, the channel and the time direction, and u is the seismic amplitude; g ₁ is the differential of the seismic amplitude along the x-direction (line), g ₂ is the differential of the seismic amplitude along the y-direction (xline), g ₃ is the differential in the z or time direction, g is the amplitude gradient vector formed by g ₁、g₂、g₃;

s7-2, construction of Structure tensor T

S7-3, carrying out eigenvalue decomposition on the structure tensor T, wherein the result is three eigenvalues lambda ₁、λ₂、λ₃ and eigenvector v, and the eigenvalue lambda ₂ has a good identification effect on the boundary of the fracture zone through screening and comparison, and the eigenvalue lambda ₂ is adopted to identify the fracture boundary.

S8, carrying out hole-fracture reservoir identification based on amplitude curvature on the seismic data subjected to fault enhancement processing in the S4.

The disconnected solution reservoirs are generally divided into three types: cracks, holes, caves. In the fracture zone enveloping surface, fracture development and fault and hole development are closely related, and generally fracture and hole development areas are often associated with denser fractures, so that a fracture type reservoir layer is generally reflected in a mess manner on an earthquake, and is often characterized by abnormal properties related to fracture and hole on the earthquake property. Through forward attribute simulation, we know that the amplitude curvature has a good identification effect on the hole-fracture type reservoir, and through the amplitude curvature attribute, we can effectively identify the hole-type reservoir and fracture zones associated with the periphery of the hole-type reservoir.

The amplitude curvature is derived from the amplitude of the seismic data by a transverse second order derivative. First, the first derivative of the main line direction and the crossline direction is calculated by using the seismic amplitude or energy, and the obtained energy gradient attribute can reflect an abnormal geologic body, which is generally called amplitude energy gradient. And then performing second-order derivation on the curved surface to obtain an amplitude curved surface, and finally calculating the curvature attribute of each amplitude according to the curve fitting. In principle, the amplitude energy gradient is a representation of the edge of the geologic body, so that the spatial distribution of the geologic body cannot be obtained by a threshold. The amplitude curvature converts the amplitude energy gradient into an attribute reflecting the envelope of the geologic body, and the spatial distribution of the geologic body can be obtained through a threshold value.

The amplitude curvature is very sensitive to strong reflection of the beads from the identification effect, holes and crack development zones around the holes can be well identified, and the amplitude curvature identification result shows that the longitudinal development difference of the disconnected control reservoir is large, the reservoir mainly develops in a shallow layer, and the development of a deep reservoir is gradually deteriorated.

Based on the space depiction method of the deep fracture system, preferably, the S8 specifically comprises the following steps:

s8-1, calculating first derivatives of the main line direction and the crossline direction by using the seismic amplitude or energy to obtain an amplitude energy gradient;

s8-2, performing second order derivation on the amplitude energy gradient to obtain an amplitude curved surface;

s8-3, calculating curvature attributes of each amplitude according to the amplitude curved surface fitting so as to identify the holes and crack development zones around the holes.

And S9, performing compressed sensing inversion based on the seismic data subjected to fault enhancement processing in the step S4 to identify the cave type reservoir.

Cave type reservoir identification is a difficulty in broken solution reservoir identification, and due to limitation of seismic resolution, the position and the size of the cave type reservoir are difficult to accurately position in seismic attribute research, so that the cave type reservoir is required to be identified by adopting a high-resolution inversion means. From the conclusion of the physical studies of the fractured-fluid reservoir rock, the difference in resistance remains an effective distinguishing parameter for carbonate reservoirs and non-reservoirs, with lower resistance generally indicating better reservoir physical properties. The cavity type reservoir is filled with shale and fluid, has obvious low impedance characteristics, and can be identified by adopting an inversion method.

The invention adopts the advanced compressed sensing high resolution inversion method in industry to identify the cave type reservoir, solves the sparse inversion problem by using the L1 norm (L2 norm is used for conventional inversion) by using the compressed sensing algorithm, and assumes that the stratum reflection coefficient can be represented by odd dipole decomposition. The algorithm is realized by introducing the wavelet matrix into wedge dictionary deconvolution, so that the accuracy and resolution of the inversion result are improved.

The objective function of the inversion is:

Min[||s-Wr||₂+λ||r||_p]

Where s is the seismic data, W is the wavelet, r is the reflection coefficient, λ is the regularization coefficient, and p is the norm. In compressed sensing inversion p=1.

The compressive sensing inversion uses a compressive sensing theory in the algorithm implementation process, and the obtained wave impedance is higher than the deterministic inversion resolution of the traditional commercial software; meanwhile, the noise immunity and the stability are better, and the elastic parameters obtained by inversion are more accurate.

Through designing a single cave and two adjacent cave forward models distributed longitudinally, inversion tests are carried out, the visible wave impedance inversion has good recognition degree on a cave type reservoir, the signal to noise ratio of a compressed sensing inversion result is higher, the resolution ratio is improved obviously, the position of the cave is more focused, and the position of the cave can be recognized more accurately.

The transverse heterogeneity of the solution reservoir is very strong, and the obtained wave impedance inversion section has certain continuity in the transverse direction, and the heterogeneity characteristic is not obvious enough, which is also determined by the characteristics of the seismic data. In order to highlight the heterogeneity of a broken solution reservoir and eliminate the influence of transverse continuous seismic reflection, the invention carries out residual impedance processing on the wave impedance inversion result.

Based on the space depiction method of the deep fracture system, preferably, the S9 specifically comprises the following steps: and carrying out residual impedance processing on the wave impedance inversion result of the compressed sensing inversion to identify a cave type reservoir.

Based on the method for spatially describing the deep fracture system, preferably, the residual impedance processing process comprises the following steps: and carrying out median filtering on the wave impedance inversion result to obtain a relatively smooth wave impedance body, and then subtracting the relatively smooth wave impedance body from the original wave impedance body to obtain the residual impedance.

And S10, fusing the results of the steps S5-S9 to finish the three-dimensional depiction of the deep fracture system.

The characterization method of the invention is characterized in that the details of the fracture control reservoir body are characterized in a step-by-step progression from the structure tensor to the amplitude curvature to the residual impedance, the structure tensor attribute is used for characterizing the fracture zone boundary and comprises all reservoirs such as caves, holes, cracks and the like, the amplitude curvature comprises the caves and the holes, cracks and the like associated with the periphery of the caves, the residual impedance is the representation of the cavern type reservoir, and the fracture control reservoir body can be precisely characterized by the comprehensive application of the three technologies.

The broken solution space carving is mainly developed on the basis of comprehensive prediction of the broken solution, and utilizes structure tensor attribute to identify the boundary of the broken solution, messy identification fracture, machine learning fracture detection, amplitude curvature identification hole-fracture type reservoir and high-resolution inversion identification cave. Determining threshold values of different attributes by well earthquake calibration, including structure tensor, messiness detection, machine learning fracture detection, amplitude curvature and compressed sensing inversion results, carrying out normalization processing on different earthquake attributes, and then carrying out spatial attribute fusion to complete spatial carving of a broken solution.

Based on the spatial characterization method of the deep fracture system, preferably, the fusing process in S10 includes: determining threshold values of different attributes by well earthquake calibration, and then carrying out normalization processing on the different earthquake attributes to carry out spatial attribute fusion; and carrying out three-dimensional visual perspective analysis on the fused broken solution attributes to complete three-dimensional depiction of the deep fracture system.

Based on the method for spatial characterization of the deep fracture system, the method is preferable to determine the threshold value of the attribute by utilizing the drilling curve and the venting leakage. When the drilling curve meets fracture, hole, cave or crack dense band, the drilling is obviously reduced, and the phenomena of leakage, emptying and the like generally occur, so that the information can be used for determining threshold values of different attributes.

Based on the space depiction method of the deep fracture system, preferably, after the three-dimensional depiction of the deep fracture system is completed in S10, a three-dimensional attribute detection technology (bodycheck) is adopted to obtain the space distribution rule of the broken solution oil and gas reservoir, namely the space distribution of different fracture-cavity units.

According to the invention, through the application of various advanced geophysical means such as seismic forward model research, seismic data fracture enhancement processing, multi-scale fracture identification, structure tensor attribute, amplitude curvature attribute, compressed sensing inversion and the like, the broken solution identification mode is effectively determined, the space spread of sliding fracture is clearly depicted, the development characteristics of different types of broken solution reservoirs are finely depicted, and a sufficient basis is provided for the space depiction and favorable target prediction of the broken solution reservoirs under the control of deep fracture.

Drawings

FIG. 1 is a schematic diagram of a disconnected solution reservoir development.

FIG. 2a is a cross-sectional view of a model one of a walk-slip fracture and a controlled reservoir seismic forward model in an embodiment.

FIG. 2b is a model migration seismic profile of a forward model one of a walk-slip fracture and a controlled-break reservoir seismic in an embodiment.

FIG. 2c is a corresponding actual seismic profile of a walk-slip fracture and a controlled-break reservoir seismic forward model, according to one embodiment.

FIG. 3a is a cross-sectional view of a model II of a walk-slip fracture and a controlled reservoir seismic forward model in an embodiment.

FIG. 3b is a model migration seismic profile of a walk-slip fracture and a controlled-break reservoir seismic forward model two in an embodiment.

FIG. 3c is a corresponding actual seismic profile of a walk-slip fracture and a controlled-break reservoir seismic forward model in an embodiment.

FIG. 4a is a model of forward modeling of a walk-slip fracture and a controlled reservoir seismic event (amplitude curvature) according to an embodiment.

FIG. 4b is a simulation result (clutter detection) of a forward model of a walk-slip fracture and a fracture-controlled reservoir seismic in an embodiment.

FIG. 4c is a model of forward modeling of a walk-slip fracture and a controlled break in reservoir seismic as a result of a seismic attribute simulation (structure tensor) in an embodiment.

FIG. 5a is a graph showing the results of two seismic attribute simulations (amplitude curvature) of a forward model of a walk-slip fracture and a controlled-break reservoir seismic in an embodiment.

FIG. 5b is a simulation result (clutter detection) of two seismic attributes of a forward model of walk-slip fracture and a controlled reservoir seismic in an embodiment.

FIG. 5c is a simulation result (structure tensor) of two seismic attributes of a forward model of a walk-slip fracture and a fracture-controlled reservoir seismic in an embodiment.

FIG. 6 is a flow chart of an anisotropic diffusion tomography enhancement process in an embodiment.

FIG. 7a is a graph showing the contrast of the effect of the anisotropic diffusion tomography in the embodiment.

FIG. 7b is a graph showing the contrast of the effect of the enhanced treatment of the anisotropic diffusion layer in the embodiment.

FIG. 8 is a schematic diagram of a fault detection technique for amplitude gradient vector clutter analysis in accordance with the present invention.

FIG. 9a is a graph showing the fracture recognition effect of the coherent analysis in the embodiment.

FIG. 9b is a graph of fault detection results of amplitude gradient vector clutter analysis in the example.

FIG. 9c is a slice along the top boundary of an Otto-series-room group in a graph of a fracture recognition effect by coherent analysis in an embodiment.

FIG. 9d is a slice along the top boundary of an Otto-chamber group in a tomographic analysis of amplitude gradient vector clutter in the example.

Fig. 10 is a diagram of a CNN network structure based on Unet of the present invention.

FIG. 11a is a plot of clutter detection versus machine learning fracture recognition result profile for an example.

FIG. 11b is a cross-section of a machine-learned test result comparing the clutter test and the machine-learned fracture recognition results in an embodiment.

FIG. 12a is a plot of the clutter detection and machine learning fracture recognition results versus plane contrast along the top boundary of an Otto series-room group in an example.

FIG. 12b is a machine-learned fault slice along the top boundary of an Otto series-room group for clutter detection and machine-learned fracture recognition result plane comparison in an embodiment.

FIG. 12c is a plot of the clutter detection along the bottom boundary of the Zhongbrian system for the clutter detection and machine learning fracture recognition results in the example.

FIG. 12d is a machine-learned fault slice along the midbrium basal boundary for clutter detection and machine-learned fracture recognition result plane comparison in the examples.

FIG. 13a is a superimposed graph of fracture and structure tensors depicting effect of cross-section of the boundary of the sliding fracture zone in an embodiment.

FIG. 13b is a superimposed graph of fracture and structure tensors versus the plane of effect of the boundary of the sliding fracture zone in an example.

FIG. 14a is a superimposed graph of fracture and amplitude curvature versus hole-fracture reservoir identification profile for an example.

FIG. 14b is a superimposed graph of fracture and amplitude curvature versus hole-fracture reservoir identification effect plane for the example.

FIG. 15a is a seismic model of a linear fracture + single cavity in an embodiment.

FIG. 15b is a seismic section view of an embodiment of a linear fracture + single cavity.

FIG. 15c is a wave impedance inversion profile of an embodiment of linear fracture + single cavity.

FIG. 15d is a cross-sectional view of the compressive sensing wave impedance inversion of a linear fracture + single cavity of an embodiment.

FIG. 15e is a seismic model of a linear fracture + dual cavity in an embodiment.

FIG. 15f is a seismic section view of an embodiment of a linear fracture + dual cavity.

FIG. 15g is a wave impedance inversion profile of the linear fracture + dual cavity of the example.

FIG. 15h is a cross-sectional view of the compressive sensing wave impedance inversion of the linear fracture + dual cavity of the embodiment.

FIG. 16a is a structure tensor attribute of compression-aware inversion and seismic attribute comparison in an embodiment.

FIG. 16b is an amplitude curvature attribute of the compression sensing rate inversion and seismic attribute contrast in an embodiment.

FIG. 16c is a compressed sensing residual impedance profile for compression sensing inversion and seismic attribute comparison in an embodiment.

FIG. 17 is a diagram of an embodiment interrupt solution property fusion effect.

FIG. 18a is a schematic illustration of an exemplary broken solution planar layout.

FIG. 18b is a three-dimensional perspective view of an embodiment of an interrupting solution.

FIG. 19a is a schematic view showing a split plane of an embodiment of an interrupted solution cavity unit.

FIG. 19b is a three-dimensional view of an embodiment interrupted solution cavity cell division.

Detailed Description

In order to more clearly illustrate the present invention, the present invention will be further described with reference to preferred embodiments. It is to be understood by persons skilled in the art that the following detailed description is illustrative and not restrictive, and that this invention is not limited to the details given herein.

The embodiment of the invention provides an optimization scheme, and adopts the optimization scheme to carry out the spatial characterization of the deep fracture system by taking a certain target area as an object.

Based on fracture development characteristics and oil gas distribution characteristics of the sliding fracture, according to actual drilling results and data, establishing a seismic geologic model which accords with the spatial distribution characteristics of a carbonate fracture-cavity type reservoir under the control of the sliding fracture, researching the kinematics and dynamics characteristics of the fracture-cavity reservoir by adopting a seismic attribute simulation analysis technology, establishing the relation between reservoir parameters and seismic attributes, determining several seismic attributes with high sensitivity to reservoir parameter changes and combinations thereof, and establishing a set of carbonate fracture-cavity type reservoir model suitable for the target area and a forward modeling method. Based on the method, the full three-dimensional carving of the fracture and the reservoir is performed by combining drilling data such as drilling, logging, geology and the like and applying various geophysical means such as seismic data fracture enhancement, multi-scale fracture identification, fracture-cavity reservoir identification and the like. The space depiction method of the deep fracture system specifically comprises the following steps:

s1, acquiring well drilling data, well logging data and seismic data of a target area;

S2, performing model forward modeling aiming at the sliding fracture-cut-off control reservoir development mode.

The seismic response of the Oregano hole seam reservoir is characterized primarily by strong amplitude, low frequency beaded or clutter reflection. In order to study the response of the hole type reservoir with different scales, geological models with different sizes and different filling media are designed by using drilling and logging data on the basis of fine interpretation of seismic data, and the model experiment method adopts a field three-dimensional observation mode and is similar to an actual field acquisition observation system. And establishing a quantitative interpretation quantity version by utilizing wave equation forward modeling, and correcting the size of an effective reservoir so as to fulfill the aim of semi-quantitatively describing the broken solution hole seam reservoir.

According to the existing research results, the inside of the sliding fracture zone develops caverns, holes and crack type reservoirs, and has the oil gas dredging function (shown in figure 1), and according to the development mode of the broken solution oil gas reservoir, different types of earthquake forward models are designed to research response characteristics of the fracture and the reservoirs on earthquakes.

In a specific application example, different sliding fracture models (as shown in fig. 2a and 3 a) are designed by taking through-well actual seismic data as a reference, vertical intervals between layers are different from 12m to 70m, a plurality of branch faults and cracks are designed from a room group to the vicinity of a sliding fracture zone of a balm group, and a plurality of caverns with widths between 10m and 30m are designed near the fracture zone. And filling different layers at corresponding speeds according to drilling results, wherein the cavity speed is filled at a low speed of 4000m/s, and the fracture zone is filled at 5500 m/s-6200 m/s.

From the forward result of the earthquake, the offset seismic section and the actual seismic section obtained from the geologic model have high similarity (as shown in fig. 2b and 2c, fig. 3b and 3 c), which shows that the designed geologic model and the actual fracture and reservoir development characteristics have high fitness, so that the answer can be used for tracing the evidence, and the earthquake response characteristics of the fracture and reservoir can be determined according to the forward result, and the conclusion is as follows: ① The response of the high and steep main section (giant slit) is mainly linear weak reflection, and the sign wave is broken; the large inclination section response is mainly linear abnormality; ② Secondary and small fracture development zones are predominantly manifested as clutter and weak reflexes; ③ The cavity is mainly formed by beads or messy strong reflection; ④ The boundary between the large surrounding rock and the stitch net is the boundary part of continuous emission and disordered reflection. Therefore, a basis is provided for the earthquake identification sliding fracture system.

And S3, screening out sensitive attributes in fracture and reservoir identification through attribute simulation of a model forward result.

After the seismic forward model is subjected to migration processing, forward seismic data with higher quality are obtained, attribute simulation calculation is performed on the forward seismic data, response characteristics of seismic attributes to caves, cracks and fractures can be analyzed, seismic sensitivity attributes are determined, and theoretical basis is provided for subsequent comprehensive identification of broken solution seismic attributes.

Because many seismic attributes are realized based on a three-dimensional algorithm at present, the first step of working is to perform pseudo three-dimensional processing on two-dimensional forward seismic data, namely, copying the two-dimensional forward seismic data, combining the two-dimensional forward seismic data according to a three-dimensional observation system to form a pseudo three-dimensional seismic data body, and performing attribute calculation on the three-dimensional seismic data body to obtain a more accurate attribute calculation result. Sensitive attributes in later fracture and reservoir identification are screened out through simulation and optimization of various attributes of forward model results, including messy detection based on amplitude gradient vectors, amplitude curvature, structure tensor and the like.

The specific characteristics of the different attributes are as follows: ① The messy detection result has better identification capability on fracture and crack development bands; ② The amplitude curvature range is larger than the actual cave range, the cave and the peripheral karst hole-crack are actually contained, the reflection of a hole-crack type reservoir is realized, and the transverse resolution capability is high; ③ Tensor properties have a better ability to delineate the boundary of the fracture zone (as shown in fig. 4a, 4b, 4c, 5a, 5b, 5 c).

S4, performing fault enhancement processing on the seismic data based on an anisotropic diffusion theory.

Background noise often exists in the pure wave superposition data or the offset data after seismic processing, interference is formed on fault identification, and especially in a low signal-to-noise ratio area, targeted filtering is needed to enhance the definition of break points and sections. The main idea of fault enhancement processing based on the anisotropic diffusion theory is to analyze and filter principal components along the construction direction of seismic data acquisition, effectively improve the signal-to-noise ratio of the seismic data, and achieve high-definition imaging of fracture surfaces on the basis of the constraint of discontinuous detection results such as coherent bodies.

The technology mainly comprises three main steps:

s4-1, construction azimuth analysis: determining the azimuth of the reflection phase axis through the earthquake inclination angle and the azimuth angle;

S4-2, breakpoint detection: calculating a coherent body and determining the breakpoint position;

S4-3, unfolding principal component analysis and smoothing (Kuwahara) filtering of a reserved boundary under the constraint of the construction and breakpoint detection results, and completing fault enhancement.

The operational flow diagram is shown in fig. 6:

And (3) calculating the inclination angle and the azimuth angle of the seismic data acquired in the step (S1), calculating a coherent body, carrying out discontinuity detection, analyzing whether the time window has strong discontinuity according to a coherent body calculation result, if the discontinuity is strong, considering the discontinuity as a fracture development position, not smoothing the seismic data, otherwise, smoothing along the structure, and adjusting the transverse channel number and the longitudinal time window in the calculation process to obtain an ideal result. This process may be iterated until the best results are obtained.

From the comparison of the fault strengthening effect of the seismic data, the signal-to-noise ratio of the processed seismic section is obviously improved, the longitudinal continuity of the fault is obviously enhanced, the multi-resolution of fault interpretation is obviously reduced, and the high-quality seismic data is provided for walk-slip fault identification (shown in fig. 7a and 7 b).

S5, detecting the seismic data subjected to fault enhancement in S4 based on amplitude gradient vector messiness so as to identify medium-small scale fracture.

In the past walk-slip fracture identification process, coherent and AFE attributes are generally adopted for carrying out fracture auxiliary depiction, because the angle of a walk-slip fracture section of a Tarim basin is large and seismic imaging is weak, the attributes are difficult to obtain a good fracture identification effect, especially the fracture longitudinal continuity identification difficulty is higher, and the longitudinal development extension degree of the walk-slip fracture is an important basis for judging the fracture source, so that a better fracture identification technology is urgently needed for carrying out higher-precision identification on the walk-slip fracture.

As shown in fig. 8, the fracture identification technology based on the amplitude gradient vector messiness detection aims at the defect of the conventional coherence technology in fracture identification, and the fracture position is found out by searching the messiness of the seismic amplitude gradient vector through each azimuth and inclination in the three-dimensional space on the assumption that the fracture surface is a surface in a local area.

Based on the space depiction method of the deep fracture system, preferably, the S5 specifically comprises the following steps:

s5-1, determining a search azimuth and an inclination angle range according to geological characteristics and data characteristics, and calculating azimuth GST messiness;

s5-2, calculating an amplitude construction tensor matrix in each searching direction, obtaining disorder degree (disorder) by solving a characteristic value of the tensor matrix, performing three-dimensional space correlation, and searching the space direction of a fracture surface, wherein the noise immunity and the stability are improved;

S5-3, searching the fracture surface azimuth, diffusing energy along the fracture surface, and determining the azimuth and the inclination corresponding to the maximum mess as the azimuth and the inclination of the fracture surface, wherein the mess value is used as the possibility of fracture (shown in figure 8).

The method starts from the seismic amplitude data body, directly searches the distribution rule of faults in the three-dimensional space, is simple and efficient, has strong interpretation of faults in horizontal slices or vertical sections, and is a good three-dimensional fault automatic tracking scheme.

As can be seen from comparing fig. 9a, 9b, 9c and 9d, the longitudinal continuity of the messy detection result is obviously enhanced, the characterization effect of the sliding fracture is greatly improved, the fault source can be accurately judged, the details of the messy detection result are rich, the fracture is more detailed than the coherent body, and the identification effect of the micro fracture and the crack development is better.

S6, performing large-scale fracture identification based on machine learning on the seismic data subjected to fault enhancement processing in the S4.

With the rapid development of machine learning technology in recent years, convolutional Neural Network (CNN) technology is widely used in seismic reservoir prediction. Theoretically, the technology is very suitable for automatic fracture tracking. The fracture detection algorithm based on CNN image segmentation converts the image classification problem into the image segmentation problem, realizes high-precision identification of fracture by using a Unet-based CNN network, and has the following main characteristics and advantages: ① The training data can be picked up interactively in a three-dimensional space, and various possible conditions capable of simulating fracture development can be generated through a random model, so that full-automatic supervised neural network learning based on big data is realized; ② The most advanced Unet CNN network is used for solving the image segmentation problem; ③ The use of a GPU solves a number of operational problems.

Fig. 10 is a schematic diagram of CNN network structure fracture recognition based on Unet, where the algorithm has high requirement on basic training data, and needs high-precision seed line picking or a large data fracture model established in advance as a support. The calculation process mainly comprises the following two steps: ① The crack development condition of the seed line is picked up interactively in the three-dimensional space, ② generates various possible conditions capable of simulating the crack development through a random model, and the fully automatic supervised neural network learning based on big data is achieved.

As shown in fig. 11 a-11 b and fig. 12 a-12 d, the machine learning technology has better effect on backbone large fracture identification in the application of identifying the sliding fracture of the cone basin, and is different from messy detection, so that the machine learning technology has more reference value on the research on the sliding fracture source. Therefore, by combining the results of messy detection and machine learning fracture identification, the fractures with different scales can be comprehensively depicted.

S7, carrying out fracture zone boundary identification based on the structure tensor on the seismic data subjected to fault enhancement processing in the S4.

A large number of well drilling results show that the broken solution oil reservoir basically develops in a sliding fracture zone, the activity of oil gas outside the fracture zone is obviously reduced, the reservoir forming capability is poor, and even if a good reservoir development exists, the oil gas discovery is difficult to obtain. Therefore, the exploration and development of the disconnected solution oil reservoir are basically carried out around the sliding fracture zone, and the characterization of the fracture zone boundary is particularly important. The early forward attribute simulation shows that the structure tensor has obvious advantages for fracture zone boundary characterization, so that the structure tensor is adopted to identify the fracture zone boundary, and the inside of the fracture zone boundary is the collection of fault, crack, hole, cave and other storage bodies.

The method firstly utilizes the seismic data to calculate the energy gradient vector, then constructs the structure tensor T, finally carries out the eigenvalue decomposition of the structure tensor to obtain three eigenvalues lambda ₁、λ₂、λ₃ and one eigenvector v, and utilizes the obtained three eigenvalues to carry out the geological anomaly identification. The structure tensor can shield the influence of stratum, so that the fracture zone boundary with strong transverse heterogeneity can be effectively identified, and a basis is provided for identification of the fracture zone boundary. The specific calculation process is as follows:

S7-1, gradient calculation:

Calculating a gradient vector (directional derivative) of energy of each point of the seismic data subjected to fault enhancement processing in the step S4:

x, y and z are respectively the line, the channel and the time direction, and u is the seismic amplitude; g ₁ is the differential of the seismic amplitude along the x-direction (line), g ₂ is the differential of the seismic amplitude along the y-direction (xline), g ₃ is the differential in the z or time direction, g is the amplitude gradient vector formed by g ₁、g₂、g₃;

s7-2, construction of Structure tensor T

S7-3, carrying out eigenvalue decomposition on the structure tensor T, wherein the result is three eigenvalues lambda ₁、λ₂、λ₃ and eigenvector v, and the eigenvalue lambda ₂ has a good identification effect on the boundary of the fracture zone through screening and comparison.

As shown in fig. 13 a-13 b, from the actual data identification effect, the structure tensor has obvious advantages in the description of the boundary of the fracture zone, has strong anti-interference capability, has relatively accurate identification effect on the spatial spreading range of the fracture zone, and has relatively good coincidence degree of fracture detection identification results. From the superimposed graph of fracture detection results and the structure tensor, the fracture is substantially within the fracture zone delineated by the structure tensor.

S8, carrying out hole-fracture reservoir identification based on amplitude curvature on the seismic data subjected to fault enhancement processing in the S4.

The disconnected solution reservoirs are generally divided into three types: cracks, holes, caves. In the fracture zone enveloping surface, fracture development and fault and hole development are closely related, and generally fracture and hole development areas are often associated with denser fractures, so that a fracture type reservoir layer is generally reflected in a mess manner on an earthquake, and is often characterized by abnormal properties related to fracture and hole on the earthquake property. Through forward attribute simulation, we know that the amplitude curvature has a good identification effect on the hole-fracture type reservoir, and through the amplitude curvature attribute, we can effectively identify the hole-type reservoir and fracture zones associated with the periphery of the hole-type reservoir.

The amplitude curvature is derived from the amplitude of the seismic data by a transverse second order derivative. First, the first derivative of the main line direction and the crossline direction is calculated by using the seismic amplitude or energy, and the obtained energy gradient attribute can reflect an abnormal geologic body, which is generally called amplitude energy gradient. And then performing second-order derivation on the curved surface to obtain an amplitude curved surface, and finally calculating the curvature attribute of each amplitude according to the curve fitting. In principle, the amplitude energy gradient is a representation of the edge of the geologic body, so that the spatial distribution of the geologic body cannot be obtained by a threshold. The amplitude curvature converts the amplitude energy gradient into an attribute reflecting the envelope of the geologic body, and the spatial distribution of the geologic body can be obtained through a threshold value.

As shown in fig. 14 a-14 b, the amplitude curvature is very sensitive to strong reflection of the beads from the identification effect, the holes and crack development zones around the holes can be well identified, the superposition of the amplitude curvature and the fracture identification result shows that the holes and the cracks have close relation with the fracture development, the favorable reservoir is generally developed in the range of the fracture zone, the amplitude curvature identification result shows that the longitudinal development difference of the controlled reservoir is large, the reservoir mainly develops in shallow layers, and the deep reservoir development gradually becomes worse.

And S9, performing compressed sensing inversion based on the seismic data subjected to fault enhancement processing in the step S4 to identify the cave type reservoir.

Cave type reservoir identification is a difficulty in broken solution reservoir identification, and due to limitation of seismic resolution, the position and the size of the cave type reservoir are difficult to accurately position in seismic attribute research, so that the cave type reservoir is required to be identified by adopting a high-resolution inversion means. From the conclusion of the physical studies of the fractured-fluid reservoir rock, the difference in resistance remains an effective distinguishing parameter for carbonate reservoirs and non-reservoirs, with lower resistance generally indicating better reservoir physical properties. The cavity type reservoir is filled with shale and fluid, has obvious low impedance characteristics, and can be identified by adopting an inversion method.

The invention adopts the advanced compressed sensing high resolution inversion method in industry to identify the cave type reservoir, solves the sparse inversion problem by using the L1 norm (L2 norm is used for conventional inversion) by using the compressed sensing algorithm, and assumes that the stratum reflection coefficient can be represented by odd dipole decomposition. The algorithm is realized by introducing the wavelet matrix into wedge dictionary deconvolution, so that the accuracy and resolution of the inversion result are improved. The objective function of the inversion is:

Min[||s-Wr||₂+λ||r||_p]

Where s is the seismic data, W is the wavelet, r is the reflection coefficient, λ is the regularization coefficient, and p is the norm. In compressed sensing inversion p=1.

The compressive sensing inversion uses a compressive sensing theory in the algorithm implementation process, and the obtained wave impedance is higher than the deterministic inversion resolution of the traditional commercial software; meanwhile, the noise immunity and the stability are better, and the elastic parameters obtained by inversion are more accurate.

15 A-15 d and 15 e-15 h, by designing a single cave and two adjacent cave forward models distributed longitudinally and performing inversion tests, the visible wave impedance inversion has better recognition degree on the cave type reservoir, the compressed sensing inversion result has higher signal-to-noise ratio and obvious resolution improvement, the cave positions are more focused, and the cave positions can be recognized more accurately.

The transverse heterogeneity of the solution reservoir is very strong, and the obtained wave impedance inversion section has certain continuity in the transverse direction, and the heterogeneity characteristic is not obvious enough, which is also determined by the characteristics of the seismic data. In order to highlight the heterogeneity of a broken solution reservoir and eliminate the influence of transverse continuous seismic reflection, the invention carries out residual impedance processing on the wave impedance inversion result. The method is that median filtering is carried out on the wave impedance inversion result to obtain a relatively smooth wave impedance body, and then the wave impedance body is subtracted from the original wave impedance body to obtain the residual impedance.

Comparing fig. 16 a-16 c, it can be seen that the structure tensor can describe the fracture contour range, the amplitude curvature can describe the strong reflection of the beads, the hole and the crack development around the hole can be well identified, the resolution is obviously improved through the compressed sensing high resolution inversion treatment, the cave type reservoir with the conventional attribute difficult to identify is clearly visible on the inversion section, and the vent leakage anastomosis degree is high.

And S10, fusing the results of the steps S5-S9 to finish the three-dimensional depiction of the deep fracture system.

The characterization method of the invention is characterized in that the details of the fracture control reservoir body are characterized in steps from the structure tensor to the amplitude curvature and then to the residual impedance, the structure tensor attribute is characterized by the boundary of fracture zones, all reservoirs such as caves, holes, cracks and the like are included, the amplitude curvature is included by the caves and holes and cracks associated with the periphery of the caves, the residual impedance is the representation of the cavern type reservoir, and the fracture control reservoir body can be accurately characterized by comprehensively applying the three technologies (as shown in fig. 17, the fracture zone contour of the fracture zone characterization by the tensor of the superimposed structure, the fracture-crack of the multi-scale fracture detection identification (including the medium-small scale fracture identification based on the amplitude gradient vector messy detection and the large-scale fracture identification based on the machine learning), the hole and crack development zone of the amplitude curvature characterization and the cavity identification by the compression perception inversion).

The space carving of the broken solution is mainly developed on the basis of comprehensive prediction of the broken solution, and the boundary of the broken solution, messy recognition fracture, amplitude curvature recognition hole-fracture type reservoir and high-resolution inversion recognition hole are recognized by utilizing the structure tensor attribute. And determining threshold values of different attributes by well earthquake calibration, carrying out normalization processing on the different earthquake attributes, and then carrying out spatial attribute fusion to complete the spatial carving of the broken solution.

In the process of attribute fusion, threshold values are determined for different attributes, and data such as drilling curves, emptying leakage points and the like are mainly referred. Statistics show that the drilling curve is sensitive to broken solution boundaries, and the drilling is faster when drilling into broken solution, so the drilling curve can be used to determine the threshold value of the attribute. When a fracture, a hole, a cave or a crack dense belt is drilled, leakage, emptying and the like are generally generated, so that the information can be used for determining threshold values of different properties. The actual calibration result shows that the predicted result and the drilling calibration result have higher correlation, the coincidence rate exceeds 80%, and the predicted result is more reliable.

And 3, carrying out three-dimensional visual perspective analysis on the fused broken solution attributes, and clearly finding out the spatial spreading characteristics (shown in fig. 18a and 18 b) of the broken solution reservoir, wherein the broken solution reservoir has obvious development characteristics along the deep fracture, and the larger the breaking activity strength is, the larger the development scale of the broken solution is.

Based on three-dimensional carving, a three-dimensional attribute detection technology (bodycheck) is adopted, so that the spatial distribution rule of the broken solution oil and gas reservoirs, namely the spatial distribution of different fracture and hole units, can be obtained. The fracture-cavity units are the fracture-control reservoir units with the same oil-water system, the same fracture-cavity units are communicated, and different fracture-cavity units are not communicated. Different fracture-cavity units can be divided according to the drilling production dynamic data by calibrating the attribute threshold value of the drilling. From the dividing result of the fracture-cavity unit, the control effect of the sliding fracture is obvious, the sliding fracture has obvious sectionality on a plane, and the fracture development modes of different parts have large differences, so that the space development and oil gas enrichment of the fracture-controlled reservoir body are determined (as shown in fig. 19a and 19 b).

It should be understood that the foregoing examples of the present invention are provided merely for clearly illustrating the present invention and are not intended to limit the embodiments of the present invention, and that various other changes and modifications may be made therein by one skilled in the art without departing from the spirit and scope of the present invention as defined by the appended claims.

Claims (7)

1. The space depiction method of the deep fracture system comprises the following steps of:

s1, acquiring well drilling data, well logging data and seismic data of a target area;

s2, performing model forward modeling aiming at a sliding fracture-cut-off control reservoir development mode;

s3, screening out sensitive attributes in fracture and reservoir identification through attribute simulation of a model forward result; the sensitive attribute comprises disorder detection, amplitude curvature and structure tensor based on an amplitude gradient vector;

S4, performing fault enhancement processing based on an anisotropic diffusion theory on the seismic data acquired in the S1;

S5, detecting the seismic data subjected to fault enhancement in the S4 based on amplitude gradient vector messiness so as to identify medium-small scale fracture; the method specifically comprises the following steps: s5-1, determining a search azimuth and an inclination angle range according to geological characteristics and data characteristics; s5-2, calculating an amplitude construction tensor matrix in each search direction, and obtaining disorder degree by solving characteristic values of the tensor matrix; s5-3, determining the azimuth and the inclination corresponding to the maximum mess degree as the azimuth and the inclination of the section, wherein the mess degree value is used as the possibility of fracture;

s6, carrying out large-scale fracture identification based on machine learning on the seismic data subjected to fault enhancement processing in the S4;

s7, carrying out fracture zone boundary identification based on the structure tensor on the seismic data subjected to fault enhancement processing in the S4;

S8, carrying out hole-crack reservoir identification based on amplitude curvature on the seismic data subjected to fault enhancement in the S4;

S9, performing compressed sensing inversion based on the seismic data subjected to fault enhancement processing in the S4, and performing residual impedance processing on a wave impedance inversion result of the compressed sensing inversion to identify a cave type reservoir; wherein the inverted objective function is: min [. S-Wr ] ₂+λ∥r∥_p, s is seismic data, W is wavelet, r is reflection coefficient, λ is regularization coefficient, p=1; wherein the residual impedance processing includes: performing median filtering on the wave impedance inversion result to obtain a relatively smooth wave impedance body, and then subtracting the relatively smooth wave impedance body from the original wave impedance body to obtain residual impedance;

s10, fusing the results of the steps S5-S9 to finish three-dimensional depiction of a deep fracture system;

The process of fusing in S10 includes: determining threshold values of different attributes by using well shock calibration and drilling time curves, wherein the threshold values comprise structure tensor, messiness detection, machine learning fracture detection, amplitude curvature and compressed sensing inversion results, and then carrying out normalization processing on different seismic attributes to carry out spatial attribute fusion; performing three-dimensional visual perspective analysis on the fused broken solution attributes to complete three-dimensional depiction of a deep fracture system; and then, a three-dimensional attribute detection technology is adopted to obtain the spatial distribution rule of the disconnected solution oil and gas reservoir.

2. The method for spatially characterizing a deep fracture system according to claim 1, wherein S2 specifically comprises:

Designing forward models of earthquakes with different sizes and different filling media on the basis of fine interpretation of the earthquakes by using drilling and logging data, performing forward simulation by using a wave equation forward method by referring to actual field three-dimensional observation parameters by forward modeling, and obtaining forward earthquake data after migration processing; performing attribute simulation calculation on the seismic data obtained by forward modeling, establishing a quantitative interpretation quantity edition, and correcting the size of an effective reservoir.

3. The method for spatially characterizing a deep fracture system according to claim 2, wherein S3 specifically comprises:

performing pseudo three-dimensional processing on the forward-modeling obtained seismic data to form a pseudo three-dimensional seismic data body, and performing attribute calculation on the pseudo three-dimensional seismic data body; sensitive properties in fracture and reservoir identification are screened out through simulation and selection of various properties of the seismic forward model results.

4. The method for spatially characterizing a deep fracture system according to claim 1, wherein S4 specifically comprises:

s4-1, construction azimuth analysis: determining the azimuth of the reflection phase axis through the earthquake inclination angle and the azimuth angle;

S4-2, breakpoint detection: calculating a coherent body and determining the breakpoint position;

S4-3, unfolding principal component analysis and preserving smooth filtering of boundaries under the constraint of the construction and breakpoint detection results, and completing fault enhancement.

5. The deep fracture system spatial characterization method of claim 1, wherein the machine learning in S6 employs a CNN image segmentation based fracture detection algorithm that implements fracture identification based on Unet CNN networks.

6. The method for spatially characterizing a deep fracture system according to claim 1, wherein S7 specifically comprises:

s7-1, gradient calculation: calculating gradient vector of energy of each point of seismic data subjected to fault enhancement processing in S4

X, y and z are respectively the line, the channel and the time direction, and u is the seismic amplitude; g ₁ is the differential of the seismic amplitude along the x-direction, g ₂ is the differential of the seismic amplitude along the y-direction, g ₃ is the differential in the z-or time-direction, g is the amplitude gradient vector formed by g ₁、g₂、g₃;

s7-2, construction of Structure tensor T

S7-3, carrying out eigenvalue decomposition on the structure tensor T, and recognizing a fracture boundary by adopting the eigenvalue lambda ₂, wherein the result is three eigenvalues lambda ₁、λ₂、λ₃ and eigenvector v.

7. The method for spatially characterizing a deep fracture system according to claim 1, wherein S8 specifically comprises:

s8-1, calculating first derivatives of the main line direction and the crossline direction by using the seismic amplitude or energy to obtain an amplitude energy gradient;

s8-2, performing second order derivation on the amplitude energy gradient to obtain an amplitude curved surface;

s8-3, calculating curvature attributes of each amplitude according to the amplitude curved surface fitting so as to identify the holes and crack development zones around the holes.