Evaluating Impact of Rain Attenuation on Space-borne GNSS Reflectometry Wind Speeds

Abstract

The novel space-borne Global Navigation Satellite System Reflectometry (GNSS-R) technique has recently shown promise in monitoring the ocean state and surface wind speed with high spatial coverage and unprecedented sampling rate. The L-band signals of GNSS are structurally able to provide a higher quality of observations from areas covered by dense clouds and under intense precipitation, compared to those signals at higher frequencies from conventional ocean scatterometers. As a result, studying the inner core of cyclones and improvement of severe weather forecasting and cyclone tracking have turned into the main objectives of GNSS-R satellite missions such as Cyclone Global Navigation Satellite System (CYGNSS). Nevertheless, the rain attenuation impact on GNSS-R wind speed products is not yet well documented. Evaluating the rain attenuation effects on this technique is significant since a small change in the GNSS-R can potentially cause a considerable bias in the resultant wind products at intense wind speeds. Based on both empirical evidence and theory, wind speed is inversely proportional to derived bistatic radar cross section with a natural logarithmic relation, which introduces high condition numbers (similar to ill-posed conditions) at the inversions to high wind speeds. This paper presents an evaluation of the rain signal attenuation impact on the bistatic radar cross section and the derived wind speed. This study is conducted simulating GNSS-R delay-Doppler maps at different rain rates and reflection geometries, considering that an empirical data analysis at extreme wind intensities and rain rates is impossible due to the insufficient number of observations from these severe conditions. Finally, the study demonstrates that at a wind speed of 30 m/s and incidence angle of 30°, rain at rates of 10, 15, and 20 mm/h might cause overestimation as large as ≈0.65 m/s (2%), 1.00 m/s (3%), and 1.3 m/s (4%), respectively, which are still smaller than the CYGNSS required uncertainty threshold. The simulations are conducted in a pessimistic condition (severe continuous rainfall below the freezing height and over the entire glistening zone) and the bias is expected to be smaller in size in real environments.

1. Introduction

The innovative Global Navigation Satellite System Reflectometry (GNSS-R) technique, offering exploitation of existing GNSS signals reflected off the Earth’s surface, has significantly grown in interest as an Earth observation technology to retrieve a variety of geophysical parameters [1,2]. Ocean surface wind monitoring using GNSS-R measurements, as a bistatic radar forward scatterometer, is one of the latest promising applications. TechDemoSat-1 (TDS-1), launched in 2014, fully demonstrated the capabilities of GNSS-R in retrieving the state of the ocean, and consequently, in obtaining high-quality winds [3]. The Cyclone Global Navigation Satellite System (CYGNSS) was launched in December 2016 and provides wind data with an unprecedented sampling rate. CYGNSS consists of eight microsatellites in the inclined orbits [4,5]. In addition, further ideas and upcoming GNSS-R missions are being developed [6,7,8].

GNSS-R satellite missions can provide a sufficient observation frequency for ocean surface monitoring, exploiting the high and growing number of GNSS transmitters and the easily employable low-cost, low-mass, and low-power GNSS receivers. In addition, it is known that GNSS L-band signals are affected by a much lower level of attenuation and scattering by clouds or raindrops in the atmosphere, compared to those from the conventional Ku- and C-band scatterometers. The general reliability of wind speed products during rainfalls is documented [9]. These advantages can be potentially a solution to overcome the limitations in capturing and forecasting severe weather events. Weather models performance degrade during extreme events due to the lack of high spatiotemporal-resolution data and the obscurity of traditional remote sensing instruments during severe weathers. Consequently, improving severe weather forecasting is the overall objective of CYGNSS and the mission studies tropical cyclone inner core process.

Rain can affect ocean monitoring with scatterometers in two ways. First, raindrops impinging the ocean surface alter the wind-induced radar signature, which, in case of GNSS-R scatterometry, is the quasi-specular reflections from multiple facets. They can create rings, stalks, and crowns [10] which alter the surface roughness, and consequently, the wind retrieval quality. Rain splash effect on GNSS-R observations is discussed [11] and a decrease in Bistatic Radar Cross Section (BRCS) at low wind speeds (≈<5 m/s) is demonstrated both with TDS-1 data analysis and based on model simulations. At high enough winds (≈>5 m/s), the intensity of forward quasi-specular scattering is controlled by surface gravity waves with lengths larger than several wavelengths of the reflected signal. As a result, the scatterometric GNSS-R measurements are insensitive to the surface modifications by raindrops impinging on the ocean at high winds.

As the second effect, rain attenuates signals passing through the atmosphere, the phenomenon which this study is focused on. Existing raindrops and water vapor in the atmosphere cause attenuation in Electromagnetic (EM) waves due to the absorption and scattering [12]. Absorption is explained by the intersection of rain cells and propagation path of EM waves and the absorbed energy is transferred into heat. Deeper absorption happens when the rain cells fill a major part of the Fresnel’s ellipsoid along the signal path. In case of large enough signal frequency, such as 12 to 18 GHz (Ku-band), as the precipitation rate increases, the signal received less at the receiver is scattered from the ocean surface and the majority of the radiation is scattered by the rain layer. These effects significantly degrade the performance of scatterometers operating at Ku-band frequencies [13].

