CN111257845B - Approximate message transfer-based non-grid target angle estimation method - Google Patents

Approximate message transfer-based non-grid target angle estimation method Download PDF

Info

Publication number
CN111257845B
CN111257845B CN202010086824.0A CN202010086824A CN111257845B CN 111257845 B CN111257845 B CN 111257845B CN 202010086824 A CN202010086824 A CN 202010086824A CN 111257845 B CN111257845 B CN 111257845B
Authority
CN
China
Prior art keywords
algorithm
value
angle
representing
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010086824.0A
Other languages
Chinese (zh)
Other versions
CN111257845A (en
Inventor
张新禹
张俊
霍凯
张双辉
刘永祥
姜卫东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202010086824.0A priority Critical patent/CN111257845B/en
Publication of CN111257845A publication Critical patent/CN111257845A/en
Application granted granted Critical
Publication of CN111257845B publication Critical patent/CN111257845B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418—Theoretical aspects

Abstract

The invention belongs to the field of array signal processing, and particularly relates to an out-of-grid target angle estimation method based on approximate message transfer, which comprises the following steps of: s1 parameter gamma to be estimatedj，j＝1，2，…，N、σ0And initializing a dictionary matrix A; s2, rapidly obtaining a signal posterior probability density function at each moment by using an AMP algorithm; s3 updating the parameter gamma to be estimated by using EM algorithmjJ is 1,2 … N, and the target number optimal solution K is obtained*Updating the noise power σ0(ii) a S4 updating dictionary matrix A by gradient descent method
Figure DDA0002382349660000011
S5 determines when the algorithm iteration converges and determines the incoming wave direction and number. The invention has the following beneficial effects: the method can automatically adjust the grid division for the targets outside the angle grid nodes at the cost of less increased calculation complexity, so that the targets are positioned on the grid nodes, the angles of the targets can be more accurately estimated, the calculation efficiency is higher, the method can be applied to a real-time multi-target high-precision angle estimation system, and the method has important engineering application value.

Description

Approximate message transfer-based non-grid target angle estimation method
Technical Field
The invention belongs to the field of array signal processing, and particularly relates to an off-grid (off-grid) target angle estimation method based on approximate message transfer.
Background
In the fields of radar, sonar, seismic remote sensing, and the like, target angle estimation (DOA) is a fundamental problem. In recent years, the mainstream idea is to introduce a sparse recovery algorithm, construct a sparse representation model of an array receiving signal according to the sparsity characteristic of the receiving signal, and complete real-time reconstruction of the signal. Compared with the traditional algorithm, the method has the remarkable advantage of improving the robustness of the signals with limited noise and sampling number and related signals.
Although solving such sparse recovery problems is difficult, accurate recovery is also possible for a strong sparse signal and well-designed dictionary matrices. Sparse Bayesian Learning (SBL) is one of the classic methods. There are two types of estimation in this approach. The first is a maximum a posteriori estimation based on sparse induced prior distribution. The second type of estimation operates in the hidden variable space with a variational representation of the hidden variable distribution that results in sparse estimates of a posteriori information beyond the mode. A series of algorithms based on the method are provided for different application scenes. These algorithms work on estimates of the hyperparameters in the Bayesian hierarchical model either by expectation-maximization iterations or by maximization iterations of evidence functions. However, each iteration of the SBL algorithm involves large-scale matrix inversion, and particularly in multi-objective solution, the calculation complexity of the SBL algorithm is very high, so that the SBL algorithm is limited to be applied to practical engineering. In order to effectively reduce the algorithm complexity, an Approximate information transfer Algorithm (AMP) is adopted in an E-Step of the SBL (see Chinese patent for details, a rapid target angle estimation method based on sparse Bayesian learning, and a patent number of ZL 2019100131299). The AMP algorithm is used for performing loop confidence coefficient propagation by using an approximate value of a central limit theorem on a bottom-layer factor graph, and effectively compensates sparse undersampling by introducing a message transfer term into an iterative threshold scheme, so that the posterior probability distribution of signals is obtained by low-complexity operation.
Consider an array of m antennas receiving signals from different directions theta1θ2… θK]TOf K targets, here [ … ]]TRepresenting the transpose of the matrix. The signal model of the SBL algorithm may be written as
y(t)＝Ax(t)+n(t) (1)
Where y (t) represents the array received signal at time t,
Figure GDA0002630746200000011
a matrix of a dictionary is represented,
Figure GDA0002630746200000012
the steering vector of the array represents a received signal vector formed by plane waves incident to the array under a far-field condition, the vector is a complex vector, the dimension of the complex vector is equal to the number of array elements of the array, and for a uniform linear array consisting of m antennas, the steering vector can be expressed as:
Figure GDA0002630746200000013
wherein d represents the distance between array elements, lambda represents the wavelength of the electromagnetic waves emitted by the array, pi is a circumferential rate constant, and theta represents any incident angle. In dictionary matrix A
Figure GDA0002630746200000014
A grid division related to a space incident angle is formed, the finer the grid division is, the higher the resolution of target angle estimation is, generally, a possible space domain interval is divided at equal intervals, and the value of N is far larger than the number m of array elements. (1) X (t) in (a) represents a received target signal vector, whose elements are in one-to-one correspondence with each column in the dictionary matrix, and when there is a target at its corresponding angle of incidence, the elements of x (t) are equal to the complex values of the signals; when the corresponding angle does not have a target, the element value of x (t) is 0. n (t) represents the additive noise of the system. Here we assume that the target signals at different times are statistically independent, an assumption that is quite common in radar signal processing. In practical engineering applications, n (t) is generally assumed to be zero-mean white Gaussian noise, i.e., white Gaussian noise
Figure GDA0002630746200000021
And is different fromThe noise between the moments is statistically independent, here
Figure GDA0002630746200000022
A complex Gaussian distribution representing the zero mean of the complex number with a variance of σ0I，σ0Representing the power of the noise, I represents the identity matrix; as used hereinafter
Figure GDA0002630746200000023
Representing a complex gaussian distribution function with a mean of α and a variance of β.
According to the signal model, the echo signals of the array can meet the distribution
Figure GDA0002630746200000024
For Bayesian derivation, the prior probability distribution of the signal vector x (t) is generally assumed to be
Figure GDA0002630746200000025
And the signals at different time instants are statistically independent, wherein
Figure GDA0002630746200000026
Is a diagonal matrix with the element gamma on the diagonaljJ is 1, and 2 … N is the variance of its corresponding element in x (t). When gamma isjWhen the angle is approximately equal to 0, the corresponding incident angle is shown
Figure GDA0002630746200000027
No target present; otherwise, the target exists. It is to be noted here that although the signals at different time instants are statistically independent, their prior probability distributions at different time instants are the same for the signal at a certain angle of incidence due to the signals from the same transmit signal source. In summary, the parameters to be estimated by the sparse bayesian algorithm include the noise power σ0And a target number K.
According to the signal model, the resolution and the angle estimation accuracy of the sparse Bayesian algorithm are determined by the dictionary matrix A reflecting the grid division density, but no matter how fine the grid is, targets outside the grid always exist, and echo signals of the targets cause mismatching of the dictionary matrix A, so that the target estimation performance of the system is reduced. Therefore, a new sparse bayesian algorithm is needed to process targets outside the grid, and the angle estimation accuracy of the algorithm is improved.
Disclosure of Invention
The invention mainly solves the technical problem that under the condition that a target is not positioned at a space angle grid node, the traditional target angle estimation algorithm based on Bayesian learning can only identify one grid node adjacent to the target, and the real angle of the target is difficult to accurately estimate.
The invention provides a grid-off target angle estimation method based on approximate message transfer (AMP), aiming at the problem that the existing Bayes learning algorithm based on approximate message transfer (AMP) (Chinese patent: a rapid target angle estimation method based on sparse Bayes learning, and patent number ZL2019100131299) can not more accurately estimate off-grid target angles.
The technical scheme adopted by the invention is as follows: an approximate message transfer based off-grid target angle estimation method comprises the following steps:
s1 parameter gamma to be estimatedj，j＝1，2，…，N、σ0And initialization of dictionary matrix A
S1.1, constructing a dictionary matrix A, A according to the required target angle estimation resolution requirement, wherein the angle corresponding to each column of the A constitutes the gridding division of the AMP algorithm on the spatial angle, and the denser the gridding division is, the higher the resolution of the angle estimation is. For the method of the invention, the equal interval division of the full angle space is adopted to initialize A, and the initial angle corresponding to the grid is recorded as
Figure GDA0002630746200000031
S1.2 in this step, the parameters needed subsequently will be initialized. The parameter to be initialized is γ ═ γ [ γ ]1γ2… γN]TAnd noise power σ0. A good initialization parameter value can greatly accelerate the convergence speed of the following algorithm and quickly obtain a correct result. Since there is no prior information about the target angle in general applications, it is initialized to γ at initialization0I.e. the signal a priori variances in the respective directions are equal. The received data from T different time points Y ═ Y (1) Y (2) … Y (T) obtained by sampling]，γ0And σ0This can be obtained from the following equation:
Figure GDA0002630746200000032
in the above formula, m is the number of array elements formed by antennas, | | … | | non-conducting phosphor2Representing the two-norm of the matrix, SNR representing the pre-estimated system signal-to-noise ratio, tr (…) representing the trace of the matrix, (…)HRepresents a conjugate transpose of the matrix;
s2 uses AMP algorithm to quickly obtain signal posterior probability density function of each time
S2.1, according to the initialization result in S1, calculates the following steps for the received data y (T), T ═ 1, and 2 … T at different times, respectively (i.e., for each T, T ═ 1, and 2 … T, steps S2.1.1-S2.1.6 are repeated until the AMP algorithm converges for the data at time T, such repeated steps are performed T times in total, for convenience of description, the time reference T is omitted in the following description of the steps, and for example, the target signal vector x (T) received at time T is abbreviated as x, and the array received signal y (T) at time T is abbreviated as y);
description of the drawings: variables occurring in Steps S2.1.1-S2.1.6
Figure GDA0002630746200000033
And
Figure GDA0002630746200000034
are all intermediate variables and have no actual physical meaning; while
Figure GDA0002630746200000035
An estimated quantity representing x estimated by the following steps; the iteration of the algorithm referred to below refers to the iteration of steps S2.1.1-S2.1.6;
S2.1.1AMP parameter initialization: for each element of x, the initial estimated parameter values are set as follows
Figure GDA0002630746200000036
Here, the first and second liquid crystal display panels are,
Figure GDA0002630746200000037
to represent
Figure GDA0002630746200000038
The (j) th element of (a),
Figure GDA0002630746200000039
to represent
Figure GDA00026307462000000310
Initial value of the jth element of (1), xjThe jth element representing the x true,
Figure GDA00026307462000000311
presentation pair
Figure GDA00026307462000000312
The initial value obtained by the estimation is obtained,
Figure GDA00026307462000000313
representing the probability density function p (x)j|γj) Expectation, where p (x)j|γj) Is shown at known gammajValue of xjIs determined by the probability density function of (a),
Figure GDA00026307462000000314
presentation pair
Figure GDA00026307462000000315
And estimating an initial value, wherein k represents the kth iteration of the algorithm, and k is 0 to represent an initialization step.
Since we generally assume the probability density function p (x)j|γj) Is a zero mean Gaussian distribution, so we can get from (5)
Figure GDA0002630746200000041
S2.1.2 linear output step: for each i-1, 2 … m, calculate
Figure GDA0002630746200000042
Figure GDA0002630746200000043
Figure GDA0002630746200000044
In the above formula
Figure GDA0002630746200000045
Representing the process of the kth iteration of the algorithm
Figure GDA0002630746200000046
Value of aijThe element representing the ith row and jth column of the dictionary matrix A, (…)iRepresents the ith element of the vector, | … | represents the modulus of the complex number,
Figure GDA0002630746200000047
representing the process of the kth iteration of the algorithm
Figure GDA0002630746200000048
The value of the one or more of,
Figure GDA0002630746200000049
representing the process of the kth iteration of the algorithm
Figure GDA00026307462000000410
The value of the one or more of,
Figure GDA00026307462000000411
representing the process of the kth iteration of the algorithm
Figure GDA00026307462000000412
The value of the one or more of,
Figure GDA00026307462000000413
representing the process of the kth iteration of the algorithm
Figure GDA00026307462000000414
The value of the one or more of,
Figure GDA00026307462000000415
representing the process of the kth iteration of the algorithm
Figure GDA00026307462000000416
The value is obtained.
S2.1.3 nonlinear output step: for each i-1, 2 … m, calculate
Figure GDA00026307462000000417
yiThe i-th element representing the received data y,
Figure GDA00026307462000000418
representing the process of the kth iteration of the algorithm
Figure GDA00026307462000000419
The value of the one or more of,
Figure GDA00026307462000000420
indicating that the k-th iteration of the algorithm is further performedNovel
Figure GDA00026307462000000421
Value of (1), function of
Figure GDA00026307462000000422
S2.1.4 Linear input step: for each j ═ 1,2 … N, calculations were made
Figure GDA00026307462000000423
Figure GDA00026307462000000424
Representing the process of the kth iteration of the algorithm
Figure GDA00026307462000000425
The value of the one or more of,
Figure GDA00026307462000000426
representing the process of the kth iteration of the algorithm
Figure GDA00026307462000000427
Value, here (…)-1Representation matrix inversion (…)*Representing the conjugate of a complex number.
S2.1.5 nonlinear input step: for each j ═ 1,2 … N, calculations were made
Figure GDA00026307462000000428
Here, the
Figure GDA00026307462000000429
Representing the (k + 1) th iteration
Figure GDA00026307462000000430
The value of the one or more of,
Figure GDA00026307462000000431
representing the (k + 1) th iteration
Figure GDA00026307462000000432
Value, function of above
Figure GDA0002630746200000051
S2.1.6 determining whether the AMP algorithm converges: computing
Figure GDA0002630746200000052
Wherein … does not calculation1A 1-norm of the matrix is represented,
Figure GDA0002630746200000053
representing the (k + 1) th iteration
Figure GDA0002630746200000054
The values, in the same way,
Figure GDA0002630746200000055
representing the kth iteration
Figure GDA0002630746200000056
The value is obtained. If the value is greater than a set threshold1Then go back to S2.1.2 for reiteration; otherwise, the loop is skipped to S2.2 to obtain p (x)jY) of the results. Threshold (THD)1Depending on factors such as the signal-to-noise ratio of the system, it needs to be adjusted according to actual conditions. Normally, the threshold is1Is between 0.1 and 0.001.
S2.2 through the steps, the posterior probability density function p (x) of the signal at the time t can be obtainedjResults of | y) are as follows
Figure GDA0002630746200000057
p(xj) Denotes xjA probability density function of; in the above formula
Figure GDA0002630746200000058
Is obtained by calculation according to the sample value y (t) at the time t in the step S2.1, and is obtained in the last iteration process after the AMP algorithm is judged to be converged
Figure GDA0002630746200000059
A value of (d); here gamma isjIs obtained from either S1 or S3, and in the first EM algorithm iteration, γ is obtainedjIs determined by the initial value in S1, and in other cases, γjIs calculated by the gamma calculated at S3 in the last EM algorithm cyclejAnd (4) determining.
S3 updating the parameter gamma to be estimated by using EM algorithmjJ is 1,2 … N, and the target number optimal solution K is obtained*Updating the noise power σ0
The posterior probability density function p (x) of the signal has been obtained in S2jY), according to the EM algorithm, this step updates the value of the parameter to be estimated one by one using the following expression
Figure GDA00026307462000000510
In the above formula, q ═ γ1… γNσ0]T，X＝[x(1) x(2) … x(T)]，N＝[n(1) n(2) … n(T)]，＜…|Y；qiDenotes that the known reception data Y ═ Y (1) Y (2) … Y (t)]And given parameter value qiIs averaged under the condition of (1), q in the above expressioniRepresents the value of q during the ith iteration of the algorithm, and qi+1Representing the q value during the (i + 1) th iteration of the algorithm. The method comprises the following steps:
s3.1 updating Gammaj,j＝1,2…N：
Because the signals at different time are independent in statistics, the parameters to be estimated are the same. Due to gammajIs updated only with
Figure GDA00026307462000000511
In this regard, the probability density function when expected can therefore become p (x)j(t)|y(t)；qi) I.e. known received data y (t)) And given parameter value qiX under the condition of (1)j(t) probability density function. The pair gamma can be obtained by formula derivationjThe updated expression of j 1,2 … N is
Figure GDA0002630746200000061
In the above expression
Figure GDA0002630746200000062
Represents gamma in the ith iteration of the algorithmj，
Figure GDA0002630746200000063
Represents gamma in the (i + 1) th iteration of the algorithmj。
Further on gammajThe partial derivative can be obtained
Figure GDA0002630746200000064
As can be seen from the above formula, γjThe update of j 1,2 … N does not need to go through matrix operation, but is simple scalar operation, so that a large amount of operation time can be saved.
S3.2 estimating target quantity optimal solution K*: different from the Chinese invention patent 'a rapid target angle estimation method based on sparse Bayesian learning', patent number ZL2019100131299, the invention needs to determine the number of targets before estimating the noise power, so as to conveniently update the noise power sigma at S3.30The intermediate result calculated in the step can be utilized to reduce the consumption of calculation resources.
This step is described in references 1, z. -m.liu, z. -t.huang, y. -y.zhou, An effective massive method for direction-of-arrival estimation of video streaming, IEEE trans. wireless Communications 11(10) (2012) 3607-3617; 2. stoica, Y.Selen, Model-order selection a review of information criteria rules, IEEESignal Process. Mag.21(4) (2004) 36-47; 3. austin, R.L.Moses, J.N.Ash, E.Ertin, on the relationship between search retrieval and parameter estimation with model order selection, IEEE J.Sel.topics Signal Process.4(3) (2010) 560-.
As can be seen from the above references, using subspace analysis, the optimal solution for the target number is
Figure GDA0002630746200000065
Wherein I is a matrix of units and I is a matrix of units,
Figure GDA0002630746200000066
Figure GDA0002630746200000067
where K is an estimate of the target number, K*Is the optimal solution of the K, and the K is the optimal solution,
Figure GDA0002630746200000068
corresponding assumed target angle in A
Figure GDA0002630746200000069
Is given a target angle
Figure GDA00026307462000000610
Refers to gammajJ is the angle corresponding to the largest K value of 1,2 … N.
S3.3 updating the noise Power σ0：
This step is described in references p.gerstoft, c.f. mecklenbranker, a.xenaki, s.nannuru, multisenshot spot basic learning for doa, IEEE Signal process. letters 23(10) (2016) 1469-.
From the above references, one can obtain
Figure GDA00026307462000000611
Wherein
Figure GDA0002630746200000071
Here, the
Figure GDA0002630746200000072
Is that
Figure GDA0002630746200000073
The covariance matrix of each column in (1), which can be obtained from the above equation
Figure GDA0002630746200000074
According to the mapping matrix P defined by S3.2, and P ═ PH＝P2Is obtained by
PRPH-σ0PPH＝ΣY-σ0I (21)
Combining tr (R-sigma)Y) 0 can deduce σ0The update formula of (2) is:
Figure GDA0002630746200000075
comparing equation (17), it can be seen that solving for K*And σ0Can be carried out simultaneously, and saves computing resources.
S4 updating dictionary matrix A by gradient descent method
Figure GDA0002630746200000076
Through the above steps, the number of targets and the grid nodes where the targets are located (i.e. exceeding the threshold) have been roughly estimated5Gamma of (2)jThe corresponding grid node, that is, the angle estimation value of the target), but since the real position of the target may be between two grid nodes, the grid node is adjusted by adopting a gradient descent algorithm in the step, so that the grid node is closer to the real position of the target, thereby improving the target angle estimation accuracy. Optimal solution K for target quantity calculated in S3.2*The target corresponding angle can be considered as γjJ is the first K in 1,2 … N*Maximum gammajTo what is providedCorresponding angle grid node
Figure GDA0002630746200000077
The angle grid nodes roughly corresponding to the targets are adjusted by using a gradient descent method, and it is noted that only the angle grid nodes currently corresponding to the targets need to be updated here
Figure GDA0002630746200000078
Without having to update the entire angle grid, to avoid high computational complexity (i.e., the angles corresponding to each target determined in S3.2)
Figure GDA0002630746200000079
And performing iterative calculation in the step until the algorithm is converged. The S4 steps collectively need to be performed K*Second). The update rule can be derived from
Figure GDA00026307462000000710
Wherein
Figure GDA00026307462000000711
A is a dictionary matrix with angle variables. From the above formula, the correlation
Figure GDA00026307462000000712
Gradient of (2)
Figure GDA00026307462000000713
Here, the
Figure GDA00026307462000000714
Representing the real part of the matrix, X is the estimate for X found by step S2,
Figure GDA00026307462000000715
is in A corresponds to a division
Figure GDA00026307462000000716
Outside angle
Figure GDA00026307462000000717
A matrix of a plurality of columns of (a),
Figure GDA00026307462000000718
corresponding to division in representation X
Figure GDA00026307462000000719
Outside angle
Figure GDA00026307462000000720
A matrix of rows of (a) is formed,
Figure GDA00026307462000000721
corresponding to angle in representation X
Figure GDA00026307462000000722
The matrix of rows of (a) is,
Figure GDA00026307462000000723
is the column vector in A, d (-) is the derivative of a (-).
Pass through (24) type pair
Figure GDA00026307462000000724
Is optimized by
Figure GDA00026307462000000725
Wherein α is the step size of the angle update, and the step size depends on the accuracy requirement of the angle estimation in practical applicationlThe value of the i-th iteration within this step is indicated. If the input of function sign (-) is a positive number then it equals 1, otherwise it equals-1.
After completing the angle update once, it needs to judge whether to converge. Considering that the initial dictionary grid may already meet the accuracy requirement, the number of iterations in this step is not too large, and when the condition is met
Figure GDA0002630746200000081
Or l is not less than4When it is, consider to
Figure GDA0002630746200000082
Has met the requirements, exits the loop,3、4is the decision threshold. Wherein the content of the first and second substances,3the method can be determined according to the precision requirement of the system on angle estimation, and the value is generally 0.1 degree;4the real-time performance of the algorithm is lower when the value is larger, which is determined according to the real-time performance requirement of the system, and the value can be 1000 under the common condition.
S5 judges when the iterative process of the algorithm converges, and determines the direction and quantity of the incoming wave
Finish γ ═ γ1γ2... γN]TAnd σ0After updating, convergence is judged by the following expression:
Figure GDA0002630746200000083
2the threshold can be set according to the system practice for determining the threshold, and is usually between 0.1 and 0.001. If the above expression is not satisfied, returning to S2 to continue the iterative operation; if the above formula is true, the loop can be exited, and the number of incoming waves is the last calculated K*In the direction of [ gamma ] - [ gamma ]1γ2… γN]TMiddle front K*Maximum gammajCorresponding direction
Figure GDA0002630746200000084
Due to the pair in S4
Figure GDA0002630746200000085
And more accurate calculation is performed, so that the angle estimation precision of the algorithm on the off-grid target is greatly improved.
The invention has the following beneficial effects: the method can automatically adjust the grid division for the targets outside the angle grid nodes at the cost of less increased calculation complexity, so that the targets are positioned on the grid nodes, the angles of the targets can be more accurately estimated, the calculation efficiency is higher, the method can be applied to a real-time multi-target high-precision angle estimation system, and the method has important engineering application value.
Drawings
FIG. 1 is a process flow diagram;
FIG. 2 is a spatial power spectrum of the method of the present invention and the conventional method when the number of array elements is 16;
FIG. 3 is a comparison of the performance of the method of the present invention with that of the conventional method with 16 array elements as a function of the signal-to-noise ratio;
FIG. 4 is a comparison of the performance of the method of the present invention with that of the conventional method with the number of samples collected when the number of array elements is 16;
(a) the target angle is 5 degrees and 30 degrees, (b) the target angle is-5 degrees and 15 degrees;
FIG. 5 is a graph comparing the operation time of the method of the present invention with that of the conventional method with the number of samples when the number of array elements is 16;
FIG. 6 is a graph showing the comparison of the operation time of the method of the present invention and the conventional method with the change of the number of samples when the number of array elements is 10;
FIG. 7 is a graph showing the comparison of the operation time of the method of the present invention and the conventional method with the number of samples when the number of array elements is 16 and the angle interval is 0.5 °.
Detailed Description
The invention is further illustrated with reference to the accompanying drawings:
FIG. 1 is a general process flow of the present invention.
The invention discloses a Bayesian learning rapid target angle estimation algorithm based on generalized approximate message transmission, which comprises the following steps of:
s1 parameter gamma to be estimatedj，j＝1，2，…，N、σ0And initializing a dictionary matrix A;
s2, rapidly obtaining a signal posterior probability density function at each moment by using an AMP algorithm;
s3 updating the parameter gamma to be estimated by using EM algorithmjJ is 1,2 … N, and the target number optimal solution K is obtained*Updating the noise power σ0；
S4 updating dictionary moment by gradient descent methodIn array A
Figure GDA0002630746200000091
S5 determines when the algorithm iteration converges and determines the incoming wave direction and number.
FIG. 2 is a spatial power spectrum of the method of the present invention and classical LASSO, RVM and Atomic Norm algorithms. The simulation is based on a uniform linear array with array elements of 16 and half-wavelength intervals. Consider two incoherent targets emitting signals incident on the array at-5 ° and 15.5 ° positions, respectively, with the signal-to-noise ratios of the target signals both being-10 dB, for a total of 50 collected samples received by the array. The four algorithms all adopt the same space grid division and dictionary matrix, namely, the division of 1-degree interval is carried out on the angle space ranging from-45 degrees to 45 degrees. As can be seen from the figure, the space universality of the four algorithms is the peak value at the target position, the LASSO algorithm has the narrowest peak value, the Atomic Norm algorithm is next to the narrow peak value, the RVM algorithm has the widest peak value, and the peak value width of the algorithm is in the middle of the Atomic Norm algorithm and the RVM algorithm. Therefore, an accurate angle estimation result can be obtained by the pole detection method.
FIG. 3 shows the variation of estimation accuracy with SNR for the Atomic Norm algorithm, the present algorithm, and the present algorithm without the step S4. The estimation accuracy is represented by the root mean square error of the angle estimation value of 50 times of simulation, and the expression is
Figure GDA0002630746200000092
Here, the
Figure GDA0002630746200000093
The angle estimate from the i-th simulation is shown. As can be seen, the algorithm of the present invention performed better than the Atomic Norm algorithm and the algorithm without the S4 step in the tests of the two sets of targets (angles of-5 °, 15.5 ° and 5 °, 30.2 °, respectively), especially in the case of low signal-to-noise ratio.
Fig. 4 shows the variation of the estimation accuracy with the number of samples collected in the four methods, where the snr is-10 dB, (a) the target angle is 5 ° and 30 °, (b) the target angle is-5 ° and 15 °. As can be seen, the RMSE exhibited a decreasing trend for all four methods as the number of samples increased. The RVM algorithm has the best estimation precision under the condition of small samples, but as the RMSE decreases slowly with the increase of the number of samples, the method disclosed by the invention has similar performance to that of the LASSO algorithm.
FIG. 5 is a graph showing comparison of the operation time with the number of samples. The calculation complexity of each algorithm is measured through the operation time, at the moment, the array element number is 16, the SNR is-5 dB, the grid interval is 1 degree, and the two target incidence angles are respectively 5 degrees and 30.2 degrees. As can be seen, the computation time increases as the number of samples increases, but the computation time of the method of the present invention is orders of magnitude longer than that of LASSO, RVM and Atomic Norm algorithms. Comparing fig. 3, it can be seen that the addition of the gradient descent algorithm only slightly increases the operation time, but the DOA estimation accuracy is significantly increased.
Fig. 6 is a comparison graph of the number of array elements 10 in addition to fig. 5, showing the change in the operation time with the number of samples. Comparing with fig. 5, it can be seen that the calculation time of each algorithm is obviously reduced as the number of array elements is reduced, and the algorithm of the present invention is superior to other algorithms.
Fig. 7 is a comparison graph of the operation time obtained by changing the angle division interval to 0.5 ° on the basis of fig. 5 as a function of the number of samples. As can be seen from comparison of FIG. 5, the grid precision has a significant influence on the operation time, the operation time of each algorithm is greatly increased, but the operation time of the algorithm of the invention is still significantly lower than that of other algorithms.
The simulation-based experimental result shows that the algorithm has strong robustness on noise, is still effective on small sample data, has obviously higher computational efficiency and angle estimation precision than the traditional method, and meets the requirement of real-time target angle estimation. The method can realize the accurate estimation of the incident angle of multiple targets under the condition that the quality of radar echo data is limited, and particularly provides technical support for missile defense and space target identification in space target monitoring under the condition of strong confrontation, and has high engineering application value.

Claims (6)

1. An estimation method of an out-of-grid target angle based on approximate message passing, the method comprising the steps of:
s1 waiting for estimationParameter gammaj，j＝1，2，…，N、σ0And initialization of dictionary matrix A
S1.1, constructing a dictionary matrix A, A according to the required target angle estimation resolution requirement, wherein the angle corresponding to each column of the A, A forms the gridding division of the AMP algorithm to the space angle, the denser the gridding division is, the higher the resolution of the angle estimation is, and the initial angle corresponding to the grid is recorded as
Figure FDA0002630746190000011
S1.2 vs. γ ═ γ1γ2…γN]TAnd noise power σ0Carrying out initialization; initialized to gamma at initialization0I.e. the signal prior variances in the respective directions are equal, and the sampled received data Y from T different times is [ Y (1) Y (2) … Y (T)]，γ0And σ0This can be obtained from the following equation:
Figure FDA0002630746190000012
in the above formula, m is the number of array elements formed by antennas, | | … | | non-conducting phosphor2Representing the two-norm of the matrix, SNR representing the pre-estimated system signal-to-noise ratio, tr (…) representing the trace of the matrix, (…)HRepresents a conjugate transpose of the matrix;
s2 uses AMP algorithm to quickly obtain signal posterior probability density function of each time
S2.1, according to the initialization result in S1, performs calculation for the received data y (T), T ═ 1, and 2 … T at different times, respectively, by repeating steps S2.1.1 to S2.1.6 for each T, T ═ 1, and 2 … T until the AMP algorithm converges for the data at time T, and such repetition steps need to be performed T times in total;
S2.1.1AMP parameter initialization: for each element of x, the initial estimated parameter values are set as follows
Figure FDA0002630746190000013
Here, the first and second liquid crystal display panels are,
Figure FDA0002630746190000014
to represent
Figure FDA0002630746190000015
The (j) th element of (a),
Figure FDA0002630746190000016
to represent
Figure FDA0002630746190000017
Initial value of the jth element of (1), xjThe jth element representing the x true,
Figure FDA0002630746190000018
presentation pair
Figure FDA0002630746190000019
The initial value obtained by the estimation is obtained,
Figure FDA00026307461900000110
representing the probability density function p (x)j|γj) To expect, p (x)j|γj) Is shown at known gammajValue of xjIs determined by the probability density function of (a),
Figure FDA00026307461900000111
presentation pair
Figure FDA00026307461900000112
Estimating an initial value, wherein k represents the kth iteration of the algorithm, and k is 0 and represents an initialization step;
suppose a probability density function p (x)j|γj) Is a zero mean Gaussian distribution, so we can get from (2)
Figure FDA00026307461900000113
S2.1.2 linear output step: for each i-1, 2 … m, calculate
Figure FDA00026307461900000114
In the above formula
Figure FDA0002630746190000021
Representing the process of the kth iteration of the algorithm
Figure FDA0002630746190000022
Value of aijThe element representing the ith row and jth column of the dictionary matrix A, (…)iRepresents the ith element of the vector, | … | represents the modulus of the complex number,
Figure FDA0002630746190000023
representing the process of the kth iteration of the algorithm
Figure FDA0002630746190000024
The value of the one or more of,
Figure FDA0002630746190000025
representing the process of the kth iteration of the algorithm
Figure FDA0002630746190000026
The value of the one or more of,
Figure FDA0002630746190000027
representing the process of the kth iteration of the algorithm
Figure FDA0002630746190000028
The value of the one or more of,
Figure FDA0002630746190000029
representing the process of the kth iteration of the algorithm
Figure FDA00026307461900000210
The value of the one or more of,
Figure FDA00026307461900000211
representing the process of the kth iteration of the algorithm
Figure FDA00026307461900000212
A value;
s2.1.3 nonlinear output step: for each i-1, 2 … m, calculate
Figure FDA00026307461900000213
yiThe i-th element representing the received data y,
Figure FDA00026307461900000214
representing the process of the kth iteration of the algorithm
Figure FDA00026307461900000215
The value of the one or more of,
Figure FDA00026307461900000216
indicating updating during the kth iteration of the algorithm
Figure FDA00026307461900000217
Value of (1), function of
Figure FDA00026307461900000218
S2.1.4 Linear input step: for each j ═ 1,2 … N, calculations were made
Figure FDA00026307461900000219
Figure FDA00026307461900000220
Representing the process of the kth iteration of the algorithm
Figure FDA00026307461900000221
The value of the one or more of,
Figure FDA00026307461900000222
representing the process of the kth iteration of the algorithm
Figure FDA00026307461900000223
Value, here (…)-1Representation matrix inversion (…)*Represents the conjugate of a complex number;
s2.1.5 nonlinear input step: for each j ═ 1,2 … N, calculations were made
Figure FDA00026307461900000224
Here, the
Figure FDA00026307461900000225
Representing the (k + 1) th iteration
Figure FDA00026307461900000226
The value of the one or more of,
Figure FDA00026307461900000227
representing the (k + 1) th iteration
Figure FDA00026307461900000228
Value, function of above
Figure FDA00026307461900000229
S2.1.6 determining whether the AMP algorithm converges: computing
Figure FDA00026307461900000230
Wherein … does not calculation1A 1-norm of the matrix is represented,
Figure FDA00026307461900000231
representing the (k + 1) th iteration
Figure FDA00026307461900000232
The values, in the same way,
Figure FDA00026307461900000233
representing the kth iteration
Figure FDA00026307461900000234
A value; if the value is greater than a set threshold1Then go back to S2.1.2 for reiteration; otherwise, the loop is skipped to S2.2 to obtain p (x)jResults of | y); threshold (THD)1The signal-to-noise ratio of the system is required to be adjusted according to actual conditions;
s2.2 through the steps, the posterior probability density function p (x) of the signal at the time t can be obtainedjResults of | y) are as follows
Figure FDA0002630746190000031
p(xj) Denotes xjA probability density function of; in the above formula
Figure FDA0002630746190000032
Is obtained by calculation according to the sample value y (t) at the time t in the step S2.1, and is obtained in the last iteration process after the AMP algorithm is judged to be converged
Figure FDA0002630746190000033
A value of (d); here gamma isjIs obtained from either S1 or S3, and in the first EM algorithm iteration, γ is obtainedjIs determined by the initial value in S1, and in other cases, γjIs calculated by S3 in the last EM algorithm loopγjDetermining;
s3 updating the parameter gamma to be estimated by using EM algorithmjJ is 1,2 … N, and the target number optimal solution K is obtained*Updating the noise power σ0
The posterior probability density function p (x) of the signal has been obtained in S2jY), according to the EM algorithm, this step updates the value of the parameter to be estimated one by one using the following expression
Figure FDA0002630746190000034
In the above formula, q ═ γ1…γNσ0]T，X＝[x(1) x(2)…x(T)]，N＝[n(1) n(2)…n(T)]，＜…|Y；qiDenotes that the known reception data Y ═ Y (1) Y (2) … Y (t)]And given parameter value qiIs averaged under the condition of (1), q in the above expressioniRepresents the value of q during the ith iteration of the algorithm, and qi+1Representing the q value in the (i + 1) th iteration of the algorithm; the method comprises the following steps:
s3.1 updating Gammaj,j＝1,2…N：
The pair gamma can be obtained by formula derivationjThe updated expression of j 1,2 … N is
Figure FDA0002630746190000035
In the above expression
Figure FDA0002630746190000036
Represents gamma in the ith iteration of the algorithmj，
Figure FDA0002630746190000037
Represents gamma in the (i + 1) th iteration of the algorithmj；
Further on gammajThe partial derivative can be obtained
Figure FDA0002630746190000038
S3.2 estimating target quantity optimal solution K*：
By using subspace analysis, the optimal solution of the target number is
Figure FDA0002630746190000039
Wherein I is a matrix of units and I is a matrix of units,
Figure FDA0002630746190000041
Figure FDA0002630746190000042
k is an estimate of the target quantity, K*Is the optimal solution of the K, and the K is the optimal solution,
Figure FDA0002630746190000043
corresponding assumed target angle in A
Figure FDA0002630746190000044
Is given a target angle
Figure FDA0002630746190000045
Refers to gammajJ is the angle corresponding to the largest K value in 1,2 … N;
s3.3 updating the noise Power σ0：
According to the formula:
Figure FDA0002630746190000046
wherein
Figure FDA0002630746190000047
Here, the
Figure FDA0002630746190000048
Is that
Figure FDA0002630746190000049
The covariance matrix of each column in (1), which can be obtained from the above equation
Figure FDA00026307461900000410
According to the mapping matrix P defined by S3.2, and P ═ PH＝P2Is obtained by
PRPH-σ0PPH＝ΣY-σ0I (18)
Combining tr (R-sigma)Y) 0 can deduce σ0The update formula of (2) is:
Figure FDA00026307461900000411
s4 updating dictionary matrix A by gradient descent method
Figure FDA00026307461900000412
Optimal solution K for target quantity calculated in S3.2*The target corresponding angle can be considered as γjJ is the first K in 1,2 … N*Maximum gammajCorresponding angle grid node
Figure FDA00026307461900000413
The angle grid nodes roughly corresponding to the targets are adjusted by using a gradient descent method, and it is noted that only the angle grid nodes currently corresponding to the targets need to be updated here
Figure FDA00026307461900000414
Without having to update the entire angle grid, to avoid high computational complexity, i.e. the angles corresponding to each target determined in S3.2
Figure FDA00026307461900000415
The iterative calculation of the step is carried out until the algorithm is converged, and the step S4 needs to execute K*Secondly; the update rule can be derived from
Figure FDA00026307461900000416
Wherein
Figure FDA00026307461900000417
A is a dictionary matrix with angle variables; from the above formula, the correlation
Figure FDA00026307461900000418
Gradient of (2)
Figure FDA00026307461900000419
Here, the
Figure FDA00026307461900000420
Representing the real part of the matrix, X is the estimate for X found by step S2,
Figure FDA00026307461900000421
is in A corresponds to a division
Figure FDA00026307461900000422
Outside angle
Figure FDA00026307461900000423
A matrix of a plurality of columns of (a),
Figure FDA00026307461900000424
corresponding to division in representation X
Figure FDA00026307461900000425
Outside angle
Figure FDA00026307461900000426
A matrix of rows of (a) is formed,
Figure FDA00026307461900000427
corresponding to angle in representation X
Figure FDA00026307461900000428
The matrix of rows of (a) is,
Figure FDA00026307461900000429
is the column vector in A, d (-) is the derivative of a (-);
pass through (21) type pair
Figure FDA00026307461900000430
Is optimized by
Figure FDA0002630746190000051
α, the step size of the angle update depends on the accuracy requirement of the angle estimation in practical application, (. degree)lThe value representing the ith iteration of this step, which is equal to 1 if the input of the function sign (·) is positive, and-1 otherwise;
after one-time angle updating is completed, whether convergence exists needs to be judged; considering that the initial dictionary grid may already meet the accuracy requirement, the number of iterations in this step is not too large, and when the condition is met
Figure FDA0002630746190000052
Or l is not less than4When it is, consider to
Figure FDA0002630746190000053
Has met the requirements, exits the loop,3、4is a decision threshold; wherein the content of the first and second substances,3the method can be determined according to the precision requirement of the system on angle estimation;4according toThe real-time requirement of the system is determined, and the larger the value is, the lower the real-time performance of the algorithm is;
s5 judges when the iterative process of the algorithm converges, and determines the direction and quantity of the incoming wave
Finish γ ═ γ1γ2...γN]TAnd σ0After updating, convergence is judged by the following expression:
Figure FDA0002630746190000054
2the threshold can be set according to the system practice for judging the threshold; if the above expression is not satisfied, returning to S2 to continue the iterative operation; if the above formula is true, the loop can be exited, and the number of incoming waves is the last calculated K*In the direction of [ gamma ] - [ gamma ]1γ2…γN]TMiddle front K*Maximum gammajCorresponding direction
Figure FDA0002630746190000055
2. The approximate messaging based off-grid target angle estimation method of claim 1, wherein: in S1.1, the space angle is divided into grids by adopting equal-interval division of a full-angle space.
3. The approximate messaging based off-grid target angle estimation method according to claim 1 or 2, wherein: s2.1.6, threshold1Is between 0.1 and 0.001.
4. The approximate messaging based off-grid target angle estimation method according to claim 1 or 2, wherein: in S4, a decision threshold is determined3The value is 0.1 deg.
5. The proximity-based messaging according to claim 1 or 2The grid target angle estimation method is characterized by comprising the following steps: in S4, a decision threshold is determined4The value is 1000.
6. The approximate messaging based off-grid target angle estimation method according to claim 1 or 2, wherein: in S5, a threshold is determined2The value is between 0.1 and 0.001.
CN202010086824.0A 2020-02-11 2020-02-11 Approximate message transfer-based non-grid target angle estimation method Active CN111257845B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010086824.0A CN111257845B (en) 2020-02-11 2020-02-11 Approximate message transfer-based non-grid target angle estimation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010086824.0A CN111257845B (en) 2020-02-11 2020-02-11 Approximate message transfer-based non-grid target angle estimation method

