CN102867330B - Spatial Complex Horizontal Reconstruction Method Based on Region Division - Google Patents

Abstract

The invention discloses a region-division-based spatial complex horizon reconstruction method. The horizon reconstruction in three-dimensional geological modeling under the complex topography conditions such as multi-fault overlapping can be realized without the needing of a fault polygon; simultaneously, the dynamic control of specified grid spacing is realized by a hybrid grid interpolation method, so that a geologic structure can be expressed finely; the reconstructed geologic horizon surface and the fault are closed strictly; and the processing, such as blocking processing, in later geological modeling is facilitated.

Description

Spatial Complex Horizontal Reconstruction Method Based on Region Division

technical field

The invention belongs to the technical field of geological modeling, in particular to a method for reconstructing complex spatial horizons based on region division.

Background technique

Geological modeling is to process limited geological data through certain technical means on the computer to reproduce the geological structure model. The essence is to transform the limited discrete space sample point data into a continuous and visible geological surface or geological body, in which layer reconstruction is the key. The most important methods of level reconstruction are the interpolation method and the model structure of the level.

Interpolation is to use the method of surface fitting to establish a continuous function through discrete input sampling points, and use this reconstructed function to find the function value at any position. Horizontal reconstruction is a method of discretely reconstructing surfaces in a limited three-dimensional space.

The commonly used spatial interpolation methods are as follows:

(1) Reciprocal distance multiplication method: The reciprocal distance multiplication method is a weighted average interpolation method. It believes that several known points with the closest distance to the interpolation point contribute the most to the interpolation point, and its contribution is inversely proportional to the distance. The power parameter controls how the weights fall off with distance from a grid point. Closer data points are given higher weights, and farther data points are given less weight. As follows:

$\mathrm{<span\; class="notranslate">\; z</span>}<span\; class="notranslate">\; =</span>\frac{\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; p</span>}}}{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; p</span>}}}}$

In the formula, z is the estimated value, z _i is the value of the i-th known point, d _i is the distance between the i-th known point and the point to be interpolated, and p is the weighting coefficient, the higher the value, the more the influence factor of the distance The larger the value, the smoother the interpolation result will be. The inverse distance weighting algorithm has a faster calculation speed, but it is easily affected by the cluster of data points. In practical applications, an isolated point is often higher than its surrounding data points. Therefore, the inverse distance weighting method needs to be restricted according to different situations.

(2) Natural adjacency point interpolation method: The natural adjacency point interpolation method is a new gridding method unique to Surfer7.0. The basic principle is that for a group of Thiessen more deformations, when a new data point is added to the data set, these Thiessen polygons will be modified, and the weight average of the neighboring points will determine the weight to be interpolated. It is homogeneous and invariant, useful for interpolating evenly spaced data, but flawed for interpolating uneven data, which is common with raw horizon data.

(3) Kriging algorithm: also known as the optimal interpolation method of spatial auto-covariance, it is an optimal interpolation method named after the South African mining engineer D.G.Krige, and is a very useful geological statistics gridding method. This method is based on variogram theory and structural analysis, and calculates the weight by introducing a variogram with distance as an independent variable, which can not only reflect the spatial structure characteristics of variables, but also reflect the random distribution characteristics of variables. In the statistical sense, it is a method of unbiased and optimal estimation of the value of regionalized variables in a limited area from the perspective of variable correlation and variability; from the perspective of interpolation, it is a method for spatial A method for linear optimal, unbiased interpolation estimation of distributed data. Since Kriging algorithm has the best effect in spatial horizon interpolation, this scheme is based on Kriging algorithm to solve the problem of horizon interpolation in complex terrain.

The above interpolation methods are all limited to solving the level fitting of the two-dimensional plane, and cannot solve the level fitting problem of the three-dimensional space. In the actual geological structure, the geological situation is extremely complex, and the usual two-dimensional plane level fitting cannot solve the fitting problem of this kind of geological level. Therefore, it is necessary to transform the three-dimensional space problem into a two-dimensional space problem, and then use the above interpolation method for level fitting.

The structure of the layer model is the organization method of layer data, and different organization methods have different characteristics for model establishment. The main level structures in the existing proposed solutions are contour lines, regular grid (GRID), and irregular triangular network (TIN). Among them, regular grid and TIN are widely used layered model structures:

(1) Regular grid (GRID)

The regular grid is a regularly arranged rectangular grid to represent the horizon model. The data structure is simple, easy to construct a network, the data storage capacity is small, and can also be compressed and stored. Various analysis and calculations are very convenient and effective. However, due to the fixed grid spacing, there is a disadvantage that it is difficult to determine a realistic representation of the appropriate grid size for complex terrain.

(2) Irregular grid (TIN)

A TIN is formed by connecting the collected spatial-level sampling points according to certain rules into many triangles that cover the entire area and do not overlap each other. TIN can express the structure and details of the terrain relatively accurately, take into account the feature points and feature lines of the layer, and represent complex terrain more accurately than the rectangular grid. But the data storage and operation of TIN is complicated. The topological relationship of triangular points and sides of TIN is relatively complex, it is difficult to perform terrain analysis on TIN, and the storage space is large, so it is generally only suitable for high-precision terrain modeling of small scale and large scale.

The prior art relevant to the present invention includes:

Regarding the reconstruction of the 3D geological level, there are quite a few implementation schemes at home and abroad. The existing schemes are divided into two types according to whether fault polygons are needed: one is the slice reconstruction scheme based on fault polygon constraints, and the other is the slice reconstruction scheme based on regional division.

(1) Layer reconstruction scheme based on fault polygon constraints

The fault polygon is the intersection of the three-dimensional fault plane and the horizon plane. The level reconstruction work needs to comprehensively analyze the fault polygons on multiple interfaces, and infer the shape of the three-dimensional spatial horizon according to the constraints of the fault polygons. Many interpretation system algorithms are also completed on the basis of each cross-sectional image, rather than directly reconstructed on the 3D original data. At present, the manual tracing method is commonly used in fault polygon combination, that is, the manual breakpoint connection method is used. This method not only has large errors and affects the accuracy of the constructed map, but also has extremely low efficiency. Wang Zhaohu, Liu Cai, etc. proposed an automatic fault polygon combination technology based on horizon interpretation, but this technology can only automatically track fault polygons when the horizon is very continuous and the data is relatively fine, and the data requirements are relatively high, so it is not applicable general situation.

(2) Level reconstruction scheme based on regional division

There are few existing schemes for 3D geological level reconstruction based on regional division. At present, only Cai Qiang and Yang Qin proposed a meshing method for overlapping subregions ^[6] . This algorithm introduces the concepts of bridge edges and intersecting ring lines, and realizes layer reconstruction in overlapping areas through a linkage subdivision algorithm. This scheme has disadvantages: First, this scheme divides overlapping sub-regions by finding intersection lines in 3D space, and the solution of intersection lines in 3D space is more complicated and slower. 2. This plan does not propose a method to solve the three-dimensional geological layer reconstruction in special cases such as multiple reverse faults, layer pinch-out, and layers rushing out of the surface and underground.

The existing schemes are divided into two categories according to the organizational structure of the reconstruction layer: the layer reconstruction based on the regular grid (GRID) and the layer reconstruction based on the triangle mesh (TIN). The regular grid is a regularly arranged rectangular grid to represent the horizon model. Due to the fixed grid spacing, it is difficult to determine a realistic representation of the appropriate grid size for complex terrain. The topological relationship of the triangle points and edges of the triangular mesh is relatively complex, it is difficult to perform terrain analysis on TIN, and the storage space is large, so it is generally only suitable for small-scale and large-scale high-precision terrain modeling.

At present, there are many schemes for geological 3D reconstruction and interpolation at home and abroad, but they all have the following limitations:

(1) Most of the current three-dimensional geological level reconstruction schemes need to obtain fault multivariate lines in advance. Manually editing fault polygons is a huge workload, and the accuracy is difficult to guarantee and the efficiency is low. However, there is no good solution for automatically generating fault polygons.

(2) Most geological three-dimensional reconstruction schemes do not consider faults. Although some schemes consider faults, they do not consider multiple reverse faults.

(3) The grid spacing of existing three-dimensional geological level reconstruction schemes based on regular grid processing is definite, which cannot dynamically adapt to the needs of geological structures, and it is difficult to express geological structures finely.

Contents of the invention

