Engee 文档
Notebook

多普勒估计

该示例显示了一个单稳态脉冲雷达，它可以检测某些范围内移动目标的径向速度。 速度由运动目标引起的多普勒频移决定。 首先，我们将确定目标在给定范围内的存在，然后，使用多普勒处理，我们将估计目标在该范围内的径向速度。

使用的功能和库

In [ ]: 
Pkg.add(["DSP"])

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`
In [ ]: 
using DSP,FFTW

function BasicMonostaticRadarExampleData()
    c = 299_792_458 # cкорость распространения сигнала
    # Создание системных объектов:
    # генератор прямоугольных импульсов
    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 приемника。请读者使用此示例探索雷达系统的设计细节。

In [ ]: 
waveform,antenna,radiator,collector,transmitter,receiver,sensormotion = BasicMonostaticRadarExampleData();

系统建模

目标

多普勒处理使用由运动目标引起的多普勒频移。 现在我们将设置三个目标，指定它们的位置，有效散射面积（EPR）和速度。

In [ ]: 
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]; # ЭПР целей
fc = radiator.OperatingFrequency;
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc);

请注意，第一个目标位于3,550，1,600和1,200m的距离处，并且以60，130和0m/s的速度移动。 第三目标不动。

分布环境

我们还需要为每个目标设置一个分发环境。 由于我们使用的是单稳态雷达，我们将定义一个双向传播模型。:

In [ ]: 
fs = waveform.SampleRate; # частота дискретизации системы
# модель канала распространения
channel = EngeePhased.FreeSpace(
    SampleRate=fs, # частота дискретизации системы
    TwoWayPropagation=true, # модель распространения
    OperatingFrequency=fc); # несущая частота сигнала
信号合成

确定了雷达系统、环境和目标后，我们将模拟从目标反射的接收回波信号。 得到的数据是每列时间快（每个脉冲内的时间）和每行时间慢（脉冲间的时间）的元素矩阵。

为了重现相同的结果，我们需要设置用于在接收器中产生噪声的初始数据。

In [ ]: 
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); # среда распространения
    
    # Отражение сигнала от целей с учетом ЭПР
    tgtsig = target(txsig);
    
    # Прием и предварительного усиление отраженного сигнала
    rxsig = collector(tgtsig,tgtang); # приемная антенна
    rxpulses[:,m] = receiver(rxsig,xor.(txstatus.>0,1)); # предусилитель
end

多普勒估计

范围检测

要估计目标的多普勒频移，首先需要使用范围检测来确定它们的位置。 由于多普勒频移将信号功率分布在I和Q通道上，因此需要使用信号能量进行检测。 这意味着必须使用非相干检测方案。

在上述示例中详细描述了检测过程，因此在这里我们将简单地执行必要的动作以估计目标范围。

In [ ]: 
# Задание вероятности ложной тревоги (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));

# Применение временного автоматического регулирования усиления (ВАРУ)
#   для компенсации потерь, зависящие от дальности.
prop_speed = radiator.PropagationSpeed; # скорость распространения сигнала
range_gates = prop_speed*fast_time_grid/2; # разрешение по дальности
lambda = prop_speed/fc; # длина волны

tvg = EngeePhased.TimeVaryingGain( # создание системного объекта ВАРУ
    RangeLoss=2*fspl.(range_gates,lambda), # потери с учетом дальности
    ReferenceLoss=2*fspl.(prop_speed/(prf*2),lambda)); # потери с учетом дальности

rxpulses = tvg(rxpulses); # применение ВАРУ

# Обнаружения пиков, привышающих порог обнаружения, после интегрирования
_,range_detect = findpeaks(pulsint(rxpulses,"noncoherent"),
    "MinPeakHeight",sqrt(threshold)); # нахождение пиков
range_estimates = round.(range_gates[range_detect]) # оценки дальности
Out[0]:
2-element Vector{Float64}:
 2000.0
 3550.0

这些估计表明存在2000m和3550m范围内的目标。

多普勒频谱

在成功估计目标范围之后，可以估计每个目标的多普勒信息。

多普勒信号估计实际上是频谱估计的过程。 因此，多普勒处理的第一步是计算接收信号的多普勒频谱。

匹配滤波器后的接收信号为矩阵，其列对应于接收脉冲。 与范围估计不同，多普勒处理使用用于脉冲累积时间（慢时间）的数据，即沿着所得矩阵的行。 由于使用了10个脉冲，因此10个样本可用于多普勒处理。 由于每个脉冲取一个样本，多普勒样本的采样频率是脉冲重复率（PRF）。

正如傅立叶理论所预测的，脉冲定位系统可以检测到的最大个位多普勒频移等于其PRF的一半。 这也对应于雷达系统可以检测到的最大个位速度。 此外，脉冲的数量确定多普勒频谱中的分辨率，其确定速度估计的分辨率。

In [ ]: 
max_speed = dop2speed(prf/2,lambda)/2; # расчет 
println("Максимальная скорость обнаружения: $(round(max_speed;sigdigits=7)) м/с");

Максимальная скорость обнаружения: 224.6888 м/с
In [ ]: 
speed_res = 2*max_speed/num_pulse_int
println("Разрешение по скорости: $(round(speed_res;sigdigits=7)) м/с");

Разрешение по скорости: 44.93776 м/с

如上面的计算所示，在此示例中，最大可检测速度为225m/s，接近（-225）或后退（+225）。 由此产生的多普勒分辨率约为45m/s，这意味着对于多普勒频谱中的分离，两个速度必须彼此相差至少45m/s，为了提高区分不同目标速度的能力，需要更 然而，可用脉冲的数量也受到目标径向速度的限制。 由于多普勒处理被限制在设定的范围内，所以在目标从一个分辨率元件切换到另一个分辨率元件之前必须接收处理中使用的所有脉冲。

由于多普勒样本的数量一般是有限的，所以重置序列以内插所得频谱是常见的。 这不会提高所获得的光谱的分辨率，但它可能改善光谱中峰值位置的估计。

多普勒频谱可以使用周期图获得。

In [ ]: 
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)

然后可以计算频谱中每个频率值对应的速度。 请注意，有必要考虑到圆周运动效果。

In [ ]: 
speed_vec = dop2speed.(fftshift(f1),lambda)/2; # пересчет сетки частот в сетку скорости
多普勒频移的估计

为了估计与每个目标相关的多普勒频移，我们需要找到每个多普勒频谱中峰值的位置。 在这个例子中，目标存在于两个不同的范围上，因此必须对每个范围重复评估过程。

首先，我们将构建对应于2000米范围的多普勒频谱。

In [ ]: 
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!("СПМ,дБВт/Гц", fontfamily="Computer Modern")
title!("Спектральная плотность мощности отраженного сигнала по дальности 2000 м", fontfamily="Computer Modern")
Out[0]:

请注意，我们只对检测峰值感兴趣，因此光谱值本身是无关紧要的。 从多普勒谱图可以看出，最大峰值以下5dB是一个很好的阈值。 因此，我们使用-5的值作为归一化多普勒频谱的阈值。

In [ ]: 
spectrum_data = spec_power1./maximum(spec_power1)
_,dop_detect1 = findpeaks(pow2db.(spectrum_data),"MinPeakHeight",-5)
sp1 = speed_vec[dop_detect1]
println("Оценки скорости:  $sp1 м/с");

Оценки скорости:  [-103.56749129975046, 3.510762416940684] м/с

结果表明，在2000米范围内有两个目标：一个速度为3.5米/秒，另一个速度为104米/秒。-104米/秒的值可以很容易地与第一个目标相关联，因为第一个目标以100米/秒的径向速度离开，鉴于本例的多普勒分辨率，它非常接近计算值。 3.5m/s的值需要更详细的解释。 由于第三目标在切向方向上运动，因此在径向方向上没有速度分量。 因此，雷达无法检测到第三目标的多普勒频移。 第三目标的真实径向速度因此为0m/s，3.5m/s的估计值非常接近真实值。

请注意，这两个目标不能单独使用范围估计来区分，因为它们的范围值是相同的。

然后将相同的操作应用于对应于3550米范围的数据。

In [ ]: 
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!("Мощность,дБВт", fontfamily="Computer Modern")
title!("Спектральная плотность мощности по дальности 3550 метров", fontfamily="Computer Modern")
Out[0]:
In [ ]: 
spectrum_data2 = spec_power2./maximum(spec_power2)
_,dop_detect2 = findpeaks(pow2db.(spectrum_data2),"MinPeakHeight",-5)
sp2 = speed_vec[dop_detect2]
println("Оценка скорости: $sp2 м/с");

Оценка скорости: [0.0] м/с

该结果显示了0m/s的估计速度，这对应于目标不在该距离处移动的事实。

结论

此示例显示了使用脉冲雷达系统估计移动目标径向速度的简单方法。 生成接收信号的多普勒频谱，并估计频谱中峰值的位置，这些峰值对应于目标的径向速度。 还讨论了由发射信号的特性和处理算法引起的多普勒处理的局限性。