A method for determining the width and intensity of tectonic stress disturbance zones on both sides of vertical faults

Abstract

The present application discloses a method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault, aiming to solve the problem in the prior art that it is difficult to accurately determine the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault. By statistically analyzing the strike and plane extension length of the vertical faults in the study area, the static rock mechanics parameters are solved; a geomechanical model and a mathematical model related to the sliding of the vertical faults along the strike are constructed, and numerical simulation experiments of the tectonic stress field are carried out to obtain the maximum horizontal principal stress intensity at each position; the average value of the maximum horizontal principal stress in each interval is statistically calculated; a contour map is drawn, and the width and intensity of the stress disturbance zone are judged based on the contour map. The present application not only gives the precise numerical value of the width of the tectonic stress disturbance zone on both sides of the fault under the specific regional tectonic stress background, but also gives the precise numerical value of the tectonic stress disturbance intensity, which significantly improves the judgment accuracy.

Description

Technical Field

[0001] The present application relates to the technical field of reservoir geomechanics in geology, and in particular to a method for determining the width and intensity of structural stress disturbance zones on both sides of a vertical fault. Background Art

[0002] Geostress is a type of stress found in rock masses (Jing Feng, Sheng Qian, Zhang Yonghui, et al., 2011. Research Progress on In-situ Geostress Measurement and Geostress Field Analysis in my country [J]. Rock and Soil Mechanics, 32(S2):51-58). It includes both paleo- and present-day geostress. Tectonic stress is crucial in both paleo- and present-day geostress and plays a crucial role in all aspects of oil and gas field exploration and development. Tectonic stress is significantly affected by faults (HUDSON JA, HARRISON JP, 2000. Engineering rock mechanics in an introduction to the principles [M]. Oxford: Elsevier). This is manifested by significant changes in the intensity and direction of tectonic stress in strata near faults compared to those in ordinary sedimentary formations, a phenomenon known as tectonic stress perturbation. Recent important oil and gas discoveries in the Shunbei area of the Tarim Basin indicate that the study of strike-slip of pre-existing near-vertical faults and their associated tectonic fracture zones is of great significance for petroleum geology. Tectonic rupture caused by paleotectonic stress is an important factor controlling the activity of vertical faults and the development of tectonic fractures (Ding, WL, Fan, TL, Yu, BS, Huang, XB & Liu, C., 2012. Ordovician carbonate reservoir fracture characteristics and fracture distribution forecasting in the Tazhong area of Tarim Basin, Northwest China, Journal of Petroleum Science and Engineering, 86-87, 62-70; Huang, L., Liu, C., He, F., Jia, H., Zhou, Y., Wang, Z., Wang, J., Liu, Y. & Li, X., 2022. Strike-slip deformation characteristics offault in craton Basin, Journal of Northwest University. Natural Science Edition,52(6),930-942).It is crucial to study the control of paleotectonic stress disturbance on the distribution characteristics of tectonic fracture zones along vertical faults. The width and intensity of the paleotectonic stress disturbance zone caused by the strike-slip faulting of vertical faults have a strong control on the width and development degree of the tectonic fracture zone (Li, Y., Ding, W., Han, J., Chen, X., Huang, C., Li, J. & Ding, S., 2024. Quantitative Prediction of the Development and Opening Sequence of Fractures in an Ultradeep Carbonate Reservoir: A Case Study of the Middle Ordovician in the Shunnan Area, Tarim Basin, China, SPE Journal, 29(6), 3091-3117; Li, Y., Ding, W., Han, J., Chen, X., Huang, C., Li, J. & Ding, S., 2024. Tectonic fractures induced by strike-slip faulting in intracratonic ultradeep carbonate rocks: Insights from the finite element method and self-adaptive constraints computational model for boundary conditions, GSA Bulletin, 136(11-12), 4512-4540). Therefore, the study of the paleo-tectonic stress disturbance zone caused by the along-strike sliding of vertical faults will have very important significance for oil and gas geology and development.

[0003] At present, the research method of tectonic stress disturbance zone is: after setting the existing fault as an abnormal zone of rock mechanical properties, the geomechanical model and mathematical model are constructed and three-dimensional finite element numerical simulation is carried out. According to the numerical simulation results, the width of the stress disturbance zone near the fault, the changes in the magnitude and direction of the ground stress and their influencing factors are qualitatively and quantitatively analyzed (Weng Jianqiao, Zeng Lianbo, Lv Wenya, Liu Qi & Zu Kewei, 2020. The width of the ground stress disturbance zone near the fault and its influencing factors, Journal of Geomechanics, 26(01), 39-47; Wang Yuanyuan, Huang Cheng, Gong Wei, et al. Analysis of Silurian fault characteristics and stress field disturbance and well location optimization in Shunbei area of Tazhong [J]. Experimental Petroleum Geology, 2024, 46(04): 674- 682; Ding Wenlong, Li Yuntao, Han Jun, et al. High-precision tectonic stress field simulation and fracture multi-parameter distribution prediction method of carbonate reservoirs and its application [J]. Petroleum & Natural Gas Geology, 2024, 45(03): 827-851; Shen Jie, Xu Hao, Deng Hucheng, et al. Study on the distribution characteristics and disturbance mechanism of ground stress field in complex fault areas - a case study of the Upper Paleozoic in Dingbei area of Ordos Basin [J / OL]. Chinese Geology, 1-19 [2025-01-22]; Lei Xianghui, Xu Dongsheng, Shi Leiting, et al. Ground stress disturbance characteristics and inversion method of fault-developed oil and gas reservoirs [J]. Journal of Liaoning University of Engineering and Technology (Natural Science Edition), 2023, 42(03): 293-300).

[0004] The second method is to combine a variety of methods (hydraulic fracturing, acoustic emission experimental measurement, wellbore collapse method, logging calculation, imaging logging method, wellbore size method and wave velocity anisotropy method, etc.) to accurately interpret the current ground stress of a single well near the fault, and analyze the disturbance effect of the ground stress based on this (Zhang Jiawei, Li Ruixue, Deng Hucheng, et al. Disturbance characteristics of the ground stress field of tight sandstone reservoirs in ultra-deep thrust-nappe structures - taking the Cretaceous reservoirs in the Bozi-Dabei area of the Tarim Basin as an example [J]. Experimental Petroleum Geology, 2024, 46(04): 760-774). Summary of the Invention

[0005] The present application provides a method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault, aiming to solve the problem in the prior art that it is difficult to accurately determine the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault based on the characteristics of the regional tectonic stress field and the vertical fault.

[0006] In a first aspect, a method for determining the width and intensity of a tectonic stress disturbance zone on both sides of a vertical fault comprises:

[0007] Step S1: determine the study area and target strata, calculate the strike and horizontal extension length of vertical faults in the study area, determine the key period of fault activity and the paleo-tectonic stress field environment; further determine the direction and intensity of the regional maximum horizontal principal stress, the intensity of the minimum horizontal principal stress, and the vertical stress intensity during the key period;

[0008] Step S2, solving static rock mechanical parameters based on core, conventional logging and array acoustic logging data of the target layer;

[0009] Step S3, constructing a geomechanical model and a mathematical model related to the sliding of the vertical fault along the strike, and conducting a numerical simulation test of the tectonic stress field to obtain the maximum horizontal principal stress intensity at each location;

[0010] Step S4, in the simulation experiment, calculate the distance from all nodes to the vertical fracture, divide the intervals according to the distance and research accuracy requirements, set the interval length or the number of intervals, and calculate the average value of the horizontal maximum principal stress in each interval; use the median of the interval to represent the overall distance from the nodes in the interval to the vertical fracture;

[0011] Step S5, selecting the nearest intervals on both sides of the vertical fault, and determining the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault;

[0012] Step S6: draw a contour map of the tectonic stress disturbance zone width and disturbance intensity as the boundary conditions change, and determine the stress disturbance zone width and disturbance intensity based on the contour map.

[0013] In the above solution, further optionally, step S2 includes:

[0014] Calculating dynamic rock mechanics parameters based on conventional logging and array acoustic logging data, wherein the dynamic rock mechanics parameters include dynamic Young's modulus and dynamic Poisson's ratio;

[0015] Obtaining static rock mechanics parameters from rock mechanics experiments on core samples, wherein the static rock mechanics parameters include static Young's modulus and static Poisson's ratio parameters;

[0016] The dynamic rock mechanics parameters and static rock mechanics parameters of the target strata in the study area are fitted to obtain the conversion model of dynamic and static rock mechanics parameters:

[0017]

[0018] Where, E represents the static Young's modulus; μ represents the static Poisson's ratio; E d represents the dynamic Young's modulus; μ d represents the dynamic Poisson's ratio.

[0019] In the above solution, further optionally, the calculating of dynamic rock mechanical parameters based on conventional logging and array acoustic logging data includes:

[0020] Based on conventional logging data, a conversion model between acoustic time difference AC and compressional time difference DTC is constructed; based on array acoustic logging data, a conversion model between compressional time difference DTC and shear time difference DTS is constructed; using the two conversion models, DTC and DTS are calculated.

