以恒定的误报概率检测(PVLT)
此示例检查具有恒定误报概率(HLT)的检测，并演示如何使用EngeePhased库中的CFARDetector和CFARDetector2D来实现HLT检测算法。
用过的功能
Pkg.add(["DSP"])
include("$(@__DIR__)/functions.jl");
导言
雷达系统执行的重要任务之一是目标检测。 检测过程本身非常简单。 它将信号与阈值进行比较。 因此，检测的真正工作是确定合适的阈值。 一般而言，阈值是正确检测概率和虚警概率两者的函数。
在许多相控阵系统中，由于与误检相关的成本，期望具有不仅使正确检测的概率最大化，而且使误检的概率保持在预设水平以下的检测阈值。
有关于如何确定检测阈值的广泛文献。 读者可能对白高斯噪声中的信号检测和来自几个样本的信号检测的示例感兴趣，其中呈现了一些众所周知的结果。 然而，所有这些经典结果都基于理论概率，并且仅限于具有已知方差（功率）的白高斯噪声。 在实际设备中，噪声通常是有色的，其功率是未知的。
CFAR算法解决了这些问题。 当需要检测特定小区，通常被称为**测试小区（CUT)**时，从邻近小区估计噪声功率。 在这种情况下，检测阈值(T)如下确定:
哪里 -评估噪音功率, -尺度（阈值）系数
从等式可以看出，阈值适应数据。 可以表明，用适当的阈值系数α，由此产生的虚警概率可以保持在一个恒定的水平，因此命名为PVLT。
小区平均。 用PVLT检测
细胞平均CFAR检测器可能是最广泛使用的CFAR检测器。 它也被用作其他CFAR方法的基础比较。 在细胞平均CFAR检测器中，从切口周围的前导和滞后细胞（所谓的训练细胞）中提取噪声样本。 噪声估计可以计算为:
哪里 $N$ -训练细胞的数目，以及 $x_m$ -在每个训练单元中选择。
如果 是二次定律检测器的输出，则 表示估计的噪声功率。 一般来说，前导和滞后训练小区的数量是相同的。 安全单元位于切口旁边，无论是在前面还是后面。 这些保护小区的目的是避免信号分量泄漏到训练小区中，这可以负面地影响噪声评估。
下图显示了案例1-D中这些单元格之间的关系。
在上述细胞平均CFAR检测器中，假设检测器中的数据来自单个脉冲，即没有脉冲积分，则阈值系数可以写为:
哪里 -虚警的所需概率。
使用自动阈值系数进行CFAR检测
在本示例的其余部分中，我们将展示如何使用EngeePhased执行cell-averaged CFAR检测。 为了简单且不失一般性，我们仍然假设噪声是白高斯。 这使我们能够将CFAR与经典检测理论进行比较。
我们可以使用以下命令创建一个CFAR探测器:
cfar = EngeePhased.CFARDetector(NumTrainingCells=20,NumGuardCells=2);
在该检测器中，我们使用20个训练和2个保护细胞。 这意味着切割的每侧有10个训练单元和1个安全单元。 如上所述，如果我们假设信号来自具有二次定律而没有脉冲积分的检测器，则可以根据训练单元的数量和虚警的期望概率来计算阈值。 假设期望的虚警概率为0.001，我们可以如下配置CFAR检测器，以便可以执行此计算。
exp_pfa = 1e-3 # вероятность ложной тревоги
cfar.ThresholdFactor = "Auto" # алгоритм определения порога (Автоматический)
cfar.ProbabilityFalseAlarm = exp_pfa;
配置的CFAR检测器如下所示。
cfar
现在我们将模拟输入数据。 由于我们想表明CFAR探测器可以将误报水平保持在某个值以下，我们将简单地模拟这些单元格中的噪声。 以下是设置:
*数据序列为23个样本长，切割为细胞12。 这在切口的每一侧留下10个训练单元和1个安全单元。
*使用100,000蒙特卡洛测试计算误报率。
rng = MersenneTwister(20)
npower = db2pow(-10); # Задание ОСШ - 10dB
Ntrials = Int(1e5);
Ncells = 23;
CUTIdx = 12;
# Модель шума после квадратичного детектора
rsamp = randn(rng,Ncells,Ntrials)+1im*randn(rng,Ncells,Ntrials);
x = abs.(sqrt.(npower/2)*rsamp).^2;
为了执行检测，将数据通过检测器。 在这个例子中，只有一个剪切，因此输出是包含所有测试的检测结果的逻辑向量。 如果结果为真，则表示目标存在于相应的测试中。 在我们的例子中，所有的检测都是假警报，因为我们只通过噪声。 结果，可以基于虚警次数和试验次数来计算虚警率。
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials;
println("act_pfa = $(act_pfa)")
结果表明，如我们所假设的，由此产生的虚警概率为0.001。
使用自定义阈值系数检测CFAR
如本示例的前面部分所解释的，CFAR检测器可以自动计算适当的阈值系数的情况仅有少数。 例如，在前面的场景中，如果我们在向检测器接收数据之前使用10脉冲非相干积分，则自动阈值将不再能够提供所需的误报率。:
npower = db2pow(-10); # значение ОСШ - 10 дБ
xn = 0;
for m = 1:10
rsamp = randn(rng,Ncells,Ntrials) .+ 1im*randn(rng,Ncells,Ntrials);
xn = xn .+ abs.(sqrt.(npower/2)*rsamp).^2 # некогерентное накопление
end
x_detected = cfar(xn,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
可能令人困惑的是，为什么我们认为0的虚警率比0.001的虚警率更差。 毕竟，0的虚警率不是一件好事吗？ 这个问题的答案在于，随着虚警概率的降低，检测概率也会降低。 在这种情况下，由于虚警的真实概率远低于可接受的值，所以检测阈值设置得太高。 相同的检测概率可以以较低的成本实现期望的虚警概率；例如，以较低的发射机功率。
在大多数情况下，必须根据具体环境和系统配置来确定阈值系数。 我们可以将CFAR检测器配置为使用自定义阈值系数，如下所示。
release!(cfar);
cfar.ThresholdFactor = "Custom";
继续查看脉冲积分的示例并使用经验数据，我们发现可以使用2.35的自定义阈值因子来实现所需的误报率。 利用该阈值，我们可以看到得到的虚警系数与期望值相对应。
cfar.CustomThresholdFactor = 2.35;
x_detected = cfar(xn,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
CFAR检测阈值
CFAR检测发生在单元中的输入信号电平超过阈值电平时。 每个小区的阈值水平取决于阈值系数和从训练小区获得的该小区中的噪声功率。 为了保持恒定的误报水平，检测阈值将与训练小区中的噪声功率成比例地增加或减少。 配置CFAR检测器以使用ThresholdOutputPort属性输出用于每次检测的阈值。 我们将为阈值系数和200个训练单元格设置自动计算模式。
release!(cfar);
cfar.ThresholdOutputPort = true;
cfar.ThresholdFactor = "Auto";
cfar.NumTrainingCells = 200;
然后，我们将生成噪声功率增加的二次输入信号。
Npoints = 10_000;
rsamp = randn(rng,Npoints).+1im*randn(rng,Npoints);
ramp = Vector(range(1,10;length=Npoints));
xRamp = abs.(sqrt.(npower*ramp./2).*rsamp).^2;
计算信号中所有像元的检测和阈值。
x_detected,th = cfar(xRamp,Vector(1:length(xRamp)));
然后将CFAR阈值与输入信号进行比较。
gr()
plot(1:length(xRamp),xRamp,label="Сигнал",line=(2,:blue,:sticks))
plot!(1:length(xRamp),th,label="Порог")
scatter!(findall(x_detected[:] .> 0),xRamp[findall(x_detected[:] .> 0)],label="Обнаружение",marker=(:circle,4,:yellow))
xlabel!("Временные отсчеты")
ylabel!("Амплитуда",legend=:topleft)
这里，阈值随着信号中噪声功率的增加而增加，以保持恒定的虚警电平。 检测发生在信号强度超过阈值的地方。
CFAR与经典Neiman-Pearson探测器的比较
在本节中，我们将使用Neumann-Pearson原理将CFAR探测器的性能与经典检测理论进行比较。 回到第一个例子，假设真正的噪声功率是已知的，理论阈值可以计算为
T_ideal = npower*db2pow(npwgnthresh(exp_pfa));
这个经典的Neiman-Pearson探测器的虚警系数可以使用这个理论阈值来计算。
act_Pfa_np = sum(x[CUTIdx,:] .> T_ideal)/Ntrials
println("act_Pfa_np = $act_Pfa_np")
由于我们知道噪声功率，经典检测理论也使我们能够获得所需的虚警系数。 CFAR检测器实现的误报率相似。
release!(cfar);
cfar.ThresholdOutputPort = false;
cfar.NumTrainingCells = 20;
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
接下来，假设两个探测器都部署在现场，噪声功率比预期高1dB。 在这种情况下，如果我们使用理论阈值，则由此产生的虚警概率将比我们想要的大四倍。
npower = db2pow(-9); # Assume 9dB SNR ratio
rsamp = randn(rng,Ncells,Ntrials).+1im*randn(rng,Ncells,Ntrials);
x = abs.(sqrt.(npower/2)*rsamp).^2;
act_Pfa_np = sum(x[CUTIdx,:] .> T_ideal)/Ntrials
println("act_Pfa_np = $act_Pfa_np")
相反，不影响CFAR检测器的操作。
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
因此，CFAR检测器能够抵抗噪声功率不确定性，更适合现场应用。
最后，我们在存在颜色噪声的情况下使用CFAR检测。 首先，我们将经典检测阈值应用于数据。
npower = db2pow(-10);
fcoeff = [-0.00057693500095478268512028119374691, 0.0004446399942713034555974438433168, 0.023749844987587823141872434007382, 0.10475951999236174372320817838045, 0.22682709001336692766770397611253, 0.2895916800267339108465591834829, 0.22682709001336692766770397611253, 0.10475951999236174372320817838045, 0.023749844987587819672425482053768, 0.0004446399942713034555974438433168, -0.00057693500095478268512028119374691]
x = abs.(sqrt(npower/2)*filt(fcoeff,1,rsamp)).^2; # colored noise
act_Pfa_np = sum(x[CUTIdx,:] .> T_ideal)/Ntrials
println("act_Pfa_np = $act_Pfa_np")
请注意，获得的虚警系数可能不符合要求。 然而，通过使用具有自定义阈值系数的CFAR检测器，我们可以获得所需的虚警系数。
release!(cfar);
cfar.ThresholdFactor = "Custom";
cfar.CustomThresholdFactor = 12.85;
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
CFAR检测范围-多普勒肖像
在前面的章节中，噪声估计是基于在一个维度上领先和落后于切割的训练单元计算的。 我们还可以对图像进行CFAR检测。 像元对应于图像中的像素，并且保护和训练像元在切口周围以条纹排列。 切周围的矩形学习带中的像元计算检测阈值。
上图中，保护条的尺寸为[2 2]，训练条的尺寸为[4 3]。 大小索引分别指示行和列中每个切割边上的单元格数。 保护条的尺寸也可以定义为2，因为它的尺寸在行和列中是相同的。
接下来，创建二维CFAR检测器。 使用1e-5虚警概率并指定5小区保护条带的大小和10小区训练条带的大小。
cfar2D = EngeePhased.CFARDetector2D(GuardBandSize=5,TrainingBandSize=10,
ProbabilityFalseAlarm=1e-5);
然后下载并构建范围多普勒图像。 图像包括来自两个静止目标和一个远离雷达的目标的返回。
hrdresp = EngeePhased.RangeDopplerResponse(
DopplerFFTLengthSource="Property",
DopplerFFTLength=RangeDopplerEx_MF_NFFTDOP,
SampleRate=RangeDopplerEx_MF_Fs);
resp,rng_grid,dop_grid = step!(hrdresp,RangeDopplerEx_MF_X,RangeDopplerEx_MF_Coeff);
#Добавление белого шума
npow = db2amp(-95);
resp = resp .+ npow/sqrt(2)*(randn(rng,size(resp)) .+ 1im*randn(rng,size(resp)))
plot_1 = heatmap(dop_grid,rng_grid,mag2db.(abs.(resp)),colorbar=false,margin=5Plots.mm,
color = :rainbow1, titlefont = font(12,"Computer Modern"),
guidefont = font(10,"Computer Modern"))
resp = abs.(resp).^2;
xlabel!("Частота Доплера, Гц")
ylabel!("Дальность, м")
title!("Дальностно-доплеровский портрет")
rangeIndx = map(x -> x[1],(findmin(abs.(rng_grid[:] .- [1000 4000]);dims=1)[2]))
dopplerIndx = map(x -> x[1],(findmin(abs.(dop_grid[:] .- [-1e4 1e4]);dims=1)[2]))
columnInds,rowInds = meshgrid(dopplerIndx[1]:dopplerIndx[2],rangeIndx[1]:rangeIndx[2]);
CUTIdx = permutedims([rowInds[:] columnInds[:]])
计算每个被测小区的检测结果。 在这个例子中，搜索区域中的每个像素是一个像元。 为范围-多普勒图像构建检测结果的图。
detections = cfar2D(resp,CUTIdx)
detectionMap = zeros(size(resp))
detectionMap[rangeIndx[1]:rangeIndx[2],dopplerIndx[1]:dopplerIndx[2]] =
reshape(float(detections),rangeIndx[2]-rangeIndx[1]+1,dopplerIndx[2]-dopplerIndx[1]+1);
heatmap(dop_grid,rng_grid,mag2db.(abs.(detectionMap .+ 1eps(1.))),colorbar=false,c = :imola,
titlefont = font(12,"Computer Modern"),guidefont = font(10,"Computer Modern"),margin=5Plots.mm)
xlabel!("Частота Доплера, Гц")
ylabel!("Дальность, м")
title!("Дальностно-доплеровский портрет после обнаружения CFAR")
已检测到三个物体。 作为cfar2D的输入信号，还可以提供随时间变化的测距仪多普勒图像数据的立方体，并且检测将在一个步骤中计算。
结论
在此示例中，我们介绍了CFAR探测器的基本概念。 特别是，我们研究了如何使用EngeePhased工具包对信号和测距仪多普勒图像执行细胞平均CFAR检测器。 对具有像元平均的CFAR检测器和具有理论计算阈值的检测器的特性的比较清楚地表明，CFAR检测器更适合于真实的现场条件。