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(["功率，dBW"],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的非波动目标 距离雷达5000米，范围分辨率为50米。 期望的性能指标是0.9检测概率(Pd)和更少的虚警概率(Pfa)。 . 由于相干检测需要相位信息，因此计算昂贵，我们使用非相干检测方案。 此外，此示例假定环境处于可用空间中。

In [ ]: 
pd = 0.9; # 正确检测的概率
pfa = 1e-6; # 虚警的概率
max_range = 5000; # [m]最大范围
range_res = 50; # [m]所需范围分辨率
tgt_rcs = 1;  # [m^2]EPR目标

单稳态雷达系统的设计

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

发电机

在这个例子中，我们选择了一个矩形波形。 所需的范围分辨率决定了波形的带宽，在矩形波形的情况下，它决定了脉冲宽度。

脉冲波形的另一个重要参数是脉冲重复率(PRF)。 PRF由最大个位数范围确定。

In [ ]: 
prop_speed = physconst("lightspeed"); # [m/s]传播速度
pulse_bw = prop_speed/(2*range_res);  # [Hz]频域脉冲带
pulse_width = 1/pulse_bw; # [c]脉冲持续时间
prf = prop_speed/(2*max_range); # [Hz]脉冲重复率
fs = 2*pulse_bw; # [Hz]系统的采样率
waveform = EngeePhased.RectangularWaveform(PulseWidth=1/pulse_bw,PRF=prf,SampleRate=fs); # 系统对象

请注意，我们已将采样率设置为带宽的两倍。

接收机噪声特性

我们假设接收器中唯一存在的噪声是热噪声，因此在该仿真中没有干扰。 热噪声的功率取决于接收器的带宽。 接收器的噪声带宽被设置为等于波形的带宽。 在实际系统中通常是这种情况。 我们还假设接收器的增益为20dB，噪声系数为0dB。

In [ ]: 
noise_bw = pulse_bw; # [Hz]接收机噪声带
receiver = EngeePhased.ReceiverPreamp(
    Gain=20, # [dB]增益
    NoiseFigure=0,# [dB]噪音系数
    SampleRate=fs, # [Hz]采样率
    EnableInputPort=true # 用于同步发射器和接收器的端口
);

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

发射器

发射机最重要的参数是峰值功率。 所需的峰值功率与许多因素有关，包括最大个位数范围、接收机中所需的SNR和脉冲宽度。 在这些因素中，接收机中所需的SNR由Pd和Pfa的设计目的以及在接收机中实现的检测方案决定。

Pd、Pfa和SNR之间的关系可以最好地用接收机的操作性能曲线(ROC)来表示。 我们可以使用以下命令构造一条曲线，其中Pd是各种Snr的Pfa的函数:

In [ ]: 
snr_db = [-Inf, 0.0, 3, 10, 13]; # Snr值（dB）
rocsnr(snr_db,SignalType="NonfluctuatingNoncoherent") # 建筑物检测特征
Out[0]:

ROC曲线表明，为了满足设计目标，Pfa=1e-6和Pd=0.9;接收信号的SNR必须超过13dB。 这是一个相当高的要求，不是很实用。 为了使雷达系统更可行，我们可以使用脉冲累积技术来降低所需的SNR。 如果我们选择对10个脉冲进行积分，那么曲线可以如下构造:

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

我们可以看到所需的功率已降至约5dB。 SNR的进一步降低可以通过积分更多脉冲来实现，但是可用于积分的脉冲数量通常由于目标移动或环境异质性而受到限制。

在上述方法中，SNR值是从曲线读取的，但通常希望仅计算所需值。 对于非相干检测方案，计算所需的SNR在理论上是相当困难的。 幸运的是，有很好的近似值，例如Albersham方程。 使用Albersham方程，所需的SNR值可以导出为

In [ ]: 
snr_min = albersheim(pd, pfa,N= num_pulse_int) # 通过求解Albrshem方程计算SNR
print("最小SNR：△(round(snr_min;sigdigits=5))dB")

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

一旦您在接收器处获得了所需的SNR值，就可以使用雷达方程计算发射器处的峰值功率。 这里我们假设发射机的增益为20dB。

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

In [ ]: 
tx_gain = 20; # [dB]发射机增益
fc = 10e9; # [Hz]信号的载波频率
lambda = prop_speed/fc; # [m]信号波长
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)kW")

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

请注意，产生的功率约为5kW，这是一个可接受的值。 为了比较，如果我们没有使用脉冲积分方法，峰值功率将是33千瓦，这是一个相当大的值。

收到所有这些信息后，我们可以设置发射器。

In [ ]: 
transmitter = EngeePhased.Transmitter(
    Gain=tx_gain, # [dB]增益
    PeakPower=peak_power, # [W]峰值功率
    InUseOutputPort=true # 启用定义的输出端口
                         # 发射器的状态（1-on，0-off）
);

由于本示例模拟的是单稳态雷达系统，因此将InUseOutputPort端口设置为true以输出发射器的状态。 该状态信号可用于接通接收器。

发送和接收装置

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

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

我们假设天线是静止的。

In [ ]: 
# 雷达天线元件
antenna = EngeePhased.IsotropicAntennaElement(
    FrequencyRange=[5e9 15e9]); # [Hz]辐射天线的频率范围
# 雷达运动模型
sensormotion = EngeePhased.Platform(
    InitialPosition=[0; 0; 0], # [m]天线的初始位置
    Velocity=[0; 0; 0]); # [m/s]天线速度

使用天线和工作频率，我们设置辐射器（Radiator）和接收天线（Collector）。

In [ ]: 
radiator = EngeePhased.Radiator(
    Sensor=antenna,             # 设置天线元件
    OperatingFrequency=fc);   # [Hz]设定载波频率

collector = EngeePhased.Collector(
    Sensor=antenna,             # 设置天线元件
    OperatingFrequency=fc);     # [Hz]设定载波频率

这样就完成了雷达系统的配置。 在以下部分中，我们将定义建模所必需的其他对象，例如目标和环境。 然后我们将模拟信号的返回，并根据模拟的信号执行范围确定。

系统建模

目标

为了测试我们的雷达探测目标的能力，我们必须首先识别目标。 让我们假设空间中有3个静止的，非波动的目标。 它们的位置和雷达横截面如下所示。

In [ ]: 
tgtpos = [[2024.66;0;0] [3518.63;0;0] [3845.04;0;0]]; # [m]3个目标初始位置的分量[[px][py][pz]]
tgtvel = [[0;0;0] [0;0;0] [0;0;0]]; # 3个目标的速度矢量的[m/s]分量[[vx][vy][vz]]
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # 目标移动模型

tgtrcs = [1.6 2.2 1.05]; # M^2EPR目标
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc); # 目标系统对象
分布环境

为了模拟信号，我们还需要确定雷达系统和每个目标之间的传播信道。

In [ ]: 
channel = EngeePhased.FreeSpace(
    SampleRate=fs, # [Hz]采样率
    TwoWayPropagation=true, # 双向传播
    OperatingFrequency=fc); # [Hz]载波频率

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

信号合成

现在我们已经准备好模拟整个系统。

合成信号是每列时间快（每个脉冲内部的时间）和每行时间慢（脉冲之间的时间）的数据矩阵。 为了可视化信号，定义快速时间网格和慢时间网格非常有用。

In [ ]: 
fast_time_grid = (0:1/fs:1/prf-1/fs); # [c]快速时间网格
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); # 分配环境的核算
    
    # 从目标的信号的反射，考虑到EPR
    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);  # [W]噪声功率
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]:

这些附图中的阈值仅用于演示目的。 请注意，来自第二和第三个目标的信号比来自第一个目标的信号弱得多，因为它们距离雷达更远。 因此，接收信号的功率取决于范围，并且阈值对位于不同距离的目标是不公平的。

一致的过滤

匹配滤波器提供增强的处理，这增加了检测阈值。 它将接收到的信号转换为发送振荡器类型的本地、时间反转和共轭副本。 因此，在创建匹配滤波器时，我们必须指定探测信号的类型。 在执行脉冲积分、阈值检测等之前，接收到的脉冲首先通过匹配滤波器以改善SNR。

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]:

在匹配滤波器级之后，SNR提高。 然而，由于接收信号的功率取决于范围，来自附近目标的信号仍然比来自位于更远距离的目标强得多。 因此，如上图所示，来自附近bean的噪声很有可能超过阈值并使目标更远。 为了确保阈值对于可检测范围内的所有目标都有效，我们可以使用时变增益因子来补偿接收回波信号中与范围相关的损失。

为了补偿与范围相关的损耗，我们首先计算对应于每个信号样本的选通脉冲，然后计算对应于每个选通脉冲的自由空间损耗。 接收到这些信息后，我们对接收到的脉冲施加时变增益，使返回看起来来自相同的参考范围（最大可检测范围）。

In [ ]: 
range_gates = prop_speed*fast_time_grid/2; # [m]范围分辨率
tvg = EngeePhased.TimeVaryingGain( # 创建VAR系统对象
    RangeLoss=2*fspl.(range_gates,lambda), # [dB]范围调整损耗
    ReferenceLoss=2*fspl.(max_range,lambda)); # [dB]最大范围内的相对损耗
rxpulses = tvg(rxpulses); # 运算后的输出信号为WARU

现在让我们在范围归一化后绘制相同的两个脉冲。

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

随时间改变增益导致噪声电平的增加。 但是，目标的返回现在独立于范围。 现在可以使用恒定阈值在整个可检测范围内进行检测。

注意，在该阶段，阈值高于每个脉冲中包含的最大功率电平。 因此，在此阶段无法检测到任何内容。 我们需要对脉冲进行积分，以便来自目标的返回回波的功率超过阈值，并且噪声水平保持在条形之下。 这是可以预料的，因为正是脉冲积分使我们能够使用一系列低功率脉冲。

不连贯的积累

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

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)m")
print("目标范围估计：$(range_estimates)m")

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

请注意，这些范围估计值仅在雷达系统可以实现的范围分辨率（50m）下准确。

结论

在这个例子中，我们设计了一个基于一组设定目标的雷达系统。 基于这些目标，计算了雷达系统的许多设计参数。 该示例还显示了如何使用开发的雷达执行范围检测任务。 在本例中，雷达使用矩形波形。