A horizontal displacement curve method for vertical fractures

Abstract

This application discloses a method for calculating the horizontal displacement curve of a vertical fracture, including: constructing a geomechanical model of the vertical fracture and calculating a corresponding mathematical model; setting different boundary conditions and constraints on the mathematical model and conducting numerical simulation of the tectonic stress field to obtain simulation results corresponding to each test; drawing the corresponding vertical displacement curve and horizontal displacement curve of the vertical fracture based on the simulation results corresponding to the different boundary conditions and constraints, and constructing a one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH, thereby calculating the linear fitting equations of the major axis and minor axis; calculating the major axis and minor axis of the axis-aligned elliptical implicit function of the target vertical fracture based on the DV curve and major axis of the target vertical fracture, and calculating the DH curve of the target vertical fracture. Through this method, the accuracy of the calculation of the DH curve of the target vertical fracture can be improved.

Description

Technical Field

[0001] The present application relates to the technical field of vertical fractures, and in particular to a method for calculating a horizontal displacement curve of a vertical fracture. Background Art

[0002] The structural characteristics of underground vertical faults have a very important impact on the formation of reservoirs and oil and gas accumulation in underground basins. At present, the research on underground vertical faults mainly focuses on the characterization of geometric characteristics. There are few literatures on systematic research, especially quantitative research on kinematic characteristics. Unlike normal faults or reverse faults that move along the dip and whose sliding direction can be determined based on the configuration relationship of the strata on both sides of the fault, vertical faults are characterized by more significant sliding along the strike. The determination of their strike-slip distance is one of the main technical difficulties at present. There are currently three methods for determining the sliding distance (horizontal displacement) along the strike of underground vertical faults:

[0003] Method 1: Determine the fault plane and the sliding direction and distance of the strike-slip fault based on the offset marks of the wavelet characteristics of the horizontal time slice. When the overlying rock layer is slightly tilted due to the foreset layer affected by early tilt and strike-slip deformation on the horizontal time slice, the sliding direction and distance of the strike-slip fault can be determined by restoring the original geometric form of the wavelet characteristics. However, the penetration points of the fault body involving linear features (such as stratigraphic marks or paleo-river channels) such as the offset marks of the wavelet characteristics of the horizontal time slice are difficult to obtain, especially in complex deformation areas. This makes the determination of the sliding direction and distance of the underground strike-slip fault have certain errors and is often speculative.

[0004] Method 2: Estimating strike-slip distance based on the offset of characteristic points on the stratigraphic thickness curve for the same set of strata on either side of a vertical fault. In this method, the offset of characteristic points on the stratigraphic thickness curve is assumed to be caused by along-strike slip of the vertical fault. Therefore, calculating the distance between them can restore the initial state of the vertical fault before along-strike slip occurred, and thus calculate the horizontal displacement of the vertical fault. However, in areas with complex tectonic activity and sedimentary system evolution patterns, the stratigraphic thickness on either side of the vertical fault is influenced by many factors and sometimes cannot fully reflect along-strike slip of the vertical fault. In addition, characteristic points on the stratigraphic thickness curve are often difficult to identify in areas of intense deformation. When multiple suspected characteristic points exist, the matching of the characteristic points on either side of the vertical fault is often poor. For example, if there are four characteristic points on one side of the vertical fault and six on the other side, it is difficult to clearly match these points, and the artifacts are significant.

[0005] Method 3: Calculate the horizontal displacement of the vertical fault based on the horizontal offset distance between pre-existing geological boundaries on either side of the vertical fault. Pre-existing geological boundaries have a very broad meaning, including pre-existing sedimentary facies, paleo-faults, paleo-drainage systems, etc. These geological boundaries undergo relative displacement due to the movement along the strike of the vertical fault. Therefore, calculating the relative displacement distance of these geological boundaries along the strike of the vertical fault can be used to restore the initial state of the vertical fault before the occurrence of along-strike slip, and then calculate the horizontal displacement of the vertical fault.

[0006] Furthermore, the three aforementioned methods only derive the horizontal displacement of a vertical fault at a specific location, typically near the center of the fault. However, a robust method for mapping the distribution of horizontal displacement along the fault's direction has yet to be developed. Obtaining horizontal displacement at different locations along the fault's direction remains a challenge. In other words, current research on horizontal displacement of vertical faults typically focuses on a one-dimensional "point" perspective, rather than a two-dimensional "line" perspective. Summary of the Invention

[0007] Based on this, in order to solve the above technical problems, a horizontal displacement curve method of vertical faults is provided to solve the problems of inaccurate calculation of horizontal displacement of vertical faults in the existing technology and inability to calculate the distribution characteristics of horizontal displacement along the direction of the vertical faults.

[0008] In a first aspect, a method for calculating a horizontal displacement curve of a vertical fracture is provided, the method comprising:

[0009] The static rock mechanical parameters of the target strata in the study area are calculated based on the core, conventional logging and array acoustic logging data of the target strata in the study area;

[0010] A geomechanical model is constructed based on the thickness of the target stratum in the study area, static rock mechanical parameters, and the planar extension length of the vertical fault, and the angle θ between the geomechanical model and any side surface of the vertical fault is changed to output multiple different geomechanical models; the geomechanical model is unitized using finite element analysis computer software to obtain a mathematical model; different boundary conditions and constraints are set for the mathematical model and a numerical simulation of the tectonic stress field is performed to obtain simulation results corresponding to each experiment; the simulation results include: the vertical displacement of the top interface of the mathematical model and the displacement in the same direction as the vertical fault; 0°<θ<90°;

[0011] According to the simulation results corresponding to different boundary conditions and constraint conditions, the corresponding vertical displacement DV curve and horizontal displacement DH curve of the vertical fracture are drawn, and a one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH is established; and the linear fitting equations of the major axis and the minor axis are calculated according to the one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH corresponding to different boundary conditions and constraint conditions;

[0012] Calculate the vertical displacement DV curve of the target vertical fracture;

[0013] The major axis and minor axis of the axis-aligned elliptical implicit function of the target vertical fracture are calculated according to the DV curve of the target vertical fracture and the linear fitting equations of the major axis and minor axis to determine the expression of the axis-aligned elliptical implicit function and calculate the DH curve of the target vertical fracture.

[0014] In the above solution, optionally, constructing a geomechanical model based on the thickness of the target stratum in the study area, static rock mechanics parameters, and the plane extension length of the vertical fracture includes:

