Симуляция загромождения моря для морской радиолокационной системы

Открыть Live Script

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

Обзор морских государств

При описании морской загроможденности в первую очередь важно установить физические свойства морской поверхности. При моделировании морского загромождения для радара существует 3 важных параметра:

$σ_{h}$ - стандартное отклонение высоты волны. Высота волны определяется как вертикальное расстояние между гребнем волны и смежным желобом волны.
$β_{0}$ - наклон волны.
$v_{w}$ - скорость ветра.

Из-за нерегулярности волн часто описываются физические свойства моря в терминах морских состояний. Число состояния моря Дугласа является широко используемой шкалой, которая представляет широкую область значений физических морских свойств, таких как высоты волн и соответствующие скорости ветра. В нижнем конце шкалы состояние моря 0 представляет собой спокойное, стекловидное состояние моря. Шкала затем протекает через состояния, которые представляют все от слегка испорченного в 1 до грубых морей с высокими высотами волн в морском состоянии 5. Высоты волн в морском состоянии 8 могут быть больше 9 м и более.

Использование searoughness функция, постройте график свойств моря для состояний моря с 1 по 5. Обратите внимание на медленное увеличение наклона волны $β_{0}$ с состоянием моря. Это является результатом увеличения длины волны и высоты волны со скоростью ветра, хотя и с различными факторами.

% Analyze for sea states 1 through 5
ss = 1:5; % Sea states 

% Initialize outputs
numSeaStates = numel(ss);
hgtsd = zeros(1,numSeaStates);
beta0 = zeros(1,numSeaStates);
vw= zeros(1,numSeaStates);

% Obtain sea state properties
for is = 1:numSeaStates
    [hgtsd(is),beta0(is),vw(is)] = searoughness(ss(is));
end

% Plot results
helperPlotSeaRoughness(ss,hgtsd,beta0,vw);

Figure contains 3 axes. Axes 1 with title Sea Wave Roughness contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line.

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

Отражательная способность

Морская поверхность сложена водой со средней минерализацией около 35 частей на тысячу. Коэффициент отражения морской воды близок к -1 для микроволновых частот и при низких углах выпаса.

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

Исследуйте отражательную способность поверхности моря в зависимости от частоты для различных морских состояний, используя seareflectivity функция. Установите угол выпаса равным 0,5 степеням и рассмотрим частоты в области значений 500 МГц до 35 ГГц.

grazAng = 0.5; % Grazing angle (deg)
freq = linspace(0.5e9,35e9,100); % Frequency (Hz)
pol = 'H'; % Horizontal polarization

% Initialize reflectivity output
numFreq = numel(freq);
nrcsH = zeros(numFreq,numSeaStates);

% Calculate reflectivity
for is = 1:numSeaStates
    nrcsH(:,is) = seareflectivity(ss(is),grazAng,freq,'Polarization',pol);
end

% Plot reflectivity
helperPlotSeaReflectivity(ss,grazAng,freq,nrcsH,pol);

$Figure contains an axes. The axes with title Sea State Reflectivity \sigma_0 contains 5 objects of type line. These objects represent SS 1, H, SS 2, H, SS 3, H, SS 4, H, SS 5, H.$

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

Поляризационные эффекты

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

pol = 'V'; % Vertical polarization

% Initialize reflectivity output
numFreq = numel(freq);
nrcsV = zeros(numFreq,numSeaStates);

% Calculate reflectivity
for is = 1:numSeaStates
    nrcsV(:,is) = seareflectivity(ss(is),grazAng,freq,'Polarization',pol);
end

% Plot reflectivity
hAxes = helperPlotSeaReflectivity(ss,grazAng,freq,nrcsH,'H');
helperPlotSeaReflectivity(ss,grazAng,freq,nrcsV,'V',hAxes);

$Figure contains an axes. The axes with title Sea State Reflectivity \sigma_0 contains 10 objects of type line. These objects represent SS 1, H, SS 2, H, SS 3, H, SS 4, H, SS 5, H, SS 1, V, SS 2, V, SS 3, V, SS 4, V, SS 5, V.$

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

Эффекты угла выпаса

Рассмотрим эффект угла выпаса. Вычислите отражательную способность моря в области значений от 0,1 до 60 степеней при частоте L-диапазона 1,5 ГГц.

grazAng = linspace(0.1,60,100); % Grazing angle (deg)
freq = 1.5e9; % L band frequency (Hz)

% Initialize reflectivity output
numGrazAng = numel(grazAng);
nrcsH = zeros(numGrazAng,numSeaStates);
nrcsV = zeros(numGrazAng,numSeaStates);

% Calculate reflectivity
for is = 1:numSeaStates
    nrcsH(:,is) = seareflectivity(ss(is),grazAng,freq,'Polarization','H');
    nrcsV(:,is) = seareflectivity(ss(is),grazAng,freq,'Polarization','V');
end

% Plot reflectivity
hAxes = helperPlotSeaReflectivity(ss,grazAng,freq,nrcsH,'H');
helperPlotSeaReflectivity(ss,grazAng,freq,nrcsV,'V',hAxes);
ylim(hAxes,[-60 -10]);

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

Пример радара морского наблюдения

Вычисление отношения загромождения к шуму

Рассмотрим горизонтально поляризованную радиолокационную систему морского наблюдения, работающую на 6 ГГц (C полоса). Определите радиолокационную систему.

% Radar parameters 
freq = 6e9;         % C-band frequency (Hz) 
anht = 20;          % Height (m) 
ppow = 200e3;       % Peak power (W)
tau  = 200e-6;      % Pulse width (sec)
prf  = 300;         % PRF (Hz)
azbw = 10;          % Half-power azimuth beamwidth (deg)
elbw = 30;          % Half-power elevation beamwidth (deg)
Gt   = 22;          % Transmit gain (dB) 
Gr   = 10;          % Receive gain (dB)
nf   = 3;           % Noise figure (dB)
Ts   = systemp(nf); % System temperature (K)

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

% Sea parameters
ss = 2;             % Sea state    

% Calculate surface state
[hgtsd,beta0] = searoughness(ss);

% Setup geometry
anht = anht + 2*hgtsd;         % Average height above clutter (m) 
surfht = 3*hgtsd;              % Surface height (m) 

% Calculate maximum range for simulation
Rua = time2range(1/prf); % Maximum unambiguous range (m)
Rhoriz = horizonrange(anht,'SurfaceHeight',surfht); % Horizon range (m)
Rmax = min(Rua,Rhoriz); % Maximum simulation range (m)

% Generate vector of ranges for simulation
Rm = linspace(100,Rmax,1000); % Range (m)
Rkm = Rm*1e-3;                % Range (km) 

% Calculate sea clutter reflectivity
grazAng = grazingang(anht,Rm,'TargetHeight',surfht); 
nrcs = seareflectivity(ss,grazAng,freq);
helperPlotSeaReflectivity(ss,grazAng,freq,nrcs,'H');

$Figure contains an axes. The axes with title Sea State Reflectivity \sigma_0 contains an object of type line. This object represents SS 2, H.$

Затем вычислите радиолокационное сечение (RCS) загромождения с помощью clutterSurfaceRCS функция. Обратите внимание на падение загромождения RCS при достижении области значений горизонта.

% Calculate clutter RCS 
rcs = clutterSurfaceRCS(nrcs,Rm,azbw,elbw,grazAng(:),tau); 
rcsdB = pow2db(rcs); % Convert to decibels for plotting
hAxes = helperPlot(Rkm,rcsdB,'RCS','Clutter RCS (dBsm)','Clutter Radar Cross Section (RCS)');
helperAddHorizLine(hAxes,Rhoriz);

Figure contains an axes. The axes with title Clutter Radar Cross Section (RCS) contains 2 objects of type line, constantline. These objects represent RCS, Horizon Range.

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

% Convert frequency to wavelength
lambda = freq2wavelen(freq); 

% Calculate and plot the clutter-to-noise ratio
cnr = radareqsnr(lambda,Rm(:),ppow,tau,... 
        'gain',[Gt Gr],'rcs',rcs,'Ts',Ts); % dB
hAxes = helperPlot(Rkm,cnr,'CNR','CNR (dB)','Clutter-to-Noise Ratio (CNR)');
ylim(hAxes,[-80 100]);
helperAddHorizLine(hAxes,Rhoriz);
helperAddBelowClutterPatch(hAxes);

Figure contains an axes. The axes with title Clutter-to-Noise Ratio (CNR) contains 3 objects of type patch, line, constantline. These objects represent Clutter Below Noise, CNR, Horizon Range.

% Range when clutter falls below noise
helperFindClutterBelowNoise(Rkm,cnr);

Range at which clutter falls below noise (km) = 18.04

Принимая во внимание путь распространения

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

% Calculate radar propagation factor for clutter 
Fc = radarpropfactor(Rm,freq,anht,surfht, ... 
    'SurfaceHeightStandardDeviation',hgtsd,...
    'SurfaceSlope',beta0,...
    'ElevationBeamwidth',elbw);
helperPlot(Rkm,Fc,'Propagation Factor', ...
    'Propagation Factor (dB)', ...
    'One-Way Clutter Propagation Factor F_C');

Figure contains an axes. The axes with title One-Way Clutter Propagation Factor F_C contains an object of type line. This object represents Propagation Factor.

В рамках вышеописанного графика видны 2 области распространения:

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

