Low-altitude Multipath Target Super-Resolution Elevation Estimation Method Based on Block Sparse Modeling

Abstract

The present invention discloses a method for super-resolution elevation estimation of low-altitude multipath targets based on block sparse modeling, including: establishing a rectangular coordinate system in the plane where the echo propagation path is located based on the positioning and orientation information of the radar and the azimuth measurement value of the target by the radar; extracting the terrain curve within a certain distance range in the direction where the target is located from the geographic information system; performing grid division in the pitch dimension for the airspace of interest, analyzing the validity of the reflection paths corresponding to each pitch angle in combination with the distance measurement value and the terrain curve, and calculating the incident pitch angle of the reflected wave; constructing a structured redundant dictionary from the direct echo steering vectors and the reflected wave steering vectors corresponding to all pitch angles of interest; establishing a structured sparse representation model for the received data in the radar element domain / subarray domain / beam domain; using a first-order fast algorithm to solve the convex combination optimization problem and recover the block sparse signal matrix; obtaining the elevation estimation value of the target based on the recovered block sparse signal, which can effectively improve the accuracy of the model.

Description

Technical Field

[0001] The present invention relates to the technical field of radar low-altitude target direction finding and positioning, and particularly relates to a super-resolution elevation angle estimation method for low-altitude multipath targets based on block sparse modeling. Background Art

[0002] Low-altitude penetration is one of the important threats faced by current air defense systems. Due to the influence of the multipath effect of radar target echoes, high-precision elevation angle estimation / altitude measurement has always been a difficult problem in the field of low-altitude detection.

[0003] In terms of the multipath propagation model, currently, the more widely studied and adopted one is to invert the height / position of the reflection point and the incident angle of the reflected wave in the digital elevation map based on the specular reflection hypothesis. However, generally, it is defaulted that the mirror plane is parallel to the ground plane, and the influence of terrain undulation on the reflection path is not considered; in recent years, some studies have also begun to consider the influence of diffuse reflection, but generally, the estimation algorithm is optimized to improve the robustness to the error / undulation of the incident angle of the reflected wave.

[0004] In terms of the estimation algorithm, with the rapid development of super-resolution direction finding technology, algorithms for estimating the elevation angle of low-altitude targets using super-resolution parameter estimation emerge in an endless stream. The more common algorithms can be divided into the following three categories: The first category is the adaptive beamforming algorithm, which can weaken the influence of diffuse reflection on the estimation performance by broadening the null of the reflected wave. However, generally, it requires a sufficient number of snapshots and non-coherent signals, and may involve multiple matrix inversion operations, with a relatively high computational complexity; the second category is the classical parameter estimation algorithm, such as the subspace algorithm or the maximum likelihood algorithm. The former deteriorates severely under the conditions of few snapshots and strongly correlated signals, while the latter requires the known or estimated number of targets; the third category is the compressed sensing or sparse recovery algorithm, which generally needs to adjust hyperparameters to ensure the estimation performance, and the existing models and methods do not consider the paired characteristics of the direct echo and the reflected wave. In addition, existing research usually estimates the multipath attenuation factor (or called the composite reflection coefficient, etc.) and the elevation angle in series. However, due to the coupling of parameters, the estimation errors of the two may affect each other, resulting in a decline in the estimation performance; to solve this problem, it is necessary to perform iteration in an alternating optimization manner, which will inevitably lead to a significant increase in computational complexity.

[0005] In summary, the existing low-altitude target elevation angle estimation methods are difficult to achieve a satisfactory compromise among model accuracy, estimation performance, and implementation complexity, and cannot guarantee the reliability under harsh conditions such as few snapshots and strongly correlated echoes. Summary of the Invention

[0006] The object of the present invention is to solve the problem that the elevation measurement performance deteriorates due to the multipath effect when the radar detects low-altitude targets in a complex geographical environment, and a super-resolution elevation estimation method for low-altitude multipath targets based on block sparse modeling is provided. On the basis of fully considering the influence of the actual terrain undulation of the position, this method accurately constructs a multipath propagation model, establishes a structured sparse representation model by combining the spatial sparsity of the target and the paired characteristics of the direct echo and the reflected wave, and proposes a block sparse recovery model and algorithm without hyperparameters and capable of rapid implementation for parameter estimation. Its model is more accurate and does not require the number of known targets and the surface reflection characteristics, has higher estimation accuracy and reliability under the conditions of strong echo correlation and few snapshots (even single snapshot), and while completing the high-precision and super-resolution elevation estimation of multiple multipath targets, it can also invert information such as the multipath attenuation factor and the inter-pulse amplitude and phase fluctuation characteristics of the target echo. Its application is not restricted by the arrangement method and system of the radar receiving array, and has a wider applicability.

[0007] In order to achieve the above object of the invention, the technical solution adopted to solve its technical problems is as follows:

[0008] A super-resolution elevation estimation method for low-altitude multipath targets based on block sparse modeling, and its specific steps are as follows:

[0009] Step 1: Based on the positioning and orientation information of the radar and the azimuth measurement value of the target by the radar, establish a rectangular coordinate system in the plane where the echo propagation path is located;

[0010] Step 2: Extract the terrain curve within a certain distance range in the direction of the target from the Geographic Information System (GIS);

[0011] Step 3: Divide the grid in the pitch dimension for the airspace of interest, analyze the validity of the reflection path corresponding to each pitch angle in combination with the distance measurement value and the terrain curve, and calculate the incident pitch angle of the reflected wave;

[0012] Step 4: Construct a structured redundant dictionary from the direct echo steering vectors and the reflected wave steering vectors corresponding to all pitch angles of interest;

[0013] Step 5: Establish a structured sparse representation model for the received data in the radar element domain / subarray domain / beam domain;

[0014] Step 6: Use the first-order fast algorithm to solve the convex combination optimization problem and recover the block sparse signal matrix;

[0015] Step 7: Obtain the elevation estimation value of the target based on the recovered block sparse signal, and further obtain the attenuation factor of the reflection path and the inter-pulse amplitude and phase fluctuation characteristics of the target echo as needed.

[0016] Furthermore, Step 1 includes the following content:

[0017] The positioning and orientation device using radar obtains information such as its own longitude, latitude, altitude, and the orientation of the antenna / array. Techniques such as monopulse angle measurement are used to estimate the azimuth angle of the target. A rectangular coordinate system is established in the plane perpendicular to the ground plane in this direction, with the radar's longitude / latitude and altitude of 0 as the coordinate origin. Denote the radar coordinates as (x R ,y R ), where x R = 0, and y R is the altitude corresponding to the radar antenna / array.

[0018] Furthermore, step 2 includes the following content:

[0019] Extract the terrain curve within the coordinate plane near the radar and within the range of the distance of interest from the Geographic Information System (GIS). Due to the discreteness of the digital map, this curve is actually a set composed of points divided at a certain resolution: {(x n ,y n )|n = 1,2,...,N}.

[0020] Furthermore, step 3 includes the following content:

[0021] Divide the airspace of interest into grids {θ l |l = 1,2,...,L} at a certain interval in the elevation dimension. Based on the measured distance value R of the radar to the target, calculate the position coordinates (x T ,y T ) = (x R + Rcosθ,y R + Rsinθ) when the target falls into each elevation grid, and further analyze the validity of the reflection path corresponding to each grid angle θ l and the incident elevation angle θ l ' of the corresponding reflected wave in combination with the terrain curve.

[0022] Specifically, the analysis methods for the reflection path validity and the incident elevation angle of the reflected wave are as follows:

[0023] ① Based on the wave path difference constraint between the direct echo of the target and the reflected wave, limit the range of potential reflection points: Given that when the wave path difference is large, the direct echo and the reflected wave can be distinguished by the distance, the length of the reflection path can be constrained to differ from the measured target distance value by no more than a certain number of distance resolution units, that is, d R + d T ≤ R + ΔR max , where d R is the distance from the reflection point to the radar, d T is the distance from the reflection point to the target, and ΔR max generally can take 1 to 3 times the distance resolution. The set of coordinates of the constrained potential reflection points Can be expressed as:

[0024]

[0025]

[0026]

[0027] ② Calculate the tangent slope at each point of the terrain curve: Here are two convenient calculation methods. One is to make the distances from the adjacent points on the left and right of this point to the tangent equal. Then, the tangent slope k n ,y n ) at the point (x n =(y n+1 -y n-1 ) / (x n+1 -x n-1 ); The other is to take the average of the connection slopes between this point and the adjacent points on the left and right. The corresponding slope expression is For general terrain undulation and map resolution, these two calculation methods are basically equivalent.

[0028] ③ Determine the position coordinates of the reflection point based on geometric relationships: Refer to Appendix Figure 3 . Under the assumption of ideal specular reflection, the necessary and sufficient condition for a point (x,y) to be a reflection point is that the angle bisector of the angles between the target to the reflection point and the reflection point to the radar coincides with the normal of the terrain curve at the reflection point. From the angle bisector theorem and the fixed-point ratio formula, the coordinates of the intersection point of the angle bisector and the line connecting the target and the radar are:

[0029]

[0030] If the connection line between this point and the reflection point (x,y) is perpendicular to the tangent of the terrain curve (with slope k), then there is:

[0031] k(y C -y)=x - x C ,

[0032] After arrangement, the necessary and sufficient condition for (x,y) to be a reflection point is:

[0033]

[0034] Therefore, the position coordinates of the reflection point corresponding to a target with elevation angle θ can be determined based on the following formula:

[0035]

[0036] Among them, represents the element in the set that makes the absolute value of f the smallest.

[0037] ④ Calculate the exact coordinates of the reflection point through interpolation: If the map resolution does not meet the requirements, the terrain curve can be interpolated at a desired high resolution within the neighborhood of the reflection point, and the above sub-steps ① to ③ are repeated. The finally obtained coordinates of the reflection point are still denoted as (x θ , y θ ).

[0038] ⑤ Judge the validity of the reflection path and calculate the incident elevation angle of the reflected wave: Considering that there may be obstructions on the reflection path, resulting in the reflected wave not being received by the radar, an effective reflection path must satisfy: the set of potential reflection point coordinates is non-empty, and for any the following equation holds.

[0039]

[0040] Therefore, it can be judged based on this whether multipath effect exists when the target elevation angle is θ. If the reflection path is effective, the incident elevation angle θ′ of the reflected wave is: θ′ = arctan[(y θ - y R ) / (x θ - x R )].[[]END]

[0041] As a special case of the above analysis, if it is roughly considered that the ground is absolutely flat (i.e., the mirror plane is parallel to the ground plane: y = y n ), the above analysis method of the reflection point coordinates can be simplified as follows: The wave path difference constraint in sub-step ① can be rewritten as Sub-step ② can be omitted (considering k n is always zero); The expression of the reflection point position coordinates in sub-step ③ is rewritten as:

[0042]

[0043] Furthermore, step 4 includes the following content:

[0044] Construct the following redundant dictionary using the direct echo steering vectors and reflected wave steering vectors corresponding to all elevation angle grid points (the number of grid points L should be much larger than the number of targets D):

[0045] A = [A1 A2 … A L

[0046]

[0047] ​Among them, \(a(\theta)\) represents the \(M\times1\) dimensional receiving array steering vector of the radar when the elevation angle of the incident signal is \(\theta\). It can be either the element domain steering vector (\(M\) corresponds to the number of antenna elements), or the subarray domain steering vector (\(M\) corresponds to the number of subarrays), or the beam domain steering vector (\(M\) corresponds to the number of beams). Denote the number of columns of \(A\) l as \(J\) l , and the number of columns of \(A\) is \(K\), then \(J\) l = 1, 2, Since the dictionary \(A\) is regularly composed of blocks \(A\) with the number of columns usually being 2 l , \(A\) can be called a structured dictionary.

[0048] Furthermore, step 5 includes the following contents:

[0049] Based on the \(M\times K\) dimensional dictionary \(A\) given in step 4, the radar (element domain / subarray domain / beam domain) received data \(X\) can be modeled as:

[0050] \(X = AS + N\)

[0051] Among them, \(X\) is composed of the range domain data after pulse compression of each channel (taking the cell where the range \(R\) is located), and its number of columns corresponds to the number of pulses \(P\) (also known as the sampling snapshot number); \(N\) is an \(M\times P\) dimensional additive white Gaussian noise; \(S\) is a \(K\times P\) dimensional signal matrix, which can be block - divided by rows according to the block - division rule of \(A\) \(S\) l is a \(J\) l \(\times P\) dimensional signal matrix, and the operator "T" represents the transpose of the matrix. The first row \(s\) l of \(S\) l1 corresponds to \(P\) samplings of the target echo with the elevation angle \(\theta\) l . If \(J\) l = 2, then the second row \(s\) l of \(S\) l2 corresponds to \(P\) samplings of its reflected wave; since the target echo and the reflected wave are from the same source, in the ideal case, they only differ by a complex constant \(\rho\) (called the multipath attenuation factor) determined by the spatial propagation characteristics and the surface reflection characteristics, that is, \(s\) l2 = \(\rho s\) l1 . It is not difficult to see that when there is a target at the grid angle \(\theta\) l , \(S\) l is composed of the target echo and its reflected wave and is a non - zero matrix; otherwise, \(S\) l = 0. Since the number of targets \(D\) is limited, and the custom - defined number of grid points \(L\) can be much larger than \(D\), \(S\) is necessarily a block - sparse matrix with most blocks being zero. Therefore, as long as the signal matrix \(S\) is recovered, the target elevation angle \(\theta\) can be determined by the positions of the non - zero blocks, and information such as the multipath attenuation factor \(\rho\) and the target inter - pulse amplitude - phase fluctuation characteristics can be obtained based on the data of the non - zero blocks.

[0052] Furthermore, step 6 includes the following contents:

[0053] To recover the block sparse signal matrix S, the following convex combination optimization problem is established:

[0054]

[0055] where, || || F represents the Frobenius norm of the matrix (degenerating to the l2 norm for vectors); the weights w of each block penalty term l are given by:

[0056]

[0057] where the operator "tr{}" represents the trace of the matrix, and the covariance matrix estimator R is given by:

[0058]

[0059] where, Diag(p) represents the diagonal matrix constructed with the vector p as the diagonal elements; the power vector p = [p1 p2 … p k T and each element in This optimization model has no hyperparameters that need to be manually adjusted or optimized, avoiding complex operations such as cross-validation.

[0060] Specifically, the fast solution of the above convex optimization problem can adopt first-order numerical optimization algorithms, such as the improved block coordinate descent method, the linearized alternating direction multiplier method (L-ADMM, or P-ADMM), the primal-dual hybrid gradient method (PDHG), or related improved algorithms. Their advantage is that the solution process does not involve complex operations such as matrix inversion, and their computational complexity will not increase explosively with the increase of the dictionary dimension and the number of snapshots. In addition, before the start of the algorithm iteration, the initial value S (0) of the matrix S to be recovered can be given by the linear minimum mean square error estimator (LMMSE): S (0) = Diag(p)·A H R -1 X.

[0061] Furthermore, step 7 includes the following content:

[0062] Based on the recovered block sparse signal matrix the spatial spectrum can be given by:

[0063] or Then the elevation angle estimation of the D targets can be determined according to the positions of the larger peaks in Z(θ l ), that is determined, i.e.,​ Furthermore, the signal matrix estimation can be corrected as follows:

[0064]

[0065] where the inter-pulse fluctuation characteristics of the d-th target echo can be given by the first row of (denoted as ); if has 2 rows (i.e., there is a reflection path for this target), denote its second row as then the corresponding multipath attenuation factor can be given by .

[0066] Due to the above technical solutions, the present invention has the following advantages and positive effects compared with the prior art:

[0067] (1) The proposed method fully considers the influence of actual terrain slope, ground object occlusion, etc. on the target echo propagation path, and can effectively improve the accuracy of the model.

[0068] (2) The proposed method transforms the elevation angle estimation into a block sparse signal recovery problem through the structured sparse representation modeling of radar received data. This model and method do not require the knowledge of the number of targets and the surface reflection characteristics, do not require manual parameter adjustment, are completely data adaptive, and have higher accuracy and reliability under the conditions of strong echo correlation and few snapshots (even single snapshot).

[0069] (3) The optimization problem designed by the proposed method can be quickly solved by a first-order numerical optimization algorithm, which can significantly reduce the computational complexity.

[0070] (4) The proposed method can not only complete the high-precision and super-resolution elevation angle estimation of multiple multipath targets, but also simultaneously obtain the attenuation factors of the corresponding reflection paths and the inter-pulse amplitude-phase fluctuation characteristics of the target echoes.

[0071] (5) The proposed method has no special requirements for the arrangement and system of the radar receiving array, and has wide applicability. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 is a flowchart of the method proposed by the present invention;

[0073] Figure 2 is a schematic diagram of the technical idea for analyzing the echo propagation path of the present invention;

[0074] Figure 3 is a schematic diagram of the geometric relationship between the direct echo and the reflection wave propagation paths of the target of the present invention;

[0075] Figure 4Analysis results of multipath effects of a certain low-altitude target by the method of the present invention under actual terrain;

[0076] Figure 5 Spatial spectrum given by the method of the present invention when the signal-to-noise ratio is 10 dB and the number of snapshots is 64;

[0077] Figure 6 Estimation of the amplitude fluctuation characteristics of the target echo given by the method of the present invention;

[0078] Figure 7 Estimation of the phase fluctuation characteristics of the target echo given by the method of the present invention;

[0079] Figure 8 Spatial spectrum given by the method of the present invention when the signal-to-noise ratio is -5 dB and the number of snapshots is 64;

[0080] Figure 9 Spatial spectrum given by the method of the present invention when the signal-to-noise ratio is 10 dB and the number of snapshots is 1. Specific implementation manners

[0081] The technical solution of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0082] Refer to Figure 1 the flowchart of, this embodiment discloses a low-altitude multipath target super-resolution elevation estimation method based on block sparse modeling, including the following steps:

[0083] Step 1: Based on the positioning and orientation information of the radar and the azimuth measurement value of the target by the radar, establish a rectangular coordinate system in the plane where the echo propagation path is located;

[0084] Step 2: Extract the terrain curve within a certain distance range in the direction where the target is located from the Geographic Information System (GIS);

[0085] Step 3: Divide the grid in the elevation dimension for the airspace of interest, analyze the validity of the reflection paths corresponding to each elevation angle in combination with the distance measurement value and the terrain curve, and calculate the incident elevation angle of the reflected wave;

[0086] Step 4: Construct a structured redundant dictionary from the direct echo steering vectors and the reflected wave steering vectors corresponding to all elevation angles of interest;

[0087] Step 5: Establish a structured sparse representation model for the received data in the radar element domain / subarray domain / beam domain;

[0088] Step 6: Solve the convex combination optimization problem using the first-order fast algorithm to recover the block-sparse signal matrix;

[0089] Step 7: Obtain the elevation angle estimate of the target based on the recovered block-sparse signal. If necessary, the attenuation factor of the reflection path and the inter-pulse amplitude-phase fluctuation characteristics of the target echo can be further obtained.

[0090] The specific descriptions of the above steps are as follows:

[0091] Furthermore, Step 1 includes the following:

[0092] Obtain information such as its own longitude, latitude, altitude, and antenna / array orientation through the radar's positioning and orientation device. Establish a rectangular coordinate system in the plane perpendicular to the ground plane at the azimuth of the target (the target azimuth angle measurement value is obtained by techniques such as monopulse angle measurement). Take the radar's longitude / latitude and altitude of 0 as the coordinate origin, the altitude of 0 as the x-axis, and the zenith as the positive y-axis direction. Denote the radar coordinates as (x R , y R ), where x R = 0, and y R is the altitude corresponding to the radar antenna / array.

[0093] Furthermore, Step 2 includes the following:

[0094] Based on the digital elevation information provided by the Geographic Information System (GIS), extract the terrain curve within the coordinate plane and within the range of interest near the radar. Due to the discreteness of the digital map, this curve is actually a set composed of points divided at a certain resolution (such as 30m, 90m): {(x n , y n ) n = 1, 2,..., N}. The digital map resolution used in this embodiment is 30 meters.

[0095] Furthermore, Step 3 includes the following:

[0096] Divide the airspace of interest in the elevation dimension (the elevation angle is defined as: 0 in the horizontal direction, positive for upward deviation, and negative for downward deviation) into grids {θ l | l = 1, 2,..., L} at a certain interval. Based on the radar's distance measurement value R of the target, calculate the position coordinates (x T , y T ) = (x R + Rcosθ, y R + Rsinθ) when the target falls into each elevation grid, and further analyze the validity of the reflection path corresponding to each grid angle θ l and its corresponding incident elevation angle θ l'. Here, the elevation grid only needs to cover the low-altitude airspace in combination with the actual beam width of the radar. It can be divided at equal angular intervals or at equal intervals in the sine space of the elevation angle. In this embodiment, the grid is divided at intervals of 1 / 80 beam width between -2.5 times the beam width and 5 times the beam width. Attached Figure 2 Intuitively introduces the technical idea of this step, that is, the curve formed by all the reflection points that satisfy the specular reflection hypothesis between the target and the radar and the actual terrain curve usually have only a very small number of intersection points, and their specific positions can be inferred according to geometric relationships.

[0097] Specifically, the analysis methods for the validity of the reflection path and the incident elevation angle of the reflected wave are as follows:

[0098] ① Based on the path difference constraint between the direct echo of the target and the reflected wave to limit the range of potential reflection points: Given that when the path difference is large, the direct echo and the reflected wave can be distinguished by distance, the length of the reflection path can be constrained to differ from the target distance measurement value by no more than a certain number of range resolution units, that is, d R +d T ≤R + ΔR max , where d R is the distance from the reflection point to the radar, and d T is the distance from the reflection point to the target; ΔR max Generally, it can take 1 to 3 times the range resolution. A certain margin needs to be reserved in its selection because there are errors in radar distance measurement, accuracy limitations in terrain curves, and possible energy leakage between adjacent range cells, etc. The set of coordinates of the potential reflection points after constraint can be expressed as:

[0099]

[0100]

[0101]

[0102] ② Calculate the tangent slope at each point of the terrain curve: Here are two convenient calculation methods. One is to make the distances from the left and right adjacent points of this point to the tangent line equal, then the tangent slope k n at (x n ) is (y n =(y n+1 -y n-1 ) / (x n+1 -x n-1 ); The other is to take the average of the connection slopes between this point and its left and right adjacent points, and the corresponding slope expression is In this embodiment, the latter method is used to calculate the tangent slope.

[0103] ③Determine the position coordinates of the reflection point based on geometric relationships: Refer to Appendix Figure 3 , under the assumption of ideal specular reflection, the necessary and sufficient condition for a point (x, y) to be a reflection point is that the angular bisector of the angle between the line connecting the target to the reflection point and the line connecting the reflection point to the radar coincides with the normal of the terrain curve at the reflection point. From the angle bisector theorem and the fixed-point ratio formula, the coordinates of the intersection point of the angular bisector and the line connecting the target and the radar are:

[0104]

[0105] If the line connecting this point and the reflection point (x, y) is perpendicular to the tangent of the terrain curve (with slope k), then there is:

[0106] k(y C - y) = x - x C

[0107] After arrangement, the necessary and sufficient condition for (x, y) to be a reflection point is:

[0108]

[0109] Therefore, the position coordinates of the reflection point corresponding to a target with an elevation angle of θ can be determined based on the following formula:

[0110]

[0111] where represents the element in the set that minimizes the absolute value of f.

[0112] ④Calculate the exact coordinates of the reflection point through interpolation: If the map resolution does not meet the requirements, then the terrain curve can be interpolated at a desired high resolution within the neighborhood of the reflection point, and the above sub-steps ① to ③ are repeated. The finally obtained reflection point coordinates are still denoted as (x θ , y θ ). In this embodiment, the selected neighborhood range is ±60 meters, the resolution of the interpolated map is 1 meter, and the interpolation method is cubic spline interpolation.

[0113] ⑤Judge the validity of the reflection path and calculate the incident elevation angle of the reflected wave: Considering that there may be obstructions on the reflection path, resulting in the reflected wave not being received by the radar, an effective reflection path must satisfy: the set of potential reflection point coordinates is non-empty, and for any x n < x θl 、 the following formula holds.

[0114]

[0115] Therefore, based on this, it can be determined whether multipath effects exist when the target pitch angle is θ. If the reflection path is valid, the incident pitch angle θ' of the reflected wave is: θ' = arctan[(y θ - y R ) / (x θ - x R )].

[0116] As a special case of the above analysis, if it is roughly considered that the ground is absolutely flat (i.e., the mirror plane is parallel to the ground plane: y = y n ), then the above analysis method of the reflection point coordinates can be simplified as follows: The path difference constraint in sub-step ① can be rewritten as Sub-step ② can be omitted (it is considered that k n is always zero); The expression of the position coordinates of the reflection point in sub-step ③ is rewritten as:

[0117]

[0118] Appendix Figure 4 shows the analysis results of the multipath effects of a target 8 km away from the radar with a pitch angle of -0.5°. Among them, the ground elevation at the radar is 252.414 m, the radar receiving array is 5 m away from the ground, and the range resolution of the radar is 3 m. It can be seen from the figure that after restricting the path difference within 6 m, the potential reflection points are limited within about 5.76 km; the final position coordinates of the reflection points obtained by interpolating the map are (813 m, 245.058 m), the incident pitch angle of the reflected wave is -0.105°, and there is no occlusion on the propagation path of the reflected wave.

[0119] Furthermore, step 4 includes the following content:

[0120] Construct the following redundant dictionary using the direct echo steering vectors and reflected wave steering vectors corresponding to all pitch angle grid points (the number of grid points L should be much larger than the number of targets D):

[0121] A = [A1 A2 … A L

[0122]

[0123] where a(θ) represents the M×1 dimensional receiving array steering vector of the radar when the incident angle of the incident signal is θ. It can be either the element domain steering vector (M corresponds to the number of antenna elements), or the sub-array domain steering vector (M corresponds to the number of sub-arrays), or the beam domain steering vector (M corresponds to the number of beams); Denote the number of columns of A l as J l , and the number of columns of A as K, then J l = 1, 2, Given that the dictionary A consists of blocks A with the number of columns usually being 2​l Regularly structured, A can be called a structured dictionary.

[0124] Actually, the basic cases of no reflection path and only one reflection path are considered here (L ≤ K ≤ 2L). If the actual terrain is more complex and there may be multiple reflection paths, only the corresponding extension of the dictionary is needed. In this embodiment, the number of grid points L is 601. After reflection path analysis, the final number of dictionary columns K is 1174; the radar receiving array is a uniformly linear array placed vertically, the number of antenna elements M is 256, the element spacing is half wavelength, and the specific expression of its steering vector is:

[0125] Further, step 5 includes the following content:

[0126] Based on the M×K-dimensional dictionary A given in step 4, the radar (element domain / subarray domain / beam domain) received data X can be modeled as:

[0127] X = AS + N

[0128] where X is composed of the range domain data after pulse compression of each channel (taking the unit where the range R is located), and its number of columns corresponds to the number of pulses P (also known as the sampling snapshot number); N is an M×P-dimensional additive white Gaussian noise; S is a K×P-dimensional signal matrix, which can be block-divided by rows according to the block rule of A as S l is J l ×P-dimensional signal matrix, and the operator "T" represents the transpose of the matrix. The first row s l of S l1 corresponds to P samplings of the target echo with the elevation angle θ l . If J l = 2, the second row s l of S l2 corresponds to P samplings of its reflected wave; since the direct echo of the target and the reflected wave are from the same source, in the ideal case, the two only differ by a complex constant ρ (called the multipath attenuation factor), that is, s l2 ≈ ρs l1 (usually |ρ| < 1). Actually, ρ is the result of the combined action of factors such as the attenuation caused by the difference in wave path and propagation medium, the surface reflection characteristics, and the receiving gain at different incident angles. It is not difficult to see that when there is a target at the grid angle θ l , S l is composed of the target echo and its reflected wave and is a non-zero matrix; otherwise, S l= 0. Since the number of targets D is finite, and the custom number of grid points L can be much larger than D, S must be a block-sparse matrix with most blocks being zero. Therefore, as long as the signal matrix S is recovered, the target elevation angle θ can be determined by the positions of the non-zero blocks, and information such as the multipath attenuation factor ρ and the inter-pulse amplitude and phase fluctuation characteristics of the target can be obtained based on the data of the non-zero blocks.

[0129] Further, step 6 includes the following content:

[0130] To recover the block-sparse signal matrix S, the following convex combination optimization problem is established:

[0131]

[0132] where, |||| F represents the Frobenius norm of the matrix (degenerating to the l2 norm for vectors); the weight w of each block penalty term l is given by the following formula:

[0133]

[0134] where the operator "tr{}" represents the trace of the matrix, and the covariance matrix estimator R is given by the following formula:

[0135]

[0136] where, Diag(p) represents a diagonal matrix constructed with the vector p as the diagonal elements; the power vector p = [p1 p2 … p k T each element in This optimization model has no hyperparameters that need to be manually adjusted or optimized, avoiding complex operations such as cross-validation.

[0137] Specifically, the fast solution of the above convex optimization problem can adopt a first-order numerical optimization algorithm, such as the improved block coordinate descent method, the linearized alternating direction multiplier method (L-ADMM, or P-ADMM), the primal-dual hybrid gradient method (PDHG), or related improved algorithms. Its advantage is that the solution process does not involve complex operations such as matrix inversion, and its computational complexity does not increase explosively with the increase of the dictionary dimension and the number of snapshots. In addition, before the algorithm iteration starts, the initial value S of the matrix S to be recovered (0) can be given by the linear minimum mean square error estimator (LMMSE): S (0) = Diag(p)·A H R -1 X. This embodiment adopts the PDHG algorithm accelerated by the modified penalty function and the Nesterov method for solution.

[0138] ​Further, step 7 includes the following content:

[0139] Based on the restored block-sparse signal matrix The spatial spectrum can be given by the following formula:

[0140] or Then the elevation angle estimation of D targets can be determined according to the positions of the larger peaks in Z(θ l ) That is Further, the signal matrix estimation can be corrected as follows:

[0141]

[0142] where Then the inter-pulse fluctuation characteristics of the echo of the d-th target can be given by the first row of (denoted as ); if has 2 rows (i.e., there is a reflection path for this target), denote its second row as then the corresponding multipath attenuation factor can be given by .

[0143] In summary, a low-altitude multipath target super-resolution elevation angle estimation method based on block sparse modeling proposed by the present invention is applicable to high-precision direction finding and positioning of radar low-altitude targets under complex terrain conditions. This method effectively utilizes terrain information to improve the accuracy of the multipath propagation model, and establishes a structured sparse representation model and proposes a fast adaptive block sparse recovery algorithm to ensure the estimation performance and reliability under the conditions of strong echo correlation and few snapshots (even single snapshot); the proposed method does not require the number of known targets and surface reflection characteristics, and while completing the high-precision and super-resolution elevation angle estimation of multiple targets, it can also provide information such as multipath attenuation factors and inter-pulse fluctuation characteristics, and its application is not restricted by the arrangement mode and system of the radar receiving array.

[0144] The effects of the present invention are verified by the following embodiments and simulation tests.

[0145] Figure 5The spatial spectrum comparison of different methods is given when the number of targets is 2, the signal-to-noise ratio is 10 dB, and the number of snapshots is 64. Among them, the targets are 8 km away from the radar, and their elevation angles are ±0.05° respectively. The receiving array is a uniform linear array with 256 elements and a half-wavelength interval, and the beam width is about 0.4°. The signal-to-noise ratio is defined as the ratio of the direct echo to the noise power of a certain target in a single receiving channel. Since the radar receiving data here is the range-domain data of specific cells after pulse compression, the number of snapshots corresponds to the number of pulses. The "unstructured model" means that all possible target elevation angles and reflection wave elevation angles are incorporated into the grid without distinction (the grid ranges of the "structured model" and the "unstructured model" are -1° to 2° and -2.87° to 2° respectively), and they are uniformly regarded as targets for elevation angle estimation. It can be seen that when the "unstructured model" is adopted, whether the multiple signal classification (MUSIC) algorithm or the sparse recovery algorithm is used for angle estimation, although there are spectral peaks at the elevation angles where the targets are located and the elevation angles where the reflection waves are located (at -0.98° and -1.35°), there may also be relatively high false peaks at other angles (especially the "noise floor" of the MUSIC algorithm is significantly elevated, which is caused by the strong correlation between the direct echo and the reflection wave of the same source), affecting the judgment of the final estimation result. In contrast, when the "structured model" proposed in this patent is adopted, both the MUSIC algorithm and the proposed algorithm can obtain a "clean" spatial spectrum and accurate super-resolution angle estimation (the target interval is about 1 / 4 of the beam width), and there are almost no false peaks in the proposed method.

[0146] Figure 6 and Figure 7 gives the amplitude-phase fluctuation characteristics of the target echo estimated by the method proposed in this patent under the simulation conditions of Figure 5 It can be seen that while obtaining the elevation angle estimation, the proposed method can also accurately invert the amplitude-phase fluctuation characteristics of the target echo. In addition, in the simulation, the amplitudes and phases of the multipath attenuation factors of the two targets are (0.5918, 103.3593°) and (0.5303, -55.2548°) respectively, and the estimated values given by the proposed method are (0.5918, 103.3104°) and (0.5314, -55.3880°) respectively. It can be seen that the proposed method can also accurately estimate the multipath attenuation factors of each target.

[0147] Figure 8The spatial spectrum of the proposed method is given when the signal-to-noise ratio is -5 dB and the number of snapshots is 64. It can be seen that under the condition of low signal-to-noise ratio, although there are spectral peaks at the target angles in the spatial spectra given by the MUSIC algorithm based on the "unstructured model" and the sparse recovery algorithm based on the "structured model", there are also a large number of false peaks at other angles, and their heights may even exceed the height of the target spectral peak, so that a reliable elevation angle estimation result cannot be obtained. In contrast, the method proposed in the present invention can still achieve reliable and accurate super-resolution elevation angle estimation.

[0148] Figure 9 The spatial spectrum of the proposed method is given when the signal-to-noise ratio is 10 dB and the number of snapshots is 1. It can be seen that under the condition of single snapshot, there is still the problem of false peaks in the sparse recovery algorithm based on the "structured model", and a reliable elevation angle estimation result cannot be provided. And the MUSIC algorithm based on the "unstructured model" fails completely due to the rank deficiency of the covariance matrix (in this case, preprocessing such as spatial smoothing needs to be used, but such operations will lead to aperture loss). In contrast, the method proposed in this patent can still provide reliable and accurate estimation results and achieve monopulse super-resolution elevation angle estimation of low-altitude targets.

[0149] As mentioned above, it is only the preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed by the present invention should be covered by the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims

1. A method for super-resolution elevation estimation of low-altitude multipath targets based on block sparse modeling, characterized in that Including: Step 1, based on the positioning and orientation information of the radar and the azimuth measurement value of the radar for the target, establish a rectangular coordinate system in the plane where the echo propagation path is located; among them, denote the radar coordinates as (x R , y R ); Step 2, extract the terrain curve within a preset distance range in the direction where the target is located from the geographic information system; wherein, the terrain curve is a set composed of points divided at a certain resolution: {(x n ,y n )|n = 1, 2, ..., N}; Step 3: Perform grid division in the pitch dimension for the airspace of interest. Combine the distance measurement values and the terrain curve to analyze the validity of the reflection paths corresponding to each pitch angle, and calculate the incident pitch angle of the reflected wave. Among them, divide the airspace of interest into grids {θ l |l = 1, 2, ..., L} in the pitch dimension at a certain interval. Based on the distance measurement value R of the radar for the target, calculate the position coordinates (x T , y T ) = (x R + Rcosθ, y R + Rsinθ) when the target falls into each pitch grid, and further analyze the validity of the reflection paths corresponding to each grid angle θ l and its corresponding incident pitch angle θ l ' of the reflected wave. The analysis methods for the validity of the reflection path and the incident pitch angle of the reflected wave include: ① The range of potential reflection points is limited based on the path difference constraint between the target direct echo and the reflected wave, where d R is the distance from the reflection point to the radar, d T is the distance from the reflection point to the target, ΔR max Take 1 to 3 times the distance resolution, and constrain the potential reflection point coordinate set It is expressed as: ② Calculate the tangent slope k at each point of the terrain curve n ; ③ Determine the position coordinates of the reflection point based on geometric relationships. Among them, the position coordinates of the reflection point corresponding to the target with an elevation angle of θ are determined based on the following formula: Among them, represents the set the element in which minimizes the absolute value of f; ④ Calculate the exact coordinates of the reflection point through interpolation, and the finally obtained coordinates of the reflection point are still denoted as (x θ , y θ ); ⑤ Determine the validity of the reflection path and calculate the incident elevation angle of the reflected wave: Considering that there may be obstructions on the reflection path, resulting in the reflected wave not being received by the radar, an effective reflection path must satisfy: the set of potential reflection point coordinates is non-empty, and for any the following equation holds: Based on this, it is determined whether multipath effect exists when the target pitch angle is θ; if the reflection path is valid, the incident pitch angle θ' of the reflected wave is: θ' = arctan[(y θ -y R ) / (x θ -x R )]; Step 4: Construct a structured redundant dictionary from the direct echo steering vectors and reflected wave steering vectors corresponding to all elevation angles of interest. Among them, the following redundant dictionary is constructed using the direct echo steering vectors and reflected wave steering vectors corresponding to all elevation angle grid points. The number of grid points L is much larger than the number of targets D: A = [A1 A2 … A L ​ Among them, \(a(\theta)\) represents the \(M\times1\) dimensional receiving array steering vector of the radar when the elevation angle of the incident signal is \(\theta\). If \(a(\theta)\) is the steering vector in the element domain, then \(M\) corresponds to the number of antenna elements; if \(a(\theta)\) is the steering vector in the subarray domain, then \(M\) corresponds to the number of subarrays; if \(a(\theta)\) is the steering vector in the beam domain, then \(M\) corresponds to the number of beams. Denote the number of columns of \(A\) l as \(J\) l , and the number of columns of \(A\) is \(K\), then \(J\) l = 1, 2, The dictionary \(A\) is regularly composed of blocks \(A\) with the number of columns usually being 2 l , and \(A\) is called a structured dictionary; Step 5: Establish a structured sparse representation model for the received data in the radar element domain, subarray domain or beam domain; Step 6: Use the first-order fast algorithm to solve the convex combination optimization problem and recover the block sparse signal matrix; Step 7: Obtain the estimated elevation angle of the target based on the recovered block sparse signal matrix. If needed, the attenuation factor of the reflection path and the inter-pulse amplitude-phase fluctuation characteristics of the target echo can be further obtained.

2. The method for estimating the super-resolution elevation angle of low-altitude multipath targets based on block sparse modeling according to claim 1, wherein Step 3 includes: If we roughly assume that the ground is absolutely flat, that is, the mirror plane is parallel to the ground plane: y = y n , then the above analysis method of the reflection point coordinates is simplified as follows: The path difference constraint in sub-step ① is rewritten as Sub-step ② is omitted, and it is considered that k n is always zero; the expression of the position coordinates of the reflection point in sub-step ③ is rewritten as:

3. The method for super-resolution elevation angle estimation of low-altitude multipath targets based on block sparse modeling according to claim 1, characterized in that, Step 5 includes: Based on the M×K-dimensional dictionary A given in Step 4, the received data X in the radar element domain, subarray domain or beam domain is modeled as: X = AS + N, Among them, X is composed of the range-domain data after pulse compression of each channel, and its number of columns corresponds to the number of pulses P; N is an M×P-dimensional additive white Gaussian noise; S is a K×P-dimensional signal matrix, which is block-divided by rows according to the block rule of A as S l is a J l ×P-dimensional signal matrix, and the operator "T" represents the transpose of the matrix; the first row s l of S l1 corresponds to P samples of the target echo with the elevation angle of θ l . If J l =2, then the second row s l of S l2 corresponds to P samples of its reflected wave; since the target echo and the reflected wave are from the same source, in the ideal case, the two only differ by a complex constant ρ determined by the spatial propagation characteristics and the surface reflection characteristics, that is, s l2 =ρs l1 ; when there is a target at the grid angle θ l , S l is composed of the target echo and its reflected wave and is a non-zero matrix; otherwise, S l =0; since the number of targets D is limited and the number of custom-defined grid points L is much larger than D, S is necessarily a block-sparse matrix with most blocks being zero. Therefore, as long as the signal matrix S is recovered, the target elevation angle θ can be determined by the positions of the non-zero blocks, and the multipath attenuation factor ρ and the information on the target's inter-pulse amplitude and phase fluctuation characteristics can be obtained based on the data of the non-zero blocks.

4. The method for estimating the super-resolution elevation angle of low-altitude multipath targets based on block sparse modeling according to claim 3, characterized in that Step 6 includes: To recover the block sparse signal matrix S, the following convex combination optimization problem is established: where, || || F represents the Frobenius norm of the matrix, which degenerates to the l2 norm for vectors; the weight w of each block penalty term l is given by the following formula: Among them, the operator "tr{}" represents the trace of the matrix, and the covariance matrix estimator R is given by the following formula: Among them, Diag(p) represents a diagonal matrix constructed with the vector p as the diagonal elements; the power vector p = [p1 p2 … p k T and the elements in ​ 5. A method for super-resolution elevation angle estimation of low-altitude multipath targets based on block sparse modeling according to claim 4, characterized in that The fast solution of the convex optimization problem uses first-order numerical optimization algorithms, including: improved block coordinate descent method, linearized alternating direction multiplier method, and primal-dual hybrid gradient method; before the start of the algorithm iteration, the initial value S of the matrix S to be restored (0) is given by the linear minimum mean square error estimator: S (0) = Diag(p)·A H R -1 X.

6. The method for estimating the super-resolution elevation angle of low-altitude multipath targets based on block sparse modeling according to claim 1, wherein, Step 7 includes: Recovery-based block sparse signal matrix The spatial spectrum is given by: or Then the elevation angle estimation of D targets can be based on the position of the larger peak in Z(θ l ), that is determined, namely The following correction is made to the signal matrix estimation: Among them, the inter-pulse fluctuation characteristic of the d-th target echo is given by the first row of denoted as ; if has 2 rows, that is, there is a reflection path for this target, and its second row is denoted as then the corresponding multipath attenuation factor is given by as follows.

Images