Target Localization in Active and Passive Radar

Localizing a target based on measurements from an array of spatially separated sensors with known position has been an important problem in radar for a long time. Time of arrival (TOA) and time difference of arrival (TDOA) are commonly used measurements for localization. This example shows how to model radar networks, configure and propagate radar waveforms, and perform TOA/TDOA estimation and localization.

Introduction

Localizing a target using radars can be realized in multiple types of radar systems. For example, the most common system is a monostatic active radar system that localizes a target by actively transmitting radar waveforms and receiving the target backscattered signals using co-located and synchronized transmitter and receiver. A bistatic passive radar system can also be used to localize a target by passively receiving target backscattered signals from a separate transmit source using synchronized and distributed radars.

TOA and TDOA are commonly used measurements that provide distance information for radar localization. For a monostatic active radar system, TOA measurements are typically adopted and estimated from round-trip propagation delays between a target and monostatic radar transceivers. For a bistatic passive radar system where the separate transmit source is asynchronous to the synchronized radar receivers, the TOA measurements at radar receivers are shifted by a common unknown delay. In this case, more suitable measurements are TDOA measurements that provide the differences between the TOA measurements such that the common unknown delay is canceled. In this example, we are going to show the applications of TOA and TDOA estimation and localization in their suitable radar systems.

TOA Estimation and Localization Using Monostatic Active Radars

In this section, we consider TOA localization of a target using a monostatic active radar system, where the radars do not need to be mutually synchronized but the transmitter and the receiver at each radar are synchronized.

Scenario Configuration

We consider using 5 millimeter wave (mmWave) radars built at road infrastructures with known locations to localize a vehicle target.

rng('default')

% RF parameters
fc = 77e9;                                      % Carrier frequency (Hz)
c = physconst('LightSpeed');                    % Speed of propagation (m/s)
bw = 200e6;                                     % Bandwidth (Hz)
sampleRate = bw;                                % Sample rate (Hz)

% Tx and Rx parameters
Pt = 1;                                         % Peak power (W) 
Gtx = 20;                                       % Tx antenna gain (dB)
Grx = 40;                                       % Rx antenna gain (dB) 
NF = 2.9;                                       % Noise figure (dB) 

% Create a sensor and sensor components to simulate waveform propagation
antenna = phased.IsotropicAntennaElement('BackBaffled',false);
transmitter = phased.Transmitter('Gain',Gtx,'PeakPower',Pt);
radiator = phased.Radiator('Sensor',antenna,'OperatingFrequency',fc);
collector = phased.Collector('Sensor',antenna,'OperatingFrequency',fc);
receiver = phased.Receiver('AddInputNoise',true,'Gain',Grx,...
         'NoiseFigure',NF,'SampleRate',sampleRate);

% Create a target
tgtrcs = 10;                                    % Target RCS (m^2)
target = phased.RadarTarget('MeanRCS',tgtrcs,'PropagationSpeed',c,...
    'OperatingFrequency',fc);

% Target platform
tgtpos = [30; 10; 5];                           % Target position (m)
tgtvel = [5; 10; 0];                            % Target velocity (m/s)
tgtplatform = phased.Platform('InitialPosition',tgtpos,'Velocity',tgtvel); 

% Radar platform
radarpos = [0 -30 100 80 -40;...
    0 50 -40 30 -20;...
    0 -5 10 5 20];                              % Radar positions (m)
numRadar = size(radarpos,2);                    % Number of radars 
radarvel = zeros(3,numRadar);                   % Radar velocities (m/s)
radarplatform = phased.Platform('InitialPosition',radarpos,'Velocity',radarvel);

Waveform Configuration

One of the most popular waveforms for radar systems is the frequency-modulated continuous waveform or FMCW Waveforms. FMCW waveform is widely adopted in automotive radar systems as this technique is mature and low-cost.

% Configure FMCW waveform
N = 1024;                                        % Number of subbands (fast-time samples)
M = 8;                                           % Number of slow-time pulses
freqSpacing = bw/N;                              % Frequency spacing (Hz)
tWave = N/sampleRate;                            % Waveform duration (s)
fmcwWaveform = phased.FMCWWaveform('SweepTime',tWave,'SweepBandwidth',bw,...
    'SampleRate',sampleRate,'SweepDirection','Up');
sig = fmcwWaveform()*ones(1,M);

Channel Configuration

For a monostatic active radar system, the propagation channel between each radar transceiver and the target is a two-way channel. We do not consider mutual radar interference in this example.

