Compact array DOA estimation method based on alternating projection and golden section search method

Abstract

The invention discloses a compact array DOA estimation method based on an alternating projection and golden section search method, and the method comprises the steps: inputting a radar array signal, initializing the total number K of search directions of arrival, an estimation direction of arrival vector theta, a current estimation direction of arrival serial number k, and a counter of k; preliminarily estimating the kth direction of arrival by using a golden section search method; dynamically updating the estimated direction-of-arrival vector theta by using an alternating projection algorithm; judging whether the current estimated direction-of-arrival serial number k is greater than the total number K of the searched direction-of-arrival; if yes, the estimation of all directions of arrival is completed, and the search is ended; otherwise, updating the current estimated direction-of-arrival serial number k and the counter, and continuously estimating the direction-of-arrival. According to the method, the accuracy, the stability and the resolution of direction of arrival estimation are improved, the required calculation time is shortened, and the requirement of real-time processing is better met.

Description

Technical Field

[0001] The invention belongs to the field of target direction estimation, and in particular relates to a compact array DOA estimation method based on alternating projection and golden section search method. Background Art

[0002] Reducing the element spacing and miniaturizing the array are of great significance in modern technology applications. The element spacing of traditional arrays is usually set to half the wavelength of the carrier, which makes the array bulky and affects its transportation and use. The element spacing of compact arrays is less than half the wavelength of the carrier, which enables the miniaturization of the array, significantly improving the flexibility and mobility of the system, making it easy to deploy in various complex environments. In addition, the use of miniaturized design can also reduce material and manufacturing costs, further improving the overall performance and reliability of the system. These factors have jointly promoted the advancement of modern communications, radar, and sensor technologies to meet the growing market demands and application challenges.

[0003] However, the reduction in the array element spacing will cause the array main lobe to become wider, thereby reducing the directivity of the signal and affecting the accuracy of the direction of arrival (DOA) estimation. This is particularly critical for applications that require high-precision positioning, such as radar and sonar. At the same time, the wide main lobe makes it impossible to accurately distinguish the contribution of each signal source in the case of similar signal sources, resulting in a decrease in resolution. This has a significant impact on applications such as high-resolution imaging or multi-target tracking, which in turn limits the performance of the system. Therefore, how to improve the accuracy and resolution of compact array DOA estimation to promote array miniaturization has become a problem that needs to be solved urgently by those skilled in the art. Summary of the invention

[0004] The purpose of the present invention is to provide a compact array DOA estimation method based on alternating projection and golden section search method to effectively improve the estimation accuracy and resolution.

[0005] In order to achieve the above tasks, the present invention adopts the following technical solutions:

[0006] A compact array DOA estimation method based on alternating projection and golden section search method, comprising:

[0007] S1, input radar array signal S, initialize the total number of search arrival directions K, estimated arrival direction vector θ, current estimated arrival direction sequence number k and k's counter Count k ;

[0008] S2, using the golden section search method, preliminarily estimates the kth wave arrival direction;

[0009] S3, using the alternating projection algorithm, dynamically updates the estimated direction of arrival vector θ;

[0010] S4, determine whether the current estimated arrival direction number k is greater than the total number of searched arrival directions K; if so, the estimation of all arrival directions has been completed and the search ends; otherwise, update the current estimated arrival direction number k and the counter Count k And return to S2.

[0011] Furthermore, in S1, θ is initialized to a zero matrix, k, Count k The initial values ​​are all set to 1, and the initialization amplitude formula is: θ = [0]; k = 1; Count k =1; the total number K of search arrival directions is known or has been estimated using existing numerical detection methods.

[0012] Furthermore, the k-th direction of arrival is preliminarily estimated by using the golden section search method in S2, including:

[0013] S2.1, calculate the first search point θ 11 and θ 12 , the specific formula is:

[0014] θ 11 =θ l1 +(1-λ)(θ u1 -θ l1 )

[0015] θ 12 =θ l1 +λ(θ u1 -θ l1 )

[0016] Among them, λ represents the search step size, θ u1 represents the upper limit of the first angle estimate, and represents θ l1 The first angle estimates the lower limit;

[0017] S2.2, determine the search interval length |θ u1 -θ l1 Is it greater than the threshold ε? 1 :

[0018] If the interval length |θ u1 -θ l1 |Less than or equal to the threshold ε 1 , the preliminary estimation of the kth wave arrival direction has been completed, the loop ends, and jumps to S2.5;

[0019] If the interval length |θ u1 -θ l1 |Greater than the threshold ε 1 , then jump to S2.3 and continue to narrow the search range;

