利用开普勒运动模型跟踪空间碎片

打开实时脚本

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

空间碎片场景

近地轨道(LEO)[1]中有3万多个大型碎片物体(直径大于10cm)和100多万个较小的碎片物体。这些碎片可能对人类在太空中的活动造成危险，破坏正在运行的卫星，并迫使时间敏感和代价高昂的规避操作。随着太空活动的增加，减少和监测太空碎片变得至关重要。

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

首先，创建一个跟踪场景，并为可重复的结果设置随机种子。

S = rng;rng (2020);scene = trackingScenario(“IsEarthCentered”,真的,“InitialAdvance”，“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下文使用该方程的四阶龙格-库塔数值积分来传播位置和速度。

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

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

空间监视雷达模型

在这个例子中，我们定义了四个用扇形雷达波束观测太空的对跖站。风扇切断了碎片物体的轨道，以最大限度地检测到物体的数量。一对监测站位于太平洋和大西洋，而另一对监测站位于两极附近。拥有四个分散的雷达可以重新探测空间碎片，以纠正其位置估计，并获得新的碎片探测。

在太平洋建立一个太空监测站Station1 =平台(场景，“位置”，[10 180 0]);在大西洋建立第二个监测站Station2 =平台(场景，“位置”，[0 -20 0]);在北极附近，在冰岛建立第三个监测站Station3 =平台(场景，“位置”，[65 -20 0]);在南极附近建立第四个监测站。Station4 =站台(场景，“位置”，[-90 0 0]);

每个站都配备了一个模拟雷达fusionRadarSensor对象。为了探测LEO范围内的碎片物体，雷达有以下要求:

探测2000公里外的10 dBsm物体
在2000公里范围内水平和垂直分辨目标，精度为100米
方位角120度，仰角30度的扇形视野
根据它的地理位置仰望太空

创建扇形单站雷达来监测空间碎片物体radar1 = fusionRadarSensor(1，.．.“UpdateRate”, 0.1,.．.10秒“ScanMode”，“没有扫描”，.．.“MountingAngles”，[0 90 0]，.．.查找“FieldOfView”(120; 30),.．.度“ReferenceRange”2 e6,.．.米“RangeLimits”， [0 2e6]，.．.“ReferenceRCS”10.．.dBsm“HasFalseAlarms”假的,.．.“HasElevation”,真的,.．.“AzimuthResolution”, 0.01,.．.度“ElevationResolution”, 0.01,.．.度“RangeResolution”, 100,.．.米“之内”,真的,.．.“DetectionCoordinates”，“场景”);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 = ufigure;查看器= trackingGlobeViewer(f，“技术”，“卫星”，“ShowDroppedTracks”、假);在全球查看器顶部添加信息框infotext = simulationInfoText(0,0,0);Infobox = uilabel(f，“文本”infotext,“FontColor”，[11 11 1]，“字形大小”11.．.“位置”，[10 20 300 70]，“可见”，“上”);显示地球上的雷达波束covcon = coverageConfig(场景);covcon plotCoverage(查看器,“ECEF”);而推进(场景)信息框。文本= simulationInfoText(场景。SimulationTime, 0, numDebris);plotPlatform(查看器碎片,“ECEF”);结束截取场景的快照快照(观众);

在虚拟地球仪上，你可以看到空间碎片用白点表示，单独的轨迹用白线表示。大多数产生的碎片物体都在接近80度的高倾角轨道上。

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

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

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

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

另外，还有一个效用函数deleteBadTracks用于更快地删除发散轨道。

定义跟踪器追踪器=追踪器“FilterInitializationFcn”@initKeplerUKF,.．.“HasDetectableTrackIDsInput”,真的,.．.“ClutterDensity”1 e-20.．.“AssignmentThreshold”1 e3,.．.“DeletionThreshold”10 [7]);重置场景、种子和全局显示重启(现场);现场。StopTime = 1800;% 30分钟明确的(观众);减少碎片的历史深度查看器。PlatformHistoryDepth = 2;%显示10西格玛协方差椭圆更好的可见性查看器。NumCovarianceSigma = 10;covcon plotCoverage(查看器,“ECEF”);初始化轨道confTracks = objectTrack.empty(0,1);而前进(场景)情节平台(观众、碎片);time = scene.SimulationTime;%生成检测检测=检测(场景);plotDetection(查看器,检测,“ECEF”);生成和更新曲目detectableInput = isdetected(跟踪器，时间，covcon);如果isLocked(跟踪)|| ~isempty(检测)[confTracks， ~， allTracks,info] =跟踪(检测，时间，detectableInput);confTracks = deleteBadTracks(跟踪器，confTracks);confTracks plotTrack(查看器,“ECEF”);结束信息框。文本= simulationInfoText(时间，数字(confTracks)， numDebris);在模拟过程中移动相机并拍摄快照开关时间情况下100 campos(viewer，[90 150 5e6]);Camorient (viewer，[0 -65 345]);绘制im1 =快照(查看器);情况下270 campos(viewer，[60 -120 2.6e6]);Camorient(观察者，[20 -45 20]);drawnow情况下380 campos(viewer，[60 -120 2.6e6]);Camorient(观察者，[20 -45 20]);绘制im2 =快照(查看器);情况下400%重置Campos (viewer，[17.3 -67.2 2.400e7]);Camorient (viewer， [360 -90 0]);drawnow情况下1500 campos(viewer，[54 2.3 6.09e6]);Camorient (viewer， [0 -73 348]);drawnow情况下1560 im3 =快照(查看器);drawnow结束结束%恢复随机种子rng(年代);