% Create a two-way free-space propagation channel
channel = phased.FreeSpace('PropagationSpeed',c,'OperatingFrequency',fc,...
    'SampleRate',sampleRate,'TwoWayPropagation',true);

Radar Data Cube

We first obtain the Radar Data Cube for each monostatic active radar. In this example, a radar data cube is defined as the complex-valued baseband samples arranged in a N-by-M matrix, where N is the number of fast-time samples and M is the number of slow-time pulses.

% Reference signal
refsig = sig;

% Estimated channel
X = cell(1,numRadar);

% Transmitted signal
txsig = transmitter(sig);

for idxRadar = 1:numRadar
    % Initialize estimated channel
    x = complex(zeros(size(sig)));

    % Transceiver chain
    for m = 1:M
        % Update radar and target positions
        [radar_pos,radar_vel] = radarplatform(tWave);
        [tgt_pos,tgt_vel] = tgtplatform(tWave);

        % Calculate the transmit angle
        [~,txang] = rangeangle(tgt_pos,radar_pos(:,idxRadar)); 

        % Radiate signal towards the receiver 
        radtxsig = radiator(txsig(:,m),txang); 

        % Propagate the signal
        chansig = channel(radtxsig,radar_pos(:,idxRadar),tgt_pos,...
            radar_vel(:,idxRadar),tgt_vel);

        % Reflect the signal off the target
        tgtsig = target(chansig);

        % Calculate the receive angle
        [~,rxang] = rangeangle(radar_pos(:,idxRadar),tgt_pos);

        % Collect signal at the receive antenna
        rxsig = collector(tgtsig,rxang);

        % Receive signal at the receiver
        x(:,m) = receiver(rxsig);
    end

    % Channel estimation via dechirping
    dechirpsig = dechirp(x,refsig(:,1));
    X{idxRadar} = conj(dechirpsig);

    % Reset platform for the next radar
    reset(radarplatform);
    reset(tgtplatform);
end

TOA Estimation and Localization

After obtaining the radar data cubes, the next step is to obtain the TOA measurements. In the code, we use phased.TOAEstimator system object to estimate TOA by configuring the 'Measurement' property to 'TOA'. The spectrum analysis method can be configured as either FFT or MUSIC. We can change the spectrum analysis method to see its impact on the TOA localization accuracy shown at the end of this section.

% Spectrum analysis method
spectrumMethod = "FFT";
 
% Configure TOA estimator
if strcmp(spectrumMethod,'FFT') 
    toaEstimator = phased.TOAEstimator('PropagationSpeed',c,...
        'Measurement','TOA','SpectrumMethod',spectrumMethod,... 
        'VarianceOutputPort',true);
else % 'MUSIC'
    toaEstimator = phased.TOAEstimator('PropagationSpeed',c,...
        'Measurement','TOA','SpectrumMethod',spectrumMethod,... 
        'VarianceOutputPort',true,...
        'ForwardBackwardAveraging',true,'SpatialSmoothing',ceil(N/2)); %#ok<UNRCH>
end

% TOA estimation
[Y,var] = toaEstimator(X,freqSpacing);

Viewing the TOA spectrum can provide us better insights on how well TOA is estimated for each radar. In the ideal case, the first peak on the TOA spectrum is picked as the TOA estimate. We can view the TOA estimation result for a radar using plotTOASpectrum.

% Plot TOA spectrum
figure
plotTOASpectrum(toaEstimator,freqSpacing,'AnchorIndex',1,'MaxDelay',500e-9);

Once we have the one-way propagation TOA estimates, the one-way propagation TOA estimation variances, and known radar positions, we can use toaposest to localize the target. As the obtained TOA estimates are the round-trip propagation delay estimates, we divide the round-trip propagation delay estimates by 2 to obtain one-way propagation TOA estimates and divide the estimated TOA variances by 4 to obtain estimated one-way propagation TOA estimation variances.

% Single-trip TOA and TOA variance
Ysingle = Y/2;
varsingle = var/4;

% Obtain TOA position estimate
tgtposest = toaposest(Ysingle,varsingle,radarpos);

We visualize the performance of the TOA localization below.

% View TOA position estimate
helperPlotTOAPositions(c,Ysingle,tgtposest,radarpos,tgtpos);

We can also check the accuracy of the TOA localization using root mean square error (RMSE).

% RMSE of the TOA position estimate
RMSE = rmse(tgtposest,tgtpos);
disp(['RMS Localization error = ', num2str(RMSE), ' meters.'])

RMS Localization error = 0.68646 meters.

TDOA Estimation and Localization Using Bistatic Passive Radars

In this section, we consider TOA-based TDOA localization of a target using a bistatic passive radar system, where the distributed radar receivers are mutually synchronized and passively receive backscattered signals from a separate asynchronous transmit source with known waveform. For a passive radar system, the separate transmit source can be another radar station, a cellular base station, a TV transmitter tower, etc.

Scenario Configuration

We consider using 5 ground radar receivers with known locations to localize an aircraft target using backscattered signals from a separate asynchronous radar transmitter tower operating at the C-band.

% RF parameters
fc = 6e9;                                       % Carrier frequency (Hz)
c = physconst('LightSpeed');                    % Speed of propagation (m/s)
bw = 100e6;                                     % Bandwidth (Hz)
sampleRate = bw;                                % Sample rate (Hz)

% Tx and Rx parameters
Pt = 10;                                        % Peak power (W) 
Gtx = 40;                                       % Tx antenna gain (dB)
Grx = 40;                                       % Rx antenna gain (dB) 
NF = 2.9;                                       % Noise figure (dB) 

% Create a sensor and sensor components to simulate waveform propagation
antenna = phased.IsotropicAntennaElement('BackBaffled',false);
transmitter = phased.Transmitter('Gain',Gtx,'PeakPower',Pt);
radiator = phased.Radiator('Sensor',antenna,'OperatingFrequency',fc);
collector = phased.Collector('Sensor',antenna,'OperatingFrequency',fc);
receiver = phased.Receiver('AddInputNoise',true,'Gain',Grx,...
         'NoiseFigure',NF,'SampleRate',sampleRate);

% Create a target
tgtrcs = 1000;                                  % Target RCS (m^2)
target = phased.RadarTarget('MeanRCS',tgtrcs,'PropagationSpeed',c,...
    'OperatingFrequency',fc);

% Target platform
tgtpos = [80; 40; 110];                         % Target position (m)
tgtvel = [50; 40; 0];                           % Target velocity (m/s)
tgtplatform = phased.Platform('InitialPosition',tgtpos,'Velocity',tgtvel); 

% Radar platform
radarpos = [0 300 100 200 150;...
    0 50 -200 300 100;...
    0 -10 10 5 20];                             % Radar positions (m)
numRadar = size(radarpos,2);                    % Number of radars 
radarvel = zeros(3,numRadar);                   % Radar velocities (m/s)
radarplatform = phased.Platform('InitialPosition',radarpos,'Velocity',radarvel);

% Separate transmitter platform
txpos = [-150; -100; 50];                       % Transmitter position (m)
txvel = [0; 0; 0];                              % Transmitter velocity (m/s)
txplatform = phased.Platform('InitialPosition',txpos,'Velocity',txvel); 

Waveform Configuration

Phase-Coded Waveforms is another type of radar waveform. Multiple phase codes can be chosen in the phase coded waveform, such as Zadoff-Chu code, Barker code, maximum length sequence (MLS), etc. In the following, we configure a phase-code waveform with MLS as its phase code.

% Configure phase-coded waveform
N = 1024;                                        % Number of subbands
M = 8;                                           % Number of slow-time pulses
freqSpacing = bw/N;                              % Frequency spacing (Hz)
numChip = 2^nextpow2(N)-1;                       % Number of chips in one phase-coded waveform
tchip = 1/bw;                                    % Chip duration (s)
tWave = numChip * tchip;                         % Modulation period (s)
prf = 1/tWave;                                   % PRF (1/s)

% Configure the phase coded waveform as the maximum length sequence
pmcwWaveform = phased.PhaseCodedWaveform('Code','Maximum Length Sequence',...
    'SampleRate',sampleRate,'NumChips',numChip,'ChipWidth',tchip,'PRF',prf);
sig = pmcwWaveform()*ones(1,M);

Channel Configuration

A bistatic radar system basically contains two one-way channels. The first channel is the channel from the separate transmit source to the target, and the second channel is the channel from the target to each radar. We specify the two channels below.

% One-way free-space propagation channel from the radar transmitter tower to the target
txchannel = phased.FreeSpace('PropagationSpeed',c,'OperatingFrequency',fc,...
    'SampleRate',sampleRate,'TwoWayPropagation',false);

% One-way free-space propagation channel from the target to each radar receiver
rxchannel = clone(txchannel);

Radar Data Cube

For a bistatic passive radar to obtain a radar data cube, a known reference signal is required. However, as the separate transmit source is asynchronous to the radar, the reference signals at the radars are not the same as local oscillator signal at the transmitter, but the local oscillator signal shifted with a common unknown delay offset.

% Common unknown delay offset in number of samples
delayoffset = 50;

% Reference signal with unknown delay offset for continuous waveform
refsig = circshift(sig,delayoffset);

In the following code, we show how to obtain the bistatic radar data cube when phase-coded waveform is used.

% Estimated channel
X = cell(1,numRadar);

% Transmitted signal
txsig = transmitter(sig);

for idxRadar = 1:numRadar
    % Initialize estimated channel
    x = complex(zeros(size(sig)));

    % Transceiver chain
    for m = 1:M
        % Update separate transmitter, radar and target positions
        [tx_pos,tx_vel] = txplatform(tWave);
        [radar_pos,radar_vel] = radarplatform(tWave);
        [tgt_pos,tgt_vel] = tgtplatform(tWave);

        % Calculate the transmit angle
        [~,txang] = rangeangle(tgt_pos,tx_pos); 

        % Radiate signal towards the receiver 
        radtxsig = radiator(txsig(:,m),txang); 

        % Propagate the signal from the transmitter to the target
        txchansig = txchannel(radtxsig,tx_pos,tgt_pos,...
            tx_vel,tgt_vel);

        % Reflect the signal off the target
        tgtsig = target(txchansig);

        % Propagate the signal from the target to each radar
        rxchansig = rxchannel(tgtsig,radar_pos(:,idxRadar),tgt_pos,...
            radar_vel(:,idxRadar),tgt_vel);

        % Calculate the receive angle
        [~,rxang] = rangeangle(radar_pos(:,idxRadar),tgt_pos);

        % Collect signal at the receive antenna
        rxsig = collector(rxchansig,rxang);

        % Receive signal at the receiver
        x(:,m) = receiver(rxsig);
    end

    % Channel estimation via frequency-domain cross-correlation
    xfreq = fft(x,N);
    xreffreq = fft(refsig,N);
    X{idxRadar} = xfreq.*conj(xreffreq);

    % Reset platform for the next radar
    reset(txplatform);
    reset(radarplatform);
    reset(tgtplatform);
end

TDOA Estimation and Localization

After obtaining the bistatic passive radar data cube, the next step is to obtain the TDOA measurements. In the following, we use phased.TOAEstimator system object to estimate TDOAs by configuring the 'Measurement' property to 'TDOA'.

% Spectrum analysis method
spectrumMethod = "MUSIC";
 
% Configure TOA estimator
if strcmp(spectrumMethod,'FFT') 
    tdoaEstimator = phased.TOAEstimator('PropagationSpeed',c,...
        'Measurement','TDOA','SpectrumMethod',spectrumMethod,... 
        'VarianceOutputPort',true); %#ok<UNRCH>
else % 'MUSIC'
    tdoaEstimator = phased.TOAEstimator('PropagationSpeed',c,...
        'Measurement','TDOA','SpectrumMethod',spectrumMethod,... 
        'VarianceOutputPort',true,...
        'ForwardBackwardAveraging',true,'SpatialSmoothing',ceil(N/2));
end

% TDOA estimation
[Y,var] = tdoaEstimator(X,freqSpacing);

Once we have the TDOA estimates, the TDOA estimation variances, and known radar positions, we can use tdoaposest to do TDOA localization.

% Obtain TDOA position estimate
tgtposest = tdoaposest(Y,var,radarpos);

We visualize the performance of the TDOA localization below.

% View TDOA position estimate
helperPlotTDOAPositions(c,Y,tgtposest,radarpos,tgtpos,txpos)

We can also check the accuracy of the TDOA localization using RMSE.

% RMSE of the TDOA position estimate
RMSE = rmse(tgtposest,tgtpos);
disp(['RMS Localization error = ', num2str(RMSE), ' meters.'])

RMS Localization error = 0.57713 meters.

Summary

In this example, we considered localization of a target in two common radar systems. The first system is a network of monostatic active radars adopting FMCW waveform. We showed how to perform TOA estimation and TOA localization using phased.TOAEstimator and toaposest. The second system is a network of bistatic passive radars adopting phase-coded waveform. We showed how to perform TOA-based TDOA estimation and TDOA localization using phased.TOAEstimator and tdoaposest.

References

[1] Reza Zekavat and R. Michael Buehrer, Handbook of Position Location: Theory, Practice, and Advances, 2019

[2] Hugh Griffiths and Christopher Baker, An Introduction to Passive Radar, Second Edition, 2022