A method for detecting abnormal ground area under the assistance of ground penetrating radar real-time interpretation

Abstract

This invention relates to the field of radar detection technology, specifically to a method for detecting abnormal ground areas with the assistance of real-time interpretation of ground-penetrating radar. The invention first acquires the local abnormal signal matrix and reference waveform set to be analyzed, removing vehicle-borne dynamic interference; further, it extracts reflection feature units, analyzes the correlation between the local waveform data corresponding to the reflection feature units and the corresponding reference waveform set, and obtains waveform feature vectors; further, based on the differences between reflection feature units in adjacent columns, it constructs a cost matrix and obtains the optimal matching result; further, based on the optimal matching result and combined with the distribution of the time points of the corresponding reflection feature units in the local abnormal signal matrix, it constructs a diagnostic feature map; finally, it archives the diagnostic feature map into a database, and based on the spatiotemporal comparison of the diagnostic feature map of the current data segment with historical records, it stably identifies diffuse anomalies within a unified spatiotemporal framework across periods, generating a diagnostic report.

Description

Technical Field

[0001] This invention relates to the field of radar detection technology, and specifically to a method for detecting abnormal ground areas with the assistance of real-time interpretation of ground penetrating radar. Background Technology

[0002] Ground-penetrating radar (GPR) is widely used for road structure condition detection due to its continuous, rapid, and non-destructive characteristics. However, in high-speed vehicle-mounted detection scenarios, road undulations, vehicle movement, and environmental noise introduce significant low-frequency interference during data acquisition, causing distortions such as overall radar signal shift and baseline fluctuations, thereby masking the true underground reflection information.

[0003] On the other hand, common anomalies in road structures include both local anomalies with clear geometric boundaries, such as faults and voids, and diffuse material property anomalies caused by changes in moisture content and density. These diffuse anomalies often do not exhibit typical hyperbolic reflections in radar images, but rather show a slow evolution of reflection energy and local waveform morphology over time. This leads to missed detections and false positives in traditional automated detection methods that rely on geometric morphology or sudden changes in local energy.

[0004] Traditional automated interpretation methods mostly focus on identifying local anomalies with clear geometric shapes. For diffuse anomalies characterized by gradual changes in attributes, there is a lack of stable detection methods, which often leads to missed detections. Summary of the Invention

[0005] To address the technical problem of unstable identification of diffuse anomalies due to dynamic interference in vehicle-mounted ground-penetrating radar (GPR) data acquisition, this invention aims to provide a method for detecting abnormal ground areas with the assistance of real-time GPR interpretation. The specific technical solution adopted is as follows:

[0006] A data matrix of a segment of data collected by ground penetrating radar along the survey line direction is obtained. Singular value decomposition is performed on the data matrix to extract the local anomaly signal matrix and reference waveform set.

[0007] Reflection feature units are extracted based on the data fluctuations in each column of the local anomaly signal matrix; the correlation between the local waveform data corresponding to the reflection feature unit and the reference waveform set extracted from the previous data segment is analyzed to obtain the waveform feature vector of each reflection feature unit in the current data segment; based on the differences in time and amplitude distribution between the reflection feature units in two adjacent columns, and combined with the differences between the waveform feature vectors, a cost matrix corresponding to the two columns of data is constructed and the optimal matching result is obtained; based on the optimal matching result, and combined with the distribution of the time points of the corresponding reflection feature units in the local anomaly signal matrix, a diagnostic feature map is constructed.

[0008] The diagnostic feature maps are archived and stored in a database. A diagnostic report is generated by comparing the diagnostic feature maps of the current data segment with historical records in a spatiotemporal manner.

[0009] Furthermore, the method for obtaining the local anomaly signal matrix includes:

[0010] Based on the first left singular vector, the first right singular vector, and the first singular value obtained from the singular value decomposition, a background matrix is ​​reconstructed; based on the background matrix, the background is removed from the data matrix to obtain a local anomaly signal matrix.

[0011] Furthermore, the method for obtaining the reference waveform set includes:

[0012] In the left singular vector matrix obtained by the singular value decomposition, a preset number of column parameters are selected starting from the second column to form a reference waveform set.

[0013] Furthermore, the method for obtaining the reflection feature unit includes:

[0014] For each column of data in the local abnormal signal matrix, a one-dimensional peak detection algorithm is used to obtain the peak, and effective peaks are selected based on a preset amplitude threshold and a preset peak spacing threshold. The time sampling point of the peak value of an effective peak, the normalized amplitude value of the peak value, the index number of the column in the local abnormal signal matrix, and the column number are used to form a reflection feature unit.

[0015] Furthermore, the method for obtaining the waveform feature vector includes:

[0016] For each of the reflection feature units, with the time sampling point corresponding to the reflection feature unit as the center, data of a preset waveform segment length is extracted from the corresponding column of the local abnormal signal matrix as the local waveform data corresponding to the reflection feature unit.

[0017] The inner product is calculated between the local waveform data and each vector in the reference waveform set extracted from the previous segment of data, and the inner product values ​​are arranged to form the waveform feature vector.

[0018] Furthermore, the method for obtaining the cost matrix includes:

[0019] The reflection feature unit is selected from each of the two adjacent columns to form a binary tuple; the time-amplitude difference cost is obtained based on the difference between the time sampling points and the difference between the amplitude values ​​of the two elements in the binary tuple; the morphological difference cost is obtained based on the difference between the waveform feature vectors of the two elements in the binary tuple; the time-amplitude difference cost and the morphological difference cost are fused to obtain the total cost.

