# An Improved Orbit Model of Space Target Radar Based on Least Square Method and Its Program Implementation

Authors: Yan Zhang; Xiao-Bin Huang; Rui Xiao; Ling Wang; Li Lu
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

Space target orbit determination module is an important component of the space target surveillance radar system. The development of this module requires very complex aerospace dynamics knowledge, which brings great difficulties to nonorbit mechanics researchers engaged in radar system design. Orbit improvement is a core content of orbit determination, and it is a necessary step to achieve high-precision orbit calculations. To this end, this paper focuses on the issue of batch processing orbit improvement of space target surveillance radar, introduces the principle of least-squares orbit improvement, the partial derivative of the model, and the main perturbation acceleration calculation methods, and gives the program design principle and implements the RadarOrbDet library. The developed library is compared and analyzed with STK and ODTK software, and the simulation results verify the effectiveness of the library. The library is also helpful for designers of space target surveillance radar systems to carry out a rapid demonstration of orbit determination indicators.


---

## Body

## 1. Introduction

The space target orbit determination module is an important component of the space target surveillance radar system. The development of this module requires a very complex knowledge of aerospace dynamics [1], which brings great difficulties to nonorbit mechanics researchers engaged in radar system design. To this end, this paper analyzes the basic orbit determination principles of space target radar and develops the RadarOrbDet open source library [2] to assist nonorbit mechanics engineers in the system design of space target surveillance radar.According to different data processing methods, the orbit determination methods can be divided into dynamic orbit determination, geometric orbit determination, and reduced dynamic orbit determination. The dynamic orbit determination method uses a dynamic model of space target to establish motion equation and determines orbit based on the constraints of the dynamic equation. It can obtain a better orbit determination effect with fewer observation data, and it is widely used in navigation satellite orbit determination [3]. Geometric orbit determination method does not consider the orbital dynamics model and obtains the orbital elements by fitting the observation data. This kind of method has been widely used in the orbit determination of low-orbit satellites [4], but the orbit extrapolation accuracy is relatively low. The reduced dynamic orbit determination method is based on orbital mechanics and combines geometric observation information for optimal weighting so as to achieve precise orbit determination. This method has also been widely used in satellite orbit determination and effectively alleviates the problem that the dynamic model is sensitive to observation errors [5, 6].Modern orbit determination technology requires the measurement model to be expressed as a function of the orbit state. Usually, the measured value is a nonlinear function of the orbital position component, and then the weighted least-squares technique is used to minimize the residuals to solve the target state vector. This solving process requires repeated iterative corrections to the initial orbital state, so it is also called orbital improvement. Various types of measurement data can be used for orbit improvement, such as angle measurement data [7, 8], GPS measurement data [9, 10], doppler measurement data [11, 12], etc. The calculation process involves various perturbation force analyses, such as the Earth’s nonspherical gravitational perturbation [13, 14], atmospheric drag perturbation [15, 16], third-body perturbation [17, 18], solar radiation perturbation [19, 20], and so on. In order to speed up data analysis tasks, many open source orbit determination software [21–25] can be used by engineers to void complicated mathematical formula derivation, but these software is mainly used to process optical angle measurement data and are not suitable for processing radar measurement data.This paper focuses on the batch orbit improvement in space target radar orbit determination, introduces the principle of least-squares orbit improvement, explains the calculation method of partial derivative and main perturbations in the model, describes the program design of the RadarOrbDet library, and finally compares the library with System Tool Kit (STK) and Orbit Determination Tool Kit (ODTK) software. The simulations verify the effectiveness of the RadarOrbDet library, which shows that the library is helpful for designers of space target surveillance radar systems to carry out a rapid demonstration of orbit determination indicators.


## 2. The Principle of Least-Squares Orbit Improvement

The basic principle of least-squares orbit determination is to find an orbit that minimizes the residual between theoretical measurements and actual measurements, which is, solvingx0 to minimize the value of the following equation:(1)Jx0=z−hx0Tz−hx0,where x0=rt0;vt0 is the state of a satellite at epoch t0, and a reference trajectory can be predicted by this orbit state; z=z1;z2⋯;zN is the actual measurement vector of radar, and each measurement contains three components, namely distance, azimuth, and elevation; hx0 is the calculated measurement obtained with reference trajectory.Sinceh is a nonlinear function, equation (1) describes a nonlinear least-squares problem, which is extremely difficult to solve. Usually, Taylor expansion is used to transform this nonlinear least-squares problem into a linear least-squares problem, and its solution can be approximately obtained using the following iterative equation:(2)x0i+1=x0i+HiTHi−1HiTz−hx0i.The initial value of the iteration isx00=x0apr, which can be obtained by initial orbit determination. The iterative process terminates until the relative change of the error for two consecutive times is less than a given threshold. H is Jacobian matrix described by the following equation:(3)Hi=∂hx0∂x0|x0=x0i.Taking into account the weight of different types of measurement, equation (2) can be modified as follows:(4)x0i+1=x0i+HiTWHi−1HiTWz−hx0i,where W=diagσρ−2,σα−2,σβ−2 is a diagonal square matrix composed of radar measurement errors.


