Main Content

Simulate a Maritime Radar PPI

This example shows how to simulate a plan position indicator (PPI) radar image for a rotating antenna array in a maritime environment. You will configure a radar scenario and spectral sea surface model, emulate a cargo ship with a cuboid extended target, generate return signals, and plot a PPI image of the result.

Configure the Radar Scenario

Set the RNG seed for repeatable results.

rng default

The scenario will consist of a large rotating uniform linear array (ULA) mounted in a fixed position above a sea surface. A large cuboid target will be placed on the sea surface.

Define the radar system parameters. Use an X-band radar at 10 GHz with a 1 degree azimuth beamwidth and 14 meter range resolution.

freq = 10e9; % Hz
azbw = 1;    % deg
rngres = 14; % m

Calculate the number of uniformly-spaced elements required to achieve the desired beamwidth, and create a phased.ULA. Set the BackBaffled property of the array element to true so that only returns from the pointing direction are included.

sinc3db = 0.8859;
numelems = round(2*sinc3db/deg2rad(azbw));
lambda = freq2wavelen(freq);
array = phased.ULA(numelems,lambda/2);
array.Element.BackBaffled = true;

Let the radar rotate at 50 rpm and express this in deg/sec. Calculate the rotation period from this rate.

rotrate = 50/60*360;     % deg/sec
rotperiod = 360/rotrate; % sec

Calculate the whole number of pulses required to minimally cover 360 degrees with the specified azimuth beamwidth. Use the rotation period to find the required pulse repetition frequency (PRF).

npulses = ceil(360/azbw);
prf = npulses/rotperiod; % Hz

The fast-time sampling rate must be a whole multiple of the PRF. Use the desired range resolution to find the required sampling rate then adjust to match this constraint with the PRF. The radar will use a short pulse with one sample per resolution cell.

sampleRate = 1/range2time(rngres); % Hz
sampleRate = prf*round(sampleRate/prf);
rngres = time2range(1/sampleRate); % m

Now create the radarTransceiver and set the required parameters.

rdr = radarTransceiver;
rdr.TransmitAntenna.OperatingFrequency = freq;
rdr.ReceiveAntenna.OperatingFrequency = freq;
rdr.TransmitAntenna.Sensor = array;
rdr.ReceiveAntenna.Sensor = array;
rdr.Waveform.PRF = prf;
rdr.Receiver.SampleRate = sampleRate;
rdr.Waveform.SampleRate = sampleRate;
rdr.Waveform.PulseWidth = 1/sampleRate;

Create a radarScenario and set the UpdateRate to 0 to let the update rate be derived from the radarTransceiver configuration.

scenario = radarScenario('UpdateRate',0);

Now configure the sea surface. The surface will be defined on a square region of size 2 km-by-2 km.

seaLength = 2e3; % m

The seaSpectrum object is used to define the spectral model for the surface. The Resolution property defines the spacing of the underlying grid in the spatial domain, and the WaveVectorSpacing defines the spacing of the grid in the frequency domain. The resolution must be a factor of the length of the surface, so specify a desired surface resolution of 1/4th of the radar's range resolution (which is sufficient to capture the shape of the waves), then adjust to match this constraint. Set the WaveVectorSpacing to 2πdivided by the length of the surface.

seaRes = rngres/4; % m
seaRes = seaLength/round(seaLength/seaRes);
spec = seaSpectrum('Resolution',seaRes,'WaveVectorSpacing',2*pi/seaLength);

The surfaceReflectivitySea object is used to specify the reflectivity model that will be associated with the surface. Use the NRL reflectivity model with a sea state of 3 for horizontal polarization.

refl = surfaceReflectivitySea('Model','NRL','SeaState',3,'Polarization','H');

Create the 2 km-by-2 km sea surface using the reflectivity and spectral models. Specify the wind speed as 10 m/s, and the wind direction as 0 degrees so the wind is blowing in the +X direction.

windSpeed = 10; % m/s
windDir = 0;    % deg
bdry = [-1 1;-1 1]*seaLength/2;

Create the radar platform. Use the kinematicTrajectory object to mount the radar 24 meters above sea level, and use the AngularVelocity parameter to specify the rotation rate about the Z axis in radians per second.

rdrHeight = 24; % m
traj = kinematicTrajectory('Position',[0 0 rdrHeight],'AngularVelocitySource','Property','AngularVelocity',[0 0 deg2rad(rotrate)]);

Now define the cuboid target. Specify the target dimensions, total RCS, position, heading, and speed.

tgtdims = struct('Length',120,'Width',18,'Height',22,'OriginOffset',[0 0 0]); % m
tgtrcs = 40;  % dBsm
tgtpos = [seaLength/3 seaLength/16 5];
tgthdg = -30; % deg
tgtspd = 8;   % m/s

Form the target velocity vector, create a kinematic trajectory, and add the target platform to the scene.

