Документация Engee
Notebook

Оценка дальности и доплера с помощью сжатия импульсов

Открыть пример в Engee

Этот пример демонстрирует эффект сжатия импульсов, когда передаваемый импульс коррелирует с принимаемым сигналом. Радарные и гидролокационные системы используют сжатие импульсов для улучшения соотношения сигнал/шум (SNR) и разрешения по дальности за счет сокращения длительности отраженных сигналов. Этот пример также демонстрирует доплеровскую обработку, при которой радиальная скорость цели определяется по доплеровскому сдвигу, возникающему при движении цели.

Определение дальности и скорости движения целей

Представленная ниже радиолокационная система передает импульсную последовательность линейных частотно-модулированных (ЛЧМ) волн на движущуюся цель и принимает эхо-сигналы. Применяя согласованную фильтрацию и доплеровскую обработку, радарная система может эффективно определять дальность и скорость цели.

Зададим требования к радиолокационной системе. В данном примере используется несущая частота fc = 3 ГГц, и частота дискретизации fs = 1 мГц.

In [ ]: 
using FFTW
fc = 3e9; # несущая частота, Гц
fs = 1e6; # частота дискретизации, Гц
c = 299_792_458; # м/с

Создайте объекты системы для моделирования радиолокационной системы. Система является моностатической. Передатчик расположен в точке (0,0,0) м и неподвижен, а цель находится в точке (5000,5000,0) м со скоростью (25,25,0) м/с и имеет эффективную площадь рассеяния (ЭПР) 1 кв. м.

In [ ]: 
antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[1e8 10e9]); # антенный элемент
transmitter = EngeePhased.Transmitter(Gain=20,InUseOutputPort=true); # передатчик
txloc = [0;0;0]; # положение передатчика
tgtloc = [5000;5000;0]; # координаты цели, дальность ~= 7071 м
tgtvel = [25;25;0]; # скорость цели,  относительная скорость ~= 35.4 м/с
target = EngeePhased.RadarTarget(Model="Nonfluctuating",MeanRCS=1,OperatingFrequency=fc);
antennaplatform = EngeePhased.Platform(InitialPosition=txloc); # объект передатчика в пространстве
targetplatform = EngeePhased.Platform(InitialPosition=tgtloc,Velocity=tgtvel); # объект цели в пространстве
radiator = EngeePhased.Radiator( # излучатель
   PropagationSpeed=c, # скорость распространения
   OperatingFrequency=fc, # несущая частота cигнала
   Sensor=antenna # антенный элемент
);
channel = EngeePhased.FreeSpace( # канал распространения
   PropagationSpeed=c, # скорость распространения
   OperatingFrequency=fc, # несущая частота
   TwoWayPropagation=false # двухнаправленное направление
);
collector = EngeePhased.Collector(
   PropagationSpeed=c, # скорость распространения
   OperatingFrequency=fc, # несущая частота cигнала
   Sensor=antenna # антенный элемент
);
receiver = EngeePhased.ReceiverPreamp( # передатчик
   NoiseFigure=0, # Фактор шума
   EnableInputPort=true, # входной порт синхронизации с передатчиком
   SeedSource="Property", # источник начального значения генератора ПСП для фазового шума
   Seed=2024 # значение генератора ПСП для фазового шума
);

Создайте ЛЧМ-сигнал с длительностью импульса 10 мкс, частотой повторения импульсов 10 кГц и шириной полосы развертки 100 кГц. Коэффициенты согласованного фильтра генерируются на основе этой сигнала.

In [ ]: 
waveform = EngeePhased.LinearFMWaveform( # генератор ЛЧМ сигнала
    PulseWidth=10e-6, # длительность импульса
    PRF=10e3, # частота следования импульсов
    OutputFormat="Pulses", # тип выхода
    NumPulses=1, # количество импульсов
    SweepBandwidth=1e5 # время наростания пилы
);
wav = waveform(); # формирование ЛЧМ-сигнала 
maxrange = c/(2*waveform.PRF); # Максимальная одназначная дальность обнаружения цели
SNR = npwgnthresh(1e-6,1,"noncoherent"); # расчет ОСШ на основе порога обнаружения
lambda = c/fc; # длина волны
tau = waveform.PulseWidth; # длительность импульса
Ts = 290; # Шумовая температура приемника

# Согласованный фильтр
Mfilter = EngeePhased.MatchedFilter( 
   Coefficients=getMatchedFilter(waveform), # коэффициенты фильтра
   GainOutputPort=true # выходной порт коэффициента передачи
);

Для улучшения доплеровского разрешения система излучает 64 импульса. Данные сигнала будем хранить в переменной rxsig, причем Матрица данных хранит быстрые временные выборки (время в пределах каждого импульса) в каждом столбце и медленные временные выборки (время между импульсами) в каждой строке.

In [ ]: 
numPulses = 64; # количество импульсов
rxsig = zeros(ComplexF64,length(wav),numPulses); # выделение памяти под выходной сигнал

