Method, device and storage medium for identifying underground anomaly interface based on ground penetrating radar

Abstract

The present invention discloses a method, device, and storage medium for identifying the interface of an underground anomaly based on ground-penetrating radar. The identification method includes: exporting data detected by ground-penetrating radar from software for preprocessing; performing a supplementary set empirical mode decomposition on the preprocessed radar data, performing a Hilbert transform on the first-order component of each intrinsic mode function to obtain an analytical form of the signal, and then obtaining the instantaneous frequency of the radar signal; combining the instantaneous frequencies of each electromagnetic echo signal to form two-dimensional radar data; and plotting an instantaneous frequency graph of the first-order component of the intrinsic mode function of the two-dimensional radar signal based on the two-dimensional radar data, and determining the position of the upper and lower interfaces of frozen soil based on the peak position in the graph. The present invention can accurately and intuitively identify the upper and lower interfaces of frozen soil, overcoming the drawback that electromagnetic wave amplitude and other characteristics are not obvious when ground-penetrating radar is detecting at a deeper location.

Description

Technical Field

[0001] The present invention belongs to the technical field of underground exploration, and relates to a method, a device and a storage medium for identifying the interface of an underground anomaly based on a ground penetrating radar. Background Art

[0002] Currently, the determination of permafrost layers in permafrost areas in engineering geological surveys relies primarily on traditional detection methods, such as landform and topographical identification and drilling. Drilling, a point-based detection method, can effectively and intuitively determine the upper and lower limits of permafrost, but it is expensive and has significant environmental impact. Furthermore, the placement of boreholes is influenced by objective factors such as topography and survey equipment, resulting in certain limitations. Determining the spatial distribution of multiple permafrost areas using traditional methods, as specified in the "Specifications for Engineering Geological Survey of Permafrost Soil" and the "Technical Specifications for the Design and Construction of Highways in Permafrost Areas," requires extensive surveys. Furthermore, if the spacing between detailed surveys is too large, permafrost sections can be missed and inaccurate, posing significant challenges and safety risks to design, construction, and operation. It is crucial to develop a non-destructive ground-penetrating radar (GPR) identification method tailored to the characteristics of permafrost distribution.

[0003] Geophysical exploration methods such as seismic, high-density resistivity, and ground-penetrating radar (GPR) are commonly used. Compared to traditional drilling methods, GPR is non-destructive, provides intuitive results, is accurate, and highly efficient. During GPR field surveys, a transmitting antenna emits pulsed electromagnetic waves at a specific excitation frequency. These pulses, when passing through media with varying electromagnetic properties within the stratum, produce distinct echoes. Engineers can determine the type and size of targets within the stratum based on the reflected wave patterns in the image. However, radar images are simply a reflection of the electric field intensity of the radar electromagnetic wave at different times, creating a "false image." Furthermore, signal clarity in traditional two-dimensional radar grayscale images is affected by numerous factors, including ground surface flatness, waterlogging potential, and radar speed. Furthermore, most targets have irregular shapes, making direct analysis of grayscale images prone to misidentification, missed detection, and inaccurate size determination. For frozen soil, the attenuation of electromagnetic waves and the influence of surrounding clutter can blur the horizon lines or make them invisible due to the clutter. Therefore, radar interpretation requires various processing techniques, such as background removal, time-effect gain, and bandpass filtering, to suppress the clutter and highlight the target. However, radar detection of frozen soil is severely affected by the upper depth and thickness, making it difficult to observe the location of the lower frozen soil interface. Therefore, the Hilbert-Huang transform method is used to process the radar signal and obtain the instantaneous frequency. The time-spectrum diagram analyzes the instantaneous frequency response characteristics when the radar detects frozen soil of different sizes, allowing the identification of the upper and lower frozen soil interfaces and the application of this method in engineering. Summary of the Invention

