format: text

text:
# Contents

## Chapter 1 Abstract

Abstract

Several experimental results show that it is possible to extract useful phase information from reflected GPS signals over the oceans. In this work we begin the development of the theoretical background to account for these results and fully understand the phenomena involved. This information will then be used to define and carry out new experiments to evaluate the feasibility of using the phase from reflected GPS signals for altimetric purposes and the advantages of using interferometric combinations of the signals at different frequencies—the PIP concept.

We focus on the coherence properties of the signals, including the PIP interferometric combination of phases in the different frequencies. In this work we will concentrate on a static, 8 m high receiver (at least in regards to the simulations), and an infinitely removed static source. As the ocean moves, the received field will pick up a random phase. We want to understand the behavior of this phase, as the goal is to carry out altimetric measurements using phase ranging. We will also keep in mind that this random phase carries geophysical information (intuitively, the bigger the significant wave height, the larger the phase excursions).

Our simulations are based on the Fresnel integral and use simulated Gaussian oceans using the Elfouhaily et al. spectrum. The simulation tool, FRESNEL, is capable of generating time series of the reflected field at different frequencies, and then analyzing their properties. This software is written in IDL.

The most important point we need to answer is whether the signal can be tracked in the aforementioned situations. As we show, the PIP combination of the signals helps clean the signal from noise—the more correlated the signals at different frequencies the more effective is the PIP mechanism. The following questions are specifically addressed:

* What is the reflected field spectrum? This is determined by the orbital ocean motion, for a static receiver. We also discuss the moving case and show that it depends on the system gain.
* Does the winding number accumulate, or does it average to zero? Winding number refers to phase accumulation. Large phase excursions (many cycles) are seen in our simulations. Our simulations show, however, that there isn’t a preferred direction for phase winding. This means that altimetric phase measurements can be accurate.
* What are the theoretical values for the average values of the interferometric fields? How do the fields correlate across frequencies? We show that the most important factor is the relation between significant wave height and the (real or synthetic) electromagnetic wavelength.
* What is the coherence time of the signals (the reflected phase)? We show it is longest for the PIP combination. For rough seas, the correlation between the fields (and therefore interferometric coherence) in different frequencies disappears, and the coherence time goes to zero even for the interferometric combinations. In calmer ocean conditions, however, our results indicate that the interferometric combination remains coherent while the individual signals lose coherence rapidly.
* What is the structure function of the reflected phase? We see a good fit with a random walk model, with a drift rate proportional to wind speed.
* In what ways is the PIP interferometric signal different and, presumably, superior to the original ones? Coherence, a basic element for altimetric purposes.

Finally, we discuss the robustness of altimetric phase measurements after low-pass filtering the PIP combination. We show that low-pass filtered PIP data can provide more robust and accurate altimetric measurements to detect slow-varying geophysical signals in the ocean.

## Chapter 2 Introduction

This is the sum of Technical Notes PIPAER-IEEC-TN-1100 and PIPAER-IEEC-TN-2100, Theoretical Feasibility Analysis. The inputs to this WP are

1. GMV/IEEC Proposal GMVSA1123/99 (GMV)
2. Relevant literature

The outputs are

1. This technical note
2. Recommendations for experiments and post-processing procedures

And the Tasks:

1. Review the relevant documentation
2. Establish the worst and best case scenarios for the concept applicability
3. Issue recommendations for experiments and post-processing procedures
4. Prepare this technical note

In addition we analyze some issues relating to the architecture of the proposed PIP instrument.

## Chapter 3 State-of-the-Art, Review

### 3.1 Historical overview

The PIP idea can be traced back to the PARIS concept \[Martín-Neira, 1993\], where the use of GNSS signals as sources of opportunity for bistatic altimetry was first proposed. The goal of the present research is to assess the possibilities of using the phase information of the reflected signal from a PARIS GNSS system for altimetric purposes. In particular, we will study the merits of using the phase from a dual frequency GNSS system in an interferometric fashion. Other uses of reflected phase data can also be envisioned (surely the temporal behavior of phase fluctuations and other characteristics contain geophysical information about the sea surface) but we will not explicitly investigate this possibility in what follows. Rather, we will, for now, content ourselves with providing the background for discussion of these ideas in future work. This is a natural byproduct of our interest in phase altimetry, as we will see.

That the GPS signal reflected from the ocean surface can retain coherence under some conditions has been demonstrated in the past. In \[Auber et al., 1994\] we have the first report of a GPS receiver “locking” onto a sea-reflected signal. Motivated by this somewhat surprising event, tracking of reflected signals was later seeked and achieved for extended periods of time in a series of experiments with flights (over the Chesapeake Bay and the Eastern Shore of Virginia) at an altitude of up to 5500 m \[Garrison et al., 1996\]. This was accomplished using an Ashtech Z-12 off-the-shelf receiver with a LHCP nadir looking antenna, and it was observed that carrier lock was obtained. This particular receiver would not have been able to track otherwise.

On the other hand, several groups have to date successfully analyzed reflected data in order to correlate the direct and reflected signals with clean code replicas. In the experiment carried out by ESA in September 1997 two receivers were employed, operating independently of each other, one tracking the direct signal with an up-looking antenna, the other recovering the reflected signal via a down-looking LHCP 9 dB helical antenna over an 18 m bridge near Rotterdam (The Netherlands), the “Zeeland Brug”. After amplification, each of the signals was sent to a different GPS receiver (GEC Plessey builder kit 2). The receivers down-converted the signals to IF (4.309 MHz, with a bandwidth of 1.9 MHz). The IF signals were sent directly to a high-performance sampling card, which sampled the data at 6.25 MHz with 2 MHz bandwidth. They were then cross-correlated with replicas of the signal adjusted for Doppler due to satellite motion and with 20 ms duration (20 C/A code periods). Using the correlation delay data from the direct and reflected signal, this group solved for the height of the bridge over mean surface (which depended on the tide) as well as a hardware bias constant via a least squares fitting procedure. In addition, the receivers recorded standard RINEX files, including SNR. The GPS signals were deemed usable if the SNR of the reflected signal was greater than 6 to 9 dB, and if the geometry was favorable (visibility and multi-path were rather complicated in some geometrical configurations, which were then excluded). Significant Wave Height during the experiment was around 1.4 meters (equivalent to a height standard deviation $`\sigma _\zeta `$ of about 35 cm), and the altimetric performance was rather poor, as only the C/A code could be used for correlation (in general, a 1% of the chip length is assumed to be the optimal precision for one measurement—this is 3 m for C/A code). A very good description and further details can be found in \[Caparrini, 1998\].

The IEEC Earth Sciences Department has been able to reproduce the analysis of the Bridge Experiment data with similar results. The analysis tools for the generation of the GPS signals and their correlation with the received ones have been developed in Matlab.

Another aspect of the Bridge Experiment analysis carried out by ESA was the attempt to extract phase information from the data. The first attempt of phase processing was an evaluation of Doppler frequency based on zero-crossings counting \[Caparrini, 1998\]. The approach was to first multiply the signal by an appropriate clear replica of the PRN code of a satellite in view. It was then expected that a peak in the Fourier transform of the resulting signal would appear at the IF plus the Doppler frequency. The direct approach was not deemed feasible, however, due to the limited length of the signals (only 10 ms long), which did not allow for enough spectral resolution. The signal multiplied by the replica was instead first filtered with a rather narrow filter in order to eliminate all but the expected Doppler-shifted component. After that, frequency was estimated using the zero-crossings count for each 10 ms stream of data, with rather good results for both the direct and the reflected signal. This suggested that the reflected signal still contained a certain degree of coherence. A second study was then devoted to the measurement of the relative tide height (that is, the difference in the height measurements at different times through the difference in the number of zero-crossings between the direct and the reflected signals). The reasoning was that if the phase could be measured at all, it would be only up to a phase bias. If the bias is assumed to be constant over some short measurement times, then the change in the phase between those times would be related directly to the change in delay between the direct and the reflected signals. This calculation, however, could not be carried through because the data was too noisy and the intervals between data measurements too far apart in time.

### 3.2 The goal of this research

As we have seen, several experimental results show that it is possible to extract useful phase information from reflected signals. In the analysis sketched below we will begin to develop the theoretical background to account for these results. The goal of this work is to continue this development and analysis of past work to fully understand the phenomena involved. This information will then be used to define and carry out new experiments to evaluate the feasibility of using the phase from reflected GPS signals for altimetric purposes.

We will focus on the coherence of the signals, including the interferometric combination of phases in the different frequencies. In this work we will concentrate on a static, 8 m high receiver (at least in regards to the simulations), and an infinitely removed static source. As the ocean moves, the received field will pick up a random phase. We want to understand the behavior of this phase, as the goal is to carry out altimetric measurements using phase ranging. We will also keep in mind that this random phase carries geophysical information (intuitively, the bigger the significant wave height, the larger the phase excursions).

The PIP concept concept consists in combining two phases from the signals at different frequencies. Imagine a very smooth and slowly changing ocean. By combining the two phases we obtain a slowly changing interferometric phase. In some sense we would expect that this interferometric phase will be less sensitive to small ripples in the ocean and thus provide a more robust ranging tool for altimetry.

