A high-precision multipath parameter fast estimation method suitable for dynamic scenes

Abstract

The application provides a high-precision multipath parameter fast estimation method suitable for a dynamic scene, which can greatly reduce the number of time-delay-Doppler units to be estimated through selection of a concerned interval; a two-dimensional focusing result is iteratively filtered by using an adaptive filter, in each iteration, a complex amplitude estimation result obtained in the last iteration is used as prior information of the current iteration, the filter of the current iteration is adaptively updated, and then the complex amplitude of each unit is estimated; and an iteration convergence judgment method is used, so that the two-dimensional channel response of the time-delay-Doppler plane can be successfully reconstructed after a limited number of iterations, and the delay and power of each path are accurately estimated.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of direct sequence spread spectrum signal processing, and particularly relates to a high-precision multipath parameter fast estimation method suitable for dynamic scenes. BACKGROUND

[0002] In a direct sequence spread spectrum measurement system, multipath effect is one of the main factors affecting measurement accuracy. The classical anti-multipath technology can be divided into four categories: antenna design-based method, receiver tracking loop design-based method, measurement value processing-based method, and multipath parameter estimation-based method. The anti-multipath antenna can produce high gain within a certain angle range and suppress incoming wave signals outside the range, but generally has a large size, high cost, and is difficult to suppress multipath interference of incoming wave signals close to the direct path signal. The receiver tracking loop design-based method is easy to apply to the existing receiver structure, has good real-time performance, but is difficult to completely eliminate multipath error in a short delay multipath scene. The measurement value processing-based technology uses the statistical characteristics of multipath error to design a filter to process the measurement data of pseudorange and carrier phase, and this type of method depends on the accuracy of the used multipath model. The multipath parameter estimation-based method estimates the multipath parameters of the received signal using statistical signal processing, and reconstructs and cancels the multipath signal according to the estimation result, and has good adaptability to different multipath scenes, and is one of the mainstream research directions of the current anti-multipath technology. However, this type of algorithm is still limited by factors such as computational complexity, storage requirements, and multipath resolution.

[0003] In the article "Multipath mitigation technique under strong multipath environment using multiple antennas" published by Kubo Nobuski et al. in Journal of Aeronautics, Astronautics and Aviation, pages 75-82, 2017, the direct signal and reflected signal received by the receiver are used to suppress multipath by taking into account the different carrier phase variation relationships when the receiving platform is stationary and moving. This method requires the receiving platform to maintain movement and is not suitable for fixed receiving scenarios. In the article "Deep neural network correlators for GNSS multipath mitigation" published by Li Haoqing et al. in IEEE Transactions on Aerospace and Electronic Systems, vol. 59, no. 2, pages 1249-1259, 2022, a deep neural network correlator method is proposed. This method can learn complex multipath characteristics, but requires a multi-correlator structure and has high computational complexity. In the article "Multipath analysis using code-minus-carrier for dynamic testing of GNSS receivers" published by Blanco-Delgado Nuria et al. in 2011 International Conference on Localization and GNSS, pages 25-30, 2011, the code-minus-carrier measurement processing method is analyzed. This method uses carrier phase measurements to smooth code phase measurements and can reduce multipath errors in the presence of thermal noise. However, this method has a large dynamic tracking error in receivers without carrier assistance. SUMMARY

[0004] Therefore, the present application provides a high-precision fast multipath parameter estimation method suitable for dynamic scenarios, which can effectively reconstruct the time delay-Doppler two-dimensional channel response and quickly estimate the code phase delay and power of each path.

[0005] A high-precision fast multipath parameter estimation method suitable for dynamic scenarios, comprising:

[0006] Step 1: Establish a received signal model under a dynamic multipath scenario;

[0007] Step 2: after slow-time windowing, Doppler compensation and code correlation, the spectrum is weighted to focus the signal in two dimensions;

[0008] Step 3: according to the power relation between the two-dimensional focusing result obtained in Step 2 and the noise floor, a concerned interval is set;

[0009] Step 4: a cross window is applied to design an adaptive filter in the concerned interval, and the accumulated result after Step 2 is filtered to estimate the two-dimensional channel response;

[0010] Step 5: Step 4 is repeated, and the two-dimensional channel response estimation value obtained in each iteration is used as the prior knowledge for the next iteration to redesign the optimal filter for each time-delay-Doppler unit in the concerned interval and perform filtering to estimate the complex amplitude of each unit.

[0011] Preferably, the method for establishing the received signal model in the dynamic multipath scenario in Step 1 comprises:

[0012] Step 1: a received signal model in a dynamic multipath scenario is established, and the specific method is as follows:

[0013] Let the expression of the transmitted waveform of the direct spread system be:

[0014] s tran (t) = C(t) exp(j2πf c t) (1)

[0015] where C(·) is the spreading pseudo-code, the code length is set as L PNS , the chip width is T chip , the pseudo-code period is T PNC = L PNS T chip , f c is the carrier frequency; t represents the time variable;

[0016] In the dynamic multipath scenario, the baseband received signal of the receiver in a coherent accumulation time interval is represented as:

[0017]

[0018] where, is the fast time, t m = mT PNC is the slow time, m = 0, 1, 2,..., M-1, M is the number of code periods of coherent accumulation, γ is the signal path number, γ = 0 represents the direct signal, Γ is the number of signal paths, A γ is the signal amplitude, is the path phase shift, f d,γ is the Doppler frequency, τ γ = τ0+Δτγ τ0is the initial delay of the direct path, Δτ γ is the relative delay of the multipath, τ m,γ = τ γ + (f d,γ / f c )t m ; without loss of generality, assume that for any γ1< γ2, there exists is an additive complex Gaussian white noise;

[0019] The sampled discrete baseband received signal is denoted as:

[0020]

[0021] where n = 0, 1, 2,, N - 1 is the fast time discrete sampling sequence, N = T PNC f s , f s is the sampling frequency, denotes the floor function, υ(n, m) is the sampled discrete baseband noise.

[0022] Preferably, the step 2 of performing slow time windowing, Doppler compensation and code correlation on the received signal and then performing spectrum weighting processing is performed according to the following method:

[0023] The sampled discrete baseband received signal is windowed in the slow time dimension and is subjected to M' point discrete Fourier transform (DFT) so that the signal energy is focused in the Doppler dimension with low sidelobes, and the result is denoted as:

[0024]

[0025] where m f = -M' / 2, -M' / 2 + 1,, M' / 2 - 1, M' ≥ M is the number of Doppler units, which is set to an integer power of 2, w M (m) is a window function with a length of M, W M′ (m f ) is the M' point DFT spectrum function thereof, denotes the DFT with respect to the variable x, m f,γ = -f d,γ T' CPI , T' CPI = M'T PNC , υ w,f (n, m f ) is the result of performing slow time windowing and M' point DFT processing on the sampled discrete baseband noise;

[0026] The Doppler compensation function is constructed as:

[0027]

[0028] s w,f (n,m f ) is Doppler compensated, and the result is denoted as:

[0029]

[0030] wherein υ dc (n,m f ) is Doppler compensated υ w,f (n,m f );

[0031] s dc (n,m f ) is fast-time dimension matched filtered, and the result is denoted as:

[0032]

[0033] wherein (·) * denotes complex conjugate:

[0034]

[0035] r C (n) is the pseudo-code autocorrelation function of C(n), and the expression of a(m f ) is:

[0036]

[0037] υ cc (n,m f ) is the result of fast-time dimension matched filtering υ dc (n,m f );

[0038] s cc (n,m f ) is approximately:

[0039]

[0040] wherein r F (m f ) = a(m f ) W M′ (m f ), * denotes convolution operation, and x(n,m f ) is the two-dimensional channel response in time-delay-Doppler plane, and the expression is:

[0041]

[0042] is a two-dimensional discrete Dirac impulse function;

[0043] The target waveform after spectrum weighting is defined as:

[0044]

[0045] where N chip = T chip f s is the average number of sampling points per chip, is a standard triangle wave function; the spectrum weighting coefficient function is constructed as:

[0046]

[0047] where n f is a discrete frequency point; the spectrum weighting is performed on s cc (n, m f ) to obtain:

[0048]

[0049] where denotes the inverse discrete Fourier transform with respect to the variable x, and υ sw (n, m f ) is the result obtained by performing spectrum weighting on υ cc (n, m f );

[0050] After the two-dimensional focusing processing of the spectrum weighting fusion after the slow-time windowing, Doppler compensation and code correlation, the two-dimensional point spread functions of each path signal component in the received signal are all unified as r sw (n) r F (m f ).

[0051] Preferably, in step 3, the method for setting the concerned interval comprises:

[0052] Let the power of the noise υ(n, m) in the baseband received signal be σ 2 , and the theoretical noise power of the signal s cc (n, m f ) after the slow-time windowing, Doppler compensation and code correlation processing be Let β dete1 denote the constant false alarm detection threshold corresponding to the false alarm rate P FA , and β dete2 denote the first peak sidelobe power of the highest power in the signal s cc (n, m f ), i.e.

[0053]