[0004] In order to solve the above problems, the present invention provides a method for identifying the interface of underground anomalies based on ground penetrating radar, which can accurately and intuitively identify the upper and lower interfaces of underground anomalies, overcoming the defect that the electromagnetic wave amplitude and other characteristics are not obvious when the ground penetrating radar detects at a deeper position.

[0005] Another object of the present invention is to provide an electronic device.

[0006] Another object of the present invention is to provide a computer storage medium.

[0007] The technical solution adopted by the present invention is a method for identifying the interface of an underground anomaly based on ground penetrating radar, which specifically includes the following steps:

[0008] S1: Export the ground penetrating radar field detection data from the software for preprocessing;

[0009] S2: Performing a supplementary ensemble empirical mode decomposition on the preprocessed radar data, performing a Hilbert transform on the first-order component of each intrinsic mode function to obtain the analytical form of the signal, and then obtaining the instantaneous frequency of the radar signal. The instantaneous frequency of each electromagnetic echo signal is combined to form two-dimensional radar data;

[0010] S3: Based on the two-dimensional radar data, a first-order component instantaneous frequency diagram of the intrinsic mode function of the two-dimensional radar signal is drawn, and the upper and lower interface positions of the underground anomaly are determined according to the peak positions in the diagram.

[0011] Furthermore, the step S1 is specifically as follows:

[0012] S1.1: Convert the data format of the radar signal obtained by detection into CSV format;

[0013] S1.2: Pre-process the radar data and take the average value of multiple data in an electromagnetic echo signal based on the number of electromagnetic echo signals in each signal.

[0014] Furthermore, the step S2 is specifically as follows:

[0015] S2.1: Add multiple pairs of positive and negative Gaussian white noise to each electromagnetic echo signal, and then decompose the electromagnetic echo signal to obtain multiple intrinsic mode functions;

[0016] S2.2: Perform a Hilbert transform on the first-order component of each intrinsic mode function to obtain the normalized frequency, converting the signal from the time domain to the frequency domain. The Hilbert transform formula is:

[0017]

[0018] is the signal after Hilbert transformation, H{c i (τ)} represents the content c in {} i (τ) performs Hilbert transform, c i (τ) represents the continuous time signal with independent variable τ, c i (t) represents the continuous time signal with independent variable t after empirical mode decomposition, that is, the i-th intrinsic mode function; the Hilbert transform of the signal can be regarded as the convolution of the signal with 1 / πt, where t represents time;

[0019] S2.3: Multiply the regularized frequency by the sampling frequency to convert it into the true instantaneous frequency;

[0020] S2.4: Combine the real instantaneous frequencies together to form two-dimensional radar data and convert it into a CSV file.

[0021] Furthermore, in step S2.1, the decomposition method for decomposing the electromagnetic echo signal to obtain multiple intrinsic mode functions is an iterative method for obtaining the components of the intrinsic mode functions, specifically:

[0022] S2.11: Find the local maxima of the time series x(t), then interpolate using a third-order spline function to obtain the upper envelope sequence of the time series x(t). Similarly, obtain the lower envelope sequence of the time series. The time series x(t) represents the relationship between time and electromagnetic wave amplitude when radar electromagnetic waves are probing underground.

[0023] S2.12: Average the sum of the upper and lower envelope values of the time series to obtain the instantaneous average value sequence m(t) of the time series envelope;

[0024] S2.13: Subtract the instantaneous mean sequence m(t) from the time series x(t) to obtain the data series h(t);

[0025] S2.14: If the first data sequence h(t) satisfies the definition of an eigenmode function, then the first eigenmode function c1(t) is defined as c1(t) = h(t). If not, then treat h(t) as the time series x(t) in step S2.11 and repeat steps S2.11 to S2.13 until the definition of an eigenmode function is satisfied.

[0026] S2.15: Subtract c1(t) from the time series x(t) to obtain the residual sequence r1(t). Then use the residual sequence r1(t) as the time series x(t) in step S2.11. Repeat steps S2.11 to S2.14 to extract the second, third, and so on, up to the nth intrinsic mode function.

