Method for reconstructing paleotectonic stress field based on vertical displacement of upright fracture

Abstract

The invention discloses a method for reconstructing a paleotectonic stress field based on vertical displacement of upright fracture, and aims to solve the problem of low precision of establishing the paleotectonic stress field in the prior art. The method comprises the following steps: solving static rock mechanical parameters of a target layer series; constructing a geomechanical model and a mathematical model related to sliding of the upright fracture along the trend; constructing a conversion model between L (DV = 0) and DV-MAX / DV-MIN; calculating L (DV = 0) and DV-MAX / DV-MIN according to a DV curve; and obtaining corresponding sigma H / sigma V, sigma h / sigma H and theta according to the L (DV = 0) and the DV-MAX / DV-MIN of the upright fracture and the conversion model, and reconstructing a regional paleotectonic stress field related to the movement of the upright fracture along the trend. Large personal errors possibly existing in typical structure recognition and matching between typical structures and theoretical modes are avoided, and the precision is remarkably improved.

Description

Technical Field

[0001] This application relates to the technical field of reservoir geomechanics in geology. Background Art

[0002] The study of the strike-slip along pre-existing nearly vertical faults and their associated tectonic fracture zones has very important significance for oil and gas geology and carbon sequestration (Bergmo, P.E.S., Grimstad, A.-A. & Lindeberg, E.G.B., 2011. Simultaneous CO2 injection and water production to optimise aquifer storage capacity, International Journal of Greenhouse Gas Control, 5, 555-564). Tectonic fractures caused by palaeo-tectonic stresses are important factors controlling the development of faults and tectonic fractures (Ding, W.L., Fan, T.L., Yu, B.S., Huang, X.B. & 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 of fault in craton Basin, Journal of Northwest University. Natural Science Edition, 52(6), 930-942). Therefore, studying the strike-slip along pre-existing nearly vertical faults and their associated tectonic fracture zones is actually studying the coupling relationship among faults - palaeo-tectonic stress fields - tectonic fractures. An important link in this is the reconstruction (determination) of the regional palaeo-tectonic stress field. Currently, the methods for reconstructing the regional palaeo-tectonic stress field based on the strike-slip characteristics of pre-existing nearly vertical fault zones mainly include:

[0003] Method 1: Use the azimuth constraint of the principal displacement zone in the positive relief and negative relief of the vertically fractured segments and the direction of the maximum horizontal principal stress (σH) in the horizontal region corresponding to the strike-slip along the strike of the vertical fracture (Li, Y.T., Ding, W.L., Zeng, T., Xia, W.Q., Cheng, X.Y., Shi, S.Y. & Ding, S.H., 2023. Structural geometry and kinematics of a strike-slip fault zone in an intracontinental thrust system: A case study of the No. 15 fault zone in the Fuling area, eastern Sichuan Basin, southwest China, Journal of Asian Earth Sciences, 242, 1 - 29; Sun, Q.Q., Fan, T.L., Holdsworth, R.E., Gao, Z.Q., Wu, J., Gao, S.C., Wang, M. & Yuan, Y.X., 2023. The spatial characterization of stepovers along deeply-buried strike-slip faults and their influence on reservoir distribution in the central Tarim Basin, NW China, Journal of Structural Geology, 170, 1 - 22);

[0004] Method 2: Use typical associated structures such as horsetail structures along the strike-slip of vertical fractures to indicate the kinematic characteristics of vertical fractures and the direction of regional horizontal stress (Deng, S., Li, H.L., Zhang, Z.P., Zhang, J.B. & Yang, X., 2019. Structural characterization of intracratonic strike-slip faults in the central Tarim Basin, AAPG Bulletin, 103(1), 109 - 137);

[0005] Method 3: The corresponding relationships between various secondary faults with different geometries and kinematics in the upright fault zone and the synthetic shear, antithetic shear, P-shear, Y-shear, T-fracture, and fold in the Riedel shear model, as well as the constraints of the corresponding strike-slip fracture patterns (including simple shear, transtensional, and transpressional) on the type and direction of regional tectonic stress (Deng, S., Zhao, R., Kong, Q. F., Li, Y. T., & Li, B., 2022. Two distinct strike-slip fault networks in the Shunbei area and its surroundings, Tarim Basin: Hydrocarbon accumulation, distribution, and controlling factors, Aapg Bulletin, 106(1), 77-102; Yao, Y. T., Zeng, L. B., Mao, Z., Han, J., Cao, D. S., & Lin, B., 2023. Differential deformation of a strike-slip fault in the Paleozoic carbonate reservoirs of the Tarim Basin, China, Journal of Structural Geology, 173, 1-12);

[0006] Method 4: Based on the fault slip data collected from outcrops in the field, computer software (such as Win_Tensor, T-Tecto, and FaultKin, etc.) and fault slip inversion methods are used to reconstruct the regional paleotectonic stress field (Nkodia, H.M.D., Boudzoumou, F., Miyouma, T., Kongota, E., Ganza, G.B., Lahogue, P. & Delvaux, D., 2024. Brittle faulting and tectonic stress history on the western margin of the Congo Basin between Kinshasa and Brazzaville: Implications for the evolution of the Malebo Pool and the Congo River, Tectonophysics, 877, 1-27; Dasgupta, S., Mukherjee, S., Vanik, N., Chatterjee, R. & Pal, S.K., 2023. Paleostress analysis and rift kinematics of the petroliferous Barmer rift basin, western Rajasthan, India, Marine and Petroleum Geology, 156, 1-29). Summary of the Invention

[0007] The present application provides a method for reconstructing the paleotectonic stress field based on the vertical displacement of vertical fractures, aiming to solve the problem of low accuracy in establishing the paleotectonic stress field in the prior art.

[0008] In a first aspect, a method for reconstructing the paleotectonic stress field based on the vertical displacement of vertical fractures is provided, including:

[0009] Step S1: Determine the study area, target formation series, and vertical fractures as the research objects, obtain the data of core, conventional logging, and array acoustic logging of the target formation series, and solve the static rock mechanical parameters of the target formation series;

[0010] Step S2: Construct a geomechanical model and a mathematical model related to the strike-slip of vertical fractures, conduct numerical simulation tests of the tectonic stress field, obtain the vertical deformation characteristics of the strata on both sides of the vertical fractures from the test results, and draw the vertical displacement D V curve of the vertical fractures corresponding to different models, and count the D VMaximum value D V-MAX and minimum value D V-MIN ratio D V-MAX / D V-MIN , as well as D V zero value position L(D V = 0);

[0011] Step S3, construct the conversion model between L(D V = 0) and D V-MAX / D V-MIN ;

[0012] Step S4, calculate L(D V = 0) and D V = 0) and D V-MAX / D V-MIN according to the D

[0013] Step S5, based on L(D V = 0) and D V-MAX / D V-MIN obtained from Step S4 for the vertical fracture, and the conversion model proposed in Step S3, obtain the corresponding σ H / σ V , σ h / σ H and θ, and reconstruct the regional paleotectonic stress field related to the strike activity of the vertical fracture accordingly.

[0014] In the above solution, optionally, the Step S1 includes:

[0015] For a study area with core data but without array acoustic logging data, conduct rock mechanics tests on core samples to obtain static rock mechanics parameters; among them, each static rock mechanics parameter takes the average value of the corresponding parameters of all samples; the static rock mechanics parameters include static Young's modulus, static Poisson's ratio, and density;

[0016] For a study area without core data, based on identifying the main lithology of the target formation, obtain the static rock mechanics parameters of the corresponding lithology in the nearby area by referring to the literature;

[0017] For a study area with core, conventional logging, and array acoustic logging data, first calculate the dynamic rock mechanics parameters according to the conventional logging and array acoustic logging, and the formula is:

[0018]

[0019] where E d represents the dynamic Young's modulus, with the unit of GPa; μ d represents the dynamic Poisson's ratio; Δt pDenote the longitudinal wave slowness as DTC, with the unit of μs / ft; Δt s Denote the shear wave slowness as DTS, with the unit of μs / ft; ρ represents the density, with the unit of g / cm 3 ; The dynamic rock mechanics parameters include the dynamic Young's modulus and the dynamic Poisson's ratio;

[0020] Then, based on the static rock mechanics parameters obtained from the rock mechanics tests of the core samples, carry out fitting work on the dynamic and static rock mechanics parameters to obtain the conversion model of the dynamic and static rock mechanics parameters of the target formation. According to the conversion model, convert the dynamic rock mechanics parameters into static rock mechanics parameters, and the conversion formula is:

[0021]

[0022]

[0023] In the formula, E d Denotes the dynamic Young's modulus, with the unit of GPa; μ d Denotes the dynamic Poisson's ratio; E represents the static Young's modulus, with the unit of GPa; μ represents the static Poisson's ratio.

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

[0025] According to the geological characteristics of the study area, construct a three-dimensional geomechanical model, which includes the position, strike, extension length of the vertical fault and the geometric characteristics of the target formation; Set the vertical fault as a discontinuous surface with a certain friction coefficient to simulate the relative sliding of the strata on both sides of the fault;

[0026] Convert the geomechanical model into a mathematical model, carry out mesh division through finite element analysis software, and divide the model into several elements and nodes;

[0027] Assign the static rock mechanics parameters obtained from step S1 to each element in the mathematical model; According to the requirements of the research accuracy, set the boundary conditions and constraint conditions of the mathematical model; Carry out numerical simulation experiments on the tectonic stress field to obtain the simulation results corresponding to each experiment;

