Shale gas reservoir space detection method and system

Abstract

The invention relates to the technical field of geological exploration, in particular to a shale gas reservoir space detection method and system, and the method comprises the following steps: collecting reservoir differentiation depth pore size and horizon point location data, judging the stability of a layer section according to a single-side average gradient value of the longitudinal change of pores, selecting an abnormal point as a transition zone, and determining the thickness of the transition zone; according to the method, through pore structure and horizon data extraction, fine division of the horizon boundary is realized, the structural differentiation identification capability is improved, physical property aliasing is avoided, wave velocity difference and local stress interference are fused, the stress response identification precision of a microcosmic unit is enhanced, and the reservoir stratification boundary division result is generated. According to the method, the internal structure difference of the reservoir is reflected, the pore axis and the fracture distribution direction are matched, the space coupling behavior is identified, the amplitude attenuation and pressure inversion distribution are combined, the reservoir unit with the stress retention and fluid wrapping characteristics is accurately positioned, and the accuracy and practicability of space identification are improved.

Description

Technical Field

[0001] The present invention relates to the field of geological exploration technology, and in particular to a shale gas reservoir space detection method and system. Background Art

[0002] The field of geological exploration technology includes methods and means for acquiring, analyzing and interpreting underground geological structures and the distribution status of mineral resources. The core content of this technical field is to use geophysics, geochemistry, remote sensing measurement and drilling technology to identify and locate underground resources, and to determine the properties of rock formations, structural characteristics and oil and gas content. Geological exploration is widely used in the preliminary investigation and development deployment of energy resources such as oil and gas, coal, metal minerals, and shale gas. Its technical means mainly include seismic wave reflection, electromagnetic response, resistivity change, ground stress distribution analysis and geochemical anomaly identification, etc., which provide a basis for subsequent resource exploitation by establishing an underground geological model.

[0003] Among them, the shale gas reservoir space detection method refers to the means used to identify and judge the reservoir space distribution and its physical characteristics in the shale layer. The technical matters targeted by the patent subject cover the identification of pore structure, fracture system and gas content in shale reservoirs. Specifically, the structural construction and parameter characterization of the reservoir space are achieved by using multi-component seismic data for wave impedance inversion, fracture anisotropy analysis, and joint inversion of elastic parameters. The method involves the construction of a wave field propagation model using pre-processed seismic waveform data, obtaining anisotropic parameters through inversion to reveal fracture density and direction, and conducting joint interpretation in combination with reservoir lithology constraints to identify spatial areas with reservoir potential.

[0004] Existing technologies often rely on macroscopic geological layers or reflection features for division when processing layers, which makes it difficult to accurately capture the continuous changes in pore structure on a vertical scale, easily leading to blurred boundaries, resulting in unclear attribute fusion and classification. In terms of stress response analysis, parameter settings are often based on a single value of wave velocity, lacking elastic adjustment under the background of local disturbances, resulting in distorted results and failure to reveal the stress difference characteristics between structures. In the trend extraction method, the methods mostly focus on the area where outliers are concentrated, and the analysis of path continuity and directional deviation is insufficient, which weakens the accurate characterization of the main control line of the deformation trend. When identifying the relationship between pores and cracks, static attribute overlap is often used to determine the distribution of contact areas, and the actual synergistic response unit cannot be identified based on the directional coupling state, which easily misjudges invalid space as a potential block. In the identification of reservoir space, attribute selection focuses on the response of a single indicator, ignoring the dual effects of stress and fluid on the reservoir capacity, resulting in incomplete identification of gas content. For example, in areas with complex fractures and pore structures interlaced, traditional analysis models are difficult to distinguish areas where synergistic behavior appears, resulting in resource deployment errors. Summary of the invention

[0005] To solve the technical problems existing in the prior art, an embodiment of the present invention provides a method for detecting shale gas reservoir spaces, including the following steps:

[0006] To achieve the above object, the present invention adopts the following technical solution: A method for detecting shale gas reservoir spaces, including the following steps: S1: Collect data on pore sizes and layer position points at different depths of the reservoir. Judge the stability of the layer segment according to the unilateral average gradient value of the longitudinal change of the pores. Select the abnormal points as the transition zone. Divide the physical property independent layer segments according to the position, and generate the result of the division of the reservoir layer boundaries; S2: Based on the result of the division of the reservoir layer boundaries, calculate the difference between the original wave velocity value and the pore mean value in the layer segment. Adjust the wave velocity response intensity according to the difference and the local stress interference data to form a set of structure response difference indexes in the reservoir; S3: Invoke the set of structure response difference indexes. Take the index offset direction in the region as the reference, delimit the change trend, set the included angle threshold, extract the group of graphics of the direction deviation interval, arrange the paths according to the gradient consistency, and extract the main control path band of the reservoir deformation trend; S4: Invoke the main control path band of the reservoir deformation trend. Select the pore-fracture contact data in the area where the path passes through. Screen the contact point group according to the included angle comparison to obtain the block set where the pore-fracture cooperative response appears; S5: Based on the block set where the pore-fracture cooperative response appears, combine the combination distribution of the wave amplitude attenuation rate and the pressure inversion value to screen the space units with stress retention and fluid inclusion characteristics to obtain the group of maps for delimiting the reservoir space blocks.

[0007] As a further solution of the present invention, the result of the division of the reservoir layer boundaries includes the distribution of the interlayer transition zone, the type of physical property independent layer segments, the level of the stability of the layer segment structure, and the level of the accuracy of the boundary division. The set of structure response difference indexes in the reservoir includes the distribution of the difference between the wave velocity and the pores, the intensity of local stress disturbance, the offset amount of stress deformation amplitude, and the wave velocity response correction coefficient. The main control path band of the reservoir deformation trend includes the consistency of the index offset direction, the continuity of the trend path, the direction included angle deviation area, and the gradient consistency section. The block where the pore-fracture cooperative response appears includes the distribution of pore-fracture contact points, the difference interval of the direction included angle, the intensity of the cooperative response appearance, and the matching relationship between the pore-fracture structure. The group of maps for delimiting the reservoir space blocks includes the wave amplitude attenuation rate interval, the combination distribution of the pressure inversion value, the identification result of the stress retention area, and the fluid inclusion characteristic unit.

[0008] As a further solution of the present invention, the specific steps for collecting data on pore sizes and layer position points at different depths of the reservoir, judging the stability of the layer segment according to the unilateral average gradient value of the longitudinal change of the pores, selecting the abnormal points as the transition zone, dividing the physical property independent layer segments according to the position, and generating the result of the division of the reservoir layer boundaries are as follows: S101: Collect the data of the pore sizes and the formation position points at different depths in the reservoir, calculate the pore size gradient values corresponding to the formation position points, calculate the unilateral average gradient value based on the pore size change in the formation, extract the unilateral average gradient values of the formation segments to form a distribution sequence, and calculate the gradient change rate of the formation segments by calculating the gradient values of adjacent formations in the distribution sequence; S102: Invoke the gradient change rate of the formation segments, screen the positions with abnormal mutations, judge the continuous abnormal points at the positions, calculate the distribution intervals of the abnormal points, delimit the range of the interlayer transition zone according to the distribution intervals of the abnormal points, and generate the position distribution of the interlayer transition zone; S103: Invoke the position distribution of the interlayer transition zone, divide the reservoir into physical property independent formation segments according to the position of the transition zone in the formation segments, adjust the boundaries of the independent formation segments, and generate the result of the reservoir layer division boundary;

