Этот пример показывает, как отслеживать маневрирующие цели с помощью различных отслеживающих фильтров. Пример показывает различие между фильтрами, которые используют одну модель движения и несколько моделей движения.
В этом примере вы задаете одну цель, которая первоначально перемещается с постоянной скоростью 200 м/с в течение 33 секунд, затем входит в постоянный поворот 10 град/с. Поворот длится 33 секунды, затем цель разгоняется по прямой со скоростью 3 м/с ^ 2.
[trueState, time, fig1] = helperGenerateTruthData; dt = diff(time(1:2)); numSteps = numel(time); figure(fig1)
Определите измерения, чтобы быть положением и добавьте нормальный случайный шум со стандартным отклонением 1 к измерениям.
% Set the RNG seed for repeatable results s = rng; rng(2018); positionSelector = [1 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 1 0]; % Position from state truePos = positionSelector * trueState; measNoise = randn(size(truePos)); measPos = truePos + measNoise;
Вы задаете trackingEKF
с моделью движения с постоянной скоростью. Вы используете первое измерение, чтобы задать начальное состояние и ковариацию состояния, и установите шум процесса как не аддитивный, чтобы задать шум процесса с точки зрения неизвестного ускорения в компонентах x, y и z. Это определение подобно тому, как функция initcvekf
работает.
initialState = positionSelector' * measPos(:,1); initialCovariance = diag([1,1e4,1,1e4,1,1e4]); % Velocity is not measured cvekf = trackingEKF(@constvel, @cvmeas, initialState, ... 'StateTransitionJacobianFcn', @constveljac, ... 'MeasurementJacobianFcn', @cvmeasjac, ... 'StateCovariance', initialCovariance, ... 'HasAdditiveProcessNoise', false, ... 'ProcessNoise', eye(3));
Для каждого измерения вы прогнозируете фильтр, вычисляете расстояние предсказанного состояния от истинного положения и корректируете фильтр с помощью измерения, чтобы получить отфильтрованную оценку положения.
dist = zeros(1,numSteps); estPos = zeros(3,numSteps); for i = 2:size(measPos,2) predict(cvekf, dt); dist(i) = distance(cvekf,truePos(:,i)); % Distance from true position estPos(:,i) = positionSelector * correct(cvekf, measPos(:,i)); end figure(fig1); plot(estPos(1,:),estPos(2,:),'.g','DisplayName','CV Low PN') title('True and Estimated Positions') axis([5000 8000 -500 2500])
Как показано на рисунке, фильтр способен очень хорошо отслеживать часть движения с постоянной скоростью, но когда цель выполняет поворот, оцененное положение фильтра расходится с истинным положением. Расстояние от оценки до истины можно увидеть на следующем графике. Во время поворота на 33-66 секундах нормированное расстояние скачет до очень высоких значений, а значит, фильтр не способен отслеживать маневрирующую цель.
fig2 = figure; hold on plot((1:numSteps)*dt,dist,'g','DisplayName', 'CV Low PN') title('Normalized Distance From Estimated Position to True Position') xlabel('Time (s)') ylabel('Normalized Distance') legend
Одним из возможных решений является увеличение шума процесса. Шум процесса представляет немоделированные условия в модели движения. Для модели постоянной скорости это неизвестные условия ускорения. Увеличивая шум процесса, вы допускаете большую неопределенность в модели движения, которая заставляет фильтр полагаться на измерения больше, чем на модель. Следующие линии создают фильтр постоянной скорости с высокими значениями шума процесса, которые соответствуют примерно 5-G повороту.
cvekf2 = trackingEKF(@constvel, @cvmeas, initialState, ... 'StateTransitionJacobianFcn', @constveljac, ... 'MeasurementJacobianFcn', @cvmeasjac, ... 'StateCovariance', initialCovariance, ... 'HasAdditiveProcessNoise', false, ... 'ProcessNoise', diag([50,50,1])); % Large uncertainty in the horizontal acceleration dist = zeros(1,numSteps); estPos = zeros(3,numSteps); for i = 2:size(measPos,2) predict(cvekf2, dt); dist(i) = distance(cvekf2,truePos(:,i)); % Distance from true position estPos(:,i) = positionSelector * correct(cvekf2, measPos(:,i)); end figure(fig1) plot(estPos(1,:),estPos(2,:),'.c','DisplayName','CV High PN') axis([5000 8000 -500 2500])
Увеличение шума процесса значительно улучшает способность фильтра отслеживать цель во время поворота. Однако есть стоимость: фильтр менее способен сглаживать измерение шум в течение периода постоянной скорости движения. Несмотря на то, что нормированное расстояние во время поворота было значительно уменьшено, нормированное расстояние увеличилось в первые 33 секунды, во время периода постоянной скорости движения.
figure(fig2) plot((1:numSteps)*dt,dist,'c','DisplayName', 'CV High PN') axis([0 100 0 50])
Другое решение - использовать фильтр, который может рассматривать все модели движения одновременно, называемый взаимодействующим фильтром с несколькими моделями (IMM). Фильтр IMM может поддерживать столько моделей движения, сколько вы хотите, но обычно используется с 2-5 моделями движения. Для этого примера достаточно трех моделей: модели постоянной скорости, модели постоянного поворота и модели постоянного ускорения.
imm = trackingIMM('TransitionProbabilities', 0.99); % The default IMM has all three models % Initialize the state and state covariance in terms of the first model initialize(imm, initialState, initialCovariance);
Фильтр IMM используется так же, как и EKF.
dist = zeros(1,numSteps); estPos = zeros(3,numSteps); modelProbs = zeros(3,numSteps); modelProbs(:,1) = imm.ModelProbabilities; for i = 2:size(measPos,2) predict(imm, dt); dist(i) = distance(imm,truePos(:,i)); % Distance from true position estPos(:,i) = positionSelector * correct(imm, measPos(:,i)); modelProbs(:,i) = imm.ModelProbabilities; end figure(fig1) plot(estPos(1,:),estPos(2,:),'.m','DisplayName','IMM')
The trackingIMM
фильтр способен отслеживать маневрирующую цель во всех трех частях движения.
Исследуя расстояние между предсказанным состоянием фильтра и истинным положением, вы видите, что фильтр IMM способен уменьшить расстояние для всех фрагментов движения. На самом деле, фильтр IMM лучше отслеживает движение, чем обе другие модели постоянной скорости, используемые ранее.
figure(fig2) hold on plot((1:numSteps)*dt,dist,'m','DisplayName', 'IMM') axis([0 100 0 50])
Чтобы лучше понять, как работает фильтр IMM, постройте график вероятностей модели как функции времени. Рисунок показывает, что фильтр инициализируется с тремя моделями, имеющими одинаковую вероятность. Когда фильтр обновляется, он быстро сходится к очень высокой вероятности, что модель является моделью постоянной скорости. После 33 секунд движения модель постоянной скорости больше не верна, и вероятность модели постоянного поворота становится очень высокой на время поворота. В последней части движения во время маневра постоянного ускорения фильтр IMM присваивает высокую вероятность того, что движение является постоянным ускорением, но фильтр менее уверен в правильности модели движения, и существует около 0,3 вероятности движения постоянной скорости.
figure plot((1:numSteps)*dt, modelProbs) title('Model Probabilities vs. Time') xlabel('Time (s)') ylabel('Model Probabilities') legend('IMM-CV','IMM-CA','IMM-CT') % Return the RNG to its previous state rng(s)
Этот пример показал вам, как отслеживать маневрирование цели с постоянным поворотом и постоянным ускорением движения. Пример показал, как можно увеличить шум процесса, чтобы захватить неизвестный маневр с помощью модели постоянной скорости. Вы также видели, как улучшить отслеживание маневрирующей цели с помощью IMM-фильтра.
helperGenerateTruthData
Эта функция генерирует основную истину траекторию.
function [Xgt, tt, figs] = helperGenerateTruthData % Generate ground truth vx = 200; % m/s omega = 10; % deg/s acc = 3; % m/s/s dt = 0.1; tt = (0:dt:floor(1000*dt)); figs = []; Xgt = NaN(9,numel(tt)); Xgt(:,1) = 0; % Constant velocity seg1 = floor(numel(tt)/3); Xgt(2,1) = vx; slct = eye(9); slct(3:3:end,:) = []; for m = 2:seg1 X0 = slct*Xgt(:,m-1); X1 = constvel(X0, dt); X1 = slct'*X1; Xgt(:,m) = X1; end % Constant turn seg2 = floor(2*numel(tt)/3); slct = eye(9); slct(3:3:end,:) = []; for m = seg1+1:seg2 X0 = slct*Xgt(:,m-1); X0 = [X0(1:4);omega]; X1 = constturn(X0, dt); X1 = X1(1:4); X1 = [X1(1:2);0;X1(3:4);0;zeros(3,1)]; Xgt(:,m) = X1; end % Constant acceleration first = true; for m = seg2+1:numel(tt) X0 = Xgt(:,m-1); if first vel = X0(2:3:end); ua = vel/norm(vel); va = acc*ua; X0(3:3:end) = va; first = false; end X1 = constacc(X0, dt); Xgt(:,m) = X1; end % Drop acceleration dimension slct = eye(9); slct(3:3:end,:) = []; Xgt = slct*Xgt; figs = [figs figure]; plot(Xgt(1,1:seg1),Xgt(3,1:seg1),'.-'); hold on; plot(Xgt(1,seg1+1:seg2),Xgt(3,seg1+1:seg2),'.-'); plot(Xgt(1,seg2+1:end),Xgt(3,seg2+1:end),'.-'); grid on; xlabel('X Position (m)'); ylabel('Y Position (m)'); title('True Position') axis equal; legend('Constant Velocity', 'Constant Turn', 'Constant Acceleration') end