CN114255188B

CN114255188B - Boundary linear smoothing method and device for drawing three-dimensional geological lithologic grid model

Info

Publication number: CN114255188B
Application number: CN202111592582.3A
Authority: CN
Inventors: 牛露佳; 王占刚; 周英
Original assignee: China University of Mining and Technology Beijing CUMTB
Current assignee: China University of Mining and Technology Beijing CUMTB
Priority date: 2021-12-23
Filing date: 2021-12-23
Publication date: 2022-06-24
Anticipated expiration: 2041-12-23
Also published as: CN114255188A

Abstract

The invention provides a boundary linear smoothing method and device for drawing a three-dimensional geological lithologic grid model. The method comprises the steps of merging grids which represent the same rock stratum attribute in the Z-axis direction in a three-dimensional geological lithology grid model by using run length codes to obtain a plurality of rock stratum sections; taking each rock stratum section in the Z-axis direction as a vertical column, and connecting the starting points of the rock stratum sections with space connection relation in each vertical column by using straight lines; the method comprises the steps that rock stratum sections with space connection relations in four adjacent vertical columns are encircled to form a hexahedron, the dividing mode of the side face is determined according to the attribute relation of the rock stratum sections in the vertical columns to which two rock stratum sections in the side face of the hexahedron belong along the direction of an X axis or a Y axis, and the dividing mode of dividing the hexahedron into tetrahedrons is obtained by combining the diagonals of the top face and the bottom face of the hexahedron; manufacturing a calculation template of the rock stratum attribute of the sampling point in the hexahedron according to the subdivision mode; and obtaining the rock stratum attribute of the sampling point by calculating the template, and determining the rock stratum attribute of each position of the hexahedron.

Description

Boundary linear smoothing method and device for drawing three-dimensional geological lithologic grid model

Technical Field

The invention relates to the technical field of engineering geological structure modeling, in particular to a boundary linear smoothing method and device for drawing a three-dimensional geological lithologic grid model.

Background

At present, the construction of a three-dimensional geological structure model is an important content in urban three-dimensional geological survey, and the fine lithological description of the three-dimensional geological structure model on the urban underground shallow layer geological structure and the geological environment of an engineering field can provide basic three-dimensional data support for solving geological problems related to engineering activities. At present, the modeling method of the three-dimensional geological structure model mainly uses an explicit modeling or implicit modeling method to solve the modeling problem with clear stratigraphic structure characteristics, and an irregular triangle table is adopted to reach the boundary of a geological object. However, in the fields of engineering geology and the like, because the stratum structure of the shallow layer fourth system is not obvious, spatial interpolation methods such as kriging indication, sequential simulation, Support Vector Machine (SVM for short) and the like are mostly adopted for a three-dimensional geological structure model, a three-dimensional regular grid is generally adopted for uniform dispersion and interpolation of an underground space in the modeling process, and the modeling result is a three-dimensional lithology regular grid model or a lithology boundary triangulation network extracted based on the regular grid.

Compared with the method for converting the three-dimensional geological structure model into the triangular net, the method for expressing the lithology model based on the three-dimensional regular grid is more convenient for spatial analysis and calculation. For example, a geological grid model is directly constructed by means of an aeroelectromagnetic method and the like, visual carving and erosion analysis of mixed terrain and geological structures, a multilayer grid model and a visual system based on a sketch and the like. In consideration of visualization angles, the drawing of the three-dimensional geological structure model is mainly realized by directly drawing mesh units or converting the mesh units into triangular meshes by using a surface drawing method, but the number of the triangular units generated by the method is more, and better conversion effect can be met by performing more simplified operations. In addition, the transformation into triangular meshes does not allow more spatial analysis in the visualization phase. For these reasons, some researchers have implemented direct volume rendering techniques for visualization and spatial analysis of geological grid models, etc., and with the development of faster volume rendering techniques and HTML5 techniques, network techniques can be used to provide interactive three-dimensional model visualization to everyone on each device.

Different from medical CT volume data, the attribute distribution of the lithological grid model has layering and definite boundary, and the display of the whole model and partial lithological geological interface (isosurface) can be easily realized by utilizing algorithms such as ray casting and the like. The ray projection algorithm needs to perform spatial sampling along the ray direction, the most common sampling mode is nearest neighbor sampling, but when the attribute value at the boundary changes, an obvious sawtooth or step phenomenon can be generated, and the drawing effect is influenced. For volume data with continuously distributed attribute values, such as medical CT volume data or physical parameter volume data, a spatial interpolation mode can be adopted to perform antialiasing processing, such as trilinear interpolation and TriCubic interpolation. However, the lithological values are not continuous and need to ensure the superposition rationality and consistency according to the spatial position of the rock stratum in the longitudinal direction, the position of the rock stratum boundary point is generally not allowed to be changed, so that the tri-linear interpolation and the tricobic interpolation cannot be used, and the global spatial interpolation method is not suitable for the real-time requirement of volume rendering.

The invention provides a boundary linear smoothing method for drawing a three-dimensional lithologic grid model according to the horizontal continuous characteristics of a rock stratum and the reasonability requirement of a space superposition relationship in the longitudinal direction.

Disclosure of Invention

An object of the first aspect of the present invention is to overcome at least one technical defect in the prior art, and to provide a boundary linear smoothing method and apparatus for three-dimensional geological lithology grid model volume rendering.

It is an object of the present invention to smooth the boundaries of a formation.

It is another object of the invention to smooth the boundaries of a three-dimensional geological lithology mesh model.

According to an aspect of the present invention, there is provided a method of linear smoothing of boundaries drawn by a three-dimensional geological lithology mesh model volume, the three-dimensional geological lithology mesh comprising meshes distributed along an X-axis, a Y-axis and a Z-axis of a coordinate system, the method comprising:

merging grids which represent the same rock stratum attribute in the Z-axis direction in the three-dimensional geological lithology grid model by using run length coding to obtain a plurality of rock stratum segments;

taking each rock stratum section in the Z-axis direction as a vertical column, and connecting the starting points of the rock stratum sections with space connection relation in each vertical column by using straight lines;

the method comprises the steps that rock stratum sections with space connection relations in four adjacent vertical columns are encircled to form a hexahedron, the dividing mode of the side face is determined according to the attribute relation of the rock stratum sections in the vertical columns to which two rock stratum sections in the side face of the hexahedron belong along the direction of an X axis or a Y axis, and the dividing mode of dividing the hexahedron into tetrahedrons is obtained by combining diagonals of the top face and the bottom face of the hexahedron; wherein, the adjacent four vertical columns refer to four vertical columns which can enclose a minimum rectangle;

manufacturing a calculation template of the rock stratum attribute of the sampling point in the hexahedron according to the subdivision mode;

and obtaining the rock stratum attribute of the sampling point through the calculation template, thereby determining the rock stratum attribute of each position in the hexahedron.

Optionally, the step of determining the segmentation mode of the side face according to the attribute relationship of the rock stratum sections in the vertical column to which the two rock stratum sections in the side face of the hexahedron belong includes:

for any one side face, judging whether the rock stratum properties of two rock stratum sections in the side face are the same or not;

if the two rock intervals are different, selecting any one of the two rock intervals as a target rock interval, and judging whether the upper rock interval or the lower rock interval connected with the starting point or the end point of the target rock interval and the other rock interval are overlapped in the Z-axis direction and whether the rock attributes are the same;

if the target rock stratum section and the upper rock stratum section are overlapped in the Z-axis direction and the rock stratum properties are the same, connecting the tail point of the target rock stratum section with the starting point of another rock stratum section;