[0028] Obtain the vertical displacement D of the vertical fault corresponding to each group of numerical simulation experiments on the tectonic stress field according to the vertical displacement of the top interface of the model V curve, and establish a one-to-one correspondence between the boundary conditions of the numerical simulation of the tectonic stress field and the vertical displacement curve D V ;

[0029] According to the vertical displacement D V curve, count the maximum value D of the vertical displacement D of each model V and the minimum value D V-MAX and minimum value DV-MIN Ratio D V-MAX / D V-MIN and vertical displacement D V zero value position L(D V = 0).

[0030] Optionally, in the above solution, step S3 includes:

[0031] When the ratio σ H / σ V of the regional horizontal maximum principal stress to the vertical stress takes different values, the zero value position L(D V = 0) and the ratio D V-MAX / D V-MIN of the maximum value to the minimum value obtained from all the models in step S2 are fitted to obtain the first conversion model of L(D V = 0) and D V-MAX / D V-MIN ;

[0032] When the angle θ between the vertical fracture and the direction of the regional horizontal maximum principal stress σ H takes different values, the effective data of L(D V = 0) and D V-MAX / D V-MIN obtained from all the models in step S2 are fitted to obtain the second conversion model of L(D V = 0) and D V-MAX / D V-MIN ;

[0033] When the ratio σ h / σ H of the regional horizontal minimum principal stress to the regional horizontal maximum principal stress takes different values, the effective data of L(D V = 0) and D V-MAX / D V-MIN obtained from all the models in step S2 are fitted to obtain the third conversion model of L(D V = 0) and D V-MAX / D V-MIN ;

[0034] In the above solution, further optionally, constructing the first, second, and third conversion models further includes:

[0035] When σ H / σ V takes different values, the effective data of L(D V = 0) and D V-MAX / D V-MIN are fitted. If multiple different σ H / σV For the numerical values, multiple corresponding fitting functions are obtained through fitting and denoted as the first conversion model;

[0036] In the first conversion model, all the fitting functions belong to the same type, but their constant coefficients are different; fit the functional relationship between the constant coefficient in each fitting function and σ H / σ V If the fitting function is not clear, then according to the one-to-one correspondence between the constant coefficient and σ H / σ V draw the first correlation curve as a discrimination chart, where the horizontal axis of the first correlation curve is σ H / σ V , and the vertical axis is each constant coefficient; for each undetermined constant coefficient, draw a curve;

[0037] When θ takes different numerical values, fit the valid data of L(D V =0) and D V-MAX / D V-MIN If multiple different θ numerical values are set, multiple corresponding fitting functions are obtained through fitting and denoted as the second conversion model;

[0038] In the second conversion model, all the fitting functions belong to the same type, but their constant coefficients are different; fit the functional relationship between the constant coefficient in each fitting function and θ. If the fitting function is not clear, then draw the second correlation curve as a discrimination chart according to the one-to-one correspondence between the constant coefficient and σ H / σ V where the horizontal axis of the second correlation curve is σ H / σ V , and the vertical axis is each constant coefficient; for each undetermined constant coefficient, draw a curve;

[0039] When σ h / σ H takes different numerical values, fit the valid data of L(D V =0) and D V-MAX / D V-MIN If multiple different σ h / σ H numerical values are set, multiple corresponding fitting functions are obtained through fitting and denoted as the third conversion model;

[0040] In the third conversion model, all the fitting functions belong to the same type, but their constant coefficients are different. Fit the functional relationship between the constant coefficient in each fitting function and σ h / σ H If the fitting function is not clear, then according to the constant coefficient and σ h / σH The one-to-one correspondence is used to draw the third correlation curve as a discrimination chart. Among them, the horizontal axis of the third correlation curve is σ h / σ H , and the vertical axis is each constant coefficient. For each undetermined constant coefficient, a curve is drawn.

[0041] In the above solution, further optionally, constructing the first, second, and third conversion models further includes:

[0042] When σ H / σ V takes different values, the effective data of L(D V =0) and D V-MAX / D V-MIN is fitted, and the fitting formula is:

[0043] Y = a - b·c X

[0044] Among them, X represents D V-MAX / D V-MIN , Y represents L(D V =0), and a, b, and c all represent constant coefficients; using the parameter σ H / σ V to specify the exponential function model between the parameters X and Y as the first conversion model:

[0045]

[0046] In the formula, X represents D V-MAX / D V-MIN , Y represents (D V =0); a 1 , b 1 and c 1 represent the constant coefficients of the first conversion model;

[0047] When θ takes different values, use the fitting formula to fit the effective data of L(D V =0) and D V-MAX / D V-MIN ;

[0048] Use the parameter θ to specify the fitting formula as the second conversion model:

[0049]

[0050] Among them, X represents D V-MAX / D V-MIN , Y represents L(D V =0); a 2 , b 2 and c 2Denote the constant coefficient of the second conversion model;

[0051] When σ h / σ H takes different numerical values, use the fitting formula to fit the L(D V =0) and D V-MAX / D V-MIN valid data;

[0052] Use the parameter σ h / σ H to specify the fitting formula as the third conversion model:

[0053]

[0054] where X represents D V-MAX / D V-MIN , Y represents L(D V =0); a 3 , b 3 and c 3 represent the constant coefficients of the third conversion model.

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

[0056] Respectively, with σ H / σ V , θ and σ h / σ H as the horizontal axis and L(D V =0) as the vertical axis, draw the corresponding discriminant curves;

[0057] According to the L(D V =0) obtained from step S4, draw a horizontal line with a constant ordinate equal to L(D V =0) on the corresponding discriminant curve, and based on the intersection point of this horizontal line and the discriminant curve, reconstruct the paleotectonic stress field related to the strike activity of the vertical fracture.

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

[0059] With σ H / σ V as the horizontal axis and L(D V =0) as the vertical axis, draw the first discriminant curve; if the horizontal line L(D V =0) intersects the first discriminant curve, record this intersection point as the first intersection point; if the horizontal line L(D V =0) does not intersect the first discriminant curve, then at L(D V =0) with respect to σ H / σ VFind the point closest to the horizontal line L(D V = 0) on the function curve graph. At this time, mark this point closest to the distance as the first intersection point;

[0060] Taking θ as the horizontal axis and L(D V = 0) as the vertical axis, draw the second discrimination curve; if this horizontal line L(D V = 0) has an intersection point with the second discrimination curve graph, mark it as the second intersection point; if this horizontal line L(D V = 0) has no intersection point with the second discrimination curve graph, then find the point closest to the horizontal line L(D V = 0) on the function curve graph of L(D V = 0) with respect to θ, and mark it as the second intersection point;

[0061] Taking σ h / σ H as the horizontal axis and L(D V = 0) as the vertical axis, draw the third discrimination curve; if this horizontal line L(D V = 0) has an intersection point with the third discrimination curve graph, mark it as the third intersection point; if this horizontal line L(D V = 0) has no intersection point with the third discrimination curve graph, then find the point closest to the horizontal line L and mark it as the third intersection point on the function curve graph of L(D V = 0) with respect to σ h / σ H ;

[0062] The abscissas of the first intersection point, the second intersection point, and the third intersection point respectively represent σ H / σ V , θ, and σ h / σ H of the regional tectonic stress field when the vertical fracture as the research object undergoes sliding along its strike, where σ H / σ V and σ h / σ H reflect the relative magnitudes of the regional horizontal maximum principal stress, the regional horizontal minimum principal stress, and the vertical stress, reflecting the type of tectonic stress; θ reflects the angle between the regional horizontal maximum principal stress and the vertical fracture. Combining with the direction of the strike-slip of the vertical fracture, the precise azimuth angle of the regional horizontal maximum principal stress can be quantitatively determined, and thus the reconstruction of the paleotectonic stress field has been completed.

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

[0064] Based on further analysis and research of the problems in the prior art, it is recognized that compared with the prior art which usually uses the azimuth constraint region of the horizontal maximum principal stress direction in the undulating section of the vertical fracture, uses typical structures such as horsetail structures to indicate the regional horizontal stress direction, and uses the corresponding relationship between each secondary fracture in the vertical fracture zone and each element in the Riedel shear model and the strike-slip fracture mode to constrain the type and direction of regional tectonic stress, the method proposed in this application does not rely on all of the above typical structures, but only uses the vertical displacement amount of the vertical fracture, with a very high degree of quantification and a very low human decision-making factor. Therefore, it avoids the large human errors that may exist in the identification of typical structures and their matching with theoretical models, and significantly improves the accuracy.

[0065] Compared with the prior art which usually can only limit the range of the horizontal maximum principal stress direction in the region and the relative strength of the regional stress, the method proposed in this application can accurately determine the azimuth angle of the horizontal maximum principal stress in the region, as well as the ratio of the horizontal maximum principal stress, the horizontal minimum principal stress and the vertical stress in the region, achieving complete quantification.

[0066] Compared with the prior art which is usually only applicable when the tectonic activity is weak after the sliding deformation along the strike of the vertical fracture, because in this case some typical structures can be well preserved, but when the typical structures related to strike-slip fractures are damaged due to strong regional compressive stress, the applicability of the prior art is greatly reduced, while the method proposed in this application has improved adaptability in this case because the deformation degree of the strata on both sides of the vertical fracture increases, and the vertical deformation amplitude and vertical displacement of the vertical fracture are more significant. Brief Description of the Drawings

