Simultaneous Mean and Covariance Correction Filter for Orbit Estimation

This paper proposes a novel filtering design, from a viewpoint of identification instead of the conventional nonlinear estimation schemes (NESs), to improve the performance of orbit state estimation for a space target. First, a nonlinear perturbation is viewed or modeled as an unknown input (UI) coupled with the orbit state, to avoid the intractable nonlinear perturbation integral (INPI) required by NESs. Then, a simultaneous mean and covariance correction filter (SMCCF), based on a two-stage expectation maximization (EM) framework, is proposed to simply and analytically fit or identify the first two moments (FTM) of the perturbation (viewed as UI), instead of directly computing such the INPI in NESs. Orbit estimation performance is greatly improved by utilizing the fit UI-FTM to simultaneously correct the state estimation and its covariance. Third, depending on whether enough information is mined, SMCCF should outperform existing NESs or the standard identification algorithms (which view the UI as a constant independent of the state and only utilize the identified UI-mean to correct the state estimation, regardless of its covariance), since it further incorporates the useful covariance information in addition to the mean of the UI. Finally, our simulations demonstrate the superior performance of SMCCF via an orbit estimation example.


Introduction
The orbit estimation problem is to obtain an accurate estimation of a space target's (e.g., satellite) position and velocity from noisy observations. The dynamic model of a space target is given below by considering J 2 perturbations [1]:r where r = [x, y, z] T is the position of the space target in the inertial coordinate frame (I-J-K), R = x 2 + y 2 + z 2 , w is the white Gaussian noise process, and F p is the instantaneous acceleration due to the J 2 perturbation [2]. J 2 is the dimensionless second zonal harmonic that quantifies the major oblateness effect of the Earth. Specially, the J 2 perturbation may cause a noticeable precession of low-earth orbit (LEO) satellites, and makes the orbital dynamics model to be a strongly nonlinear function. Thus, solving the orbit estimation problem essentially depends on designing a kind of nonlinear state estimation schemes. At present, nonlinear estimation schemes (NESs) are generally developed from a Bayesian inference framework (BIF) [3] and employed as the widely accepted methods for estimating the space target orbital state. The BIF framework structure is given in Figure 1. Obviously, to implement BIF in order to estimate the orbit state, NES requires one to compute the strongly nonlinear integral of the perturbation. As an example for Gaussian filters [4]: where p(r k−1 |Z k−1 ) is assumed to be a Gaussian distribution, and Z k−1 denotes the sequence of measurements. The numerical computation of such an integral is both complicated and intractable, such that orbit estimation accuracy may be weak and cannot meet real-time performance requirements. For orbit tracking of a space target, a data-starved and uncertain state evolution [5] requires that the tracking algorithm should be rapid and high-precision.
1. Here, data-starved does not mean a lack of measurements. Rather, it implies that due to the large number of space objectives, the collection, association, processing, and fusion of the measurements needs a lot of time. Hence, the measurement update rate and period are relatively elongated. This requires the tracking algorithm deals with data as efficiently and rapidly as possible. 2. Under a data-starved environment or long data update period, as a result of the consistent long-term (usually several orbital periods) propagation of state uncertainties in high fidelity physics models between measurement updates, a state which is initially Gaussian will inevitably become significantly non-Gaussian or uncertain if propagated over a sufficiently long time span. This requires that the tracking algorithm should be accurate when estimating the uncertain or non-Gaussian state.
The traditional NES, based on BIF, does not meet the rapidity and high-precision requirements due to the complicated and intractable nonlinear perturbation integral.
In order to avoid the nonlinear perturbation integral computation, a joint estimation and identification scheme (JEIS) simply views F p as an additional unknown input (UI) and analytically identifies F p to accurately correct the orbit state estimation [6]. Such analytical identification in JEIS decreases the algorithm complexity while improving the orbit estimation accuracy to some degree by feedback correction. In other words, this paper mainly aims to transfer the intractable nonlinear perturbation integral issue into that of perturbation (viewed as UI) identification, to simplify the orbit estimation algorithm complexity and improve its accuracy.
As a classical representative of JEIS, the expectation maximization (EM) algorithm has received continuous attention and research. The EM framework fits in with the requirements of JESI, as shown in Figure 2. However, the existing EM algorithms always model the UI as a constant variable, independent of the state [6], and only identify the mean of the UI. In the J 2 perturbation, F p (viewed as UI in JEIS) depends on the state, so it requires at least the first two moments (FTM), i.e., the mean and covariance. To further improve the accuracy of orbit estimation, it is important to research the novel EM algorithm; to identify the UI-FTM (i.e., perturbation-FTM) which simultaneously corrects the state estimation and its covariance. Unfortunately, the existing EM is incapable of dealing with the case when the UI is associated with the state. Following this idea, that the original orbit estimation is transferred into JESI, we design a novel two-stage EM algorithm (Figures 3 and 4) to deal with the key difficulty of how to simultaneously fit or identify the FTM of UI associated with the state. The main contributions are that: (1) The first EM executes joint orbit state estimation and pseudo measurement (PM) identification using networked multi-sensor observations. PMs can be understood as the indirect reflection of UI from the real observed measurements. (2) The second EM is designed for fitting the UI-FTM by UI-PMs from different sensors, and then used to simultaneously correct the orbit state estimation and its covariance in the first EM. Finally, (3) the identification and fitting computations in these two EMs are analytical, which contributes to rapid and efficient implementation. The simultaneous correction of UI-FTM to the state estimate, as shown in Figure 3, improves the orbit state estimation accuracy. By doing this, a novel simultaneous mean and covariance correction filter (SMCCF) is proposed and achieved. Indeed, vastly different from the standard EM-based filtering scheme, which only estimates the mean characteristic of UI, the novel SMCCF scheme further determines the covariance characteristic of UI. In other words, the joint properties (i.e., both mean and covariance) characterizing the UI in the new scheme provide richer information compared to the single one (i.e., only mean) of the standard approach, which contributes to improved state estimation.
In this paper, the first EM is called the EM-AM since from the signal input viewpoint it mines the actual measurements (AMs) to characterize the UI, as shown in Figure 3. Similarly, we call the second EM as the EM-PM since it utilizes the PM as input to further fit the UI-FTM. However, different from Reference [7], the expanded contributions of this paper lie in that: (1) we specifically define the pseudo measurement and explain its physical implication (see the test following Equation (19)); (2) initialization of EM-AM plays a key role on the SMCCF performance, so we refine initialization of the EM-AM by using a relatively simple forward-backward smoother, instead of the previous two-filters (see Figures 5 and 6, and Remarks 2 and 3); (3) we elaborate the relationship between EM-AM and EM-PM (see Figures 7 and 8), and the algorithm execution program of SMCCF (see Figure 9 in Section V-A), in order to facilitate application of SMCCF in practical engineering; and (4) the performance of SMCCF is thoroughly demonstrated by a space target orbit estimation problem and compared to the present classical methods (see Section VI-B). Also different from other EM-based UI identification issues, for example Reference [6], which are only concerned with the mean of UI, this paper designs the novel SMCCF which simultaneously fits the UI-FTM.

State estimation and its covariance
Joint estimation and identification  This paper is organized as follows. Section 2 formulates the problem. Sections 3 and 4 design the first and second EMs in the two-stage EM algorithm, respectively. Section 5 gives the fusion structure for combining the two EMs and achieving SMCCF. Section 6 demonstrates the superiority of SMCCF to the standard NES, including the extended Kalman filter (EKF) and the cubature Kalman filter (CKF), by an elaborate simulation of a space target orbit estimation. Finally, some conclusions are drawn in Section 7, including the future work.
Notations. The superscripts -1 and T represent the inverse and transpose operations of a matrix, respectively. If X is a positive semi-definite or positive definite matrix, we simply write X ≥ 0 or X > 0. N (x; µ, Σ) denotes the variable x obey a Gaussian distribution with mean µ and covariance Σ. p( ) is the probability, for example p(A|B) describes the conditional probability of the variable A on B. E( ) and E( | ) represent the expectation or conditional expectation, respectively. We define two operations D(x,P) =x T Px and C(x) = xx T . The symbolsˆand˜on top of a random variable, represent an estimate and its error, respectively. For example,x denotes the estimate of variable x and its estimation error isx = x −x. Tr( ) denotes the trace of matrix.

Problem Formulation
Consider a general discrete-time stochastic system with additive UIs in the dynamic and measurement models: where x ∈ R n , y ∈ R m represent the state and measurement vectors, and f (·), andh(·) are the known dynamic and measurement functions, respectively. The matrices M ∈ R n×p , and N ∈ R n×q are known. w k ∈ R n and v k ∈ R m are uncorrelated, zero-mean Gaussian white noise terms satisfying where δ kl is the Kronecker delta function. The initial state x 0 is a Gaussian vector with meanx 0 and covariance P 0 . Here, w k , v k , and x 0 are mutually independent. The parameters a k ∈ R p , and b k ∈ R q are unknown UIs which are associated with the state. Also, we set θ k = {a k , b k } . (1) is common in practical engineering. Consider the example of space object tracking, whose dynamic model is formulated as follows:

Remark 1. The model in Equation
The last term in Equation (2) is known as the J 2 perturbation, F p , which is obviously associated with the state.
If we regard the perturbation as UI, space object tracking can be modeled by Equation (1). Another example is tracking a maneuvering non-cooperative target, whose dynamic model, due to the modeling error, can be represented by: If ∆Fx k−1 is considered as the UI, which is obviously coupled with the state, we can also formulate maneuvering target tracking in terms of the model in Equation (1).
Due to the characteristic that the UI θ k is highly coupled with the system state, θ k possesses the FTM. In this case, if we only consider the UI's mean to correct the state, the estimation accuracy may be undesirable. A feasible method to improve the accuracy is to correct the state by using both the UI-FTM simultaneously. Hence, our objective is to not only identify the mean but also to fit the covariance of the UI. Unfortunately, the existing methods, such as JESI or the standard EM, always regard the UI as a constant. So they all ignore the covariance property of UI, which limits to solving the state estimation problem corresponding to Equation (1).
Combining Equations (1) and (2), we view the perturbation in space target orbit estimation or in tracking as a UI associated with the state. Our aim is to explore the novel EM algorithm for fitting the UI-FTM based on the actual measurements (AMs) Y i k−l:k = {y i k−l , · · · , y i k } (i = 1, 2, · · · , N) from multi-sensors, where l represents the sliding window length. Furthermore, the orbit state estimation is simultaneously corrected by using the fitted UI-FTM.
Different from the classical EM algorithm in Reference [6], a two-stage EM scheme is proposed, as shown in Figure 4 which is a further refinement of Figure 3. In the first stage, for each sensor with Y i k−l:k , EM carries out joint state estimation and UI-PM identification, where PM refers to pseudo measurement. N sensors need N EMs, which run in parallel. All EMs output two groups; the state estimation x 1 k ,x 2 k , · · · ,x N k and the UI-PM identification {θ 1 k ,θ 2 k , · · · ,θ N k }, where {x i k ,θ i k } corresponds to the i-th sensor or the i-th EM. Indeed, the EM in the first stage reflects or characterizes UI by the output PM from the input AM. Thus, we also call the first EM as EM-AM, from a signal input viewpoint. The state estimationx i k in EM-AM running process is only corrected by the PMθ i k . In the second stage, we further fit the UI-PM by a Gaussian mixture distribution with mean µ θ k and covariance P θ k , which can be fed back into the first EM to jointly correct the state estimation and its covariance. Similarly, we name the second EM as EM-PM. By linking these two EMs, the SMCCF is proposed by designing the fusion structure. From a sense of whether enough information is mined or utilized, thex S k by SMCCF should be superior to thex F k by the existing nonlinear filters or algorithms in accuracy, since the former further mines covariance information in addition to the mean of UI. Following Figure 4, Sections 3 and 4 derive the EM-AM and EM-PM, respectively. Section 5 further links these two kinds of EMs, and clarifies the single input-output relationship between them, by designing a new fusion structure. By doing this, the SMCCF is achieved.

EM-AM
The EM algorithm was introduced by Dempster et al. [8] to solve maximum likelihood (ML) estimation with incomplete or missing data. If we consider the state in Equation (1) as the missing data and θ k as the parameters which need to be identified, then the EM framework [9] is suitable for joint state estimation and UI-PM identification precisely.
To implement the EM-AM algorithm in each sensor in Figure 4, we first need to compute the conditional expectation, Q(θ,θ r ), of the complete data log-likelihood function L θ (X k k−l , Y k k−l ) in the E-step. Before addressing this, by first using Bayes' rule and the Markov property of the model in Equation (1), we have: Then, applying the following Gaussian approximation: Equation (4) can be rearranged to yield: where L 0 is a constant independent of θ k , and is thus omitted, and: The following subsections outline the steps to achieve EM-AM according to Figure 2; the E-step and M-step.

E-Step in EM-AM
Taking the expectation of Equation (7) with respect to the probability density function pθ r (X k k−l−1 |Y k k−l ), we can obtain: where I 0 is a constant corresponding to L 0 . Note that this posterior probability, pθ r (X k k−l−1 |Y k k−l ), corresponds to the state estimation on knowing the identification valueθ r at the r-th iteration of EM-AM. It can be computed by Kalman-based estimators, such as Gaussian filters including EKF, and CKF, among others. Thus, we have: where the computations of Equations (13) and (14) refer to: andx i|k−l:k and P i,i|k−l:k are the posterior estimates of the state x i under the measurements Y k k−l from the interval [k − l, k] with the minimum mean square error (MMSE) criterion: In Equations (8) and (12), the terms L 1 and I 1 are associated with the UI θ, but θ is hidden in the probability p θ (x k−l−1 ). Hence, we can not obtain the dominant expressions of L 1 and I 1 like we can for L 2 , L 3 , I 2 , and I 3 . In order to compute the dominant derivative of Q(θ,θ r ) with respect to θ, we have to ignore L 1 and I 1 . But ignoring L 1 and I 1 is reasonable in some sense, because: (1) if l = k − 1, then L 1 equals log p θ (x 0 ), which is independent of θ, and (2) since the UI is associated with the state in our problem, θ is time-varying. In order to quickly reflect or sketch this time-varying UI, it should gradually forget the previous information while preserving the current information in a sliding window of length l. This means that the previous measurements in interval [1, k In order to compute Q(θ,θ r ), the state estimationx i|k−l:k and its covariance P i,i|k−l:k in the interval [k−l, k] must be computed first, which is known as a smooth problem and can be implemented by using a fixed interval smoother. One example is the two-filter smoother [10], which is formulated as follows: where, as shown in Figure 5 the first filter (i.e.,x i|k−l:i , P i,i|k−l:i ) is computed by Kalman-based filtering under the model in Equation (1) and the second filter (i.e.,x i|i+1:k , P i,i|i+1:k ) is also computed with a similar Kalman-based filter but under the inverse form of Equation (1), which runs backwards in time. The unscented Kalman filter (UKS) presented in Reference [11] can be understood as an approximate implementation of this form of smoother in nonlinear systems. Another type of smoother is the forward-backward smoother [12], names as Rauch-Tung-Striebel (RTS) Smoother: Here, the forward pass applies Kalman-based filters to compute the filtering estimatesx i|k−l:i , P i,i|k−l:i and the predictive estimatesx i+1|k−l:i , P i+1,i+1|k−l:i , which are same as those in two-filter smoother. But in the backward pass, the smoothing recursion starts from last time step k and proceeds backwards in time. That is, to compute the smoothing estimation fromx i+1|k−l:k , P i+1,i+1|k−l:k tox i|k−l:k , P i,i|k−l:k , which is different from the second filter in the two-filter smoother. A derivation of the structure in Equation (18), can be found in References [13][14][15], and its computation program is shown in Figure 6. Figure 6. The forward-backward smoother.

Remark 2.
A key point of implementing Equation (17) is initialization. If we regard the interval [k − l, k] as [1, l + 1] (where the initial value isx 0 ), the first filter in Equation (17) should be initialized byx k−l−1|k−l−1:k−l−1 , which is obtained in the proposed SMCCF by knowing the identified UI-FTM in the interval [k − l − 1, k − 1]. In more detail, the first filter is easily initialized by a standard Kalman-based filter with the initial state valuex k−l−1|k−l−1:k−l−1 . The initialization of the second filter depends on that of the first filter. In other words, the second filter also operates a standard Kalman-based filter but under the inverse form of Equation (1) with the initial valuex k|k−l:k , which is obtained by the first filter. The details of the second filter are elaborated in Reference [11]. The inverse computation of Equation (1) complicates the two-filter smoother. Specifically for nonlinear systems, it is intractable and even infeasible to obtain the model's inverse form.

Remark 3.
In the forward-backward smoother, only the forward filter needs to be initialized, which is similar to the initialization of the first filter in the two-filter smoother. While the backward smoothing computation in Equation (18) proceeds backwards in time, i.e., fromx k|k−l:k tox k−l|k−l:k . Note thatx k|k−l:k is obtained by the forward filter. Hence, from the viewpoint of a simpler initialization, the forward-backward smoother outperforms the two-filter one in operation. More importantly, the former contributes to simply achieving the following M-step in EM-AM, without computing the inverse form of the nonlinear model. Due to the above two reasons, this paper implements the forward-backward smoother for the state estimation in EM-AM.

M-Step in EM-AM
So far, we have completed the computation of the expectation Q(θ,θ r ) with respect to θ by the E-step. Before maximizing Q(θ,θ r ), the pseudo measurement (PM) of UI is defined as follows: where θ is coupled to the state, but the inherent coupling characteristic is unknown. Equation (19) indeed implies the posterior identification of UI under actual measurements, i.e., that Y k k−l affects the state estimation pθ r (X k k−l |Y k k−l ) and further acts on the UI identification (θ PM ). Hence PM's definition concretely links estimation and identification in EM-AM.
The M-step requires a maximization of Q(θ,θ r ) with respect to θ. Here, we let the first derivatives of Q(θ,θ r ) with respect to θ to be zero directly: Thus, we have: where: In other words, from Equation (20) to (23), PM can be understood as an average and indirect measurement of UIs by allowing actual observations to further influence UI-PM identification. (22) and (23) indeed determines the estimated valuesâ r+1 andb r+1 at the r + 1-th iteration:

The derivatives in Equations (20) and (21) essentially treat UIs (a and b) as kinds of average values in the interval [k − l, k]. Thus UI-PM identification in Equations
Further, the second derivatives of Q(θ,θ r ) with respect to θ: are strictly negative-definite only if Q i > 0 and R i > 0, which can be always guaranteed in practice. This implies that Equations (22) and (23) uniquely maximize the convex Q(θ,θ r ). However, Equations (22) and (23) are only formally analytical and optimal since they need to approximately compute the following general integrals: where g(x) generalizes the functions f (x i−1 ) and h(x i ), while ω(x) describes pθ r x i−1 |Y k k−l or pθ r x i |Y k k−l . As mentioned in Remark 3, if we use Kalman-based filters to compute such integrals, for example the Gaussian filter, ω(x) denotes the corresponding Gaussian distribution. So numerical approximations such as EKF, and CKF can implement this integral computation (see References [16,17]). To avoid repetition, the process is omitted.
After identifying the UI-PM, how to further fit the UI-FTM with the PM will be dealt with in designing the EM-PM. According to Figure 4, every sensor runs one EM-AM to output a group of {x i k ,θ i k }. In fact, onlyx i k is directly corrected by UI-PMθ i k in the EM-AM. The following section designs the EM-PM for fitting the covariance of UI-PM to correct the estimation covariance P i k .

Remark 4.
At each recursion, the PM parameter is initialized with the previous PM identification, i.e.,θ k witĥ θ k−1 . At the beginning time,θ 0 can be initialized by prior information. For example,θ 0 is computed fromx 0 and P 0 by performing one recursion of the conventional nonlinear filtering, i.e., through the perturbation nonlinear integral to computeθ 0 :θ whose implementation is similar to that in Equation (26).

EM-PM
This section will fit the UI-FTM from the identified UI-PM in the probability domain, not the time-domain, by the EM algorithm. The UI-FTM fitting issue is to find the available probability density p(θ k ) for describing the stochastic property of UI-PM,θ k ={θ i k } N i=1 , where p(θ i k ) corresponds to the i-th sensor. It has been proven that a Gaussian mixture (GM) form can approximate any probability density function as closely as desired [16]. In order to determine the underlying probability distribution p(θ k ) of the identified UI-PM set,θ k ={θ i k } N i=1 at time k, we assume that the stochastic property of UI-PM is generated using M probability distributions where each is a Gaussian function with weight α j , mean µ j , and covariance Σ j : where ρ = α j , µ j , Σ j M j=1 are the associated parameters with UI-FTM fitting, M ∑ j=1 α j = 1, M is the number of Gaussian components in the GM, and p j θ k |µ j , Σ j = N θ k ; µ j , Σ j is a multi-dimensional Gaussian probability density. The UI-FTM fitting is further transferred into the parameter identification of ρ. A detailed discussion for GM fitting based on the EM framework has been proposed in Reference [19]. Here, we only give a brief derivation about how to apply EM for fitting UI-FTM. First, the log-likelihood function of p(θ k ) can be written as: where ϑ j = {µ j , Σ j }. Direct maximization of Equation (29) can not be achieved since we don't know which Gaussian component in the GM is the most suitable for expressing or fitting the probability p(θ i k ) of the j-th UI-PM . For dealing with this key point, a missing variable z = {z i } N i=1 is introduced, which establishes a simple relation between the UI-PM and Gaussian components. Elaborately, z i = {1, 2, · · · , M} Γ and z i = j denotes that the i-th PM,θ i k , is probably characterized by the j-th Gaussian component, i.e., p j θ i k |µ j , Σ j . By doing this, θ k , z constructs the complete-data set for further establishing the complete-data likelihood function: where the subscript z i denotes all of the Gaussian components which the i-th PM,θ i k , may correspond to. Indeed: with j ∈ {1, 2, · · · , M}. It should be emphasized thatθ i k and z i have a one-to-one correspondence, which well explains Equation (30).
Then, the E-step needs to evaluate the conditional expectation of ln p θ k , z|ρ under p z|θ i k ,ρ m : whereρ m is the fitting result of ρ in the m-th iteration. Application of the following equation: to further arrange Equation (31) obtains: with: where d is the dimension of θ, p z i = j|θ i k ,ρ m denotes the estimated probability that the j-th Gaussian component matches the i-th PM,θ i k , under knowingρ m , andα j,m andθ j,m are the m-th identification values of weight and mean/covariance corresponding to the j-th Gaussian component.
In the M-step, by making the following first derivatives as zero: the parameters of ϑ are re-estimated with the updated probabilities: Identifying α requires optimization of the following term of Equation (32): with the constraint M ∑ j=1 α j = 1. However, direct optimization with respect to α is difficult. For simplicity, we equivalently construct: Thus, making the first derivative of Equation (40) with respect to α zero, one obtains: Finally, applying the following equation: to rearrange Equation (41), gives the weight update: Summarizing the above, the fit mean, µ θ k , and covariance, P θ k , of the UI-PM set,θ k , are obtained as follows: which corresponds to the output of the EM-PM in Figure 4. Here,α j,k ,μ j,k , andΣ j,k are the terminatively optimized parameters serving for the j-th Gaussian component at time k.

Remark 5.
EM-PM is initialized by directly computing the stochastic characteristics (average and variance) ofθ k from multi-sensors in time-domain. Although it may be simple and obvious, such the initialization is in line with the stochastic essence. Moreover, the initialization has no effect on the terminative convergence of the EM-PM, since the convex optimization of Equations (36) and (37) uniquely maximize the expectation.

SMCCF
After achieving EM-AM and EM-PM, this section aims to link them and further clarify their input-output relationship with each other. Before doing this, as a comparison, we first show the execution program of directly using the state estimation from the EM-AM output in Figure 7 (as a refinement of the first EM in Figure 4). By using EKF to implement the state estimation computation in the EM-AM, we have: whereθ i k = {â i k ,b i k } denotes the identified UI-PM by the EM-AM corresponding to the i-th sensor, F k−1 and H k are the Jacobian matrices of the dynamic and measurement models. Obviously, in each EM-AM, only the state estimation is corrected by the UI-PM while its covariance computation is the same as that in standard Kalman-based filter. Besides EKF, one has a freedom to select other Kalman-based filters for implementing the state estimation in EM-AM, such as the above-mentioned nonlinear Gaussian filters including UKF, and CKF, among others.
After obtaining N groups of {x i k ,P i k }, we use the well-known Federal structure to directly fuse {x i k ,P i k } as follows: where corresponding to Figure 7, the filter weight β i = P F k (P i k ) −1 . Of course, there also exist other estimation fusion structures, which have been well studied at present. However, this paper mainly wishes to focus on demonstrating the feasibility and superiority of the simultaneous correction forms in the following (Equation (48)), as compared to the standard mean-based correction form in Equations (45) and (46). Thus, we choose the relatively simple Federal structure for trying to best eliminate the possible influence or interference of fusion on the proposed simultaneous correction here. The following continues to combine the EM-AM and EM-PM as shown in Figure 8 (which is a refinement of Figure 4). Different from Figure 7, the SMCCF adds the EM-PM between the EM-AM and Federal fusion structures to fit the UI-FTM. Before carrying out fusion, the state estimation and its covariance are simultaneously corrected by the UI-FTM. This means that if the EM-PM fits µ θ k = µ a k , µ b k and P θ k = P a k , P b k as given in Equations (43) and (44), then employing EKF again as an example, we have: Again the above-mentioned Federal structure in Equation (47) is used to fuse the state estimation and achieve the SMCCF, denoted asx S k andP S k . In fact, Figure 7 is based on the standard filter structure in Equations (45) and (46), where the UI-PM identification valuesθ i k = {â i k ,b i k } are directly used for correcting the state prediction and the measurement prediction, regardless of their covariances. But Figure 8 is based on the SMCCF structure in Equation (48), where {μ θ k , P θ k } are further mined fromθ i k by adding the second EM part to simultaneously correct the predictions and the corresponding covariances.

Remark 6.
Of course, there are other execution programs which link the EM-AM and EM-PM to achieve SMCCF. For example, at each iteration of the EM-AM, the output PM identifications by multi-sensors are fit by EM-PM, and fed back to simultaneously correct the state estimation and its covariance. But, this paper mainly aims to demonstrate the superiority of the novel JESI to the standard NES, and SMCCF to the standard nonlinear filters, including EKF, and UKF, among others. The performance analysis or comparison of different fusion structures will be explored in the future work.

Remark 7.
Maybe, one can find that the state estimation forms between Equations (17), (18), (45), (46), and (48), are slightly different from each other. At each time-recursion, the former is located in the EM-AM by applying these measurements in a sliding window interval [k − l, k] for identifying UI-PM, and further for fitting UI-FTM. Similarly, the latter is located in the fusion center by applying Kalman-based filters for recursively updating the state estimation, which is simultaneously corrected by UI-FTM. Thus, the former needs to be continuously initialized at each recursion as discussed in Remarks 2 and 3, while the latter is only initialized once byx 0 at the beginning time.

Final Algorithm of SMCCF
Referring to Figures 8 and 9, the final implementation algorithm of SMCCF is as follows.
Step 2: The first-stage EM: Joint state estimation and UI-PM identification.
1. Re-initialization-as discussed in Remarks 2 and 3. 2. Parallel EM-AM of multi-sensors for identifying UI-PMs. Given the measurements Y i k−l:k , iteratively compute the analytical UI-PM identification in Equations (22) and (23) until EM-AM terminates (r > r max ), to finally output a group of UI-PM identification valueŝ Parallel state estimation of multi-sensors. Given the full measurements Y i 1:k , implement Kalman-based filters for Equation (1) to recursively output N groups of {x i k , P i k }. As an example the EKF of Equations (45) and (46) (if only corrected by the corresponding UI-PMŝ Note that Y k k−l and Equation (26) only service for UI-PM identification.
Step 3: UI-FTM fitting. Step 4: Fuse {x i k , P i k } by the Federal structure, as given in Equation (47), or other fusion schemes.
Step 5: Repeat Step 2 to Step 4 with time proceeding.

Remark 8.
Compared with the standard NESs, the main factor, affecting the complexity of the SMCCF scheme, Figure 8, lies in the EM-AM iteration times. Owing to the analytical characteristics and to the convex optimization of identifying UI-PM, a direct method of improving EM-AM execution efficiency is to let r max ≡ 1; at each time-recursion, EM-AM is executed only once. However, to guarantee accuracy, the measurement may be extended to many samples. On the other hand, compared with the standard EM which only use the UI-mean to correct the orbit state estimation, the SMCCF scheme is slightly more computational complex because it needs EM-PM to further fit the UI-FTM. Besides the above-mentioned differences, the involved fusion structure and orbit state filtering computation in various algorithms are set to be consistent. This allows us to highlight the superiority of the SMCCF design as compared with the standard NES and EM algorithms, and further facilitates the following demonstrations of the different algorithms under a uniform benchmark by eliminating the interference of other factors.

Simulation
This section considers the space target orbit estimation with the J 2 perturbation. Let the target's position and velocity characterize the system's state, i.e., x = λ ϕ h v λ v ϕ v h T , then the dynamic equation is: where the target velocity is directly disturbed by the J 2 perturbation, which further and indirectly affects the target position. If we view the orbit perturbation as a UI coupled with the target state, then the above dynamic equations can be further simplified and rearranged as: where In Equation (50), F(t) is the transfer matrix of the orbit state, θ denotes the orbit perturbation which can be understood as an unknown input coupled with the state, and M is the multiplicative matrix serving for θ. Indeed, Equation (50) is equivalently derived from Equation (49) and the derivation process is relatively simple so we omit it and only give the result. The reason of deriving Equation (50) from Equation (49) lies in that Equation (50) has the same functional structure as Equation (1), which facilitates application of the SMCCF algorithm for dealing with the practical orbit state estimation issue.
Discretizing Equation (50) by the fourth-order Runge-Kutta yields: where T is the sampling interval. The measurement equation to observe the space target is: where (λ 0 , ϕ 0 , h 0 ) denote the sensor location. By setting different (λ 0 , ϕ 0 , h 0 ), multi-sensors observation system can be established.
The model parameters are set as follows: • The simulated tracking parameters are set as follows: the window length l = 2, five iterations are employed for EM-AM, and 30 for EM-PM, and the number of sensors and GM models are 10 and 5, respectively. Finally, our simulation results are obtained by running ten independent Monte Carlo runs. The initial orbit state estimation is: x 0 = x 0 + [ 400 400 400 0.8 0.8 0.8 ] T P 0 = diag( 400 2 400 2 400 2 0.8 2 0.8 2 0.8 2 ). For spatial scales, define the root mean square errors (RMSE)s of position and velocity at time k as: For the time scale, define the RMSEs as: where N is the independent Monte Carlo running times and K is the simulation sampling length. Figure 10 gives the space target trajectory. In Figures 11-16, all the involved algorithms are obtained by fusing the state estimation results from multiple sensors under the Federal fusion structure. The proposed SMCCF obviously outperforms the other methods, which demonstrates the feasibility and effectiveness of the two-stage EM design scheme in SMCCF, as compared to the standard EM and NES schemes. The J 2 perturbation directly affects and acts on the orbit velocity. Wonderfully, the velocity estimation performance is well improved throughout the entire sampling interval, which just demonstrates the superiority of the novel SMCCF in coping with the perturbation. The above results can also be verified by Table 1, which computes the RMSE pos and RMSE vel . The velocity estimation accuracy performance is improved by 23.4%, as compared to the standard EM, and the position performance is improved by 19.4%.

Conclusions
In the case that a UI is associated with the system state, UI estimation requires the FTM at least. For these dynamic systems with the above-considered UI forms, this paper proposes a novel SMCCF based on a two-stage EM framework, which further improves the state estimation performance by fitting the UI-FTM and using it to simultaneously correct the state estimation and its covariance. Through a space target orbit tracking example, the superiority of SMCCF, as compared to the standard NESs or EM algorithms, is demonstrated.