利用地面雷达探测和跟踪LEO卫星星座
本例展示了如何导入卫星星座的双线元素(TLE)文件，模拟星座的雷达探测，并跟踪星座。
收集和维护绕地球轨道运行的空间物体的目录是空间监视的关键任务。这项任务包括几个过程:检测和识别新物体并将其添加到目录中，更新目录中已知物体的轨道，跟踪其整个生命周期中的轨道变化，预测在大气层中的再入。在这个例子中，我们正在研究如何检测和跟踪新的卫星，并将它们添加到一个目录中。
为保证太空安全运行，防止与其他卫星或已知碎片发生碰撞，对新发射卫星进行正确探测和跟踪十分重要。航天机构通常共享发射前的信息，这些信息可用于选择搜索策略。常用的低地球轨道(LEO)卫星搜索策略是由围栏型雷达系统组成。围栏式雷达系统在空间中搜索有限体积的空间，并在卫星通过其视野时探测到它们。该策略可以快速探测和跟踪新发射的星座。［１］
从TLE文件中导入卫星星座
双线元集是保存卫星轨道信息的常用数据格式。您可以使用satelliteScenario对象导入TLE文件中定义的卫星轨道。缺省情况下，导入的卫星轨道使用SGP4轨道传播算法进行传播，该算法对LEO目标具有良好的精度。在本例中，这些轨道提供了地面真相，以测试雷达跟踪系统探测新发射卫星的能力。
创建一个卫星场景satscene = satelliteScenario;从TLE文件中添加卫星。tleFile =“leoSatelliteConstellation.tle”；constellation = satellite(satscene, tleFile);
使用卫星场景查看器来可视化星座。
玩(satscene);
模拟合成探测和轨迹星座
空间监视雷达建模
定义两个用扇形雷达波束观测太空的观测站。风扇穿过卫星轨道，以最大限度地增加探测次数。位于北美的雷达站形成了东西篱笆。
第一站在LLA内Station1 = [48 -80 0];第二站在LLA坐标Station2 = [50 -117 0];
每个站都配备了一个雷达，它是用一个fusionRadarSensor对象。为了探测低轨道范围内的卫星，雷达有以下要求:
在2000公里外探测到一个10 dBsm的物体
在2000公里范围内以100米的精度水平和垂直分辨物体
方位角120度，仰角30度的视野
仰望太空
创建扇形单台雷达fov = [40] 120;;radar1 = fusionRadarSensor (1,.．.“UpdateRate”, 0.1,.．.10秒“ScanMode”，“没有扫描”，.．.“MountingAngles”, 90 0] [0,.．.查找“FieldOfView”fov,.．.度“ReferenceRange”2000年e3,.．.米“RangeLimits”[0 2000年e3),.．.米“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定义跟踪筛选器。
%定义跟踪器追踪= trackerJPDA (“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));ne2 = dcmecef2ned(station2(1)， station2(2));covcon (1) = coverageConfig (radar1 lla2ecef (station1),四元数(ned1,“rotmat”，“帧”));covcon (2) = coverageConfig (radar2 lla2ecef (station2),四元数(ned2,“rotmat”，“帧”));covcon plotCoverage(查看器,“ECEF”）;
首先生成5小时内星座状态的全部历史。然后，模拟雷达探测并在循环中生成轨迹。
satscene。StopTime = satscene。开始时间(5)+小时;satscene。SampleTime = 10;numSteps =装天花板(秒(satscene。StopTime - satscene.StartTime) / satscene.SampleTime);在模拟过程中获取星座位置和速度。公寓= repmat (.．.结构(“PlatformID”0,“位置”(0 0 0),“速度”， [0 0 0])，.．.numSteps 40);为I =1:numel(星座)[pos, vel] = states(星座(I)，“CoordinateFrame”，“ECEF”）;为j = 1: numSteps公寓(j,我)。位置= pos (:, j) ';我的公寓(j)。速度=韦尔(:,j)”;我的公寓(j)。PlatformID =我;结束结束初始化跟踪和跟踪日志confTracks = objectTrack.empty (0,1);trackLog =细胞(1、numSteps);%初始化雷达图radarplt = helperRadarPlot (fov);为可复制的结果设置随机种子s =提高;rng (2020);一步= 0;而step < numSteps time = step*satscene.SampleTime;Step = Step + 1;生成雷达跟踪星座的探测。%处理链targets1 = assembleRadarInputs(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(detections) || isLocked(tracker) [confTracks，~，~，info] = tracker(detections,time,detectableInput);结束trackLog{一步}= confTracks;%更新查看器plotPlatform(观众、公寓(一步,:)“ECEF”，“颜色”(1 0 0),“线宽”1);plotDetection(查看器,检测,“ECEF”）;confTracks plotTrack(查看器,“ECEF”，“颜色”(0 1 0),“线宽”3);结束
上图显示了每个雷达的轨道(蓝点)和探测结果(红圈)。
恢复以前的随机种子状态rng(年代);
图;快照(观众);
经过5小时的跟踪，大约一半的星座被成功跟踪。保持部分轨道覆盖的轨道是具有挑战性的，因为卫星经常在这种配置下长时间不被发现。在这个例子中，只有两个雷达站。额外的跟踪站有望产生更好的跟踪性能。分配指标，评估跟踪性能之间的真实对象和轨迹之间的比较，如下所示。
%显示分配指标truthIdFcn = @ (x) [x.PlatformID];tam = trackAssignmentMetrics (“DistanceFunctionFormat”，“自定义”，.．.“AssignmentDistanceFcn”@distanceFcn,.．.“DivergenceDistanceFcn”@distanceFcn,.．.“TruthIdentifierFcn”truthIdFcn,.．.“AssignmentThreshold”, 1000,.．.“DivergenceThreshold”, 2000);为我= 1:numSteps在第i次更新跟踪器时提取跟踪器和实值我跟踪= trackLog {};真理=公寓(我:);根据轨迹和事实提取分配度量的摘要。[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的跟踪卫星。值为NaN的跟踪ID表示在模拟结束时未跟踪卫星。这要么意味着卫星的轨道没有穿过两个雷达中的一个的视野，要么意味着卫星的轨迹被丢弃了。跟踪器可能会由于初始检测数量不足而丢弃一个跟踪，这导致估计有很大的不确定性。另外，如果没有及时重新探测到卫星，跟踪器也会丢弃跟踪，这样就会导致卫星偏离并最终被删除。
总结
在本例中，您已经学习了如何使用satelliteScenario对象从Aerospace Toolbox™导入TLE文件中的轨道信息。您使用SGP4传播了卫星轨迹，并使用卫星场景查看器可视化了场景。您学习了如何使用传感器融合和跟踪工具箱™中的雷达和跟踪器模型来建模空间监视雷达跟踪系统。所构建的跟踪系统可以使用低保真度模型预测每颗卫星的估计轨道。
支持功能
initKeplerUKF使用开普勒运动模型初始化无气味卡尔曼滤波器。集α= 1,β= 0,卡巴= 0，保证无气味滤波器在较长预测周期内的鲁棒性。
函数filter = initKeplerUKF(检测)radarsphmeas =检测.测量;[x, y, z] = sph2cart(deg2rad(radarsphmeas(1))，deg2rad(radarsphmeas(2))，radarsphmeas(3));Radarcartmeas = [x y z];Recef2radar = detection.MeasurementParameters.Orientation;ecefmeas = detection.MeasurementParameters.OriginPosition + radarcartmeas*Recef2radar;%这相当于:% Ry90 = [0 0 -1;0 1 0;1 0 0];%帧绕y轴旋转90度% 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 ^ 2过滤器= trackingUKF (@keplerorbit @cvmeas initState,.．.“α”, 1“β”0,“卡巴”0,.．.“StateCovariance”、诊断接头(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 = kepler(state + dt*k2/2);K4 = kepler(state + dt*k3);状态=状态+ dt*(k1+2*k2+2*k3+k4)/6;函数dstate =开普勒(状态)x =状态(1:);vx =状态(2:);y =状态(3:);v =状态(4:);: z =状态(5日);: vz =状态(6日);μ= 398600.4405 * 1 e9;% s m ^ 3 ^ 2ω= 7.292115 e-5;% rad /秒R =范数([x y z]);g =μ/ r ^ 2;%坐标在非间隙坐标系中，考虑科里奥利%和向心加速度Ax = -g*x/r + 2* *vy +²*x;Ay = -g*y/r - 2* *vx +²*y;阿兹= r - g * z /;dstate = [vx;斧子;v;是的;vz; az);结束结束
isDetectable在示例中用于确定在给定时间可以检测到哪些轨迹。
函数detectInput = isDetectable(跟踪、时间covcon)如果~isLocked(跟踪器)detectInput = 0 (0,1，“uint32”）;返回结束跟踪= tracker.predictTracksToTime (“所有”、时间);如果isempty(tracks) detectInput = 0 (0,1，“uint32”）;其他的alltrackid = [tracks.TrackID];isDetectable = 0(元素个数(跟踪),元素个数(covcon),“逻辑”）;为I = 1:numel(tracks) track = tracks(I);pos_scene =。状态([1 3 5]);为J =1:numel(covcon) 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;其他的isDetectable (i, j) = false;结束结束结束detectInput = alltrackid(任何(isDetectable 2) ';结束结束
assembleRadarInput用于构造雷达体帧中的星座位姿。这是图中描述的第一步。
函数目标= assembly eradarinputs(站，星座)对于星座中的每一颗卫星，推导其相对于的姿态。雷达帧。%的输入:% - station:雷达地面站的LLA向量星座:包含ECEF位置的结构数组%和每颗卫星的ECEF速度%输出:% - targets:包含每个结构的姿态的结构数组%卫星相对于雷达，表示为本地%地面雷达帧(NED)输出的模板结构，其中包含所需的所有字段%的雷达步进法targetTemplate =结构(.．.“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, numel(constellation));为i = 1:元素个数(星座)targetPoses(我)。位置=(我).Position星座;targetPoses(我)。速度=(我).Velocity星座;targetPoses(我)。PlatformID =(我).PlatformID星座;方向和角速度为零，假设卫星到%是具有统一RCS的点目标结束然后根据地面站位置推导出ECEF中的雷达位姿。Recef2station = dcmecef2ned(station(1)， station(2));radarPose。取向=四元数(Recef2station,“rotmat”，“帧”）;radarPose。位置= lla2ecef(站);radarPose。速度= 0(1、3);radarPose。AngularVelocity = 0(1、3);最后，取差值并将每个矢量旋转到地面站。% NED轴目标= targetPoses;为i=1: numl (targetpose) thisTgt = targetpose (i);pos = Recef2station*(thisTgt.Position(:) - radarPose.Position(:));vel = Recef2station*(thisTgt.Velocity(:) - radarPose.Velocity(:)) - cross(radarPose.AngularVelocity(:)，pos(:));angVel = thisTgt.AngularVelocity(:) - radarpo . angularvelocity (:);东方= radarPose。取向的* thisTgt.Orientation;%存储到目标结构数组中目标(我).Position (:) = pos;目标(我).Velocity(:) =韦尔;目标(我).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 = radar依据{我}.MeasurementParameters。取向=侦破{我}.MeasurementParameters。取向的* Recef2station;结束结束
distanceFcn与分配指标一起使用，以评估跟踪分配。
函数d = distanceFcn(track, truth) estimate = track。状态([1 3 5 2 4 6]);true = [true . position (:);truth.Velocity (:));x =。StateCovariance([1 3 5 2 4 6]， [1 3 5 2 4 6]);D =(估计-真)' / cov *(估计-真);结束
参考
[1] Sridharan, Ramaswamy和Antonio F. Pensa编。空间监视的展望．麻省理工学院出版社,2017年。
另请参阅
对象
功能
相关的话题
您也可以从以下列表中选择网站: