Source Localization and Tracking with Passive Receivers

Since R2024b

This example uses:

A passive source localization (PSL) sensor network consists of spatially distributed passive sensors to detect, localize, and track target emitters. The passive sensors are non-transmitting receivers and hence are covert and low-cost. These benefits make PSL sensor networks attractive in many applications, such as air surveillance, acoustic source localization, etc. This example shows how to localize and track targets in a PSL sensor network.

Introduction

The PSL sensor network is different from the passive radar system described in the example Target Localization in Active and Passive Radars. While a passive radar system estimates positions of targets from their scattered signals originated from separate transmitters (like television tower, cellular base stations, navigation satellites, etc.), a PSL sensor network estimates positions of targets by intercepting the targets' emitted signals. In both PSL sensor networks and passive radar systems, the transmitters are non-cooperative to the passive receivers, and the receivers have no knowledge of the source waveforms. We show an illustration figure of a typical PSL sensor network below.

PSL tracking scenario

A common method to localize target emitters is based on time difference of arrival (TDOA) measurements from distinct passive receiver pairs. Using one receiver as the reference receiver, the other receiver can cross-correlate signals at the two receivers to estimate TDOAs between the target emitters and the receiver pair. Geometrically, each noise-free TDOA between a synchronized receiver pair defines a hyperboloid that describes the same target range difference between the two receivers in a three-dimensional (3-D) space. With multiple receiver pairs, the target can be localized at the intersection of different hyperboloids.

In the existence of noise, blockage, multipath, and multiple sources, PSL has the following challenges.

First, the TDOA measurements can be noisy. The noise and interference can make the TDOA measurements deviate from the true TDOA. In a more severe case, some TDOAs can be undetected, and some false detection of TDOAs can appear. Large TDOA measurement deviation, miss detection, and false alarms can significantly reduce the accuracy of target position estimates.
Second, the transmitters are non-cooperative, and hence their emitted signals do not contain signal IDs. When a receiver pair has multiple TDOA measurements due to false alarm and multi-source signals, there is no a priori information about the source of the TDOA measurements, if there is any. This leads to an important problem of data association. That is, assigning a TDOA measurement from a receiver pair to a potential target or to a declared false target; only the TDOA measurements assigned to the same potential target, which form a TDOA group, are used for the target's position estimation. However, data association is not always perfect. Using an incorrectly associated TDOA group can lead to a false target position estimate, which is also known as the ghost target. When the number of receiver pairs and number of TDOA measurements per receiver pair are large, the number of ghost targets can be large.

In this example, we illustrate a workflow to localize target emitters using TDOA measurements in a PSL sensor network. In the existence of multipath, we show that imperfect TDOA detections can lead to missed and ghost targets. Finally, we show a workflow to associate TDOA measurements with reduced number of ghost targets using static detection fuser and improve source position estimates via tracking.

Source Localization with Missed and False Detections

Scenario Configuration

Consider a PSL sensor network with multiple spatially distributed passive receivers, which are also called anchors in wireless localization. The targets are emitters. We configure the number of targets to be 1. You can change the configuration to multiple targets.

% Setup scenario with receiver number and target emitter number
numReceivers = 6;
numTgts = 1;

Consider each of the emitters and receivers are equipped with single isotropic antenna element.

% Antenna elements
srcElement = phased.IsotropicAntennaElement;
rxElement = phased.IsotropicAntennaElement;

Consider the scenario as a 3-D cubic space and specify its boundary.

% Scenario boundary in meters
boundLimit = [-4e2 4e2; ...                     % x_min x_max
    -4e2 4e2; ...                               % y_min y_max
    0 4e2];                                     % z_min z_max

Create a multipath channel with multiple scatterers randomly distributed in the 3-D cubic space using phased.ScatteringMIMOChannel system object. For each emitter-receiver pair, the receiver not only receives line-of-sight signal from the emitter, but also receives the emitter's non-line-of-sight signals reflected from the scatterers.

% Setup scatterer number in multipath channel
numScatterers = 5;

% Get RF parameters
[fc,bw,fs,c] = helperRFParameter;

% Generate scatter positions
scattererPos = helperScatterPosition(numScatterers,boundLimit);

% Configure multipath channel
channel = phased.ScatteringMIMOChannel('TransmitArray',srcElement, ...
    'ReceiveArray',rxElement,'PropagationSpeed',c,'CarrierFrequency',fc, ...
    'SampleRate',fs,'TransmitArrayMotionSource','Input port', ...
    'ReceiveArrayMotionSource','Input port','ScattererSpecificationSource','Property', ...
    'ScattererPosition',scattererPos,'ScattererCoefficient',ones(1,numScatterers), ...
    'SimulateDirectPath',true);

Randomly generate the receiver and target positions in a 3-D cubic space.

% Create scenario
[scenario,rxPairs,receiverIDs] = helperCreateTDOAScenario(numReceivers,numTgts,boundLimit,scattererPos);

Waveform Configuration

Consider emitter signals are generated from phase-coded waveforms. These waveforms are unknown at the passive receivers. Configure phase coded waveform with maximum-length sequence (MLS) using mlseq function and phased.PhaseCodedWaveform system object.

% Get waveform parameters
[N,tchip,prf] = helperWaveParameter(bw);

% MLS
mls = mlseq(N,1:numTgts);

% Configure emitter waveforms
wave = zeros(N,numTgts);
for s = 1:numTgts

    % Phase coded waveform with MLS
    waveform = phased.PhaseCodedWaveform('Code','Custom','CustomCode',mls(:,s), ...
        'SampleRate',fs,'ChipWidth',tchip,'PRF',prf);
    wave(:,s) = waveform();
end

Signal Propagation in Multipath Channel

Obtain received signals at multiple receivers via propagating emitter signals into the multipath channel.

% Configure transmitter and receiver
transmitter = phased.Transmitter('PeakPower',10);
receiver = phased.Receiver('NoiseFigure',10,'SampleRate',fs);

% Propagate source signals into multipath channel
[sigr,estNoisePow,anchorPos,tgtPos] = helperSignalPropagation(scenario,wave, ...
        transmitter,channel,receiver,rxPairs);

TDOA Estimation

Obtain TDOA measurements at each receiver using phased.TDOAEstimator based on the generalized cross-correlation phase transform (GCC-PHAT) algorithm.

% Configure a TDOA estimator
tdoaEstimator = phased.TDOAEstimator( ...
    'NumEstimates',4,'SampleRate',fs, ...
    'VarianceOutputPort',true,'NoisePowerSource','Input port');

% TDOA estimation
[tdoaEst,tdoaCRLB] = tdoaEstimator(sigr,estNoisePow);

For an anchor pair, plot its TDOA estimation spectrum, TDOA estimates, and true TDOAs. In the existence of noise, the estimated TDOAs can deviate from the true TDOA peaks on the TDOA spectrum. In a severe multipath environment, multiple false TDOAs can be detected.

% Plot TDOA Spectrum and TDOA estimates
figure
anchorPairIdx = 2;
[tdoaGrid,tdoaSpectrum] = plotTDOASpectrum(tdoaEstimator,'AnchorPairIndex',anchorPairIdx);

% Plot true TDOAs
hold on
helperPlotTrueTDOA(numTgts,anchorPos,tgtPos,c, ...
    tdoaGrid,tdoaSpectrum,anchorPairIdx);

Figure contains an axes object. The axes object with title TDOA Spectrum, xlabel TDOA (ns), ylabel Power (dB) contains 3 objects of type line. One or more of the lines displays its values using only markers These objects represent TDOA Spectrum, TDOA Estimate, True TDOA.

Source Localization

Before estimating position, we need to assign TDOA measurements to the correct target. However, there is no a priori information about which TDOA measurement belongs to which target. Also, the false detection of TDOAs due to multipath can make the assignment even harder. A brute-force way is to obtain all combinations of TDOA measurements to form all possible TDOA groups. This will lead to many TDOA groups for position estimation.

