# Estimation and Analysis of Powered Phase’s Orbit of the Spacecraft Based on Multiple-Point Intersection Algorithm

Authors: Honghai Ji; Yuzhou Wei; Lingling Fan; Shida Liu; Zijun Jia; Li Wang
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

To estimate and analyze the orbit of a spacecraft in the ascent phase, three algorithms are proposed in this paper: a two-satellite point-by-point intersection positioning algorithm, single-satellite positioning algorithm, and multiple-point intersection collaborative positioning algorithm. Meanwhile, two mathematic models of the spacecraft orbit are presented based on the kinematics characteristic and the least square method. Besides, by combining the results of the proposed satellite positioning algorithms with mathematic models of the spacecraft orbit, the particle filtering based on the Monte Carlo algorithm is applied to filtering of the final results of the spacecraft. Finally, the effectiveness of the two-satellite point-by-point intersection positioning algorithm and the single-satellite positioning algorithm is verified through the experiment simulation on the real data.


---

## Body

## 1. Introduction

Some countries launch spacecraft with special purposes, such as ballistic missiles and reconnaissance satellites. Therefore, it is of great strategic significance to monitor and quickly respond to hostile spacecraft launched by other countries to maintain national security [1]. Discovery of launch and detection of its orbital parameters is the first step to realize monitoring and responding. Currently, satellites are an important platform for detecting the launch and orbit parameters of the spacecraft [2].The detectors installed on satellites can be divided into active and passive ones. Active detectors use active methods (such as radar and laser) to search for targets. Passive detectors passively receive target radiation. The purpose of detection is to infer the orbital parameters of the spacecraft. Specifically, inferences are made based on observational data and through mathematical models and calculation methods. When the observation satellite flies for a period of time, the detector measures the observation data of the target relative to the moving satellite. Based on the movement model and observation model of the observation satellite and the spacecraft, the mathematical inference is made on the orbital parameters of the spacecraft (including the orbit position, initial value of velocity, and other model parameters), providing information to judge the spacecraft types and flight intentions [3, 4].Generally, the orbit of a spacecraft is divided into three sections: the ascent phase propelled by rockets, the inertial flight period in the outer space of the Earth, and the attack period after reentry into the atmosphere. Since the motion process of the active section is more complicated than the other two sections, this paper focuses on the study of the orbit estimation in the ascent phase of the spacecraft.Currently, several studies have been conducted on certain aspects of satellite positioning, including the single-satellite direction of arrival (DOA) locating methods based on the WGS-84 model [5], single-satellite locating and tracking algorithms for synchronous orbiting satellites with bearing-only measurements [6], passive satellite locating and tracking methods by single satellite [7], multisatellite coordinated locating methods [8], single-satellite orbit determination methods for navigation satellites [9–11], and autonomous orbit determination method for navigation satellites [12]. Although there have been some research achievements, more research is based on single satellite orbit determination algorithm. In addition, in order to improve the accuracy of observation and positioning, filtering algorithm can be used for estimation. However, the classical Kalman filtering is only applicable to linear systems. In order to solve the estimation problem in nonlinear systems, particle filtering algorithm [13, 14] is selected to improve the accuracy of observation. Therefore, this paper investigates the problem of orbit estimation and error analysis of the ascent phase of the spacecraft based on passive detectors of satellites. Based on different numbers of observation satellites, three algorithms for spacecraft orbit estimation are proposed in this paper: two-satellite point-by-point intersection positioning algorithm, single-satellite positioning algorithm for estimating the orbit of a single spacecraft in the ascent phase, and multiple-point intersection collaborative positioning algorithm for estimating the orbit of multiple spacecraft in the ascent phase.Through the experimental analysis of real data, the effectiveness of the proposed two-satellite point-by-point intersection positioning algorithm and single-satellite positioning algorithm for estimating the orbit of a single spacecraft in the ascent phase is verified.The main contributions of this paper are as follows:(1)
To make the time stamps of each satellite sampling point consensus, cubic spline interpolation is utilized. Furthermore, Runge–Kutta method is used to solve the differential equations, and the trajectories of the satellites are obtained.(2)
The data observed by the two satellites is converted to the basic coordinate system through coordinate transformation. The coordinates of the observed target points are calculated and three-dimensional curves are drawn in the Earth coordinate system.(3)
Two trajectory estimation models of the target spacecraft are established, using the least square method to fit the model parameters. The particle filtering algorithm is used to improve the observation positioning accuracy.(4)
The nonlinear least square method is used to solve the systematic errors in the double satellite positioning. The single-star positioning scheme, multistar positioning error estimation method and target orbit estimation method are given.


## 2. Description of the Orbit of a Spacecraft in the Ascent Phase

According to [15, 16], the schematic diagram of the orbit of a spacecraft in the ascent phase (not in actual proportion) is shown in Figure 1. The ascent phase can be divided into several sections: the vertical ascending section, the program turning section, and the gravity oblique flight section. Supposing that point A is the ground launch point, the section of AB is the vertical ascending section; the BC arc is the program turning section; the CD arc is the gravity oblique flight section, and the DE arc is the elliptical orbit. The program turning section connects the vertical ascending section and the gravity oblique flight section. Under the control of an external torque, the rocket body rotates at a certain angle. The external torque is canceled after the completion of this section, and then the rocket enters the oblique flight state. The first stage rocket is usually subject to the propulsion of “vertical section + program turning section (external torque) + front range of gravity oblique flight section” (depending on the characteristics of the rocket engine). The rear range of gravity oblique flight section is completed by the second and third stage rockets one after another. Because the gravity of the Earth and the thrust of the rocket are not in the same straight line in the oblique flight state, the trajectory of the rocket body center is a smooth curve with a certain radian [17].Figure 1
Schematic diagram of the orbit of a spacecraft in the ascent phase.To describe the movement of observation satellites and the spacecraft, an appropriate coordinate system needs to be established [18]. The basic coordinate system adopted in this study is a coordinate system that can shift with the geocenter. It takes the geocenter Oc as the origin, the rotation axis of the Earth as the z-axis, and the North Pole as the positive direction. The x-axis points from Oc to the 0 meridian at zero time. Based on this, the y-axis is determined according to the right-handed system. Finally, the rectangular coordinate system Oc-XcYcZc is established. The geocenter Oc moves in an elliptical orbit around the sun. Thus, the Oc-XcYcZc system is noninertial in theory. Since the Earth’s revolution period is much longer than the observation arc duration of the spacecraft, this system is regarded as an inertial coordinate system in a short time; that is, this basic coordinate system does not rotate with the Earth.The second coordinate system is the observation coordinate systemOs-XsYsZs that moves with the satellite (Figure 2). The satellite center Os is taken as the origin. The Xs axis is along the line of OcOs, and the direction away from the Earth is positive. The Zs axis and Xs are perpendicular to the due north. The Ys axis is determined by the right-handed system. Since the orbit of general survey satellites is not strictly through the north and south poles, the definition of this coordinate system is clear. The observation coordinate system defined in this way is called the UEN coordinate system because the three coordinate axes, respectively, point to the directions of UP, EAST, and NORTH.Figure 2
Schematic diagram of the observation coordinate system.According to the dynamics of variable mass point, the simplified motion in the ascent phase of the spacecraft in the basic coordinate system can be expressed as follows [19]:(1)r⟶¨t=F⇀e+F⇀T=−Gmr⟶t3r⇀t+v⇀rtm˙tmt.In the above equation,F⇀e represents the sum of accelerations of external forces on the spacecraft; F⇀T represents the thrust acceleration generated by the rocket; m(t) is the instantaneous mass; m˙t is the rate of mass change; r⇀t is the position vector of the spacecraft in the basic coordinate system; r⇀¨t represents the second derivative of r⇀t concerning time t, that is, acceleration; Gm is the Earth’s gravitational constant (Gm=3.986005×1014m3/s2 is taken in this paper). To more clearly express the direction of the thrust acceleration, v⟶rt is the inverse vector of the fuel injection velocity relative to the rocket tail nozzle. (1) can represent the simplified motion of the observation satellite if only the first term on the right side of this equation is retained [20]:(2)r⇀¨t=F⇀e=−Gmr⟶t3r⇀t.Given the initial value of the position and the velocity in the basic coordinate system, the numerical solution method of ordinary differential equations can be exploited to calculate the motion trajectory of the spacecraft. Different spacecraft differ fundamentally in the models ofv⇀rt and mt. Generally, mt should be strictly monotonically decreasing nonnegative functions. The direction of v⇀rt is close to or the same as the velocity direction of the spacecraft, and v⇀rt is a constant vector [21].


## 3. Spacecraft Positioning and Orbit Estimation Based on the Point-by-Point Intersection of Two Satellites

### 3.1. Two-Satellite Point-by-Point Intersection Positioning Algorithm for the Spacecraft

Based on the idea of point-by-point intersection positioning of solid geometry, a single spacecraft can be located through two satellites. Specifically, two satellites are connected to the target. Theoretically, the two straight lines intersect at one point, that is, the target point. However, since the two straight lines are calculated separately from the observation data of the two satellites, they do not intersect. As for positioning, the midpoint of the common perpendicular segment of the two straight lines is taken as the positioning result in this paper.Through simplification, the observation data of the spacecraft by the observation satellite can be determined by two dimensionless ratios in the observation coordinate system:(3)α=ysxs;β=zsxs,where xs, ys, and zs are the coordinates of the spacecraft in the observation coordinate system. By observing the target from the first satellite, current α1,β1 can be obtained, that is, the direction of the target relative to the satellite observation coordinate system. The coordinate of the first satellite in the basic coordinate system is Os1x1,y1,z1, while that of the target point in the basic coordinate system is oxs,ys,zs. The parameter equation of the line Os1−o that connects the center of the satellite observation coordinate system and the target is set to(4)x,y,zT=K1x,K1y,K1zT⋅t.Similarly, for the second satellite, it is set to(5)x,y,zT=K2x,K2y,K2zT⋅t,where t denotes the present moment.Hence,Kix,Kiy,Kiz can be represented as 1/βi,αi/βi,1, where αi=yis/xis,βi=zis/xis. xis,yis,zis is the coordinate of the target point oxs,ys,zs in the i-th observation coordinate system.The point-by-point intersection calculation is based on the transformation of the coordinatexi,yi,zi of the target point o in the two observation coordinate systems into the coordinate xiA,yiA,ziA in the basic coordinate system. Firstly, the transformation relationship between the UEN coordinate system and the basic coordinate system [22, 23] is determined:(6)xiA,yiA,ziAT=Ai⋅Kix,Kiy,KizT+x0,y0,z0,i=1,2.In (6), Ai is the rotation matrix of space coordinate system.(7)Ai=sinθi1cosθi10cosθi1sinθi2−sinθi2−sinθi1cosθi2cosθi1sinθi2cosθi2−sinθi1sinθi2,where θi1 represents the elevation angle of the i-th observation coordinate system relative to the basic coordinate system and sinθi2 represents the counterclockwise rotation angle of the i-th observation coordinate system relative to the horizontal plane of the basic coordinate system. Therefore, (4) and (5) in the observation coordinate system can be converted into the two following parameter equations in the basic coordinate system:(8)x,y,zT=x1A−x1,y1A−y1,z1A−z1T⋅t1+x1,y1,z1T,(9)x,y,zT=x2A−x2,y2A−y2,z2A−z2T⋅t2+x2,y2,z2T.In theory, the two straight lines should intersect, and the intersection point is the positioning point at which the two satellites intersect. However, it is known from the observation data that, in general, the straight linesL1 and L2 are noncoplanar straight lines. Based on this, the coordinates of the endpoints of the common perpendicular segment can be obtained. Meanwhile, through the midpoint coordinate equation, the midpoint coordinates of the common vertical line segment of the noncoplanar straight lines L1 and L2 can be obtained. Then, this midpoint is used as the positioning point after the two satellites intersect. The schematic diagram of the point-by-point intersection positioning of the two satellites is illustrated in Figure 3.Figure 3
Schematic diagram of point-by-point intersection positioning of two satellites.


### 3.2. Establishment of the Orbit Estimation Model of the Spacecraft

Section3.1 describes the positioning method of the spacecraft by the two satellites. According to the kinematic characteristics of the spacecraft in the ascent phase, the model of the spacecraft in the moving period can also be established.To estimate the orbit of spacecraft 0 under (1), it is necessary to obtain the change curve of FT over time. Based on the observation of the characteristics of the curve, two appropriate v⇀rt and mt models can be established. Within a period of time, the three-dimensional position coordinates of the spacecraft, r⇀mt, at each sampling point can be obtained by the two-satellite positioning algorithm introduced in Section 3.1. Based on this, the velocity r⇀˙t and acceleration r⇀¨t of the spacecraft during this period of time can be obtained through the central difference of r⇀t. Besides, due to the excessive noise in the process of processing acceleration, the impact of noise can be reduced by adding a nonphase delay filtering.By shifting the term in (1), the following equation can be obtained:(10)r⇀¨t+Gmr⇀t3=v⇀rtm˙tmt=F⇀T.Sincev⇀rtm˙t/mt is unknown, the exact orbital motion of the spacecraft is still not available from the model. Therefore, heuristic methods are first adopted in this paper to estimate the model structure. Then, the model parameters are estimated through real data, and the precision of the conjectural model is finally verified through data. In this way, by applying the existing real data, a line chart corresponding to each component of F⇀T and time is shown in Figure 4.Figure 4
Three-component estimation of the rocket thrust acceleration.Based on the observation of the variation trend ofF⇀T, it is inferred that the spacecraft(11)F⇀T=v⇀rtm˙tmt=r⇀¨t+Gmr⇀t3r⇀t=Kt+B.That is,(12)v⟶rtm˙tmt=c1kt+c1b,c2kt+c2b,c3kt+c3bT=Kt+B,where K=c1k,c2k,c3kT and B=c1b,c2b,c3bT are the parameters to be estimated. For the i-th data point, (12) is simplified to the following linear equation:(13)F⇀TiT=FTxi,FTyi,FTziT=Kx,Ky,KzT⋅ti+Bx,By,BzT.If there are sufficient data points, (13) is transformed into the following least square problem:(14)t11t2⋮tN1⋮1⋅KTBT=F⟶T1F⟶T2⋮F⟶TN,where N is the number of datasets. Let A=t1,t2,...,tN1,1,...,1T, and then the optimal solution of the least square problem can be obtained:(15)KTBTT=AT⋅A−1⋅AT⋅F⇀T.Denoting the optimal solution asK∗TB∗TT, they can be applied to (1) to obtain(16)r⟶¨t=−Gmr⟶t3r⟶t+K∗t+B∗.Equation (16) is a second-order ordinary differential equation. The Runge-Kutta algorithm [24, 25] can be applied to obtain the numerical solution of this equation, that is, the movement trajectory of the spacecraft.As shown in Figure4, although the three components of F⇀T have similar changes in a straight line, they still have slight bends. Besides, the linear attenuation of the spacecraft mass is realistic. When m(t) is attenuated linearly, the component of F⇀T is a hyperbola. Therefore, a linear attenuation of m(t) is considered in this study, and v⇀rt is a constant vector that should have high precision. Based on this, another conjecture is made; that is, there is a hyperbolic relationship between F⇀T and v⟶rtm˙t/mt. This relationship is described in the following equation:(17)F⟶Tt=v⟶rtm˙tmt=v⟶rt.kkt+b.LetM=b/k and N⇀=vrx,vry,vrz, and (17) can be simplified into(18)F⇀T⋅t=−F⇀T⋅M+N⇀.Similarly, following the least square method, the optimal parameter of model 2 can be obtained, that is,M∗ and N⇀∗. Since the velocity residuals and position residuals of the hyperbolic model are lower than those of the linear regression model (the specific reasons for this will be given in Section 6), the hyperbolic model will be taken as a spacecraft orbit model in the following research.


## 3.1. Two-Satellite Point-by-Point Intersection Positioning Algorithm for the Spacecraft

