A radar system's propagation path is modified by the presence of the environment in which it operates. When overlooking a surface, the free space radiation pattern is modified by the reflected wave, which results in an interference pattern with lobes and nulls. In a lobe maximum, a radar return may be increased by as much as 12 dB, whereas in a minimum the return can theoretically go to 0.

The vertical coverage pattern is a detection contour. It offers a look into the performance of a radar system at a constant signal level as specified by the free space range. The vertical coverage takes into account the interference between the direct and ground-reflected rays, as well as refraction. This example will define an L-band radar system in the presence of heavy clutter and will show you how to visualize its 3-dimensional vertical coverage over terrain.

To start, specify an L-band surveillance radar with the following parameters:

• Peak power: 3 kW

• Operating frequency: 1 GHz

• Transmit and receive antenna beamwidth: 20 degrees in azimuth and elevation

• Pulse width: 2 μs

rfs        = 50e3;           % Free-space range (m)
rdrppower  = 3e3;            % Peak power (W)
fc         = 1e9;            % Operating frequency (Hz)
hpbw       = 20;             % Half-power beamwidth (deg)
rdrpulsew  = 2e-6;           % Pulse width (s)
tiltAng    = 1;              % Radar tilt (elevation) angle (deg)
azRotation = 80;             % Radar azimuth rotation angle (deg)

% Vertical coverage parameters
minEl      = 0;              % Minimum vertical coverage angle (deg)
maxEl      = 90;             % Maximum vertical coverage angle (deg)
elStepSize = 1;              % Elevation step size (deg)
azAng      = -60:2:60;       % Azimuth angles for elevation cuts (deg)

Convert the transmit antenna's half-power beamwidth (HPBW) values to gain using the beamwidth2gain function. Assume a cosine rectangular aperture, which is a good approximation for a real-world antenna.

rdrgain    = beamwidth2gain(hpbw,'CosineRectangular'); % Transmit and receive antenna gain (dBi)

Define the Location

The next section defines the radar location. The radar is mounted on a tall tower 100 meters above the ground. The radar altitude is the sum of the ground elevation and the radar tower height referenced to the mean sea level (MSL). Visualize the location using geoplot3. The radar will appear as a red circular point in the upper right-hand corner of the image below.

rdrlat = 39.5;           % Radar latitude (deg)
rdrlon = -105.5;         % Radar longitude (deg)
rdrtowerht = 100;        % Antenna height (m)

% Import the relevant terrain data from the United States Geological
% Survey (USGS)
dtedfile = 'n39_w106_3arc_v2.dt1';
attribution = 'SRTM 3 arc-second resolution. Data available from the U.S. Geological Survey.';

% Visualize the location including terrain using the geographic globe plot
hTerrain = uifigure;
g = geoglobe(hTerrain,'Terrain','southboulder','Basemap','streets');
hold(g,'on')
h_rdrtraj = geoplot3(g,rdrlat,rdrlon,rdralt,'o','Color',[0.6350 0.0780 0.1840],'LineWidth',6,'MarkerSize',6);
campos(g,40.0194565714502,-105.564712896622,13117.1497971643);
campitch(g,-16.023558618352187); Investigate the Free Space Range

For this example, the free space range rfs is the range at which a target would have a specified signal-to-noise ratio (SNR). Assume a large target with a 10 dBsm radar cross section (RCS). At this range, the free space target SNR can be calculated as follows.

trcs       = db2pow(10); % RCS (m^2)
lambda     = freq2wavelen(fc); % Wavelength (m)
'RCS',trcs,'Gain',rdrgain)
tsnr = -0.7449

For our surveillance radar, the desired performance index is a probability of detection (Pd) of 0.8 and probability of false alarm (Pfa) below 1e-3. To make the radar system design more feasible, we can use a pulse integration technique to reduce the required SNR. For this system, we assume noncoherent integration of 64 pulses. A good approximation of the minimum SNR needed for a detection at the specified Pd and Pfa can be computed by the detectability function as follows. Note that the free space SNR satisfies the detectability requirement.

Pd = 0.8;
Pfa = 1e-3;
minsnr = detectability(Pd,Pfa,64);
isdetectable = tsnr >= minsnr
isdetectable = logical
1

The toccgh function allows one to translate receiver detection probabilities into track probabilities. Assuming the default tracker, the required Pd and Pfa translate to the following probability of target track (Pdt) and probability of false track (Pft).

[Pdt,Pft]  = toccgh(Pd,Pfa)
Pdt = 0.9608
Pft = 1.0405e-06

The Pd and Pfa requirements enable a track probability of about 96% and a false track probability on the order of 1e-6.

Plot the Vertical Coverage

The vertical coverage pattern is visualized using a Blake chart, also known as a range-height-angle chart. The range along the x-axis is the propagated range and the height along the y-axis is relative to the origin of the ray.

The vertical coverage contour is calculated using the radarvcd function. The default permittivity model in radarvcd is based on a sea permittivity model in Blake's "Machine Plotting of Radar Vertical-Plane Coverage Diagrams." Such a model is not applicable for the defined scenario, which is over land. Thus, the first step in simulating more realistic propagation is to select a more appropriate permittivity. Use the earthSurfacePermittivity function with the vegetation flag. Assume an ambient temperature of 21.1 degrees Celsius, which is about 70 degrees Fahrenheit. Assume a gravimetric water content of 0.3.

temp = 21.1; % Ambient temperature (degrees Celsius)
gwc = 0.3; % Gravimetric water content
[~,~,epsc] = earthSurfacePermittivity('vegetation',fc,temp,gwc);

Next, specify the antenna. Simulate the antenna as a theoretical sinc antenna pattern and plot.

hBeam = phased.SincAntennaElement('Beamwidth',hpbw);
pattern(hBeam,fc); Get the elevation pattern at 0 degrees azimuth.

elAng = -90:0.01:90;
pat = getVoltagePattern(hBeam,fc,0,elAng);

Specify the atmosphere and calculate the corresponding effective Earth radius. Since the latitude of the radar in this example is 39.5 degrees, a mid-latitude atmosphere model would be most appropriate. Assume that the time period is during the summer.

% Set effective Earth radius using the mid latitude atmosphere model
% in summer
[nidx,N] = refractiveidx([0 1e3], ...
'LatitudeModel','Mid','Season','Summer');

Next, calculate the vertical coverage pattern using the radarvcd function.

'SurfaceRelativePermittivity',epsc, ...
'SurfaceHeightStandardDeviation',30, ...
'Vegetation','Trees', ...
'AntennaPattern',pat,'PatternAngles',elAng, ...
'TiltAngle',tiltAng, ...
'MinElevation',minEl, ...
'ElevationStepSize',elStepSize, ...
'MaxElevation',maxEl);

Use the surface height at the site to obtain the surface refractivity from refractiveidx.

% Calculate an appropriate surface

Plot the vertical coverage using the blakechart function. The blakechart axes are formed using the Central Radio Propagation Laboratory (CRPL) reference atmosphere. The CRPL model is widely used and approximates the refractivity profile as an exponential decay versus height, which has been shown to be an excellent approximation for a normal atmosphere (i.e., an atmosphere that is not undergoing anomalous propagation like ducting). The CRPL model is tailored by setting the surface height refractivity and the refraction exponent to a value that is appropriate for the location in which the system is operating and the time of year.

Update the surface refractivity and the refraction exponent inputs in the function call. In the subsequent plot, the constant-signal-level of the contour is given by the previously calculated target SNR.

rexp = refractionexp(Ns);
blakechart(vcpkm,vcpang,'ScalePower',1,'SurfaceRefractivity',Ns,'RefractionExponent',rexp) Visualize the 3-D Vertical Coverage over Terrain

The next section calculates the vertical coverage at the specified azimuth intervals. The vertical coverage range and angle is converted to height using the range2height function using the 'CRPL' method. The range-height-angle values are then converted to Cartesian.

% Initialize Cartesian X, Y, Z outputs
numAz = numel(azAng);
numRows = numel(minEl:elStepSize:maxEl)+1;
vcpX = zeros(numRows,numAz);
vcpY = zeros(numRows,numAz);
vcpZ = zeros(numRows,numAz);
wgs84 = wgs84Ellipsoid;

% Obtain vertical coverage contour in Cartesian coordinates
for idx = 1:numAz
% Get elevation pattern at this azimuth
pat = getVoltagePattern(hBeam,fc,azAng(idx),elAng);

% Get terrain height information
[xlat,ylon,zlon] = enu2geodetic(xdir,ydir,zdir,rdrlat,rdrlon,rdralt,wgs84);
[~,~,~,htsurf] = los2(Z,R,rdrlat,rdrlon, ...
xlat,ylon,rdralt,zlon,'MSL','MSL');
htstdSurf = std(htsurf(~isnan(htsurf)));

% Calculate vertical coverage pattern
'RangeUnit','m', ...
'HeightUnit','m', ...
'SurfaceRelativePermittivity',epsc, ...
'SurfaceHeightStandardDeviation',htstdSurf, ...
'Vegetation','Trees', ...
'AntennaPattern',pat,'PatternAngles',elAng, ...
'TiltAngle',tiltAng, ...
'MinElevation',1, ...
'ElevationStepSize',1, ...
'MaxElevation',90);

