Электронное сканирование с использованием равномерного прямоугольного массива

Открыть Live Script

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

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

Сначала создадим радар с фазированной решеткой. Мы повторно используем большинство подсистем, построенных в примере Симуляция Тестовых Сигналов для Радарного Приемника. Читателям рекомендуется изучить детали разработки радарных систем с помощью этого примера. Основным различием является то, что вместо исходной одной антенны мы используем равномерный прямоугольный массив (URA) 30 на 30.

Существующий проект радара соответствует следующим спецификациям.

pd = 0.9;            % Probability of detection
pfa = 1e-6;          % Probability of false alarm
max_range = 5000;    % Maximum unambiguous range
tgt_rcs = 1;         % Required target radar cross section
int_pulsenum = 10;   % Number of pulses to integrate

Загрузите радиолокационную систему и найдите системные параметры.

load BasicMonostaticRadarExampleData;

fc = radiator.OperatingFrequency;         % Operating frequency (Hz)
v = radiator.PropagationSpeed;            % Wave propagation speed (m/s)
lambda = v/fc;                            % Wavelength (m)
fs = waveform.SampleRate;                 % Sampling frequency (Hz)
prf = waveform.PRF;                       % Pulse repetition frequency (Hz)

Далее зададим равномерный прямоугольный массив 30 на 30.

ura = phased.URA('Element',antenna,...
    'Size',[30 30],'ElementSpacing',[lambda/2, lambda/2]);
% Configure the antenna elements such that they only transmit forward
ura.Element.BackBaffled = true;

% Visualize the response pattern.
pattern(ura,fc,'PropagationSpeed',physconst('LightSpeed'),...
    'Type','powerdb');

Связать массив с излучателем и коллектором.

radiator.Sensor = ura;
collector.Sensor = ura;

% We need to set the WeightsInputPort property to true to enable it to
% accept transmit beamforming weights
radiator.WeightsInputPort = true;

Теперь нам нужно пересчитать степень передачи. Исходная степень передачи была вычислена на основе одной антенны. Для массива с 900 элементами степень, необходимая для каждого элемента, намного меньше.

% Calculate the array gain
arraygain = phased.ArrayGain('SensorArray',ura,'PropagationSpeed',v);
ag = arraygain(fc,[0;0]);

% Calculate the peak power
snr_min = albersheim(pd, pfa, int_pulsenum);
peak_power = ((4*pi)^3*noisepow(1/waveform.PulseWidth)*max_range^4*...
    db2pow(snr_min))/(db2pow(2*(transmitter.Gain+ag))*tgt_rcs*lambda^2)

peak_power = 0.0065

Новая пиковая степень составляет 0,0065 Вт.

% Set the peak power of the transmitter
transmitter.PeakPower = peak_power;

Нам также нужно разработать график сканирования фазированной решетки. Чтобы упростить пример, мы ищем только в азимутальной размерности. Мы требуем, чтобы радар искал от 45 степеней до -45 степеней по азимуту. Время повторного рассмотрения должно быть менее 1 секунды, что означает, что радар должен пересмотреть один и тот же угол азимута в течение 1 секунды.

initialAz = 45; endAz = -45; 
volumnAz = initialAz - endAz;

Чтобы определить необходимое количество сканов, мы должны знать ширину луча отклика массива. Мы используем эмпирическую формулу, чтобы оценить 3-dB ширину луча.

$G = \frac{4 π}{θ^{2}}$

где $G$ - коэффициент усиления массива, $θ$ - 3-dB лучевая ширина.

% Calculate 3-dB beamwidth
theta = radtodeg(sqrt(4*pi/db2pow(ag)))

theta = 6.7703

Ширина 3-dB луча составляет 6,77 степени. Чтобы обеспечить некоторое перекрытие луча в пространстве, мы выбираем шаг скана равным 6 степеням.

scanstep = -6;
scangrid = initialAz+scanstep/2:scanstep:endAz;
numscans = length(scangrid);    
pulsenum = int_pulsenum*numscans;

% Calculate revisit time
revisitTime = pulsenum/prf

revisitTime = 0.0050

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

Определение цели

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

tgtpos = [[3532.63; 800; 0],[2020.66; 0; 0]];
tgtvel = [[-100; 50; 0],[60; 80; 0]];
tgtmotion = phased.Platform('InitialPosition',tgtpos,'Velocity',tgtvel);

tgtrcs = [1.6 2.2];
target = phased.RadarTarget('MeanRCS',tgtrcs,'OperatingFrequency',fc);

% Calculate the range, angle, and speed of the targets
[tgtrng,tgtang] = rangeangle(tgtmotion.InitialPosition,...
    sensormotion.InitialPosition);

numtargets = length(target.MeanRCS);

Импульсный синтез

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

% Create the steering vector for transmit beamforming
steeringvec = phased.SteeringVector('SensorArray',ura,...
    'PropagationSpeed',v);

% Create the receiving beamformer
beamformer = phased.PhaseShiftBeamformer('SensorArray',ura,...
    'OperatingFrequency',fc,'PropagationSpeed',v,...
    'DirectionSource','Input port');

% Define propagation channel for each target
channel = phased.FreeSpace(...
    'SampleRate',fs,...
    'TwoWayPropagation',true,...
    'OperatingFrequency',fc);