[0027] Furthermore, the intrinsic mode function satisfies the following two conditions:

[0028] The number of extreme points in the signal is equal to the number of zero crossings or differs by at most one;

[0029] The upper envelope formed by the maximum points of the signal and the lower envelope formed by the minimum points are symmetrical about the time axis.

[0030] Furthermore, in step S2.2, the signal is converted from the time domain to the frequency domain, specifically:

[0031] According to the relevant definition of analytical signal, we can get c i (t) The analytical form of the corresponding complex signal:

[0032]

[0033] j represents the imaginary part of the function, z i (t) represents c i (t) the corresponding complex signal;

[0034] Thus we can conclude that:

[0035]

[0036] a i (t) is c i (t) is the instantaneous amplitude function;

[0037]

[0038] c i (t) instantaneous phase function;

[0039]

[0040] ω i (t) is c i (t) is a function of the instantaneous frequency.

[0041] Furthermore, the step S3 is specifically as follows:

[0042] S3.1: Filter the two-dimensional radar data to remove the instantaneous frequencies exceeding the cutoff frequency;

[0043] S3.2: Plot the processed radar data to obtain a graph of time and instantaneous frequency;

[0044] S3.3: Obtain the depth of the target based on the relationship between the ground penetrating radar wave velocity and reflection time;

[0045] S3.4: Convert the graph of time and instantaneous frequency into a graph of detection depth and instantaneous frequency, and then determine the upper and lower interfaces of the underground anomaly.

[0046] Furthermore, the depth calculation formula of the target body is:

[0047]

[0048] Where h is the depth of the target, v is the wave velocity of the ground penetrating radar, and t0 is the time it takes for the ground penetrating radar to receive the electromagnetic wave reflected by the known target.

[0049] An electronic device adopts the above method to realize the interface recognition of underground abnormal bodies.

[0050] A computer storage medium stores at least one program instruction, which is loaded and executed by a processor to implement the above-mentioned underground anomaly interface identification method based on ground penetrating radar.

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

[0052] 1. This invention is suitable for identifying the interface of underground anomalies. The instantaneous frequency curve can accurately and intuitively identify the interface of underground anomalies. It can identify underground anomalies even when the features of ground-penetrating radar images are not obvious and the location range of underground anomalies cannot be manually determined. This promotes research on the influence of electromagnetic parameters on the propagation speed and energy loss of electromagnetic waves in ground-penetrating radar.

[0053] 2. The present invention can overcome the difficulty of sampling in field environments, save manpower and material resources consumed by traditional detection methods, and reduce detection costs and time.

[0054] 3. The existing technology for detecting the thickness of road pavement surface layers is relatively mature, but the detection of deeper underground layers is still relatively vague. The present invention can accurately detect the thickness of layers that cannot be identified by existing ground-penetrating radar images. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0056] Figure 1 This is a flow chart of an embodiment of the present invention.

[0057] Figure 2 This is a forward simulation diagram of the island frozen soil model for the test example of the present invention.

[0058] Figure 3 This is a test example of the present invention, in which a two-dimensional radar image is forward simulated when a 100 MHz radar antenna frequency is used to detect 3 m thick island frozen soil.

[0059] Figure 4 This is the instantaneous frequency image of the first-order IMFs component at a 100 MHz radar antenna frequency, which is a forward simulation of the test example of the present invention.

[0060] Figure 5 This is a diagram showing the measured 100MHz two-dimensional radar detection results of the test example of the present invention.

[0061] Figure 6 This is the instantaneous frequency diagram of the first-order IMFs component of 100 MHz two-dimensional radar data measured in an embodiment of the present invention.

[0062] Figure 7 This is a diagram showing the high-density electrical detection results of an embodiment of the present invention. DETAILED DESCRIPTION

[0063] The following will be combined with the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0064] Example

[0065] A method for identifying underground anomaly interfaces based on ground penetrating radar, such as Figure 1 As shown, the specific steps include:

[0066] S1: Export the ground penetrating radar field detection data from the software for preprocessing;

