Введение в пространственно-временную адаптивную обработку

Этот пример дает краткое введение в методы пространственно-временной адаптивной обработки (STAP) и иллюстрирует, как использовать Toolbox™ Phased Array System для применения алгоритмов STAP к полученным импульсам. STAP - метод, используемый в бортовых радиолокационных системах для подавления помех от загромождения и помех от помех.

Введение

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

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

В традиционных системах MTI фильтрация загромождений часто использует тот факт, что земля не перемещается. Таким образом, загромождение занимает нулевой доплеровский интервал в допплеровском спектре. Этот принцип приводит ко многим методам фильтрации загромождения на основе Доплера, таким как подавитель импульсов. Заинтересованные читатели могут обратиться к радару Ground Clutter Mitigation with Moving Target Indication (MTI) для подробного примера подавителя импульсов. Когда сама радиолокационная платформа также движется, например, в плоскости, доплеровский компонент от возврата земли больше не равен нулю. В сложение доплеровские компоненты возвратов загромождения зависят от угла. В этом случае возврат загромождения, вероятно, будет иметь энергию через спектр Доплера. Поэтому загромождение не может быть отфильтровано только относительно доплеровской частоты.

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

Методы STAP фильтруют сигнал как в угловом, так и в Допплеровском областях (таким образом, имя «пространственно-временная адаптивная обработка»), чтобы подавить возвраты загромождения и глушителя. В следующих разделах мы моделируем возвраты от цели, загромождения и глушителя и иллюстрируем, как методы STAP фильтруют помехи от принятого сигнала.

Setup системы

Сначала зададим радиолокационную систему, начиная с системы, построенной в примере Simulation Test Signals для радиолокационного приемника.

load BasicMonostaticRadarExampleData.mat;     % Load monostatic pulse radar

Определение антенны

Предположим, что антенный элемент имеет изотропный ответ в передней полусфере и все нули в задней полусфере. Рабочая частота области значений 8-12 ГГц, чтобы соответствовать рабочей частоте системы 10 ГГц.

antenna = phased.IsotropicAntennaElement...
    ('FrequencyRange',[8e9 12e9],'BackBaffled',true); % Baffled Isotropic

Задайте ULA с 6 элементами с пользовательским шаблоном элемента. Интервал между элементами принимается равным половине длины волны формы волны.

fc = radiator.OperatingFrequency;
c = radiator.PropagationSpeed;
lambda = c/fc;
ula = phased.ULA('Element',antenna,'NumElements',6,...
    'ElementSpacing', lambda/2);

pattern(ula,fc,'PropagationSpeed',c,'Type','powerdb')
title('6-element Baffled ULA Response Pattern')
view(60,50)

Радиолокационный Setup

Затем установите антенную решетку на радиатор/коллектор. Затем задайте радиолокационное движение. Радиолокационная система установлена на плоскости, которая пролетает 1000 метров над землей. Плоскость летит вдоль оси массива ULA со скоростью, такой что она перемещается с интервалом в половину элемента матрицы в течение одного импульсного интервала. (Объяснение такой настройки приведено в следующем разделе метода DPCA.)

radiator.Sensor = ula;
collector.Sensor = ula;
sensormotion = phased.Platform('InitialPosition',[0; 0; 1000]);
arrayAxis = [0; 1; 0];
prf = waveform.PRF;
vr = ula.ElementSpacing*prf;  % in [m/s]
sensormotion.Velocity = vr/2*arrayAxis;

Цель

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

target = phased.RadarTarget('Model','Nonfluctuating','MeanRCS',1, ...
    'OperatingFrequency', fc);
tgtmotion = phased.Platform('InitialPosition',[1000; 1000; 0],...
    'Velocity',[30; 30; 0]);

Глушитель

Цель возвращает требуемый сигнал; однако в принимаемом сигнале также присутствуют несколько помех. Если Radar Toolbox доступен, пожалуйста, установите переменную hasRadarToolbox для определения простого заграждающего глушителя с эффективной излучаемой степенью 100 Вт. В противном случае в симуоляции будет использоваться сохраненный сигнал глушителя.