[0020] S2.3, substitute the first search point θ 11 and θ 12 Calculate the corresponding cost value F(θ 11 ) and F(θ 12 ), the specific formula is:

[0021] F(θ t )=tr{[IA(θ now1 )(A H (θ now1 )A(θ now1 )) -1 A(θ now1 )]R}

[0022] where θ t represents the search point, the matrix superscript H represents the conjugate transpose operation performed on the matrix, [] -1 Indicates the inverse operation of the matrix, tr{} indicates the trace of the matrix, that is, the sum of the diagonal elements of the main diagonal of the matrix, the same below; I indicates the identity matrix, R indicates the complex covariance matrix, A(θ now1 ) represents the steering matrix corresponding to the direction of arrival, and the specific formula is:

[0023]

[0024] Where S represents the array signal, e is a natural constant, and N represents the number of snapshots. represents the imaginary symbol, λ is the wavelength of the carrier, the distance vector D of each array element relative to the reference array element and the current wave arrival direction estimation vector θ now1 The specific formula is:

[0025] D=[0,d,2d,...,(M-1)d] T

[0026]

[0027] The matrix superscript T indicates that the matrix is ​​transposed, and d is the array element spacing, which is M is the number of array elements, [θ t ,θ] indicates that the estimated arrival direction vector θ has an element and θ t The vector formed;

[0028] S2.4, determine the cost value F(θ 11 ) is greater than the cost value F(θ 12 );

[0029] If the cost value F(θ 11 ) is greater than the cost value F(θ 12 ), then update the lower limit of the search interval and the first search point, and then jump to S2.2;

[0030] If the cost value F(θ 11 ) is less than or equal to the cost value F(θ 12 ), then update the upper limit of the search interval and the first search point, and then jump to S2.2;

[0031] S2.5, calculate the final search result, the specific formula is:

[0032]

[0033] Among them, θ( k ) represents the estimated result of the k-th arrival direction, which is used as the k-th element of the estimated arrival direction vector θ.

[0034] Furthermore, the specific formula for updating the lower limit of the search interval and the first search point is as follows:

[0035] θ l1 =θ 11

[0036] θ 11 =θ 12

[0037] θ 12 =θ l1 +λ(θ u1 -θ l1 )

[0038] The specific formula for updating the upper limit of the search interval and the first search point is as follows:

[0039] θ u1 =θ 12

[0040] θ 12 =θ 11

[0041] θ 11 =θ l1 +(1-λ)(θ u1 -θ l1 ).

[0042] Furthermore, the search step length λ=0.618, the first angle estimation upper limit θ u1 The initial value is 89, and the lower limit of the first angle estimation is θ l1 The initial value is -90.

[0043] Furthermore, the method of S3 dynamically updates the estimated direction of arrival vector θ by using an alternating projection algorithm, including:

[0044] S3.1, initialize the current estimated direction of arrival sum θ sn , the sum of historical wave arrival directions θso And alternately estimate the direction of arrival number j, Count j is the counter of j; the specific initialization assignment formula is:

[0045]

[0046] θ so =+∞

[0047] j=0

[0048] Count j =0

[0049] S3.2, determine the current estimated total direction of arrival θ sn and the sum of historical arrival directions θ so The difference θ sn -θ so Is it greater than the threshold ε? 2 :

[0050] If the difference |θ sn -θ so |Less than or equal to the threshold ε 2 , the arrival direction estimation of k signals has been completed, the loop ends, and jumps to S4; if the difference |θ sn -θ so |Greater than the threshold ε 2 , then jump to S3.3 and continue to accurately estimate k wave arrival directions;

[0051] S3.3, determine whether the alternate estimated direction of arrival number j is greater than the current estimated direction of arrival number k:

[0052] If the alternate estimated direction of arrival number j is greater than the current estimated direction of arrival number k, the alternate estimation of the first k signals in the array signal S has been completed, the loop ends, and the process jumps to S3.2;

[0053] If the alternate estimated direction of arrival number j is less than the current estimated direction of arrival number k, then update j = Count j +1, Count j =j, then jump to S3.4;

[0054] S3.4, using the golden section search method, estimate the j-th direction of arrival;

[0055] S3.5, update the sum of historical arrival directions θ so =θ sn and the sum of the current estimated direction of arrival Then return to S3.2;

[0056] Furthermore, the method of using the golden section search method in S3.4 to estimate the j-th direction of arrival includes:

[0057] S3.4.1, calculate the second search point θ 21 and θ 22 , the specific formula is:

[0058] θ 21 =θ l2 +(1-λ)(θ u2 -θ l2 )

[0059] θ 22 =θ l2 +λ(θ u2 -θ l2 )