[0067] S1.1: First, convert the data format of the radar signal obtained by detection into CSV format;

[0068] S1.2: On this basis, the radar data is preprocessed. According to the number of each electromagnetic echo signal, the average value of multiple data in an electromagnetic echo signal is taken. In the embodiment, the average value is taken every 10 data to obtain the processed radar data accordingly. A piece of radar data has more than 10,000 values. After processing, the data becomes less and the workload is reduced.

[0069] S2: Perform a supplementary ensemble empirical mode decomposition on the preprocessed radar data, and then obtain a new analytical form of the signal through Hilbert transform, and then obtain the instantaneous frequency of the radar signal. The instantaneous frequency of each processed electromagnetic echo signal is combined to form two-dimensional radar data.

[0070] S2.1: Perform complementary ensemble empirical mode decomposition (CEEMD): First, multiple pairs of positive and negative Gaussian white noise are added to each electromagnetic echo signal, and then the electromagnetic echo signal is decomposed;

[0071] The decomposition method is an iterative method for obtaining the components of the eigenmode function:

[0072] S2.11: Find the local maxima of the time series x(t), and then interpolate them using a third-order spline function to obtain the upper envelope sequence value of the time series x(t); similarly, obtain the lower envelope sequence value of the time series; where the time series x(t) represents the relationship between time and electromagnetic wave amplitude when radar electromagnetic waves detect underground, which is equivalent to an original signal.

[0073] S2.12: Average the sum of the upper and lower envelope values of the time series x(t) to obtain the instantaneous average value sequence m(t) of the time series envelope.

[0074] S2.13: Subtract the instantaneous mean sequence m(t) from the time series x(t) to obtain the data series h(t).

[0075] S2.14: For different data sequences, h(t) may or may not be an eigenmode function. If the first data sequence h(t) meets the definition of an eigenmode function, then the first eigenmode function c1(t) is defined as c1(t) = h(t). If not, then treat h(t) as the time series x(t) in step S2.11 and repeat steps S2.11 to S2.13 until the definition of an eigenmode function is met.

[0076] Among them, the eigenmode function meets the following two conditions:

[0077] (1) The number of extreme points (maxima and minima) in the signal is equal to the number of zero crossings or differs by at most one;

[0078] (2) The upper envelope formed by the maximum value points of the signal and the lower envelope formed by the minimum value points are symmetrical about the time axis.

[0079] S2.15: Then subtract c1(t) from the time series x(t) to obtain the residual sequence r1(t). Then use the residual sequence r1(t) as a new original sequence (i.e., the time series x(t) in step 1), and repeat the above steps S2.11 to S2.14 to extract the second, third, and so on, up to the nth intrinsic mode function; multiple intrinsic mode functions (IMFs) constitute the original data sequence.

[0080] S2.2: Perform Hilbert transform on the first-order components of each intrinsic mode function obtained in step S2.1 to obtain the regularized frequency, thereby converting the signal from the time domain to the frequency domain.

[0081] The Hilbert transform formula is:

[0082]

[0083] is the signal after Hilbert transformation, H{c i (τ)} represents the c in {} i (τ) performs Hilbert transform, c i (τ) represents the continuous time signal with independent variable τ, c i (t) represents the continuous-time signal with independent variable t after empirical mode decomposition, that is, the i-th intrinsic mode function. The Hilbert transform of the signal can be regarded as the convolution of the signal with 1 / πt, where t represents time.

[0084] Perform the Hilbert transform on the first-order component of each intrinsic mode function (IMF):

[0085]

[0086] According to the relevant definition of analytical signal, the analytical form of time series is obtained:

[0087]

[0088] j represents the imaginary part of the function, z i (t) represents c i (t) the corresponding complex signal;

[0089] Thus we can conclude that:

[0090]

[0091] a i (t) is c i (t) is the instantaneous amplitude function;

[0092]

[0093] c i (t) instantaneous phase function;

[0094]

