Multi-objective passive array direction-finding and amplitude measurement method based on iterative maximum sparse convex optimization

Abstract

The application discloses a multi-target passive array direction-finding and amplitude-measuring method based on iterative maximum sparse convex optimization, and comprises the following steps: S1, a multi-target passive array sparse signal processing model is established to obtain a received signal matrix Y; S2, the received signal matrix Y is compressed by using singular value decomposition, and a compressed signal vector y is obtained after rearrangement v ; S3, a convex optimization mathematical model with a norm constraint is constructed, and a solution x j of the convex optimization problem is obtained by solving the convex optimization problem, and the solution is recorded as S4, an IMSC optimization model is constructed, and an optimization result x opt is obtained by iterative calculation; and S5, the optimization result x opt is used to obtain multi-target angle and amplitude estimation results. The application guarantees that a target function is a convex function, maximizes non-convexity of a penalty function, improves sparsity of a solution and reduces amplitude estimation error, and realizes high-precision multi-target direction-finding and amplitude-measuring.

Description

TECHNICAL FIELD

[0001] The present application relates to passive array direction finding and amplitude measurement, in particular to a multi-target passive array direction finding and amplitude measurement method based on iterative maximum sparse convex optimization. BACKGROUND

[0002] Passive direction finding and amplitude measurement technology is widely used in radio positioning, radio astronomy, communication, navigation, radar, sonar and many other fields, and is a key technology for Internet of Things, electronic reconnaissance and electromagnetic environment perception. In the past decade, with the maturity of digital signal processing technology, passive direction finding and amplitude measurement technology has been widely developed and applied. Among them, using a passive array for multi-target direction finding and amplitude measurement has gradually become the mainstream solution.

[0003] The classical beamforming method (CBF) is a traditional passive array direction finding and amplitude measurement method. This method adjusts the weighting vector of the array antenna output signal to form beams with different directions, and the beam direction corresponding to the maximum output power is the incident direction of the target signal. This method can obtain the angle information and power information of the target at the same time. However, when there are multiple targets, the main lobe of a target is easily disturbed by the side lobes of other targets, resulting in inaccurate angle estimation and significantly elevated amplitude estimation. Therefore, the CBF method is less commonly used in multi-target direction finding and amplitude measurement.

[0004] In the early 21st century, researchers applied the sparse recovery theory in compressive sensing to the field of array signal processing, solving the above problems. A typical method among them is the l1-norm-based sparse singular value decomposition method (l1-SVD). This method introduces an l1-norm penalty to encourage sparsity in the solution, greatly improving the resolution of multiple targets. However, due to the introduction of the l1-norm penalty, the amplitude of the solution is also underestimated. Therefore, developing a new method that can simultaneously achieve high-precision direction finding and amplitude measurement for multiple targets is a problem that must be solved to further develop passive direction finding and amplitude measurement technology. SUMMARY

[0005] The present application aims to overcome the shortcomings of the prior art and provide a multi-target passive array direction finding and amplitude measurement method based on iterative maximum sparse convex optimization. The iterative maximum sparse convex optimization algorithm is used to improve the penalty function in the sparse recovery model, changing it from a traditional convex penalty function to a non-convex penalty function. While ensuring that the objective function is convex, the non-convexity of the penalty function is maximized to improve the sparsity of the solution and reduce the amplitude estimation error, achieving high-precision direction finding and amplitude measurement for multiple targets.

[0006] The purpose of the present application is achieved by the following technical solution: a multi-target passive array direction finding and amplitude measurement method based on iterative maximum sparse convex optimization, comprising the following steps:

[0007] S1. Establish a multi-target passive array sparse signal processing model to obtain a received signal matrix Y;

[0008] S2. Compress the received signal matrix Y using singular value decomposition, and rearrange to obtain a compressed signal vector y v ;

[0009] S3. Construct a convex optimization mathematical model with a norm constraint, and solve the convex optimization problem to obtain a solution x of the convex optimization problem j , and record the solution as

[0010] S4. Construct an IMS C optimization model, and iteratively calculate to obtain an optimization result x opt ;

[0011] S5. Obtain multi-target angle and amplitude estimation results by using the optimization result x opt .

