Near Real Time Satellite Event Detection, Characterization, and Operational Assessment Via the Exploitation of Remote Photoacoustic Signatures

Current active satellite maneuver detection techniques can resolve maneuvers as quickly as fifteen minutes post maneuver for large Δv when using angles-only optical tracking. Medium to small magnitude burn detection times range from 6 to 24 h or more. Small magnitude burns may be indistinguishable from natural perturbative effects if passive techniques are employed. Utilizing a photoacoustic signature detection scheme can allow for near real time maneuver detection and spacecraft parameter estimation. We define the acquisition of hypertemporal photometric data as photoacoustic sensing because the data can be played back as an acoustic signal. Studying the operational frequency spectra, profile, and aural perception of an active satellite event such as a thruster ignition or any subsystem operation can provide unique signature identifiers that support resident space object characterization efforts. A thruster ignition induces vibrations in a satellite body which can modulate reflected sunlight. If the reflected photon flux is sampled at a sufficient rate, the change in light intensity due to the propulsive event can be detected. Sensing vibrational mode changes allows for a direct timestamp of thruster ignition and shut-off events and thus makes possible the near real time estimation of spacecraft Δv and maneuver type if coupled with active observations immediately post maneuver. This research also investigates the estimation of other impulse related spacecraft parameters such as mass, specific impulse, exhaust velocity, and mass flow rate using impulse-momentum and work-energy methods. Experimental results to date have not yet demonstrated an operator-correlated detection of a propulsive event; however, the application of photoacoustic sensing has exhibited characteristics unique to hypertemporal photometry that are discussed alongside potential improvements to increase the probability of active satellite event detection. Simulations herein suggest that large, potentially destructive modal displacements are required for optical sensor detection and thus more comprehensive vibration modeling and signal-to-noise ratio improvements should be explored.


Introduction
The space environment near Earth has grown crowded since the dawn of the space age in the late 1950s. Technological advances and decreases in manufacturing costs have led to a significant increase in the space object population. Increased launch frequency by government agencies and the private sector coupled with the growing number of uniquely controlled payloads that can be included per launch demonstrate the need for a robust system to monitor and protect each near-Earth orbital regime.
A challenge inherent to the crowded space environment is how to effectively predict Resident Space Object (RSO) collision risks and use that information to ensure active satellite survival and compliance with the Inter-Agency Space Debris Coordination Committee guidelines [1]. The 18th Space Control Squadron (18 SPCS) at the Combined Space Operations Center provides conjunction data messages via a global array of sensors to operators with assets that have secondary RSOs within their screening volume. Operators will use this information to maneuver their assets into a safer orbit if required. Some operators will send a predicted ephemeris including the planned maneuver to the 18 SPCS such that they can preemptively screen it against any other secondary RSOs. The screening is to ensure the maneuver does not put the satellite into the path of another object, effectively defeating the purpose of the collision avoidance maneuver in the first place. There are, however, instances where an operator does not follow this procedure or an unexpected dynamic event occurs, and thus maneuvers and potential conjunctions are unpredictable. This scenario demonstrates the need for a method to quickly detect when an uncooperative RSO has maneuvered and estimate its new trajectory. If there were such a method, it would not only help reduce potential conjunction risks much sooner, it would also provide data that would allow satellite operational capability assessments.
In addition to the need for near real time monitoring of active satellites, there are still gaps in how both active and inactive RSOs are characterized and uniquely identified. Any source of data that would provide unique signature identifiers based on RSO shape, behavior, mass, or other parameter increases our ability to better understand the RSO population. Thus, this research presents a new method to detect satellite maneuvers in near real time and extends the method to support RSO characterization and spacecraft parameter estimation.

