Оценка дальности и доплера с помощью сжатия импульсов
Этот пример демонстрирует эффект сжатия импульсов, когда передаваемый импульс коррелирует с принимаемым сигналом. Радарные и гидролокационные системы используют сжатие импульсов для улучшения соотношения сигнал/шум (SNR) и разрешения по дальности за счет сокращения длительности отраженных сигналов. Этот пример также демонстрирует доплеровскую обработку, при которой радиальная скорость цели определяется по доплеровскому сдвигу, возникающему при движении цели.
Определение дальности и скорости движения целей
Представленная ниже радиолокационная система передает импульсную последовательность линейных частотно-модулированных (ЛЧМ) волн на движущуюся цель и принимает эхо-сигналы. Применяя согласованную фильтрацию и доплеровскую обработку, радарная система может эффективно определять дальность и скорость цели.
Зададим требования к радиолокационной системе. В данном примере используется несущая частота fc = 3 ГГц, и частота дискретизации fs = 1 мГц.
using Pkg; Pkg.add("DSP")
using FFTW,DSP
fc = 3e9; # несущая частота, Гц
fs = 1e6; # частота дискретизации, Гц
c = 299_792_458; # м/с
Создайте объекты системы для моделирования радиолокационной системы. Система является моностатической. Передатчик расположен в точке (0,0,0) м и неподвижен, а цель находится в точке (5000,5000,0) м со скоростью (25,25,0) м/с и имеет эффективную площадь рассеяния (ЭПР) 1 кв. м.
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 кГц. Коэффициенты согласованного фильтра генерируются на основе этой сигнала.
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, причем Матрица данных хранит быстрые временные выборки (время в пределах каждого импульса) в каждом столбце и медленные временные выборки (время между импульсами) в каждой строке.
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,DSP.pow2db.((abs.(resp)).^2),ylims=(0,12e3),
fontfamily=font(14,"Computer Modern"), xticks=(-200:50:200),
colorbar_title = "Мощность,дБВт",
title=("Дальностно-доплероский портрет"),
margin=5Plots.mm
)
xlabel!("Скорость, м/с")
ylabel!("Дальность, м")
Ответ показывает, что цель движется со скоростью около -40 м/с, и что цель удаляется от передатчика, поскольку скорость отрицательна.
Рассчитаем оценку дальности, соответствующие скорости распространения сигнала. График первого образца медленного времени показывает наибольший пик в районе 7000 м, что соответствует графику отклика на скорость дальности.
fasttime = unigrid(0,1/fs,1/prf,"[)"); # сетка быстрого времени
rangebins = (c*fasttime/2); # разрешение по дальности
# Построение зависимости амплитуда от дальности
plot(rangebins.*1e-3,
abs.(rxsig[:,1].*1e6),
lab="",xlabel="Дальность, км",
ylabel="Амплитуда, мкВ"
)
Определение дальности
Зададим порог обнаружения, при котором вероятность ложной тревоги будет меньше . Выполним некогерентное накопление 64 импульсов, предполагая, что сигнал находится в белом гауссовском шуме. По наибольшему пику, превышающего пороговое значение, определим оценку дальности цели.
pfa = 1e-6; # вероятность ложной тревоги
NoiseBandwidth = 5e6/2; # шумовая полоса
# Расчет мощности шума
npower = noisepow(
NoiseBandwidth, # Шумовая полоса
receiver.NoiseFigure, # фактор шума
receiver.ReferenceTemperature # шумовая температура
);
thresh = npwgnthresh(pfa,numPulses,"noncoherent"); # порог обнаружение
thresh = npower*DSP.db2pow(thresh); # чувствительность приемника
pks,range_detect = findpeaks(pulsint(rxsig,"noncoherent"),"MinPeakHeight",thresh); # обнаружение пиков
range_estimate = rangebins[range_detect[1]]; # определение оценки дальности
print("Дальность: $(round(range_estimate;sigdigits=6)), м")
Определение скорости
Выделим выборку по дальности, котора включает обнаруженную цель, и построим график оценки спектральной плотности мощности с помощью функции periodagram.
out = DSP.periodogram(rxsig[range_detect[1],:][:]; onesided=false,
nfft=256, fs=prf, window=nothing)
spec_power = fftshift(DSP.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")
Пиковая частота соответствует доплеровскому сдвигу, деленному на 2, который можно преобразовать в скорость цели. Положительная скорость означает, что цель приближается к излучателю, а отрицательная - что цель удаляется от излучателя.
Y,I = maximum(spec_power),argmax(spec_power);
lambda = c/fc;
tgtspeed = dop2speed(freq[I]/2,lambda);
print("Оценка доплеровского сдвига: $(freq[I]/2) Гц")
print("Оценка скорости: $(round(tgtspeed;sigdigits=4)) м/с")