*Article* **Improving the iGNSS-R Ocean Altimetric Precision Based on the Coherent Integration Time Optimization Model**

**Xuezhi Sun 1,†, Wei Zheng 1,2,\*,†, Fan Wu 1,† and Zongqiang Liu 1,3**


† Those authors contributed equally to this paper.

**Abstract:** Improving the altimetric precision under the requirement of ensuring the along-track resolution is of great significance to the application of iGNSS-R satellite ocean altimetry. The results obtained by using the empirical integration time need to be improved. Optimizing the integration time can suppress the noise interference from different sources to the greatest extent, thereby improving the altimetric precision. The inverse relationship between along-track resolution and signal integration time leads to the latter not being infinite. To obtain the optimal combination of integral parameters, this study first constructs an analytical model whose precision varies with coherent integration time. Second, the model is verified using airborne experimental data. The result shows that the average deviation between the model and the measured precision is about 0.16 m. The two are consistent. Third, we apply the model to obtain the optimal coherent integration time of the airborne experimental scenario. Compared with the empirical coherent integration parameters, the measured precision is improved by about 0.1 m. Fourth, the verified model is extrapolated to different spaceborne scenarios. Then, the optimal coherent integration time and the improvement of measured precision under various conditions are estimated. It was found that the optimal coherent integration time of the spaceborne scene is shorter than that of the airborne scene. Depending on the orbital altitude and the roughness of the sea surface, its value may also vary. Moreover, the model can significantly improve the precision for low signal-to-noise ratios. The coherent integration time optimization model proposed in this paper can enhance the altimetric precision. It would provide theoretical support for the signal optimization processing and sea surface height retrieval of iGNSS-R altimetry satellites with high precision and high along-track resolution in the future.

**Keywords:** coherent integration time optimization model; global navigation satellite systems reflectometry (GNSS-R); ocean altimetry precision; waveform correlation; signal optimization processing

#### **1. Introduction**

Accurately measured sea surface height (SSH) is one of the critical parameters of marine ecosystem monitoring, which is of paramount significance to applications such as fishery, oil drilling, and commercial navigation. Global navigation satellite systems reflectometry (GNSS-R) or the passive reflection and interference system (PARIS) can perform ocean altimetry as a passive remote sensing technique [1]. It has some unique advantages compared with the traditional tide gauge station and monostatic radar ocean altimetry. On the one hand, the receiver can capture multiple global navigation satellite system (GNSS) signals simultaneously, which increases the spatial coverage. On the other hand, as a novel bistatic passive remote sensing method, it has the characteristics of low cost, low power consumption, all-weather and high time revisit rate, which can make up for the shortcomings of the existing technology to a large extent. Since the technology

**Citation:** Sun, X.; Zheng, W.; Wu, F.; Liu, Z. Improving the iGNSS-R Ocean Altimetric Precision Based on the Coherent Integration Time Optimization Model. *Remote Sens.* **2021**, *13*, 4715. https://doi.org/ 10.3390/rs13224715

Academic Editors: Shuanggen Jin and Gino Dardanelli

Received: 19 October 2021 Accepted: 17 November 2021 Published: 22 November 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/).

was proposed, ground-based, air-based, and space-based altimetry experiments have been conducted successively to verify its feasibility [2–4].

Some signal processing techniques for GNSS-R ocean altimetry have been proposed. Conventionally, GNSS-R altimetry involves cross-correlating the reflected signals with the local replicas of the open navigation signals (cGNSS-R) [5]. The encrypted signal (e.g., P(Y)) has a broader bandwidth, which will bring a sharper autocorrelation function (ACF) and a higher range precision. With it, the reconstructed code GNSS-R (r-GNSS-R) altimetry improves the precision by reproducing the encrypted signals [6]. The interferometric altimetry (iGNSS-R) utilizes the full composite components of the direct signal to correlate with the reflected signal [5]. The partial interferometric GNSS-R (piGNSS-R) altimetry extracts a portion of the encrypted signal from the full composite components to obtain a sharper ACF, which is an extension of iGNSS-R [7]. Nevertheless, in the case of iGNSS-R and piGNSS-R, the loss of signal-to-noise ratio (SNR) caused by dual-channel noise is more severe than that of cGNSS-R, which leads to a two-fold improvement in precision only [8]. The SNR should be improved by increasing the signal processing time, as this will diminish the influence of thermal noise and speckle noise on altimetric precision [9]. However, a longer integration time corresponds to a longer specular point (SP) movement distance. Considering the moving speed of the sub-satellite point (km/s) and the size of the spatial footprint (~10 km), the actual inverted single height represents the average over the larger sea area which further reduces the along-track resolution [10]. The previous study has shown that in the spaceborne iGNSS-R altimetry scenario, in order to achieve an altimetric precision better than 20 cm, the signal processing time required for a single height retrieval is approximately 10 s, and the along-track spatial resolution is about 65 km [8]. If the iGNSS-R ocean altimetry with high precision and high along-track spatial resolution can be realized, it would provide critical data and information resources for small and medium scale ocean phenomena monitoring, high temporal and spatial resolution ocean gravity field model establishment, and other earth science research.

One of the data products of GNSS-R observation is the delay-Doppler maps (DDM). Fernando et al. detailed the corresponding relationship between the DDM and its spatial location [11]. The waveform corresponding to zero Doppler is extracted from the calibrated DDM waveform for GNSS-R code-delay altimetry, and the delay information of the reflected signal is extracted from this waveform. Commonly used retrieval algorithms include the peak point of the first derivative (DER), the peak of the waveform (MAX), and the half-power point (HALF) [8,12,13]. Obtaining the inversion delay requires correction of numerous parameters including the receiver and transmitter orbits, ionospheric errors, tropospheric errors, and antenna baseline offsets [14]. For the evaluation of the altimetry performance, Li et al. considered both the altimetric precision and the altimetric accuracy. Altimetric accuracy is primarily affected by systematic errors. The altimetric precision is mainly induced by the randomness of the received signal caused by thermal noise and speckle noise. Theoretically, the higher the SNR of the signal, the more accurate the SP delay information extracted from the waveform, and the better the altimetric precision [15].

It is possible to increase the altimetric precision by optimizing the payload and signal post-processing. Payload optimization includes increasing antenna gain, refining antenna pointing and receiver bandwidth, etc. [16,17]. For signal post-processing, one method to improve the altimetric precision is to increase the signal's coherent integration time and incoherent average number. Both methods inhibit the noise introduced in different ways. The former is to suppress the thermal noise introduced at the receiver end, while the latter is to suppress the speckle noise introduced in the glistening area near the SP. The total integration time of signal processing is the product of the coherent integration time and incoherent average number. In theory, high precision requires an increased integration time. Despite this, the inverse relationship between the along-track resolution and the integration time prevents the latter from increasing indefinitely. It is typically necessary to sacrifice a certain degree of precision when performing altimetry tasks requiring high spatial resolution along the orbit. To achieve an optimal altimetric precision for a given

