Satellite Hacking
Disclaimer: This guide is for educational purposes only. All content is presented to advance cybersecurity knowledge and understanding of satellite systems. Always prioritize ethical considerations and legal compliance when exploring these concepts.
This section explores how external forces (J2, atmospheric drag, solar radiation pressure, and third-body effects) cause deviations from ideal Keplerian orbits. Covers mathematical frameworks including Lagrange's equations, semi-analytical methods (DSST), and numerical integration techniques. Includes complete SGP4 propagation implementation with real ISS TLE data and practical applications for mission planning and debris tracking.
Overview: Perturbation theory is a mathematical framework that analyzes how external forces (perturbations) affect satellite orbits. Rather than treating satellite motion as a perfect Kepler orbit, perturbation theory accounts for real-world forces that cause deviations from ideal orbital mechanics. This is essential for accurate satellite propagation, mission planning, and collision avoidance.
Fundamental Concept: In perturbation theory, the satellite's motion is treated as a Keplerian orbit with small corrections applied. The approach uses the two-body problem as the base solution and adds acceleration terms representing perturbative forces. The equations of motion can be expressed as: a_total = a_kepler + a_perturbations
Key Perturbations and Their Effects:
Earth's Oblateness (J2 Perturbation): Earth is not a perfect sphere but an oblate spheroid (flattened at the poles). The J2 term, representing the quadrupole moment, is the dominant perturbation for most satellites. It causes secular (long-term) changes including apsidal precession (rotation of the line of apsides), nodal precession (rotation of the orbital plane), and mean motion rate changes. For a satellite at 400 km altitude, J2 causes a nodal precession of approximately 3.2 degrees per day. The magnitude of J2 is approximately -1.08263e-3 for Earth.
Higher-Order Harmonics (J3, J4, etc.): Beyond J2, higher zonal harmonics cause additional perturbations. J3 affects the argument of perigee precession, while J4 influences mean motion. These effects are smaller than J2 but critical for high-precision applications. Modern models use terms up to J20 or higher.
Atmospheric Drag: Satellites in low Earth orbit (LEO, typically below 2000 km) experience aerodynamic drag from the thermosphere. The drag force depends on satellite velocity (quadratically), atmospheric density (exponentially dependent on altitude and solar activity), satellite cross-sectional area, and drag coefficient (typically 2-3). The atmospheric density varies with solar flux, geomagnetic storms, and time of day. For a 1000 kg satellite at 400 km altitude, atmospheric drag can cause an orbital decay of 1-2 km per day during solar maximum.
Solar Radiation Pressure (SRP): Photons from the Sun exert momentum transfer to satellites. The force depends on solar flux (approximately 1361 W/m² at Earth's distance), satellite cross-sectional area, surface reflectivity (α), and material absorption properties. For a satellite with 10 m² cross-section and high reflectivity, SRP acceleration is on the order of 10⁻⁶ m/s², which can accumulate significantly over time. SRP is particularly significant for high-area-to-mass ratio objects (such as solar sails or large antenna structures).
Third-Body Effects (Moon and Sun): The gravitational attraction of the Moon and Sun on satellites creates perturbative forces. The magnitude depends on the satellite's distance from these bodies. For satellites in geostationary orbit (GEO), third-body effects cause oscillations in orbital elements with periods of several days. The combined effect of lunar and solar gravity can cause libration (oscillatory motion) in GEO satellites, requiring station-keeping maneuvers approximately every 2-3 days.
Tidal Forces: Differential gravitational forces from the Moon and Sun across the satellite's extent cause tidal effects. While typically small for small satellites, these forces become significant for extended bodies or large space stations. Tidal torques can cause tumbling in large structures.
Relativistic Effects: For precision applications (like atomic clock synchronization), general relativity causes time dilation effects. GPS satellites must account for relativistic frequency shifts of approximately 38 microseconds per day to maintain accuracy.
Mathematical Treatment: Perturbation theory employs several mathematical approaches:
Lagrange's Planetary Equations: These equations express the time derivatives of orbital elements (semi-major axis, eccentricity, inclination, etc.) in terms of perturbative accelerations. This approach is particularly useful for long-term propagation and understanding how perturbations affect orbital element evolution.
Delaunay Variables: These canonical variables simplify the mathematical treatment of perturbations by separating fast (orbital motion) and slow (secular) terms.
Techniques for Analysis:
Analytical Methods: Lagrange's Planetary Equations, Delaunay perturbation theory, and averaging techniques are used for simplified models and long-term trend analysis. These methods provide closed-form solutions or series expansions that reveal the structure of orbital evolution. The advantages include computational efficiency and physical insight, though accuracy decreases for complex multi-perturbation scenarios.
Semi-Analytical Methods: These combine analytical averaging with numerical integration to balance accuracy and computational efficiency. Examples include the Draper Semi-analytical Satellite Theory (DSST), which uses averaging over one orbital period while the osculating elements are integrated numerically.
Numerical Integration Methods: Direct numerical integration of the equations of motion (such as Runge-Kutta, Adams-Moulton, or symplectic integrators) using Cartesian coordinates accounts for all perturbations simultaneously and provides high accuracy. These methods handle complex perturbation scenarios with high precision but are computationally expensive for long-term propagation.
Practical Applications:
Perturbation theory has critical applications across space operations. LEO satellite mission planning benefits from accurate orbital decay prediction, allowing operators to schedule maintenance maneuvers efficiently. GEO satellite station-keeping maneuver planning requires precise understanding of third-body and solar radiation pressure effects to minimize fuel consumption. Space debris tracking and collision avoidance depend on perturbation theory to maintain accurate orbital catalogs. GPS and navigation satellite constellation management uses perturbation theory to maintain precise orbital geometry for positioning accuracy. Deep space mission orbit design leverages perturbation theory to optimize trajectories for fuel efficiency and mission objectives.
SGP4 Propagation with TLE Data (Industry Standard):
The SGP4 (Simplified General Perturbations 4) model is the industry standard for satellite propagation. It accounts for major perturbations and uses Two-Line Element (TLE) sets as input. TLE data is publicly available from NORAD for all tracked objects.
from skyfield.api import load, EarthSatellite
from datetime import datetime, timedelta
import numpy as np
# Load satellite data from TLE (International Space Station)
tle_line1 = '1 25544U 98067A 25032.85548342 .00001539 00000-0 33262-4 0 9999'
tle_line2 = '2 25544 51.6443 75.8703 0003838 32.7176 340.6461 15.49294176311377'
# Create satellite object
sat = EarthSatellite(tle_line1, tle_line2)
# Load ephemeris for Earth position
ts = load.timescale()
ephem = load('de421.bsp')
earth = ephem['earth']
# Define time range for propagation
start = ts.utc(2025, 2, 8, 0, 0, 0)
end = ts.utc(2025, 2, 8, 23, 59, 59)
# Propagate orbit (accounts for J2, atmospheric drag, SRP)
t = ts.linspace(start, end, 1440) # One point per minute
astrometric = earth.at(t).observe(sat)
ra, dec, distance = astrometric.apparent().radec()
# Extract position in ECEF coordinates
geo = earth.at(t).observe(sat).apparent()
# Analysis: Calculate orbital decay due to atmospheric drag
initial_altitude = distance[0]
final_altitude = distance[-1]
altitude_decay = initial_altitude - final_altitude
print(f"Initial Altitude: {initial_altitude:.2f} km")
print(f"Final Altitude: {final_altitude:.2f} km")
print(f"Orbital Decay (24 hours): {altitude_decay:.4f} km")
Understanding TLE Format:
A Two-Line Element set contains orbital elements in a specific format used by NORAD. Line 1 provides the satellite number, epoch year and day, drag coefficients (B* value), and checksum information. Line 2 contains inclination, Right Ascension of Ascending Node, eccentricity, argument of perigee, mean anomaly, mean motion, and revolution number at epoch.
Key Parameters in SGP4/TLE:
B* (Ballistic Coefficient): Accounts for atmospheric drag, expressed in 1/Earth radii. Higher values indicate stronger drag effect. Mean Motion: Orbital revolutions per day, approximately 15.5 for ISS. Mean Anomaly: Position of satellite in orbit, ranging from 0° to 360°.
Detailed analysis of Earth's non-uniform gravity field using spherical harmonic expansion (zonal, tesseral, and sectorial harmonics). Explains the EGM2008 gravity model with 2,190 coefficients. Covers computational methods for calculating harmonic accelerations using Legendre polynomials and recursive algorithms. Essential for meter-level precision in GPS, space debris tracking, and orbital determination for satellites across all altitude regimes.
Fundamental Concept: Earth's gravitational field is not uniform due to irregular mass distribution within the planet. Continents, ocean basins, mantle density variations, and the core create localized gravitational anomalies. Gravitational harmonics provide a mathematical framework to model these variations and their effects on satellite orbits.
Why Harmonics Matter: While the point-mass Earth approximation works reasonably well for rough calculations, it introduces significant errors (on the order of meters to kilometers) in satellite orbit predictions over extended periods. For applications like GPS positioning (requiring meter-level accuracy), space debris tracking, and orbital determination, accounting for Earth's gravitational harmonics is essential.
Spherical Harmonic Expansion: The gravitational potential V at a point can be expressed as an infinite series expansion in spherical harmonics representing Earth's non-uniform gravity field.
Zonal vs. Tesseral Harmonics:
Zonal Harmonics (m=0): These depend only on latitude and represent the axisymmetric part of Earth's gravitational field. The J2 term (related to C20) is the dominant zonal harmonic, causing secular changes in orbital elements. Tesseral Harmonics (m≠0): These have both latitudinal and longitudinal dependence and represent regional mass anomalies (like mountain ranges or basin structures), causing periodic perturbations synchronized with satellite orbital motion. Sectorial Harmonics (m=n): Special cases where order equals degree, representing large-scale regional features.
Geopotential Models:
EGM96 (Earth Gravitational Model 1996): Expanded to degree and order 360, provides approximately 55 km spatial resolution. EGM2008: Expanded to degree 2190, provides approximately 5 km resolution, incorporating gravity data from multiple satellite missions (GRACE, GOCE, etc.) and ground measurements. This is the current standard for high-precision geoid modeling. GRACE/GOCE Models: Derived from dedicated gravity measurement satellites where GOCE measured gravity gradients and GRACE measures inter-satellite distances to infer gravity variations.
Effects on Satellite Orbits and Resonance:
Harmonic perturbations cause perturbations in orbital elements, precession of the orbit plane, resonance effects when orbital frequencies match harmonic frequencies, and geosynchronous orbit stability issues. When satellite orbital characteristics satisfy certain relationships with Earth's rotation and harmonic frequencies, resonance effects can occur. Some satellites naturally maintain certain orbital parameters (frozen orbits) due to harmonic resonance.
import numpy as np
from scipy.special import eval_lpmn
from scipy.integrate import solve_ivp
"""Computing gravitational acceleration from harmonic coefficients"""
def calculate_harmonic_acceleration(r_ecef, C, S, max_degree=20):
"""Calculate acceleration due to sphere harmonics"""
# Constants
GM = 3.986004418e14 # Earth's gravitational constant (m³/s²)
ae = 6378137.0 # Earth equatorial radius (m)
# Convert to spherical coordinates
r_norm = np.linalg.norm(r_ecef)
lat = np.arctan2(r_ecef[2], np.sqrt(r_ecef[0]**2 + r_ecef[1]**2))
lon = np.arctan2(r_ecef[1], r_ecef[0])
# Calculate potential derivatives (simplified)
a = np.zeros(3)
for n in range(2, max_degree + 1):
for m in range(0, n + 1):
# Legendre polynomial calculation
P_nm = eval_lpmn(m, n, np.sin(lat))[0][-1]
# Accumulate acceleration components
factor = (GM / r_norm**2) * (ae/r_norm)**n
factor *= (n + 1) * P_nm
a[0] += factor * (C[n,m] * np.cos(m*lon))
return a
Comprehensive integration of all perturbation effects (J2, J3, atmospheric drag, SRP, third-body) into a unified satellite dynamics simulator. Demonstrates simultaneous propagation accounting for interacting perturbations and cumulative orbital evolution effects. Includes real-world ISS propagation scenarios with accurate decay rate prediction and long-term orbital analysis for mission design and constellation management.
Comprehensive Framework: Real-world satellite motion cannot be accurately predicted using a single perturbation. Advanced orbit determination requires combining all major perturbation effects simultaneously: Earth's gravitational harmonics (J2 and higher-order terms), atmospheric drag (altitude and solar activity dependent), solar radiation pressure, third-body gravity (Moon, Sun), relativistic effects, and tidal forces.
Orbit Determination and Propagation: Given initial state vectors or Two-Line Element (TLE) data, modern orbit determination uses numerical integration to propagate satellite position and velocity accounting for all perturbations. Observations from ground stations or space-based sensors are processed using least-squares fitting or Kalman filtering to refine the orbit estimate iteratively. Accuracy of ±10 meters or better is achievable for well-tracked satellites.
Real-World Mission Examples:
The International Space Station (ISS) operates as a LEO satellite at approximately 408 km altitude with 51.6° inclination, requiring severe atmospheric drag accounting with losses of approximately 100 meters altitude daily and periodic reboost maneuvers every 2-3 months, while J2 perturbations cause orbital plane precession. The GPS Constellation consists of 24-32 satellites in semi-synchronous orbits at approximately 20,200 km altitude on 12-hour periods, requiring station-keeping maneuvers every few weeks to correct for solar radiation pressure and lunar gravity perturbations, with relativistic time dilation effects (38 microseconds per day) needing correction for positioning accuracy. Geostationary Satellites including weather satellites (GOES) and communication satellites are maintained at 35,786 km altitude, requiring approximately 2-4 m/s delta-v annually for station-keeping, with libration points naturally emerging due to Earth's non-uniform gravity field and requiring active control. Sun-Synchronous Orbit satellites like Landsat and Sentinel operate at 700-800 km, exploiting J2 perturbation's nodal precession to maintain constant sun angle relative to the orbit plane, enabling consistent imaging geometry.
Comprehensive Python Implementation:
import numpy as np
from scipy.integrate import solve_ivp
from scipy.special import eval_lpmn
from datetime import datetime, timedelta
"""Advanced Satellite Dynamics Propagator with Multi-Perturbation Acceleration Model"""
class SatelliteDynamicsSimulator:
def __init__(self):
"""Initialize constants for orbital mechanics"""
# Earth parameters
self.GM = 3.986004418e14 # m³/s²
self.ae = 6378137.0 # Earth equatorial radius (m)
self.J2 = -1.08263e-3 # Second dynamic form factor
self.J3 = 2.54e-6 # Third dynamic coefficient
self.omega_earth = 7.2921150e-5 # Earth rotation rate (rad/s)
# Atmospheric model parameters (simple exponential)
self.rho_0 = 1.225 # Reference density at sea level (kg/m³)
self.H_scale = 8500 # Scale height (m) - varies 6-9 km
# Sun and Moon mass parameters
self.GM_sun = 1.32712440018e20 # m³/s²
self.GM_moon = 4.9028e12 # m³/s²
# Radiation pressure coefficient
self.c_r = 1.2 # Reflectivity (typically 1.2-1.5)
def atmospheric_density(self, r_ecef):
"""Calculate atmospheric density at altitude using exponential model"""
altitude = np.linalg.norm(r_ecef) - self.ae
rho = self.rho_0 * np.exp(-altitude / self.H_scale)
return rho
def acceleration_j2_j3(self, r_ecef):
"""J2 and J3 perturbations (dominant zonal harmonics)"""
r = np.linalg.norm(r_ecef)
z = r_ecef[2]
rho_sq = r_ecef[0]**2 + r_ecef[1]**2
factor_j2 = (3.0/2.0) * self.GM * self.J2 * (self.ae**2) / (r**5)
factor_j3 = (5.0/2.0) * self.GM * self.J3 * (self.ae**3) / (r**7)
a_j2 = factor_j2 * np.array([
5*z*r_ecef[0]/rho_sq - r_ecef[0],
5*z*r_ecef[1]/rho_sq - r_ecef[1],
5*z*z/rho_sq - 3*r_ecef[2]
])
a_j3 = factor_j3 * np.array([
5*z*(7*z*z/rho_sq - 3)*r_ecef[0]/rho_sq,
5*z*(7*z*z/rho_sq - 3)*r_ecef[1]/rho_sq,
(35*z**3/rho_sq - 15*(z**2))*r_ecef[2]/rho_sq
])
return a_j2 + a_j3
def acceleration_drag(self, r_ecef, v_ecef, mass, area):
"""Atmospheric drag acceleration"""
# Wind-relative velocity (simple model, ignoring winds)
v_wind = -self.omega_earth * np.array([-r_ecef[1], r_ecef[0], 0])
v_rel = v_ecef - v_wind
v_mag = np.linalg.norm(v_rel)
# Ballistic coefficient
rho = self.atmospheric_density(r_ecef)
C_d = 2.5 # Drag coefficient
if v_mag < 1e-6:
return np.array([0, 0, 0])
a_drag = -0.5 * rho * C_d * (area / mass) * v_mag * v_rel
return a_drag
def acceleration_srp(self, r_ecef, r_sun_ecef, mass, area, cr=None):
"""Solar Radiation Pressure acceleration"""
if cr is None:
cr = self.c_r
# Shadow function (simplified - no eclipse calculation)
r_sun_norm = np.linalg.norm(r_sun_ecef - r_ecef)
a_shadow = 1.0 # No shadow = 1, in shadow = 0
# SRP constant: ~4.56e-6 N/m²
P_rad = 4.56e-6
u_sun = (r_sun_ecef - r_ecef) / r_sun_norm
a_srp = -a_shadow * (P_rad * cr * area / mass) * u_sun
return a_srp
def acceleration_thirdbody(self, r_ecef, r_body_ecef, GM_body):
"""Third-body perturbation (Moon, Sun)"""
r_rel = r_body_ecef - r_ecef
r_mag = np.linalg.norm(r_rel)
# Direct perturbation (satellite's acceleration toward body)
a_direct = GM_body * r_rel / (r_mag**3)
# Indirect effect (Earth's acceleration toward body)
r_body_mag = np.linalg.norm(r_body_ecef)
a_indirect = -GM_body * r_body_ecef / (r_body_mag**3)
return a_direct + a_indirect
def acceleration_central(self, r_ecef):
"""Central gravity acceleration (point-mass Earth)"""
r = np.linalg.norm(r_ecef)
return -self.GM / (r**3) * r_ecef
def propagate_eci_acceleration(self, t, state, mass, area, r_sun_ecef, r_moon_ecef):
"""Total acceleration in ECI frame"""
r_ecef = state[:3]
v_ecef = state[3:]
# Central gravity
a_total = self.acceleration_central(r_ecef)
# Harmonic perturbations
a_total += self.acceleration_j2_j3(r_ecef)
# Non-gravitational perturbations
a_total += self.acceleration_drag(r_ecef, v_ecef, mass, area)
a_total += self.acceleration_srp(r_ecef, r_sun_ecef, mass, area)
# Third-body perturbations
a_total += self.acceleration_thirdbody(r_ecef, r_sun_ecef, self.GM_sun)
a_total += self.acceleration_thirdbody(r_ecef, r_moon_ecef, self.GM_moon)
return np.concatenate([v_ecef, a_total])
def propagate(self, initial_state, mass, area, duration_hours, step_minutes=1):
"""Propagate satellite orbit over time period"""
t_span = (0, duration_hours * 3600)
t_eval = np.linspace(0, duration_hours * 3600, int(duration_hours * 60 / step_minutes))
# Approximate Sun and Moon positions (simplified ephemeris)
# For production, use pyephem or astropy
r_sun = np.array([1.496e11, 0, 0]) # Approximate AU
r_moon = np.array([3.844e8, 0, 0]) # Approximate lunar distance
def derivs(t, y):
return self.propagate_eci_acceleration(
t, y, mass, area, r_sun, r_moon
)
# Integrate using RK45 method
solution = solve_ivp(
derivs, t_span, initial_state,
t_eval=t_eval, method='RK45',
dense_output=True, max_step=60
)
return solution
"""Example: ISS Orbit Propagation"""
simulator = SatelliteDynamicsSimulator()
# ISS initial state (approximate)
r_iss = np.array([6.756e6, 0, 0]) # 378 km altitude (perigee)
v_iss = np.array([0, 7680, 0]) # ~7.68 km/s circular velocity
initial_state = np.concatenate([r_iss, v_iss])
# ISS parameters
mass_iss = 420000 # kg
area_iss = 2500 # m² (effective cross-section)
# Propagate for 7 days
solution = simulator.propagate(
initial_state, mass_iss, area_iss,
duration_hours=7*24, step_minutes=5
)
# Extract final position and analyze decay
altitude_profile = np.linalg.norm(solution.y[:3], axis=0) - simulator.ae
initial_altitude = np.linalg.norm(initial_state[:3]) - simulator.ae
decay_rate = (initial_altitude - altitude_profile[-1]) / 7 # meters/day
print(f"Initial altitude: {initial_altitude:.0f} m")
print(f"Final altitude: {altitude_profile[-1]:.0f} m")
print(f"Decay rate: {decay_rate:.1f} m/day")
Conjunction Assessment & Probability of Collision (Pc) Calculation:
When debris approaches a satellite, operators must evaluate whether to perform an evasive maneuver. The decision hinges on the calculated probability of collision, which uses the Mahalanobis distance in conjunction with the covariance ellipsoid.
Probability of Collision Formula:
$$P_c = \frac{1}{\pi} \int_0^{\min(R^2, d^2)} \left(1 - \frac{r^2}{d^2}\right) dA$$
Where $R$ is the combined hard-body radius (satellite + debris), $d$ is the hard-body distance, and $A$ is the miss distance area.
For typical LEO conjunction assessment, evasive maneuvers are recommended when probability of collision exceeds $1 \times 10^{-4}$ (0.01%), become mandatory when exceeding $1 \times 10^{-3}$ (0.1%), with NASA's stringent threshold of $1 \times 10^{-5}$ (0.001%) applied for the International Space Station.
import numpy as np
from scipy.special import erf
from scipy.integrate import quad
class ConjunctionAssessment:
"""
Calculate probability of collision for debris avoidance
"""
def __init__(self, satellite_radius=2.5, debris_radius=0.1):
self.rs = satellite_radius # meters
self.rd = debris_radius
self.R = self.rs + self.rd # combined hard-body radius
def mahalanobis_distance(self, miss_vector, covariance_matrix):
"""
Calculate Mahalanobis distance: d = sqrt((x^T) * P^-1 * x)
"""
try:
cov_inv = np.linalg.inv(covariance_matrix)
d_sq = miss_vector @ cov_inv @ miss_vector
return np.sqrt(d_sq)
except:
return np.inf
def probability_of_collision(self, miss_distance, combined_covariance):
"""
Calculate probability of collision using 2D Gaussian model
"""
# Eigenvalue decomposition of covariance
eigvals, eigvecs = np.linalg.eigh(combined_covariance)
sigma1, sigma2 = np.sqrt(eigvals)
# Hard-body radius normalization
r_norm = self.R / np.sqrt(sigma1 * sigma2)
if r_norm > 3: # Beyond 3-sigma, probability is negligible
return 0.0
# Curtis & Srivastava formula
if sigma1 == 0 or sigma2 == 0:
return 0.0
# Precise formula for 2D case
pc = (1.0 / (2 * np.pi * sigma1 * sigma2)) * \
(self.R**2 / 2 +
sigma1**2 * np.exp(-self.R**2 / (2 * (sigma1**2 + sigma2**2))) +
sigma2**2 * np.exp(-self.R**2 / (2 * (sigma1**2 + sigma2**2))))
return max(0.0, min(1.0, pc))
def conjunction_analysis(self, relative_position, relative_velocity,
position_covariance, velocity_covariance,
time_to_ca):
"""
Full conjunction analysis
time_to_ca: time to closest approach (seconds)
"""
# Project covariance to time of closest approach
dt = time_to_ca
# Propagate covariance: P(t) = Phi*P(0)*Phi^T + Q
# Simplified: linear prediction
pos_cov_ca = position_covariance + \
(dt * relative_velocity[:, np.newaxis]) @ \
(dt * relative_velocity[:, np.newaxis]).T
# Miss distance at CA
miss_distance = np.linalg.norm(relative_position + relative_velocity * dt)
# Probability of collision
pc = self.probability_of_collision(miss_distance, pos_cov_ca)
return {
'pc': pc,
'miss_distance': miss_distance,
'time_to_ca': time_to_ca,
'maneuver_required': pc > 1e-4,
'maneuver_mandatory': pc > 1e-3
}
# Example
ca = ConjunctionAssessment()
# Debris approach scenario
rel_pos = np.array([500, -200, 150]) # relative position vector (meters)
rel_vel = np.array([-5, -3, 1]) # relative velocity vector (m/s)
pos_cov = np.eye(3) * (50**2) # 50m covariance
vel_cov = np.eye(3) * (0.5**2) # 0.5 m/s covariance
result = ca.conjunction_analysis(rel_pos, rel_vel, pos_cov, vel_cov, time_to_ca=600)
print(f"Probability of Collision: {result['pc']:.2e}")
print(f"Maneuver Required: {result['maneuver_required']}")
In-depth examination of CCSDS standards (TM, TC, AOS, USLP) governing spacecraft-to-ground communication. Analyzes security vulnerabilities including weak authentication, replay attacks, side-channel exploits, and signal spoofing. Presents secure CCSDS implementation with AES-256 encryption, HMAC-SHA256 authentication, sequence counters, and constant-time verification to protect against command injection and eavesdropping attacks.
CCSDS Protocols Overview:
The Consultative Committee for Space Data Systems (CCSDS) defines international standards for spacecraft communication. TM (Telemetry): Downlink data from satellite to ground station, typically encapsulated with error correction codes (Reed-Solomon). TC (Telecommand): Uplink commands from ground to satellite, requiring strong authentication and encryption. AOS (Advanced Orbiting Systems): Frame-based protocol for higher data rates, including virtual channel multiplexing. USLP (Unified Space Data Link Protocol): Modern protocol combining TM and TC with advanced security features.
Commercial Protocols:
DVB-S2/S2X: Used for satellite TV and broadband with vulnerabilities including BISS encryption weaknesses, signal piracy, and conditional access system attacks. Iridium Protocol: Proprietary TDMA system with historical vulnerabilities in early implementations that allowed unauthorized access. Starlink Protocol: Proprietary with AES-256 encryption, but limited public information due to security through obscurity.
Security Vulnerabilities and Attack Vectors:
Weak Authentication: Legacy CCSDS without HMAC-SHA256 used simple checksums instead of cryptographic signatures, allowing a 16-bit CRC to be brute-forced in microseconds. Replay Attacks: Commands without sequence numbers allow attackers to capture valid frames (such as "activate solar panels" TC frames) and replay them multiple times to cause power cycling or other harmful outcomes. Side-Channel Vulnerabilities: Timing attacks on AES implementations reveal key bits through analysis of execution time, while power analysis attacks measure current draw during cryptographic operations, both feasible against embedded satellite processors. Signal Spoofing: Transmitting false telecommand frames on the same frequency allows satellites without antenna directionality awareness to accept commands from "phantom" ground stations. Firmware Injection: Uploading malicious firmware via TC frames can reprogram the satellite's attitude control system, power management, or communication stack if integrity checking is insufficient. Denial of Service (DoS): Jamming or flooding the downlink with false telemetry frames overloads ground station processing capacity. Physical Layer Attacks: Bit-flip attacks through electromagnetic interference (EMI) during transmission exploit space environment radiation that naturally causes single-event upsets (SEUs), which sophisticated attacks can induce intentionally.
CCSDS Security Augmentation (CCSDS 352.0-B-1):
import hmac
import hashlib
import os
from cryptography.hazmat.primitives.ciphers import Cipher, algorithms, modes
from cryptography.hazmat.backends import default_backend
class CCSDSSecureProtocol:
"""
Secure CCSDS implementation with HMAC and AES
"""
def __init__(self, master_key, iv=None):
self.master_key = master_key # 32 bytes for AES-256
self.iv = iv if iv else os.urandom(16)
self.backend = default_backend()
self.sequence_counter = 0
def encrypt_telecommand(self, command_data, telecommand_code):
"""
Encrypt and authenticate a telecommand
Format: [TC_Code(2B)][Seq_Counter(4B)][IV(16B)][EncData][HMAC-SHA256(32B)]
"""
# Sequence counter prevents replay attacks
seq_bytes = self.sequence_counter.to_bytes(4, 'big')
self.sequence_counter += 1
# AES-256-CBC encryption
cipher = Cipher(
algorithms.AES(self.master_key),
modes.CBC(self.iv),
backend=self.backend
)
encryptor = cipher.encryptor()
# Pad command data to AES block size (16 bytes)
padding_len = 16 - (len(command_data) % 16)
padded_data = command_data + bytes([padding_len] * padding_len)
encrypted_data = encryptor.update(padded_data) + encryptor.finalize()
# Construct frame for HMAC
frame_for_auth = telecommand_code + seq_bytes + self.iv + encrypted_data
# HMAC-SHA256 for authentication
h = hmac.new(self.master_key, frame_for_auth, hashlib.sha256)
auth_tag = h.digest() # 32 bytes
# Final frame
tc_frame = frame_for_auth + auth_tag
return tc_frame
def decrypt_telemetry(self, tm_frame, expected_tm_code):
"""
Decrypt and verify telemetry with HMAC
"""
# Extract components
tm_code = tm_frame[0:2]
seq_counter = int.from_bytes(tm_frame[2:6], 'big')
iv = tm_frame[6:22]
encrypted_data = tm_frame[22:-32]
received_hmac = tm_frame[-32:]
# Verify code
if tm_code != expected_tm_code:
raise ValueError("Invalid telemetry code")
# Recalculate HMAC
frame_for_auth = tm_frame[:-32]
h = hmac.new(self.master_key, frame_for_auth, hashlib.sha256)
calculated_hmac = h.digest()
# Constant-time comparison (prevent timing attacks)
if not hmac.compare_digest(received_hmac, calculated_hmac):
raise ValueError("HMAC verification failed - possible tampering")
# Decrypt
cipher = Cipher(
algorithms.AES(self.master_key),
modes.CBC(iv),
backend=self.backend
)
decryptor = cipher.decryptor()
padded_data = decryptor.update(encrypted_data) + decryptor.finalize()
# Remove padding
padding_len = padded_data[-1]
data = padded_data[:-padding_len]
return {
'data': data,
'sequence': seq_counter,
'authenticated': True
}
# Example: Secure satellite command
protocol = CCSDSSecureProtocol(os.urandom(32))
# Command: Power on secondary radar (TC code 0x0A05)
command = b'\\x00\\x05' # Device ID
tc_frame = protocol.encrypt_telecommand(command, b'\\x0A\\x05')
print(f"Secure TC Frame (hex): {tc_frame.hex()}")
print(f"Frame size: {len(tc_frame)} bytes")
print(f"Authentication overhead: HMAC-SHA256 (32 bytes) + IV (16 bytes) + Seq (4 bytes) = 52 bytes")
Complete multi-layer ground station security framework covering RF reception, demodulation, protocol handling, and mission operations. Details threat models including supply chain compromise, rogue ground station attacks, and network lateral movement. Presents advanced defenses: geographic authentication, multi-signature command approval chains, ensemble anomaly detection (statistical + Isolation Forest + LSTM), microsegmentation with zero-trust architecture.
Typical Ground Station Components:
Ground stations comprise interconnected subsystems working in concert. Antenna Systems include parabolic dishes for directional communication and log-periodic antennas for wideband reception, physically isolated or shielded from interference. Receivers consist of software-defined radios (SDRs) or specialized RF receivers, often GNU Radio-based for flexibility. The Command Encoder formats and encodes uplink commands according to CCSDS standards. The Telemetry Processor decodes, decompresses, and validates downlink data. The Mission Control System serves as the primary operator interface, running on hardened, isolated networks. The Data Archive stores telemetry historical data for analysis and anomaly detection.
Security Best Practices:
Comprehensive ground station security requires implementing network segmentation with firewalls between critical and non-critical systems, preventing lateral movement of compromised components. Authentication and authorization must employ multi-factor approaches to verify operator identity and authorization level. All uplink commands require encryption with cryptographic validation to prevent command injection attacks. Intrusion detection systems (IDS) continuously monitor all network traffic for suspicious patterns and anomalies. Regular security audits and penetration testing by independent security teams identify vulnerabilities before attackers exploit them. Backup communication paths provide operational continuity if the primary system becomes compromised. Personnel security practices including comprehensive background checks for critical roles prevent insider threats and reduce human-factor vulnerability.
Advanced Ground Station Threat Model & Mitigation Strategies:
Threat actors targeting ground stations use multiple vectors simultaneously. The U.S. Space Force identified 1,200+ cyberattacks on space-related infrastructure in 2024 alone.
Attack Vector 1: Supply Chain Compromise of RF Equipment - Adversaries infiltrate manufacturers to insert firmware backdoors in receivers or software like GNU Radio. Detection: Digital signatures verification, firmware integrity checks with TPM (Trusted Platform Module), air-gapped validation systems. Mitigation: Manufacturer certificates with hardware security modules (HSM), continuous firmware attestation proofs.
Attack Vector 2: Signal Hijacking via Rogue Ground Station - Attacker establishes unauthorized ground station on satellite's uplink frequency, transmitting fraudulent telecommands. Detection: Geolocation analysis (satellite antenna pointing must match permitted ground stations), command authentication with public key infrastructure (PKI). Mitigation: Cryptographic command signing, frequency hopping protocols, multi-ground-station consensus for critical operations.
Attack Vector 3: Network Lateral Movement - Compromise of non-critical system (e.g., weather station networked with MCS) leads to privilege escalation. Detection: Behavioral anomaly detection on command injection attempts, deep packet inspection (DPI) on all outbound frames. Mitigation: Microsegmentation with zero-trust architecture, application-layer firewalls, memory tagging for exploit prevention.
import time
import logging
from enum import Enum
from dataclasses import dataclass
class ThreatLevel(Enum):
NORMAL = 0
ELEVATED = 1
CRITICAL = 2
@dataclass
class SecurityEvent:
timestamp: float
severity: ThreatLevel
event_type: str
description: str
action_taken: str
class AdvancedGroundStationDefense:
"""
Zero-Trust Ground Station Defense System
Applies Defense in Depth with continuous verification
"""
def __init__(self):
self.security_events = []
self.threat_level = ThreatLevel.NORMAL
self.command_queue = []
self.geofence_stations = {
'Primary_GSS_USA': {'lat': 35.2, 'lon': -106.5, 'radius_km': 5},
'Backup_GSS_AUS': {'lat': -35.4, 'lon': 149.2, 'radius_km': 5}
}
self.last_valid_command_time = time.time()
def geographic_authentication(self, signal_source_lat, signal_source_lon):
"""
Verify signal origin matches authorized ground station geofence
Uses TDOA (Time Difference of Arrival) cross-verification
"""
for station_name, geofence in self.geofence_stations.items():
distance_km = self.haversine_distance(
signal_source_lat, signal_source_lon,
geofence['lat'], geofence['lon']
)
if distance_km < geofence['radius_km']:
return True, station_name
return False, None
def haversine_distance(self, lat1, lon1, lat2, lon2):
"""Calculate great-circle distance between two points on Earth"""
import math
R = 6371 # Earth radius in km
dlat = math.radians(lat2 - lat1)
dlon = math.radians(lon2 - lon1)
a = (math.sin(dlat/2)**2 +
math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *
math.sin(dlon/2)**2)
c = 2 * math.asin(math.sqrt(a))
return R * c
def command_signing_verification(self, tc_frame, origin_station):
"""
Multi-signature verification for telecommand
Requires: Origin station signature + Supervisor signature + Director signature
"""
required_signatures = {
'operator': False,
'supervisor': False,
'director': False
}
# Extract signatures from TC frame
extracted_sigs = self.extract_signatures(tc_frame)
# Verify each signature against public keys
operator_pub_key = self.get_station_public_key(origin_station, 'operator')
if self.verify_rsa_signature(extracted_sigs['operator'], operator_pub_key):
required_signatures['operator'] = True
supervisor_pub_key = self.get_station_public_key(origin_station, 'supervisor')
if self.verify_rsa_signature(extracted_sigs['supervisor'], supervisor_pub_key):
required_signatures['supervisor'] = True
director_pub_key = self.get_station_public_key(origin_station, 'director')
if self.verify_rsa_signature(extracted_sigs['director'], director_pub_key):
required_signatures['director'] = True
# All three required
all_signed = all(required_signatures.values())
if not all_signed:
self.log_security_event(
ThreatLevel.ELEVATED,
'SIGNATURE_VERIFICATION_FAILURE',
f'Missing signatures: {required_signatures}'
)
return all_signed
def anomaly_detection_ensemble(self, telemetry):
"""
Ensemble anomaly detection combining multiple ML models
- Isolation Forest for multivariate anomalies
- LSTM for time-series anomalies
- Statistical methods (Z-score, Mahalanobis)
"""
anomaly_scores = []
# Statistical anomaly
z_scores = self.calculate_z_scores(telemetry)
stat_anomaly = max(z_scores) > 3.5 # p < 0.0001
anomaly_scores.append(0.95 if stat_anomaly else 0.05)
# Isolation Forest anomaly
# (Requires training on normal telemetry baseline)
iso_forest_score = self.isolation_forest_score(telemetry)
anomaly_scores.append(iso_forest_score)
# LSTM time-series anomaly
# (Detects departure from learned temporal patterns)
lstm_score = self.lstm_reconstruction_error(telemetry)
anomaly_scores.append(lstm_score)
# Ensemble voting
ensemble_anomaly_score = sum(anomaly_scores) / len(anomaly_scores)
if ensemble_anomaly_score > 0.75:
self.threat_level = ThreatLevel.CRITICAL
self.log_security_event(
ThreatLevel.CRITICAL,
'ENSEMBLE_ANOMALY_DETECTED',
f'Ensemble score: {ensemble_anomaly_score:.3f}'
)
return True
return False
def log_security_event(self, severity, event_type, description):
"""Log security event with timestamp and recommended action"""
event = SecurityEvent(
timestamp=time.time(),
severity=severity,
event_type=event_type,
description=description,
action_taken=self.determine_action(severity)
)
self.security_events.append(event)
# Real-time alerting for critical events
if severity == ThreatLevel.CRITICAL:
self.trigger_critical_alert(event)
def determine_action(self, severity):
"""Determine automated response based on threat level"""
if severity == ThreatLevel.NORMAL:
return 'MONITOR'
elif severity == ThreatLevel.ELEVATED:
return 'ISOLATE_SATELLITE'
else: # CRITICAL
return 'EMERGENCY_DISCONNECT'
def trigger_critical_alert(self, event):
"""Send alerts to all authorized personnel immediately"""
print(f"[CRITICAL ALERT] {event.event_type}: {event.description}")
# In production: send SMS, email, activate emergency bridges
# Example defense system
defense = AdvancedGroundStationDefense()
print("Advanced Ground Station Defense System initialized")
print("- Geographic authentication enabled")
print("- Multi-signature command verification active")
print("- Ensemble anomaly detection operational")
Comprehensive analysis of GNSS systems (GPS, GLONASS, Galileo, BeiDou) vulnerabilities to spoofing and jamming attacks. Covers civilian vs. military signal spoofing methodologies, multi-constellation attack coordination, and sophisticated detection algorithms. Implements pseudorange residual analysis, signal quality (SNR) monitoring, constellation consistency verification, and Mahalanobis distance-based anomaly detection for identifying spoofing attempts.
Global Navigation Satellite Systems (GNSS) Overview:
| System | Country | Satellites | Frequency | Accuracy |
|---|---|---|---|---|
| GPS | USA | 31+ | 1.2-1.6 GHz | 1-5 m (L1) |
| GLONASS | Russia | 24+ | 1.6-1.7 GHz | 5-10 m |
| Galileo | EU | 30 | 1.2-1.6 GHz | 1-4 m |
| BeiDou | China | 35+ | 1.2-1.6 GHz | 5-10 m |
Spoofing Signatures and Detection Methods:
GPS spoofing involves transmitting false satellite signals to deceive receivers. Detection methodologies include multi-antenna monitoring to detect directional inconsistencies that reveal non-satellite signal origins, signal quality monitoring to identify SNR anomalies characteristic of synthetic signals, cross-checking with alternative navigation sources (such as inertial navigation systems) to identify discrepancies that indicate spoofing, cryptographic authentication on military signals (M-Code, Y-Code) to authenticate signal origins with tamper-evident properties, and continuous monitoring of constellation geometry to detect impossible or improbable satellite positions that spoofing systems cannot replicate.
Jamming Mitigation Strategies:
Protecting navigation systems against jamming requires multi-layered technical approaches. Frequency hopping and spread spectrum techniques randomize signal frequency to prevent direct jamming. Directional antennas and null steering techniques physically reject jamming signals arriving from non-satellite directions. Spectral monitoring for anomalous signals enables detection of active jamming attempts before service degradation. Backup navigation using alternative systems provides operational continuity when GPS becomes unavailable. Anti-jamming techniques in receiver design include adaptive filtering, notch filters for specific interfering signals, and signal processing algorithms optimized for detecting weak satellite signals in high-interference environments.
Advanced GPS Spoofing Attack Methodologies:
GPS spoofing attacks have evolved significantly. In 2022, researchers demonstrated "meaconing" attacks where captured authentic GPS signals are replayed with modified timestamps, creating a false location fix up to 10 kilometers away. The attack complexity requires modest RF equipment ($5,000-$50,000) and sophisticated signal processing.
Attack Classification:
Type 1 - Civilian GPS Spoofing: Targets unauthenticated L1 (1575.42 MHz) civilian signal. No authentication, making it the easiest attack vector. Attackers generate phantom satellites in a false constellation to deceive receivers.
Type 2 - Military/M-Code Spoofing: Targets encrypted military signals (M-Code on GPS blocks III/IV). Requires prior knowledge of encryption keys. Less common but much higher impact. Government facilities and military assets remain vulnerable during key transition periods.
Type 3 - Multi-Constellation Spoofing: Simultaneously spoof GPS, GLONASS, and Galileo to overcome consistency checks. Requires coordinated RF transmission on multiple frequency bands (1.2-1.6 GHz). Detection becomes extremely difficult if spoofed constellation geometries are geometrically consistent.
import numpy as np
from scipy.signal import correlate, convolve
from scipy.optimize import least_squares
class GPSSpoofingAttackSimulator:
"""
Simulates and detects GPS spoofing attacks
Implements multiple detection algorithms
"""
def __init__(self):
# GPS satellite constellation (simplified)
self.satellites = {
'PRN_1': {'x': 20200000, 'y': 10000000, 'z': 0},
'PRN_2': {'x': 10200000, 'y': -20000000, 'z': 0},
'PRN_3': {'x': -10200000, 'y': 10000000, 'z': 20000000},
'PRN_4': {'x': 0, 'y': -10200000, 'z': -20000000},
}
self.receiver_true_position = np.array([6378137, 0, 0]) # On Earth surface
self.c = 3e8 # Speed of light (m/s)
def generate_truth_pseudoranges(self):
"""Calculate authentic pseudoranges from real satellites"""
pseudoranges = {}
for prn, sat_pos in self.satellites.items():
sat_array = np.array([sat_pos['x'], sat_pos['y'], sat_pos['z']])
range_m = np.linalg.norm(sat_array - self.receiver_true_position)
pseudoranges[prn] = range_m
return pseudoranges
def spoofing_attack_civilian(self, false_position, range_error=0):
"""
Generate spoofed civil GPS signals
false_position: [x, y, z] where attacker wants receiver to believe it is
range_error: Additional range bias (meters) to make spoofed sats appear closer
"""
spoofed_pseudoranges = {}
for prn, sat_pos in self.satellites.items():
sat_array = np.array([sat_pos['x'], sat_pos['y'], sat_pos['z']])
# Calculate range to false position
false_range = np.linalg.norm(sat_array - false_position)
# Add range error to make satellites appear stronger
spoofed_range = false_range + range_error
spoofed_pseudoranges[prn] = spoofed_range
return spoofed_pseudoranges
def detect_spoofing_residual_analysis(self, measured_pseudoranges):
"""
Detect spoofing via pseudorange residual analysis
Compares measured vs. predicted ranges based on known constellation ephemeris
"""
# Use least squares to solve: position + bias
# Measured = sqrt((x-x_sat)^2 + (y-y_sat)^2 + (z-z_sat)^2) + clock_bias
initial_guess = np.array([6378137, 0, 0, 0]) # [x, y, z, bias]
def residuals(state):
x, y, z, bias = state
receiver_pos = np.array([x, y, z])
res = []
for prn, measured_range in measured_pseudoranges.items():
sat_pos = np.array([
self.satellites[prn]['x'],
self.satellites[prn]['y'],
self.satellites[prn]['z']
])
predicted_range = np.linalg.norm(sat_pos - receiver_pos) + bias
residual = measured_range - predicted_range
res.append(residual)
return np.array(res)
# Solve for position
result = least_squares(residuals, initial_guess)
solution_position = result.x[:3]
final_residuals = result.fun
# Spoofing detection: if residuals are unusually large or systematic
residual_std = np.std(final_residuals)
residual_max = np.max(np.abs(final_residuals))
# Threshold: typical GPS residuals < 10 meters with 4+ satellites
# Spoofed signals often have larger systematic errors
spoofing_detected = residual_max > 50 or residual_std > 30
return {
'estimated_position': solution_position,
'residual_std': residual_std,
'residual_max': residual_max,
'spoofing_detected': spoofing_detected,
'confidence': min(1.0, residual_std / 20)
}
def detect_spoofing_signal_quality(self, received_signals):
"""
Detect spoofing via signal quality (SNR) anomalies
Authentic signals have realistic SNR degradation based on geometry
Spoofed signals often have artificially uniform, strong SNRs
"""
snr_values = [sig['snr_db'] for sig in received_signals]
# Analysis
snr_mean = np.mean(snr_values)
snr_std = np.std(snr_values)
snr_min = np.min(snr_values)
# Red flags for spoofing:
# 1. Unreasonably high SNR (> 50 dB is suspicious in civilian GPS)
# 2. Very low variation in SNR (all signals equally strong)
# 3. All satellites above geometric horizon but reported SNR very uniform
suspicious_factors = 0
if snr_mean > 48: # Civilian GPS typically < 45 dB
suspicious_factors += 1
if snr_std < 2: # Authentic has 3-5 dB variation typical
suspicious_factors += 1
if snr_min > 45: # Some satellites should be weaker
suspicious_factors += 1
spoofing_likelihood = suspicious_factors / 3
return {
'snr_mean': snr_mean,
'snr_std': snr_std,
'suspicious_factors': suspicious_factors,
'spoofing_likelihood': spoofing_likelihood
}
def detect_spoofing_constellation_consistency(self, ephemeris_database,
measured_positions):
"""
Detect spoofing by checking if reported satellite constellation
matches known ephemeris and orbital mechanics
"""
# Check: Do satellite positions violate known orbital laws?
# - Keplerian orbit constraints
# - Continental visibility zones
# - Orbital velocity consistency
anomalies = 0
for sat_id, measured_pos in measured_positions.items():
# Should exist in ephemeris
if sat_id not in ephemeris_database:
anomalies += 1
continue
true_pos = ephemeris_database[sat_id]
distance_error = np.linalg.norm(measured_pos - true_pos)
# Typical ephemeris uncertainty: 10-50 meters
# > 100 meters suggests spoofing
if distance_error > 100:
anomalies += 1
return {
'ephemeris_anomalies': anomalies,
'total_satellites': len(measured_positions),
'spoofing_confidence': min(1.0, anomalies / len(measured_positions))
}
# Example attack
simulator = GPSSpoofingAttackSimulator()
# Generate authentic pseudoranges
authentic_ranges = simulator.generate_truth_pseudoranges()
print(f"Authentic pseudoranges: {authentic_ranges}")
# Attacker generates spoofed signals with false position 10 km away
false_position = np.array([6378137 + 10000, 0, 0])
spoofed_ranges = simulator.spoofing_attack_civilian(false_position, range_error=100)
# Analyze with detection algorithm
detection_result = simulator.detect_spoofing_residual_analysis(spoofed_ranges)
print(f"\\nSpoofing Detection Result:")
print(f"Estimated position: {detection_result['estimated_position']}")
print(f"Spoofing Detected: {detection_result['spoofing_detected']}")
print(f"Confidence: {detection_result['confidence']:.3f}")
Strategic defense framework applying NIST Cybersecurity Model (Identify-Protect-Detect-Respond-Recover) and Zero Trust Architecture to satellite systems. Covers hardware security measures (HSMs, secure boot, tamper detection), cryptographic key management, redundancy strategies, and incident response procedures. Emphasizes continuous monitoring, immutable audit trails, and resilience against sophisticated, multi-vector attacks.
NIST Cybersecurity Framework for Space Systems:
The NIST Cybersecurity Framework provides a comprehensive approach to protecting space systems through five core functions. The Identify phase establishes foundational understanding by cataloging all system assets, thoroughly mapping system architecture, and precisely identifying critical functions that must be protected. The Protect phase implements defensive measures through access controls, comprehensive encryption schemes, multi-factor authentication mechanisms, and secure system design principles. The Detect phase establishes persistent monitoring capabilities by deploying continuous monitoring systems, implementing anomaly detection algorithms, and deploying intrusion detection systems to identify unauthorized access attempts. The Respond phase establishes operational resilience procedures including documented incident response procedures, system isolation protocols, and comprehensive forensic capability to preserve evidence and analyze attacks. The Recover phase ensures operational continuity through redundant backup systems, alternative communication paths, and automated failover mechanisms that minimize downtime during security incidents.
Zero Trust Architecture for Space:
Zero Trust Architecture represents a fundamental shift in security philosophy by rejecting the traditional perimeter-based defense model. Rather than assuming the internal network is inherently secure, Zero Trust Architecture assumes all components are potentially compromised and requires continuous validation. Every communication between systems must be authenticated and authorized, not just at the network perimeter but at every point within the system. All data in transit between satellites, ground stations, and control centers must be encrypted using cryptographically robust algorithms, and all data at rest in storage systems must be similarly protected. The architecture implements continuous verification and monitoring throughout the entire operational lifetime, continuously reassessing trust based on current threat assessments and system behavior. The principle of least privilege is rigorously applied to all access requests, ensuring that each user, process, and system component receives only the minimum permissions necessary to perform its designated function. All critical operations are recorded in immutable audit logs that provide a complete chain of custody and forensic capability, enabling comprehensive investigation of any security incidents.
Hardware Security Measures:
Protection of satellite systems requires comprehensive hardware-level security mechanisms that operate independently of software security systems. Hardware security modules (HSMs) provide dedicated physical devices for cryptographic key storage, ensuring that critical keys are protected against software-based attacks through tamper-resistant hardware designs. Secure boot mechanisms and verified firmware loading procedures ensure that satellites execute only authenticated firmware that has been cryptographically signed by authorized operators, preventing the installation and execution of malicious firmware. Physical tamper detection systems deployed on satellites and ground equipment infrastructure provide immediate alerts when unauthorized physical access or modification is attempted, enabling rapid detection of sabotage or espionage attempts. Radiation-hardened processors used throughout spacecraft are specifically designed to prevent single-event upsets (SEUs) caused by cosmic radiation and solar particle events, which could otherwise corrupt critical software or control systems and compromise spacecraft integrity in the harsh space environment.
Advanced techniques for satellite propagation, conjunction assessment, and collision probability calculation. Covers SSA fundamentals, space debris tracking using Mahalanobis distance metrics, and probability-of-collision (Pc) thresholds. Implements high-precision numerical integration with full perturbation modeling, covariance propagation for uncertainty quantification, and automated maneuver decision logic for satellite safety operations.
Space Situational Awareness (SSA): Space Situational Awareness represents the fundamental capability to detect, track, and precisely characterize all objects and events occurring in space. This capability is essential for multiple critical applications across military, commercial, and scientific domains. Collision avoidance depends on SSA to predict encounters between active satellites and space debris with sufficient accuracy to compute and execute evasive maneuvers before potential impact. Comprehensive debris tracking provides detailed catalogs of known objects in orbit, enabling operators to maintain operational awareness of the debris environment and make informed decisions about spacecraft safety. Threat detection capabilities identify potential hostile actions, spacecraft malfunctions, or anomalous activities that could compromise national security or satellite mission success. Launch planning requires accurate SSA data to identify available orbital slots that avoid known debris and active satellites, preventing tragic collisions during orbital insertion. National security intelligence applications rely on SSA to monitor foreign military and civilian space activities, providing critical early warning of potential threats to national infrastructure and space-based military systems.
Orbital Mechanics Fundamentals:
Kepler's laws describe ideal two-body motion. Modern propagators must account for all perturbations mentioned previously. State-of-the-art systems achieve position accuracy of ±100 meters or better for well-tracked objects using statistical methods and multi-sensor fusion.
Data Sources for Tracking:
Modern space surveillance leverages multiple complementary sensor types to establish comprehensive orbital awareness. Ground-based radar systems operating in the S-band (2-4 GHz) and X-band (8-12 GHz) provide continuous coverage of the space domain, tracking objects by measuring range, range-rate (doppler), and angular position with precision. Optical observation systems including ground-based telescopes and laser ranging systems provide complementary data, particularly effective for faint objects and high-altitude satellites where radar is less effective, offering excellent accuracy in angular measurements and precise distance measurements through laser ranging. Space-based sensors including the Space Force Space Surveillance System provide continuous monitoring across the entire globe without the weather dependencies or tracking gaps inherent in ground-based systems, enabling comprehensive surveillance of geostationary and high-altitude objects. Two-Line Elements (TLEs) released publicly by NORAD provide freely available orbital elements for thousands of tracked objects, serving as a convenient reference source for non-classified tracking and mission planning by commercial operators and researchers. Satellite signal reception and analysis enable independent confirmation of orbital parameters and object identity through direct communication with active satellites, providing verification of orbital state propagation accuracy and detection of unexpected orbital maneuvers or anomalies.
Conjunction Assessment: Predicting potential collisions between satellites and debris represents a critical safety function in modern space operations. Conjunction assessment algorithms integrate current orbital state information, covariance matrices representing position uncertainty, and debris catalog data to compute miss distances and collision probabilities. Common practical thresholds establish miss distances under 1 kilometer as representing significant collision risk in LEO environments requiring immediate evasive action. However, more conservative thresholds of 200 meters are commonly employed for high-value assets like the International Space Station. Modern automated systems perform continuous conjunction assessment, automatically flagging potential close approaches and computing optimal maneuver solutions to avoid collision while minimizing fuel expenditure and preserving orbital mission objectives.
Detailed exploitation methodologies for disrupting satellite command and telemetry using advanced signal jamming and GPS spoofing techniques. Covers chirp jamming, adaptive multi-tone jamming with Jamming-to-Signal Ratio (JSR) calculations, and GPS spoofing detection via SNR monitoring, cryptographic authentication, INS cross-checking, and multi-antenna direction-of-arrival analysis. Includes complete Python implementations for both attack and detection.
Satellite Signal Detection Script:
#!/usr/bin/env python3
"""
Satellite Signal Detection Script
Detects and analyzes satellite signals in S-band
Requires: GNU Radio, HackRF or USRP
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from gnuradio import gr, blocks, analog, filter
import pmt
class SatelliteSignalDetector(gr.top_block):
def __init__(self, freq=2.4e9, samp_rate=2e6):
gr.top_block.__init__(self, "Satellite Signal Detector")
# Configure SDR source
self.src = osmosdr.source()
self.src.set_sample_rate(samp_rate)
self.src.set_center_freq(freq)
self.src.set_gain_mode(False)
self.src.set_gain(40)
# Bandpass filter for target bandwidth
taps = filter.firdes.band_pass(1, samp_rate, freq-1e6, freq+1e6, 100e3)
self.bandpass = filter.fir_filter_ccc(1, taps)
# Energy detector for signal presence
self.mag = blocks.complex_to_mag()
self.threshold = blocks.threshold_ff(0.1, 0.15, 0)
# File sink for captured signals
self.sink = blocks.file_sink(gr.sizeof_gr_complex, "captured_signal.dat")
# Connect blocks
self.connect(self.src, self.bandpass)
self.connect(self.bandpass, self.mag)
self.connect(self.mag, self.threshold)
self.connect(self.bandpass, self.sink)
def detect_signal(self):
"""Run detection for specified time"""
print(f"Scanning for satellite signals at {self.src.get_center_freq()/1e9:.2f} GHz")
self.start()
self.wait()
# Analyze captured signal
data = np.fromfile("captured_signal.dat", dtype=np.complex64)
if len(data) > 0:
print(f"Signal detected: {len(data)} samples")
print(f"Signal power: {10*np.log10(np.mean(np.abs(data)**2)):.2f} dBm")
return True
return False
def estimate_bandwidth(self, signal_data, samp_rate=2e6):
"""Estimate signal bandwidth using FFT"""
fft_data = np.fft.fft(signal_data)
freqs = np.fft.fftfreq(len(signal_data), 1/samp_rate)
# Find -3dB points
psd = np.abs(fft_data)**2
peak_idx = np.argmax(psd)
peak_power = psd[peak_idx]
# Find bandwidth at half power
half_power = peak_power / 2
mask = psd >= half_power
bandwidth = np.max(freqs[mask]) - np.min(freqs[mask])
return bandwidth
# Scan common satellite frequencies
frequencies = [1.57542e9, # GPS L1
2.2e9, # S-band mobile
2.4e9, # ISM/amateur satellite
2.6e9] # Iridium
for freq in frequencies:
detector = SatelliteSignalDetector(freq=freq)
detected = detector.detect_signal()
if detected:
print(f"Signal found at {freq/1e9:.2f} GHz")
Advanced Jamming Signal Generation:
import numpy as np
from scipy import signal as sp
import matplotlib.pyplot as plt
class AdvancedJammer:
"""Generate sophisticated jamming signals"""
def __init__(self, center_freq, bandwidth, power_dbm):
self.center_freq = center_freq
self.bandwidth = bandwidth
self.power_dbm = power_dbm
def generate_chirp_jamming(self, duration, chirp_type='linear'):
"""
Generate chirp jamming signal
Types: linear, logarithmic, hyperbolic
"""
t = np.linspace(0, duration, int(duration * 100e6))
if chirp_type == 'linear':
chirp_signal = sp.chirp(t, self.center_freq - self.bandwidth/2,
duration, self.center_freq + self.bandwidth/2)
elif chirp_type == 'logarithmic':
chirp_signal = sp.chirp(t, self.center_freq - self.bandwidth/2,
duration, self.center_freq + self.bandwidth/2, method='logarithmic')
# Add noise for more effectiveness
noise = np.random.normal(0, 0.1, len(t))
jamming_signal = chirp_signal + noise
# Normalize to desired power
current_power = 10 * np.log10(np.mean(jamming_signal**2))
gain = self.power_dbm - current_power
jamming_signal *= 10**(gain/20)
return jamming_signal
def generate_adaptive_jamming(self, target_signal):
"""
Adaptive jamming that analyzes and targets specific signal characteristics
"""
# Analyze target signal
fft_signal = np.fft.fft(target_signal)
freqs = np.fft.fftfreq(len(target_signal))
# Find dominant frequencies
psd = np.abs(fft_signal)**2
dominant_freqs = freqs[np.argsort(psd)[-5:]]
# Create multi-tone jamming at dominant frequencies
t = np.linspace(0, 1, len(target_signal))
jamming_signal = np.zeros_like(target_signal)
for freq in dominant_freqs:
jamming_signal += np.sin(2 * np.pi * freq * t)
return jamming_signal
def calculate_jamming_effectiveness(self, signal_power, distance_km):
"""
Calculate Jamming-to-Signal Ratio (JSR)
Based on free-space path loss
"""
# Free-space path loss
c = 3e8 # Speed of light
wavelength = c / self.center_freq
fsp1 = 20 * np.log10(distance_km * 1000) + 20 * np.log10(self.center_freq) + 20 * np.log10(4 * np.pi / c)
# Received jamming power
p_j = self.power_dbm - fsp1
# Jamming-to-Signal Ratio
jsr = p_j - signal_power
return jsr
# Example usage
jammer = AdvancedJammer(center_freq=1.57542e9, # GPS L1
bandwidth=10e6, # 10 MHz bandwidth
power_dbm=30) # 30 dBm (1W) transmitter
# Generate different jamming signals
chirp_jam = jammer.generate_chirp_jamming(duration=0.001, chirp_type='logarithmic')
# Calculate effectiveness at 10km distance
jsr = jammer.calculate_jamming_effectiveness(signal_power=-130, distance_km=10)
print(f"Jamming-to-Signal Ratio at 10km: {jsr:.2f} dB")
if jsr > 20:
print("Effective jamming possible")
else:
print("Insufficient jamming power")
GPS Spoofing Detection:
import numpy as np
from scipy import stats
class GPSSpoofingDetector:
"""
Detect GPS spoofing using multiple techniques:
1. Signal Quality Monitoring
2. Cryptographic Authentication (SAASM)
3. Multi-antenna Analysis
4. Inertial Navigation System (INS) Cross-check
"""
def __init__(self):
self.position_history = []
self.velocity_history = []
self.snr_history = []
def monitor_signal_quality(self, current_snr, expected_snr=45):
"""
Monitor Signal-to-Noise Ratio for anomalies
Spoofed signals often have unusually consistent or high SNR
"""
self.snr_history.append(current_snr)
if len(self.snr_history) > 10:
snr_std = np.std(self.snr_history[-10:])
snr_mean = np.mean(self.snr_history[-10:])
indicators = []
# Suspicious: SNR too consistent (spoofed signals)
if snr_std < 0.5:
indicators.append("Suspiciously consistent SNR")
# Suspicious: SNR higher than expected
if snr_mean > expected_snr:
indicators.append("SNR higher than expected")
return indicators
return []
def cryptographic_authentication(self, navigation_message, key=None):
"""Verify GPS cryptographic authentication (SAASM/M-Code)"""
# Check message structure
if len(navigation_message) != 300:
return "INVALID: Wrong message length"
# Check subframe IDs (1-5)
subframe_id = (navigation_message[50] >> 2) & 0x7
if subframe_id not in [1, 2, 3, 4, 5]:
return "INVALID: Bad subframe ID"
return "AUTHENTIC"
def inertial_navigation_cross_check(self, gps_position, gps_velocity,
ins_position, ins_velocity, time_delta):
"""
Compare GPS with Inertial Navigation System
INS accumulates error over time, GPS doesn't
"""
# Calculate expected INS drift
expected_ins_error = 1.5 * time_delta # meters
# Calculate position difference
pos_diff = np.linalg.norm(np.array(gps_position) - np.array(ins_position))
vel_diff = np.linalg.norm(np.array(gps_velocity) - np.array(ins_velocity))
indicators = []
if pos_diff < expected_ins_error:
pass # Normal
elif pos_diff > 50: # Large sudden jump
indicators.append("Large position jump detected")
# Check for impossible maneuvers
if vel_diff > 100: # More than 100 m/s velocity difference
indicators.append("Impossible velocity change detected")
return indicators
def multi_antenna_analysis(self, antenna_positions, measured_phases):
"""
Use multiple antennas to detect direction of arrival
Spoofed signals come from single direction (ground)
Real satellite signals come from multiple directions
"""
# Convert phase differences to direction vectors
directions = []
for i in range(len(antenna_positions)-1):
phase_diff = measured_phases[i+1] - measured_phases[i]
direction = np.angle(phase_diff)
directions.append(direction)
# Check if all signals come from similar direction
if len(directions) > 3:
direction_std = np.std(directions)
if direction_std < 0.1:
return ["Signals from single direction - possible spoofing"]
return []
# Example usage
detector = GPSSpoofingDetector()
for i in range(100):
# Simulated GPS data
snr = 45 + np.random.normal(0, 1) # Normal variation
# Simulate spoofing attack at time 50
if i >= 50:
snr = 50 # Spoofer gives better signal
indicators = detector.monitor_signal_quality(snr)
if indicators:
print(f"Time {i}: Possible spoofing detected - {indicators}")
Systematic analysis of satellite subsystem vulnerabilities: attitude determination and control (ADCS) systems vulnerable to laser-based star tracker spoofing, reaction wheel saturation attacks inducing uncontrolled tumbling, and power system exploitation via eclipse scenarios and overcharge attacks. Includes simulations of power-based anomaly detection (SOC, temperature, load monitoring) and mathematical models of spacecraft dynamics and control system responses.
Star Tracker Spoofing Simulation:
import numpy as np
from scipy.spatial.transform import Rotation
class StarTrackerSpoofer:
"""
Simulate star tracker spoofing attack
Creates false star patterns to manipulate attitude determination
"""
def __init__(self, catalog_size=1000):
# Generate star catalog
self.star_catalog = self.generate_star_catalog(catalog_size)
def generate_star_catalog(self, n_stars):
"""Generate synthetic star positions (RA, Dec)"""
# Uniform distribution on celestial sphere
ra = np.random.uniform(0, 2*np.pi, n_stars)
dec = np.arcsin(np.random.uniform(-1, 1, n_stars))
# Assign random magnitudes
magnitude = np.random.uniform(1, 6, n_stars)
return list(zip(ra, dec, magnitude))
def spoof_star_pattern(self, actual_attitude, desired_spoof_attitude, fov=20):
"""Generate spoofed star pattern"""
# Convert attitudes to rotation matrices
R_actual = Rotation.from_euler('xyz', actual_attitude).as_matrix()
R_spoof = Rotation.from_euler('xyz', desired_spoof_attitude).as_matrix()
# Calculate rotation difference
R_diff = R_spoof @ R_actual.T
# Spoof stars in FOV
fov_rad = np.radians(fov/2)
spoofed_stars = []
for star in self.star_catalog:
if len(spoofed_stars) < 10:
spoofed_stars.append(star)
return spoofed_stars
def calculate_spoofing_power(self, target_satellite, ground_station):
"""Calculate required laser power for star tracker spoofing"""
# Parameters
satellite_altitude = 500e3 # meters
aperture_diameter = 0.1 # Star tracker aperture in meters
pixel_size = 10e-6 # Pixel size in meters
integration_time = 0.1 # Integration time in seconds
quantum_efficiency = 0.6
# Calculate spot size at satellite
wavelength = 550e-9 # Green laser
beam_divergence = 1.22 * wavelength / aperture_diameter
spot_diameter = beam_divergence * satellite_altitude
# Required photons per pixel
required_photons = 100
# Calculate required laser power
energy_per_photon = 6.626e-34 * 3e8 / wavelength
required_energy = required_photons * energy_per_photon / quantum_efficiency
required_power = required_energy / integration_time
# Account for atmospheric transmission and beam spread
atmospheric_loss = 0.8
area_ratio = (spot_diameter**2) / (pixel_size**2)
total_power = required_power / atmospheric_loss * area_ratio
return total_power
# Example usage
spoofer = StarTrackerSpoofer()
# Actual vs desired attitude
actual_attitude = [0.1, 0.05, 0.2]
desired_spoof_attitude = [0.5, 0.3, 0.8]
spoofed_stars = spoofer.spoof_star_pattern(actual_attitude, desired_spoof_attitude)
print(f"Generated {len(spoofed_stars)} spoofed stars")
print(f"Spoofing attempts to change attitude by: {np.degrees(np.array(desired_spoof_attitude) - np.array(actual_attitude))} degrees")
required_power = spoofer.calculate_spoofing_power(500, 0)
print(f"Required laser power for spoofing: {required_power:.2f} W")
Reaction Wheel Saturation Attack:
import numpy as np
from scipy.integrate import odeint
class ReactionWheelController:
"""Simulate reaction wheel control system with saturation attack"""
def __init__(self, n_wheels=4):
self.n_wheels = n_wheels
self.wheel_speeds = np.zeros(n_wheels)
self.max_speed = 6000 # RPM
self.max_torque = 0.1 # Nm
self.momentum_capacity = 50 # Nms
# Wheel configuration matrix (pyramid)
self.config_matrix = np.array([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1],
[0, 0, 1]
]).T
def apply_torque(self, desired_torque, dt=0.1):
"""Apply torque to reaction wheels"""
# Calculate wheel torques
wheel_torques = np.linalg.pinv(self.config_matrix) @ desired_torque
# Check for saturation
saturated = False
for i in range(self.n_wheels):
if abs(wheel_torques[i]) > self.max_torque:
saturated = True
wheel_torques[i] = np.sign(wheel_torques[i]) * self.max_torque
# Update wheel speeds
self.wheel_speeds += wheel_torques * dt
# Calculate actual torque delivered
actual_torque = self.config_matrix @ wheel_torques
return actual_torque, saturated
def saturation_attack(self, attack_duration=30, dt=0.1):
"""Simulate saturation attack"""
times = np.arange(0, attack_duration, dt)
saturation_events = []
# Attack pattern: oscillating torque
resonance_freq = 10 # Hz
for t in times:
# Apply torque designed to saturate wheels
attack_torque = 0.15 * np.array([
np.sin(2*np.pi*resonance_freq*t),
np.cos(2*np.pi*resonance_freq*t),
0
])
actual_torque, saturated = self.apply_torque(attack_torque, dt)
if saturated:
saturation_events.append(t)
return saturation_events
# Example usage
controller = ReactionWheelController()
print("Starting saturation attack...")
attack_events = controller.saturation_attack(attack_duration=10)
if attack_events:
print(f"Attack caused {len(attack_events)} saturation events")
print(f"Final wheel speeds: {controller.wheel_speeds}")
total_momentum = np.sum(np.abs(controller.wheel_speeds))
saturation_percentage = (total_momentum / controller.momentum_capacity) * 100
print(f"Momentum saturation: {saturation_percentage:.1f}%")
Power System Attack Detection:
import numpy as np
from dataclasses import dataclass
from typing import List
@dataclass
class PowerSystemState:
"""Current state of satellite power system"""
battery_charge: float # Wh
solar_power: float # W
load_power: float # W
battery_temp: float # °C
eclipse: bool
time: float # hours
class PowerSystemSimulator:
"""Simulate satellite power system and detect attacks"""
def __init__(self):
# System parameters
self.battery_capacity = 2000 # Wh
self.max_charge_rate = 500 # W
self.max_discharge_rate = 1000 # W
self.solar_panel_area = 10 # m²
self.panel_efficiency = 0.3
# Normal operating ranges
self.normal_battery_range = (0.2, 0.95) # 20-95% charge
self.normal_temp_range = (-10, 30) # °C
self.critical_loads = ['comms', 'computer', 'thermal']
self.state = PowerSystemState(
battery_charge=1000,
solar_power=500,
load_power=200,
battery_temp=15,
eclipse=False,
time=0
)
self.history = []
def update_state(self, dt=0.1, commands=None):
"""Update power system state"""
if commands:
self.state.load_power = commands.get('load', self.state.load_power)
# Calculate power balance
net_power = self.state.solar_power - self.state.load_power
if self.state.eclipse:
self.state.solar_power = 0
net_power = -self.state.load_power
# Update battery
if net_power > 0:
charge_increment = min(net_power * dt / 3600, self.max_charge_rate * dt / 3600)
self.state.battery_charge = min(self.battery_capacity,
self.state.battery_charge + charge_increment)
else:
discharge_increment = min(abs(net_power) * dt / 3600, self.max_discharge_rate * dt / 3600)
self.state.battery_charge = max(0, self.state.battery_charge - discharge_increment)
# Update battery temperature
if abs(net_power) > 300:
self.state.battery_temp += 0.1 * dt
else:
self.state.battery_temp -= 0.05 * dt
self.state.time += dt
self.history.append(self.state)
def detect_power_attacks(self) -> List[str]:
"""Detect potential power system attacks"""
anomalies = []
# Check battery state of charge
soc = self.state.battery_charge / self.battery_capacity
if soc < self.normal_battery_range[0]:
anomalies.append("Critical: Battery below minimum")
elif soc > self.normal_battery_range[1]:
anomalies.append("Warning: Battery overcharging")
# Check temperature
if self.state.battery_temp < self.normal_temp_range[0]:
anomalies.append("Warning: Battery temperature too low")
elif self.state.battery_temp > self.normal_temp_range[1]:
anomalies.append("Critical: Battery overheating")
# Check load anomalies
if self.state.load_power > 800:
anomalies.append("Warning: Load power excessive")
return anomalies
# Example usage
simulator = PowerSystemSimulator()
print("Initial state:")
print(f" Battery: {simulator.state.battery_charge:.0f} Wh")
print(f" Load: {simulator.state.load_power:.0f} W")
# Normal operation
for _ in range(50):
simulator.update_state(dt=0.1)
# Attack scenario
print("\nSimulating eclipse scenario...")
simulator.state.eclipse = True
for _ in range(100):
simulator.update_state(dt=0.1)
anomalies = simulator.detect_power_attacks()
if anomalies:
print(f"Anomalies detected: {anomalies}")
print(f"\nFinal battery: {simulator.state.battery_charge:.0f} Wh")
Complete spacecraft command security framework implementing CCSDS structure validation, CRC-16-CCITT checksums, command injection detection, and rate limiting. Covers semantic validation, anomalous command sequence detection, and command authorization enforcement. Includes practical Python validation engine that distinguishes legitimate commands from injection attacks, with comprehensive logging of all command receipt, parsing, and execution for forensic analysis.
CCSDS Command Structure Validator:
import re
import hashlib
from typing import List, Dict, Tuple
class SpacecraftCommandValidator:
"""
Validate spacecraft commands for potential injection attacks
Implements CCSDS command standards with security extensions
"""
def __init__(self):
# CCSDS command structure
self.command_structure = {
'primary_header': 6,
'data_field_header': 1,
'command_code': 2,
'parameters': 'variable',
'checksum': 2
}
# Allowed commands and parameters
self.allowed_commands = {
'0x01': {'name': 'POWER_ON', 'params': ['device_id', 'delay_ms']},
'0x02': {'name': 'POWER_OFF', 'params': ['device_id', 'delay_ms']},
'0x03': {'name': 'SET_MODE', 'params': ['mode_id', 'flags']},
'0x04': {'name': 'EXECUTE_SEQUENCE', 'params': ['seq_id', 'priority']},
}
self.command_history = []
self.max_history = 1000
self.command_auth_key = b'sample_key_32_bytes_for_aes_256'
def validate_command_structure(self, command_bytes: bytes) -> Tuple[bool, str]:
"""Validate command structure against CCSDS standards"""
if len(command_bytes) < 11: # Minimum CCSDS command length
return False, "Command too short"
# Extract fields
primary_header = command_bytes[0:6]
data_field_header = command_bytes[6]
command_code = command_bytes[7:9]
parameters = command_bytes[9:-2]
checksum = command_bytes[-2:]
# Verify checksum
calculated_checksum = self.calculate_checksum(command_bytes[:-2])
if calculated_checksum != checksum:
return False, "Checksum mismatch"
# Verify command code exists
cmd_code_hex = command_code.hex()
if cmd_code_hex not in self.allowed_commands:
return False, "Unknown command code"
return True, "Structure valid"
def validate_command_semantics(self, command_bytes: bytes) -> Tuple[bool, str, Dict]:
"""Validate command parameters and semantics"""
# Parse command
command_code = command_bytes[7:9].hex()
parameters = command_bytes[9:-2]
cmd_info = self.allowed_commands[command_code]
# Parse parameters
param_values = []
for i in range(0, len(parameters), 2):
param_val = int.from_bytes(parameters[i:i+2], 'big')
param_values.append(param_val)
# Validate semantics
if not self.check_rate_limit(command_code):
return False, "Rate limit exceeded", {}
if self.is_dangerous_sequence(command_code, param_values):
return False, "Dangerous command sequence detected", {}
return True, "Semantics valid", {
'command': cmd_info['name'],
'parameters': param_values
}
def detect_command_injection(self, command_bytes: bytes) -> List[str]:
"""Detect potential command injection attempts"""
anomalies = []
try:
command_str = command_bytes.decode('ascii', errors='ignore')
except:
anomalies.append("Non-ASCII data detected")
# Common injection patterns
injection_patterns = [
(r'DROP\s+TABLE', "SQL injection pattern"),
(r'EXEC\s*\(', "Execution injection pattern"),
(r'[\x00\x01\x02].*[\xff]', "Unusual byte pattern"),
]
for pattern, description in injection_patterns:
if re.search(pattern, command_str, re.IGNORECASE):
anomalies.append(f"Detected: {description}")
# Buffer overflow detection
if len(command_bytes) > 1024:
anomalies.append("Command excessively long")
return anomalies
def calculate_checksum(self, data: bytes) -> bytes:
"""Calculate CRC-16-CCITT checksum"""
crc = 0xFFFF
polynomial = 0x1021
for byte in data:
crc ^= byte << 8
for _ in range(8):
crc <<= 1
if crc & 0x10000:
crc ^= polynomial
crc &= 0xFFFF
return crc.to_bytes(2, 'big')
def is_dangerous_sequence(self, command_code: str, params: List[int]) -> bool:
"""Detect dangerous command sequences"""
# Check command history
recent_commands = self.command_history[-10:] if len(self.command_history) > 10 else self.command_history
# Example: POWER_OFF immediately after critical operation
if command_code == '0x02' and params[1] == 0: # POWER_OFF
for prev_cmd in recent_commands:
if prev_cmd['code'] == '0x04': # After EXECUTE_SEQUENCE
return True
return False
def check_rate_limit(self, command_code: str) -> bool:
"""Implement rate limiting for commands"""
import time
current_time = time.time()
# Remove old entries
self.command_history = [cmd for cmd in self.command_history
if current_time - cmd['time'] < 60]
# Count recent occurrences
recent_count = sum(1 for cmd in self.command_history
if cmd['code'] == command_code)
# Rate limit: max 10 per minute per command
return recent_count < 10
# Example usage
validator = SpacecraftCommandValidator()
# Test valid command
valid_command = bytes.fromhex('1ACFFC1D00010200030004A55A')
structure_valid, structure_msg = validator.validate_command_structure(valid_command)
print(f"Structure validation: {structure_valid} - {structure_msg}")
if structure_valid:
semantic_valid, semantic_msg, parsed_cmd = validator.validate_command_semantics(valid_command)
print(f"Semantic validation: {semantic_valid} - {semantic_msg}")
# Test injection detection
malicious_command = b'SET_MODE 1; DROP TABLE commands; --'
injection_indicators = validator.detect_command_injection(malicious_command)
if injection_indicators:
print(f"Injection detected:")
for indicator in injection_indicators:
print(f" - {indicator}")
Complete attitude control system simulation using quaternion kinematics, Euler's equation, and PD control laws. Implements Extended Kalman Filter (EKF) for attitude estimation with gyroscope bias correction and attack detection via innovation monitoring. Covers 3-sigma consistency checks, numerically stable quaternion operations, and covariance matrix updates. Essential for understanding ADCS vulnerability to sensor spoofing and demonstrating detection of malicious attitude commands.
Complete Attitude Control System Simulation:
import numpy as np
from scipy.integrate import solve_ivp
from dataclasses import dataclass
@dataclass
class SatelliteState:
"""Complete satellite state for attitude simulation"""
quaternion: np.ndarray # [q0, q1, q2, q3]
angular_velocity: np.ndarray # rad/s
time: float
control_torque: np.ndarray
class AttitudeDynamics:
"""
Complete satellite attitude dynamics simulation
Includes disturbance torques and control algorithms
"""
def __init__(self):
# Satellite inertia tensor (kg·m²)
self.J = np.diag([100, 150, 120])
# Control parameters
self.kp = 5.0 # Proportional gain
self.kd = 10.0 # Derivative gain
self.max_torque = 0.1 # N·m
# Reference attitude (target)
self.reference_quaternion = np.array([1, 0, 0, 0]) # Identity
self.state_history = []
def quaternion_kinematics(self, q: np.ndarray, omega: np.ndarray) -> np.ndarray:
"""Quaternion time derivative"""
q0, q1, q2, q3 = q
wx, wy, wz = omega
omega_quat = np.array([
[-q1, -q2, -q3],
[q0, -q3, q2],
[q3, q0, -q1],
[-q2, q1, q0]
])
return 0.5 * omega_quat @ omega
def dynamics(self, t: float, y: np.ndarray) -> np.ndarray:
"""Complete attitude dynamics"""
q = y[:4]
omega = y[4:]
# Normalize quaternion
q = q / np.linalg.norm(q)
# Calculate control torque
control_torque = self.calculate_control_torque(q, omega)
# Total torque
total_torque = control_torque
# Quaternion kinematics
dqdt = self.quaternion_kinematics(q, omega)
# Euler's equation: J·dω/dt = -ω×Jω + τ
Jomega = self.J @ omega
omega_cross_Jomega = np.cross(omega, Jomega)
dwdt = np.linalg.solve(self.J, -omega_cross_Jomega + total_torque)
return np.concatenate([dqdt, dwdt])
def calculate_control_torque(self, q: np.ndarray, omega: np.ndarray) -> np.ndarray:
"""PD control law"""
# Calculate error quaternion
q_ref_inv = self.quaternion_inverse(self.reference_quaternion)
q_error = self.quaternion_multiply(q_ref_inv, q)
# Extract vector part
q_error_vec = q_error[1:]
# PD control
control_torque = -self.kp * q_error_vec - self.kd * omega
# Limit torque
torque_norm = np.linalg.norm(control_torque)
if torque_norm > self.max_torque:
control_torque = (self.max_torque / torque_norm) * control_torque
return control_torque
def simulate(self, duration: float = 100):
"""Simulate attitude dynamics"""
# Initial conditions
q0 = np.array([1, 0, 0, 0])
omega0 = np.array([0.01, 0.005, 0.002])
y0 = np.concatenate([q0, omega0])
# Solve ODE
sol = solve_ivp(
self.dynamics, (0, duration), y0,
t_eval=np.linspace(0, duration, 1000),
method='RK45', dense_output=True, max_step=0.1
)
return sol
@staticmethod
def quaternion_multiply(q1: np.ndarray, q2: np.ndarray) -> np.ndarray:
"""Multiply two quaternions"""
a1, b1, c1, d1 = q1
a2, b2, c2, d2 = q2
return np.array([
a1*a2 - b1*b2 - c1*c2 - d1*d2,
a1*b2 + b1*a2 + c1*d2 - d1*c2,
a1*c2 - b1*d2 + c1*a2 + d1*b2,
a1*d2 + b1*c2 - c1*b2 + d1*a2
])
@staticmethod
def quaternion_inverse(q: np.ndarray) -> np.ndarray:
"""Inverse of a unit quaternion"""
return np.array([q[0], -q[1], -q[2], -q[3]])
# Example usage
dynamics = AttitudeDynamics()
print("Simulating satellite attitude control...")
sol = dynamics.simulate(duration=50)
# Calculate pointing error
q_final = sol.y[:4, -1]
error = 2 * np.arccos(min(1, abs(q_final[0])))
print(f"Final attitude error: {np.degrees(error):.2f} degrees")
Extended Kalman Filter for Attitude Estimation:
import numpy as np
from scipy.linalg import inv
from typing import List
class AttitudeEKF:
"""
Extended Kalman Filter for satellite attitude estimation
Includes attack detection capabilities
"""
def __init__(self):
# State: [q0, q1, q2, q3, ωx, ωy, ωz, bx, by, bz]
self.state_dim = 10
self.state = np.zeros(self.state_dim)
self.state[0] = 1.0 # Identity quaternion
# Covariance matrix
self.P = np.eye(self.state_dim) * 0.01
# Process noise covariance
self.Q = np.eye(self.state_dim) * 1e-6
# Measurement noise covariance
self.R_gyro = np.eye(3) * 1e-4
self.R_star = np.eye(2) * 1e-3
# Gyro parameters
self.gyro_bias = np.zeros(3)
# Attack detection
self.innovation_history = []
self.max_innovation_history = 100
self.innovation_threshold = 5.0
def predict(self, dt: float):
"""Prediction step of EKF"""
# Extract states
q = self.state[:4]
omega = self.state[4:7]
bias = self.state[7:10]
# True angular rate
omega_true = omega - bias
# Predict quaternion
omega_norm = np.linalg.norm(omega_true)
if omega_norm > 1e-10:
axis = omega_true / omega_norm
angle = omega_norm * dt
dq = np.array([
np.cos(angle/2),
axis[0] * np.sin(angle/2),
axis[1] * np.sin(angle/2),
axis[2] * np.sin(angle/2)
])
q_pred = self.quaternion_multiply(dq, q)
else:
q_pred = q
# Predict other states
omega_pred = omega
bias_pred = bias
# Update state
self.state[:4] = q_pred
self.state[4:7] = omega_pred
self.state[7:10] = bias_pred
# Calculate Jacobian
F = np.eye(self.state_dim)
# Update covariance
self.P = F @ self.P @ F.T + self.Q
# Normalize quaternion
self.state[:4] = self.state[:4] / np.linalg.norm(self.state[:4])
def update_gyro(self, gyro_measurement: np.ndarray):
"""Update step with gyro measurement"""
# Measurement model
H = np.zeros((3, self.state_dim))
H[:, 4:7] = np.eye(3)
H[:, 7:10] = np.eye(3)
# Predicted measurement
z_pred = self.state[4:7] + self.state[7:10]
# Innovation
innovation = gyro_measurement - z_pred
# Innovation covariance
S = H @ self.P @ H.T + self.R_gyro
# Kalman gain
K = self.P @ H.T @ inv(S)
# Update state
self.state = self.state + K @ innovation
# Update covariance
self.P = (np.eye(self.state_dim) - K @ H) @ self.P
# Store innovation for attack detection
self.innovation_history.append({
'innovation': innovation,
'covariance': S,
'time': len(self.innovation_history)
})
if len(self.innovation_history) > self.max_innovation_history:
self.innovation_history.pop(0)
@staticmethod
def quaternion_multiply(q1: np.ndarray, q2: np.ndarray) -> np.ndarray:
"""Multiply two quaternions"""
a1, b1, c1, d1 = q1
a2, b2, c2, d2 = q2
return np.array([
a1*a2 - b1*b2 - c1*c2 - d1*d2,
a1*b2 + b1*a2 + c1*d2 - d1*c2,
a1*c2 - b1*d2 + c1*a2 + d1*b2,
a1*d2 + b1*c2 - c1*b2 + d1*a2
])
# Example usage
ekf = AttitudeEKF()
print("Running EKF for attitude estimation...")
for i in range(100):
ekf.predict(dt=0.1)
# Simulated gyro measurement
true_omega = np.array([0.01, 0.02, 0.005])
gyro_measurement = true_omega + np.random.normal(0, 0.001, 3)
ekf.update_gyro(gyro_measurement)
if i % 20 == 0:
attitude_error = 2 * np.arccos(min(1, abs(ekf.state[0])))
print(f"Step {i}: Attitude error = {np.degrees(attitude_error):.2f}°")
print("EKF estimation complete.")
12. Attitude Dynamics & Extended Kalman Filter
Complete Attitude Control System Simulation:
Additional Resources & References:
- CCSDS (Consultative Committee for Space Data Systems) Standards: https://public.ccsds.org
- NORAD Satellite Tracking Data: https://celestrak.com/
- SGP4 Reference Implementation: github.com/skyfielddateutil/skyfield
- EGM2008 Gravity Model: http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/
- Skyfield Library Documentation: https://rhodesmill.org/skyfield/
- NASA Orbital Mechanics Course: https://www.nasa.gov/
Comments
Post a Comment