[0067] Figure 1 It is a schematic flow chart of a method for reconstructing the paleotectonic stress field based on the vertical displacement amount of a vertical fracture provided by an embodiment of this application.

[0068] Figure 2 It is a linear conversion mathematical model of AC and DTC obtained from the AC and DTC data of the target formation in the study area provided by an embodiment of this application.

[0069] Figure 3 It is a linear conversion mathematical model of DTC and DTS obtained from the DTC and DTS data of the target formation in the study area provided by an embodiment of this application.

[0070] Figure 4The size setting, vertical fracture setting, boundary condition setting of the geomechanical model provided by an embodiment of this application, and nine different geomechanical models corresponding to setting the included angle θ between the vertical fracture and one side of the geomechanical model to nine different values, and the corresponding mathematical models (the number of elements and nodes of the mathematical models are also marked below each model).

[0071] Figure 5 The vertical deformation amplitude of the nodes on the top interface of the model obtained from the numerical simulation results of the tectonic stress field provided by an embodiment of this application. fault refers to the projection of the vertical fracture on the top interface of the model, and the included angle between fault and the due north direction is set to 45°. The boundary conditions of the tectonic stress field simulation related to each figure are marked below each figure, where the regional horizontal maximum principal stress is set to 100 MPa.

[0072] Figure 6 The vertical deformation amplitude of the nodes on the top interface of the model obtained from the numerical simulation results of the tectonic stress field provided by an embodiment of this application. fault refers to the projection of the vertical fracture on the top interface of the model, and the included angle between fault and the due north direction is set to 45°. The boundary conditions of the tectonic stress field simulation related to each figure are marked below each figure, where the ratio of the regional horizontal minimum principal stress to the regional horizontal maximum principal stress is 0.3.

[0073] Figure 7 The vertical displacement amount (D V ) curve at the top of the vertical fracture F1 obtained from the numerical simulation results of the tectonic stress field provided by an embodiment of this application. The boundary conditions of the tectonic stress field numerical simulation corresponding to each curve are marked on the figure. The horizontal axis of all curves is the distance along the strike of the fracture. When the data point is located at the southernmost side of the fracture, this distance is -2500 m. When the data point is located at the center of the fracture, this distance is 0 m. When the data point is located at the northernmost side of the fracture, this distance is 2500 m.

[0074] Figure 8 The provided by an embodiment of this application when σ H / σ V are equal to 1, 2, and 3 respectively, the fitting functions of D V-MAX / D V-MIN and L(D V = 0), and all fitting functions are fitted with equation (5), where X refers to D V-MAX / D V-MIN , and Y refers to L(D V = 0).

[0075] Figure 9For an embodiment of the present application, when θ is equal to 10°, 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80° respectively, D V-MAX / D V-MIN and the fitting functions of L(D V = 0), all the fitting functions are fitted with Equation (5), where X refers to D V-MAX / D V-MIN and Y refers to L(D V = 0).

[0076] Figure 10 Curves of constants a, b and c in Equation (7) with respect to θ provided for an embodiment of the present application.

[0077] Figure 11 For an embodiment of the present application, when σ h / σ H is equal to 0.1 to 0.9 at intervals of 0.1 respectively, D V-MAX / D V-MIN and the fitting functions of L(D V = 0), all the fitting functions are fitted with Equation (5), where X refers to D V-MAX / D V-MIN and Y refers to L(D V = 0).

[0078] Figure 12 Curves of constants a, b and c in Equation (8) with respect to σ h / σ H provided for an embodiment of the present application.

[0079] Figure 13 Vertical displacement curve (A) of the upright fracture F1 as the research object provided for an embodiment of the present application, curves of functions of σ V-MAX / σ V-MIN , σ H / σ V , σ h / σ H and θ with respect to L(D V = 0) (Curve B), and the reconstruction result (C) of the paleotectonic stress field in the area related to the strike-slip of the upright fracture F1. Detailed implementation manners

[0080] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0081] In the description of the present application: Unless otherwise specified, the meaning of "a plurality of" is two or more. Terms such as "first", "second", "third", etc. in the present application are intended to distinguish the objects being referred to, and do not have special significance in terms of technical connotations (for example, they should not be understood as emphasizing importance or order, etc.). Expressions such as "including", "comprising", "having", etc. also mean "not limited to" (certain units, components, materials, steps, etc.).

[0082] There are still some problems with the methods mentioned in the background art:

[0083] Problem 1: Method 1 and Method 2 can generally only be used to limit the range of the horizontal stress direction (azimuth angle) in a defined area and cannot obtain relatively accurate values;

[0084] Problem 2: In Method 3, the research on the corresponding relationship between each fracture and the typical structures in the Riedel shear model is affected by the subjective judgment of the researcher, and there is a strong multi-solution property;

[0085] Problem 3: In Method 3 and Method 4, when the typical structures related to strike-slip fractures are damaged by strong regional tectonic movements that occurred after the strike-slip fractures, it is difficult to reconstruct the regional paleo-tectonic stress field through the Riedel shear model or rich fault-slip data.

[0086] The strong regional compressive stress causes the typical structures related to strike-slip fractures to be damaged, which is an unfavorable factor. However, it causes strong deformation of the strata on both sides of the vertical fracture, and this strong deformation is reflected as a very obvious vertical deformation amplitude and vertical displacement amount of the strata on both sides of the vertical fracture. The quantitative research on this strong deformation may be of great significance for reconstructing the regional paleo-tectonic stress field.

[0087] The purpose of the present application is to provide a method for reconstructing the paleo-tectonic stress field based on the vertical displacement of a vertical fracture, so as to solve the problem that it is difficult to quantitatively predict the regional paleo-tectonic stress field related to the strike-slip along the strike of the vertical fracture through tectonic analysis when the degree of damage to the typical structures related to the strike-slip along the strike of the vertical fracture in the prior art is relatively high, and to promote the process of oil and gas field exploration and development and the research on reservoir geomechanics.

[0088] Therefore, in an embodiment of the present application, a method for reconstructing the paleo-tectonic stress field based on the vertical displacement of a vertical fracture is provided, and this method includes the following steps S1 - step S5.

[0089] Step S1, determine the study area, the target formation, and the vertical fracture as the research object, and clarify the data conditions of the core, conventional logging, and array acoustic logging of the target formation, and solve the static rock mechanical parameters of the target formation in the study area.

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

[0091] Select the study area and target formation series according to the actual research requirements of oil and gas exploration and development, clarify the lithology of the target formation series and the depth distribution characteristics of its top and bottom interfaces in the study area. Identify the vertical faults for which paleotectonic stress field reconstruction is to be carried out. These vertical faults are required to completely penetrate the target formation series, and the intersection lines with the top and bottom interfaces of the target formation series should be continuous and clearly visible on the plane.

[0092] In an embodiment of the present application, a vertical fault F1 with a left-lateral strike-slip along its strike and a planar extension length of approximately 6.834 km and an almost straight planar trace in Basin A of China is selected as the research target. This fault F1 penetrates an important oil and gas producing formation B (as the target formation series in an embodiment of the present application). When carrying out comprehensive stratigraphic-structural interpretation based on 3D seismic data in this area, the conversion of two-way travel time and actual depth is not performed, so the depth on a time scale is used. The depth of the top interface of this oil and gas producing formation is approximately 2400 ms, the depth of the bottom interface is approximately 2600 ms, and the total thickness is approximately 200 ms. The vertical fault F1 penetrates this oil and gas producing formation.

[0093] Step S12: Understand the core, conventional logging, array acoustic logging, and other data of the target formation series in the study area.

[0094] Conduct a detailed statistics on the core, conventional logging, and array acoustic logging of each well in the study area that penetrates the target formation series. If there are no well data and core data for the target formation series in the study area, or only conventional logging data are available, it is identified as a type I study area; if the target formation series in the study area has core data but no array acoustic logging data, it is identified as a type II study area; if the target formation series in the study area has core, conventional logging, and array acoustic logging data, it is identified as a type III study area; if the target formation series in the study area has conventional logging and array acoustic logging data but no core data, it is identified as a type IV study area.

[0095] In an embodiment of the present application, 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.

[0096] Step S13: Solve the static rock mechanical parameters of the target formation series in the study area.

[0097] For the Type II study area mentioned in step S12, conduct rock mechanics tests 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, take the average of the corresponding static parameters of all samples. For the Type III study area mentioned in step S12, first calculate the dynamic rock mechanics parameters based on conventional logging and array acoustic logging, and then convert the dynamic rock mechanics parameters to static rock mechanics parameters according to parameters such as the static Young's modulus and static Poisson's ratio obtained from the rock mechanics tests of the core samples. Specifically as follows: The acoustic transit time (AC) logging data, longitudinal wave transit time (DTC), and shear wave transit time (DTS) logging data in the conventional logging data belong to transit time data and are considered to have good correlation. Therefore, first construct a conversion model between AC and DTC, and then construct a conversion model between DTC and DTS according to the array acoustic logging data; according to 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):

[0098]

[0099] Where E d is the dynamic Young's modulus, GPa; μ d is the dynamic Poisson's ratio; Δt p is the longitudinal wave transit time DTC, μs / ft; Δt s is the shear wave transit time DTS, μs / ft; ρ is the density, g / cm 3 .

