WO2003065073A1

WO2003065073A1 - A method for field calibration of system parameters in a multibeam echo sounder system

Info

Publication number: WO2003065073A1
Application number: PCT/NO2003/000012
Authority: WO
Inventors: Jan Terje BJØRKE
Original assignee: Bjoerke Jan Terje
Priority date: 2002-01-28
Filing date: 2003-01-16
Publication date: 2003-08-07
Also published as: NO20020430L; DK200401289A; NO315766B1; NO20020430D0

Abstract

A seagoing vessel carries multibeam echo sounder transducer equipment for swath coverage of the sea floor for mapping thereof, covering a survey strip by traveling a given distance. Overlapping such strips are provided by running the vessel in opposite or crossing directions. The survey strips are transformed to a common reference system by resampling recorded reflection time measurement values to establish grid models. In order to execute field calibration of system parameters like offset in vessel roll (ϕ), pitch (ω), yaw (κ), vessel travel direction (x), direction (y) in horizontal plane orthogonal to said vessel travel direction, vertical depth direction (H), scale factor (my) in a swath direction, angular scale factor (mϕ) and time delay (τ), a difference between measured survey strips and an approximation of a true sea floor relief is minimized, using the equation: X = (N1 + P2)-1 (U1 + P2L2) where formula I and X is a vector of unknown calibration parameters, A1 is a coefficient matrix, P1 is a weight matrix of the measured depth values, L1 is a vector of the differences between a mean value H*(u) and Hi(u), where formula II and n is the number of overlapping strips in grid point u and Hi(u) is the depth value of strip i in u, L2 is a vector of a priori parameter values, and P2 is a weight matrix of parameters to be computed.

Description

A method for field calibration of system parameters in a multibeam echo sounder system

The present invention relates in general to echo sounding, and more particularly to so-called multibeam echo sounder systems. Precisely, the invention con- 5 cerns a method for field calibration of system parameters in multibeam echo sounder systems for seabed mapping, in which method a vessel carries multibeam transducer equipment mounted underneath said vessel for swath coverage of plural depth values in a plane generally perpendicular to the vessel travel direction, emitting sonar energy in a number of different 10 angles to the vertical and receiving reflected sonar energy along the same angles, completing a survey strip by travelling straight ahead a given distance, overlapping survey strips are provided by running said vessel in opposite or crossing directions, and the survey strips are transformed to a common reference system by re- i5 sampling recorded reflection time measurement values to establish grid models, the system parameters being parameters in a parameter group including offset in vessel roll, pitch, yaw, vessel travel direction, direction in horizontal plane orthogonal to the vessel travel direction, vertical depth direction, scale factor in a swath direction, angular scale factor, and time delay. 20 A relatively new system in the field of underwater acoustics is the so-called multibeam system. This system can be used to obtain the geometry of the seabed below the vessel on which it is mounted. Information about the geometry of the seabed is important and can be used for many applications, e.g. dredging or making bathymetric maps. 25 In contrast to the traditional echo sounder, which uses only one vertically aimed beam, multibeam is able to measure much more data during a shorter time interval, because it transmits many beams simultaneously.

Multibeam systems are usually mounted below a ship's hull and produce a fan-shaped beam pattern. Many beams are emitted at once using beamwidths of ₃o approximately 2° angle in cross-track direction.

Many different designs of multibeam systems exist, but they are all based on the same principle. A large array of transducers emits many beams of acoustic pulses simultaneously. Depending on the water depth, the speed of the vessel and the beam angle, it may be possible to cover swaths up to 7 to 8 times the water depth below the transducer.

Due to diffraction, only a small part of the emitted energy will return to the transducer. Diffraction causes acoustic energy to be spread in all possible directi- ons as a result of irregularities on the sea bottom.

The present invention does however not deal with the aspect of recording the single echo signals correctly, so for the purposes of the present invention we suppose that single footprint reflections are recorded correctly.

Looking at the appended fig. 1 , we see, thus, that a complete measurement session, called "a strip", consists of many successive swaths executed as the vessel sails along.

The multibeam system however, has to cope with a lot of extra errors, as a consequence of the fan-shaped beam-pattern. The propagation of these errors into the resulting 3D coordinates and the precision of the multibeam measurements are aspects which have been investigated thoroughly the last few years.

In the mapping community there is a growing need for accurate digital elevation models. Such data are frequently used to calculate slope, aspect, inter visibility and perspective views. The new generation multibeam echo sounder systems measure the terrain elevation with high precision, and GPS or other geodetic systems determine the position of the sensors with high precision. Therefore, there is a pressing need for accurate calibration of the systems considered. Since the multibeam echo sounder is integrated into complex computer systems, there are many error sources. Field calibration, therefore, is a very important tool to verify that all the static error sources are properly corrected and that the random mea- surement error is lower than a given threshold.

From Andre Godin: "Field procedures for the calibration of shallow water multibeam echo-sounding systems", Canadian Hydrographic Conference, June 1996, Halifax, Nova Scotia, Canadian Hydrographic Service and from the Internet article http://crunch.tec.army.mil/information/publications/multibeam/multibeam.html is known a calibration method that is related to the present invention. However, this previous method does not consider three-dimensional movement of measured points, so that it is not able to model the interrelationship between a change of depth and change of the planar coordinates (x, y). This is a clear drawback. Another method is known from an "Operator Manual: Neptune, Bathymetric Post-Processing", Kongsberg Simrad AS, P.O. Box 111 , N-3191 Horten, Norway. This is a calibration method based on a step by step procedure. The method computes roll, pitch, yaw (heading) and time delay. The drawback of this procedure is low accuracy due to the ignorance of the correlation between the parameters to be estimated, and that the method does not utilize neither crossing strips nor more than two overlapping strips.

However, the closest prior art known to the inventor of the present invention, is believed to be the article "A study of digital correlation methods in post pro- cessing of multibeam sounder data", Jan Terje Bjørke and Ivar Maalen-Johansen, Kart og Plan no. 5, 1987, Universitetsforlaget, Oslo, Norway, pp. 473-477. This article shows how a gradient method and a least squares method can be used to match multibeam strips. However, only two strips can be matched, and correlation is only defined relative to a strip, not to a true surface or an approximation thereto. Further, this prior art method suffers from constraints like moving images (i.e. measurement sets of the strips) as short distances as possible from their initial positions, during the correlation. Further, the prior art method does not even consider the problem of gradient computation, and does not describe how to compute adjusted depth values of overlapping strips. Nor does it deal with any evaluation of the reliability of computed parameters. Finally, this prior art method presents only a very general set of calibration parameters, and not all of the parameters are closely related to the operational characteristic of field calibration. (For example the parameter dsj can only be determined if the true seabed is known, which means that it is impossible to determine this parameter from field calibration, unless the seafloor is known from other measurements.)

The present invention sets out to remedy or at least alleviate the problems of the prior art as discussed above.

In accordance with the present invention there is provided a method for field calibration of system parameters such as defined in the introduction, and the met- hod is characterized by minimizing a difference between measured survey strips and an approximation of a true seabed relief, using the equation

where

U, = A^ P, L,

and X is a vector of unknown calibration parameters, Ai is a coefficient matrix, Pi is a weight matrix of the measured depth values, l_ι is a vector of the differences between a mean value H*(u) and Hj(u), where

∑H_t{u)

H * (u) ≡ -

and n is the number of overlapping strips in grid point u and Hj(u) is the depth value of strip i in u, l_₂ is a vector of a priori parameter values, and P₂ is a weight matrix of parameters to be computed.

Preferred embodiments of the method of the invention appear from the appended dependent claims 2-6.

Thus, it appears that the present invention provides primarily a method for field calibration of system parameters like offset in roll, pitch, yaw, x, y and scale, but the invention will also provide a method for computation of an adjusted terrain model from overlapping multibeam strips, i.e. block adjustment.

The method of the present invention will accept crossing strips as well as an infinite number of strips, i.e. limited only by the computer memory. The method does not rely on objects with known shape and position like some prior art solutions. Further, the inventive method is quite efficient for computing the gradient field of the terrain relief, i.e. with regard to processing time and accuracy. In the method, measured strips are compared to an approximation of the true surface, which is another advantage. The method of the present invention also presents a comprehensive method for evaluation of the reliability of parameters that have been computed, and the set of calibration parameters is well defined.

In the following, the method of the present invention shall be explained in more detail by means of descriptions of preferred embodiments, and it is at the same time referred to the appended drawings, of which Fig. 1 is a sketch showing the general geometry of the survey situation,

Fig. 2 is an illustration of crossing and overlapping survey strips, and

Fig. 3 illustrates the use of closest neighbour grid points for gradient computation. It is first referred to fig. 1 , in which a survey vessel 10 sails along a direction x while making echo sounder measurements of the sea bottom 20. A multi-beam technique is used, whereby a number of ultrasound rays 2 are emitted simultaneously in different directions in a "fan" transverse to the x direction, hitting the seabed in small "footprint" areas (4). The collection of footprints/rays per swath usu- ally is between 80 and 260.

The different rays 2 are produced by different transducers in a transducer array 8 arranged preferably at the bottom of the vessel 10. The transducer array position defines an origin of coordinates, however a moving origin.

The transducer array 8 fires "pings" (all transducers simultaneously) to cover successive swaths 6 as the vessel advances along the x direction. Swath width is usually comparable to water depth, often 2x - 6x water depth, or even wider, and travel distance between two successive swaths will be comparable to the footprint (4) diameter. However, all these parameters can be varied within wide limits. As the vessel 10 completes a predetermined survey distance along the x direction, a "strip" 12 is generated as a collection of all successive swaths 6.

In the present invention, a set of overlapping such strips 12 of terrain are used to compute the static offset of roll, pitch, yaw, x, y, scale and time delay. Horizontal scale factor dm_y is introduced as a common parameter for all the strips. This parameter is efficient to model errors in the sound speed. The method introduces another scale factor drriφ which is associated to the angle of the beams, i.e. dψ = φdnriφ. This angle scale factor is efficient to model the angular error of the outer beams.

The strips must run in different directions; opposite directions or crossing directions, see the example in Fig. 2. In Fig. 2 the vessel run directions are indicated by means of arrows. The four strips are indicated by reference numerals 12a, 12b, 12c, 12d.

The accuracy of the computed offset corrections depends on the number of the overlapping strips, their configuration and the relief of the seabed. Resampling to a common reference system

The strips are transformed to a common reference system by resampling the measurements to grid models. The grid procedure is crucial to the calibration procedure. Therefore, accurate interpolation like kriging (see for example Hans Wackemagel, "Multivariate geostatistics", Springer, Paris, second edition, 1998), an inverse distance polynomial method or Delauney triangulation must be applied.

The computation of the approximated surface

From the gridded strips an approximation H*(u) to the true seabed relief is computed as

∑H )

H * (u) ≡ (=1

(1 )

where n i s t he n umber o f t he o veriapping s trips i n g rid p oint u and H j(u) i s t he depth value of strip i in u. If less t han two strips overlap in u, the grid point is marked no value, else H*(u) is set to the mean value of the overlapping strips according to Equation 1.

If the true terrain relief H is known, this model can be introduced by defining

H^* = H (2)

The geometry model

We assume, as previously mentioned, a coordinate system with center (Xo, Yo, Ho) in the transducer; the positive X-axis in the strip direction, the Y-axis in the swath direction, i.e. across the strip direction, and the positive H-axis pointing to the center of the earth, see Fig. 1. The coordinates of the measured points are computed as

where R_κ and R_ω are the rotation matrices for yaw and pitch, respectively. Usually, the multibeam echo sounder system is designed so that ΛΓ « 0 and ω«0. Based on the previous assumption we derive the differential corrections from Equation 3 as

dx y dκ-h dω + dX_Q cly h dφ + y dm_y + dY₀ (4) dh - ydφ+h dm_h + dH₀

where

1 ) dφ is the roll offset; 2) dω is the pitch offset; 3) di is the yaw offset; 4) dx is the X₀ offset; 5) dy is the Yo offset; 6) dHo is the H₀ offset; 7) diri_y is a scale factor in the swath direction, i.e. in the Y direction, and 8) dm_h is a scale factor in the H direction.

Usually, dm and dH₀ cannot be determined. As we will show later, we will apply the dHj factor where i is the strip number, i.e. we compute the relative vertical translation of the strips.

The computation of the gradient The computation of the gradient is very crucial to the accuracy of the calibration. Experiments have shown that appropriate values of the gradient of the terrain relief in point (x_p, y_p) can be computed from the 9-polynomial as H*(x,y) = a₀₀+ aιo(x - x ) + aoι(y - y_P) + a₁₁(x - x_p)(y - y_p) + (5) a₂o(x - x_P)² +

aι₂(x - x_p)(y - y_P)² + a₂₂(x - x_p)²(y - y_p)²

where H* is computed from Equation 1 , i.e. we base the computation of the 9-polynomial coefficients on the mean model, derived from the overlapping strips. The coefficients of the 9-polynomial is computed from the nine neighbouring grid points 16 as shown in Fig. 3. Grid point p among the nine points is the point (x_p, y_p), and the grid is designated by reference numeral 14. From Equation 5 we compute the gradient in (x_p, y_p) from the first derivative of the 9-polynomial as

5H_X = a™ and δH_y = a₀ι (6)

Least squares approximation

When the partial gradients 5H_X and δH_y are known, we can compute the H component of a translation in X and Y. Therefore, we can combine the differential corrections in Equation 4 as

dH = SH_xy dK - 5H_xh dω +

(SH_yh - y)dφ +

5H_xdXo + 5H_ydYo + (7) dH₀ + δH_yy dm_y + h dm_h

From Equation 7 we can set up the equation for least-squares adjustment as

X = (A^τ PA)^"1A^TPL (8) where X is the vector of the calibration parameters, A is the coefficient matrix, P is the weight matrix and L is the discrepancy vector.

If it should happen that some of the parameters cannot be determined by the equations considered, the computation will fail. In order to prevent numerical problems of the kind considered, we apply the method of weighted parameters, see for example Alfred Leick, "GPS: Satellite surveying", John Wiley & Sons, Inc., USA, second edition, 1995, page 126. In the case of weighted parameters the offset parameters are computed from

X = (Ni + P₂)^'1 (Ui + P₂ L₂) (9)

where

, = A^ P, A,

IX = A_x ^τ P, L, (10)

where

a) X is a vector of the unknown calibration parameters. The forthcoming sec- tion will show how X is derived from Equation 12; b) Aι is the coefficient matrix. The forthcoming section will show how Ai is derived from Table 2; c) Pi is the weight matrix of the measured depth values. The forthcoming section will show how Pi is derived from Equation 13; d) Li is a vector of the differences between the mean value H*, which is defined in Equation 1 , and the measured value Hj(u) of grid point u in strip i, i.e.

Ln(u) = H^*(u) - Hi(u) (11)

e) L₂ is a vector of the a priori parameter values. Usually, L₂ = 0; f) P₂ is the weight matrix of the parameters to be computed. The forthcoming section will show how the initial values of P₂ are derived from Table 3. The weight matrix P₂ offers great flexibility in the design of the parameters to be computed. A large value of an element in P₂, for example, will fix the corresponding offset value to its parameter value in L₂. Knowledge of the standard devi- ation σ of the different parameters to be computed can be used to determine the weights of P₂, i.e. p = 1/σ².

The calibration parameters

The unknown parameters in vector X are related to a) parameters which are individual for the strips, b) parameters which are common for all the strips, c) other parameters; for example calibration parameters related to the speed of the vessel or parameters which are individual for the beams.

For operational purposes we have selected the parameters given by

where i is the number of strips. The parameters of Equation 12 are described in Tables 1 and 2. The weight of the measured depth values

The weight matrix Pi in Equation 10 describes the weight of the observed depth values. In principle, the correlation and the variances of the measurements can be introduced into P-j, but for operational purposes we will define

Pι = l (13)

that is, the diagonal elements of Pi are equal to 1 and the other elements of P-j are equal to zero.

The weight of the parameters

The weight matrix P₂ in Equation 9 offers great flexibility in the design of the parameters. High value of the weight, for example, will fix its parameter to the initial parameter value. The values given by Table 3 can be taken as guidelines for appropriate weight values.

The variance/covariance matrix of the parameters

The variance/covariance matrix Q_x of the unknown parameters is computed from the adjustment, see for example Alfred Leick, "GPS: Satellite surveying", John Wiley & Sons, Inc., USA, second edition, 1995, page 126, as

Qx = (Nι + P₂)-X (14)

The correlation of the parameters Very often, pitch offset and X offset are correlated to a high degree. Their correlation coefficient measures how separable the two parameters are. Generally, the correlation coefficient p_x,_x is derived from Q_x in Equation 14 as

