Interpolation of Survey Coordinate Differences

This invention relates to determination of differences between ellipsoids associated with different datums in surveying activities.

A survey of locations at a survey site may involve use of location coordinates from previously-surveyed locations ("control survey points"), such as recognizable landmarks and other survey monuments, together with generation of location coordinates for newly-surveyed locations. Even where location coordinates for only newly-surveyed locations are used, these coordinates may not be fully consistent with each other when corrections are made for any errors introduced by use of a local coordinate system. The location coordinates of these locations, representing newly-surveyed locations and/or previously-surveyed locations, must be made consistent with each other in some "best fit" sense.

Survey of a chosen region often requires use of a globally defined survey ellipsoid, such as the NAD27, NAD83 or WGS84 ellipsoid, and of a locally defined ellipsoid that takes account of the local terrain. A globally defined survey ellipsoid is intended to provide a best fit, in some quantitative sense, to the terrain for the entire Earth. A locally defined ellipsoid, by contrast, is only intended to provide a best fit for a locally defined region and thus may provide a better fit over this limited region.

In a survey of a local region, location coordinates may first be obtained using a first, globally defined datum C1, after which these coordinates are re-expressed in terms of a second datum C2, which may be globally defined or locally defined, using a coordinate transformation T. Accurate coordinates for a plurality of survey control points are often available, and these survey control point coordinates are often used to "anchor" the transformation of location coordinates for other survey points. Coordinate differences between a point p(x,y,z) in the first, globally defined datum C1 and the corresponding point p'(x',y',z') in the second datum C2 are presently determined, using a point-by-point approach in which the group of survey control points chosen for "anchoring" can vary with the local terrain. FIG. 1 illustrates a vector difference .DELTA.r=(x'--x, y'--y, z'--z) between point coordinates (x,y,z) on a surface S1, defined using the datum C1, and the corresponding point coordinates (x',y',z') on a corresponding surface S2, defined using the datum C2. In the past, these coordinate differences were determined point-by-point on each of a plurality of sub-regions or patches that together make up a survey region SR, and these differences were often not continuous or consistent in moving from one sub-region to a contiguous sub-region. Thus, the coordinate differences had to be recomputed by applying the transformation T to each new point in the first datum; and parameters defining the transformation T might change in moving from one sub-region to a contiguous sub-region.

Several workers have considered the problems of representation of locations on a map or chart display and/or adjustment of location coordinates. Previous workers in this field often assume that the coordinates for locations of interest in a coordinate system, or in two or more associated coordinate systems, are consistent with each other. Further, the computations and coordinate manipulations are usually performed in a post-processing environment, rather than in a real time environment in the field at the time the survey measurements are made.

What is needed is an approach that allows real time processing, or post processing, to determine and apply a survey coordinate transformation between the first datum and the second datum in a chosen region that varies continuously and consistently with the point chosen in the region. Preferably, this transformation should be optimal in some sense over a grid of survey control points, and the transformation should depend only minimally or not at all on the particular grid chosen for such optimization. Preferably, this approach should allow use of any of a variety of generally defined coordinate transformations, each with its own set of transformation parameters that are to be optimized based on the survey control constraints imposed.

These needs are met by the invention, which provides a system for determination and application of a transformation, between a first, globally defined datum and a second datum over a chosen region, that varies continuously with the (arbitrary) choice of survey points within the region. A region is chosen, and a finite number N (N.gtoreq.2) of survey control points is identified within the region. A two-dimensional grid is defined over the region in the first datum; and for each point p=P.sub.G on the grid a transformation T, from first datum to second datum and having one or more adjustable parameters, is determined that is optimal, in the sense that the coordinate difference between first datum and second datum is "as close as possible" to the accurately known difference for each of K selected survey control points that is "closest to" the chosen grid point, where K.ltoreq.N. Each chosen grid point PG may have a different optimal transformation T associated with the grid point.

Three coordinate differences (one for each of latitude, longitude and height), computed using the transformation T at the chosen grid point P.sub.G, are associated with each grid point, thus forming three grids of difference values. For these three grids of difference values, a continuously variable interpolation function is determined that (1) substantially matches coordinate difference values (first datum versus second datum) at each grid point, either exactly or as closely as possible in some least qth power sense, and (2) extends each difference value grid to survey points that are not part of the original grid of points.

The transformation T between the first datum and the second datum may be one of the well-known four-parameter, six-parameter, seven-parameter or 12-parameter transformations used in survey work, as discussed in U.S. Pat. No. 5,614,913, issued to Nichols et al and incorporated by reference herein. Alternatively, the transformation T may be any other suitable adjustable parameter transformation between the two datums.

