Simulating Polarimetric Radar Returns for Weather Observations

Since R2021a

This example shows how to simulate polarimetric Doppler radar returns that meet the requirements of weather observations. Radar plays a critical role in weather observation, detection of hazards, classification and quantification of precipitation, and forecasting. In addition, polarimetric radar provides multiparameter measurements with unprecedented quality and information. This example shows how to simulate a polarimetric Doppler radar that scans an area of distributed weather targets. The simulation derives the radar parameters according to the well-known NEXRAD radar specifications. After collecting the received pulses using radarTransceiver, radar spectral moment and polarimetric variable estimation are performed. The estimates are compared with NEXRAD ground truth, from which error statistics are obtained and data quality is evaluated.

Radar Definition

A well-known weather radar is the Weather Surveillance Radar, 1988 Doppler (WSR-88D), also known as NEXRAD, which is operated by the US National Weather Service, FAA and DoD. For more information, see the NEXRAD Radar Operations Center website.

Radar system specifications are designed as follows.

Radar System Specifications

rng('default');                        % Set the random seed for repeatability
maxRange = 150e3;                      % Maximum unambiguous range (m)
rangeRes = 250;                        % Required range resolution (m)
numPulses = 64;                        % Number of pulses in each dwell time
fc = 3e9;                              % Frequency (Hz)   
propSpeed = 3e8;                       % Propagation speed (m/s)
lambda = propSpeed/fc;                 % Wavelength (m)

% Waveform
pulseBw = propSpeed/(2*rangeRes);      % Pulse bandwidth (Hz)
prf = propSpeed/(2*maxRange);          % Pulse repetition frequency (Hz)
fs = pulseBw;                          % Sampling rate (Hz)
pulseWidth = 1/pulseBw;                % Pulse width (s)
waveform = phased.RectangularWaveform(...
    'PulseWidth',pulseWidth,'PRF',prf,'SampleRate',fs);

% Transmitter                           
txGain = 53;                           % Transmitter gain (dB)
peakPower = 700e3;                     % Transmit peak power (W) 
transmitter = phased.Transmitter(...
    'Gain',txGain,'PeakPower',peakPower);

% Receiver
NF = 2.7;                              % Noise figure (dB)                  
rxGain = 30;                           % Receiver gain (dB)
receiver = phased.ReceiverPreamp(...
    'Gain',rxGain,'NoiseFigure',NF,'SampleRate',fs);

% Matched filter
matchedfilter = phased.MatchedFilter(...
    'Coefficients',getMatchedFilter(waveform));
matchingdelay = size(matchedfilter.Coefficients,1)-1;

Antenna Pattern

As NEXRAD is a polarimetric radar, modeling the polarimetric characteristics of the antenna is important. According to NEXRAD specifications, the antenna pattern has a beamwidth of about 1 degree and first sidelobe below -30 dB.

azang = -180:0.5:180;
elang = -90:0.5:90;
% We synthesize a pattern using isotropic antenna elements and tapering the
% amplitude distribution to make it follow NEXRAD specifications.
magpattern = load('NEXRAD_pattern.mat');
phasepattern = zeros(size(magpattern.pat));
% The polarimetric antenna is assumed to have ideally matched horizontal
% and vertical polarization pattern.
hPolAntenna = phased.CustomAntennaElement(...
    'SpecifyPolarizationPattern',true,'AzimuthAngles',azang,'ElevationAngles',elang,...
    'HorizontalMagnitudePattern',magpattern.pat,'VerticalMagnitudePattern',-inf(size(magpattern.pat)), ...
    'HorizontalPhasePattern',phasepattern,'VerticalPhasePattern',phasepattern);
vPolAntenna = phased.CustomAntennaElement(...
    'SpecifyPolarizationPattern',true,'AzimuthAngles',azang,'ElevationAngles',elang,...
    'HorizontalMagnitudePattern',-inf(size(magpattern.pat)), 'VerticalMagnitudePattern',magpattern.pat,...
    'HorizontalPhasePattern',phasepattern,'VerticalPhasePattern',phasepattern);

Plot the azimuth cut of the antenna pattern in horizontal polarization.

D = pattern(hPolAntenna,fc,azang,0);
P = polarpattern(azang,D,'TitleTop','Polar Pattern for Azimuth Cut (elevation angle = 0 degree)');
P.AntennaMetrics = 1;
removeAllCursors(P);

Associate the array with the radiator and collector for each polarization.

radiatorH = phased.Radiator(...
    'Sensor',hPolAntenna,'Polarization','Combined',...
    'OperatingFrequency',fc);
