Maritime Target Tracking Using Radar Under Rough Sea Conditions
This example demonstrates how to track targets in a maritime environment under rough sea conditions using a mechanically rotating radar. You will simulate target measurements and sea clutter in these challenging conditions with a statistical radar model. Additionally, you will learn how to define suitable models for specifying boats and marine radar for the tracking algorithm.
Introduction
The sea surface, particularly under rough conditions, reflects radar signals back strongly towards the radar, resulting in a cluttered background that can obscure the detection of actual targets. The intensity of sea clutter is significantly higher at shorter ranges due to lower path loss and relatively larger grazing angles between the sea surface and the radar beam. This clutter can easily be confused as boats with a smaller radar cross section (RCS). By adjusting the radar gain parameters, one can reduce sea clutter during detection, but this is often balanced with the risk of missing the detection of small boats at close ranges. In this example, you will learn how to use a tracker to monitor and track small boats as they approach the radar in rough sea conditions.
Setup Scenario
To simulate a maritime environment, use the radarScenario class to set up a scenario. Set the UpdateRate to 0, allowing the sensors to control the simulation rate.
rng('default'); % Reproducible results % Create a scenario scene = radarScenario('UpdateRate',0,'StopTime',200);
Next, add a sea surface with a radius of 5 kilometers is specified around the origin, simulating rough sea conditions with a sea state of 8. For further details on sea surface modeling, refer to the Maritime Radar Sea Clutter Modeling (Radar Toolbox) example.
% Add sea surface refl = surfaceReflectivitySea('SeaState',8,'Speckle','Weibull'); srf = seaSurface(scene,'RadarReflectivity',refl,'Boundary',[-5e3 5e3;-5e3 5e3]);
Add a platform to the scenario to represent the ownship or host boat, on which the radar is mounted.
ownship = platform(scene);
Define a mechanically rotating radar rotating at 20 revolutions per minute. The radar has an instantaneous field of view of 10 degrees in azimuth and 25 degrees in elevation. The radar is mounted 24 meters above the sea surface, aligned parallel to ownship. Set the FalseAlarmRate to 1e-3 to detect targets even with small signal-to-noise ratio (SNR). Set HasFalseAlarms to false to only simulate false returns from the sea clutter. For more details on the statistical radar model, refer to the algorithms section of radarDataGenerator (Radar Toolbox) class.
radar = radarDataGenerator('SensorIndex', 1, ... 'UpdateRate', 12, ... 'ScanMode','Mechanical',... 'MountingLocation', [0 0 24], ... 'MountingAngles', [0 0 0], ... 'FieldOfView', [10 25], ... 'RangeLimits', [500 5e3],... 'HasFalseAlarms',false,... 'HasNoise',true,... 'HasElevation',false,... 'HasRangeRate',false,... 'FalseAlarmRate', 1e-3, ... 'DetectionProbability',0.98,... 'ReferenceRCS',10,... 'ReferenceRange', 3e3, ... 'DetectionCoordinates', 'Body', ... 'AzimuthResolution', 4, ... 'AzimuthBiasFraction',1,... 'RangeBiasFraction',1,... 'RangeResolution', 10, ... 'MaxAzimuthScanRate',360);
Attach the radar to the ownship boat.
ownship.Sensors = radar;
Enable clutter generation from the sea surface using the clutterGenerator (Radar Toolbox) function with the radar. The UseBeam option is set to true to simulate clutter within the radar's field of view, as defined in the previous step.
clt = clutterGenerator(scene,radar,'Resolution',100,'UseBeam',true);
Simulate Clutter Measurements
To observe the clutter measurements from the sea surface reported by the radar, run the scenario for a few seconds without any targets. This helps to understand how the radar interacts with the sea to generate clutter measurements distributed across the entire radar coverage. A helper function, createDisplay, defined in the Supporting Functions section is used to visualize the radar measurements and coverage.
% Create display for plotting display = createDisplay(scene); restart(scene); % Run the scenario for 10 scans detBuffer = {}; nScans = 10; while advance(scene) && scene.SimulationTime < (nScans-1)*3 detections = detect(scene); detBuffer = [detBuffer;detections]; updateDisplay(display, scene, detBuffer); end
Notice in the image above that there is heavy clutter near the radar and the rest of the field of view seems uniformly cluttered.
Clutter Distribution
Next, quantitatively inspect this distribution using the simulated sensor data recorded in the first 10 scans of the radar.
% Compute azimuth and range of measurements
x = cellfun(@(x)x.Measurement(1),detBuffer);
y = cellfun(@(x)x.Measurement(2),detBuffer);
r = sqrt(x.^2 + y.^2);
az = atan2d(y, x);Plot the distribution of clutter measurements as a function of its range. Notice that the distribution confirms the visual interpretation that the number of clutter measurements reduce with range.
% Plot histogram of range measurements figure; % Range bins to evaluate distribution rBins = 500:100:5000; % Number of false alarms per bin nFalse = histcounts(r, rBins); % Per scan nFalse = nFalse/nScans; % Plot distribution histogram('BinCounts',nFalse,'BinEdges',rBins); xlabel('Range bins (m)'); ylabel('Number of false alarms');
Next, you fit a function to define the clutter density of the radar as a function of range. Clutter density is defined as the number of false alarms per unit measurement volume. Clutter density is an important parameter used by multi-object trackers to compute data association between potential targets and measurements from the radar. In the context of radars reporting azimuth and range measurements, the volume of the sensor is defined in azimuth-range space with units of deg-meters. Clutter density is often assumed uniform or constant in the entire measurement space, however, it can also be defined as a function of measurement to inform the tracker about high or low clutter regions. The fitted function will be reused in the next section, when you define the sensor specification for tracking. In practice, this clutter map can be obtained from prior information about the scenario or can be estimated online by collecting sensor data over a few scans.
% Compute density of clutter binVolume = diff(rBins)*360; % Volume in azimuth-range = dR*dAz % Clutter density from the recorded data Kc = nFalse./binVolume; % Density = Number per unit volume % Range at middle of bins meanRange = rBins(1:end-1)/2 + rBins(2:end)/2; % Fit a function for Kc(r) = C0 + C1/r^5; C0 = mean(Kc(end-20:end)); C1 = meanRange(1)^5*(Kc(1) - C0); % Actual density plot(meanRange, Kc, 'LineWidth',2); hold on; % Fitted density Kcest = C0 + C1./meanRange.^5; plot(meanRange, Kcest, 'LineWidth',2); xlabel('Range (m)'); ylabel('Clutter intensity (K_c)'); legend('Recording','Fitted');
Set Up the Tracker
In this section, set up a task-oriented joint probabilistic data association tracker to track targets in the presence of non-uniform clutter in the radar scan. For more details on how to use the task-oriented tracker, refer to the documentation for the JIPDATracker.
Define target specification
Use the trackerTargetSpec to create a specification for the target used by the tracker.
targetSpec = trackerTargetSpec('custom');
disp(targetSpec); CustomTarget with properties:
StateTransitionModel: []
SurvivalModel: []
Choose a constant-velocity motion model to define the motion of the targets using the targetStateTransitionModel. Set the number of motion dimensions to 2 to estimate the position, and velocity of the boats in 2-D. Modify the velocity and acceleration variance to tune the constant-velocity model assuming maximum speed of 25 m/s and max acceleration of 5 m/s^2 for the boats.
targetSpec.StateTransitionModel = targetStateTransitionModel('constant-velocity'); targetSpec.StateTransitionModel.NumMotionDimensions = 2; targetSpec.StateTransitionModel.VelocityVariance = 25^2/3*eye(2); % maxSpeed^2/3 targetSpec.StateTransitionModel.AccelerationVariance = 3^2/3*eye(2); % maxAccel^2/3
To define the survival model, assume that the targets have uniform probability to survive within the field of view between radar observations.
targetSpec.SurvivalModel = targetSurvivalModel('uniform-survival-rate');
disp(targetSpec); CustomTarget with properties:
StateTransitionModel: [1×1 ConstantVelocityModel]
SurvivalModel: [1×1 UniformSurvivalRateModel]
Define sensor specification
Use the trackerSensorSpec function to create a specification for the radar used in this example. Set the maximum number of measurements to 75 based on the average number of false alarms produced per step.
% Create a sensor spec sensorSpec = trackerSensorSpec('custom'); sensorSpec.MaxNumMeasurements = 75; disp(numel(detBuffer)/nScans);
58.6000
disp(sensorSpec);
CustomSensor with properties:
MaxNumMeasurements: 75
MeasurementModel: []
DetectabilityModel: []
ClutterModel: []
BirthModel: []
UpdateModels: 0
Next, create all the models required by the sensor specification for this sensor.
Use a azimuth-range measurement model for the sensor. In the azimuth-range space, the measurement uncertainty or noise of the detections can be specified as a constant. Use the resolution of the radar to define the uncertainty in the measurement.
% By default, the radar is mounted at the origin. sensorSpec.MeasurementModel = sensorMeasurementModel('azimuth-range'); sensorSpec.MeasurementModel.AzimuthVariance = radar.AzimuthResolution^2; sensorSpec.MeasurementModel.RangeVariance = radar.RangeResolution^2;
To create the detectability model, assume that the radar can detect targets with uniform probability of 0.98 in the entire field of view.
sensorSpec.DetectabilityModel = sensorDetectabilityModel('uniform');
sensorSpec.DetectabilityModel.DetectionProbability = radar.DetectionProbability;To create the clutter model, use a non-uniform clutter model. This model allows to specify variable clutter density as a function of measurement. The ModelData allows to use any additional information that is required for computing the clutter density. For this example, use the coefficients of the fitted model from the Clutter Distribution section as ModelData to be used during computation of clutter density.
sensorSpec.ClutterModel = sensorClutterModel('nonuniform-poisson'); sensorSpec.ClutterModel.ClutterDensityFcn = @computeClutterDensity; sensorSpec.ClutterModel.ModelData = struct('Coefficients',[C0 C1]); function Kc = computeClutterDensity(z, modelData) % z is the measurement. % z = [azimuth;range]; % Kc is the clutter density r = z(2); Co = modelData.Coefficients(1); C1 = modelData.Coefficients(2); Kc = Co + C1/r^5; end
Similar to clutter density, the tracker also needs a birth density to define the likelihood of new targets appearing in the field of view. To account for a high clutter density inside the 1 Km radius and the assuming that new targets don't get detected first time inside that region, set the birth density low enough to avoid creating new tracks inside that region. This is required to maintain the computational and memory requirements of a JIPDA tracker.
% Create uniform birth model sensorSpec.BirthModel = sensorBirthModel('uniform-poisson'); % Detection probability of the sensor Pd = sensorSpec.DetectabilityModel.DetectionProbability; % Clutter density at r=1000 Kc1000 = C0 + C1/1000^5; % 1/10 of clutter density beta = 0.1*Kc1000; % Set the constant birth density sensorSpec.BirthModel.BirthDensity = beta; disp(sensorSpec);
CustomSensor with properties:
MaxNumMeasurements: 75
MeasurementModel: [1×1 AzimuthRangeModel]
DetectabilityModel: [1×1 UniformDetectabilityModel]
ClutterModel: [1×1 NonUniformPoissonModel]
BirthModel: [1×1 UniformPoissonModel]
UpdateModels: 0
Define the JIPDA tracker
In this section, configure the JIPDA tracker to integrate target and sensor specifications defined in the previous step. The tracker uses the modeling information from the target and sensor specifications for common tasks such as prediction and data association.
tracker = multiSensorTargetTracker(targetSpec, sensorSpec, 'jipda');
tracker.MaxMahalanobisDistance = 8;
tracker.ConfirmationExistenceProbability = 0.98;
disp(tracker) fusion.tracker.JIPDATracker with properties:
TargetSpecifications: {[1×1 CustomTarget]}
SensorSpecifications: {[1×1 CustomSensor]}
MaxMahalanobisDistance: 8
ConfirmationExistenceProbability: 0.9800
DeletionExistenceProbability: 0.1000
Understand sensor data format
Use the trackerSensorSpec function to understand the format of the data required by the tracker under the defined configuration. Notice that the tracker only needs the measurements as 2-element vector defined according the the azimuth-range measurement model. In addition, it also needs the timestamp of the observations.
sensorData = dataFormat(sensorSpec); disp(sensorData)
MeasurementTime: 0
Measurements: [2×75 double]
Run Scenario
In this section, run the scenario with three boats added to the scenario. Simulate radar measurements in the presence of heavy sea clutter, and send those measurements to the tracker to track boats. The tracker is responsible to filter out clutter based on the clutter density and birth density provided by the sensor specification.
% Restart the scenario and add targets now restart(scene); addTargets(scene); % Recreate display display = createDisplay(scene); % Initialize detection buffer for each scan detBuffer = {}; % Run the scenario while advance(scene) % Detections in current field of view [dets, config] = detect(scene); % Run the tracker when the sensor completes one scan if config.IsScanDone % Assemble sensor data from detection buffer sensorData = assembleData(sensorData, detBuffer); % Update tracker tracks = tracker(sensorData); % Display updateDisplay(display, scene, detBuffer, tracks); % Empty buffer detBuffer = {}; else % Add to the buffer detBuffer = [detBuffer;dets]; updateDisplay(display, scene, detBuffer); end if abs(scene.SimulationTime - 119) < 0.05 snap = zoomAndSnap(display); end end
Results
In the image above, you can visualize the results of the tracker after the scenario has ended. Notice that the tracker is able to maintain a track on all three targets. In addition, the tracker was able to filter out all the sea clutter and only confirmed one false track during the scenario. Notice, when the boat reaches near the clutter region, the data association between measurements and target estimate is significantly more ambiguous than outside the clutter region. The JIPDA tracker handles this ambiguity by considering multiple possible data associations and representing the uncertainty in data association in the estimate of the target state.
Notice in the image below zoomed in at 1 kilometer radius that the track had a higher uncertainty in the position estimate in the clutter region, represented by the orange ellipse around the track.
figure; imshow(snap.cdata);
Summary
In this example, you learned how to simulate radar clutter from sea surfaces in rough sea conditions using a statistical radar model and observed the challenges of tracking under such conditions. Additionally, you learned how to configure a JIPDA tracker to process the radar data and use appropriate clutter modeling to assist the tracker in filtering out false alarms and maintaining tracks under ambiguous situations close to heavy clutter.
Supporting Functions
assembleData
function sensorData = assembleData(sensorData, detBuffer) % This function converts the objectDetection cell array to sensor data % format required by the tracker. x = cellfun(@(x)x.Measurement(1),detBuffer); y = cellfun(@(x)x.Measurement(2),detBuffer); az = atan2d(y,x); r = sqrt(x.^2 + y.^2); sensorData.MeasurementTime = detBuffer{end}.Time; % Time of report after the scan is complete sensorData.Measurements = [az';r']; % Measurements end
addTargets
function addTargets(scene) % Create first midsize boat with RCs = 10 Boat = platform(scene); Boat.Trajectory = waypointTrajectory( ... [1502 1951 0;1285 1480 0;1126 973 0;1126 350 0;981 -390 0;-454 -1390 0], ... 'GroundSpeed', 20*ones(6,1)); Boat.Signatures{1} = rcsSignature(Azimuth=[-180 180],Elevation=[-90;90],Pattern=10*ones(2)); % Create second large boat with RCS = 25 Boat = platform(scene); Boat.Trajectory = waypointTrajectory( ... [-2536 4062 0;-1573 4312 0;1324 4167 0;1831 3636 0], ... 'GroundSpeed', [20;20;20;20]); Boat.Signatures{1} = rcsSignature(Azimuth=[-180 180],Elevation=[-90;90],Pattern=25*ones(2)); % Create third small boat that reaches close with RCS = 0 Boat = platform(scene,'ClassID',4); Boat.Trajectory = waypointTrajectory( ... [1699 -1364 0;724 -1051 0;-354 -616 0;-2408 -616 0;-3817 0 0], ... 'GroundSpeed', 20*ones(5,1)); Boat.Signatures{1} = rcsSignature(Azimuth=[-180 180],Elevation=[-90;90],Pattern= 0*ones(2)); end
createDisplay
function display = createDisplay(scene) % This function creates the display for visualizing tracks, detections, and % ground truth. % Color order for plots clrs = [1.0000 1.0000 0.0667 0.0745 0.6235 1.0000 1.0000 0.4118 0.1608 0.3922 0.8314 0.0745 0.7176 0.2745 1.0000 0.0588 1.0000 1.0000 1.0000 0.0745 0.6510]; % Create theater plot f = figure('Units','normalized',... 'Position',[0.1 0.1 0.6 0.8],... 'Visible','on'); ax = axes(f); display.TheaterPlot = theaterPlot('XLimits',[-5e3 5e3],... 'YLimits',[-5e3 5e3],... 'ZLimits',[-5e3 5e3],... 'Parent',ax); % Create a detection plotter display.DetectionPlotter = detectionPlotter(display.TheaterPlot,... 'DisplayName','Radar Detections',... 'MarkerFaceColor',clrs(4,:),... 'MarkerEdgeColor',clrs(1,:),... 'Marker','o'); % Create a platform plotter for targets display.PlatformPlotter = platformPlotter(display.TheaterPlot,... 'DisplayName','Ground Truth',... 'MarkerFaceColor',clrs(2,:)); % Create a trajectory plotter for targets and plot their trajectories trp = trajectoryPlotter(display.TheaterPlot,... 'DisplayName','Trajectories',... 'Color',clrs(2,:),... 'LineStyle','-',... 'LineWidth',2); % Plot trajectories during construction t = linspace(0,scene.StopTime,100); pos = cell(numel(scene.Platforms)-1,1); for i = 1:numel(pos) pos{i} = lookupPose(scene.Platforms{i+1}.Trajectory,t); end trp.plotTrajectory(pos); % Create a coverage plotter display.CoveragePlotter = coveragePlotter(display.TheaterPlot,... 'DisplayName','',... 'Color',clrs(4,:),... 'Alpha',[0.1 0]); % Create a track plotter display.TrackPlotter = trackPlotter(display.TheaterPlot,... 'DisplayName','Tracks',... 'MarkerFaceColor',clrs(3,:),... 'MarkerEdgeColor',clrs(3,:),... 'ConnectHistory','on'); ax = display.TheaterPlot.Parent; hold on; % Plot range lines on the PPI az = linspace(0,360,100)'; r = 1000:1000:5000; X = zeros(0,1); Y = zeros(0,1); for i = 1:numel(r) x = r(i).*cosd(az); y = r(i).*sind(az); X = [X;nan;x]; %#ok<AGROW> Y = [Y;nan;y]; %#ok<AGROW> end plot(ax,X,Y,'g-'); % Use dark theme colors ax = display.TheaterPlot.Parent; ax.Parent.Color = 0.157*[1 1 1]; ax.Color = [0 0 0]; grid(ax,'on'); ax.GridColor = 0.68*[1 1 1]; ax.GridAlpha = 0.4; axis(ax,'equal'); ax.XColor = [1 1 1]*0.68; ax.YColor = [1 1 1]*0.68; ax.ZColor = [1 1 1]*0.68; l = legend(ax); l.Orientation = 'horizontal'; l.Location = 'northoutside'; end function updateDisplay(display, scene, detections, tracks) % This function updates the display for plotting ground truth, detections, % and tracks. % Plot ground truth poses = platformPoses(scene); display.PlatformPlotter.plotPlatform(vertcat(poses.Position)); % Plot coverage cvg = coverageConfig(scene); display.CoveragePlotter.plotCoverage(cvg); % Plot detections if nargin > 2 && ~isempty(detections) x = cellfun(@(x)x.Measurement(1),detections); y = cellfun(@(x)x.Measurement(2),detections); z = zeros(size(x)); display.DetectionPlotter.plotDetection([x y z]); end % Plot tracks if nargin > 3 && ~isempty(tracks) for tr = 1:numel(tracks) tracks(tr).State = [tracks(tr).State;0;0]; tracks(tr).StateCovariance = blkdiag(tracks(tr).StateCovariance,1,1); end [trkPos,trkPosCov] = getTrackPositions(tracks,'constvel'); display.TrackPlotter.plotTrack(trkPos,trkPosCov,string([tracks.TrackID])); end drawnow; end function snap = zoomAndSnap(display) display.TheaterPlot.Parent.XLim = [-1e3 1e3]; display.TheaterPlot.Parent.YLim = [-1e3 1e3]; snap = getframe(display.TheaterPlot.Parent.Parent); display.TheaterPlot.Parent.XLim = [-5e3 5e3]; display.TheaterPlot.Parent.YLim = [-5e3 5e3]; end
