A method for direction of arrival estimation based on golden section sparse circular array

Abstract

This invention proposes a direction-of-arrival (DOA) estimation method based on a sparse circular array using the golden ratio. The method defines the antenna position and three-dimensional coordinates using the golden ratio, constructing a sparse array. A DOA vector is established based on the azimuth and elevation angles, and the path difference and phase difference are calculated. These are then combined with multi-shot data to form a covariance matrix. To improve robustness under low snapshot and low signal-to-noise ratio conditions, adaptive diagonal loading is employed to adjust the noise eigenvalue distribution and enhance the distinguishability between the signal and noise subspaces. Addressing the issues of traditional MUSIC algorithms' accuracy depending on step size and poor real-time performance, a three-step search is designed: a coarse search to locate local maxima, a ternary trigonometric iteration to narrow the interval, and a fine search followed by quadratic surface fitting to eliminate grid errors. This method effectively suppresses sparse array grating lobes, balancing high accuracy and real-time performance.

Description

Technical Field

[0001] This invention belongs to the field of positioning, and in particular relates to the direction of arrival estimation technique for sparse antenna arrays. Background Technology

[0002] In the field of Direction of Arrival (DOA) estimation, the MUSIC (Multiple Signal Classification) algorithm has been widely used due to its super-resolution DOA estimation capabilities. Among these, sparse antenna arrays have attracted significant attention because they reduce the number of physical array elements while lowering hardware costs and system complexity. However, sparse antenna arrays face several key challenges in practical applications. Specifically, the sparse antenna distribution generates a large number of grating lobes, leading to a significant decrease in the DOA estimation accuracy of the MUSIC algorithm. Furthermore, the DOA estimation accuracy of traditional MUSIC algorithms is highly dependent on the search step size; a smaller step size results in higher DOA estimation accuracy, but simultaneously, the computational load increases exponentially, making single-snapshot DOA estimation too time-consuming and unacceptable for real-time requirements. Even more critically, in scenarios with low snapshot numbers or low signal-to-noise ratios, the estimation bias of the covariance matrix is ​​large, leading to a significant increase in DOA estimation error. Currently, although some improved solutions attempt to enhance algorithm performance by optimizing the search step size or adjusting the array layout, these solutions fail to simultaneously resolve the contradictions between real-time performance, accuracy, and robustness, and lack targeted optimization measures for complex scenarios such as low signal-to-noise ratio and low snapshot count. Therefore, developing a direction-of-arrival (DOA) estimation scheme that balances array layout rationality, algorithm real-time performance, DOA estimation accuracy, and robustness in complex scenarios has become a critical problem urgently needing to be solved in the current DOA estimation field. The DOA estimation method based on a golden section sparse circular array proposed in this invention is designed precisely to address these issues. Summary of the Invention

[0003] Based on this, this invention addresses the problems of poor accuracy and real-time performance of sparse antenna arrays by proposing a direction-of-arrival (DOA) estimation method based on a golden ratio sparse circular array. This invention achieves a synergistic improvement in "high accuracy, high real-time performance, and high robustness" by optimizing the array layout, improving the covariance matrix processing method, and refining the MUSIC search logic.

[0004] This invention first provides a model based on the golden ratio sparse antenna array:

[0005] The array optimizes the spatial distribution of array elements using the golden ratio. The antenna positions and angles are determined according to the golden ratio. Let the radius of the circular array be R, and there be a total of P antennas. A coordinate system is constructed with the center of the circular array as the origin and the direction of the first antenna as the positive x-axis. The position and angle of the p-th antenna... Based on the golden ratio design, the position coordinate vector is .

[0006] Based on the above model, the specific steps of the direction-of-arrival estimation method based on the golden section sparse circular array provided by this invention are as follows:

[0007] Step S1: Construct a golden ratio sparse circular array. The position angle and three-dimensional coordinates of each antenna are determined strictly according to the golden ratio to ensure optimal spatial distribution of array elements, reducing the interference of grating lobes on estimation accuracy from the source, while maintaining the low-cost advantage of sparse arrays.

