使用开普勒运动模型追踪太空碎片

打开直播脚本

这个例子展示了如何使用自定义运动模型来建立以地球为中心的轨迹trackingScenario，如何配置afusionRadarSensor在单基地模式下生成空间碎片的综合探测，以及如何设置多目标跟踪器来跟踪模拟目标。

空间碎片的场景

近地轨道上有3万多个大碎片(直径大于10厘米)，100多万个小碎片([1])。这些碎片可能会对人类在空间中的活动造成危险，损害操作卫星，并迫使时间敏感和昂贵的规避机动。随着空间活动的增加，减少和监测空间碎片变得至关重要。

您可以使用传感器融合和跟踪工具箱™来建模碎片轨迹，生成该碎片的合成雷达探测，并获得每个物体的位置和速度估计。

首先，创建跟踪方案并设置可重复结果的随机种子。

s =提高;rng (2020);现场= trackingScenario (“IsEarthCentered”,真的,'itsidadance'，'updateinterval'，．..“StopTime”, 3600,'updaterate'，0.1）;

您使用以地球为中心的地球固定（ECEF）参考帧。该框架的起源位于地球的中心，Z轴朝向北极点。X轴朝向赤道和格林威治经络的交点点。Y轴完成右手系统。使用该帧中的笛卡尔坐标定义了平台位置和速度。

定义碎屑运动模型

的helperMotionTrajectory此示例中使用的类定义了使用自定义运动模型函数的碎片对象轨迹。

围绕地球旋转的空间物体的轨迹可以用开普勒模型来近似，该模型假设地球是一个点质量物体，而围绕地球旋转的物体的质量可以忽略不计。地球重力场和环境扰动的高阶效应没有被考虑在内。由于运动方程是在非惯性参照系ECEF坐标系中表示的，因此考虑了科氏力和向心力。

ECEF碎片物体加速度矢量为

$\vec{一个} ＝ \frac{- μ}{r^{3.}} \vec{r} - 2 \vec{Ω} \times \frac{d}{d t} \vec{r} - \vec{Ω} \times （ \vec{Ω} \times \vec{r} ）$ ，

在哪里 $μ$ 是地球的标准引力参数， $\vec{r}$ 为ECEF碎片目标位置矢量， $r$ 是位置向量的范数吗 $\vec{Ω}$ 为地球自转矢量。

功能keplerorbit用四阶龙格-库塔数值积分将位置和速度在时间上进行传播。

首先，我们为空间碎片物体创建初始位置和速度。这是通过从随机分布中获得这些物体的传统轨道元素(半长轴、偏心度、倾角、升交点的经度、圆角和真异常角)来完成的。然后利用支撑函数将这些轨道元素转换为位置和速度矢量金宝appoe2rv．

%生成碎片的数量numDebris = 100;range = 7e6 + 1e5*randn(numDebris,1);ecc = 0.015 + 0.005*randn(numDebris,1);numDebris = 1 *rand(numDebris,1);局域网= 360 *兰德(numDebris, 1);w = 360 *兰德(numDebris, 1);ν= 360 *兰德(numDebris, 1);％转换为初始位置和速度为i = 1: numDebris [r, v] = oe2rv(范围(i), ecc (i),公司(我),局域网(i), w (i),ν(我));数据(我)。InitialPosition = r;% #好< SAGROW >数据（i）.initialvelocity = v;% #好< SAGROW >结束%创建平台并使用开普勒轨道运动模型分配它们的轨迹为i=1:numDebris碎片(i) =平台(场景);% #好< SAGROW >碎片(i)。轨迹= helperMotionTrajectory (@keplerorbit,．..“SampleRate”, 0.1,．..%积分步骤10秒“位置”，数据（i）.InitialPosition，．..“速度”、数据(我).InitialVelocity);% #好< SAGROW >结束

空间监视雷达模型

在该示例中，我们定义了四个反向站，其中扇形雷达梁调度为空间。粉丝通过碎片物体的轨道切割，以最大化物体检测的数量。一对车站位于太平洋和大西洋，而第二对监视站位于杆附近。具有四个分散的雷达允许重新检测空间碎片以校正其位置估计并获得新的碎片检测。

在太平洋上建立一个太空监测站station1 =平台(场景,“位置”, (180 0));在大西洋建立第二个监测站station2 =平台(场景,“位置”，[0 -20 0]）;靠近北极，在冰岛建立第三个监测站station3 =平台(场景,“位置”, (65 -20 0));在南极附近建立第四个监测站station4 =平台（场景，“位置”, -90 0 0);

每个观测站都配备了一个雷达模型fusionRadarSensor对象。为了探测近地轨道范围内的碎片目标，雷达有以下要求:

探测到2000公里外的10 dBsm物体
在2000公里范围内以100米的精度在水平和垂直方向分辨物体
在扇形视野中，在方位角120度和30度升高
根据其地理位置俯瞰空间

％创建扇形单体雷达，以监测空间碎片物体radar1 = fusionRadarSensor (1,．..'updaterate', 0.1,．..10秒“ScanMode”，'没有扫描'，．..“MountingAngles”, 90 0] [0,．..查找“FieldOfView”(120; 30),．..度'参考范围'2 e6,．..米“RangeLimits”，[02e6]，．..“ReferenceRCS”10，．..dBsm'hasfalsealarms'，错误的，．..“HasElevation”,真的,．..“AzimuthResolution”, 0.01,．..度“ElevationResolution”, 0.01,．..度“RangeResolution”, 100,．..米'hasins',真的,．..'检测'，“场景”）;station1。传感器= radar1;radar2 =克隆(radar1);radar2。SensorIndex = 2;station2。传感器= radar2;radar3 =克隆(radar1);radar3。SensorIndex = 3; station3.Sensors = radar3; radar4 = clone(radar1); radar4.SensorIndex = 4; station4.Sensors = radar4;

使用TrackingGlobeviewer可视化地面真相

你使用trackingGlobeViewer可视化跟踪场景中定义的所有元素:单个碎片对象及其轨迹、雷达扇形、雷达探测和轨迹。

f = uifigure;观众= trackingGlobeViewer (f,“技术”，“卫星”，“ShowDroppedTracks”，错误的）;%在全局查看器顶部添加信息框infotext = simulationInfoText (0, 0, 0);infobox = uilabel (f,“文本”，信息文本，'字体颜色'(1 1 1),“字形大小”11．..“位置”，[10 20 300 70]，“可见”，“上”）;显示雷达波束在地球上covcon = coverageConfig(现场);covcon plotCoverage(查看器,'ecref'）;尽管推进(场景)信息框。文本= simulationInfoText(场景。SimulationTime 0 numDebris);plotPlatform(查看器碎片,'ecref'）;结束%对场景进行快照快照(观众);

在虚拟地球上，你可以看到由白点代表的太空碎片，白色线条表示单独的跟踪轨迹。大多数产生的碎片物体都在倾角接近80度的轨道上。

轨迹在ECEF坐标下绘制，因此整个轨迹由于地球自转而向西旋转。经过几个轨道周期后，所有的空间碎片都会通过雷达的监视波束。

模拟合成探测和跟踪空间碎片

传感器模型使用地面真相来产生合成检测。打电话给检测方法在跟踪场景中获取场景中的所有检测。

多目标跟踪器trackerJPDA用于创建新的轨道、将检测与现有轨道关联、估计它们的状态并删除发散的轨道。设置属性hasdetableTrackidsInput.若为真，则允许跟踪器接受指示监视区域中跟踪对象是否可检测的输入。这对于不惩罚传播到雷达监视区域之外的轨道是很重要的。效用函数isDetectable计算每个模拟步骤中可检测到哪些轨道。

另外，还有一个实用函数deleteBadTracks用于更快地删除不同的曲目。

%定义跟踪器追踪= trackerJPDA (“FilterInitializationFcn”@initKeplerUKF,．..“HasDetectableTrackIDsInput”,真的,．..“ClutterDensity”，1E-20，．..“AssignmentThreshold”1 e3,．..“DeletionThreshold”，[7 10]）;重置场景，种子和球体显示重启(现场);现场。StopTime = 1800;% 30分钟清除（观众）;％减少碎片的历史深度查看器。PlatformHistoryDepth = 2;%显示10-sigma协方差省略号以提高能见度查看器。NumCovarianceSigma = 10;covcon plotCoverage(查看器,'ecref'）;％初始化曲目confTracks = objectTrack.empty (0,1);尽管前进（场景）PlotPlatform（查看者，碎片）;Time = Scene.simulationtime;%生成检测检测=检测(现场);plotDetection(查看器,检测,'ecref'）;％生成和更新曲目detectableInput = is可检测的(tracker,time, covcon);如果isLocked(tracker) || ~isempty(detections0) [confTracks， ~， allTracks,info] = tracker(detections0,time,detectableInput);confTracks = deleteBadTracks(跟踪、confTracks);confTracks plotTrack(查看器,'ecref'）;结束infobox.text = SimulationInfotext（时间，秩序（Conftracks），Numdebris）;%在模拟过程中移动相机并拍摄快照转变时间情况下100 campos(查看器，[90 150 5e6]);camorient(观众,[0 -65 345]);Drawnow im1 =快照(查看器);情况下270 campos(查看器，[60 -120 2.6e6]);Camorient(观众，[20 -45 20]);drawnow情况下380 campos(查看器，[60 -120 2.6e6]);Camorient(观众，[20 -45 20]);Drawnow im2 =快照(查看器);情况下400%重置坎波斯（观众，[17.3 -67.2 2.400E7]）;凸轮（查看器，[360 -90 0]）;drawnow情况下1500坎波（观众，[54 2.3 6.09E6]）;凸轮（查看器，[0 -73 348]）;drawnow情况下1560 im3 =快照(查看器);drawnow结束结束%恢复随机种子rng(年代);

imshow (im1);

在第一个快照中，您可以看到跟踪的对象为track T1(黄色)。该物体仅被探测两次，不足以降低轨迹的不确定性。因此，其协方差椭圆的大小比较大。您还可以观察另一个轨道T2蓝色，它被传感器检测了几次。结果表明，由于采用了更多的检测来校正状态估计，其对应的协方差椭圆要小得多。

imshow (im2);

几分钟后，如上所述，T1被删除，因为轨道的不确定性已经在没有检测的情况下变得太大。另一方面，由于额外的检测，第二轨道T2存活。

imshow（im3）

在上面的截图中，你可以看到跟踪T15(浅蓝色)即将进入雷达监视区域。Track T11(橙色)刚刚更新了一个检测，降低了其估计位置和速度的不确定性。在雷达站的配置下，经过30分钟的监视，在100个碎片中有18个轨道被初始化并确认。如果你增加模拟时间，雷达将覆盖太空360度，最终可以跟踪更多的碎片。可以探索不同的雷达站位置和配置，以增加跟踪目标的数量。

总结

在此示例中，您已经学习了如何指定自己的运动模型，以在跟踪方案中移动平台以及如何使用它们来设置跟踪器。这使您可以将该工具箱中提供的传感器融合和跟踪技术应用于更广泛的应用程序，例如在该示例中所示的地球固定坐标框架中建模和跟踪空间碎片的问题。

金宝app支持功能

本例中使用的运动模型如下所示。状态是物体的ECEF位置和速度[x;vx;y;v;z;款)．

功能国家= keplerorbit(状态,dt)%开普勒轨道执行数值积分来预测的状态%开普勒的身体。状态为[x;vx;y;vy;z;vz]%龙格-库塔4积分法: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;％m ^ 3 s ^ -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);结束结束

initkeplerukf.用你自己的运动模型初始化一个跟踪过滤器。在这个例子中，我们使用了建立地面真实的相同的运动模型，keplerorbit。

功能过滤器= initKeplerUKF(检测)%假设雷达返回[x y z]measmodel = @（x，varargin）x（[1 3 5]，:);Detnoise =检测。sigpos = 0.4;% msigvel = 0.4;% m / s ^ 2量= detection.Measurement;initState =[多边环境(1);0;量(2);0;量(3);0];过滤器= trackingUKF (@keplerorbit measmodel initState,．..'statecovariance'，诊断（Repmat（[10,10000]。^ 2,1,3）），．..“ProcessNoise”、诊断接头(repmat ([sigpos, sigvel]。^ 2,1,3)),．..'MeasurementNoise'，detnoise）;结束

oe2rv将一组6个传统的轨道元素转换为位置和速度向量。

功能(r, v) = oe2rv (a, e, i,局域网,w,ν)％参考：Bate，Mueller＆White，Astrovynamics Sec 2.5的基本原理μ= 398600.4405 * 1 e9;％m ^ 3 s ^ -2在周焦系统中表示r和vcnu = cosd(ν);snu =信德(ν);P = a*(1 - e^2);R = p/(1 + e*cnu);R_peri = [r*cnu;r * snu;0);V_peri = sqrt(mu/p)*[-snu;E + cnu;0);%转换成地心赤道坐标系家族= cosd(局域网);超人=信德(局域网);连续波= cosd (w);sw =信德(w);ci = cosd(我);如果=信德(我);R = [clan*cw-slan*sw*ci， -clan*sw-slan*cw*ci, slan*si;．..-slan* cw+clan*sw*ci， -clan*si;．..Sw *si, cw*si, ci];r = r * r_peri;v = R * v_peri;结束

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（:);[AZ，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））';结束结束

deleteBadTracks用来去除明显发散的轨迹。具体来说，本例中的发散轨迹是指当前位置落在地表上的轨迹，以及协方差过大的轨迹。

功能曲目= DELETEBADTRACKS（跟踪器，曲目）％删除不同的曲目：％ - 具有协方差> 4 * 1E8（20公里标准差）% -在LEO边界外的估计位置跟踪n = numel（轨道）;todelete = zeros（1，n，“逻辑”）;为i=1:numel(tracks) [pos, cov] = getTrackPositions(tracks(i)，[1 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 1 0]);如果规范（POS）<6500 * 1E3 ||常态（POS）> 8500 * 1E3 ||Max（COV，[]，'全部'）> 4 * 1E8 Deletetrack（跟踪器，曲目（i）.trackid）;todelete（i）= true;结束结束跟踪(删除)= [];结束

simulationInfoText用于显示模拟时间、碎片数量和场景中当前轨道数量。

功能text = simulationinfotext（时间，numtracks，numdebris）text = vertcat（String（[“经过的模拟时间:”num2str(圆(时间/ 60))“(min)”]),．..字符串([“当前曲目数量:”num2str (numTracks)]),．..字符串([“碎片总数:”num2str (numDebris))));drawnowlimitrate结束

参考文献

［1］https://www.esa.int/safety_security/space_debris/space_debris_by_ne_numbers.