A formation interface modeling method based on improved Kriging interpolation

Abstract

The present invention discloses a method for modeling the formation interface based on improved Kriging interpolation, which relates to the field of information models. Specifically, it is used to fuse the attitude information of the formation interface during the Kriging interpolation process. The method includes: obtaining formation demarcation points and formation interface attitude information; introducing a penalty factor constructed from the formation interface attitude information into the variogram fitting process, thereby improving the determination of the variogram during the Kriging interpolation process; encrypting the grid, and then determining the encrypted points by the improved Kriging interpolation; connecting the encrypted points and the formation demarcation points according to the spatial grid connection relationship to generate the formation interface. The method of the present invention effectively improves the accuracy of the geological model by integrating two different types of geological information into the geological model through a penalty factor.

Description

Technical Field

[0001] The present invention relates to the field of information models, and particularly to a method for modeling formation interfaces based on improved Kriging interpolation. Background Art

[0002] Generating spatially continuous formation interfaces from limited formation boundary point data obtained by boreholes through spatial interpolation methods to estimate underground mineral reserves and distributions is a commonly used technical means in the exploration field. Among them, the Kriging interpolation method is widely used due to its relatively high estimation accuracy.

[0003] In existing geological modeling methods based on Kriging interpolation, only the formation boundary point information obtained by boreholes is considered, and the attitude information of the formation interfaces obtained by means such as mapping is not considered. For example, in the invention patent CN109191573A published in 2019, the formation interface skeleton is generated by parabolic interpolation, and then the formation interface is generated by Kriging interpolation, improving the modeling efficiency and quality; in the invention patent CN115422830A published in 2022, the SVM support vector machine algorithm is used to replace the commonly used theoretical variogram to improve the spatial interpolation accuracy.

[0004] Existing scholars have all improved Kriging interpolation through existing technologies to improve the modeling accuracy, but have not considered combining formation boundary point data with formation interface attitude data during the modeling process to improve the modeling accuracy. Summary of the Invention

[0005] In order to overcome the defects in the above technologies, the present invention provides a method for modeling formation interfaces based on improved Kriging interpolation. By introducing a penalty factor constructed from formation interface attitude information during the variogram fitting process, a fusion modeling of formation boundary point data and formation interface attitude data is realized.

[0006] To achieve the above object, the present invention adopts the following technical solutions.

[0007] A method for modeling formation interfaces based on improved Kriging interpolation includes the following steps:

[0008] Step 1, obtaining the spatial coordinate information of formation boundary points; obtaining the attitude information of the formation interfaces and the positions where the attitude information is located;

[0009] Step 2, calculating the semi-vertical distance difference and horizontal distance between two formation boundary points, constructing a penalty factor from the attitude information, and fitting to obtain the final theoretical variogram;

[0010] Step 3, setting the spatial interpolation encryption size, calculating the plane coordinates of the encrypted projection points according to the encryption size, and then determining theZ Axis coordinates are used to determine encrypted points, and at the same time, the spatial grid connection relationship obtained during the encryption process is acquired;

[0011] Step 4: Connect the encrypted points and the formation demarcation points according to the spatial grid connection relationship to generate a formation interface.

[0012] Preferably, in Step 1,

[0013] The spatial coordinates of the formation demarcation 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 Two coordinate axes are in the horizontal plane, and X The positive direction of the Z axis is to the right, Y The positive direction of the n axis is upward, i The positive direction of the x i axis is determined according to the left-hand rule; Suppose there are y i formation demarcation points in total, and the spatial coordinates of the z i th formation demarcation point are ( i n ),

[0014] The attitude information of the formation interface includes the actual dip and actual dip angle of the actual dip section at the position where the attitude information is located; The position where the attitude information is located refers to the plane coordinates of the position where the attitude information is located; Suppose there are m pieces of attitude information in total, and the plane coordinates of the position where the k th piece of attitude information is located are ( x k , y k ), and the actual dip and actual dip angle of the corresponding actual dip section are respectively recorded as φ k and θ k , k m .

[0015] Preferably, the specific process of Step 2 is as follows:

[0016] Step 2.1: Calculate the semi-vertical distance difference i j between any two formation demarcation points r ij ​​​, calculate the horizontal distance between any two formation demarcation points according to the horizontal distance calculation formula i and j of the horizontal distance d ij ; i and j = 1, 2,..., n , there are a total of n formation demarcation points;

[0017] Step 2.2, select a theoretical variogram r (·), the theoretical variogram r (·) contains two parameters to be fitted C 0, C 1; Denote C 0 as the first parameter to be fitted, C 1 as the second parameter to be fitted;

[0018] Introduce an objective function Obj :

[0019] ;

