Determination of the Orbit of an Unknown Ultra-Small Spacecraft Based on the Circular Perturbed Motion Model and Measurements of the Doppler Frequency Shift

In this work, the determination of the orbit of an unknown ultra-small spacecraft based on Doppler measurements of the telemetry signal frequency is studied. The reception and processing of telemetry radio signals was carried out by the ground station of the Belarusian State University. In the model of perturbed circular motion, the radio signals of a small satellite were processed and the parameters of its orbit were determined. Based on a probabilistic estimate of the elevation angle and Doppler frequency shift of an ultra-small spacecraft from 10–20 measurements, a set of orbital parameters is determined for the estimated time of receiving telemetry signals. For antenna systems, the dynamics of changes in the elevation angle, azimuth, and Doppler frequency shift of telemetry radio signals for the next flights of an unknown ultra-small spacecraft was predicted. The calculated absolute errors in predicting the elevation angle, azimuth, and Doppler frequency shift did not exceed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3^{\circ}$$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3^{\circ}$$\end{document} and 250 Hz, respectively, which is sufficient for successful telemetry reception and decoding. Using the NORAD database of orbital parameters, the Chinese nanosatellite LUOJIA-1 01 of the Cubesat (6U) format (number 43485 in the NORAD system) was identified.


INTRODUCTION
The role of ultra-small spacecraft (weighing up to 10 kg) in the study, development, and use of nearearth space is constantly increasing [1][2][3]. The advantage of ultra-small spacecraft (USSC) is their speed and low cost of development due to the possibility of using commercial components, fewer ground tests, ease of execution and the low cost of the ground control and communication segment (most often one ground receiving and control station), as well as the ability to create groups of such satellites for solving commercial and scientific problems [4][5][6][7]. Therefore, many universities develop their own USSCs as educational laboratories for training specialists in aerospace specialties, conducting scientific experiments, and demonstrating technological capabilities [6][7][8].
Prediction of future positions of USSCs is important for solving problems of their control, tracking, and conducting radio communication sessions. Traditionally, the SGP4 prediction model [9,10] is used to make these predictions for one ground station. The SGP4 model uses averaged orbital parameters in the TLE (two-line elements) format of the NORAD (North American Aerospace Defense Command) [10][11][12]. As a rule, TLEs are updated daily and are available free of charge, but in the long term or in the event of military conflicts, the NORAD system has the ability to disable public access to the database of averaged orbital parameters. Less often, on board a small spacecraft, a navigation receiver is used to determine the exact coordinates and speed [8].
An alternative way to obtain initial data for predictive models of the motion of a USSC is to determine the orbit based on measurements of the characteristics of radio signals of a telemetric or command radio link, which simplifies the hardware implementation and the cost of the measuring system [13,14]. The measured parameters for a budgetary ground receiving station (GRS) are the time and the Doppler frequency shift of the received radio signal [15]. Most often, the processing of measurements consists in the refinement of the orbital parameters using the known initial data. The database of averaged orbital parameters in the TLE format contains information on 44 336 space objects (as of June 23, 2019) with dimensions over 10 cm in diameter. Several dozen unlicensed USSCs (unregistered in international organizations, that is the International Telecommunication Union (ITU)) have been launched, spacecraft with a transverse dimension of less than 10 cm, which include picosatellites and fempto satellites (weighing up to 0.1 kg), are especially dangerous Thus, for example, in January 2018, four unregistered SpaceBEE picosatellites of the American company Swarm Technologies Inc. were launched on the Indian PSLV-XL rocket. The only way to determine the orbits and track such satellites is by radio-technical trajectory measurements via telemetry or a command radio link.
The task of predicting the orbital parameters in the first hours after the launch of a USSC is especially acute, when there is no information in the NORAD system database. Therefore, the task of the initial determination of the orbital parameters of a USSC from the results of radio-technical trajectory measurements for the purpose of its further identification and cataloging is urgent.
Analysis of the orbital parameters of tracked space objects (SOs) in low orbits (an altitude up to 2000 km) using the NORAD system database (as of January 1, 2020) showed that the vast majority of spacecraft (more than 90%) move in orbits close to circular. Nano-and picosatellites developed according to the Cubesat standard can serve as an example of USSCs. These low-orbit satellites (their orbital altitude is less than 1000 km) have an almost circular orbit (with a small eccentricity value e 1). Table 1 shows the values of the eccentricities of the orbits of space objects (SOs) in low orbits. Thus, for example, the total number of SOs in low orbits launched since 1990 is 2296, while the number of SOs with an orbital eccentricity less than 0.005 is 2100 (∼91.5%).
The analytical theory of motion of artificial satellites with small eccentricities was described in [16][17][18][19]. In [20], the possibility of using the model of perturbed circular motion for the satellites of Jupiter was studied.
In this paper, we developed a method for determining the orbit of an unknown ultra-small spacecraft based on a model of circular perturbed motion and experimental measurements of the Doppler frequency.