[0095] ω i (t) is c i (t) is a function of the instantaneous frequency.

[0096] S2.3: Convert the regularized frequency of the first-order component of the intrinsic mode function into the true instantaneous frequency; the instantaneous frequency after Hilbert transform is not the true frequency and needs to be multiplied by the sampling frequency; the two-dimensional radar data is composed of each column of radar data, and the number of data points collected in each column of radar data per second is the sampling frequency.

[0097] S2.4: Combine the true instantaneous frequencies of the first-order components of the intrinsic mode functions to form two-dimensional radar data and convert them into a CSV file.

[0098] S3: The maximum instantaneous frequency of the 2D radar data is used as a characteristic parameter for identifying the frozen soil interface. The 2D radar data is plotted, and the locations of the upper and lower frozen soil interfaces are determined based on the peak positions in the instantaneous frequency plot of the first-order component of the eigenmode function of the 2D radar signal. The maximum instantaneous frequency of the 2D radar data is used as a characteristic parameter for identifying the frozen soil interface. The maximum instantaneous frequency of the first-order component of the eigenmode function corresponds to time zero, the upper frozen soil interface, and the lower frozen soil interface, respectively.

[0099] S3.1: Filter the 2D radar data (i.e., the 2D radar data in step S2.4) to remove instantaneous frequencies exceeding the cutoff frequency. The cutoff frequency is defined as the frequency at which the output signal drops to 0.707 times the maximum value while maintaining the input signal amplitude constant (the cutoff frequency is the -3dB point in frequency response, a special frequency used to describe frequency characteristics). It is typically 0.707 times the radar center frequency.

[0100] S3.2: Plot the processed radar data to obtain a graph of time and instantaneous frequency;

[0101] S3.3: Based on the relationship between the ground penetrating radar wave velocity and reflection time, the depth calculation formula of the target body is obtained:

[0102]

[0103] Where h is the depth of the target, v is the wave velocity of the ground penetrating radar, and t0 is the time it takes for the ground penetrating radar to receive the electromagnetic wave reflected by the known target.

[0104] S3.4: Convert the graph of time and instantaneous frequency into a graph of detection depth and instantaneous frequency. The upper and lower interfaces of the abnormal body (such as frozen soil) can be directly determined based on the graph of detection depth and instantaneous frequency.

[0105] Test case,

[0106] The ground penetrating radar measurement site was selected near a permafrost bridge in the Jialin Highway renovation and expansion project. According to the geotechnical engineering survey report of the Provincial Forestry Design Institute for this project, permafrost exists in this section of the road. The suitable detection depth of 100MHz ground penetrating radar is 8.7m to 10.7m, which can be used to detect permafrost in this section of the road; simulated island permafrost such as Figure 2 As shown, the specific location of the frozen soil in the model is known, which is used to verify whether the frozen soil location identified by the present invention is accurate and effective.

[0107] Step S1: Select a 100 MHz ground penetrating radar, the radar model is the CAS-S100 radar of the Chinese Academy of Sciences, and export the data of the 3 m thick island frozen soil detected by the ground penetrating radar from the software for preprocessing.

[0108] S1.1) 2D radar detection result diagram, see Figure 5 It can be seen that there is a sudden increase in amplitude within the range of 75m to 90m of the survey line, and there are downward dragging survey lines on both sides, indicating that there should be an abnormal body here. In order to avoid omission, the range of 60-90m is selected for supplementary ensemble empirical mode decomposition-Hilbert transform;

[0109] S1.2) First, convert the data format of the detected radar signal into CSV format;

[0110] S1.2) On this basis, the radar data is processed to obtain the characteristics of each electromagnetic echo signal; the average value of each electromagnetic echo signal is taken, generally every 10 data are averaged, and the processed radar data is obtained accordingly; Figure 3 As shown, the radar signature of the lower interface is barely discernible.

[0111] Step S2: Perform a supplementary set empirical mode decomposition on the data in step S1, and then perform a Hilbert transform to obtain a new analytical form of the signal, thereby obtaining the instantaneous frequency of the radar signal. The instantaneous frequency of each processed electromagnetic echo signal is combined to form two-dimensional radar data.

[0112] S2.1) Perform complementary ensemble empirical mode decomposition (CEEMD): First, multiple pairs of positive and negative Gaussian white noise are added to each electromagnetic echo signal, and then the electromagnetic echo signal is decomposed to obtain multiple intrinsic mode functions (IMFs);

[0113] S2.2) Take the first-order component of the IMFs in step S2.1 and perform Hilbert transform to obtain the regularized frequency, thereby converting the signal from the time domain to the frequency domain. The Hilbert transform formula is:

[0114]

[0115] S2.3) converting the regularized frequency of the processed first-order IMFs component into the true instantaneous frequency;

[0116] S2.4) Merge the true instantaneous frequencies of the first-order IMFs components together to form two-dimensional radar data and convert it into a CSV file.

[0117] Step S3: Plot the two-dimensional radar data, and determine the upper and lower interfaces of the frozen soil based on the peak positions in the instantaneous frequency graph of the first-order IMFs component of the two-dimensional radar signal.

[0118] S3.1) filtering the two-dimensional radar data to remove instantaneous frequencies exceeding a cutoff frequency;

[0119] S3.2) Plotting the processed radar data to obtain a graph of time and instantaneous frequency;

[0120] S3.3) Based on the relationship between the ground penetrating radar wave velocity and reflection time, the depth calculation formula of the target body is obtained:

[0121]

[0122] Where h is the depth of the target, v is the wave velocity of the ground penetrating radar, and t0 is the time it takes for the ground penetrating radar to receive the electromagnetic wave reflected by the known target.

[0123] S3.4) Convert the time vs. instantaneous frequency graph into a detection depth vs. instantaneous frequency graph, such as Figure 4 shown.

[0124] S3.5) Figure 4 It can be seen that the positions of the other two peaks below the first peak of the instantaneous frequency of each electromagnetic echo signal are the positions of the upper and lower interfaces. The position of the lower interface can be displayed more clearly, which is consistent with the Figure 2 The results are consistent. Through simulation, it is found that the method of the embodiment of the present invention can more accurately determine the location of permafrost, indicating that the permafrost interface identified by the method of the embodiment of the present invention is accurate and effective. The embodiment of the present invention can identify the location range of underground anomalies when the features of the ground penetrating radar image are not obvious and the location range of the underground anomaly cannot be manually identified.

[0125] Figure 6 The actual frozen soil identified by the method of the embodiment of the present invention is the frozen soil area measured by the high-density electrical method, such as Figure 7 As shown, darker green indicates more frozen soil. Comparison and verification revealed that the results were reliable and the shapes were similar, demonstrating the practical application of the present invention in real projects and verifying its practicability. The distribution of frozen soil in real projects is extremely uneven and irregular, and the method of the present invention can identify this characteristic.

[0126] When the electromagnetic waves emitted by the ground penetrating radar propagate downwards, when they encounter an area with a large dielectric constant value or a deep depth, the electromagnetic wave energy is severely attenuated, which is reflected in the radar image as a significant decrease in the amplitude of the electromagnetic wave, such as Figure 3 , the amplitude of electromagnetic waves at depth becomes almost a straight line. In traditional two-dimensional radar grayscale images, signal clarity is affected by a variety of factors, such as the flatness of the ground, whether the soil is prone to water accumulation, and the walking speed during radar detection. In addition, most targets do not have regular shapes, so directly viewing the grayscale image can easily lead to misjudgments, missed detections, or inaccurate size determination of the target. When using radar for interpretation, engineers need to perform various processing on the radar signal, such as background removal, time-effect gain, and bandpass filtering, to suppress clutter and highlight the target. However, radar detection of permafrost is more seriously affected by the upper depth and thickness. Due to the severe attenuation of electromagnetic wave energy deep underground, the radar amplitude characteristics of the lower interface are extremely unclear, making it difficult to identify the lower interface of the anomaly and difficult to observe the location of the upper and lower interfaces of the permafrost.

