Введение

RWR - это пассивная система поддержки радиоэлектронной борьбы [1], которая предоставляет летчику своевременную информацию об окружении его радиочастотного сигнала. RWR перехватывает ударный сигнал и использует методы обработки сигнала для извлечения информации о характеристиках перехваченного сигнала, а также о местоположении эмиттера. Эта информация может использоваться для вызова контрмер, таких как помехи, чтобы избежать обнаружения радиолокатором. Взаимодействие РЛС и самолета показано на следующей схеме.

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

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

Настройка сценария

Предположим, что наземная радиолокационная станция наблюдения работает в диапазоне L и передает чирп-сигналы $$$$длительностью 3 мкс с интервалом повторения импульсов 15 $$\mathrm{мкс}$$. Полоса пропускания передаваемого чирпа составляет 30 МГц, а несущая частота - 1,8 ГГц. РЛС наблюдения расположена у истока и неподвижна, а летательный аппарат летит с постоянной скоростью 200 м/с (~ 0,6 Маха).

% Define the transmitted waveform parameters
fs = 4e9;                     % Sampling frequency for the systems (Hz)
fc = 1.8e9;                   % Operating frequency of the surveillance radar (Hz)
T = 3e-6;                     % Chirp duration (s)
PRF = 1/(15e-6);              % Pulse repetition frequency (Hz)
BW = 30e6;                    % Chirp bandwidth (Hz)
c= physconst('LightSpeed');   % Speed of light in air (m/s)

% Assume the surveillance radar is at the origin and is stationary
radarPos= [0;0;0];          % Radar position (m)
radarVel= [0;0;0];          % Radar speed (m/s)

% Assume aircraft is moving with constant velocity
rwrPos= [-3000;1000;1000];   % Aircraft position (m)
rwrVel= [200; 0; 0];        % Aircraft speed (m/s)

% Configure objects to model ground radar and aircraft's relative motion
rwrPose = phased.Platform(rwrPos, rwrVel);
radarPose = phased.Platform(radarPos, radarVel);

Передающая антенна РЛС представляет собой однородную прямоугольную фазированную решетку 8х8, имеющую промежуток $$$$λ/2 между её элементами. Сигнал распространяется от РЛС на самолет и перехватывается и анализируется RWR. Для простоты форму сигнала выбирают в виде линейного ЧМ сигнала с пиковой мощностью 100 Вт.

% Configure the LFM waveform using the waveform parameters defined above
wavGen= phased.LinearFMWaveform('SampleRate',fs,'PulseWidth',T,'SweepBandwidth',BW,'PRF',PRF);

% Configure the Uniform Rectangular Array
antennaTx = phased.URA('ElementSpacing',repmat((c/fc)/2, 1, 2), 'Size', [8,8]);

% Configure objects for transmitting and propagating the radar signal
tx = phased.Transmitter('Gain', 5, 'PeakPower',100);
radiator = phased.Radiator( 'Sensor', antennaTx, 'OperatingFrequency', fc);
envIn = phased.FreeSpace('TwoWayPropagation',false,'SampleRate', fs,'OperatingFrequency',fc);

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

% Transmit a train of pulses
numPulses = 4;
txPulseTrain = helperRWR('simulateTransmission',numPulses, wavGen, rwrPos,...
    radarPos, rwrVel, radarVel, rwrPose, radarPose, tx, radiator, envIn,fs,fc,PRF);

% Observe the signal arriving at the RWR
pspectrum(txPulseTrain,fs,'spectrogram','FrequencyLimits',[1.7e9 1.9e9], 'Leakage',0.65)
title('Transmitted pulse train spectrogram'); caxis([-110 -90]);

RWR оснащена однородной прямоугольной матрицей 10х10 с шагом $$$$λ/2 между её элементами. Работает во всём L-диапазоне, с центральной частотой 2 ГГц. RWR прослушивает среду и непрерывно подает собранные данные в цепочку обработки.

% Configure the receive antenna
dip = phased.IsotropicAntennaElement('BackBaffled',true);
antennaRx = phased.URA('ElementSpacing',repmat((c/2e9)/2,1,2),'Size', [10,10],'Element',dip);

% Model the radar receiver chain
collector = phased.Collector('Sensor', antennaRx,'OperatingFrequency',fc);
rx = phased.ReceiverPreamp('Gain',0,'NoiseMethod','Noise power', 'NoisePower',2.5e-6,'SeedSource','Property', 'Seed',2018);

% Collect the waves at the receiver
[~, tgtAng] = rangeangle(radarPos,rwrPos);
yr = collector(txPulseTrain,tgtAng);
yr = rx(yr);

Детектор огибающей RWR

Детектор огибающей в RWR отвечает за обнаружение наличия любого сигнала. Поскольку RWR непрерывно принимает данные, цепочка приемника буферизирует и усекает принятые данные в $$$$сегменты 50 мкс.

% Truncate the received data
truncTime = 50e-6;
truncInd = round(truncTime*fs);
yr = yr(1:truncInd,:);

Поскольку RWR не знает о точной центральной частоте, используемой в сигнале передачи, он сначала использует набор фильтров, каждый из которых настроен на немного отличающуюся центральную частоту РЧ, для разделения принятых данных на поддиапазоны. Затем детектор огибающей применяется в каждой полосе, чтобы проверить, присутствует ли сигнал. В этом примере сигнал делится на поддиапазоны шириной 100 МГц. Дополнительным преимуществом такой операции является то, что вместо дискретизации всей полосы пропускания, покрываемой RWR, сигнал в каждом поддиапазоне может быть дискретизирован с понижением частоты до частоты дискретизации 100 МГц.

% Define the bandwidth of each frequency sub-band
stepFreq = 100e6;

% Calculate number of sub-bands and configure dsp.Channelizer
numChan = fs/stepFreq;
channelizer = dsp.Channelizer('NumFrequencyBands', numChan, 'StopbandAttenuation', 80);

На графике ниже показаны первые четыре полосы, созданные банком фильтров.

% Visualize the first four filters created in the filter bank of the
% channelizer
freqz(channelizer, 1:4)
title('Zoomed Channelizer response for first four filters')
xlim([0 0.2])

% Pass the received data through the channelizer
subData = channelizer(yr);

Полученные данные, subData, имеет 3 измерения. Первое измерение представляет быстрое время, второе измерение представляет поддиапазоны, и третье измерение соответствует принимающим элементам принимающего массива. Для конфигурации антенны RWR 10x10, используемой в этом примере, имеется 100 приемных элементов. Поскольку мощность передачи низкая, а шум приемника высокий, радиолокационный сигнал неотличим от шума. Поэтому принимаемая мощность суммируется по этим элементам для улучшения отношения сигнал/нос (SNR) и получения лучших оценок мощности в каждом поддиапазоне. Диапазон, который имеет максимальную мощность, используется радаром.

% Rearrange the subData to combine the antenna array channels only
incohsubData = pulsint(permute(subData,[1,3,2]),'noncoherent');
incohsubData = squeeze(incohsubData); 

% Plot power distribution
subbandPow = pow2db(rms(incohsubData,1).^2)+30;
plot(subbandPow);
xlabel('Band Index');
ylabel('Power (dBm)');

% Find the sub-band with maximum power
[~,detInd] = max(subbandPow);

Процессор сигналов RWR

Хотя мощность в выбранном диапазоне выше по сравнению с соседним диапазоном, SNR в диапазоне все еще остается низким, как показано на следующем рисунке.

subData = (subData(:,detInd,:));
subData = squeeze(subData);  %adjust the data to 2-D matrix

% Visualize the detected sub-band data
plot(mag2db(abs(sum(subData,2)))+30)
ylabel('Power (dBm)')
title('Detected sub-band from 100 channels combined incoherently')

% Find the original starting frequency of the sub-band having the detected
% signal
detfBand = fs*(detInd-1)/(fs/stepFreq);

% Update the sampling frequency to the decimated frequency
fs = stepFreq;

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

Следующим шагом для RWR является поиск направления, из которого поступают радиоволны. Эта информация об угле прихода будет использоваться для направления принимающего антенного луча в направлении излучателя и для определения местоположения излучателя на земле с использованием однолинейного подхода. RWR оценивает направление поступления, используя двухмерный модуль оценки MUSIC. Управление лучом осуществляется с использованием фазовращающего формирователя луча для достижения максимального SNR сигнала, таким образом помогая извлечению параметра формы сигнала.

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

% Configure the MUSIC Estimator to find the direction of arrival of the
% signal
doaEst = phased.MUSICEstimator2D('OperatingFrequency',fc,'PropagationSpeed',c,...
    'SensorArray',antennaRx,'DOAOutputPort',true,'AzimuthScanAngles',-50:.5:50,...
    'ElevationScanAngles',-50:.5:50, 'NumSignalsSource', 'Property','NumSignals', 1);

[mSpec,doa] = doaEst(subData);
plotSpectrum(doaEst,'Title','2-D MUSIC Spatial Spectrum Top view');
view(0,90); axis([-30 0 -30 0]);

На рисунке четко показано расположение эмиттера.

% Configure the beamformer object to steer the beam before combining the
% channels
beamformer = phased.PhaseShiftBeamformer('SensorArray',antennaRx,...
    'OperatingFrequency',fc,'DirectionSource','Input port');

% Apply the beamforming, and visualize the beam steered radiation
% pattern
mBeamf = beamformer(subData, doa);

% Find the location of the emitter
altimeterElev = rwrPos(3);
d = abs(altimeterElev/sind(doa(2)));

После применения управления лучом антенна имеет максимальный коэффициент усиления по азимуту и углу места прихода сигнала. Это дополнительно улучшает SNR перехваченного сигнала. Затем параметры сигнала извлекают в процессоре сигналов, используя один из методов частотно-временного анализа, известный как псевдо-преобразование Вигнера-Вилля, связанное с преобразованием Хафа, как описано в [2].

Сначала получают представление временной частоты перехваченного сигнала с использованием преобразования Вигнера-Вилля.

% Compute the pseudo Wigner-Ville transform
[tpwv,t,f] = helperRWR('pWignerVille',mBeamf,fs);

% Plot the pseudo Wigner-Ville transform
imagesc(f*1e-6,t*1e6,pow2db(abs(tpwv./max(tpwv(:)))));
xlabel('Frequency (MHz)'); ylabel('Time(\mus)');
caxis([-50 0]); clb = colorbar; clb.Label.String = 'Normalized Power (dB)';
title ('Pseudo Wigner-Ville Transform')

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

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

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

% Normalize the pseudo Wigner-Ville image
twvNorm = abs(tpwv)./max(abs(tpwv(:)));

% Implement a median filter to clear the noise
filImag = medfilt2(twvNorm,[7 7]);

% Use threshold to convert filtered image into binary image
BW = imbinarize(filImag./max(filImag(:)), 0.15);
imagesc(f*1e-6,t*1e6,BW); colormap('gray');
xlabel('Frequency (MHz)'); ylabel('Time(\mus)');
title ('Pseudo Wigner-Ville Transform - BW')

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

% Compute the Hough transform of the image and plot
[H,T,R] = hough(BW);
imshow(H,[],'XData',T,'YData',R,'InitialMagnification','fit');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
title('Hough transform of the image')

Пиковые положения извлекаются с помощью houghpeaks.

% Compute peaks in the transform, up to 5 peaks
P  = houghpeaks(H,5);
x = T(P(:,2)); y = R(P(:,1));
plot(x,y,'s','color','g'); xlim([-90 -50]); ylim([-5000 0])

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

lines = houghlines(BW,T,R,P,'FillGap',3e-6*fs,'MinLength',1e-6*fs);
coord = [lines(:).point1; lines(:).point2];

% Plot the detected lines superimposed on the binary image
clf;
imagesc(f*1e-6, t*1e6, BW); colormap(gray); hold on
xlabel('Frequency (MHz)')
ylabel('Time(\mus)')
title('Hough transform - detected lines')
for ii = 1:2:2*size(lines,2)
    plot(f(coord(:,ii))*1e-6, t(coord(:,ii+1))*1e6,'LineWidth',2,'Color','green');
end

% Calculate the parameters using the line co-ordinates
pulDur = t(coord(2,2)) - t(coord(1,2));         % Pulse duration
bWidth = f(coord(2,1)) - f(coord(1,1));         % Pulse Bandwidth
pulRI = abs(t(coord(1,4)) - t(coord(1,2)));     % Pulse repetition interval
detFc = detfBand + f(coord(2,1));               % Center frequency

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

helperRWR('displayParameters',pulRI, pulDur, bWidth,detFc, doa,d);

Pulse repetition interval = 14.97 microseconds
Pulse duration = 2.84 microseconds
Pulse bandwidth = 27 MHz
Center frequency = 1.8286 GHz
Azimuth angle of emitter = -18.5 degrees
Elevation angle of emitter = -17.5 degrees
Distance of the emitter = 3325.5095 m

Резюме

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

Ссылки

[1] Руководство по радиоэлектронной войне и проектированию радиолокационных систем 2013 года, Отдел вооружений Военно-морского центра воздушной войны, Пойнт-Мугу, Калифорния.

[2] Дэниел Л. Стивенс (Daniel L. Stevens), Стефани А. Шукерс (Stephanie A. Schuckers), обнаружение и выделение параметров низкой вероятности перехвата радиолокационных сигналов с использованием преобразования Хафа. Global Journal of Research in Engineering Vol 15 Выпуск 6, январь 2016 г.