[0020] In the formula, r ( d ij ) is the horizontal distance d ij substituted into the theoretical variogram r (·) to calculate the semivariance, that is, the semivariance of the formation demarcation points i and j ; is the calculated dip of the calculated dip section at the position where the k th attitude information is located; is the calculated dip angle of the calculated dip section at the position where the k th attitude information is located; φ k is the actual dip of the actual dip section at the position where the k th attitude information is located; θ k is the actual dip angle of the actual dip section at the position where the k th attitude information is located; k = 1, 2,..., m , there are a total of m attitude information; H is the penalty factor, which is a function positively correlated with and ;

[0021] Step 2.3, set the first parameter to be fitted C 0 and the second parameter to be fitted CThe initial value of 1, calculate the objective function Obj value;

[0022] Step 2.4, adjust the first parameter to be fitted C 0 and the second parameter to be fitted C 1 value, solve to obtain the minimum value of the objective function Obj and denote it as the temporary minimum value;

[0023] Step 2.5, randomly set the first fitting parameter C 0 and the second parameter to be fitted C 1 initial value in the new round of calculation, repeat steps 2.3 - 2.4, and obtain the minimum value of the objective function Obj in the new round of calculation, and compare it with the temporary minimum value; if the minimum value of the objective function Obj in the new round of calculation is less than the temporary minimum value, then replace the temporary minimum value with the minimum value of the objective function in the new round of calculation; otherwise, do not replace;

[0024] Step 2.6, loop step 2.5, and terminate when the temporary minimum value has not been replaced for consecutive w times; the w is a set parameter; take the temporary minimum value obtained at termination as the final minimum value of the objective function Obj , and the theoretical variogram corresponding to the final minimum value of this objective function Obj is the final theoretical variogram.

[0025] Preferably, the specific process of step 3 is as follows:

[0026] Step 3.1, set the encryption size to a ; obtain the spatial coordinates of n stratum boundary points, and project them onto the horizontal plane to obtain n projection points, and denote these n projection points as the basic projection points;

[0027] Step 3.2, connect the basic projection points according to triangles, ensure that each basic projection point is the vertex of a triangle, and the triangle vertices are all basic projection points, that is, form a basic triangular grid;

[0028] Step 3.3, encrypt the basic triangular grid according to the principle of equal - dividing triangles until the average side length of each triangle is less than a ; denote the points other than the basic projection points among the grid vertices of the encrypted basic triangular grid as encrypted projection points; record the connection relationship between the basic projection points and the encrypted projection points according to the encrypted basic triangular grid, and denote it as the plane grid connection relationship;

[0029] Step 3.4, select an encrypted projection point g , obtain its horizontal coordinate ( x g , y g ), calculate the horizontal distance between this encrypted projection point g and the i th formation demarcation point, and denote it as the encrypted projection horizontal distance d ig , i = 1, 2,..., n ; then replace the horizontal distance between any two formation demarcation points with this encrypted projection horizontal distance d ig and substitute it into the final theoretical variogram d ij obtained in Step 2, (·), to obtain the semivariance between this encrypted projection point r and the g th formation demarcation point, and denote it as the projection semivariance i ( r ( d ig ));

[0030] Step 3.5, calculate the horizontal distances between the encrypted projection point g and each of the n formation demarcation points in the same way as in Step 3.4, to obtain n encrypted projection horizontal distances d ig , n projection semivariances r ( d ig ); according to the weight coefficient matrix, given n projection semivariances r ( d ig ) corresponding weight groups λ 1g , λ 2g ,..., λ ig ,..., λ ng , and substitute them into the weight formula to obtain the estimated value g of the Z axis coordinate of the encrypted projection point z g ; combine the estimated value Z of the z axis coordinate of the encrypted projection point g x g , yg ), combine to form encrypted points ( x g , y g , z g );

[0031] Repeat steps 3.4 - 3.5 until each encrypted projection point generates a corresponding encrypted point;

[0032] Step 3.6, in the planar grid connection relationship, replace the corresponding encrypted projection points with encrypted points and the formation demarcation points with corresponding basic projection points to form a spatial grid connection relationship.

[0033] Preferably, in step 2.3, calculate the value of the objective function Obj as follows:

[0034] Step 2.3.1, set the initial values of the parameters in the theoretical variogram r ();

[0035] Step 2.3.2, calculate the first half part d ij of the objective function r ij according to the horizontal distance Obj and the semi-vertical offset difference Obj sub1 ,

[0036] ;

[0037] Step 2.3.3, select a dip information, obtain the position of the k th dip information and the planar coordinates of its b adjacent positions, denoted as the dip peripheral positions. A total of b + 1 dip peripheral positions are obtained; calculate the horizontal distances between the dip peripheral positions and all formation demarcation points in turn, store the horizontal distances between each dip peripheral position and all formation demarcation points in a dip horizontal distance group, and obtain b + 1 dip horizontal distance groups. Each dip horizontal distance group contains n horizontal distances d ip , p is the number of the dip peripheral position, p = 1, 2,..., b + 1; substitute each horizontal distance in each dip horizontal distance group into the theoretical variogram r () to calculate the semi-variance, and store the results calculated by substituting each dip horizontal distance group into the theoretical variogram r () separately, and obtainb +1 set of attitude semivariograms, each set of attitude semivariograms contains n semivariograms r ( d ip );

[0038] Step 2.3.4, substitute each set of attitude semivariograms corresponding to the peripheral positions of each attitude into the weight coefficient matrix, and solve to obtain the weight sets of each attitude peripheral position λ 1p , λ 2p ,..., λ ip ,..., λ np , and then substitute the weight sets of each attitude peripheral position into the weight formula respectively to calculate the Z pre-estimated value of the z-axis coordinate z p ;

[0039] Step 2.3.5, combine the horizontal coordinates of each attitude peripheral position with the Z pre-estimated value of the z-axis coordinate z p to obtain a set of spatial points, denoted as the calculated spatial points, and fit a calculated attitude section plane with this set of calculated spatial points as the k calculated attitude section plane at the position where the

[0040] According to the determination method of the actual dip and actual inclination angle of the actual attitude section plane, obtain the calculated dip k and calculated inclination angle of the calculated attitude section plane at the position where the th attitude information is located;

[0041] Step 2.3.6, repeat Steps 2.3.3 - 2.3.5, obtain the calculated dips and calculated inclination angles corresponding to all attitude information, and then given a penalty factor H , calculate the value of the second half Obj of the objective function Obj sub2 ;

[0042] .

[0043] Preferably, the determination of the actual dip and actual inclination angle of the actual attitude section plane is as follows: Define a strike line and a dip line, the strike line is the intersection line of the actual attitude section plane and the horizontal plane, and the dip line is located within the actual attitude section plane, perpendicular to the strike line and along ZThe component in the axial direction points negatively; the actual trend refers to the direction indicated by the projection of the trend line on the horizontal plane; the actual dip angle refers to the angle between the trend line and the horizontal plane.

[0044] The actual occurrence section at the position where the occurrence information is located refers to the spatial section of the formation interface at the position where the occurrence information is located.

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

[0046] ;

[0047] .

[0048] Preferably, the weight coefficient matrix and the weight formula are as follows:

[0049] ;

[0050] ;

[0051] In the formula, is the Lagrange multiplier, and the subscript v represents the number, which is the number of the encrypted projection point g or the number of the position around the occurrence p ; r ij is the semi-vertical distance difference between any two formation boundary points calculated according to the semi-vertical distance difference calculation formula i , j ; z i is the i th Z axis coordinate of the formation boundary point.

[0052] Preferably, the specific process of step 2.4 is as follows:

[0053] Step 2.4.1, taking the first parameter to be fitted C 0 and the second parameter to be fitted C 1 as the two directions of the coordinate axes, establish a plane coordinate system;

[0054] Step 2.4.2, set the initial values of the first parameter to be fitted C 0 and the second parameter to be fitted C 1 as the parameter points of this plane coordinate system;

[0055] Step 2.4.3, move the parameter points in different directions in the plane coordinate system according to a fixed length to obtain different C 0, CFor the trial parameter points of the 1 value, calculate the values of all objective functions, and select the trial parameter point with the minimum objective function value as the new parameter point until the minimum value of the objective function is obtained.

[0056] Preferably, in step 2.3.3, the b adjacent positions at the location where the attitude information is located satisfy: the horizontal distance between each adjacent position and the location where the attitude information is located is less than e , and b the points obtained by projecting the adjacent positions at the location where the attitude information is located and the location where the attitude information is located onto the horizontal plane cannot be 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 from all calculation space points to the calculated attitude section is the smallest.

[0057] The beneficial effects of the present invention are as follows:

[0058] (1) By introducing a penalty factor constructed from the attitude information of the formation interface as a constraint in the fitting process of the variogram, the fusion modeling of the formation boundary point data and the formation interface attitude data in the Kriging interpolation process is realized, effectively improving the modeling accuracy.

[0059] (2) In the fitting process of 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, ensuring the accuracy of the variogram fitting.

[0060] (3) By introducing the form of a penalty factor, the fusion of two different types of data, namely the formation boundary point and the formation interface attitude, is realized, and the algorithm is simple to implement, providing a reference idea for introducing other types of data in the geological space interpolation process. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 is a flowchart of the formation interface modeling method of the present invention.

[0062] Figure 2 is a schematic diagram of plane grid encryption.

[0063] Figure 3 is a schematic diagram of the actual attitude section and the calculated attitude section.

[0064] Figure 4 is a schematic diagram of the movement of parameter points in the process of solving the minimum value of the objective function.