where Xj and are two parameters. The reliability of the computed calibration parameters

The accuracy of the computed parameters is related to the number of strips, the amount of overlap between the strips, the run direction of the strips, the relief of the sea floor, and the sea depth. The computation of the standard deviation of the parameters from Equation 14 gives an estimate of the accuracy of the parameters, but it does not evaluate all the error sources of the computation as: (1 ) the accuracy of the resampling or (2) the accuracy of the gradient computation. Moreover, the correlation may have several local minima due to the shape of the seabed relief. However, the stability of the system of equations can be evaluated from simulation. The following example demonstrates the simulation developed for the purpose considered.

We assume the parameters given by Equation 12, see their description in Tables 1 and 2. In a first round, the calibration is computed following the ordinary calibration procedure as described by Equation 9. The first round gives vector Ci of calibration values. Before a second round is started, the measured strips are corrected according to the parameter values given by Ci + C°, where C° is a vector of offset values; appropriate values of C° are given in Table 4. The corrections of the strips are computed from the differential corrections given in Table 2. The calibration values from the second round are given by vector C₂. The difference between C₂ and C° is an expression for the reliability of the calibration. Ideally, C₂ should be equal to C°, but due to many error sources, the two may be different. If the relative value

Δ(k) = | C₂(k) - C°(k) | / C° (k) > t , (16)

the system is not able to compute reliable values of parameter k. Appropriate values of t and their interpretation are given in Table 5. The interpretation of the test above must consider the weights of the parameters. If Δ(k) « 1 , it may happen that parameter k is given high weight. If the parameter really is to be fixed to its initial value, the calibration is accepted, else the weight should be relaxed.

The test developed is very important, since it gives a comprehensive control of the calibration. Generally, it is hard to prove that a computer programme works as it should. The previous test however, covers both the implementation of the calibration procedure, the gradient model and the configuration of the measurements used in the calibration.

Block adjustment of the depth values The adjusted depth values can be computed from a similar procedure as the calibration. In a first round, the parameter values are calculated from the calibration. Then the measured strips are transformed according to the differential corrections given by Equation 4. When all the strips are corrected, i.e.,

(H_i,X_l, Y_l) = f(H_i ^' ,X_i ^' , Y: ) (17)

for = 1 ... n,

the new surface model is computed as the weighted mean model of the overlapp- ing strips. From the gridded strips the adjusted depth value Η(u) in grid point u is computed as

where n is the number of the overlapping strips in grid point u, Hj(u) is the depth value of strip I in u and pj(u) is a weight computed according to the inverse distance function as

p(u) = Ms² (19)

where s is the distance from the transducer to point u; see the geometry model in Fig. 1.

To summarize, the method of the present invention does not require knowledge of the true shape of the terrain surface, but overlapping strips must be desig- ned so that the effect of the error sources can be measured. Two totally overlapp- ing strips running in opposite directions, for example, will maximize the effect of static error on several of the system parameters.

Hence, it is possible to measure static errors in multibeam echo sounder systems from overlapping strips. The strips are matched together by the application of: (1) a gradient method and (2) least-squares adjustment with weighted parameters. All the error sources are computed simultaneously. Their variance/covariance is computed. In addition to the variance computation, the reliability of the computed parameters is evaluated from a simulation procedure. Since the measurement of the different strips does not refer to a common point pattern, the measurements are resampled to a regular grid. The resampling requires interpolation procedures as kriging or Delaunay triangulation.

It must be underlined that the description above has been centered around a certain embodiment of the method of the invention, however, several method steps can be varied. For instance, the computation of the gradient can be based on another subset of data points than the nine neighboring points, e.g., more than nine points combined with least squares computation of the coefficients in Equation 5. Further, the computation of the gradient can be based on other mathematical surfaces than the surface defined by Equation 5, e.g., other coefficients of the bi-cubic polynomial than the nine coefficients given by Equation 5.

Other calibration parameters than the parameters given by Equation 12 can be introduced, for example, parameters for the individual beams. This can be of particular interest for the outer beams.

Also, in the computation of the approximate and the adjusted surface in Equations 1 and 18, respectively, one must usually filter out unexpected depth measurements, i.e. outliers. The outliers can be removed from the data set by some kind of thresholding, e.g. by removing the measured depths which deviate more than a certain value from the computed mean value.