% Calculate associated heights
vcpht = range2height(vcp,rdrtowerht,vcpang, ...
'Method','CRPL','MaxNumIterations',2,'Tolerance',1e-2, ...
'SurfaceRefractivity',Ns,'RefractionExponent',rexp);

% Calculate true slant range and elevation angle to target
[~,vcpgeomrt,vcpelt] = ...
height2range(vcpht,rdrtowerht,vcpang, ...
'Method','CRPL', ...
'SurfaceRefractivity',Ns,'RefractionExponent',rexp);

% Set values for this iteration
vcpelt = [0;vcpelt(:);vcpang(end)];
vcpgeomrt = [0;vcpgeomrt(:);0]; % km

% Convert range, angle, height to Cartesian X, Y, Z
[vcpX(:,idx),vcpY(:,idx),vcpZ(:,idx)] = sph2cart(...
end

Convert the Cartesian vertical coverage values to geodetic and plot on the globe. The three-dimensional vertical coverage will appear as a blue mesh.

[vcpLat,vcpLon,vcpAlt] = enu2geodetic(vcpX,vcpY,vcpZ,rdrlat,rdrlon,rdralt,wgs84);
[xMesh,yMesh,zMesh] = helperCreateMesh(vcpLat,vcpLon,vcpAlt);
geoplot3(g,xMesh,yMesh,zMesh,'LineWidth',1,'Color',[0 0.4470 0.7410]); Summary

This example discussed a method to calculate a 3-dimensional vertical coverage, discussing ways to tailor the analysis to the local atmospheric conditions. The visualization created gives insight into the vertical performance of the radar system, taking into account the interference from the surface reflection and refraction.

References

1. Barton, David K. Radar Equations for Modern Radar. Norwood, MA: Artech House, 2013.

2. Blake, Lamont V. "Radio Ray (Radar) Range-Height-Angle Charts." Naval Research Laboratory, NRL Report 6650, Jan. 22, 1968.

3. Blake, Lamont V. "Ray Height Computation for a Continuous Nonlinear Atmospheric Refractive-Index Profile." Radio Science, Vol. 3 (New Series), No. 1, Jan. 1968, pp. 85–92.

4. Bean, B. R., and G. D. Thayer. CRPL Exponential Reference Atmosphere. Washington, DC: US Government Printing Office, 1959.

% Clean up by closing the geographic globe and removing the imported
% terrain data.
if isvalid(hTerrain)
close(hTerrain)
end
removeCustomTerrain("southboulder")
function pat = getVoltagePattern(hBeam,fc,azAng,elAng)
% Obtain voltage pattern from antenna element
numAz = numel(azAng);
numEl = numel(elAng);
pat = zeros(numAz,numEl);
for ia = 1:numAz
for ie = 1:numEl
pat(ia,ie) = step(hBeam,fc,[azAng(ia);elAng(ie)]);
end
end
end

function [xMesh,yMesh,zMesh] = helperCreateMesh(x,y,z)
% Organizes the inputs into a mesh, which can be plotted using geoplot3
numPts = numel(x);
xMesh = zeros(2*numPts,1);
yMesh = zeros(2*numPts,1);
zMesh = zeros(2*numPts,1);
[nrows,ncols] = size(x);
idxStart = 1;
idxEnd = nrows;
for ii = 1:ncols
if mod(ii,2) == 0
xMesh(idxStart:idxEnd) = x(:,ii);
yMesh(idxStart:idxEnd) = y(:,ii);
zMesh(idxStart:idxEnd) = z(:,ii);
else
xMesh(idxStart:idxEnd) = flipud(x(:,ii));
yMesh(idxStart:idxEnd) = flipud(y(:,ii));
zMesh(idxStart:idxEnd) = flipud(z(:,ii));
end

idxStart = idxEnd + 1;
if ii ~= ncols
idxEnd = idxStart + nrows -1;
else
idxEnd = idxStart + ncols -1;
end
end

for ii = 1:nrows
if mod(ii,2) == 0
xMesh(idxStart:idxEnd) = x(ii,:);
yMesh(idxStart:idxEnd) = y(ii,:);
zMesh(idxStart:idxEnd) = z(ii,:);
else
xMesh(idxStart:idxEnd) = fliplr(x(ii,:));
yMesh(idxStart:idxEnd) = fliplr(y(ii,:));
zMesh(idxStart:idxEnd) = fliplr(z(ii,:));
end
idxStart = idxEnd + 1;
idxEnd = idxStart + ncols -1;
end
end 