Введение

Статистическая радиолокационная модель

Статистические радиолокационные модели, такие как radarDataGenerator предложить ценные данные для ранней стадии разработки радиолокационных систем. Система моделируется с точки зрения её эксплуатационных характеристик, и минимальный набор параметров используется наряду с устоявшейся теорией, чтобы определить, как эта система будет взаимодействовать со своим окружением. Это включает спецификацию основных системных параметров, таких как рабочая частота, полоса пропускания и частота повторения импульсов (PRF), наряду с важными метриками, такими как угол и разрешение области значений, частота ложных предупреждений и вероятность обнаружения. Радиолокационная модель имеет специальные режимы для моностатической, бистатической и пассивной операции. Основным выходом этой модели являются необработанные данные обнаружения или последовательность обновлений трека. Несмотря на то, что сигналы временной области не генерируются, эмиттер может задать числовой ID для формы волны, которую он использует, так что можно учитывать пропорциональное усиление сигнала и отклонение. Обнаружения могут быть выведены в радиолокационной системе координат или в сценарии системы координат, если приемный радар знает местоположение и ориентацию передатчика. Дорожки генерируются с помощью одного из разнообразных алгоритмов отслеживания и управления дорожками. Скорость передачи данных обрабатывается спецификацией частоты обновления системы, так что обнаружения генерируются с соответствующей скоростью.

Управление сценариями

Модель работает в паре с radarScenario для генерации обнаружения или отслеживания данных с течением времени в динамическом окружении. Этот объект сценария не только управляет временем симуляции, но и может обрабатывать передачу структур данных между объектами в сцене по мере необходимости, когда запрашивается система координат обнаружения или обновления дорожки. Объект радара установлен на платформе, которая управляет динамикой, такой как положение, скорость и ориентация в системе координат сценария. Эти сценарные платформы также могут содержать информацию о сигнатуре, которая описывает, как платформа появляется на различных типах датчиков. Для наших целей, rcsSignature класс будет использоваться для придания целевым платформам радарного сечения.

Просмотр

Функциональность сканирования задается самим радаром, и может быть выполнена в механическом или электронном стиле. При механическом сканировании направление указания не может изменяться мгновенно, но ограничено максимальной скоростью вращения в азимуте и повышении. При электронном сканировании такого ограничения нет, и направление скана сбрасывается в начале каждого цикла. Радару также позволяют выполнять механическое и электронное сканирование одновременно. Рисунок ниже демонстрирует различие в шаблоне скана между этими двумя режимами. Заметьте, что в электронном режиме угол скана всегда увеличивается, а в механическом режиме положение скана всегда изменяется на смежное положение. Сканирование происходит сначала в азимуте, потому что это основное направление скана, которое является направлением с наибольшей степенью предела скана.

Сценарий: Воздушное наблюдение

Для этого сценария мы рассмотрим сканирующую радиолокационную систему с ее антенной, установленной на башне, как может быть использовано для слежения за самолетами. Мы просканируем относительно небольшую область, посмотрим на пару входящих платформ с различными профилями RCS, а также проверим обнаружение и отслеживаем данные. Мы увидим эффект ОСШ и область значений на нашу способность делать обнаружения.

Радиолокационная система

Используйте центральную частоту 1 ГГц и полосу 1,5 МГц, чтобы получить разрешение области значений около 100 метров. Чтобы смоделировать один импульс на каждое положение скана, установите частоту обновления равной желаемой PRF. Установите размер неоднозначности области значений, чтобы отразить наш выбор PRF, и установите наш верхний предел области значений в два раза больше, чем размер неоднозначности области значений, чтобы включить обнаружение объектов в неоднозначности первой области значений.

% Set the seed for repeatability
rng(2021);

% System parameters
c = physconst('lightspeed');
fc = 1e9; % Hz
bandwidth = 1.5e6; % Hz
prf = 4e3; % Hz
updateRate = prf;
rangeRes = c/(2*bandwidth); % m
rangeAmb = c/(2*prf); % m
rangeLims = [0 2*rangeAmb];

