Automatic modeling method for numerical simulation of strata movement based on three-dimensional model

Abstract

The application discloses a three-dimensional model-based numerical simulation automatic modeling method for stratum movement, and particularly relates to the technical field of three-dimensional geological modeling and numerical simulation of stratum movement, and is used for solving the problem that the existing automatic modeling of smooth interpolation smoothes the micro-undulation of a sedimentary weak surface into an ideal curved surface, so that the three-dimensional model cannot reproduce the interlayer separation and stepwise subsidence non-continuous response in the stratum movement process; the method is characterized in that: a smooth reference curved surface is generated by smooth interpolation of stratum demarcation points, a deviation scalar field is formed by calculating the normal deviation, the local autocorrelation length is calculated along the borehole path, the sedimentary weak surface influence domain is determined in combination with the lithology change interface, the positive deviation local uplift cluster is identified in the influence domain, the undulation dominant direction and the roughness amplitude are extracted, the roughness amplitude is constructed along the undulation dominant direction to form an undulation term superimposed to the smooth reference curved surface to generate a non-continuous layer surface retaining the geometric characteristics of the sedimentary weak surface, and the three-dimensional geological model is constructed by taking the non-continuous layer surface as the stratum interface.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional geological modeling and numerical simulation of rock strata movement, and more specifically, to an automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model. Background Technology

[0002] In the field of numerical simulation of mining subsidence and rock strata movement, constructing a three-dimensional model that reflects the real stratigraphic structure is a prerequisite for reliable calculations. Existing automated modeling processes typically begin by extracting the boundary points of each rock stratum from exploration data such as boreholes and profiles. Then, spatial interpolation methods are used to generate continuous three-dimensional layers. Through solid modeling and mesh generation, a numerical model that can be directly used for finite element or discrete element calculations is formed. This approach achieves automatic conversion from discrete point data to a three-dimensional geometric model, significantly reducing manual intervention, and falls entirely within the scope of three-dimensional geological modeling methods.

[0003] However, when the above process automatically interpolates and reconstructs the rock strata, the continuous smoothing interpolation process systematically smooths the inherent discontinuous geometric features of the rock strata, such as weak sedimentary surfaces and micro-undulations, into an ideal curved surface. This results in the absence of structural roughness in the constructed 3D model, which plays a controlling role in key discontinuous movement behaviors such as rock strata separation along weak surfaces and bedding slip. Consequently, the subsequent numerical simulation calculations tend to favor overall continuous deformation and cannot reproduce the discontinuous responses such as interlayer delamination and step settlement that actually occur during the movement of the rock strata. This leads to a fundamental deviation between the predicted results and the actual deformation patterns on site. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides an automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] An automated modeling method for numerical simulation of rock strata movement based on a 3D model includes the following steps:

[0007] S1. Perform smooth interpolation on the rock strata boundary points extracted from boreholes and profiles to generate smooth reference surfaces for each rock stratum.

[0008] S2. Calculate the normal deviation of each rock layer boundary point relative to the smooth reference surface to form a deviation scalar field;

[0009] S3. Calculate the local autocorrelation length of the deviation scalar field along the borehole path. Identify the adjacent boundary points at the lithological change interface where the local autocorrelation length exceeds the statistical base value and the deviation amplitude is positively correlated as the influence domain of the sedimentary weak surface.

[0010] S4. Within the influence domain of the weak depositional surface, identify the positive deviation local uplift clusters on the curved surface constructed by the normal deviation, take the extension direction of each uplift cluster as the undulation dominance direction by area weighted average, and take the maximum height difference of each uplift cluster as the roughness amplitude by area weighted average.

[0011] S5. The roughness amplitude is constructed as an undulation term along the direction of the undulation dominance, and superimposed on the smooth reference surface to generate a discontinuous layer that retains the geometric characteristics of the weak sedimentary surface. The discontinuous layer is used as the stratigraphic interface to construct a three-dimensional geological model for numerical simulation of rock strata movement.

[0012] Furthermore, smooth interpolation is performed on the rock strata boundary points extracted from boreholes and profiles to generate smooth reference surfaces for each rock stratum, including:

[0013] Outlier detection is performed on the rock strata boundary points to remove abnormal boundary points that deviate from the statistical distribution; based on the rock strata attitude revealed by the borehole, the anisotropy parameters of the strata are calculated, and Kriging interpolation that takes into account these anisotropy parameters is used to perform spatial interpolation on the removed rock strata boundary points to generate a smooth reference surface for each rock stratum.

[0014] Furthermore, the anisotropy parameters for calculation include: extracting strike azimuth and dip angles from the rock strata exposed by the borehole; determining the major axis direction of the interpolation ellipsoid using the strike azimuth; determining the stretching ratio between the major and minor axes using the tangent of the dip angle; and constructing anisotropic distance metrics to replace the Euclidean distance in Kriging interpolation for spatial interpolation.

[0015] Furthermore, the normal deviation of each rock stratum boundary point relative to the smooth reference surface is calculated to form a deviation scalar field, including:

[0016] Project the projection lines along the normal of each rock layer boundary point on the smooth reference surface, calculate the directed distance from the rock layer boundary point to the smooth reference surface, and use this directed distance as the normal deviation; perform spatial continuous interpolation on the normal deviation of each rock layer boundary point to form a deviation scalar field covering the range of the rock layer.