[0008] Step S2: Construct the received signal model and calculate its covariance matrix. Based on the golden ratio sparse circular matrix, a received signal model is constructed, representing the received signal as a superposition of the incident signal and additive white Gaussian noise. The covariance matrix of the received signal is calculated based on this model. To address the problem of large estimation deviations in the covariance matrix under low snapshot number and low signal-to-noise ratio scenarios, an adaptive diagonal loading strategy is introduced for optimization. The loading factor adopts a dynamic adaptive design, which can dynamically adjust the loading intensity according to the signal environment to ensure the stability and accuracy of the covariance matrix. Subsequently, eigenvalue decomposition is performed on the optimized covariance matrix to separate the signal subspace and noise subspace, and the MUSIC spectral function is calculated, providing a reliable spectral function basis for subsequent angle search.

[0009] Step S3: Improve the MUSIC search method. To break the contradiction between "accuracy and real-time performance" in traditional search step sizes, this invention designs a three-level progressive search mechanism of coarse search, fine search, and quadratic surface fitting. The first stage, coarse search, involves sampling within a two-dimensional range of azimuth and pitch angles. Potential peak regions are identified through local maxima detection and gradient verification. If gradient verification fails, a random selection mechanism using surrounding grid points prevents the omission of optimal solutions and rapidly narrows the search range. The second stage, fine search, uses the points found in the coarse search as centers and generates new search intervals through uniform random perturbation. The interval is iteratively narrowed using a ternary search method. After the iteration conditions are met, the interval is expanded slightly by random perturbation again. A small-step traversal using 31 precise grid points accurately locates the optimal region near the peak of the spectral function. The third stage, quadratic surface fitting, utilizes the quadratic approximation properties of the surface near the peak of the MUSIC spectrum. A quadratic surface is constructed using 9 points (8 surrounding grid points) centered on the points found in the fine search. The model coefficients are solved using the least squares method, and the extreme points are solved using Cramer's rule. Finally, a high-precision angle estimate at the sub-grid level is obtained, breaking the limitations of traditional grid search step sizes.

[0010] The beneficial effects of this invention are as follows: This invention forms a complete "hardware-algorithm-logic" optimization scheme through hardware optimization of the golden section array layout, robust covariance matrix processing with adaptive diagonal loading, and improvement of the three-stage progressive search logic. Specifically, the combination of the golden section array and diagonal loading technology effectively suppresses grating lobe interference from sparse arrays, significantly reducing angle estimation errors. The three-stage search mechanism, based on rapid range locking through coarse search, achieves a breakthrough in accuracy through fine search and quadratic surface fitting. Compared to the traditional MUSIC algorithm, it significantly improves computational efficiency while maintaining improved estimation accuracy, and maintains stable performance even in complex scenarios with low snapshot numbers and low signal-to-noise ratios. This method is applicable to multiple fields such as wireless communication, radar detection, sonar positioning, and electronic reconnaissance, meeting the accuracy and real-time requirements of DOA estimation in different scenarios, and has broad engineering application prospects and promotional value. Attached Figure Description

[0011] Figure 1 This is a model diagram of the sparse antenna array of the present invention; Figure 2 This is a diagram showing the phase difference analysis of the far-field parallel wave of the present invention. Detailed Implementation

[0012] like Figure 1 As shown, the direction-of-arrival (DOA) estimation method based on the golden section sparse circular array, specifically a method for estimating DOA based on the golden section sparse circular array, has the following detailed steps:

[0013] Step S1: Construct a golden ratio sparse circular matrix; specifically, the first... The position and angle of each antenna are determined based on the golden ratio to optimize spatial sampling.

[0014]

[0015] With the center of the circle as the origin, and the direction of the first antenna as the positive x-axis; where, Indicates the first The angle by which the straight line between the antenna and the center of the circle rotates counterclockwise relative to the positive x-axis. Indicates to Modulo operation; Indicates the total number of antennas; the first The position coordinate vector of each antenna Calculated using the golden ratio:

[0016]

[0017] in, For the first The three-dimensional coordinates of each antenna Let be the radius of the circular array.

[0018] Step S2: Construct the received signal model and calculate the covariance matrix, then perform adaptive diagonal loading processing on it. For the azimuth angle... With pitch angle directional vector for:

[0019]

[0020] like Figure 2 As shown, the path difference between different antennas is calculated using the first antenna of the array as a reference; the incident signal at infinity can be assumed to be a parallel wave, then the... Path difference between the signal received by the first antenna and the signal received by the first reference antenna for:

[0021]

[0022] in, Let be the position coordinate vector of the first antenna; after obtaining the path difference, the position coordinate vector of the second antenna can be calculated. Phase difference between the signals received by the antenna and the first reference antenna :

[0023]

[0024] The first phase difference can be used to construct the second phase difference. The steering vector of the received signal :

[0025]

[0026] in, This represents the phase difference between the signals received by the second antenna and the first reference antenna. For the 1st to the 1st The complex form of the phase difference of the received signals is given; assuming there are a total of K received signals, a received signal model is constructed. :

[0027]

[0028]

[0029]

[0030] in, This represents the steering vector matrix for K received signals; These represent the signal power of the first to the Kth signal sources, respectively. The incident signal matrix; These are additive white Gaussian noises with a mean of 0, representing the receivers of antennas 1 through P. It is an additive white Gaussian noise matrix with a mean of 0; covariance matrix The following formula is used to calculate:

[0031] The covariance matrix is ​​subjected to adaptive diagonal loading, and the adaptive diagonal loading method is as follows:

[0032]

[0033] in, As the loading factor, The identity matrix has the same dimensions as the covariance matrix; loading factor. The definition is as follows:

[0034]

[0035] in, This is a scaling factor constant term, determined a priori. Under the same conditions, different scaling factors are verified, and the scaling factor with the smallest angle measurement error is selected. Represents the trace of the covariance matrix; The number of antennas is the dimension of the covariance matrix. The signal-to-noise ratio represents the linear scale, i.e., the linear form derived from the decibel form. For the obtained covariance matrix... Eigenvalue decomposition is performed on it:

[0036]

[0037] in, , These are the signal subspace and the noise subspace, respectively. The number of antennas, The number of signal sources; and Let be the diagonal matrices representing the signal eigenvalues ​​and the noise eigenvalues, respectively; calculate the MUSIC spectral function according to the definition of the MUSIC spectral function. :

[0038]

[0039] Step S3: Improve the search method of the MUSIC algorithm. First, perform a coarse search within the two-dimensional search range. The coarse search stage is essentially a uniform sampling in the parameter space, defining the search range for azimuth and elevation angles:

[0040]

