Sparse Reconstruction-Based Inverse Scattering Imaging in a Shallow Water Environment

Bistatic sonar or multistatic sonar system can collect more scattering information of targets than a monostatic sonar system. In this paper, sparse learning via iterative minimization method (SLIM) is introduced to distinguish wave components for time-domain (TD) back-propagation (BP) inverse scattering imaging improvement. Unlike the prevailing high central frequency (>100 kHz) and wideband imaging sonar systems, a relatively low-frequency band (1–10 kHz) is considered here. Due to the low sidelobe output of SLIM, the investigated object’s surface in TD-BP image is much clearer in an ideal two-dimensional free field case. Furthermore, when the environmental information is known, this sparse reconstruction-based channel deconvolution method can be implemented to recognize, categorize the main propagating paths and then rectify their time of arrivals. Compared with the phase conjugation-based channel deconvolution method, the proposed approach’s results have fewer sidelobes and higher signal-to-background ratio in the simulation.


I. INTRODUCTION
Active acoustic imaging is the technique to locate an object's scattering strength distribution and characterize it. An acoustic color image or a tomography imaging result can be generated by actively scanning an object from various view directions or aspect angles.
An acoustic color image displays the angle-frequency spectrum of an object manifesting its scattering characteristics. Firstly, an acoustic color image is used to verify and validate the theoretic model with the experimental measurements. For example, a cylindrical shell's acoustic color has been studied, where the striations and structures in the image are associated with the elastic waves [1] or the internal structures [2]. Moreover, when the object is placed near a horizontal surface, the acoustic color images are considered, and the ray-based model [3]- [5], finite element model [6], and three-dimensional/hybrid models [7], [8] were developed to match their corresponding experimental configurations. The acoustic color image of the observed elastic waves can be generated by using fractional Fourier transform to separate each wave component [9]. After that, acoustic color images The associate editor coordinating the review of this manuscript and approving it for publication was Yunjie Yang . are also used for object classification and target identification [10], [11] because objects of different shapes or materials exhibit different scattering responses. Aside from object identification, the acoustic color images will help to determine the aspect angle of a cylinder [12].
The acoustic tomographic imaging technique often densely collects scattering information and maps the time series data to spatial scattering strength distribution. We treat the synthetic aperture sonar (SAS) imaging as one type of acoustic tomographic imaging for similar data acquisition and processing methods. An acoustic tomographic imaging system often adopts a high center frequency (>10 kHz) and a broadband/wideband transmitted pulse. The back-propagation (BP) [13] method and Fourier reconstruction method are used in a reflection tomographic imaging case [14], to reveal the geometric shape and highlight points. Acoustic tomographic imaging, or SAS imaging, can make the multiple scattering phenomena more apparent in the spatial domain [7], [15], [16]. Acoustic quasi-holographic images better reveal each wave component's features, such as propagation time and occurring position, as well as multiple scattering between objects [17], by back-propagating the time series data to a fixed plane then display the image along time and cross-range axes. Other than the mentioned methods, VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the linear sampling method (LSM) [18] is also applied in the underwater object imaging task with the partial frequency variation approach [19]. LSM, unlike acoustic tomographic methods, calculates the indicator function to mark whether a pixel/voxel is within the scatterer. The literature listed above does not involve a multipath/ multi-mode propagation effect in an underwater waveguide. It needs to deconvolve the channel impulse response (CIR) first. The acoustic color assisted by the holographic technique shows the differences among a rigid sphere, a soft sphere, and a directional source [10]. Target imaging is also numerically performed in a range-depth waveguide [20], [21] through DORT (decomposition of the time-reversal operator) and BP.
Compressive sensing or sparse reconstruction method is widely used in underwater acoustics [22], especially in the field such as direction of arrival (DOA) estimation [23]- [25], CIR estimation [26], [27], and near-field acoustic holography [28], [29]. However, seldom underwater acoustic imaging works are under the compressive sensing framework. The alternating direction method of multipliers (ADMM) has been used for distributed optimization to improve the resolution of SAS images. Similarly, ADMM is applied to solve multi-frequency far-field equations simultaneously in LSM imaging [19].
The main contributions of this paper include: (1) Through the classical back-propagation inverse scattering processing (ISCP) method, a time-domain (TD) back-propagation (BP) scheme is used here to avoid arduous work at each frequency. Furthermore, the TD-BP combined with the SLIM method outputs a better image in two-dimensional free-field space. Also, the imaging process and its corresponding physics are connected. (2) Analyzes the impact of waveguide environmental parameters on defocusing and ghosting by the conventional ISCP method. A sparse reconstruction-based channel response deconvolution method is proposed to solve these adverse effects in the underwater waveguide. This paper is organized as follows: Section II reviews the inverse scattering problem in a two-dimensional ideally free field case. The time-domain back-propagation method using the output of a matched filter (MF) or SLIM is presented and discussed. Section III studies the inverse scattering problem in a range-independent shallow water waveguide, where the multipath effect brings defocusing and ghosting into the image generated by the conventional BP algorithm. This sparse reconstruction-based channel deconvolution imaging is proposed to solve that problem. Section IV gives the conclusion.

II. ISCP IMAGING IN A TWO-DIMENSIONAL (2D) FREE FIELD SCENARIO A. A GENERAL INVERSE SCATTERING IMAGING METHOD
The inverse scattering procedure using an acoustic or electromagnetic bistatic system can be modeled as Fig.1. In Fig.1, the center of a scattering object is set as the origin O. Then, the source's location is at r inc with its distance to the origin O being R s = ||r inc ||. D is the scattering object with ∂D being its surface. The incident direction is n inc = −(cos φ inc , sin φ inc ). Similarly, when the receiver moves at r sct , its distance to the origin point is R r = ||r sct || and the scattering direction is n sct = −(cos φ sct , sin φ sct ). r t is a point in the imaging area.
A classical inverse scattering method is the backpropagation (BP) method. The BP method needs first to convert the scattered pressure data P(r sct ; r inc , ω) into a far-field pattern f (n sct , n inc , ω).
The far-field pattern f is dependent on the scattering direction n sct and the incident direction n inc . The generalized Bojarski transform builds the relation between the far-field pattern and the shape of a scattering object through physical optical (PO) approximation [13].
B g (n sct , n inc ) = f (n sct , n inc ) + f * (−n sct , −n inc ) where indicator function (r t ) only takes 1 when r t is within the scatterer. The final ∓ in (1) is ''+'' for the case that the object has the rigid/Neumann boundary condition. It clearly shows B g is the spatial Fourier transform of −∇ 2 (r t ). Therefore, the BP imaging [13] I BP (r t ) = ± k 2 8π 2 inc sct B g (n sct , n inc ) × e ik(n sct −n inc )·r t d sct d inc (2) is to find the outer and inner bounds of a scatterer's boundary.
In (2), inc is the aperture composed by different source positions, and sct is the aperture of receivers. Usually, (2) can be relaxed to where r t is a pixel in the imaging area. The disadvantage of the BP method is that it has to process data of one frequency at each time.
In practical applications, the probing signal is a wideband waveform, and the imaging result can be directly gotten from the received time-domain signal p(t) without the far-field pattern conversion: where δ(·) is the one dimensional Dirichlet function, c is the propagation speed in the medium. By reviewing Fig.1, it can be found that (4) uses the linearized approximation of ||r t − r s || ≈ R s + n inc · r t , and ||r − r t || ≈ R r − n sct · r t . This indicates (4) is originated from the time domain delay-andsum expression: Ideally, (5) should also produce peaks indicating the outer bound and the inner bound of an object's surface as (2). Hereafter, this paper uses BP's time-domain equivalence (5) to output the inverse scattering imaging results.

B. SOUND SCATTERING AND INVERSE SCATTERING IMAGING IN A 2D FREE FIELD
Here, we study the inverse scattering imaging of an infinite long rigid cylinder in a 2D free field. The scattering response S(ω) of an infinite long rigid cylinder with radius a is well studied, that where J n and H (1) n are the nth order of Bessel function and Hankel function of the first kind respectively. In the time domain, S(ω) corresponds to a specular reflection wave and creeping waves through the Sommerfeld-Watson transformation [30]. Their propagation paths are depicted as Fig.2 and their propagation time t sr (specular reflection) and t n (for the nth creeping wave) are: where is the creeping wave's propagation speed on the cylinder's surface. In (8), n ∈ Z, and n ≥ 0 represents the anticlockwise propagation direction while n < 0 is the clockwise direction. The strength of these two waves can be found in [30], and not discussed here.

C. SIMULATION IN 2D FREE FIELD SPACE
To illustrate how the inverse scattering imaging works, a simulation is built as Fig.3. In Fig.3, a 0.5 m-radius rigid cylinder is located at (x = −1 m, y = 1 m), and the source stays at (x = 20 m, y = 20 √ 3 m). The source emits a 10 ms long, 3-9 kHz up-sweep Tukey windowed linear frequency modulation (LFM) pulse, while a receiver moves circularly on a 40 m-radius circle (R = 40 m). The sound speed is 1500 m/s. The pressure field is generated by the µ-diff toolbox [31].
The creeping waves are more likely to be observed at the forward scattering zone that |φ sct − φ inc | < 90 • . Thus, Fig.4 depicts the envelopes of matched filtered received data when the receiver is at the interval of 150 • to 330 • . The simulation data perfectly matches the time of arrivals (TOAs) of specular  reflection (black line) and the creeping wave (green line). Besides, when |φ sct − φ inc | < 5 • , the scattering amplitude is far higher than other receiving angles, which displays as the brightest spot in Fig.4. The next part is to transform the received time-series data into the ISCP image. Fig.5 is the superposition of transformed images from several selected receiving angles from Fig.4. Each bright ambiguity line in Fig.5 is tangent to the cylinder's surface (the white dashed circle), which is the same as the conventional mono/bi-static sonar's image. In Fig.5, all noticeable parts are marked by white numbers. ''1'' represents the lines formed by reflection waves when the receiver is at the 150 • -200 • interval, while ''2'' corresponds to the case when the receiver is within 270 • -330 • . These nearly straight lines have neighboring points of tangency and forming a superposition region numbered as ''4''. ''5'' is the counterpart for ''2''. ''3'' an ''X''-shaped pattern, is the intersection of reflection waves' (1 and 2 in Fig.5) from selected observation angles. The creeping waves have a similar result. ''6'' is the belt formed by creeping waves accompanied by ''1''. The reason why ''6'' does not appear significant ramifications as ''1'' in Fig.5 is that the TOAs of creeping waves in Fig.4 are linearly dependent on |φ sct − φ inc | as in (8), while the reflection waves' TOAs involves the sine of |φ sct − φ inc |. The creeping waves also generate a weaker ''X''-pattern (''8'') and have two points of the tangent (''9'' and ''10'').

D. SPARSE RECONSTRUCTION BASED IMAGING METHOD
By reviewing Fig.4 and Fig.7(b), it can be found that: (1) The output of MF has a wide main lobe, which may not well separate the reflection and creeping wave (because a windowed 3-9 kHz LFM pulse, the time resolution of MF is larger than 0.167 ms); (2) TD-BP using the output of MF produces lots of unwanted ripples and sidelobes. The sparse learning via iterative minimization (SLIM) [32] algorithm is used for a higher time resolution and imaging quality enhancement, as in Fig.8, when processing the same data as Fig.4. In Fig.8, the sidelobes near the receiving angle 246 • are eliminated, and the wave structures are clearer.
It should be noticed that, when using SLIM, an implicit assumption is made that the frequency dispersion is neglectable. Here, for the 3-9 kHz bandwidth in use, according to (9), c l ranges from 1181 m/s to 1327 m/s. On the other hand, in two-dimensional free-field space, the cylindrical spread wave propagation causes a √ kR attenuation, which will reduce the impact of creeping waves' dispersion. Also, the propagation distance that waves go through will not make a big difference. Fig.9 exhibits the BP imaging output using SLIM processed data. The overall structure in Fig.9(a) is almost the same as Fig.7(a), while their zoom-in parts exist some differences. The inner and outer boundaries ''4'' and ''11'' are narrowed for the high time-delay resolution offered by SLIM. Besides, SLIM eliminates the number of sidelobes in MF output as well as lower their sidelobe levels, making the inner spot ''12'' and ripples in Fig.7(b) disappear. One can see that the sidelobe level is at least 2dB lower in Fig.8. The ambiguity belt is also alleviated in Fig.9.
Nevertheless, SLIM boosts the time-delay resolution at the cost of neglecting low amplitude components. Compared with Fig.7(a), the length of the inner/outer boundary is shorted. This is caused by some low energy peaks dropped by SLIM.

