%创建跟踪场景。场景=跟踪Cenario（“StopTime”, 30);％定义单位转换。nmi2m = 1852;%海里到米的hr2 = 3600;小时到秒KTS2MPS = NMI2M / HR2S;%节到米每秒

%创建监视雷达。传感器= fusionRadarSensor (1,'没有扫描'，......“MountingLocation”(0 0 -20),......% 20米(ASL)“MountingAngles”(45 0 0),......%[偏航纵摇“FieldOfView”30 [10],......% [z el] deg“ReferenceRange”，5e3，......% m“AzimuthResolution”2,......%度“RangeResolution”5,......% m“之内”，真的，......％报告​​INS信息“TargetReportFormat”，“检测”，......%检测不聚类“DetectionCoordinates”，“球形传感器”）;

平台(场景中,“传感器”、传感器);塔=场景。平台{1}

定义两艘小船的尺寸。昏暗的=结构(......'长度', 80,......% m“宽度”15岁的......% m“高度”5,......% m'OriginOffset'，[0 0 5/2]）;(b) i ' m％型号为小船的雷达横截面（RCS）为30 dbsm。rcs = rcsSignature (“模式”, 30);创造一个转弯轨迹。速度= 20;%结Inityaw = 130;%度initPos = [1050 790 0];半径= 200;% minitOrient =四元数([initYaw 0 0]，“eulerd”，“ZYX股票”，“帧”）;engvel =速度* kts2mps * rotatepoint（机器，[1 0 0]）';AccBody = [0（速度* kts2mps）^ 2 /半径0];Angvelbody = [0 0速度* kts2mps /半径];traj = kinematictrajectory（“位置”，initpos，“速度”，entivel，'方向'initOrient,......'AccelerationSource'，“属性”，“加速”accBody,......“AngularVelocitySource”，“属性”，“AngularVelocity”, angVelBody);％加上第一艘小船，以20节向场行驶。这是离雷达塔最近的一艘船。平台(场景中,'方面'昏暗的,“签名”rcs,“轨迹”，traj）;创造另一艘以30节速度航行的小船。这就是那艘船距雷达塔最远的％。速度= 30;%结Inityaw = 120;%度initPos = [1410 1180 0];半径= 400;% minitOrient =四元数([initYaw 0 0]，“eulerd”，“ZYX股票”，“帧”）;engvel =速度* kts2mps * rotatepoint（机器，[1 0 0]）';AccBody = [0（速度* kts2mps）^ 2 /半径0];Angvelbody = [0 0速度* kts2mps /半径];traj = kinematictrajectory（“位置”，initpos，“速度”，entivel，'方向'initOrient,......'AccelerationSource'，“属性”，“加速”accBody,......“AngularVelocitySource”，“属性”，“AngularVelocity”, angVelBody);平台(场景中,'方面'昏暗的,“签名”rcs,“轨迹”，traj）;定义大型船的尺寸。昏暗的=结构(......'长度', 400,......% m“宽度”现年60岁的......% m“高度”15岁的......% m'OriginOffset'，[0 0 15/2]）;(b) i ' m%将大型船舶的雷达截面(RCS)建模为75dbsm。rcs = rcsSignature (“模式”, 75);％创建大型船舶的轨迹，以10节的常数行进。速度= 10;%结initYaw = -135;%度initPos = [1150 1100 0];initOrient =四元数([initYaw 0 0]，“eulerd”，“ZYX股票”，“帧”）;engvel =速度* kts2mps * rotatepoint（机器，[1 0 0]）';traj = kinematictrajectory（“位置”，initpos，“速度”，entivel，'方向'initOrient,......'AccelerationSource'，“属性”，“AngularVelocitySource”，“属性”）;％将大船添加到方案。平台(场景中,'方面'昏暗的,“签名”rcs,“轨迹”，traj）;%创建一个显示来显示船只的真实、测量和跟踪位置。theaterDisplay = helperMarineSurveillanceDisplay(场景中,......“IsSea”，真的，“DistanceUnits”，“米”，......“XLim”，450 * [ -  1 1] + 1E3，“YLim”，450 * [ -  1 1] + 1E3，“ZLim”-1000年[10],......“电影”，“MarineSurveillanceExample.gif”）;slctTrkPos = 0(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);

定义雷达的测量范围和分辨率。azLimits = sensor.FieldOfView(1)/2*[-1 1];%度= [0 15e3];% msensorLimits = [azLimits; rangeLimits];sensorResolution = [sensor.AzimuthResolution; sensor.RangeResolution];定义传感器在塔上的安装位置和方向。参数(1)=结构(“帧”，“球”，......“OriginPosition”sensor.MountingLocation (:)......'OriginVelocity'，[0; 0; 0]，......'方向', rotmat(四元数(传感器。MountingAngles,“eulerd”，“zyx股票”，“帧”）,“帧”）,......“IsParentToChild”,真的,......“HasRange”，真的，“HasElevation”、传感器。HasElevation,“HasVelocity”，错误的）;％定义塔的位置，速度和方案中的方向。params（2）= struct（“帧”，'矩形的'，......“OriginPosition”tower.Trajectory.Position (:)......'OriginVelocity'，tower.traptory.velocity（:)，......'方向'，rotmat（tower.trajectory.orientation，“帧”）,......“IsParentToChild”,真的,......“HasRange”，真的，“HasElevation”，错误的，“HasVelocity”，错误的）;％创建一个跟踪等级将模拟可检测性的跟踪传感器％曲目。检测概率定义了概率从扩展对象生成至少1个检测的百分比。SENSORCONFIG =跟踪转换仪控制功能（“SensorIndex”、传感器。SensorIndex,......'sensorlimits'sensorLimits,......“SensorResolution”sensorResolution,......“DetectionProbability”, 0.99,......“SensorTransformParameters”、参数);用于更新trackingSensorConfiguration的传感器配置的%字段。configFlds = {“SensorIndex”，“IsValidTime”}；对应于雷达的分辨率单元的噪声协方差。resolutionNoise =诊断接头((sensorResolution / 2) ^ 2);sensorConfig。FilterInitializationFcn = @(变长度输入宗量)filterInitFcn(变长度输入宗量{:},params);sensorConfig。SensorTransformFcn = @ctmeas;sensorConfig。ClutterDensity = sensor.FalseAlarmRate / (sensor.AzimuthResolution * sensor.RangeResolution);使用trackingSensorConfiguration创建一个PHD跟踪器。跟踪器= trackerphd（“SensorConfigurations”sensorConfig,......“HasSensorConfigurationsInput”，真的，......“PartitioningFcn”，@（x）零件detections（x，1.5,6），......“ExtractionThreshold”, 0.75,......“DeletionThreshold”，1e-6，......“出生率”1 e-5);

初始化场景和跟踪器。重启(场景);重置(跟踪);%设置模拟以雷达的更新速度前进。场景。UpdateRate = sensor.UpdateRate;%为可重复的结果设置随机种子。rng (2019“旋风”）;%的模拟运行。snapTimes = [2 7 scenario.StopTime];%秒而推进(场景)获得当前模拟时间。Time = Scenario.simulationtime;从塔的雷达上产生探测。[侦破,~,配置]=检测(塔、时间);更新探测的测量噪声以匹配雷达分辨率。DETS = UPDATEMEASURENCENOISE（DETS，DEMOTATENOISE）;%更新追踪。trackSensorConfig = computeTrackingSensorConfig(配置、configFlds);跟踪=追踪(精细、trackSensorConfig、时间);%用当前波束位置、探测和跟踪位置更新显示。theaterDisplay(引爆器,配置,跟踪);%的快照。snapFigure (theaterDisplay,任何(时间= = snapTimes));结束writeMovie (theaterDisplay);

showSnapshot (theaterDisplay, 1)

showSnapshot(theaterDisplay,2)轴([1250 1450 1150 1350]);视图(90 [-90]);

showSnapshot(theaterDisplay,3)轴([650 850 700 900]);视图(90 [-90]);

showSnapshot(theaterDisplay,3)轴([900 1100 1250 1450]);视图(90 [-90]);

numtrks = numel（曲目）;trackid = [tracks.trackid]';长度=零（numtrks，1）;宽度=零（numtrks，1）;height = zeros（numtrks，1）;为ITRK = 1：numtrks ext = track（iTrk）。扩展;[q，d] = eig（ex）;d = 2 * sqrt（DIAG（D））;idims = 1：3;UP = [0 0 -1];[〜，iup] = max（abs（up * q））;高度（ITRK）= D（IDIMS（IUP））;IDIMS（IUP）= [];长度（iTrk）= max（d（idims））;宽度（iTrk）= min（d（idims））;结束%显示船舶的估计尺寸。dim =表(TrackID、长度、宽度、高度)

函数config = computetrackingsensorconfig（configin，flds）config = struct（）;为iFld = 1:numel(flds) thisFld = flds{iFld};配置(thisFld) = configIn。(thisFld);结束结束

函数DETS = UPDATEMEASURENCENOISE（DETS，噪声）为dets = 1:numel(dets) dets{iDet}.MeasurementNoise(:) = noise(:);结束结束

函数phd = filterinitfcn（measparam，varargin）该函数仅使用预测出生密度来模拟出生%的场景。如果输入参数个数= = 1％在FOV内均匀设置预测诞生。phdDefault = initctggiwphd;% 1。使用方位角和距离创建均匀分布状态。阿兹= 0;方位角上的所有分量。范围= Linspace（1000,10000,5）;% 5组成范围[AZ，R] = MeshGrid（AZ，范围）;%创建一个PHD过滤器来分配内存。phd = ggiwphd(zeros(7,numel(Az))，repmat(eye(7)，[1 1 numel(Az)]))，......'scalematrices'repmat(眼(3)[1 1元素个数(Az)]),......'statetransitionfcn'@constturn,“StateTransitionJacobianFcn”@constturnjac,......“MeasurementFcn”，@ ctmeas，“MeasurementJacobianFcn”@ctmeasjac,......“PositionIndex”(1 3 6),“ExtentRotationFcn”, phdDefault。ExtentRotationFcn,......'hasadditiveProcessnoise'，错误的，“ProcessNoise”, 2 *眼(4),......'instalaldecay'1 e3,“GammaForgettingFactors”，1.1 * on（1，numel（az）），......'maxnumcomponents', 10000);为i = 1:numel(Az) [sensorX,sensorY,sensorZ] = sph2cart(deg2rad(Az(i))，0,R(i));globalPos = measParam(1).Orientation'*[sensorX;sensorY;sensorZ] + measParam(1).OriginPosition(:);博士学位。States([1 3 6]，i) = globalPos;博士学位。statecovariance ([1 3 6]，[1 3 6]，i) = diag([1e5 1e5 1000]);%使用位置协方差覆盖组件之间的间隙结束% 2。你已经描述了每艘船的“运动学”状态%在视野内。接下来，添加关于它们的大小和%预期检测数。％％预计海洋中有2种类型的船只，小而且%。您可以为每种大小创建组件。phdSmall =博士;%克隆PHD过滤器用于大型船只。phdLarge =克隆(博士);%设置初始组件数。numComps = phdSmall.NumComponents;对于小型船舶，预计长度和尺寸约为100米宽度% 20米。由于方向未知，我们将创建4每个尺寸的％取向。首先，您必须将组件添加到在相同状态下，%密度。这可以通过简单地附加它来完成%小船的设置值附加（phdsmall，phdsmall）;附加（phdsmall，phdsmall）;定义形状的％自由度。大量代表%更高的尺寸确定性。景深= 1000;vx vy和的协方差。smallStateCov = diag([300 300 50]);%比例矩阵的小船= (dof - 4)*diag([100/2 20/2 10].^2);% l, w, h％彼此以45度创建4个方向。为i = 1:4 thisIndex = (i-1)*numComps + (1:numComps); / /索引R = rotmat(四元数([45*(i-1) 0 0]，“eulerd”，“ZYX股票”，“帧”）,“帧”）;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 numcoms]); / /指定一个节点phdSmall。StateCovariances([6 7]，[6 7]，thisIndex) = repmat(diag([100 100])，[1 1 numComps]);结束％小船舶产生大约10-20个检测。expNumDets = 15;不确定性= 5 ^ 2;phdSmall.Rates (,) = expNumDets /不确定性;phdSmall.Shapes (,) = expNumDets ^ 2 /不确定性;phdSmall.DegreesOfFreedom(:) =景深;对于大型船舶遵循类似的流程。追加(phdLarge phdLarge);追加(phdLarge phdLarge);largeStateCov = diag([100 5 10]);= (dof - 4)*diag([500/2 100/2 10].^2);为i = 1:4 thisIndex = (i-1)*numComps + (1:numComps); / /索引R = rotmat(四元数([45*(i-1) 0 0]，“eulerd”，“ZYX股票”，“帧”）,“帧”）;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 numcoms]); / /将数据存储在指定的节点上phdLarge。StateCovariances([6 7]，[6 7]，thisIndex) = repmat(diag([100 100])，[1 1 numComps]);结束％生成大约100-200个检测。expNumDets = 150;不确定性= 50 ^ 2;phdLarge.Rates (,) = expNumDets /不确定性;phdLarge.Shapes (,) = expNumDets ^ 2 /不确定性;phdLarge.DegreesOfFreedom(:) =景深;将大的船加入到小的船中以创造总密度。这个密度％将％添加到每个步骤的总密度。博士= phdSmall;追加(博士,phdLarge);结束用检测输入调用时％，即适应性出生密度，不要%添加任何新组件。如果nargin> 1%这将在密度中创建0个组件。博士= initctggiwphd;结束结束