[0127] The embodiment of the present invention has found through a large number of experiments and simulation data that the maximum instantaneous frequency of the first-order intrinsic mode function (IMFs) obtained after the radar data is subjected to the supplementary set empirical mode decomposition-Hilbert transform is directly related to the interface, and the upper and lower interfaces of the abnormal body (such as frozen soil) can be directly judged. The three peaks of the instantaneous frequency correspond to time zero and the upper and lower interfaces respectively; the components of other orders have no obvious relationship with the interface position. After a large number of experimental studies, the embodiment of the present invention has established the relationship between the instantaneous frequency of the first-order IMFs component of the radar signal and the frozen soil interface, and the instantaneous frequency is obtained by Hilbert transform, and then the amplitude of the instantaneous frequency is obtained. The amplitude change of the instantaneous frequency is greatly improved compared with the previous original image, and the judgment accuracy is greatly improved. The embodiment of the present invention overcomes the defect that when the ground penetrating radar detection position is deep (0-16m), the electromagnetic wave energy is severely attenuated, resulting in unclear features such as the electromagnetic wave amplitude.

[0128] The embodiments of the present invention are based on measured two-dimensional radar data. The principle of ground-penetrating radar identification is that electromagnetic waves reflect when encountering areas with significantly different dielectric constants. For frozen soil models with different shapes or other parameters, the method of the embodiments of the present invention can still be used to determine the upper and lower frozen soil interfaces, as long as the dielectric constant of the frozen soil differs from that of the surrounding soil. Ground-penetrating radar operates at a frequency of 0.1 GHz to 1.6 GHz and primarily detects the internal structure of roads, determining the presence of defects beneath the roadbed and pavement. It can also be used to identify municipal pipelines, underground cavities, groundwater, and other anomalies.

[0129] If the method for identifying an underground anomaly interface described in the embodiment of the present invention is implemented as a software functional module and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the portion that contributes to the prior art, or the portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions for causing a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the steps of the method for identifying an underground anomaly interface described in the embodiment of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as a USB flash drive, a mobile hard drive, ROM, RAM, a magnetic disk, or an optical disk.

[0130] The above description is only a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention are included in the scope of protection of the present invention.

Claims

1. A method for identifying underground anomaly interfaces based on ground penetrating radar, characterized in that: The specific steps include: S1: Export the ground penetrating radar field detection data from the software for preprocessing; S2: Performing a supplementary ensemble empirical mode decomposition on the preprocessed radar data, performing a Hilbert transform on the first-order component of each intrinsic mode function to obtain the analytical form of the signal, and then obtaining the instantaneous frequency of the radar signal. The instantaneous frequency of each electromagnetic echo signal is combined to form two-dimensional radar data; S3: Based on the two-dimensional radar data, a first-order component instantaneous frequency diagram of the intrinsic mode function of the two-dimensional radar signal is drawn, and the upper and lower interface positions of the underground anomaly are determined according to the peak positions in the diagram; The step S2 is specifically as follows: S2.1: Add multiple pairs of positive and negative Gaussian white noise to each electromagnetic echo signal, and then decompose the electromagnetic echo signal to obtain multiple intrinsic mode functions; S2.2: Perform a Hilbert transform on the first-order component of each intrinsic mode function to obtain the normalized frequency, converting the signal from the time domain to the frequency domain. The Hilbert transform formula is: is the signal after Hilbert transformation, H{c i (τ)} represents the content c in {} i (τ) performs Hilbert transform, c i (τ) represents the continuous time signal with independent variable τ, c i (t) represents the continuous time signal with independent variable t after empirical mode decomposition, that is, the i-th intrinsic mode function; the Hilbert transform of the signal can be regarded as the convolution of the signal with 1 / πt, where t represents time; S2.3: Multiply the regularized frequency by the sampling frequency to convert it into the true instantaneous frequency; S2.4: Merge the real instantaneous frequencies together to form two-dimensional radar data and convert it into a CSV file; The step S3 is specifically as follows: S3.1: Filter the two-dimensional radar data to remove the instantaneous frequencies exceeding the cutoff frequency; S3.2: Plot the processed radar data to obtain a graph of time and instantaneous frequency; S3.3: Obtain the depth of the target based on the relationship between the ground penetrating radar wave velocity and reflection time; S3.4: Convert the graph of time and instantaneous frequency into a graph of detection depth and instantaneous frequency, and then determine the upper and lower interfaces of the underground anomaly.

