ZHAI Guang,WANG Yanxin,and ZHAO Qi

1.School of Aerospace Engineering,Beijing Institute of Technology,Beijing 100081,China;2.Systems Engineering Research Institute,Beijing 100094,China

Abstract: This paper proposes partially norm-preserving filtering for a class of spacecraft in the presence of time-varying, but constant-magnitude maneuver. The augmented state Kalman filter(ASKF)is commonly used to track the maneuvering spacecraft with unknown constant propulsion;however,if the maneuver varies via time,the estimation performance will be degraded.To promote the tracking performance of the ASKF in case of time-invariant,constant-magnitude disturbance, the partially norm-preserving ASKF is developed by applying the norm constraint on the unknown maneuver. The proposed estimator, which is decomposed into two partial estimators and iteratively propagated in turns,projects the unconstrained maneuver estimation onto the Euclidian surface spanned by the norm constraint.The illustrative numerical example is provided to show the efficiency of the proposed method.

Keywords: maneuvering spacecraft, tracking, augmented state,norm-preserving,Kalman filter.

1.Introduction

Spacecraft tracking is critically important for satellite operators to maintain a clear picture of all space objects.The performance of spacecraft tracking significantly depends on an estimation of the dynamic state.The work of Kalman[1]and Luenberger[2],the techniques of state estimation are extensively studied.The Kalman filter generates an optimal solution when the motion of spacecraft is precisely modeled.However,in many practical cases,the spacecraft performs unknown maneuvers, if the maneuver sequence is not properly incorporated into the estimation, the filter evidently suffers from the performance degradation and sometimes even diverges[3].

In general, the unknown maneuvers can be considered as disturbance acted on spacecraft,and the significant challenge during the estimation is how to cope with the unknown disturbance.Over the past few decades, robust estimators for a linear system in the presence of unknown input have been developed,with a common approach being the unknown disturbance-decoupled state estimation[4].Within this approach,the disturbance vector generally satisfies the time-invariant condition [5], thus an optimal result could be obtained by minimizing the error covariance.To extend the disturbance-decoupled state estimator to the case in which the disturbance simultaneously affects the system’s state and output,Darouach&Zasadzinski[6]developed the optimal estimator filter(OEF)based on the parameterizing technique.Gillijns and De Moor[7]also proposed a three-step filter(TSF)for the joint disturbance and state estimate based on the direct feedthrough.The successful implementation of the TSF significantly depends on a full-rank feedthrough matrix of input to output.To relax the full-rank constraint on the feedthrough matrix, Cheng [8] proposed an unbiased minimum-variance estimator for the case of a feedthrough matrix of arbitrary rank. Other techniques in this area also include the minimax-optimal method[9]and the time-delayed optimal state estimator[10].

An earlier disturbance-decoupled example implemented for target tracking is the augmented state Kalman filter (ASKF), in which the constant unknown maneuver is treated as a part of the state and estimated synchronously[11,12].However,due to the increased filter dimension,the intensive computational load makes the ASKF very limited for implementation,and also,numerical problems may arise in the case of ill-conditioned systems[13].To reduce the computational load, Friedland [14] proposed a twostage Kalman filter by decomposing the ASKF into two parallel reduced order estimators. If the maneuver is constant, the two reduced order estimators propagate in turns and generate the optimal estimate of both the state and the maneuver [15]. Due to its high computational efficiency,this method is widely used to track a maneuvering target with constant propulsion[16]and achieves a better performance. To cope with the pulse-like unknown maneuvers,different multiple adaptive estimation methods have also been developed based on the ASKF.For example,the interacting multiple-model(IMM)method[17,18],which treats the unknown maneuver as part of the system state,switches the estimation among different filters based on maneuver detection [19]. When a maneuver is declared by the detector,the normal filter is augmented and the covariance is inflated to enable a re-convergenceon the maneuvering trajectories[20].Similar methods also include variable state dimension estimator [21], enhanced variable structure dimension(VSD)estimator[22]and unscented VSD estimator[23].