[0017] Furthermore, the directed distance from the rock stratum boundary point to the smooth reference surface is obtained by: taking the rock stratum boundary point as the starting point, drawing a ray along the local normal vector of the smooth reference surface, calculating the intersection point of the ray and the smooth reference surface, taking the projection length of the vector from the rock stratum boundary point to the intersection point on the normal vector as the directed distance, and taking the normal of the nearest triangular element at the discontinuity of the intersection point on the smooth reference surface.

[0018] Furthermore, the local autocorrelation length of the deviation scalar field is calculated along the borehole path. Adjacent boundary points at lithological change interfaces, whose local autocorrelation lengths exceed the statistical baseline and whose deviation amplitudes are positively correlated, are identified as the influence domains of sedimentary weak surfaces, including:

[0019] The values ​​of the deviation scalar field at each stratum boundary point are extracted along the borehole path to form a deviation sequence. The local autocorrelation length of each stratum boundary point is obtained by performing sliding autocorrelation on the deviation sequence. The location of the lithological change interface is obtained from the borehole lithology logging. The local autocorrelation lengths of all stratum boundary points are counted, and the sum of their mean and standard deviation is used as the statistical base value. The boundary points at the lithological change interface where the local autocorrelation length exceeds the statistical base value are used as seed points. The set of seed points and their spatially adjacent boundary points whose deviation amplitude is positively correlated with the seed points are determined as the influence domain of the sedimentary weak surface.

[0020] Furthermore, when calculating the local autocorrelation length of all strata boundaries, the standard deviation of the local autocorrelation length within each stratum is calculated using lithology as the stratification variable. The maximum value of the sum of the mean and the standard deviation within each stratum is used as the statistical base value. During the sliding autocorrelation calculation, the sliding window length is adaptively adjusted according to the thickness of the strata exposed by the borehole, and the window length at the interface of adjacent sedimentary cycles is taken as the smaller of the two values.

[0021] Furthermore, within the influence domain of the weak sedimentary surface, positive deviation local uplift clusters are identified on the surface constructed by the normal deviation. The extension direction of each uplift cluster is weighted by area as the dominant undulation direction, and the maximum height difference of each uplift cluster is weighted by area as the roughness amplitude, including:

[0022] Within the influence domain of the weak sedimentary surface, a deviation surface is constructed using the normal deviation as the elevation. Watershed segmentation is performed on the positive deviation region to extract spatially continuous local uplift clusters. For each local uplift cluster, a minimum bounding rectangle is fitted, and the major axis direction of the minimum bounding rectangle is taken as the extension direction of the local uplift cluster. The difference between the maximum and minimum normal deviations within the local uplift cluster is taken as the maximum elevation difference. Using the projected area of ​​each local uplift cluster as the weight, the extension directions of each local uplift cluster are weighted and averaged to obtain the undulation dominance direction. The maximum elevation difference of each local uplift cluster is weighted and averaged to obtain the roughness amplitude.

[0023] Furthermore, when using the projected area as the weight, the projected area is taken as the orthogonal projected area of ​​the local uplift cluster on the tangent plane of the smooth reference surface within the influence domain of the weak sedimentary surface; the normalized height difference is obtained by dividing the maximum height difference of each local uplift cluster by the square root of the corresponding projected area, and then the normalized height difference is weighted and averaged using the projected area as the weight. The product of the weighted average result and the square root of the total projected area of ​​the influence domain of the weak sedimentary surface is used as the roughness amplitude.

[0024] Furthermore, the roughness amplitude is constructed as an undulation term along the direction of undulation dominance and superimposed onto a smooth reference surface to generate a discontinuous layer that preserves the geometric characteristics of the weak sedimentary surface. This discontinuous layer is then used as the stratigraphic interface to construct a three-dimensional geological model for numerical simulation of strata movement, including:

[0025] Using discrete points on a smooth reference surface as carriers, an offset is added at each discrete point along the direction of undulation dominance. The magnitude of the offset is determined by the product of the roughness magnitude and the projection component of the undulation dominance direction in the tangent plane at that point, forming an undulation term. The undulation term is superimposed with the coordinates of the smooth reference surface to generate a discontinuous layer with local discontinuous undulation morphology. The discontinuous layer is used as a stratigraphic interface constraint entity to model a three-dimensional geological model containing the geometric features of rock layer boundaries.

[0026] Compared with the prior art, the present invention has the following beneficial effects:

[0027] 1. This method constructs a deviation scalar field on a smooth reference surface and extracts deviation sequences along the borehole path. By coupling the local autocorrelation length with the lithological change interface, the influence domain of sedimentary weak surfaces is separated from the deviation field. Instead of directly introducing discontinuity constraints in the interpolation stage, the method uses the anomaly of the autocorrelation structure of the deviation field to characterize the location of sedimentary weak surfaces, thus decoupling weak surface identification from smooth interpolation. This maintains the robustness of automated interpolation and captures weak surface signals masked by smoothing. Within the influence domain, the undulation dominance direction and roughness amplitude are extracted using the morphological parameters of positive deviation local uplift clusters. The discontinuous geometric features of the weak surfaces are quantified from discrete deviations into two control parameters: directionality and amplitude, providing directly superimposed undulation terms for subsequent layer reconstruction.

[0028] 2. The roughness amplitude is constructed as an undulation term along the direction of the undulation dominance and superimposed on the smooth reference surface. The generated discontinuous layer retains the micro-undulation and local uplift morphology within the influence domain of the weak sedimentary surface. This makes the stratigraphic interface of the three-dimensional geological model no longer an ideal smooth surface at the weak surface. The resulting three-dimensional geological model incorporates the structural roughness that controls the separation of rock strata along the weak surface and the bedding slip. When the model is used for numerical simulation, it can produce discontinuous deformation response at the weak sedimentary surface, reducing the overall continuous deformation deviation caused by the excessive smoothness of the bedding surface. This makes the simulation results more consistent with the interlayer delamination and step-like settlement characteristics actually presented in the rock strata movement. Attached Figure Description