[0021] The calculation formulas for the dynamic Young's modulus and dynamic Poisson's ratio of the target formation during drilling are as follows:

[0022]

[0023] Among them, E d Represents dynamic Young's modulus, in GPa; μ d is the dynamic Poisson's ratio; Δt p is the longitudinal wave time difference DTC, in μs / ft; Δt s is the shear wave time difference DTS, with the unit of μs / ft; ρ is the density, with the unit of g / cm3.

[0024] In the above solution, optionally, step S3 includes:

[0025] Based on the geological characteristics of the study area, a three-dimensional geomechanical model was constructed. The model included the location, strike, and extension of the vertical fault, as well as the geometric characteristics of the target strata. The vertical fault was set as a discontinuity with a certain friction coefficient to simulate the relative sliding of the strata on both sides of the fault.

[0026] The geomechanical model is converted into a mathematical model, and meshing is performed using finite element analysis software to divide the model into several units and nodes;

[0027] Assign the static rock mechanics parameters obtained from step S2 to each unit in the mathematical model; set the boundary conditions and constraints of the mathematical model according to the research accuracy requirements; conduct numerical simulation tests of the tectonic stress field, and obtain the horizontal maximum principal stress intensity at each location of the model from the test results;

[0028] Compare the simulation results with actual geological data to verify the accuracy and reliability of the model.

[0029] In the above solution, optionally, step S4 includes:

[0030] Adjust the coordinate system of the model to make it consistent with the direction of the regional stress field, calculate the distance from each node in the model to the vertical fracture; count the maximum value D of the distance from the node to the vertical fracture MAX and minimum value D MIN ; Set the interval length D according to the research accuracy requirementsINT ; The range of each interval is:

[0031] [D MIN +(i-1) / n*(D MAX -D MIN ), D MIN +i / n*(D MAX -D MIN )]

[0032] Where n represents the number of preset intervals; i represents the interval number;

[0033] For each interval, record the maximum horizontal principal stress of each node in the interval, and after multiple numerical simulation tests of the stress field, calculate the average value of the maximum horizontal principal stress in the interval;

[0034] The median of the interval is used to represent the overall distance between the nodes in the interval and the vertical fracture; the median of the interval is:

[0035] D MIN +(i-0.5) / m*(D MAX -D MIN ).

[0036] In the above solution, optionally, step S5 includes:

[0037] Get the median of each interval from step S4. For all the interval medians less than or equal to 0, find the interval median with the smallest absolute value, which is recorded as D LEFT-MIN The corresponding average value of the maximum horizontal principal stress is recorded as S LEFT-MIN ; For all interval medians greater than or equal to 0, find the interval median with the smallest absolute value, recorded as D RIGHT-MIN The corresponding average value of the maximum horizontal stress is recorded as S RIGHT-MIN If D LEFT-MIN =D RIGHT-MIN =0, then the two intervals are actually the same interval, and the same interval is still used in subsequent analysis;

[0038] For each interval, check the average value of the horizontal maximum principal stress S H (k,i) Whether one of the following conditions is met:

[0039] [S H (k,i)-S H (k,i-1)]·[S H (k,i+1)-S H (k,i)]<0

[0040]

[0041] Among them, S *Indicates the horizontal maximum principal stress threshold; S H represents the average value of the horizontal maximum principal stress simulation results; i represents the interval number; k represents the number of the simulation experiment;

[0042] For the median of the interval less than or equal to 0, find the interval that meets the above conditions and record the median of the interval as D LEFT-MAX (k), the corresponding average value of the maximum horizontal principal stress is recorded as S LEFT-MAX (k); For the median of the interval greater than or equal to 0, find the interval that meets the above conditions and record the median of the interval as D RIGHT-MAX (k), the corresponding stress average value is recorded as S RIGHT-MAX (k);

[0043] The intervals representing the maximum horizontal principal stress mutation on both sides of the vertical fault in all simulation tests were obtained, and the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault were calculated according to the following set of equations:

[0044]

[0045] When the median values of the interval are less than or equal to 0, the width and intensity of the tectonic stress disturbance zone are △D LEFT (k) and △S LEFT (k); When the median values of the interval are greater than or equal to 0, the width and intensity of the tectonic stress disturbance zone are △D RIGHT (k) and △S RIGHT (k);

[0046] △D LEFT (k) and ΔD RIGHT (k) are divided by the length of the vertical fault preset in the model to obtain the tectonic stress disturbance zone width coefficient △D' LEFT (k) and △D' RIGHT (k).

[0047] In the above solution, optionally, step S6 includes:

[0048] Plotting θ at different angles, △D' LEFT (k), △S LEFT (k), △D' RIGHT (k) and △S RIGHT (k) plane contour map, the horizontal axis of the plane contour map is the ratio of the horizontal maximum principal stress to the vertical stress σ H / σ V , the vertical axis is the ratio of the horizontal minimum principal stress to the horizontal maximum principal stress σ h / σ H ; where θ represents the angle between the strike of the vertical fault and the direction of the regional maximum horizontal principal stress;

[0049] According to the location of the regional stress source, the width and intensity of the tectonic stress disturbance zone on this side are determined as △D RIGHT (k) and △S RIGHT (k); the other side is △D LEFT (k) and △S LEFT (k);

[0050] Observe the trends in the contour map to determine the changing patterns of the width and intensity of the stress disturbance zone and identify the boundary of the stress disturbance zone, that is, the location of the stress mutation. Based on the contour map, determine the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault.

[0051] Calculate the disturbance band width coefficient △D' LEFT (k) and △D' RIGHT (k), and the disturbance intensity △S LEFT (k) and △S RIGHT (k).

[0052] In the above solution, further optionally, determining the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault includes:

[0053] The different θ values are numbered as α1 to α t , where θ represents the angle between the strike of the vertical fault and the direction of the regional horizontal maximum principal stress;

[0054] Get and α j The corresponding contour map, and obtain the △D' corresponding to the vertical fracture Fi from the contour map LEFT (Fi,α j ), △S LEFT (Fi,α j ),△D' RIGHT (Fi,α j ) and △S RIGHT (Fi,α j );

[0055] Get and α j+1 The corresponding contour map, and obtain the △D' corresponding to the vertical fracture Fi from the contour map LEFT (Fi,α j+1 ), △S LEFT (Fi,α j+1 ),△D' RIGHT (Fi,α j+1 ) and △S RIGHT (Fi,α j+1 );

[0056] Calculate the actual △D' corresponding to the vertical fracture Fi LEFT (Fi,θ i ), △SLEFT (Fi,θ i ),△D' RIGHT (Fi,θ i ) and △S RIGHT (Fi,θ i ), the calculation formula is:

[0057]

[0058] After calculating the actual △D of all fractures LEFT (Fi,θ i ) and △D RIGHT (Fi,θ i ), the actual width of the tectonic stress disturbance zone needs to be calculated based on the actual length of the fault. LEFT (Fi,θ i ) and △D' RIGHT (Fi,θ i ) is the disturbance band width coefficient.

[0059] Compared with the prior art, this application has at least the following beneficial effects:

[0060] Based on further analysis and research on existing technical problems, this application recognizes that compared with the existing technology that usually regards the tectonic stress disturbance zone near the fault as a whole, the method proposed in this application takes into account the width and intensity of the tectonic stress disturbance zone on different sides of the fault (determined by the regional horizontal maximum principal stress), and is therefore more advantageous in describing the heterogeneity of the tectonic stress state near the fault.

[0061] Compared with the existing technology that usually gives the range of the width of the tectonic stress disturbance zone near the fault, this application not only gives the precise numerical value of the width of the tectonic stress disturbance zone on both sides of the fault under the specific regional tectonic stress background, but also gives the precise numerical value of the tectonic stress disturbance intensity, which reduces the error caused by human judgment and significantly improves the judgment accuracy.

[0062] Compared with the existing technology that usually sets the pre-existing fault as a zone of abnormal rock mechanical properties to conduct tectonic stress disturbance zone width analysis, the method proposed in this application sets the pre-existing vertical fault as a discontinuous surface with a certain friction coefficient, resulting in the strata on both sides of the vertical fault to slide relative to each other under the background of the regional tectonic stress field, which is closer to the actual deformation characteristics of the strata on both sides when the vertical fault slides along its direction. Therefore, compared with the existing technology, the method proposed in this application is more comprehensive and in-depth in considering the causal mechanism of tectonic stress, and has more practical significance and application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1A flow chart of a method for determining the width and intensity of structural stress disturbance zones on both sides of a vertical fault provided in one embodiment of the present application.

[0064] Figure 2 An embodiment of the present application provides a linear conversion mathematical model of AC and DTC obtained from the AC and DTC data of the target strata in the study area.

[0065] Figure 3 An embodiment of the present application provides a linear conversion mathematical model of DTC and DTS obtained from the DTC and DTS data of the target strata in the study area.