Photoacoustic Sensing and Characterization
A noteworthy phenomenon becomes observable if photometric data is collected at a rate above 40 kHz. The human ear is capable of discerning frequencies from 20 Hz to 20 kHz. Thus, an acoustic signal can be generated from hypertemporal photometric data by converting the frequency content within a light curve to audio. To abide by the Nyquist Theorem, sampling at 40 kHz plus a safety margin would recover all naturally discernable sound. We define the acoustic playback of high-rate photometric data as its photoacoustic signature. The light to sound conversion is possible because photons carry an equivalent information content that an acoustic wave cannot because there is no medium for it to propagate through in the vacuum of space [2,3]. This technique has been proven robust in terrestrial tests that demonstrated the capability to accurately recover a conversation or song on the radio by collecting the light being modulated by a flexible surface nearby the acoustic source such as a plant leaf [4]. The direct analogy to applying photoacoustic sensing remotely is to imagine a satellite's solar array as the leaf in the prior example and a thruster ignition event as the acoustic source from the radio.
An example of how an acoustic interpretation of an inherently non-acoustic signal supported the characterization of a physical event is seen in the analysis of plasma wave data from the Voyager 1 spacecraft. NASA mission scientists converted the vibrations of dense ionized gas detected by the spacecraft's plasma wave instrument into an acoustic signal. While listening to the acoustic playback, the mission scientists noticed a rising tone which helped infer a continuously increasing density profile, detection of interstellar plasma, transit of the heliopause, and thus a departure from our Solar System [5]. The aural perception of similar RSO events may be one modality used to characterize satellites in a manner similar to how biometric recognition systems operate. Organizing RSO photometric event signal information into its unique frequency content, transients, pitch, aural perception, dead zones, harmonics, power level, profile, and other categories could provide the modalities needed to implement an RSO recognition system. These modalities would fulfill the universality, distinctiveness, permanence, and collectability requirements for such a system [6,7]. Fusing these modalities with other known or inferred spacecraft parameters could yield the equivalent of an n-factor authentication system that identifies and characterizes RSOs.

Resident Space Object Event Detection
The prior section discussed the potential for RSO event characterization and identification via acoustic interpretation of remote photoacoustic signatures obtained from active satellite observations. The satellite events assumed detection through either direct thruster plume observations for maneuvers, subtle changes detected in the hypertemporal photometric data frequency content due to an on-board event, or a combination of the two. This section will define specifics of the event recognition techniques from a satellite vibration mode change detection standpoint, focusing primarily on maneuvers. While the focus remains on propulsive events, these techniques could readily be applied to collision, explosive, or disintegration events, thermal snap, and abrupt changes in the space environment.
Methods to detect unmodeled dynamic events are often passive techniques that use an algorithm to sample historical ephemeris data and apply statistical methods until it can suggest an object's trajectory has deliberately shifted. Assuming an unmodeled dynamic event is a maneuver and depending on the magnitude of the burn, this sort of event detection can take anywhere from around ninety minutes up to several days to resolve. If the burn is small enough, it may be indistinguishable from the natural perturbative effects [8]. Large maneuvers may cause a complete loss of an object's trajectory and can be tedious to reacquire. An active detection technique for a Geosynchronous Earth Orbit (GEO) satellite using ground based optical tracking and sequential estimation tools showed that a Δv of 1.0 m/s could be detected as soon as fifteen minutes after a maneuver while a Δv of 0.01 m/s could take 12-24 h to discern with confidence [9]. Another promising event detection example is demonstrated in publicly available footage from ExoAnalytics which captures thruster plumes in the geosynchronous regime. 1

Maneuver Detection
Collecting hypertemporal photometric data on an active satellite has the potential to detect and characterize on-board operational events in near real time. It is assumed that the ignition and nominal operation of a thruster imparts energy into a satellite body. This energy may induce vibration in flexible components depending on the geometry, material, and inertial properties of the satellite. If any of these structures have reflective surfaces, the reflected photon flux will be modulated by the induced vibration. Further, if the relative displacement of the reflective surfaces is large enough, it should be possible to detect the change in outgoing photon flux with an optical telescope.