Помимо выбора между механическим и электронным режимами скана, поведение сканирования контролируется спецификацией пределов скана и поля зрения (FoV) по азимуту и повышению. FoV состоит из азимута и протяженности по повышению наблюдаемой области в каждом положении скана, подобной двухсторонней спецификации полуэнергетического луча. Точки скана выбираются так, чтобы общая угловая протяженность, заданная пределом скана, покрывалась неперекрывающимися сегментами с размером, равным FoV. Пределы скана заданы в системе координат радара, который может быть повернут от системы координат его платформы через MountingAngles свойство.

Мы просканируем +/-20 степени по азимуту и +/-10 степени по повышению. Установите значение FoV в 4 степени азимута и 8 степеней повышения. Наконец, укажите общее количество полных сканов, которые нужно моделировать, что будет диктовать общее время симуляции.

% Scanning parameters
azLimits = [-20 20];
elLimits = [-10 10];
fov = [4;8]; % deg
numFullScans = 20; % number of full scans to simulate

Общее количество точек скана можно найти следующим образом. Пределы скана обозначают минимальный и максимальный углы скана, привязанные к вектору boresight антенны, поэтому общая площадь, сканированная в одном направлении, больше, чем степень пределов скана, до половины FoV в этом направлении. Шаблон всегда является прямоугольной сеткой над az/el.

numScanPointsAz = floor(diff(azLimits)/fov(1)) + 1;
numScanPointsEl = floor(diff(elLimits)/fov(2)) + 1;
numScanPoints = numScanPointsAz*numScanPointsEl;

Будет использоваться спецификация отдельного от FoV углового разрешения, которая позволяет смоделировать нечто вроде оценки угла моноимпульса. Это угловое разрешение призвано определить нашу способность различать цели. и используется, чтобы вывести фактическую точность угла из нижней границы Cramer-Rao (CRLB) для оценки угла моноимпульса с этим разрешением. Точно так же разрешение области значений определяет минимальный интервал в области значений, требуемый для различения двух целей, но фактическая точность измерения области значений исходит из CRLB для оценки области значений. Точные измеряемые параметры могут использоваться путем выключения шума измерения с помощью HasNoise свойство. Как правило, ошибка измерения не переходит к 0, поскольку ОСШ увеличивается без привязки. Это может быть получено свойствами «смещения», которых по одному для каждого типа измеряемых. Для наших целей мы оставим эти по умолчанию.

Пусть наше угловое разрешение будет 1/4 FoV в каждом направлении, что позволило бы нам различать до 16 точек цели в FoV в той же области значений.

angRes = fov/4;

Ссылка по усилению цели и радиолокационного цикла

Вместо непосредственного определения таких вещей, как передающая степень, степени шума и специфика используемого алгоритма обнаружения, эта радарная модель использует концепцию радиолокации цикла усиления для преобразования RCS цели и области значений в ОСШ, который переводится непосредственно в вероятность обнаружения. Коэффициент усиления цикла радара является важным свойством любой радиолокационной системы и вычисляется ссылка целью, которая состоит из ссылки области значений и RCS и приемника рабочей характеристики (ROC) для заданных вероятностей обнаружения и ложного предупреждения.

Мы будем использовать вероятность обнаружения по умолчанию 90% для ссылки цели на 20 км и 0 дБсм. Мы будем использовать частоту ложных предупреждений по умолчанию 1e-6 ложных предупреждений на камеру разрешения за обновление.

refTgtRange = 20e3;
refTgtRCS = 0;
detProb = 0.9;
faRate = 1e-6;

Построение радара

Теперь давайте построим наш радар, используя параметры, указанные выше. Используйте нижний предел скана по повышению, чтобы задать угол установки тангажа радара с его платформы. Это указывает наше расположение так, что нижнее ребро области сканирования параллельно земле (то есть ни одна из наших точек скана не укажет радар в сторону земли).

pitch = elLimits(1); % mount the antenna rotated upwards

Включите измерения угла возвышения, генерацию ложных предупреждений и неоднозначности области значений с соответствующими Has... свойство. Чтобы смоделировать радар, который измеряет только область значений и угол, установите HasRangeRate на ложь. Чтобы получить обнаружения в координатах сценария, наш радар должен быть в курсе ориентации передатчика. Поскольку мы используем моностатический радар, этого можно достичь, включив функциональность инерционной навигационной системы (INS) с HasINS свойство.