and if the target rock interval and the lower rock interval are overlapped in the Z-axis direction and the rock stratum properties are the same, connecting the starting point of the target rock interval with the end point of the other rock interval.

Optionally, after the step of determining whether the upper rock interval or the lower rock interval connected to the start point or the end point of the target rock interval and the other rock interval overlap in the Z-axis direction and the rock attributes are the same, if the target rock interval, the upper rock interval and the lower rock interval do not satisfy the overlapping in the Z-axis direction and the rock attributes are the same, the step of determining the dividing manner of the side face according to the attribute relationship of the rock interval in the vertical column to which the two rock intervals belong in the side face of the hexahedron further includes:

and judging the height of the end points of the two rock intervals in the side surface in the Z-axis direction, selecting the end point higher in the Z-axis direction, and connecting the end point with the starting point of the other rock interval to form a cutting line.

searching a higher end point in the Z-axis direction in the end points of the two rock stratum sections in the side face, and judging whether the end point of the rock stratum section connected with the higher end point is still higher than another lower end point in the side face;

if yes, selecting the end point of the rock stratum section connected with the higher end point, and connecting the end point of the rock stratum section with the start point of another rock stratum section with the lower end point in the side face to form a cutting line.

Optionally, after the step of determining the heights of the end points of the two rock stratum sections in the side face in the Z-axis direction, if the heights of the end points of the two rock stratum sections in the side face in the Z-axis direction are the same, the step of determining the dividing manner of the side face according to the attribute relationship of the rock stratum sections in the vertical column to which the two rock stratum sections in the side face of the hexahedron belong further includes:

and calculating the lengths of two diagonal lines of the side surface, comparing the lengths of the two diagonal lines, and selecting the longer diagonal line as the cutting line.

Optionally, after the step of obtaining a splitting manner of dividing the hexahedron into tetrahedrons by combining the diagonals of the top surface and the bottom surface of the hexahedron, the method further includes:

and carrying out binary coding according to the positive and negative slope of the object line of each side face in the hexahedron as the coding of the subdivision mode of the hexahedron.

Optionally, the step of manufacturing a computation template of the rock stratum attribute of the sampling point in the hexahedron according to the subdivision mode includes:

making a section parallel to the Y-axis direction through the sampling points, wherein the section intersects with the top surface, the bottom surface and two side surfaces of the hexahedron so as to form a plurality of intersecting lines with the hexahedron and the tetrahedron;

manufacturing the calculation template according to the connection mode of each intersection line; the calculation template comprises 4 partitions along the Y-axis direction, and the head and tail intersection points of each intersection line in each partition are recorded in the partitions.

Optionally, the step of obtaining the formation property of the sampling point through the computation template includes:

determining the partition to which the sampling point belongs as a target partition;

determining coordinates of each intersection point in the target partition according to each vertex coordinate of the hexahedron, so as to determine the position relation between the sampling point and each intersection line in the target partition;

obtaining the rock stratum attribute between two adjacent intersecting lines in the target subarea according to the rock stratum attribute of each rock stratum section in the target subarea and the pinch-out mode;

and obtaining the rock stratum attribute of the sampling point according to the position relation of the sampling point and each intersection line in the target partition and the rock stratum attribute between two adjacent intersection lines in the target partition.

Optionally, after the step of obtaining the rock formation property of the sampling point through the calculation template, so as to determine the rock formation property at each position in the hexahedron, the method further includes:

judging whether a hollow rock stratum section with empty lithology exists in a hexahedron in the three-dimensional geological lithology grid model;

if the vertical column is located at the starting point or the end point of the vertical column, determining the position of the starting point or the end point of the vertical column, taking the elevation of the hexahedral hollow rock stratum as the elevation of the starting point or the end point of the vertical column where the hollow rock stratum is located, and correcting the elevation of the rock stratum with the lithology being non-empty into the elevation of the starting point or the end point of the vertical column where the rock stratum is located; and/or

Judging whether a vertical column with empty lithology exists in the vertical column where the rock stratum section in the side face of the hexahedron in the three-dimensional geological lithology grid model;

if yes, judging whether the number of the empty vertical columns is more than 1;

if the number of the vertical columns is larger than 1, the hexahedron where the hollow vertical columns are located is not processed;

if the number of the rock layer sections is equal to 1, respectively connecting the starting points and the end points of the two rock layer sections adjacent to the hollow vertical column in the hexahedron so as to divide the hexahedron into two triangular columns and reserving the triangular columns which do not contain the hollow vertical column.

According to another aspect of the present invention, there is also provided a boundary linear smoothing apparatus for volume rendering of a three-dimensional geological lithology mesh model, comprising a memory, a processor and a machine executable program stored in the memory and run on the processor, and the processor, when executing the machine executable program, implements the boundary linear smoothing method for volume rendering of a three-dimensional geological lithology mesh model according to any one of the above items.

In the boundary linear smoothing method for drawing the three-dimensional geological lithology grid model, grids which represent the same rock stratum attribute in the Z-axis direction in the three-dimensional geological lithology grid model are combined by using run length codes to obtain a plurality of rock stratum sections; taking each rock stratum section in the Z-axis direction as a vertical column, and connecting the starting points of the rock stratum sections with space connection relation in each vertical column by using straight lines; the method comprises the steps that rock stratum sections with space connection relations in four adjacent vertical columns are encircled to form a hexahedron, the dividing mode of the side face is determined according to the attribute relation of the rock stratum sections in the vertical columns to which two rock stratum sections in the side face of the hexahedron belong along the direction of an X axis or a Y axis, the diagonal lines of the top face and the bottom face of the hexahedron are combined to obtain a dividing mode for dividing the hexahedron into tetrahedrons, and a calculation template of rock stratum attributes of sampling points in the hexahedron is manufactured according to the dividing mode; the rock stratum attribute of the sampling point is obtained through calculating the template, so that the rock stratum attribute of each position in the hexahedron is determined, the problems of sawtooth and steps of a rock stratum interface are obviously solved, the linear smoothness of the boundary drawn by the three-dimensional geological lithology grid model is realized, and the drawing effect of the model is improved.

Further, judging whether a hollow rock stratum section with empty lithology exists in a hexahedron in the three-dimensional geological lithology grid model; if the vertical column is in the horizontal extension control range of the rock stratum, the connection rule and the pinch-out condition of the rock stratum are considered, the sawtooth problem caused by nearest neighbor sampling and the tail knock problem do not exist.

The above and other objects, advantages and features of the present invention will become more apparent to those skilled in the art from the following detailed description of specific embodiments thereof, taken in conjunction with the accompanying drawings.

Drawings

Some specific embodiments of the invention will be described in detail hereinafter, by way of illustration and not limitation, with reference to the accompanying drawings. The same reference numbers in the drawings identify the same or similar elements or components. Those skilled in the art will appreciate that the drawings are not necessarily to scale. In the drawings:

FIG. 1 is a flow diagram of a method for boundary linear smoothing of a three-dimensional geological lithology grid model volume rendering according to one embodiment of the invention;

FIG. 2 is a schematic diagram of initial states and states after processing by run-length encoding of a three-dimensional geological lithology grid model according to another embodiment of the present invention;

FIG. 3 is a schematic illustration of a section of rock formation in the Z-axis direction as a vertical column according to another embodiment of the present invention;

FIG. 4 is a schematic view of a cut line of one side of a hexahedron according to another embodiment of the present invention;

FIG. 5 is a schematic view of a cut line of one side surface of a hexahedron according to another embodiment of the present invention;

FIG. 6 is a schematic view of a cut line of one side surface of a hexahedron according to another embodiment of the present invention;