[0020] The total value of all pairs of the two adjacent columns is used to construct a cost matrix.

[0021] Furthermore, the method for obtaining the diagnostic feature map includes:

[0022] Based on the local anomaly signal matrix, construct a structural boundary map and a medium gradient indicator map of consistent size;

[0023] Based on the optimal matching result, according to the attribute of the unmatched reflection feature unit as the starting point or ending point of the path, the starting point mark or ending point mark is filled in at the corresponding position of the time sampling point to generate a structural boundary map.

[0024] By filling the morphological difference value of each successfully matched pair of reflection feature units into the corresponding position of the time sampling point of the reflection feature unit with the smallest temporal sequence, a medium gradient indicator map is generated.

[0025] Furthermore, the method for archiving and building a library of the diagnostic feature maps includes:

[0026] The diagnostic feature map, the geographic coordinates of the starting and ending points of the currently analyzed data segment, and the timestamp of the start time of the detection task are combined into a data record unit and stored in a structured spatiotemporal database.

[0027] Furthermore, methods for generating diagnostic reports include:

[0028] Retrieve all feature map data generated by the latest detection task and extract candidate abnormal regions from each feature map;

[0029] The candidate abnormal regions of the same type are associated across periods, and abnormal record entries are updated or created in the database.

[0030] Anomaly detection records with multiple detected attributes are retrieved. When the anomaly detection record corresponds to a structural type, univariate linear regression analysis is applied to the region area sequence and its corresponding timestamp sequence to calculate the slope of the regression line, which is used as the geometric expansion rate. The risk level is then determined based on the latest region area and geometric expansion rate. When the anomaly detection record corresponds to a physical property type, univariate linear regression analysis is applied to the maximum data value sequence and its corresponding timestamp sequence to calculate the slope of the regression line, which is used as the feature evolution rate. The risk level is then determined based on the latest maximum feature value and feature evolution rate.

[0031] A diagnostic report is generated based on the risk level classification and anomaly detection records.

[0032] Furthermore, the data matrix undergoes preprocessing, which includes:

[0033] The collected data is segmented into a data matrix by sliding window along the survey line. DC drift correction is performed on each column of data, and the time shift of surface reflection is calculated based on cross-correlation to perform surface reflection alignment.

[0034] The present invention has the following beneficial effects:

[0035] This invention first acquires the local anomaly signal matrix to be analyzed and a reference waveform set to suppress background interference and highlight local anomaly features. It then extracts reflection feature units to effectively capture possible anomaly events, providing a foundation for subsequent analysis. The correlation between the local waveform data corresponding to the reflection feature units and the corresponding reference waveform set is analyzed to obtain waveform feature vectors, characterizing the morphological properties of the reflection feature units and their similarity to the reference waveforms, facilitating more accurate anomaly classification and diagnosis. Furthermore, based on the differences in temporal and amplitude distributions between adjacent columns of reflection feature units, combined with the differences between waveform feature vectors, a cost matrix is ​​constructed for the corresponding two columns of data to assess the probability that two reflection feature units belong to the same evolutionary path and obtain the optimal matching result, providing a basis for constructing a diagnostic feature map. Based on the optimal matching result and the distribution of the corresponding reflection feature units' time points in the local anomaly signal matrix, a diagnostic feature map is constructed, providing a basis for generating a final diagnostic report. Finally, the diagnostic feature maps are archived and stored. By comparing the diagnostic feature map of the current data segment with historical records in a unified spatiotemporal framework across periods, the cumulative change trend is stably identified, generating a diagnostic report. This method performs singular value decomposition on the radar data matrix to remove vehicle dynamic interference; extracts reflection feature units column by column and completes adjacent event matching based on waveform correlation to form a diagnostic feature map; and then performs spatiotemporal comparison to generate a diagnostic report, thereby stably identifying diffuse anomalies. Attached Figure Description

[0036] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0037] Figure 1 A flowchart illustrating a method for detecting abnormal ground areas with real-time interpretation assistance from ground-penetrating radar, provided in an embodiment of the present invention;

[0038] Figure 2 This is a flowchart illustrating the process of obtaining a cost matrix according to an embodiment of the present invention. Detailed Implementation

[0039] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the following, in conjunction with the accompanying drawings and preferred embodiments, details the specific implementation, structure, features, and effects of a ground-penetrating radar real-time interpretation-assisted method for detecting abnormal ground areas according to the present invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0040] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0041] The following description, in conjunction with the accompanying drawings, details a specific scheme for an abnormal ground area detection method assisted by real-time ground-penetrating radar provided by this invention.

[0042] Please see Figure 1 The diagram illustrates a flowchart of an abnormal ground area detection method assisted by real-time ground-penetrating radar interpretation, according to an embodiment of the present invention, specifically including:

[0043] Step S1: Obtain the data matrix of a segment of data collected by the ground penetrating radar along the survey line direction, perform singular value decomposition on the data matrix, and extract the local abnormal signal matrix and reference waveform set.

[0044] In one embodiment of the present invention, a data matrix of a segment of data collected by ground penetrating radar along the survey line direction is first obtained to provide a basis for data analysis.

