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

本例介绍如何使用蒙特卡罗模拟绘制雷达系统的接收器检测曲线（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 米。同时将目标的有效散射面积 (EDA) 设为 1.5 ，工作频率设为 10 GHz。

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; # количество отсчетов в импульсе

选择探测信号 (PS) 类型及其参数

让我们组成 LFM 信号发生器的系统对象。为此需要计算
** 脉冲波段**，使用范围分辨率；
** 脉冲重复频率（PRF），以雷达最大探测距离为基础。
** 采样率 ** (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 和 Pfa，需要在目标反射信号后向接收器提供足够的信号功率。让我们利用阿尔伯森方程计算达到给定误报概率和检测概率所需的最小信噪比。

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

ОСШ = 13.11 дБ

要达到 13 dB 的信噪比，就必须向目标发射足够的功率。使用雷达方程估算出在 3000 米距离上对目标实现给定误报概率所需的峰值发射功率。接收到的信号也取决于 EPR，假定 EPR 遵循非波动模型（Swerling 0）。让我们将雷达的接收增益和发射增益设置为 20 dB。

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 的真实反射目标。该目标的反射将模拟真实的反射信号。要计算误报，请创建第二个 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属性设置为 290 K，即可产生噪声。

# Геометрия приемной антенны
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.Radiator将信号辐射到空间。
2. 使用EngeePhased.FreeSpace将信号传播到目标并返回天线。
3. 使用EngeePhased.Target将信号从目标反射出去。
4. 使用EngeePhased.ReceiverPreamp将接收到的信号通过接收放大器。
5. 使用EngeePhased.Collector在天线上接收反射信号。
6. 使用EngeePhased.MatchedFilter对放大后的信号进行匹配滤波。
7. 7. 将匹配滤波器的输出保存在目标频段的 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 和 Pfa 与理论 Pd 和 Pfa 的比较

要计算 Pd 和 Pfa，需要计算 "目标不存在 "返回值和 "目标存在 "返回值超过给定阈值的次数。这组阈值的粒度比之前模拟中用于创建直方图的粒度更细。然后，我们用脉冲数对这些计数进行归一化处理，以获得概率估计值。矢量 sim_pfa 是模拟误报概率与阈值 thresh 的函数关系。矢量 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 与 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")

使用一百万个脉冲改进建模

在之前的模拟中，低 Pfa 时的 Pd 值不会沿着平滑的曲线下降，甚至不会达到指定的工作模式。原因是在极低的 Pfa 水平下，只有极少数样本（如果有的话）超过阈值。要生成低 Pfa 时的曲线，必须使用与 Pfa 的倒数相同数量级的样本。这种建模方式非常耗时。下面的曲线使用了一百万个脉冲而不是 45,000 个脉冲。要运行此模拟，请重复示例，但将 Npulse 设为 1000000。

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

雷达的战术和技术特性

选择探测信号 (PS) 类型及其参数

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

设置目标对象

设置传播介质对象的参数

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

设置快速时间网格

探测信号的产生

计算到目标的范围

脉冲重叠

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

模拟 Pd 和 Pfa 与理论 Pd 和 Pfa 的比较

使用一百万个脉冲改进建模

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