NO315766B1 - A method for field calibration of system parameters in a multi-jet sonar system - Google Patents

Description

The present invention relates to sonar and in particular to so-called multi-beam sonar systems. The invention is more precisely about a method for field calibration of system parameters in multibeam sonar systems used for seabed mapping, where

a vessel carries multi-beam transducer equipment mounted under said vessel for swath coverage of many depth values in a plane mainly perpendicular to the vessel's direction of movement, with the emission of sonar energy at a number of different angles with the vertical and the reception of reflected sonar energy along the same angles, as a measuring strip is completed by driving straight ahead a given distance,

overlapping measuring strips are formed by driving said vessel in opposite or crossing directions, and

the measuring strips are transformed into a common reference system by resampling recorded reflection time measurement values to create lattice models,

where said system parameters are parameters in a parameter group that includes static corrections (deviation, error, offset) in roll (<J>), pitch (co), compass (k), vessel's direction of movement (x), direction (y) in the horizontal plane perpendicular on the said vessel's direction of movement, vertical depth direction (H), scaling factor (my) in the swath direction, angular scaling factor and time delay (r),

Multibeam systems are relatively new systems in underwater acoustics. Such systems are used to observe the geometry of the seabed beneath the vessel. Information about the geometry of the seabed is important in many applications, for example for dredging or for creating bathymetric maps.

In contrast to traditional single-beam echo sounders, multi-beam echo sounders are able to measure far more depth points during a given period of time, because multi-beam echo sounders operate with several beams at the same time.

Multi-beam systems are usually mounted under the boat's hull and produce a fan-shaped beam pattern. Many beams are emitted simultaneously, with a beam width of approximately 2° angle in the transverse direction.

There are many different multi-beam systems, but they are all based on the same principle. A long row of transducers/transducers sends out many beams of acoustic pulses at the same time. Depending on water depth, vessel speed and beam angle, it may be possible to cover swathes with a length of up to 7-8 times the water depth below the transducer.

Due to unevenness on the seabed, the acoustic energy will be spread in all possible directions (diffraction) so that only a small part of the emitted energy will return to the transducer.

The present invention is not concerned with issues related to registering the individual echo correctly, so we will therefore assume that the individual footprint reflection is correctly registered.

By looking at the attached figure 1, we see that a complete measurement session, which is called a measurement strip, consists of many successive swaths, carried out while the vessel is moving forward.

However, the multi-beam system is affected by many sources of error, as a result of the fan-shaped beam pattern. The effect of these error sources on the 3D coordinates the system outputs, as well as the accuracy of the multi-beam measurements, has been the subject of a number of investigations in recent years.

Within circles that carry out mapping, there is a growing need for accurate, digital elevation models. Such data is frequently used to calculate slope and outline of the aspect, perspective and see-through type. The new generation of multi-beam sonar measures the terrain height with great accuracy. GPS or other geodetic systems determine the position of the sensors with high precision. There is an urgent need for accurate calibration of multibeam sonar systems. Since such systems are integrated into complex computer systems, there are many sources of error. Field calibration is therefore an important tool for checking that all static sources of error are sufficiently corrected and that random measurement errors remain within given threshold values.

According to André Gddin: "Field procedures for the calibration of shallow water multibeam echo-sounding systems", Canadian Hydrographic Conference, June 1996, Halifax, Nova Scotia, Canadian Hydrographic Service and the Internet article: http:// crunch. tee. arrny. milfmformation/ publications/ multibeam/ multibeam. html describes a calibration method which is related to the present invention. However, the previously known method does not take into account the three-dimensional movement of the measured points. Therefore, it is not able to describe the connection between change in depth and change in the plane coordinates (x, y). This relationship represents a clear weakness of the method.

Another method is known from "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 calculates roll, pitch, compass and time delay. The disadvantage of this method is low accuracy due to the fact that the correlation between the parameters to be estimated is not taken into account, and that the method does not utilize crossing stripes or more than two overlapping stripes.

