Aviation ground penetrating radar data acquisition and inversion method based on compressed sensing technology

Abstract

The invention discloses an aviation ground penetrating radar data acquisition and inversion method based on a compressed sensing technology, and the method comprises the steps: reducing the data acquisition amount through sub-sampling, reconstructing complete radar data through a theory and actual measurement combined sparse dictionary, and carrying out the gradient optimization iteration inversion of underground medium electromagnetic parameters based on electromagnetic wave propagation equation forward modeling and an adjoint state method. And high-precision underground structure detection under the condition of low sampling rate is realized.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a method for acquiring and inverting airborne ground-penetrating radar data based on compressed sensing technology. Background Technology

[0002] Unmanned Aerial Vehicle Ground Penetrating Radar (UAV-GPR) is a geophysical exploration technology that uses a ground penetrating radar system mounted on an unmanned aerial vehicle (UAV) platform to conduct non-contact detection of shallow underground media structures.

[0003] Existing UAV-GPR systems typically acquire radar echo data along a predetermined flight path using high-density, equally spaced continuous sampling to meet the requirements of 3D imaging and Full Waveform Inversion (FWI) for data integrity and spatial sampling density. Under 3D survey conditions, the scale of acquired data increases rapidly with the number of survey lines and time sampling points, exhibiting significant high-dimensional characteristics. Although high-density sampling helps improve inversion accuracy, the following technical problems still exist in practical applications: The large number of measurement points limits flight efficiency. High-density sampling leads to a significant increase in the number of measurement points in a single mission. The data acquisition, caching, writing, and transmission processes become one of the major sources of energy consumption for the UAV platform, increasing flight time and operating costs. At the same time, it increases the load on the onboard storage system and communication modules, thereby limiting the coverage and duration of a single operation.

[0004] The data volume is enormous, resulting in significant computational and storage pressure. Full-waveform radar data demands high storage capacity and post-processing computing resources, and the strong correlation between adjacent measurement points leads to a large amount of redundant information that contributes only limited value to inversion constraints, further exacerbating the burden of data management and inversion calculations.

[0005] Inadequate adaptability to complex terrain. In high-risk areas or areas with significant topographic relief, such as landslide-prone and avalanche-prone areas, the high-density, regularly deployed measurement points are limited by flight safety, operational conditions, and time windows, making it difficult to balance detection efficiency and operational safety.

[0006] Therefore, there is an urgent need for an airborne ground-penetrating radar data acquisition method that can still meet the full waveform inversion accuracy requirements under low sampling rate conditions, so as to improve the operational efficiency and engineering applicability of the UAV-GPR system. Summary of the Invention

[0007] This invention proposes a method for acquiring and inverting airborne ground-penetrating radar data based on compressed sensing technology, including: Based on the UAV aerial ground-penetrating radar survey lines planned for the area to be measured, a complete set of measurement points is determined. Based on the sub-sampling rate, measurement points are selected from the complete set of measurement points to form a sub-sampling measurement point set, and undersampled radar data is collected. The theoretical radar response signal is obtained by forward modeling the electromagnetic wave propagation equation using underground medium parameters and then normalized to form a theoretical prior sparse dictionary; the measured radar echo signal samples are decomposed to construct a measured feature sparse dictionary; the theoretical prior sparse dictionary and the measured feature sparse dictionary are combined to form a joint sparse dictionary. Establish a sparse representation relationship between the undersampled radar data, the joint sparse dictionary, and the sparse coefficient vector. Reconstruct the complete radar data by solving the sparse coefficient vector and multiplying it with the joint sparse dictionary. The spatial distribution of the relative permittivity and conductivity of the subsurface medium is defined, and simulated radar data is obtained by forward modeling using the electromagnetic wave propagation equation. A residual objective function is constructed between the complete radar data and the simulated radar data. The gradient of the residual objective function with respect to the relative permittivity and conductivity is calculated and its spatial distribution is updated. The optimization is iterated until the residual objective function converges, and the final spatial distribution of the relative permittivity and conductivity of the subsurface medium is obtained.

[0008] Optionally, selecting measurement points from the complete set of measurement points based on the sub-sampling rate to form a sub-sampling measurement point set includes: Based on the prior information of the topographic relief and underground structure of the area to be measured, a sub-sampling rate for spatial variation is set, and measuring points are selected from the complete set of measuring points according to the sub-sampling rate to form the sub-sampling measuring point set.

[0009] Optionally, the theoretical radar response signal is obtained through forward modeling of the electromagnetic wave propagation equation using underground medium parameters, and then normalized to construct a theoretical prior sparse dictionary, including: Multiple combinations of relative permittivity and conductivity parameters of underground media are set; for each combination of parameters, the corresponding theoretical radar response signal is obtained by forward modeling of the electromagnetic wave propagation equation; the amplitude of each theoretical radar response signal is normalized; the multiple normalized theoretical radar response signals are arranged in columns to form a theoretical prior sparse dictionary.

[0010] Optionally, feature decomposition is performed on the measured radar echo signal samples to construct a sparse dictionary of measured features, including: Multiple sets of measured radar echo signal samples are selected and arranged in columns to form a measured data matrix; the measured data matrix is ​​subjected to feature decomposition to extract basis vectors; the extracted basis vectors are arranged in columns to form a measured feature sparse dictionary.

[0011] Optionally, complete radar data is reconstructed by solving the sparse coefficient vector and multiplying it by the joint sparse dictionary, including: The complete radar data to be reconstructed is represented as the product of the joint sparse dictionary and the sparse coefficient vector; a reconstruction objective function is constructed, which includes an error term between the undersampled radar data and the reconstructed data and a sparse constraint term of the sparse coefficient vector; the sparse coefficient vector is solved by a sparse optimization algorithm to minimize the reconstruction objective function; the complete radar data is obtained by multiplying the solved sparse coefficient vector with the joint sparse dictionary.

[0012] Optionally, the spatial distribution of the relative permittivity and conductivity of the underground medium is defined, and simulated radar data is obtained through forward modeling of the electromagnetic wave propagation equation, including: The spatial distribution of the relative permittivity and conductivity of the underground medium is initialized based on the arrival time and amplitude characteristics of the reflected waves in the complete radar data. The computational region is divided into multiple grid cells, and the spatial distribution is discretized and assigned to each grid cell. The propagation equation of electromagnetic waves in the grid cells is established. The electromagnetic field boundary conditions at the grid cell interface are set according to the relative permittivity and conductivity values ​​of adjacent grid cells. The time-domain waveform of the excitation source is set and substituted into the propagation equation. The electric field components of each grid cell are obtained by discretizing and solving the propagation equation and boundary conditions. The time series of electric field components corresponding to the grid cells at each measurement point in the complete set of measurement points are extracted. The extracted electric field component time series are arranged in the order of measurement points to form simulated radar data.

[0013] Optionally, a residual objective function is constructed between the complete radar data and the simulated radar data. The gradient of the residual objective function with respect to the relative permittivity and conductivity is calculated and its spatial distribution is updated. The optimization is iterated until the residual objective function converges, obtaining the final spatial distribution of the relative permittivity and conductivity of the subsurface medium, including: Calculate the difference between the complete radar data and the simulated radar data at each measurement point location, and sum the squares to construct the residual objective function; Establish the adjoint wave field equation, substitute the difference between the locations of each measuring point as the adjoint source term into the adjoint wave field equation, and obtain the distribution of the adjoint wave field in each grid cell by solving the adjoint wave field equation. The electric field components of each grid cell obtained during the electromagnetic wave forward modeling process are cross-correlated with the distribution of the accompanying wave field in the corresponding grid cell in the time domain to obtain the gradient of the residual objective function with respect to the relative permittivity and conductivity of each grid cell. The update direction of the relative permittivity and conductivity is determined according to the gradient, and the relative permittivity and conductivity values ​​of each grid cell are adjusted according to the update direction. Substitute the updated relative permittivity and conductivity spatial distribution into the electromagnetic wave propagation equation and perform forward modeling again to obtain new simulated radar data. Based on the new simulated radar data, recalculate the residual objective function. When the residual objective function satisfies the convergence condition, the current spatial distribution of relative permittivity and conductivity is taken as the final spatial distribution.

[0014] The beneficial effects of this invention are as follows: By deeply integrating compressed sensing theory and full-waveform inversion technology, inversion accuracy is maintained while significantly reducing data acquisition density, effectively solving multiple technical bottlenecks faced by existing airborne ground-penetrating radar systems in high-density sampling modes. The subsampling strategy significantly reduces the number of measurement points, directly reducing the energy consumption of data acquisition, storage, and transmission on the UAV platform, extending the endurance of a single operation and expanding the coverage area, thus improving flight efficiency and economy. The combined sparse dictionary, integrating physics-driven theoretical priors with data-driven measured characteristics, enhances the robustness and adaptability of data reconstruction, effectively suppressing information loss caused by undersampling. Full-waveform inversion is based on complete reconstructed data, avoiding inversion distortion caused by insufficient sampling in traditional methods, and ensuring the quantitative accuracy of the spatial distribution of electromagnetic parameters in underground media. The method supports flexible configuration of spatially varying subsampling rates, adapting to operational constraints in complex terrain and high-risk areas, balancing detection efficiency and safety, and significantly improving the practical value of airborne ground-penetrating radar systems in engineering scenarios such as glacier crack detection and landslide investigation. Attached Figure Description

