Обработка изображений радара с синтезированной апертурой

SAR генерирует двумерное (2D) изображение. Направление полета упоминается как направление азимута или перекрестная область значений. Направление опорного направления антенны (разворот) является ортогональным к курсу полета и упоминается как направление области значений или перекрестная дорожка. Эти два направления обеспечивают основание для размерностей, требуемых сгенерировать изображение, полученное из области в ширине луча антенны через длительность окна сбора данных. Направление перекрестной дорожки является направлением, в котором передаются импульсы. Это направление предоставляет наклонный спектр к целям вдоль курса полета. Энергия, полученная после отражения от целей для каждого импульса, должна затем быть обработана (для измерения области значений и разрешения). Направление перекрестной области значений или азимута является направлением курса полета, и это значимо, чтобы обработать ансамбль импульсов, полученных по целому курсу полета в этом направлении, чтобы достигнуть необходимого измерения и разрешения. Правильная фокусировка и в направлениях подразумевает успешную генерацию изображения в области значений и в направлениях перекрестной области значений. Это - требование для ширины луча антенны, чтобы быть достаточно широким так, чтобы цель была освещена на долгое время лучом, когда платформа проходит своя траектория. Это поможет предоставить больше информации о фазе. Ключевые условия, с которыми часто сталкиваются при работе с SAR:

Радарная настройка

Рассмотрите радар SAR, действующий в C-полосе с несущей частотой на 4 ГГц и пропускной способностью сигнала 50 МГц. Эта пропускная способность дает к разрешению области значений 3 метров. Радиолокационная система собирает данные, ортогональные к направлению движения платформы как показано в фигуре выше. Полученный сигнал является задержанной копией переданного сигнала. Задержка соответствует в целом наклонной области значений между целью и платформой. Для системы SAR наклонная область значений будет варьироваться в зависимости от времени, когда платформа пересекает путь, ортогональный к направлению луча антенны. Этот раздел ниже особого внимания при определении параметров для формы волны передачи. Пропускная способность развертки LFM может быть решена на основе желаемого разрешения области значений.

Установите физическую константу для скорости света.

c = physconst('LightSpeed');

Установите частоту центра SAR.

fc = 4e9;

Установите желаемую область значений и разрешение перекрестной области значений 3 метров.

rangeResolution = 3;  
crossRangeResolution = 3;

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

bw = c/(2*rangeResolution);

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

prf = 1000; 
aperture = 4;  
tpd = 3*10^-6; 
fs = 120*10^6;

Сконфигурируйте сигнал LFM радара.

waveform = phased.LinearFMWaveform('SampleRate',fs, 'PulseWidth', tpd, 'PRF', prf,...
    'SweepBandwidth', bw);

Примите, что скорость самолета составляет 100 м/с с длительностью рейса 4 секунд.

speed = 100;  
flightDuration = 4;
radarPlatform  = phased.Platform('InitialPosition', [0;-200;500], 'Velocity', [0; speed; 0]);
slowTime = 1/prf;
numpulses = flightDuration/slowTime +1;

maxRange = 2500;
truncrangesamples = ceil((2*maxRange/c)*fs);
fastTime = (0:1/fs:(truncrangesamples-1)/fs);
% Set the reference range for the cross-range processing.
Rc = 1000;

Сконфигурируйте передатчик SAR и получатель. Антенна смотрит в поперечном направлении, ортогональном к направлению рейса.

antenna = phased.CosineAntennaElement('FrequencyRange', [1e9 6e9]);
antennaGain = aperture2gain(aperture,c/fc); 

transmitter = phased.Transmitter('PeakPower', 50e3, 'Gain', antennaGain);
radiator = phased.Radiator('Sensor', antenna,'OperatingFrequency', fc, 'PropagationSpeed', c);

collector = phased.Collector('Sensor', antenna, 'PropagationSpeed', c,'OperatingFrequency', fc);
receiver = phased.ReceiverPreamp('SampleRate', fs, 'NoiseFigure', 30);

Сконфигурируйте канал распространения.

