Cement mixing pile hole forming uniformity detection method based on resistivity three-dimensional imaging

Abstract

The application discloses a cement mixing pile hole forming uniformity detection method based on resistivity three-dimensional imaging, and particularly relates to the field of nondestructive testing of engineering materials, and is used for solving the problem that the general resistivity inversion method is blurred in detecting the real defect boundary due to the smoothness constraint when detecting the discrete non-uniform body in the cement mixing pile, and false abnormal transition zones are prone to be generated, thereby affecting the detection accuracy; effective data is obtained by acquiring resistivity measurement data and screening, a partition constraint inversion strategy is adopted based on the effective data to weaken the smooth constraint on the core area of the pile body and to retain the conventional smooth constraint on the surrounding medium area, so as to generate a three-dimensional resistivity distribution result, the density clustering and spectral density analysis are extracted and utilized to verify the time and space of the electrical property mutation boundary, the continuous false abnormal transition zone is removed through similarity matching to generate a purified resistivity distribution result, and the uniformity detection conclusion is output based on the result, and the accuracy of the internal non-uniform defect recognition of the cement mixing pile is improved.

Description

Technical Field

[0001] This invention relates to the field of nondestructive testing technology for engineering materials, and more specifically, to a method for detecting the uniformity of borehole formation in cement mixing piles based on three-dimensional resistivity imaging. Background Technology

[0002] In the field of geotechnical engineering quality testing, resistivity three-dimensional imaging technology can be used to non-destructively evaluate the internal uniformity of artificial composite piles such as cement mixing piles. The core of this method is to measure the resistivity distribution of the pile body and the surrounding medium, and reconstruct the three-dimensional electrical structure model inside the pile body using a three-dimensional inversion algorithm. Then, the mixing quality of cement and soil can be judged based on the uniformity of the model. In the existing technology, the realization of resistivity three-dimensional imaging usually relies on a general inversion algorithm represented by the smooth constraint least squares method. By constructing and solving an optimization problem with the common goal of minimizing the data fitting difference and model roughness, a three-dimensional resistivity distribution image with a gradual spatial variation can be obtained. It is widely used when dealing with geological bodies with macroscopic electrical stratification or gradual changes.

[0003] However, when the aforementioned general three-dimensional resistivity inversion method is directly applied to the detection of borehole uniformity in cement mixing piles, the cement mixing pile body is a composite material formed by mechanically mixing multiple phase media such as cement, soil, and water. Its interior may contain discrete non-uniform bodies such as cement-rich lumps, unevenly mixed soil interlayers, or pores. The spatial smoothness constraint pursued by the general inversion algorithm is in fundamental conflict with the non-uniform characteristics of the cement-soil medium that may actually exist, which are spatially discrete and have significant electrical differences. This causes the three-dimensional resistivity model obtained by inversion to tend to blur or smooth out these real, clearly defined defects. At the same time, it may generate non-existent, continuous false anomaly transition zones in order to fit the measurement data, making it difficult to accurately and reliably distinguish between real non-uniform defects and imaging artifacts introduced by the algorithm, thus affecting the accuracy and credibility of the detection conclusions. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a method for detecting the uniformity of cement mixing pile hole formation based on resistivity three-dimensional imaging to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A method for detecting the uniformity of borehole formation in cement mixing piles based on three-dimensional resistivity imaging includes:

[0007] S1. Obtain resistivity measurement data of cement mixing piles and surrounding medium;

[0008] S2. The resistivity measurement data is filtered using an anomaly detection algorithm to obtain valid resistivity data;

[0009] S3. Based on the effective resistivity data, a partition constraint inversion strategy is adopted to divide the core area of ​​the pile into the surrounding medium area. The smooth constraint is weakened in the core area of ​​the pile, while the conventional smooth constraint is retained in the surrounding medium area to generate a three-dimensional resistivity distribution result.

[0010] S4. Extract the electrical abrupt boundary from the three-dimensional resistivity distribution results, use density clustering algorithm to analyze the spatial distribution clustering confidence of the electrical abrupt boundary, and screen out the electrical abrupt boundary with a confidence level higher than the preset confidence threshold. Then, use spectral density estimation algorithm to perform time-series spectral density analysis on the resistivity fluctuation of the screened electrical abrupt boundary to generate a set of spatiotemporally verified electrical abrupt boundaries.

[0011] S5. Based on the set of electrical abrupt change boundaries verified by time and space, a similarity matching algorithm is used to remove continuous false abnormal transition regions and generate the purified resistivity distribution results.

[0012] S6. Based on the resistivity distribution results after purification, output the test conclusion on the uniformity of the cement mixing pile hole formation.

[0013] Furthermore, S1 includes:

[0014] A multi-ring electrode array is deployed around the cement mixing pile;

[0015] A high-density electrical resistivity measurement mode combining multiple electrode spacing and multiple orientations was adopted to collect apparent resistivity data of cement mixing piles and surrounding media under different combinations of power supply electrodes and measuring electrodes, which were used as resistivity measurement data.

[0016] Furthermore, S2 includes:

[0017] The overall distribution characteristics of resistivity measurement data were analyzed based on statistical methods.

[0018] Determine the threshold for identifying abnormal data based on the overall distribution characteristics;

[0019] Measurements that deviate from the overall distribution characteristics by more than a threshold in the resistivity measurement data are considered outliers and are removed. The remaining data constitute the valid resistivity data.

[0020] Furthermore, S3 includes:

[0021] The spatial range of the core area and the surrounding medium area of ​​the cement mixing pile is defined based on the design location and geometric dimensions of the cement mixing pile;

[0022] In the inversion calculation, the part belonging to the core area of ​​the pile is subject to a minimization constraint based on the first norm or an anisotropic smoothing constraint with reduced lateral smoothing weight, while the part belonging to the surrounding medium area is subject to a minimization constraint based on the second norm.

[0023] A three-dimensional inversion calculation is performed using the effective resistivity data as the fitting target to generate a three-dimensional resistivity distribution result.

[0024] Furthermore, a three-dimensional inversion calculation is performed using the effective resistivity data as the fitting target to generate a three-dimensional resistivity distribution result. Specifically, an inversion objective function containing data fitting terms and partition smoothing constraint terms is constructed. The resistivity parameters in the three-dimensional space are adjusted through an iterative algorithm to minimize the fitting difference between the forward simulation data and the effective resistivity data until the convergence condition is met, and the final three-dimensional resistivity distribution result is output.

[0025] Furthermore, S4 includes:

[0026] Edge detection is performed on the three-dimensional resistivity distribution results to identify spatial locations where the spatial rate of resistivity change exceeds a preset gradient threshold, thus forming an initial set of electrical abrupt change boundary points.

[0027] A density-based spatial clustering algorithm is applied to the initial set of boundary points of electrical abrupt change. Different spatial clusters are divided according to the density distribution of spatial location points, and the clustering confidence of each spatial cluster is calculated.

[0028] Spatial clusters with clustering confidence scores higher than a preset confidence threshold are identified as candidate electrical mutation boundaries;

[0029] For each candidate electrical abrupt change boundary, continuous resistivity values ​​are extracted along its spatial extension direction to form a resistivity fluctuation sequence, and the power spectral density of the resistivity fluctuation sequence is estimated.

[0030] Candidate electrical abrupt change boundaries that exhibit significant spectral peak characteristics within the frequency range determined by resistivity fluctuation sequence feature analysis are included in the spatiotemporally validated set of electrical abrupt change boundaries.

[0031] Furthermore, the frequency range determined based on the resistivity fluctuation sequence characteristic analysis is specifically determined in the following way: calculate the power spectral density of the resistivity fluctuation sequence, identify the peak values ​​in the power spectrum that are significantly higher than the average background level, and determine the frequency range covered by the peak values.

[0032] Furthermore, S5 includes:

[0033] In a set of spatiotemporally validated electrical abrupt change boundaries, identify spatially adjacent pairs of electrical abrupt change boundaries;

[0034] Calculate the similarity of resistivity distribution in the regions defined by each pair of adjacent electrical abrupt boundary lines;

