Automatic feature detection and side scan sonar overlap navigation via sonar image matching
FIELD OF INVENTION 
     The invention relates to a method for matching sonar images via automatic object detection in order to implement side scan sonar overlap navigation. 
     BACKGROUND OF THE INVENTION 
     Sonar is a widely used method for determining the position of objects in the ocean. Sonar has been used for navigation and to find mines or other objects of interest on the ocean bottom. The location of objects on the ocean floor might cue an operator as to position of his vessel when he can compare the locations to a map of previously surveyed regions. Sonar scans may also reveal whether changes have occurred in the region. Also, one can develop symbolic feature maps for only data related to significant features identified by sonar. Unlike optical images, sonar images are very noisy. Thus, it is very difficult to identify objects, particularly small objects, in sonar images. Also, it is difficult to compare two sonar images to determine whether the images were taken at the same location. Some images have a very high information content while others are virtually flat or featureless. In some images subtle changes in texture are the only distinguishing feature. Another difficulty in analyzing sonar images results from the fact that the image is dependent upon the angle at which the sonar signal strikes an object. That angle, of course, is determined by the height of the sonar emitting vehicle from the ocean bottom and the distance the objects are from the source of the sonar signals. Objects which are close to the sonar source show very small shadows. Those shadows increase for objects as the distance from the sonar source increases. Thus, the sonar image of any particular object is dependent on both angle and distance from which the sonar signal is transmitted. 
     There is a need to compare sonar images to determine if they are both taken at the same location. This may be important for navigation or other purposes. After determining that two images had been taken at the same location, it may be important to know whether any changes have occurred in the position or number of objects present in one image as compared to a second image. Prior to the present invention, there has been no way to automatically compare the two images. Rather, trained, experienced technicians have made visual comparisons to determine if two prints or displays were taken of the same location. This process can take a great deal of time. Automatic detection systems developed for optical images do not work well for noisy sonar images. The noisy nature of sonar images often precludes the use of line and edge detection operators of the type used for video images. Pattern recognition techniques used on video images are not applicable to sonar images because the representation of objects on a sonar image is very dependent upon the distance and angle between the object and the source of the sonar. Thus, two images of the same object taken at a different distance and angle may appear to be quite different to the untrained eye. Connectivity analysis which has been used for video images cannot be used on noisy images such as those produced by sonar. Prior information about the scene is not available for sonar images so pattern recognition techniques cannot be used. Thus, there is a need for an automatic method of analyzing sonar images to identify objects in the image and to compare images for matching purposes. 
     An autonomous underwater vehicle requires an accurate navigational system for covert missions during which the vehicle cannot surface. Sonar has been used for obstacle avoidance while maneuvering but not for such navigation. This invention proposes side scan sonar for navigation. 
     Development of a side scan sonar overlap navigation system is complicated by the fact that objects which appear in the overlap region of one image may not be the same as objects appearing in the same overlap region on an adjacent pass. The discrepancy results from the fact that the return signal from bottom objects is direction and grazing angle dependent as well as the fact that sonar signals have a low signal-to-noise ratio which results in a somewhat random detection of objects. Nevertheless, objects in the overlap zone have a reasonable probability of being detected on both passes. There is a need for a method to correlate objects seen on one pass with objects seen on the next pass. Such a method should be able to tolerate missing objects. Such method would permit accurate navigation by overlapping and matching side scan sonar images. 
     SUMMARY OF THE INVENTION 
     We have developed a method for analyzing sonar images which involves the basic steps of normalizing the image, filtering the image and then displaying the normalized and filtered image. We may also apply a matching algorithm to the filtered and normalized image and a second image. The filtering process equalizes the contrast across the image. Because sonar images can be represented as raster images often containing 480 rows by 512 columns of pixels, one can digitize each pixel on a gray scale and then apply an algorithm to normalize the image. Then, one can convolve the image using filters. We prefer to use low pass filters having a generally round area. We further prefer to apply these filters in a sequence wherein the filter having the lowest area is applied first and subsequent filters having 2, 3 and 4 times the area of smallest low pass filter are used successively. We have found that the number of objects in an image will determine the best size filter to use. If the image has very few objects, one should use a smaller filter. If the image has a large number of objects, a larger filter is required. If one uses the wrong filter he may have either no objects detected or too many detections which leads to excessive matching computations. 
     We have found that four particular filters are useful in this process. Each filter will generate a certain number of false detections which we call false alarms depending on the conditions selected for normalizing the image. Through experience we have developed a chart of false alarm ranges which is shown in the figures. From this chart we can select the appropriate filter to be used for a particular image. In addition to the low pass filter, we also prefer to apply a spatial differentiation filter to combine the search for highlights and shadows in one step. This differentiation filter is a matrix comprised of one row and a number of columns which is equal to the number of columns in the low pass filter which has been applied to the image. After filtering objects in the image can be more reliably identified. Moreover, filtered images are well, suited for matching. We propose to match features detected twice within the common overlap regions of two adjacent sonar images and use state estimation to reduce navigation errors. By correcting both random and systematic errors in the navigation system, complete area coverage with no missed areas and without undue overlap is ensured. In particular, this capability is important for correcting drift errors associated with inertial navigations systems. Knowledge of navigation errors enables the vehicle&#39;s position and the positions of objects in the image to be estimated more accurately. 
     Other objects and advantages of the invention will become apparent as a description of certain present preferred embodiments of the invention proceeds. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a diagram showing the dependency of a sonar image upon the position of the vehicle generating the sonar signal. 
     FIG. 2 is a raw sonar image. 
     FIG. 3 is the image of FIG. 2 after being normalized in accordance with our method. 
     FIG. 4 is the image of FIG. 3 after it has been convolved using a present preferred low pass filter. 
     FIG. 5 is the image of FIG. 3 after it has been convolved with a second present preferred low pass filter. 
     FIG. 6 is the image of FIG. 5 processed to show highlight-shadow pairs in the image of FIG. 5. 
     FIG. 7 is a display of two sonar images taken at different times on which we have indicated the objects which match, the match being determined using the methods of the present invention. 
     FIG. 8 is a graph showing the rate of false detections for a 100,000 pixels image after being filtered. 
     FIG. 9 is a graph showing the rate of false detections of the same 100,000 pixels image after being spatially differentiated and filtered. 
     FIGS. 10 thru 13 are graphs of detection versus threshold for filters of various sizes. 
     FIG. 14 is a graph showing detection versus filter area for three different false alarm levels. 
     FIGS. 15 thru 18 are graphs showing how to select the constraints for threshold α. 
     FIG. 19 is a diagram based upon sonar scans showing a vehicle in the search mode without error correction. 
     FIG. 20 is a diagram similar to FIG. 19, but modified using state estimation to reduce errors in obtaining position estimates of the vehicle and features. 
     FIG. 21 is a block diagram showing the basic steps of one side scan sonar overlap navigation. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Referring to FIG. 1 we illustrate a sonar source 12 above an ocean bottom 10 having a series of objects 1 through 6 on the ocean floor. Objects 1, 3, 4 and 6 are mines of the same size. The vehicle 12 is at a known height h above the bottom floor. Objects 1 through 6 are at distances d 1  through d 6  respectively from directly below the vehicle. Assuming the vehicle emits sonar (sound, pressure) waves from its left side and then its right side, the sonar equipment will produce images. The left side image will have a scan width of s 1  and the right side image will have a scan width of s r . When the sonar wave strikes an object and is reflected back, the presence of that object will appear as a highlight and shadow in the image. The length of that shadow will depend not only on the size of the objects, but on the distance and angle at which the sonar wave is generated. In FIG. 1, lines 21, 22, 23, 24, 25 and 26 represent the shadows which would be cast by each object when insonified from the sonar source in vehicle 12. As can be seen from FIG. 1, those objects which are to the left of the vehicle will have shadows extending to the left of the object whereas objects to the right of the vehicle will have shadows extending to the right of the object. Objects 1, 3, 4 and 6 are identical in size and shape. However, the shadows cast are quite different in length depending upon the distance of the object from the sonar source. Furthermore, objects which are quite different in configuration or construction may cast an identical shadow. Shadows 22, 24 and 25 are of the same length, but result from the presence of three different objects. It should also be apparent from FIG. 1 that the relative orientation of the shadow will be in line with the path of the sonar wave which generates the shadow. The same object which has been insonified at two different angles may produce two different shadows. 
     In FIG. 2 we have shown an actual sonar image. The shadows in that image indicate the presence of an object. As one moves from right to left in the image of FIG. 2, the shadows become longer. The shadows which appear in FIG. 2 are usually reliable indicators of objects on the bottom. It is also known that mounds produce a highlight followed by a shadow while depressions produce the opposite effect of a shadow followed by a highlight. 
     A region of connected pixels of like gray level within a sonar image is called a &#34;blob&#34;. A detection of highlight blobs or shadow blobs is accomplished by growing regions of adjacent pixels whose gray levels exceed or fall below a given threshold value. We have developed two approaches to reliably detect highlight and shadow blobs. The first finds the highlight and shadow blobs individually and then uses an association algorithm. The second finds associated highlights and shadows in one step using smooth spatial differentiation. 
     A sonar image increases in contrast from near range to far range. In order to use the same threshold value at every point in the image, the image must first be normalized. Using our method the image of FIG. 2 can be normalized to the image of FIG. 3. A normalization procedure was chosen to divide the image into vertical strips such that the amount of contrast increases from one strip to the next by a constant ratio. The assumption is made that the ratio of the far ranges of the total image to the vehicle&#39;s height is 10.0, and that 8 intervals are sufficient. In the event that the common overlap region is not the total image, then the number of intervals in the partial image is less, say n, where n&lt;8. If the partial image is the total image, then n=8. The fractions of the image which determine the interval boundaries are given in equation 1. ##EQU1## Here it is assumed that the image is to the right of the vehicle&#39;s path and that the whole image is used. In general, these fractions must be adjusted depending on whether the image is to the left or right of the vehicle&#39;s path as well as whether the image is measured from the near or far edge. Once these fractions are computed, the image can be divided into intervals. Let the k th  interval go from y k-1  to y k . Let μ k  and σ k  be the the sample mean and standard deviation of the k th  strip. The normalization procedure forces the normalized image to have the relation between mean and standard deviation of a Rayleigh distribution by computing a new mean for the k th  strip by equation 2. 
     
         m.sub.k =σ.sub.k / 0.5227, for k=1 to n              (2) 
    
     Next compute the difference between the computed mean and the sample mean for each interval by equation 3. 
     
         Δμ.sub.k =mk-μ.sub.k, for k=1 to n             (3) 
    
     The following computations are required to generate the normalized image: 
     
         y.sub.k =(y.sub.k-1 +y.sub.k) / 2, for k=1 to n            (4) 
    
     
         y.sub.o =y.sub.o 
    
     
         y.sub.n+1 =y.sub.n 
    
     
         m.sub.o =m.sub.1 
    
     
         m.sub.n+1 =m.sub.n 
    
     
         Δμ.sub.o Δμ.sub.1 
    
     
         Δμ.sub.n+1 =Δμ.sub.n 
    
     The normalized sonar image is described by an intensity function IN(i, j), where i is the row index and j is the column index of the i, j th  pixel. IN(i,j) is computed as follows: 
     
         for y.sub.k-1 &lt;j&lt;y.sub.k, for k=1, . . . ,n 
    
     
         α.sub.j =(j-y.sub.k-1) / (y.sub.k -y.sub.k-1) 
    
     
         β.sub.j =1-α.sub.j 
    
     
         Δμ.sub.j =α.sub.j Δμ.sub.k +β.sub.j Δμ.sub.k-1 
    
     
         m.sub.j =α.sub.j m.sub.k +β.sub.j m.sub.k-1 
    
     
         IN(i,j)=max (0,min{255,50[(IO(i,j)+Δμ.sub.j)/m.sub.j ]}), 
    
     for all rows i where IO(i,j) is the image prior to normalization. 
     Both approaches used to associate highlights and shadows rely on a set of four low pass filter masks: ##EQU2## The corresponding areas are: 
     
         A.sub.1 =4, A.sub.2 =9, A.sub.3 =21 and A.sub.4 =37 
    
     Four different low pass images can be obtained by convolving each filter with the normalized image. FIGS. 3, 4 and 5 are examples of such images. The spatial differentiation filter size is a function of the low pass filter size. For example, if the image is to the left of the vehicle&#39;s path, then spatial differentiation is done by convolving 
     d 1  =[-1 0 0 0 0 0 0 0 1]  with the g 4  filter. The spatial differentiation filter for g 1  and g 2  has 5 elements, and the filter for g 3  and g 4  has 7 elements. The -1 goes with the shadow and the 1 with the highlight, so the array may have to be flipped depending on whether the image is to the left or right of the vehicle&#39;s path. 
     To better determine an automatic method for choosing the threshold level and smoothing filter size needed for target detection, detection studies were done on a flat, synthetic image with constant . Since this image is devoid of targets, any detections can be interpreted as false alarms. The number of highlight detections versus threshold for the various filters are shown in FIGS. 8 and 9. On a semilog plot, these curves approximate a straight line and may be modeled with 
     
         1n(FA)=K.sub.O +K.sub.1 α                            (5) 
    
     where FA is the number of false alarms per 100,000 pixels and α is the number of standard deviations that the threshold is displaced from the mean. K 0  and K 1  can be obtained by doing a least squares fit of the data. Equation 5 can be used to solve for FA in terms of α or to solve for α in terms of FA. 
     In order to determine the appropriate threshold and filter required to detect highlight and shadow blobs, methods had to be derived to determine how many blobs of a certain size are in the image. Because of the noisy nature of sonar images, this is a difficult task. Information in a sonar image is encoded by way of parameter variations in the probability density function of the returned signal. Therefore, the best one can do is obtain a noisy estimate of the number of blobs of a certain size. A detection versus threshold study for various filter sizes was done on an image of high information content after normalization. The results are graphically displayed in FIGS. 10 thru 13. FIGS. 10 thru 13 show the number of detections and false alarms for each of four filter sizes over a range of threshold values. FIG. 14 shows detections versus filter area for three different false alarm levels. This graph is the straight line fit of the false alarm data for the flat, synthetic image. On a semilog plot, these graphs suggest the, following: 
     1. Detections versus threshold for a high information content sonar image form approximately a straight line. 
     2. At high thresholds, the false alarm graph is below the real detections graph. 
     3. As the filter size increases, the false alarm graph gets closer to the real detections graph. 
     4. At a low enough threshold, the false alarm graph might cross the detections graph. 
     The log-log plot of FIG. 14 shows that the number of detections varies inversely with filter size. The graphs of FIGS. 10 thru 14 suggest some simple heuristics that are useful in developing an automatic detection algorithm: 
     1. Detections vary exponentially with threshold value. 
     2. Detections vary inversely with filter size. 
     3. The smallest filter (highest pass in frequency) yields an estimate on the total number of blobs present. 
     4. False alarm curves from the flat, synthetic image can help estimate information content at a given threshold. 
     An adaptive feature or target detection algorithm whose objective is to find suitable features for matching two different sonar images must solve two problems: 
     1. Find an appropriate filter size, 
     2. Find a threshold α that yields approximately the desired number of detections, N D . 
     The lowest frequency filter possible should be used to avoid too many detections. A large number of detections in both images can lead to false matches between subsets of the two different images. Also, a large set of points is undesirable for an initial match since the time required varies as N 4 , where N is the number of points. Finally it must solve the above two problems for images with both high information content and low information content. 
     The solution can be divided into three steps: 
     1. Find the lowest frequency filter size that yields a reasonable number of detections, 
     2. Find the range of usable information in the image, 
     3. Find a threshold that yields the desired number of detections. 
     The algorithm is easily presented using the following definitions and notation: 
     
         ______________________________________                                    
k         filter index. k = k1 is highest frequency,                      
          lowest area. k = k4 is lowest frequency,                        
          largest area                                                    
FA.sub.k (α)                                                        
          Solution of equation 5 for false alarms in                      
          terms of threshold α for k fixed.                         
α.sub.k (FA)                                                        
          Solution of equation 5 for threshold α in                 
          terms of false alarms for k fixed.                              
α.sub.FA (k)                                                        
          Solution of equation 5 for threshold α in                 
          terms of k for FA fixed.                                        
α.sub.max (k)                                                       
          Maximum threshold on a real image convolved                     
          with filter g.sub.k for which a detection will                  
          occur, obtained from a histogram.                               
D.sub.k (α)                                                         
          Number of detections on a real image                            
          convolved with filter g.sub.k threshold α,                
          obtained from region growing.                                   
D1.sub.k (FA) =                                                           
          D.sub.k (α.sub.k (FA)) Number of detections on a real     
          image, filtered with g.sub.k and with threshold                 
          set in terms of FA, obtained from region                        
          growing.                                                        
h.sub.k (α)                                                         
          Curve fit of D.sub.k (α) using two separate points        
          and the functional form,                                        
          h.sub.k (α) = exp (a.sub.k + b.sub.k · α)  
             (6)                                                          
N.sub.D   Desired number of detections, N.sub.D ≧ 8.               
N.sub.u   Usable detections.                                              
______________________________________                                    
 
    
     For purposes of dividing images into low information, medium information and high formation the following regions are defined: 
     
         ______________________________________                                    
Low information:                                                          
           4 ≦ D1.sub.k (2) &lt; 8                                    
                        k.sub.1 ≦ k ≦ k4                    
Medium information:                                                       
           4 ≦ D1.sub.k (2) ≦ 8                             
                        k = k4                                            
and                                                                       
           D1.sub.k (2) ≧ 8                                        
                        k = k1                                            
High information:                                                         
           D1.sub.k (2) ≧ 8                                        
                        k1 ≦ k ≦ k4                         
______________________________________                                    
 
    
     In the above expressions the false alarms have been set to two. 
     The function h k  (α), which gives the number of detections at a given threshold, is obtained in the following way: 
     
         ______________________________________                                    
Point 1:     D.sub.k (α.sub.max (k)) and α.sub.max (k)        
Point 2:     D1.sub.k (FA.sub.opt) and α.sub.k (FA.sub.opt)         
             If D1.sub.k (2) ≧ N.sub.u                             
             then: FA.sub.opt = 2                                         
             else: FA.sub.opt = maximum of FA subject                     
             to                                                           
              2 ≦ FA ≦ 4 and                                
              D1.sub.k (FA) ≦ N.sub.u                              
