# Airborne Integrated Navigation System Based on SINS/GPS/TAN/EOAN

Authors: Junjun Tang; Peijuan Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Considering the drawbacks that GPS signal is susceptible to obstacles and TAN becomes useless in some area when without any terrain data or with a featureless terrain field, to realize long-distance and high-precision navigation, a navigation system based on SINS/GPS/TAN/EOAN is presented. When GPS signal is available, GPS is used to correct errors of SINS; when GPS is unavailable, a terrain selection method based on the entropy weighted gray relational decision-making method is use to distinguish terrain into matchable areas and unmatchable areas; then, for the matchable areas, TAN is used to correct errors of SINS, for the unmatchable areas, EOAN is used to correct errors of SINS. The principles of SINS, GPS, TAN, and EOAN are analyzed, the mathematic models of SINS/GPS, SINS/TAN, and SINS/EOAN are constructed, and finally the federated Kalman filter is used to fuse navigation information. Simulation results show that the trajectory of SINS/GPS/TAN/EOAN is close to the ideal one in both matchable area or unmatchable area and whose navigation errors are obviously reduced, which is important for the realization of long-time and high-precision positioning.


---

## Body

## 1. Introduction

The Strapdown Inertial Navigation System (SINS) is an autonomous navigation system, which neither relies on any external information nor radiates energy to the outside world [1]. Besides, it has the advantages of high precision, good concealment, no susceptibility to interference, and so on. However, due to its navigation characteristics, the error of SINS inevitably grows with time. To realize long-distance and high-precision navigation, several aided navigation systems are used to correct its errors.Global Positioning System (GPS) is a satellite-based radio navigation system that can provide accurate velocity and position information [2, 3]. It has advantages of high precision, short observation time, and so on. However, it is susceptibility to interference, the signals might be blocked and attenuated, and when GPS information is unavailable, the error of SINS is unable to be corrected.Terrain aided navigation system is an ideal auxiliary navigation method, which is based on terrain elevation characteristics [4, 5]. The error of TAN is not accumulated while it has advantages of autonomy, concealment, continuity, and so on. However, it cannot achieve good effect in the areas where the terrain is flat or is of less terrain features. So, to make full use of TAN, it is necessary to distinguish terrain into matchable areas and unmatchable areas by measuring terrain feature parameters [6–9]. A novel selection criteria based on entropy weighted gray relational decision-making is used to divide terrain into matchable areas and unmatchable areas, in which the influences of terrain parameters such as standard deviation, roughness, entropy, and correlation coefficient are comprehensively considered.Electro-optical Aided Navigation (EOAN) System is based on image matching technology that can provide position information [10]. It has advantages of high precision and strong anti-interference ability. For the area where there is less of terrain information, EOAN can be used to correct errors of SINS [11].Integration navigation technology [12–15] has become one of the key issues because of its capability of overcoming the drawbacks of single navigation method as above [16, 17]. And it is necessary to research multisensor fusion technologies [18]. In the multisensor information fusion system, distributed filters have attracted more and more attention of researchers [19, 20]. The federated Kalman filter, which is one kind of the distributed filtering technique, is proposed by Carlson, and it has been widely used in multi-information fusion because of its flexibility and fault tolerance [15, 21–24].The ideal of the integrated navigation system is shown as Figure1. For the areas where GPS is available, GPS is used to correct errors of SINS; for the areas where GPS is unavailable but is of full terrain information, TAN is used to correct errors of SINS; and for the areas where GPS is unavailable and is of less terrain information, EOAN is used to correct errors of SINS.Figure 1
Airborne integrated navigation system.The paper is organized as follows. In Section2, terrain parameters are described and the selection method is introduced. In Section 3, the principles of SINS, GPS, TAN, and EOAN are analyzed and the mathematic models of SINS/GPS, SINS/TAN, and SINS/EOAN are constructed. In Section 4, the federated Kalman filter is proposed. In Section 5, several tests are conducted to evaluate the performance of SINS/GPS/TAN/EOAN with FKF. And conclusions are made in Section 6.


## 2. Terrain Selection Method

### 2.1. Terrain Parameters

The terrain data is generally stored in the form of grid matrix. Assume that(A)
The exemplified terrain hasm×n grids(B)
The coordinates of the grid point isi,j(C)
The height of the grid points ishi,jTerrain parameters such as standard deviation, roughness, entropy, and correlation coefficient can be defined as follows:


#### 2.1.1. Standard Deviation

The terrain standard deviation reflects the range of terrain elevation deviating from the average value. It is a macroparameter reflecting the severity of terrain overall fluctuation. And it can be defined as follows:(1)σ=1mn−1∑i=1m∑j=1nhi,j−h_ave2,h_ave=1mn∑i=1m∑j=1nhi,j,where h_ave means the average value of terrain elevation.