% Generate all TDOA combinations to form all possible TDOA groups
tdoaEstCell = num2cell(tdoaEst,1);
tdoaGroupTable = combinations(tdoaEstCell{:});
tdoaGroup = table2array(tdoaGroupTable);
tdoaGroup = tdoaGroup(sum(~isnan(tdoaGroup),2)>=3,:); % At least 3 TDOA estimates for localization

When we use all possible TDOA groups for target localization, most of the associated TDOA groups correspond to ghost targets. This issue is also known as the ghost target phenomenon. The TDOA position estimator tdoaposest can figure out part of the invalid TDOA groups using the triangle inequality principle described in the two figures below.

Triangle inequality principle

Left figure: The triangle inequalities indicate that the absolute value of the TDOA of a receiver pair is smaller than the propagation delay between the two receivers.
Right figure: The triangle inequalities indicate that the absolute value of the difference of TDOAs of two receiver pairs is smaller than the propagation delay between the two non-reference receivers.

The position estimates corresponding to the invalid TDOA groups with TDOA pairs not satisfying triangle inequalities are outputted as all-NaN column vectors from tdoaposest. The outputted non-NaN column vectors from tdoaposest are considered as potentially valid target position estimates. Although the triangle inequalities can help recognize a lot of incorrect TDOA association results, some invalid TDOA groups may still satisfy the triangle inequalities by chance and are not recognized as invalid. When the number of receiver pairs and the number of TDOA detections per receiver pair are large, many ghost targets can still be outputted from tdoaposest.

% Target position estimation
tgtPosEstAll = tdoaposest(tdoaGroup,tdoaCRLB,anchorPos);

% Display potentially valid target position estimates
tgtPosEst = tgtPosEstAll(:,~isnan(tgtPosEstAll(1,:)))

tgtPosEst = 3×2

   29.2217   29.7058
   45.6828   46.2927
   65.5990   66.3370

Comparing the potentially valid target position estimates to the true target positions, we can observe that many ghost targets in multipath environment are identified and removed.

% Display true target positions
display(tgtPos);

tgtPos = 3×1

   29.2207
   45.9492
   65.5741

Source Tracking based on Fused Detections

In this part, we present a workflow to associate more accurate TDOA groups using a static detection fuser implemented in staticDetectionFuser (Sensor Fusion and Tracking Toolbox) system object. The static detection fuser computes the likelihood of an associated TDOA group (a set of TDOA measurements from different receivers) to represent a true target versus a ghost target. To do this, for each possible TDOA group, the fuser uses tdoaposest included in helperTDOA2Pos to obtain the estimate of target position. Given the estimated position, the fuser uses helperFusedTDOA to obtain fused TDOA measurements and uses the fused TDOA measurements to compute the likelihood of each TDOA measurement originate from that target as compared to a false alarm. This likelihood is computed for each possible associated TDOA group and the most likely association is chosen using a multidimensional assignment algorithm. Compared to using triangle inequalities principle implemented tdoaposest alone, using fused TDOA detections from staticDetectionFuser can significantly reduce the number of ghost targets.

% Define a static detection fuser
tdoaFuser = staticDetectionFuser(MaxNumSensors=numReceivers-1, ...
    MeasurementFormat='custom', ...
    MeasurementFcn=@helperFusedTDOA, ...
    MeasurementFusionFcn=@helperTDOA2Pos, ...
    FalseAlarmRate=1e-10, ...
    DetectionProbability=0.99, ...
    Volume=1e-2);

To further improve the localization accuracy, we use a multi-object tracker based on the global nearest neighbor (GNN) assignment algorithm implemented in trackerGNN (Sensor Fusion and Tracking Toolbox) system object. Initialize the tracking filter as a constant-velocity extended Kalman filter (EKF) initcvekf (Sensor Fusion and Tracking Toolbox). For more information on source tracking based on fused detection, please refer to the section "Tracking Multiple Emitters with Unknown IDs" in the example Object Tracking Using Time Difference of Arrival (TDOA) (Sensor Fusion and Tracking Toolbox).