III. ISCP IMAGING IN A THREE-DIMENSIONAL (3D) RANGE-INDEPENDENT WAVEGUIDE A. SCATTERING AND INVERSE SCATTERING IMAGING
In a shallow water waveguide, we use the Ingenito model [33] to study the sound scattering by simple shape targets. Following the Ingenito model, when given a rigid sphere's scattering function S in free-field space and the coordinates of the source(r inc = (r s , φ s , z s )), target(r t = (0, 0, z t )), and receiver(r sct = (r, φ, z)), the scattered pressure P sct can be expressed as where scattering function S is S rigid (θ, φ; θ s , φ s ) = i n j n (ka) h n (ka) (−1) n (2n + 1)P n (cos β).
In (11), the solid angle β can be determined from cos β = cos θ cos θ s + sin θ sin θ s cos(φ s − φ). The reason why a rigid sphere is used here instead of a depth-extended cylinder in Sec.II, is: (1) The creeping waves' propagation speed on the sphere is a little slower than the cylinder's c l [34] for nth-order spherical Bessel function j n (x) = √ π/(2x)J n+1/2 (x), meaning their speed difference is c 2 l /(2kac+c l ); (2) In free field space, the cylinder's scattering function has a vertically standing wave term Correspondingly, in the waveguide, after plugged that into the Ingenito model, it gives The integral in (13) plays a model selection role. Next, a simulation is presented to show how multi-mode (multipath) propagation affects inverse scattering imaging. Here, the simulation's xOy horizontal configuration is similar to Fig.3 [35] and only increases R to 1 km. In the simulation, the water depth is 20 m while the source, the receiver and the rigid sphere are set to the same depth z = 10 m.

B. PHASE CONJUGATION-BASED CHANNEL DECONVOLUTION IMAGING METHOD
In Sec.III-A, it has been manifested that the multipath effect distorts the final ISCP images. However, (10) contains mode-coupling and cannot be easily resolved without vertically sampling the pressure field. According to [36], the scattering pressure of the rigid sphere can be expressed as Thus, after successful localization, phase conjugation can be used to deconvolve the CIR. However, the phase conjugation method introduces a large number of sidelobes as in Fig.12, which uses the same simulation data in Fig.11. That calls a better solution to conduct inverse scattering imaging in the underwater waveguide.  Fig.7, the image has at least lower the sidelobe with 2 dB. Also, the width of the ambiguity peaks in Fig.7 has been narrowed.

C. SPARSE RECONSTRUCTION-BASED CHANNEL DECONVOLUTION IMAGING METHOD
In this paper, a sparse reconstruction based channel responce deconvolution imaging scheme is proposed to yield a clear ISCP output. Here, SLIM is used to find the multipath  components at each receiver. The SLIM has a better time delay resolution than MF. For the same simulation data in Sec.III-A, the MF and SLIM output comparisons are displayed in Fig.13. Fig.13 shows that both MF and SLIM can find three main paths, but SLIM also gives the amplitude's polarity. When examining the receive data excited by the source at 0 • as Fig.13(c), MF's output has a strong sidelobe, while SLIM exhibits a better time and amplitude resolution. In the simulation, the time resolution is 1/6 ms because a 6 kHz bandwidth signal is used. Due to the long propagation range, the delay between two neighbor paths is larger than 1/6 ms. Fig.13 also tells that the creeping wave is hardly caught in the backward scattering direction, and indicates that creeping waves will contribute less in imaging. This naturally requires a larger far-field transmit and receive aperture to generate a complete image. All the paths found by SLIM are marked and categorized in Fig.14(a). The categorization standard includes the  polarity of amplitude and TOA relative to theoretical reflection's TOA and creeping wave's TOA, which refer to (7) and (8). Fig.14(a) agrees with Fig.13 that most categorized paths are caused by reflection waves, corresponding to the free-field case's phenomenon that creeping wave is the dominant component when the bi-static scattering angle is less VOLUME 8, 2020 than ±40 • . Thus, the output images, Fig.14(b)-(e), correspond to Path 1-5, respectively. The black dash line marks the target's contour line in the imaging plane, and a red-blue color map is used to distinguish the positive or negative amplitude path's imaging result. The first path is a positive amplitude wave, which is close to the theoretical reflection delay, resulting in the enveloping surface is on the contour line of the target (Fig.14(b)). The remaining paths have significant lags, therefore, making the enveloping surfaces in each image moving into the target's contour line. The enveloping surface in Fig.14(c) generated by the second path is still a circular arc with a smaller radius. When the time lag goes larger, the enveloping surface starts shrinking toward a point (Fig.14(d)) and then expands to a parabola (Fig.14(e)). For the paths corresponding to creeping waves, they suffer less identified parts and large time lags and only map into lines away from the target's contour (Fig.14(f)).
Group speed can help build the relation between path's TOA and normal modes propagation. The mth mode's group speed c gm approximates to [37] In our simulation, the sound speed in the water column is constant, so that (15) can be simplified as The expression above tells that the mth mode's group speed c gm can be derived from its corresponding horizontal wavenumber k rm . For the 3-9 kHz chirp signal, we take its center frequency 6 kHz as the reference and calculate the horizontal wavenumber. Thus, the paths distinguished by the SLIM algorithm has the following relation to the normal modes. Path 1 (blue '' * '' line in Fig.14(a)) corresponds to reflection wave's the first mode; Path 2 (blue ''•'' line) matches the reflection wave's the second mode; Path 3 (green '' * '' line) is the reflection wave's the third mode. Path 4 (green ''•'' line) is the reflection wave's the fifth mode; Path 5 (red ''•'' line) corresponds to the creeping wave, and its TOA matches Mode 5. Mode 4 is not identified by SLIM or MF. Mode with even number has a zero at target depth and Mode 4's derivative at the target depth is much smaller than the vertical wavenumber, making it harder to find its path. When the bi-static scattering angle is smaller than 20 • , the reflection wave becomes weaker and the creeping wave should be more possibly identified. However, in the simulation, SLIM missed the first three modes at most bi-static scattering angles less than 20 • . This phenomenon can be a further topic to study the physics behind that and will not be discussed hereafter. Actually, in the dispersive propagation media, the same mode has a different group speed value across the signal band. Both SLIM and MF ignore that effect. Inter-mode components are weaker than intra-mode components as stated in [38]. This section follows this point and treats the identified path propagated from the intra-mode coupling. The propagation time delay can be transformed from phase term exp(ik rm R) in (14). The propagation delay at frequency ω, nearly equals to group delay.
Thus, for each path component, a rectified time axis will let the path's TOA move back toward its own free-field's theoretical prediction as (7) and (8). The time axis transform of the mth mode is given as: Since SLIM neglects the dispersion, time axis scaling still takes wavenumbers at 6 kHz as reference. Fig.15(a) denotes each path's TOA value after TOA rectification. In Fig.15(a), all the identified paths move near the theoretically predicted values. All five transformed paths generate their new inverse scattering imaging as Fig.15(b) to Fig.15(f). These figures clearly show that after time rescaling, the enveloping surfaces are on the target's contour or tangent with the contour. The final output, which combines all the five paths, as shown in Fig.15(g), manifests that the sparse reconstruction based deconvolution method can generate an inverse scattering imaging result with lower sidelobe levels.

