Моделирование тестовых сигналов для радиолокационного приемника

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

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

Этот пример посвящен разработке импульсной радиолокационной системы, которая может достичь заданных технических характеристик. В нем описаны шаги по переводу проектных спецификаций, таких как вероятность обнаружения и разрешение по дальности, в параметры радарной системы - мощность передатчика и длительность импульса. Также моделируется окружающая среда и цели для получения отраженного сигнала. Наконец, к полученному сигналу применяются методы обработки сигнала для определения дальности до целей.

Используемые функции

using DSP

# Функция для построения выходного отклика радара и порогового значения
function RadarPulsePlot(ReceivePulse,threshold,t,tp,num_plot::Int64)
    thresh = sqrt(threshold)*ones(length(t),num_plot)
    ReceivePulse =  Float64.(pow2db.((abs.(ReceivePulse)).^2))
    ReceivePulse[ReceivePulse.==0] .= eps(0.)
    thresh .=  pow2db.(thresh.^2)
    t = repeat(t,1,num_plot) .+ tp[1:num_plot]'
    text_title = reshape(["Импульс № $(i)" for i in 1:num_plot],1,num_plot)
    xlabel_text = [reshape(repeat([""],1,num_plot-1),1,num_plot-1) "Время, с"]
    plot(t,ReceivePulse[:,1:num_plot],layout=(num_plot,1),label="",title = text_title,
        fontfamily="Computer Modern", titlefontsize=10, guidefontsize=10, tickfontsize=10,gridalpha=0.3,lw=1.2)
    plot!(t,thresh[:,1:num_plot],layout=(num_plot,1),xlabel=xlabel_text,lw=1.2,
        ylabel = repeat(["Мощность, дБВт"],1,num_plot),label="",ls=:dash) # 
end;

Технические характеристики устройства

Целью разработки данной импульсной радиолокационной системы является обнаружение не флуктуирующих целей с поперечным сечением радара (RCS) не менее одного квадратного метра на расстоянии до 5000 метров от радара с разрешением по дальности 50 метров. Желаемый показатель эффективности - вероятность обнаружения (Pd) 0,9 и вероятность ложной тревоги (Pfa) менее 1e-6. Поскольку когерентное обнаружение требует фазовой информации и, следовательно, требует больших вычислительных затрат, мы используем некогерентную схему обнаружения. Кроме того, в данном примере предполагается, что среда находится в свободном пространстве.

pd = 0.9; # Вероятность верного обнаружения
pfa = 1e-6; # Вероятность ложной тревоги
max_range = 5000; # Максимальная дальность
range_res = 50; # Требуемое разрешение по дальности
tgt_rcs = 1; # ЭПР цели

Проектирование моностатической радарной системы

Нам необходимо определить несколько характеристик радарной системы, таких как форма волны, приемник, передатчик и антенна, используемая для излучения и приема сигнала.

Генератор

В этом примере мы выбрали прямоугольную форму волны. Желаемое разрешение диапазона определяет полосу пропускания формы сигнала, которая, в случае прямоугольной формы сигнала, определяет ширину импульса.

Другим важным параметром импульсной формы волны является частота повторения импульсов (PRF). PRF определяется максимальным однозначным диапазоном.

prop_speed = 299_792_458; # Скорость распространения
pulse_bw = prop_speed/(2*range_res);  # Полоса импульса в частотной области
pulse_width = 1/pulse_bw; # Длительность импульса
prf = prop_speed/(2*max_range); # Частота следования импульсов
fs = 2*pulse_bw; # Частота дискретизации системы
waveform = EngeePhased.RectangularWaveform(PulseWidth=1/pulse_bw,PRF=prf,SampleRate=fs); # системый объект

Обратите внимание, что мы установили частоту дискретизации, равную удвоенной ширине полосы пропускания.

Шумовые характеристики приемника

Мы предполагаем, что единственным шумом, присутствующим в приемнике, является тепловой шум, поэтому в данном моделировании нет никаких помех. Мощность теплового шума зависит от полосы пропускания приемника. Полоса пропускания шумов приемника устанавливается равной полосе пропускания формы сигнала. Так часто бывает в реальных системах. Мы также предполагаем, что приемник имеет коэффициент усиления 20 дБ и коэффициент шума 0 дБ.

noise_bw = pulse_bw; # шумовая полоса приемника
receiver = EngeePhased.ReceiverPreamp(
    Gain=20, # усиление 
    NoiseFigure=0,# фактор шума
    SampleRate=fs, # частота дискретизации
    EnableInputPort=true # порт для синхронизации передатчика и приемника
);

Обратите внимание, что, поскольку мы моделируем моностатический радар, приемник не может быть включен, пока не выключен передатчик. Поэтому мы устанавливаем свойство EnableInputPort в true, чтобы сигнал синхронизации мог быть передан от передатчика к приемнику.

Передатчик

Наиболее важным параметром передатчика является пиковая мощность. Требуемая пиковая мощность связана со многими факторами, включая максимальную однозначную дальность, требуемый SNR в приемнике и ширину импульса. Среди этих факторов требуемый SNR в приемнике определяется целью проектирования Pd и Pfa, а также схемой обнаружения, реализованной в приемнике.

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

snr_db = [-Inf, 0.0, 3, 10, 13]; # значения ОСШ в дБ
rocsnr(snr_db,SignalType="NonfluctuatingNoncoherent") # построение характеристик обнаружения

Кривые ROC показывают, что для удовлетворения проектных целей Pfa = 1e-6 и Pd = 0,9; SNR принимаемого сигнала должен превышать 13 дБ. Это довольно высокое требование и не очень практичное. Чтобы сделать радилокационную систему более реализуемой, мы можем использовать технику накопления импульсов для снижения требуемого SNR. Если мы выберем интегрирование 10 импульсов, то кривая может быть построена следующим образом:

num_pulse_int = 10;
rocsnr([0 3 5],SignalType="NonfluctuatingNoncoherent",NumPulses=num_pulse_int)

Мы видим, что требуемая мощность снизилась примерно до 5 дБ. Дальнейшее снижение SNR может быть достигнуто за счет интеграции большего количества импульсов, но количество импульсов, доступных для интеграции, обычно ограничено из-за движения цели или неоднородности окружающей среды.

В приведенном выше подходе значение SNR считывается с кривой, но часто бывает желательно вычислить только требуемое значение. Для некогерентной схемы обнаружения расчет требуемого SNR теоретически довольно сложен. К счастью, существуют хорошие приближения, такие как уравнение Альбершема. Используя уравнение Альбершема, требуемое значение SNR можно вывести как

snr_min = albersheim(pd, pfa,N= num_pulse_int) # расчет ОСШ путем решения уравнения Альбршема
print("Минимальное ОСШ: $(round(snr_min;sigdigits=5)) дБ")

Минимальное ОСШ: 4.9904 дБ

Получив требуемое значение SNR на приемнике, можно рассчитать пиковую мощность на передатчике с помощью уравнения радара. Здесь мы предполагаем, что коэффициент усиления передатчика составляет 20 дБ.

Чтобы рассчитать пиковую мощность с помощью уравнения радара, нам также необходимо знать длину волны распространяющегося сигнала, которая связана с рабочей частотой системы. Здесь мы задаем рабочую частоту 10 ГГц.

tx_gain = 20; # Усиление передатчика
fc = 10e9; # несущая частота сигнала
lambda = prop_speed/fc; # длина волны сигнала
peak_power = ((4*pi)^3*noisepow(1/pulse_width)*max_range^4*
    db2pow(snr_min))/(db2pow(2*tx_gain)*tgt_rcs*lambda^2) # пиковая мощность
print("Пиковая мощность: $(peak_power*1e3) кВт")

Пиковая мощность: 5.22647386447291e6 кВт

Обратите внимание, что результирующая мощность составляет около 5 кВт, что является приемлемым значением. Для сравнения, если бы мы не использовали метод интеграции импульсов, то пиковая мощность составила бы 33 кВт, что является довольно большим значением.

