Underground ore body positioning and identifying method and system based on ground penetrating radar

Abstract

The application provides a ground penetrating radar-based underground ore body positioning and identification method and system, and relates to the technical field of ore body detection.The method comprises the following steps: inputting a corrected multi-dimensional point cloud feature set into a pre-trained closed-loop convolutional neural network, optimizing network parameters through transfer learning, and outputting a category confidence of an ore body target and a spatial position probability distribution point cloud; based on the spatial position probability distribution point cloud, combining electromagnetic wave double-path travel time calculation rules and Cagniard resistivity measurement results, and generating an ore body depth and horizontal positioning point cloud through a reverse projection algorithm; and performing multi-scale two-dimensional empirical mode decomposition on the reverse projection positioning point cloud, filtering noise components to extract a high-confidence target point cloud.The application improves the accuracy and stability of ore body identification.

Description

Technical Field

[0001] This invention relates to the field of ore body detection technology, and in particular to a method and system for locating and identifying underground ore bodies based on ground-penetrating radar. Background Technology

[0002] In a lead-zinc mine exploration project in the Nanling polymetallic mineral cluster, the exploration team used traditional processing methods to locate and identify underground ore bodies, but some technical defects were exposed. The geological conditions in this area are complex, with multiple sets of strike faults cutting through the area. The ore bodies are distributed in a vein-like pattern under the control of the structure, and a large number of quartz veins are developed in the surrounding rock, resulting in the mixing of ore body reflection and non-ore interference signals in the radar signal.

[0003] In practical applications, strong reflection signals from fractured zones were misidentified as thick ore bodies. This is because traditional methods fail to extract and analyze the sidelobe interference characteristics in radar signals, making it impossible to distinguish between the true reflection of the ore body and the scattering interference from the fault fractured zone, resulting in uncorrected reflection deviations. Furthermore, since the ore body depth varies between 50 and 150 meters, the electromagnetic wave propagation attenuation characteristics differ significantly at different depths. Traditional methods do not calculate reflection fluctuation values ​​and fail to capture the intensity variation patterns of reflection signals from adjacent point cloud units, leading to a lack of dynamic correction coefficients. Consequently, weak reflection signals from deep ore bodies were masked by noise, resulting in three instances of missed detection of deeply buried ore bodies, impacting the accuracy of subsequent drilling engineering design. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a method and system for locating and identifying underground ore bodies based on ground penetrating radar, which improves the accuracy and stability of ore body identification.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0006] Firstly, a method for locating and identifying underground ore bodies based on ground-penetrating radar, the method comprising:

[0007] The original radar signal is preprocessed to generate a 3D point cloud dataset containing spatial coordinates and reflection intensity. Each point cloud data unit contains a triplet of position coordinates, echo delay, and signal amplitude. Based on the 3D point cloud dataset, the goodness of fit of the reflection intensity sequence and sidelobe interference features of each point cloud unit are extracted, and the reflection deviation is calculated. Based on the 3D point cloud dataset, the spatial distribution of the binary tuple consisting of the maximum reflected signal value and orientation angle of adjacent point cloud units is extracted, and the reflection fluctuation value is calculated. The reflection deviation and reflection fluctuation value are fused to generate dynamic correction coefficients, which are used to enhance the anti-interference performance of the reflection intensity of the target point cloud subset, resulting in a corrected multi-dimensional point cloud feature set. The corrected multi-dimensional... The point cloud feature set is input into a pre-trained closed-loop convolutional neural network. The network parameters are optimized through transfer learning, and the output is a point cloud showing the category confidence and spatial location probability distribution of the ore body target. Based on the spatial location probability distribution point cloud, combined with the electromagnetic wave two-way travel time calculation rules and the Carnia resistivity measurement results, the depth and horizontal positioning point cloud of the ore body is generated through the inverse projection algorithm. Multi-scale two-dimensional empirical mode decomposition is performed on the inverse projection positioning point cloud to filter out noise components and extract high-confidence target point clouds. Based on the high-confidence target point cloud, the ore body vertex coordinates are corrected through a hyperbola fitting algorithm, and combined with the electromagnetic wave propagation speed calculation rules in the underground medium, a three-dimensional spatial morphology distribution point cloud of the ore body is generated.

[0008] Secondly, the underground ore body location and identification system based on ground-penetrating radar includes:

[0009] The preprocessing module is used to preprocess the raw radar signal to generate a three-dimensional point cloud dataset containing spatial coordinates and reflection intensity. Each point cloud data unit contains a triplet of position coordinates, echo delay, and signal amplitude.

[0010] The calculation module is used to extract the goodness of fit of the reflection intensity sequence and the sidelobe interference features of each point cloud unit based on the 3D point cloud dataset, and to calculate the reflection deviation.

[0011] The extraction module is used to extract the spatial distribution of the binary tuple consisting of the maximum value and direction angle of the reflection signal of adjacent point cloud units based on the 3D point cloud dataset, and to calculate the reflection fluctuation value.

[0012] The fusion module is used to fuse reflection deviation and reflection fluctuation values, generate dynamic correction coefficients, enhance the anti-interference of the reflection intensity of the target point cloud subset, and obtain a corrected multi-dimensional point cloud feature set.

[0013] The prediction module is used to input the corrected multi-dimensional point cloud feature set into a pre-trained closed-loop convolutional neural network, optimize the network parameters through transfer learning, and output the category confidence and spatial location probability distribution point cloud of the ore body target.

[0014] The processing module is used to generate a point cloud for ore body depth and horizontal positioning based on the spatial location probability distribution point cloud, combined with the electromagnetic wave two-way travel time calculation rules and the Carnia resistivity measurement results, through a reverse projection algorithm; perform multi-scale two-dimensional empirical mode decomposition on the reverse-projected positioning point cloud to filter out noise components and extract high-confidence target point cloud; based on the high-confidence target point cloud, correct the ore body vertex coordinates through a hyperbola fitting algorithm, and combine the electromagnetic wave propagation speed calculation rules in the underground medium to generate a three-dimensional spatial morphology distribution point cloud of the ore body.

[0015] The above-described solution of the present invention has at least the following beneficial effects:

[0016] By using a pre-trained closed-loop convolutional neural network combined with transfer learning, and leveraging the powerful feature learning and pattern recognition capabilities of neural networks, ore body targets are accurately identified from the corrected multi-dimensional point cloud feature set. The output category confidence score provides a reliable quantitative evaluation basis for the ore body identification results. By filtering out noise components through multi-scale two-dimensional empirical mode decomposition, high-confidence target point clouds are further extracted, reducing the interference of false targets on the identification results and significantly improving the accuracy and stability of ore body identification. Attached Figure Description

[0017] Figure 1 This is a schematic flowchart of the underground ore body location and identification method based on ground penetrating radar provided in an embodiment of the present invention.

[0018] Figure 2 This is a schematic diagram of an underground ore body positioning and identification system based on ground-penetrating radar provided in an embodiment of the present invention. Detailed Implementation

[0019] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0020] like Figure 1 As shown, embodiments of the present invention propose a method for locating and identifying underground ore bodies based on ground-penetrating radar, the method comprising the following steps:

[0021] Step 100: Preprocess the original radar signal to generate a three-dimensional point cloud dataset containing spatial coordinates and reflection intensity, wherein each point cloud data unit contains a triplet of position coordinates, echo delay and signal amplitude.

[0022] Step 200: Based on the 3D point cloud dataset, extract the goodness of fit of the reflection intensity sequence and the sidelobe interference features of each point cloud unit, and calculate the reflection deviation.

[0023] Step 300: Based on the three-dimensional point cloud dataset, extract the spatial distribution of the binary tuple formed by the maximum value of the reflected signal and the direction angle of the adjacent point cloud units, and calculate the reflection fluctuation value.

[0024] Step 400: Fuse reflection deviation and reflection fluctuation values ​​to generate dynamic correction coefficients, enhance the anti-interference effect on the reflection intensity of the target point cloud subset, and obtain the corrected multi-dimensional point cloud feature set;

[0025] Step 500: Input the corrected multi-dimensional point cloud feature set into the pre-trained closed-loop convolutional neural network, optimize the network parameters through transfer learning, and output the category confidence and spatial location probability distribution point cloud of the ore body target.

[0026] Step 600: Based on the spatial location probability distribution point cloud, combined with the electromagnetic wave two-way travel time calculation rules and the Carnia resistivity measurement results, the ore body depth and horizontal positioning point cloud are generated through the inverse projection algorithm.

[0027] Step 700: Perform multi-scale two-dimensional empirical mode decomposition on the inverse projection positioning point cloud to filter out noise components and extract high-confidence target point cloud.

[0028] Step 800: Based on the high-confidence target point cloud, the coordinates of the ore body vertices are corrected by the hyperbola fitting algorithm, and the three-dimensional spatial distribution point cloud of the ore body is generated by combining the calculation rules of the propagation speed of electromagnetic waves in the underground medium.

[0029] In this embodiment of the invention, by extracting reflection deviation and fluctuation values ​​and generating dynamic correction coefficients, the anti-interference ability of point cloud data is effectively enhanced and the quality of basic data is improved. By using closed-loop convolutional neural networks and transfer learning, combined with multi-scale noise filtering, the accuracy and reliability of ore body identification are significantly improved. By integrating electromagnetic wave characteristics and resistivity data and optimizing through inverse projection, hyperbola fitting, etc., high precision in ore body spatial positioning and three-dimensional morphology characterization is achieved.

[0030] In a preferred embodiment of the present invention, step 100 includes:

[0031] Step 101: Time-frequency conversion is performed on the original radar signals collected from the polymetallic mineral cluster area to generate a two-dimensional echo signal sequence image. The two-dimensional echo signal sequence image contains the temporal and amplitude information of the reflected signals caused by fault cutting and quartz vein penetration. Specifically, at the lead-zinc mine exploration site in the Nanling polymetallic mineral cluster area, the ground-penetrating radar equipment continuously emits electromagnetic waves and receives reflected signals along the preset survey line. The original radar signals are stored in the form of continuous waveform data in the time domain, which records the voltage amplitude changes of the reflected signals at different times. In order to extract the temporal patterns and amplitude differences of geological features such as fault cutting and quartz vein penetration contained in the reflected signals, the original radar signals need to be time-frequency converted. The specific process is as follows: The acquired raw time-domain signal is divided into several signal segments at fixed time intervals, each segment corresponding to the received signal of the radar antenna at a certain instant on the survey line; each time-domain signal segment is decomposed into different frequency components through time-frequency analysis, and the amplitude of each frequency component at the corresponding time point is determined; with time as the horizontal axis and frequency as the vertical axis, the amplitude of different frequency components at each time point is represented by grayscale or color intensity to generate a two-dimensional echo signal sequence image; in this image, the changes in grayscale or color brightness directly reflect the amplitude information of the reflected signal, while the time shift on the horizontal axis corresponds to the temporal process of the radar antenna moving along the survey line, completely preserving the temporal distribution characteristics of reflected signals from geological bodies such as faults and quartz veins.

[0032] Step 102: Map the temporal information of each pixel in the two-dimensional echo signal sequence image to the spatial position coordinates on the survey line, while retaining the echo time axis information of the corresponding pixel, and generating a two-dimensional data field with spatial coordinates and original time delay. Specifically, in the two-dimensional echo signal sequence image, the temporal information of the horizontal axis essentially corresponds to the movement process of the radar antenna on the survey line, and converts it into intuitive spatial position coordinates. The specific operation is as follows: First, obtain the radar antenna's movement parameters, including the coordinates of the survey line's starting point, the antenna's moving speed, and the signal acquisition time interval, and determine the antenna's moving distance corresponding to each time point. Taking the survey line's starting point as the origin and the direction along the survey line as the coordinate axis, convert each time point on the horizontal axis of the two-dimensional echo signal sequence image into the actual spatial position coordinates on the survey line according to the rule: spatial position coordinates = antenna moving speed × signal acquisition time interval × time point number + starting point coordinates. During the conversion process, retain the echo time axis information on the vertical axis of the two-dimensional echo signal sequence image, that is, the propagation time of the electromagnetic wave from transmission to reception. Use the converted spatial position coordinates as the new horizontal axis and the original echo time axis as the vertical axis, and retain the signal amplitude values ​​of each coordinate point accordingly to form a two-dimensional data field containing spatial position coordinates, original echo time, and signal amplitude. This data field realizes the conversion of temporal information into spatial position, enabling the image information to establish a preliminary correspondence with the spatial distribution of underground geological bodies.

[0033] Step 103, based on the propagation speed of electromagnetic waves in the complex medium of the mining area, converts the echo time axis information in the two-dimensional data field with spatial coordinates and original time delay into depth time delay values, specifically including:

[0034] Existing geological data of the mining area were collected, and basic information such as the stratigraphic division scheme, known ore body numbers and occurrence strata, and the occurrence of major fault structures were compiled. Simultaneously, geophysical exploration results from the past five years were retrieved, including gravity and magnetic exploration anomaly maps and interpretation reports, to determine the correspondence between physical anomalies and geological bodies. The collected data were digitized and archived, a basic data index was established, and key data sources and reliability levels were marked. A research team composed of geological engineers and geophysical engineers conducted field surveys according to the principles of grid layout and key tracking. Surface geological mapping was carried out along 10 pre-set exploration lines (500-meter spacing), recording the exposed lithology, structural fracture zones, and quartz vein outcrops. GPS was used for precise positioning (error controlled within 5 meters) to create a distribution map of surface geological phenomena. For fault fracture zone outcrops, their strike, dip, and dip angle were measured, and fault gouge and breccia samples were collected within the fracture zones, with the sampling location coordinates marked.

[0035] For the 20 boreholes already drilled in the mining area, complete core files were retrieved, and core logging was performed at 5-meter intervals according to borehole depth. The thickness of the ore body, mineral composition, and contact relationship with the surrounding rock were recorded in the cores. Microscopic identification was used to determine the lithology of the surrounding rock (such as the mineral composition and structure of limestone and sandstone). For core sections containing quartz veins, the vein thickness and penetration direction were measured, and the distribution frequency of the veins at different depths was statistically analyzed. The core data were entered into 3D geological modeling software according to the borehole coordinates to generate a preliminary 3D distribution model of the ore body, surrounding rock, and fault fracture zone. The distribution characteristics of the geological bodies were verified by combining ground magnetic and resistivity exploration data. Three geophysical profiles were re-measured at the ore body outcrop area and fault zone location to compare the spatial correspondence between physical property anomalies and geological bodies. For example, the low resistivity anomaly characteristics and extension range of the fault fracture zone were confirmed by high-density resistivity profiles, and the boundary coordinates of the geological bodies were corrected. Finally, a report on the spatial distribution characteristics of geological bodies in the mining area was generated, determining the planar distribution range, vertical extension depth, and mutual contact relationship of various geological bodies.

[0036] Thirty representative samples were selected from field-collected specimens, including ore bodies (galena-sphalerite assemblage), limestone host rocks, sandstone host rocks, fault breccia, and quartz veins. The samples were processed into standard cylinders with a diameter of 5 cm and a height of 10 cm. Their relative permittivity was measured using a dielectric constant meter in a laboratory environment (temperature 25℃, humidity 50%). Based on the conversion relationship between electromagnetic wave propagation speed and dielectric constant, the electromagnetic wave propagation speed of each sample was calculated, and the arithmetic mean of samples of the same type was taken as the reference value for the indoor experimental speed of that geological body.

[0037] Five typical test points were selected in the mining area, and a 200-meter-long ground-penetrating radar line was set up at each test point. Data acquisition was carried out using a 250MHz antenna, with a sampling rate of 200MHz and a stacking number of 16. The acquired radar data was preprocessed (DC drift removal and gain correction), and velocity analysis was performed using the common center point method (CMP). A clear reflection interface with the same phase axis was selected on the radar profile, and the two-way travel time corresponding to different offset distances was calculated. The time-distance curve was plotted and the velocity value was fitted. For complex areas such as fault fracture zones and areas with dense quartz veins, the borehole radar test method was used. The antenna was placed in a borehole (50-100 meters deep) of a known geological body, and the propagation speed of electromagnetic waves in the geological body was calculated by the radar reflection signal in the borehole.

[0038] Comparing the velocity results from indoor experiments and field tests, if the deviation exceeds 10%, five more sets of samples are added for repeated experiments, or three additional field survey lines are used for data supplementation. For key geological bodies such as ore bodies and fault fracture zones, velocity calibration is performed based on the geological conditions revealed by known boreholes. For example, at borehole locations with known ore body thickness, the accurate velocity is retrieved by comparing the radar profile reflection time with the actual thickness. The electromagnetic wave propagation velocity range and representative mean value for each geological body are ultimately determined. Based on the core logging results, the stratification boundaries of the medium are delineated in the three-dimensional geological model. For example, the surface 0-5 meters is classified as the Quaternary overburden, 5-80 meters (limestone area) as the upper surrounding rock, and depths above 80 meters as bedrock. For fault fracture zones, stratification boundaries are defined according to their actual influence range (width 5-20 meters). Quartz veins are classified into corresponding depth strata according to their penetration depth (40-100 meters). The vertical error of the stratification boundaries is controlled within ±3 meters.

[0039] Based on the stratification boundaries, the previously measured velocity values ​​of geological bodies are matched to the corresponding strata. For strata dominated by a single geological body, the average velocity of that body is directly adopted. If a stratum contains two or more geological bodies (such as a fault zone and a contact zone with the surrounding rock), the velocity value is calculated by weighting the proportion of the geological bodies within the stratum. For example, in a 50-60 meter depth stratum, 70% is limestone (velocity 65 m / ns) and 30% is quartz vein (velocity 75 m / ns), then the overall velocity of this stratum is 65 × 70% + 75 × 30% = 68 m / ns. A velocity database is constructed using relational database software (such as MySQL), and fields such as geological body type, stratum number, depth range, velocity range, average value, data source, and spatial coordinates are entered. The database is spatially linked to the mining area topographic map using GIS software to enable velocity data queries based on coordinate ranges. Three geological experts were organized to verify the database content, focusing on checking the rationality of the matching between velocity values ​​and geological body types, the accuracy of stratification boundaries, and correcting any problems found (such as abnormal velocity ranges). The final result was a database of electromagnetic wave propagation velocity in mining areas that includes spatial indexes.

[0040] For the two-dimensional data field with spatial coordinates and original time delay generated in step 102, each original echo time point in the data field is traversed. For each original echo time point, its corresponding spatial coordinates are extracted, including the horizontal coordinates of the survey line direction and the horizontal coordinates of the direction perpendicular to the survey line. Based on the extracted spatial coordinates, the previously established database of electromagnetic wave propagation velocity in the mining area is queried to determine the geological body type or medium stratification corresponding to the spatial coordinates. If the location is in a single geological body, the electromagnetic wave propagation velocity of that geological body is directly matched. If the location involves multiple medium stratifications and the electromagnetic wave propagation path passes through different stratifications, the electromagnetic wave propagation velocities of each stratification are weighted according to the actual thickness ratio of each stratification in the propagation path to obtain the comprehensive electromagnetic wave propagation velocity corresponding to the original echo time point.

[0041] For each original echo time point, the time value of its original echo time axis is determined. This time value represents the two-way propagation time of the electromagnetic wave from the surface emission point, through the underground target location, and back to the surface receiving point. Based on this original echo time value and the matched electromagnetic wave propagation speed, the two-way distance of electromagnetic wave propagation is calculated according to the inherent relationship between the distance, time, and speed of electromagnetic wave propagation in the medium. Since the depth delay value needs to reflect the underground depth, the calculated two-way propagation distance is divided by 2 to obtain the one-way distance of the electromagnetic wave propagating from the surface to the underground target location. This one-way distance is the underground depth reflected by the corresponding original echo time point, and this depth value is determined as the depth delay value for that original echo time point.

[0042] After converting the depth delay values ​​of all original echo time points in the two-dimensional data field, the vertical axis parameter of the two-dimensional data field is replaced. The original echo time axis, which is in units of time, is replaced with a depth delay value axis, which is in units of depth. The spatial coordinates and corresponding signal amplitude values ​​of each point in the data field remain unchanged; only the parameter type and value of the vertical axis are updated, thus forming a new two-dimensional data field. This new data field contains three key parameters: spatial coordinates, depth delay value, and signal amplitude. It achieves a quantitative conversion from the original echo time to underground depth, enabling the two-dimensional data field to more intuitively reflect the reflection signal characteristics of geological bodies at different underground depths.

[0043] Step 104: Extract the signal amplitude corresponding to each spatial location coordinate point, bind it with the depth delay value and spatial coordinate of the point as a triplet, and generate an anti-interference three-dimensional point cloud dataset. Specifically, in order to construct a three-dimensional data structure containing spatial location, depth and reflection intensity, it is necessary to extract the signal amplitude of the two-dimensional data field and bind it with triplet. The specific steps are as follows: Traverse each spatial coordinate point in the two-dimensional data field, and extract the signal amplitude value of the point at the corresponding depth delay value. This amplitude value directly reflects the reflection intensity of the geological body at that spatial location and depth to electromagnetic waves. The greater the reflection intensity, the more likely there is a ore body, fault, or quartz vein or other reflecting interface at that location. For each spatial coordinate point, associate and bind its spatial coordinates (including the coordinates of the survey line direction and the horizontal coordinates perpendicular to the survey line), the depth delay value obtained by step 103, and the extracted signal amplitude value to form a triplet data unit of spatial coordinates, depth delay value, and signal amplitude. Arrange and combine all triplet data units according to the distribution rules of their spatial coordinates and depth delay values ​​to construct a three-dimensional data set, namely, an anti-interference three-dimensional point cloud dataset. This dataset completely preserves the spatial location information and reflection intensity characteristics of each point underground through the triplet form.

[0044] In this embodiment of the invention, a two-dimensional echo signal sequence image is generated through time-frequency conversion, which fully preserves the temporal patterns and amplitude differences of the reflected signals from geological phenomena such as fault cutting and quartz vein interpenetration, avoiding feature loss caused by direct processing of the original signal; the temporal information is mapped to the spatial coordinates of the survey line, realizing the transformation of radar signals from the time domain to the spatial domain, establishing an intuitive correspondence between the reflected signal and the actual spatial distribution of underground geological bodies, solving the problem of the difficulty in locating spatial positions with traditional two-dimensional signals; combined with the electromagnetic wave propagation velocity database of complex media in mining areas, the echo time is accurately converted into depth delay values, eliminating the influence of differences in different geological media on depth calculation, making the signal depth information more consistent with the actual underground geological structure, and improving the accuracy of depth positioning; an anti-interference three-dimensional point cloud dataset is generated by binding triples, integrating the core information of spatial position, depth delay, and reflection intensity to form a structured three-dimensional data structure, effectively reducing the impact of noise and interference on feature extraction.

[0045] In a preferred embodiment of the present invention, step 200 includes:

[0046] Step 201: Extract the reflection intensity sequence of each point cloud unit from the 3D point cloud dataset; Step 202: Fit an ideal distribution curve based on the reflection intensity sequence, and calculate the variance of the fitting residual as a goodness-of-fit feature; Step 203: Based on the same reflection intensity sequence, locate the main lobe peak point of the reflection intensity sequence, and search for the first trough point on both sides of the main lobe peak point as the side lobe boundary; calculate the amplitude ratio based on the highest point within the main lobe peak point and the side lobe boundary as the side lobe peak attenuation rate, and measure the width of a continuous region within the side lobe boundary whose amplitude value is 10% higher than the main lobe peak value as the side lobe width feature; Step 204: Combine the goodness-of-fit feature with the side lobe peak attenuation rate and the side lobe width feature, and calculate the reflection deviation by weighted summation.

[0047] In this embodiment of the invention, step 201 may include determining each point cloud unit to be processed from the three-dimensional point cloud dataset generated in step 100, wherein each point cloud unit contains position coordinates, echo delay, and signal amplitude triplet information as defined in step 104; for each determined point cloud unit, the signal amplitude value reflecting the electromagnetic wave reflection intensity is extracted from the triplet data of that unit. Since each point cloud unit will form a series of continuous reflection intensity records under different echo delays, these signal amplitude values ​​arranged in the order of echo delay are extracted sequentially to form the reflection intensity sequence corresponding to that point cloud unit.