[0065] Figure 5 is a simplified diagram of the formation interface modeling method of the present invention.

[0066] The reference numerals are explained as follows:

[0067] 1 - Borehole, 2 - Formation boundary point, 3 - Foundation projection point, 4 - Actual occurrence section, 5 - Density projection point, 6 - Y Axis, 7 - X Axis, 8 - Calculated occurrence 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. Specific implementation mode

[0068] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention.

[0069] Figure 1 It is a flowchart of the formation interface modeling method of the present invention, Figure 5 It is a sketch of the formation interface modeling method of the present invention. As can be seen from Figure 1 And Figure 5 It can be seen that the present invention provides a formation interface modeling method based on improved Kriging interpolation. The formation interface refers to the interface between two different formations, including the following steps:

[0070] Step 1, obtain the spatial coordinate information of the formation boundary points; obtain the occurrence information of the formation interface and the position where the occurrence information is located.

[0071] The spatial coordinates of the formation 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 Two coordinate axes are in the horizontal plane, and X The positive direction of the axis is to the right, Z The positive direction of the axis is upward, Y The positive direction of the axis is determined according to the left-hand rule; assume that there are Formation boundary points, and any one of them is recorded as the i th formation boundary point, i Is the number of the formation boundary point, i = 1, 2,..., n ; The spatial coordinates of the i th formation boundary point are ( x i , y i , z i ).

[0072] The attitude information of the formation interface includes the actual dip direction and actual dip angle of the actual dip section at the position where the attitude information is located; the position where the attitude information is located refers to the planar coordinates of the position where the attitude information is located; assume there are m pieces of attitude information, and any one of them is denoted as the k th piece of attitude information. The planar coordinates of the position where the k th piece of attitude information is located are ([[]] x k , y k ), and the actual dip direction and actual dip angle of the corresponding actual dip section are respectively denoted as φ k and θ k , k = 1, 2,..., m .

[0073] In this embodiment, the determination of the actual dip direction and actual dip angle of the actual dip section is as follows: The concepts of the strike line and dip line are introduced. The strike line is the intersection line of the actual dip section and the horizontal plane. The dip line is located within the actual dip section, perpendicular to the strike line, and the component along the Z axis points in the negative direction; the actual dip direction refers to the direction indicated by the projection of the dip line on the horizontal plane, and the actual dip angle refers to the angle between the dip line and the horizontal plane.

[0074] In this embodiment, the actual dip section at the position where the attitude information is located refers to the spatial section of the formation interface at the position where the attitude information is located.

[0075] As Figure 2 shown, eight formation boundary points of a certain formation interface are obtained through drilling, and there is one formation interface attitude and its horizontal coordinate.

[0076] As Figure 3 shown, there is one formation interface attitude information in the spatial coordinate system. This actual dip section can be uniquely determined by the actual dip line and actual strike line. Among them, the recording method of the actual dip direction is the angle formed by rotating the X axis clockwise to the projection line of the actual dip line, and 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 dip section, and the component of the actual dip line along the Z axis points in the negative direction.

[0077] In Figure 2 and Figure 3In 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.

[0078] 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.

[0079] 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 , a total of The half vertical difference and Horizontal distance; i , j =1, 2, ..., n .

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

[0081] ;

[0082] .

[0083] 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.

[0084] 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.

[0085] Introducing an objective function Obj:

[0086] ;

[0087] In the formula, r ( d ij ) is the horizontal distance d ij Substitute into the theoretical variogram r ( d ) to calculate the semi-variance obtained, that is, the semi-variance of the formation boundary points i , j ; is the calculated dip of the calculated dip section at the position where the k th attitude information is located; is the calculated dip angle of the calculated dip section at the position where the k th attitude information is located; φ k is the actual dip of the actual dip section at the position where the k th attitude information is located; θ k is the actual dip angle of the actual dip section at the position where the k th attitude information is located; k = 1, 2,... m ; H is the penalty factor, which is a function positively correlated with and .

[0088] In this embodiment, the exponential function is selected as the theoretical variogram. The independent variable parameter d here is the horizontal distance between two formation boundary points.

[0089] Step 2.3, set the initial values of the first parameter to be fitted C 0 and the second parameter to be fitted C 1, and calculate the value of the objective function Obj .

[0090] Step 2.4, adjust the values of the first parameter to be fitted C 0 and the second parameter to be fitted C 1 according to the gradient descent method, solve for the minimum value of the objective function Obj , and record it as the temporary minimum value.

[0091] Step 2.5, randomly set the initial values of the first fitting parameter C 0 and the second parameter to be fitted C 1 in the new round of calculation, repeat steps 2.3 - 2.4, and obtain the objective function ObjFind the minimum value and compare it with the temporary minimum value; if the minimum value of the objective function in the new round of calculation Obj is less than the temporary minimum value, then replace the temporary minimum value with the minimum value of the objective function in the new round of calculation; otherwise, do not replace.

