# Unscented Particle Smoother and Its Application to Transfer Alignment of Airborne Distributed POS

Authors: Xiaolin Gong; Xiaorui Zheng; Xing-Gang Yan; Zhaoxing Lu
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

This paper deals with the problem of state estimation for the transfer alignment of airborne distributed position and orientation system (distributed POS). For a nonlinear system, especially with large initial attitude errors, the performance of linear estimation methods will degrade. In this paper, a nonlinear smoothing algorithm called the unscented particle smoother (UPS) is proposed and utilized in the off-line transfer alignment of airborne distributed POS. In this algorithm, the measurements are first processed by the forward unscented particle filter (UPF) and then a backward smoother is used to achieve the improved solution. The performance of this algorithm is compared with that of a similar smoother known as the unscented Rauch-Tung-Striebel smoother. The simulation results show that the UPS effectively improves the estimation accuracy and this work offers a new off-line transfer alignment approach of distributed POS for multiantenna synthetic aperture radar and other airborne earth observation tasks.


---

## Body

## 1. Introduction

Airborne synthetic aperture radar (SAR) is an important tool for earth observation. Multiantenna SAR can achieve higher accuracy of interferometric processing and ground moving target indication through multibaseline interferometry technology (i.e., interferometric SAR), and it has become an important research direction of radar remote sensing [1]. Usually, these antennae are installed on both sides of the wing of the aircraft, and the motion parameters of each antenna are needed to complete the motion compensation [2]. As a result, airborne distributed position and orientation system (distributed POS) is developed, which is used to provide the time and space reference information for all points of multiantenna SAR or other requirements [3].Airborne distributed POS is generally composed of four parts: a high-precision POS (main system), a few low-precision inertial measurement units (IMU, subsystem), a POS computer system (PCS), and postprocessing software [4, 5]. The high precision main system is usually installed in the cabin or underpart of the aircraft and performs the inertial/satellite integrated estimation with the position and velocity of global navigation satellite system (GNSS) as the measurement information. The subsystem, which is placed as close as possible to the measurement center of the load, is used to measure the acceleration and angular velocity which are sent to the PCS. In the PCS, the high precision motion information of main system is used as the reference signal to achieve the motion information of subsystem by transfer alignment technology. The nature of transfer alignment is to estimate and compensate the calculation error of motion parameters of the subsystem by using the high precision information of the main system [6]. It shows that the estimation method is crucial for distributed POS.It should be noted that distributed POS is a nonlinear system. In some emergency situations, the distributed POS is expected to try to shorten the preparation time on the ground and even start up in flight, which will bring the uncertain errors and the misalignment of initial attitude is not small. These make the nonlinearity of the system to be further increased. Thus, the performance of transfer alignment based on the linear model and linear estimation method will degrade [7].The unscented Kalman filter (UKF) is a typical nonlinear filtering method. It uses a deterministic sampling approach named as sigma points to propagate nonlinear systems and has been discussed in many literatures for the inertial/satellite integration system [8–10]. UKF is a representative filter to overcome the flaw of truncated error in extended Kalman filter (EKF) [11]. However, the UKF method is an approximation of linear minimum variance estimation and does not apply to general non-Gaussian distributions [12, 13]. Another popular solution strategy for general filtering problems is particle filter (PF) which is based on the recursive Bayesian estimation and suitable for the nonlinear and non-Gaussian distributed system [14]. The estimated accuracy of PF is higher than the accuracy of using EKF or UKF alone [7]. PF uses a set of sampling points extracted from the posterior probability to express its distribution. Choosing a reasonably recommended probability density is the core of PF. The closer the recommended density selected is to the true density, the better the filter effect is. Otherwise, it is worse or even divergent [7, 15]. The UKF is able to more accurately propagate the mean and covariance of the Gaussian approximation to the state distribution than the EKF. So PF is combined with UKF and formed a new PF called unscented particle filter (UPF) [15], in which the UKF is used for proposal distribution generation and the performance of the filtering algorithm is improved by taking into account the influence of the measured values on the state estimation.The above UKF and PF are based on the idea of forward filtering, where only the observation information of the current moment and the previous time are used. The other methods, such as smoothing estimation, can make full use of all the observation information to estimate the state of each moment, and is particularly suitable for off-line data processing. Since the smoother uses more observations than the filter, although it cannot be used in real time, the accuracy of the optimal smoothing algorithm is theoretically higher than that of the KF [16]. The unscented Rauch-Tung-Striebel smoother (URTSS) was first proposed by Särkkä in 2008 [17]. This smoother takes the benefit over unscented transformation as well as the smoothing iterative characteristics, which makes the estimation accuracy much better than that of the UKF. In ref. [18], the URTSS is applied to the off-line integrated estimation of the inertial/satellite integration system and shows an obvious accuracy advantage over extended R-T-S smoother in attitude estimation. In particle smoothing, ref. [19] proposed a particle smoothing algorithm which is a direct extension of the PF. This method has the advantage of fixed interval smoothing, but it cannot express the state variables in the past very well. Then, the forward-backward smoother and two-filter smoother based on PF are proposed [20]. This kind of smoothing algorithm needs to carry on the particle filter first and store the particles with weight selected in the approximate process and then perform the smoothing process. However, there is still the problem of particle degeneracy.For the airborne earth observation system, there are two modes of image processing: real time and off-line. In the off-line cases, there is no high requirement for fast data processing in real-time imaging mode. The smoother can be employed to provide a better solution for off-line image processing. Motivated by these, a smoother called unscented particle smoother (UPS) is proposed in this paper, which is then applied to the off-line transfer alignment of airborne distributed POS. The UPS is made by taking the advantages of the UKF, PF, and smoothing. In the UPS, UKF is used to generate the proposed distributions to overcome the shortage of particle degradation in PF. Then, the forward-backward smoother using all observations is further combined together to obtain higher estimation accuracy. The UPS can be applied to nonlinear and non-Gaussian noise systems as well. There is no need to keep still to obtain the initial motion information of the main system and subsystem before flying off.The rest of this paper is organized as follows. In Section2, the UPS are proposed. In Section 3, the design of UPS for airborne distributed POS is given in detail. In Section 4, the performance of the UPS method is demonstrated and compared with that of the URTSS in estimation accuracy by a simulation test. Finally, conclusions are given in Section 5.


## 2. Unscented Particle Smoothing Algorithm

The proposed UPS has two structures: a forward filter and a backward smoother. The forward filter is UPF used to obtain the filtering state estimation. In UPF, the proposed density is determined by UKF, which not only solves the problem of the degradation of particles but also enables particles to get the latest a posteriori information of the measurement vector when they are updated, which is helpful for particles to move toward the area with high likelihood. Then, the backward smoother is conducted after the forward filter to modify the importance weights to achieve the smoothed state estimation.Suppose the discrete form of thenx-state system equation and the observation equation is as follows:
(1)xk=fxk−1,k−1+Γk−1wk−1,yk=hxk,k+vk,where xk∈Rnx represents the state vector of the system at time k and yk∈Rny denotes the measurement vector at time tk. The functions f and h describe state and measurement models, respectively [8, 21]; Γk−1 is the state noise distribution matrix, and wk−1 and vk represent the state and measurement noises, respectively, which can be described by the corresponding probability density functions pwk−1 and pvk. In this work, wk−1 and vk are assumed to the uncorrelated white Gaussian noise. The mean value of wk−1 and vk is zero and their variance matrices are Qk−1 and Rk, respectively [8, 22].Based on (1), the process of UPS is described as follows.Step 1.
Initialization (k=0) [7, 15]. Suppose the initial state variable x0~px0, and the covariance matrix is P0.
Fori=1,2,…,N, draw the states (particles) x0i from the prior probability density px0 and set x0i~Nx^0i,P0i, where x^0i=x0i and P0i=P0. N is the number of particles.Step 2.
Importance sampling step (fori=1,2,…,N, k=1,2,3,…,T).
(1)
Update the particles with the UKF:(i)
Calculate sigma points(2)χk−10i=x^k−1i,χk−1ji=x^k−1i+nx+λPk−1ijj=1,2,…,nx,χk−1ji=x^k−1i−nx+λPk−1ijj=nx+1,nx+2,…,2nx,where the sigma point set S=χk−1ji,Wjl′;j=0,1,…,2nx,l′∈m,c is composed of the sigma points χk−1ji and their respective mean (m) and covariance (c) weights Wjl′. λ=α2nx+κ−nx denotes a scaling parameter. The parameter α determines the spread of samples that is usually set between 10−4 and 1. The constant κ is typically set to 0 or 3−nx where nx is the state dimension. The corresponding weights are given by:
(3)W0m=λnx+λ,W0c=λnx+λ+1−α2+β,Wjm=Wjc=0.5nx+λj=1,2,…,2nx,where the parameter β is used to incorporate prior knowledge of the distribution (for Gaussian distribution, β=2 is optimal).
(ii)
Time update (propagate particle into future)(4)χk∣k−1ji=fχk−1ji,k−1,x^k∣k−1i=∑j=02nxWjmχk∣k−1ji,Pk∣k−1i=∑j=02nxWjcχk∣k−1ji−x^ik∣k−1χk∣k−1ji−x^ik∣k−1T+Γk−1QkΓk−1T,γk∣k−1ji=hχk∣k−1ji,k,y^k∣k−1i=∑j=02nxWjmγk∣k−1ji,where Pk∣k−1i denotes the prior error covariance.
(iii)
Measurement update (incorporate new observation)(5)Py^ky^ki=∑j=02nxWjcγk∣k−1ji−y^k∣k−1iγk∣k−1ji−y^k∣k−1iT+Rk,Px^ky^ki=∑j=02nxWjcχk∣k−1ji−x^k∣k−1iγk∣k−1ji−y^k∣k−1iT,Kki=Px^ky^kiPy^ky^ki−1,x^ki=x^k∣k−1i+Kkiyk−y^k∣k−1i,Pki=Pk∣k−1i−KkiPy^ky^kKkiT,where y^k∣k−1i denotes the predicted measurement, Py^ky^ki denotes the covariance of predicted measurement error, and Px^ky^ki denotes the cross-covariance matrix between x^k∣k−1i and y^k∣k−1i. Kki denotes the filter gain.
(2)
Importance sampling step.
Fori=1,2,…,N, sample x~ki~qxki∣x0:k−1i,y1:k=Nx^ki,Pki and set x~ki≜x0:k−1i,x~ki and P~0:ki≜P0:k−1i,P~ki.
(3)
Update and normalize the importance weights.
Fori=1,2,…,N, evaluate the importance weights up to a normalizing constant:
(6)w~ki=wk−1ipyk∣x~kipx~ki∣xk−1iqx~ki∣x0:k−1i,y1:k.
Fori=1,2,…,N, normalize and store the importance weights:
(7)w~ki=wki∑i=1Nwki−1.
Compute the filtering state estimationx^kF:
(8)x^kF=∑i=1Nw~kix^ki.Step 3.
Using the importance weights stored in the forward filtering to conduct the backward smoothing recursion.
SetwT∣Ti=w¯T∣Ti. For k=T−1,…,1 and i=1,2,…,N, evaluate
(9)wk∣Ti=wk+1∣Ti∑j=1Nwk+1∣Tjpx~k+1j∣x~ki∑m=1Nw~kmpx~k+1j∣x~km.
(1)
Resampling
Forx~ki (i=1,2,…,N), Using sequential importance resampling (SIR) [15] to complete resampling to obtain N random samples x′0:ki, and set wki=w~ki=1/N.
(2)
Compute the smoothed state estimationx^kS
(10)x^kS=1N∑i=1Nx′ki.The principle diagram of UPS is shown in Figure1. In summary, there are two processes: a forward filtering process and a backward smoothing process. The measurements are first processed by the forward filter and the importance weight w~ki should be stored in the forward filter stage. And then, a separate backward smoothing pass is used to modify the importance weights for obtaining the smoothing solution.Figure 1
The principle diagram of UPS.


## 3. UPS Design for Airborne Distributed POS

The design of UPS for airborne distributed POS is given in this section. Firstly, the definition of coordinate frames is introduced.i and e denote the earth-center inertial frame and the earth-centered earth-fixed frame, respectively. The navigation frames of the main system and subsystem are both topocentric frames (the z-axis parallel to the upward vertical, the x-axis pointing eastward, and the y-axis pointing northward), represented by n and ns, respectively. The body frames of the main system and subsystem are denoted by b and bs. A detailed description of these coordinate frames is available in [18, 23].Figure2 show the block diagram of transfer alignment in distributed POS based on UPS.Figure 2
The block diagram of transfer alignment in distributed POS based on UPS.In Figure2, the strapdown navigator calculates the attitude (heading, pitch, and roll), velocity, and position expressed in ns-frame by using the angular rates and accelerations of subsystem. Since the strapdown navigation is an integration process, any errors of the accelerometers and gyroscopes will integrate into slowly growing errors on velocity, position, and attitude. The differences of velocity and attitude between the main system and the strapdown navigator are used in the forward filtering to estimate the growing velocity, position, and attitude errors in the strapdown navigator. The closed-loop error controller is used to correct the result of strapdown navigator by using the estimation of forward filtering. After the end of the forward filter, the backward smoother is performed to compensate the estimation of forward filtering, and the improved estimation results are used to correct the result of strapdown navigator in the feed forward error controller. The improved attitude, velocity, and position are outputted to the imaging sensor at last.


### 3.1. Nonlinear Inertial Navigation Error Model

The inertial navigation error equations are the basis of the mathematical model for transfer alignment. According to the definitions above, the nonlinear SINS error model of the subsystem based on angle error, which includes attitude error equation, velocity error equation, position error equation, and inertial senor constant error equation, is given by [24]:
(11)ϕ̇ns=I−Cnspωinsns+δωinsns−Cbsnsεbs,δV̇ns=I−Cpnsfns−2δωiens+δωensns×Vns−2δωiens+δωensns×δVns+C^bns∇bs,δL̇=δVNRM+H−VNRM+H2δH,δλ̇=secLRN+HδVE+VEsecLtanLRN+HδL−VEsecLRN+H2δH,δḢ=δVU,ε̇cbs=0,∇̇cbs=0,where ϕns=ϕEϕNϕUT denotes the attitude error vector in ns-frame, the subscripts E, N, and U denote the east, north, and up components, respectively; ωinsns denotes the rotation velocity of the ns-frame relative to the i-frame expressed in ns-frame with error δωinsns; Cbsns denotes the coordinate transformation matrix from the bs-frame to the ns-frame; εbs denotes the gyro drift of subsystem in bs-frame, which consists of random constant drift εcbs and the Gaussian white noise ωεbs with εcbs=εxbsεybsεzbsT and ωεbs=ωεxbsωεybsωεzbsT; Vns=VEVNVUT denotes the velocity vector in ns-frame with error δVns=δVEδVNδVUT; fns=fEfNfUT denotes the specific force measured by the accelerometers of subsystem expressed in ns-frame; ωiens denotes the rotation velocity of the e-frame relative to the i-frame expressed in ns-frame with error δωiens; ωensns denotes the rotation velocity of the ns-frame relative to the e-frame expressed in ns-frame with error δωensns; ∇bs denotes the accelerometer bias of subsystem in bs-frame, which consists of random constant bias ∇cbs and the Gaussian white noise ω∇bs with ∇cbs=∇xbs∇ybs∇zbsT and ω∇bs=ω∇xbsω∇ybsω∇zbsT; secL = 1/cosL. The symbols δL, δλ, and δH denote the latitude error, longitude error, and height error, respectively; and RM and RN denote the meridian and transverse radius of curvature, respectively.The nonlinear terms of attitude and velocity error equations areI−Cnspωinsns and I−Cpnsfns in (11), respectively, and the nonlinear matrix Cnsp=CpnsT is expressed as follows:
(12)Cnsp=cosϕNcosϕU−sinϕEsinϕNsinϕUcosϕNsinϕU+sinϕEsinϕNcosϕU−cosϕEsinϕN−cosϕEsinϕUcosϕEcosϕUsinϕEsinϕNcosϕU+sinϕEcosϕNsinϕUsinϕNsinϕU−sinϕEcosϕNcosϕUcosϕEcosϕN.


### 3.2. Model of Rigid Misalignment Angle and Flexure Angle

During the transfer alignment, there are flexure angles and rigid misalignment angles between the main system and the subsystem. These angles cannot be measured accurately and the flexure angles vary with the time. Therefore, the flexure angles and rigid misalignment angles should be modeled and estimated.The models of rigid misalignment angles and flexure angles are shown in the following equations:(13)ρ̇j=0,θ¨j+2βjθ̇j+βj2θj=ηjj=x,y,z,where ρj denotes the rigid misalignment angle of subsystem; θj denotes the flexure angle and described by the second-order Markov process [25]; βj=2.146/τj and τj is the correlation time; x, y, and z are the axes of bs-frame; ηj denotes the Gaussian white noise with covariance Qηj=4βj3σj2; and σj2 denotes the covariance of θj; Qη=QηxQηyQηzT.


### 3.3. Design of the Smoother for Transfer Alignment of Distributed POS

#### 3.3.1. State Equation

On the basis of models established in Subsections3.1 and 3.2, let
(14)x=x1x2T,x1=ϕEϕNϕUδVEδVNδVUδLδλδHεxbsεybsεzbs∇xbs∇ybs∇zbsT,x2=ρxρyρzθxθyθzθ̇xθ̇yθ̇zT.Then, the continuous-time system model can be given by(15)ẋ=fx,t+GtWIt,where the nonlinear function fx,t is determined by error model of transfer alignment in Subsections 3.1 and 3.2, the vector WIt represents the system noise. W1=ωεxbsωεybsωεzbsω∇xbsω∇ybsω∇zbsηxηyηzT is the zero-mean Gaussian white noise vector with covariance QI which consists of covariance Qε of gyro random drift, covariance Q∇, of accelerometer random bias, and Qη.


#### 3.3.2. Measurement Equation

The differences of velocity and attitude between the main system and the subsystem solutions are considered as the measurement vector. The measurement vector can be given as(16)yt=δψδθδγδVE′δVN′δVU′T,where δψ, δθ, and δγ denote the differences of heading angle, pitch angle, and roll angle between the main system and the subsystem; δVE′, δVN′, and δVU′ denote the differences of the east, north, and up velocities between the main system and the subsystem accounting for the lever arm between the main system and subsystem [25].The measurement vector can be expressed linearly by the state vector shown in (14). Hence, the linear measurement model can be given by
(17)yt=Hxt+vt,where the measurement noise vt=vδψvδθvδγvδVE′vδVN′vδVU′T is the white noise vector of the solution of the main system with zero mean and covariance matrix R. The measurement matrix Ht can be given as
(18)Ht=H103×303×9H2H303×303×3I3×303×903×303×303×36×24.Here, the attitude matrix of main systemCbn is denoted as Ta. Ta, H1, H2, and H3 are as follows:
(19)Cbn=cosγmcosφm−sinγmsinθmsinφm−cosθmsinφmsinγmcosφm+cosγmsinθmsinφm−cosγmsinφm+sinγmsinθmcosφmcosθmcosφmsinγmsinφm−cosγmsinθmcosφm−sinγmcosθmsinθmcosγmcosθm,where φm, θm, and γm denote the heading, pitch, and roll of the main system, respectively.
(20)H1=Ta12Ta32Ta122+Ta222Ta22Ta32Ta122+Ta222−1−Ta221−Ta322Ta121−Ta3220Ta21Ta33−Ta31Ta23Ta332+Ta312Ta31Ta13−Ta11Ta33Ta332+Ta3120,H2=H3=Ta12Ta23−Ta13Ta22Ta122+Ta2220Ta11Ta22−Ta12Ta21Ta122+Ta222Ta331−Ta3220−Ta311−Ta322−Ta31Ta32Ta332+Ta3121−Ta32Ta33Ta332+Ta312.The UPS for transfer alignment of airborne distributed POS postprocessing can be summarized into two parts: the forward filtering solution and the backward smoothing solution.Figure3 shows the data flaw of UPS used in airborne distributed POS and the relation between the state vector estimation of forward filtering x^kF and backward recursion x^kS at time k and k+1. In forward filtering solution, the filter results are applied to the nominal trajectory during each filtering step and will be reset to zero after the correction of strapdown navigation (velocity, position, and attitude). The function of the backward smoothing solution is to compensate for the estimation of the forward filter stored during the forward filtering process. The filtering and smoothing are calculated at the main system rate. And after each estimation points, a strapdown inertial navigation algorithm is performed to provide the navigation information.Figure 3
The data flaw of UPS used in airborne distributed POS.


## 3.1. Nonlinear Inertial Navigation Error Model

The inertial navigation error equations are the basis of the mathematical model for transfer alignment. According to the definitions above, the nonlinear SINS error model of the subsystem based on angle error, which includes attitude error equation, velocity error equation, position error equation, and inertial senor constant error equation, is given by [24]:
(11)ϕ̇ns=I−Cnspωinsns+δωinsns−Cbsnsεbs,δV̇ns=I−Cpnsfns−2δωiens+δωensns×Vns−2δωiens+δωensns×δVns+C^bns∇bs,δL̇=δVNRM+H−VNRM+H2δH,δλ̇=secLRN+HδVE+VEsecLtanLRN+HδL−VEsecLRN+H2δH,δḢ=δVU,ε̇cbs=0,∇̇cbs=0,where ϕns=ϕEϕNϕUT denotes the attitude error vector in ns-frame, the subscripts E, N, and U denote the east, north, and up components, respectively; ωinsns denotes the rotation velocity of the ns-frame relative to the i-frame expressed in ns-frame with error δωinsns; Cbsns denotes the coordinate transformation matrix from the bs-frame to the ns-frame; εbs denotes the gyro drift of subsystem in bs-frame, which consists of random constant drift εcbs and the Gaussian white noise ωεbs with εcbs=εxbsεybsεzbsT and ωεbs=ωεxbsωεybsωεzbsT; Vns=VEVNVUT denotes the velocity vector in ns-frame with error δVns=δVEδVNδVUT; fns=fEfNfUT denotes the specific force measured by the accelerometers of subsystem expressed in ns-frame; ωiens denotes the rotation velocity of the e-frame relative to the i-frame expressed in ns-frame with error δωiens; ωensns denotes the rotation velocity of the ns-frame relative to the e-frame expressed in ns-frame with error δωensns; ∇bs denotes the accelerometer bias of subsystem in bs-frame, which consists of random constant bias ∇cbs and the Gaussian white noise ω∇bs with ∇cbs=∇xbs∇ybs∇zbsT and ω∇bs=ω∇xbsω∇ybsω∇zbsT; secL = 1/cosL. The symbols δL, δλ, and δH denote the latitude error, longitude error, and height error, respectively; and RM and RN denote the meridian and transverse radius of curvature, respectively.The nonlinear terms of attitude and velocity error equations areI−Cnspωinsns and I−Cpnsfns in (11), respectively, and the nonlinear matrix Cnsp=CpnsT is expressed as follows:
(12)Cnsp=cosϕNcosϕU−sinϕEsinϕNsinϕUcosϕNsinϕU+sinϕEsinϕNcosϕU−cosϕEsinϕN−cosϕEsinϕUcosϕEcosϕUsinϕEsinϕNcosϕU+sinϕEcosϕNsinϕUsinϕNsinϕU−sinϕEcosϕNcosϕUcosϕEcosϕN.


## 3.2. Model of Rigid Misalignment Angle and Flexure Angle

During the transfer alignment, there are flexure angles and rigid misalignment angles between the main system and the subsystem. These angles cannot be measured accurately and the flexure angles vary with the time. Therefore, the flexure angles and rigid misalignment angles should be modeled and estimated.The models of rigid misalignment angles and flexure angles are shown in the following equations:(13)ρ̇j=0,θ¨j+2βjθ̇j+βj2θj=ηjj=x,y,z,where ρj denotes the rigid misalignment angle of subsystem; θj denotes the flexure angle and described by the second-order Markov process [25]; βj=2.146/τj and τj is the correlation time; x, y, and z are the axes of bs-frame; ηj denotes the Gaussian white noise with covariance Qηj=4βj3σj2; and σj2 denotes the covariance of θj; Qη=QηxQηyQηzT.


## 3.3. Design of the Smoother for Transfer Alignment of Distributed POS

### 3.3.1. State Equation

On the basis of models established in Subsections3.1 and 3.2, let
(14)x=x1x2T,x1=ϕEϕNϕUδVEδVNδVUδLδλδHεxbsεybsεzbs∇xbs∇ybs∇zbsT,x2=ρxρyρzθxθyθzθ̇xθ̇yθ̇zT.Then, the continuous-time system model can be given by(15)ẋ=fx,t+GtWIt,where the nonlinear function fx,t is determined by error model of transfer alignment in Subsections 3.1 and 3.2, the vector WIt represents the system noise. W1=ωεxbsωεybsωεzbsω∇xbsω∇ybsω∇zbsηxηyηzT is the zero-mean Gaussian white noise vector with covariance QI which consists of covariance Qε of gyro random drift, covariance Q∇, of accelerometer random bias, and Qη.


### 3.3.2. Measurement Equation

The differences of velocity and attitude between the main system and the subsystem solutions are considered as the measurement vector. The measurement vector can be given as(16)yt=δψδθδγδVE′δVN′δVU′T,where δψ, δθ, and δγ denote the differences of heading angle, pitch angle, and roll angle between the main system and the subsystem; δVE′, δVN′, and δVU′ denote the differences of the east, north, and up velocities between the main system and the subsystem accounting for the lever arm between the main system and subsystem [25].The measurement vector can be expressed linearly by the state vector shown in (14). Hence, the linear measurement model can be given by
(17)yt=Hxt+vt,where the measurement noise vt=vδψvδθvδγvδVE′vδVN′vδVU′T is the white noise vector of the solution of the main system with zero mean and covariance matrix R. The measurement matrix Ht can be given as
(18)Ht=H103×303×9H2H303×303×3I3×303×903×303×303×36×24.Here, the attitude matrix of main systemCbn is denoted as Ta. Ta, H1, H2, and H3 are as follows:
(19)Cbn=cosγmcosφm−sinγmsinθmsinφm−cosθmsinφmsinγmcosφm+cosγmsinθmsinφm−cosγmsinφm+sinγmsinθmcosφmcosθmcosφmsinγmsinφm−cosγmsinθmcosφm−sinγmcosθmsinθmcosγmcosθm,where φm, θm, and γm denote the heading, pitch, and roll of the main system, respectively.
(20)H1=Ta12Ta32Ta122+Ta222Ta22Ta32Ta122+Ta222−1−Ta221−Ta322Ta121−Ta3220Ta21Ta33−Ta31Ta23Ta332+Ta312Ta31Ta13−Ta11Ta33Ta332+Ta3120,H2=H3=Ta12Ta23−Ta13Ta22Ta122+Ta2220Ta11Ta22−Ta12Ta21Ta122+Ta222Ta331−Ta3220−Ta311−Ta322−Ta31Ta32Ta332+Ta3121−Ta32Ta33Ta332+Ta312.The UPS for transfer alignment of airborne distributed POS postprocessing can be summarized into two parts: the forward filtering solution and the backward smoothing solution.Figure3 shows the data flaw of UPS used in airborne distributed POS and the relation between the state vector estimation of forward filtering x^kF and backward recursion x^kS at time k and k+1. In forward filtering solution, the filter results are applied to the nominal trajectory during each filtering step and will be reset to zero after the correction of strapdown navigation (velocity, position, and attitude). The function of the backward smoothing solution is to compensate for the estimation of the forward filter stored during the forward filtering process. The filtering and smoothing are calculated at the main system rate. And after each estimation points, a strapdown inertial navigation algorithm is performed to provide the navigation information.Figure 3
The data flaw of UPS used in airborne distributed POS.


## 3.3.1. State Equation

On the basis of models established in Subsections3.1 and 3.2, let
(14)x=x1x2T,x1=ϕEϕNϕUδVEδVNδVUδLδλδHεxbsεybsεzbs∇xbs∇ybs∇zbsT,x2=ρxρyρzθxθyθzθ̇xθ̇yθ̇zT.Then, the continuous-time system model can be given by(15)ẋ=fx,t+GtWIt,where the nonlinear function fx,t is determined by error model of transfer alignment in Subsections 3.1 and 3.2, the vector WIt represents the system noise. W1=ωεxbsωεybsωεzbsω∇xbsω∇ybsω∇zbsηxηyηzT is the zero-mean Gaussian white noise vector with covariance QI which consists of covariance Qε of gyro random drift, covariance Q∇, of accelerometer random bias, and Qη.


## 3.3.2. Measurement Equation

The differences of velocity and attitude between the main system and the subsystem solutions are considered as the measurement vector. The measurement vector can be given as(16)yt=δψδθδγδVE′δVN′δVU′T,where δψ, δθ, and δγ denote the differences of heading angle, pitch angle, and roll angle between the main system and the subsystem; δVE′, δVN′, and δVU′ denote the differences of the east, north, and up velocities between the main system and the subsystem accounting for the lever arm between the main system and subsystem [25].The measurement vector can be expressed linearly by the state vector shown in (14). Hence, the linear measurement model can be given by
(17)yt=Hxt+vt,where the measurement noise vt=vδψvδθvδγvδVE′vδVN′vδVU′T is the white noise vector of the solution of the main system with zero mean and covariance matrix R. The measurement matrix Ht can be given as
(18)Ht=H103×303×9H2H303×303×3I3×303×903×303×303×36×24.Here, the attitude matrix of main systemCbn is denoted as Ta. Ta, H1, H2, and H3 are as follows:
(19)Cbn=cosγmcosφm−sinγmsinθmsinφm−cosθmsinφmsinγmcosφm+cosγmsinθmsinφm−cosγmsinφm+sinγmsinθmcosφmcosθmcosφmsinγmsinφm−cosγmsinθmcosφm−sinγmcosθmsinθmcosγmcosθm,where φm, θm, and γm denote the heading, pitch, and roll of the main system, respectively.
(20)H1=Ta12Ta32Ta122+Ta222Ta22Ta32Ta122+Ta222−1−Ta221−Ta322Ta121−Ta3220Ta21Ta33−Ta31Ta23Ta332+Ta312Ta31Ta13−Ta11Ta33Ta332+Ta3120,H2=H3=Ta12Ta23−Ta13Ta22Ta122+Ta2220Ta11Ta22−Ta12Ta21Ta122+Ta222Ta331−Ta3220−Ta311−Ta322−Ta31Ta32Ta332+Ta3121−Ta32Ta33Ta332+Ta312.The UPS for transfer alignment of airborne distributed POS postprocessing can be summarized into two parts: the forward filtering solution and the backward smoothing solution.Figure3 shows the data flaw of UPS used in airborne distributed POS and the relation between the state vector estimation of forward filtering x^kF and backward recursion x^kS at time k and k+1. In forward filtering solution, the filter results are applied to the nominal trajectory during each filtering step and will be reset to zero after the correction of strapdown navigation (velocity, position, and attitude). The function of the backward smoothing solution is to compensate for the estimation of the forward filter stored during the forward filtering process. The filtering and smoothing are calculated at the main system rate. And after each estimation points, a strapdown inertial navigation algorithm is performed to provide the navigation information.Figure 3
The data flaw of UPS used in airborne distributed POS.


## 4. Simulation and Analysis

In order to analyze the performance of the proposed UPS, numerical simulation is provided and the comparison with the URTSS using the simulated data of a distributed POS based on one flight trajectory is presented as well.


### 4.1. Design of Simulation

In this subsection, a typical “S + U” shaped flight trajectory of airborne earth observation is designed. Figure4 shows the plane trajectory and trajectory parameters, respectively. In Figure 4, the AB and CD segment can be considered as the imaging segment.Figure 4
Plane trajectory with S-shaped maneuver plus U-shaped flight.Parameters of flight trajectory are listed as follows: initial latitude is 40°, longitude is 116°, and height is 500 m; initial flight velocity is 100 m/s; and initial heading angle, pitch angle, and roll angle are 40°, 0°, and 0°, respectively. Firstly, the aircraft flies 100 seconds at constant velocity; secondly, the aircraft turns 60° clockwise (100 seconds) and then turns 60° anticlockwise (100 seconds); and thirdly, the aircraft flies 400 seconds at constant velocity, finally turns 180° clockwise (100 seconds), and continually flies 400 seconds at constant velocity. Total flight time is 1300 s.


#### 4.1.1. Data Generation

A trajectory generator is used to generate the theoretical data of the scheduled flight trajectory, which includes position, velocity, attitude, and the output data of gyros and accelerometers. The real outputs of the main system are obtained by adding the corresponding measurement noise to the theoretical position, velocity, and attitude. Then, the theoretical outputs of gyros and accelerometers are converted by rigid misalignment angles and flexure angles, and the constant noise and random noise are added to be the inertial sensor outputs of the subsystem. Meanwhile, the motion parameter benchmarks of subsystem can be obtained through transforming theoretical position, velocity, and attitude by flexure angle.


#### 4.1.2. Parameters of the Main System and Subsystem

The measurement noises of the main system at heading, pitch, roll, and velocity are0.005°(1σ), 0.0025°(1σ), 0.0025°(1σ), and 0.005m/s (1σ), respectively. Both gyro constant drift and random drift of subsystem are 0.01°/h. Both accelerometer constant bias and random bias of the subsystem are 50μg.


#### 4.1.3. Other Parameters

The rigid misalignment angle of the subsystem relative to the main system isρ=0.5∘0.5∘0.5∘T, and the lever arm between the main system and subsystem is rbm=5m0.1m0.1mT. For multiantenna SAR, these antennae are fixed under both sides of the wing. So the flexure angle rotated around the vertical axis is large, while the flexure angles around the other two axes are small. Accordingly, the correlation time of the second-order Markov processes is selected as 2 s, 5 s, and 2 s, and the covariances of flexure angles are 0.01, 0.15, and 0.01, respectively. The curves of flexure angle and flexure angle rate selected in this simulation are shown in Figures 5 and 6. All the initial attitude errors are 40∘. The particle number is 100. The data update rate of the subsystem is 100 Hz and the filter frequency is 1 Hz.Figure 5
Flexure angle onx-axis, y-axis, and z-axis.Figure 6
Flexure angle rate onx-axis, y-axis, and z-axis.


### 4.2. Simulation Results and Analysis

The differences between theoretical motion parameters of subsystem and their estimates obtained from estimation, called estimate errors, are used to assess and compare the performance of UPS and URTSS.Figures7–10 show the estimate error curves of UPS and URTSS, including the estimate errors of attitude, velocity, position, and baseline. The means of root mean square error (RMSE) and standard deviation (STD) values of motion information estimate errors in imaging segments AB and CD are calculated and shown in Table 1. The calculation formulas of RMSE and STD are given as follows:
(21)RMSE=∑i=1n′ei2n′,STD=∑i=1n′ei−e¯2n′,where ei denotes the ith estimate error, n′ denotes the number of estimate error, and e¯ denotes the mean of all estimate error.Figure 7
Attitude estimate errors of UPS and URTSS.Figure 8
Velocity estimate errors of UPS and URTSS.Figure 9
Position estimate errors of UPS and URTSS.Figure 10
Baseline estimate errors of UPS and URTSS.Table 1
Means of RMSE and STD values of estimate errors.
Parameter
URTSS
UPS
RMSE
STD
RMSE
STD
Attitude error
Heading (°)
3.7874
1.2914
0.4959
0.0217
Pitch (°)
0.5166
0.1218
0.0766
0.0781
Roll (°)
0.8238
0.0371
0.0569
0.0568
Velocity error
East (m/s)
0.0264
0.0245
0.0133
0.0134
North (m/s)
0.0403
0.0414
0.0225
0.0148
Up (m/s)
0.0125
0.0124
0.0224
0.0211
Position error
Latitude (m)
2.1444
0.0990
1.0077
0.0529
Longitude (m)
1.7280
0.0878
0.0510
0.0495
Height (m)
0.1078
0.0987
0.1294
0.0126
Baseline error
(m)
0.7039
0.0532
0.0014
0.0008For attitude estimation, from Figure7 and Table 1, we can see that there is a large and obvious trend item in the heading estimation error of URTSS, and the estimation error of URTSS on pitch and roll is prone to be affected by maneuver and their value are bigger than UPS even after 180° clockwise. Compared with URTSS, when the initial misalignment angle is 40°, the attitude estimation error of UPS is significantly reduced and the error curves are more stable both on AB and CD segments; the estimation accuracy on heading and pitch of UPS is obviously better than that of URTSS, and only the STD on roll of UPS is slightly lower than that of URTSS.For velocity estimation, Figure8 and Table 1 show that the estimation accuracy on east velocity and north velocity of UPS is better than that of URTSS and the error curves of URTSS are prone to be affected by maneuver, especially in the initial ‘S’ shape flight, the estimation error of velocity of URTSS is unstable and very large; the estimation accuracy on up velocity of UPS is lower than that of URTSS.For position estimation, Figure9 and Table 1 show that the estimation accuracy on the latitude and longitude of UPS is better than that of URTSS, and the estimation accuracy on the height of UPS is lower than that of URTSS, which corresponds to the up velocity error result. The estimation error of latitude and longitude of URTSS in the initial ‘S’ shape flight (from 0 s to 300 s) is very large due to the large velocity error in this stage.Besides, the space distance between the main system and subsystem (i.e., the baseline between two antennas of SAR) is a very important parameter for interference precision of interferometric SAR and the higher the measurement accuracy of the baseline, the higher the object’s digital elevation model accuracy of the interferometric SAR. Here, the baseline error is calculated by theoretical position and real position output of the main system and subsystem. From Figure10 and Table 1, we can see that the baseline error curve of UPS is more stable and the baseline estimation accuracy of UPS is better than that of URTSS both on RMSE and STD.From the simulation result, we can see that the estimation accuracy of motion parameters of the proposed method is higher than that of the URTSS as a whole.


## 4.1. Design of Simulation

In this subsection, a typical “S + U” shaped flight trajectory of airborne earth observation is designed. Figure4 shows the plane trajectory and trajectory parameters, respectively. In Figure 4, the AB and CD segment can be considered as the imaging segment.Figure 4
Plane trajectory with S-shaped maneuver plus U-shaped flight.Parameters of flight trajectory are listed as follows: initial latitude is 40°, longitude is 116°, and height is 500 m; initial flight velocity is 100 m/s; and initial heading angle, pitch angle, and roll angle are 40°, 0°, and 0°, respectively. Firstly, the aircraft flies 100 seconds at constant velocity; secondly, the aircraft turns 60° clockwise (100 seconds) and then turns 60° anticlockwise (100 seconds); and thirdly, the aircraft flies 400 seconds at constant velocity, finally turns 180° clockwise (100 seconds), and continually flies 400 seconds at constant velocity. Total flight time is 1300 s.


### 4.1.1. Data Generation

A trajectory generator is used to generate the theoretical data of the scheduled flight trajectory, which includes position, velocity, attitude, and the output data of gyros and accelerometers. The real outputs of the main system are obtained by adding the corresponding measurement noise to the theoretical position, velocity, and attitude. Then, the theoretical outputs of gyros and accelerometers are converted by rigid misalignment angles and flexure angles, and the constant noise and random noise are added to be the inertial sensor outputs of the subsystem. Meanwhile, the motion parameter benchmarks of subsystem can be obtained through transforming theoretical position, velocity, and attitude by flexure angle.


### 4.1.2. Parameters of the Main System and Subsystem

The measurement noises of the main system at heading, pitch, roll, and velocity are0.005°(1σ), 0.0025°(1σ), 0.0025°(1σ), and 0.005m/s (1σ), respectively. Both gyro constant drift and random drift of subsystem are 0.01°/h. Both accelerometer constant bias and random bias of the subsystem are 50μg.


### 4.1.3. Other Parameters

The rigid misalignment angle of the subsystem relative to the main system isρ=0.5∘0.5∘0.5∘T, and the lever arm between the main system and subsystem is rbm=5m0.1m0.1mT. For multiantenna SAR, these antennae are fixed under both sides of the wing. So the flexure angle rotated around the vertical axis is large, while the flexure angles around the other two axes are small. Accordingly, the correlation time of the second-order Markov processes is selected as 2 s, 5 s, and 2 s, and the covariances of flexure angles are 0.01, 0.15, and 0.01, respectively. The curves of flexure angle and flexure angle rate selected in this simulation are shown in Figures 5 and 6. All the initial attitude errors are 40∘. The particle number is 100. The data update rate of the subsystem is 100 Hz and the filter frequency is 1 Hz.Figure 5
Flexure angle onx-axis, y-axis, and z-axis.Figure 6
Flexure angle rate onx-axis, y-axis, and z-axis.


## 4.1.1. Data Generation

A trajectory generator is used to generate the theoretical data of the scheduled flight trajectory, which includes position, velocity, attitude, and the output data of gyros and accelerometers. The real outputs of the main system are obtained by adding the corresponding measurement noise to the theoretical position, velocity, and attitude. Then, the theoretical outputs of gyros and accelerometers are converted by rigid misalignment angles and flexure angles, and the constant noise and random noise are added to be the inertial sensor outputs of the subsystem. Meanwhile, the motion parameter benchmarks of subsystem can be obtained through transforming theoretical position, velocity, and attitude by flexure angle.


## 4.1.2. Parameters of the Main System and Subsystem

The measurement noises of the main system at heading, pitch, roll, and velocity are0.005°(1σ), 0.0025°(1σ), 0.0025°(1σ), and 0.005m/s (1σ), respectively. Both gyro constant drift and random drift of subsystem are 0.01°/h. Both accelerometer constant bias and random bias of the subsystem are 50μg.


## 4.1.3. Other Parameters

The rigid misalignment angle of the subsystem relative to the main system isρ=0.5∘0.5∘0.5∘T, and the lever arm between the main system and subsystem is rbm=5m0.1m0.1mT. For multiantenna SAR, these antennae are fixed under both sides of the wing. So the flexure angle rotated around the vertical axis is large, while the flexure angles around the other two axes are small. Accordingly, the correlation time of the second-order Markov processes is selected as 2 s, 5 s, and 2 s, and the covariances of flexure angles are 0.01, 0.15, and 0.01, respectively. The curves of flexure angle and flexure angle rate selected in this simulation are shown in Figures 5 and 6. All the initial attitude errors are 40∘. The particle number is 100. The data update rate of the subsystem is 100 Hz and the filter frequency is 1 Hz.Figure 5
Flexure angle onx-axis, y-axis, and z-axis.Figure 6
Flexure angle rate onx-axis, y-axis, and z-axis.


## 4.2. Simulation Results and Analysis

The differences between theoretical motion parameters of subsystem and their estimates obtained from estimation, called estimate errors, are used to assess and compare the performance of UPS and URTSS.Figures7–10 show the estimate error curves of UPS and URTSS, including the estimate errors of attitude, velocity, position, and baseline. The means of root mean square error (RMSE) and standard deviation (STD) values of motion information estimate errors in imaging segments AB and CD are calculated and shown in Table 1. The calculation formulas of RMSE and STD are given as follows:
(21)RMSE=∑i=1n′ei2n′,STD=∑i=1n′ei−e¯2n′,where ei denotes the ith estimate error, n′ denotes the number of estimate error, and e¯ denotes the mean of all estimate error.Figure 7
Attitude estimate errors of UPS and URTSS.Figure 8
Velocity estimate errors of UPS and URTSS.Figure 9
Position estimate errors of UPS and URTSS.Figure 10
Baseline estimate errors of UPS and URTSS.Table 1
Means of RMSE and STD values of estimate errors.
Parameter
URTSS
UPS
RMSE
STD
RMSE
STD
Attitude error
Heading (°)
3.7874
1.2914
0.4959
0.0217
Pitch (°)
0.5166
0.1218
0.0766
0.0781
Roll (°)
0.8238
0.0371
0.0569
0.0568
Velocity error
East (m/s)
0.0264
0.0245
0.0133
0.0134
North (m/s)
0.0403
0.0414
0.0225
0.0148
Up (m/s)
0.0125
0.0124
0.0224
0.0211
Position error
Latitude (m)
2.1444
0.0990
1.0077
0.0529
Longitude (m)
1.7280
0.0878
0.0510
0.0495
Height (m)
0.1078
0.0987
0.1294
0.0126
Baseline error
(m)
0.7039
0.0532
0.0014
0.0008For attitude estimation, from Figure7 and Table 1, we can see that there is a large and obvious trend item in the heading estimation error of URTSS, and the estimation error of URTSS on pitch and roll is prone to be affected by maneuver and their value are bigger than UPS even after 180° clockwise. Compared with URTSS, when the initial misalignment angle is 40°, the attitude estimation error of UPS is significantly reduced and the error curves are more stable both on AB and CD segments; the estimation accuracy on heading and pitch of UPS is obviously better than that of URTSS, and only the STD on roll of UPS is slightly lower than that of URTSS.For velocity estimation, Figure8 and Table 1 show that the estimation accuracy on east velocity and north velocity of UPS is better than that of URTSS and the error curves of URTSS are prone to be affected by maneuver, especially in the initial ‘S’ shape flight, the estimation error of velocity of URTSS is unstable and very large; the estimation accuracy on up velocity of UPS is lower than that of URTSS.For position estimation, Figure9 and Table 1 show that the estimation accuracy on the latitude and longitude of UPS is better than that of URTSS, and the estimation accuracy on the height of UPS is lower than that of URTSS, which corresponds to the up velocity error result. The estimation error of latitude and longitude of URTSS in the initial ‘S’ shape flight (from 0 s to 300 s) is very large due to the large velocity error in this stage.Besides, the space distance between the main system and subsystem (i.e., the baseline between two antennas of SAR) is a very important parameter for interference precision of interferometric SAR and the higher the measurement accuracy of the baseline, the higher the object’s digital elevation model accuracy of the interferometric SAR. Here, the baseline error is calculated by theoretical position and real position output of the main system and subsystem. From Figure10 and Table 1, we can see that the baseline error curve of UPS is more stable and the baseline estimation accuracy of UPS is better than that of URTSS both on RMSE and STD.From the simulation result, we can see that the estimation accuracy of motion parameters of the proposed method is higher than that of the URTSS as a whole.


## 5. Conclusion

In this paper, a nonlinear smoother called UPS has been proposed to deal with the off-line transfer alignment estimation of airborne distributed POS. In the UPS, UKF, PF, and smoothing are combined together to provide the respective advantages from each of them to improve estimation accuracy. The validity of the proposed algorithm is verified by simulation test and comparison with URTSS on the estimation accuracy of motion parameters and baseline. The simulation results show that UPS can improve estimation accuracy and has more adaptability to the maneuver than that of URTSS. This method is a better choice for off-line transfer alignment in distributed POS with large misalignment and will greatly improve the resolution of flexible baseline interferometric SAR imaging. The next work is that the practical flight experiment based on long flexible baseline will be implemented to validate the performance of this method.



---

*Source: 3898734-2018-09-30.xml*