[0060] Among them, θ u2 represents the upper limit of the second angle estimate, θ l2 represents the lower limit of the second angle estimation;

[0061] S3.4.2, determine the search interval length θ u2 -θ l2 Is it greater than the threshold ε? 1 :

[0062] If the interval length θ u2 -θ l2 |Less than or equal to the threshold ε 1 , the jth updated direction of arrival estimation has been completed, the loop ends, and jumps to S3.4.5;

[0063] If the interval length θ u2 -θ l2 |Greater than the threshold ε 1 , then jump to S3.4.3 and continue to narrow the search range;

[0064] S3.4.3, substitute the search point θ 21 and θ 22 Calculate the corresponding cost value F(θ 21 ) and F(θ 22 ), the specific formula is:

[0065] F(θ t )=tr{[IA(θ now2 )(A H (θ now2 )A(θ now2 )) -1 A(θ now2 )]R}

[0066] A(θ now2) represents the steering matrix corresponding to the current estimated direction of arrival, and the specific formula is:

[0067]

[0068] Among them, the current arrival direction estimation vector θ now2 The specific formula is:

[0069]

[0070] In the above formula, θ (-j) represents the remaining elements after removing the jth element in θ, [θ t ,θ (-j) ] represents θ t and the matrix formed by the remaining elements;

[0071] S3.4.4, determine the cost value F(θ 21 ) is greater than the cost value F(θ 22 );

[0072] If the cost value F(θ 21 ) is greater than the cost value F(θ 22 ), then update the lower limit of the search interval and the second search point, and then jump to S3.4.2;

[0073] If the cost value F(θ 21 ) is less than or equal to the cost value F(θ 22 ), then update the upper limit of the search interval and the second search point, and then jump to S2.2;

[0074] S3.4.5, calculate the final search result, the specific formula is:

[0075]

[0076] Among them, θ( j ) represents the estimated result of the j-th arrival direction, which is the j-th element of the estimated arrival direction vector θ.

[0077] Furthermore, the specific formula for updating the lower limit of the search interval and the second search point is as follows:

[0078] θ l2 =θ 21

[0079] θ 21 =θ 22

[0080] θ 22 =θ l2 +λ(θ u2 -θ l2 )

[0081] The specific formula for updating the upper limit of the search interval and the second search point is as follows:

[0082] θ u2 =θ 22

[0083] θ 22 =θ 21

[0084] θ 21 =θ l2 +(1-λ)(θ u2 -θ l2 ).

[0085] Furthermore, the second angle estimation upper limit θ u2 , the second angle estimation lower limit θ l2 The calculation formula is:

[0086]

[0087] Among them, θ (j) Represents the estimated result of the j-th direction of arrival.

[0088] A terminal device comprises a processor, a memory and a computer program stored in the memory; when the processor executes the computer program, the compact array DOA estimation method based on alternating projection and golden section search method is implemented.

[0089] A computer-readable storage medium stores a computer program; when the computer program is executed by a processor, the compact array DOA estimation method based on alternating projection and golden section search method is implemented.

[0090] Compared with the prior art, the present invention has the following technical features:

[0091] 1. The present invention uses alternating projection to optimize the wave direction vector in each iteration, so that the estimation at each step can better approximate the true value. The core of this method is to decompose the wave direction estimation process into multiple alternating projection steps, so that it can achieve efficient convergence in a complex multi-source signal environment. This step-by-step optimization strategy not only improves the stability of the estimation, but also effectively reduces the error caused by noise, ultimately achieving the purpose of improving resolution. It also takes into account the unique characteristics of compact arrays, and alternating projection further enhances the adaptability of the invention in practical applications;

[0092] 2. The present invention uses the golden section search method, based on the estimation results of the previous round of DOA, to intelligently set the search range for update, ensuring that the new round of search is more efficient. Compared with the traditional global search method, the golden section search method significantly reduces the amount of calculation while maintaining strict control of the search accuracy. This method not only improves the accuracy of DOA estimation, but also shortens the required calculation time, meeting the needs of real-time processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] Figure 1 It is a flow chart of the compact array DOA estimation method based on alternating projection and golden section search method;

[0094] Figure 2 Flow chart of searching direction of arrival for golden section search method;

[0095] Figure 3 A comparison diagram of simulation results of the present method and the MUSIC method when the incident angle is [10.46°, -5.46°] in Example 1 of the present invention;

[0096] Figure 4 A comparison diagram of simulation results of the present method and the MUSIC method when the incident angle is [10.46°, 5.46°] in Example 1 of the present invention;