#### 2.1.2. Correlation Coefficient

Terrain correlation coefficient is a characteristic parameter used to reflect terrain correlation, and it can be defined as follows:(2)R=12Rlong+Rlati,Rlong=1nm−1σ2∑i=1m−1∑j=1nhi,j−h_ave•hi+1,j−h_ave,Rlati=1nm−1σ2∑i=1m∑j=1n−1hi,j−h_ave•hi,j+1−h_ave,where Rlong means correlation coefficient of longitude direction, which reflects the correlation degree in latitude direction of terrain profile in longitude direction, and Rlati means correlation coefficient of latitude, which reflects the correlation degree in longitude direction of terrain profile in latitude direction.


#### 2.1.3. Roughness

Terrain roughness reflects the smoothness of terrain and is a characteristic parameter used to describe the local undulation, and it can be defined as follows:(3)r=12rlong+rlati,rlong=1nm−1∑i=1m−1∑j=1nhi,j−hi+1,j,rlati=1mn−1∑i=1m∑j=1n−1hi,j−hi,j+1,where rlong represents the roughness of longitude and rlati represents the roughness of latitude.


#### 2.1.4. Terrain Height Entropy

Terrain height entropy is a characteristic parameter reflecting the amount of information. The more obvious the terrain height changes, the smaller the height entropy. The height entropy can be described as follows:(4)Hf=−∑i=1m∑j=1npi,jlogpi,j,pi,j=hi,j∑i=1m∑j=1nhi,j,where pi,j means the normalized height value.


### 2.2. Entropy Weighted Gray Relational Decision-Making

Different parameters have different effects on the abundance of terrain information. Under certain conditions, the smaller the height entropy and the correlation coefficient are, the more abundant the terrain is and the greater the adaptability is. Conversely, the bigger the standard deviation and the roughness are, the more abundant the terrain is and the greater the adaptability is. Therefore, it is necessary to take all parameters into consideration and adopt a method which is able to fuse multifeatures. Assume thatT=σ,r,Hf,R is a feature vector after fused, according to which the adaptability of the terrain can be researched and a more comprehensive selection method is obtained.


#### 2.2.1. Construction of Decision Matrix

Assume that the number of terrain candidate areas ism, which are X1,X2,…,Xm, and the number of feature parameters in every candidate area is n, which are G1,G2,…Gn. The property of candidate area in parameter is xij, and the decision matrix can be obtained:(5)X=x11x12…x1nx21x22…x2n⋮⋮⋱⋮xm1xm2…xmn.


#### 2.2.2. Normalized Matrix

Regarding correlation coefficient and height entropy as cost indicators and standard deviation and roughness as benefit indictors, the decision matrixY=yijm×n can be normalized as follows:(6)yij=xij−minixijmaxixij−minixijbenefit indictors,yij=maxixij−xijmaxixij−minixijcost indictors.


#### 2.2.3. Grey Relational Matrix

Assume that the optimal solution of indictors isyj∗,j=1,2,…,n, the reference sequence consisting of yj∗, and can be presented as y0j=yj∗|j=1,2,…,n. Making yij=yij|j=1,2,…,n,i=1,2,…,m to be comparison sequence, the correlation coefficient of yi and y0 in indictor can be calculated as follows:(7)ξij=Δmin+ρΔmaxΔi0j+ρΔmax,where Δmin=miniminjy0j−yij and Δmax=maximaxjy0j−yij, respectively, means the minimum value and the maximum value of absolute differences of comparison sequence. Δi0j=y0j−yij represents the absolute differences of comparison sequence and ρ is the resolution coefficient, where 0<ρ<1. Assume that ρ=0.5.The gray relational matrix can be as follows:(8)ξ=ξ11ξ12…ξ1nξ21ξ22…ξ2n⋮⋮⋱⋮ξm1ξm2…ξmn.


#### 2.2.4. Weight Calculation of Entropy

(A)
The calculation of entropy requires logarithm, to avoid meaningless data, the data in the normalized matrix can be translated by adding 1, and the contribution of terraini to parameter j can be calculated as follows:(9)fij=yij+1∑i=1myij+1.(B)
The entropyEj can be used to represent the total contribution of all terrains to indictor j:(10)Ej=−k∑i=1mfijlnfij,j=1,2,…,n,wherek>0， and assume that k=1/lnm,Ej∈0,1. For conditions that the contributions of each terrain to indictor j is consistent, Ej tends to be 1, and the indictor is useless in the decision-making, and the weight value can be 0. Reversely, for conditions that the contributions of each terrain to indictor j are quietly different, the indictor is more useful to decision-making, and the weight value should be lager.(C)
Definegj as the difference coefficient of contribution of terrains to indictor j, and the indictor weight can be calculated as follows:(11)gj=1−Ej,j=1,2,…,n,wj=gj∑j=1ngj,j=1,2,…,n.