Получив всю эту информацию, мы можем настроить передатчик.

transmitter = EngeePhased.Transmitter(
    Gain=tx_gain, # усиление
    PeakPower=peak_power, # пиковая мощность
    InUseOutputPort=true # включение выходного порта, который определяет 
                         # состояние передатчика (1 - вкл, 0 - выкл)
);

Поскольку в данном примере моделируется моностатическая радарная система, порт InUseOutputPort устанавливается в true для вывода состояния передатчика. Этот сигнал состояния может быть использован для включения приемника.

Передающее и приемное устройства

В радарной системе сигнал распространяется в виде электромагнитной волны. Поэтому сигнал должен быть излучен и собран антенной, используемой в радарной системе. Именно здесь в дело вступают излучающая и "собирающая" антенна.

В моностатической радарной системе излучатель и коллектор используют одну и ту же антенну, поэтому сначала мы определим антенну. Для упрощения конструкции мы выбрали изотропную антенну. Обратите внимание, что антенна должна работать на рабочей частоте системы (10 ГГц), поэтому мы зададим диапазон частот антенны 5-15 ГГц.

Мы предполагаем, что антенна неподвижна.

# антенный элемент радара
antenna = EngeePhased.IsotropicAntennaElement(
    FrequencyRange=[5e9 15e9]); # диапазон частот антенны излучения
# модель движения радара
sensormotion = EngeePhased.Platform(
    InitialPosition=[0; 0; 0], # начальное положение антенны
    Velocity=[0; 0; 0]); # скорость антенны

С помощью антенны и рабочей частоты задаем излучатель (Radiator) и принимающую антенну (Collector).

radiator = EngeePhased.Radiator(
    Sensor=antenna,             # задание антенного элемента 
    OperatingFrequency=fc);   # задание несущей частоты

collector = EngeePhased.Collector(
    Sensor=antenna,             # задание  антенного элемента 
    OperatingFrequency=fc);     # задание несущей частоты

На этом конфигурация радарной системы завершена. В следующих разделах мы определим другие объекты, такие как цель и окружающая среда, которые необходимы для моделирования. Затем мы смоделируем возврат сигнала и выполним определение дальности по смоделированному сигналу.

Моделирование системы

Цели

Чтобы проверить способность нашего радара обнаруживать цели, мы должны сначала определить цели. Предположим, что в пространстве есть 3 стационарные, не флуктуирующие цели. Их положения и радиолокационные сечения приведены ниже.

tgtpos = [[2024.66;0;0] [3518.63;0;0] [3845.04;0;0]]; # компоненты начального положения 3 целей [[px] [py] [pz]]
tgtvel = [[0;0;0] [0;0;0] [0;0;0]]; # компоненты вектора скорости 3 целей [[vx] [vy] [vz]]
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # модель движения целей

tgtrcs = [1.6 2.2 1.05]; # ЭПР целей
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc); # Системный объект целей

Среда распространения

Чтобы смоделировать сигнал, нам также необходимо определить канал распространения между радарной системой и каждой целью.

channel = EngeePhased.FreeSpace(
    SampleRate=fs, # частота дискретизации
    TwoWayPropagation=true, # двунаправленное распространение
    OperatingFrequency=fc); # несущая частота

Поскольку в данном примере используется моностатическая радарная система, каналы настроены на моделирование задержек при распространении в двух направлениях.

Синтез сигналов

Теперь мы готовы к моделированию всей системы.

Синтезированный сигнал представляет собой матрицу данных с быстрым временем (время внутри каждого импульса) в каждом столбце и медленным временем (время между импульсами) в каждой строке. Для визуализации сигнала полезно определить как сетку быстрого времени, так и сетку медленного времени.

fast_time_grid = (0:1/fs:1/prf-1/fs); # сетка быстрого времени
slow_time_grid = (0:num_pulse_int-1)./prf; # сетка медленного времени

Следующий цикл моделирует 10 импульсов принимаемого сигнала.

Мы задаем затравку для генерации шума в приемнике, чтобы можно было воспроизвести те же результаты.

release!(receiver) # доступ к обновлению параметров объекта
receiver.SeedSource = "Property";  # изменение тип модели генерации шума
receiver.Seed = 2007; # задание начального состояния генератора тепловых шумов

# Выделение памяти для результирующих данных
rxpulses = zeros(ComplexF64,length(fast_time_grid),num_pulse_int);

for m = 1:num_pulse_int
    
    # Обновление пложения радара и целей
    sensorpos,sensorvel = sensormotion(1/prf);
    tgtpos,tgtvel = tgtmotion(1/prf);

    # Расчет углов по дальности от радара до целей  
    global tgtrng,tgtang = rangeangle(tgtpos,sensorpos);
    
    # Моделирование распространения сигнала от радара до целей
    pulse = waveform(); # генератор прямоугольных импульсво
    txsig,txstatus = transmitter(pulse); # усилитель передающего сигнала
    txsig = radiator(txsig,tgtang); # излучатель сгенерированного сигнала в среду
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel); # среда распространения
    
    # Отражение сигнала от целей с учетом ЭПР
    tgtsig = target(txsig);
    
    # Прием и предварительного усиление отраженного сигнала
    rxsig = collector(tgtsig,tgtang); # приемная антенна
    rxpulses[:,m] = receiver(rxsig,xor.(txstatus.>0,1)); # предусилитель
end

Определение дальности

Порого обнаружения

Детектор сравнивает мощность сигнала с заданным порогом. В радарных приложениях порог часто выбирается таким образом, чтобы Pfa была ниже определенного уровня. В данном случае мы предполагаем, что шум является белым гауссовским, а обнаружение некогерентным. Поскольку мы также используем 10 импульсов для интегрирования импульсов, пороговая мощность сигнала определяется следующим образом

npower = noisepow(noise_bw,receiver.NoiseFigure,
    receiver.ReferenceTemperature);  # мощность шума
threshold = npower * db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # порог обнаружения

Мы строим график первых двух принятых импульсов с порогом

num_pulse_plot = 2; # количество графиков
RadarPulsePlot(rxpulses,threshold,
    fast_time_grid,slow_time_grid,num_pulse_plot)

Порог на этих рисунках используется только для демонстрации. Обратите внимание, что сигналы от второй и третьей цели намного слабее, чем от первой, поскольку они находятся дальше от радара. Поэтому мощность принимаемого сигнала зависит от дальности, и порог несправедлив к целям, расположенным на разных расстояниях.

Согласованная фильтрация

Согласованный фильтр обеспечивает усиление обработки, которое повышает порог обнаружения. Он преобразует принятый сигнал в локальную, обращенную во времени и сопряженную копию переданной формы волны. Поэтому при создании согласованного фильтра мы должны указать форму передаваемого сигнала. Принятые импульсы сначала пропускаются через согласованный фильтр, чтобы улучшить SNR, прежде чем выполнять интегрирование импульсов, пороговое обнаружение и т. д.

matchingcoeff = getMatchedFilter(waveform); # расчет коэффициентов согласованного 
                                            # фильтра на основе изулченного сигнала          
matchedfilter = EngeePhased.MatchedFilter(
    Coefficients=matchingcoeff, # коэффициенты фильтра
    GainOutputPort=true); # выходной порт усиления
rxpulses, mfgain = matchedfilter(rxpulses); # отклик после согласованной фильтрации

Согласованный фильтр вносит внутреннюю задержку, поэтому местоположение пика (выходной образец с максимальным SNR) больше не совпадает с истинным местоположением цели. Чтобы компенсировать эту задержку, перенесем выход согласованного фильтра вперед и добавим нули в конце. Обратите внимание, что в реальных системах, поскольку данные накапливаются непрерывно, окончания массива не будет.

matchingdelay = size(matchingcoeff,1)-1; # длина коэффициентов согласованного фильтра
rxpulses = buffer(rxpulses[matchingdelay+1:end],200); # size(rxpulses,1) ->200

Затем порог увеличится на коэффициент усиления обработки согласованного фильтра.