Based on the idea of point-by-point intersection positioning of solid geometry, a single spacecraft can be located through two satellites. Specifically, two satellites are connected to the target. Theoretically, the two straight lines intersect at one point, that is, the target point. However, since the two straight lines are calculated separately from the observation data of the two satellites, they do not intersect. As for positioning, the midpoint of the common perpendicular segment of the two straight lines is taken as the positioning result in this paper.Through simplification, the observation data of the spacecraft by the observation satellite can be determined by two dimensionless ratios in the observation coordinate system:(3)α=ysxs;β=zsxs,where xs, ys, and zs are the coordinates of the spacecraft in the observation coordinate system. By observing the target from the first satellite, current α1,β1 can be obtained, that is, the direction of the target relative to the satellite observation coordinate system. The coordinate of the first satellite in the basic coordinate system is Os1x1,y1,z1, while that of the target point in the basic coordinate system is oxs,ys,zs. The parameter equation of the line Os1−o that connects the center of the satellite observation coordinate system and the target is set to(4)x,y,zT=K1x,K1y,K1zT⋅t.Similarly, for the second satellite, it is set to(5)x,y,zT=K2x,K2y,K2zT⋅t,where t denotes the present moment.Hence,Kix,Kiy,Kiz can be represented as 1/βi,αi/βi,1, where αi=yis/xis,βi=zis/xis. xis,yis,zis is the coordinate of the target point oxs,ys,zs in the i-th observation coordinate system.The point-by-point intersection calculation is based on the transformation of the coordinatexi,yi,zi of the target point o in the two observation coordinate systems into the coordinate xiA,yiA,ziA in the basic coordinate system. Firstly, the transformation relationship between the UEN coordinate system and the basic coordinate system [22, 23] is determined:(6)xiA,yiA,ziAT=Ai⋅Kix,Kiy,KizT+x0,y0,z0,i=1,2.In (6), Ai is the rotation matrix of space coordinate system.(7)Ai=sinθi1cosθi10cosθi1sinθi2−sinθi2−sinθi1cosθi2cosθi1sinθi2cosθi2−sinθi1sinθi2,where θi1 represents the elevation angle of the i-th observation coordinate system relative to the basic coordinate system and sinθi2 represents the counterclockwise rotation angle of the i-th observation coordinate system relative to the horizontal plane of the basic coordinate system. Therefore, (4) and (5) in the observation coordinate system can be converted into the two following parameter equations in the basic coordinate system:(8)x,y,zT=x1A−x1,y1A−y1,z1A−z1T⋅t1+x1,y1,z1T,(9)x,y,zT=x2A−x2,y2A−y2,z2A−z2T⋅t2+x2,y2,z2T.In theory, the two straight lines should intersect, and the intersection point is the positioning point at which the two satellites intersect. However, it is known from the observation data that, in general, the straight linesL1 and L2 are noncoplanar straight lines. Based on this, the coordinates of the endpoints of the common perpendicular segment can be obtained. Meanwhile, through the midpoint coordinate equation, the midpoint coordinates of the common vertical line segment of the noncoplanar straight lines L1 and L2 can be obtained. Then, this midpoint is used as the positioning point after the two satellites intersect. The schematic diagram of the point-by-point intersection positioning of the two satellites is illustrated in Figure 3.Figure 3
Schematic diagram of point-by-point intersection positioning of two satellites.


## 3.2. Establishment of the Orbit Estimation Model of the Spacecraft

Section3.1 describes the positioning method of the spacecraft by the two satellites. According to the kinematic characteristics of the spacecraft in the ascent phase, the model of the spacecraft in the moving period can also be established.To estimate the orbit of spacecraft 0 under (1), it is necessary to obtain the change curve of FT over time. Based on the observation of the characteristics of the curve, two appropriate v⇀rt and mt models can be established. Within a period of time, the three-dimensional position coordinates of the spacecraft, r⇀mt, at each sampling point can be obtained by the two-satellite positioning algorithm introduced in Section 3.1. Based on this, the velocity r⇀˙t and acceleration r⇀¨t of the spacecraft during this period of time can be obtained through the central difference of r⇀t. Besides, due to the excessive noise in the process of processing acceleration, the impact of noise can be reduced by adding a nonphase delay filtering.By shifting the term in (1), the following equation can be obtained:(10)r⇀¨t+Gmr⇀t3=v⇀rtm˙tmt=F⇀T.Sincev⇀rtm˙t/mt is unknown, the exact orbital motion of the spacecraft is still not available from the model. Therefore, heuristic methods are first adopted in this paper to estimate the model structure. Then, the model parameters are estimated through real data, and the precision of the conjectural model is finally verified through data. In this way, by applying the existing real data, a line chart corresponding to each component of F⇀T and time is shown in Figure 4.Figure 4
Three-component estimation of the rocket thrust acceleration.Based on the observation of the variation trend ofF⇀T, it is inferred that the spacecraft(11)F⇀T=v⇀rtm˙tmt=r⇀¨t+Gmr⇀t3r⇀t=Kt+B.That is,(12)v⟶rtm˙tmt=c1kt+c1b,c2kt+c2b,c3kt+c3bT=Kt+B,where K=c1k,c2k,c3kT and B=c1b,c2b,c3bT are the parameters to be estimated. For the i-th data point, (12) is simplified to the following linear equation:(13)F⇀TiT=FTxi,FTyi,FTziT=Kx,Ky,KzT⋅ti+Bx,By,BzT.If there are sufficient data points, (13) is transformed into the following least square problem:(14)t11t2⋮tN1⋮1⋅KTBT=F⟶T1F⟶T2⋮F⟶TN,where N is the number of datasets. Let A=t1,t2,...,tN1,1,...,1T, and then the optimal solution of the least square problem can be obtained:(15)KTBTT=AT⋅A−1⋅AT⋅F⇀T.Denoting the optimal solution asK∗TB∗TT, they can be applied to (1) to obtain(16)r⟶¨t=−Gmr⟶t3r⟶t+K∗t+B∗.Equation (16) is a second-order ordinary differential equation. The Runge-Kutta algorithm [24, 25] can be applied to obtain the numerical solution of this equation, that is, the movement trajectory of the spacecraft.As shown in Figure4, although the three components of F⇀T have similar changes in a straight line, they still have slight bends. Besides, the linear attenuation of the spacecraft mass is realistic. When m(t) is attenuated linearly, the component of F⇀T is a hyperbola. Therefore, a linear attenuation of m(t) is considered in this study, and v⇀rt is a constant vector that should have high precision. Based on this, another conjecture is made; that is, there is a hyperbolic relationship between F⇀T and v⟶rtm˙t/mt. This relationship is described in the following equation:(17)F⟶Tt=v⟶rtm˙tmt=v⟶rt.kkt+b.LetM=b/k and N⇀=vrx,vry,vrz, and (17) can be simplified into(18)F⇀T⋅t=−F⇀T⋅M+N⇀.Similarly, following the least square method, the optimal parameter of model 2 can be obtained, that is,M∗ and N⇀∗. Since the velocity residuals and position residuals of the hyperbolic model are lower than those of the linear regression model (the specific reasons for this will be given in Section 6), the hyperbolic model will be taken as a spacecraft orbit model in the following research.


## 4. Analysis and Solution of the Spacecraft Orbit Estimation Model considering System Errors

### 4.1. The Introduction of System Errors of the Spacecraft

Observation data inevitably contains various errors, such as random errors and system errors. In this paper, it is assumed that the random errors are white noise directly superimposed on the observation data, and the noise may be caused by background radiation interference and information processing. Meanwhile, system errors are caused by a few factors, such as satellite positioning errors, pointing mechanism errors, image calibration errors, and sensor installation errors. In this paper, only the system errors related to the satellite platform are considered. In particular, the system errors of different observation satellites are not related to each other, and the system errors of the same observation satellite for different spacecraft are the same.By applying appropriate simplified models, the system errors can be converted into the origin position errors of the observation coordinate system and three-axis pointing errors. According to engineering experiences, the origin position errors have little impact, while the three-axis pointing errors matter. The estimation of the three-axis pointing errors helps to improve the precision of the orbit estimation. Therefore, this paper only studies the three-axis pointing errors. Specifically, the error is expressed as two translation errors and one rotation error on the two-dimensional observation data plane. Based on this, it can be represented by three small constantsdα,dβ,dθ. α, β, and θ, respectively, represent the translation amount of the first observation, the translation amount of the second observation, and the rotation amount of the observation in the αβ plane.


### 4.2. Establishment and Solution of the Spacecraft Orbit Estimation Model considering System Errors