[0100] According to the wells and depths of the core samples for which the rock mechanics tests are carried out, obtain the dynamic rock mechanics parameters at the corresponding positions, and then combine the static rock mechanics parameters obtained from the rock mechanics tests to construct a conversion model between the dynamic rock mechanics parameters and the static rock mechanics parameters. According to this conversion model, the calculation results of the dynamic rock mechanics parameters of all wells in the study area in the target formation can be converted into static rock mechanics parameters. Solve the average value of the static rock mechanics parameters of all wells in the study area in the target formation.

[0101] For the Type I or Type IV study area mentioned in step S12, on the basis of identifying the main lithology of the target formation in the study area, obtain the static rock mechanics parameters of the corresponding lithology in the nearby area by referring to the literature.

[0102] In an embodiment of the present application, Basin A belongs to the Type III study area. Therefore, construct a conversion model between AC and DTC according to the AC and DTC data of all wells in the study area, as shown in the appendix Figure 2 shown, using the method of linear fitting, the goodness of fit R 2is 0.66. Then, based on the DTC and DTS data of all the wells in the study area, a conversion model between DTC and DTS is constructed. As shown in the appendix Figure 3 , a linear fitting method is also adopted, and the goodness of fit R 2 is 0.86. According to the conversion models shown in appendix Figure 2 and Figure 3 , all the AC data in the study area can be converted into DTC and DTS data. Then, according to formulas (1) and (2), the dynamic Young's modulus and dynamic Poisson's ratio can be calculated. The static rock mechanical parameters can be obtained from the rock mechanical experiments of the rock samples in the target formation in the study area. The corresponding relationships between these static rock mechanical parameters and the calculation results of the dynamic rock mechanical parameters at the corresponding wells and corresponding depths are shown in Table 1.

[0103] Table 1

[0104]

[0105] By fitting the dynamic and static rock mechanical parameters in Table 1, the conversion models of the dynamic and static rock mechanical parameters of the target formation can be obtained, as shown in formulas (3) and (4) respectively:

[0106]

[0107] The goodness of fit of formulas (3) and (4) is 0.62 and 0.63 respectively. According to formulas (3) and (4) and the calculation results of the dynamic rock mechanical parameters of the target formation in the study area, the static rock mechanical parameters of the target formation 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 .

[0108] Step S2: Construct a geomechanical model and a mathematical model related to the strike-slip along the vertical fault, and conduct a numerical simulation test of the tectonic stress field. Obtain the vertical deformation characteristics of the strata on both sides of the vertical fault from the test results and draw the vertical displacement amount (D V ) curves of the vertical fault corresponding to different models, and count the ratio of the maximum value (D V ) to the minimum value (D V-MAX ) of D V-MIN for each model, as well as the zero value position of D V (L(D V = 0)).

[0109] Step S21: Set the size of the regular quadrangular prism geomechanical model, set the vertical fault as a contact pair in the finite element analysis, and set several different regular quadrangular prism geomechanical models according to the planar extension orientation of the vertical fault.

[0110] The geomechanical model is set as a regular quadrangular prism, and the top and bottom interfaces of the regular quadrangular prism are congruent squares. To facilitate the subsequent construction of the mathematical model and result analysis, the side length of the top and bottom interfaces of the regular quadrangular prism is set to 6000 m, and the height H of the regular quadrangular prism is set to 5000 m. The vertical fracture is set as a cavity enclosed by four vertical planes with a height of H. 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 spacing of the first group of planes is set to 5 m, and the plane extension length is 5000 m; the spacing of the second group of planes is 5000 m, and the plane extension length is 5 m. The top and bottom surfaces of the pre-existing basement fracture area are both rectangles with a length of 5000 m and a width of 5 m. The height of this area is H, and the volume is 1.25×10 8 m 3 . The pre-existing basement fracture area penetrates the geomechanical model from top to bottom. The center of the intersection line of this 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 fracture and the bottom surface of the geomechanical model coincides with the center of the bottom surface of the geomechanical model.

[0111] The concept of contact pair in contact analysis for the static module in finite element analysis is introduced and used for the setting of pre-existing fractures in the three-dimensional geomechanical model. The contact pair consists of two planes (or curved surfaces), namely the contact surface and the target surface. Among them, the contact surface is the surface that undergoes passive deformation, so it is defined as the fault surface on the side that undergoes passive deformation in the regional tectonic stress field, while the target surface is the surface that undergoes active deformation, so it is defined as the fault surface on the side that undergoes active deformation in the regional tectonic stress field. According to the regional stress direction related to the target formation system and the activity of fracture F1 in the study area, the fault surface on the side close to the regional stress direction is set as the target surface, and the fault surface on the side far from the regional stress direction is set as the contact surface. The plane extension lengths of both are 5000 m, the heights are both 5000 m, and the spacing is 5 m. According to the friction coefficient of the vertical fracture in the formation close to the lithology of the target formation system, the friction coefficient between the contact surface and the target surface is set. Under the action of the regional tectonic stress, when the stress intensity at the fault surface exceeds the stress intensity threshold for relative sliding between the contact surface and the target surface, the contact surface and the target surface can undergo relative sliding, thereby realizing the relative displacement of the strata on both sides of the pre-existing fracture. After completing the size setting of the geomechanical model and the contact pair setting of the vertical fracture, the angle between the vertical fracture and any side Surface of the geomechanical model is set as θ. According to the research accuracy requirements, n values of θ are selected within 0 - 90° for θ i (1≤i≤n), and each θ i corresponds to 1 different geomechanical model. Thus, the deformation characteristic differences of the strata on both sides of the fracture can be studied when the angle between the vertical fracture F1 and the regional horizontal stress is different.

[0112] In an embodiment of the present application, the three-dimensional geomechanical model is set as a regular quadrangular prism, the length and width of the top surface are both 6000 m, and the height H of the geomechanical model is equal to 5000 m. The vertical fracture is set as a cavity enclosed by four vertical planes with a height of 5000 m. 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 spacing of the first group of planes is 5 m, and the plane extension length is 5000 m; the spacing of the second group of planes is 5000 m, and the plane extension length is 5 m. The top and bottom surfaces of the pre-existing basement fracture region are both rectangles with a length of 5000 m and a width of 5 m. The height of this region is 5000 m, and the volume is 1.25×10 8 m 3 . The pre-existing basement fracture region penetrates the geomechanical model from top to bottom. The center of the intersection line of this region 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 fracture and the bottom surface of the geomechanical model is consistent with the center of the bottom surface of the geomechanical model.

[0113] In an embodiment of the present application, the direction of the horizontal maximum principal stress is fixed from south to north. Therefore, the south side section of the vertical fracture is set as the target surface, and the north side section is set as the contact surface. The plane extension lengths of both are 5000 m, the heights are both 5000 m, and the spacing is 5 m. The friction coefficient of the pre-existing fracture is considered to be 0.15. Therefore, 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 relative displacement of the strata on both sides of the pre-existing fracture can be realized. θ is set to 10°, 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80°, corresponding to 9 different three-dimensional geomechanical models.

[0114] Step S22, in the finite element analysis software, mesh the geomechanical model constructed in step S21 with tetrahedral elements to obtain the corresponding mathematical model, where the mesh division density of the top interface of the model is higher than other positions of the model. After constructing the mathematical model, export the coordinates of all elements and nodes of the geomechanical model from the finite element analysis software.

[0115] Mesh the geomechanical model with triangular (two-dimensional) or tetrahedral (three-dimensional) elements to divide the geological model into a series of elements and nodes. The mesh division 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 in the present application, the vertical deformation amplitude (vertical displacement) of the nodes on the top surface of the model is analyzed emphatically. After the mathematical model is constructed, export the spatial coordinates (E_X, E_Y, E_Z) of all elements and the spatial coordinates (N_X, N_Y, N_Z) of the nodes of the geomechanical model from the finite element analysis software.

[0116] In one embodiment of the present application, tetrahedral elements are used to mesh a three-dimensional geomechanical model, dividing the geological model into a series of elements and nodes. The ANSYS computer software is used to mesh the geomechanical model. Based on the characteristics of the elements, the SOLID 45 element is selected to mesh the three-dimensional geomechanical model in one embodiment of the present application. The improvement of the meshing accuracy is achieved by the method of free meshing, 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 and the spatial coordinates (N_X, N_Y, N_Z) of the nodes of the geomechanical model are exported from the finite element analysis software. The total number of elements and the total number of nodes after meshing each three-dimensional geomechanical model in step S21 are shown in the appendix Figure 4 .

[0117] Step S23: Assign the static rock mechanical parameters (Young's modulus, Poisson's ratio, and density) of the target formation in the study area obtained by solving in step S13 to each element in the mathematical model, set different boundary conditions and constraint conditions for the mathematical model, and conduct numerical simulation of the tectonic stress field to obtain the simulation results corresponding to each test.

[0118] According to the calculation results of the static rock mechanical parameters in step S13, assign the static rock mechanical parameters (Young's modulus, Poisson's ratio, and density) to all elements of all geomechanical models. Define 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 ) as the key parameters affecting the numerical simulation results of the tectonic stress field. The influence of θ has been reflected in step S21, while 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, set the test data groups of σ H / σ V and σ h / σ H , where σ h / σ HLess than 1. The boundary conditions and constraint conditions of the mathematical model mainly include the regional stress application method, application intensity, and displacement constraint method on the side and top surfaces of the mathematical model. After setting reasonable boundary conditions and constraint conditions for each mathematical model, numerical simulation of the tectonic stress field is carried out to obtain the simulation results corresponding to each test. The key is to statistically analyze the vertical displacement of the nodes on the top surface of the model in each numerical simulation of the tectonic stress field.