[0054] Where PSR represents signal s cc (n,m f The point spread function r C (n)r F (m f The peak-to-sidelobe ratio (PSR) is defined as the ratio of the power of the first peak sidelobe to the peak intensity of the main lobe; when the constant false alarm rate (CFAR) detection threshold is higher than the signal s... cc (n,m f When the first peak sidelobe power reaches its highest power, the final pre-detection threshold is set to the constant false alarm rate (CFAR) detection threshold; otherwise, it is set to the signal s. cc (n,m f The first peak sidelobe power of the highest power, i.e.

[0055] β dete =max{β dete1 ,β dete2} (14)

[0056] Statistics cc (n,m f The time-delayed Doppler cells whose power exceeds the pre-detection threshold are denoted as:

[0057]

[0058] right For each point in the time delay dimension and Doppler dimension, its range is extended by several points, and the extended set of points is used as the final region of interest, as shown below:

[0059]

[0060] Preferably, the specific methods for steps 4 and 5 include:

[0061] The time-delay-Doppler two-dimensional signal obtained after spectral weighting is written in matrix form:

[0062] S sw =R sw XR F +Υ sw (17)

[0063] Where S sw R sw X, R F and Υ sw The definition is as follows:

[0064] S sw =[s sw (-M′ / 2) s sw(-M' / 2 + 1)... s sw (M' / 2 - 1)]

[0065] s sw f ) = [s sw (0, m f ) s sw (1, m f )... s sw (N-1, m f )] T

[0066]

[0067] X = [x(-M' / 2) x(-M' / 2 + 1)... x(M' / 2 - 1)]

[0068] x(m f ) = [x(0, m f ) x(1, m f )... x(N-1, m f )] T

[0069]

[0070] Y sw = [y sw (-M' / 2) y sw (-M' / 2 + 1)... y sw (M' / 2 - 1)]

[0071] y sw (m f ) = [y sw (0, m f ) y sw (1, m f )... y sw (N-1, m f )] T

[0072] where (·) T denotes the matrix transpose;

[0073] For a delay-Doppler cell (n, m f ), the expression of the data vector to be filtered is defined as:

[0074]

[0075] where

[0076] ​

[0077] δ a；(b) = [δ(-b) δ(-b+1)... δ(-b+a-1)] T , is the Dirac impulse function, K r and K d are filter size factors in time delay and Doppler dimensions, B1, B 2,n , The expressions of B4 are respectively:

[0078]

[0079]

[0080] I a denotes a x a identity matrix; substituting S sw into the equation yields:

[0081]

[0082] where,

[0083] Since r sw (n) is a periodic function, there are When the Doppler frequency of the signal path and the filter Doppler dimension factor satisfy |f d,γ |≤1 / (4T PNC ) and K d <M′ / 4, there are and At this time and are simplified as:

[0084]

[0085] where

[0086] Let k d be the waveform approximation coefficient of the Doppler dimension, and denotes the upward rounding; the matrices B5 and B6 are constructed as follows:

[0087]

[0088] When |n|≥N chip , r sw (n) = 0, there are and In addition, after windowing and accumulation in the slow time dimension, most of the energy of the Doppler dimension is focused in the main lobe, there are and In conjunction with the properties of the vectorization operator, With Simplify as:

[0089]

[0090] Where

[0091] The data vector to be filtered is represented as:

[0092]

[0093] Where

[0094]

[0095] The adaptive filter design is based on the minimum mean square error criterion, and its cost function is:

[0096]

[0097] Where, (·) H is the matrix conjugate transpose operation;

[0098] In order to minimize the cost function, the partial derivative of J(n,m f ) with respect to is solved and set to 0, so that the optimal filter coefficient vector is:

[0099]

[0100] Where

[0101]

[0102] (·) -1 represents the matrix inversion, (·) * represents the complex conjugate operation, and E(·) represents the statistical expectation;

[0103] Assuming that the responses between different time delay-Doppler units are independent of each other, and the channel response and the noise are independent of each other, then:

[0104]

[0105] Where is a diagonal matrix with the main diagonal vector as , ρ(n,m f )=E[|x(n,m f )| 2 ], d[:,a] is the a-th column of D, diag{a} represents the diagonal matrix with a as the diagonal line, and is the Hadamard product; Rυ is the noise covariance matrix, denoted as:

[0106]

[0107] For each delay-Doppler cell in the concerned interval, the calculated optimal filter coefficient vector is used to filter the constructed data vector to be filtered, i.e. traversing to obtain the estimation of the channel response in the concerned interval, thereby completing one iteration;

[0108] Step 5: Repeat step 4, and after iteration convergence, the reconstructed delay-Doppler plane two-dimensional channel response estimation and the code phase delay and power of each path are obtained.

[0109] Preferably, in step 5, after each iteration, the unit whose channel response estimation result power is higher than the theoretical noise power after step 2 processing is counted, i.e.

[0110]

[0111] wherein

[0112]

[0113] Definition

[0114]

[0115] wherein, N RoI is the total number of units in the concerned interval, N valid and N valid,lag are the number of units in the current iteration round and the last iteration round set When ΔN valid,norm is less than the preset convergence factor β iter , the iteration filtering algorithm converges, thereby effectively reconstructing the delay-Doppler plane two-dimensional channel response.

[0116] Preferably, in the first iteration round, the processing result s sw (n, m f ) obtained in step 2 is used as the prior knowledge of the channel response to design the optimal filter.

[0117] The present application has the following beneficial effects:

[0118] The present application proposes a high-precision multipath parameter fast estimation method suitable for dynamic scenes, which can solve the problem of high-precision multipath parameter estimation in dynamic scenes, and can realize fast estimation of signal path response by selecting a concerned interval and designing an adaptive filter using a cross window.

[0119] By using the slow-time windowing, the signal waveform sidelobe of the Doppler velocity can be effectively reduced, and most of the signal energy can be concentrated in the main lobe; by using the Doppler compensation, the coupling between the signal waveform of each path and the path dynamics can be removed, facilitating the spectral weighting processing and the adaptive filter design; by using the spectral weighting, the width of the signal time domain waveform can be compressed from the entire pseudo code period to several pseudo code chips, facilitating the subsequent adaptive filtering to quickly reconstruct the channel response; by selecting the attention interval, the number of time delay-Doppler units to be estimated can be greatly reduced; by using the adaptive filter to iteratively filter the two-dimensional focusing result, in each iteration, the complex amplitude estimation result obtained in the last iteration is used as the prior information of the current iteration, the filter of the current iteration is adaptively updated, and then the complex amplitude of each unit is estimated; by using the iterative convergence judgment method, the two-dimensional channel response of the time delay-Doppler plane can be successfully reconstructed after a limited number of iterations, and the delay and power of each path can be accurately estimated. BRIEF DESCRIPTION OF DRAWINGS

[0120] Figure 1 The two-dimensional focusing result (normalized) of the slow-time windowing, the Doppler compensation and the code correlation and then the spectral weighting fusion processing.

[0121] Figure 2 The iterative filtering result (normalized). DETAILED DESCRIPTION

[0122] The application provides a high-precision multipath parameter fast estimation method suitable for a dynamic scene. The application is described in detail below in combination with embodiments.

[0123] Step 1: Establish a received signal model in a dynamic multipath scene, and the specific method is as follows:

[0124] Let the expression of the transmission waveform of the direct spread system be

[0125] s tran (t) = C(t)exp(j2πf c t) (1)

[0126] Wherein C(·) is a spread spectrum pseudo code, the code length of which is L PNS , the chip width is T chip , the pseudo code period is T PNC = L PNS T chip , and f c is a carrier frequency.

[0127] Under the dynamic multipath scene, the baseband received signal of the receiver in one coherent accumulation time interval can be expressed as

[0128]

[0129] where, is fast time, t m = mT PNC is slow time, m = 0, 1, 2,..., M - 1, M is the number of code periods for coherent integration, γ is the signal path index, γ = 0 represents the direct path, Γ is the number of signal paths, A γ is the signal amplitude, is the path phase shift, f d,γ is the Doppler frequency, τ γ = τ0+ Δτ γ is the initial delay of the path, τ0is the initial delay of the direct path, Δτ γ is the relative delay of the path, τ m,γ = τ γ + (f d,γ / f c ) t m Without loss of generality, it is assumed that for any γ1< γ2, there exists τ γ1 ≤ τ γ2 , is the additive complex white Gaussian noise.

[0130] Since the Doppler frequency is much smaller than the carrier frequency, the effect of code Doppler on the discrete samples can be ignored, and the sampled discrete baseband received signal can be expressed as

[0131]

[0132] where n = 0, 1, 2,..., N - 1 is the fast time discrete sampling sequence, N = T PNC f s , f s is the sampling frequency, represents the floor function, υ(n, m) is the sampled discrete baseband noise.

[0133] Step 2: After slow time windowing, Doppler compensation, and code correlation on the received signal, the spectrum is weighted, and the specific method is as follows:

[0134] The sampled discrete baseband received signal is windowed in the slow time dimension and subjected to M' point Discrete Fourier Transform (DFT), so that the signal energy is focused in the Doppler dimension with low sidelobes. The result can be expressed as

[0135]

[0136] where m f= -M' / 2, -M' / 2 + 1,..., M' / 2 - 1, M' ≥ M is the number of Doppler bins, usually set as an integer power of 2 so that the DFT operation can be implemented using fast Fourier transform, w M (m) is a window function with length M, W M′ (m f ) is its M' point DFT spectrum function, represents the DFT with respect to variable x, m f,γ = -f d,γ T' CPI , T' CPI = M' T PNC , u w,f (n, m f ) is the result of slow time windowing and M' point DFT processing on the sampled discrete baseband noise.

[0137] The Doppler compensation function is constructed as

[0138]

[0139] The Doppler compensation is performed on s w,f (n, m f ), and the result can be represented as

[0140]

[0141] wherein u dc (n, m f ) is the result of Doppler compensation on u w,f (n, m f ).

[0142] The fast time dimension matched filtering is performed on s dc (n, m f ), and the signal after code correlation in time delay dimension can be represented as

[0143]

[0144] wherein (·) * represents complex conjugate

[0145]

[0146] r C (n) is the pseudo code autocorrelation function of C(n), and the expression of a(m f ) is

[0147]

[0148] u cc (n, m f ) is the result of Doppler compensation on udc (n,m f The result obtained by performing fast time dimension matched filtering.

[0149] s cc (n,m f It can be approximated as

[0150]

[0151] Where r F (m f )=α(m f W M′ (m f ), * indicates convolution operation, x(n,m) f The expression for the two-dimensional channel response on the time-delay-Doppler plane is as follows:

[0152]

[0153] It is a two-dimensional discrete Dirac impulse function.

[0154] Define the target waveform after spectral weighting as

[0155]

[0156] Where N chip =T chip f s It is the average number of sampling points per chip. The standard trigonometric wave function is used. The spectral weighting coefficient function is constructed as follows:

[0157]

[0158] in nf represents discrete frequency points. For s cc (n,m f ) By performing spectral weighting, we can obtain

[0159]

[0160] in υ represents the discrete inverse Fourier transform of the variable x. sw (n,m f ) for υ cc (n,m f The result is obtained by spectral weighting.

[0161] After two-dimensional focusing processing involving slow-time windowing, Doppler compensation, code correlation, and spectral weighted fusion, the two-dimensional point spread function of each path signal component in the received signal is unified to r. sw (n)r F (mf ).

[0162] Step 3: To reduce the computational load of subsequent adaptive filtering, the interested interval is set according to the power relationship between the two-dimensional focusing result after slow-time windowing, Doppler compensation and code correlation processing and the noise floor. The specific steps are as follows:

[0163] Let the power of the noise υ(n, m) in the baseband received signal be σ 2 , and the theoretical noise power of the signal s cc (n, m f ) after slow-time windowing, Doppler compensation and code correlation processing is Let β dete1 represent the constant false alarm detection threshold corresponding to the false alarm rate P FA , and β dete2 represent the first peak sidelobe power of the highest power in the signal s cc (n, m f ), that is,

[0164]

[0165] where PSR represents the peak sidelobe ratio (Peak Sidelobe Ratio, PSR) of the point spread function r cc (n)r f (m C ) of the signal s F (n, m f ), which is defined as the ratio of the first peak sidelobe power to the main lobe peak intensity. When the constant false alarm detection threshold is higher than the first peak sidelobe power of the highest power in the signal s cc (n, m f ), the final pre-detection threshold is set to the constant false alarm detection threshold, otherwise it is set to the first peak sidelobe power of the highest power in the signal s cc (n, m f ), that is,

[0166] β dete = max{β dete1 , β dete2} (14)

[0167] The time-delay-Doppler cells whose power exceeds the pre-detection threshold in s cc (n, m f ) are counted, and the set is denoted as:

[0168]

[0169] For each point in , its range is extended by several points in both time-delay and Doppler dimensions, and the extended point set is taken as the final interested interval, as shown below:

[0170]

[0171] Step 4: Apply the cross window design adaptive filter in the concerned interval to filter the results of the windowing, Doppler compensation and code correlation after spectrum weighting in the concerned interval, so as to estimate the two-dimensional channel response, the specific method is as follows:

[0172] Write the time delay-Doppler two-dimensional signal obtained after spectrum weighting as a matrix form

[0173] S sw = R sw X R F + Y sw (17)

[0174] Wherein S sw , R sw , X, R F and Y sw are defined as follows

[0175] S sw = [s sw (-M' / 2) s sw (-M' / 2+1) … s sw (M' / 2-1)]

[0176] s sw (m f ) = [s sw (0,m f ) s sw (1,m f ) … s sw (N-1,m f )] T

[0177]

[0178] X = [x(-M' / 2) x(-M' / 2+1) … x(M' / 2-1)]

[0179] x(m f ) = [x(0,m f ) x(1,m f ) … x(N-1,m f )] T

[0180]

[0181] Y sw = [υ sw (-M' / 2) υ sw(-M' / 2 + 1)... u sw (M' / 2 - 1)

[0182] u sw (m f ) = [u sw (0,m f ) u sw (1,m f )... u sw (N-1,m f )] T

[0183] where (·) T denotes matrix transposition.

[0184] For the delay-Doppler cell (n,m f ), the expression of the data vector to be filtered is defined as

[0185]

[0186] where

[0187]

[0188] δ a；(b) = [δ(-b) δ(-b+1)... δ(-b+a-1)] T , is the Dirac impulse function, K r and K d are filter size factors in delay and Doppler dimensions, respectively, B1, B 2,n , The expressions of B4are

[0189]

[0190] I a denotes a a identity matrix. Substituting S sw into the above equation, we have

[0191]

[0192] where,

[0193] Since r sw (n) is a periodic function, we have When the Doppler frequency of the signal path and the filter Doppler dimension size factor satisfy |f d,γ |≤1 / (4T PNC ) and K d <M' / 4, we have and In this case, with can be simplified as

[0194]

[0195] where

[0196] Note that k d is the waveform approximation coefficient of Doppler velocity, denotes the upward rounding. The matrices B5 and B6 are constructed as follows

[0197]

[0198] When |n|≥N chip , r sw (n) = 0, thus we have and In addition, after the windowed accumulation in the slow time dimension, most of the energy of Doppler velocity is focused in the main lobe, thus we have and Combining the properties of the vectorization operator, with can be simplified as

[0199]

[0200] where

[0201] The data vector to be filtered can be represented as

[0202]

[0203] where

[0204]

[0205]

[0206] The adaptive filter design is based on the minimum mean square error criterion, and its cost function is

[0207]

[0208] where, (·) H is the matrix conjugate transpose operation.

[0209] In order to minimize the cost function, we need to solve the partial derivative of J(n, m f ) with respect to and let it equal to 0, thus we get the optimal filter coefficient vector as

[0210]

[0211] where

[0212]

[0213] (·) -1 denotes matrix inversion, (·) * denotes complex conjugate operation, E(·) represents statistical expectation.

[0214] Assuming that the responses of different time-delay-Doppler cells are independent of each other, and the channel response and the noise are independent of each other, then

[0215]

[0216] where is the main diagonal vector, and D is a diagonal matrix with ρ(n,m f )=E[|x(n,m f )| 2 ], d[:,a] is the a-th column of D, diag{a} represents a diagonal matrix with a as the diagonal line, and is the Hadamard product. R υ is the noise covariance matrix, which can be expressed as

[0217]

[0218] For each time-delay-Doppler cell in the concerned interval, the calculated optimal filter coefficient vector is used to filter the constructed to-be-filtered data vector, that is, After traversing , the estimated value of the channel response in the concerned interval can be obtained, thereby completing one iteration.

[0219] Step 5: Repeat step 4, and the channel response estimated value obtained in each round of iteration is used as the prior knowledge in the next round of iteration, and the optimal filter is redesigned and filtered for each time-delay-Doppler cell (n,m f ) in the concerned interval. After each round of iteration, the cells whose channel response estimation results have a power higher than the theoretical noise power after the processing in step 2 are counted, that is,

[0220]

[0221] where

[0222]

[0223] Definition

[0224] ​

[0225] where N RoI is the total number of cells in the concerned interval, N valid and N valid,lag are the number of cells in the current iteration and the last iteration set respectively. When ΔN valid,norm is less than a preset convergence factor β iter , the iteration filtering algorithm converges, thus effectively reconstructing the two-dimensional channel response of the delay-Doppler plane.

[0226] In the first iteration, the processing result s sw (n,m f ) obtained in step 2 is used as the prior knowledge of the channel response to design the optimal filter, and step 4 is repeated. After the iteration converges, the reconstructed delay-Doppler plane two-dimensional channel response estimate can be effectively obtained, and the code phase delay and power of each path can be estimated.

[0227] Embodiment:

[0228] In this example, the carrier frequency is f c = 1.57542 GHz, the sampling rate is f s = 1 MHz, the spreading code uses an m sequence with a code length of L PNS = 127 and a code period of T PNC = 1 ms, the accumulation code period number is M = 16, the Doppler cell number is M' = 32, the Hamming window weighting is used, the number of signal paths in the scenario is Γ = 3, the amplitude of the direct path signal is A0 = 1, the signal multipath ratios of the two multipaths are -6.3 dB and -3.7 dB respectively, the path phase shifts are the code phase delays are τ0 = 100 chips, τ1 = 100.2 chips, and τ2 = 101.2 chips, the Doppler frequencies are f d,0 = -31.25 Hz, f d,1 = 0, and f d,2 = 31.25 Hz, the accumulation signal-to-noise ratio of the scenario is SNR in = -15 dB, the false alarm rate used when setting the concerned interval is P FA = 10 -5 , the filter size is K r = 3 and K d = 1, the Doppler waveform approximation coefficient is k d = 4, and the iteration convergence factor is β iter = 0.01.

[0229] ​Through the above five steps, the distance-Doppler unit number in the attention interval is 414, and the filtering algorithm converges after four iterations, which can effectively reconstruct the two-dimensional channel response of the estimated delay-Doppler plane. The code phase delay estimation values of the three paths are 99.95 chips, 100.20 chips and 101.22 chips, respectively, the Doppler frequency estimation values are -31.25 Hz, 0 and 31.25 Hz, respectively, and the power estimation values are -0.07 dB, 6.95 dB and 3.52 dB, respectively.

[0230] In summary, the above is only a preferred embodiment of the present application, and is not used to limit the protection scope of the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A high-precision multipath parameter fast estimation method suitable for dynamic scenes, characterized in that, The method comprises the following steps: Step 1: establishing a received signal model under a dynamic multipath scenario; Step 2: performing slow-time windowing, Doppler compensation and code correlation on the received signal, and then performing spectrum weighting processing to focus the signal in two dimensions; Step 3: setting a concerned interval according to a power relationship between a two-dimensional focusing result obtained in step 2 and a noise base; Step 4: applying a cross window to design an adaptive filter in the concerned interval, and filtering the accumulated result processed in step 2, so as to estimate a two-dimensional channel response; Step 5: repeating step 4, and taking the two-dimensional channel response estimation value obtained in each iteration as prior knowledge for the next iteration, re-designing an optimal filter for each time-delay-Doppler unit in the concerned interval and performing filtering, and then estimating the complex amplitude of each unit. The method for establishing the received signal model under the dynamic multipath scenario in step 1 comprises the following steps: Step 1: establishing a received signal model under a dynamic multipath scenario, and the specific method is as follows: Supposing that a transmission waveform expression of the direct spread system is: s tran (t) = C(t) exp(j2πf c t) (1) where C(·) is a spreading pseudo-code with a code length set to L PNS , a chip width of T chip , and a pseudo-code period of T PNC = L PNS T chip , f c is the carrier frequency; t denotes the time variable; Under the dynamic multipath scenario, a baseband received signal of the receiver in a coherent accumulation time interval is expressed as: where, is the fast time, t m = mT PNC is the slow time, m = 0, 1, 2, …, M - 1, M is the number of code periods for coherent integration, γ is the signal path index, γ = 0 represents the direct path, Γ is the number of signal paths, A γ is the signal amplitude, is the path phase shift, f d,γ is the Doppler frequency, τ γ = τ0+△τ γ is the initial delay of the path, τ0is the initial delay of the direct path,△τ γ is the relative delay of the path, τ m,γ = τ γ + (f d,γ / f c ) t m ; without loss of generality, assume that for any γ1< γ2, there exists is an additive complex Gaussian white noise; A discrete baseband received signal after sampling is expressed as: where n = 0, 1, 2,..., N - 1 is a fast time discrete sampling sequence, N = T PNC f s , f s is a sampling frequency, denotes a floor function, υ(n, m) is a sampled discrete baseband noise; The specific method for performing the slow-time windowing, Doppler compensation and code correlation on the received signal in step 2 is as follows: The discrete baseband received signal after sampling is windowed in the slow-time dimension and is subjected to M' point discrete Fourier transform (DFT), so that the signal energy is focused in the Doppler dimension with a lower sidelobe, and the result is expressed as: wherein m f = -M' / 2, -M' / 2 + 1,..., M' / 2 - 1, M' ≥ M is the number of Doppler units, set as an integer power of 2, w M (m) is a window function of length M, W M′ (m f ) is its M' point DFT spectrum function, denotes the DFT with respect to the variable x, m f,γ = -f d,γ T' CPI , T' CPI = M'T PNC , υ w,f (n, m f ) is the result of slow time windowing and M' point DFT processing on the sampled discrete baseband noise; A Doppler compensation function is constructed as: s w,f (n,m f ) are Doppler compensated, the result is represented as: wherein υ dc (n,m f ) is the result of Doppler compensation of υ w,f (n,m f ). s dc (n,m f ) are matched in the fast time dimension, and the signal after time delay dimension code correlation is represented as: wherein (·) * denotes complex conjugation: r C (n) is the autocorrelation function of the pseudo-code C(n), a(m f ) is given by the expression: υ cc (n,m f ) is the result of fast time dimension matched filtering on υ dc (n,m f ). s cc (n,m f ) is approximately: where r F (m f ) = a(m f )W M′ (m f ), * denotes convolution operation, and x(n, m f ) is a two-dimensional channel response in the delay-Doppler plane, which is expressed as It is a two-dimensional discrete Dirac impulse function; A target waveform after spectrum weighting is defined as: where N chip = T chip f s is the average number of sample points per chip, is a standard triangle wave function; the spectral weighting coefficient function is constructed as wherein n f is a discrete frequency point; and cc (n,m f ) is spectrally weighted to obtain: where ID x FT[·] denotes the inverse discrete Fourier transform with respect to the variable x, u sw (n,m f ) is the result of a spectral weighting of u cc (n,m f ). After the two-dimensional focusing processing of the slow-time windowing, Doppler compensation, code correlation and spectrum weighting fusion, the two-dimensional point spread functions of each path signal component in the received signal are unified as r sw (n)r F (m f ); The method for setting the concerned interval in step 3 comprises the following steps: Let the power of the noise υ(n,m) in the baseband received signal be σ. 2 The signal s after slow-time windowing, Doppler compensation, and code correlation processing cc (n,m f The theoretical noise power is Let β dete1 P represents the false alarm rate. FA The corresponding constant false alarm rate (CFAR) threshold, β dete2 Indicates signal s cc (n,m f The first peak sidelobe power of the highest power in ) is, i.e. where PSR represents the peak-to-sidelobe ratio of the signal s cc (n,m f ) of the point spread function r C (n)r F (m f ) is defined as the ratio of the first sidelobe power to the main lobe peak intensity; when the constant false alarm detection threshold is higher than the first sidelobe power of the highest power of the signal s cc (n,m f ), the final pre-detection threshold is set to the constant false alarm detection threshold, otherwise it is set to the first sidelobe power of the highest power of the signal s cc (n,m f ). β dete = max{β dete1 , β dete2} (14) Statistics s cc (n,m f ) in which the power exceeds a pre-detection threshold, and let this set be denoted by: For each point in , the time delay and Doppler dimension are both extended by several points, and the extended point set is taken as the final range of interest, as follows:

2. A high-precision multipath parameter fast estimation method suitable for dynamic scenes according to claim 1, characterized in that, The specific method of step 4 and step 5 comprises the following steps: The time-delay-Doppler two-dimensional signal obtained after the spectrum weighting is written in a matrix form as: S sw = R sw XR F + Y sw (17) wherein S sw , R sw , X, R F and Y sw are defined as follows: S sw = [s sw (-M' / 2) s sw (-M' / 2+1) … s sw (M' / 2-1)] s sw (m f ) = [s sw (0,m f ) s sw (1,m f ) … s sw (N-1,m f )] T X=[x(-M' / 2) x(-M' / 2+1) … x(M' / 2-1)] x(m f ) = [x(0,m f ) x(1,m f )... x(N-1,m f )] T Y sw = [υ sw (-M' / 2) υ sw (-M' / 2+1) … υ sw (M' / 2-1)] υ sw (m f )=[υ sw (0,m f ) υ sw (1,m f ) … υ sw (N-1,m f )] T wherein (·) T denotes matrix transposition; For the delay-Doppler unit (n, m f ), the expression of the data vector to be filtered is defined as: Wherein δ a；(b) = [δ(-b) δ(-b+1)... δ(-b+a-1)] T , is the Dirac impulse function, K r and K d are filter size factors in the delay and Doppler dimensions, B1, B 2,n , The expression for B4 is given by I a denotes the a x a identity matrix; and substituting S sw yields: wherein Since r sw (n) is a periodic function, we have When the Doppler frequency of the signal path and the filter Doppler size factor satisfy |f d,γ |≤1 / (4T PNC ) and K d <M′ / 4, we have and In this case and are simplified as: wherein Note k d is a waveform approximation coefficient of the Doppler, denotes a rounding up; the matrices B5 and B6 are constructed as follows: When |n|≥N chip r sw (n) = 0, there and In addition, after windowing and accumulation in the slow time dimension, most of the energy of the Doppler bins is focused in the main lobe, and there and Combining the properties of the vectorization operator, and is simplified to: wherein The data vector to be filtered is expressed as: Wherein The adaptive filter design is based on a minimum mean square error criterion, and a cost function is: wherein (·) H is the matrix conjugate transpose operation; To minimize the cost function, the partial derivative of J(n, m f ) with respect to is taken and set to zero, resulting in the optimal filter coefficient vector: Wherein (·) -1 denotes matrix inversion, (·) * denotes complex conjugation, E(·) represents statistical expectation; Supposing that the responses of different time-delay-Doppler units are independent of each other, and the channel response is independent of the noise, then: in The main diagonal vector is a diagonal matrix, ρ(n,m f )=E[|x(n,m f )| 2 ], d[:,a] is the a-th column of D, diag{a} represents a diagonal matrix with a as the diagonal, and ⊙ is the Hadamard product; R υ The noise covariance matrix is ​​expressed as: For each time-delay-Doppler cell in the concerned interval, the calculated optimal filter coefficient vector is used to filter the constructed data vector to be filtered, i.e. traversing to obtain the estimation of the channel response in the concerned interval, thereby completing one iteration; Step 5: repeating step 4, and after iterative convergence, a reconstructed time-delay-Doppler two-dimensional channel response estimation can be obtained, and the code phase delay and power of each path are estimated.

3. A high-precision multipath parameter fast estimation method suitable for dynamic scenes according to claim 2, characterized in that, In step 5, after each iteration is ended, the units whose channel response estimation result power is higher than the theoretical noise power after step 2 processing are counted, that is: Wherein It is defined as where N RoI is the total number of cells in the interval of interest, N valid and N valid,lag are the number of cells in the current iteration and the previous iteration set respectively; when ΔN valid,norm is less than a preset convergence factor β iter , the iterative filtering algorithm converges, thereby effectively reconstructing the two-dimensional channel response of the delay-Doppler plane.

4. The method for fast and high-precision multipath parameter estimation for dynamic scenarios of claim 1, wherein, In the first round of iteration, the processed result s sw (n,m f ) obtained in step 2 is used as the prior knowledge of channel response to design the optimal filter.

Images