Abstract:
A method of automatically performing a patch test for a sonar system is disclosed, where data from a plurality of overlapping 3D sonar scans of a surface as the platform is moved are used to compensate for biases in mounting the sonar system on the platform.

Description:
RELATED PATENTS AND APPLICATIONS 
       [0001]    U.S. Pat. No. 6,438,071, issued to Hansen, et al. on Aug. 20, 2002, and entitled “Method for producing a 3D image” and U.S. Patent application Ser. No. 11/504826 filed Aug. 15, 2006 by Hansen entitled “Method of constructing mathematical representations of objects from reflected sonar signals” are assigned to the assignee of the present invention and is incorporated herein by reference in their entirety including incorporated material. 
     
    
     FIELD OF THE INVENTION 
       [0002]    The field of the invention is the field of automatic calibration of three dimensional sonar imaging equipment. 
       OBJECTS OF THE INVENTION 
       [0003]    It is an object of the invention to produce an automatic method of measuring and correcting for bias in a three dimensional sonar imaging device. 
       SUMMARY OF THE INVENTION 
       [0004]    A three dimensional sonar imaging device is mounted on a device carrier or platform. Errors in mounting the device are automatically discovered and automatically compensated for by software as the device moves relative to a sonar reflecting surface. If sufficient information is not received to calculate compensation, the system automatically notifies the operator of the device carrier or the computer controlling an autonomous device carrier to move the device in an appropriate pattern to produce the necessary information. 
     
    
     DETAILED DESCRIPTION OF THE INVENTION 
       [0005]    Three dimensional sonar imaging systems use devices which are carried on device carriers or platforms such as ships, on lines attached to rigidly or loosely to ships, or on other mobile vehicles which may be tethered to a ship by a line which communicates power and information, or on free moving autonomous vehicles which operate under control of a computer system. Such systems monitor motions using sensors that translates the local sensor coordinates to the global coordinates based on bathymetric sonar. 
         [0006]    When the device is attached to the device carrier, it is difficult to mount the device accurately with respect to device carrier axis. For example, mounting the device on a ship to within a degree or two of horizontal is very difficult, and lining the device up with the directional heading of the ship in motion is worse. 
         [0007]    Without compensation, the device receiving sonar reflections from a flat, horizontal sea bottom, for example, would show the sea bottom with a few degrees of tilt. In U.S. patent application Ser. No. 11/504826, the data are fit to a plane and the slope of the seabed is returned. Moving the device carrier on to a different heading will show that the slope changes, which is impossible. 
         [0008]    Two such images can be statistically combined to give an offset or bias measurement of the roll and pitch of the installed device with respect to the device carrier. More complicated fitting of various images together will give higher order corrections to make the images converge. Such mounting biases must be detected and compensated for in a procedure called a “patch test” 
         [0009]    Manually carrying out such patch tests is very time consuming and prone to operator error, as the various biases are coupled and difficult to compensate for one at a time. The method of the invention uses a computer to compute the biases by analyzing a series of images produced as the device is moved, and mathematically computing the biases to minimize a bias factor. When the motion of the device does not provide enough images to calculate and correct for the biases, the system automatically directs the operator or controller of the device carrier to move in such a fashion that the correct images are efficiently received. The biases measured are automatically entered into the system, and used to recalculate the biases until the system converges upon biases that are of the same order as the system noise. 
         [0010]    This invention disclosure describes the use of 3D sonar data where a sonar or several sonars or echo-sounders collect simultaneous  3  sonar data. In order to use the 3D method one must be able to match two overlapping 3D data sets. This can be done by for example by the Iterative Closest Point (ICP) algorithm (presented in the patents and patent applications noted above and in references given therein) or by fitting regression planes to the data. The regression plane method will generally not provide data that can be used for estimation of linear offsets while the ICP will be useful for estimation of all offsets (linear and rotational). 
         [0011]    Note that we can use one sensor providing 3D data—a 3D sonar—or we can use multiple 2D sonars or echo sounders that together forms a 3D data set. If the environment is static, the 3D data set can be acquired over time. If the environment varies with time, the 3D data set must be acquired in a short time over which the data set can be assumed to be static—static in the sense that errors due to time variation are negligible. This 3D data set that may come from multiple descrete sensors is called the “sensor data set” in the following. They all refer to the same local frame of reference. 
         [0012]    Define vectors as follows:
       {right arrow over (p)} v : Vehicle coordinate vector   {right arrow over (p)} w : World coordinate vector   {right arrow over (p)} s : Sensor coordinate vector   S: Motion matrix—maps from sensor frame to vehicle frame. This is the offset matrix.   V: Motion matrix—maps from vehicle frame to world frame   P s : Sensor coordinate matrix   P v : Vehicle coordinate matrix   P w : World coordinate matrix
 