THEORETICAL MODEL
The mathematical model of the motion of a USSC in a circular orbit is determined by the following parameters of the state vector X: where i is the inclination of the orbit to the equatorial plane, T is the orbital period, u is the latitude argument, and Ω is the longitude of the ascending node at time t. The determination of the orbit of an unknown USSC was carried out by the University ground station based on N measurements of telemetry radio signals at several orbits for the period from Oct. 6, 2019 to Oct. 9, 2019 (UTC) where t k is the time of receiving and Δf k is the Doppler frequency shift of the received telemetry radio signals.
For the calculated point in time t 0 = 09 : 48 : 18 Oct. 10, 2019 (UTC), the USSC state vector was found X 0 = (T 0 , i 0 , u 0 , Ω 0 , ), which best meets the measurement results according to two criteria: • the elevation angle of the USSC above the ground receiving station (GRS) at the times of measurements t k must be positive (the USSC relative to the GRS should be above the horizon): • the Doppler frequency shift Δf calc k = Δf calc (t k ) of the telemetry radio signal obtained as a result of numerical simulation at the time moments of measurements t i should differ from the measured one Δf The direction to the South X The geometry of the perturbed circular motion model for determining the orbit.
less than 300 Hz (here we consider the error associated with the frequency instability of the on-board USSC transmitter): According to the results of measurement processing, the ranges of variation ΔX = (ΔT, Δi, Δu, ΔΩ) of the parameters of the state vector X 0 are determined. Next, for each obtained state vector X 0 at the measurement time t k the perturbed circular motion model is used to calculate the elevation angle el k of the USSC above the GRS and the Doppler frequency shift Δf calc k of the telemetry radio signal using the following algorithm.
1. At the time of measurements t k , we find the orbital parameters T (t k ) = T 0 , i(t k ) = i 0 , Ω(t k ) and u(t k ). When calculating the longitude of the ascending node Ω(t k ) and the argument of latitude u(t k ) in the model of disturbed circular motion, we considered the noncentrality of the Earth's gravitational field through the fundamental second zonal harmonics in the expansion for the potential of the gravitational force [20]: Here, X 0 = (T 0 , i 0 , u 0 , Ω 0 ) are the parameters of the state vector at the time t 0 , R = μT 2 (t k )/ 4π 2 1/3 is the radius of the satellite orbit, R E = 6378.137 km is the mean equatorial radius of the Earth, n = 2π/T (t k ) is the angular velocity of the satellite (mean motion), and J 2 = 0.0010826267 is the second zonal harmonic. The geometry of the perturbed circular motion model for determining the orbit is shown in Fig. 1.
2. Find the coordinates and projections of the velocity vector in the orbital coordinate system (CS) To transform coordinates and velocities from the CS OX orb Y orb Z orb to the geocentric inertial CS OXY Z, it is necessary to make two successive rotations by an angle opposite to the inclination of the i orbit around the OX orb axis and the angle opposite to the longitude of the ascending node Ω around the axis OZ orb : The rotation matrices from the CS OX orb Y orb Z orb to the rectangular coordinate system OXY Z have the form: 4. According to the coordinates of the ground receiving station (ϕ = 53 • 54 27 north latitude, λ = 27 • 33 52 east longitude, altitude H = 230 m) at the time of measurement t k we determine the radius vector of the position R GRS and the velocity vector V GRS in the geocentric inertial coordinate system: where ω E is the vector of the angular velocity of the Earth's rotation, θ LST is the local sidereal time for the GRS at the time t k , r b , r k are the projections of R GRS onto the equatorial plane and the axis perpendicular to this plane. 5. Determine the radius vector of the slant range ρ and the rate of its change in the OXYZ coordinate system: 6. Find the vector of the slant range ρ and the projection of its velocity in the topocentric rectangular CS O 1 SEU : where the rotation matrix R 2 (α) is defined as 7. Find the slant range ρ, the elevation angle el and the azimuth az, as well as the rate of their change at the time t k : Here, S, E, U are the projections of the inclined range vector ρ on the coordinate axis of the O 1 SEU system. 8. Calculate the Doppler frequency shift where c is the speed of light in vacuum, f 0 is the nominal frequency of the onboard USSC transceiver. The probability of success β j of the given set of orbital parameters (T 0 , i 0 , u 0 , Ω 0 ) for the time t 0 = 09 : 48 : 18 09. 10.2019 (UTC) by estimating only the calculated elevation angle el (j = 1) and the elevation angle and the Doppler frequency shift Δf calc (j = 2) at the moment of receiving the telemetry radio signal was found according to where N j is the number of calculated points that satisfy the criterion (1) for j = 1 and (1), (2) for j = 2 and N total = 10, 20 is the total number of measurement points at which the numerical simulation was carried out for a given set of orbital parameters (T 0 , i 0 , u 0 , Ω 0 ). The probability of success of β 2 = 50% means that the N 2 = N total /2 conditions (1) el > 0 and (2) |Δf − Δf calc | < 300 Hz are satisfied only in half of the experimental data.