FIG. 7 is a schematic view of a cut line of one side surface of a hexahedron according to another embodiment of the present invention;

FIG. 8 is a schematic view of a hexahedron being divided according to another embodiment of the present invention;

fig. 9 is a schematic view of a structure for dividing a hexahedron into tetrahedrons according to another embodiment of the present invention;

FIG. 10 is a schematic diagram of a cross-section taken through a sample point and made parallel to the Y-axis according to another embodiment of the present invention;

FIG. 11 is a schematic illustration of a plurality of solutions to the formation property of the middle quadrilateral of FIG. 10;

FIG. 12 is a schematic diagram of a computing template formed from cross-sectional intersection connections according to another embodiment of the present invention;

FIG. 13 is a schematic diagram of sample points located in a partition of a computational template according to another embodiment of the invention;

FIG. 14 is a schematic illustration of smoothing the boundaries of a three-dimensional geological lithology mesh model, according to another embodiment of the present invention;

FIG. 15 is a graphical illustration of a result of a rendering of an engineered geological three-dimensional grid phantom according to another embodiment of the present invention;

FIG. 16 is a schematic representation of a smooth before and after comparison of drawn boundaries of an engineered geological three-dimensional mesh phantom according to another embodiment of the invention;

FIG. 17 is a schematic illustration of comparison results before and after smooth rendering of YZ-plane slice boundaries according to another embodiment of the invention;

FIG. 18 is a diagram showing comparison results before and after XY-plane slice boundary smoothing rendering according to another embodiment of the present invention;

FIG. 19 is an enlarged schematic view of regions A01 and B01 of FIG. 18;

FIG. 20 is a schematic illustration of an explanation of formation interdigitation according to another embodiment of the invention;

FIG. 21 is a schematic diagram of respective calculation templates of rock formation properties in a hexahedron according to another embodiment of the present invention;

FIG. 22 is a schematic illustration of two rock intervals on one side of a hexahedron undergoing a rock pinch-out according to another embodiment of the present invention;

FIG. 23 is a lithology effect graph processed by a boundary linear smoothing method of three-dimensional geological lithology mesh model volume rendering according to another embodiment of the invention.

Detailed Description

Fig. 1 is a flowchart of a boundary linear smoothing method for volume rendering of a three-dimensional geological lithology mesh model according to an embodiment of the present invention, the three-dimensional geological lithology mesh includes meshes distributed along an X-axis, a Y-axis, and a Z-axis of a coordinate system, wherein the X-axis and the Y-axis may be understood as a transverse direction, and the Z-axis may be understood as a longitudinal direction. Of course, the X-axis direction may be regarded as the I-direction and the Y-axis direction may be regarded as the J-direction. A general stratigraphic or lithologic model is characterized by horizontal stratigraphy, i.e., longitudinal stratification and transverse boundary waviness. Referring to fig. 1, the processing method may include the following steps S102 to S110.

Step S102: and merging grids which represent the same rock stratum attribute in the Z-axis direction in the three-dimensional geological lithology grid model by using the run length codes to obtain a plurality of rock stratum sections.

In this step, in the Z-axis direction, the grids representing the same formation property that are merged together are generally adjacent grids, and if two grids representing the same formation property exist between the grids, the two grids are not merged. The X-axis and Y-axis directions in the three-dimensional geological lithology grid model are organized by adopting two-dimensional regular grids, and the grid resolution is consistent with that of the original regular grids. The start point of a rock interval may be understood as the end point of the upper portion and the end point may be understood as the end point of the lower portion. Two adjacent rock intervals up and down, the end point of the upper rock interval is also the starting point of the lower rock interval. After the grids which represent the same rock stratum attribute in the Z-axis direction are combined by the run length codes, only the starting point height and the attribute codes of each rock stratum section in the Z-axis direction can be recorded, which is equivalent to recording only the top and bottom boundaries of each rock stratum section. The attribute codes refer to codes of different rock stratum attributes, for example, the code of filling soil can be 101, the code of filling soil can be 102, the code of new loess can be 103, the code of old loess can be 104, the code of cohesive soil can be 105, and the like, which are not listed here. Run-length encoding merges grids representing the same formation property in the Z-axis direction, which may also be understood as run-length encoding compresses grids representing the same formation property in the Z-axis direction. Compared with the methods such as octree and the like, the run length coding has better compression performance on the three-dimensional geological lithology grid model. Specifically, for example, referring to fig. 2, fig. 2a on the left side of fig. 2 is an initial state of the three-dimensional geological lithology grid model, and fig. 2b on the right side is a state after processing by run-length encoding. The inactive area in fig. 2a does not belong to a rock layer and does not participate in the merging. Through the above description, those skilled in the art can understand that the specific meaning of combining the grids representing the same formation property in the Z-axis direction in the present embodiment results in multiple rock intervals.

Step S104: and taking each rock stratum section in the Z-axis direction as a vertical column, and connecting the starting points of the rock stratum sections with the space connection relation in each vertical column by using straight lines.

In this step, the spatial connection relationship may be understood as being adjacent to each other and having an intersection (overlap) in the Z-axis direction, for example, two rock stratum segments are located on adjacent vertical columns, and the two rock stratum segments have an intersection in the Z-axis direction, and the two rock stratum segments may be considered as having the spatial connection relationship. In this embodiment, the beginning of a rock interval may be understood as the upper end point and the end point may be understood as the lower end point. Two adjacent rock stratum sections up and down, the end point of the upper rock stratum section is also the starting point of the lower rock stratum section. For example, referring to fig. 3, fig. 3 includes four vertical columns including a vertical column 0, a vertical column 1, a vertical column 2, and a vertical column 3, each vertical column includes three rock strata of a rock stratum section 1, a rock stratum section 2, and a rock stratum section 3, and the starting points of the rock strata having a spatial connection relationship in each vertical column are connected by a straight line. The grid after the vertical units are combined can be regarded as a series of horizontal regularly distributed virtual drilling holes. These virtual boreholes can be referred to as vertical pillars (pilar) by the structural grid pilargrid terminology. In addition, the starting points are connected by straight lines, which belongs to a linear processing method and can reduce or eliminate the problem of visual saw teeth caused by regular grid division.

Step S106: the method comprises the steps that rock stratum sections with space connection relations in four adjacent vertical columns are encircled to form a hexahedron, the dividing mode of the side face is determined according to the attribute relation of the rock stratum sections in the vertical columns to which two rock stratum sections in the side face of the hexahedron belong along the X-axis or Y-axis direction, and the dividing mode of dividing the hexahedron into tetrahedrons is obtained by combining the diagonals of the top face and the bottom face of the hexahedron. Wherein, adjacent four vertical columns refer to four vertical columns which can enclose a minimum rectangle.

In the step, the three-dimensional geological lithologic grid model can be divided into a plurality of hexahedrons by surrounding rock stratum sections with space connection relations in four adjacent vertical columns. In the invention, for the prediction of the rock stratum property of any point in the geological space, only the hexahedron where the point is located can be considered, that is, the rock stratum property of the point can be obtained only by considering 4 vertical columns adjacent to the hexahedron. Because the hexahedrons are uniformly split according to the X-axis or Y-axis direction, the consistency of the rock stratum in the processing space of different adjacent hexahedrons can be ensured. Here, the X-axis or Y-axis direction may be understood as the coincidence in a single direction, for example, taking the X-axis direction as the positive direction and the next grid as X +1 for one grid X, or taking the Y-axis direction as the positive direction and the next grid as Y +1 for one grid Y. The property relationship of the rock interval may include a relationship of rock formation properties of the rock interval, such as whether the rock formation properties are the same, or may include a relationship between starting points and end points of the rock interval, such as a height between two starting points and a height between two end points, i.e. a positional relationship. The same formation may be understood as a formation that is laterally continuous.

Step S108: and manufacturing a calculation template of the rock stratum attribute of the sampling point in the hexahedron according to the subdivision mode.

Step S110: and obtaining the rock stratum attribute of the sampling point by calculating the template, thereby determining the rock stratum attribute of each position in the hexahedron.

In the embodiment, the grids which represent the same rock stratum attribute in the Z-axis direction in the three-dimensional geological lithology grid model are combined by using the run length codes to obtain a plurality of rock stratum segments; taking each rock stratum section in the Z-axis direction as a vertical column, and connecting the starting points of the rock stratum sections with space connection relation in each vertical column by using straight lines; the method comprises the steps that rock stratum sections with space connection relations in four adjacent vertical columns are encircled to form a hexahedron, the dividing mode of the side face is determined according to the attribute relation of the rock stratum sections in the vertical columns to which two rock stratum sections in the side face of the hexahedron belong along the direction of an X axis or a Y axis, the diagonal lines of the top face and the bottom face of the hexahedron are combined to obtain a dividing mode for dividing the hexahedron into tetrahedrons, and a calculation template of rock stratum attributes of sampling points in the hexahedron is manufactured according to the dividing mode; the rock stratum attributes of the sampling points are obtained through the calculation template, so that the rock stratum attributes of all positions in the hexahedron are determined, the problems of saw teeth and steps of a rock stratum interface are obviously solved, the linear smoothness of the boundary drawn by the three-dimensional geological lithology grid model is realized, and the drawing effect of the model is improved.

In an embodiment of the present invention, the step of determining the segmentation mode of the side face according to the attribute relationship of the rock stratum sections in the vertical column to which the two rock stratum sections belong in the side face of the hexahedron may include:

for any side face, judging whether the rock stratum properties of two rock stratum sections in the side face are the same or not;

if the target rock stratum section and the upper rock stratum section are overlapped in the Z-axis direction and the rock stratum properties are the same, connecting the tail point of the target rock stratum section with the starting point of the other rock stratum section;

and if the target rock stratum segment and the lower rock stratum segment are overlapped in the Z-axis direction and the rock stratum properties are the same, connecting the starting point of the target rock stratum segment with the end point of the other rock stratum segment.

In this embodiment, in general, the dividing manner of the side face may be determined according to the attribute relationship between the rock stratum section constituting the side face and the upper rock stratum section and the lower rock stratum section which are adjacent to the rock stratum section constituting the side face up and down. The subdivision of the side face can also be understood as connecting rock strata, namely the subdivision rule can be understood as a rock stratum connection rule and can also be understood as a rock stratum pinch-out mode or a rock stratum pinch-out rule. Specifically, for example, referring to fig. 4, in fig. 4, the two formation sections in the side may be sections B1 and B2. Assume that the formation properties of the B1 and B2 segments are different. And selecting a section B2 as a target rock interval, and judging whether the target rock interval and an upper rock interval (section A1) or a lower rock interval (section C1) connected with the starting point or the end point of the section B1 are overlapped in the Z-axis direction and whether the rock attributes are the same. If the target rock interval and the upper rock interval (A1 interval) have overlapping in the Z-axis direction and the rock stratum properties are the same, connecting the end point of the target rock interval with the starting point of the B1 interval. In this embodiment, the manner of dividing the side surface may be understood as a connected manner, and the dividing condition is a connected condition.

In an embodiment of the present invention, after the step of determining whether the upper rock interval or the lower rock interval connected to the start point or the end point of the target rock interval and the other rock interval overlap in the Z-axis direction and the rock attributes are the same, if the target rock interval, the upper rock interval and the lower rock interval do not satisfy the overlapping in the Z-axis direction and the rock attributes are the same, the step of determining the segmentation mode of the side face according to the attribute relationship of the rock interval in the vertical column to which the two rock intervals in the side face of the hexahedron belong may further include:

and judging the height of the end points of the two rock stratum sections in the side surface in the Z-axis direction, selecting the end point higher in the Z-axis direction, and connecting the end point with the starting point of the other rock stratum section to form a cutting line.

In this embodiment, the target rock interval, the upper rock stratum section and the lower rock stratum section do not satisfy the overlapping in the Z-axis direction and the same rock stratum properties, that is, the communication condition in the above embodiment is not satisfied. Specifically, for example, referring to fig. 5, in fig. 5, two rock intervals (B1 section and B2 section) in the side face may be understood as a current rock interval, a rock interval above the current rock interval may be understood as an upper rock interval (a1 section and a2 section), and a rock interval below the current rock interval may be understood as a lower rock interval (C1 section and C2 section). And taking the section B2 as a target rock stratum, and assuming that the rock stratum properties of the section B2 are different from those of the section B1, the section B2 and the section A1 do not meet the requirement of overlapping in the Z-axis direction and have the same rock stratum properties, and the section B2 and the section C1 also do not meet the requirement of overlapping in the Z-axis direction and have the same rock stratum properties. The end of segment B1 is located higher than the end of segment B2, connecting the end of segment B1 with the beginning of segment B2. And (3) rock stratum sections with different attributes on two adjacent vertical columns, wherein the tail point (the tail point can also be understood as the bottom) of the new rock stratum section is higher than the tail point of the old rock stratum section, namely the seat height of the tail point of the new rock stratum section in the Z-axis direction is higher than the coordinate of the tail point of the old rock stratum section in the Z-axis direction, the tail point higher in the Z-axis direction is selected, and the tail point is connected with the starting point of the other rock stratum section to form a cutting line, so that the transverse pinch-off of the rock stratum is realized. The new formation should be overlaid on the old formation in the order of deposition of the different formations. In this embodiment, the side face subdivision method may also be understood as subdivision according to the sequence condition of the rock interval. Two formation sections in the lateral faces may also be understood as two formation sections in a spatially connected relationship on adjacent vertical pillars. In addition, the four side surfaces of the hexahedron are divided, namely, a quadrangle formed by two rock stratum sections with space connection relations on two adjacent vertical columns is divided, then four adjacent vertical columns are selected, and starting points of the rock stratum sections with space connection relations in the four vertical columns are connected, so that the hexahedron is formed. Namely, in the invention, the hexahedron can be enclosed firstly, the side face of the hexahedron can be split, or the quadrangle formed by two rock stratum segments with space connection relation on two adjacent vertical columns can be split along the X-axis or Y-axis direction, then four adjacent vertical columns are selected, and the hexahedron is enclosed by four quadrangles formed in the four adjacent vertical columns, so that the hexahedron is obtained.

searching a higher end point in the Z-axis direction in the end points of the two rock stratum sections in the side surface, and judging whether the end point of the rock stratum section connected with the higher end point is still higher than another lower end point in the side surface; if yes, selecting the end point of the stratum section connected with the higher end point, and connecting the end point of the stratum section with the start point of another stratum section with the lower end point in the side face to form the cutting line.