Several implementations of this concept can be envisioned: a static receiver on a coast or bridge to monitor tides or floods as well as sea state. Or deployment on boats, aircraft and spacecraft for global altimetric or oceanographic measurements. Carrier altimetric signals over the oceans, even when seen by a fast moving Low Earth Orbiter, will mostly be low frequency. A strong topographic signal over the oceans would be a slope of 1 meter every 100 km. As seen by the LEO (traveling at 7 km/s, say), this translates into 0.35 Hz in L1, a rise of 5 cycles in 100/7 seconds. In L25, this becomes 0.01 Hz.

On the other hand, the roughness and motion of the sea will induce high frequency “jitter” on the received signal. This is a nuisance for altimetry, and this is where we expect the PIP concept to bring added value, as the interferometric combination should be less sensitive to these effects. It is appropriate to recall here the old ado, “What is noise to some is signal to others.” This jitter may contain very useful information.

In the author’s mind, the most important point we need to answer is whether the signal can be tracked in the aforementioned situations. As we will see, the PIP concept will help clean the signal from noise. For the most part, we will focus in this work on the low altitude, static receiver situation. The moving cases at higher altitudes will be left for future work, as we will discuss in the last Chapter.

In this report we will address, among others, the following questions:

1. What is the reflected field spectrum?
2. Does the winding number (to be described below) accumulate, or does it average to zero?
3. What are the theoretical values for the average values of the interferometric fields? That is, how do the fields correlate across frequencies?
4. What is the coherence time of the signals (the reflected phase)?
5. What is the structure function of the reflected phase?
6. In what ways is the interferometric signal different and, presumably, superior to the original ones?

## Chapter 4 Background and Applicable Models

The Rayleigh criterion defines a surface to be smooth if $`\sigma _\zeta <\lambda /(8\mathrm{cos}\theta )`$, where $`\sigma _\zeta `$ is the surface height standard deviation, $`\theta `$ is the incidence angle with respect to the normal to the surface. A more stringent condition is provided by the Fraunhofer criterion, which is used to define the far-field distance of an antenna, $`\sigma _\zeta <\lambda /(32\mathrm{cos}\theta )`$. A surface becomes smooth under two conditions: $`\sigma _\zeta 0`$ or $`\theta 90^o`$. The effective roughness of the surface is therefore $`\sigma _\zeta \mathrm{cos}\theta `$ \[Beckmann et al., 1963\].

Several aspects must be considered to understand the scattering phenomenon. We will examine two limiting situation models. In one case, we can assume that the surface is made up of many independent mirrors (specular points) whose contributions to the field add up incoherently. The other limit is that of a surface which deviates from flatness only slightly. These two cases must be treated differently. The goal in both cases will be to understand the behavior of the received phase and modulus. Let us keep in mind in the back of our heads that, at the end, we should look at the relative phase between the two available wavelengths, which should be a statistically better behaved quantity than either of the two phases.

### 4.1 Statistics of the field

In this subsection we briefly review some of the results in \[Beckmann et al., 1963\], concerning the sum

$$U=re^{i\psi }=\underset{j=1}{\overset{n}{}}e^{i\varphi _j}.$$
(4.1)
This is relevant to understand what happens when we sum the fields in a rough reflection situation. We state the result for two illuminating conditions:

1. A uniformly distributed phase from $`\pi `$ to $`\pi `$.
2. A normal phase distribution with standard deviation $`\sigma `$.

In the first case the resulting phase has again a uniformly distributed phase, but the modulus has the so-called Rayleigh distribution $`p(r)=\frac{2r}{n}\mathrm{exp}(r^2/n)`$—and therefore has a non-zero average (the average is proportional to $`n^{1/2}`$). In the second case we find that the components of the resulting phasor are normally distributed (which is a quite general result following directly from the Central Limit Theorem) and that

$$r^2=n^2e^{\sigma ^2}+n(1e^{\sigma ^2}).$$
(4.2)
In this beautiful expression we can see the coherent and incoherent contributions to the field modulus. Beckmann summarizes:

> Outside a narrow cone (or wedge) about the direction of specular reflection, the amplitude of the field scattered by a rough surface is always Rayleigh-distributed; if the surface is very rough, and grazing incidence is excluded, the amplitude of the scattered field is Rayleigh-distributed everywhere.

The specularly scattered field is composed of a coherent component and a random, Hoyt-distributed component. When the surface is very rough, the latter becomes incoherent and the former vanishes. In fact, if the surface height distribution is normal with deviation $`\sigma _\zeta `$, then

$$r^2=n^2e^{(4\pi \sigma _\zeta \mathrm{cos}\theta /\lambda )^2}+n(1e^{(4\pi \sigma _\zeta \mathrm{cos}\theta /\lambda )^2}).$$
(4.3)
Thus, we see that the magnitude of the reflected field should depend mainly on the ratio of significant wave height to wavelength. We will obtain a related result through more sophisticated analysis below.

Some verification of these ideas via simulations can be found in \[Daout et al., 1999\], where the authors, using a 1-D bistatic scattering simulator, reproduce the mentioned aspects of the phase statistical distribution. The surface is there simulated by facets, and the Kirchhoff theory for the scattering field from each facet is used to accumulate the rays with regard to their phase. The phase histograms they show verify very clearly the ideas just discussed. In particular, that the signal is very coherent at low elevations.

### 4.2 The Fresnel-Huygens-Kirchhoff integral for the field

This is basically the scalar Kirchhoff approximation, valid for surfaces with large radii of curvature compared to wavelength. Let the incoming field be described by111Here we do a scalar treatment, think of $`U`$ as a component of the field.

$$U_o(p)=\frac{e^{ikr}}{r}.$$
(4.4)
The Fresnel integral for the scattered field is (see \[Born & Wolf, 1993, p. 380, eq. 17\])

$$U(p)=\frac{i}{4\pi }\frac{e^{ik(r+s)}}{rs}(\stackrel{}{q}\widehat{n})𝑑S.$$
(4.5)
The vector $`\stackrel{}{q}=(\stackrel{}{q}_{},q_z)`$ is the scattering vector: the vector normal to the plane that would specularly reflect the rays in the direction we are looking. This vector is a function of the incoming and outgoing unit vectors $`\stackrel{}{n}_i`$ and $`\stackrel{}{n}_s`$, $`\stackrel{}{q}=k(\widehat{n}_i\widehat{n}_s)`$. The scattering vector is related to the specular angle $`\beta `$ through $`\mathrm{cos}\beta =\widehat{z}\widehat{q}/q`$. Note, as an aside, that in the nadir case all wind direction dependence disappears (surface anisotropy). Changing wind direction is akin to performing a rotation in the surface, and the integral is invariant under rotations in the nadir case. The nadir case is special because $`R[\stackrel{}{q}(\stackrel{}{x})]=\stackrel{}{q}(R[\stackrel{}{x}])`$.

### 4.3 Geometrical Optics

This is the process that dominates in specular scattering when the frequency is large (i.e., the wavelength is small compared to the wave height and to the surface correlation length). The case of L-band scattering from the ocean is at best borderline “high-frequency”, and we are carrying out some separate studies within the framework of GNSS-OPPSCAT222ESA Contract 13461/99/NL/GD, Utilization of Scatterometry Using Sources of Opportunity. to understand the relevance of this approximation in such circumstances. To this end, we have developed the speckles software tool, written in IDL (Interactive Data Language). In Geometrical Optics the surface height standard deviation is assumed to be at least of the order of the electromagnetic wavelength. According to the analysis in \[Beckmann et al., 1963\], the resulting wave will be largely incoherent. There is an important ingredient in our situation, however, and that is because we have a Woodward Ambiguity Function (WAF) at our disposal to select a given surface patch, as we will see shortly. This means that we can filter field contributions from a rather small area.

In the Geometrical Optics approximation to the Kirchhoff theory for electromagnetic scattering, the physical picture can be understood in terms of a specular point model. That is, the field at the receiver is the superposition of the fields generated by a number of ”mirrors” (not flat mirrors, thought, rather parabolic caps) on the scattering surface which are oriented in the correct manner. The radiation from each specular point is as coherent as the incoming radiation, but there is no coherence in the phase relationship between the radiation out-coming from the different mirrors if the surface height distribution has a large range. The result is that there is power at the receiver, but that the phase of the received signal is randomly distributed and the power is simply proportional to the number of scatterers. The behavior of the resulting phase will depend on the geometry, on the number and distribution of mirrors, and also on the temporal variation of these quantities. The goal of this section is to understand this relationship and quantify it.

It is useful here to stop and give a physical picture of the specular point model—one of the two ingredients in Geometrical Optics. We will treat here the radiation coming from each mirror separately, as was mentioned above. Since the size of each mirror is small compared to the Fresnel zone, we can use the Fraunhofer approximation (plane waves all the way through) to estimate the out-coming field.

Recall that the Fresnel integral for the scattered field is

$$U(p)=\frac{i}{4\pi }\frac{e^{ik(r+s)}}{rs}(\stackrel{}{q}\widehat{n})𝑑S.$$
(4.6)
In the far field (both emitter and receiver are very far compared to the size of the scatterer, i.e, the size of the scatterer is much smaller that the Fresnel zone), and since our integral is near a specular point with $`\widehat{n}\widehat{q}`$, we find that

$$U(p)=\frac{i}{4\pi }\frac{e^{ik(r^{}+s^{})}}{r^{}s^{}}e^{i\stackrel{}{q}\stackrel{}{r}}(\stackrel{}{q}\widehat{n})𝑑S,$$
(4.7)
see Figure 4.1 for a pictorial definition of the involved vectors.