[0009] As a further solution of the present invention, based on the result of the reservoir layer division boundary, the specific steps of calculating the difference between the original wave velocity value and the pore mean value in the formation segment, adjusting the wave velocity response intensity according to the difference and the local stress interference data, and forming the group of reservoir internal structure response difference indexes are as follows: S201: Based on the result of the reservoir layer division boundary, obtain the original wave velocity value in the formation segment, extract the mean value of the porosity, calculate the difference between the original value and the mean value of the porosity, compare the original wave velocity value and the pore mean value in the formation segment, obtain the difference data of each formation segment, and record the difference as the formation segment difference data; S202: Based on the formation segment difference data, monitor and collect the local stress interference data, compare the difference of each formation segment with the corresponding stress interference data, adjust the wave velocity response intensity, obtain the wave velocity correction value of the formation segment under the stress interference, and generate the wave velocity correction intensity data by comparing the differential correction values; S203: Based on the wave velocity correction intensity data, compare the stress deformation amplitudes before and after correction, calculate the stress deformation offset of the formation segment, obtain the difference value of the internal structure response of the formation segment, and form the group of reservoir internal structure response difference indexes according to the difference value.

[0010] As a further solution of the present invention, the specific calculation formula of the stress deformation offset of the formation segment is: ; wherein, represents the stress deformation offset of the th formation segment, represents the original wave velocity value before wave velocity correction of the th sampling point in the th formation segment, represents the wave velocity value after wave velocity correction of the th sampling point in the th formation segment, represents the local stress disturbance value of the th sampling point in the th layer segment, represents the stress disturbance adjustment coefficient of the th layer segment, and its value is determined by linear regression based on the average disturbance level of the previous measurement points, represents the average value of all original wave velocity sampling points in the th layer segment, represents the corrected balance factor of the th layer segment.

[0011] As a further solution of the present invention, by calling the group of structural response difference indicators, taking the index offset direction in the region as a reference, defining the change trend, setting the included angle threshold, extracting the group of direction deviation interval graphs, arranging the paths according to the gradient consistency, and extracting the main control path band of the reservoir deformation trend, the specific steps are as follows: S301: Obtain the index offset data in the region, screen according to the set offset direction threshold, extract the graph group that meets the offset direction requirements, obtain the deviation value between the direction included angle of adjacent units and the set direction threshold, calculate the included angle deviation degree, and form the offset direction distribution result of each unit; S302: Based on the offset direction distribution result of each unit, screen the graph group with an included angle deviation direction threshold between adjacent units, analyze the offset value of the graph group, arrange the distribution of the index difference according to the gradient consistency, and generate the trend path data; S303: Based on the trend path data, combined with the deformation characteristics of different units in the region, calculate the deformation amplitude and trend change of each path, compare with other path bands, and screen the path band with continuous deformation trend to obtain the main control path band of the reservoir deformation trend.

[0012] As a further solution of the present invention, the specific calculation formula of the included angle deviation degree is as follows: ; wherein, represents the included angle deviation degree between the unit and the unit in the graph group , represents the actual observed included angle value between the unit and the unit in the graph group , represents the set offset direction threshold angle, represents the average value of the observed included angles between all units in the graph group , represents the unit in the graph group and the unit The horizontal projection distance of the center point, which is calculated through the vector boundary in the geographic information system, Representative unit And the unit The boundary contact tightness interference coefficient between them, Representative graphic group The total number of all unit pairs with boundary contact relationships in it.

[0013] As a further solution of the present invention, the specific steps to call the main control path band of the reservoir deformation trend, select the pore-fracture contact data in the path crossing area, screen the contact point groups according to the included angle ratio, and obtain the pore-fracture co-response manifestation block set are as follows: S401: Call the main control path band of the reservoir deformation trend, extract the pore-fracture contact data in the area crossed by the path band, including the spatial distribution parameters of the pore axis direction and the main fracture distribution direction, combine the deformation trend of the path band, establish the spatial correlation matrix of the pore-fracture contact points, and obtain the projection distribution parameters of the pore-fracture contact points; S402: Based on the projection distribution parameters of the pore-fracture contact points, calculate the included angle between the pore axis direction and the main fracture distribution direction, set the included angle analysis interval, screen the points whose direction difference deviates from the set angle interval among all pore-fracture contact points, and count the position direction offset value of the contact points on the path band to obtain the offset angle screening contact point set; S403: Call the offset angle screening contact point set, analyze the distribution characteristics of the screening points in the area crossed by the path band, calculate the density distribution, aggregation trend and direction consistency between adjacent points of the screening points, and calibrate the blocks with direction co-response characteristics in the area to obtain the pore-fracture co-response manifestation block set.

[0014] As a further solution of the present invention, based on the pore-fracture co-response manifestation block set, combined with the combined distribution of the amplitude attenuation rate and the pressure inversion value, the specific steps to screen the spatial units with stress retention and fluid inclusion characteristics to obtain the reservoir space block delineation map set are as follows: S501: Obtain the amplitude attenuation rate and the pressure inversion value in the block, and conduct distribution analysis. Through the calculation of the amplitude attenuation rate and the pressure inversion value, obtain the amplitude attenuation and pressure inversion value intervals in the block, and respectively allocate the amplitude attenuation rate and the pressure inversion value of the block to intervals to obtain the amplitude attenuation and pressure inversion value distribution data corresponding to the block; S502: Based on the distribution of the amplitude attenuation rate and the pressure inversion value in the block, screen the spatial units with stress retention and fluid inclusion characteristics, and use the amplitude attenuation rate and the pressure inversion value in the block for analysis to obtain the stress retention characteristic block and the fluid inclusion characteristic block; S503: Based on the stress retention characteristic block and the fluid inclusion characteristic block, and in combination with the pore-fracture collaborative response information, aggregate the spatial units that meet the stress retention and fluid inclusion characteristics to obtain a set of reservoir space block delineation maps.

[0015] A shale gas reservoir space detection system, comprising: The reservoir layer section division module calculates the unilateral average gradient value of the longitudinal change of pores based on pore size and horizon point data, judges the stability of the layer section structure, selects continuous abnormal points as the interlayer transition zone, and divides physical property independent layer sections according to the position of the transition zone to obtain the reservoir layer division boundary result; The wave velocity correction analysis module calculates the difference between the original wave velocity value of the layer section and the pore mean value based on the reservoir layer section division result, adjusts the wave velocity response intensity in combination with local stress interference data, and obtains a set of reservoir internal structure response difference indexes by comparing the stress deformation amplitude offset degree before and after correction; The reservoir deformation trend analysis module, based on the set of reservoir internal structure response difference indexes, uses the index offset direction within the selected area as the judgment reference, delimits the continuous change trend, sets the direction angle threshold, and extracts the graphic group in which the direction angle between adjacent units deviates from the threshold to obtain the main control path band of the reservoir deformation trend; The pore-fracture collaborative response screening module, based on the main control path band of the reservoir deformation trend, performs single-point comparison of the angle between the pore axis direction and the main fracture distribution direction within the traversed area, and screens the contact points where the direction difference deviates from the set angle range to obtain the pore-fracture collaborative response manifestation block set; The reservoir space block identification module, based on the pore-fracture collaborative response manifestation block set, combines the amplitude attenuation rate and the combined distribution of pressure inversion value intervals within the block, screens the spatial units with stress retention and fluid inclusion characteristics, and obtains the reservoir space block delineation map.

[0016] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In the present invention, through the extraction of pore structure and horizon data, the fine division of the layer section boundary is realized, the structure differentiation recognition ability is improved, physical property mixing is avoided, the wave velocity difference and local stress interference are integrated, the stress response identification accuracy of microscopic units is enhanced, the internal structure difference of the reservoir is reflected, the pore axis and the fracture distribution direction are matched, the spatial coupling behavior is identified, and the amplitude attenuation and pressure inversion distribution are combined to accurately locate the reservoir units with stress retention and fluid inclusion characteristics, improving the accuracy and practicality of spatial identification. Description of the Drawings

[0017] To more clearly illustrate the technical solutions in the embodiments of the present invention, the following will briefly introduce the accompanying drawings required for the description of the embodiments. Obviously, the accompanying drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other accompanying drawings can be obtained based on these drawings.

[0018] Figure 1 It is a schematic diagram of the step flow of the present invention; Figure 2 It is a system module diagram of the present invention. Specific embodiments

[0019] The following will describe the technical solutions in the present invention in conjunction with the accompanying drawings.

[0020] In the embodiments of the present invention, words such as "exemplarily" and "for example" are used to represent examples, illustrations or explanations. Any embodiment or design solution described as an "example" in the present invention should not be construed as being more preferred or having more advantages than other embodiments or design solutions. Precisely, the use of the word "example" is intended to present concepts in a specific manner. In addition, in the embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one of the two can be selected.

[0021] In the embodiments of the present invention, "image" and "picture" can sometimes be used interchangeably. It should be noted that when their differences are not emphasized, the meanings they express are the same. "(of)", "corresponding", and "corresponding" can sometimes be used interchangeably. It should be noted that when their differences are not emphasized, the meanings they express are the same.

[0022] In the embodiments of the present invention, sometimes subscripts such as W1 may be written in a non-subscript form such as W1. When their differences are not emphasized, the meanings they express are the same.

[0023] To make the technical problems, technical solutions and advantages to be solved by the present invention clearer, the following will be described in detail in conjunction with the accompanying drawings and specific embodiments.

[0024] The embodiments of the present invention provide a method for detecting shale gas reservoir spaces, including the following steps: S1: Collect data on the pore sizes and layer position points at different depths of the reservoir, judge the stability of the layer structure according to the unilateral average gradient value of the longitudinal change of the pores, select continuous abnormal points as the interlayer transition zone, divide the physical property independent layer segments according to the position of the transition zone, and generate the result of dividing the reservoir stratification boundary; S2: Based on the reservoir layering boundary division results, calculate the difference between the original wave velocity value and the pore mean value in the layer segment. According to the difference and the local stress interference data, adjust the wave velocity response intensity, and compare the offset degree of the stress deformation amplitude before and after correction to form a set of reservoir internal structure response difference indicators; S3: Invoke the set of structure response difference indicators. Taking the index offset direction in the region as the judgment benchmark, delimit the continuous change trend, set the direction angle threshold, and extract the graphic group in adjacent units where the direction angle deviates from the direction threshold. Arrange the dense area of the index difference distribution in the graphic group according to the gradient consistency as the trend path, and extract the main control path band of the reservoir deformation trend; S4: Invoke the main control path band of the reservoir deformation trend, select the pore-fracture contact data in the area traversed by the path band, perform a single-point comparison on the angle between the pore axis direction and the main fracture distribution direction, and screen the contact point group whose direction difference deviates from the set angle interval to obtain the block set where the pore-fracture collaborative response appears; S5: Based on the block set where the pore-fracture collaborative response appears, combined with the combined distribution of the wave amplitude attenuation rate and the pressure inversion value interval in the block, screen the spatial units with stress retention and fluid inclusion characteristics to obtain the map set for delimiting the reservoir space block.

[0025] The reservoir layering boundary division results include the distribution of the interlayer transition zone, the type of physical property independent layer segments, the stability level of the layer segment structure, and the boundary division accuracy level. The reservoir internal structure response difference indicators include the wave velocity and pore difference distribution, the local stress disturbance intensity, the stress deformation amplitude offset, and the wave velocity response correction coefficient. The main control path band of the reservoir deformation trend includes the consistency of the index offset direction, the continuity of the trend path, the area where the direction angle deviates, and the gradient consistency section. The block where the pore-fracture collaborative response appears includes the pore-fracture contact point distribution, the direction angle difference interval, the collaborative response appearance intensity, and the pore-fracture structure matching relationship. The map set for delimiting the reservoir space block includes the wave amplitude attenuation rate interval, the combined distribution of the pressure inversion value, the identification result of the stress retention area, and the fluid inclusion characteristic unit.

[0026] The specific steps of S1 are as follows: S101: Collect the pore size and layer position data of different depths in the reservoir, calculate the pore size gradient value corresponding to the layer position, calculate the unilateral average gradient value based on the pore size change of the layer, extract the distribution sequence of the unilateral average gradient value of the layer segment, and calculate the gradient change rate of the layer segment for the gradient values of adjacent layers in the distribution sequence; When collecting reservoir differential depth pore size and horizon point position data, it is necessary to obtain the pore size values for specific well sections relying on methods such as core analysis, mercury injection experiments, and nuclear magnetic resonance logging. At the same time, combine natural gamma, resistivity curves, and imaging logging results to calibrate the depth information of each horizon point position. Each horizon point position needs to include its vertical depth and the pore size at that position. After arranging the data in depth order, compare the pore size and depth values of adjacent points pairwise. Reflect the pore size change rate, that is, the gradient value, through the relationship between the increase or decrease of the pore size and the depth change amount between two points. When executing, calculate pairwise according to the points and form a continuous sequence of the results. Then classify the point segments with the same direction in the gradient sequence, and judge whether the continuous positive or negative change segments meet the minimum number of points requirements. If they meet, extract all the gradient data in this segment for averaging to form the average gradient value of this unilateral segment. Subsequently, perform the same operation on all continuous unilateral gradient segments in the reservoir, record their starting depth position, ending depth position, and average gradient results to form a complete unilateral average gradient distribution sequence. Further, by comparing the average values of adjacent unilateral average gradient segments and combining the central depth position between the two segments, calculate the average gradient change amplitude between any adjacent segments in turn to judge the amplitude level of the gradient change. The larger the value, the more significant the change. The result of dividing the gradient difference between each pair of adjacent segments in this sequence by the center depth spacing between the two segments constitutes the layer segment gradient change rate sequence, which is used as the basic analysis data for subsequent reservoir stratification structure identification and division.

[0027] S102: Call the layer segment gradient change rate, screen the positions of abnormal mutations, judge the continuous abnormal points of the positions, calculate the distribution interval of the abnormal points, delimit the range of the interlayer transition zone based on the distribution interval of the abnormal points, and generate the position distribution of the interlayer transition zone; When calling the data of the gradient change rate of the layer section, it is necessary to screen each point of the obtained gradient change rate distribution sequence, and set a threshold range for mutation judgment. This threshold needs to be obtained through the statistical summary of the logging analysis results of multiple well sections. Taking the distribution law of the gradient change rate as the benchmark, select the numerical range at the extreme distribution position in the change rate sequence as the judgment standard. If the gradient change rate of a certain point exceeds this range, mark it as a mutation point. After screening out all mutation points, further judge its continuity in the depth direction, and detect whether the mutation points are adjacent in the vertical direction and meet the minimum continuous quantity requirement. For example, three or more adjacent mutation points are grouped into an abnormal point group, and record the depths corresponding to the starting point and the ending point in this point group as the upper and lower boundaries of this abnormal point group, forming a complete abnormal distribution section. Process all eligible abnormal point groups in turn, and finally count the depth ranges of all abnormal sections, integrate all effective continuous abnormal sections, identify the positions where there are transition zones between reservoir layers, and check again whether each transition zone section meets the minimum thickness requirement and whether the adjacent sections meet the minimum spacing setting requirement. After screening out the results that do not meet the conditions, retain the remaining effective transition zone positions and their distribution information, and output them as the distribution result of the interlayer transition zone position, providing basic support data for the subsequent division of independent reservoir layers.

[0028] S103: Call the distribution of the interlayer transition zone position, and divide the reservoir into physical property independent layers according to the position of the transition zone in the layer section, and adjust the boundaries of the independent layers to generate the result of the division of the reservoir layer boundaries; After calling the result of the distribution of the interlayer transition zone position, perform spatial matching within the original reservoir section according to the depth interval where each transition zone is located, and divide the original reservoir section according to the position of the transition zone. The execution method is: if there are two or more transition zone intervals in an original reservoir section, then use the median depth position of each transition zone as the cutting point to split the reservoir section into multiple physical property independent sub-sections. For example, the starting and ending depths remain unchanged, and the middle is segmented in turn according to the position of the transition zone. After the preliminary division, perform an adjustment operation on the boundary of each newly generated independent sub-section. When adjusting, it is necessary to call the pore size values of its boundary points and several points before and after, and judge whether there is a significant mutation by calculating the change degree of the pore size of the points near the boundary. If the difference in the average pore size of several points before and after the boundary exceeds the preset judgment value, adjust the original boundary to the position where the mutation occurs to complete the boundary correction. Each sub-section after correction should meet the requirement that the thickness is greater than the minimum division thickness. Finally, organize and summarize the starting depth, ending depth and boundary adjustment results of all eligible independent layers to form the final result dataset of the division of the reservoir layer boundaries, which is used as an important basis for constructing the geological model and reservoir classification.

[0029] The specific steps of S2 are as follows: S201: Based on the reservoir layering boundary division result, obtain the original value of the wave velocity in the layer section, extract the average value of the porosity, calculate the difference between the original value and the average value of the porosity, compare the original wave velocity value with the pore mean value in the layer section, obtain the difference data for each layer section, and record the difference as the layer section difference data; Based on the reservoir layering boundary division result, it is necessary to extract the original longitudinal wave velocity data corresponding to each physical property independent layer section that has been clarified. This data is sourced from acoustic logging data. During the acquisition process, the longitudinal wave slowness of the acoustic logging sampling points is extracted according to the corresponding depth range of the layer section, converted to the longitudinal wave velocity per unit depth, and then the porosity value is obtained by calling the conventional logging data or core analysis data within the same depth range. The porosity data within this section is averaged, that is, the porosity of all measurement points within this section is statistically analyzed and then averaged. After completing the extraction of this average value, the difference between the original wave velocity value and the average porosity value of the layer section is processed, that is, the wave velocity value of each sampling point is subtracted by the average porosity value of this section to reflect the degree of difference in the wave velocity response of the internal medium structure of the layer section. For example, if the original wave velocity value of a certain layer section is about 3500 m / s and the average porosity value is 0.17 (i.e., 17%), then the wave velocity data of each sampling point needs to be subtracted by 0.17 one by one to form the difference data. The same operation is performed successively for all layer sections to obtain the difference data for each layer section. The difference data needs to be archived according to the layer section number, start and end depths, and corresponding values, and recorded as the layer section difference data.

[0030] S202: Based on the layer section difference data, monitor and collect local stress interference data, compare the difference of each layer section with the corresponding stress interference data, adjust the wave velocity response intensity, obtain the wave velocity correction value of the layer section under stress interference, compare the differential correction values, and generate the wave velocity correction intensity data; First, micro - stress or acoustic emission monitoring means need to be carried out in the target well section to collect local stress interference data, which reflects the degree of external disturbance or tectonic stress influence on the current reservoir. During the extraction process, the measuring point positions are compared according to the depth range of the layer section, and the stress disturbance values are aggregated by layer section. At the same time, a one - to - one correspondence between the layer - section difference data and the corresponding stress disturbance values is established, and the mutual change characteristics between the difference and the stress disturbance value are judged section by section. During the comparison process, a classification processing method is adopted. The difference is divided into three grades. For example, the difference in the range of - 200 to - 100 is the low - difference section, - 100 to 100 is the medium - difference section, and above 100 is the high - difference section. At the same time, the stress disturbance values are divided into three levels according to the amplitude: weak disturbance (less than 30 MPa), medium disturbance (30 to 60 MPa), and strong disturbance (greater than 60 MPa). In each layer section, judge the combination type of its difference and stress level. For example, if the difference in a certain layer section is 150 and the disturbance value is 45 MPa, belonging to the combination of high - difference and medium - disturbance, then an adjustment action of the wave - velocity response intensity needs to be performed. The adjustment operation is to set a fixed correction ratio for the original wave - velocity value according to the stress disturbance level. The ratio setting is statistically obtained through historical seismic response data. For example, if the correction ratio for the medium - disturbance level is set to 5%, then the original wave - velocity data is reduced by 5% as the corrected value. Repeat this processing process for all layer sections in turn. After adjusting the original wave - velocity values of different layer sections under the corresponding stress disturbances, new wave - velocity correction values are obtained. By comparing the wave - velocity differences before and after correction, the response degree of each layer section to the disturbance influence can be reflected, and finally, the wave - velocity correction intensity data of each layer section is recorded.

[0031] S203: Based on the wave - velocity correction intensity data, compare the stress - deformation amplitudes before and after correction, calculate the stress - deformation offset of the layer section, obtain the difference value of the structural response within the layer section, and form a group of structural response difference indexes within the reservoir according to the difference value; The specific calculation formula for the stress - deformation offset of the layer section is: ; Where, represents the stress - deformation offset of the th layer section, represents the original wave - velocity value before wave - velocity correction of the th sampling point in the th layer section, represents the wave - velocity value after wave - velocity correction of the th sampling point in the th layer section, represents the local stress disturbance value of the th sampling point in the th layer section, represents the stress - disturbance adjustment coefficient of the th layer section, and its value is determined by linear regression based on the average disturbance level of the previous measuring points. represents the average value of all original wave velocity sampling points in the interval, represents the corrected balance factor of the interval; Suppose there are 5 sampling points in a certain interval , and the wave velocity data before and after correction and the stress disturbance values at each point are measured as follows: , , , , , , , , , , ; According to the adjustment coefficient setting: , this coefficient reflects the average adjustment level of the interval and can be determined through the correlation analysis of the average disturbance value and stress response of the previous measurement points , this value is used to smooth the influence of variance, is applicable to small sample data sets, and improves the stability of the results; Calculation process: First, calculate the part related to the wave velocity change and stress disturbance: ; ; Secondly, calculate the variance part of the difference between the original wave velocity and the average wave velocity. First, calculate the average value of the original wave velocity: ; Then calculate the variance part: ; ; Summary: ; This result shows that in the given interval, the calculation related to the wave velocity and stress disturbance shows an obvious stress deformation offset, the relationship between the wave velocity change and stress, and its influence on the overall structural response, indicating the sensitivity of the interval to external changes. This value is further used to form a group of structural response difference indicators in the reservoir, providing key data support for reservoir management.