[0119] In an embodiment of the present application, according to the calculation results of the static rock mechanics parameters in step S13, the same static rock elastic parameters are assigned to all elements of all mathematical models, specifically the static Young's modulus (39.15 GPa), static Poisson's ratio (0.232), and rock density (2.683 g / cm 3 ). When carrying out numerical simulation of the tectonic stress field for each mathematical model, σ V is set to 100 MPa to study the deformation characteristics of the model under the action of the same overlying formation gravity. σ H / σ V is set to 1, 2, and 3, namely 100, 200, and 300 MPa, 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 numerical simulation operations based on FEM, numerical simulations of the three-dimensional tectonic stress field are also supplemented when σ h / σ H is equal to 0.3 and σ H / σ V is equal to 1.25, 1.50, 1.75, 2.25, 2.50, and 2.75 to more clearly show the influence of σ H / σ V on the deformation of the model. Therefore, in an embodiment of the present application, a total of 297 three-dimensional tectonic stress field numerical simulation experiments are involved.

[0120] As shown in the appendix Figure 4 , the displacement of the nodes on the XZ plane on the south side of the three-dimensional geomechanical model is not constrained, and σ H is applied. The displacement of the nodes on the YZ plane on the east side is not constrained, and the horizontal minimum principal stress (σ h)。 The displacements of the nodes on the XZ plane on the north side and the YZ plane on the west side of the geomechanical model in the X and Y directions are set to a fixed value of 0, and the displacement in the Z direction is unconstrained. The displacements of the nodes on the XY plane on the upper side of the geomechanical model are unconstrained and a vertical stress (σ V ) is applied vertically downward. The displacement of the nodes on the XY plane on the lower side of the model in the Z direction is set to a fixed value of 0, and the displacements in the X and Y directions are unconstrained. The entire three-dimensional geomechanical model is given a vertically downward gravitational acceleration so that the simulation results are closer to the characteristic that the vertical stress increases with the increase of depth in the actual geological situation. Under the above boundary conditions and constraint conditions, numerical simulations of the tectonic stress field are carried out for each mathematical model to obtain the corresponding simulation results. Among them, the key is to statistically analyze the vertical displacement of the nodes on the top surface of the model in each numerical simulation of the tectonic stress field. The boundary conditions of several mathematical models and the vertical displacements of the corresponding nodes on the top surface of the model in an embodiment of the present application are partially shown in the appendix Figure 5 and Figure 6 .

[0121] Step S24, obtain the vertical displacement (D V ) curve of the vertical fracture corresponding to each set of numerical simulations of the tectonic stress field according to the vertical displacement of the top interface of the model in step S23, and establish the boundary conditions {θ, σ H / σ V , σ h / σ H} of the numerical simulation of the tectonic stress field and the one-to-one correspondence relationship with the D V curve.

[0122] The vertical displacement (D V ) at a certain position of the vertical fracture is defined as the difference between the displacement of the nodes on the contact surface along the Z axis and the displacement of the nodes on the target surface along the Z axis. Thus, the vertical displacement (D V ) curve corresponding to each mathematical model is calculated. The abscissa of the curve is the distance of the displacement data point from the end of one side of the fracture. In order to enhance the particularity of the central position of the vertical fracture, the distance of the displacement data point at the central position of the vertical fracture from the southern end of the fracture is set to 0 m, while the distance of the displacement data point at the southernmost end of the fracture from the southern end of the fracture is set to -2500 m, and the distance of the displacement data point at the northernmost end of the fracture from the southern end of the fracture is set to 2500 m. Therefore, the abscissa range of each D V curve is -2500 to 2500 m. When the abscissa of any point on the D V curve is less than 0, the point is located on the south side of the central part of the vertical fracture, and vice versa, it is located on the north side of the central part of the vertical fracture. Each mathematical model and the corresponding numerical simulation results of the tectonic stress field correspond to a D Vcurve, so the boundary conditions {θ,σ H / σ V ,σ h / σ H} and D V One-to-one correspondence of curves.

[0123] In one embodiment of the present application, D corresponding to each mathematical model V The curve is as follows Figure 7 As shown in the figure, the horizontal axis of each curve is the distance between the displacement data point and the southern end of the fault, with a numerical range of -2500 to 2500m. The horizontal coordinate of the displacement data point located in the center of the vertical fault is 0m, while the horizontal coordinate of the displacement data point at the southernmost end of the fault is -2500m, and the horizontal coordinate of the displacement data point at the northernmost end of the fault is 2500m. Each mathematical model and the corresponding numerical simulation results of the tectonic stress field correspond to a D V curve, so the boundary conditions {θ,σ H / σ V ,σ h / σ H} and D V In one embodiment of the present application, there are 297 models, so there are 297 D V curve.

[0124] Step S25, according to D obtained in step S24 V Curve, statistics of D for each model V Maximum value (D V-MAX ) and minimum value (D V-MIN ) ratio (D V-MAX / D V-MIN ), and D V The zero position of (L(D V =0)).

[0125] According to D obtained in step S24 V Curve, first count each D V D of the curve V Maximum value (D V-MAX ) and minimum value (D V-MIN ) ratio (D V-MAX / D V-MIN ). Vertical fracture D V Usually has positive and negative values, D V The positive or negative value represents the relative rise and fall of the strata on both sides of the vertical fault, but its absolute value can represent the strength of stratum deformation. Therefore, according to the D of each curve, VThe ratio of the maximum value to the minimum value obtained from the calculation results shall be taken as the absolute value as D in this application. V Maximum value (D V-MAX ) and minimum value (D V-MIN ) ratio (D V-MAX / D V-MIN ).

[0126] Then obtain the D of each V curve, the zero position of D V (L(D V = 0)), it should be noted that each D V curve has 3 zeros in most cases, two of which are located at the southernmost end and the northernmost end of the vertical fracture respectively, that is, when the abscissa of the displacement data points takes the minimum value of -2500m and the maximum value of 2500m respectively. However, in this application, L(D V = 0) is not the above two zeros, that is, L(D V = 0) is not equal to -2500m and 2500m, but equal to a certain value between -2500 and 2500. Divide L(D V = 0) by 2500m and then normalize L(D V = 0) to a decimal between -1 and 1. When L(D V = 0) approaches -1 infinitely, the zero position of the D V curve approaches the southernmost side of the vertical fracture infinitely. When L(D V = 0) approaches 1 infinitely, the zero position of the D V curve approaches the northernmost side of the vertical fracture infinitely. When L(D V = 0) approaches 0 infinitely, the zero position of the D V curve approaches the center of the vertical fracture infinitely. When the D V curve of a certain model has only two zeros, this curve will not be included in the analysis of the subsequent steps, because in this case, there is no third zero in the D V curve, and L(D V = 0) is meaningless. It should be noted that in the case where the D V curve has only two zeros, the parameter D V-MAX / D V-MIN is also meaningless, because the denominator D V-MIN is infinitely close to 0, so the parameter D V-MAX / D V-MIN is not included in the subsequent analysis either.

[0127] In an embodiment of this application, according to the 297 D V curves obtained in step S24, first count the D of each D V curveV Maximum value (D V-MAX ) and minimum value (D V-MIN ) of the ratio (D V-MAX / D V-MIN ). The D of the vertical fracture V usually has both positive and negative values (as shown in Attachment Figure 7 ), and the positive or negative of D V represents the relative uplift and subsidence of the strata on both sides of the vertical fracture, but its absolute value can represent the formation deformation intensity. Therefore, the ratio of the maximum value and the minimum value obtained according to the calculation result of D V of each curve should be taken as the absolute value as the D V maximum value (D V-MAX ) and minimum value (D V-MIN ) of the ratio (D V-MAX / D V-MIN ) in this application.