Although the rain attenuation can cause small-scale effects on GNSS signals, the impact on the diverse applications can be different in size. A small change in the observable cannot be always interpreted as an insignificant bias in the final products but it depends on how the product responses to the observable change. The dependency of the wind speed to a GNSS-R observable is translated by a Geophysical Model Function (GMF) which is a forward operator mapping the observable to the product. Mathematically, assuming that the wind speed GMF is a linear function, we can consider the simplified wind speed retrieval problem as:

$$\mathit{G}\overrightarrow{x}={\overrightarrow{U}}_{10}$$

(1)

where $\mathit{G}$ is a $m\times n$ GMF matrix and $\overrightarrow{x}$ and ${\overrightarrow{U}}_{10}$ are the n- and m-dimensional observable and wind speed vectors, respectively. In the case of perturbations, we have $G(\overrightarrow{x}+\overrightarrow{\delta }x)={\overrightarrow{U}}_{10}+\delta {\overrightarrow{U}}_{10}$, form which we can estimate the relative error in the wind speed as:

$$\Vert \delta {\overrightarrow{U}}_{10}\Vert /\Vert {\overrightarrow{U}}_{10}\Vert \ge \Vert \mathit{G}\Vert \Vert {\mathit{G}}^{-1}\Vert (\Vert \overrightarrow{\delta }x\Vert /\Vert \overrightarrow{x}\Vert )$$

(2)

which shows that the final error in the wind speed is a function of condition number $\mu$ which reads:

$$\mu \left(G\right)=\Vert \mathit{G}\Vert \Vert {\mathit{G}}^{-1}\Vert$$

(3)

Finally, a small change in the observable can cause large biases in the output for problems with a large condition number, so-called ill-conditioned problems.

Normally, in GNSS-R wind speed retrievals, the GMF is a nonlinear function. The condition number for nonlinear functions, unlike the discussed linear problem, is not a constant value, and varies with the point over the domain of the function. The condition number of a nonlinear but differentiable function f at point x is $x{f}^{\prime }\left(x\right)/f\left(x\right)$. Hence, a TDS-1 wind speed $U_{10}$ GMF, taking the Bistatic Radar Cross Section (BRCS) ${\sigma }^0$ as the input argument, $U_{10}=9042.24e^{-0.62{\sigma }^0}+0.99$ [9], can be considered and the associated condition number is shown in Figure 1. Accordingly, at lower values of ${\sigma }^0$, corresponding to higher wind speeds, the condition number is at higher levels. So, at a small enough ${\sigma }^0$, that is at high enough wind speeds, the GMF tends to an ill-conditioned function.

According to the above arguments, small changes in the GNSS-R observables might cause a considerable bias in the wind speed at high-value regions, which signifies the need to a characterization investigation, considering the main objective of CYGNSS. However, conducting data analysis at extreme winds and rains is impossible due to the lack of a sufficient number of observation from extreme events. This study characterizes the effect simulating the delay-Doppler Maps (DDM) and ${\sigma }^0$, as described in Section 2. Counting rain attenuation effects is explained in Section 3, which also reports the numerical results. Finally, Section 4 discusses the outcomes and provides the concluding remarks.

2. Simulating DDMs

In this study, the simulations are based on the bistatic radar equation which reads [14]:

$$\langle \mid Y(\tau ,D){\mid }^2\rangle =\frac{P_t{G}_t{\lambda }^2{T}_i^2}{{\left(4\pi \right)}^3}\int {\int }_A\frac{G_r{\mathsf{\Lambda }}^2\left(\tau \right)S^2\left(f\right)}{R_t^2{R}_r^2}{\sigma }^{0}dA$$

(4)

where $\tau$ and D are the time delay and the frequency offset, respectively, and $\langle \mid Y(\tau ,D){\mid }^2\rangle$ is the ensemble mean of the correlation as a function of these parameters. $P_t$ is the transmitter power, $G_t$ and $G_r$ are the transmitter and receiver antenna gains, respectively. $\lambda$ is the carrier wavelength, $T_i$ is coherent integration time, $R_t$ and $R_r$ are respectively the transmitter to specular point and specular point to receiver ranges, A represents the scattering area (the glistening zone), ${\sigma }^0$ is the bistatic radar cross section, ${\mathsf{\Lambda }}^2$ and $S^2$ are the triangular and sinc functions, respectively. Different studies propose DDM simulation algorithms. For instance, one can refer to [14,15,16,17] for further information. Nevertheless, the simulation algorithm conducted in this study is here explained in details for the sake of clarity. For simulating the DDMs, the sea surface is modeled firstly as a function of wind speed, which is the dominant factor controlling the surface roughness. Other oceanic phenomena such as rain splash effects and swell, which might change the roughness and, consequently, affect the DDMs, are not investigated in this study. In case of interest, one can refer to [11] which discusses the rain splash effect intensively. Then, the reflection geometry and the specular point position is determined. Finally, the power distribution in a spatial domain is computed and is then mapped to the delay-Doppler domain.

