Stratum interface modeling method based on improved Kriging interpolation

Abstract

The invention discloses a stratigraphic interface modeling method based on improved Kriging interpolation, and relates to the field of information models. Specifically, the method is used for fusing occurrence information of a stratum interface in a Kriging interpolation process. The method comprises the following steps: acquiring formation demarcation points and formation interface occurrence information; a penalty factor constructed by formation interface occurrence information is introduced into a variation function fitting process, so that determination of a variation function in a Kriging interpolation process is improved; the grids are encrypted, and then encryption points are determined through improved Kriging interpolation; and the encryption points and the stratum demarcation points are connected according to a space grid connection relation to generate a stratum interface. According to the method, two different types of geological information are fused into the geological model through the penalty factor, and the precision of the geological model is effectively improved.

Description

Technical Field

[0001] The invention relates to the field of information models, in particular to a stratum interface modeling method based on improved Kriging interpolation. Background Art

[0002] The spatial interpolation method is used to generate a spatially continuous stratigraphic interface from the limited stratigraphic boundary point data obtained from the borehole to estimate the underground mineral reserves and distribution. It is a commonly used technical means in the exploration field. Among them, the Kriging interpolation method is widely used due to its high estimation accuracy.

[0003] In the existing geological modeling methods based on Kriging interpolation, only the stratigraphic boundary point information obtained by drilling is considered, and the stratigraphic interface occurrence information obtained by means of mapping is not considered. For example, the invention patent CN109191573A published in 2019 generates the stratigraphic interface skeleton by parabolic interpolation, and then generates the stratigraphic interface by Kriging interpolation, which improves the efficiency and quality of modeling; the invention patent CN115422830A published in 2022 replaces the commonly used theoretical variation function with the SVM support vector machine algorithm to improve the accuracy of spatial interpolation.

[0004] Existing scholars have improved the Kriging interpolation by using existing technologies to improve the modeling accuracy, but have not considered combining the stratigraphic boundary point data with the stratigraphic interface occurrence data to improve the modeling accuracy during the modeling process. Summary of the invention

[0005] In order to overcome the defects in the above-mentioned technology, the present invention provides a formation interface modeling method based on improved Kriging interpolation, which realizes the fusion modeling of formation boundary point data and formation interface occurrence data by introducing a penalty factor constructed by the formation interface occurrence information in the process of variogram fitting.

[0006] To achieve the above purpose, the present invention adopts the following technical solutions: A stratum interface modeling method based on improved Kriging interpolation includes the following steps: Step 1, obtaining the spatial coordinate information of the stratigraphic boundary point; obtaining the occurrence information of the stratigraphic boundary interface and the location of the occurrence information; Step 2, calculate the semi-vertical distance difference and horizontal distance between the stratigraphic boundary points, and construct the penalty factor based on the occurrence information to obtain the final theoretical variation function by fitting; Step 3: Set the spatial interpolation encryption size, calculate the plane coordinates of the encrypted projection points according to the encrypted size, and then determine the plane coordinates of the encrypted projection points by the spatial interpolation algorithm. Z Axis coordinates, thereby determining the encryption points and obtaining the spatial grid connection relationship obtained during the encryption process; Step 4: Connect the encrypted points and the stratigraphic boundary points according to the spatial grid connection relationship to generate a stratigraphic interface.

[0007] Preferably, in step 1, The spatial coordinates of the stratigraphic boundary points are recorded in the form of a spatial rectangular coordinate system, and the recorded values ​​are X , Y , Z The values ​​corresponding to the three coordinate axis directions, where X , Y The two coordinate axes lie in the horizontal plane, and X The positive direction of the axis is right. Z The positive direction of the axis is upward. Y The positive direction of the axis is determined by the left-hand rule; assuming that the stratum boundary points have n , i The spatial coordinates of the stratigraphic boundary points are ( x i , y i , z i ), i =1, 2, ..., n ; The occurrence information of the stratum interface includes the actual inclination and the actual dip angle of the actual occurrence section at the location of the occurrence information; the location of the occurrence information refers to the plane coordinates of the location of the occurrence information; assuming that the occurrence information has m , k The plane coordinates of the location of the occurrence information are ( x k , y k ), the actual inclination and actual dip angle of the corresponding actual occurrence section are recorded as φ k and θ k , k =1, 2, ..., m .

