Geological modeling method and computer device based on the combination of geophysical exploration and drilling data

Abstract

The present invention relates to the field of geological modeling, and specifically to a geological modeling method and computer device based on the combination of geophysical and drilling data. The present invention uses point-like drilling data as a constraint condition for geophysical inversion iterative calculation, and uses the attribute model constructed by the drilling data as the initial model for geophysical inversion iterative calculation. After continuous iteration, a global optimal physical attribute model is obtained, and then the attribute model is searched through the physical attribute boundary value of the adjacent strata to obtain the interface, and then the adjacent interface is used to determine the geological classification of all grid points, and the conversion from the attribute model to the geological model is completed. This is used as a constraint condition and an initial model to achieve the fine modeling of the geological model of the texture model constructed by continuous geophysical data, thereby improving the modeling accuracy of the geological model and meeting the application requirements of digital scenes such as real three-dimensional scenes. The present invention is suitable for geological modeling of underground space.

Description

Technical Field

[0001] The present invention relates to the field of geological model building, and in particular to a geological model building method and a computer device based on the combination of geophysical exploration and drilling data. Background Art

[0002] With the continuous advancement of real-life 3D engineering construction and application, the modeling technology of the above-ground part is becoming more and more perfect, but there are still technical structural problems in the refined modeling technology of the underground geological model. Traditional geological modeling is generally modeled through drilling data, and drilling is point data, and interpolation is used between boreholes for modeling.

[0003] For example, CN117671182A discloses an improved three-dimensional geological modeling method based on Kriging interpolation, which includes: S1, obtaining the plane coordinates of the borehole and the elevations of each stratum interface corresponding to the plane coordinates; S2, determining the plane coordinates of the interpolation point, performing Kriging interpolation on the elevation of the same stratum interface, and obtaining the coordinate data of each stratum interface after interpolation; S3, generating a virtual interface, and obtaining the outer contour line of each virtual interface; S4, obtaining the bottom and top surfaces of the geological model; S5, generating a terrain entity between the bottom and top surfaces of the geological model; S6, segmenting the terrain entity based on the virtual interfaces of each stratum to obtain a layered terrain entity model; S7, setting parameters for each layer model in the layered terrain entity model to distinguish them.

[0004] However, the interpolation method used for modeling between boreholes will result in poor accuracy, especially when the boreholes are sparse and the distance between boreholes is large. The model will be seriously distorted and cannot meet the requirements of real-life 3D for the refinement of underground geological models. Summary of the invention

[0005] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a geological modeling method and a computer device based on the combination of geophysical exploration and drilling data, thereby improving the refinement of the underground geological model.

[0006] The present invention adopts the following technical solutions to achieve the above-mentioned purpose. In a first aspect, the present invention provides a geological modeling method based on the combination of geophysical exploration and drilling data, comprising:

[0007] S1. Obtain the physical property data of underground space by inversion calculation of the measured surface geophysical data;

[0008] S2, converting the physical attribute data of the underground space into the geological classification of the underground space;

[0009] S3. Construct a geological model based on the geological classification of the underground space.

[0010] Furthermore, step S1 specifically includes:

[0011] Collect the plane position data of each borehole in the survey area, the depth data of each borehole lithology layer, and the physical property data of each borehole lithology layer, and use the collected borehole data as the constraint input for participating in the geophysical inversion calculation. At the same time, use the physical property model constructed by the collected borehole data as the initial model for the geophysical inversion calculation;

[0012] Calculate the physical field of the initial model, and compare the physical field of the initial model with the measured data of geophysical exploration in the survey area and the physical property data of the lithology layer of each borehole, calculate the difference and make a judgment. If the difference is greater than the set threshold, adjust the initial model to obtain the process model;

[0013] Calculate the physical field of the process model, and compare the physical field of the process model with the measured geophysical data in the survey area and the physical property data of the lithology layers of each borehole, calculate the difference and make a judgment until the difference is less than the set threshold. Then, the process model at this time will be regarded as the global optimal physical property model.

[0014] Furthermore, step S2 specifically includes:

[0015] According to the resolution of the physical attribute model, the grid size is calculated, and the physical attribute value of each grid point is calculated using the Kriging method;

[0016] Starting from the boundary point of the borehole, according to the average values ​​of the physical properties of the upper and lower layers, search in different directions, take the grid point with the smallest difference with the average values ​​of the physical properties of the two layers as the boundary point, obtain the interface of the two lithologies, search different boundary points of different boreholes in turn, and finally obtain all the interfaces in the survey area;

