The Omega-K algorithm is widely used for creating imagery from synthetic aperture sonar and radar data. In early literature related to Omega-K beamforming, two alternative forms of the algorithm exist: one in which a critical matched-filtering stage is applied before Stolt-mapping and one in which it is applied after. The former was adopted as the standard approach in both fields. In the present study, it is demonstrated that applying the matched-filter prior to Stolt-mapping has the potential to alias acoustic energy associated with wide aperture angles, causing artifacts in imagery, whereas the alternative formulation avoids these artifacts.

The Omega-K algorithm, also known as the Range-Migration or Wavenumber algorithm, is a frequency-domain beamformer that is commonly used for focusing data captured by synthetic aperture sonar (sas) and synthetic aperture radar (SAR) systems. The principle behind the algorithm was first developed in the field of seismology1,2 and then adapted to the problem of focusing synthetic aperture radar data by Cafforio et al.3 It quickly found favor in the field of SAR and eventually sas due to its optimal point-spread function and broad-band compatibility.4–6 Today, despite the broader adoption of generalized time-domain beamforming procedures due to the proliferation of graphics processing units and high speed parallel processing,7,8 the broad-band limitations caused by the fixed integration angles of time-domain methods9 still make the Omega-K algorithm the beamformer of choice when considering ideal linear apertures.

The Omega-K algorithm converts between the domain represented by the two-dimensional Fourier transform of the raw data, i.e., the “Omega-K” domain, and the spatial spectrum of the complex image via Stolt-mapping,1 a change of variables from temporal frequency to spatial wavenumber that is usually implemented by interpolation. For signals in continuous time and space, the total process involves a two-dimensional Fourier transform, Stolt mapping, and an inverse two-dimensional Fourier transform. Application of the process to finite sized and discretely sampled data sets, however, necessitates an additional step in which a matched-filter is applied to correct for the range-migration locus at the origin, specified as the center of the imaging swath (see, e.g., the discussions in Bamler4 and chapter 4 of Hawkins10) How and when this matched filtering operation is performed is the subject of the present study.

In the earliest implementations of the algorithm described by Cafforio et al.,3 the matched filter was applied after Stolt-mapping. Bamler,4 on the other hand, proposed moving the matched-filtering step before the Stolt-mapping operation, where the relationship between the expression for the matched-filter and the system transfer function is intuitive. In fact, matched-filters of this form have found broader application in acoustics, for example, in the shadow-enhancement algorithm described by Groen et al.,11 “Quasi-holographic” imaging which preserves the focus of time-delayed, structural acoustic features and resonances,12,13 and rapid range-specific imaging used for performing autofocus.14 Descriptions of the Omega-K algorithm in primary textbooks related to synthetic aperture radar such as Carrara et al.,15 Soumeck,16 and Cumming and Wong,17 as well as several frequently cited theses related to synthetic aperture sonar such as Hawkins,10 Callow,18 and Cook,19 use the convention of Bamler.4 The adopted convention appears to be standard in acoustics with only one exception (de Paulis et al.20), which also happens to be co-authored by one of the original developers of the Omega-K algorithm for SAR.

It is worth noting that Bamler4 specifically assumes that the two variations are “equivalent from a systems theoretical point of view.” This seems to be the general assumption behind other papers discussing the algorithm and though some of the problems associated with the matched filtering process applied using the standard approach have been acknowledged (see, e.g., chapter 5 of Toleman21) to the authors' knowledge, no papers discuss any lack of equivalency between the two approaches. The present study demonstrates that the two variations, while equivalent in continuous time and space, are not equivalent upon digital implementation to sampled data captured over finite time and space. Furthermore, the variation which was more widely adopted in acoustic and radar literature can produce artifacts for near-range scatterers that are avoided using the original variation of the algorithm which places the matched-filtering step after the Stolt-mapping operation. Subsequently, for convenience “Pre-Stolt MF” will be used to describe the standard implementation in which matched filtering is applied before Stolt-mapping, and “Post-Stolt MF” will describe the case in which the matched filter is applied after Stolt-mapping.