Here all we have done is to approximate $`k(r+s)k(r^{}+s^{})\stackrel{}{q}\stackrel{}{r}`$. Now, let us focus on a single specular point: imagine that there is only one specular point on the surface. This means $`\stackrel{}{q}\widehat{n}=q`$. To perform the integral let us use a coordinate system in which the $`z`$ axis is parallel to $`\stackrel{}{q}`$. In this coordinate system the tangent plane at the specular plane is therefore parallel to the $`x`$-$`y`$ plane. Hence, $`dS=\frac{d^2x}{\stackrel{}{n}\widehat{z}}=d^2x`$. We need to compute

$$I=qe^{iqz(x,y)}d^2x,$$
(4.8)
where

$$U=\frac{i}{4\pi }\frac{e^{ik(r^{}+s^{})}}{r^{}s^{}}I.$$
(4.9)
Now we use the stationary phase approximation. The idea is that as $`q`$ gets larger the integral gets contributions only very near the specular point—the contributions farther out cancel out. Then,

$`I`$ $``$ $`I_0=q_{spec}e^{iqz_{spec}}{\displaystyle e^{i\frac{q}{2}_{ij}^2z|_{spec}x_ix_j}d^2x}`$ (4.10)
$`=`$ $`q_{spec}e^{i\stackrel{}{q}\stackrel{}{r}_{spec}}{\displaystyle \frac{2\pi }{iqdet^{1/2}\left(_{ij}^2z|_{spec}\right)}}`$ (4.11)
$`=`$ $`i_{spec}e^{i\stackrel{}{q}\stackrel{}{r}_{spec}}{\displaystyle \frac{2\pi }{det^{1/2}\left(_{ij}^2z|_{spec}\right)}}.`$ (4.12)
This result, that the square root of an inverse determinant follows from a Gaussian integral is a much used result in Quantum Field Theory (see for instance, Quantum Field Theory by L. H Ryder, Cambridge Univ Press 1996 for a derivation). It can be seen to follow from extending

$$e^{\alpha x^2}𝑑x=\sqrt{\frac{\pi }{\alpha }}$$

to higher dimensions. Now, it is not hard to show that $`det^{1/2}\left(_{ij}^2z|_{spec}\right)^1=\sqrt{r_1r_2}`$—the determinant just yields the products of the radii of curvature at the specular point: the determinant is invariant under coordinate transformations of the surface, so we can use a coordinate system along the principal directions—the surface is then described as a simple parabola. The result is then immediate (the radius of curvature is the inverse of the second derivative in such coordinates and the matrix becomes diagonal in such a coordinate system, so the computation of the determinant just yields the product of the diagonal terms). Finally,

$`U_0(p)`$ $`=`$ $`{\displaystyle \frac{i}{4\pi }}{\displaystyle \frac{e^{ik(r^{}+s^{})}}{r^{}s^{}}}i_{spec}e^{i\stackrel{}{q}\stackrel{}{r}}{\displaystyle \frac{2\pi }{det^{1/2}\left(_{ij}^2z|_{spec}\right)}}`$ (4.13)
$`=`$ $`{\displaystyle \frac{_{spec}}{2}}{\displaystyle \frac{e^{ik(r_{spec}+s_{spec})}}{r_{spec}s_{spec}}}\sqrt{r_1r_2}`$ (4.14)
The coefficient $`_{spec}`$ depends on the Fresnel coefficients and on the local geometry—we will discuss this below. The cross section is given by the ratio of the resulting field squared divided by the field squared at the surface (given above) times $`4\pi r^2`$. The result is

$$\sigma =\frac{4\pi r^2|U|^2}{|U_o|^2}=|_{spec}|^2\pi r_1r_2.$$
(4.15)
For all this to work, some requirements have to be met on top of those from the Kirchhoff approximation. Basically, this is a high frequency approximation—$`q`$ has to be large. Let us look at things in 1D, along one of the radii of curvature. The integral is essentially of the form

$`J`$ $`=`$ $`{\displaystyle _{l_x}^{l_x}}𝑑xe^{iqx^2/r_x}`$ (4.16)
$`=`$ $`\sqrt{r_x/q}{\displaystyle _{l_x\sqrt{q/r_x}}^{l_x\sqrt{q/r_x}}}e^{iu^2}𝑑u.`$ (4.17)
For the integral to become $`\sqrt{\pi }`$, we need $`l_x\sqrt{q/r_x}`$ to be large—of the order of 10. The expression for the incoherent sum of all the specular points is now immediate from Equation 4.14, given a large height deviation. This expression for the field is also useful for a not-so-incoherent sum of the fields—we will return to this point below.

#### 4.3.1 Corrections to the Geometrical Optics approximation: frequency dependence

In this section we add a comment on the next order corrections to the Geometrical Optics approximation. This is not a crucial section for the work at hand and can be taken as a small aside, but we believe it will be useful in future work.

In reality the surface is not a parabola and other terms appear. These are frequency dependent corrections to the field. In order to obtain the corrections we need to expand the integrand as a Taylor series beyond the well-behaved parabolic term. The basic idea for this type of computation is that

$$I(b)=e^{ax^2+bx}𝑑x=e^{b^2/4a}\sqrt{\pi /a},$$

and that

$$\frac{d^n}{db^n}I(b)|_{b=0}=e^{ax^2+bx}x^n𝑑x.$$

From this we can compute

$$e^{ax^2}x^4𝑑x=\frac{12}{(4a)^2}\sqrt{\pi /a}.$$

The expansion is

$$I=qe^{iqz(x,y)}d^2\stackrel{}{x}=I_0+I_4+\mathrm{},.$$
(4.18)
where

$$I_4=q_{spec}e^{iqz_{spec}}\frac{iq}{4!}_{ijkl}^4z|_{spec}x^ix^jx^kx^le^{i\frac{q}{2}_{ij}^2z|_{spec}x_ix_j}d^2\stackrel{}{x}.$$
(4.19)
Now,

$`I_4`$ $`=`$ $`q_{spec}e^{iqz_{spec}}{\displaystyle \frac{iq}{4!}}_{ijkl}^4z|_{spec}{\displaystyle x^ix^jx^kx^le^{i\frac{q}{2}_{ij}^2z|_{spec}x_ix_j}d^2\stackrel{}{x}}`$ (4.21)
$`=`$ $`q_{spec}e^{iqz_{spec}}{\displaystyle \frac{iq}{4!}}_{ijkl}^4z|_{spec}+`$
$`\left(iq_{ij}^2z|_{spec}\right)^1\left(iq_{kl}^2z|_{spec}\right)^1{\displaystyle \frac{2\pi }{iqdet^{1/2}\left(_{ij}^2z|_{spec}\right)}}`$
and

$$I_4=_{spec}e^{iqz_{spec}}\frac{2\pi \sqrt{r_1r_2}}{q\mathrm{\hspace{0.17em}4}!}_{spec}e^{iqz_{spec}}_{ijkl}^4z|_{spec}\left(_{ij}^2z|_{spec}\right)^1\left(_{kl}^2z|_{spec}\right)^1$$
(4.22)
In a future implementation of speckles (see discussion below), we would like to include these second order effects. These effects can also be accounted for theoretically using, e.g., Gaussian statistics. At the end the final expression must depend only on the usual parameters. Note that the $`I_4`$ has units of length. Thus, aside from constants, these higher order corrections may be of the form $`\frac{\sigma _\zeta }{lq}`$, say. In this last expression we see that this correction is indeed inversely proportional to frequency, as it should, since Geometric Optics is a high frequency limit.

### 4.4 The WAF zone

If the signal is filtered by cross-correlation with several C/A (or P) code periods at the delay and frequencies corresponding to the specular zone, for instance, only the surface patch in the first Delay and Doppler zone will contribute (let us call it the “WAF zone” to associate as well as distinguish this concept from the concept of Fresnel zone). This is described by the equation

$$\text{SNR}(\tau ,f_c)=\frac{1}{kT_SB_D}\frac{\lambda ^2}{(4\pi )^3}\frac{P_tG_tG_r}{R_1^2R_2^2}\sigma _{rt}^0\chi ^2(\stackrel{}{\rho },t,\delta \tau ,\delta f)𝑑A.$$

which is a slight refinement of a similar equation in \[Zavorotny et al., 1999\] and which can be found in the review \[Ruffini et al., 1999\]. The important thing to keep in mind is that the support in the integrand, the WAF zone, is the intersection of four spatial zones:

1. The receiver antenna footprint.
2. The annulus zone defined by the $`\mathrm{\Lambda }^2`$ function in $`\chi ^2`$.
3. The Doppler zone defined by the $`|S|^2`$ function in $`\chi ^2`$.
4. The scattering cross section coefficient $`\sigma _0`$.

The WAF, through its support, selects a given portion of the surface from which power is measured.

To get a feeling for these numbers, the area involved in the reflection is of the order of 10 km2 for an aircraft and 1000 km2 for a Low Earth Orbiter. For $`h>>c\tau `$, where $`\tau `$ is the chip length (1 $`\mu `$s for C/A code) and $`h`$ the receiver height), the one-WAF (or one-chip here) area can be approximated by $`A_{WAF}=2\pi hc\tau /\mathrm{cos}^2\theta `$. Note that $`A_{WAF}`$ increases only linearly with altitude.

