Engee 文档
Notebook

误报概率恒定时的检测 (FALT)

本例讨论以恒定误报概率(CFAP)进行检测,并展示如何使用 EngeePhased 库中的 CFARDetectorCFARDetector2D 实现 CFAP 检测算法。

使用的函数

In [ ]:
Pkg.add(["DSP"])
include("$(@__DIR__)/functions.jl");

简介

雷达系统执行的重要任务之一是目标探测。探测程序本身非常简单。它将信号与阈值进行比较。因此,检测的真正工作是确定适当的阈值。一般来说,阈值是正确检测概率和误报概率的函数。

在许多相控阵系统中,由于错误检测会带来一定的成本,因此检测阈值最好不仅能最大限度地提高正确检测的概率,还能将错误报警的概率保持在给定水平以下。

有关如何定义检测阈值的文献很多。读者可能会对白高斯噪声信号检测和多样本信号检测的例子感兴趣,其中给出了一些著名的结果。然而,所有这些经典结果都是基于理论概率,并局限于具有已知方差(功率)的白高斯噪声。在实际设备中,噪声通常是彩色的,其功率也是未知的。

CFAR算法解决了这些问题。当需要检测一个特定的小区(通常称为 "被测小区 "**)时,可通过邻近小区估算噪声功率。在这种情况下,检测阈值 (T) 的定义如下:

$\begin{align} T = \alpha P_n \end{align}$

其中$P_n$ 是估计噪声功率,$\alpha$ 是缩放(阈值)因子

从公式中可以看出,阈值与数据相适应。可以看出,只要有一个适当的阈值系数 α,就能保持误报概率不变,因此被称为 PVLT。

细胞平均。使用 PVLT 进行检测

细胞平均 CFAR 检测器可能是应用最广泛的 CFAR 检测器。它也被用作其他 CFAR 方法的基准比较。在单元平均 CFAR 检测器中,噪声样本是从 CUT 周围的领先和落后单元(称为训练单元)中提取的。噪声估计值的计算公式为

image_2.png 其中$N$ 是训练单元的数量,$x_m$ 是每个训练单元中的样本。

$x_m$ 是二次定律检测器的输出,$P_n$ 是估计的噪声功率。一般来说,前导和滞后训练单元的数量相同。保护单元紧邻 CUT,既在 CUT 的前方,也在 CUT 的后方。这些保护单元的目的是避免信号成分泄漏到训练单元中,从而对噪声估计产生不利影响。

下图显示了一维情况下这些单元之间的关系。

image.png

在上述细胞平均 CFAR 检测器中,假设进入检测器的数据来自单脉冲,即没有脉冲积分,则阈值因子可写成

$\begin{align} \alpha = N(P_{fa}^{-1/N}-1) \end{align}$ 其中$P_{fa}$ 是所需的误报概率。

使用自动阈值因子进行 CFAR 检测

在本示例的其余部分,我们将展示如何使用 EngeePhased 执行细胞平均 CFAR 检测。为简单起见,在不失一般性的前提下,我们仍假设噪声为白高斯噪声。这样,我们就可以将 CFAR 与经典检测理论进行比较。

我们可以使用以下命令创建 CFAR 检测器:

In [ ]:
cfar = EngeePhased.CFARDetector(NumTrainingCells=20,NumGuardCells=2);

在这个检测器中,我们使用 20 个训练单元和 2 个保护单元。这意味着在 CUT 的两侧各有 10 个训练单元和 1 个保护单元。如上所述,如果我们假设信号来自不带脉冲积分的二次定律检测器,那么阈值就可以根据训练单元的数量和所需的误报概率计算出来。假设所需的误报概率为 0.001,我们可以对 CFAR 检测器进行如下配置,以便进行计算。

In [ ]:
exp_pfa = 1e-3 # вероятность ложной тревоги
cfar.ThresholdFactor = "Auto" # алгоритм определения порога (Автоматический)
cfar.ProbabilityFalseAlarm = exp_pfa;

