# Generation and propagation of an open vortex beam

Ruoyu Zeng  1  Qi Zhao  1,2  Yuanjie Yang  1  1School of Physics, University
of Electronic Science and Technology of China, Chengdu 611731, China
2Southwest Institute of Technical Physics, Chengdu 610041, China
<EMAIL_ADDRESS>
###### Abstract

We theoretically and experimentally studied a novel class of vortex beams
named open vortex beams (OVBs). Such beams are generated using Gaussian beams
diffracted by partially blocked fork-shaped gratings (PB-FSGs).The analytical
model of OVBs in the near field and far field is given by superpositions of
Hypergeometric (HyG) modes. Unlike conventional integer and fractional vortex
beams, the OVBs can have both an open ring structure and an integer
topological charge (TC). The TC is decided by the circumference covered by the
open ring. It is also quantitatively shown that a $\pi/2$ rotation of the open
ring occurs in the propagation of an OVB due to the Gouy phase shift.
Futhermore, we demonstrate experimental generation and detection of OVBs. Our
experimental results are in very good agreement with the theory. We believe
that the OVB can be the potential candidate for numerous applications, such as
particle manipulation, quantum information and optical metrology.

††journal: osajournal††articletype: Research Article

##  1 Introduction

Optical vortices (OVs) are structured laser beams with a helical phase front
$\exp(il\varphi)$, where $\varphi$ is the azimuthal angle about the optical
axis $z$ and $l$, known as the topological charge (TC), is associated to a
well-defined orbital angular momentum (OAM) $l\hbar$ per photon[1], with
$\hbar$ denoting the reduced Plank constant. Extensive studies on
generation[2, 3], detection [4, 5] and manipulation[6, 7] of OVs gave rise to
a wide variety of emerging applications: optical communication[8], optical
imaging[9], high-resolution microscopy[10], optical metrology[11] and quantum
information[12]. In the past 20 years, the diversity and complexity of the OV
family have been extraordinarily expanded. It was found that OVs are not
limited to axially propagating singularities, but can exist and evolve in many
exotic forms, such as vortex knots[13], time-varying vortices[14], transverse
vortices[15], and ray-like trajectories[16].

OVs can also be classified from many different perspectives. A vortex beam
(VB) may possess an integer or fractional TC [17]. Therefore, VBs are divided
into two categories:integer and fractional VBs. An integer VB has an annular
intensity distribution with a dark central node produced by the phase
singularity, while a fractional VB is devoid of rotational symmetry due to a
distortion at the discontinuous phase. It is able to decompose a fractional VB
into an orthogonal set of integer modes [18]. Further considering the radial
intensity distribution, a series of different VBs were obtained theoretically
and experimentally, including Laguerre-Gaussian beams [1], Bessel-Gaussian
beams [19], Mathieu beams[20], perfect vortex beams [21], anomalous vortex
beams [22], etc. Meanwhile, in addition to the above mentioned beams with a
single vortex, VBs which contain multiple singularities may also provide novel
types of OAM. VBs embedded with specific vortex arrays were generated and
modulated, such as Hermite-Laguerre-Gaussian beams [23], helical-Ince
beams[24], SU(2) geometric modes[25], and high-power laser modes in a nearly
hemispherical cavity[26]. With such a wide range of choices on VBs,
researchers can implement related applications under distinct circumstances.

To date, most studies on VBs have focused on elucidation of beams with
cylindrical symmetry. However, when such symmetric structures are broken,
light beams can present many exceptional characteristics and reveal the nature
of basic physical laws as well. For instance, Franke-Arnold et al. used a
Gaussian beam blocked by a sector aperture to manifest the uncertainty
relation between the angular position and the OAM[27]. Similarly, Volyar et
al. showed the informational entropy made by changes in the number of the OAM
states in the azimuthally perturbed VBs [28]. The relevant non-trivial facts
indicate the necessity to thoroughly investigate the fundamental properties of
the VBs subjected to the sectorial perturbations. Besides, the truncated beam
with an open ring structure was employed several times to study the features
and effects of VBs. In [29], the Gouy phase shift on a VB was observed via a
defect in the intensity profile. In [30] and [31], the behaviors of the OVs
generated by blocked optical elements were mentioned. In [32], the focal and
optical trapping of a radially polarized VB with broken axial symmetry were
realized. In [33], the OAM spectrum of the truncated VB was used to measure
object parameters. All these explorations were mostly devoted to the
phenomenons derived from the truncated or perturbed VBs, with the basic model
of such beams remaining to be explicated.

In this article, we theoretically and experimentally study the VB with an open
ring structure, namely, open vortex beam (OVB). The OVB is generated via a
partially blocked fork-shaped grating (PB-FSG). Using the classical
diffraction theory, the OVB’s propagation in the Fresnel and Fraunhofer zone
can be described by the superposition of Hypergeometric (HyG) modes. Some
important issues are clarified with this analytical model:(1) The TC of an OVB
is determined by the range of the open ring, which is controllable by
adjusting the obstacle of the grating. In particular cases, the TC can have an
integer value. (2) The far-field solution explains the $\pi/2$ rotation of the
OVB in the propagation. We also conducted experimental generation and
observation of the OVB. A Shack-Hartmann sensor was applied to record the
wavefront of the OVB. The experimental results firmly validate the theoretical
predictions. We take that the unique properties of the OVB enable its
potential applications in optical trapping, quantum information, and optical
metrology.

##  2 Theoretical model of OVB

We first derive the general model of the OVB based on the diffractive optics.
Considering a fork-shaped grating which is blocked within a certain angular
range, namely, the PB-FSG, the transmission function in the polar coordinates
$(r,\varphi)$ is expressed as

$T(r,\varphi)=\sum_{m=-\infty}^{+\infty}t_{m}\exp\left(im\left(l\varphi-\beta
r\cos\varphi\right)\right)\times\text{rect}(\frac{2\varphi-\varphi_{0}}{2\varphi_{0}})$
(1)

where $t_{m}$ is the transmission coefficient, $\beta=2\pi/D$ with $D$ being
the period of the grating, $l$ is the singularity of the grating, rect is the
rectangle function, and $\varphi_{0}$ denotes the angular range that light can
pass through. The TC of the truncated helical phase term
$\exp(iml\varphi)\text{rect}[{(2\varphi-\varphi_{0})}/{2\varphi_{0}}]$ is
calculated using results from:

$\text{TC}=ml\times\frac{\varphi_{0}}{2\pi}$ (2)

Figure 1:  Schematic diagram of the OVB generated using a Gaussian beam
diffracted by the PB-FSG.

Then, we decompose this term into an angular Fourier series. Thus, the
transmission function of the PB-FSG is rewritten as

$T(r,\varphi)=\sum_{m=-\infty}^{+\infty}t_{m}\exp\left(-im\beta
r\cos\varphi\right)\sum_{n=-\infty}^{+\infty}a_{mn}\exp\left(in\varphi\right)$
(3)

where the coefficients $a_{mn}$ are

$a_{mn}=\frac{exp(\frac{i(ml-n)\varphi_{0}}{2})\sin(\frac{i(ml-n)\varphi_{0}}{2})}{\pi(ml-n)}$
(4)

Now using a Gaussian beam as the incident beam of the grating, the electric
field at $z=0$ is written as

$E(r,\varphi,0)=\exp(-\frac{r^{2}}{w^{2}})\times T(r,\varphi)$ (5)

where $w$ is the waist width of the Gaussian beam. When this beam illuminates
the PB-FSG, an OVB arises from the +1 order of the diffraction pattern, as
shown in Fig. 1. With the paraxial propagation of the light field in Eq. (5),
its complex amplitude at distance $z$ can be decided by the Fresnel
transformation, which takes the following form:

$\begin{split}&E(\rho,\theta,z)=\frac{-ik}{2\pi
z}\exp(ikz)exp(i\frac{k}{2z}\rho^{2})\\\
&\times\int_{0}^{2\pi}E(r,\varphi,0)\exp\biggl{\\{}\frac{ik}{2z}\big{[}r^{2}-2r\rho\cos(\varphi-\theta)\big{]}\biggl{\\}}rdrd\varphi\end{split}$
(6)

