检测曲线（ROC）第2部分：蒙特卡罗模拟

本例演示了如何使用蒙特卡罗仿真构建雷达系统的接收机检测曲线（ROC曲线）。 ROC确定系统检测目标的能力，同时在目标缺失时丢弃大的误报值（误报）。检测系统通过将接收信号值与设定的阈值进行比较来确定目标的存在或不存在。

**正确检测概率(pd)**是当目标实际存在时信号的瞬时值大于阈值的概率。

**虚警概率(pfa)**是目标缺失时信号值大于阈值的概率。在这种情况下，信号是由噪声引起的，其属性取决于噪声统计。蒙特卡罗模拟在存在和不存在目标的情况下产生非常大量的雷达信号。

ROC曲线是Pd对Pfa的依赖性的图。 ROC曲线的形状取决于接收信号的信噪比(SNR)。如果已知接收信号的SNR，则ROC曲线通过分析Pd和Pfa的值来显示系统的可靠性。如果设置Pd和Pfa，则可以确定满足此要求所需的功率。

您可以使用rocsnr函数计算理论ROC曲线。该示例显示了作为单天线雷达系统的蒙特卡罗模拟的结果获得的ROC曲线，并将该曲线与理论曲线进行了比较。

雷达的战术和技术特点

设置所需的Pd=0.9和Pfa= ，我们将设置4000米的最大雷达范围和50米的范围分辨率。将实际目标范围设置为3000米。我们还将目标的有效散射面积（ESR）设置为1.5 和10GHz的工作频率。

Pkg.add(["DSP"])
using DSP

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

using DSP

c = 299_792_458; # скорость распространения сигнала в пространстве, м/с
pd = 0.9; # вероятность правильного обнаружения
pfa = 1e-6; # вероятность ложной тревоги
max_range = 4000; # максимальная дальность, м
target_range = 3000.0; # дальность до цели, м
range_res = 50; # разрешение по дальности,м
tgt_rcs = 1.5; # ЭПР цели, м^2
fc = 10e9; # несущая частота, Гц
lambda = c/fc; # длина волны, м

任何计算Pfa和pd的仿真都需要处理多个信号。为了不占用大量内存，我们将只处理部分脉冲。我们将设置处理的脉冲数45_000和每个部分的大小10_000。

Npulse = 45_000; # количество импульсов
Npulsebuffsize = 10_000; # количество отсчетов в импульсе

探测信号类型及其参数的选择

让我们为LFM信号发生器创建一个系统对象。要做到这一点，你需要计算:
*脉冲带使用范围分辨率;
脉冲重复率（PPI），基于雷达的最大范围。
采样率（fs）**将带宽设置为两倍，因为信号是宽带的;
*脉冲持续时间，基于带宽。

pulse_bw = c/(2*range_res); # полоса импульса
prf = c/(2*max_range); # частота следования импульсов
fs = 2*pulse_bw; # частота дискретизации
pulse_duration = 10/pulse_bw; # длительность импульса

# Формирование генератора ЛЧМ-сигнала
waveform = EngeePhased.LinearFMWaveform(
    PulseWidth=pulse_duration,
    SampleRate=fs,
    SweepBandwidth=pulse_bw,
    PRF=prf
);

为了实现一定的Pd和Pfas，要求在目标反射信号后向接收器提供足够的信号功率。让我们使用Albersham方程计算实现给定虚警概率和检测概率所需的最小SNR值。

snr_min = albersheim(pd,pfa); # минимальное значение ОСШ
println("ОСШ = $(round(snr_min;sigdigits=4)) дБ")

ОСШ = 13.11 дБ

要达到SNR=13dB，必须向目标传输足够的功率。使用雷达方程估计在3000米范围内的目标实现给定SNR所需的峰值发射功率。接收信号还取决于EPR，其被假定为对应于非波动模型（Swerling0)。我们将雷达设置为具有相同的接收和发送增益系数，等于20dB。

txrx_gain = 20; # усиление передатчика и приемника

# Расчет дальности с помощью основного уравнения радилокации
peak_power = ((4*pi)^3*noisepow(1/pulse_duration)*target_range^4*
    DSP.db2pow(snr_min))/(DSP.db2pow(2*txrx_gain)*tgt_rcs*lambda^2) # 
println("peak_power = $(round(peak_power;sigdigits=5)) Вт")

peak_power = 293.18 Вт

设置发射机的系统对象的参数

让我们创建系统对象来定义仿真的发射部分：天线运动模型、天线几何形状、发射器和辐射器。

# модель движения радара
antennaplatform = EngeePhased.Platform(
    InitialPosition=[0; 0; 0], # начальное положение радара
    Velocity=[0; 0; 0] # компоненты вектора скорости радара
); 
# антенный элемент
antenna = EngeePhased.IsotropicAntennaElement(
    FrequencyRange=[5e9 15e9] # диапазон частот
);
# передатчик
transmitter = EngeePhased.Transmitter(
    Gain=txrx_gain, # усиление
    PeakPower=peak_power, # пиковая мощность
    InUseOutputPort=true # состояние радара
);
radiator = EngeePhased.Radiator(
    Sensor=antenna, # антенный элемент
    OperatingFrequency=fc # несущая частота
);

