Введение

Движущаяся цель вводит сдвиг частоты в возвращении радара из-за эффекта Доплера. Однако, поскольку большинство целей не являются твердыми телами, часто существуют другие вибрации и вращения в разных частях цели в дополнение к движению платформы. Например, когда вертолет летит, его блейды вращаются, или когда человек ходит, его руки естественно качаются. Эти движения микро шкалы вызывают дополнительные доплеровские сдвиги, называемые микро-допплеровскими эффектами, которые полезны при идентификации целевых функций. Этот пример показывает два применения, где микродопплеровские эффекты могут быть полезными. В первой заявке для определения скорости блейда вертолета используются микродопплеровские сигнатуры. Во второй заявке микродопплеровские сигнатуры используются для идентификации пешехода при возврате автомобильного радара.

Оценка скорости блейда вертолета

Рассмотрим вертолет с четырьмя блейдами ротора. Предположим, что радар расположен в источник. Укажите местоположение вертолета как (500, 0, 500), которое устанавливает его расстояние от радара в метрах и скорость (60, 0, 0) м/с.

radarpos = [0;0;0];
radarvel = [0;0;0];

tgtinitpos = [500;0;500];
tgtvel     = [60;0;0];
tgtmotion  = phased.Platform('InitialPosition',tgtinitpos,'Velocity',tgtvel);

В этой симуляции вертолет моделируется пятью рассеивателями: центром вращения и советами четырёх блейдов. Центр вращения перемещается вместе с корпусом вертолета. Каждый совет блейда отстоит на 90 степени от совета его соседних блейдов. Блейды вращаются с постоянной скоростью 4 оборота в секунду. Длина рычага каждого блейда составляет 6,5 метра.

Nblades   = 4;
bladeang  = (0:Nblades-1)*2*pi/Nblades;
bladelen  = 6.5;
bladerate = deg2rad(4*360);  % rps -> rad/sec

Все четыре советов блейдов приняты с одинаковыми отражательными способностями, в то время как отражательная способность для центра вращения сильнее.

c  = 3e8;
fc = 5e9;
helicop = phased.RadarTarget('MeanRCS',[10 .1 .1 .1 .1],'PropagationSpeed',c,...
    'OperatingFrequency',fc);

Симуляция эха вертолета

Предположим, что радар работает на 5 ГГц с простым импульсом. Частота повторения импульсов составляет 20 кГц. Для простоты предположим, что сигнал распространяется в свободном пространстве.

fs     = 1e6;
prf    = 2e4;
lambda = c/fc;

wav = phased.RectangularWaveform('SampleRate',fs,'PulseWidth',2e-6,'PRF',prf);
ura = phased.URA('Size',4,'ElementSpacing',lambda/2);
tx  = phased.Transmitter;
rx  = phased.ReceiverPreamp;
env = phased.FreeSpace('PropagationSpeed',c,'OperatingFrequency',fc,...
    'TwoWayPropagation',true,'SampleRate',fs);
txant = phased.Radiator('Sensor',ura,'PropagationSpeed',c,'OperatingFrequency',fc);
rxant = phased.Collector('Sensor',ura,'PropagationSpeed',c,'OperatingFrequency',fc);

На каждом импульсе вертолет движется по своей траектории. При этом блейды продолжают вращаться, а советы блейдов создают дополнительное перемещение и угловую скорость.

NSampPerPulse = round(fs/prf);
Niter = 1e4;
y     = complex(zeros(NSampPerPulse,Niter));
rng(2018);
for m = 1:Niter
    % update helicopter motion
    t = (m-1)/prf;
    [scatterpos,scattervel,scatterang] = helicopmotion(t,tgtmotion,bladeang,bladelen,bladerate);

    % simulate echo
    x  = txant(tx(wav()),scatterang);                    % transmit
    xt = env(x,radarpos,scatterpos,radarvel,scattervel); % propagates to/from scatterers
    xt = helicop(xt);                                    % reflect
    xr = rx(rxant(xt,scatterang));                       % receive
    y(:,m) = sum(xr,2);                                  % beamform
end

Этот рисунок показывает доплеровскую характеристику диапазона, используя первые 128 импульсов принимаемого сигнала. Вы можете увидеть отображение трёх возвратов в целевой области значений приблизительно 700 метров.