The ASKF is proved efficient only when the spacecraft is subjected to the constant unknown maneuver.However,the estimate becomes biased if the maneuver is timevarying,the faster the time-varying rate,the larger the estimation bias. Obviously, this performance degradation is mainly caused by the model mismatching on the unknown maneuvers. In general, to promote the estimation performance in the case of an unmatched model,the constrained Kalman filter (CKF), in which the state satisfies linear or nonlinear constraints,is commonly used for the estimation.By applying the constraints on the state, the estimate bias caused by model mismatching can be reduced.This can be seen,for instance,in the re-parameterization method[24],perfect measurement method [25] and projection method[24].Norm-preserving is a typical nonlinear constraint imposed on the entire or part of a system’s state. Given that the attitude quaternion satisfies the unit norm constraint,the earlier objective, to seek the norm-preserving solution of the Kalman filter, is for attitude estimation [26–28]. Zanet et al. [29] first developed the norm CKF for a time-discrete system by using the Lagrangian multiplier technique. In recent years, the norm-constrained filer is extended to the continuous-time system by considering a weighted norm [30]. The advantages of the work [29,30]have been numerically demonstrated by simulations of attitude determination.

Note that most of current earth-oriented spacecraft is propelled by the constant thruster,when describing the motion of spacecraft within the earth-centered-inertial frame,the acceleration generated by the thruster is time-varying but of constant magnitude. With a purpose of promoting the spacecraft tracking accuracies in presence of timevarying but constant magnitude maneuvers, a partially norm-preserving ASKF is proposed by applying the norm constraint on the unknown maneuvers.The novelty of this work is that we derive a partially norm-preserving ASKF with a two-stage formulation,and the ASKF is conditionally extended to a class of systems with time-varying disturbance,which is not previously considered in the literature.Additionally,we strictly prove the second order condition for optimality,based on evaluating the Hessian matrix.

2.Preliminary

Consider an earth-oriented spacecraft equipped with constant thrusters is moving in its orbit, within the earthcentered inertial frame the spacecraft motion can be represented with the discrete-time state-space equation as follows:

where xk∈ Rndenotes the system state, dk∈ Rmrepresents the maneuver vector of the well-known constant magnitude. Due to the attitude motion of spacecraft with respect to the earth centered inertial frame, the maneuver vector is always time-varying and can be described as follows:

where δk∈ Rm×mis a completely unknown time-varying matrix,because the attitude of the spacecraft moves slowly,we have

but the maneuver vector satisfies the norm constraint

where ρ > 0,Rk∈ Rp×pand yk∈ Rprepresents the measurement vector.Both the process noise wk∈Rmand the measurement uncertainties vk∈Rpare supposed to be independent,zero-mean,white random noise with wellknown covariance,denoted as Qk∈ Rm×m,and Qk>0 and Rk> 0. The time-depended matrices Ak, Gk, Γkand Ckare of appropriated dimensions that match the vectors.It also assumes that

and (Ak,Ck) is observable, while both (Ak,Gk) and(Ak,Γk) are controllable, the initial state x0is always independent of the process noise and measurement uncertainties, an unbiased estimation of initial state x0is uniquely determined by the initial covariance of

3.Unconstrained augmented Kalman filter

3.1 Augmented system

The ASKF takes the unknown maneuver as a part of state,therefore,the augmented state is formulated as follows:

where Xk∈Rn+mdenotes the augmented state vector,thus(1)and(2)are reformulated as follows:

where the over bar denotes an augmented symbol,and

3.2 Unconstrained ASKF(UASKF)

Suppose the initial maneuver is a Gaussian random variable and the unconstrained estimationis determined by

and the cross-covariance

Based on the augmented system (7), if we suppose the maneuver is invariant and take δk= 0, then the UASKF can be constructed as follows:

Note that if δk=0,the unconstrained estimators(15)–(21)work stablely and generate an optimal result.Since the dimension of the ASKF is increased to n+m, when m is compared to n, the computational load of the ASKF is also significantly increased.With a purpose of reducing the computational intensity of the estimation,the ASKF can be equivalently decomposed into two parallel reduced-order estimators. The two parallel reduced order estimators are the n-dimensional state estimator and the m-dimensional maneuver estimator,they propagate in turns,but subject to coupling equations.

Note that since δk= 0 for(1)and(2),we will demonstrate that the unconstrained estimators (15)–(21) produce a biased state estimate for(1)–(2).Define the unconstrained priori and posteriori estimate errors as follows:

Substituting(8)and(15)into(22)leads to

Similarly, substituting(19)and (20)into (23), the posteriori estimate error can be represented as follows:

Furthermore,substituting(24)into(25)yields

Decompose the posteriori estimate error and the augmented gain matirx as follows:

With(27),(26)can be reformulated as follows:

Taking the expectation of(28)and(29)yields

When time goes to infinite,i.e.,k → ∞,we have

Obviously, the state estimation is biased, which is mainly caused by the variation on the maneuver.