[0097] Figure 5 This is a comparison diagram of the root mean square error of the present method and the MUSIC method as a function of the signal-to-noise ratio when the incident angle is [10.46°] in Example 2 of the present invention;

[0098] Figure 6 This is a comparison diagram of the root mean square error of the present method and the MUSIC method as a function of the signal-to-noise ratio when the incident angle is [30.46°, -20.63°] in Example 2 of the present invention;

[0099] Figure 7 This is a comparison chart of the root mean square error of the present method and the MUSIC method versus the angle interval in Example 3 of the present invention. DETAILED DESCRIPTION

[0100] See attached Figure 1 and Figure 2 The present invention provides a compact array DOA estimation method based on alternating projection and golden section search method, comprising the following steps:

[0101] S1, input radar array signal S, initialize the total number of search arrival directions K, estimated arrival direction vector θ, current estimated arrival direction sequence number k and k's counter Count k ; Among them, θ is initialized to a 0 matrix, k, Count kThe initial values ​​are all set to 1, and the initialization amplitude formula is: θ = [0]; k = 1; Count k =1; the total number K of search arrival directions is known or has been estimated using existing numerical detection methods.

[0102] S2, using the golden section search method, preliminarily estimates the kth wave arrival direction.

[0103] S2.1, calculate the first search point θ 11 and θ 12 , the specific formula is:

[0104] θ 11 =θ l1 +(1-λ)(θ u1 -θ l1 )

[0105] θ 12 =θ l1 +λ(θ u1 -θ l1 )

[0106] Among them, the search step size λ = 0.618, the first angle estimation upper limit θ u1 The initial value is 89, and the lower limit of the first angle estimation is θ l1 The initial value is -90.

[0107] S2.2, determine the search interval length |θ u1 -θ l1 Is it greater than the threshold ε? 1 :

[0108] If the interval length |θ u1 -θ l1 |Less than or equal to the threshold ε 1 , the preliminary estimation of the kth wave arrival direction has been completed, the loop ends, and jumps to S2.5;

[0109] If the interval length |θ u1 -θ l1 |Greater than the threshold ε 1 , then jump to S2.3 and continue to narrow the search range;

[0110] In this embodiment, the threshold ε 1 The value of is 0.001.

[0111] S2.3, substitute the first search point θ 11 and θ 12 Calculate the corresponding cost value F(θ 11 ) and F(θ 12 ), the specific formula is:

[0112] F(θt )=tr{[IA(θ now1 )(A H (θ now1 )A(θ now1 )) -1 A(θ now1 )]R}

[0113] where θ t represents the search point, the matrix superscript H represents the conjugate transpose operation performed on the matrix, [] -1 Indicates the inverse operation of the matrix, tr{} indicates the trace of the matrix, that is, the sum of the diagonal elements of the main diagonal of the matrix, the same below; I indicates the identity matrix, R indicates the complex covariance matrix, A(θ now1 ) represents the steering matrix corresponding to the direction of arrival, and the specific formula is:

[0114]

[0115] Where S represents the array signal, e is a natural constant, N represents the number of snapshots, i represents the imaginary symbol, λ is the wavelength of the carrier, the distance vector D of each array element relative to the reference array element and the current direction of arrival estimation vector θ now1 The specific formula is:

[0116] D=[0,d,2d,...,(M-1)d] T

[0117]

[0118] The matrix superscript T indicates that the matrix is ​​transposed, and d is the array element spacing, which is M is the number of array elements, [θ t ,θ] indicates that the estimated arrival direction vector θ has an element and θ t The vector composed.

[0119] S2.4, determine the cost value F(θ 11 ) is greater than the cost value F(θ 12 ); if the cost value F(θ 11 ) is greater than the cost value F(θ 12 ), then update the lower limit of the search interval and the first search point, and then jump to S2.2;

[0120] The specific formulas for updating the lower limit of the search interval and the first search point are as follows:

[0121] θ l1 =θ 11

[0122] θ 11 =θ 12

[0123] θ 12 =θ l1 +λ(θ u1 -θ l1 )

[0124] In the above update and assignment formulas, the equal sign represents an assignment operation, that is, the value on the right side of the equal sign is assigned to the parameter on the left side of the equal sign.

[0125] If the cost value F(θ 11 ) is less than or equal to the cost value F(θ 12 ), then update the upper limit of the search interval and the first search point, and then jump to S2.2;

[0126] The specific formulas for updating the upper limit of the search interval and the first search point are as follows:

[0127] θ u1 =θ 12

[0128] θ 12 =θ 11

[0129] θ 11 =θ l1 +(1-λ)(θ u1 -θ l1 )

[0130] S2.5, calculate the final search result, the specific formula is:

[0131]

