Non-cooperative maneuvering spacecraft tracking via a variable structure estimator

Accepted Manuscript Non-cooperative maneuvering spacecraft tracking via a variable structure estimator

Zhai Guang, Bi Xingzi, Zhao Hanyu, Liang Bin

PII: DOI: Reference:

S1270-9638(18)30045-2 https://doi.org/10.1016/j.ast.2018.05.052 AESCTE 4605

To appear in:

Aerospace Science and Technology

Received date: Revised date: Accepted date:

6 January 2018 26 March 2018 29 May 2018

Please cite this article in press as: Z. Guang et al., Non-cooperative maneuvering spacecraft tracking via a variable structure estimator, Aerosp. Sci. Technol. (2018), https://doi.org/10.1016/j.ast.2018.05.052

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Non-cooperative Maneuvering Spacecraft Tracking via a Variable Structure Estimator Zhai Guang1, Bi Xingzi2, Zhao Hanyu3 Aerospace Engineering School, Beijing Institute of Technology. 100081, Beijing, P.R China Liang Bin4 National Laboratory for Information Science and Technology, Tinghua University. 100084, Beijing, P.R China. Abstract—A spacecraft performing unknown maneuvers significantly increases the probability of a collision between spacecraft with neighboring orbits. Traditional orbit determination filters cannot robustly cope with unknown maneuvers. To track non-cooperative, unknowingly maneuvering spacecraft, a novel variable structure estimator is proposed based on the input compensation to a normal filter. We develop an observer to estimate the unknown maneuver and then analyze the estimated maneuver performance based on error propagation. Once the maneuver detector declares the occurrence of an unknown maneuver, the estimated maneuver is fed to an extended Kalman filter as compensation, which enables the estimator to work adaptively. Finally, a series of numerical simulations are carried out to demonstrate the effectiveness of the proposed variable structure estimator, and the estimation performances in the cases of constant and time-varying maneuvers are analyzed based on the simulation results. Index Terms—Orbit determination; Non-cooperative target; Unknown maneuver; Observer; Extended Kalman filter; Variable structure estimator.

I.Introduction The techniques of tracking, predicting and estimating the orbit states of non-cooperative spacecraft are critically important for Space Situational Awareness (SSA). Due to the rapid growth of the spacecraft population, possible collisions between spacecraft in neighboring orbits are attracting increasing attention [1]. To guarantee the safety of operational space assets, satellite operators must have a clear picture of all space objects. In general, Kalman filter has been wildly used for orbit determination, it is not difficult to track a cooperative spacecraft given a well-known, planned maneuver sequence [2-4]. However, if the maneuver sequence is unknown to the tracking system, the estimation performance of a traditional filter will be significantly degraded, and sometimes the estimator even diverges [5]. In general, the divergence problem caused by unknown maneuvers could be solved by inflating the covariance [6,7]. The one-time covariance inflation method has been proved effective for coping with the impulsive maneuvers performed by noncooperative spacecraft. To prevent filter estimation divergence, this algorithm inflates the covariance to increase the weight of fresh observation during the state update stage, which allows the estimation to converge on the post-maneuver trajectory. However, this approach becomes limited if the unknown maneuvers are continuous. In contrast, the fading memory filter continuously increases the weight of the current measurement through covariance inflation [8,9], but it also increases the 1Associate Professor, Beijing Institute of Technology, Aerospace Engineering School, Post Office Box 25. No. 5, Zhongguancun South Street, Beijing. [email protected] 2 Graduated Student, Beijing Institute of Technology, Aerospace Engineering School, Post Office Box 25. No. 5 Zhongguancun South Street, Beijing. [email protected]. 3Graduated Student, Beijing Institute of Technology, Aerospace Engineering School, Post Office Box 25. No. 5 Zhongguancun South Street, Beijing. [email protected]. 4 Professor, Tsinghua University, National Laboratory for Information Science and Technology. Haidian District, Beijing. [email protected].

weight of measurement errors and sequentially degrades the tracking accuracy. Although the performance of the fading memory filter is not satisfactory for highly precise tracking, it provides a feasible way to deal with the maneuvering target tracking (MTT) problem. Many efforts have been made to address the problem of MTT in the area of orbit determination, and various adaptive methods have been proposed to prevent the divergence caused by unknown maneuvers. Previous work in this area can be broadly divided into two categories: data association and real-time tracking [10]. Data association studies focus on matching a track of observations with a spacecraft’s potential maneuvers, such as multiple hypothesis tracking filters [11] and joint probabilistic data association filters [12]; these methods are characterized by their post-processing and uncorrelated target resolution. One of the earliest real-time tracking techniques is the equivalent noise approach, within which the maneuver is modeled as white or colored equivalent orbital noise, and then the problem of MTT is converted to state estimation in presence of non-stationary process noise with unknown statistics [13,14]. Singer [15] originally proposed an augmented Kalman filter (AKF) by treating the unknown maneuver as part of a state within which the unknown maneuvers are modeled as a zeromean first-order stationary Markov process driven by white random noise. When the maneuvers occur, the system dynamics are augmented by incorporating the maneuver thrust into the system state, and then the covariance is inflated to initialize a new convergence on the maneuvering trajectories. AKE tracks the maneuvering spacecraft well if the spacecraft performs constant maneuvers, but when the unknown maneuver is time-varying, the accuracy of the orbit determination will decrease significantly. In contrast to the single-model estimator, the switching model approach (SMA) takes advantage of both standard KF and AKF, and thus it provides better tracking performance of both non-maneuvering and maneuvering spacecraft. The variable state dimension (VSD) algorithm [16], the enhanced VSD algorithm [17] and two-stage filtering [18] are in the scope of SMA. When using SMA, detectors are commonly required to determine when to switch the estimation between the standard Kalman filter and AKF [19]. As the performance of SMA is deeply influenced by maneuver detection, several detectors have been developed for maneuver detection, such as the simple chi-square test (using either residuals [20] or input estimates [21]), generalized likelihood ratio (GLR) test [22], marginalized likelihood ratio (MLR) test [23] and sequential statistical test [24]. Holzinger et al. [25] proposed an optimal control performance metric to detect and characterize orbital maneuvers. By checking the inconsistency between filtered and smoother state estimates, Kelecy [26] developed a type of finite maneuver detector based on batch least squares and the EKF; this detector has been successfully used for maneuvering target tracking. Although many efforts have been made to improve the sensitivities of the maneuver detectors, they remain limited in their ability to respond immediately to a maneuver start or conclusion [27, 28]. The input estimation (IE) approach [29-32] also provides a way to cope with unknown maneuvering spacecraft. The IE approach detects the existence of spacecraft maneuvers and directly estimates the magnitude of the unknown maneuvers. The estimated maneuver is fed into the normal filter as compensation. However, the IE approach always assumes the unknown maneuver as constant; if the actual maneuver varies via time, the accuracy of orbit determination will become worse. Gholson and Moose [33] proposed a new approach treating the unknown maneuvers as a semi-Markov process. Within this method, the possible acceleration inputs are selected according to a priori probabilities. Similarly, Kumar [34] developed an adaptive mean statistical algorithm by modeling the unknown maneuver as a modified Rayleigh-Markov process. The interacting multiple model (IMM) algorithm has recently been widely investigated for tracking unknown maneuvering spacecraft [35-38]. IMM works well when the sub-modes accurately represent the actual acceleration; it is also relatively robust to small modeling errors. However, because the computational load of IMM increases linearly with the number of maneuver modes, the hardware with high performance is a basic requirement for the IMM implementations. To enhance the robustness of the estimator, the aforementioned algorithms can also be applied integrally. For instance, Yong developed a combined estimator by integrating the VSD filter with the input estimation algorithm [39] and achieved better performance