## 3. Computation of Model’s Partial Derivatives

The key to the calculation of equation (4) lies in the computation of the matrix H. Without loss of generality, take the once radar measurement as an example and omit the superscript i.Using the derivative chain rule,H is transformed into the following:(5)H=∂h∂x⋅∂x∂x0=ΑΦ,where matrix A is the partial derivative matrix of the current measurement with respect to the current state, and Φ is the partial derivative matrix of the current state with respect to the initial state, also called the error state transition matrix.


### 3.1. Computation of Partial Derivative Matrix of Measurement with respect to the State Vector

Ignoring the optical aberration and other minor factors,, the radar measurement is only related to the current position of the satellite. Therefore, matrixA can be split into the following:(6)Α=∂h∂x=∂h∂rECI03×3,where rECI is the position vector of the satellite in the inertial coordinate system (ECI coordinate system).Radar measurement is based on the topocentric horizon coordinate system, such as SEZ (south-east- zenith) coordinate. It is difficult to directly calculate the partial derivative of the radar measurement with respect to the position vector in the ECI coordinate system. Therefore, the calculation of matrixA needs to be further decomposed by derivation chain rule(7)∂h∂rECI=∂h∂ρSEZ⋅∂ρSEZ∂ρECEF⋅∂ρECEF∂rECEF⋅∂rECEF∂rECI,where ρSEZ is the Cartesian position vector of the satellite in the SEZ coordinate system, represented by ρS,ρE,ρZ; ρECEF is a range vector from radar to satellite represented in the Earth-Fixed coordinate system (ECEF); rECEF is the satellite’s position vector in ECEF coordinate system.Notice(8)ρECEF=rECEF−rsiteECEF,where rsiteECEF is the radar position vector in ECEF coordinate system, which is a constant vector.Therefore,(9)∂ρECEF∂rECEF=∂rECEF−rsiteECEF∂rECEF=I.According to the coordinate transformation(10)ρSEZ=SEZECEFρECEF,where SEZ/ECEF is the transformation matrix from SEZ coordinate system to ECEF coordinate system.According to equation (10),(11)∂ρSEZ∂ρECEF=∂∂ρECEFSEZECEFρECEF=∂∂ρECEFSEZECEFρECEF+SEZECEF∂ρECEF∂ρECEF=SEZECEF.It can be seen from equation (11) that the partial derivative matrix is actually a coordinate transformation matrix.Similarly,(12)∂rECEF∂rECI=ECEFECI.Substituting equations (9), (11), and (12) into equation (7), we have the following:(13)∂h∂rECI=∂h∂ρSEZ⋅SEZECEF⋅ECEFECI=∂h∂ρSEZ⋅SEZECI,where SEZ/ECI is the transformation matrix from SEZ coordinate system to ECI coordinate system.In SEZ coordinate system, the relationship between Cartesian coordinate and polar coordinate is as follows:(14)ρ=ρS2+ρE2+ρZ2,α=atan−ρSρE,β=asinρZρS2+ρE2+ρZ2.So,(15)∂h∂ρSEZ=∂ρ∂ρS∂ρ∂ρE∂ρ∂ρZ∂α∂ρS∂α∂ρE∂α∂ρZ∂β∂ρS∂β∂ρE∂β∂ρZ=ρSρρEρρZρρEρS2+ρE21ρS2+ρE20−ρSρZρ2ρS2+ρE2−ρEρZρ2ρS2+ρE2ρS2+ρE2ρ2ρS2+ρE2.Substituting equation (15) into equation (13), we have the following:(16)∂h∂rECI=ρSρρEρρZρρEρS2+ρE21ρS2+ρE20−ρSρZρ2ρS2+ρE2−ρEρZρ2ρS2+ρE2ρS2+ρE2ρ2ρS2+ρE2⋅SEZECI.The calculation of the two matrices in equation (16) needs to be performed at the current observation time. Finally, substitute equation (16) into equation (6) to obtain matrix A.


### 3.2. Computation of Error State Transition Matrix

Suppose the satellite’s motion equation is as follows:(17)dxtdt=x˙=ft,x=vtat,r,v.Then,(18)Φ˙=ddt∂xt∂xt0=∂∂xt0dxtdt=∂ft,xt∂xt0=∂ft,xt∂xt⋅∂xt∂xt0=FΦ,where the matrix F is as follows:(19)F=03×313×3∂at,r,v∂r∂at,r,v∂v.Due to the iterative mode being used to solve the least-squares problem, accuracy requirements for the partial derivativesΦ and F are generally more relaxed than those for the trajectory itself. It is common to apply a simplified force model in the solution of the equation (18). The incorporation of the lowest-order zonal gravity field perturbation (C2,0) already provides an acceptable minimum model.When only the Earth’s gravity is considered, the acceleration is only related to the position, and its calculation is usually carried out in the ECEF coordinate system. The specific expression can be found in reference [26].Finally, the partial derivative matrix in the ECEF coordinate system needs to be converted to the ECI coordinate system. In the case of ignoring Coriolis force and centrifugal force, there is the following relationship:(20)∂a∂rECI=ECIECEF−1∂a∂rECEFECIECEF.At the current time, after the matrixF is calculated, the numerical differential equation method is used to solve matrix Φ.