RECEPTION AND PROCESSING OF RADIO SIGNALS OF A USSC
The reception and processing of radio signals from an unknown USSC was carried out by the ground receiving station (GRS) of the Belarusian State University on several orbits for the period from Oct. 6, 2019 to Oct. 9, 2019 (UTC). The GRS hardware part consists of: wave channel antennas with circular polarization in the range of 435-438 MHz; a receiving systems based on the IC-9100 transceiver; a receiving system based on the SDR receiver module;  a rotary device YAESU G-5500 with the control unit; a control computer. The software part of the GRS includes: the software (SW) for predicting the motion of a USSC and the characteristics of radio signals; the software for modeling and visualization of scenarios of GRS operation and carrying out express calculation of standard navigation and ballistic information; and the software for receiving and processing telemetry radio signals. The system for measuring and determining the orbit for the university ground station consists of a GPS receiver, a module for measuring the frequency and time of receiving radio signals based on a microcontroller (MC) for time processing and a two-channel digital oscillograph, and the software for processing measurements, the software for determining and refining orbits. The system for measuring and determining the orbit of the receiving ground station can measure the satellite's orbit both with initially known parameters (the tracking mode) and with unknown parameters (the omni-directional search mode). To process measurements in the tracking mode, orbital elements in the TLE format or from the own database are used as initial data. The software for predicting the motion of the satellite by the known initial parameters of the orbits calculates the radio communication sessions, the estimated parameters of the tracking systems of the rotary devices, and the frequencies of the received radio signals. The main task is to check the parameters of the satellite orbit and to refine them. The calculations use the method of differential correction of orbital parameters based on the rate of change in the slant range between the ground station and the satellite, calculated from the Doppler frequency shift of the received radio signal.
When measuring in the omni-directional search mode, the problem of the initial determination of the satellite orbit parameters is solved. The frequency of the target satellite is searched within the 435-445 MHz amateur radio band. The distinctive features of telemetry packet radio signals (packet repetition rate, presence of tags, etc.) are determined. Based on the results of measurements on several orbits, the frequency of the radio signal, the period of the satellite, the maximum duration of flight over the ground station are estimated, and the allowable limit of the measurement error is estimated. To process measurements on one satellite pass over the ground station, algorithms for determining the state vector based on the unperturbed motion model are used, and on several passes algorithms for determining the orbit are based on the simplest models of disturbed motion.

