Source code for stonesoup.sensor.radar.radar

# -*- coding: utf-8 -*-
import copy
from math import erfc
from typing import Tuple, Set, Union

import numpy as np
import scipy.constants as const

from .beam_pattern import BeamTransitionModel
from .beam_shape import BeamShape
from ...base import Property
from ...functions import cart2sphere, rotx, roty, rotz, mod_bearing
from ...models.measurement.base import MeasurementModel
from ...models.measurement.nonlinear import \
    (CartesianToBearingRange, CartesianToElevationBearingRange,
     CartesianToBearingRangeRate, CartesianToElevationBearingRangeRate)
from ...sensor.sensor import Sensor
from ...types.array import CovarianceMatrix
from ...types.detection import TrueDetection
from ...types.groundtruth import GroundTruthState
from ...types.numeric import Probability
from ...types.state import State, StateVector
from ...models.clutter.clutter import ClutterModel


[docs]class RadarBearingRange(Sensor): """A simple radar sensor that generates measurements of targets, using a :class:`~.CartesianToBearingRange` model, relative to its position. Note ---- The current implementation of this class assumes a 3D Cartesian plane. """ ndim_state: int = Property( default=2, doc="Number of state dimensions. This is utilised by (and follows in format) " "the underlying :class:`~.CartesianToBearingRange` model") position_mapping: Tuple[int, int] = Property( doc="Mapping between the targets state space and the sensors " "measurement capability") noise_covar: CovarianceMatrix = Property( doc="The sensor noise covariance matrix. This is utilised by " "(and follow in format) the underlying " ":class:`~.CartesianToBearingRange` model") clutter_model: ClutterModel = Property( default=None, doc="An optional clutter generator that adds a set of simulated " ":class:`Clutter` ojects to the measurements at each time step. " "The clutter is simulated according to the provided distribution.")
[docs] def measure(self, ground_truths: Set[GroundTruthState], noise: Union[np.ndarray, bool] = True, **kwargs) -> Set[TrueDetection]: measurement_model = CartesianToBearingRange( ndim_state=self.ndim_state, mapping=self.position_mapping, noise_covar=self.noise_covar, translation_offset=self.position, rotation_offset=self.orientation) detections = set() for truth in ground_truths: measurement_vector = measurement_model.function(truth, noise=noise, **kwargs) detection = TrueDetection(measurement_vector, measurement_model=measurement_model, timestamp=truth.timestamp, groundtruth_path=truth) detections.add(detection) # Generate clutter at this time step if self.clutter_model is not None: self.clutter_model.measurement_model = measurement_model clutter = self.clutter_model.function(ground_truths) detections = set.union(detections, clutter) return detections
[docs]class RadarRotatingBearingRange(RadarBearingRange): """A simple rotating radar, with set field-of-view (FOV) angle, range and\ rotations per minute (RPM), that generates measurements of targets, using\ a :class:`~.CartesianToBearingRange` model, relative to its\ position. Note ---- * The current implementation of this class assumes a 3D Cartesian plane. """ dwell_center: State = Property( doc="A state object, whose `state_vector` " "property describes the rotation angle of the center of the sensor's " "current FOV (i.e. the dwell center) relative to the positive x-axis " "of the sensor frame/orientation. The angle is positive if the rotation " "is in the counter-clockwise direction when viewed by an observer " "looking down the z-axis of the sensor frame, towards the origin. " "Angle units are in radians" ) rpm: float = Property(doc="The number of antenna rotations per minute (RPM)") max_range: float = Property(doc="The maximum detection range of the radar (in meters)") fov_angle: float = Property(doc="The radar field of view (FOV) angle (in radians).")
[docs] def measure(self, ground_truths: Set[GroundTruthState], noise: Union[np.ndarray, bool] = True, **kwargs) -> Set[TrueDetection]: # Read timestamp from ground truth try: timestamp = next(iter(ground_truths.copy())).timestamp except StopIteration: # No ground truths to get timestamp from return set() # Rotate the radar antenna and compute new heading self.rotate(timestamp) antenna_heading = self.orientation[2, 0] + self.dwell_center.state_vector[0, 0] # Set rotation offset of underlying measurement model rot_offset = \ StateVector( [[self.orientation[0, 0]], [self.orientation[1, 0]], [antenna_heading]]) measurement_model = CartesianToBearingRange( ndim_state=self.ndim_state, mapping=self.position_mapping, noise_covar=self.noise_covar, translation_offset=self.position, rotation_offset=rot_offset) detections = set() for truth in ground_truths: # Transform state to measurement space and generate # random noise measurement_vector = measurement_model.function(truth, noise=False, **kwargs) if noise is True: measurement_noise = measurement_model.rvs() else: measurement_noise = noise # Check if state falls within sensor's FOV fov_min = -self.fov_angle / 2 fov_max = +self.fov_angle / 2 bearing_t = measurement_vector[0, 0] range_t = measurement_vector[1, 0] # Do not measure if state not in FOV if (bearing_t > fov_max or bearing_t < fov_min or range_t > self.max_range): continue # Else add measurement measurement_vector += measurement_noise # Add noise detection = TrueDetection(measurement_vector, measurement_model=measurement_model, timestamp=truth.timestamp, groundtruth_path=truth) detections.add(detection) return detections
[docs] def rotate(self, timestamp): """Rotate the sensor's antenna This method computes and updates the sensor's `dwell_center` property. Parameters ---------- timestamp: :class:`datetime.datetime` A timestamp signifying when the rotation completes """ # Check if dwell_center has a timestamp instantiated if not sets it to incoming timestamp if self.dwell_center.timestamp is None: self.dwell_center.timestamp = timestamp # Compute duration since last rotation duration = timestamp - self.dwell_center.timestamp # Update dwell center rps = self.rpm / 60 # rotations per sec angle = self.dwell_center.state_vector[0, 0] + duration.total_seconds()*rps*2*np.pi self.dwell_center = State(StateVector([[mod_bearing(angle)]]), timestamp)
[docs]class RadarElevationBearingRange(RadarBearingRange): """A radar sensor that generates measurements of targets, using a :class:`~.CartesianToElevationBearingRange` model, relative to its position. Note ---- This implementation of this class assumes a 3D Cartesian space. """ ndim_state: int = Property( default=3, doc="Number of state dimensions. This is utilised by (and follows in format) " "the underlying :class:`~.CartesianToBearingRange` model") noise_covar: CovarianceMatrix = Property( doc="The sensor noise covariance matrix. This is utilised by " "(and follow in format) the underlying " ":class:`~.CartesianToElevationBearingRange` model")
[docs] def measure(self, ground_truths: Set[GroundTruthState], noise: Union[np.ndarray, bool] = True, **kwargs) -> Set[TrueDetection]: measurement_model = CartesianToElevationBearingRange( ndim_state=self.ndim_state, mapping=self.position_mapping, noise_covar=self.noise_covar, translation_offset=self.position, rotation_offset=self.orientation) detections = set() for truth in ground_truths: measurement_vector = measurement_model.function(truth, noise=noise, **kwargs) detection = TrueDetection(measurement_vector, measurement_model=measurement_model, timestamp=truth.timestamp, groundtruth_path=truth) detections.add(detection) # Generate clutter at this time step if self.clutter_model is not None: self.clutter_model.measurement_model = measurement_model clutter = self.clutter_model.function(ground_truths) detections = set.union(detections, clutter) return detections
[docs]class RadarBearingRangeRate(RadarBearingRange): """ A radar sensor that generates measurements of targets, using a :class:`~.CartesianToBearingRangeRate` model, relative to its position and velocity. Note ---- This class implementation assuming at 3D cartesian space, it therefore\ expects a 6D state space. """ velocity_mapping: Tuple[int, int, int] = Property( default=(1, 3, 5), doc="Mapping to the target's velocity information within its state space") ndim_state: int = Property( default=3, doc="Number of state dimensions. This is utilised by (and follows in format) " "the underlying :class:`~.CartesianToBearingRangeRate` model") noise_covar: CovarianceMatrix = Property( doc="The sensor noise covariance matrix. This is utilised by " "(and follow in format) the underlying " ":class:`~.CartesianToBearingRangeRate` model")
[docs] def measure(self, ground_truths: Set[GroundTruthState], noise: Union[np.ndarray, bool] = True, **kwargs) -> Set[TrueDetection]: measurement_model = CartesianToBearingRangeRate( ndim_state=self.ndim_state, mapping=self.position_mapping, velocity_mapping=self.velocity_mapping, noise_covar=self.noise_covar, translation_offset=self.position, velocity=self.velocity, rotation_offset=self.orientation) detections = set() for truth in ground_truths: measurement_vector = measurement_model.function(truth, noise=noise, **kwargs) detection = TrueDetection(measurement_vector, measurement_model=measurement_model, timestamp=truth.timestamp, groundtruth_path=truth) detections.add(detection) return detections
[docs]class RadarElevationBearingRangeRate(RadarBearingRangeRate): """ A radar sensor that generates measurements of targets, using a :class:`~.CartesianToElevationBearingRangeRate` model, relative to its position and velocity. Note ---- The current implementation of this class assumes a 3D Cartesian plane. """ velocity_mapping: Tuple[int, int, int] = Property( default=(1, 3, 5), doc="Mapping to the target's velocity information within its state space") ndim_state: int = Property( default=6, doc="Number of state dimensions. This is utilised by (and follows in format) " "the underlying :class:`~.CartesianToElevationBearingRangeRate` model") noise_covar: CovarianceMatrix = Property( doc="The sensor noise covariance matrix. This is utilised by " "(and follow in format) the underlying " ":class:`~.CartesianToElevationBearingRangeRate` model")
[docs] def measure(self, ground_truths: Set[GroundTruthState], noise: Union[np.ndarray, bool] = True, **kwargs) -> Set[TrueDetection]: measurement_model = CartesianToElevationBearingRangeRate( ndim_state=self.ndim_state, mapping=self.position_mapping, velocity_mapping=self.velocity_mapping, noise_covar=self.noise_covar, translation_offset=self.position, velocity=self.velocity, rotation_offset=self.orientation) detections = set() for truth in ground_truths: measurement_vector = measurement_model.function(truth, noise=noise, **kwargs) detection = TrueDetection(measurement_vector, measurement_model=measurement_model, timestamp=truth.timestamp, groundtruth_path=truth) detections.add(detection) return detections
[docs]class RadarRasterScanBearingRange(RadarRotatingBearingRange): """A simple raster scan radar, with set field-of-regard (FoR) angle, \ field-of-view (FoV) angle, range and rotations per minute (RPM), that \ generates measurements of targets, using a \ :class:`~.CartesianToBearingRange` model, relative to its position This is a simple extension of the RadarRotatingBearingRange class with \ the rotate function changed to restrict the dwell-center to within the \ field of regard. It's important to note that this only works (has been tested) in an 2D \ environment Note ---- This class implementation assuming at 3D cartesian space, it therefore\ expects a 6D state space. """ for_angle: float = Property(doc="The radar field of regard (FoR) angle (in radians).")
[docs] def rotate(self, timestamp): """Rotate the sensor's antenna This method computes and updates the sensor's `dwell_center` property. Parameters ---------- timestamp: :class:`datetime.datetime` A timestamp signifying when the rotation completes """ super().rotate(timestamp) dwell_center_max = self.for_angle/2.0 - self.fov_angle/2.0 dwell_center_min = -self.for_angle/2.0 + self.fov_angle/2.0 # If the FoV is outside of the FoR: # Correct the dwell_center # Reverse the direction of the scan pattern if self.dwell_center.state_vector[0, 0] > dwell_center_max: self.dwell_center = State( StateVector([[(2.0 * dwell_center_max) - self.dwell_center.state_vector[0, 0] ]]), timestamp) self.rpm = -self.rpm elif self.dwell_center.state_vector[0, 0] < dwell_center_min: self.dwell_center = State( StateVector([[(2.0 * dwell_center_min) - self.dwell_center.state_vector[0, 0] ]]), timestamp) self.rpm = -self.rpm
[docs]class AESARadar(Sensor): r"""An AESA (Active electronically scanned array) radar model that calculates the signal to noise ratio (SNR) of a target and the subsequent probability of detection (PD). The SNR is calculated using: .. math:: \mathit{SNR} = \dfrac{c^2\, n_p \,\beta}{64\,\pi^3 \,kT_0 \,B \,F \,f^2 \,L} \times\dfrac{\sigma\, G_a^2\, P_t}{R^4} where :math:`c` is the speed of light :math:`n_p` is the number of pulses in a burst :math:`\beta` is the duty cycle which is unitless :math:`k` is the boltzmann constant :math:`T_0` is system temperature in kelvin :math:`B` is the bandwidth in hertz :math:`F` is the noise figure unitless :math:`f` is the frequency in hertz :math:`L` is the loss which is unitless The probability of detection (:math:`P_{d}`) is calculated using the North's approximation, .. math:: P_{d} = 0.5\, erfc\left(\sqrt{-\ln{P_{fa}}}-\sqrt{ \mathit{SNR} +0.5} \right) where :math:`P_{fa}` is the probability of false alarm. In this model the AESA scan angle effects the gain by: .. math:: G_a = G_a \cos{\left(\theta\right)} \cos{\left(\phi\right)} where :math:`\theta` and :math:`\phi` are respectively the azimuth and elevation angles in respects to the boresight of the antenna. The effect of beam spoiling on the beam width is : .. math:: \Delta\theta = \dfrac{\Delta\theta}{\cos{\left(\theta\right)} \cos{\left(\phi\right)}} Note ---- The current implementation of this class assumes a 3D Cartesian plane. This model does not generate false alarms. """ rotation_offset: StateVector = Property( default=None, doc="A 3x1 array of angles (rad), specifying the radar orientation in terms of the " "counter-clockwise rotation around the :math:`x,y,z` axis. i.e Roll, Pitch and Yaw. " "Default is ``StateVector([0, 0, 0])``") position_mapping: Tuple[int, int, int] = Property( default=(0, 1, 2), doc="Mapping between or positions and state " "dimensions. [x,y,z]") measurement_model: MeasurementModel = Property( doc="The Measurement model used to generate " "measurements.") beam_shape: BeamShape = Property( doc="Object describing the shape of the beam.") beam_transition_model: BeamTransitionModel = Property( doc="Object describing the movement of the beam in azimuth and elevation from the " "perspective of the radar.") # SNR variables number_pulses: int = Property( default=1, doc="The number of pulses in the radar burst.") duty_cycle: float = Property( doc="Duty cycle is the fraction of the time the radar it transmitting.") band_width: float = Property( doc="Bandwidth of the receiver in hertz.") receiver_noise: float = Property( doc="Noise figure of the radar in decibels.") frequency: float = Property( doc="Transmitted frequency in hertz.") antenna_gain: float = Property( doc="Total Antenna gain in decibels.") beam_width: float = Property( doc="Radar beam width in radians.") loss: float = Property( default=0, doc="Loss in decibels.") swerling_on: bool = Property( default=False, doc="Is the Swerling 1 case used. If True the RCS" " of the target will change for each timestep. " "The random RCS follows the probability " "distribution of the Swerling 1 case.") rcs: float = Property( default=None, doc="The radar cross section of targets in meters squared. Used if rcs not present on " "truth. Default `None`, where 'rcs' must be present on truth.") probability_false_alarm: Probability = Property( default=1e-6, doc="Probability of false alarm used in the North's approximation") def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) if self.rotation_offset is None: self.rotation_offset = StateVector([0, 0, 0]) @property def _snr_constant(self): temp = 290 # noise reference temperature (room temperature kelvin) # convert from dB noise_figure = 10 ** (self.receiver_noise / 10) loss = 10 ** (self.loss / 10) # calculate part of snr that is independent of: # rcs, transmitted power, gain and range return (const.c ** 2 * self.number_pulses * self.duty_cycle) / \ (64 * np.pi ** 3 * const.k * temp * self.band_width * noise_figure * self.frequency ** 2 * loss) @property def _rotation_matrix(self): """_rotation_matrix getter method Calculates and returns the (3D) axis rotation matrix. Returns ------- : :class:`numpy.ndarray` of shape (3, 3) The model (3D) rotation matrix. """ theta_x = -self.rotation_offset[0, 0] # roll theta_y = -self.rotation_offset[1, 0] # pitch#elevation theta_z = -self.rotation_offset[2, 0] # yaw#azimuth return rotz(theta_z) @ roty(theta_y) @ rotx(theta_x) @staticmethod def _swerling_1(rcs): return -rcs * np.log(np.random.rand())
[docs] def gen_probability(self, truth): """Generates probability of detection of a given State. Parameters ---------- truth : The target state. Returns ------- det_prob : `float` Probability of detection. snr : `float` Signal to noise ratio. rcs : `float` RCS after the Swerling 1 case is applied. directed_power : `float` Power in the direction of the target. spoiled_gain : `float` Gain in decibels with beam spoiling applied. spoiled_width : `float` Beam width with beam spoiling applied. """ if getattr(truth, 'rcs', None) is not None: # use state's rcs if it has one rcs = truth.rcs else: rcs = self.rcs if rcs is None: raise ValueError("Truth missing 'rcs' attribute and no default 'rcs' provided") # apply swerling 1 case? if self.swerling_on: rcs = self._swerling_1(rcs) # e-scan beam steer [beam_az, beam_el] = self.beam_transition_model.move_beam( truth.timestamp) # [az,el] # effects of e-scan on gain and beam width spoiled_gain = 10 ** (self.antenna_gain / 10) * np.cos(beam_az) * np.cos(beam_el) spoiled_width = self.beam_width / (np.cos(beam_az) * np.cos(beam_el)) # state relative to radar (in cartesian space) relative_vector = truth.state_vector[self.position_mapping, :] - self.position relative_vector = self._rotation_matrix @ relative_vector # calculate target position in spherical coordinates [r, pos_az, pos_el] = cart2sphere(*relative_vector) # target position relative to beam position relative_az = pos_az - beam_az relative_el = pos_el - beam_el # calculate power directed towards target directed_power = self.beam_shape.beam_power(relative_az, relative_el, spoiled_width) # calculate signal to noise ratio snr = self._snr_constant * rcs * spoiled_gain ** 2 * directed_power / (r ** 4) # calculate probability of detection using the North's approximation det_prob = 0.5 * erfc( (-np.log(self.probability_false_alarm)) ** 0.5 - ( snr + 1 / 2) ** 0.5) return det_prob, snr, rcs, directed_power, 10 * np.log10(spoiled_gain), spoiled_width
[docs] def measure(self, ground_truths: Set[GroundTruthState], noise: Union[np.ndarray, bool] = True, **kwargs) -> Set[TrueDetection]: detections = set() measurement_model = copy.deepcopy(self.measurement_model) measurement_model.translation_offset = self.position.copy() measurement_model.rotation_offset = self.rotation_offset.copy() for truth in ground_truths: det_prob = self.gen_probability(truth)[0] # Is the state detected? if np.random.rand() <= det_prob: measured_pos = measurement_model.function(truth, noise=noise) detection = TrueDetection(measured_pos, timestamp=truth.timestamp, measurement_model=measurement_model, groundtruth_path=truth) detections.add(detection) return detections