[0092] Step 2.6, loop Step 2.5 and terminate when the temporary minimum value has not been replaced for w consecutive times; the w is a manually set parameter; take the temporary minimum value obtained at termination as the final minimum value of the objective function Obj , and the theoretical variogram corresponding to the final minimum value of this objective function Obj is the final theoretical variogram.

[0093] Step 3, set the spatial interpolation encryption size, calculate the plane coordinates of the encrypted projection points according to the encryption size, and then determine the Z axis coordinates of the encrypted projection points by the spatial interpolation algorithm, so as to determine the encrypted points, and at the same time obtain the spatial grid connection relationship obtained during the encryption process.

[0094] Step 3.1, set the encryption size to a ; obtain the spatial coordinates of n stratum boundary points, project them onto the horizontal plane to obtain n projection points, and record these n projection points as the basic projection points.

[0095] Step 3.2, connect the basic projection points according to triangles, ensuring that each basic projection point is the vertex of a triangle, and the vertices of the triangle are all basic projection points, that is, form a basic triangular grid.

[0096] As Figure 2 shown, each stratum boundary point projected onto the horizontal plane can obtain eight basic projection points, and then connect the eight basic projection points through a triangular grid to form a triangular grid.

[0097] Step 3.3, encrypt the basic triangular grid according to the principle of equal division of triangles until the average side length of each triangle is less than a ; record all points other than the basic projection points among the grid vertices of the encrypted basic triangular grid as encrypted projection points; record the connection relationship between the basic projection points and the encrypted projection points according to the encrypted basic triangular grid and denote it as the plane grid connection relationship.

[0098] In this embodiment, the principle of equal division of triangles means connecting the midpoint of the longest grid line segment of the triangular grid to all the triangle vertices corresponding to the longest grid line segment to divide the existing triangles; the triangle vertices corresponding to the longest grid line segment refer to the vertices corresponding to this grid line segment when it is used as the side of a triangle.

[0099] Step 3.4, select an encrypted projection point g , and obtain its horizontal coordinate ( x g , y g ), calculate the horizontal distance between this encrypted projection point g and the i th formation demarcation point, and denote it as the encrypted projection horizontal distance d ig ,

[0100] ;