Synthetic Light Curve Simulation
To demonstrate an event detection in this manner, a 2 kHz synthetic light curve is simulated for an active satellite with a nadir aligned, velocity constrained attitude profile in GEO assuming the basic box-wing structural model in Fig. 1 and a simplified Cook-Torrance reflectance model [10]. The light curve model was modified to operate on a flat-plate surface and relies on both diffuse and spectral bidirectional reflectance distribution functions, incidence angles, and reflected intensity values. To calculate the apparent visual magnitude as measured by an observer, the following formula can be implemented: where m v, sun is the apparent visual magnitude of the sun, ρ d is the distance between the observer and RSO, N facets is the number of facets of the box-wing model, and the sum term denotes the effect from each given facet, further defined in Eq. (2). The last term in Eq. (1) represents the additive effect of atmospheric scintillation defined as a Gaussian distribution with mean μ and σ 2 v;atm representing its variance. Unless explicitly stated otherwise, a 0.04 m v standard deviation and zero mean are assumed for the noise distribution, also denoted as N(0,0.04 2 ) as an additive effect to an object's apparent visual magnitude. 2 The intensity relative to the sun's apparent brightness generated from each facet and reflected in the direction of the observer is where d and s are the fraction of the spectral and diffuse bidirectional reflectance values, R d, i and R s, i , respectively, A i is the facet surface area, N̂i is the surface unit normal vector, L̂is the sun to RSO unit vector, and V̂is the RSO to observer unit vector. Inspecting the spectral term from [10], where D is the facet slope distribution function, G i is the geometrical attenuation factor, and F is the surface reflectance of a perfectly smooth surface, demonstrates the dependance of reflectance on both material and geometric components. The material related elements that make up the F term, namely the diffuse and spectral reflectance coefficients for each facet, C d, i and C s, i , influence the apparent brightness and how an RSO's dynamics are affected by solar radiation pressure. Spacecraft reflectance and surface area parameters implemented in the simulation are listed in Table 1. The reflected light vector and relative surface angles are depicted in Fig. 2. The 2 kHz value is chosen to ensure at least a kilohertz order of magnitude is available for acoustic analysis after application of the Nyquist Theory, to adequately capture any event frequencies, and to minimize computational requirements.
Vibration Model To simulate a thruster fire event, a fixed-free cantilever beam mode for the solar panel and a face shear mode for the satellite bus are induced as depicted in Fig. 3 and Fig. 4. A simple sinusoidal oscillation model is used for the satellite bus and solar panel surface unit normal vectors at the desired frequency to achieve the effect of a vibrating surface. It is assumed no complex deformation exists in the solar panel and thus the axis of rotation for the rigid wing model is defined at the wing-box interface. A three-dimensional surface deformation and macro model simulation should be used for an improved representation of a satellite's photometric event fingerprint. It seems reasonable that the simplistic model implemented herein and reality exhibit behavior on the same order of magnitude due to the small surface displacement magnitudes involved, especially for the satellite bus represented by the box model. A range of burn durations less than or equal to five seconds are investigateda realistic value for collision avoidance maneuvers. The event frequency is chosen somewhat arbitrarily to be 58 Hz but is loosely based on anecdotal and proprietary industry evidence that modes below 50 Hz are deliberately damped on various satellites due to their instrument vibration isolation requirements. It appears that some active satellite vibration modes exist in the 3-12 Hz range from on-orbit test data [11]. Another study that used cameras mounted on the GOSAT satellite to analyze solar panel vibrations caused by thermal snap and gas jet thruster firing events showed that a 0.215 Hz out-of-plane and 0.459 Hz in-plane mode were observed [12]. Both sets of on-orbit test data generally agree with the expected first natural frequency of a solar array described by the Table 2 parameters and fixed-free cantilever beam calculation in Eq. (4),  where K 1 is the first mode parameter constant equal to 3.52 and ω is the uniform load per unit length. The solar panel is designed to mimic the specifications in [13,14] and the material is assumed to be an aluminum composite layer with silicon solar cells. If a rectangular shaped solar array with dimensions of 6.0 m × 2.5 m were simulated, the expected first natural frequency is calculated as 0.79 Hz. If the lower natural frequencies are not intentionally damped, the 58 Hz value is likely overzealous if attempting to emulate a detectable event as the largest displacements will be generated from the first mode. The expected natural frequencies seem to suggest that collecting hypertemporal data in the kilohertz range is overkill. However, the high-rate collection remains useful in increasing event epoch definition precision and thus parameter estimation accuracy discussed further in Section 4, is also useful in supporting the subsystem characterization efforts from an acoustic standpoint in Section 2, and helps separate the effects of various noise sources from event frequencies reviewed in Section 6. Again, it is postulated that the acoustic interpretation of a satellite event can augment traditional frequency analysis like the Voyager 1 example  to identify subsystem characteristics, e.g. a momentum wheel desaturation event would sound different than an imaging instrument scanning a target. While the 58 Hz event frequency is referenced for the following analysis, the conclusions remain the same if it were run at a lower magnitude nearer the first natural frequency. A range of surface displacement values are simulated at the event frequency such that a minimum detectable displacement threshold can be determined. All referenced surface displacement values are calculated as the peak or largest deflection at the extrema of the box face and solar panel surface areas respectively as defined in Eq. (5), Eq. (6), and Figs. 3-4 as where ϕ is the angle by which the surface unit normal vectors are rotated. The face shear mode surface unit vector rotation point for the box model is defined on its area centerlines, hence the factor of four in the denominator of Eq. (5). The box-wing interface acts as the rotation point for the fixed-free beam mode in Fig. 3. The simulation also assumes a step function for the thruster ignition and shut-off events with no thrust profile effects included.
Large Displacement Test Case To study the effects of a simulated maneuver on the apparent visual magnitude, an arbitrarily large displacement test case is run with and without the additive (0,0.04 2 ) noise from the atmosphere. Results of the proof of concept are shown in Fig. 5. The simulated maneuver begins at the 1.0 s since epoch tick mark. A peak surface displacement d body and d sp on the order of tens of centimeters is needed to visually detect a change in apparent visual magnitude in a time series plot. A displacement value this large is unrealistic and would be a destructive mode. Potentially realistic surface displacements of the single-to sub-centimeter range are lost to the atmospheric noise as in Fig. 6 and require more advanced methods to detect an event, placing importance on the signal-to-noise ratio achieved by the optical equipment utilized for data collection discussed further in Section 6. The +11.1 m v average apparent visual magnitude value seems reasonable for GEO distances assuming the reflectance properties defined in Table 1.