% Create radar object
sensorIndex = 1; % a unique identifier is required
radar = radarDataGenerator(sensorIndex,'UpdateRate',updateRate,...
    'DetectionMode','Monostatic','ScanMode','Electronic','TargetReportFormat','Detections','DetectionCoordinates','scenario',...
    'HasElevation',true,'HasINS',true,'HasRangeRate',false,'HasRangeAmbiguities',true,'HasFalseAlarms',true,...
    'CenterFrequency',fc,'Bandwidth',bandwidth,...
    'RangeResolution',rangeRes,'AzimuthResolution',angRes(1),'ElevationResolution',angRes(2),...
    'ReferenceRange',refTgtRange,'ReferenceRCS',refTgtRCS,'DetectionProbability',detProb,'FalseAlarmRate',faRate,...
    'RangeLimits',rangeLims,'MaxUnambiguousRange',rangeAmb,'ElectronicAzimuthLimits',azLimits,'ElectronicElevationLimits',elLimits,...
    'FieldOfView',fov,'MountingAngles',[0 pitch 0]);

Предоставленный класс помощника может использоваться, чтобы визуализировать наш шаблон скана. Мы также будем использовать этот плоттер, чтобы анимировать шаблон скана в цикле симуляции.

scanplt = helperScanPatternDisplay(radar);
scanplt.makeOverviewPlot;

Обратите внимание, что наши точки скана в области азимута точно до заданных пределов скана азимута, в то время как точки скана повышения немного втягиваются, как необходимо, чтобы избежать перекрытия или избыточности точек скана.

Вычислите общее количество камер разрешения угла области значений для нашей системы и общее количество систем координат, которые будут использоваться для обнаружения. Затем мы можем вычислить ожидаемое количество ложных предупреждений путем умножения частоты ложных предупреждений на общее количество камер разрешения, запрашиваемых в течение симуляции.

numResolutionCells = diff(rangeLims)*prod(fov)/(rangeRes*prod(angRes));
numFrames = numFullScans*numScanPoints; % total number of frames for detection
expNumFAs = faRate*numFrames*numResolutionCells % expected number of FAs

expNumFAs = 7.9200

Сценарий и цель

Мы используем radarScenario для управления потоком симуляции верхнего уровня. Этот объект предоставляет функциональность, чтобы добавить новые платформы и быстро генерировать обнаружения от всех объектов датчика в сценарии. Он также управляет потоком времени в симуляции. Используйте общее количество систем координат, вычисленных выше, чтобы найти необходимую StopTimeи для адаптивного обновления времени симуляции на основе частоты обновления объектов в сцене установите частоту обновления сценария равной 0.

% Create scenario
stopTime = numFullScans*numScanPoints/prf;
scene = radarScenario('StopTime',stopTime,'UpdateRate',0);

Используйте platform метод для генерации новых объектов платформы, добавления их к сцене и возврата указателя для доступа к свойствам платформы. Каждая платформа имеет уникальный идентификатор, который присваивается сценарием при конструкции. Для статической платформы необходимо задать только свойство Position. Чтобы симулировать стандартную кинематическую траекторию постоянной скорости, используйте kinematicTrajectory и задайте положение и скорость.

Давайте создадим две целевые платформы в разных областях значений, как входящие с различными скоростями, так и углы.

% Create platforms
rdrPlat = platform(scene,'Position',[0 0 12]); % place 12 m above the origin
tgtPlat(1) = platform(scene,'Trajectory',kinematicTrajectory('Position',[38e3 6e3 10e3],'Velocity',[-380 -50 0]));
tgtPlat(2) = platform(scene,'Trajectory',kinematicTrajectory('Position',[25e3 -3e3 1e3],'Velocity',[-280 10 -10]));

Чтобы обеспечить удобство методов обнаружения и автоматического продвижения времени симуляции, сценарий должен быть осведомлен о радиолокационном объекте, построенном ранее. Для этого установите радар на платформу путем заполнения свойства Sensors платформы. Чтобы монтировать более одного датчика на платформе, свойство Sensors может быть массивом ячеек с объектами датчика.