[0132] Among them, θ( k ) represents the estimated result of the k-th arrival direction, which is used as the k-th element of the estimated arrival direction vector θ.

[0133] S3, using the alternating projection algorithm, dynamically updates the estimated direction of arrival vector θ.

[0134] S3.1, initialize the current estimated direction of arrival sum θ sn , the sum of historical wave arrival directions θ so And alternately estimate the direction of arrival number j, Count j is the counter of j; the specific initialization assignment formula is:

[0135]

[0136] θ so =+∞

[0137] j=0

[0138] Countj =0

[0139] S3.2, determine the current estimated total direction of arrival θ sn and the sum of historical arrival directions θ so The difference θ sn -θ so Is it greater than the threshold ε? 2 :

[0140] If the difference |θ sn -θ so |Less than or equal to the threshold ε 2 , the arrival direction estimation of k signals has been completed, the loop ends, and jumps to S4; if the difference |θ sn -θ so |Greater than the threshold ε 2 , then jump to S3.3 and continue to accurately estimate k wave arrival directions;

[0141] In this embodiment, the threshold ε 2 Set to 0.01.

[0142] S3.3, determine whether the alternate estimated direction of arrival number j is greater than the current estimated direction of arrival number k:

[0143] If the alternate estimated direction of arrival number j is greater than the current estimated direction of arrival number k, the alternate estimation of the first k signals in the array signal S has been completed, the loop ends, and the process jumps to S3.2;

[0144] If the alternate estimated direction of arrival number j is less than the current estimated direction of arrival number k, then update j = Count j +1, Count j =j, then jump to S3.4;

[0145] S3.4, using the golden section search method, estimate the j-th direction of arrival.

[0146] S3.4.1, calculate the second search point θ 21 and θ 22 , the specific formula is:

[0147] θ 21 =θ l2 +(1-λ)(θ u2 -θ l2 )

[0148] θ 22 =θ l2 +λ(θ u2 -θ l2 )

[0149] Among them, the second angle estimation upper limit and the second angle estimation lower limit are specifically formulated as follows:

[0150]

[0151] Among them, θ (j) Represents the estimated result of the j-th direction of arrival.

[0152] S3.4.2, determine the search interval length |θ u2 -θ l2 Is it greater than the threshold ε? 1 :

[0153] If the interval length |θ u2 -θ l2 |Less than or equal to the threshold ε 1 , the jth updated direction of arrival estimation has been completed, the loop ends, and jumps to S3.4.5;

[0154] If the interval length |θ u2 -θ l2 |Greater than the threshold ε 1 , then jump to S3.4.3 and continue to narrow the search range;

[0155] S3.4.3, substitute the search point θ 21 and θ 22 Calculate the corresponding cost value F(θ 21 ) and F(θ 22 ), the specific formula is:

[0156] F(θ t )=tr{[IA(θ now2 )(A H (θ now2 )A(θ now2 )) -1 A(θ now2 )]R}

[0157] A(θ now2 ) represents the steering matrix corresponding to the current estimated direction of arrival, and the specific formula is:

[0158]

[0159] Among them, the current arrival direction estimation vector θ now2 The specific formula is:

[0160]

[0161] In the above formula, θ (-j) represents the remaining elements after removing the jth element in θ, [θ t ,θ (-j) ] represents θt The matrix formed with the remaining elements.

[0162] S3.4.4, determine the cost value F(θ 21 ) is greater than the cost value F(θ 22 ); if the cost value F(θ 21 ) is greater than the cost value F(θ 22 ), then update the lower limit of the search interval and the second search point, and then jump to S3.4.2;

[0163] The specific formula for updating the lower limit of the search interval and the second search point is as follows:

[0164] θ l2 =θ 21

[0165] θ 21 =θ 22

[0166] θ 22 =θ l2 +λ(θ u2 -θ l2 )

[0167] If the cost value F(θ 21 ) is less than or equal to the cost value F(θ 22 ), then update the upper limit of the search interval and the second search point, and then jump to S2.2;

[0168] The specific formulas for updating the upper limit of the search interval and the second search point are as follows:

[0169] θ u2 =θ 22

[0170] θ 22 =θ 21

[0171] θ 21 =θ l2 +(1-λ)(θ u2 -θ l2 )

[0172] S3.4.5, calculate the final search result, the specific formula is:

[0173]

[0174] Among them, θ( j ) represents the estimated result of the j-th arrival direction, which is the j-th element of the estimated arrival direction vector θ.

[0175] S3.5, update the sum of historical arrival directions θ so=θ sn and the sum of the current estimated direction of arrival Then return to S3.2.

[0176] S4, determine whether the current estimated arrival direction number k is greater than the total number K of searched arrival directions:

[0177] If the current estimated arrival direction number k is greater than the total number of searched arrival directions K, the estimation of all arrival directions has been completed and the search ends;

[0178] If the current estimated arrival direction number k is less than or equal to the total number of searched arrival directions K, execute k = Count k +1, Count k =k, jump to S2.

[0179] Example:

[0180] In the embodiment of the present invention, the model of the simulation array signal S is:

[0181]

[0182] Among them, θ k represents the kth wave arrival angle, K r represents the number of arrival angles, A k (θ k ) represents the steering vector corresponding to the kth wave arrival angle, n is 0, and the variance is σ 2 Additive Gaussian white noise, s k represents the kth signal source, and the specific formula is:

[0183]

[0184] D=[0,d,2d,...,(M-1)d] T

[0185] s k =[s k (1),s k (2),...,s k (i),...,s k (N)]

[0186]

[0187] Where d is the array element spacing, which is M is the number of array elements, which is 4. The signal amplitude is A s =1, sampling rate F s =10kHz, f 1 is the starting frequency of the signal, f 2 is the stop frequency of the signal, is the initial phase of the signal, and the number of snapshots N=2000.

[0188] In the embodiment of the present invention, there are two signal sources. The starting frequency f of signal source 1 is 11 =1kHz, stop frequency f 12 =1kHz, initial phase The starting frequency f of signal source 2 21 =1.2kHz, end frequency f 22 =0.8kHz, initial phase

[0189] The size of the additive white Gaussian noise n is determined by the interference-to-noise ratio SNR, and the specific formula is:

[0190]

[0191] In order to accurately evaluate the effect of DOA estimation, the root mean square error of DOA is used as the evaluation index. The specific formula is:

[0192]

[0193] Among them, θ r (k) is the true direction of arrival, θ (k) To estimate the direction of arrival.

[0194] Embodiment 1:

[0195] In this embodiment, two groups of incident angles are set, namely [10.46°, -5.46°] and [10.46°, 5.46°], and the incident signals are signal source 1 and signal source 2 respectively. Other simulation parameters are the same as the initial parameters given in the example. The compact array DOA estimation method based on alternating projection and golden section search method and the MUSIC method are used to analyze the two groups of array simulation signals. The results are shown in FIG. Figure 3 and Figure 4 When the incident angle is [10.46°, -5.46°], both methods can distinguish the target angle. At this time, the root mean square error of the compact array DOA estimation method based on the alternating projection and golden section search method is 0.111, and the root mean square error of the MUSIC method is 0.502, which proves that the method of the present invention can effectively perform DOA estimation and improve the accuracy; when the incident angle is [10.46°, 5.46°], the MUSIC method has only one spectral peak and cannot distinguish two targets, while the compact array DOA estimation method based on the alternating projection and golden section search method can accurately identify the target, which proves that the method of the present invention can effectively improve the resolution.

[0196] Embodiment 2:

[0197] 3100 Monte Carlo experiments were conducted, with the signal-to-noise ratio ranging from -10dB to 20dB. Two groups of incident angles were set, namely [10.46°] (signal source 1) and [30.46°, -20.63°] (signal source 1 and signal source 2). Other simulation parameters were the same as the initial parameters given in the example. The compact array DOA estimation method based on alternating projection and golden section search method and the MUSIC method were used to analyze the simulation signal, and the root mean square error results were obtained as shown in the figure. Figure 5 and Figure 6 As shown. Figure 5 and Figure 6 Analysis shows that under the same parameters, the root mean square error of the compact array DOA estimation method based on alternating projection and golden section search method is significantly better than that of the traditional MUSIC algorithm. The above results show that the compact array DOA estimation method based on alternating projection and golden section search method proposed in the present invention can improve the accuracy.

[0198] Embodiment 3:

[0199] 1500 Monte Carlo experiments were conducted. The angle interval between signal source 1 and signal source 2 ranged from 1° to 15°. The other simulation parameters were the same as the initial parameters given in the example. The compact array DOA estimation method based on alternating projection and golden section search method and the MUSIC method were used to analyze the simulation signal. The root mean square error results were obtained as shown in the figure. Figure 7 As shown. Figure 7 The analysis shows that when the interval angle is greater than 1°, the root mean square error of the compact array DOA estimation method based on alternating projection and golden section search method decreases significantly, and reaches a lower level when the angle interval is greater than 3°, which is significantly smaller than the traditional MUSIC algorithm. Figure 3 , Figure 4 and Figure 7 ,The simulation results show that the method described in the ,presentation can improve the resolution when estimating DOA.

