simulating radar returns from moving sea surfaces -凯发k8网页登录

maritime radar systems operate in a challenging, dynamic environment. to improve detection of targets of interest and assess system performance, it is important to understand the nature of returns from sea surfaces.

in this example, you will simulate an x-band radar system used for oceanographic studies of sea states. the radar system is a fixed, offshore platform. you will first generate a moving sea surface at sea state 4 using the elfouhaily spectrum model. you will then generate iq returns from the sea surface and investigate the statistics and time-frequency behavior of the simulated sea surface signals.

define radar system parameters

define an x-band radar system with a range resolution approximately equal to 30 meters. verify the range resolution is as expected using the bw2rangeres function.

rng(2021)       % initialize random number generator
% radar parameters
freq = 9.39e9;  % center frequency (hz)
prf = 500;      % prf (hz)
tpd = 100e-6;   % pulse width (s)
azbw = 2;       % azimuth beamwidth (deg)
depang = 30;    % depression angle (deg)
antgain = 45.7; % antenna gain (db)
fs = 10e6;      % sampling frequency (hz)
bw = 5e6;       % waveform bandwidth (hz)
bw2rangeres(bw)

ans = 29.9792

next, define the sea surface as a sea state 4 with a wind speed of 19 knots and a resolution of 2 meters. set the length of the sea surface to 512 meters.

% sea surface
seastate = 4;   % sea state number 
vw = 19;        % wind speed (knots)
l = 512;        % sea surface length (m) 
ressurface = 2; % resolution of sea surface (m)
% calculate wavelength and get speed of light
[lambda,c] = freq2wavelen(freq);
% setup sensor trajectory and simulation times
rdrht = 100;                                   % radar platform height (m) 
rdrpos = [-l/2 0 rdrht];                       % radar position [x y z] (m)
numpulses = 1024;                              % number of pulses 
scenarioupdatetime = 1/prf;                    % scenario update time (sec) 
scenarioupdaterate = prf;                      % scenario update rate (hz)
simtime = scenarioupdatetime.*(numpulses - 1); % overall simulation time (sec)

model the sea surface

define a scenario with radarscenario. next, add sea motion with seaspectrum. keep the spectrumsource property default auto to use an elfouhaily spectrum. set the resolution to 2 meters. the elfouhaily model consists of a sea spectrum and an angular spreading function. the sea spectrum model gets physical properties like wind speed and fetch from the seasurface, as you will see later in this section.

you can also import your own custom model by setting the seaspectrum spectrumsource property to custom. alternative sea spectrum models in the literature include the jonswap, donelan-pierson, and pierson-moskovitz spectrums.

% create scenario
scene = radarscenario('updaterate',scenarioupdaterate, ...
    'isearthcentered',false,'stoptime',simtime);
% define elfouhaily sea spectrum 
seaspec = seaspectrum('resolution',ressurface) % see elfouhaily reference

seaspec = 
  seaspectrum with properties:
    spectrumsource: 'auto'
        resolution: 2

now, select a reflectivity model. radar toolbox™ offers 9 different reflectivity models for sea surfaces covering a wide range of frequencies, grazing angles, and sea states. the asterisk denotes the default model. type help surfacereflectivitysea or doc surfacereflectivitysea in the command window for more information on usage and applicable grazing angles for each model. the sea reflectivity models are as follows.

apl: semi-empirical model supporting low grazing angles over frequencies in the range from 1 to 100 ghz. applicable for sea states 1 through 6.
git: semi-empirical model supporting low grazing angles over frequencies in the range from 1 to 100 ghz. applicable for sea states 1 through 6.
hybrid: semi-empirical model supporting low to medium grazing angles over frequencies in the range from 500 mhz to 35 ghz. applicable for sea states 0 through 5.
masuko: empirical model supporting low grazing angles over frequencies in the range from 8 to 40 ghz. applicable for sea states 1 through 6.
nathanson: empirical model supporting low grazing angles over frequencies in the range from 300 mhz to 36 ghz. applicable for sea states 0 through 6.
nrl*: empirical model supporting low grazing angles over frequencies in the range from 500 mhz to 35 ghz. applicable for sea states 0 through 6.
rre: mathematical model supporting low grazing angles over frequencies in the range from 9 to 10 ghz. applicable for sea states 1 through 6.
sittrop: empirical model supporting low grazing angles over frequencies in the range from 8 to 12 ghz. applicable for sea states 0 through 7.
tsc: empirical model supporting low grazing angles over frequencies in the range from 500 mhz to 35 ghz. applicable for sea states 0 through 5.

for this example, set the reflectivity to the git (georgia institute of technology) model. for a general discussion about sea reflectivity models and how they are affected by frequency, grazing angle, sea state, and polarization, see the maritime radar sea clutter modeling example.

% define reflectivity model
pol = 'v'; % polarization 
reflectmodel = surfacereflectivity('sea','model','git','seastate',seastate,'polarization',pol)

reflectmodel = 
  surfacereflectivitysea with properties:
    enablepolarization: 0
                 model: 'git'
              seastate: 4
          polarization: 'v'
               speckle: 'none'

add a sea surface to the radar scenario using seasurface. assign the previously defined sea spectrum and reflectivity model to the sea surface. seasurface houses the physical properties of the surface. the windspeed, winddirection, and fetch properties are used in generating the elfouhaily sea spectrum. by thoughtful selection of these properties and the resolution of the sea surface, gravity and capillary waves with different wave ages can be produced.

capillary waves are high-frequency, small-wavelength waves that have a rounded crest with a v-shaped trough that are less than 2 cm. gravity waves are longer, low-frequency waves. the seaspectrum object cannot model breaking waves, which are waves that reach a critical amplitude resulting in wave energy being transformed into turbulent, kinetic energy. the greater the wind speed, the greater the energy transferred to the sea surface and the larger the waves.

fetch is the length of water over which wind blows without obstruction. the elfouhaily model relates the fetch to the inverse wave age $ω_{c}$ as

$ω_{c} = 0.84 \tanh {({(\frac{x}{x_{0}})}^{0.4})}^{- 0.75} .$

$x$ is calculated as

$k_{0} = \frac{g}{u_{10}}$ and $x = k_{0} fetch$ ,

where $g$ is the acceleration due to gravity in meters per second squared, $u_{10}$ is the wind speed at 10 meters above the surface in m/s, and $x_{0}$ is a constant set to 2.2e4. a fully-developed sea has an inverse wave age of 0.84, a mature sea has a value of approximately 1, and a young sea has a value between 2 and 5.

the default fetch for seasurface is infinity, which produces a fully-developed sea. more information about the elfouhaily spectrum can be found in the reference listed at the end of this example, as well as in the documentation for seaspectrum.

the elevation values of the sea surface can be obtained by using the height method on the sea surface. plot the radar and the sea surface using .

% configure sea surface
knots2mps = 0.514444; % knots to meters/sec
vw = vw*knots2mps; % wind speed (m/s)
seasurf = seasurface(scene,'spectralmodel',seaspec,'radarreflectivity',reflectmodel, ...
    'windspeed',vw,'winddirection',0,'boundary',[-l/2 l/2; -l/2 l/2])

seasurf = 
  seasurface with properties:
            windspeed: 9.7744
        winddirection: 0
                fetch: inf
        spectralmodel: [1x1 seaspectrum]
    radarreflectivity: [1x1 surfacereflectivitysea]
      reflectivitymap: 1
      referenceheight: 0
             boundary: [2x2 double]

% plot sea surface and radar 
x = -l/2:ressurface:(l/2 - 1); 
y = -l/2:ressurface:(l/2 - 1); 
[xgrid,ygrid] = meshgrid(x,y);
z = height(seasurf,[xgrid(:).'; ygrid(:).'],scene.simulationtime); 
helperseasurfaceplot(x,y,z,rdrpos)

figure contains an axes object. the axes object with title sea surface elevations, xlabel x (m), ylabel y (m) contains 2 objects of type surface, line. one or more of the lines displays its values using only markers these objects represent sea surface, radar.

investigate the statistics of the elevation values of the sea surface. first, calculate and plot the cumulative distribution function using . you can use this plot to determine the probability that an observation will be less than or equal to a certain elevation value. for example, 90% of elevations will be less than or equal to about 1 meter. 95% of elevations will be less than or equal to about 1.5 meters. if you would like to alter the statistics of the generated sea surface to increase the elevations, you can increase the wind speed to deliver more energy to the surface.

% cdf
helperseasurfacecdf(z)

figure contains an axes object. the axes object with title wave elevation cdf, xlabel wave elevation (m), ylabel probability contains an object of type line.

next, estimate the significant wave height using . the estimate is obtained by taking the mean of the highest 1/3 of the wave heights, where the wave height is defined from trough-to-crest as is visually represented in the figure.

% significant wave height
actsighgt = helperestimatesignificantwaveheight(x,y,z)

figure contains an axes object. the axes object with title wave height estimation, xlabel x (m), ylabel elevation (m) contains 3 objects of type line. one or more of the lines displays its values using only markers these objects represent wave elevation data, trough, crest.

actsighgt = 1.6340

the expected wave heights for a sea state 4 is in the range between 1.25 and 2.5 meters. note that the simulated sea surface is within the range of expected values.

expectedsighgt = [1.25 2.5]; % sea state 4
actsighgt >= expectedsighgt(1) && actsighgt <= expectedsighgt(2)

ans = logical
   1

lastly, plot a subset of the sea surface heights over time to visualize the motion due to the elfouhaily sea spectrum. note that the simulation time for the radar data collection only runs for 2 seconds, but plot over a time period of 10 seconds to get a better sense of the motion. uncomment to plot the moving sea surface.

% plot sea surface motion 
plottime = 0:0.5:10; % plot time (sec)
% helperseasurfacemotionplot(x,y,seasurf,plottime); % uncomment to plot motion

configure the radar transceiver

in this section, configure the radar system properties. define the antenna and the transmitted linear frequency modulated (lfm) waveform. assign the radar sensor to the radar platform.

% create a radar platform using the trajectory information
rdrplat = platform(scene,'position',rdrpos);
% create a radar sensor looking to the east
rdr = radartransceiver('mountingangles',[0 depang 0],'numrepetitions',1);
% configure the lfm signal of the radar
rdr.waveform = phased.linearfmwaveform('samplerate',fs,'pulsewidth',tpd, ...
    'prf',prf,'sweepbandwidth',bw);
% set receiver sample rate and noise figure
rdr.receiver.samplerate = fs;
rdr.receiver.noisefigure = 10; 
% define transmit and receive antenna and corresponding parameters
ant = phased.sincantennaelement('beamwidth',azbw);
rdr.transmitantenna.operatingfrequency = freq;
rdr.receiveantenna.operatingfrequency = freq;
rdr.transmitter.gain = antgain;
rdr.receiver.gain = antgain;
rdr.transmitantenna.sensor = ant;
rdr.receiveantenna.sensor = ant;
% add radar to radar platform
rdrplat.sensors = rdr;

generate the datacube

now that the scene and the radar system have been defined, generate returns from the sea surface with the cluttergenerator method and collect the iq data using the receive method, which both reside in the radarscenario object. by default, cluttergenerator simulates clutter returns in the mainlobe. for more information about clutter modeling, see introduction to radar scenario clutter simulation.

use the helpers and to plot the unprocessed, raw iq data.

% collect clutter returns with the cluttergenerator
cluttergenerator(scene,rdr); 
% run the scenario
numsamples = 1/prf*fs;
maxrange = 20e3; 
trng = (0:1/fs:(numsamples-1)/fs);
rnggrid = [0 time2range(trng(2:end),c)];
[~,idxtruncate] = min(abs(rnggrid - maxrange));
iqsig = zeros(idxtruncate,numpulses);
ii = 1; 
hraw = helperplotrawiq(iqsig);
while advance(scene)
    tmp = receive(scene); % nsamp x 1
    iqsig(:,ii) = tmp{1}(1:idxtruncate); 
    if mod(ii,128) == 0
        helperupdateplotrawiq(hraw,iqsig);
    end
    ii = ii   1;
end
helperupdateplotrawiq(hraw,iqsig);

