*Article* **Maritime Multiple Moving Target Detection Using Multiple-BDS-Based Radar: Doppler Phase Compensation and Resolution Improvement**

**Xiang Lan 1, Liuying Wang 2, Jinxing Li 1, Wangqiang Jiang <sup>1</sup> and Min Zhang 1,\***


**Abstract:** With the realization of global navigation satellite system (GNSS) completion, GNSS reflectometry (GNSS-R) has become increasingly popular due to the advantages of global coverage and the availability of multiple sources in terms of earth remote sensing. This paper analyzes the Beidou navigation satellite system (BDS) signal reflection detection of multiple satellites and multiple moving targets under multiple-input and multiple-output (MIMO) radar systems and proposes a series of methods to suppress multiple Doppler phase influences and improve the range detection property. The simulation results show the restored target peaks, which match the RCS data more accurately, with the GNSS-R Doppler phase influence removed, which proves the proposed method can improve target recognition and detection resolution performance.

**Keywords:** global navigation satellite system reflectometry (GNSS-R); Beidou navigation satellite system (BDS); Doppler compensation; range resolution

#### **1. Introduction**

At the moment more than 70 satellites are already in view; this brings great opportunities and challenges for both scientific and engineering applications. The different global navigation satellite system (GNSS) signals are compared and analyzed in the literature [1] in terms of detection performance and signal characteristics. A four-system positioning model for multi-satellite detection is proposed in the literature [2].

Global navigation satellite system reflectometry (GNSS-R) detection has become a popular tool for earth remote sensing, such as multiple maritime targets detection, with its gradual improvement in observation technology. Based on GNSS-R delay–Doppler map (DDM) imaging, an effective method of expression based on GNSS-R signals, the literature [3–5] shows an incoherent range walk compensation method and accurate realization methods to improve the DDM imaging effectiveness. The waveform and detection models were analyzed in [6,7] for GNSS-R detection. With the gradual maturity of the Beidou navigation satellite system (BDS), BDS satellite remote sensing studies have become increasingly widespread. The single-frequency PPP time transfer performance of BDS-2/3 is evaluated for excellent performance and BDS-3's is expected for better accuracy with the continuous development of real-time products of BDS [8]. Time–frequency-transfer and time performance technology based on BDS systems have been implemented and evaluated [9,10]. With the increase in the demand for precise positioning, an increasing number of observational reports have been proposed with constant improvement of the positioning method. The literature [11,12] presents a real-time detection and point position method based on the combination of GPS and BDS observations.

BDS-R detection still leaves some problems to be solved. Due to the low reflected signal power, it is difficult to obtain complete target information with a background of

**Citation:** Lan, X.; Wang, L.; Li, J.; Jiang, W.; Zhang, M. Maritime Multiple Moving Target Detection Using Multiple-BDS-Based Radar: Doppler Phase Compensation and Resolution Improvement. *Remote Sens.* **2021**, *13*, 4963. https://doi.org/ 10.3390/rs13244963

Academic Editors: Shuanggen Jin and Gino Dardanelli

Received: 8 November 2021 Accepted: 3 December 2021 Published: 7 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

strong clutter and noise. Reference [13] notes that for large and medium-sized targets on the sea surface, GNSS-R detection can acquire a sufficient signal-to-noise ratio (SNR) of 20~25 dB [14] after several seconds of accumulation. Next, one study in the literature [14,15] simulated a delay–Doppler map (DDM) of large-scale targets and verified the effectiveness of large-scale target detection based on short time accumulation signals using the range–Doppler Fourier method. In terms of theoretical construction and system model, reference [16] gives an extremum approximation algorithm for advanced receiver autonomous integrity monitoring, and triple-frequency combining observation models have been studied for precise point positioning [17]. Moreover, real-time, direction-constrained determination methods were proposed for point velocity detection studies [18–20]. In order to better interpret the target information of the GNSS-R signal and obtain the images that match the detected target, many signal processing methods were proposed. The incoherent range walk compensation method was proposed for spaceborne GNSS-R imaging [3], and GNSS-R-based moving target indication is studied in the literature [21]. In terms of marine target detection, ocean surface target detection and positioning were discussed [22], the delay and Doppler tracking errors were analyzed [15], and the feasibility analysis of ship detection by DDM was simulated [23]. In order to distinguish sea targets from sea clutter, the blind sea clutter suppression method was discussed [24]. Moreover, the data of TDS-1 were used to prove the effectiveness of the proposed target detection method [25], and the anomalous artifacts of TDS-1 DDM were analyzed in [26].

In this paper, multiple-BDS signal reflection detection for multiple moving targets was studied. We propose the multiple targets Doppler compensation (MTDC) method to keep target peaks from the influence of the Doppler phase and provide a target peak identification method and range inversion methods to realize the specific scatter information detection of targets. The contributions of this paper can be summarized as follows: (1) The target peak reduction problem existing in multiple-BDS signal detection was analyzed by echo formulas, and solutions are proposed for two situations: when target estimation information is acquired and when it is not. In the simulation, a difference of more than 400 Hz in the Doppler frequency was set among the detected targets, which further verifies that the proposed method is effective; (2) The target peak identification method is proposed based on 33 m sampling at the echoes, and the improved target detection peaks can express the target information more accurately in terms of amplitude and range point; (3) A range inversion method is proposed to acquire the specific scatter information in the case of low-resolution BDS signals. The results are in good agreement with the corresponding RCS data in the ship target simulation; (4) In terms of practicability, the proposed MTDC method 1 and range inverse method 1 can remove the Doppler phase influence and further acquire the targets' RCS distribution without target estimation. The target peak identification results based on 33 m echo sampling can be utilized to verify target recognition. Moreover, the results were simulated and analyzed under clutter and noise backgrounds.

The rest of the sections in this paper are organized as follows. Section 2 discusses the MIMO BDS signal detection model and lists existing problems. The Doppler compensation methods and range detection resolution design of the BDS signal based on MIMO multipletarget detection are discussed in Section 3. In Section 4, the simulation results and analysis are obtained. Finally, the paper is concluded in Section 5.

#### **2. Detection Model and Existing Problems**

#### *2.1. Detection Model*

The detection model of multiple GNSS-R signals for multiple targets was firstly outlined. The BDS signals, which are transmitted by MEO (medium earth orbit) satellites in the Beidou No. 3 system fixed at 2.1528 × 104 km above the Earth's surface, were used as detection sources. The detection simulated maritime multiple moving target detection under space-based observation. The geometry of the model is given in Figure 1, where a rectangular coordinate system is established and the XOY level coincides with the sea level. The satellites, targets, and receiver positions are marked in red.

**Figure 1.** MIMO and multiple target BDS signal detection geometric model.

#### *2.2. Signal Structure*

During the detection model, several sets of B1I and B3I signals transmitted by MEO satellites with corresponding pseudo-random noise (PRN) codes were used. A part of the signals reaches the airborne-based receiver directly, and the other part of the signals is reflected by the ship targets and arrives at the latter receiver. The received echo is the sum of all echoes transmitted by each satellite and reflected by each target. In the simulation experiment, we found the echo modulation information through the proposed geometric model with satellite ephemeris, receiver motion information and target position, altitude, and velocity information. The BDS signals' model structures are expressed as formulae with changes and states in the propagation and receivers, as follows.

In the formula expression, we set '*i*' as the satellite number and '*j*' as the target number. The transmitting BDS signal [27,28] of satellite *i* is expressed as:

$$s^i(t) = A\mathcal{C}^i D^i \cos(2\pi f\_0 t + q^i) \tag{1}$$

where *s<sup>i</sup>* (*t*) contains B1I and B3I signals. *A* is the signal amplitude, *C<sup>i</sup>* the signal ranging code, and *D<sup>i</sup>* the signal data code. During the cosine function, *f*<sup>0</sup> is the carrier frequency, and *ϕ<sup>i</sup>* is the original phase. The echoes after propagation and receiving can be expressed as follows:

$$R^i\_d = a^i\_d \cdot s^i(t - \tau^i\_{d'} f d^i\_d) \tag{2}$$

$$R\_r^{i,j} = a\_r^{i,j} \cdot s^i(t - \pi\_r^{i,j}, f d\_r^{i,j}) \tag{3}$$

where the subscripts '*d*' and '*r*' are direct signals and reflected signals, respectively. *a<sup>i</sup> <sup>d</sup>* is the propagation coefficient and scattering coefficient of direct signals, and *a i*,*j <sup>r</sup>* is the reflected signal reflected by target *j*. *τ* and *f d* are the responding transmitted delay and Doppler frequency, which are related to the carrier frequency and show different values between the B1I and B3I signals. Finally, the total received B1I and B3I signals transmitted by multiple satellites and propagated through the direct and multiple targets reflection approach is expressed as

$$R = \sum\_{i=1}^{M} (R\_{\text{B1I},d}^{i} + \sum\_{j=1}^{N} R\_{\text{B1I},r}^{i,j}) + \sum\_{i=1}^{M} (R\_{\text{B3I},d}^{i} + \sum\_{j=1}^{N} R\_{\text{B3I},r}^{i,j}) \tag{4}$$

For pulse accumulation and slow-time Doppler frequency detection, the echoes expressed in *Q* ranging code periods is (5), after data code demodulation processing.

$${}^{Q}\mathcal{R} = \sum\_{q=1}^{Q} \mu\_{Tp}(t - q \cdot T\_P) \cdot \mathcal{R} \tag{5}$$

where *QR* represents the echoes in *<sup>Q</sup>* ranging code periods, and *uTP* (*<sup>t</sup>* − *<sup>q</sup>* · *TP*) = 1, *q* · *TP* < *t* < (*q* + 1)*TP* is the window function. *TP* is the ranging code cycle period.

#### *2.3. Existing Problems*

Based on the detection model, there exist several problems which promote the development of the proposed methods:


#### **3. Methods**

#### *3.1. Multiple Targets Detection Doppler Compensation*

In order to remove the Doppler phase influence during multiple targets detection, we further studied echo processing. The B3I signal is set as an example to express the formula, and the B1I signal has the same form. Take Formulas (1) and (3) into (5) as

$$\begin{split} \mathcal{Q}\mathcal{R}\_{\text{B3I},\text{r}}^{i,j} &= \sum\_{q=1}^{\mathcal{Q}} u\_{T\_{\text{P}}}(t - q \ast T\_{\text{P}}) \cdot \mathcal{R}\_{\text{B3I},\text{r}}^{i,j} \\ &= a\_{\text{B3I},\text{r}}^{i,j} A\_{\text{B3I}}^{i} \mathcal{C}\_{\text{B3I}}^{i} \cdot \sum\_{q=1}^{\mathcal{Q}} \left\{ u \left( t - \tau\_{rj}^{i} - q \cdot T\_{\text{P}} \right) \\ &\cdot \cos(2\pi f d\_{\text{B3I},\text{r}}^{i,j} (t - \tau\_{\text{B1I},\text{r}}^{i,j} - q \cdot T\_{\text{P}}) + \varphi\_{\text{B3I}}^{i}) \right\} \end{split} \tag{6}$$

where *QR<sup>i</sup>*,*<sup>j</sup>* B3I,r is the received reflected B3I signals transmitted by satellite *i* and reflected by target *j* during *Q* ranging code periods. According to the ranging code *C<sup>i</sup>* B3I, we find the duty ratio of the signal to be 1 in a ranging code period. Therefore, the Doppler phase changes during a signal ranging code period cannot be ignored. The echoes model (6) is improved as

$$\begin{aligned} ^{QR}R\_{\rm B3I,r}^{i,j} &= a\_{\rm B3I,r}^{i,j}A\_{\rm B3I}^{i}\Gamma\_{\rm B3I}^{i} \cdot \sum\_{q=1}^{Q} \sum\_{n=1}^{N\_{\rm C}} \left\{ u\_{T\_{\rm c}} \left( t - \tau\_{\rm r}^{i,j} - qT\_{\rm P} - nT\_{\rm c} \right) \right. \\ &\cdot \cos(2\pi f d\_{\rm B3I,r}^{i,j} \left( t - \tau\_{\rm B3I,r}^{i,j} - qT\_{\rm P} - nT\_{\rm c} \right) + qt\_{\rm B3I}^{i}) \right\} \end{aligned} \tag{7}$$

where *Nc* is the number of range code elements, and *Tc* is the range code element period. In (7), the Doppler phase of echoes changes with each ranging code element, and it influences the target range peak result level of the matched filter. The normalized autocorrelation peak level changing with Doppler frequency is shown in Figure 2, where the peak level reduces slowly first and then drops quickly to lower than 0.1 when the Doppler frequency is more than 1000 Hz. Moreover, the relevant parameters are described in the simulation section. Focusing on the Doppler phase influence, we express (7) as

$$\begin{array}{c} \,^{Q}R\_{\text{B3I},\text{r}}^{i,j} = \,^{Q}R\_{\text{B3I},\text{r}}^{i,j}(\,^{h}h\_{\text{B3I},\text{r}}^{i,j})\\ ^{Q}Ph\_{\text{B3I},\text{r}}^{i,j} = \sum\_{q=1}^{Q} \sum\_{n=1}^{N\_{\text{C}}} \,^{n}\_{\text{r},\text{r}}(t-qT\_{P}-nT\_{\text{c}}) \cdot 2\pi f d\_{\text{B3I},\text{r}}^{i,j}(t-qT\_{P}-nT\_{\text{c}})\\ = \,^{Q}Ph\_{\text{B3I},\text{r}}^{i,j}(fd\_{\text{B3I},\text{r}}^{i,j}) \end{array} \tag{8}$$

where we can find that *QPh<sup>i</sup>*,*<sup>j</sup>* B3I,r is the sequence with a length of *<sup>Q</sup>* · *Nc* and is up to *f di*,*<sup>j</sup>* B3I,r. Moreover, each target reflects echo carriers with different Doppler frequencies, and we needed to reduce multiple Doppler phase influences for better detection of the target peak.

**Figure 2.** Target peak curves with the Doppler channel in the detection of BDS signals.

Before the Doppler phase compensation, the compensatory Doppler frequency of target echoes needed to be acquired. We propose two methods: method 1, without estimation and method 2, with estimation. In method 1, we acquired the Doppler frequency range of targets by satellite navigation code data and receiver data as

$$\begin{array}{l} \mathit{fd}\_{\mathrm{B1,min}}^{\mathrm{id}} = \\ \mathit{min}(\mathit{fd}\_{\mathrm{B1,sat\\_rec}}^{\mathrm{id}} + \mathit{fd}\_{\mathrm{B1l,rec\\_tar\\_max}^{\mathrm{id}}}^{\mathrm{id}} \mathit{fd}\_{\mathrm{B1l,sat\\_rec}}^{\mathrm{id}} - \mathit{fd}\_{\mathrm{B1l,rec\\_tar\\_max}^{\mathrm{id}}}^{\mathrm{id}}) \\ \mathit{fd}\_{\mathrm{B1l,max}}^{\mathrm{id}} = \\ \mathit{max}(\mathit{fd}\_{\mathrm{B1l,sat\\_rec}}^{\mathrm{id}} + \mathit{fd}\_{\mathrm{B1l,rec\\_tar\\_max}^{\mathrm{id}}}^{\mathrm{id}} \mathit{fd}\_{\mathrm{B1l,sat\\_rec}}^{\mathrm{id}} - \mathit{fd}\_{\mathrm{B1l,rec\\_tar\\_max}^{\mathrm{id}}}^{\mathrm{id}}) \end{array} \tag{9}$$

and the compensated Doppler frequencies are acquired by sampling in the range with an interval of 125 Hz as

$$\begin{aligned} \,^1f \, d\_{\text{B1I}, \text{r}}^{i, j'} &= \,^1f \, d\_{\text{B1I}, \text{min}}^i + 250 \ast (j' - 1) \\ \,^1f \, d\_{\text{B1I}, \text{min}}^i &\le \,^1f \, d\_{\text{B1I}, \text{r}}^{i, j'} \le \,^1f \, d\_{\text{B1I}, \text{max}}^i \end{aligned} \tag{10}$$

where <sup>1</sup> *f d<sup>i</sup>* B1I,*rj* is the compensation Doppler frequency of method 1 and *j* is the compensation number. Using the estimated information of targets, the compensated Doppler frequencies are set as Doppler frequency estimation values in method 2.

$$^2f d\_{\rm B1I,r}^{i,j'} = ^2f d\_{\rm B1I,e}^{i,j'} \tag{11}$$

where <sup>2</sup> *f di*,*<sup>j</sup>* B1I,r is the compensation Doppler frequency of method 2. After acquiring the corresponding Doppler frequency, the compensation phase is expressed as

$$\begin{split} ^{\mathcal{Q}}\mathbb{C}\Pr\_{\mathrm{B3I},\mathrm{r}}^{i,j'} &= \sum\_{q=1}^{\mathcal{Q}} \sum\_{n=1}^{N\_{\mathcal{C}}} u\_{T\_{\mathcal{C}}}(t - qT\_{\mathcal{P}} - nT\_{\mathcal{C}}) \cdot 2\pi f d\_{\mathrm{B3I},\mathrm{r}}^{i,j'}(t - qT\_{\mathcal{P}} - nT\_{\mathcal{C}}) \\ &= ^{\mathcal{Q}}\mathbb{C}\Pr\_{\mathrm{B3I},\mathrm{r}}^{i,j'}(f d\_{\mathrm{B3I},\mathrm{r}}^{i,j'}) \end{split} \tag{12}$$

To solve the multiple Doppler phase influences and restore the weakened target peak, we propose the multiple target Doppler compensation (MTDC) method, which compensates for the echoes with multiple Doppler phase sequences at the corresponding range code period position of targets and fetches the maximum value of the filter outputs as the compensation result. Next, (13) expresses the compensated echo.

$${}^{Q}\mathcal{C}R^{i,j,j'}\_{\text{B3I,r}} = {}^{Q}R^{i,j}\_{\text{B3I,r}} ({}^{Q}Ph^{i,j}\_{\text{B3I,r}} + {}^{Q}\mathcal{C}Ph^{i,j'}\_{\text{B3I,r}}) \tag{13}$$

where *QCR<sup>i</sup>*,*j*,*<sup>j</sup>* B3I,r denotes the echoes of satellite *i* and target *j* with compensation *j'.* The echoes from satellite *i* with multiple-target Doppler compensation can be expressed as

$${}^{\mathbb{Q}}\mathbb{C}\,R\_{\mathrm{B3I}}^{i,j'} = {}^{\mathbb{Q}}\mathbb{C}R\_{\mathrm{B3I},\mathrm{d}}^{i} + \sum\_{j=1}^{N} {}^{\mathbb{Q}}\mathbb{C}\,R\_{\mathrm{B3I},\mathrm{r}}^{i,j,j'} \tag{14}$$

Next, we matched echoes with the corresponding satellite signal as

$${}^{\mathbb{Q}}\mathbb{C}\,M\_{\mathrm{B1I}}^{i,j'} = {}^{\mathbb{Q}}\mathbb{C}R\_{\mathrm{B3I}}^{i,j'} \odot s\_{\mathrm{B3I}}^{i}(t) \tag{15}$$

Finally, (16) shows the selected maximum of each compensation as the result.

$${}^{\mathbb{Q}}\mathbb{C}M\_{\mathrm{B1I}}^{i}\mathbb{c} = \max({}^{\mathbb{Q}}\mathbb{C}M\_{\mathrm{B1I}'}^{i,1}, \dots, {}^{\mathbb{Q}}\mathbb{C}M\_{\mathrm{B1I}'}^{i,j'}, \dots, {}^{\mathbb{Q}}\mathbb{C}M\_{\mathrm{B1I}}^{i,N'})\tag{16}$$

The formula mentioned above is the one-dimensional range–Doppler compensation of multiple target reflections. In the next part, we show the Doppler and the delay two-dimensional Doppler compensation methods for Doppler delay image target peak compensation. To better express the Doppler phase changes among code element periods ranging code periods in a discrete time system, we express the Doppler phase in the form of a matrix as (17), where each row expresses a period of *TP*, and its range is from 0 to *Q* · *TP*.

$$\begin{bmatrix} ph\_{\mathrm{B1},rj}^{i} \\ 0 + 0, & \dots, & 0 + nT\_{\mathrm{C}}, & \dots, & 0 + (\mathrm{N\_{\mathrm{c}}} - 1)T\_{\mathrm{c}} \\ \vdots & & \vdots & & \vdots \\ qT\_{\mathrm{P}} + 0, & \dots, & qT\_{\mathrm{P}} + nT\_{\mathrm{c}}, & \dots, & qT\_{\mathrm{P}} + (\mathrm{N\_{\mathrm{c}}} - 1)T\_{\mathrm{c}} \\ \vdots & & \vdots & & \vdots \\ (Q-1)T\_{\mathrm{P}} + 0, & \dots, & (Q-1)T\_{\mathrm{P}} + nT\_{\mathrm{C}} & \dots, & (Q-1)T\_{\mathrm{P}} + (\mathrm{N\_{\mathrm{c}}} - 1)T\_{\mathrm{c}} \end{bmatrix} \tag{17}$$

where we can find that the number column is incremented by *Tc* in each row and is incremented by *TP* in each column. Therefore, we can take the FFT in the columns of the matrix to first acquire the Doppler information of targets. Next, we compensate the matrix in each *TP* row for the duration based on the compensation phase. Finally, we can obtain the compensated DDM results by the proposed MTDC method.

#### *3.2. Resolution Study and Peaks Identifying Methods*

In this part, we study the resolution performance of the BDS signal. As shown in Table 1 in the simulation section, the B3I signal is transmitted with a bandwidth of 20.46 MHz. Figure 3 shows the autocorrelation curves of the B3I signal and two chirp signals with bandwidths of 10.23 MHz and 20.46 MHz.



**Figure 3.** Autocorrelation curves of GNSS B3I and chirp signals.

In Figure 3, we find that the main lobe of the B3I signal is wider than that of the chirp signal at 20.46 MHz. According to Formula (18), where C is the electromagnetic wave velocity and *B* is the signal bandwidth, the range resolution of the two chirp signals can be calculated as 29.32 m. The resolution of the B3I signal is approximately 30 m, as shown in Figure 3.

$$
\Delta r = \mathbb{C} / B \tag{18}
$$

To further study the resolution of the B3I signal, we simulated two-point target detection with interval of 29 m and 33 m by using the B3I signal, setting the chirp signal at a 10.23 MHz bandwidth as a comparison. Figure 4a shows the detection peaks of the GNSS B3I signal and chirp signal reflected by two-point targets with a distance of 29 m. We found that the solid line peaks overlap completely, which means that the B3I signal cannot distinguish the two-point targets. The chirp signal, with a bandwidth of 10.23 MHz, can clearly detect two-point target peaks. Next, in Figure 4b, where the distance between two-point targets is 33 m, two signals show two target peaks, and the sidelobe of the B3I signal in the center is higher, while its sidelobes on both sides are lower than those of the chirp signal.

**Figure 4.** Two-point target detection main lobe curves of GNSS B3I and chirp signals: (**a**) The distance between the two targets is 29 m; (**b**) the distance between two targets is 33 m.

Next, we studied the sidelobe between the two target peaks of the B3I signal detection curves. Figure 5 shows the normalized sidelobe peak curves between two target peaks of the GNSS B3I signal as the distance between two targets. We found that the sidelobe peaks gradually decreased until the distance was more than 30 m, and in the range of 30~60 m, they had a linear decline as the distance increased. The results show that the GNSS B3I signal can realize a clear target peak resolution with an interval of targets near 30 m. After further simulation, we found that 33 m is the suitable interval for B3I signals with clear and smaller peak resolutions.

In BDS-R detection, the echo intervals are smaller than the resolution of the B3I signal, and in the case of bistatic detection, intervals of targets vary with the bistatic angle. To acquire a clear and stable resolution, we propose sampling the echoes with an interval of 33 m, and the RCS data for comparative verification can also be processed with an interval of 33 m.

**Figure 5.** Sidelobe between two main lobe peak ratios of the B3I and chirp signals.

#### *3.3. Resolution Study and Range Inversion Methods*

In this section, we study the peak inversion methods. For the limited resolution problem, we propose a correlation peak inversion method to determine the specific range information of targets based on one-dimensional detection data and the corresponding signal autocorrelation function. Moreover, the inversion result accuracy is related to the sampling frequency of the echoes. We propose two methods to find the approximate ranges of target positions in the time domain. In method 1, without target estimation information, we found the range as follows: (1) Acquiring matched filter results, we filtered it with a normalized amplitude > *a*0. (2) Next, we divided the results into small segments, between each pair of segments. There is a low amplitude period with a length of *b*0, where *a*<sup>0</sup> and *b*<sup>0</sup> can be utilized to adjust the target range for better adaptability of detecting targets. In method 2, we can obtain the rough target position in the echoes by rough detection. The inversion ranges can be expressed as *Range* <sup>=</sup> [*range*1,*range*2,...,*rangem*,...,*rangeM*], and the corresponding sampling points can be expressed as *N* <sup>=</sup> [*n*1, *<sup>n</sup>*2,..., *<sup>n</sup>*m,..., *<sup>n</sup>*M]. The inversion in *rangem* can be expressed as

$$\sum\_{q=1}^{n\_m} \left( s(n - n\_m + q) - s(n + q) \right) Amp(q) = M(range\_m) \tag{19}$$

where *Amp* is the modulating amplitude of the received signals. When clutter, noise, and the correlation sidelobe have little influence, the specific radar cross-section (RCS) distribution can be expressed as

$$RCS(range\_m) = Amp\tag{20}$$

According to the literature [28], clutter and noise influence can be suppressed effectively during long-term accumulation. Because the RCS of the simulation ship is large, we accumulated 10 s to effectively suppress the clutter and noise in simulation 2 and corrected

the distance and velocity migration. Therefore, we can utilize the inversion method to acquire the specific range information of targets.

#### *3.4. Target Information Inversion Method Based on GNSS-R Detection*

In signal processing, the detection results of different satellites can be separated by signal range codes. After acquiring multiple satellite detection results for the targets, the different target peaks in the range dimension can be distinguished when the targets keep a sufficiently large distance. Additionally, we can also distinguish them by the difference in beam direction. To acquire the precise coordinate information of targets in bistatic detection, we can determine the target coordinates by detecting the target range and velocity information, satellite information, and receiver information.

$$N^{i,j} \cdot T\_P \cdot \mathbb{C} + r^{i,j} = \sqrt{\frac{\left(\mathbf{x}\_s^i - \mathbf{x}\_t^j\right)^2 + \left(y\_s^i - y\_t^j\right)^2 + \left(z\_r^i - z\_t^j\right)^2}{+\sqrt{\left(\mathbf{x}\_r - \mathbf{x}\_t^j\right)^2 + \left(y\_r - y\_t^j\right)^2 + \left(z\_r - z\_t^j\right)^2}}\tag{21}$$

$$\begin{cases} (\mathbf{v}\_r \cdot \boldsymbol{\Theta}\_r^{i,1} + \mathbf{v}\_s^i \cdot \boldsymbol{\Theta}\_r^{i,2})/\lambda + (\mathbf{v}\_t^{i,j} \cdot \boldsymbol{\Theta}^{i,j,1} + \mathbf{v}\_s^i \cdot \boldsymbol{\Theta}^{i,j,2})/\lambda\\ = M^{i,j} \cdot PRF + f d^{i,j} \end{cases} \tag{22}$$

The target coordinate information can be acquired by the satellite coordinate information from the signal data and the range detection results, as shown in Formula (21). We can also obtain the target velocities through the Doppler detection results as (22), where *xi <sup>s</sup>*, *y<sup>i</sup> <sup>s</sup>*, *z<sup>i</sup> <sup>s</sup>* are the coordinates of satellite *<sup>i</sup>*, and *<sup>x</sup><sup>j</sup> <sup>t</sup>*, *y j <sup>t</sup>*, *z j <sup>t</sup>* are the coordinates of target *j*. *<sup>N</sup>i*,*<sup>j</sup>* , and *ri*,*<sup>j</sup>* are the cycle number and distance during a cycle of detected results of target *j* and satellite *i,*, respectively. In (22), **v***i*,*<sup>j</sup> <sup>t</sup>* · <sup>θ</sup>*i*,*j*,1 expresses the component of the velocity of target *i* in the direction of the bistatic angular bisector between satellite *i* and the receiver, and **v***i <sup>s</sup>* · <sup>θ</sup>*i*,*j*,2 is the component of the velocity of satellite *<sup>i</sup>* in the direction of target *j.* **<sup>v</sup>***<sup>r</sup>* · <sup>θ</sup>*i*,1 *r* and **v***<sup>i</sup> <sup>s</sup>* · <sup>θ</sup>*i*,2 *<sup>r</sup>* are the corresponding velocity components between target *i* and the receiver. *Mi*,*<sup>j</sup>* and *f di*,*<sup>j</sup>* are the Doppler cycle number and Doppler value of the detected results of target *j* and satellite *i,* respectively.

#### *3.5. Experimental Technical Scheme*

Figure 6 shows the experimental technical scheme of multiple targets and multiple-GNSS-R-signal detection. The technical route of the detection can be divided into three parts: (1) The BDS MEO satellite information can be obtained from the database [29] and the signal can simulate the parameter, as observed in the literature [27,28]. We found the ship model using the physical optics scattering method and simulated the sea clutter with the weibull distribution model; (2) After establishing the detection model, the echo was simulated and preprocessed. Where the echo was sampled, the distance migration and velocity migration were corrected, and the direct part of echo was suppressed; (3) For the third part, we simulated the DDM, 33-meter one-dimensional range image, and inverse image with multiple-target Doppler compensation, 33-meter peak identifying methods, and inverse methods. The proposed method has been bolded in the picture.

**Figure 6.** The technical process of multiple GNSS-R signals and multiple-targets detection.

#### **4. Simulation Results and Discussion**

In this section, two sets of simulation results are given to evaluate the effectiveness of the proposed MTDC method, the target peaks resolution method and the target peaks range inversion method based on the B3I signal of four BDS satellites and the MIMO system. The system parameters and target parameters are listed in Tables 1 and 2. The geometry schematic diagram of the simulation is plotted in Figure 6, where the receiver was set as airborne-based to better reflect the Doppler frequency difference of each target during GNSS-R detection. During the 10 s pulse accumulation time in simulation 2, we set the receiver plane to move horizontally and at a constant speed according to the given direction and speed, and the ship does not change its track. The influence of the Doppler and distance values instigated by the aircraft movement can be compensated for by the known aircraft movement parameters combined with the geometric model. Additionally, we set the linear delay compensation to be once per pulse period to remove the delay effect caused by ship motion.

#### *4.1. Simulation 1*

In simulation 1, we detected three sets of point targets based on the B3I signals of four BDS satellites without considering clutter and noise. The point target information is expressed in Table 2, where the points were set with an interval of 30 m, and its geometric distribution was plotted in Figure 7. In particular, to better study the detection property, the RCS ratios of target 1, target 2, and target 3 were set to 2:5:15.


#### **Table 2.** Target parameter.

**Figure 7.** Geometry schematic diagram of multi-signal and multi-target detection simulation.

To study the target peak resolution effectiveness of the sampling method at 33 m, Figures 8 and 9 show the RCS data distribution and one-dimensional range detection results with sampling at 33 m. In Figure 8, the circle line expresses the original RCS distribution of point target 3 under each satellite, and the star line represents the data with a sample at 33 m. In Figure 9, the curves are the detection peaks of point target 3 under each satellite, and its direction and length vary with the detection angle. Taking Figure 8 as a comparison, we found that the small peaks in the detection results of Figure 9 are in good agreement with the sampling data at 33 m in terms of amplitude and range point. The results can be utilized to estimate the RCS distribution of targets and verify target recognition.

To research the Doppler and range information of detection by four satellites, the simulation shows the DDM results of three sets of point targets of the B3I signal from four satellites separately in Figure 10. Because the resolution of B3I is near 30 m, each set of point targets was detected as a few light points in the DDM, and according to the number of points set in each target, we identified that the highlights are target 3, target 2, and target 1, in order from the strongest to the weakest. Among the four detection results, the Doppler frequency of each target ranged from −1000 to 1900 Hz for the high speed and different radiation angle of the satellites, but the difference in the Doppler value among targets was stable at approximately 400 Hz on account of the target points and receiver location. The motion of the receiver and the target determine the Doppler difference among targets. Moreover, the range and velocity information of the detection of four satellites can

determine the target coordinate information using the target information inversion method mentioned previously.

**Figure 8.** RCS data of point target 3 under the detection of 4 BDS satellites: (**a**) satellite 1; (**b**) satellite 2; (**c**) satellite 3; and (**d**) satellite 4.

**Figure 9.** *Cont*.

**Figure 9.** One-dimensional range detection results of point target 3 under the detection of 4 BDS satellites with 33-meter echo sampling: (**a**) satellite 1; (**b**) satellite 2; (**c**) satellite 3; and (**d**) satellite 4.

**Figure 10.** DDM of point targets under the detection of 4 BDS satellites: (**a**) satellite 1; (**b**) satellite 2; (**c**) satellite 3; and (**d**) satellite 4.

#### *4.2. Simulation 2*

In simulation 2, we designed a simulation to detect three sets of ship targets based on the B3I signals of satellite 3 considering Gaussian noise and sea clutter based on the Weibull model. The ship target information and geometric diagram are also expressed in Table 2 and Figure 7.

In Figures 11 and 12, we show the Doppler compensation effectiveness comparison of the two compensation methods. In Figure 11, we simulated the one-dimensional range detection compensated results of three ship targets without target estimation information, and the results are compensated for with single Doppler channels and the MTDC method. Figure 11a–c show the incomplete peak amplitudes of three ship targets, and Figure 11d shows the target peaks without Doppler phase influence, which shows the effectiveness of MTDC in estimating no-target information situations. Figure 12 shows the compensation results with acquired target estimation information. Figure 12a–c show that each target peak can return to the normal level under the corresponding Doppler compensation and that other target peaks lose part of their amplitudes. Figure 12d shows the MTDC method compensation results expressing all target peaks, obtaining their maximum level. The results of Figures 11 and 12 verify that the proposed MTDC method can remove the Doppler phase influence of multiple target detection.

**Figure 11.** Doppler compensation comparison of method 1 in a one-dimensional range image using B3I signals of satellite 3: (**a**) Results with Doppler channel 1 compensation; (**b**) results with Doppler channel 2 compensation; (**c**) results with Doppler channel 3 compensation; (**d**) results with the proposed MTDC method based on multiple Doppler channels.

**Figure 12.** Doppler compensation comparison of method 2 in a one-dimensional range image by B3I signals of satellite 3: (**a**) results with first target Doppler compensation; (**b**) results with second target Doppler compensation; (**c**) results with third target Doppler compensation; (**d**) results with the proposed MTDC method based on multiple estimated target Dopplers.

Based on the same sampling at 33 m, we plotted the one-dimensional range detection results of three ship targets by satellite 3 in Figure 13, where Figure 13a,c,e show the sampling detection results at 3 m and the sampling detection results at 33 m for ships 1, 2, and 3, respectively. Compared with the RCS data of Figure 13b,d,f, we found that the sampling detection at 3 m can express the general outline of the ship's target RCS without accurate amplitude and range point information. The small peaks in the sampling results at 33 m, which are in good agreement with the sampling data obtained at 33 m, express segmented scattering strong points with accurate amplitude and computable range points. The results can be more accurate in the service of target recognition.

**Figure 13.** One-dimensional range detection results of 3 ship targets under the detection of satellite 3 with 33-meter echo sampling and 3-meter echo sampling and the corresponding RCS data distribution: (**a**) ship 1 detection result; (**b**) ship 1 RCS data; (**c**) ship 2 detection result; (**d**) ship 2 RCS data; (**e**) ship 3 detection result; (**f**) ship 3 RCS data.

To further study the resolution improvement method, we simulated the detection results of satellite 3 with range inverse methods to acquire the more delicate RCS distribution of targets based on the sampling frequency *fs* = 10\**B*. Figure 14 shows the range inverse results of method 1, in which we screened the probable target information by amplitude and segment; the results are given by a certain length in the low-amplitude region. The dotted line in Figure 14 expresses the RCS data, and the circle line, which is the inverse result, which does not consider clutter, is in good agreement with the data in terms of the amplitude and range point, especially at strong scattering points. Moreover, we utilized the small square line to express the inverse results with clutter and noise backgrounds with 10 s pulse accumulation. It shows that target peaks can be recognized from clutter and noise backgrounds with the proposed inverse method.

**Figure 14.** Range inversion results of 3 ship targets under the detection of BDS satellite 3 with range inversion method 1: (**a**) ship 1; (**b**) ship 3; (**c**) left part of ship 2; (**d**) right part of ship 2.

In the other situation, since we have estimated the target information, we can acquire rough target distance ranges and reverse them using the range inverse method 2. Figure 15 shows the range inverse results based on target estimation information, where ship 2 can notably be expressed more completely and accurately than with method 1.

**Figure 15.** Range inversion results of 3 ship targets under the detection of BDS satellite 3 with target estimation information by range inversion method 2: (**a**) ship 1; (**b**) ship 2; (**c**) ship 3.

Starting from the two dimensions of Doppler and range, Figure 16 plots the DDM results of three ship targets under the detection of the B3I signal of satellite 3, where pictures (a) and (b) express whether the results consider clutter and noise. In Figure 16a, combined with target peak identification and target range inverse results, we found that the highlight groups were target 1, target 3, and target 2 from left to right. The ship targets can be detected as several highlights with a range resolution of near 30 m and a Doppler frequency resolution of approximately 8 Hz. Due to the high speed of the receiver and different received angles, the targets' Doppler frequency varies in the range of 400 Hz, and there are 8 Hz offsets in the scattered point in a target. After introducing clutter and noise backgrounds, Figure 16b shows the detection results with clutter distribution, as shown in Figure 7, with 10 s pulse accumulation. We found a clutter highlight group near the 0 Hz Doppler channel with a normalized amplitude of 0.6. In the range dimension, a few clutter peaks overlapped target 1, which protected the detection of target 1 from clutter influence. Some cluttered peaks overlapped target 2 and target 3, but the cluttered peaks are small and dispersed, which reduces the influence of the detection of the two targets. Moreover, the cluttered strong scatter point does not overlap in the simulation, supporting better detection results for the signal process. This situation has a large probability of overlapping since the scattering intensity of sea clutter is random.

**Figure 16.** DDM detection results of 3 ship targets under the detection of satellite 3: (**a**) without clutter and noise influence; (**b**) with clutter background.

Figures 12–16 show the one-dimensional range image Doppler compensation results, one-dimensional range image target peak identification results, range inverse results, and DDM results of three ship targets in simulation 2, which verify that multiple fetched ship targets can be detected with better strength using the MDTC method with noise and sea clutter backgrounds in a 10 s accumulation, and the ship RCS distribution structure can be further acquired by the proposed target peak identification method and range inversion method.

#### **5. Conclusions**

This paper studies multiple-BDS signal reflection detection for multiple moving targets in the sea. It analyzes the influence of multiple-target Doppler phases and proposes the MDTC method to suppress influence under target estimation and not under target estimation. The results show the proposed method can effectively restore target peaks under the influence of four kinds of Doppler frequencies in GNSS-R detection, where the Doppler frequency variation range is −1000–1900 Hz. For further study of target peak detection, simulations of point target and ship target detection were performed to verify the effectiveness of the proposed peak identification method and range inversion methods. The peak identification results of four satellites and three ship targets showed corresponding peaks which matched the RCS data of detected targets. The range inversion results show a better resolution of target peak detection results, which can obtain the specific RCS distribution of ships in the hundred-meter class. At last, noise and clutter backgrounds were added to simulation 2 under 10 s pulse accumulation, which shows some influence on the sidelobe of target 2 and target 3, and the three target main peaks can still be detected relatively clearly. In the future, we will study the detection simulation of small-scale targets in the sea and detection with various satellite detection angles.

**Author Contributions:** Conceptualization, methodology, formal analysis, writing—original draft preparation, X.L.; data curation, L.W. and W.J.; writing—review and editing, X.L., J.L. and M.Z.; supervision, M.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the National Natural Science Foundation of China under Grant No. 62171351, No. 6217011640, No. 62001343, No. 41901267, and No. 41701386 and the China Postdoctoral Science Foundation under Grant No. 2021T140532.

**Data Availability Statement:** The satellite coordinates based on the geodetic coordinate system were obtained from http://www.csno-tarc.cn/datacenter/ephemeris (accessed on 1 January 2021). The signal PRN parameters were obtained from literature [27,28].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

