В этом примере показано, как сгенерировать кривую рабочей характеристики приемника (ROC) радиолокационной системы с помощью симуляции Монте-Карло. Рабочая характеристика приемника определяет, как хорошо система может обнаружить цели при отклонении больших побочных значений сигналов, когда цель отсутствует (ложные предупреждения). Система обнаружения объявит присутствие или отсутствие цели путем сравнения полученного значения сигналов с предварительно установленным порогом. Вероятность обнаружения (Pd) цели является вероятностью, что мгновенное значение сигналов больше, чем порог каждый раз, когда цель на самом деле присутствует. Вероятность ложного предупреждения (PFA) является вероятностью, что значение сигналов больше, чем порог, когда цель отсутствует. В этом случае сигнал происходит из-за шума, и его свойства зависят от шумовой статистики. Симуляция Монте-Карло генерирует очень большое количество радара, возвращается с и без существующей цели. Симуляция вычисляет Pd, и PFA путем подсчета пропорции значений сигналов в каждом случае, которые превышают порог.
Кривая ПТИЦЫ РУХ строит Pd в зависимости от PFA. Форма кривой ROC зависит от полученного ОСШ сигнала. Если прибывающий ОСШ сигнала известен, то кривая ROC показывает, как хорошо система выполняет в терминах Pd и PFA. Если вы задаете Pd и PFA, то можно определить, сколько степени необходимо, чтобы достигнуть этого требования.
Можно использовать функциональный rocsnr
вычислить теоретические кривые ROC. Этот пример показывает кривую ROC, сгенерированную симуляцией Монте-Карло радиолокационной системы одно антенны, и сравнивает ту кривую с теоретической кривой.
Установите желаемую вероятность обнаружения быть 0.9 и вероятность ложного предупреждения, чтобы быть . Установите максимальную область значений радара к 4 000 метров и разрешения области значений 50 метров. Установите фактический целевой диапазон на 3 000 метров. Установите целевое радарное поперечное сечение на 1,5 квадратных метра и установите рабочую частоту на 10 ГГц. Все расчеты выполняются в основной полосе.
c = physconst('LightSpeed');
pd = 0.9;
pfa = 1e-6;
max_range = 4000;
target_range = 3000.0;
range_res = 50;
tgt_rcs = 1.5;
fc = 10e9;
lambda = c/fc;
Любая симуляция, которая вычисляет PFA и фунт, требует обработки многих сигналов. Чтобы поддержать требования к памяти на низком уровне, обработайте сигналы во фрагментах импульсов. Определите номер импульсов к процессу к 45 000 и установите размер каждого фрагмента к 10 000.
Npulse = 45000; Npulsebuffsize = 10000;
Вычислите ширину полосы пропускания импульса формы волны с помощью импульсного разрешения области значений. Вычислите импульсную частоту повторения из максимальной области значений. Поскольку сигнал является основной полосой, установите частоту дискретизации на дважды пропускную способность. Вычислите импульсную длительность от ширины полосы пропускания импульса.
pulse_bw = c/(2*range_res); prf = c/(2*max_range); fs = 2*pulse_bw; pulse_duration = 10/pulse_bw; waveform = phased.LinearFMWaveform('PulseWidth',pulse_duration,... 'SampleRate',fs,'SweepBandwidth',... pulse_bw,'PRF',prf);
Достижение конкретного Pd и PFA требует, чтобы достаточная степень сигнала прибыла в приемник после того, как цель отражает сигнал. Вычислите минимальный ОСШ, должен был достигнуть заданной вероятности ложного предупреждения и вероятности обнаружения при помощи уравнения Albersheim.
snr_min = albersheim(pd,pfa);
К достигнуть этого ОСШ, достаточная степень должна быть передана к цели. Используйте основное уравнение радиолокации, чтобы оценить пиковую степень передачи, peak_power
, требуемый достигнуть заданного ОСШ в дБ для цели в области значений 3 000 метров. Полученный сигнал также зависит от целевого радарного поперечного сечения (RCS). который принят, чтобы следовать не колеблющейся модели (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)
peak_power = 293.1830
Создание системных объектов, которые составляют часть передачи симуляции: радарная платформа, антенна, передатчик и излучатель.
antennaplatform = phased.Platform(... 'InitialPosition',[0; 0; 0],... 'Velocity',[0; 0; 0]); antenna = phased.IsotropicAntennaElement(... 'FrequencyRange',[5e9 15e9]); transmitter = phased.Transmitter(... 'Gain',txrx_gain,... 'PeakPower',peak_power,... 'InUseOutputPort',true); radiator = phased.Radiator(... 'Sensor',antenna,... 'OperatingFrequency',fc);
Создайте целевой Системный объект, соответствующий фактической цели отражения с ненулевым целевым поперечным сечением. Отражения от этой цели симулируют фактический радар, возвращается. Для того, чтобы вычислить ложные предупреждения, создайте второй целевой Системный объект с нулевым радарным сечением. Отражения от этой цели являются нулем за исключением шума.
target{1} = phased.RadarTarget(... 'MeanRCS',tgt_rcs,... 'OperatingFrequency',fc); targetplatform{1} = phased.Platform(... 'InitialPosition',[target_range; 0; 0]); target{2} = phased.RadarTarget(... 'MeanRCS',0,... 'OperatingFrequency',fc); targetplatform{2} = phased.Platform(... 'InitialPosition',[target_range; 0; 0]);
Смоделируйте среду распространения от радара до целей и назад.
channel{1} = phased.FreeSpace(... 'SampleRate',fs,... 'TwoWayPropagation',true,... 'OperatingFrequency',fc); channel{2} = phased.FreeSpace(... 'SampleRate',fs,... 'TwoWayPropagation',true,... 'OperatingFrequency',fc);
Заданный шум путем установки NoiseMethod
свойство к 'Noise temperature'
и ReferenceTemperature
свойство к 290 K.
collector = phased.Collector(... 'Sensor',antenna,... 'OperatingFrequency',fc); receiver = phased.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);
Создайте коэффициенты согласованного фильтра из Системного объекта формы волны. Затем создайте Систему согласованного фильтра object™.
MFCoeff = getMatchedFilter(waveform); matchingdelay = size(MFCoeff,1) - 1; filter = phased.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));
Создайте цикл обработки сигналов. Каждый шаг выполняется путем выполнения Системных объектов. Цикл обрабатывает импульсы дважды, однажды для целевого текущего состояния и однажды для целевого отсутствующего условия.
Излучите сигнал в пробел с помощью phased.Radiator
.
Распространите сигнал к цели и назад к антенне с помощью phased.FreeSpace
.
Отразите сигнал от цели с помощью phased.Target
.
Получите отраженные сигналы в антенне с помощью phased.Collector
.
Передайте полученный сигнал хотя получить усилитель с помощью phased.ReceiverPreamp
. Этот шаг также добавляет случайный шум в сигнал.
Фильтр соответствия усиленный сигнал с помощью phased.MatchedFilter
.
Сохраните согласованный фильтр выход в индексе интервала целевого диапазона для последующего анализа.
rcv_pulses = zeros(length(sigtrans),Npulsebuffsize); h1 = zeros(Npulse,1); h0 = zeros(Npulse,1); Nbuff = floor(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); 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 = filter(rcv_pulses); rcv_pulses = buffer(rcv_pulses(matchingdelay+1:end),size(rcv_pulses,1)); if n == 1 h1((1:Npulsebuffsize) + (k-1)*Npulsebuffsize) = rcv_pulses(rangeidx,:).'; else h0((1:Npulsebuffsize) + (k-1)*Npulsebuffsize) = rcv_pulses(rangeidx,:).'; end end if (Nrem > 0) for m = 1:Nrem rcv_pulses(:,m) = receiver(rcvsig,~(tx_status>0)); end rcv_pulses = filter(rcv_pulses); rcv_pulses = buffer(rcv_pulses(matchingdelay+1:end),size(rcv_pulses,1)); if n == 1 h1((1:Nrem) + Nbuff*Npulsebuffsize) = rcv_pulses(rangeidx,1:Nrem).'; else h0((1:Nrem) + Nbuff*Npulsebuffsize) = rcv_pulses(rangeidx,1:Nrem).'; end end end
Вычислите гистограммы целевых существующих и целевых отсутствующих возвратов. Используйте 100 интервалов, чтобы дать грубую оценку распространения значений сигналов. Установите область значений значений гистограммы от самого маленького сигнала до самого большого сигнала.
h1a = abs(h1); h0a = abs(h0); thresh_low = min([h1a;h0a]); thresh_hi = max([h1a;h0a]); nbins = 100; binedges = linspace(thresh_low,thresh_hi,nbins); figure histogram(h0a,binedges) hold on histogram(h1a,binedges) hold off title('Target-Absent Vs Target-Present Histograms') legend('Target Absent','Target Present')
Чтобы вычислить Pd и PFA, вычислите количество экземпляров, что целевой отсутствующий возврат и целевой существующий возврат превышают заданный порог. Этот набор порогов имеет более прекрасную гранулярность, чем интервалы раньше создавали гистограмму в предыдущей симуляции. Затем нормируйте эти количества на количество импульсов, чтобы получить оценку вероятностей. Векторный sim_pfa
симулированная вероятность ложного предупреждения в зависимости от порога, thresh
. Векторный sim_pd
симулированная вероятность обнаружения, также функция порога. Приемник устанавливает порог так, чтобы это могло определить, присутствует ли цель или отсутствует. Гистограмма выше предполагает, что лучший порог - приблизительно 1,8.
nbins = 1000; thresh_steps = linspace(thresh_low,thresh_hi,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 является константой по соседним индексам. Затем удалите эти значения из массивов PFA и Pd.
pfa_diff = diff(sim_pfa); idx = (pfa_diff == 0); sim_pfa(idx) = []; sim_pd(idx) = [];
Ограничьте самый маленький PFA .
minpfa = 1e-6; N = sum(sim_pfa >= minpfa); sim_pfa = fliplr(sim_pfa(1:N)).'; sim_pd = fliplr(sim_pd(1:N)).';
Вычислите теоретический PFA и значения Pd от самого маленького PFA до 1. Затем постройте теоретическую кривую PFA.
[theor_pd,theor_pfa] = rocsnr(snr_min,'SignalType',... 'NonfluctuatingNoncoherent',... 'MinPfa',minpfa,'NumPoints',N,'NumPulses',1); semilogx(theor_pfa,theor_pd) hold on semilogx(sim_pfa,sim_pd,'r.') title('Simulated and Theoretical ROC Curves') xlabel('Pfa') ylabel('Pd') grid on legend('Theoretical','Simulated','Location','SE')
В предыдущей симуляции значения Pd в низком PFA не падают вдоль плавной кривой и даже не расширяют вниз к заданному операционному режиму. Причина этого на очень низких уровнях PFA, очень немногих, если таковые имеются, выборки превышают порог. Чтобы сгенерировать кривые в низком PFA, необходимо использовать много выборок порядка инверсии PFA. Этот тип симуляции занимает много времени. Следующая кривая использует один миллион импульсов вместо 45 000. Чтобы запустить эту симуляцию, повторите пример, но установите Npulse
к 1000000.