Локализация целей на основе активных РЛС

В этом примере показано, как моделировать сценарий работы активных радиолокационных систем (РЛС) и определять оценки времени задержки распространения сигнала (TOA) для решения задачи локализацию при наличии малоразмерной цели.

Этот пример является первой частью из серии решения задачи локализации. Вторая часть посвящена методу пассивных радиолокации с применением алгоритма TDOA (Time Difference Of Arrival) для локализации самолета.

Введение

В области радиолокации задача локализации цели на основе данных, полученных от группы пространственно разнесённых датчиков с известными координатами, имеет значительную актуальность. Для решения этой задачи традиционно применяются методы, основанные на измерениях времени прихода (TOA — Time Of Arrival) и разности времён прихода (TDOA — Time Difference Of Arrival) сигналов.

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

В такой радиолокационной системе TOA определяется на основе задержек распространения сигнала между целью и радиолокационным приемопередатчиком.

Подключение вспомогательных функций

Выполним инициализацию функций, необходимых для расчета оценок позиционирования и визуализации методов локализации TOA

isdefined(Main,:init_func) || include("init.jl")

"Вспомогательные функции успешно подключены"

Рассмотрим локализацию (TOA) цели на основе системы моностатических активных РЛС. Каждый из РЛС имеют синхронизированные приемник и передатчик.

1. Формирование сценария работы РЛС

Создадим сценарий работы 5 радаров X-диапазона, имеющих известные координаты, для локализации беспилотного летательного аппарата (БПЛА) с малой эффективной площадью рассеяния (ЭПР).

# Параметры РЛС
fc = 9e9 # [Гц], несущая частота
c = physconst("LightSpeed") # [м/с], скорость распространения сигнала
bw = 200e6 # Гц, полоса сигнала
fs = bw # Гц, частота дискретизации

# Параметры приемо-передающего тракта
Pt = 1 # Вт, пиковая мощнсоть
Gtx = 40 # дБ, усиление передающей антенны
Grx = 40 # дБ, усиление приемной антенны
NF = 2.9 # дБ, коэффициент шума приемного тракта

# Формирование системных объектов (СО) функциональных узлов РЛС
antenna = EngeePhased.IsotropicAntennaElement( # СО изотропного антенного элемента
    BackBaffled=false # учет обратного рассеяния ДН антенны
) 
transmitter = EngeePhased.Transmitter( # СО - передатчик
    Gain=Gtx, # усиление передатчика
    PeakPower=Pt # Пиковая мощность 
)
radiator = EngeePhased.Radiator( # СО - передающая антенна
    Sensor=antenna, # ДН антенного элемента
    OperatingFrequency=fc # несущая частота антенны
)
collector = EngeePhased.Collector( # СО - приемная антенна
    Sensor=antenna, # ДН антенного элемента
    OperatingFrequency=fc # несущая частота антенны
)
receiver = EngeePhased.ReceiverPreamp(
    Gain=Grx, # усиление приемника
    NoiseFigure=NF, # коэффициент шума в тракте
    SampleRate=fs # частота дискретизации
)

# Создание СО - цель
tgtrcs = 1e-2 # м^2, ЭПР цели
target = EngeePhased.RadarTarget(
    MeanRCS=tgtrcs, # среднее значение ЭПР
    PropagationSpeed=c, # скорость распространения сигнала
    OperatingFrequency=fc # несущая частота сигнала
)

# Модель сценария движения цели
tgtpos = [30; 10; 15]; # [pos_x,pos_y,pos_z],м начальный вектор положения цели 
tgtvel = [5; 10; 0]; # [v_x,v_y,v_z],м/с вектор скорости цели 
tgtplatform = EngeePhased.Platform( # СО - модель движения цели
    InitialPosition=tgtpos,  # начальный вектор положения цели 
    Velocity=tgtvel # вектор скорости цели 
) 

