Доплеровская оценка

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

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

Используемые фукции и библиотеки

using DSP,FFTW

function BasicMonostaticRadarExampleData()
    c = 299_792_458 # cкорость распространения сигнала
    # Создание системных объектов:
    # генератор прямоугольных импульсов
    waveform = EngeePhased.RectangularWaveform(PulseWidth=100/c,  
        PRF=c/1e4,OutputFormat="Pulses",NumPulses=1,SampleRate= 2c/100);
    # антенный элемент (тип - изотропный)
    antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[5e9 1.5e10]);
    # излучатель узкополосного сигнала
    radiator = EngeePhased.Radiator(Sensor=antenna,
        PropagationSpeed=c,OperatingFrequency=1e10);
    # Усилитель в передающем тракте 
    transmitter = EngeePhased.Transmitter(PeakPower=5226.4791642003665401716716587543,
        Gain=20,LossFactor=0, InUseOutputPort=true,CoherentOnTransmit=true);
    # принимающая антенна 
    collector = EngeePhased.Collector(Sensor=antenna,
        PropagationSpeed=c,Wavefront="Plane",
        OperatingFrequency=1e10);
    # Предусилитель в приемной тракте 
    receiver = EngeePhased.ReceiverPreamp(Gain=20,NoiseFigure=0, 
        ReferenceTemperature=290,SampleRate=2*c/100, PhaseNoiseInputPort=false,
        EnableInputPort=true,SeedSource="Property",Seed=1000);
    # модель движения радара 
    sensormotion = EngeePhased.Platform(InitialPosition=[0;0;0],Velocity=[0;0;0]);

    return waveform, antenna, radiator, collector, transmitter, receiver, sensormotion
end;

Настройка радиолокационной системы

Сначала мы дадим определение радарной системы. Поскольку в данном примере основное внимание уделяется доплеровской обработке, мы используем радарную систему, построенную в примере Моделирование тестовых сигналов для радиолокационного приемника. Читателям предлагается изучить детали проектирования радарных систем на этом примере.

waveform,antenna,radiator,collector,transmitter,receiver,sensormotion = BasicMonostaticRadarExampleData();

Моделирование системы

Цели

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

tgtpos = [[1200; 1600; 0] [3543.63; 0; 0] [1600; 0; 1200]];  # компоненты начального положения целей
tgtvel = [[60; 80; 0] [0;0;0] [0; 100; 0]]; # компоненты скорости целей
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # модель движения целей

tgtrcs = [1.3 1.7 2.1]; # ЭПР целей
fc = radiator.OperatingFrequency;
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc);

Обратите внимание, что первая цель находятся на расстоянии 3550, 1600 и 1200 м , и движутся со скоростью 60, 130 и 0 м/с, . Разница в том, что первая цель движется в радиальном направлении, а вторая - в тангенциальном и радиальном. Третья цель не движется.

Среда распространения

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

fs = waveform.SampleRate; # частота дискретизации системы
# модель канала распространения
channel = EngeePhased.FreeSpace(
    SampleRate=fs, # частота дискретизации системы
    TwoWayPropagation=true, # модель распространения
    OperatingFrequency=fc); # несущая частота сигнала

Синтез сигналов

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

Чтобы воспроизвести те же результаты, нам нужно задать исходные данные для генерации шума в приемнике.

release!(receiver) # сброс внутреннего состояния объекта (измненеие )
receiver.SeedSource = "Property"; # изменение тип модели генерации шума
receiver.Seed = 2024; # задание начального состояния генератора тепловых шумов

prf = waveform.PRF; # частота следования импульсов
num_pulse_int = 10; # количество импульсов накопления

fast_time_grid =  (0:1/fs:1/prf-1/fs) # сетка быстрого времени
slow_time_grid = (0:num_pulse_int-1)./prf; # сетка медленного времени

# Выделение памяти для результирующих данных
rxpulses = zeros(ComplexF64,length(fast_time_grid),num_pulse_int);

for m = 1:num_pulse_int
    
    # Обновление пложения радара и целей
    sensorpos,sensorvel = sensormotion(1/prf);
    tgtpos,tgtvel = tgtmotion(1/prf);

    # Расчет углов по дальности от радара до целей  
    tgtrng,tgtang = rangeangle(tgtpos,sensorpos);
    
    # Моделирование распространения сигнала от радара до целей
    pulse = waveform(); # генератор прямоугольных импульсво
    txsig,txstatus = transmitter(pulse); # усилитель передающего сигнала
    txsig = radiator(txsig,tgtang); # излучатель сгенерированного сигнала в среду
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel); # среда распространения
    
    # Отражение сигнала от целей с учетом ЭПР
    tgtsig = target(txsig);
    
    # Прием и предварительного усиление отраженного сигнала
    rxsig = collector(tgtsig,tgtang); # приемная антенна
    rxpulses[:,m] = receiver(rxsig,xor.(txstatus.>0,1)); # предусилитель