[0012] The present application has the following beneficial effects: the present application uses an iterative maximum sparse convex optimization algorithm, improves the penalty function in the sparse recovery model, changes the traditional convex penalty function into a non-convex penalty function, maximizes the non-convexity of the penalty function while ensuring that the objective function is a convex function (the solution is a global optimal solution, i.e., the direction finding and amplitude estimation result is unique), to improve the sparsity of the solution and reduce the amplitude estimation error, make up for the shortcomings of the traditional sparse recovery method that underestimates the signal amplitude, and realize high-precision direction finding and amplitude estimation of multi-targets. BRIEF DESCRIPTION OF DRAWINGS

[0013] Figure 1 is a flow chart of the method of the present application;

[0014] Figure 2 is a flow chart of the IMS C algorithm;

[0015] Figure 3 is a schematic diagram of the direction finding and amplitude estimation results of multi-targets by the present method and the prior art;

[0016] Figure 4 is a schematic diagram of the amplitude deviation curve with the change of the signal-to-noise ratio. DETAILED DESCRIPTION

[0017] The technical solutions of the present application will be further described in detail below with reference to the accompanying drawings, but the protection scope of the present application is not limited to the following description.

[0018] As shown in Figure 1 , the multi-target passive array direction finding and amplitude estimation method based on the iterative maximum sparse convex optimization includes the following steps:

[0019] (1) Establish a multi-target passive array sparse signal processing model to obtain a received signal matrix Y;

[0020] Let Q denote the number of target radiation sources, and Q targets respectively generate Q far-field narrowband signals incident to the receiving array. Assume that the receiving array is a one-dimensional uniform linear array (ULA), the number of array elements is N (N > Q), the element spacing is half the wavelength of the narrowband signal carrier frequency, and each array element antenna is an omnidirectional antenna. The spatial electromagnetic wave signal is received by the ULA, and after pre-processing by the radio frequency front end, digital intermediate frequency quadrature sampling is performed. Let the sampling point number (also known as the frame number or snapshot number) be K, and the sampling interval be T s , then the time sampling vector is represented as t = [0, T s , 2T s ,..., (K-1)T s ]. Let the data received by the nth array element at time t be represented as y n (t), then the K data received by the nth array element can be represented as a data vector

[0021] The data vectors of array elements 1 to N are arranged in rows to obtain the received signal matrix where Y(n, k) represents the data received by the nth array element at time (k-1)T s .

[0022] In a specific embodiment, let Q = 3, N = 32, and K = 100. Three single-tone sinusoidal signal radiation sources are simulated, with a frequency of f c = 433 MHz, an angle of arrival of θ G = [85, 100, 105] degrees, and a received signal amplitude of U G = [1, 0.6, 0.85]. The received signal matrix obtained in step (1) is

[0023] (2) Use singular value decomposition to compress the received signal matrix Y, and rearrange to obtain the compressed signal vector y v .

[0024] (2a) Singular value decomposition is performed on the received signal matrix Y:

[0025] Y = U∑V T

[0026] where U represents the left singular matrix, represents the singular value matrix, and the singular values are stored in the main diagonal in reverse order, represents the right singular matrix.

[0027] (2b) Construct the singular value sampling matrix where I Q denotes a QxQ identity matrix, and 0 denotes a Qx(K-Q) all-zero matrix.

[0028] In specific embodiments,

[0029] (2c) Calculate the compressed signal matrix

[0030]

[0031] (2d) Calculate Y V Reorder by column to obtain the compressed signal (column) vector

[0032] In specific embodiments,

[0033] Step (2) extracts the singular value vector corresponding to the first Q large singular value from the received signal matrix Y by borrowing the idea of principal component analysis to obtain Y V , so that the number of signal matrix Y is reduced from K to Q, which can effectively reduce the memory occupation and calculation amount in the case of Q << K.

[0034] (3) Construct a convex optimization mathematical model with one norm constraint, and solve the convex optimization problem to obtain the solution x of the convex optimization problem j , and the solution is denoted as

[0035] (3a) Calculate the steering dictionary matrix D:

[0036] Let the number of DOA grid points be M, and divide the interval [0, 180] into M points to obtain the DOA grid Let the array element spacing vector be Calculate the steering dictionary matrix

[0037] D = exp(-1j x 2πdcos(θ))

