Engee 文档
Notebook

雷达接收器测试信号建模

本例说明如何模拟单静态脉冲雷达的接收信号，以估计目标的距离。单静态雷达的发射器位于接收器附近。发射器产生的脉冲击中目标后会产生回声，并被接收器接收。通过测量回波随时间变化的位置，可以估算出目标的距离。

本案例研究的重点是设计一个能达到给定规格的脉冲雷达系统。它描述了将探测概率和测距分辨率等设计规格转化为发射机功率和脉冲持续时间等雷达系统参数的步骤。还对环境和目标进行建模，以获得反射信号。最后，对接收到的信号采用信号处理技术，以确定目标的距离。

使用的函数

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

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

# Функция для построения выходного отклика радара и порогового значения
function RadarPulsePlot(ReceivePulse,threshold,t,tp,num_plot::Int64)
    thresh = sqrt(threshold)*ones(length(t),num_plot)
    ReceivePulse =  Float64.(pow2db.((abs.(ReceivePulse)).^2))
    ReceivePulse[ReceivePulse.==0] .= eps(0.)
    thresh .=  pow2db.(thresh.^2)
    t = repeat(t,1,num_plot) .+ tp[1:num_plot]'
    text_title = reshape(["Импульс № $(i)" for i in 1:num_plot],1,num_plot)
    xlabel_text = [reshape(repeat([""],1,num_plot-1),1,num_plot-1) "Время, с"]
    plot(t,ReceivePulse[:,1:num_plot],layout=(num_plot,1),label="",title = text_title,
        fontfamily="Computer Modern", titlefontsize=10, guidefontsize=10, tickfontsize=10,gridalpha=0.3,lw=1.2)
    plot!(t,thresh[:,1:num_plot],layout=(num_plot,1),xlabel=xlabel_text,lw=1.2,
        ylabel = repeat(["Мощность, дБВт"],1,num_plot),label="",ls=:dash) # 
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;

设备规格

该脉冲雷达系统的设计目标是探测雷达截面（RCS）至少为 1$м^2$ 的非波动目标，探测距离最远可达 5000 米，测距分辨率为 50 米。期望的性能指标是探测概率 (Pd) 为 0.9，误报概率 (Pfa) 小于$10^{-6}$ 。由于相干探测需要相位信息，计算成本较高，因此我们采用非相干探测方案。此外，本示例假定介质位于自由空间。

In [ ]: 
pd = 0.9; # Вероятность верного обнаружения
pfa = 1e-6; # Вероятность ложной тревоги
max_range = 5000; # [м] Максимальная дальность
range_res = 50; # [м] Требуемое разрешение по дальности
tgt_rcs = 1;  # [м^2] ЭПР цели

设计单静态雷达系统

我们需要确定雷达系统的几个特征，如波形、接收器、发射器和用于发射和接收信号的天线。

发电机

在本例中，我们选择了矩形波形。所需的量程分辨率决定了波形的带宽，而矩形波形的带宽决定了脉冲宽度。

脉冲波形的另一个重要参数是脉冲重复频率（PRF）。PRF 由最大个位数量程决定。

In [ ]: 
prop_speed = physconst("lightspeed"); # [м/с] Скорость распространения
pulse_bw = prop_speed/(2*range_res);  # [Гц] Полоса импульса в частотной области
pulse_width = 1/pulse_bw; # [c] Длительность импульса
prf = prop_speed/(2*max_range); # [Гц] Частота следования импульсов
fs = 2*pulse_bw; # [Гц] Частота дискретизации системы
waveform = EngeePhased.RectangularWaveform(PulseWidth=1/pulse_bw,PRF=prf,SampleRate=fs); # системый объект

请注意，我们设置的采样率等于带宽的两倍。

接收器的噪声特性

我们假定接收机中存在的唯一噪声是热噪声，因此本次模拟不存在干扰。热噪声的功率取决于接收机带宽。接收器噪声带宽等于波形带宽。实际系统中通常就是这种情况。我们还假设接收器的增益为 20 dB，噪声系数为 0 dB。

In [ ]: 
noise_bw = pulse_bw; # [Гц] шумовая полоса приемника
receiver = EngeePhased.ReceiverPreamp(
    Gain=20, # [дБ] усиление 
    NoiseFigure=0,# [дБ] фактор шума
    SampleRate=fs, # [Гц] частота дискретизации
    EnableInputPort=true # порт для синхронизации передатчика и приемника
);