tgtvel = tgtspd*[cosd(tgthdg) sind(tgthdg) 0];
tgttraj = kinematicTrajectory('Position',tgtpos,'Velocity',tgtvel,'Orientation',rotz(tgthdg).');

Finally, enable clutter generation for the radar by calling the clutterGenerator method on the scenario. The 3 dB beamwidth of the beam is used for clutter generation by default. In order to capture more of the surface return, disable this by setting UseBeam to false, and instead use a RingClutterRegion with an azimuth coverage 16 times that of the beam to encompass some sidelobes. Set the clutter resolution to half the range resolution. Returns will only be considered from a minimum ground range of 200 meters out to the maximum ground range of the surface. The azimuth center can be 0 to indicate that the region is centered about the +X direction. The ringClutterRegion method returns a handle to this region so that the parameters can be modified during the simulation.

clut = clutterGenerator(scenario,rdr,'UseBeam',false,'Resolution',rngres/2);
minGndRng = 200;
maxGndRng = seaLength/2;
azCov = 16*azbw;
azCen = 0;
reg = ringClutterRegion(clut,minGndRng,maxGndRng,azCov,azCen);

Use the provided helper function to create a visualization of the scenario.


Figure contains an axes object. The axes object with title Scenario Geometry, xlabel X (km), ylabel Y (km) contains 9 objects of type surface, line, text. One or more of the lines displays its values using only markers

Run the Simulation and Collect Returns

Each frame of the simulation will generate a range profile for one azimuth pointing direction. After receiving the raw IQ data on each frame, a phased.RangeResponse object will be used to perform matched filtering. Create the response object now, specifying the sample rate used, and get the required matched filter coefficients from the radar's waveform object.

resp = phased.RangeResponse('RangeMethod','Matched filter','SampleRate',sampleRate);
mfc = getMatchedFilter(rdr.Waveform);

Specify how much of the full 360 degrees you would like to simulate. The scan starts at 0 degrees azimuth, and 45 degrees of coverage is enough to see the target. The full 360 degrees takes about 10 minutes on a machine with 64 GB of RAM and a 3.6 GHz CPU.

azcov = 45;

Set the scenario stop time based on the desired azimuth coverage.

scenario.StopTime = azcov/360*npulses/prf;

Run the simulation. Keep track of the frame number of each loop, and get the vector of range bins from the range response object. The matrix ppi will contain the signal data formatted as range-by-azimuth.

frame = 0;
nrng = floor(time2range(1/prf)/rngres);
ppi = zeros(nrng,npulses);
while advance(scenario)
    frame = frame + 1;
    reg.AzimuthCenter = (frame-1)*rotrate/prf;
    iqsig = receive(scenario);
    [ppi(:,frame),rngbins] = resp(sum(iqsig{1},2),mfc);

Create a PPI Image

A PPI image consists of a set of range profiles arranged in radial lines to form a circular image of the scene in Cartesian space.

The IQ data covers the full range ambiguity, so start by trimming the data to the ranges of interest. Use the min and max ground range defined earlier, then find the range gate indices that correspond to those ranges.

minRange = sqrt(minGndRng^2 + rdrHeight^2);
maxRange = sqrt(maxGndRng^2 + rdrHeight^2);
minGate = find(rngbins >= minRange,1,'first');
maxGate = find(rngbins <= maxRange,1,'last');

Convert the azimuth pointing angles and range bins to rectangular coordinates, then use the surf function to plot the image. When specifying the azimuth domain for the image, use one more point than the number of pulses so that the image fully encompasses all azimuth angles.

az = rotrate*(0:frame)/prf;
gndrng = sqrt(rngbins(minGate:maxGate).^2-rdrHeight^2);
x = gndrng.*cosd(az);
y = gndrng.*sind(az);
shading flat
axis equal
axis tight
colormap winter
title('PPI Image')

Figure contains an axes object. The axes object with title PPI Image contains an object of type surface.

The gif below shows a recording of this scenario for 30 full rotations of the antenna over about 34 seconds.



In this example, you saw how to generate clutter and target returns with a rotating radar in a maritime environment. You saw how to use a spectral model to get realistic sea heights and surface shadowing, and how to emulate an extended target with a set of point targets, which allows for partial occlusion of the target by the surface. The IQ data was converted from polar format to Cartesian and plotted with the color channel of the surf function to create a simple PPI image.

Supporting Functions


function helperPlotScenarioGeometry(rdrHeight,seaLength,array,freq,tgtpos)
% Plot a visualization of the scenario

fh = figure;

% Surface
t = 0:10:360;
r = linspace(0,seaLength/2,10).';
x = r.*cosd(t);
y = r.*sind(t);
surf(x/1e3,y/1e3,zeros(size(x)),repmat(permute([0 1 1],[1 3 2]),[size(x) 1]),'EdgeColor','none')

hold on

% Antenna pattern
az = linspace(-6,6,80);
maxEl = -atand(2*rdrHeight/seaLength);
el = linspace(-90,maxEl,80);
G = pattern(array,freq,az,el);
[az,el] = meshgrid(az,el);
[x,y,~] = sph2cart(az*pi/180,el*pi/180,1);
rtg = rdrHeight./sind(-el);


% Radar
line([0 0],[0 0],[0 rdrHeight],'color','black','linestyle','--')

% Target

hold off

text(seaLength/4/1e3,-200/1e3,rdrHeight/10,'Beam Pattern')
text(-seaLength/4/1e3,seaLength/4/1e3,rdrHeight/10,'Sea Surface')

xlabel('X (km)')
ylabel('Y (km)')
zlabel('Z (m)')
title('Scenario Geometry')

cl = clim;
clim([cl(2)-30 cl(2)])

fh.Position = fh.Position + [0 0 150 150];
view([-11 58])