[0035] Regions with resistivity distribution similarity higher than a preset similarity threshold are identified as continuous pseudo-abnormal transition regions;

[0036] In the three-dimensional resistivity distribution results, the resistivity values ​​identified as continuous pseudo-abnormal transition zones are smoothed and corrected to generate the purified resistivity distribution results.

[0037] Furthermore, the resistivity distribution similarity of the regions defined by each pair of adjacent electrical abrupt boundary is calculated. Specifically, the resistivity spatial distribution data of the regions defined by each pair of adjacent electrical abrupt boundary is extracted, and the resistivity spatial distribution data is correlated or structurally similar between different directions or profiles to quantify the resistivity distribution similarity of the regions defined by each pair of adjacent electrical abrupt boundary.

[0038] Furthermore, S6 includes:

[0039] Extract the resistivity value of the core area of ​​the pile from the resistivity distribution results after purification;

[0040] Calculate the statistical coefficient of variation of resistivity values ​​in the core area of ​​the pile, and extract the spatial location and resistivity contrast of the anomalies corresponding to the set of electrical abrupt boundary verified in time and space.

[0041] The statistical coefficient of variation and the spatial location of the anomaly are compared with the resistivity contrast to the preset uniformity evaluation criteria.

[0042] Based on the comparison results, a conclusion is generated regarding the test results of the uniformity of the drilling of cement mixing piles.

[0043] Compared with the prior art, the present invention has the following beneficial effects:

[0044] 1. Significantly improves the accuracy and reliability of non-destructive testing of borehole uniformity in cement mixing piles based on resistivity 3D imaging technology. By introducing a partitioned constraint inversion strategy, differentiated smoothness constraints are applied to the core area of ​​the pile and the surrounding medium area during the inversion calculation process. The core lies in weakening or even allowing discontinuous changes in the smoothness constraint of the core area of ​​the pile. From the physical constraint level of the inversion algorithm, it changes the mathematical representation preference of the electrical structure of the pile area, so that the 3D resistivity distribution results generated by the inversion can more realistically reflect the discrete inhomogeneous bodies with significant electrical differences and relatively clear boundaries that may exist inside the cement mixing pile, such as cement-rich lumps or unevenly mixed interlayers. This effectively overcomes the inherent limitation of the general smooth constraint inversion method, which inevitably blurs or smooths out the real defect boundaries when applied to such artificial composite materials, and lays a real imaging foundation for accurate defect identification in the future.

[0045] 2. By designing a multi-step spatiotemporal verification and data purification process, the inversion imaging results are deeply analyzed and purified. Density clustering algorithm is used to screen out the true set of anomaly boundaries from the perspective of spatial distribution reliability. Then, spectral density analysis is used to verify these boundaries from the perspective of resistivity fluctuation time series characteristics. Finally, similarity matching is used to identify and remove pseudo-continuous transition regions that may be generated by the inversion algorithm to fit the data and connect real anomalies. This can actively distinguish and remove imaging artifacts introduced by the algorithm itself. As a result, the purified resistivity distribution results used to evaluate the uniformity of the pile body, as well as the defect feature parameters extracted from them, such as statistical coefficient of variation and anomaly contrast, are more representative of the true physical state of the pile material. This greatly reduces the risk of misjudgment caused by defects in the imaging method, thus making the output detection conclusions more reliable. Attached Figure Description

[0046] Figure 1 This is a flowchart of the method for detecting the uniformity of borehole formation in cement mixing piles based on three-dimensional resistivity imaging according to the present invention.

[0047] Figure 2 This is a flowchart illustrating the generation of the spatiotemporally verified electrical mutation boundary set of this invention. Detailed Implementation

[0048] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0049] Example: Figure 1 The present invention provides a method for detecting the uniformity of borehole formation in cement mixing piles based on three-dimensional resistivity imaging, comprising:

[0050] S1. Obtain resistivity measurement data of cement mixing piles and surrounding medium;

[0051] S2. The resistivity measurement data is filtered using an anomaly detection algorithm to obtain valid resistivity data;

[0052] S3. Based on the effective resistivity data, a partition constraint inversion strategy is adopted to divide the core area of ​​the pile into the surrounding medium area. The smooth constraint is weakened in the core area of ​​the pile, while the conventional smooth constraint is retained in the surrounding medium area to generate a three-dimensional resistivity distribution result.

[0053] S4. Extract the electrical abrupt boundary from the three-dimensional resistivity distribution results, use density clustering algorithm to analyze the spatial distribution clustering confidence of the electrical abrupt boundary, and screen out the electrical abrupt boundary with a confidence level higher than the preset confidence threshold. Then, use spectral density estimation algorithm to perform time-series spectral density analysis on the resistivity fluctuation of the screened electrical abrupt boundary to generate a set of spatiotemporally verified electrical abrupt boundaries.

[0054] S5. Based on the set of electrical abrupt change boundaries verified by time and space, a similarity matching algorithm is used to remove continuous false abnormal transition regions and generate the purified resistivity distribution results.

[0055] S6. Based on the resistivity distribution results after purification, output the test conclusion on the uniformity of the cement mixing pile hole formation.

[0056] S1. Obtain resistivity measurement data of cement mixing piles and surrounding medium, specifically implemented as follows:

[0057] An electrode array for resistivity measurement is deployed on the ground based on the center location of the cement mixing pile. The electrode array is arranged with the pile center as the geometric center, forming multiple concentric rings. For example, three concentric rings can be arranged, with the radii of the inner and outer rings set according to the design radius of the cement mixing pile, 1.5 times the pile diameter, and 2 times the pile diameter, respectively. Multiple measuring electrodes are distributed at equal angular intervals on each ring; for example, eight measuring electrodes in the inner ring, twelve in the middle ring, and sixteen in the outer ring. All measuring electrodes are made of metal and are driven into the soil by hammering or placing them into pre-drilled holes to ensure a stable electrical connection between the electrodes and the underground medium. The electrode deployment area needs to completely cover the projected area of ​​the cement mixing pile and the surrounding soil area of ​​interest.

[0058] After the electrodes are laid out, a high-density resistivity measuring instrument is connected for data acquisition. This instrument includes a programmable electrode switching device, a constant current source device, a high-precision voltage measuring device, and a control and recording device. The measurement employs a scanning mode combining multiple electrode spacing and multiple azimuths. Multiple electrode spacing is achieved by systematically changing the distance between the power supply electrode pair and the measuring electrode pair in the measurement sequence. For example, when using a Wenner arrangement, the electrode spacing is gradually increased starting from the minimum spacing to obtain responses at different detection depths. Multiple azimuths are achieved by changing the relative azimuth angles of the power supply electrode pair and the measuring electrode pair around the pile center multiple times at the same electrode spacing, for example, measuring every 45 degrees, to obtain information on changes in the horizontal direction.

[0059] The specific data acquisition operation is automatically executed by the control and recording device according to a pre-programmed measurement sequence. The measurement sequence control electrode switching device sequentially connects the constant current source device and the high-precision voltage measuring device to different electrode combinations. For each connection, the constant current source device injects a stable DC current of a known magnitude between the two selected power supply electrodes; for example, the current magnitude can be selected between 1 mA and 100 mA, with the specific value determined after background noise testing on-site to ensure that the generated signal voltage is much greater than the background noise. Simultaneously, the high-precision voltage measuring device measures the potential difference between the other two selected measurement electrodes. Subsequently, based on the magnitude of the injected current, the measured potential difference value, and the geometric arrangement coefficient of the quadrupole electrode arrangement, the apparent resistivity value of this measurement is calculated.

