Impact-Rebound Momentum Excitation Based Inertial Parameters and State Estimation of Defunct Space Object

: Obtaining the inertia tensors of defunct space objects is essential in on-orbit missions. When the inertia tensor of the space object is non-diagonal, the problem becomes challenging. In this case, the system does not have enough information to estimate the six independent parameters of the inertia tensor. In this paper, the problem of estimating the inertial parameters of a defunct space object with non-diagonal inertia tensor is studied. An excitation method of ejecting an impact ball from the tracking spacecraft to the object is proposed to estimate the complete inertial parameters of the object. The impact ball rebounds after colliding with the object and crashes into the atmosphere ﬁnally. After the collision, the angular momentum of the space object changes. The change is used to construct the estimation model. This paper designs an estimation model which consists of two Unscented Kalman-based estimators to estimate the inertial parameters and the motion states of the object. The observability of the estimators is proved through the observability theorem of nonlinear systems. Numerical simulations show that the estimation model is effective in estimating the complete inertial parameters of defunct objects, as well as reducing the measurement errors of the position and attitude of the object.


Introduction
In certain on-orbit operation missions such as capturing defunct objects, the attitude and orbit motions of the spacecraft will be affected by the object [1].A strong impact could damage the operating spacecraft.In this kind of mission, it is crucial to obtain knowledge of the inertial parameters of the defunct object [2,3].The identification of the inertial parameters of defunct space objects has been one of the key technologies for on-orbit service of defunct objects.Identification methods of the inertial parameters of space objects can be divided into kinematics-based methods and dynamics-based methods.These two methods have the same characteristic in building a nonlinear identification model.Due to the inevitable influence of external interference and internal factors, the measurement errors of the measurement system will directly affect the identification accuracy of the inertial parameters.So in identifying of the inertial parameters, the measured values and the inertial parameters are considered together as the state variables to be estimated.It means the identification model is constructed to estimate the motion states and inertia characteristics of the object at the same time.
There are many publications on the inertial parameters identification of defunct space objects.Philip et al. [4] proposed an algorithm to estimate the moment of inertia of the object at the assumption that the defunct object is known with some cooperative characteristic points.Lichter et al. [5] designed a Kalman filter-based estimation algorithm to estimate the standard diagonal inertia tensor of the object.Wilson et al. [6,7] proposed a least squares-based identification model to identify the inertia tensors of defunct objects on the assumption that the estimated mass characteristics of the object are known.Sheinfeld et al. [8] estimated the proportional values between the inertial parameters.Benninghoff et al. [9] proposed an algorithm to identify the center of mass and moments of inertia of defunct objects based on least square estimation.However, all parameters of this algorithm can only be obtained under the condition that the object only has nutation.Pesce et al. [10] obtained the ratio of the moment of inertia and the relative position of the object.Yuan et al. [11,12] proposed a Kalman filter-based algorithm to estimate the ratio of the moments of inertia of defunct objects.Alessia et al. [13,14] used the LIDAR-based pose measurements to estimate the moment of inertia ratios of the object.
Without direct or indirect contact, the inertia tensors of defunct space objects can not be completely identified only by visual measurement [15,16].To fully identify the inertial parameters of the object, the excitation information can be increased by actively changing the motion state of the defunct object.Murotsu et al. [17][18][19][20] proposed the estimation algorithms based on momentum conservation to identify the inertial parameters of the combination of the object and the tracking spacecraft after capturing the object.Leonard et al. [21] identified the center of mass and the inertial parameters of defunct objects by applying control torque to the object directly.Meng et al. [22,23] put forward the idea of directly contacting the object with a flexible rod or indirectly contacting the object with eddy current and magnetic field.By changing the motion state of the object through contact, the excitation information is increased.As a result, the complete inertial parameters of defunct objects are identified.Liu et al. [24] estimated the inertia parameters of the object by contacting softly with the noncooperative object.
In this paper, the problem of estimating the complete inertial parameters of defunct space objects is studied when the inertia tensor of the object is not diagonal.To estimate the inertial parameters of defunct space objects, we proposed a moment-excitation method to increase incentive information in reference [25].By ejecting an impact ball from the tracking spacecraft towards the defunct object, the motion state of the object is changed.The change of the angular momentum of the object is utilized to construct the estimation model.In [25], the impact ball is assumed to adhere to the surface of the defunct object after the collision.The impact and adhesion mechanism of the impact ball are analyzed in [26].As a further study based on references [25,26], this paper investigates the problem of estimating the complete inertial parameters of the object with non-diagonal inertia tensor.An impact-rebound-based estimation model is proposed while the impact-rebound happens between the impact ball and the object.The estimation model consists of two estimators, one is used to filter the measurement errors of the attitude and the position of the object, and the other is used to estimate the complete inertial parameters of the object.The observability of the estimation model will be proved theoretically.
This paper is organized as follows: Section 2 analyzes the problem and introduces the motion models of the object.Section 3 describes the state equation and observation equation of the estimation model.The observability of the estimation model is proved in Section 4. Numerical simulations are described in Section 5. Section 6 presents the conclusion.

