Выбор типа зондирующего сигнала с целью улучшения дальности действия радилокационной системы

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

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

Pkg.add(["DSP"])

   Resolving package versions...
   Installed CUDA_Driver_jll ──── v0.10.4+0
   Installed CUDA_Runtime_jll ─── v0.15.5+0
   Installed MLDataDevices ────── v1.6.11
   Installed MPItrampoline_jll ── v5.5.2+0
   Installed NVTX ─────────────── v0.3.5
   Installed Atomix ───────────── v1.1.0
   Installed NNlib ────────────── v0.9.27
   Installed IntervalArithmetic ─ v0.22.23
   Installed CUDA ─────────────── v5.6.1
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

using DSP
# Функция для построения выходного отклика радара и порогового значения
function RadarPulsePlot(ReceivePulse,threshold,t,tp,num_plot::Int64)
    thresh = sqrt(threshold)*ones(length(t),num_plot)
    ReceivePulse =  Float64.(DSP.db2pow.((abs.(ReceivePulse)).^2))
    ReceivePulse[ReceivePulse.==0] .= eps(0.)
    thresh .=  DSP.db2pow.(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;

function BasicMonostaticRadarExampleData()
    c = 299_792_458
    waveform = EngeePhased.RectangularWaveform(PulseWidth=100/c, 
        PRF=c/1e4,OutputFormat="Pulses",NumPulses=1,SampleRate= 2c/100);
    antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[5e9 1.5e10]);
    radiator = EngeePhased.Radiator(Sensor=antenna,
        PropagationSpeed=c,OperatingFrequency=1e10);
    collector = EngeePhased.Collector(Sensor=antenna,
        PropagationSpeed=c,Wavefront="Plane",
        OperatingFrequency=1e10);
    transmitter = EngeePhased.Transmitter(PeakPower=5226.4791642003665401716716587543,
        Gain=20,LossFactor=0, InUseOutputPort=true,CoherentOnTransmit=true);
    receiver = EngeePhased.ReceiverPreamp(Gain=20,NoiseFigure=0, 
        ReferenceTemperature=290,SampleRate=2*c/100, PhaseNoiseInputPort=false,
        EnableInputPort=true,SeedSource="Property",Seed=1000);
    sensormotion = EngeePhased.Platform(InitialPosition=[0;0;0],Velocity=[0;0;0]);

    return waveform, antenna, radiator, collector, transmitter, receiver, sensormotion
end;

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

В примере "Моделирование тестовых сигналов для радиолокационного приемника" разработан моностатический импульсный радар, имитирующий тестовые сигналы для радиолокационного приемника для удовлетворения следующий требования:

1. Минимальная эффективная площадь рассеяния (ЭПР): 1 $м^2$, не флуктуирующая;
2. Максимальная дальность однозначного определения цели: 5 км;
3. Вероятность правильного обраружения (ВПО): 0.9;
4. Вероятность ложной тревоги (ВЛТ): $10^{-6}$.

Загрузим системные объекты РЛС:

waveform,antenna,radiator,collector,transmitter,receiver,sensormotion = BasicMonostaticRadarExampleData();

Новые требования

После развертывания РЛС, добавим 2 новых требования:

1. Максимальная дальность однозначного определения цели: 5 км;
2. Модель флутуации целей: "Swerling-2"

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

Первым затронутым параметром является частота следования импульсов (ЧСИ). Его необходимо пересчитать на основе нового максимальной однозначной дальности обнаружения.

prop_speed = radiator.PropagationSpeed; # скорость распространения
max_range = 8000; # максимальная дальость
prf = prop_speed/(2*max_range) # ЧСИ
print("Новая частота следования импульсов: $(round(prf*1e-3;sigdigits=5)) кГц")

Новая частота следования импульсов: 18.737 кГц

По сравнению с 30 кГц ЧСИ существующей конструкции, новая ЧСИ стала меньше и равна 18.737 кГц. Следовательно, интервал между импульсами больше. Обратите внимание, что это тривиальное изменение в программном обеспечении радара и довольно дешевое по стоимости оборудования.

waveform.PRF = prf;

Далее, поскольку цель описывается с помощью модели Swerling-2, необходимо использовать уравнение Шнидмана (вместо уравнения Альбершема) для расчета требуемого ОСШ для достижения заданных ВПО и ВЛТ. Уравнение Шнидмана предполагает использование квадратичного детектора с некогерентным накоплением. Число импульсов для накопления равно 10.

num_pulse_int = 10; # количество импульсов накопления
pfa = 1e-6; # вероятность ложной тревоги
snr_min = shnidman(0.9,pfa,num_pulse_int,2) # расчет ОСШ
print("Минимальное ОСШ: $(round(snr_min;sigdigits=5)) дБ")

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

2. Выбор тип сигнала

Если бы мы использовали ту же прямоугольную форму сигнала в существующей конструкции, ширина импульса осталась бы прежней, поскольку она определяется разрешением по дальности. Однако, наша максимальная дальность действия увеличилась с 5 км до 8 км, а целевая модель переключилась с модели без флуктуаций на модель "Swerling-2", нам необходимо пересчитать требуемую пиковую мощность передачи.

fc = radiator.OperatingFrequency; # несущая частота
lambda = prop_speed/fc; # длина волны
peak_power = ((4*pi)^3*noisepow(1/waveform.PulseWidth)*max_range^4*
    DSP.db2pow(snr_min))/(DSP.db2pow(2*transmitter.Gain)*1*lambda^2) # расчет пиковой мощности