[0015] The size of the regular quadrangular prism geomechanical model is set according to the plane extension length of the vertical fault and the thickness of the target stratum in the study area. The cross section of the target vertical fault in the geomechanical model close to the regional stress direction is set as the target surface, and the cross section away from the regional stress direction is set as the contact surface.

[0016] In the above scheme, further optionally, the use of finite element analysis computer software to unitize the geomechanical model to obtain a mathematical model includes: meshing the geomechanical model with tetrahedral units to obtain a corresponding mathematical model, and the mesh division density of the top interface of the mathematical model is higher than that of other positions.

[0017] In the above solution, optionally, the different boundary conditions and constraint conditions include: regional stress application methods, strength application and displacement constraint methods on the side and top surfaces of the mathematical model.

[0018] In the above scheme, further optionally, the stress application intensity is mainly changed by the ratio of the regional horizontal maximum principal stress to the vertical stress (σH / σV) and the ratio of the regional horizontal minimum principal stress to the regional horizontal maximum principal stress (σh / σH).

[0019] In the above solution, optionally, the linear fitting equations for calculating the major axis and the minor axis according to the one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH corresponding to different boundary conditions and constraint conditions include:

[0020] According to the one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH corresponding to each different boundary condition and constraint condition, find the vertical displacement DV(DH-MAX) corresponding to the maximum value DH-MAX of the horizontal displacement DH, and calculate the deflection angle α and the major axis of the axis-aligned ellipse;

[0021] The one-to-one correspondence between each set of vertical displacement DV and horizontal displacement DH is rotated counterclockwise by α around the coordinate origin to obtain (DH', DV'), and the axis-aligned ellipse implicit function fitting is performed on (DH', DV'), and the fitting function expression is calculated to obtain the minor axis of the axis-aligned ellipse;

[0022] According to the major axis and minor axis of the axis-aligned ellipse corresponding to each boundary condition and constraint condition, the linear fitting equations of the major axis and minor axis are solved.

[0023] In the above solution, optionally, calculating the major axis and minor axis of the axis-aligned elliptical implicit function of the target upright fracture according to the DV curve of the target upright fracture and the linear fitting equations of the major axis and minor axis includes:

[0024] Calculating the maximum horizontal displacement DH-MAX of the target vertical fracture and the corresponding vertical displacement DV(DH-MAX) according to the DV curve of the target vertical fracture;

[0025] The deflection angle of the target vertical fracture and the major axis of the axis-aligned ellipse are solved according to the maximum horizontal displacement DH-MAX of the target vertical fracture and the corresponding vertical displacement DV(DH-MAX).

[0026] In the above solution, optionally, calculating the major axis and minor axis of the axis-aligned elliptical implicit function of the target upright fracture according to the DV curve of the target upright fracture and the linear fitting equations of the major axis and minor axis includes:

[0027] Calculating the maximum horizontal displacement DH-MAX of the target vertical fracture and the corresponding vertical displacement DV(DH-MAX) according to the DV curve of the target vertical fracture;

[0028] The deflection angle of the target vertical fracture and the major axis of the axis-aligned ellipse are solved according to the maximum horizontal displacement DH-MAX of the target vertical fracture and the corresponding vertical displacement DV(DH-MAX).

[0029] In the above solution, optionally, the calculating and obtaining the target vertical fracture DH curve includes:

[0030] According to the expression of the axis-aligned elliptical implicit function of the target vertical fracture, the elliptical implicit function corresponding to the initial data set (DH, DV) is restored after rotating the deflection angle α clockwise around the origin, and the DH curve of the target vertical fracture is solved accordingly.

[0031] In a second aspect, a computer device includes a memory and a processor, wherein the memory stores a computer program, and when the processor executes the computer program, the method steps for obtaining a horizontal displacement curve of a vertical fracture described in the first aspect are implemented.

[0032] In a third aspect, a computer program product includes a computer program / instruction, which, when executed by a processor, implements the steps of the method for obtaining a horizontal displacement curve of a vertical fracture as described in the first aspect.

[0033] This application has at least the following beneficial effects:

[0034] This application constructs a geomechanical model based on the planar extension length of the vertical fault, sets different boundary conditions and constraints, and conducts numerical simulation of the tectonic stress field. The vertical displacement of the top interface of the mathematical model under different conditions and the displacement in the same direction as the vertical fault are obtained. Based on this, a vertical displacement DV curve and a horizontal displacement DH curve are plotted, and a one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH is established. From this correspondence, the linear fitting equations for the major and minor axes are calculated. Based on the linear fitting equations for the major and minor axes, the DH curve of the target vertical fault is determined based on the vertical displacement DV curve of the actual target vertical fault. Thus, the application obtains horizontal displacements at different locations. Furthermore, the present invention proposes using numerical simulation of the tectonic stress field and contact pair settings to restore the deformation process of the vertical fault, thereby solving the linear fitting equations for the major and minor axes. Then, based on the linear fitting equations for the major and minor axes, the horizontal displacement DH curve is inverted based on the actual target vertical displacement DV curve, thereby improving the calculation accuracy of the horizontal displacement curve. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 A schematic diagram of a flow chart of a method for generating a horizontal displacement curve of a vertical fracture provided in one embodiment of the present application;

[0036] Figure 2 A linear conversion mathematical model of AC and DTC obtained from AC and DTC data of the target strata in the study area provided in one embodiment of the present application;

[0037] Figure 3 A linear conversion mathematical model of DTC and DTS obtained from DTC and DTS data of the target strata in the study area provided in one embodiment of the present application;

[0038] Figure 4 Nine different geomechanical model diagrams corresponding to the setting of geomechanical model size, vertical fracture, boundary condition, and the angle θ between the vertical fracture and one side of the geomechanical model to nine different values, and corresponding mathematical model diagrams, provided in one embodiment of the present application;

[0039] Figure 5 A diagram showing the vertical deformation amplitude and horizontal deformation amplitude of a node on the top interface of a model obtained from a numerical simulation result of a tectonic stress field provided in one embodiment of the present application;

[0040] Figure 6 A diagram showing the vertical deformation amplitude and horizontal deformation amplitude of a node on the top interface of a model obtained from a numerical simulation result of a tectonic stress field provided in one embodiment of the present application;

