Monostatic Radar Sensor POD

I'm trying to understand the probability of detection profiles associated with monostatic radar stations in the Tracking & Sensor Fusion toolbox. I've created a scenario where there is one ARSR-4 radar station located at the Statue of Liberty and several 737s start at the same location, all slowly flying westward over 800 miles. The seven aircraft maintain a constant altitude, anywhere from 0 to 30,000ft.

Below is the resulting figure which indicates how likely the radar station is to detect the aircraft at those ranges & altitudes. There were 86,400 radar updates in this simulation, so sampling error should be negligible. Can anyone explain the odd characteristics seen in these plots?

• Why are there intermitent peaks as range increases for every altitude except 5,000ft?
• Why do we see POD increasing in most plots beyond a range oif 200 miles?
• Why don't we see POD taper to zero with large ranges? For most altitudes, POD seems to be truncated rather than taper to zero.
• Why is 10,000ft worse than any other altitude, especially when 5,000ft is the best?

Here's my code in case you'd like to reproduce the plot. It took more than 2 hours of runtime on my machine.

%% 
% Wrangle flight plan, airport & radar data files

clear

%% 
% Add all radar stations to simulation

s = rng; % Save current random generator state
rng(2020); % Set random seed for predictable results

% Create scenario
scene = trackingScenario('IsEarthCentered',true,'UpdateRate',1);

% Model an ARSR-4 radar
updaterate = 1;
fov = [360;30.001];
elAccuracy = atan2d(0.9,463); % 900m accuracy @ max range
elBiasFraction = 0.1;

arsr4 = monostaticRadarSensor(1,'UpdateRate',updaterate,...
    'FieldOfView',fov,...,
    'HasElevation',true,...
    'ScanMode','Mechanical',...
    'MechanicalScanLimits',[0 360; -30 0],...
    'HasINS',true,...
    'HasRangeRate',true,...
    'HasFalseAlarms',false,...
    'ReferenceRange',463000,... = 288 miles
    'ReferenceRCS',0,...
    'AzimuthResolution',1.4,...
    'AzimuthBiasFraction',0.176/1.4,...
    'ElevationResolution',elAccuracy/elBiasFraction,...
    'ElevationBiasFraction',elBiasFraction,...
    'RangeResolution', 323,...
    'RangeBiasFraction',116/323,... Accuracy / Resolution
    'RangeRateResolution',100,...
    'DetectionCoordinates','Scenario');

% Add ARSR-4 radars at each site
radar = clone(arsr4);
radar.SensorIndex = 1;
start_ll = [40+41/60+21.7/3600,-74-02/60-42.5/3600]; % statue of liberty
platform(scene,'Position',[start_ll,0],...
    'Signatures',rcsSignature('Pattern',-50),'Sensors',{radar});
%% 
% Define all flights
m_per_ft = 0.3048;
load('737rcs.mat');
n_active_flights = 0;
platform_id_ls = [];
callsigns_ls = {};
icao_ls = {};
airplane_ls = [];
time_max = 24*60*60;
end_ll = [40+41/60+21.7/3600,-90]; % due east 834 miles
for alt=0:5000:30000
    alt = alt*m_per_ft;
    flight_route = geoTrajectory([[start_ll, alt];[end_ll, alt]],[0, time_max]);
    airplane = platform(scene,'Trajectory',flight_route);
    airplane.Signatures{1} = boeing737rcs;
    airplane_ls = [airplane_ls; airplane]; %#ok<AGROW>
end
%% 
% Run the simulation

tic;

% Declare loop variables
radar_time_vec = [];
radar_target_id_vec = [];
radar_sensor_id_vec = [];
radar_lat_vec = [];
radar_lon_vec = [];
radar_alt_vec = [];
radar_az_vec = [];
radar_rng_vec = [];
radar_rng_rate_vec = [];
radar_assign_trk_id = [];
det_radar = [];