for n = 1:numPulses
    tgtloc,tgtvel = targetplatform(1/waveform.PRF); # параметры движения цели
    tgtrng,tgtang = rangeangle(tgtloc,txloc); # параметры движения радара
    
    txsig, txstatus = transmitter(wav); # формирование передающего сигнала 
    txsig = radiator(txsig,tgtang); # излучение ЭМ-волны в пространство
    txsig = channel(txsig,txloc,tgtloc,[0;0;0],tgtvel); # канал: передатчик-цель
    txsig = target(txsig); # Отражение сигнала с учетом ЭПР цели
    txsig = channel(txsig,tgtloc,txloc,tgtvel,[0;0;0]); # канал: цель-приемник
    txsig = collector(txsig,tgtang); # формирование сигнала в приемной антенне
    rxsig[:,n] = receiver(txsig,xor.(txstatus,1)); # предусиление принятого сигнала
end

prf = waveform.PRF; # частота следования импульсов
fs = waveform.SampleRate; # частота дискретизации
response = EngeePhased.RangeDopplerResponse(
    DopplerFFTLengthSource="Property", # задание источника длины БПФ
    DopplerFFTLength=2048, # Длина БПФ по доплеру
    SampleRate=fs, # частота дискретизации
    DopplerOutput="Speed", # тип выхода 
    OperatingFrequency=fc, # несущая частота
    PRFSource="Property", # источник задания ЧСИ
    PRF=prf # частота следования импульсов
);
filt = getMatchedFilter(waveform); # коэффициенты согласованного фильтра
resp,rng_grid,dop_grid = response(rxsig,filt); # отклик по дальности и скорости
rng_grid = range(extrema(rng_grid)...,length = length(rng_grid)) # cетка дальности
dop_grid = range(extrema(dop_grid)...,length = length(dop_grid)) # cетка скорости

# построение долероского-дальностного портрета
heatmap(dop_grid,rng_grid,pow2db.((abs.(resp)).^2),ylims=(0,12e3), 
    fontfamily=font(14,"Computer Modern"), xticks=(-200:50:200),
    colorbar_title = "Мощность,дБВт",
    title=("Дальностно-доплероский портрет"),
    margin=5Plots.mm
)
xlabel!("Скорость, м/с")
ylabel!("Дальность, м")
Out[0]:

Ответ показывает, что цель движется со скоростью около -40 м/с, и что цель удаляется от передатчика, поскольку скорость отрицательна.

Рассчитаем оценку дальности, соответствующие скорости распространения сигнала. График первого образца медленного времени показывает наибольший пик в районе 7000 м, что соответствует графику отклика на скорость дальности.

In [ ]: 
fasttime = unigrid(0,1/fs,1/prf,"[)"); # сетка быстрого времени
rangebins = (c*fasttime/2); # разрешение по дальности

# Построение зависимости амплитуда от дальности
plot(rangebins.*1e-3,
    abs.(rxsig[:,1].*1e6),
    lab="",xlabel="Дальность, км",
    ylabel="Амплитуда, мкВ"
)
Out[0]:

Определение дальности

Зададим порог обнаружения, при котором вероятность ложной тревоги будет меньше $10^{-6}$. Выполним некогерентное накопление 64 импульсов, предполагая, что сигнал находится в белом гауссовском шуме. По наибольшему пику, превышающего пороговое значение, определим оценку дальности цели.

In [ ]: 
pfa = 1e-6; # вероятность ложной тревоги
NoiseBandwidth = 5e6/2; # шумовая полоса

# Расчет мощности шума
npower = noisepow(
    NoiseBandwidth, # Шумовая полоса 
    receiver.NoiseFigure, # фактор шума
    receiver.ReferenceTemperature # шумовая температура
);

thresh = npwgnthresh(pfa,numPulses,"noncoherent"); # порог обнаружение
thresh = npower*db2pow(thresh); # чувствительность приемника

pks,range_detect = findpeaks(pulsint(rxsig,"noncoherent"),"MinPeakHeight",thresh); # обнаружение пиков
range_estimate = rangebins[range_detect[1]]; # определение оценки дальности
print("Дальность: $(round(range_estimate;sigdigits=6)), м")

Дальность: 7344.92, м

Определение скорости

Выделим выборку по дальности, котора включает обнаруженную цель, и построим график оценки спектральной плотности мощности с помощью функции periodagram.

In [ ]: 
out = periodogram(rxsig[range_detect[1],:][:]; onesided=false, 
    nfft=256, fs=prf, window=nothing)
spec_power = fftshift(pow2db.(out.power))
freq = fftshift(fftfreq(256,prf))
plot(freq.*1e-3,spec_power,lab="",gridalpha=0.15)
xlabel!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("Мощность,дБВт", fontfamily="Computer Modern")
title!("Спектральная плотность мощности", fontfamily="Computer Modern")
Out[0]:

Пиковая частота соответствует доплеровскому сдвигу, деленному на 2, который можно преобразовать в скорость цели. Положительная скорость означает, что цель приближается к излучателю, а отрицательная - что цель удаляется от излучателя.

In [ ]: 
Y,I = maximum(spec_power),argmax(spec_power);
lambda = c/fc;
tgtspeed = dop2speed(freq[I]/2,lambda);
print("Оценка доплеровского сдвига: $(freq[I]/2) Гц")

Оценка доплеровского сдвига: -351.5625 Гц
In [ ]: 
print("Оценка скорости: $(round(tgtspeed;sigdigits=4)) м/с")

Оценка скорости: -35.13 м/с