Perl and Air Traffic Control

The Perl Journal September, 2003

By Richard Hogaboom

Richard Hogaboom works for Solidus Technical Solutions Inc. at MIT Lincoln Laboratory on NASA/FAA research projects. He can be reached at [email protected].

Surface surveillance of taxiing aircraft has always been important, and we at Lincoln Laboratory (in Group 42, Air Traffic Control Systems), make its study part of our mission. The public normally thinks of aircraft surveillance as an air-to-air thing— a safeguard against mid-air collisions. But some of the most devastating accidents and losses of life have occurred on the ground between taxiing aircraft and landing or departing aircraft. This article is about how Perl can be useful in improving ground- surveillance systems that guard against such disasters.

Background

Ground surveillance is complicated by several issues that are not present in air surveillance. For most air surveillance, except for low-altitude targets, the sky is the background and is not signal reflective; but for ground surveillance, any radar covering an airport surface must deal with an enormous number of reflective targets, both stationary and moving. Terminal buildings, towers of various types, baggage-handling equipment, and aircraft-towing vehicles are all signal sources that a radar must discriminate against. Seaside airports even have to filter out the ocean itself (waves) or passing ships.

This imposes an enormous burden on any ground-surveillance scheme based on radar. The signal-processing and computational resources required are considerable. No matter how well designed such systems are, they still tend to have several shortcomings: The return signal on a target has no identification; because the radar is both a single source and return receiver, it can miss small targets shadowed by large targets; and the exact velocity of targets is difficult to determine because of the slow speeds. I've seen systems represent targets by moving colored blobs of varying sizes depending on recognition confidence. It's confusing. Sometimes you just have to look out the window.

There is an alternative solution. All commercial aircraft are required to have Mode S transponders. These devices "squit" signals every half-second that contain the Mode S address of the aircraft (24 hex bits, for example, a9f7c8), which uniquely transforms into the tail number. There are two antennae situated above and below the fuselage that broadcast alternately. In one form of squitter, the long type, the GPS position is also encoded, and receivers located about the airport can forward the information for controller display. The short type of squitter contains only the Mode S address. Most aircraft, at present, only broadcast the short squitter. The scheme devised to locate these aircraft from the short-squitter signal is called "multilateration." In this article, I'll describe an algorithmic implementation of this multilateration scheme in Perl. This project is from 1995, and reflects some language features from that time. Were I to code it today, I would probably do some things differently.

Theory

The idea is actually quite simple. Several receivers are located at various positions to the side of the airport surface. These receivers are positioned on buildings, the airport surface, or on towers. The ideal configuration is to have all the receivers on the surface of the airport at exactly the same height and in positions that form a lot of equilateral triangles with other receivers. Practical considerations such as the location of terminals and obstacles and the desirability of having the full runway and taxiway surface visible from the receiver antennas make equal height positioning difficult, however.

The target (aircraft) squitter signal is received almost simultaneously by all of the receivers. The "almost" is the important part. Each receiver has to detect the leading edge of the squitter signal to within about 10 nanoseconds (ns). You then compare the absolute arrival times at a pair of receivers to get the time difference of arrival. Since the signal is traveling at the constant speed of light to both receivers (with a little fudge factor for the index of refraction of air), this means a difference of arrival time translates into a path difference that is known.

If you remember covering conic sections in calculus, then you might remember that the loci of points that have a fixed difference distance between two points is a hyperbola (or a hyperboloid) of revolution, because the target may have a height above the plane of the airport surface. You then select another third receiver and find the time difference (distance) in signal arrival of the same squitter and calculate the second hyperboloid.

The intersection of the two hyperboloids will form a hyperbola in a plane perpendicular to the plane of the three receivers. One additional piece of information is needed to fixed the target position. For surface targets, this is the known height of the airport surface; and for airborne targets, their altitude. This solution depends on there being one common receiver among the two baselines. Obviously, with a triangle, there are three possible combinations of two baselines. Any of these three may be used. Bertrand T. Fang's paper "Simple Solutions for Hyperbolic and Related Position Fixes," (IEEE Transactions on Aerospace and Electronic Systems, Vol. 26, No. 5, Sept. 1990) describes this procedure and also describes solutions with a third piece of information other than height, as well as solutions of time sums rather than differences, which result in ellipsoidal solutions.