[0128] In an embodiment of this application, each D V curve has 3 zeros in most cases (as shown in Attachment Figure 7 ), and L(D V = 0) is not the two zeros located at the ends of the vertical fracture, that is, L(D V = 0) is not equal to -2500m and 2500m. After dividing L(D V = 0) by 2500m, L(D V = 0) is normalized to a decimal between -1 and 1. When L(D V = 0) approaches -1 infinitely, the zero position of the D V curve approaches the southernmost side of the vertical fracture infinitely. When L(D V = 0) approaches 1 infinitely, the zero position of the D V curve approaches the northernmost side of the vertical fracture infinitely. When L(D V = 0) approaches 0 infinitely, the zero position of the D V curve approaches the center of the vertical fracture infinitely. When the D V curve of a certain model has only two zeros, this curve is not included in the analysis of the subsequent steps. The D V-MAX / D V-MIN and L(D V = 0) used in the embodiment of this application and their one-to-one correspondence with the numerical simulation boundary conditions of the tectonic stress field {θ, σ H / σ V , σ h / σ H} are shown in Table 2, where the data columns D V-MAX / D V-MIN and L(D VIn the case of ( = 0), the model marked by " / " and the numerical simulation results of its construction stress field are not included in the subsequent analysis because the corresponding D V curve has only two zeros.

[0129] Table 2

[0130]

[0131]

[0132]

[0133]

[0134]

[0135]

[0136]

[0137]

[0138] Step S3, construct the conversion model between L(D V = 0) and D V-MAX / D V-MIN .

[0139] Step S31, when σ H / σ V takes different values, fit the valid data of L(D V = 0) and D V-MAX / D V-MIN obtained from all models to get the conversion model 1 of L(D V = 0) and D V-MAX / D V-MIN .

[0140] According to the σ H / σ V interval and value set in step S23, when σ H / σ V takes different values, fit the valid data of L(D V = 0) and D V-MAX / D V-MIN obtained from step S25. If a total of n different σ H / σ VFor the n different numerical values, n fitting functions corresponding to these n different numerical values are obtained through fitting. The combination of these functional equations is called the conversion model 1. In the conversion model 1, all the fitting functions should belong to the same type, only the constant coefficients are different. Therefore, after constructing the conversion model 1, it is also necessary to fit the functional relationship between the constant coefficient in each fitting function and σ H / σ V If the fitting function is not clear, then according to the one-to-one correspondence between the constant coefficient and σ H / σ V A relevant curve graph is drawn as a discrimination chart. The horizontal axis of the curve graph is the value of σ H / σ V and the vertical axis is each constant coefficient. For each undetermined constant coefficient, a curve needs to be drawn.

[0141] In an embodiment of the present application, the values of σ H / σ V set in step S23 are 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, and 3. However, only when the value of σ H / σ V is 1, 2, and 3, the amount of data is relatively large, and the effective data are 62, 74, and 78 respectively. When the value of σ H / σ V takes other values, the effective data are all less than 10. Therefore, only the effective data when the value of σ H / σ V is 1, 2, and 3 are included in the analysis of this step. When σ H / σ V takes different values, the L(D V =0) and D V-MAX / D V-MIN effective data obtained from step S25 are fitted, and it is found that the function shown in formula (5) has very high applicability at this time:

[0142] Y = a - b·c X (5)

[0143] Among them, X refers to D V-MAX / D V-MIN , Y refers to L(D V =0), and a, b, and c are all constant coefficients. A total of 3 fitting function equations are obtained, and the goodness of fit is 0.96, 0.93, and 0.89 respectively, as shown in the appendix Figure 8 . The combination of these functional equations is called the conversion model 1. It is also necessary to solve the fitting functions of the constant coefficients a, b, and c in formula (5) and σ H / σ V . When σ H / σ V When it is equal to 1, 2, and 3, the constant coefficients a are 0.78, 0.75, and 0.72 respectively, the constant coefficients b are 1.14, 2.19, and 3.19 respectively, and the constant coefficients c are 0.57, 0.34, and 0.22 respectively. Therefore, the parameter σ can be used H / σ V Specify the exponential function model between parameters X and Y as the system of equations (6):

[0144]

[0145] Among them, X refers to D V-MAX / D V-MIN , and Y refers to L(D V =0). The system of equations (6) is the final expression form of transformation model 1. In this system of equations, only σ H / σ V is needed to construct the transformation model of L(D V =0) and D V-MAX / D V-MIN .

[0146] Step S32 When the boundary condition θ takes different values, fit the valid data of L(D V =0) and D V-MAX / D V-MIN obtained from all models to obtain the transformation model 2 of L(D V =0) and D V-MAX / D V-MIN .

[0147] According to the θ interval and values set in step S21, when θ takes different values, fit the valid data of L(D V =0) and D V-MAX / D V-MIN obtained in step S25. If a total of n different θ values are set, then n fitting functions corresponding to these n different values are obtained through fitting. The combination of these function equations is called the transformation model 2. In the transformation model 2, all fitting functions should belong to the same type, but only the constant coefficients are different. Therefore, after constructing the transformation model 2, it is also necessary to fit the functional relationship between the constant coefficients in each fitting function and θ. If the fitting function is not clear, then according to the one-to-one correspondence between the constant coefficients and σ H / σ V , draw the relevant curve graph as the discrimination chart. The horizontal axis of the curve graph is the value of σ H / σ V , and the vertical axis is each constant coefficient. For each undetermined constant coefficient, a curve needs to be drawn.

[0148] In one embodiment of the present application, the values of θ set in step S21 are 10°, 20°, 30°, 40°, 45°, 50°, 60°, 70° and 80°. When θ takes different values, the L(D V =0) and D V-MAX / D V-MIN valid data obtained from step S25 are fitted, and it is found that equation (5) can still be used to construct the fitting function. The results are as shown in Appendix Figure 9 . The goodness of fit of the 9 fitting functions is at least 0.88, indicating the applicability of equation (5) here. When fitting the relationship between the coefficients of each function and θ, it is found that the relationship between the constant coefficients a, b, and c in Appendix Figure 9 and θ cannot be well fitted by a conventional function. Therefore, the curves of the constant coefficients a, b, and c with respect to θ are plotted as shown in Appendix Figure 10 . For θ between 10° and 80°, the linear interpolation method is used to obtain the corresponding constant terms a, b, and c. For θ less than 10° and greater than 80°, the linear trend extrapolation method is used to obtain the corresponding constant terms a, b, and c. The conversion model 2 of L(D V =0) and D V-MAX / D V-MIN is specified as the system of equations (7) using the parameter θ:

[0149]

[0150] where X refers to D V-MAX / D V-MIN , and Y refers to L(D V =0). The curves of the constant terms a, b, and c with respect to the parameter θ are as shown in Appendix Figure 10 .

[0151] Step S33, when the boundary condition σ h / σ H takes different values, the L(D V =0) and D V-MAX / D V-MIN valid data obtained from all models are fitted to obtain the conversion model 3 of L(D V =0) and D V-MAX / D V-MIN .

[0152] According to the σ h / σ H interval and values set in step S23, when σ h / σ H takes different values, the L(D V =0) and D V-MAX / D V-MINFit the valid data. If n different values of σ are set in total h / σ H values, then n fitting functions corresponding to these n different values are obtained through fitting. The combination of these function equations is called the conversion model 3. In the conversion model 3, all fitting functions should belong to the same type, only the constant coefficients are different. Therefore, after constructing the conversion model 3, it is also necessary to fit the functional relationship between the constant coefficient in each fitting function and σ h / σ H . If the fitting function is not clear, then according to the one-to-one correspondence between the constant coefficient and σ h / σ H , draw the relevant curve graph as the discrimination chart. The horizontal axis of the curve graph is the value of σ h / σ H , and the vertical axis is each constant coefficient. For each undetermined constant coefficient, a curve needs to be drawn.

[0153] In an embodiment of the present application, the value of σ h / σ H set in step S23 is 0.1 to 0.9 at an interval of 0.1. When σ h / σ H takes different values, fit the L(D V = 0) and D V-MAX / D V-MIN valid data obtained from step S25. Equation (5) can still be used to construct the fitting function, but its applicability is lower than that in steps S21 and S22, and the result is as shown in the appendix Figure 11 . When σ h / σ H is equal to 0.5, the fitting function cannot be constructed. This may be because the data points corresponding to the higher value of the parameter L(D V = 0) are fewer and do not play a good constraining effect on the extension of the exponential function in the positive direction of the L(D V = 0) axis. When σ h / σ H is equal to 0.9, the data group (L(D V = 0), D V-MAX / D V-MIN ) is distributed irregularly, so the exponential function cannot be constructed. And when σ h / σ H is increased from 0.7 to 0.8, the goodness of fit when fitting with equation (5) has rapidly decreased from 0.93 to 0.67, which shows that one of the applicable conditions for fitting the data group (L(D V = 0), D V-MAX / D V-MIN ) with equation (5) is that σh / σ H Less than or equal to 0.7 (the goodness of fit when fitting with Equation (5) for the corresponding cases is greater than or equal to 0.92). Therefore, a total of 7 effective fitting functions were obtained in this step, which are σ h / σ H when taking 0.1, 0.2, 0.3, 0.4, 0.6, 0.7, and 0.8.

[0154] When fitting the coefficients of each function with σ h / σ H it was found that the Figure 9 constant coefficients a, b, and c in h / σ H could not be well fitted with a conventional function. Therefore, the curves of the constant coefficients a, b, and c with respect to σ h / σ H are shown in Appendix Figure 12 For σ h / σ H ranging from 0.1 to 0.9, the corresponding constant terms a, b, and c are obtained by linear interpolation. For σ h / σ H less than 0.1 and greater than 0.9, the corresponding constant terms a, b, and c are obtained by linear trend extrapolation. Using the parameter σ h / σ H to specify the conversion model 3 of L(D V = 0) and D V-MAX / D V-MIN into the system of equations (8):

[0155]

[0156] where X refers to D V-MAX / D V-MIN , and Y refers to L(D V = 0). The curves of the constant terms a, b, and c with respect to the parameter σ h / σ H are shown in Appendix Figure 12 .

[0157] Step S4, obtain the D V curve of the vertical fracture that was the research target in step S1, as well as the corresponding L(D V = 0) and D V-MAX / D V-MIN .

[0158] Step S41, obtain the D V curve of the vertical fracture that was the research object in step S1.

[0159] According to the method of calculating the vertical fracture D described in step S24, obtain the D curve of the vertical fracture under study V The horizontal axis of this curve is the distance from the displacement data point to the southernmost side of the fracture, the minimum value of the abscissa is 0, and the maximum value is the total length of the fracture V In an embodiment of the present application, the D curve of the vertical fracture F1 under study is as shown in Appendix

[0160] A. The horizontal axis of the curve is the distance from the displacement data point to the southernmost side of the fracture F1, the minimum value of the abscissa is 0, and the maximum value is the total length of the fracture F1, which is 6.834 km V Appendix Figure 13 A

[0161] Step S42: Calculate L(D = 0) and D / D according to the D curve of the vertical fracture under study V curve V =0) and D V-MAX / D V-MIN .