[0060] The measurement sequence needs to cover all valuable quadrupole combinations. One feasible sequence design method is to first fix one power supply electrode pair, then traverse all measurement electrode pairs that do not coincide with the power supply electrode to measure the potential difference, and then switch to the next power supply electrode pair, repeating the traversal measurement process until all possible valid combinations have been measured. Through this system, apparent resistivity data of cement mixing piles and the surrounding medium under a large number of different power supply and measurement electrode combinations can be collected. All these raw apparent resistivity data, corresponding electrode numbers, and current values ​​and geometric arrangement coefficients used in the calculation are stored by a control recording device to form a complete resistivity measurement dataset. The geometric parameters of the electrode layout, such as the accurate radius of each ring, the azimuth angle and planar coordinates of each measurement electrode, are accurately measured and recorded as the topological basis for building the physical model in subsequent 3D inversion. The current intensity selected during measurement, the duration of each current injection, and the voltage integration time window parameters all need to be determined through preliminary testing on site. The preliminary testing method involves selecting a few electrode arrangements for trial measurements before formal measurements. The ratio of the measured potential difference signal to background noise is observed, and the current magnitude and time parameters are adjusted until the signal stabilizes and meets preset requirements, such as a signal-to-noise ratio greater than 10:1. By deploying a multi-ring electrode array around the pile and employing a systematic multi-pole, multi-directional electrical resistivity scanning mode, a complete set of raw apparent resistivity data reflecting the three-dimensional electrical distribution of the pile and surrounding medium can be obtained. This set serves as the resistivity measurement data for subsequent processing.

[0061] S2. An anomaly detection algorithm is used to filter the resistivity measurement data to obtain valid resistivity data. The specific implementation is as follows:

[0062] After the resistivity measurement data is collected, all stored apparent resistivity data are first imported into the data processing stage for quality screening. The goal is to eliminate unreliable outlier data points and form valid resistivity data for subsequent inversion calculations. This screening process is based on statistical principles to quantitatively analyze the overall distribution characteristics of the dataset.

[0063] This analysis uses statistical methods to examine the overall distribution characteristics of resistivity measurement data. Specifically, it involves calculating several statistical characteristics of the entire resistivity measurement dataset. These include the arithmetic mean and standard deviation of all apparent resistivity data. The arithmetic mean represents the central tendency of the dataset, i.e., the average estimate of the overall resistivity level. It is calculated by summing all apparent resistivity values ​​in the dataset and then dividing the sum by the total number of data points. The standard deviation quantifies the dispersion of the data around the mean, reflecting the overall amplitude of the measurement fluctuations. It is calculated by first calculating the difference between each apparent resistivity value and the mean, then squaring each difference, summing all the squares, dividing this sum by the total number of data points, and finally taking the square root of the quotient to obtain the standard deviation. Through these calculations, two key parameters describing the overall distribution of this batch of resistivity measurement data—the mean and the standard deviation—can be obtained.

[0064] The threshold for identifying outliers is determined based on the overall distribution characteristics. The core of this step is to use the mean and standard deviation calculated in the previous step to define a reasonable range for normal data fluctuations. Measurements exceeding this range are considered potential outliers. The threshold is specifically set by using the calculated mean as the central baseline, and then expanding upwards and downwards by several times the standard deviation, thus forming a pair of numerical limits: an upper and lower threshold. For example, the threshold for identifying outliers can be set as the mean plus three times the standard deviation as the upper threshold, and the mean minus three times the standard deviation as the lower threshold. The standard deviation multiplier is a parameter that can be adjusted based on data quality. The selection of this multiplier typically considers the stringency of data screening requirements and the actual data distribution. For example, in the initial screening stage, two times the standard deviation can be used to conservatively remove obvious outliers, while in the subsequent refined processing stage pursuing higher data quality, three times the standard deviation can be used to achieve a balance between data utilization and data purity. The initial value of this standard deviation factor can be preset based on experience before data processing begins, for example, a default setting of 3. After the initial threshold application and data removal are completed, observe the proportion of removed data to the total original data. If this proportion is within an acceptable range, such as between 1% and 10%, maintain this factor setting. If the proportion is too high or too low, the factor can be fine-tuned and the threshold recalculated. The upper and lower thresholds finally determined through the above process together constitute a closed numerical range, which is statistically considered to be the normal range of the current resistivity measurement data.

[0065] Perform outlier identification and removal. Each original apparent resistivity value in the resistivity measurement data is compared one by one with the upper and lower threshold values ​​determined in the previous step. The specific comparison rule is to check whether the currently evaluated apparent resistivity value is greater than the upper threshold or less than the lower threshold. If the value meets either condition, it is determined that the apparent resistivity value deviates from the overall distribution characteristics and exceeds the threshold used to identify outliers, thus marking it as outlier data. All measurements marked as outliers are removed from the original resistivity measurement dataset. For example, suppose a specific apparent resistivity value is 50 ohm-meters, while the upper threshold calculated in the previous step is 45 ohm-meters and the lower threshold is 5 ohm-meters. Since 50 is greater than 45, exceeding the upper threshold, this value is determined as outlier data and removed.

[0066] After completing the traversal, comparison, and elimination of all data points, the remaining apparent resistivity values ​​in the original resistivity measurement dataset that were not marked as outliers were retained and reorganized into a new dataset. This new dataset eliminated obvious outliers caused by poor electrode contact, strong transient environmental interference, instrument recording errors, etc., resulting in more reliable data quality. This newly constructed dataset is defined as the effective resistivity data, which will serve as the direct input data for subsequent 3D inversion calculations. The entire screening process, through specific statistical calculations, threshold settings, and comparison logic, ensures that the data input to the inversion algorithm has consistent statistical characteristics, reducing the negative impact of outliers on the stability of the inversion results, thus laying the foundation for generating accurate 3D resistivity distribution results that reflect geological conditions. In practice, if the rejection ratio calculated based on the amount of data to be rejected after the initial screening is considered abnormally high, such as exceeding 20%, it may indicate that the threshold setting is too strict or that there are systemic problems with the quality of the original data. In this case, you can return to the step of adjusting the standard deviation factor, for example, by adjusting the factor from 3 to 2.5, recalculate the threshold, and perform the screening process again until the amount of effective resistivity data retained finally meets the basic requirements for the amount of input data for subsequent three-dimensional inversion calculations.

[0067] S3. Based on the effective resistivity data, a partitioned constraint inversion strategy is used to divide the pile core area and the surrounding medium area. The smoothing constraint is weakened in the pile core area, while the conventional smoothing constraint is retained in the surrounding medium area, generating a three-dimensional resistivity distribution result. The specific implementation is as follows:

[0068] After obtaining effective resistivity data, a three-dimensional resistivity distribution result is generated based on this data and using a partitioned constraint inversion strategy. The spatial range of the core area and the surrounding medium area of ​​the cement mixing pile is defined according to the design location and geometric dimensions of the cement mixing pile. In specific implementation, a three-dimensional spatial computational grid corresponding to the aforementioned electrode layout area is established. This grid geometrically discretizes the underground space to be inverted into multiple closely arranged regular three-dimensional cubic units. The design location of the cement mixing pile is determined by two specific parameters: the plane coordinates of its center point and the design elevation of the pile top, as clearly given on its design drawings. Its geometric dimensions mainly refer to the design pile diameter and design pile length. The spatial range of the core area of ​​the pile is precisely defined in the three-dimensional computational grid as a cylindrical region with the design pile center as the vertical axis and the design pile diameter as the cylinder diameter. This cylinder extends vertically downwards from the design pile top elevation to the design pile bottom elevation. The spatial range of the surrounding medium area is defined as all the remaining portion of the entire spatial volume covered by the three-dimensional computational grid after deducting the cylindrical volume occupied by the core area of ​​the pile. In a computer program, by calculating and determining the geometric positional relationship between the three-dimensional coordinates of the center point of each computational grid cell and the cylindrical space definition of the core area of ​​the pile, all grid cells can be clearly divided into a set of cells belonging to the core area of ​​the pile or a set of cells belonging to the surrounding medium area, thereby completing the spatial partitioning of the inversion model.