In the present embodiment, referring to fig. 6, in fig. 6, the section B2 is taken as the target rock stratum, and assuming that the rock stratum properties of the section B2 and the section B1 are different, the section B2 and the section a1 do not satisfy the overlapping in the Z-axis direction and have the same rock stratum properties, and the section B2 and the section C1 also do not satisfy the overlapping in the Z-axis direction and have the same rock stratum properties. The formation property of the current rock interval (section B1) is different from the formation properties of the current rock interval (section B2) and the lower rock interval (section C2) of the adjacent vertical columns, for example, the formation property of section B1 is B1, the formation property of section B2 is B2, and the formation property of section C2 is C2. If the top (beginning) of segment B2 is higher than the top (beginning) of segment B1 and the end point of segment B1 is lower than the end points of segments B2 and C2, the end point (bottom) of segment B2 moves down to the end point of segment C2, connecting the end point of segment C2 (end point of segment B2 after moving down) and the beginning of segment B1, so that segments B2 and C2 overlay segment B1. If the bottom of the section B2 does not move downwards, the connection mode is shown by a dotted line in FIG. 6, which can lead to filling the section B1 between the section B2 and the section C2, so that the section B2 covers the section B1, and the section C2 is wedged into the section B1, thus the form of the section B1 is too complicated. In this embodiment, the lateral dissection mode can also be understood as dissection according to the principle of continuous coverage of the rock interval. Of course, the present embodiment only deals with the continuous pinch-out of two layers (e.g., B2 and C2), and the wedging condition still occurs for the similar case of more than two layers.

In an embodiment of the present invention, after the step of determining the heights of the end points of the two rock layer sections in the side surface in the Z-axis direction, if the heights of the end points of the two rock layer sections in the side surface in the Z-axis direction are the same, the step of determining the dividing manner of the side surface according to the attribute relationship of the rock layer sections in the vertical column to which the two rock layer sections in the side surface of the hexahedron belong may further include:

and calculating the lengths of two diagonals of the side surface, comparing the lengths of the two diagonals, and selecting the longer diagonal as a cutting line.

In the present embodiment, referring to fig. 7, the section B2 is taken as the target rock formation section, and assuming that the formation properties of the section B2 and the section B1 are different, the section B2 and the section a1 do not satisfy the overlapping in the Z-axis direction and have the same formation properties, and the section B2 and the section C1 also do not satisfy the overlapping in the Z-axis direction and have the same formation properties. The heights of the end points of the B1 segment and the B2 segment in the Z-axis direction are the same, namely the coordinates in the Z-axis direction are the same. The beginning of segment B1 is higher than the beginning of segment B2, and the length of diagonal line A1C2 is greater than the length of diagonal line A2C 1. And selecting a starting point (the starting point of the section B1) higher in the Z-axis direction, and connecting the starting point with the end point of another rock interval (the end point of the section B2) to form a cutting line, namely, dividing the quadrangle according to the principle of the longest diagonal. This case is mainly for the bottom depths of the B1 section and the B2 in the adjacent vertical columns are equal (the heights of the end points in the Z-axis direction are the same), and the thin rock stratum is extinguished to the thick rock stratum by considering the rock stratum thickness weight.

When the side face of the hexahedron is divided through the above embodiment, the side face may be divided preferentially by using the communication condition, the side face may be divided by using the sequence condition if the communication condition is not satisfied, the side face may be divided by using the continuous coverage principle if the sequence condition is not applicable, and the side face may be divided by using the longest diagonal principle if the continuous coverage principle is not applicable.

For ease of understanding, the subdivision of the hexahedron is described below by way of a specific example.

Referring to fig. 8, fig. 8a thereof numbers neighboring four vertical pillars as 0(I, J), 1(I +1, J), 2(I +1, J +1), and 3(I, J +1), respectively, along the X-axis direction (I direction) and the Y-axis direction (J direction) of the grid. The current rock stratum section (B section) is encircled to form a polyhedral unit L with uneven top and bottom surfaces, namely a hexahedron. The formation properties of the B sections on each vertical column may not be exactly the same. Along the I or J direction of the grid, the starting points and the end points of the B sections in 4 adjacent vertical columns are correspondingly connected to form the side faces of 4 quadrangles. Each side surface can be divided into 2 triangles according to the division mode, namely, the triangle division is carried out, and the quadrangle triangulation can also be understood. And performing space division on the hexahedron by utilizing the top and bottom surfaces of the hexahedron and the quadrilateral triangularization construction rock stratum division surface to form a specific lithologic subdivision grid unit, wherein the division surface is shown as fig. 8b in fig. 8. A dividing plane can also be understood as a tetrahedral plane. As can be seen from the unit division result, the rock stratum attribute of the section B and the selection of the side quadrilateral diagonal jointly determine the hexahedron space division result. An example of the hexahedron being divided into tetrahedrons may be seen in fig. 9, and fig. 9 is a structural diagram of the hexahedron being divided into tetrahedrons according to another embodiment of the present invention.

Since there may be a difference in the rock formation properties of the hexahedron 4B segments, the most complicated case is that the rock formation properties of the 4B segments are all different. The side surface formed by two adjacent B sections is divided along the diagonal line with 2 possibilities. Therefore, when 4B segments are not consistent in lithologic property, 16 combinations can be obtained at most when the 4 sides are divided. And aiming at each combination condition, a feasible hexahedron subdivision mode can be given, and a calculation template of the spatial rock stratum attribute in the hexahedron is formulated based on the subdivision modes of the 16 combination conditions. In addition, adjacent hexahedrons have common faces, such as common side faces or bottom faces or top faces. The division of the faces common to adjacent hexahedrons remains uniform. For example, in order to ensure the consistency of the top and bottom common different-plane quadrilateral partitions of the upper and lower adjacent hexahedrons (the two diagonals of the top or bottom surfaces are often different planes), in view of efficiency, the invention does not use a bilinear interpolation method, but uses a diagonal plane with 1(I +1, J) and 3(I, J +1) points as default during template construction to perform unit subdivision first, and a diagonal plane in the other direction is connected to the diagonal plane, as shown in fig. 10. Generally, after the dividing mode of the four side surfaces of the hexahedron is determined, the connecting mode of the top and bottom surfaces to the corner line can be determined.

In an embodiment of the present invention, the step of preparing the calculation template of the rock property of the sampling point in the hexahedron according to the subdivision mode may include:

a section parallel to the Y-axis direction is manufactured through sampling points, and the section is intersected with the top surface, the bottom surface and two side surfaces of the hexahedron, so that a plurality of intersection lines are formed with the hexahedron and the tetrahedron;

manufacturing a calculation template according to the connection mode of each intersection line; the calculation template comprises 4 partitions along the Y-axis direction, and the head and tail intersection points of each intersection line in each partition are recorded in each partition.

In the present embodiment, a cross section parallel to the Y-axis direction (J direction) is made through the sampling point P, as shown in fig. 10. The spatial partitioning of the hexahedron may be multi-solution. The reason for the spatial partitioning into multiple solutions is the lack of formation spatial partitioning information on adjacent vertical columns or topological information constituting specific lithological units within the lithological mesh model. For example, in fig. 9, the strata properties of the sections B of 0(I, J), 1(I +1, J), 2(I +1, J +1) and 3(I, J1+1) are inconsistent, and four tetrahedrons with unknown strata properties surrounded by the dividing surfaces are formed after the diagonal lines of the side surfaces and the top and bottom surfaces of the 4 quadrangles of the hexahedron are connected. The formation properties of the middle quadrilateral of fig. 10 can be seen with a section parallel to the Y direction as uncertain. The rock formation attribute of the quadrangle can be selected from the rock formation attributes of 4B sections, so 4 possible situations exist. As shown in fig. 11, fig. 11 is a schematic diagram of a plurality of solutions of the rock formation property of the middle quadrangle in fig. 10, which can also be understood as a schematic diagram of a plurality of solutions of the hexahedral space division. Fig. 11 includes four cases, fig. 11a, fig. 11b, fig. 11c and fig. 11d, so the spatial partitioning is multi-solution. For the region with unknown rock stratum properties, the rock stratum properties can be determined according to a rock stratum property proportion principle (the relatively most one in the rock stratum property categories in the 4B sections) and a random principle (in the case that the rock stratum properties cannot be determined, one rock stratum property is randomly selected from the four B sections and allocated).

