Why do my radar detections give me the wrong position even with extremely accurate resolutions and no noise?

I am working with the sensor fusion and tracking toolbox to simulate radar detections from multiple radars that I am trying to correlate to see if it is the same target. I am using two radars to detect a single target sitting in the air between them. These radars are broken into four quadrants that each view 90 degrees, with no scanning (essentially fov is 360 degrees). When I run this scenario though, my radar detections in the X position are always off by 1.38 m in the direction the radar is looking (when the target is 5000 m away from the radars in the X direction). When I repeated this, but put the target 5000 m away in the Y direction, it was also off by 1.38 m in the Y direction, regardless of how far the target was in the X direction. Additionally, this error seems to scale with distance as at 10000 m away, the targets detected position was off by 2.75 m.

I thought this might have to do with inaccurate radar settings (resolutions/biases) and noise, but I have changed those so that the radar would basically be perfectly accurate with no noise and the error still persists. I don't mind the detections being off slightly, but having these errors despite accurate radar settings messes up my correlation, since it expects extremely accurate measurements due to the extremely accurate radar settings. Why am I not getting exact radar measurements, even though all the settings are extremely accurate and noise is turned off?

Here is my attached code. It also creates a text file at the end of all the radar detections of the target. The text file includes the time of the detection, which radar took the detection (system ID), which target it is (there is only one target in this case), the position, velocity and some measurement noise data (0 because the radar settings are extremely accurate). In this scenario, the target is sitting still at [5000 300 -1000]. Radar 1 which is positioned at [0 0 -4] always says the targets position is [5001.38 300.8 -1000.75], while Radar 2 which is positioned at [10000 0 -4] always sees the target at [4998.62 300.8 -1000.75]. Since the two radars are detecting in opposite X-directions that is why their X positions are add the 1.38 m in different directions, but they detect in the same Y-directions, so that coordinate is the same. These should be returning the actual position of [5000 300 -1000].

clear all;
close all;

%% Radar Creation
scenario = trackingScenario;
radarnum = 2;

%Platform for Radar 1
sensor1 = platform(scenario);

%Platform for Radar 2.1 
Secondsensor1 = platform(scenario);

% Radar Characteristics

updaterate = 4; 
fov = [90;88];

radar1 = monostaticRadarSensor(1, ...
    'Raster', ...
    'ScanMode', 'Mechanical', ...
    'UpdateRate', updaterate, ...           % Hz
    'FieldOfView', fov, ...                 % [az;el] deg
    'AzimuthResolution', 0.0000000001, ...        % deg
    'ElevationResolution', 0.000000001, ...
    'ReferenceRange', 375e3, ...            % m
    'ReferenceRCS', 0, ...                  % dBsm
    'RangeResolution', 0.00000000001, ...            % m
    'HasINS', true, ...
    'DetectionCoordinates', 'Scenario', ...
    'HasFalseAlarms',false, ...
    'HasNoise',false, ...
    'MaxMechanicalScanRate', [0, 0], ...
    'DetectionProbability',1, ...
    'AzimuthBiasFraction', 0, ...
    'ElevationBiasFraction', 0, ...
    'RangeBiasFraction',0);

Secondradar1 = monostaticRadarSensor(2, ...
    'Raster', ...
    'ScanMode', 'Mechanical', ...
    'UpdateRate', 5, ...           % Hz
    'FieldOfView', fov, ...                 % [az;el] deg
    'AzimuthResolution', 0.00000000001, ...        % deg
    'ElevationResolution', 0.000000001, ...
    'ReferenceRange', 375e3, ...            % m
    'ReferenceRCS', 0, ...                  % dBsm
    'RangeResolution', 0.00000000001, ...            % m
    'HasINS', true, ...
    'DetectionCoordinates', 'Scenario', ...
    'HasFalseAlarms',false, ...
    'HasNoise',false, ...
    'MaxMechanicalScanRate', [0, 0], ...
    'DetectionProbability',1, ...
    'AzimuthBiasFraction', 0, ...
    'ElevationBiasFraction', 0, ...
    'RangeBiasFraction',0);


%Get Elevation Measurement for Each Radar
radar1.HasElevation = true;
Secondradar1.HasElevation = true;

%Set Scan Limits for Each Radar
radar1.MechanicalScanLimits =[0 90;0 50]; % deg
Secondradar1.MechanicalScanLimits=[0 90; 0 50];

%Clone Radar one to create each quadrant
radar2 = clone(radar1);
radar3 = clone(radar1);
radar4 = clone(radar1);