[0048] Step 202 above may include determining the ideal distribution curve based on the reflection intensity sequence of each point cloud unit extracted in step 201. Since the distribution characteristics of the ground-penetrating radar reflection signal are determined by the physical laws of electromagnetic wave propagation in the underground medium, its ideal shape needs to be determined by combining the natural attenuation characteristics of electromagnetic wave reflection with the measured characteristics of interference-free signals. According to the physical characteristics of electromagnetic wave propagation in a homogeneous medium, when electromagnetic waves encounter a homogeneous interface and are reflected, the signal intensity is usually symmetrically distributed around the peak value, and gradually attenuates according to a fixed law as it moves away from the peak value. Based on this, the basic morphological characteristics (symmetrical attenuation trend) of the ideal reflection signal in the echo delay dimension are determined; the target is selected. Within the target detection area, for a known uniform geological section free from ore body interference, ground-penetrating radar is used to collect reflection signals from this area. Typical reflection sequences free from interference are extracted, and key characteristics such as peak position, attenuation rate, and symmetry are analyzed as a reference for the measured ideal curve shape. Combining the physical attenuation law of electromagnetic wave propagation with the characteristics of the measured interference-free signal, the specific shape of the ideal distribution curve (such as a symmetrical bell-shaped distribution centered on the peak) is determined. Based on the attenuation amplitude range of the measured interference-free signal, the core parameter range of the curve is initially set (including the echo delay interval corresponding to the peak, the maximum amplitude interval, the attenuation slope range, etc.), thereby completing the selection and determination of the ideal distribution curve.

[0049] The ideal distribution curve can be approximated by a Gaussian curve (normal distribution curve) to characterize the ideal distribution pattern of the reflection intensity sequence. The formula is as follows: ;in, This represents the theoretical value of the reflection intensity at position x on the ideal distribution curve. Here, x is the one-dimensional position coordinate of the point cloud unit in space (such as distance or depth along the detection direction), corresponding to the data point position in the reflection intensity sequence; A represents the peak amplitude of the ideal distribution curve, corresponding to the highest theoretical intensity value of the main lobe in the reflection intensity sequence, reflecting the strongest energy level of the reflected signal under ideal conditions. The peak position (mean) of the ideal distribution curve corresponds to the theoretical position coordinates of the main lobe peak point in the reflection intensity sequence, which is the spatial position where the reflected signal energy is most concentrated. It forms a corresponding reference with the actual position of the main lobe peak point that needs to be located in step 203. The width parameter (standard deviation) of the ideal distribution curve; It is a natural constant.

[0050] Based on the aforementioned determined Gaussian curve (normal distribution curve) as the ideal distribution curve, a fitting calculation is performed between the reflection intensity sequence and the Gaussian ideal distribution curve. For the reflection intensity sequence of each point cloud unit extracted in step 201, where each data point contains the correspondence between echo delay (corresponding to the position variable in the Gaussian curve) and signal amplitude (corresponding to the reflection intensity value in the Gaussian curve), this correspondence is matched with the determined Gaussian ideal distribution curve. First, the core parameter range of the Gaussian ideal distribution curve is used as the initial value, where the core parameters include the reference range of peak amplitude, the theoretical interval of peak position, and the empirical value of the width parameter. The shape of the Gaussian curve under the initial parameters is compared point by point with the actual data points of the reflection intensity sequence. During the comparison process, the focus is on the degree of overlap between the peak position of the actual data point (i.e., the echo delay corresponding to the maximum signal amplitude) and the peak position of the Gaussian curve, as well as the consistency between the attenuation trend of the actual data point on both sides of the peak with the change of echo delay and the theoretical attenuation trend of the Gaussian curve, so as to preliminarily determine the matching basis between the curve and the actual data.

[0051] Secondly, based on the results of the point-by-point comparison above, the key parameters of the Gaussian curve are adjusted successively: If there is a significant shift between the peak position of the actual reflection intensity sequence and the current peak position of the Gaussian curve, i.e., the echo delay values ​​corresponding to the two are significantly different, then the parameter representing the peak position in the Gaussian curve is adjusted so that the peak point of the curve moves towards the echo delay position of the actual peak; If the attenuation rate of the actual data points on both sides of the peak (i.e., the rate at which the signal amplitude decreases with the increase of the echo delay) does not match the theoretical attenuation trend of the Gaussian curve, for example, the actual attenuation is too fast or too slow, then the width parameter controlling the steepness of attenuation in the Gaussian curve is adjusted so that the attenuation rhythm of the curve gradually approaches the attenuation rhythm of the actual data; If the attenuation symmetry deviation of the actual reflection intensity sequence on the left and right sides of the peak is large, i.e., the signal amplitude difference under the same echo delay offset on both sides is significant, then the attenuation balance on both sides of the curve is adjusted by fine-tuning the equivalent effect strength of the width parameter on both sides of the curve, so that the curve shape gradually conforms to the actual distribution characteristics of the reflection intensity sequence in the overall trend.

[0052] After multiple rounds of parameter adjustments, when the Gaussian curve achieves optimal matching with the overall characteristics of the reflection intensity sequence in three key dimensions—that is, when the peak position of the curve basically coincides with the peak echo delay of the actual data, the attenuation trend of the curve is consistent with the attenuation rate of the actual data, the symmetry of the left and right sides of the curve matches the symmetry characteristics of the actual data, and the overlap between the curve and the actual data points is highest from a visual perspective—then parameter adjustments are stopped, and the fitting calculation between the reflection intensity sequence and the ideal Gaussian distribution curve is completed.

[0053] After the fitting calculation of the reflection intensity sequence and the Gaussian ideal distribution curve is completed in step 202, the calculation stage of the fitting residuals begins. For each actual data point in the reflection intensity sequence, since each data point contains corresponding echo delay information and signal amplitude value, it is necessary to accurately select a point with the same echo delay position as the actual data point from the already fitted ideal distribution curve, and read the corresponding amplitude value of that point on the ideal distribution curve, i.e., the fitted value of the curve under the same echo delay. Then, the difference between the signal amplitude value of the actual data point and the curve fitting value at the corresponding position is calculated. Specifically, the signal amplitude value of the actual data point is subtracted from the fitted value of the curve at that echo delay position. The resulting difference is the fitting residual of the actual data point. By performing the above operation sequentially on all actual data points in the reflection intensity sequence, the fitting residual corresponding to each data point can be obtained. These residuals completely reflect the local deviation of the actual reflected signal from the ideal distribution curve at each echo delay position.

[0054] After obtaining all the fitting residuals, statistical analysis is performed on these residuals to obtain the goodness-of-fit characteristics of the point cloud unit. The specific process is as follows:

[0055] Square all the obtained fitting residuals, that is, calculate the square value for each residual value. Squaring can amplify the influence weight of residuals with large deviations on the overall result, making the effect of significant deviations more prominent in subsequent analysis. After squaring all residuals, calculate the arithmetic mean of all residual squares. In this calculation, sum all residual squares and then divide the sum by the total number of residuals (i.e., the total number of data points in the reflection intensity sequence). The result is the arithmetic mean of all residual squares. This mean can comprehensively and quantitatively reflect the overall deviation between all actual data points and the ideal distribution curve. The arithmetic mean calculated above is determined as the variance of the fitting residuals. This variance value is determined as the goodness-of-fit feature of the point cloud unit, providing a key feature parameter for the calculation of reflection deviation in the subsequent step 204.

[0056] Step 203 above may include: based on the same reflection intensity sequence extracted in step 201, comparing all signal amplitude values ​​in the sequence, and finding the point with the largest value, which is the main lobe peak point of the reflection intensity sequence, and the corresponding signal amplitude value is the main lobe peak value; taking this main lobe peak point as the center, successively searching for adjacent data points to the left and right of the reflection intensity sequence, finding the first turning point where the signal amplitude value starts to decrease and then increase again, which is the sidelobe boundary; the one searched on the left is the left sidelobe boundary, and the one searched on the right is the right sidelobe boundary. After determining the left and right sidelobe boundaries, finding the points with the largest signal amplitude values, i.e., the highest points of the left and right sidelobes, within the range between the left sidelobe boundary and the main lobe peak point, and between the right sidelobe boundary and the main lobe peak point, respectively, and recording their corresponding amplitude values. Then, the ratio of the amplitude value of the highest point of the left sidelobe to the main lobe peak value is used as the left sidelobe peak attenuation rate, and the ratio of the amplitude value of the highest point of the right sidelobe to the main lobe peak value is used as the right sidelobe peak attenuation rate. Meanwhile, within the left sidelobe boundary range, from the left sidelobe boundary towards the main lobe peak point, the echo delay range corresponding to all continuous data points with signal amplitude values ​​higher than 10% of the main lobe peak value is measured, and the length of this range is the left sidelobe width characteristic; similarly, within the right sidelobe boundary range, the echo delay range length corresponding to the continuous region with signal amplitude values ​​higher than 10% of the main lobe peak value is measured, and this is used as the right sidelobe width characteristic.

[0057] Step 204 above may include integrating the goodness-of-fit features obtained in step 202 with the sidelobe peak attenuation rates (including left and right sidelobe peak attenuation rates) and sidelobe width features (including left and right sidelobe width features) obtained in step 203, and using these features as the basic input parameters for calculating the reflection deviation; then, according to the importance of each feature in reflecting the reflection signal deviation, assigning a corresponding weight coefficient to each feature; specifically, it is necessary to combine the detection environment of the ground penetrating radar (such as the homogeneity of the underground medium, the strength of the interference signal, etc.) and the characteristics of the ore body reflection signal (such as the stability of the reflection intensity, the typical manifestation of sidelobe interference, etc.), and pre-setting these weight coefficients through the statistical analysis of previous experimental data and theoretical analysis, and the sum of the weight coefficients of each feature is 1.0, the specific setting rules are as follows:

[0058] When the underground medium in the detection area is highly homogeneous (such as dense rock layers and low-fracture geological structures) and the interference signal is weak (no obvious electromagnetic scattering, and the proportion of clutter signal is less than 5%), the reflected signal is less affected by the medium in such an environment. The overall deviation of the signal from the ideal distribution curve is the main factor reflecting the deviation of the reflected signal. Therefore, the goodness-of-fit feature has the highest correlation with the deviation of the reflected signal, and its weight coefficient ranges from 0.4 to 0.5. The sidelobe interference has a weak impact on the signal deviation, and the weight coefficients of the left sidelobe peak attenuation rate and the right sidelobe peak attenuation rate both range from 0.1 to 0.15. The sidelobe coverage area contributes the least to the deviation, and the weight coefficients of the left sidelobe width feature and the right sidelobe width feature both range from 0.05 to 0.08.

[0059] When the underground medium in the detection area is heterogeneous (such as loose soil or fractured rock layers) and the interference signal is strong (significant electromagnetic scattering, with clutter signals accounting for 5% to 20%), in such environments, the overall signal distortion and sidelobe interference jointly affect the deviation. The correlation of the goodness-of-fit feature decreases slightly, with the weighting coefficient ranging from 0.35 to 0.45. The influence of sidelobe interference is enhanced, with the weighting coefficients for the peak attenuation rate of the left and right sidelobes ranging from 0.13 to 0.17. The contribution of sidelobe width to the deviation increases slightly, with the weighting coefficients for the left and right sidelobe width features ranging from 0.07 to 0.09.

[0060] When the ore body reflection intensity is stable (signal amplitude fluctuation coefficient <10%) but sidelobe interference is typical (e.g., high conductivity ore body causing sidelobe energy ratio >30%), sidelobe interference is the main source of signal deviation in this scenario. The sidelobe peak attenuation rate has the highest correlation with the reflection signal deviation, and the weighting coefficients for the left and right sidelobe peak attenuation rates are both in the range of 0.17~0.2. The overall signal fit is relatively high, and the weighting coefficients for the goodness-of-fit feature are in the range of 0.3~0.35. The sidelobe width has a moderate impact on the deviation, and the weighting coefficients for the left and right sidelobe width features are both in the range of 0.07~0.09.

[0061] When the sidelobes of the ore body's reflected signal are large (sidelobe coverage depth exceeds twice the main lobe width), the sidelobe width significantly interferes with signal recognition in adjacent areas. The correlation between the left and right sidelobe width features increases, with their weighting coefficients ranging from 0.09 to 0.1. The weighting coefficient for the goodness-of-fit feature ranges from 0.35 to 0.4, and the weighting coefficient for the sidelobe peak attenuation rate ranges from 0.13 to 0.17. Based on the specific conditions of the detection environment and the characteristics of the ore body's reflected signal, weighting coefficients for each feature are selected from the corresponding ranges, ensuring that the sum of all feature weighting coefficients is 1.0. Simultaneously, the weighting coefficients are fine-tuned using previously measured data from known ore body areas (e.g., comparing the deviation of the calculated reflection deviation under different weighting combinations with the deviation verified by actual borehole drilling). Finally, weighting coefficients matching the current detection scenario are determined to ensure that the contribution of each feature to the reflection deviation accurately matches its actual impact.

[0062] For example, in scenarios with tight rock layers and weak sidelobe interference, the weighting coefficients can be set as follows: goodness-of-fit feature 0.45, left sidelobe peak attenuation rate 0.15, right sidelobe peak attenuation rate 0.15, left sidelobe width feature 0.08, and right sidelobe width feature 0.07. The sum of these weighting coefficients is 1.0, which conforms to the correlation rules of each feature in this scenario. Subsequently, based on the set weighting coefficients, the weighted values ​​of each feature are calculated: the goodness-of-fit feature obtained in step 202 is multiplied by its corresponding weighting coefficient to obtain the weighted value of the goodness-of-fit feature; the left sidelobe peak attenuation rate obtained in step 203 is multiplied by its corresponding weighting coefficient to obtain the weighted value of the left sidelobe peak attenuation rate; similarly, the right sidelobe peak attenuation rate is multiplied by its corresponding weighting coefficient to obtain the weighted value of the right sidelobe peak attenuation rate; The left sidelobe width feature obtained in step 203 is multiplied by its corresponding weight coefficient to obtain the weighted value of the left sidelobe width feature. Similarly, the right sidelobe width feature is multiplied by its corresponding weight coefficient to obtain the weighted value of the right sidelobe width feature. The weighted values ​​of all features calculated above (including the weighted values ​​of the goodness-of-fit feature, the left sidelobe peak attenuation rate, the right sidelobe peak attenuation rate, the left sidelobe width feature, and the right sidelobe width feature) are summed. The sum obtained is the reflection deviation of the point cloud unit. This reflection deviation integrates the overall deviation of the signal from the ideal state reflected by the goodness-of-fit feature, as well as the influence of sidelobe interference reflected by the sidelobe peak attenuation rate and the sidelobe width feature, and fully characterizes the deviation level of the reflected signal of the point cloud unit.

[0063] In this embodiment, by extracting the goodness-of-fit features (reflecting the overall deviation of the signal from the ideal state), sidelobe peak attenuation rate (reflecting the sidelobe interference intensity), and sidelobe width features (reflecting the sidelobe interference range) of the reflection intensity sequence, a multi-dimensional characterization of the reflection signal deviation is achieved, covering key sources of deviation such as signal distortion, interference intensity, and interference range. The reflection deviation is calculated by weighted summation and fusion of multiple features, with the weighting coefficients dynamically set in combination with the detection environment and ore body signal characteristics, ensuring that the contribution of each interference factor to the deviation is accurately quantified under different scenarios, avoiding misjudgment of deviation caused by the one-sidedness of a single feature. Through the quantitative analysis of the reflection deviation, the differences between the effective signal of the ore body and interference signals such as noise and sidelobes are effectively distinguished. The feature extraction and deviation calculation process is dynamically adjusted in combination with the physical laws of electromagnetic wave propagation and the actual detection scenario, which can adapt to different geological environments (such as medium homogeneity and differences in interference intensity) and ore body reflection characteristics (such as signal stability and sidelobe performance), enhancing the robustness of the method in complex underground environments.

[0064] In a preferred embodiment of the present invention, step 300 includes:

[0065] Step 301: Extract the maximum value of the reflection intensity sequence of each point cloud unit from the three-dimensional point cloud dataset; Step 302: Calculate the direction vector of each point cloud unit and its neighboring units based on the spatial coordinates of the point cloud units, and extract the direction angle; Step 303: Combine the maximum value of the reflection intensity sequence of each point cloud unit with the direction angle to generate a binary spatial distribution dataset; Step 304: Based on the binary spatial distribution dataset, calculate the product of the standard deviation of the maximum reflection intensity between adjacent units and the rate of change of the direction angle, as the reflection fluctuation value.

[0066] In this embodiment of the invention, step 301 may include, based on the three-dimensional point cloud dataset generated in step 100, locating the reflection intensity sequence contained in each point cloud unit in the dataset. This sequence consists of signal amplitude values ​​determined in step 201 and arranged in echo delay order. The reflection intensity sequence of each point cloud unit is traversed, and the magnitudes of all signal amplitude values ​​in the sequence are compared. The signal amplitude value with the largest value is selected, which is the maximum value of the reflection intensity sequence of that point cloud unit, used to characterize the strongest energy level of the reflected signal of that unit. Through the above operations, the maximum value of the reflection intensity sequence of all point cloud units is extracted.

[0067] Step 302 above may include: based on the spatial coordinates (including X and Y plane coordinates and Z coordinates corresponding to depth delay) of each point cloud unit in the 3D point cloud dataset generated in step 100, setting a reasonable spatial neighborhood range (such as a distance threshold set based on the detection accuracy) with the spatial coordinates of the current point cloud unit as the center, and determining all other point cloud units within this range as its neighboring units; subtracting the spatial coordinates of the current point cloud unit from the spatial coordinates of the neighboring units to obtain vector data representing the relative positional relationship between the two, and the direction of this vector reflects the spatial orientation of the neighboring unit relative to the current unit; taking the X-axis of the 3D coordinate system as the reference direction, calculating the angle between the direction vector and the reference direction, and determining the angle value as the orientation angle of the neighboring unit relative to the current point cloud unit, thereby quantifying the spatial orientation relationship.

[0068] Step 303 above may include, based on the extraction of the maximum reflection intensity sequence of each point cloud unit in step 301 and the acquisition of the orientation angle in step 302, binding the maximum reflection intensity sequence of the point cloud unit with the orientation angle of each of its corresponding adjacent units to form a binary data set of maximum reflection intensity - orientation angle, wherein each binary set corresponds to a set of association information of the strongest energy of the current unit - the orientation of the adjacent units; using the spatial coordinates of the point cloud unit as an index, summarizing all the binary data corresponding to each coordinate position to form a binary spatial distribution dataset covering the entire detection area, which fully records the association features between the maximum reflection intensity and the adjacent orientation at different spatial positions.

[0069] Step 304 above may include: based on the binary spatial distribution dataset generated in step 303, extracting the maximum value of the reflection intensity sequence for each point cloud unit and all its adjacent units, calculating the standard deviation of these maximum values ​​through statistical analysis. This standard deviation reflects the dispersion of the strongest energy of reflection intensity between the current unit and its adjacent units; the larger the standard deviation, the more significant the intensity fluctuation. Calculating the rate of change of azimuth angle between adjacent units: for the azimuth angle of the current point cloud unit and its adjacent units, calculating the ratio of the difference in azimuth angle between adjacent units to the corresponding spatial distance, obtaining the rate of change of azimuth angle. This rate of change reflects the speed of change in spatial orientation; the larger the rate of change, the more obvious the orientation difference between adjacent units. Multiplying the standard deviation of the maximum reflection intensity obtained above with the rate of change of azimuth angle, the result is the reflection fluctuation value of the point cloud unit. This value comprehensively quantifies the correlation between the spatial fluctuation degree of reflection intensity and orientation change, and can effectively reflect the non-uniformity of the reflection characteristics of the underground medium.

[0070] This embodiment, by binding the maximum reflection intensity with the azimuth angle of adjacent units into a binary spatial distribution dataset, overcomes the limitations of single-unit feature analysis, achieving accurate capture of the correlation between reflection signal intensity and spatial orientation, and effectively reflecting the spatial distribution law of reflection characteristics. By calculating the reflection fluctuation value through the product of standard deviation (reflecting intensity dispersion) and azimuth angle change rate (reflecting the rate of orientation change), the synergistic influence of spatial differences in reflection intensity and orientation changes is comprehensively quantified, accurately characterizing the non-uniformity of underground medium reflection characteristics. By focusing on the relative characteristics (maximum intensity, azimuth angle) of adjacent units, the interference of single-point noise on the overall analysis is reduced, while spatial correlation is used to filter random interference signals, improving the ability to identify the spatial continuity of ore body reflection signals.

[0071] In a preferred embodiment of the present invention, step 400 includes:

[0072] Step 401: Based on the reflection deviation and reflection fluctuation value, a dynamic correction coefficient is generated through nonlinear function mapping; Step 402: According to the gradient distribution of the reflection deviation, point cloud units with reflection deviations greater than a preset threshold are selected to form a target point cloud subset; Step 403: For each unit in the target point cloud subset, its reflection intensity is multiplied by the corresponding dynamic correction coefficient to obtain the enhanced reflection intensity; Step 404: The enhanced reflection intensity is recombined with the original spatial coordinates and echo delay to generate a corrected multi-dimensional point cloud feature set.

[0073] In this embodiment of the invention, the formula corresponding to the nonlinear function in step 401 above can be:

[0074] ;

[0075] in, This represents the dynamic correction coefficient, the magnitude of which reflects the correction strength of the reflection intensity. The larger the k value, the stronger the correction strength; the closer the k value is to 1, the weaker the correction strength (close to no correction). The value of d represents the reflection deviation calculated in step 200, reflecting the degree of deviation of the point cloud unit's reflected signal from the ideal state and the level of sidelobe interference. The larger the d value, the more severe the interference to the signal, and the stronger the correction required. This represents the reflection fluctuation value calculated in step 300, reflecting the spatial non-uniformity of reflection intensity between adjacent point cloud units. The larger the value, the more drastic the change in spatial reflection characteristics, requiring targeted enhancement and correction; This represents the correction intensity adjustment parameter, used to control the maximum increase of the correction coefficient. Its value range is usually 0.5~2.0 (which can be dynamically adjusted according to the detection environment). The larger the value, the larger the maximum possible value of the correction coefficient, and the wider the adjustment range of the overall correction force; This represents a sensitivity parameter used to control the rate at which reflection deviation and fluctuation values ​​affect the correction coefficient; its value typically ranges from 0.1 to 1.0. The larger the value, the more sensitive the change in reflection deviation and fluctuation value is to the influence of the correction coefficient. The correction coefficient increases rapidly with the increase of deviation and fluctuation. It is a natural constant.

[0076] Step 402 above may include, based on the reflection deviation obtained in step 200, using the spatial coordinates of the 3D point cloud dataset as a reference, generating a gradient distribution map reflecting the spatial trend of deviation by analyzing the rate of change of reflection deviation of adjacent point cloud units. This map can intuitively present the spatial distribution characteristics of high deviation areas. A preset threshold for reflection deviation is set, which is determined based on previous detection experience or interference level assessment. Typically, a critical value that can distinguish between effective signals and strong interference signals is selected (for example, referring to the maximum deviation value of the effective signal of the ore body in historical data, the threshold is set to be slightly higher than this value). Target point cloud subsets are selected based on the gradient distribution map: all point cloud units are traversed, and point cloud units with reflection deviations greater than the preset threshold are selected. These units are usually areas that are severely interfered with and require key correction, thus constituting the target point cloud subset that needs anti-interference enhancement.

[0077] Step 403 above may include: performing a dynamic correction coefficient matching operation on the target point cloud subset (i.e., the set of point cloud units whose reflection deviation is greater than a preset threshold and requires key correction) selected in step 402; traversing each point cloud unit in the subset, associating the dynamic correction coefficient generated for it in step 401 with the unit's unique identifier (such as spatial coordinate index), ensuring that the correction coefficient matched for each point cloud unit corresponds one-to-one with its own reflection deviation and reflection fluctuation value, avoiding correction deviation caused by coefficient mismatch; performing anti-interference enhancement operation on each unit of the target point cloud subset; and for each point cloud unit, extracting from the 3D point cloud dataset generated in step 104. The original reflection intensity (i.e., the signal amplitude value in the triplet of the unit, which directly reflects the reflected energy of the electromagnetic wave at the corresponding underground location) is extracted and multiplied with the matched dynamic correction coefficient. Through this operation, the reflection intensity of the units with more severe interference (because their reflection deviation and fluctuation values ​​are large, and the corresponding dynamic correction coefficients are also large) is specifically enhanced, or the interference components are effectively suppressed. For units with more stable original signals (smaller dynamic correction coefficients), their reflection intensity is kept within a reasonable range to avoid over-correction that destroys the original characteristics of the effective signal. The enhanced reflection intensity of each point cloud unit is obtained through the above calculation.

[0078] Step 404 above may include extracting the original spatial coordinates and echo delay of each point cloud unit: obtaining the X and Y plane coordinates of the unit from the initial spatial coordinate set generated in step 102 (recording the underground plane position), and obtaining the echo delay of the unit from the depth delay value converted in step 103 (recording underground depth information). These parameters together constitute the spatial position and depth attributes of the point cloud unit. Based on the point cloud unit, the enhanced reflection intensity obtained in step 403 is associated with its original spatial coordinates (X, Y) and echo delay to form a new triplet data containing enhanced reflection intensity-spatial coordinates-echo delay. Each new triplet completely records the enhanced reflection characteristics, plane position, and depth information of the corresponding underground position. The new triplet data of all point cloud units are summarized and integrated in the order of spatial coordinate distribution to form a structured dataset.

[0079] In this embodiment, dynamic correction coefficients are generated based on reflection deviation and fluctuation values ​​to achieve precise control where the stronger the interference, the greater the correction force. By using threshold screening to focus on high-deviation areas, correction resources are concentrated on the severely interfered target point cloud subset, avoiding excessive intervention in stable signals and effectively suppressing interference components such as sidelobes and noise. Through the multiplication operation of reflection intensity and correction coefficients, the discriminability of effective signals is enhanced in a targeted manner. The reflection intensity of the interfered units is more easily distinguished from background noise after enhancement, while the signals of stable units remain unchanged, thus improving the overall signal-to-noise ratio and feature clarity of the point cloud data. The enhanced reflection intensity is recombined with the original spatial coordinates and echo delay, optimizing reflection features while completely preserving spatial position and depth information.

[0080] In a preferred embodiment of the present invention, step 500 includes:

[0081] Step 501: Convert the corrected multi-dimensional point cloud feature set into a three-dimensional feature tensor and use it as network input; Step 502: Input the three-dimensional feature tensor into the feature extraction module of a pre-trained closed-loop convolutional neural network to extract multi-scale spatial features; Step 503: Based on the multi-scale spatial features, update the parameters of the fully connected layer through a transfer learning strategy to optimize the network's adaptability to the target region, thereby obtaining optimized multi-scale spatial features; Step 504: Based on the optimized multi-scale spatial features, obtain the confidence level of the ore body category and the probability value of the presence of an ore body at each spatial location; Step 505: Fuse the confidence level and the probability value of the presence of an ore body at each spatial location to generate a point cloud with a spatial location probability distribution.

[0082] In this embodiment of the invention, step 501 may include, based on the corrected multi-dimensional point cloud feature set generated in step 404, performing format conversion on the feature set to adapt to the neural network input requirements. The corrected multi-dimensional point cloud feature set includes the enhanced reflection intensity, original spatial coordinates (X, Y), and echo delay (depth information) of each point cloud unit. These data are spatially distributed to form a structured three-dimensional data set. Subsequently, with the X-dimensional, Y-dimensional, and depth-dimensional echo delay as the three spatial axes of the tensor, the enhanced reflection intensity corresponding to each spatial position is used as the feature value of the tensor and sequentially filled into the corresponding dimensions of the tensor according to the spatial position order to form a three-dimensional feature tensor with spatial dimension, depth dimension, and feature value structure.

[0083] Step 502 above may include the following: the three-dimensional feature tensor enters the shallow convolutional layer of the feature extraction module; the shallow convolutional layer is configured with a small convolutional kernel (such as a 3×3×3 three-dimensional convolutional kernel). This type of convolutional kernel has a narrow receptive field and can focus on local small-range feature changes in the tensor; during processing, the convolutional kernel slides across the three-dimensional feature tensor according to a preset stride, performs convolution operations with the feature values ​​in the tensor point by point (i.e., multiplies the corresponding elements and then sums them), and filters out invalid information and retains valid features through an activation function (such as the ReLU function); through this process, the shallow convolutional layer extracts local subtle spatial features from the three-dimensional feature tensor, such as a sudden increase in reflection intensity in a small area (which may correspond to local reflection at the edge of the ore body), subtle differences in reflection intensity between adjacent point cloud units, and other detailed features. These features can accurately characterize the local variation law of the ore body reflection signal; finally, the shallow convolutional layer outputs multiple shallow feature maps, each feature map corresponding to a distribution of a local subtle feature, providing basic detailed features for subsequent deep processing.

[0084] After the shallow convolutional layers output shallow feature maps, the features enter the deep convolutional and pooling layers of the module for further processing. The deep convolutional layers are equipped with larger convolutional kernels (such as 5×5×5 or 7×7×7 three-dimensional convolutional kernels), which significantly expand their receptive field and can cover a wider spatial area. When the convolutional kernel slides across the shallow feature map, it integrates the correlation between multiple local subtle features and captures the feature distribution patterns in a large spatial area. For example, the extension direction of continuous areas of enhanced reflectivity and the cluster of abnormally high reflectivity points in a large area of ​​low reflectivity. These features reflect the possible global macroscopic distribution trend of the ore body. Subsequently, the pooling layer downsamples the feature map after deep convolution (such as max pooling or average pooling). It slides across the feature map through a preset pooling window and selects the maximum or average value within the window as the output. After deep convolution and pooling processing, multiple deep feature maps are output.

[0085] While shallow convolutional layers output shallow feature maps and deep convolutional and pooling layers output deep feature maps, the module fuses these two types of features through a skip connection structure, ensuring that the correlation between multi-scale features is not lost. The skip connection directly passes the shallow feature map to the processing path of the deep feature map, fusing them through feature concatenation or feature overlay. For example, it adds the local detail features at corresponding positions in the shallow feature map to the global macro features at the same positions in the deep feature map element-wise, or concatenates them along the channel dimension to form a new feature map. This fusion method can compensate for the local detail information that may be lost during deep convolution, allowing the global macro features to carry detailed annotations of local subtle features, while also associating local subtle features with the global distribution trend. The feature extraction module integrates all fused feature maps after local detail extraction by shallow convolutional layers, global trend capture by deep convolutional and pooling layers, and feature fusion by skip connections. During the integration process, the network filters and weights features of different levels and types, retaining the most critical features for ore body identification (such as matching features between local high-intensity reflection details and global continuous reflection trends). Finally, it extracts multi-scale spatial features from the three-dimensional feature tensor, which include both local subtle spatial features and global macroscopic spatial features. These multi-scale spatial features cover both small-scale reflection change details in point cloud data and large-scale spatial distribution patterns.

[0086] Step 503 above may include, based on the multi-scale spatial features extracted in step 502, using a transfer learning strategy to optimize the network's adaptability to the target detection area; inputting the extracted multi-scale spatial features into the fully connected layer of the network, using the corrected point cloud feature data of the target detection area as training samples, calculating the error between the network output and the actual ore body labeling, adjusting the weights and biases of the fully connected layer through a backpropagation algorithm to make the fully connected layer more suitable for the reflection feature distribution pattern of the target area; after multiple rounds of iterative optimization, when the error drops below a preset threshold, stopping parameter updates to obtain the optimized multi-scale spatial features, which are more in line with the ore body reflection characteristics of the target area.

[0087] Step 504 above may include calculating the confidence score and spatial location probability value of the ore body category based on the multi-scale spatial features optimized in step 503. The confidence score of the ore body category is calculated using the Softmax classifier at the end of the network, reflecting the reliability of the current feature belonging to the ore body category. The formula is as follows: C = max(Softmax( f opt )) ore Where C represents the confidence level of the ore body category, with a value range of [0, 1]. The closer the value is to 1, the higher the confidence level that the area is an ore body. fopt The optimized multi-scale spatial features output in step 503 represent feature vectors adapted to the target region through transfer learning, containing both subtle local features (such as abrupt changes in reflection intensity) and macroscopic global features (such as continuous reflection trends); Softmax( f opt ) represents the Softmax activation function, used to transform a feature vector into a class probability distribution; max(Softmax( f opt )) ore This represents the probability value corresponding to the ore body category in the Softmax output (since the classifier outputs two probabilities, ore body and non-ore body, the confidence score uses the probability value of the ore body category as the final result); the spatial location probability value reflects the possibility of an ore body existing in a single point cloud unit in three-dimensional space, as shown in the formula: P(x, y, z) = Sigmoid(W reg ·[ f opt (x, y, z), s (x, y), r (x, y, z)] + b);

[0088] Where P(x, y, z) represents the probability of a mineral body existing at the spatial location (x, y, z), with a value ranging from [0, 1]. The closer the value is to 1, the greater the probability that the point cloud unit is a mineral body; (x, y) are the planar coordinates, and z is the depth corresponding to the echo delay; Sigmoid(·) represents the Sigmoid activation function, used to map the output of the regression module to the interval [0, 1]; W reg The weight coefficients are parameters obtained through transfer learning optimization, used to quantify the influence of multi-scale features, spatial coordinates, and reflection intensity on the probability value. f opt (x, y, z) represents the optimized multi-scale spatial feature components at position (x, y, z), including local reflection details (such as signal amplitude) and global correlation features (such as reflection continuity with adjacent cells) at that position; s(x, y) represents the spatial coordinate feature of position (x, y), which is usually a normalized value of planar coordinates (such as the offset relative to the starting scan position), used to characterize the planar position of the point in the detection area; r(x, y, z) represents the reflection intensity feature of position (x, y, z), that is, the enhanced reflection intensity value in step 403, which directly reflects the strength of electromagnetic wave reflection energy at that position; b represents the bias term of the spatial regression module, which is a parameter optimized through training, used to adjust the reference offset of the regression output and improve the accuracy of probability calculation.

[0089] Step 505 above may include: based on the detection accuracy requirements and historical data verification, a confidence threshold (e.g., 0.6) is preset to distinguish between high and low confidence category judgment results; for the ore body category confidence obtained in step 504, if its value is greater than the preset threshold (e.g., 0.7 > 0.6), it is judged as high confidence, and the probability values ​​of all spatial locations in that area are multiplied by an enhancement coefficient greater than 1 (e.g., 1.2); if the confidence is less than or equal to the preset threshold (e.g., 0.5 ≤ 0.6), it is judged as low confidence, and the probability values ​​of the corresponding spatial locations are multiplied by an attenuation coefficient less than 1 (e.g., 0.8) to reduce its weight and reduce the impact of low confidence judgment on positioning. To mitigate interference with the results, for each weighted spatial location probability value, find its corresponding original spatial coordinates (X, Y plane coordinates) and echo delay (depth information); combine these three into a triplet data set of plane coordinates (X, Y) - depth (echo delay) - adjusted probability value, ensuring that each probability value accurately corresponds to a specific underground spatial location; according to the spatial coordinate distribution pattern of the detection area, arrange all triplet data sets in the order of X and Y coordinates to form a three-dimensional dataset covering the entire detection area. This dataset is the spatial location probability distribution point cloud, where each point contains location information and the corresponding probability of ore body existence, and the high-probability dense area is the potential distribution range of the ore body.

[0090] In this embodiment, by using shallow and deep convolutions and skip connections, both local subtle features (such as strong edge reflections) and global macroscopic features (such as continuous distribution trends) of the ore body's reflection are captured simultaneously, avoiding the limitations of single-scale features and improving the comprehensiveness of feature representation. By freezing the feature extraction capabilities of the pre-trained network and focusing on updating the parameters of the fully connected layers, the model can quickly adapt to the geological differences (such as medium properties and interference types) of the target detection area, solving the problems of small sample size and strong regional specificity of ground-penetrating radar data, and significantly improving the model's generalization ability and recognition efficiency in complex environments. By using the Softmax classifier and the Sigmoid activation function, quantified ore body category confidence scores and spatial location probability values ​​are generated respectively, transforming abstract features into interpretable probability indicators, providing a definite credibility reference for the positioning results, and reducing subjective judgment errors. By fusing the confidence scores and spatial probability values ​​to generate a three-dimensional probability distribution point cloud, the accuracy and reliability of ore body positioning results can be improved.

