用地面雷达探测和跟踪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毫米的物体
在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”， [0200e3]，…米“ReferenceRCS”10…dBsm“HasFalseAlarms”假的,…“HasNoise”,真的,…“HasElevation”,真的,…“AzimuthResolution”, 0.03,…度“ElevationResolution”, 0.03,…度“RangeResolution”, 2000,…M %精度~= 2000 * 0.05 (M)“DetectionCoordinates”，“球形传感器”，…“TargetReportFormat”，“检测”）;Radar2 = clone(radar1);radar2。SensorIndex = 2;
雷达处理链
在本例中，为了使雷达跟踪链正常运行，进行了多次坐标转换和轴向转换。下图说明了上面定义的输入是如何转换并传递给雷达的。
在第一步，你计算每个卫星姿态在本地雷达站NED轴。您可以通过首先获得地面站ECEF姿态并将卫星位置和速度转换为ECEF坐标来实现这一点。雷达输入是将卫星位姿与地面站位姿的差值旋转到地面站局地内NED轴上得到的。看到assembleRadarInputs金宝app实现细节的支持功能。
在第二步中，向检测对象添加所需的信息，以便跟踪器可以在ECEF状态下操作。你使用MeasurementParameters属性来实现这一目的，如addMeasurementParams金宝app支持功能。
定义跟踪器
上面定义的雷达模型输出检测。要估计卫星轨道，你需要一个跟踪器。传感器融合和跟踪工具箱™提供了各种多目标跟踪器。在本例中，您选择联合概率数据关联(Joint Probabilistic Data Association, JPDA)跟踪器，因为它提供了跟踪性能和计算成本的良好平衡。
您需要为跟踪器定义跟踪过滤器。您可以使用比SGP4更低保真度的模型，例如运动方程的开普勒积分，来跟踪卫星。通常，通过测量更新和在滤波器中加入过程噪声来补偿目标运动模型中保真度的不足。配套功能金宝appinitKeplerUKF定义跟踪过滤器。
%定义跟踪器tracker = 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。StartTime + hours(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:numel(constellation) [pos, vel] = states(constellation(I)，“CoordinateFrame”，“ECEF”）;为j = 1: numSteps公寓(j,我)。Position = pos(:，j)';我的公寓(j)。Velocity = vel(:，j)';我的公寓(j)。PlatformID = i;结束结束%初始化音轨和音轨日志confTracks = objectTrack.empty(0,1);trackLog = cell(1,numSteps);%初始化雷达图radarplt = helperRadarPlot(fov);设置随机种子以获得可重复的结果S = rng;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 = assembleRadarInputs(station2, plats(step，:));[dets2, numdets2] = radar2(targets2, time);dets2 = addMeasurementParams(dets2, numdets2, station2);检测= [dets1;dets2];updaterearplots (radarplt,targets1, targets2,dets1, dets2)%生成和更新音轨detectableInput = isDetectable(跟踪器，时间，covcon);如果~isempty(detections) || isLocked(tracker) [confTracks，~，~，info] = tracker(detections,time,detectableInput);结束trackLog{step} = 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次更新跟踪器时提取跟踪器和地面真相tracks = trackLog{i};真理= plats(i，:);根据轨迹和事实提取任务度量标准的摘要[trackAM,truthAM] = tam(tracks, truths);结束
显示每个单独记录的真值对象的累积度量results = 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表示在模拟结束时不跟踪该卫星。这要么意味着卫星的轨道没有通过两个雷达中的一个的视野，要么意味着卫星的轨迹已经被丢弃。由于初始检测数量不足，跟踪器可能会丢失跟踪，从而导致估计的不确定性很大。另一种情况是，如果没有足够快地重新探测到卫星，跟踪器可以放弃跟踪，这样缺乏更新就会导致偏离并最终删除。
总结
在本例中，您已经学习了如何使用satelliteScenario对象从TLE文件中导入轨道信息。您使用SGP4传播卫星轨迹，并使用satellite scenario Viewer可视化场景。您学习了如何使用传感器融合和跟踪工具箱™中的雷达和跟踪器模型来模拟空间监视雷达跟踪系统。所构建的跟踪系统可以利用低保真度模型预测每颗卫星的估计轨道。
金宝app支持功能
initKeplerUKF使用开普勒运动模型初始化一个Unscented卡尔曼滤波器。集α= 1,β= 0，且卡巴= 0以确保无气味过滤器在长预测周期内的鲁棒性。
函数filter = initKeplerUKF(detection) radarsphmeas = detection. measurement;[x, y, z] = sph2cart(deg2rad(radarsphmeas(1))，deg2rad(radarsphmeas(2))，radarsphmeas(3));radarcartmeans = [x y z];Recef2radar = detection. measurement . parameters . orientation;ecefmeas = detection.MeasurementParameters.OriginPosition + radarcartmeas*Recef2radar;这相当于:% Ry90 = [0 0 -1;0 0 0;[1 0 0];帧绕y轴旋转90度% nemeas (:) = r90 ' * radarcartmeas(:);% ecefmeas = lla2ecef(lla) + nedmeas * dcmecef2ned(lla(1)，lla(2));initState = [ecefmeans (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)));结束函数State = keplerorbit(State,dt)开普勒轨道执行数值积分来预测的状态%开普勒天体。状态为[x;vx;y;vy;z;vz]%龙格-库塔四阶积分法:K1 =开普勒(状态);K2 = kepler(state + 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 =kepler(state) x =state(1，:);Vx = state(2，:);y =状态(3:);Vy = state(4，:);: z =状态(5日);Vz = state(6，:);Mu = 398600.4405*1e9;% m^3 s^-2ω = 7.292115e-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;Az = -g*z/r;Dstate = [vx;ax; vz;az];结束结束
isDetectable在示例中用于确定在给定时间可检测到哪些轨迹。
函数detectInput = isDetectable(跟踪器，时间，covcon)如果~isLocked(tracker) detectInput = 0 (0,1，“uint32”）;返回结束tracks = tracker.predictTracksToTime(“所有”、时间);如果isempty(tracks) detectInput = 0 (0,1，“uint32”）;其他的alltrackid = [track . trackid];isDetectable = 0 (numel(tracks)，numel(covcon)，“逻辑”）;为I = 1:numel(tracks) track = tracks(I);Pos_scene = track。State([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 isDetectable(i,j) = true;其他的isDetectable(i,j) = false;结束结束结束detectInput = alltrackid(any(isDetectable,2))';结束结束
assembleRadarInput用于构造各雷达体框架中的星座位姿。这是图中描述的第一步。
函数目标= assembleRadarInputs(工作站，星座)对于星座中的每颗卫星，推导其相对于的位姿%雷达框架。%的输入:% - station:雷达地面站的LLA矢量星座:包含ECEF位置的结构体数组%和每颗卫星的ECEF速度%输出:% - targets:包含每个目标姿态的结构体数组%卫星相对于雷达，用本地表示%地面雷达帧(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, numel(星座));为i = 1:元素个数(星座)targetPoses(我)。Position =星座(i).Position;targetPoses(我)。速度=星座(i).速度;targetPoses(我)。PlatformID = constellation(i).PlatformID;方向和角速度为空，假设卫星为%是具有统一RCS的点目标结束然后根据地面站位置推导出ECEF中的雷达位姿Recef2station = dcmecef2ned(站(1)，站(2));radarPose。方位=四元数(接收站，“rotmat”，“帧”）;radarPose。位置= lla2ecef(站);radarPose。速度= 0 (1,3);radarPose。AngularVelocity = 0 (1,3);最后，取差值并将每个矢量旋转到地面站% NED轴target = targetpose;为i=1: numel(targetpose);pos = Recef2station*(thisTgt.Position(:) - radarPose.Position(:));level = Recef2station*(thisTgt.Velocity(:) - radarPose.Velocity(:)) - cross(radarPose.AngularVelocity(:)，pos(:)));angVel = thisTgt.AngularVelocity(:) - radarPose.AngularVelocity(:);orient = radarPose。* thisTgt.Orientation;%存储到目标结构数组目标(i).Position(:) = pos;目标(i).速度(:)=水平;目标(i).AngularVelocity(:) = angVel;目标(i)。方位=东方;结束结束
addMeasurementParam如图所示，实现雷达链过程中的步骤2。
函数dets = addMeasurementParams(dets, numdets, station)%在测量参数中增加雷达站信息Recef2station = dcmecef2ned(站(1)，站(2));为i = 1: numdets侦破{我}.MeasurementParameters。OriginPosition = lla2ecef(站);依据{我}.MeasurementParameters。IsParentToChild = true;% parent = ecef, child = radar依据{我}.MeasurementParameters。取向= dets{i}. measurementparameters。方位*接收站;结束结束
distanceFcn与分配度量一起使用，以评估跟踪分配。
函数d = distanceFcn(轨迹，真值)estimate =轨迹State([1 3 5 2 4 6]);true = [true . position (:);]truth.Velocity (:));轨道。StateCovariance([1 3 5 2 4 6]， [1 3 5 2 4 6]);D = (estimate - true)' / cov * (estimate - true);结束
参考
bbb Sridharan, Ramaswamy和Antonio F. Pensa编。空间监视的视角．麻省理工学院出版社，2017。
