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

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

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

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

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

using SpecialFunctions,DSP

function findpeaks(x,type,tresh)
    if type == "MinPeakHeight"
        value_peak = map(a -> a > tresh,x).*x
        index = (map(x -> x > tresh,x).*(1:length(x)))
        a = index[index.!=0]
        global i = 1
        out_index = Int[]
        while i <= length(a)
            global i
            if (a[i+1 > length(a) ? end : i+1] - a[i]) == 1
                len_idx = any(index[a[i]:end].==0) ? findfirst(index[a[i]:end].==0)-1 : length(a)-i+1
                push!(out_index,argmax(value_peak[a[i]:a[i+len_idx-1]])+a[i]-1)
                i += len_idx
            else
                push!(out_index,a[i])
                i += 1
            end
        end
        return x[out_index], out_index
    end
end

function noisepow(nbw::Number, nf::Number = 0, reftemp::Number = 290)
    function systemp(nf::Number,reftemp::Number = 290)
        nf >= 0 || error("Invalid input")
        nfLinear = db2pow(nf); 
        stemp = reftemp*(nfLinear); # Kelvin
        return stemp
    end
    B = 1.380649e-23
    npow =  B * systemp(nf,reftemp) * nbw
    return npow 
end

function buffer(x,n,p=0,method="delay")
    method in ["delay","nodelay"] || throw(ArgumentError("Method to be equal `delay` or `nodelay`"))
    (ndims(x) in [1,2] && length(x) == max(size(x)...)) || throw(ArgumentError("Expected x is vector "))
    p < n || throw(ArgumentError("Expected value p < n"))
    x_new = method =="delay" ? [zeros(eltype(x),p);x[:]] : x[:]
    x_buf = reshape([x_new;zeros(n-mod(length(x_new),n))],n,:)
    x_buf = vcat(x_buf,x_buf[end-p+1:end,:])
    n_delay = Int(n-(mod(length(x_buf),n)))*(Int(n-(mod(length(x_buf),n)))!=n) 
    return reshape([x_buf[:] ;zeros(eltype(x),n_delay)],n,:)
end

function fspl(R, lambda)
    (eltype(R) <: Union{Real,Array{<:Real}} &&  all(isnan.(R) .⊻ true ) &&
        all(isempty.(R) .⊻ true) && all(isinf.(R) .⊻ true)) ||
        throw(ArgumentError("Expected R to be real finite vector or scalar."))
    
    (eltype(lambda) <: Real && lambda >= 0 && !isinf(lambda) && (isempty(lambda) ⊻ true) && 
        (isnan(lambda) ⊻ true)) || throw(ArgumentError("Expected lambda to be positive scalar."))
    
    fspl_validate(L) = L < 1 ? 1 : 20*log10(L)
    return fspl_validate.(4pi*R./lambda)
end

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(["Pulse $(i)" for i in 1:num_plot],1,num_plot)
    xlabel_text = [reshape(repeat([""],1,num_plot-1),1,num_plot-1) "Time,(s)"]
    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(["Power (dBw)"],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 bandwidth
pulse_width = 1/pulse_bw; # Pulse width
prf = prop_speed/(2*max_range); # Pulse repetition frequency
fs = 2*pulse_bw; # Sampling rate
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("snr_min = $(round(snr_min;sigdigits=5))")

snr_min = 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 = $(peak_power)")

peak_power = 5226.47386447291

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

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

transmitter = EngeePhased.Transmitter(Gain=tx_gain,PeakPower=peak_power,InUseOutputPort=true);

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

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

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

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

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

antenna = EngeePhased.IsotropicAntennaElement(
    FrequencyRange=[5e9 15e9]); # диапазон частот антенны излучения

sensormotion = EngeePhased.Platform(
    InitialPosition=[0; 0; 0], # начальное положение антенны
    Velocity=[0; 0; 0]); # скорость антенны

С помощью антенны и рабочей частоты мы определяем радиатор и коллектор.

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]];
tgtvel = [[0;0;0] [0;0;0] [0;0;0]];
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;

# Pre-allocate array for improved processing speed
rxpulses = zeros(ComplexF64,length(fast_time_grid),num_pulse_int);

for m = 1:num_pulse_int

    # Update sensor and target positions
    sensorpos,sensorvel = sensormotion(1/prf);
    tgtpos,tgtvel = tgtmotion(1/prf);

    # Calculate the target angles as seen by the sensor
    global tgtrng,tgtang = rangeangle(tgtpos,sensorpos);
    # Simulate propagation of pulse in direction of targets
    pulse = waveform();
    txsig,txstatus = transmitter(pulse);
    txsig = radiator(txsig,tgtang);
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel);
   
    # Reflect pulse off of targets
    tgtsig = target(txsig);
   
    # Receive target returns at sensor
    rxsig = collector(tgtsig,tgtang);
    rxpulses[:,m] = receiver(rxsig,(txstatus.==0));
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))

([1.547656542067221e-6, 1.6136547054810927e-6, 1.2131783640036856e-6], [82, 142, 155])

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

true_range = round.(tgtrng;sigdigits=4)
range_estimates = round.(range_gates[range_detect];sigdigits=4)
println("true_range = $(true_range)")
print("range_estimates = $(range_estimates)")

true_range = [2025.0 3519.0 3845.0]
range_estimates = [2025.0, 3525.0, 3850.0]

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

Заключение

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