In order to overcome the above-mentioned shortcomings of the prior art, the present invention provides a spatially complex horizon reconstruction method based on region division. By dividing overlapping subregions, complex terrain such as multiple fault overlaps can be solved without the need for fault polygons. The level reconstruction problem in 3D geological modeling under certain conditions; at the same time, through the interpolation method of mixed grid, the dynamic control of the standard grid spacing is realized, so that the geological structure can be finely expressed, and the reconstructed geological layer plane and faults are strictly closed, which is convenient for later processing in geological modeling, such as block processing.

The technical solution adopted by the present invention to solve the technical problem is: a method for reconstructing complex spatial horizons based on regional division, comprising the following steps:

Step 1. Hybrid grid interpolation method: nest small grids inside the large grid, interpolate the large grid first, and then use the existing seed points of the large grid and the new interpolation on the large grid Point to small grid interpolation;

Step 2, defining the layer segment data structure;

Step 3. Fault level fitting: take the two-dimensional envelope of the interpretation data of the fault seed segment on the grid plane as the interpolation envelope range of the section fitting, and then use the mixed grid interpolation method to analyze the fault network within the fault envelope range. The grid points are interpolated, and the interpretation data is specifically, defining the data structure of the segment interpretation data:

Seg＝{ID,Direct,UpOrDown,Start,End,Points,LeftInfo,RightInfo},

Among them, ID is the index number of the layer to which the segment belongs,

Direct is the segment direction: the segment in the X-axis direction is X_Direct, and the segment in the Y-axis direction is Y_Direct.

UpOrDown is the relationship between a segment and a fault: Up is located on the fault hanging wall, Down is located on the fault footwall,

Start is the segment start position,

End is the segment end position,

LeftInfo is the associated information Info of the left header of the segment,

RightInfo is the associated information Info at the right end of the segment,

The section header and the section tail associated information Info structure of the section are: Info={ID, UpOrDown},

Points is a set of discrete points of the segment;

Step 4. Horizontal original data editing: delete wrong data; define horizon boundaries and faults as constraint boundaries. If there are multiple horizon segments between two constraint boundaries on a two-dimensional section, these The horizon segments are combined into one horizon segment;

Step 5. Initialize the original segment data: initialize the relationship between the segment and the constraint boundary;

Step six, sub-region division;

Step 7. Fitting of fault sub-region layer planes:

(1) Initialize the upper and lower disk seed points;

(2) Surface interpolation of the plate area above and below the fault;

(3) The interception of the actual disk surface above and below the fault;

Step 8. Horizontal plane fitting of non-fault sub-regions:

(1) Initialize the seed points: take the data points of the non-fault area as the seed points, and the data points of the boundary area of the fault sub-area as the seed points; and for the interpolated fault sub-area, all the Data points without heavy values are also used as seed points;

(2) Interpolation of non-faulted sub-regions: for the generated seed points, the hybrid grid interpolation method is used to interpolate each grid point of the non-faulted sub-regions;

Step 9: Processing of layer pinch-out and layer rushing out of the surface and underground: remove the illegal pinch-out part of the current layer and the part rushing out of the surface and underground.

Compared with the prior art, the positive effect of the present invention is that it solves the problem of stratum surface reconstruction under complex terrain in three-dimensional geological modeling, and provides new ideas for contour drawing and geological block formation. The inventive method has the following advantages:

(1) Only the original layer sampling data and fault sampling data are needed, and the layer can be reconstructed without fault polygons. Fault polygons that do not need to be manually edited, eliminating the need for cumbersome manual editing.

(2) It supports various types of faults, including normal faults and reverse faults, and can handle multiple faults at the same time. It is suitable for various complex terrain structures and can handle horizon interpolation in various complex situations. It has good adaptability .

(3) It is suitable for various applications in 3D geological modeling and 2D layer rendering, and has good versatility.

(4) The horizontal interpolation method of hybrid gridding is adopted, which not only ensures the rapid display of data, but also ensures the accuracy of fitting.

(5) Support the processing of layer pinch-out and layer rushing out of the surface and underground.

(6) A method of subdividing sub-regions by section envelope is proposed, and a specific implementation is given.

(7) Through boundary constraints, the seamless fitting of layers and faults is realized. The reconstructed geological layers and fault planes are strictly closed.

Description of drawings