Inferring Signals from Noise Structure
To investigate if any realistic displacement magnitudes can be detected in the presence of atmospheric scintillation, three techniques were implemented. The first is a visual spectrogram review, second a Fast Fourier Transform (FFT), and third a custom crosscorrelation scanning algorithm. The data collection window lasted six seconds with a simulated dynamic event occurring from 2.0-4.0 s.  Fig. 7. While the underlying technique used in the pspectrum function is an FFT, the function itself required a substantial amount of tuning to properly identify a signal with a visual inspection. Thus, a more direct analysis is implemented in the following subsection.
Fourier Transform A slightly different formulation of the short-time Fourier Transform in MATLAB using the fft function that implements a single-sided amplitude spectrum analysis was performed on the same data to determine if analyzing the raw output would allow for increased signal resolution. Unlike the spectrogram, the direct FFT power output allowed for an obvious event detection at the simulated 58 Hz event frequency as shown in Fig. 8. Decreasing the displacement magnitude d body again from 1.73 cm demonstrated a loss of signal near 1.0 cm with this technique.

Cross-Correlation Scanning Algorithm
The base form of both the Fourier Transform and the cross-correlation technique is a convolution and thus they both Thus, a custom cross-correlation signal scanning algorithm using MATLAB's xcorr function was written to scan all possible frequencies, phase lag, and burn durations in an attempt to find the highest correlated reference signal [15]. The range of possible values for vibrational mode frequency parameters is plausibly constrainable to the 0-100 Hz range based on the natural frequency analysis in Section 3.1 which significantly reduced the search space. Thus, a proof of concept test case is run using the same large displacement as Fig. 7 showing successful event resolution in Fig. 9. Implementing the scanning algorithm on prior loss of signal cases slightly improved the d body detection criteria to 7.0 mm as shown in Fig. 10. A best-case atmospheric noise distribution of N(0,0.03 2 ) allowed detection to a 5.2 mm displacement for the satellite body and 1.16 cm for the solar panel with this technique.
The signal-to-noise ratio (SNR) for all simulations was calculated using Eq. (7) in units of decibels and yielded the values in Table 3. The numerator inside the log term of Eq. (7) is the variance of the event signal while the denominator is that of the atmospheric noise source, 0.04 2 in most cases. The results suggest a lower bound on the effective SNR near −27.55 dB for satellite event detection via cross-correlation.
While this simulation assumed the atmosphere was the dominant source of apparent visual magnitude uncertainty, other sources such as instrumentation error should be considered and carefully calibrated. This research also assumes maneuvers are observable at night and the start and end epochs of the propulsive event are directly observable with a high rate, sensitive photometer. It may be possible to extend this technique into the non-visible spectrum to detect events during the daytime. Additionally, the target RSO must be acquired and tracked within the telescope field of view (FOV) throughout the event duration. For satellites in LEO this could pose a challenge, especially for long, non-impulsive maneuvers. It seems plausible to employ these tactics to satellites in GEO based on their more stationary target orbits. The logical next inquiry that arises is what the limiting detection factor becomes if there were no atmospheric turbulence to corrupt an event signal. An immediate application for such a scenario with no or less atmosphere is a satellite-to-satellite observation at a geosynchronous orbital distance.  Other noise sources such as vibration induced by fuel slosh or other spacecraft subsystem operations may need to be overcome once the atmosphere error is eliminated.
Considering the absolute best-case atmospheric scintillation assumptions and achieved SNR from simulation, a detectable d body value of 5.20 mm may be approaching the realm of a realistic satellite bus displacement although it still has the potential to remain a destructive mode at this order of magnitude. With regard to the best-case solar panel displacement demonstrating a detection, considering the tip deflection of 1.16 cm at the end of a 3.872 m fixed-free cantilever beam seems like a plausible magnitude for an on-orbit mode. However, more research using a threedimensional macro vibration model with precise material and shape selections would be required to determine realistic on-orbit displacement values given an accurate forcing function from a propulsive event.

