Low complexity rough road echo time delay estimation method based on sparse frequency sampling

Abstract

The application discloses a low-complexity rough road echo time delay estimation method based on sparse frequency sampling, and is used for solving the problems of high calculation complexity and poor estimation accuracy of existing ground penetrating radar echo time delay estimation technology. In the application, a grouping sparse frequency sampling method is adopted, the data dimension is reduced through sparse sampling, and the calculation complexity is greatly reduced. Meanwhile, cross-correlation detection is adopted to solve the time domain ambiguity caused by sparse sampling. In addition, a time delay iteration screening method is adopted, the preliminary time delay estimation value result is further iterated and screened, false time delay estimation values can be effectively avoided from being mistaken for real values, and higher estimation accuracy can be obtained.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of ground penetrating radar, time delay estimation and coherent echo recognition, and particularly relates to a low-complexity rough pavement echo time delay estimation method based on sparse frequency sampling. BACKGROUND

[0002] As a widely used non-destructive testing method in civil engineering, ground penetrating radar (GPR) is applied to pavement detection, and through detecting the reflected echo signals from underground structures, fast data acquisition is realized. In pavement detection, assuming that the pavement structure is horizontally layered, through time delay estimation and amplitude estimation of the echo signals, the vertical structure of the pavement can be inferred from the radar profile. Due to the existence of thin pavement structure, the reflected echo signals of each layer are mutually overlapped (the time delay bandwidth product of the signals is less than 1), and there is strong correlation or even coherence; and the interface roughness of each layer directly affects the frequency characteristics of the echo signals, therefore, the fast processing of the echo signals with overlapping, coherence and changing frequency characteristics has become the main challenge of GPR time delay estimation.

[0003] In order to solve the echo overlapping problem under limited bandwidth, high-resolution algorithms are widely used in the field of time delay estimation, among which the multiple signal classification (MUSIC) of the subspace class can provide asymptotically unbiased high-resolution time delay estimation under the condition that the echoes are incoherent. However, in actual engineering environment, the reflected echoes of each layer are highly correlated or even coherent, which will cause the rank deficiency of the covariance matrix of the received signals, thus leading to the failure of the MUSIC algorithm. For the echoes with linear frequency characteristics, the spatial smoothing preprocessing (SSP) technology is used to restore the rank of the data covariance matrix, so as to realize echo decorrelation and ensure the performance of time delay estimation.

[0004] However, recent studies have shown that for ultra-wideband GPR (bandwidth > 2GHz), the frequency characteristics of the backscattered echoes of rough pavement are approximately nonlinear Gaussian functions, in which case the decorrelation technology based on SSP will no longer be applicable. In order to eliminate the influence of the nonlinear frequency characteristics of the echoes, M. Sun et al. proposed to use the interpolation spatial smoothing technology combined with modified MUSIC (modified-MUSIC) for time delay estimation, and considered several possible frequency characteristics. However, this algorithm needs to perform multiple eigenvalue decomposition (EVD), and assuming that the data length is L, the computational complexity of a single EVD is O(L 3 ), therefore, the overall computational complexity of the system is extremely high, and the algorithm has poor robustness. Therefore, how to realize high-precision GPR time delay estimation with low computational complexity and considering nonlinear frequency characteristics has recently attracted great attention.

[0005] For the classical continuous frequency sampling model, in order to avoid ambiguity in the time domain search range, according to the Nyquist sampling theorem, for a fixed bandwidth B, it is necessary to ensure that the frequency interval of each sampling point Δf=B / N is not less than 1 / (2τ max ), wherein τ max represents the maximum estimable time delay, and N represents the number of frequency sampling points, therefore, the value of N cannot be too small, and this will bring higher computational load. In recent years, sparse sampling structure has been proposed and widely applied in array processing field, in order to eliminate the artifact problem caused by sparse sampling, researchers have successively proposed estimation methods based on virtual extension of differential array, Toeplitz matrix reconstruction and array interpolation, wherein the virtual extension method based on differential array cannot process coherent echoes; the Toeplitz reconstruction method requires the covariance matrix to have a joint diagonal structure; the array interpolation method will bring additional computational load, and the estimation accuracy is not high. In view of this, Huimin Pan et al. proposed a specific symmetric sparse sampling structure, combined with improved MUSIC (IMUSIC) to realize the time delay estimation of coherent echoes, but this method can only be applied to specific symmetric sampling structure, and due to the non-uniform sparse in frequency domain, it is difficult to combine with the interpolation method to process the nonlinear frequency characteristics of echoes, and it is not suitable for time delay estimation of rough road echoes.

