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 an extended target with a collection of point scatterers, 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 rectangular array (URA) mounted in a fixed position above a sea surface. A set of point targets arranged on the surface of a cuboid will be used to emulate a medium-sized container ship.

Define the radar system parameters. Use an X-band radar at 10 GHz with a 1 degree azimuth beamwidth and 10 meter range resolution. Use a wide 60 deg elevation beamwidth which will result in a URA with two rows of elements.

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

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

sinc3db = 0.8859;
numelems(1) = round(2*sinc3db/deg2rad(elbw));
numelems(2) = round(2*sinc3db/deg2rad(azbw));
lambda = freq2wavelen(freq);
array = phased.URA(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;
seaSurface(scenario,'Boundary',bdry,'RadarReflectivity',refl,'SpectralModel',spec,'WindSpeed',windSpeed,'WindDirection',windDir);

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)]);
platform(scenario,'Sensors',rdr,'Trajectory',traj);

Use the MountingAngles property on the radar to point the beam down so that it is centered between the radar nadir point and the edge of the sea surface.

depang = atand(4*rdrHeight/seaLength);
rdr.MountingAngles = [0 depang 0];

Define a cuboid target then use the provided helper function to add a set of discrete scatterers to the scene for a basic representation of the target. Specify the target dimensions, total RCS, position, heading, and speed. The specified range resolution is used to determine the spacing of the scatterers along the surface of the cuboid.

tgtdims = [120 18 22]; % length, width, height (m)
tgtrcs = 10;  % m^2
tgtpos = [seaLength/3 seaLength/16 0];
tgthdg = -30; % deg
tgtspd = 8;   % m/s
helperMakeExtendedTarget(scenario,tgtdims,tgtrcs,tgtpos,tgthdg,tgtspd,rngres)

Finally, enable clutter generation for the radar by calling the clutterGenerator method on the scenario. The clutter generator will include only the region of the surface that is within the 3 dB width of the beam. Surface shadowing is enabled by default, and can be disabled with the UseShadowing property.

clut = clutterGenerator(scenario,rdr);

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

helperPlotScenarioGeometry(rdrHeight,seaLength,depang,array,freq,tgtpos)

Figure contains an axes object. The axes object with title Scenario Geometry contains 9 objects of type surface, line, text.

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. One half is subtracted from the total number of pulses to ensure the simulation includes the exact number of pulses specified.

scenario.StopTime = azcov/360*(npulses-1/2)/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;
    iqsig = receive(scenario);
    [ppi(:,frame),rngbins] = resp(sum(iqsig{1},2),mfc);
end

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 a minimum ground range of 200 m and a maximum ground range equal to half the length of the sea surface, then find the range gate indices that correspond to those ranges.

minRange = sqrt(200^2 + rdrHeight^2);
maxRange = sqrt((seaLength/2)^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 360 degrees of azimuth.

az = rotrate*(0:frame)/prf;
gndrng = sqrt(rngbins(minGate:maxGate).^2-rdrHeight^2);
x = gndrng.*cosd(az);
y = gndrng.*sind(az);
surf(x,y,zeros(size(x)),20*log10(abs(ppi(minGate:maxGate,1:frame))))
shading flat
view(0,90)
colorbar
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.

PPI.gif

Conclusion

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

helperMakeExtendedTarget

function helperMakeExtendedTarget(scene,dims,totrcs,pos,hdg,spd,rngres)
% Add a set of point targets to a scenario to represent a cuboid target

scatSpacing = rngres*.8;
nScatsPerDim = ceil(dims/scatSpacing);

platB = rotz(hdg);
losrdr = -pos(:);
losrdr = platB.'*losrdr;

if losrdr(1) < -dims(1)/2
    if losrdr(2) < -dims(2)/2
        sideFaces = 'ws';
    elseif losrdr(2) < dims(2)/2
        sideFaces = 'w';
    else
        sideFaces = 'wn';
    end
elseif losrdr(1) < dims(1)/2
    if losrdr(2) < -dims(2)/2
        sideFaces = 's';
    elseif losrdr(2) < dims(2)/2
        error('Radar is inside target')
    else
        sideFaces = 'n';
    end
else
    if losrdr(2) < -dims(2)/2
        sideFaces = 'es';
    elseif losrdr(2) < dims(2)/2
        sideFaces = 'e';
    else
        sideFaces = 'en';
    end
end

xd = linspace(-dims(1)/2,dims(1)/2,nScatsPerDim(1));
yd = linspace(-dims(2)/2,dims(2)/2,nScatsPerDim(2));
zd = linspace(-dims(3)/2,dims(3)/2,nScatsPerDim(3));

x = []; y = []; z = [];
for ind = 1:numel(sideFaces)
    switch sideFaces(ind)
        case 'e'
            [y0,z0] = meshgrid(yd,zd);
            x0 = dims(1)/2*ones(size(y0));
        case 'n'
            [x0,z0] = meshgrid(xd,zd);
            y0 = dims(2)/2*ones(size(x0));
        case 'w'
            [y0,z0] = meshgrid(yd,zd);
            x0 = -dims(1)/2*ones(size(y0));
        case 's'
            [x0,z0] = meshgrid(xd,zd);
            y0 = -dims(2)/2*ones(size(x0));
    end
    x = [x, x0(:).'];
    y = [y, y0(:).'];
    z = [z, z0(:).'];
end

propRcsSide = 0.5; % proportion of the RCS to put on the sides of the target

rcs = totrcs*propRcsSide/numel(x)*ones(size(x));

% top face
[x0,y0] = meshgrid(xd,yd);
z0 = dims(3)/2*ones(size(x0));
x = [x, x0(:).'];
y = [y, y0(:).'];
z = [z, z0(:).'];

rcs = [rcs, totrcs*(1-propRcsSide)/numel(x0)*ones(1,numel(x0))];

% Rotate
[x,y,z] = deal(platB(1,1)*x+platB(1,2)*y+platB(1,3)*z,...
               platB(2,1)*x+platB(2,2)*y+platB(2,3)*z,...
               platB(3,1)*x+platB(3,2)*y+platB(3,3)*z);

% Translate
x = x + pos(1);
y = y + pos(2);
z = z + pos(3) + dims(3)/4;

for ind = 1:numel(x)
    traj = kinematicTrajectory('Position',[x(ind) y(ind) z(ind)],'Velocity',spd*[cosd(hdg) sind(hdg) 0]);
    platform(scene,'Trajectory',traj,'Signatures',rcsSignature('Pattern',10*log10(rcs(ind))));
end

end

helperPlotScenarioGeometry

function helperPlotScenarioGeometry(rdrHeight,seaLength,depang,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+depang);
[az,el] = meshgrid(az,el);
[x,y,~] = sph2cart(az*pi/180,el*pi/180,1);
rtg = rdrHeight./sind(-el);

surf(x.*rtg/1e3,y.*rtg/1e3,1e-2*ones(size(x)),G,'EdgeColor','none')

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

% Target
plot3(tgtpos(1)/1e3,tgtpos(2)/1e3,tgtpos(3),'*r')

hold off

text(0,0,rdrHeight*1.1,'Radar')
text(tgtpos(1)*.8/1e3,tgtpos(2)*1.1/1e3,rdrHeight/10,'Target')
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])


end