A Probabilistic Description Method of Initial Orbit Determination Results Based on Minimum Radar Admissible Region

Abstract

The invention discloses a probability description method for initial orbit determination results based on a minimum radar permissible domain, which solves the problem of how to describe the initial orbit determination results of the permissible domain method in a probabilistic form when no prior information is needed; it belongs to the field of satellite orbit control; it includes: step 1, calculating the weight coefficient of each sampling point within the minimum radar permissible domain by Bayesian formula; step 2, using a kernel density estimation method, combining the weight coefficient of each sampling point within the minimum radar permissible domain, and obtaining a minimum radar permissible domain probability density function that describes the initial orbit determination results of the minimum radar permissible domain in a probabilistic form on a radar permissible domain plane composed of right ascension change rate and declination change rate. The invention does not require prior information, limits the real orbit state within the minimum radar permissible domain, and describes the initial orbit determination results in a probabilistic form.

Description

Technical Field

[0001] The invention belongs to the technical field of satellite orbit control in the field of high-end equipment manufacturing, and relates to a probability description method for initial orbit determination results based on a minimum radar allowable domain. Background Art

[0002] Using ground-based radar to observe and monitor space targets is an important part of space situational awareness activities. Using radar observation data to determine a reliable orbit solution for space targets is a challenging task. For the initial state detection of a dynamic system, there are usually multiple solutions that can satisfy the state reflected by the observation data. This problem also exists in the above tasks.

[0003] The existing initial orbit determination algorithm takes the single minimum residual root mean square as the solution target, but in the real situation, that is, when there are measurement errors, the orbit state corresponding to this solution may not be close to the real orbit. Therefore, the results of the initial orbit determination should be described probabilistically within a reliable range to best reflect the actual situation. The orbit determination results described probabilistically will be conducive to the calculation of the correlation and collision probability between different radar observation arcs, and can also provide more reasonable initial values ​​for subsequent orbit improvements. In this process, how to determine a reliable parameter estimation range and perform uncertainty analysis within this range becomes the key to solving the problem.

[0004] For initial orbit determination and arc segment association in situational awareness missions, a popular method is the Admissible Region (AR) method, which applies the astrodynamic constraints of space targets to limit the parameters that ground equipment cannot observe, and then combines the observable parameters to estimate the possible orbital state. By finding the orbit that best matches the actual observation value within the admissible region, the dimension of the initial orbit determination problem can be reduced from six dimensions to two dimensions, which can greatly reduce the difficulty of orbit determination.

[0005] The original admissible region only considers energy constraints, which includes a large number of orbits with different semi-major axes and eccentricities. Therefore, the original admissible region contains a wide range of orbital states. Previous scholars added additional constraints to the admissible region based on some prior information, such as semi-major axis, eccentricity, perigee, apogee height, etc., to form a constrained admissible region (CAR), thereby narrowing the scope of the original admissible region. For example, DeMars et al. in DeMars, KJ and M.K. Jah, Probabilistic Initial Orbit Determination Using Gaussian Mixture Models. Journal of Guidance, Control, and Dynamics, 2013. 36 (5): p. 1324-1335. used a mixed Gaussian distribution to approximate the uniform distribution in the admissible region, and obtained a probabilistic description of the admissible region, and then obtained a probabilistic description of the initial orbital state. Psiaki et al. in Psiaki, ML, RM Weisman and MK Jah, Gaussian Mixture Approximation of Angles-Only Initial Orbit Determination Likelihood Function. Journal of Guidance, Control, and Dynamics, 2017. 40 (11): p. 2807-2819. further developed the method, using a new method to estimate Gaussian mixture components, and using the position / velocity state two-dimensional subspace instead of the distance-distance rate of change space to better describe the uncertainty of the orbit state. Hussein et al. in Islam I. Hussein, CWTR, Probabilistic admissible region for short-arcangles-only observations. 2014. formally proposed the concept of probabilistic admissible region (PAR), and proposed a method to construct the probabilistic admissible region based on the error characteristics of the observation data, and pointed out that the probabilistic admissible region is not uniformly distributed.Cai et al. used the outer probability measure theory in Cai, H., I. Hussein and M. Jah, Possibilistic admissible region using outer probability measure theory. Acta Astronautica, 2020. 177: p. 246-257. to obtain the uncertainty of the admissible region under the condition of unknown systematic errors and named it PAR+.

[0006] Existing research on the admissible domain theory is conducted in two directions: First, the original admissible domain is continuously reduced according to different constraints or by using other information. The most significant effect is: Publication No. CN115859013A; Patent name: A method and system for generating a minimum admissible threshold for optical observation arc segments, which provides a method for constructing a minimum optical admissible domain based on optical observation data. This method can greatly reduce the range of the optical admissible domain on the "distance, distance change rate" plane, and improve the solution efficiency and accuracy for subsequent arc segment associations. Second, based on different admissible domains, uncertainty analysis methods are used to describe the corresponding range probabilistically; the most typical one is the PAR+ method constructed by Cai et al. mentioned above, which is characterized by the fact that no prior information is required and it is applicable to both optical observation data and radar observation data.