Consider a synthetic aperture sonar moving at velocity v along a perfectly linear trajectory aligned with the x-axis of a Cartesian coordinate system. The sonar insonifies a set of isotropic point scatterers distributed on a plane some depth H below the sonar. Assuming a monostatic imaging system and a static position during echo reception, (“stop-and-hop”), the round trip propagation time between the sonar and a point located at (x, y, −H) on the imaging plane is
t p x , y = 2 c v τ x 2 + y 2 + H 2 ,
(1)
where c is the speed of sound and τ is “slow time,” i.e., the time the sonar has been traveling along the trajectory. The assumption of a linear aperture allows the imaging problem to be modeled in the slant range coordinate system (x, r),
t p x , r = 2 c u x 2 + r 2 ,
(2)
r = y 2 + H 2 ,
(3)
u = v τ .
(4)
Equation (2) defines the temporal locus, commonly referred to as a range-migration arc or target hyperbola, for the backscattered signal corresponding to an isotropic point scatterer located in the slant range at (x, r). The beginning assumption behind the Omega-K algorithm is that the reflectivity function for a scene can be described as a collection of point scatterers, each of which has a hyperbola embedded in the recovered data with time-delays described by Eq. (2), i.e.,
f u , t = x , r a x , r δ t 2 c u x 2 + r 2 d xdr ,
(5)
where a(x, r) is the complex scene reflectivity function modeled in the slant range and δ is the Delta function. In some formulations, δ(t) is replaced with a finite-bandwidth pulse with a defined center frequency.10 Note that the model in Eq. (5) is strictly describing the time-delays and complex amplitudes associated with the reflectivity function of the scene and terms related to absorption, spreading loss, directivity, pulse form, and bandwidth are all ignored. The function f(u, t) describes the relationship between a movement along the synthetic aperture and the delays associated with insonified scatterers during data collection, and the goal of the Omega-K algorithm is to invert Eq. (5) and recover the reflectivity function a(x, r). The Omega-K algorithm performs this inversion on F(ku, ω), the two-dimensional Fourier transform of Eq. (5), which is shown in Hawkins10 to be expressible as
F k u , ω = x , r a x , r e i r 4 ω / c 2 k u 2 + x k u d xdr ,
(6)
where a(x, r) is the spatial reflectivity function multiplied by the square-root of r. The transmission pulse term in Hawkins10 is ignored due to its replacement by the delta-function in Eq. (5) and, as in Gough and Hawkins,6 the amplitude functions are also ignored. Subsequently, a change of variables is performed via Stolt-mapping, i.e., the following substitutions are made resulting in an expression for F k x , k r,
k r 4 ω c 2 k u 2 ,
(7)
k x k u ,
(8)
F k x , k r = x , r a x , r e i x k x i r k r d xdr .
(9)
Equation (9) contains a two-dimensional Fourier transform, and it can be seen that the modified spatial reflectivity function a(x,r) can be recovered from F k x , k r by applying a two-dimensional inverse Fourier Transform. An additional step related to image origin definition is required due to the finite size of the data collection matrix. The effect of performing the Stolt-mapping process is to unify the range migration loci for each scatterer to that of the locus at the origin of the range axis which, to avoid excessive zero-padding, is generally defined as the center of the swath.4,10 To perform this origin conversion and create a focused image, the range-migration locus at the center of the swath must be deconvolved and shifted to the origin of the sampled signal matrix represented by f u , t. The deconvolution is generally performed by applying a matched-filter modelling range-migration at the center of the swath. The expression for this matched filter in the Omega-K domain is shown in both Bamler4 and Hawkins10 to be
β ω , k u = e i r 0 4 ω / c 2 k u 2 2 ( ω / c ) ,
(10)
where r0 is defined as the range to the center of the swath. Note that Eq. (10), which is a frequency-domain model of the impulse response of the measurement system at the swath center, has the same form as the various range-specific focusing filters described in the literature.11–14 Alternatively, by leveraging Eq. (8) and substituting the inverse form of Eq. (7) into the argument of Eq. (10), the matched filter can be redefined in the spatial wavenumber domain k x , k r,
β k x , k r = e i r 0 k r k r 2 + k x 2 ,
(11)
where
ω c k r 2 + k x 2 2 .
(12)
In practical implementations of the Omega-K algorithm, therefore, the options exist to focus data at the origin via Eq. (10) before performing the change of variables via Stolt-mapping, or to focus the data after Stolt-mapping by applying Eq. (11). In other words, the algorithm may use either of the two following formulations:
a x , r F 2 S F k u , ω e i r 0 4 ω / c 2 k u 2 2 ( ω / c ) ,
(13)
or
a x , r F 2 S F k u , ω e i r 0 k r k r 2 + k x 2 ,
(14)
where F 2 { } represents the two-dimensional inverse Fourier Transform and S { } represents the Stolt-mapping operation. Equation (13) represents the Pre-Stolt MF Omega-K variant, and Eq. (14) represents the Post-Stolt MF variant.

To illustrate the impact that the different algorithm variations have on image formation, data corresponding to a scan of a set of ideal points spaced from 5 to 50 m in 5-m range increments were simulated. The source is modeled as omnidirectional with a center frequency of 15 kHz and bandwidth of 20 kHz and is assumed to be co-located with the receiver. Aperture samples are spaced 1.5 cm apart along a linear aperture extending 50 m in length. Spreading loss has been applied, the data is cosine-tapered near the edges to prevent additional artifacts from abrupt signal discontinuities, and the Hilbert transform has been applied. The range migration arcs produced by the simulation and the result of applying the matched-filter defined in Eq. (11) are plotted in Fig. 1.

Figure 1(a) shows the range-migration arcs created by a set of simulated points and Fig. 1(b) shows the result of applying the matched filtering operation described in Eq. (10) to the curves in Fig. 1(a). From Fig. 1(b) it can be seen that the point at the swath center (t = 40 ms, r = 30 m) has been compressed to the point-spread function defined by its spectral coverage in the Omega-K domain, and the migration arcs to the left (i.e., earlier in time) have been inverted. Notable also are the wraparound phenomena corresponding with the close-range migration arcs: when the inverted migration arcs cross Time = 0, the left-most edge of Fig. 1(b), they reappear at the right-most edge and continue to propagate through the image in the reverse direction. Wrapping also occurs when range migration arcs cross either of the along-track boundaries visible in Fig. 1.

Related to imaging, the standard Omega-K implementation (“Pre-Stolt MF”) applies Stolt-mapping to the two-dimensional Discrete Fourier Transform (DFT) of the data shown in Fig. 1(b), whereas Post-Stolt MF beamforming applies Stolt-mapping to the two-dimensional DFT of the data shown in Fig. 1(a), which exhibits no wrapping artifacts due to being an unmodified version of the original raw data matrix. An additional point of consideration is that the range-migration arc inversion in Fig. 1(b) pushes signal energy closer to the T = 0 boundary. Figure 2 compares the images created via the two algorithm variations. Stolt-mapping was performed via spline interpolation to minimize interpolation artifacts.

In Fig. 2, the standard implementation shows a greater number and intensity of artifacts for ranges beyond 20 m. While both images exhibit some level of interpolation artifacts, which are especially visible in the targets at 5 m, the wrapped portion of the interpolation artifacts is significantly stronger in the image created via the standard implementation. The standard implementation also shows a set of artifacts between the ranges of 30 and 40 m that are not visible in Post-Stolt MF beamforming. At 40 m, the peak artifact level is −20 dB relative to the associated point scatterer level at that range in the Pre-Stolt MF variation, whereas using Post-Stolt MF the artifact level is −40 dB. Furthermore, Pre-Stolt MF shows a −6 dB loss in reconstructed peak amplitude for the near range scatterer relative to Post-Stolt MF. More differences can be demonstrated by examining the spatial spectra of individual scatterers. Figure 3 shows spectra evaluated by computing the two-dimensional DFT of a 2 × 2 m2 box around the point located at x = 0, r = 5.

A comparison of the wavenumber spectra in Figs. 3(a) and 3(b) reveals that Post-Stolt MF algorithm produces broader coverage in the wavenumber spectrum. This equates to higher resolution, a higher signal-to-noise ratio, and a wider effective synthetic aperture integration angle. Pre-Stolt MF effectively reduces the angular coverage and so may fail to reconstruct off-axis edge or structural-acoustic features. The causes of the reduction of wavenumber coverage in Fig. 3(a) are given in Sec. 5.

During the 2017 clutter experiments (“CLUTTEREX”) jointly funded by the Office of Naval Research (ONR) and the Strategic Environmental Research and Development Program (SERDP), low-frequency synthetic aperture sonar data were captured in the Gulf of Mexico using a rail-mounted synthetic aperture sonar. The sonar pinged every 0.025 m over a 40-m linear aperture, and the sonar had a bandwidth extending from 3 to 30 kHz. The altitude of the sonar above the sediment was approximately 2 m. Various unexploded ordnance (UXO) and clutter targets were arranged in a target field in front of the sonar, with the goal being to characterize the structural acoustic response of these targets as a function of angle. Figure 4 shows the results of imaging data from one scan captured using the rail system. The only difference between the images is the order of the matched-filtering step relative to the Stolt-mapping operation.

The Omega-K implementations used for creating the plots in Fig. 4 both apply approximately 20 m of zero padding to both sides of the range axis prior to beamforming. The purpose of this zero padding is to diminish the rapidity of phase oscillations in the ω-dimension and has the effect of reducing interpolation artifacts. It also somewhat reduces the severity of signal wraparound produced by the matched-filtering operation in the standard implementation; however, residual artifacts from this are still quite noticeable in Fig. 4(a).

The simulation and field experiments shown in this paper demonstrate that in a wide-beam, wideband imaging context the standard implementation of the Omega-K algorithm (“Pre-Stolt MF”), which places the matched-filtering step in the Omega-K domain, can produce excess artifacts in reconstructed imagery. The source of these artifacts is the reversal of the migration loci for near-range targets. The reversal of loci can push acoustic signals nearer to the swath boundaries creating rapid phase oscillations in the ω-domain that exacerbate interpolation artifacts. Additionally, the wide-aspect portions of the range-migration loci can alias around the range-axis, re-mapping the origins of the wrapped signal components, causing failure to focus during Stolt-mapping. Signals that alias around the range axis are no longer described by the range-loci structure defined in Eq. (2), which is leveraged by Stolt-mapping to unify the point-spread function across the sonar swath.

The artifacts produced by these effects are both identified in Fig. 2: the former is labeled in the figure annotations as “interpolation artifacts” and the latter as “MF wraparound” artifacts. As demonstrated in Fig. 3, the effect on the spatial spectra of near-range targets is to taper the coverage in wave-number space at wide aspect angles, causing an integration-angle limit and nullifying one of the innate advantages of Omega-K beamforming relative to time-domain beamforming.

The second variation of the algorithm (“Post-Stolt MF”), in which the matched filter is applied in the k x , k r domain, prevents these artifacts from forming by avoiding wrapping errors in the raw u, t domain. There does not appear to be any disadvantage from a computational standpoint to using this variation.

Though the formula of the matched-filter used in the standard implementation is intuitive and relates directly to various imaging techniques such as resonance focusing and shadow enhancement, from a practical beamforming standpoint, certain image artifacts related to signal wraparound and interpolation may be reduced or avoided by using the Post-Stolt MF beamforming algorithm for the aforementioned reasons.

The work performed in this study was sponsored by the Office of Naval Research under research Grant No. N00014-17-C-7033, and the Strategic Environmental Research and Development (SERDP) program project MR18-1051. The authors also gratefully acknowledge Steve Kargl, the PI of the CLUTTEREX experiments.

1.
R. H.
Stolt
, “
Migration by Fourier transform
,”
Geophysics
43
,
23
48
(
1978
).
2.
J.
Gazdag
and
P.
Sguazzero
, “
Migration of seismic data
,”
Proc. IEEE
72
(
10
),
1302
1315
(
1984
).
3.
C.
Cafforio
,
C.
Prati
, and
F.
Rocca
, “
SAR data focusing using seismic migration techniques
,”
IEEE Trans. Aerosp. Electron. Syst.
27
(
2)
,
194
207
(
1991
).
4.
R.
Bamler
, “
A comparison of range-Doppler and wavenumber domain SAR focusing algorithms
,”
IEEE Trans. Geosci. Remote Sens.
30
(
4
),
706
713
(
1992
).
5.
D.
Dendal
and
J. L.
Marchand
, “
Ω-k techniques advantages and weaker aspects
,” in
Proceedings IGARSS '92 International Geoscience and Remote Sensing Symposium
,
Houston, TX
(
May 26–29
,
1992
), pp.
366
368
.
6.
P. T.
Gough
and
D. W.
Hawkins
, “
Imaging algorithms for a strip-map synthetic aperture sonar: Minimizing the effects of aperture errors and aperture undersampling
,”
IEEE J. Oceanic Eng.
22
(
1
),
27
39
(
1997
).
7.
D. P.
Campbell
and
D. A.
Cook
, “
Using graphics processors to accelerate synthetic aperture sonar imaging via backpropagation
,” in
Proceedings of the High Performance Embedded Computing Workshop (HPEC'10)
,
Lexington, MA
(
September 15–16
,
2010
).
8.
I. D.
Gerg
,
D. C.
Brown
,
S. G.
Wagner
,
D.
Cook
,
B. N.
O'Donnell
,
T.
Benson
, and
T. C.
Montgomery
, “
GPU Acceleration for synthetic aperture sonar image reconstruction
,” in
Proceedings of IEEE Global Oceans 2020: Singapore–U.S. Gulf Coast
,
Virtual
(
October 5–14
,
2020
), pp.
1
9
.
9.
S. A. V.
Synnes
,
A. J.
Hunter
,
R. E.
Hansen
,
T. O.
Sæbø
,
H. J.
Callow
,
R. van
Vossen
, and
A.
Austeng
, “
Wideband synthetic aperture sonar backprojection with maximization of wave number domain support
,”
IEEE J. Oceanic Eng.
42
(
4
),
880
891
(
2017
).
10.
D. W.
Hawkins
, “
Synthetic aperture imaging algorithms: With application to wide bandwidth sonar
,” Ph.D. thesis,
University of Canterbury
,
Christchurch, New Zealand
(
1996
).
11.
J.
Groen
,
R. E.
Hansen
,
H. J.
Callow
,
J. C.
Sabel
, and
T. O.
Sabo
, “
Shadow enhancement in synthetic aperture sonar using fixed focusing
,”
IEEE J. Oceanic Eng.
34
(
3
),
269
284
(
2009
).
12.
K.
Baik
,
C.
Dudly
, and
P. L.
Marston
, “
Acoustic quasi-holographic images of scattering by vertical cylinders from one-dimensional bistatic scans
,”
J. Acoust. Soc. Am
130
(
6
),
3838
3851
(
2011
).
13.
D. J.
Zartman
,
D. S.
Plotnick
,
T. M.
Marston
, and
P. L.
Marston
, “
Quasi-holographic processing as an alternative to synthetic aperture sonar imaging
,”
Proc. Mtgs. Acoust
19
(
1
),
055011
(
2013
).
14.
T. M.
Marston
and
D. S.
Plotnick
, “
Semiparametric statistical stripmap synthetic aperture autofocusing
,”
IEEE Trans. Geosci. Remote Sens.
53
(
4
),
2086
2095
(
2015
).
15.
W. G.
Carrara
,
R. S.
Goodman
, and
R. M.
Majewski
,
Spotlight Synthetic Aperture Radar: Signal Processing Algorithms
(
Artech House
,
Norwood, MA
,
1995
).
16.
M.
Soumekh
,
Synthetic Aperture Radar Signal Processing with MATLAB Algorithms
(
Wiley
,
Hoboken, NJ
,
1999
).
17.
I. G.
Cumming
and
F. H.
Wong
,
Digital Processing of Synthetic Aperture Radar Data: Algorithms and Implementation
(
Artech House
,
Norwood, MA
,
2005
).
18.
H. J.
Callow
, “
Signal processing for synthetic aperture sonar image enhancement
,” Ph.D. thesis,
University of Canterbury
,
Christchurch, New Zealand
(
2003
).
19.
D. A.
Cook
, “
Synthetic aperture sonar motion estimation and compensation
,” Master's thesis,
Georgia Institute of Technology
,
Atlanta, GA
(
2007
).
20.
R.
De Paulis
,
C.
Prati
,
F.
Rocca
,
S.
Scirpoli
, and
S.
Tebaldini
, “
Focusing Synthetic Aperture Sonar (SAS) data with the Omega-K technique
,” in
Proceedings of the 2009 IEEE International Geoscience and Remote Sensing Symposium
,
Cape Town, South Africa
(
July 12–17
,
2009
), pp.
I-68
I-71
.
21.
M. A.
Tolman
, “
A detailed look at the omega-k algorithm for processing synthetic aperture radar data
,” Master's thesis,
Brigham Young University
,
Provo, UT
(
2008
).