The differential corrections of Equation 4 assume that yaw and pitch are approximately equal to zero. If this is not the case, the differentiation cannot be based on the simplification considered. The general formulae of the differential corrections can be of interest for multibeam echo sounders mounted on submarines or other similar crafts. Finally, it is also possible to use a multiresolution grid based on e.g. "quadtree tessellation". This kind of grid can adapt the grid to the density of the measured depth points.

Claims

C L A I M S

1. A method for field calibration of system parameters in multibeam echo sounder systems for seabed mapping, in which method - a vessel carries multibeam transducer equipment mounted underneath said vessel for swath coverage of plural depth values in a plane generally perpendicular to the vessel travel direction, emitting sonar energy in a number of different angles to the vertical and receiving reflected sonar energy along the same angles, completing a survey strip by travelling straight ahead a given distance, - overlapping survey strips are provided by running said vessel in opposite or crossing directions, and the survey strips are transformed to a common reference system by resampling recorded reflection time measurement values to establish grid models, said system parameters being parameters in a parameter group including offset in vessel roll (φ), pitch (ω), yaw ( ), vessel travel direction (x), direction (y) in horizontal plane orthogonal to said vessel travel direction, vertical depth direction (H), scale factor (m_y) in a swath direction, angular scale factor (m_φ) and time delay (τ), c h a r a c t e r i z e d b y - minimizing a difference between measured survey strips and an approximation of a true seabed relief, using the equation

X = (Ni + P₂)^"1 (Ui + P₂ L₂) (9)

where

N, = A^ P, A,

U, = A^ P, L, (10)

and X is a vector of unknown calibration parameters, Ai is a coefficient matrix, Pi is a weight matrix of the measured depth values, Li is a vector of the differences between a mean value H*(u) and Hj(u), where

and n is the number of overlapping strips in grid point u and Hj(u) is the depth value of strip i in u, L₂ is a vector of a priori parameter values, and P₂ is a weight matrix of parameters to be computed.

2. The method of claim 1 , characterized in that gradients 5H_X, δH_y of the seabed relief in a grid point u = (x_p, y_p) are computed on the basis of the closest neighbouring grid points, by means of the equations

X~~ Xp X"-Xp y=y y=y_P

and

H^*(x, y) = a₀₀ +

a₀ι(y-y_P) +

a₂₀(x-x_P)² + a₀₂(y-y_P)² + a₂ι(x-x_p)²(y-y_P) + a₁₂(x-x_p)(y-y_p)² +

where aoo ■■■■ a₂₂ are coefficients that are computed from equation (5) and equation (1) using said neighbouring grid points.

3. The method of claim 1 , characterized in that said difference is minimized also taking into consideration a further system parameter which is a tide correction (dHoi) that is individual for each survey strip (i), the previously mentioned system parameters being common for all the strips.

4. The method of claim 1 , characterized in that differential corrections of system parameters are applied as follows:

parameter description diff. corr. coefficient

dHoi tide correction, strip i 1 dφ roll offset SH_yh-y dω pitch offset -5H_xh dK yaw offset δH_xy dX₀ linear offset along strip δH_x dY₀ linear offset across strip δH_y dmφ angular scale factor (δH_yh-y) arctan (y/h] dnriv linear scale factor δH_yy dτ time delay 9H_xv

v being vessel speed.

5. The method of claim 4, characterized by executing a reliability test for determined system parameter errors, first defining

t = I C₂(k)-C°(k) | / C° (k) (16)

where C° is a vector of simulated offset values regarding all system parameters k, and C₂ is a vector of calibration values resulting from two rounds of computing using equation (9), a first round giving a vector Ci of calibration values to be added to vector C° to provide Cι+C°, whereupon a second round of computation using equation (9) is based on the differential corrections, and subsequently evaluating t, t < 0,1 signifying unquestionable reliability and t > 0,9 signifying unquestionable non-reliability.

6. The method of claim 1 , c h a r a c t e r i z e d i n that an adjusted depth value H(u) in grid point u is computed by means of inverse distance weighting, in accordance with the equation

in which pι(u) is a weight computed according to the inverse distance function as

p(u) = s^"2 (19)

s being the distance from the transducer to point u,

thereby obtaining block adjustment.