[0015] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is a schematic diagram of the data acquisition and inversion method for airborne ground-penetrating radar based on compressed sensing technology proposed in this invention. Detailed Implementation

[0016] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0017] Combination Figure 1 The flowchart illustrating the data acquisition and inversion method of airborne ground-penetrating radar based on compressed sensing technology is provided below. Based on the UAV aerial ground-penetrating radar survey lines planned for the area to be measured, a complete set of measurement points is determined. Based on the sub-sampling rate, measurement points are selected from the complete set of measurement points to form a sub-sampling measurement point set, and undersampled radar data is collected. The theoretical radar response signal is obtained by forward modeling the electromagnetic wave propagation equation using underground medium parameters and then normalized to form a theoretical prior sparse dictionary; the measured radar echo signal samples are decomposed to construct a measured feature sparse dictionary; the theoretical prior sparse dictionary and the measured feature sparse dictionary are combined to form a joint sparse dictionary. Establish a sparse representation relationship between the undersampled radar data, the joint sparse dictionary, and the sparse coefficient vector. Reconstruct the complete radar data by solving the sparse coefficient vector and multiplying it with the joint sparse dictionary. The spatial distribution of the relative permittivity and conductivity of the subsurface medium is defined, and simulated radar data is obtained by forward modeling using the electromagnetic wave propagation equation. A residual objective function is constructed between the complete radar data and the simulated radar data. The gradient of the residual objective function with respect to the relative permittivity and conductivity is calculated and its spatial distribution is updated. The optimization is iterated until the residual objective function converges, and the final spatial distribution of the relative permittivity and conductivity of the subsurface medium is obtained.

[0018] In this example, an unmanned aerial vehicle (UAV) equipped with an airborne ground-penetrating radar system conducts flight surveys along pre-set survey lines. The flight altitude is set according to the actual surface environmental conditions, typically adjusted within the range of 1.5–30m. When power lines, tall trees, or other surface obstacles are present in the detection area, the flight altitude is appropriately increased to meet obstacle avoidance and flight safety requirements; without compromising flight safety, a lower flight altitude is preferred to improve the coupling effect of electromagnetic waves on the underground medium, thereby enhancing detection depth and resolution. The survey line spacing is set within the range of 0.1–0.5m based on the target size and imaging resolution requirements.

[0019] To reduce the data acquisition density and data volume of airborne ground-penetrating radar, a sub-sampling design is implemented for the measurement points along the flight survey line, while ensuring the accuracy of subsequent inversion. The specific implementation steps are as follows.

[0020] Under standard equidistant data acquisition conditions, let the complete set of measurement points along the flight survey line be: ;in Let represent the spatial position of the i-th measurement point along the measurement line, and N represent the total number of measurement points under complete sampling conditions. The sub-sampling rate is set based on the UAV flight speed, target spatial scale, and subsequent inversion resolution requirements: Then the number of effective measurement points after subsampling is .

[0021] First, assess the target detection resolution requirements. If high-precision 3D imaging or detection of fine structures smaller than 1 meter is needed, a subsampling rate (ρ) of 0.65 to 0.7 is preferred. If the target is a large-scale stratigraphic interface or regional structure at the meter level or above, a ρ value of 0.5 to 0.6 can be selected to improve operational efficiency. Second, consider the matching relationship between the UAV's flight speed and data acquisition frequency. When the flight speed is high (close to 3 m / s), a lower subsampling rate should be selected to avoid data acquisition and storage delays due to excessively dense measurement points. When the flight speed is low (close to 1 m / s), the subsampling rate can be appropriately increased to ensure data integrity. Third, consider the sparsity requirements of the subsequent compressed sensing reconstruction algorithm. If the joint sparse dictionary has a large number of atoms (more than 600 in total) and the data sparsity is good, a lower subsampling rate can be selected. If the dictionary size is small or the data sparsity is average, the subsampling rate should be increased to ensure reconstruction accuracy.

[0022] Subsampling can be done using random subsampling: from the complete set of measurement points M measurement points are randomly selected to form a sub-sampling set: Alternatively, a block-based subsampling method can be used: the flight survey line is divided into several continuous blocks, and a subset of survey points are selected within each block to ensure the spatial uniformity of the subsampling points. Specifically, when the terrain of the detection area is highly undulating or the underground structure undergoes drastic spatial changes, random subsampling is preferred to avoid spatial aliasing effects that may result from regular sampling. When the terrain of the detection area is relatively flat and the underground structure exhibits strong spatial continuity, block-based subsampling is preferred. For complex situations where the target simultaneously possesses a large-scale layered structure and local anomalies, a combined approach can be adopted: block-based subsampling is used in the main sections of the survey line to ensure the continuity of the layered structure, while random subsampling is used in areas where local anomalies are predicted to exist to improve the detection sensitivity of anomalies.

[0023] Based on the selected set of subsampling measurement points, a subsampling operator S is constructed to describe the mapping relationship between complete sampled data and undersampled data. Its expression is: ;in, For airborne ground-penetrating radar data obtained under complete sampling conditions, This represents the undersampled data obtained after subsampling. The subsampling operator S is constructed as follows: S is represented as an M-row, N-column matrix, where M is the number of subsampled measurement points and N is the number of complete measurement points. The element in the i-th row of matrix S is 1 only at the j-th column, and all other elements are 0, where j corresponds to the index position of the i-th selected measurement point in the complete measurement point set. According to the subsampling design, data is collected only at the subsampled measurement point locations during actual flight; or, after completing the complete data collection, the original data is equivalently extracted and stored. The obtained undersampled data serves as input data for subsequent inversion and interpretation processing based on sparse constraints.

[0024] In this embodiment, a three-dimensional FDTD numerical model is constructed based on existing geological and geophysical prior information. According to the spatial distribution characteristics of electromagnetic parameters in different media, the computational domain is discretized into partitions: finer spatial grids are used in regions with high relative permittivity or drastic parameter changes to improve the accuracy of electromagnetic field numerical calculations; coarser spatial grids are used in regions with low relative permittivity or relatively gradual parameter changes to reduce computational scale. Specifically, the air layer region is uniformly discretized using coarse grids, thereby reducing the number of computational units, memory usage, and computational complexity while meeting the requirements of numerical stability and computational accuracy.

[0025] The selection of the spatial grid step size is determined based on the electromagnetic parameters of the medium and the frequency of the excitation signal, and its calculation formula is as follows: Where c is the speed of light, and f is the center frequency of the emitted excitation pulse. The relative permittivity of the corresponding medium is given by N, which is a constant, typically 10. For the air layer region, the relative permittivity is taken as 1. If the radar center frequency is 100MHz, then the spatial step size is... It can be taken as 0.3m; for a relative permittivity of The underground medium, with a spatial step size of 1. For example, underground media If the value is 9, then the grid step size for the underground space region is... The time step can be set to 0.1m. Using the subgrid FDTD method, different grids can be used for the air and subsurface regions, thus saving memory. To ensure numerical stability, the time step is determined according to the Courant stability condition; in the three-dimensional case, this means the time step ≤ the spatial step. ÷(speed of light × √3).

[0026] A sparse dictionary is constructed to sparsely represent airborne ground-penetrating radar data, thereby improving the accuracy of data reconstruction and the stability of inversion under undersampling conditions. The sparse dictionary is constructed using a combination of theoretical model-driven and measured data feature extraction methods.

[0027] A theoretical model based on underground medium parameters is used to construct theoretical response basis functions using forward modeling operators. The combination of theoretical parameters is designed according to the geological background and prior information of the area to be measured. For snow and ice media, the relative permittivity is typically taken in the range of 2 to 4, and the electrical conductivity in the range of 0.001 to 0.01 S / m; for soil media, the relative permittivity is taken in the range of 4 to 20, and the electrical conductivity in the range of 0.01 to 0.1 S / m. Specific sampling methods and quantity settings for parameter combinations can be found in other embodiments of this patent. The i-th theoretical response vector is represented as... ,in is the forward operator, Let be the vector of the i-th group of subsurface medium parameters. The forward modeling operator is implemented using the aforementioned three-dimensional subgrid FDTD method. A corresponding numerical model is established for each parameter combination. The source is set as a Gaussian modulated pulse or a Ricker wavelet, and the receiver is set at a fixed distance above the source, typically 0.5 to 2 meters. The time series of the electric field components is recorded as the theoretical radar response signal. After amplitude normalization of the theoretical response vector, a theoretical prior sparse dictionary is constructed. The specific steps of the normalization process are as follows: for each theoretical response vector, firstly, the DC component is removed, i.e., the average value of all elements of the vector is subtracted; then, the maximum absolute value of the vector is calculated, and each element of the vector is divided by this maximum absolute value, so that the amplitude of the normalized vector is between -1 and +1. A theoretical prior sparse dictionary is then established. .