rdresp  = phased.RangeDopplerResponse('PropagationSpeed',c,'SampleRate',fs,...
    'DopplerFFTLengthSource','Property','DopplerFFTLength',128,'DopplerOutput','Speed',...
    'OperatingFrequency',fc);
mfcoeff = getMatchedFilter(wav);
plotResponse(rdresp,y(:,1:128),mfcoeff);
ylim([0 3000])

Хотя возвраты выглядят так, как будто они от разных целей, они на самом деле все от одной цели. Возврат центра происходит из центра вращения и намного сильнее по сравнению с двумя другими возвратами. Эта интенсивность обусловлена тем, что отражение от корпуса вертолета сильнее по сравнению с советами блейдов. На график показана скорость -40 м/с для центра вращения. Это значение соответствует истинной радиальной скорости цели.

tgtpos = scatterpos(:,1);
tgtvel = scattervel(:,1);
tgtvel_truth = radialspeed(tgtpos,tgtvel,radarpos,radarvel)

tgtvel_truth =

  -43.6435

Другие два возврата из советов блейдов, когда они приближаются или отходят от цели на максимальной скорости. По графику скорости, соответствующие этим двум обнаружениям сближения и отхода, составляют около 75 м/с и -160 м/с соответственно.

maxbladetipvel = [bladelen*bladerate;0;0];
vtp = radialspeed(tgtpos,-maxbladetipvel+tgtvel,radarpos,radarvel)
vtn = radialspeed(tgtpos,maxbladetipvel+tgtvel,radarpos,radarvel)

vtp =

   75.1853


vtn =

 -162.4723

Можно связать все три обнаружения с одной и той же целью посредством дальнейшей обработки, но эта тема выходит за возможности этого примера.

Обратный микро-допплеровский анализ блейда

Частотно-временное представление микродопплеровских эффектов может выявить больше информации. Этот код создает представление частота-время в обнаруженном целевом интервале области значений.

mf  = phased.MatchedFilter('Coefficients',mfcoeff);
ymf = mf(y);
[~,ridx] = max(sum(abs(ymf),2)); % detection via peak finding along range
pspectrum(ymf(ridx,:),prf,'spectrogram')

Рисунок показывает микродопплеровскую модуляцию, вызванную советами блейда вокруг постоянного доплеровского сдвига. Изображение предполагает, что каждый совет блейда вводит синусоидально-подобную доплеровскую модуляцию. Как отмечается на рисунке ниже, в пределах каждого периода синусоиды присутствуют три дополнительные синусоиды, появляющиеся на равном расстоянии. Этот внешний вид говорит о том, что вертолет оборудован четырьмя равноудаленными блейдами.

hanno = helperAnnotateMicroDopplerSpectrogram(gcf);

В дополнение к количеству блейдов, изображение также показывает, что период каждой синусоиды, Tr, составляет около 250 мс. Это значение означает, что блейд возвращается в свое исходное положение после 250 мс. В этом случае угловая скорость вертолета составляет около 4 оборотов в секунду, что соответствует параметру симуляции.

Tp = 250e-3;
bladerate_est = 1/Tp

bladerate_est =

     4

Это изображение также показывает скорость совет, которая может быть получена из максимального Доплера. Максимальный Доплер находится примерно на 4 кГц от постоянного Доплера, введенного массовым движением. Вычислите обнаруженную максимальную скорость совета.

Vt_detect = dop2speed(4e3,c/fc)/2

Vt_detect =

   120

Это значение является максимальной скоростью совета вдоль радиального направления. Чтобы получить правильную максимальную скорость совета, должна быть принята во фактор относительная ориентация. Поскольку блейды вращаются по кругу, на обнаружение не влияет угол азимута. Исправьте только угол возвышения для результата максимальной скорости совета.

doa = phased.MUSICEstimator2D('SensorArray',ura,'OperatingFrequency',fc,...
    'PropagationSpeed',c,'DOAOutputPort',true,'ElevationScanAngles',-90:90);
[~,ang_est] = doa(xr);
Vt_est = Vt_detect/cosd(ang_est(2))

Vt_est =

  164.0793

На основе скорректированной максимальной скорости совета и скорости вращения лезвия вычислите длину блейда.

bladelen_est = Vt_est/(bladerate_est*2*pi)

bladelen_est =

    6.5285

Обратите внимание, что результат соответствует параметру симуляции 6,5 метров. Информация, такая как количество блейдов, длина блейда и скорость вращения блейда, может помочь идентифицировать модель вертолета.

Идентификация пешеходов в автомобильном радаре

Рассматривая ego автомобиль с автомобильной радиолокационной системой FMCW, ширина полосы пропускания которой 250 МГц и работает на 24 ГГц.

bw = 250e6;
fs = bw;
fc = 24e9;
tm = 1e-6;
wav = phased.FMCWWaveform('SampleRate',fs,'SweepTime',tm,...
    'SweepBandwidth',bw);

Эго-машина едет по дороге. По пути стоит припаркованный со стороны улицы автомобиль и за автомобилем выходит человек. Сцена проиллюстрирована на следующей схеме

Основываясь на этой настройке, если эго-автомобиль не может идентифицировать присутствие пешехода, может произойти авария.

egocar_pos = [0;0;0];
egocar_vel = [30*1600/3600;0;0];
egocar = phased.Platform('InitialPosition',egocar_pos,'Velocity',egocar_vel,...
    'OrientationAxesOutputPort',true);

parkedcar_pos = [39;-4;0];
parkedcar_vel = [0;0;0];
parkedcar = phased.Platform('InitialPosition',parkedcar_pos,'Velocity',parkedcar_vel,...
    'OrientationAxesOutputPort',true);
parkedcar_tgt = phased.RadarTarget('PropagationSpeed',c,'OperatingFrequency',fc,'MeanRCS',10);

ped_pos = [40;-3;0];
ped_vel = [0;1;0];
ped_heading = 90;
ped_height = 1.8;
ped = phased.BackscatterPedestrian('InitialPosition',ped_pos,'InitialHeading',ped_heading,...
    'PropagationSpeed',c,'OperatingFrequency',fc,'Height',1.6,'WalkingSpeed',1);

chan_ped = phased.FreeSpace('PropagationSpeed',c,'OperatingFrequency',fc,...
    'TwoWayPropagation',true,'SampleRate',fs);
chan_pcar = phased.FreeSpace('PropagationSpeed',c,'OperatingFrequency',fc,...
    'TwoWayPropagation',true,'SampleRate',fs);

tx = phased.Transmitter('PeakPower',1,'Gain',25);
rx = phased.ReceiverPreamp('Gain',25,'NoiseFigure',10);

Экстракция микро-допплера пешеходов

Следующий рисунок показывает карту Доплера диапазона, сгенерированную с радара эго-автомобиля с течением времени. Поскольку припаркованный автомобиль является гораздо более сильной мишенью, чем пешеходный, пешехода легко затеняет припаркованный автомобиль на карте range-Doppler. В результате карта всегда показывает одну цель.

Это означает, что обычная обработка не может удовлетворить наши потребности в этой ситуации.

Эффект Микро-Допплера во времени частотного диапазона может быть хорошим кандидатом, чтобы идентифицировать, есть ли пешеходная сигнатура, встроенная в радиолокационный сигнал. В качестве примера в следующем разделе моделируется возврат радара в течение 2,5 секунд.

Tsamp = 0.001;
npulse = 2500;
xr = complex(zeros(round(fs*tm),npulse));
xr_ped = complex(zeros(round(fs*tm),npulse));
for m = 1:npulse
    [pos_ego,vel_ego,ax_ego] = egocar(Tsamp);
    [pos_pcar,vel_pcar,ax_pcar] = parkedcar(Tsamp);
    [pos_ped,vel_ped,ax_ped] = move(ped,Tsamp,ped_heading);
    [~,angrt_ped] = rangeangle(pos_ego,pos_ped,ax_ped);
    [~,angrt_pcar] = rangeangle(pos_ego,pos_pcar,ax_pcar);

    x = tx(wav());
    xt_ped = chan_ped(repmat(x,1,size(pos_ped,2)),pos_ego,pos_ped,vel_ego,vel_ped);
    xt_pcar = chan_pcar(x,pos_ego,pos_pcar,vel_ego,vel_pcar);
    xt_ped = reflect(ped,xt_ped,angrt_ped);
    xt_pcar = parkedcar_tgt(xt_pcar);
    xr_ped(:,m) = rx(xt_ped);
    xr(:,m) = rx(xt_ped+xt_pcar);
end

xd_ped = conj(dechirp(xr_ped,x));
xd = conj(dechirp(xr,x));

В имитированном сигнале xd_ped содержит только возврат пешехода в то время как xd содержит возврат как пешехода, так и припаркованного автомобиля. Если сгенерировать спектрограмму, используя только возврат пешехода, получаем график, показанный ниже.