[0066] Figure 4 An embodiment of the present application provides a geomechanical model size setting, a vertical fracture setting, a boundary condition setting, and 9 different geomechanical models corresponding to 9 different values of the angle θ between the vertical fracture and one side of the geomechanical model, and a corresponding mathematical model.

[0067] Figure 5 This is the numerical simulation result of the horizontal maximum principal stress of the top interface of the model provided in one embodiment of the present application, with the boundary condition being σ H / σ V Equal to 3, σ h / σ H equal to 0.2, θ equal to 30°.

[0068] Figure 6 This is the numerical simulation result of the horizontal maximum principal stress of the top interface of the model provided in one embodiment of the present application, with the boundary condition being σ H / σ V Equal to 2, σ h / σ H equal to 0.6, θ equal to 60°.

[0069] Figure 7 The calculation results of the distance between the node and the vertical fracture of the top interface of the model provided in one embodiment of the present application are as follows: the boundary condition is that θ is equal to 30°.

[0070] Figure 8 The calculation result of the distance between the node and the vertical fracture of the top interface of the model provided in one embodiment of the present application, the boundary condition is that θ is equal to 60°.

[0071] Figure 9 The relationship between the average value of the maximum horizontal principal stress in each interval and the median value of the interval is provided in an embodiment of the present application, and the boundary condition is σ H / σ V is equal to 1 and θ is equal to 45°, where curves 1 to 9 correspond to σ h / σ HEqual to 0.1~0.9 with an interval of 0.1.

[0072] Figure 10 The relationship between the average value of the maximum horizontal principal stress in each interval and the median value of the interval is provided in an embodiment of the present application, and the boundary condition is σ h / σ H is equal to 0.3 and θ is equal to 45°, where curves 1 to 9 correspond to σ H / σ V Equal to 1 to 3 with an interval of 0.25.

[0073] Figure 11 When θ is equal to 50°, ΔD is provided in one embodiment of the present application. LEFT About σ H / σ V and σ h / σ H Plane contour map of .

[0074] Figure 12 When θ is equal to 50°, ΔS is provided in one embodiment of the present application. LEFT About σ H / σ V and σ h / σ H Plane contour map of .

[0075] Figure 13 When θ is equal to 50°, ΔD is provided in one embodiment of the present application. RIGHT About σ H / σ V and σ h / σ H Plane contour map of .

[0076] Figure 14 When θ is equal to 50°, ΔS is provided in one embodiment of the present application. RIGHT About σ H / σ V and σ h / σ H Plane contour map of .

[0077] Figure 15 When θ is equal to 60°, ΔD is provided in one embodiment of the present application. LEFT About σ H / σ V and σ h / σ H Plane contour map of .

[0078] Figure 16 When θ is equal to 60°, ΔS is provided in one embodiment of the present application. LEFT About σ H / σ V and σ h / σ H Plane contour map of .

[0079] Figure 17 When θ is equal to 60°, ΔD is provided in one embodiment of the present application. RIGHT About σ H / σ V and σ h / σ H Plane contour map of .

[0080] Figure 18 When θ is equal to 60°, ΔS is provided in one embodiment of the present application. RIGHT About σ H / σ V and σ h / σ H Plane contour map of . DETAILED DESCRIPTION

[0081] In order to make the purpose, technical solutions and advantages of this application more clear, the following further describes this application in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0082] In the description of this application: unless otherwise specified, "a plurality of" means two or more. Expressions such as "include", "comprising", "having" and the like also mean "not limited to" (certain units, components, materials, steps, etc.).

[0083] Based on further analysis and research of existing technical problems, this application recognizes that the methods mentioned in the background technology have the following three problems:

[0084] Question 1: Usually, the pre-existing fault is set as the zone of abnormal rock mechanical properties to carry out the tectonic stress disturbance zone width analysis. Such a setting is still applicable in the analysis of current ground stress disturbance characteristics, but it is not applicable in the analysis of ancient tectonic stress disturbance zones because the strata on both sides of the fault do not produce relative sliding, which does not conform to the actual geological conditions.

[0085] Question 2: The tectonic stress disturbance zone near the fault is usually compared as a whole, but the differences in the ancient tectonic stress disturbance characteristics on both sides of the fault are not clearly portrayed. In the context of regional tectonic movement, there must be significant differences between the tectonic stress disturbance characteristics on the side of the fault closer to the regional stress source and the tectonic stress disturbance characteristics on the side of the fault farther from the regional stress source.

[0086] Question 3: The width of the tectonic stress disturbance zone near a fault is usually given as a range, rather than a precise coefficient or ratio. This increases the error caused by human judgment and reduces the research precision. Therefore, quantitative research on the width and intensity of the paleotectonic stress disturbance zone associated with the along-strike slip of the vertical fault is urgently needed. In this process, it is also necessary to pay attention to the differences in the width and intensity of the paleotectonic stress disturbance zone on both sides of the vertical fault, and the strata near the fault cannot be simply considered as a whole.

[0087] The purpose of this application is to provide a method for determining the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault, so as to solve the problem in the existing technology that it is difficult to quickly determine the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault based on the characteristics of the regional tectonic stress field and the vertical fault, and to promote the exploration and development process of oil and gas fields and reservoir geomechanics research.

[0088] In one embodiment of the present application, a method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault is provided, comprising:

[0089] Step S1, determine the study area and target strata, count the strike and plane extension length of the vertical faults that cut through the target strata in the study area, and determine the key tectonic geological history period when the vertical faults slipped along the strike and the corresponding paleo-tectonic stress field environment.

[0090] Step S11, according to the actual research requirements of oil and gas exploration and development, determine the study area and target strata, and count the strike and plane extension length of vertical faults in the study area that penetrate the target strata.

[0091] Based on the actual needs of oil and gas exploration and development research, the study area and target strata are selected, and the lithology of the target strata and the depth distribution characteristics of the top and bottom interfaces within the study area are clarified. Based on detailed interpretation of 2D or 3D seismic data, all vertical faults in the study area that intersect the target strata are identified, numbered, and their strike (expressed by the fracture azimuth, ranging from 0° to 180°) and horizontal extension are determined. These vertical faults must completely intersect the target strata, and their intersections with the top and bottom interfaces of the target strata must be continuous and clearly visible in the plane.

[0092] In one embodiment of the present application, China's A Basin is selected as the study area, and the important oil and gas producing layer B is selected as the target stratum. When conducting a comprehensive stratigraphic-structural interpretation based on three-dimensional seismic data in Basin A, the conversion between two-way travel time and actual depth was not performed. Therefore, the depth based on time was used, and the top interface depth of the producing layer B was obtained to be approximately 2400ms, the bottom interface depth to be approximately 2600ms, and the total thickness to be approximately 200ms. Based on the detailed interpretation of two-dimensional or three-dimensional seismic data, the planar distribution characteristics of a total of 13 vertical faults that interrupted the producing layer B in Basin A were determined. These faults all interrupted the producing layer B, and the planar traces were nearly straight and clearly identifiable. The azimuths and planar extension lengths of these faults are shown in Table 1.

[0093] Table 1

[0094]

[0095]

[0096] Step S12, determining the key tectonic geological historical periods and corresponding paleo-tectonic stress field environments during which each vertical fault in the study area slipped along its strike.

[0097] Based on the relevant literature on the tectonic evolution analysis of the study area and the tectonic evolution history mapping technology, the key tectonic geological history period related to the sliding of each vertical fault in the study area was determined. Each fault should be analyzed separately because vertical faults that penetrate the same target layer are not necessarily formed in the same tectonic geological history period. After determining the key tectonic geological history period related to the sliding of each vertical fault in the study area, the regional horizontal maximum principal stress (σ H ) direction, σ H Strength, regional horizontal minimum principal stress (σ h ) strength and vertical stress (σ V ) intensity. It should be noted that σ H The direction of is from the stress source to the fracture direction, for example, σ H Pointing from the due east to the due west, the stress source is on the due east side of the vertical fracture, then σ H The direction is -90°; if σ H From the southwest to the northeast, the stress source is on the southwest side of the vertical fracture, then σ H The direction is 45°; if σ H From the southeast direction to the northwest direction, the stress source is on the southeast side of the vertical fault, then σ H The direction should be between -45° and -90°.

[0098] In one embodiment of the present application, the key structural geological history periods associated with the sliding of faults F1 to F6 along the strike in Table 1 are consistent and are designated as period I; the key structural geological history periods associated with the sliding of faults F7 to F9 along the strike are consistent and are designated as period II; the key structural geological history periods associated with the sliding of faults F10 to F13 along the strike are consistent and are designated as period III. The σ corresponding to periods I, II, and III is H direction, σ H Strength, σ h Intensity and σ V The strength is shown in Table 2.

[0099] Table 2

[0100]

[0101] Step S2: Ascertain the core, conventional logging, array acoustic logging and other data of the target stratum in the study area, and solve the static rock mechanical parameters of the target stratum in the study area.