[0028] From actual airborne ground-penetrating radar data or historical measurement data, priority should be given to selecting measurement point data with high signal-to-noise ratio and clear echo characteristics; ensuring that the samples are spatially representative and cover different terrain conditions and underground structure types; the sample size N2 is determined based on the total number of measurement points along the survey line and data diversity, typically taking 20% ​​to 50% of the complete measurement points, but not less than 100 and not more than 300. Sample signals need to be preprocessed before selection, including removing DC offset, filtering high-frequency noise, and aligning the time zero point. An actual measurement sample set should be established. Each column vector represents a radar echo signal at a typical measurement point or measurement line location.

[0029] We perform feature decomposition on the measured sample set, extract a set of basis vectors that can effectively represent the main structural features of the measured data, and construct a sparse dictionary of measured features. The eigenvalue decomposition employs singular value decomposition, decomposing matrix X into the product of a left singular vector matrix U, a singular value diagonal matrix, and a right singular vector matrix. The column vectors corresponding to the first K largest singular values ​​in matrix U are extracted as basis vectors. The value of K is determined based on the cumulative energy ratio, typically set to ensure that the sum of squares of the first K singular values ​​accounts for 85% to 95% of the total sum of squares of singular values. Specific feature extraction algorithms and parameter settings can be found in other embodiments of this patent. Its objective function can be expressed as follows: ,in, It is the corresponding sparse representation coefficient matrix, and satisfies the preset sparse constraint conditions.

[0030] The theoretical prior sparse dictionary and the measured feature sparse dictionary are combined to form a joint sparse dictionary. The joint dictionary is constructed by concatenating D1 and D2 column by column. To ensure the consistency of the two types of dictionary atoms on the numerical scale, D1 and D2 are column normalized before concatenation to ensure that the L2 norm of each atomic vector is equal to 1. .

[0031] The complete data vector x to be reconstructed is represented as a linear combination of atomic vectors in the joint sparse dictionary, i.e. ,in, This represents a sparse coefficient vector, where the number of non-zero elements is less than a preset threshold. Under undersampling conditions, the corresponding observed data is represented as follows: ,in, This represents the data acquisition mapping operator. This is the noise term. The relationship between the data acquisition mapping operator G and the aforementioned subsampling operator S is that G equals S. The noise term n is assumed to be additive white Gaussian noise with a mean of 0 and a variance determined based on the signal-to-noise ratio estimation of the measured data. The sparsity threshold is typically set to 5% to 15% of the total number of atoms in the joint dictionary.

[0032] Different weight constraints are applied to the components of the sparse coefficient vector corresponding to the theoretical prior dictionary and the measured feature dictionary, respectively, to balance the influence of the two types of information on the data reconstruction result. The objective function can be summarized as follows: ,in, and Representing the theoretical prior sparse dictionary respectively sparse dictionary of measured features The corresponding sparse coefficient subvector, and Here, λ1 and λ2 are the weighting coefficients. The sparse constraint function is defined as the L1 norm of the coefficient vector. The weighting coefficients λ1 and λ2 are determined as follows: the optimal combination is selected from a set of candidate values ​​using cross-validation, with the candidate values ​​typically ranging from 0.001 to 0.1. In the absence of validation data, empirical settings can be used. When the theoretical prior information is relatively reliable, λ1 is set to 0.01 to 0.02 and λ2 to 0.02 to 0.05; when the measured data quality is good and the sample is sufficient, λ1 is set to 0.02 to 0.05 and λ2 to 0.01 to 0.02.

[0033] After completing the construction of the joint sparse dictionary, the undersampled acquisition data is reconstructed using the above sparse representation relationship to recover the complete full waveform data. The specific implementation method is as follows: (1) Recovery of complete full waveform data. Using the obtained sparse coefficient vector, according to Reconstructing the complete data, including To recover the complete full waveform data. (2) Constraints and output of the reconstruction results. In some embodiments, amplitude constraints or energy constraints can be applied to the reconstruction results to suppress noise amplification effects, ensure the rationality of the reconstructed waveform in a physical sense, and finally output for subsequent inversion or imaging processing.

[0034] After completing the construction of the joint sparse dictionary and establishing the sparse representation relationship, the airborne ground-penetrating radar data acquired under undersampling conditions is reconstructed based on the sparse representation model to recover the complete full waveform data.

[0035] Solve the corresponding sparse coefficient vector based on the undersampled observation data. The solution method employs the Fast Iterative Threshold Soft Recovery (FISTA) algorithm or other sparse optimization algorithms. The specific implementation steps of the FISTA algorithm are as follows: initialize the sparse coefficient vector and auxiliary vector to zero, and initialize the momentum parameter to 1; in each iteration, first calculate the gradient direction, then perform gradient descent updates, with the step size determined by the reciprocal of the largest eigenvalue of matrix GD; next, perform a soft threshold shrinkage operation, setting the soft threshold to the step size multiplied by the weight coefficient; update the auxiliary vector and momentum parameter; terminate the iteration when the relative change of the objective function is less than a preset threshold (e.g., 0.0001) or the number of iterations reaches an upper limit (e.g., 200 times). Detailed implementation of the sparse optimization algorithm can be found in the descriptions of other embodiments of this patent. Furthermore, a joint sparse dictionary is used to reconstruct the complete data, the expression of which is: ,in, This represents the complete waveform data obtained from the reconstruction.

[0036] To suppress noise amplification or non-physical oscillations that may occur under undersampling conditions, amplitude constraints, energy constraints, or smoothing constraints are applied to the reconstruction results to ensure the physical rationality of the reconstructed waveform in time and space. The amplitude constraint is implemented by determining reasonable upper and lower limits based on the amplitude statistical characteristics of the measured data, typically set to 0.8 to 1.2 times the maximum absolute value of the measured data. Elements in the reconstructed data exceeding this range are truncated. Specific implementation methods for energy constraints and smoothing constraints can be found in other embodiments of this patent. The constrained reconstructed data is then used as input data for subsequent full waveform inversion or imaging processing.

[0037] Based on the reconstructed complete waveform data, the electromagnetic parameters of the underground medium are iteratively solved using the full waveform inversion method. The spatial distribution of the relative permittivity and conductivity of the underground medium is obtained through inversion, and a corresponding three-dimensional visualization model is generated. The specific implementation steps are as follows.

[0038] First, based on the subgrid finite difference time-domain (FDTD) method, the control equations describing the propagation of electromagnetic waves from airborne ground-penetrating radar are numerically solved, and the model parameters in the k-th iteration are obtained. Under the given conditions, calculate the corresponding forward-modeling simulated radar data. The forward modeling calculation employs the same three-dimensional subgrid FDTD method as the theoretical prior dictionary construction method. A numerical model is established based on the spatial distribution of relative permittivity and conductivity in the current iteration. A transmitter is set up at each measurement point and the electric field response at the receiver point is calculated.

[0039] Secondly, construct the objective function. Used to characterize observation data With forward simulation data The difference between them is expressed as follows: , in, and These represent the locations of the radar transmitter and receiver, respectively. Represents a time variable.

[0040] To improve the efficiency of inversion calculation, the adjoint state method is used to calculate the gradient of the objective function with respect to the model parameters. Based on this, the model parameters are iteratively updated using an optimization algorithm. The optimization algorithm can be the conjugate gradient method, the L-BFGS method, or other quasi-Newton algorithms. The model update formula for the next iteration is: ,in, The step size factor is determined using a line search method, employing the Armijo criterion. If the condition is not met, the step size is reduced to 0.5 times the original value, and the test is repeated until the condition is met or the step size is less than the minimum threshold. Detailed implementation of the line search can be found in the descriptions of other embodiments of this patent. Let be the quasi-Hessian matrix or its approximate form constructed by the optimization algorithm.

[0041] Repeat the above forward modeling, gradient calculation, and model update process until the objective function is achieved. The process converges to a preset threshold or the number of iterations reaches a predetermined upper limit, ultimately yielding the inversion results of the electromagnetic parameters of the subsurface medium. Based on this, a three-dimensional inversion model of the underground structure is generated for subsequent interpretation and analysis.

[0042] After completing the full waveform data reconstruction and subsequent full waveform inversion processing based on compressed sensing, the effectiveness of the method in underground structural exploration is verified by constructing a numerical verification model, and the corresponding geological interpretation results are obtained.

[0043] The high-precision relative permittivity distribution of the underground medium obtained from the full waveform inversion is used as a verification reference model, denoted as... UAV-GPR forward modeling data simulation and measurement point subsampling: Based on the reference model, numerical forward modeling methods are used to generate airborne ground-penetrating radar echo data under complete measurement point conditions. .

[0044] Based on this, a subsampling operator S is introduced to perform subsampling processing on the complete measurement point data to obtain undersampled observation data. This was used to simulate the undersampling data acquisition process under actual UAV-GPR operating conditions.

[0045] For undersampled observation data, the aforementioned compressed sensing reconstruction method is used to reconstruct the full waveform data, thus obtaining complete waveform data. ,in, This represents the compressed sensing reconstruction operator. Further, using the reconstructed waveform data as input, full waveform inversion processing is performed to obtain the inversion results of the subsurface medium parameters. ,in, This represents the full waveform inversion operator.

