# A New Method to Improve the Measurement Accuracy of Autonomous Astronomical Navigation

Authors: Yu-Feng Yang; Xiang Han; Ning-Ning Song; Zi-Chen Wang
Journal: Journal of Mathematics
 (2022)

---

## Abstract

Based on the theory of atmospheric refraction, combined with the atmospheric parameter data of NCEP (National Centers for environmental prediction), the Fourier interpolation fitting algorithm is used to model and analyze the parameters affecting atmospheric refraction on a global scale. The atmospheric temperature and density model with space-time variation is constructed. The spacecraft state equation and the measurement equation with the starlight apparent height as the observation quantity are established. Moreover, the Unscented Kalman filter is applied to the indirect sensitive horizon autonomous astronomical navigation of starlight refraction. The relative error of fitting the measured data with the spatiotemporal atmospheric temperature model established in this paper is less than 2%. The position estimation error of the navigation system is 94 m, and the velocity estimation error is 0.16 m/s. Compared with the traditional model, the navigation and positioning considering complex atmospheric changes are more accurate.


---

## Body

## 1. Introduction

Astronautical navigation is a positioning navigation method that uses the information of objects (moon, sun, other planets, and stars) measured by optical sensors to calculate the position of the carrier. Astronomical navigation and inertial navigation both belong to autonomous navigation technologies. Astronautical navigation is one of the widely used autonomous positioning technologies in the fields of space, aviation, and navigation. With the development of satellite application technology, the number of various spacecraft in space operation has increased dramatically. The amount of information transmission has increased dramatically too. Relying entirely on the ground station for measurement and control will cause the overload of the ground station, which will result in the information jam of the measurement and control system [1, 2]. With the development of manned spaceflight and deep space exploration technology, the requirement for spacecraft autonomous navigation capability has become more urgent. Astronomical navigation has become the most effective autonomous navigation method for spacecraft because of its high navigation accuracy, no error accumulation with time, strong anti-interference ability, and the ability to provide position and attitude information at the same time [3]. Autonomous astronavigation technology mainly includes direct and indirect sensitive horizons. The starlight refraction indirectly sensitive horizon method utilizes the optical properties of the atmosphere. For the spacecraft, when starlight passes through the Earth’s atmosphere, it refracts and bends toward the center of the Earth due to uneven atmospheric density, which will result in the star’s viewing position up from its actual position. The refraction angle depends on starlight frequency and atmospheric density [4, 5]. If the refraction angle of a known star near the horizon is measured, the direction of the Earth’s horizon in the geocentric inertial coordinate system can be obtained. The measurement equation will be established which is about the relationship between the apparent height of the refracted light relative to the Earth and the refraction angle, and then the position and velocity information of the spacecraft can be calculated [6].In the early 1906, before the Apollo program was completed, the Charles Stark Draper Laboratory (CSDL) conducted in-depth research on the use of navigation techniques such as starlight atmospheric refraction, starlight atmospheric attenuation, and other navigation technology methods [7, 8]. In the early 1980s, the French National Space Science Research Center (CNES) studied LEO satellites as navigation objects and carried out several experiments on the atmospheric refraction characteristics of the stratosphere and predicted that the precision of the navigation technology of star atmospheric refraction could reach 300 m [9]. Until the 1990s, when starlight atmospheric refraction navigation technology was put into practice in the MADAN (Multimission Attitude Determination and Autonomous Navigation) system in the United States, researchers found the effect of atmospheric density on navigation [10, 11]. In 2013, Ning et al. analyzed the atmospheric environment when they studied the method of autonomous astronomical navigation of the Earth’s direct sensitive horizon. They proposed that atmospheric density is the main factor affecting atmospheric refraction. At the same time, based on the atmospheric transmission characteristics, the stratospheric atmospheric refraction model was improved. However, due to the complex changes in the atmosphere, the accuracy of the model still needs to be improved [12]. In 2015, Ning et al. analyzed the impact of star sensors on navigation positioning accuracy and proposed a method to improve navigation accuracy using IUKF (iterated sampling Kalman filter) [13]. In 2017, Han et al. proposed a new data processing optimization method for star sensors, which improves the stability and accuracy of the star sensor’s attitude output [14]. In 2019, He et al. proposed a new stellar map centroid extraction algorithm to improve the accuracy and anti-interference performance of stellar centroid extraction in order to improve navigation accuracy. However, due to the refraction of stellar light, few researchers have studied the effect of atmospheric components on navigation [15]. An accurate atmospheric refraction model is required when applying starlight refraction theory to spacecraft navigation. Therefore, based on NCEP global atmospheric parameter data, this paper establishes an atmospheric parameter model considering time and space changes using the Fourier interpolation algorithm, studies the changes of atmospheric density with latitude, modifies the traditional atmospheric parameter model combined with the characteristics of star light transmission in the atmosphere, and uses unscented Kalman filter for indirectly sensitive horizon autonomous astronomical navigation of star light refraction.