[0102] Step S21: Ascertain the core, conventional logging, and array acoustic logging data for the target formation within the study area. Detailed statistics are collected for the core, conventional logging, and array acoustic logging data for each well drilled through the target formation within the study area. If the target formation within the study area has no well or core data, or only conventional logging data, it is considered a Type I study area. If the target formation within the study area has core data but no array acoustic logging data, it is considered a Type II study area. If the target formation within the study area has core, conventional logging, and array acoustic logging data, it is considered a Type III study area. If the target formation within the study area has conventional logging and array acoustic logging data but no core data, it is considered a Type IV study area.

[0103] In one embodiment of the present application, formation B in basin A has core data, conventional logging data, and array acoustic logging data, and therefore belongs to a type III study area.

[0104] S22 solves the static rock mechanics parameters of the target strata in the study area.

[0105] For the Type II study area mentioned in step S21, rock mechanics tests are carried out on core samples to obtain parameters such as the static Young's modulus, static Poisson's ratio and density of the rock. For each parameter, the average value of the corresponding static parameters of all samples is taken. For the type III study area mentioned in step S21, the dynamic rock mechanical parameters are first calculated based on conventional logging and array acoustic logging, and then the dynamic rock mechanical parameters are converted into static rock mechanical parameters based on the static Young's modulus and static Poisson's ratio obtained from the rock mechanical test of the core sample. Specifically, the acoustic time difference (AC) logging data and the longitudinal time difference (DTC) and shear time difference (DTS) logging data in the conventional logging data are time difference data and are considered to have good correlation. Therefore, the conversion model of AC and DTC is first constructed, and then the conversion model of DTC and DTS is constructed based on the array acoustic logging data; Based on these two conversion models, DTC and DTS can be calculated, and the dynamic Young's modulus and dynamic Poisson's ratio of all wells in the target formation are calculated according to the following formulas (1) and (2):

[0106]

[0107] Among them, E d is the dynamic Young's modulus, in GPa; μ d is the dynamic Poisson's ratio; Δt p is the longitudinal wave time difference DTC, in μs / ft; Δt s is the shear time difference DTS, μs / ft; ρ is the density, unit is g / cm 3 .

[0108] Dynamic rock mechanical parameters at corresponding locations are obtained based on the drilling location and depth of the core samples used for rock mechanics testing. Combined with the static rock mechanical parameters obtained from the rock mechanics tests, a conversion model between dynamic and static rock mechanical parameters is constructed. This conversion model can be used to convert the dynamic rock mechanical parameter calculations for all wells drilled in the study area into static rock mechanical parameters in the target formation. The average static rock mechanical parameters for all wells drilled in the study area in the target formation are then calculated.

[0109] For the Class I or Class IV study area mentioned in step S21, based on the identification of the main lithology of the target strata in the study area, the static rock mechanics parameters of the corresponding lithology in the nearby area are obtained by consulting the literature.

[0110] In one embodiment of the present application, Basin A belongs to the Type III study area, so the AC and DTC conversion model is constructed based on the AC and DTC data of all wells in the study area, such as Figure 2 As shown, the linear fitting method is used, and the goodness of fit R 2is 0.66. Then, a DTC and DTS conversion model is constructed based on the DTC and DTS data of all wells in the study area, as shown in Figure 3 As shown, linear fitting is also used, and the goodness of fit R 2 is 0.86. Figure 2 and Figure 3 The conversion model shown can convert all AC data in the study area into DTC and DTS data. Dynamic Young's modulus and dynamic Poisson's ratio can then be calculated according to equations (1) and (2). Static rock mechanical parameters can be obtained from rock mechanical experiments on rock samples from the target strata in the study area. The corresponding relationship between these static rock mechanical parameters and the dynamic rock mechanical parameters calculated for the corresponding wells and depths is shown in Table 3.

[0111] Table 3

[0112] Sample number Dynamic Young's modulus / GPa Static Young's modulus / GPa Dynamic Poisson's ratio Static Poisson's ratio 1 75.807 47.151 0.341 0.238 2 75.997 46.749 0.340 0.212 3 76.288 43.695 0.344 0.210 4 73.650 39.600 0.345 0.240 5 74.662 41.265 0.342 0.226 6 64.259 34.800 0.307 0.205 7 67.457 34.791 0.302 0.217 8 70.906 42.865 0.321 0.238 9 71.300 41.329 0.309 0.205 10 67.752 42.304 0.320 0.268

[0113] Fitting the dynamic and static rock mechanics parameters in Table 3 can obtain the conversion models of the dynamic and static rock mechanics parameters of the target layer, as shown in formulas (3) and (4), respectively:

[0114]

[0115] The goodness of fit of formulas (3) and (4) is 0.62 and 0.63, respectively. Based on formulas (3) and (4) and the calculation results of the dynamic rock mechanical parameters of the target strata in the study area, the static rock mechanical parameters of the target strata in the study area can be calculated. The calculation results show that the average value of the static Young's modulus is 39.15 GPa, the average value of the static Poisson's ratio is 0.232, and the average value of the rock density is 2.683 g / cm 3 .

[0116] Step S3: constructing a geomechanical model and a mathematical model related to the strike-side sliding of the vertical fault, and conducting a numerical simulation test of the tectonic stress field to obtain the maximum horizontal principal stress intensity at each position of the model from the test results.

[0117] Step S31 , setting the size of the regular tetrahedron geomechanical model, setting the vertical fracture as a contact pair in the finite element analysis, and setting a number of different regular tetrahedron geomechanical models according to the plane extension orientation of the vertical fracture.

[0118] The geomechanical model is set as a regular square prism, and the top and bottom interfaces of the regular square prism are congruent squares. In order to facilitate the subsequent mathematical model construction and result analysis, the side lengths of the top and bottom interfaces of the regular square prism are set to 6000m, and the height H of the regular square prism is set to 5000m. The vertical fault is set as a cavity enclosed by four vertical planes with a height of H. These four planes can be divided into two groups. The planes belonging to different groups are perpendicular to each other, and each group consists of two planes: the spacing between the first group of planes is set to 5m, and the plane extension length is 5000m; the spacing between the second group of planes is 5000m, and the plane extension length is 5m. The top and bottom surfaces of the existing basement fault area are both rectangles with a length of 5000m and a width of 5m. The height of the area is H, and the volume is 1.25×10 8 m 3 The pre-existing basement fault area runs through the geomechanical model from top to bottom, and the center of the intersection line of the area and the top surface of the geomechanical model coincides with the center of the top surface of the geomechanical model, and the center of the intersection line of the fault and the bottom surface of the geomechanical model coincides with the center of the bottom surface of the geomechanical model.

[0119] The contact pair concept, used in contact analysis within the statics module of finite element analysis, is introduced and applied to the configuration of pre-existing fractures in 3D geomechanics models. A contact pair consists of two planes (or surfaces): a contact surface and a target surface. The contact surface is the surface undergoing passive deformation and is therefore defined as the fault plane on the side undergoing passive deformation in the regional tectonic stress field. The target surface is the surface undergoing active deformation and is therefore defined as the fault plane on the side undergoing active deformation in the regional tectonic stress field. Based on the regional stress direction associated with the target strata and the activity of fault F1 in the study area, the section closest to the regional stress direction was set as the target surface, and the section away from the regional stress direction was set as the contact surface. Both sections had a planar extension length of 5000 m, a height of 5000 m, and a spacing of 5 m. The friction coefficient between the contact surface and the target surface was set based on the friction coefficient of the vertical fault in the strata with similar lithology to the target strata. Under the influence of regional tectonic stress, when the stress intensity at the section exceeds the stress intensity threshold for relative sliding between the contact surface and the target surface, relative sliding can occur between the contact surface and the target surface, thereby achieving relative displacement of the strata on both sides of the pre-existing fault. After completing the geomechanical model size setting and the vertical fault contact pair setting, the angle between the vertical fault and any side surface of the geomechanical model was set to θ. Based on the research accuracy requirements, θ was selected from n values within the range of 0 to 90°. i (1≤i≤n), each θ i Each corresponds to a different geomechanical model, which can be used to study the differences in deformation characteristics of the strata on both sides of the fault when the angle between the vertical fault F1 and the regional horizontal stress is different.

[0120] In one embodiment of the present application, the three-dimensional geomechanical model is set as a regular square prism, the length and width of the top surface are both 6000m, and the height H of the geomechanical model is equal to 5000m. The vertical fault is set as a cavity enclosed by four vertical planes with a height of 5000m. These four planes can be divided into two groups, and the planes belonging to different groups are perpendicular to each other. Each group consists of two planes: the first group of planes has a spacing of 5m and a plane extension length of 5000m; the second group of planes has a spacing of 5000m and a plane extension length of 5m. The top and bottom surfaces of the pre-existing basement fault area are both rectangles with a length of 5000m and a width of 5m. The height of the area is 5000m and the volume is 1.25×10 8 m 3 The pre-existing basement fault area runs through the geomechanical model from top to bottom, and the center of the intersection line of the area and the top surface of the geomechanical model coincides with the center of the top surface of the geomechanical model, and the center of the intersection line of the fault and the bottom surface of the geomechanical model coincides with the center of the bottom surface of the geomechanical model.

