Simulating Radar Signals with Atmospheric Refraction Effects

The performance of a radar system is closely tied to the environment in which it operates. This example will discuss one of the environmental factors that is outside of the system designer's control, namely the atmosphere. In particular, this example will focus on the refraction effects due to the atmosphere and its impact on target parameter estimation.

This example demonstrates errors due to refraction for an L-band surveillance radar operating at a high latitude during winter. It first configures a radar scenario and reviews built-in atmospheric models available in Radar Toolbox™. Next it presents a simple exponential approximation for the refractivity versus height. It then uses these models to generate IQ using radarScenario and processes the IQ to generate detections. Lastly, the detections are then used in an analysis showing the refraction errors for the specified geometry.

Define Radar Parameters

Define an L-band monostatic surveillance radar that has coherent processing intervals (CPIs) that are 24 pulses long. The sensor is a stationary system mounted to a tower 10 meters above the surface.

% Set random number generator to ensure results are repeatable
rng('default')

% Sensor parameters
freq       = 2e9;  % L-band center frequency (Hz)
numPulses  = 24;   % Number of pulses in a CPI
prf        = 1e3;  % PRF (Hz)
tpd        = 1e-5; % Pulse width (s)
fs         = 3e7;  % Sampling frequency (Hz)
bw         = fs/2; % Bandwidth (Hz)
sensorHgt  = 10;   % Sensor height above the surface (m)

Given the bandwidth, determine the range resolution using the bw2rangeres function.

bw2rangeres(bw)    % Range resolution (m)

ans = 9.9931

Configure the Radar Scenario

First configure the target positions for the analysis. Within this scenario, 21 different target positions will be considered. Set the target elevation angles to a constant 2 degrees with true slant ranges (also called the geometric range) varying from 1 to 30 km. Set the target radar cross section (RCS) to 20 dBsm, which is representative of a large aircraft.

% Target parameters
numTgtPos  = 21; % Target positions
tgtRCSdBsm = 20; % Target RCS (dBsm)
tgtEl      = 2;  % Target elevation (deg)
tgtSR      = linspace(1,30,numTgtPos)*1e3; % Target slant range (m)

Calculate the true target heights over a curved Earth using the range2height function with an effective Earth radius set to the true Earth radius. Use the target trajectory to define an appropriate update time for the radar scenario.

% Calculate target heights
Rearth   = physconst('EarthRadius'); % True Earth radius (m)
tgtHgt   = range2height(tgtSR,sensorHgt,tgtEl, ...
    'EffectiveEarthRadius',Rearth);  % Target heights (m)
tgtSpeed = 90;                       % Target speed (m/s) 
updateTime = (tgtSR(2) - tgtSR(1))/tgtSpeed; 

Calculate the flight duration and create a radar scenario.

flightDuration = (numTgtPos - 1).*updateTime;
updateRate = 1/updateTime; 
scene = radarScenario('UpdateRate',updateRate, ...
    'IsEarthCentered',true,'StopTime',flightDuration);

Add the Atmosphere

Simulating the refraction phenomenon requires a model of the atmosphere. For this example, the built-in atmospheres available in the refractiveidx function will be used. The refractiveidx function offers six ITU reference atmospheric models, namely:

• Standard atmospheric model, also known as the Mean Annual Global Reference Atmosphere (MAGRA)

• Summer high latitude model (higher than 45 degrees)

• Winter high latitude model (higher than 45 degrees)

• Summer mid latitude model (between 22 and 45 degrees)

• Winter mid latitude model (between 22 and 45 degrees)

• Low latitude model (smaller than 22 degrees)

Select the high latitude model in winter. More information about these atmospheric models can be found in the Modeling Target Position Errors Due to Refraction example.

The Central Radio Propagation Laboratory developed a widely used reference atmosphere model that approximates the refractivity profile as an exponential decay versus height. In general, the CRPL model is a good approximation for refractivity profiles under normal propagation conditions. This exponential decay model underlies the CRPL method in the following functions:

• height2range

• height2grndrange

• range2height

• slant2range

• refractionexp

It also exists as an atmosphere option within a radar scenario.

The CRPL exponential model is given by

$$\mathit{N}={\mathit{N}}_{\mathit{S}}\exp \left(-{\mathit{c}}_{\mathit{e}}\mathit{h}\right),$$