## 3.1. Computation of Partial Derivative Matrix of Measurement with respect to the State Vector

Ignoring the optical aberration and other minor factors,, the radar measurement is only related to the current position of the satellite. Therefore, matrixA can be split into the following:(6)Α=∂h∂x=∂h∂rECI03×3,where rECI is the position vector of the satellite in the inertial coordinate system (ECI coordinate system).Radar measurement is based on the topocentric horizon coordinate system, such as SEZ (south-east- zenith) coordinate. It is difficult to directly calculate the partial derivative of the radar measurement with respect to the position vector in the ECI coordinate system. Therefore, the calculation of matrixA needs to be further decomposed by derivation chain rule(7)∂h∂rECI=∂h∂ρSEZ⋅∂ρSEZ∂ρECEF⋅∂ρECEF∂rECEF⋅∂rECEF∂rECI,where ρSEZ is the Cartesian position vector of the satellite in the SEZ coordinate system, represented by ρS,ρE,ρZ; ρECEF is a range vector from radar to satellite represented in the Earth-Fixed coordinate system (ECEF); rECEF is the satellite’s position vector in ECEF coordinate system.Notice(8)ρECEF=rECEF−rsiteECEF,where rsiteECEF is the radar position vector in ECEF coordinate system, which is a constant vector.Therefore,(9)∂ρECEF∂rECEF=∂rECEF−rsiteECEF∂rECEF=I.According to the coordinate transformation(10)ρSEZ=SEZECEFρECEF,where SEZ/ECEF is the transformation matrix from SEZ coordinate system to ECEF coordinate system.According to equation (10),(11)∂ρSEZ∂ρECEF=∂∂ρECEFSEZECEFρECEF=∂∂ρECEFSEZECEFρECEF+SEZECEF∂ρECEF∂ρECEF=SEZECEF.It can be seen from equation (11) that the partial derivative matrix is actually a coordinate transformation matrix.Similarly,(12)∂rECEF∂rECI=ECEFECI.Substituting equations (9), (11), and (12) into equation (7), we have the following:(13)∂h∂rECI=∂h∂ρSEZ⋅SEZECEF⋅ECEFECI=∂h∂ρSEZ⋅SEZECI,where SEZ/ECI is the transformation matrix from SEZ coordinate system to ECI coordinate system.In SEZ coordinate system, the relationship between Cartesian coordinate and polar coordinate is as follows:(14)ρ=ρS2+ρE2+ρZ2,α=atan−ρSρE,β=asinρZρS2+ρE2+ρZ2.So,(15)∂h∂ρSEZ=∂ρ∂ρS∂ρ∂ρE∂ρ∂ρZ∂α∂ρS∂α∂ρE∂α∂ρZ∂β∂ρS∂β∂ρE∂β∂ρZ=ρSρρEρρZρρEρS2+ρE21ρS2+ρE20−ρSρZρ2ρS2+ρE2−ρEρZρ2ρS2+ρE2ρS2+ρE2ρ2ρS2+ρE2.Substituting equation (15) into equation (13), we have the following:(16)∂h∂rECI=ρSρρEρρZρρEρS2+ρE21ρS2+ρE20−ρSρZρ2ρS2+ρE2−ρEρZρ2ρS2+ρE2ρS2+ρE2ρ2ρS2+ρE2⋅SEZECI.The calculation of the two matrices in equation (16) needs to be performed at the current observation time. Finally, substitute equation (16) into equation (6) to obtain matrix A.


## 3.2. Computation of Error State Transition Matrix

Suppose the satellite’s motion equation is as follows:(17)dxtdt=x˙=ft,x=vtat,r,v.Then,(18)Φ˙=ddt∂xt∂xt0=∂∂xt0dxtdt=∂ft,xt∂xt0=∂ft,xt∂xt⋅∂xt∂xt0=FΦ,where the matrix F is as follows:(19)F=03×313×3∂at,r,v∂r∂at,r,v∂v.Due to the iterative mode being used to solve the least-squares problem, accuracy requirements for the partial derivativesΦ and F are generally more relaxed than those for the trajectory itself. It is common to apply a simplified force model in the solution of the equation (18). The incorporation of the lowest-order zonal gravity field perturbation (C2,0) already provides an acceptable minimum model.When only the Earth’s gravity is considered, the acceleration is only related to the position, and its calculation is usually carried out in the ECEF coordinate system. The specific expression can be found in reference [26].Finally, the partial derivative matrix in the ECEF coordinate system needs to be converted to the ECI coordinate system. In the case of ignoring Coriolis force and centrifugal force, there is the following relationship:(20)∂a∂rECI=ECIECEF−1∂a∂rECEFECIECEF.At the current time, after the matrixF is calculated, the numerical differential equation method is used to solve matrix Φ.


## 4. Calculation of Main Perturbation Acceleration

In solving least-squares, a reference trajectory of the satellite needs to be calculated. So, it is necessary to appropriately select the perturbation force of the satellite according to the problem and the accuracy requirement. For LEO satellites, the Earth’s nonspherical gravity and atmospheric drag perturbation are usually considered. The following describes the satellite acceleration caused by these two perturbation forces.


### 4.1. Gravitational Acceleration of the Earth

UserECI and rECEF to denote the position vectors of the satellite under the ECI and ECEF coordinate system, respectively. The acceleration can be obtained by calculating the gradient of the potential function U, namely,(21)aECI=∇UrECEF=∂rECI∂rECEF∂U∂rECEF∂rECEF∂rECEFT+∂U∂φ∂φ∂rECEFT+∂U∂λ∂λ∂rECEFT=ECIECEFaECEF.where aECI and aECEF are the accelerations of the satellite in the ECI and ECEF coordinates, respectively, rECEF=xECEF2+yECEF2+zECEF2, (rECEF, φ, λ) is the altitude, latitude, and longitude of the satellite in the ECEF coordinate system; ECI/ECEF is the coordinate transformation matrix from ECI to ECEF coordinate system.The three components ofaECEF can be calculated by the following equation:(22)axECEF=xECEFrECEF∂U∂rECEF−xECEFzECEFrECEF2xECEF2+yECEF2∂U∂φ−yECEFxECEF2+yECEF2∂U∂λ,ayECEF=yECEFrECEF∂U∂rECEF−ybzrECEF2xECEF2+yECEF2∂U∂φ+xbxECEF2+yECEF2∂U∂λ,azECEF=zECEFrECEF∂U∂rECEF+xECEF2+yECEF2rECEF2∂U∂φ.The partial derivative of the potential functionU to (rECEF, φ, λ) is as follows:(23)∂U∂rECEF=−GMerECEF2∑n=0∞∑m=0nn+1RerECEFnP¯nmsinφC¯nmcosmλ+S¯nmsinmλ,∂U∂φ=GMerECEF2∑n=0∞∑m=0nRerECEFnPnm+1sinφ−mtanφPnmsinφC¯nmcosmλ+S¯nmsinmλ,∂U∂λ=GMerECEF2∑n=0∞∑m=0nRerECEFnmP¯nmsinφS¯nmcosmλ−C¯nmsinmλ,where G is the gravitational constant, Me is the mass of the Earth, Re is the average radius of the Earth’s equator, P¯nmx is the normalized Legendre polynomial, and C¯nm and S¯nm are the normalized gravitational potential coefficients, which can be read out directly from the Earth’s gravity field model file.


### 4.2. Atmospheric Drag Acceleration

The satellite acceleration caused by atmospheric drag can be written as follows:(24)adrag=−12CDAmρvr2ev,where m is the satellite mass, CD is the atmospheric drag coefficient, and its typical value ranges from 1.5 to 3.0. The relative velocity vr can be written approximately as follows:(25)vr=v−ω⊕×r,with the inertial satellite velocity vector v, the position vector r, and the Earth’s angular velocity vector ω⊕ of size 0.7292 × 10−4 rad/s. Note that the calculation of equations (24) and (25) is carried out in the true-of-date (TOD) coordinate system, and the corresponding value in the ECI coordinate system can be obtained after a coordinate transformation.The drag depends on the atmospheric densityρ at the location of the satellite, and the atmospheric density is usually calculated by the semiempirical atmospheric density model. At present, the commonly used atmospheric model in orbit determination is NRLMSISE-00 atmospheric model, which plays an important role in satellite orbit determination and prediction.The NRLMSISE-00 atmospheric model has 8 input items: The number of days from January 1 of the current year to the current day, the number of seconds from 00:00:00 to the calculating time, geographical longitude, latitude, altitude, the solar radiation flux of 10.7 cm on the previous day (F10.7), the average F10.7 for 81 days (3 solar rotation cycles, with the current day as the midpoint), average geomagnetic index (Ap) of current day and the twenty 3 h average values of Ap before the calculating time. The outputs include the number densities of N2, O2, He, Ar, N, H, O, and O+, the neutral atmospheric temperature, and the atmospheric density.These space environment parameters are provided by the SpaceWeather.txt file. Each line in the file represents a spatial environment parameter record for the specified date. The record has 32 fields, as shown in Figure1. The meanings of the fields related to atmospheric drag are shown in Table 1.Figure 1
Parameter format of spatial environment parameter record.Table 1
The meaning of space environment parameters.
ColumnNameDescription001–004yyyyYear006–007mmMonth009–010ddDay048–050Ap0Ap index (00:00–03:00)052–054Ap3Ap index (03:00–06:00)056–058Ap6Ap index (06:00–09:00)060–062Ap9Ap index (09:00–12:00)064–066Ap12Ap index (12:00–15:00)068–070Ap15Ap index (15:00–18:00)072–074Ap18Ap index (18:00–21:00)076–078Ap21Ap index (21:00–24:00)114–118ObsF10.7F10.7, previous day (observed)126–130Obs Lst81F10.7, arithmetic mean for the previous 81 days (observed)