where $k$ is the wave number and $(\rho,\theta)$ are the polar coordinates in
a plane offset by distance $z$ from the initial plane. By substituting Eq.(3)
into Eq.(5) and Eq.(5) into Eq.(6), one can separate an arbitrary order
diffracted beam from the diffraction integral. The $m$-th order diffracted
beam is given by:

$\begin{split}&E_{m}(\rho,\theta,z)=\frac{-ik}{2\pi
z}\exp(ikz)exp(i\frac{k}{2z}\rho^{2})\times\int_{0}^{+\infty}\int_{0}^{2\pi}\sum_{n=-\infty}^{+\infty}a_{mn}\exp\left(-\frac{r^{2}}{w^{2}}\right)\\\
&\times\exp\left(in\varphi-im\beta
r\cos\varphi\right)\times\exp\biggl{\\{}\frac{ik}{2z}\big{[}r^{2}-2r\rho\cos(\varphi-\theta)\big{]}\biggl{\\}}rdrd\varphi\end{split}$
(7)

One can obtain the following formula by solving Eq. (7):

$\begin{split}&E_{m}(\rho_{m},\theta_{m},z)=\frac{-ik}{2\pi
z}\exp(ikz)\exp(i\frac{k}{2z}\rho^{2})\\\
&\times\sum_{n=-\infty}^{+\infty}a_{mn}\frac{\rho_{m}^{\lvert
n\rvert}}{2^{\lvert n\rvert+1}\sqrt{(w^{-2}-ik/2z)^{\lvert
n\rvert+2}}}\frac{\Gamma(\frac{1}{2}\lvert n\rvert+1)}{\Gamma(\lvert
n\rvert+1)}\\\ &{{}_{1}}{F}{{}_{1}}\biggl{[}\frac{1}{2}\lvert n\rvert+1,\lvert
n\rvert+1;-{\rho_{m}^{2}}/{(4w^{-2}-2ik/z)}\biggr{]}\exp\bigl{[}in(\theta_{m})\bigr{]}\end{split}$
(8)

where $\Gamma(x)$ is the Gamma function, ${{}_{1}}F{{}_{1}}(a,b;x)$ is a
confluent hypergeometric function, and $(\rho_{m},\theta_{m})$ are transformed
from $(\rho,\theta)$:

$\begin{cases}&\rho_{m}=\sqrt{m^{2}\beta^{2}+k^{2}\rho^{2}/z^{2}+2m\beta
k\rho\cos\theta/z}\\\ &\theta_{m}=\mp\tan^{-1}\bigl{[}\cot\theta+m\beta
z/(k\rho\sin\theta)\bigr{]}\\\ \end{cases}$ (9)

Particularly, Eq.(7), in the case that the diffracted beam is +1 order,
describes the OVB in the near-field regime. It can be obviously seen that the
OVB is reprensented as the superposition of HyG vortex modes. VBs depicted by
a similar mathematical formula are found in Ref.[34] and [35]. We also notice
that the sign of $\theta_{m}$ are inverse to that of $ml$, which implies the
orientation of the energy flow in the OVB.

In the far-field regime, the diffraction integral is written as:

$\begin{split}&E(\rho,\theta,z)=\frac{-ik}{2\pi
z}\exp(ikz)exp(i\frac{k}{2z}\rho^{2})\\\
&\times\int_{0}^{+\infty}\int_{0}^{2\pi}E(r,\varphi,0)\times\exp\big{[}-i\frac{kr\rho}{z}\cos(\varphi-\theta)\big{]}rdrd\varphi\end{split}$
(10)

By the fact that the light field in the focal plane of a Fourier lens is
equivalent to that in the Fraunhofer zone, we can employ an alternative form
of Eq. (12):

$\begin{split}&E(\rho,\theta,z)=\frac{-ik}{2\pi
z}\exp(ikz)exp(i\frac{k}{2z}\rho^{2})\\\
&\times\sum_{m=-\infty}^{+\infty}\biggl{\\{}t_{m}\mathscr{F}\bigl{(}\exp\left(-im\beta
r\cos\varphi\right)\bigr{)}*\sum_{n=-\infty}^{+\infty}a_{mn}\mathscr{F}\bigl{[}{\exp\left(in\varphi\right)\times
exp(-\frac{r^{2}}{w^{2}})}\bigr{]}\biggl{\\}}\end{split}$ (11)

where $\mathscr{F}$ is the two-dimensional Fourier transformation (FT) and $*$
is the convolution. The first FT characterizes the position of the $m$-th
order diffracted beam:

$\begin{split}\mathscr{F}\bigl{(}\exp\left(-im\beta
r\cos\varphi\right)\bigr{)}=2\pi\times\delta(\rho\cos{\theta}+m\beta)\delta(\rho\sin{\theta})\end{split}$
(12)

where $\delta(x)$ is the Dirac delta function. The second FT gives the
intensity and phase distribution of the $m$-th order diffracted beam:

$\begin{split}&\mathscr{F}\bigl{(}{\exp\left(in\varphi\right)\times
exp(-\frac{r^{2}}{w^{2}})}\bigr{)}=\pi\bigl{(}\frac{k\rho}{2z}\bigr{)}^{\lvert
n\rvert}w^{\lvert n\rvert+2}\frac{\Gamma(\frac{1}{2}\lvert
n\rvert+1)}{\Gamma(\lvert n\rvert+1)}\\\
&{{}_{1}}{F}{{}_{1}}(\frac{1}{2}\lvert n\rvert+1,\lvert
n\rvert+1;-\frac{k^{2}\rho^{2}w^{2}}{4z^{2}})\exp\bigl{(}in(\theta\mp\frac{\pi}{2})\bigr{)}\end{split}$
(13)

Therefore, the OVB in the far field is expressed as

$\begin{split}&E_{+1}(\rho,\theta,z)=\frac{-ik}{2\pi
z}\exp(ikz)\exp(i\frac{k}{2z}\rho^{2})\\\ &\times\sum_{n=-\infty}^{+\infty}\pi
a_{1n}\bigl{(}\frac{k\rho}{2z}\bigr{)}^{\lvert n\rvert}w^{\lvert
n\rvert+2}\frac{\Gamma(\frac{1}{2}\lvert n\rvert+1)}{\Gamma(\lvert
n\rvert+1)}\\\ &{{}_{1}}{F}{{}_{1}}(\frac{1}{2}\lvert n\rvert+1;\lvert
n\rvert+1;-\frac{k^{2}\rho^{2}w^{2}}{4z^{2}})\exp\bigl{(}in(\theta\mp\frac{\pi}{2})\bigr{)}\end{split}$
(14)

Compared with the initial phase, each HyG mode of Eq.(16) adds a $\pi/2$ phase
rotation, the direction of which depends on whether the energy flow of the OVB
circulates in the clockwise or counterclockwise, namely, the sign of the TC.
The open ring structure makes this rotation visible in the OVB’s propagation.
Visualizations of this effect are shown in the upcoming discussion.

This model gives a precise description of the OVB, and also lays the
foundation for studying the OVB by numerical simulations, of which we extend
details in the next section.

##  3 Simulation results

The desired PB-FSGs defined by Eq. (1) are shown in Fig. 2 (a)-(c). The
essential parameters are $\varphi_{0}=\pi/2\;\&\;\pi\;\&\;3\pi/2$ and $l=4$.
The far-field distributions of the OVBs generated by these gratings can be
approximately fitted using Eq. (16), as shown in Fig. 2 (d)-(f). The OVBs have
a quarter ring, a half ring and a three-quarter ring respectively. The
corresponding phase profiles are also illustrated in Fig. 3 (g)-(i). As can be
seen, in the original angle range that light cannot be transmitted, the phase
temporarily changes and then falls back, which does not contribute to the
total TC. However, the local phase variation in the open ring is no different
from those of any conventional VBs with closed ring structures. Thus, the TC
in Eq. (3) is conserved in the propagation of an OVB. In addition, if the
circumference of the open ring is appropriate, the phase increment in a cycle
can be a multiple of $2\pi$. The "cut-off" of the phase profiles in (g)-(i)
result to TCs of 1,2 and 3. This indicates that an OVB can possess an integer
TC.

