Введение

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

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

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

Setup сценария

Предположим, что наземный радар наблюдения работает в полосе L и передает щебетание сигналов 3 $$\mu s$$ длительность в интервале повторения импульса 15 $$\mu s$$. Полоса пропускания переданного щебета составляет 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);

Передающая антенна радара является равномерной прямоугольной фазированной решеткой 8x8, имеющей интервал $$\lambda$$/ 2 между его элементами. Сигнал переходит от радара к самолету и перехватывается и анализируется RWR. Для простоты форма волны выбрана в виде линейной FM волны с пиковой степенью 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 оборудован равномерным прямоугольным массивом 10x10 с интервалом $$\lambda$$/ 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 $$\mu s$$ сегменты.

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

Поскольку RWR не знает точной центральной частоты, используемой в форме волны передачи, он сначала использует банк фильтров, каждый из которых настроен на слегка разную частоту RF-центра, чтобы разделить принятые данные на поддиапазоны. Затем огибающий детектор прикладывается в каждой полосе, чтобы проверить, присутствует ли сигнал. В этом примере сигнал разделяется на поддиапазоны пропускной способности 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 приемных элементов. Поскольку степень передачи низкая, а шум приемника высок, радиолокационный сигнал неотличим от шума. Поэтому мощность приемника суммируется между этими элементами, чтобы улучшить отношение сигнал/нос (ОСШ) и получить лучшие оценки степени в каждом поддиапазоне. Полоса, которая имеет максимальную степень, является той, которая используется радаром.

% 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

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

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;

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

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

Предположим, что плоскость заземления плоская и параллельна плоскости 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)));

После применения рулевого управления лучом антенна имеет максимальное усиление в азимуте и угле возвышения прихода сигнала. Это дополнительно улучшает ОСШ перехваченного сигнала. Затем параметры сигнала извлекают в сигнальном процессоре с помощью одного из методов частотно-временного анализа, известных как псевдослучайное преобразование Вигнера-Вилля, связанное с преобразованием Хафа, как описано в [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')

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

% 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] Electronic Warfare and Радиолокационные Системы Engineering Handbook 2013, Naval Air Warfare Center Weapons Division, Пойнт-Мугу, Калифорния.

[2] Дэниел Л. Стивенс, Стефани А. Шукерс, обнаружение и экстракция параметров малой вероятности точки пересечения радиолокационных сигналов с помощью преобразования Хафа. Global Journal of Research In Engineering Vol 15 Выпуск 6, Январь 2016