[0046] Comparison, evaluation, and validity determination of inversion results: The subsurface medium parameter model obtained by inversion... Compared with the reference model Comparative analysis is conducted, and the inversion accuracy is evaluated using error indices or structural consistency indices. The error can be expressed as... . Simultaneously, the consistency of subsurface interface location, stratum depth, and anomaly spatial distribution can be combined to perform qualitative and quantitative comprehensive verification of the inversion results. Quantitative evaluation indicators also include: root mean square error of relative permittivity; and interface location error, defined as the difference between the depth of the inverted stratigraphic interface and the interface depth of the reference model. Qualitative evaluation involves comparing the inversion results with two-dimensional cross-sectional images of the reference model to observe the consistency of major stratigraphic structures, interface morphology, and anomaly locations. When the error indicators are below preset thresholds, the relative error E is less than 0.1, the root mean square error is less than 10% of the range of dielectric constant variation in the reference model, the average interface location error is less than one grid step, and the subsurface structural features obtained through inversion maintain a high degree of consistency with the reference model in spatial distribution and geometric morphology, the judgment method demonstrates good application performance in subsurface structural detection.

[0047] This invention combines compressed sensing reconstruction and full waveform inversion to maintain inversion accuracy while increasing data acquisition volume, thus significantly improving the operational efficiency of airborne ground-penetrating radar.

[0048] Optionally, selecting measurement points from the complete set of measurement points based on the sub-sampling rate to form a sub-sampling measurement point set includes: One of the following methods can be used: selecting measurement points from the complete measurement point set according to a pseudo-random sequence to form the sub-sampling measurement point set; dividing the flight survey line into multiple continuous blocks, and selecting measurement points in each block according to the sub-sampling rate to form the sub-sampling measurement point set; setting a spatial variation sub-sampling rate based on the prior information of the terrain undulation and underground structure of the area to be measured, and selecting measurement points from the complete measurement point set according to the spatial variation sub-sampling rate to form the sub-sampling measurement point set.

[0049] In this embodiment, when selecting measurement points from the complete set of measurement points based on the sub-sampling rate to form a sub-sampling measurement point set, any one of the following methods can be used: pseudo-random sequence selection, block selection, or spatially varied sub-sampling rate selection.

[0050] When using a pseudo-random sequence selection method, a linear congruence generator or Mason tween algorithm is used to generate a pseudo-random number sequence. The seed value of the pseudo-random number sequence is set according to the survey line number or timestamp to ensure the difference in subsampling patterns between different survey lines. Each element in the pseudo-random number sequence is modulo the total number of complete survey points N to obtain the corresponding survey point index. When the survey point corresponding to this index has not yet been selected, it is added to the subsampling survey point set. This process is repeated until the number of survey points in the subsampling survey point set reaches M. To avoid local over-density or over-sparseness caused by pseudo-random selection, a minimum survey point spacing constraint is set. When the minimum distance between a newly selected survey point and an already selected survey point is less than a set threshold, the survey point is discarded and the next pseudo-random index is generated. The minimum survey point spacing threshold is usually set to 0.3 to 0.5 times the complete sampling interval.

[0051] When using a block-based selection method, the flight survey line is divided into multiple continuous blocks of equal length along the survey line direction. The number of blocks is set to 5 to 20, and the length of each block is equal to the total length of the survey line divided by the number of blocks. Within each block, the number of measurement points to be selected is determined according to the sub-sampling rate ρ. The number of measurement points in a block is equal to the total number of complete measurement points contained in the block multiplied by the sub-sampling rate ρ and rounded down. Measurement points are selected within each block at uniform intervals, with the uniform interval distance equal to the block length divided by the number of measurement points to be selected within the block. When there is a deviation between the location of a complete measurement point within a block and the theoretical uniform selection location, the actual measurement point closest to the theoretical location is selected and added to the sub-sampling measurement point set. The block-based selection method ensures the uniform spatial distribution of sub-sampling measurement points throughout the entire survey line, avoiding spatial clustering that may occur due to random selection.

[0052] When using the spatial variation subsampling rate selection method, different subsampling rates are set for different spatial locations based on the topographic relief and prior information on the underground structure of the area to be measured. In areas with dramatic topographic relief or complex underground structures, a higher subsampling rate is set, ranging from 0.7 to 0.9, to ensure data acquisition density in that area. In areas with flat topography or simple underground structures, a lower subsampling rate is set, ranging from 0.4 to 0.6, to reduce the amount of data collected. The degree of topographic relief is quantified by extracting the rate of elevation change along the survey line using a digital elevation model. Areas with dramatic topographic relief are identified when the elevation difference between adjacent measuring points exceeds 2 meters. Prior information on underground structures comes from existing geological data or pre-scanning results of the exploration area. Areas with complex structures are identified when the number of reflected waves exceeds 5 or the rate of change of reflected wave amplitude exceeds 30% based on the pre-scanning radar data. The survey line is divided into several spatial segments, each with a length of 5 to 10 meters. A corresponding subsampling rate is assigned to each spatial segment based on its topographic relief and structural complexity. Within each spatial segment, measurement points are selected from the complete measurement points contained in that segment according to the sub-sampling rate of that segment. The selection method can be uniform interval or pseudo-random method. The selected measurement points are summarized to form the final sub-sampling measurement point set.

[0053] This invention optimizes the distribution of measurement points according to actual needs by flexibly combining various subsampling selection methods, while ensuring the accuracy of data reconstruction, thereby improving data acquisition efficiency.

[0054] Optionally, the theoretical radar response signal is obtained through forward modeling of the electromagnetic wave propagation equation using underground medium parameters, and then normalized to construct a theoretical prior sparse dictionary, including: Multiple combinations of relative permittivity and conductivity parameters of underground media are set; for each combination of parameters, the corresponding theoretical radar response signal is obtained by forward modeling of the electromagnetic wave propagation equation; the amplitude of each theoretical radar response signal is normalized; the multiple normalized theoretical radar response signals are arranged in columns to form a theoretical prior sparse dictionary.

[0055] In this embodiment, when constructing the theoretical prior sparse dictionary, multiple different combinations of relative permittivity and conductivity parameters for underground media are set as inputs to the theoretical model. The range of relative permittivity values ​​is determined according to the geological type of the area to be measured: 2 to 4 for snow and ice media, 4 to 20 for soil media, and 5 to 15 for rock media. Within the set range, the relative permittivity is discretized and sampled using equal or logarithmic intervals. The interval step size for equal intervals is set to 0.1 to 0.5, while the logarithmic interval uses denser intervals in low-value areas and sparser intervals in high-value areas. The conductivity value range is set to 0.001 S / m to 0.1 S / m, and within this range, it is sampled using logarithmic intervals. The base of the logarithmic interval is set to 1.5 to 3, meaning each value is equal to the previous value multiplied by the base. Different relative permittivity values ​​are combined with different conductivity values ​​by Cartesian product to generate multiple sets of parameter combinations. The total number of parameter combinations is usually set to 200 to 500 sets to cover the range of possible changes in the electromagnetic parameters of the medium in the area under test.

[0056] For each parameter combination, the corresponding theoretical radar response signal is obtained through forward modeling of the electromagnetic wave propagation equation. The forward modeling calculation solves Maxwell's equations using the aforementioned finite-difference time-domain method, establishing a simplified layered medium model or a medium model containing a single target body. The layered medium model includes an air layer and a homogeneous subsurface medium layer, while the target body model embeds a spherical or cylindrical anomaly within the homogeneous medium. The relative permittivity and conductivity of the current parameter combination are assigned to the subsurface medium layer or anomaly region, with the relative permittivity of the air layer fixed at 1 and the conductivity fixed at 0. The radar transmitter and receiver positions are set, with the horizontal distance between the transmitter and receiver set to 0.2 to 1 meter and the vertical height set to 1 to 5 meters above ground level. The transmitter excitation signal uses a Ricker wavelet or Gaussian derivative pulse, with the center frequency consistent with the actual radar system. The calculation is performed step-by-step using the finite-difference time-domain method, extracting the time series of the electric field component at the receiver position as the theoretical radar response signal corresponding to that parameter combination. The number of time sampling points for the theoretical radar response signal is consistent with the single-channel recording duration and sampling rate of the actual radar system.

[0057] Each theoretical radar response signal undergoes amplitude normalization. The normalization method involves iterating through all time sampling points of the signal, finding the maximum absolute value of the signal amplitude, and dividing the value of each sampling point by this maximum absolute value. The normalized signal amplitude ranges from -1 to +1, eliminating amplitude differences between signals corresponding to different parameter combinations while preserving the signal's waveform characteristics and phase information. When the maximum absolute value of the theoretical radar response signal is less than a preset noise threshold, the signal corresponding to that parameter combination is deemed too weak and discarded from the dictionary. The noise threshold is typically set to 10 to 100 times the numerical accuracy of the forward modeling calculation.

[0058] Multiple normalized theoretical radar response signals are arranged column-wise to form a theoretical prior sparse dictionary D1. Each theoretical radar response signal is a column of the dictionary, with the number of time sampling points corresponding to the number of rows and the number of parameter combinations corresponding to the number of columns. Dictionary D1 is stored as a two-dimensional array or matrix data structure, using double-precision floating-point numbers to ensure numerical accuracy. After the dictionary is constructed, the 2-norm of each column is calculated. If the 2-norm of a column is less than 0.1, that column is considered redundant or invalid and is deleted from the dictionary. The retained column vectors are orthogonalized using Gram-Schmidt orthogonalization or QR decomposition to reduce the cross-correlation coefficient between different column vectors. The threshold for the cross-correlation coefficient is set to 0.95; if the cross-correlation coefficient between two column vectors exceeds this threshold, one column is retained and the other is deleted.