Publications (2)

Publication Number Publication Date
CN111257845A CN111257845A (en) 2020-06-09
CN111257845B true CN111257845B (en) 2020-09-22

Family

ID=70945657

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010086824.0A Active CN111257845B (en) 2020-02-11 2020-02-11 Approximate message transfer-based non-grid target angle estimation method

Country Status (1)

Country Link
CN (1) CN111257845B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112346123B (en) * 2020-11-06 2023-02-10 中国地震灾害防御中心 VIA (visual analysis of seismic data) double-parameter analysis method
CN113281702B (en) * 2021-04-30 2024-02-09 中国人民解放军战略支援部队信息工程大学 Method for directly positioning beyond-view-range target by cooperating short-wave multi-station angle with satellite time frequency
CN116879862B (en) * 2023-09-08 2023-12-01 西安电子科技大学 Single snapshot sparse array space angle super-resolution method based on hierarchical sparse iteration

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104375133A (en) * 2014-11-11 2015-02-25 西北大学 Estimation method for space two-dimensional DOA
KR101783777B1 (en) * 2016-05-17 2017-10-10 국방과학연구소 Method for estimating direction of arrival using covariance fitting
CN109116293A (en) * 2018-08-22 2019-01-01 上海师范大学 A kind of Wave arrival direction estimating method based on sparse Bayesian out of place
CN109490819A (en) * 2018-11-16 2019-03-19 南京邮电大学 A kind of Wave arrival direction estimating method out of place based on management loading
CN109658467A (en) * 2018-12-12 2019-04-19 浙江工业大学 A kind of endoscopic images sensing reconstructing method based on multiword allusion quotation modified compressed sensing framework
CN109752710A (en) * 2019-01-07 2019-05-14 中国人民解放军国防科技大学 Rapid target angle estimation method based on sparse Bayesian learning
CN110208735A (en) * 2019-06-12 2019-09-06 西北工业大学 A kind of DOA Estimation in Coherent Signal method based on management loading

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104375133A (en) * 2014-11-11 2015-02-25 西北大学 Estimation method for space two-dimensional DOA
KR101783777B1 (en) * 2016-05-17 2017-10-10 국방과학연구소 Method for estimating direction of arrival using covariance fitting
CN109116293A (en) * 2018-08-22 2019-01-01 上海师范大学 A kind of Wave arrival direction estimating method based on sparse Bayesian out of place
CN109490819A (en) * 2018-11-16 2019-03-19 南京邮电大学 A kind of Wave arrival direction estimating method out of place based on management loading
CN109658467A (en) * 2018-12-12 2019-04-19 浙江工业大学 A kind of endoscopic images sensing reconstructing method based on multiword allusion quotation modified compressed sensing framework
CN109752710A (en) * 2019-01-07 2019-05-14 中国人民解放军国防科技大学 Rapid target angle estimation method based on sparse Bayesian learning
CN110208735A (en) * 2019-06-12 2019-09-06 西北工业大学 A kind of DOA Estimation in Coherent Signal method based on management loading

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
An effeicient maximum likelihood method for direction of arrival estimation via sparse bayesian learning;Z.M.Liu等;《IEEE Trans. Wireless Communications》;20121231;第11卷(第10期);全文 *
快速块稀疏贝叶斯学习算法的理论与应用;刘本源;《中国博士学位论文全文数据库 信息科技辑》;20170215;全文 *

