多普勒估计
该示例显示了一个单稳态脉冲雷达，它可以检测某些范围内移动目标的径向速度。 速度由运动目标引起的多普勒频移决定。 首先，我们将确定目标在给定范围内的存在，然后，使用多普勒处理，我们将估计目标在该范围内的径向速度。
使用的功能和库
Pkg.add(["DSP"])
using DSP,FFTW
function BasicMonostaticRadarExampleData()
c = 299_792_458 # 信号传播速度
# 创建系统对象:
# 矩形脉冲发生器
waveform = EngeePhased.RectangularWaveform(PulseWidth=100/c,
PRF=c/1e4,OutputFormat="Pulses",NumPulses=1,SampleRate= 2c/100);
# 天线元件(各向同性型)
antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[5e9 1.5e10]);
# 窄带信号发射器
radiator = EngeePhased.Radiator(Sensor=antenna,
PropagationSpeed=c,OperatingFrequency=1e10);
# 传输路径中的放大器
transmitter = EngeePhased.Transmitter(PeakPower=5226.4791642003665401716716587543,
Gain=20,LossFactor=0, InUseOutputPort=true,CoherentOnTransmit=true);
# 接收天线
collector = EngeePhased.Collector(Sensor=antenna,
PropagationSpeed=c,Wavefront="Plane",
OperatingFrequency=1e10);
# 前置放大器在接收路径
receiver = EngeePhased.ReceiverPreamp(Gain=20,NoiseFigure=0,
ReferenceTemperature=290,SampleRate=2*c/100, PhaseNoiseInputPort=false,
EnableInputPort=true,SeedSource="Property",Seed=1000);
# 雷达运动模型
sensormotion = EngeePhased.Platform(InitialPosition=[0;0;0],Velocity=[0;0;0]);
return waveform, antenna, radiator, collector, transmitter, receiver, sensormotion
end;
# 计算噪声功率的函数
function noisepow1(nbw::Number, nf::Number = 0, reftemp::Number = 290)
function systemp(nf::Number,reftemp::Number = 290)
nf >= 0 || error("Invalid input")
nfLinear = db2pow(nf);
stemp = reftemp*(nfLinear); # Kelvin
return stemp
end
B = 1.380649e-23
npow = B * systemp(nf,reftemp) * nbw
return npow
end;
配置雷达系统
首先，我们将定义雷达系统。 由于本示例侧重于多普勒处理，我们使用示例中构建的雷达系统Simulation of test signals for radar приемника。请读者使用此示例探索雷达系统的设计细节。
waveform,antenna,radiator,collector,transmitter,receiver,sensormotion = BasicMonostaticRadarExampleData();
系统建模
目标
多普勒处理使用由运动目标引起的多普勒频移。 现在我们将设置三个目标，指定它们的位置，有效散射面积（EPR）和速度。
tgtpos = [[1200; 1600; 0] [3543.63; 0; 0] [1600; 0; 1200]]; # 目标初始位置的组成部分
tgtvel = [[60; 80; 0] [0;0;0] [0; 100; 0]]; # 目标的速度组件
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # 目标移动模型
tgtrcs = [1.3 1.7 2.1]; # EPR目标
fc = radiator.OperatingFrequency;
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc);
请注意，第一个目标位于3,550，1,600和1,200m的距离处，并且以60，130和0m/s的速度移动。 第三目标不动。
分布环境
我们还需要为每个目标设置一个分发环境。 由于我们使用的是单稳态雷达，我们将定义一个双向传播模型。:
fs = waveform.SampleRate; # 系统的采样率
# 分销渠道模型
channel = EngeePhased.FreeSpace(
SampleRate=fs, # 系统的采样率
TwoWayPropagation=true, # 分布模式
OperatingFrequency=fc); # 信号的载波频率
信号合成
在确定了雷达系统、环境和目标后，我们将模拟从目标反射的接收回波信号。 得到的数据是每列时间快（每个脉冲内的时间）和每行时间慢（脉冲间的时间）的元素矩阵。
为了重现相同的结果，我们需要设置用于在接收器中产生噪声的初始数据。
release!(receiver) # 重置对象的内部状态（修改）
receiver.SeedSource = "Property"; # 更改噪声生成模型的类型
receiver.Seed = 2024; # 设置热噪声发生器的初始状态
prf = waveform.PRF; # 脉冲重复率
num_pulse_int = 10; # 累积脉冲数
fast_time_grid = (0:1/fs:1/prf-1/fs) # 快速时间网格
slow_time_grid = (0:num_pulse_int-1)./prf; # 慢时间网格
# 为生成的数据分配内存
rxpulses = zeros(ComplexF64,length(fast_time_grid),num_pulse_int);
for m = 1:num_pulse_int
# 更新雷达位置和目标
sensorpos,sensorvel = sensormotion(1/prf);
tgtpos,tgtvel = tgtmotion(1/prf);
# 根据从雷达到目标的范围计算角度
tgtrng,tgtang = rangeangle(tgtpos,sensorpos);
# 从雷达到目标的信号传播仿真
pulse = waveform(); # 矩形脉冲发生器
txsig,txstatus = transmitter(pulse); # 发射信号放大器
txsig = radiator(txsig,tgtang); # 周三生成信号的发射器
txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel); # 分布环境
# 从目标的信号的反射，考虑到EPR
tgtsig = target(txsig);
# 反射信号的接收和预放大
rxsig = collector(tgtsig,tgtang); # 接收天线
rxpulses[:,m] = receiver(rxsig,xor.(txstatus.>0,1)); # 前置放大器
end
多普勒估计
范围检测
要估计目标的多普勒频移，首先需要使用范围检测来确定它们的位置。 由于多普勒频移将信号功率分布在I和Q通道上，因此需要使用信号能量进行检测。 这意味着必须使用非相干检测方案。
在上面的例子中详细描述了检测过程，因此在这里我们将简单地执行必要的动作来估计目标范围。
# 设置虚警概率(pfa)
pfa = 1e-6;
# 在所考虑的系统中，噪声带等于采样频率的一半。
noise_bw = receiver.SampleRate/2; # 噪音波段
npower = noisepow1(noise_bw,
receiver.NoiseFigure,receiver.ReferenceTemperature); # 噪音功率
threshold = npower * db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # 检测阈值
# 应用一致过滤并更新检测阈值
matchingcoeff = getMatchedFilter(waveform); # 议定的系数的计算
# 基于细化信号的滤波器
matchedfilter = EngeePhased.MatchedFilter( # 创建一致的筛选器对象
Coefficients=matchingcoeff, # 滤波器系数
GainOutputPort=true); # 增益输出端口
rxpulses, mfgain = matchedfilter(rxpulses); # 一致过滤后的响应
threshold = threshold * db2pow(mfgain); # 新的检测阈值的计算
# 匹配滤波器中的延迟补偿
matchingdelay = size(matchingcoeff,1)-1; # 匹配滤波器的长度
rxpulses = buffer(rxpulses[(matchingdelay+1):end],size(rxpulses,1));
# 临时自动增益控制(VAR)的应用
# 以补偿依赖于范围的损失。
prop_speed = radiator.PropagationSpeed; # 信号传播速度
range_gates = prop_speed*fast_time_grid/2; # 范围分辨率
lambda = prop_speed/fc; # 波长
tvg = EngeePhased.TimeVaryingGain( # 创建VAR系统对象
RangeLoss=2*fspl.(range_gates,lambda), # 范围调整损失
ReferenceLoss=2*fspl.(prop_speed/(prf*2),lambda)); # 范围调整损失
rxpulses = tvg(rxpulses); # WARU应用
# 积分后检测超过检测阈值的峰值
_,range_detect = findpeaks(pulsint(rxpulses,"noncoherent"),
"MinPeakHeight",sqrt(threshold)); # 寻找山峰
range_estimates = round.(range_gates[range_detect]) # 范围估计
这些估计表明存在2000m和3550m范围内的目标。
多普勒频谱
在成功估计目标范围之后，可以估计每个目标的多普勒信息。
多普勒信号估计实际上是频谱估计的过程。 因此，多普勒处理的第一步是计算接收信号的多普勒频谱。
匹配滤波器后的接收信号为矩阵，其列对应于接收脉冲。 与范围估计不同，多普勒处理使用用于脉冲累积时间（慢时间）的数据，即沿着所得矩阵的行。 由于使用了10个脉冲，因此10个样本可用于多普勒处理。 由于每个脉冲取一个样本，多普勒样本的采样频率是脉冲重复率（PRF）。
正如傅立叶理论所预测的，脉冲定位系统可以检测到的最大个位多普勒频移等于其PRF的一半。 这也对应于雷达系统可以检测到的最大个位速度。 此外，脉冲的数量确定多普勒频谱中的分辨率，其确定速度估计的分辨率。
max_speed = dop2speed(prf/2,lambda)/2; # 计算方法
println("最大检测速度：△(round(max_speed;sigdigits=7))m/s");
speed_res = 2*max_speed/num_pulse_int
println("速度分辨率：△(round(speed_res;sigdigits=7))m/s");
如上面的计算所示，在此示例中，最大可检测速度为225m/s，接近（-225）或后退（+225）。 由此产生的多普勒分辨率约为45m/s，这意味着对于多普勒频谱中的分离，两个速度必须彼此相差至少45m/s，为了提高区分不同目标速度的能力，需要更 然而，可用脉冲的数量也受到目标径向速度的限制。 由于多普勒处理被限制在设定的范围内，所以必须在目标从一个分辨率元件切换到另一个分辨率元件之前接收处理中使用的所有脉冲。
由于多普勒样本的数量一般是有限的，所以重置序列以内插所得频谱是常见的。 这不会提高所获得的光谱的分辨率，但它可能改善光谱中峰值位置的估计。
多普勒频谱可以使用周期图获得。
num_range_detected = length(range_estimates);
# 计算2000米范围内的围图
period_out_1 = DSP.periodogram(rxpulses[range_detect[1],:][:],nfft=256,fs=prf)
p1, f1 = period_out_1.power,period_out_1.freq # 提取功率(p1)和频率(f1)值
# 计算3550米范围内的围图
period_out_2 = DSP.periodogram(rxpulses[range_detect[2],:][:],nfft=256,fs=prf)
p2, f2 = period_out_2.power,period_out_2.freq; # 提取功率(p2)和频率(f2)值
然后可以计算频谱中每个频率值对应的速度。 请注意，有必要考虑到圆周运动效果。
speed_vec = dop2speed.(fftshift(f1),lambda)/2; # 将频率网格重新计算到速度网格中
多普勒频移的估计
为了估计与每个目标相关的多普勒频移，我们需要找到每个多普勒频谱中峰值的位置。 在这个例子中，目标存在于两个不同的范围上，因此必须对每个范围重复评估过程。
首先，我们将构建对应于2000米范围的多普勒频谱。
spec_sig1 = fft([rxpulses[range_detect[1],:][:];zeros(256-length(rxpulses[range_detect[1],:][:]))])
spec_power1 = fftshift(real(spec_sig1.*conj.(spec_sig1))./100)
plot(fftshift(fftfreq(256,prf)).*1e-3,pow2db.(spec_power1),lab="",gridalpha=0.15,margin=5*Plots.mm)
xlabel!("频率，千赫", fontfamily="Computer Modern")
ylabel!("SPM，dBW/Hz", fontfamily="Computer Modern")
title!("在2000m范围内反射信号的光谱功率密度", fontfamily="Computer Modern")
请注意，我们只对峰值检测感兴趣，因此频谱值本身是无关紧要的。 从多普勒谱图可以看出，最大峰值以下5dB是一个很好的阈值。 因此，我们使用-5的值作为归一化多普勒频谱的阈值。
spectrum_data = spec_power1./maximum(spec_power1)
_,dop_detect1 = findpeaks(pow2db.(spectrum_data),"MinPeakHeight",-5)
sp1 = speed_vec[dop_detect1]
println("速度估计：$sp1m/s");
结果表明，在2000米范围内有两个目标：一个速度为3.5米/秒，另一个速度为104米/秒。-104米/秒的值可以很容易地与第一个目标相关联，因为第一个目标以100米/秒的径向速度离开，鉴于本例的多普勒分辨率，它非常接近计算值。 3.5m/s的值需要更详细的解释。 由于第三目标在切向方向上运动，因此在径向方向上没有速度分量。 因此，雷达无法检测到第三目标的多普勒频移。 第三目标的真实径向速度因此为0m/s，3.5m/s的估计值非常接近真实值。
请注意，仅使用范围估计值无法区分这两个目标，因为它们的范围值相同。
然后将相同的操作应用于对应于3550米范围的数据。
spec_sig2 = fft([rxpulses[range_detect[2],:][:];zeros(256-length(rxpulses[range_detect[2],:][:]))])
spec_power2 = fftshift(real(spec_sig2.*conj.(spec_sig2))./100)
plot(fftshift(fftfreq(256,prf)).*1e-3,pow2db.(spec_power2),lab="",gridalpha=0.15,margin=5*Plots.mm)
xlabel!("频率，千赫", fontfamily="Computer Modern")
ylabel!("功率，dBW", fontfamily="Computer Modern")
title!("光谱功率密度在3550米范围内", fontfamily="Computer Modern")
spectrum_data2 = spec_power2./maximum(spec_power2)
_,dop_detect2 = findpeaks(pow2db.(spectrum_data2),"MinPeakHeight",-5)
sp2 = speed_vec[dop_detect2]
println("速度等级：＜sp2m/s");
该结果显示了0m/s的估计速度，这对应于目标不在该距离处移动的事实。
结论
此示例显示了使用脉冲雷达系统估计移动目标径向速度的简单方法。 生成接收信号的多普勒频谱，并估计频谱上峰值的位置，这些峰值对应于目标的径向速度。 还讨论了由发射信号的特性和处理算法引起的多普勒处理的局限性。