[0045] In a preferred embodiment of the present invention, considering that the original ground-penetrating radar data is acquired in the form of a continuous data stream and includes DC drift introduced by the instrument and misalignment of the signals on the time axis caused by small changes in road surface elevation, these factors can interfere with the accuracy of subsequent analysis. Therefore, it is necessary to perform segmentation and standardization preprocessing on the data stream.

[0046] Based on this, the collected data is segmented into a data matrix by sliding window along the survey line direction, DC drift correction is performed on each column of data, and the time shift of surface reflection is calculated based on cross correlation to perform surface reflection alignment.

[0047] Specifically, a sliding window process is applied to the collected data to generate a series of data segments. Each data segment contains N consecutive A-scan data points, forming a data set of size [missing information]. Two-dimensional data matrix M is the number of time sampling points for each A-scan, N is the number of A-scan channels in each data segment, and i represents the sequence number of the data segment;

[0048] As an example, the sliding occurs in the "survey line direction". The width of the sliding window is N, which is the number of consecutive A-scans N contained in each data segment. N is set to 256, and the movement step is also N. The A-scan spacing is 4 centimeters, corresponding to a road length of 10.24 meters (about 10 meters).

[0049] Furthermore, regarding the matrix For each column (i.e. each A-scan data point), the average value of all sampled data points in that column is subtracted to perform DC drift correction;

[0050] Finally, let's take the first A-scan data (matrix) from the current segment of data being analyzed. Using the surface reflection location of the first column as a baseline, the cross-correlation function between each subsequent A-scan data point and this baseline A-scan data point is calculated to find the time shift corresponding to the peak value of each A-scan data point. Then, the entire A-scan data column is shifted to obtain a preprocessed data point whose size remains unchanged. The data matrix is ​​denoted as the data matrix corresponding to the segment of data being analyzed.

[0051] It should be noted that the calculation of the cross-correlation function is a well-known technique in the field. The specific value of M is determined by the time window length and sampling frequency of the radar system. N can be configured according to the system memory and real-time requirements. The moving step size (or overlap rate) of the sliding window can be adjusted by the implementer. This is common knowledge in the field and will not be elaborated further.

[0052] Considering the significant differences in signal components within the data matrix: the background signal, composed of reflections from stable layered structures and the macroscopic trend of vehicle bumps, has strong energy and varies continuously in space (along the survey line direction); while the disturbance caused by local underground anomalies (such as loosening and micro-cracks) has weaker energy and does not possess universal continuity in space.

[0053] Since Singular Value Decomposition (SVD) can sort and decompose matrix components according to the energy contribution and mode stability of the signal in the matrix, performing SVD on the data matrix to extract the local anomalous signal matrix and reference waveform set can effectively suppress background interference, highlight local anomalous features, and provide a reliable foundation for subsequent reflection feature unit extraction and waveform matching.

[0054] Preferably, in one embodiment of the present invention, considering that the macroscopic background signal composed of stable layered structure reflection and vehicle bumps has the characteristics of high energy and continuous spatial variation, the dominant singular component of the matrix corresponding to this background signal is the first left singular vector obtained based on singular value decomposition. First right singular vector and the first singular value Reconstruct to obtain the background matrix This characterizes macroscopic background interference;

[0055] Background is removed from the data matrix based on the background matrix to obtain the local anomaly signal matrix. .

[0056] Specifically, ,matrix Macroeconomic background interference has been eliminated, and the current data segment has been retained. Details of all signals related to local anomalies.

[0057] Considering that each column of the left singular vector matrix is ​​an orthogonal basis vector, capable of representing the main reflection waveform change patterns within the current data segment, a predetermined number of column parameters are selected from the second column of the left singular vector matrix obtained by singular value decomposition to form a reference waveform set. .

[0058] As an example, the preset column number parameter K can be an integer between 5 and 10. In this example, it is set to 5, which corresponds to columns 2 to 6 (column K+1). Column 1 is skipped to avoid treating the macro background as a reference waveform.

[0059] It should be noted that using singular value decomposition to decompose a matrix is ​​a technique well-known to those skilled in the art, and will not be elaborated further.

[0060] It should be noted that one embodiment of the present invention employs a hysteresis adaptive mechanism, where processing the current data segment relies on some parameters generated from the previous data segment (reference waveform set and time distribution variance). and amplitude distribution variance Therefore, for the first data segment (index) Because a parameter-bootstrapping initialization procedure must be executed in the preceding data segment, this procedure is executed only once when processing the first data segment.

[0061] Specifically: the reference waveform set obtained based on the first segment of data. , , It is used directly in the subsequent first stage of data analysis and stored for use in the next stage of data analysis.

[0062] The calculation methods for the time distribution variance and the amplitude distribution variance will be explained in subsequent steps.

[0063] Step S2: Extract reflection feature units based on the data fluctuations in each column of the local abnormal signal matrix; analyze the correlation between the local waveform data corresponding to the reflection feature unit and the reference waveform set extracted from the previous data segment to obtain the waveform feature vector of each reflection feature unit in the current data segment; construct the cost matrix corresponding to the two columns of data and obtain the optimal matching result based on the differences in time distribution and amplitude distribution between reflection feature units in adjacent columns and the differences between waveform feature vectors; construct a diagnostic feature map based on the optimal matching result and the distribution of the time points of the corresponding reflection feature units in the local abnormal signal matrix.