Clocks are not perfect, and this type of application requires very accurate clocks. It turns out that maintaining exactly synchronized accuracy is not necessary. All clocks drift to some extent, even atomic clocks. The problem is not to determine the time of day, but only the difference in times of arrival measured in nanoseconds. The solution to this is to provide a reference clock signal from a reference transponder to all the receivers. This transponder is situated in the center of the airport and visible to all the receivers. The reference transponder has its own clock. This clock drifts. It is not important how the reference transponder clock drifts, but how the baseline receiver clocks drift relative to it.

By utilizing the fixed known distances from the reference transponder to each of the receivers (and thus, the known time differences of the propagation of light through air and the times with which each receiver stamps the reference transponder signal), it is possible to do a clock fit and calculate the clock drift of each receiver clock. This drift is then taken into account when doing the time-difference calculations from the actual target signals. The process of position calculation using a triangle of receivers is referred to as "trilateration."

The several receivers are needed for a variety of reasons:

1. To provide full surface coverage.
2. To provide redundant receivers, because a squitter might be blocked by another target or its own aircraft wing or fuselage.
3. To provide redundancy against multipath signal reflections off the runway or other aircraft that delay and scramble the receiver input.
4. To provide several possible triangle combinations, some of which might provide better geometric characteristics than others.

MLDA (Multilateration Data Analysis)

The choice for the data analysis was Perl. Its impressive capabilities, in the form of built-in functions, made it superior to nawk or any shell. Everything about the language seemed well integrated. And I realized that the calculations, even for hundreds of aircraft, were not all that significant for the average workstation. This made the interpreted, quick-turn-around modify/run cycle of Perl attractive.

Atlanta Airport

The data collection part of the project was conducted at Atlanta International Airport in 1995. The receivers, built by Cardion Corp., were called Receiver/Transmitter (RTs) because they could receive squitters and transmit Mode S interrogations as well. There were five RTs and the single reference transponder in the center of the airport surface. Data was taken during normal airport operations and the results sent to us for analysis. All data was in ASCII text.

Because Atlanta is a busy airport, data from the RTs was collected for many aircraft at once. Each time a squitter was received at an RT, its leading edge was time-tagged with an integer 16-bit time with a granularity of 10 ns. A data record with this time, the Mode S aircraft ID, the RT ID, and a confidence indicator was transmitted to a central workstation. All RTs would also continuously receive squitters from the reference transponder. The central workstation would save all this, and the data would be shipped to us.

The mlda Script

The main Perl script is called "mlda" (located in the ml/da/ directory in the package of downloadable source code for this article, available at http://www.tpj.com/source/). I won't try to describe every line of code, but I will go through the script in enough detail to outline its structure and function, and where Perl comes to the fore and where I had to maneuver around it. In this narrative, I am assuming that the reader is following along in the mlda code.

Input Options

The script has only three options; -s to skip input file processing if you have a reference transponder file already; -i ModeSID if you want to trilaterate on a single target; and -t c2/c3/i2 if you want to change the solution method to two- or three-dimensional analytic or two-dimensional iterative. The single argument is the raw RT data. One of the first things I liked about Perl was the easy argument handling; no getopt()/getsubopt() complexities.

Initialization

BEGIN block processing is next (actually first, since this is done on program init, but it's next in the code). Various tunable parameters are set. These are used to give maximums or minimums for data file processing, such as a minimum number of records needed to do a reasonable solution or a maximum time difference between RT records in order to consider them from the same squitter from the same target. Some physics constants needed for the trilateration algorithm like the speed of light in air are set:

# c = 2.997925e8 m/s or 2.997925 m/10ns
$c = 2.997925;

# index of refraction - est. based on air, 20 deg. C, 75% hum.
$ir = 1.000350;

# adjust speed of light by index of refraction
$c = $c/$ir;

Next are the RT pairs for which differences in time of arrival are calculated, and the trilateration RT combinations used to do the final position fix. Hashes are used to store the index of the RT combination subscripted by the RT IDs of the pair. Finally, we see the RT locations and the known signal-propagation delays from the reference transponder (based on distance and the speed of light).

Input Data File

The RT data file is input. Perl, of course, makes short work of this. Unfortunately, I do not have any actual old data files from this project. I do, however, have another program, simpts, that is designed to generate a fake input file for a single linearly moving target with simulated squitter receptions from each RT with the start point, the direction, the number of points, the point separation, and a little jitter added. The output from simpts (in the simpts directory in the downloadable source code) is converted into Cardion input file format ready for the mlda script.

Each squitter received by an RT generates a record. These have to be merged so that all the data from a single squitter and a single aircraft is together for trilateration analysis. Generally, since the differences in time of arrival are on the order of hundreds to thousands of nanoseconds, and the squitters are broadcast only every half second, then the relevant records will be close together in the input file. Because reception is sometimes scrambled, there are often bit errors in the Mode S addresses. This makes a simple string comparison of two addresses unreliable. We took the approach of allowing two addresses to differ by a tunable parameter (a number of bits) and still be declared equal. Most aircraft have addresses that differ by a lot of bits. We chose 6 bits (out of 24) as a reasonable cutoff. Anything less, and the addresses were considered equal. The merging of low-confidence data into high-confidence data was a goal. The cntbitd() subroutine (See Listing 1) does this bit comparison. Perl is often thought of as a text processing language, but I've found that it can handle binary files or binary calculations with equal ease.

After confidence merging, the merging or blocking of same-squitter, same-aircraft data is done. First, targets are broken out into entries in a $targets{} hash with the Mode S address (six hex digits) as the key. Then all data records have their data join()ed into a string, which is pushed onto the array @$hex by the name of the address. I used:

@$hex[++$#$hex] = join(' ', $idh[0], $wt[0], $rt[0], $cf[0], $st[0]);

Why I didn't just use a push(), I don't know. If I were writing it today, I probably would. The Perl feature of symbolic named arrays is one of the most useful syntactic constructs imaginable. Next, loop over keys(%targets) and build a new record of all the merged data and write to a file with the address as the name.

Clock Smoothing

The next problem is to do clock smoothing. What's that? Well, no matter how accurate, all the RT clocks will have a slightly different time, and even more important, a slightly different drift. In order to measure differences in times of arrival to the granularity of 10 ns, it's important to take this into account. Other practical considerations arise. At times the reference transponder signal is blocked from an RT, sometimes a target squitter is received by only some of the RTs, and sometimes reception is scrambled so that a squitter is missed (if it's more than 6 bits off).

Another real difficulty stemmed from the specific clocks selected for the RTs: There were jumps in clock readings that were way outside the normal jitter. I refer to these here and in the code as outliers. The clocks were oven-controlled crystals. This was a big mistake. These clocks had the necessary accuracy, but they did not drift linearly, or they exhibited quantum jitter within the time periods of the several seconds over which it was necessary to clock fit. The more expensive rubidium clocks are not only more accurate, but they drift very linearly over several minutes, and would have been a much better choice.

I cannot tell exactly how much rubidiums would have improved the clock fits and, thus, the solutions, but I suspect that something of the order of a factor of three to a full order of magnitude in the average root mean square deviation could have been achieved. I had to deal with these outliers somehow. I chose to do two clock fits, one long (over the whole target file), and one short (over the several seconds necessary to eliminate jitter and missed data).

Back to the code. The first thing to do is delete the %targets entry for the reference transponder, a9f7c8, since its processing will be separate from all other targets; the necessary data is in the file a9f7c8. (See line 412 of the mlda script.) I then read in a9f7c8, which contains the arrival time stamps at each of the RTs for every squitter from the transponder, and calculate the difference of this time and the propagation time of light between the RT and the reference transponder. I then difference this time difference with the propagation difference times to other RTs (all mod 64K since the clocks roll over at 16 bits). Got that? The idea is to see how the clocks on RT pairs are drifting relative to each other from what the difference time of arrival should be from the fixed reference. The relevant line is:

$refd{$reft[$rt_cnt], $i, $j} = ( ( $tmp[$i*2+1] - $rtd[$i] ) 
                                - ( $tmp[$j*2+1] - $rtd[$j] ) ) % $k64;

This saves the difference of differences in a hash indexed by the reference clock time and the RT IDs. If any data is missing, a -1 is entered.

Now comes the C part. I needed to do some kind of smoothing fit on this data, both to get a more accurate measure of clock drift over time and to fill in missing reference transponder data from missed squitters at the RTs. I started out with the idea of a least squares fit, but the outliers made this difficult. I needed something more tolerant of really bad data that would give me some way of obtaining a reasonable first fit that would help eliminate the outlier data and be ready for a second short-term fit.

Numerical Recipes in C (5) (Cambridge University Press, 1993, ISBN: 0521431085) provided a good solution in medfit(), or Least Absolute Deviation Median Fit. This was much more tolerant of outliers. However, it was not in Perl, so I had to either recode or call it from Perl. I decided to call it. This was where I had to maneuver around Perl the most. I wrote the reference clock time $reft[$rt_cnt] and the $refd{$reft[$rt_cnt], $i, $j} to the clk_cor_raw file with prt_refd("clk-cor-raw"); and called lad-medfit with:

# do linear clock fit - use Numerical Recipes void medfit() routine
open LSF, "lad-medfit 0|" || die "mlda: lad-medfit open failure: $!";
( @clk_par = <LSF> ) == 10 || die "mlda: lad-medfit completion failure: $!";
close LSF;

This construct runs lad-medfit (with a 0 arg which says to long-term fit) and opens a pipe from which the program output can be read. These clock-fit parameters are written to the lad-par file for visual analysis. Outlier elimination is then done if any point is too far from the long-term linear estimate. Bad points are written to lad-out for visual analysis and are deleted by setting $refd{$reft[$i], $l, $m} = -1;. The new, corrected $refd{} is written to the clk-cor-out file. Then the script does a short-term clock fit with the number of points to fit over ($st_clk):

# do short term linear clock fit - use Numerical Recipes void medfit() routine
open LSF, "lad-medfit $st_clk|" || die "mlda: lad-medfit open failure: $!";
@clk_par = <LSF>;
close LSF;

I then append the clock-smoothed points to the lad-par (again for visual analysis), overwrite $refd{} with the smoothed data, thus eliminating any -1s, and write to the clk-cor-smt file. I won't get into the details of lad-medfit.c because this is all in C. It's in the lsf directory in the source code if you want to see the details.

Solution

Now it's possible to actually do trilateration. We trilaterate for each target file in keys(%targets). These files have the same data format as the reference transponder file, a9f7c8, except that they are for ephemeral targets moving into and out of coverage.

Trilaterations are performed on all triangles that have data, and the results are written to *.tls (trilateration solution) files. Three different types of solution are calculated: a two-dimensional analytical, a three-dimensional analytical, and a two-dimensional iterative solution (1,2,3,4 in "References"). All three are output to the *.tls files. The routine that calculates the solutions is ml($pt, @rtt);, where $pt is the master workstation time at which the clock corrections are calculated, and @rtt contains the RT time stamps for the squitter in question. The *.tls files will look like:

c3a     544292 0 1 4   1859.990    292.129
c3b     544292 0 1 4   1721.927    494.637
c3      544292 0 3 4   1667.606    571.296
c3      544292 1 3 4   1682.204    556.787
trk     5442921678.509    558.461  3 TK
c3      545342 0 2 3   1669.138    553.378
trk     5453421669.138    553.378  1 TK
c3      545743 0 1 2 Solution not obtainable - abs(R1p-R2p-R12) > 
                                      e and abs(R1n-R2n-R12) > e.
c3      545743 0 1 3 Solution not obtainable - abs(R1p-R2p-R12) > 
                                      e and abs(R1n-R2n-R12) > e.
c3      545743 0 1 4 Solution does not exist - time differences inconsistent.
c3      545743 0 2 3 Solution not obtainable - abs(R1p-R2p-R12) > 
                                      e and abs(R1n-R2n-R12) > e.
c3      545743 0 2 4 Solution does not exist - time differences inconsistent.
c3      545743 0 3 4 Solution does not exist - time differences inconsistent.

A normal output with valid solution will show the solution type (c3), the time stamp, the three RTs used in the solution, and the x, y-coordinates (North, South). If the solution fails, a diagnostic is given in place of x,y. If two solutions are admitted (sometimes, the solution is ambiguous) then both are given (c3a and c3b).

For c3 solutions, a tracker is run. A tracker is an algorithm that attempts to eliminate bad data points and smooth over good points. Most trackers take some weighted average over time and attempt to do a reasonable estimate of target movement. This tracker only looks at the solution space of a single squitter solution; that is, all the trilaterations over any triangle of RTs that had good enough data for a solution, and attempts to eliminate outliers and average over the remaining solutions.

The tracker will drop a track if too much time passes between points, and start a new track. The tracker will eliminate points that show a velocity jump that is too large. The BEGIN block has some tunable tracker parameters. After bad data is eliminated, then an average of the good data is taken and output to the tracker records. If all five RTs receive the squitter and all 10 (5 combination 3) triangle combinations have a good solution, then the average is over 10 points or solutions for the same squitter. Most of the time, this is not true, and tracker solutions average two to five points.

The geometry of the triangles and the target location relative to that triangle is important in determining position confidence. An equilateral triangle with the target in the middle is best; an elongated triangle with the target to the side is less accurate. This has to do with differences in time of arrival being smaller, the ambiguity in position being larger in the odd geometries, and the intersection of two hyperbolas being out on their wings rather than at their inflection points. These *.tls files are then parsed by other simple Perl scripts such as trkx and trix (also in the ml/da/ directory in the downloadable source for this issue) to get data and do plots.

Figure 1 gives you some idea of how all this works. There are three RTs located at x-, y-coordinates (0,0), (1,0), and (1,1). Hyperbola one is a result of difference time measurements between RTs (1,0) and (1,1). Hyperbola two is a result of a different difference time measurement between RTs (0,0) and (1,0). (Actually, their reflection in the lower half plane is also a valid solution, but gnuplot—with which this image was generated—will not plot multiple-valued functions. The intersection of the two curves is the solution. This 2D figure is a simplification of the real 3D case, however, the intersection of hyperboloids of revolution are difficult to visualize, and for a target in the plane, this representation is accurate. Also note the existence of multiple intersection points and thus, multiple solutions.

Some reasonable criteria for selecting the single true location are needed. Solutions are often eliminated because they are off the runway surface. If several trilateration solutions with different RT combinations are available, then the common solution would yield the true target location. And finally, the previous target location a second or so before the current solution would be very close. The targets cannot jump suddenly (with very high velocity).

Other Code

There are several other directories to the project. Look at the README file in the project main directory. The most important are the da (Data Analysis), lsf (Least Squares Fit), simpts (SIMulation PoinTS), and anal_sol (ANALytical SOLution). A lot of the code in these other directories is in C and represents experimental attempts at better solutions to parts of the problem.

Click here to download the updated MultiLateration.tar file.

Perl Advantages

I came up with a list of 10 Perl-language features that made it superior to its competitors for this project, at least at the time. These are:

1. Option handling.
2. Dynamic allocation, both initial and expansion of arrays and hashes.
3. Simple I/O.
4. Ease of handling both binary and text files.
5. Symbolic declarations (the @$hex[++$#$hex] stuff).
6. Hashes and their multisubscript comma-separated lists.
7. Ease of calling other processes (and getting the data back).
8. Rapid development cycle.
9. Wide variety of useful built-in functions that play well together.
10. Context-sensitive function returns.

This effort was my first significant development with Perl. Before this, I had used Perl for many trivial tasks and got to know the language well enough to feel comfortable using it in a more complex project like this. No doubt the code could be improved with more modern Perl syntax, but the underlying algorithms are sound.

References

1. Fang, Bertrand T. Simple Solutions for Hyperbolic and Related Position Fixes, IEEE Transactions on Aerospace and Electronic Systems, Vol. 26, No. 5, Sept. 1990.
2. Ho, K.C., Chan Y.T. Solution and Performance Analysis of Geolocation by TDOA, IEEE Transactions on Aerospace and Electronic Systems, Vol. 29, No. 4, Oct. 1993.
3. Friedland, B., Hutton, M.F. Journal of the Institute of Navigation, Vol. 20, No. 2, Summer 1973.
4. Boisvert, R. Dr. MIT Lincoln Laboratory, Unpublished Notes and Communications, 1995.
5. William H. Press, et al. Numerical Recipes in C, Cambridge University Press, 1993, ISBN: 0521431085.

TPJ

Listing 1

sub cntbitd
{
   local($len, $h1, $h2) = @_;
   local($xo, $i, $cnt, @pos);
   local($hex) = "0123456789abcdefABCDEF";

   return -1 if $len > 8 || length($h1) < $len || length($h2) < $len;
   return -2 if substr($h1, -$len, $len) =~ /[^$hex]/ || substr($h2, -$len, 
                                                         $len) =~ /[^$hex]/;
   $xo = hex($h1) ^ hex($h2);

   @pos = ();
   for ( $i=0,$cnt=0 ; $i<$len*4 ; $i++ )
   {
      if ( $xo & 1 )
      {
         $cnt++;
         push(@pos, $i);
      }
      $xo >>= 1;
   }
   unshift(@pos, $cnt);
   return @pos;
}