## 2. Principle of Star Refraction Autonomous Navigation Method

It is observed from the spacecraft that the apparent height of the refracted light relative to the Earth isha, but in fact, it is at a slightly lower height hgfrom the ground, that is, the refraction height (as shown in Figure 1). Starlight atmospheric refraction angle is the angle between the two light rays. Star sensor captures the stars, establishes the inertial coordinate system through the geometric relationship between navigation objects to obtain the state parameters, takes the refraction apparent height as the observation of navigation model, and uses the filtering algorithm to achieve high-precision autonomous navigation.Figure 1
Principle diagram of atmospheric refraction of starlight.The accuracy of indirectly sensitive horizons for starlight refraction depends on the accuracy of atmospheric density. However, the actual changes in the Earth’s atmosphere are very complex and depend on many factors such as time, latitude, season, and day of solar activity. Most of the density data are calculated by the indirect measurement of temperature or pressure, so the accuracy is not high, which limits the density to be better estimated and modeled. The change of atmospheric density mainly depends on the change of temperature. In this paper, by creating a new temperature model, the actual change of atmospheric density can be analyzed accurately; thereby, the accuracy of navigation measurement equation will be improved.


### 2.1. Global Atmospheric Temperature and Density Model

#### 2.1.1. Atmospheric Temperature

Based on the NCEP FNL (final operational global analysis) global operation analysis data, the temperature data are analyzed and modeled. Given any time (day) and space (longitude and latitude), the atmospheric temperature model is expressed as(1)T=a0+∑i=1Naicosiwt+bisiniwt,where T is the atmospheric temperature, t is the time, ω is the time period of the model, a0, ai, and bi are the parameters to be determined. When N is 4 in the process of fitting the measured data, the model fits well with the true atmospheric change trend. Equation (1) can be interpolated for temperature variations with latitude and height.Using the model to fit temperature data for latitudes of 0 and 10 degrees, the error between the fitting and the measured temperature data is expressed by relative error1/M∑i=1Myi−y/y2, where yi is the value obtained by model fitting, yis the measured data, and M is the data volume. Taking the latitude range of 0∼10° as an example, when the latitude is 0°, the relative error between the fitting results and the measured results is 2.0%. When the latitude is 10°, the relative error between the fitting results and the measured results is 1.83%. The temperature variation curve with time in the range of 0∼10° is obtained by interpolation, as shown by the solid line in Figure 2(a). Figure 2(b)shows the variation of atmospheric temperature with time and altitude by interpolating altitude at fixed latitude based on the fitting results of Figure 2(a).Figure 2
Temperature model fitting curves. (a) Temperature variation with time and latitude at the same altitude; (b) temperature variation with time and altitude at the same latitude.
(a)(b)In order to verify the accuracy of model fitting, equation (1) is used to fit the temperature data from NCEP between 11.37 and13.06 km, interpolate different latitudes, and fit to obtain the change of temperature with time at different heights. When the interpolation height is 12 km, the change result of temperature with time is shown as the solid line in Figure 3.Figure 3
Variation of temperature with time at typical global latitudes. (a) Latitude is 5°; (b) latitude is 55°; (c) latitude is 95°; (d) latitude is 165°.
(a)(b)(c)(d)


#### 2.1.2. Atmospheric Density

The calculation model of atmospheric pressure with respect to altitude is(2)P=1.273×105e−0.1633h.The gas equation of state isρ=P/RT(R = 287, R is the gas constant). The atmospheric density can be obtained from the analysis of equations (1) and (2). The variation of atmospheric density with latitude at different altitudes is fitted by the model, as shown in Figure 4.Figure 4
Variation of atmospheric density with latitude at different altitudes. (a) Height range from 1 to 13 km; (b) height range from 13 km to 25 km; (c) height range from 25 km to 37 km; (d) height range from 37 km to 50 km.
(a)(b)(c)(d)As shown in Figure4, under the action of the Earth’s gravity, the atmosphere is dense to the lower layer, so the air density near the ground is the largest, and the density will continue to decrease with the increase of height. From the surface, the peak value of density with latitude decreases rapidly. Therefore, the model conforms to the vertical distribution law of atmospheric density.