% Mount radar to platform
rdrPlat.Sensors = radar;

The rcsSignature класс может использоваться, чтобы задать RCS платформы как функцию угла аспекта. Для простой сигнатуры constant-RCS просто установите Pattern свойство скалярному значению. Пусть одна платформа будет 20 дБсм, а другая 4 дБсм, чтобы продемонстрировать различие в обнаруживаемости.

% Set target platform RCS profiles
tgtRCS = [20, 4]; % dBsm
tgtPlat(1).Signatures = rcsSignature('Pattern',tgtRCS(1));
tgtPlat(2).Signatures = rcsSignature('Pattern',tgtRCS(2));

Обнаружительная способность

Давайте рассмотрим теоретический ОСШ, который можно ожидать для вне целей. Поскольку мы используем простую кинематическую траекторию, мы можем вычислить положение истинности и диапазон как таковой:

time = (0:numFrames-1).'/updateRate;
tgtPosTruth(:,:,1) = tgtPlat(1).Position + time*tgtPlat(1).Trajectory.Velocity;
tgtPosTruth(:,:,2) = tgtPlat(2).Position + time*tgtPlat(2).Trajectory.Velocity;
tgtRangeTruth(:,1) = sqrt(sum((tgtPosTruth(:,:,1)-rdrPlat.Position).^2,2));
tgtRangeTruth(:,2) = sqrt(sum((tgtPosTruth(:,:,2)-rdrPlat.Position).^2,2));

Затем мы можем использовать известные RCS цели и коэффициент усиления радиолокационного цикла, вычисленный радарной моделью, чтобы получить ОСШ истинности для каждой цели:

tgtSnr = tgtRCS - 40*log10(tgtRangeTruth) + radar.RadarLoopGain;

Минимальный ОСШ, требуемый для обнаружения, учитывая схему обнаружения только по величине, такую как CFAR, учитывая нашу вероятность обнаружения и частоту ложных предупреждений,

minSnr = 20*log10(erfcinv(2*radar.FalseAlarmRate) - erfcinv(2*radar.DetectionProbability)); % dB

Идя всего на бит дальше, мы можем вычислить нашу минимальную область значений для обнаружения:

Rd = 10.^((tgtRCS - minSnr + radar.RadarLoopGain)/40); % meters

Давайте сравним среднюю область значений наших целей и ОСШ с минимальной областью значений и порогами ОСШ. Общее расстояние перемещения наших целей достаточно мало, чтобы ОСШ на основе дальности не изменялся значительно.

snrMargin = mean(tgtSnr,1) - minSnr

snrMargin = 1×2

    8.0827    0.0019

rangeMargin = (Rd - mean(tgtRangeTruth,1))/1e3 % km

rangeMargin = 1×2

   23.5299    0.0028

Первая цель примерно на 8 дБ ярче, чем необходимо для обнаружения с заданными параметрами CFAR, в то время как вторая цель считается едва обнаруживаемой.

Моделирование обнаружений

Основной цикл симуляции начинается с вызова advance на объекте сценария. Этот метод переводит сценарий в следующий раз, когда объект в сцене определяется, требует обновления и возвращает false, когда достигается заданное время остановки, выходя из цикла. Обратите внимание, что первый вызов advance не шагает время симуляции, так что первая система координат набора данных может произойти на 0 с. В качестве альтернативы использованию advance, scene.SimulationStatus можно проверить, чтобы определить, выполнялся ли сеанс sim до своего завершения.

Внутри цикла методы обнаружения и генерации дорожек могут быть вызваны для любого объекта датчика индивидуально (см. Датчик step метод), но существуют удобные методы для генерации обнаружений и треков для всех датчиков и от всех целей в сцене с одним вызовом функции. Потому что у нас есть только один датчик, detect(scene) является хорошим решением, чтобы получить список обнаружений для всех целей в сцене. Инициируя генерацию обнаружений на верхнем уровне, подобном этому, сцена сама будет обрабатывать передачу требуемых временных параметров, INS и других строений данных датчику. Класс визуализации, используемый ранее, также будет использоваться здесь, чтобы анимировать шаблон скана.

allDets = []; % initialize list of all detections generated
lookAngs = zeros(2,numFrames);
simTime = zeros(1,numFrames);
frame = 0; % current frame
while advance(scene) % advance until StopTime is reached
    
    % Increment frame index and record simulation time
    frame = frame + 1;
    simTime(frame) = scene.SimulationTime;
    
    % Generate new detections
    dets = detect(scene);
    
    % Record look angle for inspection
    lookAngs(:,frame) = radar.LookAngle;
    
    % Update coverage plot
    scanplt.updatePlot(radar.LookAngle);
    
    % Compile any new detections
    allDets = [allDets;dets];
    
end

Поведение скана

Постройте график зарегистрированных углов обзора для первых двух полных сканов, чтобы просмотреть шаблон. Радар сначала сканирует по азимуту и поднимается по повышению в конце каждого азимутального скана, мгновенно сбрасывая угол обзора в начале каждого нового скана.

figure;
plot(simTime(1:2*numScanPoints)*1e3,lookAngs(:,1:2*numScanPoints));
xlabel('Time (ms)');
ylabel('Look Angle (deg)');
legend('Azimuth','Elevation');
title('Look Angles');

Обнаружения

Давайте посмотрим содержимое objectDetection как вывод нашего радара.

allDets{1} % look at the first detection

ans = 
  objectDetection with properties:

                     Time: 7.5000e-04
              Measurement: [3x1 double]
         MeasurementNoise: [3x3 double]
              SensorIndex: 1
            ObjectClassID: 0
    MeasurementParameters: [1x1 struct]
         ObjectAttributes: {[1x1 struct]}

Time - время симуляции, в которое было сгенерировано обнаружение. SensorIndex сообщает нам, каким датчиком сгенерировано это обнаружение (что важно, когда у нас есть несколько датчиков и мы используем метод обнаружения уровня сцены). Measurement - измеряемое значение, сопоставленное с обнаружением. Формат зависит от выбора нашего режима обнаружения и выходных координат. MeasurementNoise приводит отклонение или ковариацию измерения и используется в моделях трекера. Поскольку мы выводим обнаружения в координатах сценария, поле измерения является просто предполагаемым вектором положения цели, и шум измерения дает ковариацию этой оценки положения.

The ObjectAttributes поле содержит некоторые полезные метаданные:

allDets{1}.ObjectAttributes{1}

ans = struct with fields:
    TargetIndex: 3
            SNR: 12.5894

Это говорит нам, что обнаружение произошло с платформы с индексом 2, который является уникальным идентификатором, присвоенным сценарием нашей первой целевой платформе (второй платформе, которая будет добавлена к сцене). Мы также видим, что ОСШ этого обнаружения о том, что мы ожидали для нашей первой цели. Если обнаружение является ложным предупреждением, TargetIndex будет -1, что является недопустимым идентификатором платформы. Скомпилируйте целевой индекс, ОСШ и время каждого обнаружения.

detTgtIdx = cellfun(@(t) t.ObjectAttributes{1}.TargetIndex, allDets);
detTgtSnr = cellfun(@(t) t.ObjectAttributes{1}.SNR, allDets);
detTime   = cellfun(@(t) t.Time, allDets);
firstTgtDet  = find(detTgtIdx == tgtPlat(1).PlatformID);
secondTgtDet = find(detTgtIdx == tgtPlat(2).PlatformID);

Мы можем найти общее количество обнаружений FA с:

numFAs = numel(find(detTgtIdx==-1))

numFAs = 8

который близок к ожидаемому количеству FA, вычисленных ранее.

Соберите данные измерений между обнаружениями и постройте график вместе с положениями истинности для двух целевых платформ.