than the single approach. Other similar estimators include the VSD-IMM algorithm, the IMM-FT (filtering-through) algorithm and the thrust-Fourier-coefficient event representation with the filtering-through algorithm [40]. Motivated by the problem of tracking spacecraft engaged in unknown, time-varying maneuvers, we develop a novel variable structure estimator based on the input compensation to the EKF in this work. Based on observations from groundbased radar, a maneuver observer is first constructed to estimate the unknown maneuver; the error propagation during maneuver estimation is also analyzed. The maneuver estimation is then fed to the EKF once the maneuver has been declared by the detector. Rather than the constant maneuver assumption within the aforementioned IE approaches, the variable structure estimator developed in this work can adaptively cope with unknown, time-varying maneuvers by regulating the observer matrix, and the negative effect from the maneuver-varying rate on estimation performance can also be controlled. Finally, a series of simulations is carried out to demonstrate the performance of the proposed estimators, and the advantages and disadvantages are discussed on the basis of the simulation results.

II.Orbit Determination Using Extended Kalman Filter A. Two-Body Dynamics Generally, the orbit motion of spacecraft is described 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 dr =v dt

dv μ =− 2 dt r

(1)

2 2 º r 2 z °½ °­ r 3J 2 Re ª 5z − 1» − Zˆ ¾ + u + w ® − 2 « 2 2 r r r ¬ ¼ r r ¿° ¯°

(2)

