[0064] First define some basic geological structures and terms of this plan:
[0065] Horizon: refers to a specific position in the stratigraphic sequence. The stratum level can be the boundary of the stratigraphic unit, or it can be a marker layer belonging to a specific era.
[0066] Fault: A structure in which the crustal rock layer is broken due to a certain strength and has obvious relative movement along the fracture surface is called a fault.
[0067] Vertical fault: A fault with a small distance.
[0068] Gridization: The discrete point data is logically divided to form a regular logical grid, which is convenient for layer interpolation.
[0069] Interpolation: The process of using known points to calculate unknown points.
[0070] Fitting: A process of using the data after the level interpolation is completed to form the level.
[0071] Such as image 3 As shown, a surface reconstruction method in a complex space includes the following steps:
[0072] The first step, raw fault data preprocessing
[0073] From a geological point of view, a fault is a fracture structure in which rock blocks on both sides are significantly displaced along the fracture surface after the rock mass is broken by force. Therefore, the fault performance in the seismic interpretation system is generally a relatively steep curved structure. It is usually divided into: normal faults with a relatively descending hanging wall and reverse faults with a relatively rising hanging wall. In actual geological structures, some faults have very small distances. When we generate structural maps, we usually use a line to describe this fault. In computer processing, it can be described as a vertical fault, and its projection on the horizontal plane is a broken line. .
[0074] The preprocessing of the original fault data includes three aspects: fault interpolation, fitting fault planes and generating fault polygons related to the horizon based on the original three-dimensional data of the horizon. Because the fault has no heavy value in the depth domain (or time domain), the fault interpolation is relatively simple. You only need to determine a fault envelope range, you can use the kriging algorithm to directly interpolate, and then fit the section according to the interpolation data . Due to its special structure, vertical faults do not need to be fitted. A fault polygon is a polygonal structure where a fault and a horizon intersect spatially. Its boundary is composed of fault lines. The fault lines are divided into hanging wall lines and lower wall lines. The upper and lower wall lines form a closed polygon in a three-dimensional space, such as Figure 4 Shown.
[0075] A fault may have multiple fault polygons, but each fault polygon is associated with a unique horizon. The relationship between fault polygon and horizon is as Figure 5 As shown, here we define the fault set F, the polygon set Fi of the i-th fault, and the horizon set L, which have:
[0076] F={F 1 , F 2 , L F m }
[0077] F i = P 1 , P 2 , L P k i , 0 ≤ k i ≤ n Formula (2-1)
[0078] L={L 1 , L 2 , L L n }
[0079] Each fault polygon is composed of an upper inventory set and a lower inventory set (for vertical faults, the upper inventory and the lower inventory completely overlap):
[0080] P i ={up i , Lp i } Formula (2-2)
[0081] We use F i P j Representing the j-th polygon of the i-th fault, the set of fault polygons associated with the horizon is expressed as follows:
[0082] L i = F i 1 P j 1 , F i 2 P j 2 L F i a P j a , 0 ≤ a ≤ m Formula (2-3)
[0083] The second step, the original layer data preprocessing
[0084] For some horizon data, some of the original data may violate the actual geological structure, such as Image 6 As shown, on the upper or lower wall of the fault, there is one and only one horizon, and the original horizon data may not meet this condition, which requires preprocessing of the original horizon data to remove Illegal point data. The preprocessing of original horizon data is a necessary process to ensure the correctness of horizon interpolation.
[0085] The third step, layer interpolation
[0086] Horizon interpolation is a necessary process to fit the discrete points of stratum into layers. In this scheme, a stratum grid (Row×Col) is calculated in advance through the specified interpolation accuracy, and stratum interpolation is to calculate each grid The horizon point value at the point. In general, horizon interpolation is the process of using known points to calculate unknown points. Such as Figure 7 As shown, when interpolating point P, it is necessary to search for known points around to perform interpolation calculation. The search range is generally a circle with P as the center and R as the radius. The radius R can be a fraction of the work area, or it can be manually designated according to actual needs.
[0087] Due to the existence of various faults, horizon interpolation is restricted by faults, and not all the seed points found in the search range can be used for calculation. In this case, the traditional point-finding scheme is not suitable for the horizon under such complex conditions. Bit interpolation. Such as Figure 8 As shown, the definition requires interpolation point set I={I 1 , I 2 , L, I k , L, I max }, and call the known points the seed point, and define the seed point set S={S 1 , S 2 , L, S k , L, S N }. If we set the search radius R=8 (that is, search for eight nearby grids), then in the interpolation I 1 , The seed point set found is {S 1 , S 2 , S 3 , S 4 , S 5 , S 7 }, but only seed points {S 1 , S 2 , S 3 } Can be used to calculate I 1 The Z value of, which can be used to calculate the seed point of a certain interpolation point is called the legal point of the interpolation point.
[0088] Such as Picture 9 As shown, the process of horizon interpolation is as follows:
[0089] Initialize the associated fault polygons: get all the fault polygons associated with the horizon to be interpolated, initialize the constrained fault boundary set for horizon interpolation, and discretize the associated fault polygon boundary into three-dimensional control points according to the grid. These three-dimensional control points will be The boundary control during horizon interpolation plays an important role.
[0090] Initialize the grid properties: determine the grid parameters (number of rows, columns, etc.), project the associated fault polygon to the grid, and determine the attributes of each grid (the grid is outside the fault, boundary or inside).
[0091] Initialize seed point data: Obtain the original data of the horizon to be interpolated and the three-dimensional control points of the associated fault polygon boundary, convert them into a discrete three-dimensional seed point data structure, obtain the fault-related attributes of each seed point, and then convert the seed points according to (x , Y) coordinates are projected into each grid. Define the seed point structure as follows (where type is the seed point type, cp is the three-dimensional point coordinate, uF is the set of faults associated with the hanging wall, 1F is the set of faults associated with the bottom wall, and the faults in uF and 1F are from small to large according to the depth value. Sort):
[0092] S i ={type i , Cp i , UF i , LF i } Formula (2-4)
[0093] Such as Picture 10 Shown:
[0094] Initialize seed point S 1 Attributes: the type is horizon point data, the coordinates are (x, y, z), the hanging wall associated fault set is {null}, and the foot wall associated fault set is {null};
[0095] Initialize seed point S 2 Attributes of: the type is the upper fault inventory data, the coordinates are (x, y, z), the hanging wall associated fault set is {null}, and the foot wall associated fault set is {F1, F2};
[0096] Initialize seed point S 4 Attributes of: the type is horizon point data, the coordinates are (x, y, z), the hanging wall associated fault set is {F1}, and the foot wall associated fault set is {F2};
[0097] Initialize seed point S 6 Attributes of: Type is fault count data, coordinates are (x, y, z), hanging wall related fault set is {F1, F2}, bottom wall related fault set is {null};
[0098] Initialize seed point S 8 Attributes of: Type is upper fault inventory data, coordinates are (x, y, z), hanging wall related fault set is {F4}, foot wall related fault set is {F3}.
[0099] Initialize the interpolation point data: before interpolation, each channel of the grid is initialized to determine the number of horizon points that need to be interpolated for each channel of the grid (e.g. grid channel 1 has 1 interpolation point and grid channel 3 has 3 interpolation points Point, grid lane 8 has no interpolation points due to the existence of vertical faults), and the fault-related attributes of each interpolation point. Define the interpolation point structure as follows:
[0100] I i ={zValue i , UF i , LF i } Formula (2-5)
[0101] Such as Picture 11 Shown:
[0102] Initialize interpolation point I 1 Attributes of: z value to be calculated, hanging wall associated fault set is {null}, bottom wall associated fault set is {null};
[0103] Initialize interpolation point I 4 Attributes of: z value to be calculated, hanging wall associated fault set is {F1}, bottom wall associated fault set is {F3};
[0104] Initialize interpolation point I 7 The attribute of: z value to be calculated, hanging wall associated fault set is {F1, F2, F3}, bottom wall associated fault set is {null};
[0105] Initialize interpolation point I 8 Attributes of: z value to be calculated, hanging wall associated fault set is {F1}, bottom wall associated fault set is {F3}.
[0106] Horizon interpolation: According to the interpolation point data obtained in the previous step, the interpolation point is calculated 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. After all the interpolation points that can rely on the seed point for interpolation calculation are processed, the interpolation point is regarded as the new seed point for the horizon compensation interpolation, until All interpolation points are processed. In this way, we can generate a complete geological layer surface through layer fitting.
[0107] In the level interpolation algorithm, searching for legal seed points within the allowed search range is one of the most important parts. Let the interpolation point be I i , The associated fault set uF i And lF i; The searched seed point is S j , The associated fault set uF j And lF j. Define the horizontal line between the interpolation point and the seed point I i S j Intersecting fault polygon boundary set Boundary (where uP is the boundary set of the fault hanging wall, lP is the boundary set of the fault bottom wall, up k Is fault F k The boundary of the polygonal hanging wall of the associated interpolation level, lp k Is fault F k The boundary of the polygon bottom plate of the associated interpolation level, m is the total number of fault polygons associated with the interpolation level):
[0108] Boundary = {uP, lP}
[0109] uP = { up i 1 , up i 2 , L , up i M } , 0 ≤ i M ≤ m Formula (2-6)
[0110] lP = { lp k 1 , lp k 2 , L , lp k N } , 0 ≤ k N ≤ m
[0111] In order to get the correct seed point for interpolation calculation, we define the following basic judgment rules:
[0112] A. If And Seed point S j Effective (or legal).
[0113] B. If there is fault F k ∈uF i And F k ∈lF j (Or F k ∈lF i And F k ∈uF j ), then the seed point S j Invalid (or illegal).
[0114] C. If there is a vertical fault boundary in Boundary, the seed point S j invalid.
[0115] D. If the hanging wall boundary and the lower wall boundary of a certain fault exist in Boundary at the same time, the seed point S j invalid.
[0116] E. If uF j The last fault or lF j The first fault in Boundary is a normal fault, and the boundary of this normal fault is included in Boundary, then the seed point S j invalid.
[0117] F if the Boundary set does not satisfy the conditions in D and E, and there is uF i = UF j And lF i =lF j , Then the seed point S j effective.
[0118] G. If there is fault F k ∈uF i With lp k ∈ F k , Lp k ∈lP (or F k ∈lF i , And there is up k ∈ F k , Up k ∈uP), then the seed point S j invalid.
[0119] H. If there is fault F k ∈uF j With lp k ∈ F k , Lp k ∈lP (or F k ∈lF j , And there is up k ∈ F k , Up k ∈uP), then the seed point S j invalid.
[0120] I. In other cases, it is judged as seed point S j effective.
[0121] For the layer interpolation under the boundary constraint of any complex polygon, the judging of the legitimacy of its seed points can be based on the above basic judgment criteria. Such as Picture 12 As shown, according to the above judgment criteria, there are:
[0122] (1) If the interpolation I 1 , S 1 Meet criterion A, valid; for S 5 , There is S 5 ={uF={F 1 }, lF={F 2 }}, uP={null}, lP={lp 1 , Lp 2 }, lp 1 ∈ F 1 , Lp 1 ∈lP satisfies the criterion H and is invalid; for S 2 , In line with criterion I, valid; for S 10 And S 11 , According to criterion D, it is an invalid seed point.
[0123] (2) If the interpolation I 2 , According to criteria F, S 5 Effective; according to criteria B, S 2 Invalid; for seed point S 1 With I 2 ={uF={F 1 }, lF={F 2 }}, uP={null}, lP={lp 1 , Lp 2 }, lp 1 ∈ F 1 , Lp 1 ∈lP satisfies criterion G and is invalid.
[0124] (3) If the interpolation I 4 , According to criterion E, line segment I 4 S 16 Crossing normal fault F 3 So the seed point S 16 invalid.
[0125] (4) If the interpolation I 5 , According to criteria C, S 19 invalid.
[0126] For a certain interpolation point I x For example, if the set of legal seed points found is S′={S 1 , S 2 , L, S k }, and k≥2, then kriging can be used for interpolation. The coordinate of the point set corresponding to S′ is SP={p 1 , P 2 , L, p k }, where p i (1≤i≤k) is the three-dimensional coordinate point (x i , Y i ,z i ), I x The coordinates are (x 0 , Y 0 , Z), z is unknown. Kriging horizon interpolation is to use the known point set SP and the xy coordinates of the interpolation point (x 0 , Y 0 ) To calculate the z value of the interpolation point.
[0127] For some special applications, such as block formation, the horizon boundary is required to be relatively smooth, and too sparse horizon grid will cause obvious boundary jagged. On the other hand, in the 3D seismic interpretation system, if the horizon grid is too fine, the amount of data processing will increase sharply. In order to solve this contradiction, a boundary optimization scheme based on grid nesting is proposed here, that is, in the original grid, if a certain grid is exactly on the boundary of the fault polygon, the grid will be re-meshed and refined. Sub-interpolation, such as Figure 13 Shown.
[0128] The fourth step, the horizon fits the surface
[0129] When the horizon interpolation is completed, the data obtained is still a bunch of discrete three-dimensional points. To obtain the horizon surface, it is also necessary to fit the interpolated data, that is, connect the interpolated grid point data according to the determined rules Formed into a triangle surface, several triangle surfaces are seamlessly connected to become a geological layer surface. In this scheme, because the boundary constraint conditions are determined in advance, in the process of connecting triangles, as long as the boundary interpolation points adjacent to the horizon and the fault are connected with the upper and lower disk line segments of the fault polygon, a complete layer and fault can be obtained. Seamlessly connected surface structure, such as Figure 14 As shown, the problem that the horizon and the fault cannot be completely intersected in the traditional gridding horizon fitting method is solved.