Engee 文档
Notebook

干扰下的频率适应

本例说明了如何对雷达系统进行干扰下的频率适应建模。

连接库和辅助功能

In [ ]: 
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
In [ ]: 
using DSP
In [ ]: 
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;

导言

有源相控电子扫描天线阵列可以使用相同的天线阵列硬件帮助解决许多应用问题。它们可用于雷达、REB 和通信系统。然而，这些系统的运行环境非常复杂，可能会引入不必要的干扰。例如，中继器上的干扰装置可能会重复接收雷达信号并重新发射，从而混淆雷达。频率适应是抵消干扰源产生的信号的有效方法，有助于这些系统的有效运行。

本例模拟了静止的单静态雷达和移动的飞机目标的情景（类似情景在在干扰情况下调谐中心频率一例中讨论）。飞机随后会产生一个欺骗信号，使雷达产生混淆。一旦雷达探测到干扰源，就可以使用中心跳频技术使雷达绕过干扰。

1.无干扰环境下的系统建模

假设原点有一个单静态 X 波段雷达。

In [ ]: 
fc = 10e9; # несущая частота работы РЛС X-диапазона, Гц 
fs = 2e6; # частота дискретизации модели, Гц
c = 3e8; # скорость распросранения сигнала, м/с
lambda = c/fc; # длина волны, м

radar_pos = [0;0;0]; # положение РЛС
radar_vel = [0;0;0]; # скорость РЛС

接收天线由 64 个等距矩形天线阵元件（8x8）组成，元件间距为半波长：

In [ ]: 
antenna = EngeePhased.URA(
    Element=EngeePhased.CosineAntennaElement(), # косинусный элемент АР 
    Size=[8 8], # размер антенной решетки
    ElementSpacing=[lambda/2 lambda/2],  # расстояние между элементами
)
Out[0]:
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

为了更好地描述 AR 的指向性，我们形成了一个二维泰勒窗函数：

In [ ]: 
taper = TaylorSpectrumWindow(8); # взвешанная фукция окно Тэйлора
taperURA = taper.*taper' # формирование 2д взвешанной функции
antenna.Taper = taperURA # взвешивание АР

# Визуализация оконной функции
surface(taperURA,title="2-мерная оконная фукция Тейлора")
Out[0]:

天线阵列指向性图可以使用函数 pattern 绘制：

In [ ]: 
pattern(antenna,fc,Type="powerdb")
plot!(title="ДН прямоугольной АР",colorbar_title="Мощность, дБВт")
Out[0]:

在所考虑的雷达系统中，信号源将是线性频率调制（LFM）信号：

In [ ]: 
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,  # частота дискретизации модели
);

接下来，让我们创建目标和分发环境对象：

In [ ]: 
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); # дублирование канала обратного распространения

利用之前定义的系统对象，我们可以计算出目标在无干扰情况下的反射信号：

In [ ]: 
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

In [ ]: 
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

In [ ]: 
mfcoeff1 = getMatchedFilter(wav) # Коэффициентов согласованнго фильтра
mf1 = EngeePhased.MatchedFilter(Coefficients=mfcoeff1);

对生成的光束进行匹配滤波，并显示雷达操作的结果：

In [ ]: 
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="Выход СФ")
Out[0]:

2.添加干扰成分

现在让我们通过模拟有用信号和时延来模拟干扰信号：

In [ ]: 
 # имитация помехового сигнала путемдублирование генератор сигнала
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); # фильтрация суммарного сигнала

让我们看看模型的输出图是如何变化的：

In [ ]: 
plot(r/1000,abs.(y1j),xlabel="Дальность, км",ylabel="Амплитуда",
    lab="",xticks=-1e6:5:1e6,title="Выход СФ при воздействии помехи")
Out[0]:

我们可以注意到，在加入干扰后，匹配滤波器的输出端形成了两个相关峰，其中一个峰主要来自干扰，并遮挡了目标的真实位置。

3.调整中心频率

为了减少干扰成分的影响，可以采用频率调谐方法，即通过几个信号波段的值来改变信号的中心频率：

In [ ]: 
deltaf = 2.5 * bw # смещение центральной частоты на 250 кГц
xh = deepcopy(wav(deltaf)); # сигнал со смещенной центральной частотой

让我们来构建两个信号（无偏置和有偏置的 LFM 信号）之和的频谱图：

In [ ]: 
calc_spectrogramm(x.+xh,1/fs;down_lim = -40)
Out[0]:

该频谱图显示，两个信号的频谱相对偏移了约 250 kHz。

让我们针对移位信号重复模拟雷达运行情况：

In [ ]: 
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]; # формирование луча

需要一个调谐到偏移频率、带宽为 WL 的带通滤波器来提取有用成分并去除干扰成分。

将低通巴特沃斯滤波器的系数设为 9 阶，并将中心带宽移至 250 kHz：

In [ ]: 
include("$(@__DIR__)/data_filter_coeff.jl") # считывание коэффициентов из jl файла
# смещение центральной частоты на 250 кГц
bf2 = buttercoef.*exp.(1im*2*pi*deltaf*Vector(0:length(buttercoef)-1)/fs);

在 WF 频段实现滤波：

In [ ]: 
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!("Выход СФ после перестройки")
Out[0]:

从 NF 输出端的信号图可以看出，噪声成分的相关峰值已被成功抑制。 因此，中心频率调整的应用允许

结论

在案例研究中，使用 EngeePhased 库系统对象对频率可调雷达的建模进行了研究。这种方法减少了干扰成分的影响，并成功隔离了有用信号。