[0064] Existing methods cannot simultaneously distinguish, within a unified framework, the geometrical interruption of reflected signals caused by structural damage (such as cracks or voids) and the morphological evolution of reflected signals caused by material degradation (such as gradual changes in moisture content or decreased density). Therefore, step S2 constructs a diagnostic feature map to provide an analytical basis for subsequent synchronous identification.

[0065] Considering that the fluctuation characteristics of local anomalous signals and background or macro trend signals are different in a series of data, we first extract reflection feature units based on the data fluctuations in each column of data in the local anomalous signal matrix. This can effectively capture possible anomalous events and provide a basis for subsequent analysis.

[0066] Preferably, in one embodiment of the present invention, the local abnormal signal is manifested as a sudden change in amplitude or a change in local shape in the time-space two-dimensional matrix, while the background signal and macro trend are continuous and smooth in space, so the abnormal event presents as a peak shape in a series of data;

[0067] Based on this, for each column of data in the local abnormal signal matrix, a one-dimensional peak detection algorithm is used to obtain the peak. At the same time, to reduce false peaks caused by random noise, a preset amplitude threshold is also used. and preset peak spacing threshold Valid peaks were selected;

[0068] As an example, It can be set to 1.5 times the signal standard deviation. It can be set to one-quarter of the wavelength corresponding to the center frequency. In the first round of screening, when the peak value of the detected wave is less than... In the first round, peaks are filtered out; in the second round, the peaks retained in the first round are arranged in chronological order, and the distance between the peak values ​​of two adjacent peaks is checked. If it is less than... Then, only the one with the larger amplitude (peak value) is retained, and the other one is removed; after two rounds of screening, the effective peak value is obtained.

[0069] It should be noted that the center frequency refers to the main frequency of the electromagnetic pulse emitted by the radar system, which is usually determined by the antenna and the transmitted waveform. The methods for obtaining the signal standard deviation and the center frequency of the ground penetrating radar system are well-known technologies in the art and will not be elaborated here.

[0070] A reflection feature unit is formed by taking the time sampling point t of the peak value of an effective wave peak, the normalized amplitude value a of the peak value, the column index x and the column number j in the local anomalous signal matrix.

[0071] The normalized amplitude value is the peak value after linear normalization, which is normalized in the data dimension composed of the peak values ​​of all valid peaks; the column numbers are sorted according to the time order of the valid peaks (peak points) in a column.

[0072] Considering reflections from different causes, even with similar time and amplitude, the morphology (including phase, broadening, and waveform details) of the reflected wavelets differs, and these morphological differences are directly related to the medium properties and structural state of the reflector. To analyze and distinguish these morphological differences, the correlation between the local waveform data corresponding to the reflection feature unit and the reference waveform set extracted from the previous data segment is analyzed. This yields the waveform feature vector for each reflection feature unit in the current data segment, characterizing its morphological properties and its similarity to the reference waveform. This not only enhances the ability to identify minute or diffuse anomalies but also distinguishes different media or damage types when time and amplitude are similar, achieving more accurate anomaly classification and diagnosis.

[0073] Preferably, in one embodiment of the present invention, for each reflection feature unit, data of a preset waveform segment length is extracted from the corresponding column of the local abnormal signal matrix, centered on the time sampling point corresponding to the reflection feature unit, as the local waveform data corresponding to the reflection feature unit; this method of obtaining local waveform data can ensure that the extracted segment completely covers the main wavelet components of the reflection event, avoid the loss of waveform features due to insufficient extraction, and at the same time preserve its local morphological features, which is convenient for subsequent analysis.

[0074] As an example, the preset waveform segment length L is 1.5 times the wavelength corresponding to the center frequency, in order to cover a complete sub-wavelength and capture the main waveform details;

[0075] Furthermore, considering that the morphological differences of different reflectors directly reflect the medium properties and structural state, the local waveform data is used to calculate the inner product with each vector in the reference waveform set extracted from the previous data segment, and the inner product values ​​are arranged to form the waveform feature vector.

[0076] It should be noted that each column of data in the reference waveform set corresponds to a vector. When calculating the inner product, the longest vector among the two vectors is used as the reference, and the beginning ends are aligned. Zeros are padded at the end of the shortest vector to align the lengths at both ends. The inner product is used to show the correlation between the local waveform data and each vector in the reference waveform set. After arrangement, the overall correlation is shown.

[0077] In this process, ground-penetrating radar continuously collects data along the survey line. The reflected signals of adjacent data segments usually originate from the same type of underground structure or material properties, and the waveform morphology has a certain continuity. By using the reference waveform set of the previous segment as a benchmark, the local waveform of the current segment can be compared with the typical waveform of the previous segment. Using the waveform set of the previous segment as a stable reference helps to eliminate the interference of single-segment data anomalies and highlight the characteristics of local anomalies, corresponding to the "hysteresis adaptive" core mechanism adopted.

[0078] Since the reflection events generated by a continuous underground anomaly in two adjacent columns (adjacent channels, sampling points) should not only remain continuous in spatial location, but their reflection waveform characteristics should also remain highly similar or evolve smoothly, a cost matrix corresponding to the two columns of data is constructed based on the differences in temporal and amplitude distributions between the reflection feature units in two adjacent columns, combined with the differences between waveform feature vectors. This is to assess the possibility that two events (reflection feature units) belong to the same evolution path and obtain the optimal matching result, providing a basis for subsequent path tracing and the construction of diagnostic feature maps.

[0079] Preferably, in one embodiment of the present invention, please refer to Figure 2 The diagram illustrates a flowchart of a cost matrix acquisition method provided by an embodiment of the present invention, specifically including:

[0080] Step S201: Select one reflection feature unit from each of the two adjacent columns to form a binary tuple.

[0081] First, arbitrarily select two adjacent data columns in the local anomalous signal matrix. Then, take one reflection feature unit from each of the two adjacent columns to form a binary pair for matching analysis. The analysis process is the same for all binary pairs in all adjacent columns in the local anomalous signal matrix. Only one example is described here and will not be repeated.

[0082] Step S202: Obtain the time-amplitude difference cost based on the difference between the time sampling points and the difference between the amplitude values ​​of the two elements in the tuple.

[0083] For the two elements in a binary tuple, the differences are analyzed along two dimensions: time sampling point and amplitude value. As an example, a lag adaptive mechanism is also used, with the index of the current data segment being i, and the variance of the time distribution obtained when processing the previous data segment is retrieved. and amplitude distribution variance ;

[0084] Using the squared Euclidean distance between the time sampling points of the two elements as the numerator, As the denominator, the ratio of the fractions serves as the time-dependent value, characterizing the differences in temporal distribution among the reflective feature units; the squared Euclidean distance between the normalized amplitude values ​​of the two elements is used as the numerator. As the denominator, the fractional ratio serves as the amplitude cost, characterizing the differences in amplitude distribution among the reflection feature units; the sum of the time cost and the amplitude cost is used as the time-amplitude difference cost, which allows the matching cost to adapt to the overall fluctuation level of signals in different road segments, avoiding the problem of insufficient or oversensitive sensitivity of fixed thresholds under different geological conditions.

[0085] It should be noted that the variance of the time sampling points of all reflection feature units in the currently analyzed data segment is taken as the time distribution variance, and the variance of the normalized amplitude values ​​of all reflection feature units is taken as the amplitude distribution variance. These two variance values ​​are stored for use in processing the next data segment. The first data segment uses the two variance values ​​calculated from its own data.

[0086] In one embodiment of the present invention, a preset positive parameter for division by zero, such as 0.01, can be added to the denominator. The unit of the parameter is the same as the unit of the denominator to avoid the denominator being zero.

[0087] Step S203: Obtain the morphological difference value based on the difference between the waveform feature vectors of the two elements in the binary tuple.

[0088] Further analysis of the differences between waveform feature vectors reveals that, as an example, the squared Euclidean distance between two waveform feature vectors is used as the morphological difference denominator to characterize the differences between waveform feature vectors.

[0089] Step S204: Integrate the time-amplitude difference value and the morphological difference value to obtain the total value; construct a cost matrix from the total values ​​of all pairs of adjacent columns.

[0090] Finally, the generation value quantified from different dimensions is integrated to obtain the total generation value, which characterizes the associated generation value of the two reflective feature units of the analyzed binary belonging to the same evolutionary path, reflecting the possibility of belonging to the same evolutionary path.

[0091] As an example, the product of the preset first weight and the time-amplitude difference cost value, plus the product of the preset second weight and the morphological difference cost value, is used as the total value of the analyzed binary pair; wherein, the preset first weight and the preset second weight are used to balance the relative importance of the two costs, and in this example, both can be set to 1; in other embodiments of the present invention, the implementer can adjust them as needed.

[0092] Furthermore, the total value of all pairs of corresponding binary pairs in two adjacent columns is used to construct a cost matrix. The size is , where m is the number of reflection feature units (events) in the xth column (channel) and n is the number of reflection feature units (events) in the (x+1)th column (channel).

[0093] Finally, the optimal matching result for the globally optimal assignment of the cost matrix is ​​obtained using the Hungarian algorithm. The Hungarian algorithm is a well-known combinatorial optimization algorithm for solving assignment problems. When the number of events in two paths is not equal, the cost matrix is ​​expanded into a square matrix by supplementing it with virtual rows or columns with maximum cost values ​​to ensure that the algorithm can handle unbalanced matching problems.

[0094] The optimal matching result contains a set of matching pairs, which includes all successfully matched event pairs (reflection feature unit pairs); it also includes a set of unmatched events, which includes all unmatched events (reflection feature units).

[0095] Since different types of underground anomalies behave differently during the matching process, in order to separate and visualize the differences in behavior, a diagnostic feature map is constructed based on the optimal matching result and the distribution of the time points of the corresponding reflection feature units in the local anomaly signal matrix, providing a basis for the final generation of a diagnostic report.

[0096] Preferably, in one embodiment of the present invention, a structural boundary map and a medium gradient indicator map of the same size are constructed based on the local anomaly signal matrix. When these two feature maps are created, the initial value of each position is set to zero, and the time attributes of each position of each map are the same as those of the same position in the local anomaly signal matrix, corresponding to the same time sampling point.

[0097] Based on the optimal matching result, according to the attribute of the unmatched reflection feature unit as the starting point or ending point of the path, the starting point mark or ending point mark is filled in at the corresponding position of the time sampling point to generate a structural boundary map.

[0098] If the reflection feature unit j in column x fails to find a matching reflection feature unit in column x+1, then the reflection feature unit j in column x is the end point of an evolution path, and the end point marker -1 is filled in at the corresponding position of the end point time sampling point; if the reflection feature unit k in column x+1 fails to find a matching reflection feature unit in column x, then the reflection feature unit k in column x+1 is the start point of an evolution path, and the start point marker +1 is filled in at the corresponding position of the start point time sampling point.

[0099] By marking the start and end points of the path, the outlines of anomalous bodies with clear physical boundaries, such as cavities and pipelines, are directly indicated.

[0100] By filling the morphological difference value of each successfully matched pair of reflective feature units into the corresponding position of the time sampling point of the reflective feature unit with the smallest temporal sequence, a medium gradient indicator map is generated. The greater the morphological difference, the larger the data value is filled in, thus appearing as a highlighted area in the map, indicating the presence of a gradient in medium properties at the corresponding position, thereby achieving stable detection of such anomalies.

[0101] It should be noted that when a reflective feature unit serves as both the start and end point, based on the probability of other reflective feature units in the eight neighborhoods of this reflective feature unit serving as both the start and end point in the generated medium gradient indicator map, the one with the highest probability is selected as the final attribute.

[0102] It should be noted that, in another embodiment of the present invention, a path tracing cost map of consistent size is also constructed based on the local abnormal signal matrix; based on the optimal matching result, the path tracing cost map is generated by filling the total value of each pair of successfully matched reflection feature units into the corresponding position of the time sampling point of the reflection feature unit with the smallest time sequence; the higher the data value of this map, the more it indicates that the event evolution path is discontinuous, abrupt, or drastic, providing a basis for abnormal warning.

[0103] Generate diagnostic feature maps with clear engineering implications, complete the "real-time interpretation" of data collected by ground penetrating radar, and provide a basis for the final generation of diagnostic reports.

[0104] Step S3: Archive the diagnostic feature maps into a database, and generate a diagnostic report based on the spatiotemporal comparison of the diagnostic feature maps of the current data segment with the historical records.

[0105] Considering that diffuse anomalies are usually characterized by slow evolution, weak reflection, and unstable morphology, their cumulative change trend can only be stably identified within a unified spatiotemporal framework across cycles. Therefore, diagnostic feature maps are archived and stored in a database.

[0106] Considering that the evolution of underground structures has obvious temporal and spatial dependence, the diagnostic feature maps of the same location show different characteristics in different detection cycles, reflecting the evolution trend of the same reflection event in different cycles. Thus, a specific underground disease can be tracked as an independent entity. Therefore, by comparing the diagnostic feature map of the current data with the historical data in time and space, a diagnostic report can be generated to achieve a quantitative and trend assessment of the development of structural damage, the cumulative effect of material deterioration, and the emergence of new anomalies. This helps to identify potential risks in advance and realize long-term dynamic monitoring and early warning of the health status of underground structures.

[0107] Preferably, in one embodiment of the present invention, a data recording unit is constructed by the diagnostic feature map, the geographic coordinates of the starting and ending points corresponding to the currently analyzed segment of data, and the timestamp of the starting time of the detection task, and stored in a structured spatiotemporal database.

[0108] Preferably, in one embodiment of the present invention, all feature map data generated by the latest detection task are retrieved, and candidate abnormal regions are extracted from each feature map;

[0109] Specifically, on each structural boundary map, a connected component labeling algorithm, a well-known technique in image processing, is applied to aggregate spatially adjacent non-zero pixels. Spatially adjacent +1 and -1 labels are grouped together to form one or more discrete regions. Each discrete region is defined as a "candidate structural anomaly region," and its corresponding geographic bounding box and area are calculated.

[0110] Among them, the geographic location bounding box: obtain the column range and row range covered by this region in the data matrix, and combine the geographic coordinates associated with the column index and the depth information associated with the row index to determine the minimum enclosing rectangle of this region in geographic space, thus obtaining the geographic location bounding box.

[0111] Region area: The approximate physical area of ​​the abnormal region is obtained by multiplying the total number of pixels in the discrete region by the physical size represented by a single pixel (determined by the product of the channel spacing and the depth sampling interval).

[0112] On each medium gradient indicator map, the 80th percentile of all non-zero pixel values ​​is taken as the anomaly threshold. A connected component labeling algorithm is applied to aggregate all pixels greater than the anomaly threshold that are spatially connected to form one or more discrete regions. Each discrete region is defined as a "candidate property-type anomaly region," and the corresponding geographic bounding box, region area, maximum value of feature values ​​within the region, and average value of feature values ​​are calculated.

[0113] Considering that for diffuse anomalies, single-period data often fails to form clear boundaries, and their subtle evolution in multi-period detection is more diagnostically significant, it is necessary to rely on cross-period spatiotemporal matching to enhance their identification stability. Therefore, further cross-period association of similar types is carried out: for each newly identified candidate anomaly region (including structural and physical property types), based on the geographic location bounding box, all historical anomaly records of the same type that are geographically overlapping (historical candidate anomaly regions) are retrieved in the spatiotemporal database.

[0114] Specifically, based on their respective geographic location bounding boxes, the intersection-union ratio (IUR) between newly identified candidate anomaly regions and historical anomaly records is calculated. An overlap threshold of 0.5 is set. When the IUR is greater than 0.5, it is considered an overlap, and the newly identified candidate anomaly region is used as the latest observation record in the corresponding historical record. If there is no overlap, a new anomaly record entry is created in the database.

[0115] A unique anomaly record entry is updated or created in the database. Each entry contains its type (structural / physical property) and its time-series attributes of each detection, forming an anomaly detection record.

