模拟和跟踪在地球中心的场景中飞行的飞机

打开实时脚本

这个例子向您展示了如何使用地球中心trackingScenario和一个geoTrajectory物体来模拟跨越数千公里的飞行轨迹。您使用两种不同的模型来生成飞机的综合探测:单站雷达和ADS-B报告。您使用多目标跟踪器来估计飞机轨迹，比较跟踪性能，并探索ADS-B对整体跟踪质量的影响。

在美国，联邦航空管理局(FAA)负责管理整个国家空域每天数以千计的航班。商业航班通常从起飞机场到抵达机场，全程都被跟踪。空中交通管制系统是一个复杂的多层次系统。机场控制塔负责监测和处理机场附近的情况，而空中交通管制中心(ARTCC)负责跨越构成国家空域的各个区域的远程航路监视。

在过去的几十年里，空中交通/监视雷达的能力和复杂性已经显著增加。飞机上增加的应答器增加了雷达设施和飞机之间的双向通信，允许非常准确的位置估计和有利于控制中心的决策。到2020年，所有飞行在10,000英尺以上的飞机都必须配备自动相关监视广播(ADS-B)应答器，以广播其机载估计位置。该信息由空中交通管制中心接收和处理。

创建一个空中交通场景

首先创建一个以地球为中心的场景。

%保存当前随机发生器状态S = rng;为可预测的结果设置随机种子rng (2020);%创建场景scene = trackingScenario()“IsEarthCentered”,真的,“UpdateRate”1);

定义飞机模型和轨迹

的matfileflightwaypoints本例所附包含从威奇托到芝加哥的飞行轨迹的合成坐标和时间信息。你用ageoTrajectory对象来创建飞行轨迹。

负载(“flightwaypoints.mat”) flightRoute = geoTrajectory(lla，到达时间);飞机=平台(场景，“轨迹”, flightRoute);

如今，商用飞机都配备了GPS接收器。作为ADS-B的骨干，星载GPS的精度可以设定为符合ADS-B的要求。ADS-B中用于定位和速度的导航精度类别分别称为NACp和NACv。根据FAA规定[1]，NACp必须小于0.05海里，NACv必须小于10米/秒。在本例中，您使用agpsSensor模型以50 m的位置精度和10 m/s的速度精度配置adsbTransponder模型。受波音737的启发，您还为飞机使用了更逼真的RCS签名。

