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

打开生活的脚本

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

空间碎片的场景

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

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

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

s =提高;rng (2020);现场= 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用四阶龙格-库塔数值积分将位置和速度在时间上进行传播。

首先，我们为空间碎片物体创建初始位置和速度。这是通过从随机分布中获得这些物体的传统轨道元素(半长轴、偏心度、倾角、升交点的经度、圆角和真异常角)来完成的。然后利用支撑函数将这些轨道元素转换为位置和速度矢量金宝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 >数据(我)。InitialVelocity = v;% #好< SAGROW >结束%创建平台并使用开普勒轨道运动模型分配它们的轨迹为i=1:numDebris碎片(i) =平台(场景);% #好< SAGROW >碎片(i)。轨迹= helperMotionTrajectory (@keplerorbit,…“SampleRate”, 0.1,…%积分步骤10秒“位置”、数据(我)。InitialPosition,…“速度”、数据(我).InitialVelocity);% #好< SAGROW >结束

空间监视雷达模型

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

在太平洋上建立一个太空监测站station1 =平台(场景,“位置”, (180 0));在大西洋建立第二个监测站station2 =平台(场景,“位置”-20年[0 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),…度“ReferenceRange”2 e6,…米“RangeLimits”[0 2 e6],…“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 = uifigure;观众= trackingGlobeViewer (f,“技术”，“卫星”，“ShowDroppedTracks”、假);%在全局查看器顶部添加信息框infotext = simulationInfoText (0, 0, 0);infobox = uilabel (f,“文本”infotext,“FontColor”(1 1 1),“字形大小”11…“位置”，[10 20 300 70]，“可见”，“上”）;显示雷达波束在地球上covcon = coverageConfig(现场);covcon plotCoverage(查看器,“ECEF”）;而推进(场景)信息框。文本= simulationInfoText(场景。SimulationTime 0 numDebris);plotPlatform(查看器碎片,“ECEF”）;结束%对场景进行快照快照(观众);

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

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

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

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

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

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

%定义跟踪器追踪= trackerJPDA (“FilterInitializationFcn”@initKeplerUKF,…“HasDetectableTrackIDsInput”,真的,…“ClutterDensity”1 e-20…“AssignmentThreshold”1 e3,…“DeletionThreshold”10 [7]);重置场景，种子和球体显示重启(现场);现场。StopTime = 1800;% 30分钟明确的(观众);减少碎片的历史深度查看器。PlatformHistoryDepth = 2;%显示10-sigma协方差省略号以提高能见度查看器。NumCovarianceSigma = 10;covcon plotCoverage(查看器,“ECEF”）;%初始化跟踪confTracks = objectTrack.empty (0,1);而推进(场景)plotPlatform(观众,碎片);时间= scene.SimulationTime;%生成检测检测=检测(现场);plotDetection(查看器,检测,“ECEF”）;%生成和更新轨道detectableInput = is可检测的(tracker,time, covcon);如果isLocked(tracker) || ~isempty(detections0) [confTracks， ~， allTracks,info] = tracker(detections0,time,detectableInput);confTracks = deleteBadTracks(跟踪、confTracks);confTracks plotTrack(查看器,“ECEF”）;结束信息框。Text = simulationInfoText(time, nummel (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.400 e7));Camorient(查看器，[360 -90 0]);drawnow情况下1500 campos(查看器，[54 2.3 6.09e6]);Camorient(查看器，[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;% 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);结束结束

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

函数过滤器= initKeplerUKF(检测)%假设雷达返回[x y z]测量模型= @(x,varargin) x([1 3 5]，:);detNoise = detection.MeasurementNoise;sigpos = 0.4;% msigvel = 0.4;% m / s ^ 2量= detection.Measurement;initState =[多边环境(1);0;量(2);0;量(3);0];过滤器= trackingUKF (@keplerorbit measmodel initState,…“StateCovariance”、诊断接头(repmat([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，《天体动力学基础》第2.5章μ= 398600.4405 * 1 e9;% s m ^ 3 ^ 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 (:);(阿兹,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 =元素个数(跟踪);删除= 0 (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]);如果Norm (pos) < 6500*1e3 || Norm (pos) > 8500*1e3 || max(cov，[]，“所有”) > 4*1e8 deleteTrack(tracker,tracks(i).TrackID);删除(i) = true;结束结束跟踪(删除)= [];结束

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