[0007] Although the minimum admissible region theory can significantly reduce the original admissible region and limit the most likely orbital state to a very small region, it is only applicable to optical observation data and not to radar observation data due to the influence of different observation parameters.

[0008] Although the probability admissible domain can be used to perform uncertainty analysis on the basis of the constraint admissible domain, the constraint admissible domain can only reduce the admissible domain to a limited extent. Large-scale uncertainty analysis has no practical significance for the results of initial orbit determination, and cannot meet the technical requirements for subsequent radar observation arc segment association.

[0009] In addition, prior information, as a necessary condition for constructing the constraint admissible domain and the probability admissible domain, cannot be obtained in advance in some specific situation awareness tasks. Summary of the invention

[0010] In order to solve the technical problem of how to describe the initial orbit determination result of the admissible domain method in a probabilistic form when there is no prior information, the present invention provides a probabilistic description method for the initial orbit determination result based on the minimum radar admissible domain. The method calculates the weight coefficient of each sampling point within the minimum radar admissible domain by the Bayesian formula; the kernel density estimation method is used, combined with the weight coefficient of each sampling point within the minimum radar admissible domain, to obtain the minimum radar admissible domain probability density function that describes the initial orbit determination result of the minimum radar admissible domain in a probabilistic form on the radar admissible domain plane composed of the right ascension change rate and the declination change rate. In the minimum radar admissible domain, the present invention describes the result of the initial orbit determination in a probabilistic form, with the purpose of giving the initial orbit determination result in a more complete way and limiting the real orbit state to the minimum radar admissible domain as much as possible. The possibility of each sampling point in the minimum admissible domain becoming a real orbit is described by probability, which is convenient for subsequent arc segment association.

[0011] The purpose of the present invention is specifically achieved through the following technical solutions:

[0012] The present invention discloses a method for describing the probability of initial orbit determination results based on a minimum radar permissible region, the method comprising:

[0013] Step 1: Calculate the weight coefficient of each sampling point within the minimum radar allowable domain using the Bayesian formula;

[0014] Step 2: Using the kernel density estimation method and combining the weight coefficients of each sampling point within the minimum radar allowable domain, the minimum radar allowable domain probability density function is obtained on the radar allowable domain plane composed of the right ascension change rate and the declination change rate, which describes the initial orbit determination result of the minimum radar allowable domain in a probabilistic form.

[0015] In step 1, the method of using the Bayesian formula to calculate the weight coefficient of each sampling point within the minimum radar allowable domain includes:

[0016]

[0017] in, is the weight coefficient; η is the adjustment parameter; m is the total number of observation data; i is the sequence number of observation data;

[0018] is the actual observation distance; is the sampling observation distance; is the distance observation error; is the Jacobian matrix of distance propagation; Σ B is the uncertainty matrix;

[0019] is the actual observed right ascension; For sampling observation right ascension; is the right ascension observation error; is the Jacobian matrix of right ascension propagation;

[0020] is the actual observed declination; is the sampling observation declination; σ 2 is the declination observation error; is the Jacobian matrix of the declination propagation.

[0021] Furthermore, the uncertainty matrix Σ B for:

[0022]

[0023]

[0024] in, are the accuracy information of distance, distance change rate, right ascension and declination respectively; σ ρ is the radar ranging accuracy, σ is the radar angle measurement accuracy, T is the time span of the observation arc, and N is the number of observation points.

[0025] Furthermore, the Jacobian matrix of distance propagation is Jacobian matrix for right ascension propagation Jacobian matrix for declination propagation

[0026] in, is the propagation function of right ascension; is the propagation function of declination; is the propagation function of distance; is the actual arc attribute; B is the constructed arc attribute,

[0027] Where ρ is the distance, is the rate of change of distance, α is the right ascension, and δ is the declination; is a positive real number; is the total real number.

[0028] In step 2, the probability density function of the minimum radar allowable region is:

[0029]

[0030] in, is the probability density function of the minimum radar admissible region; K(·) is the kernel density function; is the radar allowable domain plane; is the rate of change of right ascension, is the rate of change of declination; is the i-th radar admissible domain plane; is the rate of change of the ith right ascension, is the i-th declination change rate; i is the sequence number of the observation data; h=(h1,h2,…,h d ) is a window width vector, d≧1, where h1 is the first window width; h2 is the second window width; is the i-th weight coefficient; e is the natural base.

[0031] In step 1, the method for constructing the minimum radar allowable domain range includes:

[0032] S1. Convert the received radar observation data from the station coordinate system to the geocentric inertial coordinate system through the conversion matrix to obtain the radar allowable domain arc segment attributes;

[0033] S2. Based on the radar admissible domain arc segment attributes, the original radar admissible domain is constructed on the radar admissible domain plane, and the least squares solution of the sampling points in the original radar admissible domain is obtained by using the Levenberg-Marquardt algorithm;

[0034] S3. Based on the least squares solution, a loss function is set to describe the standardized deviation of the theoretical observation value relative to the actual observation value. The loss function value of each sampling point in the original radar allowable domain is calculated through the loss function to determine the minimum radar allowable domain range.

[0035] The conversion matrix in S1 converts the azimuth and elevation angles in the received radar observation data into right ascension and declination. The conversion matrix is:

[0036]

[0037] in, is the radar allowed domain arc segment attribute, that is, the radar observation data in the geocentric inertial coordinate system, ρ is the distance, is the distance change rate, α is the right ascension, δ is the declination; λ S is the station longitude, Ω G It is the angle between the geocentric inertial coordinate system and the X-axis of the geocentric coordinate system, which can be calculated by looking up the astronomical almanac; is the latitude of the station, B rad It is the radar observation data in the station coordinate system, including distance, distance change rate, azimuth and elevation angle;

[0038] M2 is the second elementary transformation matrix,

[0039] M3 is the third elementary transformation matrix,

[0040] Among them, θ is the rotation angle, M2[θ] is the second elementary transformation matrix about the angle θ; M3[θ] is the third elementary transformation matrix about the angle θ.

[0041] In S2, based on the radar admissible domain arc segment attributes, the method of constructing the original radar admissible domain on the radar admissible domain plane includes:

[0042]

[0043] in, is the original radar allowed domain, is the radar allowable domain plane; is the rate of change of right ascension, is the rate of change of declination;

[0044]

[0045] Among them, Z is the coefficient matrix, which is a real symmetric matrix; 11 is the eleventh coefficient; z 12 is the twelfth coefficient; z 13 is the thirteenth coefficient; z 21 is the twenty-first coefficient; z 22 is the 22nd coefficient; z 23 is the twenty-third coefficient; z 31 is the thirty-first coefficient; z 32 is the thirty-second coefficient; z 33 is the thirty-third coefficient;

[0046] ρ is the distance between the spacecraft and the observation station in the radar observation data; δ is the declination; ω 3,α is the seventh sub-coefficient; ω 3,δ is the eighth sub-coefficient; is the distance change rate; ω1 is the second sub-coefficient; ω4 is the fifth sub-coefficient; μ E is the earth's gravitational constant; α is the right ascension; S(ρ) is the coefficient polynomial: S(ρ)=ρ 2 +ω5ρ+ω0; where ω0 is the first sub-coefficient; ω5 is the sixth sub-coefficient:

[0047] ω0=||R|| 2 ;

[0048]

[0049] ω5=2R·L;

[0050]

[0051] Where R is the position vector of the observation platform; is the velocity vector of the observation platform; L is the unit vector of the satellite relative to the observation platform; L α is the partial derivative of the satellite's unit vector relative to the observation platform with respect to right ascension; L δis the partial derivative of the satellite's unit vector with respect to the observation platform with respect to declination.

[0052] In S2, the method of using the Levenberg-Marquardt algorithm to obtain the least squares solution of the sampling points in the original radar admissible domain includes:

[0053] Using the least squares idea, uniform sampling is performed in the original radar admissible domain, and the first loss function of each sampling point is calculated. Use the Levenberg-Marquardt algorithm to find the first loss function The sampling point with the smallest value is the least square solution in the original radar admissible domain; among which, the first loss function is:

[0054]

[0055] in, is the first loss function; m is the total number of radar observation data; i is the sequence number of radar observation data;

[0056] is the actual observation distance; is the observation distance to be estimated; is the distance observation error;

[0057] is the actual observed right ascension; is the observed right ascension to be estimated; is the right ascension observation error;

[0058] is the actual observed declination; is the observed declination to be estimated; σ 2 is the declination observation error.

[0059] In S3, the specific steps include:

[0060] Based on the least squares solution, set the loss function to describe the standardized deviation of the theoretical observations from the actual observations:

[0061]

[0062] The loss function is used to calculate the loss function value of each sampling point in the original radar permissible domain, and the threshold β ≥ 3 is set; the loss function value that meets the loss function value is eliminated. The sampling points that meet The sampling points are included in the minimum radar allowable domain, so as to determine the range of the minimum radar allowable domain;

[0063] in, is the loss function; m is the total number of observations; i is the sequence number of the observations;

[0064] is the actual observation distance; is the theoretical observation distance; is the distance observation error; is the Jacobian matrix of distance propagation; Σ B is the uncertainty matrix;

[0065] is the actual observed right ascension; is the theoretical observed right ascension; is the right ascension observation error; is the Jacobian matrix of right ascension propagation;

[0066] is the actual observed declination; is the theoretical observed declination; σ 2 is the declination observation error; is the Jacobian matrix of the declination propagation.

[0067] The beneficial effects of the present invention are:

[0068] Using the Bayesian formula to calculate the weight coefficient of each sampling point within the minimum radar permissible domain can objectively reflect the probability of the track corresponding to each sampling point becoming the real track. The calculated probability is proportional to the weight coefficient, and the subsequent probability density estimation will be more accurate.

[0069] The probability density is estimated using a weighted kernel density estimation method, and the probability density function of the minimum radar allowable domain is obtained; compared with the more common parameter estimation method and semi-parameter estimation method, the weighted kernel density estimation method disclosed in the present invention can solve the strong nonlinear problem of initial orbit determination.

