# Source Localization and Tracking with Passive Receivers

*Since R2024b*

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.

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

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

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.

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 = 'Source Tracking Using TDOA Measurements';

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

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