[0091] In a preferred embodiment of the present invention, step 600 includes:

[0092] Step 601: Construct a three-dimensional spatial grid coordinate system based on the spatial location probability distribution point cloud; Step 602: Calculate the theoretical two-way travel time from each grid node to the radar object in the three-dimensional spatial grid coordinate system, combining the electromagnetic wave two-way travel time calculation rules; Step 603: Integrate the Carnia resistivity measurement results and map the Carnia resistivity data to the corresponding nodes in the three-dimensional spatial grid coordinate system; convert the resistivity value of each node into a dielectric constant value according to the resistivity-dielectric constant mapping relationship; convert the dielectric constant value into the electromagnetic wave propagation speed value of each grid node based on the electromagnetic wave propagation speed calculation rules; Step 604: Based on the theoretical two-way travel time and electromagnetic wave propagation speed value, construct a spherical projection equation system with the radar antenna position as the origin; substitute the theoretical two-way travel time and electromagnetic wave propagation speed value into the spherical projection equation system; solve the spherical projection equation system to derive the three-dimensional spatial coordinates of the reflected signal source; aggregate the spatial coordinates of all reflected signal sources to generate a point cloud for the ore body depth and horizontal positioning.

[0093] In this embodiment of the invention, step 601 may include, based on the spatial location probability distribution point cloud output in step 505, first determining the three-dimensional spatial boundary of the target detection area: by extracting the X coordinate (horizontal horizontal direction), Y coordinate (horizontal vertical direction), and depth value (Z direction) corresponding to the echo delay of all points in the probability distribution point cloud, taking their maximum and minimum values ​​respectively, and delineating the X-axis range [X] in the three-dimensional space. min X max ], Y-axis range [Y min Y max ] and Z-axis range [Z min Z max Based on the preset detection accuracy requirements (e.g., lateral positioning error ≤ 0.5m, depth positioning error ≤ 0.2m), the grid resolution parameters △X, △Y, and △Z are set (where △X is the grid spacing in the X direction, △Y is the grid spacing in the Y direction, and △Z is the grid spacing in the Z direction). Based on the three-dimensional spatial boundary and resolution, the target area is divided into several regular cubic grids. The vertex of each grid is defined as a grid node, and each node is assigned a unique three-dimensional coordinate (X, Y, △Z). k Y k Z k ), where k is the node index, X k ∈[X min X max ], Y k ∈[Y min Y max ], Z k ∈[Z min Z max This allows for the construction of a three-dimensional spatial grid coordinate system covering the entire detection area, enabling the discretization of underground space.

[0094] Step 602 above may include the total time from when the electromagnetic wave is emitted from the ground-penetrating radar antenna, to when it propagates to the underground target point (grid node), is reflected, and then returns to the antenna for reception; first, the surface position coordinates of the radar antenna are determined as (X0, Y0, 0), where Z = 0 represents the surface reference plane; for any node (X0, Y0, 0) in the three-dimensional spatial grid coordinate system... k Y k Z k ), calculate the spatial straight-line distance D between the node and the antenna. k Based on spatial geometric relationships, D k The horizontal distance is the composite distance of the horizontal distance and the vertical depth. The horizontal distance and vertical depth are calculated as Z. k Using the initial electromagnetic wave propagation velocity v0 (based on the average velocity of the medium obtained from previous geological surveys), the one-way propagation time t = D at this node is calculated. k / v0; thus, the theoretical two-way travel time t is obtained. d = 2 × t, that is, t d =2D k / v0, traverse all grid nodes in this way to obtain the theoretical two-way travel time corresponding to each node, which will be used as the time constraint parameter for subsequent inversion.

[0095] Step 603 above may include collecting measured resistivity data of the target detection area, which records the planar coordinates (such as the horizontal X coordinate and the vertical Y coordinate) and the corresponding resistivity value (in ohms-meters) of each measurement point; matching the measured data with the three-dimensional spatial grid coordinate system: for each measurement point, based on its X and Y coordinates, finding the grid node with the closest planar position in the constructed three-dimensional grid (for example, if a measurement point has X = 2.3 meters and Y = 5.7 meters, and the grid node X interval is 0.5 meters, then a node with X = 2.5 meters and Y = 5.5 meters is matched), and directly assigning the resistivity value of the measurement point to this matched grid node; For nodes in the 3D mesh that are not covered by measured data (i.e., nodes for which no corresponding measured point was found), the resistivity value is estimated using inverse distance weighted interpolation: Taking the node as the center, at least three mesh nodes with known resistivity values ​​are selected within a certain range; the 3D spatial distance from the node to each known node is calculated (considering both horizontal and depth distances); the closer the known node, the greater its influence on the target node (closer nodes have higher weights, farther nodes have lower weights); the resistivity values ​​of each known node are summed according to their weights to obtain the estimated resistivity value of the target node, ensuring that all mesh nodes have corresponding resistivity attributes.

[0096] Through preliminary core experiments on similar geological areas (such as rock strata with the same lithology as the exploration area), a large amount of experimental data was obtained. The resistivity and dielectric constant of different core samples were measured, and a resistivity-dielectric constant correspondence table was established (for example, when the resistivity is 100 ohm-meters, the dielectric constant is 8; when the resistivity is 200 ohm-meters, the dielectric constant is 6, etc.). For each node in the three-dimensional grid, the corresponding dielectric constant was looked up in the above experimental correspondence table based on its determined resistivity value. If the node resistivity value is exactly in the table, the corresponding dielectric constant is directly used. If the resistivity value is between two values ​​in the table (for example, if the node resistivity is 150 ohm-meters, 100 corresponds to 8 and 200 corresponds to 6 in the table), linear interpolation is used to calculate (150 is between 100 and 200, so the dielectric constant is 7). If the resistivity value exceeds the range in the table, the trend in the table is used for extrapolation (for example, when the resistivity is greater than 200 ohm-meters, the value is estimated according to the value change pattern in the table). Finally, the corresponding dielectric constant value is determined for each node.

[0097] The propagation speed of electromagnetic waves in a medium is inversely proportional to the square root of the dielectric constant of the medium, and the speed of light in a vacuum is a known fixed value (approximately 300,000 km / s). For each grid node, the square root is first calculated based on its dielectric constant value (for example, when the dielectric constant is 9, the square root is 3). Then, the speed of light in a vacuum is divided by this square root (300,000 km / s ÷ 3 = 100,000 km / s) to obtain the actual propagation speed of electromagnetic waves at that node. This node-specific propagation speed replaces the regional average speed used in step 602, so that the propagation speed of each grid node can accurately reflect the medium characteristics at its location, achieving a fine distribution of propagation speed in three-dimensional space.

[0098] Step 604 above may include taking the actual position of the ground-penetrating radar antenna on the ground surface as the origin (the coordinates of this position have been predetermined and are denoted as the ground surface reference point), and treating each node in the three-dimensional spatial grid as a potential signal source that may reflect electromagnetic waves (i.e., the location of a possible ore body reflection interface); according to the basic law of electromagnetic wave propagation, that is, the total propagation distance of the electromagnetic wave from the antenna to the node and then back to the antenna is equal to half of the product of the electromagnetic wave propagation speed at that node and the round-trip travel time (because the round-trip travel time is the total time for the electromagnetic wave to travel to and from the node, and the one-way distance needs to be divided by 2).

[0099] For each grid node, measure its three-dimensional straight-line distance to the radar antenna. The calculation must consider both horizontal and depth positional differences: first, calculate the distance between the node and the antenna on the horizontal plane (using the difference between their X and Y coordinates), then combine this with the node's depth value (Z coordinate) to calculate the straight-line distance from the antenna to the node (i.e., the actual one-way distance of electromagnetic wave propagation) through spatial geometry. Based on the node's electromagnetic wave propagation speed (the node-specific speed determined in step 603) and theoretical two-way travel time (the node's two-way travel time calculated in step 602), calculate the theoretical one-way propagation distance: multiply the propagation speed by the two-way travel time, then divide the result by 2 (because the two-way travel time includes the round-trip time of the electromagnetic wave, the one-way distance needs to be half). Correlate the node's actual spatial distance with the theoretical propagation distance to form a correspondence between actual distance and theoretical propagation distance.

[0100] The above operations are performed on each grid node one by one. The relationship between the actual distance and the theoretical propagation distance of each node constitutes an independent equation. The equations of all nodes are combined to form a spherical projection equation system. Each equation corresponds to a sphere with the antenna as the center and the theoretical propagation distance as the radius. If a node is located on this sphere, it means that its position conforms to the electromagnetic wave propagation law and it may be an effective reflection source. If it deviates from the sphere, it does not conform to the law. In this way, the spatial position of each node is bound to the electromagnetic wave propagation characteristics, providing a quantitative basis for subsequent screening of effective reflection sources. The least squares iterative method is used to solve the spherical projection equation system: For each grid node, the difference between its actual spatial distance and the propagation speed × two-way travel time ÷ 2 (i.e., the residual) is calculated. The smaller the residual, the more the position of the node conforms to the electromagnetic wave propagation law. A residual threshold is set. (e.g., 0.1 meters, determined according to detection accuracy requirements), compare the residuals of each node one by one: if the node residual is less than or equal to the threshold, the node is determined to be a valid reflection signal source, that is, the location of the node is very likely to be the ore body reflection interface; if the residual is greater than the threshold, it means that the node does not conform to the electromagnetic wave propagation law, and is determined to be an invalid node and excluded; through this screening, valid reflection signal sources that meet the conditions are accurately extracted from all grid nodes. These nodes constitute the core data for ore body positioning; the X and Y coordinates of each node reflect its horizontal planar position, and the Z coordinate reflects its underground depth position; these coordinates are organized and integrated in spatial distribution order to form a complete three-dimensional dataset, which is the ore body depth and horizontal positioning point cloud, in which each point accurately corresponds to the location of the ore body that may exist underground.

[0101] This invention integrates spatial location probability distribution point clouds, electromagnetic wave two-way travel time, and Cania resistivity data to construct multi-dimensional constraints. This overcomes the limitation of single radar signals being easily interfered with in complex media, enabling ore body positioning to simultaneously consider signal probability characteristics, propagation time patterns, and differences in medium electrical properties, significantly reducing positioning errors. By constructing a three-dimensional spatial grid coordinate system, the underground space is discretized into regular nodes. Combining inverse distance weighted interpolation and physical quantity conversion, each node is assigned precise resistivity, dielectric constant, and propagation velocity attributes, achieving spatial fine-grained allocation of electromagnetic wave propagation parameters and avoiding overall errors caused by regional average parameters. Based on the spherical projection equations and residual screening mechanism, effective reflected signal sources conforming to the electromagnetic wave propagation laws are accurately extracted from massive grid nodes. Invalid interference nodes are eliminated by quantizing residual thresholds, ensuring that the positioning point cloud focuses on the real ore body reflection interface and improving the reliability of the results.

[0102] In a preferred embodiment of the present invention, step 700 includes:

[0103] Step 701: Based on the inverse projection positioning point cloud, establish a two-dimensional grid coordinate system on the horizontal plane; Step 702: Use the point cloud depth value in each grid cell as elevation data to generate a two-dimensional depth distribution matrix; Step 703: Perform multi-scale two-dimensional empirical mode decomposition on the two-dimensional depth distribution matrix to obtain intrinsic mode function components at different scales; Step 704: Filter out high-frequency noise components in the intrinsic mode function components according to a preset noise frequency threshold to obtain the retained intrinsic mode function components; Step 705: Reconstruct the target depth distribution matrix using the retained intrinsic mode function components; Step 706: Map the grid cells in the reconstructed target depth distribution matrix with elevation values ​​greater than the confidence threshold to high-confidence target point clouds.

[0104] In this embodiment of the invention, step 701 may include: based on the reverse projection positioning point cloud generated in step 604, extracting the horizontal coordinates (X coordinates and Y coordinates) of all point clouds, and determining the spatial range within the horizontal plane: defining the X-axis range by statistically analyzing the minimum and maximum values ​​of the X coordinates, and defining the Y-axis range by statistically analyzing the minimum and maximum values ​​of the Y coordinates; setting the resolution of the two-dimensional grid (e.g., the grid interval in both the X and Y directions is 0.3 meters) according to the density of the target point cloud and the noise filtering accuracy requirements, and uniformly dividing the defined horizontal range into several square grid units according to this resolution, with each grid unit uniquely identified by the X and Y coordinates of its upper left corner or center, thereby constructing a two-dimensional grid coordinate system covering the entire horizontal detection area.

[0105] Step 702 above may include: traversing all points in the reverse projection positioning point cloud, assigning each point to the corresponding grid cell in the two-dimensional grid coordinate system constructed in step 701 according to its horizontal coordinates (X, Y) (i.e. finding the grid cell to which the X and Y coordinates of the point belong); for each grid cell, collecting the depth values ​​(Z values) of all points falling into it, and obtaining the representative depth value of the grid cell by calculating the average value (or taking the median to avoid the influence of extreme values), which is used as the elevation data of the cell; according to the arrangement order of the grid cells in the two-dimensional coordinate system (such as from left to right, from top to bottom), filling the elevation data of all grid cells into a matrix in sequence to form a two-dimensional depth distribution matrix. The rows of the matrix correspond to the grid number in the Y direction, and the columns correspond to the grid number in the X direction. Each element value in the matrix is ​​the depth value of the corresponding grid cell, which intuitively presents the distribution of the depth of the underground ore body on the horizontal plane.

[0106] Step 703 above may include traversing each grid cell in the two-dimensional depth distribution matrix and identifying local extrema: for each cell, comparing its depth value with that of its surrounding neighboring cells (e.g., cells in the top, bottom, left, right, and diagonal directions, a total of 8 directions); if the depth value of the cell is greater than that of all neighboring cells, it is marked as a local maximum point; if it is less than that of all neighboring cells, it is marked as a local minimum point; collecting all local maximum points and connecting these points using an interpolation algorithm (e.g., cubic spline interpolation) to form an upper envelope covering the entire matrix range (reflecting the local maximum trend of depth values); similarly, collecting all local minimum points and connecting them using interpolation... The lower envelope (reflecting the local minimum trend of depth values) is formed. The upper and lower envelopes must completely cover the edge area of ​​the matrix to ensure no data gaps. The average values ​​of the upper and lower envelopes at each grid cell are calculated: for each cell in the matrix, the arithmetic mean of the upper and lower envelope values ​​at that location is taken to obtain a smooth trend surface (reflecting the overall trend of change at that scale). The original depth value of each grid cell is subtracted from the trend surface value at the corresponding location to obtain a new set of values, forming the first component. This component mainly contains small-scale details that fluctuate rapidly in the original data (such as random noise or local interference), i.e., the high-frequency component.

[0107] The residual matrix is ​​formed by subtracting the first high-frequency component from the original matrix (this matrix has removed small-scale noise and retains smoother features). The above steps are repeated on the residual matrix to identify local maxima and minima again, reconstruct the upper and lower envelopes, and calculate a new trend surface. The new trend surface is then subtracted from the residual matrix to obtain the second component. This component has a lower fluctuation frequency than the first high-frequency component, corresponding to medium-scale features (such as local boundaries of the ore body or continuously distributed details), i.e., the mid-frequency component. The decomposition is repeated on the new residual matrix (the value after subtracting the mid-frequency component from the previous residual): Each round... During the iteration, the fluctuations of the envelope become smoother and smoother, and the frequencies of the separated components gradually decrease. After multiple iterations, when there are no more obvious local extrema in the residual matrix (or the fluctuation amplitude is less than the preset threshold), the decomposition stops. At this point, the final components are low-frequency components, which reflect the overall distribution trend of the ore body depth in the entire exploration area (such as the large-scale depth variation pattern). They are arranged from high frequency to low frequency. High-frequency components correspond to small-scale noise or subtle interference, mid-frequency components correspond to the local features of the ore body, and low-frequency components correspond to the overall distribution trend. These components completely cover all scale features in the original depth distribution matrix.

[0108] Step 704 above may include setting a noise frequency threshold by combining previous exploration data in similar geological areas or by analyzing the noise characteristics of the current exploration area (such as the typical fluctuation frequency of interference signals); for example, based on the highest fluctuation frequency of noise signals in historical data, the threshold may be set as a critical value that can distinguish between rapid, irregular fluctuations and slow, regular fluctuations; examining each intrinsic mode function component obtained from step 703 one by one, i.e., observing high-frequency components, whose values ​​change drastically between adjacent grid cells, exhibiting rapid, irregular fluctuations (such as sudden increases or decreases in depth values ​​of adjacent cells without a continuous trend), which is consistent with the characteristics of noise signals; observing mid- and low-frequency components, whose values ​​change gently, and whose depth values ​​of adjacent grid cells show a continuous, gradual trend (such as gradually deepening along a certain direction or remaining stable), which is consistent with the regular distribution characteristics of ore body reflection signals; screening and eliminating noise components, comparing the fluctuation frequency of each component with the preset threshold, and if the fluctuation frequency of a component is higher than the threshold (i.e., it belongs to rapid, irregular fluctuations), it is determined to be a noise component and eliminated from the component set; retaining all mid- and low-frequency components with fluctuation frequencies lower than the threshold, as these components contain the effective characteristics of the ore body.

[0109] Step 705 above may include collecting all the low- and medium-frequency intrinsic mode function components selected in step 704, ensuring that these components cover the medium-scale characteristics (such as local boundaries) and overall distribution trends (such as large-scale depth variations) of the ore body; for each retained low- and medium-frequency component, numerical superposition is performed according to the grid cell position; for each cell in the two-dimensional grid, the corresponding value of the cell in each retained component is found, and these values ​​are added one by one (for example, if a grid cell has a value of 5 in the first mid-frequency component, a value of 3 in the second mid-frequency component, and a value of 10 in the low-frequency component, then the superimposed value is 5 + 3 + 10 = 18). The superposition results of all grid cells are arranged according to their position (row and column order) in the two-dimensional coordinate system to form a new two-dimensional depth distribution matrix.

[0110] Step 706 above, extracting high-confidence target point cloud, may include: according to the reliability requirements of ore body detection, setting a preset elevation confidence threshold (such as setting a minimum depth threshold or a difference threshold with the background value based on the typical range of ore body depth in historical data); for each grid cell, if its elevation value (depth value) is greater than the preset confidence threshold, it is determined that the cell contains a high-confidence ore body signal; extracting the horizontal coordinates (X, Y) of these cells and their corresponding depth values ​​(Z), combining them into three-dimensional point cloud data, and finally generating a high-confidence target point cloud that only contains points that meet the depth characteristics and whose noise has been filtered out, accurately focusing on the core distribution area of ​​the ore body.

[0111] In a preferred embodiment of the present invention, step 800 includes:

[0112] Step 801: Based on the high-confidence target point cloud, extract the peak reflection intensity points as candidate point sets for ore body vertices; Step 802: Fit the candidate point sets for ore body vertices to a hyperbolic equation and calculate the fitting residual; Step 803: Based on the fitting residual, select candidate points whose residuals are less than a preset error threshold as valid ore body vertices; Step 804: Based on the calculation rules for the propagation speed of electromagnetic waves in underground media, correct the depth coordinates of the valid ore body vertices; Step 805: Using the corrected ore body vertices as control points, generate a three-dimensional spatial surface of the ore body through radial basis function interpolation; Step 806: Uniformly sample on the three-dimensional spatial surface to generate a three-dimensional spatial morphology distribution point cloud of the ore body.

[0113] In this embodiment of the invention, step 801, which extracts the candidate point set of the ore body vertex, may include: traversing the high-confidence target point cloud generated in step 706, retrieving the original reflection intensity data corresponding to each point (this data comes from the reflection energy value recorded by the radar signal receiver); selecting the top 20% of points with reflection intensity values ​​that are significantly higher than the surrounding points (e.g., the intensity of a certain point is more than 30% higher than the average intensity of all points within a 5-meter radius); and collecting the three-dimensional coordinates (X, Y horizontal coordinates and Z depth coordinates) of the points that meet the conditions to form the candidate point set of the ore body vertex. These points are considered as potential marker points of the top boundary of the ore body due to their prominent reflection energy.

[0114] Step 802 above, polynomial curve fitting and residual calculation may include: visualizing the distribution of point cloud on the horizontal plane to observe whether there is an obvious bending trend (such as a single-peak curve, gentle undulations or multi-segment gradual change characteristics); if the candidate points show a simple gradual change trend (such as an approximate straight line or slight bending), a second-order polynomial is selected; if there are obvious undulations or inflection points, a third-order polynomial is selected (to avoid overfitting of the curve due to excessively high order, resulting in meaningless fluctuations). At the same time, the fitting principal axis is determined according to the extension direction of the ore body (such as X as the independent variable if it extends along the X-axis, and Y as the independent variable if it extends along the Y-axis).

[0115] Based on the horizontal coordinates of the candidate points, a polynomial fitting model is established; for example, if fitting along the X-axis, the model form is Y coordinate = aX 3 +bX 2 +cX+d (third-order polynomial), where a, b, c, and d are coefficients to be determined; iterate through all candidate points, substitute the X coordinate of each point into the model, and calculate the deviation between the model's predicted Y coordinate and the actual Y coordinate of that point (i.e., predicted value - actual value); by adjusting the values ​​of a, b, c, and d, minimize the sum of squared deviations of all points (the smaller the sum of squared deviations, the better the model fits the actual distribution). After multiple rounds of iterative adjustments, finally determine a set of optimal coefficients to obtain a polynomial fitting curve that fits the distribution trend of the candidate points.

[0116] For each candidate point, the horizontal distance from it to the fitted curve is calculated as the residual. For a curve fitted along the X-axis, the predicted Y-coordinate corresponding to the X-coordinate of the point on the fitted curve is first found. The difference between the actual Y-coordinate and the predicted Y-coordinate of the point (i.e., the horizontal deviation distance) is calculated. The shortest horizontal distance from the point to the fitted curve is calculated using the vertical distance formula (to ensure that the residual reflects the true degree of deviation). The residual result directly reflects the consistency between the candidate point and the overall distribution trend. The smaller the value, the more the point conforms to the natural distribution law of the ore body apex.

[0117] Step 803 above, screening effective orebody vertices may include: based on the accuracy requirements of on-site detection, setting a residual error threshold (e.g., referring to the experience value of similar projects, setting the threshold to 0.4 meters); judging the residuals of candidate points one by one; if the residual is less than or equal to the threshold, it means that the distribution of the point conforms to the overall trend of the orebody vertices and is judged as an effective orebody vertices; if the residual exceeds the threshold, it is regarded as an abnormal point affected by local interference and is excluded; summarizing the three-dimensional coordinates of all effective vertices to form an effective vertex set, which can accurately reflect the main distribution characteristics of the top of the orebody.

[0118] Step 804 above, correcting the depth coordinates of valid vertices, may include: for each valid orebody vertex, obtaining the electromagnetic wave propagation velocity of its corresponding grid node (this velocity comes from the refined velocity value based on resistivity conversion in step 603, taking into account the propagation differences of different media); combining this with the two-way travel time of the radar signal corresponding to that vertex (the total time from signal transmission to reception), recalculating the depth coordinates: multiplying the propagation velocity by the two-way travel time, and then dividing by 2 (i.e., the one-way propagation distance) to obtain the corrected depth value. This method eliminates the error caused by velocity approximation in the initial depth calculation, making the depth coordinates closer to the actual burial conditions.

[0119] Step 805 above, generating the three-dimensional spatial surface of the ore body, may include: using the corrected effective ore body vertices as key control points, constructing the three-dimensional spatial surface using the inverse distance weighted interpolation method; first determining the spatial range of the interpolation to ensure coverage of all effective vertices and the area where the ore body may extend; calculating the three-dimensional distance from each blank position in the area to be interpolated to each control point, with the closer the control point, the greater its influence weight on that position; and weighting the depth values ​​of the control points according to the weight ratio to obtain the depth values ​​of the blank positions, gradually filling the entire area to form a continuous three-dimensional surface that can fully present the undulating shape of the top of the ore body.

[0120] Step 806 above, generating a three-dimensional spatial morphology distribution point cloud of the ore body, may include performing regular sampling on the generated three-dimensional spatial surface: setting the sampling interval according to the required level of detail (e.g., every 0.3 meters in the horizontal direction, and naturally undulating with the surface in the depth direction), uniformly selecting sampling points within the X and Y horizontal range of the surface; recording the three-dimensional coordinates (X and Y horizontal positions and corresponding Z depth) of each sampling point, and integrating these coordinates into a three-dimensional dataset, i.e., the three-dimensional spatial morphology distribution point cloud of the ore body, which clearly shows the overall morphology, horizontal distribution range, and depth variation of the ore body underground.

[0121] This embodiment fully considers the impact of the heterogeneity of the underground medium on the propagation speed, corrects the depth coordinates of the effective vertices, avoids depth deviation caused by velocity approximation, and significantly improves the calculation accuracy of the ore body burial depth. Using the corrected vertices as control points, a continuous and smooth three-dimensional spatial surface is generated through radial basis function interpolation, which can effectively fill the blank areas between vertices and completely restore the undulating shape and spatial distribution trend of the top of the ore body. The morphological distribution point cloud generated by uniform sampling on the three-dimensional surface not only retains the key morphological features of the ore body, but also ensures the data density balance through regular sampling.

[0122] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for locating and identifying underground ore bodies based on ground-penetrating radar, characterized in that, The method includes: The original radar signal is preprocessed to generate a three-dimensional point cloud dataset containing spatial coordinates and reflection intensity. Each point cloud data unit contains a triplet of position coordinates, echo delay, and signal amplitude. Based on the 3D point cloud dataset, the goodness of fit of the reflection intensity sequence and the side lobe interference features of each point cloud unit are extracted, and the reflection deviation is calculated. Based on a 3D point cloud dataset, the spatial distribution of binary pairs consisting of the maximum value and orientation angle of the reflected signal of adjacent point cloud units is extracted, and the reflected wave value is calculated. By fusing reflection deviation and reflection fluctuation values, dynamic correction coefficients are generated to enhance the anti-interference effect on the reflection intensity of a subset of the target point cloud, resulting in a corrected multi-dimensional point cloud feature set. The corrected multi-dimensional point cloud feature set is input into a pre-trained closed-loop convolutional neural network. The network parameters are optimized through transfer learning, and the point cloud of the ore body target's category confidence and spatial location probability distribution is output. Based on the spatial location probability distribution point cloud, combined with the electromagnetic wave two-way travel time calculation rules and the Carnia resistivity measurement results, the ore body depth and horizontal positioning point cloud are generated through the inverse projection algorithm. Multi-scale two-dimensional empirical mode decomposition is performed on the inverse projection positioning point cloud to filter out noise components and extract high-confidence target point cloud; Based on the high-confidence target point cloud, the coordinates of the ore body vertex are corrected by the hyperbola fitting algorithm, and combined with the calculation rules of the propagation speed of electromagnetic waves in the underground medium, a three-dimensional spatial distribution point cloud of the ore body is generated.

2. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 1, characterized in that, The raw radar signal is preprocessed to generate a 3D point cloud dataset containing spatial coordinates and reflection intensity. Each point cloud data unit contains a triplet of position coordinates, echo delay, and signal amplitude, including: The original radar signals collected from the polymetallic mineral cluster area are converted into time and frequency to generate a two-dimensional echo signal sequence image. The two-dimensional echo signal sequence image contains the timing and amplitude information of the reflected signals caused by fault cutting and quartz vein interpenetration. The temporal information of each pixel in the two-dimensional echo signal sequence image is mapped to the spatial position coordinates on the survey line, while the echo time axis information of the corresponding pixel is retained, generating a two-dimensional data field with spatial coordinates and original time delay. Based on the propagation speed of electromagnetic waves in the complex medium of the mining area, the echo time axis information in the two-dimensional data field with spatial coordinates and original time delay is converted into depth time delay value. Extract the signal amplitude corresponding to each spatial location coordinate point, bind it with the depth delay value and spatial coordinates of that point as a triple, and generate an anti-interference 3D point cloud dataset.

3. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 2, characterized in that, Based on a 3D point cloud dataset, the goodness of fit of the reflection intensity sequence and sidelobe interference features of each point cloud unit are extracted, and the reflection deviation is calculated, including: Extract the reflection intensity sequence of each point cloud unit from the three-dimensional point cloud dataset; Based on the reflection intensity sequence, an ideal distribution curve is fitted, and the variance of the fitting residual is calculated as a goodness-of-fit feature. Based on the same reflection intensity sequence, locate the main lobe peak point of the reflection intensity sequence, and search for the first trough point on both sides of the main lobe peak point as the center as the side lobe boundary; calculate the amplitude ratio based on the main lobe peak point and the highest point within the side lobe boundary as the side lobe peak attenuation rate; within the side lobe boundary range, measure the width of the continuous region where the amplitude value is 10% higher than the main lobe peak value as the side lobe width feature. The reflection deviation is calculated by weighted summation by integrating the goodness-of-fit feature with the sidelobe peak attenuation rate and sidelobe width feature.

4. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 3, characterized in that, Based on a 3D point cloud dataset, the spatial distribution of binary pairs consisting of the maximum value and orientation angle of the reflected signals from adjacent point cloud units is extracted, and the reflection fluctuation value is calculated, including: Extract the maximum value of the reflection intensity sequence for each point cloud unit from the three-dimensional point cloud dataset; Based on the spatial coordinates of point cloud units, calculate the direction vector of each point cloud unit and its neighboring units, and extract the direction angle; The maximum value of the reflection intensity sequence of each point cloud unit is combined with the orientation angle to generate a binary spatial distribution dataset; Based on the aforementioned binary spatial distribution dataset, the product of the standard deviation of the maximum reflection intensity between adjacent units and the rate of change of the direction angle is calculated as the reflection fluctuation value.

5. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 4, characterized in that, By fusing reflection deviation and reflection fluctuation values, dynamic correction coefficients are generated to enhance the anti-interference effect of the reflection intensity of a subset of the target point cloud, resulting in a corrected multi-dimensional point cloud feature set, including: Based on reflection deviation and reflection fluctuation values, dynamic correction coefficients are generated through nonlinear function mapping. Based on the gradient distribution of the reflection deviation, point cloud units with reflection deviations greater than a preset threshold are selected to form a subset of the target point cloud; For each cell in the target point cloud subset, its reflection intensity is multiplied by the corresponding dynamic correction coefficient to obtain the enhanced reflection intensity; The enhanced reflection intensity is recombined with the original spatial coordinates and echo delay to generate a corrected multi-dimensional point cloud feature set.

6. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 5, characterized in that, The corrected multi-dimensional point cloud feature set is input into a pre-trained closed-loop convolutional neural network. The network parameters are optimized through transfer learning, and the output is a point cloud containing the category confidence score and spatial location probability distribution of the ore body target, including: The corrected multi-dimensional point cloud feature set is converted into a three-dimensional feature tensor and used as network input. The three-dimensional feature tensor is input into the feature extraction module of a pre-trained closed-loop convolutional neural network to extract multi-scale spatial features; Based on the multi-scale spatial features, the parameters of the fully connected layer are updated through a transfer learning strategy to optimize the network's adaptability to the target region, thereby obtaining the optimized multi-scale spatial features. Based on the optimized multi-scale spatial features, the confidence level of the ore body category and the probability value of the existence of an ore body at each spatial location are obtained; By combining the confidence level and the probability value of the presence of a ore body at each spatial location, a point cloud of spatial location probability distribution is generated.

7. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 6, characterized in that, Based on the spatial location probability distribution point cloud, combined with the electromagnetic wave two-way travel time calculation rules and Carnia resistivity measurement results, a point cloud for ore body depth and horizontal positioning is generated through an inverse projection algorithm, including: A three-dimensional spatial grid coordinate system is constructed based on the point cloud of spatial location probability distribution; Combining the electromagnetic wave two-way travel time calculation rules, the theoretical two-way travel time from each grid node to the radar object is calculated in a three-dimensional spatial grid coordinate system; By integrating the Carnia resistivity measurement results, the Carnia resistivity data are mapped to the corresponding nodes in the three-dimensional spatial grid coordinate system; based on the resistivity-dielectric constant mapping relationship, the resistivity values ​​of each node are converted into dielectric constant values; based on the calculation rules of electromagnetic wave propagation speed in the medium, the dielectric constant values ​​are converted into electromagnetic wave propagation speed values ​​of each grid node. Based on the theoretical two-way travel time and electromagnetic wave propagation speed, a spherical projection equation system is constructed with the radar antenna position as the origin. The theoretical two-way travel time and electromagnetic wave propagation speed are substituted into the spherical projection equation system. The spherical projection equation system is solved to derive the three-dimensional spatial coordinates of the reflected signal source. The spatial coordinates of all reflected signal sources are aggregated to generate a point cloud of the ore body depth and horizontal positioning.

8. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 7, characterized in that, Multi-scale two-dimensional empirical mode decomposition is performed on the inverse projection localization point cloud to filter out noise components and extract high-confidence target point cloud, including: Based on the reverse projection positioning point cloud, a two-dimensional grid coordinate system is established on the horizontal plane; The point cloud depth value within each grid cell is used as elevation data to generate a two-dimensional depth distribution matrix. The two-dimensional depth distribution matrix is ​​subjected to multi-scale two-dimensional empirical mode decomposition to obtain intrinsic mode function components at different scales; Based on a preset noise frequency threshold, high-frequency noise components in the intrinsic mode function components are filtered out to obtain the retained intrinsic mode function components. The target depth distribution matrix is ​​reconstructed using the retained intrinsic mode function components; In the reconstructed target depth distribution matrix, grid cells with elevation values ​​greater than the confidence threshold are mapped to high-confidence target point clouds.

9. The method for locating and identifying underground ore bodies based on ground-penetrating radar according to claim 8, characterized in that, Based on a high-confidence target point cloud, the coordinates of the ore body vertices are corrected using a hyperbola fitting algorithm. Combined with the rules for calculating the propagation speed of electromagnetic waves in underground media, a three-dimensional spatial distribution point cloud of the ore body is generated, including: Based on the high-confidence target point cloud, the peak points of reflection intensity are extracted as the candidate point set of ore body vertices; Fit a hyperbolic equation to the candidate vertex point set of the ore body and calculate the fitting residual; Candidate points with residuals less than a preset error threshold are selected based on the fitting residuals and used as valid orebody vertices. Based on the calculation rules for the propagation speed of electromagnetic waves in underground media, the depth coordinates of the effective ore body apex are corrected. Using the corrected orebody vertices as control points, a three-dimensional spatial surface of the orebody is generated by radial basis function interpolation. The distribution point cloud of the three-dimensional spatial morphology of the ore body is generated by uniformly sampling on a three-dimensional curved surface.

10. A ground-penetrating radar-based underground ore body location and identification system, wherein the system implements the method as described in any one of claims 1 to 9, characterized in that, include: The preprocessing module is used to preprocess the raw radar signal to generate a three-dimensional point cloud dataset containing spatial coordinates and reflection intensity. Each point cloud data unit contains a triplet of position coordinates, echo delay, and signal amplitude. The calculation module is used to extract the goodness of fit of the reflection intensity sequence and the sidelobe interference features of each point cloud unit based on the 3D point cloud dataset, and to calculate the reflection deviation. The extraction module is used to extract the spatial distribution of the binary tuple consisting of the maximum value and direction angle of the reflection signal of adjacent point cloud units based on the 3D point cloud dataset, and to calculate the reflection fluctuation value. The fusion module is used to fuse reflection deviation and reflection fluctuation values, generate dynamic correction coefficients, enhance the anti-interference of the reflection intensity of the target point cloud subset, and obtain a corrected multi-dimensional point cloud feature set. The prediction module is used to input the corrected multi-dimensional point cloud feature set into a pre-trained closed-loop convolutional neural network, optimize the network parameters through transfer learning, and output the category confidence and spatial location probability distribution point cloud of the ore body target. The processing module is used to generate ore body depth and horizontal positioning point clouds based on spatial location probability distribution point clouds, combined with electromagnetic wave two-way travel time calculation rules and Carnia resistivity measurement results, through the inverse projection algorithm; and to perform multi-scale two-dimensional empirical mode decomposition on the inverse projection positioning point clouds to filter out noise components and extract high-confidence target point clouds. Based on the high-confidence target point cloud, the coordinates of the ore body vertex are corrected by the hyperbola fitting algorithm, and combined with the calculation rules of the propagation speed of electromagnetic waves in the underground medium, a three-dimensional spatial distribution point cloud of the ore body is generated.

Images