collectorH = phased.Collector(...
    'Sensor',hPolAntenna,'Polarization','Combined',...
    'OperatingFrequency',fc);
radiatorV = phased.Radiator(...
    'Sensor',vPolAntenna,'Polarization','Combined',...
    'OperatingFrequency',fc);
collectorV = phased.Collector(...
    'Sensor',vPolAntenna,'Polarization','Combined',...
    'OperatingFrequency',fc);

To simulate a polarimetric radarTransceiver, define radarTransceiver for each polarization.

xcvrHH = radarTransceiver('Waveform',waveform, ...
    'Transmitter',transmitter,'TransmitAntenna',radiatorH, ...
    'Receiver',receiver,'ReceiveAntenna',collectorH,...
    'MechanicalScanMode','Circular','MechanicalScanRate',0);
xcvrVV = radarTransceiver('Waveform',waveform, ...
    'Transmitter',transmitter,'TransmitAntenna',radiatorV, ...
    'Receiver',receiver,'ReceiveAntenna',collectorV,...
    'MechanicalScanMode','Circular','MechanicalScanRate',0); 

Weather Target

Generally, weather radar data is categorized into three levels. Level-I data is raw time series I/Q data as input to the signal processor in the Radar Data Acquisition unit. Level-II data consists of the radar spectral moments (reflectivity, radial velocity, and spectrum width) and polarimetric variables (differential reflectivity, correlation coefficient, and differential phase) output from the signal processor. Level-III data is the output product data of the radar product generator, such as hydrometeor classification, storm total precipitation, and tornadic vortex signature.

In this example, Level-II data from KTLX NEXRAD radar at 04:33:45 UTC on May 19th, 2013 is used. This data comes from a convective precipitation that occurred in Oklahoma and is used to generate scattering matrix of weather targets. The data is available via FTP download. It represents a volume scan that includes a series of 360-degree sweeps of the antenna at predetermined elevation angles completed in a specified period of time. In the data file name KTLX20130519_043345_V06, KTLX refers to the radar site name, 20130519_043345 refers to the date and time when the data was collected, and V06 refers to the data format of version 6. In this simulation, the lowest elevation cut (0.5 degree) is extracted from the volume scan data as an example.

Read the Level-II data into the workspace. Store it in the nexrad structure array, which contains all the radar moments as well as an azimuth field that specifies the azimuth angle for each radial data point in the Cartesian coordinate system. For simplicity, load NEXRAD data that was transformed from a compressed file to a MAT-file and round azimuth to the nearest integer.

load NEXRAD_data.mat;
nexrad.azimuth = round(nexrad.azimuth);

Define an area of interest (AOI) in terms of azimuth and range in Cartesian coordinates.

az1 = 42;           % Starting azimuth angle (degree) 
az2 = 48;           % Ending azimuth angle (degree) 
rg1 = 20000;        % Starting range (m) 
rg2 = 30000;        % Ending range (m)  
blindRng = 2000;    % NEXRAD's blind range
b0 = blindRng/rangeRes;         % Number of range bins within the blind range
numAz = az2-az1+1;              % Number of azimuth angles in the AOI
numBins = (rg2-rg1)/rangeRes+1; % Number of range bins in each azimuth
% Select AOI data and store it in _nexrad_aoi_ structure array, which
% contains all the radar moments, as well as starting and ending azimuth
% and range indices. 
nexrad_aoi = helperSelectAOI(nexrad,az1,az2,rg1,rg2,blindRng,rangeRes);

Weather Measurements Collection

To model weather target comprised of distributed volume scatterers, each radar resolution volume is regarded as an equivalent scattering center whose scattering matrix is updated between pulse to pulse. In this example, the Monte Carlo method is used to simulate the randomly scattered wave field and generate scattering matrix from the NEXRAD Level-II data.

