Factor graph based distributed spaceborne radar networked height estimation of air targets
By using a factor graph model and a nonparametric confidence propagation algorithm, the problem of large elevation angle measurement error in spaceborne radar networking was solved, achieving high-precision and highly robust aerial target altitude estimation.
Patent Information
- Application Number
- CN202310220001.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2023-03-09
- Publication Date
- 2026-01-30
- Estimated Expiration
- 2043-03-09
AI Technical Summary
Single-satellite radars have large pitch angle measurement errors, making it difficult to accurately track aerial targets and estimate altitude. Traditional centralized nonlinear filtering algorithms have poor robustness.
The distributed spaceborne radar networking method based on factor graphs constructs an airborne target motion model and a spaceborne radar measurement model, introduces local variable coupling function relationships, uses a nonparametric confidence propagation algorithm to calculate the posterior marginal distribution, and combines iterative approximation method to estimate the target altitude.
It improves the accuracy and robustness of aerial target altitude estimation. The algorithm has a fast convergence speed, high robustness, and higher accuracy than traditional algorithms.
Smart Images
Figure CN116243299B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of spaceborne radar target tracking and relates to a method for estimating the altitude of airborne targets in a distributed spaceborne radar network based on factor graphs. Technical Background
[0002] Spaceborne radar possesses all-weather, all-time strategic and tactical early warning capabilities, and enjoys advantages such as being unaffected by the Earth's curvature and less vulnerable to attack, making it crucial in early warning and defense systems. When a single spaceborne radar tracks and locates aerial targets, the target elevation angle measurement error is generally larger than the radial distance and azimuth angle measurement errors, resulting in a lack of aerial target altitude estimation capabilities. This lack of aerial target altitude information significantly amplifies target location errors and hinders target attribute identification and threat assessment. Therefore, networking multiple spaceborne radars and jointly utilizing their radial distance and azimuth angle measurements, combined with target motion models, and employing nonlinear filtering methods to effectively estimate aerial target altitude is of great significance.
[0003] Spaceborne radar target tracking utilizes nonlinear filtering algorithms to achieve high-precision target tracking based on the target's motion model and radar measurements. However, spaceborne radar target tracking involves complex coordinate transformations and nonlinearities in radar measurements, requiring high tracking accuracy. Furthermore, due to the unique characteristics of the radar platform, traditional centralized nonlinear filtering algorithms exhibit poor robustness. Therefore, addressing the challenges of accurate target tracking by a single spaceborne radar and the nonlinear issues involved in the tracking process, it is necessary to research high-precision, robust distributed nonlinear filtering methods for networked spaceborne radar systems to achieve accurate and stable altitude estimation of aerial targets. Summary of the Invention
[0004] The purpose of this invention is to provide a method for estimating the altitude of airborne targets in a distributed spaceborne radar network based on factor graphs, in order to solve the problems in the prior art of large elevation angle measurement errors for single spaceborne radars, difficulty in accurately tracking airborne targets, and inaccurate altitude estimation.
[0005] The technical solution adopted in this invention is a method for estimating the altitude of airborne targets in a distributed spaceborne radar network based on factor graphs, comprising the following steps:
[0006] Step 1: In the geocentric fourth equatorial coordinate system, construct an airborne target motion model and a spaceborne radar measurement model. Combine the airborne target motion model and the spaceborne radar measurement model to construct a spaceborne radar target tracking model based on a factor graph, i.e., a factor graph.
[0007] Step 2: Based on the factor graph, introduce the coupling function relationship between local variables of the air target and obtain the posterior marginal distribution of the local variables;
[0008] The posterior marginal distribution of the local variable is calculated using a nonparametric confidence propagation algorithm, and the estimated value of the local variable is obtained based on the calculation result of the posterior marginal distribution.
[0009] Step 3: Based on the estimated values of the local variables, the target height of the aerial target in the geodetic coordinate system is obtained using an iterative approximation method.
[0010] Furthermore, the specific content of step 1 is as follows:
[0011] S1.1 Constructing an aerial target motion model in the geocentric fourth equatorial coordinate system:
[0012] x k =Fx k-1 +ω k-1 ,
[0013] Where, x k Let be the target's position, velocity, and state vector in the fourth geocentric equatorial coordinate system at time k.
[0014] x k-1 Let be the target's position, velocity, and state vector in the fourth geocentric equatorial coordinate system at time k-1. ω k-1 Assuming process noise
[0015] F is the target state transition matrix:
[0016]
[0017] in, It is the matrix direct product operator, where I3 is the three-dimensional identity matrix and T is the sampling interval;
[0018] S1.2 Constructing a spaceborne radar measurement model in the fourth geocentric equatorial coordinate system:
[0019]
[0020] in, Let be the target's position vector in the coordinate system of the nth radar array. v k ~N(0,R), Let N be the radial range and azimuth measurement errors of the N radars; the radial range and azimuth measurement errors of the nth radar are... Then the measurements of N radars at time k are
[0021] S1.3, Combine the target motion model and the radar measurement model to perform factor graph modeling:
[0022] For the aforementioned airborne target motion model and the aforementioned spaceborne radar measurement model, the variables... and The joint probability density function can be decomposed as:
[0023]
[0024] Where p(x) k |x k-1 Let p(y) be the state transition probability density function. k |x k Let be the likelihood probability density function, and then describe the decomposition form of the joint probability density function according to the factor graph to obtain the factor graph model.
[0025] Furthermore, the specific process of obtaining the posterior marginal distribution of the local variables of the aerial target by introducing the coupling function relationship between local variables in step 2 is as follows:
[0026] In the factor graph, a coupling factor node g is introduced. ji To represent the global variable x k Two replicated state variables and The relationship between adjacent nodes j and i
[0027] Then, at time k, the local variables The posterior marginal distribution is:
[0028]
[0029] in For state transition messages, For measuring messages, For coupled messages, N j Let j be the set of all neighboring nodes of node j.
[0030] Furthermore, in step 2, the local variables at time k are calculated based on the factor graph model. The specific process of the posterior marginal distribution is as follows:
[0031] At time k, compute the particle-like form of the state transition message. and the particle form of measurement messages
[0032] Calculate coupled messages using message iteration. In the l-th iteration, kernel density estimation is used to calculate the message product of all neighboring nodes of node j: Update based on coupling distribution relationship After a total of L iterations, the coupled message is obtained.
[0033] The posterior distribution is calculated using the kernel density estimation method:
[0034] Then, the local variables are obtained based on the posterior distribution. The estimated value.
[0035] The beneficial effects of this invention are as follows: This invention describes the motion model of an aerial target and the measurement model of a spaceborne radar in the fourth geocentric equatorial coordinate system. It also models the target tracking of a networked spaceborne radar system based on factor graphs, introducing coupling parameters between local variables to derive a local variable factor graph model based on coupling function relationships. Considering the large linearization error in solving the Jacobian matrix using the traditional extended Kalman filter algorithm, a particle-based message passing algorithm based on importance sampling, i.e., a nonparametric confidence propagation algorithm, is proposed by combining confidence propagation and particle filtering methods. Finally, the target altitude is calculated based on the state estimate of the aerial target in the fourth geocentric equatorial coordinate system. Compared with traditional algorithms, this invention introduces coupling function relationships, resulting in faster convergence, higher robustness, and higher altitude estimation accuracy. Attached Figure Description
[0036] Figure 1 The message passing factor diagram for target tracking of spaceborne radar in the global state in this invention;
[0037] Figure 2 The message passing factor graph between local variables in this invention;
[0038] Figure 3 The message passing factor graph between all nodes is introduced in this invention after the coupling parameter is introduced;
[0039] Figure 4 This is a comparison chart of the root mean square error of target height estimation under four algorithms: EKF, DCEKF, NBP, and geometric method in the embodiment.
[0040] Figure 5-(a) is a comparison of the root mean square error of target height estimation under different coupling parameters when the number of samples is 400 in the NBP algorithm of the embodiment;
[0041] Figure 5-(b) is a comparison of the root mean square error of target height estimation under different coupling parameters when the number of samples is 800 in the NBP algorithm of the embodiment;
[0042] Figure 6 The example uses the NBP algorithm to estimate the target height under different measurement errors, and the root mean square error is plotted. Detailed Implementation
[0043] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.
[0044] This invention provides a method for estimating the altitude of airborne targets in a distributed spaceborne radar network based on factor graphs, including the following:
[0045] Step 1: In the geocentric fourth equatorial coordinate system, construct the airborne target motion model and the spaceborne radar measurement model, and combine the airborne target motion model and the spaceborne radar measurement model to construct a spaceborne radar target tracking model based on the factor graph, i.e., the factor graph.
[0046] Step 2: In the factor graph model, the coupling function relationship between local variables is introduced to obtain the posterior marginal distribution of the local variables of the aerial target;
[0047] The posterior marginal distribution of the local variables of the aerial target is calculated using a nonparametric confidence propagation algorithm, and the estimated value of the variable node is obtained based on the calculation result of the posterior marginal distribution.
[0048] Step 3: Based on the estimated values of the variable nodes, the target height of the aerial target in the geodetic coordinate system is obtained by using an iterative approximation method.
[0049] In some embodiments, the specific content of step 1 is as follows:
[0050] S1.1 Assuming the target is flying at a constant speed, the motion model of the aerial target constructed in the geocentric fourth equatorial coordinate system can be described as follows:
[0051] x k =Fx k-1 +ω k-1 (1),
[0052] Where, x k Let x be the target's position, velocity, and state vector in the fourth geocentric equatorial coordinate system at time k. k-1 Let be the target's position, velocity, and state vector in the fourth geocentric equatorial coordinate system at time k-1, and F is the target state transition matrix:
[0053]
[0054] in, It is the matrix direct product operator; I3 is the three-dimensional identity matrix; T is the sampling interval. ω k-1 Assuming process noise
[0055] S1.2 Constructing a spaceborne radar measurement model in the fourth geocentric equatorial coordinate system:
[0056]
[0057] in, Let be the target's position vector in the coordinate system of the nth radar array. v k ~N(0,R), Let N be the radial range and azimuth measurement errors of the N radars; the radial range and azimuth measurement errors of the nth radar are... Then the measurements of N radars at time k are
[0058] The specific process of constructing the spaceborne radar measurement model is as follows: Assume that at time k, a total of N radars measure the target, and the coordinates of each radar are represented as (X... n ,Y n Z n ), n=1,2,…,N. The target measurement information obtained by each radar is radial range and azimuth At the nth radar node, the target's motion model can be approximated as:
[0059]
[0060] in Let be the local position and velocity state vector of the nth radar node target. For the local process noise of the nth radar node, assume
[0061] Then the measurement equation for the nth radar at time k can be expressed as:
[0062]
[0063] in Let be the position vector of the target in the coordinate system of the nth radar array, and Where g(·) represents the transformation relationship from the geocentric fourth orbital coordinate system to the radar array coordinate system. and Let the measurement error of the nth radar at time k be assumed.
[0064] Define the measurements of all radars at time k as The global measurement model can then be expressed as:
[0065]
[0066] Among them, v k ~N(0,R),
[0067] S1.3, Combine the target motion model and the radar measurement model to perform factor graph modeling:
[0068] For the aforementioned airborne target motion model and the aforementioned spaceborne radar measurement model, the variables... and The joint probability density function can be decomposed as:
[0069]
[0070] Where p(x) k |x k-1 Let p(y) be the state transition probability density function. k |x k Let be the likelihood probability density function. Then, based on the factor graph description of the joint probability density function, we obtain the factor graph model as follows: Figure 1 As shown, where f k and h k These are the factor nodes representing the state transition function and the likelihood function.
[0071] In step 1, it should be noted that since the target measurement information (radial distance, azimuth, and elevation angle) provided by the spaceborne radar is located in the radar array coordinate system, while the satellite platform and target position information are in the geocentric fourth equatorial coordinate system, and this invention selects the geocentric fourth equatorial coordinate system as the filtering state space (refer to GB / T 32296-2015), the distributed spaceborne radar network air target altitude estimation method based on factor graphs in this invention involves coordinate transformation from the radar array measurement coordinate system to the geocentric fourth coordinate system. The specific transformation process is as follows:
[0072] (1) Transformation of radar measurement coordinate system to antenna array coordinate system:
[0073]
[0074] Where, r z A z E and E represent the radial distance, azimuth, and elevation angle of the aerial target as measured by the radar, respectively; x a y a , z a This represents the position of the aerial target in the antenna array coordinate system.
[0075] (2) Transformation of the antenna array coordinate system to the satellite body coordinate system:
[0076]
[0077] Where, γ c , and ψ c M1, M2, and M3 are the roll angle, elevation angle, and yaw angle of the radar antenna array, respectively; M1, M2, and M3 are rotation matrices that rotate counterclockwise around the x, y, and z axes by a certain angle according to the right-hand rule, respectively.b y b , z b This refers to the position of an aerial target in the satellite's coordinate system.
[0078] (3) Transformation from the satellite body coordinate system to the satellite orbit coordinate system:
[0079]
[0080] Where, γ s , and ψ s X1, Y1, and Z1 represent the roll angle, pitch angle, and yaw angle of the satellite attitude, respectively; X1, Y1, and Z1 represent the positions of the aerial target in the satellite orbital coordinate system.
[0081] (4) Transformation of satellite orbital coordinate system to geocentric second orbital coordinate system:
[0082] r0 is the distance from the satellite to the Earth's center; f is the satellite's true anomaly angle; x2, y2, z2 are the positions of the target in the second geocentric orbital coordinate system.
[0083] (5) From the second geocentric orbital coordinate system to the first geocentric equatorial coordinate system:
[0084]
[0085] ω, Ω and η i X1, Y1, and Z1 represent the satellite's perigee argument, right ascension of the ascending node, and orbital inclination, respectively; X1, Y1, and Z1 represent the target's position in the geocentric first equatorial coordinate system.
[0086] (6) From the first geocentric equatorial coordinate system to the second geocentric equatorial coordinate system:
[0087]
[0088] Where [A], [B], [C], and [D] are the polar motion, rotation, nutation, and precession matrices, respectively. X2, Y2, and Z2 represent the target's position in the geocentric second equatorial coordinate system.
[0089] (7) From the second geocentric equatorial coordinate system to the fourth geocentric equatorial coordinate system:
[0090]
[0091] Where, ω e Let t be the Earth's rotational angular velocity, calculated from the epoch beginning at t0. X4, Y4, and Z4 represent the target's position in the fourth geocentric equatorial coordinate system.
[0092] In some embodiments, a coupling factor node g is introduced into the factor graph. jiTo represent the global variable x k Two replicated state variables and The relationship between adjacent nodes j and i
[0093] Then, at time k, the local variables The posterior marginal distribution is:
[0094]
[0095] in, For state transition messages, For measuring messages, For coupled messages, N j Let j be the set of all neighboring nodes of node j.
[0096] The specific implementation process of this step is as follows:
[0097] like Figure 1 In the factor graph shown, the global variable x is obtained at time k through recursive calculation on the factor graph. k The posterior marginal distribution and the global variable x k+1 Predicted distribution:
[0098]
[0099] Among them, the definition The update obtained from the posterior distribution at time k is the same as the state estimate obtained by the traditional KF method.
[0100] Equations (16) and (17) are for updating and predicting the target global state variables.
[0101] Based on the spaceborne radar measurement model obtained in step 1, a coupling factor node g is introduced. ji To represent the global variable x k Two replicated state variables and The relationship between adjacent nodes j and i, where N j Let j represent the set of neighbors.
[0102] g ij Defined as an exponential function:
[0103]
[0104] Where β is the coupling parameter, and β > 0.
[0105] Figure 2 Let u,i∈N be the message passing factor graph between local variables. jAs β increases, neighboring variables... and The correlation between them becomes increasingly stronger, thus promoting the consistency of the state distribution in the factor graph. As β→∞, we can obtain... At this point, the distribution of the local state variables is the same throughout the factor graph, that is... And it is consistent with the global variable distribution p(x) k |y k )same.
[0106] Therefore, at time k, the local variables The posterior marginal distribution is:
[0107]
[0108] In some embodiments, in step 2, the local variables at time k are calculated based on the factor graph model. The process of the posterior marginal distribution is as follows:
[0109] At time k, compute the particle-like form of the state transition message. and the particle form of measurement messages
[0110] Calculate coupled messages using message iteration. In the l-th iteration, kernel density estimation is used to calculate the message product of all neighboring nodes of node j: Update based on coupling distribution relationship After a total of L iterations, the coupled message is obtained.
[0111] The posterior distribution is calculated using the kernel density estimation method:
[0112] Then, the local variables are obtained based on the posterior distribution. The estimated value.
[0113] The specific implementation process of this step is as follows:
[0114] The traditional extended Kalman filter algorithm makes a quasi-linear assumption on the nonlinearity of the measurement model to represent messages using mean and variance. However, in strongly nonlinear problems, the estimation performance of this assumption cannot be guaranteed. Therefore, for the nonlinear, non-Gaussian probability distribution in the factor graph, nonparametric belief propagation (NBP) is used to calculate message passing in the factor graph, combining belief propagation and particle filtering. After adding the coupling function in step S2.1, at time k, the connection form of all factor nodes and variable nodes is as follows: Figure 3 As shown.
[0115] From factor node f k-1 Pass to variable node The message can be written as:
[0116] A particle-based message passing algorithm is used, and an importance sampling method is employed to construct messages. The particle-based representation. Without loss of generality, the proposed distribution is chosen. Then the s-th particle sampled from the proposed distribution is Then message The particle representation is as follows:
[0117]
[0118] in,
[0119] and
[0120]
[0121] From factor node g ij Pass to variable node The message can be written as:
[0122] Importance sampling recommendation distribution selection coupling function Then message The particle representation is as follows:
[0123]
[0124] in,
[0125] and
[0126]
[0127] From factor nodes Pass to variable node The message is Recommendational distribution selection of importance sampling likelihood function However, due to the observation equation Only with the location of the target Relevant, therefore, suggested distribution It can be simplified to The location of the target It can be represented as:
[0128]
[0129] Where g -1(·) represents the transformation relationship from the radar array coordinate system to the geocentric fourth orbit coordinate system, where d and α are the target radial distance and azimuth angle obtained by the radar, respectively, and θ is the elevation angle, and θ~U(0,2π).
[0130] Therefore, from the factor node Pass to variable node The message is The particle form is:
[0131]
[0132] in,
[0133] and
[0134]
[0135] In equation (29), and This represents the probability distribution function of the observation noise. Because... The particles are directly sampled from equation (28), therefore, the weights of the S particles are equal, all being 1 / S.
[0136] From variable node Passed to factor node f k The message can be written as:
[0137] For this form of multiple message product, particle sampling is performed using kernel density estimation.
[0138] Let p(x) represent the product of D input messages.
[0139]
[0140] To simplify the description, the representation of time variable k and node j is removed from equation (31). Based on the transmission of factor node messages above, the message m of the i-th input... i (x) can be represented as Therefore, the input message m i (x) can be represented nonparametrically using kernel density estimation:
[0141]
[0142] Wherein, N(x; x i,s ,Λ i Let x represent the normalized Gaussian density function, with mean and covariance x. i,s and Λ i In equation (32), the S Gaussian density functions choose the same covariance.
[0143]
[0144] in,
[0145]
[0146] Therefore, the product of the D input messages can be expressed as:
[0147]
[0148] It is easy to see that the product of D Gaussian density functions is still a Gaussian density function.
[0149]
[0150] in
[0151]
[0152] Therefore, the Gaussian density function shown in equation (34) weight It can be represented as:
[0153]
[0154] in, The weights of the Gaussian components in the Gaussian mixture shown by the D input messages are denoted as .
[0155] From equation (33), we can see that Let represent the product of D Gaussian mixtures, where each Gaussian mixture contains N Gaussian components. It contains a total of N D The computational complexity of directly sampling Gaussian components is O(N). D Therefore, importance sampling is used for the product of D Gaussian mixtures, and auxiliary variables are introduced to identify the Gaussian components of the sampled data.
[0156] θ i ∈{1,2,…,N} i=1,…,N (37), and θ 1:D ={θ1,θ2,…,θ D Discrete random variable θ i =l i This indicates the l-th term in the i-th Gaussian mixture. i One sample is sampled from each Gaussian component. Then equation (33) can be written as:
[0157]
[0158] The suggested distribution for importance sampling is p(θ). 1:D),and
[0159]
[0160] From the advisory distribution p(θ) 1:D Obtain the auxiliary variable θ 1:D ∈{1,2,…,N} D A set of samples l 1:D Then, the D weighted Gaussian density function products shown in equation (34) are calculated using equations (35) and (36). mean covariance and weight Then from the Gaussian density function Sample a particle Calculate weights Importance weights Normalization is performed to make the sum of the weights equal to 1. Therefore, the particle form of the product of D input messages is: It is worth noting that the computational complexity of this sampling method is only O(DN).
[0161] From variable node Passed to factor node g ji The message can be written as:
[0162] It is also in the form of a product of multiple messages, and particle sampling is performed using the kernel density estimation method described above.
[0163] It is worth noting the particle-based messages When entering a message Just about The function, therefore the variable node Divided into the first three items and the last three items Then, respectively, the messages... and Particle-based processing, and finally the particles and particles Merge to form the target state posterior probability distribution information Particleization method and same.
[0164] The specific algorithm flow for this step is shown in Table 2:
[0165]
[0166]
[0167]
[0168] In some embodiments, step 3 is specifically performed as follows:
[0169] The target's position and velocity vector in the fourth geocentric equatorial coordinate system is obtained using the nonparametric confidence propagation algorithm in step 2. Then, based on this position vector [x, y, z], the target's height in the geodetic coordinate system is calculated using an iterative approximation method. Considering the Earth's curvature, its equatorial radius is R = 6378.137 km, its polar radius is r = 6356.752 km, and f = (Rr) / r. The iteration accuracy is δ, and the maximum number of iterations is K.
[0170] First calculate z0 = z;
[0171] Step 2: Calculate A = rz0, B = 2Rx0, C = 2(R 2 -r 2 );
[0172] Step 3: Calculate the initial iteration value when K=0
[0173] Step 4: Calculation
[0174] Step 5: Calculate the iteration value t when K = k k =t k-1 -f(t k-1 ) / f′(t k-1 );
[0175] Step 6: When |t k -t k-1 |≤δ, or after reaching the maximum number of iterations K, the final t is obtained;
[0176] The final target altitude is:
[0177]
[0178] in:
[0179]
[0180] Example
[0181] Aerial single target altitude estimation under multi-radar network
[0182] Experimental background:
[0183] Consider a target flying at a constant speed in a straight line at an altitude of 9 km, tracked by three space-based radars. The orbital parameters of the three satellites are shown in Table 3. Other parameters are: satellite attitude angles are all 0°; radar antenna array attitude angles are all 0°; and roll angle θ is 0°. c =30°, pitch angle Yaw angle φ c =0°; the standard deviation of radar measurement noise is the same, radial distance 0.09km, azimuth angle 0.03°.
[0184] Table 3. Number of orbital elements for the three satellites.
[0185]
[0186] Analysis of Experiment 1:
[0187] Four algorithms were employed to estimate the target height: Extended Kalman Filter (EKF), Distributed Consistent Extended Kalman Filter (DCEKF), the Nonparametric Belief Propagation (NBP) algorithm used in this invention, and the geometric method. The DCEKF algorithm is a message calculation method based on extended Kalman filtering on a factor graph. Its specific steps are as follows:
[0188] Based on formula (5) for messages By performing a linear approximation, we can see from formulas (4) and (5) that:
[0189]
[0190] in for Given the Jacobian matrix, then:
[0191] Variable Node To its neighbor factor node g ji The message transmission method is as follows:
[0192]
[0193] Factor node g ji To variable node The message is:
[0194] Substituting equations (18) and (45) into equation (44), we initialize... and An iterative equation can be obtained, then the l-th iteration is:
[0195] After L iterations, the node Accepting neighbor nodes and To calculate the mean of its own posterior marginal distribution With variance
[0196]
[0197] The specific algorithm flow of the DCEKF step is shown in Table 1:
[0198]
[0199]
[0200] In the parameters of the two factor graph-based DCEKF and NBP algorithms, the coupling parameter β = 500. In the NBP algorithm, the number of sampled particles S = 800 and the number of iterations L = 400. For the EKF, DCEKF, and NBP algorithms, the system noise covariance matrix of these three algorithms is:
[0201] Q = diag{10 -3 10 -3 10 -3 10 -3 / T,10 -3 / T,10 -3 / T}; T is the sampling time, T = 15s.
[0202] The simulation yields the position of the aerial target in the fourth geocentric orbital coordinate system. According to step S3, the position is transformed to the geodetic coordinate system to obtain the target altitude estimate.
[0203] The algorithm performance evaluation metric is the root mean square error of the altitude of the aerial target, i.e.:
[0204]
[0205] Where N represents the number of Monte Carlo simulations, taken as N = 100, H k and h k This represents the target's true altitude (9km) and estimated altitude at time k.
[0206] Simulation results are shown below Figure 4 , Figure 5-(a) and 5-(b) as well as Figure 6 .
[0207] like Figure 4 This is a comparison chart of the root mean square errors (RMSEs) of aerial target altitude estimation using four methods: EKF, DCEKF, NBP, and the geometric method. From... Figure 4As can be seen, the RMSE of target altitude estimation using the traditional centralized EKF algorithm is 0.5km. However, the space environment of space-based radar is complex, and the robustness and flexibility of the traditional centralized filter are poor. The RMSE of target altitude estimation using the geometric method is between 1.5km and 3km, and since it only uses measurement information, the error is relatively large. In contrast, the RMSEs of altitude estimation using the DCEKF algorithm and the distributed NBP algorithm are 0.85km and 0.55km, respectively. The NBP algorithm improves the altitude estimation performance by 35.3% compared to the DCEKF algorithm. Therefore, the NBP algorithm not only possesses the flexibility and robustness of distributed algorithms, but also achieves an estimation accuracy that is basically comparable to that of the centralized EKF algorithm.
[0208] Analysis of Experiment 2:
[0209] Using the nonparametric confidence propagation algorithm (NBP) from this invention, the target height estimation results under different particle types, different coupling parameters, and different measurement errors were compared. Figure 5-(a) and 5-(b) It can be seen that the root mean square error of target height estimation under different coupling parameters is as follows when the number of sampled particles is S=400 and S=800 in the NBP algorithm. The estimation results are listed in Table 4.
[0210] Table 4. Target height estimation results of the NBP algorithm under different parameters.
[0211] Particle number S = 400 Number of particles S = 800 Coupling parameter β = 50 0.08196km 0.6596km Coupling parameter β = 500 0.7511km 0.5506km Coupling parameter β = 5000 0.6389km 0.4661km
[0212] As shown in Table 4, the height estimation accuracy increases with the increase of the coupling parameter β and the number of sampling particles S. This is because as β increases, the correlation between adjacent variables becomes stronger, and the measurement information of neighboring nodes is utilized more fully, resulting in better final height estimation accuracy. However, according to step S2.2, when β→∞, the distribution of local variables is the same as that of global variables. Therefore, as β increases, the height estimation accuracy will not continue to improve but will approach a certain constant value. A larger number of particles results in a smaller approximation error for nonlinear distributions, leading to better final estimation accuracy. Of course, increasing the number of particles also increases the computational complexity and communication load of the algorithm.
[0213] Analysis of Experiment 3:
[0214] like Figure 6 This is a comparison of the root mean square error (RMSE) of target altitude estimation using the NBP algorithm under different measurement errors. The three radars have the same standard deviation and are set at the following values: radial distance 0.07 km, azimuth 0.02°; radial distance 0.09 km, azimuth 0.03°; radial distance 0.11 km, azimuth 0.04°. Figure 6 As shown, the closer the radar is to the target, the smaller the measurement error, and the better the final target height estimation accuracy.
[0215] To address the problems of excessive elevation angle error in target tracking by single-satellite radars, lack of aerial target altitude estimation capability, strong measurement nonlinearity of satellite-borne radars, and high sensor robustness requirements, this invention proposes a distributed satellite-borne radar network aerial target altitude estimation method based on factor graphs. This method establishes a satellite-borne radar target tracking model based on factor graphs by constructing an aerial target motion model and a satellite-borne radar measurement model (excluding elevation angle measurement). Considering the nonlinearity of radar measurements, a DCEKF algorithm based on factor graph message passing is derived. To address the large approximation error of extended Kalman filtering for highly nonlinear problems, a sampling-based nonparametric confidence propagation algorithm (NBP) is proposed to calculate message passing in the factor graph. Simulation results demonstrate that the NBP algorithm achieves better altitude estimation accuracy than DCEKF, EKF, and geometric methods. The influence of NBP algorithm parameters on estimation performance is analyzed: the larger the coupling parameters between local variables and the number of sampled particles, the higher the estimation accuracy of the NBP algorithm; simultaneously, the smaller the radar measurement error, the higher the target altitude estimation accuracy.
Claims
1. A factor graph based distributed spaceborne radar netted air target height estimation method, characterized in that, The method comprises the following steps: Step 1, in the fourth ecliptic coordinate system of the earth center, a motion model of an air target and a measurement model of a spaceborne radar are constructed, and a factor graph-based spaceborne radar target tracking model, i.e. a factor graph, is constructed by combining the motion model of the air target and the measurement model of the spaceborne radar; Step 2, based on the factor graph, a coupling function relationship between local variables of the air target is introduced to obtain a posteriori marginal distribution of the local variables; A non-parametric credibility propagation algorithm is used to calculate the posteriori marginal distribution of the local variables, and an estimated value of the local variables is obtained according to a calculation result of the posteriori marginal distribution; Step 3, according to the estimated value of the local variables, an iterative approximation method is used to obtain a target height of the air target in a geodetic coordinate system. The specific content of the step 1 is as follows: S1.1, a motion model of an air target is constructed in the fourth ecliptic coordinate system of the earth center: x k = Fx k-1 + ω k-1 , where x k is the position and velocity state vector of the target in the fourth ECEF coordinate system at time k, x k-1 is the position and velocity state vector of the target in the fourth ECEF coordinate system at time k-1, ω k-1 is the process noise, which is assumed to be F is a target state transition matrix: wherein is the matrix direct product operator symbol, I3is the three-dimensional identity matrix, and T is the sampling interval; S1.2, a measurement model of a spaceborne radar is constructed in the fourth ecliptic coordinate system of the earth center: wherein, is the position vector of the target in the nth radar array coordinate system, wherein g(·) is the transformation relationship from the geocentric fourth orbital coordinate system to the radar array coordinate system; v k ~ N(0, R), is the radial range and azimuth measurement error of the N radars; the radial range and azimuth measurement of the nth radar is then the measurement of the N radars at time k is S1.3, factor graph modeling is performed by combining the target motion model and the radar measurement model: For the air target motion model and the space-borne radar measurement model, the joint probability density function of variables and can be decomposed as: where p(x k ) is the state transition probability density function, p(y k-1 |x k ) is the likelihood probability density function, and p(x k ) is the prior probability density function. And then, a decomposition form of the joint probability density function is described according to the factor graph, i.e. a factor graph model is obtained.
2. The factor graph based distributed space-borne radar netted airborne target height estimation method of claim 1, wherein, The specific process of introducing the coupling function relationship between the local variables to obtain the posteriori marginal distribution of the local variables of the air target in the step 2 is as follows: In the factor graph, a coupling factor node g ji is introduced to represent the global variable x k with two replicated state variables and on neighboring nodes j and i, Then, the posterior marginal distribution of the local variable at time k is given by: wherein is a status transition message, is a measurement message, is a coupling message, N j is the set of all neighbor nodes of node j.
3. The factor graph based distributed space-borne radar netted airborne target height estimation method of claim 2, wherein, In step 2, the posterior marginal distribution of the local variable at time k is calculated according to the factor graph model The specific process is as follows: kth instant, the granular form of the state transition message is calculated and the granular form of the metrology message The coupling messages are calculated by message iteration In the first iteration, the message product of all the neighbor nodes of node j is calculated by kernel density estimation: The coupling distribution is updated according to the coupling distribution relationship After L iterations, the coupling message is obtained The posterior distribution is calculated according to the kernel density estimation method: Then, the local variables are obtained based on the posterior distribution. The estimated value.