Кривые обнаружения (ROC) часть 2: Моделирование методом Монте-Карло

Открыть пример в Engee

В этом примере показано, как с помощью моделирования Монте-Карло построить кривую обнаружения приемника (ROC-кривая) радиолокационной системы. ROC определяет, насколько хорошо система может обнаруживать цели, отбрасывая при этом большие значения ложных сигналов, когда цель отсутствует (ложные тревоги). Система обнаружения определяет наличие или отсутствие цели, сравнивая полученное значение сигнала с заданным порогом.

Вероятность правильного обнаружения (Pd) - это вероятность того, что мгновенное значение сигнала больше порогового значения, когда цель действительно присутствует.

Вероятность ложной тревоги (Pfa) - это вероятность того, что значение сигнала больше порога, когда цель отсутствует. В этом случае сигнал обусловлен шумом, и его свойства зависят от статистики шума. При моделировании методом Монте-Карло генерируется очень большое количество радиолокационных сигналов в условиях наличия и отсутствия цели.

ROC-кривая представляет собой график зависимости Pd от Pfa. Форма ROC-кривой зависит от отношения сигнал шум (ОСШ) принимаемого сигнала. Если ОСШ принимаемого сигнала известен, то кривая ROC показывает насколько надежна система, анализируя значения Pd и Pfa. Если задать Pd и Pfa, то можно определить, какая мощность необходима для выполнения этого требования.

Вы можете использовать функцию rocsnr для расчета теоретических ROC-кривых. В этом примере показана ROC-кривая, полученная в результате моделирования методом Монте-Карло одноантенной радарной системы, и проведено сравнение этой кривой с теоретической кривой.

Тактико-технические характеристики радара

Установим желаемую Pd = 0,9 и Pfa = $10^{-6}$, зададим максимальную дальность действия радара 4000 метров и разрешение по дальности 50 метров. Установите фактическую дальность до цели 3000 метров. Также установим эффективную площадь рассеяния (ЭПР) цели в 1.5 $м^2$ и рабочую частоту в 10 ГГц.

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

Выбор тип зондирующего сигнала (ЗС) и его параметры

Сформируем системный объект генератора ЛЧМ-сигнала. Для этого необходимо вычислить:

• Полосу импульсов, используя разрешение по дальности;
• Частоту повторения импульсов (ЧПИ), исходя из максимальной дальности работы радара.
• Частоту дискретизации (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, необходимое для достижения заданной вероятности ложной тревоги и вероятности обнаружения, используя уравнение Альбершема.

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

ОСШ = 13.11 дБ

Чтобы достичь ОСШ = 13 дБ, на цель должна быть передана достаточная мощность. Используйте уравнение радара, чтобы оценить пиковую мощность передачи, необходимую для достижения заданного ОСШ для цели на дальности 3000 метров. Принимаемый сигнал также зависит от ЭПР, которое, как предполагается, соответствует нефлутуирующей модели (Swerling 0). Установим для радара одинаковые коэффициенты усиления по приему и передаче, равные 20 дБ.

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

# Расчет дальности с помощью основного уравнения радилокации
peak_power = ((4*pi)^3*noisepow(1/pulse_duration)*target_range^4*
    db2pow(snr_min))/(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 # несущая частота
);

Настройка объекта цели

Создадим объект Target, соответствующий реальной отражающей цели с заданным ЭПР. Отражения от этой цели будут имитировать реальные отраженные сигналы. Чтобы рассчитать ложные тревоги, создадим вторую цель с нулевой ЭПР. Отражения от этой цели равны нулю, за исключением шума.

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 значение "Noise temperature", а для свойства шумовой температуры 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);

Создайте коэффициенты согласованного фильтра из объекта waveform . Затем создайте системный объект согласованного фильтра.

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];

Перекрытие импульсов

Создадим цикл обработки сигнала. Каждый шаг выполняется путем выполнения системных объектов. Цикл обработки вызывается дважды: один раз для условия "цель присутствует" и один раз для условия "цель отсутствует".

1. Излучение сигнала в пространство с помощью EngeePhased.Radiator.
2. Распространение сигнала к цели и обратно к антенне с помощью EngeePhased.FreeSpace
3. Отражение сигнала от цели с помощью EngeePhased.Target.
4. Прохождение принятого сигнала через усилитель приема с помощью EngeePhased.ReceiverPreamp.Этот шаг также добавляет случайный шум к сигналу.
5. Прием отраженного сигнала на антенну с помощью EngeePhased.Collector.
6. Согласованный фильтр усиленного сигнала с помощью EngeePhased.MatchedFilter.
7. Сохраните выход согласованного фильтра в индексе бина целевого диапазона для дальнейшего анализа.

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, подсчитайте количество случаев, когда возврат "цель отсутствует" и возврат "цель присутствует" превышают заданный порог. Этот набор пороговых значений имеет более тонкую гранулярность, чем бины, использованные для создания гистограммы в предыдущем моделировании. Затем нормализуем эти подсчеты по количеству импульсов, чтобы получить оценку вероятности. Вектор 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 как строго монотонно убывающую функцию от Thresh. Чтобы выразить 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 до $10^{-6}$ .

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

Вычислим теоретические значения Pfa и Pd от наименьшего Pfa до 1. Затем постройте теоретическую кривую 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")

Улучшение моделирования с помощью одного миллиона импульсов

В предыдущем моделировании значения Pd при низком Pfa не падают по плавной кривой и даже не доходят до заданного рабочего режима. Причина этого в том, что при очень низких уровнях Pfa очень немногие выборки, если таковые вообще имеются, превышают пороговое значение. Для генерации кривых при низком Pfa необходимо использовать количество выборок порядка, обратного Pfa. Этот тип моделирования занимает много времени. Следующая кривая использует миллион импульсов вместо 45 000. Чтобы запустить это моделирование, повторите пример, но установите Npulse равным 1000000.