[0038] Where 1j represents a unit imaginary number.

[0039] In specific embodiments, let M = 181, then θ = [0, 1, 2,..., 180], The received signal matrix Y in step (1) is obtained by the following formula in the simulation case:

[0040] Y = D G X G + W

[0041] Where, represents the steering matrix, and the calculation method is as follows:

[0042] D G = exp(-1j x 2πdcos(θ G ))

[0043] represents the target signal matrix, each row of which represents the signal sampling data of a radiation source. In specific embodiments, there are:

[0044]

[0045] W represents an additive complex Gaussian white noise matrix, in embodiments SNR = 0 dB (based on the maximum signal amplitude of 1 V) to simulate the thermal noise in actual equipment, while the random seed is fixed using the rng(0) code on the MATLAB platform to repeat the test;

[0046] (3b) Repeat the steering dictionary matrix D Q times to obtain the block diagonal dictionary matrix

[0047]

[0048] In specific embodiments, stored in the form of a sparse matrix, which can reduce memory consumption and computational overhead.

[0049] (3c) Solve the following convex optimization problem to obtain

[0050]

[0051] where Φ1(x) = λ|x| represents a norm penalty function, λ > 0 represents a regularization parameter, which can be calculated as follows:

[0052]

[0053] where norminv represents the inverse cumulative distribution function of the normal distribution, the above formula represents the inverse function value of 0.999 probability in the normal cumulative distribution function with a mean of 0 and a standard deviation of , where σ represents the standard deviation of the real and imaginary parts of the additive complex Gaussian white noise matrix W. In specific embodiments, σ = 0.7071 and λ = 12.3609.

[0054] The convex optimization problem in step (3c) can be solved by the max-min (MM) algorithm, with the following steps:

[0055] (3d) Set the iteration number to j, and the maximum iteration number to N J , the optimization variable of the jth iteration to x j , and let j = 1,

[0056] (3e) Construct the diagonal matrix W j-1 -1 = diag(|x j-1|. / Φ′(|x j-1 |)), and optimize variable x using this diagonal matrix. j-1 Update:

[0057]

[0058] Where diag(·) represents the diagonal matrix construction operation, . / represents point division, and Φ′(·) represents the derivative of the penalty function Φ(·). I represents the identity matrix, G -1 Let G be the inverse matrix, specifically, for the norm penalty function Φ. l1 (x)=λ|x|, there is Φ′(|x|)=λ;

[0059] (3f) Determine if j = N J If the condition is true, then the solution x to the convex optimization problem is obtained. j Let the solution be denoted as Otherwise, let j = j + 1 and execute step (3e).

[0060] In a specific embodiment, the maximum number of iterations used to solve the convex optimization problem in step (3c) is N. J =50.

[0061] (4) Construct the IMSC optimization model and iteratively calculate the optimization result x. opt The IMSC algorithm process is as follows: Figure 2 As shown:

[0062] (4a) Let the iteration variable i = 1,

[0063] (4b) will The median value is greater than the threshold T i The elements are placed into a set, resulting in a set of non-zero indices. The elements in the set are arranged in ascending order, and the number of elements is expressed as follows:

[0064] Specifically, the threshold can be set Where std(·) represents taking the vector standard deviation;

[0065] (4c) If Then the iteration ends, let Proceed to step (5); otherwise proceed to step (4d).

[0066] In the first iteration of the specific embodiment, the threshold T1 = 5.8632. In the second iteration, T2 = 5.9065. At this time there is The iteration is ended after one round. The reason why the iteration converges fast is that the initial value is selected properly. If the iteration number of the initial value problem is low, the iteration number here will increase.

[0067] (4d) Take the l column of , According to the original order of the l column in , form a sub-dictionary matrix

[0068] (4e) Solve the following semi-definite programming (SDP) problem to get a positive definite lower bound vector r opt

[0069]

[0070] S·Lr≥α min

[0071]

[0072] where the optimization variable α min represents the minimum eigenvalue of , diag(r) represents the vector r into a diagonal matrix, ≥0 represents the matrix positive definite.

[0073] In specific embodiments, the semi-definite programming problem in step (4e) is solved by the CVX toolbox of MATLAB platform. In the first round of iteration,