[0162] According to the D curve of the vertical fracture under study obtained from step S41, and then solve the corresponding L(D V =0) and D V =0) and D V-MAX / D V-MIN .

[0163] In an embodiment of the present application, the D curve of the vertical fracture F1 obtained from step S41 is as shown by the red line in Appendix V Appendix Figure 13 A. According to the method mentioned in step S25, the corresponding L(D V =0) is equal to -0.736, and the obtained D V-MAX / D V-MIN is equal to 0.308, which are all marked in Appendix Figure 13 A

[0164] Step S5: According to L(D V =0) and D V-MAX / D V-MIN of the vertical fracture under study obtained from step S42, and the three conversion models proposed in step S3, obtain the corresponding σ H / σ V 、σ h / σ H and θ, and accordingly reconstruct the regional paleotectonic stress field related to the strike-slip activity of the fracture F1

[0165] Step S51: Respectively use σ H / σ V 、σh / σ H With σ / σ and θ as the horizontal axes and L(D V = 0) as the vertical axis, three function curves are plotted.

[0166] First, with σ / σ H / σ V as the horizontal axis and L(D V = 0) as the vertical axis, a discrimination curve is plotted. Equation (6) gives the conversion model 1 of L(D H / σ V = 0) and D / D V = 0) when σ / σ V-MAX / D V-MIN takes any value. Substituting the D / D V-MAX / D V-MIN of the vertical fracture obtained as the research object from step S42 into the variable X of equation (6), the functional relationship of L(D V = 0) with respect to σ / σ H / σ V can be obtained. According to this functional relationship, the corresponding function curve can be plotted. Similarly, the function curves of L(D V = 0) with respect to σ / σ h / σ H and θ can be plotted successively according to equations (7) and (8).

[0167] In an embodiment of the present application, when D / D V-MAX / D V-MIN is equal to 0.308, it is respectively substituted into equations (6) to (8), and the functional relationships of L(D V = 0) with respect to σ / σ H / σ V , σ / σ h / σ H and θ can be obtained, and the corresponding function curves can be plotted accordingly, as shown by the thick black, red, and blue solid lines in Appendix Figure 13 B respectively.

[0168] Step S52: According to the L(D V = 0) of the vertical fracture obtained as the research object from step S42, a horizontal line with a constant ordinate equal to this L(D V = 0) is plotted on the function curve graph obtained from step S51. According to the intersection points of this horizontal line and the function curve graph obtained from step S51, the regional paleotectonic stress field related to the strike activity of fracture F1 is reconstructed.

[0169] According to the L(D V = 0) of the vertical fracture obtained as the research object from step S42, a horizontal line with a constant ordinate equal to this L(DV a horizontal line L where (D = 0), and perform the following operations:

[0170] Operation 1: If the horizontal line L intersects the function curve of L(D V = 0) with respect to σ H / σ V , then consider the abscissa of this intersection point as the σ H / σ V value, which represents the σ H / σ V value related to the strike - slip of the vertical fracture under study. If the horizontal line L does not intersect the function curve of L(D V = 0) with respect to σ H / σ V , then find the point on the function curve of L(D V = 0) with respect to σ H / σ V that is closest to the horizontal line L, and identify this as intersection point 1 here.

[0171] Operation 2: If the horizontal line L intersects the function curve of L(D V = 0) with respect to σ h / σ H , then consider the abscissa of this intersection point as the σ h / σ H value, which represents the σ h / σ H value related to the strike - slip of the vertical fracture under study. If the horizontal line L does not intersect the function curve of L(D V = 0) with respect to σ h / σ H , then find the point on the function curve of L(D V = 0) with respect to σ h / σ H that is closest to the horizontal line L, and identify this as intersection point 2 here.

[0172] Operation 3: If the horizontal line L intersects the function curve of L(D V = 0) with respect to θ, then consider the abscissa of this intersection point as the θ value, which represents the θ value related to the strike - slip of the vertical fracture under study. If the horizontal line L does not intersect the function curve of L(D V = 0) with respect to θ, then find the point on the function curve of L(D V = 0) with respect to θ that is closest to the horizontal line L, and identify this as intersection point 3 here.

[0173] The abscissas of the above intersection points 1, 2, and 3 respectively represent the σ of the corresponding regional tectonic stress field when the vertical fault under study undergoes strike-slip along its strike. H / σ V 、σ h / σ H and θ, where σ H / σ V and σ h / σ H can reflect the relative magnitudes of the regional horizontal maximum principal stress, the regional horizontal minimum principal stress, and the vertical stress, and can be used to reflect the type of tectonic stress (Anderson, E.M., 1951. The Dynamics of Faulting and Dyke Formation with Application to Britain, The Journal of Geology, 51(2), 140 - 140). The parameter θ reflects the angle between the regional horizontal maximum principal stress and the vertical fault. Combining with the direction of the strike-slip of the vertical fault, the precise azimuth of the regional horizontal maximum principal stress can be quantitatively determined. Thus, the reconstruction of the paleotectonic stress field has been completed.

[0174] In an embodiment of the present application, from step S42, it is solved that L(D V =0) of the vertical fault F1 is equal to -0.736. A horizontal line L with a vertical coordinate equal to -0.736 is drawn in Appendix Figure 13 B, and the following operations are performed:

[0175] Operation 1: The horizontal line L and the function curve of L(D Figure 13 =0) in Appendix V B with respect to σ H / σ V have 1 intersection point. The abscissa of this intersection point is 1.863. Then, it is considered that σ H / σ V related to the strike-slip of the vertical fault F1 is 1.863.

[0176] Operation 2: The horizontal line L and the function curve of L(D Figure 13 =0) in Appendix V B with respect to σ h / σ H have 1 intersection point. The abscissa of this intersection point is 0.065. Then, it is considered that σ h / σ H related to the strike-slip of the vertical fault F1 is 0.065.

[0177] Operation 3: The horizontal line L and the L(D Figure 13 =0) in AppendixV When the function curve of = 0 with respect to θ has two intersection points, and the abscissas are 27.14° and 35.98° respectively, it is considered that θ related to the strike-slip of the vertical fault F1 is equal to 27.14° or 35.98°.

[0178] σ related to the strike-slip of the vertical fault F1 H / σ V and σ h / σ H are equal to 1.863 and 0.065 respectively, corresponding to a type III stress field (σ H > σ V > σ h ). θ is equal to 27.14° or 35.98°, which indicates that the included angle between the regional horizontal maximum principal stress and the vertical fault is 27.14° or 35.98°. It can be seen from step S1 that the vertical fault F1 undergoes a left-lateral strike-slip along its strike. Therefore, after the regional horizontal maximum principal stress direction rotates clockwise by 27.14° or 35.98°, it is parallel to the trace of the vertical fault F1 on the plane. Attached Figure 13 Figure C shows the plane trace of the vertical fault F1, and its strike is NW - 56.68°. So the regional horizontal maximum principal stress direction is NWW - 83.82° or SWW - 92.66°, as shown in attached Figure 13 Figure C. Thus, the reconstruction of the paleotectonic stress field has been completed.

[0179] The technical features of the above embodiments can be combined arbitrarily. For the sake of concise description, not all possible combinations of the technical features in the above embodiments are described. However, as long as the combinations of these technical features do not conflict, they should all be considered as the scope recorded in this specification.

Claims

1. A method for reconstructing paleo-tectonic stress field based on the vertical displacement of vertical faults, characterized in that: include: Step S1, determining the study area, target layer and vertical fault as the study object, and obtaining the core, conventional logging and array acoustic logging data of the target layer to solve the static rock mechanics parameters of the target layer; Step S2, constructing geomechanical models and mathematical models related to the sliding of the vertical fault along the strike, and conducting numerical simulation tests of the tectonic stress field, obtaining the vertical deformation characteristics of the strata on both sides of the vertical fault from the test results, and drawing the vertical displacement D of the vertical fault corresponding to different models V curve, and calculate the D of each model V Maximum value D V-MAX and minimum value D V-MIN The ratio D V-MAX / D V-MIN , and D V The zero position of L(D V =0); Step S3, construct L(D V =0) and D V-MAX / D V-MIN The conversion model between Step S4, according to the D V The curve calculates L(D V =0) and D V-MAX / D V-MIN ; Step S5, according to the vertical fracture L(D V =0) and D V-MAX / D V-MIN , and the conversion model proposed in step S3, obtain the corresponding σ H / σ V , σ h / σ H and θ, and based on this, the regional paleotectonic stress field associated with the along-strike activity of the vertical faults is reconstructed.

2. The method for reconstructing the paleo-tectonic stress field based on the vertical displacement of vertical faults according to claim 1, characterized in that: The step S1 comprises: For the study area with core data but without array acoustic logging data, rock mechanics tests on core samples were carried out to obtain static rock mechanics parameters; each static rock mechanics parameter was taken as the average value of the corresponding parameters of all samples; the static rock mechanics parameters included static Young's modulus, static Poisson's ratio and density; For the study area without core data, the static rock mechanical parameters of the corresponding lithology in the nearby area were obtained by consulting the literature on the basis of identifying the main lithology of the target strata; For the study area with core, conventional logging and array acoustic logging data, the dynamic rock mechanical parameters are first calculated based on conventional logging and array acoustic logging. The formula is: Among them, E d Represents dynamic Young's modulus, in GPa; μ d represents the dynamic Poisson's ratio; Δt p Indicates the longitudinal wave time difference DTC, the unit is μs / ft; Δt s represents the shear wave time difference DTS, in μs / ft; ρ represents the density, in g / cm 3 ; The dynamic rock mechanics parameters include dynamic Young's modulus and dynamic Poisson's ratio; Then, according to the static rock mechanics parameters obtained from the rock mechanics test of the core sample, the dynamic and static rock mechanics parameters are fitted to obtain the conversion model of the dynamic and static rock mechanics parameters of the target layer. According to the conversion model, the dynamic rock mechanics parameters are converted into static rock mechanics parameters. The conversion formula is: In the formula, E d Represents dynamic Young's modulus, in GPa; μ d represents dynamic Poisson's ratio; E represents static Young's modulus, in GPa; μ represents static Poisson's ratio.

3. The method for reconstructing the ancient tectonic stress field based on the vertical displacement of vertical faults according to claim 1 is characterized in that: The step S2 comprises: According to the geological characteristics of the study area, a three-dimensional geomechanical model is constructed, which includes the location, strike, extension length of the vertical fault and the geometric characteristics of the target strata; the vertical fault is 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 S1 to each unit in the mathematical model; set the boundary conditions and constraint conditions of the mathematical model according to the requirements of research accuracy; carry out numerical simulation tests of tectonic stress fields to obtain simulation results corresponding to each test; According to the vertical displacement of the top interface of the model, the vertical displacement D of the vertical fault corresponding to each set of tectonic stress field numerical simulation experiments is obtained. V Curve, establish the boundary conditions and vertical displacement curve D for numerical simulation of tectonic stress field V One-to-one correspondence between ; According to the vertical displacement D V Curve, calculate the vertical displacement D of each model V The maximum value D V-MAX and minimum value D V-MIN The ratio D V-MAX / D V-MIN , and the vertical displacement D V The zero position of L(D V =0).

4. The method for reconstructing the ancient tectonic stress field based on the vertical displacement of vertical faults according to claim 1 is characterized in that: The step S3 comprises: When the ratio of the regional maximum horizontal principal stress to the vertical stress is σ H / σ V When different values ​​are taken, the zero position L(D V =0) and the ratio of the maximum value to the minimum value D V-MAX / D V-MIN The effective data is fitted to obtain L(D V =0) and D V-MAX / D V-MIN The first conversion model of When the vertical fracture and the regional maximum horizontal principal stress σ H When the angle θ of the direction takes different values, the L(D V =0) and D V-MAX / D V-MIN The effective data is fitted to obtain L(D V =0) and D V-MAX / D V-MIN The second conversion model of When the ratio of the regional minimum horizontal principal stress to the regional maximum horizontal principal stress is σ h / σ H When different values ​​are taken, the L(D V =0) and D V-MAX / D V-MIN The effective data is fitted to obtain L(D V =0) and D V-MAX / D V-MIN The third conversion model.

5. The method for reconstructing the ancient tectonic stress field based on the vertical displacement of vertical faults according to claim 4 is characterized in that: Building the first, second, and third transition models also includes: When σ H / σ V When different values ​​are taken, L(D V =0) and D V-MAX / D V-MIN If multiple different σ are set, H / σ V numerical value, then multiple corresponding fitting functions are obtained by fitting and recorded as the first conversion model; In the first transformation model, all fitting functions are of the same type, but their constant coefficients are different; the constant coefficient in each fitting function is related to σ H / σ V If the fitting function is not clear, then according to the constant coefficient and σ H / σ V , a first correlation curve diagram is drawn as a discriminant diagram, wherein the horizontal axis of the first correlation curve diagram is σ H / σ V , the vertical axis is the constant coefficients of each item, and a curve is drawn for each undetermined constant coefficient; When θ takes different values, L(D V =0) and D V-MAX / D V-MIN The valid data is fitted. If multiple different θ values ​​are set, multiple corresponding fitting functions are obtained through fitting and recorded as the second conversion model; In the second conversion model, all fitting functions belong to the same type, but their constant coefficients are different; the functional relationship between the constant coefficient and θ in each fitting function is fitted. If the fitting function is unclear, a second correlation curve diagram is drawn as a discriminant diagram according to the one-to-one correspondence between the constant coefficient and θ, wherein the horizontal axis of the second correlation curve diagram is θ and the vertical axis is the constant coefficients of each item; a curve is drawn for each undetermined constant coefficient; When σ h / σ H When different values ​​are taken, L(D V =0) and D V-MAX / D V-MIN If multiple different σ are set, h / σ H If the value is , then multiple corresponding fitting functions are obtained by fitting and recorded as the third conversion model; In the third transformation model, all fitting functions belong to the same type, but their constant coefficients are different. The constant coefficients in each fitting function are related to σ h / σ H If the fitting function is not clear, then according to the constant coefficient and σ h / σ H , a third correlation curve diagram is drawn as a discriminant diagram, wherein the horizontal axis of the third correlation curve diagram is σ h / σ H , the vertical axis is the constant coefficients, and a curve is drawn for each undetermined constant coefficient.

6. The method for reconstructing the ancient tectonic stress field based on the vertical displacement of vertical faults according to claim 5, characterized in that: Building the first, second, and third transition models also includes: When σ H / σ V When different values ​​are taken, L(D V =0) and D V-MAX / D V-MIN The effective data is fitted, and the fitting formula is: Y=a-b·c X Where X represents D V-MAX / D V-MIN , Y represents L(D V =0), a, b and c all represent constant coefficients; using the parameter σ H / σ V The exponential function model between parameters X and Y is concretized as the first conversion model: In the formula, X represents D V-MAX / D V-MIN , Y represents (D V =0); a1, b1 and c1 represent constant coefficients of the first conversion model; When θ takes different values, the fitting formula is used to calculate L(D V =0) and D V-MAX / D V-MIN Fitting the valid data; The fitting formula is concretized into the second conversion model using parameter θ: Where X represents D V-MAX / D V-MIN , Y represents L(D V =0); a2, b2 and c2 represent constant coefficients of the second conversion model; When σ h / σ H When different values ​​are taken, the fitting formula is used to calculate L(D V =0) and D V-MAX / D V-MIN Fitting the valid data; Using the parameter σ h / σ H The fitting formula is concretized into the third conversion model: Where X represents D V-MAX / D V-MIN , Y represents L(D V =0); a3, b3 and c3 represent constant coefficients of the third conversion model.

7. The method for reconstructing paleo-tectonic stress field based on vertical displacement of vertical faults according to claim 1, characterized in that: The step S5 comprises: σ H / σ V , θ and σ h / σ H As the horizontal axis, L(D V =0) as the vertical axis to draw the corresponding discriminant curve; According to L(D V =0), draw a vertical coordinate on the corresponding discriminant curve that is always equal to L(D V =0), and based on the intersection of the horizontal line and the discriminant curve, the regional paleo-tectonic stress field related to the along-strike activity of the vertical fault is reconstructed.

8. The method for reconstructing the ancient tectonic stress field based on the vertical displacement of vertical faults according to claim 7, characterized in that: The step S5 further comprises: σ H / σ V As the horizontal axis, L(D V =0) as the vertical axis to draw the first discriminant curve; if the horizontal line L(D V =0) and the first discriminant curve graph, the intersection is recorded as the first intersection; if the horizontal line L(D V =0) has no intersection with the first discriminant curve, then in L(D V =0) about σ H / σ V Find the horizontal line L(D V =0) the point closest to the first intersection point; With θ as the horizontal axis and L(D V =0) as the vertical axis to draw the second discriminant curve; if the horizontal line L(D V =0) and the second discriminant curve graph have an intersection, which is recorded as the second intersection; if the horizontal line L(D V =0) has no intersection with the second discriminant curve graph, then in L(D V =0) on the function curve of θ and the horizontal line L(D V =0) is the point closest to the second intersection point; σ h / σ H As the horizontal axis, L(D V =0) as the vertical axis to draw the third discriminant curve; if the horizontal line L(D V =0) and the third discriminant curve graph, it is recorded as the third intersection point; if the horizontal line L(D V =0) has no intersection with the third discriminant curve graph, then in L(D V =0) about σ h / σ H Find the point on the function curve graph that is closest to the horizontal line L and record it as the third intersection point; The abscissas of the first, second and third intersections represent the σ of the corresponding regional tectonic stress field when the vertical fault as the research object slides along its strike. H / σ V , θ and σ h / σ H , where σ H / σ V and σ h / σ H It reflects the relative sizes of the regional maximum horizontal principal stress, the regional minimum horizontal principal stress and the vertical stress, and reflects the type of tectonic stress; θ reflects the angle between the regional maximum horizontal principal stress and the vertical fault. Combined with the direction of sliding along the vertical fault, the precise azimuth of the regional maximum horizontal principal stress can be quantitatively determined. At this point, the reconstruction of the paleo-tectonic stress field has been completed.

Images