[0070] At the same time, the present invention also narrows the scope of uncertainty analysis from the original radar allowable domain to the minimum radar allowable domain. The minimum radar allowable domain greatly narrows the scope of the original radar allowable domain, includes the least squares solution of the original radar allowable domain, and includes the true orbit state with a very high probability, and can more objectively and specifically describe the initial orbit determination result probabilistically. The minimum radar allowable domain can include the true orbit state with a probability greater than 99.7%, which means that the uncertainty analysis under the minimum radar allowable domain has more practical application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] The present invention is further described in detail below based on the accompanying drawings and embodiments.

[0072] Figure 1 This is a schematic diagram of the location of the least squares solution in the original radar admissible domain.

[0073] Figure 2It is a schematic diagram of the positional relationship between the minimum radar allowable domain range and the true value of the least squares solution.

[0074] Figure 3 It is a schematic diagram of the weight coefficients of each sampling point.

[0075] Figure 4 It is a three-dimensional schematic diagram of the probability density distribution within the minimum radar allowable domain.

[0076] Figure 5 It is a two-dimensional schematic diagram of the probability density distribution within the minimum radar allowable domain. DETAILED DESCRIPTION

[0077] Embodiment 1

[0078] Embodiment 1 of the present invention provides a method for describing the probability of an initial orbit determination result based on a minimum radar permissible region, the method comprising:

[0079] Step 1: Calculate the weight coefficient of each sampling point within the minimum radar allowable domain using the Bayesian formula;

[0080] Step 2: Using the kernel density estimation method and combining the weight coefficients of each sampling point within the minimum radar allowable domain, the minimum radar allowable domain probability density function is obtained on the radar allowable domain plane composed of the right ascension change rate and the declination change rate, which describes the initial orbit determination result of the minimum radar allowable domain in a probabilistic form.

[0081] In step 1, the method of using the Bayesian formula to calculate the weight coefficient of each sampling point within the minimum radar allowable domain includes:

[0082]

[0083] in, is the weight coefficient; η is the adjustment parameter; m is the total number of observation data; i is the sequence number of observation data;

[0084] is the actual observation distance; is the sampling observation distance; is the distance observation error; is the Jacobian matrix of distance propagation; Σ B is the uncertainty matrix;

[0085] is the actual observed right ascension; For sampling observation right ascension; is the right ascension observation error; is the Jacobian matrix of right ascension propagation;

[0086] is the actual observed declination; is the sampling observation declination; σ2 is the declination observation error; is the Jacobian matrix of the declination propagation.

[0087] Among them, the minimum radar allowable domain of step one is constructed on An irregular area on a plane, the least squares solution is at the geometric center of the area. Every point in the minimum radar allowable domain may become a real orbit state, but the probability of each point is different. In order to accurately describe the probability of each point becoming a real orbit, it is also necessary to estimate the probability density of the minimum radar allowable domain.

[0088] To estimate the probability density, we must first determine the weight of each sampling point within the minimum radar allowable domain. They all correspond to a set of orbital parameters x0. Under the assumption of Gaussian noise, the probability density of x0 can be calculated using the Bayesian formula:

[0089]

[0090] The present invention solves the problem of describing the probability of the initial orbit determination result without prior information, which means that x0 obeys a uniform distribution. The expression of

[0091]

[0092] The above formula expresses a positive proportional relationship and defines the weight coefficient. Among them, η is the adjustment parameter.

[0093]

[0094] Mapped on the radar allowed domain plane, we have

[0095]

[0096] This coefficient is proportional to the probability density and can be used to characterize the probability that the orbit corresponding to the sampling point becomes the real orbit. From the weight coefficient expression, it can be seen that the least squares orbit is also the orbit with the largest weight.

[0097] Uncertainty matrixΣ B for:

[0098]

[0099]

[0100] in, are the accuracy information of distance, distance change rate, right ascension and declination respectively; σ ρ is the radar ranging accuracy, σ is the radar angle measurement accuracy, T is the time span of the observation arc, and N is the number of observation points.

[0101] Jacobian matrix for distance propagation Jacobian matrix for right ascension propagation Jacobian matrix for declination propagation

[0102] in, is the propagation function of right ascension; is the propagation function of declination; is the propagation function of distance; is the actual arc attribute; B is the constructed arc attribute,

[0103] Where ρ is the distance, is the rate of change of distance, α is the right ascension, and δ is the declination; is a positive real number; is the total real number.

[0104] In step 2, the probability density function of the minimum radar allowable region is:

[0105]

[0106] in, is the probability density function of the minimum radar admissible region; K(·) is the kernel density function; is the radar allowable domain plane; is the rate of change of right ascension, is the rate of change of declination; is the i-th radar admissible domain plane; is the rate of change of the ith right ascension, is the i-th declination change rate; i is the sequence number of the observation data; h=(h1,h2,…,h d ) is a window width vector, d≥1, where h1 is the first window width; h2 is the second window width; is the i-th weight coefficient; e is the natural base.

[0107] There are many methods for estimating the probability density function of the population through samples, which are usually divided into parameter estimation methods, semi-parametric estimation methods and non-parametric estimation methods. The parameter estimation method is to assume that the variable obeys a certain form of probability distribution (such as normal distribution), and estimate a small number of parameters that can fully describe the distribution based on the sample; the non-parametric estimation method does not need to assume the probability distribution form obeyed by the population, and can directly give the probability density function based on the sample; the semi-parametric estimation method takes into account the uncertainty of the model and is a combination of parameter estimation and non-parametric estimation.