2. The method for identifying underground anomaly interfaces based on ground penetrating radar according to claim 1, characterized in that: The step S1 is specifically as follows: S1.1: Convert the data format of the radar signal obtained by detection into CSV format; S1.2: Pre-process the radar data and take the average value of multiple data in an electromagnetic echo signal based on the number of electromagnetic echo signals in each signal.

3. The method for identifying underground anomaly interfaces based on ground penetrating radar according to claim 1, characterized in that: In step S2.1, the decomposition method for decomposing the electromagnetic echo signal to obtain multiple eigenmode functions is an iterative method for obtaining eigenmode function components, specifically: S2.11: Find the local maxima of the time series x(t), then interpolate using a third-order spline function to obtain the upper envelope sequence of the time series x(t). Similarly, obtain the lower envelope sequence of the time series. The time series x(t) represents the relationship between time and electromagnetic wave amplitude when radar electromagnetic waves are probing underground. S2.12: Average the sum of the upper and lower envelope values of the time series to obtain the instantaneous average value sequence m(t) of the time series envelope; S2.13: Subtract the instantaneous mean sequence m(t) from the time series x(t) to obtain the data series h(t); S2.14: If the first data sequence h(t) satisfies the definition of an eigenmode function, then the first eigenmode function c1(t) is defined as c1(t) = h(t). If not, then treat h(t) as the time series x(t) in step S2.11 and repeat steps S2.11 to S2.13 until the definition of an eigenmode function is satisfied. S2.15: Subtract c1(t) from the time series x(t) to obtain the residual sequence r1(t). Then use the residual sequence r1(t) as the time series x(t) in step S2.

11. Repeat steps S2.11 to S2.14 to extract the second, third, and so on, up to the nth intrinsic mode function.

4. The method for identifying underground anomaly interfaces based on ground penetrating radar according to claim 3, characterized in that: Among them, The eigenmode function satisfies the following two conditions: The number of extreme points in the signal is equal to the number of zero crossings or differs by at most one; The upper envelope formed by the maximum points of the signal and the lower envelope formed by the minimum points are symmetrical about the time axis.

5. The method for identifying underground anomaly interfaces based on ground penetrating radar according to claim 1, characterized in that: In step S2.2, the signal is converted from the time domain to the frequency domain, specifically: According to the relevant definition of analytical signal, we can get c i (t) The analytical form of the corresponding complex signal: j represents the imaginary part of the function, z i (t) represents c i (t) the corresponding complex signal; Thus we can conclude that: a i (t) is c i (t) is the instantaneous amplitude function; c i (t) instantaneous phase function; ω i (t) is c i (t) is a function of the instantaneous frequency.

6. The method for identifying underground anomaly interfaces based on ground penetrating radar according to claim 5, characterized in that: The depth calculation formula of the target body is: Where h is the depth of the target, v is the wave velocity of the ground penetrating radar, and t0 is the time it takes for the ground penetrating radar to receive the electromagnetic wave reflected by the known target.

7. An electronic device, characterized in that: The method according to any one of claims 1 to 6 is used to realize the interface identification of underground abnormal bodies.

8. A computer storage medium, characterized in that The storage medium stores at least one program instruction, and the at least one program instruction is loaded and executed by the processor to implement the underground anomaly interface identification method based on ground penetrating radar according to any one of claims 1 to 6.

Images