[0200] The above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. These modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application, and should all be included in the protection scope of the present application.

Claims

1. A compact array DOA estimation method based on alternating projection and golden section search method, characterized in that: include: S1, input radar array signal S, initialize the total number of search arrival directions K, estimated arrival direction vector θ, current estimated arrival direction sequence number k and k's counter Count k ; S2, using the golden section search method, preliminarily estimates the kth wave arrival direction; S3, using the alternating projection algorithm, dynamically updates the estimated direction of arrival vector θ; S4, determine whether the current estimated arrival direction number k is greater than the total number of searched arrival directions K; if so, the estimation of all arrival directions has been completed and the search ends; otherwise, update the current estimated arrival direction number k and the counter Count k And return to S2.

2. The compact array DOA estimation method based on alternating projection and golden section search method according to claim 1, characterized in that: In S1, θ is initialized to a zero matrix, k, Count k The initial values ​​are all set to 1, and the initialization amplitude formula is: θ = [0]; k = 1; Count k =1; the total number K of search arrival directions is known or has been estimated using existing numerical detection methods.

3. The compact array DOA estimation method based on alternating projection and golden section search method according to claim 1, characterized in that: S2 uses the golden section search method to preliminarily estimate the k-th wave arrival direction, including: S2.1, calculate the first search point θ 11 and θ 12 , the specific formula is: i 11 =θ l1 +(1-λ)(θ) u1 -θ l1 ) i 12 =θ l1 +λ(θ u1 -θ l1 ) Among them, λ represents the search step size, θ u1 represents the upper limit of the first angle estimate, and represents θ l1 The first angle estimates the lower limit; S2.2, determine the search interval length |θ u1 -θ l1 |Is it greater than the threshold ε1: If the interval length |θ u1 -θ l1 | is less than or equal to the threshold ε1, the preliminary estimation of the kth wave arrival direction has been completed, the loop ends, and jumps to S2.5; If the interval length θ u1 -θ l1 If it is greater than the threshold ε1, then jump to S2.3 and continue to narrow the search range; S2.3, substitute the first search point θ 11 and θ 12 Calculate the corresponding cost value F(θ 11 ) and F(θ 12 ), the specific formula is: F(θ t )=tr{[IA(θ now1 )(A H (i now1 )A(θ now1 )) -1 A(θ now1 )]R} where θ t represents the search point, the matrix superscript H represents the conjugate transpose operation performed on the matrix, [] -1 Indicates the inverse operation of the matrix, tr{} indicates the trace of the matrix, that is, the sum of the diagonal elements of the main diagonal of the matrix, the same below; I indicates the identity matrix, R indicates the complex covariance matrix, A(θ now1 ) represents the steering matrix corresponding to the direction of arrival, and the specific formula is: Where S represents the array signal, e is a natural constant, N represents the number of snapshots, i represents the imaginary symbol, λ is the wavelength of the carrier, the distance vector D of each array element relative to the reference array element and the current direction of arrival estimation vector θ now1 The specific formula is: D=[0,d,2d,...,(M-1)d] T The matrix superscript T indicates that the matrix is ​​transposed, and d is the array element spacing, which is M is the number of array elements, [θ t ,θ] indicates that the estimated arrival direction vector θ has an element and θ t The vector formed; S2.4, determine the cost value F(θ 11 ) is greater than the cost value F(θ 12 ); If the cost value F(θ 11 ) is greater than the cost value F(θ 12 ), then update the lower limit of the search interval and the first search point, and then jump to S2.2; If the cost value F(θ 11 ) is less than or equal to the cost value F(θ 12 ), then update the upper limit of the search interval and the first search point, and then jump to S2.2; S2.5, calculate the final search result, the specific formula is: Among them, θ( k ) represents the estimated result of the k-th arrival direction, which is used as the k-th element of the estimated arrival direction vector θ.