We have
       
 
         [0000]      P s ={{right arrow over (p)} s } 
         [0000]    and correspondingly for the vehicle and world systems. P s  is (4×N) where N is the number of points. 
         [0021]    For a given data set, the vehicle coordinates are given by: 
         [0000]      P v =SP s    (1) 
         [0000]    and the world coordinates are given by: 
         [0000]      P w =VP v =VSP s    (2) 
         [0022]    Assume two data sets—data set 1 global frame of reference and data set 2 global frame of reference. The two data sets cover the same area so that 
         [0000]      P w2,1 =WP w1    (3) 
         [0023]    where P w1  is data set 1 expressed in data set 1 global frame of reference and P w2,1  is data set 1 expressed in data set 2 global frame of reference. 
         [0024]    Likewise 
         [0000]        P   w1,2   =W   −1   P   w2    (4) 
         [0025]    where P w2  is data set 2 expressed in data set 2 global frame of reference and P w1,2  is data set 2 expressed in data set 1 global frame of reference. 
         [0026]    W is thus mapping the data from one frame of reference to the other and vice versa. 
         [0027]    Assume that the two data sets represent different geometries, so that the problem has a solution other than the trivial. Also assume the the alignment is not perfect, thus 
         [0000]      W≠I 
         [0028]    Where I is the identity matrix. 
         [0029]    W can be found in a number of ways: 
         [0030]    By using Fast ICP noted in the above referenced patents and patent applications. 
         [0031]    By using methods shown in the following references 
         [0032]    P.Besl and N.McKay. A method for registration of 3-D shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence. 14/2):239-256. February 1992 
         [0033]    S.Rusinkiewicz, O. Hall-Holt and M. Levoy. Realtime 3d model acquisition. In Proceedings of Siggraph (SIGGRAPH 2002), pages 438-446, 2002. 
         [0034]    By using multiple regression analysis noted in the above referenced patents and patent applications 
         [0035]    Let us assume that W has been found. 
         [0036]    We want to minimize 
         [0000]        E=P   w2,1   −P   w1    (5) 
         [0037]    as we should get the same coordinates from the two different frames of reference if all the coordinate transformations were correct—i.e. the offset matrix S is correct. 
         [0038]    Let us assume that the coordinate transformations are not correct, that V is known with high accuracy (i.e. we obtain position, pitch, roll and heading with high accuracy) and that S contains errors (i.e. pitch, roll, yaw and position offset errors). 
         [0039]    Combining (2), (3) and (5) yields 
         [0000]        E =( W− 1) V   1   SP   S1    (6) 
         [0040]    Where index 1 denotes data set number 1. Note that V varies from data set to data set while S is constant. This means that we assume that the sensor is fixed with respect to the vehicle. 
         [0041]    Let 
         [0000]        A =( W− 1) V   1   S    (7) 
         [0042]    In order to minimise E, we need to minimise A. To minimise A means to maximise A −1 . 
         [0043]    A may not be invertible, but its pseudoinverse A +  can be written 
         [0000]        A   + =( AA   H ) +   A   H    (8) 
         [0044]    Using the spectral theorem, we can write 
         [0000]    
       
         
           
             
               
                 
                   
                     AA 
                     H 
                   
                   = 
                   
                     
                       ∑ 
                       i 
                     
                      
                     
                       
                         λ 
                         i 
                       
                        
                       
                         
                           q 
                           → 
                         
                         i 
                       
                        
                       
                         
                           q 
                           → 
                         
                         i 
                         H 
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
         [0045]    As AA H  is Hermitian, its eigenvectors {right arrow over (q)} i  are orthonormal. 
         [0046]    If the relationship between data sets 1 and 2 is such that one eigenvectors {right arrow over (q)} n  dominates—it&#39;s eigenvalue is larger then the others 
         [0000]      λ α =max{λ i },   (10) 
         [0000]    then we need to correct that dominant degree of freedom, correct S correspondingly and continue the iterations until W≈I. The criterion will simply be the sum of eigenvalues: 
         [0000]    
       
         
           
             
               
                 
                   σ 
                   = 
                   
                     
                       
                         ∑ 
                         i 
                       
                        
                       
                         λ 
                         i 
                       
                     
                     &lt; 
                     k 
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
     
         [0047]    Where the limit k has to be according to the quality of V. If (11) is computed for every iteration, one will continue to iterate as long as the sum of eigenvalues decreases. When it converges, we have reached the noise level of the system. 
         [0048]    In each step of the iteration, correct each position with 
         [0000]    
       
         
           
             
               
                 
                   
                     d 
                     → 
                   
                   = 
                   
                     
                       1 
                       2 
                     
                      
                     
                       
                         λ 
                         α 
                       
                     
                      
                     
                       
                         q 
                         → 
                       
                       α 
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
         [0049]    where λ α  is defined by (10) and {right arrow over (q)} α  is the corresponding eigenvector. 
         [0050]    Using (2), this yields for data set 1: 
         [0000]        V   1   S   d   {right arrow over (p)}   s1i   ={right arrow over (d)}   (13) 
         [0051]    and for data set 2: 
         [0000]        V   2   S   d   {right arrow over (p)}   s2i   ={right arrow over (d)}   (14) 
         [0052]    where indexes 1 and 2 represent data sets 1 and 2 respectively, S d  is the correction of Sand {right arrow over (p)} s1i  and {right arrow over (p)} s2i  denote all the sensor coordinate vectors. 
         [0053]    From (13) and (14) respectively, we get 
         [0000]        S   d   =V   1   −1   {right arrow over (d)}{right arrow over (p)}   s1i   H ( {right arrow over (p)}   s1i   {right arrow over (p)}   s1i   H ) +   (15) 
         [0000]      and 
         [0000]        S   d   =V   2   −1   {right arrow over (d)}{right arrow over (p)}   s2i   H ({right arrow over (p)} s2i   {right arrow over (p)}   s2i   H ) +   (16) 
         [0054]    If data set 1 includes N sensor coordinate vectors and data set 2 includes M sensor coordinate vectors, there are (M+N) equations that define the 12 entries in the motion difference matrix S d . 
         [0055]    This can also be written 
         [0000]        S   d   =V   1   −1   {{right arrow over (d)}}P   s1   H ( P   s1   P   s1   H ) +   (17) 
         [0000]      and 
         [0000]        S   d   =V   2   −1   {{right arrow over (d)}}P   s2   H ( P   s2   P   s2   H ) +   (18) 
         [0000]    respectively. Here, {{right arrow over (d)}} is the 4 by N or M matrix that is formed when we use {right arrow over (d)} as column vector N or M times respectively. 
         [0056]    Use of the pseudoinverse of (P P H ) means that this represents the least square solution. 
         [0057]    The offset matrix S is found automatically as the device carrier is moving and producing overlapping images by the following algorithm. 
         [0058]    Based on overlapping data from two images, find A from (7) 
         [0059]    Compute AA H . 
         [0060]    Find the eigensets of AA H    
         [0061]    Find σ from (11). 
         [0062]    Compare σ with the results from previous iterations and see if it has converged. 
         [0063]    Find λ n =max{λ i } and compared to σ. If this ratio is above a factor (for example 0.9), set the correction vector 
         [0000]    
       
         
           
             
               d 
               → 
             
             = 
             
               
                 1 
                 2 
               
                
               
                 
                   λ 
                   n 
                 
               
                
               
                 
                   
                     q 
                     → 
                   
                   n 
                 
                 . 
               
             
           
         
       
     
         [0000]    Compute S d  from (17) and (18), ignoring the noise terms in the spectral decomposition of (P s2 P s2   H ) + . 
         [0064]    Compute the corrected offset matrix: 
         [0000]    
       
      
       S 
       c 
       =S−S 
       d  
      
     
         [0065]    Continue with item 1) again. 
         [0066]    With a 3D sonar this process can be automated and overlapping data sets will provide large data sets with high reliability. 
         [0067]    If there is insufficient data, the system automatically notifies the operator to move the device carrier to a location where sufficient data may be taken. 
         [0068]    In general the device carrier must carry the device over the same area several times and the area should have features so that errors in position of an object will be revealed. If the sea bed is completely flat, pitch and roll errors will be detected, but there is no information in the sea bed data that allows for detection of heading or linear offset errors. The automatic method described above will show this as the eigensets corresponding to these parameters will never dominate. Hence the method will provide input to how the vehicle should move in order to obtain a full set of parameters at a given accuracy. 
         [0069]    For a device carrier operated by an operator, the system sends a signal to proceed to a location where there is a sea bed feature and when such features are detected, the system directs the operator to traverse several paths over the feature on a different heading. The new set of images can then be automatically registered to give the heading and linear offset errors. 
         [0070]    In an autonomous vehicle the computer directing the vehicle is instructed to do the same thing. 
         [0071]    Obviously, many modifications and variations of the present invention are possible in light of the above teachings. It is therefore to be understood that, within the scope of the appended claims, the invention may be practiced otherwise than as specifically described.