Обнаружение при постоянной вероятности ложной тревоги (ПВЛТ)¶
В этом примере рассматривается обнаружение c постоянной вероятностью ложной тревоги (ПВЛТ) и показано, как использовать CFARDetector и CFARDetector2D из библиотеки EngeePhased для реализации алгоритма обнаружения с ПВЛТ
Используемые фукции¶
Pkg.add(["DSP"])
include("$(@__DIR__)/functions.jl");
Введение¶
Одна из важных задач, которую выполняет радарная система - обнаружение цели. Сама процедура обнаружения довольно проста. Он сравнивает сигнал с пороговым значением. Поэтому настоящая работа по обнаружению заключается в определении соответствующего порога. В общем случае, порог является функцией как вероятности правильного обнаружения, так и вероятности ложной тревоги.
Во многих системах с фазированными решетками из-за затрат, связанных с ложным обнаружением, желательно иметь порог обнаружения, который не только максимизирует вероятность правильного обнаружения, но и удерживает вероятность ложной тревоги ниже заданного уровня.
Существует обширная литература о том, как определить порог обнаружения. Читателям могут быть интересны примеры обнаружения сигнала в белом гауссовском шуме и обнаружения сигнала по нескольким образцам, где приведены некоторые хорошо известные результаты. Однако все эти классические результаты основаны на теоретических вероятностях и ограничены белым гауссовским шумом с известной дисперсией (мощностью). В реальных устройствах шум часто бывает цветным, а его мощность неизвестна.
Алгоритм CFAR решает эти проблемы. Когда требуется обнаружить определенную ячейку, часто называемую тестируемой ячейкой (CUT), мощность шума оценивается по соседним ячейкам. В таком случае порог обнаружения (T) определяется следующим образом:
$\begin{align} T = \alpha P_n \end{align}$
где $P_n$ - мощность шума оценки, $\alpha$ - масштабный (пороговый) коэффициент
Из уравнения видно, что порог адаптируется к данным. Можно показать, что при соответствующем пороговом коэффициенте α результирующая вероятность ложной тревоги может поддерживаться на постоянном уровне, отсюда и название ПВЛТ.
Усреднение по ячейкам. Обнаружение с ПВЛТ¶
Детектор CFAR с усреднением по ячейкам, вероятно, является наиболее широко используемым детектором CFAR. Он также используется в качестве базового сравнения для других методов CFAR. В детекторе CFAR с усреднением по ячейкам образцы шума извлекаются из опережающих и отстающих ячеек (так называемых обучающих ячеек) вокруг CUT. Оценка шума может быть вычислена как:
где $N$ - количество обучающих ячеек, а $x_m$ - выборка в каждой обучающей ячейке.
Если $x_m$ является выходом детектора квадратичного закона, то $P_n$ представляет собой оценочную мощность шума. В общем случае количество опережающих и отстающих обучающих ячеек одинаково. Защитные ячейки располагаются рядом с CUT, как впереди, так и позади нее. Цель этих защитных ячеек - избежать просачивания компонентов сигнала в обучающую ячейку, что может негативно повлиять на оценку шума.
На следующем рисунке показана связь между этими ячейками для случая 1-D.
В приведенном выше детекторе CFAR с усреднением по ячейкам, предполагая, что данные в детектор поступают от одного импульса, т.е. без интегрирования импульсов, пороговый коэффициент можно записать в виде:
$\begin{align} \alpha = N(P_{fa}^{-1/N}-1) \end{align}$ где $P_{fa}$ - требуемая вероятность ложной тревоги.
Обнаружение CFAR с помощью автоматического порогового коэффициента¶
В остальной части этого примера мы покажем, как использовать EngeePhased для выполнения усредненного по ячейкам обнаружения CFAR. Для простоты и без потери общности мы по-прежнему предполагаем, что шум является белым гауссовским. Это позволяет сравнить CFAR с классической теорией обнаружения.
Мы можем создать детектор CFAR с помощью следующей команды:
cfar = EngeePhased.CFARDetector(NumTrainingCells=20,NumGuardCells=2);
В этом детекторе мы используем 20 обучающих и 2 защитных ячейки. Это означает, что с каждой стороны CUT находится 10 обучающих и 1 защитная ячейка. Как упоминалось выше, если мы предполагаем, что сигнал поступает от детектора с квадратичным законом без интегрирования импульсов, порог может быть рассчитан на основе количества обучающих ячеек и желаемой вероятности ложной тревоги. Предполагая, что желаемая вероятность ложной тревоги составляет 0,001, мы можем настроить детектор CFAR следующим образом, чтобы можно было выполнить этот расчет.
exp_pfa = 1e-3 # вероятность ложной тревоги
cfar.ThresholdFactor = "Auto" # алгоритм определения порога (Автоматический)
cfar.ProbabilityFalseAlarm = exp_pfa;
Настроенный детектор CFAR показан ниже.
cfar
Теперь мы смоделируем входные данные. Поскольку мы хотим показать, что детектор CFAR может поддерживать уровень ложных тревог ниже определенного значения, мы просто смоделируем шум в этих ячейках. Вот настройки:
- Последовательность данных имеет длину 23 образца, а CUT - это ячейка 12. Таким образом, остается 10 обучающих ячеек и 1 защитная ячейка с каждой стороны CUT.
- Коэффициент ложной тревоги рассчитывается с помощью 100 тысяч испытаний Монте-Карло.
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}$
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 с классическим детектором Неймана-Пирсона¶
В этом разделе мы сравним производительность детектора CFAR с классической теорией обнаружения, используя принцип Неймана-Пирсона. Возвращаясь к первому примеру и предполагая, что истинная мощность шума известна, теоретический порог может быть рассчитан как
T_ideal = npower*db2pow(npwgnthresh(exp_pfa));
Коэффициент ложной тревоги этого классического детектора Неймана-Пирсона может быть рассчитан с использованием этого теоретического порога.
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")
Далее предположим, что оба детектора развернуты в поле и мощность шума на 1 дБ больше, чем ожидалось. В этом случае, если мы используем теоретический порог, результирующая вероятность ложной тревоги будет в четыре раза больше, чем мы хотим.
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 по дальностно-доплеровскому портрету¶
В предыдущих разделах оценка шума вычислялась на основе обучающих ячеек, опережающих и отстающих от CUT в одном измерении. Мы также можем выполнять обнаружение CFAR на изображениях. Ячейки соответствуют пикселям на изображениях, а защитные и обучающие ячейки располагаются в полосах вокруг CUT. Порог обнаружения вычисляется по ячейкам в прямоугольной полосе обучения вокруг CUT.
На рисунке выше размер защитной полосы составляет [2 2], а размер обучающей полосы - [4 3]. Индексы размеров означают количество ячеек на каждой стороне CUT в строке и столбце соответственно. Размер защитной полосы можно также определить как 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 больше подходит для реальных полевых условий.