## 4.1. Gravitational Acceleration of the Earth

UserECI and rECEF to denote the position vectors of the satellite under the ECI and ECEF coordinate system, respectively. The acceleration can be obtained by calculating the gradient of the potential function U, namely,(21)aECI=∇UrECEF=∂rECI∂rECEF∂U∂rECEF∂rECEF∂rECEFT+∂U∂φ∂φ∂rECEFT+∂U∂λ∂λ∂rECEFT=ECIECEFaECEF.where aECI and aECEF are the accelerations of the satellite in the ECI and ECEF coordinates, respectively, rECEF=xECEF2+yECEF2+zECEF2, (rECEF, φ, λ) is the altitude, latitude, and longitude of the satellite in the ECEF coordinate system; ECI/ECEF is the coordinate transformation matrix from ECI to ECEF coordinate system.The three components ofaECEF can be calculated by the following equation:(22)axECEF=xECEFrECEF∂U∂rECEF−xECEFzECEFrECEF2xECEF2+yECEF2∂U∂φ−yECEFxECEF2+yECEF2∂U∂λ,ayECEF=yECEFrECEF∂U∂rECEF−ybzrECEF2xECEF2+yECEF2∂U∂φ+xbxECEF2+yECEF2∂U∂λ,azECEF=zECEFrECEF∂U∂rECEF+xECEF2+yECEF2rECEF2∂U∂φ.The partial derivative of the potential functionU to (rECEF, φ, λ) is as follows:(23)∂U∂rECEF=−GMerECEF2∑n=0∞∑m=0nn+1RerECEFnP¯nmsinφC¯nmcosmλ+S¯nmsinmλ,∂U∂φ=GMerECEF2∑n=0∞∑m=0nRerECEFnPnm+1sinφ−mtanφPnmsinφC¯nmcosmλ+S¯nmsinmλ,∂U∂λ=GMerECEF2∑n=0∞∑m=0nRerECEFnmP¯nmsinφS¯nmcosmλ−C¯nmsinmλ,where G is the gravitational constant, Me is the mass of the Earth, Re is the average radius of the Earth’s equator, P¯nmx is the normalized Legendre polynomial, and C¯nm and S¯nm are the normalized gravitational potential coefficients, which can be read out directly from the Earth’s gravity field model file.


## 4.2. Atmospheric Drag Acceleration

The satellite acceleration caused by atmospheric drag can be written as follows:(24)adrag=−12CDAmρvr2ev,where m is the satellite mass, CD is the atmospheric drag coefficient, and its typical value ranges from 1.5 to 3.0. The relative velocity vr can be written approximately as follows:(25)vr=v−ω⊕×r,with the inertial satellite velocity vector v, the position vector r, and the Earth’s angular velocity vector ω⊕ of size 0.7292 × 10−4 rad/s. Note that the calculation of equations (24) and (25) is carried out in the true-of-date (TOD) coordinate system, and the corresponding value in the ECI coordinate system can be obtained after a coordinate transformation.The drag depends on the atmospheric densityρ at the location of the satellite, and the atmospheric density is usually calculated by the semiempirical atmospheric density model. At present, the commonly used atmospheric model in orbit determination is NRLMSISE-00 atmospheric model, which plays an important role in satellite orbit determination and prediction.The NRLMSISE-00 atmospheric model has 8 input items: The number of days from January 1 of the current year to the current day, the number of seconds from 00:00:00 to the calculating time, geographical longitude, latitude, altitude, the solar radiation flux of 10.7 cm on the previous day (F10.7), the average F10.7 for 81 days (3 solar rotation cycles, with the current day as the midpoint), average geomagnetic index (Ap) of current day and the twenty 3 h average values of Ap before the calculating time. The outputs include the number densities of N2, O2, He, Ar, N, H, O, and O+, the neutral atmospheric temperature, and the atmospheric density.These space environment parameters are provided by the SpaceWeather.txt file. Each line in the file represents a spatial environment parameter record for the specified date. The record has 32 fields, as shown in Figure1. The meanings of the fields related to atmospheric drag are shown in Table 1.Figure 1
Parameter format of spatial environment parameter record.Table 1
The meaning of space environment parameters.
ColumnNameDescription001–004yyyyYear006–007mmMonth009–010ddDay048–050Ap0Ap index (00:00–03:00)052–054Ap3Ap index (03:00–06:00)056–058Ap6Ap index (06:00–09:00)060–062Ap9Ap index (09:00–12:00)064–066Ap12Ap index (12:00–15:00)068–070Ap15Ap index (15:00–18:00)072–074Ap18Ap index (18:00–21:00)076–078Ap21Ap index (21:00–24:00)114–118ObsF10.7F10.7, previous day (observed)126–130Obs Lst81F10.7, arithmetic mean for the previous 81 days (observed)


## 5. Programming Program Design

This section introduces the program design idea of perturbation acceleration calculation and least-squares orbit improvement.Figure2 shows the program design flow chart of the calculation of the Earth’s gravitational acceleration. The input of the program is the position vector rECI, the coordinate transformation matrix from ECI to ECEF coordinate system, and the order of the Earth’s nonspherical gravity m, n. The calculation process is as follows: firstly, calculate the satellite’s coordinate rECEF, and then convert it into geodetic coordinate (rECEF, φ, λ), use the geodetic coordinate and the input nonspherical gravity order to calculate the value of the associated Legendre polynomial, then use the gravity potential coefficient file GGM03C.grv to calculate partial derivative according to equation (23), and finally use equation (21) to output aECI.Figure 2
The flow chart of the program design for the calculation of the Earth’s gravitational acceleration.Figure3 is a program flow chart for calculating the perturbation acceleration of atmospheric drag. The main calculation process in Figure 3 is as follows: the mean solar time, the number of seconds from 00:00:00 to the current time of the day, and the number of days from January 1 to the current day of the year are calculated using the input UTC time; the required Apgeo magnetic datum is obtained and calculated using the input space weather data document; the longitude and latitude of the satellite are calculated by the input satellite position vector. According to the above information, NRLMSISE-00 model is used to calculate the atmospheric density. Based on the position and velocity vectors of the satellite, the velocity of the satellite relative to the atmosphere is calculated according to the Earth’s angular velocity vector. Finally, the atmospheric drag acceleration can be calculated from the relative velocity, atmospheric density, input area-to-mass ratio, and drag coefficient.Figure 3
The flow chart of the program design for calculating the acceleration of atmospheric drag.Figure4 is a program design flow chart of the least-squares orbit improvement algorithm. The radar measurement data yio of each point in Figure 4 includes three components: range, azimuth, and elevation, namely (ρ,α,β); when using the numerical method to solve the reference trajectory, the perturbation factors should be considered as comprehensively as possible according to the accuracy requirement of the problem; while solving the error state transition matrix, only the center gravity of Earth and low-order nonspherical perturbation should be considered. The nonspherical perturbation of the Earth is introduced when solving the error state transition matrix, and the order is 4 × 4. When calculating the reference orbit, the order is 21 × 21.Figure 4
The flow chart of the program design for least-squares orbit improvement.The program is divided into three loops:The first loop is the internal loop for solving the numerical differential equation. The program’s inputs include Initial momentt0; Initial satellite orbit state x00=x0apr; N-point radar measurement data y10,y20,…,yN0; least-squares convergence threshold ε; number of iterations k = 0; measurement data index number i = 1; other relevant parameters, such as EOP files, atmospheric environment parameters, satellite surface-to-mass ratio, etc. The output results are the reference trajectory and error state transfer matrix at time i. The second loop is to traverse each measurement value and solve the single orbit improvement according to the least-squares principle. After iteration, update the orbital state at t0: x0k+1=x0k+HTWH−1HTWΔz. The third loop is to perform multiple least-squares orbit improvements to make the final result close to the actual orbit state. The output is the improvement orbit elements x0lsq of t0 time.Getting the time interval ∆T and current orbit state of the satellite xtc=r,v from the input parameters of the differential equation in the inner loop, we can follow these steps to solve the numerical differential equation:(1)
Current timetc = t0+∆t < ti(2)
According to the force model and related parameters, calculate accelerations of two-body, nonspherical perturbation, planets, atmospheric drag, etc., so as to construct the vectorx˙tc=v;a(3)
Calculate∂atwo−body+anon−sphericalperturbation/∂r to construct matrix F, then solve the differential equation Φ˙tc=FΦtc


## 6. Software Interface Call and Validity Verification

### 6.1. Earth’s Nonspherical Gravitational Acceleration