hasRadarToolbox = false;
Fs = waveform.SampleRate;
rngbin = c/2*(0:1/Fs:1/prf-1/Fs).';
if hasRadarToolbox
    jammer = barrageJammer('ERP',100);
    jammer.SamplesPerFrame = numel(rngbin);
    jammermotion = phased.Platform('InitialPosition',[1000; 1732; 1000]);
end

Беспорядок

В этом примере мы симулируем загромождение, используя постоянную гамма- модель с гамма- значением -15 дБ. Литература показывает, что такое значение гаммы может использоваться для моделирования местности, покрытой лесом. Для каждой области значений возврат загромождения может рассматриваться как комбинация возвратов из многих небольших загроможденных закрашенных фигур на этой области значений звонка. Поскольку антенна находится в заслонке, вклад загромождения только спереди. Чтобы упростить расчет, используйте ширину азимута 10 степеней для каждой закрашенной фигуры. Снова, если Radar Toolbox недоступен, в симуляции будет использован сохраненный сигнал загромождения.

if hasRadarToolbox
    mountingAng = [30,0,0];
    Rmax = 5000;
    Azcov = 120;
    clutter = constantGammaClutter('Sensor',ula,'SampleRate',Fs,...
        'Gamma',-15,'PlatformHeight',1000,...
        'OperatingFrequency',fc,...
        'PropagationSpeed',c,...
        'PRF',prf,...
        'TransmitERP',transmitter.PeakPower*db2pow(transmitter.Gain),...
        'PlatformSpeed',norm(sensormotion.Velocity),...
        'PlatformDirection',[90;0],...
        'MountingAngles',mountingAng,'ClutterAzimuthCenter', ...
        'ClutterMaxRange', Rmax,...
        'ClutterAzimuthSpan',Azcov,...
        'PatchAzimuthSpan',10,...
        'OutputFormat','Pulses');
end

Пути распространения

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

tgtchannel = phased.FreeSpace('TwoWayPropagation',true,'SampleRate',Fs,...
    'OperatingFrequency', fc); 
jammerchannel = phased.FreeSpace('TwoWayPropagation',false,...
    'SampleRate',Fs,'OperatingFrequency', fc);

Область Цикла

Теперь мы готовы моделировать возвраты. Собирайте 10 импульсов перед обработкой. Seed генератора случайных чисел из модели jammer устанавливается на константу, чтобы получить воспроизводимые результаты.

numpulse = 10; % Number of pulses
tsig = zeros(size(rngbin,1), ula.NumElements, numpulse);
jsig = tsig; tjcsig = tsig; tcsig = tsig; csig = tsig; 

if hasRadarToolbox
    jammer.SeedSource = 'Property';
    jammer.Seed = 5;
    clutter.SeedSource = 'Property';
    clutter.Seed = 5;
else
    load STAPIntroExampleData;
end

for m = 1:numpulse
    
    % Update sensor, target and calculate target angle as seen by the sensor
    [sensorpos,sensorvel] = sensormotion(1/prf);          
    [tgtpos,tgtvel] = tgtmotion(1/prf);          
    [~,tgtang] = rangeangle(tgtpos,sensorpos);       

    % Update jammer and calculate the jammer angles as seen by the sensor
    if hasRadarToolbox
        [jampos,jamvel] = jammermotion(1/prf);  
        [~,jamang] = rangeangle(jampos,sensorpos);    
    end
    
    % Simulate propagation of pulse in direction of targets
    pulse = waveform();
    [pulse,txstatus] = transmitter(pulse);
    pulse = radiator(pulse,tgtang);
    pulse = tgtchannel(pulse,sensorpos,tgtpos,sensorvel,tgtvel);
    
    % Collect target returns at sensor
    pulse = target(pulse);
    tsig(:,:,m) = collector(pulse,tgtang);
    
    % Collect jammer and clutter signal at sensor
    if hasRadarToolbox
        jamsig = jammer();
        jamsig = jammerchannel(jamsig,jampos,sensorpos,jamvel,sensorvel);
        jsig(:,:,m) = collector(jamsig,jamang);
        
        csig(:,:,m) = clutter();
    end
    
    % Receive collected signals
    tjcsig(:,:,m) = receiver(tsig(:,:,m)+jsig(:,:,m)+csig(:,:,m),...
        ~(txstatus>0));                         % Target + jammer + clutter
    tcsig(:,:,m) = receiver(tsig(:,:,m)+csig(:,:,m),...
        ~(txstatus>0));                         % Target + clutter
    tsig(:,:,m) = receiver(tsig(:,:,m),...
        ~(txstatus>0));                         % Target echo only
end

Истинная целевая область значений, угол и допплер

Целевой угол азимута составляет 45 степени, а угол возвышения - около -35,27 степени.

tgtLocation = global2localcoord(tgtpos,'rs',sensorpos);
tgtAzAngle = tgtLocation(1)
tgtAzAngle = 44.9981
tgtElAngle = tgtLocation(2)
tgtElAngle = -35.2651

Целевая область значений - 1732 м.

tgtRng = tgtLocation(3)
tgtRng = 1.7320e+03

Целевая доплеровская нормированная частота составляет около 0,21.

sp = radialspeed(tgtpos, tgtmotion.Velocity, ...
                sensorpos, sensormotion.Velocity);
tgtDp = 2*speed2dop(sp,lambda);  % Round trip Doppler
tgtDp/prf
ans = 0.2116

Общий полученный сигнал содержит возвраты от цели, загромождения и глушителя вместе взятых. Сигнал является кубом данных с тремя размерностями (области значений x количество элементов x количество импульсов). Заметьте, что возврат загромождения доминирует над общим возвращением и маскирует возврат целевого значения. Мы не можем обнаружить цель (синяя вертикальная линия) без дальнейшей обработки на данном этапе.

ReceivePulse = tjcsig;
plot([tgtRng tgtRng],[0 0.01],rngbin,abs(ReceivePulse(:,:,1)));
xlabel('Range (m)'), ylabel('Magnitude');
title('Signals collected by the ULA within the first pulse interval')

Figure contains an axes. The axes with title Signals collected by the ULA within the first pulse interval contains 7 objects of type line.

Теперь исследуйте возвраты в 2-D угловой Доплеровской (или пространственно-временной) области. В целом, реакция генерируется путем сканирования всех областей значений и углов азимута для заданного угла возвышения. Поскольку мы точно знаем, где находится цель, мы можем вычислить ее область значений и угол возвышения относительно антенной решетки.

tgtCellIdx = val2ind(tgtRng,c/(2*Fs));
snapshot = shiftdim(ReceivePulse(tgtCellIdx,:,:));  % Remove singleton dim
angdopresp = phased.AngleDopplerResponse('SensorArray',ula,...
              'OperatingFrequency',fc, 'PropagationSpeed',c,...
              'PRF',prf, 'ElevationAngle',tgtElAngle);
plotResponse(angdopresp,snapshot,'NormalizeDoppler',true);
text(tgtAzAngle,tgtDp/prf,'+ Target')

Figure contains an axes. The axes with title Angle-Doppler Response Pattern contains 2 objects of type image, text.

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

Подавление загромождения с помощью компенсатора DPCA

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

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

Предположим, что N является количеством элементов ULA. Возврат загромождения, принимаемый в антенне 1 через антенну N-1 во время первого импульса, аналогичен возврату загромождения, принимаемому в антенне 2 через антенну N во время второго импульса. Путем вычитания импульсов, принятых в этих двух подрешетках, в течение двух импульсных интервалов, загромождение может быть отменено. Стоимость этого метода является апертурой, которая на один элемент меньше исходного массива.

Теперь определите отменщик DPCA. Алгоритму может потребоваться поиск по всем комбинациям угла и Допплера, чтобы найти цель, но для примера, потому что мы точно знаем, где цель, мы можем направить наш процессор в эту точку.

rxmainlobedir = [0; 0];
stapdpca = phased.DPCACanceller('SensorArray',ula,'PRF',prf,...
    'PropagationSpeed',c,'OperatingFrequency',fc,...
    'Direction',rxmainlobedir,'Doppler',tgtDp,...
    'WeightsOutputPort',true)
stapdpca = 
  phased.DPCACanceller with properties:

            SensorArray: [1x1 phased.ULA]
       PropagationSpeed: 299792458
     OperatingFrequency: 1.0000e+10
              PRFSource: 'Property'
                    PRF: 2.9979e+04
        DirectionSource: 'Property'
              Direction: [2x1 double]
    NumPhaseShifterBits: 0
          DopplerSource: 'Property'
                Doppler: 6.3429e+03
      WeightsOutputPort: true
       PreDopplerOutput: false

Сначала примените отменщик DPCA как к целевому возврату, так и к возврату загромождения.

ReceivePulse = tcsig;
[y,w] = stapdpca(ReceivePulse,tgtCellIdx);

Обработанные данные объединяют всю информацию в пространстве и через импульсы, чтобы стать одним импульсом. Затем исследуйте обработанный сигнал во временном интервале.

plot([tgtRng tgtRng],[0 1.2e-5],rngbin,abs(y));
xlabel('Range (m)'), ylabel('Magnitude');
title('DPCA Canceller Output (no Jammer)')

Figure contains an axes. The axes with title DPCA Canceller Output (no Jammer) contains 2 objects of type line.

Сигнал теперь четко отличается от шума, и загромождение отфильтровано. Из угла Доплеровской характеристики весов процессора DPCA ниже, мы также видим, что веса дают глубокое ядро вдоль хребта загромождения.

angdopresp.ElevationAngle = 0;
plotResponse(angdopresp,w,'NormalizeDoppler',true);
title('DPCA Weights Angle Doppler Response at 0 degrees Elevation')

Figure contains an axes. The axes with title DPCA Weights Angle Doppler Response at 0 degrees Elevation contains an object of type image.

Хотя результаты, полученные DPCA, очень хороши, радиолокационная платформа должна удовлетворять очень строгим требованиям в своем движении, чтобы использовать этот метод. Кроме того, метод DPCA не может подавить помехи джаммера.

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

ReceivePulse = tjcsig;
[y,w] = stapdpca(ReceivePulse,tgtCellIdx);
plot([tgtRng tgtRng],[0 8e-4],rngbin,abs(y));
xlabel('Range (m)'), ylabel('Magnitude');
title('DPCA Canceller Output (with Jammer)')

Figure contains an axes. The axes with title DPCA Canceller Output (with Jammer) contains 2 objects of type line.

plotResponse(angdopresp,w,'NormalizeDoppler',true);
title('DPCA Weights Angle Doppler Response at 0 degrees Elevation')

Figure contains an axes. The axes with title DPCA Weights Angle Doppler Response at 0 degrees Elevation contains an object of type image.

Загромождение и подавление помех с помощью SMI Beamformer

Чтобы подавить загромождение и глушитель одновременно, нам нужен более сложный алгоритм. Оптимальные веса приемника, когда интерференция распределена Гауссом, заданы [1]

w=kR-1s

где k - скалярный коэффициент, R - пространственно-временная ковариационная матрица интерференционного сигнала, и s - желаемый пространственно-временной вектор управления. Точная информация R часто недоступна, поэтому мы будем использовать алгоритм sample matrix inversion (SMI). Алгоритм оценивает R из выборок обучающих ячеек и затем использует его в вышеупомянутом уравнении.

Теперь задайте SMI-формирователь луча и примените его к сигналу. В дополнение к информации, необходимой в DPCA, светоформер SMI должен знать количество защитных камер и количество обучающих камер. Алгоритм использует выборки в обучающих камерах, чтобы оценить интерференцию. Таким образом, мы не должны использовать камеры, которые близки к целевой камере, для оценок, потому что они могут содержать некоторую информацию о цели, то есть мы должны задать защитные камеры. Количество защитных камер должно быть четным числом, которое должно быть разделено поровну перед и позади целевой камеры. Количество обучающих камер также должно быть четным числом и разделяться поровну перед и позади цели. Обычно, чем больше количество обучающих камер, тем лучше оценка интерференции.

tgtAngle = [tgtAzAngle; tgtElAngle];
stapsmi = phased.STAPSMIBeamformer('SensorArray', ula, 'PRF', prf, ...
    'PropagationSpeed', c, 'OperatingFrequency', fc, ...
    'Direction', tgtAngle, 'Doppler', tgtDp, ...
    'WeightsOutputPort', true,...
    'NumGuardCells', 4, 'NumTrainingCells', 100)
stapsmi = 
  phased.STAPSMIBeamformer with properties:

            SensorArray: [1x1 phased.ULA]
       PropagationSpeed: 299792458
     OperatingFrequency: 1.0000e+10
              PRFSource: 'Property'
                    PRF: 2.9979e+04
        DirectionSource: 'Property'
              Direction: [2x1 double]
    NumPhaseShifterBits: 0
          DopplerSource: 'Property'
                Doppler: 6.3429e+03
          NumGuardCells: 4
       NumTrainingCells: 100
      WeightsOutputPort: true

[y,w] = stapsmi(ReceivePulse,tgtCellIdx);
plot([tgtRng tgtRng],[0 2e-6],rngbin,abs(y));
xlabel('Range (m)'), ylabel('Magnitude');
title('SMI Beamformer Output (with Jammer)')

Figure contains an axes. The axes with title SMI Beamformer Output (with Jammer) contains 2 objects of type line.

plotResponse(angdopresp,w,'NormalizeDoppler',true);
title('SMI Weights Angle Doppler Response at 0 degrees Elevation')

Figure contains an axes. The axes with title SMI Weights Angle Doppler Response at 0 degrees Elevation contains an object of type image.

Результат показывает, что SMI-формирователь луча может отличать сигналы как от загромождителя, так и от глушителя. Угловой допплеровский шаблон весов SMI показывает глубокое ядро вдоль направления глушителя.

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

Сокращение вычислительных затрат с помощью Canceller ADPCA

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

Адаптивный компенсатор DPCA (ADPCA) отфильтровывает загромождение так же, как и DPCA, но он также имеет возможность подавлять глушитель, так как он оценивает ковариационную матрицу интерференции с помощью двух последовательных импульсов. Поскольку задействованы только два импульса, расчет значительно уменьшается. В сложение, поскольку алгоритм адаптирован к интерференции, он также может терпеть некоторые нарушения порядка движения.

Теперь определите компенсатор ADPCA, а затем примените его к полученному сигналу.

stapadpca = phased.ADPCACanceller('SensorArray', ula, 'PRF', prf, ...
    'PropagationSpeed', c, 'OperatingFrequency', fc, ...
    'Direction', rxmainlobedir, 'Doppler', tgtDp, ...
    'WeightsOutputPort', true,...
    'NumGuardCells', 4, 'NumTrainingCells', 100)
stapadpca = 
  phased.ADPCACanceller with properties:

            SensorArray: [1x1 phased.ULA]
       PropagationSpeed: 299792458
     OperatingFrequency: 1.0000e+10
              PRFSource: 'Property'
                    PRF: 2.9979e+04
        DirectionSource: 'Property'
              Direction: [2x1 double]
    NumPhaseShifterBits: 0
          DopplerSource: 'Property'
                Doppler: 6.3429e+03
          NumGuardCells: 4
       NumTrainingCells: 100
      WeightsOutputPort: true
       PreDopplerOutput: false

[y,w] = stapadpca(ReceivePulse,tgtCellIdx);
plot([tgtRng tgtRng],[0 2e-6],rngbin,abs(y));
xlabel('Range (m)'), ylabel('Magnitude');
title('ADPCA Canceller Output (with Jammer)')

Figure contains an axes. The axes with title ADPCA Canceller Output (with Jammer) contains 2 objects of type line.

plotResponse(angdopresp,w,'NormalizeDoppler',true);
title('ADPCA Weights Angle Doppler Response at 0 degrees Elevation')

Figure contains an axes. The axes with title ADPCA Weights Angle Doppler Response at 0 degrees Elevation contains an object of type image.

График временного интервала показывает, что сигнал успешно восстановлен. Угловая доплеровская характеристика весов ADPCA аналогична той, которая производится весами SMI.

Сводные данные

Этот пример представил краткое введение в пространственно-временную адаптивную обработку и проиллюстрировал, как использовать различные алгоритмы STAP, а именно SMI, DPCA и ADPCA, для подавления помех загромождения и помех помех в принятых импульсах.

Ссылка

[1] Дж. Р. Герси, пространственно-временная адаптивная обработка для радара, Artech House, 2003