[0029] Figure 1 This is a flowchart of the automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model, as described in this invention.

[0030] Figure 2 This is a flowchart illustrating the determination of the influence domain of the weak deposition surface in this invention. Detailed Implementation

[0031] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0032] Example: Figure 1 This invention presents an automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model, which includes the following steps:

[0033] S1. Perform smooth interpolation on the rock strata boundary points extracted from boreholes and profiles to generate smooth reference surfaces for each rock stratum.

[0034] S2. Calculate the normal deviation of each rock layer boundary point relative to the smooth reference surface to form a deviation scalar field;

[0035] S3. Calculate the local autocorrelation length of the deviation scalar field along the borehole path. Identify the adjacent boundary points at the lithological change interface where the local autocorrelation length exceeds the statistical base value and the deviation amplitude is positively correlated as the influence domain of the sedimentary weak surface.

[0036] S4. Within the influence domain of the weak depositional surface, identify the positive deviation local uplift clusters on the curved surface constructed by the normal deviation, take the extension direction of each uplift cluster as the undulation dominance direction by area weighted average, and take the maximum height difference of each uplift cluster as the roughness amplitude by area weighted average.

[0037] S5. The roughness amplitude is constructed as an undulation term along the direction of the undulation dominance, and superimposed on the smooth reference surface to generate a discontinuous layer that retains the geometric characteristics of the weak sedimentary surface. The discontinuous layer is used as the stratigraphic interface to construct a three-dimensional geological model for numerical simulation of rock strata movement.

[0038] When implementing S1, outlier detection is performed on the strata boundary points extracted from boreholes and profiles to remove abnormal boundary points that deviate from the statistical distribution. Outlier detection uses a single stratum bedding plane as the processing unit, grouping strata boundary points belonging to the same stratum bedding plane into a point set. The median of the spatial coordinates of each stratum boundary point in the point set is calculated in the strike, dip, and vertical directions. The Mahalanobis distance from each stratum boundary point to the median is then calculated using the median as the benchmark. The covariance matrix for Mahalanobis distance calculation is composed of the sample covariance of the coordinates of the stratum boundary points within the point set. Stratum boundary points whose Mahalanobis distance exceeds the chi-square distribution quantile are marked as outliers and removed. The degrees of freedom of the chi-square distribution quantile are taken from the degrees of freedom of the three-dimensional spatial coordinates. The quantile threshold is set according to the total number of rock stratum boundary points and the expected misclassification rate. The expected misclassification rate is the reciprocal of the total number of rock stratum boundary points. The quantile threshold is the chi-square distribution quantile corresponding to the upper probability of the chi-square distribution being equal to the expected misclassification rate. When the Mahalanobis distance of a rock stratum boundary point is greater than the quantile threshold, the rock stratum boundary point is judged as an abnormal boundary point.

[0039] Anisotropy parameters of the bedding plane are calculated based on the attitude of the rock strata revealed by boreholes. The strike azimuth and dip angles are extracted from the borehole-revealed rock strata attitude. The rock strata attitude is obtained from measured records of the core axis angle and dip direction during core sampling, or from fitting bedding plane traces interpreted from borehole imaging. The major axis direction of the interpolated ellipsoid is determined by the strike azimuth, which extends along the strike azimuth in the horizontal plane. The stretching ratio between the major and minor axes is determined by the tangent of the dip angle. The stretching ratio is the ratio of the radius of the major axis to the radius of the minor axis of the interpolated ellipsoid; a larger tangent value indicates a larger stretching ratio, suggesting stronger spatial continuity of the bedding plane in the strike direction than in the dip direction. The minor axis direction is orthogonal to the major axis direction and lies within the dip plane of the bedding plane; the radius of the minor axis is determined by dividing the radius of the major axis by the stretching ratio.

[0040] An anisotropic distance metric is constructed to replace the Euclidean distance in Kriging interpolation for spatial interpolation. The anisotropic distance metric projects the coordinate difference vector between the point to be interpolated and the known stratum boundary point onto the three principal axes of the interpolation ellipsoid. Each projected component is divided by its corresponding principal axis radius, and the sum of the squares is taken as the square root. The three principal axes of the interpolation ellipsoid are the major axis, minor axis, and normal axis, with the normal axis parallel to the normal vector of the stratum surface. The major axis radius is a range parameter determined by the hysteresis distance at which the spatial autocorrelation of the stratum boundary point along the major axis reaches a stationary state. The minor axis radius is determined by dividing the major axis radius by the stretching ratio. The normal radius is determined by dividing the major axis radius by the normal compression ratio, which is the square of the dip cotangent.

[0041] Kriging interpolation, taking anisotropic parameters into account, is used to spatially interpolate the strata boundaries after outlier removal, generating smooth reference surfaces for each stratum. In the variogram of the Kriging interpolation, the distance parameter is replaced by anisotropic distance metric instead of Euclidean distance. The variogram model is selected from spherical, exponential, and Gaussian models based on the minimum cross-validation residual criterion. During the Kriging interpolation solution, a system of Kriging equations is established for each interpolated grid node, consisting of known strata boundaries in its neighborhood. The coefficient matrix of the Kriging equations is jointly constructed by the variogram values ​​under the anisotropic distance metric and Lagrange multiplier constraints, yielding elevation estimates for the interpolated grid nodes. The smooth reference surface is constructed by organizing the elevation estimates of all interpolated grid nodes into a regular or triangular grid, with the grid resolution set to half the average spacing between strata boundaries.