[0017] Assign geological classification to the grid points of non-boundary points between adjacent interfaces, randomly select non-boundary points, query the geological classification of the boundary points in each direction with the smallest distance from the selected non-boundary points, if the boundary numbers of at least two boundary points in each direction are different, determine the geological classification of the non-boundary point according to the common geological classification of boundary points with different boundary numbers, and for the remaining closed or semi-closed non-boundary points that cannot be selected, determine their geological classification values ​​according to another geological classification of the boundary points on the interface;

[0018] After all the grid points are assigned geological classifications, the isolated points of the grid points are reassigned geological classifications, and the isolated points are assigned geological classifications of the adjacent grid points. The size of the isolated points is determined according to the resolution.

[0019] Furthermore, step S3 specifically includes:

[0020] Adjacent grid points of the same geological classification are rendered in a way that multiple grid points form a tetrahedron or hexahedron to form the smallest unit of the model and finally form a geological model.

[0021] In a second aspect, the present invention provides a computer device, including a memory, wherein the memory stores program instructions, and when the program instructions are executed, the above-mentioned geological modeling method based on the combination of geophysical exploration and drilling data is executed.

[0022] The beneficial effects of the present invention are:

[0023] The present invention uses point-like drilling data as constraints for geophysical inversion iterative calculations, and uses the attribute model constructed by the drilling data as the initial model for geophysical inversion iterative calculations. After continuous iterations, a linear and unbiased global optimal physical attribute model is finally obtained. The attribute model is then searched through the physical attribute boundary values ​​(average values) of adjacent strata to obtain the interface. Finally, the adjacent interface is used to determine the geological classification of all grid points, and the conversion from the attribute model to the geological model is completed, thereby realizing the construction of a skeleton model of the geological model through point-like drilling and trenching data, and using this as a constraint and initial model to realize the refined modeling technology of the geological model that uses continuous geophysical data to construct a texture model. Through this technical invention, the modeling accuracy of the geological model can be improved to meet the application requirements of digital scenes such as real three-dimensional scenes. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 It is a flow chart of a geological model building method based on the combination of geophysical exploration and drilling data provided by an embodiment of the present invention;

[0025] Figure 2 is a schematic diagram of a demarcation point provided by an embodiment of the present invention;

[0026] Figure 3 It is a schematic diagram of the interface of different lithologies between boreholes provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0027] To make the purpose, technical solution and advantages of the embodiments of the present invention more clear, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0028] The present invention uses drilling data to construct a skeleton model of the geological model, and the geophysical data of inter-hole continuity to construct a texture model. The skeleton model constructed by drilling is used as a constraint condition for geophysical inversion calculation, and the underground model is obtained by inversion calculation. The obtained underground model is then used for forward calculation to obtain the model physical property response. The model response is compared with the initial geophysical measured data. If the preset threshold is not met, the model is modified to calculate the model physical property response. After repeated iterative calculations, the preset threshold is reached. The underground model finally obtained by the inversion calculation is the global optimal underground model that meets the linear unbiased estimation, thereby realizing the refined modeling of the geological model.

[0029] Specifically, Figure 1 As shown, the present invention provides a geological modeling method based on the combination of geophysical exploration and drilling data, comprising:

[0030] S1. The physical property data of underground space is obtained by inversion calculation of the measured surface geophysical data.

[0031] Collect the plane position of each borehole in the survey area, the depth data of each borehole lithology layer, and the physical property data of each borehole lithology layer as the constraint input for participating in the geophysical inversion calculation. At the same time, the physical property model constructed by the borehole data is used as the initial model for the geophysical inversion calculation;

[0032] The physical field of the initial model is calculated by the finite difference method, finite element method or other simulation methods, and the physical field of the initial model is compared with the geophysical survey data in the survey area and the physical property data of the lithology layers of each borehole to calculate the difference;

[0033] The difference is judged. If the difference is greater than the set threshold, the initial model is adjusted to obtain the process model, the physical field of the process model is calculated, and the physical field of the process model is compared with the geophysical measured data in the survey area and the physical property data of the lithology stratification of each borehole. This step is repeated until the difference is less than the set threshold. At this time, the process model is the global optimal physical property model.