where ${\mathit{N}}_{\mathit{S}}$ is the surface refractivity (N-units), $\mathit{h}$ is the height (km), and ${\mathit{c}}_{\mathit{e}}$ is a decay constant (1/km) calculated as

${\mathit{c}}_{\mathit{e}}=\ln \frac{{\mathit{N}}_{\mathit{S}}}{{\mathit{N}}_{\mathit{S}}-\Delta \mathit{N}}$.

In the above equation, $\Delta \mathit{N}$ is the difference between the surface refractivity ${\mathit{N}}_{\mathit{S}}$ and the refractivity at a height of 1 km. $\Delta \mathit{N}$ can be approximated by an exponential function expressed as

$\Delta \mathit{N}=7.32\exp \left(0.005577{\mathit{N}}_{\mathit{S}}\right)$.

The CRPL exponential model is easily modified to match local conditions by setting appropriate surface refractivity and refractive exponent values. Use the surface refractivity output from the refractiveidx function as an input to the CRPL model. Calculate the decay constant (refraction exponent) using the refractionexp function.

% Define CRPL model parameters
[~,N0] = refractiveidx(0,'LatitudeModel','High','Season','Winter');  % Surface refractivity (N-units)
rexp   = refractionexp(N0);                                          % Corresponding refraction exponent (1/km)

Define the atmosphere model of the radar scenario as CRPL. Set the surface refractivity and refraction exponent to the previously obtained values.

atmosphere(scene,'CRPL','SurfaceRefractivity',N0,'RefractionExponent',rexp)

ans = 
  AtmosphereCRPL with properties:

    SurfaceRefractivity: 313.4066
     RefractionExponent: 0.1440
       MaxNumIterations: 10
              Tolerance: 0.0100

Create a Radar Transceiver

First define a sensor trajectory using geoTrajectory and add the trajectory to a platform.

% Surface parameters
surfHgt    = 52;  % Surface height at sensor (m)

% Setup sensor trajectory
sensoralt  = sensorHgt + surfHgt; % Surface altitude (m) 
sensorpos1 = [42.30141195164364,-71.376098592854,sensoralt]; % Natick, MA USA [Lat (deg) Lon (deg) Alt (m)] 
sensorpos2 = sensorpos1;
sensorTraj = geoTrajectory([sensorpos1;sensorpos2],[0;flightDuration]);

% Add sensor platform to scenario
sensorplat = platform(scene,'Trajectory',sensorTraj);

Next configure the monostatic radar sensor using radarTransceiver. Define the antenna as a rectangular phased array of size 4-by-4 with half-wavelength spacing. Set the waveform modulation to be an LFM.

% Configure radarTransceiver
[lambda,c] = freq2wavelen(freq);
elementSpacing = lambda/2;
uraSize = [4 4];
ant = phased.URA('Size',[4 4],'ElementSpacing',elementSpacing);
sensor = radarTransceiver('MountingAngles',[-90 0 0],'NumRepetitions',numPulses);
sensor.Waveform = phased.LinearFMWaveform('SampleRate',fs,'PulseWidth',tpd, ...
    'PRF',prf,'SweepBandwidth',bw);
sensor.Transmitter.Gain = 30;
sensor.Receiver.Gain = 30;
sensor.Receiver.SampleRate = fs;
sensor.Receiver.NoiseFigure = 0;
sensor.TransmitAntenna.Sensor = ant;
sensor.ReceiveAntenna.Sensor = ant;
sensor.TransmitAntenna.OperatingFrequency = freq;
sensor.ReceiveAntenna.OperatingFrequency = freq;

Attach the sensor to the platform.

sensorplat.Sensors = sensor;

Define Target Trajectory

Convert target range and angle to geodetic coordinates using aer2geodetic.

% Define a structure for the target truth
tgtTruth.Range = tgtSR;
tgtTruth.Elevation = tgtEl*ones(1,numTgtPos);
tgtTruth.Height = tgtHgt;

% Convert azimuth, elevation, and range to geodetic positions
tgtpos     = zeros(numTgtPos,3);       % Target [lat lon alt]
tgtposTime = zeros(numTgtPos,1);       % Position times
refS       = referenceSphere('earth'); % Earth model 
for m = 1:numTgtPos   
    [tgtLat,tgtLon,tgtalts] = aer2geodetic(90,tgtEl,tgtSR(m), ...
        sensorpos1(1),sensorpos1(2),sensoralt,refS);
    tgtpos(m,:) = [tgtLat tgtLon tgtalts];
    tgtposTime(m) = (m - 1).*updateTime;
end

Define the target trajectory using geoTrajectory and attach to a target platform.

% Add target platform to scenario
tgtTraj = geoTrajectory(tgtpos,tgtposTime);
platform(scene,'Trajectory',tgtTraj);

Visualize the radar and target geometry by using the helper helperRadarGlobeViewer. The radar platform is shown as the blue dot and the target platform is shown as the red dot. The white line indicates the constant elevation target trajectory.

% Visualize the scenario
gv = helperRadarGlobeViewer; 
showScenario(gv,scene)
positionCamera(gv,[42.1417 -71.3920 1.4668e+03],[359.9981 5.3302 31.7642])

Setup a Signal Processing Chain

In this section, the objects for a signal processing chain will be configured. At a high-level, the signal processing chain proceeds as depicted.

1. First, the simulated IQ is range and Doppler processed.

2. Next, the data is beamformed to broadside.

3. The beamformed data is then processed by a 2-dimensional constant false alarm rate detector (CFAR) to generate detections. Only the highest SNR detection is retained.

4. The azimuth and elevation angles of the target are then estimated using a sum/difference beamformer.

5. Lastly, the range of the target is further refined by quadratic interpolation. The output of this chain is the target detection.

First, configure phased.RangeDopplerResponse to create an object to perform range and Doppler processing on the simulated datacube.

% Initialize range/Doppler processing objects
matchingCoeff = getMatchedFilter(sensor.Waveform);
rngdopresp = phased.RangeDopplerResponse('SampleRate',fs, ...
    'DopplerFFTLengthSource','Property','DopplerFFTLength',numPulses, ...
    'DopplerOutput','Speed','PRFSource','Property','PRF',prf);

Define a steering vector object that will be used to create a range-angle plot, as well as to beamform prior to performing CFAR detection.

% Beamforming 
steeringvec = phased.SteeringVector('SensorArray',ant);

% Create steering vector for initial range-angle visualization
elVec = -90:90;
sv = steeringvec(freq,[zeros(1,numel(elVec)); elVec]);

Next, define a 2-D CFAR with 3 guard cells and 15 training cells on each side of the cell under test in range and 1 guard cell and 3 training cells in Doppler. Set the false alarm rate to 1e-6.

% 2-D CFAR
cfar = phased.CFARDetector2D('GuardBandSize',[3 1], ...
    'TrainingBandSize',[15 3], ...
    'ProbabilityFalseAlarm',1e-6,'OutputFormat','Detection index','NoisePowerOutputPort',true);

Define a sum/difference beamformer for angle estimation and create a range estimator object for improved range estimation.

% Angle estimation
sumdiff = phased.SumDifferenceMonopulseTracker2D('SensorArray',ant,'OperatingFrequency',freq);

% Range estimation object
rangeestimator = phased.RangeEstimator;

Simulate IQ and Create Detections

Simulating the IQ is as simple as calling the receive method on the radar scenario scene. Advance the scenario and collect IQ from radarTransceiver.

% Initialize simulation
numSamples = 1/prf*fs;
dets = struct('Range',zeros(1,numTgtPos), ...
    'Angles',zeros(2,numTgtPos), ...
    'SNRdB',zeros(1,numTgtPos));
for m = 1:numTgtPos
    % Advance scene
    advance(scene);
    
    % Update geometry visualization
    updateDisplay(gv,scene.SimulationTime,[scene.Platforms{:}]);  
    positionCamera(gv,[42.1417 -71.3920 1.4668e+03],[359.9981 5.3302 31.7642]);
    drawnow; 
    pause(0.1);

    % Collect IQ
    tmp = receive(scene); 
    iqsig = tmp{1}; 

Range and Doppler process the simulated IQ.

    % Range-Doppler process
    [resp,rngVec,dopVec] = rngdopresp(iqsig,matchingCoeff);

Next, beamform the data to get a sense of the range-angle performance.

    % Beamform across all angles
    respAng = (sv'*squeeze(max(resp,[],3)).');

