The TPL Triangulation Works Well Even With Relatively Large Errors

  I added a random error of 30% to the calculated intensity values to simulate measurement errors and increased the number of measurement per side to 101 to reduce the statistical error for the estimated position. The lines of position did not always meat at a point so for the "best estimated position" I chose the point, z_b, for which the sum of its distances from the lines was a minimum. It was easier to do this with Cartesian coordinates so I converted the complex numbers that I was using with z' being the position of a minimum for the inverse intensity along a side of the triangle and e' the direction of the corresponding line of position. P_i is a projection operator that gives the distance of point z_b from the line of position. The formulas for the calculating z_b are rather simple.

With random errors the curves for the inverse intensity are less well defined and I used a log scale for the inverse intensity to better show the spread of the data for the smaller values.

The example shown is for an unusual event with a relatively large error for the estimated position but shows how z_b fits inside a small triangle formed by the lines of position. A error of 100 meters in the estimated position would narrow down the search area considerably.

More data gives better results but if there are obstructions near the black box the data along a constant course may not be complete due to the TPL passing through shadows.

Using a TPL to Locate a "Black Box"

  The towed pinger locator used to search for the black boxes of Flight 370 was the U.S. Navy's TPL-25. To locate the position of a pinger it is towed behind a vessel along a constant heading and detects ping intensity along a traverse. What one can do is look for a minimum intensity along track. At this point the pinger position will be on a line normal to the track. Repeating this process for a number of different headings one can use the intersection of the lines of position as an estimate of the pinger's location. To illustrate how this works I assumed that the pinger was at a point within a triangle determined by three vertices. If one plots the inverse intensity of the detected pings for the points along a track one will get curves similar to the ones below that one can fitted to find their minimum points.

The inverse intensity was a better match to the quadratic functions used for the fit and as a result gave a better estimate of the lines of position. Ten points along each side of the triangle were used to make estimate of the minimum position. The three normal lines cross at a point.

If one is not very close to the pinger one may still get a line of position from the triangle and a rough estimate of how far away it is.

Sound Can Play Tricks On Us

  The formula for depth as a function of assumed position was derived for a model for sound that used the inverse square law for its propagation. This may not always be the case. When the distance from the source is large compared to the depth the model for propagation is closer to that a slab of material. The sound energy is confined to the region between the floor and ceiling defined by the ocean bottom and surface so the rate at which the intensity drops off will be closer to an inverse distance law. This is probably why whale songs travel such long distances. Another case is when the sound energy is confined by a tubular channel and the intensity can travel long distances with little change.

  So one has to be on guard against a pinger's environment channeling the ping's energy. There may be whisper points that give false positions for its location. And the sonar operators were probably correct in their decision not to claim that they had located the Flight 370 black boxes but rather chose to wait until they had more reliable physical evidence to go by.

Progress on Finding a Pinger Location

  I made a little progress yesterday on finding a pinger's location based on the ping intensity data in the simplified problem previously mentioned. The first result is that the pinger's depth for an assumed position, z₀, can be determined from the intensity data, a.

A grid size of 4 km with 3 points on each side was assumed. The second result is that the gradient of the intensity data can be found using a linear fit and that gives a direction for the source.

The matrix M just contains information on the relative positions for the data points. The gradient for 1/a was found since its function is a little simpler. There is a total of 9 data points in the grid used. The computed direction was off by about 1.6 degrees. Doing a quadratic fit allows one to get a better estimate of the direction relative to the corner used as the reference point.

So in theory it is possible to locate a pinger using just the intensity data. The reference point and direction gives one line of position for the location the source of the pings. A second set of data for an area some distance away would be needed for another line of position and their intersection would be an estimate of the pinger's location. In the real world the data is corrupted by errors but using a large number of data points can improve results. If enough ping data was collected it still may be possible to determine location of the black boxes of Flight 370.