[0108] Due to the strong nonlinear characteristics of the initial orbit determination problem, the orbit determination results do not necessarily obey the normal distribution or other specific forms of distribution. Therefore, both the parameter estimation method and the semi-parametric estimation method are no longer applicable. Kernel density estimation is one of the non-parametric estimation methods, also known as the Parzen window, which is currently the most effective and widely used non-parametric probability density estimation algorithm. The present invention uses the kernel density method to estimate the probability density function, and the key lies in the selection of the kernel function and the determination of the window width.

[0109] Let x1, x2, …, x n is a sample obtained by sampling a one-dimensional random variable x (independent and identically distributed), then the probability density function f(x) of x can be estimated as:

[0110]

[0111] Where K(·) is the kernel function, including Gaussian function, Epanechnikov function, Biweight function, etc.; h is the window width. When the weights of the variables are different, 1 / n needs to be corrected to the weight value w of each variable. i ,Right now

[0112]

[0113] For a d-dimensional random variable X=(x1,x2,…,x d ) T , let the sample be X i =(x i1 ,x i2 ,…x id ) T (i=1,2,…n), let h=(h1,h2,…,h d ) is a window width vector, then the kernel density estimate of X is:

[0114]

[0115] The present invention selects Gaussian kernel function as kernel density function, and the expression is

[0116]

[0117] The above formula is the estimated probability density function of the minimum radar allowable region. At this point, the initial orbit determination result within the minimum radar allowable region can be described in the form of probability.

[0118] In step 1, the method for constructing the minimum radar allowable domain range includes:

[0119] S1. Convert the received radar observation data from the station coordinate system to the geocentric inertial coordinate system through the conversion matrix to obtain the radar allowable domain arc segment attributes;

[0120] S2. Based on the radar admissible domain arc segment attributes, the original radar admissible domain is constructed on the radar admissible domain plane, and the least squares solution of the sampling points in the original radar admissible domain is obtained by using the Levenberg-Marquardt algorithm;

[0121] S3. Based on the least squares solution, a loss function is set to describe the standardized deviation of the theoretical observation value relative to the actual observation value. The loss function value of each sampling point in the original radar allowable domain is calculated through the loss function to determine the minimum radar allowable domain range.

[0122] The conversion matrix in S1 converts the azimuth and elevation angles in the received radar observation data into right ascension and declination. The conversion matrix is:

[0123]

[0124] in, is the radar allowed domain arc segment attribute, that is, the radar observation data in the geocentric inertial coordinate system, ρ is the distance, is the distance change rate, α is the right ascension, δ is the declination; λ S is the station longitude, Ω G It is the angle between the geocentric inertial coordinate system and the X-axis of the geocentric coordinate system, which can be calculated by looking up the astronomical almanac; is the latitude of the station, B rad It is the radar observation data in the station coordinate system, including distance, distance change rate, azimuth and elevation angle;

[0125] M2 is the second elementary transformation matrix,

[0126] M3 is the third elementary transformation matrix,

[0127] Among them, θ is the rotation angle, M2[θ] is the second elementary transformation matrix about the angle θ; M3[θ] is the third elementary transformation matrix about the angle θ.

[0128] In order to facilitate the description and calculation of the orbital state of space targets in the observation space, the radar allowable domain is usually constructed on the (right ascension rate, declination rate) plane. To facilitate subsequent calculations, it is necessary to map the actual observations to the radar allowable domain plane, that is, to convert the radar observation data (including: distance, distance change rate, azimuth, pitch angle) into radar allowable domain arc attributes (including: distance, distance change rate, right ascension, declination).

[0129] In S2, based on the radar admissible domain arc segment attributes, the method of constructing the original radar admissible domain on the radar admissible domain plane includes:

[0130]

[0131] in, is the original radar allowed domain, is the radar allowable domain plane; is the rate of change of right ascension, is the rate of change of declination;

[0132]

[0133] Among them, Z is the coefficient matrix, which is a real symmetric matrix; 11 is the eleventh coefficient; z 12 is the twelfth coefficient; z 13 is the thirteenth coefficient; z 21 is the twenty-first coefficient; z 22 is the 22nd coefficient; z 23 is the twenty-third coefficient; z 31 is the thirty-first coefficient; z 32 is the thirty-second coefficient; z 33 is the thirty-third coefficient;

[0134] ρ is the distance between the spacecraft and the observation station in the radar observation data; δ is the declination; ω 3,α is the seventh sub-coefficient; ω 3,δ is the eighth sub-coefficient; is the distance change rate; ω1 is the second sub-coefficient; ω4 is the fifth sub-coefficient; μ E is the earth's gravitational constant; α is the right ascension; S(ρ) is the coefficient polynomial: S(ρ)=ρ 2 +ω5ρ+ω0; where ω0 is the first sub-coefficient; ω5 is the sixth sub-coefficient:

[0135] ω0=||R|| 2 ;

[0136]

[0137] ω5=2R·L;

[0138]

[0139] Where R is the position vector of the observation platform; is the velocity vector of the observation platform; L is the unit vector of the satellite relative to the observation platform; L α is the partial derivative of the satellite's unit vector relative to the observation platform with respect to right ascension; L δ is the partial derivative of the satellite's unit vector with respect to the observation platform with respect to declination.

[0140] Among them, S1 constructs the radar allowable domain arc attribute When constructing the original radar allowed domain, only the energy constraint of the orbit needs to be considered. For near-Earth space targets, the spacecraft runs on an elliptical orbit, and the orbit energy is less than 0.

[0141]

[0142] Observe the above formula, the key to solving the problem is to find the number of orbital roots in the allowed domain and Cartesian orbital parameters Analyzing the positional relationship between the observation radar equipment, space targets and the center of the earth, we have

[0143]

[0144] In the formula, R, are the position and velocity vector of the observation platform respectively, where

[0145]

[0146] The orbital energy available is about The expression is:

[0147]

[0148] in,

[0149]

[0150] The orbital energy of a space target can be expressed as

[0151]

[0152] So, in The original radar allowable domain range on the plane is

[0153]

[0154] In S2, the method of using the Levenberg-Marquardt algorithm to obtain the least squares solution of the sampling points in the original radar admissible domain includes:

[0155] Using the least squares idea, uniform sampling is performed in the original radar admissible domain, and the first loss function of each sampling point is calculated. Use the Levenberg-Marquardt algorithm to find the first loss function The sampling point with the smallest value is the least square solution in the original radar admissible domain; among which, the first loss function is:

[0156]

[0157] in, is the first loss function; m is the total number of radar observation data; i is the sequence number of radar observation data;

[0158] is the actual observation distance; is the observation distance to be estimated; is the distance observation error;

[0159] is the actual observed right ascension; is the observed right ascension to be estimated; is the right ascension observation error;

[0160] is the actual observed declination; is the observed declination to be estimated; σ 2 is the declination observation error.

[0161] Among them, the first loss function The fitting effect of the orbital state represented by the sampling points and the radar observation data is characterized. The smaller the value of means the better the fit between the orbit state and the radar observation data. The Levenberg-Marquardt algorithm is used to find the sampling point with the best fit, which is the least squares solution within the original radar admissible domain.

[0162] In S3, the specific steps include:

[0163] Based on the least squares solution, set the loss function to describe the standardized deviation of the theoretical observations from the actual observations:

[0164]

[0165] The loss function is used to calculate the loss function value of each sampling point in the original radar permissible domain, and the threshold β ≥ 3 is set; the loss function value that meets the loss function value is eliminated. The sampling points that meet The sampling points are included in the minimum radar allowable domain, so as to determine the range of the minimum radar allowable domain;

[0166] in, is the loss function; m is the total number of observations; i is the sequence number of the observations;

[0167] is the actual observation distance; is the theoretical observation distance; is the distance observation error; is the Jacobian matrix of distance propagation; Σ B is the uncertainty matrix;

[0168] is the actual observed right ascension; is the theoretical observed right ascension; is the right ascension observation error; is the Jacobian matrix of right ascension propagation;

[0169] is the actual observed declination; is the theoretical observed declination; σ 2 is the declination observation error; is the Jacobian matrix of the declination propagation.

[0170] Among them, the sources of uncertainty that lead to the radar allowable domain can be divided into accidental uncertainty and epistemic uncertainty. Accidental uncertainty refers to the inherent randomness of the system, which can be characterized by probability distribution, such as the distance observation error σ of the observation equipment ρ and angle observation error σ; cognitive uncertainty refers to the information that causes errors under the influence of the system, which is not a random phenomenon. In the radar allowable domain, the latter can be analyzed by deriving the propagation law of the uncertainty of the arc segment attributes to the observation space. Therefore, the theoretical observation value errors of distance, right ascension, and declination are jointly affected by two parts: the first part is the observation error of the sensor σ 2 , the second part is the allowable domain uncertainty error have

[0171]

[0172] According to the characteristics of the standard normal distribution, we have

[0173]

[0174] The above analysis will be used to design the minimum allowable domain boundary based on radar observation data. The purpose of constructing the minimum radar allowable domain is to reduce its scope based on the original radar allowable domain, but to ensure that the true orbit state can appear in the minimum radar allowable domain with the greatest possible probability. The loss function is defined to describe the standardized deviation of the theoretical observation value relative to the actual observation value (absolute error divided by the standard deviation).

[0175]

[0176] In the formula, JMinAR represents the average angular deviation of the track. It can be seen that the loss function strictly conforms to the normal distribution. Referring to its 3σ criterion, the threshold is set to β. For any sampling point in the original radar allowable domain Determine the loss function Is the value of less than β? The corresponding sampling point is considered to be very likely to be the real orbit state and is included in the minimum radar allowable domain. It is not included in the minimum radar tolerance region. When setting the parameter β, let β≥3. In this way, the actual orbit state corresponds to The probability of being included in the minimum radar tolerance domain is greater than 99.7%. The range setting method of the minimum radar tolerance domain excludes sampling points with large average angle deviation, which is also in line with objective laws.

[0177] Verify the instance

[0178] The embodiment of the present invention takes a simulated radar observation arc of a LEO target as an example to verify the above scheme. Table 1 shows the specific information and observation error information of the selected space target.

[0179] Table 1

[0180] Target Information Numeric Epoch time (UTCG) 20Mar 2020 15:04:37.000 a / km 6923.161730 e 0.000832 i / ° 53.149 Ω / ° 341.470 w / ° 154.344 f / ° 346.534 Observation arc duration / s 61 Angle observation error / ° 0.1 Distance observation error / m 30

[0181] Based on the information disclosed in Table 1, a probability description of the initial orbit determination result based on the minimum radar allowable region is given.

[0182] The radar observation data (range, range change rate, azimuth, pitch angle) are converted into radar allowable domain arc attributes (range, range change rate, right ascension, declination); the original radar allowable domain is constructed on the (right ascension change rate, declination change rate) plane, and the least squares solution is obtained using the Levenberg-Marquardt algorithm. Figure 1 Indicates the position of the least squares solution (the "*" shape in the figure) in the original radar admissible domain.

[0183] Based on the least squares solution obtained, the uncertainty propagation law of the accuracy information contained in the observation data is derived, the loss function is set and the loss function value of each sampling point is calculated to determine the minimum radar allowable domain range. Figure 2 It represents the range of the minimum radar admissible domain and the positional relationship between the least squares solution and the true value.

[0184] Precise sampling is performed in the minimum radar allowable domain, and the weight coefficient of each sampling point is calculated using the Bayesian formula. Figure 3 Indicates the weight coefficient of each sampling point.

[0185] Combined with the weight coefficient of each sampling point, the weighted two-dimensional kernel density estimation method is used to finally obtain the probability density function of the minimum radar admissible domain. Figure 4 and Figure 5 are the probability density functions obtained after estimation. Figure 4 is the three-dimensional probability density distribution, Figure 5The two-dimensional probability density distribution of . It can be seen from the two figures that the probability of the initial orbit determination result becoming the true orbit state is mainly concentrated around the least squares solution. Usually, the least squares solution has the greatest probability of becoming the true orbit state, which is also in line with the normal uncertainty analysis logic.

[0186] The beneficial effects of the present invention are:

[0187] Using the Bayesian formula to calculate the weight coefficient of each sampling point within the minimum radar permissible domain can objectively reflect the probability of the track corresponding to each sampling point becoming the real track. The calculated probability is proportional to the weight coefficient, and the subsequent probability density estimation will be more accurate.

[0188] The probability density is estimated using a weighted kernel density estimation method, and the probability density function of the minimum radar allowable domain is obtained; compared with the more common parameter estimation method and semi-parameter estimation method, the weighted kernel density estimation method disclosed in the present invention can solve the strong nonlinear problem of initial orbit determination.

[0189] At the same time, the present invention also narrows the scope of uncertainty analysis from the original radar allowable domain to the minimum radar allowable domain. The minimum radar allowable domain greatly narrows the scope of the original radar allowable domain, includes the least squares solution of the original radar allowable domain, and includes the true orbit state with a very high probability, and can more objectively and specifically describe the initial orbit determination result probabilistically. The minimum radar allowable domain can include the true orbit state with a probability greater than 99.7%, which means that the uncertainty analysis under the minimum radar allowable domain has more practical application value.

[0190] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art can easily think of changes or substitutions within the technical scope disclosed by the present invention, which should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. A method for describing the probability of initial orbit determination results based on the minimum radar allowable region, characterized in that: The method includes: Step 1: Calculate the weight coefficient of each sampling point within the minimum radar allowable domain using the Bayesian formula; Step 2: Using the kernel density estimation method and combining the weight coefficients of each sampling point within the minimum radar allowable domain, a probability density function of the minimum radar allowable domain is obtained on the radar allowable domain plane composed of the right ascension change rate and the declination change rate, which describes the initial orbit determination result of the minimum radar allowable domain in a probabilistic form; In step 1, the method of using the Bayesian formula to calculate the weight coefficient of each sampling point within the minimum radar allowable domain includes: ;in, is the weight coefficient; To adjust the parameters; m is the total number of observations; i is the sequence number of the observation data; is the actual observation distance; is the sampling observation distance; is the distance observation error; is the Jacobian matrix of distance propagation; is the uncertainty matrix; is the actual observed right ascension; For sampling observation right ascension; is the right ascension observation error; is the Jacobian matrix of right ascension propagation; is the actual observed declination; To sample the observed declination; is the declination observation error; is the Jacobian matrix of declination propagation; The uncertainty matrix for: ； ； in, They are the accuracy information of distance, distance change rate, right ascension and declination; is the radar ranging accuracy, is the radar angle measurement accuracy, T is the time span of the observation arc, N is the number of observation points; In step 2, the probability density function of the minimum radar allowable region is: ； in, is the probability density function of the minimum radar admissible domain; K (·) is the kernel density function; is the radar allowable domain plane; is the rate of change of right ascension, is the rate of change of declination; For the i A radar admissible domain plane; For the i The rate of change of right ascension, For the i Declination change rate; i is the sequence number of the observation data; is a window width vector, d ≧1, where is the first window width; is the second window width; For the i Weight coefficients; Is the natural base.

2. The method according to claim 1, characterized in that Jacobian matrix for distance propagation ; Jacobian matrix of right ascension propagation ; Jacobian matrix for declination propagation ; in, is the propagation function of right ascension; is the propagation function of declination; is the propagation function of distance; is the actual arc attribute; is the constructed arc segment attribute, ; in, For distance, is the distance change rate, is right ascension, is declination; is a positive real number; is the total real number.

3. The method according to claim 1, characterized in that In step 1, the method for constructing the minimum radar allowable domain range includes: S1. Convert the received radar observation data from the station coordinate system to the geocentric inertial coordinate system through the conversion matrix to obtain the radar allowable domain arc segment attributes; S2. Based on the radar admissible domain arc segment attributes, the original radar admissible domain is constructed on the radar admissible domain plane, and the least squares solution of the sampling points in the original radar admissible domain is obtained by using the Levenberg-Marquardt algorithm; S3. Based on the least squares solution, a loss function is set to describe the standardized deviation of the theoretical observation value relative to the actual observation value. The loss function value of each sampling point in the original radar allowable domain is calculated through the loss function to determine the minimum radar allowable domain range.

4. The method according to claim 3, characterized in that The conversion matrix in S1 converts the azimuth and elevation angles in the received radar observation data into right ascension and declination. The conversion matrix is: ； in, is the radar allowed domain arc segment attribute, that is, the radar observation data in the geocentric inertial coordinate system, , For distance, is the distance change rate, is right ascension, is declination; is the station longitude, It is the angle between the geocentric inertial coordinate system and the X-axis of the geocentric coordinate system, which is calculated by looking up the astronomical almanac; is the station latitude, It is the radar observation data in the station coordinate system, including distance, distance change rate, azimuth and elevation angle; is the second elementary transformation matrix, ; is the third elementary transformation matrix, ; in, θ is the rotation angle, M 2[ θ ] is about angle θ The second elementary transformation matrix of ; M 3[ θ ] is about angle θ The third elementary transformation matrix.

5. The method according to claim 3, characterized in that In S2, based on the radar admissible domain arc segment attributes, the method of constructing the original radar admissible domain on the radar admissible domain plane includes: ； ； in, is the original radar allowed domain, is the radar allowable domain plane; is the rate of change of right ascension, is the rate of change of declination; ； Among them, Z is the coefficient matrix, which is a real symmetric matrix; is the eleventh coefficient; is the twelfth coefficient; is the thirteenth coefficient; is the twenty-first coefficient; is the twenty-second coefficient; is the twenty-third coefficient; is the thirty-first coefficient; is the thirty-second coefficient; is the thirty-third coefficient; is the distance between the spacecraft and the observation station in the radar observation data; is declination; is the seventh sub-coefficient; is the eighth sub-coefficient; is the distance change rate; is the second sub-coefficient; is the fifth sub-coefficient; is the Earth's gravitational constant; is right ascension; The coefficients are polynomials: ;in, is the first sub-coefficient; is the sixth sub-coefficient: in, is the position vector of the observation platform; is the velocity vector of the observation platform; is the unit vector of the satellite relative to the observation platform; is the partial derivative of the satellite's unit vector relative to the observation platform with respect to right ascension; is the partial derivative of the satellite's unit vector with respect to the observation platform with respect to declination.

6. The method according to claim 3, characterized in that In S2, the method of using the Levenberg-Marquardt algorithm to obtain the least squares solution of the sampling points in the original radar admissible domain includes: Using the least squares idea, uniform sampling is performed in the original radar admissible domain, and the first loss function of each sampling point is calculated. , use the Levenberg-Marquardt algorithm to find the first loss function The sampling point with the smallest value is the least square solution in the original radar admissible domain; among which, the first loss function is: ； in, is the first loss function; m is the total number of radar observation data; i is the serial number of radar observation data; is the actual observation distance; is the observation distance to be estimated; is the distance observation error; is the actual observed right ascension; is the observed right ascension to be estimated; is the right ascension observation error; is the actual observed declination; is the observed declination to be estimated; is the declination observation error.

7. The method according to claim 3, characterized in that In S3, the specific steps include: Based on the least squares solution, set the loss function to describe the standardized deviation of the theoretical observations from the actual observations: ； The loss function is used to calculate the loss function value of each sampling point in the original radar permissible domain and set the threshold ; Eliminate sampling points that meet the loss function value , only the The sampling points are included in the minimum radar allowable domain, so as to determine the range of the minimum radar allowable domain; in, is the loss function; m is the total number of observations; i is the sequence number of the observation data; is the actual observation distance; is the theoretical observation distance; is the distance observation error; is the Jacobian matrix of distance propagation; is the uncertainty matrix; is the actual observed right ascension; is the theoretical observed right ascension; is the right ascension observation error; is the Jacobian matrix of right ascension propagation; is the actual observed declination; is the theoretical observed declination; is the declination observation error; is the Jacobian matrix of the declination propagation.

Images