[0041] in, This indicates the left boundary of the azimuth search range. Indicates the right boundary of the azimuth search range; This indicates the lower boundary of the pitch angle search range. This represents the upper boundary of the pitch angle search range; here we take... Generate discrete grid points within the search range. and{ Perform a coarse search using the following search method:

[0042]

[0043]

[0044] in, and Let be the search step size for the grid. To ensure the accuracy of the search range boundaries, the step size for the coarse search is . Set all angles to 5. Set to 1 within the angle range; and The number of grid points; compare the MUSIC spectral values ​​of all grid points. Then, assuming the searched grid points are If the grid points near a point satisfy the following conditions, then the point... Local maxima:

[0045]

[0046]

[0047] in, and They are respectively The previous and next azimuth grid points, and They are respectively The previous and next pitch angle grid points; then for Perform gradient validation:

[0048]

[0049]

[0050] in, and Indicates the spectral function respectively Calculate the gradients in the azimuth and pitch directions. If the gradient at a point satisfies the true peak characteristics of the two-dimensional MUSIC spectral function (i.e., the grid points to the left of the peak have a gradient greater than zero, the grid points to the right of the peak have a gradient less than zero, the grid points below the peak have a gradient greater than zero, and the grid points above the peak have a gradient less than zero), then that point can be used as a coarse search point. If the gradient verification at that point does not meet the conditions, or the gradient is less than a preset threshold... Using this point as the center, eight grid points are randomly selected around it; the randomly selected point is directly used as the point for coarse search.

[0051] Then, using the grid points obtained from the coarse search as the center... By employing random perturbation, a new search interval is formed in its vicinity:

[0052]

[0053]

[0054]

[0055] in, Indicates uniform distribution. , To meet the interval Uniform distribution and These are the left and right boundaries of the new azimuth search interval, respectively. and These represent the left and right boundaries of the new pitch angle search interval, respectively; the probability of a perturbation to the left and right of the coarse search point is completely equal, ensuring an unbiased search of the parameter space; in the azimuth angle... With pitch angle New search zones were formed in the vicinity. and .

[0056] In the fine search phase, the search interval is iteratively narrowed using the ternary method; taking the azimuth angle as an example, the left dividing point is calculated. right dividing point :

[0057]

[0058] Compare the MUSIC spectral values ​​at two points and ;like > Then the new interval for the next iteration is ;like Then the new interval for the next iteration is ;like = Then the new interval for the next iteration is Until the maximum number of iterations (30) is reached or the change occurs between two adjacent iteration intervals. Less than the convergence threshold When the interval change occurs, the iteration stops. Defined as

[0059]

[0060] in, Indicates the first iteration and For the first The left and right intervals of the next iteration; when the iteration converges, the final convergent interval is defined as... At this point, the convergence interval is expanded again using random perturbation:

[0061]

[0062]

[0063] in, and To comply with The uniform distribution; the interval obtained after perturbation This serves as the next search interval; similarly, we can obtain... The search range is the pitch angle range; within the search range and A fine search is performed using a small grid step size much smaller than that of the coarse search. Regardless of the size of the search interval at this time, the search intervals in both azimuth and elevation angles are divided into 31 more precise grid points. Then, by comparing the MUSIC spectral values ​​of all grid points, the grid point with the maximum MUSIC spectral value is taken as the center point for the next stage.

[0064] In two-dimensional music spectrum In the model, the surface near the peak can be approximated as a quadratic surface, which breaks the limitation of the traditional search step size and further improves the positioning accuracy. Therefore, quadratic surface fitting is performed on the two-dimensional MUSIC spectral function:

[0065]

[0066] in, The coefficients to be determined are: This represents the extreme point to be found, i.e., the final angle estimate; assuming the angle found in the local fine search stage is... , express The MUSIC spectral values; Centered on the grid, select 8 surrounding grid points for a total of 9 points. Constructing a quadratic surface, we can list the equations at each point:

[0067]

[0068] in The MUSIC spectral values ​​are for the selected 8 grid points; the equation is then transformed into matrix form:

[0069]

[0070] in,

[0071]

[0072] Solve the coefficient matrix using the least squares method. Least squares objective function for:

[0073]

[0074] right Taking the derivative and setting it to zero, we obtain the normal equation:

[0075]

[0076] Simplifying, we get:

[0077]

[0078] After obtaining the coefficient matrix, based on the condition that the first derivative of a quadratic surface is zero, we can obtain the system of equations:

[0079]

[0080] in, This represents the extreme point to be found, i.e., the final estimated value; and Two-dimensional MUSIC spectra right and Find the partial derivative; applying Cramer's rule to this equation directly yields the expression for the extreme points, from which the final estimate can be directly calculated. :

[0081]

[0082] Thus, direction-of-arrival estimation based on a sparse circular array with the golden ratio has been achieved.

[0083] The above-described embodiments are merely preferred embodiments of the present invention and are not intended to limit the present invention. The scope of protection of the present invention is not limited to the above embodiments. For those skilled in the art, all equivalent substitutions, modifications, alterations, etc., made according to the disclosure of the present invention should be included within the scope of protection of the present invention.

Claims

1. A direction-of-arrival estimation method based on a golden ratio sparse circular array, characterized in that: This invention includes a method for constructing a geometric model of a sparse circular array based on the golden ratio, an adaptive diagonal loading method for processing the covariance matrix of the received signal, and a three-step angle search method. The method for constructing the geometric model of the sparse circular array based on the golden ratio uses the center of the circle as the origin and the direction of the first antenna as the positive x-axis. Starting from the positive x-axis, the position angles of each antenna in the sparse circular array are determined in a counter-clockwise direction. The coordinate vectors of the other antenna positions are calculated using the center of the circle as a reference. The method for processing the covariance matrix of the received signal uses the first antenna along the positive x-axis of the antenna array as a reference, calculates the phase difference between different antennas through incident wave vector projection, obtains the received signal model, and calculates the covariance matrix. The method then applies the covariance matrix... The difference matrix undergoes adaptive diagonal loading, with the loading factor determined by a combination of scaling factor, covariance matrix trace, number of antennas, and signal-to-noise ratio. The three-step angle search method utilizes the MUSIC algorithm to perform eigenvalue decomposition on the processed covariance matrix, obtaining the signal subspace and noise subspace. A three-step search approach is employed: a coarse search stage uniformly samples and verifies local maxima points in the parameter space; a small-interval search iteratively narrows the search interval using a ternary search method; a local fine search performs refined sampling within a small range; and finally, quadratic surface fitting is used to solve for the surface coefficients using the least squares method to determine the final estimated angle, i.e., the angle corresponding to the highest peak of the MUSIC spectral function is selected as the estimated target direction of arrival.

2. The direction-of-arrival estimation method based on a golden ratio sparse circular array as described in claim 1, characterized in that, The method for constructing the geometric model of the golden ratio sparse circular matrix includes constructing a golden ratio sparse circular matrix; specifically, the first... The position and angle of each antenna are determined based on the golden ratio to optimize spatial sampling. ； With the center of the circle as the origin, and the direction of the first antenna as the positive x-axis; where, Indicates the first The angle by which the straight line between the antenna and the center of the circle rotates counterclockwise relative to the positive x-axis. Indicates to Modulo operation; Indicates the total number of antennas; the first The position coordinate vector of each antenna Calculated using the golden ratio: ； in, For the first The three-dimensional coordinates of each antenna Let be the radius of the circular array.

3. The direction-of-arrival estimation method based on a golden ratio sparse circular array as described in claim 1, characterized in that, The adaptive diagonal loading processing method for the received signal covariance matrix includes the following steps: Step S1, for azimuth angle With pitch angle directional vector for: ； The path difference between different antennas is calculated using the first antenna of the array as a reference; the incident signal at infinity can be assumed to be a parallel wave, then the... Path difference between the signal received by the first antenna and the signal received by the first reference antenna for: ; in, Let be the position coordinate vector of the first antenna; after obtaining the path difference, the position coordinate vector of the second antenna can be calculated. Phase difference between the signals received by the antenna and the first reference antenna : ； The first phase difference can be used to construct the second phase difference. The steering vector matrix of the received signal : ； in, This represents the phase difference between the signals received by the second antenna and the first reference antenna. For the 1st to the 1st The complex form of the phase difference of the received signals; assuming there are a total of K received signals, the received signal model matrix can be constructed. : ； ； ； in, This represents the steering vector matrix for K received signals; These represent the signal power of the first to the Kth signal sources, respectively. The incident signal matrix; These are additive white Gaussian noises with a mean of 0, representing the receivers of antennas 1 through P. It is an additive white Gaussian noise matrix with a mean of 0; the covariance matrix can be derived from it. : ; Where M is the number of snapshots. Indicates to The conjugate transpose of a matrix; Step S2: Perform adaptive diagonal loading on the covariance matrix. The adaptive diagonal loading method is as follows: ； in, The loaded covariance matrix, As the loading factor, The identity matrix has the same dimensions as the covariance matrix; loading factor. The definition is as follows: ； in, This is a scaling factor constant term, determined a priori. Under the same conditions, different scaling factors are verified, and the scaling factor with the smallest angle measurement error is selected. Represents the trace of the covariance matrix; The number of antennas is the dimension of the covariance matrix. The signal-to-noise ratio represents the linear scale, which is the linear form derived from the decibel form.

4. The direction-of-arrival estimation algorithm based on a golden section sparse circular array as described in claim 1, characterized in that, The three-step angle search method includes the following steps: Step S1: For the obtained covariance matrix Eigenvalue decomposition is performed on it: ； in, , These are the signal subspace and the noise subspace, respectively. The number of antennas, The number of signal sources; and Let be the diagonal matrices representing the signal eigenvalues ​​and the noise eigenvalues, respectively; calculate the MUSIC spectral function according to the definition of the MUSIC spectral function. : ； Step S2: Perform a coarse search within the two-dimensional search range; the coarse search stage essentially involves uniform sampling in the parameter space to define the search range for azimuth and elevation angles. ； in, This indicates the left boundary of the azimuth search range. Indicates the right boundary of the azimuth search range; This indicates the lower boundary of the pitch angle search range. This represents the upper boundary of the pitch angle search range; discrete grid points are generated within the search range. and{ Perform a coarse search using the following search method: ； ； in, and This represents the search step size for the grid. and The number of grid points; compare the MUSIC spectral values ​​of all grid points. Then, assuming the searched grid points are If the grid points near this point satisfy the following conditions: ； ； Then point This is a local maximum point; where, and They are respectively The previous and next azimuth grid points, and They are respectively The previous and next pitch angle grid points; then for Perform gradient validation: ； ； in, and Indicates the spectral function respectively Calculate the gradients in the azimuth and pitch directions. If the gradient at a point satisfies the true peak characteristics of the two-dimensional MUSIC spectral function (i.e., the grid points to the left of the peak have a gradient greater than zero, the grid points to the right of the peak have a gradient less than zero, the grid points below the peak have a gradient greater than zero, and the grid points above the peak have a gradient less than zero), then that point can be used as a coarse search point. If the gradient verification at that point does not meet the conditions, or the gradient is less than a preset threshold... Using this point as the center, eight grid points are randomly selected around it; the randomly selected points are directly used as the points for the coarse search; then the grid points obtained from the coarse search are used... Centered on a target area, a new search range is formed in its vicinity using random perturbation. ； ； ； in, Indicates uniform distribution. , To meet the interval Uniform distribution and These represent the left and right boundaries of the new azimuth search interval, respectively. and These represent the left and right boundaries of the new pitch angle search interval, respectively; the probability of a perturbation to the left and right of the coarse search point is completely equal, ensuring an unbiased search of the parameter space; in the azimuth angle... With pitch angle New search zones were formed in the vicinity. and ; Step S3, the fine search stage, involves iteratively narrowing the search interval using the ternary method; taking the azimuth angle as an example, the left dividing point is calculated. right dividing point : ； Compare the MUSIC spectral values ​​at two points and ;like > Then the new interval for the next iteration is ;like Then the new interval for the next iteration is ;like = Then the new interval for the next iteration is Until the maximum number of iterations is reached or the change between two adjacent iteration intervals is reached. Less than the convergence threshold When the interval change occurs, the iteration stops. Defined as ； in, Indicates the first iteration and For the first The left and right boundaries of the interval in the next iteration; when the iteration converges, the interval of final convergence is defined as... At this point, the convergence interval is expanded again using random perturbation: ； ； in, and To comply with Uniform distribution; and These represent the left and right boundaries of the azimuth interval after the second random perturbation; and the interval obtained after the perturbation. This serves as the next search interval; similarly, we can obtain... The pitch angle search range, and These represent the left and right intervals of the pitch angle after the second random perturbation; then within the search interval... and A fine search is performed with a grid step size much smaller than that of the coarse search. Regardless of the size of the search interval at this time, the search intervals in the azimuth and elevation angles are divided into more precise grid points. Then, by comparing the MUSIC spectral values ​​of all grid points, the grid point with the maximum MUSIC spectral value is taken as the center point of the next stage. Step S4: In the two-dimensional music spectrum In the process, a quadratic surface fitting is performed, and the surface near the peak can be approximated as a quadratic surface: ； in, The coefficients to be determined are: This represents the extreme point to be found, i.e., the final angle estimate; assuming the angle found in the local fine search stage is... , express The MUSIC spectral values; Centered on the grid, select 8 surrounding grid points for a total of 9 points. Constructing a quadratic surface, we can list the equations at each point: ； in The MUSIC spectral values ​​are for the selected 8 grid points; the equation is then transformed into matrix form: ； ； in, Represents the MUSIC spectral value matrix. Represents a quadratic term matrix. The coefficient matrix is ​​solved using the least squares method. Least squares objective function for: ； right Taking the derivative and setting it to zero, we obtain the normal equation: ； Simplifying, we get: ； After obtaining the coefficient matrix, based on the condition that the first derivative of a quadratic surface is zero, we can obtain the system of equations: ； in, This represents the extreme point to be found, i.e., the final estimated value; and Two-dimensional MUSIC spectra right and Find the partial derivative; applying Cramer's rule to this equation directly yields the expression for the extreme points, from which the final estimate can be directly calculated. : ； Ultimately, a direction-of-arrival (DOA) estimation based on a sparse circular array with the golden ratio is achieved.

Images