[0041] Figure 7 A graph showing the vertical displacement (DV) of the top of the vertical fault F1 obtained from the numerical simulation results of the tectonic stress field provided in one embodiment of the present application;

[0042] Figure 8 A graph showing the horizontal displacement (DH) of the top of a vertical fault F1 obtained from numerical simulation results of a tectonic stress field provided in one embodiment of the present application;

[0043] Figure 9 A schematic diagram of obtaining the deflection angle α provided in one embodiment of the present application;

[0044] Figure 10 This is an example diagram of the results of performing axis-aligned elliptical implicit function fitting based on the first (DH', DV') data set provided in one embodiment of the present application;

[0045] Figure 11 This is an example diagram of the results of performing axis-aligned ellipse implicit function fitting based on the second (DH', DV') data set provided in one embodiment of the present application;

[0046] Figure 12 This is an example diagram of the results of performing axis-aligned ellipse implicit function fitting based on the third (DH', DV') data set provided in one embodiment of the present application;

[0047] Figure 13 A DV curve diagram of a target vertical fracture provided in one embodiment of the present application;

[0048] Figure 14 A schematic diagram of the maximum horizontal displacement of a target fracture and the position of the maximum horizontal displacement provided in one embodiment of the present application;

[0049] Figure 15 This is a diagram showing the calculation results of the DH curve of fracture F1 provided in one embodiment of the present application. DETAILED DESCRIPTION

[0050] 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.

[0051] In one embodiment, Figure 1As shown, a method for horizontal displacement curve of vertical fracture is provided, the method comprising:

[0052] Step S1: Calculate the static rock mechanical parameters of the target strata in the study area based on the core, conventional logging and array acoustic logging data of the target strata in the study area.

[0053] Specifically, step S1 comprises the following steps:

[0054] Step S11: According to the actual research requirements of oil and gas exploration and development, determine the study area, target strata and vertical faults as the research objects.

[0055] Select the study area and target strata based on the actual needs of oil and gas exploration and development research. Clarify the lithology of the target strata and the depth distribution characteristics of the top and bottom interfaces within the study area. Identify the vertical faults for quantitative vertical displacement research and quantitative horizontal displacement prediction. These vertical faults must completely penetrate the target strata, and the intersection lines with the top and bottom interfaces of the target strata must be continuous and clearly visible in the plane.

[0056] For example, a vertical fault F1 in Basin A, China, with a horizontal extension of approximately 5,000 meters and a nearly straight horizontal trajectory, was selected as a research target. F1 intersects a key oil and gas producing layer B (the target stratum in one embodiment of the present invention). The top interface of the oil and gas producing layer is approximately 5,000 meters deep, the bottom interface is approximately 10,000 meters deep, and the total thickness is approximately 5,000 meters. F1 intersects the oil and gas producing layer.

[0057] Step S12: Clarify the core, conventional logging, array acoustic logging and other data of the target strata in the study area.

[0058] A detailed statistical analysis of the core, conventional logging, and array acoustic logging data for each well drilled through the target formation in the study area was performed. If no well or core data, or only conventional logging data, were available for the target formation within the study area, the area was considered a Type I study area. If core data but no array acoustic logging data were available for the target formation within the study area, the area was considered a Type II study area. If core, conventional logging, and array acoustic logging data were available for the target formation within the study area, the area was considered a Type III study area. If conventional logging and array acoustic logging data but no core data were available for the target formation within the study area, the area was considered a Type IV study area.

[0059] For example, formation B in basin A has core data, conventional logging data, and array acoustic logging data, so it belongs to a type III study area.

[0060] Step S13: Solve the static rock mechanics parameters of the target strata in the study area.

[0061] For the Type II study area mentioned in step S12, rock mechanics tests are carried out on core samples to obtain rock parameters such as static Young's modulus, static Poisson's ratio and density. 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 S12, 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):

[0062]

[0063]

[0064] Where Ed is the dynamic Young's modulus, GPa; μd is the dynamic Poisson's ratio; Δtp is the longitudinal wave time difference DTC, μs / ft; Δts is the shear wave time difference DTS, μs / ft; ρ is the density, g / cm3.

[0065] 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.

[0066] For the Class I or Class IV study area mentioned in step S12, 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.