where μ , Re and J 2 denote the Earth gravitational coefficient, equatorial radius and second-order Earth oblateness

 y,z   ]T are the orbital position and velocity vectors, r = x 2 + y 2 + z 2 is the coefficient, respectively. r = [x, y,z ]T and v = [x, orbital radius, u ∈ R3×1 represents the unknown maneuver acceleration and Zˆ represents the unit vector of the Earth-centered inertial (ECI) frame. For simplicity, only the J2 perturbation term is taken into account; the other disturbance terms are neglected. The vector of w denotes the normally distributed process noises acting on the spacecraft. All of the vectors in Eq. (1) are represented with respect to the ECI frame. Define the state vector x as the combination of the position and velocity, then reformulate Eqs. (1) and (2) using space-state equations. This yields dx =f (x ) + Du + Fw dt

(3)

where f ( x ) = [ x , y , z, − μ

x x y y z z + μ 5 , − μ 3 + μ 5 , − μ 3 + μ 5 ]T 3 r r r r r r

(4)

and x = [r T , vT ] T is the six-dimensional state vector. D ∈ R 6 × 3 and F ∈ R 6 × 3 are the input and noise matrix, respectively. By using the Taylor series, Eq. (4) can be linearized in the vicinity of the current state xk : f ( xk +1 ) ≈ f ( xk ) +

∂f ( x ) ∂x

( xk +1 − xk ) .

(5)

x = xk

With Eqs. (4) and (5), the discretization can be performed, and the dynamics can be represented by the state transfer function as xk +1 = ĭk +1,k xk + Bk uk + Gk wk 6×6 where ĭk +1,k ∈ R is the state transfer matrix

(6)

ĭk +1,k = I 6×6 + Δt

∂f ( x ) ∂x x = xk +1

(7)

and satisfies ĭ(tk +1 , tk )ĭ(tk , tk +1 ) = I6×6

(8)

2.2 Measurement from Ground-based Radar In this work, it is assumed that the spacecraft is continuously tracked by ground-based radar. Phased-array radar is widely used to measure the three-dimensional motion of spacecraft within radar-fixed range and angle coordinates. Because groundbased radars are fixed on the Earth’s surface and move by following the Earth’s rotation, measurements within the radarfixed reference frame should be transformed to ECI coordinates based on radar navigation and attitude reference data, which are commonly provided by the strap-down inertial navigation system (INS) mounted on the radar face.

ȡ

Oz Ozradar

r

Ozradar

Oxradar Oy

ȡ

ϕ Oyradar

γ

O

Oxradar

Ox Oyradar

Fig. 1 Geometric Illustration of Ground-based Radar Measurements Generally, INS is erected local vertical and pointed north after the gyro-compassing alignment procedure. After the tracking begins, the radar measures the spacecraft motion within the South-East-Zenith (SEZ) frame, and finally generates the output as follows: zk = h( xSEZ (k )) + ȣk

(9)

and 2 2 2 ª xSEZ º + ySEZ + zSEZ « » 2 2 2 « x SEZ + y SEZ + zSEZ » « » ySEZ » h( xSEZ ( k )) = « a tan( ) « » ρ « » « » zSEZ a sin( ) « » ρ ¬ ¼

(10)

where state vector x ESZ = [ ȡT , ȡ T ]T is defined within the SEZ frame, ȡ = [ xSEZ , y SEZ , zSEZ ]T denotes the spacecraft position, zk = [ ρ k , ρ k ,ϕ k , γ k ]T is the measurement output that includes the range, range rate, azimuth, elevation , while ȣk indicates the

measurement uncertainties. The vector ȡ could be obtained by rotating the vector r from the ECI frame into the EarthCentered Earth-fixed (ECEF) frame as ECEF ρ ECEF = ȉ ECI → ECEF r − pRadar

(11)

ECEF represents the location of where ȉ ECI →ECEF denotes the rotation matrix from the ECI frame to the ECEF frame, while pRadar

the radar within the ECEF frame. Next, rotate the vector ȡ ECEF into the SEZ frame as ECEF ρ = ȉ ECEF → SEZ (ȉ ECI → ECEF r − pRadar )

(12)

where ȉ ECEF →SEZ denotes the rotation matrix from the ECEF frame into the ESZ frame. With a similar rotation, the velocity vector of ȡ can also be expressed in term of r and r as ρ = ȉ ECEF → SEZ ȉ ECI → ECEF r + ȉ ECEF → SEZ ȉ ECI → ECEF r

(13)

ECEF x ECEF = Θx + ΛpRadar

(14)

With Eqs.(12) and (13), then we have

where ȉ ªȉ Θ = « ECEF → SEZ  ECI → ECEF ȉ ¬ ECEF → SEZ ȉ ECI → ECEF

0 º ȉ ECEF → SEZ ȉ ECI → ECEF »¼

ªȉ º Λ = « ECEF → SEZ » 0 ¬ ¼

(15)

(16)

Finally, the measurement in Eq. (14) can be linearized at the vicinity of the current state prediction with the Taylor series as (17)

zk +1 = h( x k ) + H k ( xk+1/ k − x k ) + ȣk +1

where H k is the measurement Jacobian matrix. It can be calculated as Hk =

∂h( xk ) ∂x k

= x = x k

∂h( x SEZ ( k )) ∂x SEZ ( k )

x = x k

∂x SEZ ∂x

(18) x = x k

B. Traditional Extended Kalman Filter The extended Kalman filter (EKF) has been wildly used for orbit determination based on ground observations. The EKF can be viewed as a method for optimal estimation in the minimum mean square error (MMSE) sense. If the system inputs are well known, EKF generates the near-optimal estimation in real time. In general, the propagation of the EKF can be divided into two stages. The first stage is one-step prediction based on the dynamics: xˆ k+1/k = ĭk+1,k x k + Buk

(19)

T Pk+1/k = ĭk+1,k Pk ĭk+1,k + GQG T

(20)

where x k and Pk are the state estimate and the posteriori covariance, respectively, while xˆk+1/k and Pk+1/k are the one-step predictions. Q is the covariance of the dynamic noise. If the dynamic noise is uncorrelated, zero-mean and Gaussian distributed, the covariance of w(k ) can be presented as E( wi w j ) = Qδij

(21)

The symbol δ ij is the Kronecker function, and Q is the diagonal matrix Q = diag ([ qx2 ; q 2y ; qz2 ;])

(22)

where q x , q y and q z are the noise standard deviations in three directions. The second stage is the one-step update, propagated as K k +1 = PkT+1,k H kT+1 ( H k +1 Pk +1,k H kT+1 + Rk +1 ) −1

(23)

Pk +1 = ( I − K k +1 H k +1 ) Pk +1,k

(24)

x k+1 = xˆk+1 + Kk +1Șk +1

(25)

where K k +1 is the Kalman gain matrix, x k+1 indicates the state estimation and Șk +1 denotes the measurement residual, which can be calculated as Șk +1 = zk +1 − zˆk +1

(26)

R denotes the covariance of measurement uncertainties that are also uncorrelated, zero-mean and Gaussian distributed: E ( ȣi ȣ j ) = Rδ ij

(27)

R = diag ([r12 , r22 , r32 , r42 ])

(28)

where ri 2 is the standard deviation for each measurement component.

III.Unknown Maneuver Observer As stated earlier, if the orbital maneuvers are well known, the EKF works effectively to generate a near-optimal estimation of orbit trajectory. However, if the spacecraft performs continuous maneuvers of unknown amplitude at an unknown time, the estimation diverges and the spacecraft is likely to be lost. When the EKF reaches a relatively steady state, both the Kalman gain and covariance matrix remain nearly constant; therefore, even if the measurement residual is increased, neither the Kalman gain nor the covariance matrix can automatically respond to the dynamic uncertainties caused by the unknown maneuvers. The aim of this work is to develop a novel Variable Structure Estimator (VSE) used for orbit determination in the presence of time-varying unknown maneuvers. Similar to the VSD estimator, the proposed VSE also varies its structure when the maneuver start/conclusion is detected, but it never assumes the unknown maneuver to be constant. The main idea of the VSE is to use a first-order observer to estimate the unknown maneuver based on observations. To respond to unknown maneuvers, the estimated maneuvers are fed into the standard KF; hence, if the observer is stable, the EKF could ensure re-convergence on the maneuvering trajectory. A. Basic Maneuver Observer In this section, the first-order observer used to estimate the unknown maneuver acceleration is proposed. To estimate the maneuver acceleration, the straightforward idea is to construct the basic observer as follows: uˆ k+1 = ( I − Ȍ k ) uˆ k + Ȍ k uk

(29)

where uˆ k denotes the estimation of the maneuver and Ȍ k is the observer gain matrix. Without loss of generality, the matrix Ȍ is chosen as

Ȍ k = diag[λ1 , λ2 , λ3 ]

(30)

λi − 1 < 1 i = 1,2,3

(31)

Naturally, if

is satisfied, then the basic observer in Eq. (29) is globally asymptotically stable. Furthermore, reformulate Eq. (29) and express the basic observer equivalently as Δuˆ k +1 = Ȍ k ( uk − uˆ k )

(32)

As Eq. (32) shows, if the matrix Ȍk is appropriately selected, it ensures that uˆ k decreases when uˆk > uk , while uˆ k increases when uˆk < uk . Define the error of estimation as e k +1 = uk +1 − uˆ k +1

(33)

Δuˆ k +1 = uˆ k +1 − uˆ k − uk +1 +uk +1 +uk − uk = uˆ k +1 − uk +1 − uˆ k +uk +uk +1 − uk

(34)

Then, rewrite the right side of Eq. (32), which yields

= − ek +1 + ek +Δuk +1

Substitute Eq. (34) into Eq. (32), which leads to ek+1 = ( I − Ȍ k )ek + Δuk+1

(35)

where Δuk+1 represents the time-varying rate of the unknown maneuver. Obviously, Eq. (35) is the error dynamics when using the basic observer to estimate the unknown maneuvers. Because the observer matrix Ȍk satisfies the conditions in Eq. (31), if Ȍk is properly selected, then the error of the steady state depends only on the term of Δuk+1 . Generally, it is

reasonable to assume that the rate of the maneuver is always bounded; that is, at any time tk, there exists the constant τ i ≥ 0 that always satisfies [ Δux (k) , Δu y (k) , Δuz (k) ] ≤ T[τ1 ,τ 2 ,τ 3 ]

(36)

where T represents the time step. With the assumption in Eq. (36), the estimate error in Eq. (35) can also be proved to be bounded. B. Order Reduction Next, according to Eq. (6), the maneuver acceleration can be represented as uk = BK+ ( xk +1 − ĭk +1,k xk )

(37)

where B + represents the left-invertible matrix of Bk and satisfies B + Bk = I . Substitute Eq. (37) into Eq. (32); then, the K

K

unknown maneuver acceleration can be observed as follows: Δuˆ k +1 = − Ȍ k ( uˆ k − B + ( xk +1 − ĭk +1,k xk ))

(38)

Naturally, if all of the terms related to x k +1 and x k are available, the basic observer in Eq. (38) can be used to estimate the unknown maneuver accelerations through the sequential time update. However, because the term of x k +1 is always unavailable at the current time of t = t k , the basic observer in Eq. (38) should be reconstructed by eliminating the terms associated with x k +1 . Thus, an auxiliary variable

ȟˆ k

is introduced as ȟˆ k = uˆ k − īxk

(39)

where ī is the 3 × 6 matrix to be determined. By substituting Eq. (39) into Eq. (38), the basic observer can be reformulated as Δȟˆ k +1 = − Ȍ k ( uˆk − BK+ ( xk +1 − ĭk +1, k xk )) − ī (xk +1 − xk ) = − Ȍ k ((ȟˆ k + īxk ) − BK+ ( xk +1 − ĭk +1, k xk )) − ī (xk +1 − xk )

(40)

= − Ȍ k (ȟˆ k + īxk ) + ( Ȍ k BK+ − ī ) xk +1 − ( Ȍ k BK+ ĭk +1, k − ī ) xk

Obviously, by selecting

ī = Ȍk B+ ,

the term of x k +1 can be eliminated as ȟˆ k +1 = ( I − Ȍ k ) ȟˆ k − ( Ȍ k ī + Ȍ k B K+ ĭ k +1, k − ī ) x k

(41)

Eqs. (39) and (41) provide a type of reduced-order observer that depends only on the current state. Because x k can be calculated by using the current observations, the unknown maneuver acceleration can be estimated based on the proposed observer. C. Noise Propagations Most importantly, according to the structure of the VSE, the estimation of the maneuvers is fed into the EKF as input compensation, and it is therefore necessary to evaluate the estimation performance. Here, we use the notation xk to represent the noise-corrupted state obtained from the observation. According to Eq. (17), xk can be computed as follows: h( xk-1 ) + H k ( xk − xk-1 ) = h( xk-1 ) + H k ( xk − xk-1 ) + ȣ k +1

(42)

Then, we have xk = xk − Hk+ ȣk

where H k+ represents the left-invertible matrix of Hk . Substituting Eqs. (39) and (6) into Eq. (33) yields

(43)

ek +1 = uk +1 − (ȟˆ k +1 + īxk +1 ) = uk +1 − ī (ĭk +1, k xk + Bk uk )) − ȟˆ k +1 = uk +1 − ī (ĭk +1, k xk + Bk uk )) − (( I − Ȍ k )ȟˆ k − ( Ȍ k ī + Ȍ k BK+ ĭk +1, k − ī ) xk ) = uk +1 − (ȟˆ k + īxk ) + Ȍ k (ȟˆ k + īxk ) − ( īĭk +1, k − Ȍ k BK+ ĭk +1, k ) xk − īBk uk = uk +1 − (ȟˆ k + īxk ) + Ȍ k (ȟˆ k + īxk ) − ( ī − Ȍ k BK+ )ĭk +1,k xk − īBk uk

Because we have

Ȍk B+ = ī

(44)

, ek +1 = uk +1 − (ȟˆ k + īxk ) + Ȍ k (ȟˆ k + īxk ) − īBk uk

(45)

Substituting Eq. (43) into Eq. (45) leads to ek +1 = uk +1 − (ȟˆ k + īxk ) + Ȍk (ȟˆ k + īxk ) − īBk uk + ( ī − Ȍk ) H k+ ȣ k = uk +1 − uˆk + Ȍ k uˆk − īBk uk + ( ī − Ȍ k ) H k+ ȣ k = uk +1 − uk + uk − uˆk + Ȍk uˆk − īBk uk + ( ī − Ȍk ) H k+ ȣk

(46)

+ k

= (uk − uˆk ) − Ȍ k (uk − uˆk ) + (uk +1 − uk ) + ( ī − Ȍ k ) H ȣ k ek

ek

Δuk +1

= ( I − Ȍk )ek + Δuk+1 + ( I − Ȍ k ) īH k+ ȣ k

As the time goes to infinity, k → ∞ , we have ek = ( I − Ȍ k )ek + Δuk + ( Ȍ k − I ) īH k+ ȣk

(47)

ek = Ȍk −1 (Δuk + (Ȍk − I ) īHk+ȣk )

(48)

Reformulating Eq. (47) leads to

As stated earlier, because Ȍk B+ = ī and Bk = Gk , we have k

ek = AΔuk − Dȣk

(49)

where A = diag (

1 1 1 , , )

(50)

λ1 λ2 λ3

D = diag (λ1 − 1, λ2 − 1, λ3 − 1) Bk+ H k+

(51)

Then, with Eq. (49), the covariance in the x direction can be evaluated as follows: E(ekxekx ) ≤

1

λ

2 1

4

2 T2 ⋅τ x2 + (λ1 − 1)2 (¦b1i ri )2 = qmx

(52)

i =1

2 represents the covariance of the estimated maneuver in x direction. where b1i denotes the element of matrix B + H k+ and qmx k

Similarly, the covariance in the y and z directions has the same expression as in Eq. (52). From Eq.(52), the estimation performance of maneuver observer is mainly effected by observer gain parameters λi , the rate of acceleration change τ x , observed period T and measurement noise ȣk . From Eq. (52), one can see that for the case of constant unknown maneuvers, that is, τ x = 0 , the error of the estimated maneuver depends only on the measurement uncertainties. Because the estimated maneuver is corrupted by the uncertainties, when it is fed to the EKF, the estimation error must be taken into account as additional dynamic noise. Therefore, the original dynamic covariance Q should be replaced by Q ′ as Q′ = Qm + Q

(53)

2 2 2 Q m = diag ( qmx ,qmy ,qmz )

(54)

where

IV.Variable Structure Estimator Based on Input Compensation A. Framework of the Variable Structure Estimator

Based on the proposed unknown maneuver observer, the Variable Structure Estimator (VSE) is proposed to adaptively cope with the unknown maneuver. Within the framework of the VSE, if the spacecraft moves freely, the state of the spacecraft is estimated routinely by using the EKF. However, if the unknown maneuver is detected, the structure of the estimator will be changed by incorporating the maneuver observer into the propagations, where the observer is used to estimate the unknown maneuver, and then the maneuver estimation is fed to the EKF as an input compensation. Additionally, if the end of the maneuver is declared, the estimator will switch back to the EKF by switching off the observer. The structure of the VSE is shown in Fig. 2. Using the VSE, the estimation convergence could be guaranteed even if the unknown maneuver occurs, and the unknown maneuver and spacecraft state could be estimated synchronously with high quality.

Fig. 2. Structure of VSE B. Maneuver Detector According to the structure of the VSE, it is significantly important to change the estimator structure exactly at the time of maneuver start or end; otherwise, the estimation error will be built up and the estimation performance will be degraded. Therefore, the maneuver detectors should be developed to declare the maneuver start and end as accurately as possible. In general, if the unknown maneuver occurs, the estimation covariance remains nearly constant when EKF is used to process observations, but the observational residuals in Eq. (26) will increase to reveal the maneuver. To construct the maneuver start detector, a straightforward method is to detect the maneuver occurrence depending on the sensitivity of the residual on the maneuvers; then, the detector can be represented as ϒ ks = ȘTk ( H k Pk HTk + R k ) −1 Șk

(55)

Obviously, ϒ is a scalar function in terms of the observational residual; it is updated during the estimation propagations. s k

If the maneuver has occurred, ϒ ks will grow with the increase of Șk +1 . Moreover, with additional consideration of the inherent noise of the residual, Eq. (54) can be applied to declare the maneuver start if ϒ ks exceeds a predefined threshold. If the maneuver start is declared, then the observer is initialized to estimate the maneuver. Similarly, another detector is required to detect the end of the maneuver. When the spacecraft is subjected to the maneuver, the unknown maneuver is estimated by the observer; then, the end detector can be developed by calculating the magnitude of the estimated maneuver as follows: ϒ ck = uTk ( Pkth ) −1 u k

(56)

Based on Eq. (56), it is easy to understand that once the magnitude of the maneuver estimation falls to predefined thresholds, the maneuver could be declared to the end; then, the observer will be terminated and the estimator will switch back to the EKF. C. Initialization of the Observer

When the observer is initialized to compensate the EKF, the estimator unavoidably suffers from transient oscillation and degrades the estimation performance. To suppress the undesirable transient oscillation, the observer should be initialized with a given initial condition that is much closer to the true state. Generally, the observer can be initialized in different ways. If the amplitude of the maneuver is roughly known as u , then

ȟˆ 0

can be set as ȟˆ 0 = u − īx0 , which allows the observer to

estimate the unknown maneuver with an smaller initial error. The possible oscillation caused by the inaccuracy of the maneuver estimation can thus be avoided. Additionally, if the amplitude of the maneuver is completely unknown, u can be set as zero and the maneuver can be estimated initially as − īx0 . The proposed initializing method is advantageous for smoothing the output of EKF compensated by the maneuver observer.

V.

Numerical Results

To verify the performance of the VSE developed in this work, the orbit determination problem of a spacecraft subjected to unknown maneuvers is simulated. The spacecraft moves in a circular orbit, and it is observed continuously by ground radars equivalently spaced on the Earth’s surface. Within the simulations, the observational and dynamic uncertainties are assumed to be zero-mean Gaussian distributed with constant covariance, and all of the orbital perturbations are ignored within the simulations except the J2 term. The transformation matrix associated with the measurement representations are treated as deterministic. The estimator and the measurement work at a frequency of 5Hz. To demonstrate the estimation performance in case of different types of unknown maneuvers, we apply the constant and time-varying unknown accelerations on the spacecraft. Some important parameters for the simulation are listed in Table 1. Table 1. Simulation Parameters parameters Orbital elements

a, i, e, Ω,ω, M0 Measurement uncertainties ( [δρ k , δρ k , δϕ k , δγ k ] )

Value 7200km,30deg,0.02,5deg,45deg,0deg

[300m,50m / s,0.01°,0.01°]

Dynamic noises (m/s2)

10−12

Initial position error(m)

[200, 200, 200]

Initial velocity error(m/s)

[50,50,50]

A. EKF Results In this section, to illustrate the estimation performance of the EKF, two different scenarios are simulated for the problem of orbit determination. The first scenario considers the case in which the spacecraft orbits the Earth freely; therefore, the motion of the spacecraft is uniquely determined by the initial condition of the orbital dynamics. In the second scenario, the constant unknown maneuver with amplitude of 0.1m / s 2 is applied along the three directions, as stated earlier. Because the nonzero maneuver is completely unknown to the EKF, the estimation performance will be seriously degraded. The position estimation errors of the first scenario are presented in Fig. 3, where the solid lines represent the position estimation errors and the dashed line denotes the ı-boundaries of the estimate errors. The figure shows that at the beginning of the estimation, the estimation error is greater than 100 m, but the estimation accuracy gradually improves over time, and after about 300 s the estimation reaches a stable level and achieves accuracy better than 25 m (ı) along three directions. The result in Fig. 3 also indicates that when the spacecraft is free of unknown maneuvers, the EKF could work accurately for orbital determination.

X-Error(m)

100 Estimate Error σ-Boundary

50 0

0

500

1000 Time(s)

1500

2000

0

500

1000 Time(s)

1500

2000

0

500

1000 Time(s)

1500

2000

Y-Error(m)

100 50 0

Z-Error(m)

100 50 0

Fig. 3. Position estimation errors of the EKF when the unknown maneuver is zero

X-Error(m)

400 Estimate Error σ-Boundary

200 0

0

500

1000 Time(s)

1500

2000

0

500

1000 Time(s)

1500

2000

0

500

1000 Time(s)

1500

2000

Y-Error(m)

400 200 0

Z-Error(m)

400 200 0

Fig. 4. Position estimation errors of the EKF when the unknown maneuver is nonzero When we apply the unknown constant maneuver at the time of t = 500s along each direction within the second scenario, the EKF diverges because it cannot robustly cope with the unknown maneuvers. As shown in Fig. 4, compared with the results in Fig. 3, the estimation errors are at the same level before t = 500s . However, once the unknown maneuver occurs in each direction, the estimation error begins to increase, and finally the estimation diverges. During filter iteration, the Kalman gain and covariance are updated independently of the measurement residual and finally approach a relatively constant value. Even if the unknown maneuver has occurred and the observation residual has been increased, the Kalman gain and covariance still remain relatively constant, and the weight of fresh observation used for estimation correction cannot be improved; hence, the filter diverges. B. VSE Results

1)

Constant unknown maneuver ( u = 0 ) In this section, we demonstrate the estimation performance of VSE when the unknown maneuver occurs. Because the

maneuver observer proposed in Section 3 is sensitive to the time-varying rate of the maneuver, we simulate two different scenarios for comparison. The first scenario considers the spacecraft accelerated by a constant unknown maneuver, while the

second considers a time-varying maneuver. Within the framework of VSE, the maneuver detector plays an important role in changing the structure of the filter. To highlight detector-based switching of the VSE, the pulse-like acceleration, which includes both the start and end during the simulation, is implemented as an unknown maneuver. Without loss of generality, the pulse-like acceleration along the x direction is represented as °­ 0 t ∉[txk , tx ( k +n ) ] ux = ® °¯β x t ∈[txk , tx ( k +n ) ]

(57)

where the subscript x denotes the x direction, β x is the scale factor to regulate the amplitude of the constant maneuvers, t xk and t x ( k + n ) are the start and end time of the maneuver, respectively, and Δt x = t x ( k + n ) − t xk represents the width of the thrust. The constant unknown maneuvers along the y and z directions also have the same expression as Eq. (57). The parameters of the unknown constant maneuver are given in Table 2. Table 2. Parameters of the Pulse-like Acceleration direction

X

y

z

Amplitude(m/s2)

0.01

0.01

0.01

Start time(s)

400

800

1000

Conclusion time(s)

2800

3200

3400

Threshold of start detector

0.1

0.1

0.12

Threshold of end detector

0.7

0.7

0.7

During the simulation, the observer always works and generates the maneuver estimation, even when no maneuver occurs. However, it is only when the maneuver detector declares the maneuver occurrence that the maneuver estimation is fed to the EKF as compensation. The unknown maneuver estimation and actual maneuver are as shown in Fig. 5, where the solid lines represent the actual maneuvers and the dashed lines represent the estimated maneuvers. In this simulation, because the maneuver is supposed to be completely unknown, the initial condition of the observer is assigned according to Eq. (36) by

Z-Maneuver(m/2)

Y-Maneuver(m/2)

X-Maneuver(m/2)

setting u = 0 . 0.2 Maneuver end declaration

0.1 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0.2 Maneuver end declaration

0.1 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0.2 Maneuver end declaration

0.1 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Fig. 5. Real and estimated maneuvers along each direction At the beginning, the estimations follow the actual maneuver approximately, but the estimation still suffers from time delays with respect to the actual maneuver. Moreover, the maneuver estimation is used as the end detector. When the maneuver estimation falls below 0.7, as shown in Fig. 5, the end detector declares the conclusion of the maneuver, and the input compensation along the corresponding direction is switched off. Fig. 6 shows the maneuver estimation errors along three directions. The result clearly indicates that during the estimation, the estimation accuracy is always better than

0.005m / s 2 ;

however, once the unknown maneuver occurs and disappears, the estimation error increases remarkably, which

X-Error(m/s2)

could be treated as a response to the external disturbance caused by the unknown maneuvers. 0.015

0.005

Y-Error(m/s2)

0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0.015 0.01 0.005 0

Z-Error(m/s2)

Estimate Error σ-Boundary

0.01

0.015 0.01 0.005 0

Fig. 6. Estimation errors of the unknown maneuver along each direction During the estimation, the maneuver start detectors are continuously calculated according to Eq. (55). The time history is plotted in Fig. 7, which shows that the detector is relatively stable at a level of less than 0.05 along each direction, but when the unknown maneuver is suddenly applied to the spacecraft, the observational residuals continuously increase. Therefore, the value of the start detector immediately jumps to a peak value greater than 0.2 and consequently declares the occurrence of the maneuver. Then, the structure of the filter is changed and the maneuver estimation is fed into the EKF as input compensation.

x-detector)

0.4 Maneuver Start Declaration

0.2

Threshold

0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

y-detector

0.4 Maneuver Start Declaration

0.2

Threshold

0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

z-detector

0.4 Maneuver Start Declaration

0.2 0

Threshold

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Fig. 7. Maneuver start detectors Once the EKF is compensated by the estimated maneuver, the observational residuals and start detector fall to stable levels, as shown in Fig. 7. The results in Fig. 7 also reveal that the detector proposed in Eq. (55) is adequately sensitive when used to detect unknown maneuvers. The detector also experiences a remarkable jump at the time of maneuver conclusion; this jump can be treated as a response to the maneuver end but not the maneuver occurrence. Figs. 8 and 9 show the VSE estimation errors of the orbital position and velocity results. Before the time of the start maneuver, the orbit of the spacecraft is continuously estimated by the EKF, and an accuracy similar to that represented in Fig.

3 is obtained. When the first maneuver is applied in the x direction (t = 400 s), the detector successfully declares the start of the unknown maneuver. Then, the maneuver estimation generated by the observer is fed into the EKF. According to the VSE framework, when the estimation switching is completed, the covariance of the dynamic uncertainties is immediately reinitialized according to Eq. (53). As a response to the EKF compensation, transient oscillation appears in the estimation result. As shown in Fig. 8, a remarkable transient oscillation occurs at times t = 400, 600 and 800 s along each direction, which is consistent with the time of the unknown maneuver start. The accuracy of the orbit determination, which is represented by the sigma boundary, suffers from degradation after filter switching, but still obtains a better estimate accuracy of 80 m (ı) when it converges to a stable level. Correspondingly, the transient oscillations on the position estimation also occur when the maneuver in terminated at the times of t = 2800, 3200 and 3400 s. The estimation results of the orbital velocity are presented in Fig. 9. Similar to the orbital position estimation, the orbit velocity estimation is also sensitive to the unknown maneuver. Due to the transient response to the maneuver occurrence and conclusion, the average performance is significantly degraded, as shown by the ı boundaries in the plot; however, the maneuver-affected velocity estimation error, and that of the position estimation, are nearly at the same levels as in the non-maneuver-affected stage when the EKF compensation is applied.

X-Error(m)

200

Estimate Error σ-Boundary

100 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Y-Error(m)

200 100 0

Z-Error(m)

200 100 0

Z-Vel-Error(m/s)

Y-Vel-Error(m/s)

X-Vel-Error(m/s)

Fig. 8. Orbit position estimation error in case of constant unknown maneuver 20 Estimate Error σ-Boundary

10 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

20 10 0

20 10 0

Fig. 9. Orbit velocity estimation error in case of constant unknown maneuver

Based on the preceding simulation results, it can be confirmed that compared with the traditional EKF, the proposed VSE can adaptively cope with the unknown constant maneuvers by incorporating the maneuver observer into the estimation. If the transient response of the filter on the unknown maneuver is neglected, the performance of the orbit determination remains nearly the same as the case free of maneuvers.

2)

Time-varying unknown maneuver ( u ≠ 0 ) In this section, we apply the time-varying unknown maneuvers, which are selected as sinusoidal acceleration, to the

spacecraft. The sinusoidal acceleration affects the spacecraft during the whole simulation. The sinusoidal acceleration along the x direction is described as follows: u x = ψ x sin(ωx t + ϕ x 0 )

(58)

where ϕ x 0 is the initial phase, ψ x represents the amplitude of thrust, ω x is used to regulate the time-varying rate of the maneuver thrust and increasing ω x will speed up the time-varying rate of maneuver. Given a sinusoidal acceleration via time, the EKF is compensated only when the amplitude of the maneuver is greater than the predefined threshold. The sinusoidal acceleration along y and z is similar to Eq. (58). The parameters of the sinusoidal thrust are given in Table 3, the eigenvalue of the observer gain matrix is selected as λi = 1.8 . Noted that with a purpose of comparison, the Variable State Dimension estimator (proposed in reference 27) is also simulated. Because the VSD estimator always assumes the unknown maneuver as constant, when the unknown maneuver is time-varying, the estimation performance of VSD estimator will be degraded. Table 3. Parameters of the sinusoidal acceleration Direction

x

y

z

Amplitude (m/s2)

0.01

0.02

0.01

Start time (s)

0

0

0

End time (s)

4000

4000

4000

Thrust frequency (rad/s)

0.001

0.001

0.001

Bias (m/s2)

0

0

0

Initial phase (rad)

0

0

0

The maneuver estimation is plotted together with the actual maneuvers in Fig. 10. Although the maneuver varies during the estimation, the estimated acceleration follows the actual accelerations closely, and no remarkable deviation appears until the simulation is terminated. The estimation error is presented in Fig. 11. As analyzed earlier, when the unknown maneuver is time-varying, the estimation error on the maneuver amplitude is affected by its time-varying rate: the faster the maneuver varies, the greater the estimation errors. Because the sinusoidal acceleration along the y direction varies with the same period but the maximum amplitude, the time-varying rate is the maximum, which leads to the biggest estimation error on the y direction.

X-Maneuver(m/2) Y-Maneuver(m/2) Z-Maneuver(m/2)

0.2 0 Real maneuver Estimated maneuver -0.2 0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0.2 0 -0.2 0 0.2 0 -0.2 0

X-Error(m/s2)

Fig. 10. Actual and estimated maneuvers 0.015

0.005

Y-Error(m/s2)

0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0.01 0.005 0

Z-Error(m/s2)

Estimate Error σ-Boundary

0.01

0.01 0.005 0

Fig. 11. Maneuver estimation error for sinusoidal unknown maneuver The time history of the maneuver start detectors are also presented in Fig. 12. At the beginning of the simulation, the initial value of the applied acceleration is zero; as such, the observational residuals gradually increase together with the amplitude of the acceleration over time. Therefore, the detector also increases gradually to the peak value. Compared with the sharp peaks within the case of constant unknown maneuvers, the detector in this case takes nearly 60 s to reach its peak values, and therefore the maneuvers are declared with an obvious time delay along each direction. Moreover, at the end of one sinusoidal period, the sinusoidal acceleration approaches zero, and accordingly the observer-based EKF compensation is also terminated based on the end detectors. As a response to the estimator structure change, a start detector peak occurs near the end of the sinusoidal period, as shown in the figure.

X-detector

0.2 Threshold

0.1 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Y-detector

0.2 Threshold

0.1 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Z-detector

0.2 Threshold

0.1 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Fig. 12. History of maneuver start detectors Figs. 13 and 14 shows the estimation errors by using the Variable Structure Estimator proposed in this work. At the beginning of the simulation, the observations are processed by the EKF without input compensation. With the growth of the unknown maneuver, the estimation error on the orbit position and velocity rapidly grow via time. However, when the start detectors are greater than the predefined threshold, the estimated maneuver from the observer is fed to the EKF. Thus, the errors on the position and velocity estimation fall rapidly, and after approximately 100 s the error reaches a stable level of less than 100 m and 20m/s, respectively.

X-Error(m)

600 Estimate Error σ-Boundary

400 200 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Y-Error(m)

600 400 200 0

Z-Error(m)

600 400 200 0

Fig. 13. Orbit position estimation error of Variable Structure Estimator As a response to the input compensation, the error jump also appears in the result. For example, the error jumps at the times of t = 100 and 3100 s along the x direction, which is consistent with the time of the start and conclusion declaration. Correspondingly, the performance of the estimation, which is represented by the sigma boundary, is also affected by the filter structure changes, but when the filter converges to the stable level, the average accuracies on the position and velocity estimation are still better than 50 m (ı) and 5 m/s (ı).

X-Vel-Error(m/s) Y-Vel-Error(m/s) Z-Vel-Error(m/s)

40

Estimate Error σ-Boundary

20 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

40 20 0

40 20 0

Fig. 14. Orbit velocity estimation error of Variable Structure Estimator Figs. 15 and 16 represent the results of Variable State Dimension estimator. As stated earlier, the Variable State Dimension estimator only works well when the unknown maneuver is constant. From the results, one can see that the sigma boundary on the position error is always greater than 100m, while the one on the velocity is always greater than 10m/s. Obviously, it can be confirmed that for a time-varying maneuver, the estimation performances on orbit position and velocity are both degraded, and finally achieves a worse result than Variable Structure estimator developed in this work.

X-Error(m)

600 Estimate Error σ-Boundary

400 200 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Y-Error(m)

600 400 200 0

Z-Error(m)

600 400 200 0

Fig. 15. Orbit Position estimation of Variable State Dimension estimator

X-Vel-Error(m/s) Y-Vel-Error(m/s) Z-Vel-Error(m/s)

100 Estimate Error σ-Boundary

50 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

100 50 0

100 50 0

Fig. 16. Orbit velocity estimation error of Variable State Dimension estimator To demonstrate the effect of the time-varying rate of the unknown maneuver, we increase the time-varying rate by improving the frequency of the maneuver from 0.001 rad/s to 0.004 rad/s. The position estimation errors, which are presented in Fig. 17, clearly show that when increasing the time-varying rate of the maneuver, the average estimation accuracy and peaks caused by the filter structure changes become worse, which agrees well with the conclusion from Eq. (52).

X-Error(m)

300 Estimate Error σ-Boundary

200 100 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Y-Error(m)

300 200 100 0

Z-Error(m)

300 200 100 0

Fig. 17. Orbit position estimation error of Variable Structure Estimator ( ω x , y , z = 0.004rad / s ) Moreover, to verify the effect from the eigenvalue of the observer matrix, we decrease λi from 1.8 to 1.2, and the estimation error of the orbit position is given in Fig. 18. Obviously, when decreasing λi , according to Eq. (52), the estimation errors caused by the maneuver varying rate are increased, and compared with the result in Fig. 13 the estimation performance is degraded.

X-Error(m)

300 Estimate Error σ-Boundary

200 100 0

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

0

500

1000

1500

2000 Time(s)

2500

3000

3500

4000

Y-Error(m)

300 200 100 0

Z-Error(m)

300 200 100 0

Fig. 18. Orbit position estimation error of Variable Structure Estimator ( λi = 1.2 )

VI.

Conclusion

Overall, to address the orbit determination problem of spacecraft engaged in unknown, time-varying maneuvers, the VES filter is developed based on the maneuver observer to improve estimation performance. Within the framework of the proposed VSE filter, the EKF could be compensated by the maneuver estimation generated by the observer, which allows estimation convergence even if unknown maneuvers occur. The simulation result indicates that the proposed filter in this work is effective for unknown maneuvers, but with different performances for different cases. When the maneuver start or end is declared, remarkable estimation error jumps appear and degrade the average estimation performance. Furthermore, when we increase the time-varying rate of the unknown maneuver or decrease the eigenvalues of the observer matrix, the estimation performance is also degraded. In the future work, we will incorporate the proposed estimator into the formation maneuvering control, where the leader satellite performs an unknown maneuver, while the follower satellites are required to maintain the configuration of the formation. Due to the robustness to the unknown maneuver, the proposed estimator is expected to enable a closed-loop control with better performance.

Acknowledgment This work was supported by National Natural Science Foundation of China (11102018). The authors thank the Associate Editor and anonymous reviewers for their many helpful comments.

References [1]Baird, M. A., “Maintaining Space Situational Awareness and Taking It to the Next Level,” Air and Space Power Journal, Vol. 27, No. 5, 2013, pp. 50–72 [2]Kiani M, Pourtakdoust S H. Adaptive square-root cubature–quadrature Kalman particle filter via KLD-sampling for orbit determination. Aerospace Science and Technology. 2015, 46:159-167. [3]Calabia A, Jin S. Assessment of conservative force models from GRACE accelerometers and precise orbit determination[J]. Aerospace Science and Technology, 2016, 49:80-87. [4]Kim B Y, Yoon J C, Choi K H, et al. Real time numerical dynamic orbit determination of geostationary satellite for time synchronization service[J]. Aerospace Science & Technology, 2003, 7(5):385-395.

[5]Dahmani M, Meche A, Keche M, et al. An improved fuzzy alpha-beta filter for tracking a highly maneuvering target[J]. Aerospace Science and Technology, 2016, 58:298-305. [6]Kim W H, Han S, Lee J G. “An Optimal FIR Filter With Fading Memory,” IEEE Signal Processing Letters, Vol. 18, No. 5, 2011, pp. 327-330. doi:10.1109/LSP.2011.2131648 [7]Amirzadeh A, Karimpour A. “An interacting Fuzzy-Fading-Memory-based Augmented Kalman Filtering method for maneuvering target tracking,” Digital Signal Processing, Vol. 23, No. 5, 2013, pp. 1678-1685. doi:10.1016/j.dsp.2013.05.002 [8]Liu M, Lai J, Li Z, et al. An adaptive cubature Kalman filter algorithm for inertial and land-based navigation system[J]. Aerospace Science & Technology, 2016, 51(2):52-60. [9] Sacks J., Sorenson H., “Nonlinear Extensions of the Fading Memory Filter,” IEEE Transactions on Automatic Control, Vol. 16, No. 5, 2003, pp. 506–507. [10] Bar-Shalom. “Tracking and Data Association, “Journal of the Acoustical Society of America, Vol. 87, No. 2, 1988, pp. 918-919. [11] Kelecy, T., Jah, M., DeMars, K., “Application of a Multiple Hypothesis Filter to Near GEO High Area-to-mass Ratio Space Objects State Estimation,” Acta Astronautica, Vol. 81, No. 2, 2012, pp. 435–444. [12] Blackman, S. S., “Multiple Hypothesis Tracking for Multiple Target Tracking,” IEEE Aerospace & Electronic Systems Magazine, Vol. 19, No. 1, 2009, pp. 5–18. [13]P. R. Belanger. “Estimation of Noise Covariance Matrices for a Linear Time-Varying Stochastic Process,” Automatica, Vol. 10, No. 3, 1974, pp. 267–275. [14]Hough M E. “Improved Performance of Recursive Tracking Filters Using Batch Initialization and Process Noise Adaptation,” Journal of Guidance Control & Dynamics, Vol. 22, No. 5, 2012, pp. 675-681. [15]Singer R A. “Estimating optimal tracking filter performance for manned maneuvering targets,” IEEE Transactions on Aerospace & Electronic Systems, Vol. AES-6, No. 4, 1970, pp. 473-483. [16]Cloutier, J. R., Lin, C. F., and Yang, C., “Enhanced Variable Dimension Filter for Maneuvering Target Tracking,” IEEE Transactions on Aerospace and Electronics System, Vol. 29, No. 3, 1993, pp. 786–797. [17]J. R. Cloutier, C. F. Lin, and C. Yang. “Maneuvering Target Tracking via Smoothing and Filtering Through Measurement Concatenation,” AIAA Journal of Guidance, Control, and Dynamics, Vol. 16, No. 2, 1993, pp. 377–384. [18]W. D. Blair. “Fixed-Gain Two-Stage Estimators for Tracking Maneuvering Targets,” IEEE Trans. Aerospace and Electronic Systems, Vol.AES-29, No. 3, 1993, pp. 1004–1014. [19]Bar-Shalom, Y., and Birmiwal, K., “Variable Dimension Filter for Maneuvering Target Tracking,” IEEE Transactions on Aerospace and Electronics System, Vol. AES-18, No. 5, 1982, pp. 621–629. [20]H. Lee and M.-J. Tahk. “Generalized input-estimation technique for tracking maneuvering targets,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 35, No. 4, 1999, pp.1388–1402. [21]Lee S C, Liu C Y. “Trajectory Estimation of Reentry Vehicles by Use of On-Line Input Estimator,” Journal of Guidance Control & Dynamics, Vol. 22, No.6, 1999, pp. 808-815. [22]F. Gustafsson. “The Marginalized Likelihood Test for Detecting Abrupt Changes,” IEEE Trans. Automatic Control, Vol. AC-41, 1996, pp. 66–78. [23]A. S. Willsky. “A Survey of Design Methods for Failure Detection in Dynamic Systems,” Automatica, Vol. 12, No. 6, 1976, pp. 601– 611. [24]Ru J, Jilkov V P, Li X R, et al. “Detection of Target Maneuver Onset,” IEEE Transactions on Aerospace & Electronic Systems, Vol. 45, No. 2, 2009, pp. 536-554. [25]Holzinger, M. J., Scheeres, D. J., Alfriend, K. T., “Object Correlation, Maneuver Detection, and Characterization Using ControlDistance Metrics,” Journal of Guidance, Control, and Dynamics, Vol. 35, No. 4, 2012, pp. 1312–1325. [26] Kelecy, T., Jah, M., “Detection and Orbit Determination of a Satellite Executing Low Thrust Maneuvers,” Acta Astronautica, Vol. 66, Nos. 5–6, 2010, pp. 798–809. [27]Gary M. Goff, Johnathan T. Black, and Joseph A. Beck. "Orbit Estimation of a Continuously Thrusting Spacecraft Using Variable Dimension Filters", Journal of Guidance, Control, and Dynamics, Vol. 38, No. 12, 2015, pp. 2407-2420. [28]Li X R, Jilkov V P. “A Survey of Maneuvering Target Tracking—Part IV: Decision-Based Methods,” Proceedings of SPIE - The International Society for Optical Engineering, Vol. 4728, No. 28, 2002, pp. 154-158. [29]Chan Y T, Hu A G C, Plant J B. “A Kalman Filter Based Tracking Scheme with Input Estimation,” IEEE Transactions on Aerospace

& Electronic Systems, Vol. AES-15, No. 2, 1979, pp.237-244. [30]Bogler P L. “Tracking a Maneuvering Target Using Input Estimation,” IEEE Transactions on Aerospace & Electronic Systems, Vol. AES–23, No. 3, 1987, pp. 298-310. [31]Li S, Jiang X, Liu Y. Innovative Mars entry integrated navigation using modified multiple model adaptive estimation[J]. Aerospace Science and Technology, 2014, 39:403-413. [32]Bahari M H, Pariz N. “High maneuvering target tracking using an input estimation technique associated with fuzzy forgetting factor,” Scientific Research & Essays, Vol. 4, No. 10, 2009, pp.936-945 [33]Gholson N H, Moose R L. “Maneuvering Target Tracking Using Adaptive State Estimation,” IEEE Transactions on Aerospace & Electronic Systems, Vol. AES-13, No. 3, 1977, pp. 310-317. [34]Kumar K S P, Zhou H. “A 'current' statistical model and adaptive algorithm for estimating maneuvering targets,” Journal of Guidance Control & Dynamics, Vol. 7, No. 5, 1984, pp. 596-602. [35]Gary M.Goff, JonathanT. Black, and Joseph A. Beck. “Tracking Maneuvering Spacecraft with filter-through approaches using interacting multiple models,” Acta Astronautica, Vol. 114, 2015, pp. 152–163. [36]Xiong K, Wei C L, Liu L D. Robust multiple model adaptive estimation for spacecraft autonomous navigation[J]. Aerospace Science and Technology, 2015, 42:249-258. [37]Jeon D, Eun Y. Distributed asynchronous multiple sensor fusion with nonlinear multiple models[J]. Aerospace Science and Technology, 2014, 39:692-704. [38]Yu M, Chen W H, Chambers J. State Dependent Multiple Model-Based Particle Filtering for Ballistic Missile Tracking in a LowObservable Environment[J]. Aerospace Science and Technology, 2017, 67. [39]Yong H P, Seo J H, Lee J G. “Tracking using the variable-dimension filter with input estimation,” IEEE Transactions on Aerospace & Electronic Systems, Vol. 31, No. 1, 1995, pp. 399-408. [40]Ko H C, Scheeres D J. “Tracking Maneuvering Satellite Using Thrust-Fourier-Coefficient Event Representation,” Journal of Guidance Control Dynamics, Vol. 39, No. 11, 2016, pp. 1-9.