Расширенное отслеживание объекта с радаром для морского наблюдения

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

Морской сценарий наблюдения

Симулируйте морской радар наблюдения, смонтированный наверху башни, пропускающей поставки в гавани. Симуляция местоположения башни и движение поставок в сценарии управляемы trackingScenario.

% Create tracking scenario.
scenario = trackingScenario('StopTime',30);

% Define unit conversions.
nmi2m = 1852;           % nautical miles to meters
hr2s = 3600;            % hours to seconds
kts2mps = nmi2m/hr2s;   % knots to meters per second

Морской радар наблюдения

Добавьте морской радар наблюдения в башню. Радар смонтирован 20 метров над уровнем моря (ASL). Радар смотрит в гавань, рассматривая 30 секторов азимута степени. Перечислены общие технические требования для морского радара наблюдения:

  • Чувствительность: 0 dBsm 5 км

  • Поле зрения: азимут на 30 градусов, вертикальное изменение на 10 градусов

  • Разрешение азимута: 2 градуса

  • Разрешение области значений: 5 м

Смоделируйте морской радар с указанными выше техническими условиями с помощью fusionRadarSensor.

% Create surveillance radar.
sensor = fusionRadarSensor(1,'No scanning', ...
    'MountingLocation',[0 0 -20], ...      % 20 meters (ASL)
    'MountingAngles',[45 0 0], ...         % [yaw pitch roll] deg
    'FieldOfView',[30 10], ...             % [az el] deg
    'ReferenceRange',5e3, ...              % m
    'AzimuthResolution',2, ...             % deg
    'RangeResolution',5, ...               % m
    'HasINS',true,...                      % Report INS information
    'TargetReportFormat','Detections',...  % Detections without clustering
    'DetectionCoordinates','Sensor spherical');

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

platform(scenario,'Sensors',sensor);
tower = scenario.Platforms{1}
tower = 

  Platform with properties:

       PlatformID: 1
          ClassID: 0
         Position: [0 0 0]
      Orientation: [0 0 0]
       Dimensions: [1x1 struct]
             Mesh: [1x1 extendedObjectMesh]
       Trajectory: [1x1 kinematicTrajectory]
    PoseEstimator: [1x1 insSensor]
         Emitters: {}
          Sensors: {[1x1 fusionRadarSensor]}
       Signatures: {[1x1 rcsSignature]  [1x1 irSignature]  [1x1 tsSignature]}

Добавьте три поставки в гавани в секторе наблюдения радара. Две меньших поставки поворачиваются на уровне 20 и 30 узлов, большая поставка перемещается в постоянном заголовке на уровне 10 узлов.

% Define the dimensions for the two small ships.
dim = struct( ...
    'Length',80, ... % m
    'Width',15, ...  % m
    'Height',5, ...  % m
    'OriginOffset', [0 0 5/2]); % [x y z] m

% Model the radar cross section (RCS) of the small ships as 30 dBsm.
rcs = rcsSignature('Pattern',30);

% Create a turning trajectory.
speed = 20;     % knots
initYaw = 130;  % deg
initPos = [1050 790 0];
radius = 200;   % m

initOrient = quaternion([initYaw 0 0], 'eulerd', 'ZYX', 'frame');
initVel = speed*kts2mps*rotatepoint(initOrient,[1 0 0])';
accBody = [0 (speed*kts2mps)^2/radius 0];
angVelBody = [0 0 speed*kts2mps/radius];
traj = kinematicTrajectory('Position',initPos,'Velocity',initVel,'Orientation',initOrient, ...
    'AccelerationSource','Property','Acceleration',accBody, ...
    'AngularVelocitySource','Property','AngularVelocity', angVelBody);

% Add the first small ship, traveling at 20 knots to the scenario. This is
% the closest ship to the radar tower.
platform(scenario,'Dimensions',dim,'Signatures',rcs,'Trajectory',traj);

% Create the other small ship, traveling at 30 knots. This is the ship
% which is farthest from the radar tower.
speed = 30;     % knots
initYaw = 120;  % deg
initPos = [1410 1180 0];
radius = 400;   % m

initOrient = quaternion([initYaw 0 0],'eulerd','ZYX','frame');
initVel = speed*kts2mps*rotatepoint(initOrient,[1 0 0])';
accBody = [0 (speed*kts2mps)^2/radius 0];
angVelBody = [0 0 speed*kts2mps/radius];
traj = kinematicTrajectory('Position',initPos,'Velocity',initVel,'Orientation',initOrient, ...
    'AccelerationSource','Property','Acceleration',accBody, ...
    'AngularVelocitySource','Property','AngularVelocity', angVelBody);

platform(scenario,'Dimensions',dim,'Signatures',rcs,'Trajectory',traj);

% Define the dimensions for the large ship.
dim = struct( ...
    'Length',400, ... % m
    'Width',60, ...   % m
    'Height',15, ...  % m
    'OriginOffset', [0 0 15/2]); % [x y z] m

% Model the radar cross section (RCS) of the large ship as 75 dBsm.
rcs = rcsSignature('Pattern',75);

% Create the large ship's trajectory, traveling at a constant heading at 10 knots.
speed = 10;     % knots
initYaw = -135; % deg
initPos = [1150 1100 0];

initOrient = quaternion([initYaw 0 0],'eulerd','ZYX','frame');
initVel = speed*kts2mps*rotatepoint(initOrient,[1 0 0])';
traj = kinematicTrajectory('Position',initPos,'Velocity',initVel,'Orientation',initOrient, ...
    'AccelerationSource','Property','AngularVelocitySource','Property');

% Add the large ship to the scenario.
platform(scenario,'Dimensions',dim,'Signatures',rcs,'Trajectory',traj);

% Create a display to show the true, measured, and tracked positions of the ships.
theaterDisplay = helperMarineSurveillanceDisplay(scenario, ...
    'IsSea',true,'DistanceUnits','m', ...
    'XLim',450*[-1 1]+1e3,'YLim',450*[-1 1]+1e3,'ZLim',[-1000 10], ...
    'Movie','MarineSurveillanceExample.gif');
slctTrkPos = zeros(3,7); slctTrkPos(1,1) = 1; slctTrkPos(2,3) = 1; slctTrkPos(3,6) = 1;
slctTrkVel = circshift(slctTrkPos,[0 1]);
theaterDisplay.TrackPositionSelector = slctTrkPos;
theaterDisplay.TrackVelocitySelector = slctTrkVel;
theaterDisplay();
snapnow(theaterDisplay);

Мультипредназначайтесь для средства отслеживания GGIW-PHD

Создайте trackerPHD сформировать дорожки из радарных обнаружений, сгенерированных от трех поставок в гавани. Средство отслеживания PHD включает оценку размера поставок, позволяя нескольким обнаружениям быть сопоставленным к отдельному объекту. Это важно в ситуациях, таких как морское наблюдение, где размер объектов, обнаруженных датчиком, больше разрешения датчика, приводящего к нескольким обнаружениям, сгенерированным вдоль поверхностей поставок.

Средство отслеживания использует filterInitFcn поддерживание функции, чтобы инициализировать постоянную угловую скорость вращения фильтр Gamma Gaussian Inverse Wishart (GGIW) PHD. filterInitFcn добавляют компоненты рождения к ИНТЕНСИВНОСТИ ДОКТОРА ФИЛОСОФИИ на каждом временном шаге. Эти компоненты рождения добавляются однородно в поле зрения датчика. Их размеры и ожидаемое количество обнаружений указаны с помощью предшествующей информации о типах поставок, ожидаемых в гавани.

Средство отслеживания использует гамма распределение компонентов GGIW-PHD, чтобы оценить, сколько обнаружений должно быть сгенерировано от объекта. Средство отслеживания также вычисляет обнаружительную способность каждого компонента в плотности с помощью пределов датчика. Используйте trackingSensorConfiguration смоделировать настройку датчика для trackerPHD.

% Define radar's measurement limits and resolution.
azLimits = sensor.FieldOfView(1)/2*[-1 1]; % deg
rangeLimits = [0 15e3];                    % m
sensorLimits = [azLimits;rangeLimits];
sensorResolution = [sensor.AzimuthResolution;sensor.RangeResolution];

% Define the sensor's mounting location and orientation on the tower.
params(1) = struct('Frame','Spherical', ...
    'OriginPosition',sensor.MountingLocation(:), ...
    'OriginVelocity',[0;0;0], ...
    'Orientation',rotmat(quaternion(sensor.MountingAngles,'eulerd','zyx','frame'),'frame'),...
    'IsParentToChild', true, ...
    'HasRange',true,'HasElevation',sensor.HasElevation,'HasVelocity',false);

% Define the tower's location, velocity, and orientation in the scenario.
params(2) = struct('Frame','Rectangular', ...
    'OriginPosition',tower.Trajectory.Position(:), ...
    'OriginVelocity',tower.Trajectory.Velocity(:), ...
    'Orientation',rotmat(tower.Trajectory.Orientation,'frame'), ...
    'IsParentToChild', true, ...
    'HasRange',true,'HasElevation',false,'HasVelocity',false);

% Create a trackingSensorConfiguration to model the detectability of the
% tracks by the sensor. The detection probability defines the probability
% of generating at least 1 detection from the extended object.
sensorConfig = trackingSensorConfiguration('SensorIndex',sensor.SensorIndex, ...
    'SensorLimits',sensorLimits,...
    'SensorResolution',sensorResolution,...
    'DetectionProbability',0.99,...
    'SensorTransformParameters',params);

% Fields on the sensor configuration used to update trackingSensorConfiguration.
configFlds = {'SensorIndex','IsValidTime'};

% Noise covariance corresponding to a resolution cell of the radar.
resolutionNoise = diag((sensorResolution/2).^2);

sensorConfig.FilterInitializationFcn = @(varargin)filterInitFcn(varargin{:},params);
sensorConfig.SensorTransformFcn = @ctmeas;
sensorConfig.ClutterDensity = sensor.FalseAlarmRate/(sensor.AzimuthResolution*sensor.RangeResolution);

% Create a PHD tracker using the trackingSensorConfiguration.
tracker = trackerPHD('SensorConfigurations',sensorConfig, ...
    'HasSensorConfigurationsInput',true, ...
    'PartitioningFcn',@(x)partitionDetections(x,1.5,6), ...
    'DeletionThreshold',1e-6,...
    'BirthRate',1e-5);

Симулируйте и отследите корабли

Следующий цикл совершенствует положения поставок до конца сценария. Для каждого шага вперед в сценарии средство отслеживания обновляется с обнаружениями от поставок в поле зрения радара.

% Initialize scenario and tracker.
restart(scenario);
reset(tracker);

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

% Set random seed for repeatable results.
rng(2019);

% Run simulation.
snapTimes = [2 7 scenario.StopTime]; % seconds
while advance(scenario)
    % Get current simulation time.
    time = scenario.SimulationTime;

    % Generate detections from the tower's radar.
    [dets,~,config] = detect(tower,time);
    config.SensorIndex = sensor.SensorIndex;

    % Update measurement noise of detections to match radar's resolution.
    dets = updateMeasurementNoise(dets,resolutionNoise);

    % Update tracker.
    trackSensorConfig = computeTrackingSensorConfig(config,configFlds);
    tracks = tracker(dets,trackSensorConfig,time);

    % Update display with current beam position, detections, and track positions.
    theaterDisplay(dets,config,tracks);

    % Take snapshot.
    snapFigure(theaterDisplay,any(time==snapTimes));
end
writeMovie(theaterDisplay);

Следующий рисунок показывает радарные обнаружения, показавшие красными точками и предполагаемыми местоположениями дорожки, показавшими желтыми квадратами, аннотируемыми ID дорожки и степенью предполагаемого отслеживаемого объекта, показавшей желтым эллипсом. Радарная башня расположена в начале координат, (0,0), который не показан на рисунке. Поле зрения радара обозначается этими двумя красными линиями, пересекающими верх и низ фигуры. Все поставки лежат в поле зрения радара, и потому что размер поставок намного больше, чем область значений радара и разрешение азимута, несколько обнаружений сделаны вдоль поверхностей поставок, видимых к радару.

showSnapshot(theaterDisplay,1)

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

showSnapshot(theaterDisplay,2)
axis([1250 1450 1150 1350]); view([-90 90]);

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

showSnapshot(theaterDisplay,3)
axis([650 850 700 900]); view([-90 90]);

Следующий рисунок также показывает, что средство отслеживания PHD оценило и местоположение, размер и заголовок другой маленькой поставки в сценарии. Это - поставка, которая была ранее закрыта другими двумя поставками. Несмотря на поглощение газов, предполагаемое местоположение, размер и ориентацию тесно совпадает с поставкой.

showSnapshot(theaterDisplay,3)
axis([900 1100 1250 1450]); view([-90 90]);

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

numTrks = numel(tracks);
TrackID = [tracks.TrackID]';
Length = zeros(numTrks,1);
Width = zeros(numTrks,1);
Height = zeros(numTrks,1);
for iTrk = 1:numTrks
    ext = tracks(iTrk).Extent;
    [Q,D] = eig(ext);
    d = 2*sqrt(diag(D));
    iDims = 1:3;

    up = [0 0 -1];
    [~,iUp] = max(abs(up*Q));
    Height(iTrk) = d(iDims(iUp));
    iDims(iUp) = [];

    Length(iTrk) = max(d(iDims));
    Width(iTrk) = min(d(iDims));
end

% Display a table of the estimated dimensions of the ships.
dims = table(TrackID,Length,Width,Height)
dims =

  3x4 table

    TrackID    Length    Width     Height
    _______    ______    ______    ______

       1        102.7    18.213    16.543
       2       472.05     55.85     8.953
       3       97.196     18.57    17.626

Вспомните, что истинными размерностями поставок дают:

Большая поставка

  • Длина: 400 м

  • Width: 60 м

  • Высота: 15 м

Маленькая поставка

  • Длина: 80 м

  • Width: 15 м

  • Высота: 5 м

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

Сводные данные

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

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

computeTrackingSensorConfig

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

function config = computeTrackingSensorConfig(configIn,flds)
config = struct();
for iFld = 1:numel(flds)
    thisFld = flds{iFld};
    config.(thisFld) = configIn.(thisFld);
end
end

updateMeasurementNoise

Устанавливает шум измерения обнаружений на основе заданной шумовой ковариации.

function dets = updateMeasurementNoise(dets,noise)
    for iDet = 1:numel(dets)
        dets{iDet}.MeasurementNoise(:) = noise(:);
    end
end

filterInitFcn

Изменяет фильтр, возвращенный initctggiwphd совпадать со скоростями и ожидаемым количеством обнаружений для отслеженных кораблей.

function phd = filterInitFcn(measParam,varargin)
% This function uses the predictive birth density only to simulate birth in
% the scenario.

if nargin == 1 % Set predictive birth uniformly inside the FOV.
    phdDefault = initctggiwphd;
    % 1. Create uniformly distributed states using Azimuth and Range.
    az = 0; % All components at azimuth.
    ranges = linspace(1000,10000,5); % 5 components for range
    [Az,R] = meshgrid(az,ranges);

    % create a PHD filter to allocate memory.
    phd = ggiwphd(zeros(7,numel(Az)),repmat(eye(7),[1 1 numel(Az)]),...
        'ScaleMatrices',repmat(eye(3),[1 1 numel(Az)]),...
        'StateTransitionFcn',@constturn,'StateTransitionJacobianFcn',@constturnjac,...
        'MeasurementFcn',@ctmeas,'MeasurementJacobianFcn',@ctmeasjac,...
        'PositionIndex',[1 3 6],'ExtentRotationFcn',phdDefault.ExtentRotationFcn,...
        'HasAdditiveProcessNoise',false,'ProcessNoise',2*eye(4),...
        'TemporalDecay',1e3,'GammaForgettingFactors',1.1*ones(1,numel(Az)),...
        'MaxNumComponents',10000);

    for i = 1:numel(Az)
        [sensorX,sensorY,sensorZ] = sph2cart(deg2rad(Az(i)),0,R(i));
        globalPos = measParam(1).Orientation'*[sensorX;sensorY;sensorZ] + measParam(1).OriginPosition(:);
        phd.States([1 3 6],i) = globalPos;
        phd.StateCovariances([1 3 6],[1 3 6],i) = diag([1e5 1e5 1000]); % Cover gaps between components using position covariance
    end

    % 2. You have described the "kinematic" states of each of the ships
    % inside the field of view. Next, add information about their sizes and
    % expected number of detections.
    %
    % It is expected that there are 2 types of ships in the sea, small and
    % large. You create components for each size.
    phdSmall = phd;

    % Clone the PHD filter for large ships.
    phdLarge = clone(phd);

    % Set initial number of components.
    numComps = phdSmall.NumComponents;

    % For small ships, the expected size is about 100 meters in length and
    % 20 meters in width. As the orientation is unknown, we will create 4
    % orientations for each size. First, you must add components to the
    % density at same states. This can be done by simply appending it
    % Setup values for small boats
    append(phdSmall,phdSmall);
    append(phdSmall,phdSmall);

    % Degrees of freedom for defining shape. A large number represents
    % higher certainty in dimensions.
    dof = 1000;

    % Covariance in vx, vy and omega.
    smallStateCov = diag([300 300 50]);

    % Scale matrix for small boats
    smallShape = (dof - 4)*diag([100/2 20/2 10].^2); % l, w and h

    % Create 4 orientations at 45 degrees from each other.
    for i = 1:4
        thisIndex = (i-1)*numComps + (1:numComps);
        R = rotmat(quaternion([45*(i-1) 0 0],'eulerd','ZYX','frame'),'frame');
        phdSmall.ScaleMatrices(:,:,thisIndex) = repmat(R*smallShape*R',[1 1 numComps]);
        phdSmall.StateCovariances([2 4 5],[2 4 5],thisIndex) = repmat(R*smallStateCov*R',[1 1 numComps]);
        phdSmall.StateCovariances([6 7],[6 7],thisIndex) = repmat(diag([100 100]),[1 1 numComps]);
    end

    % Small ships generate approximately 10-20 detections.
    expNumDets = 15;
    uncertainty = 5^2;
    phdSmall.Rates(:) = expNumDets/uncertainty;
    phdSmall.Shapes(:) = expNumDets^2/uncertainty;
    phdSmall.DegreesOfFreedom(:) = dof;

    % Follow similar process for large ships.
    append(phdLarge,phdLarge);
    append(phdLarge,phdLarge);
    largeStateCov = diag([100 5 10]);
    largeShape = (dof - 4)*diag([500/2 100/2 10].^2);
    for i = 1:4
        thisIndex = (i-1)*numComps + (1:numComps);
        R = rotmat(quaternion([45*(i-1) 0 0],'eulerd','ZYX','frame'),'frame');
        phdLarge.ScaleMatrices(:,:,thisIndex) = repmat(R*largeShape*R',[1 1 numComps]);
        phdLarge.StateCovariances([2 4 5],[2 4 5],thisIndex) = repmat(R*largeStateCov*R',[1 1 numComps]);
        phdLarge.StateCovariances([6 7],[6 7],thisIndex) = repmat(diag([100 100]),[1 1 numComps]);
    end

    % Generate approximately 100-200 detections.
    expNumDets = 150;
    uncertainty = 50^2;
    phdLarge.Rates(:) = expNumDets/uncertainty;
    phdLarge.Shapes(:) = expNumDets^2/uncertainty;
    phdLarge.DegreesOfFreedom(:) = dof;

    % Add large ships to small ships to create total density. This density
    % is added to the total density every step.
    phd = phdSmall;
    append(phd,phdLarge);
end

% When called with detection input i.e. the adaptive birth density, do not
% add any new components.
if nargin > 1
    % This creates 0 components in the density.
    phd = initctggiwphd;
end

end
Для просмотра документации необходимо авторизоваться на сайте