fast_time_grid = unigrid(0, 1/fs, 1/prf, '[)');

% Pre-allocate array for improved processing speed
rxpulses = zeros(numel(fast_time_grid),pulsenum);

for m = 1:pulsenum

    % Update sensor and target positions
    [sensorpos,sensorvel] = sensormotion(1/prf);
    [tgtpos,tgtvel] = tgtmotion(1/prf);

    % Calculate the target angles as seen by the sensor
    [tgtrng,tgtang] = rangeangle(tgtpos,sensorpos);

    % Calculate steering vector for current scan angle
    scanid = floor((m-1)/int_pulsenum) + 1;
    sv = steeringvec(fc,scangrid(scanid));
    w = conj(sv);
    
    % Form transmit beam for this scan angle and simulate propagation
    pulse = waveform();
    [txsig,txstatus] = transmitter(pulse);
    txsig = radiator(txsig,tgtang,w);
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel);
    
    % Reflect pulse off of targets
    tgtsig = target(txsig);
    
    % Beamform the target returns received at the sensor
    rxsig = collector(tgtsig,tgtang);
    rxsig = receiver(rxsig,~(txstatus>0));    
    rxpulses(:,m) = beamformer(rxsig,[scangrid(scanid);0]);
end

Согласованный фильтр

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

% Matched filtering
matchingcoeff = getMatchedFilter(waveform);
matchedfilter = phased.MatchedFilter(...
    'Coefficients',matchingcoeff,...
    'GainOutputPort',true);
[mf_pulses, mfgain] = matchedfilter(rxpulses);
mf_pulses = reshape(mf_pulses,[],int_pulsenum,numscans);

matchingdelay = size(matchingcoeff,1)-1;
sz_mfpulses = size(mf_pulses);
mf_pulses = [mf_pulses(matchingdelay+1:end) zeros(1,matchingdelay)];
mf_pulses = reshape(mf_pulses,sz_mfpulses);

% Pulse integration
int_pulses = pulsint(mf_pulses,'noncoherent');
int_pulses = squeeze(int_pulses);

% Visualize
r = v*fast_time_grid/2;
X = r'*cosd(scangrid); Y = r'*sind(scangrid);
clf;
pcolor(X,Y,pow2db(abs(int_pulses).^2));
axis equal tight 
shading interp
axis off
text(-800,0,'Array');
text((max(r)+10)*cosd(initialAz),(max(r)+10)*sind(initialAz),...
    [num2str(initialAz) '^o']);
text((max(r)+10)*cosd(endAz),(max(r)+10)*sind(endAz),...
    [num2str(endAz) '^o']);
text((max(r)+10)*cosd(0),(max(r)+10)*sind(0),[num2str(0) '^o']);
colorbar;

На карте скана мы четко видим два пика. Близкий находится около 0 степеней азимута, удаленный - около 10 степени азимута.

Обнаружение и оценка области значений

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

range_gates = v*fast_time_grid/2;
tvg = phased.TimeVaryingGain(...
    'RangeLoss',2*fspl(range_gates,lambda),...
    'ReferenceLoss',2*fspl(max(range_gates),lambda));
tvg_pulses = tvg(mf_pulses);

% Pulse integration
int_pulses = pulsint(tvg_pulses,'noncoherent');
int_pulses = squeeze(int_pulses);

% Calculate the detection threshold

% sample rate is twice the noise bandwidth in the design system
noise_bw = receiver.SampleRate/2;
npower = noisepow(noise_bw,...
    receiver.NoiseFigure,receiver.ReferenceTemperature);
threshold = npower * db2pow(npwgnthresh(pfa,int_pulsenum,'noncoherent'));
% Increase the threshold by the matched filter processing gain
threshold = threshold * db2pow(mfgain);

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

N = 51;
clf;
surf(X(N:end,:),Y(N:end,:),...
    pow2db(abs(int_pulses(N:end,:)).^2));
hold on;
mesh(X(N:end,:),Y(N:end,:),...
    pow2db(threshold*ones(size(X(N:end,:)))),'FaceAlpha',0.8);
view(0,56);
axis off

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

[~,peakInd] = findpeaks(int_pulses(:),'MinPeakHeight',sqrt(threshold));
[rngInd,angInd] = ind2sub(size(int_pulses),peakInd);
est_range = range_gates(rngInd); % Estimated range
est_angle = scangrid(angInd);    % Estimated direction

Доплеровская оценка

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

for m = numtargets:-1:1
    [p, f] = periodogram(mf_pulses(rngInd(m),:,angInd(m)),[],256,prf, ...
                'power','centered');
    speed_vec = dop2speed(f,lambda)/2;   

    spectrum_data = p/max(p);
    [~,dop_detect1] = findpeaks(pow2db(spectrum_data),'MinPeakHeight',-5);
    sp(m) = speed_vec(dop_detect1); % Estimated Doppler speed
end

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

------------------------------------------------------------------------
               Estimated (true) target parameters
------------------------------------------------------------------------
                 Range (m)       Azimuth (deg)   Speed (m/s)
Target 1:    3625.00 (3622.08)   12.00 (12.76)   86.01 (86.49)
Target 2:    2025.00 (2020.66)    0.00 (0.00)   -59.68 (-60.00)

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

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

Документация