[0006] In summary, the existing problems of the prior art are: 1) considering the influence of interface roughness, but only applicable to non-sparse uniform sampling, high complexity; 2) sparse sampling reduces a certain complexity, but the array is limited and difficult to combine with interpolation method, and cannot process rough road echoes. SUMMARY

[0007] The purpose of the present application is to solve the problems of high computational complexity and poor estimation accuracy of the existing time delay estimation method considering interface roughness, and to provide a low-complexity rough road echo time delay estimation method based on sparse frequency sampling, so as to improve the performance of ground penetrating radar time delay estimation.

[0008] The purpose of the present application can be achieved by adopting the following technical solutions:

[0009] A low-complexity rough road echo time delay estimation method based on sparse frequency sampling, the estimation method comprising the following steps:

[0010] The transmitting end transmits a ground penetrating radar signal;

[0011] The ground penetrating radar signal reaches the multi-layer rough road and is reflected by each layer interface to form multiple backscattering echoes;

[0012] The receiving end detects the backscattering echoes to obtain initial observation data;

[0013] The initial observation data is whitened, and the received data is obtained after eliminating the radar pulse from the initial observation data;

[0014] The received data is grouped and sparsely frequency sampled to obtain a plurality of low-dimensional received data sub-sample sets;

[0015] The covariance matrix of each sub-sample set is calculated;

[0016] The linear interpolation is performed on the covariance matrix of each sub-sample set;

[0017] Each linearly interpolated covariance matrix is subtracted by a respective noise matrix;

[0018] The spatial smoothing preprocessing is performed on each covariance matrix after subtracting the noise matrix;

[0019] The modified-MUSIC is applied to each spatially smoothed covariance matrix to obtain a plurality of time delay estimation spectra;

[0020] The minimum time delay estimation spectrum period T is calculated;

[0021] The spectrum peak search is performed on each time delay estimation spectrum to obtain a plurality of time delay estimation values;

[0022] The time delay of the first echo is calculated: the cross-correlation operation is performed on the plurality of time delay estimation spectra to obtain a cross-correlation spectrum, and the time delay τ1 of the first echo is obtained according to the time corresponding to the maximum value of the cross-correlation spectrum;

[0023] Each time delay estimation value obtained by the spectrum peak search is selected, and the time delay value in each time delay estimation value belonging to the interval [τ1, τ1+T) is selected and recorded in an array;

[0024] The row sum of the received data is calculated and averaged to obtain a column vector X;

[0025] The time delay values in the array are constructed into a complete dictionary D, and K correct time delay values are screened out through K iterations for a given number of echoes K, wherein the iteration screening process is as shown in Figure 1 and specifically includes the following steps:

[0026] T1, initialize the iteration number i=1, the initial residual b=X, the initial selected matrix S=[], the initial amplitude matrix a=[], and the initial time delay screening result t3=[];

[0027] T2, the contribution of each column vector in the complete dictionary D to the residual b is calculated by least square, the column vector with the maximum contribution value is selected as u, and the time delay value corresponding to u is recorded as τ;

[0028] T3, update the amplitude matrix a = u\X, select the matrix S = [Su];

[0029] T4, update the residual b = X-aS, update t3 = [t3τ], update i = i+1, delete the element u from the complete dictionary D;

[0030] Repeat the above T2-T4 operations until i>K, end the screening, and output the time delay screening result t3 as an array containing K elements.

[0031] Further, in the grouped sparse frequency sampling, m integers that are prime to each other are selected, and the i-th integer is represented as M i i = 1…m, respectively, obtain the i-th set of sub-sampling sets r i of the i-th group of received data with a frequency sampling interval of M i i = 1…m, wherein Δf represents the minimum frequency interval that satisfies the Nyquist sampling theorem, and taking m = 2 as an example, the grouped sparse frequency sampling structure is as shown in Figure 2 The grouped sparse frequency sampling method is simple to implement, and each sub-sampling set is uniformly sparse sampling, which can be directly combined with subsequent linear interpolation to solve the nonlinear frequency characteristics of the echo.

[0032] Further, in the linear interpolation of the covariance matrix of each set of sub-sampling sets, the i-th set of covariance matrices R i i = 1…m constructed from the i-th set of received data sub-sampling sets r are linearly interpolated, wherein the superscript H represents the complex conjugate transpose operation, and m linearly interpolated covariance matrices are obtained, and the i-th linearly interpolated covariance matrix is specifically represented as:

[0033]

[0034] wherein B represents an interpolation transformation matrix;

[0035] The covariance matrix after linear interpolation eliminates the influence of the nonlinear frequency characteristics of the echo, and satisfies the Vandermonde structure on the mathematical model, and can be combined with subsequent spatial smoothing preprocessing to solve the coherence problem of the echo.

[0036] Further, in the noise matrix subtraction of each sub-sampling set after linear interpolation, since the covariance matrix linear interpolation step introduces a colored noise term, in order to ensure the subsequent step expansion, the noise component of the covariance matrix needs to be removed, and the noise variance of the i-th linearly interpolated covariance matrix R is estimated to obtain the noise variance of the i-th linearly interpolated covariance matrix The i-th noise-free covariance matrix is calculated as follows:

[0037]

[0038] where Σ = Λ -1 Λ -H represents the diagonal matrix of the frequency domain radar pulse.

[0039] Further, the noise variance estimation is calculated by a propagation operator PM method, which estimates the noise variance without eigenvalue decomposition of the covariance matrix, has a smaller calculation complexity, and has the following process:

[0040] A given data covariance matrix R xx is divided into the following form:

[0041]

[0042] where G1, G2, H1 and H2 are (KxK), ((N-K)xK), (Kx(N-K)) and ((N-K)x(N-K)) dimensional sub-matrices, N and K represent the dimensions of the data covariance matrix R xx and the number of echoes, respectively, and the noise variance σ 2 of the data covariance matrix R xx is derived as follows:

[0043]

[0044] where, tr{·} represents the trace operator of the matrix, the superscript + represents the pseudo-inverse of the matrix, I N-k represents the (N-k)x(N-k) dimensional unit matrix.

[0045] Further, in the spatial smoothing preprocessing, in order to eliminate the influence of the coherence of the echoes, the following modification is made to the ith noiseless covariance matrix to obtain the ith spatially smoothed preprocessed covariance matrix as follows:

[0046]

[0047] where, represents the lth sub-band of the ith noiseless covariance matrix L i is an integer greater than 1 and less than half the length;

[0048] After spatial smoothing preprocessing, the rank of the covariance matrix can be restored, thereby realizing de-coherence of the echoes; the dimension υ i of the ith covariance matrix after spatial smoothing preprocessing is N i -Li +1, where N i is the length of the i-th subsample set r i while the dimension of the i-th spatially smoothed preprocessed covariance matrix υ i is not less than 2K.

[0049] Further, the modified MUSIC is used to calculate the delay estimation spectrum because the echo frequency characteristics are unknown, theoretically, the original MUSIC cannot be directly used, a noise component is introduced in the expression of the original modified MUSIC delay estimation spectrum, the expression of the original modified MUSIC delay estimation spectrum is modified, and the modified MUSIC is used for delay estimation, the process is as follows:

[0050] The i-th spatially smoothed preprocessed covariance matrix R is subjected to eigenvalue decomposition to obtain the corresponding i-th noise subspace U ni , i = 1 … m.

[0051] Subsequently, the time domain search range t = t1: Δt: t2 of the modified MUSIC is set, where t1 and t2 are the minimum and maximum delays respectively, and Δt is the time step;

[0052] The orthogonality of the signal subspace and the noise subspace is utilized to construct the modified MUSIC delay estimation spectrum, the expression of the modified MUSIC delay estimation spectrum does not need to distinguish between the i-th spatially smoothed preprocessed covariance matrix υ i is odd and even, and the introduction of the noise component is avoided, so that more accurate delay estimation can be obtained;

[0053] The i-th modified MUSIC delay estimation spectrum expression is as follows:

[0054]

[0055] Where λ ik (t), i = 1 … m; k = 1, 2 … 2K represents the k-th smaller eigenvalue of R , real{·} represents the real part operation, K is the number of echoes, is the diagonal matrix of the steering vector of the i-th spatially smoothed preprocessed covariance matrix R .

[0056] Substituting calculation obtains the i-th modified delay estimation spectrum P i (t), i = 1 … m.

[0057] Further, the minimum time delay estimation spectrum period calculation method is as follows:

[0058] For the m integers defined in the above steps, the i integer is represented as M i , i = 1…m, and M i is the maximum integer in M max ;

[0059] Calculate the minimum time delay estimation spectrum period:

[0060] T = (t2-t1) / (2(M max )-1)

[0061] Wherein, t1 and t2 are the minimum and maximum time delays in the modified-MUSIC time domain search range respectively.

[0062] Further, the time delay of the first echo utilizes the characteristics of the rough road backscattering echo, and according to the reflection coefficient law of each layer interface, the amplitude of each echo is attenuated in turn, so that the time corresponding to the maximum value of the cross-correlation spectrum is detected as the time delay of the first echo.

[0063] Further, the ground penetrating radar signal transmitted by the transmitting end adopts an ultra-wideband stepped frequency radar, and the ultra-wideband stepped frequency radar has rich high frequency information in the frequency spectrum, and can obtain high time and space resolution.

[0064] The present application has the following advantages and effects relative to the prior art:

[0065] (1) The present application groups the received data for sparse frequency sampling, obtains multiple low-dimensional sub-sampling sets, greatly reduces the computational complexity of subsequent time delay estimation; each sub-sampling set is a sparse uniform sampling, simple to implement, free in structure, and can directly combine linear interpolation to solve the nonlinear frequency characteristics of the echo, and the sub-sampling set covariance matrix after linear interpolation satisfies the Vandermonde structure, which can directly combine spatial smoothing preprocessing to solve the coherence problem of the echo;

[0066] (2) The present application modifies the time delay estimation spectrum expression of modified-MUSIC, applies modified-MUSIC to the covariance matrix of each sub-sampling set, and calculates their respective time delay estimation spectrum;

[0067] (3) The present application proposes to perform cross-correlation detection on multiple time delay estimation spectra, obtain the time delay τ of the first echo from the cross-correlation spectrum, and then calculate the minimum time delay estimation spectrum period T to obtain a correct time delay estimation interval [τ, τ+T); the method utilizes the attenuation characteristics of the rough road echo signal and the periodic repetition characteristics of the time delay estimation spectrum caused by sparse frequency sampling, and through cross-correlation operation on multiple time delay estimation spectra, the correct time delay estimation interval and the time delay estimation value within the interval are obtained, successfully avoiding the time domain ambiguity problem caused by sparse frequency sampling.

[0068] (4) The present application proposes an iterative screening method, which performs iterative screening on the time delay estimation value in the time delay estimation interval [τ, τ+T), can distinguish the false time delay estimation value and the true time delay estimation value in the correct interval range, and finally screens out the correct time delay estimation value, solving the problem of poor estimation accuracy of the original modified-MUSIC time delay estimation, and greatly improving the overall estimation performance. BRIEF DESCRIPTION OF DRAWINGS

[0069] The drawings described herein are used to provide further understanding of the present application, and form a part of the present application. The illustrative embodiments of the present application and their descriptions serve to explain the present application, and do not constitute an improper limitation on the present application. In the drawings:

[0070] Figure 1 is a flow chart of the time delay iterative screening method disclosed by the present application;

[0071] Figure 2 is a schematic diagram of the grouping sparse frequency sampling structure disclosed by the present application;

[0072] Figure 3 is a block diagram of the specific implementation in the embodiment of the present application;

[0073] Figure 4 is a schematic diagram of the rough road structure in the embodiment of the present application;

[0074] Figure 5 is a simulation result diagram under different signal-to-noise ratio conditions in the embodiment of the present application. DETAILED DESCRIPTION

[0075] In order to make the purpose, technical scheme and advantages of the embodiments of the present application more clear, the technical scheme in the embodiments of the present application will be described clearly and completely below with reference to the drawings of the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0076] EMBODIMENT

[0077] The embodiment specifically discloses a low-complexity rough road echo time delay estimation method based on sparse frequency sampling, and a specific implementation manner is as shown in the figure Figure 3 , and comprises the following steps:

[0078] S1, defining a received signal model: considering a three-layer rough road model, the received signal is the superposition of three backscattering echoes, and the road structure is given in Figure 4 Table 1:

[0079] Table 1. Road structure simulation parameters

[0080]

[0081] The radar data model is used to represent the backscattering echo of the medium separated by the rough interface, and the frequency domain expression is as follows:

[0082]

[0083] Wherein, f represents the frequency, e(f) represents the radar pulse at the frequency f, K is the number of backscattering echoes, K=3 in the embodiment; t k and s k represent the propagation time delay and reflection coefficient of the kth echo respectively, w k (f) represents the frequency characteristics of the kth backscattering echo at the frequency f, and n(f) is the frequency domain additive white Gaussian noise.

[0084] For N=72 non-sparse discrete frequency sampling points in the entire bandwidth B=3GHz, the above radar data model is rewritten into a vector matrix form as follows:

[0085] r=ΛAs+n

[0086] Wherein, r represents an N×1-dimensional received signal vector obtained by step frequency radar measurement, r=[r(f1)r(f2)…r(f N )] T , wherein the i-th frequency sampling point f i =f1+(i-1)Δf, i=1, 2…N, the radar starting frequency f1=0.5GHz, is the frequency step, and the superscript T represents the transposition operation; Λ is an N×N diagonal matrix, and the diagonal elements are the frequency domain form of the radar pulse; A represents a mode matrix, A=[W1a(t1)W2a(t2)…W K a(tk)], wherein W k a(t k ) represents a mode vector of the kth echo, wherein a(t k ) is a direction vector, which contains the time delay information of the echo, W k is an N x N diagonal matrix whose diagonal elements are the frequency characteristics of the kth echo, W k = diag{w k (f1)w k (f2)…w k (f N ), where diag{·} denotes a diagonal matrix; s represents a K x 1 dimensional vector composed of medium reflection coefficients, s = [s1 s2 s3] T ; n is an N x 1 dimensional noise vector with a mean of zero and a covariance matrix of σ 2 I, I represents an N x N dimensional unit matrix.

[0087] In order to use the feature structure-based algorithm in subsequent use, the radar signal needs to be whitened, and the received signal vector r is divided by the diagonal matrix Λ of the radar pulse, and finally the whitened received signal vector r w is obtained as follows:

[0088] r w = Λ -1 r = As + b

[0089] where b is a new noise vector obtained after the data is whitened.

[0090] S2, obtain the grouping sparse frequency sampling signal covariance matrix: select two prime integers M1 = 5 and M2 = 7, and obtain sparse frequency sub-sampling sets r1 and r2 with frequency sampling intervals of M1Δf and M2Δf, and the frequency sampling point numbers of r1 and r2 are N1 = 11 and N2 = 15, respectively.

[0091] Calculate the covariance matrix of the ith sub-sampling set r i , i = 1, 2 , where the superscript H represents the complex conjugate transpose operation.

[0092] S3, construct an interpolation transformation matrix: the signal bandwidth B = 3 GHz, and the frequency characteristics of the echo are approximately Gaussian functions, i.e. w k (f) = exp(-b k f 2 ), where b k represents the roughness of the kth echo reflection interface, which is related to the root mean square height σ h and the correlation length L h of the interface, and C(σ h , L h ) is the diagonal matrix of the frequency characteristics of the backscattering echo of a certain interface.

[0093] Define a set of G = 5 root mean square heights σ = {σh1 ,σ h2 …σ hG},and a set of P=5 possible correlation lengths L={L h1 ,L h2 …L hP},the model vector associated with the set σ and L is calculated and denoted as a matrix as follows:

[0094] C r =[C(σ h1 ,L h1 )C(σ h1 ,L h2 )…C(σ h1 ,L hP )C(σ h2 ,L h1 )…C(σ hG ,L hP )]

[0095] C r contains the echo frequency characteristics in all the considered ranges, then a virtual linear frequency characteristics matrix C r is defined, which has one-to-one correspondence with the elements in C v . Different from C r , each element in C v has a uniform linear frequency characteristic, and the matrix is as follows:

[0096]

[0097] where, i=1,2,…,G;j=1,2,…,P represents the linear frequency characteristic diagonal matrix of the backscattering echo of the rough interface with root mean square height σ hi and correlation length L hj , which corresponds to C(σ hi ,L hj ).

[0098] The least square solution of BC r =C v is calculated to obtain the value of the optimal interpolation transformation matrix B as follows:

[0099]

[0100] S4, Linear interpolation of covariance matrix: based on the obtained transformation matrix B, the i-th sub-sampling set covariance matrix R i ,i=1,2 is linearly interpolated to obtain two new covariance matrices:

[0101]

[0102] The ith interpolated covariance matrix is estimated by the propagator method (PM) The noise variance of the ith interpolated covariance matrix Then, let Subtract the corresponding noise matrix to obtain the noise-free covariance matrix of the ith subsample set:

[0103]

[0104] Where, ∑=Λ -1 Λ -H , Λ is the diagonal matrix of the frequency domain radar pulse in S1.

[0105] S5, spatial smoothing preprocessing: the ith noise-free covariance matrix is subjected to spatial smoothing preprocessing (SSP) to eliminate the coherence of the echo, and the SSP operation makes the following modification to the ith noise-free covariance matrix to obtain the ith coherence-eliminated covariance matrix As follows:

[0106]

[0107] Where, represents the lth subband of the ith noise-free covariance matrix , and let L1=L2=K=3, L i represents the dimension of the ith data covariance matrix after SSP.

[0108] S6, application of modified MUSIC:

[0109] First, the eigenvalue decomposition is performed on the ith SSP covariance matrix to obtain the ith noise subspace U ni ;

[0110] Then, set the time domain search range of the MUSIC algorithm t=t1:Δt:t2, where the minimum search delay t1=0, the maximum search delay t2=12.5nS, and the time step Δt=10 -3 nS;

[0111] Further, the steering vector of the ith SSP covariance matrix is calculated Where, a i (t) represents the direction vector of the ith subsample set r i .

[0112] The modified MUSIC time delay estimation spectrum expression is modified to obtain the modified time delay estimation spectrum expression as follows:

[0113]

[0114] Where, λ ik (t), i = 1…m; k = 1, 2…2K represents The kth smallest eigenvalue;

[0115] By utilizing the orthogonality of the signal subspace and the noise subspace, a modified-MUSIC time delay estimation spectrum is constructed, and two modified time delay estimation spectra, P1(t) and P2(t), are calculated.

[0116] S7. Time Delay Estimation Spectral Peak Search: Subsequently, after normalizing P1(t) and P2(t) respectively, a spectral peak search is performed to find the time delay positions corresponding to the peaks with amplitudes greater than 0.3, and these positions are recorded as an array t1 = [tde1, tde2] without duplicate elements, where tde... i P represents i (t) An array consisting of the estimated peak search delay values;

[0117] Because sparse frequency sampling is used, t1 contains false delay estimates. In order to find the true value, a further fusion delay estimation method is needed to obtain the correct estimation interval and the delay estimate within the interval.

[0118] S8. Subsample set fusion delay estimation: In order to eliminate false delay estimates in the delay set t1 of S7, a fusion delay estimation method is adopted to find a common solution by combining the delay estimation spectra of the two subsample sets.

[0119] First, regarding P i Perform cross-correlation operations on (t), i = 1, 2 to obtain P crx (t)=P1(t)P2(t);

[0120] Based on the characteristics of backscattered echoes from rough road surfaces, the amplitude of each path decreases sequentially. This characteristic can be used to detect P. crx The time corresponding to the maximum cross-correlation peak in (t) is the time delay estimate τ1 of the first echo;

[0121] Calculate the minimum time delay to estimate the spectral period The time delay estimation interval is obtained as [τ1, τ1+T]. Elements in t that belong to the interval [τ1, τ1+T] are selected and recorded as t2. At this time, t2 is an array containing d time delay estimates (d≥K).

[0122] S9. Iterative selection of candidate delay sets:

[0123] If the number of elements d in t2 is equal to K, the estimation process is completed, and the time delay estimation result is output as t2, and the step ends;

[0124] If the number of elements d in t2 is greater than K, the candidate time delay set screening based on the iterative least square is started;

[0125] First, the observation signal X is constructed. Since the frequency domain whitening received signal r w is a multi-snapshot data, the snapshot number snp = 100, and thus r w is an (N x snp) dimensional matrix. The sum of the rows of r w is taken and averaged to form an N x 1 dimensional column vector, which is the observation signal X.

[0126] Subsequently, the d elements in t2 are taken as time delay candidate values (d > K), and an N x d dimensional complete dictionary D is constructed, wherein each column vector in the complete dictionary D corresponds to each time delay in t2.

[0127] Given the iteration number K, the complete dictionary D and the observation signal X, the iterative screening is started, and the specific process is as follows:

[0128] T1, initialize the iteration number i = 1, the initial residual b = X, the initial selected matrix S = [], the initial amplitude matrix a = [], and the initial time delay screening result t3 = [];

[0129] T2, calculate the contribution of each column vector in the complete dictionary D to the residual b by least square, select the column vector with the largest contribution value as u, and record the time delay value corresponding to u as τ;

[0130] T3, update the amplitude matrix a = u\X, and the selected matrix S = [Su];

[0131] T4, update the residual b = X-aS, update t3 = [t3τ], update i = i+1, and delete the element u from the complete dictionary D;

[0132] T5, repeat the above T2-T4 operations until i > K, end the screening, and output the time delay screening result t3 as an array containing K elements.

[0133] After K iterations of the above steps, K accurate time delay estimations are finally obtained, recorded as a (1 x K) dimensional array t3, the estimation process is completed, and the time delay estimation result is output as t3.

[0134] As a contrast, in the modified-MUSIC algorithm proposed by M. Sun, the spectral peak size of the time delay estimation spectrum P(t) is directly taken as a basis to select the time corresponding to the K largest peaks as the estimation value, which is extremely easy to select false time delay estimation value, thereby serious estimation deviation occurs; and the K estimation values obtained through iteration and screening effectively avoid misjudging false time delay estimation value as true value, thereby obtaining higher estimation accuracy.

[0135] The application effect of the present application will be described in detail in combination with simulation results.

[0136] (1) Complexity analysis:

[0137] In the MUSIC type time delay estimation algorithm, the eigenvalue decomposition and peak search contribute most of the calculation amount, in the present embodiment, the calculation complexity of single eigenvalue decomposition and spectrum peak search of the rough road echo time delay estimation method based on sparse frequency sampling adopted in the present embodiment is Wherein the first term and the second term are the calculation complexity sum brought by the eigenvalue decomposition and spectrum peak search of the two parts of the sub-sampling sets r1 and r2, wherein n is the grid number of the spectrum peak search, m1=N1-L1+1=9, m2=N2-L2+1=13;

[0138] And under the same grid number of the spectrum peak search, the calculation complexity of single eigenvalue decomposition and spectrum peak search of the original modified-MUSIC algorithm based on non-sparse frequency sampling is O{N 3 +n(N-K)(2(N-K+1))}, under the parameters given in the present embodiment, the calculation complexity is reduced by Since Therefore, it can be seen that the calculation complexity is significantly reduced.

[0139] (2) Estimation accuracy analysis:

[0140] In order to evaluate the performance of the present application, based on MATLAB, 100 times of Monte Carlo experiments are carried out under the parameters given in the above steps S1 to S9, and the RMSE of the estimation result is given in Figure 5 .

[0141] From Figure 5It can be seen that in the whole signal-to-noise ratio SNR=0:4:20dB range, the RMSE of time delay estimation of the original modified-MUSIC algorithm does not decrease obviously with the increase of signal-to-noise ratio, while the RMSE of the method proposed in the application is less than that of the original modified-MUSIC algorithm in the whole SNR range, and it can be seen that the RMSE of time delay estimation of the algorithm proposed in the application gradually decreases with the increase of signal-to-noise ratio, although the downward trend of the RMSE curve tends to be flat due to the limitation of the spectral peak search grid, but the RMSE at SNR=20dB is still nearly two orders of magnitude smaller than that of the original modified-MUSIC algorithm, and the accuracy of time delay estimation is obviously improved.

[0142] In summary, the embodiment can realize rough road echo time delay estimation with low computational complexity and high estimation accuracy, the method is simple to realize, has wide applicability and strong robustness, and can be used for ground penetrating radar echo signal time delay estimation of various road structures.

[0143] The above embodiment is a preferred embodiment of the application, but the embodiments of the application are not limited by the above embodiment, and any changes, modifications, substitutions, combinations, simplifications made without departing from the spirit and principles of the application shall be equivalent replacement methods, and all shall be included in the protection scope of the application.

Claims

1. A low-complexity rough road echo time delay estimation method based on sparse frequency sampling, characterized in that, The estimation method comprises the following steps: The transmitting end transmits a ground penetrating radar signal; The ground penetrating radar signal reaches the multi-layer rough road surface, is reflected by the interfaces of each layer, and forms multiple backscattering echoes; The receiving end detects the backscattering echoes to obtain initial observation data; The initial observation data is subjected to white processing, and the radar pulse is removed from the initial observation data to obtain receiving data; The receiving data is subjected to grouping and sparse frequency sampling to obtain multiple groups of low-dimensional receiving data sub-sampling sets; The covariance matrix of each group of sub-sampling sets is calculated; The covariance matrix of each group of sub-sampling sets is subjected to linear interpolation; Each covariance matrix after linear interpolation is respectively subtracted by a noise matrix; The covariance matrix after subtracting the noise matrix is subjected to spatial smoothing preprocessing; A modified MUSIC is applied to each covariance matrix after spatial smoothing preprocessing to obtain multiple time delay estimation spectra; The minimum time delay estimation spectrum period T is calculated; Spectrum peak searching is performed on each time delay estimation spectrum to obtain multiple groups of time delay estimation values; The time delay of the first echo is calculated, that is, cross-correlation operation is performed on the multiple time delay estimation spectra to obtain a cross-correlation spectrum, and the time when the maximum value of the cross-correlation spectrum is obtained is taken as the time delay value τ1 of the first echo; Each group of time delay estimation values obtained by spectrum peak searching is selected, and the time delay value belonging to the interval [τ1, τ1+T) in each group of time delay estimation values is selected respectively and recorded in an array; Row summation and averaging are performed on the receiving data to obtain a column vector X; The time delay values in the array are constructed into a complete dictionary D, and for a given number of echoes K, K correct time delay values are screened out through K iterations, wherein the iteration screening process is as follows: T1, initialize the iteration number i = 1, the initial residual b = X, the initial selected matrix S = [], the initial amplitude matrix a = [], and the initial time delay screening result t3 = []; T2, the contribution of each column vector in the complete dictionary D to the residual b is calculated by least square, the column vector with the maximum contribution value is selected and recorded as u, and the time delay value corresponding to u is recorded as τ; T3, update the amplitude matrix a = u\X, and the selected matrix S = [S u]; T4, update the residual b = X-aS, update t3 = [t3 τ], update i = i + 1, and delete the element u from the complete dictionary D; Repeat the operations of T2-T4 until i > K, end the screening, and output the time delay screening result t3 as an array containing K elements.

2. The low complexity rough road echo time delay estimation method based on sparse frequency sampling according to claim 1, characterized in that, The m integers are selected from the sparse frequency sampling of the group, and the i-th integer is represented as M i , i = 1…m, respectively, to obtain a frequency sampling interval of M i Δf, the i-th group of received data sub-sampling set r i , i = 1…m, wherein Δf represents a minimum frequency interval satisfying the Nyquist sampling theorem.

3. The low complexity rough road echo time delay estimation method based on sparse frequency sampling of claim 1, wherein, In the linear interpolation of the covariance matrix of each group of sub-sampling sets, the i-th group of sub-sampling sets r i The i-th group of covariance matrices is constructed according to the following formula The linear interpolation is performed, where the superscript H represents the complex conjugate transpose operation, to obtain m linearly interpolated covariance matrices, and the i-th linearly interpolated covariance matrix is specifically represented as: Wherein, B represents an interpolation transformation matrix.

4. The low complexity rough road echo time delay estimation method based on sparse frequency sampling of claim 1, wherein, In the noise matrix subtracted from each linear-interpolated sub-sampling set of covariance matrices, the noise variance of the i-th linear-interpolated covariance matrix is estimated In the noise matrix subtracted from each linear-interpolated sub-sampling set of covariance matrices, the noise variance of the i-th linear-interpolated covariance matrix is estimated The i-th noise-free covariance matrix is calculated as follows: where ∑ = Λ -1 Λ -H Λ represents the diagonal matrix of the frequency domain radar pulse.

5. The low-complexity rough road echo time delay estimation method based on sparse frequency sampling according to claim 4, characterized in that, The noise variance estimation is calculated by a propagation operator PM method, and the process is as follows: Given a data covariance matrix R xx is partitioned into the form: where G1, G2, H1and H2are (KxK), ((N-K)xK), (Kx(N-K)) and ((N-K)x(N-K)) dimensional sub-matrices, respectively, N and K denote the dimension of the data covariance matrix R xx and the number of echoes, respectively, the noise variance σ xx of the data covariance matrix R 2 is derived as follows: wherein tr{•} denotes the trace operator of a matrix, the superscript + denotes the pseudo-inverse of a matrix, I N-k denotes the (N-k) x (N-k) identity matrix.

6. The low complexity rough road echo time delay estimation method based on sparse frequency sampling of claim 1, wherein, In the spatial smoothing preprocessing, the ith noiseless covariance matrix With the following modification, the ith spatially smoothed preprocessed covariance matrix is given by wherein represents the lth subband of the ith noiseless covariance matrix i L is an integer greater than 1 and less than half the length.​ 7. The low complexity rough road echo time delay estimation method based on sparse frequency sampling of claim 2, wherein, The modified-MUSIC is modified-MUSIC of the ith spatially smoothed preprocessed covariance matrix Eigenvalue decomposition is performed to obtain the corresponding ith noise subspace U ni ; Define the dimension of the ith spatially smoothed preconditioned covariance matrix υ i = N i - L i + 1, where N i is the length of the sub-sampled set r i . Subsequently, the time domain search range t of the MUSIC algorithm is set as t = t1:Δt:t2, wherein t1 and t2 are the minimum and maximum time delays respectively, and Δt is a time step; By utilizing the orthogonality of the signal subspace and the noise subspace, a modified MUSIC time delay estimation spectrum is constructed, and the expression of the i th modified MUSIC time delay estimation spectrum is as follows: where λ ik (t),i = 1...m; k = 1,2...2K represent the kth smallest eigenvalue of , real{·} denotes the real part operation, and K is the number of echoes, denotes the steering vector of the ith spatially smoothed pre-processed covariance matrix of the diagonal matrix; Substituting the calculation of the i-th modified time delay estimation spectrum P i (t).

8. The low-complexity rough road echo time delay estimation method based on sparse frequency sampling according to claim 7, characterized in that, The minimum time delay estimation spectrum period is calculated as follows: For m integers that are prime to each other, the i-th integer is denoted as M i , i = 1...m, let M i be the largest integer among them, and let M max be the smallest integer among them. The minimum time delay estimation spectrum period is calculated as follows: T = (t2 - tl) / (2(M max )-1) Wherein, t1 and t2 are the minimum and maximum time delays respectively.

9. The sparse frequency-based low-complexity rough road echo time delay estimation method according to claim 1, characterized in that, The ground penetrating radar signal transmitted by the transmitting end is an ultra-wideband stepped frequency radar.

Images