All these determinations of transformations and interpolation mappings between corresponding survey locations in a first datum and in a second datum can be made in the field in near-real time, or these determinations can be made by post-processing after the survey crew returns to the office. The invention will be most useful for in-the-field determinations, where errors in measurement, or in assumptions concerning a local coordinate system used for such measurements, can be promptly corrected if necessary.

FIG. 1 illustrates coordinate differences between first and second datums.

FIG. 2 illustrates application of the invention in one embodiment.

FIG. 3 illustrates a regional grid applied to a survey region.

FIG. 4 illustrates a grid of transformation differences between first and second datum coordinates, constructed according to the invention, to a plurality of coordinate grids, corresponding to the regional grid in FIG. 2.

FIG. 5 is a flow chart illustrating a method for practicing the invention.

FIG. 6 illustrates apparatus suitable for practicing the invention.

FIG. 1 (not drawn to scale) illustrates use of the invention in one environment. A user 11 employs survey apparatus 13 to conduct a survey at selected points in a chosen survey region SR. The survey apparatus preferably uses timed location determination (LD) signals received from a plurality of LD signal sources 15A, 15B, 15C, 15D, which may be part of a Global Positioning System (GPS), a Global Orbiting Navigational Satellite System (GLONASS), a modified low earth orbit (LEO) system, a Loran system or any other suitable LD system. However, the user apparatus 13 may also be conventional manual surveying apparatus (transit, survey pole, etc.), if desired. The survey is performed using a first, globally defined datum and associated first (global) ellipsoid 17 (and associated coordinate system C1) that is chosen as a best fit to the entire Earth's surface. The survey location coordinates are then to be transformed to a second datum (global or local) and associated second ellipsoid 19 (and associated coordinate system C2) that either (1) provides a better fit to at least part of the chosen region than does the first ellipsoid or (2) is more convenient for reporting the survey results. The survey location coordinates, referred to the first datum and referred to the second datum, will generally not coincide so that coordinate first datum-to-second datum differences are usually extant.

The survey region SR preferably includes N survey control points P.sub.SC-n (n=1, 2, . . . , N; N.gtoreq.3) (shown in FIG. 3) whose location coordinates (x.sub.SC-n,y.sub.SC-n,z.sub.SC-n) are accurately known for each datum, from performance of one or more previous surveys. The survey control point coordinates are often used to "anchor" or evaluate the accuracy of the transformation between the first datum and the second datum in terms of the accuracy of the transformation at each survey control point.

A grid G of spaced apart points P.sub.m,G, with corresponding location coordinates (x.sub.m,G,y.sub.m,G,z.sub.m,G) (m=1, 2, . . . , M), is imposed on the survey region SR in the first datum, as illustrated in FIG. 3. For a selected positive integer K (K.ltoreq.N), which may but need not vary with the choice of grid point P.sub.m,G, a set NSC.sub.m of K selected survey control points p'.sub.SC-k (k=1, 2, . . . , K) that are nearest or closest to the grid point p.sub.m,G is identified for each point on the grid G. Two nearest survey control sets, NSC.sub.m1 and NSC.sub.m2, may contain different points drawn from the survey control point set, based on the locations of the corresponding grid points p.sub.m1,G and P.sub.m2,G relative to the survey control points. For example, if K=3, the two sets NSC.sub.6 and NSC.sub.8 for the grid points p.sub.6,G and p.sub.8,G in FIG. 3 contain the survey control points {p.sub.SC-1, p.sub.SC-2, p.sub.SC-3 } and {p.sub.SC-1, p.sub.SC-3, p.sub.SC-N }, respectively.

The criterion used for "closeness" between a grid point p.sub.m,G and a survey control point p.sub.SC-n may be the conventional Euclidean distance ##EQU1## More generally, closeness may be defined as a qth power distance, such as ##EQU2## where the coefficients c.sub.x, c.sub.y and c.sub.z are selected non-negative real numbers, not all zero, and q is a selected positive number. With the choices q=2 and c.sub.x =c.sub.y =c.sub.z =1, the conventional Euclidean distance is recovered.

For each nearest control survey point set NSC.sub.m, corresponding to a grid point p.sub.m,G, a transformation T[p;p.sub.m,G ] of location coordinates of a point p in the first datum into location coordinates in the second datum is determined that minimizes a generalized distance between the coordinates of each survey control point p'.sub.SC-k in the set NSC.sub.m and the image under the transformation T of the survey control point p'.sub.SC-k, summed over all survey control points in the set NSC.sub.m. This generalized distance between NSC.sub.m and its image under T may be defined, for example, as a sum, ##EQU3## where the individual quantities d(p'.sub.SC-k,T[p'.sub.SC-k ;p.sub.m,G ];r) are defined in Eq. (1) or Eq. (2), for a particular choice of the positive exponent r, which need not coincide with the choice of the exponent q in Eq. (2).

The transformation T[p;p.sub.m,G ] may, for example, be any of the four-parameter, six-parameter, seven-parameter or 12-parameter coordinate transformations discussed in U.S. Pat. No. 5,614,913, incorporated by reference herein. Alternatively, the transformation T[p;p.sub.m,G ] may be any other suitable transformation of location coordinates.

A general form of linear transformation T on a three-dimensional coordinate system is one which carries Cartesian coordinates (x,y,z) in the first system into coordinates (x',y',z') in the second system, where the two sets of coordinates are related by ##EQU4## and where at least one of the coefficients h.sub.ij (i=1,2,3; j=1,2,3,4) is an adjustable parameter. In the most general transformation T of this sort, all 12 of the coefficients h.sub.ij are independently adjustable. In a well known special transformation, due to Molodensky, only seven parameters are independently adjustable, and the transformation T is expressed in the form ##EQU5## where L is a selected nonzero overall scale factor, .theta.12, .theta.13 and .theta.23 are selected angles of rotation about the third, second and first rectangular coordinates axes (or other selected axes), respectively, and .DELTA.x, .DELTA.y, and .DELTA.z are selected translation distances parallel to the three coordinates axes. Any other suitable three-dimensional transformation can also be used here.

A general form of linear transformation T on a two-dimensional coordinate system is one which carries Cartesian coordinates (x,y) in the first system into coordinates (x',y') in the second system, where the two sets of coordinates are related by ##EQU6## and where at least one of the coefficients h.sub.ij (i=1, 2; j=1, 2, 4) is an adjustable parameter. In the most general transformation T of this sort, all six of the coefficients h.sub.ij are independently adjustable. A special transformation of interest is given by ##EQU7## where L is a selected nonzero overall scale factor, .theta. is a selected rotation angle about a selected axis, and .DELTA.x and .DELTA.y are selected translation distances parallel to the two coordinate axes. Any other suitable two-dimensional transformation can also be used here.

The transformation T[p;p.sub.m,G ] is now applied to the grid point p.sub.m,G (p=p.sub.m,G), and coordinate differences

are computed for each grid point p.sub.m,G, where T[p.sub.m,G ].sub.w (w=x, y, z) indicates the value of the corresponding coordinate w under the transformation T[p;p.sub.m,G ] applied to the grid point P.sub.m,G. Here, s is another positive exponent that may, but need not, coincide with either of the exponents q and r. Recall that each grid point P.sub.m,G may have a different associated transformation T[p;p.sub.m,G ], or may have the same transformation form but with different transform parameters, or may have the same transformation.

A difference interpolation function IF(x,y,z) is now defined on one or more of the location coordinates x.sub.m,G, y.sub.m,G and/or z.sub.m,G (m=1, 2, . . . , M) of the grid points that provides a best fit to the function value(s) .DELTA.(x.sub.m,G), .DELTA.(y.sub.m,G) and/or .DELTA.(z.sub.m,G). The interpolation function IF is then extended to all coordinate values in the survey region SR and provides an estimate of the first datum-to-second datum coordinate differences that is continuous over the survey region SR.

The original grid G may be re-expressed as first, second and third coordinate grids, one of which is illustrated in FIG. 4 with w=x, y or z, corresponding to the respective x-, y- or z-coordinate values of the individual grid points P.sub.m,G. A separate difference interpolation function IF is optionally defined on each of the grids, corresponding to w=x, y and z. Alternatively, the same form for the interpolation function IF, with different parameters corresponding to each of the grids for w=x, y and z, can be used here.

Where a set of two-dimensional difference numbers d(x.sub.i,y.sub.i) (i=1, 2, . . . , M) is to be approximated or interpolated by a function IF(x,y), exact and approximate approaches are available to determine a suitable function f. The quantities d(x.sub.i,y.sub.i) may be the differences d(p'.sub.SC-k,T[p'.sub.SC-k ;p.sub.m,G ];s) between the coordinates of each survey control point p'.sub.SC-k in the set NSC.sub.m and the image under the transformation T[p;p.sub.m,G ] of the survey control point p'.sub.SC-k, as developed in the preceding. If the individual coordinates in the coordinate constraint set CS={x.sub.i,y.sub.i)} satisfy