Maneuver Estimation
Using the dynamic event detection techniques of the prior section, photoacoustic sensing provides a means to directly timestamp maneuvers and other operational events in space. With event epoch error bounds as large as a full orbital period demonstrated in some of the prior literature review, it becomes impossible to accurately estimate spacecraft parameters that involve any sort of time dependency with this technique. However, direct event observations via hypertemporal photometry provide the ability to estimate impulse-related spacecraft parameters that were not possible with legacy techniques due to the large error incurred by generally comparing orbital changes before and after an estimated event epoch. While tracking the RSO throughout the entire event duration would support photoacoustic characterization efforts, only the initial and final event epochs require observability to allow for estimation of the following parameters. For propulsive events, the ability to directly define the impulse duration allows for estimation of Δv magnitude, direction, and maneuver type once the RSO state is resolved to a desired level of uncertainty post maneuver. We define "near real time" estimation as the time duration between maneuver end and post-maneuver state resolution. To study how sensitive the aforementioned parameters are to direct event start and end epoch observation precision assumed via photoacoustic sensing, a Δv simulation is run. Satellite operators at Optus, an Australian telecommunications company, were kind enough to provide both the states before and after stationkeeping maneuvers as well as exact maneuver start epochs and burn durations for various satellites in their constellation in GEO as a comparison to truth.
The dynamic model used for the propagation of the Optus satellite included a 20 × 20 EGM-96 gravity model, luni-solar perturbations, and an area-averaged exponential drag and solar radiation pressure model. It is assumed that the ballistic coefficient and equivalent solar coefficient-area-mass ratio were estimated via observations made on the RSO prior to observing a maneuver. Propagating the true Optus states provided before and after the maneuver to the given maneuver end epoch and computing their differences yielded the results in Table 4. Using a ±0.50 s uncertainty on the maneuver end epoch incurred the errors also listed in Table 4. Thus, for photoacoustic sensing to act as a reliable parameter estimation scheme, event start and end epoch observations must perform to a tenth-second or less uncertainty range to bound the estimation error to low single digit percentages. In the ExoAnalytics data referenced in Section 3, it is Orbit Determination The prior Δv simulation assumed operator provided states were available post-maneuveran unrealistic assumption when studying uncooperative or non-operational RSOs. To investigate how quickly a post-maneuver RSO state can be resolved via observation, an Extended Kalman Filter (EKF) is implemented using simulated range and range-rate measurements from three near-equatorial stations for an RSO in a 791 km near-circular orbit. A State Noise Compensation (SNC) algorithm was included in the EKF to prevent filter saturation and account for any dynamic model uncertainty. The ground station locations were chosen to mimic the Kwajalein Atoll, Diego Garcia, and Arecibo sites currently in the U.S. Space Surveillance Network, although the Arecibo site was not visible during the time span shown in Figs. 11-12. A complex Earth model using an FK5 precession, nutation, and polar motion correction Fig. 11 Post-fit state 3σ standard deviation trace output from the EKF starting from maneuver end epoch to solution convergence. The apparent discontinuity at 27 min is caused by obtaining measurements from the Kwajalein Atoll station acquired after a measurement gap from 10 to 27 minthe new measurement geometry reduces uncertainty and the effects of light time were accounted for in the measurement simulation. Measurement uncertainties based on demonstrated values by LeoLabs [16,17], observation frequency, and a priori covariance values are given in Table 5. To inject operational realism into the scenario, a seventeen-minute observation gap is included in the time study between the two visible ground stations. As seen in Figs. 11-12, the EKF required eleven observations at a sixty second interval over a span of thirty minutes to achieve a 13-m 3σ state uncertainty post maneuver. Back propagation to the maneuver end epoch with a Rausch-Tung-Striebel smoother maintained an uncertainty well within the referenced 100-m 3σ state uncertainty bound referenced previously.
With more consistent measurements, priority tasking, and ideal station coverage, the uncertainty convergence time can be improved. While the 3σ state uncertainty is a common metric used to confidently resolve an RSO state, it is not the only way to allow near real time estimation. For the same simulation, the measurement residuals and innovations covariance in Fig. 12 can give clues to solution stability. After only two measurements, the range residual has stabilized and the 3σ innovations covariance has achieved a steady 20-m value. Thus, it seems generally plausible to estimate Tables 4 and 6 parameters within five minutes of an active satellite maneuver using this technique and given assumptions. Studying residual behavior to indicate state stability should be used with caution as residuals alone do not guarantee an accurate orbit and should be later substantiated with the traditional state covariance convergence or other metrics.

Spacecraft Operational Assessment
Near real time Δv and maneuver type estimation defined in Section 4 allow for rapid response to uncooperative satellite operations and anomalous unmodeled dynamic events. Knowledge of these parameters combined with the now directly observed event epochs thanks to the photometric event detection techniques of Section 3 allow for a deeper understanding of spacecraft operational parameters specific to propulsive events. Through the application of orbital momentum and energy conservation equations about the known maneuver impulse, the prior techniques can be extended to estimate thruster mass flow rate, specific impulse, and fuel consumed given an a priori mass estimate. While the origins of the specific orbital momentum and energy integrals of motion in Eq. (8) and Eq. (9) are well established, it is important to understand their derivation and application as it has direct implications into the ability to estimate certain spacecraft mass properties. The base form of the specific orbital momentum ĥis defined as r̂Â v̂¼ ĥð8Þ where r̂is the satellite position vector, v̂is its velocity vector, and the specific orbital energy ε is defined by where v and r are the magnitudes of the state vectors from Eq. (8) and μ is Earth's gravitational constant. To infer operational information about an active satellite, the maneuver specific parameters must be related to known orbital quantities. In its most basic form, a maneuver is simply an impulse on this system and the initial and final orbital momentum and energy are calculable quantities from the research of Section 4. Thus, it is plausible to expand the specific orbital momentum and energy equations to their mechanical forms and include the varying work and impulse terms imparted by any dynamic events. The general form of the conservation of momentum equation is defined in Eq. (10) as where F is the force imparted on an object from time t i to t f , Δp is the total change in momentum, and m i , m f and v i , v f are the initial and final masses and scalar velocities of the object, respectively. The work W imparted on an object due to a force F is defined in Eq. (11) as where C is the trajectory from initial position x(t i ) to final position x(t f ), v is the velocity, and E i , E f are the initial and final energies. Extending these basic equations to the orbit problem and converting from scalar to vector quantities while including any rotational terms will yield the expanded forms in Eq. (12) and Eq. (13). These forms assume the only force contributing to the respective work and impulse integrals is due to the maneuver and that the object mass, position, and velocity are time varying. For the full orbital impulse-momentum conservation equation, where the parameter at h t ð Þ denotes the time-varying acceleration vector due to thrust, m : is the thruster mass flow rate, τ̂t ð Þ is the time-varying rotational torque imparted on the spacecraft by the maneuver, Δt is the duration of the impulse, I i and I f are the initial and final mass moments of inertia, ω̂i and ω̂f are the initial and final angular velocity, and finally h i and h f are the initial and final momentum wheel terms. Regarding the orbital work-energy conservation equation, where the U E terms describe the gravitational potential energy at a given position and L̂t ð Þ is the angular work performed over the maneuver. Momentum or energy losses due to imperfect thrust vectoring, momentum wheel saturation, or the change in moment of inertia due to decreasing spacecraft mass are assumed to be negligible for a satellite with a nadir aligned, velocity constrained attitude profile over a short maneuver duration. Thus, the rotational energy, work, torque, and momentum terms cancel over the time span of the maneuver and can be removed from the equations. The remaining form consists of two equations and five unknownsmass, mass flow rate, and the position, velocity, and acceleration due to thrust as functions of time during the impulse. To solve for the time varying functions, one can divide the Δv magnitude by the burn duration to get the average acceleration magnitude over the maneuver. The direction of the acceleration is the same as the Δv pointing vector. Thus, the force model can be augmented to include the acceleration due to the thrust and the satellite position and velocity as functions of time can be discretized. A fourth order polynomial is then used to fit each component of the discretized state over the maneuver, thus reducing the unknowns from five to two.
An analytical form of the acceleration due to thrust as a function of time can also be formulated based on the specific impulse in Eq. (14) and rocket equations defined in Eq. (15) as where F th is the overall force imparted on the satellite due to thrust, I sp is the specific impulse of the thruster, g 0 is the standard acceleration due to gravity at sea level, a th is the acceleration due to thrust, and v e is the exhaust velocity. Continuing with the rocket equation, where Δm is the change in mass and Δv is the change in velocity imparted on the satellite, or simply the delta-v. Solving for the scalar acceleration due to thrust, vectorizing, and injecting time dependency yields the form in Eq. (16). The acceleration unit vector as a function of time aû t ð Þ can either use the polynomial form mentioned previously or the constant Δv relative pointing vector derived in Section 4.
Substitution into the simplified mechanical momentum and energy equations yields their final forms in Eq. (17)(18)(19) below.
With two equations and mass and mass flow rate remaining as the final two unknown parameters, it seems as if it would be possible to solve the system of equations. However, the prior conversion from specific to mechanical energy and momentum in Eq. (12) and Eq. (13) sheds the linear independence with respect to the mass and mass flow rate terms. If one attempts to solve for the two unknown parameters, the resulting output will show no unique solution. One potential workaround for the linear dependence problem is to determine if there are any numerical mass and mass flow rate pairs that satisfy Eq. (17)(18)(19). Unfortunately, inputting a grid of mass parameter pairs in a Monte Carlo fashion produces another ambiguous solution space with a range of pairs that satisfy the system balance. Thus, application of momentum and energy to the orbital impulse problem does not allow for standalone estimation of mass and massflow rate and an a priori estimate of one or the other is necessary to solve either equation.
Utilizing the true spacecraft mass provided by Optus in the same simulation from Section 4 to solve for the mass flow rate included in Eq. (17) yields an estimate within 2 % of the true value. Use of Eq. (14) and Eq. (15) then allows for calculation of the thruster specific impulse and exhaust velocity with results tabulated in Table 6.
The results in this section reiterate the importance of efforts in the Space Situational Awareness (SSA) community to estimate RSO mass and the resulting operational assessments that can be gained through estimates of the Table 4 and Table 6 parameters. One mass estimation method that seems promising uses astrometric and photometric data fusion [18]. Utilizing this method to get a good a priori mass estimate would subvert the linear dependence problem discussed previously.
A potential impact of estimating fuel consumption over time is the ability to monitor if an uncooperative satellite is saving enough fuel for a proper post-mission disposal or graveyard orbit. An estimate of specific impulse also helps characterize the operational capability and mission of an active satellite. Assuming the Table 6 parameters are known in addition to the operational frequency and potentially a derived thrust profile, it is conceivable to compare these propulsion system values to known model specifications and thus constrain a specific class or even unique model to an observed event. Identifying a thruster model may help with other RSO identification parameters such as launch date or country of origin.

