В этом примере вводится основная концепция микродоплеровского эффекта в радиолокационном возврате цели из-за вращения этой цели. Вы можете использовать микро-доплеровскую сигнатуру для идентификации цели.
Движущаяся цель вносит сдвиг частоты в возврат радара из-за эффекта Доплера. Однако, поскольку большинство целей не являются жесткими телами, в разных частях цели в дополнение к движению платформы часто присутствуют другие вибрации и вращения. Например, когда вертолет летает, его лопасти вращаются или когда человек ходит, его руки качаются естественным образом. Эти микромасштабные движения производят дополнительные доплеровские сдвиги, называемые микродоплеровскими эффектами, которые полезны для идентификации целевых признаков. В этом примере показаны два применения, в которых могут быть полезны микродоплеровские эффекты. В первой заявке микродоплеровские сигнатуры используются для определения скорости лопасти вертолета. Во второй заявке микродоплеровские сигнатуры используются для идентификации пешехода при возврате автомобильного радара.
Рассмотрим вертолет с четырьмя лопастями ротора. Предположим, что радар расположен в начале координат. Укажите местоположение вертолета как (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
Это изображение также показывает скорость Vt наконечника, которая может быть получена из максимального доплеровского значения. Максимальная доплеровская величина находится на расстоянии около 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 метра. Такая информация, как количество лопастей, длина лопастей и скорость вращения лопастей, может помочь идентифицировать модель вертолета.
Рассматриваем эго-автомобиль с автомобильной радиолокационной системой 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);
На следующем рисунке показана доплеровская карта дальности, сформированная с помощью радара эго-автомобиля с течением времени. Поскольку припаркованный автомобиль является гораздо более сильной мишенью, чем пешеход, пешеход легко затеняется припаркованным автомобилем на карте дальности - доплеровской. В результате на карте всегда отображается одна цель.
Это означает, что традиционная обработка не может удовлетворить наши потребности в этой ситуации.
Микро-доплеровский эффект во временной частотной области может быть хорошим кандидатом для идентификации наличия сигнатуры пешехода, встроенной в радиолокационный сигнал. В качестве примера, следующий раздел моделирует возврат радара в течение 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');
По кривой видно, что есть примерно четыре области. Регион А представляет собой наиболее значительный вклад в сигнал, которым является припаркованный автомобиль. Область D представляет шум. Поэтому области В и С обусловлены смешиванием возвращения припаркованного автомобиля и возвращения пешехода. Потому что возвращение от пешехода гораздо слабее, чем возвращение от припаркованного автомобиля. В области В его еще можно замаскировать остатком возврата из припаркованного автомобиля. Поэтому мы выбираем область 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] Чен, В.К., Микро-доплеровский эффект в радаре, Artech House, 2011
[2] Чен, V. C., Ф. Ли, судно. Хо и Х. Вехслер, «Микро-доплеровский эффект в радаре: исследование феномена, модели и моделирования», IEEE Transactions on Aerospace and Electronic Systems, Vol 42, No. 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