### 2.2. Relationship between Apparent Height and Atmospheric Density and Refraction Angle

The apparent heightha can be used to confirm the relationship between the refraction angle R with the position of the spacecraft. The approximate mathematical relationship between haand Ris given in [16]:(3)haR,ρ=h0−HlnR+Hlnkλρ02πReH1/2+RHRe2π1/2,where kλ is the scattering parameter; it is only related to the wavelength of light λ. ρ0 is the density at the reference height h0. H is the height of the density scaler at height h0. Re is the mean radius of the Earth’s equator. The traditional model only considers the variation of atmospheric density with height, but the density is inhomogeneous. The newly established atmospheric density model can obtain the atmospheric density data reflecting the variation of latitude, time, and height, which makes the refraction model more accurate.From the geometric relations in Figure1, it can be seen that(4)ha=r2−u2+utanR−Re−a.In the above equation,u=rs⋅us, rs is the position vector of spacecraft; usis the direction vector of starlight before refraction; a is a very small amount, which can be ignored. Equations (3) and (4) establish the relationship between refraction measurement and spacecraft position, which is the key to apply starlight refraction to spacecraft autonomous astronomical navigation.


## 2.1. Global Atmospheric Temperature and Density Model

### 2.1.1. Atmospheric Temperature

Based on the NCEP FNL (final operational global analysis) global operation analysis data, the temperature data are analyzed and modeled. Given any time (day) and space (longitude and latitude), the atmospheric temperature model is expressed as(1)T=a0+∑i=1Naicosiwt+bisiniwt,where T is the atmospheric temperature, t is the time, ω is the time period of the model, a0, ai, and bi are the parameters to be determined. When N is 4 in the process of fitting the measured data, the model fits well with the true atmospheric change trend. Equation (1) can be interpolated for temperature variations with latitude and height.Using the model to fit temperature data for latitudes of 0 and 10 degrees, the error between the fitting and the measured temperature data is expressed by relative error1/M∑i=1Myi−y/y2, where yi is the value obtained by model fitting, yis the measured data, and M is the data volume. Taking the latitude range of 0∼10° as an example, when the latitude is 0°, the relative error between the fitting results and the measured results is 2.0%. When the latitude is 10°, the relative error between the fitting results and the measured results is 1.83%. The temperature variation curve with time in the range of 0∼10° is obtained by interpolation, as shown by the solid line in Figure 2(a). Figure 2(b)shows the variation of atmospheric temperature with time and altitude by interpolating altitude at fixed latitude based on the fitting results of Figure 2(a).Figure 2
Temperature model fitting curves. (a) Temperature variation with time and latitude at the same altitude; (b) temperature variation with time and altitude at the same latitude.
(a)(b)In order to verify the accuracy of model fitting, equation (1) is used to fit the temperature data from NCEP between 11.37 and13.06 km, interpolate different latitudes, and fit to obtain the change of temperature with time at different heights. When the interpolation height is 12 km, the change result of temperature with time is shown as the solid line in Figure 3.Figure 3
Variation of temperature with time at typical global latitudes. (a) Latitude is 5°; (b) latitude is 55°; (c) latitude is 95°; (d) latitude is 165°.
(a)(b)(c)(d)


### 2.1.2. Atmospheric Density

The calculation model of atmospheric pressure with respect to altitude is(2)P=1.273×105e−0.1633h.The gas equation of state isρ=P/RT(R = 287, R is the gas constant). The atmospheric density can be obtained from the analysis of equations (1) and (2). The variation of atmospheric density with latitude at different altitudes is fitted by the model, as shown in Figure 4.Figure 4
Variation of atmospheric density with latitude at different altitudes. (a) Height range from 1 to 13 km; (b) height range from 13 km to 25 km; (c) height range from 25 km to 37 km; (d) height range from 37 km to 50 km.
(a)(b)(c)(d)As shown in Figure4, under the action of the Earth’s gravity, the atmosphere is dense to the lower layer, so the air density near the ground is the largest, and the density will continue to decrease with the increase of height. From the surface, the peak value of density with latitude decreases rapidly. Therefore, the model conforms to the vertical distribution law of atmospheric density.


## 2.1.1. Atmospheric Temperature

Based on the NCEP FNL (final operational global analysis) global operation analysis data, the temperature data are analyzed and modeled. Given any time (day) and space (longitude and latitude), the atmospheric temperature model is expressed as(1)T=a0+∑i=1Naicosiwt+bisiniwt,where T is the atmospheric temperature, t is the time, ω is the time period of the model, a0, ai, and bi are the parameters to be determined. When N is 4 in the process of fitting the measured data, the model fits well with the true atmospheric change trend. Equation (1) can be interpolated for temperature variations with latitude and height.Using the model to fit temperature data for latitudes of 0 and 10 degrees, the error between the fitting and the measured temperature data is expressed by relative error1/M∑i=1Myi−y/y2, where yi is the value obtained by model fitting, yis the measured data, and M is the data volume. Taking the latitude range of 0∼10° as an example, when the latitude is 0°, the relative error between the fitting results and the measured results is 2.0%. When the latitude is 10°, the relative error between the fitting results and the measured results is 1.83%. The temperature variation curve with time in the range of 0∼10° is obtained by interpolation, as shown by the solid line in Figure 2(a). Figure 2(b)shows the variation of atmospheric temperature with time and altitude by interpolating altitude at fixed latitude based on the fitting results of Figure 2(a).Figure 2
Temperature model fitting curves. (a) Temperature variation with time and latitude at the same altitude; (b) temperature variation with time and altitude at the same latitude.
(a)(b)In order to verify the accuracy of model fitting, equation (1) is used to fit the temperature data from NCEP between 11.37 and13.06 km, interpolate different latitudes, and fit to obtain the change of temperature with time at different heights. When the interpolation height is 12 km, the change result of temperature with time is shown as the solid line in Figure 3.Figure 3
Variation of temperature with time at typical global latitudes. (a) Latitude is 5°; (b) latitude is 55°; (c) latitude is 95°; (d) latitude is 165°.
(a)(b)(c)(d)


## 2.1.2. Atmospheric Density

The calculation model of atmospheric pressure with respect to altitude is(2)P=1.273×105e−0.1633h.The gas equation of state isρ=P/RT(R = 287, R is the gas constant). The atmospheric density can be obtained from the analysis of equations (1) and (2). The variation of atmospheric density with latitude at different altitudes is fitted by the model, as shown in Figure 4.Figure 4
Variation of atmospheric density with latitude at different altitudes. (a) Height range from 1 to 13 km; (b) height range from 13 km to 25 km; (c) height range from 25 km to 37 km; (d) height range from 37 km to 50 km.
(a)(b)(c)(d)As shown in Figure4, under the action of the Earth’s gravity, the atmosphere is dense to the lower layer, so the air density near the ground is the largest, and the density will continue to decrease with the increase of height. From the surface, the peak value of density with latitude decreases rapidly. Therefore, the model conforms to the vertical distribution law of atmospheric density.


## 2.2. Relationship between Apparent Height and Atmospheric Density and Refraction Angle

The apparent heightha can be used to confirm the relationship between the refraction angle R with the position of the spacecraft. The approximate mathematical relationship between haand Ris given in [16]:(3)haR,ρ=h0−HlnR+Hlnkλρ02πReH1/2+RHRe2π1/2,where kλ is the scattering parameter; it is only related to the wavelength of light λ. ρ0 is the density at the reference height h0. H is the height of the density scaler at height h0. Re is the mean radius of the Earth’s equator. The traditional model only considers the variation of atmospheric density with height, but the density is inhomogeneous. The newly established atmospheric density model can obtain the atmospheric density data reflecting the variation of latitude, time, and height, which makes the refraction model more accurate.From the geometric relations in Figure1, it can be seen that(4)ha=r2−u2+utanR−Re−a.In the above equation,u=rs⋅us, rs is the position vector of spacecraft; usis the direction vector of starlight before refraction; a is a very small amount, which can be ignored. Equations (3) and (4) establish the relationship between refraction measurement and spacecraft position, which is the key to apply starlight refraction to spacecraft autonomous astronomical navigation.


## 3. Navigation System Simulation Verification

### 3.1. Unscented Kalman Filter

Assume a discrete nonlinear system(5)xk+1=Fxk,uk,k+ωk,(6)zk=Hxk,uk,k+vk,where xk is the system state vector. uk is the input control vector. wk is the system noise vector. zk is the observation vector. vk is the measurement noise vector. Here, a series of sampling points around x^k are selected whose mean and covariance are x^k and Pk. A series of sampling points around x^k|k are selected whose mean and covariance are x^k|k and Pk|k. These sampling points generate corresponding transform sampling points through the nonlinear system. Let the state variable be N×1 dimension, then 2n+1 sampling points and their weights are as follows:(7)χ0,k=x^kW0=τn+τχi,k=x^k+n+τPk|kiWi=12n+τχi+n,k=x^k−n+τPk|kiWi+n=12n+τ.In the above equation,τ∈R, when Pk|k=ATA, Pk|ki takes line i of A.When Pk|k=AAT, Pk|ki takes column i of A. The standard unscented Kalman filter algorithm is as follows:(1)
Initializationk≥1:(8)x^0=Ex0,(9)P0=Ex0−x^0x0−x^0T.(2)
Calculating sampling points:(10)χk−1=x^k−1x^k−1+n+τPk−1ix^k−1−n+τPk−1i,i=1,2,…n.(3)
Time update:(11)χk|k−1=Fχk−1,uk−1,k−1,(12)x^k¯=∑i=02nWiχi,k|k−1,(13)Pk¯=∑Wiχi,k|k−1−x^k¯χi,k|k−1−x^k¯T+Qk,(14)yk|k−1=Hχk|k−1,uk,k,(15)y^k¯=∑i=02nWiyi,k|k−1.(4)
Measurement update:(16)Py^ky^k=∑i=02nWiyi,k|k−1−y^k¯yi,k|k−1−y^k¯T+Rk,(17)Pxkyk=∑i=02nWiχi,k|k−1−x¯k¯yi,k|k−1−y^k¯T,(18)Kk=PxkykPy^ky^k−1,(19)x^k=x^k¯+Kkyk−y^k¯,(20)Pk=Pk¯−KkPy^ky^kKkT.Qk and Rk are the covariance of system noise and measurement noise, respectively; when xkis assumed to be a Gaussian distribution, n+τ is equal to 3.


### 3.2. State Equation and Measurement Equation of Navigation System

When establishing the system equation of state, only the second-order harmonic terms perturbed by the gravity of the Earth’s mass center and the gravitational field are considered, while other perturbation factors are equivalent to Gaussian white noise. The state model of the satellite navigation system in epoch (J2000.0) geocentric equatorial coordinate system is(21)dxdt=vxdydt=vydzdt=vzdvxdt=−μxr31−J2Rer7.5z2r2−1.5+ΔFxdvydt=−μyr31−J2Rer7.5z2r2−1.5+ΔFydvzdt=−μzr31−J2Rer7.5z2r2−1.5+ΔFz,where r=x2+y2+z2; then, the equation of state can be expressed as(22)X˙t=fX,t+wt.In equations (21) and (22), the state vectors X=xyzvxvyvzT,x,y,z,vx,vy, and vzare the position and velocity of the satellite in the X, Y, and Zdirections, respectively. μis the gravitational constant of the Earth, μ=3.986×1014m3/s2. ris the size of the satellite position vector. J2is the gravitational coefficient of the Earth, J2=1.08263×10−3. ΔFx,ΔFy,ΔFzare equivalent to Gaussian white noise.In the astronomical navigation of starlight refraction indirectly sensitive horizon, the refraction apparent heighthais usually selected as the observation measurement, but there is no direct relationship between the apparent height haand the refractive angle R. The apparent height needs to be obtained indirectly through the relationship between the refractive height hgand R. Therefore, equation (4) is transformed into the measurement equation by adding the Gaussian measurement noise v. Equation (23) is combined with the orbital dynamics equation to calculate the position information of the aircraft:(23)ha=rs2−u2+utanR−Re−a+v.


### 3.3. Simulation Result

Computer simulation is carried out based on the state and measurement equations of an autonomous astronomical navigation system with indirectly sensitive horizon for Earth satellites. Assuming the initial orbital root number of the satellite is shown in Table1.Table 1
Number of satellite initial orbits.
Semimajor axis7136.635 kmEccentricity0.001809Inclination65°Right ascension of the ascending node30.00°Argument of periapsis30.00°Time past perigee0After establishing the system state equation and measurement equation, according to the noise characteristics,n=6,τ=−3are selected. The unscented Kalman filter is used for the computer simulation of autonomous navigation with indirect sensitive horizon. The noise variance matrix of discrete system is(24)q1=1×10−3m/s2,q2=2×10−3m/s2,Q=diagq1,q2,q1,q1,q2,q1.Initial state value is(25)X0=4.590×106m,4.388×106m,3.228×106m,−4.612×103m/s,5.014×103m/s,5.876×103m/sT.Initial estimation error matrix is(26)P0=diag600m2,600m2,500m2,2m/s2,1.6m/s2,1.5m/s2.Assuming the starlight wavelength to be 550 nm, atmospheric density model error is 1%1σ, and the measured noise variance R=EvvT=80m2. Assuming that the system noise is not related to the measured noise, the sampling period of the simulation system is 3 seconds. The simulation results of position and speed according to the above simulation conditions are shown in Figure 5.Figure 5
UKF simulation results. (a) Position estimation error; (b) velocity estimation error.
(a)(b)Considering the variation of density with latitude and altitude, the simulation results are shown in Figure5. The position estimation error of unscented Kalman filter indirect sensitive horizon navigation system is 94 m, and the velocity estimation error is 0.16 m/s. It can be seen from Figure 5that the navigation system with improved atmospheric density model has stability. Because the traditional Kalman filter can only be applied to linear systems, the estimation accuracy is reduced. However, the state equation and measurement equation are nonlinear. Unscented Kalman filter can filter the nonlinear system directly without calculating the Jacobian matrix of state equation and measurement equation, so it is higher than Kalman filter in estimation accuracy.


## 3.1. Unscented Kalman Filter

Assume a discrete nonlinear system(5)xk+1=Fxk,uk,k+ωk,(6)zk=Hxk,uk,k+vk,where xk is the system state vector. uk is the input control vector. wk is the system noise vector. zk is the observation vector. vk is the measurement noise vector. Here, a series of sampling points around x^k are selected whose mean and covariance are x^k and Pk. A series of sampling points around x^k|k are selected whose mean and covariance are x^k|k and Pk|k. These sampling points generate corresponding transform sampling points through the nonlinear system. Let the state variable be N×1 dimension, then 2n+1 sampling points and their weights are as follows:(7)χ0,k=x^kW0=τn+τχi,k=x^k+n+τPk|kiWi=12n+τχi+n,k=x^k−n+τPk|kiWi+n=12n+τ.In the above equation,τ∈R, when Pk|k=ATA, Pk|ki takes line i of A.When Pk|k=AAT, Pk|ki takes column i of A. The standard unscented Kalman filter algorithm is as follows:(1)
Initializationk≥1:(8)x^0=Ex0,(9)P0=Ex0−x^0x0−x^0T.(2)
Calculating sampling points:(10)χk−1=x^k−1x^k−1+n+τPk−1ix^k−1−n+τPk−1i,i=1,2,…n.(3)
Time update:(11)χk|k−1=Fχk−1,uk−1,k−1,(12)x^k¯=∑i=02nWiχi,k|k−1,(13)Pk¯=∑Wiχi,k|k−1−x^k¯χi,k|k−1−x^k¯T+Qk,(14)yk|k−1=Hχk|k−1,uk,k,(15)y^k¯=∑i=02nWiyi,k|k−1.(4)
Measurement update:(16)Py^ky^k=∑i=02nWiyi,k|k−1−y^k¯yi,k|k−1−y^k¯T+Rk,(17)Pxkyk=∑i=02nWiχi,k|k−1−x¯k¯yi,k|k−1−y^k¯T,(18)Kk=PxkykPy^ky^k−1,(19)x^k=x^k¯+Kkyk−y^k¯,(20)Pk=Pk¯−KkPy^ky^kKkT.Qk and Rk are the covariance of system noise and measurement noise, respectively; when xkis assumed to be a Gaussian distribution, n+τ is equal to 3.


## 3.2. State Equation and Measurement Equation of Navigation System

