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;

导言

电子扫描有源相控阵天线可以帮助解决使用相同天线阵列设备的各种任务。 它们可用于雷达，电子战和通信系统。 然而，这些类型的系统在其中操作的环境是复杂的并且可能引入不需要的干扰。 例如，中继器上的干扰器可能会重复接收到的雷达信号并重新传输它以混淆雷达。 频率适应可以是对抗由干扰源产生的信号并有助于这些系统的高效操作的有效方法。

在此示例中，对具有静止单稳态雷达和移动飞机目标的场景进行建模（在示例Center frequency tuning under impact中考虑了类似的场景помех)。然后，飞机将产生混淆雷达的替换信号。 一旦雷达检测到干扰源，就可以使用中心频率调谐技术来允许雷达绕过干扰。

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的方向性，我们将形成一个2维泰勒窗函数。:

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

我们可以看到，添加干扰后，在匹配滤波器的输出端形成2个相关峰，其中1个峰来自干扰，模糊了目标的真实位置。

3. 调谐中心频率

为了减少干扰分量的影响，可以使用频率调谐方法，该方法包括将信号的中心频率改变为几个信号频带的值。:

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

让我们构建一个2个信号之和的频谱图-一个无偏和一个偏移FM信号。:

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

频谱图显示信号频谱相对于彼此偏移约250kHz。

让我们重复一个偏移信号的雷达操作场景的模拟:

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

为了隔离有用分量并去除干扰分量，需要一个带通滤波器，调谐到具有ZS频带的偏移频率。

让我们设置9阶低通巴特沃斯滤波器的系数，并将中心带移位250kHz:

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

我们在SW波段实现滤波:

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

从SF输出端的信号的曲线图可以看出，干扰分量的相关峰值已经被成功地抑制。
因此，允许使用中心频率调谐

结论

在示例中，我们考虑了使用EngeePhased库的系统对象对频率可调谐雷达的仿真。 这种方法使得能够减小干扰分量的影响并成功地隔离有用信号。