arsrupdatetime  = 1;
sim_start = datetime('now');
progress = [];
pred_runtime = [];
scene_spherical = clone(scene);
scene_spherical.Platforms{1,1}.Sensors{1,1}.DetectionCoordinates = 'Sensor spherical';
f = waitbar(0,'Simulation progress');
while advance(scene)
    advance(scene_spherical);
    time = scene.SimulationTime;
    
    % Generate radar detections at the defined rate
    if mod(time,arsrupdatetime) == 0 
        % Generate synthetic radar detections
        rng(2020); % Set random seed for predictable results
        dets = detect(scene);
        n_radar_dets = length(dets);
        rng(2020); % Set random seed for predictable results
        dets_spherical = detect(scene_spherical);
        % dets = removeBelowGround(dets);
        det_radar = [det_radar; dets]; %#ok<AGROW>
        if ~isempty(dets)
            for i=1:length(dets)
                position = dets{i,1}.Measurement([1 2 3]);
                position_spherical = dets_spherical{i,1}.Measurement([1 2 3]);
                velocity = dets{i,1}.Measurement([4 5 6]);
                velocity_spherical = dets_spherical{i,1}.Measurement([4]);
                sensor_position = scene.Platforms{1,2}.Position;
                lla = fusion.internal.frames.ecef2lla(position');
                
                radar_time_vec = [radar_time_vec; time]; %#ok<AGROW>
                radar_target_id_vec = [radar_target_id_vec; dets{i,1}.ObjectAttributes{1,1}.TargetIndex]; %#ok<AGROW>
                radar_sensor_id_vec = [radar_sensor_id_vec; dets{i,1}.SensorIndex]; %#ok<AGROW>
                radar_lat_vec = [radar_lat_vec; lla(1)]; %#ok<AGROW>
                radar_lon_vec = [radar_lon_vec; lla(2)]; %#ok<AGROW>
                radar_alt_vec = [radar_alt_vec; lla(3)]; %#ok<AGROW>
                
                radar_az_vec = [radar_az_vec; position_spherical(1)]; %#ok<AGROW>
                radar_rng_vec = [radar_rng_vec; position_spherical(3)]; %#ok<AGROW>
                radar_rng_rate_vec = [radar_rng_rate_vec; velocity_spherical(1)]; %#ok<AGROW>
            end
        end
    end
    
    progress = [progress; [toc  time_max-time]];
    if (mod(time,36) == 0) && (time > 0)
        p = polyfit(progress(:,1),progress(:,2),2);
        r = roots(p);
        r_real = real(r);
        r_pos_real = r_real(r_real > 0);
        if isempty(r_pos_real)
            waitbar(time/time_max,f,sprintf('Simulation progress is %.0f%% (%d/%d)\nMinutes remaining: unknown',100*time/time_max,time,time_max))
        else
            pred_runtime = [pred_runtime; min(r_pos_real)];
            etc = min(r_pos_real)/60 - toc/60;
            waitbar(time/time_max,f,sprintf('Simulation progress is %.0f%% (%d/%d)\nMinutes remaining: %.1f',100*time/time_max,time,time_max,etc))
        end
    end
end
radar_df = timetable(sim_start + seconds(radar_time_vec + .001),radar_target_id_vec,radar_sensor_id_vec,radar_lat_vec,...
    radar_lon_vec,radar_alt_vec,radar_az_vec,radar_rng_vec,radar_rng_rate_vec,...
    'VariableNames',{'TargetID','SensorID','Lat','Lon','Alt','Azimuth','Range','RangeRate'});
head(radar_df)
waitbar(1,f,sprintf('Simulation took %.1f minutes to complete',toc/60))
%%
% Visualize results
radar_df.Loc = [radar_df.Lon radar_df.Lat radar_df.Alt];

alt_rng = my_range(radar_df.Alt);
alt_rng(1) = 0;
tp = theaterPlot('XLim',my_range(radar_df.Lon),'YLim',my_range(radar_df.Lat),'ZLim',alt_rng);
% radarCovPlotter = coveragePlotter(tp,'DisplayName','Radar Coverage');
% covcon = coverageConfig(scene);
% covcon.Position = fusion.internal.frames.ecef2lla(covcon.Position);
% plotCoverage(radarCovPlotter,covcon);
radarPlotter = detectionPlotter(tp,'DisplayName','Radar Detections','MarkerEdgeColor','blue','Marker','o');
plotDetection(radarPlotter, radar_df.Loc);
grid on
%%
% Visualize ground truth
gl = helperGlobeViewer;
setCamera(gl,[28.9176  -95.3388  5.8e5],[0 -30 10]);

% Show radar sites
plotPlatform(gl,scene.Platforms(1:5),'d');

% Show radar coverage
covcon = coverageConfig(scene);
plotCoverage(gl,covcon);

% Show flight routes
restart(scene);
for i=2:8
    plotTrajectory(gl,scene.Platforms{1,i});
end

% Show some detections
plotDetection(gl,[det_radar(1:500)]);

m_per_mi = 1609;
target_ls = unique(radar_df.TargetID);
radar_df.Range = radar_df.Range/m_per_mi;
range_max = norm(fusion.internal.frames.lla2ecef([start_ll,0]) -...
                    fusion.internal.frames.lla2ecef([end_ll,30000*m_per_ft]))/m_per_mi;
edge_ls = 0:range_max/50:range_max;
n_expected_obs = 24*60*60/(length(edge_ls) - 1);
figure
for i=1:length(target_ls)
    subplot(2,4,i)
    N = histcounts(table2array(radar_df(radar_df.TargetID == target_ls(i), "Range")), edge_ls);
    histogram('BinEdges',edge_ls,'BinCounts',N/n_expected_obs)
    title(sprintf("%dft alt",(i-1)*5000))
    xlabel("Range (miles)")
    ylabel("POD")
end

figure
plot(progress(:,1),progress(:,2))
hold on
plot(progress(1,1):progress(end,1), polyval(polyfit(progress(:,1),progress(:,2),2),progress(1,1):progress(end,1)))

[Chris Raper]

It seems the strange qualities of the plot above were caused by falsely attributing detections to various altitudes and/or ranges for reasons that aren't clear to me. I modified my script to simulate 1 airplane at a time, increasing the altitude with each successive simulation. The following figure gives my new results, which are much cleaner.

However, I still don't fully understand the shape of these POD curves. An ARSR-4 is supposed to have a POD of 0.9 at a range of 288 miles. Why is the simulation not showing that? I suppose it is because the RCS for the 737 is much larger than the reference RCS? If so, could someone suggest where I could find radar cross section data for smaller aircraft besides the 737?

Here's my updated code:

s = rng; % Save current random generator state
rng(2020); % Set random seed for predictable results

% Create scenario
scene = trackingScenario('IsEarthCentered',true,'UpdateRate',1);

% Model an ARSR-4 radar
updaterate = 5;
fov = [360;30.001];
elAccuracy = atan2d(0.9,463); % 900m accuracy @ max range
elBiasFraction = 0.1;

arsr4 = monostaticRadarSensor(1,'UpdateRate',updaterate,...
    'FieldOfView',fov,...,
    'HasElevation',true,...
    'ScanMode','Mechanical',...
    'MechanicalScanLimits',[0 360; -30 0],...
    'HasINS',true,...
    'HasRangeRate',true,...
    'HasFalseAlarms',false,...
    'ReferenceRange',463000,... = 288 miles
    'ReferenceRCS',0,...
    'AzimuthResolution',1.4,...
    'AzimuthBiasFraction',0.176/1.4,...
    'ElevationResolution',elAccuracy/elBiasFraction,...
    'ElevationBiasFraction',elBiasFraction,...
    'RangeResolution', 323,...
    'RangeBiasFraction',116/323,... Accuracy / Resolution
    'RangeRateResolution',100,...
    'DetectionCoordinates','Scenario');

% Add ARSR-4 radars at each site
radar = clone(arsr4);
radar.SensorIndex = 1;
start_ll = [40+41/60+21.7/3600,-74-02/60-42.5/3600]; % statue of liberty
platform(scene,'Position',[start_ll,0],...
    'Signatures',rcsSignature('Pattern',-50),'Sensors',{radar});
%% 
% Define all flights
m_per_ft = 0.3048;
load('737rcs.mat');
platform_id_ls = [];
airplane_ls = [];
time_max = 24*60*60;
end_ll = [40+41/60+21.7/3600,-90]; % due east 834 miles
alt_ft_ls = 0:5000:40000;
n_active_flights = length(alt_ft_ls);
flight_route = geoTrajectory([[start_ll, 0];[end_ll, 0]],[0, time_max]);
airplane = platform(scene,'Trajectory',flight_route);
airplane.Signatures{1} = boeing737rcs;
airplane_ls = [airplane_ls; airplane]; %#ok<AGROW>
platform_id_ls = [platform_id_ls; scene.Platforms{end}.PlatformID]; %#ok<AGROW>
%% 
% Run the simulation

tic;

% Declare loop variables
radar_time_vec = [];
radar_target_id_vec = [];
radar_sensor_id_vec = [];
radar_lat_vec = [];
radar_lon_vec = [];
radar_alt_vec = [];
radar_az_vec = [];
radar_rng_vec = [];
radar_rng_rate_vec = [];
radar_assign_trk_id = [];
det_radar = [];

arsrupdatetime  = 5;
sim_start = datetime('now');
progress = [];
pred_runtime = [];

for j=1:n_active_flights
    restart(scene)
    alt = alt_ft_ls(j)*m_per_ft;
    scene.Platforms{1, 2}.Trajectory = geoTrajectory([[start_ll, alt];[end_ll, alt]],[0, time_max]);
    scene_spherical = clone(scene);
    scene_spherical.Platforms{1,1}.Sensors{1,1}.DetectionCoordinates = 'Sensor spherical';
    f = waitbar(0,'Simulation progress');
    while advance(scene)
        advance(scene_spherical);
        time = scene.SimulationTime;

        % Generate radar detections at the defined rate
        if mod(time,arsrupdatetime) == 0 
            % Generate synthetic radar detections
            rng(2020); % Set random seed for predictable results
            dets = detect(scene);
            n_radar_dets = length(dets);
            rng(2020); % Set random seed for predictable results
            dets_spherical = detect(scene_spherical);
            % dets = removeBelowGround(dets);
            det_radar = [det_radar; dets]; %#ok<AGROW>
            if ~isempty(dets)
                for i=1:length(dets)
                    position = dets{i,1}.Measurement([1 2 3]);
                    position_spherical = dets_spherical{i,1}.Measurement([1 2 3]);
                    velocity = dets{i,1}.Measurement([4 5 6]);
                    velocity_spherical = dets_spherical{i,1}.Measurement([4]);
                    sensor_position = scene.Platforms{1,2}.Position;
                    lla = fusion.internal.frames.ecef2lla(position');

                    radar_time_vec = [radar_time_vec; time]; %#ok<AGROW>
                    radar_target_id_vec = [radar_target_id_vec; j]; %#ok<AGROW>
                    radar_sensor_id_vec = [radar_sensor_id_vec; dets{i,1}.SensorIndex]; %#ok<AGROW>
                    radar_lat_vec = [radar_lat_vec; lla(1)]; %#ok<AGROW>
                    radar_lon_vec = [radar_lon_vec; lla(2)]; %#ok<AGROW>
                    radar_alt_vec = [radar_alt_vec; lla(3)]; %#ok<AGROW>

                    radar_az_vec = [radar_az_vec; position_spherical(1)]; %#ok<AGROW>
                    radar_rng_vec = [radar_rng_vec; position_spherical(3)]; %#ok<AGROW>
                    radar_rng_rate_vec = [radar_rng_rate_vec; velocity_spherical(1)]; %#ok<AGROW>
                end
            end
        end

        steps_remaining = (n_active_flights - j + 1)*time_max - time;
        progress = [progress; [toc  steps_remaining]];
        if (mod(time,120) == 0) && (time > 0)
            p = polyfit(progress(:,1),progress(:,2),2);
            r = roots(p);
            r_real = real(r);
            r_pos_real = r_real(r_real > 0);
            completion = ((j-1)*time_max+time)/(n_active_flights*time_max);
            if isempty(r_pos_real)
                waitbar(completion,f,sprintf('Simulation progress is %.0f%% (%d/%d)\nMinutes remaining: unknown',100*time/time_max,time,time_max))
            else
                pred_runtime = [pred_runtime; min(r_pos_real)];
                etc = min(r_pos_real)/60 - toc/60;
                waitbar(completion,f,sprintf('Simulation progress is %.0f%% (%d/%d)\nMinutes remaining: %.1f',100*time/time_max,time,time_max,etc))
            end
        end
    end
end
radar_df = timetable(sim_start + seconds(radar_time_vec + .001),radar_target_id_vec,radar_sensor_id_vec,radar_lat_vec,...
    radar_lon_vec,radar_alt_vec,radar_az_vec,radar_rng_vec,radar_rng_rate_vec,...
    'VariableNames',{'TargetID','SensorID','Lat','Lon','Alt','Azimuth','Range','RangeRate'});
m_per_mi = 1609;
radar_df.Range = radar_df.Range/m_per_mi;
% head(radar_df)
waitbar(1,f,sprintf('Simulation took %.1f minutes to complete',toc/60))
%%
radar_df.Loc = [radar_df.Lon radar_df.Lat radar_df.Alt];
target_ls = unique(radar_df.TargetID);
range_max = norm(fusion.internal.frames.lla2ecef([start_ll,0]) -...
                    fusion.internal.frames.lla2ecef([end_ll,30000*m_per_ft]))/m_per_mi;
edge_ls = 0:range_max/50:range_max;
range_coords = (edge_ls(1:end-1) + edge_ls(2:end))/2;
n_expected_obs = 24*60*60/(length(edge_ls) - 1)/arsrupdatetime;
color_ls = colormap('turbo');
color_ls = color_ls(1:floor(length(color_ls)/length(target_ls)):length(color_ls),:);
figure
hold on
legend()
for i=1:length(target_ls)
    % subplot(2,5,i)
    N = histcounts(table2array(radar_df(radar_df.TargetID == target_ls(i), "Range")), edge_ls);
    % histogram('BinEdges',edge_ls,'BinCounts',N/n_expected_obs)
    plot(range_coords,N/n_expected_obs,'-','DisplayName',sprintf("%dft alt",alt_ft_ls(i)),'Color',color_ls(i,:))
end
title("POD of ARSR-4 for 737 based on Range & Altitude")
xlabel("Range (miles)")
ylabel("POD")

figure
plot(progress(:,1),progress(:,2))
hold on
plot(progress(1,1):progress(end,1), polyval(polyfit(progress(:,1),progress(:,2),2),progress(1,1):progress(end,1)))

%%
% Support functions
function x_rng = my_range(x)
    x_min = min(x);
    x_max = max(x);
    x_rng = [x_min x_max];
end