配置后的 CFAR 检测器如下所示。

In [ ]:
cfar
Out[0]:
CFARDetector:
    Method=CA
    Rank=1
    NumGuardCells=2
    NumTrainingCells=20
    ThresholdFactor=Auto
    ProbabilityFalseAlarm=0.001
    CustomThresholdFactor=1
    OutputFormat=CUT result
    ThresholdOutputPort=false
    NoisePowerOutputPort=false
    NumDetectionsSource=Auto
    NumDetections=1

现在我们将模拟输入数据。由于我们想证明 CFAR 检测器可以将误报率控制在一定值以下,因此我们将对这些单元中的噪声进行简单建模。设置如下

  • 数据序列长度为 23 个样本,CUT 为第 12 个单元格。这就在 CUT 的两侧各留下了 10 个训练单元和 1 个保护单元。
  • 误报率通过 100,000 次蒙特卡罗测试计算得出。
In [ ]:
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;

要进行检测,可将数据通过检测器运行。在本例中,只有一个 CUT,因此输出是一个逻辑向量,包含所有试验的检测结果。如果结果为真,则表示目标出现在相应的试验中。在我们的示例中,所有检测结果都是误报,因为我们只通过了噪声。因此,我们可以根据误报次数和试验次数来计算误报率。

$\begin{align} \alpha = N(P_{fa}^{-1/N}-1) \end{align}$

In [ ]:
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials;
println("act_pfa = $(act_pfa)")
act_pfa = 0.001

结果显示,误报概率为 0.001,符合我们的预期。

使用自定义阈值因子进行 CFAR 检测

如本示例前半部分所述,只有在少数情况下 CFAR 检测器可以自动计算适当的阈值因子。例如,在前面的示例中,如果我们在数据到达检测器之前使用 10 脉冲非相干积分,自动阈值将无法提供所需的误报率:

In [ ]:
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")
act_pfa = 0.0

也许有人会不解,为什么我们认为误报率为 0 比误报率为 0.001 更糟糕。毕竟,误报率等于 0 不是一件好事吗?这个问题的答案在于,当误报率降低时,检测到的概率也会降低。在这种情况下,由于真正的误报概率远低于可接受值,因此检测阈值设置得过高。可以用较低的成本(例如较小的发射机功率)实现相同的检测概率和所需的误报概率。

在大多数情况下,阈值系数需要根据具体环境和系统配置来确定。如下图所示,我们可以将 CFAR 检测器配置为使用自定义阈值因子。

In [ ]:
release!(cfar);
cfar.ThresholdFactor = "Custom";

继续以脉冲积分为例,利用经验数据,我们发现可以使用 2.35 的自定义阈值因数来达到所需的误报率。使用该阈值后,我们发现误报率与预期值相符。

In [ ]:
cfar.CustomThresholdFactor = 2.35;
x_detected = cfar(xn,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
act_pfa = 0.00091

CFAR 检测阈值

当单元中的输入信号电平超过阈值电平时,就会进行 CFAR 检测。每个单元的阈值电平取决于阈值因子和该单元中来自训练单元的噪声功率。为了保持误报率不变,检测阈值的增减将与训练单元中的噪声功率成正比。配置 CFAR 检测器,使用 ThresholdOutputPort 属性输出用于每次检测的阈值。让我们设置自动计算阈值系数的模式和 200 个训练单元。

In [ ]:
release!(cfar);
cfar.ThresholdOutputPort = true;
cfar.ThresholdFactor = "Auto";
cfar.NumTrainingCells = 200;

然后生成一个噪声功率不断增大的二次输入信号。

In [ ]:
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;

计算信号中所有单元的检测值和阈值。

In [ ]:
x_detected,th = cfar(xRamp,Vector(1:length(xRamp)));

然后将 CFAR 门限值与输入信号进行比较。

In [ ]:
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)
Out[0]:

在这里,阈值会随着信号中噪声功率的增加而增加,以保持误报水平不变。当信号电平超过阈值时,就会发生检测。

