Initial Orbit Determination Method for Low Earth Orbit objects using Too-Short Arc based on Bistatic Radar

The problem of initial orbit determination (IOD) for low Earth orbit (LEO) objects using bistatic radar too-short arc (TSA) observations is addressed. For TSA observations, the traditional IOD methods suffer low accuracy. For LEO objects with stable attitude, the high order kinematic parameters can be obtained from the time derivatives of the radar echo phase. In this paper, we propose an analytical IOD method using bistatic radar TSA observations, which contain the position measurements (bistatic range, azimuth angle, and elevation angle) and the high order kinematic measurements (bistatic velocity, acceleration, and jerk). As the undetermined target state variables constitute a complex system of equations that can only be solved iteratively, we define an auxiliary coordinate system based on the bistatic geometry to help reduce the equations to one unary quartic equation. Further, we derive closed-form expressions of the orbital state. The performance of the proposed method is evaluated using linearization approximations. Numerical simulations are carried out for several typical LEO observation scenarios to demonstrate the performance of the proposed method.


I. INTRODUCTION
A S the development of space activities, the number of low Earth orbit (LEO) objects is increasing rapidly [1]. To maintain a catalogue of these LEO objects and detect new objects, the radar systems are widely used in space surveillance tasks. Because of the large number of objects and relatively limited sensors resources, a radar surveillance mode has been developed to improve the observing efficiency. Unlike in the traditional tracking mode, radar systems operating in the surveillance mode can produce a large number of too-short arc (TSA) observations, where the arc is usually several seconds [2]. When the detected object is new, the initial orbit determination (IOD) needs to be done to provide the basis for further data association and cataloguing [3], [4]. Using TSA observations, the classical IOD methods, such as the Laplace's, Gauss', or Double-r iteration, are unreliable and give poor results [5]. Therefore, it is highly important to develop a more accurate IOD method using TSAs.
The IOD problem has been studied extensively using TSA. The TSA measurements usually contain limited information to obtain a full orbit with valid accuracy, which is the intrinsic weakness of IOD problem [2]. Milani et al. developed a method using TSA angle and angle-rate measurements, and then the range and range rate were constrained to the admissible region [6]. The true orbit is obtained by searching parameters in the region. After identifying two TSAs from the same object, the orbit can be determined [7]. Based on the work of Milani, DeMars et al. proposed an IOD method using a discretization of the admissible region [8], and then developed a solution of IOD using Gaussian mixture models [9]. Ansalone and Hinagawa developed the method from optical data by exploiting genetic algorithms [10], [11]. They determined the optimal IOD solution by minimizing the observation errors. These methods apply to the angle-only measurements. As radars provide range information, the IOD method using two positions and the times, which is referred to as the Lambert problem, has been widely developed [12]- [14]. These methods are not applicable for TSAs. Tommei et al. developed the admissible region for radar data [15]. Considering the admissible region using radar observations, DeMars et al. applied it to the IOD from data association [16]. The methods solve the problem iteratively. The analytical method using three stations was proposed in [17]. However, the combination for three stations would increase the complexity of IOD process. The Herrick-Gibbs method used three position vectors and epochs to obtain the middle velocity vector. The method used the Taylor-series expansion to approximate the velocity [5]. Since the high accuracy doppler measurements can be obtained by the phase of radar echoes [18], IOD methods using doppler measurements were studied. Shang et al. proposed an IOD method using measurements from two stations in [19]. The method derived an analytical solution of IOD using the ranges, velocities, and accelerations from two stations simultaneously. Zhang et al. proposed a method using single-site higher order doppler measurements, where the radial acceleration and jerk are utilized [20]. Bistatic radars offer several advantages over monostatic radars, such as the configuration flexibility [21], [22]. Thus, developing an IOD method using bistatic radar is considered, and the crux of the problem is the complex relationship between the target state and doppler measurements.
In this paper, we propose a new IOD method based on bistatic radar TSA measurements. The method uses six measurements from a single arc observation to analytically determine the target state, where the bistatic range, azimuth and elevation angles are used for position determination, and the bistatic velocity, acceleration, and jerk are used for velocity determination. To solve the problem that the estimated parameters are coupling in a high order algebraic system, a new coordinate system is defined for decoupling, and finally the analytical expression is derived by coordinate transformation. The RMSE (root-mean-square error) is used to evaluate the IOD performance, and singular cases with low accuracy of velocity determination are considered. This paper is organized as follows. Section II introduces the bistatic observation model and the relationship between observations and the object state, then the IOD method is presented by solving the equations formed by the measurements and state, including position determination and velocity determination. Section III presents the linearized error analysis using the RMSE to evaluate the performance of IOD. Section IV presents Monte Carlo simulations and the comparison between the theoretical result and the simulation result. Observation scenarios of typical LEO objects are considered. Some singular cases are discussed to avoid bad IOD performance. Finally, Section V concludes the paper.

A. OBSERVATION MODEL
The orbital state is described in the Earth-centered Inertial (ECI) system. Let the position vector be r and the velocity vector be v. Let X = r T , v T T denote the state vector. As shown in Figure 1, the two stations are located at S t and S r , and the station position vectors are S t and S r , respectively. The target is located at P, and S t S r P forms the bistatic plane. To simplify the analysis, the Earth is regarded as an ideal sphere. In the absence of perturbations, the target motion is regarded as the two-body problem. For a single short arc observation, the station is assumed to be fixed, while the station velocity, acceleration, and jerk caused by the Earth's rotation are included.
For a single short arc, the bistatic radar produces six measurements, which are designated the observation vector Y = [ρ, az, el,ρ,ρ, ... ρ ] T , where ρ is the bistatic range. Based on the receiving station, the topocentric coordinate system (SEZ) can be established. Then, the angle measurements az and el are defined in the SEZ system.ρ,ρ, and ... ρ are the bistatic velocity, acceleration, and jerk, respectively.
The bistatic range ρ is the sum of the bistatic transmitterto-target range ρ t and the bistatic receiver-to-target range ρ r : The SEZ system is defined in [23] and shown in Figure 2, where the origin is located at the receiving station, the Saxis points the South, the E-axis points the East, and SEZ composes a right-hand coordinate system. The azimuth (az) is the angle measured from north, clockwise to the location beneath the object, and the elevation (el) is measured from the local horizon positive up to the object. Hence, the angle measurements are: The azimuth angle is usually reconciled to range from 0 to 360°, and the value of elevation angle usually ranges from 0 to 90°.
The bistatic velocityρ is defined as the sum of the two velocities of the target relative to the two stations, and the acceleration and jerk are obtained by differentiating kinematic measurements with respect to time: To establish the observation equations of the high order kinematic measurements, the method of obtaining the measurements is considered. The bistatic velocity, acceleration, and jerk can be estimated from the phase of radar echoes. Assuming that the radar system uses the single carrier frequency coherent pulse train, the transmitted coherent pulse train can be expressed as where N is the number of the pulses, T p is the pulse repetition period, f c is the carrier frequency, and p t is the waveform function of the pulse. The received baseband signal is: where w (t) is the additive white Gaussian noise, and τ n is the time delay. τ n is related to the bistatic range: where c is the speed of light. The relative ranges are modelled as a polynomial. Let t 0 denote the reference time and ∆t = t − t 0 denote the relative time. The relative range can be described as: In the case of short arc observations, the target motion is approximately a third-order polynomial. Then, the bistatic doppler measurements can be derived as: The maximum likelihood (ML) estimation is used to extract the coefficients in the equation. As the range is usually derived from the time delay, only the three high order kinematic measurements are obtained by the ML estimation. Let θ denote the estimated parameter vector. The ML estimator can be obtained as:θ where s r is the discrete-time signal sampled from the baseband received signal. The three high order kinematic measurements can be estimated from the estimation process. In practice, most of radars use the linear frequency modulated (LFM) pulse train for space surveillance. The errors of measurements are analyzed based on the results of references [24]. When the reference time t 0 is the middle time of the integration time interval, the CRLB of the high order kinematic measurements are derived in [25]. In fact, the measurement vector Y can be considered to be obtained from the instantaneous state r T 0 , v T 0 T at the reference time t 0 . For simplicity, the subscripts are omitted. As the bistatic high order kinematic measurements are the sum of those relative to both stations, without loss of generality, the terms of the target relative to the receiving station are considered. The velocity can be formulated in the ECI system:ρ where ρ r is the relative vector from the station to the target, and˙ S r is the velocity vector of receiving station caused by the Earth's rotation. The accelerationρ r can be generated by differentiatingρ r with respect to time: where v is the norm of v, and¨ r is the acceleration vector of the target in the two-body model.¨ S r is the acceleration vector of receiving station caused by the Earth's rotation. The motion vectors¨ r and ... r can be derived from the target state: where r is the norm of r and µ is the gravitational parameter of the Earth. Since the jerk vector of the target can be represented by the position vector r and the velocity vector v, the jerk ... ρ r can be expressed as follows: ...
... S r the jerk vector of receiving station caused by the Earth's rotation. Similarly, the terms relative to the transmitting station can also be formulated.