Observations
The experimental goal of this research was to observe an active satellite thruster fire event and apply the estimation techniques of prior sections. Satellite operators at Optus and Iridium provided their maneuver schedule such that a precise orbital location and epoch were known for each station-keeping maneuver. The optical telescopes used are courtesy of the SERC 0.7 and 1.8m geotrackers at the Mt. Stromlo Observatory in Australia with collaboration from EOS Space Systems. The hypertemporal sampling rate detector is based on a Hamamatsu PMT sensor (H11901-20) that is sensitive over the entire visible spectrum. The photometric data is sampled at a rate of 50 kHz (up to 100 kHz), time stamped in UTC via GPS signal, and stored in binary files. The detector is built on the Beaglebone Black PC board with the real time operation software written in C++. The H11901 type PMTs operate with a highly sensitive photocathode that makes them suitable for even the most demanding photon counting applications. During the GEO observations at the 1.8 m telescope, the PMT detector was operated at the maximum gain of 4.0 × 10 6 and our estimations suggested that the system should respond to the low intensity of the photon flux reflected off the GEO targets.
The high sampling rate of 50 kHz is useful in that it helps separate out noise sources from any event signals. As predicted in Section 3.1, any detected events will likely manifest in the sub-100 Hz range. The effect of atmospheric scintillation is unstable and changes from 10 Hz to 100 Hz which is near the predicted event spectra. The 50 kHz sampling rate should be able to easily separate any dynamic event signatures from atmospheric effects or other noise sources. A potential limiting factor in the collection efforts is likely to be the SNR achieved by each telescope. The sensitivity of the 0.7-m telescope was defined to approximately a +9 stellar magnitude with an FOV of 0.5 deg. No sensitivity test was available for the 1.8 m due to required mirror maintenance.
Thus far, observations on these constellations have not yielded a maneuver detection yet. However, at least one anomalous event is included for further analysis. In addition to collaboration with Optus and Iridium, over twenty hours of data on 29 unique active satellites including the GLONASS constellation, Ajisai, Cosmos 2527, and Sentinel 1a & 1b were collected without any knowledge of their operational schedules. While some of the phenomena included in this section do not require hypertemporal photometric data nor photoacoustic conversion to be observed, the data are still included as a capability demonstration. Secondary measurements from LeoLabs also helped determine whether unexplained signals could be correlated to any unmodeled dynamic events.
Iridium flare Data collected on Iridium 14 revealed a classic flare, or abrupt peak in visual magnitude, due to favorable observation geometry and the relative attitude progression throughout its orbit. The flare included two auxiliary peaks on either side of the main peak shown in Fig. 13. It is likely the auxiliary peaks were caused by door-sized antennas offset by forty degrees from the panel antenna that is known to cause the main flare in light intensity. The photoacoustic playback of the flare allows for detection of the main intensity rise and fall as well as the auxiliary peaks. As this is a common satellite event that is easily detected in the photometric domain, the audio playback does not reveal any unique insights besides it being a demonstration of the acoustic domain conversion.
Ajisai tumble chirp The Experimental Geodetic Payload (EGP) satellite Ajisai was a useful calibration object for the detection system. Ajisai has a known spin rate of 1.25 Hz which can be seen in Fig. 14 plotted with the FFTW tool [19]. The photoacoustic playback of Ajisai's signal can be described as a partially muted tumble chirp at the aforementioned rate, or three times the base  frequency. Like the flare, deriving spin rates from photometry is a relatively well-known phenomenon and thus the acoustic analysis acts more as a confirmation of known RSO behavior.
SLR ticking In many of the collected data arcs, a 10 Hz and sometimes a 60 Hz signal would show up briefly and then disappear with no apparent correlation to satellite events. The signal itself sounded like a reverberating tick that was easily identifiable whenever it appeared. Eventually, it was discovered that the nearby satellite laser ranging (SLR) equipment operation correlated with the signal epochs. While the SLR operated at 60 Hz, the photometric system became blind during the laser pulses approximately every 100 ms due to what is likely a time domain aliasing effect. Thus, the 10 Hz signal can be acquired in an FFT. One example of the 60 Hz signal detection is shown in Fig. 15.
Australian power grid Another mystery signal is explained by the Australian power grid operating at 50 Hz. A commonly observed phenomenon is the flickering of exterior lights in close proximity to the observatory at two times the grid frequency as shown in Fig. 16. The lights are motion sensor activated, likely initiated by some nearby kangaroos, and thus can show up sporadically in the data. This particular signal was indistinguishable from background noise in the acoustic domain.
Uncorrelated burst-like transient A yet unexplained signal that appeared similarly at least twice across two unique satellite collections can be acoustically described as a short burst-like transient. Each time the profile appeared in the time-series data there were two distinguishable peaks. There were many other cases where a star passed the FOV which caused a gradual rise and fall in light intensity. This transient appeared much too short to be explained by a star passing. Acoustic playback of the potential event sounds somewhat noisy but does have a detectable dual-burst or knocking characteristic to it. The frequency content is rich in the sub-10 Hz range as shown in Fig. 17. It was confirmed  by Iridium that no maneuvers were planned while these transients were observed. Further analysis is needed to determine if the events originated from the spacecraft or if a measurement system aberration may have been the source.
Many factors may have contributed to the experimental null finding. The propulsion system on the observed Optus satellites may have simply not induced enough vibration in the satellite body or any modes that were activated were damped. The effective SNR from the optical system also may have been insufficient to distinguish any signals from atmospheric or other noise sources. The 1.8 m telescope was unavailable for most of the collection efforts and the 0.7 m geotracker was only sensitive to +9 m v and thus the 1.8 m was likely required to detect dynamic events near the predicted +11.1 m v . Also, the atmospheric condition, high clouds, water vapor, or dust could have attenuated any event signal and thus the overall intensity of light focused on the sensor may not have exceeded the noise level.