posAccuracy = 50;%米velAccuracy = 10;% m / sgps = gpsSensor(“PositionInputFormat”，“大地”，“HorizontalPositionAccuracy”posAccuracy,…“VerticalPositionAccuracy”posAccuracy,“VelocityAccuracy”10);adsbTx = adsbtranssponder (“MW2020”，“UpdateRate”， 1,“全球定位系统”全球定位系统(gps),“类别”, adsbCategory.Large);负载(737 rcs.mat）;飞机。签名{1}=波音737rcs;

沿途增加监测站

美国联邦航空局使用的远程监视雷达有几种型号。空中航线监视雷达4 (ARSR-4)是在1990年代推出的一种雷达，可以在250海里(463公里)的远程范围内提供任何1平方米物体的3D返回。大多数ARSR-4雷达位于美国大陆边界，而射程稍短的雷达大多位于美国联邦航空局在美国大陆的雷达站。在本例中，按照ARSR-4的常用规格对单个雷达类型进行建模，如下所示:

更新速度:12秒
最大射程(1平方米目标):463公里
距离分辨率:323米
射程精度:116米
方位视野:360度
方位分辨率:1.4度
方位精度:0.176度
高度精度:900米

每个雷达站场景增加一个平台。这些平台的RCS签名设置为 $-$ 50分贝以避免产生不必要的雷达回波。

默认情况下，雷达检测报告在雷达安装平台主体框架中，在本例中为每个雷达站位置的局部东北下架框架。但是，在本例中，您将设置DetectionCoordinates财产场景以便在以地球为中心的地球固定帧(ECEF)中输出检测结果，这使得跟踪器可以在一个共同帧中处理来自不同雷达站点的所有检测结果。

模拟ARSR-4雷达update = 1/12;Fov = [360;30.001];elAccuracy = atan2d(0.9463);% 900m精度@最大范围elBiasFraction = 0.1;arsr4 = fusionRadarSensor(1，“UpdateRate”updaterate,…“FieldOfView”, 360;30.001),…“HasElevation”,真的,…“ScanMode”，“机械”，…“MechanicalAzimuthLimits”360年[0],…“MechanicalElevationLimits”(-30 0),…“之内”,真的,…“HasRangeRate”,真的,…“HasFalseAlarms”假的,…“ReferenceRange”, 463000,…“ReferenceRCS”0,…“RangeLimits”， [0463000]，…“AzimuthResolution”, 1.4,…“AzimuthBiasFraction”, 0.176/1.4,…“ElevationResolution”elAccuracy / elBiasFraction…“ElevationBiasFraction”elBiasFraction,…“RangeResolution”, 323,…“RangeBiasFraction”, 116/323,…精度/分辨率“RangeRateResolution”, 100,…“DetectionCoordinates”，“场景”）;%在每个ARTCC站点增加ARSR-4雷达radarsitesLLA = [41.4228 -88.0583] 0;…40.6989 -89.8258 0;…[39.2219 -95.2461];为I =1:3 radar = clone(arsr4);雷达。SensorIndex = i;平台(场景,“位置”radarsitesLLA(我,:)…“签名”rcsSignature (“模式”, -50),“传感器”,{雷达});结束

你使用adsbReceiver来模拟ADS-B信息的接收。ADS-B信息包含飞机自带GPS仪器测量的位置。该信息通常被编码并在1090mhz信道上广播，以供附近的ADS-B接收器接收。每个监测站的接收范围为200公里。在这个例子中，你假设监控站之间有完美的通信。因此，中央接收器在至少一个电台的范围内接收广播的ADS-B消息。

%定义中心adsbReceiveradsbRx = adsbReceiver(“ReceiverIndex”2);adsbRange = 2e5;

想象场景

你使用trackingGlobeViewer显示地球上的平台、轨迹、探测和轨迹。

查看器= trackingGlobeViewer(“技术”，“streets-dark”，…“TrackLabelScale”, 1.3,“TrackHistoryDepth”, 4000,…“CoverageMode”，“报道”）;%显示雷达站点plotPlatform(观众,scene.Platforms(2:结束);显示雷达覆盖covcon = coverageConfig(场景);covcon plotCoverage(查看器,“ECEF”，“α”, 0.1);%显示飞行路线plotTrajectory(观众、flightRoute);

监视雷达有一个天线盲锥，有时被称为“沉默锥”。它是雷达正上方的一大块空间，由于天线扫描的限制而无法被监视。雷达网络重叠覆盖是这一盲区的常见缓解策略。然而，使用重叠策略，仍然可能存在网络未完全覆盖的区域。在这个例子中，最南端的雷达的沉默锥(上图中橙色显示)仅部分被网络中相邻的雷达覆盖。这就造成了一个盲点，飞机将不会被任何雷达发现。

定义一个中央雷达跟踪器和一个跟踪引信。

通常，一个ARTCC维护其监视区域内所有物体的轨迹，并在物体进入新区域时将跟踪信息传递给下一个ARTCC。在本例中，您为所有雷达定义了一个集中式跟踪器。你用atrackerGNN目标融合多部雷达对飞机的雷达探测。

radarGNN = trackerGNN(“TrackerIndex”, 1…“MaxNumTracks”, 50岁,…“FilterInitializationFcn”@initfilter,…“ConfirmationThreshold”[3 - 5],…“DeletionThreshold”, 5 [5],…“AssignmentThreshold”2000年[250]);

还将雷达航迹与ADS-B接收机获得的ADS-B航迹融合。为此，您需要配置一个中心trackFuser对象。你设置ConfirmationThreshold和DeletionThreshold以解释ADS-B接收器速率为1 Hz和雷达跟踪器速率为1/12 Hz之间的更新速率差异。为了允许至少两个雷达航迹被分配到中央航迹，命中次数必须至少为 $2 \times 12$ ．

fuser = trackFuser“FuserIndex”3,“MaxNumSources”2,…“AssignmentThreshold”(250 500),…“StateFusion”，“十字路口”，…“StateFusionParameters”，“跟踪”，…“ConfirmationThreshold”(2 3 * 12),…“DeletionThreshold”4 [4] * 12);

用雷达和ADS-B跟踪飞机

在本节中，您将模拟该场景并对雷达跟踪器、跟踪引信器、ADS-B应答器和接收机进行步进。

%声明一些变量detBuffer = {};radarTrackLog = {};fusedTrackLog = {};adsbTrackLog = {};图像= {};snapTimes = [84,564,1128, 2083];跟踪标签和颜色adsblabel =“ads - b”;radarlabel =“雷达”;Fusedlabel = string(“% s \ n”，”“，“融合”));Adsbclr = [183 70 255]/255;Radarclr = [255 255 17]/255;Fusedclr = [255 105 41]/255;而advance(scene) time = scene. simulationtime;更新飞机在地球上的位置plotPlatform(观众、飞机、“TrajectoryMode”，“没有”）;生成和绘制合成雷达探测[dets, configs] =检测(场景);dets = postProcessDetection(dets);detBuffer = [detBuffer;依据);% #好< AGROW >detBuffer plotDetection(查看器,“ECEF”）;%的踪迹adsbTracks = objectTrack.empty;radarTracks = objectTrack.empty;fusedTracks = objectTrack.empty;%生成ADS-B消息truePose =姿态(飞机，“真正的”，“CoordinateSystem”，“大地”）;adsbMessage = adsbTx(truePose。位置,truePose.Velocity);在每个监测站的范围内可以接收到信息如果isWithinRange (radarsitesLLA truePose。位置，adsbRange) adsbTracks = adsbRx(adsbMessage，时间);结束在扫描结束时更新雷达跟踪器如果all([configs.IsValidTime]) && (~isempty(detBuffer) || isLocked(radarGNN)) radarTracks = radarGNN(detBuffer,time);detBuffer = {};结束%熔断器轨迹如果isLocked(fuser) || ~isempty([adsbTracks;radarTracks]) fusedTracks = fuser([adsbTracks;radarTracks]，time);结束%绘图只跟踪一次雷达扫描如果all([configs.IsValidTime]) labels = [repmat(adsblabel,1,numel(adsbTracks))，…repmat (radarlabel 1元素个数(radarTracks)),…repmat (fusedlabel 1元素个数(fusedTracks)));color = [repmat(adsbclr, numel(adsbTracks)， 1);…repmat(radarclr, numel(radarTracks)， 1);…repmat(fusedclr, numel(fusedTracks)， 1)];plotTrack(观众[adsbTracks;radarTracks;fusedTracks),“ECEF”，…“LabelStyle”，“自定义”，“CustomLabel”、标签“颜色”、颜色、…“线宽”3);结束记录估计的飞机数据fusedTrackLog = [fusedTrackLog， {fusedTracks}];% #好< AGROW >radarTrackLog = [radarTrackLog， {radarTracks}];% #好< AGROW >adsbTrackLog = [adsbTrackLog， {adsbTracks}];% #好< AGROW >移动相机，拍快照images = moveCamera(查看器，飞机，时间，snapTimes，图像);结束

图imshow(图片{1});

在场景开始时，飞机远离最南端的监测站，没有发送ADS-B信息。因此，飞机只能被雷达跟踪。请注意，ARSR在高度上的探测相对不准确，这通常是可以接受的，因为空中交通管制员将飞机水平分开，而不是垂直分开。最不精确的探测是由位于较远距离的雷达站点产生的。尽管如此，这些探测还是与轨道有关。图中，白线为真实轨迹，黄线为雷达跟踪器估计的轨迹。飞行的第一段距离远离任何雷达，相应的探测具有很高的测量噪声。此外，等速运动模型不能很好地模拟起飞后初始飞行转弯时的运动。雷达航迹被传递到引信器，引信器输出一个引信航迹，显示为橙色，紧跟雷达航迹。融合航迹不同于雷达航迹，因为融合器在航迹状态中加入了自身的过程噪声。

图像图imshow ({2});

在上面的快照中，飞机处于ADS-B通信范围内，并且建立了新的ADS-B航迹。熔断器对新轨迹进行了处理，提高了熔断轨迹的精度。

{3}图imshow(图像);

在上面的快照中，飞机进入了沉默的锥体。雷达跟踪器在几次更新后删除航迹，没有任何新的探测。此时，引信只依赖ADS-B轨迹来估计飞机的位置。

{4}图imshow(图像);

当飞机进入第二和第三监视雷达站所覆盖的区域时，建立了新的雷达航迹。来自两个雷达站的探测由雷达跟踪器融合，航迹引信器将新的雷达航迹与ADS-B航迹融合。

分析结果

你比较雷达和ADS-B记录的轨迹数据。真实的位置和速度可从geoTrajectory。您可以使用OSPA度量来比较仅从雷达、仅从ADS-B和从与ADS-B融合的雷达跟踪的质量。的默认设置trackOSPAMetric目标比较NEES(归一化估计误差平方)的轨道位置和检测分配。

ospa = trackOSPAMetric;radarospa = 0 (ceil(numel(radarTrackLog)/12)，1);%计算雷达跟踪器OSPA为i=1:12:numel(radarTrackLog) tracks = radarTrackLog{i};(真理。位置,~,真理。= lookupPose(flightRoute, i-1，“ECEF”）;Radarospa (i) = ospa(轨迹，真相);结束%计算ADS-B OSPAadsbospa = 0 (numel(adsbTrackLog)，1);重置(ospa);为i=1:numel(adsbTrackLog) tracks = adsbTrackLog{i};(真理。位置,~,真理。= lookupPose(flightRoute, i-1，“ECEF”）;Adsbospa (i) = ospa(轨迹，真相);结束%计算融合OSPAfusedTrackLog = 0 (numel(fusedTrackLog)，1);重置(ospa);为i=1: nummel (fusedTrackLog) tracks = fusedTrackLog{i};(真理。位置,~,真理。= lookupPose(flightRoute, i-1，“ECEF”）;Fusedospa (i) = ospa(轨迹，真相);结束

% Plot OSPA图保存在情节(0:12:(元素个数(radarTrackLog) 1)) / 60, radarospa(1:12:结束),“颜色”, radarclr);情节((0:(元素个数(adsbTrackLog) 1)) / 60, adsbospa,“颜色”adsbclr,“线宽”2);情节((0:(元素个数(fusedTrackLog) 1)) / 60, fusedospa,“颜色”, fusedclr);l =传奇(“雷达”，“ads - b”，“雷达+ ADS-B”）;l.Color = [0.1 0.1 0.1];l.TextColor = [1 1 1];包含(的时间(分钟)) ylabel (“OSPA”) ax = gca;网格在;盒子在;斧子。Color = [0.1 0.1 0.1];斧子。GridColor = [1 1 1];

OSPA度量显示了ADS-B航迹与雷达航迹融合所获得的改进。从模拟时间19分钟到25分钟，由于飞机在雷达网络的盲点上空飞行，雷达仅OSPA较高。ADS-B在该区域的可用性显著提高了跟踪性能，如融合OSPA所示。此外，该指标在开始时出现两个峰值，这可归因于在轨迹初始转弯期间恒速滤波器的性能不佳以及ADS-B的不可用性。大约在40分钟左右，由于区域内ADS-B可用性的丧失，仅ADS-B的OSPA被降级。在模拟的后面部分，雷达和ADS-B都可用。只有雷达的OSPA总体上比只有ADS-B的差。这是因为与GPS相比，雷达的垂直精度较差。

总结

在本例中，您学习了如何创建以地球为中心的场景，并使用大地坐标定义轨迹。您还学习了如何建模空中路线监视雷达和生成合成检测。您将这些检测结果提供给多目标跟踪器，并估计飞机的位置、速度和航向。通过加入和融合ADS-B信息，提高了跟踪性能。您对ADS-B报告进行建模，并将其集成到跟踪解决方案中。在本例中，只对一次飞行进行了建模。在空中交通管制员必须保持安全隔离距离的情况下，对多个航班进行建模时，ADS-B的好处可以得到进一步体现。

万博1manbetx支持功能

initfilter定义所使用的扩展卡尔曼滤波器trackerGNN．飞机运动可以用等速运动模型很好地近似。因此，使用相对较小的过程噪声来为动态分配更多的权重，而不是在长范围内预计会有相当大的噪声的测量。

函数Filter = initfilter(detection) Filter = initcvekf(detection);过滤器。StateCovariance = 100*filter.StateCovariance;% initcvekf使用测量噪声作为默认值过滤器。ProcessNoise = 0.1;结束

postProcessDetection通过两个操作对检测结果进行后处理:

它消除了位于地球表面以下的雷达探测，因为这些不能由真正的雷达产生。
它提高了垂直于径向分量的测速分量的检测噪声水平。

函数detsout = postProcessDetection(detsin) n = numel(detsin);Keep = 0 (1,n，“逻辑”）;为i=1:n meas = detsin{i}.Measurement(1:3)';Lla = fusion.internal.frames.ecef2lla(meas);如果Lla (1) = true;其他的Keep (i) = false;结束结束Detsout = detsin(keep);velCovTransversal = 100^2;为i=1:numel(detsout) velCov = detsout{i}. measuring noise (4:6,4:6);[P,D] = e (v);[m,ind] = min(diag(D));D = velCovTransversal * eye(3);D(ind,ind) = m;correctedCov = P*D*P';detsout{i}. measuring noise (4:6,4:6) = correctedCov;结束结束

isWithinRange如果飞机位置在任何监视站的ADS-B接收机范围内，则返回true。

函数flag = isWithinRange(stationsLLA, pos, range) cartDistance = fusion.internal.frames.lla2ecef(stationsLLA) - fusion.internal.frames.lla2ecef(pos);flag = any(vecnorm(cartDistance,2,2) < range);结束

moveCamera指定跟踪飞机并拍摄快照的新相机位置。

函数images = moveCamera(查看器，飞机，时间，snapTimes，图像)如果mod(time,12) == 0 pos = aircraft . position;如果Time == 0 campos(viewer, pos(1)，pos(2)， 2e4);elseif时间== 1100%缩小坎波斯(观众、pos (1), pos (2), 1 e5);elseif时间== 2000%放大坎波斯(观众、pos (1), pos (2), 2.6 e4);其他的保持原来的高度，但要跟飞机坎波斯(观众、pos (1), pos (2));结束结束%的快照如果任何(时间== snapTimes) drawnow img =快照(查看器);Images = [Images， {img}];结束结束