Based on the three-axis pointing errorsdα,dβ,dθ, the original observation coordinate system is corrected. According to the rotation matrix, the coordinate conversion relationship can be obtained as follows:(19)α′,β′,1T=cosθsinθdαsinθcosθdβ001⋅α,β,1T.α′,β′T is corrected α,βT. Applying α′,β′T to (6), subsequent observations can be obtained. At this time, under the basic coordinate system, it can be assumed that the numbers of the two observation satellites are 01 and 02, respectively. Two observation satellites have formed stereoscopic overlapping observation of spacecraft 0. Then the observation equation of the two satellites to the target point is modified as(20)x−x6,y−y6,z−z6T=1βi′,ai′βi′,1T⋅t,x−x9,y−y9,z−z9T=1βi′,αi′βi′,1T⋅t.In an ideal state, the lines connecting the two observation satellites to the target spacecraft should intersect at the target point. Thus, by making (8) and (9) approximately equal, the following equation holds:(21)x1A−x1y1A−y1z1A−z1⋅t1+x1y1z1=x2A−x2y2A−y2z2A−z2⋅t2+x2y2z2+e⇀t,where t1 and t2 are two parameters of the linear parameter equation before correction. It is suggested that the result of error correction has little impact on these parameters, so they can be adopted as known variables. e⇀t represents the errors in the three directions. Because the operation of the rotation matrix introduces nonlinearity, the estimation of the system error becomes a nonlinear least square problem.At each sampling moment, the observation equations of the observation satellites 01 and 02 to the target point in the basic coordinate system contain three unknown parameters. The nonlinear least square method is directly applied to estimate parametersdα1,dβ1,θ1,dα2,dβ2,θ2 in the static nonlinear model to minimize the sum of squared errors of the elements in e⇀t. Meanwhile, the shortest length of the common perpendicular segment is obtained, which satisfies the constraint condition of (21).After the calculated parametersdα1,dβ1,θ1,dα2,dβ2,θ2 are, respectively, applied to the observation equation of the two observation satellites to the target, the corrected observation equation is obtained. Then, following the point-by-point intersection method introduced in Section 3.1, the new observation point of the spacecraft can be calculated.


## 4.1. The Introduction of System Errors of the Spacecraft

Observation data inevitably contains various errors, such as random errors and system errors. In this paper, it is assumed that the random errors are white noise directly superimposed on the observation data, and the noise may be caused by background radiation interference and information processing. Meanwhile, system errors are caused by a few factors, such as satellite positioning errors, pointing mechanism errors, image calibration errors, and sensor installation errors. In this paper, only the system errors related to the satellite platform are considered. In particular, the system errors of different observation satellites are not related to each other, and the system errors of the same observation satellite for different spacecraft are the same.By applying appropriate simplified models, the system errors can be converted into the origin position errors of the observation coordinate system and three-axis pointing errors. According to engineering experiences, the origin position errors have little impact, while the three-axis pointing errors matter. The estimation of the three-axis pointing errors helps to improve the precision of the orbit estimation. Therefore, this paper only studies the three-axis pointing errors. Specifically, the error is expressed as two translation errors and one rotation error on the two-dimensional observation data plane. Based on this, it can be represented by three small constantsdα,dβ,dθ. α, β, and θ, respectively, represent the translation amount of the first observation, the translation amount of the second observation, and the rotation amount of the observation in the αβ plane.


## 4.2. Establishment and Solution of the Spacecraft Orbit Estimation Model considering System Errors

Based on the three-axis pointing errorsdα,dβ,dθ, the original observation coordinate system is corrected. According to the rotation matrix, the coordinate conversion relationship can be obtained as follows:(19)α′,β′,1T=cosθsinθdαsinθcosθdβ001⋅α,β,1T.α′,β′T is corrected α,βT. Applying α′,β′T to (6), subsequent observations can be obtained. At this time, under the basic coordinate system, it can be assumed that the numbers of the two observation satellites are 01 and 02, respectively. Two observation satellites have formed stereoscopic overlapping observation of spacecraft 0. Then the observation equation of the two satellites to the target point is modified as(20)x−x6,y−y6,z−z6T=1βi′,ai′βi′,1T⋅t,x−x9,y−y9,z−z9T=1βi′,αi′βi′,1T⋅t.In an ideal state, the lines connecting the two observation satellites to the target spacecraft should intersect at the target point. Thus, by making (8) and (9) approximately equal, the following equation holds:(21)x1A−x1y1A−y1z1A−z1⋅t1+x1y1z1=x2A−x2y2A−y2z2A−z2⋅t2+x2y2z2+e⇀t,where t1 and t2 are two parameters of the linear parameter equation before correction. It is suggested that the result of error correction has little impact on these parameters, so they can be adopted as known variables. e⇀t represents the errors in the three directions. Because the operation of the rotation matrix introduces nonlinearity, the estimation of the system error becomes a nonlinear least square problem.At each sampling moment, the observation equations of the observation satellites 01 and 02 to the target point in the basic coordinate system contain three unknown parameters. The nonlinear least square method is directly applied to estimate parametersdα1,dβ1,θ1,dα2,dβ2,θ2 in the static nonlinear model to minimize the sum of squared errors of the elements in e⇀t. Meanwhile, the shortest length of the common perpendicular segment is obtained, which satisfies the constraint condition of (21).After the calculated parametersdα1,dβ1,θ1,dα2,dβ2,θ2 are, respectively, applied to the observation equation of the two observation satellites to the target, the corrected observation equation is obtained. Then, following the point-by-point intersection method introduced in Section 3.1, the new observation point of the spacecraft can be calculated.


## 5. Discussion on Single-Satellite Single-Spacecraft Positioning and Multisatellite, Multispacecraft Cooperative Positioning

### 5.1. Analysis and Solution of the Single-Satellite Positioning Problem

In certain cases, this study may only meet the conditions for single-satellite and single-spacecraft positioning, which makes the consideration of the single-satellite positioning method necessary. The single satellite under study is recorded as No. 02, and the observed spacecraft is recorded as No. 01. Then, the orbit estimation of spacecraft No. 01 by satellite No. 02 is considered in this section.Assumption 1.
Due to the influence of atmosphere forces, it is difficult to observe the tail flame of the rocket in the vertical ascent phase. Moreover, the motion equation in the turning stage of the program is complicated. Therefore, this paper focuses on the backstage of the gravity oblique flight stage, and the simulation data in this paper also focus on this stage.
According to the analysis in Section4.1, it can be considered that the system errors of the observation satellite are known. Thus, the observation linear equation l3 of the satellite to the target point is known. At the sampling time ti, the observation point xi9,yi9,zi9 is obtained from the observation satellite, and it conforms to the linear equation l3:(22)xi9,yi9,zi9T=x1Ai−x1i,y1Ai−y1i,z1Ai−z1iT⋅ti+x1i,y1i,z1iT.
Supposing that the orbit estimation model of the single spacecraft is the hyperbolic regression model established in Section3.2, (17) can be applied to obtain the observation point (xi9, yi9, zi9):(23)Gm⋅xi9xi92+yi92+zi923/2,Gm⋅yi9xi92+yi92+zi923/2,Gm⋅zi9xi92+yi92+zi923/2T+x¨i9,y¨i9,z¨i9T=N⇀1t+M,where N⇀=Nx,Ny,Nz and M are unknown parameters to be estimated and ti is the i-th sampling time. Combining (23) and (24), given a certain amount of sampled data, the position of spacecraft No. 0 at each sampling time can be calculated. The details of the calculation will be introduced in Section 6. The schematic diagram of the single-satellite spacecraft positioning is shown in Figure 5.
It is worth mentioning that the above problems can be solved by some numerical optimization algorithms, such as filtering algorithm [26–28], genetic algorithm [29, 30], and particle swarm algorithm [31]. Through the above methods, the system errors of spacecraft No. 02, the trajectory equation of spacecraft No. 01, and the three-dimensional position coordinates of spacecraft No. 01 can be obtained through the observation of satellite No. 02 in the basic coordinate system at any sampling time in the interval [50 s, 170 s] with a step length of 0.2 s.Figure 5
Schematic diagram of the single-satellite positioning.


### 5.2. Analysis of the Multisatellite, Multispacecraft Positioning Problem

In a few cases, there are many spacecraft and observation satellites. This section presents a multisatellite and multispacecraft collaborative positioning scheme.It is assumed that there arem observation satellites and n spacecraft. Firstly, each spacecraft is located. For spacecraft nj, the observations of every two of the m observation satellites can be exploited to locate the spacecraft separately. The obtained Cm2 observation coordinates are averaged, and the average value is defined as the coordinate of the multisatellite collaborative positioning. The estimated coordinate of the j-th spacecraft at time k is(24)xjk,yjk,zjkT=∑i=1Cm2xij,yij,zijCm2.Letdαi,dβi,θi be the system errors of the i-th satellite, and let αijk,βijk be the observation of the i-th satellite to the j-th spacecraft at time k.For the same spacecraftj, it is observed by a satellite with the other (m-1) satellites. Set dαlij,dβlij,θlijT as the i-th observation satellite to collaborate with the l-th observation satellite. As for the estimated system errors of the j-th spacecraft observation, each satellite can obtain (m-1) groups of dαlij,dβij,dθlij,T, i = 1, 2, …, m-1. Then, for every other spacecraft, there are (m-1) groups of dαlij,dβlij,dθlijT, j = 1, 2, …, n. Therefore, by averaging dαlij,dβlij,dθlij,T of all the groups, the system errors of the i-th observation satellite can be obtained:(25)dαi,dβi,dθi=∑l=1m−1∑j=1ndαlij,dβlij,dθlij,Tn⋅m−1.According to (24) and (25), the orbit of each spacecraft and the system errors of each observation satellite can be estimated.It should be noted that this model is based on the idea of point-by-point intersection and positioning of two satellites. Therefore, the conditions should be retained. Also, since there is no sufficient data at present, the model only provides an algorithm in theory.


## 5.1. Analysis and Solution of the Single-Satellite Positioning Problem

In certain cases, this study may only meet the conditions for single-satellite and single-spacecraft positioning, which makes the consideration of the single-satellite positioning method necessary. The single satellite under study is recorded as No. 02, and the observed spacecraft is recorded as No. 01. Then, the orbit estimation of spacecraft No. 01 by satellite No. 02 is considered in this section.Assumption 1.
Due to the influence of atmosphere forces, it is difficult to observe the tail flame of the rocket in the vertical ascent phase. Moreover, the motion equation in the turning stage of the program is complicated. Therefore, this paper focuses on the backstage of the gravity oblique flight stage, and the simulation data in this paper also focus on this stage.
According to the analysis in Section4.1, it can be considered that the system errors of the observation satellite are known. Thus, the observation linear equation l3 of the satellite to the target point is known. At the sampling time ti, the observation point xi9,yi9,zi9 is obtained from the observation satellite, and it conforms to the linear equation l3:(22)xi9,yi9,zi9T=x1Ai−x1i,y1Ai−y1i,z1Ai−z1iT⋅ti+x1i,y1i,z1iT.
Supposing that the orbit estimation model of the single spacecraft is the hyperbolic regression model established in Section3.2, (17) can be applied to obtain the observation point (xi9, yi9, zi9):(23)Gm⋅xi9xi92+yi92+zi923/2,Gm⋅yi9xi92+yi92+zi923/2,Gm⋅zi9xi92+yi92+zi923/2T+x¨i9,y¨i9,z¨i9T=N⇀1t+M,where N⇀=Nx,Ny,Nz and M are unknown parameters to be estimated and ti is the i-th sampling time. Combining (23) and (24), given a certain amount of sampled data, the position of spacecraft No. 0 at each sampling time can be calculated. The details of the calculation will be introduced in Section 6. The schematic diagram of the single-satellite spacecraft positioning is shown in Figure 5.
It is worth mentioning that the above problems can be solved by some numerical optimization algorithms, such as filtering algorithm [26–28], genetic algorithm [29, 30], and particle swarm algorithm [31]. Through the above methods, the system errors of spacecraft No. 02, the trajectory equation of spacecraft No. 01, and the three-dimensional position coordinates of spacecraft No. 01 can be obtained through the observation of satellite No. 02 in the basic coordinate system at any sampling time in the interval [50 s, 170 s] with a step length of 0.2 s.Figure 5
Schematic diagram of the single-satellite positioning.


## 5.2. Analysis of the Multisatellite, Multispacecraft Positioning Problem

In a few cases, there are many spacecraft and observation satellites. This section presents a multisatellite and multispacecraft collaborative positioning scheme.It is assumed that there arem observation satellites and n spacecraft. Firstly, each spacecraft is located. For spacecraft nj, the observations of every two of the m observation satellites can be exploited to locate the spacecraft separately. The obtained Cm2 observation coordinates are averaged, and the average value is defined as the coordinate of the multisatellite collaborative positioning. The estimated coordinate of the j-th spacecraft at time k is(24)xjk,yjk,zjkT=∑i=1Cm2xij,yij,zijCm2.Letdαi,dβi,θi be the system errors of the i-th satellite, and let αijk,βijk be the observation of the i-th satellite to the j-th spacecraft at time k.For the same spacecraftj, it is observed by a satellite with the other (m-1) satellites. Set dαlij,dβlij,θlijT as the i-th observation satellite to collaborate with the l-th observation satellite. As for the estimated system errors of the j-th spacecraft observation, each satellite can obtain (m-1) groups of dαlij,dβij,dθlij,T, i = 1, 2, …, m-1. Then, for every other spacecraft, there are (m-1) groups of dαlij,dβlij,dθlijT, j = 1, 2, …, n. Therefore, by averaging dαlij,dβlij,dθlij,T of all the groups, the system errors of the i-th observation satellite can be obtained:(25)dαi,dβi,dθi=∑l=1m−1∑j=1ndαlij,dβlij,dθlij,Tn⋅m−1.According to (24) and (25), the orbit of each spacecraft and the system errors of each observation satellite can be estimated.It should be noted that this model is based on the idea of point-by-point intersection and positioning of two satellites. Therefore, the conditions should be retained. Also, since there is no sufficient data at present, the model only provides an algorithm in theory.


## 6. Simulation Analysis

To verify the effectiveness of the proposed algorithm, this paper uses real data to construct a simulation platform on MATLAB for verification. In the simulation, some data from spacecraft No. 0 and the observation satellites No. 01 and No. 02 are adopted.


### 6.1. Satellite Motion Trajectory and Satellite Observation Trajectory Simulation

The position of the observation satellite at any time is the prerequisite for estimating the spacecraft information. However, this position is not known and needs to be calculated through existing data. By applying the known initial positionr⇀0 and initial acceleration r⇀¨0 of satellite No. 02 into (2), the motion trajectory r⇀t of satellite No. 02 can be obtained. Since the problem to be solved is a second-order nonlinear ordinary differential equation, it is difficult to obtain the analytical solution. To address this issue, the Runge-Kutta algorithm is adopted in this paper to obtain the numerical solution. The three-dimensional positions of observation satellite 02 at 50.0 s, 100.0 s, 150.0 s, 200.0 s, and 250.0 s are calculated using ODE45 function of MATLAB shown in Table 1.Table 1
The three-dimensional position of observation satellite 02 at the sampling point.
tsxmymzm50.01.773806e + 068.161384e + 064.516700e + 06100.01.501626e + 068.126764e + 064.684680e + 06150.01.227700e + 068.082699e + 064.847216e + 06200.09.523491e + 058.029253e + 065.004127e + 06250.06.956792e + 057.971293e + 065.144636e + 06In addition to the above points, the three-dimensional position of satellite can be observed in the basic coordinate system at any time. The simulation results are shown in Figure6. Among them, the running period is 2.5 hours.Figure 6
The trajectory of satellite No. 02 orbiting the Earth in 2.5 hours (m).
(a)(b)(c)(d)The simulation of the trajectory (meter) of satellite No. 02 orbiting the Earth in 2.5 hours is shown in Figure6. It can be seen from the trajectory that the observation satellite flies close to the ground surface and is operating well.


### 6.2. The Two-Satellite Point-by-Point Intersection Algorithm for the Prediction of the Spacecraft Trajectory

#### 6.2.1. Algorithm Implementation

To verify the effectiveness of the two-satellite point-by-point intersection positioning algorithm introduced in Section3.1, the observation satellites No. 01 and No. 02 form a three-dimensional overlap observation of spacecraft No. 0. In this case, the motion trajectory of spacecraft No. 0 is observed through the two-satellite point-by-point intersection positioning algorithm.Double-star point-by-point positioning requires data at the same time. Firstly, the simulation data of satellites No. 01 and No. 02 should be preprocessed for 3 times’ spline interpolation, so that the data are sampled at an interval of 0.2 s in the closed interval of [50 s, 170 s], so that the data time points are consistent. The algorithm flow is shown in Figure7.Figure 7
Flow chart of double-star location algorithm.In Figure8, the purple curve shows the observation position curve of spacecraft No. 0 in the interval [50 s, 170 s]. The red and blue curves are, respectively, the motion trajectories of observation satellites No. 01 and No. 02.Figure 8
Observation position of spacecraft No. 0 in the interval [50 s, 170 s].


#### 6.2.2. Orbit Estimation Model

Two kinds of spacecraft trajectory estimation models are introduced. Residual analysis is conducted on the two models to verify and compare the effectiveness of the models. In Section6.2.1, through the two-satellite point-by-point intersection algorithm, the motion trajectory of spacecraft No. 0 within [50 s, 170 s] is calculated. We select 10 s as a sampling interval and a total of 13 sampling points. After the corresponding acceleration F⇀Tt is obtained based on the second-order difference, the data calculated from 13 points are applied to the two spacecraft orbital models introduced in Section 3.2. In this way, the coefficients in Model 1 (linear regression model) are obtained as follows:(26)K=0.150694,0.0267504,0.185861T,B=−7.85566,15.9406,8.91922T.ParametersK and B are applied to (1), and the Runge-Kutta algorithm is exploited to find the motion trajectory of the spacecraft, according to the residual formula:(27)r=∑Si−Si′2N.In the above formula, N is 8.46286e + 02, Si is the spacecraft position obtained by the point-by-point intersection positioning algorithm, and Si′ is the result estimated by Model 1. The position residuals and the velocity residuals of spacecraft No. 0 in [50 s, 170 s] are listed in Table 2. Similarly, the same data is input into Model 2, and the parameters of Model 2 are solved and listed in Table 2.(28)M=−261.323753,N⇀=3753.699997,−1995.763805,−4290.903058T.Table 2
Position and velocity residuals of the two models.
Position residuals (m)Velocity residuals (m/s)Model 14.09361e + 048.46286e + 02Model 22.67343e + 046.39409e + 02It can be seen that Model 2 has smaller position and velocity residuals, and it achieves higher prediction precision. Thus, Model 2 is used as the final application model to calculate the estimated value of the displacement and corresponding acceleration of spacecraft No. 0.Model 1 (linear model) is an estimation model with compoundv⟶rtm˙t/mt, which contains the inverse vector v⟶rt and instantaneous mass mt of the jet velocity at the nozzle of the rocket tail; it has the characteristics of simple calculation. Model 2 (hyperbolic model) makes a specific estimate of v⟶rt and gives the corresponding approximate value. mt is generally a nonnegative strictly monotonically decreasing function. According to image feature of the thrust acceleration F⟶T, it shows a weak radian. In addition, the combustion rate m˙t can generally be set as a constant, so mt can be assumed to be a linear function. From Table 2, compared with Model 1, Model 2 has higher accuracy.


### 6.3. The Two-Satellite Point-by-Point Intersection Algorithm for the Prediction of the Spacecraft Orbit

It is mentioned in this section that the observation position through the particle filtering algorithm on the calculated results can effectively reduce the residuals. Therefore, the observation position obtained by the system model is filtered to obtain a more accurate observation position. Meanwhile, it is mentioned in this section that the parameter that needs to be adjusted for the particle filtering algorithm is mainly the number of particles. Based on the experience, the total number of particles is set to 1000.The final results are shown in the following table.According to the spacecraft observation model, the residuals before and after particle filtering can be calculated as in Table3.Table 3
Position and velocity comparison.
Position residuals (m)Velocity residuals (m/s)Without particle filtering2.67343e + 046.39409e + 02With particle filtering1.44085e + 046.39352e + 02The variation curves of the observed position and estimated trajectory versus errors are illustrated in Figure9. Compared with no particle filtering (Figure 9(a)), the position residuals after particle filtering (Figure 9(b)) are significantly reduced. Therefore, the effectiveness of particle filtering is verified. During the measurement process, the length of the common perpendicular line of the two trajectories is 100 meters, indicating that the position observation results are accurate. However, after two differences, a lot of noise appears in the obtained acceleration. Thus, the acceleration is smoothed. However, compared with the observed data, the motion model parameters obtained in theory suffer from larger errors. Compared with the situation without particle filtering (Figure 9(a)), the maximum distance error with particle filtering (Figure 9(b)) is reduced from 6×104 to 3.5×104. To sum up, it can be verified that the effect of particle filtering does improve the observation results.Figure 9
Variation curves of the observed position and estimated trajectory versus errors. (a) The estimation error without particle filtering. (b) The estimation error with particle filtering.
(a)(b)To further verify the effectiveness of the entire two-satellite positioning algorithm, three position components (t-x, t-y, and t-y) and three velocity components (t-vx, t-vy, and t-vz) of spacecraft No. 0 after particle filtering are shown in Tables 4–6 and Figure 10.Table 4
Observation of spacecraft 0.
Time (s)Observation position (m)Observation velocity (m/s)xyzvxvyvz50−1.1045785E + 066.2026543E + 061.1340530E + 06−7.4903883E + 027.8251039E + 029.5203116E + 0260−1.1131004E + 066.2103568E + 061.1442698E + 06−9.3028087E + 027.8195722E + 021.1023444E + 0370−1.1233125E + 066.2184561E + 061.1564104E + 06−1.0912921E + 038.3312097E + 021.3019919E + 0380−1.1352505E + 066.2270647E + 061.1706173E + 06−1.2918543E + 038.7182984E + 021.5207654E + 0390−1.1490173E + 066.2361108E + 061.1870506E + 06−1.4794826E + 039.3487786E + 021.7629851E + 03100−1.1646975E + 066.2455991E + 061.2058790E + 06−1.6627190E + 039.6992430E + 021.9973859E + 03110−1.1824947E + 066.2554542E + 061.2271690E + 06−1.8905454E + 031.0071563E + 032.2702085E + 03120−1.2025020E + 066.2656861E + 061.2511673E + 06−2.1194352E + 031.0520625E + 032.5351808E + 03130−1.2248703E + 066.2762167E + 061.2780546E + 06−2.3455157E + 031.0697871E + 032.8286792E + 03140−1.2497717E + 066.2871388E + 061.3080340E + 06−2.6280271E + 031.1075056E + 033.1616176E + 03150−1.2773908E + 066.2983536E + 061.3413011E + 06−2.8949518E + 031.1338597E + 033.4867523E + 03160−1.3079042E + 066.3099141E + 061.3781158E + 06−3.2136372E + 031.1845338E + 033.8839258E + 03170−1.3416105E + 066.3218029E + 061.4187374E + 06−3.5230289E + 031.1994513E + 034.2283819E + 03Table 5
Estimations of spacecraft 0.
Time (s)Estimation position (m)Estimation velocity (m/s)xyzvxvyvz50−1.1053335E + 066.2034360E + 061.1350124E + 06−7.6829756E + 027.7969235E + 039.7559569E + 0260−1.1135553E + 066.2111612E + 061.1454274E + 06−8.8420484E + 027.6487188E + 031.1173389E + 0370−1.1230364E + 066.2187452E + 061.1573817E + 06−1.0130099E + 037.5255725E + 031.2747035E + 0380−1.1338359E + 066.2262251E + 061.1709458E + 06−1.1479464E + 037.4408492E + 031.4393986E + 0390−1.1460175E + 066.2336408E + 061.1861961E + 06−1.2895491E + 037.3978234E + 031.6120649E + 03100−1.1596510E + 066.2410358E + 061.2032160E + 06−1.4384300E + 037.4002235E + 031.7934354E + 03110−1.1748126E + 066.2484577E + 061.2220966E + 06−1.5952926E + 037.4523179E + 031.9843530E + 03120−1.1915861E + 066.2559586E + 061.2429380E + 06−1.7609510E + 037.5590229E + 032.1857932E + 03130−1.2100640E + 066.2635958E + 061.2658512E + 06−1.9363541E + 037.7260386E + 032.3988928E + 03140−1.2303493E + 066.2714329E + 061.2909590E + 06−2.1226160E + 037.9600243E + 032.6249870E + 03150−1.2525568E + 066.2795407E + 061.3183993E + 06−2.3210565E + 038.2688256E + 032.8656573E + 03160−1.2768161E + 066.2879985E + 061.3483269E + 06−2.5332536E + 038.6617746E + 033.1227952E + 03170−1.3032739E + 066.2968959E + 061.3809175E + 06−2.7469620E + 039.1178683E + 033.3815571E + 03Table 6
Position and velocity comparison.
Position residuals (m)Velocity residuals (m/s)Spacecraft No. 02.67343e + 046.39409e + 02Figure 10
The three-component change curves of the velocity and displacement of spacecraft No. 0 in the case of Model 2.
(a)(b)It can be seen that the trajectory model of spacecraft No. 0 is consistent with the three-component change curve of velocity and displacement in the case of Model 2. Since the vertical axis in Figure10(a) has a large scale, the three-component position curve is approximately a straight line. Besides, in Figure 10(b), the three components of velocity change in different directions with time, resulting in changes in the direction and the velocity of the spacecraft. This conforms to the actual operation law of the spacecraft.


## 6.1. Satellite Motion Trajectory and Satellite Observation Trajectory Simulation

The position of the observation satellite at any time is the prerequisite for estimating the spacecraft information. However, this position is not known and needs to be calculated through existing data. By applying the known initial positionr⇀0 and initial acceleration r⇀¨0 of satellite No. 02 into (2), the motion trajectory r⇀t of satellite No. 02 can be obtained. Since the problem to be solved is a second-order nonlinear ordinary differential equation, it is difficult to obtain the analytical solution. To address this issue, the Runge-Kutta algorithm is adopted in this paper to obtain the numerical solution. The three-dimensional positions of observation satellite 02 at 50.0 s, 100.0 s, 150.0 s, 200.0 s, and 250.0 s are calculated using ODE45 function of MATLAB shown in Table 1.Table 1
The three-dimensional position of observation satellite 02 at the sampling point.
tsxmymzm50.01.773806e + 068.161384e + 064.516700e + 06100.01.501626e + 068.126764e + 064.684680e + 06150.01.227700e + 068.082699e + 064.847216e + 06200.09.523491e + 058.029253e + 065.004127e + 06250.06.956792e + 057.971293e + 065.144636e + 06In addition to the above points, the three-dimensional position of satellite can be observed in the basic coordinate system at any time. The simulation results are shown in Figure6. Among them, the running period is 2.5 hours.Figure 6
The trajectory of satellite No. 02 orbiting the Earth in 2.5 hours (m).
(a)(b)(c)(d)The simulation of the trajectory (meter) of satellite No. 02 orbiting the Earth in 2.5 hours is shown in Figure6. It can be seen from the trajectory that the observation satellite flies close to the ground surface and is operating well.


## 6.2. The Two-Satellite Point-by-Point Intersection Algorithm for the Prediction of the Spacecraft Trajectory

### 6.2.1. Algorithm Implementation

To verify the effectiveness of the two-satellite point-by-point intersection positioning algorithm introduced in Section3.1, the observation satellites No. 01 and No. 02 form a three-dimensional overlap observation of spacecraft No. 0. In this case, the motion trajectory of spacecraft No. 0 is observed through the two-satellite point-by-point intersection positioning algorithm.Double-star point-by-point positioning requires data at the same time. Firstly, the simulation data of satellites No. 01 and No. 02 should be preprocessed for 3 times’ spline interpolation, so that the data are sampled at an interval of 0.2 s in the closed interval of [50 s, 170 s], so that the data time points are consistent. The algorithm flow is shown in Figure7.Figure 7
Flow chart of double-star location algorithm.In Figure8, the purple curve shows the observation position curve of spacecraft No. 0 in the interval [50 s, 170 s]. The red and blue curves are, respectively, the motion trajectories of observation satellites No. 01 and No. 02.Figure 8
Observation position of spacecraft No. 0 in the interval [50 s, 170 s].


### 6.2.2. Orbit Estimation Model

Two kinds of spacecraft trajectory estimation models are introduced. Residual analysis is conducted on the two models to verify and compare the effectiveness of the models. In Section6.2.1, through the two-satellite point-by-point intersection algorithm, the motion trajectory of spacecraft No. 0 within [50 s, 170 s] is calculated. We select 10 s as a sampling interval and a total of 13 sampling points. After the corresponding acceleration F⇀Tt is obtained based on the second-order difference, the data calculated from 13 points are applied to the two spacecraft orbital models introduced in Section 3.2. In this way, the coefficients in Model 1 (linear regression model) are obtained as follows:(26)K=0.150694,0.0267504,0.185861T,B=−7.85566,15.9406,8.91922T.ParametersK and B are applied to (1), and the Runge-Kutta algorithm is exploited to find the motion trajectory of the spacecraft, according to the residual formula:(27)r=∑Si−Si′2N.In the above formula, N is 8.46286e + 02, Si is the spacecraft position obtained by the point-by-point intersection positioning algorithm, and Si′ is the result estimated by Model 1. The position residuals and the velocity residuals of spacecraft No. 0 in [50 s, 170 s] are listed in Table 2. Similarly, the same data is input into Model 2, and the parameters of Model 2 are solved and listed in Table 2.(28)M=−261.323753,N⇀=3753.699997,−1995.763805,−4290.903058T.Table 2
Position and velocity residuals of the two models.
Position residuals (m)Velocity residuals (m/s)Model 14.09361e + 048.46286e + 02Model 22.67343e + 046.39409e + 02It can be seen that Model 2 has smaller position and velocity residuals, and it achieves higher prediction precision. Thus, Model 2 is used as the final application model to calculate the estimated value of the displacement and corresponding acceleration of spacecraft No. 0.Model 1 (linear model) is an estimation model with compoundv⟶rtm˙t/mt, which contains the inverse vector v⟶rt and instantaneous mass mt of the jet velocity at the nozzle of the rocket tail; it has the characteristics of simple calculation. Model 2 (hyperbolic model) makes a specific estimate of v⟶rt and gives the corresponding approximate value. mt is generally a nonnegative strictly monotonically decreasing function. According to image feature of the thrust acceleration F⟶T, it shows a weak radian. In addition, the combustion rate m˙t can generally be set as a constant, so mt can be assumed to be a linear function. From Table 2, compared with Model 1, Model 2 has higher accuracy.


## 6.2.1. Algorithm Implementation

To verify the effectiveness of the two-satellite point-by-point intersection positioning algorithm introduced in Section3.1, the observation satellites No. 01 and No. 02 form a three-dimensional overlap observation of spacecraft No. 0. In this case, the motion trajectory of spacecraft No. 0 is observed through the two-satellite point-by-point intersection positioning algorithm.Double-star point-by-point positioning requires data at the same time. Firstly, the simulation data of satellites No. 01 and No. 02 should be preprocessed for 3 times’ spline interpolation, so that the data are sampled at an interval of 0.2 s in the closed interval of [50 s, 170 s], so that the data time points are consistent. The algorithm flow is shown in Figure7.Figure 7
Flow chart of double-star location algorithm.In Figure8, the purple curve shows the observation position curve of spacecraft No. 0 in the interval [50 s, 170 s]. The red and blue curves are, respectively, the motion trajectories of observation satellites No. 01 and No. 02.Figure 8
Observation position of spacecraft No. 0 in the interval [50 s, 170 s].


## 6.2.2. Orbit Estimation Model

Two kinds of spacecraft trajectory estimation models are introduced. Residual analysis is conducted on the two models to verify and compare the effectiveness of the models. In Section6.2.1, through the two-satellite point-by-point intersection algorithm, the motion trajectory of spacecraft No. 0 within [50 s, 170 s] is calculated. We select 10 s as a sampling interval and a total of 13 sampling points. After the corresponding acceleration F⇀Tt is obtained based on the second-order difference, the data calculated from 13 points are applied to the two spacecraft orbital models introduced in Section 3.2. In this way, the coefficients in Model 1 (linear regression model) are obtained as follows:(26)K=0.150694,0.0267504,0.185861T,B=−7.85566,15.9406,8.91922T.ParametersK and B are applied to (1), and the Runge-Kutta algorithm is exploited to find the motion trajectory of the spacecraft, according to the residual formula:(27)r=∑Si−Si′2N.In the above formula, N is 8.46286e + 02, Si is the spacecraft position obtained by the point-by-point intersection positioning algorithm, and Si′ is the result estimated by Model 1. The position residuals and the velocity residuals of spacecraft No. 0 in [50 s, 170 s] are listed in Table 2. Similarly, the same data is input into Model 2, and the parameters of Model 2 are solved and listed in Table 2.(28)M=−261.323753,N⇀=3753.699997,−1995.763805,−4290.903058T.Table 2
Position and velocity residuals of the two models.
Position residuals (m)Velocity residuals (m/s)Model 14.09361e + 048.46286e + 02Model 22.67343e + 046.39409e + 02It can be seen that Model 2 has smaller position and velocity residuals, and it achieves higher prediction precision. Thus, Model 2 is used as the final application model to calculate the estimated value of the displacement and corresponding acceleration of spacecraft No. 0.Model 1 (linear model) is an estimation model with compoundv⟶rtm˙t/mt, which contains the inverse vector v⟶rt and instantaneous mass mt of the jet velocity at the nozzle of the rocket tail; it has the characteristics of simple calculation. Model 2 (hyperbolic model) makes a specific estimate of v⟶rt and gives the corresponding approximate value. mt is generally a nonnegative strictly monotonically decreasing function. According to image feature of the thrust acceleration F⟶T, it shows a weak radian. In addition, the combustion rate m˙t can generally be set as a constant, so mt can be assumed to be a linear function. From Table 2, compared with Model 1, Model 2 has higher accuracy.


## 6.3. The Two-Satellite Point-by-Point Intersection Algorithm for the Prediction of the Spacecraft Orbit

It is mentioned in this section that the observation position through the particle filtering algorithm on the calculated results can effectively reduce the residuals. Therefore, the observation position obtained by the system model is filtered to obtain a more accurate observation position. Meanwhile, it is mentioned in this section that the parameter that needs to be adjusted for the particle filtering algorithm is mainly the number of particles. Based on the experience, the total number of particles is set to 1000.The final results are shown in the following table.According to the spacecraft observation model, the residuals before and after particle filtering can be calculated as in Table3.Table 3
Position and velocity comparison.
Position residuals (m)Velocity residuals (m/s)Without particle filtering2.67343e + 046.39409e + 02With particle filtering1.44085e + 046.39352e + 02The variation curves of the observed position and estimated trajectory versus errors are illustrated in Figure9. Compared with no particle filtering (Figure 9(a)), the position residuals after particle filtering (Figure 9(b)) are significantly reduced. Therefore, the effectiveness of particle filtering is verified. During the measurement process, the length of the common perpendicular line of the two trajectories is 100 meters, indicating that the position observation results are accurate. However, after two differences, a lot of noise appears in the obtained acceleration. Thus, the acceleration is smoothed. However, compared with the observed data, the motion model parameters obtained in theory suffer from larger errors. Compared with the situation without particle filtering (Figure 9(a)), the maximum distance error with particle filtering (Figure 9(b)) is reduced from 6×104 to 3.5×104. To sum up, it can be verified that the effect of particle filtering does improve the observation results.Figure 9
Variation curves of the observed position and estimated trajectory versus errors. (a) The estimation error without particle filtering. (b) The estimation error with particle filtering.
(a)(b)To further verify the effectiveness of the entire two-satellite positioning algorithm, three position components (t-x, t-y, and t-y) and three velocity components (t-vx, t-vy, and t-vz) of spacecraft No. 0 after particle filtering are shown in Tables 4–6 and Figure 10.Table 4
Observation of spacecraft 0.
Time (s)Observation position (m)Observation velocity (m/s)xyzvxvyvz50−1.1045785E + 066.2026543E + 061.1340530E + 06−7.4903883E + 027.8251039E + 029.5203116E + 0260−1.1131004E + 066.2103568E + 061.1442698E + 06−9.3028087E + 027.8195722E + 021.1023444E + 0370−1.1233125E + 066.2184561E + 061.1564104E + 06−1.0912921E + 038.3312097E + 021.3019919E + 0380−1.1352505E + 066.2270647E + 061.1706173E + 06−1.2918543E + 038.7182984E + 021.5207654E + 0390−1.1490173E + 066.2361108E + 061.1870506E + 06−1.4794826E + 039.3487786E + 021.7629851E + 03100−1.1646975E + 066.2455991E + 061.2058790E + 06−1.6627190E + 039.6992430E + 021.9973859E + 03110−1.1824947E + 066.2554542E + 061.2271690E + 06−1.8905454E + 031.0071563E + 032.2702085E + 03120−1.2025020E + 066.2656861E + 061.2511673E + 06−2.1194352E + 031.0520625E + 032.5351808E + 03130−1.2248703E + 066.2762167E + 061.2780546E + 06−2.3455157E + 031.0697871E + 032.8286792E + 03140−1.2497717E + 066.2871388E + 061.3080340E + 06−2.6280271E + 031.1075056E + 033.1616176E + 03150−1.2773908E + 066.2983536E + 061.3413011E + 06−2.8949518E + 031.1338597E + 033.4867523E + 03160−1.3079042E + 066.3099141E + 061.3781158E + 06−3.2136372E + 031.1845338E + 033.8839258E + 03170−1.3416105E + 066.3218029E + 061.4187374E + 06−3.5230289E + 031.1994513E + 034.2283819E + 03Table 5
Estimations of spacecraft 0.
Time (s)Estimation position (m)Estimation velocity (m/s)xyzvxvyvz50−1.1053335E + 066.2034360E + 061.1350124E + 06−7.6829756E + 027.7969235E + 039.7559569E + 0260−1.1135553E + 066.2111612E + 061.1454274E + 06−8.8420484E + 027.6487188E + 031.1173389E + 0370−1.1230364E + 066.2187452E + 061.1573817E + 06−1.0130099E + 037.5255725E + 031.2747035E + 0380−1.1338359E + 066.2262251E + 061.1709458E + 06−1.1479464E + 037.4408492E + 031.4393986E + 0390−1.1460175E + 066.2336408E + 061.1861961E + 06−1.2895491E + 037.3978234E + 031.6120649E + 03100−1.1596510E + 066.2410358E + 061.2032160E + 06−1.4384300E + 037.4002235E + 031.7934354E + 03110−1.1748126E + 066.2484577E + 061.2220966E + 06−1.5952926E + 037.4523179E + 031.9843530E + 03120−1.1915861E + 066.2559586E + 061.2429380E + 06−1.7609510E + 037.5590229E + 032.1857932E + 03130−1.2100640E + 066.2635958E + 061.2658512E + 06−1.9363541E + 037.7260386E + 032.3988928E + 03140−1.2303493E + 066.2714329E + 061.2909590E + 06−2.1226160E + 037.9600243E + 032.6249870E + 03150−1.2525568E + 066.2795407E + 061.3183993E + 06−2.3210565E + 038.2688256E + 032.8656573E + 03160−1.2768161E + 066.2879985E + 061.3483269E + 06−2.5332536E + 038.6617746E + 033.1227952E + 03170−1.3032739E + 066.2968959E + 061.3809175E + 06−2.7469620E + 039.1178683E + 033.3815571E + 03Table 6
Position and velocity comparison.
Position residuals (m)Velocity residuals (m/s)Spacecraft No. 02.67343e + 046.39409e + 02Figure 10
The three-component change curves of the velocity and displacement of spacecraft No. 0 in the case of Model 2.
(a)(b)It can be seen that the trajectory model of spacecraft No. 0 is consistent with the three-component change curve of velocity and displacement in the case of Model 2. Since the vertical axis in Figure10(a) has a large scale, the three-component position curve is approximately a straight line. Besides, in Figure 10(b), the three components of velocity change in different directions with time, resulting in changes in the direction and the velocity of the spacecraft. This conforms to the actual operation law of the spacecraft.


## 7. Conclusions

Aiming at passive detection satellites, three algorithms are proposed in this paper to estimate the orbit of the spacecraft in the ascent phase: the two-satellite point-by-point intersection positioning algorithm, the single-satellite positioning algorithm for estimating the orbit of the single spacecraft in the ascent phase, and the multiple-point intersection collaborative positioning algorithm for estimating the orbit of multiple spacecraft in the ascent phase.The three proposed algorithms are based on the least square method, nonlinear least square method, Monte Carlo algorithm, genetic algorithm, Runge-Kutta algorithm, point-by-point intersection method of spatial solid geometry, and coordinate system conversion method in image processing. The experimental analysis of real data verifies the effectiveness of the two-satellite point-by-point intersection algorithm and single-satellite positioning algorithm for estimating the orbit of the single spacecraft in the ascent phase. [32].



---

*Source: 6219698-2022-02-24.xml*