Also Published As

Publication number Publication date
CN111257845A (en) 2020-06-09

Similar Documents

Publication Publication Date Title
CN111257845B (en) Approximate message transfer-based non-grid target angle estimation method
CN106646344B (en) A kind of Wave arrival direction estimating method using relatively prime battle array
WO2018094565A1 (en) Method and device for beamforming under pulse noise
CN109490819B (en) Sparse Bayesian learning-based method for estimating direction of arrival of wave in a lattice
CN108872926B (en) Amplitude-phase error correction and DOA estimation method based on convex optimization
CN107450045B (en) DOA estimation method based on FOCUSS secondary weighting algorithm
CN109239649B (en) Novel co-prime array DOA estimation method under array error condition
CN111337893A (en) Off-grid DOA estimation method based on real-value sparse Bayesian learning
CN111767791A (en) Arrival angle estimation method based on anti-regularization deep neural network
CN109239646B (en) Two-dimensional dynamic direction finding method for continuous quantum water evaporation in impact noise environment
CN111337873A (en) DOA estimation method based on sparse array
Xu et al. Spatial information theory of sensor array and its application in performance evaluation
CN115236584A (en) Meter-wave radar low elevation angle estimation method based on deep learning
CN109783960B (en) Direction-of-arrival estimation method based on grid part refinement
CN114720938A (en) Large-scale antenna array single-bit sampling DOA estimation method based on depth expansion
CN113759303A (en) Non-grid DOA (angle of arrival) estimation method based on particle swarm optimization
Yang et al. A correlation-aware sparse Bayesian perspective for DOA estimation with off-grid sources
CN112731273B (en) Low-complexity signal direction-of-arrival estimation method based on sparse Bayesian
CN115015832A (en) Large-scale array amplitude-phase error and target direction joint estimation method under non-uniform noise
CN110471026B (en) Phase-enhanced meter-wave radar target low elevation DOA estimation method
CN105303009A (en) Super-resolution spectrum estimation method based on compressed sensing and regular MFOCUSS
CN109298384B (en) Non-uniform linear array direction of arrival angle estimation method based on variational Bayes inference
CN109100679B (en) Near-field sound source parameter estimation method based on multi-output support vector regression machine
CN117092585B (en) Single-bit quantized DoA estimation method, system and intelligent terminal
CN114280533B (en) Sparse Bayesian DOA estimation method based on l0 norm constraint

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant