Многоканальная широкополосная система пеленгации
Многоканальная широкополосная система пеленгации малоразмерных объектов
В примере рассматривается моделирование многоканальной радиолокационной системы для решения задачи пеленгации множественных целей. В качестве целевых малоразмерных объектов выступают беспилотные летательные аппараты (БПЛА) - "квадрокоптер".
Перед выполнением скрипта рекомендуется очистить рабочую область с помощью команды engee.clear(). Если это не требуется можете удалить или закомментировать данную команду. Также, предоставляется выбор типа визуализации графика: интерактивный или статический. Управление реализуется с помощью переменной is_dinamic_plot.
engee.clear() # очистка рабочей области
is_dinamic_plot = true # выбор типа графика:true (plotlyjs()) - интерактивный, false (gr()) - статический
is_dinamic_plot ? plotlyjs() : gr()
Подключение вспомогательных файлов и библиотек
Для реализации рассматриваемой системы необходимо подключить сторонние библиотеки и вспомогательные функции:
let
installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["JLD2","Hungarian"]
for pack in list_packages
pack in installed_packages || Pkg.add(pack)
end
end
using LinearAlgebra,Hungarian,DataFrames
import JLD2
include("utils/plot_sig_and_spec.jl") # фунция построения осциллограммы и спектрограммы
include("utils/viz_scene.jl") # фунция построения сценария работы РЛС
include("utils/calc_plot_rcs.jl") # фунция построения диаграммы обратного рассеяния
include("utils/sort_ang_mse.jl") # фунция сортироки оценок пеленга
include("utils/plot_spec_azel.jl") # функция построения пространственного спектра
include("utils/calc_table_ang_analysis.jl") # функция расчета таблицы для анализа точности
1. Разработка системной модели сценария работы РЛС
Рассмотрим сценарий работы многоканальной широкополосной системы пеленгации при наличии нескольких целей.
- Считается, что РЛС расположена в начале координат, ориентирована вдоль оси x. Для решения задачи определения пеленга используется прямоугольная фазированная антенная решетка (ФАР) и алгоритм цифровой обработки (MUSIC, MVDR, BeamScan)
- По направлению к РЛС против оси x движутся четыре малоразмерных объекта БПЛА типа "Квадракоптер". Каждый из объектов имеет собственную траекторию движения. Дальность до целей составляет 1000-1500 метров, скорость 130-150 м/с. Азимутальные направления на цели лежат в пределах [-20, 20] град, угломестные - [0, 30] град.
Примерный облик сценария приведен на рисунке ниже.

1.1 Описание структуры РЛС
Опишем более детально структурную схему системной модели РЛС:
-
Генератор ЗС: в качестве зондирующего сигнала используется импульсный сигнал с линейной частотной модуляцией (ЛЧМ);
-
Передатчик: в передатчике усиливается сигнал с учетом потерей в передающем тракте;
-
Передающая ФАР: формируется диаграмма направленности (ДН) с помощью прямоугольной ФАР;
-
Среда распространения: учитывает затухание на распространение, частотную дисперсию, влияние погоды;
-
Цель: описывает ДОР с учетом широкополосности спектра сигнала и траекторию движения цели;
-
Передающая ФАР: формируется ДН с помощью приемной широкополосной прямоугольной ФАР,;
-
Приемник: реализует усиление принятого сигнала с добавлением собственных тепловых шумов (модель МШУ).
-
Алгоритм пеленга: представляет собой один из цифровых алгоритмов (MUSIC,BeamScan,MVDR), определяющий направление на цель в азимутальной и угломестной плоскости;
Схема работы модели представлена на рисунке ниже.

1.2 Параметры РЛС
Ключевым параметров РЛС является значение несущей частоты. При выборе частотного диапазона необходимо учитывать задачи РЛС поскольку чем выше значение несущей частоты тем выше затухание и ослабление при распространении, при этому растет чувствительность к малоразмерных объектам.
При решении задачи пеленгации малоразмерных объектов целесообразно использовать диапазоны соизмеримые с размерами целей (20-30 см): L- или S-диапазоны. Выберем L-диапазон (1-2 ГГц) радиолокационной системы, и несущее значение частоты примем равной 1 ГГц (λ ≈ 0,3 м).
# Фундаментальные параметры РЛС
fc = 1e9 # [Гц], несущая частота
c = physconst("LightSpeed") # [м/с], скорость распространения сигнала
lambda = c / fc; # длина волны, м
println("Несущая частота: $(fc*1e-9) ГГц")
println("Скорость света: $c м/с")
println("Длина волны: $lambda м")
1.3 Параметры зондируемого сигнала
В качестве зондирующего импульса выберем сигнала с внутриимпульсной модуляцией - импульсный ЛЧМ-сигнал, с полосой bw - 100 МГц, длительность T - 10 мкс и скважностью Q - 2. Такой сигнал является широкополосный поскольку отношение полосы сигнала к несущей частоте больше 5% (в примере - 10%). Формирование сигнала реализуется с помощью системного объекта EngeePhased.LinearFMWaveform . С помощью параметра n_pulses , n_groups_imps можно задавать количество импульсов накопления в пачке и количество пачек соотвественно.
Важно, для широкополосной системы необходимо задание количества частотных поддиапазонов NBand, в пределах которых сигнал будет считаться узкополосным. Общая полоса определяется частотой дискретизации Ширина поддиапазона вычисляется как (1)
тогда центральные частоты поддиапазонов рассчитываются по формуле (2)
# --- Параметры сигнала ---
bw = 100e6; # Полоса ЛЧМ, Гц
T = 10e-6; # Длительность импульса, с
Q = 2 # скважность импульса
prf = round(1 / (Q * T);sigdigits = 15); # Частота повторения импульсов (округление при необходимости)
type_dir = "Up" # направление изменения частоты ["Up","Down"]
type_sweep = "Positive" # тип изменения частоты ["Positive","Symmetric"]
n_pulses = 2 # количество импульсов накопления
n_groups_imps = 2 # количество пачек импульсов
fs = 2*bw; # Частота дискретизации, Гц
NBand = 64; # количество частотных поддиапазонов
T_coh = n_pulses/prf # время корегентного накопления
ΔR = c/(2*bw) # разрешение по дальности, м
ΔV = lambda/(2*T_coh) # разрешение по скорости, м/с
Rmin = c*T/2 # минимальная одназначная дальность, м
Rmax = c/(2*prf) # максимальная одназначная дальность, м
Vmax = lambda*prf/4 # максимальная одназначная скорость, м/с
# формирование систмного объекта генератора ЛЧМ-сигнала
waveform = EngeePhased.LinearFMWaveform(
SampleRate = fs, # Частота дискретизации, Гц
PulseWidth = T, # Длительность импульса
PRF = prf, # Частота повторения импульсов
SweepBandwidth = bw, # Полоса ЛЧМ
SweepDirection = "Up", # направление изменения частоты
SweepInterval = "Positive", # тип пилы, нессиметричный или симметричный
OutputFormat = "Pulses", # тип формата выхода []
NumPulses = n_pulses # количество импульсов
);
Визуализируем осциллограмму и спектрограмму сигнала для заданного количество импульсов в 1 пачке.
# расчет ЛЧМ-сигнала
sig_lfm = waveform();
# построение осциллограммы и спектрограммы
plot_sig_and_spec(sig_lfm;fs=fs,name_sig = "ЛЧМ");
Тактико-технические характеристики сформированного ЗС приведены ниже:
println("Частота дискретизации: $(round(fs*1e-6;sigdigits=5)) МГц")
println("Длительность импульса: $(round(T*1e6;sigdigits=5)) мкс")
println("Период следования: $(round(1e6/prf;sigdigits=5)) мкс")
println("Частота следования: $(round(prf*1e-3;sigdigits=5)) кГц")
println("Скважность: $(round(Q;sigdigits=5))")
println("Время когерентного накопления: $(round(T_coh*1e3;sigdigits=5)) мс")
println("Количество импульсов в пачке: $n_pulses")
println("Количество пачек: $n_groups_imps")
println("Разрешение по дальности : $(round(ΔR;sigdigits=3)) м")
println("Разрешение по скорости : $(round(ΔV;sigdigits=3)) м/с")
println("Максимальная одназначная дальность : $(round(Rmax;sigdigits=3)) м")
println("Минимальная одназначная дальность : $(round(Rmin;sigdigits=3)) м")
println("Максимальная одназначная скорость: $(round(Vmax;sigdigits=3)) м/с")
1.4 Формирование антенны
Для пеленгации по двум направлениям - азимут и угол места требуется двухмерная антенная решетка. В примере используется наиболее распространенный тип - прямоугольная эквидистантная ФАР, задаваемая с помощью СО EngeePhased.URA. В качестве элемента ФАР выступает изотропный излучатель - EngeePhased.IsotropicAntennaElement.
# параметры антенны
d = lambda /2 # расстояние между элементами
freq_range = [fc-bw fc+bw] # частотный диапазон излучения
Nx = 24 # Число элементов по оси X (разрешение по азимуту)
Ny = 24 # Число элементов по оси Y (разрешение по углу места)
Δaz = rad2deg(lambda/(d*(Nx-1))) # разрешение по азимуту, град
Δel = rad2deg(lambda/(d*(Ny-1))) # разрешение по углу места, град
# Общее число элементов N = Nx * Ny
# СО изотропного антенного элемента
ant_elem = EngeePhased.IsotropicAntennaElement(
FrequencyRange = freq_range, # частотный диапазон
BackBaffled = true # удаление обратного рассеяния ДН антенны
)
# создание прямоугольной ФАР
array = EngeePhased.URA(
Element = ant_elem, # элемент антенны
Size = [Nx Ny], # размер решетки
ElementSpacing = [d d], # расстояние между элементами
);
Визуализируем диаграмму направленности ФАР с помощью функции pattern.
pattern(array,fc) |> display;
Параметры и характеристики прямоугольной ФАР приведены ниже:
println("Расстояние между элементами : $(round(d;sigdigits=3)) м")
println("Размеры прямоугольной ФАР: $(Nx)x$Ny [Nx, Ny]")
println("Частотный диапазон излучения: $(freq_range*1e-9) ГГц")
println("Разрешение по азимуту: $(round(Δaz;sigdigits=3)) °")
println("Разрешение по углу места: $(round(Δel;sigdigits=3)) °")
1.5 Приемо-передающий тракт
Следующий компонент РЛС - приемо-передающий тракт. Он выполняет 2 основные задачи:
- Излучение зондирующего сигнала в пространство в виде электромагнитной волны (ЭМВ);
- Прием отраженного сигнала от цели.
Передающий тракт включает 2 компонента: передатчик, реализующий предварительное усиление и ФАР. Аналогично сформирован и приемный тракт, включающий приемную ФАР и предусилитель.
# Параметры приемо-передающего тракта
Pt = 1e3 # Вт, пиковая мощность
Gtx = 20 # дБ, усиление передающей антенны
LFtx = 1.5 # дБ, потери в передающем тракте
Grx = 20 # дБ, усиление приемной антенны
NF = 2 # дБ, коэффициент шума приемного тракта
LFrx = 1.5 # дБ, потери в приемном тракте
TN = 290 # шумовая температура
NP = physconst("boltzmann") * bw * TN; # Вт, мощность шума
# создания СО передатчик
transmitter = EngeePhased.Transmitter(
Gain = Gtx, # усиление передатчика
PeakPower = Pt, # Пиковая мощность
LossFactor=LFtx # потери в передающем тракте
)
# создания СО широкополосный передающей ФАР
radiator = EngeePhased.WidebandRadiator(
Sensor = array, # ДН антенного элемента
SampleRate = fs, # частота дискретизации
CarrierFrequency = fc, # несущая частота
NumSubbands = NBand # количество частотных поддиапазонов
)
# создания СО широкополосный приемной ФАР
collector = EngeePhased.WidebandCollector(
Sensor = array, # ДН антенного элемента
SampleRate = fs, # частота дискретизации
CarrierFrequency = fc, # центральная частота
NumSubbands = NBand # количество частотных поддиапазонов
)
# создания СО приемного предусилителя (МШУ)
receiver = EngeePhased.ReceiverPreamp(
NoiseMethod="Noise power", # способ задания шума ["Noise power", "Noise temperature"]
Gain = Grx, # усиление приемника, дБ
NoisePower = NP, # мощность шума
LossFactor=LFrx, # потери в приемном тракте
NoiseFigure = NF, # коэффициент шума в тракте
SampleRate = fs # частота дискретизации
);
Параметры и характеристики приемопередающего тракта приведены ниже:
println("-----------------Передатчик-----------------------------")
println("Мощность передатчика: $(round(Pt*1e-3;sigdigits=5)) кВт")
println("Усиление передатчика: $(round(Gtx;sigdigits=5)) дБ")
println("Потери в передающем тракте: $(round(LFtx;sigdigits=5)) дБ")
println("--------------------------------------------------------")
println("-----------------Приемник-----------------------------")
println("Усиление приемника: $(round(Gtx;sigdigits=5)) дБ")
println("Потери в приемном тракте: $(round(LFtx;sigdigits=5)) дБ")
println("Шумовая температура: $(round(LFrx;sigdigits=5)) дБ")
println("Коэффициент шума: $(round(NF;sigdigits=5)) дБ")
println("Шумовая температура: $(round(TN;sigdigits=5)) °К")
println("Мощность шума: $(round(NP;sigdigits=5)) Вт")
1.6 Блок цифровой обработки (пеленгатор)
Последний блок РЛС - алгоритм пеленгации. Такой блок решает задачу определения направления прихода сигнала и, как следствие, направление целевого объекта по отношению к РЛС. Рассмотрим наиболее используемые: "BeamScan", "MVDR", "MUSIC"
az_scan = -20:0.5:20 # диапазон сканирования по азимуту
el_scan = 0:0.5:30 # диапазон сканирования по углу места
n_sig = 4 # количество оценок (количество целей)
algorithm = "MUSIC" # @param ["Beamscan","MVDR","MUSIC"]
beamscan = EngeePhased.BeamscanEstimator2D(
SensorArray = array, # антенна
OperatingFrequency = fc, # несущая частота
AzimuthScanAngles = az_scan, # диапазон сканирования по азимуту
ElevationScanAngles = el_scan, # диапазон сканирования по углу места
DOAOutputPort = true, # включение выхода оценок пеленга
NumSignals = n_sig # количество оценок
);
mvdr = EngeePhased.MVDREstimator2D(
SensorArray = array, # антенна
OperatingFrequency = fc, # несущая частота
AzimuthScanAngles = az_scan, # диапазон сканирования по азимуту
ElevationScanAngles = el_scan, # диапазон сканирования по углу места
DOAOutputPort = true, # включение выхода оценок пеленга
NumSignals = n_sig # количество оценок
);
music = EngeePhased.MUSICEstimator2D(
SensorArray = array, # антенна
OperatingFrequency = fc, # несущая частота
AzimuthScanAngles = az_scan, # диапазон сканирования по азимуту
ElevationScanAngles = el_scan, # диапазон сканирования по углу места
DOAOutputPort = true, # включение выхода оценок пеленга
NumSignalsSource="Property",
NumSignals = n_sig, # количество оценок
ForwardBackwardAveraging=false # учет прямого и обратного усреднения
);
doa_alg = if algorithm == "Beamscan"
println("Алгоритм: Beamscan")
deepcopy(beamscan)
elseif algorithm == "MVDR"
println("Алгоритм: MVDR");
deepcopy(mvdr)
elseif algorithm == "MUSIC"
println("Алгоритм: MUSIC");
deepcopy(music)
end;
1.7 Фоноцелевая обстановка
Блок фоноцелевой обстановки включается 2 компонента: среду распространения и целевые объекты. Для моделирования среды с учетом широкополосности сигнала используется СО EngeePhased.WidebandLOSChannel, он позволяет учесть частотной дисперсии при распространении объекта и влияния погодных условий (гидрометеоры, град, изменение давления и т.п.).
# параметры среды распространения
is_atmosphere = true # учет атмоферного влияния на распространения ЭМВ
T_ch = 15 # температура атмосферы,°С
press_air = 101325.0 # давление атмосферы, ПА
vap_density = 7.5 # плотность влажного пара
liq_density = 0.0 # плотность конденсата
rain_rate = 0.0 # интенсивность осадков
channel_fwd = EngeePhased.WidebandLOSChannel(
PropagationSpeed = c,
SampleRate = fs,
OperatingFrequency = fc,
NumSubbands = NBand,
TwoWayPropagation = false, # учет двунаправленного распространения сигнала
SpecifyAtmosphere = is_atmosphere,
Temperature = T_ch,
DryAirPressure = press_air,
WaterVapourDensity = vap_density,
LiquidWaterDensity = liq_density,
RainRate = rain_rate
);
channel_bwd = deepcopy(channel_fwd); # канал распространения от целей до РЛС (идентичен)
Для формирования цели используется СО EngeePhased.WidebandBackscatterRadarTarget. Он позволяет описать пользовательскую ДОР цели с заданной точностью сетки пространственного разрешения (азимут и углу места). В примере подгружается матрица ДОР (файл quad_rcs_pattern.jld2) для аналитически рассчитанной модели БПЛА типа квадрокоптер. Более подробное описание расчета ДОР рассматривается в [2,3].
data = JLD2.load("quad_rcs_pattern.jld2")
num_tgts = n_sig # количество целей (зададим равное количество оценок пеленгатора)
rcs_az = data["az"] # диапазон учета ЭПР по азимуту (-180:180), град
rcs_el = data["el"] # диапазон учета ЭПР по углу места (-90:90), град
rcs_pattern = data["rcs"] # матрица ДОР цели [181x361], м^2
type_tgt = "Nonfluctuating" # тип цели ["Nonfluctuating","Swerling1","Swerling2","Swerling3","Swerling4"]
# формирование СО ЭПР широкополосных целей
rcs_targets = EngeePhased.WidebandBackscatterRadarTarget(
AzimuthAngles = rcs_az, # диапазон учета ЭПР по азимуту
ElevationAngles = rcs_el, # диапазон учета ЭПР по углу места
RCSPattern = rcs_pattern, # ДОР цели
PropagationSpeed = c, # скорость распространения, м/с
OperatingFrequency = fc, # несущая частота, Гц
NumSubbands = NBand, # количество частотных поддиапазонов
Model = type_tgt # тип флуктуации ЭПР цели
);
Для визуализации объекта БПЛА и его ДОР воспользуемся функцией calc_plot_rcs:
text_title = "Визуализация и ДОР цели"
k_norm = 0.6 # коэффициент масштабирования размеров ДОР
calc_plot_rcs(rcs_pattern,rcs_az,rcs_el;title = text_title,k = k_norm)
Общий вид ДОР БПЛА имеет изрезанный, лепестковый характер из-за геометрических особенности цели - наличия неоднородностей в виде винтов . Можно заметить: максимальное значение ЭПР достигается при облучении БПЛА сверху и снизу, что обусловлено максимальной облучаемой площадью, в горизонтальной плоскости ЭПР уменьшается.
1.8 Сценарий движения
Последним этапом формирования модели является задания траектории движения цели и РЛС. В сценарии считается, что РЛС расположена в начале системы координат и стационарна. Начальное положение и скорости целей задаются в переменных targets_pos и targets_vel, модель движения выбрана равномерная
# Модель движения цели
targets_pos = [
[1000; 100; 200];; # 1 цель
[1200; -10; 200];; # 2 цель
[1100; -150; 150];; # 3 цель
[950; 50; 320];; # 4 цель
] # начальное положение цели, м
targets_vel = [
[-130; 0; 0];; # скорость 1-ой цели (на РЛС параллельно оси x)
[-150; 0; 0];; # скорость 2-ой цели (на РЛС параллельно оси x)
[-200; 0; 0];; # скорость 3-ой цели (на РЛС параллельно оси x)
[-230; 0; 0];; # скорость 4-ой цели (на РЛС параллельно оси x)
] # скорость цели, м/с
type_traj = "Velocity" # модель движения - равномерная ["Velocity","Acceleration","Custom"]
targets_platform = EngeePhased.Platform(
MotionModel = type_traj,
InitialPosition = targets_pos, # начальное положение
Velocity = targets_vel # скорость цели
)
# Модель движения РЛС
radar_pos = [0; 0; 0] # Положение РЛС
radar_vel = [0; 0; 0] # скорость РЛС
radar_platform = EngeePhased.Platform(
InitialPosition = radar_pos, # положение РЛС
Velocity = radar_vel # скорость РЛС
);
abs_dist,ang_true = rangeangle(targets_pos,radar_pos)
for i in 1:num_tgts
println("----------------Цель № $(i)----------------------------")
println("азимут: $(round(ang_true[1,i];sigdigits=3))°, угол места: $(round(ang_true[2,i];sigdigits=3))°")
println("Дальность от РЛС до цели №$(i): $(round(abs_dist[i]*1e-3;sigdigits=5)) км")
end
С помощью функции viz_scene отобразим начальное состояния сценария с учетом направлений визирования целей:
viz_scene(targets_pos, radar_pos) |> display;
В результате, были сформированы следующие тактико-технические характеристики РЛС:
println("-----------------РЛС-----------------------------")
println("Несущая частота: $(round(fc*1e-9;sigdigits=3)) ГГц")
println("Скорость света: $c м/с")
println("Длина волны: $(round(lambda;sigdigits=3)) м")
println("Количество частотных поддиапазонов: $NBand м")
println("Положение РЛС: $radar_pos [Px,Py,Pz], м")
println("скорость РЛС: $radar_vel [Vx,Vyc,Vz],м/с")
println("--------------------------------------------------------")
println("-----------------Генератор-----------------------------")
println("Частота дискретизации: $(round(fs*1e-6;sigdigits=5)) МГц")
println("Длительность импульса: $(round(T*1e6;sigdigits=5)) мкс")
println("Период следования: $(round(1e6/prf;sigdigits=5)) мкс")
println("Частота следования: $(round(prf*1e-3;sigdigits=5)) кГц")
println("Скважность: $(round(Q;sigdigits=5))")
println("Время когерентного накопления: $(round(T_coh*1e3;sigdigits=5)) мс")
println("Количество импульсов в пачке: $n_pulses")
println("Количество пачек: $n_groups_imps")
println("Разрешение по дальности : $(round(ΔR;sigdigits=3)) м")
println("Разрешение по скорости : $(round(ΔV;sigdigits=3)) м/с")
println("Максимальная одназначная дальность : $(round(Rmax;sigdigits=3)) м")
println("Минимальная одназначная дальность : $(round(Rmin;sigdigits=3)) м")
println("Максимальная одназначная скорость: $(round(Vmax;sigdigits=3)) м/с")
println("--------------------------------------------------------")
println("-----------------Передатчик-----------------------------")
println("Мощность передатчика: $(round(Pt*1e-3;sigdigits=5)) кВт")
println("Усиление передатчика: $(round(Gtx;sigdigits=5)) дБ")
println("Потери в передающем тракте: $(round(LFtx;sigdigits=5)) дБ")
println("--------------------------------------------------------")
println("-----------------Приемник-----------------------------")
println("Усиление приемника: $(round(Gtx;sigdigits=5)) дБ")
println("Потери в приемном тракте: $(round(LFtx;sigdigits=5)) дБ")
println("Шумовая температура: $(round(LFrx;sigdigits=5)) дБ")
println("Коэффициент шума: $(round(NF;sigdigits=5)) дБ")
println("Шумовая температура: $(round(TN;sigdigits=5)) °К")
println("Мощность шума: $(round(NP;sigdigits=5)) Вт")
println("--------------------------------------------------------")
println("-----------------Антенна-------------------------------")
println("Расстояние между элементами : $(round(d;sigdigits=3)) м")
println("Размеры прямоугольной ФАР: $(Nx)x$Ny [Nx, Ny]")
println("Частотный диапазон излучения: $(freq_range*1e-9) ГГц")
println("Разрешение по азимуту: $(round(Δaz;sigdigits=3)) °")
println("Разрешение по углу места: $(round(Δel;sigdigits=3)) °")
println("--------------------------------------------------------")
println("-----------------Пеленгатор-----------------------------")
println("Диапазон сканирования по азимуту: $(az_scan) град.")
println("Диапазон сканирования по углу места: $(el_scan) град.")
println("Количество оценок: $n_sig")
println("--------------------------------------------------------")
2. Расчет сценария работы РЛС
2.1 Прогон модели
С применением ранее инициализированных СО сформируем сценарий работы РЛС. После выполнения расчета будут доступны следующие переменные:
-
data_y- отраженный сигнал на выходе приемного устройства -
data_ang_true- истинные значения пеленга для N циклов (пачек импульсов) -
data_ang_est- оценки пеленга для N циклов (пачек импульсов) -
data_spec- оценка спектра мощности на выходе алгоритма пеленгации для выбранного алгоритма
Важно отметить, что положение целей обновляется от импульса к импульса. Поэтому, в цикле работы системы на вход передатчика необходимо подавать 1 импульс ЛЧМ-сигнала. Для этого используются вспомогательные переменные sig_lfm_one_imps и sig_buffer_one_groups.
# Выделение памяти для выходного сигнала и оценок пеленга
data_y = zeros(ComplexF64,length(sig_lfm),Nx*Ny,n_groups_imps)
data_ang_true = zeros(2, num_tgts, n_groups_imps); # истинные направления на цели
data_ang_est = zeros(2, num_tgts, n_groups_imps); # истинные направления на цели
data_spec = zeros(length(el_scan),length(az_scan), n_groups_imps); # данные спектра от углов визирования на выходе пеленгатора
len_one_imps = round(Int,length(sig_lfm)/n_pulses) # количество отсчетов 1 импульса
sig_lfm_one_imps = sig_lfm[1:len_one_imps] # выделение 1 импульса
sig_buffer_one_groups = similar(data_y[:,:,1]) # буффер для накопления пачки импульсов
println("Расчет сценария запущен...")
@inbounds for i in 1:n_groups_imps
println("Расчет пачки №$(i)...")
@inbounds for k in 1:n_pulses
println("Расчет импульса №$(k)...")
# Обновление положения РЛС и целей с шагом периода импульса 1/prf
radar_pos, radar_vel = radar_platform(1/prf)
tgts_pos, tgts_vel = targets_platform(1/prf)
# Расчет и сохранение пеленга целей
_,ang_doa = rangeangle(tgts_pos,radar_pos)
data_ang_true[:,:,i] .= deepcopy(ang_doa)
# Излучение широкополосного сигнала в среду
sig_tx = transmitter(sig_lfm_one_imps) # прохождение передающего тракта
sig_rad = radiator(sig_tx,ang_doa) # излучение сигнала в пространство
# Распространение от РЛС до целей
sig_fwd = channel_fwd(sig_rad,radar_pos,tgts_pos,radar_vel,tgts_vel)
# Отражение от целей
sig_refl = rcs_targets(sig_fwd,ang_doa)
# Распространение от целей к РЛС
sig_ret = channel_bwd(sig_refl,tgts_pos,radar_pos,tgts_vel,radar_vel)
# приемный тракт
y = collector(sig_ret,ang_doa) # приемная ФАР
sig_buffer_one_groups[(1:len_one_imps).+(k-1)*len_one_imps,:] .= receiver(y) # предварительное усиление с учетом теплового шума приемника (МШУ)
end
# расчет пеленга
data_spec[:,:,i],data_ang_est[:,:,i] = doa_alg(sig_buffer_one_groups); # оценка спектра и направлений прихода целей
data_y[:,:,i] .= sig_buffer_one_groups; # сохранение принятого сигнала
end
println("Расчет сценария завершен!")
2.2 Анализ точности системы
Построим график спектра мощности сигнала на выходе алгоритма пеленгации от азимуту и угла места. Сетка угол наследуется из пункта 1.6, где задавался алгоритм пеленгации и точность сетки.
n_imp = 1 # номер анализируемого импульса
plot_size = (600,600) # разрешение
plot_spec_azel(az_scan,el_scan,data_spec[:,:,n_imp];
title = "Cпектр мощности на выходе пеленгатора (алгоритм $(algorithm))",
size = plot_size
)
На графике видны 4 характерных пика, сигнализирующие о наличии целей в данных направлениях. Ниже представлено сравнение истинных значений направлений и оценок пеленга для соответствующий целей.
На выходе пеленгатора оценки выдаются в порядке детектирования, поэтому предварительно выполним сортировку по критерию минимума Евклидова расстояния (L2-норма) оценок с помощью функции sort_ang_mse:
sort_ang_est = sort_ang_mse(data_ang_true,data_ang_est) # сортированные данные оценок углов
for k in 1:n_groups_imps
println("--------------- Пачка импульсов №$(k) -----------------------------")
for i in 1:num_tgts
println(" ---------------- Цель №$(i) ----------------------------")
println(" Цель: азимут = $(round(data_ang_true[1,i,k];sigdigits=4)) °, угол места = $(round(data_ang_true[2,i,k];sigdigits=4))°");
println(" Оценка: азимут = $(round(sort_ang_est[1,i,k];sigdigits=4))°, угол места = $(round(sort_ang_est[2,i,k];sigdigits=4))°");
println(" Ошибка: азимут = $(round(abs(sort_ang_est[1,i,k] - data_ang_true[1,i,k]);sigdigits=4))°, угол места = $(round(abs(sort_ang_est[2,i,k] - data_ang_true[2,i,k]);sigdigits=4))°");
end
println("Максимальная погрешность: $( round(maximum(abs.(data_ang_true[:,:,k].-sort_ang_est[:,:,k])); sigdigits=4) )°")
println("-------------------------------------------------------------")
end
Максимальная погрешность пеленга составила 0.538°, что сопоставимо с ценой деления сетки (0.5°). Как следствие используемый алгоритм применим для решаемой задачи. Для увеличения точности можно увеличить количество элементов ФАР и уменьшить анализируемую сетку углов
3. Анализ и сравнение алгоритмов пеленгации
В заключительном разделе исследуем точность алгоритмов пеленгации (Beamscan,MVDR,Music) при различных отношении сигнал-шум.
Для этого оценим текущее ОСШ с помощью основного уравнения радиолокации radareqsnr. Автоматически вычислить ОСШ и сгенерировать jl-файл расчета можно с применением приложения "Расчет уравнений радиолокации". Ниже представлен пример использования приложения для расчета ОСШ для дальности 1.5 км, длине 0.3 м, длительность импульса 10 мкс и пиковой мощности 1 кВт.

В результате, для текущей конфигурации значение ОСШ составило 34.1 дБ. После генерации в файловом браузере появился файла расчета RadarEquationScript_2026-05-29_21:46:25.749.jl. Результирующий код с параметризацией представлен ниже:
tgtrng = maximum(abs_dist) # дальность до цели
# расчет ОСШ с точностью до 3 знаков
snr = round(radareqsnr(
lambda, # длина волны
tgtrng, # дальность до цели
Pt, # импульсная мощность
T, # длительность импульса
RCS = 0.01, # ЭПР цели (выберем среднее значение)
Gain = Grx, # усиление антенны
Loss = LFtx+LFrx, # потери в тракте
Ts = TN # шумовая температура
)[1];sigdigits = 3)
println("Текущее ОСШ: $snr дб")
Для анализа эффективности алгоритма выделим последний импульс сигнала на выходе приемника и сформируем вектор увеличения мощности шума noise_vec и значений ОСШ snr_vec:
y = deepcopy(data_y[:,:,end]) # копируем данные сигнала для последнего импульса
ang_trues_end = data_ang_true[:,:,end]; # копируем истинные направления целей
noise_vec = collect(20:10:70) # вектор увеличения мощности шума от 20 до 70 дБ относительно NP
snr_vec = snr .- noise_vec # вектор значений ОСШ
println("Вектор ОСШ: $(round.(snr_vec;sigdigits=3)) дБ")
3.1 Алгоритм MUSIC
Рассмотрим эффективность работы алгоритма MUSIC. Выполним 6 вызовов СО MUSIC с уменьшением ОСШ в 10 раз от вызова к вызову:
# выделением памяти под простравнственный спектр и оценки углов для MUSIC
spec_music = zeros(length(el_scan),length(az_scan), length(noise_vec))
ang_music = zeros(2, num_tgts, length(noise_vec)) # истинные направления на цели
@inbounds for i in 1:length(noise_vec)
println("Расчет оценки для ОСШ $(round(snr_vec[i];sigdigits=3)) дБ")
release!(music) # сбрасываем системный объекта к начальному состоянию
noise_sig = sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # новые даные для шумового сигнала
spec_music[:,:,i],ang_music[:,:,i] = music(y.+noise_sig); # оценка спектра и направлений прихода целей
end
Сопоставим результаты пеленгации для заданных ОСШ:
fig_spec_music_vec = [] #
@inbounds for i in 1:length(noise_vec)
push!(fig_spec_music_vec,
plot_spec_azel(az_scan,el_scan,spec_music[:,:,i];
title = "MUSIC (ОСШ = $(round(snr_vec[i];sigdigits=3)) дб)"
)
)
end
plot(fig_spec_music_vec...)
При 24 дБ наблюдается уверенный прием всех 4 целей. По мере уменьшения ОСШ количество детектируемых целей падает: для -15.9 дБ видно 2 цели, для -25.9 дБ детектирование целей теряется .
Далее, рассмотрим максимальную ошибку позиционирования для всех целей:
sort_ang_music = similar(ang_music) # сортированные данные оценок углов
[sort_ang_music[:,:,i] .= sort_ang_mse(ang_trues_end,ang_music[:,:,i]) for i in 1:length(snr_vec)]
# построение таблицы
calc_table_ang_analysis(ang_trues_end,sort_ang_music,snr_vec)
При уменьшении ОСШ до -16 дБ точность позиционирования падает для слабых целей, по мере дальнейшего ослабления сигнала ошибка приобретает случайный характер.
3.2 Алгоритм MVDR
Далее, рассмотрим эффективность работы алгоритма MVDR.
# выделением памяти под простравнственный спектр и оценки углов для MVDR
spec_mvdr = zeros(length(el_scan),length(az_scan), length(noise_vec))
ang_mvdr = zeros(2, num_tgts, length(noise_vec)) # истинные направления на цели
@inbounds for i in 1:length(noise_vec)
println("Расчет оценки для ОСШ $(round(snr_vec[i];sigdigits=3)) дБ")
release!(mvdr) # сбрасываем системный объекта к начальному состоянию
noise_sig = sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # новые даные для шумового сигнала
spec_mvdr[:,:,i],ang_mvdr[:,:,i] = mvdr(y.+noise_sig); # оценка спектра и направлений прихода целей
end
Сопоставим результаты пеленгации для заданных ОСШ:
fig_spec_mvdr_vec = [] #
@inbounds for i in 1:length(noise_vec)
push!(fig_spec_mvdr_vec,
plot_spec_azel(az_scan,el_scan,spec_mvdr[:,:,i];
title = "MVDR (ОСШ = $(round(snr_vec[i];sigdigits=3)) дб)"
)
)
end
plot(fig_spec_mvdr_vec...)
При 14 дБ наблюдается уверенный прием всех 4 целей. По мере уменьшения ОСШ количество детектируемых целей падает: для -16 дБ видно 2 цели, для -20 дБ уже не видно ни одной.
Далее, рассмотрим максимальную ошибку позиционирования для всех целей:
sort_ang_mvdr = similar(ang_mvdr) # сортированные данные оценок углов
[sort_ang_mvdr[:,:,i] .= sort_ang_mse(ang_trues_end,ang_mvdr[:,:,i]) for i in 1:length(snr_vec)]
# построение таблицы
calc_table_ang_analysis(ang_trues_end,sort_ang_mvdr,snr_vec)
При уменьшении ОСШ до 0 дБ точность позиционирования падает для слабых целей, по мере дальнейшего ослабления сигнала ошибка приобретает случайный характер
3.3 Алгоритм Beamscan
Далее, рассмотрим эффективность работы алгоритма BeamScan.
# выделением памяти под простравнственный спектр и оценки углов для beamscan
spec_beamscan = zeros(length(el_scan),length(az_scan), length(noise_vec))
ang_beamscan = zeros(2, num_tgts, length(noise_vec)) # истинные направления на цели
@inbounds for i in 1:length(noise_vec)
println("Расчет оценки для ОСШ $(round(snr_vec[i];sigdigits=3)) дБ")
release!(beamscan) # сбрасываем системный объекта к начальному состоянию
noise_sig = sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # новые даные для шумового сигнала
spec_beamscan[:,:,i],ang_beamscan[:,:,i] = beamscan(y.+noise_sig); # оценка спектра и направлений прихода целей
end
Сопоставим результаты пеленгации для заданных ОСШ:
fig_spec_beamscan_vec = [] #
@inbounds for i in 1:length(noise_vec)
push!(fig_spec_beamscan_vec,
plot_spec_azel(az_scan,el_scan,spec_beamscan[:,:,i];
title = "beamscan (ОСШ = $(round(snr_vec[i];sigdigits=3)) дб)"
)
)
end
plot(fig_spec_beamscan_vec...)
Для любого значения ОСШ (от -20 до 30 дБ) наблюдает устойчивое детектирование целей. При этом, ширина главных лепесткой значительно больше, чем для предыдущих алгоритм, что говорит о плохой разрешающей способности и невозможности детектирования 2 близких целей.
Далее, рассмотрим максимальную ошибку позиционирования для всех целей:
sort_ang_beamscan = similar(ang_beamscan) # сортированные данные оценок углов
[sort_ang_beamscan[:,:,i] .= sort_ang_mse(ang_trues_end,ang_beamscan[:,:,i]) for i in 1:length(snr_vec)]
# построение таблицы
calc_table_ang_analysis(ang_trues_end,sort_ang_beamscan,snr_vec)
Как мы видим, даже при малом ОСШ цели верно детектируются, что говорит о высокой помехоустойчивости алгоритма
3.4 Сравнение алгоритмов
Выполним сравнение результатов для выбранного значения ОСШ snr_i :
plot_spec_2d = [];
for snr_i in eachindex(snr_vec)
spec_beamscan_snr_0db_norm = maximum(spec_beamscan[:,:,snr_i];dims=1)[:]./maximum(spec_beamscan[:,:,snr_i])
spec_mvdr_snr_0db_norm = maximum(spec_mvdr[:,:,snr_i];dims=1)[:]./maximum(spec_mvdr[:,:,snr_i])
spec_music_snr_0db_norm = maximum(spec_music[:,:,snr_i];dims=1)[:]./maximum(spec_music[:,:,snr_i])
fig = plot(az_scan,spec_beamscan_snr_0db_norm,lab="BeamScan",xlab="Азимут, град",
ylab="Нормированная мощность",lw=2,
title="Сравнение алгоритмов пеленгации для ОСШ = $(round(snr_vec[snr_i];sigdigits = 3)) дБ"
)
plot!(az_scan,spec_mvdr_snr_0db_norm,lab="MVDR",lw=2)
plot!(az_scan,spec_music_snr_0db_norm,lab="MUSIC",lw=2)
push!(plot_spec_2d,fig)
end
plot(plot_spec_2d...;layout = (length(plot_spec_2d),1),size = (600,300*length(plot_spec_2d)))
Как было ранее замечено, алгоритм BeamScan имеет худшее разрешение, но лучшее помехоустойчивость: при -26 дБ алгоритм смог детектировать цель по направлению 6.5 градусов. Алгоритм MVDR имеет отличное разрешение при достаточном ОСШ (от 10 дБ), но при этом сильно подвержен шумовому влияния. Наилучшими характеристиками по разрешению и помехоустойчивости обладает алгоритм MUSIC: даже при ОСШ = -6 дБ на пространственном спектре наблюдается уверенное детектирование 4 целей***.***
Заключение
Таким образом, в примере был рассмотрен подход моделирование многоканальной широкополосной системы пеленгации с множественными целями в среде моделирования Engee.
Также, было произведен сравнительный анализ алгоритмов пеленгации (MUSIC,BeamScan,MVDR) при различных значениях отношения сигнал шум (в диапазоне от -36 до 14 дБ):
- алгоритм MVDR : показан высокую эффективность - высокую разрешающую способность при значениях ОСШ больше нуля, но сильнее всех подвержен шумовому воздействию;
- алгоритм Beamscan: при досточном ОСШ имеет наибольшую ширину главного лепестка и как следствие, наихудшую разрешающую способность пространственного спектра. При этом, обладает наилучшую помехоустойчивость и способен обнаружить цель даже при ОСШ< -20 дБ;
- алгоритм Music: показан высокую эффективность в большом динамическом диапазоне (от -20 до 14 дБ) и обладает наилучшей разрешающей способностью среди всех представленных алгоритмов.
Источники
- Ван Трис, Х. Оптимальная обработка сигналов. Нью-Йорк: Wiley-Interscience, 2002.
- Семичастнов, А. Е., и Д. А. Балакин. «Моделирование сценариев работы радиолокационной системы для классификации беспилотных летательных аппаратов и птиц на основе микродоплеровских сигнатур в среде Engee». Электронные библиотеки, т. 28, вып. 4, ноябрь 2025 г., сс. 943-52, doi:10.26907/1562-5419-2025-28-4-943-952.
- A. R. Gorbunov, M. R. Akmalov, A. E. Semichastnov and D. A. Balakin, "Algorithm for Drone Recognition from Birds Based on a Neural Network Approach Using Radar Data," 2025 7th International Youth Conference on Radio Electronics, Electrical and Power Engineering (REEPE), Moscow, Russian Federation, 2025, pp. 1-8, doi: 10.1109/REEPE63962.2025.10971050.