4.Partially norm constrained ASKF

4.1 Augmented performance index

In this section,in order to reduce the negative effects from the approximation in(4),major effort is dedicated to reconstructing (15)–(21) by incorporating the norm constraint on the maneuver.Define the augmented posteriori covariance matrix as follows:

where the partial covariance is defined as follows:

Additionally,consider the normalized maneuver estimation,satisfying the constraint in(5),

Discompose(19)as

we have

Substituting(39)into(27)leads to

Thus with the consideration of the norm constraint,the optimal estimation can be obtained by minimizing the augmented objective function:

Examine the right side of(41),it follows:

where tr{·} denotes the trace of a matrix. For notational convenience, define the partial performance index as follows:

where the scalar λk+1denotes the Lagrangian multiplier.Obviously,sinceandonly depend on the term ofandrespectively, therefore, the minimization on Jk+1can be performed separately onandthus we have

4.2 Partially constrained estimator for system state

(i)Prediction

(ii)Correction

where

4.3 Partially constrained estimator for maneuver

Next,we focus our attention to construct the partial estimator for maneuver by minimizing the performance indexwhich can be treated as an extremum problem with a constraint. Based on (15)–(21),can be computed as follows:

Consider the fact that

The minimization of trace is equivalent to one of the covariances, then substitute (55) into (44), take the partial derivative of(44)with respect toand λk+1and set the result to zero, then the first-order condition for the minimization can be represented as follows:

Expanding(56)and(57)leads to

and also,according to the principle of the matrix operation,we have

Substituting(61)into(60)leads to

Furthermore, substitute (62) into (59), a second order equation of λk+1is obtained as

where

By using Vieta theorem with(63),the optimal Lagrange multiplier can be obtained as follows:

where

Equation (66) indicates that there are two possible solutions for the optimal Lagrange multiplier, to determine the sign for the minimization, the second order condition should be examined. In general, the second order condition can be obtained by taking the second order derivative on the performance index,however,it will lead to the m×n dimensional matrix equation and significantly increases the computational complexities. Without the loss of generality, we consider the case that the gainis a column vector, where the measurement is reduced to a scalar,but the obtained results are still identical for the vector measurement.

Generally, the second order condition of an extremal problem can be founded by evaluating its Hessian matrix.The performance index achieves the maximum when the Hessian matrix is negative definite, while the minimum when the Hessian matrix is positive definite.In the following, we will show the performance indexis minimized when choosing the sign as plus in (66), while is maximized when choosing minus. When the gainis a column vector,the Hessian matrix of the performance indexcan be obtained by differentiating(56)again

Equation(68)shows that the Hessian matrix is diagonal,with(66),we have

With(67),(69)can be rewritten as follows:

Return to(52),it is not difficult to get

where Rk+1is the covariance of the scalar measurement noise. Furthermore, with (70) and (71), it clearly reveals that when the plus sign is chosen, the Hessian matrix is positive definite,and the minimum performance index occurs. Conversely, if the minus sign is chosen, the performance index is maximized.

In order to prove that the matrixis positive definite when choosing the plus sign in (66),rewrite the termas follows:

Multiply ηk+1on both sides and we have

