**4. Imaging Algorithm**

The key issue of the SAS imagery is the decoupling. Inspecting (10), the coupling phase is a function of the range, instantaneous frequency and Doppler frequency. Due to the range variance, we cannot perform the decoupling operation in the range Doppler domain. Based on the characteristic of the coupling phase, the sub-block processing method is used to perform the decoupling operation. In this section, important steps of the presented method are introduced in detail.

#### *4.1. 2D FT*

In this step, each receiver's data is transformed into 2D frequency domain based on the FT. Each receiver data is undersampling. The energy of the frequency *f <sup>t</sup>* ∈ [−*PRF*/2, *PRF*/2] is a combination of all the energy related to the frequency points (··· , *f <sup>t</sup>* − *PRF*, *f <sup>t</sup>* , *f <sup>t</sup>* + *PRF*, ···) ∈ [−*M* · *PRF*/2, *M* · *PRF*/2]. Here, *f <sup>t</sup>* ∈ [−*PRF*/2, *PRF*/2] denotes the Doppler frequency related to the sampling rate of each receiver's data. In other words, the spectrum in the Doppler domain is aliased *M* times. Here the pulse repetition frequency is denoted by *PRF*, which is also the sampling frequency of

each receiver's data in the azimuth dimension. The alias can be suppressed by the coherent processing of multireceiver data.

#### *4.2. Decoupling*

The coupling phase shown as (10) is characterized by range variance. Although the coupling phase is accurate, it is hard to compensate the coupling phase completely. Here, the sub-block processing method is exploited to perform the decoupling operation. Since it is hoped that the echo signal of a target is in the same range cell (usually at *r*), the range deviation from its desired position denotes the RCM. Based on (10), the RCM in the 2D frequency domain is expressed as:

$$
\varphi\_{\rm dc\\_i}(f\_{\sf T}, f\_{\sf t}; r) = -\varphi\_i(f\_{\sf T}, f\_{\sf t}; r) - 4\pi f\_{\sf T} \frac{r}{c} \tag{11}
$$

From (11), we see that the RCM is a function of the range, instantaneous frequency and Doppler frequency. Based on this characteristic, the sub-block processing method rather than interpolation is exploited to carry out the range cell migration correction (RCMC). The decoupling is decomposed into two steps. One is the bulk decoupling, and the other is the differential decoupling. Based on (11), the filtering function of the bulk decoupling is given by:

$$H\_{\rm bulk\\_j} = \text{conj}\{P(f\_\tau)\} \cdot \exp\{j\varphi\_{\rm de\\_j}(f\_{\tau\prime}f\_t; r\_\mathbf{c})\}\tag{12}$$

where *r*<sup>c</sup> represents the center range of the mapping swath; conj(·) represents the complex conjugate.

The bulk decoupling simultaneously performs the range matched filtering. Considering targets at reference range, the coupling is completely removed after this step. However, other targets suffer from the residual error, which is *ϕ*de\_*i*(*fτ*, *ft*;*r*) − *ϕ*de\_*i*(*fτ*, *ft*;*r*c). Based on the sub-block processing method, the differential decoupling is used to solve this problem. Before performing the differential decoupling, the whole swath is virtually segmented into *N* sub-blocks in the range direction. Based on (11) and (12), the filtering function of the differential decoupling is written as

$$H\_{\rm i\\_n} = \exp\{j\varphi\_{\rm de}(f\_{\tau\,\prime}f\_{t\,\prime}; r\_{\rm ref\\_n}) - j\varphi\_{\rm de\\_i}(f\_{\tau\,\prime}f\_{t\,\prime}; r\_{\rm c})\}\tag{13}$$

where the center range of the *n*-th sub-block is used as the reference range, which is denoted by *r*ref\_*n*. The variable *n* ∈ [1, *N*] denotes the sub-block index.

Based on (13), we remove the coupling between the range and azimuth dimensions for the data of the *n*-th sub-block. Then, the data is transformed into the time domain via the inverse Fourier transformation (IFT) in the range direction. For each sub-block, the differential decoupling is carried out using (13). The processed sub-blocks are extracted and stored in the range Doppler domain. The coupling between the range and azimuth dimensions is cancelled when *N* sub-blocks are reckoned into a new signal matrix.

### *4.3. Azimuth Compression*

The sonar transmits and receives signals and then stores the received signal. This process is conducted at a set of locations along the moving trajectory. The signal collected at these locations forms different units of synthetic aperture called azimuth sample serials. By processing the data collected within the synthetic aperture time, an equivalent large sonar array can be obtained. By coherently processing the echo data, we obtain the high resolution in the azimuth dimension.

In practice, the azimuth echo of the SAS system can also be regarded as a wideband signal. The matched filter can still be used to compress the azimuth signal. According to the theory of the matched filter [9], it is necessary to set parameters of the matched filter to be the same as Doppler parameters of the azimuth signal. As a result, parameters should be adjusted at different ranges in order to get high resolution in azimuth across the whole swath. This is the main difference between the range compression and azimuth compression.

Inspecting (8), the azimuth modulation only depends on the range *r* and Doppler frequency *ft*. It is not a function of the instantaneous frequency *f<sup>τ</sup>* any more. This term is responsible for the azimuth matched filtering after performing the decoupling between the range and azimuth dimensions. Considering the new signal matrix in the range-Doppler domain, the azimuth compression is directly carried out based on the azimuth modulation. The filtering function is given by:

$$H\_{\text{ac\\_i}} = \exp\{-j\varphi\_{\text{ac\\_i}}(f\_t; r)\}\tag{14}$$

After this step, the signal is compressed in azimuth. However, the azimuth spectrum of each receiver's data is aliased *M* times due to the (1/*M*)-th undersampling.

#### *4.4. Coherent Superposition*

Monostatic SAS system has a transducer, which is used as the transmitter and receiver at different times. The wide swath in range requires a low rate of *PRF*. However, the high resolution in azimuth requires high rate of *PRF*. In other words, the different requirements of the pulse repetition frequency in range and azimuth dimensions lead to a trade-off between swath width and azimuth resolution.

To solve this issue, the multireceiver SAS including a transmitter and receiver array is used. When the sampling of each individual receiver occurs at a rate of *PRF*, the effective sampling of the equivalent monostatic system is *M* · *PRF*. After transmitting a pulse, *M* samples can be recorded. Based on the PCA method [12], the unambiguous spectrum satisfying the Nyquist rate is recovered before the SAS imagery. With the PCA method, the PTRS of multireceiver SAS is decomposed into two parts. One depends on *di*. The other is independent of *di*, and it is similar to the PTRS of monostatic SAS. When it comes to recover the unambiguous spectrum, the phase related to *di* must be simultaneously compensated. The subsequent processing can be done with the help of imaging algorithms based on monostatic SAS system. With the presented method, the PTRS cannot be decomposed into aforementioned two parts. According to the theory of the linear time invariant system [9], we can exchange the processing order between the SAS imagery and recovery of the unambiguous spectrum. In other words, each receiver's data can be focused before recovering the unambiguous spectrum. Considering each receiver's data, we obtain focusing results in the range Doppler domain by using the decoupling and azimuth compression. Since each receiver's data is sampled by the pulse repetition frequency *PRF*, the spectrum of a single receiver data must be aliased. Fortunately, the coherent processing of multireceiver data can suppress this alias. The spectrums of all receiver data are coherently superposed in the range Doppler domain. The resultant data would satisfy the Nyquist rate in azimuth, and the equivalent sampling rate is *M* · *PRF*. The high resolution image is obtained after an azimuth IFT.

According to the presented steps, the block diagram of the proposed algorithm is shown in Figure 2.

**Figure 2.** Block diagram of the presented processor. See the text (sections: PTRS, azimuth modulation and coupling term, and imaging algorithm) for all terms and full names of abbreviations used in the figure.

#### **5. Comparison with Traditional Methods**