The invention will be illustrated by way of example with reference to the accompanying drawings, in which:

Fig. 1 is a flow chart of the inventive method;

Fig. 2 is a schematic diagram of deleting erroneous data;

Fig. 3 is a schematic diagram of initializing original segment data.

Detailed ways

First define some basic geological structures and scheme terms:

Horizon: Refers to a specific position in a stratigraphic sequence.

Fault: After the rock mass is fractured by force, the rock blocks on both sides have significant displacement along the fracture surface.

Segment data: A data set consisting of a continuous set of discrete points on a section of a horizon or fault.

Hanging wall of a fault: Located above the fault plane.

Fault footwall: located below the fault plane.

Stratigraphic pinch-out: refers to the sedimentary layer towards the edge of the sedimentary basin, and its thickness gradually becomes thinner until there is no deposition.

Gridding: Logically divide the discrete point data to form a regular logical grid, which is convenient for horizon interpolation.

Interpolation: The process of using known points to calculate unknown points.

Fitting: A process of using the data after horizon interpolation to form horizons.

A method for reconstructing complex spatial horizons based on regional division, as shown in Figure 1, includes the following steps:

Step 1. Hybrid grid interpolation method

In the existing commonly used grid interpolation methods, if the grid is too sparse, the accuracy will be insufficient, and if the grid is too dense, the interpolation speed will be affected. This scheme proposes a hybrid grid interpolation method with two-level grid nesting. The first-level grid ensures the speed of interpolation, while the second-level grid ensures the accuracy of interpolation.

The hybrid grid is a secondary grid, that is, a small grid is nested inside a large grid. When interpolating, first interpolate the points on the large grid, and then use the existing seed points of the large grid and new interpolation points on the large grid to interpolate the small grid. Due to the ratio of the size grid, the range of finding the seed point and the number of seed points can be manually input, so that the interpolation can be set differently for different situations. This not only solves the problem of insufficient interpolation accuracy due to too few seeds supported by direct Kriging, but also improves the interpolation speed.

(1) Large grid interpolation

For the large grid point to be interpolated, the interpolation calculation is performed on the interpolation point by searching the legal seed points around the interpolation point. In some horizons where the seed points are sparse, it is impossible to calculate all the interpolation points only by relying on the seed points for interpolation. In this case, keep the interpolation points that cannot calculate the Z value, and after all the interpolation points that can rely on the seed point for interpolation calculation are processed, use the interpolated point as the new seed point for horizon compensation interpolation until All interpolation points are processed. The specific method of large grid interpolation is as follows:

1) Scan large grid points one by one to determine whether interpolation is required for large grid points. If there are already data points in this large grid point, interpolation is not required, otherwise interpolation is required.

2) For the large grid points that need to be interpolated, according to the manually input search range, that is, the number of large grids to be searched again, search for the seed point data of the large grids around the current interpolation point.

3) Sector classification is performed on the seed points according to the manually input sectors.

4) Arrange the seed points found in each sector according to their distance from the current large grid point from small to large.

5) Select the threshold value according to the number of seed points in each sector input, and select the seed point data within the threshold range closest to the grid point to be interpolated within the sector range and add it to the current large grid interpolation seed point queue.

6) If the current large grid point does not find effective seed point data for interpolation, the current large grid point does not interpolate, and turns to the next large grid point, otherwise, it turns to step 7).

7) Use the seed points in the current large grid interpolation seed point queue to perform Kriging interpolation on the current large grid points.

8) If all the large grid points have been processed once, check whether there is a large grid point interpolation failure, if there is a large grid point interpolation failure, then use the large grid point that has been successfully interpolated as the seed point, and return Step 1) Perform re-interpolation until all grid points are successfully interpolated.

(2) Small grid point interpolation

After large grid interpolation, all large grid points on the horizon grid have been interpolated with corresponding data. Then use the existing large grid point data to perform subdivision and interpolation of small grids. In order to improve the interpolation speed, the interpolation of small grid points is processed in units of large grids. For a large grid, search the seed points in the current large grid and the seed point data in the eight surrounding large grids as the original seed points, and then use Kriging interpolation to fit these original seed points into a surface.

All small grid points in the large grid can be interpolated directly by taking the value on this surface, avoiding multiple Kriging interpolation. That is, through one piece of data and one Kriging interpolation fitting, the interpolation of small grid points in the entire large grid is realized, which greatly improves the efficiency of interpolation. For a certain interpolation point I _x , if the legal seed point set found is S'={S ₁ , S ₂ ,...,S _k }, and k≥2 is satisfied, kriging can be used for interpolation. The coordinates of the point set corresponding to S' are SP={p ₁ ,p ₂ ,…,p _k }, where p _i (1≤i≤k) is the three-dimensional coordinate point (x _i , y _i , z _i ), I _x The coordinates of are (x ₀ ,y ₀ ,z), and z is unknown. Kriging horizon interpolation is the process of using the known point set SP and the xy coordinates (x ₀ , y ₀ ) of the interpolation point to calculate the z value of the interpolation point.

Step 2. Define the layer segment data structure

Horizontal segment data is a collection of continuous discrete points. Since a segment is a collection of discrete points, the attributes of a segment are uniform for all seed points on the segment, and the relationship with faults is logical.

In this solution, the data structure of the segment interpretation data is defined:

Seg＝{ID,Direct,UpOrDown,Start,End,Points,LeftInfo,RightInfo}

ID: The index number of the layer to which the segment belongs.

Direct: segment direction: the segment in the X-axis direction is X_Direct, and the segment in the Y-axis direction is Y_Direct.

UpOrDown: The relationship between the segment and the fault: Up is located on the fault hanging wall, and Down is located on the fault footwall.

Start: segment start position.

End: The end position of the segment.

The associated information structure of the header and tail of the definition section is: Info＝{ID,UpOrDown}

LeftInfo: the associated information Info of the left section header of the section.

RightInfo: the related information Info of the right end of the segment.

Points: A collection of discrete points for the segment.

The horizon set is defined as: S={S ₁ , S ₂ . . . S _i }, where S _i is the i-th horizon set. The data structure of layer S _i is S _i ={ID,XSegments,YSegments}, ID is the layer number of layer S _i , XSegments is the data set of all segments belonging to the X-axis direction of layer S _i , and YSegments is Y A collection of all segment data in the axis direction.

Step 3. Fault plane fitting

A fault is a structure in which the crustal rock layer is ruptured due to the force reaching a certain strength, and there is obvious relative movement along the rupture surface. The staggered planes of fault fractures are called fault planes. The rock blocks on both sides of the fault plane are called the fault plate, the rock blocks above the fault plane are called the hanging wall, and the rock blocks below the fault plane are called the footwall. According to the unique nature of the faults, they can be divided into normal faults with relatively lower hanging wall and reverse faults with relatively higher hanging wall.

Section fitting needs to determine the interpolation range of section fitting. Since this scheme is based on the grid, the two-dimensional envelope of the interpretation data of the fault seed segment on the grid plane is taken as the range of the interpolation envelope for section fitting. Then, the interpolation method of mixed grid is used to interpolate the fault grid points within the fault envelope.

Step 4: Layer raw data editing

For some horizon data, the original data may partly violate the actual geological structure, mainly as follows: As shown in Figure 2, there is and can only be one horizon plane in the hanging wall or footwall of a fault , and when interpreting seismic data, erroneous data may appear. Therefore, it is necessary to automatically crop and manually edit the original data of the horizon, delete the wrong data for interpolation of the horizon space, and fit the correct horizon.

If horizon boundaries and faults are defined as constraint boundaries, then all horizon segments between the two constraint boundaries have the same attributes. If there are multiple horizon segments between two constraint boundaries on a 2D section, these horizon segments need to be combined into one horizon segment. The specific combination method is as follows:

(1) Sort the horizon segments on the current profile in ascending order according to the start position of the segment.

(2) Count the faults on the current profile.

(3) Process each fault. Hanging wall processing: Find the horizon segment Seg adjacent to the hanging wall of the current fault, and then use this horizon segment as the initial segment to find the next horizon segment that can be merged with it for merging. Seg update points to the merged segment, recursively merges the next segment until a certain constraint boundary; footwall processing: find the horizon segment Seg adjacent to the current fault footwall, and then use this horizon segment as the initial segment to find the next one that can Merge with the merged layer segment, Seg update points to the merged segment, and recursively merge the next segment until a certain constraint boundary.