[0008] Preferably, the specific process of step 2 is as follows: Step 2.1: Calculate any two stratum boundary points according to the semi-vertical distance difference calculation formula i , j Half vertical difference r ij , according to the horizontal distance calculation formula, we can calculate any two stratum boundary points i , j Horizontal distance d ij ; i , j =1, 2, ...,n , the stratigraphic boundary points are n indivual; Step 2.2, select a theoretical variation function r (·), the theoretical variation function r (·) contains two parameters that need to be fitted C 0. C 1. C 0 is recorded as the first parameter to be fitted, C 1 is recorded as the second parameter to be fitted; Introducing an objective function Obj : ; In the formula, r ( d ij ) is the horizontal distance d ij Substitute the theoretical variogram r The semivariance calculated in (·) is the stratigraphic boundary point i , j The semivariance of For the k The calculated tendency of the calculated occurrence section at the location of each occurrence information; For the k The calculated inclination angle of the calculated occurrence section at the location of the occurrence information; φ k For the k The actual inclination of the actual occurrence section where the occurrence information is located; θ k For the k The actual inclination angle of the actual occurrence section at the location of the occurrence information; k =1, 2, ..., m , the production information is m indivual; H is the penalty factor, which is and Positively correlated functions; Step 2.3, set the first parameter to be fitted C 0 and the second parameter to be fitted C The initial value of 1 is used to calculate the objective function Obj The value of Step 2.4: Adjust the first parameter to be fitted according to the gradient descent method C 0 and the second parameter to be fitted C 1, and solve to get the objective function Obj The minimum value of is recorded as the temporary minimum value; Step 2.5, randomly set the first fitting parameters in the new round of calculationsC 0 and the second parameter to be fitted C 1, repeat steps 2.3-2.4 to obtain the objective function in the new round of calculation Obj The minimum value of is calculated and compared with the temporary minimum value; if the objective function in the new round of calculation Obj If the minimum value of is less than the temporary minimum value, the temporary minimum value is replaced with the minimum value of the objective function in the new round of calculation; otherwise, it is not replaced; Step 2.6, repeat step 2.5, in succession w The temporary minimum value is replaced when no replacement occurs; w is a set parameter; the temporary minimum value obtained at the end is the objective function Obj The final minimum value of the objective function Obj The theoretical variation function corresponding to the final minimum value of is the final theoretical variation function.

[0009] Preferably, the specific process of step 3 is as follows: Step 3.1, set the encryption size to a ; Get n The spatial coordinates of the stratigraphic boundary points are projected onto the horizontal plane to obtain n projection point, the n The projection point is recorded as the basic projection point; Step 3.2, connect the basic projection points according to the triangle, ensure that each basic projection point is a vertex of the triangle, and the triangle vertices are all basic projection points, that is, form a basic triangle mesh; Step 3.3: According to the principle of equally divided triangles, encrypt the basic triangle mesh until the average side length of each triangle is less than a ; All points other than the basic projection points in the mesh vertices of the encrypted basic triangular mesh are recorded as the encrypted projection points; the connection relationship between the basic projection points and the encrypted projection points is recorded according to the encrypted basic triangular mesh, and recorded as the plane mesh connection relationship; Step 3.4, select an encrypted projection point g , get its horizontal coordinates ( x g , y g ), calculate the encrypted projection point g With i The horizontal distance between the stratigraphic boundary points is recorded as the encrypted projection horizontal distance d ig , i =1, 2, ..., n ; Then use the encrypted projection horizontal distance d ig Replace the horizontal distance between any two stratigraphic boundary pointsd ij Substitute the final theoretical variance function obtained in step 2 into r (·), get the encrypted projection point g With i The semivariance of the stratigraphic boundary points is recorded as the projected semivariance r ( d ig ); Step 3.5: Calculate the encrypted projection points in the same way as step 3.4 g Respectively n The horizontal distance between the stratigraphic boundary points is n Encrypted projection horizontal distance d ig , n Projected semivariance r ( d ig ); according to the weight coefficient matrix given n Projected semivariance r ( d ig ) The corresponding weight group λ 1g , λ 2g , ..., λ ig , ..., λ ng , and substitute it into the weight formula to get the encrypted projection point g of Z Axis coordinate estimation z g ; Encrypt the projected points Z Axis coordinate estimation z g With horizontal coordinates ( x g , y g ) to form an encryption point ( x g , y g , z g ); Repeat steps 3.4-3.5 until all encrypted projection points generate a corresponding encrypted point; Step 3.6, in the plane grid connection relationship, the corresponding encrypted projection points are replaced by the encrypted points, and the corresponding basic projection points are replaced by the stratum boundary points to form a spatial grid connection relationship.

[0010] Preferably, in step 2.3, the objective function is calculated Obj The specific process is as follows: Step 2.3.1, set the theoretical variation function r (·) The initial values ​​of the parameters; Step 2.3.2, according to the horizontal distance d ij , half vertical distance difference r ij Calculate the objective function Obj The first half of Obj sub1 The value of ; Step 2.3.3, select a production information and obtain the k The location of the production information and its b The plane coordinates of the adjacent positions are recorded as the surrounding positions of the occurrence, and a total of b +1 occurrence peripheral position; calculate the horizontal distances between the occurrence peripheral position and all stratum boundary points in sequence, store the horizontal distances between each occurrence peripheral position and all stratum boundary points in an occurrence horizontal distance group, and obtain b +1 occurrence horizontal distance group, each occurrence horizontal distance group has n Horizontal distance d ip , p is the number of the position around the occurrence, p =1, 2, ..., b +1; each horizontal distance in each occurrence horizontal distance group is sequentially brought into the theoretical variation function r Calculate the semivariance in (·) and substitute each occurrence horizontal distance group into the theoretical variogram r (·) The calculated results are stored separately, and we get b +1 occurrence semivariance group, each with n Semivariance r ( d ip ); Step 2.3.4, substitute the semivariance group corresponding to each occurrence peripheral position into the weight coefficient matrix, and solve to obtain the weight group of each occurrence peripheral position λ 1p , λ 2p , ..., λ ip , ..., λ np , and then the weight group of each position around the occurrence is respectively brought into the weight formula to calculate the weight of each position around the occurrence Z Axis coordinate estimation z p ; Step 2.3.5: Compare the horizontal coordinates of each occurrence peripheral position with Z Axis coordinate estimation z p Combining them, we get a set of space points, recorded as calculation space points, and use this set of calculation space points to fit a calculation occurrence section as the first k The calculated occurrence section where each occurrence information is located; According to the determination method of the actual inclination and actual dip of the actual occurrence section, the first k The calculated tendency of the calculated occurrence section at the location of the occurrence information and calculate the inclination ; Step 2.3.6, repeat steps 2.3.3-2.3.5 to obtain the calculated inclination and calculated dip angle corresponding to all occurrence information, and then give the penalty factor H , calculate the objective function Obj The second half of Obj sub2 The value of .