[0067] For example, 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 in the figure, a linear fitting method is used, and the goodness of fit R2 is 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 the figure. Figure 3As shown, the linear fitting method is also used, and the goodness of fit R2 is 0.86. Figure 2 and 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 formulas (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 1.

[0068] Table 1

[0069]

[0070]

[0071] Fitting the dynamic and static rock mechanics parameters in Table 1 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:

[0072]

[0073]

[0074] The goodness of fit for 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 / cm3.

[0075] Step S2: construct a geomechanical model based on the thickness of the target stratum in the study area, static rock mechanical parameters and the planar extension length of the vertical fault, and change the angle θ between the geomechanical model and any side surface of the vertical fault to output multiple different geomechanical models; and use finite element analysis computer software to unitize the geomechanical model to obtain a mathematical model; and set different boundary conditions and constraints for the mathematical model and carry out numerical simulation of the tectonic stress field to obtain simulation results corresponding to each experiment; the simulation results include: the vertical displacement of the top interface of the mathematical model and the displacement in the same direction as the strike of the vertical fault; 0°<θ<90°.

[0076] Step S2 specifically includes the following steps:

[0077] Step S21: The size of the regular tetrahedron geomechanical model is set according to the plane extension length of the vertical fault in the study area and the thickness of the target stratum. The vertical fault is set as a contact pair in the finite element analysis. Several different regular tetrahedron geomechanical models are set according to the plane extension orientation of the vertical fault.

[0078] The planar extension length of the vertical fault F1, obtained from step S11 and used as the research object, is defined as Lfault. The geomechanical model is set as a regular square prism, with the top and bottom interfaces being congruent squares with a side length equal to 1.2*Lfault. The height H of the regular square prism is equal to the thickness of the target stratum obtained from step S11. The vertical fault is defined as a cavity enclosed by four vertical planes of height H. These four planes can be divided into two groups, with planes belonging to different groups perpendicular to each other. Each group consists of two planes: the first group has planes spaced (Lfault) / 1000 apart and a planar extension length of Lfault; the second group has planes spaced (Lfault) / 1000 apart and a planar extension length of (Lfault) / 1000. The top and bottom surfaces of the existing basement fault region are both rectangles with a length of Lfault and a width of (Lfault) / 1000. The height of this region is H, and the volume is H*(Lfault)2 / 1000. The center of the intersection line of the pre-existing basement fault area and the top surface of the geomechanical model is consistent 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 is consistent with the center of the bottom surface of the geomechanical model.

[0079] The concept of contact pairs, used in contact analysis within the statics module of finite element analysis, was introduced and applied to the configuration of pre-existing faults in a three-dimensional geomechanical model. 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 passive deformation side of the regional tectonic stress field. The target surface is the surface undergoing active deformation and is therefore defined as the fault plane on the active deformation side of the regional tectonic stress field. Based on the regional stress direction associated with the target strata in the study area and the activity of fault F1, the target surface is set as the section closest to the regional stress direction, and the contact surface is set as the section away from the regional stress direction. Both planar extensions have a length of Lfault, a height of H, and a spacing of (Lfault) / 1000. The friction coefficient between the contact surface and the target surface is set based on the friction coefficient of vertical faults in strata with similar lithology to the target strata. Under the influence of regional tectonic stress, when the stress intensity at the fracture exceeds the stress intensity threshold for relative sliding between the contact and target surfaces, relative sliding occurs between the contact and target surfaces, thereby achieving relative displacement of the strata on either side of the pre-existing fault. After setting the geomechanical model dimensions and the vertical fault contact pair, the angle between the vertical fault and any side surface of the geomechanical model is set to θ. Based on the required accuracy, n values ​​θi (1≤i≤n) are selected for θ within the range of 0 to 90°, each corresponding to a different geomechanical model. This allows the study of the differences in deformation characteristics of the strata on either side of the fault when the angle between the vertical fault F1 and the regional horizontal stress varies.

[0080] In one embodiment of the present application, the three-dimensional geomechanical model is set as a regular quadrangular prism, and the length and width of the top surface are both 1.2*Lfault. Since Lfault is equal to 5000m obtained from step S11, the side lengths of the geomechanical model are all 6000m. The height H of the geomechanical model is equal to the thickness of the target stratum, which is 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 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×108m3. The center of the intersection line of the pre-existing basement fault area and the top surface of the geomechanical model is consistent 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 is consistent with the center of the bottom surface of the geomechanical model.

[0081] In one embodiment of the present invention, the direction of the horizontal maximum principal stress is fixed from south to north. Therefore, the southern section of the vertical fault is set as the target surface, and the northern section is set as the contact surface. Both have a planar extension length of 5000m, a height of 5000m, and a spacing of 5m. The friction coefficient of the pre-existing fault is assumed to be 0.15, so the friction coefficient between the contact surface and the target surface is set to 0.15. 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, and 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 nine different three-dimensional geomechanical models.

[0082] Step S22: The geomechanical model constructed in step S21 is meshed using tetrahedral elements in finite element analysis software to obtain a corresponding mathematical model, wherein 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 derived from the finite element analysis software.

[0083] The geomechanical model is meshed using triangular (2D) or tetrahedral (3D) elements, dividing the model into a series of elements and nodes. The meshing accuracy of the top surface of the geomechanical model is higher than that of the bottom surface and the four side surfaces. This is because, during the numerical simulation of the tectonic stress field, the displacement of the nodes on the top surface of the model along the X, Y, and Z coordinate axes is analyzed to demonstrate the model's vertical displacement and along the fault strike. After the mathematical model is constructed, the spatial coordinates (E_X, E_Y, E_Z) of all elements and nodes of the geomechanical model (N_X, N_Y, N_Z) are derived from the finite element analysis software.

[0084] In one embodiment of the present invention, 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 invention. 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 .

[0085] Step S23: Assign the static rock mechanical parameters (Young's modulus, Poisson's ratio and density) of the target stratum in the study area obtained from step S13 to each unit in the mathematical model, set different boundary conditions and constraints for the mathematical model, and carry out numerical simulation of the tectonic stress field to obtain the simulation results corresponding to each test.

[0086] According to the calculation results of the static rock mechanics parameters in step S13, static rock mechanics parameters (Young's modulus, Poisson's ratio and density) are assigned to all units of all geomechanical models. The three parameters of θ, the ratio of the regional horizontal maximum principal stress to the vertical stress (σH / σV), and the ratio of the regional horizontal minimum principal stress to the regional horizontal maximum principal stress (σh / σH) 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 S21, and the influence of σH / σV and σh / σH needs to be achieved by changing the boundary conditions of the mathematical model. According to the research accuracy requirements, the test data group of σH / σV and σh / σH is set, where σh / σH is less than 1. The boundary conditions and constraints of the mathematical model mainly include the regional stress application method, applied strength 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, numerical simulation of the tectonic stress field was carried out to obtain the simulation results corresponding to each experiment. The focus was on counting the vertical displacement of the top surface node of the model and the displacement in the same direction as the vertical fault in each numerical simulation of the tectonic stress field.

[0087] In one embodiment of the present invention, according to the calculation results of the static rock mechanical parameters in step S13, all units of all mathematical models are assigned the same static rock elastic parameters, specifically static Young's modulus (39.15GPa), static Poisson's ratio (0.232) and rock density (2.683g / cm3). When carrying out numerical simulation of the tectonic stress field for each mathematical model, σV is set to 100MPa 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 300MPa, to study the deformation characteristics of the model when σH / σV is different. σ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 the deformation characteristics of the model when σh / σH is different. Considering the comprehensiveness of this study and the efficiency of FEM-based numerical simulations, we also conducted additional 3D tectonic stress field simulations for σh / σH = 0.3 and σH / σV = 1.25, 1.50, 1.75, 2.25, 2.50, and 2.75 to more clearly demonstrate the effect of σH / σV on model deformation. Therefore, in one embodiment of the present invention, a total of 297 3D tectonic stress field numerical simulation experiments were involved.

[0088] like Figure 4 As shown, the displacement of the nodes on the XZ surface on the south side of the three-dimensional geomechanical model is unconstrained and is applied with σH, and the displacement of the nodes on the YZ surface on the east side is unconstrained and is applied with a horizontal minimum principal stress (σh). The displacement 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 is set to a constant value of 0, and the displacement in the Z direction is unconstrained. The displacement of the nodes on the XY surface on the upper side of the geomechanical model is unconstrained and a vertical stress (σV) is applied vertically downward. The displacement of the nodes on the XY surface on the lower side of the model in the Z direction is set to a constant value of 0, and the displacement in the X and Y directions is unconstrained. The entire three-dimensional geomechanical model is given a vertical downward gravity acceleration to make the simulation results closer to the characteristic of increasing vertical stress with increasing depth in actual geological conditions. Under the above boundary conditions and constraints, a numerical simulation of the tectonic stress field is carried out for each mathematical model to obtain the simulation results corresponding to each model, wherein the focus is on statistically analyzing the vertical displacement of the model top surface nodes and the displacement in the same direction as the vertical fault in each numerical simulation of the tectonic stress field. The boundary conditions of the mathematical model of an embodiment of the present invention and the corresponding vertical displacement of the model top surface nodes and the displacement in the same direction as the vertical fault are as follows: Figure 5 and 6 As shown. Among them, Figure 5 Figures (A) through (I) show the vertical deformation amplitude, and Figures (J) through (R) show the vertical deformation amplitude. Fault refers to the projection of a vertical fault on the top surface of the model, and the angle between fault and true north is set to 45°. The boundary conditions for the tectonic stress field simulation associated with each figure are indicated below each figure, where the regional maximum horizontal principal stress is set to 100 MPa. Figure 6 Figures (A) through (I) show the vertical deformation amplitude, and Figures (J) through (R) show the vertical deformation amplitude. Fault refers to the projection of a vertical fault on the top surface of the model, and the angle between fault and true north is set to 45°. The boundary conditions for the tectonic stress field simulation associated with each figure are indicated below each figure, where the ratio of the regional minimum horizontal principal stress to the regional maximum horizontal principal stress is always 0.3.

[0089] Step S3: Draw the vertical displacement DV curve and the horizontal displacement DH curve of the vertical fracture according to the simulation results corresponding to different boundary conditions and constraints, and construct a one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH; and calculate the linear fitting equations of the major axis and the minor axis according to the one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH corresponding to different boundary conditions and constraints.

[0090] Specifically, step S3 includes the following steps:

[0091] Step S31: Calculate the vertical displacement (DV) and horizontal displacement (DH) of the vertical fault corresponding to each set of tectonic stress field numerical simulation experiments based on the vertical displacement of the model top interface in step S23 and the displacement in the same direction as the vertical fault, and obtain a one-to-one correspondence between DH and DV, i.e., (DH, DV).

[0092] The vertical displacement (DV) of a certain position of the vertical fault is defined as the difference between the displacement of the node on the contact surface along the Z axis and the displacement of the node on the target surface along the Z axis. The horizontal displacement (DH) of a certain position of the vertical fault is defined as the difference between the displacement of the node on the contact surface along the fault direction (when the displacement of the node along the Y direction is positive, its displacement along the fault direction is defined as a positive value, otherwise it is defined as a negative value) and the displacement of the node on the target surface along the fault direction. When DH is a negative value, the vertical basement fault has a left-sliding characteristic along the direction, otherwise it has a right-sliding characteristic. The positive and negative values ​​of DH only represent the sliding direction of the fault along its direction, while the absolute value of DH represents the activity intensity of the fault along its direction. Therefore, in the present invention, the high value of DH is considered to be the high value of its absolute value, rather than the high value of the true value of DH. For example, for a left-sliding fault, a DH of -100m is a high value relative to a DH of -50m, and -100<-50 cannot be used to determine that the high value of DH is -50m. The vertical displacement (DV) and horizontal displacement (DH) curves of the vertical fracture corresponding to each mathematical model are calculated. The horizontal axis of the curve is the distance of the displacement data point from the end of one side of the fracture. Each position on the vertical fracture corresponds to a unique DV and DH, so a one-to-one correspondence between DH and DV corresponding to each mathematical model can be obtained, that is, (DH, DV).

[0093] In one embodiment of the present invention, the DV and DH curves corresponding to each mathematical model are as follows: Figure 7 and 8 As shown, the horizontal axis of each curve is the distance between the displacement data point and the southern end of the fault. However, in order to show the distance between the position where DV=0 and the center of the vertical fault, the distance between the displacement data point at the center of the vertical fault and the southern end of the fault is set to 0m, while the distance between the displacement data point at the southernmost end of the fault and the southern end of the fault is set to -Lfault / 2, i.e., -2500m, and the distance between the displacement data point at the northernmost end of the fault and the southern end of the fault is set to Lfault / 2, i.e., 2500m. Therefore, the horizontal axis range of each DV and DH curve is -2500 to 2500m. Figure 7In the figure, the boundary conditions of the numerical simulation of the tectonic stress field corresponding to each curve are marked on the figure. The horizontal axis of all curves is the distance along the fault strike. When the data point is located at the southernmost side of the fault, the distance is -2500m, when the data point is located at the center of the fault, the distance is 0m, and when the data point is located at the northernmost side of the fault, the distance is 2500m. Figure 8 The boundary conditions for the numerical simulation of the tectonic stress field corresponding to each curve are marked on the figure. The horizontal axis of all curves is the distance along the fault strike. When the data point is located at the southernmost side of the fault, the distance is -2500m, when the data point is located at the center of the fault, the distance is 0m, and when the data point is located at the northernmost side of the fault, the distance is 2500m.

[0094] Step S32: According to (DH, DV) obtained in step S31, find the vertical displacement DV(DH-MAX) corresponding to the maximum value DH-MAX of DH, and calculate the deflection angle α and the major axis of the axis-aligned ellipse accordingly.

[0095] According to (DH, DV) obtained in step S31, first find the maximum value of DH DH-MAX, then find the vertical displacement DV (DH-MAX) corresponding to this value, and calculate the deflection angle α according to formula (5):

[0096]

[0097] Where: DV(|DH-MAX|) is the DV corresponding to DH-MAX, m.

[0098] In the subsequent steps, (DH, DV) will be rotated counterclockwise around the coordinate origin by an angle α and then fitted with an axis-aligned ellipse implicit function. Here, the major axis of the axis-aligned ellipse can be calculated based on DH-MAX and DV(DH-MAX), that is:

[0099]

[0100] where γ is half the major axis of the axis-aligned elliptical implicit function, m.

[0101] In one embodiment of the present invention, taking the three-dimensional geomechanical model with σH / σV equal to 1, σh / σH equal to 0.3, and θ equal to 45° as an example, a plane coordinate system is constructed with DH as the horizontal axis and DV as the vertical axis. The "DH-DV" diagram obtained based on the three-dimensional geomechanical modeling and the numerical simulation of the tectonic stress field is presented as follows: Figure 9The set shown is composed of scattered points located in the second and third quadrants. Among them, there are always two scattered points near the origin of the coordinate system, representing the northeastern end and the southwestern end of the strike-slip fault (DH and DV are both 0). Starting from a scattered point near the origin of the coordinate system, the curve obtained by connecting all the scattered points clockwise resembles an ellipse whose major and minor axes are not parallel (perpendicular) to the coordinate axes. Because each independent variable (DH) does not correspond to a unique dependent variable (DV), DH and DV present an elliptical implicit function relationship closely related to the curve. According to formulas (5) and (6), the angle α between the major axis of the ellipse and the horizontal axis of the "DH-DV" coordinate system and the major axis 2γ of the ellipse can be calculated respectively.

[0102] Step S33: Rotate (DH, DV) obtained in step S31 counterclockwise by α around the coordinate origin to obtain (DH', DV'), and perform axis-aligned elliptical implicit function fitting on (DH', DV') to obtain the function expression.

[0103] After the data set (DH, DV) obtained in step S31 is rotated counterclockwise by α degrees around the coordinate origin, the data set (DH, DV) is converted into the data set (DH', DV'):

[0104] D' H =D H ·cosα-D V ·sinα (7)

[0105] D' V =D H ·sinα+D V ·cosα (8)

[0106] By fitting the implicit function of the data set (DH', DV'), we can obtain the axis-aligned ellipse implicit function with the major axis on the horizontal axis of the coordinate system, the minor axis parallel to the vertical axis of the coordinate system, and the axis passing through the origin of the coordinate system. Its expression is:

[0107]

[0108] Where γ and ω are the key constant terms of the axis-aligned elliptic implicit function. γ is equal to half of the absolute value of the minimum value of DH' in the data set (DH', DV'), which can also be calculated from formula (6), m; ω is equal to the maximum value of DV', m.

[0109] In one embodiment of the present invention, the data set (DH, DV) obtained in step S31 is rotated counterclockwise by α degrees around the coordinate origin to obtain the data set (DH', DV'), and the partial axis-aligned ellipse implicit function obtained after fitting is as shown in the attached figure. Figures 10-12 shown. Figure 10The numerical simulations of the tectonic stress fields corresponding to all figures are carried out when θ is equal to 30°, and other boundary conditions are marked in the upper left corner of each figure. Figure 11 The numerical simulations of the tectonic stress fields corresponding to all figures are carried out when θ is equal to 45°, and other boundary conditions are marked in the upper left corner of each figure. Figure 12 The numerical simulations of the tectonic stress fields corresponding to all figures are carried out when θ is equal to 60°, and other boundary conditions are marked in the upper left corner of each figure.

[0110] Step S34: Count the major axis and minor axis of the axis-aligned ellipse, and solve the linear fitting equations of the major axis and the minor axis.

[0111] The shape and position of the axis-aligned elliptical implicit function obtained in step S33 are determined by the parameters γ and ω. Each set of boundary conditions for the numerical simulation of the tectonic stress field corresponds to an axis-aligned elliptical implicit function. The γ and ω of all axis-aligned elliptical implicit functions are counted and linear fit is performed on these two parameters to obtain the linear fitting equations of γ and ω.

[0112] In one embodiment of the present invention, γ and ω of all axis-aligned elliptical implicit functions show an obvious linear positive correlation. Specific data are shown in Table 2. The fitting function is:

[0113] ω=0.4007·γ-0.6663 (10)

[0114] The goodness of fit of this equation as high as 0.9894 proves its applicability.

[0115] Table 2

[0116]

[0117]

[0118]

[0119]

[0120] Step S4: Calculate the vertical displacement DV curve of the target vertical fracture.

[0121] Specifically, step S4 includes: obtaining a DV curve of the vertical fracture as a research target according to the structural analysis results.

[0122] According to the method for calculating the vertical fracture DV curve described in step S31, the DV curve of the vertical fracture as the research object is obtained. The horizontal axis of the curve is the distance between the displacement data point and the southernmost side of the fracture. The minimum value of the horizontal axis is 0, and the maximum value is the total length of the fracture.

[0123] In one embodiment of the present invention, the DV curve of the vertical fracture F1 as the research object is shown in the attached figure. Figure 13 As shown in the attached figure, the plane trace of the vertical fault F1 is shown in the attached figure. Figure 14 As shown, the horizontal axis of the curve is the distance between the displacement data point and the southernmost side of the fault F1. The minimum value of the horizontal axis is 0, and the maximum value is the total length of the fault F1, 5000m.

[0124] Step S5: Calculate the major axis and minor axis of the axis-aligned elliptical implicit function of the target vertical fracture according to the DV curve of the target vertical fracture and the linear fitting equations of the major axis and minor axis to determine the expression of the axis-aligned elliptical implicit function and calculate the DH curve of the target vertical fracture.

[0125] Specifically, step S5 includes the following steps:

[0126] Step S51: Calculate the maximum horizontal displacement DH-MAX of the vertical fault as the research target and the corresponding DV(DH-MAX) based on methods such as horizontal offset of geological boundaries.

[0127] The maximum horizontal displacement DH-MAX of the vertical fault being studied is calculated using methods such as horizontal offset of geological boundaries. Examples of geological boundaries include abnormally high and low values ​​of stratum thickness, paleo-drainage systems, and paleo-faults. The calculated maximum horizontal displacement DH-MAX of the vertical fault should generally be located approximately at the center of the fault, as the displacement distribution pattern of vertical faults along their strike indicates that DH-MAX typically occurs at the center of the fault. The specific location of DH-MAX is the midpoint between two breakpoints on the vertical fault, resulting from the offset of the geological boundary. After determining DH-MAX and its corresponding location, the DV at the same location as DH-MAX, namely DV(DH-MAX), can be calculated from the DV curve of the vertical fault.

[0128] In one embodiment of the present invention, the vertical fault F1, the target of study, displaces the axis of a nearby syncline (serving as the displaced geological boundary). The axis of the displaced syncline intersects F1 at points P1 and P2, respectively. The distance between these two points is 490 m, representing the maximum left-lateral strike-slip distance DH-MAX of F1, which is 490 m. The distance from the midpoint of points P1 and P2 to the southernmost side of the fault is 2 km. Therefore, the position of the displacement data point corresponding to DH-MAX is 2 km, and the corresponding DV(DH-MAX) is 42.40 m. Figure 13 shown.

[0129] Step S52: Calculate the deflection angle and the major axis of the axis-aligned ellipse based on DH-MAX and DV(DH-MAX) of the vertical fracture as the research target.

[0130] According to formulas (5) and (6), and DH-MAX and DV(DH-MAX) calculated in step S51, the deflection angle α and the major axis of the axis-aligned ellipse can be calculated.

[0131] In one embodiment of the present invention, according to formulas (7) and (8), and DH-MAX and DV(DH-MAX) calculated in step S42, the tangent value of the deflection angle α is calculated to be 42.40 / 490, that is, 0.0865, and the deflection angle α is 4.95°. The major axis 2γ of the axis-aligned ellipse is equal to 491.831 m, and γ is equal to 241.915 m.

[0132] Step S53: Solve the expression of the axis-aligned elliptical implicit function applicable to the vertical fracture as the research target.

[0133] According to the linear fitting equations of the major axis and minor axis of the axis-aligned elliptical implicit function in step S34, the major axis and minor axis of the axis-aligned elliptical implicit function applicable to the vertical fracture as the research target are solved, and then the expression of the axis-aligned elliptical implicit function is determined.

[0134] In one embodiment of the present invention, according to formula (12) in step S34 and the major axis half-length γ of the axis-aligned elliptic implicit function calculated in step S51, the minor axis half-length ω of the axis-aligned elliptic implicit function is calculated to be 98.69 m. Then, according to formula (11), the expression of the axis-aligned elliptic implicit function applicable to fracture F1 can be obtained as follows:

[0135]

[0136] Step S54: The axis-aligned elliptical implicit function obtained in step S53 is rotated clockwise around the origin by a deflection angle α to restore it to the elliptical implicit function corresponding to the initial data set (DH, DV), and the DH curve of the vertical fracture as the research target is solved accordingly.

[0137] Substituting formulas (7) and (8) into the axis-aligned elliptic implicit function equation (9) yields:

[0138]

[0139] Here, (DH, DV) represents the data set before rotation by angle α, i.e., the original data set. Substituting the various parameters obtained from steps S52 and S53 into formula (12) yields a conversion model between DH and DV for the fracture being studied. The fracture's DH curve, obtained from step S41, can then be determined. Formula (12) requires the following conditions to ensure the unique solution of the DH curve: ① Both ends of the DH and DV curves are zero; ② When DH equals DH-MAX, DV equals DV(DH-MAX).

[0140] In one embodiment of the present invention, formulas (9) and (10) and the deflection angle α obtained in step S52 are substituted into formula (11) to obtain:

[0141]

[0142] For formula (13), the following two conditions are used to ensure that the DH curve of fault F1 has a unique solution: ① Both ends of the DH and DV curves are zero; ② When DH is equal to 490m, DV is equal to 42.4m.

[0143] According to the DV curve of the fracture F1 obtained in step S41, the DH curve of the fracture can be determined. The result is shown in the attached figure. Figure 15 shown.

[0144] In the above-mentioned method for calculating the horizontal displacement curve of a vertical fault, a geomechanical model is constructed based on the planar extension length of the vertical fault. Different boundary conditions and constraints are set, and numerical simulation of the tectonic stress field is performed. The vertical displacement of the top interface of the mathematical model under different conditions and the displacement in the same direction as the vertical fault are obtained. Based on this, a vertical displacement DV curve and a horizontal displacement DH curve are plotted, and a one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH is established. From this correspondence, the linear fitting equations of the major and minor axes are calculated. Based on the linear fitting equations of the major and minor axes, the DH curve of the target vertical fault is determined based on the vertical displacement DV curve of the actual target vertical fault. Thus, the present application obtains horizontal displacements at different locations. Furthermore, the present invention proposes to use numerical simulation of the tectonic stress field and contact pair settings to restore the deformation process of the vertical fault, thereby solving the linear fitting equations of the major and minor axes. Then, based on the linear fitting equations of the major and minor axes, the horizontal displacement DH curve is inverted based on the actual target vertical displacement DV curve, thereby improving the calculation accuracy of the horizontal displacement curve.

[0145] The advantages of this application include the following:

[0146] (1) Upgrading from one dimension to two dimensions: Compared with the prior art which usually calculates the horizontal displacement of a one-dimensional single point position of a vertical fault based on the horizontal offset distance of the geological boundaries on both sides of the vertical fault, the present invention proposes a method for determining the horizontal displacement curve of a vertical fault based on the vertical displacement curve of the vertical fault (easy to obtain) and the horizontal displacement of a single point position (relatively easy to obtain), thereby upgrading the method from only being able to obtain data of a one-dimensional "point" position to being able to obtain data of a two-dimensional "line" position.

[0147] (2) Reduction of human factors: Compared with the prior art which usually determines the horizontal displacement of vertical faults based on the human determination of the horizontal offset distance of geological boundaries, the method proposed in the present invention has fewer human factors in each link. Whether it is calculating the total length of the fault, calculating the thickness of the target layer, geomechanical modeling, mathematical modeling, numerical simulation of the tectonic stress field, or fitting the functional relationship of important parameters, all are completely quantitative. Different operators will obtain the same results when using the method proposed in the present invention.

[0148] (3) Innovatively using a numerical simulation method to reproduce the deformation of the strata on both sides of the vertical fault to carry out research on the horizontal displacement of the vertical fault: Compared with the existing technology that usually inverts the initial state of the strata on both sides of the vertical fault based on the current structural characteristics to obtain the horizontal displacement of the vertical fault, the method proposed in the present invention uses numerical simulation of the tectonic stress field and contact pair setting to restore the deformation process of the vertical fault. The vertical fault produces vertical displacement and displacement along its strike under the background of the regional tectonic stress field, which are reflected as vertical displacement and horizontal displacement respectively in the current structural form. The two have a close internal structural connection, and the method proposed in the present invention reproduces this process. This forward research method has not yet appeared in the study of the horizontal displacement of the vertical fault.

[0149] In one embodiment, a computer device is provided, which may be a server. The computer device includes a processor, memory, and a network interface connected via a system bus. The processor of the computer device is configured to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program, and a database. The internal memory provides an environment for the operation of the operating system and computer program in the non-volatile storage medium. The network interface of the computer device is configured to communicate with an external terminal via a network connection. When executed by the processor, the computer program implements the aforementioned method for calculating a horizontal displacement curve for a vertical fracture.

[0150] In one embodiment, a computer program product is further provided, including a computer program / instruction. When the computer program / instruction is executed by a processor, the computer program / instruction involves all or part of the process in the above embodiment method.

[0151] Those skilled in the art will appreciate that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to memory, storage, database or other media used in the embodiments provided in this application may include at least one of non-volatile and volatile memory. Non-volatile memory may include read-only memory (ROM), magnetic tape, floppy disk, flash memory or optical memory, etc. Volatile memory may include random access memory (RAM) or external cache memory. As an illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0152] 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.