[0069] In the inversion calculation, different mathematical constraint strategies are applied to the cell sets of different partitions. For the part belonging to the core area of ​​the pile, a minimum constraint based on the first norm or an anisotropic smoothing constraint with reduced lateral smoothing weight is applied. The specific implementation of the minimum constraint based on the first norm is as follows: when constructing the mathematical expression of the inversion objective function, for the resistivity parameters of all grid cells in the core area of ​​the pile, the change of these resistivity parameters relative to a certain initial reference model or relative to the corresponding resistivity value of the model in the previous iteration step is calculated one by one. Then, the absolute value of the change of each cell is taken, and the absolute values ​​of all grid cells are summed to obtain a sum. Finally, this sum is multiplied by a constraint weight coefficient specifically for the core area of ​​the pile and added as a constraint term to the overall objective function. The specific value of the constraint weight coefficient for the core area of ​​the pile is used to control the relative strength of the first norm constraint term in the entire objective function. The determination of its value can be accomplished through systematic numerical experiments. For example, start with a small empirical value as the initial value, and then gradually increase the value. After each adjustment, run a simplified inversion calculation. Observe whether the resistivity distribution image of the core area of ​​the pile in the output three-dimensional resistivity model shows that obvious abrupt changes and sharp boundary features are allowed, while maintaining overall stability. This will determine whether the weight coefficient is appropriate. When the model features meet the requirements, the adjustment is stopped and the weight coefficient value is fixed. The anisotropic smoothing constraint with reduced lateral smoothing weight is implemented by independently calculating the square of the difference between the resistivity values ​​of each grid cell in the core area of ​​the pile and its adjacent grid cells in two mutually perpendicular directions in the horizontal plane and in the vertical direction orthogonal to the horizontal plane. However, when calculating the contribution of the square difference in the horizontal direction, a pre-set lateral smoothing weight factor less than 1 is multiplied, while when calculating the contribution of the square difference in the vertical direction, a weight factor with a default value of 1 is used. Then, the weighted square differences of all cells in all directions are summed and multiplied by the constraint weight coefficient of the core area of ​​the pile to form a constraint term that is added to the objective function. The constraint term constructed in this way mathematically allows for more drastic resistivity changes in the horizontal direction, thereby physically adapting to the possible vertical stratification or local blocky anomalies in the electrical structure of the cement mixing pile.

[0070] For the portion belonging to the surrounding medium region, a minimization constraint based on the L2 norm is applied. This is achieved by calculating the change in resistivity parameters of all grid cells within the surrounding medium region relative to the resistivity values ​​of a predefined initial uniform model or relative to the resistivity values ​​of adjacent model cells when constructing the mathematical expression for the inversion objective function. Then, the change in resistivity of each cell is squared, and the squared values ​​of all grid cells are summed to obtain a total. Finally, this sum is multiplied by a constraint weighting coefficient specifically for the surrounding medium region and added as a constraint term to the overall objective function. The setting of the constraint weight coefficients for the surrounding medium region is usually aimed at achieving a smooth spatial transition of the model. Their values ​​can be determined by setting an initial reference value based on general geophysical inversion experience, and then manually fine-tuned by observing the balance between the data fit difference and the overall model roughness during the inversion iteration process. Alternatively, the L-curve method can be used to automatically determine the values. The L-curve method involves calculating and plotting a series of curves showing the relationship between the data fit difference norm and the model roughness norm on logarithmic coordinates for different weight coefficient values. This curve typically presents an L-shape. The weight coefficient value corresponding to the inflection point of the L-shape is selected as the final value, thus achieving a reasonable balance between data fit and model smoothness. The smoothing constraint of the surrounding medium region is usually applied isotropically, meaning that the same weight factor with a value of 1 is used in calculations in all spatial directions, including both horizontal and vertical directions.

[0071] A three-dimensional inversion calculation is performed using effective resistivity data as the fitting target to generate a three-dimensional resistivity distribution result. Specifically, this involves constructing an inversion objective function that includes a data fitting term and a partitioned smoothing constraint term. The data fitting term measures the overall difference between the forward simulation data and the effective resistivity data. It is typically calculated using the squared form of the L2 norm of the difference between corresponding data points. Specifically, it calculates the difference between the apparent resistivity value of the forward simulation corresponding to each electrode arrangement and the apparent resistivity value of the corresponding measurement point in the effective resistivity data. Then, it squares all the differences and sums all the squared values ​​to obtain a total. The forward simulation data is a sequence of theoretical apparent resistivity values ​​calculated by solving the steady current field distribution of a point current source in a three-dimensional non-uniform conductive medium, based on the resistivity value of each grid cell in the current iteration of the three-dimensional resistivity model and the known electrode arrangement geometric parameters. The partitioned smoothing constraint term is a linear sum of the aforementioned core area constraint term and the surrounding medium area constraint term, weighted according to their respective determined constraint weight coefficients. The final inversion objective function used for iterative solution is a linear combination of the data fitting term and the partition smoothing constraint term multiplied by a global regularization coefficient. This global regularization coefficient is used to globally control the relative importance between the data fitting term and the model constraint term. Its value can be determined by trying several values ​​of different orders of magnitude and running the inversion calculation separately. Then, based on the complexity of the output model and the size of the data fitting residuals, a value that can produce reasonable geological interpretation results and whose fitting residuals are acceptable is selected as the final value.

[0072] The resistivity parameter values ​​of all grid cells in 3D space are iteratively adjusted using an algorithm to minimize the fitting difference between the forward modeling data and the effective resistivity data. The iterative algorithm starts with a pre-defined initial 3D resistivity model, for example, setting the resistivity values ​​of all grid cells to a uniform background resistivity value estimated based on site experience or previous surveys. In each iteration, forward modeling data is first obtained using the current 3D resistivity model and known electrode arrangement parameters. Then, the objective function value under the current model is calculated using the method described above. Next, the gradient vector of the objective function with respect to the resistivity parameter of each grid cell is calculated. This gradient vector indicates the direction of the steepest decrease in each resistivity parameter to reduce the objective function value. The optimal step size for updating the model parameters along the gradient direction is then determined using a linear search method or by approximating the Hessian matrix constructed using a quasi-Newton method, thereby calculating the specific correction amount for the resistivity parameter of each grid cell. The current resistivity values ​​of all grid cells are then added to the corresponding correction amount to obtain the updated 3D resistivity model. This calculation process is repeated continuously, forming an iterative loop.

[0073] The iterative process continues until a preset convergence condition is met. The convergence condition is typically set to determine if the number of iterations has reached a preset maximum limit, such as 50 iterations. When any convergence condition is met, the iteration loop terminates, and the resistivity values ​​of all grid cells in the currently obtained 3D resistivity model are output as the final result. This final result is a 3D data volume containing the specific resistivity values ​​of each grid cell, i.e., the 3D resistivity distribution result. This result reflects a 3D electrical structure that, under the partitioning constraint strategy, fits effective resistivity data while allowing for electrical abrupt changes in the core area of ​​the pile and requiring smooth changes in the surrounding medium. The entire inversion calculation process ensures the physical rationality of the inversion imaging for the cement mixing pile detection scenario.

[0074] Figure 2 The flowchart of generating the spatiotemporally verified set of electrical abrupt boundary points according to the present invention is given. S4: Extract the electrical abrupt boundary points from the three-dimensional resistivity distribution results, analyze the spatial distribution clustering confidence of the electrical abrupt boundary points using a density clustering algorithm, and screen out the electrical abrupt boundary points with confidence scores higher than a preset confidence threshold. Then, perform time-series spectral density analysis on the resistivity fluctuations of the screened electrical abrupt boundary points using a spectral density estimation algorithm to generate the spatiotemporally verified set of electrical abrupt boundary points. The specific implementation is as follows:

[0075] After obtaining the three-dimensional resistivity distribution results, electrical abrupt change boundaries that may reflect the inhomogeneity within the cement mixing pile are extracted, and these boundaries are spatiotemporally verified to form a reliable boundary set. This process first performs edge detection on the three-dimensional resistivity distribution results to identify spatial locations where the resistivity spatial change rate exceeds a preset gradient threshold. Specifically, for the three-dimensional resistivity distribution data volume, the resistivity change gradient of each grid cell in three orthogonal directions is calculated. The gradient calculation uses the central difference method, which divides the resistivity difference between adjacent grid cells by the center-point spacing of the grid cells. For each grid cell, the square root of the sum of the squares of its gradient components in the three directions is taken to obtain the total amplitude of the resistivity spatial change rate at that cell. A preset gradient threshold is used to determine whether the change is significant enough to constitute a boundary. This threshold is set based on a statistical analysis of the gradient amplitude values ​​of the entire 3D resistivity distribution data volume. First, the average and standard deviation of the gradient amplitude values ​​of all grid cells are calculated. Then, the preset gradient threshold is set as the average value plus an adjustable multiple multiplied by the standard deviation. The specific value of this multiple is determined by analyzing the statistical distribution of the gradient amplitude values. For example, the cumulative frequency distribution curve of the gradient amplitude values ​​can be observed, and a multiple such that the preset gradient threshold corresponds to the high percentile value of the curve, such as 80% to 90%, can be selected. In this way, the preset gradient threshold is set at a level that highlights significant changes while suppressing background noise. The total amplitude value of the calculated resistivity spatial change rate in all grid cells is compared with the preset gradient threshold. The 3D coordinates of the center points of grid cells whose total amplitude value is greater than the preset gradient threshold are extracted, forming a set containing numerous spatial coordinate points. This set is the initial electrical abrupt change boundary point set.

[0076] A density-based spatial clustering algorithm is applied to the initial set of electrical mutation boundary points to divide the area into different spatial clusters based on the density distribution of spatial locations. Then, for each spatial cluster, the following operations are performed: the clustering confidence score is calculated; if the clustering confidence score is higher than a preset confidence threshold, it is identified as a candidate electrical mutation boundary and further validated by spectral density analysis; otherwise, it is ignored. After processing all spatial clusters, the finally validated boundaries are added to the set. The implementation of the spatial clustering algorithm requires setting two key parameters: the neighborhood search radius and the minimum number of points threshold. The neighborhood search radius defines the radius of the neighborhood search around a point in three-dimensional space. This radius is set based on the spatial statistical characteristics of the points in the initial electrical mutation boundary point set. It is calculated by averaging the distance from each point in the set to its nearest k-th neighbor, and then setting the neighborhood search radius to 2 to 3 times this average distance, where k can be an integer between 5 and 10. The minimum point count threshold specifies the minimum number of points a dense region should contain. Its setting considers the overall size of the initial electrical mutation boundary point set; for example, the minimum point count threshold can be set between 1% and 5% of the total number of points, but not less than a minimum cardinality to ensure statistical significance, such as 5 points. During algorithm execution, an unvisited point is randomly selected from the initial electrical mutation boundary point set. Then, using this point as the center and a set neighborhood search radius as the radius, all points within that neighborhood are searched. If the number of points in the neighborhood is greater than or equal to the minimum point count threshold, a new spatial cluster is created, and the core point and all points in its neighborhood are added to this cluster. The same neighborhood search and expansion process is then recursively performed on each point in the neighborhood until no more points can be added, thus forming a complete spatial cluster. If the number of points in the neighborhood is less than the minimum point count threshold, the point is temporarily marked as a noise point. After traversing all points in the point set, points not assigned to any cluster are ultimately considered noise points and removed. Through this process, the initial set of electrical abrupt boundary points is divided into several spatially dense sets of points, each set being called a spatial cluster. The algorithm also records the number of points contained in each spatial cluster and the spatial distribution range of these points.

[0077] Calculate the cluster confidence score for each spatial cluster. Cluster confidence score is a comprehensive quantification of the compactness and significance of a spatial cluster's spatial distribution. Its calculation is based on two attributes of the spatial cluster: the first attribute is the density of points within the cluster, calculated by dividing the total number of points in the cluster by its approximate volume in space. This approximate volume is estimated by multiplying the ranges of the three-dimensional coordinates of the points within the cluster in each direction. The second attribute is the regularity of the spatial distribution shape of the points within the cluster, measured by calculating the eigenvalues ​​of the covariance matrix of the point coordinates within the cluster and the ratio of the largest to the smallest eigenvalue to measure the degree of anisotropy. A specific calculation method involves multiplying the normalized density value (the cluster density divided by the largest density among all clusters) by a factor reflecting regularity, such as subtracting the reciprocal of the eigenvalue ratio from 2 to limit it to between 1 and 2, and then multiplying by a scaling factor that adjusts the result to the range of 0 to 1. The final result is a value between 0 and 1, which is then used as the cluster confidence score for the spatial cluster. The specific value of the scaling factor can be determined empirically by testing on typical datasets, so that the confidence value has good discriminative power.

[0078] Spatial clusters with cluster confidence scores higher than a preset confidence threshold are identified as candidate electrical mutation boundaries. The preset confidence threshold is used to filter out clusters that are spatially significant and reliable. This threshold is set based on an analysis of the distribution of cluster confidence scores for all calculated spatial clusters. For example, the cluster confidence scores of all spatial clusters are sorted in ascending order, their mean and standard deviation are calculated, and the preset confidence threshold is set to the mean plus 0.5 times the standard deviation, or a fixed value such as 0.6 can be set empirically. Through this comparison, only spatial clusters with calculated cluster confidence scores greater than the preset confidence threshold are retained. Each retained cluster represents a spatially continuous and concentrated potential electrical mutation boundary, i.e., a candidate electrical mutation boundary.

[0079] For each candidate electrical abrupt change boundary, continuous resistivity values ​​are extracted along its spatial extension direction to form a resistivity fluctuation sequence. Specifically, for a spatial cluster corresponding to a candidate electrical abrupt change boundary, a principal axis direction that best represents its extension direction is first determined based on the spatial distribution of its points. Principal component analysis of the three-dimensional coordinates of points within the cluster is used to determine the first principal component direction as the principal axis direction. Then, along this principal axis direction, a series of sampling points are set at spatial intervals smaller than the size of a three-dimensional grid cell, for example, an interval set to half the size of the grid cell. For each sampling point, the resistivity value at that point is calculated using a three-dimensional linear interpolation method based on the three-dimensional resistivity distribution data volume. The resistivity values ​​of all sampling points arranged sequentially along the principal axis direction are connected according to the sampling order to form a discrete sequence with the sampling point sequence number as the independent variable and the resistivity value as the dependent variable. This sequence is the resistivity fluctuation sequence.

[0080] Power spectral density estimation is performed on a resistivity fluctuation sequence. Power spectral density estimation is used to analyze the intensity of different frequency components in the resistivity fluctuation sequence. In practice, the resistivity fluctuation sequence is first detrended, i.e., the linear trend term is removed, making its mean zero. Then, the periodogram method is used for estimation. The sequence is treated as a discrete-time signal, and its discrete Fourier transform is calculated to obtain a series of complex Fourier coefficients. The square of the modulus of each Fourier coefficient is then calculated and divided by the sequence length to obtain the power spectral density estimate for the corresponding frequency point. These frequency points are determined by the sampling interval and the sequence length, ranging from zero to half the sampling frequency.

[0081] Candidate electrical abrupt change boundaries exhibiting significant spectral peaks within the frequency range determined by resistivity fluctuation sequence feature analysis are included in a spatiotemporally validated set of electrical abrupt change boundaries. The frequency range determined by resistivity fluctuation sequence feature analysis is specifically defined as follows: First, the estimated power spectral density of the resistivity fluctuation sequence is calculated. Then, the average value of the entire estimated power spectral density is calculated as the average background level. Next, local maxima points with spectral density values ​​significantly higher than the average background level are identified on the power spectrum; these points are the peak values. A peak judgment threshold is set for identifying significant peaks, for example, the average background level plus twice its standard deviation. Frequency points with power spectral density values ​​greater than this peak judgment threshold are identified, and it is examined whether these frequency points form continuous intervals on the frequency axis. These continuous intervals are used as the frequency range determined by resistivity fluctuation sequence feature analysis. For a candidate electrical abrupt change boundary, the power spectral density estimate of its resistivity fluctuation sequence is examined to determine whether there is at least one obvious peak within the frequency range defined by the resistivity fluctuation sequence feature analysis. This peak value must not only be higher than the peak judgment threshold but also significantly higher than the spectral density values ​​of its adjacent frequency points. Typically, the peak height is required to exceed 50% of the average height of the valleys on both sides. Candidate electrical abrupt change boundaries that meet this condition are considered to have passed both spatial clustering and sequence spectral feature tests and are thus included in the final spatiotemporally validated set of electrical abrupt change boundaries. This set contains electrical boundaries that are spatially clustered and exhibit significant abrupt spectral changes in resistivity along the boundary extension direction, providing reliable input for distinguishing between real defects and imaging artifacts.

