A Novel Imaging Algorithm for Wide-Beam Multiple-Receiver Synthetic Aperture Sonar Systems

Abstract

In existing imaging algorithms for wide-beam multiple-receiver synthetic aperture sonar (SAS) systems, the double-square-root (DSR) range history of each receiver is generally converted into the sum of a single-square-root (SSR) range history and an error term using displaced phase center aperture (DPCA) approximation. Therefore, before imaging, each receiver’s error term needs to be individually compensated in the azimuth frequency domain, which is computationally expensive. As a result, a novel wide-beam multiple-receiver SAS system algorithm with low complexity and high precision is suggested. First, the translation relationship between the range histories of the reference receiver and other receivers is used to derive an SSR approximation range history that takes into account the azimuth variance of the non-stop-hop-stop time while ignoring differential range curvature (DRC) between the range histories from different receivers. Then, using the principle of stationary phase (POSP), the two-dimensional (2-D) spectrum of the point target is obtained. Finally, the multiple-receiver data are transformed into monostatic SAS-equivalent data for imaging after phase correction, time delay correction, and azimuth reconstruction. The range-Doppler (RD) algorithm is used as an example to explain the specific steps of the proposed approach. Simulation data and ChinSAS data experiments verify that the proposed algorithm achieves an imaging performance that is comparable to that of the existing wide-beam algorithm, but with much higher computational efficiency, making it suitable for real-time imaging.

1. Introduction

Synthetic aperture sonar (SAS) [1,2,3,4,5,6,7] utilizes a small-sized aperture that moves uniformly in a straight line in space to synthesize a large-sized aperture, which can obtain high-azimuth-resolution images independent of range and wavelength. The azimuth movement of the platform during pulse transmission and reception cannot be neglected since the SAS speed is comparable to the speed of sound [8]. The situation where the stop-hop-stop assumption is not valid is called the non-stop-hop-stop mode [9,10], and the round-trip propagation time of the signal is called the non-stop-hop-stop time [11,12].

The range history in the non-stop-hop-stop mode is a complex double-square-root (DSR) form concerning the non-stop-hop-stop time. Although Callow [11] obtained an analytical expression of the non-stop-hop-stop time of each receiver, it is challenging to employ it directly in the development of a successful imaging method. The existing line-by-line imaging algorithms always simplify the DSR range history to obtain the two-dimensional (2-D) spectrum by approximating the non-stop-hop-stop time. The non-stop-hop-stop time may be seen as the signal propagation time in the beam center in the narrow-beam case, but this approximation is inaccurate in the wide-beam case [10]. Although a time-domain beamformer can reconstruct SAS images with an arbitrary beam width [13], this requires expensive processing resources [14,15]. The main challenge for line-by-line wide-beam imaging algorithms is to derive a sufficiently accurate range history and 2-D spectrum by approximating the non-stop-hop-stop time. For multiple-receiver SAS systems, there are generally three major methods to obtain the 2-D spectrum.

The first method is to obtain the numerical solution of the 2-D spectrum of each receiver individually [16]. This method may provide an accurate image in a wide-beam case, but it requires multiple iterations to solve the point of stationary phase (PSP) for each receiver individually, resulting in low computational efficiency compared with other analytical methods.

The second method is to mathematically transform the DSR range history into a form that is easy to solve for the 2-D spectrum. Loffeld’s bistatic formula (LBF) [17] expanded the transmitting phase history and the receiving phase history with the Taylor series, and then combined the two individual PSPs to obtain the approximate bistatic PSP. Based on the LBF, Zhang developed corresponding range-Doppler [18], chirp scaling (CS) [19], and WK [20] algorithms to achieve monostatic conversion and imaging. In [21], by defining a new variable called the half quasi-bistatic angle (HQBA), an analytic formula of the point target response in the spectral domain was developed for BiSAR. Furthermore, considering the variability of the baseline in a multiple-receiver SAS, Tian [10] solved the fourth-order polynomial equation concerning the HQBA to obtain the quasi-analytical 2-D spectrum of each receiver. With the method of series reversion (MSR) [22], the DSR range history is expanded into a power series of the azimuth time, and the 2-D spectrum is obtained through MSR [23,24]. The main advantage of this method is that the accuracy of the 2-D spectrum can be adjusted by the order of the power series. In [25], an indirect range-Doppler algorithm (IDRDA) that avoids a complicated Taylor series expansion of the 2-D spectrum and just requires an explicit PSP was proposed. The IDRDA simplifies the design of the imaging processor and avoids the substantial computation load caused by interpolation in the traditional RD algorithm. However, in the wide-beam case, especially at far ranges, using the beam-center approximation of the non-stop-hop-stop time may degrade the image quality of the IDRDA significantly, which is also a common drawback of all the SAS imaging algorithms mentioned above. Recently, Ning [26] expanded the analytical expression of the non-stop-hop-stop time into a Taylor series directly, and then obtained the 2-D spectrum of the wide-beam signal through MSR. However, due to the coupling between the baseline length of the receiver and Doppler, fusion imaging is required after preprocessing each receiver individually, which greatly reduces the execution efficiency of the algorithm.

The third method is to process the multiple-receiver SAS data into monostatic SAS equivalents [19] using displaced phase center aperture (DPCA) technology [27,28]; then, the monostatic SAS imaging algorithms can be employed directly. The DSR range history can be converted into the sum of a single-square-root (SSR) range history and an error term using DPCA. Bonifant [27] first compensated the stop-hop-stop approximation by multiplying a fixed phase concerning the center of the swath. Subsequently, Yang [9] proposed a modified DPCA (MDPCA) to approximate the round-trip distance of the signal to twice the nearest range of the target, i.e., the non-hop-stop-hop time is regarded as the signal propagation time in the center of the transmitting beam. Zhang [29] compared the approximation errors in [9,27], and concluded that the MDPCA would provide a more accurate image than the algorithm in [27] for the narrow-beam case, while these two algorithms exhibited significant performance degradation in the wide-beam case. To improve the imaging performance of the MDPCA on wide-band signals, [30,31] additionally compensated for the point target reference spectral error and the residual quadratic coupling error in the RD and CS algorithms, respectively. However, this does not improve the imaging performance of the MDPCA in the wide-beam case. Zhang [32] improved the MDPCA algorithm utilizing the approximate relationship between the non-stop-hop-stop time and the azimuth angle within the beam, and obtained the approximate analytical solution for the 2-D spectrum. Tian [33] expanded the analytical expression of the non-stop-hop-stop time into a Taylor series, and then substituted it into the approximate DPCA range history to obtain the analytical two-dimensional spectrum. The algorithms in [32,33], to some extent, take into account the azimuth variance of the non-stop-hop-stop time, improving the imaging performance under wide-beam conditions. However, due to the coupling of the baseline length of receivers and Doppler frequency in the DPCA error term, the two algorithms both need to preprocess each receiver’s data separately, and then fuse them into a high-precision image; thus, the two algorithms are extremely inefficient.

To achieve low-complexity and high-precision imaging in the wide-beam case, we developed a novel imaging algorithm suitable for multiple-receiver SAS systems in this manuscript. First, we note that when the baseline length of the receiver is zero, the DSR range history will degenerate into an SSR form that considers the azimuth variance of the non-stop-hop-stop time. Therefore, according to the translation relationship of range histories between the reference receiver (baseline length is zero) and other receivers in the range–azimuth plane, and ignoring the differential range curvature [34] (DRC) between range histories of different receivers, the SSR range history of each receiver is derived. Then, the 2-D spectrum can be derived easily through POSP. Finally, phase correction, time delay correction, and azimuth reconstruction are performed, causing the undersampled data to be converted into a monostatic SAS equivalent that satisfies the sampling theorem, and the image is focused through the monostatic range-Doppler (RD) algorithm. The proposed algorithm considers the azimuth variance of the non-stop-hop-stop time and does not require a complex preprocess for each receiver individually used in existing algorithms; thus, it is suitable for wide-beam imaging, and has high computational efficiency. The imaging time for an actual area that is 200 m × 320 m is only one-tenth of that of existing wide-beam algorithms. The approximate error of the algorithm mainly comes from the DRC between different receivers, which increases with the array length and beam width. However, the experimental results in the following sections indicate that for existing wide-beam systems (such as HISAS 1030, HISAS 2040, muscle-SAS, etc.), the approximate error of the proposed algorithm fully meets the imaging requirements of less than one-eighth of the wavelength.

The remainder of this manuscript is organized as follows: In Section 2, the approximate range history based on relative range shifts is derived, the approximate error is analyzed, and the applicable conditions of the approximate range history are analyzed. In Section 3, the 2-D spectrum, the imaging process, and the computational efficiency of the proposed algorithm are shown. Section 4 compares the imaging performance of the proposed algorithm with the MDPCA algorithm, the wide-beam algorithm presented in [32], and the back projection (BP) algorithm, via simulation data and real data experiments. The experimental results verify the superiority of the proposed algorithm. Section 5 serves as a summary of this research.

2. The Approximate Range History Based on Relative Range Shift

2.1. The Exact Range History

The range history geometry of the multiple-receiver SAS is shown in Figure 1. The SAS platform moves uniformly in a straight line with velocity $v$ along the azimuth direction, and the point target $P$ is located at $\left(0,r\right)$. When the transmitter is located at $\left(vt,0\right)$, the signal is transmitted and $R_T\left(t,r\right)$ is the distance between the center of the transmitter and the point $P$. When the receiver moves to $\left(vt+d_i+v{\tau }_i^{*},0\right)$, the signal is received and $R_{R,i}\left(t,r\right)$ is the distance between the center of the receiver and the point $P$, where $d_i$ is the distance between the center of the transmitter and the center of the receiver $i$, $v{\tau }_i^{*}$ is the azimuthal movements of the SAS platform during pulse transmission and reception, and ${\tau }_i^{*}$ is the non-stop-hop-stop time.

The exact range history can be expressed as follows:

$$\begin{array}{ll}{R}_i^{*}\left(t,r\right)& =R_T\left(t,r\right)+R_{R,i}\left(t,r\right)\\ & =\sqrt{r^2+{\left(vt\right)}^2}+\sqrt{r^2+{\left(vt+v{\tau }_i^{*}+d_i\right)}^2}\end{array}$$

(1)

The exact range history can also be expressed as the product of the velocity of sound $c$ and the non-stop-hop-stop time ${\tau }^{*}$. Therefore, the equation can be expressed as follows:

$$\sqrt{r^2+{\left(vt\right)}^2}+\sqrt{r^2+{\left(vt+v{\tau }_i^{*}+d_i\right)}^2}=c{\tau }_i^{*}$$

(2)

By solving Equation (2), we obtain the following:

$${\tau }_i^{*}=\frac{v{d}_i+v^{2}t+c\sqrt{v^2{t}^2+r^2}}{c^2-v^2}+\frac{\sqrt{{\left(v{d}_i+v^{2}t+c\sqrt{v^2{t}^2+r^2}\right)}^2+\left(c^2-v^2\right)\left(2vt{d}_i+d_i^2\right)}}{c^2-v^2}$$

(3)

Thus, the exact range history [25] of the receiver $i$ is the following:

$$\begin{array}{ll}{R}_i^{*}\left(t,r\right)& =c{\tau }_i^{*}\\ & =\frac{v{d}_i+v^{2}t+c\sqrt{v^2{t}^2+r^2}}{c^2-v^2}c+\frac{\sqrt{{\left(v{d}_i+v^{2}t+c\sqrt{v^2{t}^2+r^2}\right)}^2+\left(c^2-v^2\right)\left(2vt{d}_i+d_i^2\right)}}{c^2-v^2}c\end{array}$$

(4)

where $c$ is the velocity of sound in water.

The exact range history in Equation (4) is difficult to directly use in the derivation of an effective imaging algorithm. A common approach is to transform the DSR range history in Equation (1) into a form that is easy to solve the 2-D spectrum, which requires approximating the non-stop-hop-stop time ${\tau }^{*}$ in Equation (3). The MDPCA [9] approximates the round-trip propagation distance of the signal to twice the nearest range of the target, i.e., the non-stop-hop-stop time is approximated to $2r/c$. This approximation is effective in the narrow-beam case, but the error increases in the wide-beam case, thereby reducing image quality. The algorithm in [32] approximates the round-trip propagation distance of the signal to $2r/\cos \theta$, where $\theta =\arctan \left[\left(vt+d_i/2\right)/r\right]$, as shown in Figure 1. This approximation improves the accuracy of the DSR range history in Equation (1) in the wide-beam case, thereby improving image quality. However, the algorithm in [32] requires a complex preprocess, resulting in low computation efficiency.

The MDPCA and the algorithm in [32] both adopt DPCA technology to transform the DSR range history into the SSR range history. Next, we will show how we obtain an SSR range history that is suitable for wide-beam imaging based on the relative range shift rather than on DPCA technology.

2.2. The Approximate Range History

The reference receiver is defined as the receiver where $d_i=0$. According to Equation (4), the range history of the reference receiver can be expressed as follows:

$$\begin{array}{ll}{R}_{ref}^{*}\left(t,r\right)& =\frac{2v^{2}t+2c\sqrt{v^2{t}^2+r^2}}{c^2-v^2}c\\ & \approx 2\beta ct+2\sqrt{v^2{t}^2+r^2}\end{array}$$

(5)

where $\beta =v^2/c^2$. Note that the range history of the reference receiver ($d_i=0$) shown in Equation (5) is a simple SSR form, rather than the complex DSR form of the range history of other receivers ($d_i\ne 0$) shown in Equation (4); therefore, we use the translation relationship of range histories between the reference receiver and other receivers to approximate the range histories of other receivers.

Figure 2 shows the range history diagram of the relative positions of the receiver $i$ and the reference receiver in the range–azimuth plane. The range history of the reference receiver is shown by the dashed line. When the reference receiver is at $t=t_{ref}$, the range history is the minimum value. If the derivative of $R_{ref}^{*}\left(t,r\right)$ with respect to $t$ in Equation (5) is 0, the azimuth time for obtaining the minimum value of $R_{ref}^{*}\left(t,r\right)$ is the following:

$$t_{ref}=-\frac{r}{c}$$

(6)

The non-stop-hop-stop times of different receivers at the minimum range history are different because there is a slight difference in the squint angle between the target and different receivers. However, due to this difference being very small, the non-stop-hop-stop time of different receivers at the minimum range history can be reasonably replaced by $-r/c$. Therefore, by shifting the reference receiver range history $R_{ref}^{*}\left(t,r\right)$ of Equation (5) to the left by $d_i/\left(2v\right)$ along the azimuth time axis, we can obtain the azimuth time for obtaining the minimum range history of the receiver $i$ as follows:

$$t_i=-\frac{r}{c}-\frac{d_i}{2v}$$

(7)

Substituting $t_i$ and $t_{ref}$ into Equations (4) and (5), respectively, the range translation amount between the receiver $i$ and the reference receiver is obtained is obtained:

$$\begin{array}{l}\Delta R_{offset,i}\left(r\right)=R_i^{*}\left(t_i,r\right)-R_{ref}^{*}\left(t_{ref},r\right)\\ =\frac{\frac{v{d}_i}{2}-\frac{v^{2}r}{c}+c\sqrt{v^2{\left(\frac{r}{c}+\frac{d_i}{2v}\right)}^2+r^2}}{c^2-v^2}c\\ +\frac{\sqrt{{\left(\frac{v{d}_i}{2}-\frac{v^{2}r}{c}+c\sqrt{v^2{\left(\frac{r}{c}+\frac{d_i}{2v}\right)}^2+r^2}\right)}^2-\left(c^2-v^2\right)\frac{2v{d}_{i}r}{c}}}{c^2-v^2}c+2\beta r-2r\sqrt{\beta +1}\end{array}$$

(8)

Considering that $\Delta R_{offset,i}\left(r\right)$ is too complex, it will not be expanded in the following section as it is unnecessary.

Figure 3 shows the schematic diagram of the translated range history of the receiver $i$, where $\delta$ is the DRC. After translation, the lowest point of the range history between the receiver $i$ and the reference receiver is aligned at $t_i$. Ignoring the DRC between other receivers and the reference receiver, the approximate range history of the receiver $i$ after translation is as follows:

$$R_i\left(t,r\right)=2\beta c\left(t+\frac{d_i}{2v}\right)+2\sqrt{r^2+v^2{\left(t+\frac{d_i}{2v}\right)}^2}+\Delta R_{offset,i}\left(r\right)$$

(9)

where $\Delta R_{offset,i}\left(r\right)$ is shown in Equation (8).

2.3. The Approximate Error of The Range History

Next, we analyze the approximate error under both narrow- and wide-beam situations, and compare it with the range history approximation error of the MDPCA algorithm, and the wide-beam algorithm presented in [32]. For convenience, we henceforth refer to the algorithm in [32] as the wide-beam algorithm based on DPC approximate slant range (DASR-WBA).

The approximate error $\Delta R_i\left(t,r\right)$ of the range history in Equation (9) can be expressed as the following:

$$\Delta R_i\left(t,r\right)=R_i^{*}\left(t,r\right)-R_i\left(t,r\right)$$

(10)

where $R_i^{*}\left(t,r\right)$ and $R_i(t,r)$ are shown in Equations (4) and (9), respectively. This approximate error mainly comes from ignoring the DRC between the reference receiver and other receivers.

The range history of the MDPCA [9] algorithm is the following:

$$\begin{array}{ll}{R}_{\mathrm{MDPCA},i}(t,r)& =R_T\left(t,r\right)+R_{R,i}\left(t,r\right)\\ & =\sqrt{r^2+{\left(vt\right)}^2}+\sqrt{r^2+{\left(vt+v{\tau }_i^{*}+d_i\right)}^2}\\ & \approx 2\sqrt{r^2+{\left(vt+v\underset{\raisebox{1ex}{${\tau }_i^{*}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\underbrace{\frac{r}{c}}}+\frac{d_i}{2}\right)}^2}+\underset{\mathrm{DPCA}\mathrm{error}\mathrm{term}}{\underbrace{\frac{{\left(v\frac{2r}{c}+d_i\right)}^2}{4r}}}\end{array}$$

(11)

Note that the non-stop-hop-stop time ${\tau }_i^{*}$ in Equation (11) is approximated to be $2r/c$ without considering the squint angle of the target; thus, the range history in Equation (11) is only suitable for the narrow-beam case. The narrow-beam assumption is the main source of the approximation error for the MDPCA algorithm. By combining Equations (4) and (11), we obtain the approximation error $\Delta R_{\mathrm{MDPCA},i}\left(t,r\right)$ of the MDPCA algorithm:

$$\Delta R_{\mathrm{MDPCA},i}\left(t,r\right)=R_i^{*}\left(t,r\right)-R_{\mathrm{MDPCA},i}\left(t,r\right)$$

(12)

The range history of the DASR-WBA [32] is as follows:

$$\begin{array}{ll}{R}_{\mathrm{DASR}-\mathrm{WBA},i}\left(t,r\right)& =R_T\left(t,r\right)+R_{R,i}\left(t,r\right)\\ & =\sqrt{r^2+{\left(vt\right)}^2}+\sqrt{r^2+{\left(vt+v{\tau }_i^{*}+d_i\right)}^2}\\ & \approx 2\sqrt{r^2+{\left(vt+v\underset{\raisebox{1ex}{${\tau }_i^{*}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\underbrace{\frac{r}{c\cdot \cos \theta }}}+\frac{d_i}{2}\right)}^2}+\underset{\mathrm{DPCA}\mathrm{error}\mathrm{term}}{\underbrace{\frac{{\left(v\frac{2r}{c\cdot \cos \theta }+\stackrel{\mathrm{coupling}\mathrm{term}}{\overbrace{d_i\cdot \cos \theta }}\right)}^2}{4r}}}\end{array}$$

(13)

where $\theta$ is the target squint angle.

Compared to the MDPCA algorithm, the improvement in the range history in Equation (13) is that the non-stop-hop-stop time is represented as $2r/\left(c\cdot \cos \theta \right)$. To some extent, this approximation makes the range history suitable for the wide-beam case. Note that the baseline length of the receivers and Doppler frequency are coupled together in the DPCA error term in Equation (13). Therefore, the DPCA error term of each receiver requires individual compensation in the azimuth frequency domain, and then fuses the preprocessing data of each receiver to provide a high-precision image, causing a lower computational efficiency compared with that of the MDPCA and the proposed algorithm. The computational efficiencies of the MDPCA [9], the DASR-WBA [32], and the proposed algorithm are discussed in Section 4.

By combining Equations (4) and (13), we obtain the approximation error $\Delta R_{\mathrm{DASR}-\mathrm{WBA},i}\left(t,r\right)$ of the DASR-WBA

$$\Delta R_{\mathrm{DASR}-\mathrm{WBA},i}\left(t,r\right)=R_i^{*}\left(t,r\right)-R_{\mathrm{WB},i}\left(t,r\right)$$

(14)

Typical system parameters are shown in

Table 1. System parameters.

Parameters	Beam Situations
	Narrow Beam	Wide Beam
Carrier Frequency, kHz	150	100
Bandwidth, kHz	20	20
PRI, s	0.4	0.4
Transmitter Length, m	0.08	0.04
Receiver Length, m	0.04	0.02
Receiver Number	48	96
Beam Width, °	6.3	19
SAS Velocity, m/s	2.4	2.4

. To study the approximate error over long distances and arrays, we appropriately increased the PRI and array length. The approximate error of the edge receiver is the maximum value; thus, we selected the edge receiver for analysis. The approximate errors in Equations (12), (14) and (10) are shown in Figure 4.

As can be seen from Figure 4a, in the case of the narrow beam, the approximate error of the MDPCA algorithm increases rapidly as the azimuth increases, but the maximum approximate error still meets the imaging requirements of less than one-eighth of the wavelength. In Figure 4b,c, the range histories of the DASR-WBA [32] and the proposed algorithm take into account the azimuth variance of the non-stop-hop-stop time; therefore, their approximate errors are smaller than that of the MDPCA algorithm, and the maximum value is far less than one-eighth of the wavelength. Generally speaking, the three algorithms all meet the imaging requirements in the case of the narrow beam. In the case of the wide beam, Figure 4d shows that the maximum approximate error of the MDPCA algorithm exceeds one-eighth of the wavelength, which does not meet the imaging requirements. Although the approximate error of the algorithm in Figure 4f is slightly larger than that of the DASR-WBA in Figure 4e, the approximate errors of the two algorithms are far less than one-eighth of the wavelength, which meets the imaging requirements. The results of Figure 4f verify the effectiveness of the proposed range history in the case of the wide beam. In addition, because the DASR-WBA needs to preprocess the data of each receiver and then fuse them for imaging, the computational efficiency is much lower than that of the proposed algorithm, which is discussed in Section 4.

The approximate error shown in Equation (10) mainly comes from the DRC between other receivers and the reference receiver, which increases as the array length and beam width increase. To study the applicable conditions for the range history shown in Equation (9), we analyzed the approximate error shown in Equation (10) under different array lengths and beam widths, as shown in Figure 5. When the beam width is 19°, the system parameters are the same as those shown in Table 1. When the beam width is 25°, the transmitter length is 0.03 m, the receiver length is 0.015 m, and the other parameters are the same as those in Table 1.

From Figure 5a, it can be seen that the approximate errors for all array length values (less than 4.92 m) meet the imaging requirements of less than one-eighth of the wavelength. Figure 4b shows that, when the array lengths are 1.92 m, 2.88 m, and 3.96 m, the approximate errors are less than one-eighth of the wavelength. When the array length increases to 4.92 m, the approximate error exceeds one-eighth of the wavelength at the near range, but still meets the imaging requirements at the far range. The beam width of existing wide-beam systems (such as HISAS 1030, HISAS 2040, muscle-SAS, etc.) is usually less than 25°, and the array length is far less than the maximum allowable array length of 3.96 m shown in Figure 5. Therefore, the approximate range history shown in Equation (9) can fully meet the requirements of wide-beam imaging.

3. The Wide-Beam Imaging Algorithm

3.1. Two-Dimensional Spectrum

In this section, we take the RD algorithm as an example to provide detailed steps for the proposed imaging algorithm. Assuming that the carrier frequency of the transmitting signal is $f_0$, the frequency modulation slope is $K_r$, the range fast time is $\tau$, the azimuth slow time is $t$, the pulse envelope is ${\omega }_r\left(\tau \right)$, and the beampattern determined by the receiver $i$ and the transmitter is ${\omega }_a\left(t\right)$, the baseband signal of the receiver $i$ in the 2-D time domain can be expressed as follows:

$$ss\left(\tau ,t,r\right)={\omega }_r\left[\tau -\frac{R_i^{}\left(t,r\right)}{c}\right]{\omega }_a\left(t\right)\exp \left[-j2\pi f_0\frac{R_i^{}\left(t,r\right)}{c}\right]\exp \left\{j\pi K_r{\left[\tau -\frac{R_i^{}\left(t,r\right)}{c}\right]}^2\right\}$$

(15)

By performing the range Fourier transform (FT) on Equation (15), we then can obtain the following:

$$Ss\left(f_r,t,r\right)={\displaystyle {\int }_{-\infty }^{\infty }ss\left(\tau ,t,r\right)\exp \left(-j2\pi f_r\tau \right)d\tau }$$

(16)

where $f_r$ is the range frequency.

The phase in the integral in Equation (16) is as follows:

$$\theta \left(\tau \right)=-\frac{2\pi f_0{R}_i\left(t,r\right)}{c}+\pi K_r{\left[\tau -\frac{R_i\left(t,r\right)}{c}\right]}^2-2\pi f_r\tau$$

(17)

According to POSP, by making the derivative of $\theta \left(\tau \right)$ zero with respect to $\tau$, we obtain the following:

$$\tau =\frac{f_r}{K_r}+\frac{R_i\left(t,r\right)}{c}$$

(18)

Substituting Equation (18) into Equation (17), we obtain the integral result of Equation (16)

$$Ss\left(f_r,t,r\right)=W_r\left(f_r\right){\omega }_a\left(t\right)\exp \left[-j\frac{2\pi \left(f_0+f_r\right)}{c}{R}_i\left(t,r\right)\right]\exp \left(-j\frac{\pi f_r{}^2}{K_r}\right)$$

(19)

where $W_r\left(f_r\right)$ is the range spectrum envelope.

Then, we perform the azimuth FT on Equation (19)

$$SS\left(f_r,f_a,r\right)={\displaystyle {\int }_{-\infty }^{\infty }Ss\left(f_r,t,r\right)\exp \left(-j2\pi f_{a}t\right)dt}$$

(20)

where $f_a$ is the azimuth frequency.

Substituting Equation (19) into Equation (20), the phase in the integral in Equation (20) is as follows:

$$\theta \left(t\right)=-\frac{2\pi \left(f_0+f_r\right)}{c}{R}_i\left(t,r\right)-\frac{\pi f_r^2}{K_r}-2\pi f_{a}t$$

(21)

Next, by substituting the range history shown in Equation (9) into Equation (21), and taking the derivative of $\theta \left(t\right)$ with respect to $t$, we obtain the following:

$$\begin{array}{ll}\frac{\mathrm{d}}{\mathrm{d}t}\theta \left(t\right)& =-\frac{2\pi \left(f_0+f_r\right)}{c}\frac{\mathrm{d}{R}_i\left(t,r\right)}{\mathrm{d}t}-2\pi f_a\\ & =-\frac{2\pi \left(f_0+f_r\right)}{c}\left[2\beta c+\frac{2v^2\left(t+\frac{d_i}{2v}\right)}{\sqrt{r^2+v^2{\left(t+\frac{d_i}{2v}\right)}^2}}\right]-2\pi f_a\end{array}$$

(22)

Since $\Delta R_{offset,i}\left(r\right)$ in Equation (9) is independent of $t$, the derivative of $\Delta R_{offset,i}\left(r\right)$ with respect to $t$ is zero. According to POSP, if Equation (22) is zero, then the solution is as follows:

$$t=-\frac{cr\left(2f_0\beta +2f_r\beta +f_a\right)}{v\sqrt{4v^2{\left(f_0+f_r\right)}^2-c^2{\left(2f_0\beta +2f_r\beta +f_a\right)}^2}}-\frac{d_i}{2v}$$

(23)

Substituting Equation (23) into Equation (21), we obtain the 2-D spectrum shown in Equation (20)

$$SS\left(f_r,f_a,r\right)=W_r\left(f_r\right)W_a\left(f_a\right)\exp \left[j{\varphi }_i\left(f_r,f_a,r\right)\right]$$

(24)

where $W_a\left(f_a\right)$ is the azimuth spectrum envelope, and

$$\begin{array}{ll}{\varphi }_i\left(f_r,f_a,r\right)=& -\frac{4\pi r{f}_0}{c}\sqrt{1-\beta }\sqrt{\frac{f_r^2}{f_0^2}+\frac{2f_r}{f_0}B\left(f_a\right)+D{\left(f_a\right)}^2}-\frac{\pi f_r^2}{K_r}\\ & -\frac{2\pi f_0}{c}\Delta R_{offset,i}\left(r\right)-\frac{2\pi f_r}{c}\Delta R_{offset,i}\left(r\right)+2\pi f_a\frac{d_i}{2v}\end{array}$$

(25)

where $\beta =v^2/c^2$, $B\left(f_a\right)=1-\frac{f_a}{2f_0\left(1-\beta \right)}$, and $D\left(f_a\right)=\sqrt{1-\frac{f_a\left(4\beta f_0+f_a\right)}{4\beta \left(1-\beta \right)f_0^2}}$.

Then, we expand the square root in Equation (25) into a power series of $f_r$, and retain it as the quadratic term to obtain the following:

$$\begin{array}{ll}{\varphi }_i\left(f_r,f_a,r\right)\approx & -\frac{2\pi r\sqrt{1-\beta }\left[D{\left(f_a\right)}^2-B{\left(f_a\right)}^2\right]}{cD{\left(f_a\right)}^3{f}_0^{}}{f}_r^2-\frac{\pi f_r^2}{K_r}-\frac{4\pi r\sqrt{1-\beta }B\left(f_a\right)}{cD\left(f_a\right)}{f}_r\\ & -\frac{4\pi r{f}_{0}D\left(f_a\right)\sqrt{1-\beta }}{c}-\frac{2\pi f_0}{c}\Delta R_{offset,i}\left(r\right)-\frac{2\pi f_r}{c}\Delta R_{offset,i}\left(r\right)+\pi f_a\frac{d_i}{v}\end{array}$$

(26)

In Equation (26), the first two terms are range compression terms, which can be compensated for in the 2-D frequency domain. The third term is the range migration term, which is compensated for through interpolation in the RD domain. The fourth term is the azimuth compression term, which can be compensated for in the RD domain. The fifth term and the sixth term are the phase correction term and the time delay correction term, respectively, which are compensated for in the range time domain and range frequency domain, respectively. The time delay correction term is weakly range-dependent; thus, it can be uniformly compensated for using the time delay correction term at the center of the swath in the range frequency domain. The seventh term is the azimuth reconstruction term, which can be compensated for by arranging the data of each receiver in order in the azimuth direction.

3.2. Imaging Process

The detailed flow of the proposed algorithm is shown in Figure 6. Firstly, phase correction and time delay correction are performed to compensate for the fifth and sixth terms in Equation (26), respectively. Phase correction can be performed in the 2-D time domain. Due to $\Delta R_{offset,i}\left(r\right)$ being weakly range-dependent, the time delay correction can be performed in the range frequency domain using the time delay at the center of the swath. The phase correction and time delay correction functions are, respectively, as follows:

$$H_{\mathit{phase},i}\left(r\right)=\exp \left[j\frac{2\pi f_0}{c}\Delta R_{\mathit{offset},i}\left(r\right)\right]$$

(27)

$$H_{\mathit{delay},i}\left(f_r\right)=\exp \left[j\frac{2\pi f_r}{c}\Delta R_{\mathit{offset},i}\left(r_{\mathit{ref}}\right)\right]$$

(28)

After correction, the data of each receiver are arranged in order along the azimuth time axis with the interval of $d_i/\left(2v\right)$ to complete the azimuth reconstruction. Then, the range compression in the 2-D frequency domain is performed. The range compression function is as follows:

$$H_{rc}\left(f_r,f_a,r_{ref}\right)=\exp \left\{j\frac{2\pi r_{ref}\sqrt{1-\beta }\left[D{\left(f_a\right)}^2-B{\left(f_a\right)}^2\right]}{cD{\left(f_a\right)}^3{f}_0^{}}{f}_r^2\right\}\exp \left(j\frac{\pi f_r^2}{K_r^2}\right)$$

(29)

The signal after range compression can be expressed as the following:

$$SS\left(f_r,t,r\right)=W_r\left(f_r\right)W_a\left(f_a\right)\exp \left[-j\frac{4\pi r\sqrt{1-\beta }B\left(f_a\right)}{cD\left(f_a\right)}{f}_r\right]\exp \left[-j\frac{4\pi r{f}_{0}D\left(f_a\right)\sqrt{1-\beta }}{c}\right]$$

(30)

where the first phase term is the range migration correction term, and the second phase term is the azimuth compression term. The signal is transformed into the RD domain via inverse Fourier transform (IFT) for range migration correction and azimuth compression, and the azimuth compression function is as follows:

$$H_a\left(f_a,r\right)=\exp \left[j\frac{4\pi r{f}_{0}D\left(f_a\right)\sqrt{1-\beta }}{c}\right]$$

(31)

Finally, the azimuth IFT is performed to obtain the focused image.

After the time delay correction shown in Equation (28) is completed, the 2-D spectrum shown in Equation (26) is converted into an equivalent monostatic 2-D spectrum. Next, in addition to the RD algorithm introduced in this study, other equivalent monostatic algorithms (such as CS, WK, etc.) can also be applied to imaging, which are not elaborated on in this manuscript.

3.3. Computational Efficiency

In this section, we discuss the computational efficiency of the proposed algorithm, the MDPCA algorithm, and the DASR-WBA.

We analyze the computational efficiency of these algorithms from the perspectives of Fast Fourier Transform (FFT) and interpolation operations, as these are the main sources of the computation. Assuming that the number of receivers is $N_{re}$, the number of pulses is $N_p$, and the number of range samples of the data matrix is $N_r$, then the number of azimuth samples of the data matrix is $N_a=N_{re}\cdot N_p$, and the total number of samples of the data matrix is $N_{data}=N_{re}\cdot N_p\cdot N_r$.

Figure 6 shows that the proposed algorithm includes $2N_a$ range FFTs (or IFFTs), $2N_r$ azimuth FFTs (or IFFTs), and $N_{data}$ sinc interpolations. Since multiplication operations on computers require much more time than addition operations, we only focus on multiplication operations. The number of complex multiplications for one range FFT (or IFFT) is $\frac{N_r}{2}{\log }_2{N}_r$; the number of complex multiplications for one azimuth FFT (or IFFT) is $\frac{N_a}{2}{\log }_2{N}_a$; and the number of complex multiplications for one sinc interpolation is 8 (8-point sinc interpolation). Therefore, the total number of complex multiplications for the proposed algorithm is the following:

$$\begin{array}{ll}{N}_{cm}& =2N_a\frac{N_r}{2}{\log }_2{N}_r+2N_r\frac{N_a}{2}{\log }_2{N}_a+8N_a{N}_r\\ & =N_a{N}_r\left[{\log }_2\left(N_a{N}_r\right)+8\right]\end{array}$$

(32)

According to [9], the MDPCA algorithm has the same number of range FFTs (or IFFTs), azimuth FFTs (or IFFTs), and sinc interpolations compared to the proposed algorithm; thus, the total number of complex multiplications is also shown in Equation (32).

The DASR-WBA includes two procedures, which are preprocessing and the imaging process. In preprocessing, we extract the data of the receiver $i$ from the original data matrix first. The number of azimuth samples of the extracted data is $N_p$. Secondly, in order to satisfy the sampling theorem, we perform an azimuth FFT on the extracted data of the receiver $i$, and then extend the extracted data via azimuth spectrum extension. The number of azimuth samples of the extended data is $N_a$, which is the same as the number of azimuth samples of the original data. Thirdly, the range compression, the Doppler phase compensation, the micro-range migration correction, and the bistatic phase compensation are performed successively. Finally, the extended data of all of the receivers are coherently superposed in the RD domain. After the above preprocessing, the multiple-receiver data are converted into monostatic SAS-equivalent data, and then the RD algorithm is used for imaging. Therefore, the total number of complex multiplications of this algorithm is the following:

$$\begin{array}{ll}{N}_{cm}^{*}& =N_{re}\left(N_r\frac{N_p}{2}{\log }_2{N}_p+2N_a{N}_r{\log }_2{N}_r\right)+N_a{N}_r{\log }_2{N}_r+8N_a{N}_r+N_r\frac{N_a}{2}{\log }_2{N}_a\\ & =N_a{N}_r\left[\left(2N_{re}+1\right){\log }_2{N}_r+\frac{1}{2}{\log }_2\left(N_p{N}_a\right)+8\right]\end{array}$$

(33)

According to Equations (32) and (33), the ratio of the number of complex multiplications between the DASR-WBA and the proposed algorithm is given by the following

$$\begin{array}{ll}{R}_{cm}^{}& =\frac{N_{cm}^{*}}{N_{cm}}\\ & =\frac{N_a{N}_r\left[\left(2N_{re}+1\right){\log }_2{N}_r+\frac{1}{2}{\log }_2\left(N_p{N}_a\right)+8\right]}{N_a{N}_r\left[{\log }_2\left(N_a{N}_r\right)+8\right]}\\ & =\frac{\left(2N_{re}+1\right){\log }_2{N}_r+\frac{1}{2}{\log }_2\left(N_p{N}_a\right)+8}{{\log }_2\left(N_a{N}_r\right)+8}\end{array}$$

(34)

$R_{cm}^{}$ can be used as a rough comparison of the calculation time between the proposed algorithm and the DASR-WBA. In general, $R_{cm}^{}$ is much greater than 1. For example, assuming $N_{re}$ is 48, $N_p$ is 50, and $N_r$ is 2048, then $R_{cm}^{}$ is 35.8, which shows that the number of complex multiplications of the DASR-WBA is 35.8 times that of the proposed algorithm in this research.

4. Experiments and Results

4.1. Simulation Experiments

In this section, the MDPCA algorithm, the DASR-WBA, the BP algorithm, and the proposed algorithm in this study were used to image in narrow-beam and wide-beam cases. The system parameters are shown in

Table 2. System parameters.

Parameters	Beam Situations
	Narrow Beam	Wide Beam
Carrier Frequency, kHz	150	100
Bandwidth, kHz	20	20
PRI, s	0.4	0.4
Transmitter Length, m	0.08	0.04
Receiver Length, m	0.04	0.02
Receiver Number	48	96
Beam Width, °	6.3	19
Pulse Number	50	100
SAS Velocity, m/s	2.4	2.4

. Two point targets were located at the near range and far range, which were $P_1\left(3\mathrm{m},80\mathrm{m}\right)$ and $P_2\left(0\mathrm{m},250\mathrm{m}\right)$, respectively. To ensure that the farthest target experienced a completely synthetic aperture length, the numbers of pulses in the case of the narrow beam and wide beam were set to 50 and 100, respectively. The receiver number and pulse number in the wide-beam case were twice as large as that in the narrow-beam case, resulting in the quantity of data in the wide-beam case being four times as much as that in the narrow-beam case. The computing operating system used was Microsoft Windows 11 (64 bit), the CPU was an Inter(R) Core (TM) quad-core [email protected] GHz, the memory was 16 GB, and the MATLAB version used was R2021a. Note that there was no weighting adopted for both the range and azimuth processing in the subsequent simulations and real data imaging.

Firstly, the imaging performance of the MDPCA algorithm, the DASR-WBA, the BP algorithm, and the proposed algorithm in the case of the narrow beam was studied. Figure 7 shows the azimuth profiles of the two point targets, $P_1\left(3\mathrm{m},80\mathrm{m}\right)$ and $P_2\left(0\mathrm{m},250\mathrm{m}\right)$.

From Figure 7, it can be seen that for the near-range target $P_1\left(3\mathrm{m},80\mathrm{m}\right)$, the azimuth profiles of the four algorithms are close. For the far-range target $P_2\left(0\mathrm{m},250\mathrm{m}\right)$, the sidelobe of the MDPCA algorithm is only slightly higher than those of other algorithms. Figure 7 indicates that the imaging degradation caused by the azimuth variance of the non-stop-hop-stop time in the narrow-beam case is slight; thus, the improvements of the proposed algorithm, the DASR-WBA, and the BP algorithm are small, which corresponds to the experimental results in Figure 4a–c in Section 2.

Next, we studied the imaging performance in the case of the wide beam. Figure 8 shows the azimuth profiles of $P_1\left(3\mathrm{m},80\mathrm{m}\right)$ and $P_2\left(0\mathrm{m},250\mathrm{m}\right)$.

Figure 8 shows that for the near-range target, the sidelobe of the proposed algorithm is about 2 dB lower than that of the MDPCA algorithm, which is close to those of the DASR-WBA and the BP algorithm. For the far-range target, the advantages of the proposed algorithm are more obvious. The MDPCA algorithm has a 4 dB increase in the sidelobe for the far-range target, which may result in ghost targets. The MDPCA algorithm shows worse performance in the wide-beam case compared to the proposed algorithm in this research, indicating that the imaging degradation caused by the azimuth variance of the non-stop-hop-stop time in the wide-beam case cannot be ignored; this is consistent with the experimental results of Figure 4d–f in Section 2.

The imaging performance of the proposed algorithm in this study is relatively similar to that of the DASR-WBA and the BP algorithm, regardless of whether it is under narrow- or wide-beam conditions. This is because the DASR-WBA and the BP algorithm also consider the azimuthal variance of the non-stop-hop-stop time. However, the DASR-WBA and the BP algorithm require complex processing, resulting in significant computational complexity.

To quantitatively analyze the imaging results of the four algorithms, we compared the imaging results through the impulse response width (IRW), peak sidelobe level ratio (PSLR), and integrated sidelobe level ratio (ISLR) of the azimuth profiles, as shown in

Table 3. Imaging parameters of simulated data.

Beam Situations	Algorithm	Calculation Time, s	Target	IRW, cm	PSLR, dB	ISLR, dB
Narrow beam	MDPCA	3.79	P1	4.83	−18.26	−15.83
			P2	4.78	−19.51	−19.42
	DASR-WBA	35.64	P1	4.83	−18.38	−15.86
			P2	4.78	−19.96	−19.43
	Proposed algorithm	4.17	P1	4.83	−18.38	−15.83
			P2	4.76	−19.97	−19.43
	BP	—	P1	4.83	−18.89	−15.86
			P2	4.54	−20.02	−20.64
Wide beam	MDPCA	16.90	P1	2.53	−15.89	−15.07
			P2	2.76	−14.77	−15.99
	DASR-WBA	311.23	P1	2.44	−18.14	−15.13
			P2	2.40	−19.58	−17.45
	Proposed algorithm	16.65	P1	2.44	−18.15	−15.14
			P2	2.60	−19.31	−17.41
	BP	—	P1	2.42	−17.86	−16.18
			P2	2.36	−19.03	−20.31

. Table 3 also provides a comparison of the calculation times. Since the computational complexity of the BP algorithm is much higher than those of other algorithms; thus, the calculation time of the BP algorithm is not listed. The bold font shows the parameters of the proposed algorithm in this study.

According to Table 3, we can see that in the case of the narrow beam, the IRW, PSLR, and ISLR of the four algorithms are very close for the near-range target (P₁). For the far-range target (P₂), the PSLRs of the proposed algorithm and the DASR-WBA are about 0.4 dB lower than the PSLR of the MDPCA algorithm. In the case of the wide beam, for the near-range target (P₁), the IRWs of the proposed algorithm, the DASR-WBA, and the BP algorithm are about 0.1 cm smaller than that of the MDPCA algorithm. The PSLRs of the proposed algorithm, the DASR-WBA and the BP algorithm are about 2 dB lower than that of the MDPCA algorithm. For the far-range target (P₂), compared with the MDPCA algorithm, the IRW of the proposed algorithm is reduced by 0.16 cm, the PSLR is reduced by 4.54 dB, and the ISLR is reduced by 1.42 dB. The IRW, PSLR, and ISLR of the proposed algorithm are close to those of the DASR-WBA and the BP algorithm, indicating that the imaging performance of these algorithms is comparable.

In terms of calculation times, the proposed algorithm is much shorter than that of the DASR-WBA, with a reduction of 89.4% and 94.6% in the narrow-beam and wide-beam situations, respectively. Although the DASR-WBA can reduce the calculation time through parallel computing, this requires complex hardware configurations for the SAS. According to Equation (34), in the narrow-beam and wide-beam situations, the calculation times of the proposed algorithm should be reduced by 97.2% and 98.5%, respectively. The main reason for the calculation time reduction being smaller than the theoretical values is the efficient matrix FFT operations in MATLAB and the time-consuming loop operations in the sinc interpolation.

The above simulation experimental results indicate that the algorithm proposed in this paper has good imaging performance in both narrow-beam and wide-beam cases, and the computational efficiency is much higher than that of the existing wide-beam imaging algorithm.

4.2. Real Data Experiments

In this section, we conducted experiments using the data collected from a sea trial in 2017 by ChinSAS, an interference SAS system of the Naval University of Engineering with a beam width of 6.3°. The MDPCA algorithm, the DASR-WBA, the proposed algorithm, and the BP algorithm were used to image a swath with an area of 200 m × 320 m. The imaging results of some of the areas are shown in Figure 9.

As can be seen from Figure 9, the imaging results of the four algorithms are similar, but there are some differences. To quantitatively analyze the imaging results of the four algorithms, we compared the imaging results of different algorithms using the azimuth profiles of the two point-like targets in Figure 9 (marked in yellow boxes), as well as the IRW, PSLR, and ISLR of the azimuth profiles, as shown in Figure 10 and

Table 4. Imaging parameters of algorithms.

Target	Algorithm	IRW, cm	PSLR, dB	ISLR, dB
1	MDPCA	16.64	−3.79	−2.28
	DASR-WBA	8.30	−6.96	−0.72
	Proposed algorithm	8.95	−6.09	−1.24
	BP	24.81	−8.50	−5.18
2	MDPCA	15.95	−6.56	0.95
	DASR-WBA	24.05	−10.63	−1.02
	Proposed algorithm	24.51	−11.85	−1.32
	BP	18.30	−12.96	−2.56

, respectively.

Figure 10 clearly shows the advantages of the proposed algorithm. The main lobe width of the proposed algorithm is narrower than that of the MDPCA algorithm, and the sidelobe height of the proposed algorithm is much lower than that of the MDPCA algorithm. The imaging result of the DASR-WBA is similar to that of the proposed algorithm, while the results of the BP algorithm are slightly better than those of other algorithms.

According to Table 4, it can be seen that for Target 1, the parameters of the proposed algorithm and the DASR-WBA are comparable, and the IRW and PSLR of both algorithms are superior to those of the MDPCA algorithm. Although the ISLR of the MDPCA algorithm is superior to those of the proposed algorithm and the DASR-WBA, this is due to the wide pulse width of the MDPCA algorithm. Additionally, the IRW of the BP algorithm is higher than those of other algorithms, while PSLR and ISLR are much lower than those of other algorithms.

For Target 2, the parameters of the proposed algorithm are relatively close to those of the DASR-WBA, and both are far superior to those of the MDPCA algorithm. Additionally, all parameters of the BP algorithm are superior to other algorithms. Compared with other algorithms, the parameters of the proposed algorithm are closer to those of the BP algorithm.

In addition, for an area with a size of 200 m × 320 m, the calculation times of the MDPCA algorithm, the DASR-WBA, and the proposed algorithm are 50.05 s, 517.70 s, and 55.82 s, respectively. The calculation time of the proposed algorithm is similar to that of the MDPCA algorithm, but much shorter than that of the DASR-WBA. The reason is that the proposed algorithm does not require complex preprocessing for each receiver. Due to the calculation complexity of the BP algorithm being much higher than those of other algorithms, the calculation time of the BP algorithm is not shown.

5. Conclusions

In this study, a novel low-complexity and high-precision imaging algorithm based on relative range shifts for wide-beam multiple-receiver SAS was proposed. The algorithm obtained an SSR range history of each receiver by simply translating the range history of the reference receiver in the range–azimuth plane. The azimuth variance of the non-stop-hop-stop time is considered in the SSR range history, through which the 2-D spectrum of each receiver can be conveniently derived for imaging. The proposed algorithm may provide comparable imaging performance to other wide-beam algorithms, but with a much higher computational efficiency. Compared with other wide-beam algorithms, the proposed algorithm is more suitable for real-time imaging of wide-beam multiple-receiver SAS systems.