[0011] Preferably, the actual inclination and the actual dip angle of the actual strike section are determined as follows: a strike line and a dip line are defined, wherein the strike line is the intersection line of the actual strike section and the horizontal plane, and the dip line is located in the actual strike section, perpendicular to the strike line and along the horizontal plane. Z The component direction of the axis direction is negative; the actual inclination refers to the direction pointed by the projection of the inclination line on the horizontal plane; the actual inclination angle refers to the angle between the inclination line and the horizontal plane; The actual strike section at the location of the strike information refers to the spatial section of the stratigraphic interface at the location of the strike information.

[0012] Preferably, the semi-vertical distance difference calculation formula and the horizontal distance calculation formula are respectively: ; .

[0013] Preferably, the weight coefficient matrix and weight formula are as follows: ; ; In the formula, is the Lagrange multiplier, subscript v Indicates the number, which is the number of the encrypted projection point g Or it is the number of the position around the occurrence p ; r ij To calculate the boundary points of any two strata according to the semi-vertical distance difference calculation formula i ,j The half vertical difference; z i For the i Stratigraphic boundary Z Axis coordinates.

[0014] Preferably, the specific process of step 2.4 is as follows: Step 2.4.1, take the first parameter to be fitted C 0 and the second parameter to be fitted C 1 is the two directions of the coordinate axis, and a plane coordinate system is established; Step 2.4.2, set the first parameter to be fitted C 0 and the second parameter to be fitted C The initial value of 1 is the parameter point of the plane coordinate system; Step 2.4.3: Move the parameter point in different directions according to a fixed length in the plane coordinate system to obtain different C 0. C 1, calculate the values ​​of all objective functions, and select the trial parameter point with the smallest objective function value as the new parameter point until the minimum value of the objective function is obtained.

[0015] Preferably, in step 2.3.3, the location of the production information b The adjacent positions satisfy the following conditions: The horizontal distance between each adjacent position and the position of the occurrence information is less than e ,and b The points obtained by projecting the adjacent position and the position of the occurrence information onto the horizontal plane simultaneously cannot be located on the same straight line. b ≥2, e is a set parameter; in step 2.3.5, the calculated attitude section satisfies that the sum of the squares of the vertical distances of all calculated spatial points from the calculated attitude section is minimal.

[0016] The beneficial effects of the present invention are: (1) By introducing a penalty factor constructed based on the formation interface strike information as a constraint in the fitting process of the variogram, the fusion modeling of the formation boundary point data and the formation interface strike data in the Kriging interpolation process is realized, which effectively improves the modeling accuracy.

[0017] (2) In the process of fitting the variogram, the gradient descent method is used to solve the minimum value of the objective function under different initial values, and the minimum value of the objective function is determined by comparing the minimum values, thus ensuring the accuracy of the variogram fitting.

[0018] (3) By introducing the penalty factor, the fusion of two different types of data, stratigraphic boundary points and stratigraphic interface attitude, is realized. The algorithm is simple to implement, which provides a reference idea for introducing other types of data in the process of geological space interpolation. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 The present invention is a flow chart of the formation interface modeling method.

[0020] Figure 2 Schematic diagram of plane mesh encryption.

[0021] Figure 3 Schematic diagram of the actual occurrence section and the calculated occurrence section.

[0022] Figure 4 Schematic diagram of parameter point movement in the process of solving the minimum value of the objective function.

[0023] Figure 5 A simplified diagram of the formation interface modeling method of the present invention.

[0024] The following are the descriptions of the reference numerals: 1-drill hole, 2-stratum boundary point, 3-basic projection point, 4-actual occurrence section, 5-encrypted projection point, 6- Y Axis, 7- X axis, 8-calculated attitude section, 9-actual strike line, 10-calculated strike line, 11-actual dip line, 12-calculated dip line, 13-actual dip, 14-calculated dip, 15-actual dip angle, 16-calculated dip angle, 17-parameter point, 18- Z axis. DETAILED DESCRIPTION

[0025] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0026] Figure 1 is a flow chart of the formation interface modeling method of the present invention, Figure 5 FIG. 1 is a simplified diagram of the formation interface modeling method of the present invention. Figure 1 and Figure 5 It can be seen that the present invention provides a stratum interface modeling method based on improved Kriging interpolation, wherein the stratum interface refers to the interface between two different strata, comprising the following steps: Step 1, obtain the spatial coordinate information of the stratigraphic boundary point; obtain the occurrence information of the stratigraphic boundary interface and the location of the occurrence information.

[0027] The spatial coordinates of the stratigraphic boundary points are recorded in the form of a spatial rectangular coordinate system, and the recorded values ​​are X , Y, Z The values ​​corresponding to the three coordinate axis directions, where X , Y The two coordinate axes lie in the horizontal plane, and X The positive direction of the axis is right. Z The positive direction of the axis is upward. Y The positive direction of the axis is determined by the left-hand rule; assuming that the stratum boundary point is , and record any one of them as i A stratigraphic dividing point, i is the number of the stratigraphic boundary point. i =1, 2, ..., n ;No. i The spatial coordinates of the stratigraphic boundary points are ( x i , y i , z i ).

[0028] The occurrence information of the stratum interface includes the actual inclination and actual dip angle of the actual occurrence section at the location of the occurrence information; the location of the occurrence information refers to the plane coordinates of the location of the occurrence information; assuming that the occurrence information has m , and record any one of them as k Birth information, k The plane coordinates of the location of the occurrence information are ( x k , y k ), the actual inclination and actual dip angle of the corresponding actual occurrence section are recorded as φ k and θ k , k =1, 2, ..., m .

[0029] In this embodiment, the actual inclination and actual dip of the actual strike section are determined as follows: the concepts of strike line and dip line are introduced, the strike line is the intersection line of the actual strike section and the horizontal plane, and the dip line is located in the actual strike section, perpendicular to the strike line and along the horizontal plane. Z The component in the axial direction is negative; the actual inclination refers to the direction of the projection of the inclination line on the horizontal plane, and the actual inclination angle refers to the angle between the inclination line and the horizontal plane.

[0030] In this embodiment, the actual strike section at the position where the strike information is located refers to the spatial section of the stratum interface at the position where the strike information is located.

[0031] like Figure 2As shown, there are eight stratigraphic boundary points of a stratigraphic interface obtained by drilling, and one stratigraphic interface occurrence and its horizontal coordinate.

[0032] like Figure 3 As shown in the figure, there is a formation interface occurrence information in the spatial coordinate system. The actual occurrence section can be uniquely determined by the actual dip line and the actual strike line. The actual dip is recorded as X The axis rotates clockwise to the angle between the projection line of the actual dip line and the actual dip angle. The actual dip angle is the angle between the actual dip line and the horizontal plane. The actual dip line is perpendicular to the actual strike line and is located within the actual strike section. Z The axial direction component is negative.

[0033] exist Figure 2 and Figure 3 In the figure, the reference numerals are explained as follows: 1 is a borehole, 2 is a stratigraphic boundary point, 3 is a basic projection point, 4 is an actual occurrence section, 5 is an encrypted projection point, 6 is a Y-axis, 7 is an X-axis, 8 is a calculated occurrence section, 9 is an actual strike line, 10 is a calculated strike line, 11 is an actual inclination line, 12 is a calculated inclination line, 13 is an actual inclination, 14 is a calculated inclination, 15 is an actual dip angle, 16 is a calculated dip angle, and 18 is a Z-axis.

[0034] Step 2: Calculate the semi-vertical distance difference and horizontal distance between the stratigraphic boundary points, construct the penalty factor based on the occurrence information, and fit the final theoretical variation function.

[0035] Step 2.1: Calculate any two stratum boundary points according to the semi-vertical distance difference calculation formula i , j Half vertical difference r ij , according to the horizontal distance calculation formula, we can calculate any two stratum boundary points i , j Horizontal distance d ij , a total of The half vertical difference and Horizontal distance; i , j =1, 2, ..., n .

[0036] In this embodiment, the semi-vertical distance difference calculation formula and the horizontal distance calculation formula are respectively: ; .

[0037] like Figure 2 As shown, in this embodiment, eight stratigraphic boundary points are obtained, so thirty-two semi-vertical distance differences and horizontal distances can be obtained.

[0038] Step 2.2, select a theoretical variation function r ( d ), the theoretical variation function r ( d ) contains two parameters that need to be fitted C 0. C 1. A manually set parameter δ and an independent variable parameter d ;Will C 0 is recorded as the first parameter to be fitted, C 1 is recorded as the second parameter to be fitted.

[0039] Introducing an objective function Obj : ; In the formula, r ( d ij ) is the horizontal distance d ij Substitute the theoretical variogram r ( d ), i.e., the semivariance calculated from the stratigraphic boundary point i , j The semivariance of For the k The calculated tendency of the calculated occurrence section at the location of each occurrence information; For the k The calculated inclination angle of the calculated occurrence section at the location of the occurrence information; φ k For the k The actual inclination of the actual occurrence section where the occurrence information is located; θ k For the k The actual inclination angle of the actual occurrence section at the location of the occurrence information; k =1, 2, ..., m ; H is the penalty factor, which is and Positively correlated function.

[0040] In this embodiment, the exponential function is selected As the theoretical variogram. The independent variable parameter here is d It is the horizontal distance between the boundary points of two strata.

[0041] Step 2.3, set the first parameter to be fitted C 0 and the second parameter to be fitted C The initial value of 1 is used to calculate the objective function ObjThe value of .

[0042] Step 2.4: Adjust the first parameter to be fitted according to the gradient descent method C 0 and the second parameter to be fitted C 1, and solve to get the objective function Obj and record it as the temporary minimum.

[0043] Step 2.5: Randomly set the first fitting parameter in the new round of calculation C 0 and the second parameter to be fitted C 1, repeat steps 2.3-2.4 to obtain the objective function in the new round of calculation Obj The minimum value of is calculated and compared with the temporary minimum value; if the objective function in the new round of calculation Obj If the minimum value of is less than the temporary minimum value, the temporary minimum value is replaced with the minimum value of the objective function in the new round of calculation; otherwise, it is not replaced.

[0044] Step 2.6, repeat step 2.5, in succession w The temporary minimum value is replaced when no replacement occurs; w is a manually set parameter; the temporary minimum value obtained at the end is the objective function Obj The final minimum value of the objective function Obj The theoretical variation function corresponding to the final minimum value of is the final theoretical variation function.

[0045] Step 3: Set the spatial interpolation encryption size, calculate the plane coordinates of the encrypted projection points according to the encrypted size, and then determine the plane coordinates of the encrypted projection points by the spatial interpolation algorithm. Z Axis coordinates are used to determine the encryption points and obtain the spatial grid connection relationship obtained during the encryption process.

[0046] Step 3.1, set the encryption size to a ; Get n The spatial coordinates of the stratigraphic boundary points are projected onto the horizontal plane to obtain n projection point, the n The projection point is recorded as the basic projection point.

[0047] Step 3.2, connect the basic projection points according to the triangle, ensure that each basic projection point is a vertex of the triangle, and the triangle vertices are all basic projection points, that is, form a basic triangle mesh.

[0048] like Figure 2 As shown, each stratigraphic boundary point is projected into the horizontal plane to obtain eight basic projection points, and then the eight basic projection points are connected through a triangular mesh to form a triangular mesh.

[0049] Step 3.3: According to the principle of equally divided triangles, encrypt the basic triangle mesh until the average side length of each triangle is less than a ; All points other than the basic projection points in the mesh vertices of the encrypted basic triangular mesh are recorded as encrypted projection points; the connection relationship between the basic projection points and the encrypted projection points is recorded according to the encrypted basic triangular mesh, and recorded as the plane mesh connection relationship.

[0050] In this embodiment, the principle of equally dividing triangles refers to connecting the midpoint of the longest grid line segment of the triangular mesh with all the triangle vertices corresponding to the longest grid line segment to divide the existing triangle; the triangle vertex corresponding to the longest grid line segment refers to the vertex corresponding to the grid line segment when it serves as the side of the triangle.

[0051] Step 3.4, select an encrypted projection point g , get its horizontal coordinates ( x g , y g ), calculate the encrypted projection point g With i The horizontal distance between the stratigraphic boundary points is recorded as the encrypted projection horizontal distance d ig , ; Then use the encrypted projection horizontal distance d ig Replace the horizontal distance between any two stratigraphic boundary points d ij Substitute the final theoretical variation function obtained in step 2 to obtain the encrypted projection point g With i The semivariance of the stratigraphic boundary points is recorded as the projected semivariance r ( d ig ).

[0052] Step 3.5: Calculate the encrypted projection points in the same way as step 3.4 g Respectively n The horizontal distance between the stratigraphic boundary points is n Encrypted projection horizontal distance d ig , n Projected semivariance r ( d ig ); according to the weight coefficient matrix given n Projected semivariance r ( d ig ) The corresponding weight group λ 1g, λ 2g , ..., λ ig , ..., λ ng , and substitute it into the weight formula to get the encrypted projection point g of Z Axis coordinate estimation z g ; Encrypt the projection points Z Axis coordinate estimation z g With horizontal coordinates ( x g , y g ) to form an encryption point ( x g , y g , z g ).

[0053] Repeat steps 3.4-3.5 until all encrypted projection points generate a corresponding encrypted point.

[0054] Step 3.6, in the plane grid connection relationship, the corresponding encrypted projection points are replaced by the encrypted points, and the corresponding basic projection points are replaced by the stratum boundary points to form a spatial grid connection relationship.

[0055] Step 4, obtaining the encrypted points, the spatial grid connection relationship and the stratum boundary points, and connecting the encrypted points and the stratum boundary points according to the spatial grid connection relationship to generate a stratum interface.

[0056] At this point, the modeling of the stratigraphic interface is completed.

[0057] In this embodiment, the calculation objective function in step 2.3 is Obj The specific process is as follows: Step 2.3.1, set the theoretical variation function r ( d ) in the initial value of the parameter.

[0058] Step 2.3.2, according to the horizontal distance d ij , half vertical distance difference r ij Calculate the objective function Obj The first half of Obj sub1 The value of .