[0116] Retrieve anomaly detection records that have been detected multiple times:

[0117] When the anomaly detection record corresponds to a structured record, the database is used to extract the area sequence and timestamp sequence of each detected region. The timestamp sequence is a sequence consisting of the timestamps of the start time of the detection task.

[0118] By applying univariate linear regression analysis to the regional area sequence and its corresponding timestamp sequence, the slope of the regression line is calculated as the geometric expansion rate, characterizing the evolution trend of the extent of structural anomalies (such as cavities or voids). For example, a positive rate value indicates that the extent of structural anomalies (such as cavities or voids) is expanding.

[0119] Risk levels are determined based on the latest regional area and geometric expansion rate. As an example, area and expansion rate thresholds are set. When the latest regional area is greater than the area threshold and the geometric expansion rate is less than the expansion rate threshold, the risk level is "Stable - Attention"; when the latest regional area is less than or equal to the area threshold and the geometric expansion rate is greater than or equal to the expansion rate threshold, the risk level is "Developing - Early Warning".

[0120] The area threshold can be set in the range of 0.2–1.0 m² according to the safety control standards of tunnel lining or roadbed structure, and is 0.6 m² in this example; the expansion rate threshold can be set in the range of 0.01–0.05 m² / month, and is 0.03 m² / month in this example.

[0121] When the anomaly detection record corresponds to a physical property type, the database is used to extract the sequence of maximum eigenvalues ​​and timestamps for each detected region. The maximum eigenvalue refers to the maximum value among the morphological differences within the corresponding region.

[0122] By applying univariate linear regression analysis to the sequence of maximum data values ​​and their corresponding timestamp sequences, the slope of the regression line is calculated and used as the characteristic evolution rate to characterize the evolution of material properties within the corresponding region. For example, a positive rate value indicates that the material properties (such as moisture content and porosity) in this region are continuously deteriorating.

[0123] Risk levels are determined based on the latest maximum eigenvalue and eigenvalue evolution rate. As an example, eigenvalue thresholds and evolution rate thresholds are set. When the latest maximum eigenvalue is greater than the eigenvalue threshold and the eigenvalue evolution rate is less than the evolution rate threshold, the risk level is "Stable - Attention". When the latest maximum eigenvalue is less than or equal to the eigenvalue threshold and the eigenvalue evolution rate is greater than or equal to the evolution rate threshold, the risk level is "Developing - Early Warning".

[0124] In this process, the morphological difference cost value is linearly normalized in the corresponding data dimension. The eigenvalue threshold can be set to 0.3–0.6 based on the normalization range of the morphological difference cost value, and is 0.45 in this example. The evolution rate threshold can be set to 0.01–0.05 / month, and is 0.03 / month in this example.

[0125] At this point, a risk level label is attached to the anomaly detection record in the database.

[0126] Finally, based on the risk levels and anomaly detection records, a diagnostic report is generated. Specifically, a tabular diagnostic report is created based on the anomaly detection records in the database that have risk level labels. Each row in the report corresponds to an anomaly detection record, listing its core information: unique ID (naming rules set by the implementer), geographical location, anomaly type (structural / physical property type), current state quantitative indicators (such as area or maximum eigenvalue), historical evolution rate (geometric expansion rate or eigenvalue evolution rate), and the final comprehensive risk level.

[0127] After generating the diagnostic report, the locations of all identified anomaly detection records can be marked with different colors or symbols on a Geographic Information System (GIS) layer. The colors or symbols directly correspond to the risk levels assessed in the diagnostic report. This report and map enable managers to identify high-risk road sections and specific defects, and allocate maintenance resources to the areas most in need. The colors or symbols can be chosen by the implementer and are not limited here.

[0128] It should be noted that, in other embodiments of the present invention, regarding the risk level classification, implementers may further refine it into a multi-level risk system beyond "stable-concern" and "development-early warning," for example, including at least levels such as "safe," "minor deterioration," "significant deterioration," and "serious early warning." The boundaries of different levels can be reset and combined according to threshold parameters. In one embodiment of the present invention, targeting the urban roadbed detection scenario, since geological conditions, equipment parameters, and safety standards vary in different regions, implementers may also adjust the threshold parameters according to actual needs; the present invention does not limit this.

[0129] In summary, to address the technical problem of unstable identification of diffuse anomalies in vehicle-mounted ground-penetrating radar (GPR) data acquisition due to dynamic interference, this invention provides a method for detecting anomalous ground areas with real-time GPR interpretation assistance. This invention first acquires the local anomalous signal matrix to be analyzed and a reference waveform set; then, it extracts reflection feature units, analyzes the correlation between the local waveform data corresponding to the reflection feature units and the corresponding reference waveform set, and obtains waveform feature vectors; further, based on the differences in time and amplitude distribution between adjacent columns of reflection feature units, combined with the differences between waveform feature vectors, it constructs a cost matrix for the corresponding two columns of data and obtains the optimal matching result; further, based on the optimal matching result and the distribution of the corresponding reflection feature units' time points in the local anomalous signal matrix, it constructs a diagnostic feature map; finally, it archives the diagnostic feature map and generates a diagnostic report by comparing the diagnostic feature map of the current data segment with historical records in a spatiotemporal manner. This method performs singular value decomposition on the radar data matrix to remove vehicle-mounted dynamic interference; extracts reflection feature units column by column and completes adjacent event matching based on waveform correlation to form a diagnostic feature map; and then performs spatiotemporal comparison to generate a diagnostic report, thereby stably identifying diffuse anomalies.