channel = phased.FreeSpace('PropagationSpeed', c, 'OperatingFrequency', fc,'SampleRate', fs,...
    'TwoWayPropagation', true);

Настройка сцены

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

targetpos= [800,0,0;1000,0,0; 1300,0,0]'; 

targetvel = [0,0,0;0,0,0; 0,0,0]';

target = phased.RadarTarget('OperatingFrequency', fc, 'MeanRCS', [1,1,1]);
pointTargets = phased.Platform('InitialPosition', targetpos,'Velocity',targetvel);
% The figure below describes the ground truth based on the target
% locations.
figure(1);h = axes;plot(targetpos(2,1),targetpos(1,1),'*g');hold all;plot(targetpos(2,2),targetpos(1,2),'*r');hold all;plot(targetpos(2,3),targetpos(1,3),'*b');hold off;
set(h,'Ydir','reverse');xlim([-10 10]);ylim([700 1500]);
title('Ground Truth');ylabel('Range');xlabel('Cross-Range');

Симуляция сигнала SAR

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

% Define the broadside angle
refangle = zeros(1,size(targetpos,2));
rxsig = zeros(truncrangesamples,numpulses);
for ii = 1:numpulses
    % Update radar platform and target position
    [radarpos, radarvel] = radarPlatform(slowTime);
    [targetpos,targetvel] = pointTargets(slowTime);
    
    % Get the range and angle to the point targets
    [targetRange, targetAngle] = rangeangle(targetpos, radarpos);
    
    % Generate the LFM pulse
    sig = waveform();
    % Use only the pulse length that will cover the targets.
    sig = sig(1:truncrangesamples);
    
    % Transmit the pulse
    sig = transmitter(sig);
    
    % Define no tilting of beam in azimuth direction
    targetAngle(1,:) = refangle;
    
    % Radiate the pulse towards the targets
    sig = radiator(sig, targetAngle);
    
    % Propagate the pulse to the point targets in free space
    sig = channel(sig, radarpos, targetpos, radarvel, targetvel);
    
    % Reflect the pulse off the targets
    sig = target(sig);
    
    % Collect the reflected pulses at the antenna
    sig = collector(sig, targetAngle);
    
    % Receive the signal  
    rxsig(:,ii) = receiver(sig);
    
end

Визуализируйте полученный сигнал.

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

imagesc(real(rxsig));title('SAR Raw Data')
xlabel('Cross-Range Samples')
ylabel('Range Samples')

Выполните сжатие области значений.

Каждая строка полученного сигнала, который содержит всю информацию от каждого импульса, может быть соответствующей отфильтрованная, чтобы получить de-chirped/range сжатый сигнал.

pulseCompression = phased.RangeResponse('RangeMethod', 'Matched filter', 'PropagationSpeed', c, 'SampleRate', fs);
matchingCoeff = getMatchedFilter(waveform);
[cdata, rnggrid] = pulseCompression(rxsig, matchingCoeff);

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

imagesc(real(cdata));title('SAR Range Compressed Data')
xlabel('Cross-Range Samples')
ylabel('Range Samples')

Сжатие азимута

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

rma_processed = helperRangeMigration(cdata,fastTime,fc,fs,prf,speed,numpulses,c,Rc);
bpa_processed = helperBackProjection(cdata,rnggrid,fastTime,fc,fs,prf,speed,crossRangeResolution,c);

Визуализируйте итоговое изображение SAR.

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

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

figure(1);
imagesc((abs((rma_processed(1700:2300,600:1400).'))));
title('SAR Data focused using Range Migration algorithm ')
xlabel('Cross-Range Samples')
ylabel('Range Samples')

figure(2)
imagesc((abs(bpa_processed(600:1400,1700:2300))));
title('SAR Data focused using Back-Projection algorithm ')
xlabel('Cross-Range Samples')
ylabel('Range Samples')

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

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

Ссылки

Приложение

Алгоритм миграции области значений

function azcompresseddata = helperRangeMigration(sigData,fastTime,fc,fs,prf,speed,numPulses,c,Rc)

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

Установите промежуток частоты области значений.

frequencyRange = linspace(fc-fs/2,fc+fs/2,length(fastTime));
krange = 2*(2*pi*frequencyRange)/c;

Установите перекрестную область значений wavenumber.

kaz = 2*pi*linspace(-prf/2,prf/2,numPulses)./speed;

Сгенерируйте матрицу перекрестной области значений wavenumbers, чтобы совпадать с размером полученного двумерного сигнала SAR

kazimuth = kaz.';
kx = krange.^2-kazimuth.^2;

Установите финал wavenumber достигать фокусировки азимута.

kx = sqrt(kx.*(kx > 0));
kFinal = exp(1i*kx.*Rc);

Выполните двумерный БПФ на области значений сжатый сигнал.

sdata =fftshift(fft(fftshift(fft(sigData,[],1),1),[],2),2);

Выполните объемное сжатие, чтобы получить сжатие азимута в диапазоне ссылки. Выполните фильтрацию 2D сигнала БПФ с новой перекрестной областью значений wavenumber, чтобы достигнуть полной фокусировки в диапазоне ссылки и как побочный продукт, частичная фокусировка целей, не лежащих в диапазоне ссылки.

fsmPol = (sdata.').*kFinal;

Выполните Интерполяцию Stolt, чтобы достигнуть фокусировки для целей, которые не лежат в диапазоне ссылки

stoltPol = fsmPol;
for i = 1:size((fsmPol),1)
    stoltPol(i,:) = interp1(kx(i,:),fsmPol(i,:),krange(1,:));
end
stoltPol(isnan(stoltPol)) = 1e-30;
stoltPol = stoltPol.*exp(-1i*krange.*Rc);
azcompresseddata = ifft2(stoltPol);
end

Алгоритм задней проекции

function data = helperBackProjection(sigdata,rnggrid,fastTime,fc,fs,prf,speed,crossRangeResolution,c)

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

data = zeros(size(sigdata));
azimuthDist = -200:speed/prf:200;%azimuth distance

Ограничьте область значений и пиксели перекрестной области значений, обрабатываемые, чтобы уменьшать время вычислений.

rangelims = [700 1400];
crossrangelims = [-10 10];

Индексируйте сетку области значений в соответствии с пределами области значений.

rangeIdx =  [find(rnggrid>rangelims(1), 1) find(rnggrid<rangelims(2),1,'last')];

Индексируйте расстояние азимута в соответствии с пределами перекрестной области значений.

crossrangeIdxStart = find(azimuthDist>crossrangelims(1),1);
crossrangeIdxStop = find(azimuthDist<crossrangelims(2),1,'last');
for i= rangeIdx(1):rangeIdx(2)  % Iterate over the range indices
    
    % Using desired cross-range resolution, compute the synthetic aperture
    % length
    lsynth= (c/fc)* (c*fastTime(i)/2)/(2*crossRangeResolution);
    lsar = round(lsynth*length(azimuthDist)/azimuthDist(end)) ;
    % Ensure lsar is an odd number
    lsar = lsar + mod(lsar,2);    
    
    % Construct hanning window for cross-range processing, to suppress the
    % azimuthal side lobes
    hn= hanning(lsar).';
    % Iterate over the cross-range indices
    for j= crossrangeIdxStart:crossrangeIdxStop 
        % azimuth distance in x direction over cross-range indices
        posx= azimuthDist(j);
        % range in y-direction over range indices
        posy= c*fastTime(i)/2;
        % initializing count to zero
        count= 0;
        % Iterate over the synthetic aperture
        for k= j-lsar/2 +1:j+ lsar/2 
            % Time delay for each of range and cross-range indices
            td= sqrt((azimuthDist(k)- posx)^2 + posy^2)*2/c;
            cell= round(td*fs) +1 ;
            signal = sigdata(cell,k);
            count= count + hn(k -(j-lsar/2))*signal *exp(1j*2*pi*fc*(td));
        end
        % Processed data at each of range and cross-range indices
        data(i,j)= count;
    end
    
end
end