[0034] S2. Convert the physical property data of the underground space into the geological classification of the underground space.

[0035] According to the resolution of the physical attribute model, the grid size is calculated, and the physical attribute value of each grid point is calculated using the Kriging method;

[0036] Starting from the boundary point of the borehole, search in different directions with the average value of the physical properties of the upper and lower layers of lithology, and take the grid point closest to the average value of the physical properties of the two layers of lithology as the boundary point to obtain the interface between the two lithologies. Search different boundary points of different boreholes in turn, and finally obtain all the interfaces in the survey area;

[0037] Then, assign geological classifications to other grid points that are not boundary points between adjacent interfaces. Randomly select a non-boundary point and query the geological classifications of the six directional boundary points closest to it. If the boundary numbers of at least two directional boundary points are different, the geological classification of the non-boundary point can be determined based on the common geological classification of boundary points with different boundary numbers. For the remaining closed or semi-closed non-boundary points that cannot be selected, their geological classification values ​​can be determined based on another geological classification of the geological classification of the boundary point on the interface;

[0038] After all grid points are assigned geological classifications, the isolated points of the grid points are re-assigned to the geological classifications of the adjacent grid points. The size of the isolated points is determined by the resolution.

[0039] S3. Construct a geological model based on the geological classification of the underground space.

[0040] Adjacent grid points of the same geological classification are rendered in a manner that every eight grid points form a hexahedron or in a manner that every four grid points form a tetrahedron, to form the smallest unit of the model, and thus form a geological model.

[0041] The following takes electrical prospecting as an example to explain the solution of the present invention in detail.

[0042] like Figure 3 As shown in the figure, between the two boreholes ZK01 and ZK02, the lithology from top to bottom is limestone, slate and mudstone, and the resistivity is R1, R2 and R3 respectively. An electrical depth survey line WT1 is deployed, and the 12 measuring points are WT1-1 to WT1-12. The conventional modeling method is to use interpolation to build the model between boreholes. Figure 2 The dotted line is the interface between different lithologies between boreholes. In fact, the lithology interface between boreholes is a solid line.

[0043] Specifically, the inversion calculation process includes:

[0044] The underground space is divided into regular grids of a certain size, and each grid is regarded as a resistor. By assigning different resistivity values ​​to the grids, the resistivity property model of the underground space can be approximately simulated. When measuring by direct current method, the measuring device is divided into positive and negative electrodes AB at the transmitting end and positive and negative electrodes MN at the measuring end. The measured resistivity is equal to the voltage at the measuring end MN divided by the AB circuit at the transmitting end multiplied by the device coefficient to obtain the observed resistivity value, which is equivalent to the series-parallel resistance of the grid resistance within the electric field influence range under the measuring end MN (the electric field influence range is controlled by the length of the AB pole at the transmitting end. The longer the length, the deeper the electric field depth), that is, the total resistance of all grids within the box range between MN, multiplied by the unit grid length divided by the cross-sectional area to obtain the theoretical calculated resistivity value of the measuring end MN.

[0045] Establish the objective function, R i * is the observed resistivity data, R i is the theoretically calculated resistivity data, N is the number of geological model subdivision units, and the goal is to minimize the difference between the observed resistivity and the theoretically calculated resistivity. At the same time, it is necessary to ensure that the theoretically calculated resistivity at the locations x1 and x2 of boreholes ZK01 and ZK02 is consistent with the measured resistivity of the borehole lithology specimens, thereby reconstructing the resistivity attribute model of the underground space:

[0046]

[0047] The above equations are nonlinear problems and are generally solved by local optimization methods. Starting from the initial model, the global optimal resistivity attribute model is approached through multiple iterations. The initial model is based on the resistivity interpolation of different rock samples measured in boreholes ZK01 and ZK02 to obtain the resistivity distribution of each grid unit between the boreholes, that is, the resistivity of the grid unit above the dotted line L1 is R1, the resistivity of the grid unit between L1 and L2 is R2, and the resistivity of the grid unit below L2 is R3.

[0048] The iterative update form of the attribute model is as follows:

[0049] , m k is the attribute model at the kth iteration, and are the step size and update direction of the attribute model at the kth iteration respectively. When k=1, m1 is the initial model.

[0050] In order to find the optimal update direction, the objective function can be Perform a second-order Taylor expansion:

[0051]

[0052] Compare the two sides of the above equation Taking the derivative, we get:

[0053]

[0054] when When will reach an extreme value, when:

[0055]

[0056] Rewrite the above formula as: ,in, is the gradient of the objective function, is the Hessian matrix, then: .

[0057] Update the step size based on the solution , gradient G k The updated model can be obtained by calculating the objective function of the updated model and iterating continuously until the preset conditions are met. The final updated model is the global optimal underground space attribute model.

[0058] Specifically, the geological classification process includes:

[0059] Take the center point of the grid cell in the previous inversion calculation process, calculate the coordinates of the center point (X, Y, Z), and assign the resistivity value R of the grid cell to the corresponding center point to form a dot matrix, each point contains four values ​​(X, Y, Z, R).

[0060] According to the resolution of the final geological model, the appropriate model grid size and grid point coordinates (x, y, z) are calculated. The Kriging method is used to perform grid interpolation calculations on the point matrix of the previous step to obtain the resistivity value r of the grid point, thereby forming a new point matrix, each point containing (x, y, z, r).

[0061] Starting from the dividing point P0 on the lithologic interface L1 of borehole ZK01, the resistivity of P0 is the average value (R1+R2) / 2 of the lithologic resistivities R1 and R2 of the upper and lower layers of the interface L1. Searches are performed in three different directions, x, y, and z, and the grid point with the resistivity closest to (R1+R2) / 2 is selected as the dividing point.

[0062] For example, the resistivities of adjacent P1 and P2 in the x direction are R P1 and R P2 , if (R1+R2) / 2>R P1 >R P2 or (R1+R2) / 2 <R P1 <R P2 , then select P1 as the dividing point, if (R1+R2) / 2>R P2 >R P1 or (R1+R2) / 2 <R P2 <R P1 , then select P2 as the demarcation point. If R P1 >(R1+R2) / 2>R P2 or R P1 >(R1+R2) / 2>R P2 , then compare the resistivity difference between it and the two adjacent points, and take the smaller difference as the dividing point.

[0063] Compare the adjacent points in the y and z directions according to the above method, such as Figure 2As shown in the figure, we get the adjacent demarcation points P1, P3, and P6 in the three directions of P0 respectively, and mark whether they are demarcation points as "1" (1=yes, 0=no), mark the demarcation sequence as 1 (0=boundary, 1=L1, 2=L2), mark the geological classification as 1 and 2 (0=air / boundary, 1=limestone, 2=slate, 3=mudstone), and mark whether P0 has been searched as 1 (1=searched, 0=not searched). Then the grid points P0, P1, P3, and P6 are (x0, y 0,z0,r0,1,1,1 / 2,1), (x1,y1,z1,r1,1,1 / 2,0), (x3,y3,z3,r3,1,1,1 / 2,0), (x6,y6,z6,r6,1,1,1 / 2,0), and then search for the adjacent boundary points of P1, P3, and P6 respectively and mark them. No longer mark the grid points that are boundary points, and no longer search the grid points that have been searched. Do this and eventually search out all the boundary points of the boundary interface L1.

[0064] Similarly, all the boundary points of the lithologic interface L2 of borehole ZK01 can be obtained, and finally all the interfaces passing through borehole ZK01 can be obtained. Then, different boundary points of other boreholes (such as ZK02) are searched in turn.

[0065] For grid points on the boundary, the mark whether it is a boundary point is 1, the mark boundary sequence number is 0, the mark geological classification is 0, and the mark whether it has been searched is 1.

[0066] For the trap interfaces that do not pass through the boreholes, according to the resistivity values ​​of the interfaces of each borehole, the corresponding interfaces are searched and marked in turn according to the above method. Finally, all the interfaces in the survey area are obtained.

[0067] Randomly select a non-dividing point, query the geological classification of the six directional dividing points closest to it, at least two directional dividing points must have different boundary numbers (except for the boundary point with number 0), otherwise, randomly select a non-dividing point again. After selecting a non-dividing point (such as point Ps in the figure, the boundary numbers of its positive and negative Z directions are 1 and 2 respectively), the geological classification of the non-dividing point can be determined according to the common geological classification of the boundary points with different boundary numbers (the geological classification of the boundary point of point Ps in the positive Z direction is 1 / 2, the geological classification of the boundary point in the negative Z direction is 2 / 3, the common geological classification is 2, and the common geological classification is slate), and its geological classification is marked. Therefore, the geological classification of point Ps is 2 (slate), and all non-dividing points connected to it are marked with the same geological classification (stop marking when encountering a dividing point). For the remaining non-boundary points of the closed or semi-closed traps that cannot be selected, their geological classification values ​​can be determined based on another geological classification of the geological classification of the boundary points on the interface. For example, the geological classification of all connected non-boundary points in the area above the interface L1 is 1 (limestone), the geological classification of all connected non-boundary points in the area below the interface L2 is 3 (mudstone), and the geological classification of all connected non-boundary points in the closed area Dn below the interface L2 is 2 (limestone).

[0068] After all the grid points are assigned geological classifications, the isolated points of the grid point geological classification are reassigned and assigned to the geological classification of the adjacent grid points. The size of the isolated point is determined according to the resolution. For example, Dl in the figure is an area formed by one or several points, and the geological classification is 2 (slate), which is called an isolated point. It is generally caused by human error and external interference. Therefore, the boundary points and non-boundary points of the area need to be marked and reset, and re-marked as the geological classification 3 (mudstone) of the adjacent grid points outside the area.

[0069] The above is only a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the form disclosed herein, and should not be regarded as excluding other embodiments, but can be used in various other combinations, modifications and environments, and can be modified within the scope of the concept described herein through the above teachings or the technology or knowledge of the relevant field. The changes and modifications made by those skilled in the art shall not deviate from the spirit and scope of the present invention, and shall be within the scope of protection of the claims attached to the present invention.

Claims

1. A geological modeling method based on the combination of geophysical exploration and drilling data, characterized in that: include: S1. Obtain the physical property data of underground space by inversion calculation of the measured surface geophysical data; S2, converting the physical attribute data of the underground space into the geological classification of the underground space; S3. Construct a geological model based on the geological classification of underground space; Step S1 specifically includes: Collect the plane position data of each borehole in the survey area, the depth data of each borehole lithology layer, and the physical property data of each borehole lithology layer, and use the collected borehole data as the constraint input for participating in the geophysical inversion calculation. At the same time, use the physical property model constructed by the collected borehole data as the initial model for the geophysical inversion calculation; Calculate the physical field of the initial model, and compare the physical field of the initial model with the measured data of geophysical exploration in the survey area and the physical property data of the lithology layer of each borehole, calculate the difference and make a judgment. If the difference is greater than the set threshold, adjust the initial model to obtain the process model; Calculate the physical field of the process model, and compare the physical field of the process model with the geophysical measured data in the survey area and the physical property data of the lithology layers of each borehole, calculate the difference and make a judgment, until the difference is less than the set threshold, then the process model at this time is regarded as the global optimal physical property model; Step S2 specifically includes: According to the resolution of the physical attribute model, the grid size is calculated, and the physical attribute value of each grid point is calculated using the Kriging method; Starting from the boundary point of the borehole, according to the average values ​​of the physical properties of the upper and lower layers, search in different directions, take the grid point with the smallest difference with the average values ​​of the physical properties of the two layers as the boundary point, obtain the interface of the two lithologies, search different boundary points of different boreholes in turn, and finally obtain all the interfaces in the survey area; Assign geological classification to the grid points of non-boundary points between adjacent interfaces, randomly select non-boundary points, query the geological classification of the boundary points in each direction with the smallest distance from the selected non-boundary points, if the boundary numbers of at least two boundary points in each direction are different, determine the geological classification of the non-boundary point according to the common geological classification of boundary points with different boundary numbers, and for the remaining closed or semi-closed non-boundary points that cannot be selected, determine their geological classification values ​​according to another geological classification of the boundary points on the interface; After all the grid points are assigned geological classifications, the isolated points of the grid points are reassigned geological classifications, and the isolated points are assigned geological classifications of the adjacent grid points. The size of the isolated points is determined according to the resolution.

2. The geological modeling method based on the combination of geophysical exploration and drilling data according to claim 1 is characterized in that: Step S3 specifically includes: Adjacent grid points of the same geological classification are rendered in a way that multiple grid points form a tetrahedron or hexahedron to form the smallest unit of the model and finally form a geological model.

3. A computer device, comprising a memory, wherein the memory stores program instructions, characterized in that: When the program instructions are executed, the geological model building method based on the combination of geophysical exploration and drilling data as described in any one of claims 1-2 is executed.

Images