In this section, we begin with the comparison of the PTRS, which is the basis of fast imaging algorithms. To directly use the method of stationary phase [9], the series expansion of the two-way slant range is often used by traditional methods [12,16,19]. From Figure 1, the two-way slant range can be formulated as:

$$R\_i(t;r) = c\tau\_i = \sqrt{r^2 + \left(vt\right)^2} + \sqrt{r^2 + \left(vt + v\tau\_i + d\_i\right)^2} \tag{15}$$

where *r*<sup>2</sup> + (*vt*) <sup>2</sup> denotes the instantaneous range between the target and transmitter, and *r*<sup>2</sup> + (*vt* + *vτ<sup>i</sup>* + *di*) <sup>2</sup> is the instantaneous range between the *i*-th receiver and target. Here, *<sup>v</sup>τ<sup>i</sup>* represents the forward distance during the signal reception. Considering the complex expression of the accurate time delay, *v* · *τ<sup>i</sup>* is often approximated by 2*vr*/*c* for simplicity.

(15) is the sum of two square roots, which make it difficult to acquire the PTRS based on the method of stationary phase [9]. To solve this problem, (15) is expanded into a power series, which is given by:

$$\mathcal{R}\_i(t;r) \approx \sum\_{q=0}^{Q} k\_{i\\_q} t^q \tag{16}$$

where *k*-coefficients in (16) can be calculated by the rule of series expansion [9]. To obtain the analytical PSP, the parameter *Q* often satisfies 2 ≤ *Q* ≤ 4.

Inspecting (16), the series expansion would lead to the approximation error. Besides, the error of stop-and-hop approximation is not completely compensated. Since the presented method is based on the accurate time delay, both issues are successfully avoided.

We now discuss the link of the SAS imagery between the presented method and traditional methods [12,16,19]. Expanding the PTRS as a Taylor series with respect to the instantaneous frequency, various imaging algorithms are developed. The second order series approximation is exploited by the range-Doppler (R-D) algorithm [24] and chirp scaling (CS) algorithm [25]. The nonlinear CS algorithm [26,27] and quartic phase algorithm [28] are based on the third and fourth order approximation, respectively. The range migration algorithm (RMA) [29,30] requires that the PTRS is a linear function of the range. Generally, the series approximation of the PTRS is unavoidable if traditional fast algorithms were still exploited. Based on traditional methods, we conclude that the azimuth modulation is independent of the instantaneous frequency. With this characteristic, we first derive the azimuth modulation by setting *f<sup>τ</sup>* = 0 in (7). Then, the difference between the PTRS and azimuth modulation denotes the coupling term. Consequently, the Taylor expansion of the PTRS is avoided. Based on (7), (8) and (10), the block diagram related to the deduction of the PTRS, coupling term and azimuth modulation is shown in Figure 3a. The dashed rectangle in Figure 3a denotes important terms deduced by the presented method. Figure 3b shows the deduction of traditional methods. It can be found that traditional methods are very tedious. In this paper, we provide a novel aspect for the accurate deduction of the PTRS, azimuth modulation and coupling term.

**Figure 3.** Deduction of three important functions. (**a**) Presented method; (**b**) traditional methods. See the text (section: Introduction) for all terms and full names of abbreviations used in the figure.

In general, our method owns two major advantages. Traditional fast imaging methods cannot directly use the accurate time delay, which is often exploited by the BP algorithm. Approximations are usually used to deduce the PTRS, azimuth modulation and coupling term. Unfortunately, the approximations degrade the imaging performance. In this paper, the PTRS, azimuth modulation and coupling term are deduced based on the accurate time delay. Consequently, these functions avoid approximations. It is the first advantage of the presented method.

The second advantage is that our imaging scheme can be simply extended to any other PTRS. The presented imaging scheme does not require the series expansion of the PTRS with respect to the instantaneous frequency. With the presented method, the key step is to deduce the azimuth modulation. Considering the analytic PTRS, the azimuth modulation can be directly obtained by setting *f<sup>τ</sup>* = 0 in the PTRS. When the PTRS is complicated, the azimuth modulation is deduced by using two steps. The first step should calculate the azimuth PSP by setting *f<sup>τ</sup>* = 0 in PSP. The subsequent step sets *f<sup>τ</sup>* = 0 in the PTRS and substitutes the azimuth PSP into the PTRS. Based on the PTRS and azimuth modulation, the coupling term can be calculated. After carrying out the decoupling operation based on the sub-block processing method, the imaging process is decomposed into two separate filtering processes in the range and azimuth dimensions.

#### **6. Simulations and Real Data Processing**

#### *6.1. Simulation Results*

In this section, simulations are exploited to validate the presented method. The SAS parameters are listed in Table 1.


**Table 1.** The SAS System Parameters.

#### 6.1.1. Processing Results of Presented Method

To understand the presented method, the processing results of the main steps are discussed in detail. For clarity, we suppose that there is a point target located at coordinates (141 m, 17 m). Considering the first receiver's data shown in Figure 4a, we carry out the bulk decoupling in the 2D frequency domain. The resultant signal is shown in Figure 4b. Then, we perform the differential decoupling, and the result is shown in Figure 4c. After this step, the azimuth compression is conducted based on the azimuth modulation. Figure 4d depicts the result after the azimuth compression. Each receiver's data is undersampled using the pulse repetition frequency, which cannot satisfy the Nyquist rate in azimuth. Due to this reason, all results shown in Figure 4 are aliased in the azimuth dimension. Fortunately, the alias can be suppressed by processing the multireceiver data coherently.

Inspecting Figure 4c,d, we find that both results are visually indistinguishable. In fact, the result shown in Figure 4c is considered to be the input of the subsequent filter, i.e., azimuth compression. Since the azimuth compression performs the phase compensation in the frequency domain, the signal magnitudes are visually indistinguishable in the frequency domain. However, the major difference can be found in the space domain. Applying the azimuth IFT to Figure 4c,d, Figure 5 shows the results in the 2D space domain. From Figure 5b, the signal shown in Figure 5a is compressed in the azimuth dimension, and the circled part represents the recovered target. Each receiver's data sampled by the

pulse repetition frequency is undersampling. Due to this reason, the ghost targets are introduced. The coherent processing of multireceiver data can solve this problem.

**Figure 4.** Processing results of a single receiver data. (**a**) Single receiver data; (**b**) bulk decoupling; (**c**) differential decoupling; (**d**) azimuth compression.

**Figure 5.** Results in the two-dimensional (2D) space domain. (**a**) Differential decoupling; (**b**) azimuth compression.

Based on the steps in Section 4, *M* results corresponding to *M* receiver data are obtained. Each result is similar to Figure 4d. We coherently superpose *M* results, and the signal is shown in Figure 6a. The coherent superposition is equivalent to the improvement of the sampling rate in azimuth. Therefore, the data shown in Figure 6a satisfies the Nyquist rate, which is increased to *M* · *PRF*. Performing an azimuth IFT, we obtain the high resolution image, which is shown in Figure 6b. To visually examine the focusing performance, we depict the azimuth slice corresponding to Figure 6b. The azimuth slice is shown in Figure 6c. For comparison, the data shown in Figure 4a is directly processed by the presented method. The azimuth slice corresponding to Figure 5b is also depicted in Figure 6c. From Figure 6c, the azimuth slice related to a single receiver data is aliased,

as each receiver data is undersampled by the pulse repetition frequency. By coherently processing multireceiver SAS data, we obtain the high resolution image. Therefore, we conclude that the presented method successfully focuses the point target.

**Figure 6.** Processing results of all receiver data. (**a**) Coherent superposition; (**b**) focused target; (**c**) azimuth slice.

#### 6.1.2. Influence of Sub-Block Width on Imagery

In general, the imaging algorithm can be decomposed into two steps. The first step is to derive the PTRS, coupling term and azimuth modulation. With the presented method, we accurately deduce three terms. The second step is to design the imaging algorithm based on the PTRS, coupling term and azimuth modulation. Since the coupling term in the 2D frequency domain is range variant, we perform the decoupling based on the sub-block processing method. However, it is hard to compensate the coupling term completely. In other words, the sub-block processing method results in the residual error. Here, we discuss the influence of the residual error on the SAS imagery. The imaging scenario consisting of 18-point targets is shown in Figure 7. The targets are marked by T1, T2, ... , and T18, respectively. T1, T2, ... , and T6 are located at close range. T7, T8, ... , and T12 are located at medium range. The remaining targets are supposed to locate at far range. The coupling term is characterized by space variance, which may lead to the space variance of the optimal sub-block width. When the performance of the presented method is not inferior to that of traditional method, the sub-block width used by the presented method is defined as the optimal width of the sub-block. In Figure 7, we depict three sub-blocks, which are denoted by the red, blue and pink rectangles. Each sub-block consists of six targets.

**Figure 7.** Simulated scenario with 18-point targets.

We first focus on the targets circled by the red rectangle in Figure 7. In this sub-block, the reference range used by the differential decoupling is 53 m. The difference between the target range and reference range denotes the half width of the sub-block. The BP algorithm [18] based on (1) is viewed as the precise method. Here, the results of the BP algorithm [22] are used as the criteria. For multireceiver SAS systems, the PCA method [12] is widely used. At this point, we mainly conduct the comparison between the presented method and PCA based R-D algorithm. Figure 8 shows the azimuth slices of T1, . . . , and T6.

**Figure 8.** *Cont*.

**Figure 8.** Azimuth slices of close targets. (**a**) Sub-block width with 20 m; (**b**) sub-block width with 14 m; (**c**) sub-block width with 12 m; (**d**) sub-block width with 10 m; (**e**) sub-block width with 8 m; and (**f**) sub-block width with 2 m.

From Figure 8, the performance of the presented method is improved by decreasing the sub-block width. The large width of the sub-block generates great residual error, which noticeably degrades the imaging performance. When the sub-block widths are 14, 12 and 10 m, the performance of the presented method is inferior to that of the PCA method. The performance of the presented method is nearly close to that of the PCA method when the sub-block width is decreased to 8 m. In other words, this sub-block width can satisfy the imagery with high performance at close range. The improvement is still obtained when we choose much narrower sub-block width. However, the improvement is not noticeable. Figure 8f enhances this conclusion.

Next, we use the peak sidelobe level ratio (PSLR) and integral sidelobe level ratio (ISLR) to evaluate the imaging performance. The quality parameters are shown in Table 2.


**Table 2.** Quality parameters for targets at close range. PCA: phase center approximation; BP: back projection; PSLR: peak sidelobe level ratio; ISLR: integral sidelobe level ratio.

From Table 2, we see that the presented method obtains a low resolution image with a large sub-block width. The PSLR and ISLR related to T1, T2 and T3 enhances this conclusion. For the target T4, the focusing performance of the presented method is mostly close to that of the PCA method. Considering the target T5, the focusing performance of the presented method is superior to that of the PCA method. In this case, the residual error introduced by the sub-block processing method can be negligible. Inspecting quality parameters of T6, the focusing performance is improved by decreasing the sub-block width. However, the improvement is slight, as the residual error does not dramatically influence the imaging performance. Therefore, we conclude that the sub-block width with 8 m can satisfy the high performance imagery at close range.

We now concentrate on the targets at medium range. In this case, the reference range used for the differential decoupling is 143 m. After the data processing, the azimuth slices are shown in Figure 9.

**Figure 9.** Azimuth slices of medium range targets. (**a**) Sub-block width with 20 m; (**b**) sub-block width with 14 m; (**c**) sub-block width with 12 m; (**d**) sub-block width with 10 m; (**e**) sub-block width with 8 m; and (**f**) sub-block width with 2 m.

Inspecting Figure 9, we nearly obtain the same conclusions drawn from Figure 8. When the targets are at medium range, the optimal width of the sub-block is still 8 m. Figure 10b also strengths this conclusion. Since the stop-and-hop error is not completely compensated, the performance of the PCA method is slightly degraded. Based on Figures 8 and 9, the optimal sub-block width of both cases is mostly identical. In other words, the optimal width of the sub-block is nearly range invariant. Table 3 lists quality parameters for targets at medium range.

**Figure 10.** Azimuth slices of far targets. (**a**) Sub-block width with 2 m; (**b**) sub-block width with 8 m; (**c**) sub-block width with 10 m; (**d**) sub-block width with 12 m; (**e**) sub-block width with 14 m; and (**f**) sub-block width with 20 m.


**Table 3.** Quality parameters for targets at medium range.

Based on Table 3, the performance of the presented method is improved with the decreasing of the sub-block width. When the sub-block width is decreased to 8 m, we obtain the image which outperforms that of the PCA method. Due to the residual error of the stop-and-hop approximation, the quality parameters of the PCA method are slightly lowered.

The last experiment focuses on the imaging performance at far range. The reference range used for the differential decoupling in this sub-block is 194 m. The azimuth slices are shown in Figure 10. The residual error of the stop-and-hop assumption increases with the range. Consequently, the PCA method suffers from large residual error when the targets are located at far range. With the presented method, the imaging performance is nearly close to that of the PCA method when the sub-block width is 14 m. Figure 10e enhances this conclusion. By decreasing the sub-block width, the imaging performance of the presented method is improved. From Figure 10b, the sub-block width with 8 m satisfies the high performance imagery. In this case, we obtain the image which is mostly identical to that of the BP algorithm. However, the performance of the PCA method is inferior to that of the presented method and BP algorithm, because the residual error of the stop-and-hop approximation is not completely compensated. Inspecting Figures 8 and 9, the optimal sub-block width is about 8 m for three cases. In practice, the large sub-block width can be used at far range without loss of the imaging performance.

For targets at far range, the PSLR and ISLR are listed in Table 4. When the sub-block width is large, the error introduced by the sub-block processing method noticeably degrades the imaging performance of the presented method. The focusing performance of T15, T16, T17 and T18 also enhances this conclusion. The result of the presented method can be improved by decreasing the sub-block width. With the presented method, we obtain the high resolution results when the sub-block width is less than 8 m. The PSLR and ISLR related to T13 and T14 further strength this conclusion. However, the performance of the PCA method is greatly affected by the residual error of the stop-and-hop approximation at far range.


**Table 4.** Quality parameters for targets at far range.

Since the coupling term is space variant, the sub-block processing method is exploited to perform the decoupling operation. However, the sub-block method introduces the residual phase error, which is expressed as |*ϕde*\_*i*(*fτ*, *ft*;*r*ref\_*<sup>n</sup>* ± Δ*r*/2) − *ϕde*\_*i*(*fτ*, *ft*;*r*ref\_*n*)|. Here, Δ*r* denotes the sub-block width. Generally speaking, the imaging performance of the presented method highly depends on the sub-block width. In each sub-block, the residual phase error introduced by the decoupling operation should be limited within *π*/4 [31]. Under this condition, the influence of the phase error can be neglected.