%Clone Radar two to create each quadrant
Secondradar2 = clone(Secondradar1);
Secondradar3 = clone(Secondradar1);
Secondradar4 = clone(Secondradar1);

% Mount radar at orgin
radar1.MountingLocation = [0.1 0 -4];
radar2.MountingLocation = [0 0.1 -4];
radar3.MountingLocation = [-0.1 0 -4];
radar4.MountingLocation = [0 -0.1 -4];
radarpos(1,:) = [radar2.MountingLocation(1),radar3.MountingLocation(2),radar2.MountingLocation(3)];

% Mount second radar nearby
Secondradar1.MountingLocation = [10000.1 0 -4];
Secondradar2.MountingLocation = [10000 0.1 -4];
Secondradar3.MountingLocation = [9999.9 0 -4];
Secondradar4.MountingLocation = [10000 -0.1 -4];
radarpos(2,:) = [Secondradar2.MountingLocation(1),Secondradar3.MountingLocation(2),Secondradar2.MountingLocation(3)];

% Make sure each part of radar one is facing a different direction
radar1.MountingAngles = [0 0 0];
radar2.MountingAngles = [90 0 0];
radar3.MountingAngles = [180 0 0];
radar4.MountingAngles = [270 0 0];

Secondradar1.MountingAngles = [0 0 0];
Secondradar2.MountingAngles = [90 0 0];
Secondradar3.MountingAngles = [180 0 0];
Secondradar4.MountingAngles = [270 0 0];

%Attach Sensors to Radars - Only change if you don't want four quadrants
sensor1.Sensors = {radar1,radar2,radar3,radar4};
Secondsensor1.Sensors = {Secondradar1,Secondradar2,Secondradar3,Secondradar4};
  
%% Target Creation

ht = 1000;
wp = [5e3, 300, -ht;5e3, 300, -ht; 5e3, 300, -ht]; %Sit in the same spot
traj = waypointTrajectory('Waypoints',wp,'TimeOfArrival',[0 10 20]);    
platform(scenario,'Trajectory', traj);

%% Run Sim
% Set simulation to advance at the update rate of the radar
scenario.UpdateRate = 20;

velos = [];
tracks = [];
detLog = {};
timeLog = [];

