В этом примере показано, как сгенерировать морской сценарий, симулируйте радарные обнаружения от морского радара наблюдения и сконфигурируйте мультицелевое средство отслеживания Плотности гипотезы вероятности (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 м
Смоделируйте морской радар с указанными выше техническими условиями с помощью monostaticRadarSensor.
% Create surveillance radar. sensor = monostaticRadarSensor(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,'DetectionCoordinates','Sensor spherical');
Добавьте башню в сценарий как стационарная платформа с радаром, смонтированным сверху его.
platform(scenario,'Sensors',sensor);
tower = scenario.Platforms{1}
tower =
Platform with properties:
PlatformID: 1
ClassID: 0
Dimensions: [1x1 struct]
Trajectory: [1x1 kinematicTrajectory]
PoseEstimator: [1x1 insSensor]
Emitters: {}
Sensors: {[1x1 monostaticRadarSensor]}
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);
Создайте 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. sensorConfig = trackingSensorConfiguration('SensorIndex',sensor.SensorIndex, ... 'SensorLimits',sensorLimits,... 'SensorResolution',sensorResolution,... '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); % 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 99.38 17.597 16.464
2 474.46 52.653 8.6582
3 102.33 20.182 17.635
Вспомните, что истинными размерностями поставок дают:
Большая поставка
Длина: 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 100]); % 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