利用地面雷达探测和跟踪LEO卫星星座
本示例演示如何导入卫星星座的TLE (Two-Line Element)文件，模拟星座的雷达探测，并跟踪星座。
填充和维护围绕地球运行的空间物体目录的任务在空间监视中至关重要。这项任务包括几个过程:探测和识别新天体并将其添加到目录中，更新目录中已知天体的轨道，跟踪它们整个生命周期的轨道变化，以及预测重返大气层。在本例中，我们将研究如何探测和跟踪新卫星，并将其添加到目录中。
为保证太空安全运行，防止与其他卫星或已知碎片碰撞，正确探测和跟踪新发射的卫星非常重要。太空机构通常分享发射前的信息，这些信息可用于选择搜索策略。低地球轨道(LEO)卫星搜索策略由围栏式雷达系统组成。围栏式雷达系统在空间中搜索有限体积，并在卫星通过其视野时进行探测。该策略可以快速探测和跟踪新发射的星座。[1]
从TLE文件中导入卫星星座
双线元集是保存卫星轨道信息的常用数据格式。您可以使用satelliteScenario对象导入TLE文件中定义的卫星轨道。缺省情况下，导入的卫星轨道采用SGP4轨道传播算法进行传播，该算法为LEO目标提供了较好的精度。在这个例子中，这些轨道提供了地面真相，以测试雷达跟踪系统探测新发射卫星的能力。
创建一个卫星场景卫星场景=卫星场景;从TLE文件中添加卫星。tleFile =“leoSatelliteConstellation.tle”；星座=卫星(satscene, tleFile);
使用卫星场景查看器来可视化星座。
玩(satscene);
模拟合成探测和轨道星座
空间监视雷达建模
定义两个用扇形雷达波束观测太空的站。风扇穿过卫星轨道，以最大限度地增加探测到的数量。位于北美的雷达站形成了一个东西篱笆。
% LLA第一站坐标Station1 = [48 -80 0];% LLA中第二站坐标Station2 = [50 -117 0];
每个站都配备了一个雷达，该雷达是用fusionRadarSensor对象。为了探测LEO范围内的卫星，雷达有以下要求:
探测2000公里外的10 dBsm物体
在2000公里范围内水平和垂直分辨目标，精度为100米
有120度方位角和30度仰角的视野
仰望太空
%创建扇形单站雷达Fov = [120;40];radar1 = fusionRadarSensor(1，.．.“UpdateRate”, 0.1,.．.10秒“ScanMode”，“没有扫描”，.．.“MountingAngles”，[0 90 0]，.．.查找“FieldOfView”fov,.．.度“ReferenceRange”2000年e3,.．.米“RangeLimits”， [0 2000e3]，.．.米“ReferenceRCS”10.．.dBsm“HasFalseAlarms”假的,.．.“HasNoise”,真的,.．.“HasElevation”,真的,.．.“AzimuthResolution”, 0.03,.．.度“ElevationResolution”, 0.03,.．.度“RangeResolution”, 2000,.．.M %精度~= 2000 * 0.05 (M)“DetectionCoordinates”，“球形传感器”，.．.“TargetReportFormat”，“检测”）;Radar2 =克隆(radar1);radar2。SensorIndex = 2;
雷达处理链
在本例中，执行了几个坐标转换和轴转换以正确运行雷达跟踪链。下图说明了上述定义的输入是如何转换并传递给雷达的。
在第一步中，计算每个卫星在本地雷达站NED轴上的姿态。首先要获得地面站ECEF姿态，并将卫星位置和速度转换为ECEF坐标。雷达输入是将卫星位姿与地面站位姿的差值旋转成地面站局部NED轴。看到assembleRadarInputs实现细节的支持功能。
在第二步中，向检测对象添加所需的信息，以便跟踪器可以使用ECEF状态进行操作。你可以使用MeasurementParameters属性在每个对象检测中实现该目的，如图所示addMeasurementParams支持功能。
定义跟踪器
您在上面定义的雷达模型输出探测。为了估计卫星轨道，你要使用跟踪器。传感器融合和跟踪工具箱™提供了多种多目标跟踪器。在本例中，您选择联合概率数据关联(JPDA)跟踪器，因为它能很好地平衡跟踪性能和计算成本。
您需要为跟踪器定义一个跟踪过滤器。您可以使用比SGP4保真度低的模型，例如运动方程的开普勒积分来跟踪卫星。通常，目标运动模型中保真度的缺乏通过测量更新和在滤波器中加入过程噪声来补偿。支持函数initKeplerUKF定义跟踪筛选器。
定义跟踪器追踪器=追踪器“FilterInitializationFcn”@initKeplerUKF,.．.“HasDetectableTrackIDsInput”,真的,.．.“ClutterDensity”1 e-40.．.“AssignmentThreshold”1 e4,.．.“DeletionThreshold”, 8 [5],.．.“ConfirmationThreshold”8 [5]);
运行模拟
在本例的其余部分中，您将逐步完成该场景来模拟雷达探测和跟踪卫星。本节使用trackingGlobeViewer可视化。您可以使用这个类显示带有不确定性椭圆的传感器和跟踪数据，并显示每个卫星的真实位置。
查看器= trackingGlobeViewer(“ShowDroppedTracks”假的,“PlatformHistoryDepth”, 700);定义每个雷达的覆盖配置，并在地球上可视化。Ned1 = dcmecef2ned(station1(1)， station1(2));Ned2 = dcmecef2ned(station2(1)， station2(2));covcon(1) = coverageConfig(radar1,lla2ecef(station1)，四元数(ned1，“rotmat”，“帧”));covcon(2) = coverageConfig(radar2,lla2ecef(station2)，四元数(ne2，“rotmat”，“帧”));covcon plotCoverage(查看器,“ECEF”）;
您首先生成超过5小时的星座状态的整个历史。然后，模拟雷达探测并在循环中生成轨迹。
satscene。停止时间=卫星场景。开始时间+小时(5);satscene。SampleTime = 10;numSteps = ceil(seconds(satscene.)- satscene.StartTime)/satscene.SampleTime);在模拟过程中获得星座位置和速度Plats = repmat(.．.结构(“PlatformID”0,“位置”，[0 0 0]，“速度”， [0 0 0])，.．.numSteps 40);为I =1:数字(星座)[pos, vel] = state(星座(I)，“CoordinateFrame”，“ECEF”）;为j = 1: numSteps公寓(j,我)。Position = pos(:，j)';我的公寓(j)。速度= vel(:，j)';我的公寓(j)。PlatformID = i;结束结束初始化跟踪和跟踪日志confTracks = objectTrack.empty(0,1);trackLog = cell(1,numSteps);%初始化雷达图radarplt = helperRadarPlot(fov);为可重复的结果设置随机种子S = rng;rng (2020);Step = 0;而step < numSteps time = step*satscene.SampleTime;Step = Step + 1;生成雷达跟踪星座的探测。%加工链targets1 = assemblerarinputs (station1, plats(step，:));[dets1,numdets1] = radar1(targets1, time);dets1 = addMeasurementParams(dets1,numdets1,station1);targets2 = assemblerarinputs (station2, plats(step，:));[dets2, numdets2] = radar2(targets2, time);dets2 = addMeasurementParams(dets2, numdets2, station2);检测= [dets1;dets2];updateRadarPlots(radarplt,targets1, targets2,dets1, dets2)生成和更新曲目detectableInput = isdetected(跟踪器，时间，covcon);如果~isempty(detects) || isLocked(tracker) [confTracks，~，~，info] = tracker(detects,time,detectableInput);结束trackLog{step} = confTracks;%更新查看器plotPlatform(观众、公寓(一步,:)“ECEF”，“颜色”，[1 0 0]，“线宽”1);plotDetection(查看器,检测,“ECEF”）;confTracks plotTrack(查看器,“ECEF”，“颜色”，[0 10 0]，“线宽”3);结束
上图显示了每个雷达视角下的轨道(蓝点)和探测(红圈)。
恢复之前的随机种子状态rng(年代);
图;快照(观众);
经过5个小时的跟踪，大约一半的星座被成功跟踪。保持部分轨道覆盖的轨道具有挑战性，因为卫星通常在这种配置下长时间不被发现。在这个例子中，只有两个雷达站。更多的跟踪站有望产生更好的跟踪性能。分配指标，通过比较真实对象和轨道来评估跟踪性能，如下所示。
%显示分配指标truthIdFcn = @(x)[x. platformid];tam = trackAssignmentMetrics(“DistanceFunctionFormat”，“自定义”，.．.“AssignmentDistanceFcn”@distanceFcn,.．.“DivergenceDistanceFcn”@distanceFcn,.．.“TruthIdentifierFcn”truthIdFcn,.．.“AssignmentThreshold”, 1000,.．.“DivergenceThreshold”, 2000);为我= 1:numSteps在第i次更新时提取跟踪器和ground truthtracks = trackLog{i};真理= plats(i，:);根据轨道和真相提取分配指标的摘要。[trackAM,truthAM] = tam(轨道，真理);结束
显示每个单独记录的真值对象的累积度量结果= truthMetricsTable(tam);结果(:,{“TruthID”，“AssociatedTrackID”，“BreakLength”，“EstablishmentLength”})
ans =40×4表TruthID AssociatedTrackID BreakLength EstablishmentLength _______ _________________ ___________ ___________________ 1 52 5 232 24 594 3 4 0 2 56 437 55 11 492 5 18 0 6 48 513 811 7 21 0 8 27 661 9 0 1221 10 0 0 1339 37 12 0 1056 1504 11 43 13南0 1800 1800 15南南0 0 1800 16南⋮0 1800
上表列出了发射星座中的40颗卫星，并显示了具有相关轨道id的跟踪卫星。轨道ID为NaN表示模拟结束时卫星未被跟踪。这意味着卫星的轨道没有穿过两个雷达中的一个的视野，或者卫星的轨道已经下落。跟踪器可能会由于初始探测数量不足而掉线，这导致估计有很大的不确定性。另外，如果卫星没有很快被重新检测到，跟踪器可以丢弃跟踪，这样缺乏更新就会导致偏离并最终删除。
总结
在本例中，您已经学习了如何使用satelliteScenario对象，以从TLE文件中导入轨道信息。您使用SGP4传播了卫星轨迹，并使用卫星场景查看器可视化了场景。您学习了如何使用传感器融合和跟踪工具箱™中的雷达和跟踪器模型来建模空间监视雷达跟踪系统。所构建的跟踪系统可以使用低保真模型预测每个卫星的估计轨道。
支持功能
initKeplerUKF初始化Unscented卡尔曼滤波器使用开普勒运动模型。集α= 1,β= 0，和卡巴= 0，确保无气味滤波器在较长的预测周期内具有鲁棒性。
函数filter = initKeplerUKF(detection) radarsphmeas = detection. measurement;[x, y, z] = sph2cart(deg2rad(radarsphmeas(1))，deg2rad(radarsphmeas(2))，radarsphmeas(3));雷达图像= [x y z];Recef2radar = detection.MeasurementParameters.Orientation;ecefmeas = detection.MeasurementParameters.OriginPosition + radarcartmeas*Recef2radar;%这相当于:% Ry90 = [0 0 -1;0 10 0;10 0 0];%帧旋转90度绕y轴% nedmeas(:) = Ry90' * radarcartmeas(:);% ecefmeas = lla2ecef(lla) + nedmeas * dcmecef2ned(lla(1)，lla(2));initState = [ecefmeas(1);0;ecefmeas (2);0;ecefmeas (3);0);Sigpos = 2;% mSigvel = 0.5;% m / s ^ 2filter = trackingUKF(@keplerorbit，@cvmeas,initState，.．.“α”, 1“β”0,“卡巴”0,.．.“StateCovariance”， diag(repmat([1000, 10000].^2,1,3))，.．.“ProcessNoise”、诊断接头(repmat ([sigpos, sigvel]。^ 2,1,3)));结束函数状态= keplerorbit(状态，dt)% keplerorbit执行数值积分来预测的状态%开普勒体。状态是[x;vx;y;vy;z;vz]%龙格-库塔四阶积分法:K1 =开普勒(状态);K2 =开普勒(状态+ dt*k1/2);K3 =开普勒(状态+ dt*k2/2);K4 =开普勒(状态+ dt*k3);状态=状态+ dt*(k1+2*k2+2*k3+k4)/6;函数Dstate =开普勒(状态)x =状态(1，:);Vx = state(2，:);y =状态(3:);Vy =状态(4，:);: z =状态(5日);Vz = state(6，:);Mu = 398600.4405*1e9;% m^3 s^ 2Omega = 7.292115e-5;% rad /秒R =范数([x y z]);G = /r^2;%坐标在非区间坐标系中，考虑科里奥利%和向心加速度Ax = -g*x/r + 2* *vy + ^2*x;Ay = -g*y/r - 2* *vx + ^2*y;Az = -g*z/r;Dstate = [vx;ax;vy;ay;vz;az];结束结束
isDetectable在示例中用于确定在给定时间可以检测到哪些轨道。
函数detectInput = isdetected(跟踪器，时间，covcon)如果~isLocked(跟踪器)detectInput = 0 (0,1，“uint32”）;返回结束tracks = tracker.predictTracksToTime(“所有”、时间);如果isempty(轨道)detectInput =零(0,1，“uint32”）;其他的alltrackid = [tracks.TrackID];isdetected =零(数字(轨道)，数字(covcon)，“逻辑”）;为I = 1: number (tracks) track = tracks(I);Pos_scene = track。状态([1 3 5]);为J =1: config = covcon(J);%到传感器框的旋转位置:d_scene = pos_scene(:) - config.Position(:);Scene2sens = rotmat(配置。取向,“帧”）;D_sens = scene2sens*d_scene(:);(阿兹,el) = cart2sph (d_sens (1) d_sens (2), d_sens (3));如果abs(rad2deg(az)) <= config.FieldOfView(1)/2 && abs(rad2deg(el)) < config.FieldOfView(2)/2 isdetected (i,j) = true;其他的isdetected (i,j) = false;结束结束结束detectInput = alltrackid(any(isdetected,2))';结束结束
assembleRadarInput用于构造每个雷达体帧中的星座位姿。这是图中描述的第一步。
函数目标= assemblerarinputs(站，星座)对于星座中的每一颗卫星，推导其相对于的姿态雷达框架。%的输入:% - station:雷达地面站LLA矢量星座:包含ECEF位置的结构数组%和每颗卫星的ECEF速度%输出:目标:包含每个结构的姿态的结构数组%卫星相对于雷达，以本地表示地面雷达帧(NED)输出的模板结构，包含所需的所有字段%采用雷达步进法targetTemplate = struct(.．.“PlatformID”0,.．.“ClassID”0,.．.“位置”, 0(1、3).．.“速度”, 0(1、3).．.“加速”, 0(1、3).．.“定位”四元数(0,0,0),.．.“AngularVelocity”, 0(1、3).．.“维度”结构(.．.“长度”0,.．.“宽度”0,.．.“高度”0,.．.“OriginOffset”， [0 0 0])，.．.“签名”{{rcsSignature}});首先填充当前卫星ECEF姿态targetpose = repmat(targetTemplate, 1，数字(星座));为i = 1:元素个数(星座)targetPoses(我)。位置=星座(i).位置;targetPoses(我)。速度=星座(i).速度;targetPoses(我)。PlatformID =星座(i).PlatformID;%方向和角速度为空，假设卫星到%是具有统一RCS的点目标结束然后根据地面站位置，推导出ECEF中的雷达位姿Recef2station = dcmecef2ned(station(1)， station(2));radarPose。方位=四元数(Recef2station，“rotmat”，“帧”）;radarPose。位置= lla2ecef(站);radarPose。速度= 0 (1,3);radarPose。AngularVelocity = 0 (1,3);最后，取差值并将每个矢量旋转到地面站% NED轴目标=目标姿势;为i=1: nummel (targetpose) thisTgt = targetpose (i);pos = Recef2station*(thisTgt.Position(:) - radarPose.Position(:));vel = Recef2station*(thisTgt.Velocity(:) - radarPose.Velocity(:)) - cross(radarPose.AngularVelocity(:)，pos(:));angVel = thisTgt.AngularVelocity(:) - radarPose.AngularVelocity(:);orient = radarPose。Orientation' * thisTgt.Orientation;存储到目标结构数组目标(i).Position(:) = pos;目标(i).Velocity(:) = vel;目标(i).AngularVelocity(:) = angVel;目标(i)。方位=东方;结束结束
addMeasurementParam实现雷达链过程中的步骤2，如图所示。
函数dets = addMeasurementParams(dets, numdets, station)在测量参数中添加雷达站信息Recef2station = dcmecef2ned(station(1)， station(2));为i = 1: numdets侦破{我}.MeasurementParameters。OriginPosition = lla2ecef(站);依据{我}.MeasurementParameters。IsParentToChild = true;% parent = ecef, child =雷达依据{我}.MeasurementParameters。Orientation = dets{i}. measurementparameters。定向' * Recef2station;结束结束
distanceFcn与分配指标一起使用，以评估跟踪分配。
函数d = distanceFcn(轨道，真相)估计=轨道。状态([1 3 5 2 4 6]);true = [true . position (:);truth.Velocity (:));Cov =轨道。StateCovariance([1 3 5 2 4 6]， [1 3 5 2 4 6]);D =(估计-正确)' / cov *(估计-正确);结束
参考
[1]斯里达兰、拉马斯瓦米、安东尼奥·f·彭萨编。空间监测展望．麻省理工学院出版社，2017年。
您也可以从以下列表中选择一个网站: