GPS-based satellite formation flight simulation and applications to ionospheric remote sensing ============================================================================================== * YuXiang Peng * Wayne A. Scales * Thom R. Edwards ## Abstract The Virginia Tech Formation Flying Testbed (VTFFTB), a GNSS-based hardware-in-the-loop (HIL) simulation testbed for spacecraft formation flight, is developed and applied to ionospheric remote sensing. The current VTFFTB consists of GNSS RF hardware signal simulators, multi-constellation multifrequency GNSS receivers, a navigation and control system, an STK visualization system, and an ionospheric remote sensing system. GPS signals are emulated using GNSS simulator scenarios that include ionospheric phenomena. A formation of two spacecraft (“chief” and “deputy”) is considered. GNSS receiver data are used to produce space-based Total Electron Content (TEC) and scintillation measurements. A reference low Earth orbit (LEO) scenario is benchmarked with past simulation results to validate functionality. A LEO formation flying mission is designed to probe two Equatorial Spread F (ESF) scenarios with plasma bubbles. The results investigate the structure of ionospheric irregularities and demonstrate that the GPS-based satellite formation is able to measure vertical electron density by differencing 1D GPS vertical TEC. ## 1 INTRODUCTION Spacecraft formation flying, a novel space mission architecture of multiple satellites orbiting together in a specific configuration, is significantly advantageous to space science and engineering. The concept can be dated back to the early 1960s, when spaceflight pioneers investigated the spacecraft rendezvous and docking problem. A critical technique of spacecraft formation flying is being able to precisely control the relative motion of collaborative distributed space systems. Compared to traditional single-spacecraft missions, formation flight missions can be more cost-efficient, agile, flexible, scalable, and sustainable by deploying a group of collaborative small satellites. A number of CubeSat (a class of small satellites) formation flying missions for various applications including Earth science, astronomy and astrophysics, planetary science, heliophysics, and technology demonstrations are presented by Bandyopadhyay et al.1 In principle, the observation dimension and resolution can be increased by deploying scientific instruments on a cluster of satellites with the ability of rearranging their formation geometry. The European Space Agency (ESA)’s Cluster II mission is one of the pioneering spacecraft formation flying missions for space science. In 2000, a constellation of four identical Cluster II spacecraft was launched into highly elliptical polar orbits around Earth and formed a tetrahedral formation to observe the Sun-Earth electromagnetic interactions by making three-dimensional in situ measurements.2 After the success of Cluster II, NASA took a step forward in 2015 and launched the Magnetospheric Multiscale (MMS) mission, which consists of four identical satellites to probe Earth’s magnetosphere in higher spatial and temporal resolution with a concentration on magnetic reconnection.3 The Gravity Recovery and Climate Experiment (GRACE), led by NASA and the German Aerospace Center (DLR), is another renowned formation flight mission. In March 2002, two GRACE twin satellites were launched into low Earth orbit (LEO) with an in-track separation of 220 km to map Earth’s gravity field with a high level of accuracy.4 Eleven years later, ESA launched a constellation of three SWARM satellites into polar orbits to precisely survey Earth’s magnetic field, ionosphere, and other geoscience problems.5 When applying formation flying concepts to those satellite imaging or space-borne interferometry systems, the aperture becomes freely adjustable by altering the relative distances among a team of satellites. This opens up new design opportunities for more powerful remote sensing systems suitable for next-generation space exploration missions, such as the Stellar Imager (SI), the Milli-Arc-Second Structure Imager (MASSIM), and Black Hole Imager (BHI) missions. All of these future NASA formation flying astrophysical missions are planned to be launched in the 2030s.6 Spacecraft positioning, navigation, and timing (PNT) play an important role in satellite missions. The Global Navigation Satellite System (GNSS) is the most ubiquitous modern real-time PNT technology and is commonly utilized in navigation for spacecraft formation flight. In LEO, GNSS can provide absolute navigation for each individual satellite and highly accurate relative navigation based on differential GNSS measurement techniques. Besides the well-known American Global Positioning System (GPS) and the Russian GLONASS, the European Union’ Galileo and the Chinese BeiDou (or COMPASS) are two other rapidly growing global GNSS constellations. GNSS also provides space weather assessment capabilities. The ionosphere impacts the propagation of GNSS radio-frequency (RF) signals between the GNSS constellations and receivers. As a dominant source of GNSS positioning errors, ionospheric plasma can cause time delays, scintillation, Faraday rotation, or even loss of lock of GNSS signals and thus significantly impacts receiver performance. Consequently, GNSS is widely applied to ionospheric remote sensing. For example, multifrequency GNSS receivers have been used broadly for ground-based observations of Total Electron Content (TEC) and ionospheric scintillation. In addition, space-based satellite missions with onboard GPS receivers are designed for ionospheric investigations. For example, the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) mission satellites utilize radio occultation (RO) techniques for three-dimensional ionospheric tomography.7 TEC and scintillation measurements are able to reveal or verify many ionospheric phenomena such as electron density irregularities, spread F bubbles, plasma plumes, mid-latitude troughs, storm-enhanced densities (SED), and sub-auroral polarization streams (SAPS). The observational capabilities of traditional ground-based (eg, radar) or space-based single-satellite techniques (eg, in situ, tomography) cannot always meet the specific temporal and spatial demands to detect these ionospheric phenomena, where GNSS-based multi-satellite remote sensing applications present the possibility of fulfilling these demands. Plasma bubbles, polar cap patches, tongues of ionization, etc may have spatial scales from hundreds of km down to meter-scale irregularities driven by secondary instabilities.8 Ionospheric irregularities with a few hundreds of meter spatial scales corresponding to the Fresnel length of GNSS signals can be a potential source causing amplitude scintillation.9 It is difficult for a single satellite to probe the boundary of complicated plasma plumes that vary with time and space. GPS-RO with a satellite constellation (eg, COSMIC) cannot pinpoint the exact locations of ionospheric irregularities that cause scintillation as well. To combine the observational advantages of satellite formation flying and the growing GNSS network, GNSS-based spacecraft formation flying in LEO can be applied to ionospheric remote sensing, irregularity tracking, and morphology. The robustness of spacecraft formation flying missions highly depends on the performance of the Guidance, Navigation, and Control (GNC) system. A platform to test the functionality and reliability of all the system components is needed to demonstrate mission feasibility and implementation. To effectively develop and validate GNC algorithms, assess performance of onboard GNSS receiver(s), and emulate various high-fidelity mission scenarios, hardware-in-the-loop (HIL) simulations are essential to the mission designs of GNSS-based spacecraft formation flying. Since the early 2000s, NASA’s Goddard Space Flight Center (GSFC) had used the Formation Flying TestBed (FFTB) to support their satellite formation flying mission (eg, MMS) designs, validations, and developments. Leitner first established the FFTB, a GPS-based HIL simulation platform for formation flying clusters and constellations of satellites.10 Busse used the FFTB to demonstrate centimeter-level precision LEO formation flying using single-differential carrier-phase GPS based on an adaptive extended Kalman filter (EKF).11 Burns upgraded the FFTB and then performed both pure software and HIL formation control simulations, which then showed that the mean range error of HIL simulation is about 2.5 times worse than the error of pure software simulation.12 Mohiuddin used the FFTB to test and evaluate a satellite relative navigation estimator based on a double-differential carrier-phase GPS technique.13 The German Aerospace Center (DLR) has also actively involved HIL simulation of spacecraft formation flying. Gill used the FFTB to demonstrate real-time HIL simulation for GPS-based autonomous formation flying of two near-polar LEO micro-satellites.14 Leung conducted closed-loop GPS-based HIL spacecraft formation flying simulations to validate the real-time performance of a GPS relative navigation algorithm in a four-spacecraft scenario.15 Yamamoto presented a closed-loop testbed for GPS-based HIL autonomous formation flying of LEO satellites developed at DLR and used it to test and validate the GPS-based GNC flight software for PRISMA, a Swedish technology demonstration mission of spacecraft formation flying and rendezvous.16 Besides NASA and DLR, several other universities have also developed GNSS-based HIL testbeds to support relative navigation and control study or mission designs for satellite formation flying, such as Marji at University of Calgary,17 Eyer at University of Toronto,18 and the D’Amico research group at Stanford University.19 Park et al20 at Yonsei University established a GPS-based HIL satellite formation flying simulation testbed and developed algorithms for relative navigation and formation control. These testbeds were primarily developed to verify real-time HIL algorithms and/or test space-based GNSS navigation systems however, none of them were prototyped for ionospheric remote sensing applications. Because of the similarity of infrastructure, Park et al20 is chosen as the paradigm for testbed development and benchmarking for validations in this paper. The main purpose of this work is to develop a GNSS-based HIL simulation testbed for small-satellite formation flying with ionospheric remote sensing capability that offers multidimensional observation flexibility and scalability for future cost-effective missions. The first step is to develop a GPS-based HIL simulation testbed for spacecraft formation flying (the Virginia Tech Formation Flying TestBed [VTFFTB]) and validate its performance by benchmarking with past closed-loop real-time reference simulation results. The second step is to propose, demonstrate, and verify an original ionospheric remote sensing technique using GPS-based satellite formation flight for specific scenarios on the VTFFTB. All elements of the testbed are designed for multi-constellation GNSS use; however, the example results shown in this paper are limited to GPS only. The VTFFTB will ultimately become a new mission incubator for novel ionospheric remote sensing techniques or morphology observations using GNSS-based spacecraft formation flying. The remainder of this paper is organized as follows. Section 2 gives a system overview of the VTFFTB and some details on each component predominantly from the hardware standpoint. Section 3 describes the principle and methodology of the navigation and control system and the ionospheric remote sensing system. Section 4 demonstrates the application of the VTFFTB to ionospheric remote sensing from scenario designs to measurement result analysis. Finally, Section 5 summarizes the current accomplishments, issues, and future work. Appendix A presents and analyzes the software and HIL simulation results, as well as validates the performance of the developed VTFFTB by comparing to the results of Park et al.20 ## 2 OVERVIEW OF THE VIRGINIA TECH FORMATION FLYING TESTBED The VTFFTB is a HIL testbed developed for GNSS-based spacecraft formation flying mission designs and closed-loop simulations. The current VTFFTB mainly consists of GPS and Galileo RF hardware signal simulators, a multiconstellation multi-frequency GNSS receiver, a navigation and control system, an Analytical Graphics, Inc. (AGI) Systems Tool Kit (STK) visualization system, and an ionospheric remote sensing system. The system configuration of the VTFFTB is shown in Figure 1. As the core software framework, the navigation and control system plays the central role in managing and executing all simulations. The navigation and control system is primarily composed of a GNSS measurement system, a flight control system, and a remote control system. These consist of a package of scripts and functions written in MATLAB and executed on a laptop computer during simulation. The STK visualization system can be run on the same computer together with the navigation and control system, or on another separate workstation for spacecraft trajectory real-time display or post-simulation replay. Two LEO satellites (hereon referred to as a “chief” and a “deputy”) were simulated on the VTFFTB in this work. Due to only one GNSS receiver being used for spacecraft navigation during this study, the GNSS data and navigation solutions of the chief satellite during the scenario were recorded before the closed-loop real-time simulation, similar to the circumstance in Park et al.20 ![FIGURE 1](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F1.medium.gif) [FIGURE 1](https://navi.ion.org/content/67/1/3/F1) FIGURE 1 System configuration of the VTFFTB The GNSS signal simulators emulate and output the GPS and/or Galileo signals that would be received by the spacecraft in a customized scenario that can simulate the signal delay caused by TEC and amplitude scintillation effects on the RF signals received by a receiver. The RF signals from the RF generator(s) are fed into the GNSS receiver through a coaxial cable. The receiver can constantly produce raw GPS measurements including time, pseudorange, carrier-phase, carrier-to-noise density (C/N), etc and report ephemeris of all tracked GPS satellites. The GNSS measurement system continuously acquires data with respect to the deputy spacecraft and processes them in real time along with the pre-recorded GPS data from the chief spacecraft’s receiver. An EKF is implemented in a fixed rate to estimate the relative state (eg, position and velocity) between the chief and the deputy spacecraft by combining both model-based predictions and differential carrier-phase measurements. The estimations of the relative state are sent to the flight control system to calculate the required thrust of the deputy spacecraft maneuver. Next, the updated dynamics information is transmitted to the remote control system to generate motion commands. Finally, the remote control system delivers the spacecraft’s motion command through TCP/IP back to the GNSS simulator(s) and propagates the spacecraft in the scenario for the upcoming time step, thereby closing the loop of the HIL spacecraft formation flying simulation. The closed-loop configuration (shown in Figure 1 inside the green dashed box) described previously involves real-time maneuver control on the deputy spacecraft during the simulation, which mimics the data flow with feedback in real-world circumstances. After each closed-loop simulation, the collected GPS measurement data are utilized in the ionospheric remote sensing system to further analyze the TEC and scintillation observations for a specific ionospheric scenario. ### 2.1 GNSS simulators The Spirent GSS8000 GNSS RF signal generator and its operational computer with proprietary software are used to emulate and output GPS signals (L1, L2, and L5). The primary simulation software on the operational computer, known as “SimGEN” or “PosApp,” is synchronized with the hardware simulator via TCP/IP to manage and customize a simulation scenario through a graphical user interface (GUI). Users can configure a wide range of scenario features in SimGEN such as vehicle type, GNSS constellation, RF output power level, and model ionospheric delays (via TEC), scintillations, tropospheric delays, etc in the scenario. The “SimGEN Remote” capability of SimGEN allows real-time external motion controls (eg, via TCP/IP sockets) for the vehicle in a scenario, which enables closed-loop HIL simulation. After each simulation, each vehicle’ motion data are recorded to give the “true” states, which will be compared with GNSS measurements or filter estimations to evaluate errors. The GSS8000 is chosen to develop the VTFFTB because its capabilities meet the objectives of the work very well. Also, the Spirent GNSS simulators have been used by other groups for spacecraft formation flying testbed development (eg, NASA,10 DLR,14 University of Calgary,17 Yonsei University,20 etc). ### 2.2 GNSS receiver A NovAtel multi-frequency multi-constellation GNSS receiver is used in the VTFFTB. The OEM628 is a 120-channel receiver with the capability of tracking all current and upcoming GNSS constellations satellite signals including GPS, GLONASS, Galileo, BeiDou, and QZSS. By default, L1C/A, L2P (codeless/semi-codeless), and L5 GPS signals from the GNSS simulators are tracked during HIL simulations. With a maximum data logging rate of 20 Hz, this receiver is able to provide GNSS observables and navigation solutions for LEO spacecraft. Another NovAtel multi-frequency multi-constellation GNSS receiver, GPStation-6, is used for benchmarking purposes to support the development of the ionospheric remote sensing system as it is equipped with commercial TEC and scintillation measurement capabilities. Note that NovAtel OEM receivers have also been used by other groups for HIL formation flying simulation (eg, University of Calgary,17 University of Toronto,18 Stanford University,19 etc). ### 2.3 STK visualization subsystem Developed by Analytical Graphics, Inc., STK is a physics-based software package widely used for aerospace simulation, modeling, and analysis. STK offers a user-friendly GUI interface to design and customize a scenario that includes satellites, aircraft, and many other facilities or objects. After each HIL simulation is completed, the motion data of chief and deputy satellites are loaded into STK to reproduce the spacecraft’s simulated scenario trajectories. In addition, a scripting interface in STK called Connect allows the integration between MATLAB and STK via TCP/IP using the STK’s Object model. The motion data of both chief and deputy satellites can be delivered from MATLAB to STK in real time, which enables satellite trajectory display on STK in real time. ## 3 SOFTWARE ALGORITHMS AND DESIGN METHODOLOGY ### 3.1 Navigation and control system overview The flowchart in Figure 2 highlights the algorithms of the navigation and control system. After the start and initialization, all the communications between different components in the HIL system are established, and the simulated scenario is activated to run. In order to execute the task loop (inside the red dashed box in Figure 2) in real time and synchronize the software with the hardware systems, a timer class (object) is implemented to schedule the dynamics propagation, remote control, GNSS measurement, state estimation, and flight control system in a specified constant rate. The default fixed rate is 1 Hz during the simulations in this work. The first task inside the loop is satellite orbit propagation, which is followed by remote control of the deputy spacecraft motion in the SimGEN scenario. Next, a GNSS measurement system constantly fetches data from the receiver. If GNSS measurements are available, ie, at least one GPS space vehicle (SV or PRN) is tracked, another criterion (at least one common GPS satellite is tracked with qualified data) will be judged to tell if EKF estimation is currently available or not. As long as both judgements are valid, the differential GNSS measurements will be performed and used to optimize the state estimation based on EKF. When there is no available measurement or estimation, the predicted relative states will be directly used as the estimated state to bypass the EKF. The estimated relative state will then be fed into the flight controller, and afterwards, the thrust to drive the deputy spacecraft is output according to the desired relative orbit. If the scenario time is less than the preset total time, the thrust will be saved in a buffer to perform orbital maneuver of the deputy satellite in the next time step. When the scenario is completed, the task loop will be terminated and necessary data will be returned and saved. The design methodologies of each critical subsystem are described as follows. ![FIGURE 2](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F2.medium.gif) [FIGURE 2](https://navi.ion.org/content/67/1/3/F2) FIGURE 2 Key algorithm flowchart of the navigation and control system ### 3.2 Dynamics propagator The dynamics model for satellite trajectory propagation in LEO can be expressed as21 ![1][1] 1 where ![Graphic][2] is the position vector of the satellite in the Earth-centered inertial (ECI) frame, μ is Earth’s gravitational parameter (3.986 × 1014 m3 · s−2), and ![Graphic][3] is the additional acceleration due to all perturbing forces. The first term on the right hand side of Equation (1) is based on Newton’s law of universal gravitation (ie, 2-body model), while the second perturbation term can be expanded as ![2][4] 2 where ![Graphic][5] is the acceleration due to Earth’s non-spherical gravitational potential, ![Graphic][6] is the atmospheric drag acceleration, ![Graphic][7] is the acceleration due to solar radiation pressure, ![Graphic][8] is the acceleration due to third-body (eg, Sun or Moon), and ![Graphic][9] is the thrust acceleration due to an active onboard control system. An 8 × 8 subset of the joint gravity model 3 (JGM-3) was selected to model Earth’s gravitational field, which is written into a MATLAB function based on an eighth-order spherical harmonics function.22 It can be flexibly extended to a higher order gravity model with the cost of demanding more computation power. Testing showed that the run time of an 8 × 8 JGM-3 model propagator is around 0.02 seconds (on a Dell quad-core laptop), while the run time of a 2-body model propagator is only about 0.002 seconds, which is 10 times faster compared to the 8 × 8 JGM-3 model. Atmospheric drag, solar radiation pressure, and third-body effects were modeled during some software simulations, but not adopted during any HIL simulations in order to optimize the computation speed of the dynamics propagator in real time. When ignoring these perturbations, the orbit propagation fidelity is reduced compared to real-world scenarios. However, these perturbation impacts are relatively small during a simulation lasting only a few hours and do not significantly affect the overall results. A 1-Hz pulsed plasma thruster (PPT) system with a scalable maximum unit thrust cap is used to maneuver the deputy spacecraft. Note that no active control is applied on the chief satellite during the simulation scenarios. It should be noted that reliable propagation of satellite orbits for the current version of SimGEN requires that the motion be consistent with the Spirent’s use of a first-order Euler method in the numerical motion propagation.23 This was discovered when we first implemented a Runge-Kutta fourth-order (RK4) integrator, which resulted in propagation inconsistencies flagged by a “kinetic error” in SimGEN. Coordinate transformations between the ECI (J2000) and the ECEF (Earth-Centered, Earth-Fixed) frames are repetitively performed during the dynamics propagation. The formation flight controller is executed under a body-frame coordinate system (ie, Local-Vertical, Local-Horizontal) that requires coordinate transformations between the ECI and body-frame systems in each time step. The body frame is a coordinate system originating at the center of mass of a spacecraft (usually the chief) that is often used to describe the relative dynamics of spacecraft formation flying. As shown in Figure 3, the unit vector ![Graphic][10] is the radial direction and is directed from the center of Earth to the chief satellite. The unit vector ![Graphic][11] is referred to as the “in-track” or “along-track” direction and is directed along the velocity vector of the chief satellite. The unit vector ![Graphic][12] is the “cross-track” direction and is perpendicular to the orbital plane of the chief satellite with the positive direction along with the angular momentum vector. ![FIGURE 3](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F3.medium.gif) [FIGURE 3](https://navi.ion.org/content/67/1/3/F3) FIGURE 3 Body frame of satellite formation flying ### 3.3 Differential GPS measurement The single-differential carrier-phase measurement technique used by Busse et al11 is applied in the VTFFTB GNSS measurement system to estimate the relative state between the chief and deputy satellite. For small-baseline formation flying, some GNSS ranging errors (eg, ionospheric delays, troposphere delays, and GPS satellite clock bias) could be canceled out or neglected by this approach. The carrier-phase of GPS satellite *m* with respect to vehicle *i* is expressed as ![3][13] 3 Here, *CP* is the carrier-phase measured in number of cycles, *λ* is the wavelength of the current band of GNSS signals, *R* is the geometric range to the satellite (true distance without error) in meters, *c* is the speed of light in a vacuum, *δT* is the satellite clock bias in seconds, *δt* is the receiver clock bias in seconds, *I* denotes ionospheric signal delay in meters, *T* is the troposphere signal delay in meters, *M* is the multipath effect delay in meters, *D* denotes errors caused by the geometric dilution effect in meters, *W* denotes other noise effects (including receiver noise, relativity effects, etc), and *N* is the carrier-phase measurement (integer) ambiguity. After consolidating the multipath delay and geometric dilution effect into the noise term *W*, the single difference *CP* with respect to the GPS satellite *m* between the chief spacecraft *A* and deputy spacecraft *B* can be approximated as ![4][14] 4 The relative position is defined as the difference between the absolute position of chief and deputy satellite, written as ![Graphic][15]. Equation (4) can then be rewritten into the measurement matrix form shown in Equation (5), after changing the units from number of cycles into meters and assuming *n* satellites are available at the time of measurement: ![5][16] 5 On the right hand side of Equation (5), the positions of all visible GPS satellites ![Graphic][17] are computed by the GPS ephemeris from a navigation file. The chief satellite position ![Graphic][18] is modeled by a dynamics propagator. The relative position ![Graphic][19] and differential receiver clock bias Δδ*t**AB* are included in the state matrix to be estimated. The differential noise term ![Graphic][20] is modeled as white Gaussian noise (WGN). Based on first available observables, initial estimation is required to determine ![Graphic][21] (assuming all of ![Graphic][22] are constant integers). ### 3.4 EKF-based estimation An EKF is used for relative state estimation in the VTFFTB navigation and control system. A detailed explanation of this method can be found in Yaakov et al24; therefore, an abbreviated description will be provided here. The differential GPS measurements are incorporated into the EKF together with the dynamics propagations. The state to be estimated can be defined as an 8 × 1 matrix: ![6][23] 6 Initialization of the state to be estimated *X*(t), the state covariance *P*(t), the process noise covariance *Q* (assumed as a constant), and the measurement noise covariance *R*(t) are required to begin the filtering. After initialization, the relative state is predicted at a fixed time step (Δ*t*) forward via the dynamics propagator *f*: ![Graphic][24], where the superscript “*P*” denotes “prediction.” Then the nonlinear relative dynamics is linearized, which yields the state prediction Jacobian ![Graphic][25]. The continuous linearized Jacobian matrix can then be discretized as ![Graphic][26]. Next, the state prediction covariance can be determined using *P**P*(*t*1)= Φ(*t**o*)*P*(*t**o*)Φ*T*(*t**o*)+*Q*. If the estimation judgement is valid, the measurement matrix of differential carrier-phase *Z*(*t*1) can be formed, and the value of *n* can be determined at this moment. Meanwhile, the modeling measurement matrix can be predicted by equating ![Graphic][27] to the right hand side of Equation (5). Now that both predicted measurement and actual measurement are ready, the measurement residual can be calculated as ![Graphic][28]. Subsequently, the nonlinear observation is linearized by evaluating the measurement prediction Jacobian ![Graphic][29]. Following that, the residual covariance can be determined using *S*(*t*1) = *R*(*t*1)+*H*(*t*1)*P**P*(*t*1)*H**T*(*t*1). The measurement noise covariance matrix *R* is updated in every step based on raw GPS measurements through receiver logging. Standard derivations of carrier-phase were used to construct the diagonal matrix *R* by taking a square and then multiplying by 2. Here, the standard deviations are a function of the geometry (ie, dilution of precision) and the estimated precision of the measurements. The values are the square roots of the diagonal terms of the covariance matrix that relate to the position under estimation. At this point, the Kalman Filter gain can be given as *K*(*t*1) = *P**P*(*t*1)*H**T*(*t*1)*S*−1(*t*1). Finally, the state estimation can be updated as *X*(*t*1) = *X**P*(*t*1)+*K*(*t*1)Δ*Z*(*t*1). The state covariance is updated using the Joseph form covariance: *P*(*t*1) = [*I*8 × 8 − *K*(*t*1)*H*(*t*1)]*P**P*(*t*1)[*I*8 × 8 − *K*(*t*1)*H* (*t*1)]*T*+*K*(*t*1)*R*(*t*1)*K**T*(*t*1), where *I*8 × 8 is an 8-by-8 identity matrix. ### 3.5 Flight controller In order to achieve the desired satellite formation, a flight controller is designed and implemented to maneuver the deputy spacecraft using only thrust of three unit-vector directions (radial, in-track, or cross-track). To solve the required thrust given the intended relative state, the VTFFTB’s flight controller is based on the State Dependent Riccati Equation (SDRE) technique25 and implemented by Park et al20 for relative orbit control. The SDRE-technique-solved equations can be expressed in the state-dependent coefficient (SDC) form ![Graphic][30], where the state matrix *A* in this work corresponds to the Hill-Clohessy-Wiltshire equations of relative motion (![Graphic][31]; orbital eccentricity ≈ 0), *B* is the control input matrix, ![Graphic][32] is the state of relative position and velocity expressed in the body frame ![Graphic][33], and ![Graphic][34] is the thrust expressed in the body frame to be solved. The solution to this equation can be cast as an optimization problem with a performance index *J*, written as ![7][35] 7 where ![Graphic][36] is the desired relative state vector. The weighted control matrices *Q**ctr* and *R**ctr* are written as the identity matrix multiplied by a constant. This is where the tradeoff between the energy required to control the formation flying and the errors associated with the current formation flying configuration is considered. Using the SDRE control law, the optimization problem here can be solved by finding *G* in the following expression: ![Graphic][37]. *G* here is the positive-definite solution of the Riccati equation, ![Graphic][38], which can be found by using the “care” function in MATLAB. As long as the Riccati equation is solved for the linear quadratic regulators, the cost function (*J*) can be minimized. Further, the thrust limit of this real-time feedback control system in each unit-vector direction can be customized. A thrust cap (ΔVmax) of 1 m/s2 in all unit directions was set as default. ### 3.6 Ionospheric remote sensing system A multi-frequency GNSS ionospheric remote sensing system is developed to take advantage of the growing GNSS network for space weather studies. By logging and processing raw GNSS measurements through the OEM628 receiver, a package of MATLAB software has been developed to estimate the GNSS TEC and amplitude scintillation index S4. Only GPS is considered in this study, however, the framework allows multi-constellation TEC and S4 measurements. The software algorithm of the system is outlined in Figure 4. Raw pseudorange, carrier-phase, ephemeris, and carrier to noise density (C/N) data with a maximum logging rate of 20 Hz are extracted from the observation files generated by the receiver. Next, raw pseudorange relative slant TECs and raw carrier-phase relative slant TECs along with the elevation data are computed and saved. Here, “relative” means a bias due to “differential-code-biases (DCB).” Due to receiver firmware limitations, S4 is calculated by taking the ratio of the standard deviation of C/N to the average of C/N over a 1-minute period. This yields a “pseudo” S4, which is highly correlated to the S4 computed via S/N instead of C/N0. This does not compromise the goal of identifying scintillation and detecting ionospheric irregularities.9 The fitted relative slant TECs can be estimated based on a least-squares fitting method. To get the bias free absolute slant TECs, the satellite and receiver DCB must be properly estimated and eliminated. The TEC bias estimation is based on the least squares method developed by Gaposchkin and Coster at MIT Lincoln Laboratory,26 and this method has been implemented for DCB estimations in the MIT Madrigal GPS TEC database.27 Finally, a mapping function is applied on the elevation data in order to convert the absolute slant TECs into vertical TECs. To validate the algorithm for TEC and S4 measurements using raw data taken from the OEM628 receiver, the GPStation-6 receiver is utilized to benchmark the measurement results. As the OEM628 and the GPStation-6 were both connected to either the same RF output of the GSS-8000 simulator or the rooftop GNSS antenna located at Blacksburg, VA (37.2°N, 80.4°W), the GPStation-6 receiver was used to directly generate TECs and S4 to benchmark with the TECs and S4 measured and processed from OEM628. More details of the system validation and benchmarking work can be found in Peng.28 ![FIGURE 4](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F4.medium.gif) [FIGURE 4](https://navi.ion.org/content/67/1/3/F4) FIGURE 4 Flowchart of the ionospheric remote sensing system ### 3.7 Testing and validation of real-time formation flying A reference formation flying scenario similar to the one described by Park et al20 has been established to test and validate the functionality of the VTFFTB for closed-loop real-time formation flying. The initial classical orbital elements and equivalent ECEF epoch state vectors are provided by Park et al.20 The software and hardware simulations were both run for 1 hour with the same set of initial conditions for benchmarking. The reference scenario is composed of two LEO satellites with an initial in-track separation of 1000 m and an altitude of approximately 550 km. This scenario is designed to perform the formation acquisition required to maneuver the deputy spacecraft into an in-track offset of 100 m with no other radial or cross-track separation and maintain such formation flying until the end of the simulation. Only the deputy satellite was under active control and the observation file of the chief satellite was pre-recorded. This is described in more details in the Appendix. ## 4 IONOSPHERIC REMOTE SENSING DEMONSTRATION ### 4.1 Scenario design In this section, an Equatorial Spread F (ESF) scenario is designed to demonstrate the application of GPS-based spacecraft formation flying to observation of ionospheric irregularities. ESF is an ionospheric irregularity that forms in the post-sunset equatorial region of the ionosphere, which is often observed as plasma bubbles or plumes.8 A bottomside spread F event between 250 and 350 km in altitude before 23:00 local time (LT) was observed by Rodrigues et al.29 A simplified version of this ESF observation was simulated in a SimGEN scenario on the VTFFTB. Due to limitations in SimGEN’s current capabilities, horizontal variations in TEC are not supported by default and are instead emulated by modeling amplitude scintillation in a grid with minimum 10°latitude versus 1-hour LT. Vertical ionospheric irregularities can be simulated by user-defined vertical TEC profiles. A sufficient number of data points can be programmed to ensure appropriate resolution of the ionospheric profile. The ionospheric model is setup/programmed into SimGEN before running each HIL simulation. A homogeneous irregularity region of a plasma bubble was created in a cuboid. This region ranges from 290 to 350 km vertically, 20:00 LT to 21:00 LT longitudinally, and 0° to 10° S latitudinally. For the purpose of this study, it is assumed that this irregularity region appears near the Jicamarca Radio Observatory (JRO) [11°38′S, 77°02′W] in Peru starting at 01:00:00 (GPS time) on March 21, 2016. Using the TIE-GCM model,30 the local vertical electron density profile and the corresponding TEC profile of this region were designed as shown in Figure 5. The depletion of electron density between 290 and 350 km emulates a simplified plasma bubble vertically. The horizontal boundary of the irregularity region was simulated by setting up a 10° latitude × 15° longitude region with an S4 index of 0.4. A static cuboid region with a fixed S4 value cannot represent realistic scenarios of scintillation effects caused by ESF; however, this is sufficient for the concept demonstration described here. Other robust scintillation software simulators (such as Xu et al31) as well as hardware simulation techniques32 will be implemented to increase the VTFFTB’s simulation capability of ionospheric scintillations impacting LEO satellites in the future. ![FIGURE 5](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F5.medium.gif) [FIGURE 5](https://navi.ion.org/content/67/1/3/F5) FIGURE 5 (A) Vertical electron density profile; (B) corresponding vertical TEC profile Taking advantage of the GNSS receivers on a pair of spacecraft flying in formation, TEC, electron density, and ionospheric scintillation can be probed in a novel way with a high level of resolution. STK and the dynamics propagator in MATLAB are utilized to design the orbits of this formation flying scenario. A low-inclination, small-eccentricity, elliptical LEO is selected to let the spacecraft constantly fly above and monitor the equatorial region for ESF observations. The elliptical feature allows the change of spacecraft altitude, which is critical to measuring the vertical electron density profile of the ionosphere. The initial classical orbit elements of the chief and deputy spacecraft are listed in Table 1. The orbit is designed to achieve and maintain a 1-km negative radial offset and a 15-km negative in-track offset after a short formation acquisition process. This ionospheric remote sensing scenario of two satellites passing through the irregularity region in formation is summarized and illustrated in Figure 6A. View this table: [TABLE 1](https://navi.ion.org/content/67/1/3/T1) TABLE 1 Initial orbit conditions of the remote sensing scenario ![FIGURE 6](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F6.medium.gif) [FIGURE 6](https://navi.ion.org/content/67/1/3/F6) FIGURE 6 (A) Illustration of the ionospheric remote sensing scenario; (B) in-track and radial offset change of the relative orbit HIL simulation of this ionospheric remote sensing scenario has been performed using the VTFFTB. As shown in Figure 7A, after 98 seconds, the desired radial offset is reached within ±1 m error bound; after about 57 seconds, the desired in-track offset is reached within ±1 m error bound. The initial/final in-track and radial offsets during a HIL simulation are shown in Figure 6B. After a 1-hour HIL simulation with ΔVmax of 1.0 m/s2, the total thrust consumed in radial, in-track, and cross-track direction are 27.95 m/s, 6.95 m/s, and 0.61 m/s, respectively. As shown in Figure 7B, thrust was applied to correct the relative orbit during the early phase in both radial and in-track directions. A small amount of radial thrust was used to maintain the radial offset throughout the rest of the scenario. ![FIGURE 7](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F7.medium.gif) [FIGURE 7](https://navi.ion.org/content/67/1/3/F7) FIGURE 7 HIL simulation results of the remote sensing scenario: (A) relative position history; (B) thrust history Both chief and deputy spacecraft can continuously monitor the TEC and S4 above them. By using the vertical TEC measured by two spacecraft at different altitudes in formation with onboard GNSS receivers, the vertical electron density between the spacecraft can be estimated, and then the vertical boundary of the plasma bubble can be detected. When the spacecraft cross the horizontal edge of the ESF region, a gradient of TEC or electron density can be expected in the observation data. ### 4.2 Ionospheric measurements from HIL simulations The ionospheric remote sensing system was used to process the TEC and S4 measured by both the chief and deputy spacecraft during the scenario. GPS PRN 7 and PRN 9 were tracked by both the chief and deputy spacecraft, and their vertical TEC (using GPS L1 and L2) are plotted in Figure 8. A nearly-constant TEC region can be seen in the chief measurements from about 30.17 TECU (at 1957 seconds) to 30.10 TECU (at 2165 seconds) with respect to both GPS PRN 7 and PRN 9, and a nearly-constant TEC region can also be seen in the deputy measurements from about 30.17 TECU (at 1961 seconds) to 30.10 TECU (at 2169 seconds) relative to both GPS PRN 7 and PRN 9. Both chief and deputy satellites were rising in altitude during that period, and these results indicate both spacecraft crossed the vertical boundaries of the plasma bubble region. This, given a vertical speed of approximately 295.79 m/s for each spacecraft, yielded a 1183.16 m vertical separation, which is fairly consistent with the true vertical difference of 984.90 m. Note that there is a data gap in Figure 8C,D due to a loss of lock of PRN 9. ![FIGURE 8](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F8.medium.gif) [FIGURE 8](https://navi.ion.org/content/67/1/3/F8) FIGURE 8 Vertical TEC measurements: (A) chief using PRN 7; (B) deputy using PRN 7; (C) chief using PRN 9; (D) deputy using PRN 9 As shown in Figure 9, vertical electron density (Ne) profiles are derived based on a differential-TEC approach by using the TEC measurements in Figure 8. The modeled Ne is plotted in red, and the measured Ne is plotted in blue. In Figure 9A, measurements from GPS PRN 7 are consistent with the model, but the measured Ne of the plasma bubble region is noisy. This could be caused by signal scintillation with respect to the particular communication link from the formation satellite(s) to the GPS PRN 7. In Figure 9B, measurements from GPS PRN 9 are consistent with the model and less noisy. To avoid potential measurement inaccuracy, TEC biases must be carefully estimated. Note that the number and geometry of GPS PRNs being tracked do not substantively affect the Ne retrieval accuracy due to the simplifying assumption of 1D spatial variations in the ionospheric model. ![FIGURE 9](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F9.medium.gif) [FIGURE 9](https://navi.ion.org/content/67/1/3/F9) FIGURE 9 Retrieved vertical electron density profiles using (A) PRN 7 and (B) PRN 9 Figure 10A shows the S4 indices calculated by the 20 Hz C/N (GPS L1 band) measurements of the chief satellite versus altitude. Measurements by all 27 tracked SVs are shown in the same plot, which reveals the amplitude scintillation is more severe between 290 and 350 km. As the time resolution of S4 computation is 1 minute, the corresponding vertical resolution is about 17.7 km. The vertical electron density is more appropriate to be used to look at the exact location of the vertical boundary in terms of the spatial resolution. However, the S4 index measurement can increase the confidence for the vertical boundary detection as a sanity check. The horizontal boundary can be estimated by analyzing the geographic distribution of GPS L1 S4 magnitude as shown in Figure 10B. The chief satellite path is projected on the latitude-longitude plane. S4 measured every minute by different PRNs is plotted vertically on each location along the path. The scintillation events (when S4 > 0.2) are concentrated in a region starting at [7.1°S, 82.7°W] and ending at [8.8°S, 67.6°W], consistent with the S4 cell region modeled in SimGEN. The 1 Hz C/N measurements can be examined to further confirm the horizontal boundaries. Data analysis showed that the fluctuation of C/N observed by the chief satellite is about 2 seconds ahead of that observed by the deputy satellite. This is consistent with the 15 km in-track separation during the formation keeping and a horizontal speed of 7-8 km/s. ![FIGURE 10](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F10.medium.gif) [FIGURE 10](https://navi.ion.org/content/67/1/3/F10) FIGURE 10 GPS L1 S4 measurements by the chief satellite: (A) vertical distribution; (B) horizontal distribution ### 4.3 Validation using a sounding rocket measurement profile To evaluate the resolution capability of the vertical Ne measurement technique described earlier, an ESF Ne profile measured from a past sounding rocket (SR) campaign33 was simulated in the VTFFTB. The PLUMEX I sounding rocket was launched on July 17, 1979, from Roi Namur Island of the Kwajalein Atoll (9.40°N, 167.48°E, ~9° magnetic dip), where an intense ESF event was taking place. The upleg Ne profile measured by the fixed-bias Langmuir probe on PLUMEX I is used to derive the vertical TEC profile emulated in a scenario. The 3-hour simulation of this scenario was conducted with the same initial orbit conditions shown in Table 1. Figure 11 shows the sounding rocket Ne measurement plotted with symbol “×” and labeled as “SR Ne.” The VTFFTB derived vertical Ne using GPS L1 and L2 TEC measurements from PRN 24 are plotted with symbol “ ” and labelled as “VTFFTB Ne.” The discrepancies between “SR Ne” and “VTFFTB Ne” at each altitude are plotted with symbol “+” and labeled as “Difference.” The average absolute discrepancy across the whole altitude range (from ~236 km to ~567 km) is 6.56 × 104 cm−3. Relatively good agreement is seen for the Region 1 bubble (between ~480 km and ~520 km) and the Region 3 bubble (between ~340 km and ~380 km), with the average absolute discrepancy of 5.28 × 104 cm−3 and 6.15 × 104 cm−3, respectively. The Region 2 bubble (between ~395 km and ~435 km), given a relatively larger average absolute discrepancy of 7.86 × 104 cm−3, is still detectable. Overall, the results show that the differential GPS TEC technique with satellite formation flight can be used to retrieve the vertical Ne profile from realistic ESF observations. This is a stepping stone to the fully integrated formation flying observation of ionospheric structures, which further demonstrates the practicability of applying this remote sensing technique for future ionospheric probing missions. Work to consider the ability to resolve ionospheric irregularities with much smaller density fluctuations is ongoing. ![FIGURE 11](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F11.medium.gif) [FIGURE 11](https://navi.ion.org/content/67/1/3/F11) FIGURE 11 Comparison between VTFFTB simulated and sounding rocket33 electron density Ne measurement profiles. ESF plasma density bubbles are indicated by 1, 2, and 3 ### 4.4 Formation configuration sensitivity study The VTFFTB can serve as a mission incubator for different formation configurations. The formation configuration with the 1-km radial offset and 15-km in-track offset demonstrates its practicality to sense an ESF region and generate localized vertical Ne retrieval. Three more formation configurations with different radial offsets were simulated in order to study the sensitivity of vertical Ne retrieval in comparison to the baseline configuration described in the previous subsection. In addition to the baseline 1-km radial offset, formations keeping with 100-m, 3-km, and 10-km radial offsets were simulated on the VTFFTB. Using the PRN7 results as an example comparison, the vertical Ne measurements from the four different radial offsets are presented in Figure 12. The formation configuration with a 100-m radial offset, as shown in Figure 12A, exhibits the most scattering behavior in the Ne measurements. The overall shape is consistent with the trend of the SR vertical Ne profile. However, the plasma bubbles are difficult to identify without further applying any signal processing techniques (eg, low-pass filter). As a comparison, the plasma bubbles can be clearly recovered from the vertical Ne retrieved with the configurations of 1-km, 3-km, or 10-km radial offset. The measurements in Figure 12B-D exhibit much less scattering compared to the 100-m case (Figure 12A). It can be observed that as the radial offset increases, the fine scale spatial structures are not well resolved. It appears from Figure 12 that the 1-km and 3-km radial offset cases (Figure 12B,C) produce the most optimal vertical Ne retrieval profile. ![FIGURE 12](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F12.medium.gif) [FIGURE 12](https://navi.ion.org/content/67/1/3/F12) FIGURE 12 GPS Ne measurement from HIL simulations using PRN 7 with different radial offsets: (A) 100 m; (B) 1 km; (C) 3 km; (D) 10 km Electron density wavenumber spectra were calculated to investigate, in more detail, the resolution of various spatial scales depending on the altitude separation between the chief and the deputy. The wavenumber *k* is related to the spatial scale λ by *k* = 2π/λ. Spatial Fourier Transforms of Ne, ie, Ne(*k*), were computed for the four radial offset results and the SR Ne model in Figure 12. Figure 13 shows the wavenumber spectra, which are defined here as the square of the absolute value of Ne(*k*), ie, |Ne(*k*)|2. The equatorial plasma bubbles (EPBs) can be clearly seen in the spectra for 0.1 < *k* < 0.2, which is equivalent to the spatial scale range 30 km < λ < 60 km. Therefore, the separations of 1-km, 3-km, and 10-km resolve the EPBs quite well. As expected, the EPBs can at least be detected in the raw Ne measurement using a radial offset of 100 m; however, the noise level is too high to accurately retrieve the smaller scale Ne structures. The 1-km and 3-km radial offsets can resolve the smaller scale Ne structures fairly well. However, it appears that 1 km is the more optimal configuration in this scenario. Note that for smaller scale Ne structures, the 3-km spectrum starts to deviate from the SR spectrum for *k* > 1, which corresponds to a spatial scale λ < 6 km. The 10-km spectrum starts to deviate from the SR spectrum for *k* > 0.3, which corresponds to a spatial scale λ < 20 km. ![FIGURE 13](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F13.medium.gif) [FIGURE 13](https://navi.ion.org/content/67/1/3/F13) FIGURE 13 Wavenumber spectra comparison of Ne retrieval with different formation radial offsets This comparison indicates that vertical Ne retrieval characteristics are sensitive to the altitude separation or radial offset of the two-satellite formation configuration. When the altitude separation between two LEO satellite/receivers is too small, the measured Ne results become noisier, and the measurement accuracy is reduced. This is due to the measured VTEC noise level becoming close to the corresponding VTEC difference between the two different altitudes. When the altitude separation between two LEO receivers becomes sufficiently large, this effect is reduced, the measurement becomes less noisy, and the Ne profile is more easily resolvable. However, for too large a separation, the Ne resolution is expected to be reduced for those ionospheric structures with a spatial scale less than twice the altitude separation (as seen in Figure 13). An analysis based on a simplified model corroborates the above stated concept on how the Ne measurement sensitivity depends on the formation height separation and the TEC noise level. The model adds WGN with a signal-to-noise ratio of 75 dB to the discrete VTEC data derived from the SR Ne profile. Given a fixed height separation between two satellites, the altitude of both satellites is slid across the SR vertical VTEC profile, which mimics sampling VTEC at different altitudes during a formation flying scenario. By applying the vertical Ne retrieval algorithm with the four different height separations, comparable Ne retrieval characteristics are obtained as in Figure 12. By increasing the WGN level, the Ne measurement scattering level is increased accordingly. Therefore, the reduction in resolution at large vertical separation is a lack of spatial sampling according to the Nyquist sampling theorem. The reduction in resolution at small vertical separation is essentially a measurement noise effect. In summary, the sensitivity study in this subsection demonstrates that the VTFFTB can be utilized to determine optimal formation configurations for vertical Ne retrieval depending on the TEC noise level and spatial resolution requirement. ## 5 SUMMARY AND FUTURE WORK The goal of this work was to develop the VTFFTB, a GNSS-based HIL simulation testbed for spacecraft formation flying with the capability of ionospheric remote sensing mission designs and demonstrations. The current infrastructure of VTFFTB is composed of GNSS RF hardware simulators, two multi-constellation multi-frequency GNSS receivers, a navigation and control system, an STK visualization system, and an ionospheric remote sensing system. A GNSS receiver was utilized to receive emulated GPS (L1, L2, and L5) signals from a GPS simulator and provide real-time measurements for navigation and ionospheric remote sensing of a pair of LEO satellites. The navigation and control system was developed to perform real-time dynamics propagation, EKF-based relative orbit estimations using single-differential carrier-phase measurements, and maneuver commands for the deputy spacecraft using the SDRE control technique. The mission visualization system using STK allows real-time trajectory display or replay of spacecraft formation flying simulations. A reference test scenario of two LEO spacecraft was established with the initial in-track separation of 1000 m and simulated successfully to achieve a desired leader-follower configuration of a 100-m along-track offset (see Appendix A). Both software and HIL simulations were performed and analyzed for this 1-hour reference scenario, and the overall results are consistent with past reference simulation results (see Appendix A). Additional software simulation shows the thrust cap is sensitive to the relative orbit history and total thrust needed to achieve the desired formation (see Appendix A). This has validated the functionalities of the VTFFTB for GPS-based HIL spacecraft formation flight simulations. Together with the demonstrated HIL simulation capability, the ionospheric remote sensing system has been applied to ionospheric tomography in special scenarios simulated on the VTFFTB that include specific ionospheric phenomenon. A LEO spacecraft formation flying mission to observe two ESF (plasma bubbles) scenarios has been designed, illustrated, and simulated on VTFFTB. Measurement results from the HIL simulation demonstrated that GPS-based formation flight of two satellites is able to probe the vertical electron density and study the irregularity structure. Specifically, a vertical boundary can be determined by examining the measured vertical electron density profile or S4, while a horizontal boundary can be investigated by TEC gradient or analyzing the direction of scintillation event(s). It has been shown that the VTFFTB can be used to optimize formation configurations for vertical Ne measurements given TEC noise level and ionospheric structure spatial scales. During the HIL simulations, a small amount of relative orbit deviation was caused by poor initial EKF estimations. Additional work is required to improve the initialization of the EKF estimation. A higher-fidelity dynamics model will be implemented to propagate the satellite orbit. The carrier-phase GPS measurement model can be upgraded from single-differential to double-differential. Further, the Spirent GSS7800 Galileo simulator will be utilized to incorporate the rapidly growing Galileo constellation, which can potentially increase both the navigation accuracy and remote sensing resolution. The ionosphere environment modeled in the current HIL simulation is simplified to a 1D spatial variation. Higher fidelity ionospheric modeling will be incorporated in the future. Additional GNSS receivers will be incorporated to track both the chief and deputy satellites in real time and simulate multi-spacecraft formation flying scenarios for higher observational robustness. Autonomous maneuvers in response to real-time ionospheric conditions are significantly advantageous to observe multi-scale ionospheric structures with higher temporal and spatial resolution. The limited fuel budget of small satellites requires development of appropriate control algorithms to enable such ionospheric-based feedback maneuvering. Additional space weather information (eg, seasonal or real-time ionospheric profiles) can facilitate the execution of such autonomous ionospheric tracking; however, this is beyond the scope of current work and will be considered in the near future. The flexible infrastructure of VTFFTB is capable of supporting the design and development of more novel upper atmospheric missions using spacecraft formation flying. ## HOW TO CITE THIS ARTICLE Peng YX, Scales WA, Edwards TR. GPS-based satellite formation flight simulation and applications to ionospheric remote sensing. *NAVIGATION*. 2020;67:3–21. [https://doi.org/10.1002/navi.354](https://doi.org/10.1002/navi.354) ## ACKNOWLEDGMENTS Thanks to Gregory Earle, Dylan Thomas, and Michael Esswein for providing useful suggestions on this work. ## APPENDIX A #### A.1. Software simulation To test the flight control system, pure software simulations were performed before running the real-time HIL simulation. In the software simulation, motions of the chief and deputy spacecraft were driven by the dynamics propagator and the flight controller. The GNSS measurement system was deactivated, and the EKF was not applied. Figure A1 presents the relative position and thrust history of software simulation based on the 2-body gravity model and a maximum cap ΔVmax of 1.4 m/s2. These results demonstrate that the flight controller is able to maneuver the relative orbit into the target state. At 124 seconds, a 100-m in-track formation is achieved, and the formation keeping phase begins at 125 seconds, as shown in Figure A1A. Total thrust applied for the 1-hour simulation in radial, in-track, and cross-track direction are 2.43 m/s, 75.48 m/s, and 1.07 × 10−5 m/s, respectively. Due to an overshoot, the closest in-track separation during this simulation is 38.17 m (at 56 seconds) during the formation acquisition phase. ![FIGURE A1](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F14.medium.gif) [FIGURE A1](https://navi.ion.org/content/67/1/3/F14) FIGURE A1 Software simulation results: (A) relative position history; (B) thrust history After the deputy overshoots the chief, the controller commands the deputy to apply thrust in the opposite direction of motion, as shown in Figure A1B. During the early portions of the acquisition, the controller commands a large amount of thrust which causes the deputy to require additional thrust in order to counteract the momentum applied during earlier portions of the maneuver. This makes the approach overzealous and a pattern of oscillating relative distance appears in these situations. From the classic feedback control perspective, the overshoot can be understood as an underdamped feedback system that is correctable and sensitive to the thrust limit. Correlations between the thrust limit (cap) and total thrust or the overshoot severity can be expected. In a real scenario, the amount of thrust that the satellite can apply is limited. It is imperative to understand how this affects the control and total amount of thrust applied by the spacecraft. An analysis of the thrust cap is made with changes of the thrust limit being applied during the maneuver for the reference formation-flying scenario. Variable thrust caps (0-10 m/s2) have been applied on the same controller, and the analysis results are shown in Figure A2. Each thrust command was applied once per second. In Figure A2A, the total thrust over the entire simulation versus the thrust limit are plotted in blue, red, and green for the radial, in-track, and cross-track directions, respectively. The minimum in-track offset over the entire simulation against the thrust limit is shown in Figure A2B. The inoperable portion in Figure A2A is a very low thrust cap region (ΔVmax < 0.02 m/s2), which is unrealistic and would not be used for formation keeping and control. The thrust in this inoperable region is too low, and the desired formation is never reached. When the formation acquisition can be reached, the in-track total thrust required to get to the desired location has three distinct regions. Region 1 (0.02 m/s2 < ΔVmax < 1.20 m/s2) is the low damping ratio region due to the relatively weak thruster, where the relative state oscillation would lead to a negative in-track offset as shown in Figure A2B. This could potentially cause collision hazards for the spacecraft and as such is not recommended. Region 2 (1.20 m/s2 < ΔVmax < 1.95 m/s2) is where the local minimum of fuel consumption for the in-track direction can be found for this particular scenario, reaching a local minimum when ΔVmax = 1.4 m/s2. The largest minimum in-track offset (69.02 m) can be found at the right edge of this region. Region 3 (ΔVmax > 1.95 m/s2) is the high damping ratio region due to the relatively high thrust cap. When the thrust cap is very high, the total thrust needed approaches a constant value. Additional analysis shows that most of the thrust is required during the period of formation acquisition, and the thrust needed to keep in formation is quite low and stable. Therefore, extra attention is required to design an optimized controller for optimal fuel consumption and satellite maneuverability. ![FIGURE A2](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F15.medium.gif) [FIGURE A2](https://navi.ion.org/content/67/1/3/F15) FIGURE A2 Sensitivity analysis on thrust cap: (A) total thrust versus thrust cap; (B) minimum in-track offset versus thrust cap #### A.2. HIL simulation A real-time closed-loop HIL simulation has been conducted for the same reference scenario, and the results are highly consistent with the software simulation result in Figure A1. The discrepancies between the real-time HIL simulation results and software simulation results are computed and plotted in Figure A3. During the formation acquisition, a minimum in-track distance of 38.17 m was reached at 56 seconds due to an overshoot. This minimum separation is smaller than in Figure 7 in Park et al20 due to the utilization of different controller parameters. Again, the in-track separation of 100 m is achieved at 124 seconds, and the formation keeping phase begins at 125 seconds. Here, the formation acquisition is achieved faster than in Park et al.20 Given a maximum thrust cap of 1.4 m/s2 in each unit direction, the total thrust applied during the 1-hour simulation in the radial, in-track, and cross-track direction are 5.63 m/s, 76.88 m/s, and 1.12 m/s, respectively. Along the in-track direction, a total thrust of 75.45 m/s is required for the formation acquisition phase, whereas only 1.44 m/s is used during the formation keeping phase. During the formation keeping, thrust in all directions became low and relatively stable. This indicates that the desired state of formation flying is fairly fuel-efficient during the formation keeping process. Compared with the software simulation, the HIL simulation required more thrust in all three directions to account for more “real-world” effects in the HIL environment, including the measurement noise, processing noise, estimation errors, etc. Compared to Figure 10 in Park et al,20 more thrust is consumed in this HIL simulation for formation acquisition, which is consistent with a more rapid convergence time. Further analysis revealing the characteristics of the formation keeping process in each body frame unit-direction is summarized in Table A1 (VTFFTB). Compared to the results from Park et al20 as shown in Table A1 here, the formation keeping accuracy of VTFFTB for this scenario is comparable at a submeter level. The maximum error along the in-track direction is 2.8 m, and the mean error is approximately 1 cm. The largest formation keeping deviations in this simulation are caused by the relatively initial estimates of the EKF in the system that led to the sudden thrust at 702 seconds in Figure A3B. In summary, the real-time closed-loop HIL simulation results of a reference scenario GNSS-based spacecraft formation flying on the VTFFTB have been presented and discussed. Relative orbit history shows a successful navigation and control process, and the results are sensitive to the dynamics propagator, state estimations, and controller. As discussed in the software simulation section, thrust history is sensitive to the scenario orbit and the control algorithm, including the thrust cap value. Due to the flawed initial state estimations of the system, the EKF implementation still needs some improvement. Regardless, the overall results have validated the feasibility and performance of spacecraft formation flying simulation on the VTFFTB via benchmarking with the results in Park et al.20 ![FIGURE A3](https://navi.ion.org/https://navi.ion.org/content/navi/67/1/3/F16.medium.gif) [FIGURE A3](https://navi.ion.org/content/67/1/3/F16) FIGURE A3 Discrepancy between HIL and software simulation results in Figure A1: (A) relative position history; (B) thrust history View this table: [TABLE A1](https://navi.ion.org/content/67/1/3/T2) TABLE A1 Position difference between the chief and deputy satellites during formation keeping ## Footnotes * **SUPPORTING INFORMATION** Additional supporting information may be found online in the Supporting Information section at the end of this article. * Received October 3, 2018. * Revision received December 30, 2019. * Accepted January 14, 2020. * © 2020 Institute of Navigation This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. ## REFERENCES 1. 1.Bandyopadhyay S, Foust R, Subramanian GP, Chung SJ, Hadaegh FY. Review of formation flying and constellation missions using nanosatellites. Journal of Spacecraft and Rockets. 2016;567-578. 2. 2.Escoubet CP, Schmidt R, Goldstein ML. Cluster-science and mission overview. In: The Cluster and Phoenix Missions. Dordrecht: Springer; 1997:11-32. 3. 3.Curtis S. The magnetospheric multiscale mission … resolving fundamental processes in space plasmas. Technical Report. 1999:1-56. 4. 4.Tapley BD, Bettadpur S, Ries JC, Thompson PF, Watkins M. GRACE measurements of mass variability in the Earth system. Science. 2004;305(5683):503-505. [Abstract/FREE Full Text](https://navi.ion.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzMDUvNTY4My81MDMiO3M6NDoiYXRvbSI7czoxNzoiL25hdmkvNjcvMS8zLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 5. 5.Friis-Christensen E, Lühr H, Knudsen D, Haagmans R. Swarm—an Earth observation mission investigating geospacer. Advances in Space Research. 2008;41(1):210-216. [CrossRef](https://navi.ion.org/lookup/external-ref?access_num=10.1016/j.asr.2006.10.008&link_type=DOI) [Web of Science](https://navi.ion.org/lookup/external-ref?access_num=000253590400027&link_type=ISI) 6. 6.Alfriend KT, Vadali SR, Gurfil P, How JP, Berger LS. Spacecraft Formation Flying: Dynamics, Control and Navigation. UK: Elsevier; 2010. 7. 7.Anthes RA, Bernhardt PA, Chen Y, et al. The COSMIC/-FORMOSAT-3 mission: early results. Bulletin of the American Meteorological Society. 2008;89(3):313-333. [CrossRef](https://navi.ion.org/lookup/external-ref?access_num=10.1175/BAMS-89-3-313&link_type=DOI) [Web of Science](https://navi.ion.org/lookup/external-ref?access_num=000254848800013&link_type=ISI) 8. 8.Kelley MC. The Earth’s Ionosphere: Plasma Physics and Electrodynamics. 96 Academic Press; 2009. 9. 9.Kintner PM, Psiaki M, Humphreys T, et al. A Beginner’s Guide to Space Weather and GPS. Lecture Notes. 2008. 10. 10.Leitner J. A hardware-in-the-loop testbed for spacecraft formation flying applications. IEEE Aerospace Conference Big Sky, MT; 2001;2:615. 11. 11.Busse FD, How JP, Simpson J. Demonstration of adaptive extended Kalman filter for low earth orbit formation estimation using CDGPS. NAVIGATION. 2003;50:79-93. 12. 12.Burns R, Naasz B, Gaylor D, Higinbotham J. An environment for hardware-in-the-loop formation navigation and control. AIAA/AAS Astrodynamics Specialist Conference and Exhibit. 2004; Providence, USA:4735. 13. 13.Mohiuddin S, Psiaki ML. Satellite relative navigation using carrier-phase Differential GPS with integer ambiguities. AIAA Guidance, Navigation, and Control Conference and Exhibit. 2005; San Francisco, CA:6054. 14. 14.Gill E, Naasz B, Ebinuma T. First results from a hardware-in-the-loop demonstration of closed-loop autonomous formation flying. 26th Annual AAS Guidance and Control Conference. 2003; Breckenridge, CO. 15. 15.Leung S, Montenbruck O. Real-time navigation of formation-flying spacecraft using global positioning system measurements. Journal of Guidance, Control, and Dynamics. 2005;28(2):226-235. [CrossRef](https://navi.ion.org/lookup/external-ref?access_num=10.2514/1.7474&link_type=DOI) 16. 16.Yamamoto T, D’Amico S. Hardware-in-the-loop demonstration of GPS-based autonomous formation flying. NAVITEC. 2008; Noordwijk, Netherlands. 17. 17.Marji Q. Precise Relative Navigation for Satellite Formation Flying Using GPS. Master Thesis. University of Calgary, Canada; 2008. 18. 18.Eyer JK. A Dynamics and Control Algorithm for Low Earth Orbit Precision Formation Flying Satellites. Ph.D. Dissertation. University of Toronto, Canada; 2009. 19. 19.Eddy D, Giralo V, D’Amico S. Development and verification of the Stanford GNSS navigation testbed for spacecraft formation-flying. Technical Note. 2017. 20. 20.Park JI, Park HE, Park SY, Choi KH. Hardware-in-the-loop simulations of GPS-based navigation and control for satellite formation flying. Advances in Space Research. 2010;46(11): 1451-1465. 21. 21.Vallado DA. Fundamentals of Astrodynamics and Applications. 4th ed. Springer Science & Business Media; 2001. 22. 22.Montenbruck O, Gill E. Satellite Orbits. Heidelberg, Germany: Springer-Verlag; 2001. 23. 23.Jones KG. Vehicle reference frames and vehicle motion handling in SimGEN (PC) software. Spirent Communication plc. Technical Report. DGP01062AAA. 2011. 24. 24.Yaakov BS, Li XR, Thiagalingam K. Estimation with Applications to Tracking and Navigation. New York; Johh Wiley and Sons; 2001. 25. 25.Cloutier JR. State-dependent Riccati equation techniques: an overview. IEEE American Control Conference. 1997;2:923-936. 26. 26.Gaposchkin EM, Coster AJ. GPS L1-L2 bias determination. Technical Report 971. MIT Lincoln Laboratory; 1993. 27. 27.Rideout W, Coster A. Automated GPS processing for global total electron content data. GPS Solutions. 2006;10(3):219-228. 28. 28.Peng Y. GNSS-based Spacecraft Formation Flying Simulation and Ionospheric Remote Sensing Applications. Master Thesis. Virginia Polytechnic Institute and State University; 2017. 29. 29.Rodrigues FS, De Paula ER, Abdu MA, et al. Equatorial spread F irregularity characteristics over Sao Luis, Brazil, using VHF radar and GPS scintillation techniques. Radio Science. 2004;39 (1):1-9. 30. 30.Richmond AD, Ridley EC, Roble RG. A thermosphere/ionosphere general circulation model with coupled electrodynamics. Geophysical Research Letters. 1992;19(6):601-604. 31. 31.Xu D, Morton Y, Jiao Y, Rino C. Simulation and tracking algorithm evaluation for scintillation signals on LEO satellites traveling inside the ionosphere. IEEE/ION PLANS 2018. Monterey, CA; April 2018:1143-1150. 32. 32.Forte B, Aquino M, Elmas ZG, Vadakke SV. On the use of the of the Cornell Scintillation Model within a Spirent signal simulator for the simulation of ionospheric scintillation events. International Ionospheric Effects Symposium (IES2011). 2011; Alexandria, VA. 33. 33.Rino CL, Tsunoda RT, Petriceks J, Livingston RC, Kelley MC, Baker KD. Simultaneous rocket-borne beacon and in situ measurements of equatorial spread F—intermediate wavelength results. Journal of Geophysical Research: Space Physics. 1981;86 (A4):2411-2420. [1]: /embed/graphic-3.gif [2]: /embed/inline-graphic-1.gif [3]: /embed/inline-graphic-2.gif [4]: /embed/graphic-4.gif [5]: /embed/inline-graphic-3.gif [6]: /embed/inline-graphic-4.gif [7]: /embed/inline-graphic-5.gif [8]: /embed/inline-graphic-6.gif [9]: /embed/inline-graphic-7.gif [10]: /embed/inline-graphic-8.gif [11]: /embed/inline-graphic-9.gif [12]: /embed/inline-graphic-10.gif [13]: /embed/graphic-6.gif [14]: /embed/graphic-7.gif [15]: /embed/inline-graphic-11.gif [16]: /embed/graphic-8.gif [17]: /embed/inline-graphic-12.gif [18]: /embed/inline-graphic-13.gif [19]: /embed/inline-graphic-14.gif [20]: /embed/inline-graphic-15.gif [21]: /embed/inline-graphic-16.gif [22]: /embed/inline-graphic-17.gif [23]: /embed/graphic-9.gif [24]: /embed/inline-graphic-18.gif [25]: /embed/inline-graphic-19.gif [26]: /embed/inline-graphic-20.gif [27]: /embed/inline-graphic-21.gif [28]: /embed/inline-graphic-22.gif [29]: /embed/inline-graphic-23.gif [30]: /embed/inline-graphic-24.gif [31]: /embed/inline-graphic-25.gif [32]: /embed/inline-graphic-26.gif [33]: /embed/inline-graphic-27.gif [34]: /embed/inline-graphic-28.gif [35]: /embed/graphic-10.gif [36]: /embed/inline-graphic-29.gif [37]: /embed/inline-graphic-30.gif [38]: /embed/inline-graphic-31.gif