figure contains an axes object. the axes object with title raw iq, xlabel pulses, ylabel samples contains an object of type surface.

process sea surface returns

use phased.rangeresponse to pulse compress the returned iq data. visualize the range-time behavior using . if the sea surface was static, you would see straight, horizontal lines in the plot, but the sea returns exhibit dynamism. at a given range, the magnitude changes with time. note that the returns occur in a small subset of range due to the geometry (radar height, beamwidth, and depression angle).

% pulse compress
matchingcoeff = getmatchedfilter(rdr.waveform);
rngresp = phased.rangeresponse('rangemethod', 'matched filter', ...
    'propagationspeed',c,'samplerate',fs);
[pcresp,rnggrid] = rngresp(iqsig,matchingcoeff); 
% plot
pcsigmagdb = mag2db(abs(pcresp));
[maxval,maxidx] = max(pcsigmagdb(:));
pcsigmagdb = pcsigmagdb - maxval;
helperrangetimeplot(rnggrid,prf,pcsigmagdb);

figure contains an axes object. the axes object with title range versus time, xlabel time (sec), ylabel range (m) contains an object of type surface.

use the range index of the maximum value and visualize the magnitude versus time at this range bin using .

% plot magnitude versus time
[idxrange,~] = ind2sub(size(pcsigmagdb),maxidx); 
helpermagtimeplot(pcsigmagdb(idxrange,:),numpulses,prf,rnggrid(idxrange));

figure contains an axes object. the axes object with title magnitude versus time at range 194.87 (m), xlabel time (sec), ylabel magnitude (db) contains an object of type line.

next, generate a time-frequency plot of that same range bin using . generate the plot using the short-time fourier transform function stft.

% stft
[s,f,t] = stft(pcresp(idxrange,:),scenarioupdaterate);
helperstftplot(s,f,t,lambda,rnggrid(idxrange));

figure contains an axes object. the axes object with title stft at range 194.87 (m), xlabel time (sec), ylabel speed (m/s) contains an object of type surface.

notice that the detected speed of the range bin changes over time due to the motion of the sea surface. over longer simulation times, periodicities are detectable.

lastly, form a histogram from the magnitude values of the data in the ranges between 180 and 210 meters using . note the shape of the histogram resembles a weibull distribution. if you have the statistics and machine learning toolbox™, you can use the fitdist function to get the parameters for the weibull distribution from the magnitude of the data.

% look at a subset of data in range and convert to decibel scale
[~,idxmin] = min(abs(rnggrid - 180)); 
[~,idxmax] = min(abs(rnggrid - 210)); 
pcsigmagdb = mag2db(abs(pcresp(idxmin:idxmax,:)));
% remove any inf values
pcsigmagdb(isinf(pcsigmagdb(:))) = []; 
% shift values to be positive
pcsigmagdb = pcsigmagdb(:) - min(pcsigmagdb(:))   eps; 
% weibull parameters
% note: these values can be obtained if you use fitdist(pcsigmagdb,'weibull')
scale = 37.8589;
shape = 7.80291;
% plot histogram with theoretical pdf overlay
helperhistplot(pcsigmagdb,scale,shape);

figure contains an axes object. the axes object with title amplitude histogram, xlabel normalized magnitude (db), ylabel pdf contains 2 objects of type histogram, line.

summary

in this example, you learned how to:

build a radar scenario,
model a moving sea surface,
investigate sea surface statistics and behavior,
generate iq, and
perform time-frequency processing on simulated iq returns.

this example can easily be extended to include targets to support applications ranging such as maritime surveillance and radar imaging.

references

apel, john r. principles of ocean physics. international geophysics series, v. 38. london; orlando: academic press, 1987.
elfouhaily, t., b. chapron, k. katsaros, and d. vandemark. "a unified directional spectrum for long and short wind-driven waves." journal of geophysical research: oceans 102, no. c7 (july 15, 1997): 15781-96. .
doerry, armin w. "ship dynamics for maritime isar imaging ." sandia report. albuquerque, new mexico: sandia national laboratories, february 2008.
ward, keith d., simon watts, and robert j. a. tough. sea clutter: scattering, the k-distribution and radar performance. iet radar, sonar, navigation and avionics series 20. london: institution of engineering and technology, 2006.

supporting functions

helperseacolormap

function cmap = helperseacolormap(n)
% color map for sea elevation plot
c = hsv2rgb([2/3 1 0.2; 2/3 1 1; 0.5 1 1]);
cmap = zeros(n,3);
cmap(:,1) = interp1(1:3,c(:,1),linspace(1,3,n)); 
cmap(:,2) = interp1(1:3,c(:,2),linspace(1,3,n));
cmap(:,3) = interp1(1:3,c(:,3),linspace(1,3,n)); 
end

helperseasurfaceplot

function helperseasurfaceplot(x,y,z,rdrpos)
% plot sea elevations
% color map for sea elevation plot
seacolormap = helperseacolormap(256);
% plot
figure
z = reshape(z,numel(x),numel(y)); 
surf(x,y,z)
hold on
plot3(rdrpos(1),rdrpos(2),rdrpos(3),'ok','linewidth',2,'markerfacecolor','k')
legend('sea surface','radar','location','best')
shading interp
axis equal
xlabel('x (m)')
ylabel('y (m)')
zlabel('elevations (m)')
stdz = std(z(:));
minc = -4*stdz;
maxc = 4*stdz;
minz = min(minc,rdrpos(3)); 
maxz = max(maxc,rdrpos(3)); 
title('sea surface elevations')
axis tight
zlim([minz maxz])
hc = colorbar('southoutside');
colormap(seacolormap)
hc.label.string = 'elevations (m)';
hc.limits = [minc maxc];
drawnow
pause(0.25)
end

helperseasurfacecdf

function helperseasurfacecdf(x)
% calculate and plot cdf 
% form cdf
x = x(~isnan(x));
n = length(x);
x = sort(x(:));
yy = (1:n)' / n;
notdup = ([diff(x(:)); 1] > 0);
xx = x(notdup);
yy = [0; yy(notdup)];
% create vectors for plotting
k = length(xx);
n = reshape(repmat(1:k, 2, 1),2*k,1);
xcdf = [-inf; xx(n); inf];
ycdf = [0; 0; yy(1 n)];
% plot
figure
plot(xcdf,ycdf,'linewidth',2)
grid on
title('wave elevation cdf') 
xlabel('wave elevation (m)')
ylabel('probability')
drawnow
pause(0.25)
end

helperestimatesignificantwaveheight

function sighgt = helperestimatesignificantwaveheight(x,y,z)
% plot an example showing how the wave heights are estimated and estimate
% the significant wave height 
% plot wave height estimation
figure 
numx = numel(x);
z = reshape(z,numx,numel(y));
zest = z(numx/2   1,:);
plot(x,zest,'linewidth',2)
grid on
hold on
xlabel('x (m)')
ylabel('elevation (m)')
title('wave height estimation')
axis tight
idxmin = islocalmin(z(numel(x)/2   1,:));
idxmax = islocalmax(z(numel(x)/2   1,:));
co = colororder; 
plot(x(idxmin),zest(idxmin),'v', ...
    'markerfacecolor',co(2,:),'markeredgecolor',co(2,:))
plot(x(idxmax),zest(idxmax),'^', ...
    'markerfacecolor',co(3,:),'markeredgecolor',co(3,:))
legend('wave elevation data','trough','crest','location','southwest')
% estimate significant wave height
wavehgts = [];
for ii = 1:numx
    zest = z(ii,:);
    idxmin = islocalmin(zest);
    troughs = zest(idxmin);
    numtroughs = sum(double(idxmin)); 
    idxmax = islocalmax(zest);
    crests = zest(idxmax);
    numcrests = sum(double(idxmax)); 
    numwaves = min(numtroughs,numcrests); 
    wavehgts = [wavehgts ...
        abs(crests(1:numwaves) - troughs(1:numwaves))]; %#ok 
end
wavehgts = sort(wavehgts); 
idxtopthird = floor(numel(wavehgts)*2/3);
sighgt = mean(wavehgts(idxtopthird:end));
drawnow
pause(0.25)
end

helperseasurfacemotionplot