B. IOD PROCESS
The target position vector can be obtained using three position measurements, [ρ, az, el]. According to the relationship between the SEZ system and the measurements, a transformation is performed to yield The bistatic receiver-to-target range ρ r , according to the geometry of bistatic stations and the target, can be written as: where L is the baseline range, and θ is the angle between the vectors − − → S r P and − − → S r S t . Let e LOS represent the unit vector of the line-of-sight (LOS) direction, which is written as: Finally, the position state of the target in ECI can be expressed by: where the transformation matrices are provided in [23]. Note that M ECF←SEZ depends only on the latitude and longitude of the receiving station.
To get the velocity vector, the three high order kinematic measurements are required. The three equations in Equation (3) are used to solve the three velocity components in the ECI system. Compared with the monostatic radar, the bistatic radar produces a measured velocityρ that is combined by the velocities relative to both stations, and the velocity of the target relative to a single station is unknown. Therefore, the three velocity components in the ECI system are coupled in three equations of Equation (3), which is not feasible to directly calculate an analytical solution. Thus, an auxiliary coordinate systemXŶẐ is considered as shown in Figure 3, where the origin coincides with the target and the axesX andŶ lie in the bistatic plane. AxisŶ is collinear with the bisector of the bistatic angle, axisẐ is vertical to the bistatic plane, andXŶẐ composes a right-hand coordinate system. The unit vectors of this system are obtained: The transformation matrices between the ECI system and XŶẐ system are: The stations lie in the bistatic plane, which means that the position vectors between stations and the target can be described in this system: Based on the fact that eŶ is collinear with the bisector, according to the angular bisector properties, the coefficients of (20) have the following relationship: Because the position vector has been determined, the components in Equation (20) can be calculated. The velocity vector is described inXŶẐ system aŝ v = (vx, vŷ, vẑ) Combining the definitions of theXŶẐ system and the expression ofρ in Equations (3) and (10), vŷ can be derived as: Therefore, the components vx and vẑ are solutions to the two other equations. To solve this problem, the high order doppler measurements are described in theXŶẐ system.
coordinates of station velocities in theXŶẐ system. The bistatic velocityρ equalsρ r +ρ t , wherė After substituting Equation (24) into Equation (3), the acceleration can be obtained: To simplify the form, a substitution for the coefficients in former equations is done. The specific substitution is provided in Appendix A. Then the single velocity expressed in Equation (25) is simplified as: Considering that the bistatic accelerationρ equalsρ r +ρ t , a further substitution is formed: Thenρ can be expressed as: To separate the two unknowns, a further variable substitution is made: Then (28) can be expressed as: The two unknowns can be separated, and v ẑ can be written by vx: Similarly, substituting Equations (24) and (25) into Equation (3), the bistatic jerk ... ρ can be simplified and expressed as a binary cubic equation of vx and vẑ: .
Then v ẑ is used to replace vẑ and (32) is updated as: The specific substitution is provided in Appendix A. Let m 2 equal m 1 −ρ. Combining with (30), an equation of vx is obtained: According to the angular bisector properties shown in Equation (21), the coefficient of the highest term v 3 x in Equation (34), which is c 1 − c4k1 l1 , equals 0. Therefore, Equation (34) can be simplified into a quartic equation of vx: Additionally, the specific substitution is provided in Appendix A. To solve Equation (35), the method for finding roots for quartic equations, which is provided in Appendix B. After finding vx, vẑ can be derived using (31). Noticing that there exist two ambiguity solutions, the unique solution can be determined by combining with (32). Once the velocity in theXŶẐ system is expressed, the velocity vector v can be calculated by using the transformation matrix shown in (19).
Regardless of the measured error, the system of equations could produce four complex values of vx, which lead to four sets ofˆ v, and there is at least one real set that matches the true value of target state. All these sets of solutions meet the existing conditions, and some additional information could be used to solve this problem. The measured rate of angle can be used, but in order to pursue better orbit determination accuracy, it is not used for deriving the analytical expression.
In actual observation scenarios, due to the disturbance of the measured error, the equation (35) will produce at most four complex solutions, which means that the unique true value solution may have an non-zero imaginary part. In this case, the imaginary part of the unique root is discarded and the error between the real part and the true value are considered.
In summary, by using the bistatic radar observation including six measurements: bistatic range, azimuth, elevation, bistatic velocity, acceleration, and jerk, an analytical expression of target state is obtained.