threshold = threshold * db2pow(mfgain);

На следующем графике показаны те же два импульса после прохождения через согласованный фильтр.

RadarPulsePlot(rxpulses,threshold,
    fast_time_grid,slow_time_grid,num_pulse_plot)

После этапа согласованного фильтра SNR улучшается. Однако, поскольку мощность принимаемого сигнала зависит от дальности, сигнал от близкой цели все равно намного сильнее, чем от цели, находящейся на большем расстоянии. Поэтому, как показано на рисунке выше, шум от близкого бина имеет значительный шанс превысить порог и затмить цель, находящуюся дальше. Чтобы пороговое значение было справедливым для всех целей в пределах обнаруживаемой дальности, мы можем использовать изменяющийся во времени коэффициент усиления, чтобы компенсировать зависящие от дальности потери в принимаемом эхосигнале.

Чтобы компенсировать потери, зависящие от дальности, мы сначала вычисляем импульс строба, соответствующие каждому образцу сигнала, а затем рассчитываем потери в свободном пространстве, соответствующие каждому импульсу строба . После получения этой информации мы применяем изменяющееся во времени усиление к принимаемому импульсу, чтобы возвраты были как будто из одного и того же эталонного диапазона (максимальная обнаруживаемая дальность).

range_gates = prop_speed*fast_time_grid/2; # разрешение по дальности
tvg = EngeePhased.TimeVaryingGain( # создание системного объекта ВАРУ
    RangeLoss=2*fspl(range_gates,lambda), # потери с учетом дальности
    ReferenceLoss=2*fspl(max_range,lambda)); # относительные потери при максимальной дальности
rxpulses = tvg(rxpulses); # выходной сигнал после работы ВАРУ

Теперь построим график тех же двух импульсов после нормализации диапазона

RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,num_pulse_plot)

Изменение коэффициента усиления во времени приводит к увеличению уровня шума. Однако возврат цели теперь не зависит от диапазона. Теперь для обнаружения во всем обнаруживаемом диапазоне можно использовать постоянный порог.

Обратите внимание, что на этом этапе порог находится выше максимального уровня мощности, содержащегося в каждом импульсе. Поэтому на данном этапе пока ничего нельзя обнаружить. Нам необходимо выполнить интегрирование импульсов, чтобы мощность возвращаемых эхо-сигналов от целей превысила порог, а уровень шума остался ниже планки. Это ожидаемо, поскольку именно интеграция импульсов позволяет нам использовать последовательность импульсов меньшей мощности.

Некогерентное накопление

Мы можем еще больше улучшить SNR, некогерентно интегрируя (видеоинтегрируя) полученные импульсы.

rxpulses = pulsint(rxpulses,"noncoherent"); # некогерентное накопление
RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,1)

После этапа интеграции видео данные готовы к финальному этапу обнаружения. На рисунке видно, что все три эхо-сигнала от целей превышают пороговое значение, а значит, могут быть обнаружены.

Обнаружение дальности

Наконец, для интегрированных импульсов выполняется пороговое детектирование. Схема обнаружения идентифицирует пики, а затем переводит их положение в диапазоны целей.

_,range_detect = findpeaks(rxpulses,"MinPeakHeight",sqrt(threshold)); # обнаружение пиков

Истинные и обнаруженные дальности целей показаны ниже:

true_range = round.(tgtrng;sigdigits=4)
range_estimates = round.(range_gates[range_detect];sigdigits=4)
println("Истинная дальность целей: $(true_range) м")
print("Оценка дальности целей: $(range_estimates) м")

Истинная дальность целей: [2025.0 3519.0 3845.0] м
Оценка дальности целей: [2025.0, 3525.0, 3850.0] м

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

Заключение

В этом примере мы спроектировали радарную систему на основе набора заданных целей. Исходя из этих целей, были рассчитаны многие конструктивные параметры радиолокационной системы. В примере также показано, как использовать разработанный радар для выполнения задачи обнаружения дальности. В этом примере радар использовал прямоугольную форму волны.