function helperseasurfacemotionplot(x,y,seasurf,plottime) %#ok 
% color map for sea elevation plot
seacolormap = helperseacolormap(256);
% get initial height 
[xgrid,ygrid] = meshgrid(x,y);
z = height(seasurf,[xgrid(:).'; ygrid(:).'],plottime(1)); 
% plot
hfig = figure;
hs = surf(x,y,reshape(z,numel(x),numel(y)));
hold on
shading interp
axis equal
xlabel('x (m)')
ylabel('y (m)')
zlabel('elevations (m)')
stdz = std(z(:));
minz = -4*stdz; 
maxz = 4*stdz; 
title('sea surface motion plot')
xlim([-50 50])
ylim([-50 50])
zlim([minz maxz])
hc = colorbar('southoutside');
colormap(seacolormap)
hc.label.string = 'elevations (m)';
hc.limits = [minz maxz];
view([-45 12])
% change figure size
ppos0 = get(hfig,'position');
figwidth = 700;
figheight = 300;
set(gcf,'position',[ppos0(1),ppos0(2),figwidth,figheight])
numtimes = numel(plottime);
for ii = 2:numtimes
    % update plot
    z = height(seasurf,[xgrid(:).'; ygrid(:).'],plottime(ii));
    hs.zdata = reshape(z,numel(x),numel(y));
    drawnow
    pause(0.1)
end
pause(0.25)
end

helperplotrawiq

function hraw = helperplotrawiq(iqsig)
% plot raw iq magnitude
figure()
hraw = pcolor(mag2db(abs(iqsig)));
hraw.edgecolor = 'none';
title('raw iq')
xlabel('pulses')
ylabel('samples')
hc = colorbar;
hc.label.string = 'magnitude (db)'; 
drawnow
pause(0.25)
end

helperupdateplotrawiq

function helperupdateplotrawiq(hraw,iqsig)
% update the raw iq magnitude plot
hraw.cdata = mag2db(abs(iqsig)); 
drawnow
pause(0.25)
end

helperrangetimeplot

function helperrangetimeplot(rnggrid,prf,pcsigmagdb)
% range-time plot
figure() 
numpulses = size(pcsigmagdb,2); 
himage = pcolor((1:numpulses)*1/prf,rnggrid,pcsigmagdb);
himage.edgecolor = 'none';
shading interp
xlabel('time (sec)')
ylabel('range (m)')
hc = colorbar;
hc.label.string = 'magnitude (db)'; 
axis tight
title('range versus time')
caxis([-20 0])
ylim([140 260])
drawnow
pause(0.25)
end

helpermagtimeplot

function helpermagtimeplot(magvals,numpulses,prf,rngval)
% magnitude vs. time plot
figure() 
plot((1:numpulses)*1/prf,magvals,'linewidth',2)
grid on
xlabel('time (sec)')
ylabel('magnitude (db)')
axis tight
title(sprintf('magnitude versus time at range %.2f (m)',rngval))
ylim([-20 0])
drawnow
pause(0.25)
end

helperstftplot

function helperstftplot(s,f,t,lambda,rngval)
% time-frequency plot 
figure()
s = mag2db(abs(s)); % convert to db
s = s - max(s(:)); % normalize
speed = dop2speed(f,lambda)/2; % speed (m/s)
hstft = pcolor(t,speed,s);
hold on
hstft.edgecolor = 'none';
colorbar
xlabel('time (sec)')
ylabel('speed (m/s)')
title(sprintf('stft at range %.2f (m)',rngval))
caxis([-20 0])
ylim([-1 1])
drawnow
pause(0.25)
end

helperhistplot

function helperhistplot(pcsigmagdb,scale,shape)
% amplitude distribution
% histogram
figure
hhist = histogram(pcsigmagdb,'normalization','pdf'); % create histogram
grid on
hold on
xlabel('normalized magnitude (db)')
ylabel('pdf')
title('amplitude histogram')
% lognormal overlay
edges = hhist.binedges;
x = (edges(2:end)   edges(1:end-1))./2;
z = x./scale;
w = exp(-(z.^shape));
y = z.^(shape - 1).*w .*shape./scale;
plot(x,y,'-r')
end