[0121] In one embodiment of the present application, the direction of the horizontal maximum principal stress is fixed from south to north, so the south section of the vertical fault is set as the target surface, and the north section is set as the contact surface. The plane extension length of both is 5000m, the height is 5000m, and the spacing is 5m. The friction coefficient of the pre-existing fault is considered to be 0.15, so the friction coefficient between the contact surface and the target surface is set to 0.15. Under the action of regional tectonic stress, when the stress intensity at the section exceeds the stress intensity threshold for relative sliding between the contact surface and the target surface, relative sliding can occur between the contact surface and the target surface, and the relative displacement of the strata on both sides of the pre-existing fault can be achieved. θ is set to 10°, 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80°, corresponding to 9 different three-dimensional geomechanical models.

[0122] In step S32, the geomechanical model constructed in step S31 is meshed using tetrahedral elements in finite element analysis software to obtain a corresponding mathematical model. The mesh density at the top interface of the model is higher than that at other locations in the model. After the mathematical model is constructed, the coordinates of all elements and nodes of the geomechanical model are exported from the finite element analysis software.

[0123] The geomechanical model is meshed using triangular (2D) or tetrahedral (3D) elements, dividing the model into a series of elements and nodes. The top surface of the geomechanical model has higher meshing accuracy than the bottom and four side surfaces. This is because the top surface of the model is horizontal and serves as the outer surface of the model, which can more intuitively reflect the planar heterogeneity of horizontal stress. This approach significantly improves accuracy if contour maps of the horizontal maximum principal stress at the top surface of the model are required. After the mathematical model is constructed, the spatial coordinates (E_X, E_Y, E_Z) of all elements and the spatial coordinates (N_X, N_Y, N_Z) of all nodes in the geomechanical model are exported from the finite element analysis software.

[0124] In one embodiment of the present application, tetrahedral elements are used to mesh the three-dimensional geomechanical model, dividing the geological model into a series of elements and nodes. ANSYS computer software is used to mesh the geomechanical model, and SOLID 45 elements are selected based on the characteristics of the elements to mesh the three-dimensional geomechanical model in one embodiment of the present application. The improvement of meshing accuracy is achieved by using a free meshing method, and the meshing accuracy of the top surface of the geomechanical model is higher than that of the bottom surface and the four side surfaces. After the mathematical model is constructed, the spatial coordinates (E_X, E_Y, E_Z) of all elements of the geomechanical model and the spatial coordinates (N_X, N_Y, N_Z) of the nodes are derived from the finite element analysis software. The total number of elements and the total number of nodes of each three-dimensional geomechanical model after meshing in step S21 are shown in FIG. Figure 4 .

[0125] In step S33, the static rock mechanical parameters (Young's modulus, Poisson's ratio and density) of the target stratum in the study area obtained from step S22 are assigned to each unit in the mathematical model, and different boundary conditions and constraints are set for the mathematical model, and numerical simulation of the tectonic stress field is carried out to obtain the simulation results corresponding to each test.

[0126] According to the calculation results of the static rock mechanical parameters in step S22, all elements of all geomechanical models are assigned static rock mechanical parameters (Young's modulus, Poisson's ratio and density). H / σ V ), the ratio of the regional minimum horizontal principal stress to the regional maximum horizontal principal stress (σ h / σ H ) These three parameters are defined as the key parameters affecting the numerical simulation results of the tectonic stress field, among which the influence of θ has been reflected in step S31, and σ H / σ V and σ h / σ HThe influence of needs to be realized by changing the boundary conditions of the mathematical model. According to the research accuracy requirements, set σ H / σ V and σ h / σ H The test data set, where σ h / σ H Less than 1. The boundary conditions and constraints of the mathematical model mainly include the regional stress application method, strength application and displacement constraint method on the side and top surfaces of the mathematical model. After setting reasonable boundary conditions and constraints for each mathematical model, a numerical simulation of the tectonic stress field is carried out to obtain the simulation results corresponding to each test. The focus is on statistically analyzing the horizontal maximum principal stress intensity at each position of the model in each numerical simulation of the tectonic stress field. For any simulation test i (where 1≤i≤n, n is the total number of tectonic stress field numerical simulation tests carried out), a one-to-one correspondence between the node position and the horizontal maximum principal stress can be obtained, that is, {σ H} i =F i (N_X,N_Y,N_Z).

[0127] In one embodiment of the present application, according to the calculation results of the static rock mechanical parameters in step S22, all units of all mathematical models are assigned the same static rock elastic parameters, specifically the static Young's modulus (39.15 GPa), the static Poisson's ratio (0.232) and the rock density (2.683 g / cm 3 ). When carrying out numerical simulation of tectonic stress field for each mathematical model, σ V Both are set to 100 MPa to study the deformation characteristics of the model under the same overlying stratum gravity. H / σ V is set to 1, 2, and 3, i.e., 100, 200, and 300 MPa, to investigate σ H / σ V Deformation characteristics of the model at different times. h / σ H is set to 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9 to study σ h / σ H Considering the comprehensiveness of this study and the efficiency of numerical simulation based on FEM, it is also added that when σ h / σ H When σ is equal to 0.3, H / σ V Numerical simulation of the three-dimensional tectonic stress field when σ is equal to 1.25, 1.50, 1.75, 2.25, 2.50 and 2.75 to more clearly show H / σ VTherefore, in one embodiment of the present application, a total of 297 three-dimensional tectonic stress field numerical simulation experiments are involved.