(4) If there is no fault, all the horizon segments only need to be merged into one segment according to the starting and ending positions.

Step 5. Initialize the original segment data

Since the horizon segment is constrained by the constraint boundaries defined by horizon boundaries and faults, the segment head and segment tail must be associated with a certain constraint boundary. Initializing the original segment data means initializing the relationship between the segment and the constraint boundary, as shown in Figure 3.

Step 6. Sub-area division

The purpose of sub-region division is to divide the entire horizon into several small sub-regions, and these sub-regions have no weight value, so that the interpolation of more complicated heavy faults and reverse faults becomes interpolation without weight value. The division of sub-regions is determined according to the envelope range of the section, and each fault corresponds to two sub-regions, namely the upper wall sub-region and the lower wall sub-region.

Define the envelope range of the section Fi as Pi, and Pi is the envelope range of the two-dimensional projection of the section Fi on the grid. Define the upper plate area corresponding to the fault Fi as Ui, and the lower plate area as Di. Then Ui and Di are the current fault envelope and the set of a series of fault envelopes associated with the current fault Fi.

The specific method of sub-region division is as follows:

(1) If the envelope Pi of the fault Fi does not overlap with the envelopes of other faults, the upper and lower plates of the fault Fi are the envelope Si of the fault Fi, that is, Di={Pi}, Ui={Pi}.

(2) If the envelopes {Ps...Pt} of multiple faults {Fs...Ft} overlap with the envelopes of fault Fi and are directly adjacent to the upper (lower) wall of fault Fi, then The upper (lower) plate area of the fault Fi is the set of envelopes of all related faults, ie Ui(Di)={Pi,Ps...Pt}.

(3) After the upper and lower plate regions of all faults are divided, the remaining non-fault horizon region is regarded as a single non-fault sub-region LH.

Then the entire area of horizon H can be expressed as: H={U1, D1...Ui, Di, LH}. Ui and Di are the upper and lower wall regions of fault Fi associated with horizon H, respectively, and LH is the non-fault region.

Step 7: Fitting of fault sub-region horizon planes

(1) Initialize the upper and lower disk seed points

To interpolate a certain sub-area, first determine the seed point data for this sub-area interpolation. These seed point data have no heavy value on the two-dimensional plane projection.

Let the data segment structure of the i-th layer be S _i ={ID,XSegments,YSegments}, where ID is the layer number of the layer, XSegments is the data set of all segments in the X-axis direction, and YSegments is all the segments in the Y-axis direction collection of data. Define UpSegs as the seed segment of the hanging wall of the fault, and DownSegs as the seed segment of the footwall of the fault. Then the collection form of UpSegs or DownSegs is

UpSegs={seg1,seg2...segi}

DownSegs={seg1,seg2...segj}

Among them, the segment segi or segj is a certain segment in the segment set of horizon S _i , and has the following attributes: the segment head of segi or segj is associated with the hanging wall of the current fault, and the segment tail of segi or segj is associated with the footwall of the current fault . Let the set of seed points in the upper plate area of the fault be UpPoints, and the set of seed points in the lower plate area be DownPoints. Then the expressions of UpPoints and DownPoints are:

UpPoints={p1,p2...pi}

DownPoints={p1,p2...pj}

In the above formula, pi and pj are points in the Ui and Di segments, respectively, and the two-dimensional projections of the two on the grid coordinates are respectively in the upper and lower subregions (UpPoints and DownPoints) of the fault.

(2) Surface interpolation of the upper and lower plates of the fault

Using the interpolation method of mixed grids, taking UpPoints and DownPoints as interpolation seed points respectively, interpolation is performed on each grid point in the range of the upper and lower plates of the fault.

(3) Interception of the actual disk surface above and below the fault

The layer plane data of the upper and lower plate areas of the fault after interpolation is defined as the initial upper and lower plate areas, and the initial upper and lower plate data are within the range of the entire fault upper and lower plate areas, which are larger than the actual upper and lower plate plane data, passing through the fault plane . It is necessary to intercept the initial upper and lower disk data according to the constraints of the fault. In order to strictly seal the upper and lower walls and faults, we seal the actual lower walls at the same time as the section.