请注意，由于我们模拟的是单静态雷达，因此在发射机关闭之前，接收机无法打开。因此，我们设置标志EnableInputPort=true ，以便同步信号从发射机传输到接收机。

发射机

最重要的发射机参数是峰值功率。所需的峰值功率与许多因素有关，包括最大个位数量程、接收器所需的信噪比和脉冲宽度。其中，接收器所需的信噪比由设计目标 Pd 和 Pfa 以及接收器中实施的检测方案决定。

Pd、Pfa 和信噪比之间的关系可以用接收机工作特性曲线（ROC）来表示。我们可以使用以下命令绘制曲线，其中 Pd 是不同信噪比下 Pfa 的函数：

In [ ]: 
snr_db = [-Inf, 0.0, 3, 10, 13]; # значения ОСШ в дБ
rocsnr(snr_db,SignalType="NonfluctuatingNoncoherent") # построение характеристик обнаружения
Out[0]:

ROC 曲线显示，要达到 Pfa = 1e-6 和 Pd = 0.9 的设计目标，接收信号的信噪比必须超过 13 dB。这是一个相当高的要求，而且不太实用。为了使无线电定位系统更加可行，我们可以使用脉冲累积技术来降低所需的信噪比。如果我们选择对 10 个脉冲进行积分，可以绘制出如下曲线：

In [ ]: 
num_pulse_int = 10;
rocsnr([0 3 5],SignalType="NonfluctuatingNoncoherent",NumPulses=num_pulse_int)
Out[0]:

我们可以看到，所需功率已降至约 5 dB。通过整合更多脉冲可以进一步降低信噪比，但由于目标移动或环境不均匀，可供整合的脉冲数量通常有限。

在上述方法中，信噪比值是从曲线中读取的，但通常只需要计算所需值即可。对于非相干检测方案，计算所需的信噪比在理论上相当困难。幸运的是，有一些很好的近似值，如阿尔伯沙姆方程。利用阿尔伯沙方程，所需的信噪比值可推导为

In [ ]: 
snr_min = albersheim(pd, pfa,N= num_pulse_int) # расчет ОСШ путем решения уравнения Альбршема
print("Минимальное ОСШ: $(round(snr_min;sigdigits=5)) дБ")

Минимальное ОСШ: 4.9904 дБ

得到接收器所需的信噪比值后，我们就可以利用雷达方程计算发射器的峰值功率。这里我们假设发射机增益为 20 dB。

要利用雷达方程计算峰值功率，我们还需要知道传播信号的波长，这与系统的工作频率有关。这里我们将工作频率设为 10 GHz。

In [ ]: 
tx_gain = 20; # [дБ] Усиление передатчика
fc = 10e9; # [Гц] несущая частота сигнала
lambda = prop_speed/fc; # [м] длина волны сигнала
peak_power = ((4*pi)^3*noisepow1(1/pulse_width)*max_range^4*
    DSP.db2pow(snr_min))/(DSP.db2pow(2*tx_gain)*tgt_rcs*lambda^2) # пиковая мощность
print("Пиковая мощность: $(peak_power*1e3) кВт")

Пиковая мощность: 5.22647386447291e6 кВт

请注意，得出的功率约为 5 千瓦，这是一个可以接受的值。相比之下，如果我们不使用脉冲积分法，峰值功率将达到 33 千瓦，这是一个相当大的数值。

获得所有这些信息后，我们就可以配置发射机了。

In [ ]: 
transmitter = EngeePhased.Transmitter(
    Gain=tx_gain, # [дБ] усиление
    PeakPower=peak_power, # [Вт] пиковая мощность
    InUseOutputPort=true # включение выходного порта, который определяет 
                         # состояние передатчика (1 - вкл, 0 - выкл)
);

由于本例模拟的是单静态雷达系统，因此 InUseOutputPort 设置为 true，以输出发射机状态。该状态信号可用于打开接收机。

发射机和接收机

在雷达系统中，信号以电磁波的形式传播。因此，信号必须由雷达系统中使用的天线辐射和收集。这就是辐射天线和 "收集 "天线发挥作用的地方。

在单静态雷达系统中，发射器和收集器使用同一天线，因此我们首先要确定天线。为了简化设计，我们选择了各向同性天线。请注意，天线必须在系统工作频率（10 GHz）下工作，因此我们将天线频率范围设定为 5-15 GHz。

我们假设天线是静止的。

In [ ]: 
# антенный элемент радара
antenna = EngeePhased.IsotropicAntennaElement(
    FrequencyRange=[5e9 15e9]); # [Гц] диапазон частот антенны излучения