2.1. Ocean Surface Modeling

The BRCS describes the strength of the signal after reflection from a particular point on the rough surface and scattered in the direction of the receiver’s antenna. The intensity of the propagated signal is a function of roughness and, as a result, the BRCS describes the glistening zone being rough. In these numerical simulations, we use the equation derived from the well-known geometric-optics limit of the Kirchhoff approximation [14,18]:

$${\sigma }_0=\mid {\mathfrak{R}}^2\mid {(\mid \overrightarrow{q}\mid /q_z)}^{4}P({\overrightarrow{q}}_{\perp }/q_z)$$

(5)

where $\mathfrak{R}$ is the complex Fresnel coefficient which can be determined by the complex dielectric constant of sea water, signal polarization, and local incidence angle [14]. $\overrightarrow{q}={\overrightarrow{q}}_{\perp }+q_z\widehat{z}$ is the scattering vector, which can be obtained from the locations of the transmitter and receiver satellites and the specular point. $P({\overrightarrow{q}}_{\perp }/q_z)$ is the Probability Density Function (PDF) of the slopes of the large-scale sea surface components at wavenumbers much higher than $k_{*}=(2\pi /3\lambda )cos\left({\theta }_{inc}\right)$, where ${\theta }_{inc}$ is the incidence angle.

For numerical simulations, the PDF for a given slope $s={\overrightarrow{q}}_{\perp }/q_z$, can be obtained using the Gaussian statistics of anisotropic slopes [14]:

$$P\left(s\right)=\frac{1}{2\pi {\sigma }_{s_x}{\sigma }_{s_y}\sqrt{1-b_{x,y}^2}}exp\left(-\frac{1}{2(1-b_{x,y}^2)}\left(\frac{s_x^2}{{\sigma }_{s_x}^2}-2b_{x,y}\frac{s_x{s}_y}{{\sigma }_{s_x}{\sigma }_{s_y}}+\frac{s_y^2}{{\sigma }_{s_y}^2}\right)\right)$$

(6)

where $s_x$ and $s_y$ are the upwind and crosswind ocean slope components and ${\sigma }_{s_x}$ and ${\sigma }_{s_y}$ denote the upwind and crosswind Mean Square Slope (MSS) components, respectively. The slope variances and correlations can be determined from a ocean surface elevation spectrum $W\left(\overrightarrow{k}\right)$ with an integration over wavenumbers smaller than $k_{*}$, which can be in turn modeled using different spectra models such as the one proposed in [19]. One can refer to [14] for the integrals deriving the wind-dependent slope variances and correlations which are valid for any arbitrary spectrum $W\left(\overrightarrow{k}\right)$. Instead, the direct equations are here used, deriving the MSS components [20]:

$${\sigma }_{s_x}=3.16\times {10}^{-3}{U}_{10}$$

(7)

$${\sigma }_{s_y}=0.003+1.92\times {10}^{-3}{U}_{10}$$

(8)

where $U_{10}$ is the wind speed at 10 m above the ocean surface. For simplicity, the spectrum is here supposed symmetrical with respect to a wind direction which is in turn assumed along the major x- or y-axis ($b_{x,y}=0$). Finally, ${\sigma }^0$ is computed for different scenarios leading to simulated DDMs as described in the following subsections.

The type of scattering considered here and modeled by Equation (4) is the quasi-specular scattering, which takes place on oceans induced with intense enough winds ($\ge \approx$4 m/s), in other words, for rough surfaces with a high enough (≫1) Rayleigh parameter, $R_a=(2\pi /\lambda )hcos\left({\theta }_{inc}\right)$, where h is the root-mean-square of surface heights. For the case of scattering from insufficiently rough oceans (induced by weak winds), the mechanism changes the rough surface scattering dominated by or combined with the coherent specular reflection component which is absent at higher winds (at a large Rayleigh parameter). This type of scattering and its transition to the regime of strong diffuse scattering is discussed and modeled in [21]. In this study, the scattering mechanism change at low winds is neglected for simplicity since the main concern is the attenuation effect at high winds where the wind speed GMF has a high condition number and might be potentially sensitive to BRCS small changes, as discussed in the introduction.

2.2. Reflection Geometry

The reflection geometry is defined with locations of transmitter and receiver satellites and the specular point on the ocean surface. For more realistic simulation, eight real TDS-1 reflectometry events, at different incidence angles of 0°, 10°, 20°, 30°, 40°, 50°, 60°, and 70° are collected as the case scenarios. Then, the direct information on the position of TDS-1 and GPS satellites are used as the coordinates of the receiver and transmitter, respectively. Appendix A reports on the position and velocity of the satellites in the Earth-centered Earth-Fixed (ECEF) reference frame, in each of eight case scenarios, in

Table A1. Receiver and transmitter satellite coordinates (m) in the Earth-centered Earth-fixed (ECEF) system.

Incidence Angle	Receiver			Transmitter
	$\mathit{x}$	$\mathit{y}$	$\mathit{z}$	$\mathit{x}$	$\mathit{y}$	$\mathit{z}$