III. ERROR ANALYSIS.
To analyze the performance of the proposed IOD method, the covariance matrix of the target state is required. As the transition from the target state to radar observation is nonlinear, the state estimate covariance does not have a compact analytical expression. The linearization method is used to approximately generate the covariance. Let C X denote the state covariance matrix and C Y denote the observation covariance matrix. According to the principle of the linearization method, C X can be derived as: The root-mean-square error is frequently used to express the differences between the observed values and the estimated values. The initial orbit can be regarded as the estimated state derived from the radar observations. Hence, we define the position root-mean-square error (PRMSE) and the velocity root-mean-square error (VRMSE) to represent the position accuracy and the velocity accuracy [26]. They can be generated as: To simplify the analysis, the covariance matrix is derived under the assumption of the two-body motion model, where the oblateness is neglected. Suppose that the measurements from the receiving station are disturbed only by random noise, which means that the system errors caused by sensors are removed. Let [σ ρ , σ az , σ el , σρ, σρ, σ ... ρ ] represent the standard deviations of the six measurements.
Considering the observation covariance matrix, some studies for LFM signal radar measurements have been done in [18] and [25]. The position measurements [ρ, az, el] and kinematic measurements [ρ,ρ, ... ρ ] are mutually independent. The matrix can be written as where C position is a diagonal matrix of measurement variances. The high order kinematic measurements are estimated by the maximum likelihood method using the high precision phase characteristic from radar echoes. The covariance matrix C doppler is the inverse of the Fisher information matrix.

A. POSITION ROOT-MEAN-SQUARE ERROR (PRMSE).
Let C pos ECI denote the covariance matrix of the position vector in ECI system, and C pos SEZ denote that in the SEZ system. As the transformation matrix is provided in [26], C pos ECI can be obtained: Since M ECI←SEZ is related to the Local Sidereal Time (LST), the station longitude, and the station latitude, it is orthogonal with no relationship with measurements. The PRMSE is the square root of the trace of C pos ECI , which equals the trace of C pos SEZ . Let Y pos = [ρ, az, el] T denote the observation vector that is used for position determination. Then C pos SEZ can be derived using linearization approach: Based on the chain rule, the partial derivatives can be written as: (41) where cos (θ) = − − → S r S t · e LOS /L. Let C θ represent cos (θ) for simplicity. Other terms can be derived using derivative rules. After substituting Equation (41) into Equation (40) and using Equation (37), the PRMSE can be obtained:

B. VELOCITY ROOT-MEAN-SQUARE ERROR (VRMSE).
Similarly, let C vel ECI denote the covariance matrix of the velocity vector in the ECI system. The VRMSE is the square root of the trace of C vel ECI , which is derived by: The velocity vector is determined using six measurements, and the transformation relationship is: Because the matrix is relative to position measurements, the partial derivatives can be derived: Note that the matrix is determined only by the position measurements Y pos . The partial derivatives of the matrix to kinematic measurements are 0, and that to position measurements is derived: As the unit vectors of the three axes in theXŶẐ coordinate system are provided in Equation (18), let − − → vec denote the vector ρr ρr + ρt ρt , the specific partial derivatives can be obtained as:: Because the position vector is known, the terms can be calculated by the derivation rule.
Consider the partial derivatives of the three velocity components in theXŶẐ coordinate system. The derivative of vŷ is simply calculated by Equation (23): The expressions of vx and vẑ are too complicated to directly analyze the partial derivations of measurements. Thus the binary system of equations composed of (30) and (32) is used. Taking the partial derivative of the two equations, two linear equations about ∂vx ∂ Y T and ∂vẑ ∂ Y T can be obtained: where ob i is the ith element of Y (i = 1, ..., 6). To simplify the expression, the coefficients of each term are substituted as follows: Since the intermediate variables in the process of velocity determination are given, the coefficients can be obtained by using the chain rule. Then by solving Equation (49), the partial derivatives can be derived as: It is observed through simulations that the denominator of the partial derivative greatly affects the error of velocity determination. Define the singular discriminant delta=p 2 q 1 − p 1 q 2 . When delta equals 0, the performance of velocity determination is bad and this case is regarded as singular. Note that when the orbit plane is coincident with the bistatic plane, the variances rẑ and vẑ are equal to zero, which results in the infinitesimal coefficients of q 2 and p 2 . Thus, the coplanar case is considered the known singular case. Now, by combining Equations (37) and (43), the VRMSE can be obtained.

IV. SIMULATION RESULTS
To verify the performance of the proposed IOD method, some numerical simulations are provided in this section. Similar to the derivation of the method described in the former text, the J2 perturbations are neglected for simplification. The simulation scenarios are based on low Earth orbit (LEO) satellites. The orbital elements are shown in Table 1. To verify the performance of the method for different orbits, Case 4 considers the orbit with a high efficiency of 0.7, which is usually regarded as a HEO object. For each case, the extrapolation from the time is performed based on the two-body model. The extrapolation step length is set as 20s. Therefore, several observations along the whole passing arc are obtained. Each observation epoch is regarded as the short arc reference time, and then the measurement vector  Table 2. The target is visible to both stations for all simulation scenarios.  In Section III, the RMSEs of the position and velocity vectors are derived, which represents the effect of the measurement error on the estimated state error. However, the procedure is based on a linearization approximation approach while the observation model is nonlinear.The model is approximately described using the first-order partial derivative of the state with respect to the observation. Monte Carlo experiments with 100,000 runs are performed to verify the reliability of the linearization. Each experiment can produce the deviations of position and velocity vectors, respectively, and the mean square error of these deviations can be approximated as the estimate error under the real observation model, as the number of experiments is large enough. In each figure, the legend 'MC' means the Monte Carlo simulation result, and the legend 'THEO' means the theoretical result. he radar parameters are provided in Table 3 for simulations. Then, the covariance matrices related to the parameters can be derived [25]:  (42), the PRMSE corresponds to the range error and two angle errors. Two observation cases are considered, one is that the target orbit passes over the receiving station right on the top, which means that the orbit coincides with the SZ plane, and the other is that the target passes over at a random angle. In Cases 1 and 2, the position measurements and the PRMSE results are provided in Figures 4 and Figure 5, respectively. The figures depict that the theoretical PRMSE is consistent with the Monte Carlo results. Analyzing the magnitude of each item of PRMSE in Equation (42) and comparing with the simulation results, it can be found that the PRMSE is mainly affected by ρ r and el. For a single pass, when the elevation reaches the maximum and the range ρ r reaches the minimum, the PRMSE is attained. We conclude that the best performance of position determination appears when the target passes over the receiving station.

B. VRMSE.
Since the velocity determination does not have a simple analytical expression and is derived from a quartic equation, the performance analysis is complicated. After simplifying the covariance matrix in Equation (43), it is influenced by many factors including position geometry and velocity components. Simulations of different aspects are considered. First, two non-singular observation cases are considered, and theoretical results with Monte Carlo results are compared. Then, as mentioned before, the denominator of partial derivatives in Equation (51) would cause singular cases when it nearly equals 0. Simulations are carried out for this case. Finally, by traversing the ascending node and changing the orbit inclination, and in all visible arcs, a contour map of the velocity error is provided as a reference to choose the better observation scenarios.

1) Nonsingular cases
The two Cases 2 and 3 have different inclinations of 90°a nd 75°, and the ascending node is set at 55°. Case 4 is the orbit with high eccentricity. The results are shown in Figures 6, 7, and 8, where the subfigures have the following representations: (a) As the sign to judge singular cases, the denominator in Equation (51) is given, which is denoted as delta; (b) the comparison between theoretical VRMSE and Monte Carlo VRMSE is done; (c) The errors of velocity components in ECI system are provided; (d) the comparison between theoretical total VRMSE and the partial VRMSE, when the measurements are only perturbed by the jerk and angle errors, is provided.
As depicted in Figs. (a), the arcs have no cases in which delta equals 0. Obviously, the error is inversely proportional to the absolute value of delta. Figs. (b) depict that the theoretical VRMSE values are consistent with the Monte Carlo results. Note that the minimum VRMSE appears when the target nearly passes over the receiving station, and the VRMSE can reach several meters per second when the elevation angle is larger than 25°. Combined with the analysis in PRMSE, it shows that the best performance of the proposed IOD method appears when the elevation angle reaches the largest value for a nonsingular pass. The results shown in Figs. (c) demonstrate that the linearization approximation is reliable. The velocity error components in the three directions of the ECI system are consistent with the simulations. Experiments show that the total VRMSE is mainly influenced by the partial derivative of the estimated velocity with respect to the jerk and the angles. Therefore, the comparison is provided in Figs. (d). It has been verified that for the whole arc, this partial VRMSE value is almost the same as the total VRMSE value, which means that the error of velocity determination is determined by the three measurements: az, el, and jerk.