[0074] (4f) Calculate the penalty parameter vector where β∈[0, 1] represents a non-convex factor.

[0075] The larger the β is, the stronger the non-convexity of the penalty function is. When β = 0, the penalty function degenerates into a norm function (at this time, the penalty function is convex), and when β = 1, it can guarantee that the objective function is convex while the non-convexity of the penalty function is the strongest. In specific embodiments, β = 1.

[0076] (4g) Solve the following convex optimization problem to get the optimal sub-vector

[0077]

[0078] where represents a non-convex penalty function. Specifically, an inverse tangent non-convex penalty function is as follows:

[0079]

[0080] Note that the convex optimization problem described in step (4g) has the same format as the convex optimization problem described in step (3c), therefore in specific embodiments, the MM algorithm is also used to optimize the convex optimization problem in step (4g), with the maximum number of iterations being N J = 50. At this time, the derivative function Φ'(x) of the penalty function in step (3e) is the derivative function of a non-convex penalty function For the arctangent non-convex penalty function, its derivative function is

[0081]

[0082] Note that the calculation of the above formula should use scalar operation instead of vector operation.

[0083] (4h) Update the optimization result

[0084]

[0085] wherein represents the position sequence number of m in the set .

[0086] (4i) Let i = i + 1, and execute step (4b);

[0087] (5) Obtain the multi-target angle and amplitude estimation results by using the optimization result x opt .

[0088] (5a) Rearrange the optimization result vector x opt by column to obtain the optimization result matrix

[0089] (5b) Calculate to obtain the estimation matrix

[0090]

[0091] (5c) Average each row element of the estimation matrix to obtain the time-averaged estimation vector

[0092] In specific embodiments, the true value of the estimation vector is

[0093]

[0094] In specific embodiments, the obtained by the above steps of the present application is as shown in the attached Figure 3 .

[0095] (5d) From Select the element with the largest value Let the corresponding element position be m, one radiation source target is detected, the estimated wave arrival angle is θ(m), and the estimated wave arrival amplitude is

[0096] (5e) If Q targets have been detected, the algorithm ends, otherwise let Go to step (5d).

[0097] In the embodiments of the present application, the technical effects of the present application are described in detail in combination with simulation experiments.

[0098] 1. Simulation conditions:

[0099] Simulation 1: Multi-target direction finding and amplitude measurement simulation experiment

[0100] To test the effectiveness of the method for multi-target direction finding and amplitude measurement, two existing technologies are added as a comparison on the basis of the specific embodiments. The existing technology 1 selects the classical beam forming method (CBF), and the existing technology 2 selects the sparse singular value decomposition method based on one norm (l1-SVD). The rest of the simulation parameters are set the same as in the specific embodiments.

[0101] Simulation 2: Single target amplitude error comparison simulation experiment

[0102] To compare the deviation of the target signal amplitude estimation of the present method and the existing technology, a single target amplitude error simulation is performed. A single tone sinusoidal signal source is generated in the simulation, with a wave arrival angle of 80 degrees and a received signal amplitude of 1.5V. Let Q=1, N=8, K=1400, SNR from 0dB with 0.5dB interval to 10dB, and the MM algorithm iteration number N J =5 when the initial value is solved. The rest of the simulation parameters are set the same as in the specific embodiments. The amplitude deviation is represented as the difference between the target estimated amplitude obtained by each method and the actual target amplitude.

[0103] Simulation process in software and hardware environment:

[0104] Hardware environment: CPU is Inter Core i7-8750H, main frequency is 2.2Ghz, RAM is 24GB, frequency is 2666MHz.

[0105] Software environment: The system environment is Windows 11 Education 21H2, and MATLAB R2022b simulation software is used.

[0106] 2. Simulation content and result analysis:

[0107] Simulation 1: Multi-target direction finding and amplitude measurement simulation experiment

[0108] The results of the direction finding and amplitude measurement of multiple targets by the present method and the existing technology are shown in the accompanying Figure 3As shown, the direction finding results are consistent with the actual results, and the amplitude measurement results are shown in Table 1, wherein the average error of the amplitude measurement is the average value of the difference between the estimated amplitude of each target and the true amplitude after taking the modulus.

[0109] Table 1. Comparison of amplitude direction finding results of three technologies in a specific embodiment