The intersection of the above cross section with the hexahedron and tetrahedron (which can also be understood as lithologic splitting plane) presents three features: and in the characteristic 1, the end points of all the intersecting lines can be obtained by the coordinates of 8 vertexes of the hexahedron and the coordinates of the sampling point in the X direction. And 2, the number and the connection mode of all intersecting lines are not changed along with the change of the position of the cross section in the hexahedron. And 3, all intersecting lines can be assigned to 4 partitions. In general, the division into 4 partitions may be based on the coordinates of the end points of the intersecting lines and the coordinates of the diagonal lines of the top and bottom surfaces of the hexahedron. Therefore, the intersection line end points of all 16 kinds of unit L space division conditions are classified and summarized, and a calculation template is made according to the cross section intersection line connection mode, as shown in FIG. 12. In fig. 12, all the intersection endpoints can be normalized to 13 nodes and divided into 4 partitions along the Y direction (J direction), and the coordinates of the 13 nodes are only illustrative and do not limit the present invention. The sampling point P (point P) is located within one of the partitions, for example, see fig. 13, with point P located in the first partition. Specifically, for example, according to the hexahedral division manner (each division manner may be represented by one code), each partition records each intersection line (head and tail nodes) in the partition from top to bottom in turn. Fig. 13 may be a calculation template of the hexahedron division manner in fig. 9, indicating the connection manner of the respective intersecting lines. The formation properties between adjacent intersecting lines are known from the formation properties of the region where each intersecting line is located, and in general, the region where each intersecting line is located can be understood as a region above the intersecting line and a region below the intersecting line.

In one embodiment of the present invention, the step of obtaining the formation property of the sampling point by calculating the template may include:

determining a partition to which the sampling point belongs as a target partition;

determining coordinates of each intersection point in the target partition according to each vertex coordinate of the hexahedron, and accordingly determining the position relation between the sampling point and each intersection line in the target partition;

In this embodiment, the coordinates of each intersection point in the target partition are determined according to each vertex coordinate of the hexahedron, so as to determine the position relationship between the sampling point and each intersection line in the target partition, that is, the sampling point can be determined to be located between which two adjacent intersection lines, and further the attribute of the sampling point is determined according to the rock stratum attribute between the two intersection lines. The pinch-out mode can refer to the division mode of the side surface in the above embodiment.

In an embodiment of the present invention, after the step of determining the rock stratum property at each position in the hexahedron by calculating the rock stratum property of the sampling point through the template, the method may further include:

and if the vertical column is present, determining the position of the starting point or the end point of the vertical column of the hollow rock layer, taking the elevation of the hexahedral hollow rock layer as the elevation of the starting point or the end point of the vertical column of the hollow rock layer according to the position of the starting point or the end point of the vertical column of the hollow rock layer, and correcting the elevation of the rock layer with the lithology being non-hollow to the elevation of the starting point or the end point of the vertical column of the rock layer.

In this embodiment, if a hollow rock interval exists in the hexahedron, it is necessary to determine whether the hollow rock interval is located at the starting point or the end point of the vertical column where the hollow rock interval is located, if the starting point is located, the elevation of the hollow rock interval is used as the elevation of the starting point of the vertical column where the hollow rock interval is located, and the elevation of the non-hollow rock interval in the hexahedron is used as the elevation of the starting point of the vertical column where the hollow rock interval is located, if the end point is located, the elevation of the hollow rock interval is used as the elevation of the end point of the vertical column where the hollow rock interval is located, and the elevation of the non-hollow rock interval in the hexahedron is used as the elevation of the end point of the vertical column where the hollow rock interval is located. Typically, elevation differences at the top or bottom boundaries of a three-dimensional geological lithology grid model will result in formation voids on certain formation segments that make up the sides of the hexahedron, e.g., see FIG. 14, where the left half of the top formation D of FIG. 14a in FIG. 14 is a void region, which may also be understood as a void region, which may result in non-void formations that do not completely overlie the underburden. If the top of formation D is connected to the top of formation B along a connecting line, which may also be understood as a dashed line in fig. 14B, formation D is not in full contact with the overburden (formation C). Therefore, it is necessary to find a hollow rock interval with empty lithology in the hexahedron in the three-dimensional geological lithology grid model. The rock stratum segments forming the hexahedron are assumed to be four vertical-column B-rock intervals, namely 4B segments. If the section B with the empty lithology is located at the starting point of the vertical column, the elevation of the section B with the empty lithology is set to be the elevation of the starting point of the vertical column, namely the thickness of the section B is zero, the elevation of the section B with the non-empty lithology is set to be the elevation of the starting point of the vertical column, and the processing result is as shown in fig. 14B, so that the corrected rock stratum D can be completely contacted with the overburden and completely covered on the underburden.

In the present embodiment, assuming the engineering geological formation level, in the prior art, since only the adjacent intervals in the horizontal direction of the input point are processed, the attribute calculation process is performed within the control range of the horizontal extension of the formation (as shown in fig. 14 a). The layer connection rule provided by the invention actually provides a connection mode and a rock layer pinch-out mode of rock layer segments on adjacent vertical columns, a rock layer structure model organized according to a PillarGrid grid can be formed (as shown in figure 14c), and due to the huge difference of the top and bottom elevations between rock layers, the direct connection pinch-out can cause a long and narrow obtuse triangle, namely the tail knock problem. Therefore, in the processing of the embodiment, the control range of the horizontal extension of the rock formation is considered, and the connection rule and the pinch-out condition of the rock formation section are considered, and the processing effect is as shown in fig. 14b, and the problem of jaggy caused by nearest neighbor sampling shown in fig. 14a and the problem of tail knock shown in fig. 14c do not exist.

In one embodiment of the present invention, the step of smoothing the boundaries of the three-dimensional geological lithology mesh model may comprise:

if the number of the vertical columns is more than 1, the hexahedron where the empty vertical columns are located is not processed;

In this embodiment, if the number of the empty vertical pillars is greater than 1, it indicates that the hexahedron where the empty vertical pillars are located is outside the three-dimensional geological lithology grid model, and therefore the hexahedron where the empty vertical pillars are located may not be processed. In the present embodiment, the hollow vertical column is understood to be in a hexahedron as long as a part of the hollow vertical column is located in the hexahedron. In addition, the rock formation property of any vertical column adjacent to the empty vertical column can be used as the property of the empty vertical column.

For ease of understanding, the following detailed description of the implementations and testing of the present invention is provided.

In the invention, a testing environment can be built by using ThreeJS, and a volume drawing and anti-aliasing boundary smoothing flow of the run-length compressed three-dimensional geological lithology grid is realized by using WebGL 2.0 GLSL. Since the compressed three-dimensional grid may have unevenness, a Graciano method can be adopted to perform light transmission algorithm realization and accelerated optimization on data organization by using two-dimensional textures.

The test data is a geological three-dimensional grid model of engineering in a certain area, as shown in fig. 15, each formation property code on the right side in fig. 15 represents one formation property, such as a code 101, a code 102, a code 103, a code 104, a code 105, a code 106, a code 107, a code 108, a code 109, a code 110 and a code 201, and the total number of lithological layers is 11, and each formation layer is horizontally deposited and transversely continuous. Dividing the geological three-dimensional grid model according to the scale of 64 × 100 grid along the X axis, the Y axis and the Z axis, wherein the grid interval is 10 meters × 1 meter. In order to ensure the browsing effect in all directions, the grid can be stretched 10 times along the Z-axis direction.

And comparing the integral effect of the integral model volume rendering boundary smoothing treatment with the integral effect of the integral model volume rendering non-smoothing treatment, wherein the ray projection sampling of the boundary non-smoothing volume rendering algorithm is carried out according to the nearest neighbor attribute value, and the boundary smoothing volume rendering algorithm adopts the method of the invention to carry out attribute calculation. Comparative effects as shown in fig. 16, fig. 16a is a graph of the effects without smoothing, and fig. 16b is a graph of the effects after smoothing. Obviously, after the smoothing treatment, the sawteeth on the rock stratum interface disappear, and the drawing effect is improved. On the other hand, since the boundary smoothing process requires more texture reading operations and attribute calculation operations than the nearest neighbor sampling, the frame rate is obviously reduced, and the interactivity of the model is affected. For example, the frame rate after the non-smoothing process is 45 and the frame rate after the smoothing process is 30.

Comparing the rendering effect of the entire model rendering boundary smoothed with that of the non-smoothed longitudinal slices (YZ plane), as shown in fig. 17, fig. 17a is a graph of the non-smoothed effect and fig. 17b is a graph of the smoothed effect. Obviously, the sawteeth of the vertical slice rock layer interface disappear after the boundary smoothing treatment, and the drawing effect is improved. Meanwhile, the rendering frame rate after the boundary smoothing processing is also significantly decreased, for example, the frame rate after the boundary smoothing processing is 66, and the frame rate after the boundary smoothing processing is 42.

Comparing the rendering effect of the horizontal slices (XY-plane) without smoothing the rendering boundary of the entire model volume, as shown in fig. 18, fig. 18a is an effect diagram without smoothing and fig. 18b is an effect diagram after smoothing. It is clear that there is some improvement in the jaggy of the formation interface after the boundary smoothing process. Meanwhile, the rendering frame rate after the boundary smoothing processing is also significantly reduced, for example, the frame rate after the non-smoothing processing is 80, and the frame rate after the smoothing processing is 64. However, there are some unreasonable points in the boundary smoothing processing of the lateral slice, and specifically, for example, the a01 region of fig. 18a and the B01 region of fig. 18B are enlarged, and the enlargement effect is as shown in fig. 19. In fig. 19, the large serration pattern in fig. B01 is improved as compared to fig. a 01. However, as shown in the circle of fig. 19, there are still some problems. For example, a "burred" boundary still exists at the portion shown by circle g3, and unreasonable formation penetration still exists at the portions shown by circles g1 and g 2. These problems are related to the spatial overlap of local strata with each other caused by the hexahedral internal strata connection treatment method.

Specifically, the phenomenon in circle g1 is explained taking the phenomenon in circle g1 as an example. Referring to fig. 20, fig. 20 is a schematic diagram illustrating a principle of explaining formation interpenetrating according to another embodiment of the present invention, where fig. 20a is a schematic diagram of a slice of an XY plane, dots in fig. 20a are nodes of a formation property, and the slice of the XY plane includes a cell 1, a cell 2, and a cell 3, where the cell is a hexahedron. Wherein the dark parts represent the formation encoded as 105 and the light parts represent the formation encoded as 106. The formation 106 in cell 1 in fig. 20b is in the form of a spiral twist because the pinch-off occurs at the upper and lower points of the same vertical column, and there is some twist pinch-off of the two formations in cell 2, which causes the two formations to spatially intersect. The reason for this is that the formation connection rule only considers the space contact characteristics of the adjacent vertical columns, but does not consider the local and overall space characteristics of the rock formation, so that the two rock formations are not overlapped up and down after the cell is divided along the diagonal, as in the cell 1 of fig. 20b, the two adjacent sides 105 rock formations and 106 rock formations of the 105 rock formation vertical column are overlapped up and down, which is contradictory. According to the subdivision rule of the present invention, cell 1 and cell 2 can perform cell division and attribute calculation using the templates No. 2 and No. 13 in fig. 21, respectively, and the phenomenon shown in fig. 19 occurs. Therefore, if this phenomenon is eliminated, the mesh model can only be optimized to supplement the spatial sequence stacking relationship of the global rock strata and the connection relationship of the rock strata in the local cells.

In conclusion, the processing method of the three-dimensional geological lithology grid model provided by the invention improves the sawtooth problem of the volume drawing lithology boundary. According to the space sampling characteristics of the volume rendering ray projection algorithm, the surrounding attributes of sampling points are required to be used for calculation. The invention relates to a method for computing lithological attributes of any point in a model space by defining rock stratum connection rules on adjacent columns and 16 hexahedron subdivision templates, which aims to highlight rock stratum boundaries and ensure smoothness and continuity of the boundaries, longitudinally compressing a three-dimensional geological lithological grid by using run length codes to form vertical columns which are transversely and regularly distributed. The method is realized through WebGL, and the method is verified. The 16 hexahedral split templates can be shown in fig. 21, and fig. 21 totally comprises 16 templates of 0-15. After the step of obtaining the subdivision mode of dividing the hexahedron into tetrahedrons by combining the diagonals of the top surface and the bottom surface of the hexahedron, the method may further include: and performing binary coding according to the positive and negative slope of the object line of each side face in the hexahedron as the coding of the subdivision mode of the hexahedron. The object line can also be understood as a diagonal line. The 16 hexahedron subdivision templates sequentially form a side quadrangle anticlockwise according to the number of the vertical column in fig. 9, and binary coding is carried out according to the positive and negative slope of the object line of the quadrangle as the template number, and the corresponding code in fig. 9 is 0110 (the 6 th template in fig. 21). The slope binarization processing is convenient for template query, and meanwhile, the object line codes of the quadrangles are consistent in adjacent hexahedrons, so that the spatial consistency of the rock stratum in the processing of the transverse adjacent units is ensured.

By verification, the invention can obviously smoothen the rock stratum interface in the longitudinal direction, improve the sawtooth problem of the volume drawing lithologic boundary and improve the drawing effect. However, in the transverse direction, the problem of rock stratum boundary line sawtooth is improved to a certain extent, but because the reasonable space overlapping relation of mutual coverage between rock stratums in the subdivision units cannot be obtained globally, the rock stratum penetrating problem obviously and unreasonably exists at the rock stratum boundary, and the smoothing effect is influenced. Meanwhile, since the display frame rate of volume rendering is reduced by smoothing, real-time interactive rendering is not suitable.

The core of the invention is that 16 unit subdivision templates are designed for attribute calculation according to the space sampling characteristics and the continuous assumption of the engineering geological rock level, so that the condition of different attribute combinations can be uniformly processed, and the GPU parallel processing is convenient to perform.

Because the hexahedron where the sampling point is located mainly utilizes the local space connection rule (the rule of quadrangle triangulation) and the horizontal extension characteristic of rock stratum attributes, rock strata can be extinguished inside the hexahedron (unit), which is different from the rock stratum connection mode of the traditional Pillargrid grid, and the problem of tail warping caused by overlarge top-bottom position difference of adjacent rock strata can be avoided. Of course, the unit subdivision of the invention is more complex, and the subdivided fine areas are more, which is not suitable for the organization management of the grids. The subdivision method is suitable for the requirements of online visualization and rapid smoothing of the grid model, on one hand, hexahedron units do not need to be stored, only the position of a sampling point is needed to dynamically construct a hexahedron for attribute calculation, on the other hand, the volume rendering process based on the non-uniform grid does not need to be used, but the light projection volume rendering algorithm process based on the regular grid is completely utilized.

Because horizontal continuous assumption is carried out aiming at the characteristics of engineering geology, and clear sequence relations do not exist among rock stratums, the invention cannot process a fault model and cannot process the condition that the rock stratums with large inclination angles are not communicated transversely. On one hand, the subsequent research can consider increasing rock stratum connection constraint and improving the horizontal rock stratum boundary smoothing effect, and on the other hand, the method can be applied to lithologic three-dimensional grid multi-resolution drawing.

In some other embodiments, for any one side face, if the rock formation properties of two rock intervals in the side face are the same, and the two rock intervals have an intersection in the Z-axis direction, the start points of the two rock intervals may be connected, and the end points of the two rock intervals may be connected, so as to implement pinch-out of the rock formation. Specifically, referring to fig. 22, the formation properties of the B1 segment and the B2 segment are the same, and may both be denoted as B properties, and there is an intersection in the Z-axis direction, and the start point of the B1 segment and the start point of the B2 segment may be connected, and the end point of the B1 segment and the end point of the B2 segment may be connected, so as to implement rock pinch-out. In addition, when the side face of the present embodiment is triangularly split, because the rock formation properties of two rock intervals in the side face are the same, and the two rock intervals have an intersection in the Z-axis direction, the start point of any one of the two rock intervals may be connected to the end point of the other rock interval, and taking this diagram as an example, the start point of B1 may be connected to the end point of B2, or the end point of B1 may be connected to the start point of B2.

Referring to fig. 23, fig. 23 is a lithology effect diagram processed by a boundary linear smoothing method of three-dimensional geological lithology mesh model volume rendering according to another embodiment of the present invention. As can be seen from the figure, the problems of sawtooth and steps of the rock stratum interface are obviously improved, and the drawing effect of the model is improved.

Based on the same concept, the invention also provides a boundary linear smoothing device for the three-dimensional geological lithology grid model volume rendering, which comprises a memory, a processor and a machine executable program stored in the memory and operated on the processor, wherein when the processor executes the machine executable program, the boundary linear smoothing device realizes the boundary linear smoothing method for the three-dimensional geological lithology grid model volume rendering according to any one of the embodiments.

The above embodiments can be combined at will, and according to any one of the above preferred embodiments or a combination of multiple preferred embodiments, the embodiments of the present invention can achieve the following beneficial effects:

in the boundary linear smoothing method for drawing the three-dimensional geological lithology grid model, grids which represent the same rock stratum attribute in the Z-axis direction in the three-dimensional geological lithology grid model are combined by using run length codes to obtain a plurality of rock stratum sections; taking each rock stratum section in the Z-axis direction as a vertical column, and connecting the starting points of the rock stratum sections with space connection relation in each vertical column by using straight lines; the rock stratum segments representing the same rock stratum in four adjacent vertical columns are encircled into a hexahedron, and a subdivision mode of dividing the hexahedron into tetrahedrons is determined according to the attribute relationship of two rock stratum segments in each side of the hexahedron along the X-axis or Y-axis direction; manufacturing a calculation template of the rock stratum attribute of the sampling point in the hexahedron according to the subdivision mode; the rock stratum attributes of the sampling points are obtained through the calculation template, so that the rock stratum attributes of all positions in the hexahedron are determined, the problems of saw teeth and steps of a rock stratum interface are obviously solved, the linear smoothness of the boundary drawn by the three-dimensional geological lithology grid model is realized, and the drawing effect of the model is improved.

Thus, it should be appreciated by those skilled in the art that while a number of exemplary embodiments of the invention have been illustrated and described in detail herein, many other variations or modifications consistent with the principles of the invention may be directly determined or derived from the disclosure of the present invention without departing from the spirit and scope of the invention. Accordingly, the scope of the invention should be understood and interpreted to cover all such other variations or modifications.

Claims

1. A method of linear smoothing of boundaries drawn by a three-dimensional geological lithology mesh model volume, the three-dimensional geological lithology mesh comprising meshes distributed along an X-axis, a Y-axis and a Z-axis of a coordinate system, the method comprising:

merging grids which represent the same rock stratum attribute in the Z-axis direction in the three-dimensional geological lithology grid model by using run length codes to obtain a plurality of rock stratum sections;

obtaining rock stratum attributes of the sampling points through the calculation template, and accordingly determining rock stratum attributes of all positions in the hexahedron;

the step of determining the segmentation mode of the side face according to the attribute relationship of the rock stratum sections in the vertical column to which the two rock stratum sections in the side face of the hexahedron belong comprises the following steps:

if the target rock interval and the lower rock interval are overlapped in the Z-axis direction and the rock stratum properties are the same, connecting the starting point of the target rock interval with the end point of another rock interval;

if the target rock stratum section, the upper rock stratum section and the lower rock stratum section do not meet the condition that the target rock stratum section, the upper rock stratum section and the lower rock stratum section are overlapped in the Z-axis direction and have the same rock stratum property, judging the height of the tail points of the two rock stratum sections in the side face in the Z-axis direction, selecting the tail point higher in the Z-axis direction, and connecting the tail point with the starting point of the other rock stratum section to form a cutting line; or searching the higher end point in the Z-axis direction in the end points of the two rock stratum sections in the side face, and judging whether the end point of the rock stratum section connected with the higher end point is still higher than the other lower end point in the side face; if so, selecting the end point of the rock stratum section connected with the higher end point, and connecting the end point of the rock stratum section with the start point of another rock stratum section with the lower end point in the side face to form a cutting line;

and if the heights of the tail points of the two rock stratum sections in the side surface in the Z-axis direction are the same, calculating the lengths of the two diagonal lines of the side surface, comparing the lengths of the two diagonal lines, and selecting the longer diagonal line as the cutting line.

2. The method for linear smoothing of boundaries rendered by a three-dimensional geological lithology grid model according to claim 1, wherein,

after the step of obtaining a splitting manner of dividing the hexahedron into tetrahedrons by combining the diagonals of the top surface and the bottom surface of the hexahedron, the method further includes:

3. The method for linear smoothing of boundaries rendered by a three-dimensional geological lithology grid model according to claim 1, wherein,

the step of manufacturing a calculation template of the rock stratum attribute of the sampling point in the hexahedron according to the subdivision mode comprises the following steps of:

4. The method for linear smoothing of boundaries rendered by a three-dimensional geological lithology grid model according to claim 3, wherein,

the step of obtaining the rock stratum attribute of the sampling point through the calculation template comprises the following steps:

5. The method for linear smoothing of boundaries rendered by a three-dimensional geological lithology grid model according to claim 1, wherein,

after the step of obtaining the rock stratum attribute of the sampling point through the calculation template so as to determine the rock stratum attribute of each position in the hexahedron, the method further comprises the following steps:

if the vertical column is present, determining the position of the starting point or the end point of the vertical column of the hollow rock layer, taking the elevation of the hexahedral hollow rock layer as the elevation of the starting point or the end point of the vertical column of the hollow rock layer according to the position of the starting point or the end point of the vertical column of the hollow rock layer, and correcting the elevation of the rock layer with the lithology being non-hollow as the elevation of the starting point or the end point of the vertical column of the rock layer; and/or

6. An apparatus for linear smoothing of boundaries of a three-dimensional geological lithologic mesh model volume rendering comprising a memory, a processor and a machine executable program stored on the memory and run on the processor, and the processor when executing the machine executable program implements the method for linear smoothing of boundaries of a three-dimensional geological lithologic mesh model volume rendering according to any one of claims 1 to 5.