DISCUSSION
The total number of Cubesat satellites in the NORAD database is ∼180. To narrow the search ranges for the parameters of an unknown satellite, we preliminarily analyzed the possible orbital parameters of satellites based on the averaged orbital parameters in the TLE format as of October 1, 2019. Calculations show that in Fig. 2 there are two types of characteristic Cubesats orbits. The first type of orbits is associated with nano-and picosatellites launched from the International Space Station (ISS). They have an orbital inclination of about 51.6 • , an orbital period from 90 to 94 min, and an orbital altitude from 330 km to 430 km. The second type of orbit is associated with nano-and picosatellites launched into a sun-synchronous orbit by a passing launch jointly with a main satellite (usually a loworbit remote sensing satellite). They have an orbital inclination from 98 • to 99 • , an orbital period from 94 to 102 min, and an orbital altitude from 500 km to 800 km.
The ground receiving station of the Belarusian State University received and processed telemetry radio signals from an unknown spacecraft in the range of radio amateur frequencies 435-445 MHz. Since the orbital parameters of the spacecraft were unknown, it was impossible to predict the exact pointing angles of the antenna systems and the Doppler frequency shift of the received telemetry radio signal. Due to the low power of the radio signal, telemetry packets could not be decoded.
In the calculations, the intervals between radio signals Δt 2 on adjacent loops were 90, 96 (mostly), and 102 min. This made it possible to assume that this unknown small satellite belongs to the satellites with orbits of the second type. For the determination algorithm, the ranges of variation of the parameters of the circular orbit were set: the orbital period T : 94-102 min with a step of 1 s, the inclination of the orbit Table 2. The ranges of variation of orbital parameters with a probability of success β 1 higher than 50% based on an estimate of only the elevation angle of an unknown ultra-small spacecraft based on N = 10, 15, 20 measurements    Table 1 shows the ranges of variation of the orbital parameters for the estimated time t 2 with the probability of success β 1 above 50% based on the analysis of the elevation angle of an unknown ultrasmall spacecraft by N total = 10, 15, 20 measurements. It was found that the ranges of variation of the orbital period T and orbital inclination i with an increase in the number of measurements from N total = 10 to 20 are reduced by 24% from 120 to 91 s and from 0.38 • to 0.29 • , respectively, while the ranges for the latitude u and longitude of the orbit's ascending node Ω remain unchanged. Table 2 shows the ranges of variation of the orbital parameters for the estimated time t 2 with the probability of success β 2 above 50% based on the estimate of the elevation angle and the Doppler frequency shift of the unknown ultra-small spacecraft by N total = 10, 15 and 20 measurements. It was found that the ranges of variation of the orbital period Т and orbital inclination i with an increase in the number of measurements from N total = 10 to 20 decrease from 44 to 7 s (by 84%) and from 0.14 • to 0.02 • (by 86%), respectively. With an increase in the number of measurements from N total = 10 to 20 the ranges of variation of the argument of latitude u and longitude of the ascending node of the orbit Ω decrease from 21 • to 4 • (by 81%) and from 10 • to 6 • (by 40%), respectively. Figure 3 shows the dependence of the number of possible sets of orbital parameters X on the success probabilities of the given set of orbital parameters β 1 and β 2 based on the measurement data for 10 and 20 points. During analysis at the calculated points only by the values of the elevation angle el > 0 (Fig. 3a), there are 16 106 (N total = 10) and 20 805 (N total = 20) sets of orbital parameters N X with a probability of success in the range of 50-60%. With a 100% probability, the number of possible orbital parameters is 5065 (N total = 10) and 1099 (N total = 20), which does not allow its unique identification.
Analysis by two parameters (elevation el and Doppler frequency shift Δf calc ) allows reducing the number of possible state vectors X by 2 orders of magnitude. It was found that there are 187 (N total = 10) and 42 (N total = 20) sets of orbital parameters N X with the probability of success in the interval of 50-60% and one set of orbital parameters X = (5855с, 97.98 • , 115 • , 360 • ) and (5855 с, 97.98 • , 115 • , 359 • ) with the probability of success of β 2 = 95% and 100% respectively. Thus, a statistical analysis of the orbital parameters of an unknown ultra-small spacecraft based on an estimate of the elevation angle and the Doppler frequency shift makes it possible to uniquely determine the orbit.
For the given set of orbital parameters for the calculated time t 2 = 9 : 48 : 18 09. 10   angle el, azimuth az and Doppler frequency shift of telemetry radio signals Δf calc were numerically simulated in the range from 0 : 0 : 0 to 23 : 59 : 59 10.10.2019 (UTC) of next flights. Based on the simulated data, the telemetry packets of the small satellite LUOJIA-1 01 (№ 43485 in the NORAD system) were successfully received and decoded at the flight interval from 10 : 07 : 50 to 10 : 20 : 50 10.10.2019 (UTC), which is confirmed by good agreement between the calculated Doppler curves and the measurements of the BSU GRS, as shown in Fig. 4. The results were published on the DK3WN SatBlog amateur radio site. Figure 5 shows the results of estimating the accuracy of predicting the elevation angle el, azimuth az and the Doppler frequency shift of radio telemetry signals Δf D in the interval of successful reception and decoding of telemetry packets of the small LUOJIA-1 01 satellite in the perturbed circular motion model in comparison with the SGP 4 model. As can be seen from the graphs, the absolute errors in predicting the elevation and azimuth did not exceed 3 • and the absolute error in predicting the Doppler frequency shift of telemetry radio signals did not exceed 250 Hz, which is sufficient for the successful reception of telemetry radio signals and their decoding.

CONCLUSIONS
In this work, a method for calculating satellite orbit was developed on the basis of a model of circular disturbed motion and measurements of the Doppler frequency shift. A method for determining an orbit has been developed, considering a preliminary analysis of the possible orbital parameters of satellites based on averaged orbital parameters in the TLE format. The proposed calculation algorithm made it possible to identify the Chinese nanosatellite of the Cubesat (6U) format LUOJIA-1 01 (No. 43485 in the NORAD system), which is confirmed by the successful reception of telemetric radio signals and their decoding, as well as good agreement between the calculated and experimental Doppler curves.