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

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

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

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

using DSP,FFTW

using SpecialFunctions

function dop2speed(dp,lambda) # пересчет из доплеровской частоты в скорость
    (eltype(dp) <: Union{Real,Array{<:Real}} &&  all(isnan.(dp) .⊻ true ) &&
        all(isempty.(dp) .⊻ true) && all(isinf.(dp) .⊻ true)) ||
        throw(ArgumentError("Expected sp to be real value."))
    (eltype(lambda) <: Real && lambda >= 0 && !isinf(lambda) && (isempty(lambda) ⊻ true) && 
        (isnan(lambda) ⊻ true)) || throw(ArgumentError("Expected lambda to be positive scalar."))
    return dp.*lambda
end

function findpeaks(x,type,tresh)
    if type == "MinPeakHeight"
        value_peak = map(a -> a > tresh,x).*x
        index = (map(x -> x > tresh,x).*(1:length(x)))
        a = index[index.!=0]
        global i = 1
        out_index = Int[]
        while i <= length(a)
            global i
            if (a[i+1 > length(a) ? end : i+1] - a[i]) == 1
                len_idx = any(index[a[i]:end].==0) ? findfirst(index[a[i]:end].==0)-1 : length(a)-i+1
                push!(out_index,argmax(value_peak[a[i]:a[i+len_idx-1]])+a[i]-1)
                i += len_idx
            else
                push!(out_index,a[i])
                i += 1
            end
        end
        return x[out_index], out_index
    end
end

function noisepow(nbw::Number, nf::Number = 0, reftemp::Number = 290)
    function systemp(nf::Number,reftemp::Number = 290)
        nf >= 0 || error("Invalid input")
        nfLinear = db2pow(nf); 
        stemp = reftemp*(nfLinear); # Kelvin
        return stemp
    end
    B = 1.380649e-23
    npow =  B * systemp(nf,reftemp) * nbw
    return npow 
end

function buffer(x,n,p=0,method="delay")
    method in ["delay","nodelay"] || throw(ArgumentError("Method to be equal `delay` or `nodelay`"))
    (ndims(x) in [1,2] && length(x) == max(size(x)...)) || throw(ArgumentError("Expected x is vector "))
    p < n || throw(ArgumentError("Expected value p < n"))
    x_new = method =="delay" ? [zeros(eltype(x),p);x[:]] : x[:]
    x_buf = reshape([x_new;zeros(n-mod(length(x_new),n))],n,:)
    x_buf = vcat(x_buf,x_buf[end-p+1:end,:])
    n_delay = Int(n-(mod(length(x_buf),n)))*(Int(n-(mod(length(x_buf),n)))!=n) 
    return reshape([x_buf[:] ;zeros(eltype(x),n_delay)],n,:)
end

function fspl(R, lambda)
    (eltype(R) <: Union{Real,Array{<:Real}} &&  all(isnan.(R) .⊻ true ) &&
        all(isempty.(R) .⊻ true) && all(isinf.(R) .⊻ true)) ||
        throw(ArgumentError("Expected R to be real finite vector or scalar."))
    
    (eltype(lambda) <: Real && lambda >= 0 && !isinf(lambda) && (isempty(lambda) ⊻ true) && 
        (isnan(lambda) ⊻ true)) || throw(ArgumentError("Expected lambda to be positive scalar."))
    
    fspl_validate(L) = L < 1 ? 1 : 20*log10(L)
    return fspl_validate.(4pi*R./lambda)
end

function BasicMonostaticRadarExampleData()
    c = 299_792_458
    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);
    collector = EngeePhased.Collector(Sensor=antenna,
        PropagationSpeed=c,Wavefront="Plane",
        OperatingFrequency=1e10);
    transmitter = EngeePhased.Transmitter(PeakPower=5226.4791642003665401716716587543,
        Gain=20,LossFactor=0, InUseOutputPort=true,CoherentOnTransmit=true);
    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;

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

Сначала мы дадим определение радарной системы. Поскольку в данном примере основное внимание уделяется доплеровской обработке, мы используем радарную систему, построенную в примере Simulating Test Signals for a Radar Receiver. Читателям предлагается изучить детали проектирования радарных систем на этом примере.

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

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

Цели

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

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);

Обратите внимание, что первая и третья цели находятся на расстоянии 2000 м и движутся со скоростью 100 м/с. Разница в том, что первая цель движется в радиальном направлении, а третья - в тангенциальном. Вторая мишень не движется.

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

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

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;

# Pre-allocate array for improved processing speed
rxpulses = zeros(ComplexF64,length(fast_time_grid),num_pulse_int);

for m = 1:num_pulse_int
    
    # Update sensor and target positions
    sensorpos,sensorvel = sensormotion(1/prf);
    tgtpos,tgtvel = tgtmotion(1/prf);

    # Calculate the target angles as seen by the sensor
    tgtrng,tgtang = rangeangle(tgtpos,sensorpos);
    
    # Simulate propagation of pulse in direction of targets
    pulse = waveform();
    txsig,txstatus = transmitter(pulse);
    txsig = radiator(txsig,tgtang);
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel);
    
    # Reflect pulse off of targets
    tgtsig = target(txsig);
    
    # Receive target returns at sensor
    rxsig = collector(tgtsig,tgtang);
    rxpulses[:,m] = receiver(rxsig,xor.(txstatus.>0,1));
end

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

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

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

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

# calculate initial threshold
pfa = 1e-6;
# in loaded system, noise bandwidth is half of the sample rate
noise_bw = receiver.SampleRate/2; 
npower = noisepow(noise_bw,
    receiver.NoiseFigure,receiver.ReferenceTemperature);
threshold = npower * db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent"));

# apply matched filter and update the threshold
matchingcoeff = getMatchedFilter(waveform);
matchedfilter = EngeePhased.MatchedFilter(
    Coefficients=matchingcoeff,
    GainOutputPort=true);
rxpulses, mfgain = matchedfilter(rxpulses);
threshold = threshold * db2pow(mfgain);

# compensate the matched filter delay
matchingdelay = size(matchingcoeff,1)-1;
rxpulses = buffer(rxpulses[(matchingdelay+1):end],size(rxpulses,1));

# apply time varying gain to compensate the range dependent loss
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);

# detect peaks from the integrated pulse
_,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("max_speed = $(round(max_speed;sigdigits=7))");

max_speed = 224.6888

speed_res = 2*max_speed/num_pulse_int
println("speed_res = $(round(speed_res;sigdigits=7))");

speed_res = 44.93776

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

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

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

num_range_detected = length(range_estimates);
period_out_1 = periodogram(rxpulses[range_detect[1],:][:],nfft=256,fs=prf)
p1, f1 = period_out_1.power,period_out_1.freq
period_out_2 = periodogram(rxpulses[range_detect[2],:][:],nfft=256,fs=prf)
p2, f2 = period_out_2.power,period_out_2.freq;

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

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)
xlabel!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("СПМ,дБВт/Гц", fontfamily="Computer Modern")
title!("Спектральная плотность мощности отраженного сигнала", 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 = $sp1");

sp1 = [-103.56749129975047, 3.5107624169406937]

Результаты показывают, что на дальности 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)
xlabel!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("Мощность,дБВт", fontfamily="Computer Modern")
title!("Спектральная плотность мощности", 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 = $sp2");

sp2 = [0.0]

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

Заключение

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