The method found to be closest to the invention is believed to be the article "A study of digital correlation methods in post processing of multibeam sounder data", Jan Terje Bjørke and Ivar Maalen-Johansen, Kart og Plan no. 5, 1987, University Press, Oslo, Norway, pp. 473-477. This article shows how a gradient method and least squares method can be used to match multibeam stripes. However, here it is allowed that only two strips can be matched, and correlation is defined only in relation to a strip, not in relation to the true terrain surface or an approximation thereof. This earlier method also suffers from the fact that conditions of a type have been introduced in the calculation where images (ie measurement sets of stripes) are moved the shortest possible distance from their starting positions during the correlation. Furthermore, the previous method does not even describe how the gradients or the equalized depth values between overlapping stripes should be calculated. Nor does it describe how the credibility of the calculated parameters should be derived. Finally, it is the case that the former method specifies a very general set of calibration parameters where not all parameters are equally well related to the operational conditions a field calibration is subject to (for example, the parameter ds; can only be determined if the true seabed is known, something which means that it is impossible to determine this parameter from field calibration unless the seabed is known from other surveys).

The present invention will improve or at least remedy some

of the weaknesses of the previous method discussed above.

In accordance with the present invention, a method for field calibration of system parameters as defined in the introduction has been provided. The procedure is characterized in particular by: minimizing the difference between measured stripes and an approximation to a true seabed topography by using the equation and X is a vector of unknown calibration parameters, A1 is a coefficient matrix, Pi is a weight matrix for the measured depth values, Li is a vector of the differences between a mean value H<*>(u) and H,(u), where

and n is the number of overlapping stripes at grid point u and Hj(u) is stripe in its depth value at point u, L2 is a vector with a priori (given in advance) parameter values, and P? is a weight matrix for parameters to be calculated.

Preferred embodiments of the invention appear from the attached non-independent patent claims 2-6.

It thus appears that the present invention mainly provides a method for field calibration of system parameters such as static deviations (offset) in roll, pitch, compass, x, y and scaling, but the invention also offers a method for calculating a leveled terrain model from overlapping multi-beam strips, i.e. .block smoothing.

The method of the present invention accepts intersecting stripes and an infinite number of stripes, i.e. it is limited only by the computer's storage capacity. The method does not depend on objects with a known shape and position, like some other previous similar solutions. Furthermore, the method according to the invention is effective when it comes to calculating the gradient field of the terrain topography, i.e. with regard to processing time and accuracy. In the method, the measured stripes are compared with an approximation to the real surface, and this is a further advantage.

The method according to the present invention also exhibits a comprehensive method for evaluating the credibility of the calculated parameters, and it uses a well-defined group of calibration parameters.

In what follows, the method according to the present invention will be explained in more detail by means of descriptions of preferred embodiments, and reference is also made to the attached drawings, where

figure 1 is a sketch illustrating general geometry in a surveying situation,

figure 2 is an illustration of crossing and overlapping measuring strips and

Figure 3 illustrates the use of nearest neighbor grid points for gradient calculation.

Reference is first made to Figure 1, where a survey vessel 10 sails along direction x while sonar measurements of the seabed 20 are made.

multi-beam sonar principle is used where a number of ultrasound beams are sent out simultaneously in different directions in a "fan" perpendicular to the x-direction. The aforementioned rays hit the seabed in small "footprint areas" (4). The number of footprints or rays in a swath usually varies between 80 and 260.

The different beams 2 are formed by different transducers (transducers) in a transducer row 8 preferably mounted in the bottom of the vessel 10. The position of the transducer row defines a reference point for the coordinate system, but note a moving reference point.

The transducer array 8 fires "pings" (all transducers simultaneously) to cover successive swaths 6 as the vessel moves in the x-direction. The swath length is usually proportional to the water depth, often 2x - 6x the water depth, or even further swaths. The distance between two consecutive swathes is often related to the size of the footprint (4). However, all these parameters can be varied within wide limits.

As the vessel 10 completes a predefined measurement line in the x-direction, a strip 12 is formed as a collection of all successive swaths 6.

In the present invention, a number of such strips 12 which overlap each other are used to calculate static deviations (offset) in roll, pitch, compass, x, y, scale and time delay. Horizontal scale factor dmy is introduced as a common parameter for all strips. This parameter is used to model errors in the speed of sound. The method also introduces another scaling factor dm$ which is associated with the angle of the rays, i.e. d+= §dm$. This angular scaling factor is effective in modeling the angular error of the outermost rays.

The stripes must be run in different directions; opposite directions or crossing directions, see the example in figure 2.1 figure 2 the vessel's driving directions are shown with arrows. The four strips are indicated by reference numbers 12a, 12b, 12c and 12d.

The accuracy of the calculated offset corrections depends on the number of overlapping stripes, their configuration and the relief of the seabed.

Conversion to a common reference system

The stripes are transformed into a common reference system by recalculating (resampling) the measurements to grid models. The grid procedure is crucial for the calibration calculation. Therefore, it is necessary to rely on accurate interpolation such as kriging (see for example Hans Wackernagel, "Multivariate geostatistics", Springer, Paris, second edition, 1998), polynomial method with inverse distance weight or Delauney triangulation.

Calculation of the approximate terrain surface

From the grid models of the strips, an approximation H<*>(u) to the true seabed surface can be calculated by

where n is the number of overlapping stripes in grid point u and Hi(u) is the depth value of u in stripe i. If less than two stripes overlap each other in u, the grid point is marked as not measured, otherwise H<*>(u) is set to the mean value of the overlapping stripes using equation 1.

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

The geometry model

We assume, as previously mentioned, a coordinate system with center (X0, Y0, H0) 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 points towards the center of the earth, see figure 1. The coordinates of the measured points are calculated by

where RK and Ra are the rotation matrices for compass and pitch respectively. Usually the multibeam sonar system is designed so that k«0 and<g>>*0. Based on the previous assumptions, we can derive the differential corrections from equation 3 by

where

1) d<{> is roll error; 2) deo is pitch error; 3) dK is compass error; 4) dx is X0 error; 5) dy is Y0wrong; 6) dH0er H0fault;

7) dmy is the scaling factor in the shear direction, i.e. in the Y direction, and

8) dmher scaling factor in the H direction.

Usually it is not possible to determine dmhand dHo. As we will see later, we will use the factor dHj (where i is the stripe number) to calculate the relative displacement of the stripes in the vertical direction.

Calculation of the gradient

Calculation of the gradient is very critical to the accuracy of the calibration. Experiments have shown that applicable values for the gradient of the terrain relief at the point (xp, yp) can be calculated from the 9-polynomial by

where H<*>is found by equation 1, i.e. we base the calculation of the 9-polynomial coefficients on the mean model derived from overlapping strips. The coefficients of the 9-polynomial are calculated from the nine neighboring grid points 16 as shown in figure 3. Grid point p among the nine points is given by (xp, yp), and the grid is indicated by reference number 14. From equation 5 we find the gradient in (xp, yp) from the first derivative of the 9-polynomial at

Least squares approximation

When the partial gradients dHx and dHy are known, we can calculate the H component of a displacement in X and Y. Consequently, we can combine the differential corrections in equation 4 as

From equation 7, we can set up the equation for least squares smoothing where X is the vector of 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 listed, the calculation will break down. To prevent numerical problems of this type, we introduce the method of weighted parameters, see for example Alfred Leick, "GPS: Satellite surveying", John Wiley & Sons, Inc., USA, second edition, 1995, page 126.1 the case of weighting is performed the calculation of the awiks parameters from

there

a) X is a vector with the unknown calibration parameters. In the following section, it will be shown how X is derived from equation 12; b) A-\ is the coefficient matrix. In the following section, it will be shown how At is derived from table 2; c) Pi is the weight matrix for the measured depth values. In the following section, it will be shown how Pi is derived from equation 13; d) Li is a vector of the differences between the mean value H<*> (as defined in equation 1) and the measured value Hj(u) for grid point u in strip i, i.e. e) L2 is a vector of a priori parameter values. Usually l_2 = 0 is set; f) P2 is the weight matrix for the parameters to be calculated. In the following section, it will be shown how the initial values in P2 are derived from

table 3.

The weighting matrix P2 offers great flexibility in the design of the parameters to be calculated. A large value on an element in P2, for example, will lock the associated awiks value to its parameter value in L2. Knowledge of the standard deviation a of the various parameters to be calculated can be used to determine the weights in P2, i.e. p = 1/cT<2>.

The calibration parameters

The unknown parameters in vector X are related to

a) parameters that are individual to the strips,

b) parameters that are common to all the strips,

c) other parameters such as calibration parameters related to the speed of the vessel or parameters that are individual to the beams.

For operational purposes, we have chosen parameters given by

where i is the number of stripes. The parameters in equation 12 are described in tables 1 and 2.

The weight of the measured depth values

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

ie, the diagonal elements of Pi are equal to 1, and the other elements of Pi are equal to zero.

The weight of the parameters

The weighting matrix P2i equation 9 provides great flexibility in the design of the parameters. A large weight value, for example, will lock its parameter to the initial parameter value. The values given in table 3 are examples of appropriate weight values.

The variance/covariance matrix of the parameters

The variance/covariance matrix Qxtil the unknown parameters is calculated from the smoothing, see for example Alfred Leik, "GPS: Satellite surveying", John Wiley & Sons, Inc., USA, second edition, 1995, page 126, as

The correlation of the parameters

It often happens that pitch errors and X errors are highly correlated. Their correlation coefficient measures how separable the two parameters are. The correlation coefficient px<x can generally be derived from Qxi equation 14 by

where X, and X are two parameters.

The reliability of the calculated calibration parameters

The accuracy of the calculated parameters is related to the number of stripes, the degree of overlap between the stripes, the direction of travel of the stripes, the relief of the seabed and the sea depth. Calculation of the standard deviation of the parameters from equation 14 gives an estimate of the accuracy of the parameters, but the calculation does not evaluate all sources of error such as: (1) the accuracy of the conversion (resampling) or (2) the accuracy of the gradient calculation. Furthermore, the correlation can have several local minima depending on the relief of the seabed. However, the stability of the system of equations can be evaluated by simulation. The following example demonstrates the simulation developed for that purpose.

We assume the parameters given by equation 12, see their description in tables 1 and 2.1 first the calibration is calculated in the ordinary way described in equation 9. This first calculation round gives vector C, of calibration values. Before a second calculation round starts, the measured stripes are corrected according to the parameter values given by Ci+ C°, where C° is a vector of deviation values; appropriate values for C° are given in table 4. The correction to the strips is calculated from the differential corrections given in table 2. The calibration values from the second round are given by vector C2. The difference between C2 and C° is an expression of the credibility of the calibration. Ideally, C2 should be equal to C°, but due to many sources of error, the two can be different. If the relative value

is the system unable to calculate credible values for parameter k. Suitable values for t and their interpretation are given in Table 5. The interpretation of the above test must take into account the weight of the parameters. If A(k) « 1, it can be thought that parameter k is given a lot of weight. If the parameter is really to be locked to its initial value, the calibration is accepted, otherwise the weight must be reduced.

The test that has been developed is very important, because it provides a comprehensive check on the calibration. It is generally difficult to prove that a computer program works as it should. However, the preceding test covers both the implementation of the calibration procedure, the gradient model and the configuration of the measurements used in the calibration.

Block smoothing of depth values

The equalized depth values can be calculated from a corresponding procedure used in the calibration. In the first instance, the parameter values are calculated from the calibration. The measured strips are then transformed according to the differential corrections given in equation 4. When all the strips have been corrected, i.e.

the new relief model is calculated as the weighted average model of the overlapping stripes. From the grid model of the strips, the equalized depth value H(u) at grid point u is calculated by where n is the number of overlapping stripes at grid point u, Htfu) is the strip at its depth value at point u and pj(u) is a weight calculated according to the inverse distance function provided by

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

In summary, the method according to the present invention does not require knowledge of the true terrain surface, but overlapping strips must be configured so that the effect of the error sources can be measured. For example, two completely overlapping stripes running in opposite directions will maximize the impact of static errors on several of the system parameters.

Consequently, it is possible to measure the static error of the multibeam sonar systems from overlapping streaks.

The strips are fitted using: (1) a gradient method and (2) least squares smoothing with weighted parameters. All error sources are calculated simultaneously. Their variance/covariance is derived. In addition to the variance calculation, the credibility of the calculated parameters is evaluated in a simulation.

Since the measurement of the various stripes does not refer to a common point pattern, the measurements are converted to a regular grid. The resampling requires interpolation procedures such as kriging or Delaunay triangulation.

It must be emphasized that the above description is centered on a certain embodiment of the method according to the invention. However, several of the steps in the procedure can be varied. For example, the calculation of the gradient can be based on a different subset of data points than the nine neighboring points, for example more than nine points can be used in combination with the least squares method for deriving the coefficients in equation 5. Furthermore, the calculation of the gradient can be based on other mathematical surfaces than the one defined in equation 5, for example other coefficients in the bicubic polynomial than the nine coefficients given by equation 5.,

Other calibration parameters than the parameters given by equation 12 can be introduced, such as parameters for the individual beams, which can be of particular interest for the outer beams.

In the calculation of the approximated and the leveled surface in equations 1 and 18, respectively, you usually have to filter out unexpected depth measurements, i.e. atypical values/gross errors. The gross errors can be removed from the data set by certain thresholding techniques, for example by removing the depth measurements that deviate more than a given value from the calculated mean value.

The differential corrections in equation 4 assume that compass and pitch are approximately equal to zero. If this is not the case, the differential corrections cannot be based on the aforementioned simplification. The general formulas for the differential corrections may be of interest for multibeam sonar mounted on submarines or similar vessels.

Finally, it can be stated that it is possible to use a multi-resolution grid based on, for example, a quadtree division ("quadtree tesselation"). Such grids are able to adapt to the density of the measured depth points.

Claims (6)

1. Procedure for field calibration of system parameters in multibeam sonar systems for seabed mapping, where a vessel carries multi-beam transducer equipment mounted under said vessel for swath coverage of many depth values in a plane mainly perpendicular to the vessel's direction of movement, with the emission of sonar energy at a number of different angles with the vertical and the reception of reflected sonar energy along the same angles, as a measuring strip is completed by driving straight ahead a given distance, overlapping measuring strips are formed by driving said vessel in opposite or crossing directions, and the measuring strips are transformed into a common reference system by resampling recorded reflection time measurement values to create grid models, where said system parameters are parameters in a parameter group that includes static corrections in roll (<f>), pitch (co), compass (k), the vessel's direction of movement (x), direction (y) in the horizontal plane perpendicular to said vessel's direction of movement, vertical depth direction (H), scaling factor (my) in the swath direction, angular scaling factor (rn*) and time delay (x), characterized by the difference between measured stripes and an approximation to a true seabed topography is minimized using the following equation: and X is a vector of unknown calibration parameters, Ai is a coefficient matrix, Pi is a weight matrix for the measured depth values, Li is a vector of the differences between a mean value H<*>(u) and Hi(u), where and n is the number of overlapping stripes in grid point u and Hj(u) is stripe in its depth value in u, l_2 is a vector of a priori parameter values, and P2 is a weight matrix for parameters to be calculated. 2. Method according to claim 1 characterized in that the gradients dHx, dHy to the relief of the seabed at a grid point u = (xp, yp) are calculated on the basis of the nearest neighboring grid points using the equations where a0o■ ■ a22 are coefficients calculated from equation (5) and equation (1) using the aforementioned neighboring grid points. 3. Method according to claim 1, characterized in that said difference is minimized by also taking into account a further system parameter which is a tidal correction (dH0j) which is individual for each measuring strip (i), the previously specified system parameters being common to all strips. 4. Method according to claim 1, characterized in that differential corrections to the system parameters are defined as follows: where v is the vessel's speed. 5. Method according to claim 4, characterized in that a credibility test is performed for system parameter errors determined by first defining where C° is a vector of simulated correction values for all system parameters k, and C2 is a vector of calibration values determined from two calculation rounds where equation (9) is used, where the first round gives a vector Ci with calibration values which is added to vector C°' to give Ci+ C°, after which a second round of calculation using equation (9) is based on the differential corrections, and then t is evaluated, with t < 0.1 indicating clearly credible and t > 0.9 indicating clearly unreliable. 6. Method according to claim 1, characterized in that an equalized depth value H(u) at grid point u is calculated using inverse distance weighting according to the equation where pi(u) is a weight calculated according to the inverse distance function which where s is the distance from the transducer to point u, thereby achieving block smoothing.