用地面雷达探测和跟踪LEO卫星星座
该示例演示了如何导入一个卫星星座的双线元素(Two-Line Element, TLE)文件，模拟该星座的雷达探测，并跟踪该星座。
在空间监视中，对环绕地球的空间物体进行统计和维护是至关重要的任务。这项任务包括几个过程:探测和识别新的物体并将它们添加到目录中，更新目录中已知物体的轨道，跟踪它们在整个生命周期中的轨道变化，并预测在大气层中的再入。在这个例子中，我们将学习如何探测和跟踪新的卫星，并将它们添加到目录中。
为了保证太空中的安全运行，防止与其他卫星或已知碎片发生碰撞，正确探测和跟踪新发射的卫星非常重要。航天机构通常会共享发射前的信息，这些信息可以用来选择搜索策略。通常采用的是由栅栏型雷达系统组成的近地轨道卫星搜索策略。围栏型雷达系统在空间中搜索有限的空间，并在卫星通过其视场时进行探测。该策略可以快速检测和跟踪新发射的星座。［1］
从TLE文件导入卫星星座
双线元组是保存卫星轨道信息的常用数据格式。你可以使用satelliteScenario对象导入在TLE文件中定义的卫星轨道。默认情况下，导入的卫星轨道使用SGP4轨道传播算法进行传播，该算法为LEO目标提供了良好的精度。在这个例子中，这些轨道提供了地面真相，以测试雷达跟踪系统检测新发射卫星的能力。
%创建一个卫星场景satscene = satelliteScenario;从TLE文件中添加卫星。tleFile =“leoSatelliteConstellation.tle”；星座=卫星(satscene, tleFile);
使用卫星场景查看器来可视化星座。
玩(satscene);
模拟合成探测和轨道星座
模拟空间监视雷达
用扇形雷达波束定义两个空间站。风扇穿过卫星轨道，以最大限度地增加探测次数。位于北美的雷达站形成了一个东西栅栏。
LLA第一站坐标[48 -80 0];LLA第二站坐标[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金宝app支持功能的实现细节。
在第二步中，将所需信息添加到检测对象，以便跟踪器可以使用ECEF状态进行操作。您使用MeasurementParameters属性来实现该目的，如addMeasurementParams金宝app支持功能。
定义一个跟踪器
上面定义的雷达模型输出检测。要估计卫星轨道，你需要一个跟踪器。传感器融合和跟踪工具箱™提供了多种多目标跟踪器。在本例中，您选择联合概率数据关联(Joint Probabilistic Data Association, JPDA)跟踪器，因为它提供了跟踪性能和计算成本的良好平衡。
您需要为跟踪器定义一个跟踪过滤器。你可以使用一个比SGP4低保真度的模型，比如运动方程的开普勒积分，来跟踪卫星。通常，目标运动模型保真度的不足是通过测量更新和在滤波器中加入过程噪声来补偿的。支持函数金宝appinitKeplerUKF定义跟踪过滤器。
%定义跟踪器追踪= 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(constellation) [pos, vel] = states(constellation(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，时间);dets1 = addMeasurementParams (dets1 numdets1 station1);targets2 = assembleRadarInputs(station2, plats(step，:))); / /设定目标[dets2, numdets2] = radar2(targets2，时间);dets2 = addMeasurementParams(dets2, numdets2, station2);检测= [dets1;dets2];updateRadarPlots(radarplt,targets1, targets2,dets1, dets2)%生成和更新轨道detectableInput = is可检测的(tracker,time, covcon);如果~isempty(detecttions) || isLocked(tracker) [confTracks，~，~，info] = tracker(detecttions,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 57 5 232 45 680 598 3 29 0 664 4 61 11 492 5 18 0 436 6 54 5 807 513 22日0 8 47 6 675 9 42 0 1221 10 56 0 1500 11 49 1339 12 40 0 1056 13南0 1800 1800 15南南0 0 1800 16南⋮0 1800
上表列出了发射星座中的40颗卫星，并显示了具有相关轨道id的跟踪卫星。轨道ID值NaN表示在模拟结束时卫星没有被跟踪。这要么意味着卫星的轨道没有通过两个雷达中的一个的视场，要么意味着卫星的轨道已经被丢弃。由于初始检测次数不足，跟踪器可能会丢弃跟踪，这导致估计的不确定性很大。另一种情况是，如果卫星没有很快被重新探测到，跟踪器就会放弃跟踪，这样缺乏更新就会导致偏离并最终被删除。
总结
在这个例子中，你已经学会了如何使用satelliteScenario对象，以从TLE文件导入轨道信息。您使用SGP4传播卫星轨迹，并使用卫星场景查看器可视化场景。您学习了如何使用传感器融合和跟踪工具箱™中的雷达和跟踪器模型来建模空间监视雷达跟踪系统。所构建的跟踪系统可以利用低保真度模型预测每颗卫星的估计轨道。
金宝app支持功能
initKeplerUKF使用开普勒运动模型初始化Unscented卡尔曼滤波器。
函数filter = initKeplerUKF(detection) radarsphmeas = detection. measurement;[x, y, z] = sph2cart(deg2rad(radarsphmeas(1))，deg2rad(radarsphmeas(2))，radarsphmeas(3));雷达carmeas = [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,．．.“StateCovariance”、诊断接头(repmat((1000、10000)。^ 2,1,3)),．．.“ProcessNoise”、诊断接头(repmat ([sigpos, sigvel]。^ 2,1,3)));结束函数国家= keplerorbit(状态,dt)%开普勒轨道执行数值积分来预测的状态%开普勒的身体。状态为[x;vx;y;vy;z;vz]%龙格-库塔四阶积分法:k1 =开普勒(状态);K2 =开普勒(状态+ dt*k1/2);K3 = kepler(state + dt*k2/2);K4 = kepler(state + dt*k3);State = State + 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 = norm([x y z]);g =μ/ r ^ 2;坐标在非惯性坐标系中，考虑科里奥利和向心加速度Ax = -g*x/r + 2* vy + ^2*x;Ay = -g*y/r - 2* vx + ^2*y;阿兹= r - g * z /;dstate = [vx;斧子;v;是的;vz; az);结束结束
isDetectable用于确定在给定时间内哪些轨迹是可检测的。
函数detectInput = isDetectable(跟踪、时间covcon)如果~isLocked(tracker) 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(el)) <= config.FieldOfView(1)/2 && abs(rad2deg(el)) <= config.FieldOfView(2)/2 isdetected (i,j) = true;其他的isDetectable (i, j) = false;结束结束结束detectInput = alltrackid(任何(isDetectable 2) ';结束结束
assembleRadarInput用于构造每个雷达机体框架中的星座位姿。这是图中描述的第一步。
函数targets = assembleRadarInputs(车站，星座)对星座中的每一颗卫星，求其相对于的位姿%雷达框架。%的输入:% - 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: numel(targetpose) thisTgt = targetpose (i);pos = Recef2station*(thisTgt.Position(:) - radarPose.Position(:)));vel = Recef2station*(thisTgt.Velocity(:) - radarPose.Velocity(:)) - cross(radarPose.AngularVelocity(:));angVel = thisTgt.AngularVelocity(:) - radarPose.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)状态([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 =(估计值- true)' / cov *(估计值- true);结束
参考
Sridharan, Ramaswamy和Antonio F. Pensa编。空间监视展望．麻省理工学院出版社,2017年。
你也可以从以下列表中选择一个网站:
