Range-doppler algorithm of multireceiver synthetic aperture sonar using nonuniform signal

TL;DR

This paper introduces a new imaging algorithm for multireceiver synthetic aperture sonar that improves resolution in nonuniform conditions.

Contribution

A robust imaging algorithm for multireceiver SAS that handles nonuniform and redundant receiver cases effectively.

Findings

The proposed method models multireceiver SAS spectra using monostatic and bistatic components.

Sub-block filtering and interpolation techniques improve image resolution in nonuniform scenarios.

Simulations and experiments confirm the method's superior performance over traditional approaches.

Abstract

When the structures of multireceiver synthetic aperture sonar (SAS) are applied, imaging techniques are not always able to acquire images with outstanding resolution in nonuniform situations. In order to cope with this issue, we present a robust imaging algorithm in this paper. Based on our method, the receivers adopted by imaging algorithms are firstly determined. Then, we model the multireceiver SAS spectrum made up of the monostatic one and bistatic one. The bistatic portion is resolved via the application of sub-block filtering in the spectral domain. In the range-Doppler domain, the spectrums of determined receivers are coherently added. In the spectral domain, the compression in the range direction relying on multiplications is handled. With the help of interpolations, the coupling between both directions is removed in the range-Doppler domain, where the subsequent step is to do the compression in the azimuth direction. Our method can well focus the SAS datasets with the slight nonuniform case and redundant receiver case. However, traditional method just works well with the slight nonuniform case. In comparison with traditional method, our method can get high-resolution image all the time. The results based on simulations and experiments fully validate our method.

Figures

Geometric configuration of SAS system.

Phase error.(A) the range is 30m. (B) the range is 120m. (C) the range is 210m.

Flowchart of proposed imaging processor.

Imaging results with slight deviation.(A) FB approach, (B) our approach, (C) BP algorithm.

Azimuth profiles with slight deviation.(A) T1, (B) T2, (C) T3, (D) T4.

Imaging results with major deviation.(A) FB approach, (B) our approach, (C) BP algorithm.

Azimuth profiles with major deviation.(A) T1, (B) T2, (C) T3, (D) T4.

Imaging results of real SAS with major deviation.(A) FB approach, (B) Presented approach, (C) BP algorithm.

Equations

Full text

1. Introduction

As a mature underwater imaging technology [1,2], synthetic aperture sonar (SAS) [3–5] can carry out the imaging of the seafloor, and it is very important for the ocean observation [6–8] and underwater imagery. With the increasing demand of underwater imagery [9–11], SAS is viewed to be the good choice of small target detection [12–15] and recognition [16–20]. To achieve the large-scale swath, the emitter should emit signal with low pulse repetition frequency (PRF). To prevent spectrum ambiguity and retain outstanding resolution in the azimuth direction, it is necessary to apply large PRF which must be much broader than the Doppler bandwidth. As a result, the requirement for outstanding resolution and large-scale imaging area cannot be met by relying on traditional monostatic SAS systems. In order to cope with this contradiction, multireceiver synthetic aperture sonar (SAS) [21–24] indicating an emitter and a receiver array are mounted on the same sonar carrier is proposed. Using multiple receivers in the azimuth dimension, echoed signal from multiple receivers can be obtained at the same time. The demand of instantaneous sampling in range is alleviated by the enhancement of spatial sampling in azimuth. Compared to monostatic SAS, it simultaneously owns the ability of large-scale swath and high resolution [25,26].

For multireceiver SAS, the imagery is a core job. Based on the range-Doppler (RD) algorithm [27–29], the SAS imagery can be decomposed into three main steps, i.e., range compression, decoupling and azimuth compression. Since RD algorithm is simple and easily understood, it is popularly used by various systems. In [30], the RD algorithm based on hyperbolic slant range is presented. However, the higher order error which would have an impact on the imaging quality is ignored by this technique. The 2D spectrum and RD algorithm are designed in [31] via applying the stationary phase approach and the expansion of Taylor series with respect to the SAS overall range. In general, this method owns many approximations resulting in residual phase error. In [32], Zhongs et al develop the RD algorithm based on the lower-order Taylor series expansion of SAS overall range. With this method, the influence of approximation error on the migration curve is not corrected, and the imaging performance would be introduced by the higher sidelobes. In [33], Yangs et al present the RD algorithm based on the numerical computation technique, which is utilized to acquire stationary point and range migration function. In practice, the range migration on the basis of the calculated stationary point disregards the Doppler impacts giving rise to further phase error.

In general, the above methods need that the SAS system works in the ideal case [34,35], where the sonar moves precisely the half aperture of receiver array. Usually, we call it the uniform sampling case. Unluckily, the ocean is disturbed by the wind, current, etc. The moving speed of SAS system cannot be always met the constant requirement [36–38]. Therefore, the strict requirement cannot be satisfied, and the nonuniform sampling data would be caused in the azimuth dimension. When the SAS system suffers from nonuniform sampling, the above imaging methods cannot generate the high-resolution results. In order to cope with this issue, nonuniform data should be employed to recover the uniform data. In [39], the spectrum of nonuniform data is considered to be the weighted spectrum of uniform signal. Additionally, the weighting function [40] is barely affected by the azimuth and range frequencies. The main work of this method focuses on the spectrum calculation of nonuniform data and the weighting function. However, both works are time-consuming. Based on the filter bank (FB) method [41,42], each receiver data is considered to be the filtering results of nonuniform signal, and each receiver signal is uniformly sampled. From the filter inversion concept, the uniform signal can be efficiently reconstructed based on each receiver spectrum and inversion of filter bank. Some researchers have demonstrated that both methods are equivalent [43]. To some degree, both methods relax the SAS moving requirement between two adjacent pulses compared to traditional methods. However, this relaxation is just limited to a single receiver aperture in the azimuth dimension. When the deviation from the uniform sampling case is more than a receiver aperture in the azimuth dimension, the overlapped receivers and redundant signal between two adjacent pulses are caused. When the SAS system works in this case, neither of two uniform signal reconstruction techniques could provide images with excellent resolution.

Regarding our work of the paper, an imagery method which can simultaneously process the SAS data in the uniform and nonuniform cases is offered. With our method, we firstly choose the receivers adopted by our imaging algorithm. Based on this, we model the multireceiver SAS spectrum consisting of monostatic one and bistatic one, which is compensated for each receiver data in the spectrum domain. In the Doppler domain, the spectrums associated with chosen receivers are coherently superposed. We carry out the compression in the spectral domain of range dimension. The data is then converted into the range-time and azimuth-Doppler domain, the range migration is hereafter removed by correcting the decoupling between both dimensions with interpolations. After this operation, the compression and correction of shifting in azimuth is carried out. Compared to traditional method, our solution can greatly enhance the imaging capability.

The construction of our work is presented below. Section 2 discloses the geometric configuration and model of SAS system. Hereafter, we discuss our method in detail in the next section. The suggested approach is validated in section 4 with the help of tests and simulations. Lastly, the work in this paper is concluded.

2. Geometric configuration and model

Fig 1 illustrates the geometry of multireceiver SAS system consisting of one transmitter and M equally distributed receiver sensors. The transmitter is filled with black while no colors are filled by the receivers. The system [44] moves in the azimuth dimension with an ideal velocity V. We suppose that the position of transmitter in the azimuth dimension is 0 when the system plans to work. After some working time $[eqn]$ , the azimuth position of transmitter is $[eqn]$ . The range history from transmitter location to an assumed object (r, 0) can be represented as

Geometric configuration of SAS system.

$$[eqn]$$

In Eq. (1), the slow time and range are denoted by $[eqn]$ and r, individually. The receiver location is $[eqn]$ in the azimuth dimension, and $[eqn]$ denotes the limited distance of the subsystem made up of the sender and the mth receiver sensors. The distance from this receiver location to the assumed object (r, 0) is expressed by

$$[eqn]$$

We assume that the chirp waveform [45,46] is emitted, and the echo collected by the mth receiver is provided by

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the acoustic velocity and carrier frequency, individually. The fast time is denoted by the symbol $[eqn]$ .

3. Multireceiver SAS RD algorithm

3.1. 2D Spectrum

The main job of this section is to deduce the 2D spectrum of echo shown in Eq. (3).

Before discussing our imaging algorithm, we firstly present 2D spectrum [47–49], which is an indispensable work of fast focusing algorithms. Depending on the Fourier transform (FT) with respect to the fast time, Eq. (3) is converted into the spectral domain in the range direction and the result is written as

$$[eqn]$$

where the instantaneous frequency related to the fast time is shown by $[eqn]$ . $[eqn]$ is the spectrum of $[eqn]$ in Eq. (3).

In azimuth direction, exploiting the FT generates the 2D spectrum, and it is indicated by

$$[eqn]$$

where the Doppler frequency is shown by $[eqn]$ in the azimuth dimension.

For ease of clarification, we define the emitter phase and receiver phase $[eqn]$ and $[eqn]$ . Based on the stationary phase principle [50], two stationary phase points corresponding to emitter phase and receiver phase are developed, and they are indicated by $[eqn]$ and $[eqn]$ , respectively. With the quadratic expansion, the emitter and receiver phases are given by

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the first-order derivative with respect to the slow time. The first-order derivative is zero according to the stationary phase principle, i.e., $[eqn]$ and $[eqn]$ . $[eqn]$ and $[eqn]$ are the second-order derivative.

Based on Eq. (6) and Eq. (7), the emitter and receiver phases in Eq. (5) are replaced by the approximated terms, and the resultant integral shown in Eq. (5) is reformulated to

$$[eqn]$$

Exploiting the stationary phase principle again, Eq. (8) would be solved and the analytical expression of 2D spectrum is provided by

$$[eqn]$$

Examining Eq. (9), the initial component represents the chirp spectrum. The subsequent two items are the bistatic phase arising from the displaced separation between emitter and receiver. The monostatic SAS spectrum is the precondition of RD algorithm, and its expression is very close to the fourth factor in(9). The stop-and-hop assumption usually applied in the radar field [51–53] results in the deviation [54], and it is shown by the final item.

The SAS system owns 50 receivers, and the total aperture of receiver array is 1.5 m. The aperture of emitter is 0.06 m. The carrier frequency of sonar system uses 150000 Hz. 20000 Hz is used by the chirp waveform bandwidth. The interval of sent chirp is 300 ms. The sonar is operated with 2.5 m/s. The system works with uniform case. The phase error of presented spectrum shown in Eq. (9) increases with the aperture. Based on the accurate spectrum, the phase error of largest aperture is displayed in Fig 2. When the SAS focusing is performed, the error not exceeding $[eqn]$ would not noticeably affect the SAS quality. The maximum error of our spectrum is not exceed to $[eqn]$ . Therefore, the preseted spectrum shown in Eq. (9) has the ability to generate high quality results.

Phase error.(A) the range is 30m. (B) the range is 120m. (C) the range is 210m.

3.2. Imaging processor

Based on the 2D spectrum shown in Eq. (9), the work of this part is devoted to the development of RD processor [28] for multireceiver configuration. According to the sonar moving velocity and the pulse repetition interval (PRI), the moving distance between two adjacent pulses is given by $[eqn]$ . In practice, the azimuth sampling should be adhering to the Nyquist sampling like the time-domain sampling. The azimuth sampling is strictly carried out based on the receiver, and each receiver aperture is the sampling interval. That is to say, the moving distance per pulse should be strictly half length of the array. Since the spatial sampling interval is too small, the deviation caused by the nonuniform sampling can be neglected. This is the key of our approach. Using the real aperture of each receiver, the number of receivers exploited by the imaging processor is conformed, and it is shown as

$$[eqn]$$

where d indicates the aperture of receiver. $[eqn]$ shows that the nearest integer is selected by the imaging processor.

Based on Eq. (9), the matched filter is used to perform the compression in range, and the filter item is written by

$$[eqn]$$

where $[eqn]$ shows the conjugate operation.

Following the aforementioned filtering, every receiver data is divided into Q blocks. Every block data is converted into spectral domain, and we aim to correct the second item in Eq. (9). For the qth block, the correction term is given by

$$[eqn]$$

where $[eqn]$ denotes the center of each block.

Actually, the large error after correction lies in the edge of this block. When the SAS focusing is performed, the error not exceeding $[eqn]$ would not noticeably affect the SAS quality. We define the block length $[eqn]$ . Then, this relationship is shown as

$$[eqn]$$

Eq. (13) is further reformulated as

$$[eqn]$$

Based on Eq. (14), the block length is got. Consequently, the total block number Q is got.

All block datasets in the spectral domain are then transformed into fast time domain, where all blocks are integrated together to generate an entire receiver signal.

Subsequently, we plan to eliminate the impacts of the bistatic phase expressed by the third term in Eq. (9), and the compensation function is expressed as

$$[eqn]$$

After this compensation, the coherent summation of all sensor spectrums is operated, and this operation is shown as

$$[eqn]$$

When we complete the coherent summation, all receiver datasets are coerced into the uniform equivalent signal. To compensate the fourth term in Eq. (9), it is expanded up to the quadratic term, and the expression is shown as

$$[eqn]$$

where a temporary symbol $[eqn]$ is defined. The primary item in Eq. (17) indicates the modulation in azimuth. The subsequent item represents the range and azimuth coupling. The final factor in Eq. (17) shows the chirp characteristic resulted from the quadratic coupling of both dimensions.

In the spectral domain, the approximation error that arises from Eq. (17) is compensated. Besides, the quadratic coupling in Eq. (17) is corrected. To perform this step, the filter is supplied by

$$[eqn]$$

where $[eqn]$ denotes the range of particular target in the imagery region.

Following the preceding action, the inverse Fourier transform (IFT) is applied to datasets in the range direction, and the datasets are converted into the range time domain from range frequency domain. Then, the linear coupling between both dimensions in Eq. (17) is eliminated via making use of interpolation technique, and the delay offset arising from the range cell migration is represented by

$$[eqn]$$

The term shown in Eq. (19) is not the same for each range bin. Therefore, it cannot be corrected by using the multiplex multiplication at the same time. Considering this factor, the linear coupling between both directions in Eq. (17) is cancelled via the application of interpolation technique. The physical meaning of this step eliminates the coupling between two directions.

After the former handling stage, we aim to do the matched filtering and deviation correction in azimuth. The physical meaning of this step performs the compression in the azimuth direction. We provide the handling filter shown as

$$[eqn]$$

In the azimuth dimension, we do the IFT, and the image with outstanding quality is created. Based on the processing mentioned above, Fig 3 shows the processing structure of described solution.

Flowchart of proposed imaging processor.

3.3. Computation complexity

The basic operations used by our method mainly include the FT/IFT, multiplex multiplication and interpolation. In [55], the computation complexity for these operations has been discussed. In this section, the size of original data matrix is supposed to be $[eqn]$ . On the basis of this, we directly present the computation complexity of our method, and it is shown as

$$[eqn]$$

where N is the receiver numbers exploited by our imaging processor. The interpolation length is represented by $[eqn]$ .

It is apparently visible from Fig 2 that our solution mainly exploits the FT/IFT and multiplex multiplication. Therefore, the presented method has the characteristics of high efficiency.

4. Simulations and experiments

4.1. Simulation discussion

This part mainly concentrates on the validations of presented method based on simulations. In the imaging area, four targets are at 50m, 100m, 150m and 200m are focused, and the targets are labeled T1, T2, T3 and T4, individually. The corresponding azimuth positions are 5m, 8m, 12m and 15m, individually. The SAS system owns 50 receivers, and the total aperture of receiver array is 1.5 m. The aperture of emitter is 0.06 m. The carrier frequency of sonar system uses 150000 Hz. 20000 Hz is used by the chirp waveform bandwidth. The interval of sent chirp is 300 ms. Here, the sonar is operated with 2.475 m/s. Based on this speed, the system should need 49.5 receivers. Traditional nonuniform reconstruction method can well work in this case. The performance of presented method can be compared to the high performance of traditional method. Based on these configurations, we discover that the SAS system is nonuniformly sampled by the multiple receiver sensors installed along the moving direction. With this case, the forward deviation of consecutive two pulses is shorter than the real aperture of a single receiver, and there are no redundant receivers. Due to this reason, we call it the case of slight deviation.

Using the FB method [41], presented approach in this paper and back projection (BP) algorithm [56], the handling images are acquired. Here, the images using BP algorithm are the benchmark. Fig 4A demonstrates the results using FB technique, and Fig 4B displays images via the presented approach. Fig 4C shows the results of BP method. Comparing imaging outcomes of three methods, it is apparent that three methods can precisely focus all targets. It shows that our approach in this paper can get the high-resolution image when the moving deviation between two adjacent pulses is not less than the real aperture of a single receiver.

Imaging results with slight deviation.(A) FB approach, (B) our approach, (C) BP algorithm.

Fig 5 illustrates the azimuth profiles of four targets. As seen in the profiles of all targets, the provided approach mostly achieves identical profiles relying on the FB method [41] and BP algorithm [56]. This still proves the efficacy of provided approach.

Azimuth profiles with slight deviation.(A) T1, (B) T2, (C) T3, (D) T4.

The quality metrics can assist us to assess the performance of various approaches quantitatively. We would compute the quality metrics to further assess the quality of both methods. In Table 1, the quality metrics such as the peak sidelobe ratio (PSLR), integrated sidelobe ratio (ISLR) and azimuth resolution (AR) are adopted to compare the focusing quality of both approaches. Based on the quality metrics in Table 1, the quality metrics of presented approach is slightly inferior to those of FB method [41] and BP algorithm [56]. These distinctions are generally so insignificant that these differences can be disregarded. The discussed technique is capable of obtaining a fine image. That is to say, the presented method can get the high-resolution image. The findings gathered from Fig 4 and Fig 5 are generally in agreement with that based on Table 1.

Table 1: Imaging quality with slight deviation.

The second simulation concentrates on the major deviation between two adjacent pulses. In this simulation, the towed speed of sonar system is changed to 2.2 m/s. Based on this speed, the system needs 44 receivers, and 6 redundant receivers are exited. With this case, traditional nonuniform reconstruction method cannot well reconstruct the signal. However, our method can still work with this case. The other parameters are identical to parameters of previous case. Based on the parameters of this case, it should be noted that there are 6 redundant receivers. Using the FB method [41], discussed method and BP algorithm [56], four targets are focused. Fig 6 exhibits the handling outcomes. The outcomes related to FB approach [41] are illustrated in Fig 6A. We can clearly see that the targets are failed to be focused. Actually, the redundant receivers make the failure of reconstruction filters, which is based on the matrix inversion. Due to this reason, the targets cannot be successfully focused. Based on the discussed approach, all targets can be still focused with the high-resolution. The performance between our method and BP algorithm is nearly the same. This shows that our approach is better than FB method [41] when the SAS system works in the case of major deviation.

Imaging results with major deviation.(A) FB approach, (B) our approach, (C) BP algorithm.

Here, we still pay attention to the profiles in the azimuth dimension. The profiles related to focused targets in the azimuth dimension are presented in Fig 7. Table 2 shows the quality metrics including PSLR, ISLR and AR. Based on Fig 7 and Table 2, the FB method cannot focus the targets any more. With the presented approach, the targets are well focused. The quality parameters based on our method are slightly inferior to those of BP algorithm. In conclusion, the suggested technique is capable of significantly enhancing the handing quality even if the system operates with the major deviation case. As a result, we figure out that the suggested solution is capable of achieving outstanding outcomes in any cases.

Table 2: Imaging quality with major deviation.

Azimuth profiles with major deviation.(A) T1, (B) T2, (C) T3, (D) T4.

4.2. Experiment and discussions

Currently, the focusing in the case of major deviation is discussed. With our system, 48 receivers are exploited and the real aperture of each receiver in the azimuth dimension is 0.04 m. The carrier frequency and bandwidth of the chirp waveform are 150000 Hz and 20000 Hz, respectively. The towed speed of sonar system is 2.5 m/s, and PRI is 0.32 s. With these parameters, it is easily found that the system suffers from 8 redundant receivers. Exploiting the FB approach [41], proposed approach and BP [56], the real datasets are handled. Fig 8A shows the outcome of FB approach. The outcome of suggested approach can be seen in Fig 8B. The outcome based on BP algorithm is presented in Fig 8C. Since discussed SAS system suffers from many redundant receivers, the reconstruction filters of FB method [41] based on the computation of matrix inversion would be not a success job. Consequently, the targets cannot be focused by exploiting the FB approach [41], and Fig 8A confirms this conclusion. The presented method well solves this issue. Therefore, the target can be successfully focused, and Fig 8B confirms this conclusion. Furthermore, the performance based on our method in Fig 8B is close to that based on BP algorithm in Fig 8C. This further verifies the efficacy of our approach. The handling time of this experiment is recorded. The FB approach and our approach need 848 s and 145 s, respectively. The BP algorithm needs 7531 s. Based on the experiment in this part, we can get the conclusion that our method is much better than FB method in the case of major deviation.

Imaging results of real SAS with major deviation.(A) FB approach, (B) Presented approach, (C) BP algorithm.

In Table 3, we still present the entropy and contrast [57] corresponding to images in Fig 8. Based on Table 3, the entropy and contrast of FB method dramatically deviate from corresponding parameters of our method and BP algorithm. That is to say, traditional method fails to recover the SAS image. Comparing the entropy and contrast of our method and those of BP algorithm, the parameter differences between our method and BP algorithm can be disregarded. It means that, our method can recover the SAS image with high quality.

Table 3: Imaging quality with the major deviation case.

5. Conclusion

This article offers the RD algorithm for SAS made up of many receiver sensors. With our approach, the receivers applied to our imaging processor are determined. Then, the 2D spectrum including the bistatic term and monostatic equivalent term is deduced based on the quadratic expansion of emitter and receiver phases. With the help of the block correction approach, tackling the bistatic term would yield the monostatic equivalent term, which is further compensated based on the quadratic decoupling, linear decoupling and azimuth compression. Using our approach, the multireceiver SAS datasets can well be focused when the SAS works in the case of slight deviation. In the case of major deviation, our method can still get the high-resolution image while FB method suffers from serious distortion. The focused results based on the simulated datasets and real SAS datasets verify the merits of our approach.

Based on our method, we can process the SAS datasets with high performance all the time in any cases. However, traditional reconstruction method just works in slight nonuniform case. Our work does not consider the motion compensation, which is our future work.

Supporting information

S1 FileThis is the raw data file of experiment.(XLSX)