#### 2.2.5. Comprehensive Evaluation Value of Gray Correlation Degree

The comprehensive evaluation value of gray correlation degree of terrainXi is as follows:(12)Di=∑j=1nξijwj.The comprehensive evaluation values can be sorted; the bigger the evaluation value is, the higher the adaptability of the terrain is.


## 2.1. Terrain Parameters

The terrain data is generally stored in the form of grid matrix. Assume that(A)
The exemplified terrain hasm×n grids(B)
The coordinates of the grid point isi,j(C)
The height of the grid points ishi,jTerrain parameters such as standard deviation, roughness, entropy, and correlation coefficient can be defined as follows:


### 2.1.1. Standard Deviation

The terrain standard deviation reflects the range of terrain elevation deviating from the average value. It is a macroparameter reflecting the severity of terrain overall fluctuation. And it can be defined as follows:(1)σ=1mn−1∑i=1m∑j=1nhi,j−h_ave2,h_ave=1mn∑i=1m∑j=1nhi,j,where h_ave means the average value of terrain elevation.


### 2.1.2. Correlation Coefficient

Terrain correlation coefficient is a characteristic parameter used to reflect terrain correlation, and it can be defined as follows:(2)R=12Rlong+Rlati,Rlong=1nm−1σ2∑i=1m−1∑j=1nhi,j−h_ave•hi+1,j−h_ave,Rlati=1nm−1σ2∑i=1m∑j=1n−1hi,j−h_ave•hi,j+1−h_ave,where Rlong means correlation coefficient of longitude direction, which reflects the correlation degree in latitude direction of terrain profile in longitude direction, and Rlati means correlation coefficient of latitude, which reflects the correlation degree in longitude direction of terrain profile in latitude direction.


### 2.1.3. Roughness

Terrain roughness reflects the smoothness of terrain and is a characteristic parameter used to describe the local undulation, and it can be defined as follows:(3)r=12rlong+rlati,rlong=1nm−1∑i=1m−1∑j=1nhi,j−hi+1,j,rlati=1mn−1∑i=1m∑j=1n−1hi,j−hi,j+1,where rlong represents the roughness of longitude and rlati represents the roughness of latitude.


### 2.1.4. Terrain Height Entropy

Terrain height entropy is a characteristic parameter reflecting the amount of information. The more obvious the terrain height changes, the smaller the height entropy. The height entropy can be described as follows:(4)Hf=−∑i=1m∑j=1npi,jlogpi,j,pi,j=hi,j∑i=1m∑j=1nhi,j,where pi,j means the normalized height value.


## 2.1.1. Standard Deviation

The terrain standard deviation reflects the range of terrain elevation deviating from the average value. It is a macroparameter reflecting the severity of terrain overall fluctuation. And it can be defined as follows:(1)σ=1mn−1∑i=1m∑j=1nhi,j−h_ave2,h_ave=1mn∑i=1m∑j=1nhi,j,where h_ave means the average value of terrain elevation.


## 2.1.2. Correlation Coefficient

Terrain correlation coefficient is a characteristic parameter used to reflect terrain correlation, and it can be defined as follows:(2)R=12Rlong+Rlati,Rlong=1nm−1σ2∑i=1m−1∑j=1nhi,j−h_ave•hi+1,j−h_ave,Rlati=1nm−1σ2∑i=1m∑j=1n−1hi,j−h_ave•hi,j+1−h_ave,where Rlong means correlation coefficient of longitude direction, which reflects the correlation degree in latitude direction of terrain profile in longitude direction, and Rlati means correlation coefficient of latitude, which reflects the correlation degree in longitude direction of terrain profile in latitude direction.


## 2.1.3. Roughness

Terrain roughness reflects the smoothness of terrain and is a characteristic parameter used to describe the local undulation, and it can be defined as follows:(3)r=12rlong+rlati,rlong=1nm−1∑i=1m−1∑j=1nhi,j−hi+1,j,rlati=1mn−1∑i=1m∑j=1n−1hi,j−hi,j+1,where rlong represents the roughness of longitude and rlati represents the roughness of latitude.


## 2.1.4. Terrain Height Entropy

Terrain height entropy is a characteristic parameter reflecting the amount of information. The more obvious the terrain height changes, the smaller the height entropy. The height entropy can be described as follows:(4)Hf=−∑i=1m∑j=1npi,jlogpi,j,pi,j=hi,j∑i=1m∑j=1nhi,j,where pi,j means the normalized height value.


## 2.2. Entropy Weighted Gray Relational Decision-Making

