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;

设置雷达系统

首先，我们定义一个雷达系统。由于本例的重点是多普勒处理，因此我们使用雷达接收机测试信号模拟 示例中构建的雷达系统。我们鼓励读者使用该示例探索雷达系统设计的细节。

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);

请注意，第一个目标的距离分别为 3550 米、1600 米和 1200 米，速度分别为 60 米/秒、130 米/秒和 0 米/秒。不同的是，第一个目标沿径向移动，第二个目标沿切向和径向移动。第三个目标没有移动。

传播环境

我们还需要为每个目标配置传播环境。由于我们使用的是单静态雷达，因此将设置双向传播模型：

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

这些估算结果表明，目标距离在 2000 米和 3550 米之间。

多普勒频谱

一旦成功估算出目标距离，就可以估算出每个目标的多普勒信息。

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

经过匹配滤波器后的接收信号是一个矩阵，其列与接收到的脉冲相对应。与测距估算不同，多普勒处理使用的是脉冲累积时间（慢速时间）内的数据，即所得矩阵的各行数据。由于使用了 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 м/с

如上计算所示，在本例中，可探测到的最大速度为 225 米/秒，要么正在接近（-225），要么正在远离（+225）。由此得出的多普勒分辨率约为 45 m/s，这意味着两个速度必须相差至少 45 m/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]:

请注意，我们只对峰值检测感兴趣，因此频谱值本身并不重要。从多普勒频谱图中我们可以看出，低于最大峰值 5 dB 是一个很好的阈值。因此，我们使用 -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 m/s 的值很容易与第一个目标联系起来，因为第一个目标的径向速度为 100 m/s，考虑到本例的多普勒分辨率，这个值与计算值非常接近。至于 3.5 米/秒的数值，则需要更详细的解释。由于第三个目标沿切线方向移动，因此径向上没有速度分量。因此，雷达无法探测到第三个目标的多普勒位移。因此，第三个目标的真实径向速度为 0 m/s，而 3.5 m/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] м/с

结果显示估计速度为 0 米/秒，这与目标在此范围内没有移动这一事实相符。

结论

本示例展示了使用脉冲雷达系统估算移动目标径向速度的简单方法。我们生成了接收信号的多普勒频谱，并估算了频谱上与目标径向速度相对应的峰值位置。我们还讨论了辐射信号的特性和处理算法对多普勒处理造成的限制。