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

打开实时脚本

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

空间碎片场景

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

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

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

S = rng;rng (2020);scene = trackingScenario()“IsEarthCentered”,真正的);

您使用地球中心-地球固定(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);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 >结束创建平台并使用开普勒轨道运动模型分配它们的轨迹为i=1:numDebris碎片(i) =平台(场景);% #好< SAGROW >碎片(i)。轨迹= helperMotionTrajectory(@keplerorbit，…“SampleRate”, 0.1,…%集成步骤10秒“位置”、数据(我)。InitialPosition,…“速度”、数据(我).InitialVelocity);% #好< SAGROW >结束

空间监视雷达模型

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

在太平洋建立一个太空监控站车站1 =站台(场景、“位置”，[10 180 0]);%在大西洋建立第二个监测站车站2 =站台(场景、“位置”，[0 -20 0]);在北极附近，在冰岛建立第三个监测站车站3 =站台(场景、“位置”，[65 -20]);在南极附近建立第四个监测站车站4 =站台(场景、“位置”，[-90 0 0]);

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

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

%制造扇形单站雷达以监测空间碎片物体radar1 = monostaticRadarSensor(1，…“UpdateRate”, 0.1,…10秒“ScanMode”，“没有扫描”，…“MountingAngles”，[0 90 0]，…查找“FieldOfView”(120; 30),…度“ReferenceRange”, 2000000,…米“ReferenceRCS”10…dBsm“HasFalseAlarms”假的,…“HasElevation”,真的,…“AzimuthResolution”, 0.01,…度“ElevationResolution”, 0.01,…度“RangeResolution”, 100,…米“之内”,真的,…“DetectionCoordinates”，“场景”);station1。传感器= radar1;Radar2 = clone(radar1);radar2。SensorIndex = 2;station2。传感器= radar2;Radar3 = clone(雷达1);radar3。SensorIndex = 3; station3.Sensors = radar3; radar4 = clone(radar1); radar4.SensorIndex = 4; station4.Sensors = radar4;

在虚拟的地球仪上想象地面的真相

在这个例子中，helperScenarioGlobeViewer提供一个虚拟地球仪，用于可视化跟踪场景中定义的所有元素:单个碎片物体及其轨迹、雷达风扇、雷达探测和轨迹。

globeDisplay = helperscenariglobeviewer;显示地球上的雷达波束covcon = coverageConfig(场景);plotCoverage (globeDisplay covcon);设置TargetHistoryLength以显示碎片物体的完整轨迹globeDisplay。targethhistorylength = 1000;现场。停止时间= 3600;现场。UpdateRate = 0.1;而advance(scene) time = scene. simulationtime;updateDisplay (globeDisplay、时间碎片);结束提前(globeDisplay);

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

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

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

传感器模型使用地面真值生成合成检测。调用检测方法上的跟踪场景，以获得所有的检测场景。

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

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

%定义跟踪器tracker = trackerJPDA(“FilterInitializationFcn”@initKeplerUKF,…“HasDetectableTrackIDsInput”,真的,…“ClutterDensity”1 e-20…“AssignmentThreshold”1 e4,…“DeletionThreshold”10 [7]);%复位场景，种子，和全球显示重启(现场);现场。StopTime = 1800;% 30分钟明确(globeDisplay);globeDisplay。targethhistorylength = 2;plotCoverage (globeDisplay covcon);%初始化音轨confTracks = objectTrack.empty(0,1);而advance(scene) time = scene. simulationtime;%生成检测Detections =检测(场景);%生成和更新音轨detectableInput = isDetectable(跟踪器，时间，covcon);如果~(isempty(detections) && ~isLocked(tracker)) [confTracks， ~， allTracks,info] = tracker(detections,time,detectableInput);confTracks = deleteBadTracks(跟踪器，confTracks);结束%更新全局显示updateDisplay (globeDisplay、时间碎片,检测,[],confTracks);在模拟过程中移动相机并拍摄快照开关时间情况下100 setCamera(globeDisplay，[90 150 5e6]，[0 -65 345]);im1 = snap(globeDisplay);情况下270 setCamera(globeDisplay，[60 -120 2.6e6]，[20 -45 20]);情况下380 setCamera(globeDisplay，[60 -120 2.6e6]，[20 -45 20]);im2 = snap(globeDisplay);情况下400%重置setCamera(globeDisplay，[17.3 -67.2 2.400e7]， [360 -90]);情况下setCamera(globeDisplay，[54 2.3 6.09e6]， [0 -73 348]);情况下1560 im3 = snap(globeDisplay);结束结束%恢复随机种子rng(年代);

imshow (im1);

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

imshow (im2);

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

imshow (im3)

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

总结

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

金宝app支持功能

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

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

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

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

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

函数[r,v] = 1 (a,e,i,l,w,nu)参考:贝特，穆勒和怀特，天体动力学基础第2.5节Mu = 398600.4405*1e9;% m^3 s^-2%表示局部系统中的r和vCnu = cosd(nu);Snu = sind(nu);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);转换成地心赤道坐标系Clan = cosd(lan);Slan = sind(lan);Cw = cosd(w);Sw = sind(w);Ci = cosd(i);Si = sind(i);R = [clan*sw-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 = 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))';结束结束

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

函数铁轨= deleteBadTracks(铁轨，铁轨)%移除发散轨迹:% -协方差bbbb4 *1e8的航迹(20公里标准差)% -在LEO边界之外估计位置的轨道N = numel(tracks);toDelete = 0 (1,n，“逻辑”);为i=1:numel(tracks) [pos, cov] = getTrackPositions(tracks(i)，[1 0 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;结束结束tracks(toDelete) = [];结束

参考文献

[1]https://www.esa.int/Safety_Security/Space_Debris/Space_debris_by_the_numbers