# 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 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 ```

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);```

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}\text{\hspace{0.17em}}}\mathrm{exp}\left(-{\mathit{c}}_{\mathit{e}\text{\hspace{0.17em}}}\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}\text{\hspace{0.17em}}}=\mathrm{ln}\text{\hspace{0.17em}}\frac{{\mathit{N}}_{\mathit{S}\text{\hspace{0.17em}}}}{{\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\mathrm{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 ```

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

```% 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```