[0101] Then replace the horizontal distance between any two formation demarcation points with this encrypted projection horizontal distance d ig and substitute it into the final theoretical variogram obtained in Step 2 to obtain the semi-variance between this encrypted projection point d ij and the g th formation demarcation point, and denote it as the projection semi-variance i ( r ( d ig ).

[0102] Step 3.5, calculate the horizontal distances between the encrypted projection point g and the n formation demarcation points respectively in the way of Step 3.4 to obtain n encrypted projection horizontal distances d ig , n projection semi-variances r ( d ig ); according to the weight coefficient matrix, given n projection semi-variances r ( d ig ), corresponding weight groups λ 1g , λ 2g ,..., λ ig ,..., λ ng , and substitute them into the weight formula to obtain the g axis coordinate estimated value Z of the encrypted projection point z g ; the Z axis coordinate estimated value z g of the encrypted projection point and the horizontal coordinate ( xg , y g ), combined to form encrypted points ( x g , y g , z g ).

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

[0104] Step 3.6, in the planar grid connection relationship, replace the corresponding encrypted projection points with encrypted points and the formation boundary points with corresponding basic projection points to form a spatial grid connection relationship.

[0105] Step 4, obtain the encrypted points, spatial grid connection relationship, and formation boundary points, and connect the encrypted points and formation boundary points according to the spatial grid connection relationship to generate the formation interface.

[0106] So far, the modeling of the formation interface is completed.

[0107] In this embodiment, the calculation of the objective function Obj in step 2.3 is as follows:

[0108] Step 2.3.1, set the initial values of the parameters in the theoretical variogram r ( d ).

[0109] Step 2.3.2, calculate the first half d ij , half vertical distance difference r ij of the objective function Obj , Obj sub1 value.

[0110] .

[0111] Step 2.3.3, select a dip information, obtain the position of the k th dip information and the planar coordinates of its b adjacent positions, denoted as the dip surrounding positions. A total of b + 1 dip surrounding positions are obtained; calculate the horizontal distances between the dip surrounding positions and all formation boundary points in turn, and store the horizontal distances between each dip surrounding position and all formation boundary points in a dip horizontal distance group to obtain b + 1 dip horizontal distance groups, and each dip horizontal distance group contains n horizontal distances d ip, p is the number of the positions around the attitude, p = 1, 2, ..., b + 1, and each group of horizontal distances of the attitude corresponds to the positions around the attitude participating in the calculation of this group; successively substitute each horizontal distance in each group of horizontal distances of the attitude into the theoretical variogram r ( d ) to calculate the semi-variance, and substitute each group of horizontal distances of the attitude into the theoretical variogram r ( d ) separately store the calculated results, and obtain b + 1 groups of attitude semi-variances, and each group of attitude semi-variances has n semi-variances r ( d ip ), and each group of attitude semi-variances corresponds to the positions around the attitude corresponding to the group of horizontal distances of the attitude participating in the calculation of this group.

[0112] In this embodiment, the adjacent positions at the b where the attitude information is located satisfy: the horizontal distance between each adjacent position and the position where the attitude information is located is less than e , and b the points obtained by projecting the adjacent positions at the b and the position where the attitude information is located onto the horizontal plane cannot be on the same straight line; the e ≥ 2; the

[0113] λ 1p 2p , λ 2p , ..., λ ip , ..., λ np , and then substitute the weight groups of each position around the attitude into the weight formula respectively to calculate the Z axial coordinate estimated value z p of each position around the attitude.

[0114] Z z axial coordinate estimated value p p to obtain a set of spatial points, denoted as the calculated spatial points, and fit a calculated attitude section with this set of calculated spatial points as the calculated attitude section at the position where the k th attitude information is located.

[0115] According to the method for determining the actual dip direction and actual dip angle of the actual occurrence section described in Step 1, the calculated dip direction of the calculated occurrence section is obtained. and the calculated dip angle ; the calculated occurrence section satisfies that the sum of the squares of the vertical distances from all calculated spatial points to the calculated occurrence section is the smallest.

[0116] Step 2.3.6, repeat Steps 2.3.3 - 2.3.5, obtain the calculated dip direction and calculated dip angle corresponding to all occurrence information, and then given a penalty factor H , calculate the latter half Obj of the objective function Obj sub2 value.

[0117] .

[0118] In this embodiment, in Step 2.4, adjusting the values of the first parameter to be fitted C 0 and the second parameter to be fitted C 1 according to the gradient descent method, and solving for the minimum value of the objective function Obj is as follows:

[0119] Step 2.4.1, taking the two directions of the first parameter to be fitted C 0 and the second parameter to be fitted C 1 as the two axes directions, establish a plane coordinate system.

[0120] Step 2.4.2, set the initial values of the first parameter to be fitted C 0 and the second parameter to be fitted C 1 as the parameter points of this plane coordinate system.

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

[0122] Figure 4 is the schematic diagram of the movement of the parameter points C 0, C 1 during the solution process of the minimum value of the objective function. Figure 4 In the figure, the reference numeral 17 is the parameter point.

[0123] In this embodiment, the weight coefficient matrix and the weight formula are as follows respectively:

[0124] ;

[0125] ;

[0126] In the formula, is the Lagrange multiplier, and the subscript v represents the number, which is the number of the encrypted projection points g or the number of the positions around the occurrence p .

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

Claims

1. A method for modeling the formation interface based on improved Kriging interpolation, characterized in that, It includes the following steps: Step 1, obtain the spatial coordinate information of the formation demarcation points; obtain the attitude information of the formation interface and the position where the attitude information is located; Step 2, calculate the semi-vertical distance difference and horizontal distance between each pair of formation demarcation points, construct a penalty factor from the attitude information, and fit to obtain the final theoretical variogram; Step 3, set the spatial interpolation encryption size, calculate the plane coordinates of the encrypted projection points according to the encryption size, and then determine the Z-axis coordinates of the encrypted projection points by the spatial interpolation algorithm, so as to determine the encrypted points, and at the same time obtain the spatial grid connection relationship obtained during the encryption process; Step 4, connect the encrypted points and the formation demarcation points according to the spatial grid connection relationship to generate the formation interface; The specific process of Step 2 is as follows: Step 2.1, calculate the half vertical distance difference r between any two formation boundary points i and j according to the half vertical distance difference calculation formula ij , calculate the horizontal distance d between any two formation boundary points i and j according to the horizontal distance calculation formula ij ; i, j = 1, 2,..., n, and there are n formation boundary points in total; Step 2.2, select a theoretical variogram r(·), and this theoretical variogram r(·) contains two parameters C0 and C1 to be fitted; record C0 as the first parameter to be fitted, and C1 as the second parameter to be fitted; Introduce an objective function Obj: where r(d ij ) is the semi-variance calculated by substituting the horizontal distance d ij into the theoretical variogram r(·), that is, the semi-variance of the stratigraphic demarcation points i and j; is the calculated dip of the calculated attitude section at the position where the k-th attitude information is located; is the calculated dip angle of the calculated attitude section at the position where the k-th attitude information is located; is the actual dip of the actual attitude section at the position where the k-th attitude information is located; θ k is the actual dip angle of the actual attitude section at the position where the k-th attitude information is located; k = 1, 2,..., m, and there are m attitude information in total; H is a penalty factor, which is a function and that is positively correlated; Step 2.3, set the initial values of the first parameter to be fitted C0 and the second parameter to be fitted C1, and calculate the value of the objective function Obj; Step 2.4, adjust the values of the first parameter to be fitted C0 and the second parameter to be fitted C1 according to the gradient descent method, solve for the minimum value of the objective function Obj, and record it as the temporary minimum value; Step 2.5, randomly set the initial values of the first fitting parameter C0 and the second fitting parameter C1 in the new round of calculation, repeat Steps 2.3 - 2.4, obtain the minimum value of the objective function Obj in the new round of calculation, and compare it with the temporary minimum value; if the minimum value of the objective function Obj in the new round of calculation is less than the temporary minimum value, then replace the temporary minimum value with the minimum value of the objective function Obj in the new round of calculation; Otherwise, do not replace; Step 2.6, loop Step 2.5, and terminate when the temporary minimum value is not replaced continuously for w times; the w is a set parameter; take the temporary minimum value obtained at the termination as the final minimum value of the objective function Obj, and the theoretical variogram corresponding to the final minimum value of this objective function Obj is the final theoretical variogram.

2. The formation interface modeling method based on improved Kriging interpolation according to claim 1, wherein, In Step 1, The spatial coordinates of the formation boundary points are recorded in the form of a spatial rectangular coordinate system, and the recorded values are the corresponding values in the directions of the three coordinate axes X, Y, and Z. Among them, the two coordinate axes X and Y are located in the horizontal plane, the positive direction of the X-axis is to the right, the positive direction of the Z-axis is upward, and the positive direction of the Y-axis is determined according to the left-hand rule. Suppose there are n formation boundary points in total, and the spatial coordinates of the i-th formation boundary point are (x i , y i , z i ), where i = 1, 2,..., n; The attitude information of the formation interface includes the actual dip direction and actual dip angle of the actual attitude section at the position where the attitude information is located; the position where the attitude information is located refers to the planar coordinates of the position where the attitude information is located; assume there are m pieces of attitude information, and the planar coordinates of the position of the k-th attitude information are (x k , y k ), and the actual dip direction and actual dip angle of the corresponding actual attitude section are respectively denoted as and θ k , where k = 1, 2,..., m.

3. A method for modeling the formation 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 as a; obtain the spatial coordinates of n formation demarcation points, project them onto the horizontal plane to obtain n projection points, and record these n projection points as the basic projection points; Step 3.2, connect the basic projection points according to triangles, ensure that each basic projection point is the vertex of a triangle, and the triangle vertices are all basic projection points, that is, form a basic triangular grid; Step 3.3, encrypt the basic triangular grid according to the principle of equal division of triangles until the average side length of each triangle is less than a; record all the points other than the basic projection points among the grid vertices of the encrypted basic triangular grid as encrypted projection points; record the connection relationship between the basic projection points and the encrypted projection points according to the encrypted basic triangular grid, and record it as the plane grid connection relationship; Step 3.4, select an encrypted projection point g, and obtain its horizontal coordinates (x g , y g ). Calculate the horizontal distance between the encrypted projection point g and the i-th formation boundary point, and denote it as the encrypted projection horizontal distance d ig , where i = 1, 2,..., n; then use the encrypted projection horizontal distance d ig to replace the horizontal distance d ij between any two formation boundary points, and substitute it into the final theoretical variogram r(·) obtained in Step 2 to obtain the semivariance between the encrypted projection point g and the i-th formation boundary point, and denote it as the projection semivariance r(d ig ); Step 3.5: Calculate the horizontal distances between the encrypted projection point g and the n formation demarcation points in the same way as in Step 3.4 to obtain n encrypted projection horizontal distances d ig , and n projection semi-variances r(d ig ); According to the weight coefficient matrix, given the weight groups λ ig corresponding to the n projection semi-variances r(d 1g ), λ 2g ,..., λ ig ,..., λ ng , and substitute them into the weight formula to obtain the estimated value z g of the Z-axis coordinate of the encrypted projection point g; Combine the estimated value z g of the Z-axis coordinate of the encrypted projection point with the horizontal coordinates (x g , y g ) to form an encrypted 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 planar grid connection relationship, replace the corresponding encrypted projection points with encrypted points and replace the corresponding basic projection points with formation demarcation points to form a spatial grid connection relationship.

4. A method for modeling the formation interface based on improved Kriging interpolation according to claim 1, characterized in that, In Step 2.3, calculate the value of the objective function Obj. The specific process is as follows: Step 2.3.

1. Set the initial values of the parameters in the theoretical variogram r(·). Step 2.3.2, according to the horizontal distance d ij , the semi-vertical distance difference r ij calculate the first half part Obj of the objective function Obj sub1 value, Obj sub1 = ∑[r(d ij ) - r ij 2 ;​ Step 2.3.

3. Select a dip information, obtain the planar coordinates of the position where the k-th dip information is located and its b adjacent positions at b, denoted as the positions around the dip, and a total of b + 1 positions around the dip are obtained. Calculate the horizontal distances between the positions around the attitude and all the formation demarcation points in sequence, store the horizontal distances between each position around the attitude and all the formation demarcation points in a group of attitude horizontal distances, and obtain b + 1 groups of attitude horizontal distances. There are n horizontal distances d in each group of attitude horizontal distances ip , where p is the number of the position around the attitude, p = 1, 2,..., b + 1; Substitute each horizontal distance in each group of attitude horizontal distances into the theoretical variogram r(·) to calculate the semi-variance in sequence, and store the results calculated by substituting each group of attitude horizontal distances into the theoretical variogram r(·) separately, obtaining b + 1 groups of attitude semi-variances. There are n semi-variances r(d ip ); Step 2.3.4, substitute the semi-variance groups of attitudes corresponding to the positions around each attitude into the weight coefficient matrix respectively, and solve to obtain the weight groups λ 1p , λ 2p ,..., λ ip ,..., λ np , and then substitute the weight groups of each position around the attitude into the weight formula respectively to calculate the estimated value z of the Z-axis coordinate of each position around the attitude p ; Step 2.3.5, combine the horizontal coordinates of the positions around each attitude with the estimated value z of the Z-axis coordinate p to obtain a set of spatial points, denoted as the calculated spatial points, and fit a calculated attitude section with this set of calculated spatial points as the calculated attitude section at the position where the k-th attitude information is located; According to the determination method of the actual dip direction and actual dip angle of the actual occurrence section, the calculated dip direction of the calculated occurrence section at the position where the k-th occurrence information is located is obtained and the calculated dip angle Step 2.3.6, repeat Steps 2.3.3 - 2.3.5 to obtain the calculated dip direction and calculated dip angle corresponding to all occurrence information, and then given a penalty factor H, calculate the latter half Obj of the objective function Obj sub2 value, 5. A method for modeling the formation interface based on improved Kriging interpolation according to any one of claims 1-4, characterized in that, The determination of the actual dip and actual dip angle of the actual dip section is as follows: Define the strike line and dip line. The strike line is the intersection line of the actual dip section and the horizontal plane. The dip line is located within the actual dip section, perpendicular to the strike line and the component along the Z-axis direction points downward; the actual dip refers to the direction indicated by the projection of the dip line on the horizontal plane; the actual dip angle refers to the angle between the dip line and the horizontal plane. The actual dip section at the position where the dip information is located refers to the spatial section of the formation interface at the position where the dip information is located.

6. A method for modeling the formation interface based on improved Kriging interpolation according to any one of claims 1-4, characterized in that The calculation formulas for the semi-vertical distance difference and the horizontal distance are as follows:

7. A method for modeling the formation interface based on improved Kriging interpolation according to claim 3 or 4, characterized in that The weight coefficient matrix and the weight formula are as follows: Wherein, φ is the Lagrange multiplier, the subscript v represents the number, which is the number g of the encrypted projection point or the number p of the position around the attitude; r ij is the semi-vertical distance difference between any two formation boundary points i and j calculated according to the semi-vertical distance difference calculation formula; z i is the Z-axis coordinate of the i-th formation boundary point.

8. A method for modeling the formation interface based on improved Kriging interpolation according to claim 1, characterized in that The specific process of Step 2.4 is as follows: Step 2.4.

1. Establish a planar coordinate system with two directions using the first parameter to be fitted C0 and the second parameter to be fitted C1 as the coordinate axes. Step 2.4.

2. Set the initial values of the first parameter to be fitted C0 and the second parameter to be fitted C1 as the parameter points of this planar coordinate system. Step 2.4.

3. Move the parameter points in different directions at a fixed length within the planar coordinate system to obtain trial parameter points with different C0 and C1 values, calculate the values of all objective functions, and select the trial parameter point with the minimum objective function value as the new parameter point until the minimum value of the objective function is obtained.

9. A method for modeling the formation interface based on improved Kriging interpolation according to claim 4, wherein, In Step 2.3.3, the b adjacent positions at the position where the dip information is located satisfy: the horizontal distance between each adjacent position and the position where the dip information is located is less than e, and the points obtained by projecting the b adjacent positions and the position where the dip information is located onto the horizontal plane cannot be on the same straight line, b ≥ 2, and e is a set parameter. In Step 2.3.5, the calculated dip section satisfies that the sum of the squares of the vertical distances from all calculated spatial points to the dip section is the smallest.

Images