If the line connecting two points does not cross the fault, the two points are said to be connected. We take a point that must be on the actual upper and lower disk data as the initial source point, and recursively search for the actual upper and lower disks through connectivity. In fact, the recursive search method of the disk is as follows:

1) Traversing the grid of the initial upper wall to find a point p that must be on the grid of the fault plane, then the point p must be on the actual upper wall, so point p is taken as the initial source point.

2) Set up a data stack STACK to store the source point. First push the source point p into the data stack STACK.

3) Judging whether the data station STACK is empty, if it is empty, the recursive search ends. Otherwise, take out a little q from the data stack STACK. Traverse the four points pi (i=1, 2, 3, 4) of the point q on the grid coordinates.

4) If pi can be connected with q, then add pi to the actual disk, and push pi into the data stack STACK at the same time. If pi cannot be connected with q, it means that q is actually the boundary point of the disk surface, then change the z value of point q to the z value of the fault plane corresponding to the grid point.

The same is true for the recursive search of the actual lower disk and the processing of the upper disk.

Step 8. Horizontal plane fitting of non-fault sub-regions

(1) Initialize the seed point

If only the data points located in the non-fault area are selected as seed points for horizon interpolation fitting, the transition between the horizon plane of the non-fault sub-area and the horizon plane of the fault sub-area will not be smooth. The invention not only uses the data points of the non-fault region as the seed points, but also uses the data points of the boundary regions of the fault sub-regions as the seed points, thereby solving the problem of unsmooth transition between the sub-regions. For the interpolated fault sub-region, all the data points without repeated values on the two-dimensional plane are also used as seed points.

(2) Interpolation of non-faulted sub-regions

For the seed points that have been generated, use the interpolation method of the mixed grid to interpolate each grid point in the non-fault sub-region. After the interpolation is completed, the surface fitting of the entire horizon is completed.

Step 9: Treatment of layer pinch-out and layer rushing out of the surface and underground

The situation of layer pinch-out and layer rushing out of the surface and the bottom is essentially the same, and both pinch-out layers and the surface and bottom are used to constrain the current layer. Therefore, it is only necessary to remove the illegal pinch-out part of the current layer and the part rushing out of the surface and underground.

Define the current interpolation horizon as the main horizon S, the pinch-out horizon with the current horizon as H _i , the ground surface as UH, and the underground as DH. Defining H _i , UH, and DH as constrained horizons, there is a relative relationship between constrained horizons and main horizons, that is, which horizon is above and which is below. Then, the treatment methods for layer pinch-out and layer rushing out of the surface and underground are as follows:

(1) If the current constrained horizon is the surface, intercept all points that go out of the surface (that is, the Z value on the grid point is smaller than the Z value of the corresponding grid point on the ground surface).

(2) If the current constrained horizon is underground, intercept all the points that go out of the underground (that is, the Z value on the grid point is greater than the Z value of the corresponding grid point on the underground surface).

(3) If the current constrained horizon is a pinch-out horizon. If the pinch-out horizon is on the main horizon S, intercept the point whose Z value on the grid point is smaller than the Z value of the grid point corresponding to the constraint horizon. If the pinch-out horizon is under the main horizon S, intercept the point whose Z value on the grid point is greater than the Z value of the grid point corresponding to the constraint horizon.

Claims (3)