4. The compact array DOA estimation method based on alternating projection and golden section search method according to claim 3, characterized in that: The specific formula for updating the lower limit of the search interval and the first search point is as follows: i l1 =θ 11 i 11 =θ 12 i 12 =θ l1 +λ(θ u1 -θ l1 ) The specific formula for updating the upper limit of the search interval and the first search point is as follows: i u1 =θ 12 i 12 =θ 11 i 11 =θ l1 +(1-λ)(θ) u1 -θ l1 )。 5. The compact array DOA estimation method based on alternating projection and golden section search method according to claim 1, characterized in that: S3 uses the alternating projection algorithm to dynamically update the estimated arrival direction vector θ, including: S3.1, initialize the current estimated direction of arrival sum θ sn , the sum of historical wave arrival directions θ so And alternately estimate the direction of arrival number j, Count j is the counter of j; the specific initialization assignment formula is: i so =+∞ j＝0 Count j =0 S3.2, determine the current estimated total direction of arrival θ sn and the sum of historical arrival directions θ so The difference θ sn -θ so |Is it greater than the threshold ε2: If the difference |θ sn -θ so | is less than or equal to the threshold ε2, the arrival direction estimation of k signals has been completed, the loop ends, and jumps to S4; if the difference |θ sn -θ so | is greater than the threshold ε2, jump to S3.3 and continue to accurately estimate k wave arrival directions; S3.3, determine whether the alternate estimated direction of arrival number j is greater than the current estimated direction of arrival number k: If the alternate estimated direction of arrival number j is greater than the current estimated direction of arrival number k, the alternate estimation of the first k signals in the array signal S has been completed, the loop ends, and the process jumps to S3.2; If the alternate estimated direction of arrival number j is less than the current estimated direction of arrival number k, then update j = Count j +1, Count j =j, then jump to S3.4; S3.4, using the golden section search method, estimate the j-th direction of arrival; S3.5, update the sum of historical arrival directions θ so =θ sn and the sum of the current estimated direction of arrival Then return to S3.

2.

6. The compact array DOA estimation method based on alternating projection and golden section search method according to claim 5, characterized in that: S3.4 uses the golden section search method to estimate the j-th direction of arrival, including: S3.4.1, calculate the second search point θ 21 and θ 22 , the specific formula is: i 21 =θ l2 +(1-λ)(θ) u2 -θ l2 ) i 22 =θ l2 +λ(θ u2 -θ l2 ) Among them, θ u2 represents the upper limit of the second angle estimation, θ l2 represents the lower limit of the second angle estimation; S3.4.2, determine the search interval length |θ u2 -θ l2 |Is it greater than the threshold ε1: If the interval length |θ u2 -θ l2 | is less than or equal to the threshold ε1, the estimation of the jth updated direction of arrival has been completed, the loop ends, and jumps to S3.4.5; If the interval length |θ u2 -θ l2 | is greater than the threshold ε1, jump to S3.4.3 and continue to narrow the search range; S3.4.3, substitute the search point θ 21 and θ 22 Calculate the corresponding cost value F(θ 21 ) and F(θ 22 ), the specific formula is: F(θ t )=tr{[IA(θ now2 )(A H (i now2 )A(θ now2 )) -1 A(θ now2 )]R} A(θ now2 ) represents the steering matrix corresponding to the current estimated direction of arrival, and the specific formula is: Among them, the current arrival direction estimation vector θ now2 The specific formula is: In the above formula, θ (-j) represents the remaining elements after removing the jth element in θ, [θ t ,θ (-j) ] represents θ t and the matrix formed by the remaining elements; S3.4.4, determine the cost value F(θ 21 ) is greater than the cost value F(θ 22 ); If the cost value F(θ 21 ) is greater than the cost value F(θ 22 ), then update the lower limit of the search interval and the second search point, and then jump to S3.4.2; If the cost value F(θ 21 ) is less than or equal to the cost value F(θ 22 ), then update the upper limit of the search interval and the second search point, and then jump to S2.2; S3.4.5, calculate the final search result, the specific formula is: Among them, θ( j ) represents the estimated result of the j-th arrival direction, which is the j-th element of the estimated arrival direction vector θ.

7. The compact array DOA estimation method based on alternating projection and golden section search method according to claim 6, characterized in that: The specific formula for updating the lower limit of the search interval and the second search point is as follows: i l2 =θ 21 i 21 =θ 22 i 22 =θ l2 +λ(θ u2 -θ l2 ) The specific formula for updating the upper limit of the search interval and the second search point is as follows: i u2 =θ 22 i 22 =θ 21 i 21 =θ l2 +(1-λ)(θ) u2 -θ l2 )。 8. The compact array DOA estimation method based on alternating projection and golden section search method according to claim 6, characterized in that: The second angle estimation upper limit θ u2 , the second angle estimation lower limit θ l2 The calculation formula is: Among them, θ (j) Represents the estimated result of the j-th direction of arrival.

9. A terminal device comprising a processor, a memory and a computer program stored in the memory; characterized in that: When the processor executes the computer program, the compact array DOA estimation method based on the alternating projection and golden section search method according to any one of claims 1 to 9 is implemented.

10. A computer-readable storage medium, wherein a computer program is stored in the medium; characterized in that: When the computer program is executed by a processor, the compact array DOA estimation method based on alternating projection and golden section search method according to any one of claims 1 to 9 is implemented.

Images