Problem Analysis and System Modeling
In this paper, the problem of estimating the inertial parameters of defunct space objects is studied when the inertia tensor of the object is not diagonal.At this time, there are six inertial parameters of defunct objects.According to references [15,16], the six inertial parameters cannot be completely identified without additional excitation.In this paper, an impact-rebound based estimation model is proposed to estimate the inertial parameters of defunct objects based on references [25,26].
Figure 1 shows the schematic diagram of the collision between the impact ball and the defunct space object.The impact ball is ejected from the tracking spacecraft towards the defunct object.Then the ball collides with the object.As proved in reference [26], the relative velocity of the ball should be controlled to less than 82.27 m/s to make sure that the object will not be punctured by the collision between the ball and the object.In addition, the impact ball should collide with the object eccentrically.In this situation, the angular momentum of the space object is changed by the collision.As the tracking spacecraft is close to the object, the influence of gravity gradient can be ignored.In addition, it is assumed that the defunct space object experiences no external torque except the collision between the impact ball and the object.Thus the angular momentum of the system consisting of the impact ball and the object is conserved.The motion states of the object and the impact ball are obtained through the measurement system installed on the tracking spacecraft.It should be noted that this paper focuses on the parameter estimation of the defunct object.The attitude and the position of the defunct object are provided by the measurement system.The problem in this paper can be summarized as follows.Based on the measurements of the state of the space object and the impact ball, an estimation model is designed to estimate the inertial parameters of the defunct object, as well as filtering the errors of the motion state measurements of the object.

Coordinate Frames
At first, the coordinate frames used in the paper are introduced.This paper concerns three objects: the tracking spacecraft, the defunct space object and the impact ball.Three coordinate frames are introduced to describe the motion state: the geocentric equatorial inertial coordinate frame {N − x n y n z n }, the local-vertical, local-horizontal Euler-Hill (LVLH) coordinate frame{L − x l y l z l }, and the body coordinate frame of the defunct object {U − x u y u z u }.The coordinate frames are shown in Figure 1.The origin of the inertial frame {N} is defined at the center of Earth.The axis x n points to the direction of the vernal equinox.The axis z n is the spin axis of the Earth.The axis y n , axis x n , and axis z n form a right-handed orthogonal coordinate frame.The origin of the LVLH frame {L} is defined at the center of mass of the tracking spacecraft.The axis x l coincides with the position vector of the spacecraft and points from the center of mass of Earth to the spacecraft.y l is perpendicular to the axis x l , and is defined in the orbit plane of the tracking spacecraft.The positive direction of y l is along with the direction of the velocity of the tracking spacecraft.The axis z l forms a right-handed system with the x l and y l .The frame {U} is originated at the center of mass of the defunct object.When the three axes of the frame {U} are aligned along the principal inertial axes, the inertia tensor of the object is diagonal.In this paper, the parameter estimation model is designed for the space objects with non-diagonal inertia tensors.The three axes of {U} are not aligned along the principal inertial axes.
The attitude motion equation and the orbit motion equation of the space object are presented in the following context.

Attitude Motion
This section introduces the attitude motion equation of the object.In this paper, the (3-1-2) successive rotation Euler angle is selected to describe the attitude of the object relative to the inertial frame {N}.The Euler angles of the object is defined as: T is introduced to represent the angular velocity of the object in the body coordinate frame {U} relative to the inertial frame {N}.The attitude kinematics equations are given by [27]: where ψu , φu , θu represent the acceleration of the Euler angles ψ u , ϕ u , θ u .
Euler's dynamics equation is written as: where [I u ] is the inertia tensor of the object.[I u ] is defined as: where I = I 1 , • • •, I 6 are the six independent parameters of the inertial tensor [I u ].Assuming that the object is rigid, the derivatives of the inertial parameters are zeros: and M is the total external torque applied in the object.It is assumed that the defunct space object experiences no control force or torque, so M = 0 in this paper.
A state vector composed of the Euler angles and the angular velocity of the object is defined as: The state equation of x q is obtained according to Equations (1) and (2): where f α (x q ) and f ω (x q ) are defined as: Equations ( 1) and (2) represent the attitude motion equations of the object in inertial frame {N}.The attitude of the object is updated by the equations presented in Equations ( 6) and (7).

Orbit Motion
In this section, the relative motion equations of the object in LVLH frame {L} are introduced.The schematic diagram of the relative motion of the tracking spacecraft and the defunct object is shown in Figure 2. The positions of the tracking spacecraft and the object in the inertial frame {N} are noted as R s and R u .The relative position of the object in LVLH frame {L} is defined as r = [x, y, z] T .Assuming that the orbit of the tracking spacecraft is circular, the acceleration of the angular velocity of the tracking spacecraft is zero.ω s , the angular velocity of frame {L} with respect to {N}, is written as: where n 0 is the average orbital angular velocity of the tracking spacecraft.n 0 is given by: where R s is the modulus of R s , and µ represents the gravitational constant of the Earth.The Clohessy-Wiltshire equations are used to describe a simplified model of the relative motion state of the object in LVLH frame {L} [27]: where v = [ ẋ, ẏ, ż] T is the relative velocity of the object in LVLH frame {L}, and a = [ ẍ, ÿ, z] T is the relative acceleration of the object in LVLH frame {L}.
A state vector x p composed of the relative position and the relative velocity of the object is defined as: According to Equation (10), the state equation of x p is: where f r (x p ) and f v (x p ) are defined as: As a conclusion, Equation (10) represents the relative motion equation of the object in LVLH frame {L}.The relative position and the relative velocity of the object are updated by Equations ( 12) and (13).
As seen from the above theoretical analysis, the attitude and position of the object are updated by the attitude motion equation and the orbit motion equation.The core problem of updating the motion state is that the inertia tensor of the object in attitude dynamics equation (Equation ( 2)) is unknown.The inertial parameters of the object cannot be completely identified solely relying on vision.In this paper, an impact-rebound momentum excitation method is designed to estimate the six inertial parameters of the object.Since the system consisting of the impact ball and the object is free from any external force, the system follows the law of angular momentum conservation after the collision.After ejecting the impact ball, an autonomous control mechanism is introduced to compensate the reaction force of the impact ball on the tracking spacecraft.The relative motion of the spacecraft is still in the original orbit after ejecting the impact ball.