[0059] Step 2.3.3, select a production information and obtain thek The location of the production information and its b The plane coordinates of the adjacent positions are recorded as the surrounding positions of the occurrence, and a total of b +1 occurrence peripheral position; calculate the horizontal distances between the occurrence peripheral position and all stratum boundary points in sequence, store the horizontal distances between each occurrence peripheral position and all stratum boundary points in an occurrence horizontal distance group, and obtain b +1 occurrence horizontal distance group, each occurrence horizontal distance group has n Horizontal distance d ip , p is the number of the position around the occurrence, p =1, 2, ..., b +1, each attitude horizontal distance group corresponds to the attitude surrounding position participating in the calculation of this group; each horizontal distance in each attitude horizontal distance group is sequentially brought into the theoretical variation function r ( d ) and substitute each occurrence horizontal distance group into the theoretical variogram r ( d ) The results of the calculations are stored separately, and the b +1 occurrence semivariance group, each with n Semivariance r ( d ip ), each occurrence semivariance group corresponds to the occurrence peripheral position corresponding to the occurrence horizontal distance group participating in the calculation of this group.

[0060] In this embodiment, the location of the production status information b The adjacent positions satisfy the following conditions: The horizontal distance between each adjacent position and the position of the occurrence information is less than e ,and b The points obtained by simultaneously projecting the adjacent position and the position of the occurrence information onto the horizontal plane cannot be located on the same straight line; b ≥2; e It is a manually set parameter.

[0061] Step 2.3.4, substitute the semivariance group corresponding to each occurrence peripheral position into the weight coefficient matrix, and solve to obtain the weight group of each occurrence peripheral position λ 1p , λ 2p , ..., λ ip , ..., λ np , and then the weight group of each position around the occurrence is respectively brought into the weight formula to calculate the weight of each position around the occurrence ZAxis coordinate estimation z p .

[0062] Step 2.3.5: Compare the horizontal coordinates of each occurrence peripheral position with Z Axis coordinate estimation z p Combining them, we get a set of space points, recorded as calculation space points, and use this set of calculation space points to fit a calculation occurrence section as the first k The calculated occurrence section at the location of the occurrence information.

[0063] According to the method for determining the actual inclination and actual dip of the actual occurrence section described in step 1, the calculated inclination of the calculated occurrence section is obtained. and calculate the inclination ; The calculated attitude section satisfies that the sum of the squares of the vertical distances of all calculated spatial points from the calculated attitude section is minimum.

[0064] Step 2.3.6, repeat steps 2.3.3-2.3.5 to obtain the calculated inclination and calculated dip angle corresponding to all occurrence information, and then give the penalty factor H , calculate the objective function Obj The second half of Obj sub2 The value of .

[0065] In this embodiment, the first parameter to be fitted is adjusted according to the gradient descent method in step 2.4. C 0 and the second parameter to be fitted C 1, and solve for the objective function Obj The specific process is as follows: Step 2.4.1, take the first parameter to be fitted C 0 and the second parameter to be fitted C 1 is the two directions of the coordinate axis, and a plane coordinate system is established.

[0066] Step 2.4.2, set the first parameter to be fitted C 0 and the second parameter to be fitted C The initial value of 1 is the parameter point of the plane coordinate system.

[0067] Step 2.4.3: Move the parameter point in different directions according to a fixed length in the plane coordinate system to obtain different C 0. C 1, execute steps 2.3.2-2.3.6, calculate the values ​​of all objective functions, and select the trial parameter point with the smallest objective function value as the new parameter point until the minimum value of the objective function is obtained.

[0068] Figure 4 Solving the process parameter points for the minimum value of the objective function C 0. C 1 is a moving diagram. Figure 4 Reference numeral 17 is a parameter point.

[0069] In this embodiment, the weight coefficient matrix and weight formula are as follows: ; ; In the formula, is the Lagrange multiplier, subscript v Indicates the number, which is the number of the encrypted projection point g Or it is the number of the position around the occurrence p .

[0070] The above are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A stratum interface modeling method based on improved Kriging interpolation, characterized in that: The following steps are involved: Step 1, obtaining the spatial coordinate information of the stratigraphic boundary point; obtaining the occurrence information of the stratigraphic boundary interface and the location of the occurrence information; Step 2, calculate the semi-vertical distance difference and horizontal distance between the stratigraphic boundary points, and construct the penalty factor based on the occurrence information to obtain the final theoretical variation function by fitting; Step 3: Set the spatial interpolation encryption size, calculate the plane coordinates of the encrypted projection points according to the encrypted size, and then determine the plane coordinates of the encrypted projection points by the spatial interpolation algorithm. Z Axis coordinates, thereby determining the encryption points and obtaining the spatial grid connection relationship obtained during the encryption process; Step 4: Connect the encrypted points and the stratigraphic boundary points according to the spatial grid connection relationship to generate a stratigraphic interface.

2. A stratum interface modeling method based on improved Kriging interpolation according to claim 1, characterized in that: In step 1, The spatial coordinates of the stratigraphic boundary points are recorded in the form of a spatial rectangular coordinate system, and the recorded values ​​are X , Y , Z The values ​​corresponding to the three coordinate axis directions, where X , Y The two coordinate axes lie in the horizontal plane, and X The positive direction of the axis is right. Z The positive direction of the axis is upward. Y The positive direction of the axis is determined by the left-hand rule; assuming that the stratigraphic boundary points have n , i The spatial coordinates of the stratigraphic boundary points are ( x i , y i , z i ), i =1, 2, ..., n ; The occurrence information of the stratum interface includes the actual inclination and actual dip angle of the actual occurrence section at the location of the occurrence information; the location of the occurrence information refers to the plane coordinates of the location of the occurrence information; assuming that the occurrence information has m , k The plane coordinates of the location of the occurrence information are ( x k , y k ), the actual inclination and actual dip angle of the corresponding actual occurrence section are recorded as φ k and θ k , k =1, 2, ..., m .

3. The method for modeling stratum interface based on improved Kriging interpolation according to claim 1 is characterized in that: The specific process of step 2 is as follows: Step 2.1: Calculate any two stratum boundary points according to the semi-vertical distance difference calculation formula i , j Half vertical distance difference r ij , according to the horizontal distance calculation formula, we can calculate any two stratum boundary points i , j Horizontal distance d ij ; i , j =1, 2, ..., n , the stratigraphic boundary points are n indivual; Step 2.2, select a theoretical variation function r (·), the theoretical variation function r (·) contains two parameters that need to be fitted C 0. C 1. C 0 is recorded as the first parameter to be fitted, C 1 is recorded as the second parameter to be fitted; Introducing an objective function Obj : ； In the formula, r ( d ij ) is the horizontal distance d ij Substitute the theoretical variogram r The semivariance calculated in (·) is the stratigraphic boundary point i , j The semivariance of For the k The calculated tendency of the calculated occurrence section at the location of each occurrence information; For the k The calculated inclination angle of the calculated occurrence section at the location of the occurrence information; φ k For the k The actual inclination of the actual occurrence section where the occurrence information is located; θ k For the k The actual inclination angle of the actual occurrence section at the location of the occurrence information; k =1, 2, ..., m , the production information is m indivual; H is the penalty factor, which is and Positively correlated functions; Step 2.3, set the first parameter to be fitted C 0 and the second parameter to be fitted C The initial value of 1 is used to calculate the objective function Obj The value of Step 2.4: Adjust the first parameter to be fitted according to the gradient descent method C 0 and the second parameter to be fitted C 1, and solve to get the objective function Obj The minimum value of is recorded as the temporary minimum value; Step 2.5, randomly set the first fitting parameters in the new round of calculations C 0 and the second parameter to be fitted C 1, repeat steps 2.3-2.4 to obtain the objective function in the new round of calculation Obj The minimum value of is calculated and compared with the temporary minimum value; if the objective function in the new round of calculation Obj If the minimum value of is less than the temporary minimum value, the temporary minimum value is replaced as the minimum value of the objective function in the new round of calculation; Otherwise, no replacement; Step 2.6, repeat step 2.5, in succession w The temporary minimum value is replaced when no replacement occurs; w is a set parameter; The temporary minimum obtained at the end is the objective function Obj The final minimum value of the objective function Obj The theoretical variation function corresponding to the final minimum value of is the final theoretical variation function.

4. The method for modeling stratum interface based on improved Kriging interpolation according to claim 1, characterized in that: The specific process of step 3 is as follows: Step 3.1, set the encryption size to a ; Get n The spatial coordinates of the stratigraphic boundary points are projected onto the horizontal plane to obtain n projection point, the n The projection point is recorded as the basic projection point; Step 3.2, connect the basic projection points according to the triangle, ensure that each basic projection point is a vertex of the triangle, and the triangle vertices are all basic projection points, that is, form a basic triangle mesh; Step 3.3: According to the principle of equally divided triangles, encrypt the basic triangle mesh until the average side length of each triangle is less than a ; All points other than the basic projection points in the mesh vertices of the encrypted basic triangular mesh are recorded as encrypted projection points; the connection relationship between the basic projection points and the encrypted projection points is recorded according to the encrypted basic triangular mesh, and recorded as the plane mesh connection relationship; Step 3.4, select an encrypted projection point g , get its horizontal coordinates ( x g , y g ), calculate the encrypted projection point g With i The horizontal distance between the stratigraphic boundary points is recorded as the encrypted projection horizontal distance d ig , i =1, 2, ..., n ; Then use the encrypted projection horizontal distance d ig Replace the horizontal distance between any two stratigraphic boundary points d ij Substitute the final theoretical variance function obtained in step 2 into r (·), get the encrypted projection point g With i The semivariance of the stratigraphic boundary points is recorded as the projected semivariance r ( d ig ); Step 3.5: Calculate the encrypted projection points in the same way as step 3.4 g Respectively n The horizontal distance between the stratigraphic boundary points is n Encrypted projection horizontal distance d ig , n Projected semivariance r ( d ig ); according to the weight coefficient matrix given n Projected semivariance r ( d ig ) The corresponding weight group λ 1g , λ 2g , ..., λ ig , ..., λ ng , and substitute it into the weight formula to get the encrypted projection point g of Z Axis coordinate estimation z g ; Encrypt the projection points Z Axis coordinate estimation z g With the horizontal coordinate ( x g , y g ) to form an encryption point ( x g , y g , z g ); Repeat steps 3.4-3.5 until all encrypted projection points generate a corresponding encrypted point; Step 3.6, in the plane grid connection relationship, the corresponding encrypted projection points are replaced by the encrypted points, and the corresponding basic projection points are replaced by the stratum boundary points to form a spatial grid connection relationship.

5. The method for modeling stratum interface based on improved Kriging interpolation according to claim 3 is characterized in that: In step 2.3, calculate the objective function Obj The specific process is as follows: Step 2.3.1, set the theoretical variation function r (·) The initial values ​​of the parameters; Step 2.3.2, according to the horizontal distance d ij , half vertical distance difference r ij Calculate the objective function Obj The first half of Obj sub1 The value of ； Step 2.3.3, select a production information and obtain the k The location of the production information and its b The plane coordinates of the adjacent positions are recorded as the surrounding positions of the occurrence, and a total of b +1 occurrence peripheral position; Calculate the horizontal distances between the occurrence peripheral position and all the stratum boundary points in sequence, store the horizontal distances between each occurrence peripheral position and all the stratum boundary points in an occurrence horizontal distance group, and obtain b +1 occurrence horizontal distance group, each occurrence horizontal distance group has n Horizontal distance d ip , p is the number of the position around the occurrence, p =1, 2, ..., b +1; each horizontal distance in each occurrence horizontal distance group is sequentially brought into the theoretical variation function r Calculate the semivariance in (·) and substitute each occurrence horizontal distance group into the theoretical variogram r (·) The calculated results are stored separately, and we get b +1 occurrence semivariance group, each with n Semivariance r ( d ip ); Step 2.3.4, substitute the semivariance group corresponding to each occurrence peripheral position into the weight coefficient matrix, and solve to obtain the weight group of each occurrence peripheral position λ 1p , λ 2p , ..., λ ip , ..., λ np , and then the weight group of each position around the occurrence is respectively brought into the weight formula to calculate the weight of each position around the occurrence Z Axis coordinate estimation z p ; Step 2.3.5: Compare the horizontal coordinates of each occurrence peripheral position with Z Axis coordinate estimation z p Combine them to get a set of space points, record them as calculation space points, and use this set of calculation space points to fit a calculation occurrence section as the first k The calculated occurrence section where each occurrence information is located; According to the determination method of the actual inclination and actual dip of the actual occurrence section, the first k The calculated tendency of the calculated occurrence section at the location of the occurrence information and calculate the inclination ; Step 2.3.6, repeat steps 2.3.3-2.3.5 to obtain the calculated inclination and calculated dip angle corresponding to all occurrence information, and then give the penalty factor H , calculate the objective function Obj The second half of Obj sub2 The value of 。 6. A stratum interface modeling method based on improved Kriging interpolation according to any one of claims 2 to 5, characterized in that: The actual inclination and actual dip of the actual strike section are determined as follows: a strike line and a dip line are defined, wherein the strike line is the intersection of the actual strike section and the horizontal plane, and the dip line is located in the actual strike section, perpendicular to the strike line and along the horizontal plane. Z The component direction of the axis direction is negative; the actual inclination refers to the direction pointed by the projection of the inclination line on the horizontal plane; the actual inclination angle refers to the angle between the inclination line and the horizontal plane; The actual strike section at the location of the strike information refers to the spatial section of the stratigraphic interface at the location of the strike information.

7. A stratum interface modeling method based on improved Kriging interpolation according to any one of claims 3 to 5, characterized in that: The semi-vertical distance difference calculation formula and the horizontal distance calculation formula are respectively: ； 。 8. A stratum interface modeling method based on improved Kriging interpolation according to claim 4 or 5, characterized in that: The weight coefficient matrix and weight formula are as follows: ； ； In the formula, is the Lagrange multiplier, subscript v Indicates the number, which is the number of the encrypted projection point g Or it is the number of the position around the occurrence p ; r ij To calculate the boundary points of any two strata according to the semi-vertical distance difference calculation formula i , j The half vertical difference; z i For the i Stratigraphic boundary Z Axis coordinates.

9. The method for modeling stratum interface based on improved Kriging interpolation according to claim 3, characterized in that: The specific process of step 2.4 is as follows: Step 2.4.1, take the first parameter to be fitted C 0 and the second parameter to be fitted C 1 is the two directions of the coordinate axis, and a plane coordinate system is established; Step 2.4.2, set the first parameter to be fitted C 0 and the second parameter to be fitted C The initial value of 1 is the parameter point of the plane coordinate system; Step 2.4.3: Move the parameter point in different directions according to a fixed length in the plane coordinate system to obtain different C 0. C 1, calculate the values ​​of all objective functions, and select the trial parameter point with the smallest objective function value as the new parameter point until the minimum value of the objective function is obtained.

10. The method for modeling stratum interface based on improved Kriging interpolation according to claim 5, characterized in that: In step 2.3.3, the location of the occurrence information b The adjacent positions satisfy the following conditions: The horizontal distance between each adjacent position and the position of the occurrence information is less than e ,and b The points obtained by projecting the adjacent position and the position of the occurrence information onto the horizontal plane simultaneously cannot be located on the same straight line. b ≥2, e is a set parameter; In step 2.3.5, the calculated occurrence section satisfies the requirement that the sum of the squares of the vertical distances of all calculated spatial points from the calculated occurrence section is minimal.

Images