A multi-scale fracture detection method based on high-precision spectral decomposition

Abstract

The application discloses a multiscale fracture detection method based on high-precision spectral decomposition, comprising the following steps: S10, inputting a seismic data cube and performing super-resolution spectral decomposition; S20, performing adaptive multiscale frequency band division; S30, performing PDV feature extraction on each frequency band according to the divided multiscale frequency bands, and completing multiscale PDV feature body construction; S40, performing fusion of features of each frequency band through double-channel 3DCNN; and S50, obtaining a fracture probability body through spatial attention weighting. The application breaks through the resolution limit and improves the noise immunity through high-precision time-frequency spectral decomposition, multiscale deep feature fusion and phase-enhanced fracture response.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of oil and gas reservoir exploration, and particularly relates to a multi-scale fracture detection method based on high-precision spectral decomposition. BACKGROUND

[0002] With the deepening of unconventional oil and gas reservoir exploration, the accurate identification of micro-fracture systems (fracture distance <10m) has become a core problem of reservoir evaluation. The traditional seismic fracture detection technology faces three bottlenecks:

[0003] Resolution limit of spectral decomposition method: Theoretical level: Conventional short-time Fourier transform (STFT) and continuous wavelet transform (CWT) are restricted by the Heisenberg uncertainty principle, and there is an inherent contradiction in time-frequency resolution. When the target layer thickness is lower than the tuning thickness (λ / 4, usually 15-30ms), the time-frequency spectrum energy is severely diffused (according to the SEG 2023 report, the frequency band of the thin layer interference zone is widened by more than 45%). Industrial practice defects: Commercial software (such as Petrel TM The signal-to-noise ratio of the spectral decomposition module is deteriorated sharply (SNR<1.5dB) at high frequencies above 45Hz, resulting in micro-fracture responses being overwhelmed by noise.

[0004] Fragmentation problem of multi-scale analysis: Scale fragmentation: Although existing methods (such as the multi-scale coherence cube proposed by Wang et al. 2021) use multiple frequency band analysis, the scale results are only fused by linear weighting, and no frequency band correlation model is established. The measured data of a shale gas block showed that the detection error of the crack system with changing direction was 32%. Adaptive loss: Fixed scale division strategy (such as 10 / 20 / 40Hz three frequency bands) is difficult to adapt to complex geological background. In the igneous rock intrusion area, high-frequency components are strongly absorbed, resulting in a crack detection missing rate of more than 60%.

[0005] Insufficient physical characterization of fracture response characteristics: Amplitude dependence: The mainstream coherence algorithm (C1 / C3) and curvature attribute are based on amplitude discontinuity. However, in the area where gypsum rock develops, the amplitude mutation caused by lithology change is confused with the fracture response, and the proportion of false anomalies is more than 50%. Insufficient use of phase field: Although the phase derivative variance (PDV) theory has been proposed, its application in the multi-scale domain is still limited to a single frequency band. A case study of a carbonate oil field showed that the boundary positioning error of the conventional PDV method for dissolved fractures was ±15m. SUMMARY

[0006] To solve the above problems, the application provides a multi-scale fracture detection method based on high-precision spectral decomposition, which breaks through the resolution limit and improves the noise immunity through high-precision time-frequency spectrum decomposition, multi-scale depth feature fusion and phase-enhanced fracture response.

[0007] To achieve the above object, the technical scheme adopted by the present application is: a multi-scale fracture detection method based on high-precision spectral decomposition, comprising the steps of:

[0008] S10, inputting a seismic data cube, performing super-resolution spectral decomposition;

[0009] S20, performing adaptive multi-scale frequency band division;

[0010] S30, according to the divided multi-scale frequency bands, respectively extracting PDV features of each frequency band, and completing multi-scale PDV feature body construction;

[0011] S40, dual-channel 3DCNN fusion, fusing the features of each frequency band;

[0012] S50, obtaining a fracture probability body through spatial attention weighting.

[0013] Further, the super-resolution spectral decomposition comprises the steps of:

[0014] S101: initialization of time-frequency representation of seismic signals;

[0015] input preprocessing: reading a post-stack seismic data body, performing amplitude normalization processing on each gather signal, and eliminating energy deviation caused by differences in acquisition equipment; basic time-frequency spectrum generation: using a complex Morlet wavelet as a mother wavelet, performing continuous wavelet transform CWT, discretizing in the scale space, calculating wavelet coefficients, and generating an initial time-frequency spectrum;

[0016] S102: accurate reconstruction of instantaneous frequency field:

[0017] phase derivative calculation: calculating the time direction derivative of the phase component of the wavelet coefficient to obtain an instantaneous angular frequency estimate, which re-maps the spread energy to the real physical frequency coordinate;

[0018] energy redistribution mechanism: compressing the wavelet coefficient along the frequency axis to the instantaneous frequency position;

[0019] S103: time-frequency spectrum sharpening and noise suppression:

[0020] Gaussian focusing enhancement: applying an adaptive Gaussian window on the synchronized compressed time-frequency spectrum, a narrow window in the high-frequency region enhances the transient feature, and a wide window in the low-frequency region suppresses the aliasing;

[0021] noise adaptive filtering: using median filtering in the time-frequency domain: calculating the energy median in the local time window, setting the energy below 20% of the median to zero, and retaining the burst high-frequency component caused by the fracture;

[0022] S104: time-frequency data body reconstruction:

[0023] Four-dimensional data volume construction: for each sample point of seismic data volume, calculate its synchronous compressive time-frequency spectrum, reorganize into four-dimensional tensor, get super-resolution time-frequency spectrum cube.

[0024] Further, the adaptive multi-scale frequency band division includes the steps of:

[0025] S201: three-dimensional spectrum entropy field construction; input data: super-resolution time-frequency spectrum cube; entropy value calculation: for each spatial point (x, y) of time-frequency spectrum T s (t, ω), calculate time-frequency information entropy H:

[0026]

[0027] wherein ω is a frequency variable, P(ω) is a normalized power spectrum probability density, ω min ~ ω max is a value range of the frequency variable, and t is a time variable;

[0028] Geological significance mapping: low entropy area (H < 2.0): effective signal, high entropy area (H > 3.5): noise or geologically irrelevant signal;

[0029] S202: scale space generation guided by geological targets:

[0030] Scale parameter definition: microscopic scale: scale factor s = 10, target geological body is microcrack, frequency coverage range is 60-120Hz; mesoscopic scale: scale factor s = 20, target geological body is small and medium-sized fault, frequency coverage range is 30-60Hz; macroscopic scale: scale factor s = 40, target geological body is large fault zone, frequency coverage range is 10-30Hz;

[0031] Dynamic scale activation: if the maximum entropy value < 3.0 and the dominant frequency > 45Hz, the microscopic scale is activated; if the stratum dip angle > 30°, the mesoscopic scale is forcibly activated; if there is igneous rock intrusion body, the microscopic scale is disabled;

[0032] S203: dominant frequency band extraction:

[0033] Peak energy search: in the low entropy area (H < 2.5), search for the local energy maximum value point ω peak along the frequency axis, which satisfies:

[0034] |T s (ω peak )|>1.5×median(|T s (ω)|);

[0035] wherein T s (ωpeak) is the synchronous compressive time-frequency spectrum at ω peakThe amplitude of the band, the median is the median of the band energy;

[0036] Geophysical constraint frequency width: the band width Δω is dynamically adjusted according to the formation quality factor Q:

[0037]

[0038] Where, α is the lithology coefficient, Vp / Vs is the ratio of longitudinal and transverse wave velocity.

[0039] Further, the multi-scale PDV feature body construction includes:

[0040] Input data: adaptive division of frequency band;

[0041] Phase derivative variance PDV calculation: for each scale of time-frequency spectrum data:

[0042] (1) Extract the instantaneous phase field: φ ω (x, y, t);

[0043] (2) Calculate the amplitude of the phase gradient:

[0044] (3) Calculate the variance of the gradient amplitude in the local window to generate the PDV feature cube.

[0045] At the fracture boundary Sharp increase, representing a sharp increase in PDV value; weak PDV response in lithology change area.

[0046] Further, the dual-channel 3DCNN architecture includes:

[0047] Channel 1, input feature is high-frequency PDV, network module is densely connected convolution block, physical constraint mechanism is convolution kernel size ≤3 3 , target detection object is micro crack;

[0048] Channel 2, input feature is low-frequency PDV, network module is residual dilated convolution block, physical constraint mechanism is dilated rate = 2, target detection object is large fault;

[0049] The feature fusion method adopts cross-scale feature splicing: the dual-channel outputs are spliced in the feature dimension to form a multi-scale joint feature body, channel 1 micro feature locates fracture details, and channel 2 macro feature maintains structural continuity; Physical constraint rule: when the spatial distance of dual-channel fracture labels is <10m, forced fusion is carried out to avoid fracture fragmentation.

[0050] Further, the spatial attention enhancement includes:

[0051] Input dual-channel fusion feature body F cat , point multiplication attention weight A to get enhanced feature Faat：

[0052] F att = A o F cat ;

[0053] where A is the attention weight, used to enhance the fracture edge and suppress the homogeneous background.

[0054] Further, the attention weight calculation formula is:

[0055] A(x, y, t) = 1 / (1 + e^[-k-(G(x, y, t)-μ)])

[0056] k is the response sharpening factor; μ is the gradient activation threshold, and G(x, y, t) is the F cat gradient enhanced change value.

[0057] Further, the F cat gradient enhanced change value calculation formula is:

[0058]

[0059] where arctan is the inverse tangent function, is the energy body space gradient, and γ is the gradient sensitivity coefficient.

[0060] Further, the double-channel fusion feature body Fcat is compressed, the L2 norm is calculated along the channel dimension, the saliency energy body E(x, y, t) = ||Fcat||2 is generated, and the high energy area corresponding to the fracture or stratigraphic boundary is obtained.

[0061] Further, the fracture probability is predicted, and the enhanced feature Fatt obtained through spatial attention enhancement is compressed to a single channel through a 1x1 convolution to obtain the fracture probability body P f ;

[0062] The geological rationality constraints include: vertical continuity: 3D conditional random field (CRF) constraint is adopted to ensure that the fracture extends in the time direction; dip angle control: when the predicted fracture dip angle > 75°, it needs to satisfy to avoid the straight-up false appearance;

[0063] The loss function is:

[0064]

[0065] where BCE represents the main loss, Y label is the fracture label explained by experts, N is the total number of training samples, i is the ith training sample, L edge is the edge loss term, is the direction loss term, and n is the stratigraphic normal vector, For the vertical gradient of the fracture probability, lambda1 and lambda2 are weight coefficients.

[0066] The beneficial effects of the technical solution are:

[0067] The present application realizes resolution limit reconstruction, realizes 1.5-meter-level micro-fault accurate detection (traditional method >=8 meters), and solves the problem that "seismic cannot identify fault distance less than tuning thickness".

[0068] The present application realizes resolution limit reconstruction, realizes 1.5-meter-level micro-fault accurate detection (traditional method >=8 meters), and solves the problem that "seismic cannot identify fault distance less than tuning thickness".

[0069] The present application realizes resolution limit reconstruction, realizes 1.5-meter-level micro-fault accurate detection (traditional method >=8 meters), and solves the problem that "seismic cannot identify fault distance less than tuning thickness".

[0070] The present application realizes resolution limit reconstruction, realizes 1.5-meter-level micro-fault accurate detection (traditional method >=8 meters), and solves the problem that "seismic cannot identify fault distance less than tuning thickness".

[0071] The present application realizes resolution limit reconstruction, realizes 1.5-meter-level micro-fault accurate detection (traditional method >=8 meters), and solves the problem that "seismic cannot identify fault distance less than tuning thickness". BRIEF DESCRIPTION OF DRAWINGS

[0072] Figure 1 A multi-scale fault detection method based on high-precision spectrum decomposition of the present application is shown in the flowchart;

[0073] Figure 2 A three-dimensional convolution network diagram in the embodiment of the present application is shown in the flowchart;

[0074] Figure 3 A three-dimensional convolution network diagram in the embodiment of the present application is shown in the flowchart;

[0075] Figure 4 A three-dimensional convolution network diagram in the embodiment of the present application is shown in the flowchart; DETAILED DESCRIPTION

[0076] In order to make the purpose, technical scheme and advantages of the present application clearer, the present application will be further described below with reference to the drawings.

[0077] The present application proposes a multi-scale fault detection method based on high-precision spectrum decomposition, as shown in Figure 1 The steps include:

[0078] S10, input the seismic data cube, and perform super-resolution spectrum decomposition;

[0079] S20, performing adaptive multi-scale frequency band division;

[0080] S30, according to the divided multi-scale frequency bands, respectively extracting PDV features of each frequency band, and completing multi-scale PDV feature body construction;

[0081] S40, dual-channel 3DCNN fusion, fusing the features of each frequency band;

[0082] S50, obtaining the fracture probability body through spatial attention weighting.

[0083] As an optimization scheme of the above embodiment, the super-resolution spectral decomposition comprises the steps of:

[0084] S101: initialization of time-frequency representation of seismic signal;

[0085] Input preprocessing: reading the post-stack seismic data body, performing amplitude normalization processing on each gather signal, and eliminating the energy deviation caused by the difference of acquisition equipment; basic time-frequency spectrum generation: using complex Morlet wavelet as mother wavelet, performing continuous wavelet transform CWT, discretizing in scale space, calculating wavelet coefficients, and generating initial time-frequency spectrum;

[0086] S102: accurate reconstruction of instantaneous frequency field:

[0087] Phase derivative calculation: calculating the time direction partial derivative of the phase component of the wavelet coefficient to obtain the instantaneous angular frequency estimation, which remaps the diffused energy to the real physical frequency coordinate;

[0088] Energy redistribution mechanism: compressing the wavelet coefficient along the frequency axis to the instantaneous frequency position;

[0089] S103: time-frequency spectrum sharpening and noise suppression:

[0090] Gaussian focusing enhancement: applying an adaptive Gaussian window on the synchronous compressed time-frequency spectrum, a narrow window in the high-frequency region enhances the transient feature, and a wide window in the low-frequency region suppresses the aliasing;

[0091] Noise adaptive filtering: using median filtering in time-frequency domain: calculating the energy median in the local time window, setting the energy below 20% of the median to zero, and retaining the burst high-frequency component caused by fracture;

[0092] S104: time-frequency data body reconstruction:

[0093] Four-dimensional data body construction: calculating the synchronous compressed time-frequency spectrum of each sampling point of the seismic data body, reorganizing it into a four-dimensional tensor, and obtaining a super-resolution time-frequency spectrum cube;

[0094] Physical mechanism breakthrough: through the instantaneous frequency field reconstruction, the dispersion energy is refocused to the real frequency coordinates, solving the problem of spectrum aliasing in thin layer interference zone

[0095] Geological adaptability design: time-varying Gaussian window adapts to stratum absorption and attenuation, and entropy-driven frequency band extraction avoids noise interference.

[0096] Industrialization implementation: parallel computing framework makes the processing time of 100GB data body less than 4 hours, and the output four-dimensional time-frequency body can be directly input into the multi-scale feature extraction module.

[0097] This technology lays a resolution foundation for micro-fracture detection, and its core is to deeply couple the mathematical frequency redistribution theory and the physical law of seismic wave propagation, forming a high-precision spectral decomposition method with geological interpretability.

[0098] As shown in Table 1, the above-mentioned synchronous compression change method of the present application is superior to the traditional STFT method.

[0099] Table 1

[0100] Parameter Conventional STFT Synchronous compressive transform Thin layer resolution ≥ 15 ms ≤ 5 ms Window width at 90 Hz 12 ms 3 ms Break response SNR 1.8 dB 6.5 dB Computational complexity O(N log N) O(N 2 )→GPU parallelization

[0101] As an optimization scheme of the above-mentioned embodiment, the information entropy theory and the geological physical constraint are fused, and the traditional fixed frequency band division is broken through, and the adaptive multi-scale frequency band division includes the following steps:

[0102] S201: three-dimensional spectral entropy field construction; input data: super-resolution time-frequency spectrum cube (four-dimensional tensor: space x time x frequency); entropy value calculation: for each spatial point (x, y) of the time-frequency spectrum T s (t,ω), the time-frequency information entropy H is calculated:

[0103]

[0104] Wherein, ω is a frequency variable, P(ω) is a normalized power spectrum probability density, ω min ~ω max is the value range of the frequency variable, and t is a time variable;

[0105] Geological significance mapping: low entropy area (H<2.0): effective signal (such as fracture response, stratum interface), high entropy area (H>3.5): noise or geologically irrelevant signal (such as multiple wave);

[0106] S202: scale space generation guided by geological targets:

[0107] Scale parameter definition: Micro scale: scale factor s = 10, target geological body is micro crack (<5m throw), frequency coverage range is 60-120Hz; meso scale: scale factor s = 20, target geological body is small and medium-sized fault (5-20m), frequency coverage range is 30-60Hz; macro scale: scale factor s = 40, target geological body is large fault zone (>20m), frequency coverage range is 10-30Hz;

[0108] Dynamic scale activation: if the maximum entropy value <3.0 and the main frequency >45Hz, the micro scale is activated; if the stratum dip angle >30°, the meso scale is activated forcibly; if there is igneous rock intrusion, the micro scale is disabled;

[0109] S203: dominant frequency band extraction:

[0110] Peak energy search: in the low entropy area (H <2.5), the local energy maximum point ω is searched along the frequency axis peak , which satisfies:

[0111] |T s (ω peak )|>1.5×median(|T s (ω)|);

[0112] Wherein, T s (ωpeak) is the amplitude of the spectrum at ω peak , and median is the median of the band energy;

[0113] Geophysical constraint frequency width: the frequency band width Δω is dynamically adjusted according to the stratum quality factor Q:

[0114]

[0115] Wherein, α is the lithology coefficient, Vp / Vs is the ratio of longitudinal wave velocity to transverse wave velocity

[0116] Shale: Q = 8-12 Q = 8-12→ narrow frequency band (focus on thin layer);

[0117] Sandstone: Q = 15-20→ wide frequency band (adapt to rapid phase change).

[0118] The information entropy theory, rock physics law and signal significance analysis are deeply coupled in the application, the frequency band division of "geological scene self-perception" is realized, and the industry pain point that the fixed scale model cannot adapt to complex geological conditions is fundamentally solved. As shown in Table 2, the advantages of the adaptive method proposed in the application compared with the traditional fixed frequency band division method.

[0119] Table 2

[0120]

[0121] As an optimization scheme of the above embodiment, the multi-scale PDV feature body construction comprises:

[0122] Input data: adaptive divided frequency bands; (micro / meso / macro scale corresponding frequency band range)

[0123] Phase derivative variance PDV calculation:

[0124] For each scale of time-frequency spectrum data:

[0125] (1) Extract the instantaneous phase field: φ ω (x, y, t);

[0126] (2) Calculate the phase gradient amplitude:

[0127] (3) Calculate the variance of the gradient amplitude in the local window to generate the PDV feature cube.

[0128] At the fracture boundary The sharp increase represents a sudden increase in PDV value; the PDV response in the lithology change area is weak.

[0129] As an optimization scheme of the above embodiment, as shown in Figure 2 , the dual-channel 3DCNN architecture comprises:

[0130] Channel 1, the input feature is high-frequency PDV (60-120Hz), the network module is a densely connected convolution block, the physical constraint mechanism is that the convolution kernel size is less than or equal to 3 3 , and the target detection object is microcrack (less than 5m fracture distance);

[0131] Channel 2, the input feature is low-frequency PDV (10-30Hz), the network module is a residual dilated convolution block, the physical constraint mechanism is that the dilated rate is 2, and the target detection object is a large fault (more than 20m);

[0132] The feature fusion method adopts cross-scale feature splicing: the dual-channel outputs are spliced in the feature dimension to form a multi-scale joint feature body, the micro-scale feature of channel 1 locates the fracture details, and the macro-scale feature of channel 2 maintains the structural continuity; the physical constraint rule is that when the spatial distance of the dual-channel fracture labels is less than 10m, it is forced to be fused to avoid fracture fragmentation.

[0133] As an optimization scheme of the above embodiment, the spatial attention enhancement comprises:

[0134] Input dual-channel fusion feature body F cat , multiply the attention weight A to obtain the enhanced feature F aat：

[0135] F att =A⊙F cat ;

[0136] where A is the attention weight, used to enhance the fracture edge and suppress the homogeneous background, as shown in Figure 3 ;

[0137] Specifically, the attention weight calculation formula is:

[0138] A(x, y, t) = 1 / (1 + e^[-k · (G(x, y, t) - μ)])

[0139] k is the response sharpening factor, and the typical value is 8, the greater the fracture boundary is sharper; μ is the gradient activation threshold, and the typical value is 0.6, G(x, y, t) is the F cat gradient enhanced change value.

[0140] Specifically, the F cat gradient enhanced change value calculation formula is:

[0141]

[0142] where arctan is the inverse tangent function, is the energy body space gradient, γ is the gradient sensitivity coefficient, and the empirical value is 0.05, low γ value → enhance weak edge response (suitable for micro cracks), high γ value → suppress noise gradient (suitable for igneous rock area).

[0143] Specifically, the double-channel fusion feature body Fcat is compressed, and the L2 norm is calculated along the channel dimension to generate the saliency energy body E(x, y, t) = ||Fcat||2, and the high energy area corresponds to the fracture or stratigraphic boundary.

[0144] As an optimization scheme of the above embodiment, the fracture probability prediction compresses the enhanced feature Fatt obtained by spatial attention enhancement to a single channel through 1x1 convolution to obtain the fracture probability body P f ;

[0145] The geological rationality constraints include: vertical continuity: 3D conditional random field (CRF) constraint is adopted to ensure that the fracture extends in the time direction; dip angle control: when the predicted fracture dip angle > 75°, it needs to meet to avoid the upright false impression;

[0146] The loss function is:

[0147]

[0148] where BCE represents the main loss, Y label is the fracture label explained by experts, N is the total number of training samples, i is the i-th training sample, L edge is the edge loss term, For the directional loss term, n is the formation normal vector, For the fracture probability vertical gradient, λ1 and λ2 are weight coefficients, as shown in Figure 4

[0149] Industrial output: fracture probability volume: 0-1 continuous value;

[0150] Binary fracture model:

[0151]

[0152] τ(t): time-varying threshold, deep tolerance increases, such as t>2s τ=0.6.

[0153] The present application realizes the leap of seismic fracture detection from "black box model" to "physical interpretation AI" through the triple technology closed loop of physical feature (PDV) driving + network architecture geological constraint + attention focusing edge, and provides a reliable solution for complex geological target identification.

[0154] As shown in Table 3, the advantages of the present method compared with traditional AI are illustrated.

[0155] Table 3

[0156]

[0157]

[0158] The above shows and describes the basic principles and main features of the present application and the advantages of the present application. It should be understood by those skilled in the art that the present application is not limited by the above examples, and the above examples and descriptions in the specification are only to illustrate the principles of the present application. Without departing from the spirit and scope of the present application, various changes and improvements can be made to the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.​

Claims

1. A multi-scale fracture detection method based on high-precision spectral decomposition, characterized in that, The method comprises the steps of: S10, inputting a seismic data cube, super-resolution spectral decomposition, comprising the steps of: S101: initialization of time-frequency representation of seismic signals; Input preprocessing: reading the post-stack seismic data volume, performing amplitude normalization processing on each gather signal, and eliminating energy deviation caused by acquisition equipment differences; basic time-frequency spectrum generation: using a complex Morlet wavelet as a mother wavelet, performing continuous wavelet transform CWT, discretizing in the scale space, calculating wavelet coefficients, and generating an initial time-frequency spectrum; S102: accurate reconstruction of instantaneous frequency field: Phase derivative calculation: calculating the time direction partial derivative of the phase component of the wavelet coefficient to obtain an instantaneous angular frequency estimate, which remaps the dispersed energy to the real physical frequency coordinate; Energy redistribution mechanism: compressing the wavelet coefficient along the frequency axis to the instantaneous frequency position; S103: time-frequency spectrum sharpening and noise suppression: Gaussian focusing enhancement: applying an adaptive Gaussian window to the synchronously compressed time-frequency spectrum, a narrow window in the high-frequency region enhances transient characteristics, and a wide window in the low-frequency region suppresses aliasing; Noise adaptive filtering: using median filtering in the time-frequency domain: calculating the energy median in the local time window, setting the energy below 20% of the median to zero, and retaining the burst high-frequency component caused by the fracture; S104: time-frequency data volume reconstruction: Four-dimensional data volume construction: for each sampling point of the seismic data volume, calculate its synchronously compressed time-frequency spectrum, reorganize it into a four-dimensional tensor, and obtain a super-resolution time-frequency spectrum cube; S20, adaptive multi-scale frequency band division is performed; S30, according to the divided multi-scale frequency bands, respectively, the PDV feature extraction of each frequency band is carried out, and the multi-scale PDV feature volume construction is completed; S40, double-channel 3DCNN fusion, the features of each frequency band are fused; S50, through spatial attention weighting, a fracture probability volume is obtained.

2. The method of claim 1, wherein, The adaptive multi-scale frequency band division comprises the steps of: S201: Construction of a 3D Spectral Entropy Field; Input Data: Super-resolution time-frequency spectrum cube; Entropy Calculation: Time-frequency spectrum T for each spatial point (x,y) s (t,ω), calculate the time-frequency information entropy H: where ω is the frequency variable, P(ω) is the normalized power spectral probability density, ω min ~ ω max is the range of values of the frequency variable, t is the time variable; Geological significance mapping: low entropy H<2.0: effective signal, high entropy H>3.5: noise or geologically irrelevant signal; S202: scale space generation oriented to geological targets: Scale parameter definition: microscopic scale: scale factor s=10, target geological body is microcrack, frequency coverage range is 60-120Hz; mesoscopic scale: scale factor s=20, target geological body is small fault, frequency coverage range is 30-60Hz; macroscopic scale: scale factor s=40, target geological body is large fault zone, frequency coverage range is 10-30Hz; Dynamic scale activation: if the maximum entropy value is less than 3.0 and the dominant frequency is greater than 45Hz, the microscopic scale is activated; if the stratum dip angle is greater than 30°, the mesoscopic scale is forcibly activated; if there is an igneous rock intrusion body, the microscopic scale is disabled; S203: dominant frequency band extraction: Peak energy search: Within the low-entropy region H < 2.5, search for local energy maxima along the frequency axis ω peak that satisfy: |T s (ω peak )| > 1.5 x median(|T s (ω)|); where T s (ω peak ) is the amplitude of the spectrum at ω peak , median is the median of the band energy. Geophysical constraint frequency width: the frequency band width Δω is dynamically adjusted according to the stratum quality factor Q: where a is the lithology coefficient, V p / V s is the ratio of P-wave to S-wave velocity.

3. The multi-scale fracture detection method based on high-precision spectral decomposition of claim 1 or 2, wherein, The multi-scale PDV feature volume construction comprises: Input data: adaptively divided frequency bands; Phase derivative variance PDV calculation: for time-frequency spectrum data of each scale: (1) Extract the instantaneous phase field: φ ω (x, y, t); (2) Calculate the phase gradient magnitude: (3) calculating the variance of the gradient amplitude in the local window to generate a PDV feature cube; Phase gradient amplitude at break boundary The PDV value increases sharply, which represents a sharp increase; the PDV response in the lithology change area is weak.

4. The method of claim 1, wherein, The double-channel 3DCNN architecture comprises: Channel 1, input feature is high-frequency PDV, network module is dense connection convolution block, physical constraint mechanism is convolution kernel size ≤ 3 3 , target detection object is micro crack; Channel 2, input feature is low-frequency PDV, network module is residual dilated convolution block, physical constraint mechanism is dilated rate = 2, and the target detection object is a large fault; The feature fusion mode adopts cross-scale feature splicing: the outputs of the two channels are spliced in the feature dimension to form a multi-scale joint feature body, the micro features of channel 1 locate the fracture details, and the macro features of channel 2 maintain the structural continuity; the physical constraint rule is: when the spatial distance of the fracture labels of the two channels is less than 10 m, it is forced to be fused to avoid fragmentation of the fracture.

5. The method of claim 1, wherein, The spatial attention enhancement includes: input dual-channel fusion feature body F cat dot product attention weight A to obtain enhanced feature F aat： F att = A O F cat ; Wherein, A is the attention weight, used to strengthen the fracture edge and suppress the homogeneous background.

6. The method of claim 5, wherein, The attention weight calculation formula is: A(x,y,t)=1 / (1+e^[-k·(G(x,y,t)-μ)]) k is a response sharpening factor; μ is a gradient activation threshold, G(x, y, t) is F cat gradient intensification change value.

7. The method of claim 6, wherein, F cat The gradient reinforcement change value calculation formula is: where arctan is the inverse tangent function, is the energy body spatial gradient, and γ is the gradient sensitivity coefficient.

8. The method of claim 7, wherein, To the double channel fusion feature body F cat The feature compression is carried out, the L2 norm is calculated along the channel dimension, and the saliency energy body E(x, y, t) = ||F cat ||2, and the high-energy area corresponds to the fracture or the formation boundary.

9. The multi-scale fracture detection method based on high-precision spectral decomposition of claim 1 or any one of claims 5-8, wherein, Fracture probability prediction, enhanced features F obtained by spatial attention enhancement are compressed to single channel to obtain fracture probability body P att Fracture probability prediction, enhanced features F obtained by spatial attention enhancement are compressed to single channel to obtain fracture probability body P f ; Geological constraints include: vertical continuity: 3D conditional random field (CRF) is used to ensure the continuity of faults in the time direction; dip control: when the predicted fault dip is greater than 75°, the following conditions need to be met to avoid the upright false appearance; The loss function is: where BCE represents the main loss, Y label is the fracture label explained by experts, N is the total number of training samples, i is the ith training sample, L edge is the edge loss term, is the direction loss term, n is the formation normal vector, is the vertical gradient of fracture probability, λ1 and λ2 are weight coefficients.

Images