# Модель сценария движения РЛС
radarpos = [   # м, начальные положения РЛС    
    0 -30. 100 80 -40; 
    0 50 -40 30 -20; 
    0 -5 7 5 2
]                        
numRadar = size(radarpos,2) # Число РЛС 
radarvel = zeros(3,numRadar) # [v_x,v_y,v_z] [м/с], скорости РЛС
radarplatform = EngeePhased.Platform( # СО - модель движения РЛС
    InitialPosition=radarpos,  # начальный вектор положения  
    Velocity=radarvel # вектор скорости РЛС 
);

2. Формирование зондирующего сигнала (ЗС)

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

# Формирование непрерывного сигнала с ЛЧМ
N = 1024 # число отсчетов за время нарастания (отсчеты быстрого времени )
M = 8 # количество импульсов накопления (отсчеты медленного времени)
freqSpacing = bw/N # Гц, разрешение по частоте
tWave = N/fs; # с, время нарастания пилы (длительность импульса)
fmcwWaveform = EngeePhased.FMCWWaveform(  # системный объект FMCWWaveform
    SweepTime=tWave, # длительность импульса
    SweepBandwidth=bw, # полоса сигнала (ширина спектра)
    SampleRate=fs, # частота дискретизации
    SweepDirection="Up" # направление изменения частоты
)

# Моделирование непрерывного ЛЧМ сигнала для М-каналов
sig_lfm = ComplexF64.(fmcwWaveform()*ones(1,M))
println("Размерность ЗС: $(size(sig_lfm))")

Размерность ЗС: (1024, 8)

Применим функцию plot_sig_and_spec для построение графика осциллограммы и спектрограммы непрерывного сигнала с ЛЧМ:

plot_sig_and_spec(sig_lfm[:,1];fs=fs,name_sig = "непрерывного ЛЧМ")

На спектрограмме можно пронаблюдать изменение частоты от 0 до 200 МГц за 5 мкс, что соответствует заданным требованиям по девиации частоты `bw` и длительности импульса `tWave` генератора

3. Канал распространения

В рассматриваемом сценарии моделью канала распространения сигнала является свободное пространство между каждым РЛС и целью в условиях отсутствия взаимных, пассивных и активных радиолокационных помех. Для реализации данной модели воспользуемся системным объектом EngeePhased.FreeSpace, позволяющий учесть затухание при распространении на несущей частоте, доплеровское смещение при наличии относительного движения цели и двунаправленное распространение (от цели и обратно).

# Формирование канала распространения сигнала с учетом двунаправленного распространения
channel = EngeePhased.FreeSpace( # СО - модель канала (свободное пространство)
    PropagationSpeed=c, # скорость распространения
    OperatingFrequency=fc, # несущая частота
    SampleRate=fs, # частота дискретизации
    TwoWayPropagation=true # учет двунаправленного распространения
)

FreeSpace:
    PropagationSpeed=299792458
    OperatingFrequency=9.0e9
    SampleRate=2.0e8
    TwoWayPropagation=true
    MaximumDistance=10000.0

4. Расчет сценария работы активной системы РЛС

После инициализации системных объектов, необходимых для расчета сценария работы системы, выполним расчет отраженного сигнала от БПЛА для каждой РЛС. Суммарный отраженный сигнал будем хранить в переменной X с размерностью [N,M,R]:

N — количество отсчётов в быстром времени,

M — количество импульсов в медленном времени.

R — количество РЛС в активной системе.

# Создание опорного сигнала
refsig = deepcopy(sig_lfm[:,1])

# Выделение памяти под отраженный сигнал для всех РЛС
X = zeros(ComplexF64,size(sig_lfm)...,numRadar)

# Сигнал после передатчика
txsig = transmitter.(sig_lfm)

for rad_i in 1:numRadar
    # Инициализация входного сигнала 
    x = zeros(ComplexF64,size(sig_lfm))

    # Приемо-передающий тракт
    for m in 1:M
        # Обновление положения цели и РЛС
        radar_pos,radar_vel = radarplatform(tWave)
        tgt_pos,tgt_vel = tgtplatform(tWave)

        # Расчет углов пеленга передающего устройства
        _,txang = rangeangle(tgt_pos,radar_pos[:,rad_i]) 

        # Излучение ЗС в пространство
        radtxsig = radiator(txsig[:,m],txang) 

        # Распространение сигнала в пространстве
        chansig = channel(
            radtxsig,radar_pos[:,rad_i],
            tgt_pos,radar_vel[:,rad_i],tgt_vel
        )

        # Отражение сигнала от цели
        tgtsig = target(chansig)

        # Расчет углов пеленга приемного устройства
        _,rxang = rangeangle(radar_pos[:,rad_i],tgt_pos)

        # Прием отраженного сигнала антенной
        rxsig = collector(tgtsig,rxang)

        # Предварительное усиление приемника
        x[:,m] .= receiver(rxsig)
    end

    dechirpsig = dechirp(x,refsig)
    X[:,:,rad_i] = conj.(dechirpsig)

    # Сброс положения РЛС и цели для расчета отраженного для следующей  РЛС
    reset_plt!(radarplatform,radarpos)
    reset_plt!(tgtplatform,tgtpos)
end

5. Расчет оценок задержки сигнала и локализация цели

После расчета отраженного сигнала для каждой РЛС следующим шагом является получение измерений времени прихода сигнала. Для этого воспользуемся объект TOAEstimator для оценки времени прихода сигнала, настроив Measurement="TOA". Метод спектрального анализа можно настроить как на FFT, так и на MUSIC.

spectrumMethod = "FFT"; # @param ["FFT", "MUSIC"]

# Формирование объекта TOAEstimator
toaEstimator = TOAEstimator(
    PropagationSpeed=c, # скорость распространения
    Measurement="TOA", # метод локального позиционирования 
    SpectrumMethod=spectrumMethod,  # спектральный метод
    VarianceOutputPort=true, # включение выхода СКО ошибки задержки
    SpatialSmoothing=ceil(Int64,N/2) # для алгоритм MUSIC
);

Выполним вызов объекта "toaEstimator" для расчета вектора задержек Y1 и вектора дисперсии var1

# Расчет оценок задержек и дисперсии
Y1,toa_var1 = toaEstimator(X,freqSpacing)
println("Оценки задержек для каждой РЛС: $(round.(Y1.*1e9;sigdigits=6)) нc")

Оценки задержек для каждой РЛС: [235.0, 500.0, 575.0, 365.0, 515.0] нc

# Построение спектра отраженного сигнала с пересчетом в задержку распространения
plotTDOASpectrum(
    toaEstimator, # объект расчета задержек
    anchorIdx=1:5, # Номера РЛС для отображения
    MaxDelay=700e-9, # Предел по оси задержек
    DinRange=30, # Динамический диапазон СПМ
    OneSidedSpectrum=true # построение одностороннего спектра
)

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

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

# Корректировка с учетом двунаправленного распространения
Y_single = Y1 ./ 2 # задержка, с
var_single = toa_var1 ./ 4 # дисперсия, с^2

# Расчет оценки положения цели
tgtposest = toaposest(Y_single,var_single,radarpos)

println("Оценка положений $(round.(tgtposest;sigdigits=6)) м")
println("Истинное положение $(tgtpos) м")

Оценка положений [30.1159, 9.97408, 15.2852] м
Истинное положение [30, 10, 15] м

Ниже приведена визуализация процесса локализации методом TOA:

helperPlotTOAPositions(c,Y_single,tgtposest,radarpos,tgtpos)

Мы также можем проверить точность определения местоположения с помощью среднеквадратичной ошибки (RMSE).

# Расчет СКО позиционирования
RMSE = rmse(tgtposest,tgtpos)
println("СКО позиционирования: $(round(RMSE;sigdigits=6)) м")

СКО позиционирования: 0.178334 м

Заключение

В примере продемонстрировано решение задачи локализации малоразмерной цели с применением активной многопозиционной РЛС в среде моделирования Engee. .

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

Локализация целей на основе активных РЛС