[0032] The specific steps of S3 are as follows: S301: Obtain the index offset data in the area, screen according to the set offset direction threshold, extract the graphic group that meets the offset direction requirements, obtain the deviation value of the included angle between adjacent units and the set direction threshold, calculate the included angle deviation degree, and form the offset direction distribution result of each unit; The specific calculation formula of the included angle deviation degree is: ; Among them, represents the angular deviation of the unit in the graphic group from the unit and the unit ; represents the actual observed angular value of the unit in the graphic group from the unit and the unit ; represents the set threshold angle of the offset direction represents the average value of the observed angles between all units in the graphic group ; represents the horizontal projection distance of the unit in the graphic group from the center point of the unit and the unit , calculated through the vector boundary in the geographic information system represents the boundary contact tightness interference coefficient between the unit and the unit ; represents the total number of pairs of units with boundary contact relationships in the graphic group ; Considering a specific example, assume that the graphic group contains four units , and the specific observed angular value has been obtained from GIS data: , , The set threshold angle of the offset direction is calculated The horizontal projection distance of the center point m, and the boundary contact tightness interference coefficient is calculated through GIS data. It is known from the interface data that it is a function of the ratio of the unit boundary length to the area and the total number of adjacent unit pairs in the graphic group ; Calculation process: ; ; ; This result shows that the unit and The angular deviation is 141.5°, far exceeding the set threshold, indicating that the structural changes between the units in this area are relatively significant. The larger the value of the angular deviation, the more drastic the actual observed direction change is relative to the preset direction threshold. This is crucial for analyzing the dynamic changes of geological structures and their impact on reservoir performance, and further guiding reservoir management and development decisions. The calculation of this value takes into account the direct deviation of the angle and the influence after adjusting for distance and contact tightness, providing a method that comprehensively considers spatial relationships and geometric configurations to evaluate the interaction between units.

[0033] S302: Based on the offset direction distribution results of each unit, screen the graphic groups of the angular deviation direction thresholds between adjacent units, analyze the offset values of the graphic groups, arrange the distribution of the index differences according to gradient consistency, and generate the strike path data; First, spatial coding processing needs to be carried out on each geological unit in the divided area of the target reservoir, and the index offset values formed by each unit under different tectonic stages or different stress states are extracted. This offset value can be from the data results generated by previous modules such as wave velocity correction intensity and stress deformation offset. After data collection for all units in sequence, the offset direction threshold is set as the screening basis. This threshold needs to be set according to the main direction of tectonic stress in this area and the distribution of offset angles to set the upper and lower boundary limits. For example, the offset direction threshold is set to 45° ± 10°, that is, the unit data with the offset direction between 35° and 55° is screened. When performing the screening operation, it is judged whether the offset direction of each unit falls within this interval. If so, the unit is extracted to form a graphic group. Subsequently, the spatial relationship between adjacent units in each graphic group is analyzed, and their direction angles are extracted pair by pair. The direction angle value is obtained by calculating the difference between the main offset direction angles of the two units, and then the difference operation is performed between this angle value and the set offset direction threshold to obtain its deviation degree in the direction. According to the size of the deviation value, it can be classified into different levels. For example, the deviation value less than 5° is the first level, 5° to 10° is the second level, and greater than 10° is the third level. In this way, the offset direction distribution results of each unit in the entire area are formed, which includes content such as graphic number, position, main direction angle, angular deviation degree, etc.

[0034] S303: Based on the strike path data, combined with the deformation characteristics of different units in the area, calculate the deformation amplitude and trend change of each path, and compare with other path bands to screen the path bands with continuous deformation trends to obtain the main control path band of reservoir deformation trends; First, it is necessary to traverse and extract the node sequences inside each path, and call three types of data: the index offset direction, offset amount, and included angle deviation degree of the nodes on each path. By sorting and segmenting the direction change trends in the node sequence, it is judged whether there are obvious turning points or angle mutation points inside the path. If the angle difference between the offset direction of a certain node and the previous node exceeds the set included angle threshold, then this point is defined as the direction mutation point inside the path, and the path segmentation is recorded. At the same time, the cumulative difference of the offset amount trend on the path is judged. If the total offset amount is always in the same change direction, it means that this path segment is a stable offset path segment. After performing this type of analysis on all paths, select the path segments with continuous path length, stable direction offset angle, and large cumulative offset amount, and define them as the main control path band of the deformation trend. The screening criteria are set as follows: the path length is not less than 5 units, the total change of the offset direction included angle does not exceed 20°, and the total change value of the offset amount is greater than 100 units. The path segments that meet this screening condition are finally defined as the main control path band of the reservoir deformation trend. Each node in this path band records elements such as its spatial position, direction, offset amount, and deviation degree level, and is output in the form of a graphic sequence for regional reservoir structure analysis.

[0035] The specific steps of S4 are as follows: S401: Call the main control path band of the reservoir deformation trend, extract the pore-fracture contact data within the area traversed by the path band, including the spatial distribution parameters of the pore axis direction and the main fracture distribution direction, and combine the deformation trend of the path band to establish the spatial correlation matrix of the pore-fracture contact points and obtain the projection distribution parameters of the pore-fracture contact points; First, according to the fault, fracture and tectonic surface deformation characteristics in the 3D geological model of the reservoir, the target area of the reservoir is divided into a clear path zone area. For example, a spatial range from coordinates (X = 200m, Y = 800m) to coordinates (X = 500m, Y = 1200m) is determined as the path zone crossing area. Subsequently, in the reservoir physical property database, all pore-fracture contact point data within this area are retrieved. For example, the axis direction angle of pore point P1 is 35°, and the corresponding fracture distribution direction angle is 50°; the axis direction angle of pore point P2 is 40°, and the corresponding fracture distribution direction angle is 55°. Such data are read one by one. After collecting all pore-fracture point data within the path zone, the pore axis direction and the main fracture distribution direction are projected onto the main control direction of the path zone (for example, the overall trend direction of the path zone is taken as 45° northeast). Taking the starting point of the path zone as the coordinate zero point, the projected positions of the pore-fracture contact points are measured and recorded along the path zone direction. For example, the projection of pore-fracture contact point P1 is 120m from the starting point of the path zone, and the projection of pore-fracture contact point P2 is 165m from the starting point of the path zone. Then, the pairwise relative position distances and the corresponding axis direction angle differences between the pore-fracture contact points within the path zone are calculated. For example, the axis direction angle difference between point P1 and P2 is 5°, and the spatial projection distance difference is 45m. Similar calculations are carried out one by one to form a spatial correlation matrix, where each cell in the matrix records the position distance and the axis direction angle difference between the corresponding two points. Finally, through matrix analysis, the projection position distribution parameters, axis angle difference parameters and spatial distance parameters of the pore-fracture contact points are extracted to form a set of projection distribution parameters of the pore-fracture contact points.