Different parameters have different effects on the abundance of terrain information. Under certain conditions, the smaller the height entropy and the correlation coefficient are, the more abundant the terrain is and the greater the adaptability is. Conversely, the bigger the standard deviation and the roughness are, the more abundant the terrain is and the greater the adaptability is. Therefore, it is necessary to take all parameters into consideration and adopt a method which is able to fuse multifeatures. Assume thatT=σ,r,Hf,R is a feature vector after fused, according to which the adaptability of the terrain can be researched and a more comprehensive selection method is obtained.


### 2.2.1. Construction of Decision Matrix

Assume that the number of terrain candidate areas ism, which are X1,X2,…,Xm, and the number of feature parameters in every candidate area is n, which are G1,G2,…Gn. The property of candidate area in parameter is xij, and the decision matrix can be obtained:(5)X=x11x12…x1nx21x22…x2n⋮⋮⋱⋮xm1xm2…xmn.


### 2.2.2. Normalized Matrix

Regarding correlation coefficient and height entropy as cost indicators and standard deviation and roughness as benefit indictors, the decision matrixY=yijm×n can be normalized as follows:(6)yij=xij−minixijmaxixij−minixijbenefit indictors,yij=maxixij−xijmaxixij−minixijcost indictors.


### 2.2.3. Grey Relational Matrix

Assume that the optimal solution of indictors isyj∗,j=1,2,…,n, the reference sequence consisting of yj∗, and can be presented as y0j=yj∗|j=1,2,…,n. Making yij=yij|j=1,2,…,n,i=1,2,…,m to be comparison sequence, the correlation coefficient of yi and y0 in indictor can be calculated as follows:(7)ξij=Δmin+ρΔmaxΔi0j+ρΔmax,where Δmin=miniminjy0j−yij and Δmax=maximaxjy0j−yij, respectively, means the minimum value and the maximum value of absolute differences of comparison sequence. Δi0j=y0j−yij represents the absolute differences of comparison sequence and ρ is the resolution coefficient, where 0<ρ<1. Assume that ρ=0.5.The gray relational matrix can be as follows:(8)ξ=ξ11ξ12…ξ1nξ21ξ22…ξ2n⋮⋮⋱⋮ξm1ξm2…ξmn.


### 2.2.4. Weight Calculation of Entropy

(A)
The calculation of entropy requires logarithm, to avoid meaningless data, the data in the normalized matrix can be translated by adding 1, and the contribution of terraini to parameter j can be calculated as follows:(9)fij=yij+1∑i=1myij+1.(B)
The entropyEj can be used to represent the total contribution of all terrains to indictor j:(10)Ej=−k∑i=1mfijlnfij,j=1,2,…,n,wherek>0， and assume that k=1/lnm,Ej∈0,1. For conditions that the contributions of each terrain to indictor j is consistent, Ej tends to be 1, and the indictor is useless in the decision-making, and the weight value can be 0. Reversely, for conditions that the contributions of each terrain to indictor j are quietly different, the indictor is more useful to decision-making, and the weight value should be lager.(C)
Definegj as the difference coefficient of contribution of terrains to indictor j, and the indictor weight can be calculated as follows:(11)gj=1−Ej,j=1,2,…,n,wj=gj∑j=1ngj,j=1,2,…,n.


### 2.2.5. Comprehensive Evaluation Value of Gray Correlation Degree

The comprehensive evaluation value of gray correlation degree of terrainXi is as follows:(12)Di=∑j=1nξijwj.The comprehensive evaluation values can be sorted; the bigger the evaluation value is, the higher the adaptability of the terrain is.


## 2.2.1. Construction of Decision Matrix

Assume that the number of terrain candidate areas ism, which are X1,X2,…,Xm, and the number of feature parameters in every candidate area is n, which are G1,G2,…Gn. The property of candidate area in parameter is xij, and the decision matrix can be obtained:(5)X=x11x12…x1nx21x22…x2n⋮⋮⋱⋮xm1xm2…xmn.


## 2.2.2. Normalized Matrix

Regarding correlation coefficient and height entropy as cost indicators and standard deviation and roughness as benefit indictors, the decision matrixY=yijm×n can be normalized as follows:(6)yij=xij−minixijmaxixij−minixijbenefit indictors,yij=maxixij−xijmaxixij−minixijcost indictors.


## 2.2.3. Grey Relational Matrix

Assume that the optimal solution of indictors isyj∗,j=1,2,…,n, the reference sequence consisting of yj∗, and can be presented as y0j=yj∗|j=1,2,…,n. Making yij=yij|j=1,2,…,n,i=1,2,…,m to be comparison sequence, the correlation coefficient of yi and y0 in indictor can be calculated as follows:(7)ξij=Δmin+ρΔmaxΔi0j+ρΔmax,where Δmin=miniminjy0j−yij and Δmax=maximaxjy0j−yij, respectively, means the minimum value and the maximum value of absolute differences of comparison sequence. Δi0j=y0j−yij represents the absolute differences of comparison sequence and ρ is the resolution coefficient, where 0<ρ<1. Assume that ρ=0.5.The gray relational matrix can be as follows:(8)ξ=ξ11ξ12…ξ1nξ21ξ22…ξ2n⋮⋮⋱⋮ξm1ξm2…ξmn.


## 2.2.4. Weight Calculation of Entropy

(A)
The calculation of entropy requires logarithm, to avoid meaningless data, the data in the normalized matrix can be translated by adding 1, and the contribution of terraini to parameter j can be calculated as follows:(9)fij=yij+1∑i=1myij+1.(B)
The entropyEj can be used to represent the total contribution of all terrains to indictor j:(10)Ej=−k∑i=1mfijlnfij,j=1,2,…,n,wherek>0， and assume that k=1/lnm,Ej∈0,1. For conditions that the contributions of each terrain to indictor j is consistent, Ej tends to be 1, and the indictor is useless in the decision-making, and the weight value can be 0. Reversely, for conditions that the contributions of each terrain to indictor j are quietly different, the indictor is more useful to decision-making, and the weight value should be lager.(C)
Definegj as the difference coefficient of contribution of terrains to indictor j, and the indictor weight can be calculated as follows:(11)gj=1−Ej,j=1,2,…,n,wj=gj∑j=1ngj,j=1,2,…,n.


## 2.2.5. Comprehensive Evaluation Value of Gray Correlation Degree

The comprehensive evaluation value of gray correlation degree of terrainXi is as follows:(12)Di=∑j=1nξijwj.The comprehensive evaluation values can be sorted; the bigger the evaluation value is, the higher the adaptability of the terrain is.


## 3. Mathematical Model

### 3.1. Error Model of SINS

Here, we define the navigation coordinates frames as an east-north-vertical geography coordinates frames. The position errors, velocity errors, attitude angle errors, gyro drift, and accelerometer bias are formed as the state vector. The state vector is determined by(13)XSINS=δLδλδVEδVNΦEΦNΦU∇bx∇byεbxεbyεbzT,where δL and δλ represent the position errors, δVE and δVN are east and north velocity errors, φE, φN, andφU are misalignment angles of pitch, roll, and yaw, ∇bx and ∇by represent accelerometer bias, and εbx, εby, andεbz are gyro drifts along x, y, and z axes.The system state equation can be constructed as(14)X˙SINS=FSINSXSINS+WSINS,where(15)WSINS=ωδLωδλωVEωVNωϕEωϕNωϕU00000,F21=VERsec2LsinL,F31=2ωiecosL+VERsec2LVN,F41=−2ωiecosL+VERsec2LVE,F34=2ωiesinL+VERtanL,F43=−2ωiesinL+VERtanL,F57=−F75=−ωiecosL+VER,F56=−F65=ωiesinL+VERtanL,F61=−ωiesinL,F71=ωiecosL+VERsec2L.Cij is element of Cnb, which is the direction cosine matrix from n-frame to b-frame:(16)FSINS=0001R00000000F210secLL000000000F310VNtanLRF340−fUfNC11C21000F410F430fU0−fEC12C22000000−1R0F56F5700−C11−C21−C31F6101R0F650−VNR00−C12−C22−C32F710tanLR0F75VNR000−C13−C23−C33000000000000000000000000000000000000000000000000000000000000.


### 3.2. Mathematical Model of SINS/GPS

In the areas where GPS is available, GPS is used to correct errors of SINS. Take the difference between the position and velocity information of the INS and the GPS as the observation measurement and establish the system measurement equation as follows:(17)Zk_GPS=LSINS−LGPSλSINS−λGPSVE_SINS−VE_GPSVN_SINS−VN_GPS=Hk_GPSXSINSXGPS+Vk_GPS,where LSINS, λSINS, VE_SINS, and VN_SINS are position and velocity information of SINS, LGPS, λGPS, VE_GPS, and VN_GPS are position and velocity information of GPS, Vk_GPS is the measurement noise matrix, and Hk_GPS is the measurement matrix, and the measurement matrix can be expressed as(18)Hk_GPS=10000000000000−1000010000000000000−1000010000000000000−1000010000000000000−1.


### 3.3. Mathematical Model of SINS/TAN

In the areas where GPS is unavailable but is of enough terrain information, TAN is used to correct errors of SINS. The aircraft is assumed to be equipped with a radio altimeter to collect real-time height data during the navigation. According to the height map stored in advance, the matching algorithm is able to find the best nearest point, whose longitude and latitude are the output of TAN. Take the difference between the position information of the INS and the TAN as the observation measurement and establish the system measurement equation as follows:(19)Zk_TAN=LSINS−LTANλSINS−λTAN=Hk_TANXSINSXTAN+Vk_TAN,where LTAN and λTAN are position information outputted by TAN. And Vk_TAN is the measurement noise matrix and Hk_TAN is the measurement matrix, which can be expressed as(20)Hk_TAN=100000000000−100100000000000−1.


### 3.4. Mathematical Model of SINS/EOAN

In the areas where GPS is unavailable and is of less terrain information, EOAN is used to correct errors of SINS [11]. The aircraft is assumed to be equipped with a video camera to collect real-time image data during the navigation. Compared with the image database stored in advance and according to the hypothesis that the gray value of each point in the above image is equal, the distance between the target and aircraft can be calculated, and the real-time position of the aircraft can be converted. Take the difference between the position information of the INS and the EOAN as the observation measurement and establish the system measurement equation as follows:(21)Zk_EOAN=LSINS−LEOANλSINS−λEOAN=Hk_EOANXSINSXEOAN+Vk_EOAN,where LEOAN and λEOAN are position information obtained by EOAN, Vk_EOAN is the measurement noise matrix, and Hk_EOAN is the measurement matrix, which can be expressed as(22)Hk_EOAN=100000000000−100100000000000−1.


## 3.1. Error Model of SINS

Here, we define the navigation coordinates frames as an east-north-vertical geography coordinates frames. The position errors, velocity errors, attitude angle errors, gyro drift, and accelerometer bias are formed as the state vector. The state vector is determined by(13)XSINS=δLδλδVEδVNΦEΦNΦU∇bx∇byεbxεbyεbzT,where δL and δλ represent the position errors, δVE and δVN are east and north velocity errors, φE, φN, andφU are misalignment angles of pitch, roll, and yaw, ∇bx and ∇by represent accelerometer bias, and εbx, εby, andεbz are gyro drifts along x, y, and z axes.The system state equation can be constructed as(14)X˙SINS=FSINSXSINS+WSINS,where(15)WSINS=ωδLωδλωVEωVNωϕEωϕNωϕU00000,F21=VERsec2LsinL,F31=2ωiecosL+VERsec2LVN,F41=−2ωiecosL+VERsec2LVE,F34=2ωiesinL+VERtanL,F43=−2ωiesinL+VERtanL,F57=−F75=−ωiecosL+VER,F56=−F65=ωiesinL+VERtanL,F61=−ωiesinL,F71=ωiecosL+VERsec2L.Cij is element of Cnb, which is the direction cosine matrix from n-frame to b-frame:(16)FSINS=0001R00000000F210secLL000000000F310VNtanLRF340−fUfNC11C21000F410F430fU0−fEC12C22000000−1R0F56F5700−C11−C21−C31F6101R0F650−VNR00−C12−C22−C32F710tanLR0F75VNR000−C13−C23−C33000000000000000000000000000000000000000000000000000000000000.


## 3.2. Mathematical Model of SINS/GPS

In the areas where GPS is available, GPS is used to correct errors of SINS. Take the difference between the position and velocity information of the INS and the GPS as the observation measurement and establish the system measurement equation as follows:(17)Zk_GPS=LSINS−LGPSλSINS−λGPSVE_SINS−VE_GPSVN_SINS−VN_GPS=Hk_GPSXSINSXGPS+Vk_GPS,where LSINS, λSINS, VE_SINS, and VN_SINS are position and velocity information of SINS, LGPS, λGPS, VE_GPS, and VN_GPS are position and velocity information of GPS, Vk_GPS is the measurement noise matrix, and Hk_GPS is the measurement matrix, and the measurement matrix can be expressed as(18)Hk_GPS=10000000000000−1000010000000000000−1000010000000000000−1000010000000000000−1.


## 3.3. Mathematical Model of SINS/TAN

In the areas where GPS is unavailable but is of enough terrain information, TAN is used to correct errors of SINS. The aircraft is assumed to be equipped with a radio altimeter to collect real-time height data during the navigation. According to the height map stored in advance, the matching algorithm is able to find the best nearest point, whose longitude and latitude are the output of TAN. Take the difference between the position information of the INS and the TAN as the observation measurement and establish the system measurement equation as follows:(19)Zk_TAN=LSINS−LTANλSINS−λTAN=Hk_TANXSINSXTAN+Vk_TAN,where LTAN and λTAN are position information outputted by TAN. And Vk_TAN is the measurement noise matrix and Hk_TAN is the measurement matrix, which can be expressed as(20)Hk_TAN=100000000000−100100000000000−1.


## 3.4. Mathematical Model of SINS/EOAN

In the areas where GPS is unavailable and is of less terrain information, EOAN is used to correct errors of SINS [11]. The aircraft is assumed to be equipped with a video camera to collect real-time image data during the navigation. Compared with the image database stored in advance and according to the hypothesis that the gray value of each point in the above image is equal, the distance between the target and aircraft can be calculated, and the real-time position of the aircraft can be converted. Take the difference between the position information of the INS and the EOAN as the observation measurement and establish the system measurement equation as follows:(21)Zk_EOAN=LSINS−LEOANλSINS−λEOAN=Hk_EOANXSINSXEOAN+Vk_EOAN,where LEOAN and λEOAN are position information obtained by EOAN, Vk_EOAN is the measurement noise matrix, and Hk_EOAN is the measurement matrix, which can be expressed as(22)Hk_EOAN=100000000000−100100000000000−1.


## 4. Federated Kalman Filter

### 4.1. System Scheme

The system scheme is shown in Figure2. SINS is the public system, which provides position, velocity, and attitude angle information. The first subfilter is composed of SINS and GPS, the second subfilter is composed of SINS and TAN, and the third subfilter is composed of SINS and EOAN. All the above three subfilters are Kalman filters. The federal filter is composed of the three subfilters and a main filter, where the subfilter can independently update the time and measurement, and the result will be sent to the main filter. And the main filter will feed back to the subfilters after data fusion. It is the initial value of the next cycle. It is important that the local estimations are not related to each other, so the system has good fault tolerance and can meet the requirement of long-distance and high-precision navigation.Figure 2
Design structure of the filter.


### 4.2. The Traditional Kalman Filter

It is assumed that the system equation and measurement equation of the linear discrete system can be written as follows:(23)Xik=Фikk−1Xik−1+Wik,Zik=HikXik+Vik,where Xik is the estimated state of subfilter i, Zik is the measurement of subfilter i, Фik/k−1 is the transfer matrix from tk−1 to tk, Wik is noise sequence of system incentive, Hik is the measurement matrix, Vik is noise sequence of measurement, and EVikViTk=Rik.So, the subfilter equation is shown as follows.(1)
Time update is(24)X^ikk−1=Фikk−1X^ik−1,Pikk−1=Фikk−1Pik−1ФiTkk−1+Γik−1Qik−1ΓiTk−1.(2)
Measurement update is(25)Kik=Pikk−1HiTkHiTkPikk−1HiTk+Rik−1,Pikk=I−KikHikPikk−1I−KikHikT+KikRikKiTk,X^ikk=X^ikk−1+KikZik−HikkX^ikk−1.


### 4.3. The Federal Filter

The fusion algorithm of the federal filter [18] is as follows:(26)Pg=∑i=13Pi−1kk−1,X^g=Pg∑i=13Pi−1kkX^ikk,where X^g is global state estimation after fusion, X^ik/k is the local optimum estimation of subfilter, Pg is the variance matrix of error, and Pik/k is the local optimum estimation.


### 4.4. Information Distribution

The principle of information distribution in the federated Kalman filter is(27)Pi−1k=βiPg−1kX^ik=X^gkQi−1k=βiQg−1ki=1,2,3,∑i=13βi+βm=1,where βi represents the information distribution factor of subfilter and βm represents the information distribution factor of main filter, choosing βm=0.The estimation accuracy of the Federal filter is equivalent to that of the Kalman filter; however, the Federal filter has a small amount of calculation and is of strong fault tolerance. And distribution coefficient is the key to the performance of the federal filter. Considering that the estimation accuracy of each subfilter is decided by the error covariance matrix, the distribution coefficient can be determined be the eigenvalue of the error covariance matrix; thus,(28)βik=trPik−1−1∑i=13trPik−1−1.


## 4.1. System Scheme

The system scheme is shown in Figure2. SINS is the public system, which provides position, velocity, and attitude angle information. The first subfilter is composed of SINS and GPS, the second subfilter is composed of SINS and TAN, and the third subfilter is composed of SINS and EOAN. All the above three subfilters are Kalman filters. The federal filter is composed of the three subfilters and a main filter, where the subfilter can independently update the time and measurement, and the result will be sent to the main filter. And the main filter will feed back to the subfilters after data fusion. It is the initial value of the next cycle. It is important that the local estimations are not related to each other, so the system has good fault tolerance and can meet the requirement of long-distance and high-precision navigation.Figure 2
Design structure of the filter.


## 4.2. The Traditional Kalman Filter

It is assumed that the system equation and measurement equation of the linear discrete system can be written as follows:(23)Xik=Фikk−1Xik−1+Wik,Zik=HikXik+Vik,where Xik is the estimated state of subfilter i, Zik is the measurement of subfilter i, Фik/k−1 is the transfer matrix from tk−1 to tk, Wik is noise sequence of system incentive, Hik is the measurement matrix, Vik is noise sequence of measurement, and EVikViTk=Rik.So, the subfilter equation is shown as follows.(1)
Time update is(24)X^ikk−1=Фikk−1X^ik−1,Pikk−1=Фikk−1Pik−1ФiTkk−1+Γik−1Qik−1ΓiTk−1.(2)
Measurement update is(25)Kik=Pikk−1HiTkHiTkPikk−1HiTk+Rik−1,Pikk=I−KikHikPikk−1I−KikHikT+KikRikKiTk,X^ikk=X^ikk−1+KikZik−HikkX^ikk−1.