[0042] When implementing S2, projection lines are drawn along the normals of the smooth reference surface at each rock stratum boundary point. The directed distance from the rock stratum boundary point to the smooth reference surface is calculated, and this directed distance is used as the normal deviation. The smooth reference surface is organized into a regular mesh or a triangular mesh. When determining the local normal vector of the smooth reference surface at each rock stratum boundary point, if the smooth reference surface is organized into a regular mesh, the normal vector of the face element containing the mesh node closest to the rock stratum boundary point is taken as the local normal vector. If the smooth reference surface is organized into a triangular mesh, the normal vector of the triangular face element containing the projection of the rock stratum boundary point is taken as the local normal vector. Starting from the rock stratum boundary point, a ray is drawn along the local normal vector, and the intersection point of the ray and the smooth reference surface is calculated. The intersection point of the ray and the smooth reference surface is obtained by solving the simultaneous equations of the ray and the plane equations of the smooth reference surface face elements. The ray equation takes the coordinates of the rock stratum boundary point as the starting point of the ray and the local normal vector as the direction of the ray. When the smooth reference surface has local discontinuities, i.e., when the projection of the rock strata boundary point falls on the shared boundary or shared vertex of two or more triangular facets, the normal vector of the triangular facet closest to the rock strata boundary point is taken as the local normal vector for ray intersection. The projection length of the vector from the rock strata boundary point to the intersection point onto the local normal vector is the directed distance. The sign of the directed distance is determined by the spatial position of the rock strata boundary point relative to the smooth reference surface. When the rock strata boundary point is located on the side pointed to by the local normal vector, the directed distance takes a positive value; when the rock strata boundary point is located on the opposite side pointed to by the local normal vector, the directed distance takes a negative value.

[0043] Spatially continuous interpolation is performed on the normal deviations at each stratum boundary point to form a deviation scalar field covering the entire stratum surface. The spatial continuous interpolation uses the same interpolation framework as the Kriging interpolation, which considers anisotropic parameters, used when generating the smooth reference surface, with the normal deviations at each stratum boundary point as the interpolation object. A variogram of the normal deviation is calculated at each stratum boundary point. The distance parameter in the variogram uses anisotropic distance metrics. The three principal axes of the interpolation ellipsoid in the anisotropic distance metrics follow the major axis, minor axis, and normal directions determined when generating the smooth reference surface. The three principal axis radii of the interpolation ellipsoid also follow the major axis, minor axis, and normal radius determined when generating the smooth reference surface. The variogram model is selected from the spherical model, exponential model, and Gaussian model according to the minimum cross-validation residual criterion. The minimum cross-validation residual criterion involves calculating the sum of squared residuals for leave-one-out cross-validation for each of the three variogram models and selecting the variogram model with the smallest sum of squared residuals as the variogram model used for interpolation. During Kriging interpolation, a system of Kriging equations is established for each grid node in the deviation scalar field, consisting of the normal deviations of known rock strata boundaries within its neighborhood. Solving this system yields the deviation estimates for each grid node. The deviation scalar field is composed of the deviation estimates of all grid nodes organized in a regular grid format, with the grid resolution matching that of the smooth reference surface. The coverage area of ​​the rock strata by the deviation scalar field means that the spatial extent of the deviation scalar field coincides with the circumscribed rectangle of the spatial extent of the smooth reference surface projected onto the horizontal plane. This ensures that each grid node within the area covered by the smooth reference surface has a deviation estimate, providing complete input data for subsequent extraction of the deviation sequence along the borehole path.

[0044] By calculating the directed distance from the rock strata boundary point to the smooth reference surface, the geometric difference between the discrete rock strata boundary point and the smooth reference surface is quantified as a normal deviation. The normal deviation is then spatially interpolated to form a deviation scalar field. The deviation scalar field records the degree and direction of deviation of the rock strata boundary point relative to the smooth reference surface in the normal direction in the form of a continuous field. The discontinuous deviation caused by the weak sedimentary surface is manifested as a region in the deviation scalar field with an increased local spatial autocorrelation length, while random deviation is manifested as a region in the deviation scalar field with a shorter spatial autocorrelation length. The spatial continuity of the deviation scalar field provides a data basis for the subsequent analysis of the local autocorrelation length of the deviation sequence along the borehole path.

[0045] Figure 2A flowchart for determining the influence domain of the weak sedimentary surface according to the present invention is provided. In implementation S3, the values ​​of the deviation scalar field at each stratum boundary point are extracted along the borehole path to form a deviation sequence. The deviation scalar field is stored in a regular grid format, and the spatial coordinates of each stratum boundary point correspond to a grid node position in the deviation scalar field. The deviation estimate at that grid node is taken as the value of the deviation scalar field at the stratum boundary point. Extraction along the borehole path refers to sorting the stratum boundary points traversed by a single borehole sequentially according to the borehole depth from smallest to largest, and forming a one-dimensional deviation sequence from the deviation scalar field values ​​corresponding to each sorted stratum boundary point. The arrangement order of the deviation scalar field values ​​in the deviation sequence is consistent with the vertical stacking relationship of the stratum boundary points exposed by the borehole, and the spacing between adjacent deviation scalar field values ​​in the deviation sequence is the depth spacing of the corresponding stratum boundary points along the borehole.