设置目标对象

创建与具有指定EPR的真实反射目标相对应的Target对象。来自该目标的反射将模拟真实的反射信号。要计算假警报，请创建具有零EPR的第二个目标。除噪声外，来自该目标的反射为零。

target = [EngeePhased.RadarTarget() EngeePhased.RadarTarget()] # вектор целей
targetplatform = [ EngeePhased.Platform()  EngeePhased.Platform()] # вектор моделей движения целей
# задание параметров целей
target[1] = EngeePhased.RadarTarget(
    MeanRCS=tgt_rcs, # ЭПР 
    OperatingFrequency=fc # несущая частота
);
target[2] = EngeePhased.RadarTarget(
    MeanRCS=1e-16,
    OperatingFrequency=fc
)
# задание параметров движения целей
targetplatform[1] = EngeePhased.Platform(
    InitialPosition=[target_range; 0; 0] # начальное положение
)

targetplatform[2] = EngeePhased.Platform(
    InitialPosition=[target_range; 0; 0] # начальное положение
);

设置分发环境对象的参数

让我们设置信号传播的条件-从雷达到目标和返回的传播通道。

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

设置发射机的系统对象的参数

让我们通过将NoiseMethod属性设置为**"噪声温度"，将噪声温度属性设置为ReferenceTemperature**设置为290K来产生噪声。

# Геометрия приемной антенны
collector = EngeePhased.Collector(
    Sensor=antenna, # антенна
    OperatingFrequency=fc # Несущая частота
);
# формирование передатчика
receiver = EngeePhased.ReceiverPreamp(
    Gain=txrx_gain, # усиление
    NoiseMethod="Noise temperature", # источник шумов - тепловые шумы
    ReferenceTemperature=290.0, # шумовая температура
    NoiseFigure=0, # фактор шума
    SampleRate=fs, # частота дискретизации
    EnableInputPort=true # синхронизация с передатчиком
);
# Задание начального значения источника шума
receiver.SeedSource = "Property";
receiver.Seed = 2010;

设置快速时间网格

快速时间网格是脉冲持续时间内的一组时间样本。每个计数对应于一个范围分辨率元素。

fast_time_grid = unigrid(0,1/fs,1/prf,"[)"); # сетка быстрого времени
rangebins = c*fast_time_grid/2; # элементы разрешения по дальности

探测信号的生成

我们将生成您要传输的探测信号的脉冲。

wavfrm = waveform(); # формирования ЛЧМ-сигнала

接下来，我们将创建一个包括发射天线放大的发射信号。

sigtrans,tx_status = transmitter(wavfrm);

从波形对象创建一致的滤波器系数。然后创建一个一致的过滤器系统对象。

MFCoeff = getMatchedFilter(waveform); # коэффициенты фильтра
matchingdelay = size(MFCoeff,1) - 1; # заддержка на согласованную фильтрацию
mfilter = EngeePhased.MatchedFilter( # согласованный фильтр
    Coefficients=MFCoeff,
    GainOutputPort=false 
);

计算到目标的范围

计算到目标的范围，然后计算范围分辨率向量中的元素编号。由于目标和雷达是静止的，因此在整个仿真周期中使用位置和速度值。可以假设分辨率向量中的数字在整个模拟中是恒定的。

# параметры движения радара и целей
ant_pos = antennaplatform.InitialPosition;
ant_vel = antennaplatform.Velocity;
tgt_pos = targetplatform[1].InitialPosition;
tgt_vel = targetplatform[1].Velocity;
# Расчет дальности и азимутальный угол от радара до целей
tgt_rng,tgt_ang = rangeangle(tgt_pos,ant_pos);
# поиск номера элемента разрешения по дальности
rangeidx = val2ind(tgt_rng,rangebins[2]-rangebins[1],rangebins[1])[1];

脉冲重叠

让我们创建一个信号处理周期。每个步骤都通过执行系统对象来执行。处理周期被调用两次：一次用于"目标存在"条件，一次用于"目标缺失"条件。

使用EngeePhased将信号辐射到空间。散热器。
使用EngeePhased将信号传播到目标并返回天线。自由空间
使用EngeePhased反射来自目标的信号。目标。
使用EngeePhased将接收信号通过接收放大器。ReceiverPreamp。该步骤还将随机噪声添加到信号中。
使用EngeePhased接收反射信号到天线。收藏家。
匹配放大信号滤波器使用EngeePhased。MatchedFilter。
将匹配滤波器的输出保存在目标范围的bin索引中以供进一步分析。

rcv_pulses = zeros(ComplexF64,length(sigtrans),Npulsebuffsize);
h1 = zeros(ComplexF64,Npulse,1);
h0 = zeros(ComplexF64,Npulse,1);
Nbuff = floor(Int64,Npulse/Npulsebuffsize);
Nrem = Npulse - Nbuff*Npulsebuffsize;
for n = 1:2 # H1 and H0 Hypothesis
    trsig = radiator(sigtrans,tgt_ang);
    trsig = channel[n](trsig,
        ant_pos,tgt_pos,
        ant_vel,tgt_vel);
    global trsig  = trsig
    rcvsig = target[n](trsig);
    rcvsig = collector(rcvsig,tgt_ang);
    
    for k = 1:Nbuff
        for m = 1:Npulsebuffsize
            rcv_pulses[:,m] = receiver(rcvsig, .!(tx_status.>0));
        end
        rcv_pulses = mfilter(rcv_pulses);
        rcv_pulses = buffer(rcv_pulses[(matchingdelay+1):end][:,:],size(rcv_pulses,1));
        if n == 1
            h1[(1:Npulsebuffsize) .+ (k-1)*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,:]);
        else
            h0[(1:Npulsebuffsize) .+ (k-1)*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,:]);
        end
    end
    if (Nrem > 0)
        for m = 1:Nrem
            rcv_pulses[:,m] = receiver(rcvsig, .!(tx_status.>0));
        end
        rcv_pulses = mfilter(rcv_pulses);
        rcv_pulses = buffer(rcv_pulses[(matchingdelay+1):end][:,:],size(rcv_pulses,1));
         if n == 1
            h1[(1:Nrem) .+ Nbuff*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,1:Nrem]);
        else
            h0[(1:Nrem) .+ Nbuff*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,1:Nrem]);
         end
    end
end;

绘制匹配滤波器输出的直方图

让我们计算"目标存在"和"目标不存在"信号的直方图。使用100个样本来获得信号值的扩展的粗略估计。让我们设置直方图值的限制，从最小信号到最大信号。

h1a = abs.(h1);
h0a = abs.(h0);
thresh_low = minimum([h1a;h0a]);
thresh_hi  = maximum([h1a;h0a]);
nbins = 100;
binedges = range(thresh_low,thresh_hi;length=nbins);
histogram(h0a,label="Target Absent",)
histogram!(h1a,label="Target Present")
title!("Target-Absent Vs Target-Present Histograms")

模拟和理论Pd和Pfas的比较

要计算Pd和Pfa，请统计"目标缺失"返回和"目标存在"返回超过指定阈值的情况数。这组阈值比先前模拟中用于创建直方图的区间具有更精细的粒度。然后我们通过脉冲数对这些计算进行归一化以获得概率估计。 Sim_pfa向量是作为阈值函数的虚警的建模概率。 Sim_pd向量是检测的建模概率，它也是阈值的函数。接收器设置阈值，以便它可以确定目标的存在或不存在。上面的直方图显示最佳阈值约为1.8。

nbins = 1000;
thresh_steps = range(thresh_low,thresh_hi;length=nbins)
sim_pd = zeros(1,nbins)
sim_pfa = zeros(1,nbins)
for k = 1:nbins
    thresh = thresh_steps[k]
    sim_pd[k] = sum(h1a .>= thresh)
    sim_pfa[k] = sum(h0a .>= thresh)
end
sim_pd = sim_pd./Npulse;
sim_pfa = sim_pfa./Npulse;

为了构建实验ROC曲线，需要反转Pfa曲线，以便绘制Pd vs Pfa图。只有当您可以将Pfa表示为阈值的严格单调递减函数时，才能反转Pfa曲线。若要以这种方式表示Pfa，请在相邻索引中查找pfa为常量的数组的所有索引。然后从Pd和Pfa阵列中删除这些值。

pfa_diff = diff(sim_pfa[:])
idx = (pfa_diff.!=0).*(1:length(pfa_diff))
sim_pfa = sim_pfa[idx[idx .!= 0]]
sim_pd = sim_pd[idx[idx .!= 0]];

让我们将Pfa的最小值限制为 .

minpfa = 1e-6;
N = sum(sim_pfa .>= minpfa)
sim_pfa = reverse(sim_pfa[1:N])
sim_pd = reverse(sim_pd[1:N]);

计算从最低Pfa到1的pfa和Pd的理论值。然后绘制理论Pfa曲线。

theor_pd,theor_pfa = rocsnr(snr_min,"Pd&Pfa",
    SignalType="NonfluctuatingNoncoherent",
    MinPfa=minpfa,NumPoints=N,NumPulses=1);
plot(theor_pfa,theor_pd,label="Theoretical",legend_position=:bottomright,
    legend_background_color =:transparent)
scatter!(sim_pfa,sim_pd,lab = "Simulated",xscale= :log10,ms=2,
    markerstrokecolor = :red)
title!("Simulated and Theoretical ROC Curves")
xlabel!("Pfa")
ylabel!("Pd")

改进了100万个脉冲的仿真

在之前的模拟中，具有低Pfa的Pd值不会沿着平滑曲线下降，甚至没有达到设定的工作模式。这样做的原因是，在非常低的Pfa水平下，很少（如果有的话）样本超过阈值。要生成具有低Pfa的曲线，必须使用逆Pfa阶的样本数。这种类型的建模需要大量的时间。下面的曲线使用一百万个脉冲而不是45,000个。要运行此模拟，请重复该示例，但将Npulse设置为1000000。

![图像。png]（附件：图片。png)