Shallow water height mapping with interferometric synthetic aperture sonar

Shallow Water Height Mapping With Interferometric Synthetic Aperture Sonar Sérgio Rui Silva1, Sérgio Cunha1,2, Aníbal Matos1, Nuno Cruz1 1Porto University, Faculty of Engineering (www.fe.up.pt), 2CIMAR (www.cimar.org) Email: {srui, sergio, anibal, nacruz}@fe.up.pt AbstractHeight mapping of shallow water areas is an important task for many commercial and scientific applications like river navigability, infrastructure maintenance or natural resource monitoring. The use of an autonomous boat presents several advantages that ease the use of synthetic aperture images to create three-dimensional topographic maps through interferometric techniques. Sample data obtained during test trials illustrate how synthetic aperture can be used to generate imagery and bathymetry data. KeywordsHeight mapping, interferometry, synthetic aperture, sonar, autonomous boat I. INTRODUCTION This article describes an interferometric synthetic aperture sonar mounted on an autonomous boat that has been developed at the University of Porto ([4,5]). High quality synthetic aperture imagery is obtained through the use of this platform and its ability to gather high precision navigation data that is then used to synthesize the sonar images in a motion error free manner. The images are then processed to obtain the bottom topography maps of rivers, estuaries, dams, harbors, etc. Topography mapping of the river beds is of high importance for a large set of reasons. Very few river beds have been scanned with submeter accuracy. Furthermore, these are dynamic environments and periodic studies are valuable for river navigability assessment, for instance. Also, dynamic sand dunes that form in rivers and overgrowth of submerged plants have an impact in the ecosystem and so, should be identified and mapped. These procedures must be done in a regular basis and can be efficiently done with an interferometric synthetic aperture sonar. Surveillance of harbors is another important application for this system. Here, it is important to map objects like fallen containers, rocks and sand dunes that can be hazardous to ships using the docks. Using an autonomous boat presents several advantages as the sonar platform. Precise navigation is possible using a RTKGPS that is further enhanced with the data from an inertial navigation unit taking the navigation error down to the less than 1 centimeter. The synthetic aperture algorithms used (based on back-projection) employ this data to synthesize motion error free images. This enables low phase errors and the assembly of high quality interferograms. The height uncertainty is therefore controlled and the process of phase unwrapping for the height estimation is less complex. Furthermore the position of the boat and its velocity is controlled in a closed loop to follow predefined paths (this is possible with small surface crafts). As so, dual pass interferometry is possible with this sonar system, enabling the use of longer baselines, reducing the height estimation errors. As further advantage, the autonomous boat enables low operation costs and fast deployment, which is crucial in many surveillance applications. In comparison to real aperture or multi-beam system systems, this approach uses lower cost hardware to obtain the same level of performance. Figure 1: Autonomous with submerged transducer arrays. The autonomous boat carries a double array of transducers arranged in the vertical direction. They operate at 200kHz, with a beam width of 18º. After pulse compression and synthetic aperture image formation, a resolution of 1.5cm in range and 1cm in the along-track dimension is obtained. The theoretical height estimation error is within centimeter level. The nominal operation speed of the system is 0.5m/s, but can be extended to 1m/s. Results obtained during test missions in the Douro River, Portugal, are presented as examples of the system capabilities. This include height mapping through one-pass and dual-pass interferometry. This system enables gathering of both reflectivity and height maps. Combined together, they can enable better knowledge of the underwater territory demonstrating that an interferometric synthetic aperture sonar operated from the surface is a powerful tool for shallow water remote sensing. II. INTERFEORMETRIC HEIGHT MAPPING A synthetic aperture sonar system provides high resolution reflectivity maps of underwater areas by coherently combining echoes and forming a large virtual array. The obtained resolution is in the order of the dimension of the transducer used and more importantly is independent of the range ([2, 7]). Nevertheless this images provide only information in twodimensions, which is not sufficient for most underwater imaging applications. Figure 2: Transducer array. If the area observed is imaged by two synthetic apertures and a suitable looking angle, the information included in the phases can be used to obtain the third dimension and thus supply height estimation ([1, 3]). This is a know technique in synthetic aperture radar systems, but still a relatively unexplored concept in synthetic aperture sonar systems. The reason for this is linked to the difficulties that motion and medium phase errors induced in the obtention of good quality interferograms. Nevertheless there have been some very successful approaches to interferometric synthetic aperture sonar, as described in [13] and [14]. Interferometric height mapping is not restricted to be used with synthetic aperture ([15]). It is possible to create interferograms from real aperture sonar images. Nevertheless, using synthetic aperture has several advantages inherent to the resolution abilities of the synthetic aperture sonar, especially when there is the need to cover large areas as the resolution is independent of the range. Moreover interferometric synthetic aperture sonar enables coregistered imagery and bathymetry, which is of great interests for mapping applications. III. SYSTEM DESCRIPTION The described sonar system is built considering the operation in an autonomous boat as it platform (Figure 1). The boat is based on a modular architecture, comprising two pontoons, two independent thrusters for high maneuverability, batteries for energy, a CPU module with navigation sensors, a high speed communications link to shore for boat control and real-time data visualization from the sonar system itself. The navigation system is constituted by a RTK/DGPS, an inertial unit and digital compass unit. The autonomous boat is capable of following pre-programmed paths at a specified reference velocity ([4, 5, 6]). The sonar system is constituted by two vertically separated transducer arrays (Figure 2). Each individual array looking angle is adjustable and the separation between both arrays is 0.5m. This separation distance is a trade-off between fringe spacing and height resolution. Figure 3: Sonar system. Typically the two arrays are submerged at depths of 1m and 1.5m. The elements have 5cm of diameter and work at a frequency of 200kHz proving a maximum bandwidth of 40kHz and a beam-width of 18 degrees. The transducer array is dismountable from the boat for ease of transport, but has a mechanical support that enables precise and repetitive attachment. This mechanical support also guarantees a rigid bond to the boat while minimizing oscillation from the transducer array arm. The sonar is built around a FPGA system that generates the transmitting signals (chirps), controls the power and low-noise amplifiers, synchronizes the data to GPS time and manages the data recording thought the analog-to-digital channels ([6]). IV. AUTONOMOUS BOAT CONTROL The motion of the autonomous boat influences considerably the performance of the sonar system. To achieve the best results, the boat should move in straight lines at constant speed over ground and with minimal roll and pitch motions. Furthermore, it is desirable that the intended boat trajectories are described with minimal deviation to allow multiple sweeps in close parallel lines which are mandatory for dual passage interferometric height mapping. To accomplish these goals, a control system that automatically drives the vehicle along user specified trajectories at given speeds was developed. The control system is organized in two independent control loops, one for the velocity and the other for the horizontal position of the vehicle. The velocity loop determines the common mode actuation of the boat while the horizontal plane loop determines the differential mode actuation. These values are then combined to produce the commands for the starboard and port thrusters. The velocity loop is based on a proportional plus integral controller that assures that the velocity of the boat is, in steady state, the one defined by the user. The controller parameters where tuned to assure a smooth motion by rejecting high frequency noise from the navigation sensors. Different controller parameters can be used to obtain the best behavior at different velocity ranges. 00.05 0.1 0.15 0.2 0.25 0.3 0.35 0 10 20 30 40 50 60 70 80 cross track error [m] occurrences Figure 4: Line tracking performance. The horizontal plane loop implements a line tracking algorithm. It is composed by an outer controller that computes a heading reference based on the cross track error (distance of the boat to the desired straight line) and an inner controller that drives the vehicle heading to the given reference ([5]). This two stage control loop assures zero cross track error in steady state regardless of the water current and the tuning of the controller parameters took into account a dynamic model of the boat. Relevant to performance of the automatic control system of the boat, and in general to its behavior, is the quality of the real time estimates of the vehicle position, attitude and velocities. These are produced by a navigation system that relies on data provided by a magnetic compass, an inertial measurement unit, a pair of L1 GPS receivers and a L1+L2 RTK GPS receiver. When receiving DGPS corrections, the navigation system produces position estimates with accuracy better than 1.5cm for an update rate of 10 Hz. The whole system has already been tested in operational scenarios with great success. In particular, it has been observed that the vehicle describes intended trajectories with average cross tracking errors below 10cm. For illustration purposes, Figure 4 presents a histogram of the cross tracking error for a 50 m straight line described by the autonomous boat while collecting SAS data. V. INTERFEROMETRIC HEIGHT ESTIMATION PROBLEM Consider two synthetic aperture images a swath formed from two separate apertures. These images can be obtained from two different transducers mounted in the same system at different positions (single pass system) or from two different passages through the same swath (dual pass system). Phase Map 2D Phase Unwrapping Height Calc Height Estimation Unwrapped Phase Map Synthetic Aperture Processor Auto−Focus Navigation Solution Cross−track Compressed Data Along/Cross−track Compressed Data Estimate Rough Bottom Height Raw Echo Data 1 Pulse Compression Image Registration Synthetic Aperture Processor Auto−Focus Navigation Solution Cross−track Compressed Data Along/Cross−track Compressed Data Estimate Rough Bottom Height Reflectivity Map 2 Raw Echo Data 2 Pulse Compression Reflectivity Map 1 Interferogram Phase Filter Figure 5: Interferometric processing flow diagram. Each synthetic aperture image is processed independently following the normal steps for image generation (Figure 5): pulse compression; synthetic aperture processing using known navigation data and auto-focus. The synthetic aperture images are formed using a back-projection algorithm ([10]) that uses the available precise navigation data and a prior rough estimate of the bottom height to synthesize motion error free images ([4, 8, 9, 12]). The images are also submitted through an auto-focus step to reduce phase errors due to motion uncertainties and medium fluctuations. Prior to the creation of the interferogram, the images have to be carefully aligned through a registration step. This step identifies several control points in both images, computes the geometric displacement relationship between each pair of control points and finally warps through interpolation one image so it overlays the other ([17]). It is also useful to filter the phase data due to the specular nature of the reflectivity maps. The geometric model of the interferometric height estimation can be seen in Figure 6. The relation in the ranges to the same point given by the two images is obtained by: () 222 12 2 () 2cosRRRBBR ψβ +Δ = + − + (1.1) Where B is the separation arm between the receiving transducers, β is the angle this arm makes with the horizontal line that intersects the transducer T1 position and ψ is the angle of 1 R with also this horizontal line. R Δcan be given, in the single pass case, by the phase difference, φ , between the images: () 2 2 n R φπλ π + Δ= (1.2) β h R1 R2 1 Tψ z h0 y T 2 B ΔR Figure 6: Interferometric height mapping geometry. In the two pass case there is a factor of 2 in the phase that accounts for the fact that the images are formed by different transmitting transducer position. The height estimation is given by: 01 sin( )hh R ψ =− (1.3) In both cases, the obtained images can be modeled by the following equations: 1 1 2 2 4 11 4 22 2 (, ) , 2 (, ) , R j sb R j sb f fxy ss t Re ddt c f fxy ss t Re ddt c πλ πλ ττ ττ − − ⎧⎛⎞ =− ⎪⎜⎟ ⎪⎝⎠ ⎨⎛⎞ ⎪=− ⎜⎟ ⎪⎝⎠ ⎩ ∫∫ ∫∫ (1.4) Where ( , ) b s st τ is the pulse compressed echo data. The interferogram which gives origin to the phase differences is given by: 21 f fff φ =∠ ⋅ (1.5) In the particular case of this sonar system we can simplify the height estimation equation given two situations: single-pass and two-pass. In the single pass the transducers are arranged in a vertical manner and so 90º β =(Figure 7). In this case, we get: () 1 sin sin R RB B ψψ −Δ ⎛⎞ Δ≈ ⇒=⎜⎟ ⎝⎠ (1.6) This approximation is valid for range differences ( R Δ) much smaller than the range itself ( 1 R ) which is the case for smooth topographies. R h R1 R2 1 Tψ z h0 y T 2 B Δ Figure 7: Interferometric height mapping geometry (1-pass). For the dual passage case, then 0º β =(Figure 8): () 1 cos cos R RB B ψψ −Δ ⎛⎞ Δ≈− ⇒=− ⎜⎟ ⎝⎠ (1.7) Δ h R1 ψ z h0 y 2 R B T1T2 R Figure 8: Interferometric height mapping geometry (2-pass). To get an estimate of the height error we can differentiate the height function in regard to phase: 1dh dh d R dR dd d d φϕ φ ϕ Δ =⋅ ⋅ Δ (1.8) And get: 0cos( ) sin( ) 2 h R B φ ϕλ σσ ϕβ π =− + (1.9) Here, φ σ is the phase error which can be derived from the signal-to-noise ratio and 0 R is the medium range to the scene. From this equation we see that the longer the arm between the transducers, the lower the height estimation error. For the arm variation (which is an important factor in dual-passage operation), the height estimation error can be given by: 1dh dh d R dR dd dB d φϕ ϕ Δ =⋅ ⋅ Δ (1.10) Which in turn gives: 0cos( ) cot( ) hB R B ϕ σσ ϕβ =+ (1.11) So also in this case the longer the arm, the more insensible is the height error to the variation of the arm. 00.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.02 0.04 0.06 0.08 Phase error (deg) Height error (m) Height error vs Phase error 00.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0.5 1 1.5 2 2.5 3 Arm error (m) Height error (m) Height error vs Arm error B=0.5m (single-pass) B=1m (dual-pass) B=5m (dual-pass) B=0.5m (single-pass) B=1m (dual-pass) B=5m (dual-pass) Figure 9: Height error sensibility. Figure 9 show the height error estimate for increasing errors of phase and arm measurement error. It is possible to see that a longer arm produces lower height error for the same phase or arm error. The arm can’t be arbitrary long also due to decorrelation between the pair of synthetic aperture images. The maximum arm after which the images have an inappropriate degree of correlation is given by the following equation in the case of single-pass ([18]): 0 sin cv XT R B λ ρψ = (1.12) And for the case of dual-pass: 0 2 cos 2sin ch XT R B λ ψ ρψ = (1.13) XT ρ is the cross-track resolution. For the parameters of the sonar in consideration this results in 27 and 38m, respectively. Nevertheless, this will also produce a greater phase variation, which is problematic since the phase is a wrapped measurement. In equation (1.2), the phase difference φ is only known modulo 2 π . This means that the absolute phase ϕ is related to φ through a wrapping operator: () ( ) mod , 2W φϕ ϕ ππ π == +−(1.14) Therefore to obtain the height map estimate it is necessary to know the integer value of cycles (since the beginning of the scene) and retrieve the absolute phase. In the absence of noise, this is simply done by integrating the phase difference and adding 2 π if the difference of value at a discontinuity is smaller than π −or subtracting 2 π if the value is larger than π . However, in the presence of noise this phase integration will easily fail. Several algorithms have been derived to effectively solve this problem, such as path following methods or minimum norm methods ([16]). In this paper we will focus on a particular minimum norm method know as the least square phase unwrapping. If the adjacent phase differences do not differ by more than π ±, then the wrapped gradient of the measured wrapped phase differences, φ ∇, will equal the gradient of the absolute phase differences ϕ ∇, so that the solution that minimizes the following equation, () () 2 22 1, , , 00 2 12 2 ,1 , , 00 MN x i j ij ij ij MN y ij ij ij ij φφ φ φφ φ ε −− + == −− + == −−∇ + +−−∇= ∑∑ ∑∑ (1.15) is the minimum least squares solution to the phase unwrapping problem. Here, , x ij φ ∇is the discrete phase gradient in the xdirection and , y ij φ ∇the gradient in the y direction. Figure 10: Height map, DouroRiver (single-pass). Putting this equation in the form: ()() ()() 1, , 1, , 1 , , 1 , ,1, ,,1, 22 i j ij i j ij ij ij ij xy yy ij i j ij ij ij φφφ φφφρ φφ φφ ρ +−+− −− ⎧−+ + −+ = ⎪ ⎨∇−∇ +∇−∇ = ⎪ ⎩ (1.16) This equation can be identified as the discrete Poisson’s equation. After manipulation of this equation, it is possible to efficiently calculate a solution through the discrete cosine transform (DCT): { } ()() () , 1 2 cos / cos / 2 ij DCT DCT iM jN ρ ϕππ −⎧⎫ ⎪⎪ =⎨⎬ +− ⎪⎪ ⎩⎭ (1.17) VI. RESULTS The system was tested in the Douro River, Portugal. During these tests it was possible to obtain some experimental data. The test site has a muddy river bed with smooth height variations caused my accumulation of sand and mud. The results are still preliminary, due to difficulties in obtaining comparison data and the early development stage of the sonar hardware. The autonomous boat was programmed to make passages near the shore line and in close proximity to a harbor entrance. Figure 10 show the height map of an area near the shore. To test the dual-passage possibility the boat was programmed to make a parallel pass on the same area, with a 1m of distance between passes (Figure 11). This type of control would be very difficult to obtain if not for the use an autonomous boat. The arm variation is depicted in Figure 12. -2.709 -2.7085 -2.708 -2.7075 -2.707 -2.7065 -2.706 x 10 4 1.546 1.5461 1.5461 1.5462 1.5462 1.5463 x 10 5 Boat position (m, N) (m, E) 1st pass 2nd pass Figure 11: Trajectory described by the autonomous boat. The relative arm variation is small, but nevertheless it has to be accounted for during the height calculation stage. The resulting height map obtained with the two passes is presented in Figure 8. -10 -5 0 5 10 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Along-Track (m) Cross-Track (m) Arm Figure 12: Arm in the dual-passage mode. For this system an arm estimation error in the centimeter level will produce a height error on the same accuracy level. Figure 13: Height map, Douro River (dual-pass). At the entrance of the harbor there is a visible sand dune formed by accumulation of river sand (Figure 14). Figure 14: Height map of a harbor entrance in the Douro River. The coregistration of swath reflectivity imagery with bathymetric data, enabled by interferometric synthetic aperture sonar, renders the results much more readable. VII. CONCLUSION An approach to shallow water height mapping was presented in this paper. The use of synthetic aperture sonar imagery to produce interferometric height mapping presents several advantages, but also shortcomings. An autonomous, the control of which is made automatically presents ease of use and deployment in areas suitable for its operation, like harbors, rivers, dams, etc. This system can be used successfully in many applications. Dual pass interferometry is an important technique that might help reduce the height estimation errors through the use of a longer separation between images. This is specially the case when longer wavelengths are used. Nevertheless, care must be made to the extension of this separation due to the decreasing of fringe spacing, image decorrelation and specially motion errors. Also one must guarantee that the swath maintains its coherence during the time needed to obtain the echo data. ACKNOWLEDGMENT This project has been possible thanks to the contributions of several institutions. The first author is supported by the FCT through the Ph.D. scholarship SFRH/BD/19976/2004/5R38. The authors would like to emphasize the support given by CIIMAR in financing part of the described system. REFERENCES [1] Richards, M.A., "A Beginner's Guide to Interferometric SAR Concepts and Signal Processing [AESS Tutorial IV]," Aerospace and Electronic Systems Magazine, IEEE , vol.22, no.9, pp.5-29, Sept. 2007. [2] Silkaitis, J.M.; Douglas, B.L.; Hua Lee, "Synthetic-aperture sonar imaging: system analysis, image formation, and motion compensation," Signals, Systems and Computers, 1995. 1995 Conference Record of the Twenty-Ninth Asilomar Conference on , vol.1, no., pp.423-427 vol.1, 30 Oct-1 Nov 1995. [3] Griffiths, H.D.; Rafik, T.A.; Meng, Z.; Cowan, C.F.N.; Shafeeu, H.; Anthony, D.K., "Interferometric synthetic aperture sonar for high resolution 3-D mapping of the seabed," Radar, Sonar and Navigation, IEE Proceedings - , vol.144, no.2, pp.96-103, Apr 1997. [4] Silva, S.R.; Cunha, S.; Matos, A.; Cruz, N., "An Autonomous Boat Based Synthetic Aperture Sonar," Oceans 2007 , vol., no., pp.1-7, Sept. 29 2007-Oct. 4 2007. [5] N. Cruz, A. Matos, S. Cunha, S. Silva, "ZARCO - An Autonomous Craft for Underwater Surveys", Proceedings of the 7th Geomatic Week, Barcelona, Spain, Feb 2007. [6] S. Silva, S. Cunha, A. Matos, N. Cruz, “An In-SAS System For Shallow Water Surveying”, Proceedings of the 7th Geomatic Week, Barcelona, Spain, Feb 2007. [7] P. T. Gough, “Unified Framework for Modern Synthetic Aperture Imaging Algorithms”. The International Journal of Imaging Systems and Technology, Vol. 8, pp. 343-358. [8] Silkaitis, J.M.; Douglas, B.L.; Hua Lee, "Synthetic-aperture sonar imaging: system analysis, image formation, and motion compensation," Signals, Systems and Computers, 1995. 1995 Conference Record of the Twenty-Ninth Asilomar Conference on , vol.1, no., pp.423-427 vol.1, 30 Oct-1 Nov 1995. [9] Fortune, S.A.; Hayes, M.P.; Gough, P.T., "Statistical autofocus of synthetic aperture sonar images using image contrast optimisation," OCEANS, 2001. MTS/IEEE Conference and Exhibition , vol.1, no., pp.163-169 vol.1, 2001. [10] Shippey, G.; Banks, S.; Pihl, J., "SAS image reconstruction using Fast Polar Back Projection: comparisons with Fast Factored Back Projection and Fourier-domain imaging," Oceans 2005 - Europe , vol.1, no., pp. 96101 Vol. 1, 20-23 June 2005. [11] Tomiyasu, K., "Tutorial review of synthetic-aperture radar (SAR) with applications to imaging of the ocean surface," Proceedings of the IEEE , vol.66, no.5, pp. 563-583, May 1978. [12] Fornaro, G., "Trajectory deviations in airborne SAR: analysis and compensation ," Aerospace and Electronic Systems, IEEE Transactions on , vol.35, no.3, pp.997-1009, Jul 1999. [13] Chatillon, J.; Adams, A.E.; Lawlor, M.A.; Zakharia, M.E.; “SAMI: a low-frequency prototype for mapping and imaging of the seabed by means of synthetic aperture”, Oceanic Engineering, IEEE Journal of Volume 24, Issue 1, Jan. 1999 Page(s):4-15, Digital Object Identifier 10.1109/48.740152. [14] Pinto, M.A.; Fohanno, F.; “Interferometric transmission synthetic aperture sonar”, Electronic Engineering in Oceanography, 1997. Technology Transfer from Research to Industry., Seventh International Conference on 23-25 June 1997 Page(s):113 – 119. [15] D.E. Pryor; “Theory and test of bathymetric side scan sonar”,OCEANS '88. 'A Partnership of Marine Interests'. Proceedings 31 Oct.-2 Nov. 1988 Page(s):379 - 384 vol.2 Digital Object Identifier 10.1109/OCEANS.1988.23531. [16] Ghiglia, Dennis C.; Pritt, Mark D., “Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software”, ISBN: 978-0-47124935-1, May 1998. [17] Wen Yang; Chuanzhao Han; Hong Sun; Yongfeng Cao, “Registration of high resolution SAR and optical images based on multiple features”, Geoscience and Remote Sensing Symposium, 2005. IGARSS apos; 05. Proceedings. 2005 IEEE International, Volume 5, Issue , 25-29 July 2005 Page(s): 3542 – 3544. [18] Zebker, H.A.; Villasenor, J., "Decorrelation in interferometric radar echoes," Geoscience and Remote Sensing, IEEE Transactions on , vol.30, no.5, pp.950-959, Sep 1992.