このページは機械翻訳を使用して翻訳されました。最新版の英語を参照するには、ここをクリックします。
地上レーダーによる LEO コンスタレーションの検出と追跡
この例では、衛星コンスタレーションの 2 行要素 (TLE) ファイルをインポートし、コンスタレーションのレーダー検出をシミュレートし、コンスタレーションを追跡する方法を示します。
地球を周回する宇宙物体のカタログを作成し、維持する作業は、宇宙監視において極めて重要です。このタスクは、新しいオブジェクトの検出と識別、カタログへの追加、カタログ内の既知のオブジェクトの軌道の更新、オブジェクトの寿命全体にわたる軌道の変化の追跡、大気圏への再突入の予測など、いくつかのプロセスで構成されます。この例では、新しい衛星を検出して追跡し、カタログに追加する方法を学びます。
宇宙での安全な運用を保証し、他の衛星や既知のデブリとの衝突を防ぐためには、新しく打ち上げられた衛星を正しく検出し追跡することが重要です。宇宙機関は通常、打ち上げ前の情報を共有しており、それを使用して検索戦略を選択できます。フェンス型レーダー システムで構成される低軌道 (LEO) 衛星検索戦略が一般的に使用されています。フェンス型レーダーシステムは、宇宙の有限の空間を探索し、視野を通過する衛星を検出します。この戦略により、新たに打ち上げられたコンスタレーションを迅速に検出し追跡することができる。[1]
TLE ファイルから衛星コンスタレーションをインポートする
2 行要素セットは、衛星の軌道情報を保存するための一般的なデータ形式です。satelliteScenario オブジェクトを使用して、TLE ファイルで定義された衛星軌道をインポートできます。デフォルトでは、インポートされた衛星軌道は、LEO オブジェクトに対して優れた精度を提供する SGP4 軌道伝播アルゴリズムを使用して伝播されます。この例では、これらの軌道は、新しく打ち上げられた衛星を検出するレーダー追跡システムの機能をテストするためのグラウンドトゥルースを提供します。
% Create a satellite scenario satscene = satelliteScenario; % Add satellites from TLE file. tleFile = "leoSatelliteConstellation.tle"; constellation = satellite(satscene, tleFile);
衛星シナリオ ビューアーを使用して、コンスタレーションを視覚化します。
play(satscene);
合成検出とトラックコンスタレーションのシミュレーション
宇宙監視レーダーのモデリング
宇宙を見つめる扇形のレーダービームを持つ 2 つのステーションを定義します。ファンは衛星の軌道を横切って検出数を最大化します。北米に位置するレーダー基地は東西のフェンスを形成します。
% First station coordinates in LLA station1 = [48 -80 0]; % Second station coordinates in LLA station2 = [50 -117 0];
各ステーションにはレーダーが装備されており、fusionRadarSensor オブジェクトを使用してモデル化されています。LEO 範囲の衛星を検出するために、レーダーには次の要件があります。
最大2000km離れた10dBsmの物体を検出
2000 km の範囲で 100 m の精度で水平方向と垂直方向の物体を解像します。
方位角120度、仰角30度の視野を持つ
宇宙を見上げる
% Create fan-shaped monostatic radars fov = [120;40]; radar1 = fusionRadarSensor(1,... 'UpdateRate',0.1,... 10 sec 'ScanMode','No scanning',... 'MountingAngles',[0 90 0],... look up 'FieldOfView',fov,... degrees 'ReferenceRange',2000e3,... m 'RangeLimits', [0 2000e3], ... m 'ReferenceRCS', 10,... dBsm 'HasFalseAlarms',false,... 'HasNoise', true,... 'HasElevation',true,... 'AzimuthResolution',0.03,... degrees 'ElevationResolution',0.03,... degrees 'RangeResolution',2000, ... m % accuracy ~= 2000 * 0.05 (m) 'DetectionCoordinates','Sensor Spherical',... 'TargetReportFormat','Detections'); radar2 = clone(radar1); radar2.SensorIndex = 2;
レーダー処理チェーン
この例では、レーダー追跡チェーンを適切に実行するために、いくつかの座標変換と軸変換が実行されます。以下の図は、上で定義した入力がどのように変換され、レーダーに渡されるかを示しています。
最初のステップでは、ローカル レーダー ステーションの NED 軸で各衛星の姿勢を計算します。これを実現するには、まず地上局の ECEF 姿勢を取得し、衛星の位置と速度を ECEF 座標に変換します。レーダー入力は、衛星の姿勢と地上局の姿勢の差を取得し、その差を地上局のローカル NED 軸に回転させることによって取得されます。実装の詳細については、assembleRadarInputs サポート関数を参照してください。
2 番目のステップでは、トラッカーが ECEF 状態で動作できるように、検出オブジェクトに必要な情報を追加します。この目的を達成するには、addMeasurementParams サポート関数に示されているように、各オブジェクト検出で MeasurementParameters プロパティを使用します。
トラッカーの定義
上記で定義したレーダー モデルは検出結果を出力します。衛星の軌道を推定するには、トラッカーを使用します。Sensor Fusion and Tracking Toolbox ™ は、さまざまなマルチオブジェクト トラッカーを提供します。この例では、追跡パフォーマンスと計算コストのバランスが優れているため、Joint Probabilistic Data Association (JPDA) トラッカーを選択します。
トラッカーの追跡フィルターを定義する必要があります。衛星を追跡するには、運動方程式のケプラー積分など、SGP4 よりも忠実度の低いモデルを使用できます。多くの場合、ターゲットのモーション モデルの忠実度の欠如は、測定値の更新とフィルターへのプロセス ノイズの組み込みによって補われます。サポート関数 initKeplerUKF は追跡フィルターを定義します。
% Define the tracker tracker = trackerJPDA('FilterInitializationFcn',@initKeplerUKF,... 'HasDetectableTrackIDsInput',true,... 'ClutterDensity',1e-40,... 'AssignmentThreshold',1e4,... 'DeletionThreshold',[5 8],... 'ConfirmationThreshold',[5 8]);
シミュレーションの実行
この例の残りの部分では、シナリオをステップ実行して、レーダー検出をシミュレートし、衛星を追跡します。このセクションでは視覚化のために trackingGlobeViewer を使用します。このクラスを使用して、不確実性の楕円とともにセンサーと追跡データを表示し、各衛星の実際の位置を示します。
viewer = trackingGlobeViewer('ShowDroppedTracks',false, 'PlatformHistoryDepth',700); % Define coverage configuration of each radar and visualize it on the globe ned1 = dcmecef2ned(station1(1), station1(2)); ned2 = dcmecef2ned(station2(1), station2(2)); covcon(1) = coverageConfig(radar1,lla2ecef(station1),quaternion(ned1,'rotmat','frame')); covcon(2) = coverageConfig(radar2,lla2ecef(station2),quaternion(ned2, 'rotmat','frame')); plotCoverage(viewer, covcon, 'ECEF');
まず、5 時間にわたるコンスタレーションの状態の履歴全体を生成します。次に、レーダー検出をシミュレートし、ループでトラックを生成します。
satscene.StopTime = satscene.StartTime + hours(5); satscene.SampleTime = 10; numSteps = ceil(seconds(satscene.StopTime - satscene.StartTime)/satscene.SampleTime); % Get constellation positions and velocity over the course of the simulation plats = repmat(... struct('PlatformID',0,'Position',[0 0 0], 'Velocity', [0 0 0]),... numSteps, 40); for i=1:numel(constellation) [pos, vel] = states(constellation(i),'CoordinateFrame',"ECEF"); for j=1:numSteps plats(j,i).Position = pos(:,j)'; plats(j,i).Velocity = vel(:,j)'; plats(j,i).PlatformID = i; end end % Initialize tracks and tracks log confTracks = objectTrack.empty(0,1); trackLog = cell(1,numSteps); % Initialize radar plots radarplt = helperRadarPlot(fov); % Set random seed for reproducible results s = rng; rng(2020); step = 0; while step < numSteps time = step*satscene.SampleTime; step = step + 1; % Generate detections of the constellation following the radar % processing chain 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); detections = [dets1; dets2]; updateRadarPlots(radarplt,targets1, targets2 ,dets1, dets2) % Generate and update tracks detectableInput = isDetectable(tracker,time, covcon); if ~isempty(detections) || isLocked(tracker) [confTracks,~,~,info] = tracker(detections,time,detectableInput); end trackLog{step} = confTracks; % Update viewer plotPlatform(viewer, plats(step,:),'ECEF', 'Color',[1 0 0],'LineWidth',1); plotDetection(viewer, detections,'ECEF'); plotTrack(viewer, confTracks, 'ECEF','Color',[0 1 0], 'LineWidth',3); end
![{"String":"Figure contains 2 axes objects. Axes object 1 with title Radar 1 Field of View contains 772 objects of type line. Axes object 2 with title Radar 2 Field of View contains 1052 objects of type line.","Tex":["Radar 1 Field of View","Radar 2 Field of View"],"LaTex":[]}](detectandtrackaleosatellithwithgroundradarsexample_02_ja_JP.png)
上のグラフは、各レーダーの観点から見た軌道 (青い点) と検出 (赤い円) を示しています。
% Restore previous random seed state
rng(s);figure; snapshot(viewer);
5時間の追跡後、コンスタレーションの約半分の追跡に成功しました。この構成では衛星が長期間検出されないことが多いため、部分的な軌道カバレッジで追跡を維持するのは困難です。この例では、レーダーステーションは 2 つだけです。追跡ステーションの追加により、追跡パフォーマンスが向上することが期待されます。実際のオブジェクトとトラックを比較して追跡パフォーマンスを評価する割り当てメトリックを以下に示します。
% Show Assignment metrics truthIdFcn = @(x)[x.PlatformID]; tam = trackAssignmentMetrics('DistanceFunctionFormat','custom',... 'AssignmentDistanceFcn',@distanceFcn,... 'DivergenceDistanceFcn',@distanceFcn,... 'TruthIdentifierFcn',truthIdFcn,... 'AssignmentThreshold',1000,... 'DivergenceThreshold',2000); for i=1:numSteps % Extract the tracker and ground truth at the i-th tracker update tracks = trackLog{i}; truths = plats(i,:); % Extract summary of assignment metrics against tracks and truths [trackAM,truthAM] = tam(tracks, truths); end
% Show cumulative metrics for each individual recorded truth object results = truthMetricsTable(tam); results(:,{'TruthID','AssociatedTrackID','BreakLength','EstablishmentLength'})
ans=40×4 table
TruthID AssociatedTrackID BreakLength EstablishmentLength
_______ _________________ ___________ ___________________
1 52 5 232
2 24 0 594
3 4 0 56
4 55 11 492
5 18 0 437
6 48 8 811
7 21 0 513
8 27 0 661
9 39 0 1221
10 50 0 1504
11 43 0 1339
12 37 0 1056
13 NaN 0 1800
14 NaN 0 1800
15 NaN 0 1800
16 NaN 0 1800
⋮
上記の表には、打ち上げられたコンスタレーションの 40 個の衛星がリストされており、追跡された衛星と関連する追跡 ID が表示されています。トラック ID の値が NaN の場合、シミュレーションの終了までに衛星が追跡されていないことを示します。これは、衛星の軌道が 2 つのレーダーのいずれかの視野を通過しなかったか、衛星の軌道が外れたことを意味します。トラッカーは、初期検出数が不十分なためにトラックをドロップする可能性があり、その結果、推定に大きな不確実性が生じます。あるいは、衛星がすぐに再検出されない場合、トラッカーはトラックをドロップする可能性があり、更新が不足すると、相違が生じ、最終的には削除されます。
まとめ
この例では、Aerospace Toolbox ™ の satelliteScenario オブジェクトを使用して、TLE ファイルから軌道情報をインポートする方法を学習しました。SGP4 を使用して衛星の軌道を伝播し、衛星シナリオ ビューアーを使用してシナリオを視覚化しました。Sensor Fusion and Tracking Toolbox ™ のレーダーおよびトラッカー モデルを使用して、宇宙監視レーダー追跡システムをモデル化する方法を学びました。構築された追跡システムは、低忠実度モデルを使用して各衛星の推定軌道を予測できます。
サポート機能
initKeplerUKF は、ケプラー運動モデルを使用して、Unscented Kalman フィルターを初期化します。長い予測期間にわたって無香料フィルタの堅牢性を確保するには、Alpha = 1、Beta = 0、Kappa = 0 に設定します。
function filter = initKeplerUKF(detection) radarsphmeas = detection.Measurement; [x, y, z] = sph2cart(deg2rad(radarsphmeas(1)),deg2rad(radarsphmeas(2)),radarsphmeas(3)); radarcartmeas = [x y z]; Recef2radar = detection.MeasurementParameters.Orientation; ecefmeas = detection.MeasurementParameters.OriginPosition + radarcartmeas*Recef2radar; % This is equivalent to: % Ry90 = [0 0 -1 ; 0 1 0; 1 0 0]; % frame rotation of 90 deg around y axis % nedmeas(:) = Ry90' * radarcartmeas(:); % ecefmeas = lla2ecef(lla) + nedmeas * dcmecef2ned(lla(1),lla(2)); initState = [ecefmeas(1); 0; ecefmeas(2); 0; ecefmeas(3); 0]; sigpos = 2;% m sigvel = 0.5;% m/s^2 filter = trackingUKF(@keplerorbit,@cvmeas,initState,... 'Alpha', 1, 'Beta', 0, 'Kappa', 0, ... 'StateCovariance', diag(repmat([1000, 10000].^2,1,3)),... 'ProcessNoise',diag(repmat([sigpos, sigvel].^2,1,3))); end function state = keplerorbit(state,dt) % keplerorbit performs numerical integration to predict the state of % Keplerian bodies. The state is [x;vx;y;vy;z;vz] % Runge-Kutta 4th order integration method: k1 = kepler(state); 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; function dstate=kepler(state) x =state(1,:); vx = state(2,:); y=state(3,:); vy = state(4,:); z=state(5,:); vz = state(6,:); mu = 398600.4405*1e9; % m^3 s^-2 omega = 7.292115e-5; % rad/s r = norm([x y z]); g = mu/r^2; % Coordinates are in a non-intertial frame, account for Coriolis % and centripetal acceleration ax = -g*x/r + 2*omega*vy + omega^2*x; ay = -g*y/r - 2*omega*vx + omega^2*y; az = -g*z/r; dstate = [vx;ax;vy;ay;vz;az]; end end
この例では、isDetectable を使用して、特定の時間にどのトラックが検出可能かを判断します。
function detectInput = isDetectable(tracker,time,covcon) if ~isLocked(tracker) detectInput = zeros(0,1,'uint32'); return end tracks = tracker.predictTracksToTime('all',time); if isempty(tracks) detectInput = zeros(0,1,'uint32'); else alltrackid = [tracks.TrackID]; isDetectable = zeros(numel(tracks),numel(covcon),'logical'); for i = 1:numel(tracks) track = tracks(i); pos_scene = track.State([1 3 5]); for j=1:numel(covcon) config = covcon(j); % rotate position to sensor frame: d_scene = pos_scene(:) - config.Position(:); scene2sens = rotmat(config.Orientation,'frame'); d_sens = scene2sens*d_scene(:); [az,el] = cart2sph(d_sens(1),d_sens(2),d_sens(3)); if abs(rad2deg(az)) <= config.FieldOfView(1)/2 && abs(rad2deg(el)) < config.FieldOfView(2)/2 isDetectable(i,j) = true; else isDetectable(i,j) = false; end end end detectInput = alltrackid(any(isDetectable,2))'; end end
assembleRadarInput は、各レーダー本体フレームのコンスタレーションポーズを構築するために使用されます。これは図で説明されている最初のステップです。
function targets = assembleRadarInputs(station, constellation) % For each satellite in the constellation, derive its pose with respect to % the radar frame. % inputs: % - station : LLA vector of the radar ground station % - constellation : Array of structs containing the ECEF position % and ECEF velocity of each satellite % outputs: % - targets : Array of structs containing the pose of each % satellite with respect to the radar, expressed in the local % ground radar frame (NED) % Template structure for the outputs which contains all the field required % by the radar step method targetTemplate = struct( ... 'PlatformID', 0, ... 'ClassID', 0, ... 'Position', zeros(1,3), ... 'Velocity', zeros(1,3), ... 'Acceleration', zeros(1,3), ... 'Orientation', quaternion(1,0,0,0), ... 'AngularVelocity', zeros(1,3),... 'Dimensions', struct( ... 'Length', 0, ... 'Width', 0, ... 'Height', 0, ... 'OriginOffset', [0 0 0]), ... 'Signatures', {{rcsSignature}}); % First fill in the current satellite ECEF pose targetPoses = repmat(targetTemplate, 1, numel(constellation)); for i=1:numel(constellation) targetPoses(i).Position = constellation(i).Position; targetPoses(i).Velocity = constellation(i).Velocity; targetPoses(i).PlatformID = constellation(i).PlatformID; % Orientation and angular velocity are left null, assuming satellite to % be point targets with a uniform rcs end % Then derive the radar pose in ECEF based on the ground station location Recef2station = dcmecef2ned(station(1), station(2)); radarPose.Orientation = quaternion(Recef2station,'rotmat','frame'); radarPose.Position = lla2ecef(station); radarPose.Velocity = zeros(1,3); radarPose.AngularVelocity = zeros(1,3); % Finally, take the difference and rotate each vector to the ground station % NED axes targets = targetPoses; for i=1: numel(targetPoses) thisTgt = targetPoses(i); pos = Recef2station*(thisTgt.Position(:) - radarPose.Position(:)); vel = Recef2station*(thisTgt.Velocity(:) - radarPose.Velocity(:)) - cross(radarPose.AngularVelocity(:),pos(:)); angVel = thisTgt.AngularVelocity(:) - radarPose.AngularVelocity(:); orient = radarPose.Orientation' * thisTgt.Orientation; % Store into target structure array targets(i).Position(:) = pos; targets(i).Velocity(:) = vel; targets(i).AngularVelocity(:) = angVel; targets(i).Orientation = orient; end end
addMeasurementParam は、図に示されているように、レーダー チェーン プロセスのステップ 2 を実装します。
function dets = addMeasurementParams(dets, numdets, station) % Add radar station information to the measurement parameters Recef2station = dcmecef2ned(station(1), station(2)); for i=1:numdets dets{i}.MeasurementParameters.OriginPosition = lla2ecef(station); dets{i}.MeasurementParameters.IsParentToChild = true; % parent = ecef, child = radar dets{i}.MeasurementParameters.Orientation = dets{i}.MeasurementParameters.Orientation' * Recef2station; end end
distanceFcn は、追跡割り当てを評価するために割り当てメトリックとともに使用されます。
function d = distanceFcn(track, truth) estimate = track.State([1 3 5 2 4 6]); true = [truth.Position(:) ; truth.Velocity(:)]; cov = track.StateCovariance([1 3 5 2 4 6], [1 3 5 2 4 6]); d = (estimate - true)' / cov * (estimate - true); end
リファレンス
[1] Sridharan、Ramaswamy、Antonio F. Pensa 編。宇宙監視の展望。MITプレス、2017年。