[0130] It should be noted that the order of the above embodiments of the present invention is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0131] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

Claims

1. A method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation, characterized in that, The method includes: A data matrix of a segment of data collected by ground penetrating radar along the survey line direction is obtained, and singular value decomposition is performed on the data matrix to extract the local anomaly signal matrix and reference waveform set; Reflection feature units are extracted based on the data fluctuations in each column of the local anomaly signal matrix; the correlation between the local waveform data corresponding to the reflection feature unit and the reference waveform set extracted from the previous data segment is analyzed to obtain the waveform feature vector of each reflection feature unit in the current data segment; based on the differences in time and amplitude distribution between the reflection feature units in two adjacent columns, and combined with the differences between the waveform feature vectors, a cost matrix corresponding to the two columns of data is constructed and the optimal matching result is obtained; based on the optimal matching result, and combined with the distribution of the time points of the corresponding reflection feature units in the local anomaly signal matrix, a diagnostic feature map is constructed. The diagnostic feature maps are archived and stored in a database. A diagnostic report is generated based on a spatiotemporal comparison of the diagnostic feature maps of the current data segment with historical records. The method for obtaining the cost matrix includes: The process involves selecting one reflection feature unit from each of two adjacent columns to form a binary tuple; obtaining the time-amplitude difference cost based on the difference between the time sampling points and the difference between the amplitude values ​​of the two elements within the binary tuple; obtaining the morphological difference cost based on the difference between the waveform feature vectors of the two elements within the binary tuple; fusing the time-amplitude difference cost and the morphological difference cost to obtain the total cost; and constructing a cost matrix from the total cost of all binary tuples corresponding to two adjacent columns.

2. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 1, characterized in that, The method for obtaining the local anomaly signal matrix includes: Based on the first left singular vector, the first right singular vector, and the first singular value obtained from the singular value decomposition, a background matrix is ​​reconstructed; based on the background matrix, the background is removed from the data matrix to obtain a local anomaly signal matrix.

3. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 1, characterized in that, The method for obtaining the reference waveform set includes: In the left singular vector matrix obtained by the singular value decomposition, a preset number of column parameters are selected starting from the second column to form a reference waveform set.

4. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 1, characterized in that, The method for obtaining the reflection feature unit includes: For each column of data in the local abnormal signal matrix, a one-dimensional peak detection algorithm is used to obtain the peak, and effective peaks are selected based on a preset amplitude threshold and a preset peak spacing threshold. The time sampling point of the peak value of an effective peak, the normalized amplitude value of the peak value, the index number of the column in the local abnormal signal matrix, and the column number are used to form a reflection feature unit.

5. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 4, characterized in that, The method for obtaining the waveform feature vector includes: For each of the reflection feature units, with the time sampling point corresponding to the reflection feature unit as the center, data of a preset waveform segment length is extracted from the corresponding column of the local abnormal signal matrix as the local waveform data corresponding to the reflection feature unit. The inner product is calculated between the local waveform data and each vector in the reference waveform set extracted from the previous segment of data, and the inner product values ​​are arranged to form the waveform feature vector.

6. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 1, characterized in that, The method for obtaining the diagnostic feature map includes: Based on the local anomaly signal matrix, construct a structural boundary map and a medium gradient indicator map of consistent size; Based on the optimal matching result, according to the attribute of the unmatched reflection feature unit as the starting point or ending point of the path, the starting point mark or ending point mark is filled in at the corresponding position of the time sampling point to generate a structural boundary map. By filling the morphological difference value of each successfully matched pair of reflection feature units into the corresponding position of the time sampling point of the reflection feature unit with the smallest temporal sequence, a medium gradient indicator map is generated.

7. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 1, characterized in that, The method for archiving and building a library of the diagnostic feature maps includes: The diagnostic feature map, the geographic coordinates of the starting and ending points of the currently analyzed data segment, and the timestamp of the start time of the detection task are combined into a data record unit and stored in a structured spatiotemporal database.

8. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 1, characterized in that, Methods for generating diagnostic reports include: Retrieve all feature map data generated by the latest detection task and extract candidate abnormal regions from each feature map; The candidate abnormal regions of the same type are associated across periods, and abnormal record entries are updated or created in the database. Anomaly detection records with multiple detected attributes are retrieved. When the anomaly detection record corresponds to a structural type, univariate linear regression analysis is applied to the region area sequence and its corresponding timestamp sequence to calculate the slope of the regression line, which is used as the geometric expansion rate. The risk level is then determined based on the latest region area and geometric expansion rate. When the anomaly detection record corresponds to a physical property type, univariate linear regression analysis is applied to the maximum data value sequence and its corresponding timestamp sequence to calculate the slope of the regression line, which is used as the feature evolution rate. The risk level is then determined based on the latest maximum feature value and feature evolution rate. A diagnostic report is generated based on the risk level classification and anomaly detection records.

9. The method for detecting abnormal ground areas with the assistance of real-time ground-penetrating radar interpretation according to claim 1, characterized in that, The data matrix has undergone preprocessing, including the following steps: The collected data is segmented into a data matrix by sliding window along the survey line. DC drift correction is performed on each column of data, and the time shift of surface reflection is calculated based on cross-correlation to perform surface reflection alignment.

Images