IV. CONCLUSION
This paper focuses on the inverse scattering imaging problems in a bistatic active sonar imaging. The time-domain back-propagation method is used for broadband signal probing cases. Moreover, combined with SLIM in two-dimensional free-field space, TD-BP not only reduces more sidelobe and raises at least a 2 dB peak-to-background ratio. The commonly used inverse scattering imaging methods will encounter defocusing and multiple ghosting introduced by a multipath effect in a shallow water environment. This problem is solved by using the SLIM algorithm, which is a sparse reconstruction method to identify the significant propagating paths. The identified paths are tagged into their possible categories, then rectified in the time domain. The proposed method achieves the effect of CIR deconvolution and yields a final low sidelobe inverse scattering image. Simulations in two-dimensional free-field space and three-dimensional range-independent waveguide are done to validate the methods developed in this paper. He is currently an Associate Professor with Zhejiang University. His research interests include statistical signal processing, speech processing, pattern recognition, and image processing.
T. C. YANG received the Ph.D. degree in high-energy physics from the University of Rochester, Rochester, NY, USA, in 1971. He was a Pao Yu-Kong Chair Professor with Zhejiang University, Hangzhou, China. From 2012 to 2014, he was a National Science Counsel Chair Professor with National Sun Yat-Sen University, Kaohsiung, Taiwan. He spent a period of 32 years with the Dispersive Wave Guide Effects Group, Naval Research Laboratory, Washington, DC, USA, as the Head of the Arctic Section, the Head of the Acoustic Signal Processing Branch, and a Consultant to the division on research proposals. He is currently a Professor with Zhejiang University. He has pioneered matched mode processing for a vertical line array and matched-beam processing for a horizontal line array. His current research interests include environmental impacts on underwater acoustic communications and networking, exploiting the channel physics to characterize and improve performance, environmental acoustic sensing and signal processing using distributed networked sensors, methods for improved channel tracking and data-based source localization, geoacoustic inversions, waveguide invariants, effects of internal waves on sound propagation in shallow water, arctic acoustics, and so on. He is a Fellow of the Acoustical Society of America.