% Create a GNN tracker
tracker = trackerGNN(FilterInitializationFcn=@initcvekf, ...
                     AssignmentThreshold=20); 

% Configure tracking display
trackingDisplay = HelperTDOATrackingExampleDisplay(XLimits=boundLimit(1,:), ...
                                           YLimits=boundLimit(2,:), ...
                                           ZLimits=boundLimit(3,:));
trackingDisplay.Title = ''; % records new GIF if title is specified

Advance the tracking scenario to fuse and track target emitters.

while advance(scenario)
    % Current elapsed time
    time = scenario.SimulationTime;

Propagate waveform from each emitter to the passive receivers through multipath channels. Obtain TDOA estimates and TDOA estimation Cramer-Rao lower bound (CRLB) using phased.TDOAEstimator.

    % Simulate TDOA measurements
    [sigr,estNoisePow] = helperSignalPropagation(scenario,wave, ...
        transmitter,channel,receiver,rxPairs);

    % TDOA estimates and TDOA estimation CRLB
    [tdoaEst,tdoaCRLB] = tdoaEstimator(sigr,estNoisePow);

The TDOA estimation CRLB obtained above is the lower bound of the TDOA estimation variance in a single-target free-space environment. In multipath environment with GCC-PHAT TDOA estimation algorithm used, the TDOA estimation variance is larger than TDOA estimation CRLB. In this example, assume TDOA estimation variance is 100 times larger than TDOA estimation CRLB.

    % TDOA estimation variance
    tdoaVar = tdoaCRLB*100;

Organize TDOA estimates and TDOA variance into objectDetection (Sensor Fusion and Tracking Toolbox) format, which captures the TDOA estimates as the measurement and capture the TDOA estimation variance as the measurement variance.

    % TDOA detections
    tdoaDets = helperTDOADetection(scenario,tdoaEst,tdoaVar,rxPairs);

Use staticDetectionFuser to associate and fuse TDOA measurements into position detections.

    % Position detections
    posDets = tdoaFuser(tdoaDets);

Track fused source positions using trackerGNN.

    if ~isempty(posDets)
        % Update tracker with position detections
        tracks = tracker(posDets, time);
    end 

    % Update tracking display
    trackingDisplay(scenario, receiverIDs, tdoaDets, posDets, tracks);
end

Figure contains an axes object. The axes object with xlabel X (m), ylabel Y (m) contains 15 objects of type line, patch, text. One or more of the lines displays its values using only markers These objects represent Emitters (Targets), Passive Receivers (Anchors), Static Fused Position, Tracks, (history).

% Write new GIF if title is specified
writeGIF(trackingDisplay);

In the above animation result, the passive receivers are depicted with the upward-pointing triangles, and the target emitter is depicted with the downward-pointing triangle. Each hyperboloid describes the surface with the same range difference (the same TDOA) between the target emitter and a receiver pair. The dark solid trajectory is the true trajectory of the target emitter. The growing blue bold curve shows how the target emitter is tracked by the passive receivers. We can see that in the multipath environment, the PSL sensor network can localize and track the target emitter using fused detections. The number of ghost targets is small due to the help of the static detection fuser. A few false tracks are observed, but they do not exist for a long time.

Summary

This example shows how to localize and track target emitters in a 3-D space using a PSL sensor network. In particular, this example shows: 1) how to propagate unknown phase-coded waveforms into a multipath propagation channel; 2) how to estimate TDOAs and 3-D positions of target emitters in a PSL sensor network; 3) what the ghost target phenomenon is, and how to use the triangle inequality principle and a static detection fuser to reduce the number of ghost targets; 4) how to use a multi-object tracker to improve TDOA localization accuracy via tracking.

Reference

[1] Daniel E. Hack, Lee K. Patton and Braham Himed, Passive MIMO Radar Networks, in Chapter 19 of Novel Radar Techniques and Applications Volume 1: Real Aperture Array Radar, Imaging Radar, and Passive and Multistatic Radar, SciTech, 2017

[2] Mateusz Malanowski, Target Localization, in Chapter 8 of Signal Processing for Passive Bistatic Radar, Artech House, 2019

[3] T. Sathyan, A. Sinha and T. Kirubarajan, "Passive Geolocation and Tracking of an Unknown Number of Emitters," in IEEE Transactions on Aerospace and Electronic Systems, vol. 42, no. 2, pp. 740-750, April 2006

[4] Marco Compagnoni et al., "A Geometrical-Statistical Approach to Outlier Removal for TDOA Measurements," in IEEE Transactions on Signal Processing, vol. 65, no. 15, pp. 3960-3975, Aug. 2017

Helper Functions

helperRFParameter Function

function [fc,bw,fs,c] = helperRFParameter
fc = 10e9;                                       % Carrier frequency (Hz)
bw = 500e6;                                      % Bandwidth (Hz)
fs = bw;                                         % Sample rate (Hz)
c = physconst('LightSpeed');                     % Emission speed (m/s)
end

helperScatterPosition Function

function scatterPos = helperScatterPosition(numScatters,boundLimit)
xScatter = boundLimit(1,1) + (boundLimit(1,2)-boundLimit(1,1))*rand(1,numScatters);
yScatter = boundLimit(2,1) + (boundLimit(2,2)-boundLimit(2,1))*rand(1,numScatters);
zScatter = boundLimit(3,1) + (boundLimit(3,2)-boundLimit(3,1))*rand(1,numScatters);
scatterPos = [xScatter;yScatter;zScatter];
end

helperWaveParameter Function

function [N,tchip,prf] = helperWaveParameter(bw)
N = 2^nextpow2(1024)-1;                          % Number of chips in one phase-coded waveform
tchip = 1/bw;                                    % Chip duration (s)
tWave = N * tchip;                               % Modulation period (s)
prf = 1/tWave;                                   % PRF (1/s)
end

helperSignalPropagation Function

function [sigr,estNoisePow,anchorpos,targetpos] = helperSignalPropagation(scenario,wave, ...
        transmitter,channel,receiver,rxPairs)
% Find targets
platIDs = cellfun(@(x)x.PlatformID, scenario.Platforms);
anchors = scenario.Platforms(ismember(platIDs,rxPairs(:)));
targets = scenario.Platforms(~ismember(platIDs,rxPairs(:)));

% Anchor position and velocity
L = numel(anchors);
anchorpos = zeros(3,L);
anchorvel = zeros(3,L);
for l = 1:L
    anchorpos(:,l) = anchors{l}.Trajectory.Position;
    anchorvel(:,l) = anchors{l}.Trajectory.Velocity;
end

% Target position and velocity
K = numel(targets);
targetpos = zeros(3,K);
targetvel = zeros(3,K);
for k = 1:K
    targetpos(:,k) = targets{k}.Trajectory.Position;
    targetvel(:,k) = targets{k}.Trajectory.Velocity;
end

% Signal propagation
signal = transmitter(wave);
N = size(signal,1);
sigr = zeros(N,1,L);
% Loop for different passive receivers
for l = 1:L
    % Signal from target emitters to a passive receiver
    sigp = zeros(N,K);
    
    % Loop for different target emitters
    for k = 1:K
        % Propagate signal from each emitter to each passive receiver
        sigp(:,k) = channel(signal(:,k),targetpos(:,k),targetvel(:,k),eye(3), ...
            anchorpos(:,l),anchorvel(:,l),eye(3));
    end

    % Collect and receive signals from target emitters
    sigc = sum(sigp,2);
    sigr(:,1,l) = receiver(sigc);
end

% Noise power estimation
estNoisePow = reshape(vecnorm(fft(sigr)).^2/N,[1,L,1]);
end

helperPlotTrueTDOA Function

function helperPlotTrueTDOA(numTgts,anchorPos,tgtPos,c, ...
    tdoaGrid,tdoaSpectrum,anchorPairIdx)
% TDOA truth
tdoaTruth = zeros(1,numTgts);
for k = 1:numTgts
    rangeThis = norm(anchorPos(:,anchorPairIdx+1)-tgtPos(:,k)) ...
        -norm(anchorPos(:,1)-tgtPos(:,k));
    tdoaTruth(k) = rangeThis/c;
end

% Plot TDOA truth
[~,tdoaTruthIdx] = min(abs(tdoaGrid-tdoaTruth));
scaledTdoaTruth = engunits(tdoaTruth);
plot(scaledTdoaTruth,pow2db(tdoaSpectrum(tdoaTruthIdx)),'go', ...
    'LineWidth',1,'MarkerSize',10,'DisplayName','True TDOA')
legend('Location','southwest')
end

helperTDOADetection Function

function tdoaDets = helperTDOADetection(scenario,tdoaEst,tdoaVar,rxPairs)
% TDOA detection
tdoaDets = cell(0,1);
measParams = repmat(struct('OriginPosition',zeros(3,1)),2,1);
sampleDetection = objectDetection(scenario.SimulationTime,0,'MeasurementParameters',measParams);
platIDs = cellfun(@(x)x.PlatformID, scenario.Platforms);
for i = 1:size(rxPairs,1)
    sampleDetection.SensorIndex = i;
    r1Plat = rxPairs(i,1);
    r2Plat = rxPairs(i,2);
    plat1Pose = pose(scenario.Platforms{platIDs == r1Plat},'true');
    plat2Pose = pose(scenario.Platforms{platIDs == r2Plat},'true');

    % Fill detections
    tdoaEstThis = tdoaEst(:,i);
    tdoaEstThis = tdoaEstThis(~isnan(tdoaEstThis));
    effNumEst = numel(tdoaEstThis);
    thisTDOADetections = repmat({sampleDetection},effNumEst,1);
    for j = 1:numel(thisTDOADetections)
        thisTDOADetections{j}.Measurement = tdoaEstThis(j);
        thisTDOADetections{j}.MeasurementNoise = tdoaVar(i);
        thisTDOADetections{j}.MeasurementParameters(1).OriginPosition = plat1Pose.Position(:);
        thisTDOADetections{j}.MeasurementParameters(2).OriginPosition = plat2Pose.Position(:);
    end

    % From all receiver pairs
    tdoaDets = [tdoaDets;thisTDOADetections]; %#ok<AGROW>
end
end

helperFuseTDOA Function

function tdoa = helperFusedTDOA(pos, params)
% Recalculate TDOA based on fused position
r1 = norm(pos(1:3) - params(1).OriginPosition(1:3));
r2 = norm(pos(1:3) - params(2).OriginPosition(1:3));
tdoa = (r1 - r2)/physconst('LightSpeed');
end

helperTDOA2Pos Function

function varargout = helperTDOA2Pos(tdoaDets, reportDetection)
if nargin < 2
    reportDetection = false;
end

% Use at least 4 TDOA measurements for accurate localization
if numel(tdoaDets)>=4
    % Location of the known receivers
    anchorPos = zeros(3,numel(tdoaDets)+1);
    % Reference location
    anchorPos(:,1) = tdoaDets{1}.MeasurementParameters(2).OriginPosition(:);
    tdoaEst = zeros(1,numel(tdoaDets));
    tdoaVar = zeros(1,numel(tdoaDets));
    for i = 1:numel(tdoaDets)
        % Non-reference location
        anchorPos(:,i+1) = tdoaDets{i}.MeasurementParameters(1).OriginPosition(:);
        tdoaEst(i) = tdoaDets{i}.Measurement;
        tdoaVar(i) = tdoaDets{i}.MeasurementNoise;
    end

    % TDOA position estimate
    [pos,posCov] = tdoaposest(tdoaEst,tdoaVar,anchorPos);

    % Set position and position covariance of an invalid TDOA group
    if any(isnan(pos))
        pos = zeros(3,1);
        posCov = 1e10*eye(3);
    end
else
    pos = zeros(3,1);
    posCov = 1e10*eye(3);
end

if reportDetection
    varargout{1} = objectDetection(tdoaDets{1}.Time,pos,'MeasurementNoise',posCov);
else
    varargout{1} = pos;
    if nargout > 1
        varargout{2} = posCov;
    end
end
end