[0110]

[0111] As shown in the accompanying Figure 3 Table 1, the present application can measure the direction and amplitude of multiple targets, and the average error of the amplitude measurement is significantly reduced compared with the prior art.

[0112] Simulation 2: Single target amplitude error comparison simulation experiment

[0113] Under the conditions of simulation 2, a single target amplitude error comparison simulation experiment was performed, and the amplitude deviation curves of the prior art 1, the prior art 2 and the present application for the target under different SNRs were obtained, as shown in the accompanying Figure 4 The amplitude error is the average value obtained after 100 Monte Carlo experiments at each SNR.

[0114] As shown in the accompanying Figure 4 , the absolute value of the amplitude deviation of the three methods decreases with the increase of SNR, and the absolute value of the amplitude deviation of the present application is the lowest, which indicates that the present application has better amplitude measurement accuracy. As shown in the accompanying Figure 4 , the amplitude deviation of the prior art 1 and the present application is negative, indicating that the estimated amplitude is slightly lower than the actual amplitude, while the prior art 2 is the opposite, and its estimated amplitude is higher than the actual amplitude, which provides a direction for subsequent error compensation.

[0115] The above description shows and describes one preferred embodiment of the present application, but as mentioned above, it should be understood that the present application is not limited to the form disclosed herein, and should not be considered as excluding other embodiments, but can be used in various other combinations, modifications and environments, and can be modified within the scope of the inventive concept described herein by the above-mentioned teaching or related art or knowledge. Changes and modifications made by those skilled in the art without departing from the spirit and scope of the present application shall be within the scope of the appended claims of the present application.

Claims

1. A multi-objective passive array direction-finding and amplitude-measuring method based on iterative maximum sparse convex optimization, characterized in that: The method comprises the following steps: S1. establishing a multi-target passive array sparse signal processing model to obtain a received signal matrix Y; S2. Compress the received signal matrix Y using singular value decomposition, and rearrange to obtain a compressed signal vector y v ; S3. Construct a convex optimization mathematical model with norm constraint, and solve the convex optimization problem to obtain the solution of the convex optimization problem, and record the solution as S4. Construct the IMS C optimization model, and iteratively calculate the optimization result x opt ; The step S4 comprises: S401. Let an iteration variable i = 1, Q represents the number of target radiation sources, and M represents the number of DOA grid points. S402. put elements with middle value greater than threshold T i into a set, get non-zero index set: The elements in the set are arranged in ascending order, and the number is represented as S403. Determine whether a condition is met If the condition is satisfied, the iteration ends, and the value of and goes to step S5. If the condition is not satisfied, go to step S404; S404. Take of the l column, in the l column original order to form a sub-dictionary matrix of the l column N represents the number of array elements; S405. Solve the following semi-definite programming problem to obtain a positive definite lower bound vector r opt s.t.r > a min where the optimization variables α min denotes the minimum eigenvalue of A, diag(r) denotes the vector r vectorized as a diagonal matrix, denotes that the matrix is positive definite; S406. Compute the penalty parameter vector where β ∈ [0, 1] represents a non-convex factor; The greater the β is, the stronger the non-convexity of the penalty function is. When β = 0, the penalty function degenerates into a norm function, and the penalty function is a convex function. When β = 1, the non-convexity of the penalty function is the strongest while the objective function is a convex function; S407. Solve the following convex optimization problem to obtain the optimal sub-vector wherein represents a non-convex penalty function, specifically an inverse tangent non-convex penalty function as follows: The derivative of the penalty function Φ'(x) in step S304 at this time is the derivative of the non-convex penalty function For the arctangent non-convex penalty function, the derivative is The calculation of the above formula uses scalar operation; S408. Update the optimization result wherein m represents the position number of m in the set {0, 1, 2, 3, 4, 5, 6, S409. Let i = i + 1, and return to step S402; S5. Utilize the optimization result x opt Obtain multi-target angle and amplitude estimation results.

2. The multi-objective passive array direction finding and amplitude measurement method based on iterative maximum sparse convex optimization according to claim 1, characterized in that: The step S1 comprises: Supposing that Q represents the number of target radiation sources, Q targets respectively generate Q far-field narrowband signals to be incident to a receiving array, and the receiving array is a one-dimensional uniform linear array (ULA), the number of array elements is N, N > Q, the element spacing is half the wavelength of the carrier frequency of the narrowband signal, and each array element antenna is an omnidirectional antenna; The spatial electromagnetic wave signal is received by the ULA, preprocessed by a radio frequency front end, and then subjected to digital intermediate frequency quadrature sampling, where the sampling point number is K and the sampling interval is T s The time sampling vector is represented as t = [0, T s , 2T s , …, (K-1)T s ] Let the data received by the ntharray element at time t be denoted as y n (t), then the K data received by the ntharray element is denoted as a data vector The data vectors of the 1st to Nth elements are arranged in rows to obtain a received signal matrix where Y(n, k) represents the data received by the nth element at the (k-1)T s time.

3. The multi-objective passive array direction finding and amplitude measurement method based on iterative maximum sparse convex optimization according to claim 2, characterized in that: The step S2 comprises: S201. Singular value decomposition is performed on the received signal matrix Y: Y = U∑V T wherein denotes the left-singular matrix, denotes the singular-value matrix, whose main diagonal, in reverse order, holds the singular values, denotes the right-singular matrix; S202. Constructing a singular value sampling matrix where I Q denotes a QxQ identity matrix and 0 denotes a Qx(K-Q) zero matrix. S203. Calculate the compressed signal matrix S204. Y V rearranged by column, resulting in a compressed signal column vector 4. The multi-objective passive array direction finding and amplitude measurement method based on iterative maximum sparse convex optimization according to claim 2, characterized in that: The step S3 comprises: S301. A steering dictionary matrix D is calculated: Let the number of DOA grid points be M, and the interval [0, 180] is evenly divided into M points to obtain the DOA grid Let the array element spacing vector be The steering dictionary matrix is calculated as D = exp (-1j x 2pi d cos (theta)) Wherein, 1j represents a virtual unit; S302. Repeat the steering dictionary matrix D for Q times to obtain a block diagonal dictionary matrix S303. Solve the following convex optimization problem to obtain Wherein, Phi1 (x) = lambda x represents a norm penalty function, and lambda > 0 represents a regularization parameter, and the calculation is as follows: where norminv denotes a normal inverse cumulative distribution function, the above expression denotes an inverse function value of a normal cumulative distribution function with a mean of 0 and a standard deviation of 0.999 probability, and σ denotes a standard deviation of real and imaginary parts of the additive complex Gaussian white noise. where norminv denotes a normal inverse cumulative distribution function, the above expression denotes an inverse function value of a normal cumulative distribution function with a mean of 0 and a standard deviation of 0.999 probability, and σ denotes a standard deviation of real and imaginary parts of the additive complex Gaussian white noise. S304. The convex optimization problem is solved by a maximum minimum method, and the steps are as follows: A1, let the iteration number be j, and the maximum iteration number be N J , the optimization variable of the jth iteration is x j , and let j = 1, A2, construct the diagonal matrix W j-1 -1 = diag(|x j-1 |. / Φ'( |x j-1 | ), and update the optimization variable x j-1 by this diagonal matrix where diag(·) denotes the operation of constructing a diagonal matrix,. / denotes point division, and Φ'(·) denotes the derivative of the penalty function Φ(·), I denotes the identity matrix, G -1 denotes the inverse matrix of G, for a norm penalty function Φ l1 (x) = λ |x|, with Φ'( |x| ) = λ; A3. If j = N, then go to step A4. J If yes, then the solution x of the convex optimization problem is obtained j Let this solution be denoted by Otherwise, let j = j + 1 and go to step A2.

5. The multi-objective passive array direction finding and amplitude measurement method based on iterative maximum sparse convex optimization according to claim 4, characterized in that: The step S5 comprises: S501. Per-column rearranging the optimization result vector x opt , resulting in an optimization result matrix S502. Compute the estimation matrix S503. Average each row element of the estimate matrix to obtain a time-averaged estimate vector S504. From the element with the largest value Suppose the corresponding element position is m, then a radiation source target is detected, the estimated angle of arrival is θ(m), and the estimated angle of arrival amplitude is S505. If Q targets have been detected, the step ends, otherwise let Proceeding to step S504.

Images