A low-complexity root-music-based drone DOA estimation method

Abstract

The application discloses a low-complexity Root-MUSIC based UAV DOA estimation method, and belongs to the field of array signal processing; specifically, first, a uniform linear array (ULA) signal model is built, single-snapshot data vector is calculated, FFT transformation is performed on received multiple-snapshot sample data, an average power spectrum diagram is solved, and a target DOA approximate interval is determined according to the number of peak points and corresponding angle information. Then, the determined DOA approximate interval is divided into several parts, an iteration initial value in each part interval is solved according to a principle that an average distance from a unit circle arc is the shortest, and the iteration initial value is updated and sorted; according to a traditional Root-MUSIC algorithm, multiple-snapshot data is subjected to forward and backward smoothing decorrelation, a sample covariance matrix is solved, and characteristic decomposition is performed, so that a noise subspace is solved; then, polynomial coefficients are solved by using the noise subspace. Finally, the polynomial is solved, and a signal source DOA is obtained; the application narrows a selection interval of a root iteration initial value, accelerates an algorithm convergence speed, and reduces an algorithm complexity.

Description

Technical Field

[0001] This invention belongs to the field of array signal processing, specifically relating to a method for estimating the direction-of-arrival (DOA) of unmanned aerial vehicles (UAVs) based on the low-complexity Root-MUSIC (Root Multiple Signal Classification) algorithm. Background Technology

[0002] Unmanned aerial vehicles (UAVs) have broad application prospects in military operations, aerial photography and videography, surveying and mapping, search and rescue, environmental monitoring, and precision agriculture, among other military and civilian industries. Beyond their beneficial uses, the hazards and safety risks posed by UAVs have gradually attracted public attention. Uncontrolled and undetected UAVs can pose serious threats to military operations, air traffic, human life, and privacy. Therefore, DOA (Depth of Account) estimation for UAVs has significant social value.

[0003] In applications involving real-time monitoring of drones, high latency requirements are necessary, while devices such as FPGAs (Field-Programmable Gate Arrays) have limited computing power. Therefore, researching low-complexity DOA estimation algorithms is of great significance.

[0004] DOA estimation algorithms for super-resolution, represented by Root-MUSIC algorithms in the subspace class, can overcome the Rayleigh limit. However, due to their large computational load, they are difficult to implement in engineering. The computational load is mainly concentrated in two relatively complex computational units: eigenvalue decomposition of the covariance matrix and polynomial root finding.

[0005] For the first unit, existing algorithms mostly reduce the complexity of feature decomposition by addressing real-valued processing, dimensionality reduction, and avoiding feature decomposition, but this results in a significant loss of accuracy. For the second unit, existing improved algorithms can reduce algorithm complexity while maintaining accuracy, but there is still room for improvement in terms of complexity. Summary of the Invention

[0006] To address the high complexity and engineering difficulty of the classic Root-MUSIC algorithm for DOA estimation, this invention proposes a low-complexity Root-MUSIC-based method for UAV DOA estimation. 1) By combining FFT (Fast Fourier Transform) operations and the LADTUC (Least Average Distance to Unit Circle) principle, the optimal initial value for root-finding iterations and the initial value update strategy are determined, narrowing the initial value selection interval and avoiding randomness in initial value updates; 2) The Muller root-finding algorithm based on factor cancellation is used to prevent the initial values ​​from converging to the same root, reducing the number of initial value updates; 3) A preset threshold is used to determine the authenticity of the roots. When the number of true roots equals the number of sources, the algorithm terminates, avoiding the need to find all roots before selecting true roots. Ultimately, this achieves low complexity.

[0007] The UAV DOA estimation method based on low-complexity Root-MUSIC described in this invention includes the following steps:

[0008] Step 1: Build a uniform linear array ULA signal model. A far-field narrowband signal hits the array, and the data vector of a single snapshot received by the array is calculated.

[0009] Data vector x of the i-th snapshot i (t)

[0010] x i (t)=A(θ)s i (t)+n i (t)

[0011] s i (t) represents a D×1 dimensional signal vector, where D is the number of far-field narrowband signals;

[0012] n i (t) represents M×1 dimensional white noise, where M is the number of elements in a uniform linear array;

[0013] A(θ)=[a(θ1),…,a(θ j ),…,a(θ D )] is the M×D dimensional array steering matrix; θ j It is the angle between the incident direction of the far-field narrowband signal and the normal direction of the array.

[0014] Step 2: Using the l-th snapshot data vector, perform a 2M-point Discrete Fourier Transform (DFT) to calculate the average power spectrum P of the l-th snapshot data. m (k), and normalize it;

[0015] Average power spectrum Pm (k):

[0016]

[0017] X i (k) is the data vector x of the i-th (i = 0, 1, ..., l-1) snapshot. i (t) Perform a 2M-point DFT to obtain the frequency domain sampled signal;

[0018]

[0019] x i (n) represents the received data of the nth array element in the i-th snapshot;

[0020] Step 3: Count all peak points of the normalized power spectrum, filter and remove peak points with normalized power values ​​less than -3dB, and calculate the frequency of each peak point after filtering.

[0021] The frequency of the i'th peak point after filtering is

[0022] N is the number of peak points whose normalized power value is greater than or equal to -3dB after screening, k i' Let i' be the frequency domain sampling point number corresponding to the i'th peak point;

[0023] Step 4: Determine the approximate range of the target DOA by the number of peak points after power spectrum screening;

[0024] Specifically:

[0025] When there is only one peak point, the frequency f corresponding to the peak point is used as the center, and the frequency range corresponding to one beamwidth is extended to the left and right respectively. The expanded frequency range is converted into an angle range, which serves as the approximate range where the target DOA is located.

[0026] When there are multiple peak points, the frequency f0 corresponding to the leftmost peak point and the frequency f corresponding to the rightmost peak point are used. N-1 Centered on, respectively to the left of frequency f0 and at frequency f N-1 Each of the two sides extends by one beamwidth corresponding to a frequency range, i.e. The expanded frequency range is converted into an angle range, which serves as the approximate range where the target DOA is located.

[0027] Step 5: Divide the approximate interval of the target DOA into several parts. Based on the principle of minimizing the average distance from the unit arc, calculate and update the initial values ​​for each interval.

[0028] First, the general rule for dividing the DOA intervals is: the first set of intervals is 2 times the DOA estimation interval. 1Divide into equal parts, the second group of intervals is 2 times the DOA estimation interval. 2 Divide into equal parts, the third group of intervals is 2 times the DOA estimation interval. 3 Divide into equal parts, and so on; a total of 2 parts are divided. i The number of groups is i = 1, 2, ..., Q, where Q is the approximate number of groups into which the target DOA is divided.

[0029] Then, for the current interval U ij Select a point q from the sector region containing this interval. ij The point q is such that its expected distance from all points on the arc of the sector is minimized. ij That is, the interval U ij The initial value point of the iteration;

[0030] Similarly, the initial points for iteration in all intervals are obtained;

[0031] Finally, the initial values ​​of all intervals are updated in ascending order of their indices, that is, according to the group order of the divided intervals, and within each group, they are updated in ascending order of the DOA interval range.

[0032] Step 6: According to the Root-MUSIC algorithm, perform forward and backward smoothing decoherence on the l snapshot data, calculate the covariance matrix, and further perform eigenvalue decomposition on the covariance matrix to obtain the noise subspace and polynomial coefficients.

[0033] Specifically:

[0034] Step 601: Divide the uniform linear array into p interleaved subarrays, each subarray having m array elements;

[0035] Step 602, for single snapshot data x i (t), the data model of the p subarrays corresponding to its forward and backward smoothing is as follows:

[0036]

[0037] k' = 1, 2, ..., p

[0038] A k' and B k' They are m×M data matrices; x ik' (t) represents the single snapshot data x i (t) A forward-smoothed data model of p subarrays; For single snapshot data x i (t) A backward-smoothed data model of p subarrays;

[0039] Step 603: Using the data model of p subarrays, calculate the single snapshot data x. i The data covariance matrix (t) for:

[0040]

[0041]

[0042] k' = 1, 2, ..., p

[0043] in

[0044] Step 604, using the covariance matrix of a single snapshot, calculate the covariance matrix R of the data from snapshot l as follows:

[0045]

[0046] Step 605: Perform eigenvalue decomposition on the covariance matrix R to obtain the noise subspace and polynomial coefficients;

[0047]

[0048] Among them U S U N ,Σ S ,Σ N Let represent the signal subspace, the noise subspace, the diagonal matrix composed of D large eigenvalues, and the diagonal matrix composed of MD small eigenvalues, respectively.

[0049] Depend on The coefficients c of the 2(m-1)th order polynomial are obtained. i i = 0, 1, ..., 2m-2:

[0050]

[0051] Where e ij represent The element in the i-th row and j-th column, where i,j∈[0,m-1];

[0052] Ultimately, the root-mUSIC polynomial is...

[0053] x i Let p(x) represent the i-th power of the independent variable x in the root-finding polynomial.

[0054] Step 7: Solve the polynomial using the updated initial value points to obtain the source DOA.

[0055] Specifically:

[0056] Step 701: Select initial values ​​for iteration according to the initial value update strategy, and iterate to find the root using the Muller root-finding algorithm. When the tolerance of the root between two iterations is less than the preset tolerance, the iteration ends, and the obtained root x is calculated.temp If the iteration fails to converge after reaching the maximum number of iterations, discard the initial value and select a new initial value.

[0057] Step 702, calculate the obtained root x temp The distance from the unit arc Δd = ||x temp The algorithm calculates |-1| and compares this distance with a preset threshold thres to determine the authenticity of the root. If it is a false root, the algorithm selects the initial value for the next iteration according to the initial value update strategy and returns to step 701 to continue iterating. Otherwise, if the root is found and the number of true roots is less than the number of sources, the root factor obtained this time is eliminated in the polynomial, the initial value is updated, and the algorithm continues to iterate to find the root. The algorithm ends when all true roots are found.

[0058] The process of comparing the truth values ​​of roots is as follows:

[0059] If Δd > thres, then it is determined to be a false root; otherwise, it is determined to be a true root. Continue to determine whether the number of true roots is equal to the number of sources D. If it is less than D, update the polynomial; otherwise, solve for the target DOA.

[0060] Specifically:

[0061] The number of proper roots is num, and each root is represented as r. i If i = 1, 2, ..., num, (num < D), then the polynomial is updated as follows:

[0062]

[0063] All proper roots are represented as r. i ,i=1,2,…,D, calculate the angle of the target DOA:

[0064]

[0065] Where λ is the signal wavelength, arg{r i} represents the proper root r i The phase angle.

[0066] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0067] (1) The present invention provides a method for estimating the DOA of unmanned aerial vehicles based on low-complexity Root-MUSIC. It uses FFT operation to determine the approximate range of DOA, narrows the selection range of the initial value of the root-finding iteration, and accelerates the convergence speed of the algorithm.

[0068] (2) The present invention provides a method for DOA estimation of unmanned aerial vehicles based on low-complexity Root-MUSIC. It utilizes the Muller root-finding algorithm based on factor cancellation method to avoid the initial values ​​of iterations converging to the same root, reduce the number of initial value updates, and effectively improve the DOA solution speed.

[0069] (3) The present invention provides a low-complexity Root-MUSIC-based method for estimating the DOA of unmanned aerial vehicles (UAVs). This method uses a preset threshold to determine the authenticity of the roots and terminates the algorithm when the number of true roots equals the number of DOAs of the source. This avoids the need to find the roots corresponding to all initial values ​​of the iterations and then determine the distance between the roots and the unit circle to select the true roots, thus reducing the complexity of the algorithm. Attached Figure Description