Figures 8a, 9a and 10f are based on the sub-block width with 20 m. We find that the slices of the presented method are inferior to the slices of the PCA algorithm and BP method. The reason behind this is that the residual coupling error is not completely compensated. Hence, there is still the coupling between the range and azimuth dimensions. T1, T7 and T18 are far away from the sub-block center. The corresponding sub-blocks have a large width. Based on the decoupling phase of reference targets at the sub-block center, the decoupling operation across the whole sub-block is carried out. When the targets are at the sub-block edge, this operation introduces large residual coupling error, which is not limited within *π*/4. Due to this reason, the azimuth focusing performance such as the azimuth slice, PSLR and ISLR is seriously degraded. The focusing results of T1, T7 and T18 enhance this conclusion. When the targets are close to the sub-block center, the residual coupling error does not noticeably

affect the SAS focusing performance. Under this condition, the residual coupling error limited within *π*/4 can be negligible. Considering the trade-off between the imaging efficiency and performance, the optimal width of the sub-block is often exploited. In this case, the imaging performance of the presented method is nearly close to that of the BP method. The processing results of T5, T11 and T14 enhance this conclusion. Generally, the presented method can obtain the high resolution image across the whole swath based on the optimal width of the sub-block. Besides, it is very suitable for the SAS imagery with wide swath.

#### 6.1.3. Imaging Performance at Scenario Edge

Since the error of the stop-and-hop approximation is space variant, it is difficult to compensate this error using traditional methods. Fortunately, the presented method successfully solves this problem based on the accurate time delay of the signal. We now concentrate on the focusing performance at the scenario edge. Figure 11 shows the imaging scenario including three ideal targets.

**Figure 11.** Simulated scenario with three-point targets.

Based on the presented method, the PCA method and BP algorithm, the azimuth slices of focused targets are shown in Figure 12.

**Figure 12.** Azimuth slices of T19, T20 and T21. (**a**) T19; (**b**) T20; (**c**) T21.

With the presented method and PCA method, the target coordinates in Figure 12a are mostly accordant with coordinates shown in Figure 11. Considering the PCA method, there is a slight deviation of the azimuth coordinate in Figure 12b. Generally, this deviation can be negligible at close range. Inspecting Figure 12c, the PCA method suffers from noticeable deviation of the azimuth coordinate, and the deviation is about 0.01 m. Since the PCA method does not completely compensate the error of the stop-and-hop approximation, the residual error is introduced. It increases with the range and slow time. Due to this reason, the deviation can be negligible at close range. However, the deviation leads to the distortion when there are distributed targets at far range. Using the presented method, the targets across the whole swath are well focused.

#### *6.2. Real Data Processing*

We tested the presented method based on the real data. The data has 4800 sampling points in the range dimension and 3200 spatial sampling points in the azimuth dimension. For the transmitted signal, the center frequency and bandwidth are 150 kHz and 20 kHz, respectively. The receiver array, including 40 uniformly spaced receivers, is 1.6 m in azimuth. The velocity of the sonar platform is 2.5 m/s. When it comes to the operation of the differential decoupling with the presented method, two cases including four sub-blocks and eight sub-blocks in the range dimension are considered. The corresponding results are shown in Figure 13a,b, respectively. From Figure 13, it can be seen that the processing results of both cases are mostly identical. Therefore, we can draw a conclusion that the requirement of the sub-block segmentation in the range dimension can be relaxed. In practice, we can use large sub-block width without loss of performance when the real data is processed based on the presented method. For comparison, the real data is still processed by the BP algorithm. Figure 14 shows the processing result of the BP algorithm. Inspecting Figures 13 and 14, the presented method provide the high resolution result, which is mostly identical to that of the BP algorithm.

**Figure 13.** Processing results of the presented method. (**a**) Four sub-blocks; (**b**) eight sub-blocks.

**Figure 14.** Processing results of the BP algorithm.

The processing time with both methods are listed in Table 5.


**Table 5.** Processing time of imaging methods.

Both algorithms are developed based on Matlab 2012a. The processing time of the BP algorithm based on the sinc interpolation is 35,234 s. There are two schemes to develop the program of the presented method. With the first scheme, the calculation of the PSP and azimuth PSP are integrated with the imaging algorithm. It is time consuming due to the numerical evaluation of PSPs. In practice, the numerical calculation of PSPs can be carried out ahead of the focusing. With stored PSPs, we run the imaging algorithm. This is the second scheme. With this scheme, the processing time with four sub-blocks is decreased to 608 s. With the second scheme, the efficiency is dramatically improved compared with the first scheme. At this point, the second scheme is used with the real data processing. From Table 5, the presented algorithm with eight sub-blocks costs 1149 s. Overall, the processing time of the presented method increases with sub-blocks in range, because more time is needed to perform the differential decoupling. In comparison with the BP algorithm, the efficiency of presented method has been improved 30.7 times at least. Nowadays, the parallel algorithm and tools such as graphics processing unit (GPU), faster Fourier transform in the west (FFTW) and intel@ math kernel library (MKL) can be used to improve the efficiency of the presented method dramatically. Our future work is to optimize our method using parallel algorithms.