[0059] This invention constructs a theoretical prior dictionary covering the medium parameter space through systematic parameter combination design and forward modeling, providing a physical constraint basis for subsequent data reconstruction.

[0060] Optionally, feature decomposition is performed on the measured radar echo signal samples to construct a sparse dictionary of measured features, including: Multiple sets of measured radar echo signal samples are selected and arranged in columns to form a measured data matrix; the measured data matrix is ​​subjected to feature decomposition to extract basis vectors; the extracted basis vectors are arranged in columns to form a measured feature sparse dictionary.

[0061] In this example, when constructing a sparse dictionary of measured features, multiple sets of measured radar echo signal samples are selected as input data. The sources of these samples include pre-acquisition data of the area to be measured, historical measurement data of adjacent areas, or existing detection data under similar geological conditions. Sample selection follows the principle of representativeness, covering radar echo signals from different measurement points, different terrain conditions, and different underground reflection characteristics. The sample size N² needs to balance dictionary construction quality and computational complexity, and is typically set within the range of 100 to 300. The upper limit of the sample size range is used when the geological conditions of the area to be measured are complex or the underground structure changes drastically; the lower limit is used when the geological conditions are simple or the structure is uniform. The number of time sampling points for the sample signals is consistent with the single-channel recording duration and sampling rate of the actual radar system, ensuring that the samples and the data to be reconstructed have the same time dimension. Before selecting the sample signals, they need to be preprocessed, including removing DC offset, filtering high-frequency noise, and aligning the time zero point. The DC offset is achieved by calculating the signal mean and subtracting the mean from each sampling point. High-frequency noise is filtered out by using a low-pass filter with a cutoff frequency set to 2 to 3 times the radar center frequency. The time zero point is aligned by detecting the arrival time of the direct wave and aligning the peak values ​​of the direct waves of all signals to the same time position.

[0062] The selected N² measured radar echo signal samples are arranged column-wise to form a measured data matrix X. Each radar echo signal is a column of the matrix, the number of time sampling points of the signal corresponds to the number of rows L of the matrix, and the number of samples corresponds to the number of columns N². Matrix X has L rows and N² columns. Each element of matrix X represents the amplitude of a sample signal in a certain column at a certain time sampling point in a certain row. After the matrix is ​​constructed, it is centered by calculating the mean of each row and subtracting the mean of that row from all elements. This centering process eliminates amplitude deviations between different samples and highlights the signal variation characteristics.

[0063] The measured data matrix X is subjected to eigenvalue decomposition to extract basis vectors. The eigenvalue decomposition method is singular value decomposition (SVD). Following the aforementioned SVD method, matrix X is decomposed into the product of a left singular vector matrix U, a singular value diagonal matrix S, and the transpose VT of the right singular vector matrix V. The column vectors of the left singular vector matrix U are the extracted basis vectors. Each basis vector corresponds to a singular value, and the magnitude of the singular value reflects the contribution of that basis vector to the original data matrix. When extracting basis vectors, the singular values ​​are arranged in descending order, and the basis vectors corresponding to the first K largest singular values ​​are selected. The value of K is determined using the cumulative energy proportion criterion. The sum of squares of the first K singular values ​​is calculated, and this sum is divided by the sum of squares of all singular values ​​to obtain the cumulative energy proportion. Extraction stops when the cumulative energy proportion reaches 95% to 99%. This threshold is adjusted according to the data reconstruction accuracy requirements; it is the upper limit of the threshold range when higher accuracy is required, and the lower limit of the threshold range when lower accuracy is required. The current K basis vectors can represent the main information of the original data matrix; subsequent basis vectors mainly correspond to noise or secondary features.

[0064] The extracted K basis vectors are arranged column-wise to form the measured feature sparse dictionary D2. Each basis vector serves as a column of the dictionary. The dimension of the basis vectors is equal to the number of rows in the dictionary corresponding to the number of time sampling points L, and the number of extracted basis vectors K corresponds to the number of columns in the dictionary. D2 is an L-row, K-column matrix. Since the basis vectors already satisfy the orthogonality condition during extraction (i.e., the dot product between different basis vectors is zero), the column vectors of the measured feature sparse dictionary D2 are orthogonal and require no additional orthogonalization. D2 is stored as a two-dimensional array or matrix data structure, using double-precision floating-point numbers. After the dictionary is constructed, each column is normalized by calculating its 2-norm and dividing all elements of that column by the 2-norm. After normalization, the 2-norm of each column equals 1, ensuring that different atomic vectors in the dictionary have a uniform energy scale.

[0065] This invention enhances the dictionary's adaptability to real data by performing feature decomposition on measured data and extracting basis vectors that can effectively characterize the main structural features of actual radar echo signals.

[0066] Optionally, complete radar data is reconstructed by solving the sparse coefficient vector and multiplying it by the joint sparse dictionary, including: The complete radar data to be reconstructed is represented as the product of the joint sparse dictionary and the sparse coefficient vector; a reconstruction objective function is constructed, which includes an error term between the undersampled radar data and the reconstructed data and a sparse constraint term of the sparse coefficient vector; the sparse coefficient vector is solved by a sparse optimization algorithm to minimize the reconstruction objective function; the complete radar data is obtained by multiplying the solved sparse coefficient vector with the joint sparse dictionary.

[0067] In this embodiment, the complete radar data to be reconstructed is represented as the product of a joint sparse dictionary and a sparse coefficient vector, and a sparse representation relationship is established as described above. Complete radar data contains radar echo signals at all measurement point locations in the complete measurement point set, while the undersampled radar data actually acquired only contains signals at a subset of measurement point locations corresponding to the subsampled measurement point set. The goal of the reconstruction process is to utilize the prior information from the undersampled radar data and the joint sparse dictionary to recover the radar echo signals at the unacquired measurement point locations, thereby obtaining the complete radar data.

[0068] A reconstruction objective function is constructed to achieve data reconstruction. This objective function comprises two main components: a data fitting term and a sparsity constraint term. The data fitting term measures the degree of matching between the undersampled radar data and the reconstructed data at the subsampled measurement point locations. It is expressed as the undersampled radar data minus the result of the subsampling operator acting on the reconstructed data, followed by the squared 2-norm of this difference. A smaller data fitting term indicates that the reconstructed data is closer to the measured data at the acquired measurement point locations. The sparsity constraint term constrains the sparsity of the sparse coefficient vector. It is expressed as the 1-norm of the sparse coefficient vector, i.e., the sum of the absolute values ​​of all elements in the sparse coefficient vector. A smaller sparsity constraint term indicates fewer non-zero elements in the sparse coefficient vector, and the reconstruction result better matches the signal's sparse prior. The reconstruction objective function is a weighted sum of the data fitting term and the sparsity constraint term. The sparsity regularization parameter λ controls the trade-off between the two terms; a larger λ value emphasizes sparsity, while a smaller λ value emphasizes data fitting accuracy.

[0069] The sparse coefficient vector is solved using a sparse optimization algorithm to minimize the reconstruction objective function. Any one of the following sparse optimization algorithms can be used: orthogonal matching pursuit, basis pursuit, iterative soft thresholding, or alternating direction multiplier method. When using the orthogonal matching pursuit algorithm, the algorithm employs an iterative greedy strategy to progressively select the atom vector in the dictionary that is most relevant to the current residual. In each iteration, the inner product of the current residual and all column vectors in the dictionary is calculated, and the column vector with the largest absolute value of the inner product is added to the selected atom set. The values ​​of the sparse coefficient vector on the selected atom are updated using the least squares method, the new residual is calculated, and the next iteration begins. The upper limit of the number of iterations is set to 2 to 3 times the sparsity, which is estimated based on the signal complexity and is typically set to 5% to 20% of the number of columns in the dictionary. The algorithm terminates when the 2-norm of the residual is less than a preset threshold or when the upper limit of the number of iterations is reached. The residual threshold is set to 0.01 to 0.1 times the 2-norm of the undersampled data. When using the basis pursuit algorithm, the algorithm transforms the sparse optimization problem into a linear programming problem. An auxiliary variable is introduced to split each element of the sparse coefficient vector into positive and negative parts. The objective function and constraints of the linear programming problem are then constructed, and the linear programming solver is called to obtain the sparse coefficient vector. The convergence tolerance of the basis pursuit algorithm is set to 1×10⁻⁶. -6 Up to 1×10 -8The maximum number of iterations is set to 1000 to 5000. When using the iterative soft thresholding algorithm, the algorithm updates the sparse coefficient vector iteratively. Each iteration includes a gradient descent step and a soft thresholding shrinking step. The gradient descent step updates the sparse coefficient vector based on the gradient of the data fitting term. The soft thresholding shrinking step performs element-wise soft thresholding on the updated vector elements. If the absolute value of an element is less than the threshold, it is set to zero; otherwise, the element is retained and its absolute value is reduced. The soft threshold is equal to the sparse regularization parameter λ multiplied by the gradient descent step size, which is set to 0.001 to 0.01. The number of iterations is set to 100 to 500. When using the alternating direction multiplier method, the algorithm introduces auxiliary variables and Lagrange multipliers, decomposing the original optimization problem into multiple subproblems that are solved alternately. Each subproblem has a closed-form solution or is easily solved numerically. The convergence speed of this algorithm is faster than the iterative soft thresholding algorithm. The penalty parameter is set to 1 to 10, and the maximum number of iterations is set to 50 to 200.