imshow (im1);

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

imshow (im2);

几分钟后，如上图所示，T1被删除了，因为轨道的不确定性已经变得太大了，没有检测到。另一方面，由于额外的检测，第二条轨道T2幸存下来。

imshow (im3)

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

总结

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

金宝app支持功能

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

函数状态= keplerorbit(状态，dt)% keplerorbit执行数值积分来预测的状态%开普勒体。状态是[x;vx;y;vy;z;vz]%龙格-库塔4积分法: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];结束结束

initKeplerUKF初始化跟踪过滤器与您自己的运动模型。在这个例子中，我们使用与建立基本真理相同的运动模型，keplerorbit。

函数filter = initKeplerUKF(检测)%假设雷达返回[x y z]Measmodel = @(x,varargin) x([1 3 5]，:);detNoise =检测。测量噪声;Sigpos = 0.4;% mSigvel = 0.4;% m / s ^ 2meas =检测。initState = [meas(1);0;量(2);0;量(3);0];filter = trackingUKF(@keplerorbit,measmodel,initState，.．.“StateCovariance”， diag(repmat([10,10000].^2,1,3))，.．.“ProcessNoise”、诊断接头(repmat ([sigpos, sigvel]。^ 2,1,3)),.．.“MeasurementNoise”, detNoise);结束

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

函数[r,v] = oe2rv(a,e,i,lan,w,nu)参考资料:Bate, Mueller & White，《天体动力学基础》第2.5章Mu = 398600.4405*1e9;% m^3 s^ 2%在焦周系统中表达r和vCnu = cosd(nu);Snu = sind(nu);P = a*(1 - e²);R = p/(1 + e*cnu);R_peri = [r*cnu;r * snu;0);V_peri =√(mu/p)*[-snu;E + cnu;0);%转换成地心赤道坐标系Clan = cosd(lan);Slan = sind(lan);Cw = cosd(w);Sw = sind(w);Ci = cosd(i);Si = sind(i);R = [clan*cw-slan*sw*ci， -clan*sw-slan*cw*ci, slan*si;.．.Slan *cw+clan*sw*ci， -slan*sw+clan*cw*ci， -clan*si;.．.Sw *si, cw*si, ci];r = r *r_peri;v = R*v_peri;结束

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))';结束结束

deleteBadTracks用于去除明显偏离的轨迹。具体来说，本例中的发散轨迹是当前位置落在地球表面的轨迹和协方差过大的轨迹。

函数tracks = deleteBadTracks(跟踪器，轨道)%去除发散轨迹:% -轨迹协方差> 4*1e8(20公里标准差)% -在近地轨道范围外估计位置的轨迹N =数字(轨道);toDelete = 0 (1,n，“逻辑”);为i=1: number (tracks) [pos, cov] = getTrackPositions(tracks(i)，[1 0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0]);如果Norm (pos) < 6500*1e3 || Norm (pos) > 8500*1e3 || max(cov，[]，“所有”) > 4*1e8 deleteTrack(跟踪器，跟踪(i).TrackID);删除(i) = true;结束结束tracks(toDelete) = [];结束

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

函数text = simulationInfoText(time,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_the_numbers