Using Frame Theoretic Convolutional Gridding for Robust Synthetic Aperture Sonar Imaging

06/26/2017 ∙ by John McKay, et al. ∙ 0

Recent progress in synthetic aperture sonar (SAS) technology and processing has led to significant advances in underwater imaging, outperforming previously common approaches in both accuracy and efficiency. There are, however, inherent limitations to current SAS reconstruction methodology. In particular, popular and efficient Fourier domain SAS methods require a 2D interpolation which is often ill conditioned and inaccurate, inevitably reducing robustness with regard to speckle and inaccurate sound-speed estimation. To overcome these issues, we propose using the frame theoretic convolution gridding (FTCG) algorithm to handle the non-uniform Fourier data. FTCG extends upon non-uniform fast Fourier transform (NUFFT) algorithms by casting the NUFFT as an approximation problem given Fourier frame data. The FTCG has been show to yield improved accuracy at little more computational cost. Using simulated data, we outline how the FTCG can be used to enhance current SAS processing.



There are no comments yet.


page 6

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Synthetic aperture Sonar has proven to be one of the most important developments underwater imaging. Indeed, SAS can produce high resolution images of much larger areas than was previously possible [1]. Most reconstruction algorithms developed for Sonar fall into two categories: back-projection and Fourier domain (FD). Back-projection methods have shown promise in terms of resolution but requires significant computational effort. On the other hand, FD methods “can reconstruct images in real time on standard workstation computers” [2] but with varying quality [3].

Many different methods fall into the FD category, including wavenumber and range Doppler algorithms in SAS [4] and polar-format algorithms (PFAs) in SAR [5] [6] [7]. In each, a interpolation step emerges as the main issue limiting resolution. From a numerical analysis point of view, this problem is even more difficult since no one optimal mathematical framework exists for two dimensional interpolation [2]. Of the ways this interpolation step has been addressed, Strolt mapping [8]

, chirp z-transforms

[6], and non-uniform fast Fourier Transforms (NUFFTs) [9, 5] have had some success in SAS. Although chirp z-transforms seemlingly yield the fastest and most accurate algorithms, they are susceptible to irregular vehicle movements and therefore suitable in only the most ideal cases [6]. On the other hand, NUFFTs - namely fast Gaussian gridding (FGG) NUFFTs - have demonstrated a robustness to noise and jitter in the Fourier space in synthetic aperture radar applications while still retaining reasonable computational efficiency [10, 9, 7]. It is with this in mind that we present a similar but further refined algorithm for tomographic imaging: frame theoretic convolutional gridding (FTCG). The FTCG algorithm was developed in [11] where the authors demonstrate how by cleverly incorporating frame theoretic concepts into NUFFTs, high quality image reconstructions can be had at comparable computational efficiencies and with properties guaranteeing bounded error.

In this paper, we will look to briefly introduce SAS imaging processing, acquaint the reader to the FTCG algorithm in the one dimensional manner presented in [11], show how it can be adapted to two dimensions as to be applied to SAS, and then detail applications of this method to simulated data.

Ii Synthetic Aperture Sonar

SAS is as powerful as it is because it is able to overcome unavoidable hardware constraints. The along track resolution of Sonar, Radar, and other tomographic systems is proportional to the width of the aperture of the signal-emitting device [2]. Unfortunately, though, there is a trade-off between narrow apertures and device size: the smaller the aperture, the larger the necessary array of sensors. This factor has limited Sonar imaging in the past, but as illustrated in Figure 1, SAS processing has been able to work by utilizing a series of successive pings and combining them as to mimic a smaller aperture. Of the two methods we mentioned previously, note that in this setting back-projection approaches essentially perform several costly inverse Radon transforms and FD methods leverage the Fourier slice projection theorem (FSPT) to attempt to reconstruct scenes from frequency data quickly [12]. Herein lies the problem: the FSPT will rarely (if ever) yield data points in some sort of equispaced arrangement. This means that one must employ some sort of correction or use a flexible strategy like a NUFFT [12, 13]. If a correction (i.e. interpolation) is used, it may struggle with data collected under irregular vehicle movements and/or incorrect sound-speed estimation which can incur unsightly artifacts like blurring [14].

This all means that a well designed, robust NUFFT algorithm can be quite valuable for SAS imaging, let alone other tomographic fields. Some success has been had in synthetic aperture radar in applying NUFFTs to their similar FD algorithms and tend to offer a far greater robustness than interpolation-based approaches [12]. [4] provided an introductory view of NUFFTs in SAS, but this is not a well-studied problem.

Fig. 1: SAS imaging procedure with focus on Fourier domain methods for reconstruction.

Iii Image Reconstruction from Non-Uniform Fourier Data

Below we present a description of standard NUFFT and the FTCG algorithm that is used to enhance its performance.

Iii-a The NUFFT

The NUFFT essentially maps irregular Fourier data to a uniform grid by attempting to quantify the “contribution” of Fourier basis element to each data point. For example, if we consider an image as a (one-dimensional) reflectivity function from sound waves, then its NUFFT approximation is given by


where is assumed to be piece-wise smooth on , , are its Fourier coefficients on some raster , and is a truncation threshold determined by the support of . The calculation performed in (1) can be viewed as where is the th order Fourier partial sum approximation of . Here is a positive window function that is essentially compactly supported in the domain while also having the property that its Fourier coefficients, are essentially compactly supported. This is necessary to reduce the impact of aliasing while still maintaining computational efficiency. The support in the Fourier domain of determines the threshold in (1). Finally, should have an analytical expression so as not to require calculation on point values . Gaussian window functions have shown promise in SAR [10], though there are several other options.

Much work has been done in determining the quadrature weights employed in approximating the convolutional integral . In particular, iterative methods have been demonstrated to be effective [15].

For the purpose of our discussion, we will rewrite in condensed form. It is convenient to replace the sum limits with which does not affect the solution (only zeroes are added in) but makes comparison to the new method more straightforward. We write




and , which will be especially relevant later when we discuss frames.

As we have alluded to before, NUFFTs are known to have manageable computational cost. While the exact cost depends on the window function, the idea is that has small support in the Fourier domain. In this way, the convolutional expense is small (which is represented by the matrix ) so that can be calculated quickly and the subsequent reconstruction can be done by the weighted sums represented in (2).

Iii-B Frame Approximation

The frame theoretic FFT was designed in [16] to improve the accuracy of by incorporating the properties of frames into the design of the quadrature weights . A brief review of the method is provided below.

Definition III.1.

A frame is a set in a Hilbert space such that


for some

. Any such set of vectors that satisfy (

4) must span even if the set contains redundancies.

Observe that a frame is similar to a basis for but with a relaxation on linear independence.

Although the set , is typically not a basis if , for some sampling patterns forms a frame [17]. It was also demonstrated in [11] that even if this is not the case, the assumption that is a frame can still be quite useful.

The frame operator is defined as


where are known as the frame coefficients. Observe that, since is both linear and invertible [18], we can say


Typically the dual frame is not known in closed form, so an approximation is needed. Moreover, in applications we only know a finite number of frame coefficients, e.g.  , . Hence we are interested in the partial sum frame approximation


The method developed in [16] approximates the dual frame so that (7) converges to (6). The method is based on projecting onto a so called admissible frame for . The convergenge of (7) is then dependent on the properties of the admissible frame.

Definition III.2.

An admissible frame of is a frame chosen so that

  1. for some , , and all .

  2. for some , , and all .

Given an admissible frame, an approximation to the dual frame can then be written as


In [16] it was shown that the coefficients in (8) can be solved for using with


where is the Moore Penrose pseudo-inverse. From the discussion above, we can now approximate any as


To link these ideas to (2), we define , which was also used in [11] with respect to this data frame. Thus, we modify the finite frame approximation in (7) to our problem of recovering reflectivity function as111We note that (7) is not equivalent to (11), but we use the same variable so as to avoid cumbersome notation.


or, in matrix form


The frame approximation has several advantages over NUFFTs, not the least of which is better reconstruction as well as a known upper error bound for certain rasters, even when corrupted with noise [11]. Thus, the frame approximation provides a certain degree of confidence that other methods cannot. On the other hand, this comes with considerable cost; requires operations, in addition to the (off-line) calculation of , making it infeasible for many machines.

Iii-C Frame Theoretic Convolutional Gridding (FTCG)

The FTCG strategy bridges the gap between and that takes advantage of the efficiency of the NUFFT but provides higher accuracy by replacing the traditional vector of quadrature weights with a matrix of quadrature weights informed by frame theory. Specifically, the FTCG approximation is given by


where is an -banded matrix. The comparison to the traditional convolutional gridding method is easily seen in the matrix formula for in (2) and in particular in (3). Here , derived below, is banded as compared to diagonal .222As discussed in [11], for diagonal, the coefficients of are determined using standard numerical quadrature. In this case, the convergence of inherently depends upon the density of data points. Specifically, the density of points must be inversely proportional to the number of data points. However, in SAS imaging, this is not necessarily the case, as increasing can lengthen the spectrum. Hence in order to guarantee convergence, a different convolution approximation is needed.

We first describe how to design , the diagonal matrix of quadrature weights . Observe that the Frobenious norm difference between and is


where is a constant stemming from frame bounds, and and are defined in (12) and (2) respectively. Furthermore, we see that


which can be simplified to show


Therefore, solving


provides a straight-forward way to find optimal diagonal elements to solve (17). Let us now remove the constraint for to be diagonal, then we see that the minimum is achieved for


Inserting this into (2) yields coefficients


so that by (12) we end up with

! Therefore, we have a link between the two methods and we can formally designate a trade-off between computational efficiency and accuracy based on the band we impose for the inverted matrix. In other words, the coefficient matrix

in FTCG is such that


where is the Hadamard product and is a binary matrix with ones along a band and zeros elsewhere. This isolation of the elements of along a band provides the key aspect to the FTCG method: increasing improves the accuracy of the method while reducing its efficiency. In [11] it was suggested to choose so that the cost of the convolutional sums mirrors that of subsequent FFTs, making for a minimal increase in computation over but with the promise of a better reconstruction.

Fig. 2: error of 3 1-D methods over varied raster sizes.

Another important advantage of the FTCG over the standard CG method is improved convegence properties. As Figure 2 shows in a experiment over a jittered raster, simple convolutional gridding garners little benefit from more data. In fact, as was observed in [19], increasing the number of Fourier samples potentially reduces the accuracy since it amounts to increasing the size of the spectral domain. Hence the density between points actually grows. As discussed previously, the convergence of numerical quadrature methods is linked to the density of data points, i.e. , for . On the other hand, the convergence of the frame approximation is based on the truncation error, that is, the accuracy improves with increased . However, it requires computing the inverse of an full matrix. The FTCG provides a balance between the two. The accuracy increases with bandwidth , allowing practitioners with large datasets to adjust the band for both efficiency and accuracy.

Iv Two Dimensional FTCG

The authors of [20] present a frame approximation for non-uniform Fourier data points. The method is based on an extension of the one-dimensional admissible frame methodology into dimensions. We briefly describe the method below for .

Suppose we are given a raster of data where each element, , is such that for some


In sequel we will refer to defined similarly for some .

The data frames are defined as


As in the case, we define an admissible frame for .

Definition IV.1.

, , is an admissible frame with respect to in two dimensions given that the following are satisfied

  1. for , ,

  2. for , , and

For example,


satisfies the admissibility requirements.

Although the admissible frame construction yields accurate reconstruction for a variety of different data arrangements with bounded maximal error, it is computationally inefficient. Further, as suggested in (23) is not readily applicable to a NUFFT scheme. Hence we seek to adapt the method in [20] towards an accurate and efficient FTCG algorithm.

Fig. 3: Example Fourier data raster arrangements of noisy grids (left), asterisk patterns (center), and SAS-like set-ups (right).

As in the case, we use the convolutional gridding method for construction. In this regard, let us define as a bounded window function with essentially compact support over (that is, decays rapidly outside of that region) and be such that


For as defined in (23), we can show that is also an admissible frame for . To do so, consider the two criteria for admissibility:


using sufficient , from the admissibility of and bound for , and

again utilizing and from the admissibility of . Thus, having satisfied the criteria a) and b), we see that, indeed, is an admissible frame for our data frame.

The FTCG is analogous to the version. Assuming that our original function can be reconstructed from a separable Fourier space, we augment (2) to333For comparative purposes we write the sum to include all values . As before, when computing the sum, only a subset of those values are used since is compactly supported.


with truncation term

A more compact form for (25) is expressed as


where and each component defined as


The frame approximation with is given by


such that for we have


As in the case, the relationship between and can be understood by the optimization problem expressed by (17) with the analogous terms. In particular, we can construct the FTCG approximation of band as


where .

We note here that this construction assumes separability of the Fourier space. Future investigations will consider the non-separable case.

V Experiments

To show the effectiveness of the FTCG algorithm, we present the following results using synthetic data mimicking Sonar image capture. Future investigations will include real SAS data.

We look to show the viability of FTCG in Sonar by demonstrating its ability on a noisy grid, an asterisk pattern, and then a more challenging side-scan-SAS-like arrangement. Each of these arrangements were chosen to reflect an inherent robustness to data sampling patterns. Indeed, as demonstrated below, the FTCG is readily used for a variety of data environments.

For each case, we looked to reconstruct the function


which is similar to that used in [20].

V-a Noisy Grid Data Arrangement

Fig. 4: Function reconstruction from noisy data arrangement.
Fig. 5: Function reconstruction from asterisk data arrangement.
Fig. 6: Function reconstruction from SAS data arrangement.

In , it was proved in [11] that is an admissible frame for if the raster is such that

Building on a configuration that extends this optimal setting, as a first experiment we consider

In this case we used a band of 15 () for a matrix which we define as


of dimension meaning that we used approximately 1.7% of the information provided (and, therefore, 1.7% of the computations).

Figure 6 compares the frame and FTCG approximations for the example in (33). In terms of PSNR, the frame approximation was able to achieve a value of 39.79dB while FTCG gave a value of 25.91dB. It is important to note that no post processing was employed to reduce the effects of Fourier approximation on the the underlying function.444More precisely, it is the the Fourier partial sum divided by the window function, . Specifically, while the periodic function in (33) requires a significant number of Fourier coefficients to achieve a high resolution recovery, [21], for piecewise smooth functions, the error approximation should be measured against the Fourier reconstruction (divided by the window function ) of the underlying function.

This experiment serves as a baseline proof of concept that the FTCG can produce satisfactory reconstructions at substantially less cost. To test out more complex patterns that strain many of the frame theoretic assumptions, we present the next two experiments: the asterisk and SAS-like arrangements.

V-B Asterisk Data Arrangement

While one could think of the asterisk pattern given in Figure 3 as useful with regards to other tomographic imaging settings such as magnetic resonance imaging, a similar sampling could be found in circular SAS settings wherein a vehicle circles around a target, avoiding the shadowing present in side-scan [22]. For this reason, and the fact that it is vastly different than the noisy grid, we sought out how the frame approximation and our 2D FTCG algorithm would perform. We used with a band of 23 (), representing 10% of given in (34).

Both methods performed better in our asterisk pattern case than the noisy grid (of which we had sound analytical inferences towards accuracy). This is likely due to the increased density in the low frequency domain, since our test function can be suitably resolved. In particular, observe the formation of matrix given in (34) in Figure 7 (left). There is a high intensity region along the band (especially in the top left region) and decay off the center. Given how well our images appear, the information contained on the off banded areas has minimal information and is, therefore, not really needed. Figure 6 shows the reconstructions and their error. The frame approximation yielded a PSNR value of 40.06dB and FTCG produced a value of 37.48dB meaning that, despite the banded structure and limited information, FTCG performed only marginally worse.

Fig. 7: Intensity maps of the magnitude values matrices given by (34) used for the asterisk sampling pattern (left) and side-scan SAS arrangement (right) experiments.

V-C SAS Data Arrangement

In our last experiment, we show a case where FTCG out performs the frame approximation. Here we use a side-scan-like arrangement (similar to the polar format common in synthetic aperture radar as well [7]). The band was 15 which represented only of .

While both reconstructions were challenged, the FTCG method was able to provide a clearer picture of the original function. This fact is reflected by their PSNR values of 21.04 and 23.27 for the frame and FTCG approximations, respectively. The point-wise error did lend some merit to the frame approximation over FTCG though this may reveal some issue with that metric.

A key reason why the FTCG yielded a higher PSNR may be shown in Figure 7 (right). Observe the thin high intensity bands with large swaths of low magnitude information. The diagonal itself is of the highest intensity and our method captures this while disregarding the rest. In particular we note that the condition number while . Therefore we conclude that the FTCG serves as a pre-conditioner for the standard frame approximation. We will further pursue using the FTCG as a pre-conditioner in future investigations.

Vi Discussion

Understanding how to overcome the issues of non-uniform Fourier data is essential to improving SAS. While this approach to Sonar has dramatically improved underwater imaging, we see that most high-quality SAS reconstructions come from back-projection methods that are infeasible for computers without GPUs or other very fast processors. To achieve higher resolutions in real time by automated underwater vehicles, methods have to be developed that yield better and more consistent performance than the standard interpolations that typically cannot handle irregular motion and noise due to ill conditioning. Hence we have presented the frame approximation as well as its less expensive surrogate, the FTCG algorithm, as alternatives. Both methods demonstrated improved resolution properties on synthetic data arranged on several rasters without any preprocessing to correct for imperfections. Further work into adapting FTCG towards SAS should yield even better results with actual SAS data.


  • [1] Roy Edgar Hansen, “Introduction to synthetic aperture sonar,” InTech, 2011.
  • [2] Michael P Hayes and Peter T Gough, “Synthetic aperture sonar: a review of current status,” Oceanic Engineering, IEEE Journal of, vol. 34, no. 3, pp. 207–224, 2009.
  • [3] Charles V Jakowatz Jr and Neall Doren, “Comparison of polar formatting and back-projection algorithms for spotlight-mode sar image formation,” in Defense and Security Symposium. International Society for Optics and Photonics, 2006, pp. 62370H–62370H.
  • [4] Peter T Gough and Michael P Hayes, “Fast fourier techniques for sas imagery,” in Oceans 2005-Europe. IEEE, 2005, vol. 1, pp. 563–568.
  • [5] Bo Fan, Yuliang Qin, Peng You, and Hongqiang Wang, “An improved pfa with aperture accommodation for widefield spotlight sar imaging,” Geoscience and Remote Sensing Letters, IEEE, vol. 12, no. 1, pp. 3–7, 2015.
  • [6] Daiyin Zhu and Zhaoda Zhu, “Range resampling in the polar format algorithm for spotlight sar image formation using the chirp z-transform,” Signal Processing, IEEE Transactions on, vol. 55, no. 3, pp. 1011–1023, 2007.
  • [7] Bo Fan, Jiantao Wang, Yuliang Qin, Hongqiang Wang, and Huaitie Xiao, “Polar format algorithm based on fast gaussian grid non-uniform fast fourier transform for spotlight synthetic aperture radar imaging,” IET Radar, Sonar & Navigation, vol. 8, no. 5, pp. 513–524, 2014.
  • [8] Peter T Gough and David W Hawkins, “Unified framework for modern synthetic aperture imaging algorithms,” .
  • [9] Fredrik Andersson, Randolph Moses, and Frank Natterer, “Fast fourier methods for synthetic aperture radar imaging,” Aerospace and Electronic Systems, IEEE Transactions on, vol. 48, no. 1, 2012.
  • [10] Leslie Greengard and June-Yub Lee, “Accelerating the nonuniform fast fourier transform,” SIAM review, vol. 46, no. 3, pp. 443–454, 2004.
  • [11] Anne Gelb and Guohui Song, “A frame theoretic approach to the nonuniform fast fourier transform,” SIAM Journal on Numerical Analysis, vol. 52, no. 3, pp. 1222–1242, 2014.
  • [12] Mehrdad Soumekh, Synthetic aperture radar signal processing, New York: Wiley, 1999.
  • [13] BR Breed, “k-omega beamforming on non-equally spaced line arrays,” The Journal of the Acoustical Society of America, vol. 110, no. 2, pp. 1199–1202, 2001.
  • [14] D. A. Cook and D. C. Brown, “Analysis of phase error effects on stripmap sas,” IEEE Journal of Oceanic Engineering, vol. 34, no. 3, pp. 250–261, July 2009.
  • [15] James G Pipe and Padmanabhan Menon, “Sampling density compensation in mri: rationale and an iterative numerical solution,” Magnetic resonance in medicine, vol. 41, no. 1, pp. 179–186, 1999.
  • [16] Guohui Song and Anne Gelb, “Approximating the inverse frame operator from localized frames,” Applied and Computational Harmonic Analysis, vol. 35, no. 1, pp. 94 – 110, 2013.
  • [17] John J Benedetto, Alexander M Powell, and Hui-Chuan Wu, “Mri signal reconstruction by fourier frames on interleaving spirals,” in Biomedical Imaging, 2002. Proceedings. 2002 IEEE International Symposium on. IEEE, 2002, pp. 717–720.
  • [18] Ole Christensen, An introduction to frames and Riesz bases, vol. 7, Springer, 2003.
  • [19] Adityavikram Viswanathan, Anne Gelb, Douglas Cochran, and Rosemary Renaut, “On reconstruction from non-uniform spectral data,” Journal of Scientific Computing, vol. 45, no. 1-3, pp. 487–513, 2010.
  • [20] Guohui Song, Jacqueline Davis, and Anne Gelb, “A two-dimensional inverse frame operator approximation technique,” arXiv preprint arXiv:1511.02922, 2015.
  • [21] Jan S Hesthaven, Sigal Gottlieb, and David Gottlieb, Spectral methods for time-dependent problems, vol. 21, Cambridge University Press, 2007.
  • [22] HJ Callow, RE Hansen, S Synnes, and TO Sæbø, “Circular synthetic aperture sonar without a beacon,” in Proc. Underwater Acoustic Measurements (UAM), 2009.