[0070] The complete radar data is obtained by multiplying the sparse coefficient vector obtained from the solution with the joint sparse dictionary. Matrix multiplication is performed according to standard matrix multiplication rules. Each column of the joint sparse dictionary is multiplied by the corresponding element of the sparse coefficient vector, and the sum is obtained to get the value of the reconstructed data at each time sampling point. The reconstructed complete radar data is an N-row L-column matrix containing the radar echo signals of all N measurement points in the complete measurement point set at L time sampling points. The signals of the reconstructed data at the sub-sampled measurement point locations should be highly consistent with the undersampled radar data. The signals at the unsampled measurement point locations are generated by linear combination interpolation using the sparse dictionary. After reconstruction, the quality of the reconstructed data is evaluated by calculating the root mean square error (RMSE) between the reconstructed data at the sub-sampled measurement point locations and the undersampled radar data. The RMS error should be less than 10% of the amplitude of the undersampled data. If the error exceeds this threshold, the sparse regularization parameter λ is adjusted or the number of dictionary atoms is increased, and reconstruction is performed again.

[0071] This invention efficiently solves the sparse coefficient vector using a sparse optimization algorithm and utilizes the dual prior constraints of physical and data in the joint dictionary to achieve high-precision data reconstruction under undersampling conditions.

[0072] Optionally, the spatial distribution of the relative permittivity and conductivity of the underground medium is defined, and simulated radar data is obtained through forward modeling of the electromagnetic wave propagation equation, including: The spatial distribution of the relative permittivity and conductivity of the underground medium is initialized based on the arrival time and amplitude characteristics of the reflected waves in the complete radar data. The computational region is divided into multiple grid cells, and the spatial distribution is discretized and assigned to each grid cell. The propagation equation of electromagnetic waves in the grid cells is established. The electromagnetic field boundary conditions at the grid cell interface are set according to the relative permittivity and conductivity values ​​of adjacent grid cells. The time-domain waveform of the excitation source is set and substituted into the propagation equation. The electric field components of each grid cell are obtained by discretizing and solving the propagation equation and boundary conditions. The time series of electric field components corresponding to the grid cells at each measurement point in the complete set of measurement points are extracted. The extracted electric field component time series are arranged in the order of measurement points to form simulated radar data.

[0073] In this embodiment, when initializing the spatial distribution of the relative permittivity and conductivity of the underground medium based on the arrival time and amplitude characteristics of reflected waves in the complete radar data, the arrival times of the main reflection events at each measuring point are extracted from the reconstructed complete radar data matrix. The arrival time of reflected waves is picked up using an automatic detection algorithm, traversing the time series of each radar signal. When the signal amplitude exceeds 3 to 5 times the standard deviation of the background noise, it is determined to be a valid reflection event, and this moment is recorded as the arrival time of the reflected wave. For the arrival times of reflected waves at different measuring points on the same reflecting interface, continuous time contour lines are established using time-domain correlation matching or phase tracking methods. The spatial sampling interval of the time contour lines is consistent with the measuring point interval. According to the two-way travel time relationship, the arrival time of reflected waves is equal to the total time it takes for the electromagnetic wave to propagate from the transmitting antenna to the reflecting interface and back to the receiving antenna. This time is related to the depth of the reflecting interface and the electromagnetic wave velocity in the medium. The electromagnetic wave velocity is inversely proportional to the relative permittivity; the velocity is equal to the speed of light divided by the square root of the relative permittivity. Assuming the initial model is a horizontally layered structure, the relative permittivity is estimated layer by layer from the surface downwards. The relative permittivity of the first layer is inferred from the arrival time of the first reflecting interface and the assumed depth. The assumed depth is initially set to 5% to 10% of the survey line length as a rough estimate, and the relative permittivity is calculated using two-way travel time relationships. The relative permittivity of subsequent layers is calculated sequentially based on the time difference between adjacent reflecting interfaces and the layer thickness differences. When no obvious reflecting interface can be identified, the initial relative permittivity is set to a uniform value based on the geological background of the area to be measured: 3 for snow and ice media, 9 for soil media, and 7 for rock media.

[0074] The initial spatial distribution of conductivity is estimated based on the amplitude attenuation characteristics of reflected waves. Amplitude values ​​at different propagation distances are extracted from the same reflecting interface. The attenuation rate of amplitude with propagation distance reflects the conductivity of the medium. The amplitude attenuation coefficient is obtained by fitting a linear relationship between the logarithmic amplitude and the propagation distance. Conductivity and the attenuation coefficient are positively correlated; the larger the attenuation coefficient, the higher the conductivity. When the amplitude attenuation characteristics are not obvious or the signal-to-noise ratio is low, the initial conductivity value is uniformly set to 0.01 S / m as the default configuration. The initialized spatial distribution of relative permittivity and conductivity is stored in a two-dimensional array. The row dimension of the array corresponds to the spatial depth direction, and the column dimension corresponds to the horizontal direction of the survey line. The range of array element values ​​is physically constrained: the relative permittivity is not less than 1 and not greater than 30, and the conductivity is not less than 0 and not greater than 0.5 S / m.

[0075] The computational domain is divided into multiple grid cells. The horizontal range of the computational domain covers the spatial span of the complete set of measurement points, while the vertical range extends from the ground surface to the maximum detection depth. The maximum detection depth is estimated based on the maximum recording time of the radar signal and the average wave velocity of the medium, and is calculated as the recording time multiplied by the average wave velocity and then divided by 2. A regular rectangular grid is used for grid division. The horizontal grid step size is consistent with the measurement point interval or is an integer fraction of the measurement point interval. The vertical grid step size is determined based on the center frequency of the transmitted pulse and the relative permittivity of the medium, set according to the method for determining the spatial step size in the aforementioned forward modeling calculation to ensure numerical calculation accuracy. The total number of grid cells is equal to the number of horizontal cells multiplied by the number of vertical cells. When the computational domain is large enough that the total number of grid cells exceeds one million, a variable mesh technique is used, employing a finer mesh in the target area and a coarser mesh in the boundary area. The step size in the transition area of ​​the mesh size varies geometrically, with a scaling factor set to 1.1 to 1.3.

[0076] The spatial distribution is discretized and assigned to each grid cell. All grid cells are traversed, and the corresponding parameter values ​​are obtained by interpolation from the initialized spatial distribution arrays of relative permittivity and conductivity based on the spatial coordinates of the grid cell center point. Bilinear interpolation is used; when the grid cell center point is located between known data points, the interpolation result is calculated based on the parameter values ​​of four adjacent known data points and distance weights. After assignment, each grid cell stores its corresponding relative permittivity and conductivity values. The data structure uses a structure array or parallel array to support fast indexing and access.

[0077] The propagation equations of electromagnetic waves in the grid cells are established, based on the time-domain form of Maxwell's equations, and solved discretly using the finite-difference time-domain method. Electric and magnetic field components are arranged alternately in a Yee grid, with the electric field component located at the midpoint of the grid edge and the magnetic field component at the center of the grid surface. Electromagnetic boundary conditions at the grid cell interfaces are set based on the relative permittivity and conductivity values ​​of adjacent grid cells. The permittivity at the interface is taken as the arithmetic mean or harmonic mean of the permittivity of the two adjacent cells. The arithmetic mean is applicable when the field component is perpendicular to the interface, while the harmonic mean is applicable when the field component is parallel to the interface. The conductivity is treated in the same way as the permittivity. A perfectly matched layer absorbing boundary condition is used at the outer boundary of the computational domain. The thickness of the perfectly matched layer is set to 10 to 20 grid cells, and the conductivity within the layer increases according to a geometric series or polynomial law, with the maximum conductivity value set as the theoretical optimum.

[0078] The time-domain waveform of the excitation source is defined, and the excitation source location is set above the ground at the same height as the actual radar transmitting antenna. The excitation source type is a current source or an electric field source. The time-domain waveform uses a Ricker wavelet or a Gaussian derivative pulse. The center frequency of the Ricker wavelet is consistent with that of the actual radar system, and the time delay is usually set to 2 to 3 times the pulse duration to ensure that the main lobe of the excitation signal is completely contained within the calculation time window. The time-domain characteristics of the Gaussian derivative pulse are similar to those of the Ricker wavelet, making it suitable for broadband detection scenarios. The discrete sampling point interval of the time-domain waveform is equal to the time step, which is determined according to the Coulomb stability condition and set according to the method for determining the time step in the aforementioned forward modeling calculation.

[0079] The excitation source's time-domain waveform is substituted into the propagation equation, and the solution is obtained by progressively solving the problem according to the time step. The calculation for each time step is divided into two half-steps: the first half-step updates the magnetic field components, and the second half-step updates the electric field components. The magnetic field component update is calculated based on the spatial curl of the electric field component from the previous time step, while the electric field component update is calculated based on the spatial curl of the magnetic field component from the current time step, as well as the relative permittivity and conductivity of the grid cells. The update method uses the same difference scheme as the aforementioned forward modeling calculation. During the calculation, the electric field component values ​​of each grid cell at each time step need to be stored. The storage strategy is to choose full storage or thinned storage based on the memory capacity. Thinned storage is performed every 5 to 10 time steps.

[0080] The electric field component time series of the corresponding grid cells for each measuring point in the complete measuring point set is extracted. The nearest neighbor grid cell is determined based on the spatial coordinates of the measuring point, and the electric field component values ​​of that grid cell are extracted from the stored electric field component data at all time steps. Since the electric field components are vectors, the vertical component or the horizontal component along the measuring line direction is typically extracted as the radar received signal. The length of the extracted time series is equal to the total number of time steps, and the sampling interval is equal to the time step length.

[0081] The extracted electric field component time series are arranged according to the measurement point order to form simulated radar data. The measurement point order is consistent with the spatial numbering of the measurement points in the complete measurement point set. The arrangement is a two-dimensional matrix of N rows and L columns, where N is the total number of complete measurement points and L is the number of time sampling points for a single radar signal. The simulated radar data is stored as a double-precision floating-point array, consistent with the data format of the reconstructed complete radar data, which facilitates subsequent residual calculation and gradient solving.

[0082] This invention generates simulated data through forward modeling based on a physical model, providing an iteratively updated data benchmark for subsequent full-waveform inversion and ensuring the physical consistency of the inversion process.

[0083] Optionally, a residual objective function is constructed between the complete radar data and the simulated radar data. The gradient of the residual objective function with respect to the relative permittivity and conductivity is calculated and its spatial distribution is updated. The optimization is iterated until the residual objective function converges, obtaining the final spatial distribution of the relative permittivity and conductivity of the subsurface medium, including: Calculate the difference between the complete radar data and the simulated radar data at each measurement point location, and sum the squares to construct the residual objective function; Establish the adjoint wave field equation, substitute the difference between the locations of each measuring point as the adjoint source term into the adjoint wave field equation, and obtain the distribution of the adjoint wave field in each grid cell by solving the adjoint wave field equation. The electric field components of each grid cell obtained during the electromagnetic wave forward modeling process are cross-correlated with the distribution of the accompanying wave field in the corresponding grid cell in the time domain to obtain the gradient of the residual objective function with respect to the relative permittivity and conductivity of each grid cell. The update direction of the relative permittivity and conductivity is determined according to the gradient, and the relative permittivity and conductivity values ​​of each grid cell are adjusted according to the update direction. Substitute the updated relative permittivity and conductivity spatial distribution into the electromagnetic wave propagation equation and perform forward modeling again to obtain new simulated radar data. Based on the new simulated radar data, recalculate the residual objective function. When the residual objective function satisfies the convergence condition, the current spatial distribution of relative permittivity and conductivity is taken as the final spatial distribution; if it does not satisfy the condition, the steps of solving the adjoint wave field, calculating the gradient, updating the parameters, and performing forward modeling are repeated until the convergence condition is satisfied.

[0084] In this embodiment, when calculating the difference between complete radar data and simulated radar data at each measurement point location and summing the squares to construct the residual objective function, all measurement point locations in the complete measurement point set are traversed. For each measurement point location, the time series corresponding to the complete radar data and simulated radar data at that measurement point are extracted. The complete radar data is an N-row L-column matrix obtained through compressed sensing reconstruction, and the simulated radar data is an N-row L-column matrix obtained through forward modeling, where N is the total number of complete measurement points and L is the number of time sampling points for a single-channel radar signal. For the i-th measurement point, the difference between the complete radar data value and the simulated radar data value at all time sampling points j is calculated. The difference is equal to the element in the i-th row and j-th column of the complete radar data minus the element in the i-th row and j-th column of the simulated radar data. This difference is squared, and by traversing all measurement point locations i from 1 to N and all time sampling points j from 1 to L, the squares of all differences are summed to obtain the residual objective function value. The smaller the residual objective function value, the higher the degree of fit between the simulated radar data and the complete radar data, and the closer the corresponding subsurface medium parameter model is to the real situation. The calculation result of the residual objective function is stored as a single scalar value, with the data type being a double-precision floating-point number. This value serves as the core indicator for judging the convergence of the iteration.

[0085] When establishing the adjoint wave field equation, it is derived from the forward Maxwell equations by variational derivation of the residual objective function. The electric and magnetic field components of the adjoint equation satisfy a propagation relationship similar to, but in reverse time, as in the forward equations. The temporal evolution of the adjoint electric field is jointly controlled by the spatial curl of the adjoint magnetic field, the relative permittivity of the medium, and the conductivity, while the temporal evolution of the adjoint magnetic field is controlled by the spatial curl of the adjoint electric field. The difference between the complete radar data and the simulated radar data at each measurement point is substituted into the adjoint wave field equation as the adjoint source term. The adjoint source term is spatially located at the grid cell corresponding to each measurement point and temporally corresponds to each time sampling point of the radar signal at that measurement point. The adjoint source term is introduced by adding a source term contribution to the right side of the adjoint electric field equation. The magnitude of the source term contribution is equal to the data difference at that measurement point at that time. The spatial distribution maps the measurement point location to the nearest neighbor grid cell through the Dirac function or interpolation. The adjoint wave field equations are solved using a time-reverse approach. Starting from the maximum recorded moment of the radar signal, the time step gradually decreases towards the zero point, maintaining the same size as in the forward modeling but in the opposite direction. During the time-reverse solution, the values ​​of the adjoint electric and magnetic fields in each grid cell are updated using the same time-domain finite-difference scheme as in the forward modeling, but in the order of updating the adjoint electric field first, followed by the adjoint magnetic field. The solution process requires storing the adjoint electric field components of each grid cell at each time step. The storage strategy is the same as in the forward modeling, allowing for either full storage or thinned storage, with the thinned storage interval set to 5 to 10 time steps. The adjoint wave field also employs a fully matched layer absorbing boundary condition at the boundary of the computational domain to ensure no spurious reflections occur at the boundary.

[0086] After obtaining the distribution of the adjoint wave field in each grid cell by solving the adjoint wave field equation, the electric field components of each grid cell obtained during the electromagnetic wave forward modeling are cross-correlated with the distribution of the adjoint wave field in the corresponding grid cells in the time domain to obtain the gradient of the residual objective function with respect to the relative permittivity and conductivity of each grid cell. The time domain cross-correlation operation is performed independently for each grid cell. For a grid cell at a spatial location r, the electric field component values ​​of that grid cell stored in the forward modeling calculation at all time steps are extracted, and the adjoint electric field component values ​​of that grid cell stored in the adjoint wave field calculation at all time steps are also extracted. For the gradient calculation of the relative permittivity, the time derivative of the forward modeling electric field is multiplied by the adjoint electric field and the dot product operation is performed step by step. The dot product result is the sum of the corresponding components of the two vectors. The dot product results of all time steps are integrated over time. The integration is achieved by discrete summation, that is, the dot product results of each time step are multiplied by the time step size and then accumulated. The integrated result is multiplied by the negative vacuum permittivity to obtain the gradient of the relative permittivity of that grid cell. For calculating the conductivity gradient, the forward electric field and the accompanying electric field are directly multiplied step-by-step. The dot product results for all time steps are then integrated over time, and the negative value of the integral result is used to obtain the conductivity gradient of that grid cell. After the gradient calculation is completed, each grid cell corresponds to a relative permittivity gradient value and a conductivity gradient value. The gradient values ​​are stored as a two-dimensional array corresponding to the grid cell array, and the array elements are double-precision floating-point numbers.

[0087] When determining the update direction for relative permittivity and conductivity based on the calculated gradient, the steepest descent method, the conjugate gradient method, or the quasi-Newton method can be used. With the steepest descent method, the update direction is directly taken as the negative direction of the gradient, i.e., the update direction equals the gradient multiplied by -1. With the conjugate gradient method, the update direction is a weighted combination of the negative direction of the current gradient and the update direction of the previous iteration. The weighting coefficients are calculated according to the Fletcher-Reeves formula or the Polak-Ribiere formula, and are equal to the square of the 2-norm of the current gradient divided by the square of the 2-norm of the previous gradient. With the L-BFGS quasi-Newton method, the update direction is obtained by multiplying the approximate Hessian inverse matrix by the negative direction of the gradient. The approximation of the Hessian inverse matrix uses limited memory to store the gradient differences and parameter differences of the most recent 5 to 10 iterations, avoiding the memory overhead caused by storing the complete Hessian matrix. After determining the update direction, the update step size is determined using a line search method. The goal of the line search is to find the step size that maximizes the decrease in the residual objective function. The line search employs either backtracking line search or exact line search. Backtracking line search starts with a large initial step size and gradually reduces it until the Armijo or Wolfe condition is satisfied. The initial step size is set to 0.1, the reduction factor is set to 0.5, and the maximum number of backtracking iterations is set to 20. Exact line search calculates the objective function values ​​corresponding to multiple candidate step sizes in the update direction and selects the step size that minimizes the objective function. Candidate step sizes are sampled logarithmically within the range of 0.001 to 0.5.

[0088] The relative permittivity and conductivity values ​​of each grid cell are adjusted according to the update direction and determined step size. For each grid cell, the new relative permittivity value is equal to the old value minus the step size multiplied by the update direction component of the relative permittivity of that grid cell, and the new conductivity value is equal to the old value minus the step size multiplied by the update direction component of the conductivity of that grid cell. The updated parameter values ​​must be constrained within a physically reasonable range. The constraint range for the relative permittivity is set according to the geological type of the area to be measured: 1.5 to 5 for snow and ice media, 3 to 25 for soil media, and 4 to 15 for rock media. If the updated value is less than the lower limit, it is forcibly assigned to the lower limit; if it is greater than the upper limit, it is forcibly assigned to the upper limit. The constraint range for conductivity is 0 to 0.2 S / m; negative values ​​are forcibly assigned to 0, and values ​​exceeding 0.2 S / m are forcibly assigned to 0.2 S / m. The parameter update process adopts a grid cell-by-grid method. After the update, the relative permittivity and conductivity of all grid cells constitute a new spatial distribution.

[0089] The updated relative permittivity and conductivity spatial distribution are substituted into the electromagnetic wave propagation equations to perform a new forward modeling calculation to obtain new simulated radar data. The forward modeling method is exactly the same as the previous forward modeling steps, using the finite-difference time-domain method to solve Maxwell's equations, while keeping the mesh generation, time step, boundary conditions, and excitation source settings unchanged. After the forward modeling is completed, the electric field component time series of the corresponding mesh elements at each measurement point in the complete measurement point set are extracted and arranged in the order of the measurement points to form a new simulated radar data matrix. Based on the new simulated radar data, the residual objective function is recalculated using the same method as the initial calculation: traversing all measurement points and time sampling points, the sum of squared differences between the complete radar data and the new simulated radar data is calculated to obtain the new residual objective function value.

[0090] When determining whether the residual objective function meets the convergence criteria, the convergence criteria include the relative change criterion and the absolute change criterion. The relative change criterion is calculated by dividing the absolute difference between the residual objective function values ​​of two consecutive iterations by the previous residual objective function value. Convergence is determined when this relative change is less than a set threshold, which is usually set to 1 × 10⁻⁶. -4 Up to 1×10 -3 The absolute change criterion directly calculates the absolute difference between the residual objective function values ​​of two adjacent iterations. Convergence is determined when this absolute difference is less than a set threshold, which is set based on the data amplitude and is typically 1 × 10⁻¹⁰ of the square of the complete radar data amplitude. -6 Up to 1×10 -5 In addition to the convergence criterion for the objective function, a maximum number of iterations is set as the termination condition, ranging from 30 to 50 iterations. The iteration terminates regardless of convergence when the maximum number of iterations is reached. When the residual objective function satisfies the convergence condition, the current spatial distribution of relative permittivity and conductivity is output as the final spatial distribution. The output data format is a two-dimensional array, with rows corresponding to the depth direction and columns corresponding to the horizontal direction. The array elements are the parameter values ​​of each grid cell. If the convergence condition is not met, the steps of solving the adjoint wave field, calculating the gradient, updating parameters, and performing forward modeling are repeated to enter the next iteration cycle. The iteration counter is incremented by 1, and the iteration history information is updated for direction calculation using the conjugate gradient method or the quasi-Newton method.

[0091] This invention achieves high-precision quantitative inversion of electromagnetic parameters of underground media by efficiently calculating gradients using the adjoint state method and iteratively updating parameters using an optimization algorithm, thereby improving detection accuracy.

Claims

1. A method for acquiring and inverting airborne ground-penetrating radar data based on compressed sensing technology, characterized in that, include: Based on the UAV aerial ground-penetrating radar survey lines planned for the area to be measured, a complete set of measurement points is determined. Based on the sub-sampling rate, measurement points are selected from the complete set of measurement points to form a sub-sampling measurement point set, and undersampled radar data is collected. The theoretical radar response signal is obtained by forward modeling the electromagnetic wave propagation equation using underground medium parameters and then normalized to form a theoretical prior sparse dictionary. A sparse dictionary of measured features is constructed by performing feature decomposition on measured radar echo signal samples. The theoretical prior sparse dictionary is combined with the measured feature sparse dictionary to form a joint sparse dictionary; Establish a sparse representation relationship between the undersampled radar data, the joint sparse dictionary, and the sparse coefficient vector. Reconstruct the complete radar data by solving the sparse coefficient vector and multiplying it with the joint sparse dictionary. The spatial distribution of the relative permittivity and conductivity of the subsurface medium is defined, and simulated radar data is obtained by forward modeling using the electromagnetic wave propagation equation. A residual objective function is constructed between the complete radar data and the simulated radar data. The gradient of the residual objective function with respect to the relative permittivity and conductivity is calculated and its spatial distribution is updated. The optimization is iterated until the residual objective function converges, and the final spatial distribution of the relative permittivity and conductivity of the subsurface medium is obtained.

2. The method according to claim 1, characterized in that, Based on the sub-sampling rate, measurement points are selected from the complete measurement point set to form a sub-sampling measurement point set, including: Based on the prior information of the topographic relief and underground structure of the area to be measured, a sub-sampling rate for spatial variation is set, and measuring points are selected from the complete set of measuring points according to the sub-sampling rate to form a sub-sampling measuring point set.

3. The method according to claim 1, characterized in that, The theoretical radar response signal is obtained through forward modeling of the electromagnetic wave propagation equation using underground medium parameters, and then normalized to construct a theoretical prior sparse dictionary, including: Multiple combinations of relative permittivity and conductivity parameters of underground media are set; for each combination of parameters, the corresponding theoretical radar response signal is obtained by forward modeling of the electromagnetic wave propagation equation; the amplitude of each theoretical radar response signal is normalized; the multiple normalized theoretical radar response signals are arranged in columns to form a theoretical prior sparse dictionary.

4. The method according to claim 1, characterized in that, The measured radar echo signal samples are subjected to feature decomposition to construct a sparse dictionary of measured features, including: Multiple sets of measured radar echo signal samples are selected and arranged in columns to form a measured data matrix; the measured data matrix is ​​subjected to feature decomposition to extract basis vectors; the extracted basis vectors are arranged in columns to form a measured feature sparse dictionary.

5. The method according to claim 1, characterized in that, Complete radar data is reconstructed by solving the sparse coefficient vector and multiplying it by the joint sparse dictionary, including: The complete radar data to be reconstructed is represented as the product of the joint sparse dictionary and the sparse coefficient vector; a reconstruction objective function is constructed, which includes an error term between the undersampled radar data and the reconstructed data and a sparse constraint term of the sparse coefficient vector; the sparse coefficient vector is solved by a sparse optimization algorithm to minimize the reconstruction objective function; the complete radar data is obtained by multiplying the solved sparse coefficient vector with the joint sparse dictionary.

6. The method according to claim 1, characterized in that, The spatial distribution of the relative permittivity and conductivity of the underground medium is defined, and simulated radar data is obtained through forward modeling of the electromagnetic wave propagation equation, including: The spatial distribution of the relative permittivity and conductivity of the underground medium is initialized based on the arrival time and amplitude characteristics of the reflected waves in the complete radar data. The computational region is divided into multiple grid cells, and the spatial distribution is discretized and assigned to each grid cell. The propagation equation of electromagnetic waves in the grid cells is established. The electromagnetic field boundary conditions at the grid cell interface are set according to the relative permittivity and conductivity values ​​of adjacent grid cells. The time-domain waveform of the excitation source is set and substituted into the propagation equation. The electric field components of each grid cell are obtained by discretizing and solving the propagation equation and boundary conditions. The time series of electric field components corresponding to the grid cells at each measurement point in the complete set of measurement points are extracted. The extracted electric field component time series are arranged in the order of measurement points to form simulated radar data.

7. The method according to claim 1, characterized in that, Construct a residual objective function between the complete radar data and the simulated radar data, calculate the gradient of the residual objective function with respect to the relative permittivity and conductivity, update its spatial distribution, and iterate until the residual objective function converges to obtain the final spatial distribution of the relative permittivity and conductivity of the subsurface medium, including: Calculate the difference between the complete radar data and the simulated radar data at each measurement point location, and sum the squares to construct the residual objective function; Establish the adjoint wave field equation, substitute the difference between the locations of each measuring point as the adjoint source term into the adjoint wave field equation, and obtain the distribution of the adjoint wave field in each grid cell by solving the adjoint wave field equation. The electric field components of each grid cell obtained during the electromagnetic wave forward modeling process are cross-correlated with the distribution of the accompanying wave field in the corresponding grid cell in the time domain to obtain the gradient of the residual objective function with respect to the relative permittivity and conductivity of each grid cell. The update direction of the relative permittivity and conductivity is determined according to the gradient, and the relative permittivity and conductivity values ​​of each grid cell are adjusted according to the update direction. Substitute the updated relative permittivity and conductivity spatial distribution into the electromagnetic wave propagation equation and perform forward modeling again to obtain new simulated radar data. Based on the new simulated radar data, recalculate the residual objective function. When the residual objective function satisfies the convergence condition, the current spatial distribution of relative permittivity and conductivity is taken as the final spatial distribution.

Images