[0153] The above-described embodiments merely represent several implementation methods of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that a person skilled in the art could make various modifications and improvements without departing from the spirit of the present application, all of which fall within the scope of protection of the present application. Therefore, the scope of protection of the present patent application shall be determined by the appended claims.

Claims

1. A method for calculating the horizontal displacement curve of a vertical fracture, characterized in that: The method comprises: The static rock mechanical parameters of the target strata in the study area are calculated based on the core, conventional logging and array acoustic logging data of the target strata in the study area; A geomechanical model is constructed based on the thickness of the target stratum in the study area, static rock mechanical parameters, and the planar extension length of the vertical fault, and the angle θ between the geomechanical model and any side surface of the vertical fault is changed to output multiple different geomechanical models; the geomechanical model is unitized using finite element analysis computer software to obtain a mathematical model; different boundary conditions and constraints are set for the mathematical model and a numerical simulation of the tectonic stress field is performed to obtain simulation results corresponding to each experiment; the simulation results include: the vertical displacement of the top interface of the mathematical model and the displacement in the same direction as the vertical fault; 0°<θ<90°; According to the simulation results corresponding to different boundary conditions and constraint conditions, the corresponding vertical displacement DV curve and horizontal displacement DH curve of the vertical fracture are drawn, and a one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH is established; and the linear fitting equations of the major axis and the minor axis are calculated according to the one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH corresponding to different boundary conditions and constraint conditions; Calculate the vertical displacement DV curve of the target vertical fracture; The major axis and minor axis of the axis-aligned elliptical implicit function of the target vertical fracture are calculated according to the DV curve of the target vertical fracture and the linear fitting equations of the major axis and minor axis to determine the expression of the axis-aligned elliptical implicit function and calculate the DH curve of the target vertical fracture.

2. The horizontal displacement curve method of vertical fracture according to claim 1, characterized in that: The geomechanical model is constructed based on the thickness of the target strata in the study area, static rock mechanical parameters and the plane extension length of the vertical fault, including: The size of the regular quadrangular prism geomechanical model is set according to the plane extension length of the vertical fault and the thickness of the target stratum in the study area. The cross section of the target vertical fault in the geomechanical model close to the regional stress direction is set as the target surface, and the cross section away from the regional stress direction is set as the contact surface.

3. The horizontal displacement curve method of vertical fracture according to claim 2, characterized in that: The method of unitizing the geomechanical model to obtain a mathematical model using finite element analysis computer software includes: meshing the geomechanical model using tetrahedral units to obtain a corresponding mathematical model, wherein the meshing density of the top interface of the mathematical model is higher than that of other positions.

4. The horizontal displacement curve method of vertical fracture according to claim 1, characterized in that: The different boundary conditions and constraint conditions include: regional stress application methods on the side and top surfaces of the mathematical model, and application strength and displacement constraint methods.

5. The horizontal displacement curve method of vertical fracture according to claim 4, characterized in that: The stress application intensity is mainly changed by the ratio of the regional horizontal maximum principal stress to the vertical stress (σH / σV) and the ratio of the regional horizontal minimum principal stress to the regional horizontal maximum principal stress (σh / σH).

6. The method for calculating the horizontal displacement curve of vertical fracture according to claim 1, characterized in that: The linear fitting equations for calculating the major axis and the minor axis according to the one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH corresponding to different boundary conditions and constraint conditions include: According to the one-to-one correspondence between the vertical displacement DV and the horizontal displacement DH corresponding to each different boundary condition and constraint condition, find the vertical displacement DV(DH-MAX) corresponding to the maximum value DH-MAX of the horizontal displacement DH, and calculate the deflection angle α and the major axis of the axis-aligned ellipse; The one-to-one correspondence between each set of vertical displacement DV and horizontal displacement DH is rotated counterclockwise by α around the coordinate origin to obtain (DH', DV'), and the axis-aligned ellipse implicit function fitting is performed on (DH', DV'), and the fitting function expression is calculated to obtain the minor axis of the axis-aligned ellipse; According to the major axis and minor axis of the axis-aligned ellipse corresponding to each boundary condition and constraint condition, the linear fitting equations of the major axis and minor axis are solved.

7. The horizontal displacement curve method of vertical fracture according to claim 1, characterized in that: Calculating the major axis and minor axis of the axis-aligned elliptical implicit function of the target upright fracture according to the DV curve of the target upright fracture and the linear fitting equations of the major axis and minor axis includes: Calculating the maximum horizontal displacement DH-MAX of the target vertical fracture and the corresponding vertical displacement DV(DH-MAX) according to the DV curve of the target vertical fracture; The deflection angle of the target vertical fracture and the major axis of the axis-aligned ellipse are solved according to the maximum horizontal displacement DH-MAX of the target vertical fracture and the corresponding vertical displacement DV(DH-MAX).

8. The method for calculating the horizontal displacement curve of vertical fracture according to claim 1, characterized in that: The DH curve for obtaining the target vertical fracture by calculation includes: According to the expression of the axis-aligned elliptical implicit function of the target vertical fracture, the elliptical implicit function corresponding to the initial data set (DH, DV) is restored after rotating the deflection angle α clockwise around the origin, and the DH curve of the target vertical fracture is solved accordingly.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of the method according to any one of claims 1 to 8 are implemented.

10. A computer program product comprising a computer program / instructions, characterized in that When the computer program / instructions are executed by a processor, the steps of the method according to any one of claims 1 to 8 are implemented.

Images