Parameter Estimation Filter Design
The attitude motion equations and the orbit motion equations of the defunct object are introduced in the last section.In this section, a parameter estimation model based on the Kalman filter algorithm is designed to estimate the motion state and the inertial parameters of the object.It should be noted that this paper focuses on the problem of estimating the inertial parameters of the object based on the measurements of the motion state.The Euler angles and the relative position of the object are provided by the measurement system installed on the tracking spacecraft.In the simulations, the measurements are set as the sum of the simulated values and noise signals.
The state vector x of the parameter estimation model includes the Euler angles, the angular velocity, the relative position, the relative velocity and the inertial parameters of the object.x is written as: x = [x q , x p , I] T where x q = [α T , ω T ] T represents the Euler angles and the angular velocity of the object, x p = [r T , v T ] T represents the relative position and the relative velocity of the object, T is a vector composed of the inertial parameters of the object.
According to the expression of the state vector in Equation ( 14), the state equation of the estimation model is expressed as: where f q (x(k)) is defined by Equation ( 5), and f p (x(k)) is defined by Equation (11).f I (x(k)) is the derivatives of the inertial parameters of the object.Under the assumption that the object is rigid, the derivatives of the inertial parameters are zeros.The expression of f I (x(k)) is presented as: The general observation equation is written as: where z(k) is the vector of the measurement data.h(x(k)) describes the relationship between the state vector and the measurement data.ξ(k) represents the noise in the measurements.To completely estimate the inertial parameters of the defunct object, three sets of observation equations are established in this paper: two observation equations based on measurements and an observation equation based on angular momentum conservation model.The former two observation equations provide the measurement information about the Euler angles and the relative position of the object, which are widely used in filtering the motion state of the object.The latter observation equation is the core observation equation in this paper, which is obtained by the law of angular momentum conservation.
The specific expressions and meanings of the three observation equations will be given in the following.

Observation Equations Based on Measurements
As mentioned in the previous section, the observation equations are composed of two parts: the measurements-based equations and the angular momentum conservationbased equations.This section introduce the measurements based equations.The measurements based equations include the Euler angles' measurements and the relative position measurements.
The measurements of the Euler angles of the object are noted as α(k).Define z 1 (k) as: Substituting Equation ( 18) into Equation ( 17), yields the expression of h 1 (k): where E 3×3 means a 3 × 3 identity matrix, and O 3×15 represents a 3 × 15 zero matrix.Now the first observation equation is obtained: where z 1 (k) represents the Euler angles of the object obtained from the measurement system, h 1 (x(k)) is calculated through Equation ( 19), and ξ 1 (k) is the measurement error.The measurement of the relative position of the object is noted as r(k).r(k) is chosen as the observation variable z 2 (k): Substituting Equation (21) into Equation ( 17), yields the expression of h 2 (k): Now the second observation equation is obtained: where z 2 (k) represents the measurement of the relation position of the object, h 2 (x(k)) is calculated through Equation ( 22), and ξ 2 (k) is the measurement error.Note that the estimator consisting of Equations ( 14), ( 15), ( 20) and ( 23), which is noted as estimator Σ 1 , is effective in filtering the attitude, the relative position and the relative velocity of the object.However, the inertia tensor of the object with six parameters can not be completely estimated without other excitation information.To increase the amount of independent measurements, an approach of excitation that the tracking spacecraft send out a ball to change the angular momentum of the space object is proposed.An impact ball is ejected from the tracking spacecraft to the defunct object and collides with the object.The inertial parameters of the object become observable after obtaining the change of the state of the object.In this paper, estimator Σ 1 is utilized to filter the motion state of the object before the collision.Another estimator Σ 2 is designed to estimate the complete inertial parameters of the object after the collision, as well as filtering the motion state of the object.As mentioned above, estimator Σ 2 concludes three observation equations.The former two observation equations provide the measurement information about the Euler angles and the relative position of the object, which are expressed in Equations ( 20) and (23).Next, the angular momentum conservation-based equation is introduced.

Observation Equation Based on Angular Momentum Conservation Model
As mentioned above, the impact ball rebounds after the collision between the impact ball and the object.The positions of the impact ball and the object are identified by the measurement system installed on the tracking spacecraft.Due to the assumptions that the gravity gradient torque and atmospheric resistance torque are ignored in the scenario of the close-range tracking missions and there is no active torque applied on the object, the system composed of the ball and the defunct object follows the law of the angular momentum conservation.
The epoch of the collision is noted as t c .The angular momentum of the ball and the object when t ≤ t c are noted as L b (t − c ) and L u (t − c ), respectively.The angular momentum of the ball and the object when t > t c are noted as L b (t + c ) and L u (t + c ), respectively.The angular momentum of the system which consists of the impact ball and the object is noted as L s with respect to the center of mass of the object.In the following, L s is calculated in two phases: t ≤ t c and t > t c .