[0070] Figure 1 This is a schematic diagram of a drone DOA estimation method based on low-complexity Root-MUSIC according to the present invention.

[0071] Figure 2 This is a flowchart of a low-complexity Root-MUSIC-based method for estimating the DOA of unmanned aerial vehicles (UAVs) according to the present invention.

[0072] Figure 3 This is a schematic diagram illustrating the calculation of a single snapshot data vector using an array signal model in an embodiment of the present invention. Detailed Implementation

[0073] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments of the present invention and the features thereof can be combined with each other.

[0074] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0075] This invention discloses a low-complexity root-mUSIC-based method for estimating the target DOA of unmanned aerial vehicles (UAVs). It determines the optimal initial value and initial value update strategy for the root-finding iteration based on the FFT and LADTUC principles, and uses the Muller root-finding algorithm and a true / false root threshold to solve for the target DOA. Figure 1 As shown, it includes the following steps:

[0076] S1. Approximate estimate of the interval where the target DOA is located;

[0077] The target DOA range is determined by FFT operation: First, FFT is performed on the received multi-shot sample data to obtain the average power spectrum; then, the number of peak points of the power spectrum and the corresponding angle information are selected to determine the approximate range of the target DOA.

[0078] S2. Selection of the optimal initial value for the root-finding iteration and determination of the initial value update strategy;

[0079] Specifically, the determined approximate DOA interval is divided into several parts, and the initial iteration value in each interval is obtained according to the principle of minimizing the average distance from the unit arc: a point is selected within the sector interval such that the expected distance to all points on the sector arc is minimized; the initial iteration value is updated according to the group order of the divided intervals (i.e., in the order of the first group, the second group, and the third group), and the initial value in each group interval is selected in order of increasing DOA interval range.

[0080] S3. Calculation of coefficients of the root-finding polynomial;

[0081] Specifically, following the traditional Root-MUSIC algorithm, the coherence is first smoothed in both directions and the sample covariance matrix is ​​obtained. Then, the covariance matrix is ​​decomposed into eigenvalues ​​to obtain the noise subspace. Finally, the polynomial coefficients are obtained using the noise subspace.

[0082] S4. Solve the polynomial to obtain the source DOA;

[0083] First, an initial value is selected for iteration based on the initial value update strategy. The Muller root-finding algorithm is then used to iterate and find roots. The iteration ends when the tolerance between two iterations is less than a preset tolerance. Second, the distance between the obtained root and the unit circle is compared with a preset threshold to determine the authenticity of the root. If it is a false root, the initial value for the next iteration is selected according to the initial value update strategy. If it is a true root and the number of true roots is less than the number of sources, the root factor obtained in this iteration is eliminated from the polynomial, the initial value is updated, and the iteration continues. This process continues until all true roots are found, at which point the algorithm ends.

[0084] The UAV DOA estimation method based on low-complexity Root-MUSIC described in this invention, such as... Figure 2 As shown, it includes the following steps:

[0085] Step 1: Construct a Uniform Linear Array (ULA) signal model. A far-field narrowband signal hits the array, and the data vector of a single snapshot received by the array is calculated.

[0086] like Figure 3 As shown, consider an M-element uniform linear array with an element spacing of d = λ / 2. D far-field narrowband signals with wavelength λ strike the array, and the angle between the incident direction and the array normal direction is defined as the incident angle θ.i (i = 1, 2, ..., D), then the data vector x of the i-th snapshot received by the array i (t)

[0087] x i (t)=A(θ)s i (t)+n i (t)

[0088] s i (t) represents a D×1 dimensional signal vector, where D is the number of far-field narrowband signals;

[0089] n i (t) is M×1 dimensional white noise with zero mean and variance σ. 2 The array has a complex Gaussian distribution, and the output noise of each element is statistically independent; M is the number of elements in the uniform linear array.

[0090] A(θ)=[a(θ1),…,a(θ j ),…,a(θ D [)] is an M×D dimensional array steering matrix, composed of linearly independent array steering vectors a(θ) i The array consists of the following components, and the array steering vector is defined as follows:

[0091]

[0092] θ j It is the angle between the incident direction of the far-field narrowband signal and the normal direction of the array.

[0093] Step 2: Using the l-th snapshot data vector, perform a 2M-point Discrete Fourier Transform (DFT) to calculate the average power spectrum P of the l-th snapshot data. m (k), and normalize it;

[0094] First, for the data vector x of the i-th (i = 0, 1, ..., l-1) snapshot... i (t) Performing a 2M-point DFT (Discrete Fourier Transform) yields the frequency domain sampled signal X. i (k), X i (k) can be quickly implemented using FFT; where X i (k) is defined as:

[0095]

[0096] This represents the received data of the nth element in the i-th snapshot;

[0097] Then, based on the 2M-point DFT of the l-th snapshot data, the average power spectrum P of the l-th snapshot data can be obtained. m(k):

[0098]

[0099] After normalization, we get:

[0100]

[0101] Step 3: Count all peak points of the normalized power spectrum, filter and remove peak points with normalized power values ​​less than -3dB, and calculate the frequency of each peak point after filtering.

[0102] The frequency of the i'th peak point after filtering is

[0103] N is the number of peak points whose normalized power value is greater than or equal to -3dB after screening, k i ' represents the frequency domain sampling point number corresponding to the i-th peak point;

[0104] Step 4: Determine the approximate range of the target DOA by the number of peak points after power spectrum screening;

[0105] Specifically:

[0106] When there is only one peak point, the frequency f corresponding to the peak point is used as the center, and the frequency range corresponding to one beamwidth is extended to the left and right respectively. The expanded frequency range is converted into an angle range, which serves as the approximate range where the target DOA is located.

[0107] The formula for converting to an angle interval is: θ = sin -1 (fλ / d), where d is the element spacing and λ is the signal wavelength.

[0108] When there are multiple peak points, the frequency f0 corresponding to the leftmost peak point and the frequency f corresponding to the rightmost peak point are used. N-1 Centered on, respectively to the left of frequency f0 and at frequency f N-1 Each of the two sides extends by one beamwidth corresponding to a frequency range, i.e. The expanded frequency range is converted into an angle range, which serves as the approximate range where the target DOA is located.

[0109] Step 5: Divide the approximate interval of the target DOA into several parts. Based on the principle of minimizing the average distance from the unit arc, calculate and update the initial values ​​for each interval.

[0110] First, the general rule for dividing the DOA intervals is: the first set of intervals is 2 times the DOA estimation interval. 1 Divide into equal parts, the second group of intervals is 2 times the DOA estimation interval. 2 Divide into equal parts, the third group of intervals is 2 times the DOA estimation interval.3 Divide into equal parts, and so on; a total of 2 parts are divided. i The number of groups is i = 1, 2, ..., Q, where Q is the approximate number of groups into which the target DOA is divided.

[0111] The approximate range of DOA is Z. f Divide into Q groups, with each group having an average of 2... i If each portion is a range U, then each range U is a range U. ij The scope is:

[0112] U ij =[2πf est_min +jΔw i ,2πf est_min +(j+1)Δw i )

[0113]

[0114] i∈[1,Q],j∈[0,2 i -1]

[0115] f est_min The minimum frequency corresponding to the approximate interval of the target DOA divided by multiple peak points; f est_max The maximum frequency corresponding to the approximate interval of the target DOA divided by multiple peak points;

[0116] Then, the initial value point for each interval is selected as follows:

[0117] For the current interval U ij Select a point q from the sector region containing this interval. ij The point q is such that its expected distance from all points on the arc of the sector is minimized. ij That is, the interval U ij The initial value point of the iteration;

[0118] The specific solution method is as follows:

[0119] Let R(cosθ,sinθ) be the interval U ij Let θ be any point on the arc of the sector region, where θ is the angle between the line connecting the point to the center of the circle and the x-axis, and θ∈[θ]. min ,θ max ), θ min ,θ max They are intervals U ij The left and right boundaries are defined by θ; θ follows a uniform distribution with the following probability density function:

[0120] Select point q ij The expected value of the squared distance to R(cosθ,sinθ) is:

[0121]

[0122] For L 2 (x ij ,y ij Find the partial derivative and set it to zero:

[0123]

[0124] Easy to obtain:

[0125] Similarly, the initial points for iteration in all intervals are obtained;

[0126] Finally, the initial values ​​of all intervals are updated in ascending order of their indices, i.e., according to the group order of the divided intervals. Update: Within each interval, update in ascending order of DOA interval range.

[0127] Step 6: According to the Root-MUSIC algorithm, perform forward and backward smoothing decoherence on the l snapshot data, calculate the covariance matrix, and further perform eigenvalue decomposition on the covariance matrix to obtain the noise subspace and polynomial coefficients.

[0128] Specifically:

[0129] Step 601: Divide the uniform linear array into p interleaved subarrays, each subarray having m elements; define the following two m×M data matrices:

[0130]

[0131] Among them I m×m J m×m Let k' be the identity matrix and the permutation matrix with anti-diagonal line equal to 1, respectively.

[0132] Their data covariance matrix can be obtained as follows:

[0133] k' = 1, 2, ..., p

[0134] in

[0135] Step 602, for single snapshot data x i (t), the data model of the p subarrays corresponding to its forward and backward smoothing is as follows:

[0136]

[0137] k' = 1, 2, ..., p

[0138] x ik'(t) represents the single snapshot data x i (t) A forward-smoothed data model of p subarrays; For single snapshot data x i (t) A backward-smoothed data model of p subarrays;

[0139] Step 603: Using the data model of p subarrays, calculate the single snapshot data x. i The data covariance matrix (t) for:

[0140]

[0141] Step 604, using the covariance matrix of a single snapshot, calculate the covariance matrix R of the data from snapshot l as follows:

[0142]

[0143] Step 605: Perform eigenvalue decomposition on the covariance matrix R to obtain the noise subspace and polynomial coefficients;

[0144]

[0145] Among them U S U N ,Σ S ,Σ N Let represent the signal subspace, the noise subspace, the diagonal matrix composed of D large eigenvalues, and the diagonal matrix composed of MD small eigenvalues, respectively.

[0146] Depend on The coefficients c of the 2(m-1)th order polynomial are obtained. i i = 0, 1, ..., 2m-2:

[0147]

[0148] Where e ij represent The element in the i-th row and j-th column, where i,j∈[0,m-1];

[0149] Ultimately, the root-mUSIC polynomial is...

[0150] x i Let p(x) represent the i-th power of the independent variable x in the root-finding polynomial.

[0151] Step 7: Solve the polynomial using the updated initial value points to obtain the source DOA.

[0152] Specifically:

[0153] Step 701: Select initial values ​​for iteration according to the initial value update strategy, and iterate to find the root using the Muller root-finding algorithm. When the tolerance of the root between two iterations is less than the preset tolerance, the iteration ends, and the obtained root x is calculated. temp If the iteration fails to converge after reaching the maximum number of iterations, discard the initial value and select a new initial value.

[0154] The specific process of Muller's iterative root finding is as follows:

[0155] The three initial values ​​for the roots of the Muller root-finding algorithm are obtained from the initial value q during iteration:

[0156] x k-1 =q,

[0157] Where w q Let Δw be the phase angle of q, and Δw be the interval U in which q lies. ij Divide the x into 10 equal parts, and the quadratic polynomial function value corresponding to the three points is p(x k-2 ), p(x k-1 ), p(x k );

[0158] By (x) k-2 ,p(x k-2 )),(x k-1 ,p(x k-1 )),(x k ,p(x k Construct the parabola y(x) = ax 2 +bx+c, introducing a new variable simplifies the calculation:

[0159] Substituting the variables into the parabola equation, we get:

[0160] in:

[0161] The two zeros of y(λ) are:

[0162] Take the zero point of y(λ) that is closer to x k At that zero point, we have:

[0163] The iterative formula for the Muller method is:

[0164] x k+1 =x k +λ4(x k -x k-1 )

[0165] Continuously update the value of the root until |x k+1 -x k If | < γ (γ is the tolerance) or the number of iterations equals the preset maximum number of iterations, the iteration ends.

[0166] Step 702, calculate the obtained root x temp The distance from the unit arc Δd = ||x temp The algorithm calculates |-1| and compares this distance with a preset threshold thres to determine the authenticity of the root. If it is a false root, the algorithm selects the initial value for the next iteration according to the initial value update strategy and returns to step 701 to continue iterating. Otherwise, if the root is found and the number of true roots is less than the number of sources, the root factor obtained this time is eliminated in the polynomial, the initial value is updated, and the algorithm continues to iterate to find the root. The algorithm ends when all true roots are found.

[0167] The process of comparing the truth values ​​of roots is as follows:

[0168] If Δd > thres, then it is determined to be a false root; otherwise, it is determined to be a true root. Continue to determine whether the number of true roots is equal to the number of sources D. If it is less than D, update the polynomial; otherwise, solve for the target DOA.

[0169] Specifically:

[0170] The number of proper roots is num, and each root is represented as r. i If i = 1, 2, ..., num, (num < D), then the polynomial is updated as follows:

[0171]

[0172] All proper roots are represented as r. i ,i=1,2,…,D, calculate the angle of the target DOA:

[0173]

[0174] Where λ is the signal wavelength, arg{r i} represents the proper root r i The phase angle.

[0175] After the above steps, the target DOA can be estimated using the low-complexity Root-MUSIC algorithm.

[0176] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A low-complexity Root-MUSIC based UAV DOA estimation method, characterized in that, The specific steps are as follows: First, a uniform linear array ULA signal model is constructed. A far-field narrowband signal is applied to this array, and the vector of a single snapshot data received by the array is calculated. Then, using... The single snapshot data vectors are processed separately. After point-based Discrete Fourier Transform (DFT), calculate the average power spectrum. and normalization; In the single snapshot data vector, the first... Secondary snapshot data vector The calculation is as follows: express 3D signal vector The number of far-field narrowband signals; for 3D white noise, The number of elements in a uniform linear array; yes 3D array steering matrix; The angle between the incident direction of the far-field narrowband signal and the normal direction of the array; Then, all peak points of the normalized power spectrum are counted, and peak points with normalized power values ​​less than -3dB are filtered out and removed. The frequency of each peak point after filtering is calculated. The approximate range of the target DOA is determined by the number of peak points after filtering the power spectrum. Specifically: When there is only one peak point, use the frequency corresponding to the peak point. Centered on the center, extend to the left and right by a frequency range corresponding to one beamwidth, that is... The expanded frequency range is converted into an angle range, which serves as the approximate range where the target DOA is located. When there are multiple peak points, the frequency corresponding to the leftmost peak point is used. and the frequency corresponding to the rightmost peak point Centered on, at different frequencies Left side and frequency Each of the two sides extends by one beamwidth corresponding to a frequency range, i.e. The expanded frequency range is converted into an angle range, which serves as the approximate range where the target DOA is located. Next, the approximate interval of the target DOA is divided into several parts. Based on the principle of minimizing the average distance from the unit arc, the initial values ​​for iteration within each interval are obtained and updated. The general rule for dividing the DOA intervals is as follows: the first set of intervals is the DOA estimation interval. Divide into equal parts, the second group of intervals is the DOA estimation interval. Divide into equal parts, the third group of intervals is the DOA estimation interval. Divide into equal parts, and so on; the total area is divided into 100 parts. share, , The number of groups to divide the approximate range of the target DOA; For the current interval Select a point from the sector area containing that interval. The point is such that the expected distance to all points on the arc of the sector is minimized. That is, the interval The initial value point of the iteration; Similarly, we obtain the initial points for iteration in all intervals; Update the initial values ​​of all intervals in ascending order of their indices, that is, update them according to the group order of the divided intervals, and update them in ascending order of the DOA interval range within each group. Finally, according to the Root-MUSIC algorithm, The data from the next snapshot is smoothed and decohered in both directions, and the covariance matrix is ​​calculated. The covariance matrix is ​​then subjected to eigenvalue decomposition to obtain the noise subspace and polynomial coefficients. The polynomial is solved using the updated initial values ​​of the iteration to obtain the source DOA.

2. The UAV DOA estimation method based on low-complexity Root-MUSIC as described in claim 1, characterized in that, The Average power spectrum of the second snapshot data The calculation formula is: For the first Secondary snapshot data vector conduct Point DFT yields the frequency domain sampled signal; Indicates the first The second quick shot The received data of each array element.

3. The method for estimating the DOA of a UAV based on low-complexity Root-MUSIC as described in claim 1, characterized in that, The filtered first The frequency of each peak point is ; This represents the number of peak points whose normalized power values ​​are greater than or equal to -3 dB after screening. For the first The frequency domain sampling point sequence number corresponding to each peak point.

4. The UAV DOA estimation method based on low-complexity Root-MUSIC as described in claim 1, characterized in that, The pair The snapshot data is subjected to forward and backward smoothing decoherence, the covariance matrix is ​​calculated, and eigenvalue decomposition is used to obtain the noise subspace and polynomial coefficients; specifically: Step 601: Divide the uniform linear array into interleaved... There are several subarrays, and the number of array elements in each subarray is... ; Step 602, for single snapshot data Its forward smoothing and backward smoothing correspond to The data model for each subarray is: and They are respectively The data matrix; For single snapshot data Forward smooth Data model for each subarray; For single snapshot data Backward smoothing Data model for each subarray; Step 603, using The data model of each subarray is used to calculate single snapshot data. Data covariance matrix for: in ; Step 604: Calculate the covariance matrix using the data from a single snapshot. Covariance matrix of the second snapshot data for: Step 605, for the covariance matrix Perform eigenvalue decomposition to obtain the noise subspace and polynomial coefficients; in Representing the signal subspace, noise subspace, respectively A diagonal matrix composed of large eigenvalues A diagonal matrix composed of small eigenvalues; Depend on get coefficients of the order polynomial : in represent The Middle Line 1 Column elements, ; Ultimately, the root-mUSIC polynomial is... ; Represents the root-finding polynomial independent variable of Power of 1.

5. The method for estimating the DOA of a UAV based on low-complexity Root-MUSIC as described in claim 1, characterized in that, The polynomial coefficients are solved using the updated iterative initial value points to obtain the source DOA; Specifically: Step 701: Select initial values ​​for iteration according to the initial value update strategy, and iteratively find roots based on the Muller root-finding algorithm. When the tolerance of the root between two iterations is less than the preset tolerance, end the iteration and record the obtained roots. Remain in step 702; Otherwise, continue iterating. If convergence is not achieved after reaching the maximum number of iterations, discard the initial value for this iteration and select a new initial value. Step 702, calculate the obtained root Distance from the unit arc and compare this distance with a preset threshold. The algorithm compares the roots to determine their authenticity: if a root is false, the algorithm selects the initial value for the next iteration based on the initial value update strategy and returns to step 701 to continue iterating; otherwise, if the root is true and the number of true roots is less than the number of sources, the algorithm eliminates the root factor obtained this time in the polynomial, updates the initial value, and continues iterating to find roots; the algorithm ends when all true roots are found. The process of comparing the truth values ​​of roots is as follows: if If it is true, then it is determined to be a false root; otherwise, it is determined to be a true root. Continue to determine whether the number of true roots equals the number of sources. If less than If the polynomial is updated, then the target DOA is solved; otherwise, the target DOA is solved. Specifically: The number of true roots is Each root is represented as The polynomial is then updated to: All proper roots are represented as Calculate the angle of the target DOA: in, For the signal wavelength, Indicates true root The phase angle.

Images