0°	3,401,792.25	5,401,767.00	2,892,517.25	12,928,387.67	20,596,486.63	11,085,709.07
10°	5,432,997.00	−2,902,608.00	−3,353,593.75	19,042,175.45	−9,244,237.63	−15,715,969.40
20°	6,467,244.50	−2,193,562.50	1,560,833.37	21,924,510.86	−7,603,784.91	13,351,620.42
30°	−5,847,129.50	3,453,629.50	−1,752,086.87	−25,527,012.40	6,709,641.30	2,884,300.00
40°	−2,831,923.25	3,264,313.00	−5,519,979.00	−19,885,837.37	15,623,426.78	−8,200,966.13
50°	−6,043,847.00	3,411,350.00	985,065.06	−15,004,963.22	16,008,724.98	−14,947,308.25
60°	5,958,662.00	−2,340,375.75	2,830,701.50	11,365,226.82	10,898,201.04	21,366,130.62
70°	−5,196,265.50	−2,030,612.37	−4,253,424.50	12,082,232.34	−9,818,137.67	−21,484,129.2

and

Table A2. Receiver and transmitter satellite velocities (m/s) in the Earth-centered Earth-fixed (ECEF) system.

Incidence Angle	Receiver			Transmitter
	$\dot{\mathit{x}}$	$\dot{\mathit{y}}$	$\dot{\mathit{z}}$	$\dot{\mathit{x}}$	$\dot{\mathit{y}}$	$\dot{\mathit{z}}$
0°	3056.24	1719.31	−6776.05	−1186.06	−710.54	2746.57
10°	−3928.48	169.99	−6531.11	−963.524	1811.42	−2215.26
20°	1064.68	−2079.62	−7270.98	1641.96	261.90	−2546.16
30°	−764.41	2345.64	7218.34	−412.20	−208.37	−3149.80
40°	5362.44	−2980.10	−4512.57	−1080.42	87.56	2862.90
50°	1712.27	889.13	7385.07	−1836.21	465.19	2367.17
60°	2208.92	−2655.07	−6812.81	−1438.98	2318.73	−362.66
70°	−4829.46	27.42	5889.65	2312.36	1455.58	639.18

.

2.2.1. Calculating the Specular Point Position

Once the positions of the satellites are known, the specular point position can be determined. To this aim, the algorithm proposed in [22] is used. The specular point position is calculated so that three conditions are met: the angle between the incoming and reflected signal with respect to the surface normal is equal, the total path between the transmitter, specular point, and the receiver is the minimum, and finally, the specular point relies on the WGS84 Earth geoid. For a more precise simulation, when needed, local variations in the geoid height can be taken into account (one can refer to [23] for example). Then, the signal path magnitude as a function of the specular point position reads [22]:

$$P(\overrightarrow{S})=\mid (\overrightarrow{T}-\overrightarrow{S})+(\overrightarrow{R}-\overrightarrow{S})\mid$$

(9)

where $\overrightarrow{S}$, $\overrightarrow{T}$, and $\overrightarrow{R}$ denote the specular point, transmitter and receiver position vectors. Since the transmitter and receiver locations are known, the specular point location can be obtained with an iterative minimizing approach as the path equation is non-linear. Taking the partial differential of the specular point $\overrightarrow{S}$ with respect to x, y, and z for the purpose of minimizing this path:

$$\begin{array}{r}dS=\left(\frac{T_x-S_x}{\sqrt{{(T_x-S_x)}^2+{(T_y-S_y)}^2+{(T_z-S_z)}^2}}+\frac{R_x-S_x}{\sqrt{{(R_x-S_x)}^2+{(R_y-S_y)}^2+{(R_z-S_z)}^2}}\right)dx\\ +\left(\frac{T_y-S_y}{\sqrt{{(T_x-S_x)}^2+{(T_y-S_y)}^2+{(T_z-S_z)}^2}}+\frac{R_y-S_y}{\sqrt{{(R_x-S_x)}^2+{(R_y-S_y)}^2+{(R_z-S_z)}^2}}\right)dy\\ +\left(\frac{T_z-S_z}{\sqrt{{(T_x-S_x)}^2+{(T_y-S_y)}^2+{(T_z-S_z)}^2}}+\frac{R_z-S_z}{\sqrt{{(R_x-S_x)}^2+{(R_y-S_y)}^2+{(R_z-S_z)}^2}}\right)dz\end{array}$$

(10)

We have the change in $\overrightarrow{S}$ converging to the minimum path step as:

$$\mathsf{\Delta }\overrightarrow{S}=\frac{(\overrightarrow{T}-\overrightarrow{S})}{\mid (\overrightarrow{T}-\overrightarrow{S})\mid }+\frac{(\overrightarrow{R}-\overrightarrow{S})}{\mid (\overrightarrow{R}-\overrightarrow{S})\mid }$$

(11)

After the initial guess, which is simply the point below the receiver satellite and on WGS84 surface reference, the iteration is carried out as we have ${\overrightarrow{S}}_{temp,i}={\overrightarrow{S}}_{temp,i-1}+\mathsf{\Delta }{\overrightarrow{S}}_i$, where i is the iteration step and $temp$ indicate the intermediate being of the value. To ensure the specular point lies on the ocean surface (on the WGS84), ${\overrightarrow{S}}_{temp}$ is scaled by the radius of the Earth, $r/\mid {\overrightarrow{S}}_{temp}\mid$, where:

$$r=a_{WGS84}\sqrt{\frac{1-e_{WGS84}^2}{1-(e_{WGS84}^2+cos\left(\lambda \right))}}$$

(12)

in which $a_{WGS84}$ and $e_{WGS84}$ are the semimajor axis and eccentricity of the WGS84 reference ellipsoid. So, the specular point can be obtained as:

$$\overrightarrow{S}=r\frac{{\overrightarrow{S}}_{temp}}{\mid {\overrightarrow{S}}_{temp}\mid }$$

(13)

when $d\overrightarrow{S}$ is smaller than a specified tolerance after several iterations, it is considered determined.

2.2.2. Transformation of the Vectors to the Model Reference Frame

To ease the simulation process, all of the position and velocity vectors are transformed to a model reference frame. The model reference frame is defined so that its origin is at the center of the Earth. The specular point is along the z-axis, and the x–y plane is parallel to the plane tangent to the surface at the specular point. The system is right-handed and the x-axis is toward the receiver satellite trajectory on the x–y plane. This reference frame is shown in Figure 2.

To obtain the vectors in the model reference frame, they are first transferred to a local East, North, Up (ENU) coordinate system centered at the specular point. This transformations is as follows:

$$\left[\begin{array}{c}{x}_{enu}\\ y_{enu}\\ z_{enu}\end{array}\right]=\left[\begin{array}{ccc}-sin{\lambda }_s& cos{\lambda }_s& 0\\ -sin{\varphi }_{s}cos{\lambda }_s& -sin{\varphi }_{s}sin{\lambda }_s& cos{\varphi }_s\\ cos{\varphi }_{s}cos{\lambda }_s& cos{\varphi }_{s}sin{\lambda }_s& sin{\varphi }_s\end{array}\right]\left[\begin{array}{c}{x}_{ecef}\\ y_{ecef}\\ z_{ecef}\end{array}\right]$$

(14)

where ${\lambda }_s$ and ${\varphi }_s$ denote the geodetic longitude and latitude of the specular point, respectively. The ENU and the model refrence frame have equivalent z components, while the the horizontal orientation can be different. With a rotation of the coordinate system with respect to z-axis, ENU coordinates are transformed to the model reference frame. Then, the transformation of the horizontal components are first conducted with calculating the horizontal angular difference between the two systems:

$$\mathsf{\Delta }\theta ={\theta }_h-{\theta }_r$$

(15)

where:

$${\theta }_h={tan}^{-1}(x_{enu}/y_{enu})$$

(16)

$${\theta }_h={tan}^{-1}\left(\frac{x_{t,enu}-x_{r,enu}}{y_{t,enu}-y_{r,enu}}\right)$$

(17)

denoting the coordinates components of the reciever and trasnmitter satellite with r and t subscripts, respectively. Finally, the vectors in the model reference frame can be easily obtained:

$$\left[\begin{array}{c}{x}_m\\ y_m\\ z_m\end{array}\right]=\left[\begin{array}{c}\sqrt{x_{enu}^2+y_{enu}^2}cos\mathsf{\Delta }\theta \\ \sqrt{x_{enu}^2+y_{enu}^2}sin\mathsf{\Delta }\theta \\ z_{enu}\end{array}\right]$$

(18)

and in the following, all of the vectors are designated in the model reference frame.

2.2.3. Establishment of the Glistening Zone

The glistening zone is considered as the area surrounding the specular point. So, a spherical plane is considered so that the curvature of Earth is taken into account. To this aim, a grid is etablished whose elements are designated in a 2D angular frame $({\gamma }_1,{\gamma }_2)$. The components ${\gamma }_1$ and ${\gamma }_2$ are the angular differences with respect to the z-axis of the model reference frame (intersecting the specular point), while the endpoints of the rays are at its origin. The boundary of the glistening zone can be considered as $\gamma =R_g/r$, where $R_g$ is the radius of the circular glistening zone. Finally, the grid elements can be obtained in the model reference frame by rotation of the specular point position ${\overrightarrow{S}}_m=[0,0,1]$ with respect to y-axis and x-axis as large as ${\gamma }_1$ and ${\gamma }_2$, respectively:

$${\overrightarrow{g}}_m({\gamma }_1,{\gamma }_2)={\overrightarrow{S}}_m{M}_1\left({\gamma }_1\right)M_2\left({\gamma }_2\right)$$

(19)

where,

$$M_1\left({\gamma }_1\right)=\left[\begin{array}{ccc}-cos{\gamma }_1& 0& sin{\gamma }_1\\ 0& 1& 0\\ sin{\gamma }_1& 0& cos{\gamma }_1\end{array}\right]$$

(20)

$$M_2\left({\gamma }_2\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos{\gamma }_2& -sin{\gamma }_2\\ 0& sin{\gamma }_2& cos{\gamma }_2\end{array}\right]$$

(21)

2.3. Power Distribution Calculation

In this step, the received power is calculated using Equation (4). However, the direct calculation of the power map in the delay-Doppler domain is not possible unless it is first obtained in the spatial domain (over the established grid over ocean surface) and then mapped to the delay-Doppler domain afterward.

2.3.1. Power Distribution in the Spatial Domain

For each grid element we need the corresponding delay and Doppler shift values to calculate the received power scattered form that surface element. The delay refers to the time delays of the signal scattered from different grid elements. The delay value at every grid point on the ocean surface can be obtained using Equation (9) and the path delay at the specular point (the minimum path delay) can be used to reference other paths across the grid. To be more specific, the C/A delay for each grind position vector $\overrightarrow{g}$ is obtained as:

$$\tau \left(\overrightarrow{g}\right)=\frac{f_{C/A}}{c}P\left(\overrightarrow{g}\right)$$

(22)

where $f_{C/A}=1.025e^6$ Hz and $c=2.99792458e^8$ are the frequency of C/A code and the speed of light in a vacuum, respectively, and P is given in Equation (9). Finally, $\mathsf{\Delta }\tau =\tau \left(\overrightarrow{g}\right)-\tau (\overrightarrow{S})$ is the delay with respect to the specular point.

The relative motion between the transmitter, receiver and grid elements cause the Doppler shift. In an arbitrary case, a signals Doppler frequency between two points $\overrightarrow{A}$ and $\overrightarrow{B}$ with velocity vectors of ${\overrightarrow{A}}_v$ and ${\overrightarrow{B}}_v$, respectively, is as follows:

$$D_{\overrightarrow{A}\overrightarrow{B}}=\frac{\mid ({\overrightarrow{A}}_v-{\overrightarrow{B}}_v)·\frac{(\overrightarrow{A}-\overrightarrow{B})}{\mid \overrightarrow{A}-\overrightarrow{B}\mid }\mid f}{c}$$

(23)

where f is the signal frequency. In GNSS-R, the signal is received after being reflected and the total Doppler frequency at grid element $\overrightarrow{g}$, after adding the rate of the receiver clock drift $d_{clk}$, reads:

$$D\left(\overrightarrow{g}\right)=D_{Tg}+D_{gR}+d_{clk}$$

(24)

which can be rewritten as:

$$D\left(\overrightarrow{g}\right)=\frac{\mid ({\overrightarrow{g}}_v-{\overrightarrow{T}}_v)·{\overrightarrow{u}}_{Tg}\mid f_{L1}}{c}+\frac{\mid ({\overrightarrow{R}}_v-{\overrightarrow{g}}_v)·{\overrightarrow{u}}_{gR}\mid f_{L1}}{c}+d_{clk}$$

(25)

where $\overrightarrow{u}$ is the unit vector between the two points shown as the subscripts and $f_{L1}=1575.42e6$ Hz is GPS L1 carrier signal frequency. The receiver clock drift is usually determined as an unknown parameter in positioning solutions. Nevertheless, this value is assumed zero in simulations of this study. In addition, the velocity of sea surface (grid elements) can be assumed equivalent to zero since its magnitude is much smaller than the transmitter and receiver satellites. This calculation is carried out at every grid element to map the scattered signal frequency over the glistening zone. Similar to delay computations, the Doppler frequency at the specular point is used to reference other frequencies across the grid as $\mathsf{\Delta }D=D\left(\overrightarrow{g}\right)-D(\overrightarrow{S})$.

2.3.2. Power Distribution in Delay-Doppler Domain

Once the power distribution is transformed from the spatial domain to the delay-Doppler domain, the simulated DDM is generated. To this aim, the relationship between the two domains must be first considered. Figure 3 schematically illustrates an idealized DDM in both the spatial and delay-Doppler domain as well as the equi-Doppler curves, shown in blue, and the equi-range lines, ellipses in black. The right panel shows the corresponding DDM with the well-known horseshoe shape. The surface elements in the spatial coordinate domain, specified at the intersection of equi-range and equi-Doppler curves, and DDM pixels in the delay-Doppler domain are connected with a relation which is not one-to-one. In other words, every pixel of the DDM is not mapped to exactly one element of the spatial domain. Each DDM pixel is proportional to scattered power from the pair of pixels located symmetrically with respect to the intersection of the incidence plane with the ocean surface, which is the y-axis in the illustrated case. This is valid for every DDM pixel, except the pixel with maximum value, which normally corresponds to the specular point and the pixels with zero value (or a value close to the noise-floor), the dark blue area of the DDM in the right panel, which correspond to non-intersecting equi-range and equi-Doppler lines in the spatial domain. This non-one-to-one relation causes the so-called ambiguity problem.

Dealing with the quasi-specular scattering from the glistening zone and the discussed ambiguity, the received power can be spatially filtered with the help of the Woodward Ambiguity Function (WAF), $\chi (\mathsf{\Delta }f,\mathsf{\Delta }\tau )\approx S\left(\mathsf{\Delta }f\right)\mathsf{\Lambda }\left(\mathsf{\Delta }\tau \right)$, applied as follows:

$$\langle \mid Y(\mathsf{\Delta }\tau ,\mathsf{\Delta }D){\mid }^2\rangle ={\chi }^2(\mathsf{\Delta }\tau ,\mathsf{\Delta }D)*\mathsf{\Sigma }(\mathsf{\Delta }\tau ,\mathsf{\Delta }D)$$

(26)

where * represents the two dimensional convolution operation and:

$$\mathsf{\Sigma }\left(\overrightarrow{g}\right)=\frac{T_i^2{D}^2\left(\overrightarrow{g}\right){\sigma }^0\left(\overrightarrow{g}\right)ds}{4\pi R_t^2\left(\overrightarrow{g}\right)R_r^2\left(\overrightarrow{g}\right)}$$

(27)

in which $ds$ is the area of the grid element. To accelerate the numerical simulations, Equation (26), with the help of Fourier transformation $\mathcal{F}$, can be rewritten as [15]:

$$\mathcal{F}\left[\langle \mid Y(\mathsf{\Delta }\tau ,\mathsf{\Delta }D){\mid }^2\rangle \right]=\mathcal{F}\left[\mathsf{\Sigma }(\mathsf{\Delta }\tau ,\mathsf{\Delta }D)\right].\mathcal{F}\left[{\chi }^2(\mathsf{\Delta }\tau ,\mathsf{\Delta }D)\right]$$

(28)

Finally, using the procedure described in this section, the DDMs are simulated for each scenario reported in Table A1 and Table A2. Figure 4 shows a simulated and measured DDM onboard TDS-1.

3. Rain Attenuation

After simulating the DDMs, rain attenuation effects on them is applied as described in the following subsections.

3.1. Attenuation Model

For practical purposes, the specific attenuation ${\gamma }_R$, for a given rain rate R, can be calculated simply using the power-law relationship:

$${\gamma }_R=k{R}^{\alpha }$$

(29)

The coeficient values k and $\alpha$ are a function of the frequency f. They are obtainable from equations proposed by the International Telecommunication Union-Radio communation sector (ITU-R), specifically ITU-R P.838-3, which have been developed from curve-fitting to power-law coefficients derived from scattering calculations [25]:

$${log}_{10}{k}_{H/V}=\sum _{j=1}^4\left(a_{j}exp\left[-{\left(\frac{{log}_{10}f-b_j}{c_j}\right)}^2\right]\right)+m_k{log}_{10}f+c_k$$

(30)

$${\alpha }_{H/V}=\sum _{j=1}^5\left(a_{j}exp\left[-{\left(\frac{{log}_{10}f-b_j}{c_j}\right)}^2\right]\right)+m_{\alpha }{log}_{10}f+c_{\alpha }$$

(31)

The coefficient are depended on the signal polarization as shown with subscripts H or V for horizontal and and vertical polarization, respectively. Values of $a_j$, $b_j$, $c_j$, $m_{k/\alpha }$, and $c_{k/\alpha }$ for each polarization are given in Appendix B from ITU-R P.838-3. Coefficients k and $\alpha$, for a circular polarization, can be obtained from those for horizontal and vertical polarization as follows:

$$k=[k_H+k_V+(k_H-k_V){cos}^2\theta cos2{\tau }_t]/2$$

(32)

$$\alpha =[k_H{\alpha }_H+k_V{\alpha }_V+(k_H{\alpha }_H-k_V{\alpha }_V){cos}^2\theta cos2{\tau }_t]/2k$$

(33)

where $\theta$ is the path elevation angle and ${\tau }_t$ is the polarization tilt angle relative to the horizon (=45° for circular polarization). At GPS L1 frequency (f = 1575 MHz), the values for k and $\alpha$ are $24.312e^{-5}$ and $0.9567$, respectively [26]. Finally, the rain attenuation $\left(G_{rain}\right)$ is computed with a forward model as follows [26]:

$$G_{rain}=exp(-{\gamma }_R{h}_f(csc{\theta }_t+csc{\theta }_r))$$

(34)

where $h_f$ is the freezing height in km (assumed 6 km in this study), and ${\theta }_t$ and ${\theta }_r$ are the elevation angles to the transmitter and receiver satellites, respectively. The elevation angles of both satellites are assumed equal, as would be the case at the secular point. In addition, the rain rate is assumed constant from the surface up to the freezing height. In Equation (34), $G_{rain}$ is the power out as percentage of the power in, which is simply convertible to attenuation in dB ($=10{log}_{10}\left(G_{rain}\right)$).

3.2. The Effect on GNSS-R Winds

Finally, the rain effects are calculated as described previously and summarized in Figure 5. The simulation is conducted for rain rates up to 20 mm/h effecting the observation at eight different TDS-1 case scenarios with incidence angles of 0°, 10°, 20°, 30°, 40°, 50°, 60°, and 70° as reported in Table A1 and Table A2.

As a means to extract the wind information, retrieval algorithms are developed based on different quantities obtained from the DDMs as the input into the empirical GMFs. The GMF is in turn derived from a large collocation study between the observed DDM measurement together with wind data from other sources such as scatterometers, situ buoys and Numerical Weather Prediction (NWP) model data.

Among diverse observables, the BRCS computed from the signal-to-noise ratio (SNR) of the correlated power around the specular point is used to retrieve wind speed form TDS-1 measurements [3] and the average normalized BRCS around the specular point is used in the CYGNSS wind speed algorithm [27]. As a result, the rain attenuation effects on this parameter and, consequently, the derived wind speed is investigated. Figure 6 demonstrates the rain attenuation effects on the BRCS at different rain rates and incidence angles, which is followed by Figure 7 showing the eventual modification on the derived BRCS and wind speeds. The differences indicated by Figure 7 are at the incidence angle of 30° which has the highest occurrence probability in TDS-1 dataset.

4. Discussion and Conclusions

Rain attenuation, depending on the signal frequency, can affect the satellite scatterometric wind retrievals. This effect on the novel space-borne GNSS-R, providing wind information over oceans, needs to be characterized, since one of the main objectives of this technique is providing unaffected data during severe weather. Although the L-band GNSS signals are not attenuated at high levels by the raindrops in the atmosphere, the high condition number of GMFs at high wind speed ranges may result in a considerable bias in the resultant wind speed. In this study, the signal attenuation impact on the BRCS ${\sigma }^0$ at the specular point as one of the important quantities derived from DDMs, as well as the bias in the wind speed, is evaluated. Empirical GNSS-R data analysis is impossible due to the lack of a sufficient number of observations at extreme wind speeds and rain rates. Consequently, this study is conducted simulating the DDMs from a satellite with orbital properties similar to TDS-1.

Figure 6 shows that the modification due to the rain attenuation is slight in at all of the incidence angles. As visualized in this figure (and as reported in Figure 1 in terms of the condition number), at wind speeds below 5–10 m/s, the BRCS change over the wind speed is very gradual (d${\sigma }^0/$d$U_{10}$ is small) and the wind speed is barely responsive to BRCS changes. As a result, although the modifications in BRCS are relatively larger at low winds, the resultant bias is negligibly small (as shown in Figure 7). At higher wind speeds, the slope of the curves in Figure 6 becomes very steep and small changes in BRCS correspond to relatively large changes in the derived wind speed. So, the component of wind speed error due to measurement error will be larger here. At typical precipitations (<5 mm/h), rain attenuation is negligibly ineffective so that at a wind speed of 30 m/s and incidence angle of 30°, the bias is smaller than ≈0.35 m/s (1%).

The numerical results demonstrate that rain rates of 10, 15 and 20 mm/h, might cause overestimation as large as ≈ 0.65 m/s (2%), 1.00 m/s (3%), and 1.3 m/s (4%), respectively, at a wind speed of 30 m/s and incidence angle of 30°. These statistics can be considerable if they exceed the required uncertainty threshold. This threshold is established considering the intrinsic errors in the wind speed retrieval algorithms, which represent the errors that would be still present even if the measurements were perfect. In practice, perfectly modeling the dependence of measured DDMs on non-wind variables is impossible and the algorithms are at a level of simplification. Such unavoidable error sources are combined with those from other phenomena, including the rain attenuation effect. The uncertainty threshold for CYGNSS mission is 2 m/s at wind speeds below 20 m/s and ±10% above 20 m/s. The errors shown in Figure 7 are always much smaller than this threshold, and thereby might not be statistically significant, even at the incidence angle of 70°, which is the worst case scenario shown in Figure 6.

In addition, according to [28], the average rain rate during hurricanes is 3.9 mm/h and majority of the area (>50 %) are covered with a rainfall rate of 0.25–6.25 mm/h, estimated from composites of rain gauge data measuring rain rate on Pacific islands during passage of typhoons over a 21-year period. Only over 4% of the total covered area the rain rate is larger than 19 mm/h, which may be for a very short time (few minutes). In hurricane events, extremely large rain rates happen in connective cells which are highly localized and moving and evolving rapidly. So, the severe rainfalls happen for few minutes in fast-moving cells. However, the simulation here counts the rain attenuation effects very pessimistically. It considers constant and continuous precipitation along the entire path of the direct and reflected signal below the freezing height and over the entire glistening zone. The ambiguity function (explained in Section 2.3.2) smooths the effect of rain attenuation over an area of ≈25 km. So, the effect of extreme rainfalls over small local regions can be significantly reduced. As a result, the effect of rainfall during hurricanes can be noticeably smaller in reality than the reported statistics.

It is recently demonstrated that rain splash altering ocean roughness can be considerable at low wind speeds. The characterization of the signal attenuation effect in this study complements the previous discussions and helps scientists obtain a better understanding in the development of models correcting rain effects on wind products and potentially in detecting rainfalls over oceans with GNSS-R measurements. Generally, this study demonstrates that rain attenuation effects can be ignored, considering the current technical and theoretical level of space-borne GNSS-R maturity, and the technique shows the required insensitivity to atmospheric rain attenuation for providing high-quality data from tropical storms.