tgtPosObs = cell2mat(cellfun(@(t) t.Measurement,allDets,'UniformOutput',0).');
subplot(1,2,1);
helperPlotPositionError( tgtPosTruth(:,:,1),time,tgtPosObs(:,firstTgtDet),detTime(firstTgtDet),scene.SimulationTime );
title('Target 1 Position Error');
subplot(1,2,2);
helperPlotPositionError( tgtPosTruth(:,:,2),time,tgtPosObs(:,secondTgtDet),detTime(secondTgtDet),scene.SimulationTime );
title('Target 2 Position Error');
set(gcf,'Position',get(gcf,'Position')+[0 0 600 0]);

Большинство наших обнаружений происходит из первой цели, что неудивительно, учитывая анализ обнаруживаемости ранее. Однако он не был обнаружен в каждой системе координат, и вторая цель сгенерировала много обнаружений, несмотря на низкий ОСШ. Несмотря на то, что первая цель легко обнаруживается, ошибка измерения велика, потому что цель находится в нашей второй области значений неоднозначности, и никакие значения не были выполнены, в то время как наша вторая цель, когда обнаружена, показывает ошибку положения, которая согласуется с нашим разрешением областью значений 100 метров и относительно низким разрешением угла.

Давайте также рассмотрим общее отклонение обнаружений по ОСШ. Общая дисперсия является суммой предельных отклонений в каждом декартовом направлении (X, Y и Z). Это отклонение включает эффекты оценки в пространстве угол-диапазон наряду с преобразованием этой статистики в координаты сценария.

detTotVar = cellfun(@(t) trace(t.MeasurementNoise),allDets);
figure;
subplot(1,2,1);
helperPlotPositionTotVar( detTgtSnr(firstTgtDet),detTotVar(firstTgtDet) );
title('First Target');
subplot(1,2,2);
helperPlotPositionTotVar( detTgtSnr(secondTgtDet),detTotVar(secondTgtDet) );
title('Second Target');
set(gcf,'Position',get(gcf,'Position')+[0 0 600 0]);

Обратите внимание, что, хотя абсолютная ошибка положения для первой цели была больше из-за неоднозначности области значений, более низкое отклонение измерения отражает как больший ОСШ первой цели, так и ее, по-видимому, малый диапазон.

Симулируйте дорожки

The radarDataGenerator также способен выполнять отслеживание без разомкнутого контура. Вместо вывода объектов обнаружения, он может выводить обновления отслеживания.

Давайте еще раз запустим симуляцию, но на этот раз мы позволим нашей радарной модели выводить треки непосредственно. Во-первых, переустановите некоторые свойства модели радара, чтобы сгенерировать дорожки в системе координат сценария.

release(radar); % unlock for re-simulation
radar.TargetReportFormat = 'Tracks';
radar.TrackCoordinates = 'Scenario';

Отслеживание может быть выполнено несколькими различными алгоритмами. Существуют простые альфа-бета-фильтры и фильтры Калмана линейного, расширенного и неароматического разнообразия, наряду с моделями движения с постоянной скоростью, постоянным ускорением и постоянным поворотом. Алгоритм, который будет использоваться, задан как FilterInitializationFcn свойство. Это указатель на функцию или имя функции в вектор символов, т.е. @initcvekf или 'initcvekf', который принимает начальное обнаружение как вход и возвращает инициализированный объект трекера. Может использоваться любая определяемая пользователем функция с той же подписью. Мы будем использовать поставляемый фильтр Калмана с постоянной скоростью.

radar.FilterInitializationFcn = @initcvekf;

Общая структура цикла sim одинаковая, однако мы должны назвать step радара функцию вручную. Необходимыми входами являются struct целевых положений, полученная вызовом targetPoses на нашей радиолокационной платформе, struct INS, приобретенной вызовом pose на радиолокационной платформе и время симуляции тока.

restart(scene); % restart scenario

allTracks = []; % initialize list of all track data
while advance(scene)
    
    % Generate new track data
    tgtPoses = targetPoses(rdrPlat);
    ins = pose(rdrPlat);
    tracks = radar(tgtPoses,ins,scene.SimulationTime);
    
    % Compile any new track data
    allTracks = [allTracks;tracks];
    
end

Осмотр объекта дорожки:

allTracks(1) % look at first track update

ans = 
  objectTrack with properties:

             TrackID: 1
            BranchID: 0
         SourceIndex: 1
          UpdateTime: 0.0093
                 Age: 34
               State: [6x1 double]
     StateCovariance: [6x6 double]
     StateParameters: [1x1 struct]
       ObjectClassID: 0
          TrackLogic: 'History'
     TrackLogicState: [1 1 0 0 0]
         IsConfirmed: 1
           IsCoasted: 0
      IsSelfReported: 1
    ObjectAttributes: [1x1 struct]

TrackID является уникальным идентификатором файла дорожки и SourceIndex - идентификатор объекта трекера, сообщающий о дорожке, который полезен, когда в сцене есть несколько трекеров. Когда используется трекер с несколькими гипотезами, BranchID приводит индекс используемой гипотезы. The State поле содержит информацию о состоянии этого обновления. Потому что мы используем модель постоянной скорости и координаты дорожки кадра сценария, State состоит из векторов положения и скорости. UpdateTime задает время, в которое было сгенерировано это обновление дорожки. Другое важное свойство, IsCoasted, сообщает вам, использовалось ли это обновление дорожки новое обнаружение цели для обновления фильтра или просто распространялось вперед во времени от последнего обнаружения цели (дорожка «покрыта»).

Структура ObjectAttributes задает идентификатор платформы, сигнатура которой обнаруживается, и ОСШ этой подписи во время обновления.

allTracks(1).ObjectAttributes

ans = struct with fields:
    TargetIndex: 3
            SNR: 18.5910

Соберите целевой индекс и время обновления. Также проверьте, сколько обнаружений было использовано для обновления наших треков.

trackTgtIdx = arrayfun(@(t) t.TargetIndex,[allTracks.ObjectAttributes]);
updateTime = [allTracks.UpdateTime];
firstTgtTrack = find(trackTgtIdx == tgtPlat(1).PlatformID);
secondTgtTrack = find(trackTgtIdx == tgtPlat(2).PlatformID);

numUpdates = numel(find(~[allTracks.IsCoasted])) % total number of track updates with new detection data

numUpdates = 31

Это примерно количество не-FA обнаружений, сгенерированных в первой части.

Давайте рассмотрим отчетные позиции в плоскости XY и их ковариации, а также данные о положении истины. Поскольку первая цель находится в первой области значений неоднозначности и никакие значения не выполняются, давайте сначала найдем «истинное» неоднозначное положение для первой цели как таковой:

los = tgtPosTruth(:,:,1) - rdrPlat.Position;
R = sqrt(sum(los.^2,2));
los = los./R;
Ramb = mod(R,rangeAmb);
tgtPosTruthAmb = rdrPlat.Position + los.*Ramb;

Теперь используйте предоставленную функцию helper, чтобы построить визуализацию дорожки положения.

figure;
subplot(1,2,1);
helperPlotPositionTrack( tgtPosTruthAmb,allTracks(firstTgtTrack) );
title('First Target Track');
subplot(1,2,2);
helperPlotPositionTrack( tgtPosTruth(:,:,2),allTracks(secondTgtTrack) );
title('Second Target Track');
set(gcf,'Position',get(gcf,'Position')+[0 0 600 0]);

В обоих случаях оценки положения начинаются с сходимости по траектории истинности (или неоднозначной области значений). Для первой цели у нас есть больше обновлений треков, и окончательное состояние выглядит ближе к установившейся сходимости. Глядя на ковариационные контуры, мы видим ожидаемое начальное уменьшение размера, когда отслеживающий фильтр прогревается.

Чтобы присмотреться к нашей эффективности отслеживания, постройте график величины ошибки расчета положения для обеих целей и укажите, какие выборки были обновлены обнаружениями целей.

figure;
subplot(1,2,1);
helperPlotPositionTrackError( tgtPosTruthAmb,time,allTracks(firstTgtTrack),updateTime(firstTgtTrack) );
title('First Target Ambiguous Position Error');
subplot(1,2,2);
helperPlotPositionTrackError( tgtPosTruth(:,:,2),time,allTracks(secondTgtTrack),updateTime(secondTgtTrack) );
title('Second Target Position Error');
set(gcf,'Position',get(gcf,'Position')+[0 0 600 0]);

Для первой цели ошибка положения уменьшается на выборках, у которых было новое обнаружение для обновления файла дорожки (обозначено красным кругом). Поскольку вторая цель имела намного более низкий ОСШ, ошибка в некоторых отдельных обнаружениях была достаточно большой, чтобы увеличить ошибку положения дорожки. Несмотря на это, вторая ошибка положения целевой дорожки имеет небольшой тренд вниз, что указывает на то, что предполагаемое состояние, вероятно, сходится. Однако мы видим проблему с нашей способностью отслеживать первую цель. Поскольку он находится в первой области значений неоднозначности и имеет компонент движения, перпендикулярный линии зрения, кажущаяся скорость фактически меняется, и фильтр Калмана с постоянной скоростью не может сходиться.

Заключение

В этом примере мы сконфигурировали статистическую радиолокационную систему для использования в симуляции и анализе обнаруживаемости и отслеживаемости целей. Мы видели, как использовать пределы скана и свойства поля зрения для определения шаблона скана, как запустить радарную модель с помощью инструмента управления сценариями и как просмотреть сгенерированные обнаружения и дорожки.

Вспомогательные функции

function helperPlotPositionError( tgtPosTruth,time,tgtPosObs,detTime,T )

tgtPosTruthX = interp1(time,tgtPosTruth(:,1),detTime).';
tgtPosTruthY = interp1(time,tgtPosTruth(:,2),detTime).';
tgtPosTruthZ = interp1(time,tgtPosTruth(:,3),detTime).';

err = sqrt((tgtPosTruthX - tgtPosObs(1,:)).^2+(tgtPosTruthY - tgtPosObs(2,:)).^2+(tgtPosTruthZ - tgtPosObs(3,:)).^2);

plot(detTime*1e3,err/1e3,'--o');
grid on;
xlim([0 T*1e3]);
ylabel('RMS Error (km)');
xlabel('Time (ms)');

end

function helperPlotPositionTotVar( detTgtSnr,detTotVar )

plot(detTgtSnr,detTotVar,'--o');
grid on;
xlabel('SNR (dB)');
ylabel('Total Var (m^2)');

end

function helperPlotPositionTrack( tgtPosTruth,tracks )

plot(tgtPosTruth(:,1)/1e3,tgtPosTruth(:,2)/1e3);
hold on;

updateIdx = find(~[tracks.IsCoasted]);
theta = 0:0.01:2*pi;

for ind = 1:2:numel(updateIdx)
    
    % plot update
    T = tracks(updateIdx(ind));
    plot(T.State(1)/1e3,T.State(3)/1e3,'*black');
    sigma = T.StateCovariance([1 3],[1 3]);
    [evec,eval] = eigs(sigma);
    mag = sqrt(diag(eval)).';
    u = evec(:,1)*mag(1);
    v = evec(:,2)*mag(2);
    x = cos(theta)*u(1) + sin(theta)*v(1);
    y = cos(theta)*u(2) + sin(theta)*v(2);
    plot(x/1e3 + T.State(1)/1e3,y/1e3 + T.State(3)/1e3,'magenta');
    
end

hold off;
grid on;
xlabel('X (km)');
ylabel('Y (km)');
legend('Truth Position','Position Estimates','Covariance','Location','northwest');

end

function helperPlotPositionTrackError( tgtPosTruth,time,tracks,updateTime )

state = [tracks.State];

sx = state(1,:);
sy = state(3,:);
sz = state(5,:);

x = interp1(time,tgtPosTruth(:,1),updateTime);
y = interp1(time,tgtPosTruth(:,2),updateTime);
z = interp1(time,tgtPosTruth(:,3),updateTime);

err = sqrt((x-sx).^2+(y-sy).^2+(z-sz).^2);

isUpdate = ~[tracks.IsCoasted];

plot(updateTime*1e3,err);
hold on;
plot(updateTime(isUpdate)*1e3,err(isUpdate),'or');
hold off;
grid on;
xlabel('Update Time (ms)');
ylabel('Error (m)')
legend('Coasted','Update','Location','northwest');

end