2) Singular cases
An arc with the singular observation scenario, regarded as Case 5, is used for simulation, where the discriminant delta has zero points. The orbit inclination is 90°and the ascending node is 62°. As shown in Figure 9(a), delta has four zero points. Since the denominator equal to zero makes no practical sense, the near-zero interval is observed. Figure 9(b) depicts the results of VRMSE containing singular cases. When delta is nearly 0 points, the nonlinear effects cannot be ignored. Although the two results are not consistent, the Monte Carlo results show that the proposed solution of velocity does have poor performance when delta equals 0, and the lowest accuracy of velocity determination could reach kilometers per second. This is intolerable for the absolute velocity of LEO targets. Therefore, in actual scenarios, a larger absolute magnitude discriminant delta is the premise for better IOD accuracy.

3) Scenario traversal
The discriminant delta has a great influence on the VRMSE and is jointly determined by the target position and velocity. Simulations for some common scenarios are carried out to help avoid singular cases. By traversing the ascending node, the target appears at the position in the entire visible range. By changing the orbit inclination, the effect of different velocity directions on the VRMSE can also be represented.
The first cases have a constant orbit inclination of 90°, and the results are shown in Figure 10. The relationship between the VRMSE and the ascending node Omega is provided. The numbers indicated in the figure are VRMSE contours, where the unit is meters per second. The outermost line distinguishes the visibility. Because the station longitudes are 60°and 70°, the traversal range of Omega is from 40°to 70°. The second cases have a constant orbit inclination of 75°and are shown in Figure 11. The outermost line is at an angle to the central axis, and the whole graph is symmetric about the center point. The largest VRMSE, which reaches kilometers per second, appears when the target is nearly located in the plane, about which the two station positions are symmetric.
The best VRMSE appears when the target passes over the top of one station and can reach several meters per second. As a consequence, for a single pass, the optimal observation scenario, in which the error of the IOD method can satisfy the needs for the initial orbit, is when the target passes over the top of stations, and the target is not located near the central vertical plane of the two-station baseline.

V. CONCLUSION
In this paper, we present a new initial orbit determination method using a single too-short-arc observation based on bistatic radar. By using measurements including the bistatic range, azimuth and elevation angles, the bistatic velocity, acceleration, and jerk, a closed-form solution for the space object state is derived, which contains the position and velocity vectors. As the relationship between the state and the bistatic high order kinematic measurements is more complex than that of a monostatic radar, the undetermined parameters are coupling in the observation equations. By defining an auxiliary coordinate system based on bistatic geometry, the   parameters are separated and the equations are transformed into a binary system that the analytical solution can be derived. Finally, we use coordinate transformations to obtain the closed-form expressions of the orbital state. We evaluate the performance of the proposed method by RMSE, and the position and velocity errors are represented by the linearization approach. Simulations for theoretical and Monte Carlo results are performed for verification. Additionally, the singular cases with poor IOD performance are considered, where the errors are not acceptable for normal LEO objects.
To avoid singular cases, simulations of some LEO scenarios are provided as a reference for better IOD performance.
The method uses the measurements obtained from a single too-short arc observation to derive the complete orbital state vector. No additional information or multiple observations are required. The process provides the setup for new detected targets that must be catalogued by the space surveillance system. When more observations are gathered, the data association or precise orbit determination can be performed based on the derived initial orbit. .

APPENDIX A
For simplifying the expression in Equation (26), the specifitic substitutions are provided as follows: For simplifying the expression in Equation (32), the details of the variable substitution are provided in the text: (60) where¨ S r =[a rx , a ry , a rz ] T and¨ S t =[a tx , a ty , a tz ] T are the coordinates of station accelerations in theXŶẐ system.
And the further substitutions from (32) to (33) are:

APPENDIX B
There are methods to solve the quartic equations, and we provide the procedure of Ferrari's solution for finding vx in Equation (35). Given a general form of a quartic equation as We can get the roots following this process: