Экстракция параметра формы волны от полученного импульса

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

Этот пример требует Image Processing Toolbox™

Введение

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

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

  1. Импульсный интервал повторения

  2. Центральная частота

  3. Пропускная способность

  4. Импульсная длительность

  5. Направление прибытия

  6. Положение эмиттера

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

Setup сценария

Примите, что наземный радар наблюдения действует в полосе L и передает сигналы щебета 3 μs длительность в импульсном интервале повторения 15 μ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 универсальная прямоугольная фазированная решетка, имея интервал λ/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]);

Figure contains an axes. The axes with title Transmitted pulse train spectrogram contains an object of type image.

RWR оборудован 10x10 универсальный прямоугольный массив с интервалом λ/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 μ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])

Figure contains an axes. The axes with title Zoomed Channelizer response for first four filters contains 4 objects of type line.

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

Полученные данные, subData, имеет 3 размерности. Первая размерность представляет быстро-разовое, второе измерение представляет поддиапазоны, и третья размерность соответствует элементам получения массива получения. Для RWR's 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)');

Figure contains an axes. The axes contains an object of type line.

% 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')

Figure contains an axes. The axes with title Detected sub-band from 100 channels combined incoherently contains an object of type line.

% 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. Регулирование луча сделано с помощью формирователя луча сдвига фазы, чтобы достигнуть максимального ОСШ сигнала, таким образом помочь экстракции параметра формы волны.

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

Figure contains an axes. The axes with title 2-D MUSIC Spatial Spectrum Top view contains an object of type surface.

Фигура ясно показывает местоположение эмиттера.

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

После применения регулирования луча антенна имеет максимальное усиление в азимуте и угле возвышения прибытия сигнала. Это далее улучшает ОСШ прерванного сигнала. Затем параметры сигнала извлечены в сигнальном процессоре с помощью одного из методов частотно-временного анализа, известных как псевдо Wigner-Ville, преобразовывают вместе с Преобразованием Хафа как описано в [2].

Во-первых, выведите представление частоты времени прерванного сигнала с помощью Wigner-Ville, преобразовывают.

% 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')

Figure contains an axes. The axes with title Pseudo Wigner-Ville Transform contains an object of type image.

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

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

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

Figure contains an axes. The axes with title Pseudo Wigner-Ville Transform - BW contains an object of type image.

Используя Преобразование Хафа, бинарное псевдо изображение Wigner-Ville сначала преобразовывается к 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])

Figure contains an axes. The axes with title Hough transform of the image contains 2 objects of type image, line.

Используя эти положения, 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

Figure contains an axes. The axes with title Hough transform - detected lines contains 4 objects of type image, line.

% 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] Дэниел Л. Стивенс, Стефани А. Шукерс, Обнаружение и Экстракция Параметра Низкой Вероятности Радарных Сигналов Точки пересечения с помощью Преобразования Хафа. Глобальный Журнал Исследования В Техническом выпуске 6 Vol 15, январь 2016