Figure 2:  (a)-(c): PB-FSGs with the parameters: $l=4$;(a)
$\varphi_{0}=\pi/2$, (b) $\varphi_{0}=\pi$, (c) $\varphi_{0}=3\pi/2$. (d)-(f):
Corresponding far-field intensity distributions. (g)-(i): Corresponding far-
field phase distributions.

More examples of simulated OVBs are demonstrated in Fig. 3. The applied
parameters are as follows: $l=3-9$;
(a):$\varphi_{0}=\pi$,(b)$\varphi_{0}=3\pi/2$. Within our expectations, the
radius of the OVB increases with the singularity index $l$ of the PB-FSG.

Figure 3:  Simulated OVBs generated by PB-FSGs with parameters: $l=3-9$;
(a)$\varphi_{0}=3\pi/2$, (b)$\varphi_{0}=\pi$.

Another vital effect of the OVB is its rotation. A fast Fresnel algorithm [36]
can be employed to simulate the OVB’s propagation after it passing through a
thin lens (focal length $f=5$m). In Fig. 4, two OVBs with opposite-sign TCs
propagate along the optical axis $z$ and rotate in the vicinity of the Fourier
plane at $z=5$m. The three-quarter ring with TC=+3 rotates clockwise, while
the ring with TC=-3 rotates counterclockwise. This indicates that the
simulation results are consistent with the theoretical model.

The simulation results roughly visualize propagation characteristics of the
OVB. The open ring structure and tunable integer TC of the OVB can be similar
to but different from the fractional VB’s role in inducing rotation of
particles [37]. The rotation near the Fourier plane may also provide novel
optical manipulations.

Figure 4:  The OVBs rotate near the Fourier plane: (a) TC=+3,
$\varphi_{0}=3\pi/2$ (b) TC=-3, $\varphi_{0}=\pi$.

##  4 Experimental results

To validate theoretical studies presented above, we also carried out an
experimental corroboration. Normally, phase gratings can be encoded with
computer generated holograms (CFGs). In our experiments we used PB-FSGs, which
we generated with an spatial light modulator (SLM). A sketch of the
experimental setup is shown in Fig.5. A He-Ne laser beam expanded and
collimated by a beam expander (BE), illuminates a reflective SLM (Holoeye
PLUTO VIS), where the PB-FSG is displayed. The diffraction pattern of the PB-
FSG consists of multiple order diffracted beams. After a 4$f$ system with an
iris at the focal plane ($L_{1}-L_{1},f=50$mm) the first-order diffracted
light is subsequently directed to a Fourier lens (${L_{2}},f=75$mm) and
finally imaged onto a charge-coupled device (CCD). Polarizers placed ahead and
behind the SLM are used to filter the polarization and attenuate the energy of
the beam, respectively.

Figure 5: Schematic of the experimental set-up. NA, neutral attenuators; BE,
beam expander; AS, aperture slot; SLM, spatial light modulator; CCD, charge-
coupled device; $L_{1}-L_{1}$, $4f$-system.  Figure 6: OVBs generated by PB-
FSGs with parameters: $l=3-9$; (a)$\varphi_{0}=3\pi/2$, (b)$\varphi_{0}=\pi$.

The experimental results of the intensity profiles recorded on the CCD are
shown in Fig. 6, which corresponds to the simulations shown in Fig. 3 in the
previous section. We can see from the comparison that the experimental results
are in high concordance with the fitting results based on the theory.

Moreover, we observed fundamental properties of the OVB, as illustrated in
Fig. 6. Firstly, the $\pi/2$ rotation of the OVB was probed by adjusting the
position of the CCD within the focal length of $L_{2}$.

Figure 7:  (a): The experimental observation of the rotation of an OVB.
(b)-(c): wavefronts of OVBs generated using specific grating parameters: (b)
$l=7,\varphi_{0}=\pi$, (c) $l=7,\varphi_{0}=3\pi/2$. (d)-(e): the relationship
between the TC and the singularity index of the PB-FSG with
(d)$\varphi_{0}=3\pi/2$, (e)$\varphi_{0}=\pi$.

This rotation was not obvious when the beam initially passed ${L_{2}}$, but
could be clearly recognized near the focal plane, as outlined in Fig.6(a). The
open ring was ultimately rotated by 90 degrees compared with the input, which
agrees well with the theory.

Secondly, to detect the TC of the OVB, a Shack-Hartmann wavefront sensor was
used to record the phase profile of the beam. Fig. 6 (b)-(c) give two examples
of measured phase distributions with specific grating parameters: (b)
$l=7,\varphi_{0}=\pi$, (c) $l=7,\varphi_{0}=3\pi/2$. The TCs of different OVBs
were calculated via phase circulations. Wih such a scheme, the relationship
between the TC and the parameters of the PB-FSG were obtained, as shown in
Fig. 6 (d)-(e). We can see that the experimental results are in good agreement
with the theory, which indicates that the OVB can possess both an integer TC
and an open ring structure.

##  5 Conclusions

In conclusion, here we reported the OVB with novel properties. The analytical
model of the OVB is given by superpositions of HyG modes. Some unique effects
such as $\pi/2$ rotations and "cut" TCs can be validated using this model and
numerical simulations. Experimental results are in high concordance with the
theoretical predictions. Our work provides a theoretical foundation for the
previously studied VBs with open rings as well as future researchs in similar
conditions. The open ring structures and the phase circulations of OVBs are
completely different from those of integer and fractional OVs, which we expect
to be useful in particle manipulation, optical metrology, and quantum
information. We are going to extend this study by using optical forces of OVBs
to guide particles, or measuring angles of objects by OVBs.

## Funding

This work is supported by The National Natural Science Foundation of China
(11874102), Sichuan Province Science and Technology Support Program
(2020JDRC0006), and Fundamental Research Funds for the Central Universities
under grant no. ZYGX2019J102.

## Disclosures

The authors declare no conflicts of interest.

