Частотная адаптация при воздействии помехи
В этом примере показано, как моделировать радилокационные системы с частотной адаптацией при воздействии помех.
Подключение библиотек и вспомогательных функций
let
intalled_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["DSP"]
for pack in list_packages
pack in intalled_packages || Pkg.add(pack)
end
end
using DSP
function TaylorSpectrumWindow(N, nbar = 4, sll = 30)
function calcFm(m, sp2, A, NBAR)
n = Array(1:NBAR-1)
p = [Array(1:m-1); Array((m+1):(NBAR-1))]
Num = prod((1 .- (m^2/sp2)./(A^2 .+(n .-0.5).^2)));
Den = prod((1 .- m.^2 ./p.^2));
return ((-1)^(m+1).*Num)./(2 .*Den)
end
A = acosh((10^(sll/20)))/π
sp2 = nbar^2 / (A^2 + (nbar - 0.5)^2)
Fm = zeros(nbar - 1)
sum = zeros(N)
k = Array(0:N-1)
xi = @. (k-0.5*N+0.5)/N
for m = 1:(nbar-1)
Fm[m] = calcFm(m,sp2,A,nbar);
@. sum = Fm[m]*cos(2*pi*m*xi)+sum;
end
@. sum = 2sum + 1
return sum
end
function calc_spectrogramm(x,t_step;down_lim = -40)
plotlyjs()
out,f,t = EngeePhased.Functions.spectrogram(x;
window = DSP.kaiser(128,10),
noverlap=120,nfft=256,fs=1/t_step,
centered=true,freqloc = "yaxis",out = :data
)
out_sig = map(x -> x < down_lim ? down_lim : x,DSP.pow2db.(abs.(out)))
fig = Plots.heatmap(
t[:].*1e6,f*1e-6,out_sig,color= :jet,
gridalpha=0.5,colorbar_title = "Мощность, дБВт",
xticks=0:50:250,yticks=(-1:0.2:1),
)
Plots.xlabel!("Время, мкс")
Plots.ylabel!("Частота Доплера, МГц")
Plots.title!("Cпектрограмма сигнала")
return fig
end;
Введение
Активные фазированные антенные решетки с электронным сканированием могут помогать решать множество задач, использующих одно и то же оборудование антенной решетки. Их могут использовать в радарах, РЭБ и системах связи. Однако окружающая среда, в которой работают системы такого типа, сложна и может вносить нежелательные помехи. Например, устройство постановки помех на ретрансляторе может повторить принятый сигнал радара и повторно передать его, чтобы сбить радар с толку. Адаптация к частоте может быть эффективным методом противодействия сигналам, генерируемым источниками помех, и способствовать эффективной работе этих систем.
В этом примере моделируется сценарий со стационарным моностатическим радаром и движущейся авиационной целью (похожий сценарий рассматривается в примере Перестройка центральной частоты при воздействии помех). Затем самолет будет генерировать подменяющий сигнал, который сбивает радар с толку. Как только радар обнаружит источник помех, могут быть использованы методы перестойки центральной частоты, позволяющие РЛС обойти помехи.
1. Моделирование системы в среде без помех
Предположим, что в начале координат находится моностатический радар X-диапазона.
fc = 10e9; # несущая частота работы РЛС X-диапазона, Гц
fs = 2e6; # частота дискретизации модели, Гц
c = 3e8; # скорость распросранения сигнала, м/с
lambda = c/fc; # длина волны, м
radar_pos = [0;0;0]; # положение РЛС
radar_vel = [0;0;0]; # скорость РЛС
Приемная антенна состоит из 64 элементов (8х8) эквидистатной прямоугольной антенной решетки с шагом половины длины волны между элементами:
antenna = EngeePhased.URA(
Element=EngeePhased.CosineAntennaElement(), # косинусный элемент АР
Size=[8 8], # размер антенной решетки
ElementSpacing=[lambda/2 lambda/2], # расстояние между элементами
)
Для лучшей характеристики направленности АР, сформируем 2-мерную оконную функцию Тейлора:
taper = TaylorSpectrumWindow(8); # взвешанная фукция окно Тэйлора
taperURA = taper.*taper' # формирование 2д взвешанной функции
antenna.Taper = taperURA # взвешивание АР
# Визуализация оконной функции
surface(taperURA,title="2-мерная оконная фукция Тейлора")
Диаграмму направленности антенной решетки можно построить с помощью функции pattern
:
pattern(antenna,fc,Type="powerdb")
plot!(title="ДН прямоугольной АР",colorbar_title="Мощность, дБВт")
В рассматриваемом РЛС источником сигнала будет являться сигнал с линейной частотной модуляцией (ЛЧМ):
release!(antenna)
# Формирование системных объектов элементов РЛС
bw = 1e5 # полоса сигнала в частотной области
pw = 1e-4 # длительной импульса
wav = EngeePhased.LinearFMWaveform( # генератор ЛЧМ
SampleRate=fs, # частота дискретизации модели
PulseWidth= pw, # длительность импульса
SweepBandwidth= bw, # полоса сигнала в частотной области
PRF=4e3, # частота следования импульсов
FrequencyOffsetSource="Input port" # включение входа для настройки центральной частоты
)
tx = EngeePhased.Transmitter( # передатчик
Gain=20, # усиление, дБ
PeakPower=500 # пиковая мощность, Вт
)
txArray = EngeePhased.WidebandRadiator( # Широкополосная передающая АР
SampleRate=fs, # частота дискретизации модели
Sensor=antenna, # задание геометрии АР
CarrierFrequency=fc # несущая частота
)
rxArray = EngeePhased.WidebandCollector(
SampleRate=fs, # частота дискретизации модели
Sensor=antenna, # задание геометрии АР
CarrierFrequency=fc # несущая частота
)
rxPreamp = EngeePhased.ReceiverPreamp( # предусилитель
Gain=10, # усиление, дБ
NoiseFigure=5, # фактор шума, дБ
SampleRate=fs, # частота дискретизации модели
);
Далее, создадим объекты цели и среды распространения:
target = EngeePhased.RadarTarget(
MeanRCS=100, # среднее значение ЭПР цели
OperatingFrequency=fc # несущая частота
)
target_pos = [8000;1000;1000] # начальное положение цели
target_vel = [100;0;0] # вектор скорости цели
envout = EngeePhased.WidebandFreeSpace(
TwoWayPropagation=false, # Учет двунаправленного распространения
SampleRate = fs, # частота дискретизации
OperatingFrequency=fc, # несущая частота
PropagationSpeed=c # скорость распространения
)
envin = deepcopy(envout); # дублирование канала обратного распространения
С помощью ранее заданных системных объектов системы вычислим отраженный сигнал от цели в условиях отсуствия помехи:
tgtRng, tgtAng = rangeangle(target_pos, radar_pos)
x = deepcopy(wav(0)) # Генератор
xt = tx(x) # передатчик
xtarray = txArray(xt,tgtAng) # передающая антенная решетка
yp = envout(xtarray,radar_pos,target_pos,radar_vel,target_vel) # канал прямого распространения
yr = target(yp) # цель
ye = envin(yr,target_pos,radar_pos,target_vel,radar_vel) # среда распространения
yt = rxArray(ye,tgtAng) # приемная антенная решетка
yt = rxPreamp(yt); # приемник
Вычислим направление и реализуем формирование луча с помощью BeamscanEstimator2D
и SubbandPhaseShiftBeamformer
estimator = EngeePhased.BeamscanEstimator2D(
SensorArray=antenna, # геометрия АР
DOAOutputPort=true, # Включение выхода уголов
OperatingFrequency=fc, # несущая частота
NumSignals=1, # количество сигналов
AzimuthScanAngles=-40:40, # диапазон сканирования по азимуту
ElevationScanAngles=-60:60 # диапазон сканирования по углу места
)
doa = estimator(yt)[2] # расчет оценки пеленга с помощью СО BeamscanEstimator2D
beamformer = EngeePhased.SubbandPhaseShiftBeamformer(
SensorArray=antenna, # геометрия АР
OperatingFrequency=fc, # несущая частота
DirectionSource="Input port", # способ задания углов визирования
SampleRate=fs, # частота дискретизации
WeightsOutputPort=true # весовые коэффициенты
)
ybf = beamformer(yt,doa)[1]; # применения алгоритма формирования СО SubbandPhaseShiftBeamformer
Выполним согласованную фильтрацию, для этого вычислим коэффициенты фильтра по средствам функции getMatchedFilter
mfcoeff1 = getMatchedFilter(wav) # Коэффициентов согласованнго фильтра
mf1 = EngeePhased.MatchedFilter(Coefficients=mfcoeff1);
Применим согласованную фильтрацию к сформированному лучу и визуализируем результат работы РЛС:
y1 = mf1(ybf) # применение согласованного фильтра
nSamples = wav.SampleRate/wav.PRF # количество отсчетов в периоде импульса
t = ((0:nSamples-1).-(length(mfcoeff1)-1))./fs # временная сетка задержек
r = t*c/2 # пересчет из задержки в дальность
# Результат работы РЛС
plot(r/1000,abs.(y1),xlabel="Дальность, км",ylabel="Амплитуда",lab="",xticks=-1e6:5:1e6,title="Выход СФ")
2. Добавление помеховой компоненты
Теперь смоделируем помеховый сигнал по средством имитации полезного сигнала и задержки во времени:
# имитация помехового сигнала путемдублирование генератор сигнала
jwav = deepcopy(wav)
xj = jwav(0) # формирование помеховой составляющей
Npad = ceil(Int64,3500/(c/fs)) # количество элементов задержки (1 отсчет - 150 метров)
xj = circshift(xj,Npad) # сдвиг сигнала на распространения
txjam = EngeePhased.Transmitter(
Gain=10, # усиление помехи
PeakPower=5 # мощность помехи
)
# Расчет отраженного сигнала + помеха
xj = txjam(xj) # излучение помехи
ye = envin(yr.+xj,target_pos,radar_pos,target_vel,radar_vel) # канал прямого распространения
yt = rxArray(ye,tgtAng) # приемная АР
yt = rxPreamp(yt) # предусилитель
ybfj = beamformer(yt,doa)[1] # формирователь луча
y1j = mf1(ybfj); # фильтрация суммарного сигнала
Посмотрим каким образом преобразился выходной график модели:
plot(r/1000,abs.(y1j),xlabel="Дальность, км",ylabel="Амплитуда",
lab="",xticks=-1e6:5:1e6,title="Выход СФ при воздействии помехи")
Можем заметить, что после добавление помехи на выходе согласованного фильтра сформировалось 2 корреляционных пика, причем 1 пик от помехи преобладает и затеняет истинное положение цели.
3. Перестройка центральной частоты
Для уменьшения влияния помеховой составляющей можно воспользоваться подходом перестраивания частоты, заключающийся в изменение центральной частоты сигнала на значение нескольких полос сигнала:
deltaf = 2.5 * bw # смещение центральной частоты на 250 кГц
xh = deepcopy(wav(deltaf)); # сигнал со смещенной центральной частотой
Построим спектрограмму суммы 2 сигналов - несмещенного и смещенного ЛЧМ сигнала:
calc_spectrogramm(x.+xh,1/fs;down_lim = -40)
На спектрограмме видно, что спектры сигналов смещены относительно друг друга примерно на 250 кГц.
Повторим моделирование сценария работы РЛС для смещенного сигнала:
xth = tx(xh) # передатчик
xtharray = txArray(xth, tgtAng) # передающая ФАР
yph = envout(xtharray,radar_pos,target_pos,radar_vel,target_vel) # канал прямого распространения
yrh = target(yph) # цель
yeh = envin(yrh.+xj,target_pos,radar_pos,target_vel,radar_vel) # канал обратного распространения
yth = rxArray(yeh,tgtAng) # приемная ФАР
yth = rxPreamp(yth) # предусилитель
ybfh = beamformer(yth,doa)[1]; # формирование луча
Для выделения полезной составляющей и удаления помеховой компоненты требуется полосовой фильтр, настроенный на смещенную частоту с полосой ЗС.
Зададим коэффициенты фильтра НЧ Баттерворта 9 порядка и сместим центральную полосу на 250 кГц:
include("$(@__DIR__)/data_filter_coeff.jl") # считывание коэффициентов из jl файла
# смещение центральной частоты на 250 кГц
bf2 = buttercoef.*exp.(1im*2*pi*deltaf*Vector(0:length(buttercoef)-1)/fs);
Реализуем фильтрация в полосе ЗС:
mfcoeff2 = getMatchedFilter(wav,deltaf)
mf2 = EngeePhased.MatchedFilter(Coefficients=mfcoeff2)
# выделение полосым фильтром и согласованная фильтрация
yb2 = mf2(DSP.filt(bf2,1,ybfh))
plot(r./1000,abs.(yb2),lab="")
xlabel!("Дальность, км")
ylabel!("Амплитуда")
title!("Выход СФ после перестройки")
Из графика сигнала на выходе СФ можно заметить, что корреляционный пик помеховой составляющей успешно подавлен.
Таким образом, применение перестройки центральной частоты позволило
Заключение
В примере был рассмотрено моделирование РЛС, перестраиваемого по частоте, при использовании системных объектов библиотеки EngeePhased. Такой подход позволил уменьшить воздействие помехой компоненты и успешно выделить полезный сигнал.