end

Доплеровская оценка

Обнаружение дальности

Чтобы оценить доплеровский сдвиг целей, сначала необходимо определить их местоположение с помощью обнаружения по дальности. Поскольку доплеровский сдвиг распределяет мощность сигнала как по I, так и по Q каналам, для обнаружения необходимо использовать энергию сигнала. Это означает, что необходимо применять некогерентные схемы обнаружения.

Процесс обнаружения подробно описан в вышеупомянутом примере, поэтому здесь мы просто выполним необходимые действия для оценки дальности цели.

# Задание вероятности ложной тревоги (pfa)
pfa = 1e-6;
# В рассматриваемой системе, шумовая полоса равна половине частоты дискретизации
noise_bw = receiver.SampleRate/2; # шумовая полоса
npower = noisepow(noise_bw,
    receiver.NoiseFigure,receiver.ReferenceTemperature); # мощность шума
threshold = npower * db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # порог обнаружения

# Применение согласованной фильтрации и обновление порога обнаружения
matchingcoeff = getMatchedFilter(waveform); # расчет коэффициентов согласованного 
                                            # фильтра на основе изулченного сигнала
matchedfilter = EngeePhased.MatchedFilter( # создание объекта согласованного фильтра
    Coefficients=matchingcoeff, # коэффициенты фильтра
    GainOutputPort=true); # выходной порт усиления
rxpulses, mfgain = matchedfilter(rxpulses); # отклик после согласованной фильтрации
threshold = threshold * db2pow(mfgain); # расчет нового значения порога обнаружения

# Компенсация задержки в согласованном фильтре
matchingdelay = size(matchingcoeff,1)-1; # длина согласованного фильтра
rxpulses = buffer(rxpulses[(matchingdelay+1):end],size(rxpulses,1));

# Применение временного автоматического регулирования усиления (ВАРУ)
#   для компенсации потерь, зависящие от дальности.
prop_speed = radiator.PropagationSpeed; # скорость распространения сигнала
range_gates = prop_speed*fast_time_grid/2; # разрешение по дальности
lambda = prop_speed/fc; # длина волны

tvg = EngeePhased.TimeVaryingGain( # создание системного объекта ВАРУ
    RangeLoss=2*fspl.(range_gates,lambda), # потери с учетом дальности
    ReferenceLoss=2*fspl.(prop_speed/(prf*2),lambda)); # потери с учетом дальности

rxpulses = tvg(rxpulses); # применение ВАРУ

# Обнаружения пиков, привышающих порог обнаружения, после интегрирования
_,range_detect = findpeaks(pulsint(rxpulses,"noncoherent"),
    "MinPeakHeight",sqrt(threshold)); # нахождение пиков
range_estimates = round.(range_gates[range_detect]) # оценки дальности

2-element Vector{Float64}:
 2000.0
 3550.0

Эти оценки позволяют предположить наличие целей в диапазоне 2000 м и 3550 м.

Доплеровский спектр

После того как успешно оценены дальности до целей, можно оценить доплеровскую информацию для каждой цели.

Оценка доплеровского сигнала - это, фактически, процесс оценки спектра. Поэтому, первым шагом в доплеровской обработке является расчет доплеровского спектра принимаемого сигнала.

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

Как предсказывает теория Фурье, максимальный однозначный доплеровский сдвиг, который может обнаружить импульсная локационная система, равен половине ее PRF. Это также соответствует максимальной однозначной скорости, которую может обнаружить радилокационная система. Кроме того, количество импульсов определяет разрешение в доплеровском спектре, от которого зависит разрешение оценки скорости.

max_speed = dop2speed(prf/2,lambda)/2; # расчет 
println("Максимальная скорость обнаружения: $(round(max_speed;sigdigits=7)) м/с");

Максимальная скорость обнаружения: 224.6888 м/с

speed_res = 2*max_speed/num_pulse_int
println("Разрешение по скорости: $(round(speed_res;sigdigits=7)) м/с");

Разрешение по скорости: 44.93776 м/с

Как показано в расчетах выше, в данном примере максимальная обнаруживаемая скорость составляет 225 м/с, либо приближающаяся (-225), либо удаляющаяся (+225). Результирующее доплеровское разрешение составляет около 45 м/с, что означает, что для разделения в доплеровском спектре две скорости должны отличаться друг от друга не менее чем на 45 м/с. Чтобы улучшить способность различать разные скорости цели, необходимо большее количество импульсов. Однако количество доступных импульсов также ограничено радиальной скоростью цели. Поскольку доплеровская обработка ограничена заданным диапазоном, все импульсы, используемые в обработке, должны быть приняты до того, как цель перейдет из одного элемента разрешения в другой.

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

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

