Abstract
The differential geometric filter is implemented to estimate the absolute and relative positions of the spacecraft in a formation. The extended Kalman Filter is also implemented as a benchmark for the differential geometric estimation. Only relative positions between the spacecraft are measured. Relative positions are measured using a wireless local positioning system (WLPS) installed in all spacecraft. Two different scenarios are studied: (1) the observations include WLPS measurements only and (2) the observations include WLPS measurements in addition to measurements for the absolute position of one spacecraft made by a radar that takes measurements from the earth’s surface. Results show that the differential geometric estimation has better stability performance and a faster convergence rate compared to the extended Kalman filter.
1. Introduction
The relative and absolute positions estimation of spacecraft formations is a fundamental task in many space missions. Relative position estimation plays an important role in Spacecraft Formation Flying (SFF) missions, a subject that has been the focus of many researches during last decade. Some SFF missions require that multiple spacecraft, in different orbits, communicate without interruption, for example, satellites in the CITRIS-COSMIC system are required to communicate with each other to monitor the ionospheric irregularities [1]. Relative positions between satellites, such as Cluster and Cluster-II satellites launched by the European Space Agency, are estimated and controlled to support many collaborative tasks where satellites are required to maintain a specific formation in a continuous manner within the mission period [2].
Several relative position estimation methods have been developed. A GPS-like technology that can be applied to SFF have been introduced [3, 4], where each spacecraft is equipped with a communication system to localize other spacecraft in the formation. This system provides high precision estimates for relative positions. Yet, complex hardware is required on each spacecraft. This system does not stand alone and requires a GPS system.
A relative position and attitude estimation through a Vision-Based Navigation System (VISNAV) has been addressed extensively in the literature [5–8]. The VISNAV consists of a Position Sensing Diode (PSD) sensor and an array of LED source beacons installed on a plate. The PSD detects the energy of the light source emitted by the beacons. Then the PSD generates a current flow to four terminals that are installed on the PSD plate. At that point, the centroid of the current flow is determined to measure the direction of the energy source. The covariance study shows that the performance could be affected by the relative distance between spacecraft. The accuracy factors of the estimation depend on the number of beacons installed on the spacecraft. Three or more beacons are suggested to ensure the observability of the system [5].
Psiaki proved the feasibility of estimating the orbits of two spacecraft through measuring their relative distance, azimuth, and elevation angles [9]. The orbital elements estimation accuracy depends on the relative distances between the spacecraft. Only two spacecraft formations were considered. As shown in [9], the relative position between spacecraft impacts the estimation performance. In addition, the system becomes unobservable in a few cases, such as zero inclination.
The extended Kalman filter (EKF) has been widely implemented into the nonlinear system to estimate the relative and absolute positions in SFF [5, 10–12]. The system nonlinearity causes the EKF estimation to be sensitive to initial condition errors, and the linearization step in the EKF may result in estimation instability. The Unscented Kalman Filter (UKF) has been implemented for relative attitude and position estimation in SFF [13, 14]. The UKF has a better robustness to initial condition errors compared to the EKF. However, its computational complexity is higher than the standard EKF. If compared to the standard EKF, the UKF has a faster convergence rate, but the computational time required for each update is longer than that of the EKF [15]. The Batch filter [16], an offline filtering method, has been implemented in [9] for absolute positions estimation using the relative positions measurements. Both studies in [9, 12] show that the configuration of the SFF affects the stability and the accuracy of estimation.
The differential geometric (DG) theory has been widely implemented in the control and guidance research area [17, 18]. In calculus, the differential geometry represents the linear approximation of a smooth curve [19]. References [20, 21] introduced the differential geometric approach into the missile tracking applications. In [20], the system model is transformed from the arc length domain into the time domain; then the guidance law is applied. Both papers have shown the feasibility of implementing the DG for missile guidance.
The differential geometric approach for nonlinear systems has been extended recently for the nonlinear estimation purpose [22]. The DG estimation avoids the linearization step in the EKF. In DG estimation, the nonlinear dynamics are mapped to a linear domain, where a linear estimator can be implemented. Then, the inverse transformation is applied to the estimator [22]. Reference [22] details how the DG estimation can be implemented when the number of available measurements is too few to handle all the system nonlinearities and shows that both the optimal control and the linear filtering approaches can be applied in the mapped linear domain.
In this paper, the DG filter is implemented in estimating the spacecraft relative and absolute positions in formation. While most of the developments in the literature focus on the estimation of relative positions [3, 5], this paper estimates relative and absolute positions of all spacecraft. The differential geometric (DG) estimation and the extended Kalman filter (EKF) are implemented and compared for estimation using wireless local positioning system (WLPS) measurements. The WLPS installed on each spacecraft enables that spacecraft to determine the relative positions of other spacecraft located in its coverage area via time-of-arrival (TOA) and directional-of-arrival (DOA) measurements [12, 23, 24].
In this paper, two scenarios are considered: (1) observations include WLPS measurements only and (2) observations include WLPS measurements in addition to the absolute position of one spacecraft measured by radar systems installed on the earth. Section 2 presents an overview on the WLPS. Section 3 presents the orbit dynamics model. The derivation of the DG Estimation equations for the problem of relative and absolute positions estimation are presented in Section 4. Section 5 presents the EKF implementation. Section 6 discusses the simulation results and presents a complexity analysis that compares the computational costs between the DG filter and the EKF.
2. Wireless Local Positioning System (WLPS)
The WLPS consists of two basic components [12, 23, 24]: A dynamic base station (DBS) and a transponder (TRX). Each DBS is capable of localizing TRXs located in its coverage area via TOA and DOA measurements, as shown in Figure 1. The DBS periodically broadcasts an identification (IDR) signal once every ID Request Repetition Time (IRT) as shown in Figure 2. A TRX within the DBS coverage area receives the IDR signal and transmits a response signal that includes its own ID back to the DBS within the IRT period. The ID of each TRX allows the DBS to distinguish one TRX from another. It also allows the DBS to easily track multiple TRXs located in its coverage area.
As shown in Figure 2, the range of a TRX is measured by comparing the TOA of the signal from the TRX at the DBS receiver and the time of the transmission of the signal from the DBS transmitter. Here, it is assumed that the processing time delay in the TRX is known. The processing time estimate can be included in the signal packet transmitted from the TRX to the DBS in order to allow the DBS to correctly measure the range. The DBS, equipped with antenna arrays, allows DOA estimation and beamforming. In addition, beamforming enhances the performance of the DBS by reducing the interference effects [25]. The DOA is measured by each spacecraft relative to its body fixed coordinate system. In this paper, we assume that the attitudes of all spacecraft are known. Hence, all the DOA measurements can be computed relative to a fixed reference frame.
Thus, a WLPS allows single node positioning. In other words, each node equipped with a DBS can independently localize the TRXs located in its coverage area and its field of view (FOV). Now, if all spacecraft are equipped with both a DBS and TRX, each spacecraft can find the position of other spacecraft located in its FOV and coverage area. The position information across multiple spacecraft can be fused to improve the localization performance [26, 27].
3. System Model
In this section, the state and measurement models are derived. Two sets of models are derived for two scenarios: (1) only WLPS measurements are available and (2) WLPS measurements and those taken by a tracking system (such as a radar installed on the ground) are available. For the first scenario, a four-spacecraft formation is considered. Each spacecraft is equipped with both a DBS and a TRX. For the second scenario, a two-spacecraft formation is considered; one spacecraft is equipped with a DBS and the other spacecraft is equipped with a TRX. The radar measures the absolute range, the absolute azimuth, and elevation angles of one spacecraft with respect to a ground station. Since the position of the ground station is known, we assume the absolute position measurement is expressed with respect to earth center [28].
The estimated states are considered to be the spacecrafts' absolute and relative positions and their velocities. If the spacecrafts' orientations are known, the WLPS measurements can be expressed in the inertial reference frame.
3.1. State Model
We assume the spacecraft's absolute position in the Earth Centered Inertial (ECI) frame is , the velocity vector is , and the acceleration vector is . The dynamic model of each spacecraft with respect to the earth center, can be expressed using a two-body model, which is given as [29] where, is the Earth gravitational constant, and is a zero-mean gaussian acceleration uncertainty, and ; here, refers to the statistical expectation operation. The vectors and are the absolute position and acceleration of the spacecraft, respectively, and denotes the absolute distance of the spacecraft from the center of the Earth.
3.2. WLPS Measurements Model
In the first scenario, we incorporate only WLPS measurements for localization. In Figure 3, and represent the th and the th spacecraft absolute position vectors expressed in an ECI frame. The relative position vector of the th spacecraft observed by the th spacecraft equipped with DBS, expressed in ECI frame, is
In general, there are a maximum of 12 possible measurement sets for a four-spacecraft configuration. Here, we consider only four measurement sets in the paper: , , , and . If the orientation of the spacecraft (-reference frame in Figure 3) is known, the relative position vector, , expressed in the th spacecraft reference frame would be where is the Directional Cosine Matrix (also known as Attitude Matrix [16]) of the th spacecraft relative to the ECI frame. The relative distance between spacecraft and is a scalar linear function of TOA, , that is, where is the speed of light, and refers to the magnitude of vector. The DOA between the two spacecraft is given by: Here, is the relative azimuth angle, and is the relative elevation angle between the two spacecraft and . Using (6) and (7), we can express the WLPS measurement between the th and the th spacecraft as where denotes WLPS measurement noise, which is assumed to be zero-mean gaussian with .
Using (A.2) in Appendix A, the WLPS measurement vector in Cartesian coordinates is The new measurement noise is a function of , and are defined in (8).
3.3. Radar System Model
In the second scenario, an additional radar system that measures the spacecraft range, azimuth, and elevation angles is included. It is assumed that the th spacecraft that is observed by the radar is also equipped with a DBS. The absolute position measurement made by the radar is where is radar system noise, which is assumed to be zero-mean gaussian, and . The is the absolute distance of the th spacecraft from the radar, is the azimuth angle, and is the elevation angle of th spacecraft with respect to a radar-fixed coordinate system. The radar position is known in the ECI frame. So, using radar measurements, it is straightforward to calculate the absolute distance of the th spacecraft from the earth's center, and the azimuth, , and elevation, angles with respect to the ECI frame. In this preliminary analysis, we will assume for simplicity that we measure directly , , and . Azimuth and elevation angles are related to the coordinates of the spacecraft through where , , and are defined in (3).
3.4. Measurement Models for Simulation
In the first scenario, only the linear filtering method presented in Section 4 will be implemented. Four-spacecraft formation is required, as discussed in Section 4.3. Thus, the measurement vector, is where is the measurement noise vector. The corresponding measurement noise covariance is given as
The which denotes the relative position measurement noise covariance in the ECI frame in (14) is given as where is the noise covariance of the relative position measurement, expressed in cartesian coordinates, and is the attitude matrix [16]. is given as a function of [28]. The transformation from to is shown in Appendix B.
For the second scenario, both the linear filtering method and pole placement, which are presented in Section 4, will be implemented for estimation purposes. Only two-spacecraft formation is considered. The measurement vector is where . The corresponding measurement covariance matrix is where , which the “diag” represents the diagonal matrix, and , , and are the measurement noise variances for absolute distance, azimuth, and elevation angle, respectively.
4. Differential Geometric Filter
DG estimation was recently extended to deal with nonlinear dynamic systems with fewer measurements than required to handle all nonlinearities [22]. DG estimation avoids the linearization step encountered in the EKF. In the DG estimation, the nonlinear system is transformed into a linear system. The transformation requires the state vector, , to be expressed as a function of the measurements vector, , that is .
The DG estimation developement is explained in [22] and is briefed here for completness of presentation. The system model for DG estimation is constructed by defining a new state vector such that the measurement vector can be expressed as a linear function of Here, , , , and are linear matrices, is the process noise, and is the measurement noise. The state equation (18) includes a linear system, , and a nonlinear input function, .
When the available measurements are not enough to handle the nonlinearities in the system, a set of pseudomeasurements, , is required. The pseudo-measurements can be obtained by taking the derivatives of , that is, .
The pseudo-errors are defined as
Then, the estimated state vector, , that includes the pseudomeasurements is
Letting be the total number of measurements, (21) shows that the size of the estimated state vector, , is twice the size of the measurement vector, . In addition, it can be shown that the time rate of change of the pseudoerrors is [22] where is the gain matrix which will be presented in the next section, and and denote th to th row of the and matrices, respectively.
In [22], two gain computation methods are presented, which are the pole placement method and Kalman Filter (or Linear Filtering) method.
4.1. Linear Filtering Gain
The state (18) is linear. Thus, the Kalman Filter theory can be applied. Note that the input function, , is not used in the updating process in both DG estimation and EKF. The DG estimation gain matrix, , is determined by where is the state covariance matrix corresponding to the vector, and is the measurement noise covariance matrix. Then, the state estimates, , and its covariance are updated as follows: where is the identity matrix, the superscript − denotes preupdate estimates, and the superscript + denotes postupdate estimates. The denotes any estimated pseudomeasurements vectors and the estimated state vectors that are not observed by measurements. denotes th to th row of the gain matrix, .
The DG estimation filtering process works as follows [22]. The gain matrix, , is determined using (23). Then the estimated states and state covariance matrix are updated using (24). The updated states and the pseudoerrors are propagated by (18) and (22) with their respective input function. The state covariance matrix is propagated using (25) where is the process noise covariance, which is defined as .
Both the and matrices in (18) and (25) are linear time invariant matrices; thus, (25) can be expressed in time discrete representation, which is given as where, subscript and denote the current and next time step, respectively, is the state transition matrix, and is the discrete-time process noise covariance. Both and are calculated as follows [30]. Let where is the time interval between two successive measurements, is the exponential matrix of matrix, and is the matrix representation of the system dynamic model. For the nonlinear dynamic model case, for example, (1), is the first order Taylor series expansion (or Jacobian matrix) of the dynamic model. Both and are given as
4.2. Pole Placement Method
Another filtering method presented in [22] is the pole placement method. The pole placement method has been widely used in controls and estimation. Unlike the Kalman filter, the gain, , in pole placement, does not change dynamically. Pole placement also requires full system observability. However, the pole placement method guarantees the linearity of the measurement model, and hence guarantees stability.
The pole placement gain matrix is computed as follows. Given that the matrix has dimensional and canonical form as the following: and matrix has dimensional. The system is considered to be fully observable if the rank of matrix is full rank [31] as the following:
If the full rank observability condition in (30) is fulfilled, the gain matrix is calculated using the following equation: where is the eigenvalue vector, and eig represents the eigenvalue.
If a set of eigenvalues are given, the gain matrix can be computed using the Ackermann's formula [16]. The selection of the eigenvalues impacts both stability and accuracy of estimation. Therefore, negative eigenvalues are always selected to ensure the stability of the estimation process.
Here, our goal is to set up the DG estimation filter for the two scenarios introduced in Section 3. For the first scenario, only the WLPS measurement is available while we intend to estimate the spacecraft absolute position. Thus, it is required to derive the expression of absolute position of spacecraft in terms of relative position. For the second scenario, an additional absolute position measurement of spacecraft is included.
4.3. Scenario One—Relative Position Estimation with WLPS-Only Measurement
The absolute distance of a spacecraft from the Earth's center is where and denote absolute and relative positions, respectively, defined in (2) to (4).
The orbital equation of spacecraft is [29]
In (33), represents the distance of spacecraft to the earth's center, is the semimajor axis, is the eccentricity, and is the true anomaly. In this scenario, it is assumed that all spacecraft in the formation have the same semimajor axis, eccentricity, and true anomaly. Hence, (33) shows that the distances from all spacecraft to the earth's center would be the same. In this case, , leads to
If we consider a four-spacecraft configuration, then (34) can provide a closed form solution for absolute positions using relative position measurements. For example, the absolute position of spacecraft 1 can be obtained from Here, , , and and , , and are the absolute and relative position elements introduced in (2), respectively. The right hand side of (35) is the square of relative distances between the spacecraft. Similarly, we can determine the absolute positions of the other spacecraft. Note that the solution in (35) requires the computation of inverse of the relative position measurements matrix. Therefore, the closed form solution needs an accurate estimation or unbiased relative position measurements.
Consider the four sets of WLPS measurements mentioned in Section 3 where the total number of measurements, , in (13) is 12. However, there are 24 states to be estimated; thus, at least 12 pseudomeasurements are required. The pseudomeasurements are defined as the relative velocities. Using the dynamic motion of spacecraft defined in (1), the relative accelerations between spacecraft are
When all spacecraft in the formation have the same semimajor axis, , eccentricity, , and true anomaly, , at all times, (36) can be simplified to.
The nonlinear input function, , in (18) is defined as a function of measurements, , and pseudomeasurement, . The only nonlinear function of the dynamic system in this scenario is the relative acceleration introduced in (37). Therefore, the input function consists of the relative acceleration between the spacecraft, which is where is the pseudomeasurements vector defined in (20), which can be written as .
Equation (38) depicts that the absolute range of spacecraft that is required for the relative dynamic model. Equation (35) determines the absolute position using relative position measurements. However, (35) is vulnerable to the error in , which could result in a large estimation error in the absolute range, . Then, we consider two additional pseudomeasurements which are the absolute range and eccentric anomaly, . The rate of change of absolute range is [29] where and were defined in (33). The rate of change of eccentric anomaly is:
Substituting (40) into (39), we obtain the expression of in terms of eccentric anomaly and range only. Adding (39) and (40) into the input function, in (38), we obtain
Then, the estimate state vector, is
The linear matrices , , and in (18) associated with state models in (41) and (42) correspond to
For the WLPS-only measurement, the matrix in (19) associated with the measurement model in (13) is
4.4. Scenario Two—Relative Position Estimation with WLPS and Radar Measurement
Consider a radar tracks one of the spacecraft in the two-spacecraft formation, for example, th spacecraft where represents either spacecraft 1 or spacecraft 2. The absolute distance to the spacecraft is measured at all times. Then, the range and eccentric anomaly pseudomeasurements in (42) can be ommitted. Therefore, the state vector for the second scenario becomes where, is the absolute range, is the absolute azimuth angle, is the absolute elevation angle, and , , are their time derivatives, respectively. Based on (47), the nonlinear function are the relative acceleration and the acceleration in terms of , , and . Then, the modified input function, is where , , and are the spacecraft polar accelerations, which are the second order time derivatives of the absolute range, absolute azimuth, and elevation shown in (10)-(11) respectively. Their corresponding equations are where the first order time derivation of polar coordinate and cartesian coordinate are presented in Appendix A.
In this case, the , , , and matrices for DG estimation's state and measurement models in (18) and (19), are The process noises in the matrix in (54) are assumed along the acceleration axes only.
5. Implementation of Extended Kalman Filter
The standard EKF has been implemented for both scenarios. The state model, (1), and the measurement models, (13) and (16), correspond to the first and second scenarios, respectively. In [12], the EKF has been derived for the WLPS-measurements expressed in polar coordinates.
For -number of spacecraft in the formation, the state estimate vector for EKF is
For the first scenario, the measurement model is:
For the second scenario, the linearized measurement model is where is the jacobian of the three radar measurement components introduced in (10) with respect to the th spacecraft absolute positions and velocities [28], which is given by where and , , and are defined in (10).
Using the measurement noise covariance obtained in (14) and (17) together with the linearized state and measurement models, the EKF for relative position estimation is implemented.
6. Simulation Results and Discussions
Simulations are conducted to study and compare the accuracy and convergence performance between differential geometric estimation and EKF. Two differential geometric estimation methods discussed in Section 4 are implemented in the simulation, which are the pole placement and the linear filter. The accuracy performance measured is their root mean square error (RMSE). Here, the mean square error (MSE) is the average of the square of estimation error over , , and axes of spacecraft position. Thus, where is the estimated absolute position of the th spacecraft, and is the true absolute position of the th spacecraft.
The convergence rate is determined by the amount of time required by the estimator's RMSE to fall within a RMSE threshold, . The is determined as the average RMSE at the estimator's RMSE steady state. In addition, the complexity of DG Estimation and EKF are compared through the number of multiplications required by each estimator in an iteration.
In the simulation, all spacecraft are in low earth orbit (LEO), with a semimajor-axis of 7000 km, the eccentricity of 0.04, and 0 degree of initial true anomaly. The argument of perigee (ARGPER), the Right ascension of the ascending node (RAAN) and the inclination (INC) of each spacecraft are specified in each scenario in order to meet the formation configuration requirements.
The initial attitude of each spacecraft is constructed using a euler angle with 3-1-3 orientation at periapsis [32]. The first rotation is the spacecraft RAAN angle, followed by a second rotation (inclination angle). Then, the third rotation angle is 90 degrees for spacecraft 1, 3, and 4 and −90 degrees for spacecraft 2. To ensure that one of the spacecraft reference frame axes points toward the center of the Earth at all times, all spacecraft are set to rotate at a constant angular velocity, which is, . Here, is the angular velocity vector expressed in the th spacecraft reference frame. Then, the angular velocity of spacecraft expressed in ECI is . The mean motion is ; where is the earth grativational constant, and is the semimajor axis.
The simulation assumptions are as follows: (i) process and measurement noises for all spacecraft are zero-mean Gaussian; (ii) Initial conditions of spacecraft are well known; (iii) The attitude of each spacecraft is well known; (iv) All spacecraft are equipped with a DBS and TRX; thus, each one of them can localize others; (v) The reference frame of the DBS is aligned with the corresponding spacecraft's attitude reference frame; (vi) The radar system measurement is available at all times for the second scenario; (vii) no signal transmission delay; (viii) The process noise variance is km/s2 along each axis for all scenarios.
6.1. Observability Analysis
The pole placement method presented in Section 4.2 requires the system to be fully observable. In this section, the system observability analysis of these two scenarios is presented. For simplicity, we consider a two-spacecraft formation case for the first scenario in this observability analysis. The estimated states for DG estimation, after the transformation, is . The matrix for the dynamic model in (18) is
We assume that only relative position is measured in cartesian coordinate. The matrix is given as:
The rank of the observability matrix using (30) is three, which is not full rank. Therefore, the pole placement method presented in Section 4.2 is not applicable.
In the second scenario, a radar measures one spacecraft's absolute position (see (16) and (47)). Both the and matrix are given in the (52) and (55), respectively. The rank of the observability matrix is twelve, which is full rank. Then, the pole placement method will be implemented in this scenario and compared with other estimation methods.
6.2. Scenario One
Consider only WLPS measurements are available in the four-spacecraft formation. The measurements vector, state vector, and its corresponding nonlinear model and matrices are given in (13), (41) to (45), respectively. Three formation configurations that represent the short range (≃0.25 km), medium range (≃60 km), and long range (≃1200 km) are considered. The spacecraft formation configurations for different formation size are shown in Table 1. We study the impact of formation size on both EKF and DG estimation performance.
The measurement noise of WLPS is assumed to be 1 meter in distance (computed based on TOA measurement) and 0.001 degree in DOA. Then, the measurement noise covariance matrix is constructed using (14).
The process noise parameter, in (1) and (18), is assumed to be on each absolute and relative acceleration axis for both EKF and DG Estimation. The initial states standard deviations for EKF are 1 km and km/sec in absolute position and velocity, respectively. The initial states standard deviations for DG Estimation are 1 km in relative position, 0.5 km/sec in relative velocity, 2 km in absolute range, 0.001 degree in eccentric anomaly, and 0.001 km/sec for the pseudoerrors.
6.3. Discussion—Scenario One
Table 2 shows the mean RMSE in absolute position for DG estimation and EKF with respect to different mean relative distances. The results show that the performance of DG estimation improves as the mean relative distance increases when the measurements consist of WLPS only. In Figure 4, the signal-to-noise ratio (SNR) between estimation error and relative distance decreases when the spacecraft in the formation get closer (or 1/SNR increases), in which case the errors in the transformation from relative positions to absolute positions using (35) become significant.
Table 2 depicts that the EKF has a better accuracy compared to DG Estimation, when only WLPS measurements are used. The EKF does not require the inverse transformation (34). However, when the relative distance between spacecraft increases, the effect of DOA estimation error on the positioning error increases (see Figure 5). Accordingly, the performance of EKF drops slightly as the mean relative distance increases.
If only one DBS is available in a three-spacecraft configuration [12], the EKF would be unstable when the mean relative distance between spacecraft falls within a certain range (≤10 km). EKF stability can be improved by increasing the total number of DBSs installed in the formation, which is shown in Table 2. In addition, increasing the number of DBSs improves the estimation accuracy.
6.4. Scenario Two
Consider WLPS and radar measurements are available. Here, the two-spacecraft formation with only the medium range formation configuration is considered. The radar is measuring the absolute range, azimuth and elevation angle of spacecraft 1 at all time. The configuration of both spacecraft is similar to the configuration of S/C 1 and S/C 2 shown in Table 1. The measurement vector, state vector, and its corresponding nonlinear model and matrices are given in (16), (47) to (55), respectively. The performance of the EKF and the DG estimation for a given range of measurement noise is studied. The measurement noises of WLPS and radar are depicted in Table 3. The measurement noise covariance matrix is constructed using (17).
The process noise parameter, , is assumed to be for all axes in EKF. For DG Estimation, the standard deviation of process noises for the relative position acceleration is , the absolute distance's acceleration is , and the absolute azimuth and elevation angles' acceleration is . The initial states standard deviations for EKF are 2 km in absolute position and 0.1 km/sec in absolute velocity. The initial states standard deviations for DG Estimation are 0.1 km in relative position, 0.05 km/sec in relative velocity, 1 km in absolute range, and 0.05 degrees in azimuth and elevation angles. In addition, the initial variance for the rate of change of the range is 0.1 km/sec, and the rate of change of the absolute azimuth and elevation angles is 0.001 degree.
The simulation for each set of corresponding measurement noise is repeated for 10 times to compare the consistency of results produced by both filters. Then the maximum, minimum, and mean RMS of each set of measurement noise for both EKF and DG Estimation are compared.
6.4.1. Linear Filtering Gain
First, the DG Estimation is compared with the EKF performance using the gain computation presented in Section 4.1. Figures 6 and 7 show the RMSE of DG Estimation and EKF for measurement noise standard deviation of 0.001 km in relative range, 0.001 degree in relative azimuth and elevation angle, 0.01 km in absolute range, and 0.01 degree in absolute azimuth and elevation angle. Figure 6 shows that the DG Estimation converges faster than EKF, because the DG Estimation avoids the linearization steps in the estimation.
However, Figure 7 shows that EKF is more accurate than DG Estimation. The nonlinear dynamic model in DG estimation is expressed in terms of measurements and pseudomeasurements (see (41) and (48)). While the linearization step is eliminated in the DG estimator, the propagation step is dependent on the measurements noises. On the other hand, the nonlinear dynamic model of EKF is expressed in terms of the estimated state vector with additional linearization steps. The EKF is able to achieve a better accuracy if the estimation converges.
Figure 8 compares the minimum and maximum RMSE of DG estimation and EKF for a given set of measurement noise standard deviation. The DG gain is computed using the method shown in Section 4.1. With additional radar measurement, the DG estimation accuracy is improved. In addition, the minimum and maximum RMSE of DG Estimation are very close compared to EKF because the DG estimation does not require linearization steps. Therefore, a stable estimation can be always achieved.
Figure 8 also shows that the EKF has better accuracy performance than the DG estimation. However, the stability of the EKF is not always guaranteed. In Figure 8, we observe large RMSEs in the EKF estimation at higher measurement noise levels. The large RMSE is caused by the unstable estimation of EKF, which is due to the linearization of the nonlinear radar measurement and state models in the EKF algorithm. On the other hand, the DG Estimation does not require any linearization step.
6.4.2. Pole Placement Method
Figures 9 and 10 compare the RMSE performance between the DG Estimation and the EKF. The EKF's configuration is the same as in the previous case. On the other hand, the pole placement method is as follows. The poles for the relative position are , for the relative velocity vector are , and for the absolute position and velocity in polar coordinates are and , respectively. The DG estimation gain is computed using the Ackermann's formula [16].
Figure 9 shows that when the EKF estimation is stable, the EKF has better accuracy performance than the DG estimation. However, the linear system characteristic of the DG estimation allows it to converge faster than the EKF. Figure 10 shows that the EKF may suffer instability due to the linearization steps in the algorithm. In Figure 10, the EKF is able to converge to lower RMSE than the DG estimation, but the RMSE diverges after a certain time period. However, the DG estimation does not suffer any instability. Therefore, the pole placement method in DG estimation does guarantee the fast convergence and stability of estimation algorithm.
6.5. Complexity Analysis
In [33], a formalism of function for complexity computational of matrices operation has been introduced, which are: where, , , and are the dimensions of matrices that perform these operations. Here, is the number of states in DG estimation, is the number of states in EKF, is the number of measurements, and and are the number of process noise in DG estimation and EKF, respectively.
We assume the computational complexity of the jacobian matrices is in the order of . In addition, the propagation of the state models is considered to have the complexity of in both filters. The EKF requires an order of for measurement model computation, while the DG estimation only requires an order of computation complexity. Therefore, the EKF requires higher computational complexity to compute the measurement.
The EKF requires the order of and computation complexity for gain matrix. There are two different computation complexities in DG estimation for gain matrix computation. For the pole placement method, there is no computational complexity because the gain matrix is predetermined. For the linear filtering method, the gain computational complexity is in the order of and . Both EKF and DG estimation may have same computation complexity if all states are observed. However, the DG estimation has higher computational complexity if not all states are observed ().
The covariance update and propagation's computational complexity for the EKF are in the order of , , and . There is no computational complexity for covariance update and propagation for DG estimation if the pole placement method is used. However, for the linear filtering method, the computational complexities are in the order of , , and .
Therefore, for fully observable case, the DG estimation has lower computational complexity. In addition, among the two DG estimation approaches implemented in this paper, the pole placement method has lower computational complexity compared to the linear filtering method. But the EKF has lower computation complexity if not all states are fully observed.
7. Conclusion
This paper presents an implementation of the differential geometric estimation in relative and absolute positions estimation for spacecraft formations. In the transformed linear domain, two approaches have been implemented in this paper: the pole placement and the linear filtering approaches. Simulations are conducted to compare the estimation performance of the differential geometric estimation and the extended Kalman filter. The results show that the differential geometric estimation has a faster convergence rate, and has better stability, compared to the extended Kalman filter. The complexity analysis shows that the differential geometric estimation has lower complexity when the system is fully observable.
Appendices
A. Time Derivative of State Vectors in Polar and Cartersian Coordinate
For simplicity, in this section, we assume , and . (1)First order time derivatives of polar coordinates is (2)The Cartesian representation in polar coordinates is (3)First order time derivatives of cartesian coordinates is
B. Conversion of Measurement Noise Covariance from Polar to Cartesian Coordinate
Reference [28] shows the expression cartesian coordinate measurement noise in terms of polar coordinate vector. Given that is the relative range, and are the relative elevation and azimuth angles at the current time step, respectively. The corresponding measurement noises in standard deviation are , and , respectively. Then, the measurement noise covariance, , in cartesian coordinate is given as where the “diag” represents the diagonal matrix, and is given as