An Instrument Error Correlation Model for GNSS Reflectometry

1. Introduction

Launched in 2016, the Cyclone Global Navigation Satellite System (CYGNSS) mission provides a suite of land and ocean surface sensing products through a technique known as Global Navigation Satellite System reflectometry (GNSS-R). With GNSS-R, the eight CYGNSS satellites act as receivers in a bistatic radar system, where GPS satellites are the transmitters, and the Earth surface is the target of interest.

The CYGNSS mission aims to improve tropical cyclone prediction by measuring ocean surface wind speeds over the tropics [1]. Other satellite observations are impaired by difficulty in collecting windspeeds underneath thick cloud layers in tropical cyclones, and in-situ measurements often are challenging and dangerous to collect due to extreme conditions within tropical cyclones. By exploiting GPS broadcasts in the L-band where there is limited attenuation from thick rain and cloud layers, CYGNSS can collect observations through mature tropical cyclones.

With eight satellites in equatorial low-Earth orbit, CYGNSS can collect global observations with a median revisit time of 3.8 hours. However, the sampling characteristics of CYGNSS differ from many other space-based weather observations. Instead of imaging, CYGNSS data are comprised of tracks of specular points that streak across the Earth surface. Images are reconstructed from composites of multiple satellites. This unique sampling mechanism poses challenges for assimilation into numerical weather prediction (NWP) models.

Data assimilation (DA) is the component of a NWP system that incorporates observations with evolving forecast models. A common data DA technique poses the problem as variational equation that minimizes a cost function [2]:

$$\arg \underset{\mathbf{x}}{\min }\mathbf{J}\left(\mathbf{x}\right)={\left(\mathbf{x}-{\mathbf{x}}_{\mathbf{b}}\right)}^T{\mathbf{B}}^{-1}\left(\mathbf{x}-{\mathbf{x}}_{\mathbf{b}}\right)+{\left(\mathbf{y}-H\left[\mathbf{x}\right]\right)}^T{\mathbf{R}}^{-1}\left(\mathbf{y}-H\left[\mathbf{x}\right]\right)$$

(1)

where $\mathbf{J}\left(\mathbf{x}\right)$ represents the cost function, $\mathbf{x}$ is the variational argument, ${\mathbf{x}}_{\mathbf{b}}$ is the background model state of dimension n, $\mathbf{B}$ is the n-by-n background error covariance matrix, $\mathbf{y}$ is the observation vector of dimension p, $H$ is the forward model operator for transforming from model space to observation space, and $\mathbf{R}$ is the p-by-p observation error covariance matrix.

In practice, a great deal of focus is placed on estimating B [3,4] because incorrectly specified B can have significant impacts to the performance of the model. Famously, introducing perfect observations to a perfect model can still lead to forecast degradations if B is poorly constructed [5].

On the other hand, $\mathbf{R}$ is comparatively overlooked. NWP data assimilation (DA) schemes often assume each observation error is completely uncorrelated [6,7,8]. This is performed for three main reasons: to simplify the assumptions needed for minimizing the cost function in Equation 1 and reduce the computational resources required to solve the forecast problem; to reduce the data handling and storage requirements for observation data at low latencies; and because frequently the information required to specify correlated observation errors was unavailable or difficult to ascertain.

There are some major recent advancements on this front. Discussions of correlated instrument errors are generally limited to explorations for interchannel correlated errors in infrared [9,10] and microwave [11] sounders, with some physical models of correlation proposed [12]. Spatial observation error correlation is discussed in [11,13] but generally restricted to estimates of total observation error correlation as estimated from a weather model [14]. More recently, instrument error inventories for microwave sounders have been established [15].

Nevertheless, the assumption that observation errors are independent and uncorrelated is not valid for any observatory. In practice, observation error correlations are minimized by thinning the dataset or ‘super-obbed’ samples, which improves overall model skill by mitigating the effects of correlated error, which may be interpreted by the model as a real signal [16,17,18].

The premise of not fully exploiting observation data is unappealing, considering that the process of acquiring these data usually requires the considerable capital expense of building and flying remote sensing satellites. The highest utility of dense observations in space and time occurs when the application (i.e., the model) resolution matches the observation resolution. For many imagers and sounders, there may be a reasonable case that thinning is an appropriate way to scale dense observation data to match model gridding.

For CYGNSS, however, thinning the data is an especially unpalatable solution, as the peak utility of CYGNSS data is during the relatively uncommon occurrence that a specular point passes through the eyewall of a tropical cyclone. Thinning may miss this observation or misrepresent the structure of the storm in the model. By a similar token, ‘super-obbing’ is not a practicable remedy, because tracks are one dimensional and collocated observations are potentially hours apart. Blending data from a wide temporal window risks misrepresenting the dynamics of the tropical cyclone and negates the fast-revisit utility of the observatory.

Estimating the correlated error structure for observations can improve model skill in NWP [6,13,19,20,21]. A number of methods have been demonstrated to estimate the observation error correlation matrix [22,23], most notably Desroziers’ diagnostic [24].

Our work does not estimate the full observation error correlation matrix that is typical of these prior works. Instead, we provide a first-principles, bottoms-up, tunable engineering model for how the CYGNSS instrumentation itself can cause observations to contain correlated errors. This correlation model is validated and tuned against empirical observation data during a period when two CYGNSS assets were in a specific orbital geometry where observations nearly coincided in space in time. We further discuss the challenges and opportunities that result from this correlated error model, as well as highlight potential for applicability to future assimilation investigations.

2. Materials and Methods

2.1. The CYGNSS Observatory

The CYGNSS main observable is normalized bistatic radar cross section (NBRCS), ${\sigma }^o$. This is derived from the radar range equation [25]:

$${\mathit{\sigma }}^{\mathit{o}}=\frac{{\mathit{P}}^{\mathit{g}}{\left(4\mathit{\pi }\right)}^3{\mathit{L}}^{\mathit{a}\mathit{t}\mathit{m}}}{{\mathit{P}}^{\mathit{T}}{\mathit{G}}^{\mathit{T}}{\mathit{\lambda }}^2{\mathit{G}}^{\mathit{R}}{\mathit{R}}^{\mathit{t}\mathit{o}\mathit{t}}\overline{\mathit{A}}}$$

(2)

where $P^g$ is the calibrated power received by the nadir science antenna; $L^{atm}$ are the losses due to attenuation in the atmosphere; $P^T$ is the power of the transmitter, which in this case is a GPS satellite; $G^T$ is the antenna gain of the GPS satellite in the direction of the Earth specular point; $\lambda$ is the wavelength of the GPS signal, which for CYGNSS is the GPS L1 signal at 1575 MHz; $G^R$ is the gain of the CYGNSS nadir science antenna in the direction of the Earth specular point; $R^{tot}$ is the total range loss of the transmission via the Earth specular point; and $\overline{A}$ is the effective scattering area of the specular point.

NBRCS is directly related to the mean square slope (MSS) of ocean wave spectra via [25,26]:

$${\mathit{\sigma }}^{\mathit{o}}\left(\mathit{\theta }\right)=\frac{{\left|\mathfrak{R}\left(\mathit{\theta }\right)\right|}^2}{\mathit{M}\mathit{S}\mathit{S}}$$

(3)

where $\theta$ is the angle of incidence of the scattering geometry and $\mathfrak{R}\left(\theta \right)$ is the Fresnel electric field reflection coefficient for the ocean surface at the specular point. CYGNSS retrieves windspeed from ${\sigma }^o$ via an empirical Geophysical Model Function (GMF) [27] that is processed operationally via a stored look-up table (LUT).

In practice, un-normalized bistatic radar cross section $\sigma$ is calculated for every delay and Doppler, the native coordinates of the delay-Doppler Mapping Instrument (DDMI). ${\sigma }^o$ is calculated from the average of $\sigma$ over a 3-by-5 bin window centered at the specular point, and divided by the effective area $\overline{A}$ corresponding to the surface area bounded by the 3-by-5 window [$m^2/m^2$].

After launch of the CYGNSS mission, it became clear that uncertainty in the GPS effective isotropic radiated power (EIRP) would significantly impact calibration quality, [28,29,30]. The EIRP is simply the product $P^T{G}^T$ from Equation 2. Errors are present in both terms: not only is there general uncertainty with the actual (versus published) transmit power of GPS satellites, newer generations of the constellation operate with a flexible power mode and dynamically change transmit power levels across their respective orbits [31]. Further, only a subset of the GPS antenna patterns were published [32] and they were not sufficiently detailed to constrain the error in EIRP for CYGNSS calibration.

To remedy this, Wang et. al. developed a novel calibration technique that utilizes the onboard zenith CYGNSS antenna and empirically derived relationship from measured GPS antenna patterns to estimate GPS EIRP in the direction of the Earth specular point at the time of observation [33]. This new technique has been adopted in the CYGNSS Level 1 calibration algorithm since version 3.0, and modifies the radar equation to become the following [34]:

$${\mathit{\sigma }}^{\mathit{o}}=\frac{{\mathit{P}}^{\mathit{g}}4\mathit{\pi }{\mathit{L}}^{\mathit{a}\mathit{t}\mathit{m}}{\mathit{L}}^{\mathit{Z}}{\mathit{G}}^{\mathit{Z}}}{{\mathit{P}}^{\mathit{Z}}{\mathit{\zeta }}^{\mathit{E}}{\mathit{G}}^{\mathit{R}}{\mathit{R}}^{\mathit{t}\mathit{o}\mathit{t}}\overline{\mathit{A}}}$$

(4)

where all the terms are the same as in Equation 2, with the additions of $L^Z$, the transmission range loss from the GPS source to the CYGNSS zenith antenna (written as a loss in the numerator to distinguish from the $R^{tot}$ range losses); the CYGNSS zenith antenna gain $G^Z$ in the direction of the GPS satellite; the power received from the CYGNSS zenith receiver $P^Z$; and the specular-zenith ratio ${\zeta }^E$ for the GPS EIRP as described in [33].

Each of the terms in Equation 4 can potentially be a source of correlated error, but for practical reasons this study only evaluates those sources with the largest error magnitude and therefore the largest potential utility for future data assimilation users. The five largest sources of error are explored as shown in Table 1, representing a 1-sigma error magnitude for a reference windspeed at 10 m/s. The following sections evaluate each of these terms in depth.

2.2. The Bottoms-up Correlated Error Model

The CYGNSS Level 1 correlated error model is constructed to represent the major error sources based on the physical and operational characteristics of the CYGNSS observatories. These errors are correlated over both space and time, and depend on several variables relating to the operation of the CYGNSS constellation.

Construction of the error model is intended to be realistic enough to represent measurable and plausible correlated errors, yet flexible enough to allow for tuning and parameterization. The general modeled structure of correlated error takes the form:

$$R_{mod}\left(i,j\right)=\frac{1}{\mathcal{N}}\sum _n{K}_n(i,j)$$

(5)

where $R_{mod}\left(i,j\right)$ is the normalized correlated error between any two samples $i$ and $j$ with a domain $-1\le R_{mod}\le 1$. There are $n$ component terms that are added together, and each of the $n$th term corresponds to a unique source of error. $K_n\left(i,j\right)$ represents the error covariance between any two samples $i$ and $j$ for a specific component $n$. $\mathcal{N}$ is a normalization constant to ensure that the diagonal terms of the $R_{mod}$ matrix is 1.

Each of the $K_n$ terms can be further generalized as:

$$K_n\left(i,j\right)={\left(E\left(n\right)\right)}^2\cdot Ecor{r}_n\left(i,j\right)$$

(6)

where $E\left(n\right)$ represents the magnitude of each error component $n$ as calculated in [33,34,36], and can be thought of as the variance of the error; and $Ecor{r}_n\left(i,j\right)$ is a derived function with a domain $-1\le Ecorr\le 1$ that represents the correlation in error between samples $i$ and $j$ for each component $n$.

Each $Ecorr$ function depends on different arguments, depending on the nature of the error source. The individual errors are explored in depth in the following sections. Further, several tuning parameters have been added to assist with validation. The initial values of the tuning parameters are 1 and essentially leave the model output unmodified. However, this model assumes that the relative magnitude of the error components.

2.2.1. Model Assumptions

Several assumptions are made regarding the overall model and each individual $Ecorr$ term. The first is that the decorrelation scales of interest are on the order of seconds to minutes. This restricts the problem to errors that decorrelate on timescales that would be of relevance to numerical weather prediction users. CYGNSS observations may have correlated error structure that evolves at longer timescales, say, daily or seasonally. However, this structure would probably be better resolved via other calibration and processing techniques and would not add significant value to a correlated error matrix. As a result, all errors between CYGNSS samples $i$ and $j$ are assumed to be 100% uncorrelated if there is more than 10 minutes of separation in observation time. This corresponds to roughly 4500 km of distance, where observations could reasonably be treated as wholly uncorrelated.

Further, all error terms $Ecor{r}_n$ are assumed to be uncorrelated with each other. This assumption simplifies the implementation and calculation of the model, but is less supported theoretically. There are a number of plausible rationales for why many of these error sources are correlated. However, this model hopes to correct for these empirically with the implementation of tuning parameters.

Finally, there is a general assumption throughout that all the error sources in the model exhibit wide-sense stationarity during the timescales of interest. Therefore, the error magnitude values $E\left(n\right)$ are treated as constants. In practice, the magnitude of errors and their correlation time and space scale sizes may well vary, e.g. seasonally and regionally. A fully operational implementation of the methodology developed here would take these dependencies into account, e.g. by an appropriate parameterization of the properties of the errors. Here, we only consider stationary error properties in order to illustrate how they are determined and how they would be combined to be used by a DA scheme.

2.2.1. The Full Correlated Error Model

Considering the terms in Table 1 and applying Equation 5, the Level 1 correlated error model is constructed:

$$R_{mod}\left(i,j\right)=\frac{1}{\mathcal{N}}\left(\begin{array}{c} R_{P^g}{\left(i,j\right)}^{\hspace{0.17em}}+\\ R_{P^Z}{\left(i,j\right)}^{\hspace{0.17em}}+\\ R_{G^R}{\left(i,j\right)}^{\hspace{0.17em}}+\\ R_{G^Z}{\left(i,j\right)}^{\hspace{0.17em}}+\\ R_{\zeta }{\left(i,j\right)}^{\hspace{0.17em}}\end{array}\right)$$

(7)

where $R_{P^g}$ is the correlated error component from the calibrated nadir receiver power and discussed in Appendix B; $R_{P^Z}$ is the correlated error component from the zenith receiver and discussed in Appendix C; $R_{G^R}$ and $R_{G^Z}$ are the correlated error components from the nadir and zenith CYGNSS antenna patterns, respectively, and discussed in Appendix D; and $R_{\zeta }$ is the correlated error component due to the zenith-specular ratio used for dynamic EIRP estimation and discussed in Appendix E.

2.3. Verification Techniques

One of the primary challenges in constructing the CYGNSS Level 1 correlated error model $R_{mod}$ is identifying a plausible validation scenario. In practice, it can be difficult to disentangle the various sources of correlated observation error structure, e.g. those caused by the inherent behavior and calibration of the instrument and, by the geophysical retrieval and inversion process, and by the representativity errors imposed when observations are gridded and ingested into models. Further, the choice of ground truth may impose an additional source of correlated error, such as when using reanalysis data or another observation source.

In an effort to disentangle these correlated errors and isolate only those due to instrumental sources, this work validates the correlated error structure $R_{mod}$ by matching up near- simultaneous collocated observations made by two different observatories at nearly identical geometries, generating modeled NBRCS ${\sigma }_{mod}^o$ from reanalysis data, and single- and double-differencing the results.

2.3.1. Curating Matchup Observations

As a constellation of eight small satellites in the same orbital plane, CYGNSS is generally unable to collect collocated, near-simultaneous observations from different satellites. If the CYGNSS observatories were equally distributed across the orbital plane, each asset would follow the next at approximately a 10 minute lag. However, in that time both of the uniquely defining features of the CYGNSS observation changes: (1) the surface of the Earth rotates underneath the observatories, changing the ground track, and (2) the GPS satellites that serve as the source of radar signal advance in their orbit. Therefore in 10 minutes, it can be quite challenging to develop 1-to-1 matchup conditions suitable for investigating the correlated error structure at timescales of seconds to minutes.

CYGNSS has no onboard thrusters, and orbit phasing is controlled solely through differental drag. At several junctures during the CYGNSS mission, one of the observatories has advanced within the plane to be nearly overlapping with another, yet at a slightly different altitude. Generally, the greater the altitude difference between the observatories, the faster the relative precession within the orbit plane.

An exhaustive review of each observatory’s ephemerides for the life of the mission identified a few matchup opportunities. Matchups near the beginning of the mission were preferred, as the orbit planes of the satellites tend to drift apart over the lifetime of the mission. After further filtering by operational status of each observatory to ensure they were in similar attitude configurations and operating in similar science modes, a 24 hour period starting 0Z on 11 SEP 2019 was chosen, where FM1 and FM5 were in a nearly identical orbital phase for a sustained period. A representation of the ground sampling during this period is shown in Figure 1.

Each “track” of observations (when a series of observations are made in close succession sharing a CYGNSS receiver and GPS transmitter) was matched sample-for-sample between FM1 and FM5 by first minimizing the distance between observations, and then quality controlling for several factors:

• Matched tracks must both have greater than 300 samples;
• Individual sample matchups are valid if samples are within 0.5 degrees (great circle distance);
• Individual samples are screened to ensure no quality control flags apply;
• The matched track is only valid if 60% of the data remains after all other matchup criteria applies.

The resulting matchup conditions produced 103 matched tracks within the 24 hour period of near-simultaneous observations.

2.3.2. Generating Model NBRCS

A forward model for the CYGNSS observatory (generating ${\mathsf{\sigma }}_{\mathrm{m}\mathrm{o}\mathrm{d}}^{\mathrm{o}}$ from wind data) can be challenging. Operationally, CYGNSS uses a GMF LUT that maps between the two quantities. However, for the purposes of this analysis, we prefer to use a physically-representative forward operator that represents the physical dynamics of the roughening of the ocean surface due to locally driven winds.

For each sample matched up, ${\mathsf{\sigma }}_{\mathrm{m}\mathrm{o}\mathrm{d}}^{\mathrm{o}}$ is generated by computing a spectrally-corrected ocean surface MSS from an ERA5-forced [37] WaveWatch III wave model as described in [38]. Because the resolution of ERA5 is much coarser than CYGNSS Level 1 observations in both space and time, the ${\mathsf{\sigma }}_{\mathrm{m}\mathrm{o}\mathrm{d}}^{\mathrm{o}}$ is matched up with ${\mathsf{\sigma }}_{\hspace{0.17em}}^{\mathrm{o}}$ via a tri-linear interpolation across three dimensions (latitude, longitude, and time).

Therefore, for every track of data, there are observations from two different CYGNSS observatories, and two modeled NBRCS from reanalysis data.

2.3.3. Estimating Total Correlated Error

Single- and double-differencing are common techniques to calibrate remote sensing instruments, especially radiometers [39,40] and radars [41]. Both are useful for quantifying the correlated error structures in CYGNSS observations.

For each track of data, there are two single-differenced datasets,

$$\mathit{S}{\mathit{D}}_{\mathit{o}\mathit{b}\mathit{s}}={\mathit{\sigma }}_{\mathbf{1}}^{\mathit{o}}-{\mathit{\sigma }}_{\mathbf{5}}^{\mathit{o}}$$

(8a)

$$\mathit{S}{\mathit{D}}_{\mathit{m}\mathit{o}\mathit{d}}={\mathit{\sigma }}_{\mathit{m}\mathit{o}\mathit{d},\mathbf{1}}^{\mathit{o}}-{\mathit{\sigma }}_{\mathit{m}\mathit{o}\mathit{d},\mathbf{5}}^{\mathit{o}}$$

(8b)

where $S{D}_{obs}$ is the single-differenced data from the two matched-up observations from FM1 and FM5, and $S{D}_{mod}$ is the single-differenced data of the model-generated NBRCS for the specific coordinates of FM1 and FM5. The double-difference is computed by differencing these two quantities:

$$\mathit{D}\mathit{D}=\mathit{S}{\mathit{D}}_{\mathit{o}\mathit{b}\mathit{s}}-\mathit{S}{\mathit{D}}_{\mathit{m}\mathit{o}\mathit{d}}$$

(9)

Both the SDs and DD are useful for our analysis. The $S{D}_{obs}$ term represents the difference in ${\sigma }^o$measured by two different CYGNSS satellites. Despite the strict matchup conditions, these assets are still measuring different fields of view at different times. These differences may result in some residual correlated structure. In addition, any correlated structure in $S{D}_{obs}$ may not reveal structure imposed from fundamental properties of the earth system. Therefore, $S{D}_{mod}$ allows us to identify the correlated structure of any systematic differences in observation target. An example of the utility of these differencing techniques is shown in Figure 2.

Since we are primarily interested in correlated error, and not absolute error, our primary investigative tool is the autocorrelation function. The correlated error for each track can be computed by autocorrelating the DD using the standard formulation:

$$\mathit{\rho }\left(\mathit{\tau }\right)=\frac{1}{(\mathit{N}-\mathit{\tau })}\frac{1}{{{\mathit{\sigma }}_{\mathit{x}}}_{\mathit{i}}}\frac{1}{{{\mathit{\sigma }}_{\mathit{x}}}_{\mathit{i}+\mathit{\tau }}}\sum _{\mathit{i}=1}^{\mathit{N}-\mathit{\tau }}{(\mathit{X}}_{\mathit{i}}-\overline{\mathit{X}}){(\mathit{X}}_{\mathit{i}+\mathit{\tau }}-\overline{\mathit{X}})$$

(10)

where $\tau$ is the lag, $N$ is the total number of lags observed, $X_i$ is the value of timeseries $X$ at time $i$, $\overline{X}$ is the mean of timeseries $X$, and the standard deviation is formulated as normal, ${\sigma }_i=\sqrt{\frac{\sum _{i=1}^{N-\tau }{(X}_i-\overline{X}{)}^2}{N-\tau }}$.

Because the error correlation behavior can vary from track-to-track, we also construct a bulk autocorrelation which approximates the autocorrelation behavior for all tracks sampled. For this we consolidate $M$ tracks, the bulk autocorrelation is:

$$\mathcal{R}\left(\mathit{\tau }\right)=\frac{1}{\mathit{M}\left(\mathit{N}-\mathit{\tau }\right)}\frac{1}{{{\mathit{\sigma }}_{\mathit{y}}}_{\mathit{i}}}\frac{1}{{{\mathit{\sigma }}_{\mathit{y}}}_{\mathit{i}+\mathit{\tau }}}\sum _{\mathit{i}=1}^{\mathit{N}-\mathit{\tau }}{(\mathit{Y}}_{\mathit{i}}-\overline{\mathit{Y}}){(\mathit{Y}}_{\mathit{i}+\mathit{\tau }}-\overline{\mathit{Y}})$$

(11)

where $Y_i=\sum _{j=1}^M{X}_{i,j}$, $j$ is the index for track number, and the standard deviation and mean are calculated as before, but with the consolidated track series $Y$.

2.3.4. Model Tuning

The constructed correlated error model $R_{mod}\left(i,j\right)$ is designed to estimate the correlated error between arbitrary samples $i$ and $j$. Validating the correlated error with empirical matches is challenging since the empirical error correlation can only be estimated in a broader, statistical sense. To create an appropriate comparison, we introduce the modeled error autocorrelation ${\tilde{\rho }}_{mod}$, which is a function of lag $\tau$:

$${\tilde{\rho }}_{mod}\left(\tau \right)=\frac{1}{N-\tau }\sum _{i=1}^{N-\tau }{R}_{mod}\left(i,i+\tau \right)$$

(12)

where $R_{mod}$ is calculated in Equation 7, $\tau$ is the lag, $N$ is the total number of lags observed. The quantity ${\tilde{\rho }}_{mod}$ is reasonably comparable to the autocorrelation $\rho \left(\tau \right)$ of observed error for a single track of samples. A similar analog can be made for bulk modeled error autocorrelation ${\tilde{\mathcal{R}}}_{mod}\left(\tau \right)$, which will estimate the autocorrelation behavior of the error model across M tracks:

$${\tilde{\mathcal{R}}}_{mod}\left(\tau \right)=\frac{1}{M}\sum _{j=1}^M{\tilde{\rho }}_{mod,j} \left(\tau \right)$$

(13)

where ${\tilde{\rho }}_{mod,j} \left(\tau \right)$ is the modeled autocorrelation for track $j$ at lag $\tau$.

A parametrized version of $R_{mod}$ is constructed of the form

$$R_{mod}\left(\alpha ,\beta ,\gamma ,\delta \right){\left.\right|}_{i,j}=\frac{1}{\mathcal{N}}\left[R_{P^g}\left(\alpha ,\beta \right)+R_{P^Z}\left(\alpha ,\beta \right)+R_{G^R}\left(\gamma ,\delta \right)+R_{G^Z}\left(\gamma ,\delta \right)+R_{\zeta }\left(\beta \right) \right]$$

(14)

where

• $\alpha$ represents the relative magnitude of the white noise component of the error, which decorrelates at $\tau =1$;
• $\beta$ represents the relative magnitude long-decay pedestal, or any residual correlated errors at the edge of our timescales of interest;
• $\gamma$ represents the relative magnitude of the correlated caused by terms $R_{G^R}$ and $R_{G^Z}$, which exhibit smooth decay as samples spread apart when projected through the nadir and zenith antenna coordinates, respectively [see Appendix D for an in-depth discussion]; and
• $\delta$ represents the relative decorrelation roll-off in terms $R_{G^R}$ and $R_{G^Z}$.

With an appropriate benchmark, the tuning parameters $\alpha$, $\beta$, $\gamma$, and $\delta$ are iterated such that ${\tilde{\mathcal{R}}}_{mod}$ matches a target signal. The tuning parameters are applied such that they can be modified to change specific characteristics of the modeled behavior:

The parameters modify $R_{mod}$ directly and are described in depth in Appendix B, Appendix C，Appendix D and Appendix E. The target for matching the modeled autocorrelation ${\tilde{\mathcal{R}}}_{mod}$ is the bulk autocorrelation of the double differenced data DD:

$$\mathcal{R}\left[DD\right]\left(\tau \right)\cong {\tilde{\mathcal{R}}}_{mod}\left(\tau \right)$$

(15)

Because this is an under-determined system, there is no unique solution to optimizing the tuning parameters. Further, these parameters all relate to one another, and changing one will impact the others. Instead, $R_{mod}$ is tuned heuristically such that the modeled behavior matches the empirical data at key points: to match the white noise component at lag $\tau =1$; to match the rolloff at lags $\tau =5$ and $\tau =30$; and to match the endpoint behavior at lag $\tau =100$.

3. Results

3.1. Bulk Behavior

The constructed model $R_{mod}$ is first tuned to match the overall behavior of $\mathcal{R}\left[DD\right]\left(\tau \right)$ as described in Section 2.3.4. With the appropriate tuning parameters, the bulk modeled autocorrelation ${\tilde{\mathcal{R}}}_{mod}$ closely resembles $\mathcal{R}\left[DD\right]\left(\tau \right)$ in most important aspects, including relative magnitude of white noise component, relative rate of decorrelation roll-off, and endpoint behavior at large lags. The final tuned ${\tilde{\mathcal{R}}}_{mod}$ is shown in Figure 3 with $\mathcal{R}\left[S{D}_{obs}\right]$, $\mathcal{R}\left[S{D}_{mod}\right]$, and $\mathcal{R}\left[DD\right]$.

The behavior in Figure 3 is worth discussing in detail. The single-differenced model generated ${\sigma }_{mod}^o$ decays at a much slower rate than the single-differenced observations and double-differenced data suggests that the errors of the slight mismatch in observation location and time are generating errors on a fundamentally different scale than the instrument errors. Further, the fact that the double-differenced decorrelation is almost identical to the observation single-differencing suggests that single-differencing near-simultaneous observations are a reasonable approximation for bulk error correlation.

Evaluating the double difference trace in Figure 3 also reveals important qualities about the CYGNSS error correlation structure. First, the uncorrelated error component accounts for roughly 5% of the total system error. Therefore, assumptions that samples be treated as independent with uncorrelated error are empirically refuted. Further the error decorrelation rolloff is swift: about 50% of the error has decorrelated within 7 seconds of observation. Using a flat-plane approximation of Earth, this corresponds to a distance of approximately 50 km, which is generally the scale of the gridding of modern global weather models, but much coarser than the resolution of state-of-the-art regional models. Finally, the endpoint behavior at large lag times suggest a small, positive correlated error (~2%) at larger timescales. The statistical properties of the lag correlation become less stable at larger lags, but the existence of a long-timescale correlation is not surprising, considering all observations from a single receiver share a common electronics and processing chain.

The fact that our tuned model $R_{mod}$ can approximate bulk error correlation structure ${\tilde{\mathcal{R}}}_{mod}$ with similar features as $\mathcal{R}\left[DD\right]$ validates the assumption that the overall instrument correlated error can be modeled from the fundamental components of the instrument observable, in our case, the radar-range equation (Equation 2). Further, that $R_{mod}$ can be generated between arbitrary samples suggests that an observation error correlation matrix R can be generated dynamically from first-principles given appropriate knowledge about the instrument and retrieval.

The value of the chosen tuning parameters is also worth investigating. The untuned model is defined by having parameters $\alpha =\beta =\gamma =\delta =1$. Figure 3 demonstrates that the modification of tuning parameters can change the overall model behavior significantly to match observed behavior. Figure 3 further suggests that the untuned $R_{mod}$ generally overestimates both the relative magnitudes of uncorrelated error (i.e., white noise, tuned by $\alpha$) and totally-correlated error (i.e., endpoint correlation, tuned by $\beta$). The chosen parameters for tuning are shown in Table 2.

The tuning parameter $\gamma$ is used to vary the relative magnitude of the near-sample correlated error roll-off. In our tuned model, this parameter remains at 1. The tuning parameter $\delta$ is used to enhance or reduce the steepness of the rolloff in correlation, caused by two observations moving farther apart in the nadir and zenith coordinate systems. If $\delta >1$, it enhances the steepness of decay. The fact that our optimized tuning maintains $\delta =1$ suggests that the filtering theory discussed in Lemma C2 is a reasonable model of decay.

3.2. Single-Track Comparisons

The tuned model can be compared to single track autocorrelations $\rho \left(\tau \right)$ of matched-sample double differences via the modeled error autocorrelation ${\tilde{\rho }}_{mod}$ as described in Equation 12. We compare three exemplar tracks in Figure 4.

The single-track autocorrelation behavior shown in Figure 4 is also revealing of the limitations of this work. The single-track error autocorrelation for double-differenced data is highly variable. This may be due to both artifacts in the data (processing, quality control, insufficient data) as well as real behavior of the observation. It is worth articulating that autocorrelation of data is generally less stable with smaller datasets and at longer lags. With our quality control parameters, we flag out significant quantities of data, which decreases the stability of autocorrelations from track to track. We further note the difficulty in exploring the dynamics of this behavior in a statistical sense: our aperture of observation where two CYGNSS assets are in a near-overlap condition is quite rare, and as the mission progresses, the orbit planes of the satellites have drifted, making further matchup scenarios less representative. These behaviors may evolve differently day-to-day or as a function of observed windspeed, but the paucity of data in this configuration makes it challenging to establish sufficient baselines to test for significance.

3.3. Dynamic Correlated Error Estimation and Impact of Tuning

The model $R_{mod}$ can produce plausible dynamics suggesting broader importance of a realistic dynamic correlated error model. Figure 5 plots two nearby tracks captured nearly simultaneously by the same receiver (but different GPS transmitters).

Because correlated error $R_{mod}$ can be generated for any arbitrary observation using the bottoms-up model, the dynamics of how instrument errors decorrelate can be explored. If the observations were treated as completely independent without any correlated error, the matrix in Figure 5a would contain non-zero elements exclusively along the main diagonal. However, we observe several structural elements. The ‘pixelation’ pattern is largely a result of the $R_{\zeta }$ term and originates from the coarse mapping of the estimated GPS antenna gain pattern as a function of scattering incidence ${\theta }_{inc}$. This phenomenon is explored in depth in Appendix D. The smooth decorrelation roll-off is from the $R_{G^R}$ and $R_{G^Z}$ components of the model, and represent the direct and reflected signals moving about the nadir and zenith the antenna patterns in CYGNSS receivers as the observation is collected along the track. We note that the model allows for cross-track correlation where the two different tracks share a CYGNSS receiver but have a different GPS transmitter. In both the untuned and tuned cases, these tracks appear to have uncorrelated cross-track error. There are residual correlations between the tracks when the scattering geometry is such that the two tracks share similar incidence angles or are near similar antenna coordinates, in addition to the shared receiver noise at lag $\tau =0$.

We can also determine $R_{mod}$ for our matchup tracks where CYGNSS observes similar track geometries with the same GPS transmitter but with two different receivers. The correlated structure for one of the tracks in the previous example (but matched by difference receivers) is shown in Figure 6.

The model $R_{mod}$ allows for correlated error between two different receivers that share a GPS transmitter. We note that observations sharing the same GPS transmitter may contain correlated error due to correlated misestimation of GPS transmit power. This may be an incomplete articulation of the full cross-receiver error correlation. For example, CYGNSS assets in similar orbital regimes may experience similar thermal and environmental conditions that cause correlated errors during our timescales of interest. We assume that correlated error due to misestimation of GPS transmit power is nearly constant for the timescales of interest.

As demonstrated in Figure 6, the cross-track correlation virtually disappears after tuning the model $R_{mod}$. Taken with the data presented in Figure 5, there appears to be negligible cross-track error correlation, either from when two tracks share a CYGNSS receiver or when two tracks share a GPS transmitter. This is a novel result, but challenging to validate in practice, as the double-differencing may eliminate any shared error correlation between two tracks.

4. Discussion and Conclusions

This work produces a first-principles estimate of correlated instrument error with results that approximate observed statistical behavior. The correlated error model $R_{mod}$ exhibits complex dynamics that are generally plausible, but difficult to validate outside of the large-scale statistical estimation discussed in Section 3.1. Nevertheless, we believe that this model presents a significant advancement in the estimation of the spatial and temporal correlation structure of instrument error for remote sensing systems, and in particular, for the unique considerations of GNSS-R measurements by CYGNSS.

In essence, this work answers a theoretical exercise for enumerating the plausible engineering reasons why two data points from the same observing constellation can share correlated sources of error. The instrument-originating sources of correlated error is likely to be a small component of the overall correlated structures of observation error. A companion work will explore how the retrieval that maps from NBRCS to windspeed produces correlated error structures in space and time.

The correlated error model $R_{mod}$ is designed to compute the estimated correlated error from arbitrary CYGNSS samples, but in practice, this may prove burdensome to implement. It is likely that a given sufficient assumptions about the stationarity of the instrument correlated error, and the overall magnitude of the instrument-originating sources of correlated error, this information can be conveyed as a static LUT.

The fact that the bulk statistical behavior between the estimated autocorrelation model averaged over all tracks is consistent with the observed double difference autocorrelation (c.f. Figure 3) is reassuring, but the disparity between the individual track model estimated autocorrelations and the single-track double difference estimations (c.f. Figure 4) suggests that either the model is insufficient at capturing the dynamics track-to-track or that the correlated error may not be as stationary as assumed.

Further, the fact that the single-differenced observations generally drive the double-differenced autocorrelation behavior suggests that Powell et. al.’s metric of simultaneity and collocation [42] generally applies for CYGNSS when satellites are in a near overlap condition. This interesting result suggests that double-differencing may not be required for statistical estimations of GNSS-R errors given a sufficiently large dataset of observations that nearly overlap.

The specific values of the tuning parameters recovered suggest that the correlated error model $R_{mod}$ overestimates both the uncorrelated and highly-correlated components of instrument error. This is useful information for future studies, because as an opportunistic measurement, GNSS-R has limited insight to the error structures in the source signals from GPS. The fact that the error correlation decays relatively rapidly provides evidence that nearby samples in space and time can add significant new information to a forecast. Further, the fact that the tuning parameters $\gamma$ and $\delta =1$ suggest that the overall structure of the error correlation decay is consistent with theory posited in Appendix D (that the primary source of error correlation is from the application of smoothing filters in the production of antenna gain patterns) and that the application of the Filtering Lemma (Lemma D2) is reasonable.