Before the Collision
The expression of L s is given by: Note that the mass of the impact ball is smaller than the mass of the object.Therefore, the impact ball is regarded as a point when calculating the angular momentum of the ball.
where m b is the mass of the impact ball, r bu (t − c ) is the relative position of the impact ball relative to the object when t ≤ t c , and v bu (t − c ) is the relative velocity of the impact ball relative to the object when t ≤ t c .According to Equation (25), the reference point of the angular momentum of the impact ball is the center of mass of the object.In the inertial frame {N}, the position of the object is variable, which leads to the slight variation of L b (t − c ) within a short phase.In this paper, L b (t − c ) is considered to be a constant vector.The influence of the variation of L b (t − c ) on the estimation model is considered as an error term.L u (t − c ) is the angular momentum of the object before the collision.L u (t − c ) is calculated by: where C nu (t − c ) is the rotation matrix from {U} to {N} when t ≤ t c , and ω u (t − c ) is the angular velocity of the object when t ≤ t c .C nu is calculated by: where c and s are the abbreviations of cos and sin.In the inertial frame {N}, L u (t − c ) remains unchanged before the collision.
The time when the spacecraft sends the ball is noted as t b .Because L b and L u are considered to be invariant in t ∈ [t b , t c ], the measurements of the attitude and the position at any moment can be selected to calculate the angular momentum of the system.Note that although the accurate estimations of the inertial parameters of the object cannot be obtained, the estimator Σ 1 is effective in reducing the measurement errors of the attitude and the position of the object, which will be proved in Section 4. So in the simulations, the values of L b (t − c ) and L u (t − c ) are expressed by the angular momentum of the ball and the object at t = t c .

After the Collision
At the epoch t = t c , the ball collides with the object.The angular velocity of the object exert sudden changes.The measurement system provides the measured values of the motion state of the object.The variation of the angular momentum of the object is used to estimate the six inertial parameters of the object.After the collision, the ball rebounds.The collision between the ball and the defunct object follows the law of angular momentum conservation.Another observation equation is expected based on the rebound momentum excitation.
The angular momentum L s after the collision is calculated by: where where r bu (t + c ) and v bu (t + c ) are the relative position and relative velocity of the impact ball relative to the object when t > t c .Similarly, L b (t + c ) is considered to be a constant value within a short time.The variation of L b (t + c ) is considered as an error term in the simulations.The angular momentum of the object L u (t + c ) is calculated by: where C nu (t + c ) is the rotation matrix from {U} to {N} when t > t c , and ω u (t + c ) is the angular velocity of the object when t > t c .L u (t + c ) remains unchanged after the collision.Substituting Equation ( 24) into (28), yields: The change of the angular momentum of the impact ball caused by the collision is noted as ∆L b : According to Equations ( 25) and ( 29), ∆L b is relevant to the mass, the position and the velocity of the impact ball and the object.In this paper, the filtering of the position and the velocity of the impact ball is not considered.The measurement error of ∆L b is packaged as an error term.In addition, the influence of the measurement error of ∆L b will be discussed in the simulations.Substituting Equations ( 26), ( 30) and (32) into Equation (31), yields: Equation ( 33) is based on the law of angular momentum conservation.∆L b is chosen as the observation variable z 3 (k): Substituting Equations ( 33) and (34) into Equation (17), yields the expression of h 3 : Now the third observation equation is obtained: where z 3 (k) represents the change of the angular momentum of the impact ball obtained from the measurement system, h 3 (x(k)) is calculated through Equation (35), and ξ 3 (k) is the measurement error.

Filter Design
In the previous two subsections, the state equation and observation equation of the estimation model are introduced.This subsection presents the filter process of iteratively updating the state vector of the estimation model.The estimation model established in this paper is based on the Unscented Kalman filter (UKF) model.Compared with the Kalman filter (KF) and the Extended Kalman filter (EKF), the UKF has higher accuracy for systems with strong nonlinearities, and has been widely used to estimate the attitude of spacecrafts [28,29].
The state vector of the estimation model is presented in Equation (14).At first, the UKF transforms the state vector to 2n + 1 Sigma sampling points: where n is the length of the state vector, X(k) = [X 0 , X 1 , • • •, X 2n ] is the vector composed of Sigma points, P is variance of x, and λ is a scaling parameter used to reduce the total prediction errors.
The state equation is given by Equation (15).The Sigma points presented in Equation (38) are propagated through the nonlinear state function f (x): where X i (k) is the i-th column of X(k), and X i (k + 1) is the predicted value of X i (k).
Then the predicted value x (k + 1) and the covariance P (k + 1) of the state vector x is calculated by: where are the weights of the Sigma points and the weights of the covariance matrix P. w m and w c are calculated by: where a controls the sampling state, and β ≥ 0 is a weight coefficient.Then the predicted value x (k + 1) is transformed to a new Sigma point X (k + 1) according to Equation (38).The observation equation is presented in Equation (37).Substituting X (k + 1) into the observation equation, the predicted observation ẑ (k + 1) of x (k + 1) is obtained: where Z i (k + 1) is the predicted observation of X i (k + 1), and X i (k + 1) is the i-th column of X (k + 1).The predicted state vector x(k + 1) is updated by the new observation z(k + 1) and the predicted observation ẑ (k + 1): where K(k + 1) is the gain matrix and is given by: where P z k z k is the covariance of ẑ (k + 1), and P x k z k is the cross covariance between x (k + 1) and ẑ (k + 1).Finally, the posterior covariance P is updated by: In the estimation model proposed in this paper, two estimators are designed: the estimator before the collision Σ 1 is used to estimate the motion state of the object, and the estimator after the collision Σ 2 is designed to estimate the six inertial parameters as well as the motion state.When t ≤ t c , only the measured Euler angles and the measured relative position of the object are provided.Only the scaled values of the moments of inertia can be identified if based solely on vision.The attitude and the position of the object can be estimated according to the motion equations, which will be proved in Section 4. The state vector of Σ 1 contains the Euler angles, the angular velocity, the relative position, the relative velocity and the inertial parameters of the object.The state equation of Σ 1 is Equation (15) and the observation equation is the visual-based observation equation presented in Equations ( 20) and (23).As a summary, Σ 1 consists of Equations ( 14), ( 15), ( 20) and (23).Estimator Σ 1 is designed to estimate the motion state when t ≤ t c .When t > t c , the third observation equation presented in Equation ( 36) is introduced.Equation (36), which is based on the law of angular momentum conservation, is the core of estimating the non-diagonal inertia tensor with six independent parameters.The estimator Σ 2 is composed of three parts: the state vector in Equation ( 14), the state equation in Equation (15), and three observation equations written together in Equation (37).In summary, the measurement system starts working when t > 0. The estimator Σ 1 works when 0 < t ≤ t c , and the estimator Σ 2 starts when t > t c .The algorithm flowchart of the estimators is shown in Figure 3.The initial state vector x 0 consists of α 0 , ω 0 , r 0 , v 0 , and I 0 .P 0 is the initial covariance of the state vector.

Observability Analysis
This Section investigates the observability of the estimation model contained the estimator Σ 1 and the estimator Σ 2 .In this paper, estimator Σ 2 is the core of estimating the complete inertial parameters of the object, as well as estimating the attitude and relative position.Estimator Σ 1 is utilized to filter the measurements of the Euler angles and the relative position of the object.So the estimator Σ 1 is proved to be observable in estimating the motion state which consists of the Euler angles, the angular velocity, the relative position, and the relative velocity of the space object.The estimator Σ 2 is proved to be observable in estimating the motion state and the six inertial parameters.At first, a theorem in references [30,31], which discusses the observability of nonlinear systems, is presented.

Theorem 1.
The following form is used to model the behavior of physical, biological, or social systems, where x ∈ M ⊆ R m is the state vector, y ∈ R n is the vector of observed system outputs, and u ∈ R p is the vector of control inputs.The Lie derivative of h with respect to f at x ∈ M is: where ∂h(x) ∂x is the partial derivative of h(x) with respect to x.Let R be the observability matrix for the nonlinear system, whose rows are formed from the gradients of the Lie derivatives of h(x).The system is locally weakly observable at x 0 if R has full column rank at x 0 .We say that the system satisfies the observability rank condition at x 0 .
According to Theorem 1, the nonlinear system satisfying the observability rank condition is locally weakly observable at x 0 .R 1 and R 2 are introduced to express the observability matrix of estimator Σ 1 and estimator Σ 2 : Each column in the observability matrix R 1 and R 2 represents the partial derivative of each parameter of the state vector x: 1th-3th columns are the partial derivative of α, 4th-6th columns are the partial derivative of ω, 7th-9th columns are the partial derivative of r, 10th-12th columns are the partial derivative of v, and 13th-18th columns are the partial derivative of I.In the following, the observability of the designed estimators is proved.

The Observability of Estimator after the Collision
Note that the dimension of the state vector of the estimator Σ 2 is 18.A sufficient condition of the estimator to be observable is that the observability matrix R 2 is full-ranked as: According to Equation (37), the observation equations of the estimator Σ 2 are presented as: To simplify the expressions of formulas, C, I, ω, B, and b are introduced to substitute 56) in this subsection.Then Equation ( 56) is rewritten as: Taking partial derivatives of h 1 (x), h 2 (x), and h 3 (x) with respect to x, yields: where E 3×3 is a 3 × 3 identity matrix.O m×n is a m × n all-zero matrix.The matrixes D 0 , D 1 , and D 2 represent the partial derivatives of the observation equation h 2 (x) to the Euler angles, the angular velocity and the inertial parameters, respectively.The expressions of D 0 , D 1 , and D 2 are given by: where C 2 , C 1 , and C 3 in Equation (61) represent the transformation matrix corresponding to each rotation of the Euler angles θ, ϕ, ψ, respectively.The expressions of C 2 , C 1 , and C 3 are given by: The expressions of ω * and b * in Equation ( 63) are defined as: where The Lie differentiation of h(x) with respect to the function f (x) is defined by: where D 3 , D 4 , and D 5 represent the Lie differentiation of h 1 (x), h 2 (x), and h 3 (x) with respect to f (x), respectively.The expressions of D 3 , D 4 , and D 5 are presented as: Then the partial derivative of L f h(x) with respect to the state vector x is calculated as: where the symbol " " represents the matrix component which will not be utilized in proving the observability of the system.D 6 -D 10 are calculated by: In Equations ( 73) and (75), D 6 and D 8 are not simplified.This is because D 6 and D 8 contain no inertial parameters.The partial derivative of the terms contained D 6 with respect to the inertial parameters are zeros.∂D 0 f α /∂I and ∂D 1 f ω /∂I in Equation (77) are calculated by: where D 7,i , i = 1, 2, 3 represents the i-th row in D 7 .ω × in Equation ( 79) is defined as: The Lie differentiation of L f h(x) with respect to the function f (x) is defined by: where D 11 represents the Lie differentiation of L f h 1 (x) with respect to f (x).The expression of D 11 is given by: Then the partial derivative of L 2 f h 1 with respect to x is obtained: where D 12 is the partial derivative of L 2 f h 1 with respect to I. D 12 is calculated by: A 21 × 18 matrix F is introduced to represent a submatrix of the observability matrix R 2 : Note that F is composed of the beginning 21 rows of the observability matrix R 2 .If rank(F) = 18 is proved, then Equation ( 53) is proved to be true.In other words, rank(F) = 18 is a sufficient condition for the estimator Σ 2 to be observable.The rank of the identity matrix E 3×3 is 3.For det D 7 = sec x 1 = 0, the rank of D 7 is 3.A 9 × 6 matrix N is introduced: The rank of N is solved by programming in MATLAB.The program is shown in Figure A1.According to the programming results, the rank of N is 5 when ω = b, and the rank of N is 6 when ω = b.It means that if the angular velocity of the object is changed by the ball after the collision, the condition expressed as ω = b is valid.As a result, the rank of N is 6.As E 3×3 , D 7 and N are column full rank matrix, the rank of F is calculated by: In other words, as long as the angular velocity of the object is affected by the collision between the ball, the rank of the observability matrix R 2 is equal to 18.The nonlinear estimator in Equation ( 37) is observable under the circumstance after the collision between the impact ball and the object.

The Observability of Estimator before the Collision
In this subsection, the observability of the estimator Σ 1 is proved.Note that estimator Σ 1 is utilized to filter the Euler angles, the angular velocity, the relative position, and the relative velocity.A sufficient condition of the estimator Σ 1 to be observable is that the observability matrix R 1 is full-ranked as: rank(R 1 ) = 12 (88) A 12 × 12 matrix H is introduced to represent the beginning 12 rows of matrix R 1 .If rank(H) = 12 is proved, then Equation ( 88) is proved to be true.Note that H is the beginning 12 columns in the beginning 12 rows of matrix F. According to Equation (87), the rank of H is calculated by: As a result, estimator Σ 2 is observable in estimating the Euler angles, the angular velocity, the relative position and the relative velocity of the object.
This section proves the observability of the estimation model.Estimator Σ 1 is proved to be observable in estimating the motion state of the defunct object.Estimator Σ 2 is proved to be observable in estimating the motion state and the inertial parameters of the object when the inertia tensor of the defunct object is non-diagonal.Estimator Σ 1 is unobservable for the complete six inertial parameters.When the inertia tensor is diagonal, the inertial parameters only contains the moment of inertial parameters: I = [I 1 , I 2 , I 3 ] T .In this time, estimator Σ 1 is effective in estimating the three principal axes of the inertial parameters.

Numerical Simulations
This section presents numerical simulation examples to demonstrate the performance of the estimation model.The estimation model is composed of two estimators: estimator Σ 1 works when t ≤ t c and estimator Σ 2 works when t > t c .When t < t c , the space object performs an undisturbed attitude motion.At the epoch t = t c , the impact ball collides with the object.When t > t c , the ball rebounds and crashes into the atmosphere finally.Due to the collision, extra information is introduced into the problem.As a result, the inertia tensor of the object becomes observable, as proved in Section 4. The initial conditions and simulation parameters of the impact ball and the space object are shown in Table 1.In addition, the initial value of R s is set as [7178, 0, 0] T km.The initial value of R u is set as [7178.01,0, 0] T km.The sampling frequency of the measurement system is set as 10 Hz.The initial time is set as t 0 = 0 s.The ball is ejected from the tracking spacecraft at t = t b > 0. In the simulations, t b is chosen as 200 s.The ball hits the defunct object and adheres to the object at t = t c , where t c > t b .The value of t c is calculated from the initial parameters, which results in t c = 201 s.
The errors of the estimations of the state variables are defined as: where e α represents the estimation errors of the Euler angles, e ω represents the estimation error of the angular velocity, e r represents the estimation errors of the relative position, e v represents the estimation error of the relative velocity, and e I = [e I 1 , e I 2 , • • •, e I 6 ] T represents the relative errors of the estimations of the inertial parameters.As shown in Equations ( 90) and (91), the absolute errors are utilized to describe the errors of the Euler angles and the angular velocity to avoid zero-crossing issues, while the relative errors are utilized to describe the errors of the inertial parameters as shown in Equation (92).
Table 1.Initial conditions and parameters of the impact ball and the object.

Impact Ball Object
Euler angles (

Estimation Results
As introduced in Section 3, the estimation model proposed in this paper contains two estimators: estimator Σ 1 before the collision and estimator Σ 2 after the collision.This subsection presents the simulation results of the estimation model and analyzes the effectiveness of the estimation model.In the simulation, the measured Euler angles are set as the sum of the true values and zero-mean Gaussian white noise signals with the standard deviation of 1.15 • .The measured relative position is set as the sum of the true value and zero-mean Gaussian white noise with the standard deviation of 0.63 m.The algorithm is initialized as follows.The initial values of the Euler angles are set as the measured values of the measurement system, the initial value of angular velocity is set as [0, 0, 0] T • /s, the initial value of relative position is set as the measured value obtained from the measurement system, the initial value of relative velocity is set as [0, 0, 0] T m/s, and the initial values of the inertial parameters are set as [18,18,18,3,3,3] T kgm 2 /s.
Figure 4 shows the true values of the object's states.Figure 4a,b represent the true values of the Euler angles and the relative position of the object.Figure 4c shows that the angular momentum of the object exerts a change when the collision happens at t = 201 s.The amount of the change of the object is equal to the change of the angular momentum of the ball after the collision.

Before the Collision
This section shows the simulation results of the estimator Σ 1 before the collision.In this phase, the measurements of the Euler angles and the relative position of the object are known.
Figure 5 shows the simulation results of the Euler angles, the relative position, and the inertial parameters of the object.As shown in Figure 5a, the range of the Euler angles is limited to [−180 • , 180 • ]. Figure 5b shows that the errors of the estimated Euler angles are significantly declined within 10 s.At t = 200 s, e α is less than 0.2 • .Figure 5c shows that the estimated errors of the relative position declines from 1.5 m to 0.1 m within 20 s. Figure 5d shows that the average value of e r is less than 0.05 m when t = 200 s.However, the estimator has poor performance in estimating the inertial parameters.The dotted lines and the solid lines indicate the true values and the estimations of the inertial parameters in Figure 5e,f.As shown in Figure 5e,f, the relative errors of the estimations of I j , j = 1, 2, • • •, 6 are more than 100% when t = 200 s.The curves of the estimations of the inertial parameters are divergent in the simulation time.This result demonstrates that the non-diagonal inertia tensor of a defunct object is not observable if only the measured value of the attitude is provided.

After the Collision
This subsection presents the simulation results of the estimator Σ 2 expressed in Equation (37) when t > t c .The initial value of the state vector x in the estimator Σ 2 is set as the estimation of the state vector of the estimator Σ 1 when t = t c .In the numerical simulations in this section, the error of ∆L b given in Equation ( 32) is set as 5 percent of the true value of ∆L b .The estimation results of state parameters are shown in Figures 6-8.
Figure 6 shows the estimation results of the Euler angles and the angular velocity after the collision.Figure 6a shows the estimations of the Euler angles.The solid lines and the dotted lines represent the estimations and the true values of the Euler angles, respectively.As shown in Figure 6b, the estimation of the angular velocity has slight fluctuations in the beginning 200 s.It can be seen from Figure 6c,d that the estimations of the Euler angles converge within 100 s.Note that in the simulations the standard deviation of the errors of the measured Euler angles is set to be 1.15 • .When t > 800, the errors of the estimated Euler angles converge to the region [−0.15 • , 0.15 • ], and the error of the angular velocity converges to [−0.005, 0.005] • /s. Figure 7 shows the estimation results of the relative position and the relative velocity after the collision.As shown in Figure 7a,b, the initial error of the the relative position is around 0.5 m, and the initial error of the relative velocity is more than 0.5 m/s.The estimations of the relative position and the relative velocity converge rapidly within 20 s. Figure 7c,d show the estimation errors of the relative position and the relative velocity.It can be seen from these figures that the relative position converges to 0.05 m within 100 s.When t = 1400 s, the average of the errors of the relative position estimations is less than 0.03 m.The average of the errors of the relative velocity estimations is less than 2 × 10 −4 m/s.This shows that the estimator Σ 2 is effective in filtering the position and the velocity of the object.Figure 8a shows the estimations of the inertial parameters of the object.It can be seen that the estimations of the inertial parameters stabilize within 600 s. Figure 8b corresponds to the errors of the estimations of the inertial parameters.The relative errors of the inertial parameters are reduced from 150% to 8% within 800 s.When t > 1200 s, the estimations of the inertial parameters tend to be stable.When t = 1400 s, the average of the errors of the inertial parameters estimations is within 5%.This result demonstrates that the estimator Σ 2 is effective in estimating the complete inertial parameters of the object.

The Influence of Initial Values of Inertial Parameters
The last subsection presents the simulation results of the estimation model under normal circumstances.This subsection analyzes the influence of the initial values of the inertial parameters I 0 on the estimations of the estimator Σ 2 .In the above simulations, the initial inertial parameters I 0 are set as [18,18,18,3,3,3] T kgm 2 /s.According to Table 1, the true values of the inertial parameters are [8, 7, 6, 1.3, 1.2, 1.1] T kgm 2 /s.The average of the relative errors of the initial inertial parameters is around 145%.In this subsection, the initial inertial parameters I 0 are set as [50, 50, 50, 10, 10, 10] T kgm 2 /s.At this time, the relative error of the initial inertial parameters I 0 is more than 600%.The total simulation time is set as 2000 s.
Figure 9 shows the errors of the estimations after the collision with the large initial error of I 0 .Figure 9a,b show the errors of the Euler angles and the angular velocity estimations of the object.As shown in Figure 9a, the errors of the Euler angles estimations decline from 2 • to 0.2 • within 1000 s.When t > 1000 s, the estimation errors converge to the region [−0.2 • , 0.2 • ].When t = 2000 s, e α is around 0.12 • .Figure 9b shows that e ω converges to the region [−2 × 10 −3 , 2 × 10 −3 ] • /s when t > 1000 s. Figure 9c,d show that the errors of the relative position and the relative velocity estimations decline quickly and converge within 500 s. Figure 9e,f are the errors of the inertial parameters estimations of the object.As shown in Figure 9e,f, the relative errors of the estimations of the inertial parameters decline quickly.Within the simulation time of 1000 s, the curves tend to be converged.The average of the relative error of the inertial parameters is less than 8% when t = 1000 s.The estimation errors of the state variables are tabulated in Table 2.It can be seen from the second line of Table 2 that the estimation model is effective in estimating the state vector of the object when the average of errors of the initial inertial parameters is 600%.The estimation errors of the motion states are basically unchanged compared with the estimations with lower errors of the initial inertial parameters shown in Figures 6 and 7.This demonstrates that the estimation model performs robustly in estimating the motion state and the inertial parameters of the defunct space object.

The Influence of Relative Distance
This subsection analyzes the influence of the relative distance ||r|| on the estimations of the state variables.In the above simulations, the relative distance is set as 10 m.The measurement error of the relative distance is in proportion to the true value of the relative distance.It means that the larger the distance, the larger the measurement errors of the state variables.It can be seen from Equation (32) that, in the estimation model proposed in this paper, the measurement error of the relative distance is reflected in the error of the angular momentum of ∆L b .When the relative distance is set as 10 m, the relative error of ∆L b is assumed as 5%.In this section, the farer relative distance between the object and the tracking spacecraft is considered.The relative distance is set as 100 m.Correspondingly, the relative error of ∆L b is assumed as 25%.
Figure 10 shows the errors of the state estimations when the relative distance r is 100 m. Figure 10a,b are the errors of the Euler angles and the angular velocity estimations of the object.As shown in Figure 10a, the errors of the estimated Euler angles gradually decline within 200 s.e α is around 0.2 • when t = 2000 s.As shown in Figure 10b, e ω is less than 2 × 10 −3 • /s when t = 2000 s. Figure 10c,d are the errors of the relative position and the relative velocity estimations.e r declines from 1 m to 0.05 m and e v declines from 0.1 m/s to 10 −3 m/s within 1000 s.This demonstrates that the estimator after the collision is effective in estimating the motion states when the relative distance between the tracking spacecraft and the object is somewhat large.Compared with the estimations shown in Figures 6 and 7, the errors of the motion state are little larger.However, the relative distance has a slight influence on the estimations of the inertial parameters.Figure 10e,f are the errors of the inertial parameters estimations of the object.As shown in Figure 10e,f, the relative errors of the initial estimations of I j , j = 1, 2, • • •, 6 are more than 100% when t = 200 s.The estimations of the inertial parameters converge to the region [−25%, 25%] within 1000 s.When t = 2000 s, e I j is less than 23%.The estimation errors of the state variables are tabulated in Table 2.It can be seen from the third line of Table 2 that the estimation model is effective in estimating the state vector when e ∆L b is 25%.To further analyze the influence of the relative distance on the estimation model, the relative distance r between the tracking spacecraft and the object is considered when ||r|| ∈ [10,30,50,70, 100] m.Correspondingly, the error of the angular momentum measurement e ∆L b is set as: e ∆L b ∈ [5%, 10%, 15%, 20%, 25%].In this section, only the errors of the attitude and the inertial parameters estimations are considered.In order to show the results conveniently, the root mean square (RMS) of the state variables are defined.Different from e α presented in Equation ( 90), e α,RMS is defined as the relative RMS of the average values of e α in the last N estimations of the Euler angles in the simulation time.Similarly, the statistical RMS error of the angular velocity is defined as e ω,RMS .The statistical RMS error of the inertial parameters is defined as e I,RMS .e α,RMS , e ω,RMS , and e I,RMS are calculated by: In the simulations, N is chosen as 100. Figure 11 is the statistic results of the RMS relative errors of the states with different relative distance.As shown in Figure 11, e α,RMS and e ω,RMS increase slowly with the increase of e ∆L b .Compared with the true values of the Euler angles and the angular velocity, the relative RMS error of the Euler angles is less than 1% and the relative RMS error of the angular velocity is less than 3%.The relative RMS error of the inertial parameters e I,RMS is 5% when the relative distance ||r|| is 10 m, while e I,RMS is up to 25% when ||r|| is 100 m.This result demonstrates that the inertial parameters are more sensitive to the relative distance than the attitude and the position.

Conclusions
This paper proposes a momentum excitation approach which enables the observation of the inertia tensor of a defunct space object with six independent variables.Based on this idea, an impact-rebound based estimation model composed of two estimators is developed.Before the collision happens, the first estimator is engaged to provide a robust estimation of the attitude and the position of the object, while the accurate values of the inertial parameters cannot be obtained.After the collision happens, the second estimator is activated to estimate the attitude, the position and the inertial parameters of the object.The observability of the two estimators is proved theoretically.Numerical simulations illustrate the effectiveness of the parameter estimation model.In addition, the simulations also show that the inertial parameters are more sensitive to the relative distance than the attitude and the position.The estimation model is robust in estimating the state vector under the condition of large initial errors and large relative distance.

Figure 1 .
Figure 1.The schematic diagram of the collision between the object and the impact ball.

Figure 2 .
Figure 2. The relative motion of the tracking spacecraft and the object.

Figure 3 .
Figure 3.The flowchart of the estimator after the collision.

Figure 4 .
Figure 4.The true values of the Euler angles, the relative position, and the angular momentum of the defunct space object: (a) true values of the Euler angles; (b) true value of the relative position; (c) true value of the angular momentum.

Figure 5 .
Figure 5.The estimations of the state variables of the object before the collision: (a) Euler angles estimations; (b) errors of the Euler angles estimations; (c) relative position estimations; (d) errors of the relative position estimations; (e) the moment of inertia estimations; (f) the product of inertia estimations.

Figure 6 .Figure 7 .
Figure 6.The estimations of the attitude of the object after the collision: (a) estimations of the Euler angles; (b) estimations of the angular velocity components; (c) errors of the Euler angle estimations; (d) error of the angular velocity estimation.

Figure 8 .
Figure 8.The estimations of the inertial parameters of the object after the collision: (a) errors of the estimations of the moment of inertia; (b) errors of the estimations of the product of inertia.

Figure 9 .
Figure 9.The errors of the estimations after the collision with large errors of initial inertial parameters: (a) errors of the Euler angle estimations; (b) errors of the angular velocity estimations; (c) errors of the relative position estimations; (d) errors of the relative velocity estimations; (e) errors of the moment of inertia estimations; (f) errors of the product of inertia estimations.

Figure 10 .
Figure 10.The errors of the estimations after the collision at large distance: (a) errors of the Euler angle estimations; (b) errors of the angular velocity estimations; (c) errors of the relative position estimations; (d) errors of the relative velocity estimations; (e) errors of the moment of inertia estimations; (f) errors of the product of inertia estimations.
where x i represents the i-th estimation values among the N estimations of the state variables.

Figure 11 .
Figure 11.The statistical RMS relative errors of the attitude and the inertial parameters of the object under different relative distance ||r||.
(23)ξ 1 (k) are presented together in Equation(20).The expressions of z 2 (k), h 2 (k), and ξ 2 (k) are presented together in Equation(23).The expressions of z 3 (k), h 3 (k), and ξ 3 (k) are presented together in Equation (36).Equation (36) depends on the change of the angular momentum of the object caused by the collision and constitutes the core observation equation proposed in this paper.

Table 2 .
Estimation errors of the state variables when t = 2000 s.