# модель движения радара
sensormotion = EngeePhased.Platform(
    InitialPosition=[0; 0; 0], # [м] начальное положение антенны
    Velocity=[0; 0; 0]); # [м/с] скорость антенны

根据天线和工作频率，我们定义了辐射器（Radiator）和接收天线（Collector）。

In [ ]: 
radiator = EngeePhased.Radiator(
    Sensor=antenna,             # задание антенного элемента 
    OperatingFrequency=fc);   # [Гц] задание несущей частоты

collector = EngeePhased.Collector(
    Sensor=antenna,             # задание  антенного элемента 
    OperatingFrequency=fc);     # [Гц] задание несущей частоты

至此，雷达系统配置完成。在接下来的章节中，我们将定义模拟所需的目标和环境等其他对象。然后，我们将模拟信号返回，并根据模拟信号进行测距。

系统建模

目标

要测试雷达探测目标的能力，我们必须首先定义目标。假设太空中有 3 个静止、无波动的目标。它们的位置和雷达截面如下。

In [ ]: 
tgtpos = [[2024.66;0;0] [3518.63;0;0] [3845.04;0;0]]; # [м] компоненты начального положения 3 целей [[px] [py] [pz]]
tgtvel = [[0;0;0] [0;0;0] [0;0;0]]; # [м/с] компоненты вектора скорости 3 целей [[vx] [vy] [vz]]
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # модель движения целей

tgtrcs = [1.6 2.2 1.05]; # м^2 ЭПР целей
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc); # Системный объект целей

传播环境

为了建立信号模型，我们还需要定义雷达系统与每个目标之间的传播信道。

In [ ]: 
channel = EngeePhased.FreeSpace(
    SampleRate=fs, # [Гц] частота дискретизации
    TwoWayPropagation=true, # двунаправленное распространение
    OperatingFrequency=fc); # [Гц] несущая частота

由于本例使用的是单静态雷达系统，因此通道的设置是为了模拟两个方向的传播延迟。

信号合成

现在，我们可以为整个系统建模了。

合成后的信号是一个数据矩阵，每列为快速时间（每个脉冲内的时间），每行为慢速时间（脉冲之间的时间）。为了使信号可视化，我们需要定义快速时间网格和慢速时间网格。

In [ ]: 
fast_time_grid = (0:1/fs:1/prf-1/fs); # [с] сетка быстрого времени
slow_time_grid = (0:num_pulse_int-1)./prf; # [c] сетка медленного времени

下一个周期模拟接收信号的 10 个脉冲。

设置接收机中噪声发生器的初始状态，以便重现相同的结果。

In [ ]: 
release!(receiver) # доступ к обновлению параметров объекта
receiver.SeedSource = "Property";  # изменение тип модели генерации шума
receiver.Seed = 2007; # задание начального состояния генератора тепловых шумов

# Выделение памяти для результирующих данных
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);

    # Расчет углов по дальности от радара до целей  
    global 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

范围确定

检测阈值

探测器将信号强度与预先设定的阈值进行比较。在雷达应用中，阈值的选择通常是使 Pfa 低于某一水平。在这种情况下，我们假设噪声是白高斯噪声，检测是不连贯的。由于我们也使用 10 个脉冲进行脉冲积分，阈值信号功率定义如下

In [ ]: 
npower = noisepow1(noise_bw,receiver.NoiseFigure,
    receiver.ReferenceTemperature);  # [Вт] мощность шума
threshold = npower * DSP.db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # порог обнаружения

我们将前两个接收到的脉冲与阈值对比绘制成图

In [ ]: 
num_pulse_plot = 2; # количество графиков
RadarPulsePlot(rxpulses,threshold,
    fast_time_grid,slow_time_grid,num_pulse_plot)
Out[0]:

这些图中的阈值仅用于演示目的。请注意，第二和第三个目标的信号比第一个目标弱得多，因为它们距离雷达更远。因此，接收到的信号强度与距离有关，阈值对不同距离的目标是不公平的。

匹配过滤

匹配滤波器可提供处理增益，提高检测阈值。它将接收到的信号转换成发射振荡器类型的本地、时间反转和连续副本。因此，在创建匹配滤波器时，我们必须指定探测信号的类型。接收到的脉冲首先要通过匹配滤波器，以提高信噪比，然后再进行脉冲积分、阈值检测等。

In [ ]: 
matchingcoeff = getMatchedFilter(waveform); # расчет коэффициентов согласованного 
                                            # фильтра на основе изулченного сигнала          
matchedfilter = EngeePhased.MatchedFilter(
    Coefficients=matchingcoeff, # коэффициенты фильтра
    GainOutputPort=true); # выходной порт усиления
rxpulses, mfgain = matchedfilter(rxpulses); # отклик после согласованной фильтрации

匹配滤波器会引入内部延迟，因此峰值位置不再与真正的目标位置相匹配。为了补偿这一延迟，我们将匹配滤波器的输出前移，并在末端添加零。需要注意的是，在实际系统中，由于数据是连续累积的，因此阵列不会有终点。

In [ ]: 
matchingdelay = size(matchingcoeff,1)-1; # длина коэффициентов согласованного фильтра
rxpulses = buffer(rxpulses[matchingdelay+1:end],200); # size(rxpulses,1) ->200

这样，阈值就会因匹配滤波器的处理增益而增加。

In [ ]: 
threshold = threshold * db2pow(mfgain)

下图显示的是通过匹配滤波器后的两个脉冲。

In [ ]: 
RadarPulsePlot(rxpulses,threshold,
    fast_time_grid,slow_time_grid,num_pulse_plot)
Out[0]:

经过匹配滤波器阶段后，信噪比有所改善。然而，由于接收到的信号功率与距离有关，来自近距离目标的信号仍然比来自远距离目标的信号要强得多。因此，如上图所示，来自近处垃圾箱的噪声有很大几率超过阈值，并掩盖远处的目标。为了确保阈值对可探测范围内的所有目标都有效，我们可以使用时变增益来补偿接收到的回声中与范围相关的损失。

为了补偿随距离变化的损耗，我们首先计算每个信号样本对应的选通脉冲，然后计算每个选通脉冲对应的自由空间损耗。获得这些信息后，我们对接收到的脉冲应用时变增益，使回波与来自同一参考范围（最大可探测范围）的回波相同。

In [ ]: 
range_gates = prop_speed*fast_time_grid/2; # [м] разрешение по дальности
tvg = EngeePhased.TimeVaryingGain( # создание системного объекта ВАРУ
    RangeLoss=2*fspl.(range_gates,lambda), # [дБ] потери с учетом дальности
    ReferenceLoss=2*fspl.(max_range,lambda)); # [дБ] относительные потери при максимальной дальности
rxpulses = tvg(rxpulses); # выходной сигнал после работы ВАРУ

现在绘制范围归一化后的两个脉冲图

In [ ]: 
RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,num_pulse_plot)
Out[0]:

增益随时间的变化导致噪声增加。不过，目标返回现在与量程无关。现在可以在整个可探测范围内使用恒定阈值进行探测。

请注意，在此阶段，阈值高于每个脉冲所含的最大功率水平。因此，在这个阶段还无法检测到任何东西。我们需要进行脉冲积分，使目标返回的回波功率超过阈值，而噪声水平保持在条形以下。这是意料之中的，因为正是脉冲积分使我们能够使用较低功率的脉冲序列。

不连贯累积

通过对接收到的脉冲进行非相干积分（视频积分），我们可以进一步提高信噪比。

In [ ]: 
rxpulses = pulsint(rxpulses,"noncoherent"); # некогерентное накопление
RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,1)
Out[0]:

在视频积分步骤之后，数据就可以进行最后的检测步骤了。从图中我们可以看到，三个目标回波都超过了阈值，因此可以被检测到。

范围检测

最后，对集成脉冲进行阈值检测。检测电路可识别峰值，然后将其位置转换为目标范围。

In [ ]: 
_,range_detect = findpeaks(rxpulses,"MinPeakHeight",sqrt(threshold)); # обнаружение пиков

真实的目标范围和检测到的目标范围如下所示：

In [ ]: 
true_range = round.(tgtrng;sigdigits=4)
range_estimates = round.(range_gates[range_detect];sigdigits=4)
println("Истинная дальность целей: $(true_range) м")
print("Оценка дальности целей: $(range_estimates) м")

Истинная дальность целей: [2025.0 3519.0 3845.0] м
Оценка дальности целей: [2025.0, 3525.0, 3850.0] м

请注意，这些测距估计值的精确度只能达到雷达系统所能达到的测距分辨率（50 米）。

结论

在本案例研究中，我们根据一组指定目标设计了一个雷达系统。根据这些目标，计算了雷达系统的许多设计参数。本例还展示了如何使用所设计的雷达执行测距任务。在本例中，雷达使用了矩形波形。