Обычно коэффициент распространения загромождения и отражательная способность моря объединяются как продукт $σ_{C} F_{C}^{4}$ , потому что измерения отражательной способности поверхности, как правило, являются измерениями продукта, а не просто отражательной способности $σ_{C}$ . Вычислите этот продукт и постройте график результатов.

% Combine clutter reflectivity and clutter propagation factor
FcLinear = db2mag(Fc); % Convert to linear units
combinedFactor = nrcs.*FcLinear.^2;
combinedFactordB = pow2db(combinedFactor);
helperPlot(Rkm,combinedFactordB,'\sigma_CF_C', ...
    '\sigma_CF_C (dB)', ...
    'One-Way Sea Clutter Propagation Factor and Reflectivity');

$Figure contains an axes. The axes with title One-Way Sea Clutter Propagation Factor and Reflectivity contains an object of type line. This object represents \sigma_CF_C.$

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

% Calculate one-way loss associated with rain
elAng = height2el(surfht,anht,Rm); % Elevation angle (deg) 
numEl = numel(elAng);
Latmos = zeros(numEl,1);
for ie = 1:numEl
    Latmos(ie,:) = tropopl(Rm(ie),freq,anht,elAng(ie));
end
helperPlot(Rkm,Latmos,'Atmospheric Loss','Loss (dB)','One-Way Atmospheric Loss');

Figure contains an axes. The axes with title One-Way Atmospheric Loss contains an object of type line. This object represents Atmospheric Loss.

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

% Re-calculate CNR including radar propagation factor 
cnr = radareqsnr(lambda,Rm(:),ppow,tau,... 
        'gain',[Gt Gr],'rcs',rcs,'Ts',Ts, ...
        'PropagationFactor',Fc,...
        'AtmosphericLoss',Latmos); % dB
helperAddPlot(Rkm,cnr,'CNR + Propagation Factor + Atmospheric Loss',hAxes);

Figure contains an axes. The axes with title Clutter-to-Noise Ratio (CNR) contains 4 objects of type patch, line, constantline. These objects represent Clutter Below Noise, CNR, Horizon Range, CNR + Propagation Factor + Atmospheric Loss.

% Range when clutter falls below noise
helperFindClutterBelowNoise(Rkm,cnr);

Range at which clutter falls below noise (km) = 10.44

Понимание погодных эффектов

Так же, как погода влияет на обнаружение цели, погода также влияет на обнаружение загромождения. Рассмотрим эффект дождя в моделируемых областях значений. Сначала вычислите ослабление дождя.

% Calculate one-way loss associated with rain
rr = 50;                           % Rain rate (mm/h)
polAng = 0;                        % Polarization tilt angle (0 degrees for horizontal) 
elAng = height2el(surfht,anht,Rm); % Elevation angle (deg) 
numEl = numel(elAng);
Lrain = zeros(numEl,1);
for ie = 1:numEl
    Lrain(ie,:) = cranerainpl(Rm(ie),freq,rr,elAng(ie),polAng);
end
helperPlot(Rkm,Lrain,'Rain Loss','Loss (dB)','One-Way Rain Loss');

Figure contains an axes. The axes with title One-Way Rain Loss contains an object of type line. This object represents Rain Loss.

Пересчет CNR. Включите путь распространения и потери дождя. Обратите внимание, что из-за наличия дождя происходит лишь незначительное снижение CNR.

% Re-calculate CNR including radar propagation factor and rain loss
cnr = radareqsnr(lambda,Rm(:),ppow,tau,... 
        'gain',[Gt Gr],'rcs',rcs,'Ts',Ts, ...
        'PropagationFactor',Fc, ...
        'AtmosphericLoss',Latmos + Lrain); % dB
helperAddPlot(Rkm,cnr,'CNR + Propagation Factor + Atmospheric Loss + Rain',hAxes);

Figure contains an axes. The axes with title Clutter-to-Noise Ratio (CNR) contains 5 objects of type patch, line, constantline. These objects represent Clutter Below Noise, CNR, Horizon Range, CNR + Propagation Factor + Atmospheric Loss, CNR + Propagation Factor + Atmospheric Loss + Rain.

% Range when clutter falls below noise
helperFindClutterBelowNoise(Rkm,cnr);

Range at which clutter falls below noise (km) = 9.61

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

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

Сильная зависимость от морского состояния,
Пропорциональная зависимость от частоты,
Зависимость от поляризации, которая уменьшается с увеличением частоты, и
Сильная зависимость от угла выпаса при низких углах выпаса

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

Ссылки

Бартон, Дэвид К. Основные уравнения радиолокации для современного радара. 1-е издание. Norwood, MA: Artech House, 2013.
Блейк, Л.В. Построение графических изображений радиолокационных вертикально-плоскостных покрытий. Отчет военно-морской исследовательской лаборатории 7098, 1970.
Gregers-Hansen, V. and Mittal, R. «Улучшенная эмпирическая модель для радиолокационной отражательной способности морского загромождения». NRL/MR/5310-12-9346, 27 апреля 2012 года.
Ричардс, М. А., Джим Шеер и Уильям А. Холм. Принципы современного радара. Raleigh, NC: SciTech Pub., 2010.

function helperPlotSeaRoughness(ss,hgtsd,beta0,vw)
% Creates 3x1 plot of sea roughness outputs 

% Create figure
figure

% Plot standard deviation of sea wave height
subplot(3,1,1)
plot(ss,hgtsd,'-o','LineWidth',1.5)
ylabel([sprintf('Wave\nHeight') '\sigma (m)'])
title('Sea Wave Roughness')
grid on; 

% Plot sea wave slope
subplot(3,1,2)
plot(ss,beta0,'-o','LineWidth',1.5)
ylabel(sprintf('Wave\nSlope (deg)'))
grid on; 

% Plot wind velocity 
subplot(3,1,3)
plot(ss,vw,'-o','LineWidth',1.5)
xlabel('Sea State')
ylabel(sprintf('Wind\nVelocity (m/s)'))
grid on; 
end

function hAxes = helperPlotSeaReflectivity(ss,grazAng,freq,nrcs,pol,hAxes)
% Plot sea reflectivities 

% Create figure and new axes if axes are not passed in
newFigure = false; 
if nargin < 6
    figure();
    hAxes = gca; 
    newFigure = true;
end

% Get polarization string
switch lower(pol)
    case 'h'
        lineStyle = '-';
    otherwise
        lineStyle = '--';
end

% Plot
if numel(grazAng) == 1
    hLine = semilogx(hAxes,freq(:).*1e-9,pow2db(nrcs),lineStyle,'LineWidth',1.5);
    xlabel('Frequency (GHz)')
else
    hLine = plot(hAxes,grazAng(:),pow2db(nrcs),lineStyle,'LineWidth',1.5);
    xlabel('Grazing Angle (deg)')
end

% Set display names
numLines = size(nrcs,2);
for ii = 1:numLines
    hLine(ii).DisplayName = sprintf('SS %d, %s',ss(ii),pol);
    if newFigure
        hLine(ii).Color = brighten(hLine(ii).Color,0.5);
    end
end

% Update labels and axes 
ylabel('Reflectivity \sigma_0 (dB)')
title('Sea State Reflectivity \sigma_0')
grid on
axis tight
hold on; 

% Add legend
legend('Location','southoutside','NumColumns',5,'Orientation','Horizontal');
end

function varargout = helperPlot(Rkm,y,displayName,ylabelStr,titleName)
% Used in CNR analysis

% Create figure 
hFig = figure;
hAxes = axes(hFig); 

% Plot
plot(hAxes,Rkm,y,'LineWidth',1.5,'DisplayName',displayName);
grid(hAxes,'on');
hold(hAxes,'on');
xlabel(hAxes,'Range (km)')
ylabel(hAxes,ylabelStr);
title(hAxes,titleName);
axis(hAxes,'tight');

% Add legend
legend(hAxes,'Location','Best')

% Output axes
if nargout ~= 0 
    varargout{1} = hAxes;
end
end

function helperAddPlot(Rkm,y,displayName,hAxes)
% Used in CNR analysis

% Plot
ylimsIn = get(hAxes,'Ylim');
plot(hAxes,Rkm,y,'LineWidth',1.5,'DisplayName',displayName);
axis(hAxes,'tight');
ylimsNew = get(hAxes,'Ylim');
set(hAxes,'Ylim',[ylimsIn(1) ylimsNew(2)]); 
end

function helperAddHorizLine(hAxes,Rhoriz)
% Add vertical line indication horizon range

xline(Rhoriz.*1e-3,'--','DisplayName','Horizon Range','LineWidth',1.5);
xlims = get(hAxes,'XLim');
xlim([xlims(1) Rhoriz.*1e-3*(1.05)]);
end

function helperAddBelowClutterPatch(hAxes)
% Add patch indicating when clutter falls below the noise

xlims = get(hAxes,'Xlim');
ylims = get(hAxes,'Ylim');
x = [xlims(1) xlims(1) xlims(2) xlims(2) xlims(1)];
y = [ylims(1) 0 0 ylims(1) ylims(1)];
hP = patch(hAxes,x,y,[0.8 0.8 0.8], ...
    'FaceAlpha',0.3,'EdgeColor','none','DisplayName','Clutter Below Noise');
uistack(hP,'bottom');
end

function helperFindClutterBelowNoise(Rkm,cnr)
% Find the point at which the clutter falls below the noise
idxNotNegInf = ~isinf(cnr); 
Rclutterbelow = interp1(cnr(idxNotNegInf),Rkm(idxNotNegInf),0); 
fprintf('Range at which clutter falls below noise (km) = %.2f\n',Rclutterbelow)
end

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