[0082] S5. Based on the spatiotemporally verified set of electrical abrupt change boundaries, a similarity matching algorithm is used to remove continuous pseudo-abnormal transition regions, generating the purified resistivity distribution results. The specific implementation is as follows:

[0083] After obtaining a spatiotemporally validated set of electrical abrupt change boundaries, a similarity matching algorithm is used to eliminate continuous false anomaly transition zones that may be generated during the inversion process, thereby generating a purified resistivity distribution result that more closely approximates the actual geological conditions. This process first identifies all spatially adjacent electrical abrupt change boundary pairs within the spatiotemporally validated set. Specifically, for each electrical abrupt change boundary in the set, the three-dimensional coordinates of its spatial geometric center point are calculated. This center point is obtained by calculating the average coordinates of all points within the spatial cluster corresponding to the boundary, i.e., calculating the arithmetic mean of the X, Y, and Z coordinates of all points to obtain the coordinates representing the center point of the boundary. Then, the Euclidean distance between the center points of any two different electrical abrupt change boundaries in the set is calculated, i.e., the sum of the squares of the coordinate differences between the two points in the X, Y, and Z directions is calculated, and then the square root of the sum of squares is taken to obtain the distance value. To determine whether two boundaries are spatially adjacent, a distance threshold for adjacency judgment needs to be set. This threshold is based on an analysis of the statistical distribution of distances between all possible boundary pairs in a spatiotemporally validated set of electrically abrupt boundary pairs. For example, the average distance of all boundary pairs can be calculated, and the distance threshold for adjacency judgment can be set between 1.5 and 2 times this average. The specific multiple is determined by observing the histogram of distance value distribution, selecting a multiple that keeps the number of adjacent boundary pairs within a reasonable range. The calculated distance between the center points of any two boundaries is compared with the set distance threshold for adjacency judgment. If the distance is less than or equal to the distance threshold, the two electrically abrupt boundaries are determined to form a pair of spatially adjacent electrically abrupt boundary pairs. By traversing all possible boundary combinations, all boundary pairs that satisfy the adjacency condition are identified, and each pair of boundaries is recorded.

[0084] The resistivity distribution similarity between the regions defined by each pair of adjacent abrupt electrical boundary changes is calculated. First, the spatial distribution data of resistivity within the regions defined by each pair of adjacent abrupt electrical boundary changes needs to be extracted. Specifically, for each pair of adjacent abrupt electrical boundary changes, a spatial region between the two boundaries is defined in the 3D resistivity distribution results. This region can be approximated as a 3D space enclosed by the minimum outer envelope cuboids of the two boundaries. The resistivity values ​​of all grid cells within this 3D space are extracted to form a 3D resistivity data sub-block. The resistivity distribution similarity of this region is calculated by calculating the correlation coefficient or structural similarity index between the spatial resistivity distribution data in different directions or profiles. If the correlation coefficient method is used, multiple parallel two-dimensional profiles are cut along multiple different spatial directions or cut into the three-dimensional data sub-block, such as along the direction connecting the center points of two boundaries, and several profiles perpendicular to this direction. A series of continuous resistivity values ​​are extracted from each profile to form a sequence. Then, the Pearson correlation coefficient between any two different direction or profile sequences is calculated. The Pearson correlation coefficient is calculated by dividing the covariance of the two sequences by the product of their respective standard deviations. Finally, the average of all the calculated correlation coefficient values ​​is taken as the resistivity distribution similarity of the region defined between the pair of adjacent boundaries. If the structural similarity index method is used, resistivity value sequences from different directions or cross-sections are obtained. Each pair of sequences is treated as two one-dimensional "images". The product of their brightness comparison value, contrast comparison value, and structural comparison value is calculated. The brightness comparison value is a function of the mean of the two sequences, the contrast comparison value is a function of the standard deviation of the two sequences, and the structural comparison value is a function of the product of the covariance and standard deviation of the two sequences. Then, these three comparison values ​​are combined, usually by multiplying the three, to obtain an index between -1 and 1. This index value is the structural similarity index of the two sequences. Finally, the arithmetic mean of the index values ​​of all paired sequences is calculated as the resistivity distribution similarity of the region.

[0085] Regions with resistivity distribution similarity exceeding a preset similarity threshold are identified as continuous pseudo-anomaly transition zones. The preset similarity threshold is used to determine whether the region between two real anomaly boundaries exhibits overly similar, smoothly transitioning electrical characteristics, potentially indicating artificial generation by the inversion algorithm rather than genuine geological structures. This preset similarity threshold is set based on the analysis of numerous forward simulation cases. For example, multiple theoretical models containing known discrete anomalies are constructed for inversion, calculating the similarity value distribution between boundary pairs in these real cases. Simultaneously, a homogeneous medium model is constructed and subjected to strong smoothing constraints for inversion to generate typical pseudo-transition zones, calculating their similarity value distribution. The overlapping regions of these two types of distributions are analyzed, and a value located at the edge of the overlapping region that can distinguish most pseudo-transition zone cases is selected as the preset similarity threshold. For example, 0.75 can be chosen as the initial reference value based on analysis. The resistivity distribution similarity of each pair of adjacent boundaries obtained in the previous calculation step is compared with the preset similarity threshold. If the calculated similarity value of a certain pair of boundaries is higher than the preset similarity threshold, the three-dimensional spatial region between that pair of boundaries is determined to be a continuous pseudo-anomaly transition zone.

[0086] In the 3D resistivity distribution results, the resistivity values ​​identified as continuous pseudo-anomaly transition zones are smoothed and corrected to generate a purified resistivity distribution result. The smoothing correction operation is performed on every grid cell within the 3D spatial region identified as a continuous pseudo-anomaly transition zone. One smoothing correction method is the region mean replacement method, which calculates the average resistivity of a normal background region immediately surrounding the continuous pseudo-anomaly transition zone. This background region is defined as a shell region extending outward from the outer surface of the transition zone with a certain thickness, for example, three times the size of the 3D grid cell. Then, the calculated average resistivity of the background region is used to replace the original resistivity values ​​of all grid cells within the transition zone. Another method is spatial filtering, such as using a 3D Gaussian filter, which filters only the resistivity values ​​within the transition zone. The size of the Gaussian filter kernel is set according to the spatial scale of the transition zone; for example, the standard deviation of the kernel in the three directions is set to one-sixth of the width of the transition zone in that direction. The smoothing correction process needs to be performed cell-by-cell, updating the resistivity values ​​at the corresponding positions in the 3D resistivity distribution results. After smoothing and correcting all identified continuous pseudo-anomaly transition zones, artificially created regions with gradual electrical transitions that might have connected two real anomalies in the original 3D resistivity distribution were corrected to resistivity distributions more consistent with their surrounding background, resulting in a purified resistivity distribution. This result preserves the spatiotemporally validated boundaries of real electrical abrupt changes and the anomalies they enclose, while suppressing or removing pseudo-transition zones that might be introduced by the inversion algorithm due to smoothness constraints. This makes the characterization of the internal inhomogeneity of cement mixing piles more accurate and reliable, providing a cleaner data foundation for the final quality evaluation.

[0087] S6. Based on the resistivity distribution results after purification, output the conclusion on the uniformity of cement mixing pile hole formation. The specific implementation is as follows:

[0088] After obtaining the resistivity distribution results after purification, the uniformity of the cement mixing pile borehole is quantitatively evaluated based on these results, and the final test conclusion is output. This process first extracts the resistivity values ​​of the core area of ​​the pile from the purified resistivity distribution results. Specifically, based on the spatial range of the core area of ​​the pile predefined in the inversion calculation, this range is a cylindrical space extending from the designed pile top elevation to the designed pile bottom elevation, with the designed pile center as the axis and the designed pile diameter as the diameter. In the three-dimensional mesh data volume corresponding to the purified resistivity distribution results, all mesh cells whose center point coordinates lie within this cylindrical space are identified. All the resistivity values ​​stored in these identified mesh cells are extracted to form a set containing multiple resistivity values; this set represents the resistivity values ​​of the core area of ​​the pile. These resistivity values ​​represent the true estimates reflecting the electrical distribution of the pile material obtained after the aforementioned series of purification treatments.

[0089] Calculate the statistical coefficient of variation (COP) of resistivity values ​​in the core area of ​​the pile. The COP is a dimensionless statistic used to quantify the relative dispersion of resistivity values ​​within the core area, thus characterizing overall uniformity. The calculation process involves three steps: First, calculate the arithmetic mean of the resistivity values ​​in the core area. This involves summing all resistivity values ​​in the set and then dividing the sum by the total number of resistivity values. Second, calculate the standard deviation of the resistivity values ​​in the core area. This is done by first calculating the difference between each resistivity value and the arithmetic mean obtained in the previous step, then squaring each difference, summing all the squares, dividing this sum by the total number of resistivity values, and finally taking the square root of the quotient to obtain the standard deviation. Third, divide the calculated standard deviation by the arithmetic mean; the quotient is the COP. A larger COP indicates greater relative fluctuation in resistivity within the core area and poorer uniformity; conversely, a smaller COP indicates better uniformity.

[0090] The spatial location and resistivity contrast of the anomalous bodies corresponding to the spatiotemporally verified set of electrically abrupt boundary changes are extracted. Here, an anomalous body refers to a region enclosed or marked by the boundaries in the spatiotemporally verified set of electrically abrupt boundary changes, exhibiting a significant electrical difference from the surrounding medium. For each electrically abrupt boundary in the set of spatiotemporally verified electrically abrupt boundary changes, the spatial location of its corresponding anomalous body is determined by calculating the three-dimensional coordinates of the geometric center point of the region enclosed by that boundary. These center point coordinates can be obtained by calculating the average of the coordinates of the center points of all grid cells within the boundary. Resistivity contrast is used to quantify the degree of electrical difference between the anomalous body and the surrounding background medium. When calculating resistivity contrast, a representative resistivity value inside the anomalous body is first determined; for example, the arithmetic mean of the resistivity values ​​of all grid cells within the anomalous body is taken as the anomalous body resistivity. Next, representative resistivity values ​​of the region immediately adjacent to the anomaly are determined. In the purified resistivity distribution results, the background region is defined as a shell region extending outward from the anomaly boundary with a specific thickness. The thickness of this shell region is set according to the size of the three-dimensional mesh cells, for example, three times the average size of the mesh cells. The resistivity values ​​of all mesh cells within this shell region are extracted, and their arithmetic mean is calculated as the background resistivity. Finally, the resistivity contrast is calculated as the absolute value of the difference between the anomaly resistivity and the background resistivity, divided by the background resistivity to obtain a dimensionless ratio representing the relative difference. For the spatiotemporally verified set of electrical abrupt boundary conditions, the spatial coordinates and resistivity contrast values ​​of the anomaly corresponding to each boundary are calculated one by one, forming a set of quantitative data on the spatial distribution and intensity characteristics of the anomaly.

[0091] The statistical coefficient of variation (SCVA) and the spatial location and resistivity contrast of anomalies are compared with a pre-defined uniformity evaluation standard. This standard is a pre-established set of quantitative and qualitative criteria for judging the uniformity of cement-mixed piles, including multiple specific numerical thresholds and spatial constraints. The thresholds are set based on statistical analysis of a large number of historical engineering case data, theoretical understanding of the physical properties of cement-soil materials, and the combined results of forward numerical simulation experiments. For example, the SCVA threshold for judging overall uniformity is set as follows: a batch of cement-mixed pile cases verified as having acceptable mixing uniformity through traditional core sampling and mechanical testing are collected. The SCVA of each case is calculated using this testing method. Then, the distribution range of the SCVA values ​​for all these acceptable cases is calculated, such as their mean and standard deviation. The SCVA threshold for judging acceptable uniformity is set to the mean of the SCVA of acceptable cases plus twice the standard deviation, and the SCVA threshold for judging excellent uniformity is set to the mean of the SCVA of acceptable cases plus one standard deviation. To determine the threshold for resistivity contrast of anomalous bodies, a theoretical geological model containing anomalous bodies with different contrasts was constructed for forward simulation and inversion analysis using this method. The ability of anomalous bodies to be reliably detected and identified under different contrast levels was observed, and the lower limit of contrast that could be stably identified and corresponded to unacceptable defects in engineering was used as the determination threshold. Regarding the constraints on the spatial location of anomalous bodies, key limiting areas, such as the pile top loading area and the area of ​​maximum pile bending moment, were determined based on pile stress analysis and construction quality control specifications, using engineering experience.

[0092] The pre-defined uniformity evaluation criteria typically include several specific criteria. The first criterion concerns the statistical coefficient of variation (SCV), which includes a threshold for acceptable SCV and a threshold for excellent SCV. For example, a SCV of less than or equal to 0.15 indicates excellent overall uniformity of the pile; a SCV greater than 0.15 and less than or equal to 0.25 indicates acceptable uniformity; and a SCV greater than 0.25 indicates unacceptable uniformity. The second criterion concerns the spatial location of anomalous bodies, including a series of spatial constraints. For example, no spatiotemporally verified anomalous bodies are allowed within the upper third of the pile depth; and within the middle third of the pile depth, the lateral distribution of anomalous bodies must not exceed 50% of the central area of ​​the pile diameter. The third criterion concerns the resistivity contrast of anomalous bodies, including a resistivity contrast threshold for excessively strong defects. For example, the resistivity contrast of a single anomalous body must not exceed 1.5; otherwise, it is considered to have excessively strong local defects. The comparison operation involves comparing the actual calculated statistical coefficient of variation (SCV) value with the corresponding SCV thresholds for determining compliance and excellence in the standard. It also checks whether the spatial coordinates of the actual anomaly conform to the spatial constraints specified in the standard and whether the resistivity contrast value of each anomaly exceeds the resistivity contrast threshold specified in the standard for determining excessive defects.

[0093] Based on the comparison results, a conclusion is generated regarding the uniformity of the cement mixing pile borehole. The conclusion is generated based on a comprehensive logical judgment of the above comparison results. If the calculated statistical coefficient of variation (SCV) is less than or equal to the threshold for excellent performance, and all detected anomalies fully meet all the limiting clauses in the preset uniformity evaluation criteria in terms of spatial location and resistivity contrast, the conclusion is that the uniformity of the cement mixing pile borehole is excellent. If the SCV is greater than the threshold for excellent performance but less than or equal to the threshold for acceptable performance, and all anomalies meet the criteria, the conclusion is that the uniformity is acceptable. If the SCV is greater than the threshold for acceptable performance, or if one or more anomalies violate spatial constraints or their resistivity contrast exceeds the threshold for excessive defects, the conclusion is that the uniformity of the cement mixing pile borehole is unacceptable. For boundary cases, such as a SCV slightly exceeding the threshold for acceptable performance, or individual anomaly parameters slightly exceeding limits, the conclusion is that there are minor defects in uniformity, requiring comprehensive judgment based on the specific engineering conditions. The test results are output in the form of a written report. The report must list the specific values ​​of the calculated statistical coefficient of variation and their comparison with the thresholds used to determine whether the test is acceptable or excellent; list the spatial coordinates of all detected anomalies and their compliance with spatial constraints; list the resistivity contrast values ​​of all anomalies and their comparison with the resistivity contrast threshold used to determine excessive defects; and, based on a comprehensive analysis of the above comparisons, provide a clear qualitative conclusion. This conclusion, based on the purified resistivity data and verified anomaly boundary information, significantly improves the accuracy and reliability of judging the internal mixing uniformity of cement mixing piles.

