Частотная адаптация при воздействии помехи

В этом примере показано, как моделировать радилокационные системы с частотной адаптацией при воздействии помех.

Подключение библиотек и вспомогательных функций

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],  # расстояние между элементами
)

URA:
    Element=CosineAntennaElement: FrequencyRange=[0.0 1.0e20] CosinePower=[1.5 1.5]
    Size=[8 8]
    ElementSpacing=[0.015 0.015]
    Lattice=Rectangular
    ArrayNormal=x
    Taper=1

Для лучшей характеристики направленности АР, сформируем 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. Такой подход позволил уменьшить воздействие помехой компоненты и успешно выделить полезный сигнал.