[0128] like Figure 4 As shown, the displacement of the nodes on the XZ surface on the south side of the 3D geomechanical model is not constrained and is imposed with σ H The displacement of the nodes on the east side of the YZ plane is not constrained, and the horizontal minimum principal stress (σ h The displacements of the nodes on the XZ surface on the north side and the YZ surface on the west side of the geomechanical model in the X and Y directions are set to a constant value of 0, and the displacement in the Z direction is not constrained. The displacements of the nodes on the XY surface on the upper side of the geomechanical model are not constrained, and a vertical stress (σ V ), the displacement in the Z direction of the nodes on the XY plane on the lower side of the model is set to a constant value of 0, and the displacements in the X and Y directions are not constrained. The entire three-dimensional geomechanical model is given a vertical downward gravitational acceleration to make the simulation results closer to the characteristics of the actual geological situation where vertical stress increases with increasing depth. Under the above boundary conditions and constraints, a numerical simulation of the tectonic stress field is carried out on each mathematical model to obtain the simulation results corresponding to each model, where the focus is on statistically analyzing the horizontal maximum principal stress intensity of the nodes on the top surface of the model in each numerical simulation of the tectonic stress field. In one embodiment of the present application, a plane contour map of the horizontal maximum principal stress intensity of the top interface of the model drawn according to the results of the numerical simulation of the tectonic stress field is partially shown in Figure 5 and Figure 6 For any simulation test i (where 1≤i≤297) in the embodiment of the present application, a one-to-one correspondence between the node position and the horizontal maximum principal stress can be obtained, that is, {σ H} i =F i (N_X,N_Y,N_Z).

[0129] Step S4: For each simulation test, the distances of all nodes from the vertical fracture are calculated, and the interval length or number of intervals is set according to the distance and the requirements of research accuracy. When the distance of the node from the vertical fracture is in a certain interval, the average value of the horizontal maximum principal stress at the corresponding node is counted.

[0130] Step S41: For each simulation test i, calculate the distance D of all nodes from the vertical fracture. i (N_X,N_Y,N_Z).

[0131] For the i-th simulation test, the coordinates (N_X, N_Y, N_Z) of each node in the corresponding mathematical model have been obtained in step S32. First, all nodes are rotated counterclockwise around the axis (N_X = 0, N_Y = 0) by an angle θ (the definition of θ is shown in step S31). The new coordinates of the nodes (N_X', N_Y', N_Z') are:

[0132]

[0133] The distance D between the node and the vertical fracture i The expression for (N_X, N_Y, N_Z) is:

[0134] D i (N_X, N_Y, N_Z) = N_X' (6)

[0135] After calculating the distance D1(N_X, N_Y, N_Z) of all nodes from the vertical fracture in the first simulation test by using equations and equation groups (5) and (6), this operation can be repeated until the distances of nodes from the vertical fracture in all simulation tests are calculated.

[0136] In one embodiment of the present application, the distances between nodes and vertical fractures in 297 simulation tests were calculated according to equations (5) and (6). The plane contour map of the distances between the model top interface nodes and the vertical fractures in some simulation tests is shown as follows: Figure 7 and Figure 8 shown.

[0137] Step S42: setting the interval length or number of intervals according to the distance between the node and the vertical fault obtained in step S41 and the research accuracy requirement.

[0138] According to the distance between the node and the vertical fracture obtained in step S41, the interval length or number of intervals is set according to the research accuracy. If the interval length is set, for any simulation test, the minimum and maximum values of the distance between the node and the vertical fracture are first counted and set as D MIN and D MAX , then according to D MAX and D MIN The difference and research accuracy require, set the interval length D INT , the higher the research accuracy requirement, D INT The smaller the D INT If the number of intervals is set, a uniform number of intervals n can be set for all simulation tests, and then the minimum value D of the distance between the node and the vertical fracture in each simulation test can be used. MIN and the maximum value D MAX Divide into n intervals, where the distance D between the node of the i-th interval (1≤i≤n and a positive integer) and the vertical fracture is iThe value range of (N_X, N_Y, N_Z) is [D MIN +(i-1) / n*(D MAX -D MIN ),D MIN +i / n*(D MAX -D MIN )].

[0139] In one embodiment of the present application, the number of intervals is set to 30 according to the distance between the node and the vertical fracture obtained in step S41. Then, for each simulation test, the distance D between the node and the vertical fracture is i (N_X, N_Y, N_Z) are divided into 30 corresponding numerical intervals, and the distance D between the node in the i-th interval (1≤i≤30 and a positive integer) and the vertical fracture is i The value range of (N_X, N_Y, N_Z) is [D MIN +(i-1) / 30*(D MAX -D MIN ),D MIN +i / 30*(D MAX -D MIN )].

[0140] Step S43: According to the interval length or number of intervals set in step S42, for each simulation test, the average value of the horizontal maximum principal stress at the corresponding node when the distance from the node to the vertical fracture is in a certain interval is calculated.

[0141] According to the interval length or number of intervals set in step S42, for the kth (1≤k≤n and a positive integer, n is the total number of simulation tests) simulation test, if the distance D between a node and the vertical fracture is k (N_X, N_Y, N_Z) is in the interval [D MIN +(i-1) / m*(D MAX -D MIN ),D MIN +i / m*(D MAX -D MIN )], where i is the interval number, m is the total number of intervals, 1≤i≤m and is a positive integer, then the horizontal maximum principal stress simulation result corresponding to the node is marked in the corresponding interval. For the kth simulation test, after counting the distance D of all nodes from the vertical fracture, k After (N_X, N_Y, N_Z) belongs to which interval, the corresponding horizontal maximum principal stress simulation results are also marked in each interval. For the i-th interval, calculate the average value S of all horizontal maximum principal stress simulation results marked as belonging to this interval H(k,i). Repeat this operation until the average value of the horizontal maximum principal stress simulation results of all intervals in all simulation tests is calculated. In subsequent analysis, the interval number is not used because it does not directly reflect the position of the corresponding interval. Instead, the median value of the interval is used, that is, D MIN +(i-0.5) / m*(D MAX -D MIN ), which directly represents the overall distance between the nodes in the interval and the vertical fracture.

[0142] In one embodiment of the present application, the number of intervals set in step S42 is 30. For the kth (1≤k≤297 and a positive integer) simulation test, if the distance D between a node and the vertical fracture is k (N_X, N_Y, N_Z) is in the interval [D MIN +(i-1) / 30*(D MAX -D MIN ),D MIN +i / 30*(D MAX -D MIN )], where i is the interval number, the horizontal maximum principal stress simulation result corresponding to the node is marked in the corresponding interval. For the kth simulation test, after counting the distance D from all nodes to the vertical fracture, k After (N_X, N_Y, N_Z) belongs to which interval, the corresponding horizontal maximum principal stress simulation results are also marked in each interval. For the i-th interval, calculate the average value S of all horizontal maximum principal stress simulation results marked as belonging to this interval H (k, i). Repeat this operation until the average value of the horizontal maximum principal stress simulation results of all intervals in all simulation tests is calculated. H (k,i) and the median D of the corresponding interval i MIN +(i-0.5) / 30*(D MAX -D MIN ) Figure 9 and Figure 10 shown.

[0143] Step S5: For each simulation test, select one "node-to-vertical-fault distance" interval closest to the vertical fault on each side of the vertical fault, and determine the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault.

[0144] In step S51 , one “distance between a node and a vertical fracture” interval closest to the vertical fracture is selected on each side of the vertical fracture.

[0145] For each simulation test, the corresponding “distance of the node from the vertical fracture” intervals have been obtained from step S43 and the median values D of these intervals have been calculated. MIN +(i-0.5) / m*(D MAX -D MIN ), where i is the interval number, m is the total number of intervals, 1≤i≤m and is a positive integer. Then for all interval medians less than or equal to 0, there is an interval median with the smallest absolute value. This interval is closest to the vertical fracture, and the median of this interval is defined as D LEFT-MIN The average value of the maximum horizontal principal stress of the nodes in this interval is defined as S LEFT-MIN For all interval medians greater than or equal to 0, there is also an interval median with the smallest absolute value. This interval is closest to the vertical fracture. The median of this interval is defined as D RIGHT-MIN The average value of the maximum horizontal principal stress of the nodes in this interval is defined as S RIGHT-MIN When D LEFT-MIN Equal to D RIGHT-MIN And when it is equal to 0, the two intervals mentioned above are actually the same interval. In the subsequent analysis, the two intervals are still analyzed using the same interval. The above steps are applied to all simulation tests until D LEFT-MIN (k), S LEFT-MIN (k), D RIGHT-MIN (k) and S RIGHT-MIN (k), where 1≤k≤n is a positive integer, and n is the total number of simulation trials.

[0146] In one embodiment of the present application, for each simulation test, the corresponding “distance between the node and the vertical fracture” intervals have been obtained from step S43 and the median values D of these intervals have been calculated. MIN +(i-0.5) / 30*(D MAX -D MIN ), where i is the interval number, 1≤i≤30 and is a positive integer. Then for all interval medians less than or equal to 0, there is an interval median with the smallest absolute value. This interval is closest to the vertical fracture. The median of this interval is defined as D LEFT-MIN The average value of the maximum horizontal principal stress of the nodes in this interval is defined as S LEFT-MIN For all interval medians greater than or equal to 0, there is also an interval median with the smallest absolute value. This interval is closest to the vertical fracture. The median of this interval is defined as D RIGHT-MIN The average value of the maximum horizontal principal stress of the nodes in this interval is defined as S RIGHT-MIN The above steps are applied to all simulation tests until D LEFT-MIN (k), S LEFT-MIN (k), DRIGHT-MIN (k) and S RIGHT-MIN (k), where 1≤k≤297 and is a positive integer.

[0147] Step S52 , respectively determining the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault.

[0148] For each simulation test, D has been obtained from step S51. LEFT-MIN (k), S LEFT-MIN (k), D RIGHT-MIN (k) and S RIGHT-MIN (k), where 1≤k≤n is a positive integer, and n is the total number of simulation tests. Now we need to find a position on both sides of the vertical fault that can represent the maximum horizontal principal stress mutation, which will serve as the boundary of the tectonic stress disturbance zone on both sides of the vertical fault. H (k,i), where i is the interval number, m is the total number of intervals, 1≤i≤m and is a positive integer. First, find the position representing the maximum horizontal principal stress mutation when the median of the interval is less than or equal to 0, specifically:

[0149] When S H When (k,i) satisfies the following inequality, the i-th interval is identified as the interval representing the location of the horizontal maximum principal stress mutation:

[0150] [S H (k,i)-S H (k,i-1)]·[S H (k,i+1)-S H (k,i)]<0(7)

[0151] If there is no S that satisfies inequality (9) H (k,i), if S H (k,i) satisfies the following inequality group, and the i-th interval can also be considered to represent the location of the horizontal maximum principal stress mutation:

[0152]

[0153] Among them S * The threshold of the maximum horizontal principal stress is exceeded. When this threshold is exceeded, it indicates that the maximum horizontal principal stress between adjacent intervals has changed significantly. The threshold needs to be determined according to the research accuracy and is usually set to a smaller value, such as 1MPa. The median of the interval representing the sudden change of the maximum horizontal principal stress when the median of the interval is less than or equal to 0 is defined as D LEFT-MAX (k), the average value of the horizontal maximum principal stress corresponding to this interval is defined as S LEFT-MAX (k).

[0154] When the median values of the interval are all greater than or equal to 0, find the position representing the maximum horizontal principal stress mutation, specifically:

[0155] When S H When (k,i) satisfies inequality (7), the i-th interval is identified as the interval representing the location of the horizontal maximum principal stress mutation. If there is no S that satisfies inequality (7) H (k,i), if S H (k,i) satisfies the following inequality group, and the i-th interval can also be considered to represent the location of the horizontal maximum principal stress mutation:

[0156]

[0157] The median of the interval representing the maximum principal stress mutation when the median of the interval is greater than or equal to 0 is defined as D RIGHT-MAX (k), the average value of the horizontal maximum principal stress corresponding to this interval is defined as S RIGHT-MAX (k).

[0158] By repeating the above steps, the intervals representing the maximum horizontal principal stress mutation on both sides of the vertical fault in all simulation tests can be obtained. The width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault can be calculated according to the following set of equations:

[0159]

[0160] Among them, when the median values of the interval are less than or equal to 0, the width and intensity of the tectonic stress disturbance zone are △D LEFT (k) and △S LEFT (k), when the median values of the interval are greater than or equal to 0, the width and intensity of the tectonic stress disturbance zone are △D RIGHT (k) and △S RIGHT (k) Since the length of vertical faults in the model setting is 5000m, △D LEFT (k) and ΔD RIGHT (k) are divided by 5000m to obtain the tectonic stress disturbance zone width coefficient △D' LEFT (k) and △D' RIGHT (k).

[0161] In one embodiment of the present application, the horizontal maximum principal stress threshold S * is set to 1 MPa, and the △D' obtained according to the above steps LEFT (k), △S LEFT (k), △D' RIGHT (k) and △S RIGHT (k) As shown in Table 4.

[0162] Table 4

[0163]

[0164]

[0165]

[0166]

[0167]

[0168]

[0169]

[0170] Step S6, draw the contour map of the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault with respect to the boundary conditions of the numerical simulation of the tectonic stress field, and judge the width and intensity of the stress disturbance zone of the vertical fault that breaks through the target layer in the study area based on the contour map.

[0171] Step S61: Draw ΔD' when θ is equal to 10°, 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80° in sequence. LEFT , △S LEFT ,△D' RIGHT and △S RIGHT The plane contour map of the graph is σ H / σ V , the vertical axis is σ h / σ H .

[0172] When θ is equal to 10°, ΔD' is plotted according to the result in step S52. LEFT , △S LEFT ,△D' RIGHT and △S RIGHT About σ H / σ V and σ h / σ H Repeat this operation when θ is equal to 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80° in sequence, and 36 contour maps can be obtained.

[0173] In one embodiment of the present application, when θ is equal to 10°, ΔD′ is plotted according to Table 4 in step S52. LEFT , △S LEFT ,△D' RIGHT and △S RIGHT About σ H / σ V and σ h / σH Repeat this operation when θ is equal to 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80° in sequence to obtain 36 contour maps.

[0174] Step S62, judging the width and intensity of the stress disturbance zone of the vertical fault that cuts through the target stratum in the study area based on a series of contour maps obtained from step S61.

[0175] According to steps S11 and S12, the strike direction and σ of the vertical faults Fi (where 1≤i≤n, n is the total number of vertical faults) that penetrate the target strata in the study area can be obtained respectively. H The azimuth of the vertical fracture Fi and σ can be calculated. H The angle θ i , and we know which side of the vertical fault the regional stress source is located on (see the definition of the direction of the regional horizontal maximum principal stress in step S12). According to the analysis of steps S3 to S5, when the regional stress source is located on one side of the vertical fault, the width and intensity of the tectonic stress disturbance zone on this side are respectively defined as △D RIGHT and △S RIGHT , and the width and intensity of the tectonic stress disturbance zone on the side of the vertical fault away from the regional stress source are defined as △D LEFT and △S LEFT Therefore, the width and intensity of the tectonic stress disturbance zone on each vertical fault are determined by the location of the regional stress source on that side of the vertical fault. RIGHT and △S RIGHT , and the other side is △D LEFT and △S LEFT 10°, 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80° are numbered as α1 to α9. For vertical fracture Fi, if α j ≤θ i <α j+1 , where 1≤j≤8, then according to and α j and α j+1 Corresponding △D' LEFT , △S LEFT ,△D' RIGHT and △S RIGHT The contour map (obtained from step S61) obtains the △D' of the vertical fracture Fi LEFT , △S LEFT ,△D' RIGHT and △S RIGHT , the specific steps are as follows:

[0176] The first step is to obtain the regional maximum principal stress intensity, regional minimum principal stress intensity and vertical stress intensity during the tectonic movement period related to the sliding along the vertical fault Fi from step S12, and calculate the σ corresponding to the vertical fault Fi based on this. H / σ V and σ h / σ H . Obtained from step S61 and α j Corresponding △D' LEFT , △S LEFT ,△D' RIGHT and △S RIGHT Contour map, according to the σ corresponding to the vertical fracture Fi H / σ V and σ h / σ H From and α j Corresponding △D' LEFT , △S LEFT ,△D' RIGHT and △S RIGHT Obtain the △D' corresponding to the vertical fracture Fi in the contour map LEFT (Fi,α j ), △S LEFT (Fi,α j ),△D' RIGHT (Fi,α j ) and △S RIGHT (Fi,α j );

[0177] The second step is to obtain the sum α from step S61. j+1 Corresponding △D' LEFT , △S LEFT ,△D' RIGHT and △S RIGHT Contour map, according to the σ corresponding to the vertical fracture Fi H / σ V and σ h / σ H From and α j+1 Corresponding △D' LEFT , △S LEFT ,△D' RIGHT and △S RIGHT Obtain the △D' corresponding to the vertical fracture Fi in the contour map LEFT (Fi,α j+1 ), △S LEFT (Fi,α j+1 ),△D' RIGHT (Fi,α j+1 ) and △S RIGHT (Fi,α j+1 );

[0178] The third step is to calculate the actual △D' corresponding to the vertical fracture Fi according to the following equation group (11): LEFT (Fi,θ i ), △S LEFT (Fi,θ i ),△D' RIGHT (Fi,θ i ) and △S RIGHT (Fi,θ i ):

[0179]

[0180] Repeat the above three steps until the actual △D' of all vertical fractures is calculated. LEFT (Fi,θ i ), △S LEFT (Fi,θ i ),△D' RIGHT (Fi,θ i ) and △S RIGHT (Fi,θ i ). It should be noted that if θ i ≤α1, then set j in the above three steps to 1 and j+1 to 2; if θ i ≥α9, then set j in the above three steps to 8 and j+1 to 9. LEFT (Fi,θ i ) and △D' RIGHT (Fi,θ i ) It is also necessary to calculate the actual width of the tectonic stress disturbance zone based on the actual length of the fault, because △D' LEFT (Fi,θ i ) and △D' RIGHT (Fi,θ i ) is a ratio rather than actual width, see step S52 for details.

[0181] In one embodiment of the present application, the strike and σ of the vertical fault Fi (where 1≤i≤13) that cuts through the target strata in the study area can be obtained according to Table 1 and S12, respectively. H The azimuth of the vertical fracture Fi and σ H The angle θ i The angles 10°, 20°, 30°, 40°, 45°, 50°, 60°, 70°, and 80° are numbered α1 to α9. The following uses the vertical fracture F1 listed in Table 1 as an example to illustrate how to calculate △D. LEFT (Fi,θ i ), △S LEFT (Fi,θ i),△D RIGHT (Fi,θ i ) and △S RIGHT (Fi,θ i ).

[0182] For the vertical fault F1, its strike is 59.2° (Table 1), and its slip along the strike occurred during the tectonic activity period I (Table 2). H The azimuth angle is 0°, σ H / σ V and σ h / σ H are equal to 2.4 and 0.33 respectively, Fi and σ H The angle θ i The angle is 59.2°, and the regional stress source is located on the southeast side of the vertical fault F1. Therefore, the width and intensity of the tectonic stress disturbance zone calculated on the southeast side are △D' RIGHT and △S RIGHT The width and intensity of the tectonic stress disturbance zone calculated on the northwest side are △D' LEFT and △S LEFT Angle θ i Between α6 and α7, so according to the △D' corresponding to α6 and α7 obtained from step S61 LEFT , △S LEFT ,△D' RIGHT and △S RIGHT Analyze the contour map and the △D' corresponding to α6 LEFT , △S LEFT ,△D' RIGHT and △S RIGHT The contour maps are as follows Figures 11 to 14 As shown, △D' corresponding to α7 LEFT , △S LEFT ,△D' RIGHT and △S RIGHT The contour maps are as follows Figures 15 to 18 The σ related to the strike slip of the vertical fault F1 is shown in Figure 2. H / σ V and σ h / σ H are equal to 2.4 and 0.33 respectively, so according to Figures 11 to 14 It can be obtained when σ H / σ V and σ h / σ H When △D' is equal to 2.4 and 0.33 respectively LEFT , △S LEFT ,△D' RIGHT and △S RIGHTare equal to 0.1400, 64.88MPa, 0.1363 and 75.83MPa respectively. Figures 15 to 18 It can be obtained when σ H / σ V and σ h / σ H When △D' is equal to 2.4 and 0.33 respectively LEFT , △S LEFT ,△D' RIGHT and △S RIGHT These values are 0.0820, 53.70 MPa, 0.1368, and 56.00 MPa, respectively. Substituting these values into formula (12) yields:

[0183]

[0184] From this, we can calculate the △D' of the vertical fracture F1 LEFT , △S LEFT ,△D' RIGHT and △S RIGHT The values are 0.0867, 54.59 MPa, 0.1368, and 57.59 MPa, respectively. Based on the 2036 m extension length of vertical fault F1 (see Table 1), it can be calculated that the tectonic stress disturbance zone on the southeast side of vertical fault F1 is 278.52 m wide and 57.59 MPa strong, while the tectonic stress disturbance zone on the northwest side of vertical fault F1 is 176.52 m wide and 54.59 MPa strong.

[0185] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

Claims

1. A method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault, characterized in that: include: Step S1: determine the study area and target strata, calculate the strike and plane extension length of vertical faults in the study area, determine the key period of fault activity and the paleo-tectonic stress field environment; further determine the direction and intensity of the regional maximum horizontal principal stress, the intensity of the minimum horizontal principal stress, and the vertical stress intensity during the key period; Step S2, solving static rock mechanical parameters based on core, conventional logging and array acoustic logging data of the target layer; Step S3, constructing a geomechanical model and a mathematical model related to the sliding of the vertical fault along the strike, and conducting a numerical simulation test of the tectonic stress field to obtain the intensity of the maximum horizontal principal stress at each location; Step S4, in the numerical simulation experiment, calculate the distances from all nodes to the vertical fracture, divide the intervals according to the distance and research accuracy requirements, set the interval length or the number of intervals, and calculate the average value of the horizontal maximum principal stress in each interval; use the median value of the interval to represent the overall distance of the nodes in the interval from the vertical fracture; Step S5, selecting the nearest intervals on both sides of the vertical fault, and determining the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault; Step S6, drawing a contour map of the tectonic stress disturbance zone width and disturbance intensity as a function of boundary conditions, and determining the stress disturbance zone width and disturbance intensity based on the contour map; Wherein, the step S3 specifically includes: Based on the geological characteristics of the study area, a three-dimensional geomechanical model was constructed. The model included the location, strike, and extension of the vertical fault, as well as the geometric characteristics of the target strata. The vertical fault was set as a discontinuity with a certain friction coefficient to simulate the relative sliding of the strata on both sides of the fault. The geomechanical model is converted into a mathematical model, and meshing is performed using finite element analysis software to divide the model into several units and nodes; Assign the static rock mechanics parameters obtained from step S2 to each unit in the mathematical model; set the boundary conditions and constraints of the mathematical model according to the research accuracy requirements; conduct numerical simulation tests of the tectonic stress field, and obtain the horizontal maximum principal stress intensity at each location of the model from the test results; Compare the simulation results with actual geological data to verify the accuracy and reliability of the model.

2. The method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault according to claim 1, characterized in that: The step S2 comprises: Calculating dynamic rock mechanics parameters based on conventional logging and array acoustic logging data, wherein the dynamic rock mechanics parameters include dynamic Young's modulus and dynamic Poisson's ratio; Obtaining static rock mechanics parameters from rock mechanics experiments on core samples, wherein the static rock mechanics parameters include static Young's modulus and static Poisson's ratio parameters; The dynamic rock mechanics parameters and static rock mechanics parameters of the target strata in the study area are fitted to obtain the conversion model of dynamic and static rock mechanics parameters: Where, E represents the static Young's modulus; μ represents the static Poisson's ratio; E d represents the dynamic Young's modulus; μ d represents the dynamic Poisson's ratio.

3. The method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault according to claim 2, characterized in that: The calculation of dynamic rock mechanical parameters based on conventional logging and array acoustic logging data includes: Based on conventional logging data, a conversion model between acoustic time difference AC and compressional time difference DTC is constructed; based on array acoustic logging data, a conversion model between compressional time difference DTC and shear time difference DTS is constructed; DTC and DTS are calculated using the conversion model. The calculation formulas for the dynamic Young's modulus and dynamic Poisson's ratio of the target formation during drilling are as follows: Among them, E d Represents dynamic Young's modulus, in GPa; μ d is the dynamic Poisson's ratio; Δt p is the longitudinal wave time difference DTC, in μs / ft; Δt s is the shear wave time difference DTS, with the unit of μs / ft; ρ is the density, with the unit of g / cm3.

4. The method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault according to claim 1, characterized in that: The step S4 comprises: Adjust the coordinate system of the model to make it consistent with the direction of the regional stress field, calculate the distance from each node in the model to the vertical fracture; count the maximum value D of the distance from the node to the vertical fracture MAX and minimum value D MIN ; Set the interval length D according to the research accuracy requirements INT ; The range of each interval is: [D MIN +(i-1) / n*(D MAX -D MIN ),D MIN +i / n*(D MAX -D MIN )] Where n represents the number of preset intervals; i represents the interval number; For each interval, record the maximum horizontal principal stress of each node in the interval, and after multiple numerical simulation tests of the stress field, calculate the average value of the maximum horizontal principal stress in the interval; The median of the interval is used to represent the overall distance between the nodes in the interval and the vertical fracture; the median of the interval is: D MIN +(i-0.5) / m*(D MAX -D MIN )。 5. The method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault according to claim 1, characterized in that: The step S5 comprises: Get the median of each interval from step S4. For all the interval medians less than or equal to 0, find the interval median with the smallest absolute value, which is recorded as D LEFT-MIN The corresponding average value of the maximum horizontal principal stress is recorded as S LEFT-MIN ; For all interval medians greater than or equal to 0, find the interval median with the smallest absolute value, recorded as D RIGHT-MIN The corresponding average value of the maximum horizontal stress is recorded as S RIGHT-MIN If D LEFT-MIN =D RIGHT-MIN =0, then the two intervals are actually the same interval, and the same interval is still used in subsequent analysis; For each interval, check the average value of the horizontal maximum principal stress S H (k,i) Whether one of the following conditions is met: [S H (k,i)-S H (k,i-1)]·[S H (k,i+1)-S H (k,i)]<0 Among them, S * Indicates the horizontal maximum principal stress threshold; S H represents the average value of the horizontal maximum principal stress simulation results; i represents the interval number; k represents the number of the simulation experiment; For the median of the interval less than or equal to 0, find the interval that meets the above conditions and record the median of the interval as D LEFT-MAX (k), the corresponding average value of the maximum horizontal principal stress is recorded as S LEFT-MAX (k); For the median of the interval greater than or equal to 0, find the interval that meets the above conditions and record the median of the interval as D RIGHT-MAX (k), the corresponding stress average value is recorded as S RIGHT-MAX (k); The intervals representing the maximum horizontal principal stress mutation on both sides of the vertical fault in all simulation tests were obtained, and the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault were calculated according to the following set of equations: When the median values of the interval are less than or equal to 0, the width and intensity of the tectonic stress disturbance zone are △D LEFT (k) and △S LEFT (k); When the median values of the interval are greater than or equal to 0, the width and intensity of the tectonic stress disturbance zone are △D RIGHT (k) and △S RIGHT (k); △D LEFT (k) and ΔD RIGHT (k) are divided by the length of the vertical fault preset in the model to obtain the tectonic stress disturbance zone width coefficient △D' LEFT (k) and △D' RIGHT (k).

6. The method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault according to claim 1, characterized in that: The step S6 comprises: Plotting θ at different angles, △D' LEFT (k), △S LEFT (k), △D' RIGHT (k) and △S RIGHT (k) plane contour map, the horizontal axis of the plane contour map is the ratio of the horizontal maximum principal stress to the vertical stress σ H / σ V , the vertical axis is the ratio of the horizontal minimum principal stress to the horizontal maximum principal stress σ h / σ H ; where θ represents the angle between the strike of the vertical fault and the direction of the regional maximum horizontal principal stress; According to the location of the regional stress source, the width and intensity of the tectonic stress disturbance zone on this side are determined as △D RIGHT (k) and △S RIGHT (k); the other side is △D LEFT (k) and △S LEFT (k); Observe the trends in the contour map to determine the changing patterns of the width and intensity of the stress disturbance zone and identify the boundary of the stress disturbance zone, that is, the location of the stress mutation. Based on the contour map, determine the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault. Calculate the disturbance band width coefficient △D' LEFT (k) and △D' RIGHT (k), and the disturbance intensity △S LEFT (k) and △S RIGHT (k).

7. The method for determining the width and intensity of the tectonic stress disturbance zone on both sides of a vertical fault according to claim 6, characterized in that: Determining the width and intensity of the tectonic stress disturbance zone on both sides of the vertical fault includes: The different θ values are numbered as α1 to α t , where θ represents the angle between the strike of the vertical fault and the direction of the regional horizontal maximum principal stress; Get and α j The corresponding contour map, and obtain the △D' corresponding to the vertical fracture Fi from the contour map LEFT (Fi,α j ), △S LEFT (Fi,α j ),△D' RIGHT (Fi,α j ) and △S RIGHT (Fi,α j ); Get and α j+1 The corresponding contour map, and obtain the △D' corresponding to the vertical fracture Fi from the contour map LEFT (Fi,α j+1 ), △S LEFT (Fi,α j+1 ),△D' RIGHT (Fi,α j+1 ) and △S RIGHT (Fi,α j+1 ); Calculate the actual △D' corresponding to the vertical fracture Fi LEFT (Fi,θ i ), △S LEFT (Fi,θ i ),△D' RIGHT (Fi,θ i ) and △S RIGHT (Fi,θ i ), the calculation formula is: After calculating the actual △D of all fractures LEFT (Fi,θ i ) and △D RIGHT (Fi,θ i ), the actual width of the tectonic stress disturbance zone needs to be calculated based on the actual length of the fault. LEFT (Fi,θ i ) and △D' RIGHT (Fi,θ i ) is the disturbance band width coefficient.

Images