CFAR 与经典的奈曼-皮尔逊检测器的比较

在本节中,我们将利用奈曼-皮尔逊原理比较 CFAR 检测器与经典检测理论的性能。回到第一个例子,假设真实噪声功率已知,理论阈值可计算为

In [ ]:
T_ideal = npower*db2pow(npwgnthresh(exp_pfa));

利用这个理论阈值可以计算出经典奈曼-皮尔逊检测器的误报率。

In [ ]:
act_Pfa_np = sum(x[CUTIdx,:] .> T_ideal)/Ntrials
println("act_Pfa_np = $act_Pfa_np")
act_Pfa_np = 0.00109

由于我们知道噪声功率,因此经典检测理论也能让我们获得所需的误报率。CFAR 检测器实现的误报率与此类似。

In [ ]:
release!(cfar);
cfar.ThresholdOutputPort = false;
cfar.NumTrainingCells = 20;
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
act_pfa = 0.001

接下来,假设两个探测器都部署在现场,噪声功率比预期高 1 dB。在这种情况下,如果我们使用理论阈值,所产生的误报率将是我们想要的四倍。

In [ ]:
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")
act_Pfa_np = 0.00407

相比之下,CFAR 探测器的性能不会受到影响。

In [ ]:
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
act_pfa = 0.00105

因此,CFAR 探测器对噪声功率不确定性具有很强的鲁棒性,更适合现场应用。

最后,我们在存在彩色噪声的情况下使用 CFAR 检测器。我们首先对数据应用经典的检测阈值。

In [ ]:
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")
act_Pfa_np = 0.0

请注意,这样得出的误报率无法满足要求。不过,通过使用带有自定义阈值因子的 CFAR 检测器,我们可以获得所需的误报率。

In [ ]:
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")
act_pfa = 0.00089

通过测距-多普勒肖像进行 CFAR 检测

在前面的章节中,噪声估计是根据在一个维度上领先和落后于 CUT 的训练单元计算出来的。我们也可以在图像上进行 CFAR 检测。单元对应于图像中的像素,保护单元和训练单元放置在 CUT 周围的带状区域中。检测阈值通过 CUT 周围矩形训练带中的单元计算得出。

在上图中,保护带的大小为 [2 2],训练带的大小为 [4 3]。大小指数分别表示 CUT 两侧各行和各列的单元数。保护带的大小也可以定义为 2,因为它在行和列中的大小相同。

接下来,创建一个二维 CFAR 检测器。使用 1e-5 的误报概率,并指定 5 个单元的保护条带大小和 10 个单元的训练条带大小。

In [ ]:
cfar2D = EngeePhased.CFARDetector2D(GuardBandSize=5,TrainingBandSize=10,
  ProbabilityFalseAlarm=1e-5);

然后加载并构建一个测距多普勒图像。图像包括两个静止目标和一个远离雷达的目标的回波。

In [ ]:
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!("Дальностно-доплеровский портрет")
Out[0]:
In [ ]:
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[:]])
Out[0]:
2×82885 Matrix{Int64}:
  41   42   43   44   45   46   47   48  …  155  156  157  158  159  160  161
 171  171  171  171  171  171  171  171     855  855  855  855  855  855  855

计算每个被测单元的探测结果。在本例中,单元是指搜索区域内的每个像素。为测距-多普勒图像绘制检测结果图。

In [ ]:
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")
Out[0]:

检测到三个物体。也可将随时间变化的测距多普勒图像数据立方体作为 cfar2D 的输入,这样就能一步计算出探测结果。

结论

在本例中,我们介绍了 CFAR 探测器背后的基本概念。特别是,我们研究了如何使用EngeePhased工具包在信号和测距多普勒图像上执行细胞平均 CFAR 检测器。比较细胞平均 CFAR 检测器和具有理论计算阈值的检测器的性能,可以清楚地看出 CFAR 检测器更适合实际现场条件。