Finding a Position Based on an Amplitude Function

  The search for the black boxes of Flight 370 managed to pick up some pings on a number of occasions but were unable to get good location for them. The depth and decaying batteries were complicating factors. Under favorable conditions one can use the amplitude of the detected pings to estimate a source location.

  On April 14th I came up with a method that located a source using just four sets of observed amplitudes and positions. In lieu of measured amplitudes computed values for an assumed amplitude function were used instead but the actual location of the source was not used to find a best fit. The bottom of the water was assumed to be constant so that the depth would be known and a value was assumed for the source strength. The situation is defined by the figure below. Complex numbers were used to minimize the number of terms in the amplitude function with d being the depth and s the distance of a surface vessel from the point directly over the source. The amplitude, a, was assumed to follow an inverse square law.

A bar above a complex number such as z indicates its complex conjugate which is found by replacing the imaginary component with its negative value. For an ordinary least squares fit we need the amplitude deviations, δ,  or the differences between the computed values and the "measured" values but that didn't work very well. So instead I tried fitting the relative amplitude deviations, ρ, found by dividing a calculated deviation by its measured amplitude at a position.  The formulas can also be simplified by using relative positions, ζ, and the source strength at the surface directly above the source, Φ₀.

The objective function which is what one tries to minimize is the sum of the squares of the relative amplitude deviations. If the measurements are taken for a small number of points on a grid one can estimate the center, ζ₀, of the amplitude function. One does not have to be directly over the source to locate its position. One of the data points was used as a reference position.

The fit found the position that was assumed for the source. The assumed value below are unprimed while the estimates have a prime on them. From the fit position one can compute heading and range from the reference position to the source location.

The geometry of the ocean bottom can act to focus the sound waves of the pings and alter their amplitudes. Changes in the density of sea water with depth can deflect the path of the sound waves much like the change in density of the air with altitude can alter the index of refraction at a given height and alter the observed position of a star. Under more complicated conditions a least squares fit would be more difficult but it still may be possible if one can collect more information of the ocean floor. One also needs to know the hazards in the neighborhood of the pingers if one wants to move about and collect objects from the ocean floor for investigation.

If the source strength is unknown the range of the source will be uncertain but one might still be able to use a grid search at well separated locations to get the relative headings and use intersecting arcs to estimate an actual position for the source.

Using Spherical Trigonometry To Work From Knows To Unknowns

  Spherical trigonometry uses the known sides and angles of a spherical triangle formed by great circle arcs to find those quantities that are unknown. In the figure below for the English Channel problem the labels for the station points are also used for the interior angles, M, R and X, of the spherical triangle and the sides opposite to these angles, x, r and m, are also indicated.

The unit vector tangent to a great circle arc at a point can be used to find the bearing to some other point on the arc. The trick is to eliminate the component of the second vector in the direction of the first. The two simple procedures shown below are very useful. They determine the direction of the tangent and the heading from the first point to the second. The functions E_r, E_θ and E_φ compute the radial and angular unit vectors from the angular coordinates of a point.

We started with the coordinates and headings for stations M and R. We need the heading of R relative to M to help determine the interior angles at M and R. The dot product of the two position vectors, e_M and e_R, are used to find the length of the connecting arc, x.

We now have enough information to calculate the interior angle at X using the law of cosines for angles. The values of X and x allow us to calculate the common ratio, ρ, for use with the law of sines and we can then find the lengths of sides r and m.

Either r or m can be used to find the position vector for station X. For m the formulas are,

e_n = N cos(hdg_MX) + E sin(hdg_MX)
e_X = e_M cos(m) + e_n sin(m).
E and N are the tangents to the sphere at M that give the direction of increase for φ and θ. The directions of the North Pole, the intersection of the Prime Meridian and the Equator and their easterly normal vector can be used to find the coordinates of station X.

The search that was done to find the intersection point wasn't really necessary. But it shows more directly that the arcs and their intersection are what determines the position for X.

Good Intersection Point for a Tangent Plane

  While trying to minimize the map distortions that result from treating the lines of latitude and longitude as straight lines in a plane I decided to do the 2D intersect in the plane tangent to the unit sphere at the coordinates of the War Memorial. The directions of east, E, north, N, and zenith, Z, there define this plane and the projections onto it. The deviation of intersection point for the tangent plane from that found for the intersecting arcs on the surface of the sphere was measurable in seconds of arc. Finding the required coordinate rotation can be a little tricky. Rotating the radial line through M to the z-axis in the plane defined by these two directions and then rotating the line in the direction of E to the x-axis while keeping the z-axis fixed works fairly well.

Supplemental (Apr 22): The positions of x'X and y'X were found for the intersection of the two straight lines determined by the heading angles for the two known positions. The third coordinate, z'X can be computed with the Pythagorean Theorem. The following calculation allows us to compare the components of X' with those of X and also the new coordinates. The basis B is the set of coordinate axes found for station M and acting on it with rotation R shows that they map onto the xyz-axes. D gives the 3D coordinates of the known positions, M and R, and the calculated positions, X and X', for comparison.

 A good reference for finding the rotation R is Kuipers, Quaternions and Rotation Sequences (1999), Sect. 4.6 Great Circle Navigation (p. 91).

Comparison of Spherical Triangulation with Plane Triangulation

  We can make a comparison of the results for triangulation on a sphere with triangulation in a plane. To do so we need to convert from angular distances to ordinary distances and take into account the change in the length of a degree of longitude with changes in latitude. The angles and distances are easier to convert and compare if we assume a unit sphere.

The unknown position X' for the plane triangulation is shifted about 1.1 minutes of arc from the previous spherical triangulation position X. The arc lengths, α₁ and α₂, and the distances from the known points, Δs_M and Δs_R, are similar in value but also differ slightly too.

So we see that due to the map distortions plane triangulation does not work as well at a distance of about 30 kilometers.

Edit (Apr 21): corrected the error on deg of longitude which used sine instead of cosine and got a better agreement of the results. Also made minor changes to the text.

Triangulating on a Sphere

  In Google Earth I selected two positions on the coast of the English Channel in the UK near Dover and another in France for a spherical geometry triangulation check. The positions and the arcs of the spherical triangle connecting them can be seen in the following figure (click to enlarge).

Setting placemarks in Google Earth allows to get the coordinates for the positions chosen to a fraction of a second of arc and the ruler can be used to read the headings of the connecting arcs to a small fraction of a degree. An estimate of the position of the Semaphore Station on the French side was estimated by determining the point of intersection for great circle arcs through the assumed known points on the English side and along the measured headings. A 3D plot of the spherical triangle is shown in the figure below. The arc lengths, α₁ and α₂, were found by doing a search for the values which minimized the distance between the two arcs.

The resulting position for X was within a few seconds of arc of the observed coordinates for the Semaphore station.

Geodesics on a Sphere

  As the bounds of a map increase the curvature of the Earth starts to deviate from a plane and this will distort the positions on the map even if we use a projection onto a plane tangent to the surface at the center of the map. The distances, angles from north and the direction of lines of sight are altered from those of plane geometry. The simplest solution is to use a sphere to represent the surface of the Earth and work with spherical geometry.

  The first question that we need to ask is how do we measure distances on the surface of a sphere? If we stretch a string between two points on the surface it becomes something like a straight line and we can used the string length as a measure of the distance. The string's arc is the arc common to the sphere's surface and a plane passing through the two points and the center of the sphere. This is a section of a great circle. The equator and meridian arcs that are used in spherical coordinates are great circles. The other lines of latitude are not but are arcs drawn at equal distances from the pole of the coordinate system.

  To show that the great circle is the minimum distance between two points we can compare the great circle length  with the lengths of other arcs connecting them. These arcs will have centers on a great circle passing through the center of the great circle connecting the two points and normal to it. To compare the arcs we need their lengths, Δs, and the maximum distance, β, of the arc normal to the great circle connecting the points.

When the lengths are computed and plotted we see that β = 0 corresponds to the minimum length which is that of the great circle.

So we see that great circles play the role of straight lines in spherical geometry. The great circle is a geodesic.

A Third Practical Positioning Example

  When one can bounce a signal off an object and measure the time it takes for the signal to travel to and from the object one can get an estimate of the range of the object. One can do this with the "echo" of a radar or sonar pulse and the strength of the returned pulse is improved if the object has a transponder. The strength of a pulse inside a cone follows the inverse square law and this also applies to a reflected pulse so the strength of the returned pulse follows an inverse 4th power law. The intensity with a transponder gives a stronger return signal but one needs power to generate the return pulse.

  Again one can triangulate the unknown position, X, given two positions, P₁ and P₂, and their ranges, r₁ and r₂. To find X one needs the distance of X, c, from the line connecting P₁ and P₂. If the closest point on the line to X is A then the distance, a, is determined as follows.

Using the formulas derived above one can easily find the distances a and c and consequently the positions A and X.

A Couple of Practical Positioning Examples

  The problem presented in the last blog was somewhat contrived but it illustrates two methods that can be used to find a line of position in a plane. The straight line assumed that we can determine the direction of north at the unknown position in order to measure the angle. If we can do this for one point we can also do it for another and we can find the position by determining the point of intersection for two lines.

If the direction of north is not readily available one can use three known positions and two angular separations to determine the unknown position. These are shown in the following figure.

For each angle we can find the center of a circle of position using the procedure described in the last blog then we can draw arcs from the common point to determine the location of X.

The line P₂X turns out to be perpendicular to the line between the centers Z₁₂ and Z₂₃ and the point X is the same distance from this line as P₂.

  One uses two circles to determine a position in Celestial Navigation. The latitude is a circle for which the distance from the North Pole is a constant. Latitude can be found by measuring the altitude of the Pole. Measuring the altitude of a star gives a second circle of position. The point on the Earth's surface directly beneath the star at the time of the altitude's measurement is needed and can be found using the star's right ascension and declination and sidereal time. The altitude of the star determines the distance of the second circle from the point beneath the star. And the intersection of these two circles gives an estimate of position. One needs the altitude of at least two points of the celestial sphere and the time of their measurement to determine a position on the surface of the Earth. A third star measurement or a computed position based on the heading, speed and time from the last known position is needed to determine which of the two possible points of intersection to use for the estimated position.

Determining One's Position on a Map

  With the recent difficulties encountered over the borders of Ukraine and the search for the black boxes of Flight 370 it might be useful to review some surveying problems. How can one determine one's position on a map? First of all one needs to determine the positions of locations relative known positions in a coordinate system to draw the map. This can be done by using triangulation in which one measures the angles of the unknown locations relative to pairs of known positions then computes the unknown positions using geometry or trigonometry. One can do something similar to determine one's position, P, on the map. Only two angles are needed. One can use the angle of one, A, relative to north and the angle between it and a second point, B.

The angle, φ, of the line PA relative to north is the same for all points on it including the point A so we have one line of position for P. A second line of position uses the fact that all points for which the angle between points A and B is the same angle θ lie on a circle. This is the inscribed angle theorem which is Prop. XXI in Bk I of Euclid's Elements.The trick is to find the center and radius of this circle. The center is equidistant from points A and B so it lies on a line bisecting the one between them.

The position for P on the line bisecting the line AB forms an isosceles triangle and we can assume that the center of the circle, Z, is at a distance x from the line AB. The exterior angle theorem (Euclid, Bk III, Prop. XXXII) tells us that the angle of A relative to the perpendicular bisector of line AB is equal to the sum of the two opposite angle of triangle APZ which is θ. From plane trigonometry we know that b/x = 2 tan(θ).

 We can also determine P geometrically using the following construction.

The point D is the intersection of two lines, one perpendicular to the line AB and a second drawn through the midpoint, C, of line AB at an angle θ relative to the perpendicular bisector, CZ. It gives us the distance x from the line AB. The intersection of the circle about Z through A and B with the line at angle φ from north through A is the required position.