Plot the range-angle and range-Doppler data.

    % Convert to dB
    angRespMagdB = mag2db(abs(respAng));
    respMagdB = mag2db(abs(resp));

    if m == 1
        % Initialize range-angle plot
        figure
        tiledlayout(1,2)
        nexttile
        hRAMaxes = gca;
        hRAMplot = pcolor(rngVec*1e-3,elVec,angRespMagdB);
        hRAMplot.EdgeColor = 'none';
        hRAMplot.FaceColor = 'interp';
        ylabel('Elevation Angles (deg)')
        xlabel('Range (km)')
        title('Range-Angle')
        hC = colorbar;
        hC.Label.String = 'Magnitude (dB)';

        % Initialize range-Doppler plot
        nexttile
        hRDMaxes = gca;
        hRDMplot = pcolor(rngVec*1e-3,dopVec,squeeze(max(respMagdB,[],2)).');
        hRDMplot.EdgeColor = 'none';
        hRDMplot.FaceColor = 'interp';
        ylabel('Speed (m/s)')
        xlabel('Range (km)')
        title('Range-Doppler')
        hC = colorbar;
        hC.Label.String = 'Magnitude (dB)';
    else
        % Update range-angle plot
        hRAMplot.CData = angRespMagdB;

        % Update range-Doppler plot
        hRDMplot.CData = squeeze(max(respMagdB,[],2)).';
    end
    xlim(hRAMaxes,[max((tgtSR(m) - 2e3),rngVec(1)) min((tgtSR(m) + 2e3),rngVec(end))]*1e-3)
    xlim(hRDMaxes,[max((tgtSR(m) - 2e3),rngVec(1)) min((tgtSR(m) + 2e3),rngVec(end))]*1e-3)
    
    drawnow limitrate
    pause(0.25)

Finally, detect and estimate the target range and angle using the helper function helperDetectAndEstimate. This helper function contains the broadside beamforming, CFAR detection, and angle and range estimation processing steps.

    % Detect targets and estimate range and angle
    dets = helperDetectAndEstimate(cfar,sumdiff,steeringvec,rangeestimator, ...
        m,rngVec,dopVec,resp,dets);
end

Assess the Refraction Errors

Predict the range and elevation that will be output from the detector including refraction effects. For the calculation of the predicted range and elevation angle, use the slant2range function with the CRPL method selected.

% Estimate propagated range using CRPL
[predicted.Range,predicted.Elevation] = slant2range( ...
    tgtTruth.Range,sensorHgt,tgtTruth.Height, ...
    'Method','CRPL');

Now calculate the range error. The range error is the difference between the propagated range and the true slant range. The detected errors in blue follow the predicted range errors, shown as the dashed line in the figure. Note that the IQ processing introduced an additional source of bias. This bias can be attributed to factors like:

• The range resolution of the system, which is fairly large at about 10 m and

• Range bin straddling.

% Calculate the range error (m)
co = colororder;
figure('Name','Range Error (m)')
plot(tgtSR*1e-3,dets.Range - tgtTruth.Range,'o','MarkerFaceColor',co(1,:),'LineWidth',1.5)
hold on
grid on
plot(tgtSR*1e-3,predicted.Range - tgtTruth.Range,'--k','LineWidth',1.5)
xlabel('Slant Range (km)')
ylabel('Range Error (m)')
title('Range Error (m)')
legend('Detected Range Error','Predicted Range Error','Location','Best')

Next calculate the elevation error. The elevation error is the difference between the detected elevation angle and the true elevation angle. Note that the detected errors match the predicted elevation errors well. As the target moves farther out in range, the elevation errors increase.

% Calculate the elevation angle error (deg)
figure('Name','Elevation Error (deg)')
plot(tgtSR*1e-3,dets.Angles(2,:) - tgtTruth.Elevation,'o','MarkerFaceColor',co(1,:),'LineWidth',1.5)
hold on
grid on
plot(tgtSR*1e-3,predicted.Elevation - tgtTruth.Elevation,'--k','LineWidth',1.5)
xlabel('Slant Range (km)')
ylabel('Elevation Error (deg)')
title('Elevation Error (deg)')
legend('Detected Elevation Error','Predicted Elevation Error','Location','Best')

Next, calculate the height error. The height error is the difference between the apparent height, which is the height calculated assuming no refraction, and the true target height. To calculate the apparent height, use the range2height function with the true Earth radius. Note that the detected height errors match the predicted height errors well. As the target moves farther out in range, the height errors slowly increase. At the farthest range, height errors of about 20 m can be seen.

% Calculate apparent height from detections
dets.Height = range2height(dets.Range,sensorHgt,dets.Angles(2,:), ...
    'EffectiveEarthRadius',Rearth);

% Calculate apparent height from predicted path
predicted.Height = range2height(predicted.Range,sensorHgt,predicted.Elevation, ...
    'EffectiveEarthRadius',Rearth);

% Plot height errors (m)
figure('Name','Height Error (m)')
plot(tgtSR*1e-3,dets.Height - tgtTruth.Height,'o','MarkerFaceColor',co(1,:),'LineWidth',1.5)
hold on
grid on
plot(tgtSR*1e-3,predicted.Height - tgtTruth.Height,'--k','LineWidth',1.5)
xlabel('Slant Range (km)')
ylabel('Height Error (m)')
legend('Detected Height Error','Predicted Height Error','Location','Best')

Summary

This example demonstrated how to model refraction for an L-band surveillance radar operating at a high latitude during winter. It also demonstrated how to configure and implement a signal processing chain to detect the target. Lastly, it discussed the target detection errors due to refraction.

Helpers

helperDetectAndEstimate

function dets = helperDetectAndEstimate(cfar,sumdiff,steeringvec,rangeestimator, ...
    tgtIdx,rngVec,dopVec,iqPCDP,dets)
% Detection processing. Performs the following: 
% - Broadside beamforming
% - 2D CFAR
% - Angle estimation
% - Range estimation

% Number of bins
numRngBins = numel(rngVec);
numDopBins = numel(dopVec);

% Convert datacube to power
freq = sumdiff.OperatingFrequency; % Operating frequency (Hz)
sv = steeringvec(freq,[0;0]); % Broadside
iqPCDPtmp = permute(iqPCDP,[2 1 3]); % <elements x range x Doppler>
iqPCDPBF = (sv'*reshape(iqPCDPtmp,[],numRngBins*numDopBins)); % Beamform <beams x range x Doppler>
iqPow = abs(squeeze(iqPCDPBF)).^2; % Power (Watts)
iqPow = reshape(iqPow,numel(rngVec),numel(dopVec)); % Power <range x Doppler>

% Run 2D CFAR. The indices of the cells under test for the CFAR must
% permit the stencil size.
numRngStencil = cfar.GuardBandSize(1) + cfar.TrainingBandSize(1);
numDopStencil = cfar.GuardBandSize(2) + cfar.TrainingBandSize(2);
idxRngTest = (numRngStencil + 1):(numel(rngVec) - numRngStencil);
idxDopTest = (numDopStencil + 1):(numel(dopVec) - numDopStencil);
[idxRngGrid,idxDopGrid] = meshgrid(idxRngTest,idxDopTest);
[idxDet,noisePwr] = cfar(iqPow,[idxRngGrid(:).'; idxDopGrid(:).']);

% Assume maximum detection SNR is target
idxDetLinear = sub2ind(size(iqPow),idxDet(1,:),idxDet(2,:));
snrDet = iqPow(idxDetLinear)./noisePwr;
[~,idxMax] = max(snrDet);
idxRng = idxDet(1,idxMax);
idxDop = idxDet(2,idxMax);

% Estimate true angle. Assume sensor is aligned with target.
dets.Angles(:,tgtIdx) = sumdiff(squeeze(iqPCDP(idxRng,:,idxDop)), ...
    [0; 0]);
dets.Angles(2,tgtIdx) = -dets.Angles(2,tgtIdx); 

% Beamform to detected angle
sv = steeringvec(freq,dets.Angles(:,tgtIdx));
iqPCDPBF = (sv'*squeeze(iqPCDP(:,:,idxDop)).');

% Estimate SNR and range
dets.SNRdB(tgtIdx) = pow2db(abs(squeeze(iqPCDPBF(idxRng))).^2./noisePwr(idxMax));
dets.Range(tgtIdx) = rangeestimator(iqPCDPBF(:),rngVec,idxRng);
end

×

Abrir ejemplo

Tiene una versión modificada de este ejemplo. ¿Desea abrir este ejemplo con sus modificaciones?

No, sobrescribir la versión modificada Sí