clf;
spectrogram(sum(xd_ped),kaiser(128,10),120,256,1/Tsamp,'centered','yaxis');
clim = get(gca,'CLim');
set(gca,'CLim',clim(2)+[-50 0])

Обратите внимание, что качели рук и ног производят много параболических кривых во временном частотном диапазоне вдоль пути. Поэтому такие функции могут использоваться, чтобы определить, существует ли в сцене пешехода.

Однако, когда мы генерируем спектрограмму непосредственно от общего возврата, мы получаем следующий график.

spectrogram(sum(xd),kaiser(128,10),120,256,1/Tsamp,'centered','yaxis');
clim = get(gca,'CLim');
set(gca,'CLim',clim(2)+[-50 0])

Мы замечаем, что возврат припаркованного автомобиля продолжает доминировать в возврате, даже в частотный диапазон времени. Поэтому частотная характеристика показывает только Допплер относительно припаркованного автомобиля. Падение частоты Допплера происходит из-за того, что эго-автомобиль приближается к припаркованному автомобилю и относительная скорость падает к 0.

Чтобы увидеть, есть ли возврат, скрытый за сильной отдачей, мы можем использовать сингулярное разложение. Следующий график показывает распределение сингулярного значения дешифрованных импульсов.

[uxd,sxd,vxd] = svd(xd);
clf
plot(10*log10(diag(sxd)));
xlabel('Rank');
ylabel('Singular Values');
hold on;
plot([56 56],[-40 10],'r--');
plot([100 100],[-40 10],'r--');
plot([110 110],[-40 10],'r--');
text(25,-10,'A');
text(75,-10,'B');
text(105,-10,'C');
text(175,-10,'D');

Из кривой ясно, что существует примерно четыре области. Самый значительный вклад в сигнал представляет области A, которым является припаркованный автомобиль. Область D представляет шум. Поэтому области B и C обусловлены смешением возврата припаркованного автомобиля и возвращения пешехода. Потому что возврат от пешехода намного слабее, чем возврат от припаркованной машины. В области B его еще можно замаскировать остатком возврата из припаркованной машины. Поэтому мы выбираем область C, чтобы восстановить сигнал, и затем снова строим график частотной характеристики времени.

rk = 100:110;
xdr = uxd(:,rk)*sxd(rk,:)*vxd';
clf
spectrogram(sum(xdr),kaiser(128,10),120,256,1/Tsamp,'centered','yaxis');
clim = get(gca,'CLim');
set(gca,'CLim',clim(2)+[-50 0])

При успешно отфильтрованном возврате из автомобиля появляется микродопплеровская подпись от пешехода. Поэтому можно сделать вывод, что в месте происшествия есть пешеходы и действовать соответствующим образом, чтобы избежать ДТП.

Сводные данные

Этот пример вводит базовую концепцию эффекта микродопплера и показывает его влияние на целевой возврат. Это также показывает, как извлечь микродопплеровскую сигнатуру из принятого сигнала I/Q, а затем вывести соответствующие целевые параметры из информации микродопплера.

Ссылки

[1] Chen, V. C., The Micro-Doppler Effect in Radar, Artech House, 2011

[2] Chen, V. C., F. Li, S.-S. Ho, and H. Wechsler, «Micro-Doppler Effect in Radar: Phenomenon, модель Simulation Study», Транзакции IEEE по аэрокосмическим и электронным системам, том 42, № 1, январь 2006 г.

Служебные функции

Функциональные helicopmotion моделирует движение нескольких рассеивателей вертолета.

function [scatterpos,scattervel,scatterang] = helicopmotion(...
    t,tgtmotion,BladeAng,ArmLength,BladeRate)

prf = 2e4;
radarpos = [0;0;0];
Nblades  = size(BladeAng,2);

[tgtpos,tgtvel] = tgtmotion(1/prf);

RotAng     = BladeRate*t;
scatterpos = [0 ArmLength*cos(RotAng+BladeAng);0 ArmLength*sin(RotAng+BladeAng);zeros(1,Nblades+1)]+tgtpos;
scattervel = [0 -BladeRate*ArmLength*sin(RotAng+BladeAng);...
    0 BladeRate*ArmLength*cos(RotAng+BladeAng);zeros(1,Nblades+1)]+tgtvel;

[~,scatterang] = rangeangle(scatterpos,radarpos);

end