We can add another element at this point. In the above discussion, the WAF is assumed to result from choosing a given delay and Doppler in the reflected signal. However, as time passes the signal reflected by a given surface element cannot be characterized by a fixed Delay and Doppler: both will change as the situation evolves. It is clear, however, that a more sophisticated WAF algorithm can be devised that focuses on that specific surface patch. This observation opens the door to longer integration times, which can be useful for altimetric purposes.

## Chapter 5 Simulation Tools: fresnel and speckles

These are two modeling tools we have developed to understand, via a computer simulation, some of the characteristics of the reflected field. fresnel models the reflected field (just the carrier part after removing the source time-dependence) using an ocean model and direct integration of the Fresnel-Huygens-Kirchhoff integral described above. We have carried out several types of simulations: a slowly vertically moving receiver, a static receiver, a realistic ocean and a so-called chaotic ocean. We will focus here on the static realistic ocean case.

speckles has been developed with two goals in mind. We will eventually attempt to compute the reflected field using a GO method: finding first the specular points and integrating just around these to sum up the resulting field. We will just outline here some of its characteristics and initial results. We have found that this is very difficult to do (in general, comparisons with the Fresnel integral are poor). We have seen, however that the number of specular points correlates reasonably well with the reflected Fresnel field, specially in rough ocean simulated conditions. We are attempting to understand this phenomenon better, and to study its implications: GO, in some form of another, is presently the leading model used for the analysis of reflected power in GPSR. It is a simple method, and experimental results to date seem to validate it. We would like to understand its somewhat puzzling high performance, since L-band is not really high-frequency in the ocean surface case. That is, GO is not, strictly speaking applicable in the whole range of ocean conditions with an L-band instrument—the wavelength in neither large of short enough. Nonetheless, some semi-empirical modifications of this model are currently being successfully used. We would like to understand this success and, more importantly, its limitations. We will briefly mention some of this work, although it is not central to the study.

### 5.1 fresnel

This routine ( written in Interactive Data Language—IDL) integrates the field from the ocean surface below, using nadir incidence. The receiver is hovering a few meters over the water. The ocean model is Gaussian (we use a routine from B. Chapron using the Elfouhaily spectrum) \[Elfouhaily et al., 1997\], “chaotic” (random heights, moving) and a non-Gaussian simulations are planned for future work. The software simulates the reflected complex field at L1, L2 and L5 frequencies. The goals of this part of the work were to show that the simulation is possible, to compare the behavior of the field in the different cases, and to assess the potential of the PIP idea.

#### 5.1.1 Random ocean

The random ocean is generated by choosing a height standard deviation. The code assigns a random height with the desired characteristic to each 10 cm by 10 cm square of the reflecting surface. In addition, a velocity random field is chosen, and each square moves up and down with a chosen period.

#### 5.1.2 Elfouhaily spectrum

The procedure to compute the ocean state, based on the IDL routine kindly supplied by B. Chapron, is, starting from the energy spectrum, to generate a plausible random spectrum. To do this, the square root of the energy spectrum is taken, and a uniformly random phase appended to it. The appropriate time dependence is also added to this random phase—using the dispersion relation for gravity ocean waves. Then, the inverse Fourier transform is taken. Using this procedure, a moving, random ocean with the appropriate spectral characteristics is generated. Using a uniform random phase distribution in the different frequencies means that the resulting ocean is Gaussian. This can be changed, of course, but we have not attempted to analyze the subtleties of non-uniform phase distributions in this work.

Generally the sea state is characterized by its power spectrum. As the sea evolves in a random process, the spectral density of sea elevations is obtained by the Fourier transform of the autocorrelation function of the elevations $`\zeta (x,y)`$. The spectral density can also be estimated by the Fourier transform of one realization of $`\zeta (x,y)`$ \[Blackman and Tukey, 1958\]:

$$F(k_x,k_y)=|TF[\zeta (x,y)]|^2$$

If we determine the matrix $`\varphi `$ of random phases uniformly distributed between 0 and $`2\pi `$, the sea height at the point $`r=(x,y)`$ is:

$$\zeta (r)=TF^1[\sqrt{F(k_x,k_y)}e^{i\varphi }]$$

Thus the probability density functions of heights and slopes are Gaussian. For this study we use the unified spectrum of \[Elfouhaily et al., 1997\]. This is just a consequence of the Central Limit Theorem applied to sums of harmonics.

#### 5.1.3 Field calculation

The field is calculated as a complex number using the Fresnel integral expression:

$$U(p)=\frac{i}{4\pi }\frac{e^{ik(r+s)}}{rs}(\stackrel{}{q}\widehat{n})𝑑S.$$
(5.1)
To be precise, we Fraunhofer-expand the part of the integrand corresponding to $`s`$ the distance to the infinitely far transmitter. That is, we write $`s=ks^{}\stackrel{}{k}_{in}\stackrel{}{x}`$ in the phase part of the integrand, where $`\stackrel{}{x}`$ denotes the position of the scattering point. Thus, the integrand becomes,

$$U(p)=\frac{i}{4\pi }\frac{e^{iks^{}}}{s^{}}\frac{e^{ikr\stackrel{}{k}_{in}\stackrel{}{x}}}{r}(\stackrel{}{q}\widehat{n})𝑑S.$$
(5.2)
We actually throw out all the $`s^{}`$ dependence (in effect assuming an incoming plane wave, as well as on the Fresnel coefficient ) and obtain

$$U(p)=\frac{i}{4\pi }\frac{e^{ikr\stackrel{}{k}_{in}\stackrel{}{x}}}{r}(\stackrel{}{q}\widehat{n})𝑑S.$$
(5.3)
This means that the outgoing field would be a plane wave if the surface were flat, with unit modulus.

The ocean surface is divided into squares 10 cm wide (or less), and an area of up to 200x200 meters is integrated. For the most part we will concentrate on the first chip zone, with radius $`r=\sqrt{2h\tau _c}`$.

The resulting field is a complex number in this scalar treatment, $`U(p(t))=re^{i\varphi (t)}`$. In order to evaluate the phase as the receiver moves up we take

$$\delta \varphi _i=\mathrm{ln}(w(t_i)/w(t_{i1}))$$
(5.4)
with $`w(t)=e^{i\varphi (t)}`$. Summing the $`\delta \varphi _i`$ yields the overall phase. This method is related to the definition of winding number of a curve around the origin in the complex plane:

$$I=_\gamma \frac{dz}{z}=_\gamma d(\mathrm{ln}z)=\underset{i}{}\mathrm{ln}\frac{z_i}{z_{i1}}.$$
(5.5)
A full cycle, 360 degrees, is equivalent to 19 (L1), 24 (L2), 25 (L5), 86 (L12) or 568 (L25) cm. We will return to this definition shortly.

### 5.2 speckles

As mentioned above, this code has been developed to study the use of Geometric Optics for the study of ocean surface reflections in L-band. It will also be used to study specular point statistics, using both Gaussian and, eventually, non-Gaussian models. Taking a Gaussian ocean model, we determine the specular points distribution and compute the reflected field using geometrical optics. We are to compare this to straight integration of Fresnel integral. Then the analysis will be compared with a non-linear ocean model.

#### 5.2.1 Specular point determination

We consider the incident radiation as a plane wave normal to a Gaussian surface and a receiver at $`(x_r,y_r,z_r)`$. The point $`(x,y,\zeta )`$ of the surface is considered to be a specular point if the slopes follow rather clear geometrical conditions (see for instance \[D.E.Freund 1997\]):

$$\frac{d\zeta }{dx}=\frac{xx_r}{z_r\zeta }$$
(5.6)
$$\frac{d\zeta }{dy}=\frac{yy_r}{z_r\zeta }.$$
(5.7)
Our approach consists in the detection of slope variations around the expected specular slopes defined above. We compute at each point P of the surface, and for both directions $`x`$ and $`y`$, the expected specular slope S—i.e. the slope needed to reflect the incident wave to the receiver—and the slopes of the facets located just before (slope S\[before\]) and after (slope S\[after\]) the point P. The condition for P to be a speckle is that (S\[before\]-S) and (S -S\[after\]) must be of the same sign. This determination is attractive because we doesn’t need to define any tolerance on the specular slopes. We are neglecting saddle point by this approch (saddle points are not well defined in discrete spaces).

Then an analysis is made on the radius of curvature at the specular point. The Kirchhoff approximation leads to a condition on the radius of curvature at the specular facet,

$$\sqrt{q\sqrt{r_1r_2}}>10.$$
(5.8)
The product of the principal radii of curvature can be computed using the derivatives of $`\zeta (r)`$, see \[Barrick, 1968\]:

$$|r_1r_2|=\frac{\left(1+\zeta _x^2+\zeta _y^2\right)^2}{\left|\zeta _{xx}\zeta _{yy}\zeta _{xy}^2\right|}$$

The following Figures (5.2 to 5.4) present the specular point position (circles) for three different times on the moving surface. The receiver is located at 8 meters above the surface, which has a 5 m by 5 m dimension and 1 cm of resolution. The stars correspond to the specular points that are under the curvature condition. The specular distribution is greatly modified within 10 ms.

#### 5.2.2 Field computation

The idea is to compute the Fresnel field with the contribution of scatterers only. We integrate on the vicinity of each mirror. The result is then compared to the straight Fresnel computation. To date these comparisons have not been very successful. It appears that in order for the Geometric Optics approach to work at this level (not after substantial averaging) we need to impose more stringent frequency or surface scale characteristics and/or take into account higher order effects. If the main contribution is from specular points, as we suspect, the key may be to calculate the field from each scatterer to higher order than the stationary phase approximation. Work in this area is planned for the future.

## Chapter 6 Coherence and Structure Functions

### 6.1 The coherence and structure functions, and the coherence time

Another question that needs to be answered is how the sea surface motion and varying geometry affect the coherence of the reflected signal. There are several models for the sea spectrum, and these can be used to try to extract this information. Consider the toy model of one-dimensional scattering, with static, lined up, receiver and transmitter, with the scatterer (a single facet) also in line but now moving up and down with average speed $`v`$. It is easy to see that the correlation of the direct and the scattered signal is zero unless the coherent integration time is less than the characteristic time of the surface at GPS wavelengths. The coherent integration time is roughly that time for which the RMS phase error is 1 radian \[Thomson et al., 1990\]. More precisely, let the coherence function be defined by

$$c(T)=\left|\frac{1}{T}_0^Te^{i\varphi (t)}𝑑t\right|,$$
(6.1)
where $`\varphi (t)`$ is the error phase. Then, the coherence time is the time $`T`$ it takes for $`C(T)^2`$ to drop to 0.5—say \[Thomson et al., 1990\]. Note that for long times the coherence function can also be related to the magnitude of the zero-frequency spectral component of the (windowed) normalized field associated to the error phase (that is, in our simulations, the Fourier transform of the normalized field). We will get back to this point below.

The coherence time of the GPS signal is about 2–3 hours. What about the reflected signal? If the WAF zone and the sea surface are fixed, we should expect the same. If we include time dependence (and we ignore all motions), we can rewrite the equation on the statistics of the signal by

$$U(t)=re^{i\psi }=\underset{j}{}e^{i\varphi _ji\omega t}=e^{i\omega t}\underset{j}{}e^{i\varphi _j},$$
(6.2)
so the resulting field is as coherent as the incident one. But in real situations neither of this is the case. Both the ocean surface and the WAF patch are moving (although this second effect can be corrected in principle by the use of SAR techniques). The characteristic time of the surface can roughly be defined as the time it takes the surface to move a wavelength—in the present case $`\lambda /v`$. This thought experiment is useful to understand the limitation in coherent integration time for a moving surface like the ocean. An advantage of using a synthetic wavelength using two GNSS frequencies is that this coherent integration time should be longer. We will verify in a moment that this is the case.

An associated concept to the coherence function is the structure function of the phase. This is defined by

$$s(T)=\left(\varphi (t+T)\varphi (t)\right)^2_t.$$
(6.3)
It gives us a measure of how tied up are two points in the phase as we vary their temporal separation. If this phase drift can be approximated by a random walk stochastic process \[Herring, 1990\],

$$s(T)=dT$$
(6.4)
we can use the associated “Geophysically Induced Phase drift” (GIP) rate $`\sqrt{d}`$ as a geophysical parameter if it correlates well to U10 in our models. This parameters has units of cycle per square root second. We will check this below (but let us anticipate that it does).

The coherence time and the structure function are related. In fact,

$$C(T)^2=\frac{1}{T^2}_0^T_0^Te^{i(\varphi (t)\varphi (t^{}))}𝑑t𝑑t^{}.$$
(6.5)
If we assume Gaussian statistics for the phase fluctuations \[Thomson et al., 1990\] we can simplify things considerably,

$`C(T)^2`$ $`=`$ $`{\displaystyle \frac{1}{T^2}}{\displaystyle _0^T}{\displaystyle _0^T}e^{s(t^{}t)/2}𝑑t𝑑t^{}.`$ (6.6)
$`=`$ $`{\displaystyle \frac{2}{T}}{\displaystyle _0^T}\left(1{\displaystyle \frac{\tau }{T}}\right)e^{s(\tau )/2}𝑑\tau .`$ (6.7)
For a random walk process, then,

$$C(T)^2=\frac{2}{T}_0^T\left(1\frac{\tau }{T}\right)e^{d\tau /2}𝑑\tau .$$
(6.8)
See Figure 6.1 for a graph of this function for various drift rates.

Note that in this model $`C(T)^2`$ decreases to zero as time increases. This fact alerts us that there is something amiss: we see in our simulations large phase excursions (certainly much larger than a cycle), but there is coherence left at large times—see Figures 6.2, 6.3, 6.4 and 6.5. How is this possible? The Gaussian model for phase fluctuations does not really apply to our situation, however. Gaussian means that phase excursions probabilities are bell-shaped, of course. Are they? In our simulations the phase of the reflected signal gathers about a point in the complex plane—the field does seem to have a Gaussian distribution about an average point. The cumulative phase does not have to have such a distribution, however. In fact, it could conceivably indulge in arbitrarily large excursions while spending more time in that average field point: this would still result in a Gaussian field distribution. However, the probability distribution for such a phase history would not be a Gaussian distribution (centered about 0, say). It would consist of a series of Gaussian-like humps spaced one cycle apart, with the one about zero the largest. In this situation we could well have a very long coherence time (since the phase modulo 2$`\pi `$ could be a rather narrow Gaussian distribution about a point) while the structure function may show high drift. This type of behavior is evident in our simulations when plotting the field points (see Figure 6.6). Note that the interferometric signal has a stronger average field value. See Figure 6.7 for a histogram of the L25 cumulative phase illustrating the multiple-hump distribution that we alluded to.

Notice that the above discussion is relevant to understand the coherence properties that were found in the Bridge Experiment reflected signal (integration times where of 10 ms, short enough for the static approximation to work in that geometry). In longer time scales the accumulated phase will pick up a random component, but if there is an overall phase change (due for instance, to a tide) it should be distinguishable from the random part (which should not add nor subtract to the overall phase on the average).

### 6.2 Doppler spread of the reflected signal

Another thought experiment is useful. Imagine that the surface height is varying as a function of time as

$$\zeta (t)=\alpha t+Ar(t),$$

where $`\alpha `$ is a constant and the second term adds a random, noise-like displacement between $`A`$ and $`A`$. We can think of this as the effect of a tide superimposed on the wave motion on the sea. The scattered phase will not be coherent if A is large enough, but the accumulated phase will certainly contain information about $`\alpha `$. What happens if there are many scatterers, how fast is the residual phase varying? This will give us a feel for the Phase Lock Loop (PLL) bandwidth required for tracking. It is not hard to see that the bandwidth required will be less than 100 Hz (assuming that the sea surface vertical speed is not greater than 10 m/s), in a static geometry (or if a static-WAF is used). The largest frequency present in the signal cannot be greater than the contribution due to the fastest moving surface patch—the result is always the linear combination of such contributions, $`U(t)re^{i\psi (t)}=_{j=1}^ne^{i\varphi _j(t)}`$. Now, if the maximal orbital velocity for a wave is 10 m/s (see Figure 6.8, this implies a Doppler of $`2\times 10/\lambda 100`$ Hz (remember that the signal is bouncing off the surface). Thus, a sampling rate of at least 200 Hz is required to retrieve the phase of the reflected signal (even if the SNR is high enough). In our simulations we have seen that it is safe to sample the signal at such rate.

Here is another useful question: given that a receiving LEO is moving at, say, 7 km/s, would it be possible to integrate coherently (correlate) the direct and reflected signals for 1 ms, if the LEO antenna has a footprint of 10 meters? Note that this “eye on the sea” is moving very fast over the surface, and this is introducing a random, time-dependent phase to the signal. The answer is no. But if the footprint is enlarged enough, yes. This may be an important problem for the PIP, but there are ways through which it could be fixed. One obvious one is to use a WAF that sticks to a given surface patch, or that it moves slowly over the sea surface. Some of these ideas are discussed within the context of the GNSS-OPPSCAT project (Report WP3320). To get an idea of the induced Doppler, we can reason that the scattered field is just a sum of the type

$$F=R(t)e^{i\mathrm{\Phi }(t)}\underset{i=1}{\overset{N}{}}r_i(t)e^{i\varphi _i(t)},$$

which already looks like a Fourier sum if we approximate $`r_i(t)`$ as constant and the exponent by its first derivative. Thus, we can guess that the maximum frequency in the spectrum is simply Doppler due to the fastest moving patch on the ocean. As an upper bound we can use 10 m/s—see Figure 6.9. Multiplying times 2 and dividing by $`\lambda `$ we obtain a scattering Doppler of about 100 Hz.

To be a bit more precise, let us write

$$F=_{A_{WAF}}e^{iq_zz}d^2x=_{A_{WAF}}e^{iq_z(z_0+\dot{z}t)}d^2x,$$

using the simpler backscattering expression in the Fraunhofer zone. Notice the restriction of the integral to the (static) WAF zone. Let us now rewrite the transform of the field:

$`\stackrel{~}{F}(\omega )`$ $`=`$ $`{\displaystyle Fe^{i\omega t}𝑑t}={\displaystyle e^{iq_zz_o}e^{i(\omega tq_z\dot{z})}𝑑t𝑑x}`$
$``$ $`{\displaystyle _{A(q\dot{z}\omega )}}e^{iq_zz_o}d^2x,`$
if we assume that vertical motions take over a “long enough” time compared to the time scale associated to the frequency we are interested in. This is saying that the spectral component at a given frequency is proportional to the area on the surface moving at the right speed to produce that frequency.

We can extend this reasoning to the moving-receiver case. There are two limiting situations we can think about. One is the “high gain” situation. Let us first discuss this situation—we will get back to the low gain situation at the end. Assume for now the ocean is frozen. In this situation the WAF zone is very small (due to antenna or processing gain) and coherence time is definitely influenced by the fact that as the receiver moves the contributing surface is changing. If the WAF zone is smaller than the correlation length of the ocean (certainly an extreme case, but never mind), the coherence time is limited by $`l/v`$, where $`v`$ is the receiver speed and $`l`$ is the correlation length of the surface—assuming large enough height deviations. To be more precise, if we use a replica with fixed Doppler and Delay, the “eye on the sea” will move over the ocean at a speed similar to the receiver’s. This will introduce additional bandwidth in the carrier if the ocean is rough—even if is frozen. This bandwidth will be mainly proportional to ocean surface rugosity, not to ocean surface motion. It can be very large (kHz), depending on the receiver speed and the WAF area size. Following the above reasoning, and keeping in mind that the WAF zone is now a moving filter on the ocean surface,

$`F(t)`$ $`=`$ $`{\displaystyle _{A_{WAF}(t)}}e^{iq_zz(\stackrel{}{x})}d^2x={\displaystyle \mathrm{\Xi }^{WAF}(\stackrel{}{x}\stackrel{}{v}t)e^{iq_zz(\stackrel{}{x})}d^2x}`$
$`=`$ $`{\displaystyle \mathrm{\Xi }^{WAF}(\stackrel{}{x})e^{iq_zz(\stackrel{}{x}+\stackrel{}{v}t)}d^2x}`$
$`=`$ $`{\displaystyle _{A_{WAF}}}e^{iq_zz(\stackrel{}{x}+\stackrel{}{v}t)}d^2x`$
Now let us look at $`\dot{F}(t)`$:

$$\dot{F}(t)=\stackrel{}{v}_{A_{WAF}}\stackrel{}{}\left(e^{iq_zz(\stackrel{}{x}+\stackrel{}{v}t)}\right)d^2x=\stackrel{}{v}_{A_{WAF}}e^{iq_zz(\stackrel{}{x}+\stackrel{}{v}t)}\widehat{n}𝑑l$$
(6.9)
where the last is an integral over the boundary of the WAF area and $`\stackrel{}{n}`$ is the normal to the boundary. The result illustrates the “edge effect”, i.e., that the change is just due to the change at the boundary of the WAF area. We have used a 2D version of Green’s theorem in the last step. This equation says that the change in the field is proportional to the variation of the field contribution at the edges of the WAF area (properly mapped in the velocity direction). The faster the velocity and the larger the difference, the larger the rate of change of the field. To relate this to the spectrum, note that this will contribute a high frequency componene to the spectrum. Whether this is a relatively large or small contribution depends on the total field, which is proportional to the total area. Thus, for large WAF zones this high frequency effect contributes a small portion of the total spectrum.

If we allow for a moving surface, $`z=z(\stackrel{}{x},t)`$ it is readily seen that the complete result is the sum of two distinct effects:

$$\dot{F}(t)=\stackrel{}{v}_{A_{WAF}}e^{iq_zz(\stackrel{}{x}+\stackrel{}{v}t,t)}\widehat{n}𝑑l+\mathrm{\Xi }^{WAF}(\stackrel{}{x}\stackrel{}{v}t)_t\left(e^{iq_zz(\stackrel{}{x},t)}\right)d^2x.$$
(6.10)
Let us now return to

$$F(t)=_{A_{WAF}(t)}e^{iq_zz(\stackrel{}{x})}d^2x=\mathrm{\Xi }^{WAF}(\stackrel{}{x}\stackrel{}{v}t)e^{iq_zz(\stackrel{}{x})}d^2x$$
(6.11)
and compute the Fourier transform:

$$\stackrel{~}{F}(\omega )=_{\mathrm{}}^{\mathrm{}}𝑑te^{i\omega t}\mathrm{\Xi }^{WAF}(\stackrel{}{x}\stackrel{}{v}t)e^{iq_zz(\stackrel{}{x})}d^2x.$$
(6.12)
To evaluate this, let us assume for simplicity and without loss of generality that $`\mathrm{\Xi }^{WAF}`$ is a box of size $`S`$ moving along in the $`x`$ direction with velocity $`v_x`$:

$$\mathrm{\Xi }^{WAF}(\stackrel{}{x}\stackrel{}{v}t)=\mathrm{\Xi }_x(xv_xt)\mathrm{\Xi }_y(y).$$
(6.13)
Then

$$_{\mathrm{}}^{\mathrm{}}𝑑te^{i\omega t}\mathrm{\Xi }^{WAF}(\stackrel{}{x}\stackrel{}{v}t)=2\frac{\mathrm{\Xi }_y(y)e^{i\omega x/v_x}}{w}\mathrm{sin}(\omega S/v_x).$$
(6.14)
Hence,

$$\stackrel{~}{F}(\omega )=\frac{2\mathrm{sin}(\omega S/v_x)}{\omega }\mathrm{\Xi }_y(y)e^{i\omega x/v_x}e^{iq_zz(\stackrel{}{x})}d^2x.$$
(6.15)
The characteristic time is roughly given by $`S/v_x`$, as can be seen from the multiplying sinc. Frequencies higher than $`v_x/S`$ are supressed by this factor. The power in frequencies smaller than this are modulated by the horizontal roughness of $`\mathrm{exp}iq_zz(\stackrel{}{x})`$ in the scale defined by $`v_x/\omega `$ in the direction of motion: if the surface does not vary in that scale, there will be little power at the frequency $`v_x/S`$.

In this approximation we have implicitly assumed a very small WAF in relation to the geometry, since we have worked in the Fraunhofer approximation. This means that all the points in the reflecting patch have the same geometry and generate the same geometrically induced Doppler. In recent experimental conditions this is hardly the case, and the Doppler of the reflected signal is dominated by the Doppler amplitude of the WAF and glistening zones. If the glistening zone is large enough, then radiation will be received from the entire WAF. The Doppler span of the WAF zone determines then the Doppler spread of the received field. This is the basic mechanism for Delay-Doppler mapping and the basis for the relationship betweeen waveforms and the characteristics of the sea surface. Even if we restrict the signal to the first chip in the C/A code, there will be quite a bit of geometric Doppler. This second source of Doppler is in principle removable by many means: high antenna gain or processing gain. By any of these means we can reduce the patch to the Fraunhofer zone mantaining SNR. It will then become increasingly important to use SAR techniques to focus and eliminate the edge effect we just discussed. We plan to revisit these calculations extending them to the Fresnel zone in future work.

The ultimate limit to the Doppler width and therefore to the coherent integration time of the signal is ocean motion–that cannot be anticipated and compensated for.

The size of the first Fresnel zone for a receiver at 350 km is about $`2\sqrt{2h\lambda }`$ or about 1 km. In this narrow area, the geometric Doppler spans about 200 Hz edge to edge (the receiver looking down and is moving along at 7 km/s) and the corresponding integration time is 5 ms. This Doppler is of the same order of magnitude of ocean induced spread. Integration times longer than 10 ms will certainly begin to be sensitive to ocean motion. The edge effect we discussed earlier will also start to play a role at longer integration times : for a 1 km WAF (roughly corresponding to 5 ms integration time) the corresponding frequency is of 7 Hz— still too small to make spotlight processing necessary. This may be the best way to retrieve ocean induced Doppler spread. See Table-6.1 for a summary of these ideas.

We have two different new ways to extract sea-state information from the Carrier bandwidth, by playing with the eye location. One the one hand, we can use a static WAF to measure vertical velocities—this will entail using SAR techniques from moving platforms. On the other, we can measure surface roughness by letting the WAF drift over the surface at a speed of choice. This is just a matter of tuning the matched filter appropriately, and it could even be done from a static platform.

## Chapter 7 Statisitical Properties of the Reflected Fields

### 7.1 Analysis of field correlations

In this section we address the question of the correlation between the fields at different frequencies (say $`L_1`$ and $`L_2`$ and the future $`L_5`$) for different sea conditions. An approach to this problem is to rewrite the Fresnel integral using the “zone” concept. That is, we classify areas in the surface according to their distance to the receiver. We assume we are in the Fraunhofer-emitter zone, so the distance to the emitter plays no role in this discussion. The idea is then to rewrite the field as a sum of field contributions from equal delay ($`r`$)zones:

$$U_q=q_{r_{min}}^{r_{max}}e^{iqr}𝒜^{}(r)𝑑r.$$

This idea can be found in \[Berry, 1972\]. The function $`𝒜^{}(r)`$ is an “areal” density function that takes into account how much area contributes to each zone. It’s exact expression is not of immediate concern. In the case of a flat surface it just becomes one. In order to extend the integration over the whole $`r`$-axis, let us simply extend the definition of $`𝒜^{}(r)`$ to be zero outside $`r_{min}`$ and $`r_{max}`$. We then have,

$$U_q=q_{\mathrm{}}^{\mathrm{}}e^{iqr}𝒜^{}(r)𝑑r.$$

We can now read: the field is the Fourier transform of $`𝒜^{}(r)`$. And the question about the relationship of the field between frequencies becomes a question on the relationship between the Fourier components of this areal function at different frequencies.

Let us forget for now the meaning of the areal function, and just think in terms of the Fourier transform of a function with limited support, since we are looking at the Fourier transform of a function with support $`Lr_{min}r_{max}`$. The larger the support of the function, the more small-scale structure can be found in the Fourier transform—everything else the same. We are accostumed to thinking about this fact in the opposite domain: a function with a lot of small time-scale structure will have a large bandwidth—a large support in the frequncy domain.

Intuitively, the larger the support two given nearby frequency components will end up sampling the pulse at more separated regions. It is useful to think about two sine functions of slightly different frequencies, running side by side. As a rule of thumb, we expect that the correlation between these Fourier components will be sensitive to the correlation function of the (areal) function at a distance $`(\mathrm{\Delta }\lambda )L/\lambda _{mean}`$. This is because $`\mathrm{\Delta }\lambda `$ is the sampling distance difference gained per cycle, and $`L/\lambda _{mean}`$ is the number of cycles available in a region of size $`L`$. Therefore, before carrying out any calculation we can say that if the areal function is correlated at the distance dictated by this (maximal possible) frequency separation, then we will see correlation in the fields. It helps to look at the separation between the equi-delay contours in the surface of integration. If this separation is smaller than the correlation length of the surface, then high correlation between the fields is possible. There is an extra consideration, however. We have been assuming that the separation between contours is constant, but this is not the case. That is, we have been implicitly assuming that the contours are fixed on the surface, and all we have considered is the correlation between this contours. But in fact, as the ocean moves, the contours move, especially with large seas. This introduces an additional source of decorrelation, which gets worse with larger seas. The time-correlation between the fields is given by

$`U_{q_1}U_{q_2}^{}_t`$ $`=`$ $`q_1q_2{\displaystyle 𝒜_1^{}(r)𝒜_2^{}(r^{})_te^{iq_1riq_2r^{}}𝑑r𝑑r^{}}.`$
This is a difficult beast to deal with, so we will change strategy.

The correlation between the fields can also be calculated in another manner. It will be useful to approximate things entirely in the Fraunhofer zone, i.e., with both transmitter and receiver in the far field, assuming Gaussian statistics. The field in the nadir case is given, in 1-D and up to a constant, by

$$U=q_Ae^{iqz(x)}𝑑x.$$
(7.1)
Now,

$$|U|^2=q^2_A_Ae^{iq\left(z(x)z(x^{})\right)}𝑑x𝑑x^{},$$
(7.2)
and

$$|U|^2=q^2_A_Ae^{iq\left(z(x)z(x^{})\right)}𝑑x𝑑x^{}.$$
(7.3)
Now, if we assume Gaussian statistics, with

$$P(z(x)z(x^{}))=\frac{1}{2\pi \sigma ^2\sqrt{1\rho ^2(xx^{})}}\mathrm{exp}\left[\frac{z(x)^22\rho (xx^{})z(x)z(x^{})+z^2(x^{})}{2\sigma ^2(1\rho ^2(xx^{}))}\right],$$
(7.4)
we can carry out these calculations explicitly. Here $`\sigma `$ is the height standard deviation from the (zero) mean, and $`\rho `$ is the correlation function of the surface. It is healthy to keep in mind that in the gaussian case, there are the only parameters, together with the frequencies, that can appear in the final expressions. If we further assume a Gaussian correlation function with correlation length $`l`$, the answer to any question must be expressed in terms of dimensionally meaningful expressions contantaining, $`\sigma `$, $`l`$, $`\lambda _1`$ and $`\lambda _2`$.

For the case at hand, the key result is (\[Beckmann et al., 1963\], p. 190)

$$e^{iq_1z(x)+q_2z(x^{})}=\mathrm{exp}\left[\frac{1}{2}\sigma ^2(q_1^2+2\rho (xx^{})q_1q_2+q_2^2)\right].$$
(7.5)
This implies, for instance,

$$U_1U_2^{}=q_1q_2_A_A\mathrm{exp}\left[\frac{1}{2}\sigma ^2(q_1^22\rho (xx^{})q_1q_2+q_2^2)\right]𝑑x𝑑x^{}.$$
(7.6)
In the following we use a simple approximation,

$$\rho (u)=1|u|/l,|u|l,$$
(7.7)
else zero. In this simple case we model a Gaussian correlation function by a linear approximation. For simplicity we will also work in the 1-D case. The first step now is to define $`u=xx^{}`$, and $`v=x+x^{}`$. The Jacobian of this change of variables is $`1/2`$. The integral becomes (let $`A`$ be a 1-D area from $`L`$ to $`L`$)

$`U_1U_2^{}`$ $`=`$ $`q_1q_2{\displaystyle _{\sqrt{2}L}^{\sqrt{2}L}}{\displaystyle _{\sqrt{2}L+u}^{\sqrt{2}Lu}}{\displaystyle \frac{dudv}{2}}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\sigma ^2(q_1^22\rho (u)q_1q_2+q_2^2)\right]`$ (7.8)
$`=`$ $`\sqrt{2}Lq_1q_2{\displaystyle _{\sqrt{2}L}^{\sqrt{2}L}}𝑑u\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\sigma ^2(q_1^2+2\rho (u)q_1q_2+q_2^2)\right].`$ (7.9)
Now, using our simplified correlation function model (and we will henceforth assume that $`L>l`$), we find

$$U_1U_2^{}=2\sqrt{2}L\left(\frac{l}{\sigma ^2}\left(1e^{\sigma ^2q_1q_2}\right)e^{\sigma ^2\mathrm{\Delta }q^2/2}+q_1q_2(\sqrt{2}Ll)e^{\sigma ^2(q_1^2+q_2^2)/2}\right).$$
(7.10)
For L-band, $`q=2k`$ is about 60 per meter. If $`\sigma `$ is greater than 0.1 m, $`\sigma ^2q^2`$ is greater than 40, and all the exponential terms in this expression vanish, except for $`\mathrm{exp}[\sigma ^2\mathrm{\Delta }q^2/2]`$. That is, for $`\sigma >0.1`$ m,

$$U_1U_2^{}\frac{2\sqrt{2}Ll}{\sigma ^2}e^{\sigma ^2\mathrm{\Delta }q^2/2}.$$
(7.11)
It also follows that

$$C(U_1,U_2)=\frac{U_1U_2^{}}{\sqrt{U_1U_1^{}U_2U_2^{}}}e^{\sigma ^2\mathrm{\Delta }q^2/2}.$$
(7.12)
According to this result, the critical parameter in any reasonable ocean state is simply the significant wave height. Imagine for instance an ocean with a very large correlation length and large height standard deviation. This is mirror-like ocean with a global up-down displacement. It is quite clear that ih this example $`U_1U_2^{}`$ will be zero, since this product will be a number in the complex plane with a reasonable magnitude and a rather random phase. The average of such a complex number is zero.

The general result, relevant for smaller significant wave heights, is

$$C(U_1,U_2)=\frac{\frac{l}{\sigma ^2}\left(1e^{\sigma ^2q_1q_2}\right)e^{\sigma ^2\mathrm{\Delta }q^2/2}+q_1q_2(\sqrt{2}Ll)e^{\sigma ^2(q_1^2+q_2^2)/2}}{\sqrt{\left(\frac{l}{\sigma ^2}\left(1e^{\sigma ^2q_1^2}\right)+q_1^2(\sqrt{2}Ll)e^{\sigma ^2q_1^2}\right)\left(\frac{l}{\sigma ^2}\left(1e^{\sigma ^2q_2^2}\right)+q_2^2(\sqrt{2}Ll)e^{\sigma ^2q_2^2}\right)}}.$$
(7.13)
We show a plot for the case $`l=5`$ and $`L=20`$ in Figure 7.1.

Although we did not discuss it, it is rather immediate that the average field for the pure fields are governed again by $`\mathrm{exp}[\sigma ^2q^2/2]`$—this can be read off Equation 7.5 assuming a zero correlation function to infer $`\mathrm{exp}[iqz(x)]=\mathrm{exp}[\sigma ^2q^2/2]`$.

We have checked this trend in our data. In general, the interferometric combination average fields are larger. For instance, for a $`\sigma _\zeta `$= 79 cm simulation (50 seconds at 0.0025 temporal resolution, 44 m size, spatial resolution 10cm, U10=1100cmpers)—see Figure 7.4, the L1, L2, L5, L12, L25 average fields were:

> * $`(+0.03,0.01)`$
> * $`(+0.00,0.06)`$
> * $`(0.04,+0.03)`$
> * $`(0.02,+0.09)`$
> * $`(0.12,+0.07)`$

while for U10=4 ms, $`\sigma _\zeta `$= 12 cm simulation, 19 m side (see Figures 6.5 , 6.6 and 6.7 for more details about this simulation)

> * $`(+0.05,0.03)`$
> * $`(0.06,0.12)`$
> * $`(0.04,+0.10)`$
> * $`(0.09,+0.12)`$
> * $`(0.63,+0.20)`$

This is where the real strength of PIP becomes evident. Moreover, low-pass filtering (see next section) does not change these numbers much—the average field is robust.

The Fraunhofer approximation has given us a useful means to understand the correlation properties of the signals in different frequencies. We should note, however, that in the GPS case we are not in the Fraunhofer zone—we are in the WAF zone, which is is usually substantially larger. Understanding all the implications of the extension to the WAF zone involves the analysis of the “beast” in Equation 7.1 which is substantially more complicated by the factors inthe areal function. This is left for future work.

### 7.2 Filtering

Another way to understand and take advantadge of the coherence proporties is to focus on the low frequency aspects of the reflected signals for the static case. As we mentioned above, in the static situation the coherence function of the scattered phase is related to the zeroth Fourier component of the field—see Equation 6.1. Filtering the signal to retrieve slow varying geophysical signals is thus a plausible approach. As an example, we show in Figure 7.3 the fields with and without filtering at 0.5 Hz (relevant for the case of a LEO). The result is quite good, as can be seen. This behavior is also expected from the analysis above of the average field above, of course. The interferometric combination is again superior, especially L25. For the satellite of aircraft case, good modeling of the receiver position and filtering can be used to remove all but the slowly varying geophysical signals from the phase drift.

Important points: filtering doesn not seem to affect the field mean value (a very desirable result). It does significantly reduce the scatter, however. For instance, in Figure 7.2, we have the fields for a 3 m/s wind. The standard deviations of the field go from 0.87, 0.82, 0.79, 0.73, 0.68 to 0.51, 0.51, 0.49, 0.48, 0.29 (for L1, L2, L5, L12, L25) after 0.5 Hz filtering of the fields. See also the results for higher wind speed is Figures 7.3 and 7.4. The mean field remains virtually unchanged after this process.

## Chapter 8 Simulation experiments with winding number and GIP

Consider the following scenario. The receiver is static. Initially the ocean is calm, totally flat. Then a disturbance gradually appears, peaks, and slowly disappears. The ocean final state is the same as the original state: total flatness. We have run several simulations with this scenario. The result is that the winding number of the field history is non-zero. That is, the field can loop several times arond zero.

The first example is in Figure 8.1 and accompanying Figure 8.2.

How is this possible? As can be seen in the example in Figure 8.3, the sum of zero-winding-number curves can result in a winding number 1 curve. Since the resulting total field is a linear combination of such fields, all that we need to show is how the moving ocean can generate such winding number zero curves. This is easy: the only requirements are that a) the initial and final phase be the same, with the phase moving in between, and b) that the field magnitude increase or decrease before returning to its original value. It is possible to imagine examples based on GO reasoning (a usefull approximation, at any rate) to obtain this behaviour. Suppose, for instance, that the distance from a specular point to the receiver decreases and then increases, while the radius of curvature decreases then increases. It is easy to see that this yields a clockwise loop in phasor space (with the usual angular convention used in the complex plane).

We have seen the phase dispersion effect, very much related to winding number (which is a special case of the first when initial and final states are the same), in all wind conditions, with U10=2 m/s and up. At first sight this poses a significant problem in the use of phase for altimetry. How can we use a “defiting” measure for altimetry? It is now very important to characterize the statistical and geophysical properties of this phase drift (Geophysically Induced Phase drift, or GIP for short).

### 8.1 The uses of winding number

Winding number, or more generally, phase dispersion under a moving ocean is certainly related to sea state characteristics. An interesting Gedanken111Gedanken (thought) experiments are little mental exercises for exploring hypothetical experiments. The terminology comes from Einstein’s work, who is said to have engaged in them even when he was very young.: imagine a satellite tracking a reflection over a smooth region, then rough, then smooth again. The winding number induced will be related to the roughness encountered. We have already seen that, in general, the phase drifts without a specific chirality in our simulations. It remains to be seen if non-gaussian effects in the surface can lead to chirality, or preferred rotation direction. If this is the case, this phenomenon could conceivably be used for geophysical measurements. We have already seen that GIP is related to wind speed in our simulations.

## Chapter 9 Recommendations for Post-processing Procedures

Based on the work carried out so far, we are prepared to outline some of the possible incarnations of the PIP concept. What we must emphasize is that the strength lies in the noise cancellation in the interferometric combination. Filtering is best carried out after this combination.