spatial resolution along track, it is necessary to optimize the combination of coherent and incoherent parameters. Accordingly, Martin-Neira et al. derived an altimetric precision prediction model with coherent integration time as a variable [18], but this model overlooked the influence of the correlation between the waveforms. You et al. constructed a waveform correlation model from the time domain and frequency domain to predict the upper limit of the waveform coherence time [19,20]. Li et al., reconstructed the altimetric precision model by considering the correlation of the waveforms which had reasonable consistency with the measured results [21], but the influence of the coherent integration time on the precision was not discussed in detail. Currently, traditional intermediate frequency (IF) data processing uses empirical parameters to approximate an optimal precision. The coherent integration time normally takes 10 ms in airborne scenarios [21] and 1 ms in spaceborne scenarios [22]. However, the actual optimal coherent integration processing depends on factors such as the position and relative motion of the transmitter and receiver. This needs to be estimated through accurate modelling.

Different from previous studies, in order to improve the precision of iGNSS-R ocean altimetry, this study constructs a coherent integration time optimization model by deducing the conversion relationship among coherent integration time, waveform correlation, and altimetric precision. Furthermore, the model can more accurately estimate the variation of precision with coherent integration time in different iGNSS-R altimetry applications so as to optimize the final precision result.

#### **2. Signal Processing and Height Inversion**

In order to examine the accuracy of the coherent integration optimization model, we use the data obtained by the Institute for Space Studies of Catalonia (IEEC) through an airborne experiment on the Baltic Sea on the 3 December 2015, for verification. The delay difference between direct and reflected signals is calculated using interference processing. The raw IF datasets are processed on the ground with the software receiver, including data acquisition, processing, height retrieval, and precision calculation.

The aircraft's altitude was about 3 km during the experiment, and the velocity was about 50 m/s. The direct and reflected signals are captured by the 8-element phased array antennas of RHCP (right-handed circular polarization) and LHCP (left HCP) respectively. The radiofrequency (RF) signals received by the antenna elements are filtered, amplified, and down-converted into IF signals. The IF signals are quantized by a comparator, and the quantized signals are connected to a field-programmable gate array (FPGA) through D-type flip-flops whose parallel processing capability is used for sampling. Analogue signals received by an element are quantized into one-bit in-phase components and one-bit quadrature components per sample. The sampling frequency is 80 MHz. Finally, the sampled digital signal is transmitted into the ground receiving station through the PCIe bus [23]. The acquired data is processed by a software-defined receiver. A simplified diagram of data processing and altimetry retrieval is depicted in Figure 1.

**Figure 1.** Simplified flowchart of airborne IF data processing.

#### *2.1. Signal Processing*

#### 2.1.1. Data Fetch

The IF data is read in 4 bytes at a time, which correspond to the in-phase and quadrature components of the 16 antenna elements sampled at a time. The binary data is converted from 0, 1 code to 1, −1 code. That means the logic level is converted to a polarity nonreturn-to-zero level. The direct signal and the reflected signal sampled each time can be expressed as [24]:

$$\begin{array}{l} S\_{\text{up}}^{nT\_{\mathcal{E}}}(k) = \sum\_{i=1}^{8} s\_{i,\text{up\\_I}}(k) + j \sum\_{i=1}^{8} s\_{i,\text{up\\_Q}}(k) \\ S\_{\text{down}}^{nT\_{\mathcal{E}}}(k) = \sum\_{i=9}^{16} s\_{i,\text{down\\_I}}(k) + j \sum\_{i=9}^{16} s\_{i,\text{down\\_Q}}(k) \end{array} \tag{1}$$

where *SnTc* up (*k*) and *SnTc* down(*k*) are the amplitudes of the kth sample of the direct signal and the reflected signal during the nth coherent integration, *s* is the signal component of each antenna element, *i* is the antenna element number, up\_I and up\_Q respectively represent the in-phase and quadrature components of the zenith pointing antenna, and down\_I and down\_Q denote the in-phase and quadrature components of the nadir pointing antenna respectively.

#### 2.1.2. Coherent Integration

The coherent integration time is set to *Tc* and the sampling frequency *f*<sup>s</sup> is 80 MHz. In this case, *f*<sup>s</sup> × *Tc* samples constitute the sequence of each coherent integration, which is expressed as *<sup>S</sup>nTc* up [*K*]=[*SnTc* up (1), *<sup>S</sup>nTc* up (2)... *<sup>S</sup>nTc* up (*f*<sup>s</sup> <sup>×</sup> *Tc*)]. Removing single-frequency interference is essential before coherent integration can be performed to prevent it from influencing the measurement result. Through fast Fourier transform (FFT), the time domain signal *SnTc* up [*K*] and *SnTc* down[*K*] are transformed into the frequency domain *<sup>S</sup>nTc* up [ *<sup>f</sup>* ] and *SnTc* down[ *f* ], and the amplitude anomalies in the frequency spectrum are identified and filtered out. The filtered direct signal and the reflected signal are cross-correlated in the frequency domain to obtain the complex waveform *<sup>y</sup>*(*nTc*, *<sup>τ</sup>*) = *<sup>F</sup>*−1(*S*∗*nTc* up [ *<sup>f</sup>* ] <sup>×</sup> *<sup>S</sup>nTc* down[ *f* ]). It should

be declared that the phased array antenna is used in this experiment, and beamforming is usually employed to improve the SNR in the signal post-processing. However, the effect of beamforming may submerge the influence of the coherent integration time on the SNR of power waveforms. To reveal the role of coherent integration time, this research abandons beamforming.

#### 2.1.3. Retracking and Incoherent Average

The complex waveforms after coherent integration need to be incoherently averaged to degrade the influence of speckle noise. With the vertical drift of the aircraft during this period, the waveform series need to be compensated with respect to the first one (retracking) [25]. Figure 2A shows a power waveform obtained after incoherent averaging.

**Figure 2.** Examples of measured and simulated power waveforms. (**A**) Normalized measured power waveform processed from the IF data. The *τ*obs DER is the maximum point of the first derivative. (**B**) Normalized simulated power waveform generated based on the Z-V model. The *τ*<sup>m</sup> DER is the maximum point of the first derivative. The *<sup>τ</sup>*<sup>m</sup> sp is the nominal SP calculated based on the WGS-84 reference ellipsoid.

#### *2.2. Height Inversion*

2.2.1. Delay Estimation and Error Correction

The delay of the reflected signal through the SP can be estimated from the power waveform. The delay estimation can be divided into two types: fixed-point tracking and model fitting. Fixed-point tracking involves tracking a given point of the waveform, such as DER, MAX, and HALF. The DER delay estimation is applied in this study, which can be expressed as [15]:

$$
\tau\_{\rm DEP}^{\rm obs,corr} = \tau\_{\rm DEM}^{\rm obs} + (\tau\_{\rm sp}^{\rm m} - \tau\_{\rm DEM}^{\rm m}) \tag{2}
$$

where *τ*obs DER is the delay corresponding to the maximum of the measured waveform's first derivative (Figure 2A), and *τ*<sup>m</sup> sp and *τ*<sup>m</sup> DER correspond to the nominal SP delay and the maximum of the simulated waveform's first derivative. The nominal SP *τ*<sup>m</sup> sp is the minimum of the reflection path calculated from the WGS-84 reference coordinate when the positions of the transmitter and the receiver are known. The simulated waveform (Figure 2B) is generated according to the Zavorotny and Voronovich (Z-V) model, which takes into account the influence of geometry, instrument configuration, and ocean state [26]. The difference between *τ*<sup>m</sup> sp and *τ*<sup>m</sup> DER obtained from the simulated waveform is used as the deviation correction for DER tracking.

The estimated bistatic delay information of the SP also needs error correction, including tropospheric error *ρ*trop, ionospheric error *ρ*iono, and antenna baseline error *ρ*bl [12,14]. As the ionosphere is located above 60 km, both direct and reflected signals received by the airborne platform go through the same data path for downlink transmission, so the

delay deviation can be ignored. The tropospheric error correction is derived from the Saastamoinen model [21]:

$$
\rho\_{\text{trop}} = \frac{4.6}{\sin e} (1 - e^{-H\_{\text{R}}/H\_{\text{trop}}}) \tag{3}
$$

where *e* is the satellite elevation angle, *H*<sup>R</sup> is the receiver height, and *H*trop is the average tropospheric height. The antenna baseline offset is calculated from the path difference between the zenith pointing antenna and the nadir pointing antenna relative to the SP under the known airborne position and attitude information. The corrected SP delay can be converted to the actual height from the receiver to the sea surface [15].

$$H\_{\varepsilon}^{\text{obs}} = \frac{r\_{\text{DER}}^{\text{obs,corr}} - (p\_{\text{iono}} + p\_{\text{trop}} + p\_{\text{bl}})}{2 \times \sin(\varepsilon)} \tag{4}$$

2.2.2. Height Retrieval and Precision Calculation

After delay estimation and error correction, the vertical distance from the sea surface to the reference ellipsoid can be computed as:

$$\text{SSH} = H\_{\text{c}}^{\text{obs}} - H\_{\text{c}}^{\text{m}} = H\_{\text{c}}^{\text{obs}} - \frac{\mathbf{r}\_{\text{sp}}^{\text{m}}}{2 \times \sin(\varepsilon)} \tag{5}$$

where *H*<sup>m</sup> *<sup>e</sup>* is the height of the aircraft on the WGS84 reference ellipsoid obtained by using bistatic geometry. A fitted piecewise linear function is subtracted from the measured SSH sequence to generate zero mean, near-white noise residuals. The precision is obtained by calculating the standard deviation of the SSH residuals for each trajectory.

$$\sigma\_{\rm h}^{\rm obs}(T\_{\rm c}) = \sqrt{\left\langle \left| \text{SSH}\_{\rm residual}^{\rm k} - \left\langle \text{SSH}\_{\rm residual}^{\rm k} \right\rangle \right|^2 \right\rangle} \tag{6}$$

It is primarily due to the following considerations that linear fitting is used instead of the geoid model. The altimetry precision is affected by zero-mean random error, which is mainly due to the random nature of the received signals caused by thermal noise and speckle noise. The SSH residual after subtracting the linear fit can be used to evaluate randomness. In contrast, precision will be affected by the errors in the geoid model [27].

#### **3. Construction of Coherent Integration Time Optimization Model**

The altimetric precision is related to the power uncertainty of the incoherent average waveform at the SP. The uncertainty value of the incoherent average waveform will be affected by the correlation between the waveforms. To more accurately reflect the variation of precision with the coherent integration time, the waveform correlation is considered in the precision prediction model. Since the correlation between the waveforms varies with the coherent integration parameters, variable conversion is performed on the reconstructed precision model to obtain the coherent integration time optimization model whose precision varies with the coherent integration time in this section.

#### *3.1. The Reconstruction of Altimetric Precision Model*

According to [18], the estimated altimetric precision can be converted by the uncertainty of the incoherent average power waveform at SP.

$$\sigma\_h^{\rm m}(\tau) = \frac{c}{2 \times \cos(i)} \frac{\sigma\_Z(\tau)}{\overline{Z}'(\tau)} = \frac{1}{2 \times \cos(i)} \frac{1}{S\_h(\tau)} \frac{\sigma\_Z(\tau)}{\overline{Z}(\tau)}\tag{7}$$

where *i* is the incident angle, *c* is the speed of light in vacuum, *Z*(*τ*) = *Z*(*t*, *τ*) is the average of the power waveform, and the derivative of which corresponds to *Z* (*τ*). *Sh* = *Z* (*τ*)/*cZ*(*τ*) indicates the logarithmic derivative of the power waveform at SP which

is defined as the altimetric sensitivity. *σZ*(*τ*)/*Z*(*τ*) is the ratio of the standard deviation to power amplitude which is defined as the effective independent incoherent average number [21]. Since *σZ*(*τ*) depends on the correlation between the waveforms, the ratio is derived analytically in this section.

As introduced in Section 2, the complex waveform after coherent integration of direct and reflected signals can be represented by a discrete array *y*(*nTc*, *τ*):

$$y(nT\_{c\prime}, \mathbf{r}) = \left[ y(nT\_{c\prime}, \mathbf{r}\_1), y(nT\_{c\prime}, \mathbf{r}\_2) \dots \dots y(nT\_{c\prime}, \mathbf{r}\_n) \right] \tag{8}$$

where *nTc* is the nth coherent integration process, *Tc* is the coherent integration time which is usually in milliseconds, and *τ<sup>n</sup>* is the code delay whose resolution is inversely proportional to the signal sampling rate. The square of each complex waveform results in a one-shot power waveform expressed as [21]:

$$z(nT\_{c\prime}, \mathbf{r}) = y(nT\_{c\prime}, \mathbf{r})y^\*(nT\_{c\prime}, \mathbf{r})\tag{9}$$

In the case of an incoherent average of *N*<sup>I</sup> one-shot power waveforms, the power waveform can be expressed as:

$$Z(\mathbf{r}) = \frac{1}{N\_{\rm I}} \sum\_{n=0}^{N\_{\rm I}-1} z(nT\_{\rm c}, \mathbf{r}) = \langle z(nT\_{\rm c}, \mathbf{r}) \rangle \tag{10}$$

where  stands for the ensemble average. A power waveform is usually processed in seconds, which can be expressed as:

$$T\_{\rm I} = T\_{\rm c} \times N\_{\rm I} \tag{11}$$

Due to the relative motion of the transmitter and the receiver, the correlation between two complex waveforms separated by *nT <sup>c</sup>* in time is determined by the coherence of the signals from two ocean contribution areas. These regions are elliptical, and the length of the semi-axis is related to *τn*. The correlation between complex waveforms is defined as [19]:

$$\mathcal{C}\_{\mathcal{Y}}(\bar{n}T\_{\mathfrak{c}},\boldsymbol{\tau}) = \left\langle \boldsymbol{y}(nT\_{\mathfrak{c}},\boldsymbol{\tau})\boldsymbol{y}^\*(nT\_{\mathfrak{c}}+\bar{n}T\_{\mathfrak{c}},\boldsymbol{\tau})\right\rangle \tag{12}$$

where *nT <sup>c</sup>* is the time interval between complex waveforms. The spatial geometry of correlation function is illustrated in Figure 3. It can be seen that the correlation between the two sea surface scattering signals gradually decreases as *nT <sup>c</sup>* increases. An analogous correlation between one-shot power waveforms can also be derived as:

$$\mathbb{C}\_{z}(\overline{n}T\_{c\prime}\tau) = \langle \left[ z(nT\_{c\prime}\tau) - \langle z(nT\_{c\prime}\tau) \rangle \right] \left[ z(nT\_{c\prime} + \overline{n}T\_{c\prime}\tau) - \langle z(nT\_{c\prime} + \overline{n}T\_{c\prime}\tau) \rangle \right] \rangle \tag{13}$$

**Figure 3.** Sea surface contribution area corresponding to *<sup>y</sup>*(*nTc*, *<sup>τ</sup>*)*y*∗(*nTc* <sup>+</sup> *nT <sup>c</sup>*, *<sup>τ</sup>*) .

The correlation time of the signal is usually measured in milliseconds, while the time of incoherent averaging can reach the order of seconds. Thus, it can be assumed that there is no correlation between the incoherent average power waveforms, i.e., *CZ*(*t*, *τ*) = 0*t* = 0. The variance of the power waveform can be expressed as [21]:

$$
\sigma\_Z^2(\mathbf{r}) = \mathbb{C}\_Z(0, \mathbf{r}) = \left\langle \left[ Z(\mathbf{t}, \mathbf{r}) - \langle Z(\mathbf{t}, \mathbf{r}) \rangle \right]^2 \right\rangle \tag{14}
$$

The effective incoherent average number can be reconstructed following Appendix A as [21]:

$$\frac{1}{N\_{\text{eff}}} = \frac{\sigma\_Z^2(\tau)}{\overline{Z}^2(\tau)} = \frac{\frac{1}{N\_{\text{I}}} \sum\_{\overline{n}=-(N\_{\text{I}}-1)}^{N\_{\text{I}}-1} \left| \mathbb{C}\_y(\overline{n}T\_{c\text{\text{\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textdegree\textmathbf$$

#### *3.2. The Relationship between Model Parameters and Coherent Integration Time*

As discussed in Section 3.1, the altimetric precision is determined mainly by two parameters, *Sh* and *N*eff. In this section, the variables of these two parameters are derived to be expressed in terms of coherent integration time. The re-derived parameters are brought into Equation (7) to construct the coherent integration time optimization model.

#### 3.2.1. Altimetric Sensitivity

As discussed in [28], the complex waveform is composed of useful signal and noise components

$$y(nT\_{\varepsilon}, \tau) = y\_{\text{s}}(nT\_{\varepsilon}, \tau) + y\_{\text{nd}}(nT\_{\varepsilon}, \tau) + y\_{\text{nr}}(nT\_{\varepsilon}, \tau) + y\_{\text{ndr}}(nT\_{\varepsilon}, \tau) \tag{16}$$

where *y*s(*t*, *τ*) is the cross-correlation value for the useful signal term, and *y*nd(*t*, *τ*),*y*nr(*t*, *τ*) and *y*ndr(*t*, *τ*) are the cross-correlation values of the direct and reflected noise terms.

Assuming that the signal component has no correlation with the noise components, the power waveform can be expressed as [29]:

$$Z(\tau) = \left\langle \left| y\_s(nT\_{\varepsilon\prime}\tau) \right|^2 \right\rangle + \left\langle \left| y\_{\rm nr}(nT\_{\varepsilon\prime}\tau) \right|^2 \right\rangle + \left\langle \left| y\_{\rm nd}(nT\_{\varepsilon\prime}\tau) \right|^2 \right\rangle + \left\langle \left| y\_{\rm ndr}(nT\_{\varepsilon\prime}\tau) \right|^2 \right\rangle \tag{17}$$

The components above can be interpreted as the expression for the coherent integration time [18,20]:

$$\begin{cases} \left|y\_{s}(nT\_{\rm c},\tau)\right|^{2}\right> = 2P\_{\rm d} \iint \frac{2\eta\_{\rm c}G\_{\rm r}(\overset{\rightarrow}{\rho})\rho\_{\rm d}}{4\pi R\_{\rm c}^{2}(\overset{\rightarrow}{\rho})R\_{\rm c}^{2}(\overset{\rightarrow}{\rho})} \times \Lambda^{2}(\Delta\tau) \times \text{sinc}^{2}(\Delta fT\_{\rm c})d^{2}p\\ \left|y\_{\rm nr}(nT\_{\rm c},\tau)\right|^{2}\right> = \frac{2kT\_{\rm nc,r}}{T\_{\rm c}}2P\_{\rm d}\\ \left|y\_{\rm nd}(nT\_{\rm c},\tau)\right|^{2}\right> = \frac{2kT\_{\rm nc,d}}{T\_{\rm c}}2P\_{\rm r}\\ \left|y\_{\rm ndr}(nT\_{\rm c},\tau)\right|^{2}\right> = \frac{2kT\_{\rm nc,r}}{T\_{\rm c}}2kT\_{\rm c\rm c\,d}B \end{cases} \tag{18}$$

where *B* represents the equivalent noise bandwidth of the receiver, *T*rec\_d and *T*rec\_r represent the equivalent input noise temperature of the up-looking and down-looking chains, *P*r represents the total power of the reflected signal at the input of the correlator, *P*<sup>d</sup> represents the total power of the direct signal at the input of the correlator, *P*t represents the power of the transmitted signal, and *G*r represents the antenna gain [30,31]. Based on Equations (17) and (18), the altimetric sensitivity can be represented by the coherent integration time [5].

#### 3.2.2. Effective Incoherent Average Number

As discussed in Equation (15), the effective incoherent average at the SP is closely related to the waveform correlation [32]. Therefore, when the processing time of the incoherent average power waveform is *T*I, the variation of *N*eff with the coherent integration time can be approximated following Appendix B as:

$$\frac{1}{N\_{\text{eff}}} = \frac{T\_{\text{c}} \sum\_{\overline{n}=-\left(\overline{\lambda}\_{\text{f}}-1\right)}^{\overline{\lambda}\_{\text{f}}-1} \left(2P\_{\text{d}} \iint \frac{2R\_{\text{f}}G\_{\text{f}}(\overline{\overline{p}})\overline{\eta}}{4\pi R\_{\text{f}}^{2}(\overline{p})R\_{\text{f}}^{2}(\overline{p})}\Lambda^{2}(\Delta\tau)\text{sinc}^{2}(\Delta f T\_{\text{c}})\exp^{-j2\pi\Delta f(\overline{n}T\_{\text{c}})}d^{2}p + \frac{2P\_{\text{d}}\mathfrak{u}\_{0}^{\text{d}}\,\text{tr}\,(\overline{\mu}T\_{\text{c}})}{\overline{\lambda}\_{\text{c}}}\Delta\overline{\eta} + \frac{R\_{\text{f}}\mathfrak{u}\_{0}^{\text{d}}\,\text{tr}\,(\overline{\mu}T\_{\text{c}})}{4\pi R\_{\text{f}}^{2}(\overline{p})R\_{\text{f}}^{2}(\overline{p})}\Lambda^{2}(\Delta\tau)\text{sinc}^{2}(\Delta f T\_{\text{c}})d^{2}p + \frac{2P\_{\text{d}}\mathfrak{u}\_{0}^{\text{d}}}{\overline{\lambda}\_{\text{c}}} + \frac{R\mathfrak{u}\_{0}^{\text{d}}\mathfrak{u}\_{0}^{\text{d}}}{\overline{\lambda}\_{\text{c}}}\right)}\tag{19}$$

where *n*<sup>d</sup> <sup>0</sup> = <sup>2</sup>*kT*<sup>d</sup> and *<sup>n</sup>*<sup>u</sup> <sup>0</sup> = 2*kT*<sup>u</sup> are the power spectral density of thermal noise from the up-looking and down-looking chains respectively, tri() is the triangular function, Δ*τ* is the difference between *τ<sup>n</sup>* and the scattering point delay *τ*( → *p* ), and Δ*f* is the difference between 0 Doppler and the scattering point Doppler *f*( → *p* ).

Accordingly, the coherent integration time optimization model can be derived by simultaneously applying Equations (7) and (17)–(19). In accordance with Equation (17), a larger coherent integration time is needed to improve the altimetric sensitivity and precision. Nevertheless, the limitation of *N*eff on the altimetric precision leads to a relatively smaller expected value of the coherent integration time. As a consequence, the best precision result is achieved by balancing the time required for the coherent integration of both aspects.

#### **4. Results and Application**

#### *4.1. Validation of Coherent Integration Time Optimization Model*

The proposed model precision results are compared with the experimental results to verify the effectiveness of the coherent integration time optimization model in this section. The measured precision is calculated as described in Section 2. With the power waveform processing time *T*<sup>I</sup> = 1 s unchanged, the precision variation *σ*obs *<sup>h</sup>* (*Tc*) with coherent integration time can be determined by varying the value of *Tc*. To avoid the error interference caused by aircraft turning, two straight flight trajectories are selected for precision calculations. Figure 4B illustrates the inverted SSH sequence when *Tc* is 10 ms during the experiment. The corresponding SSH residual sequence is shown in Figure 4C.

**Figure 4.** Experimental area and IF data processing results. (**A**) The specular points tracks of the raw IF data. (**B**) Measured SSH relative to the WGS84 reference ellipsoid. (**C**) SSH residual after subtracting the fitted value.

The coherent integration time optimization model has been acquired in the following steps as described in Section 3. Obtain auxiliary data during the experiment, including the position and velocity information of the receiver after dual-frequency orbit determination, the GPS satellite position and velocity using the IGS precision ephemeris interpolation method, wind speed (~7 m/s) and payload parameters. These simulation parameters are listed in Table 1. The elevation angle is computed with the bistatic geometry. The computation of the variations in the effective incoherent average number and altimeric sensitivity with coherent integration time can be obtained by applying Equations (17)–(19). By substituting the intermediate results into Equation (7), the coherent integration time optimization model can be calculated.

**Table 1.** Air-based simulation parameters that are consistent with those for the experiment described in Section 2.


A comparison between the variance of the model and the measured result is shown in Figure 5. To evaluate the effectiveness of the model, the mean deviation of precision ' *σ*<sup>m</sup> *<sup>h</sup>* (*Tc*) − *<sup>σ</sup>*obs *<sup>h</sup>* (*Tc*) ( is used. *σ*<sup>m</sup> *<sup>h</sup>* (*Tc*) is the precision calculated by the model, whereas *σ*obs *<sup>h</sup>* (*Tc*) is the measured precision. The results show that the average deviation between the red and blue curves is 0.85 m, and the average deviation between the red and black curves is 0.16 m. Therefore, the proposed model considering the correlation of the waveforms is in good agreement with the measured result.

**Figure 5.** The variation of the precision with coherent integration time under different conditions. Red indicates the measured result. Black indicates the model result considering the correlation between the waveforms. Blue indicates the model result without considering the correlation between the waveforms.

The model curve and the measured curve differ by a certain amount. On the one hand, the measured precision result is generally worse than that of the model, with an average of 0.16 m. This may be due to the lack of beamforming in the signal processing, which causes the SNR to decrease. On the other hand, the trend of the measured curve fluctuates. As opposed to the simulated results, the actual measurement results may contain other stochastic factors such as aircraft mechanical vibrations, attitude changes, and antenna pointing, etc. These random factors can be used to optimize the model by obtaining recorded data of aircraft attitude and vibration in the future.

#### 4.1.1. Altimetric Sensitivity

The reciprocal of the altimetric sensitivity simulation is compared with the measured result, as shown in Figure 6. It can be appreciated that the measured and the simulated curves are in good agreement. As the coherent integration time increases, the reciprocal of the sensitivity approaches a limit value. The difference between the curves may be caused by unaccounted noises [33,34].

**Figure 6.** A comparison of the reciprocal value of altimetric sensitivity between measured (red) and simulated data (black).

#### 4.1.2. Effective Incoherent Average Number

It can be seen from Equation (19) that the effective incoherent average is a nonlinear function of the coherent integration time. The variations of *N*<sup>I</sup> and *N*eff with the coherent integration time are compared, as shown in Figure 7. It reveals that when *T*<sup>I</sup> is constant, *N*<sup>I</sup> and *Tc* are inversely proportional. However, *N*eff and *Tc* are non-linearly related. The difference between *N*eff and *N*<sup>I</sup> will gradually decrease as the coherent integration time increases. If there were no correlation between the waveforms i.e., *Cy*(*nT <sup>c</sup>*, *<sup>τ</sup>*) = 0, *<sup>n</sup>* <sup>=</sup> 0, then *N*eff would be equal to *N*I. *N*eff also varies with the length delay *cτ<sup>n</sup>* as shown in Figure 8. It is a hybrid effect of SNR and Doppler bandwidth at different delays [35].

**Figure 7.** Comparison of the effective incoherent average number (black) and the incoherent average number (red).

**Figure 8.** Effective incoherent average number relative to SP.

#### *4.2. Application of Coherent Integration Time Optimization Model*

The conformity between the model and the measured results indicates that the proposed model can better reflect the influence of coherent integration time on the altimetric precision. Consequently, the model can be used to determine the optimal coherent integration time to improve experimental data processing. In addition, the model can be applied to a variety of altimetry situations to provide a theoretical reference for data optimization processing.

#### 4.2.1. Model Application: Airborne Experiment Scenario

In the airborne altimetry scene described in Section 4.1, the traditional integration parameters cannot yield the optimal precision. In constructing a coherent integration time optimization model based on the same experimental parameters, the optimal coherent integration time is calculated to be 7.5 ms, as shown in Figure 5. By applying this integral parameter to the processing of the measured data, the precision is 0.81 m. As compared to the empirical coherent integration (10 ms) [21], the altimetric precision is improved by 0.09 m.

#### 4.2.2. Model Application: Extrapolation to Spaceborne Scenario

Given the reasonable levels of agreement between the simulated model and the airborne experimental data, the proposed coherent integration time optimization model has been implemented to simulate spaceborne iGNSS-R data. Since there is no dedicated iGNSS-R altimetry satellite at present, the results of model optimization can provide a theoretical reference for improving the in-orbit performance of an iGNSS-R-based ocean altimeter in the future. The altimetric precision is mainly determined by the system instrument parameters and observation geometry. These factors will lead to different model optimization results. According to the model's expression, the main parameters that affect the precision are respectively analyzed, including the receiver's orbital altitude, sea surface roughness, elevation angle, and antenna gain. Previously, these parameters were discussed, and this section re-evaluates them in light of the novel model [31,36]. It needs to be clarified that the Z-V model does not account for the influence of coherent scattering under low wind speed conditions; the minimum wind speed is set to 4 m/s to avoid errors caused by the simulation waveform. The fixed system parameters are summarized in Table 2. The dependence of the precision variation and the optimal coherent integration time on different parameters is analyzed as follows.


**Table 2.** The common simulation parameters for spaceborne altimetry scenarios.

As can be seen from Figure 9:

**Figure 9.** The simulation variation curves of the estimated altimetric precision with the coherent integration time under different conditions. (**A**) When the wind speed is 6 m/s, the elevation is 60◦, the antenna gain is 15 dBi, and the orbit altitude varies. (**B**) When the orbit altitude is 600 km, the elevation is 60◦, the antenna gain is 15 dBi, and the wind speed varies. (**C**) When the orbit altitude is 600 km, the wind speed is 6 m/s, the antenna gain is 15 dBi, and the elevation varies. (**D**) When the orbit altitude is 600 km, the wind speed is 6 m/s, the elevation is 60◦, and the antenna gain varies.

(1) As shown in Figure 9A, the precision decreases with the altitude of the orbit. When the orbit altitude increases from 200 km to 1000 km, the optimal precision decreases from 0.43 m to 0.94 m. This is because, according to Equation (17), the signal power decreases squarely with the transmission distance, causing the SNR to decrease. On the other hand, the receiver velocity decreases with altitude, resulting in increased waveform correlation and decreased effective incoherent average.

As the orbital altitude increases, the optimal coherent integration time increases from 1.5 ms to 3 ms. Higher altitude results in greater energy loss, which calls for a longer time for coherent integration, which improves the SNR and therefore the altimetric sensitivity. The higher the orbital altitude of the receiver, the more noticeable the improvement in precision of the optimal coherent integration time when compared to the empirical spaceborne coherent integration parameters (1 ms). Therefore, the coherent integration time optimization model is preferred for improving the altimetric precision under the low SNR conditions.

(2) As shown in Figure 9B, the precision decreases with wind speed. With an increase in wind speed from 4 m/s to 10 m/s, the optimal precision decreases from 0.57 m to 0.85 m. At low wind speeds, the sea surface is smooth, the waveform coherence is strong, and the effective incoherent average is small. However, the SNR is high, which enhances precision.

The optimal coherent integration time will increase from 1.6 ms to 3 ms with wind speed. This is because high wind speed increases ocean scattering. Simultaneously, the higher the wind speed, the worse the SNR, and the greater the precision improved by the optimal coherent integration time.

(3) In Figure 9C, the optimal precision increases from 0.94 m to 0.63 m as the elevation angle changes from 40◦ to 80◦. On the one hand, increasing the elevation angle shortens the propagation path of the navigation signal, which reduces the signal energy loss. On the other hand, according to Figure 3, the iso-delay area decreases with the elevation angle, thereby increasing the effective incoherent average.

Since both the altimetric sensitivity and the effective incoherent average decrease with the elevation angle, the balance between the two leads to little change in the optimal coherent integration time. Accordingly, as the elevation angle increases, the SNR of the reflected signal upgrades and the improvement of the precision by the optimal coherent integration time becomes insignificant.

(4) As illustrated in Figure 9D, when the antenna gain is increased from 15 dBi to 30 dBi, the optimal precision improves from 0.72 m to 0.31 m. The high gain antenna can increase signal energy and upgrade the SNR. When the antenna gain is increased from 25 dBi to 30 dBi, the optimal precision is only improved from 0.34 m to 0.31 m. This indicates that an increase in antenna gain will not be sufficient to improve precision at high SNRs.

Similarly, the optimal coherent integration time does not greatly improve the precision in the case of high gain. Since the change of the antenna gain does not affect the effective incoherent average, the optimal coherent integration time does not change much with the antenna gain.

(5) In general, the optimal coherent integration time calculated in the spaceborne simulation is shorter than that for the airborne scenario. This is because the satellite receiver has a high velocity, making the Doppler bandwidth and the effective incoherent average large. The shorter coherent integration time can effectively reduce the influence of speckle noise.

The simulation analysis of the above parameters indicates that the coherent integration time optimization model improves the precision more significantly when the SNR of the reflected waveform is relatively weak. Limited by the code bandwidth, improving the SNR of the signal does not increase the precision infinitely but approaches the limit value. Additionally, the optimal coherent integration time varies with the simulation parameters, and the spaceborne scenario tends to be shorter than the airborne scenario. Increasing the altitude of the orbit and the wind speed will also increase the optimal coherent integration time.

#### **5. Conclusions**

This paper constructs a coherent integration time optimization model from the perspective of signal processing to improve iGNSS-R ocean altimetric precision. The research mainly involves three aspects: the processing of IF data and the extraction of precision information; the derivation and verification of the coherent integration time optimization model; and the application of the validated model to airborne and spaceborne altimetry mission scenarios to predict the optimal solution for precision.

(1) To estimate the SP delay for airborne data processing, we use the DER method. The altimetric precision is evaluated based on the random characteristics of the SSH. (2) For the purpose of optimizing the precision performance, we consider the influence of the correlation between waveforms on the covariance of the power waveform, and derive the coherent integration time optimization model from the statistical properties of the waveform. The model is verified by airborne measurements. Results indicate that the average deviation between model and measurement precisions is 0.16 m, which is acceptable. (3) Based on the optimal coherent integration time of the model solution, we processed the experimental data, and the precision is improved by about 0.1 m relative to that obtained using the empirical parameters. The positive validation of the coherent integration time optimization model provides a tool for evaluating the optimal coherent integration time and its precision under different orbital altitudes, sea surface roughness, elevation angles, and antenna gains in spaceborne iGNSS-R altimetry scenarios. The results show that the optimal coherent integration time for the spaceborne scene is shorter than that for the airborne scene. In the case of low SNR, the model can improve the precision to the decimeter level. Due to the combined effect of the SNR and the Doppler bandwidth, the optimal coherent integration time increases with the orbit altitude and the sea surface roughness.

The coherent integration time optimization model developed in this paper is capable of optimizing the waveform processing parameters and improving the altimetric precision. It can be applied to the optimization of data processing and high-precision retrieval for future spaceborne iGNSS-R altimetry missions. In future research, we will also analyze and evaluate the influence of the calibration errors of the instrument, the receiver orbit, and the atmosphere on the precision.

**Author Contributions:** Writing—original draft, X.S.; Writing—review and editing, W.Z., F.W. and Z.L.; X.S., W.Z. and F.W. contributed equally to this paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China under Grant (41774014, 41574014), the Liaoning Revitalization Talents Program under Grant (XLYC2002082), the Frontier Science and Technology Innovation Project and the Innovation Workstation Project of Science and Technology Commission of the Central Military Commission under Grant (085015), the Independent Research and Development Start-up Fund of Qian Xuesen Laboratory of Space Technology (Y-KC-WY-99-ZY-000-025), and the Outstanding Youth Fund of the China Academy of Space Technology.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We would like to thank the Institute of Space Sciences (ICE, CSIC) and the Institute for Space Studies of Catalonia (IEEC) for providing the raw data processing results of the airborne experiment. We would also like to thank IGS for providing SP3 precise orbit documents.

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

#### **Appendix A**

Simultaneous Equations (8) and (13), *Cz*(*nT <sup>c</sup>*, *<sup>τ</sup>*) can be derived as [21]:

$$\begin{array}{lcl} \mathbb{C}\_{\varepsilon}(\overline{n}\Gamma\_{\varepsilon},\mathbf{r}) &= & \langle [z(nT\_{\varepsilon\tau}) - \langle z(nT\_{\varepsilon\tau}),\tau\rangle] \rangle [z(nT\_{\varepsilon\tau} + \overline{n}T\_{\varepsilon\tau},\tau) - \langle z(nT\_{\varepsilon\tau} + \overline{n}T\_{\varepsilon\tau},\tau)\rangle] \rangle \\ &= & \langle [y(nT\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau},\tau) - \langle y(nT\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau},\tau)\rangle] \rangle \\ & \times [y(nT + \overline{n}T\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau} + \overline{n}T\_{\varepsilon\tau},\tau) - \langle y(nT\_{\varepsilon\tau} + \overline{n}T\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau} + \overline{n}T\_{\varepsilon\tau},\tau)\rangle] \rangle \\ &= & \langle y(nT\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau},\tau)y(nT + \overline{n}T\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau} + \overline{n}T\_{\varepsilon\tau},\tau) \rangle \\ & - \langle y(nT\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau},\tau) \rangle \langle y(nT + \overline{n}T\_{\varepsilon\tau},\tau)y^\*(nT\_{\varepsilon\tau} + \overline{n}T\_{\varepsilon\tau},\tau) \rangle \end{array} \tag{A1}$$

In general, the complex waveforms follow the circular complex Gaussian statistics. The correlation function of the one-shot power waveforms can be simplified by using the complex Gaussian moment theorem [37]:

$$\begin{array}{l} \langle y(nT\_{\mathsf{c}},\mathsf{r})y^\*(nT\_{\mathsf{c}},\mathsf{r})y(nT+\tilde{n}T\_{\mathsf{c}},\mathsf{r})y^\*(nT\_{\mathsf{c}}+\tilde{n}T\_{\mathsf{c}},\mathsf{r})\rangle\\ =\langle y(nT\_{\mathsf{c}},\mathsf{r})y^\*(nT\_{\mathsf{c}},\mathsf{r})\rangle\langle y(nT+\tilde{n}T\_{\mathsf{c}},\mathsf{r})y^\*(nT\_{\mathsf{c}}+\tilde{n}T\_{\mathsf{c}},\mathsf{r})\rangle\\ +\langle y(nT\_{\mathsf{c}},\mathsf{r})y^\*(nT\_{\mathsf{c}}+\tilde{n}T\_{\mathsf{c}},\mathsf{r})\rangle\langle y^\*(nT\_{\mathsf{c}},\mathsf{r})y(nT+\tilde{n}T\_{\mathsf{c}},\mathsf{r})\rangle\end{array} \tag{A2}$$

Substituting Equation (A2) into Equation (A1), we can get:

$$\begin{array}{rcl} \mathbb{C}z(\tilde{n}T\_{\mathfrak{c}},\mathsf{r}) &=& \langle y(nT\_{\mathfrak{c}},\mathsf{r})y^\*(nT\_{\mathfrak{c}}+\tilde{n}T\_{\mathfrak{c}},\mathsf{r})\rangle \langle y^\*(nT\_{\mathfrak{c}},\mathsf{r})y(nT+\tilde{n}T\_{\mathfrak{c}},\mathsf{r})\rangle\\ &=& \left|\mathbb{C}y(\tilde{n}T\_{\mathfrak{c}},\mathsf{r})\right|^2 \end{array} \tag{A3}$$

$$\begin{array}{rcl} \mathsf{C}\_{\mathsf{z}}(0,\mathsf{r}) & = & \left| \mathsf{C}\_{\mathsf{y}}(0,\mathsf{r}) \right|^{2} \\ & = & \left\langle y(nT\_{\mathsf{c}},\mathsf{r}) y^{\*}(nT\_{\mathsf{c}},\mathsf{r}) \right\rangle \left\langle y^{\*}(nT\_{\mathsf{c}},\mathsf{r}) y(nT\_{\mathsf{c}},\mathsf{r}) \right\rangle \\ & = & \left\langle y(nT\_{\mathsf{c}},\mathsf{r}) y^{\*}(nT\_{\mathsf{c}},\mathsf{r}) \right\rangle \left\langle y^{\*}(nT\_{\mathsf{c}},\mathsf{r}) y(nT\_{\mathsf{c}},\mathsf{r}) \right\rangle \\ & = & \left\langle z(nT\_{\mathsf{c}},\mathsf{r}) \right\rangle \left\langle z(nT\_{\mathsf{c}},\mathsf{r}) \right\rangle \end{array} \tag{A4}$$

According to Equations (10) and (14), the variance of the power waveforms after incoherent average can be expressed as [21]:

$$\begin{array}{rcl} \sigma\_Z^2(\tau) &=& \mathcal{C}\_Z(0,\tau) \\ &=& \langle [Z(t,\tau) - \langle Z(t,\tau) \rangle][Z(t,\tau) - \langle Z(t,\tau) \rangle] \rangle \\ &=& \langle Z(t,\tau)Z(t,\tau) \rangle - \langle Z(t,\tau) \rangle \langle Z(t,\tau) \rangle \\ &=& \langle Z(t,\tau)Z(t,\tau) \rangle - \langle z(t,\tau) \rangle \langle z(t,\tau) \rangle \end{array} \tag{A5}$$

Substitute Equation (10) into the first term of Equation (A5), we can get:

$$\begin{aligned} \langle Z(\mathbf{t}, \mathbf{r}) Z(\mathbf{t}, \mathbf{r}) \rangle &= \begin{array}{c} \frac{1}{N\_{\mathrm{I}}^{2}} \left\langle \sum\_{i=0}^{N\_{\mathrm{I}}-1} z(iT\_{\mathbf{c}}, \mathbf{r}) \sum\_{j=0}^{N\_{\mathrm{I}}-1} z(jT\_{\mathbf{c}}, \mathbf{r}) \right\rangle \\ &= \frac{1}{N\_{\mathrm{I}}^{2}} \sum\_{i=0}^{N\_{\mathrm{I}}-1} \sum\_{j=0}^{N\_{\mathrm{I}}-1} \langle z(iT\_{\mathbf{c}}, \mathbf{r}) z(jT\_{\mathbf{c}}, \mathbf{r}) \rangle \end{aligned} \tag{A6}$$

According to Equation (A1), *z*(*iTc*, *τ*)*z*(*jTc*, *τ*) can be represented as:

$$
\langle z(iT\_c, \tau)z(jT\_c, \tau) \rangle = \mathbb{C}\_z((i-j)T\_c, \tau) + \langle z(nT\_c, \tau) \rangle \langle z(nT\_c, \tau) \rangle \tag{A7}
$$

$$
\begin{aligned}
\sigma\_Z^2(\tau) &= \mathbb{C}\_Z(0, \tau) \\
&= \frac{1}{N\_1^2} \sum\_{i=0}^{N\_1 - 1} \sum\_{j=0}^{N\_1 - 1} \mathbb{C}\_z((i-j)T\_{c\_\tau}\tau) \\
&= \sum\_{k=-(N\_1 - 1)}^{N\_1 - 1} \frac{N\_1 - k}{N\_1^2} \mathbb{C}\_z(kT\_{c\_\tau}\tau) \\
&= \frac{1}{N\_1} \sum\_{\tilde{n}=-(N\_1 - 1)}^{N\_1 - 1} \mathbb{C}\_z(\tilde{n}T\_{c\_\tau}\tau)
\end{aligned} \tag{A8}
$$

Simultaneous Equations (A3), (A4) and (A8), we can get:

$$\begin{array}{rcl} \frac{1}{N\_{\text{eff}}} & = & \frac{\sigma\_Z^2(\tau)}{Z^2(\tau)}\\ & = & \frac{\sigma\_Z^2(\tau)}{\langle z(nT\_\ell, \tau) \rangle \langle z(nT\_\ell, \tau) \rangle}\\ & = & \frac{\sum\_{\substack{N\_1 = 1\\N\_1 \ge 0}} \left| \mathbb{C}\_{\mathcal{Y}}(\bar{n}\boldsymbol{T}\_\ell, \tau) \right|^2}{N\_1 \left| \mathbb{C}\_{\mathcal{Y}}(0, \tau) \right|^2} \end{array} \tag{A9}$$

**B.**

Assuming that there is no correlation between the signal term and noise term, the complex waveforms of iGNSS-R can include the following four items [18]:

$$\mathcal{C}\_{\mathcal{Y}}(\tilde{n}T\_{\mathfrak{c}},\mathsf{r}) = \mathcal{C}\_{\mathcal{Y},\mathsf{s}}(\tilde{n}T\_{\mathfrak{c}},\mathsf{r}) + \mathcal{C}\_{\mathcal{Y},\mathsf{nd}}(\tilde{n}T\_{\mathfrak{c}},\mathsf{r}) + \mathcal{C}\_{\mathcal{Y},\mathsf{nu}}(\tilde{n}T\_{\mathfrak{c}},\mathsf{r}) + \mathcal{C}\_{\mathcal{Y},\mathsf{nd}\mathfrak{u}}(\tilde{n}T\_{\mathfrak{c}},\mathsf{r}) \tag{A10}$$

Since the noise of the reflected signal is much larger than that of the direct signal, can be assumed. According to [20], the covariance of signal components can be simplified as:

$$\begin{array}{rcl} \mathbb{C}\_{\mathsf{y},\mathsf{s}}(\overline{\mathsf{u}}T\_{\mathsf{c}},\mathsf{r}) &=& \langle y\_{\mathsf{s}}(n\mathsf{T}\_{\mathsf{c}},\mathsf{r})y\_{\mathsf{s}}^{\*}(n\mathsf{T}\_{\mathsf{c}}+\overline{n}\mathsf{T}\_{\mathsf{c}},\mathsf{r})\rangle \\ &=& 2P\_{\mathsf{d}} \| \int \frac{2P\_{\mathsf{f}}\mathsf{G}\_{\mathsf{f}}(\overrightarrow{\mathsf{p}})\mathsf{c}\_{\mathsf{0}}}{4\pi R\_{\mathsf{t}}^{2}(\overrightarrow{\mathsf{p}})R\_{\mathsf{t}}^{2}(\overrightarrow{\mathsf{p}})} \Lambda^{2}(\mathsf{A}\mathsf{r})\mathrm{sinc}^{2}(\Delta f T\_{\mathsf{c}})\exp^{-j2\pi\Lambda f(\overline{n}\mathsf{T}\_{\mathsf{c}})}d^{2}p \end{array} \tag{A11}$$

According to [21], the covariance of the noise component can be simplified as:

$$\begin{array}{l} \mathsf{C}\_{\mathsf{y,nd}}(\vec{n}T\_{\mathsf{c}},\mathsf{r}) = \frac{2P\_{\mathsf{d}}n\_{0}^{\mathsf{d}}}{T\_{\mathsf{c}}} \mathsf{tri}(\frac{\vec{n}T\_{\mathsf{c}}}{T\_{\mathsf{c}}})\\ \mathsf{C}\_{\mathsf{y,ndu}}(\vec{n}T\_{\mathsf{c}},\mathsf{r}) = \frac{Bn\_{0}^{\mathsf{d}}n\_{0}^{\mathsf{u}}}{T\_{\mathsf{c}}} \mathsf{tri}(\frac{\vec{n}T\_{\mathsf{c}}}{T\_{\mathsf{c}}}) \end{array} \tag{A12}$$

Simultaneous Equations (A10)–(A12), *Cy*(*nT <sup>c</sup>*, *<sup>τ</sup>*) can be converted as:

$$\begin{array}{rcl} \mathbb{C}\_{\mathcal{Y}}(\overline{\boldsymbol{n}}\boldsymbol{T}\_{\boldsymbol{c}},\boldsymbol{\tau}) &=& 2\mathcal{P}\_{\mathbf{d}} \iint \frac{2\mathcal{P}\_{\mathbb{C}}\mathcal{C}\_{\mathbf{r}}(\overline{\boldsymbol{p}})\boldsymbol{\tau}\_{0}}{4\pi\mathcal{R}\_{\mathbf{r}}^{2}(\overline{\boldsymbol{p}})\mathcal{R}\_{\mathbf{r}}^{2}(\overline{\boldsymbol{p}})} \Lambda^{2}(\boldsymbol{\Delta}\boldsymbol{\tau}) \times \operatorname{sinc}^{2}(\boldsymbol{\Delta}\boldsymbol{f}\boldsymbol{T}\_{\mathbf{c}}) \exp^{-j2\pi\Lambda f(\overline{\boldsymbol{n}}\boldsymbol{T}\_{\mathbf{c}})}d^{2}\boldsymbol{p} \\ &+ \frac{2\mathcal{P}\_{\mathbb{d}}\mathcal{U}\_{0}^{\mathbb{d}}}{\mathcal{T}\_{\mathbf{c}}} \operatorname{tr}(\frac{\overline{\boldsymbol{n}}\boldsymbol{T}\_{\mathbf{c}}}{\mathcal{T}\_{\mathbf{c}}}) + \frac{Bn\_{0}^{\mathsf{d}}n\_{0}^{\mathsf{u}}}{\overline{\boldsymbol{T}\_{\mathbf{c}}}} \operatorname{tr}(\frac{\overline{\boldsymbol{n}}\boldsymbol{T}\_{\mathbf{c}}}{\mathcal{T}\_{\mathbf{c}}}) \end{array} \tag{A13}$$

Simultaneous Equations (A9) and (A13), the relationship between the reciprocal of the effective incoherent average number and the coherent integration time can be derived as Equation (19).

#### **References**