the multilinear function ##EQU8## satisfies

for each value in the coordinate constraint set. If the inequalities (11) and (12) are not all satisfied, the multilinear function in Eq. (13) is modified to ##EQU9## where the avoidance requirements {j.noteq.i, x.sub.j .noteq.x.sub.i, y.sub.j .noteq.y.sub.i } for the product .PI. are interpreted as follows. Where two x-coordinate values x.sub.i and x.sub.j in the constraint set CS satisfy x.sub.i =x.sub.j with i.noteq.j, the ratio (x-x.sub.j)/(x.sub.i -x.sub.j) is formally removed (absent) from the product .PI.. Similarly, where two y-coordinate values y.sub.i and y.sub.j in the constraint set CS satisfy y.sub.i =y.sub.j with i.noteq.j, the ratio (y-y.sub.j)/(y.sub.i -y.sub.j) is formally removed (absent) from the product .PI., and similarly for the ration ((x-x.sub.j)/(x.sub.i -x.sub.j) with x.sub.i =x.sub.j and x.sub.i .noteq.x.sub.j). This insures that no factor in the denominator in the product .PI. is identically zero. In practice, if two x-coordinate values x.sub.i and x.sub.j (or two y-coordinate values y.sub.i and y.sub.j) in the constraint set CS are not identical but are sufficiently close (e.g., within a few centimeters or decimeters of each other), the ratio (x-x.sub.j)/(x.sub.i -x.sub.j) (or (y-y.sub.j)/(y.sub.i -y.sub.j)) is preferably formally removed from the product .PI..

A three-dimensional multilinear function, corresponding to the function IF'(x,y) in Eq. (15) may be defined by

The multilinear functions IF'(x,y) and IF"(x,y,z) defined in Eqs. (15) and (16) may be polynomials of degree up to M-1 in x and up to degree M-1 in y (and up to degree M-1 in z). If this polynomial degree is too large for the memory available, or for other reasons, the multilinear function IF'(x,y) or IF"(x,y,z) in Eq. (15) or Eq. (16) may be replaced by a lower degree polynomial, such as ##EQU11## where M.sub.x, M.sub.y and M.sub.z are suitably chosen non-negative integers. The coefficients g.sub.k,l,m in Eq. (17) may be determined by requiring a least squares fit for the coordinate pairs (x.sub.i,y.sub.i,z.sub.i) in the constraint set CS. That is, these coefficients are chosen so that a mean square error term ##EQU12## is minimized. This requires that the coefficients g.sub.k,l satisfy the relations ##EQU13## from which the coefficients g.sub.k,l,m can be determined by standard linear algebra techniques.

If interpolation is desired only over a two-dimensional region or surface, a suitable analogue of the polynomial IG"(x,y,z) in Eq. (17) is ##EQU14## and the analogue of the minimization relation, Eq. (19), becomes ##EQU15##

Another suitable interpolation function, defined on a two-dimensional surface, is a "bilinear" interpolation function IB(x,y;i,j) that is defined on each grid "rectangle" {x.sub.i .ltoreq.x.ltoreq.x.sub.i+1, y.sub.j .ltoreq.x.ltoreq.y.sub.j+1 } as

where the coefficients a.sub.i,j, b.sub.i,j, c.sub.i,j, and d.sub.i,j, are chosen so that

The coefficients a.sub.i,j, etc. are found by standard linear algebra techniques, using Eqs. (22) and (23). By matching the values at each corner of a rectangle, using a bilinear function, one insures that two bilinear functions, defined on contiguous grid rectangles, will substantially agree on the boundary line that the rectangles share.

Another suitable interpolation function, defined on a two-dimensional surface, is a "bicubic" interpolation function IC(x,y;i,j) that is defined on each grid "rectangle" {x.sub.i .ltoreq.x.ltoreq.x.sub.i+2, y.sub.j .ltoreq.x.ltoreq.y.sub.j+2 } as ##EQU16## where the coefficients a.sub.i,j, . . . , k.sub.i,j are chosen so that

The coefficients a.sub.i,j, etc. are found by standard linear algebra techniques, using Eqs. (24) and (25).

Another suitable method uses spline functions, such as a cubic spline, extended to two- or three-dimensional coordinates. Spline methods are discussed by M. J. D. Powell, Approximation theory and methods, Cambridge University Press, 1981, pp. 212-240, 261-265, 277-295. Any other suitable interpolation function can also be used here.

Each of the interpolation functions provides an interpolation function discussed or mentioned on the preceding for the differences d(p'.sub.SC-k,T[p'.sub.SC-k ;p.sub.m,G ];r), extended to all points in the survey region SR. Other interpolation functions can also be defined and used here.

FIG. 5 is a flow chart illustrating one method for practicing the invention. In step 41, the system selects or imposes a grid G of spaced apart points P.sub.m,G (m=1, 2, . . . , M) on the survey region SR. In step 43, for each grid point P.sub.m,G, a "nearest set" NSC.sub.m of (K) survey control points, from among all survey control points, is determined.

In step 45, a transformation T[p;p.sub.m,G ], depending upon the particular grid point P.sub.m,G selected, is determined, mapping a first coordinate system onto a second coordinate system, that (optionally) minimizes the generalized differences between the survey control points in NSC.sub.m and their images under the transformation T[p;p.sub.m,G ]. In step 47, for each grid point P.sub.m,G, a difference d(p.sub.m,G,T[p.sub.m,G ;p.sub.m,G ]) is determined between at least one location coordinate of P.sub.m,G and the corresponding location coordinate of the image under the transformation T[p;p.sub.m,G ]. In step 49, an interpolation function IF'(x,y) (or IF"(x,y,z))is determined that interpolates coordinate differences d(p,T[p;p.sub.m,G ]) over all points p in the survey region.

The interpolation function, IF'(x,y) or IF"(x,y,z) or IG'(x,y) or IG"(x,y,z) or IB(x,y;i;j) or IC(x,y;i;j) or other, is continuous in the coordinates x and y (and z) and is employed to approximate the actual distances between any point in the survey region in the first coordinate system and the corresponding point in the second coordinate system. The interpolation function thus obtained can be used to provide continuous and consistent interpolation values for the coordinate differences .DELTA.r=(.DELTA.x,.DELTA.y,.DELTA.z) between the first and second datums in a survey project.

FIG. 6 illustrates apparatus 61 suitable for practicing the invention. The apparatus 61 includes an LD module 63, including and LD signal antenna 65 and LD signal receiver/processor 67 that receive and process LD signals and determine the present location coordinates of the LD module 63. At a sequence of selected times, the present location coordinates are presented to a computer 69 that includes a microprocessor 71 and a memory unit 73 containing coordinates for one or more survey control points. Optionally, the LD signal receiver/processor 67 and the microprocessor 71 are part of the same data processing unit. The microprocessor 71 is programmed: (1) to obtain a coordinate grid {P.sub.m,G } for the survey region; (2) to obtain a coordinate transformation T[p;p.sub.m,G ] in the first datum between the first and second datums, using the survey control point(s) coordinates as discussed above; (3) to determine one, two or three coordinate differences d(p.sub.m,G,T[p.sub.m,G ;p.sub.m,G ]) between the first and second datums; and (4) to determine an interpolation function for the coordinate differences that is continuous across the survey region. Optionally, the computer 69 also includes a data/command entry system 75, such as a keyboard or light pen or touch-based screen, and/or a visually perceptible or audibly perceptible display module 77.

Optionally, the display module 77 is arranged to display at least one of the following: (1) a selected number (K) of survey control locations that are closest to a selected location; (2) at least one transformation parameter that describes the transformation T; (3) the difference between at least one coordinate of the selected location in the first coordinate system and the corresponding coordinate of the image under T in the second coordinate system of the selected location; and (4) an estimate, using the interpolation function, of the difference between at least one coordinate of the selected location in the first coordinate system and the corresponding coordinate of the image under T in the second coordinate system of the selected location. Optionally, the computer 69 can be programmed so that a selected number (K) of survey control locations that is closest to the present location is automatically displayed, in one mode of computer operation.

As an example illustrating application of the invention, assume that a valuable object was buried or otherwise concealed, intentionally or unintentionally, at some time in the past, using a local (second) coordinate system to determine its location at that time. A user who wishes to retrieve the object at a later time may have access to the burial location coordinates in the second system. However, the user may not have access to coordinate difference values of the second coordinate system, relative to a first coordinate system that is now relied upon by the user. Alternatively, the user may only have access to coordinate differences between the first and second coordinate systems at widely spaced apart coordinate reconciliation locations, and the burial location may have second system coordinates that are far from any of these reconciliation locations. If a coordinate difference interpolation function is available that reconciles the corresponding coordinates of the first and second systems over the region of interest, the user may determine the location coordinates of the buried object in the first system for any reasonable burial location coordinates in the second system.