num_range_detected = length(range_estimates);
# Расчет периограммы для дальности 2000 м
period_out_1 = periodogram(rxpulses[range_detect[1],:][:],nfft=256,fs=prf) 
p1, f1 = period_out_1.power,period_out_1.freq # излечение значений мощности (p1) и частот (f1)
# Расчет периограммы для дальности 3550 м
period_out_2 = periodogram(rxpulses[range_detect[2],:][:],nfft=256,fs=prf)
p2, f2 = period_out_2.power,period_out_2.freq; # излечение значений мощности (p2) и частот (f2)

Затем можно рассчитать скорость, соответствующую каждому значению частоты в спектре. Обратите внимание на то, что необходимо учесть эффект кругового движения.

speed_vec = dop2speed.(fftshift(f1),lambda)/2; # пересчет сетки частот в сетку скорости

Оценка доплеровского сдвига

Чтобы оценить доплеровский сдвиг, связанный с каждой целью, нам нужно найти местоположение пиков в каждом доплеровском спектре. В данном примере цели присутствуют на двух разных диапазонах, поэтому процесс оценки необходимо повторить для каждого диапазона.

Сначала построим доплеровский спектр, соответствующий дальности 2000 метров.

spec_sig1 = fft([rxpulses[range_detect[1],:][:];zeros(256-length(rxpulses[range_detect[1],:][:]))]) 
spec_power1 = fftshift(real(spec_sig1.*conj.(spec_sig1))./100)
plot(fftshift(fftfreq(256,prf)).*1e-3,pow2db.(spec_power1),lab="",gridalpha=0.15,margin=5*Plots.mm)
xlabel!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("СПМ,дБВт/Гц", fontfamily="Computer Modern")
title!("Спектральная плотность мощности отраженного сигнала по дальности 2000 м", fontfamily="Computer Modern")

Обратите внимание, что нас интересует только обнаружение пиков, поэтому сами значения спектра не имеют значения. Из графика доплеровского спектра видно, что 5 дБ ниже максимального пика - это хороший порог. Поэтому мы используем значение -5 в качестве порога для нормализованного доплеровского спектра.

spectrum_data = spec_power1./maximum(spec_power1)
_,dop_detect1 = findpeaks(pow2db.(spectrum_data),"MinPeakHeight",-5)
sp1 = speed_vec[dop_detect1]
println("Оценки скорости:  $sp1 м/с");

Оценки скорости:  [-103.56749129975046, 3.510762416940684] м/с

Результаты показывают, что на дальности 2000 м есть две цели: одна со скоростью 3,5 м/с, а другая -104 м/с. Значение -104 м/с можно легко связать с первой целью, поскольку первая цель уходит с радиальной скоростью 100 м/с, что, учитывая доплеровское разрешение данного примера, очень близко к расчетному значению. Значение 3,5 м/с требует более подробного объяснения. Поскольку третья цель движется по тангенциальному направлению, компонент скорости в радиальном направлении отсутствует. Поэтому радар не может обнаружить доплеровское смещение третьей цели. Истинная радиальная скорость третьей цели, следовательно, равна 0 м/с, а оценка в 3,5 м/с очень близка к истинному значению.

Обратите внимание, что эти две цели невозможно различить, используя только оценку дальности, поскольку их значения дальности одинаковы.

Затем те же операции применяются к данным, соответствующим дальности 3550 метров.

spec_sig2 = fft([rxpulses[range_detect[2],:][:];zeros(256-length(rxpulses[range_detect[2],:][:]))]) 
spec_power2 = fftshift(real(spec_sig2.*conj.(spec_sig2))./100)
plot(fftshift(fftfreq(256,prf)).*1e-3,pow2db.(spec_power2),lab="",gridalpha=0.15,margin=5*Plots.mm)
xlabel!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("Мощность,дБВт", fontfamily="Computer Modern")
title!("Спектральная плотность мощности по дальности 3550 метров", fontfamily="Computer Modern")

spectrum_data2 = spec_power2./maximum(spec_power2)
_,dop_detect2 = findpeaks(pow2db.(spectrum_data2),"MinPeakHeight",-5)
sp2 = speed_vec[dop_detect2]
println("Оценка скорости: $sp2 м/с");

Оценка скорости: [0.0] м/с

Этот результат показывает расчетную скорость 0 м/с, что соответствует тому факту, что цель на этой дистанции не движется.

Заключение

В этом примере показан простой способ оценки радиальной скорости движущихся целей с помощью импульсной радиолокационной системы. Сгенерировали доплеровский спектр принимаемого сигнала и оценили расположение пиков по спектру, которые соответствуют радиальной скорости цели. Также обсудиили ограничения доплеровской обработки вызванные характеристиками излученного сигнала и алгоритмом обработки.