The rodAccelHarmonic interface function for calculating the Earth’s gravitational acceleration is realized by C language according to the programming principle of Part 5, and it is integrated into the RadarOrbDet library. The input and output of this function are shown in Figure2.A satellite is simulated in STK11.2, and its orbital elements are set as shown in Table2. In order to compare STK with rodAccelHarmonic function under the same perturbation conditions, only the nonspherical gravitational perturbation of the Earth (21st order) is set in STK. The simulation time period is set from 2020-10-18 04:00:00 UTCG to 2020-10-19 04:00:00 UTCG. UTCG stands for Gregorian Universal Time Coordinated. Use STK’s report function to output the satellite epoch time every 1 minute and the position and acceleration under ECI coordinate system. The epoch time is used to calculate the coordinate transformation matrix. The calculated matrix and the position in ECI coordinate system are used as the input parameters of the rodAccelHarmonic function. The acceleration output from the report is taken as the standard value and compared with the acceleration calculated by the rodAccelHarmonic function. The relative error described in equation (26) is used for quantitative analysis.(26)relative error=asoftware−aSTKaSTK.Table 2
Satellite orbit elements.
EpochSemimajor axis (km)EccentricityInclinationArgument of perigeeAscension of ascending nodeTrue anomaly2020-10-08 04:00:006678.14028.5°0°0°0°There are 1441 data points in total. Figure5 plots the acceleration (solid line) and relative error curve (dashed line) calculated by the rodAccelHarmonic function. It can be seen from the figure that the relative error is at the level of 10−6, which shows that the software algorithm is effective. The main source of error is the loss of accuracy due to the truncation of the text output of the STK report. In addition, the acceleration curve in the figure shows periodic changes because when the only nonspherical gravitational perturbation is considered, the semimajor axis of the satellite will change periodically.Figure 5
Earth’s gravitational acceleration and its error.


### 6.2. Atmospheric Drag Acceleration

According to the program design principle in Part 5, the rodAccelAirDragPerturbation interface function for calculating the atmospheric drag perturbation acceleration is realized by using C language and integrated into RadaOrbDet library. The input and output of this function are shown in Figure3.Two satellites are simulated in STK11.2, and their initial orbital elements are set to be the same, as shown in Table2, except that the epoch time is changed to 2016-06-16 04:00:00. The first satellite is set with a 21st order gravity model pulse atmospheric drag perturbation, and the second satellite is only set with a 21st order gravity model. Using the report function of STK, the accelerations of the two satellites at 2016-06-16 04:00:00 (STK11 does not contain the latest atmospheric model data, so use this time to ensure that our software and STK use the same atmospheric model data) are obtained respectively, and the atmospheric drag acceleration of the second satellite can be obtained by subtracting the two accelerations. In addition, another atmospheric drag acceleration can be calculated by calling the rodAccelAirDragPerturbation interface function. The calculation results of STK and our software are shown in Table 3, and the relative error described in equation (26) is used for quantitative comparison.Table 3
Comparison of atmospheric drag acceleration (STK and out software).
CategoryItemAcceleration (m/s2)xyzSTKSatellite 1−8.9509556686−1.123596 × 10−4−1.37593 × 10−5Satellite 2−8.9509556684−1.037988 × 10−4−8.7521 × 10−6Atmospheric drag perturbation−2 × 10−10−8.5608 × 10−6−5.0072 × 10−6SoftwareAtmospheric drag perturbation−6.6174 × 10−24−8.4633 × 10−6−4.9502 × 10−6Relative error1.14%It can be seen from Table3 that the error is 1.14%, indicating that the software algorithm in this paper is effective. The main sources of error are the coordinate and time conversion of the calculation program, and the calculation truncation error, etc.


### 6.3. Least Squares Orbit Improvement

Using the programming principle in part 5, the least-squares orbit improvement interface function rodLeadsquare is realized and integrated into the RadarOrbDet library.  The input and output  of this function are shown in Figure4.A satellite is simulated in STK11.2 and its propagator mod is set to two-body. Its orbital elements settings are shown in Table2. Only the epoch time is changed to 2021-09-07 04:00:00. The longitude, latitude, and altitude of the radar site are (120°, 30°, 0). Select 10 data points with an interval of 1 second from 19:15:01 to 19:15:10, and add measurement errors of 100 m, 0.1°, and 0.1° to the distance, azimuth, and elevation, respectively.The software in this paper and the ODTK software are used to compare the orbit determination, respectively. The calculation results of ODTK and the software are shown in Table4. The relative error of position and velocity described in equation (27) are used for quantitative comparison.(27)relative position error=rsoftware−rODTKrODTK,relative velocity error=vsoftware−vODTKvODTK.Table 4
Comparison of orbit determination results (ODTK and our software).
CategoryOrbit determination results (epoch: 19:15:00; unit: km, km/s)xyzvxvyvzODTK5187.983690.772002.66−4.845.292.88Software5190.663692.612004.97−4.865.282.87Relative error0.06%0.32%It can be seen from Table4 that the position error of orbit determination is 0.06%, and the velocity error is 0.32%, indicating that the software algorithm in this paper is effective.


## 6.1. Earth’s Nonspherical Gravitational Acceleration

The rodAccelHarmonic interface function for calculating the Earth’s gravitational acceleration is realized by C language according to the programming principle of Part 5, and it is integrated into the RadarOrbDet library. The input and output of this function are shown in Figure2.A satellite is simulated in STK11.2, and its orbital elements are set as shown in Table2. In order to compare STK with rodAccelHarmonic function under the same perturbation conditions, only the nonspherical gravitational perturbation of the Earth (21st order) is set in STK. The simulation time period is set from 2020-10-18 04:00:00 UTCG to 2020-10-19 04:00:00 UTCG. UTCG stands for Gregorian Universal Time Coordinated. Use STK’s report function to output the satellite epoch time every 1 minute and the position and acceleration under ECI coordinate system. The epoch time is used to calculate the coordinate transformation matrix. The calculated matrix and the position in ECI coordinate system are used as the input parameters of the rodAccelHarmonic function. The acceleration output from the report is taken as the standard value and compared with the acceleration calculated by the rodAccelHarmonic function. The relative error described in equation (26) is used for quantitative analysis.(26)relative error=asoftware−aSTKaSTK.Table 2
Satellite orbit elements.
EpochSemimajor axis (km)EccentricityInclinationArgument of perigeeAscension of ascending nodeTrue anomaly2020-10-08 04:00:006678.14028.5°0°0°0°There are 1441 data points in total. Figure5 plots the acceleration (solid line) and relative error curve (dashed line) calculated by the rodAccelHarmonic function. It can be seen from the figure that the relative error is at the level of 10−6, which shows that the software algorithm is effective. The main source of error is the loss of accuracy due to the truncation of the text output of the STK report. In addition, the acceleration curve in the figure shows periodic changes because when the only nonspherical gravitational perturbation is considered, the semimajor axis of the satellite will change periodically.Figure 5
Earth’s gravitational acceleration and its error.


## 6.2. Atmospheric Drag Acceleration

According to the program design principle in Part 5, the rodAccelAirDragPerturbation interface function for calculating the atmospheric drag perturbation acceleration is realized by using C language and integrated into RadaOrbDet library. The input and output of this function are shown in Figure3.Two satellites are simulated in STK11.2, and their initial orbital elements are set to be the same, as shown in Table2, except that the epoch time is changed to 2016-06-16 04:00:00. The first satellite is set with a 21st order gravity model pulse atmospheric drag perturbation, and the second satellite is only set with a 21st order gravity model. Using the report function of STK, the accelerations of the two satellites at 2016-06-16 04:00:00 (STK11 does not contain the latest atmospheric model data, so use this time to ensure that our software and STK use the same atmospheric model data) are obtained respectively, and the atmospheric drag acceleration of the second satellite can be obtained by subtracting the two accelerations. In addition, another atmospheric drag acceleration can be calculated by calling the rodAccelAirDragPerturbation interface function. The calculation results of STK and our software are shown in Table 3, and the relative error described in equation (26) is used for quantitative comparison.Table 3
Comparison of atmospheric drag acceleration (STK and out software).
CategoryItemAcceleration (m/s2)xyzSTKSatellite 1−8.9509556686−1.123596 × 10−4−1.37593 × 10−5Satellite 2−8.9509556684−1.037988 × 10−4−8.7521 × 10−6Atmospheric drag perturbation−2 × 10−10−8.5608 × 10−6−5.0072 × 10−6SoftwareAtmospheric drag perturbation−6.6174 × 10−24−8.4633 × 10−6−4.9502 × 10−6Relative error1.14%It can be seen from Table3 that the error is 1.14%, indicating that the software algorithm in this paper is effective. The main sources of error are the coordinate and time conversion of the calculation program, and the calculation truncation error, etc.


## 6.3. Least Squares Orbit Improvement

Using the programming principle in part 5, the least-squares orbit improvement interface function rodLeadsquare is realized and integrated into the RadarOrbDet library.  The input and output  of this function are shown in Figure4.A satellite is simulated in STK11.2 and its propagator mod is set to two-body. Its orbital elements settings are shown in Table2. Only the epoch time is changed to 2021-09-07 04:00:00. The longitude, latitude, and altitude of the radar site are (120°, 30°, 0). Select 10 data points with an interval of 1 second from 19:15:01 to 19:15:10, and add measurement errors of 100 m, 0.1°, and 0.1° to the distance, azimuth, and elevation, respectively.The software in this paper and the ODTK software are used to compare the orbit determination, respectively. The calculation results of ODTK and the software are shown in Table4. The relative error of position and velocity described in equation (27) are used for quantitative comparison.(27)relative position error=rsoftware−rODTKrODTK,relative velocity error=vsoftware−vODTKvODTK.Table 4
Comparison of orbit determination results (ODTK and our software).
CategoryOrbit determination results (epoch: 19:15:00; unit: km, km/s)xyzvxvyvzODTK5187.983690.772002.66−4.845.292.88Software5190.663692.612004.97−4.865.282.87Relative error0.06%0.32%It can be seen from Table4 that the position error of orbit determination is 0.06%, and the velocity error is 0.32%, indicating that the software algorithm in this paper is effective.


## 7. Conclusion

To the issue of batch orbit improvement in space target radar orbit determination, this paper introduces the principle of least-squares orbit improvement, describes the calculation method of partial derivatives in the model, and explains the Earth’s nonspherical gravity and atmospheric drag perturbations. Program design ideas and simulations show the effectiveness of software design. In addition, the RadarOrbDet library developed in this paper can well assist the designers of space target surveillance radar systems in the rapid design and verification of orbit determination indicators.



---

*Source: 2341442-2022-02-24.xml*