______________________________________                                    
 
    
     If D1 k  (FA opt ) is less than or equal to D k  (α max ), then h k  (α) is constant and equals D1 k  (FA opt ). Otherwise the two points above are used to define h k  (α) according to equation 6. 
     The algorithm defined below will find a suitable number of detections for high information content images without exceeding N D . For low information content images, it will adjust N D  downward so that an excessive number of false alarms are not included in the detections found by region growing. 
     1. First, find the maximum filter size that yields the minimum number of detections suitable for a match. Specifically, set FA =2 and k=k4. Search down in k until D1 k  (FA)≧4 or k&lt;k1. If k goes below k1 declare insufficient data and exit. Else save k and D=D1 k  (FA). 
     2. Determine if sufficient information exists to decrease the filter size to allow for more detections than just the minimum. Explicitly, if D1 k1  (FA)≧8 and D&lt;8 then search down in k until D1 k  (FA)≧8 and save k; otherwise, go directly to step 3. 
     3. Estimate the maximum usable detections subject to an upper bound, by searching up in FA until D1 k1  (FA)≧N D  or FA=4. Set N u  =Minimum {N D , D1 k1  (FA)}. 
     4. Bound α from below by determining which of a number of constraints is active. Evaluate the following, which are shown graphically in FIGS. 15 thru 18. 
     a) α A  =α k  (Nu) as shown in FIG. 15 
     b) If h k  (α) is not constant, then α B  is the solution to: N u  =h k  (α B ) which is shown in FIG. 16 
     As can be seen in FIG. 17 α C  is the solution to: 
     
         h.sub.k (α.sub.C)=FA.sub.k (α.sub.C) 
    
     
         α.sub.B =Maximum of {α.sub.B, α.sub.C } 
    
     c) Finally from FIG. 18 we see that if h k  (α) is constant then 
     
         α.sub.B =α.sub.k (h.sub.k (α)) 
    
       Consequently, α=Maximum {α A , Δ B  } 
     5. Use threshold α found in step 4 to obtain detections to be used for image matching. 
     Having thus selected a filter and threshold value, the sonar image of FIG. 2 can be normalized to produce the sonar image of FIG. 3. Then, that image is convolved to produce the sonar image of FIG. 4 or FIG. 5 depending upon the filter which is used. That image can be thresholded and region grown to form the image of FIG. 6 or used for matching. 
     In FIG. 7 we show two images which have been processed with our method and whose features have been matched by computer. In FIG. 7 lines between the images connect matching features which are labeled M0 thru M10 in each image. Our method using a computer program can quick match features which would be difficult or impossible for humans alone to do. Furthermore, feature-based matching is required for side scan sonar. Therefore, our image matching method is particularly useful for side scan overlap navigation as described below. 
     Although sonar images have not previously been used for navigation purposes the present method permits such usage. We assign coordinates or map positions to objects in the images. These coordinates can be used to determine the position of the vessel, to steer the vessel or for other navigation purposes. 
     A side scan sonar creates a sound reflectance image of the ocean bottom to the port and starboard side of the sonar bearing vehicle by active insonification of the water volume using a very narrow azimuth beam pattern. In search mode, the vehicle travels in a straight line and creates a sonar image whose length is proportional to the distance traveled and whose width is proportional to the range of the sonar. In order to do feature based navigation some man-made or natural objects must be detected in this first pass. The detected objects are assigned coordinates in the vehicle&#39;s navigation system. 
     The vehicle then makes a U-turn, and follows a translated path parallel to the previous path so as to create another sonar image with a fixed percentage of overlap with the previous image. This procedure is repeated with alternating left and right U-turns at the end of each pass until the desired area is completely covered. In the region of overlap, some of the objects detected in the first pass must be detected again. In each of these sightings, the position of each object will differ. Therefore, different coordinates will be assigned to these objects. Matching techniques are used to pair the same objects. 
     After one has determined by our matching techniques that the same object is in two images, he can check to see if that object has been recorded twice under different coordinates. If so, he can update the navigation system using information from our matching algorithms to reflect this fact. If these errors are not corrected, they will accumulate during the mission thereby making the navigation system less and less accurate. 
     FIG. 19 shows simulated results of side scan sonar overlap navigation using compass and RPM meter without error correction. The solid line 50 shows the vehicle&#39;s path in a search mode. The solid vehicle icon 51 shows the vehicle&#39;s true position at the time each feature is detected. The outline of the vehicle 52 shows the perceived vehicle location at the same time. Filled circles 53 indicate true positions of features; unfilled circles 54 indicate perceived feature positions. Arrows show sonar range. Features exist in the overlap region. 
     FIG. 20 shows simulated results of side scan sonar overlap navigation using compass and RPM meter with error correction. Objects are depicted and numbered as in FIG. 19. When a feature is detected a second time, state estimation is used to correct the perceived vehicle and perceived feature locations. The result is that by the time the smallest feature is detected twice, the corrections are such that the perceived vehicle locations are nearly identical to the actual vehicle locations. Likewise, the perceived feature locations nearly coincide with the actual feature locations. 
     In FIG. 21 we provide a block diagram of our matching technique. We begin with a first sonar image or set of images 60 and a second sonar image or set of images 62. Next we apply a matching algorithm to images 60 and 62 as indicated by block 64. Then we correct bias errors by state estimation as indicated by block 66. Finally, at block 68 we update positions of the vehicle and features. 
     Once the matching algorithm pairs the detections of one pass with the detections of an adjacent pass, this information is acted upon by the state estimation algorithm. The state estimation algorithm works with position data, the x and y coordinates of the centroids of the detected objects. As discussed above, the position estimates of the centroids of the same object as seen on separate passes will generally differ because of errors. With an appropriate mathematical model describing the statistical behavior of the errors, the state estimation algorithm uses the position data to refine the position estimates of the objects as well as to update the position estimate of the vehicle given by the navigation system. 
     The requirements for the matching of two side scan sonar images with known overlap are fairly simple because of the following assumptions: 
     1. Since the vehicle will be approximately the same height on each pass, the two adjacent passes will have the same scaling. 
     2. Readings from the navigation system can be used to correct orientation errors to within one degree. 
     3. In the deep ocean, the number of important features on each pass will be relatively small. 
     A procedure derived from &#34;Point Pattern Matching by Relaxation&#34; by Ranade and Rosenfeld was sufficient to do the initial matching. In addition to modifications of this approach, we have developed two additional matching algorithms. All three of these procedures must be used to find successful matches required to implement side scan sonar overlap navigation. Each of these procedures will be discussed in turn. 
     In the initial matching algorithm, navigational data is used to provide an upper bound on the expected error between two sets of feature points {P} and {Q} to be matched. If the distance between any two points P i  and Q j  exceeds the upper bound, then they cannot be matched and are assigned scores of zero. Use of this initial test eliminates impossible matches and greatly reduces the computational cost, which is proportional to the fourth power of the number of points. Those points which are possible matches are scored for every possible pair (P i , Q j ) in the following way. First set {P} is translatcd such that point P i  is placed on top of point Q j . Then for every point P h  in {P} and h≈i, the closest point Q k  (in {Q} with k≈j) to P h  is found and its distance squared, d 2 , is computed. The average errors in the points {P} and {Q} can be approximated with a scale factor r. With 
     
         x=d.sup.2 /r.sup.2                                         (7) 
    
     
         and ##EQU3## the score of pair (P.sub.i, Q.sub.j) is incremented by the corresponding Δs for each P.sub.h examined. A matrix [s.sub.ij ] is formed from the scores with s.sub.ij the score of pair (P.sub.i, Q.sub.j). The matching process is completed by choosing pairs whose score is both the row and column maximum. A data file of the matched pairs and their scores is then created for further processing. 
    
     In images containing many highlight/shadow blobs (regions of connected pixels of like gray level), it is possible to find subsets of points from one image that match subsets from the other image even though the matches are incorrect from a global point of view. A post-processing matching algorithm is used to eliminate these incorrect matches in the following way. 
     1. Rank all pairs according to their score starting with the highest score. 
     2. Compute an average of the vector translations defined by the pairs appearing in the ranked list. The first value of the average is the translation defined by the highest scored pair or first pair. 
     3. Find a new score for each pair based on the difference between their translation average and the average translation. Obtain a new score by using equations (7) and (8) to score the difference distance squared. 
     4. Accept the pair (or match) as good if its new score is greater than the threshold. 
     5. We prefer to update the average of the translation by the translation of a newly accepted pair. 
     The result of post matching is a data file of matched pairs of points all of which determine a common translation. Pairs of points which do not share this common translation are eliminated even if they had high scores from the initial matching. 
     We then prefer to use a third matching procedure confirm matches on sets of points too large to initially match in a reasonable amount of time. Input to this procedure consists of: 
     1. Data file of matched points {P} and {Q} from the initial match or from the postprocessing match. 
     2. Two data files of points {R} and {S} yet to be matched. 
     This procedure computes a match between sets {R} and {S} in which set {R} and set {P} are from one image and set {S} and set {Q} are from the other image. The method is as follows: 
     1. Compute the translation needed to align set {P} with the set {Q} using their centroids. 
     2. Apply this translation to set {R}. 
     3. Match {R} with {S} by finding points from {R} and {S} that are the closest to each other. Use equations (7) and (8) to score each matched pair. 
     The time required by this algorithm goes as N 2 , whereas the time required by the initial matching algorithm goes as N 4 , where N is the number of points. 
     After matching detections from adjacent passes, state estimation is used to reduce errors in both vehicle and feature positions. The state estimation procedure seeks a maximum likelihood estimate to the state vector μ of the system based upon the measurement vector z corrupted by the error vector e. The vehicle path geometry induces a model of the system, 
     
         z=f(u)+e                                                   (9) 
    
     with error statistics 
     
         E{e}=O                                                     (10) 
    
     
         and 
    
     
         E{ee.sup.T }=R,                                            (11) 
    
     where R is the covariance matrix. The well known maximum likelihood estimate is given by that u which minimizes the functional, 
     
         J(u)=[z-f(u)].sup.T R.sup.-1 [z-f(u)].                     (12) 
    
     This mildly nonlinear problem is quickly solved by the iterative procedure. 
     
         Δu.sub.k =H.sup.T R.sup.-1 H(z-f(u.sub.k))           (13) 
    
     
         u.sub.k+1 =u.sub.k +Δu.sub.k                         (14) 
    
     
         H=δf(u) /δu.sup.T                              (15) 
    
     The state variables which make up the vector u are: 
     xo j , yo j  the coordinates of the centroid of the j th  object, 
     xv i , yv i  the coordinates of the vehicle when its sonar &#34;sees&#34; object j, and 
     b k  the constant but unknown biases in the navigation data. 
     The measurements which make up the vector z are: 
     α j  the heading angle of the vehicle when its sonar &#34;sees&#34; object j 
     d j  the range of object j measured by the sonar, and 
     x i , y i  the change in position of the vehicle from point i-1 to point i as measured by the navigation system. 
     Typical equations making up the vector function f(u) are: ##EQU4## 
     
         d.sub.j =[(xo.sub.j -xv.sub.i).sup.2 +(yo-yv.sub.i).sup.2 ].sup.1/2 +ed.sub.j                                                 (17) 
    
     
         Δx.sub.i =xv.sub.i=xv.sub.i-1 +ab.sub.k +ex.sub.i    (18) 
    
     
         Δy.sub.i =yv.sub.i=yv.sub.i-1 +ab.sub.s +ey.sub.i    (19) 
    
     The convariance matrix R is a function of the type of navigation system, the geometry of the vehicle path, and the time interval. The choice of unknown biases represented by the terms ab k  and ab s  will also depend upon the type of navigation system used. 
     State estimation is able to reduce the errors in all of the state variables whenever the position coordinates (xo j , yo j ) of an object appear in more than one measurement equation. When an object is sighted twice, the graph connecting the position points of the vehicle and the object contains a loop. Such a loop contributes to redundancy, allowing errors to be averaged and thus reduced. 
     Our method reduces navigation errors during covert underwater operations during which the vehicle cannot surface to get a GPS fix. It enables overlap navigation with no missed areas and without undue overlap. 
     We are also able to identify unknown but constant bias errors enabling the vehicle&#39;s position and the positions of objects in the image to be estimated more accurately. 
     Although we have disclosed certain present preferred embodiments of the invention, the invention is not limited thereto, but may be variously embodied within the scope of the following claims.