[0094] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0095] It should be noted that this invention can be deployed on the device itself to realize embedded applications, or it can run on a PC or other terminal with a user interface, thereby meeting various hardware environments and usage requirements.

[0096] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. A computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, all or part of the processes or functions according to the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. Computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wireless or wired transmission; wired transmission methods include optical fiber, twisted pair, coaxial cable, etc.; wireless transmission includes infrared, microwave, etc. Computer-readable storage media can be any available medium that a computer can access or a data storage device such as a server or data center that contains one or more sets of available media. Available media can be magnetic media (e.g., floppy disks, hard disks, magnetic tapes), optical media (e.g., DVDs), or semiconductor media. Semiconductor media can be solid-state drives.

[0097] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and modules described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0098] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0099] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical modules; they may be located in one place or distributed across multiple network modules. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs.

[0100] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0101] If a function is implemented as a software module and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0102] The above are merely specific embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0103] In conclusion, the above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for detecting the uniformity of borehole formation in cement mixing piles based on three-dimensional resistivity imaging, characterized in that, include: S1. Obtain resistivity measurement data of cement mixing piles and surrounding medium; S2. The resistivity measurement data is filtered using an anomaly detection algorithm to obtain valid resistivity data; S3. Based on the effective resistivity data, a partition constraint inversion strategy is adopted to divide the core area of ​​the pile into the surrounding medium area. The smooth constraint is weakened in the core area of ​​the pile, while the conventional smooth constraint is retained in the surrounding medium area to generate a three-dimensional resistivity distribution result. S4. Extract the electrical abrupt boundary from the three-dimensional resistivity distribution results, use density clustering algorithm to analyze the spatial distribution clustering confidence of the electrical abrupt boundary, and screen out the electrical abrupt boundary with a confidence level higher than the preset confidence threshold. Then, use spectral density estimation algorithm to perform time-series spectral density analysis on the resistivity fluctuation of the screened electrical abrupt boundary to generate a set of spatiotemporally verified electrical abrupt boundaries. S5. Based on the set of electrical abrupt change boundaries verified by time and space, a similarity matching algorithm is used to remove continuous false abnormal transition regions and generate the purified resistivity distribution results. S6. Based on the resistivity distribution results after purification, output the test conclusion on the uniformity of the cement mixing pile hole formation.

2. The method for detecting the uniformity of cement mixing pile borehole formation based on three-dimensional resistivity imaging according to claim 1, characterized in that, S1 includes: A multi-ring electrode array is deployed around the cement mixing pile; A high-density electrical resistivity measurement mode combining multiple electrode spacing and multiple orientations was adopted to collect apparent resistivity data of cement mixing piles and surrounding media under different combinations of power supply electrodes and measuring electrodes, which were used as resistivity measurement data.

3. The method for detecting the uniformity of cement mixing pile borehole formation based on resistivity three-dimensional imaging according to claim 1, characterized in that, S2 include: The overall distribution characteristics of resistivity measurement data were analyzed based on statistical methods. Determine the threshold for identifying abnormal data based on the overall distribution characteristics; Measurements that deviate from the overall distribution characteristics by more than a threshold in the resistivity measurement data are considered outliers and are removed. The remaining data constitute the valid resistivity data.

4. The method for detecting the uniformity of cement mixing pile borehole formation based on three-dimensional resistivity imaging according to claim 1, characterized in that, S3 includes: The spatial range of the core area and the surrounding medium area of ​​the cement mixing pile is defined based on the design location and geometric dimensions of the cement mixing pile; In the inversion calculation, the part belonging to the core area of ​​the pile is subject to a minimization constraint based on the first norm or an anisotropic smoothing constraint with reduced lateral smoothing weight, while the part belonging to the surrounding medium area is subject to a minimization constraint based on the second norm. A three-dimensional inversion calculation is performed using the effective resistivity data as the fitting target to generate a three-dimensional resistivity distribution result.

5. The method for detecting the uniformity of cement mixing pile borehole formation based on three-dimensional resistivity imaging according to claim 4, characterized in that, Using the effective resistivity data as the fitting target, a three-dimensional inversion calculation is performed to generate the three-dimensional resistivity distribution result. Specifically, an inversion objective function containing data fitting terms and partition smoothing constraint terms is constructed. The resistivity parameters in the three-dimensional space are adjusted through an iterative algorithm to minimize the fitting difference between the forward simulation data and the effective resistivity data until the convergence condition is met, and the final three-dimensional resistivity distribution result is output.

6. The method for detecting the uniformity of cement mixing pile borehole formation based on resistivity three-dimensional imaging according to claim 1, characterized in that, S4 include: Edge detection is performed on the three-dimensional resistivity distribution results to identify spatial locations where the spatial rate of resistivity change exceeds a preset gradient threshold, thus forming an initial set of electrical abrupt change boundary points. A density-based spatial clustering algorithm is applied to the initial set of boundary points of electrical abrupt change. Different spatial clusters are divided according to the density distribution of spatial location points, and the clustering confidence of each spatial cluster is calculated. Spatial clusters with clustering confidence scores higher than a preset confidence threshold are identified as candidate electrical mutation boundaries; For each candidate electrical abrupt change boundary, continuous resistivity values ​​are extracted along its spatial extension direction to form a resistivity fluctuation sequence, and the power spectral density of the resistivity fluctuation sequence is estimated. Candidate electrical abrupt change boundaries that exhibit significant spectral peak characteristics within the frequency range determined by resistivity fluctuation sequence feature analysis are included in the spatiotemporally validated set of electrical abrupt change boundaries.

7. The method for detecting the uniformity of cement mixing pile borehole formation based on three-dimensional resistivity imaging according to claim 6, characterized in that, The frequency range determined based on the resistivity fluctuation sequence characteristic analysis is specifically determined in the following way: calculate the power spectral density of the resistivity fluctuation sequence, identify the peaks in the power spectrum whose spectral density values ​​are significantly higher than the average background level, and determine the frequency range covered by the peaks.

8. The method for detecting the uniformity of cement mixing pile borehole formation based on three-dimensional resistivity imaging according to claim 1, characterized in that, S5 include: In a set of spatiotemporally validated electrical abrupt change boundaries, identify spatially adjacent pairs of electrical abrupt change boundaries; Calculate the similarity of resistivity distribution in the regions defined by each pair of adjacent electrical abrupt boundary lines; Regions with resistivity distribution similarity higher than a preset similarity threshold are identified as continuous pseudo-abnormal transition regions; In the three-dimensional resistivity distribution results, the resistivity values ​​identified as continuous pseudo-abnormal transition zones are smoothed and corrected to generate the purified resistivity distribution results.

9. The method for detecting the uniformity of cement mixing pile borehole formation based on three-dimensional resistivity imaging according to claim 8, characterized in that, The resistivity distribution similarity between the regions defined by each pair of adjacent electrical abrupt boundary is calculated by extracting the resistivity spatial distribution data within the regions defined by each pair of adjacent electrical abrupt boundary, and quantifying the resistivity distribution similarity between each pair of adjacent electrical abrupt boundary by calculating the correlation coefficient or structural similarity index of the resistivity spatial distribution data in different directions or profiles.

10. The method for detecting the uniformity of cement mixing pile borehole formation based on resistivity three-dimensional imaging according to claim 1, characterized in that, S6 include: Extract the resistivity value of the core area of ​​the pile from the resistivity distribution results after purification; Calculate the statistical coefficient of variation of resistivity values ​​in the core area of ​​the pile, and extract the spatial location and resistivity contrast of the anomalies corresponding to the set of electrical abrupt boundary verified in time and space. The statistical coefficient of variation and the spatial location of the anomaly are compared with the resistivity contrast to the preset uniformity evaluation criteria. Based on the comparison results, a conclusion is generated regarding the test results of the uniformity of the drilling of cement mixing piles.

Images