## References

  * [1] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of laguerre-gaussian laser modes,” Phys. Rev. A 45, 8185–8189 (1992).
  * [2] N. R. Heckenberg, R. McDuff, C. P. Smith, and A. G. White, “Generation of optical phase singularities by computer-generated holograms,” Opt. Lett. 17, 221–223 (1992).
  * [3] V. V. Kotlyar, A. A. Almazov, S. N. Khonina, V. A. Soifer, H. Elfstrom, and J. Turunen, “Generation of phase singularity through diffracting a plane or gaussian beam by a spiral phase plate,” J. Opt. Soc. Am. A 22, 849–861 (2005).
  * [4] G. C. G. Berkhout and M. W. Beijersbergen, “Method for probing the orbital angular momentum of optical vortices in electromagnetic waves from astronomical objects,” Phys. Rev. Lett. 101, 100801 (2008).
  * [5] Q. Zhao, M. Dong, Y. Bai, and Y. Yang, “Measuring high orbital angular momentum of vortex beams with an improved multipoint interferometer,” Photon. Res. 8, 745–749 (2020).
  * [6] Y. Yang, Q. Zhao, L. Liu, Y. Liu, C. Rosales-Guzmán, and C.-w. Qiu, “Manipulation of orbital-angular-momentum spectrum using pinhole plates,” Phys. Rev. Applied 12, 064007 (2019).
  * [7] Y. Wen, I. Chremmos, Y. Chen, Y. Zhang, and S. Yu, “Arbitrary multiplication and division of the orbital angular momentum of light,” Phys. Rev. Lett. 124, 213901 (2020).
  * [8] J. Wang, “Advances in communications using optical vortices,” Photon. Res. 4, B14–B28 (2016).
  * [9] L. Torner, J. P. Torres, and S. Carrasco, “Digital spiral imaging,” Opt. Express 13, 873–881 (2005).
  * [10] P. Bouchal, L. Štrbková, Z. Dostál, and Z. Bouchal, “Vortex topographic microscopy for full-field reference-free imaging and testing,” Opt. Express 25, 21428–21443 (2017).
  * [11] V. D’Ambrosio, N. Spagnolo, L. Del Re, S. Slussarenko, Y. Li, L. C. Kwek, L. Marrucci, S. P. Walborn, L. Aolita, and F. Sciarrino, “Photonic polarization gears for ultra-sensitive angular measurements,” Nature Communications 4, 2432 (2013).
  * [12] A. Mair, A. Vaziri, G. Weihs, and A. Zeilinger, “Entanglement of the orbital angular momentum states of photons,” Nature 412, 313–316 (2001).
  * [13] J. Leach, M. R. Dennis, J. Courtial, and M. J. Padgett, “Vortex knots in light,” New Journal of Physics 7, 55–55 (2005).
  * [14] L. Rego, K. M. Dorney, N. J. Brooks, Q. L. Nguyen, C.-T. Liao, J. San Román, D. E. Couch, A. Liu, E. Pisanty, M. Lewenstein, L. Plaja, H. C. Kapteyn, M. M. Murnane, and C. Hernández-García, “Generation of extreme-ultraviolet beams with time-varying orbital angular momentum,” Science 364 (2019).
  * [15] A. Chong, C. Wan, J. Chen, and Q. Zhan, “Generation of spatiotemporal optical vortices with controllable transverse orbital angular momentum,” Nature Photonics 14, 350–354 (2020).
  * [16] Y. Shen, X. Yang, D. Naidoo, X. Fu, and A. Forbes, “Structured ray-wave vector vortex beams in multiple degrees of freedom from a laser,” Optica 7, 820–831 (2020).
  * [17] J. Leach, E. Yao, and M. J. Padgett, “Observation of the vortex structure of a non-integer vortex beam,” New Journal of Physics 6, 71–71 (2004).
  * [18] D. Deng, M. Lin, Y. Li, and H. Zhao, “Precision measurement of fractional orbital angular momentum,” Phys. Rev. Applied 12, 014048 (2019).
  * [19] F. Gori, G. Guattari, and C. Padovani, “Bessel-gauss beams,” Optics Communications 64, 491 – 495 (1987).
  * [20] C. Lóxpez-Mariscal, J. C. Gutiérrez-Vega, G. Milne, and K. Dholakia, “Orbital angular momentum transfer in helical mathieu beams,” Opt. Express 14, 4182–4187 (2006).
  * [21] P. Vaity and L. Rusch, “Perfect vortex beam: Fourier transformation of a bessel beam,” Opt. Lett. 40, 597–600 (2015).
  * [22] Y. Yang, Y. Dong, C. Zhao, and Y. Cai, “Generation and propagation of an anomalous vortex beam,” Opt. Lett. 38, 5418–5421 (2013).
  * [23] E. G. Abramochkin and V. G. Volostnikov, “Generalized hermite-laguerre-gauss beams,” Physics of Wave Phenomena 18, 14–22 (2010).
  * [24] J. B. Bentley, J. A. Davis, M. A. Bandres, and J. C. Gutiérrez-Vega, “Generation of helical ince-gaussian beams with a liquid-crystal display,” Opt. Lett. 31, 649–651 (2006).
  * [25] Y. Shen, Z. Wan, X. Fu, Q. Liu, and M. Gong, “Vortex lattices with transverse-mode-locking states switching in a large-aperture off-axis-pumped solid-state laser,” J. Opt. Soc. Am. B 35, 2940–2944 (2018).
  * [26] Y. F. Chen, H. T. Ke, Y. C. Tseng, M. X. Hsieh, J. C. Tung, Y. H. Hsieh, H. C. Liang, and K. F. Huang, “High-power structured laser modes: direct generation of a vortex array,” Opt. Lett. 45, 4096–4099 (2020).
  * [27] S. Franke-Arnold, S. M. Barnett, E. Yao, J. Leach, J. Courtial, and M. Padgett, “Uncertainty principle for angular position and angular momentum,” New Journal of Physics 6, 103–103 (2004).
  * [28] A. Volyar, M. Bretsko, Y. Akimova, and Y. Egorov, “Orbital angular momentum and informational entropy in perturbed vortex beams,” Opt. Lett. 44, 5687–5690 (2019).
  * [29] J. Hamazaki, Y. Mineta, K. Oka, and R. Morita, “Direct observation of gouy phase shift in a propagating optical vortex,” Opt. Express 14, 8382–8392 (2006).
  * [30] D. M. Cottrell, J. A. Davis, and T. J. Hernandez, “Fraunhofer diffraction of a partially blocked spiral phase plate,” Opt. Express 19, 12873–12878 (2011).
  * [31] J. A. Davis and J. B. Bentley, “Azimuthal prism effect with partially blocked vortex-producing lenses,” Opt. Lett. 30, 3204–3206 (2005).
  * [32] Z. Man, L. Du, Y. Zhang, C. Min, S. Fu, and X. Yuan, “Focal and optical trapping behaviors of radially polarized vortex beam with broken axial symmetry,” AIP Advances 7, 065109 (2017).
  * [33] G. Xie, H. Song, Z. Zhao, G. Milione, Y. Ren, C. Liu, R. Zhang, C. Bao, L. Li, Z. Wang, K. Pang, D. Starodubov, B. Lynn, M. Tur, and A. E. Willner, “Using a complex optical orbital-angular-momentum spectrum to measure object parameters,” Opt. Lett. 42, 4482–4485 (2017).
  * [34] E. Karimi, G. Zito, B. Piccirillo, L. Marrucci, and E. Santamato, “Hypergeometric-gaussian modes,” Opt. Lett. 32, 3053–3055 (2007).
  * [35] V. V. Kotlyar, R. V. Skidanov, S. N. Khonina, and V. A. Soifer, “Hypergeometric modes,” Opt. Lett. 32, 742–744 (2007).
  * [36] J. A. Davis and D. M. Cottrell, “Ray matrix analysis of the fast fresnel transform with applications towards liquid crystal displays,” Appl. Opt. 51, 644–650 (2012).
  * [37] S. H. Tao, X.-C. Yuan, J. Lin, X. Peng, and H. B. Niu, “Fractional optical vortex beam induced rotation of particles,” Opt. Express 13, 7726–7731 (2005).


# Comparison of CYGNSS and Jason-3 Wind Speed Measurements via Gaussian
Processes

William Bekerman and Joseph Guinness This work was supported in part by the
National Science Foundation, Division of Mathematical Statistics under grant
numbers 1916208 and 1953088. (Corresponding author: Joe Guinness).The authors
are with Department of Statistics and Data Science, Cornell University,
Ithaca, NY 14853 USA (e-mail: guinness@cornell.edu).

###### Abstract

Wind is a critical component of the Earth system and has unmistakable impacts
on everyday life. The CYGNSS satellite mission improves observational coverage
of ocean winds via a fleet of eight micro-satellites that use reflected GNSS
signals to infer surface wind speed. We present analyses characterizing
variability in wind speed measurements among the eight CYGNSS satellites and
between antennas. In particular, we use a carefully constructed Gaussian
process model that leverages comparisons between CYGNSS and Jason-3 during a
one-year period from September 2019 to September 2020. The CYGNSS sensors
exhibit a range of biases, most of them between -1.0 m/s and +0.2 m/s with
respect to Jason-3, indicating that some CYGNSS sensors are biased with
respect to one another and with respect to Jason-3. The biases between the
starboard and port antennas within a CYGNSS satellite are smaller. Our results
are consistent with, yet sharper than, a more traditional paired comparison
analysis. We also explore the possibility that the bias depends on wind speed,
finding some evidence that CYGNSS satellites have positive biases with respect
to Jason-3 at low wind speeds. However, we argue that there are subtle issues
associated with estimating wind speed-dependent biases, so additional careful
statistical modeling and analysis is warranted.

##  I Introduction

Wind is a crucial component of the atmosphere and climate, having significant
implications in numerous areas of daily life, from safety and transportation
to industry and science. Recording accurate wind measurements provides
critical information for precisely defining weather hazards [1], building
skyscrapers, and landing aircraft [2]. Reliable wind speed observations are
also consequential in allowing us to efficiently conduct crop spraying [3],
monitor the global climate [4], and avoid the obstruction of essential global
shipping routes [5].

Scientists and engineers have developed a diverse suite of tools for measuring
wind speeds. Weather stations and buoys are often equipped with anemometers to
measure wind speed directly. Many earth-observing satellites carry sensors
capable of inferring wind speeds. A thorough review of satellite-based methods
for measuring wind speeds is provided in [6], which includes radiometers,
scatterometers, and altimeters, which are usually attached to low-earth-
orbiting satellites. In addition, geostationary satellites are capable of
inferring upper-air wind speeds by detecting movements in clouds via derived
motion winds algorithms [7].

The Cyclone Global Navigation Satellite System (CYGNSS) is a fleet of eight
micro-satellites that use the scattered signals from existing GNSS satellites
to infer wind speeds at the ocean surface [8, 9]. CYGNSS is a relatively new
and low-cost system that, due to its ability to distribute its sensing effort
over eight satellites, has the advantage of greater spatial-temporal coverage
of the oceans relative to a single-satellite system, an important feature for
its mission of monitoring tropical cyclones.

Our primary goal is to study the internal variability in wind speed
measurements among the eight CYGNSS satellites and across each satellite’s
starboard and port antennas. As noted in [10], this variability is still under
study, and further calibrations are a possibility: “it is more likely the
differences in bias, both between antennas and between spacecraft, are caused
by residual errors in the engineering calibration, which is performed
individually for each spacecraft and antenna. This is an ongoing area of
investigation by the CYGNSS project team . . .” Our secondary goal is to study
differences between CYGNSS and Jason-3 wind speed measurements.

Several recent articles study statistical properties of CYGNSS wind speed
measurements. CYGNSS was compared with weather model forecast winds and found
to be positively biased at low wind speeds and negatively biased at high wind
speeds, but no comparisons were made among individual CYGNSS satellites [11].
In the tropics, a similar pattern of positive bias for low wind speeds and
negative bias for high wind speeds, relative to hourly-averaged buoy data, was
detected [10]. CYGNSS biases with respect to modeled and sensed wind speeds
were assessed in the context of building a complex bias correction algorithm
[12].

The spatial-temporal patterns of low-earth-orbiting satellite observations,
coupled with the inherent spatial-temporal variability in wind speeds, present
a data-analytic challenge for conducting the desired comparisons in our study.
By design, the eight CYGNSS satellites do not measure winds concurrently at
the same locations, so it is difficult to draw direct comparisons between
observations from most pairs of these satellites. Thus, we rely on repeated
approximate crossings between CYGNSS and Jason-3 for indirect comparisons.
Jason-3 uses reflected signals from its radar altimeter to infer wind speed
and other ocean surface parameters. We believe that Jason-3 is a suitable
comparison because, like CYGNSS, it retrieves a snapshot of surface wind
speed, as opposed to an hourly average or a model forecast, both of which can
be overly smooth. To maximize statistical power, it is important to make
judicious use of the observations arising from the limited number of nearby
CYGNSS-Jason-3 passes, while avoiding a subjective determination of which
passes constitute a close-enough match.

To address these issues, we analyze the data using carefully constructed
Gaussian process models. Gaussian process models have become an indispensable
tool for analyzing and interpolating scattered remote sensing data. Recent
examples include the analysis of Argo float data [13], Orbiting Carbon
Observatory-2 data [14, 15], surface temperatures [16], Microwave Atmospheric
Satellite data [17], and Jason-3 wind speeds [18].

Our models for CYGNSS and Jason-3 wind speeds contain bias and variance
parameters that are directly related to our study goals. In particular, each
model has parameters that are interpreted as the expected difference between
CYGNSS starboard- and port-measurements and Jason-3 measurements if they had
measured wind speed at the exact same location and time. These parameters are
estimated via maximum likelihood and generalized least squares, which uses a
variance-minimizing linear combination of observations to make efficient use
of the available data. Computational challenges often associated with Gaussian
process models are overcome by downsampling across time and using a state-of-
the-art Gaussian process approximation implemented in the publicly available
GpGp R package [19]. The Gaussian process model and associated computational
techniques are the main methodological novelties of this work.

We find that there are significant and persistent differences among some pairs
of the CYGNSS sensors of a magnitude up to 1.11 m/s. There are smaller
differences between the starboard and port sensors from the same satellite.
Five of the eight CYGNSS satellites have a negative bias with respect to
Jason-3 measurements. No substantial differences in variances among the eight
satellites were detected. These results are successfully validated against a
traditional empirical analysis, which shows similar trends but higher
uncertainty. In addition, we investigate the possibility that bias depends on
wind speed, finding some evidence that CYGNSS measurements are larger than
Jason-3 at low wind speeds, though we argue that more careful analysis is
needed. We conclude the paper with a discussion of the results and suggestions
for how to modify the models to study variation in the biases. All of the code
necessary for reproducing our results is available in a GitHub repository at
https://github.com/WillBekerman/satellite-wind-speeds.

##  II Datasets and Data Processing

We compile one year of measurements recorded by CYGNSS and Jason-3 between
September 28, 2019 and September 25, 2020, specifically, CYGNSS Level 2
Science Data Record Version 3.0 and Jason-3 Level-2 X-GDR Data. Due to missing
records in the Jason-3 data during the weeks of February 1, 2020 to February
14, 2020 and June 13, 2020 to June 19, 2020, we omit these periods from our
data collection, yielding 49 total weeks of satellite measurements for our
analysis.

After acquiring the data in NetCDF format, we process the data in R, retaining
information about spatial location, time of measurement, wind speed, CYGNSS
satellite number, and CYGNSS sensor (port vs. starboard). We omit any
observations with missing data and standardize the times to seconds since
2020-01-01 00:00 UTC. Since CYGNSS does not record any measurements over land,
we keep Jason-3 observations only over open oceans and semi-enclosed seas. The
data are saved in standard R Data format.

In Figure 1, we compare the wind speed measurements taken over the same
latitudes by CYGNSS 1, CYGNSS 4, and Jason-3 during the week of November 30 to
December 6, 2019. The value in each pixel is the average of all measurements
taken within the pixel over the week. While the spatial patterns of wind
speeds are similar among the three satellites, there are subtle differences.
CYGNSS 4 appears to record larger wind speeds than CYGNSS 1, and Jason-3 tends
to have less-smooth wind fields with more of the largest values.

Figure 1: Wind speed (m/s) measurements recorded by CYGNSS 1, CYGNSS 4, and
Jason-3 between November 30, 2019 and December 6, 2019. The value in each
pixel is the sample average of all measurements falling within the pixel
during the week.  Figure 2: Distances (km) between each pair of CYGNSS micro-
satellites and between Jason-3 and CYGNSS on September 29, 2019. Gray lines
represent distances between pairs of CYGNSS micro-satellites while colored
lines depict distances between Jason-3 and CYGNSS. One black line is shown to
represent the distance between CYGNSS 1 and CYGNSS 4.

The pixel-wise comparisons in Figure 1 can be misleading because even though
the satellites have complete coverage in a one-week period, the timing of the
measurements will differ among satellites, so differences could be
attributable to the interaction between locally changing wind conditions and
the observational times, rather than bias. In seeking to quantify these
differences between satellite measurements, the simplest approach is to
directly compare measurements recorded by CYGNSS and Jason-3 taken within
small space-time windows. To explore the feasibility of the space-time
matching analysis, we plot the time-varying distances between the eight CYGNSS
satellites and Jason-3 over the course of one full day in Figure 2. All of the
CYGNSS satellites come reasonably close to Jason-3 at one or more points
during the day, with some variability in the number and proximity of such
occurrences. By contrast, while some pairs of CYGNSS satellites nearly always
measure winds at nearby locations, some pairs never do, like CYGNSS 1 and
CYGNSS 4. In the next section, we propose a model designed to address the
challenge of making efficient use of the CYGNSS and Jason-3 comparisons.

##  III Analysis

###  III-A Model Description

We first describe the statistical model used in our analysis in mathematical
notation, and then provide interpretations for the model and its statistical
parameters. Our analysis strategy is to fit separate models to datasets
consisting of Jason-3 data and one CYGNSS satellite, and repeat the analysis
for each week and each CYGNSS satellite, which means we will have many
separate fits of the same general model formulation. This allows us to study
whether biases are consistent over time. In the interest of keeping the number
of symbols manageable, we do not provide notation for each individual model
fit; the notation below is a model for an arbitrary CYGNSS satellite in an
arbitrary week. Our results will contain a re-estimation of the parameters for
each dataset.

We label the $n$ observations from both CYGNSS and Jason-3 with
$i=1,\ldots,n$. We use $Y_{i}$ for the wind speed associated with observation
$i$ and model it as follows:

$\displaystyle Y_{i}$
$\displaystyle=\mu_{i}+a_{k(i)}+Z(x_{i},t_{i})+\varepsilon_{i}$ (1)
$\displaystyle\mu_{i}$
$\displaystyle=b_{0}+b_{1}t_{i}+b_{2}(\mbox{lat})_{i}+b_{3}(\mbox{lat})^{2}_{i}+b_{4}(\mbox{lat})^{3}_{i}$
(2) $\displaystyle Z$ $\displaystyle\sim GP(0,K)$ (3)
$\displaystyle\varepsilon_{1},\ldots,\varepsilon_{n}$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}N(0,\sigma^{2}).$ (4)

The term $\mu_{i}$ is intended to capture broad-scale variation in wind speed
over time $t_{i}$ and latitude, independent of satellite. The mapping $k(i)$
indicates which sensor produced observation $i$, with $k=1$ indicating
Jason-3, $k=2$ indicating CYGNSS starboard, and $k=3$ indicating CYGNSS port.
We model space-time variation in wind speeds with the Gaussian process (GP)
$Z$, with inputs spatial location $x_{i}$ and time $t_{i}$, assumed to have
mean zero and space-time Matérn covariance function
$\mbox{Cov}(Z(x_{i},t_{i}),Z(x_{j},t_{j}))=K((x_{i},t_{i}),(x_{j},t_{j}))$,

$\displaystyle K((x_{i},t_{i}),(x_{j},t_{j}))$
$\displaystyle=\mbox{Mat\'{e}rn}(d_{ij};\theta_{1},\theta_{2})$ (5)
$\displaystyle d_{ij}$
$\displaystyle=\Bigg{(}\frac{\|x_{i}-x_{j}\|^{2}}{\theta_{3}^{2}}+\frac{(t_{i}-t_{j})^{2}}{\theta_{4}^{2}}\Bigg{)}^{1/2},$
(6)

where $\theta_{1}$ controls the variance of $Z$, $\theta_{2}$ is the Matérn
smoothness parameter, $\theta_{3}$ controls the spatial-decay of the
covariances, and $\theta_{4}$ the temporal decay.

Since $Z$ is a mean-zero process that depends only on space-time location, and
not the satellite or random independent error, we interpret it as a wind speed
anomaly, with the caveat that the anomaly is relative to a specific linear-in-
time, cubic-in-latitude mean field. The mean field contains an intercept
$b_{0}$, which means that $a_{1}$, $a_{2}$, and $a_{3}$ are not separately
identifiable, but differences such as $a_{2}-a_{1}$ and $a_{3}$ \- $a_{1}$
are. As a consequence, without additional outside information, our analysis is
not able to determine whether CYGNSS or Jason-3 is biased with respect to the
true wind field, but it is capable of assessing whether CYGNSS and Jason-3 are
biased with respect to one another via estimation of differences such as
$a_{2}-a_{1}$.

Figure 3: Model-based estimates of CYGNSS starboard (S) vs Jason-3 bias
($a_{2}-a_{1}$) and CYGNSS port (P) vs Jason-3 bias ($a_{3}-a_{1}$) for each
CYGNSS satellite and each of the 49 weeks in our study. Starboard and port
biases from the same week are connected with a black line. Discussed in more
detail in Section III-C.

In terms of the model, we are principally interested in how CYGNSS starboard,
CYGNSS port, and Jason-3 measurements would differ if they had measured wind
speed at the same location at the same time. To see how this quantity relates
to our model parameters, suppose that measurement $i$ was taken by CYGNSS
starboard and measurement $j$ was from Jason-3, and those two measurements
were recorded at the same time and location. Then

$\displaystyle Y_{i}-Y_{j}=a_{2}-a_{1}+\varepsilon_{i}-\varepsilon_{j}\sim
N(a_{2}-a_{1},2\sigma^{2}),$ (7)

meaning that $a_{2}-a_{1}$ is the bias and $2\sigma^{2}$ is the variance of
the difference. These are the parameters of interest in our study. By
comparing the estimates of $a_{2}-a_{1}$ and $a_{3}-a_{1}$ across CYGNSS
satellites, we can achieve our primary goal of understanding variability among
the CYGNSS measurements.

###  III-B Model Estimation

We fit our model separately to each of the 49 weeks and each of the eight
CYGNSS satellites, producing a total of 392 model fits. To reduce
computational burden, we fit each model using $20{,}000$ randomly selected
observations from CYGNSS and $20{,}000$ from Jason-3, and we employ a popular
computationally efficient Gaussian process approximation proposed by [20],
implemented in the R package GpGp [19]. Each model fit delivers estimates of
the model parameters via maximization of the approximate likelihood function,
as well as standard errors for the mean parameters.

###  III-C Model-Based Results

In Figure 3, for each week and each CYGNSS satellite, we plot the CYGNSS vs.
Jason-3 bias estimates $a_{2}-a_{1}$ (starboard) and $a_{3}-a_{1}$ (port).
Estimates for starboard and port from the same week and CYGNSS satellite are
connected by a black line. The bias estimates vary by week, satellite, and
antenna. Five of the eight CYGNSS satellites (1, 2, 3, 5, and 6) produce
negative biases with respect to Jason-3 for every week and for both antennas.
The three other satellites (4, 7, and 8) have a mix of negative and positive
biases across weeks. CYGNSS 1 has the largest negative biases, on average
approximately -0.83 m/s for starboard and -0.94 m/s for the port antenna.
Nearly all of the CYGNSS 1 biases are more negative than the most negative
CYGNSS 4 bias, suggesting that CYGNSS 1 and CYGNSS 4 were persistently biased
with respect to one another during our study period. Within each satellite and
antenna, the bias estimates vary by roughly 0.5 to 1.0 m/s from week to week.
The differences across weeks in the biases could be due to uncertainty in the
parameter estimates, rather than a bias that truly varies over time.

By inspecting the black lines in Figure 3, we observe differences between the
estimates of the starboard and port biases from the same satellite within the
same week. Figure 4 explores these differences in more detail by directly
plotting estimates of $a_{3}-a_{2}$, which measure the bias between the
starboard and port antennas. Most of the starboard vs. port bias estimates are
smaller than the CYGNSS vs. Jason-3 biases, with most magnitudes less than
0.25 m/s, and show a mix of both negative and positive biases, though there
are notably more negative than positive biases. CYGNSS 5 is the most lopsided
with 46 of the 49 biases being negative.

Figure 5 displays estimates of $\sqrt{2}\sigma$, which we recall from Equation
(7) is the model’s standard deviation of the difference between two
observations taken by different sensors at the same location and time. The
estimates generally fall between 0.6 and 0.9 m/s, meaning that the size of the
noise is roughly equal to the largest CYGNSS vs. Jason-3 biases. There is some
variation of the estimates of $\sqrt{2}\sigma$ across weeks but no substantial
differences among the eight CYGNSS satellites.

Figure 4: Estimates of $a_{3}-a_{2}$ for each week and each CYGNSS satellite.
Within each satellite, estimates are sorted and spaced horizontally to
visually depict the empirical quantile function of the estimates. Figure 5:
Estimates of $\sqrt{2}\sigma$ for each week and each CYGNSS satellite. Within
each satellite, estimates are sorted and spaced horizontally to visually
depict the empirical quantile function of the estimates.

###  III-D Empirical Explorations

To validate our model-based results, we conduct additional analyses based on
simple averages of differences between CYGNSS and Jason-3 wind speeds that
fall within small space-time windows. For each of the eight CYGNSS satellites,
each antenna, and each week between September 28, 2019 and September 25, 2020,
we divide the week into two-hour windows and find the pair of CYGNSS and
Jason-3 observations that are closest in distance within the two-hour window,
ignoring any windows that do not have a pair that fall within 25km. We then
take the average of the differences between the selected pairs of CYGNSS and
Jason-3 wind speeds. We refer to the averages of these differences as our
empirical bias estimates.

Analogously to Figure 3, we plot in Figure 6 the empirical bias estimates over
all CYGNSS satellites, antennas, and weeks. The general patterns of the
empirical and model-based bias estimates are quite similar. The same five
CYGNSS satellites have largely negative biases with respect to Jason-3, while
the other three have a mix of negative and positive biases. The average size
of the biases are similar as well, ranging roughly from -1.0 to 0.05 m/s. As
in the model-based analysis, the port biases in CYGNSS 5 are, on average, more
negative than the starboard biases. The empirical biases differ in that the
variation across weeks is larger than in the model-based biases, which
produces more overlap among the eight CYGNSS satellites and two antennas. For
instance, whereas there was essentially no overlap between the CYGNSS 1 and
CYGNSS 4 model-based biases, the empirical biases show more substantial
overlap. In addition, the difference between the CYGNSS 5 starboard and port
model-based biases is more clear than the difference between the empirical
biases.

Figure 6: Empirical estimates of CYGNSS starboard (S) vs Jason-3 bias and
CYGNSS port (P) vs Jason-3 bias for each CYGNSS satellite and each of the 49
weeks in our study. Starboard and port biases from the same week are connected
with a black line. Discussed in more detail in Section III-D.

To further compare the model-based and empirical bias estimates, we average
the estimates over the 49 weeks and plot them in Figure 7. The estimates
generally follow the 45 degree line, with the model-based estimates being
slightly more positive than the empirical estimates. The points tend to
cluster by satellite, suggesting that the difference between starboard and
port within a CYGNSS satellite is generally smaller than the differences among
the eight CYGNSS satellites. Interestingly, every starboard point is northeast
of its corresponding port point, indicating that the average empirical and
model-based starboard biases are more positive than the corresponding port
biases for every CYGNSS satellite.

Figure 7: Average CYGNSS satellite bias estimates from our model-based
approach and our empirical strategy.

Previous studies have explored whether bias depends on the magnitude of the
wind speed. These analyses are complicated by the fact that we never have
access to the “true” wind speed. One could take the more accurate measurement
as the “true” wind speed and estimate bias as a function of the more accurate
measurement. This approach is not without its drawbacks; when one measurement
is high, the other is likely to be lower, due to standard regression-to-the-
mean. To partially circumvent this issue, we plot in Figure 8 the difference
between CYGNSS and Jason-3 (CYGNSS minus Jason-3) against their average. As
before, we break the year into two-hour intervals and within each interval, we
extract the closest pair of observations, provided that the closest distance
is less than 25km. We see that among the port antennas, CYGNSS measurements
are usually larger than Jason-3 for small average wind speed, but for larger
average wind speeds, Jason-3 records tend to be larger. The pattern is similar
for all satellites. The overall negative bias for satellites 1, 2, 3, 5, and 6
is also evident from the plots. The patterns for the starboard antennas are
similar (not shown).

Figure 8: CYGNSS minus Jason-3 vs. the average of CYGNSS and Jason-3 for pairs
of observations. Each pair represents the two closest observations from each
two hour window over the entire year, provided that the distance is less than
25 km.

##  IV Conclusions

Our main finding is that during our study period of September 2019 to
September 2020, persistent biases existed among the wind speed measurements
recorded by the eight CYGNSS satellites and between some CYGNSS satellites and
Jason-3. Considering the averages of the model-based parameter estimates over
the study period, the largest bias between pairs of CYGNSS sensors was 1.11
m/s (CYGNSS 4 starboard minus CYGNSS 1 port), and the largest CYGNSS vs.
Jason-3 bias was -0.94 m/s (CYGNSS 1 port minus Jason-3). We discovered
smaller biases between the starboard and port antennas within a satellite,
with the largest average bias being 0.25 m/s (CYGNSS 5 starboard - CYGNSS 5
port).

It is not surprising to us that the two sensors on the same CYGNSS satellite
would be reasonably well-calibrated with respect to one another. Since they
tend to measure wind speeds at fairly close locations, direct comparison
across antennas is easier. However, direct comparisons between CYGNSS
satellites are more difficult due to the fact that some pairs nearly always
measure wind speeds at disparate locations. Similar to other studies that make
indirect comparisons with forecast winds or buoy data, we use Jason-3 as an
intermediary to achieve indirect comparisons among every pair of CYGNSS
satellites. We believe that Jason-3 data is appropriate because both CYGNSS
and Jason-3 attempt to measure snapshots of wind speed within small space-time
windows and they pass each other somewhat regularly.

These findings were facilitated by the use of a Gaussian process model that
contained parameters directly related to the expected difference between
measurements from different instruments taken at the same time and location.
In addition to the bias parameters, the models contained a parameter related
to the variance of the difference between two observations taken by different
sensors at the same time and location. We did not find substantial differences
in the estimates of noise parameters among the models for the eight CYGNSS
satellites. The size of the noise averaged about 0.75 m/s, meaning that the
noise was of roughly equal magnitude to the largest biases, implying that
about half of the mean square errors between measurements from different
sensors was due to bias, half due to noise. In other words, eliminating the
bias could potentially reduce the mean squared errors by a factor of two.

We validated the model-based findings with a more traditional empirical
analysis that searched for matching pairs of observations from different
sensors within small space-time windows. The general pattern of the estimates
in the empirical analysis was similar to that in the model-based analysis. The
analyses differed in that the week-to-week variation in the empirical
estimates of the bias was larger. This is expected because the model-based
estimates use generalized least squares, which provides variance-minimizing
parameter estimates under the assumed model.

Some studies have suggested that the size and direction of the bias may depend
on the magnitude of the true wind speed. We caution against over-interpreting
these results because we do not have access to the true wind speed. If a
particular sensor is chosen as the reference, we expect to see positive bias
at low wind speeds and negative bias at high wind speeds due to regression-to-
the-mean effects, even if the bias does not vary with wind speed. We attempted
to mitigate this effect by comparing paired differences against paired
averages, finding that generally the Jason-3 wind speeds increase relative to
CYGNSS as the average of their two measurements increases. This aspect is
certainly worthy of more exploration, with consideration of the aforementioned
statistical issues. Due to the differing accuracies of CYGNSS and Jason-3, the
average may not be the best estimate of the wind speed. We could also seek out
a third wind speed measurement to serve as the baseline. One could also pursue
model-based estimates of biases that depend on true wind speed. To this end,
consider the following extension of our model:

$\displaystyle Y_{i}$
$\displaystyle=\mu_{i}+a_{k(i)}+b_{k(i)}Z(x_{i},t_{i})+\varepsilon_{i}$

which contains a sensor-dependent slope multiplying the wind speed anomaly
$Z(x_{i},t_{i})$. Inferences about biases that depend on wind speed could be
obtained via estimation of $a_{k(i)}$ and $b_{k(i)}$. This is still a Gaussian
process model, so we could use the same methodology to fit the model, though
one would have to be careful about non-identifiability of parameters; for
example, the variance of $Z(x_{i},t_{i})$ is not identifiable separately from
the $b_{k(i)}$ parameters.

One could imagine that the bias depends on various other factors, such as
latitude, time, or GNSS satellite. This sort of variation could be handled
within our model framework by adding interactions between the bias and the
desired factor. To capture biases that vary in space, we could extend our
model as

$\displaystyle Y_{i}$
$\displaystyle=\mu_{i}+a_{k(i)}+Z(x_{i},t_{i})+W_{k(i)}(x_{i})\varepsilon_{i}$
(8)

where $W_{1},W_{2}$, and $W_{3}$ are independent spatial Gaussian processes.
Then we can interpret $a_{2}-a_{1}+W_{2}(x_{i})-W_{1}(x_{i})$ to be the
spatially-varying starboard bias, and $a_{3}-a_{1}+W_{3}(x_{i})-W_{1}(x_{i})$
as the spatially-varying port bias. We suspect that we would need to use more
than one week of data at a time to accurately estimate a spatially-varying
bias.

##  V Supplementary Materials

We run our analysis using R 4.0.5 [21] on platform x86_64-w64-mingw32/x64
running under Windows >= 8 x64. We make frequent usage of R packages "fields"
[22], "maps" [23], and "ggplot2" [24] for creating visualizations, and "GpGp"
[19] for modelling.

We maintain a Github repository at https://github.com/WillBekerman/satellite-
wind-speeds which contains all data, as well as R scripts to replicate our
analysis.

## Acknowledgment

The authors would like to thank David Moroni for technical support downloading
CYGNSS data.

## References

  * [1] Ibidun O Adelekan.  A survey of rainstorms as weather hazards in southern nigeria.  Environmentalist, 20(1):33–39, 2000.
  * [2] M Smith and L Chow.  Prediction method for aerodynamic noise from aircraft landing gear.  In 4th AIAA/CEAS Aeroacoustics Conference, page 2228, 1998.
  * [3] A Melese Endalew, C Debaer, N Rutten, J Vercammen, Mulugeta Admasu Delele, Herman Ramon, BM Nicolaï, and Pieter Verboven.  A new integrated cfd modelling approach towards air-assisted orchard spraying. part i. model development and effect of wind speed and direction on sprayer airflow.  Computers and Electronics in Agriculture, 71(2):128–136, 2010.
  * [4] Scott Eichelberger, James McCaa, Bart Nijssen, and Andrew Wood.  Climate change effects on wind speed.  North American Windpower, 7:68–72, 2008.
  * [5] Liliana Rusu, Alina Beatrice Raileanu, and Florin Onea.  A comparative analysis of the wind and wave climate in the black sea along the shipping routes.  Water, 10(7):924, 2018.
  * [6] IR Young, E Sanina, and AV Babanin.  Calibration and cross validation of a global wind and wave database of altimeter, radiometer, and scatterometer measurements.  Journal of Atmospheric and Oceanic Technology, 34(6):1285–1306, 2017.
  * [7] Jaime Daniels, Wayne Bresky, Steve Wanzong, Chris Velden, and Howard Berger.  GOES-R advanced baseline imager (ABI) algorithm theoretical basis document for derived motion winds, February 2019.
  * [8] Christopher S Ruf, Scott Gleason, Zorana Jelenak, Stephen Katzberg, Aaron Ridley, Randall Rose, John Scherrer, and Valery Zavorotny.  The CYGNSS nanosatellite constellation hurricane mission.  In 2012 IEEE International Geoscience and Remote Sensing Symposium, pages 214–216. IEEE, 2012.
  * [9] C Ruf, M Unwin, J Dickinson, R Rose, D Rose, M Vincent, and A Lyons.  CYGNSS: Enabling the future of hurricane prediction [remote sensing satellites].  IEEE Geoscience and Remote Sensing Magazine, 1(2):52–67, 2013.
  * [10] Shakeel Asharaf, Duane E Waliser, Derek J Posselt, Christopher S Ruf, Chidong Zhang, and Agie W Putra.  CYGNSS ocean surface wind validation in the tropics.  Journal of Atmospheric and Oceanic Technology, 38(4):711–724, 2021\.
  * [11] Daniel Pascual, Maria Paola Clarizia, and Christopher S Ruf.  Improved CYGNSS wind speed retrieval using significant wave height correction.  Remote Sensing, 13(21):4313, 2021.
  * [12] Faozi Saïd, Zorana Jelenak, Jeonghwan Park, and Paul S Chang.  The noaa track-wise wind retrieval algorithm and product assessment for CYGNSS.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–24, 2021\.
  * [13] Mikael Kuusela and Michael L Stein.  Locally stationary spatio-temporal interpolation of argo profiling float data.  Proceedings of the Royal Society A, 474(2220):20180400, 2018.
  * [14] Jouni Susiluoto, Alessio Spantini, Heikki Haario, Teemu Härkönen, and Youssef Marzouk.  Efficient multi-scale Gaussian process regression for massive remote sensing data with satgp v0. 1.2.  Geoscientific Model Development, 13(7):3439–3463, 2020.
  * [15] Matthias Katzfuss, Joseph Guinness, Wenlong Gong, and Daniel Zilber.  Vecchia approximations of Gaussian-process predictions.  Journal of Agricultural, Biological and Environmental Statistics, 25(3):383–414, 2020.
  * [16] Nick A Rayner, Renate Auchmann, Janette Bessembinder, Stefan Brönnimann, Yuri Brugnara, Francesco Capponi, Laura Carrea, Emma MA Dodd, Darren Ghent, Elizabeth Good, et al.  The EUSTACE project: delivering global, daily information on surface air temperature.  Bulletin of the American Meteorological Society, 101(11):E1924–E1947, 2020.
  * [17] Weitong Ruan, Adam B. Milstein, William Blackwell, and Eric L. Miller.  Multiple output Gaussian process regression algorithm for multi-frequency scattered data interpolation.  In 2017 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), pages 3992–3995, 2017.
  * [18] Joseph Guinness.  Permutation and grouping methods for sharpening Gaussian process approximations.  Technometrics, 60(4):415–429, 2018.
  * [19] J Guinness and M Katzfuss.  GpGp: Fast Gaussian process computation using Vecchia’s approximation.  R package version 0.4.0, 2018.
  * [20] Aldo V Vecchia.  Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society: Series B (Methodological), 50(2):297–312, 1988.
  * [21] R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2013.
  * [22] Douglas Nychka, Reinhard Furrer, John Paige, and Stephan Sain.  fields: Tools for spatial data, 2017.  R package version 11.6.
  * [23] Original S code by Richard A. Becker, Allan R. Wilks. R version by Ray Brownrigg. Enhancements by Thomas P Minka, and Alex Deckmyn.  maps: Draw Geographical Maps, 2018.  R package version 3.3.0.
  * [24] Hadley Wickham.  ggplot2: Elegant Graphics for Data Analysis.  Springer-Verlag New York, 2016.

|  William Bekerman is currently pursuing the B.S. degree in biometry and
statistics from Cornell University, Ithaca, NY, USA. He is involved in
research with the National Biomedical Center for Advanced ESR Technology and
the Department of Statistics and Data Science at Cornell University, as well
as the Department of Biostatistics at Columbia University. In 2019, he was
involved in summer research with the Department of Applied Statistics at New
York University. His research interests include statistical learning theory,
high-dimensional statistical inference, spatial statistics, signal and image
processing, along with applications to genomics, demography, geosciences, and
other disciplines. Mr. Bekerman is also a Student Member of the American
Statistical Association and the Biophysical Society and a member of the Hunter
R. Rawlings III Cornell Presidential Research Scholars program.
---|---
|  Joseph Guinness received his undergraduate degree in Math and Physics from
Washington University in St. Louis in 2007 and his Ph.D. in Statistics from
University of Chicago in 2012. He worked as a postdoctoral scholar and then
Assistant Professor at North Carolina State University in the Department of
Statistics from 2012 to 2017 and is now Associate Professor in the Department
of Statistics and Data Science at Cornell University in Ithaca, New York. His
research interests include computational statistics for spatial-temporal data,
with applications in the Earth sciences and remote sensing.
---|---