### 9.1 The PIP instrument: architecture and data products

Based on our results, we have the following suggestions.

1. Use two wideband (200 Hz for the static receiver, possibly more for the moving case) PLL’s to extract the phase in each frequency.
2. Combine the phases and low-pass filter to recover altimetric trends. The normalized interferometric field is more coherent, so it can be used to track slow-varying altimetric changes (such as a tide). We have seen that the interferometric fiels is a more robust source, with a better defined phase.
3. Use phase drift as a geophysical parameter More generally, the phase and power time series are probably rich information source of sea state.

To be more precise, the following recipes for PIP implementation can be conceived: One, the most natural one perhaps, is the Single PLL PIP (SPIP)

1. Multiply the two signals, form the interferometric signal.
2. Filter from 0 to 0.5 Hz
3. Use a single PLL to extract the phase.

The disadvantadge is that the field magnitude will oscillate quite a bit. This may not be important, given that a filter will be used. An alternative would be to used two PLLs (MPIP)

1. Extract the phase in L1 and in L2
2. Combine the two phases
3. Filter from 0 to 0.5 Hz.

This is probably similar to the previous case, but it takes the field magnitude out of the picture from the very beginning—probably not a good thing. Finally, a non-contender is to

1. Filter from 0 to 0.5 Hz
2. Extract the phase in L1 and in L2
3. Combine the two phases

This is a bad option, it makes no use of the jitter cancellation across frequencies.

We have focused on the static case. For other situations we just need to work with excess phase and carry out the same filtering process.

## Chapter 10 Conclusions and Future Work

### 10.1 What we have found: PIP’s superior performance

1. The mean field is largest for the PIP case L25 for all sea conditions. In general the mean field is proportional to $`\mathrm{exp}[\sigma ^2q^2/2]`$, where $`q=4\pi /\lambda `$ and $`\lambda `$ is the pure or interferometric wavelength.
2. The coherence integration time is consequently substantially larger for the interferometric combination, especially L25 (if the mean sea level scatter is very large the effect disappears, of course). The interferometric phase is coherent in a reasonable range of sea conditions, unlike the single frequency phases. This implies that the PIP interferometric combination is superior to extract altimetric low frequency trends in the phase, as we discussed at the beginning.
3. The phase in our simulations behaves like a random walk, and the drift rate is directly related to sea state: in general, the larger the wind speed, the higher the drift rate. This implies that a system capable of tracking the phase of the reflected signals can provide, aside from altimetric measurements, sea state information.
4. Filtering the received field does not affect its mean value, while decreasing its scatter. It can probably be used for altimetric purposes, as it removes noise but leaves the slow-varying geophysical signals alone.
5. No long-term trends in phase drift chirality have been detected in our simulations with Gaussian ocean models based on \[Elfouhaily et al., 1997\]. We have performed simulations up to 2 minutes long. This is good news for our altimetric efforts. It seems that despite the phase drift present, and the winding number it can lead to, in practice this effect should not introduce systematic effects on altimetric determination. It is important to check that the simulations are correctly imitated by nature, of course. Among other things, non-Gaussian effects in the sea could conceivably lead to real drifting.

The single most important parameter in the correlation between the fields at different frequencies is the significant wave height. This has to be compared to the synthetic wavelength. For rough seas, the correlation bethween the fields (and therefore interferometric coherence) in different frequencies disappears, and the coherence time goes to zero even for the interferometric combinations. In calmer ocean conditions, however, our results indicate that the interferometric combination remains coherent.

### 10.2 What we’d now like to know: future work

Due to the limited scope of this study we were not able to cover all the interesting aspects of this problem. Future theoretical work includes the extension of the present results to the case in which the receiver is far away from the ground and is moving. We would also like to understand the impact of changing the ocean statistics from gaussian to non-gaussian. Another important aspect that deserves further research is to extend the analysis to a bistatic situation. So far, our simulations have only considered the monostatic case. We intend to continue this work within the scope of PARIS-$`\alpha `$:

1. Examine more realistic situations (higher, faster). This includes, in particular, understanding the issues that will arise in the aircraft and LEO scenario, and it will require larger simulations—in fact it may impossible to simulate the ocean with the required size and resolution, as the WAF zone are increases linearly with height.
2. Understand the theoretical issues involving the correlations and Doppler spread of the reflected fields without using the Fruanhofer approximation for the receiver.
3. Extend the analysis to the bistatic situation.
4. Carry out a detailed study of phase drift versus wind speed.
5. Look deeper at the advantages (which we have seen and are expected) of using even closer frequencies. This may result in recomendations for GALILEO, for instance.
6. Understand the impact of non-Gaussian ocean effects.
7. Understand better the possible sources of chirality in GID. This is important to assure altimetric accuracy.
8. Obtain and analyze data for the static situation.
9. Understand the effectiveness of GO or its higher order corrections.

It is rather clear that experimental work is needed in this field. The most important experiment for the concept at this point, is a redo of the bridge experiment: a static receiver over the moving ocean, using both L1 and L2. The desired product from such an experiment from the point of view of the present work, would time-series’ of the fields. MEATEX campaigns will also provide aircraft data which will be very useful to understand some of these issues.

## Chapter 11 Acknowlegments

This is a document produced for ESA under ESTEC Contract No. 14071/99/NL/MM. The author’s are grateful to ESA for funding and permission to freely distribute this document, and are especially grateful to Manuel Martin-Neira of ESTEC, the technical officer in charge of this contract, for very valuable comments and suggestions during the course of this research. .