Naturally, (75) indicates thatis positive definite, as a result of which, det(Θ +

(ii)Constrained Kalman gain and posterior covariance

The optimal Lagrange multiplier for the minimum performance index can be represented as follows:

Then substitute(76)into(62)and lead to

According to the Kalman filter algorithm, the constrained maneuver estimation can be obtained by correcting the one-step prediction as follows:

Furthermore,

where

Substitute(82)into(51)and lead to

Because the measurement noise is always uncorrelated with theandthus(83)can be reformulated by canceling the term related toand

Finally,we have

Equation (85) indicates that the constrained posterior covariance can be updated according to (54), but replacing the unconstrained gain as the constrained ones.

4.4 Performance analysis

Similarly, the constrained co-variance is updated as follows:

With(80),we have

and also

According to Cauchy-Schwarz inequality, we have

The result indicates that ifthe constrained estimator produces a result more accurate than the unconstrained estimator.

Rewrite(89)as follows:

Obviously,κkis also a representation of the estimation error,the smaller κk,the larger the error.Thus by selecting appropriate a threshold scale factor κT,it can be confirmed that only if κk> 1 orthe constrained estimator works with a better performance than the unconstrained one,otherwise the unconstrained estimator maintains an advantage over the constrained one.Moreover,it is easily to understand that the constrained estimate bias of the system state will be reduced when the maneuver is estimated with a better performance,as a result of which,the estimation accuracy of the system state is consequently improved.

4.5 Summary

The constrained estimation can be performed in the centralized way,however,the increased dimension of the centralized estimator will increase the computational load.With a purpose of reducing the computational complexity,it is recommended to propagate by the parallel way,in which,the centrialized estimator is seperated into a set of partial estimators of lower dimensions,the summary of the constrained estimator is shown in Table 1.

Table 1 Summary of the constrained estimator

5.Numerical results

As an illustrative example, the orbit determination of a non-cooperative spacecraft by using the proposed estimator will be simulated. Moreover, with a purpose of comparison,the unconstrained estimation results are also presented. The spacecraft is supposed to be three-axis stabilized and equipped with a constant magnitude apogee engine.During orbit maneuvering,the apogee engine always generates thrust acceleration in the along-track direction.For simplicity,we suppose the spacecraft initially works in a circular orbit, the orbit angular velocity, which is completely unknown to the ground station, is denoted as Ω,then with the earth centered inertial frame, the thrust acceleration can be described as follows:

where T is the well-known constant thrust force, ?0is the initial phase, while M is the mass of the spacecraft.Obviously,the thrust acceleration unknowingly varies via time, but always with constant magnitude.The following assumptions also hold: (i) The true orbit dynamics are only perturbed process noise.(ii)The attitude of spacecraft is properly controlled.(iii)The spacecraft is continuously tracked by ground stations.

The estimation performance,which is adopted to evaluate the relationship between constrained estimation accuracy and unconstrained ones quantitatively, is defined as follows:

5.1 Dynamics and observations

5.1.1 Orbit dynamics

Generally, the orbit motion of the spacecraft is described by using the two-body model.When the spacecraft undergoes an unknown maneuver,the motion can be represented with the linear perturbing acceleration within its dynamics

where μedenotes the earth g ravitational coefficient. r =are the orbital position and velocity vectors. The vector of w denotes the normally distributed process noises acting on the spacecraft.All of the vectors in(97)are represented with respect to the ECI frame. Define the state vector x as the combination of the position and velocity, then reformulate (97) by using space-state equations and yield

where

where D = F = [0 I]Tare the input and noise matrix,respectively.By using the Taylor series, (99) can be linearized in the vicinity of the current state as follows:

with(99)and(100),a discretization can be performed on the orbit dynamics,finally it leads to

where Φk+1,kis the state transfer matrix

where Δt is the time step for update.

5.1.2 Measurement

Suppose the spacecraft is tracked by ground-based phasedarray radar in a time-discrete fashion.Phased-array radar,which is fixed on the earth’s surface and follows the spinning of the earth, measures the spacecraft motion within radar-fixed range and angle coordinates. The nonlinear measurement model is represented as follows:

where

and

Generally,the vector

could be obtained by rotating the vector r from the ECI frame into the earth-centered earth-fixed(ECEF),and then into the SEZ frame as follows:

where TECI→ECEFis the rotation matrix from the ECI frame to the ECEF frame, while TECEF→SEZis the one from the ECEF frame to the SEZ frame,whilerepresents the location of the radar within the ECEF frame.With a similar rotation,the velocity vector can also be expressed as

We then have

where

Substitute (109) into (106), then the measurement can be represented in terms ofx(k).Furthermore,linearize the measurement in the vicinity of state prediction and yield

whereHkis the measurement Jacobian matrix. It can be calculated as follows:

5.2 Numerical results

Consider a spacecraft in an equatorial orbit initially. The orbit eccentricity and altitude are 0 and 600 km, respectively.The mass of the spacecraft is 1 000 kg,and the constant thruster with a magnitude of 20 N is fixed on the spacecraft. During the entire simulation, the thruster always works to generate the sinusoidal maneuvering acceleration within thex-yplane.The initial conditions for the estimation are listed as follows:the initial orbit statex0=[6 978 000 m,0 m,0 m,0 m/s,7 557.9 m/s,0 m/s], the initial maneuverd0= [0 m/s2,0 m/s2,0 m/s2], the initial covariance for the state error is ofwhilefor the unknown maneuver acceleration. The covariance of the dynamic noise is set toQ=10?12×I3×3,whileR=diag([40 000,400,7.62×10?9,7.62×10?9]) for measurement uncertainties. The estimator works at a frequency of 1 Hz,and the simulation totally takes one period of the initial circular orbit.

The unconstrained maneuver estimations along both directions are represented in Fig.1,where the dash lines denote the true maneuvers, while the solid lines are the estimated results. From the left plots, one can see that the unconstrained estimator suffers from dramatic oscillations due to the initial errors,after the initial oscillations,the estimated maneuvers follow the true ones but with a remarkable time delay, which can be explained by the fact that the unconstrained estimator always suppose the maneuvering acceleration to be constant,therefore the unconstrained estimator cannot reflect the variation of the maneuvers in time. The right plots show the maneuver estimation error on the maneuvers, the results reveal that the errors along both directions liberate throughout the simulations,and the maximum errors are all greater than 0.005 m/s2.

Fig.1 Unconstrained estimation results of maneuver

The most interesting results are shown in Fig.2. Compared with the results in Fig.1,it clearly reveals that once the norm constraint is applied to the maneuver estimation,the unconstrained estimations are normalized to preserve the magnitude of the maneuver, and the constrained estimator becomes more sensitive to the variation of the maneuver,as a result of which,the drastic oscillations at the beginning are significantly suppressed,moreover,the time delays of the estimation, as well as the estimation errors,are both reduced.

Fig.2 Norm constrained estimation results of maneuver

Fig.3 shows the magnitudes of the estimated maneuvering accelerations.It clearly shows in the case of the unconstrained estimation, the estimated maneuver dramatically oscillates within the first 100 s, but when the estimation reaches a relatively stable stage, the estimated magnitude is smaller than 0.02 for most of the time.However,for the constrained estimation, because we select the parameterκT= 0.975, the magnitude of the maneuvering acceleration is projected into a range of [0.019 5, 0.02], which makes the estimated maneuvering acceleration closer to the true ones.

Fig.3 Norm of estimated maneuver

The estimation errors on the orbit position and velocity are represented in Fig. 4 and Fig. 5, where the solid lines represent the estimation errors, while the dash lines represent theσboundary. Due to the time delays of the maneuvering acceleration estimation, the estimation errors of both constrained and unconstrained cases suffer from the initial librations.After reaching a relatively stable stage,the unconstrained estimator approximately achieves an accuracy of 30 m (σ) along both directions, while the constrained estimator achieves a better accuracy of about 20 m(σ),and also,the accuracies on the velocity are promoted after applying the norm constraint on the maneuver,as shown in the right plots within each figure.The results clearly reveal that the estimation accuracy is improved after incorporating the constant norm constraint on the unknown maneuvering acceleration.

Fig.4 Position errors of unconstrained estimation

Fig.5 Position errors of constrained estimation

We also increase the maneuvering acceleration from 0.02 m/s2to 0.05 m/s2, the estimation errors on the orbit position and velocity are represented in Fig.6 and Fig.7.It is obvious that we get the similar results but with the increased estimation error.After reaching a relatively stable stage,the unconstrained estimator approximately achieves an accuracy of 55 m (σ) along both directions, while the constrained estimator achieves a better accuracy of about 40 m(σ).The results clearly reveal that the estimation accuracy will be degraded after increasing the magnitude of the unknown maneuvering acceleration.

Fig. 6 Position errors of unconstrained estimation after increasing maneuvering acceleration

Fig. 8 shows the estimation performance of the constrained estimator relative to the unconstrained ones.One can see that at the beginning of the simulation,the estima-tion performance along both directions dramatically oscillates, which are mainly caused by the initial convergence of the both estimators, but after about 10 s, the performances are always greater than 1,and then greater than 1.2 after about 1 500 s,finally it achieves the maximum values of 1.5 and 1.4 alongxandydirections,respectively.Based on the results in Fig. 8, it can be confirmed that the constrained estimator obtains an enhanced performance than the unconstrained ones.

Fig.7 Position errors of constrained estimation after increasing maneuvering acceleration

Fig.8 Estimation performance of constrained estimator

6.Conclusions

The primary contribution of this work is the development of the partially norm-preserving augmented estimator.The proposed estimator is implemented to track spacecraft subjected to the time-varying maneuver but of constant magnitude.Based on the framework of the traditional augmented Kalman filter,the norm constraint is incorporated into the estimation, the estimation errors are minimized by using the Lagrangian multiplier technique, and with the constrained Kalman gain matrix, the estimation is projected to the Euclidian surface spanned by the norm constraint.Finally, to highlight the advantages of the proposed estimator, both the unconstrained and constrained estimators are simulated to track a maneuvering spacecraft.Numerical results indicate that the proposed estimator achieves an enhanced performance.