# Relative Status Determination for Spacecraft Relative Motion Based on Dual Quaternion

Authors: Jun Sun; Shijie Zhang; Xiande Wu; Fengzhi Guo; Yaen Xie
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

For the two-satellite formation, the relative motion and attitude determination algorithm is a key component that affects the flight quality and mission efficiency. The relative status determination algorithm is proposed based on the Extended Kalman Filter (EKF) and the system state optimal estimate linearization. Aiming at the relative motion of the spacecraft formation navigation problem, the spacecraft relative kinematics and dynamics model are derived from the dual quaternion in the algorithm. Then taking advantage of EKF technique, combining with the dual quaternion integrated dynamic models, considering the navigation algorithm using the fusion measurement by the gyroscope and star sensors, the relative status determination algorithm is designed. At last the simulation is done to verify the feasibility of the algorithm. The simulation results show that the EKF algorithm has faster convergence speed and higher accuracy.


---

## Body

## 1. Introduction

There are higher demands for the satellite navigation technology and control accuracy in space missions than before. The accuracy of satellite attitude determination is the key factor for the performance of the attitude control. Hence, how to effectively improve the accuracy of satellite attitude determination has received increased attention during the last years.With rapid multitarget acquisition and pointing and tracking capabilities, agile spacecraft has great practical value and wide prospect for application. Moreover, agile spacecraft has many advantages such as lighter mass, smaller size, cheaper price, and shorter development cycle, which makes great differences to space technology compared with classical spacecraft.However, it is difficult to establish relative kinematics and dynamics model for double agile satellites’ collaborated earth observing mission because of the rapid attitude maneuverability of agile satellites. Classical modeling approaches separately address the orbit and attitude problems, which complicates models and the spacecraft attitude control. To describe spiral motion, dual quaternion was proved to be the most compact and efficient tool among other approaches. Considering the advantages of dual quaternion in describing spiral and relative motion, building double agile satellites integrated orbit and attitude dynamics equations not only can enhance the computing efficiency but also may avoid the effect of orbit and attitude decoupling.The Kalman Filter which has many advantages such as fast calculation, small memory resource, and real-time is a recursive algorithm based on minimum mean-square error [1, 2]. The basic Kalman Filter is limited to a linear assumption. The spiral motion of spacecraft in practical applications, however, can be nonlinear. As a result, spacecraft kinematics and dynamics equations are nonlinear [3–5]. In view of the effect of nonlinearity, this paper, using EKF (Extended Kalman Filter) technology, combined with the dual quaternion integrated dynamics modeling, designs a navigation algorithm of gyroscope/star sensor combined measurement and carries on the simulation to verify the feasibility of the algorithm.The paper is organized as follows. Section2 describes the relative kinematic and dynamic equations based on dual quaternion. Section 3 designs the measurement models of gyroscope/star sensor, because satellite attitude determination mostly adopts gyroscope/star sensor cooperative determination. Based on the models designed, the state and measurement equations which are linearized in Section 4 are derived. Section 5 shows simulation results. Finally, Section 6 contains the main conclusions.


## 2. The Relative Dynamic Equations Based on Dual Quaternion

### 2.1. Kinematic Equation Based on Dual Quaternion

Figure1 shows two rigid-body spacecraft orbiting the earth: one is spacecraft 1 and the other is spacecraft 2. Their body fixed frames, respectively, are Oa-xayaza and Ob-xbybzb. We can see the relative motion of the two spacecraft as the spiral motion of the Ob-xbybzb relative to the Oa-xayaza. The relative motion is described by dual quaternion; the expression is
(1)q^ba=qba+ε12rbaa∘qba=qba+ε12qba∘rbab,
where qba is attitude quaternion of the spacecraft 1 relative to spacecraft 2 and rba=rb-ra is relative position vector of two-spacecraft center mass.Figure 1
The orbit of the relative motion.The kinematic equation for the relative motion of two spacecraft expressed by dual quaternion is governed by(2)2q^˙ba=ω^baa∘q^ba=q^ba∘ω^bab,
where
(3)ω^baa=ωbaa+εvbaa=ωbaa+ε(r˙baa+rbaa×ωbaa),ω^bab=ωbab+εvbab=ωbab+ε(r˙bab+ωbab×rbab).ωbaa, ωbab are attitude velocity of the spacecraft 2 relative to spacecraft 1 expressed in Oa-xayaza and Ob-xbybzb. rbaa, rbab are relative position vectors of the center mass of the spacecraft expressed in Oa-xayaza and Ob-xbybzb.ω
^
b
a
b is relative velocity motor of two spacecraft expressed in Ob-xbybzb, it can be expressed as
(4)ω^bab=ω^bb-q^ba*∘ω^aa∘q^ba,
where ω^bb, ω^aa are velocity motors of two spacecraft expressed in their own body fixed frames, q^ba is spiral motion of Ob-xbybzb with respect to Oa-xayaza, and q^ba* is conjugation of q^ba.


### 2.2. Dynamic Equation Based on Dual Quaternion

For rigid-body spacecraft, the dual inertia operator is defined as(5)M^=m^E+I^=mddεE+εI=[mddε+εIxxεIxyεIxzεIyxmddε+εIyyεIyzεIzxεIzymddε+εIzz],
where m is the mass of spacecraft, I is the inertia matrix, and E is a 3 × 3 identity matrix:
(6)M^-1=ddεI-1+ε1mE.Based on single spacecraft dynamic equationF^c=(d/dt)H^c=M^(d/dt)ω^c+ω^c×M^ω^c=M^ω^˙cc+ωc×M^ω^c, differentiating (4), and combining with the single spacecraft kinematic equation expressed by dual quaternion q^˙=(1/2)ω^bi∘q^=(1/2)q^∘ω^bb and q^˙*=-(1/2)ω^bb∘q^*, assume that [ω^bab×]=(1/2)([ω^bab+]-[ω^bab–]), the equation can be written as
(7)ω^˙bab=ω^˙bb-q^ba*∘ω^˙aa∘q^ba-q^˙ba*∘ω^aa∘q^ba-q^ba*∘ω^aa∘q^˙ba=ω^˙bb-q^ba*∘ω^˙aa∘q^ba-q^ba*∘ω^aa∘12q^ba∘ω^bab+12ω^bab∘q^ba*∘ω^aa∘q^ba=ω^˙bb-q^ba*∘ω^˙aa∘q^ba+12([ω^bab+]-[ω^bab–])q^ba*∘ω^aa∘q^ba.
The following equation is derived from the single spacecraft dynamic equation:
(8)ω^˙bab=-M^b-1(ω^bb×M^bω^bb)+M^b-1F^bb.
Substituting (8) into (7) yields
(9)ω^˙bab=-M^b-1(ω^bb×M^bω^bb)+M^b-1F^bb-q^ba*∘ω^˙aa∘q^ba+[ω^bab×]q^ba*∘ω^aa∘q^ba.
Equation (8) is relative dynamic equation based on dual quaternion.ω
^
b
b is obtained through (4) and expressed as
(10)ω^bb=ω^bab+q^ba*∘ω^aa∘q^ba.ω^˙bab is the derivative of ω^bab in Ob-xbybzb, M^b is dual inertia operator of spacecraft 2, and F^bb is dual force acting on the spacecraft 2 at the center of mass.According to the demands, we do simplified assumptions: assume that the orbit of spacecraft 1 is circular orbit, and spacecraft 1 without maneuver; thusω^aa=[00-ω]T+ε[ωra00]T.The following equation is obtained:(11)ω^˙aa=0,
where ω is orbit angular velocity of spacecraft 1 and ra is orbit radius of spacecraft 1. Relative dynamic equation can be simplified as
(12)ω^˙bab=-M^b-1(ω^bb×M^bω^bb)+M^b-1F^bb+[ω^bab×]q^ba*∘ω^aa∘q^ba.


## 2.1. Kinematic Equation Based on Dual Quaternion

Figure1 shows two rigid-body spacecraft orbiting the earth: one is spacecraft 1 and the other is spacecraft 2. Their body fixed frames, respectively, are Oa-xayaza and Ob-xbybzb. We can see the relative motion of the two spacecraft as the spiral motion of the Ob-xbybzb relative to the Oa-xayaza. The relative motion is described by dual quaternion; the expression is
(1)q^ba=qba+ε12rbaa∘qba=qba+ε12qba∘rbab,
where qba is attitude quaternion of the spacecraft 1 relative to spacecraft 2 and rba=rb-ra is relative position vector of two-spacecraft center mass.Figure 1
The orbit of the relative motion.The kinematic equation for the relative motion of two spacecraft expressed by dual quaternion is governed by(2)2q^˙ba=ω^baa∘q^ba=q^ba∘ω^bab,
where
(3)ω^baa=ωbaa+εvbaa=ωbaa+ε(r˙baa+rbaa×ωbaa),ω^bab=ωbab+εvbab=ωbab+ε(r˙bab+ωbab×rbab).ωbaa, ωbab are attitude velocity of the spacecraft 2 relative to spacecraft 1 expressed in Oa-xayaza and Ob-xbybzb. rbaa, rbab are relative position vectors of the center mass of the spacecraft expressed in Oa-xayaza and Ob-xbybzb.ω
^
b
a
b is relative velocity motor of two spacecraft expressed in Ob-xbybzb, it can be expressed as
(4)ω^bab=ω^bb-q^ba*∘ω^aa∘q^ba,
where ω^bb, ω^aa are velocity motors of two spacecraft expressed in their own body fixed frames, q^ba is spiral motion of Ob-xbybzb with respect to Oa-xayaza, and q^ba* is conjugation of q^ba.


## 2.2. Dynamic Equation Based on Dual Quaternion

For rigid-body spacecraft, the dual inertia operator is defined as(5)M^=m^E+I^=mddεE+εI=[mddε+εIxxεIxyεIxzεIyxmddε+εIyyεIyzεIzxεIzymddε+εIzz],
where m is the mass of spacecraft, I is the inertia matrix, and E is a 3 × 3 identity matrix:
(6)M^-1=ddεI-1+ε1mE.Based on single spacecraft dynamic equationF^c=(d/dt)H^c=M^(d/dt)ω^c+ω^c×M^ω^c=M^ω^˙cc+ωc×M^ω^c, differentiating (4), and combining with the single spacecraft kinematic equation expressed by dual quaternion q^˙=(1/2)ω^bi∘q^=(1/2)q^∘ω^bb and q^˙*=-(1/2)ω^bb∘q^*, assume that [ω^bab×]=(1/2)([ω^bab+]-[ω^bab–]), the equation can be written as
(7)ω^˙bab=ω^˙bb-q^ba*∘ω^˙aa∘q^ba-q^˙ba*∘ω^aa∘q^ba-q^ba*∘ω^aa∘q^˙ba=ω^˙bb-q^ba*∘ω^˙aa∘q^ba-q^ba*∘ω^aa∘12q^ba∘ω^bab+12ω^bab∘q^ba*∘ω^aa∘q^ba=ω^˙bb-q^ba*∘ω^˙aa∘q^ba+12([ω^bab+]-[ω^bab–])q^ba*∘ω^aa∘q^ba.
The following equation is derived from the single spacecraft dynamic equation:
(8)ω^˙bab=-M^b-1(ω^bb×M^bω^bb)+M^b-1F^bb.
Substituting (8) into (7) yields
(9)ω^˙bab=-M^b-1(ω^bb×M^bω^bb)+M^b-1F^bb-q^ba*∘ω^˙aa∘q^ba+[ω^bab×]q^ba*∘ω^aa∘q^ba.
Equation (8) is relative dynamic equation based on dual quaternion.ω
^
b
b is obtained through (4) and expressed as
(10)ω^bb=ω^bab+q^ba*∘ω^aa∘q^ba.ω^˙bab is the derivative of ω^bab in Ob-xbybzb, M^b is dual inertia operator of spacecraft 2, and F^bb is dual force acting on the spacecraft 2 at the center of mass.According to the demands, we do simplified assumptions: assume that the orbit of spacecraft 1 is circular orbit, and spacecraft 1 without maneuver; thusω^aa=[00-ω]T+ε[ωra00]T.The following equation is obtained:(11)ω^˙aa=0,
where ω is orbit angular velocity of spacecraft 1 and ra is orbit radius of spacecraft 1. Relative dynamic equation can be simplified as
(12)ω^˙bab=-M^b-1(ω^bb×M^bω^bb)+M^b-1F^bb+[ω^bab×]q^ba*∘ω^aa∘q^ba.


## 3. Gyroscope and Star Sensor Measurement Model

### 3.1. Gyroscope Measurement Model

Gyroscope is one of attitude measurement sensors which are the most widely used in modern navigation control areas, generally output data in the form of angular velocity or angular velocity increment with respect to inertial space. Combined with other sensors, such as accelerometers, it can navigate and posit the carrier.Generally, directing output of the gyroscope is projected in the gyroscope coordinate system of angular velocity in inertial system. For simplicity, assuming gyroscope measurement coordinate system coincides with the body coordinate system of spacecraft. Gyroscope mathematical model can be represented by this:(13)ω~(t)=ω(t)+β(t)+ηv(t),β˙(t)=ηu(t).ω
~ is continuous measurement output value of the gyroscope, β(t) is gyroscope drift, and ηv(t) and ηu(t) are Gaussian white noise process with uncorrelated zero mean, and they conformed:
(14)E{ηv(t)ηvT(t)}=I3×3σv2δ(t-τ),E{ηu(t)ηuT(t)}=I3×3σu2δ(t-τ).


### 3.2. Star Sensor Measurement Model

Star sensor is currently of the highest accuracy, the most widely used attitude sensor; stellar is its reference source of attitude measurement, and it can output vector direction of stars in star tracker coordinates. It provides high-accuracy measurements of spacecraft attitude control and navigation.The output of star sensor can be a reference axis vector and also be Euler angles or quaternions. Euler angles obtained through the star sensor data processing unit and the star sensor output is:(15)ψm=ψ+vψ,θm=θ+vθ,φm=φ+vφ.ψ
m, θm, φm are measurements of star sensor, ψ, θ, and φ are true attitude angles, and vs=[vψvθvφ]T is star sensor measurement noise.


## 3.1. Gyroscope Measurement Model

Gyroscope is one of attitude measurement sensors which are the most widely used in modern navigation control areas, generally output data in the form of angular velocity or angular velocity increment with respect to inertial space. Combined with other sensors, such as accelerometers, it can navigate and posit the carrier.Generally, directing output of the gyroscope is projected in the gyroscope coordinate system of angular velocity in inertial system. For simplicity, assuming gyroscope measurement coordinate system coincides with the body coordinate system of spacecraft. Gyroscope mathematical model can be represented by this:(13)ω~(t)=ω(t)+β(t)+ηv(t),β˙(t)=ηu(t).ω
~ is continuous measurement output value of the gyroscope, β(t) is gyroscope drift, and ηv(t) and ηu(t) are Gaussian white noise process with uncorrelated zero mean, and they conformed:
(14)E{ηv(t)ηvT(t)}=I3×3σv2δ(t-τ),E{ηu(t)ηuT(t)}=I3×3σu2δ(t-τ).


## 3.2. Star Sensor Measurement Model

Star sensor is currently of the highest accuracy, the most widely used attitude sensor; stellar is its reference source of attitude measurement, and it can output vector direction of stars in star tracker coordinates. It provides high-accuracy measurements of spacecraft attitude control and navigation.The output of star sensor can be a reference axis vector and also be Euler angles or quaternions. Euler angles obtained through the star sensor data processing unit and the star sensor output is:(15)ψm=ψ+vψ,θm=θ+vθ,φm=φ+vφ.ψ
m, θm, φm are measurements of star sensor, ψ, θ, and φ are true attitude angles, and vs=[vψvθvφ]T is star sensor measurement noise.


## 4. Relative Attitude Determination Algorithm Based on Extended Kalman Filter (EKF)

### 4.1. State Equation and Linearization

State variableX(t)=[q^baω^ba]T. Nonlinear continuity equation is expressed as
(16)X˙(t)=f(X(t),t)+w(t),
where(17)f(X(t),t)=[q^baω^ba]T=[12q^ba∘ω^ba-M^b-1(ω^bb×M^bω^bb)+Mb-1F^bb-q^ba*∘(ω^aa)∘q^ba+[ω^bab×]q^ba*∘(ω^aa)∘q^ba].Make the formula (16) to be linearization and discretization; we can get
(18)Xk=Φk,k-1Xk,k-1+Wk-1,
where Φk,k-1=I+Ck-1T, T is sampling time:
(19)Ck-1=∂f∂XT|X=Xk-1=[∂q^˙ba∂q^baT∂q^˙ba∂ω^baT∂ω^˙ba∂q^baT∂ω^˙ba∂ω^baT].


### 4.2. Measurement Equation and Linearization

The measurement equation is as follows:(20)Z(t)=h(X(t),t)+v(t),
where v(t) is measurement noise. This is a nonlinear equation of variable X. Making it to be linearization and discretization, so
(21)Zk=HkXk+Vk.


### 4.3. EKF Filtering Algorithm Process

State single-step forecast:(22)x^k,k-1=x^k-1+f(x^k-1)T.
Single-step forecast error of mean square:(23)Pk,k-1=Φk,k-1Pk-1Φk,k-1T+Qk-1,
whereΦk,k-1 is single-step transfer matrix of system.
Filtration transmission gain:(24)Kk=Pk,k-1HkT[HkPk·k-1HkT+Rk]-1.
H
k is measurement matrix and Rk is measurement noise matrix.
Status estimated value:(25)x^k=x^k,k-1+Kk[zk-h(x^k,k-1)].
Filtering error of mean square:(26)Pk=[I-KkHk]Pk,k-1.


## 4.1. State Equation and Linearization

State variableX(t)=[q^baω^ba]T. Nonlinear continuity equation is expressed as
(16)X˙(t)=f(X(t),t)+w(t),
where(17)f(X(t),t)=[q^baω^ba]T=[12q^ba∘ω^ba-M^b-1(ω^bb×M^bω^bb)+Mb-1F^bb-q^ba*∘(ω^aa)∘q^ba+[ω^bab×]q^ba*∘(ω^aa)∘q^ba].Make the formula (16) to be linearization and discretization; we can get
(18)Xk=Φk,k-1Xk,k-1+Wk-1,
where Φk,k-1=I+Ck-1T, T is sampling time:
(19)Ck-1=∂f∂XT|X=Xk-1=[∂q^˙ba∂q^baT∂q^˙ba∂ω^baT∂ω^˙ba∂q^baT∂ω^˙ba∂ω^baT].


## 4.2. Measurement Equation and Linearization

The measurement equation is as follows:(20)Z(t)=h(X(t),t)+v(t),
where v(t) is measurement noise. This is a nonlinear equation of variable X. Making it to be linearization and discretization, so
(21)Zk=HkXk+Vk.


## 4.3. EKF Filtering Algorithm Process

State single-step forecast:(22)x^k,k-1=x^k-1+f(x^k-1)T.
Single-step forecast error of mean square:(23)Pk,k-1=Φk,k-1Pk-1Φk,k-1T+Qk-1,
whereΦk,k-1 is single-step transfer matrix of system.
Filtration transmission gain:(24)Kk=Pk,k-1HkT[HkPk·k-1HkT+Rk]-1.
H
k is measurement matrix and Rk is measurement noise matrix.
Status estimated value:(25)x^k=x^k,k-1+Kk[zk-h(x^k,k-1)].
Filtering error of mean square:(26)Pk=[I-KkHk]Pk,k-1.


## 5. Simulation and Result Analysis

Under the background of the final collaborative period of two-spacecraft earth observing (relative distance is 1000-0 meters), we carry on the simulation based on kinematics and dynamics equations. Simulation verifies the two spacecraft’s kinematics and dynamics equations, exerts perturbative force and moment to spacecraft 2, and verifies that the relative velocity and angular velocity of the spacecraft change along with the variation of force and moment. Integrating the two spacecraft’s kinematics and dynamics equation, integration step is 0.1 seconds.Suppose that the spacecraft 1 moves in circular orbit whose height is 500 kilometers. Because spacecraft 1 is directional to earth, its body coordinate system coincides with orbital coordinate system.The three axes of the spacecraft 2’s body coordinate system coincide with its inertia principal axis, its mass is 1000 kg, and its moment of inertia is(27)J=[722000820000876]kg·m2.
The attitude and position of spacecraft 2 relative to the spacecraft 1 at initial time are(28)qba(0)=[0.3772-0.4329-0.66450.4783]T,rbab(0)=[-200-100100]Tm,ωbab(0)=[0.000960.001400.00105]Trad/s,r˙bab(0)=[0.250.3820.074]Tm/s.
The attitude and position of the goal are(29)qba(t)=[1000]T,rbab(t)=[000]Tm,ωbab(t)=[000]Trad/s,r˙bab(t)=[000]Tm/s.The measurement noise mean square error of star sensorσs=1/3600, initial gyro constant drift b0=[555]T deg/h, mean square error of measurement noise σg1=0.005 deg/h, and mean square error of slope white noise σg1=0.1 deg/h.Relative motion of two spacecraft changes in coordinate systemOb-xbybzb is as shown.Figures2~3, respectively, are perturbative force and moment to spacecraft 2Figure 2
Control force.Figure 3
Control moment.Figures4, 5, 6, 7, 8, and 9, respectively, are attitude quaternion, attitude angular, attitude angular velocity, relative position, and relative velocity in the case of dynamics, star sensor, and filtering. They are parameters estimation of EKF algorithm.Figure 4
Attitude quaternion.Figure 5
Quaternion.Figure 6
Attitude angular.Figure 7
Attitude angular velocity.Figure 8
Relative position.Figure 9
Relative velocity.


## 6. Conclusion

Based on Extended Kalman Filter, combined with the spacecraft relative dynamics equation which is described by dual quaternion, this paper designs the relative attitude algorithm of gyroscope/star sensor combined measurement. The simulation results show that the EKF algorithm has faster convergence speed and higher accuracy.



---

*Source: 602724-2014-05-08.xml*