Conclusions
The research herein proposed a new method to support characterization of resident space objects by exploiting remote photoacoustic signatures. It was postulated that interpreting hypertemporal photometric data as a photoacoustic signal can support characterization of active satellites and distinguish between spacecraft subsystem events. With a maneuver start and end epoch actively observed via hypertemporal photometric data collection, it was shown that the burn duration, event frequency, Δv magnitude, direction, and maneuver type could be accurately estimated within five minutes of a completed maneuver. This near real time estimation scheme supports operational risk reduction by significantly shrinking the maneuver detection lag. It was also shown that with an a priori estimate of spacecraft mass, the thruster mass flow rate, specific impulse, exhaust velocity, and fuel consumed by a maneuver could be accurately estimated via impulse-momentum and work-energy methods. These parameters help characterize the operational capability of active satellites and can also support RSO identification. The ability to estimate fuel consumption supports efforts to monitor uncooperative satellite behavior for adherence to proper post-mission disposal actions. This research provides further motivation to experimentally detect an active satellite maneuver event via direct plume observation or vibration analysis. Current modeling suggests that relatively large, potentially destructive modal displacements are required for optical sensor detection and thus further refinements in vibration modeling and SNR improvements should be explored. The resulting knowledge that can be attained from such an event detection shows promise in characterizing RSO behavior and maintaining a sustainable space environment.