## 4.3. The Federal Filter

The fusion algorithm of the federal filter [18] is as follows:(26)Pg=∑i=13Pi−1kk−1,X^g=Pg∑i=13Pi−1kkX^ikk,where X^g is global state estimation after fusion, X^ik/k is the local optimum estimation of subfilter, Pg is the variance matrix of error, and Pik/k is the local optimum estimation.


## 4.4. Information Distribution

The principle of information distribution in the federated Kalman filter is(27)Pi−1k=βiPg−1kX^ik=X^gkQi−1k=βiQg−1ki=1,2,3,∑i=13βi+βm=1,where βi represents the information distribution factor of subfilter and βm represents the information distribution factor of main filter, choosing βm=0.The estimation accuracy of the Federal filter is equivalent to that of the Kalman filter; however, the Federal filter has a small amount of calculation and is of strong fault tolerance. And distribution coefficient is the key to the performance of the federal filter. Considering that the estimation accuracy of each subfilter is decided by the error covariance matrix, the distribution coefficient can be determined be the eigenvalue of the error covariance matrix; thus,(28)βik=trPik−1−1∑i=13trPik−1−1.


## 5. Simulation Experiment

The simulation terrain can be divided into matchable areas and unmatchable areas using the terrain selection method based on entropy weighted gray relational decision-making, and the map of divided terrain is shown as Figure3.Figure 3
Map of divided terrain.The simulation conditions and parameters can be set as Table1.Table 1
Simulation conditions and parameters.
TermsConditions and parametersLatitude initial error of SINS (“)5Longitude initial error of SINS (“)5East velocity error (m/s)0.1North velocity error (m/s)0.1Accelerometer biases (ug)20Gyroscope constant (°/h)0.1Positioning error of TAN (m)50Positioning error of EOAN (m)50In the matchable area, the ideal trajectory starts from (31.0325 N, 118.0300 E) and ends at (31.0435 N, 118.0390 E). And the simulated navigation trajectories of SINS, SINS/TAN, and the proposed method are shown in Figure4.Figure 4
Navigation trajectory in the matchable area.And errors of trajectories in the matchable area are calculated as Table2.Table 2
Errors of trajectories in the matchable area.
TrajectoryErrors of latitudeErrors of longitudeMaximumAverageMaximumAverageSINS186.5326151.0221344.6527337.2541SINS/TAN95.712545.632457.892135.2196SINS/GPS/TAN/EOAN88.306942.678155.632634.1945In the unmatchable area, the ideal trajectory starts from (31.0724 N, 118.0750 E) and ends at (31.0834 N, 118.0840 E). And the simulated navigation trajectories of SINS, SINS/TAN, and the proposed method are shown in Figure5.Figure 5
Navigation trajectory in the unmatchable area.And errors of trajectories in the unmatchable area are calculated as Table3.Table 3
Errors of trajectories in the unmatchable area.
TrajectoryErrors of latitudeErrors of longitudeMaximumAverageMaximumAverageSINS186.5326151.0221344.6527337.2541SINS/TAN350.1225177.5213465.2315215.2196SINS/GPS/TAN/EOAN90.174341.569859.324936.3312In Figures3 and 4, the red line represents ideal trajectory, the blue line represents SINS trajectory, the green line represents SINS/TAN trajectory, and the orange line represents trajectory of the proposed method. Obviously, the presented navigation method can obtain a trajectory which is much closer to the ideal one.It can be seen from Table2 that SINS/TAN can effectively reduce the positioning error in the matchable area, while the positioning accuracy of SINS/GPS/TAN/EOAN is slightly higher than that of SINS/TAN.From Table3, it can be seen that SINS/TAN becomes useless in the unmatchable area, while SINS/GPS/TAN/EOAN can still significantly improve the positioning accuracy.


## 6. Conclusion

The airborne integrated navigation system proposed in this paper is a combination of SINS, GPS, TAN, and OEAN. The federal Kalman filter is used to fuse the navigation information. Compared with the simulation results, it can be seen that there is a lager deviation between the trajectory of pure SINS and that of ideal one; the trajectory of TAN is close to the ideal one in only the matchable area, while the trajectory of SINS/GPS/TAN/EOAN is close to the ideal one in both matchable area or unmatchable area, and whose navigation errors are obviously reduced, which is important for the realization of long-time and high-precision positioning.



---

*Source: 3218942-2020-11-12.xml*