1. A method for spatially complex horizon reconstruction based on regional division, characterized in that: comprise the steps: Step 1. Hybrid grid interpolation method: nest small grids inside the large grid, interpolate the large grid first, and then use the existing seed points of the large grid and the new interpolation on the large grid Points are interpolated on a small grid, The method for interpolating a large grid is as follows: 1) Scan the large grid points one by one to determine whether interpolation is required for the large grid points. If there are already data points in the large grid points, no interpolation is required, otherwise interpolation is required; 2) For the large grid points that need to be interpolated, according to the manual input search range, search for the seed point data of the large grid around the current interpolation point; 3) Carry out sector classification to the seed point according to the manually input sector; 4) For the seed points found in each sector, arrange them from small to large according to their distance from the current large grid point; 5) Select the threshold according to the number of seed points in each sector input, and select the seed point data within the threshold range closest to the grid point to be interpolated in this sector to add to the seed point queue of the current large grid interpolation; 6) If the current large grid point does not find effective seed point data for interpolation, the current large grid does not interpolate, and turns to the next large grid point, otherwise it turns to step 7); 7) Carry out Kriging interpolation to the current large grid point with the seed point in the current large grid interpolation seed point queue; 8) If all the large grid points have been processed once, check whether there is a large grid point interpolation failure, if there is a large grid point interpolation failure, then use the large grid point that has been successfully interpolated as the seed point, and return Step 1) performs re-interpolation until all grid points are interpolated successfully; Step 2, defining the layer segment data structure; Step 3. Fault level fitting: take the two-dimensional envelope of the interpretation data of the fault seed segment on the grid plane as the interpolation envelope range of the section fitting, and then use the mixed grid interpolation method to analyze the fault network within the fault envelope range. The grid points are interpolated, and the interpretation data is specifically, defining the data structure of the segment interpretation data: Seg＝{ID,Direct,UpOrDown,Start,End,Points,LeftInfo,RightInfo}, Among them, ID is the index number of the layer to which the segment belongs, Direct is the segment direction: the segment in the X-axis direction is X_Direct, and the segment in the Y-axis direction is Y_Direct. UpOrDown is the relationship between a segment and a fault: Up is located on the fault hanging wall, Down is located on the fault footwall, Start is the segment start position, End is the segment end position, LeftInfo is the associated information Info of the left header of the segment, RightInfo is the associated information Info at the right end of the segment, The section header and the section tail associated information Info structure of the section are: Info={ID, UpOrDown}, Points is a set of discrete points of the segment; Step 4. Horizontal original data editing: delete wrong data; define horizon boundaries and faults as constraint boundaries. If there are multiple horizon segments between two constraint boundaries on a two-dimensional section, these The horizon segments are combined into one horizon segment; Step 5. Initialize the original segment data: initialize the relationship between the segment and the constraint boundary; Step six, sub-region division; Step 7. Fitting of fault sub-region layer planes: (1) Initialize the upper and lower disk seed points; (2) Surface interpolation of the plate area above and below the fault; (3) The interception of the actual disk surface above and below the fault; Step 8. Horizontal plane fitting of non-fault sub-regions: (1) Initialize the seed points: take the data points of the non-fault area as the seed points, and the data points of the boundary area of the fault sub-area as the seed points; and for the interpolated fault sub-area, all the Data points without heavy values are also used as seed points, and the heavy values are specifically that the envelopes between faults overlap with each other; (2) Interpolation of non-faulted sub-regions: for the generated seed points, the hybrid grid interpolation method is used to interpolate each grid point of the non-faulted sub-regions; Step 9: Processing of layer pinch-out and layer rushing out of the surface and underground: remove the illegal pinch-out part of the current layer and the part rushing out of the surface and underground. 2. the spatially complex horizon reconstruction method based on regional division according to claim 1, is characterized in that: the method for combining a plurality of horizon segments into one horizon segment is: (1) Sort the horizon segments on the current profile in ascending order according to the start position of the segment; (2) Count the faults on the current profile; (3) Each fault is processed, and the described processing of each fault is specifically, hanging wall processing: find the horizon segment Seg adjacent to the hanging wall of the current fault, and then use this horizon segment as the initial segment to search for the next A horizon segment that can be merged with it is merged, the Seg is updated to point to the merged segment, and the next segment is merged recursively until a certain constraint boundary is reached. Footwall processing: find the horizon segment Seg adjacent to the footwall of the current fault, and then use This horizon segment is the initial segment, find the next horizon segment that can be merged with it and merge, Seg update points to the merged segment, recursively merge the next segment until a certain constraint boundary; (4) If there is no fault, it is only necessary to merge all horizon segments into one segment according to the starting and ending positions. 3. the spatial complex horizon reconstruction method based on region division according to claim 1, is characterized in that: the method for described subregion division is: (1) If the envelope of a certain fault does not overlap with the envelope of other faults, the upper and lower plates of the fault are the envelope of the fault; (2) If the envelope range of multiple faults overlaps with that of a certain fault and is directly adjacent to the hanging wall or footwall of the fault, the hanging wall or footwall subarea of the fault is defined as all relevant faults. The set of envelope ranges; (3) After the upper and lower plate areas of all faults are divided, the remaining non-faulted horizon area is regarded as a single non-faulted sub-area.