Monostatic Radar Sensor POD

3 次查看（过去 30 天）
Chris Raper
Chris Raper 2021-6-28
编辑： Chris Raper 2021-7-2
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)))
  1 个评论
Chris Raper
Chris Raper 2021-7-2
编辑：Chris Raper 2021-7-2
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Tracking and Sensor Fusion 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by