[0046] The local autocorrelation length of each stratum boundary point is obtained by performing a sliding autocorrelation operation on the deviation sequence. The sliding autocorrelation operation involves setting a sliding window on the deviation sequence, and for each stratum boundary point, calculating the autocorrelation function by taking a segment of the deviation sequence within the sliding window that includes that stratum boundary point. The autocorrelation function is calculated as a sequence of correlation coefficients between the deviation sequence segment and itself at different lag times. The lag time at which the autocorrelation function first falls below the autocorrelation threshold is the local autocorrelation length of that stratum boundary point. The autocorrelation threshold is the value of the autocorrelation function at zero lag multiplied by a decay coefficient, which is set empirically, for example, 1 / e, where e is a natural constant. The sliding window length is adaptively adjusted according to the thickness of the strata exposed by the borehole. The adjustment method is to take the sliding window length as the ratio of the stratum thickness to the average spacing between stratum boundaries, rounded up, and then multiplied by the average spacing between stratum boundaries. The stratum thickness is determined by the depth difference between the top and bottom interfaces of the strata in the borehole lithology log. When the sliding window length at the interface of adjacent sedimentary cycles is taken as the smaller of the two values, the sedimentary cycle interface refers to the boundary between adjacent sedimentary cycles revealed by the borehole. The sedimentary cycle is divided by the cyclical change of lithology in the lithological sequence, which is either finer or coarser upwards. The difference in the thickness of the rock layers on both sides of the cycle interface is significant. Taking the sliding window length corresponding to the smaller thickness avoids the sliding window crossing the cycle interface and mixing in the bias signals of different sedimentary environments, thus ensuring the sensitivity of the local autocorrelation length to the weak sedimentary surface signal within a single cycle.

[0047] The location of lithological change interfaces is obtained from the borehole lithology log. The borehole lithology log records lithological segmentation information along the borehole depth. A lithological change interface refers to the boundary between two adjacent lithological segments. The depth coordinates of the lithological change interface are directly read from the borehole lithology log. The stratum boundary point located at the lithological change interface is the stratum boundary point whose spatial coordinates are consistent with the depth of the lithological change interface. The depth of the lithological change interface is compared with the depth coordinates of each stratum boundary point. The stratum boundary point whose depth coordinates are less than the depth of the lithological change interface is marked as located at the lithological change interface. The preset tolerance is set according to the accuracy of the borehole logging depth sampling, for example, 1 / 2 of the borehole logging sampling interval.

[0048] The local autocorrelation lengths of all strata boundaries are statistically analyzed. Using lithology as the stratification variable, the standard deviation of the local autocorrelation lengths within each lithological layer is calculated. The maximum sum of the mean and standard deviation within each lithological layer is used as the statistical baseline. Using lithology as the stratification variable means grouping all strata boundaries according to their respective lithologies, with each group corresponding to a specific lithology determined by borehole lithology logging. For each lithology-specific strata boundary group, the mean and standard deviation of the local autocorrelation lengths are calculated, with the standard deviation representing the sample standard deviation. The sum of the mean and standard deviation of the local autocorrelation lengths represents the threshold at which the statistically high local autocorrelation length within that lithological layer is considered. The maximum sum of the mean and standard deviation of the local autocorrelation lengths for each group is used as the statistical baseline, ensuring that the statistical baseline reflects the highest upper limit of statistical fluctuations in local autocorrelation lengths within each lithological layer, thus avoiding the underestimation of weak sedimentary surface signals due to background differences in local autocorrelation lengths caused by lithological variations.

[0049] Points at lithological change interfaces where the local autocorrelation length exceeds the statistical baseline are designated as seed points. The local autocorrelation lengths of each boundary point at the lithological change interface are compared with the statistical baseline value one by one; boundary points with local autocorrelation lengths greater than the statistical baseline value are marked as seed points. The set of seed points and their spatially adjacent boundary points whose deviation amplitudes are positively correlated with the seed points are defined as the influence domain of the sedimentary weak surface. Spatially adjacent means that the seed point is located within the same lithological layer in three-dimensional space, and the Euclidean distance from the seed point is less than the neighborhood search radius. The neighborhood search radius is twice the average local autocorrelation length within the lithological layer where the seed point is located, ensuring that the search radius is adapted to the spatial extent of the sedimentary weak surface. Deviation amplitude refers to the absolute value of the deviation scalar field's estimated deviation at the lithological layer boundary point. A positive correlation between deviation amplitude and seed point means that the deviation scalar field value of the stratum boundary point and the deviation scalar field value of the seed point have the same sign, i.e., both are positive or both are negative. This indicates that the stratum boundary point and the seed point are in the same direction of deviation relative to the smooth reference surface. This consistency indicates that the deviation is caused by the same sedimentary weak surface structure rather than random fluctuations. Starting from the seed point, spatially adjacent boundary points satisfying the positive correlation of deviation amplitude are searched within the neighborhood search radius and included in the sedimentary weak surface influence domain. For the boundary points included in the sedimentary weak surface influence domain, the search continues to find spatially adjacent boundary points that satisfy the positive correlation of deviation amplitude until the search boundary no longer expands, forming a spatially connected sedimentary weak surface influence domain. The connectivity of the sedimentary weak surface influence domain allows the influence range of the same sedimentary weak surface to be completely extracted in three-dimensional space, providing a complete input boundary for subsequent identification of positive deviation local uplift clusters and quantitative description of the geometric features of sedimentary weak surfaces within the sedimentary weak surface influence domain.

[0050] When implementing S4, a deviation surface is constructed within the influence domain of the sedimentary weak surface, using the normal deviation as the elevation. The influence domain of the sedimentary weak surface, determined by S3, is the set of spatially connected boundary points on the strata. Deviation estimates are extracted from the deviation scalar field at each grid node within the influence domain of the sedimentary weak surface; these deviation estimates are the numerical values ​​of the normal deviation. A deviation surface in the form of a regular grid is constructed using the planar coordinates of each grid node within the influence domain of the sedimentary weak surface as the horizontal and vertical coordinates, and the corresponding deviation estimates as the vertical coordinate. A positive deviation region refers to a continuous region on the deviation surface where the deviation estimate is greater than zero; the positive deviation corresponds to a stratum boundary point located on the side pointed to by the normal vector of the smooth reference surface.

[0051] Watershed segmentation is performed on regions with positive deviations to extract spatially continuous local uplift clusters. Watershed segmentation uses the deviation estimate as the topographic height, treating the positive deviation region as the topographic surface. Separate catchment basins are identified through simulated flooding, with each catchment basin representing a spatially continuous local uplift cluster. The seed point for watershed segmentation is the local maximum of the deviation estimate within the positive deviation region. The water accumulation threshold is set based on the standard deviation of the deviation estimate within the influence zone of the weak sedimentary surface. The water accumulation threshold is a preset multiple of the standard deviation of the deviation estimate, selected to balance over-segmentation and under-segmentation in watershed segmentation; for example, 0.5 times the standard deviation of the deviation estimate. Adjacent catchment basins are merged when the ratio of the merging depth of adjacent catchment basins to the maximum deviation estimate is less than the water accumulation threshold. The resulting independent catchment basins are the extracted positive deviation local uplift clusters.

[0052] For each local uplift cluster, a minimum bounding rectangle is fitted. Each grid node within the local uplift cluster is projected onto a tangent plane of a smooth reference surface within the influence domain of the weak depositional surface. The tangent plane is taken as the tangent plane of the smooth reference surface at the center point of the local uplift cluster. The center point of the local uplift cluster is the projection point of the grid node containing the maximum estimated deviation value within the local uplift cluster onto the smooth reference surface. The minimum bounding rectangle is calculated on the tangent plane using the rotating caliper method for the set of projected points. The minimum bounding rectangle is the rectangle with the smallest area, containing all projected points. The major axis of the minimum bounding rectangle is taken as the extension direction of the local uplift cluster. The extension direction is an angle, representing the dominant orientation of the local uplift cluster within the tangent plane. The difference between the maximum and minimum normal deviation values ​​within the local uplift cluster is taken as the maximum elevation difference. The maximum elevation difference has the dimension of length, characterizing the undulation amplitude of the local uplift cluster in the normal direction.

[0053] The undulation dominance direction is obtained by weighting the extension directions of each local uplift cluster using the projected area as the weight, and the roughness amplitude is obtained by weighting the maximum height difference of each local uplift cluster. The projected area is taken as the orthographic projection area of ​​the local uplift cluster on the tangent plane of the smooth reference surface within the influence domain of the weak deposition surface. The orthographic projection area is determined by the convex hull area of ​​the projection point of the local uplift cluster. In the weighted averaging process of the undulation dominance direction, the extension directions of each local uplift cluster are converted into unit direction vectors. The projected area is used as the weight to weight each direction vector, and the weighted average is obtained to obtain the weighted average resultant vector. The angle calculated from the weighted average resultant vector is then used as the undulation dominance direction. This process avoids the influence of the periodicity of the angle on the calculation of the average direction. The initial weighted average of the roughness amplitude is obtained by directly weighting the maximum height difference of each local uplift cluster using the projected area as the weight.

[0054] To further eliminate the influence of local hump cluster area differences on roughness amplitude, the normalized height difference is obtained by dividing the maximum height difference of each local hump cluster by the square root of its corresponding projected area. The normalized height difference characterizes the degree of undulation per unit area. Then, a weighted average of the normalized height differences is calculated using the projected area as the weight, yielding the normalized roughness. The roughness amplitude is calculated as follows: Where R represents the roughness amplitude, rn represents the normalized roughness, and Attotal represents the total projected area of ​​the sedimentary weak surface influence domain. The total projected area of ​​the sedimentary weak surface influence domain is the orthographic projection area of ​​all grid nodes within the influence domain onto the tangent plane of the smooth reference surface. The roughness amplitude has the dimension of length, reflecting both the undulation intensity of the positive deviation local uplift clusters within the sedimentary weak surface influence domain and matching the overall spatial scale of the influence domain. This allows the roughness amplitude to be used as the amplitude parameter for the undulation term subsequently superimposed onto the smooth reference surface, directly embedded in the 3D geological model construction process.

[0055] In implementing S5, discrete points on the smooth reference surface are used as carriers. An offset is added at each discrete point along the undulation-dominant direction. The magnitude of the offset is determined by the product of the roughness magnitude and the projection component of the undulation-dominant direction onto the tangent plane of the discrete point, forming the undulation term. The discrete points on the smooth reference surface are all grid nodes of the smooth reference surface mesh, and each grid node carries three-dimensional coordinates and normal vector information. The undulation-dominant direction is determined by S4, and is a unit direction vector embedded in the tangent plane of the smooth reference surface where the depositional weak surface influence domain is located. For each grid node on the smooth reference surface, it is determined whether the grid node falls within the spatial range of the depositional weak surface influence domain. If the horizontal projection coordinates of the grid node are within the projection range of the depositional weak surface influence domain on the horizontal plane, then the grid node belongs to the depositional weak surface influence domain and requires an offset; if the grid node does not belong to the depositional weak surface influence domain, the offset is zero. This conditional judgment ensures that the undulation term is only applied within the depositional weak surface influence domain, and the rest of the smooth reference surface remains unchanged, achieving the embedding of local discontinuous geometric features.

[0056] For grid nodes belonging to the influence domain of weak depositional surfaces, the offset is calculated as: D = R × cosα; where D represents the offset, R represents the roughness amplitude, α represents the angle between the dominant undulation direction and the reference direction in the tangent plane at the grid node, and cosα represents the projection component of the dominant undulation direction in the tangent plane of the grid node. The roughness amplitude is determined by S4, and is a scalar with length dimensions. The projection component cosα reflects the geometric matching relationship between the dominant undulation direction and the tangent plane of the grid node. When the dominant undulation direction is parallel to the tangent plane of the grid node, cosα reaches its maximum value, and the complete roughness amplitude is added in the offset direction; when there is an angle between the dominant undulation direction and the tangent plane of the grid node, only the effective component of the roughness amplitude in the tangent plane is added, so that the magnitude of the offset adaptively matches the local tangent plane direction.

[0057] The undulation term consists of the offsets at all grid nodes in the smooth reference surface grid that belong to the influence domain of the weak deposition surface. The undulation term of each grid node is a vector with the direction along the local normal vector direction of the smooth reference surface at the grid node and the magnitude equal to the offset D.

[0058] The undulation term is superimposed on the coordinates of the smooth reference surface to generate a discontinuous layer with a locally discontinuous undulation morphology. The superposition process involves vector-adding the 3D coordinate vectors of each grid node in the smooth reference surface mesh with the corresponding undulation term vector to obtain new 3D coordinates for the discontinuous layer grid nodes. The topological connectivity of the discontinuous layer mesh is completely consistent with that of the smooth reference surface mesh; only the grid node coordinates undergo normal offset within the influence domain of the weak depositional surface. Within the influence domain of the weak depositional surface, the discontinuous layer exhibits a local bulge or depression morphology along the direction of undulation dominance. The sign of the bulge is determined by the sign of the corresponding normal deviation. A positive normal deviation results in a bulge formed by an additional offset along the positive normal direction, while a negative normal deviation results in a depression formed by an additional offset along the negative normal direction, preserving the discontinuous geometric characteristics of the weak depositional surface.

[0059] Using discontinuous strata as stratigraphic interfaces to constrain solid modeling, a 3D geological model incorporating the geometric features of strata boundaries is generated. During solid modeling, two adjacent discontinuous strata are used as the top and bottom interfaces of a stratigraphic unit, respectively. Triangular prisms or tetrahedral solid units are constructed by connecting corresponding mesh nodes on the top and bottom interfaces, forming closed stratigraphic unit solids. If multiple strata exist, the solid modeling process is repeated for each pair of adjacent discontinuous strata to obtain solid models of each stratigraphic unit. The solid models of all stratigraphic units are combined to form the 3D geological model. In the 3D geological model, the node coordinates of each stratigraphic unit mesh inherit the local geometric features of the discontinuous strata. At the influence zone of weak sedimentary surfaces, the strata exhibit non-smooth undulations, fully preserving the geometric features of the weak sedimentary surfaces.

[0060] The generated 3D geological model serves as the input model for numerical simulation of strata movement. Discontinuous bedding provides geometric constraints for simulating strata separation along sedimentary weak surfaces and bedding-parallel slip. The stratigraphic interfaces constructed with discontinuous bedding exhibit realistic micro-undulations in the influence domain of sedimentary weak surfaces. This overcomes the problems caused by using continuous smooth interpolation for 3D modeling of bedding, which results in the smooth filling of the geometric features of sedimentary weak surfaces and the inability to simulate the overall continuous deformation of strata during simulation, thus preventing the generation of slippage and separation of the hanging wall along sedimentary weak surfaces. This allows the numerical simulation to reproduce the discontinuous responses of interlayer delamination and step-like subsidence during strata movement.

[0061] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0062] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0063] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and inventive constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

Claims

1. An automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model, characterized in that, Includes the following steps: S1. Perform smooth interpolation on the rock strata boundary points extracted from boreholes and profiles to generate smooth reference surfaces for each rock stratum. S2. Calculate the normal deviation of each rock layer boundary point relative to the smooth reference surface to form a deviation scalar field; S3. Calculate the local autocorrelation length of the deviation scalar field along the borehole path. Identify the adjacent boundary points at the lithological change interface where the local autocorrelation length exceeds the statistical base value and the deviation amplitude is positively correlated as the influence domain of the sedimentary weak surface. S4. Within the influence domain of the weak depositional surface, identify the positive deviation local uplift clusters on the curved surface constructed by the normal deviation, take the extension direction of each uplift cluster as the undulation dominance direction by area weighted average, and take the maximum height difference of each uplift cluster as the roughness amplitude by area weighted average. S5. The roughness amplitude is constructed as an undulation term along the direction of the undulation dominance, and superimposed on the smooth reference surface to generate a discontinuous layer that retains the geometric characteristics of the weak sedimentary surface. The discontinuous layer is used as the stratigraphic interface to construct a three-dimensional geological model for numerical simulation of rock strata movement.

2. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 1, characterized in that, Smooth interpolation is performed on the rock strata boundary points extracted from boreholes and profiles to generate smooth reference surfaces for each rock stratum, including: Outlier detection is performed on the rock strata boundary points to remove abnormal boundary points that deviate from the statistical distribution; based on the rock strata attitude revealed by the borehole, the anisotropy parameters of the strata are calculated, and Kriging interpolation that takes into account these anisotropy parameters is used to perform spatial interpolation on the removed rock strata boundary points to generate a smooth reference surface for each rock stratum.

3. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 2, characterized in that, The anisotropy parameters for calculation include: extracting strike azimuth and dip angles from the rock strata exposed by the borehole; determining the major axis direction of the interpolation ellipsoid using the strike azimuth; determining the stretching ratio of the major and minor axes using the tangent of the dip angle; and constructing anisotropic distance metrics to replace the Euclidean distance in Kriging interpolation for spatial interpolation.

4. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 1, characterized in that, Calculate the normal deviation of each rock stratum boundary point relative to the smooth reference surface to form a deviation scalar field, including: Project the projection lines along the normal of each rock layer boundary point on the smooth reference surface, calculate the directed distance from the rock layer boundary point to the smooth reference surface, and use this directed distance as the normal deviation; perform spatial continuous interpolation on the normal deviation of each rock layer boundary point to form a deviation scalar field covering the range of the rock layer.

5. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 4, characterized in that, The directed distance from the rock stratum boundary point to the smooth reference surface is obtained by: taking the rock stratum boundary point as the starting point, drawing a ray along the local normal vector of the smooth reference surface, calculating the intersection point of the ray and the smooth reference surface, taking the projection length of the vector from the rock stratum boundary point to the intersection point onto the normal vector as the directed distance, and taking the normal of the nearest triangular element at the discontinuity of the intersection point on the smooth reference surface.

6. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 1, characterized in that, The local autocorrelation length of the deviation scalar field is calculated along the borehole path. Adjacent boundary points at lithological change interfaces, whose local autocorrelation lengths exceed the statistical baseline and whose deviation amplitudes are positively correlated, are identified as the influence domains of sedimentary weak surfaces, including: The values ​​of the deviation scalar field at each stratum boundary point are extracted along the borehole path to form a deviation sequence. The local autocorrelation length of each stratum boundary point is obtained by performing sliding autocorrelation on the deviation sequence. The location of the lithological change interface is obtained from the borehole lithology logging. The local autocorrelation lengths of all stratum boundary points are counted, and the sum of their mean and standard deviation is used as the statistical base value. The boundary points at the lithological change interface where the local autocorrelation length exceeds the statistical base value are used as seed points. The set of seed points and their spatially adjacent boundary points whose deviation amplitude is positively correlated with the seed points are determined as the influence domain of the sedimentary weak surface.

7. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 6, characterized in that, When calculating the local autocorrelation length of all strata boundaries, the standard deviation of the local autocorrelation length within each stratum is calculated using lithology as the stratification variable. The maximum value of the sum of the mean and the standard deviation within each stratum is used as the statistical base value. During the sliding autocorrelation calculation, the sliding window length is adaptively adjusted according to the thickness of the strata exposed by the borehole, and the window length at the interface of adjacent sedimentary cycles is taken as the smaller of the two values.

8. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 1, characterized in that, Within the influence domain of the weak depositional surface, local uplift clusters with positive deviations are identified on the surface constructed by the normal deviation. The extension direction of each uplift cluster is weighted by area as the dominant undulation direction, and the maximum height difference of each uplift cluster is weighted by area as the roughness amplitude, including: Within the influence domain of the weak sedimentary surface, a deviation surface is constructed using the normal deviation as the elevation. Watershed segmentation is performed on the positive deviation region to extract spatially continuous local uplift clusters. For each local uplift cluster, a minimum bounding rectangle is fitted, and the major axis direction of the minimum bounding rectangle is taken as the extension direction of the local uplift cluster. The difference between the maximum and minimum normal deviations within the local uplift cluster is taken as the maximum elevation difference. Using the projected area of ​​each local uplift cluster as the weight, the extension directions of each local uplift cluster are weighted and averaged to obtain the undulation dominance direction. The maximum elevation difference of each local uplift cluster is weighted and averaged to obtain the roughness amplitude.

9. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 8, characterized in that, When the projected area is used as the weight, the projected area is the orthogonal projected area of ​​the local uplift cluster on the tangent plane of the smooth reference surface within the influence domain of the weak sedimentary surface. The normalized height difference is obtained by dividing the maximum height difference of each local uplift cluster by the square root of the corresponding projected area. Then, the normalized height difference is weighted and averaged using the projected area as the weight. The product of the weighted average result and the square root of the total projected area of ​​the influence domain of the weak sedimentary surface is used as the roughness amplitude.

10. The automated modeling method for numerical simulation of rock strata movement based on a three-dimensional model according to claim 1, characterized in that, Roughness amplitude is constructed as a undulation term along the direction of undulation dominance and superimposed onto a smooth reference surface to generate a discontinuous layer that preserves the geometric characteristics of the weak sedimentary surface. This discontinuous layer is then used as the stratigraphic interface to construct a three-dimensional geological model for numerical simulation of strata movement, including: Using discrete points on a smooth reference surface as carriers, an offset is added at each discrete point along the direction of undulation dominance. The magnitude of the offset is determined by the product of the roughness magnitude and the projection component of the undulation dominance direction in the tangent plane at that point, forming an undulation term. The undulation term is superimposed with the coordinates of the smooth reference surface to generate a discontinuous layer with local discontinuous undulation morphology. The discontinuous layer is used as a stratigraphic interface constraint entity to model a three-dimensional geological model containing the geometric features of rock layer boundaries.

Images