2. DYNAMIC AND MEASUREMENT MODEL

The orbit dynamic model used in this paper is:

(1)$$\dot{{\bi x}}={\bi F}\lpar {\bi x}\comma \; t\rpar +B \cdot {\bi d}\lpar t\rpar $$

where ${\bi x}=\left[\matrix{{\bi r}\lpar t\rpar & \dot{{\bi r}}\lpar t\rpar }\right]^{\rm T}$ is the state of the satellite, including the position and velocity vectors. F(x, t) is the a priori dynamic model of the satellite. d(t) is the unknown manoeuvre acceleration and $B=\left[\matrix{0_{3\times 3} & I_{3\times3}}\right]^{\rm T}$ is the auxiliary coefficient matrix of the manoeuvre acceleration. In this paper, only an impulsive manoeuvre is considered:

(2)$${\bi d}\lpar t\rpar =\left\{\matrix{{\bf 0}\comma \; \hfill & t\, \ne\, t^\ast \hfill \cr \Delta {\bi v}\cdot \delta \lpar t\rpar \comma \; \hfill & t=t^\ast}\right. $$

where Δv is the velocity vector of the impulsive manoeuvre and t* is the occurrence time of the impulsive manoeuvre.

Proper measurements are needed for tracking the manoeuvring satellite. Three radar measurements of range and range-rate are employed in this paper. The measurement data is modelled as:

(3)$$\eqalign{y_{\rho i} & = \Vert {\bi r}-R_i\Vert +\omega_{\rho i} \cr y_{vi} & = \displaystyle{{\lpar \dot{{\bi r}}-\dot{{\bi R}}_i \rpar \cdot \lpar {\bi r}-{\bi R}_i\rpar }\over{\Vert {\bi r}-{\bi R}_i\Vert }} +\omega_{vi}} $$

where r and R_i are the positions of the satellite and the i-th radar station. $\dot{{\bi r}}$ and $\dot{{\bi R}}_{i}$ are the velocities of the satellite and the i-th radar station. y _{ρ i} and y _vi are the range and range-rate measurements of the i-th radar station. ω_ρi and ω_vi are the measurement noises which are considered as zero-mean white Gaussian noises in this paper.

A linear discrete-time system is derived from Equations (1) and (3):

(4)$${\bi x}_k ={\brPhi}_{k\comma k-1} {\bi x}_{k-1} +G_{k-1} {\bi d}_{k-1} $$

(5)$${\bi y}_k =H_k {\bi x}_k +{\bi v}_k $$

where x_k ∈ ℝⁿ is the state of the satellite at the discrete time t _k and ${\brPhi}_{k\comma k-1} \approx I_{6\times 6} \left.+\displaystyle{{\partial {\bi F}}\over{\partial {\bi x}}}\right\vert _{{\bi x}={\bi x}_{k-1}} \cdot \lpar t_{k} -t_{k-1}\rpar $ is the state transition matrix from time t _k−1 to t _k. y_k ∈ ℝ^m is the measurement data, $H_{k} =\left.\displaystyle{{\partial {\bi y}}\over{\partial {\bi x}}}\right\vert _{{\bi x}={\bi x}_{k}}$ is the measurement mapping matrix and v_k is the corresponding measurement error. d_k−1 is the manoeuvre acceleration at time t _k−1. The value of d _k−1 is modelled as:

(6)$${\bi d}_{k-1} =\left\{\matrix{{\bf 0}\comma \; & t^{\ast} \notin \lsqb t_{k-1} \comma \; t_k \rpar \cr \displaystyle{{\Delta {\bi v}}\over{t_k -t_{k-1}}}\comma \; & t^{\ast} \in \lsqb t_{k-1} \comma \; t_k \rpar }\right. $$

The coefficient matrix G _k−1 is defined from Equation (1) as:

(7)$$G_{k-1} =B\cdot \lpar t_k -t_{k-1} \rpar =\left[\matrix{0_{3\times3} \cr I_{3\times3}}\right]\cdot \Delta t_k $$

3. TRACKING MANOEUVRING SATELLITE WITH EKF/IDE METHOD

3.1. Unbiased minimum-variance method for the input and state estimation

According to Gillijns and De Moor (Reference Gillijns and De Moor2007a; Reference Gillijns and De Moor2007b), the unbiased minimum-variance input and state estimation of a linear discrete-time system described in Equations (4) and (5) are proposed as:

(8)$$\hat{{\bi x}}_{k\vert k-1} ={\brPhi}_{k\comma k-1} \hat{{\bi x}}_{k-1\vert k-1} $$

(9)$$\hat{{\bi d}}_{k-1} =M_k \lpar {\bi y}_k -H_k \hat{{\bi x}}_{k\vert k-1}\rpar $$

(10)$$\hat{{\bi x}}_{k\vert k}^{\ast} =\hat{{\bi x}}_{k\vert k-1} +G_{k-1} \hat{{\bi d}}_{k-1} $$

(11)$$\hat{{\bi x}}_{k\vert k} =\hat{{\bi x}}_{k\vert k}^\ast +K_k \lpar {\bi y}_k -H_k \hat{{\bi x}}_{k\vert k}^{\ast}\rpar $$

where $\hat{{\bi d}}_{k-1}$ and $\hat{{\bi x}}_{k\vert k}$ are the estimations of the input and state of the system. The derivation of the gain matrices K _k and M _k is presented in Gillijns and De Moor (Reference Gillijns and De Moor2007a). In this case, it is proven that the estimations for unknown input and state given in Equations (8)~(11) are unbiased and have minimum variance. Therefore, this unbiased minimum-variance input and state estimation method is robust to unknown inputs.

However, it should be noted that this estimation method brings an additional observable requirement compared with the classical EKF method. The gain matrix M _k of the unknown input estimation in Equation (9) is presented as:

(12)$$M_k =\lsqb \lpar H_k G_{k-1} \rpar ^T{\tilde R}_k^{-1} \lpar H_k G_{k-1} \rpar \rsqb ^{-1}\lpar H_k G_{k-1}\rpar ^T{\tilde R}_k^{-1} $$

where ${\tilde R}_{k} =H_{k} \Phi_{k\comma k-1} P_{k-1} \Phi_{k\comma k-1}^{T} H_{k}^{T} +R_{k}$. It can be observed that the estimation of the unknown input is workable only when the Three-Dimension (3D) square matrix $\lpar H_{k} G_{k-1} \rpar ^{T}{\tilde R}_{k}^{-1} \lpar H_{k} G_{k-1} \rpar $ is invertible. Thus, a necessary condition is that the matrix H _k G _k−1 is column full rank. The measurement vector y is divided into two types: the velocity-related measurement y_r and the velocity-independent measurement y_v, where:

(13)$${\bi y}=\left[\matrix{{\bi y}_r & {\bi y}_v}\right]^{\rm T} $$

Then the measurement matrix H _k can be presented as:

(14)$$H_k =\left[\matrix{\displaystyle{{\partial {\bi y}_r}\over{\partial {\bi r}}} & \displaystyle{{\partial {\bi y}_r}\over{\partial {\bi v}}} \cr \displaystyle{{\partial {\bi y}_v}\over{\partial {\bi r}}} & \displaystyle{{\partial {\bi y}_r}\over{\partial {\bi v}}}}\right]=\left[\matrix{\displaystyle{{\partial {\bi y}_r}\over{\partial {\bi r}}} & {\bf 0} \cr \displaystyle{{\partial {\bi y}_v }\over{\partial {\bi r}}} & \displaystyle{{\partial {\bi y}_r}\over{\partial {\bi v}}}}\right]$$

Combining Equations (7) and (14), the matrix H _k G _k−1 can be presented as:

(15)$$H_k G_{k-1} =\left[\matrix{{\bf 0} \cr \displaystyle{{\partial {\bi y}_r}\over{\partial {\bi v}}} \cr}\right]_{m\times 3} \cdot \Delta t_k $$

where m is the dimension of the measurement. It can be observed in Equation (15) that there must be at least three independent velocity measurements to make the matrix H _k G _k−1 column full rank. To satisfy this necessary condition, three ground radar stations are employed in this paper, and the range and range-rate measurements of the ground stations are applied for tracking the manoeuvring satellite.

3.2. EKF/IDE method for manoeuvre tracking

3.2.1. Manoeuvre detection and state estimation

Manoeuvre detection plays an important role in tracking a manoeuvring target. A manoeuvre detection method based on input estimation is proposed in this section.

Combining Equations (4), (5) and (9), the estimation of the unknown manoeuvre acceleration can be presented as:

(16)$$\hat{{\bi d}}_{k-1} =M_k \lpar H_k G_{k-1} {\bi d}_{k-1} +H_k {\brPhi}_{k\comma k-1} \lpar x_{k-1} -\hat{{\bi x}}_{k-1} \rpar +v_k\rpar $$

The unbiased condition which can be seen in Gillijns and De Moor (Reference Gillijns and De Moor2007a) is:

(17)$$M_k H_k G_{k-1} =I_{3\times 3} $$

and yields:

(18)$$\hat{{\bi d}}_{k-1} -{\bi d}_{k-1} =M_k H_k {\brPhi}_{k\comma k-1} \lpar {\bi x}_{k-1} -\hat{{\bi x}}_{k-1}\rpar +M_k {\bi v}_k $$

Equation (18) represents the estimate error of the unknown acceleration, which is affected by the state estimate error ${\bi x}_{k-1} -\hat{{\bi x}}_{k-1}$ and the measurement error v_k. The variance of the acceleration estimate error is derived as:

(19)$$E_{k-1} =E\lsqb \lpar \hat{{\bi d}}_{k-1} -{\bi d}_{k-1} \rpar \lpar \hat{{\bi d}}_{k-1} -{\bi d}_{k-1} \rpar ^T\rsqb =M_k \lpar H_k {\brPhi}_{k\comma k-1} P_{k-1} \Phi_{k\comma k-1}^T H_k^T +R_k \rpar M_k^T $$

To detect the impulsive manoeuvre, a manoeuvre indicator is proposed as:

(20)$$MD_{k-1} =\sqrt {\hat{{\bi d}}_{k-1}^T E_{k-1}^{-1} \hat{{\bi d}}_{k-1}} $$

When the target is not manoeuvring at t _k−1, the true value of the manoeuvre is d_k−1 = 0, therefore Equation (20) can be rewritten as:

(21)$$MD_{k-1} =\sqrt {\lpar \hat{{\bi d}}_{k-1}^ -{\bi d}_{k-1}\rpar ^TE_{k-1}^{-1} \lpar \hat{{\bi d}}_{k-1}-{\bi d}_{k-1}\rpar } $$

which is the Mahalanobis distance of the acceleration estimate error, thus the value of the manoeuvre indicator stays in a certain range. However, when the impulsive manoeuvre occurs at t _k−1 where d_k−1 ≠ 0, the manoeuvre indicator MD _k−1 diverges from the original value and goes beyond the boundary of the Mahalanobis distance. It can be seen that this boundary plays an important role in distinguishing impulsive manoeuvres from the non-manoeuvre period. If the threshold is too high, the sensitivity of manoeuvre detection will reduce. If the threshold is too low, a false-manoeuvre-alarm will easily occur. To simplify the analysis of the threshold, the variance of the acceleration estimate error is considered as a diagonal matrix, meaning that each dimension of the acceleration estimate error is independent of each other.

(22)$$E_{k-1}^\ast =\left[\matrix{\sigma_{k-1.1}^2 & & \cr & \sigma_{k-1.2}^2 & \cr & & \sigma_{k-1.3}^2 } \right]$$

where σ_{k−1, 1}, σ_{k−1, 2}, σ_{k−1, 3} are the main diagonal elements of the matrix E _k−1, representing the standard deviations for each dimension.

Further, in order to obtain the threshold value of manoeuvre detection, its statistical characteristics need to be considered. Figure 1 shows the probability distribution of the acceleration estimate given by Equation (9) in a non-manoeuvre simulation lasting 500 seconds. The red curve in Figure 1 is a fitted Gaussian distribution. It is shown that the probability distribution of acceleration estimate error is in good agreement with a Gaussian distribution. Therefore, the acceleration estimate error in a non-manoeuvre process can be approximated as a random error following the Gaussian distribution.

Figure 1.

The probability distribution of acceleration estimation error in a non-manoeuvre process.

The maximum boundary of the acceleration estimate error in a non-manoeuvre process can be determined as the triple standard deviation as:

(23)$$\hat{{\bi d}}_{k-1}^\ast =\left[\matrix{3\sigma_{k-1\comma 1} & 3\sigma_{k-1\comma 2} & 3\sigma_{k-1\comma 3} }\right]^{\rm T} $$

where $\hat{{\bi d}}_{k-1}^{\ast}$ is the assumed maximum acceleration estimate error. Figure 2 shows the acceleration estimate error and the 3σ boundary in the non-manoeuvre process. It is shown that the 3σ boundary does represent the range that the acceleration estimation error can reach, indicating the validity of Equation (23).

Figure 2.

The acceleration estimation and the 3σ boundary in three axes.

Substituting Equations (22) and (23) into Equation (20), the threshold value of the manoeuvre detection is presented as:

(24)$$MD_T =\sqrt {\lpar \hat{{\bi d}}_{k-1}^\ast \rpar ^{\rm T}\lpar E_{k-1}^\ast\rpar ^{-1}\hat{{\bi d}}_{k-1}^\ast } =3\sqrt 3 $$

MD _T represents the largest value of the manoeuvre indicator MD _k−1 during the non-manoeuvre process. If the value of MD _k−1 is greater than MD _T, it means that an unknown manoeuvre has occurred, causing the manoeuvre indicator to deviate from the original range. Therefore, the manoeuvre estimation of the proposed method can be adjusted as follows:

(25)$$\hat{{\bi d}}_{k-1} = \left\{\matrix{ {\bf 0}\comma \; \hfill & MD_{k-1} \le MD_T \cr \hat{{\bi d}}_{k-1} \comma \; & MD_{k-1} \gt MD_T }\right. $$

When the manoeuvre indicator MD _k−1 is below the threshold MD _T, a manoeuvre does not occur, thus the acceleration estimate $\hat{{\bi d}}_{k-1}$ can be directly set to zero. Conversely, when the manoeuvre indicator MD _k−1 exceeds the threshold MD _T, the acceleration estimate remains the result of the IDE method proposed in Equation (9). With the detection and estimation of impulsive manoeuvres, the state estimation of the target can be proposed. When the manoeuvre is judged as not occurring, a classical EKF method is applied for state estimation. When the manoeuvre is detected as occurring, the state estimation provided by Equation (11) is taken for the manoeuvring target. With the combination of EKF and unbiased minimum-variance input and state estimation method, the manoeuvre tracking method proposed in this paper can provide a good estimation accuracy both in the non-manoeuvre and manoeuvre processes.

A brief summary of the proposed method is shown as follows:

At the k-th time step, for the given state estimation and the covariance at the previous time step $\hat{{\bi x}}_{k-1\vert k-1}$ and P _k−1, and the measurement data at the current time step y_k:

1. 1. Calculate the manoeuvre acceleration estimation of the previous time step $\hat{{\bi d}}_{k-1}$ according to Equation (9).

2. 2. Calculate the manoeuvre indicator MD _k−1 according to Equation (20).

3. 3. Correct the acceleration estimation $\hat{{\bi d}}_{k-1}$ based on the manoeuvre detection as proposed in Equation (25).

4. 4. Obtain the state estimation and the covariance at the current time step $\hat{{\bi x}}_{k\vert k}$ and P _k. If the manoeuvre is detected as not occurring at the previous step, calculate $\hat{{\bi x}}_{k\vert k}$ and P _k by the EKF method. Otherwise, calculate them using Equations (8)~(11).

The computational expense of the proposed method also needs to be considered. If it is assumed that the number of impulsive manoeuvres is far less than the total number of intervals in the entire tracking process, which is quite natural in practice, the computational expense of the fourth step can be considered to be the same as for the EKF method. Therefore, the manoeuvre detection and estimation of the proposed method from the first step to the third step results in additional computational costs compared to the EKF method. Table 1 shows a comparison of the computational time for the EKF method, the ARUKF method and the proposed method over 1,000 time steps in the J2 model. It is shown that the EKF/IDE method has a higher computation time. The computation time of the EKF/IDE method is about 2·5 times that of the EKF method and 1·7 times that of the ARUKF method. However, the computation time of EKF/IDE method is still in the same order of magnitude as the EKF method, and it can meet the needs of real-time manoeuvre tracking.

Table 1.

The computation time comparison for EKF, ARUKF and proposed method

3.2.2. Covariance inflation for the unobserved impulsive manoeuvre

In the previous sections, the measurement data is assumed to exist throughout the entire tracking process. However, the target is not always observable in practice due to the constraints of the ground station. Assume that from the epoch k ₁ to the epoch k ₃, the target is not observed. Then the state estimation of the target can only be obtained by orbit extrapolation due to the lack of measurement data:

(26)$$\hat{{\bi x}}_{k_3 \vert k_3 } ={\brPhi}_{k_1 \comma k_3 } \hat{{\bi x}}_{k_1 \vert k_1 } $$

(27)$$P_{k_3 } ={\brPhi}_{k_1 \comma k_3 } P_{k_1 } {\brPhi}_{k_1 \comma k_3 }^{\rm T} $$

If the target does not manoeuvre in the unobserved arc, state estimation based on the orbit extrapolation is feasible, and the tracking performance of the proposed EKF/IDE method after the unobserved arc is still good. However, if the target manoeuvres during the unobserved arc, there will be a mismatch between the state estimation error and the covariance matrix. For example, assuming that an impulsive manoeuvre ${\bi d}_{k_{2}}$ occurs at the epoch of k ₂ during the unobserved arc, where k ₁ < k ₂ < k ₃, the state estimation error at the end of the unobserved arc is:

(28)$$\hat{{\bi x}}_{k_3 \vert k_3 } -{\bi x}_{k_3 \vert k_3 } ={\brPhi}_{k_1 \comma k_3 } \lpar \hat{{\bi x}}_{k_1 \vert k_1 } -{\bi x}_{k_1 \vert k_1 } \rpar -{\brPhi}_{k_2 +1\comma k_3 } G_{k_2 }{\bi d}_{k_2 } $$

Therefore, there is a mismatch between the state estimation error and the nominal covariance matrix $P_{k_{3} }$:

(29)$${\rm E}\lsqb \lpar \hat{{\bi x}}_{k_3 \vert k_3 } -{\bi x}_{k_3 \vert k_3 }\rpar \lpar \hat{{\bi x}}_{k_3 \vert k_3 } -{\bi x}_{k_3 \vert k_3 } \rpar ^{\rm T}\rsqb -P_{k_3 } ={\brPhi}_{k_2 +1\comma k_3 }G_{k_2 } {\bi d}_{k_2 } {\bi d}_{k_2 }^{\rm T} G_{k_2 }^{\rm T} {\brPhi}_{k_2 +1\comma k_3}^{\rm T} $$

The mismatch will reduce the tracking performance of the proposed method and lead to the errors in manoeuvre detection after the unobserved arc. Thus, an approach of covariance inflation is proposed in this section to correct this mismatch.

At the epoch k ₃ + 1, the unobserved arc ends and the measurement data ${\bi y}_{k_{3} +1}$ is available for tracking. The approach of covariance inflation is applied at this time step as follows:

(30)$${\bf \varepsilon}_{k_3 +1} ={\bi y}_{k_3 +1} -H_{k_3 +1} {\brPhi}_{k_3 +1\comma k_3 } \hat{{\bi x}}_{k_3 \vert k_3} $$

(31)$$\gamma =\displaystyle{{tr\lpar {\bf \varepsilon}_{k_3 +1} {\bf \varepsilon}_{k_3 +1}^{\rm T} -\beta R_k\rpar }\over{tr\lpar H_{k_3 +1} \Phi_{k_3 +1\comma k_3 } P_{k_3 } \Phi_{k_3 +1\comma k_3 }^T H_{k_3+1}^{\rm T} \rpar }} $$

(32)$$P_{k_3} = \left\{\matrix{ P_{k_3}\comma \; & \gamma \lt 1 \cr \gamma P_{k_3} &\comma \; \gamma \gt 1 }\right. $$

γ in Equation (30) represents the covariance inflation factor. ${\bf \varepsilon}_{k_{3} +1}$ is the measurement residual of the orbit prediction and R _k is the variance matrix of the measurement error. β is the softening factor to eliminate the interference of the measurement error. In Equation (30), the numerator shows the covariance matrix of the actual measurement residual and the denominator shows the nominal variance of the measurement residual based on orbit extrapolation. Therefore, γ shows the ratio of the actual measurement residual to the nominal covariance. If the ratio is greater than one, it indicates that there is a mismatch between the state estimation error and the covariance matrix at epoch k ₃. The covariance matrix is scaled by the factor γ, and can then be used in the state estimation proposed in the previous section. With the inflation of the covariance, the mismatch caused by the unobserved manoeuvre is tolerated, and the tracking performance of the proposed method after the unobserved arc is improved, as demonstrated by the simulations in the next section.

4. SIMULATIONS AND VALIDATION

To illustrate the feasibility of the proposed method, several simulations are presented in this section. In the simulations. the initial true position and velocity of the satellite are [−5,726·21, −3,472·06, 1,156·51] km and [−1·6, 0·15, −7·49] km/s. The initial state estimation errors are [1, 1, 1] km for the position and [1, 1, 1] m/s for the velocity. The standard deviations of the range and the range-rate measurement error are 10 m and 0.1 m/s. The tracking process lasts for 1,000 s in cases 1 to 5 and lasts for 2,000 s in case 6. The measurement sampling step is 1 s. It is assumed that the manoeuvring target has been identified and the tracking association has been done properly. The results of the adaptively robust unscented Kalman filter (ARUKF) method proposed by Wang et al. (Reference Wang, Sun and Li2014) are presented for reference.

In cases 1, 2 and 3, the single impulsive manoeuvre situation is demonstrated. The impulsive manoeuvre occurs at epoch 500 s. The characteristic velocities of the manoeuvres are 100 m/s, 10 m/s and 1 m/s, respectively, in those three cases, and the manoeuvre vectors are along the direction of the satellite's velocity.

Figure 3 illustrates a comparison between the proposed EKF/IDE method and the ARUKF method in case 1. Figures 3(a) and 3(b) show the position and velocity estimation errors of the methods. Figure 3(c) shows the manoeuvre indicator proposed in Equation (20), and Figure 3(d) shows the manoeuvre estimation provided by the proposed EKF/ IDE method. In Figures 3(a) and 3(b), it can be seen that the proposed EKF/IDE method has almost the same tracking performance as the ARUKF method. After the manoeuvring instance, the EKF/IDE method maintains the accuracy of the position estimation, and the velocity estimation error shows a small jump. The state estimation error of the ARUKF method shows a huge increase at the manoeuvre moment. The largest position and velocity estimation errors of the ARUKF method are 1,200 m and 100 m/s, while the largest estimation error of EKF/IDE method is 12 m and 0·5 m/s after the manoeuvre. However, the state estimation error of the ARUKF method decreases rapidly in this case, and there is almost no difference between the results of these two methods after 550 s. Figure 3(c) shows that the manoeuvre indicator of the EKF/IDE method is always below the threshold during the non-manoeuvre arc and it is beyond the threshold at 500 s which is exactly the moment of occurrence of the impulsive manoeuvre. At this moment, the value of the manoeuvre indicator is about 300 which is much larger than the threshold, and this indicates the occurrence of the impulsive manoeuvre. As shown in Figure 3(d), the proposed EKF/IDE method provides an estimation of the impulsive manoeuvre which is not provided by the ARUKF method.

Figure 3.

Tracking performance of EKF/IDE and ARUKF in case 1.

A characteristic velocity of 100 m/s is applied in case 1. In this condition, the advantage of the EKF/IDE method compared to the ARUKF method is not obvious. Therefore, the smaller characteristic velocities are applied in cases 2 and 3. The results are illustrated in Figures 4 and 5.

Figure 4.

Estimation performance of EKF/IDE and ARUKF in case 2.

Figure 5.

Estimation performance of EKF/IDE and ARUKF in case 3.

In case 2, the characteristic velocity of the manoeuvre is 10 m/s. The largest position and velocity estimation errors of ARUKF method are 175 m and 10 m/s, while the largest estimation errors of EKF/IDE method are 10 m and 0·4 m/s. The proposed EKF/IDE method has a better performance than the ARUKF method from the occurrence of the manoeuvre to the end of the tracking process, which is shown in Figures 4(a) and 4(b). The impulsive manoeuvre is also detected and determined accurately as shown in Figures 4(c) and 4(d).

In case 3, the characteristic velocity of the impulsive manoeuvre is 1 m/s. It is shown in Figures 5(a) and 5(b) that the largest position and velocity estimation errors of the ARUKF method are 100 m and 1 m/s while the largest errors of EKF/IDE method are 18·7 m and 1 m/s during the tracking process. The ARUKF method has a large state estimation error from the occurrence of the impulsive manoeuvre, and the convergence speed of the ARUKF method is quite slow. The proposed EKF/IDE method shows a better state estimation result and a faster convergence speed than the ARUKF method. With the decrease of the characteristic velocity of the manoeuvre, the difficulty of manoeuvre detection increases. The ARUKF method detects the impulsive manoeuvre at 848 s, which means the manoeuvre detection delay of the ARUKF method is 348 s in this situation. Meanwhile, the manoeuvre detection delay of EKF/IDE method is 18 s as shown in Figures 5(c) and 5(d), which shows the better manoeuvre detection ability of the proposed EKF/IDE method. It should be noted that the maximum value of the manoeuvre indicator during the whole tracking process is just slightly larger than the threshold value, as shown in Figure 5(c). This shows that a single impulsive manoeuvre with a characteristic velocity smaller than 1 m/s is quite difficult to detect accurately by the proposed EKF/IDE method with the given measurement error. A series of Monte Carlo experiments with different characteristic velocities are implemented to show the manoeuvre detection ability of the proposed method more clearly. The results are shown in Table 2.

Table 2.

Monte Carlo experiments of different characteristic velocities

In Table 2, six groups of Monte Carlo experiments are presented, and the characteristic velocities of the impulsive manoeuvre are set as 10, 5, 2, 1, 0·8 and 0·5 m/s respectively. With the decrease of the characteristic velocity, the average error and the standard deviation of the manoeuvre estimation and the detection delay increase. When the characteristic velocity is smaller than 1 m/s, the manoeuvre estimation has quite a large relative error, and the manoeuvre detection delay rises significantly, which shows that it is difficult to obtain an accurate estimation and detection of an impulsive manoeuvre by the proposed method under this situation. The rate of the missing detection is also a problem when the characteristic velocity is too small. As shown in Table 2, there are three instances of missing detections over 1,000 experiments when the characteristic velocity is 0·8 m/s, and the number of missing detections increases to 115 when the characteristic velocity is 0·5 m/s. This shows that manoeuvre detection easily fails when the characteristic velocity of the impulsive manoeuvre is smaller than a minimum limit of 0·5 m/s under the measurement noise level used in this paper.

In the first three cases, the situations of the single impulsive manoeuvre are tested. After that, the situations of the multiple impulsive manoeuvres are presented.

In case 4, five impulsive manoeuvres occur at 300, 400, 500, 600, 700 s, and the characteristic velocity of each manoeuvre is 5 m/s. Figure 6 demonstrates the tracking performance of the EKF/IDE method and the ARUKF method in this case. Figures 6(a) and 6(b) show that the proposed EKF/IDE method has a better tracking performance than the ARUKF method during the manoeuvre process. The position estimation error of the ARUKF method shows a significant fluctuation with a maximum value of 120 m, and the velocity estimation error of the ARUKF method keeps above 1 m/s until the end of the manoeuvre process, while the EKF/IDE method provides a stable tracking result with a smaller estimation error. Figure 6(c) shows the manoeuvre detection result of the EKF/IDE method with which five impulsive manoeuvres are detected accurately without detection delay. Figure 6(d) shows the manoeuvre estimation error of the EKF/IDE method. The manoeuvre estimation error is basically consistent with the result of the single impulsive manoeuvre situation as shown in Table 2.

Figure 6.

Estimation result of EKF/IDE and ARUKF in case 4.

When the interval and the characteristic velocity of the impulsive manoeuvres are as large as we set in case 4, the estimation result of the EKF/IDE method is similar to a single manoeuvre situation because the proposed method can be fully convergent between two adjacent manoeuvres. Thus, the situation with a smaller interval and characteristic velocity is tested in the next simulation.

In case 5, there are ten impulsive manoeuvres in the period 500–509 s. The interval between two adjacent manoeuvres is 1 s and the characteristic velocity of each manoeuvre is 1 m/s. The sum of the manoeuvres is 10 m/s, which is similar to case 2. Thus, the state estimation result of EKF/IDE method in case 5 is similar to the result in case 2, as shown in Figures 7(a) and 7(b). However, the result of the manoeuvre estimation and detection shows some differences in case 5. Figures 7(c) and 7(d) demonstrate the manoeuvre detection and estimation result in case 5. There are ten impulsive manoeuvres during the whole tracking process, but only three manoeuvres were detected at 501, 504 and 508 s with estimated magnitudes of 2·5248 m/s, 2·1130 m/s and 2·3813 m/s. The last manoeuvre at 509 s was not detected because the velocity covariance at this moment was large enough to cover the measurement residual caused by this last manoeuvre. Although the manoeuvre estimations were not identical to the true manoeuvres, the proposed EKF/IDE method still provides state estimation with an acceptable precision, which shows its good manoeuvre tracking ability in this situation.

Figure 7.

Estimation result of EKF/IDE and ARUKF in case 5.

In cases 1 to 5, all the manoeuvres occur when the measurement is available, which is actually an ideal situation. Thus, in case 6, the manoeuvre is set to occur in an unobserved arc. The tracking process in case 6 lasted for 2,000 s, and the measurement data from 600 s to 1,200 s was not available. Three manoeuvres occurred at 300 s, 900 s and 1500 s, and the characteristic velocity of each manoeuvre was 10 m/s. Figure 8 demonstrates the result of this simulation. At 300 s and 1,500 s, the tracking performance of the proposed method was similar to case 2 because those two manoeuvres occurred when the measurement was available, as shown in Figures 8(a) and 8(b). However, the second manoeuvre occurred in an unobserved arc at 900 s. Therefore, the state estimation error of the proposed method increased rapidly from 900 s to 1,200 s because of the lack of the measurement. At 1,201 s, just after the end of the unobserved arc, a manoeuvre was detected by the proposed method, as shown in Figure 8(c). The manoeuvre estimation at this moment was close to 10 m/s, as shown in Figure 8(d), which indicates the manoeuvre in the unobserved arc. After this moment, the state estimation error of the proposed EKF/IDE method decreased rapidly, and the tracking process worked effectively. It should be noted that the manoeuvre estimation after the unobserved arc is not always close to the true unobserved manoeuvre. In case 6, the length of the unobserved arc was only 10 minutes. If the length was much longer than that, for example, 1 or 2 hours, the estimation accuracy of the unobserved manoeuvre will be influenced significantly because of the extrapolation of the state estimation error during the unobserved arc. However, the state estimation after the unobserved arc is always effective even in that type of situation, which relies on the proposed covariance inflating approach.

Figure 8.

Estimation result of EKF/IDE and ARUKF in case 6.

Finally, in case 7, an orbit determination process based on actual measurement data was demonstrated to validate the proposed method. The tracking period started from 10:00:00 on 1 June 2016. The entire measurement period lasted 900 seconds and the measurement worked at a frequency of 10 Hz. The position of ground stations were: h ₁ = 1·1 km, ϕ₁ = 38 · 6^°, λ₁ = 87 · 2^°; h ₂ = 1.2 km, ϕ₂ = 40 · 1^°, λ₂ = 98 · 2^°; h ₃ = 1·0 km, ϕ₃ = 35 · 4^°, λ₃ = 97 · 6^°, where h _i, ϕ _i, λ_i are the height, latitude and longitude of the i-th ground station respectively. The covariance matrix of measurement noise was set as: R = diag[(10 m) ², (10 m) ², (10 m) ², (0 · 1 m/s) ², (0 · 1 m/s) ², (0 · 1 m/s) ²]. In this actual measurement case, unlike in the previous simulation cases, due to the lack of actual position and velocity data of the target, the measurement residual Δ y was used to verify the accuracy of the estimation:

(33)$$\Delta {\bi y}_k ={\bi y}_k -h\lpar \hat{{\bi x}}_k \comma \; t_k\rpar $$

where y_k is the actual measurement, and $h\lpar \hat{{\bi x}}_{k} \comma \; t_{k} \rpar $ is the nominal measurement obtained from the state estimation $\hat{{\bi x}}_{k}$. When the measurement residual remains within triple standard deviation of the measurement error, the state estimation can be considered to be accurate. If the measurement residual exceeds triple the standard deviation of the measurement error, it indicates that the state estimation deviates from the actual value.

Figure 9 shows measurement residuals of the EKF, ARUKF and proposed EKF/IDE methods. In order to make the image more concise and show the trend of measurement residuals with time more clearly, the sampling interval of measurement residuals in Figure 9 is set as 5 seconds.

Figure 9.

Measurement residuals of EKF, ARUKF and EKF/IDE in case 7.

In Figure 9, The residuals of EKF and ARUKF exceeded the 3σ boundary at 385 s, which indicates that an unknown manoeuvre occurred. The residuals of the EKF method continued to increase, indicating that its filtering diverged. At 480 s, the ARUKF method detected the manoeuvre and inflated the covariance, therefore its measurement residuals stopped growing. However, it can be observed in Figures 9(b), 9(d) and 9(f) that the range-rate measurement residuals of the ARUKF method did not return to near zero, but remained at a higher value, and still exceeded the 3σ boundary many times from 600 s to the end of the tracking process. This indicates that the ARUKF method has a long manoeuvre detection delay and takes a long time to re-converge the filter after detecting the manoeuvre. The proposed EKF/IDE method has good robustness against unknown manoeuvres, as shown in Figure 9. The measurement residuals of the proposed method were well maintained within the 3σ boundary throughout the entire tracking process, except for a few points, which were caused by measurement errors of more than three standard deviations. Figure 10 shows the manoeuvre indicator of the proposed EKF/IDE method. The manoeuvre indicator exceeded the threshold for the first time at 368·7 s, indicating the unknown manoeuvre, while the ARUKF method detected it at about 480 s. Then the manoeuvre indicator again exceeded the threshold at 548·8 s, which was caused by the unknown thrust during the period from 368·7 s to 548·8 s. Then, the indicator remained below the threshold due to the end of the manoeuvre. Overall, the proposed EKF/IDE method demonstrates its feasibility and validity in this example with real measurement data. The manoeuvre detection delay was lower than that of the ARUKF method, and the state estimation accuracy during and after the manoeuvres was also better.

Figure 10.

Manoeuvre indicator of the EKF/IDE method.