% Preallocate datacube
rxh_aoi = complex(zeros(numBins,numPulses));    % H-pol 2D datacube
rxv_aoi = complex(zeros(numBins,numPulses));    % V-pol 2D datacube
iqhh = complex(zeros(numAz,numBins,numPulses)); % H-pol 3D datacube
iqvv = complex(zeros(numAz,numBins,numPulses)); % V-pol 3D datacube
rdrpos = [0;0;0];   % Radar position (m)
el = 0;             % Elevation of target plane relative to radar
simTime = 0;        % Simulation time  
sceneAxes = eye(3); % Scenario axes
% Mechanical scanning
for ii = 1:numAz
    az = nexrad.azimuth(nexrad_aoi.r1+ii-1);   % current azimuth
    axes = azelaxes(wrapTo180(az),el);         % current axes
    % Preallocate target position and target poses
    tgtpos = zeros(3,numBins);      
    tgtPoses = repmat(struct('Position',zeros(3,1),...
        'Signatures',{}),1,numBins); 
    for jj = 1:numBins
        rpos = (nexrad_aoi.b1+jj-2)*rangeRes + blindRng;             
        tpos = [rpos*cosd(el)*cosd(az);rpos*cosd(el)*sind(az);rpos*sind(el)];
        tgtpos(:,jj) = tpos; 
    end
    % Calculate the target range
    [tgtrng,~] = rangeangle(tgtpos,rdrpos,axes); 
    % Initial mechanical scan angle for each dwell time
    release(xcvrHH);
    xcvrHH.InitialMechanicalScanAngle = az; 
    release(xcvrVV);
    xcvrVV.InitialMechanicalScanAngle = az; 
    for mm = 1:numPulses  
        randStream = RandStream('twister','Seed',0);
        for jj = 1:numBins
            % Update scattering matrix of weather target every PRT
            [shhpat,svvpat] = helperScatteringMatrix(lambda,nexrad_aoi.ZH(ii,jj),nexrad_aoi.Vr(ii,jj),...
                nexrad_aoi.SW(ii,jj),nexrad_aoi.ZDR(ii,jj),nexrad_aoi.Rhohv(ii,jj),nexrad_aoi.Phidp(ii,jj),simTime,randStream);
            tgtPoses(jj).Signatures = rcsSignature('EnablePolarization',true, ...
                'ShhPattern',shhpat*tgtrng(jj),...
                'SvvPattern',svvpat*tgtrng(jj));
            tgtPoses(jj).Position = tgtpos(:,jj); 
        end
        % Step H-pol radarTransceiver
        rxh = xcvrHH(tgtPoses,simTime,sceneAxes); 
        % Step V-pol radarTransceiver
        rxv = xcvrVV(tgtPoses,simTime,sceneAxes); 
        % Matched filtering
        rxh = matchedfilter(rxh); 
        rxv = matchedfilter(rxv);
        rxh = [rxh(matchingdelay+1:end);zeros(matchingdelay,1)];
        rxv = [rxv(matchingdelay+1:end);zeros(matchingdelay,1)];
        % Discard blind range data and select AOI data
        rxh_aoi(:,mm) = rxh(nexrad_aoi.b1+b0:nexrad_aoi.b2+b0);
        rxv_aoi(:,mm) = rxv(nexrad_aoi.b1+b0:nexrad_aoi.b2+b0);
        % Update simulation time
        simTime = simTime + 1/prf; 
    end
    iqhh(ii,:,:) = rxh_aoi;
    iqvv(ii,:,:) = rxv_aoi;
end

Weather Radar Moment Estimation

Using pulse pair processing, calculate all the radar moments from estimates of correlations, including reflectivity, radial velocity, spectrum width, differential reflectivity, correlation coefficient, and differential phase.

azimuth = nexrad.azimuth(nexrad_aoi.r1:nexrad_aoi.r2);
range = (nexrad_aoi.b1-1:nexrad_aoi.b2-1)*rangeRes + blindRng; 
for ii = 1:numAz
    for jj = 1:numBins
        moment0 = helperPulsePairProcessing(squeeze(iqhh(ii,jj,:)),squeeze(iqvv(ii,jj,:)),lambda,prf,range(jj));
        moment.ZH(ii,jj) = moment0.ZH;
        moment.Vr(ii,jj) = moment0.Vr;
        moment.SW(ii,jj) = moment0.SW;
        moment.ZDR(ii,jj) = moment0.ZDR;
        moment.Rhohv(ii,jj) = moment0.Rhohv;
        moment.Phidp(ii,jj) = moment0.Phidp;    
    end
end

Simulation Result

Compare the simulation result with the NEXRAD ground truth. Evaluate the simulated data quality using error statistics, a sector image, a range profile, and a scatter plot. Error statistics are expressed as the bias and standard deviation of the estimated radar moments compared to the NEXRAD Level-II data (truth fields).

Reflectivity

Reflectivity, $$Z$$, is the zeroth moment of the Doppler spectrum and is related to liquid water content or precipitation rate in the resolution volume. Because values of $$Z$$ that are commonly encountered in weather observations span many orders of magnitude, radar meteorologists use a logarithmic scale given by $$10lo{g}_{10}Z$$ as dBZ, where $$Z$$ is in units of mm^6/m^3.

[ZHBias,ZHStd] = helperDataQuality(nexrad_aoi,moment,range,azimuth,'ZH');

Radial Velocity

Radial velocity, $$V_r$$, is the first moment of the power-normalized spectra, which reflects the air motion toward or away from the radar.

[VrBias,VrStd] = helperDataQuality(nexrad_aoi,moment,range,azimuth,'Vr');

Spectrum Width

Spectrum width, $${\sigma }_v$$, is the square root of the second moment of the normalized spectrum. The spectrum width is a measure of the velocity dispersion, that is, shear or turbulence within the resolution volume.

[SWBias,SWStd] = helperDataQuality(nexrad_aoi,moment,range,azimuth,'SW');

Differential Reflectivity

Differential reflectivity, $$Z_{DR}$$, is estimated from the ratio of the returned power for the horizontal and vertical polarization signals. The differential reflectivity is useful in hydrometeor classification.

[ZDRBias,ZDRStd] = helperDataQuality(nexrad_aoi,moment,range,azimuth,'ZDR');

Correlation Coefficient

The correlation coefficient, $${\rho }_{hv}$$, represents the consistency of the horizontal and vertical returned power and phase for each pulse. The correlation coefficient plays an important role in determining polarimetric radar system performance and classifying radar echo types.

[RhohvBias,RhohvStd] = helperDataQuality(nexrad_aoi,moment,range,azimuth,'Rhohv');

Differential Phase

The differential phase, $${\varphi }_{DP}$$, is the difference in the phase delay of the returned pulse between the horizontal and vertical polarizations. The differential phase provides information on the nature of the scatterers that are being sampled.

[PhidpBias,PhidpStd] = helperDataQuality(nexrad_aoi,moment,range,azimuth,'Phidp');

Error Statistics

Figures in previous sections provide a visual qualitative measure of the simulation quality. This section of the example shows the quantitative comparison of the estimates with NEXRAD ground truth as error statistics including bias and standard deviation (STD).

MomentName = {'ZH';'Vr';'SW';'ZDR';'Rhohv';'Phidp'};
BiasEstimate = [round(mean(ZHBias),2);round(mean(VrBias),2);round(mean(SWBias),2);...
    round(mean(ZDRBias),2);round(mean(RhohvBias),3);round(mean(PhidpBias),2)];
StdEstimate = [round(mean(ZHStd),2);round(mean(VrStd),2);round(mean(SWStd),2);...
    round(mean(ZDRStd),2);round(mean(RhohvStd),3);round(mean(PhidpStd),2)];
StdTruth = helperTruthSTD(nexrad_aoi.SW,nexrad_aoi.Rhohv,numPulses,1/prf,lambda);
StdTheory = [round(StdTruth.ZH,2);round(StdTruth.Vr,2);round(StdTruth.SW,2);...
    round(StdTruth.ZDR,2);round(StdTruth.Rhohv,3);round(StdTruth.Phidp,2)];
Unit = {'dB';'m/s';'m/s';'dB';'Unitless';'degree'};
T = table(MomentName,BiasEstimate,StdEstimate,StdTheory,Unit);
disp(T);

    MomentName    BiasEstimate    StdEstimate    StdTheory        Unit    
    __________    ____________    ___________    _________    ____________

    {'ZH'   }            0            1.37          1.25      {'dB'      }
    {'Vr'   }        -0.09            0.57          0.48      {'m/s'     }
    {'SW'   }        -0.06            0.56          0.47      {'m/s'     }
    {'ZDR'  }         0.03            0.36          0.32      {'dB'      }
    {'Rhohv'}        0.001           0.008         0.006      {'Unitless'}
    {'Phidp'}         0.03            2.15          2.11      {'degree'  }

By comparison, the mean bias for all the radar estimates is close to zero, and the STD of all the radar estimates is consistent with the STD of theory. It is expected that the STD of estimates is a little bit larger than the STD of theory because of sampling error. If numPulses is increased, the STD of estimates will be closer to the STD of theory, while it will take a longer time to run the simulation.

Summary

This example showed how to simulate the polarimetric Doppler radar returns from an area of distributed weather targets. Visual comparison and error statistics showed the estimated radar moments are consistent with NEXRAD ground truth. With this example, you can further explore the simulated time series data in other applications such as waveform design, clutter filtering, and data quality evaluation for weather radar.

References

[1] Doviak, R and D. Zrnic. Doppler Radar and Weather Observations, 2nd Ed. New York: Dover, 2006.

[2] Zhang, G. Weather Radar Polarimetry. Boca Raton: CRC Press, 2016.

[3] Li, Z, et al. "Polarimetric phased array weather radar data quality evaluation through combined analysis, simulation, and measurements," IEEE Geosci. Remote Sens. Lett., vol. 18, no. 6, pp. 1029–1033, Jun. 2021.