print("Пиковая мощность передатчика: $(round(peak_power*1e-3;sigdigits=5)) кВт")

Пиковая мощность передатчика: 44.821 кВт

Пиковая мощность примерно в восемь раз превышает предыдущее требование. Такое улучшение оборудования не так просто реализовать по 2 причинам:

1. Существующее оборудование радара рассчитано на получение импульса с пиковой мощностью около 5200 Вт. Хотя большинство конструкций оставляют некоторый запас выше требуемой мощности, маловероятно, что существующая система может обеспечить в восемь раз большую мощность;
2. Замена оборудования для вывода такой высокой мощности обходится очень дорого. Следовательно, текущая конструкция должна быть изменена в соответствии с новой целью за счет использования более сложных методов обработки сигналов.

Сигнал с линейной частотной модуляцией (ЛЧМ)

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

Давайте рассмотрим детали конфигурации сигнала с ЛЧМ.

Желаемое разрешение диапазона определяет полосу пропускания сигнала. Для ЛЧМ-сигнала полоса пропускания равна его ширине развертки. Однако ширина импульса больше не ограничена величиной, обратной ширине полосы пропускания импульса, поэтому можно использовать гораздо большую ширину импульса. Мы используем ширину импульса, которая в 20 раз больше, и устанавливаем частоту дискретизации, которая в два раза превышает ширину полосы пропускания импульса.

range_res = 50; # разрешение по дальности
pulse_bw = prop_speed/(2*range_res); # ширина спектра ЛЧМ-сигнала
pulse_width = 20/pulse_bw; # длительность импульса ЛЧМ-сигнала
fs = 2*pulse_bw; # частота дискретизации

waveform = EngeePhased.LinearFMWaveform(
    SweepBandwidth=pulse_bw,
    PulseWidth=pulse_width,
    PRF=prf,
    SampleRate=fs
);

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

peak_power = ((4*pi)^3*noisepow(1/waveform.PulseWidth)*max_range^4*
    DSP.db2pow(snr_min))/(DSP.db2pow(2*transmitter.Gain)*1*lambda^2)
    
print("Пиковая мощность передатчика: $(round(peak_power*1e-3;sigdigits=5)) кВт")

Пиковая мощность передатчика: 2.2411 кВт

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

release!(transmitter)
transmitter.PeakPower = peak_power;

3. Системное моделирование

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

Цели

Как и в случае с вышеупомянутым примером, мы предполагаем, что в среде свободного пространства есть 3 цели. Однако теперь моделью цели является модель флуктуации цели - "Swerling-2". Зададим положения радара и целей, и ЭПР целей

tgtpos = [[2024.66; 0; 0] [6518.63; 0; 0] [6845.04; 0; 0]];  # компоненты начального положения целей
tgtvel = [[0;0;0] [0;0;0] [0;0;0]]; # компоненты скорости целей
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # модель движения целей

tgtrcs = [2.2 1.1 1.05]; # ЭПР целей
fc = radiator.OperatingFrequency; # несущая частота
target = EngeePhased.RadarTarget(
    Model="Swerling2", # модель флуктуации
    MeanRCS=tgtrcs, # среднее значение ЭПР цели
    OperatingFrequency=fc # несущая частота
);

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

Далее, настроим канал распространения между радаром и каждой целью.

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

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

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

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

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

# Выделение памяти для результирующих данных
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,true);
    
    # Прием и предварительного усиление отраженного сигнала
    rxsig = collector(tgtsig,tgtang); # приемная антенна
    rxpulses[:,m] = receiver(rxsig, .!(txstatus.>0) ); # предусилитель
end

4. Определение оценки дальности

Расчет порога обнаружения

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

noise_bw = receiver.SampleRate/2;
npower = noisepow(noise_bw,
    receiver.NoiseFigure,receiver.ReferenceTemperature);
threshold = npower * db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent"));

pulseplotnum = 2;
RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,pulseplotnum)

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

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

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

matchingcoeff = getMatchedFilter(waveform);
matchedfilter = EngeePhased.MatchedFilter(
    CoefficientsSource="Property",
    Coefficients=matchingcoeff,
    GainOutputPort=true
);
rxpulses, mfgain = matchedfilter(rxpulses);
threshold = threshold * DSP.db2pow(mfgain);

Для компенсации задержки в согласованном фильтре, воспользуеся фукнкцией buffer.

matchingdelay = size(matchingcoeff,1)-1;
rxpulses = buffer(rxpulses[(matchingdelay+1):end][:,:],size(rxpulses,1));

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

range_gates = prop_speed*fast_time_grid/2; 
lambda = prop_speed/fc;

tvg = EngeePhased.TimeVaryingGain(
    RangeLoss=2*fspl(Vector(range_gates),lambda),
    ReferenceLoss=2*fspl(max_range,lambda));
rxpulses = tvg(rxpulses);

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

Теперь,необходимо накапливать принимаемые импульсы некогерентно для дальнейшего улучшения ОСШ.

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)
println("Истинная дальность: $(true_range) м")

Истинная дальность: [2025.0 6519.0 6845.0] м

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

Оценка дальности: [2025.0, 6525.0, 6850.0] м

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

Заключение

В этом примере мы использовали тип сигнала с ЛЧМ для определения дальности действия. Используя такой тип сигнала, получилось снизить требуемую пиковую мощность передачи, достигнув большей дальности обнаружения (в 8 км) для модели флуктуации цели - Swerling-2.