When establishing the system equation of state, only the second-order harmonic terms perturbed by the gravity of the Earth’s mass center and the gravitational field are considered, while other perturbation factors are equivalent to Gaussian white noise. The state model of the satellite navigation system in epoch (J2000.0) geocentric equatorial coordinate system is(21)dxdt=vxdydt=vydzdt=vzdvxdt=−μxr31−J2Rer7.5z2r2−1.5+ΔFxdvydt=−μyr31−J2Rer7.5z2r2−1.5+ΔFydvzdt=−μzr31−J2Rer7.5z2r2−1.5+ΔFz,where r=x2+y2+z2; then, the equation of state can be expressed as(22)X˙t=fX,t+wt.In equations (21) and (22), the state vectors X=xyzvxvyvzT,x,y,z,vx,vy, and vzare the position and velocity of the satellite in the X, Y, and Zdirections, respectively. μis the gravitational constant of the Earth, μ=3.986×1014m3/s2. ris the size of the satellite position vector. J2is the gravitational coefficient of the Earth, J2=1.08263×10−3. ΔFx,ΔFy,ΔFzare equivalent to Gaussian white noise.In the astronomical navigation of starlight refraction indirectly sensitive horizon, the refraction apparent heighthais usually selected as the observation measurement, but there is no direct relationship between the apparent height haand the refractive angle R. The apparent height needs to be obtained indirectly through the relationship between the refractive height hgand R. Therefore, equation (4) is transformed into the measurement equation by adding the Gaussian measurement noise v. Equation (23) is combined with the orbital dynamics equation to calculate the position information of the aircraft:(23)ha=rs2−u2+utanR−Re−a+v.


## 3.3. Simulation Result

Computer simulation is carried out based on the state and measurement equations of an autonomous astronomical navigation system with indirectly sensitive horizon for Earth satellites. Assuming the initial orbital root number of the satellite is shown in Table1.Table 1
Number of satellite initial orbits.
Semimajor axis7136.635 kmEccentricity0.001809Inclination65°Right ascension of the ascending node30.00°Argument of periapsis30.00°Time past perigee0After establishing the system state equation and measurement equation, according to the noise characteristics,n=6,τ=−3are selected. The unscented Kalman filter is used for the computer simulation of autonomous navigation with indirect sensitive horizon. The noise variance matrix of discrete system is(24)q1=1×10−3m/s2,q2=2×10−3m/s2,Q=diagq1,q2,q1,q1,q2,q1.Initial state value is(25)X0=4.590×106m,4.388×106m,3.228×106m,−4.612×103m/s,5.014×103m/s,5.876×103m/sT.Initial estimation error matrix is(26)P0=diag600m2,600m2,500m2,2m/s2,1.6m/s2,1.5m/s2.Assuming the starlight wavelength to be 550 nm, atmospheric density model error is 1%1σ, and the measured noise variance R=EvvT=80m2. Assuming that the system noise is not related to the measured noise, the sampling period of the simulation system is 3 seconds. The simulation results of position and speed according to the above simulation conditions are shown in Figure 5.Figure 5
UKF simulation results. (a) Position estimation error; (b) velocity estimation error.
(a)(b)Considering the variation of density with latitude and altitude, the simulation results are shown in Figure5. The position estimation error of unscented Kalman filter indirect sensitive horizon navigation system is 94 m, and the velocity estimation error is 0.16 m/s. It can be seen from Figure 5that the navigation system with improved atmospheric density model has stability. Because the traditional Kalman filter can only be applied to linear systems, the estimation accuracy is reduced. However, the state equation and measurement equation are nonlinear. Unscented Kalman filter can filter the nonlinear system directly without calculating the Jacobian matrix of state equation and measurement equation, so it is higher than Kalman filter in estimation accuracy.


## 4. Conclusion

Aiming at the inaccuracy of the traditional atmospheric parameter model, based on NCEP atmospheric parameter data, this paper establishes a spatiotemporal atmospheric temperature model by using the Fourier interpolation algorithm. By using the new atmospheric temperature model, the variation of atmospheric density with latitude is analyzed, and the relationship between apparent altitude and atmospheric is given. Then, the spacecraft state equation and the measurement equation taking the starlight apparent height as the observation are established, and the Unscented Kalman filter is applied to the starlight refraction indirectly sensitive horizon autonomous astronomical navigation. The relative error of fitting the measured data with the spatiotemporal variation atmospheric temperature model established is less than 2%. The result of filtering is that the position estimation error is 94 m, and the velocity estimation error is 0.16 m/s. Compared with the traditional model, the navigation and positioning considering the variation of density with altitude and latitude are more accurate.



---

*Source: 3649662-2022-02-03.xml*