A Rapid Method For Orbital Coverage Statistics With \mathbf{J_2} Using Ergodic Theory
AAAS 21-222A RAPID METHOD FOR ORBITAL COVERAGE STATISTICS WITH J USING ERGODIC THEORY
Andrew J. Graven * , Alan H. Barr † , and Martin W. Lo ‡ Quantifying long-term statistical properties of satellite trajectories typically entailstime-consuming trajectory propagation. We present a fast, ergodic method of ana-lytically estimating these for J − perturbed elliptical orbits, broadly agreeing withtrajectory propagation-derived results. We extend the approach in Graven and Lo(2019) to estimate: (1) Satellite-ground station coverage with limited satellitefield of view and ground station elevation angle with numerically optimized for-mulae, and (2) long-term averages of general functions of satellite position. Thismethod is fast enough to facilitate real-time, interactive tools for satellite constel-lation and network design, with an approximate × GPU speedup.
Figure 1(a)
Space tracks of a satellite orbit (red at apogeeand yellow at perigee) intersecting the visibility cones(blue) of two ground stations (white). Computing the por-tion of time the satellite spends in view of its ground sta-tion via trajectory propagation is relatively slow. How-ever, using ergodic theory, we’re able to rapidly estimatethis quantity by integrating over the interior of the conewith respect to a carefully chosen probability measure.See
Figure 4 for a definition of the cone.
Figure 1(b)
An example from a key class of applications:real time data generation & visualization. This is a 3-dimensional heat map of percent visibility, from ρ , definedin Equation 1 and computed using Equation 14, for , distinct satellite orbits, generated in 4.9 seconds on a dualcore 2.7GHz laptop CPU, or 2 milliseconds on the TitanV GPU. The plot represents , orbits generated froma × × grid of semimajor axis, eccentricity andinclination values. * Department of Mathematics, Cornell University, Ithaca, NY, 14853, United States † Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA, 91125, UnitedStates ‡ Principal Engineer, Mission Design and Navigation Section, Jet Propulsion Laboratory, California Institute of Technol-ogy, Pasadena, CA, 91109, United States a r X i v : . [ a s t r o - ph . E P ] F e b NTRODUCTION
Dynamical Systems theory can be divided into two areas; the first consists of methods of solvingdifferential equations and is well known. The second area is ergodic theory, which is concerned withthe coverage and transport properties of dynamical systems and is less well known, perhaps due toits difficulty. Poincar´e made fundamental contributions to both areas, creating the geometric theoryof differential equations and pioneering the study of deterministic chaos in dynamical systems.Examples of problems and results in the domain of ergodic theory are the Poincar´e recurrencetheorem, the ergodic hypothesis of statistical mechanics, and the coverage of ground stations bysatellites in orbit about a central body. This last example is the primary focus of this paper.Advanced satellite and satellite constellation planning often requires an understanding of thelong-term behavior of the proposed orbit or constellation. Quantities such as average satellite-ground station visibility, atmospheric drag and sun β − angle may inform the likelihood that missionrequirements are satisfied. Estimating these by trajectory propagation can be costly and time con-suming due to the time scales and small step sizes necessary for accurate estimates. This is com-pounded by the high-dimensionality of the design space, adding a combinatorial challenge findingacceptable (or, what’s more, optimal) designs. Thus in many cases, such an approach may require asparse sampling of the design space, the use of low-fidelity simulations, or significant computationalresources.In this paper we present a fast analytical approach to estimating a wide range of long-term statis-tics of aperiodic J -perturbed circular and elliptical orbits. We provide formulae for estimating: (1)Satellite-ground station view period ratios, with limited satellite field of view (FOV) and groundstation elevation angle taken into account. (2) Averages of general functions of satellite position:drag force, gravity gradient, radiation, etc. Applying the Birkhoff-Kinchin Theorem of Ergodic the-ory, we express these quantities in terms of a definite integral. In certain cases, symmetries may beexploited to reduce the dimension of the integral, thereby further accelerating numerical evaluation.The evaluation of these formulae is sufficiently fast that it’s feasible to use these in real-time con-current engineering applications. These formulae turn out to be independent of the value of J aslong as J (cid:54) = 0 . Thus the results presented in this paper can be applied to any body with J (cid:54) = 0 .The value of J is only necessary to verify the aperiodicity of the orbit to guarantee its ergodicity.The mathematical symbols used throughout the paper are explained in-line and collected in the“Notation” section at the end of the paper. For the sake of compactness of notation, we define asymmetric truncation function, Trunc b ( x ) sending x to the nearest point in [ − b, b ] :Trunc b ( x ) := min { max { x, − b } , b } Throughout the paper it will also be convenient to extend the domain of cos − and sin − : sin − ( x ) := sin − ( Trunc ( x ))cos − ( x ) := cos − ( Trunc ( x )) The View Period Ratio ρ The view period ratio, ρ (also referred to as the “view period percentage”) for a given satellite-ground station pair is the asymptotic fraction of time the satellite is able to communicate with itsground station. Specifically, if T is the total flight time, and P ( T ) is the portion of flight time forwhich the satellite and the ground station are able to communicate, then the view period ratio, ρ , is2runc ( x ) sin − ( x ) cos − ( x ) given by the limit: ρ = lim T →∞ P ( T ) T (1)The existence of this limit implies that the approximation ρT ≈ P ( T ) improves * as T grows.And, under very mild assumptions, one can show that this limit indeed exists. Rapid computationof this quantity in various contexts is a central focus of this paper. The J Model
As a consequence of the Earth’s rotation, it isn’t a perfect sphere, but an oblate spheroid. Thisbreaks the assumption of radial symmetry of the central body in the standard 2 Body Problem,resulting in a perturbed potential and equations of motion. This perturbation can be quantified witha spherical harmonic model of the central body, which allows us to describe the perturbed bodyby a sequence of coefficients J , J , ... , the normalized zonal harmonic gravitational coefficients.For many applications, the J term effectively dominates the others. The first two terms for Earthare: J = 1 . · − and J = − . · − , with magnitude continuing to drop off for higherorder harmonics. Thus, most of this perturbation can be captured by the effect of J alone. The J − perturbed equations of motion are given in Equation 2. ¨ x = − µxr (1 − J R B r ( 5 z r − y = − µyr (1 − J R B r ( 5 z r − z = − µzr (1 − J R B r ( 5 z r − r = ( x + y + z ) (2) The Mean Linear J Model
Vallado shows that the average secular motion of this system is well-approximated by a linearprecession of the orbital elements: M = mean anomaly, Ω =
Longitude of the Ascending Node, ω = Argument of Periapsis. In particular, much of the long-term behavior of the system is accuratelycaptured by the linear flow: M ( t ) = M + ˙ M t Ω( t ) = Ω + ˙Ω t , mod πω ( t ) = ω + ˙ ωt , mod π (3) * In terms of relative error, not necessarily absolute error. ˙ M = (cid:114) µa [1 + 3 J R B a (1 − e ) (3 cos ( i ) − − (cid:114) µa R B a (1 − e ) J cos( i ) + Ω B ˙ ω = (cid:114) µa R B a (1 − e ) J (4 − ( i )) (4)This model is used throughout the paper. The Invariant Measure µ A key result in Ergodic theory, the Birkhoff–Khinchin Theorem, asserts a time mean-space meanequivelence for a certain class of “Ergodic” dynamical systems, see Arnold 1989 and Sinai 1976. If φ t ( x ) : S × R → S is the trajectory of an an ergodic system starting at x ∈ S , then the Birkhoff-Kinchin Theorem asserts that, for almost every x ∈ S , there exists a probability measure µ suchthat for any measurable function, f : S → R : lim T →∞ T T (cid:90) t =0 f ( φ t ( x )) dt = (cid:90) S f ( x ) dµ (5)The probability measure, µ , (also referred to an an invariant * measure) can be interpreted as theinfinitesimal proportion of time the state of the system spends at any given point in its state space.For example, if µ were the uniform distribution: µ = Vol ( S ) − , we could conclude that in thelong-term, (cid:126)x ( t ) spends the same amount of time in each region of its state-space.Here, the dynamical system of interest is the J − perturbed 2 Body Problem, with φ t ( x ) theposition of the satellite over time. We can’t expect µ to be uniform in this case, however. This isclear from the fact that the ground-tracks of satellite orbits are biased to extreme latitudes, as in Figure 2(a) . In addition, for elliptical orbits one should expect a bias of the distribution towardslarger radii due to the inverse relationship between satellite velocity and radial position
Figure 3(a) .Equation 6 from Graven and Lo 2019 provides the invariant measure for elliptical orbits. Thisextends the invariant measure for circular orbits from Lo 1994, provided in Equation 7. The Invariant Measure for Elliptical Satellite Orbits: µ e ( r, λ, L ) = r cos( λ )2 π a (cid:112) sin ( i ) − sin ( λ ) (cid:112) a e − ( a − r ) (6) The Invariant Measure for Circular Satellite Orbits: µ c ( λ, L ) = cos( λ )2 π (cid:112) sin ( i ) − sin ( λ ) (7) SATELLITE-GROUND STATION VIEW PERIOD RATIOS
The original motivation for this work was accelerating the computation of the view period ratio, ρ . Lo 1994 provided this result for circular orbits. In the following sections, we demonstrate theutility of the result for elliptical orbits. * invariant in the sense that if we define the flow of the system: F : S × R + → S by F ( x, t ) = (cid:126)x ( t ) s.t. (cid:126)x (0) = x ,then ∀ U ⊆ S, ∀ t > , µ ( F ( U, t )) = µ ( U ) * Asymptotes on the boundary of the distribution remove most of the detail from the standard heat map, so we normal-ize the distribution to lie in [0 , igure 2(a) Ground tracks of a circular orbit. The tracksexhibit a higher density at extreme latitudes, which isin agreement the well-known fact that ground stations athigher latitudes tend to have a higher level of connectivitywith their satellites. Note that the long-term density of theground tracks is longitude-independent.
Figure 2(b)
Normalized ∗ probability density µ c of circularsatellite orbit position, from Equation 7. µ c is a measureof the infinitesimal fraction of time the satellite spends atany given point in its state-space. Darker blue indicates ahigher probability. The probability density closely corre-sponds with the distribution of ground tracks in Figure 2(a) . Figure 3(a)
The space tracks of an elliptical orbit. Thetracks vary from yellow at perigee, to red at apogee. Whileit’s more difficult to see the distribution of space tracks herethan it is the for the ground tracks in
Figure 2(a) , theyare very similar. We observe higher space track densityat extreme latitudes just like before, with the main distinc-tion being that we also observe a greater relative density ofspace tracks at perigee and apogee than in-between.
Figure 3(b)
The normalized ergodic probability density µ e , from Equation 6, of an elliptical orbit. Analogousto Figure 2(b) , darker blue indicates a higher probabilitydensity. The analogy of the probability density in
Fig-ure 2(b) to the density of space tracks in
Figure 3(a) isless apparent on account of the higher-dimensionality andthe variation in velocity with respect to altitude distortingthe apparent distribution of space tracks. he General Problem Geometry Figure 4(a)
The elevation angle constraint on the regionof visibility. The elevation angle (cid:15) is the minimum com-munication angle above the horizon. The ground stationmask angle, θ , is the maximum communication angle be-tween it and the satellite. Figure 4(b)
Visualization of the satellite FOV constrainton the region of visibility from Equation 8. The FOV angle β is the largest angle off-nadir at which the satellite cancommunicate. A formula is provided in Equation 8. Let S be a satellite at some radius r , and suppose that its ground station has its line of sightconstrained by an elevation angle (cid:15) , such that it can only communicate with satellites at least (cid:15) radians above the horizon. Also suppose S has a nadir-pointing communication instrument withan FOV angle constraint, β , such that S can only send/receive within a cone of angular radius β about nadir. Then, the communication geometry is as in Figure 4. Applying standard trigonometricidentities, we can describe each of these constraints by the parameter θ , shown in Figure 4. θ , elev ( r ) = cos − (cid:18) R B r cos( (cid:15) ) (cid:19) − (cid:15)θ , FOV ( r ) = sin − (cid:18) rR B sin( β ) (cid:19) − β, sin( β ) < R B r cos − (cid:18) R B r (cid:19) , sin( β ) ≥ R B r (8)And if both constraints are in effect, one can compute the effective ground station mask angle bytaking the minimum of the two quantities from Equation 8. θ ( r ) = min { θ , elev ( r ) , θ , FOV ( r ) } (9)Finally, given a ground station mask angular radius θ , the ground station mask itself can be de-6cribed in coordinates suitable for integration, as in Equation 10. ( Ground Station Mask )( r ) = (cid:8) ( r, λ, L ) (cid:12)(cid:12) λ ∈ [ λ min ( r ) , λ max ( r )] L ∈ [ − L bound ( λ, r ) , L bound ( λ, r )] (cid:9) λ min ( r ) = Trunc i ( g − θ ( r )) λ max ( r ) = Trunc i ( g + θ ( r )) L bound ( λ, r ) = cos − (cid:18) cos( θ ( r )) − sin( λ ) sin( g )cos( λ ) cos( g ) (cid:19) (10)Of course, the proportion of time the satellite spends in the ground station mask will be exactly thesatellite-ground station view period ratio. Moreover, the invariant measure µ provides the infinites-imal proportion of time the satellite spends in each region of its state-space. Thus, the view periodratio can be computed as the integral over the ground station mask with respect to µ . We will bereferring to these formulae extensively in the following sections. One Satellite, One Ground Station
First, we will consider the most basic case of a single satellite in a circular orbit communicat-ing with a single ground station. This can then straightforwardly be extended to more complexgeometries.
The Circular Orbit View Period Formula
Suppose we have a satellite in a circular orbit withsemimajor axis a and inclination i . Then, recalling equations 7 and 10, we can compute the satellite-ground station view period ratio, ρ as: ρ = (cid:90) GroundStation Mask dµ c = λ max (cid:90) λ min L bound ( λ ) (cid:90) − L bound ( λ ) µ c ( λ, L ) dLdλ Plugging in the formulae for λ min , λ max , L bound and µ c , and evaluating the inner integral yields: ρ = 1 π Trunc i ( g + θ ) (cid:90) Trunc i ( g − θ ) cos( λ ) cos − ( cos( θ ) − sin( λ ) sin( g )cos( λ ) cos( g ) ) (cid:112) sin ( i ) − sin ( λ ) dλ Now, note that as the latitude approaches ± i , the integrand blows up, and that singularities in theintegrand tend to slow down numerical integration. To avoid this issue, we take the change ofvariables: λ ( α ) = sin − (sin( i ) sin( α )) (11)This yields a formula for the satellite-ground station view period ratio for circular orbits in terms ofa single definite integral, given in Equation 12. Accurate evaluation of this formula can be achievedwith only 10 Gaussian quadrature nodes, and enables ≈ , ratios/sec on a dual core 2.7GHzlaptop CPU, or ≈ . · ratios/sec on the core, . GHz
Titan V GPU.7 = 1 π α (cid:90) α cos − (cid:32) cos( θ ) − sin( i ) sin( α ) sin( g ) (cid:112) − sin ( i ) sin ( α ) cos( g ) (cid:33) dα,α = sin − (cid:0) sin ( Trunc i ( g − θ )) sin( i ) − (cid:1) ,α = sin − (cid:0) sin ( Trunc i ( g + θ )) sin( i ) − (cid:1) ,θ = min cos − (cid:18) R B a cos( (cid:15) ) (cid:19) − (cid:15), sin − (cid:18) aR B sin( β ) (cid:19) − β, sin( β ) ≤ R B a cos − (cid:18) R B a (cid:19) , sin( β ) ≥ R B a (12) Equation 12:
The simplified satellite-ground station view period ratio, ρ formula for circular orbits. a, i, (cid:15), β are the semimajor axis (km), orbital inclination (rad), ground station elevation angle (rad)and satellite FOV angle (rad) respectively. The Elliptical Orbit View Period Formula
Now, suppose we have a satellite in an elliptical orbitwith semimajor axis a , eccentricity e and inclination i . Then, noticing that for each radius r ∈ [ a (1 − e ) , a (1 + e )] , the cross-section of the cone of visibility is exactly described by Equation 10.Thus, we can compute the view period ratio as the integral over these masks with respect to µ e : ρ = a (1+ e ) (cid:90) a (1 − e ) λ ( r ) (cid:90) λ ( r ) L ( r,λ ) (cid:90) − L ( r,λ ) µ e ( r, λ, L ) dLdλ dr Applying the same approach as in the circular case, but extending the change of variables: r ( θ ) = a (1 − e sin( θ )) , λ ( α ) = sin − (sin( i ) sin( α )) (13)The simplified view period ratio formula for the elliptical case is given in Equation 14. ρ = 1 π π (cid:90) − π α ( θ ) (cid:90) α ( θ ) (1 − e sin( θ )) cos − (cid:32) cos( θ ( θ )) − sin( α ) sin( i ) sin( g ) (cid:112) − sin ( α ) sin ( i ) cos( g ) (cid:33) dαdθ,α ( θ ) = sin − (cid:0) sin ( Trunc i ( g − θ ( θ ))) sin( i ) − (cid:1) ,α ( θ ) = sin − (cid:0) sin ( Trunc i ( g + θ ( θ ))) sin( i ) − (cid:1) ,θ ( θ ) = min cos − ( P ( θ ) cos( (cid:15) )) − (cid:15), sin − (cid:18) sin( β ) P ( θ ) (cid:19) − β, sin( β ) ≤ P ( θ )cos − ( P ( θ )) , sin( β ) ≥ P ( θ ) P ( θ ) = R B a (1 − e sin( θ )) (14) Equation 14:
The simplified satellite-ground station view period ratio formula for elliptical orbits. a, e, i, (cid:15), β are the semimajor axis (km), eccentricity (unitless), orbital inclination (rad), groundstation elevation angle (rad) and satellite FOV angle (rad) respectively.8 ne Satellite, Many Ground Stations
The introduction of additional, potentially overlapping, ground stations turns out to non-additivelyincrease the complexity of computing the view period ratio. We conceptualize this problem as fol-lows: Given a set of N ground stations and a satellite trajectory, we can represent each groundstation by a latitude, longitude, ground station mask angle triple: ( λ , L , θ ) , ( λ , L , θ ) , . . . , ( λ N , L N , θ N ) and the satellite trajectory by its orbital elements ( a, e, i ) . Using this information we’d like todetermine the expected total visibility time for the satellite with the complete set of ground stations.We will begin with the circular case. The Circular Case
Here, we are given the information outlined above, except the satellite tra-jectory is circular, so we only need its semimajor axis a and inclination i . For example, the situationfor N = 5 may appear as in Figure 5. Figure 5 : Example set of ground station masks for multiple ground stationsNote that we can’t directly compute the view period ratio for each ground station and sum thembecause the ground station masks may overlap, resulting in the double counting of some regions.With a little work, it is possible to take this into account. The view period ratio formula in this caseis given by Equation 15. ρ = 12 π α max (cid:90) α min Length ( I ( λ ( α ))) dαα min = sin − sin Trunc i (cid:18) min k { λ k − θ k } (cid:19) sin( i ) α max = sin − sin Trunc i (cid:18) min k { λ k + θ k } (cid:19) sin( i ) (15)Dealing with the complexities introduced by the overlapping ground station masks is a lengthy pro-cess. Thus, the explanation of Equation 15, including the definition of definition of Length ( I ( λ ( α ))) ( I ( λ ( α ))) is at worst O ( N log( N )) if N is the number of ground stations. Thus, we only pick up a logarithmic term in the computationalcomplexity, in comparison to the case of non-overlapping ground stations. The Elliptical Case
Taking a similar approach, the formula in the elliptical case is given by: ρ = 12 π π (cid:90) − π α max ( r ) (cid:90) α min ( r ) (1 − e sin( θ )) Length ( I ( λ ( α ) , r ( θ ))) dαdθα min ( r ) = sin − sin Trunc i (cid:18) min k { λ k − θ k ( r ) } (cid:19) sin( i ) α max ( r ) = sin − sin Trunc i (cid:18) min k { λ k + θ k ( r ) } (cid:19) sin( i ) (16) Many Satellites, One (or Many) Ground Station(s)
The natural next step after considering the single satellite-single ground station and single satellite-multiple ground station cases is considering what can be said of the coverage if we’re working withan ensemble satellites and one (or several) ground station(s). If the satellite orbits are statisticallydependent, then this appears to be a difficult problem. However under the assumption of indepen-dence, a lot can be said using only the view period formulae from the previous sections.
Two Satellites, One Ground Station
Suppose we have two satellites, S and S , and a groundstation, G . Also suppose that the view period ratios for S and S are ρ and ρ respectively. Then,note that we can view ρ , ρ ∈ [0 , as the instantaneous probability that either of the satellites willbe visible to the ground station at any point in time. Then, because their trajectories are independent: P ( S and S visible to G ) = P ( S visible to G ) P ( S visible to G ) = ρ ρ (17)Similarly, the probability that exactly one of the satellites is in view of the ground station are ρ (1 − ρ ) and ρ (1 − ρ ) respectively. Moreover, we can compute the total coverage ratio: P ( S or S visible to G ) = ρ + ρ − ρ ρ (18)These extended coverage ratios have a very natural interpretation: They represent the expectedproportion of flight time for which the given condition holds. (e.g. ρ ρ is the expected proportionof flight time for which we can expect both satellites to have line-of-sight with G . Similarly, ρ + ρ − ρ ρ is the expected proportion of flight time for which we can expect at least one satellite tohave line-of-sight with G .) N Satellites, One Ground Station
From here, it’s straightforward to generalize to N satellites S , S , . . . , S N in communication with a single ground station, G , each satellite with a view periodratio, ρ i . Firstly, we can extend Equation 17 to compute the probability that any subset S { i ,...,i N (cid:48) } = { S i , S i , . . . , S i N (cid:48) } of the satellites is visible to G : P (cid:0) S { i ,...,i N (cid:48) } visible to G (cid:1) = P (cid:92) i ∈{ i ,...,i N (cid:48) } S i visible to G = (cid:89) i ∈{ i ,...,i N (cid:48) } ρ i (19)10ext, using P ( A ) = 1 − P ( ¬ A ) , we can see the overall coverage ratio of the N satellites: P ( ≥ of the S i visible to G ) = 1 − N (cid:89) i =1 (1 − ρ i ) (20)Clearly these formulae can be applied analogously when there are multiple ground stations simplyby replacing the single ground station view period ratios with the multiple ground station view pe-riod ratios. In addition, in the multiple ground station case, these formulae can be applied to analyzeindividual ground stations and, in general, any subset of the ground stations under consideration. THE LONG-TERM MEAN VALUE OF FUNCTIONS OF SATELLITE POSITION
Recall that the Birkhoff-Kinchin Theorem asserts a time mean-space mean equivalence for er-godic dynamical systems. This relation is given explicitly in Equation 5. In the case of computingview period ratios, the formulae provided in the previous sections can be thought of as applicationsof the theorem to f = Vis (the indicator function on the region of visibility of the ground sta-tion(s)), and (cid:126)x ( t ) the position of the satellite as a function of time. But this formula actually appliesin much greater generality. In particular, f : S → R n can be any measurable scalar or vector-valued(and even possibly tensor-valued) function on the state space of the satellite. For example:• Satellite speed: v : S → R , ie. v ( r, λ, L ) = (cid:115) µ (cid:18) r − a (cid:19) .• Atmospheric density: ρ : S → R .• Drag force per unit area: D : S → R , ie. D ( r, λ, L ) = Cd ρ ( r, λ, L ) v ( r, λ, L ) .• The magnetic field strength or direction: | B | : S → R , B : S → R .• The gravity gradient tensor: Γ : S → R × .It’s important to note that the standard form of the Birkhoff-Kinchin Theorem assumes that the func-tion being averaged can be made to depend exclusively on the spatial state of the system. However,many functions of interest are not necessarily time independent. For example, both the Earth’s at-mospheric density and magnetic field fluctuate over time. In cases such as this, mean values for eachlocation in S can be considered as an approximation - although this will be application dependentand requires further validation. Circular Orbits
Applying the Birkhoff-Kinchin theorem, we can straightforwardly write down an integral formulafor the long term time mean of any measurable function f of the position of the satellite: E µ c ( f ) = i (cid:90) − i π (cid:90) − π µ c ( λ, L ) f ( r, λ, L ) dLdλ (21)It then becomes a triviality to compute higher order quantities such as variance:Var µ c ( f ) = i (cid:90) − i π (cid:90) − π µ c ( λ, L ) ( f ( r, λ, L ) − E µ c ( f )) dLdλ (22)11With the square taken entry-wise if f isn’t scalar-valued)Also note that because the probability measure has no dependence on the longitude, if f does aswell, then we can directly evaluate the inner-most integral in each case, reducing these to singleintegrals. A similar approach can be used for latitude in some cases as well. Elliptical Orbits
Doing the same as in the previous section, but for the elliptical measure: E µ e ( f ) = a (1+ e ) (cid:90) a (1 − e ) i (cid:90) − i π (cid:90) − π µ ( r, λ, L ) f ( r, λ, L ) dLdλdr (23)We can compute the variance in this case as well:Var µ e ( f ) = a (1+ e ) (cid:90) a (1 − e ) i (cid:90) − i π (cid:90) − π µ ( r, λ, L ) ( f ( r, λ, L ) − E µ ( f )) dLdλdr (24)Just as in the previous section, it will often be possible to reduce these volume integrals to single ordouble integrals by directly evaluating one or more of the integrals in the expression. There can besignificant benefit to doing so, because accurate evaluation of volume integrals can be a numericallycostly operation.In many cases, it’s possible to evaluate these formulae in closed form. For example, we can applyEquation 23 to easily compute the mean radius of of an elliptical orbit with respect to time: E µ ( r ) = a (1+ e ) (cid:90) a (1 − e ) i (cid:90) − i π (cid:90) − π µ ( r, λ, L ) rdLdλdr = aπ π (cid:90) θ = − π (1 − e sin( θ )) dθ = a (cid:18) e (cid:19) Which is in agreement with the standard formula. NUMERICAL RESULTS
In this section, numerical accuracy results for some of the preceding formulae are provided. Ingeneral, the “ground truth” we compare these numerical results against is the value calculated viadirect trajectory propagation on an RK-78 integrator, using the J -perturbed force model providedin Equation 2. The absolute error and percent error metrics used here are defined in equation 25.Absolute Error := | True Value − Estimate | Percent Error := | True Value − Estimate | True Value ·
100 =
Absolute ErrorTrue Value · (25)Both of these metrics are utilized throughout this section, but primarily absolute error. To illustratethe reason for this, we consider the “One Satellite, One Ground Station” view period ratio case.Here, note that each view period ratio will be a number ρ ∈ [0 , . . If ρ = . , for example,this would imply that the satellite and ground station can communicate approximately of thetime. Now, suppose that the estimate given by Equation 14 were . . Then, the absolute error = | . − . | = . has a very natural interpretation: the view period ratio estimate from theformula is off by of total flight time. Similarly, the percent error = | . − . | . ·
100 = 4% hasthe interpretation: the view period ratio estimate from the formula deviates by from the trueview period ratio. In this case, each of these work as effective metrics. However, consider a slightly12ifferent case, where ρ true = . , and ρ estimate = . . Then, the absolute error = | . − . | = . , and the percent error = . . ·
100 = 100% . In this scenario, we can still interpret these valuesin the same way as before, however while the true and estimate values are far closer than in the firstcase, the percent error metric represents ρ estimate as a very poor estimate. This issue only worsensas ρ true → , and if ρ true = 0 , then the percent error isn’t well-defined. Thus, although there may becases where percent error is a useful metric, absolute error is preferred here due to its robustness. One Satellite, One Ground Station
The numerical results given below demonstrate the distribution of the absolute error of the viewperiod ratio formulae given by equations 12 and 14 when compared against the method of directorbit propagation for 6000 (integration) days, using the J -perturbed Earth model. Each of the50,000 cases given were randomly sampled from the following sample space:Semimajor Axis a ∈ [6371 . , (km)Eccentricity e ∈ [0 , . (unitless)Orbital Inclination i ∈ [0 , (deg)Ground Station Latitude g ∈ [0 , (deg)Elevation Angle (cid:15) ∈ [0 , (deg)Field of View Angle β ∈ [0 , (deg)In addition, we require that the apogee is at least 160km above the surface. Figure 6(a)
A histogram of the absolute error, defined inEquation 25, for the 50,000 cases assigned above, with ver-tical lines denoting the values of absolute error such that , , and of the cases have error less thanthat value. In addition, cases with error ≥ . account for . of the total of cases. Figure 6(b)
The same as
Figure 6(a) , the difference beingthat the cases with orbital inclination within . ◦ are fil-tered out (amounting to 3.408% of the , cases). No-tice that with these critically (or nearly critically) inclinedorbits filtered out, each of the percentiles of cases occursearlier. These differences are most pronounced for the and thresholds. igure 7 : A scatter plot of the absolute error, defined in Equation 25, for each of the 50,000 cases asa function of the inclination of the orbit. Note the significant increase in error at critical inclination,which occurs at inclination ≈ . ◦ . Critical inclination is the orbit inclination at which thesatellite orbit experiences zero apogee drift, which constitutes a degenerate case for the view periodformula, so this was to be expected. Figure 6 illustrates the effect removing these cases has on theerror profile. Figure 8 : Similar to the plots in Figure 6, showing percent error, defined in Equation 25. Caseswhere ρ true = 0 , but ρ estimate (cid:54) = 0 were handled by setting them to error. This occurred in of the , cases (or .68%). 14 ONCLUSION & FUTURE WORK
In this work, we provided an analytical method of estimating the long-term mean value of anyfunction of satellite position, under the assumption of a J − perturbed and aperiodic orbit. Spe-cial emphasis was placed on applying this approach to rapid ground coverage assessment, withintegral formulae optimized for numerical computation provided in Equations 12, 14, 15 and 16.In particular, formulae for computing coverage while taking into account ground station elevationangle, satellite field of view, and multiple ground stations for both circular and elliptical satelliteorbits were provided. Note that these formulae are really just a special case of the general methodgiven in Equations 21, 23, 22 and 24 for estimating the time mean and variance of any function ofthe satellite trajectory. We anticipate that this will enable accelerated evaluation of many relevantflight parameters such as velocity, drag, magnetic field strength/direction, gravity gradient and sun β − angle, among other quantities of interest. It’s worth noting that each of these formulae are com-pletely independent of any physical constants, other than the radius of the body R B , and (implicitly) J , which is necessary to determine the periodicity (or the lack thereof) of a given orbit.The feasibility of extending this approach to quantifying the dynamics of multiple satellites wasalso investigated. The main barrier to this extension was the requirement that the satellite trajectoriesbe statistically independent - a property which is unlikely to hold for most satellite constellationsand formations, in which the orbits are often commensurate by design. The authors are activelyworking on extending the underlying theory to more effectively handle these cases. It’s anticipatedthat this can be accomplished without significant increases in complexity or computational cost.As was alluded to in the “Many Satellites, One (or Many) Ground Station(s)” section, it can beinstructive to view, µ c and µ e as probability densities on the instantaneous position of the satellite.On this interpretation, it becomes natural to consider quantities such as expectation and variance. Inaddition, it suggests applications to orbit determination via maximum likelihood estimation (e.g. forsituational awareness or exoplanet TLE estimation) and the utilization of Chebyshev-like bounds onthe distribution of values of functions of satellite position. However further research will be neces-sary to accurately assess the efficacy of these and the limitations of the probabilistic interpretationof µ c and µ e .Throughout this paper, a linearized J model was assumed. The J − perturbed two body problemmay be well-approximated by a linear precession of the orbital elements Mean Anomaly, Argumentof Periapsis, and Longitude of the Ascending Node. This linearization is a key step in the derivationof the results provided here. However, because this linearization only approximately captures thebehavior of the system, with bias towards the secular effects, there are classes of orbits for whichour method performs rather poorly. The most notable of which is the class of critically inclinedorbits. As we saw in the numerical results section, there’s a sharp increase in view period ratioestimation error for orbits which are near critical inclination ( ≈ . ◦ ). Critically inclined orbitsoften exhibit pathological behavior in their own right, but our model introduces additional inaccu-racies. In particular, in order for the ergodicity assumption on an orbit to be satisfied, it must beaperiodic which, in this context, means that the three linear precession rates ( ˙ M , ˙ ω, ˙Ω ) are rationallyindependent * . While such orbits constitute a set of measure zero, near-periodic orbits tend to exhibitreduced accuracy. Critical inclination is a special case of periodicity in which ˙Ω = 0 , resulting inparticularly poor accuracy in its vicinity. The use of a higher order (non-linearized) model wouldlikely ameliorate many of these accuracy issues by capturing the effect of critical inclination. * Specifically, (cid:64) α ∈ Q such that ˙ M ˙ ω = α or ˙ ω ˙Ω = α or ˙Ω˙ M = α CKNOWLEDGMENTS
This research was carried out in part at the Jet Propulsion Laboratory, California Institute of Tech-nology under a contract with the National Aeronautics and Space Administration (80NM0018D0004).This work was sponsored in part by the Caltech Summer Undergraduate Research Fellowship Pro-gram. This work was also supported in part by the Hummer-Tuttle gift to Professor Al Barr throughthe Caltech Division of Engineering and Applied Science.
NOTATION ρ = The View Period Ratio, ρ ∈ [0 , µ = The Standard Gravitational Parameter µ c = The Invariant Measure for Circular Orbits µ e = The Invariant Measure for Elliptical Orbits a = Orbit Semimajor Axis e = Orbit Eccentricity i = Orbit Inclination Angle M = Mean Anomaly Ω = Longitude of the Ascending Node ω = Argument of Periapsis r = The Radial Component of Satellite Position, r ∈ ( a (1 − e ) , a (1 + e )) λ = The Latitude Component of Position, λ ∈ [ − i, i ] L = The Longitude Component of Position L ∈ [ − π, π ] rad or L ∈ [ − , ◦ R B = The Radius of the Central Body J = The Second Zonal Harmonic g = The Latitude of the Ground Station g ∈ ( − π , π ) (cid:15) = The Elevation Angle of the Ground Station, (cid:15) ∈ (0 , π ) β = The Field of View (FOV) Angle of the Satellite, β ∈ (0 , π ) θ = The Ground Station Mask Radius S = The State Space of the Satellite, S = [ a (1 − e ) , a (1 + e )] × [ − i, i ] × [ − π, π ]Ω B = The Precession Rate of the Central BodyE µ c ( f ) = The Expected Value of f With Respect to the Probability Measure µ c For Circular OrbitsE µ e ( f ) = The Expected Value of f With Respect to the Probability Measure µ e For Elliptical OrbitsVar µ c ( f ) = The Variance of f With Respect to the Probability Measure µ c For Circular OrbitsVar µ e ( f ) = The Variance of f With Respect to the Probability Measure µ e For Elliptical OrbitsTrunc b ( x ) = The Truncation Function, Trunc b ( x ) = min { max { x, − b } , b } APPENDIX A
Note that we can’t just compute the naive view period ratio for each ground station and sum thembecause the ground station masks may overlap, as in the Figure 5, resulting in double counting someregions. We will need to do a bit of additional work to avoid that issue. Firstly, note that the naiveview period integral for ground station k with ground station mask radius θ k is given by: ρ = (cid:90) Trunc i ( λ k + θ k ) λ = Trunc i ( λ k − θ k ) (cid:90) L k +cos − ( cos( θk ) − sin( λ ) sin( λk )cos( λ ) cos( λk ) ) L = L k − cos − ( cos( θk ) − sin( λ ) sin( λk )cos( λ ) cos( λk ) ) cos( λ )2 π (cid:112) sin ( i ) − sin ( λ ) dLdλ
16n particular for each value of λ in range, we want to integrate over the values of L in the interval: [ L k, min ( λ ) , L k, max ( λ )] := (cid:20) L k − cos − (cid:18) cos( θ k ) − sin( λ ) sin( λ k )cos( λ ) cos( λ k ) (cid:19) , L k + cos − (cid:18) cos( θ k ) − sin( λ ) sin( λ k )cos( λ ) cos( λ k ) (cid:19)(cid:21) Thus, when there are multiple ground stations, we want to integrate over the union of the intervals: L range ( λ ) = N (cid:91) k =1 [ L k, min ( λ ) , L k, max ( λ )] The simplest way to work with this union of intervals programmatically is to combine them into aunion of disjoint intervals. There is a simple algorithm for doing so.
Algorithm: Interval Merge Given a set of N intervals: [ x , y ] , [ x , y ] , . . . , [ x N , y N ] :1. Reorder the intervals in increasing order based on the lower bounds of the intervals. That is: [ x , y ] , [ x , y ] , . . . , [ x N , y N ] → [ x (cid:48) , y (cid:48) ] , [ x (cid:48) , y (cid:48) ] , . . . , [ x (cid:48) N , y (cid:48) N ] , s.t: x (cid:48) ≤ x (cid:48) ≤ . . . ≤ x (cid:48) N
2. Push the first interval onto the stack.3. For each interval in the ordered list:(a) If the current interval does not overlap with the interval on the top of the stack, push itonto the top.(b) If the current interval overlaps with stack top and ending time of current interval is morethan that of stack top, update stack top with the ending time of current interval.4. Return the new list of intervals [ a , b ] , [ a , b ] , . . . , [ a M , b M ] , M ≤ N .So, we can write the pseudocode for the approach as follows:We are given a satellite with semimajor axis a and orbit inclination i , and a set of ground stationswith latitude-longitude coordinates ( λ , L ) , ( λ , L ) , . . . , ( λ N , L N ) and ground station mask radii θ , θ , . . . , θ N .Then, for each λ , we can write down the interval giving the range of values of L for each groundstation: [ L , min ( λ ) , L , max ( λ )] , [ L , min ( λ ) , L , max ( λ )] , . . . , [ L N, min ( λ ) , L N, max ( λ )] Now, there are two steps of preprocessing we need to do:1. Remove all empty intervals (ie. intervals corresponding to ground stations which have novisibility at the latitude λ .2. Some intervals may wrap around past π . Break these up into two intervals, wrapping theinterval around: [ L k, min ( λ ) , L k, max ( λ )] → [ L k, min ( λ ) , π ] , [ − π, L k, max ( λ ) − π ] (also do the analogous for intervals wrapping past − π )17abel this new set of intervals: [ x , y ] , [ x , y ] , . . . , [ x N (cid:48) , y N (cid:48) ] Then run the interval merge algo-rithm on these to yield the new list of ordered disjoint intervals: I ( λ ) := { [ a , b ] , [ a , b ] , . . . , [ a M , b M ] } Furthermore, let: λ min = Trunc i (min k { λ k − θ k } ) , λ max = Trunc i (min k { λ k + θ k } ) .Then, we can write the view period integral as: ρ = (cid:90) λ max λ = λ min (cid:88) [ a k ,b k ] ∈ I ( λ ) (cid:90) b k L = a k cos( λ )2 π (cid:112) sin ( i ) − sin ( λ ) dLdλ = (cid:90) λ max λ = λ min (cid:88) [ a k ,b k ] ∈ I ( λ ) cos( λ )( b k − a k )2 π (cid:112) sin ( i ) − sin ( λ ) dλ = (cid:90) λ max λ = λ min cos( λ )2 π (cid:112) sin ( i ) − sin ( λ ) (cid:88) [ a k ,b k ] ∈ I ( λ ) ( b k − a k ) dλ = 12 π (cid:90) λ max λ = λ min cos( λ ) Length ( I ( λ )) (cid:112) sin ( i ) − sin ( λ ) dλ (where Length ( I ( λ )) := (cid:80) [ a k ,b k ] ∈ I ( λ ) ( b k − a k ) , the total length of the union of the intervals)Finally, we take the change of variables from Equation 11 to remove any singularities from theintegrand, yielding Equation 15. REFERENCES [1] V. Arnold,
Ergodic Problems of Classical Mechanics . New York: Addison-Wesley, 1989.[2] A. Graven and M. W. Lo, “The Long-Term Forecast of Station View Periods for Elliptical Orbits,”
AASAstrodynamics Specialist Conference , 2019. https://arxiv.org/abs/2010.06021.[3] W. McClain and D. Vallado,
Fundamentals of Astrodynamics and Applications . Space Technology Li-brary, Springer Netherlands, 2001.[4] Y. G. Sinai,
Introduction to Ergodic Theory . Princeton University Press, Princeton N.J., 1976.[5] M. W. Lo, “The Long-Term Forecast of Station View Periods,” tech. rep., Jet Propulsion Laboratory,Pasadena, California, 1994.[6] S. K. Stein, “”Mean Distance” in Kepler’s Third Law,”
Mathematics Magazine , Vol. 50, No. 3, 1977,pp. 160–162.[7] I. Gkolias, J. Daquin, F. Gachet, and A. J. Rosengren, “From Order to Chaos in Earth Satellite Orbits,”
The Astronomical Journal , Vol. 152, oct 2016, p. 119, 10.3847/0004-6256/152/5/119.[8] R. C. Enaganti,
Merge Overlapping Intervals , 2020 (accessed 2020-07-15). ..