A method for processing electrical prospecting data
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- HUA HUI ENGINEERING DESIGN GROUP CO LTD
- Filing Date
- 2023-09-18
- Publication Date
- 2026-05-12
Smart Images
Figure CN117148459B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of data processing technology, and specifically relates to a method for processing electrical exploration data. Background Technology
[0002] my country is currently in a period of rapid economic development, and its demand for resources and energy is increasing day by day. Therefore, it is urgent to conduct further exploration to find resources and energy supplies for the sustainable development of the national economy.
[0003] Electrical resistivity tomography (EPT) is an important geophysical exploration technique used to obtain information about the properties of subsurface rocks and soils. The main principle of existing EPT is that electric current is transmitted through the subsurface, and the voltage gradient is measured to estimate the resistivity distribution of the subsurface material. This relies on electrode arrangement, current injection, and voltage measurement, and then empirical formulas or forward and inverse models are used to estimate the resistivity distribution.
[0004] However, existing electrical exploration methods have some limitations and challenges.
[0005] First, the amount of data that can be actually obtained in the field is limited, and the number of unknown model parameters is often far greater than the amount of data, making it difficult to obtain an accurate resistivity distribution during forward and inverse modeling.
[0006] Secondly, when further analyzing the exploration results, it usually relies on experienced experts to conduct some professional assessments based on resistivity distribution maps. This includes assessing the degree of damage to the target rock. The assessment results are easily affected by subjective factors such as the expert's knowledge level, work experience, and emotional state, resulting in poor consistency and accuracy of the assessment results. Summary of the Invention
[0007] To address the limitations of existing technologies where the amount of data obtained from actual field exploration is limited, and the number of unknown model parameters often far exceeds the data volume, making it difficult to obtain accurate resistivity distributions during forward and inverse modeling, further analysis of exploration results typically relies on experienced experts to conduct professional assessments based on resistivity distribution maps. This includes assessing the degree of damage to the target rock. However, the assessment results are easily influenced by subjective factors such as the expert's knowledge level, work experience, and emotional state, leading to poor consistency and accuracy of the assessment results. This invention provides a method for processing electrical exploration data.
[0008] This invention provides a method for processing electrical resistivity tomography (OTT) data, comprising:
[0009] S101: Distributed electrodes are arranged around the target rock;
[0010] S102: Obtain the voltage value of each electrode and use the voltage value of each electrode as a boundary condition;
[0011] S103: Determine the Laplace equations that satisfy the boundary conditions for each point within the target rock.
[0012] S104: Solve for the voltage distribution within the target rock using the finite element method;
[0013] S105: Based on the voltage distribution within the target rock, the resistivity distribution within the target rock is calculated by inversion;
[0014] S106: Calculate the porosity of the target rock based on the resistivity distribution within the target rock;
[0015] S107: Determine the degree of damage to the target rock based on the porosity within the target rock.
[0016] Compared with the prior art, the present invention has at least the following beneficial technical effects:
[0017] (1) In this invention, in the face of limited data obtained from actual field exploration, the Laplace equation is determined by using the actual measured voltage value as the boundary condition. Then, the Laplace equation is solved by the finite element method, and the voltage distribution and resistivity distribution problems are decomposed into discrete elements. The voltage distribution is more accurately solved by using the limited data obtained from actual measurement. Then, the resistivity distribution in the target rock is calculated by inversion based on the voltage distribution in the target rock, thereby improving the accuracy of resistivity distribution calculation.
[0018] (2) In this invention, the porosity of the target rock is automatically calculated based on the resistivity distribution within the target rock, and then the degree of damage to the target rock is automatically determined based on the porosity. This eliminates the need for human intervention, reduces the influence of subjective factors, and improves the consistency and accuracy of rock damage assessment. Attached Figure Description
[0019] The preferred embodiments will now be described in a clear and easy-to-understand manner, in conjunction with the accompanying drawings, to further explain the above-mentioned characteristics, technical features, advantages, and implementation methods of the present invention.
[0020] Figure 1 This is a flowchart illustrating a method for processing electrical exploration data provided by the present invention;
[0021] Figure 2 This is a schematic diagram of the structure of an electrical exploration data processing method provided by the present invention;
[0022] Figure 3 This is a schematic diagram of a contour map provided by the present invention. Detailed Implementation
[0023] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the specific implementation methods of the present invention will be described below with reference to the accompanying drawings. Obviously, the drawings described below are merely some embodiments of the present invention. For those skilled in the art, other drawings and other implementation methods can be obtained based on these drawings without any creative effort.
[0024] To keep the drawings concise, each figure only schematically shows the parts relevant to the invention, and these do not represent the actual structure of the product. Furthermore, to facilitate understanding, in some figures, only one of components with the same structure or function is schematically depicted, or only one is labeled. In this document, "one" not only means "only one," but can also mean "more than one."
[0025] It should also be further understood that the term "and / or" as used in this specification and the appended claims refers to any combination of one or more of the associated listed items and all possible combinations, and includes such combinations.
[0026] In this document, it should be noted that, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to fixed connections, detachable connections, or integral connections. They can refer to mechanical connections or electrical connections. They can refer to direct connections or indirect connections through an intermediate medium, or internal connections between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.
[0027] Furthermore, in the description of this invention, the terms "first," "second," etc., are used only for distinguishing descriptions and should not be construed as indicating or implying relative importance.
[0028] In one embodiment, refer to the appendix to the specification. Figure 1 This diagram illustrates a flowchart of the electrical resistivity tomography data processing method provided by the present invention. (See attached specification.) Figure 2 The diagram shows a schematic representation of the electrical exploration data processing method provided by the present invention.
[0029] This invention provides a method for processing electrical resistivity tomography (OTT) exploration data, comprising:
[0030] S101: Arrange distributed electrodes around the target rock.
[0031] In one possible implementation, S101 specifically involves: arranging distributed electrodes at different thicknesses along the thickness direction of the target rock.
[0032] It should be noted that if the target rock has a large vertical thickness, a single horizontal electrode arrangement may be limited by depth. By arranging electrodes at different thicknesses, more dimensional resistivity data can be obtained, thus providing information on the resistivity distribution at different depths or thicknesses, rather than just a horizontal profile underground, which helps to understand the vertical resistivity distribution of the target rock.
[0033] S102: Obtain the voltage value of each electrode and use the voltage value of each electrode as the boundary condition.
[0034] The voltage values of each electrode can also be referred to as edge voltage values.
[0035] Boundary conditions, in mathematical models or physical problems, refer to the conditions specified at the system boundaries or specific locations to limit the solution. In the context of rock resistivity measurement, boundary conditions play a crucial role in determining the boundary characteristics of the rock, thus enabling a reasonable solution to the resistivity measurement problem.
[0036] S103: Determine the Laplace equations for each point within the target rock that satisfy the boundary conditions.
[0037] The Laplace equation is an important partial differential equation that describes the spatial distribution of a scalar field (typically electric potential, temperature, pressure, etc.). In this invention, the Laplace equation is used to describe voltage distribution. Voltage distribution is often influenced by boundary conditions, which can include the voltage applied to the electrodes, current density, charge density, etc. By introducing appropriate boundary conditions into the Laplace equation, the potential distribution within the rock can be determined, thereby aiding in the interpretation of electrical exploration data and understanding the resistivity distribution of the rock.
[0038] In one possible implementation, the present invention proposes a novel approach to addressing the shortcomings of the Laplace equation, wherein S103 specifically includes sub-steps S1031 and S1032:
[0039] S1031: Determine the relationship between resistivity and voltage at various points within the target rock:
[0040]
[0041] in, Represents the vector differential operator. Indicates electric field strength. Represents electric flux density. Indicates magnetic field strength. Indicates magnetic flux density. q represents current density, and q represents charge density.
[0042] S1032: Construct the Laplace equations that satisfy the boundary conditions at various points within the target rock:
[0043]
[0044] Among them, J n This represents the current density of the excitation current applied to the boundary of the target rock. Let n represent the boundary voltage of the target rock, and n represent the normal vector within the target rock.
[0045] In this invention, the Laplace equation is used as the mathematical model to describe the electric field distribution, ensuring mathematical consistency during the solution and analysis process. This helps ensure that the obtained electric field and resistivity distributions satisfy fundamental electromagnetic physics principles, providing a strong theoretical foundation for data processing.
[0046] S104: Solve for the voltage distribution within the target rock using the finite element method.
[0047] The Finite Element Method (FEM) is a numerical analysis method used to solve partial differential equations and boundary value problems. It can divide a complex physical domain into small elements, establish a mathematical model, and approximate the solution to the physical equations by solving the linear system to obtain the numerical solution to the problem.
[0048] In one possible implementation, the present invention proposes a novel finite element method for solving the problem, wherein S104 specifically includes sub-steps S1041 to S10410:
[0049] S1041: Applying the variational principle to the Laplace equation satisfied at various points within the target rock, the voltage distribution functional is obtained. :
[0050]
[0051] Where (x, y) represents the coordinates of each point, where x represents the x-coordinate and y represents the y-coordinate. This represents the voltage distribution within the target rock, where σ represents resistivity. Represents the vector differential operator. Indicates the spatial extent of the target rock.
[0052] Among them, the voltage distribution functional is used to describe the distribution of electric potential in space.
[0053] S1042: According to the requirements of finite element mesh generation, select the number, density and distribution of meshes, divide the target rock into multiple triangular elements, and mark the elements and nodes to obtain M elements and N nodes.
[0054] It should be noted that the process of dividing the target rock into triangular units is a process of discretization, which divides a continuous physical space into a finite number of discrete regions for numerical analysis. By dividing the target rock into multiple triangular units and applying an interpolation function, this method is applicable to various geological conditions and irregular terrains, making it flexible in handling complex geological situations.
[0055] S1043: Each triangular element has a sub-functionality. The sum of the integrals of the sub-functionalities of all triangular elements constitutes the functionality of the entire field. The functionality can then be expressed as:
[0056]
[0057] in, Let F represent the field of the e-th triangular unit. e (φ) represents the sub-functionality of the e-th triangular unit. m represents the total number of triangular units.
[0058] S1044: The voltage distribution function within each cell is approximated using an interpolation function.
[0059]
[0060] in, Let represent the voltage value at point (x, y), and let a, b, and c represent the coefficients to be solved.
[0061] The three nodes of the triangular element are numbered counterclockwise as I(x1,y1), J(x2,y2), and K(x3,y3). Substituting I(x1,y1), J(x2,y2), and K(x3,y3) into the voltage distribution function, we obtain:
[0062]
[0063] in, Represents the interpolation function matrix. Represents the voltage distribution matrix. This represents the interpolation function term.
[0064] It should be noted that the voltage distribution function within the cell is determined by considering the potential values at the three nodes.
[0065] S1045: Solve for the partial derivatives of the voltage distribution function, and obtain the following:
[0066]
[0067] in, This represents the coefficient matrix of the e-th triangular unit. This represents the voltage distribution matrix of the e-th triangular unit.
[0068] S1046: Calculate the subfunctional F representing the e-th triangular unit based on the partial derivative of the voltage distribution function. e (φ):
[0069]
[0070] in, This indicates the matrix transpose.
[0071] S1047: Calculate the functional of current density in each triangular element:
[0072]
[0073] in, The generalized function representing the current density of the e-th triangular element. Q represents the voltage distribution within the target rock. e This represents the current density of the e-th triangular unit.
[0074] S1048: Subfunctional F for each triangular unit e (φ) and the generalized function of current density By accumulating these values, we obtain the generalized function F(φ) for the entire target rock:
[0075]
[0076] in, This represents the current density matrix of the e-th triangular element. Let M denote the matrix transpose, M denote the total coefficient matrix, and Q denote the total current density matrix.
[0077] It should be noted that summing the integrals of the sub-functionalities of all triangular elements to construct the functionality of the entire field allows for integrated modeling of the entire rock, thus providing a more comprehensive description of the electric field behavior. Simultaneously, it ensures the consistency of voltage distribution and current density throughout the rock, thereby avoiding the accumulation of local errors.
[0078] S1049: Taking the extreme values of the generalized function over the entire target rock, we obtain the finite element equations:
[0079]
[0080] It should be noted that stable numerical solutions can be obtained by taking the extreme values of the functional. Taking the extreme values is equivalent to finding the minimum or maximum points of the functional, which can effectively avoid numerical instability and oscillations. Taking the extreme values of the functional is equivalent to optimizing the voltage distribution within the sensitive field. By finding the extreme points of the functional, the optimal voltage distribution under given boundary conditions can be obtained, thus yielding the optimal solution for the voltage distribution in the sensitive field.
[0081] S10410: Solve the finite element equations using the Gaussian partial principal element elimination method to obtain the voltage distribution within the target rock. .
[0082] It should be noted that the Gaussian elimination method can effectively solve large-scale finite element equations, which means that the voltage distribution within the target rock can be obtained in a relatively short time, thus improving computational efficiency.
[0083] This invention provides a comprehensive, flexible, and highly accurate solution method using a novel finite element method. It innovatively employs the calculation of a generalized function and optimizes the solution by taking the extreme values of the function, ensuring global consistency and numerical stability, and achieving optimal solution. Simultaneously, it decomposes the voltage and resistivity distribution problems into discrete elements, better utilizing the limited, actually measured data to accurately solve for the voltage distribution.
[0084] S105: Based on the voltage distribution within the target rock, the resistivity distribution within the target rock is calculated by inversion.
[0085] Inversion calculation refers to the process of inferring or estimating unknown parameters or scenarios from a set of known observation data or measurement results through mathematical modeling and calculation methods. In electrical exploration, inversion calculation estimates the resistivity distribution within the target rock by using observed voltage distribution data.
[0086] In one possible implementation, the present invention proposes a novel inversion calculation method, wherein S105 specifically includes sub-steps S1051 to S1055:
[0087] S1051: Introducing regularization parameters to construct the objective function for resistivity calculation:
[0088]
[0089] Where minF(ρ) represents the objective function, f d f represents the data fitting error term. m λ represents the regularization constraint term, and λ represents the regularization factor.
[0090] The data fitting error term can be expressed as:
[0091]
[0092] The target rock contains M triangular elements, where ρ represents an M-dimensional resistivity distribution vector, and the elements in ρ correspond to the resistivity values of each element. V(ρ) represents the boundary voltage value when the resistivity distribution is ρ, and U represents the measured boundary voltage value. Represents the modulo operation of a matrix.
[0093] The regularity constraint term can be expressed as:
[0094]
[0095] Where L represents the regularization matrix and ρ0 represents the initial resistivity value.
[0096] In this invention, a regularization constraint term is introduced to impose additional constraints on the objective function, such as smoothness and sparsity, thereby limiting the solution space and reducing the instability of the problem. By controlling the magnitude of the regularization parameter, the degree of constraint can be adjusted, balancing the relationship between data fitting and model constraints, making the solution more stable and reliable, and thus improving ill-conditioned problems.
[0097] S1052: Taking the derivative of the objective function for resistivity calculation and setting the derivative value to 0, we have:
[0098]
[0099] in, The Jacobian matrix can be represented as: .
[0100] S1053: With respect to the derivative Perform a Taylor expansion, retaining only the linear terms of the Taylor expansion, to obtain the increment at the k-th iteration. :
[0101]
[0102] in, This represents the resistivity distribution in the k-th iteration. Let Jacobian matrix be the value of the k-th iteration. Represents the resistivity distribution as The boundary voltage value at time, where U represents the measured boundary voltage value. Let represent the regularization factor, L represent the regularization matrix, and ρ0 represent the initial resistivity value. This indicates the matrix transpose.
[0103] Select I represents the identity matrix. Then the increment It can be simplified to:
[0104]
[0105] S1054: Determine whether the convergence condition is met during the k-th iteration:
[0106]
[0107] Among them, f d This represents the data fitting error term. Represents the resistivity distribution as The boundary voltage value at the time, where U represents the measured boundary voltage value and ε represents the preset error value.
[0108] Those skilled in the art can set the preset error value ε according to the actual situation, and the present invention does not limit it.
[0109] It should be noted that a convergence condition is introduced, which is determined by comparing the error between the current voltage distribution and the measured boundary voltage value to see if the convergence condition is met. This controls the error, ensuring that the iterative process is carried out within the allowable error range and preventing infinite loops and excessive iterations.
[0110] Specifically, by introducing a preset error value ε, the convergence condition can be flexibly adjusted during the iteration process. When the error value reaches the preset range, the convergence condition is met, the iteration can be terminated, and the final resistivity distribution can be output. This allows for adjustment of the preset error value based on specific problems and data conditions, thus providing better control over the iteration process.
[0111] S1055: When the convergence condition is met, terminate the iteration and output the resistivity distribution of the k-th iteration. When the convergence condition is not met, Continue the next iteration until the convergence condition is met.
[0112] It should be noted that an iterative solution can gradually approach the optimal solution. In each iteration, by calculating the derivative and increment of the objective function, the optimization direction of the current resistivity distribution can be found, causing the objective function to gradually decrease. This can accelerate the convergence speed, especially when the objective function is complex and highly nonlinear; iterative optimization can approach the optimal solution more quickly.
[0113] In this invention, by solving the objective function after introducing regularization constraints, additional constraints such as smoothness and sparsity can be introduced into the objective function, thereby limiting the solution space, reducing the instability of the problem, and thus improving ill-conditioned problems. Introducing an iterative solution method, by gradually adjusting parameters and optimizing the solution, can better approximate the optimal solution and improve the success rate and accuracy of the solution, thereby improving the efficiency and stability of solving inverse problems.
[0114] In one possible implementation, after S105, the following is also included:
[0115] Based on the resistivity distribution within the target rock, the resistivity at each point is obtained.
[0116] Calculate the rate of change of resistance at each point:
[0117]
[0118] Where c represents the rate of change of resistance, and ρ represents resistivity. This represents the average resistivity.
[0119] If the rate of change of resistance at the current point is greater than the preset value, it is determined that the resistivity at the current point is abnormal.
[0120] In this invention, by calculating the rate of change of resistivity and comparing it with a preset value, abnormal points in the resistivity distribution can be identified, which helps to eliminate inaccurate or abnormal data caused by noise, equipment problems or other interference factors, thereby improving the quality and reliability of the data.
[0121] When abnormal resistivity is removed, select the four points closest to the current point, take the average resistivity of the four points, and replace the abnormal resistivity with the average resistivity of the four points.
[0122] In this invention, replacing the abnormal resistivity value with the average of the nearest few points helps to smooth the resistivity distribution, reduces local discontinuities, makes the resistivity image more coherent, easier to interpret and understand, helps to reduce the impact of local noise on the data, and helps to improve the stability and consistency of the resistivity distribution.
[0123] S106: Calculate the porosity of the target rock based on the resistivity distribution within the target rock.
[0124] In one possible implementation, the present invention proposes a novel method for calculating porosity, wherein S106 specifically includes sub-steps S1061 and S1062:
[0125] S1061: When the target rock is in a dry region, the porosity within the target rock is calculated using the following formula:
[0126]
[0127] Where ρ represents the resistivity of the rock, ρ a Indicates the resistivity of air. This indicates porosity.
[0128] It should be noted that the resistivity of the air at the geographical location of the target rock is a known quantity that can be found. Specifically, you can refer to the reference values of atmospheric or air resistivity provided in relevant literature and data on geography or environmental science, or you can use a professional resistivity measuring instrument to measure the resistivity of the air.
[0129] S1062: When the target rock is in a humid area, the porosity within the target rock is calculated using the following formula:
[0130]
[0131] Where, ρ w The resistivity of water in pores, S r It indicates the saturation of water in the pores.
[0132] It should be noted that the resistivity and saturation of water in the pores of the target rock are known quantities that can be found. Specifically, you can consult relevant literature and data on geography or environmental science to find reference values for the resistivity and saturation of water in the pores, or you can use professional measuring instruments to measure the resistivity and saturation of water in the pores on-site.
[0133] In this invention, since different moisture conditions affect the resistivity measurement results, different formulas are used in dry and humid areas respectively, which can more accurately estimate the porosity in the target rock.
[0134] S107: Determine the degree of damage to the target rock based on the porosity within the target rock.
[0135] In one possible implementation, the present invention proposes a novel method for determining the degree of rock damage, wherein S107 specifically includes sub-steps S1071 and S1072:
[0136] S1071: Map the porosity within the target rock to a probability value within the interval [0,1] using the sigmoid function.
[0137]
[0138] in, This represents the sigmoid function, where e represents the natural logarithm. This indicates porosity.
[0139] S1072: Converts the damage level of the target rock using a step function and outputs the result.
[0140]
[0141] in, The function represents a step function, where 1 represents intact rock, 2 represents slightly damaged rock, 3 represents moderately damaged rock, and 4 represents severely damaged rock.
[0142] It should be noted that a mathematical relationship has been established between porosity and the degree of rock damage. By mapping probability values to damage levels, the degree of damage can be quantified into specific numerical values. This helps to more accurately describe the physical state of rocks.
[0143] In this invention, the porosity of the target rock is automatically calculated based on its resistivity distribution, and then the degree of damage to the target rock is automatically determined based on this porosity. This eliminates the need for human intervention, reduces the influence of subjective factors, and improves the consistency and accuracy of rock damage assessment.
[0144] In one possible implementation, after S105, the following is also included:
[0145] S108: Construct a resistivity distribution map based on the resistivity distribution at various points within the target rock.
[0146] A resistivity map is a geophysical image used to visualize the resistivity distribution of underground rocks or soil. These images are typically displayed in two or three dimensions to show resistivity values at different locations.
[0147] In this invention, a resistivity distribution map is constructed based on the resistivity distribution at various points within the target rock. This map graphically presents the resistivity distribution within the target rock, making the data more intuitive and easier to understand. Furthermore, the resistivity distribution map can display resistivity values at different locations, thereby helping to identify spatial variations and trends within the rock. This is crucial for understanding the geological characteristics and distribution of rocks.
[0148] In one possible implementation, the present invention proposes a novel method for constructing resistivity distribution maps, wherein S108 specifically includes sub-steps S1081 to S1083:
[0149] S1081: Perform triangular subdivision of each point within the target rock.
[0150] In one possible implementation, sub-step S1081 specifically includes grandchild steps S10811 to S10814:
[0151] S10811: Construct an initial large triangle containing all the points.
[0152] S10812: Select a point to be inserted. If the point to be inserted is inside an existing triangle, connect the point to be inserted to the three vertices of the triangle to form three new triangles. If the point to be inserted is on the common side of two existing triangles, connect the point to be inserted to the remaining vertices of the two triangles to form four new triangles.
[0153] It should be noted that by inserting new points within existing triangles or on common edges, it is possible to ensure that the generated triangles closely match the points within the rock, thereby improving the accuracy of the model.
[0154] In this invention, triangles can be adaptively generated based on the distribution of actual points, making it suitable for point data with different densities and distributions, thereby improving the applicability of the model.
[0155] S10813: When two triangles share a common circumcircle, swap the diagonals of the two triangles.
[0156] It should be noted that when two triangles share a common circumcircle, one of them is an obtuse triangle. When an angle of a triangle is very close to 180 degrees, it may lead to numerical instability or increased error. By exchanging diagonals, this situation can be reduced, improving the robustness of the numerical model. Furthermore, exchanging diagonals of a common circumcircle helps reduce the scalene property of triangles (i.e., large differences in side lengths), thereby improving numerical stability.
[0157] S10814: Select other points in sequence, repeat S1082 and S1083 until all points are inserted, and perform triangulation on each point in the target rock.
[0158] In this invention, by performing fine triangular subdivision of each point within the target rock, high-resolution data representation can be achieved, which helps to capture the details of underground resistivity changes.
[0159] S1082: Generate contour maps based on the resistivity distribution at various points within the target rock.
[0160] In one possible implementation, sub-step S1082 specifically includes grandchild steps S10821 to S10825:
[0161] S10821: Based on the resistivity distribution at various points within the target rock, the coordinates of each point are obtained, where the point coordinates are expressed as... x represents the horizontal axis, y represents the vertical axis, and ρ represents the resistivity.
[0162] S10822: Determine the positional relationship between the contour lines with resistivity ρ* and the target side of the triangle, where the coordinates of the two endpoints of the target side are respectively... and ,but:
[0163] when At that time, the contour line with resistivity ρ* passes through the target edge.
[0164] when At that time, the contour lines with resistivity ρ* are separated from the target edge.
[0165] when When the resistivity is ρ*, the contour line passes through the endpoint of the target side.
[0166] S10823: When a contour line with resistivity ρ* crosses the target edge, determine the coordinates of the intersection point:
[0167]
[0168] Where x* represents the x-coordinate of the intersection point and y* represents the y-coordinate of the intersection point.
[0169] It should be noted that the intersection points determine the position of the contour lines on the graph, which can help depict the boundary between different resistivity values.
[0170] S10824: Connect the points with resistivity ρ* to form contour lines.
[0171] S10825: Smooth the contour lines.
[0172] Specifically, contour lines can be smoothed using methods such as Bézier curve fitting, spline interpolation, and smoothing filtering.
[0173] In this invention, generating contour maps helps to visualize abstract resistivity data in an intuitive way. Contour lines make it easier to understand and interpret the distribution of resistivity within a target rock, thus providing a better understanding of subsurface structures. Furthermore, contour maps show the spatial trends of resistivity variation. Users can identify rapid changes and gradients in resistivity within a region by observing the density, spacing, and shape of the contour lines, which is crucial for studying subsurface geological features and resource exploration.
[0174] S1083: Color fill the contour map to construct a resistivity distribution map.
[0175] In one possible implementation, S1083 specifically includes steps S10831 to S10833:
[0176] S10831: Get the number of contour lines n.
[0177] S10832: n contour lines divide the contour map into n+1 types of regions.
[0178] Reference manual attached Figure 3 This diagram illustrates a contour map provided by the present invention. Specifically, taking three contour lines as an example, the three contour lines are ρ1, ρ2, and ρ3, and... Regions smaller than ρ1 are classified as the first type of region. Figure 3 In the I), the region between ρ1 and ρ2 is considered as the second type of region ( Figure 3 In II), the region between ρ2 and ρ3 is considered as the third type of region ( Figure 3 Regions greater than ρ3 (in the third category) are classified as the fourth category region. Figure 3 (IV in the text). Therefore, the three contour lines divide the contour map into four types of regions.
[0179] S10833: Select n+1 colors and fill the n+1 types of regions respectively to construct a resistivity distribution map.
[0180] Specifically, a color gradient can be used to select a start color and an end color, and then create a smooth color transition between them. This can be achieved by gradually adjusting the channel values in the RGB color space. For example, a transition from blue to red represents a transition from low resistivity to high resistivity.
[0181] It should be noted that each resistivity range has a corresponding fill color, which makes the resistivity distribution map easier to understand and interpret.
[0182] In this invention, a novel resistivity distribution map construction method is used. Through triangulation, contour mapping, and color filling, the constructed resistivity distribution map can be presented in an intuitive and clear manner, making the underground resistivity distribution readily apparent.
[0183] Compared with the prior art, the present invention has at least the following beneficial technical effects:
[0184] (1) In this invention, in the face of limited data obtained from actual field exploration, the Laplace equation is determined by using the actual measured voltage value as the boundary condition. Then, the Laplace equation is solved by the finite element method, and the voltage distribution and resistivity distribution problems are decomposed into discrete elements. The voltage distribution is more accurately solved by using the limited data obtained from actual measurement. Then, the resistivity distribution in the target rock is calculated by inversion based on the voltage distribution in the target rock, thereby improving the accuracy of resistivity distribution calculation.
[0185] (2) In this invention, the porosity of the target rock is automatically calculated based on the resistivity distribution within the target rock, and then the degree of damage to the target rock is automatically determined based on the porosity. This eliminates the need for human intervention, reduces the influence of subjective factors, and improves the consistency and accuracy of rock damage assessment.
[0186] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.
[0187] The above embodiments merely illustrate several implementation methods of the present invention, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these all fall within the protection scope of the present invention. Therefore, the protection scope of this invention patent should be determined by the appended claims.
Claims
1. A method for processing electrical resistivity tomography (OTT) exploration data, characterized in that, include: S101: Distributed electrodes are arranged around the target rock; S102: Obtain the voltage value of each electrode and use the voltage value of each electrode as a boundary condition; S103: Determine the Laplace equations that satisfy the boundary conditions for each point within the target rock. S104: The voltage distribution within the target rock is determined using the finite element method; S105: Based on the voltage distribution within the target rock, the resistivity distribution within the target rock is calculated by inversion; S106: Calculate the porosity of the target rock based on the resistivity distribution within the target rock; S107: Determine the degree of damage to the target rock based on the porosity within the target rock.
2. The electrical resistivity tomography data processing method according to claim 1, characterized in that, S103 specifically includes: S1031: Determine the relationship between resistivity and voltage at various points within the target rock: in, Represents the vector differential operator. Indicates electric field strength. Represents electric flux density. Indicates magnetic field strength. Indicates magnetic flux density. q represents current density, and q represents charge density; S1032: Construct the Laplace equations for each point within the target rock that satisfy the boundary conditions: Among them, J n This represents the current density of the excitation current applied to the boundary of the target rock. Let n represent the boundary voltage of the target rock, and n represent the normal vector within the target rock. This represents the voltage distribution within the target rock, where σ represents the resistivity.
3. The electrical resistivity tomography data processing method according to claim 2, characterized in that, S104 specifically includes: S1041: Applying the variational principle to the Laplace equation satisfied at each point within the target rock, the voltage distribution functional is obtained. : Where (x, y) represents the coordinates of each point, where x represents the x-coordinate and y represents the y-coordinate. This represents the voltage distribution within the target rock, where σ represents resistivity. Represents the vector differential operator. Indicates the spatial extent of the target rock; S1042: According to the requirements of finite element mesh generation, select the number, density and distribution of meshes, divide the target rock into multiple triangular elements, and mark the elements and nodes to obtain M elements and N nodes; S1043: Each triangular element has a sub-functionality. The sum of the integrals of the sub-functionalities of all triangular elements constitutes the functionality of the entire field. The functionality can then be expressed as: in, Let F represent the field of the e-th triangular unit. e (φ) represents the sub-functionality of the e-th triangular unit, and m represents the total number of triangular units; S1044: The voltage distribution function within each cell is approximated using an interpolation function. in, Let (x, y) represent the voltage value at point (x, y), and let a, b, and c represent the coefficients to be solved. The three nodes of the triangular element are numbered counterclockwise as I(x1,y1), J(x2,y2), and K(x3,y3). Substituting I(x1,y1), J(x2,y2), and K(x3,y3) into the voltage distribution function, we obtain: in, Represents the interpolation function matrix. Represents the voltage distribution matrix. Indicates the interpolation function term; S1045: Solve for the partial derivatives of the voltage distribution function, and obtain the following: in, This represents the coefficient matrix of the e-th triangular unit. This represents the voltage distribution matrix of the e-th triangular element; S1046: Calculate the subfunctional F representing the e-th triangular unit based on the partial derivative of the voltage distribution function. e (φ): in, Indicates matrix transpose; S1047: Calculate the functional of current density in each triangular element: in, The generalized function representing the current density of the e-th triangular element. Q represents the voltage distribution within the target rock. e This represents the current density of the e-th triangular unit; S1048: Subfunctional F for each triangular unit e (φ) and the generalized function of current density By accumulating the results, we obtain the generalized function F(φ) for the entire target rock: in, This represents the current density matrix of the e-th triangular element. denoted by matrix transpose, M represents the total coefficient matrix, and Q represents the total current density matrix; S1049: Taking the extreme values of the generalized function of the entire target rock, we obtain the finite element equations: S10410: Solve the finite element equations using the Gaussian partial principal element elimination method to obtain the voltage distribution within the target rock. .
4. The electrical resistivity tomography data processing method according to claim 1, characterized in that, Specifically, S105 includes: S1051: Introducing regularization parameters to construct the objective function for resistivity calculation: Where minF(ρ) represents the objective function, f d f represents the data fitting error term. m λ represents the regularization constraint term, and λ represents the regularization factor. The data fitting error term can be expressed as: The target rock contains M triangular units, where ρ represents an M-dimensional resistivity distribution vector, the elements in ρ correspond to the resistivity values of each unit, V(ρ) represents the boundary voltage value when the resistivity distribution is ρ, and U represents the measured boundary voltage value. Represents the modulo operation of a matrix; The regularity constraint term can be expressed as: Where L represents the regularization matrix and ρ0 represents the initial resistivity value; S1052: Taking the derivative of the objective function for resistivity calculation and setting the derivative value to 0, we have: in, The Jacobian matrix can be represented as: ; S1053: With respect to the derivative Perform a Taylor expansion, retaining only the linear terms of the Taylor expansion, to obtain the increment at the k-th iteration. : in, This represents the resistivity distribution in the k-th iteration. Let Jacobian matrix be the value of the k-th iteration. Represents the resistivity distribution as The boundary voltage value at time, where U represents the measured boundary voltage value. Let represent the regularization factor, L represent the regularization matrix, and ρ0 represent the initial resistivity value. Indicates matrix transpose; Select I represents the identity matrix. Then the increment It can be simplified to: S1054: Determine whether the convergence condition is met during the k-th iteration: Among them, f d This represents the data fitting error term. Represents the resistivity distribution as The boundary voltage value at the time, where U represents the measured boundary voltage value and ε represents the preset error value; S1055: When the convergence condition is met, the iteration ends, and the resistivity distribution of the k-th iteration is output. When the convergence condition is not met, Continue the next iteration until the convergence condition is met.
5. The electrical resistivity tomography data processing method according to claim 1, characterized in that, S106 specifically includes: S1061: When the target rock is in a dry region, the porosity within the target rock is calculated using the following formula: Where ρ represents the resistivity of the rock, ρ a Indicates the resistivity of air. Indicates porosity; S1062: When the target rock is in a humid area, the porosity of the target rock is calculated using the following formula: Where, ρ w The resistivity of water in pores, S r It indicates the saturation of water in the pores.
6. The electrical resistivity tomography data processing method according to claim 1, characterized in that, S107 specifically includes: S1071: The porosity within the target rock is mapped to a probability value within the interval [0,1] using the sigmoid function. in, This represents the sigmoid function, where e represents the natural logarithm. Indicates porosity; S1072: Convert the damage level of the target rock using a step function and output: in, The function represents a step function, where 1 represents intact rock, 2 represents slightly damaged rock, 3 represents moderately damaged rock, and 4 represents severely damaged rock.
7. The electrical resistivity tomography data processing method according to claim 1, characterized in that, Following S105, the following is also included: S108: Construct a resistivity distribution map based on the resistivity distribution at various points within the target rock.
8. The electrical resistivity tomography data processing method according to claim 7, characterized in that, S108 specifically includes: S1081: Perform triangular subdivision of each point within the target rock; S1082: Generate contour maps based on the resistivity distribution at various points within the target rock; S1083: Fill the contour map with color to construct the resistivity distribution map.
9. The electrical resistivity tomography data processing method according to claim 8, characterized in that, S1081 specifically includes: S10811: Construct an initial large triangle containing all the points; S10812: Select a point to be inserted. When the point to be inserted is inside an existing triangle, connect the point to be inserted to the three vertices of the triangle to form three new triangles. When the point to be inserted is on the common side of two existing triangles, connect the point to be inserted to the remaining vertices of the two triangles to form four new triangles. S10813: When two triangles share a common circumcircle, swap the diagonals of the two triangles. S10814: Select other points in sequence, repeat S1082 and S1083 until all points are inserted, and perform triangulation on each point in the target rock.
10. The electrical resistivity tomography data processing method according to claim 9, characterized in that, S1082 specifically includes: S10821: Based on the resistivity distribution at various points within the target rock, the coordinates of each point are obtained, where the point coordinates are expressed as... x represents the horizontal axis, y represents the vertical axis, and ρ represents the resistivity; S10822: Determine the positional relationship between the contour lines with resistivity ρ* and the target side of the triangle, where the coordinates of the two endpoints of the target side are respectively... and ,but: when At that time, contour lines with resistivity ρ* pass through the target edge; when At that time, the contour lines with resistivity ρ* are separated from the target edge; when At that time, contour lines with resistivity ρ* pass through the endpoints of the target edge; S10823: When a contour line with resistivity ρ* crosses the target edge, determine the coordinates of the intersection point: Where x* represents the x-coordinate of the intersection point and y* represents the y-coordinate of the intersection point; S10824: Connect the points with resistivity ρ* to form contour lines; S10825: Smooth the contour lines.