while advance(scenario)
    
    % Current simulation time
    simTime = scenario.SimulationTime;
    dets = {};
    
    % Target poses in the Radar's coordinate frame
    targets1 = targetPoses(Secondsensor1);
 
    [dets,configs] = detect(scenario);
    
    if(isempty(dets) ~= 1)
       for i = 1:size(dets,1)
           dist = sqrt(sum((dets{i,1}.Measurement' - radarpos(dets{i,1}.SensorIndex,:)) .^2));
           %Set Radars minimum and maximum distances
           if(abs(dets{i,1}.Measurement(3)) > 20 && dist + randi([-3000 3000]) < 14000 && dist + randi([-100 100]) > 1000)
               %Record Data for printing later
               detLog = [detLog;dets(i,1)];
               velos = [velos;targets1(dets{i,1}.ObjectAttributes{1,1}.TargetIndex - 1).Velocity];
           end
       end
    end
end


%Create simulation data text files

fileID = fopen('Testing_Scenario.txt','w'); %Name of file created
sd = 0.0; %Standard deviation of Velocity Noise in m
normal_noise = sd*randn(size(velos,1),size(velos,2)); %Create Noise for Velocity
velosNoise = velos + normal_noise; %Apply Noise to Velocity
format long g
start_time = detLog{1,1}.Time * 1000;

%Create File
for i = 1:size(detLog,1)
fprintf(fileID,'Message 2 \n',1);
if(size(detLog{i,1}.MeasurementNoise,1) == 3)
   fprintf(fileID,'Message Length 60 \n');
else
    fprintf(fileID,'Message Length 120 \n');
end
fprintf(fileID,'Time %.0f \n',detLog{i,1}.Time * 1000 - start_time);
fprintf(fileID,'System ID %i \n', detLog{i,1}.SensorIndex);
fprintf(fileID,'Target ID %i \n', detLog{i,1}.ObjectAttributes{1,1}.TargetIndex - radarnum);
fprintf(fileID,'Track Classification 205 \n');
fprintf(fileID,'Track ECEF X %.2f \n', round(detLog{i,1}.Measurement(1),2));
fprintf(fileID,'Track ECEF Y %.2f \n', detLog{i,1}.Measurement(2));
fprintf(fileID,'Track ECEF Z %.2f \n', detLog{i,1}.Measurement(3));
fprintf(fileID,'Track ECEF Vx %.2f\n', velosNoise(i,1));
fprintf(fileID,'Track ECEF Vy %.2f\n', velosNoise(i,2));
fprintf(fileID,'Track ECEF Vz %.2f\n', velosNoise(i,3));
if(size(detLog{i,1}.MeasurementNoise,1) == 3)
    fprintf(fileID,'P Elements 6 \n');
    fprintf(fileID,'P Matrix (1,1) %f\n', detLog{i,1}.MeasurementNoise(1,1));
    fprintf(fileID,'P Matrix (1,2) %f\n', detLog{i,1}.MeasurementNoise(1,2));
    fprintf(fileID,'P Matrix (1,3) %f\n', detLog{i,1}.MeasurementNoise(1,3));
    fprintf(fileID,'P Matrix (2,2) %f\n', detLog{i,1}.MeasurementNoise(2,2));
    fprintf(fileID,'P Matrix (2,3) %f\n', detLog{i,1}.MeasurementNoise(2,3));
    fprintf(fileID,'P Matrix (3,3) %f\n', detLog{i,1}.MeasurementNoise(3,3));
else
    fprintf(fileID,'P Elements 21 \n');
    fprintf(fileID,'P Matrix (1,1) %f\n', detLog{i,1}.MeasurementNoise(1,1));
    fprintf(fileID,'P Matrix (1,2) %f\n', detLog{i,1}.MeasurementNoise(1,2));
    fprintf(fileID,'P Matrix (1,3) %f\n', detLog{i,1}.MeasurementNoise(1,3));
    fprintf(fileID,'P Matrix (1,4) %f\n', detLog{i,1}.MeasurementNoise(1,4));
    fprintf(fileID,'P Matrix (1,5) %f\n', detLog{i,1}.MeasurementNoise(1,5));
    fprintf(fileID,'P Matrix (1,6) %f\n', detLog{i,1}.MeasurementNoise(1,6));
    fprintf(fileID,'P Matrix (2,2) %f\n', detLog{i,1}.MeasurementNoise(2,2));
    fprintf(fileID,'P Matrix (2,3) %f\n', detLog{i,1}.MeasurementNoise(2,3));
    fprintf(fileID,'P Matrix (2,4) %f\n', detLog{i,1}.MeasurementNoise(2,4));
    fprintf(fileID,'P Matrix (2,5) %f\n', detLog{i,1}.MeasurementNoise(2,5));
    fprintf(fileID,'P Matrix (2,6) %f\n', detLog{i,1}.MeasurementNoise(2,6));
    fprintf(fileID,'P Matrix (3,3) %f\n', detLog{i,1}.MeasurementNoise(3,3));
    fprintf(fileID,'P Matrix (3,4) %f\n', detLog{i,1}.MeasurementNoise(3,4));
    fprintf(fileID,'P Matrix (3,5) %f\n', detLog{i,1}.MeasurementNoise(3,5));
    fprintf(fileID,'P Matrix (3,6) %f\n', detLog{i,1}.MeasurementNoise(3,6));
    fprintf(fileID,'P Matrix (4,4) %f\n', detLog{i,1}.MeasurementNoise(4,4));
    fprintf(fileID,'P Matrix (4,5) %f\n', detLog{i,1}.MeasurementNoise(4,5));
    fprintf(fileID,'P Matrix (4,6) %f\n', detLog{i,1}.MeasurementNoise(4,6));
    fprintf(fileID,'P Matrix (5,5) %f\n', detLog{i,1}.MeasurementNoise(5,5));
    fprintf(fileID,'P Matrix (5,6) %f\n', detLog{i,1}.MeasurementNoise(5,6));
    fprintf(fileID,'P Matrix (6,6) %f\n', detLog{i,1}.MeasurementNoise(6,6));
end
fprintf(fileID, '\n');
end

Answer 1

[Peter Khomchuk]

The radar model used in monostaticRadarSensor takes into account the effects of atmospheric refraction. This is the reason that monostaticRadarSensor does not return ‘ideal’ measurements even when the noise variances and the biases are zero. The atmospheric refraction influences the range and the angle measurements which in turn produces the observed discrepancy with the true values. Note that this effect is different from random measurement noise. It is caused by variation of the air density with height which causes bending of the electromagnetic waves. monostaticRadarSensor does not provide a way to turn the atmospheric refraction off, so it is always present in a simulation.