微多普勒效应建模¶
该示例介绍了微多普勒效应的概念，并演示了如何构建一个能够从反射信号中提取微多普勒的雷达系统。这一特征可进一步用于解决目标分类问题。
连接库和函数¶
In [ ]:
let
installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["DSP","FFTW","LinearAlgebra"]
for pack in list_packages
pack in installed_packages || Pkg.add(pack)
end
end
In [ ]:
using DSP,FFTW,LinearAlgebra
In [ ]:
gr()
default(titlefontsize=12,top_margin=5Plots.px,guidefont=10,
fontfamily = "Computer Modern",colorbar_titlefontsize=8,size=(800,400),
margin = 2Plots.mm
)
function 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)
return scatterpos,scattervel,scatterang
end
function spectrogramm_ped(x,t_step;down_lim = -35)
out,f,t = EngeePhased.Functions.spectrogram(sum(x;dims=1)[:];
window = DSP.kaiser(128,10),
noverlap=120,
nfft=256,fs=1/t_step,
centered=true,freqloc = "yaxis",out = :data
)
out_sig = map(x -> x < down_lim ? down_lim : x,DSP.pow2db.(abs.(out)))
fig = Plots.heatmap(
t[:],f,out_sig,color= :jet,
gridalpha=0.5,
colorbar_title = "Мощность, дБВт",
xticks=0:0.5:3,
yticks=-500:100:500,
)
Plots.xlabel!("Время, с")
Plots.ylabel!("Частота Доплера, Гц")
return fig
end
function calc_spectrogram_quad(Data,ts,T;limit = -150)
x = Data./maximum(abs.(Data))
dT = T/length(x)
F = 1/dT
spec_quad = DSP.spectrogram(x,128,120;window=DSP.kaiser(128,3.95),nfft=512,fs=F)
power_quad = map(x -> x < limit ? limit : x,pow2db.(spec_quad.power))
power_quad_quad_new = zeros(size(power_quad))
power_quad_quad_new[1:round(Int64,size(power_quad,1)/2),:] .= power_quad[round(Int64,size(power_quad,1)/2)+1:end,:]
power_quad_quad_new[round(Int64,size(power_quad,1)/2)+1:end,:] .= power_quad[1:round(Int64,size(power_quad,1)/2),:]
fig = Plots.heatmap(spec_quad.time*1e3, fftshift(spec_quad.freq)*1e-3, power_quad_quad_new,color= :jet,
gridalpha=0.5,ylims=(-10,10),xlims=(0,500),colorbar_title = "Мощность, дБВт",
xticks=0:50:500,yticks=-10:2:10)
Plots.xlabel!("Время, мс")
Plots.ylabel!("Частота Доплера, кГц")
Plots.title!("Спектрограмма отраженного сигнала от вертолета")
return fig
end
function helperAnnotateSpectrogram(fig)
Plots.plot!(fig,[20;50],[4;7],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(50, 7.5, ("Лопасть 1", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[70;120],[4;5.5],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(120, 6, ("Лопасть 2", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[20;50],[-4;-7],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(50, -7.5, ("Лопасть 3", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[70;120],[-4;-5.5],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(120, -6, ("Лопасть 4", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[0;250],[0;-0.1],linewidth=3,c=:white,lab="",arrow=Plots.Arrow(:close, :both, 0.5, 1))
Plots.annotate!(120, -1, ("Период", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[124;125],[0.2;4],linewidth=3,c=:white,lab="",arrow=Plots.Arrow(:close, :both, 0.5, 1))
Plots.annotate!(137, 2, ("Скорость вращения", 10, :white, :left),annotationfontfamily="Roboto")
end
function plotting_svd(x)
Plots.plot(10*log10.(x),lab="",ylabel="Сингулярные значения",xlabel="ранг")
Plots.plot!([55], seriestype="vline", label="",ls=:dashdot,color="red",lw=2)
Plots.plot!([101], seriestype="vline", label="",ls=:dashdot,color="red",lw=2)
Plots.plot!([109], seriestype="vline", label="",ls=:dashdot,color="red",lw=2)
Plots.annotate!(25, -10, ("A", 10, :black, :left),annotationfontfamily="Roboto")
Plots.annotate!(75, -10, ("B", 10, :black, :left),annotationfontfamily="Roboto")
Plots.annotate!(102, -10, ("C", 10, :black, :left),annotationfontfamily="Roboto")
Plots.annotate!(175, -10, ("D", 10, :black, :left),annotationfontfamily="Roboto")
end
Out[0]:
1.微多普勒效应的概念¶
微运动*的概念与微多普勒效应密不可分：
微动是指物体某一组成部分的微小运动或偏差，是整个物体主要体积运动之外的额外运动[1]。微动的来源可以是振动的表面、直升机旋翼的旋转叶片、行走的人做周期性的手脚运动、鸟类拍打翅膀或其他原因。
微多普勒效应是一种与反射信号频谱成分丰富有关的现象，在雷达研究过程中由于物体的微动而产生。
考虑物体运动的两种情况：简单平移运动（上图）和复杂平移-旋转运动（下图）。
在这种情况下，物体是匀速向我们运动的，因此多普勒频率具有恒定值，与时间无关。
第二种情况则不同：旋转运动（如直升机叶片的旋转）在反射信号的频谱中引入了周期性变化的成分，这就是微多普勒效应的特征。
频谱中包含微多普勒成分的反射信号频谱图称为微多普勒特征（或缩写为MDS）。鸟类 MDS 的示例如下：
2.行人 MDS 建模¶
在下面的模拟中，我们考虑开发一种基于连续 LFM 信号的汽车雷达，以识别道路上的人行横道。
2.1 模型场景描述¶
模型考虑了以下情景：
- 一辆装有雷达的汽车在内侧车道上以每小时 30 公里的速度从起点驶向一辆站立的汽车和一名行人；
- 站立的汽车距离行人 39.4 米。
- 行人的速度为 1 米/秒。
该场景的直观图如下所示：
2.2 模型参数¶
穿越识别需要较高的范围和速度分辨率。汽车雷达使用连续线性频率调制（FM）信号来满足这一要求。
In [ ]:
bw = 250e6 # полоса сигнала, Гц
fs = bw # частота дискретизации, Гц
fc = 24e9 # несущая частота, Гц
tm = 1e-6 # время развертки, с
c = 3e8 # скорость света, м/с
egocar_pos = [0;0;0] # м, вектор положения автомобиля с радаром
egocar_vel = [30*1600/3600;0;0] # м/с, вектор скорости автомобиля с радаром
parkedcar_pos = [39;-4;0] # м, вектор положения впереди стоящего автомобиля
parkedcar_vel = [0;0;0] # м/с, вектор скорости впереди стоящего автомобиля
ped_pos = [40;-3;0] # м, вектор положения пешехода
ped_vel = [0;1;0] # м, вектор скорости пешехода
ped_heading = 90 # град, начальное направление движения пешехода
ped_height = 1.8; # рост пешехода
2.3 系统目标形成¶
为了形成物体的反射信号--行人使用系统对象EngeeRadars.BackscatterPedestrian
In [ ]:
wav_fmcw = EngeePhased.FMCWWaveform(
SampleRate=fs, # частота дискретизации
SweepTime=tm, # время развертки
SweepBandwidth=bw # полоса сигнала
)
# СО автомобиля с радаром
egocar = EngeePhased.Platform(
InitialPosition=egocar_pos, # начальное положение
Velocity=egocar_vel, # вектор скорости
OrientationAxesOutputPort=true # включение выходного порта ориентации в пространстве
)
# Формирование объекта впереди стоящего автомобиля
parkedcar = EngeePhased.Platform(
InitialPosition=parkedcar_pos, # начальное положение
Velocity=parkedcar_vel, # скорость
OrientationAxesOutputPort=true # включение выходного порта ориентации в пространстве
)
parkedcar_tgt = EngeePhased.RadarTarget(PropagationSpeed=c,OperatingFrequency=fc,MeanRCS=10)
# Формирование объекта - пешеход
ped = EngeeRadars.BackscatterPedestrian(
InitPosition=ped_pos, # начальное положение
InitHeading=ped_heading, # начальное направление движения
PropagationSpeed=c, # скорость распространения сигнала
OperatingFrequency=fc, # несущая частота
Height=ped_height, # рост
WalkingSpeed=1 # скорость движения
)
# каналы распространения: до пешехода и ближайшего автомобиля
chan_ped = EngeePhased.FreeSpace(
PropagationSpeed=c, # скорость распространения
OperatingFrequency=fc, # несущая частота
TwoWayPropagation=true, # учет двунаправленного распространения
SampleRate=fs # частота дискретизации
)
chan_pcar = EngeePhased.FreeSpace(
PropagationSpeed=c, # скорость распространения
OperatingFrequency=fc, # несущая частота
TwoWayPropagation=true, # учет двунаправленного распространения
SampleRate=fs # частота дискретизации
)
# приемник и передатчик
tx = EngeePhased.Transmitter(PeakPower=1,Gain=25)
rx = EngeePhased.ReceiverPreamp(Gain=25,NoiseFigure=10);
2.4 模型模拟¶
要正确表示行人 MDS，就必须累积至少一个人体运动周期内的反射信号。行人的速度为 1 米/秒，因此一个完整的周期需要 2 秒钟。让我们以 1 毫秒为单位，模拟 2.5 秒的行人运动场景：
In [ ]:
Tsim = 2.5 # время симуляции
Tsamp = 0.001 # шаг расчета симуляции, с
npulse = round(Int64,Tsim/Tsamp) # количество импульсов
# Выделение памяти под:
xr = complex(zeros(round(Int64,fs*tm),npulse)) # суммарный отраженный сигнал
xr_ped = complex(zeros(round(Int64,fs*tm),npulse)) # отраженный сигнал от пешехода
@inbounds 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)
# расчет отраженного сигнала
global x = tx(wav_fmcw()) # зондирующий сигнал
# прохождение сигналов в канале прямого и обратного распространения
xt_ped = chan_ped(repeat(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));
2.5 处理接收到的信号¶
在模拟信号中，xd_ped 只包含行人的响应，而xd 则包含行人和停放汽车的响应。如果我们仅使用反射的行人信号来创建频谱图，则可得到以下 MDS：
In [ ]:
spectrogramm_ped(xd_ped,Tsamp;down_lim = -35)
title!("Спектрограмма отраженного сигнала от пешехода")
Out[0]:
行人和汽车总反射信号的频谱图如下：
In [ ]:
spectrogramm_ped(xd,Tsamp;down_lim = -35)
title!("Спектрограмма отраженного сигнала от пешехода и автомобиля")
Out[0]:
可以看出，汽车的 MDS 在行人的背景 MDS 中占主导地位。因此，如果不进行额外处理，就无法检测到行人。
2.6 利用 svd 分解进行处理¶
要从较强的信号中检测出较弱的信号，我们可以使用奇异值分解，这样就能突出感兴趣的区域。
In [ ]:
out_svd = svd(xd,full=true) # расчет компонента разложения
uxd,sxd,vxd = out_svd.U,out_svd.S,out_svd.V; # считывание матриц S,U,V
使用函数plotting_svd 显示主对角奇异矩阵S ：
In [ ]:
plotting_svd(sxd)
Out[0]:
在奇异值与矩阵秩的关系图上，我们通常可以区分出 4 个区域：
- A - 包括对频谱图有重要贡献的主要信号；
- B - 行人和汽车信号的组合，以第二个信号为主
- C - 行人信号和汽车信号的组合，以前者为主；
- D - 接收设备的噪声。
因此，要有效提取行人 MDS，就必须从整个频谱图中提取位于等级值 103-110 内的 C 区域：
In [ ]:
rk = 103:110 # диапазон рангов матрицы сигнулярных значений
# расчет обработанной спектрограммы
sxd_new = [sxd;zeros(size(vxd,1)-length(sxd))]
xdr = uxd[:,rk]*Diagonal(sxd_new)[rk,:]*vxd';
我们再次绘制经过额外处理后的行人和汽车 MDS 图（为正确表示，我们将下限降至 -40 dB，因为主要功率已被去除）：
In [ ]:
spectrogramm_ped(xdr,Tsamp;down_lim = -40)
title!("МДС пешехода и автомобиля после svd преобразования")
Out[0]:
滤除汽车信号后，行人的 MDS 清晰可见，而汽车的 MDS 水平则明显下降。
因此，在额外处理（svd 分解）的帮助下，可以在附近汽车干扰信号的背景下突出行人的光谱成分。
3 直升机 MDS 建模¶
3.1 模型情景描述¶
考虑以下情况
- 一架直升机沿 x 轴以 100 米/秒的速度在 x-z 平面上以 [-25,500,500] mesres 的距离匀速飞行。
- 雷达静止不动，位于原点。
在这种情况下，直升机在运动开始时靠近雷达，在过渡到 x 轴正区域后开始远离。这种情况的可视化效果如下图所示：
3.2 模型参数¶
让我们来设置雷达和目标（直升机）移动场景的参数：
In [ ]:
radarpos = [0;0;0] # вектор положения РЛС, м
radarvel = [0;0;0] # вектор скорости РЛС, м/с
tgtinitpos = [-25;500;500] # начальное положение цели, м
tgtvel = [100;0;0]; # скорость цели, м/с
接下来，绘制直升机的几何图形，指定其主要参数（叶片数量、叶片长度、螺旋桨速度）：
In [ ]:
f_rot = 4 # частота вращения винта, c
Nblades = 4 # количество лопастей
bladeang = Matrix((0:Nblades-1)')*2*pi/Nblades # вектор углов лопастей
bladelen = 6.5 # длина лопасти, м
bladerate = deg2rad(360*f_rot); # скорость вращения винта рад/с
Out[0]:
下一阶段是组建雷达系统。设置系统的基本参数：
In [ ]:
c = 3e8 # скорость распространения сигнала, м/с
fc = 5e9 # несущая частота, Гц
fs = 1e6 # частота дискретизации, Гц
prf = 2e4 # частота следования импульсов, Гц
lambda = c/fc; # длина волны, м
3.3 创建系统对象¶
现在，利用目标场景和雷达系统的参数，让我们使用库EngeePhased 创建模型的系统对象：
In [ ]:
# Формирование цели
helicop = EngeePhased.RadarTarget(
MeanRCS=[10 .1 .1 .1 .1], # ЭПР элементов вертолета (фюзеляжа и 4 лопастей)
PropagationSpeed=c, # скорость распространения
OperatingFrequency=fc # несущая частота сигнала
)
tgtmotion = EngeePhased.Platform(InitialPosition=tgtinitpos,Velocity=tgtvel)
# Генератор зондирующего сигнала
wav = EngeePhased.RectangularWaveform(
SampleRate=fs, # частота дискретизации
PulseWidth=2/fs, # длительность импульса
PRF=prf # частота следования импульсов
)
# Формирование антенной решетки (прямоугольная 4х4)
ura = EngeePhased.URA(
Size=[4 4], # размер антенной решетки
ElementSpacing=[lambda/2 lambda/2] # расстояние между элементами
)
tx = EngeePhased.Transmitter() # передатчик
rx = EngeePhased.ReceiverPreamp() # предусилитель
# Cреда распространения
env = EngeePhased.FreeSpace(
PropagationSpeed=c, # скорость распространения, м/с
OperatingFrequency=fc, # несущая частота, Гц
TwoWayPropagation=true, # учет двунаправленного распространения
SampleRate=fs # частота дискретизации
)
# Передающая АР
txant = EngeePhased.Radiator(
Sensor=ura, # геометрия АР
PropagationSpeed=c, # скорость распространения, м/с
OperatingFrequency=fc # несущая частота, Гц
)
# Принимающая АР
rxant = EngeePhased.Collector(
Sensor=ura, # геометрия АР
PropagationSpeed=c, # скорость распространения, м/с
OperatingFrequency=fc # несущая частота, Гц
);
3.4 运行模型模拟¶
让我们通过循环调用系统对象来实现直升机移动时雷达运行场景的计算：
In [ ]:
NSampPerPulse = round(Int64,fs/prf) # количество отсчетов в 1 импульсе
Niter = Int64(1e4) # количество итераций (число излученных импульсов)
y = complex(zeros(NSampPerPulse,Niter)) # выделение памяти для отраженного сигнала
@inbounds for m in 1:Niter
# обновление положение и скорости вертолета
t = (m-1)/(prf) # шаг симуляции по времени
global scatterpos,scattervel,scatterang = helicopmotion(t,tgtmotion,bladeang,bladelen,bladerate,prf,radarpos)
# расчет отраженного сигнала
x = txant(tx(wav()),scatterang) # зондирующий сигнал
xt = env(x,radarpos,scatterpos,radarvel,scattervel) # сигнал, прошедший среду
xt = helicop(xt) # отраженный сигнал от вертолета
global xr = rx(rxant(xt,scatterang)) # принятый сигнал
y[:,m] = sum(xr;dims=2) # формирование луча
end
3.5 信号处理¶
利用系统对象RangeDopplerResponse ，我们实现了测距-多普勒处理：匹配测距滤波和基于速度 FFT 的频谱处理。
In [ ]:
rdresp = EngeePhased.RangeDopplerResponse(
PropagationSpeed=c, # скорость распространения сигнала
SampleRate=fs, # частота дискретизации
DopplerFFTLengthSource="Property", # Cпособ задания БПФ (в параметрах СО)
DopplerFFTLength=128, # длина БПФ
DopplerOutput="Speed", # Доплероский выход (скорость)
OperatingFrequency=fc # несущая частота
)
mfcoeff = getMatchedFilter(wav)
fig1 = plotResponse(rdresp,y[:,1:128],mfcoeff)
Plots.ylims!(0,3000)
Plots.plot!(fig1,title="Дальностно-доплероский портрет",xlabel="Скорость, м/с",
ylabel="Дальность, м",titlefontsize=12,top_margin=5Plots.px,guidefont=10,
colorbar_title="Мощность, дБВт",fontfamily = "Computer Modern",colorbar_titlefontsize=8)
Out[0]:
可以看出，测距-多普勒肖像捕捉到了 3 个分量：中央分量是直升机机身，侧面分量是桨叶的旋转。
让我们估算一下直升机的径向速度：
In [ ]:
tgtpos = scatterpos[:,1] # положение вертолета на момент окончания симуляции
tgtvel = scattervel[:,1] # скорость вертолета на момент окончания симуляции
tgtvel_truth = radialspeed(tgtpos,tgtvel,radarpos,radarvel)
println("радиальная скорость фюзеляжа вертолета $(round(tgtvel_truth;sigdigits=3)) м/с")
我们还可以看到螺旋桨径向速度的最大值和最小值：
In [ ]:
maxbladetipvel = [bladelen*bladerate;0;0] # скорость лопастей винта
vtp = radialspeed(tgtpos,-maxbladetipvel+tgtvel,radarpos,radarvel)
vtn = radialspeed(tgtpos,maxbladetipvel+tgtvel,radarpos,radarvel)
println("максимальное значение радиальной скорости $(round(vtp;sigdigits=3)) м/с")
println("минимальное значение радиальной скорости $(round(vtn;sigdigits=3)) м/с")
让我们执行匹配滤波并绘制反射信号的频谱图：
In [ ]:
# Cогласованная фильтрация
mf = EngeePhased.MatchedFilter(Coefficients=mfcoeff) # согласованный фильтр
ymf = mf(y) # обработка
ridx = findmax(sum(abs.(ymf);dims=2)[:])[2]; # обнаружение пика по дальности
# Построение МДС
T = Niter/prf # длительность спектрограммы
fig1 = calc_spectrogram_quad(ymf[ridx,:],1/prf,T;limit=-100)
Out[0]:
中心成分表示直升机机身的反射信号，周期成分表示叶片的旋转。使用函数helperAnnotateSpectrogram 可视化 MDS 分量：
In [ ]:
helperAnnotateSpectrogram(fig1)
Out[0]:
从频谱图中我们可以估算出振荡周期，大约等于 250 毫秒。在这种情况下，叶片数量的估计值等于：
In [ ]:
Tp = 250e-3 # период МДС
bladerate_est = 1/Tp # оценка количества лопастей
println("Оценка числа лопастей $(round(bladerate_est;sigdigits=2))")
多普勒频率的最大偏差约等于 4 kHz，则直升机叶片的径向速度估计为
In [ ]:
f_doppler = 4e3 # амплитуда доплероской частоты, Гц
Vt_detect = dop2speed(f_doppler,lambda)/2
println("Радиальная скорость лопастей $(Vt_detect) м/с")
为了估算真实速度，我们需要确定桨叶位置角问题的解。让我们使用在系统对象MUSICEstimator2D 中实现的算法 "音乐"：
In [ ]:
doa = EngeePhased.MUSICEstimator2D(
SensorArray=ura, #
OperatingFrequency=fc, # несущая частота
PropagationSpeed=c, # скорость распространения сигнала
DOAOutputPort=true, # включение выхода углового направления
ElevationScanAngles=-90:90 # диапазон сканирования по углу места
)
_,ang_est = doa(xr) # расчет углов визирования
Vt_est = Vt_detect/cosd(ang_est[2]) # оценка скорости движения лопастей
println("Оценка скорости лопастей $(round(Vt_est;sigdigits=3)) м/с")
在估算出速度和叶片数量后，让我们来确定叶片的大致长度：
In [ ]:
bladelen_est = Vt_est/(bladerate_est*2*pi)
println("Оценка длины лопасти $(round(bladelen_est;sigdigits=3)) м")
长度的真实值为 6.5 米，因此误差为 4%。
结论¶
本示例考虑了 "微多普勒效应 "和微多普勒信号的概念。我们使用EngeePhased系统对象模拟了两种雷达情况下反射信号的微多普勒信号：存在空中目标（直升机）的地面雷达和存在行人和附近汽车的汽车雷达。