[0036] S402: Based on the projection distribution parameters of the pore-fracture contact points, calculate the angle between the pore axis direction and the main fracture distribution direction, set an angle analysis interval, screen out the points whose direction differences deviate from the set angle interval among all pore-fracture contact points, and count the position direction offset values of the contact points on the path zone to obtain a set of offset angle screening contact points; First, extract the projection angle between the pore axis direction and the main crack distribution direction of the pore-crack contact points one by one from the obtained parameter set. For example, if the axis direction of the pore-crack contact point P3 is 25° and the corresponding main crack distribution direction is 40°, the calculated angle difference is 15°. Subsequently, set the angle analysis interval. Based on the reservoir deformation trend and actual production experience, it is stipulated that the angle difference within the range of 10° to 20° is the effective angle analysis interval. Then, calculate the projection angles of all pore-crack contact points within the path zone one by one. Suppose the calculated angle of the pore-crack contact point P4 is 25°. Compared with the set angle analysis interval, this angle exceeds the interval range, so the pore-crack contact point P4 is marked as a direction difference deviation point. Next, using the main control direction angle of the path zone (such as 45° in the northeast direction) as a reference, calculate the relative offset angles between the axis directions of each direction difference deviation point and the main control direction of the path zone. For example, if the axis direction of point P4 is 70°, the direction offset angle of point P4 is 70° - 45° = 25°. Calculate the offset angles of all direction difference deviation points within the path zone area in sequence, and statistically obtain the offset angle set of all direction difference deviation pore-crack contact points, that is, obtain the offset angle screening contact point set.

[0037] S403: Call the offset angle screening contact point set, analyze the distribution characteristics of the screening points within the path zone crossing area, calculate the density distribution, aggregation trend, and direction consistency between adjacent points of the screening points, calibrate the blocks with direction coordination characteristics within the area, and obtain the pore-crack coordination response manifestation block set; First, extract the position coordinates one by one from the contact point set for distribution feature analysis. For example, within the range of the path band crossing area (such as from (X = 200m, Y = 800m) to (X = 500m, Y = 1200m)), statistically analyze the spatial position coordinates of the direction difference deviating from the contact point. For example, point P5 is located at (X = 300m, Y = 900m), and point P6 is located at (X = 320m, Y = 920m). Complete the marking of all screened point coordinates one by one, and then conduct a density distribution statistics of the screened points. The specific method is to divide the path band area into regular grid cells (such as using a 50m×50m square), and count the number of direction deviations from the contact point within each grid cell. If there are 6 screened points in a grid cell, the density is 0.0024 points / m². Then, determine the aggregation trend, that is, calculate the distance between each screened point and its nearest neighbor point. For example, calculate the distance between point P5 and its nearest neighbor point P6, which is approximately 28m. Set the distance threshold between adjacent points to 30m. When the distance between a point and its nearest point is less than the threshold of 30m, it is determined to be in an aggregated state. Next, calculate the consistency of the axis directions of each point within each aggregated state point set, specifically by calculating the average value and standard deviation of the direction angles of the axis directions of the points in the aggregated state point set. For example, in the aggregated state point set, there are point P5 (axis direction 60°), point P6 (axis direction 62°), and point P7 (axis direction 61°). Calculate the average value of their axis directions to be 61°, and the standard deviation to be 0.82°. Set the standard deviation determination threshold for direction consistency to 5°. If the standard deviation is less than 5°, it is determined that the direction consistency of the aggregated point set is strong. Complete the analysis of the direction consistency of the screened points within the entire path band area. Finally, calibrate the position of the area with strong direction consistency through coordinate recording (such as the area range is from 300m to 350m in the X direction and from 900m to 950m in the Y direction) to obtain the block set of the co - response manifestation of pore - fracture.

[0038] The specific steps of S5 are as follows: S501: Obtain the amplitude attenuation rate and pressure inversion value within the block and conduct distribution analysis. Through the calculation of the amplitude attenuation rate and pressure inversion value, obtain the amplitude attenuation and pressure inversion value intervals within the block. Allocate the amplitude attenuation rate and pressure inversion value of the block into intervals respectively to obtain the amplitude attenuation and pressure inversion value distribution data corresponding to the block; First, retrieve the initial and end values of the wave amplitude at each measuring point in the seismic acquisition database for this block. For example, at measuring point A, the initial wave amplitude value is 200 and the end wave amplitude value is 120. The wave amplitude attenuation rate of 0.2 / m is obtained by dividing the difference between the initial and end wave amplitude values (200 - 120 = 80) by the propagation distance (e.g., 400m). After repeating similar calculations for other measuring points, all the wave amplitude attenuation rate data within the block are obtained. At the same time, retrieve the pressure inversion numerical values at each measuring point in the block from the pressure inversion dataset. For example, the pressure inversion value at measuring point A is 35 MPa and at measuring point B is 42 MPa. Call and record the pressure inversion value data of all measuring points in the block one by one. Then, based on the calculated data of the wave amplitude attenuation rate of all measuring points, determine the maximum and minimum values of the attenuation rate (e.g., 0.1 / m to 0.5 / m). Divide this interval with a step size of 0.1 / m, which are four sub - intervals: 0.1 - 0.2 / m, 0.2 - 0.3 / m, 0.3 - 0.4 / m, and 0.4 - 0.5 / m. And according to the wave amplitude attenuation rate numerical values of the measuring points, classify each measuring point into the corresponding interval. For example, the attenuation rate of measuring point A, 0.2 / m, belongs to the 0.1 - 0.2 / m interval. Then, use all the pressure inversion value data within the block to determine the maximum and minimum values of the pressure inversion value. For example, in the interval of 25 MPa to 50 MPa, divide it into sub - intervals according to 5 MPa as one interval, such as 25 - 30 MPa, 30 - 35 MPa, 35 - 40 MPa, 40 - 45 MPa, 45 - 50 MPa. Then, allocate each measuring point to the corresponding interval according to the pressure inversion value. For example, the pressure inversion value of measuring point A is 35 MPa, which belongs to the 35 - 40 MPa interval. Through the above steps, complete the interval allocation of the wave amplitude attenuation rate and pressure inversion value of all measuring points in the block one by one, count the number of measuring points corresponding to each sub - interval, and record the distribution position coordinate data, so as to obtain the wave amplitude attenuation and pressure inversion value distribution data corresponding to the block.

[0039] S502: Based on the distribution of the wave amplitude attenuation rate and pressure inversion value within the block, screen the spatial units with stress retention and fluid inclusion characteristics, and use the wave amplitude attenuation rate and pressure inversion value within the block for analysis to obtain the stress retention characteristic block and the fluid inclusion characteristic block; First, the measured point data of the allocated intervals are called one by one to determine the spatial units with stress retention characteristics. The specific operation is to simultaneously call and numerically compare the wave amplitude attenuation rate and the pressure inversion value. The combined condition that the wave amplitude attenuation rate is greater than 0.3 / m and the pressure inversion value is higher than 40 MPa is set as the unit with stress retention characteristics. For example, the wave amplitude attenuation rate of measuring point C is 0.35 / m and the pressure inversion value is 42 MPa. After judgment, it meets the above conditions and is marked as the stress retention characteristic point. Then, another set of condition judgments is made on the wave amplitude attenuation rate and the pressure inversion value simultaneously. The combined condition that the wave amplitude attenuation rate is lower than 0.2 / m and the pressure inversion value is lower than 35 MPa is set as the fluid inclusion characteristic unit. For example, the wave amplitude attenuation rate of measuring point D is 0.15 / m and the pressure inversion value is 32 MPa. After this numerical judgment, if it meets the conditions, it is marked as the fluid inclusion characteristic point. Next, the coordinate data of all the measuring points with stress retention characteristics are centralized, and their spatial position coordinates in the block are recorded and the distribution quantity is counted to form the stress retention characteristic block. The same operation is also performed on the coordinate data of all the measuring points with fluid inclusion characteristics, centralizing the coordinate positions and counting the distribution to obtain the fluid inclusion characteristic block. Through the above operations, the stress retention characteristic block and the fluid inclusion characteristic block are obtained respectively.

[0040] S503: Based on the stress retention characteristic block and the fluid inclusion characteristic block, combined with the pore-fracture collaborative response information, aggregate the spatial units that meet the stress retention and fluid inclusion characteristics to obtain the mapping group for delineating the reservoir space block; First, call the pore-fracture collaborative response information, specifically, call the coordinate range of the block with the pore-fracture collaborative response manifestation characteristics obtained previously. For example, the coordinate range of the pore-fracture collaborative response manifestation block is 500m to 600m in the X coordinate and 1000m to 1100m in the Y coordinate. Compare the spatial range of this block with the spatial coordinate range of the stress retention characteristic block one by one to determine the spatial position relationship between the two blocks. If there is an overlap or spatial adjacency (the position offset does not exceed 50m) between the two blocks, record the overlapping area or adjacent area as the stress retention characteristic space unit with pore-fracture collaborative response. For example, the spatial coordinate range of the stress retention characteristic block is 550m to 650m in the X coordinate and 1050m to 1150m in the Y coordinate. After comparing and judging the numerical values of the coordinate position ranges, it is found that the two blocks overlap in the area of 550m to 600m in the X coordinate and 1050m to 1100m in the Y coordinate. Subsequently, perform the same call and coordinate range comparison operation on the spatial coordinates of the fluid inclusion characteristic block. For example, the coordinates of the fluid inclusion characteristic block are 480m to 560m in the X coordinate and 980m to 1060m in the Y coordinate. After comparing the coordinate ranges, there is an overlapping range of 500m to 560m in the X coordinate and 1000m to 1060m in the Y coordinate with the pore-fracture collaborative response block. Record this area as the fluid inclusion characteristic space unit with pore-fracture collaborative response. Finally, integrate the specific coordinate range data of the above two types of space units, aggregate the overlapping or adjacent areas of the spatial coordinates, and draw a schematic diagram of the spatial coordinate distribution to show the delineation result of the reservoir space block in a graphical manner, and finally obtain the reservoir space block delineation map set.

[0041] Please refer to Figure 2 , a shale gas reservoir space detection system, comprising: The reservoir layer section division module calculates the unilateral average gradient value of the longitudinal change of pores based on the pore size and layer position point data, judges the stability of the layer section structure, selects continuous abnormal points as the interlayer transition zone, and divides the physical property independent layer sections according to the position of the transition zone to obtain the division result of the reservoir layer boundary; The wave velocity correction analysis module calculates the difference between the original wave velocity value of the layer section and the pore mean value based on the reservoir layer section division result, adjusts the wave velocity response intensity in combination with the local stress interference data, and obtains the reservoir internal structure response difference index group by comparing the offset degree of the stress deformation amplitude before and after correction; The reservoir deformation trend analysis module selects the index offset direction in the selected area as the judgment benchmark based on the reservoir internal structure response difference index group, delimits the continuous change trend, sets the direction angle threshold, and extracts the graphic group in which the direction angle of adjacent units deviates from the threshold to obtain the main control path band of the reservoir deformation trend; The pore-fracture co-response screening module conducts single-point comparison of the angle between the pore axis direction and the main fracture distribution direction within the traversing area based on the main control path zone of reservoir deformation trend, and screens out the contact points where the direction difference deviates from the set angle range, obtaining the pore-fracture co-response manifestation block set; The reservoir space block identification module, based on the pore-fracture co-response manifestation block set, combines the amplitude attenuation rate and the combined distribution of pressure inversion value intervals within the block, and screens out the spatial units with stress retention and fluid inclusion characteristics, obtaining the reservoir space block delineation map.

[0042] As mentioned above, it is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should all be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claimed rights.

Claims

1. A method for detecting shale gas reservoir spaces, characterized in that, The system includes: S1: Collect data on pore sizes and formation position points at different depths in the reservoir. Determine the stability of the formation intervals based on the unilateral average gradient value of the longitudinal pore variation. Select the abnormal points as the transition zones, divide the physical property independent formation intervals according to their positions, and generate the result of the reservoir layer boundary division. S2: Based on the result of the reservoir layer boundary division, calculate the difference between the original wave velocity value and the pore mean value in the formation interval. Adjust the wave velocity response intensity according to the difference and the local stress interference data to form a set of reservoir internal structure response difference indicators. S3: Invoke the set of structure response difference indicators. Taking the index offset direction within the region as the reference, delimit the change trend, set the angle threshold, extract the graph set of the direction deviation intervals, arrange the paths according to the gradient consistency, and extract the main control path band of the reservoir deformation trend. S4: Invoke the main control path band of the reservoir deformation trend. Select the pore-fracture contact data in the area traversed by the path, screen the set of contact points according to the angle comparison, and obtain the block set of the pore-fracture collaborative response appearance. S5: Based on the block set of the pore-fracture collaborative response appearance, combine the combined distribution of the wave amplitude attenuation rate and the pressure inversion value, screen the spatial units with stress retention and fluid inclusion characteristics, and obtain the map set of the reservoir space block delimitation.

2. The shale gas reservoir space detection method according to claim 1, characterized in that: The result of the reservoir layer boundary division includes the distribution of the interlayer transition zone, the type of physical property independent formation intervals, the stability level of the formation interval structure, and the accuracy level of the boundary division. The set of reservoir internal structure response difference indicators includes the distribution of the difference between wave velocity and pores, the intensity of local stress disturbance, the offset amount of stress deformation amplitude, and the wave velocity response correction coefficient. The main control path band of the reservoir deformation trend includes the consistency of the index offset direction, the continuity of the trend path, the direction angle deviation area, and the gradient consistency section. The block of the pore-fracture collaborative response appearance includes the distribution of pore-fracture contact points, the difference interval of the direction angle, the appearance intensity of the collaborative response, and the matching relationship between pore-fracture structures. The map set of the reservoir space block delimitation includes the wave amplitude attenuation rate interval, the combined distribution of the pressure inversion value, the identification result of the stress retention area, and the fluid inclusion characteristic unit.

3. The shale gas reservoir space detection method according to claim 1, characterized in that: The specific steps for collecting data on pore sizes and formation position points at different depths in the reservoir, determining the stability of the formation intervals based on the unilateral average gradient value of the longitudinal pore variation, selecting the abnormal points as the transition zones, dividing the physical property independent formation intervals according to their positions, and generating the result of the reservoir layer boundary division are as follows: S101: Collect data on pore sizes and formation position points at different depths in the reservoir. Calculate the pore size gradient value corresponding to the formation position point. Based on the pore size change situation of the formation, calculate the unilateral average gradient value, extract the distribution sequence of the formation interval unilateral average gradient value, and calculate the gradient change rate of the formation interval for the gradient values of adjacent formations in the distribution sequence. S102: Invoke the gradient change rate of the formation interval, screen the abnormally mutated positions, judge the continuous abnormal points of the positions, calculate the distribution interval of the abnormal points, delimit the range of the interlayer transition zone according to the distribution interval of the abnormal points, and generate the position distribution of the interlayer transition zone. S103: Invoke the position distribution of the interlayer transition zone. Divide the reservoir into physical property independent formation intervals according to the position of the transition zone in the formation interval, and adjust the boundaries of the independent formation intervals to generate the result of the reservoir layer boundary division.

4. The shale gas reservoir space detection method according to claim 1, wherein: Based on the reservoir layering boundary division result, calculate the difference between the original wave velocity value and the pore mean value in the layer segment. According to the difference and the local stress interference data, adjust the wave velocity response intensity. The specific steps to form the structural response difference index group within the reservoir are as follows: S201: Based on the reservoir layering boundary division result, obtain the original wave velocity value in the layer segment, extract the mean porosity, calculate the difference between the original value and the mean porosity, compare the original wave velocity value and the pore mean value in the layer segment, obtain the difference data for each layer segment, and record the difference as the layer segment difference data; S202: Based on the layer segment difference data, monitor and collect local stress interference data, compare the difference of each layer segment with the corresponding stress interference data, adjust the wave velocity response intensity, obtain the wave velocity correction value of the layer segment under stress interference, compare the differential correction values, and generate wave velocity correction intensity data; S203: Based on the wave velocity correction intensity data, compare the stress deformation amplitudes before and after correction, calculate the stress deformation offset of the layer segment, obtain the difference value of the structural response within the layer segment, and form the structural response difference index group within the reservoir according to the difference value.

5. The shale gas reservoir space detection method according to claim 1, wherein: The specific calculation formula for the stress deformation offset of the layer segment is: ; Among them, represents the stress and deformation offset of the th layer segment, represents the original wave velocity value before wave velocity correction at the th sampling point in the th layer segment, represents the wave velocity value after wave velocity correction at the th sampling point in the th layer segment, represents the local stress disturbance value at the th sampling point in the th layer segment, represents the stress disturbance adjustment coefficient of the th layer segment, and its value is determined by linear regression based on the average disturbance level of the previous measurement points, represents the average value of all original wave velocity sampling points in the th layer segment, represents the correction balance factor of the th layer segment.

6. The shale gas reservoir space detection method according to claim 1, wherein: The specific steps to call the structural response difference index group, delimit the change trend based on the index offset direction within the region, set the included angle threshold, extract the graphic group of the direction deviation interval, and arrange the paths according to the gradient consistency to extract the main control path band of the reservoir deformation trend are as follows: S301: Obtain the index offset data within the region, screen according to the set offset direction threshold, extract the graphic group that meets the offset direction requirements, obtain the deviation value between the included angle of adjacent units and the set direction threshold, calculate the included angle deviation degree, and form the offset direction distribution result of each unit; S302: Based on the offset direction distribution result of each unit, screen the graphic group with an included angle deviation from the direction threshold between adjacent units, analyze the offset value of the graphic group, arrange the distribution of the index difference according to the gradient consistency, and generate the trend path data; S303: Based on the trend path data, combine the deformation characteristics of different units within the region, calculate the deformation amplitude and trend change of each path, compare with other path bands, and screen the path band with continuous deformation trend to obtain the main control path band of the reservoir deformation trend.

7. The shale gas reservoir space detection method according to claim 1, characterized in that: The specific calculation formula for the included angle deviation degree is: ； Among them, represents the graphic group in the unit and the deviation angle of the unit The included angle deviation, represents the graphic group in the unit and the unit The actual observed included angle value, represents the set offset direction threshold angle, represents the graphic group The average value of the observed included angles between all units in, represents the graphic group in the unit and the unit The horizontal projection distance of the center point, calculated by the vector boundary in the geographic information system, represents the unit and the unit The boundary contact tightness interference coefficient between, represents the graphic group The total number of pairs of units with boundary contact relationships in all.

8. The shale gas reservoir space detection method according to claim 1, characterized in that: The specific steps to call the main control path band of the reservoir deformation trend, select the pore-fracture contact data in the area where the path passes through, screen the contact point group according to the included angle comparison, and obtain the block set where the pore-fracture co-response appears are as follows: S401: Call the main control path band of the reservoir deformation trend, extract the pore-fracture contact data in the area where the path band passes through, including the spatial distribution parameters of the pore axis direction and the main fracture distribution direction, combine the deformation trend of the path band, establish the spatial correlation matrix of the pore-fracture contact points, and obtain the projection distribution parameters of the pore-fracture contact points; S402: Calculate the angle between the pore axis direction and the main fracture distribution direction based on the projection distribution parameters of the pore-fracture contact points, set the angle analysis interval, screen the points whose direction differences deviate from the set angle interval among all the pore-fracture contact points, and statistically calculate the position direction offset value of the contact points on the path band to obtain the offset angle screening contact point set; S403: Call the offset angle screening contact point set, analyze the distribution characteristics of the screened points within the area where the path band passes through, calculate the density distribution, aggregation trend and direction consistency between adjacent points of the screened points, and demarcate the blocks with direction coordination characteristics within the area to obtain the pore-fracture coordination response apparent block set.

9. The shale gas reservoir space detection method according to claim 1, wherein: Based on the pore-fracture coordination response apparent block set, combined with the combined distribution of the amplitude attenuation rate and the pressure inversion value, the specific steps for screening the spatial units with stress retention and fluid inclusion characteristics to obtain the reservoir space block delineation map set are as follows: S501: Obtain the amplitude attenuation rate and the pressure inversion value within the block, and conduct distribution analysis. Through the calculation of the amplitude attenuation rate and the pressure inversion value, obtain the amplitude attenuation and pressure inversion value intervals within the block, and respectively allocate the amplitude attenuation rate and the pressure inversion value of the block into intervals to obtain the amplitude attenuation and pressure inversion value distribution data corresponding to the block; S502: According to the distribution of the amplitude attenuation rate and the pressure inversion value within the block, screen the spatial units with stress retention and fluid inclusion characteristics, and use the amplitude attenuation rate and the pressure inversion value within the block for analysis to obtain the stress retention characteristic block and the fluid inclusion characteristic block; S503: Based on the stress retention characteristic block and the fluid inclusion characteristic block, combined with the pore-fracture coordination response information, aggregate the spatial units that meet the stress retention and fluid inclusion characteristics to obtain the reservoir space block delineation map set.

10. A shale gas reservoir space detection system, characterized in that, According to any one of claims 1-9, it includes the following steps: The reservoir layer section division module calculates the unilateral average gradient value of the longitudinal change of the pores based on the pore size and the layer position point data, judges the stability of the layer section structure, selects the continuous abnormal points as the interlayer transition zone, and divides the physical property independent layer sections according to the position of the transition zone to obtain the reservoir layer division boundary division result; The wave velocity correction analysis module calculates the difference between the original wave velocity value of the layer section and the pore mean value based on the reservoir layer section division result, adjusts the wave velocity response intensity in combination with the local stress interference data, and obtains the in-reservoir structure response difference index group by comparing the offset degree of the stress deformation amplitude before and after correction; The reservoir deformation trend analysis module selects the index offset direction within the selected area as the judgment reference based on the in-reservoir structure response difference index group, demarcates the continuous change trend, sets the direction angle threshold, and extracts the graphic group in which the direction angle between adjacent units deviates from the threshold to obtain the main control path band of the reservoir deformation trend; The pore-fracture coordination response screening module conducts single-point comparison of the angle between the pore axis direction and the main fracture distribution direction within the area where the main control path band of the reservoir deformation trend passes through, screens the contact points whose direction differences deviate from the set angle interval, and obtains the pore-fracture coordination response apparent block set; Based on the pore-fracture collaborative response manifestation block set, the reservoir space block identification module combines the amplitude attenuation rate and the combined distribution of pressure inversion value intervals within the block to screen out spatial units with stress retention and fluid inclusion characteristics, and obtains a reservoir space block delineation map.

Images