\documentclass[letter,12pt,dissertation]{OUdissertation}
\usepackage{graphicx, subfigure} %allows .eps and .epsi graphics to be inserted
\usepackage{epic} %allows use of latex graphics
\usepackage{eepic} %allows use of latex graphics
\usepackage{amssymb,amsmath}
\usepackage{amsfonts}
\usepackage[amssymb]{SIunits}
\usepackage{multirow}
\usepackage{psfrag}
\usepackage{float}
\usepackage{bar}
\usepackage[centerfoot]{pageno}
\usepackage{rotating}
\usepackage{xcolor}
\usepackage{url}
\setlength{\footskip}{0.365in}

\let\autoref\ref

\usepackage{natbib} %with amermeteorsoc.bst
\bibpunct{(}{)}{;}{a}{}{,}

\newcommand{\be}{\begin{equation}} %these definitions save typing
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\pd}[2]{\frac{\partial#1}{\partial#2}}
\newcommand{\pdd}[3]{\frac{\partial^{2}#1}{\partial#2 \partial#3}}
\newcommand{\bse}{\begin{subequations}}
\newcommand{\ese}{\end{subequations}}
\def\bal#1\eal{\begin{align}#1\end{align}}%allows for shortcut to align

\begin{document}
\title{RADAR POLARIMETRY FOR BIOLOGICAL APPLICATIONS}
\author{PHILLIP MICHAEL STEPANIAN}
\depositdate{2015}
\majorfield{SCHOOL OF METEOROLOGY}
\memberone{Dr. Phillip Chilson, Chair}
\membertwo{Dr. Robert Palmer}
\memberthree{Dr. Guifu Zhang}
%next two members are used in dissertation
\memberfour{Dr. Terry Schuur}
\memberfive{Dr. Jeffrey Kelly}
%%%%
\begin{preface}
\prefacesectionX{\normalfont{\large Dedication}}
To my family

%\prefacesection{Acknowledgements}



%must leave a blank line next to avoid bug:

\tableofcontents\pagestyle{plain}
%\listoftables
\listoffigures\pagestyle{plain}
\prefacesection{Abstract}

Radar aeroecology has steadily progressed since its inception in the 1950's, and additional advancements in radar technology, networks, modeling, and validation continue to drive the science further.  As human impacts steadily increase within the airspace and across the landscape, conflicts with wildlife will motivate new solutions for mitigating these negative interactions.  As one of the few sources of widespread surveillance of the aerosphere, radar will continue to play a valuable role in observing and quantifying animal life aloft.

The United States' next-generation weather radar network (NEXRAD) provides national coverage of the airspace at ten-minute resolution.  These data are quality-controlled, archived, and freely available for download, resulting in an efficient source of animal observations.  The NEXRAD weather surveillance radars have recently been upgraded to dual-polarizations, yielding three additional routine data products: differential reflectivity ($Z_{DR}$), differential phase ($\psi_{DP}$), and co-polar cross-correlation coefficient ($\rho_{HV}$).  While much speculation has been generated over the application of these products to biological studies, little work has been done on the subject.

The topic of this dissertation is the application of NEXRAD to biological studies, specifically focusing on the interpretation and use of polarimetric radar products. An overview of the polarimetric products is presented, and their biological interpretation is outlined.  Laboratory radio scattering measurements are detailed, providing a technique for characterizing the radar characteristics of individuals.  Typical manifestations of biological scatter are shown, including the roost emergences of bats and birds, as well as widespread nocturnal migration.  From these polarimetric signatures in widespread migration, a method is developed that extracts profiles of migrant flight orientations by exploiting morphological patterns in $\rho_{HV}$.  Finally, two validation methods are described using radar simulation and passive acoustic localization of nocturnal flight calls.

Radar polarimetry---especially in biological applications---is a young area of research, with much work still to be done.  As scientists and engineers continue to probe the polarimetric products, future developments in our understanding will enhance the use of polarimetric information as a tool in the study of airborne organisms.  This multidisciplinary effort will require the longterm cooperation of radar specialists and biologists to create and apply useful tools for addressing challenges in ecology and conservation.

\end{preface}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Introduction}
\label{intro}

\section{Motivation}
\label{intro:motivation}
The science of meteorology was born out of a need to describe and predict the impact of atmospheric processes on a specific group of terrestrial animals---namely, humans.  Whether forecasting travel conditions, anticipating agricultural yields, or mitigating risk to life and property, attempts in understanding the atmosphere have been largely motivated by human interests.  While efforts within the pure atmospheric sciences have increased our mechanistic knowledge of the atmosphere, resulting advances are often cast into the context of implications on human life and wellbeing.  

Growing along a different branch of science was a desire to understand the other life covering the planet.  Ecologists sought to characterize the interactions among plants, animals, and the underlying habitats that they occupy.  Scientists probed forests, prairies, deserts, and tundras, analyzing the living organisms within.  Beyond the terrestrial landscape, oceanographers studied marine ecosystems---the ebb and flow of tidal regions, the persistence of currents, the turbulent mixing of estuaries, and the many other complexities of fluid habitats.  As ecology advanced, understanding of the intricate web of terrestrial and oceanic systems continued to take shape.

Despite the developments in ecological science, a major habitat remained largely overlooked.  While land covers approximately 28\% of Earth's surface and the oceans cover 72\%, the atmosphere envelopes 100\% of Earth's surface.  Considering that some animals can be observed flying at altitudes in excess of ten kilometers, the atmosphere is by far the largest expanse of habitat on Earth.  But unlike oceanography, meteorology rarely acknowledged the wildlife living within the dynamic fluid body of air.  This natural overlap between atmospheric science and ecology has spawned several cross disciplines such as aerobiology, biometeorology, and aeroecology. 

It is undeniable that the airspace is home to a wide range of atmospheric flora and fauna; however, it is only relatively recently that humans have developed technology to occupy this habitat. As the impact of humans on the ecosphere continues to proliferate, there is an increasing need to mitigate human-wildlife conflicts on land and within the aerosphere.  Such conflicts include anthropogenic inclusions or occupancy within the airspace by aircraft, high-rise structures, towers and antennae, wind turbines, and transmission lines.  Other effects include the modification of land surfaces for agriculture, transportation, or human habitation, all of which contribute to wildlife habitat changes and fragmentation.  Perhaps the most notable influence of humans is the anthropogenic modification of atmospheric composition, which has lead to shifting climatic zones and animal ranges.  

The combination of these effects has created an urgent demand for the application of weather and climate science to the greater Earth ecological system.  One facet of this endeavor is the development of high spatiotemporal resolution observational tools that can provide long-term surveillance of both the atmosphere itself, as well as the wildlife within.  These observations must include measurements of atmospheric fluid dynamics; radiative transfer and thermodynamics; and the contribution and distribution of atmospheric constituents---both chemical (e.g., water vapor, carbon dioxide, methane, and other trace gases) and biological (e.g., pollen, spores, insects, birds, and bats).  It is with these observations that scientists can characterize the ecology of the atmosphere and the changes that are ongoing.  This knowledge can ultimately provide guidance for future use of the airspace and policies protecting the wildlife living in it. 

\section{Goals}
\label{intro:goals}
Radar is a proven source of both atmospheric and biological observations.  Efforts in adapting radar technology to weather surveillance have resulted in advanced systems and sophisticated techniques for atmospheric monitoring; however, biological radar applications are comparatively rudimentary.  The purpose of this dissertation is the development and introduction of several ecological monitoring methods to be applied to operational weather surveillance radars.  While some dedicated biological radar systems are currently available on the market, the utilization of serendipitous animal observations by operational systems provides the greatest benefit at the lowest cost to ecologists.  The scope primarily focuses on birds, however, some applications for bats and insects are included.  The ultimate goal of this dissertation is to provide novel observational methods for monitoring airborne organisms in support of research goals within ornithology, chiropterology, and entomology.



\section{Structure of the Dissertation}
\label{intro:structure}
The remaining dissertation chapters are organized as follows.  Chapter~\ref{background} provides background material on the past and present use of radar for biological surveillance. Additionally, a brief overview of the theory and products provided by polarimetric pulsed Doppler radar are presented.  Chapter~\ref{interp} introduces the biological interpretation of polarimetric radar products, as well as their deviations from the traditional meteorological interpretations. Starting with the simplest case, Chapter \ref{measure} describes radio wave scatter from isolated individuals by considering laboratory measurements of a bird and bat specimen. Extending this scattering theory, Chapter~\ref{NEXRAD} explores the resulting features of bulk biological scatter from operation weather radars. Chapter~\ref{align} describes a technique for extracting flight orientation profiles of migrating organisms.  Chapter~\ref{sim} combines the theoretical aspects of biological scattering and radar operation to formulate a Lagrangian radar simulator, and compares synthesized radar signals to observations to assess the validity of the underlying theory and assumptions.  Finally, Chapter~\ref{local} introduces bioacoustics as a source of ground truth observations that can link radar measurements to the underlying bird species composition. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Background}
\label{background}

\section{Historical Perspectives}
\label{background:history}
Modern radar was developed during the second world war as a military countermeasure for advanced warning of incoming enemy aircraft \citep{Buderi1998}.  An unintended result of probing the atmosphere with radio waves was the detection of non-aircraft signals, which later were identified as either precipitation or so-called `angels' \citep{Atlas1959}.  The value of observations of weather systems motivated the adoption of radar technology in meteorological studies, and has contributed to major operational and research applications.  The source of this latter class of `angels' was a subject of controversy, with equally many believing in meteorological, inanimate particulate, and volant biological causes \citep{Plank1956, Tolbert++1958, Kocurek+1967, Chadwick+1983}.  There was unquestionable proof that airborne organisms could be---and were indeed---detected in radar measurements \citep{Lack+1945, Crawford1949}, but the extent to which this occurred was often subject to speculation in the absence of routine validation \citep{Chadwick+1983}.  This uncertainty is likely the reason that radar technology was not readily adopted into the biological sciences, as it was in meteorology.

Subsequent radar studies on biological subjects were primarily performed as a novelty topic in electrical engineering. Laboratory measurements performed by \cite{Edwards+1959} investigated the dependence of radar cross section on polar look angle for a pigeon (\emph{Columba ivia}), starling (\emph{Sturnus vulgaris}), and sparrow (\emph{Passer domesticus}). Additionally, the scattering contribution of feathers was investigated by measuring a pigeon before and after plucking, and effects of wing position were investigated for a rook (\emph{Corvus frugilegus}).  \cite{Blacksmith+1965} performed crude lab measurements of ducks and chickens, but as with several similar studies of the time, emphasis was on novelty of the topic rather than serious analysis or biological application.  Several studies attempted radar cross section measurements by releasing live insects and birds from aircraft within a radar sampling volume \citep{Glover++1966,Konrad++1968}, others performed measurements on dead specimens tethered to lofted balloons \citep{Richter+1973}, and some attempted such characterizations on wild birds flying within the airspace \citep{Eastwood+1966}. 

Many of the other studies on biological radar scatter during this period were performed by meteorologists within the context of weather research. Atmospheric scientists used small insects as passive tracers of clear-air atmospheric motions \citep{Lhermitte1966, Richter++1973}, while radar meteorologists performed studies to characterize and censor biological clutter.  Through the latter half of this period several biologists began applying radar measurements to pure biological studies, effectively founding the modern science of radar aeroecology \citep{Eastwood1967, Bruderer1969, Gauthreaux1970, Able1970, Alerstam1972, Williams++1972, Alerstam+1973}.

The following decades generated a growth in biological radar studies, primarily consisting of an overlap among three distinct areas:
(1) the mechanistic exploration of biological radio wave scatter by electrical engineers and radio physicists, 
(2) the application of radar measurements to biological studies by ecologists, and 
(3) the study of biological radar echoes by meteorologists for clear-air studies or clutter suppression.  Some highlights of this period included the use of tracking-radar baseband signals for wing-beat characterization \citep{Schaefer1968, Bruderer+1972, Vaughn1974}, multi-radar networks for larger scale migratory studies \citep{Gauthreaux1971}, Doppler observations of flying animals \citep{Martinson1973}, and polarimetric measurements \citep{Mueller+1985, Riley1985}.  It was also during this period that the concept of radar-based taxonomic classification began to gain traction as the future of biological radar \citep{Martinson1973,Williams+1980, Larkin1980}.

Developments through the 1990's established the contemporary landscape of radar aeroecology. These areas included the use of the newly upgraded weather surveillance radar network in the United States \citep{Gauthreaux+1998, Russell+1998}, the ornithological tracking and identification efforts by the Swiss \citep{Bruderer1994, Liechti++1995, Bruderer1997a, Bruderer1997b} and Swedes \citep{Alerstam+1999}, British entomological studies \citep{Riley+1990}, deployments of portable marine radars for ecological studies \citep{Riley+1990, Cooper++1991}, and radar-based aircraft birdstrike mitigation \citep{Haykin++1991}.  Additionally, studies focusing on meteorological data quality and echo classification continued to gain importance \citep{Wilzak++1995, Russell+1997, Zrnic+1998}.


\section{State of the Art}
\label{background:sota}
Current advances in radar aeroecology are progressing across several research fronts.  Radio wave scatter by animals is still being approached as a measurement and modeling topic in electromagnetics \citep{Lang++2004, Bachmann+2007, Melnikov++2012, Melnikov++2014, Melnikov++2015}.  Some studies are focused on development of practical radar-based retrieval methods and method development \citep{Larkin++2002,Martin+2007,Schmaljohann++2008,Nebuloni++2008,Horn+2008,Zaugg++2008,Buler+2009,Taylor++2010,Cruz++2013}.  Others are finding ways to exploit larger radar networks for ecological applications \citep{Dokter++2011, Chilson++2012}.  Currently, weather radar networks across Europe (OPERA, Fig.~\ref{back:opera}) and the United States (NEXRAD, Fig.~\ref{back:NEXRAD}) are being used for studies in biological sciences.


\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/back_opera}
  \caption{The OPERA network. Image from http://www.eumetnet.eu/radar-network.}
\label{back:opera}
\end{figure}



A large volume of work is applying radar to basic ecology questions in migration \citep{Alerstam++2001,Dinevich++2003,Hedenstrom++2009, Dokter++2013,ONeal++2014}, animal behavior \citep{ShamounBaranes++2010,Dokter++2013b,Diehl2013, VanDoren++2015}, and phenology \citep{Kelly++2012}.  Additionally, radar is playing a significant role in animal conservation, including land management \citep{Bonter++2009, Buler+2014} and wind farm siting \citep{Huppop++2006,Plonczkier+2012}.  Efforts directed toward human health and well-being include agricultural crop pest monitoring \citep{Westbrook2008, Bell++2013} and aircraft birdstrike prevention \citep{Zakrajsek+2001,Nohara++2011}.  Radar studies of bioscatter is still an active topic in meteorology, with applications related to clear air studies \citep{Contreras+2008}, hydrometeor classification \citep{Park++2009, Chandrasekar++2012, AlSakka++2013}, and data quality \citep{Martiner+2001, Zhang++2005, Lakshmanan++2010})

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/back_nexrad}
  \caption{The NEXRAD network covering the contiguous United States. Image from http://virtual.clemson.edu/birdrad/com1a.htm.}
  \label{back:NEXRAD}
\end{figure}

Much of the increase in radar ecological studies can be attributed to an increase in radar technology and data availability. Several aeroecology labs have developed specialized biological radar units, such as the mobile marine radar system of \cite{Cooper++2001}, the vertically-looking entomological radar of \cite{Chapman++2003}, and the Swiss Superfledermaus \citep{Bruderer++2010}.  

The applications related to to conservation and air safety have motivated several commercially available dedicated biological radar systems \citep{Nohara++2007}.  These units and software are produced by DeTect, Inc. (MERLIN\texttrademark, \cite{Kelly++2007}), Accipiter, Inc. (Avian Radar, \cite{Nohara++2007}), SRC  Inc. (BSTAR\texttrademark~and SR Hawk\texttrademark, \cite{SRC}), Versar (MARS$^\circledR$, \cite{VERSAR}), Robin (3D FLEX, \cite{ROBIN}), and others.

While specialized biological radar systems are becoming more common, much of the research is still conducted using atmospheric radars, including the polarimetric NEXRAD network in the U.S. \citep{Chilson++2012}, the OPERA network in Europe \citep{ShamounBaranes++2014}, and other serendipitous sites.


\section{Weather Surveillance Radar}
\label{background:wsr}
Radar samples the atmosphere by directing radiation through an antenna to the region of the airspace that is of interest and recording the backscattered signal \citep{DoviakZrnic}. These resulting signals are transformed into data products, sometimes generally referred to as \emph{measurables} or \emph{observables}, that correspond to specific sampling volumes located radially along the antenna beam direction. That is, the retrieved information is a function of antenna pointing direction in azimuth and elevation, as well as range from the radar. Depending on the position and motion of the antenna, different spatial and temporal information may be extracted from the radar data. This discussion will focus on radar systems that operate in surveillance modes, as opposed to those that actively move the antenna to follow a defined object, that is, tracking radar. A number of different sampling techniques have been developed and optimized for surveillance radar systems to extract specific information in time and space based on the application needs. In particular, three sampling strategies commonly utilized on radar platforms include fixed beam, azimuth scanning and elevation scanning.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/back_scan}
  \caption{Physical sampling technique (left column), data raster image (center column) and physical coordinate image (right column) for (a) fixed-beam
sampling, (b) elevational scanning and (c) azimuthal scanning. \citep[reproduced from][]{Stepanian++2014}}
\label{back:scan}
\end{figure}

Fixed-beam, or spotlight, sampling is the most basic technique, requiring no antenna motion. In this case, the radar is pointed in a fixed direction, often vertically or near vertically, and repeatedly samples the airspace (Fig.~\ref{back:scan}a, left). As organisms pass through the sampled volumes, each pulse returns a profile of the radar measurables as a function of range along the beam (Fig.~\ref{back:scan}a, center). While each sample returns a one-dimensional profile of measurables, consecutive samples are typically displayed in series as a function of time to form a range-time indicator, or RTI display (Fig.~\ref{back:scan}a, right). In addition, it is common that the results of several consecutive pulses are averaged together to reduce the contamination by random noise, thereby improving the data product quality. It is important to note that the RTI display does not directly contain any two-dimensional spatial information, but only the temporal evolution of the sampled radial. Under specialized circumstances, however, it may be possible to infer spatial information by making assumptions regarding the motion of the individuals in the radar space. For example, if a large flock of birds passes over the radar with their formation unchanging in time, the temporal changes in the RTI plot may be used to infer spatial changes, that is, the flock formation. This assumption may be valid in some ecological situations in which the heading, speed and altitude of individual organisms remain fairly constant as they travel over large horizontal extents (i.e. large scale migrations). In many cases, an organism may remain within the beam for several consecutive pulses, resulting in a horizontal line of increased intensity across the RTI display (e.g. Schmaljohann et al. 2008, Fig. 2). Additional examples of fixed-beam sampling in ecological applications have been demonstrated by \cite{Moran++2000}, \cite{Chapman++2003}, \cite{Martin+2007}, and Dokter et al. (2013; 2013b).

In the second technique, the antenna scans in elevation at a fixed horizontal azimuth angle. In some cases, the antenna will scan from one horizon, through the vertical, to the opposite horizon (e.g. Fig.~\ref{back:scan}b, left), while other times, only a subset of elevation angles are covered. These resulting radar data products are a function of elevation angle and range from the radar (sometimes referred to as an E-Scope, Fig.~\ref{back:scan}b, center). Unlike fixed-beam sampling, elevation scanning contains explicit spatial information with the surveyed airspace encompassing a vertical cross section and is typically displayed as such in the range-height indicator (RHI) format (Fig.~\ref{back:scan}b, right). \cite{Gauthreaux1991}, \cite{Martin+2007} and \cite{vanGasteren++2008} provide ecological demonstrations of elevational scanning. 

Perhaps the most common sampling method in radar applications is the azimuthal rotation of the antenna (Fig.~\ref{back:scan}c, left). If equipped with an antenna that produces a fan beam (with the broadest portion of the beam aligned vertically), the antenna is typically rotated at a fixed elevation angle, providing instantaneous surveillance across a range of heights. Common to marine radar systems and airport surveillance radars, this scanning strategy provides rapid frame-to-frame updates, but cannot provide the altitude of the detected objects \citep{Larkin+2012}. In the case of weather radars, a parabolic antenna creates a highly focused beam, which is swept azimuthally through a number of elevation angles. A typical beam width for weather radars is about one degree. While these volume coverage patterns lead to slower update times, they provide three-dimensional coverage of the airspace. In both cases, these methods yield information about the horizontal distribution of organisms, with radar products being a function of azimuth angle and range from the radar (sometimes referred to as a B-Scope, Fig.~\ref{back:scan}c, center). In the case of weather radar, the three-dimensional data products, often referred to as volume data, can be thought of as a collection of two-dimensional scans (as in Fig.~\ref{back:scan}c, center), each at a different elevation angle. These two-dimensional, constant elevation angle data products � often referred to as sweeps, tilts or cuts � are typically displayed as viewed from above in the plan position indicator (PPI) (Fig.~\ref{back:scan}c, right). Azimuthal scanning has been applied widely within the ecology literature, with \cite{Horn+2008}, \cite{Buler+2009}, and Dokter et al. (2010; 2013b) demonstrating several applications.


\section{Dual-polarization pulsed Doppler Radar}
\label{background:dualpol}
Many conventional radar systems transmit and receive radiation in a single polarization, while some use polarization diversity to gain additional information on the characteristics of scatterers and the propagation medium \cite{DoviakZrnic}.  There are several schemes for making polarimetric measurements, but this dissertation will be limited to orthogonal dual-polarizations.  That is, waves are transmitted in the horizontal polarization plane with oscillations parallel to the horizon, as well as the orthogonal vertical plane.  Additionally, two modes of orthogonal dual-polarization will be considered: Alternating Transmit and Receive (ATAR) and Simultaneous Transmit and Receive (STAR).  In ATAR, a horizontally polarized pulse is transmitted, scattered, and received.  Subsequently, a vertical pulse is transmitted, scattered, and received, and the process repeats.  Radars operating in STAR configuration transmit a single pulse containing both horizontal and vertical components, which is scattered and received as a mixed-polarization packet.  These mixed polarization signals are finally sequestered upon reception. The formal definitions of polarimetric observables are given in Section~\ref{background:products} .


\section{Radio Scattering of Isolated subjects}
\label{background:radio}
When a dual-polarized incident electromagnetic field ($\vec{\textbf{E}}_i$) impinges on a scatterer in space, the backscattered field ($\vec{\textbf{E}}_b$) can be described by the complex backscattering matrix for the scatterer
\begin{equation}
\label{back:smatrix}
\textbf{S}= 
\begin{bmatrix}
s_{hh} & s_{hv} \\
s_{vh} & s_{vv}
 \end{bmatrix}
 \end{equation}
 with $h$ and $v$ denoting horizontal and vertical polarizations, respectively  \cite{DoviakZrnic}.  The first subscript index specifies the backscattered polarization, and the second subscript index specifies the incident polarization.  For example, the complex function $s_{vh}$ represents the backscattered vertical polarization component of an incident horizontally polarized wave. Thus, matching subscript indices (i.e., $s_{hh}$ and $s_{vv}$) denote copolar scattering functions, while different subscript indices are cross-polar scatter functions.  For a scatterer located a range, $r$, from the antenna, the backscattered field at the receiver can be calculated as
 \begin{equation}
 \label{eq:waveScatter}
 \vec{\textbf{E}}_b= \textbf{S}\vec{\textbf{E}}_i\frac{\exp{(-jkr)}}{r}
 \end{equation} 
with $k$ denoting the probing radar wavenumber, $\frac{2\pi}{\lambda}$ \citep{DoviakZrnic}.  When describing the relationship between the incident and backscattered fields for a given scatterer, two quantities are typically reported.  The backscatter radar cross-sections
 \begin{subequations}
 \begin{align}
 \sigma_{b,hh} &= 4\pi |s_{hh}|^2\\
 \sigma_{b,vv} &= 4\pi |s_{vv}|^2\\
 \sigma_{b,hv} &= 4\pi |s_{hv}|^2\\
 \sigma_{b,vh} &= 4\pi |s_{vh}|^2
 \end{align}
 \end{subequations} 
describe the amplitude of backscattered polarized energy provided for an incident polarization, relative to an isotropic scatterer \citep{DoviakZrnic}.  The phase shifts introduced upon scattering are simply the arguments of the scatter functions for the incident and backscattered polarizations,
 \begin{subequations}
 \begin{align}
 \gamma_{hh} &= \arg(s_{hh})\\
 \gamma_{vv} &= \arg(s_{vv})\\
 \gamma_{hv} &= \arg(s_{hv})\\
 \gamma_{vh} &= \arg(s_{vh})
 \end{align}
 \end{subequations} 
 with the differential phase shifts due to scattering, $\delta$, being the difference between two scatter phases.

\section{Definitions of Radar Moments and Polarimetric Products}
\label{background:products}
\subsection{Radar Reflectivity ($\eta$), Reflectivity Factor ($Z$), and \\Equivilant Reflectivity Factor ($Z_e$)}
\label{background:products:z}
Consider a single-polarization radar transmitting a pulse with power, $P_t$, and wavelength, $\lambda$, through an antenna with gain, $G$.  The returned power, $P_r$, depends on the backscattering cross-section of the scatterer $\sigma_b$, its range from the radar, $r$, and its position within the antenna beam pattern, $\vartheta$, and nominal range sample, $r_o$.  The result of this relation is the radar equation for point scatterers,
\begin{equation} 
P_{r} = \frac{P_t G^2 \lambda^2 \sigma_{b} f^4(\vartheta) |W^2(r,r_o)| } {(4 \pi)^3 r^4}\;\;\;\;\;\;\;\;\;\;\;\; \textrm{[W]}
 \label{back:radarEq}
\end{equation}
in which all units are in MKS and weighting functions are linear fractions between zero and unity \citep{DoviakZrnic}.  The weighting functions, $W$ and $f$, represent the position of the scatterer within the beam, and account for variations in sensitivity within the range sample ($W$) and beam width ($f$).  When multiple scatterers are in a single sample volume, their collective backscattering cross-section can be related to the resolution volume, yielding the radar reflectivity,
\begin{equation} 
\eta = \frac{1}{V}\sum_{i=1}^{vol} \sigma_{b,i}\;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{cm}^2~\textrm{km}^{-3}]~\textrm{or}~[\textrm{m}^{-1}]
 \label{back:eta}
\end{equation}
which summarizes the cumulative backscattering area within the volume, $V$, \citep{Russell+1997, Chilson++2012_2Z}.  In the case that scatterers are spherical dielectrics that are small compared to the probing wavelength, \cite{Mie1908} shows that the backscattering cross-section is related to the scatterer size and composition as
\begin{equation} 
\sigma_b = \frac{\pi^5}{\lambda^4}|K|^2D^6 \;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{m}^2]
 \label{back:mie}
\end{equation}
with $D$ being the scatterer diameter, and $|K|=|\frac{m^2-1}{m^2+2}$| with $m$ defined as the complex refractive index of the scatterer.  Within this framework, size of the scatterer is explicitly related to the amplitude of the radio scatter
\begin{equation} 
D = \left(\frac{4\pi\lambda^4}{\pi^5 |K|^2}|s|^2\right)^{1/6}\;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{m}]
\end{equation}
such that scatter amplitudes can unambiguously relate to diameter.

Substituting (\ref{back:mie}) into (\ref{back:eta}) yields the so-called Rayleigh approximation to volume-filling radio scatter
\begin{equation} 
\eta = \frac{\pi^5|K|^2}{\lambda^4V}\sum_{i=1}^{vol} D_i^6\;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{mm}^6~\textrm{m}^{-7}]~\textrm{or}~[\textrm{m}^{-1}]
 \label{back:rayleigh}
\end{equation}
with the Radar Reflectivity Factor defined as 
\begin{equation} 
z = \frac{1}{V}\sum_{i=1}^{vol} D_i^6\;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{mm}^6~\textrm{m}^{-3}]
 \label{back:z}
\end{equation}
\citep{DoviakZrnic}.  The power received by the radar is often calibrated to relate to scatterers within the Rayleigh assumption.  Additionally, because backscattered power can span many orders of magnitude, it is often reported in decibel units with respect to one unit of Z.  The process of range correction and scaling to $Z$ is computed as
\begin{equation} 
Z = 10\textrm{log}_{10}(P_r)+20\textrm{log}_{10}(r)+C\;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{dBZ}]
 \label{back:calZ}
\end{equation}
with a calibration constant, $C$ \citep{DoviakZrnic}.  In many cases, the scatterers sampled by the radar may not fulfill the requirements for the Rayleigh approximation to be valid.  Because these scatterers can still be related to a value of $Z$, a more general term of Equivalent Reflectivity Factor ($Z_e$) is sometimes used to describe the Mie scatterers that would correspond to these signals. 

Through the duration of this dissertation, a lowercase $z$ will refer to linar radar reflectivity factor (in mm$^6$ m$^{-3}$), while a capital letter $Z$ will refer to the decibel reflectivity factor in units of dBZ.  Additionally, subscripts will denote the polarization of reflectivity factor (i.e., $Z_H$ is horizontal reflectivity factor), with non-subscripted $Z$ implying the horizontal polarization.  

\subsection{Mean Radial Velocity ($v_r$) and Spectrum Width ($\sigma_v$)}
\label{background:products:v}
For the Doppler systems considered herein, a coherent receiver compares the signal phase upon reception with the transmit phase to detect phase shifts caused by the motion of scatterers \citep{DoviakZrnic}.  The motions contributing to these phase shifts are the components along the radial vector of the radar, or in other words, motions toward or away from the radar.  As signals are received at the radar, the contributions of backscattered power are inherently tied to the motion of the scattering objects.  That is, the total power can be decomposed by its frequency content into the contributions from scatterers of different velocities along the radial vector.  The result of this decomposition is the Doppler velocity spectrum, which describes the underlying frequency shifts of received power contributions \citep{DoviakZrnic}.  

The Doppler spectrum is the convolution of the frequency modulating effects from several sources.  The first source of frequency content is from the scatterers themselves; the velocity components along the radial vector determine the mean frequency shift, while the variation in velocity broadens the peak.  Other spectral variations are caused by sampling effects such as scanning motion and width of the beam, both of which broaden the overall spectrum.  Finally, the spectrum is also affected by computational artifacts like the choice of windowing functions and discrete Fourier transform length.  

In weather radar applications, the Doppler spectrum is assumed to follow a Gaussian distribution, and is characterized by the moments of the fitted distribution.   The zeroth moment, or integrated area under the curve, is the signal power, and is related to $Z$ and $\eta$.  The first moment is the peak of the Gaussian fit, and indicated the mean radial velocity of scatterers ($v_r$).  A typical convention is that negative values denote motions toward the radar.  The second central moment of the distribution is the variance of the Gaussian fit, and is reported as the spectrum width of the mean Doppler radial velocity ($\sigma_v$).  In general, spectrum width indicates the diversity of scatterer velocities within a sampling volume. 

As with any discrete sampling problem, the maximum frequency that can be resolved is determined by the sampling rate.  This maximum, Nyquist frequency relates to the maximum resolvable radial velocity magnitude and Pulse Repetition Frequency (PRF) as
\begin{equation} 
v_a = \pm\frac{\textrm{PRF} \lambda}{4}\;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{m~s}^{-1}]
 \label{back:va}
\end{equation}
in which velocities above $v_a$ or below $v_a$ will be aliased, or wrapped, to the opposite side of the Doppler spectrum \citep{DoviakZrnic}.

\subsection{Differential Reflectivity ($Z_{DR}$)}
\label{background:products:zdr}
 The differential reflectivity is the linear ratio (or logarithmic difference) between radar reflectivity factors of orthogonal polarizations.  For the case of weather radar,
 \begin{equation} 
Z_{DR}=Z_H-Z_V  \;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{dBZ}]
 \end{equation} 
with subscripts indicating the received polarization.  In meteorological applications, scatterers are generally assumed to follow the Rayleigh conditions, linking $Z_{DR}$ to the bulk aspect ratio of the scatterers \cite{DoviakZrnic}.  For larger scatterers, the conditions of the Rayleigh assumption are not valid, and $Z_{DR}$ has a more complicated interpretation as the amplitude from orthogonal polarizations oscillate independently.  Additional complications arise when scatterers are non-spherical and aligned \citep{Ryzhkov+2007}, and will be discussed in Section \ref{interp:zdr}.
 
 
\subsection{Differential Phase Shifts ($\phi_{DP}$, $\Phi_{DP}$, and $K_{DP}$)}
\label{background:products:phidp}
As radio wave propagates through a medium, the phase is continuously modified along the path.  One way of conceptualizing these phase shifts is by considering many forward-scattering events within the medium that each introduce a slight phase shift upon scattering \citep{Bohren+2006}.  When the medium is stratified, these forward scattering phase shifts are described in terms of ray refraction \citep{Bean+1966}.  When a medium is homogenous along the dimension perpendicular to wave propagation, these phase shifts are a continuous modification to the original phase.  The magnitude of these phase shifts is dependent on the refractive index, size, shape, and orientation of the scattering medium at the polarization of the propagating wave.  

Just as $Z_{DR}$ is useful for comparing amplitudes at dual polarizations, the phase difference at dual polarizations ($\phi_{DP}$) can be measured.  As discussed in Section~\ref{background:radio}, the modification of incident, dual-polarized radio waves upon scattering is described by (\ref{eq:waveScatter}), with the \textbf{S} matrix describing the phase shifts occurring due to scattering ($\delta$).  When considering the full radar process of transmission, scattering, and reception, two additional terms are added,
 \begin{equation}
 \vec{\textbf{E}}_b=  \frac{\exp{(-jkr)}}{r}\textbf{R}\textbf{S} \textbf{T}\vec{\textbf{E}}_i
 \end{equation} 
with $\textbf{T}$ describing the propagation of the transmitted wave of initial phase offset, $\psi_t$, along the path to the scatter as,
 \begin{equation}
 \textbf{T}=  
 \begin{bmatrix}
1 & 0 \\
0 & \exp(\psi_t+\frac{1}{2}\Phi_{DP})
 \end{bmatrix}
 \end{equation} 
 and \textbf{R} describing the propagation of the scattered wave to the radar receiver
  \begin{equation}
 \textbf{R}=  
 \begin{bmatrix}
1 & 0 \\
0 & \exp(\psi_r+\frac{1}{2}\Phi_{DP})
 \end{bmatrix}
 \end{equation}
 with $\psi_r$ denoting the phase offset introduced within the receiver \citep{Melnikov++2015}. In both terms, $\Phi_{DP}$ is the differential phase offset resulting from propagation through the medium between the radar and the scatterer. In total, the signal at the radar receiver is a combination of the initial transmitted field ($\psi_t$), the propagation to the scatterer ($\frac{1}{2}\Phi_{DP}$), the effects of scattering ($S$; or $\sigma_b$ and $\delta$), the propagation back to the radar ($\frac{1}{2}\Phi_{DP}$), and the effects within the receiver chain ($\psi_r$) \citep{Melnikov++2015}.  In most meteorological applications, the propagation medium is described by the incremental phase shifts in propagation, $K_{DP}$, in degrees per unit propagation distance such that
  \begin{equation} 
\Phi_{DP}=\int_{distance} K_{DP} \;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{deg}]
 \end{equation} 
in which $distance$ is the round trip propagation through the medium.  Additionally, the system differential phases on transmission and reception are not generally known, but their combined effect ($\psi_{sys}$) can be obtained on a case-by-case basis \citep{Melnikov++2015}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Co-polar Correlation Coefficient ($\rho_{HV}$)}
\label{background:products:rhohv}
When comparing orthogonal polarizations, $Z_{DR}$ considers differences in signal amplitude while $\phi_{DP}$ describes differences in signal phase---both of which are taken in the ensemble sense.  The correlation coefficient between co-polar channels, $\rho_{HV}$, indicates the overall similarity between the orthogonal, dual-polarization, co-polar signals,
   \begin{equation} 
\rho_{HV}=\frac{<s_{vv}s_{hh}^*>}{<|s_{hh}|^2>^{1/2}<|s_{vv}|^2>^{1/2}}
 \end{equation} 
 with $<\cdot>$ indicating ensemble averages \citep{DoviakZrnic}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Derived Products}
\label{background:devProd}

\subsection{Spatial Textures}
\label{background:devProd:texture}
When viewing images of radar products, it is often not only the magnitude of the variables that reveals information on scatterers, but also their spatial variation.  In some cases, the volume-to-volume continuity can indicate how homogenously or densely scatterers are distributed in space.  To this end, the spatial standard deviation of fields can provide additional information. \cite{Park++2009} uses a running 5-gate standard deviation of reflectivity factor to determine how variable values are along radials for hydrometeor classification.  Additionally, both 5-gate and 25-gate running standard deviations of $K_{DP}$ are calculated to describe local scatterer variability.  Texture parameters can be taken on any field across one or two dimensions.  In the latter case, a two dimensional structuring element is applied across both range and azimuth to obtain 2D spatial statistics. 

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/back_vpr.pdf}
  \caption{Vertical stratification of scatterers and the resulting PPIs for (a-b) near-ground enhanced scatterers, and (c-d) an elevated layer of enhanced scatterers.}
\label{back:vpr}
\end{figure}

\subsection{Vertical Profile of Reflectivity (VPR)}
\label{background:devProd:VPR}
Many phenomena in the atmosphere vary vertically, or in response to the inherent vertical stratification of pressure, temperature, and humidity.  When considering the three scanning techniques in Fig~\ref{back:scan}, it is clear that the altitude of the radar beam increases with range from the radar.  As a result, different ranges typically sample different vertical regions of the atmosphere. In some applications, such as fixed-beam vertical sampling, these vertical changes are the desired measurements.  In other applications, however, the desired measurements are meant to capture horizontal changes only.  For example, the azimuthal samples taken by weather radars are typically meant to measure horizontal features such as storms and fronts.  Figure~\ref{back:vpr} demonstrates two common situations in which vertical heterogeneity manifests as false horizontal structures. In the first case, the concentration of scatterers is highest near the ground and decreases in altitude (Fig.~\ref{back:vpr}a), resulting in a PPI with the greatest scatter intensity near the radar (Fig.~\ref{back:vpr}b).  The second case shows an elevated layer of enhanced scatter (Fig.~\ref{back:vpr}c), which results in a ring of higher reflectivity around the radar (Fig.~\ref{back:vpr}d).  In both situations, the distribution of scatterers is horizontally homogenous and vertically stratified, but the resulting PPIs imply horizontally-variable reflectivity structures. 

Sometimes it is necessary to account for these changes in the vertical beam location. For example, when estimating precipitation rates it is often the case that precipitation can be vertically variable.  In some situations, shallow precipitation may be well-represented at close ranges, but the beam may overshoot it farther away (as in Fig. \ref{back:vpr}a-b).  The opposite case will sometimes occur when high-based precipitation evaporates before reaching the ground, resulting in virga, and an apparent lack of precipitation near the radar (the opposite of Fig. \ref{back:vpr}a-b).  Finally, it is common that high elevation angles will result in such high beam propagation that the 0 $^\circ$C boundary is crossed and different precipitation types are recorded (e.g., rain at lower levels and ice aloft).  These cases will result in PPIs with discontinuities in range, and the occurrence of the \emph{radar bright band} \citep{DoviakZrnic}. 

Under the assumption that the atmosphere is horizontally homogenous along the beam, profiles of radar variables can be obtained by relating the height of the beam with range from the radar.  When the variable is radar reflectivity factor, these vertical profiles of reflectivity can be related to the density, composition, and distribution of scatterers in height.  Alone, these profiles are valuable in that they describe the vertical structure of the atmosphere; but they may also be used to correct or normalize for vertical variability such that the resulting PPIs do not show these artificial effects \citep{Koistinen1991}.  


\subsection{The Velocity Azimuth Display (VAD)}
\label{background:devProd:VAD}
As discussed in Section~\ref{background:products:v}, the velocity component resolved by Doppler radar is that along the radial vector, i.e., toward or away from the radar.  As a result, each ray only resolves one component of the overall wind field.  One way for weather surveillance radars to obtain more information on wind structure while scanning in azimuth is by applying the Velocity Azimuth Display (VAD) technique to piece together the flow field from wind components in different azimuths and elevations \citep{Browning+1968}.  Because the beam rises in altitude with range from the radar (as discussed in the previous section), the resulting retrievals yield a vertical profile of wind speed and direction.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Biological Interpretation of Radar Products}
\label{interp}
\section{Overview}
\label{interp:overview}
Doppler weather surveillance radars, such as those comprising the NEXRAD network in the United States, provide a wealth of information on the large-scale movements of airborne organisms, including birds, bats, and insects \citep{Chilson++2012}.  The increased prevalence of polarimetric radars---most notably, the recent NEXRAD upgrade to dual-polarizations---has resulted in enhanced capabilities for observing the physical and behavioral characteristics of airborne animals \citep{Zrnic+1999}.  The strong dependence of polarimetric measurements on the size, shape, and orientation of scatterers can provide unique information on the characteristics of scatterers that has been otherwise unattainable, although creating the means for utilizing these data in biological applications is still a developing research front \citep{Gauthreaux++2008, Bridge++2011}.  In meteorological applications,  polarimetric radar measurements are used in classification algorithms to characterize hydrometeor types \citep{Zrnic+1999, Park++2009, Chandrasekar++2012}.  The parallel for biological taxonomic classification has been widely speculated \citep{Bridge++2011}; however, research on this topic has been limited to a few studies \citep{Mueller+1985,Zrnic+1998,Bachmann+2007, vanDenBroeke2013, Melnikov++2014}. 

Weather radar measurements are typically interpreted in meteorological applications, and as a result, much of the common knowledge of these radar observations falls within this context.  When measuring signals from biological scatterers, however, many of these assumptions can be invalid, leading to incorrect interpretations.  The following sections describe some deviations from these common assumptions, and outline the biological interpretation of NEXRAD products. 

\section{Deviations from Meteorological Interpretations}
\label{interp:meteo}
\subsection{Departures from the Rayleigh Scattering Regime}
\label{interp:meteo:rayleigh}
In the case of biological scatter, the basic assumptions put forth by \cite{Mie1908} are never true; not surprisingly, very few animals are spherical.  Furthermore, it is typically rare that the Rayleigh conditions are good approximations for biological scatterers.  In some cases, insects can be small enough to roughly act as Rayleigh scatterers, especially at S-band, but most birds and bats are too large with respect to the radar wavelength.  Similarly, the aerial densities of insects can satisfy the volume-filling assumption, but birds and bats do not always fill the beam uniformly.  This is especially true at range gates close to the radar, where smaller resolution volumes may only capture one or two individuals.


\subsection{Axes of Rotational Symmetry}
\label{interp:meteo:symmetry}
One difficulty in characterizing polarimetric quantities for biological scatterers is their lack of rotational symmetry in azimuth \citep{Gauthreaux++2008}.  For example, a single value characterizes the differential reflectivity of a 1 mm water droplet with no additional qualifications.  The value for oblate rain drops can vary depending on the elevation look angle, but are still rotationally symmetric in azimuth.  Even highly irregular hydrometeors such as ice aggregates can be described independent from azimuth due to high number concentrations and random orientations.  With few exceptions, bulk meteorological scatterers are characterized without azimuthal view angle qualifications.  Conversely, one does not simply provide bulk polarimetric values for biological scatterers.  

Consider the common question,``What is a typical $Z_{DR}$ for a bird?"  The first challenge is that a wide variety of bird sizes and shapes exist (Fig.~\ref{interp:birds}).  Additionally, the body positions of birds can change drastically while in flight (e.g., wings out, wings in, soaring, climbing, diving, etc.), which can all affect the overall aspect ratio.  Even in the ideal case of densely-spaced, single-species ensembles, polarimetric quantities can span a large portion of their available dynamic range depending on the azimuthal orientation. This case is illustrated by \cite{vanDenBroeke2013}, in which the differential reflectivity of a Purple Martin (\emph{Progne subis}) colony ranges from $-4$ to $6$ dB and correlation coefficient spans $0.5$ through $1$.  Similar cases are presented in the following sections.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/interp_birds.pdf}
\caption{Birds of Oklahoma and their relative scales. (clockwise from upper-left) Ruby-throated Hummingbird, Great Blue Heron, Yellow Warbler, Northern Cardinal, Scissor-tailed Flycatcher, and Bald Eagle.}
\label{interp:birds}
  \end{figure}

\subsection{Significance of Depolarization}
\label{interp:meteo:depolar}
Another consequence of the non-spherical body shapes of animals is the non-negligible cross-polar scattering terms, $s_{vh}$ and $s_{hv}$.  As incident electromagnetic fields induce currents within and along an animal's body, the irregular shape enables cross coupling and depolarization \citep{Melnikov++2015}.  These contributions are almost always negligible in meteorological scatterers, with very few exceptions \citep{Ryzhkov+2007}.  The effects of depolarization are especially significant within the STAR configuration, when cross-polar terms contribute to the reported measurements \citep{Melnikov++2015}.  The result of cross-polar contributions to $Z_{DR}$ are presented in Section~\ref{interp:zdr}.

\subsection{Behavioral Considerations}
\label{interp:meteo:alignment}
The motion, orientation, and distribution of meteorological scatterers are typically dictated by kinematics and thermodynamics---most commonly, gravity and air motions \citep{DoviakZrnic}.  In some specialized cases, such as electrically charged storms, hydrometeors can assume more unique configurations \citep{Ryzhkov+2007}.  As sentient, volant creatures, biological scatterers have a wider range of configurations beyond those driven by simple kinematics.  The behaviorally-driven motions and configurations create signals uncommon in meteorological observations.


\section{Radar Reflectivity ($\eta$) and Reflectivity Factor ($Z$)}
\label{interp:z}
Much as in meteorological applications, radar reflectivity describes the combined effects of the size and density of scatterers.  Because organisms are non-spherical with no rotational symmetry around a vertical axis, their radar cross-sections are highly variable. Thus, reflectivity factors must be interpreted with respect to the anticipated animal orientations.  In some situations, animal orientations within the resolution volume can be treated as random, allowing an ensemble averaged radar cross-section to represent the scattering area of each organism.  In other cases, clear alignment of organisms exists, and cross-sections must account for orientation to yield correct interpretations.  Cases illustrating both situations are presented in Section~\ref{NEXRAD}.

Once an appropriate radar cross-section has been determined, the radar reflectivity factor  value can be related to the approximate number of organisms within the sampling volume \citep{Chilson++2012_2Z}.  By relating equations \ref{back:rayleigh} and \ref{back:z}, the value of $Z$ can be transformed into $\eta$ by adding a radar-dependent constant, $10\log_{10}(10000\pi^5|K|^2\lambda^{-4})$, with $\lambda$ being the radar wavelength in cm \citep{Chilson++2012_2Z}.  Using equation \ref{back:eta}, the reflectivity can be decomposed into the contributions from multiple scatterers of a given radar cross-section.  For some applications, the ensemble of scatterers within a volume will be a single species (see Section~\ref{NEXRAD:ex}).  In these cases, equation \ref{back:eta} is solved by simply dividing the total reflectivity by the expected radar cross-section of an individual.  In other cases, the resolution volume may contain several species of birds, bats, or insects.  When this is the case and organismal compositions can be estimated, the number of specific species can be solved for by considering the proportional contributions from each species.

For example, consider an experimental setup in which a vertically-pointing microphone records calling birds in flight within a radar sampling volume.  From these acoustic data, the relative proportion of birds can be estimated.  If only duck and sparrow calls are recorded, then the reflectivity contributions can be described as
\begin{equation} 
\eta = \eta_{duck}+\eta_{sparrow} = \frac{1}{V}\sum_{i=1}^{vol} \sigma_{b,i}^{duck} + \frac{1}{V}\sum_{i=1}^{vol} \sigma_{b,i}^{sparrow}\;\;\;\;\;\;\;\;\;\;\;\;  [\textrm{cm}^2~\textrm{km}^{-3}]
 \label{interp:eta_species}
\end{equation}
with superscripts indicating the radar cross-section of the specified bird type. The number of each bird type can be estimated as  
\begin{subequations}
\begin{align}
n_{duck} &= \textrm{floor}\left(\frac{V \eta_{duck}}{\sigma_b^{duck}}\right)\\
n_{sparrow} &= \textrm{floor}\left(\frac{V \eta_{sparrow}}{\sigma_b^{sparrow}}\right)
 \label{interp:eta_species2}
 \end{align}
\end{subequations}
with floor($\cdot$) indicating truncation to integer values. The decomposed values of $\eta_{duck}$ and $\eta_{sparrow}$ are not known, but acoustic recordings can provide proxies for relative animal abundance \citep{Farnsworth++2004, Horton++2015}.  If the recordings indicate a duck to sparrow ratio of 1:10, the total number of each can be solved for by the addition of the relation,
\begin{equation} 
n_{duck} = \frac{1}{10}n_{sparrow}
 \label{interp:eta_species}
\end{equation}
providing an estimate of the total number of each in the radar sample volume.  An acoustic recording and processing method is presented in Chapter~\ref{local} that can provide species-specific information of calling animals in the airspace. 

Another major factor when considering the reflectivity of multiple species ensembles is their relative size distribution.  In most applications, the relevant biological measures are absolute number, number concentration, or biomass---all of which are generally dominated by the smaller, more prevalent organisms.  While the small organisms, e.g., insects, often contribute most to these metrics, radar measurements are strongly biased toward larger body sizes.  As a simple example, consider a 1-mm gnat and a 1-cm housefly.  Under the Rayleigh assumption, reflectivity is proportional to the product of the insect diameter to the sixth power and the number concentration (equation \ref{back:z}). As a result, 
\begin{subequations}
\begin{align}
z_{gnat} &= N_{c,gnat}(10^{-3})^6 = N_{c,gnat}10^{-18}\\
z_{fly} &= N_{c,fly}(10^{-2})^6 = N_{c,fly}10^{-12}
 \label{interp:gnat}
 \end{align}
\end{subequations}
with $N_c$ denoting the insect number concentration.  Based on the diameter component of reflectivity factor, the housefly contribution is six orders of magnitude greater than the gnat; or in other words, the reflectivity contribution from size is 1,000,000 times greater for the fly.  Similarly, when a large size discrepancy exists between populations of scatterers, reflectivity-weighted measurements will be heavily biased toward larger organisms.  In this scenario, gnats could only have a comparable reflectivity contribution to houseflies if their number concentration exceed that of houseflies by a factor of one million.  This analogy can be extended further by considering a volume of 1-cm houseflies and 10-cm warblers, although the effects of resonance may lessen the disparity. In every case, the only way smaller organisms will have a significant contribution to the final signal is if their number concentrations are many orders of magnitude higher than the larger organisms.  In the extreme case of organisms on the order of a meter (e.g., cranes, raptors, waterfowl, etc.), it may be unlikely that any other organisms could reasonably produce notable signal contributions.   

\section{Mean Radial Velocity ($v_r$) and Spectrum Width ($\sigma_v$)}
\label{interp:v}
A number of motion contributions can be recorded within the Doppler spectrum, including the movement of organisms as well as their periodic body motions, e.g., wing flapping \citep{Bruderer1997a}.  For typical spectral processing, the reduction to a Gaussian fit only retains the average movement vector of the sampled organisms and a measure of the velocity diversity along the beam direction within a sampling volume. Because these two measures are only sensitive to motions along the beam radial, they are highly dependent on the pointing direction of the beam. For example, a beam pointed vertically at a flock of migrating geese only resolves the velocity from changes in altitude, not horizontal motions.

In addition to the purposeful movements of animals, atmospheric motions also contribute to the velocity of volant animals \citep{Green+2002, Chapman++2011}.  From a reference frame fixed at a point on the ground, an insect with negligible flight speed will be embedded within the bulk flow and move as a passive tracer of the surrounding wind field \citep{Lhermitte1966}. For faster flying insects, birds, and bats, this velocity with respect to the ground is the vector sum of the flight velocity of the organism and the wind velocity, forming the so-called velocity triangle \citep{Chapman++2011}.  Within this framework, the animal flight vector defines the \emph{air speed} and \emph{heading} while the resultant is the \emph{ground speed} and \emph{track} \citep{Larkin1980,Green+2002}. Figure~\ref{interp:triangles} demonstrates the vector contributions to overall movement for (a) passive drift, (b) flight along the wind direction, and (c) flight across the wind direction.  In all cases, the radial velocity would be the radial component of the resultant (i.e., ground speed) vector.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/interp_triangles.pdf}
\caption{Velocity vector contributions for passive drift, flight along the wind, and flight across the wind.  Dotted red arrows indicate wind vectors, dotted blue arrows indicate animal flight vectors (air speeds), and black arrows indicate the resultant (ground speed).}
\label{interp:triangles}
  \end{figure}

When multiple organisms are captured within a single resolution volume, the corresponding Doppler spectrum will be multimodal with peaks associated with the radial component of each individual.  If individuals are traveling in highly variable directions, the Gaussian fit will reduce these details to a broad distribution with a mean velocity of approximately zero.  In such cases, the radial velocity will be indicative of the prevailing winds. Following the discussion on gnats and flies in the previous section, the Doppler spectrum will be strongly weighted toward larger organisms. In a scenario with similar concentrations of 1-cm insects flying southward and 10-cm birds flying northward, it is likely that the Gaussian fit of the Doppler spectrum would indicate exclusively northward motions.  In this respect, the topic of bird-insect delineation is likely a problem of determining bird presence or absence.

Presently, knowledge of the speed of the organisms of interest can be used to eliminate signals from other organisms, and is commonly used for bird-insect delineation \citep{Larkin1991, Gauthreaux+2003, Liu++2005, Zhang++2005, Cruz++2013}. Additionally, the diverse flapping of collections of organisms often results in enhanced texture in the velocity product across adjacent resolution bins, yielding a possible indicator of biological scatter \citep{Dokter++2011}.  Furthermore, these variations in velocities lead to higher values of the spectrum width of radial velocity in Doppler radars \citep{Koistinen2000, vanGasteren++2008, Holleman++2008}.

\section{Differential Reflectivity ($Z_{DR}$)}
\label{interp:zdr}
Differential reflectivity is perhaps the most explored polarimetric quantity in biological applications \citep{Mueller+1985,Zrnic+1998,Melnikov++2012}, although its physical interpretation can be quite complicated.  Unlike meteorological scatterers, most organisms fall in the resonance scattering regime, which results in oscillating radar cross-sections with slight changes in size \citep{Bruderer1997a}.  As a result, the typical connection of $Z_{DR}$ to physical aspect ratio is generally invalid, with the possible exception of small insects.  

Within the ATAR framework, $Z_{DR}$ can be simplified to 
\begin{equation} 
Z_{DR} = 10\textrm{log}_{10}\left(\left<\frac{|s_{hh}|^2}{|s_{vv}|^2}\right>\right)
 \label{interp:zdrATAR}
\end{equation}
with $<\cdot>$ indicating the expected value \citep{DoviakZrnic}.  For STAR configuration, orthogonal polarizations are transmitted with uncontrolled initial phases and an arbitrary differential transmit phase, $\psi_t$.  Upon scattering, components from both polarizations are depolarized into the orthogonal channels, with contributions summing coherently.  That is, the signal received on the horizontal channel will be the scattered co-polar horizontal wave combined with the incident vertical wave that has been depolarized into the horizontal.  The way that these contributions combine---constructively, destructively, or otherwise---is determined by their phase differential at the receiver.  In the absence of propagation effects (a good assumption for most biological, clear air scatter), this differential phase will be a result of the differential phase upon scattering (during depolarization), and equally, the system differential phase on transmission \citep{Melnikov++2015}.

As a consequence of the high depolarization of biological scatterers, values of $Z_{DR}$ are just as dependent on the arbitrary transmit phase differential as the scatterers themselves.  Furthermore, changes in the differential transmit phase will result in different $Z_{DR}$ values.  In summary, if a radar engineer inadvertently hits his head on a low-hanging segment of transmission waveguide, the $Z_{DR}$ of a Great Blue Heron will instantly change.   

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/interp_zdr.pdf}
\caption{View from above of scattered waves from a bird oriented (a) northward and (b) southward. Blue vectors indicate the scattered co-polar horizontal wave component ($S_{HH}$) and the red vectors indicate the scattered cross-polar horizontal wave component ($S_{HV}$).}
\label{interp:zdrcross}
  \end{figure}

A second consequence of this depolarization process is the lack of axisymmetric scatter for individuals and aligned organisms \citep{Melnikov++2015}.  During the depolarization process, surface currents are excited along the organism's body.  The body will act  as an antenna, radiating the scattered waves outward.  If the phase of the depolarized wave is determined by the animal anatomy, $Z_{DR}$ will be different when looking at the left side than the right side, for example.  Consider the idealized schematic plan configuration in Fig~\ref{interp:zdrcross}a in which a radar located at the left of the image emits a dual polarized wave toward a bird.  The blue vectors indicate the scattered co-polar horizontal wave component.  Imagine the incident vertically polarized wave component reaches the bird and is depolarized such that the scattered wave has an initial phase toward the bird's head (red).  The resulting received horizontal wave would be the coherent wave addition with constructive interference.  Now consider the exact same configuration except the bird is oriented southward (Fig~\ref{interp:zdrcross}b).  While the incident waves are identical in every way, the depolarization of the vertical wave component will still excite surface currents with respect to anatomy, resulting in a cross-polar scattered initial phase again directed headward.  In this case, the waves will add destructively.  The end result is that the $Z_{DR}$ of bird will be different when viewing from one side and the other.  More formally, 
\begin{subequations}
\begin{align} 
\arg(s_{hv}(\vartheta)) &= -\arg(s_{hv}(-\vartheta))\\
\arg(s_{vh}(\vartheta)) &= -\arg(s_{vh}(-\vartheta))\\
\arg(s_{hv}(\vartheta)) &= \arg(s_{vh}(\vartheta))\\
 \label{interp:phaseFlop}
 \end{align}
\end{subequations}
with $\vartheta$ denoting the azimuthal look angle and $-\vartheta$ denoting the axisymmetric azimuth.  Currently this process is speculation, but the observations presented in Section~\ref{NEXRAD} demonstrate an absence of axes of symmetry in $Z_{DR}$, supporting some orientation-dependent variation.  


\section{Differential Phase Shifts ($\psi_{DP}$ and $\delta$)}
\label{interp:phidp}
The phase measurements provided by NEXRAD can take on different interpretations depending on the propagation medium and scattering bodies.  From \cite{Zrnic+1998}, the measured differential phase ($\phi_{DP}$) is the sum of the contributions from propagation through a medium ($\Phi_{DP}$) and the differential scattering phase ($\delta$).  While the former study presents the measured phase shift as a proxy for animal $\delta$ values,  \cite{Ryzhkov+2007} and \cite{Melnikov++2015} show that an additional parameter, the system differential phase ($\psi_{sys}$), contributes to the measured differential phase and varies among NEXRAD sites.  Furthermore, this system differential phase can be decomposed into two components: the differential phase upon transmission and the differential phase shift in reception.  Considering system phase effects motivates a more general form of measured differential phase, hereafter denoted as
\begin{equation} 
\psi_{DP} = \psi_t+\frac{1}{2}\Phi_{DP}+\delta+\frac{1}{2}\Phi_{DP}+\psi_r
 \label{interp:psi}
\end{equation}
with
\begin{equation} 
\psi_{sys} = \psi_t+\psi_r
 \label{interp:psiSys}
\end{equation}
and the standard definition of 
\begin{equation} 
\phi_{DP} = \Phi_{DP}+\delta
 \label{interp:phiDP}
\end{equation}
in which $\psi_{DP}$ is the value recorded and reported by the radar \citep{Melnikov++2015}.

Values of $\psi_{sys}$ can be obtained by observing $\psi_{DP}$ on the outer edges of rain showers in otherwise clear air.  In such cases, there is no propagation differential phase along the clear-air path from the radar to the edge of the storm.  Additionally, the rain does not contribute any differential phase upon scattering.  Therefore, the differential phase recorded is the result of phase offsets introduced during the combination of transmission and reception, although these two effects cannot be separated.  Once the value of $\psi_{sys}$ has been obtained for a given radar, it can be subtracted from the $\psi_{DP}$ measurements to yield true $\phi_{DP}$ values.  It should again be noted that the stability of $\psi_{DP}$ in time is a largely unexplored topic, and it is not clear how long the validity of $\psi_{sys}$ values remain, or to what extent they vary with diurnal environmental variability (e.g., temperature of waveguide).

In many meteorological applications, differential scatter phase is negligible, and the range derivative of the differential phase, $K_{DP}$, is calculated for practical uses \citep{DoviakZrnic}.  Because $\psi_{DP}$ introduces a constant phase offset, it does not affect $K_{DP}$ values, and is likely overlooked by many radar meteorologists.  Similarly, the negligible depolarization of most meteorological scatterers inhibits the phase-controlled cross-polar effects on $Z_{DR}$, again resulting in little impact on meteorological applications.


\section{Co-polar Correlation Coefficient ($\rho_{HV}$)}
\label{interp:rhohv}
The cross-correlation coefficient between co-polar channels ($\rho_{HV}$) is often used for discriminating between meteorological and biological signals \citep{Park++2009,vanDenBroeke2013}.  While the magnitude of $\rho_{HV}$ from biological scatter can regularly overlap that of meteorological echoes, it does provide some consistent morphological features. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Radio Scattering Measurements of Biota}
\label{measure}
\section{Background}
\label{measure:background}
The radar variables described in the previous chapter are the expected values of ensembles of individuals.  To interpret these data at a taxonomic level, radial scattering measurements of individuals need to be performed.  Several studies have measured the effect of azimuth look angle on radar cross section \citep{Edwards+1959, Blacksmith+1965}, but few have made measurements at orthogonal polarizations.  The goal of this study was to develop a methodology of making precision radio measurements of dead animal subjects.  Two major challenges complicated the radar cross-section measurements of an animal carcass.  First, the size and dielectric composition of the subjects yield relatively weak echo signals that are easily overcome by surrounding clutter or environmental interference.  Second, the bodies of dead birds and bats are quite frail and flaccid, making exact and consistent positioning difficult.  Coupled with the low radar cross-section, many typical methods for supporting the body (e.g., foam pedestals and supports) were simply too obstructive for this application.  The following describes the hardware implementation used for these measurements, and discusses our methods for overcoming these technical challenges.

A number of initial efforts to conduct radar scattering measurements of live and dead animal specimens were attempted using a scatterometer (detailed in \citealt{Kong2014}) inside an anechoic chamber (Fig.~\ref{measure:fail}).  In these cases, small echoes from imperfections in the chamber---such as exposed power outlets, lights, and extension cords---contaminated the measurements.  Later attempts to isolate the animal signals were made by swinging a suspended animal to generate a non-zero Doppler frequency, and removing zero-Doppler clutter.  It is unclear to what extent the radial motion contributed to variation in the boresite position of the animal, and how that affected the quality of the measurements.  Furthermore, when attempting measurements of different body positions, it is possible that the boresite location changed.   

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/measure_fail.JPG}
\caption{A previous experimental setup for measuring the radar cross-section of a Brazilian free-tailed bat in an anechoic chamber.}
\label{measure:fail}
  \end{figure}

To avoid these complications, a new method for making precision measurements was developed.  The major methodological changes were threefold:
\begin{enumerate}
\item A network analyzer was used to make dual-polarization measurements across a large frequency sweep.  
\item The animals were suspended by an apparatus that allowed azimuthal rotation while keeping the center of mass within the antenna boresite.
\item Measurements were taken outdoors in a skyward configuration to avoid clutter and second-trip echoes.
\end{enumerate}

Network analyzers provide the ability to measure the pattern of antennas---a process technically similar to radar cross-section measurements.  Additionally, these devices are small, lightweight, and portable enough to enable convenient onsite measurements.  The network analyzer used for these measurements was the Agilent E8364B, a general-purpose network analyzer equipped with four receiver ports.  To implement direct polarimetric measurement, we decoupled the instrument ports as proposed by \cite{Agilent2004}, using them to sample the four polarization return channels from the scatterer.  In this configuration, the received signal components are exactly to co-polar and cross-polar contributions of the \textbf{S} matrix from equation \ref{back:smatrix}.  To ensure reliable measurements, the transmit and receive antennas require precise control of the cross-polarization coupling and phase center position. The antennas used are dual-polarized, ETS Lindgren quad-ridged horns (model 3164-05).  Confirmed by our laboratory measurements, these antennas provide over 24 dB of isolation between orthogonal polarizations and more than 25 dB of backfire suppression \citep{Lindgren2005}.  

Obtaining precise and reproducible body positions when measuring such a small, limp, and fragile subject is a challenging task.  It is also imperative that the apparatus suspending the animal does not interact with the measurement, either as obtrusive clutter or a source of electromagnetic coupling.  Many previous measurement efforts have used polyfoam supports or enclosures, citing their negligible electromagnetic contributions \citep{Blacksmith+1965}.  While such a rig may be suitable for larger (in the RCS sense) objects, our experience has shown that the weak scattering characteristics of birds and bats are often of the same magnitude as otherwise negligible clutter.  

To eliminate these clutter sources, we opted to perform outdoor sky-looking measurements and constructed a planar positioning apparatus that minimized interference (Fig.~\ref{measure:cal}). The planar positioner consists of two identical wooden structures separated by four meters, with each structure being made up of a rotating disk atop a 2.3-meter stationary post.  A 360-degree protractor is fastened to each disk, providing rotations in 1-degree increments. Segments of 0.25 mm diameter low-permittivity monofilament line spanned between the disks, suspending the scatterer in the antenna boresite.  The antenna was placed at ground level facing upward with cables running to the network analyzer several meters away.  The network analyzer was then connected to the control computer even farther away to minimize any sidelobe clutter from the operators (Fig.~\ref{measure:cal}). 

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/measure_cal.pdf}
\caption{Experimental setup for calibration using metal sphere. Annotated close-ups of components are provided in Fig.~\ref{measure:batSetup}}
\label{measure:cal}
  \end{figure}
  
Post processing of the observations allowed us to remove unwanted artifacts in the data, thereby improving our RCS measurements. For example, the total returned signal during the measurement consists of reflections from the scatterer under investigation, along with reflections from other unwanted background signals. The background signals were characterized by conducting measurements in the absence of the object being investigated. In this way the noise floor can be set by the averaged value of the time dependent signals from the background (here, -103dBm).  Moreover, to mitigate statistical errors in amplitude measurements and reduce noise in the samples, 10 consecutive samples of the received signal were averaged in time. The temporal averaging was implemented when characterizing both the background and the scattering properties of the bat. Signal post-processing using time gating was implemented to further improve the quality of the RCS measurements by better suppressing second trip reflections, especially those from the ground.
     
Various calibration techniques exist for radar cross-section measurements, but all rely on measurements of objects having well-determined RCS characteristics, such as spheres, cylinders or plates composed of conducting material. Here we used a conductive sphere with a diameter of 15.82 mm (Fig.~\ref{measure:cal}).  The size of the sphere was chosen to be comparable with that of a typical bird or bat. The true (calibration) RCS for a sphere was calculated using a precision computational electromagnetics tools available in WIPL-D \citep{Kolundzija+2002}. By performing the calculation of the calibration standard using the WIPL-D software we removed the possibility of inherited errors in approximate methods, which could affect the quality of calibration. 



\section{Measurements of a Brazilian Free-tailed Bat (\emph{Tadarida brasiliensis})}
\label{measure:tadarida}
The first specimen for measurement was a Brazilian free-tailed bat carcass that had been collected (IACUC R09-29), frozen for temporary storage, and thawed prior to measurements.  Four 0.25 mm diameter low-permittivity monofilament lines spanned between the two counterfacing disks, piercing the bat in several places, and suspending it in the middle of the rig within the antenna boresite. Upon applying slight tension to all lines, the bat was suspended in a rigid, flight-like body position above the antennas (Fig.~\ref{measure:batSetup}).  Initially, the bat faced upward such that the antenna recorded at tail-on view.  The measurement process was conducted in two steps, repeated 360 times.
\begin{figure}[h!]
\centering
\includegraphics[width=6in]
  {figures/measure_batSetup.pdf}
  \caption{Laboratory equipment setup for field measurements. (a) The control area located several meters away from the antenna. (b) Suspended Brazilian free-tailed bat. (c) Network analyzer with dual transmit and dual receive channels. (d) Quad-ridged dual-polarization antennas.}
\label{measure:batSetup}
\end{figure}
The network analyzer acquired 10 samples of the measurement, swept over a range of frequencies. Each sample was made up of the transmission of a single-polarization wave and the subsequent reception of both co-polar and cross-polar components, followed by the transmission and reception of the orthogonal polarization (i.e., ATAR with simultaneous co- and cross-polar reception).  Following the measurement, the disks were rotated by one degree, resulting in a one-degree change in azimuth.  Measurements for the new azimuth were acquired, the disks were rotated, and the process was repeated.  Care was taken to ensure that the measurements were taken while the bat was completely stationary. Any movement resulting from wind gusts was assumed to corrupt the measurement, and the data acquisition was repeated for that azimuth.  Overall, once the setup and calibration was complete, the 360-degree measurements took around 4 hours to finish.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/measure_RCS.pdf}
\caption{Measured radar cross-sections for a Brazilian free-tailed bat for (a) co-polar and (b) cross-polar components.}
\label{measure:RCS}
  \end{figure}
  
Following post processing and calibration, the result was the S-parameter matrix as a function of frequency and azimuth.  Theoretically, the sequestration of polarizations by the network analyzer should result in axisymmetric radar cross-sections across the bat body (as with ATAR radar operation).  Furthermore, measurement cross-polar terms ($s_{hv}$ and $s_{vh}$) should be equal.  With this in mind, the side of the bat with the most similar cross-polar terms was taken as the highest-quality measurement, and reflected across the body axis.  From these final data, the radar cross-sections were computed using equations 2.3(a-d).  The resulting polar plots of radar cross-sections at S-band are shown in Fig.~\ref{measure:RCS}.  Additionally, the corresponding differential radar cross-section of the co-polar channels was computed (Fig.~\ref{measure:DRCS}).

  \begin{figure}[t]
\centering
\includegraphics[width=3.5in]
  {figures/measure_DRCS.pdf}
\caption{Measured differential radar cross-sections for a Brazilian free-tailed bat at X-band.}
\label{measure:DRCS}
  \end{figure}
  
 Comparing the magnitude of the co-polar scattering components, it is clear that neither of the terms are consistently larger than the other, resulting in both positive and negative differential radar-cross-sections.  In fact, all four scattering magnitudes are of similar size (around -40 dBsm).  This measurement supports the significance of the cross-polar contributions to STAR radar signals.  Furthermore, because the interference of the cross-polar contributions is determined by the transmission phases, these components can effectively double or negate the co-polar signals, depending totally on the phase.  This effect can be demonstrated by using the measured \textbf{S} matrix to add the co-polar and  cross-polar contributions with varying transmit phase offsets to find the final STAR radar cross-sections,
\begin{subequations} 
\begin{align}
\sigma_{b,H} &= 4\pi|s_{hh}+s_{hv}\exp(\j\psi_t)|^2\\
\sigma_{b,V} &= 4\pi|s_{vv}+s_{vh}\exp(-\j\psi_t)|^2
 \label{measure:star}
 \end{align}
\end{subequations}
in which positive values of $\psi_t$ are defined as horizontal transmit phases preceding the vertical.  For the first calculation, the symmetric scatter phases were used, resulting in symmetric differential radar cross-sections (Fig.~\ref{measure:DRCSstarSym}).  From these results, the difference between radar cross-sections can range from -48.3 to 49.7 dB with azimuthal morphology depending on the transmit differential phase.


  \begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/measure_DRCSstarSym.pdf}
\caption{X-band differential radar cross-sections as a function of transmit phase offset for a Brazilian free-tailed bat using symmetric cross-polar scatter phases.}
\label{measure:DRCSstarSym}
  \end{figure}

  \begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/measure_DRCSstarAnti.pdf}
\caption{X-band differential radar cross-sections as a function of transmit phase offset for a Brazilian free-tailed bat using anti-symmetric cross-polar scatter phases.}
\label{measure:DRCSstarAnti}
  \end{figure}

For the second calculation, anti-symmetric cross-polar scatter phases were used (as in equations 3.5a-d), resulting in non-symmetric differential radar cross-sections (Fig.~\ref{measure:DRCSstarAnti}).  From these results, the difference between radar cross-sections can range from -46.3 to 56.1 dB, again with azimuthal morphology depending on the transmit differential phase.  In both cases, whether symmetric or not, the magnitude and sign of the differential radar cross-section at a fixed view azimuth is highly variable and dependent on the transmit phase differential. 


\section{Measurements of a Brown-headed Cowbird (\emph{Molothrus ater})}
\label{measure:BHCO}
The second specimen was a Brown-headed Cowbird carcass that had been collected (IACUC R09-29), frozen for temporary storage, and thawed prior to measurements.  Colloquially, it is common knowledge that the contributions of a bird's wings are negligible to it's overall radar-cross-section.  However, this oft-quoted fact has never been confirmed by supporting measurements.  To explore this subject, the process used to measure the bat was repeated for the cowbird.  To expedite the measurements, azimuth increments were changed from 1-deg to 5-deg, resulting in slightly coarser resolution.  Additionally, the 360-degree measurements were repeated for three variations: First the bird was suspend with wings stretched outward; second, the wings were tied against the bird's body; third, the wings were removed from the carcass (Fig.~\ref{measure:birdSetup}).  In each case, great care was taken to keep the body in a fixed position while only manipulating the wings.  


\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/measure_birdSetup.pdf}
  \caption{Brown-headed Cowbird configuration for measurements with (a) wings outstretched, (b) wings tucked, and (c) wings removed.}
\label{measure:birdSetup}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/measure_bird.pdf}
  \caption{Brown-headed Cowbird RCS measurements with (a) wings outstretched, (b) wings tucked, and (c) wings removed. Diagrams are on a fixed, relative decibel scale.}
\label{measure:bird}
\end{figure}


Overall, there is a noticeable difference in the X- band morphology of the cross-section in azimuth between the measurements taken with wings out and wings in (Fig.~\ref{measure:bird}a, b).  This would support that wing flapping patterns are observable from modulations in radar amplitude returns.  Considering the great similarity of the measurements with wings removed to wings in (Fig.~\ref{measure:bird}b,c), it is clear that wings themselves are not contributing to the wing-flapping signature.  Instead, it is possible that the repetitive flexing and relaxing motions of the chest muscles are contributing to these changes, resulting in the flapping signatures.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Observation and Interpretation of Biological Signals in NEXRAD}
\label{NEXRAD}

\section{Background}
\label{NEXRAD:back}
The NEXRAD network provides continuous, continental scale surveillance of the airspace over the United States and the animals within \citep{Chilson++2012}.  These data are quality controlled, archived, downloadable, and freely available over the internet, making them ideal for broad scale or long range aeroecological studies.  The weather surveillance radars in the NEXRAD network operate at S-band and have Doppler and polarimetric capabilities. 



\section{Diurnal Biological Signals}
\label{NEXRAD:day}

The signals typical of the daytime hours are dominated by features of the convective boundary layer.  Organisms in the airspace include insects that become trapped in thermals as well as soaring birds such as hawks, eagles, kites, and swallows \citep{Russell+1997}. Depending on the time of year, these signals may have no organized movements, apart from those stemming from thermally-driven circulations, or may contain daytime migrants such as raptors \citep{Diehl2013}.  Figure~\ref{NEXRAD:DFXday} shows the six NEXRAD products for a typical summer day (03 Aug. 2014, 19:05 UTC) in Laughlin, Texas (KDFX).  Reflectivity values are generally low (5-15 dBZ), velocity is weak (approximately 5 m s$^{-1}$), and $Z_{DR}$ is positive throughout the region (4-6 dB).  The differential phase has been calibrated to remove the system differential phase contributions (i.e., conversion from $\psi_{DP}$ to $\phi_{DP}$), as evidenced by the fact that the leading edge of precipitation has a value of zero, and $\rho_{HV}$ is generally low with a slight difference in magnitude longitudinally across the domain.  Overall, the combination of these patterns indicates a boundary layer dominated by thermally-lofted insects that are very slightly preferentially oriented.

The radar products for the same date and time to the north-east in Austin, Texas (KEWX) are shown in Fig.~\ref{NEXRAD:EWXday}.  Overall the patterns are quite similar, but the velocity and spectrum width indicate greater boundary layer development, including horizontal convective rolls running north to south.  Additionally, both $\rho_{HV}$ and $\phi_{DP}$ reveal greater azimuthal asymmetry across the region, indicating a higher degree of preferential orientation.

\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_dayDFX.pdf}
  \caption{NEXRAD products from Laughlin AFB, Texas (KDFX) on 03 August 2014. Products are (clockwise from upper left) radar reflectivity factor ($Z$, [dBZ]), radial velocity ($v_r$, [m s$^{-1}$]), spectrum width ($\sigma_v$, [m s$^{-1}$]), differential phase ($\phi_{DP}$, [deg]), co-polar correlation coefficient ($\rho_{HV}$, [~~]), and differential reflectivity ($Z_{DR}$, [dB]).}
  \label{NEXRAD:DFXday}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_dayEWX.pdf}
  \caption{NEXRAD products from Austin, Texas (KEWX) on 03 August 2014. Products are (clockwise from upper left) radar reflectivity factor ($Z$, [dBZ]), radial velocity ($v_r$, [m s$^{-1}$]), spectrum width ($\sigma_v$, [m s$^{-1}$]), differential phase ($\phi_{DP}$, [deg]), co-polar correlation coefficient ($\rho_{HV}$, [~~]), and differential reflectivity ($Z_{DR}$, [dB]).}
  \label{NEXRAD:EWXday}
\end{figure}



\section{Roost Exidus}
\label{NEXRAD:ex}
\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_cave.jpg}
  \caption{A bat emergence at Frio Cave in Rio Frio, Texas (23 June 2011).}
  \label{NEXRAD:cave}
\end{figure}
\subsection{Brazilian Free-tailed Bat (\emph{Tadarida brasiliensis})}
\label{NEXRAD:ex:tadarida}
In some cases, large aggregations of single-species communities form identifiable signatures on radar.  Such is the case for the Brazilian free-tailed bats (\emph{Tadarida brasiliensis}) that inhabit caves and bridges in the south central plains.  These bats roost communally, and immediately following sunset, exit the roost in a mass exodus that can be observed routinely on radar \citep{Frick++2012}.  A photograph of such an emergence at Frio Cave is shown in Fig. \ref{NEXRAD:cave}.  To decrease competition for food among individuals, the bats spread throughout the region to forage for insects.  The result of this behavior with respect to radar sensing is twofold.  First, the emergence signatures at these known cave sites are predominantly a single, known species.  Second, because of the divergence behavior, the radar views ensembles of bats from most or all azimuth look angles. That is, some bats are flying toward the radar, yielding a head-on view; some are flying away, giving a tail-on view; and others are flying in each perpendicular direction.  

\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_DFX_all.pdf}
  \caption{NEXRAD products from Laughlin AFB, Texas (KDFX) on 03 August 2014. Products are (clockwise from upper left) radar reflectivity factor ($Z$, [dBZ]), radial velocity ($v_r$, [m s$^{-1}$]), spectrum width ($\sigma_v$, [m s$^{-1}$]), differential phase ($\psi_{DP}$, [deg]), co-polar correlation coefficient ($\rho_{HV}$, [~~]), and differential reflectivity ($Z_{DR}$, [dB]).}
  \label{NEXRAD:DFX}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_EWX_all.pdf}
  \caption{NEXRAD products from Austin, Texas (KEWX) on 03 August 2014. Products are (clockwise from upper left) radar reflectivity factor ($Z$, [dBZ]), radial velocity ($v_r$, [m s$^{-1}$]), spectrum width ($\sigma_v$, [m s$^{-1}$]), differential phase ($\psi_{DP}$, [deg]), co-polar correlation coefficient ($\rho_{HV}$, [~~]), and differential reflectivity ($Z_{DR}$, [dB]).}
  \label{NEXRAD:EWX}
\end{figure}

Figure~\ref{NEXRAD:DFX} shows a typical evening bat emergence as seen by KDFX.  Two known bat caves have been annotated in the reflectivity image.  The divergence signatures can be identified on radial velocity by a transition from inward to outward velocities along a radial.  For comparison, Figure~\ref{NEXRAD:EWX} shows the same emergence as seen from the Austin, Texas radar.  In both figures, Ney Cave and a second unidentified colony to the south are visible.  From both radar perspectives, the emergences surrounding Ney Cave have highly variable quantities in all of the radar products.  In these cases, spectrum width shows high velocity variance in the head-on and tail-on measurements with low values at cross-beam flight directions.   The physical interpretation is that the divergent velocity (i.e., the component directed radially away from the cave) is highly variable in magnitude, but quite consistent in direction.  For example, if the individual bat motions were predominantly isotropic in the horizontal with a slight divergent component, the spectrum width would be approximately equal around the circumference of the divergent ring.  The lower values at cross-beam flight angles means that the bats are actually flying perpendicular to the beam with little variance in the radial direction.  In total, the behavior must be that the bats are flying strictly away from the cave at variable speeds. 

To compare the differential phases between the two radars, it is necessary to remove the system specific differential phase component.  Figure \ref{NEXRAD:DFXcal} demonstrates this process for KDFX.  It is first necessary to find a case in which precipitation is within the range of the radar, but with a clear-air path in between.  Starting with the raw differential phase measurement ($\psi_{DP}$), the values along the leading edge of the precipitation must be isolated (Fig. \ref{NEXRAD:DFXcal}a).  Because these volumes contain only precipitation, their $\delta$ values are zero.  Furthermore, the clear air path to and away from this volume will accumulate no phase differential from propagation.  As a result, the values of differential phase in these volumes are due only to the combined effects of the system on transmission and reception.  After finding this system value, it can be subtracted from $\psi_{DP}$ to yield $\phi_{DP}$  (Fig. \ref{NEXRAD:DFXcal}b).  Finally, because the typical density of airborne organisms does not contribute significant propagation phase shifts \citep{Zrnic+1998}, these values are representative of the phase shifts due to scattering.  With this in mind, the values above 270$^\circ$ are likely negative phase shifts, so a constant $360^\circ$ is subtracted to yield negative values (Fig. \ref{NEXRAD:DFXcal}c). This process was repeated for KWEX, resulting in $\delta$ measurements of the same bat emergence from two perspectives (Fig. \ref{NEXRAD:Ney}). 


\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_DFX_cal.pdf}
  \caption{Differential phase calibration for KDFX.}
  \label{NEXRAD:DFXcal}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_Ney.pdf}
  \caption{Simultaneous $\delta$ measurements of Ney Cave from KDFX and KEWX.}
  \label{NEXRAD:Ney}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_EWX_P.pdf}
  \caption{Detail of the radial leading edges of $\phi_{DP}$ in the emergence signatures.}
  \label{NEXRAD:PEWX}
\end{figure}

\clearpage
One difficulty in differential phase calibration is the ambiguity of the direction of phase shifts.  The decision to use the range (-90$^\circ$, 270$^\circ$) was determined subjectively, and it is quite difficult to determine whether phase shifts of $180^\circ$ represent forward or backward $\delta$ values.  It appears that a discontinuity exists in the phase measurements of KEWX, with $\delta$ values jumping from +50$^\circ$ to +250$^\circ$ over a short segment of the roost ring (Fig. \ref{NEXRAD:Ney}).  It is possible that these values should be -110$^\circ$, indicating a transition from positive to negative $\delta$ values.   Figure \ref{NEXRAD:PEWX} details the emergence differential phase signatures for KEWX, showing the point at which the ray impinges on the leading edge of each emergence.  The points are marked by $\delta$ values around 0-50$^\circ$, which do decrease outward hitting a discontinuity at the jump from -90$^\circ$ to 270$^\circ$.  Assuming that the red pixels should, in fact, be negative values exceeding -90$^\circ$, the differential phase decreases continuously to a minimum value of nearly -270$^\circ$.

 To help with interpretation, the differential reflectivity and correlation coefficient from both radars can be considered (Fig.~\ref{NEXRAD:NeyZDR}, Fig.~\ref{NEXRAD:NeyCC}).  Unlike the X-band measurements in the previous chapter, the NEXRAD S-band $Z_{DR}$ observations are generally positive, likely implying that the bats are \emph{almost} acting as rayleigh scatterers.  For most of the emergence signatures, areas of lower $Z_{DR}$ appear to be coincident with higher positive magnitudes of $\phi_{DP}$.


\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_NeyZDR.pdf}
  \caption{Simultaneous $Z_{DR}$ measurements of Ney Cave from KDFX and KEWX.}
  \label{NEXRAD:NeyZDR}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_NeyCC.pdf}
  \caption{Simultaneous $\rho_{HV}$ measurements of Ney Cave from KDFX and KEWX.}
  \label{NEXRAD:NeyCC}
\end{figure}

\clearpage

\subsection{Purple Martin (\emph{Progne subis})}
\label{NEXRAD:ex:PUMA}
The behavioral drivers of Purple Martins (\emph{Progne subis}) are quite similar to those of bats. Both roost communally, taking part in mass exoduses to forage for insects.  But while bats emerge at night, Purple Martins do so in the morning.  Purple Martins have been a popular candidate for radar studies for their unambiguous signatures as they emerge into the airspace \citep{Gauthreaux+1998, Russell+1998, Kelly++2012, vanDenBroeke2013}.  As insectivores, they also play an important role in regulating populations of aerial insects \citep{Kelly++2013}.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_MartinsOKC.JPG}
\caption{Purple Martin roost in Oklahoma City, Oklahoma in August 2013.}
\label{NEXRAD:OKCPUMA}
  \end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_PUMA.pdf}
  \caption{A Purple Martin exodus as seen by three local weather radars. The radar and roost locations are annotated on the map.}
\label{NEXRAD:PUMA}
\end{figure}
The first step toward characterizing radar polarimetric signatures of a phenomenon is physically understanding the phenomenon.  To investigate the behavior of Purple Martins at their late summer roosts, a mobile radar was deployed at the site of a roost in Garland, Texas in August 2011.  This roost emergence is regularly observed by the local NEXRAD site (KDFW) as well as two Terminal Doppler Weather Radars (TDWRs, TDAL and TDFW).  The simultaneous emergent signatures in reflectivity for the three radars (Fig.~\ref{NEXRAD:rax}) reveal the different spatial resolutions between the TDWRs (beam width of 0.5$^\circ$) and NEXRAD  (beam width of 1$^\circ$). Through the remainder of the day, the mobile radar (RaXpol, \cite{Pazmany+2009}) was stationed at the roost site scanning continuously in RHI mode for nine hours, generating RHI images every thirteen seconds and resulting in nearly 2500 RHI images. Figure \ref{NEXRAD:rax}a shows one such image, with the region outlined by the black box shown in Fig.~\ref{NEXRAD:rax}b. From manual inspection, the high intensity signals below 0.5 km in altitude are ground and electronic clutter. Similarly, the low intensity arc at a range of 0.6 km is 63 due to the radar receiver. Confirmed by visual observations at the time of the measurements, the remaining moderate intensity blobs consist entirely (or at least predominantly) of purple martins. The pixels above the white dashed line in Fig. \ref{NEXRAD:rax}a are considered noise, and used to calculate the noise statistics \citep{Stepanian++2014}. A threshold of 0.3 standard deviations above the mean noise level was chosen, resulting in the binary image in Fig. \ref{NEXRAD:rax}c. We determined that birds could be adequately identified by their size alone, resulting in the filtering of blobs above 10 pixels and below 2 pixels in size. Figure \ref{NEXRAD:rax}d shows the resulting filtered image containing only purple martin-like signals. Having confirmed that the identifications were correct on a selection of images, the algorithm could then be applied to the full nine hours of RHI images.  In each image, the number of birds as a function of height were counted, resulting in a time-altitude-count plot (Fig.~\ref{NEXRAD:Garland})

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_rax.png}
\caption{a) Backscattered power RHI image (no range correction) from the RaXpol with white dashed line separating noise-only heights. b) Blowup of the boxed region in (a), showing the region of the aerosphere containing birds. c) Binary image resulting from the thresholding of (b). d) Filtered power image with only bird signals.}
\label{NEXRAD:rax}
  \end{figure}


\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_Garland.pdf}
\caption{The time evolution of Purple Martin height distributions at the Garland, Texas roost on 01 August 2011.}
\label{NEXRAD:Garland}
  \end{figure}

The patterns in Fig.~\ref{NEXRAD:Garland} indicate that Purple Martins initially exit the roost and remain concentrated in the lowest kilometer of the atmosphere.  As the sun drives boundary layer development, the birds use thermals to gain altitude up to 3 km.  This behavior is likely driven by the passive lofting of insects by thermals as the boundary layer deepens.  In response to the accumulation of this food source aloft, the birds are willing to expend the extra energy to climb to these altitudes to feed.  

Unfortunately, KFWX was not polarimetric at the time of this study.  To explore the polarimetric signatures of a Purple Martin exodus, we revisit KEWX on 03 August 2014 (i.e., the same radar and date as the diurnal and bat examples).  Figure~\ref{NEXRAD:EWXpuma} shows the radar products at the time of the morning Purple Martin emergence.  Two roosts can be identified by their divergence signatures in radial velocity: one to the north-east and one to the south-west.  An additional area of enhanced bioscatter is west of the radar, and is likely the morning reentry of bats at Bracken Cave.  The remaining background bioscatter is made up of insects, as evidenced by the low velocity and high $Z_{DR}$.  


\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_EWX_puma.pdf}
  \caption{NEXRAD products from Austin, Texas (KEWX) on 03 August 2014. Products are (clockwise from upper left) radar reflectivity factor ($Z$, [dBZ]), radial velocity ($v_r$, [m s$^{-1}$]), spectrum width ($\sigma_v$, [m s$^{-1}$]), differential phase ($\phi_{DP}$, [deg]), co-polar correlation coefficient ($\rho_{HV}$, [~~]), and differential reflectivity ($Z_{DR}$, [dB]).}
  \label{NEXRAD:EWXpuma}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_EWX_pumaZoom.pdf}
  \caption{NEXRAD products from Austin, Texas (KEWX) on 03 August 2014, zooming in on the north-east roost emergence.  Products are (clockwise from upper left) radar reflectivity factor ($Z$, [dBZ]), radial velocity ($v_r$, [m s$^{-1}$]), spectrum width ($\sigma_v$, [m s$^{-1}$]), differential phase ($\phi_{DP}$, [deg]), co-polar correlation coefficient ($\rho_{HV}$, [~~]), and differential reflectivity ($Z_{DR}$, [dB]).}
  \label{NEXRAD:EWXpumaZoom}
\end{figure}

Focusing in on the roost to the north-east, it is apparent that the signature is made up of several concentric rings of enhanced reflectivity (Fig.~\ref{NEXRAD:EWXpumaZoom}).  Having witnessed the morning exodus of Purple Martins, it is often the case that the colony departs in several discrete stages or pulses, resulting in several expanding rings.  Unlike the bat emergences, which did not surpass $-1$ dB, the Martin signature has large regions of negative differential reflectivity, with areas below $-4$ dB.  Again, due to the effect of cross-polar contributions, a physical interpretation of $Z_{DR}$ is difficult, but because both the bat and martin cases are observed by the same radar on the same day, we can assume that the differential transmit phase is constant between the two. 

Similar to the spectrum width in bat emergences, low values at cross-beam orientations indicate flight in strict divergence from the cave with a variety of airspeeds.  Correlation coefficient, is low in regions of positive differential reflectivity, and high in regions of negative differential reflectivity.  Furthermore, the regions of lowest correlation coefficient are coincident with areas of high spatial variability in differential phase. \cite{Ryzhkov2007} demonstrates that the effects of nonuniform beam filling can be diagnosed by lowering of correlation coefficient and azimuthal variability in differential phase.  With this in mind, it is likely that the aerial distribution of Purple Martins in regions of especially low correlation coefficient are non beam filling, resulting in these signatures.

\section{Widespread Nocturnal Migration}
\label{NEXRAD:mig}
\begin{figure}[h!]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_seasons.pdf}
  \caption{Seasonal variations in the polarimetric fields of the $0.5^\circ$ PPI from KTLX (Oklahoma City, Oklahoma). Tiles cover a 300 km by 300 km square domain, centered on the radar site, and are the first volume scan following midnight local time (i.e., 05 UTC in Spring, Summer, and Autumn, and 06 UTC in Winter). Cases showing widespread migratory alignment (i.e., Spring and Autumn) are annotated with an approximate orientation vector.}
\label{NEXRAD:seasons}
\end{figure}
\sloppy
The phenomenon of avian migration is marked by the widespread, directional, purposeful movements of millions of individuals.  This shared behavior results in aligned individuals across large regional expanses.  As a seasonal phenomenon, not all widespread bioscatter is migration; however, polarimetry can indicate the presence of migration by azimuthal signatures resulting from mass alignment.  For many birds, migration is a nocturnal activity, initiating at sunset and continuing until sunrise.  A typical case of nocturnal bird migration is shown in the Spring case at Oklahoma City, Oklahoma (KTLX, Fig.~\ref{NEXRAD:seasons}) in which the radar sweeps are the $0.5^\circ$ elevation angle immediately following midnight local time.  This spring migration is characterized by high reflectivity, exceeding 25 dBZ, indicating a high number of birds in the airspace.  Radial velocity magnitudes reach 30 m s$^{-1}$, indicating either strong winds or active, directional flight.  Differential reflectivity ranges from $-4$ through $4$ dBZ with strong azimuthal variability, forming apparent quadrants.  Both differential phase and correlation coefficient display strong symmetry across the flight orientation vector (solid arrow), and correlation coefficient shows consistent antisymmetry across the perpendicular axis (i.e., head-tail axis, dashed line).  


As spring transitions to summer, migrating birds leave Oklahoma for their northern breeding grounds.  The summer resident birds of Oklahoma have no need to make high altitude flights at night, and as a result, there is an absence of birds in the nocturnal airspace during the summer.  Nonetheless, nighttime bioscatter is still quite high in the summer, reaching 35 dBZ (Fig.~\ref{NEXRAD:seasons}).  These signals are made up primarily of insects, most of which have been pumped aloft by thermals in the convective boundary layer, and remain suspended in the nocturnal residual layer.  Additionally, the airspace in some cases may contain bats that are foraging on these insects.  These processes result in low velocities, likely representative of the bulk winds.  Differential reflectivity is high, often exceeding 7 dBZ, again indicating insects.  Differential phase and correlation coefficient have less azimuthal structure and are comparatively unorganized.  

The onset of fall migration yields nearly identical patterns as the spring, only oriented southward (Fig.~\ref{NEXRAD:seasons}).  Birds from Canada and the Northern United States again fill the airspace in nocturnal migration over central Oklahoma, resulting in high reflectivity and velocity magnitudes.  Again, differential phase and correlation coefficient are both highly symmetric along the flight orientation of the organisms.  Differential reflectivity is similar as the spring case, but has developed a secondary lobe of negative values that is symmetric with the original lobe.  Through the winter, bioscatter is limited.  The migratory birds have moved south to their wintering grounds, and resident birds have no reason to fly at altitude, especially at night.  The result is negligible reflectivity and velocity, and unorganized polarimetric variables.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_syn.png}
  \caption{Synoptic map for United States on 03 May 2013 at 05:00 UTC (local midnight). Prepared by the National Centers for Environmental Prediction, Weather Prediction Center.}
\label{NEXRAD:syn}
\end{figure}

These migratory patterns over Oklahoma City are not unique; Similar polarimetric morphology can be observed at many (if not all) radar sites when bird migration is ongoing.  The following sample case consists of a spring migration event in the northeastern United States in which migrants (primarily passerines) were utilizing the southerly winds preceding the passage of a synoptic front (Figs.~\ref{NEXRAD:syn} and \ref{NEXRAD:NE}). Such migration episodes occur on an approximately weekly basis during the spring and autumn months, and form similar morphological features.  The location of seven radar sites is annotated in this non-quality controlled, composite reflectivity (CREF) image.  The square domains cover a 300 km by 300 km tile centered on each radar, and serve as the analysis regions for the following figures.  

Figure~\ref{NEXRAD:Z} displays the underlying radial (i.e., Plan Position Indicator, PPI) data for the domains annotated in Fig.~\ref{NEXRAD:NE} that were used to form the CREF.  The radar operating frequency for each site is listed with the site location and station identifier.  All PPIs are at the 0.5$^\circ$ elevation angle.  Columns span the thirty minutes preceding and following 05 UTC, at ten-minute resolution.  The corresponding PPI images of Doppler radial velocity ($v_r$, Fig.~\ref{NEXRAD:V}), differential phase ($\psi_{DP}$,  Fig.~\ref{NEXRAD:P}), co-polar correlation coefficient ($\rho_{HV}$,  Fig.~\ref{NEXRAD:C}), and differential reflectivity ($Z_{DR}$,  Fig.~\ref{NEXRAD:ZDR}) are included for comparison.

\begin{figure}[t]
\centering
\includegraphics[width=3.5in]
  {figures/NEXRAD_NE.pdf}
  \caption{Composite reflectivity for a region of the northeastern United States on 03 May 2013 at 05:00 UTC (local midnight).  Circles indicate the location of the seven sample radar sites.  The square domains surrounding each site are 300 km by 300 km, and correspond to the tiles presented in all following figures. Image constructed at http://soar.ou.edu/legacy.html}
\label{NEXRAD:NE}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_Z.pdf}
  \caption{NEXRAD 0.5$^\circ$ PPIs of radar reflectivity factor ($Z$) for the seven radar sites for the 1-hour period surrounding 05 UTC.}
\label{NEXRAD:Z}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_V.pdf}
  \caption{NEXRAD 0.5$^\circ$ PPIs of Doppler radial velocity  ($v_r$)for the seven radar sites for the 1-hour period surrounding 05 UTC.}
\label{NEXRAD:V}
\end{figure}

Radial velocity indicates the widespread northward motions typical of spring migrants (Fig.~\ref{NEXRAD:V}).  The morphological features of $\psi_{DP}$ images in migration are quite consistent across time, with strong azimuthal symmetry across the orientation axis (Fig.~\ref{NEXRAD:P}).  Slight variations exist at the off-tail view angles, which can be highlighted by wrapping in the colormap; for example, consider the dark blue colors that symmetrically correspond with reds in KBGM.  As with most of the radars, these regions have slight discrepancies in $\psi_{DP}$ magnitudes, although the color scale accentuates it most in the KBGM figure.  Across sites, the contributions of the system differential phase result in different magnitudes with site-specific shifts.  Again, because biological scatterers can have potentially large and highly variable $\delta$ values, combined with the continuous cyclic nature of phase, there is an ambiguity of the direction of phase shifts near 180$^\circ$. From a visual standpoint, $\psi_{DP}$ serves as an effective means of deducing orientation, with major morphological features diagnosing symmetry axes.

Correlation coefficient is similar across sites, with symmetry across the flight orientation axis, and antisymmetry across the perpendicular head-tail axis (Fig.~\ref{NEXRAD:C}).  In all cases, the half of the domain viewing tailward angles are lower, with values around 0.3, while headward view angles are high and can approach 1.   

In the case of NEXRAD, differential reflectivity is the least morphologically consistent radar product.  Some studies utilizing other radar systems have presented highly symmetric $Z_{DR}$ fields that could have the potential of yielding orientation information (e.g., \cite{Lang++2004}).  But even in strong migration with widespread alignment, NEXRAD $Z_{DR}$ fields are generally asymmetric (Fig.~\ref{NEXRAD:ZDR}).  As discussed in the previous sections, this asymmetry in biological echoes is a result of the Simultaneous Transmit And Receive (STAR) configuration employed by NEXRAD \citep{Melnikov++2015}.  As described in \cite{Ryzhkov+2007}, one consequence of STAR configuration is the potential for cross-coupling of the orthogonally polarized wave components.  This phenomenon is negligible for most hydrometers (with some exceptions presented in \cite{Ryzhkov+2007}); however, the irregular shape of biological scatterers results in non-negligible depolarization (i.e., cross-polar) contributions \citep{Melnikov++2015}.  The resulting signal received at the radar antenna is the coherent sum of the original polarization component and the depolarized wave, with the phase offset between the two contributions determining the resultant wave.  A dominant factor influencing this phase offset at the receiver is the system differential phase upon transmission---a quantity that is not recorded or controlled for NEXRAD systems \citep{Ryzhkov+2007,Melnikov++2015}. Because of this variation in system differential transmit phase, the morphology of azimuthal asymmetries can vary across NEXRAD sites. Furthermore, it has been shown that variations in radar frequency within S-band, such as those employed among NEXRAD sites, can result in inherent variation in $Z_{DR}$ measurements \citep{Melnikov++2012}.  From inspecting many $Z_{DR}$ images, we note that the morphology of these asymmetries seem to remain stable through a migratory season, but they vary substantially among radar sites, further supporting a cause linked to variations in radar system transmit phase or frequency.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_P.pdf}
  \caption{NEXRAD 0.5$^\circ$ PPIs of differential phase ($\psi_{DP}$) for the seven radar sites for the 1-hour period surrounding 05 UTC.}
\label{NEXRAD:P}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_C.pdf}
  \caption{NEXRAD 0.5$^\circ$ PPIs of co-polar correlation coefficient ($\rho_{HV}$) for the seven radar sites for the 1-hour period surrounding 05 UTC.}
\label{NEXRAD:C}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_ZDR.pdf}
  \caption{NEXRAD 0.5$^\circ$ PPIs of differential reflectivity ($Z_{DR}$) for the seven radar sites for the 1-hour period surrounding 05 UTC.}
\label{NEXRAD:ZDR}
\end{figure}

\clearpage
\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_ZR98.pdf}
  \caption{The relationship between differential scatter phase and differential reflectivity for birds and insects, reproduced from \cite{Zrnic+1998}.  Boxes have been annotated to distinguish approximate regions for birds (red) and insects (blue).  These same annotations are included in Figs. \ref{NEXRAD:scatter1} and \ref{NEXRAD:scatter2}.}
\label{NEXRAD:ZR98}
\end{figure}

The relationship between the differential scatter phase and differential reflectivity has been shown to be diagnostic of the type of biological scatterers filling the resolution volumes \citep[][see Fig.~\ref{NEXRAD:ZR98}]{Zrnic+1998}.  While biologically-driven reasoning alone can support our assumption of bird-dominated scattering in this case, it is worth checking whether these clustering boundaries are maintained.  Figures \ref{NEXRAD:scatter1} and \ref{NEXRAD:scatter2} show the scatter for $\psi_{DP}-Z_{DR}$ relationships for the seven radars from the previous case, using the PPIs from the first scan following 05 UTC (i.e., local midnight).  Values of $\psi_{DP}$ can be thought of as $\delta$ with a constant offset, shifting the scatter up or down the vertical axis.  For these cases of migration, bird scatter spans much of the dynamic range in both differential reflectivity and phase.  Furthermore, of the hundreds of cases analyzed that displayed similar morphology, all produced similar $\psi_{DP}-Z_{DR}$ scatter plots.


\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_scatter1.pdf}
  \caption{Scatter plots of differential phase and differential reflectivity for four of the seven radars denoted in Fig. \ref{NEXRAD:NE} for the first scan following 05:00 UTC.  Boxes distinguishing birds (red) from insects (blue) from Fig. \ref{NEXRAD:ZR98} have been included for comparison.}
\label{NEXRAD:scatter1}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=6in]
  {figures/NEXRAD_scatter2.pdf}
  \caption{Scatter plots of differential phase and differential reflectivity for three of the seven radars denoted in Fig. \ref{NEXRAD:NE} for the first scan following 05:00 UTC.  Boxes distinguishing birds (red) from insects (blue) from Fig. \ref{NEXRAD:ZR98} have been included for comparison.}
\label{NEXRAD:scatter2}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Extracting Flight Orientation Profiles from Polarimetric Radar}
\label{align}

\section{Background}
\label{align:background}
Much like the VAD, polarimetric quantities can be presented as a function of azimuth as a way to visualize angular variations in scattering characteristics \citep{Bachmann+2007}.  While the Spectral VAD (SVAD) methods demonstrated in \citep{Bachmann+2007} provide insight into the fundamental scattering physics of birds and insects, many practical SVAD applications are limited by the general lack of spectral data from operational radars.  As a result, there is a continuing lack of biological, polarimetric methods that can be implemented on readily-available, i.e., level II, radar data.  Several studies have illustrated the presence of symmetric and quasi-symmetric morphological features in radar images of uniform mass migrations, which under the assumption of horizontal homogeneity, can be diagnostic of flight orientation \citep{Lang++2004,Diehl+2005, Westbrook2008, vanDenBroeke2013}.  

Provided that the axis of azimuthal symmetry in $\rho_{HV}$ is coincident with the axis of mean animal orientation, radar polarimetry can provide a basis for diagnosing migratory flight direction. The consistency of $\rho_{HV}$ morphology over time and across NEXRAD sites makes it a good candidate for automated analysis. Furthermore, the symmetry and oscillatory nature of $\rho_{HV}$ provide a convenient parallel to the VAD technique for pure translation \citep{Browning+1968}. As illustrated in the schematics in Fig.~\ref{align:fig2}, the flight orientation corresponds with the global $\rho_{HV}$ minimum in the top two cases and the local $\rho_{HV}$ maximum in the bottom case.  While the azimuthal variation of $\rho_{HV}$ is different in each of the three morphological cases, the flight orientation in all three can be obtained by a common method. As illustrated in Fig.~\ref{align:fig2}, a sinusoidal fit to $\rho_{HV}$ that cycles over 360 degrees provides a global minimum that is collocated with the axis of symmetry and flight orientation angle. 

\begin{figure}[t]
\centering
\includegraphics[width=4in]
  {figures/align_fig2.pdf}
  \caption{Conceptual schematics of the three common $\rho_{HV}$ morphologies. (left) Idealized PPIs of symmetric $\rho_{HV}$ patterns. (right) The corresponding azimuthal functions of $\rho_{HV}$ and their idealized sinusoidal fits.}
\label{align:fig2}
\end{figure}

While quite simplistic, this technique yields results that are consistent with those obtained by human inspection of the radar images, provided several conditions are met:
\begin{enumerate}
\item\textbf{The radar domain is well-filled with organisms.}  Because this technique relies on fitting a model to $\rho_{HV}$ measurements across all azimuths, the measurements must provide enough data points to enable this model fit. This condition breaks down in weak or localized animal movements, or when movements are spatially constrained or poorly sampled (e.g., lack of organisms over ocean at coastal radar sites). The effect of this condition is investigated in Section~\ref{align:sensitivity}.
\item\textbf{Flight orientation is horizontally homogeneous across the domain.} Much like the VAD technique, any variation in the mean flight orientation is assumed to exist only in height.  The validity of this assumption is typically maintained for large scale nocturnal migration, provided the absence of significant meteorological inhomogeneities (e.g., synoptic fronts or precipitation) or topographic obstacles.
\item\textbf{Organismal composition is horizontally homogeneous across the domain.} As above, the condition of horizontal homogeneity ensures azimuthal symmetry in the $\rho_{HV}$ field, and enables the model fit. 
\end{enumerate}
In many cases, the validity of these conditions cannot be known \emph{a priori}; however, several indicators of retrieval quality are presented in Sections \ref{align:method} and \ref{align:demo}.

\section{Computational Method}
\label{align:method}
Let $\rho_{HV}(\phi_i,r_j)$ denote the volume measurement of $\rho_{HV}$ at the azimuth, $\phi_i$, and range gate, $r_j$. Starting with a $0.5^\circ$ PPI of $\rho_{HV}$, flight orientation extractions are conducted independently at each range gate.  Following \cite{Browning+1968}, the azimuthal oscillation of a sinusoid with a wavelength of 360 degrees can be decomposed into the Fourier series of the form
\begin{equation}
\widehat{\rho_{HV}}(\phi,r_j) = c_{0,j} + c_{1,j} \cos(\frac{\pi}{180^\circ}\phi) + c_{2,j} \sin(\frac{\pi}{180^\circ}\phi)
\end{equation}
in which $\widehat{\rho_{HV}}$ denotes the sinusoidal fit to the $\rho_{HV}$ measurements, azimuth angle ($\phi$) is in degrees clockwise from north, and $c_{0,j}$,  $c_{1,j}$, and $c_{2,j}$ are the Fourier coefficients defining the sinusoidal fit at the $j^{th}$ range gate.  Using the $\rho_{HV}$ measurements, the Fourier coefficients for the $j^{th}$ range gate can be obtained by optimizing the system of equations
\begin{multline}
\begin{bmatrix}
1 & \cos(\frac{\pi\phi_1}{180^\circ}) & \sin(\frac{\pi\phi_1}{180^\circ}) \\
1 & \cos(\frac{\pi\phi_2}{180^\circ}) & \sin(\frac{\pi\phi_2}{180^\circ}) \\
\vdots & \vdots & \vdots \\
1 & \cos(\frac{\pi\phi_n}{180^\circ}) & \sin(\frac{\pi\phi_n}{180^\circ}) \\
 \end{bmatrix}
 \begin{bmatrix}
c_{0,j} \\
c_{1,j} \\
c_{2,j}
 \end{bmatrix}
 =
\begin{bmatrix}
\rho_{HV}({\phi_1,r_j}) \\
\rho_{HV}({\phi_2,r_j}) \\
 \vdots \\
\rho_{HV}({\phi_n,r_j})
 \end{bmatrix}
\end{multline}
using least-squares or a similar method.  Once the Fourier coefficients have been obtained, they can be inserted in (6.1) to find the functional form of the sinusoidal fit, $\widehat{\rho_{HV}}(\phi,r_j)$.  The flight orientation at the $j^{th}$ range gate is the azimuth corresponding to the minimum in $\widehat{\rho_{HV}}(\phi,r_j)$.  Assuming standard refraction, these range-dependent flight orientations can be mapped to height above ground level ($h$) using 
\begin{equation}
h_j = \sqrt{r_j^2+(\frac{4}{3} r_0)^2 + \frac{8}{3}r_0r_j\textrm{sin}(\theta)} -\frac{4}{3} r_0\\
\end{equation}
with the elevation angle, $\theta$, and the Earth radius, $r_0$ \citep{DoviakZrnic}. 

\begin{figure*}[b]
\centering
\includegraphics[width=6in]
  {figures/align_fig3.pdf}
  \caption{Automated orientation profile extraction method: (left) Sample $\rho_{HV}$ PPI with 25 km and 75 km range rings. (center) Scatter plots of $\rho_{HV}$ verses azimuth and sine fits for the 25 km and 75 km ranges. (right) Extracted flight orientation profile.}
\label{align:fig3}
\end{figure*}

Figure~\ref{align:fig3} demonstrates this process for the 25-km and 75-km range gates.  Beginning with a PPI of $\rho_{HV}$ (Fig.~\ref{align:fig3}, left), the measured values at each range gate are presented as a function of azimuth (Fig.~\ref{align:fig3}, center).  Using (6.1) and (6.2) on these $\rho_{HV}$-Azimuth Displays produces the sine fits that are overlaid, and the minima indicate the flight orientation.  The ranges corresponding with these minima are substituted in (6.3) to yield flight orientation as a function of height (Fig.~\ref{align:fig3}, right).

The three conditions listed in Section~\ref{align:background} provide the ideal case requirements for applying this orientation extraction technique; however, the model fit will still return results if conditions are not met.   Several techniques can be applied before, during, or after the orientation retrievals to reveal the validity of these conditions. 

Condition 1 requires the domain to be filled with organisms such that all of the azimuths at a given range contain data.  Missing data points result in less information to guide the model fit, contributing to retrieval errors.  The amount of missing data can easily be obtained prior to the analysis, and may be dealt with in several ways.  One option is to only perform retrievals on ranges within a minimum threshold of missing data.  For example, any ranges missing more than 25\% of azimuth data points can be flagged as `retrieval not possible'.  An alternative that is sometimes used in VAD applications is to limit analysis to some range annulus which should be well-filled with organisms.  In this case, all ranges corresponding with altitudes under 100 meters or over 1.5 kilometers, for example, could be ignored to avoid clutter or organism-free regions.  Finally, it is possible to perform retrievals on all ranges and construct orientation error bounds with respect to the amount of missing data points.  This process is the topic of Section~\ref{align:sensitivity}.
 
Condition 2 requires that organisms are in constant alignment horizontally.  In the extreme opposite case of random orientation, the azimuthal variation in $\rho_{HV}$ will be random and the Fourier fit will return a sinusoid with minimal amplitude, i.e., a nearly flat line through the measurements.  In this case, the orientation angle associated with this meaningless minimum will be returned as the flight orientation.  This situation can be avoided by setting a minimum magnitude of sinusoid amplitude that must be exceeded for the retrievals to be valid.  Another possibility is to calculate a measure of goodness of fit between the $\rho_{HV}$ measurements and the sine fit.  Because the sine model was chosen for convenience rather than upon a physical basis, most goodness of fit measures will yield poor values even in ideal cases.  Nonetheless, these metrics can still potentially identify non-migration cases by their exceptionally poor fits.
 
Condition 3 requires that organismal composition is horizontally homogenous, resulting in symmetric $\rho_{HV}$ fields.  As above, goodness of fit metrics can provide a rough assessment of this condition.  An alternative technique can provide an overall quality check once all retrievals have been performed, and can help to evaluate whether the three conditions have been met.  In general, the orientation of migrants should not vary substantially between consecutive altitude bins or time steps; or in other words, orientation retrievals should display spatiotemporal continuity.  Measures of 2D variability, such as standard deviation within a 3$\times$3 square structuring element, can reveal anomalous results indicative of poor retrievals. A demonstration of these techniques is presented in Section~\ref{align:demo}


\section{Comparisons to VAD Techniques}
\label{align:compare}
The current standard for obtaining flight velocity profiles from radar is the Velocity Azimuth Display (VAD, \cite{Browning+1968,Gauthreaux+1998}).  The result of the VAD is the \emph{ground speed} of migrants, or their velocity with respect to the stationary Earth surface. As such, these ground velocities are the vector sum of the active flight velocity, or \emph{air speed}, of the organisms and the ambient atmospheric motions. The orientation provided by polarimetric analysis is associated with the intentional orientation of organisms in flight, corresponding with their air speed velocity. Thus, obtaining the corresponding orientation vector requires subtracting an ambient wind measurement from the VAD result \citep{Larkin1980}. Several sources of routine wind data are available for such calculations. These data sources can either be wind measurements from surface anemometers or balloon-borne soundings \citep{Buler+2014}, or the output wind fields from atmospheric models \citep{Jiang++2013}.  

Balloon soundings are launched twice daily from National Weather Service forecast offices at 00 UTC and 12 UTC. Vertical profiles of wind speed and direction are recorded incrementally as a function of pressure level, and can be downloaded from the National Climatic Data Center (NCDC, \cite{NCDC}) or the University of Wyoming's upper air data portal \citep{UW}.  In general, these balloon launch sites are located near NEXRAD installations, providing local wind information that can be used for correcting VADs. While this close proximity is the ideal case, many radar locations do not have collocated soundings.  Of the 143 NEXRAD sites in the contiguous US, only $41\%$ have soundings within a 50 km range, and $53\%$ within a 100 km range.  In these cases, the nearest neighboring sounding is typically used.

As an alternative, the wind fields from an atmospheric model may be used to correct VADs  \citep{Jiang++2013}.  The North American Regional Reanalysis (NARR, \cite{NARR}) is a combined model and assimilated dataset containing wind vectors at 29 geopotential height levels. NARR model output can be downloaded from the NCDC, with a model resolution of 3-hours and 32-km grid spacing \cite{NCDC}.  Another model option available at the NCDC is the Rapid Refresh (RAP), which provides 1-hour updates at 50 vertical levels and 13-km horizontal grid spacing \citep{RAP}. While the concept of airspeed retrieval is simple---the vector subtraction of wind contributions---the process is quite sensitive to the accuracy of wind data.

The following comparison focuses on the Albany, New York NEXRAD site (KENX, 2865 MHz) during the Spring 2013 migratory season. Level II radar data spanning February through June 2013 were downloaded from the NCDC, and the nightly analysis period was limited to include all PPIs taken between civil twilights (approximately 01 UTC through 09 UTC). The PPIs containing weather signals were removed by manual inspection. The resulting dataset contained a total of 982 volume scans.  Velocity dealiasing was performed as needed using \citep{Sheldon++2013}, and a VAD was calculated for each $0.5^\circ$ PPI following \citep{Browning+1968}.


\begin{figure}[t]
\centering
\includegraphics[width=4in]
  {figures/align_map.pdf}
  \caption{Analysis domain for the Albany, New York (KENX) orientation retrieval comparisons. Radiosonde launch sites were located at the Albany (ALB) and Upton (OKX) weather forecast offices.  Two range rings are included for reference.}
\label{align:map}
\end{figure}


To evaluate the use of soundings for VAD correction, balloon-borne wind profiles were downloaded for two radiosonde launch sites.  The first site is located at the Albany weather forecast office [Fig.~\ref{align:map}, ALB], approximately 22 km from the radar location.  To represent radar sites without collocated soundings, the nearest neighboring sounding---taken 218 km away in Upton, NY [Fig.~\ref{align:map}, OKX]---was used. To test model-based VAD corrections, model fields from the NARR and RAP were downloaded from the NCDC. For both models, the wind profiles for the grid cells encompassing KENX were extracted.  Finally, all sounding and model wind profiles were cast onto a common 10-meter vertical coordinate system using nearest-neighbor interpolation. To test the temporal variation in sounding-based corrections, VADs were corrected using the 00 UTC, 12 UTC, and temporally-nearest soundings.  Both model-based corrections used the temporally-nearest model output.  In all, a total of eight corrections were tested.


Following the method described in Section~\ref{align:method}, an orientation profile was extracted from each of the 982 PPIs. Additionally, the percent of missing azimuth data and coefficient of determination, $R^2$, were calculated for each retrieval.  To eliminate data biased by ground clutter, sparse animal distributions, or non-migration, range gates with $R^2<0.3$ were removed.  The resulting orientation profiles were cast onto the same 10-meter vertical coordinate system as VADs and winds using the mean orientation value in each height bin.  These orientation profiles were subtracted pointwise from the ground speed vector profiles from each corresponding VAD, as well as the eight air speed vector profiles from each wind-corrected VAD.  The resulting absolute offsets between polarimetric orientation and velocity-based techniques were averaged within bins of one hour and fifty meters in height (Fig.~\ref{align:comp}).


\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/align_comp.pdf}
  \caption{Comparisons of nine velocity-based orientation retrievals with polarimetric retrievals. Subscripts list the wind correction data source, with `nearest' denoting the temporally-nearest measurement. Each pixel value is the mean absolute offset of the 982 cases, binned by 1 hour and 50 meters.} 
\label{align:comp}
\end{figure}

Figure~\ref{align:comp}a shows the offset between the airspeed orientation retrieved from the $\rho_{HV}$ field and the ground speed orientation directly from the VAD.  These offsets are the result of three possible contributions. First, while migrants typically attempt to use winds to aid their movements, they often are not perfectly aligned with the wind.  In this case, offsets represent a physical deviation between airspeed and ground speed, indicating the degree of migrant wind drift or compensation against the wind. The second source of offsets results from a poor model fit to the $\rho_{HV}$ data.  This happens most commonly at high altitudes where azimuthal data loss contributes to large retrieval errors.  Additionally, in some cases of two-minima morphology (as in Fig.~\ref{align:fig2}, bottom), the magnitude of  the low $\rho_{HV}$ lobes is slightly different, resulting in a biased retrieval of several degrees toward the deeper lobe. The final source of offsets results from the poor performance of the VAD.  These cases are especially problematic in high-wind cases when velocity aliasing (and subsequent dealiasing) can yield errors in the velocity fields, biasing VAD retrievals.  Overall, the orientation retrievals from $\rho_{HV}$ and VADs are not expected to match, and so offsets in Fig.~\ref{align:comp}a can either be physical (i.e., due to migrant wind drift or compensation) or retrieval errors in either method. Conversely, the following wind-corrected VADs should correspond to those orientations retrieved from $\rho_{HV}$, and all offsets must be the combined result of retrieval errors in both methods.


The offsets in Fig.~\ref{align:comp}b-d show the deviations between polarimetric orientations and VADs corrected by the nearby ALB soundings. The VAD in (b) was corrected using the 00 UTC sounding, taken at 7:00pm local time.  Through much of the spring in Albany, this sounding is taken before local sunset, with winds characterizing those of the convective boundary layer.  In many cases, these wind profiles are almost immediately invalid following the evening transition into the nocturnal residual layer.  The errors associated with subtracting an invalid wind profile can often exceed the magnitude of uncorrected VAD offsets.  This effect can be seen in Fig.~\ref{align:comp}b as offsets increase steadily through the night, as the boundary layer stabilizes and diverges from the convective sounding at 00 UTC. The reverse effect is observed when using the 12 UTC ALB sounding (Fig.~\ref{align:comp}d), with offsets steadily decreasing through the night as the nocturnal boundary layer wind fields approach the 12 UTC measurement.  Figure~\ref{align:comp}c uses the temporally-nearest ALB sounding to correct the VADs (i.e., 00 UTC for VADs taken before 06 UTC, and 12 UTC for VADs taken after).  In many cases this technique may help mitigate errors due to the temporal variation in wind fields when applying sounding-based corrections; however, caution should be used to avoid incorporating springtime 00 UTC convective soundings to represent nocturnal winds.  This caveat is especially significant in regions of prominent nocturnal wind variations such as coastal land breezes, katabatic winds, or nocturnal low-level jets.  Figure~\ref{align:comp}e-g shows the corresponding offsets resulting from wind corrections from the UPT soundings.  At 218 kilometers away, the Upton wind profiles are affected by spatial variability in horizontal winds, resulting in overall higher offsets as compared to the ALB correction.

The offsets associated with model-based corrections from the NARR (Fig.~\ref{align:comp}h) and RAP (Fig.~\ref{align:comp}i) have significantly different results.  On average the NARR provides minor improvements over the original VADs.  Conversely, the RAP increases offsets overall---despite its higher temporal and spatial resolution.  Two possible causes may be the sources of these discrepancies.  First, the winds used for VAD correction are within the lowest 3 km above the ground level, well within the effects of surface features such as small scale topography, land cover, and other roughness effects that are not modeled. These effects can lead to slight deviations in both models. The second, and most significant, source of discrepancy is the data that is assimilated into the models.  The NARR wind assimilation includes several measurements including those from rawinsondes, pibals, and aircraft \citep{NARR}.  The RAP assimilation also includes wind measures from NEXRAD VADs and wind profilers---both of which can be biased by migrating animals, particularly during time periods such as those presented in this analysis \citep{RAP}.  As a result, it is expected that the NARR would better reflect true air motions, while the RAP could include contributions from the very animals that VAD corrections are meant to characterize.
 
While this example is not meant to be a comprehensive analysis of VAD corrections by soundings and models, it does illustrate the need for such future work.  There is a clear disparity between polarimetric orientation retrievals and those obtained by velocity corrections.  Step-by-step inspection of the underlying retrieval processes has revealed that some contribution to these offsets is a result of inaccurate VAD wind corrections, which we attribute primarily to spatiotemporal inaccuracies in underlying wind data. The extent to which each method contributes to the final offsets cannot currently be quantified within this study.  It is possible that visual ground truth via vertically-pointing thermal video recording may be the only true source of accuracy assessments for these two radar methods, and should be the topic of a followup study.


\section{Sensitivity Analysis}
\label{align:sensitivity}
The method described in Section~\ref{align:method} relies on fitting a function to a discrete set of data points. As a result, the availability of these data points will influence the model fit and the resulting orientation retrieval. In the ideal case, the full radar domain would be filled with organisms such that all azimuths contain a data point. In realistic cases, however, this may not be true. Losses in data availability can be split into two types. The first type is the absence of sparse data points. This scattered lack of data can be the result of a localized absence of organisms or censoring of small clutter sources. The second type is the loss of contiguous regions of data. These losses can be the result of a large-scale absence of organisms, radar beam occultation, or anomalous propagation. The following investigates the method precision subject to data losses of these two types.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/align_sensitivity2.pdf}
  \caption{(a) The frequency of retrieval error with respect to the 0\% missing data retrieval for increasing azimuthal data losses. Each column is associated with a 1-deg reduction in data, with cases summing to 1,345,800.  The 95\% confidence interval of retrieval errors is included for reference. (b) The mean retrieval error with respect to no data loss as a function of missing data wedge width and position. (c) The standard deviation of retrieval error with respect to no data loss as a function of missing data wedge width and position.}
\label{align:sens}
\end{figure}

\subsection{Sensitivity to Sparse Data Loss}
To quantify the precision of the method as random data is lost, the same Spring 2013 KENX dataset (as in Section~\ref{align:compare}) was analyzed.  For each of the 982 PPIs, only range gates with data in at least 99.0\% of azimuths were considered, resulting in 13,458 distinct range cases.  The orientations associated with these well-filled ranges served as baseline cases to which degraded retrievals were compared.  For each of the 13,458 cases, random azimuths were sequentially removed one at a time until only 2\% of azimuths remained, and the resulting orientation retrievals from each reduction were compared to the original.  To avoid anomalous results due to `lucky' or `unlucky' random seeds, the process of removing random azimuths was repeated 100 times for each case, yielding a total of 1,345,800 results for each reduction in data.  The resulting errors from these 1,345,800 cases as a function of data reduction are presented in Fig.~\ref{align:sens}a.  The frequencies in each column of Fig.~\ref{align:sens}a sum to 1,345,800 (i.e., the total number of cases), and can be thought of as discrete probability density functions of expected error for a given percent of missing data.  As a result, confidence intervals can be computed to indicate the error associated with an indicated percent of cases. For example, when a random 25\% of azimuth data is removed, 95\% of cases had errors between $\pm3.1$ degrees. When azimuth data was reduced to 50\%, 95\% of cases had errors between $\pm5.4$ degrees.  Even in the extreme case of a 75\% loss of data, the 95\% error bounds are only $\pm9.3$ degrees (Fig.~\ref{align:sens}a). Overall, although a sine fit is a relatively rudimentary model for azimuthal variations in $\rho_{HV}$, the fit proves quite precise in cases of moderate random data loss.



\subsection{Sensitivity to Contiguous Regional Data Loss }
In some cases, such as beam occultation, contiguous azimuth data points are lost.  To quantify the errors associated with these losses, the same 13,458 range cases were systematically degraded and compared with the original retrievals.  For each case, contiguous wedges ranging from 15 to 315 degrees wide were removed in 15-degree intervals.  Additionally, the position of these missing wedges were varied with respect to the original orientation angle (i.e., sine minimum).  Each pixel in Fig.~\ref{align:sens}b represents the average error across the 13,458 cases for a given sized wedge of missing data, and offset of the wedge center with respect to the original sine minimum.  Figure~\ref{align:sens}c shows the standard deviations corresponding to the mean errors shown in Fig.~\ref{align:sens}b. As expected, larger data losses generally correspond to larger retrieval errors and more deviation in the errors.  Interestingly, the position of the missing wedge with respect to the original sine minimum has a significant, systematic impact on the sign of the error (Fig.~\ref{align:sens}b).  For common sources of beam occultation, azimuthal losses typically do not exceed 15 degrees, yielding an average error magnitude of 2 degrees. For cases along a coastline where 180 contiguous degrees are missing, mean errors can range between $\pm25$ degrees, with the actual error varying erratically depending on the migrant orientation with respect to the coastline.


\begin{figure*}[t]
\centering
\includegraphics[width=6in]
  {figures/align_demoCase.pdf}
  \caption{Retrieval results over central Oklahoma (KTLX) for 29 April 2013, including the time evolution of vertical profiles of (a) reflectivity factor, (b) retrieved flight orientation, (c) percent of missing azimuth data, (d) coefficient of determination, (e) amplitude of the sinusoidal fit, and (f) 2D standard deviation over a $3\times3$ structuring element. The 0 dBZ contour is included for reference (b, dashed).  Vertical lines denote local sunset (01:13 UTC) and sunrise (11:41 UTC).} 
\label{align:fig4}
\end{figure*}


\section{Method Demonstration}
\label{align:demo}
As a simple demonstration, this method is applied to all radar scans for the Oklahoma City, Oklahoma NEXRAD (KTLX) on 29 April 2013, thereby extending the case presented in Fig.~\ref{align:fig3} through the full 24-hour period.  Figure~\ref{align:fig4} shows the time evolution of several quantities relevant to flight orientation extraction.  To aid inspection, na\"{i}ve Vertical Profiles of Reflectivity (VPRs) were constructed by linearly averaging reflectivity factor PPIs across azimuths, and mapping the resulting range profiles to height using (3).  The time evolution of VPRs (Fig.~\ref{align:fig4}a) help indicate the density and distribution of the sampled organisms.  The corresponding orientation retrievals are shown in Fig.~\ref{align:fig4}b, with the 0 dBZ isopleth included for reference. As described in Section~\ref{align:method}, several metrics for data quality are presented in Fig.~\ref{align:fig4} including (c) percent of missing azimuth data, (d) coefficient of determination, $R^2$, (e) amplitude of sine fit, and (f) 2D circular standard deviation of orientation retrievals using a $3\times3$ square structuring element.  The anomalous vertical lines in each image are cause by sun spurs, and denote local sunset (01:13 UTC) and sunrise (11:41 UTC).

The large scale ascent of nocturnal migrants following sunset is characterized by spatially and temporally continuous orientation retrievals, with organisms oriented northward.  The corresponding large $R^2$ values and sine amplitudes indicate high-quality orientation retrievals throughout the night up to approximately 2 km above ground level. At higher altitudes, the high spatiotemporal deviations in orientation indicate poor retrievals---a result that is supported by the low $R^2$ values.  Figure~\ref{align:fig4}c reveals that the inconsistent results at above 2.5 km correspond to large amounts of missing azimuth data, likely due to sparse animal density at high altitudes.  Additionally, inconsistent retrievals can be observed near ground level, possibly due to non-migratory movements or obstacle avoidance. 

Following sunrise, the airspace begins to fill with the diurnal occupants---likely insects and soaring birds.  Figure~\ref{align:fig4}b seems to indicate a coherent region of south-westward migration; however, the low $R^2$ and sine amplitudes reveal that this is not an organized movement, but rather a consistent bad retrieval.  Subsequent inspection of the underlying PPIs for this period confirm the lack of organization in the polarimetric fields, and the persistence of a sine minimum not related to coherent $\rho_{HV}$ morphology or migratory movements.  While manual inspection of images such as Fig.~\ref{align:fig4} can reveal a unique perspective on animal movements, the combination of these retrieval metrics can aid automated analysis of flight orientation by providing a method of quality control.  Within an automated framework, meteorological measurements such as winds, turbulence, temperature, boundary layer height, or cloud cover could also be compared to orientation profiles to deduce migratory flight preferences.  Overall, the ability to measure the change in orientation in time and altitude provides a wealth of new data to explore.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Radar Simulation for Biological Applications}
\label{sim}

\section{Background}
\label{sim:back}
A major difficulty in applying radar-based methods to ecology is a persistent lack of validation datasets.  Generally, the underlying biological organisms or behavioral processes that produce the resulting radar signatures are not known.  In some limited cases, optical or thermal imaging can yield information on the position and identity of organisms, but these data are generally difficult to obtain.  In an effort to understand the interactions among individual airborne organisms and the changes in their aerial and terrestrial ecosystems, Lagrangian ecological models (known as agent-based or individual-based models) have been developed to emulate the behavioral reactions to these changes.  While both ecological modeling and radar measurements have been applied in biological studies, no framework currently exists for comparing model results with this common source of biological observations.  The value of such a comparison technique has been demonstrated many times in meteorological applications in which atmospheric model output are related to weather radar measurements, providing model validation \citep{Cheong++2008,Ryzhkov++2011} or a physical interpretation of radar measurements \citep{Muschinski++1999,Jung++2008a}.  Beyond the value in validating ecological models, simulating radar signals from a set of known biological scatterers can enhance the use of radar in ecological applications.  This can include the ability to test different scanning strategies on a static distribution of scatterers, or the development biological radar retrieval algorithms using the known model state as truth. To satisfy this wide range of applications, a computational framework for synthesizing polarimetric, pulsed-Doppler, baseband radar signals has been developed.  In any of these applications, the goal of radar simulation is to emulate the forward process of radar propagation, scattering, and processing to convert a known field of scatterers into a corresponding radar dataset (Fig. \ref{sim:block}).  Following this process, the synthesized radar data can be used to infer the underlying physical scatterers. 
\begin{figure*}[t]
\centering
\includegraphics[width=4in]
  {figures/sim_block.pdf}
  \caption{Conceptual diagram of the forward (modeling) and backward (retrieval) radar measurement processes for biological sampling.} 
\label{sim:block}
\end{figure*}

\section{Notation and Coordinate Transformations}
\label{sim:notation}
For the following discussion, the term `agent' will be used to describe a modeled biological point scatterer (i.e., bird, bat, or insect) within a three-dimensional aerial domain.  As such, the instantaneous state of an agent can be described by its position, velocity, and orientation within the modeled airspace.  It is assumed that some technique (e.g., an ecological model) has produced a series of such information for a total of $N$ agents across $M$ time steps.  The ultimate goal is to relate these native model output fields to quantities that are more relevant to radar measurements.  Though conceptually simple, defining an agent's state with respect to a radar requires computations within and across several coordinate systems, often leading to notational difficulties.  In an effort to avoid ambiguities, the following defines the notation of these coordinate systems, and the sequential transformations through them.

\begin{figure*}[t]
\centering
\includegraphics[width=6in]
  {figures/sim_coords.pdf}
  \caption{Coordinate system definitions: a) The radar-centered cartesian coordinate system. b) Blowup of agent location in (a) showing velocity components. c) The radar-centered spherical coordinate system. d) Blowup of agent location in (c) showing the projection of the velocity vector onto the `ray' coordinate basis. e) Details of the agent position in (d), showing orientation angle definitions. } 
\label{sim:coords}
\end{figure*}

As a general starting point, it is assumed that the format of the ecological model output describes the geographical position of the $i^{th}$ agent at the $t^{th}$ time step as
%
\begin{equation}
\label{eq:e1}
\textbf{X}_{geo}(i,t) = \langle lon(i,t), lat(i,t), alt(i,t)\rangle
\end{equation}
%
using degrees longitude ($lon$), degrees latitude ($lat$), and altitude in meters above mean sea level ($alt$).  Similarly, the instantaneous motion of an agent is described as
%
\begin{equation}
\label{eq:e2}
\textbf{V}(i,t) = \langle u(i,t), v(i,t), w(i,t)\rangle \;\;\;\;\;\;\;\;\; [\textrm{m s}^{-1}]
\end{equation}
%
where $u$, $v$, and $w$ represent the agent's zonal, meridional, and vertical velocity components, respectively (Fig. \ref{sim:coords}b).  Thus, the ecological model output consists of two arrays--position and velocity--each with size ($N\times M\times 3$).  Starting from these native output arrays, it is first necessary to determine the desired location of the simulated radar:
%
\begin{equation}
\label{eq:e3}
\textbf{X}_{geo,rad} = \langle lon_{rad}, lat_{rad}, alt_{rad} \rangle.
\end{equation}
%
Once this radar location has been decided, a Cartesian coordinate system is defined with the radar at the origin, the positive $\hat{X}$-axis pointing east, the positive $\hat{Y}$-axis pointing north, and the positive $\hat{Z}$-axis pointing opposite the force of gravity (Fig. \ref{sim:coords}a).  The agent position array is transformed onto this radar-centered coordinate system
\begin{equation}
\label{eq:e4}
\textbf{X}_{cart}(i,t) = \langle x(i,t), y(i,t), z(i,t)\rangle\\
\end{equation}
through the geodetic relations
%
\begin{subequations}
\begin{align}
\label{eq:e5}
x &= (lon-lon_{rad}) \frac{a\pi\cos(lat)}{180^\circ[1-e^2\sin^2(lat)]^{1/2}} &[\textrm{m}]\\
y &= 111,200 (lat-lat_{rad}) &[\textrm{m}]\\
z &= alt - alt_{rad} &[\textrm{m}]
\end{align}
\end{subequations}
%
in which $a$ is the ellipsoidal Earth's major axis radius in meters and $e$ is the eccentricity \citep{Rapp1991}.  In the 1984 World Geodetic System, these values are $a = 6,378,137.0$ m and $e = 0.081819$.

Following this coordinate transformation, the Cartesian agent location array $\textbf{X}_{cart}(i,t)$ implicitly assumes that the Earth is flat.  In other words, the agent positions take the x-y plane as the Earth surface, with $z$ representing height above the radar. In this system, the radar beam path is subject to two distorting effects: an artificial upward propagation due to the curvature of the Earth, and radio refraction due to inhomogenieties in fields of pressure, temperature, and humidity \citep{DoviakZrnic}.  Within the simulation framework, it is more computationally efficient to consider radar beams that follow simple ray geometry (i.e., straight-line propagation paths).  To allow this geometry, the combined effect of these two sources of beam deformation are calculated and accounted for by `pre-distorting' the position of all agents in space with respect to the radar.    In this distorted Cartesian system, the radar beam is subject to ray geometry across a flat Earth, while retaining the same relative positions of scatterers with respect to the beam.  The process of pre-distortion is achieved in two steps, first moving the beam and each agent to remove the curvature of the Earth, and second adjusting agent positions to straightening the path of the beam (Fig.~\ref{sim:refrac}).  Hereafter the Cartesian system $\textbf{X}_{cart}(i,t)$ will refer to agent positions that have been pre-distorted to account for Earth curvature and beam refraction.

\begin{figure*}[t]
\centering
\includegraphics[width=4in]
  {figures/sim_refrac.pdf}
  \caption{Distortion of agent locations to provide ray beam geometry. (a) Original agent locations relative to a round Earth and refracting beam. (b) Adjustment of the beam path and agent location to create `flat Earth' coordinates.  (c) Adjustment of agent location to straighten beam to a ray.} 
\label{sim:refrac}
\end{figure*}

The resulting Cartesian agent location array $\textbf{X}_{cart}(i,t)$ is converted into radar-centered spherical coordinates 
%
\begin{equation}
\label{eq:e6}
\textbf{X}_{sph}(i,t) = \langle r(i,t), \phi(i,t), \theta(i,t)\rangle
\end{equation}
%
using
%
\begin{subequations}
\begin{align}
r &= \sqrt{x^2 + y^2 + z^2} &[\textrm{m}]\label{eq:e7a}\\
\phi &=  \textrm{arcsin}\left(\frac{z}{r}\right) &[\textrm{deg}]\label{eq:e7b}\\
\theta &= \textrm{arctan2}(x,y) &[\textrm{deg}]\label{eq:e7c}
\end{align}
\end{subequations}
%
with $r$ denoting the slant range between the agent and the radar in meters, $\phi$ denoting the elevation angle in degrees above the horizon, $\theta$ denoting the azimuth angle in degrees clockwise from north, and $\textrm{arctan2}(\cdot)$ denoting the four-quadrant inverse tangent.  The result of these transformations is an array of three-dimensional agent positions in the radar-relative spherical coordinates that are common in radar applications (Fig. \ref{sim:coords}c).

A similar procedure is needed to produce the final two quantities relevant to polarimetric Doppler radar measurements, namely, radial velocity and radar-relative orientation.  Because both of these quantities are dependent on the position and velocity of the agent with respect to the radar beam, it is convenient to introduce two additional coordinate bases to describe these quantities.  Following the notation presented in \cite{Bringi+2001}, the first system is centered on the agent with the positive $\hat{K}$-axis pointing radially away from the radar, the positive $\hat{H}$-axis pointing to the left of the radial vector and parallel to the ground, and the $\hat{V}$-axis completing the orthogonal, right-hand coordinate basis (Fig. \ref{sim:coords}d).  Describing this `ray' system in terms of common polarimetric radar language, $\hat{H}$ is the horizontal polarization axis, $\hat{V}$ is the vertical polarization axis, and $\hat{K}$ is the radial axis of an incident beam directed at the agent.  

The second system is centered on the agent's body, in constant alignment with the agent's orientation.  Some agent-based models may explicitly output flight orientations as agents bank, climb, and turn throughout the domain.  In this case, however, orientation is defined based on the agent's motion under the assumption that each is oriented head-first along the horizontal component of it's velocity vector with the body and both wings parallel with the horizon.  In other words, while an agent may climb or turn, the inflight orientation will not include pitch or roll---only yaw.  As a result, this Cartesian `body' system is oriented with the positive $\hat{F}$-axis pointing forward along the agent's horizontal velocity vector, the positive $\hat{L}$-axis pointing to the agent's left wing and parallel to the ground, and the $\hat{B}$-axis emanating out of the agent's back, completing the orthogonal, right-hand coordinate basis.


Following these definitions, an agent's radial velocity, $V_r$, is the projection of the native velocity vector, \textbf{V}$(i,t)$, onto the $\hat{K}$ axis, and can be calculated using the rotation relation
%
\begin{equation}
\label{eq:e8}
V_r(i,t) = -u\sin(\theta)\sin(\phi)-v\cos(\theta)\sin(\phi)+w\cos(\phi)
\end{equation}
 %
based on the agent position angles found in (\ref{eq:e7b},c).  The radar-relative orientation can be described in terms of the angular differences between the `ray' and `body' coordinate systems.  This orientation with respect to the radar beam is found using three Euler angles that define three rotations of the agent's `body' frame to a final alignment with the radar `ray' frame.  In successive order, these rotations are:

\hangindent=0.7cm $\alpha \equiv$  rotation around the $\hat{B}$-axis (i.e., heading or yaw), turning to the agent's left, and yielding the new `prime body' coordinates, $\hat{F}' \hat{L}' \hat{B}'$.

\hangindent=0.7cm $\beta \equiv$  rotation around the $\hat{L}'$-axis (i.e., elevation or pitch), with the head rising and tail dropping, and yielding the new `double-prime body' coordinates, $\hat{F}'' \hat{L}'' \hat{B}''$.

\hangindent=0.7cm $\gamma \equiv$  rotation around the $\hat{F}''$-axis (i.e., bank or roll), with the right wing dropping and left wing rising, and resulting in the alignment with the `ray' coordinates, $\hat{K} \hat{H} \hat{V}$.

\noindent
By imposing level inflight orientation, the radar-relative orientation will only require $\alpha$ and $\beta$, with $\gamma$ always assuming a value of zero.  Following the conventions shown in Fig. 1e, the these angles are calculated using
%
\begin{subequations}
\begin{align}
\alpha &= \textrm{sgn}(-u \textrm{cos}(\theta)+v \textrm{sin}(\theta))\ldots\notag\\
&\;\;\;\;\;\;\;\;\textrm{arccos}\left(\frac{xu+yv}{\sqrt{x^2+y^2}\sqrt{u^2+v^2}}\right) & \textrm{[deg]} \label{eq:e9a}\\
\beta &= -\phi & \textrm{[deg]}\label{eq:e9b}
\end{align}
\end{subequations}
%
with $\textrm{sgn}(\cdot)$ denoting the signum function.
The result of these transformations are three new arrays that describe the position $\langle r, \phi, \theta \rangle$, velocity $\langle V_r\rangle$, and orientation $\langle \alpha, \beta\rangle$ of agents in terms that are most relevant to polarimetric Doppler radar. 

\section{Radar Signal Synthesis for a Single Pulse}
\label{sim:single}
To simplify this initial formulation of the radar signal, only a single transmitted pulse will be considered.  This will eliminate the need to consider the pulse-to-pulse motion of agents and the mechanical scanning of the radar beam.  Building upon these results, Section \ref{sim:real} will formulate the full algorithm that considers beam and agent motions. 
\subsection{Defining the Radar System}
With only the position of the radar currently defined, it is necessary to specify the characteristics of the radar system.  For a single transmitted pulse, eight specifications are necessary: the radar wavelength ($\lambda$), transmit power ($P_t$), pulse width ($\tau$), receiver sample time ($\tau_s$), receiver bandwidth ($B$), system differential phase on transmit ($\psi_t$), antenna beam pattern ($f$) and gain ($G$).  Additionally, the radar pointing direction (i.e., boresight) must be specified by the antenna elevation angle ($\phi_{pt}$) and azimuth ($\theta_{pt}$).

When an existing radar system is to be emulated, it is often possible to obtain the precise specifications of the system (e.g., a measured antenna beam pattern).  However, when such information are not available, or when testing hypothetical radar designs, generalized expressions may be used.  For example, \cite{DoviakZrnic} provides a functional expression for the one-way beam weighting pattern for a circularly-symmetric beam of a given beamwidth ($\vartheta_b$) as 
%
\begin{equation}
\label{eq:e10}
f^2(\vartheta)=  \left\{ \frac{8J_2[(1.27 \pi \sin\vartheta)/\vartheta_b]} {[(1.27 \pi \sin\vartheta)/\vartheta_b]^2} \right\}^2
\end{equation}
%
where $\vartheta$ is the angular distance off of boresight in radians, and $J_2(\cdot)$ denotes a second order Bessel function.  Similarly, the range weighting function ($|W(r,r_o)|$) can be represented, assuming a Gaussian transfer function, as 
%
\begin{equation} 
\displaystyle |W(r,r_o)| = [\mbox{erf}(x+b)-\mbox{erf}(x-b)]/2
 \label{eq:e11}
\end{equation}
%
in which
%
\begin{equation} 
\begin{array}{l}
\displaystyle b = B\tau\pi/4 \sqrt{\ln2}\\
\displaystyle x = 2aB/c (r_o-r)\\
\displaystyle a = \pi/2 \sqrt{\ln2}
\end{array} 
 \label{eq:e11}
\end{equation}
%
with $r_o$ being the location of the center of the nominal range gate [m], and erf($\cdot$) denoting the error function \citep{DoviakZrnic}.  

The combination of the beam and range weighting functions defines the resolution (i.e., size) of the pulse sampling volumes, but their location in space is determined by the ranges corresponding with the receiver sample times.  In other words, immediately following the transmitted pulse, a sample will be taken every $\tau_s$ seconds, placing the center of a resolution volume every $c \tau_s/2$ meters along the ray.  This mapping between range and time provides the link between the time-series data stream that will be synthesized and the location of each nominal range gate in space, such that the $n^{th}$ sample following the transmitted pulse corresponds to a distance
%
\begin{equation} 
r_o(n) = \frac{c \; n \; \tau_s}{2}\;\;\;\;\;\;\;\;\;\;\; \textrm{[m]}
 \label{eq:e12}
\end{equation}
%
with $c$ representing the speed of light [m s$^{-1}$]. Combined with the two boresight angles ($\phi_{pt},\theta_{pt}$), $r_o(n)$ provides the spatial coordinates for the center of each resolution volume for the $n$ time-series samples following the transmitted pulse. 

\subsection{Calculating Echo Amplitude and Phase}
The defining characteristic of polarimetric radar is the transmission and reception of radiation of diverse polarizations---in this case, horizontal and vertical.  The method of sequestering these signals upon transmission and reception can be conducted simultaneously or alternating \citep{DoviakZrnic}.  The following formulation will consider Simultaneous Transmit And Receive (STAR) operation, but can be extended to Alternating Transmit and Simultaneous Receive (ATSR) or Alternating Transmit and Alternating Receive (ATAR) modes by modifying the transmission and reception calculations for each polarization.

Using the radar system definitions above, the power contribution from the $i^{th}$ agent for the $n^{th}$ sample following a single pulse is calculated using the radar range equation for point scatterers,
%
\begin{equation} 
P_{r,i}(n) = \frac{P_t G^2 \lambda^2 \sigma_{b,i} f^4(\vartheta_i) |W^2(r_i,r_o)| } {(4 \pi)^3 r_i^4}\;\;\;\;\;\; \textrm{[W]}
 \label{eq:e13}
\end{equation}
%
in which $\vartheta_i$ is the angular distance of the $i^{th}$ scatterer off the center of the beam axis, $r_i$ is the range to the $i^{th}$ scatterer, and $\sigma_{b,i}(\alpha,\beta)$ is the radar cross-section at the given polarization, corresponding with the agent's current orientation (discussed in Section \ref{sim:single:scattermod}) ~\citep{Skolnik}.  The echo phase from the $i^{th}$ agent is defined as
%
\begin{equation} 
\psi_i = \frac{4 \pi r_i}{\lambda} + \frac{4 \pi T_s}{\lambda}V_{r,i} + \psi_{s,i} + \psi_t \;\;\;\;\;\; \textrm{[rad]}
 \label{eq:e14}
\end{equation}
%
using the initial transmitted phase ($\psi_t$), as well as the range ($r_i$), radial velocity ($V_{r,i}$), and scatter phase ($\psi_{s,i}(\alpha,\beta)$) of the $i^{th}$ agent ~\citep{DoviakZrnic}.  Some potential treatments for calculating or omitting $\psi_{s,i}(\alpha,\beta)$ are discussed in Section \ref{sim:single:scattermod}.

With each agent's contribution to the echo power and phase calculated for each sample, the complex echo voltage contribution from the $i^{th}$ agent for the $n^{th}$ sample is calculated using the in-phase and quadrature-phase components,
%
\begin{equation}
\begin{array}{l}
\displaystyle I_i(n) = \sqrt{\frac{P_{r,i}(n)}{2}}\cos{\psi_i} \\
\displaystyle Q_i(n) = \sqrt{\frac{P_{r,i}(n)}{2}}\sin{\psi_i} \\
\displaystyle V_i(n) = I_i(n) + \jmath Q_i(n)
\end{array} 
\label{eq:iq}
\end{equation}
%
with $\jmath$ representing the imaginary unit (i.e., $\sqrt{-1}$).
%
Next, the time-series components from each of the $N$ agents are coherently summed to obtain the final radar time-series echo voltage for the given polarization:
\begin{equation}
 V(n) = \sum_{i=1}^{N} V_i(n).
 \label{eq:sum}
\end{equation}
%
For ATAR operation, this process is calculated twice---once using $\sigma_{b,hh}$ and once with $\sigma_{b,vv}$.  For STAR, these calculations must be computed for cross-polar contributions as well.  The final dual-polarization time-series for STAR operation is defined as the coherent sum of co- and cross-polar contributions,
\begin{subequations}
\begin{align}
V_{h}(n) = V_{hh}(n) + V_{hv}(n) \\
V_{v}(n) = V_{vv}(n) + V_{vh}(n)
\label{eq:star}
\end{align}
\end{subequations}
in which the phase offsets have been determined by the value of $\psi_t$ and $\psi_{s,i}$ in equation \ref{eq:e14}.  Finally, if desired, white noise can be added to the final time-series at a level defined by the user.  The resulting synthesized time-series segment represents the range-time samples from a single pulse along a ray.  

\subsection{Scattering Models}
\label{sim:single:scattermod}
Although the agent behavioral model may tailor the motion of agents to mimic a desired species, a method is still needed to describe the physical characteristics of these point scatterers. More specifically, a method for determining the radar cross section (RCS) and scatter phase for each agent is required. In general, the RCS of an agent will vary as a function of size, shape, orientation, and composition. It is possible that it may even change in time (e.g. the modulation in RCS as a bird flaps its wings). Depending on the application, it may not be necessary to include all of these considerations. For example, many studies using only radar reflectivity factor may assume agents to be dielectric spheres of a constant size, while a study of wing beat patterns would require that these spheres oscillate in size over time. Furthermore, a study on biological polarimetric signatures would certainly require non-spherical RCS definitions.

One challenge within this computational framework is determining the scattering characteristics of agents as a function of look angle and polarization.  Namely, the radar cross section and scatter phases must be defined at the wavelength that is being simulated.  One possibility is to use methods such as T-matrix or Rayleigh-Gans to calculate approximations for these values \citep{Bringi+2001}.  Other techniques can use measurements (such as this in Chapter \ref{measure}) to define scattering characteristics as a lookup table of radar cross sections and phases.  In either case, the quality of the synthesized radar moments and polarimetric products is only as good as the underlying scattering models. With this modular design, it is possible to test the same biological behavioral simulations with different scattering models (e.g., spheres versus spheroids) to investigate the differences in the observed radar moments.  


\section{Radar Signals for Realistic Sampling}
\label{sim:real}
Section \ref{sim:single} described the synthesis of a single pulse, but typical scanning strategies require the transmission and reception of many pulses, often directed at different regions of space.  In this case, several more radar parameters must be specified that define the Volume Coverage Pattern (VCP).  These parameters include the set of azimuth and elevation angles that the antenna will cover, the Pulse Repetition Time ($PRT$), and the antenna scanning rate ($\omega_r$).  Additionally, the pulse-to-pulse motion of the agents must be considered as the radar samples the airspace.  The agent-based model developed for the following examples has a template resolution of 1 sec, and is linearly interpolated to result in unique data at each PRT.  The resulting effects of agent motions within the beam, coupled with the pulse-to-pulse beam motions, are samples that vary in time at the PRT yielding realistic Doppler spectra.

The process of synthesizing realistic radar data within this framework is primarily an exercise in bookkeeping---solving a relatively simple set of equations for several hundred-thousand agents at four unique polarization combinations, and coherently combining the results in a specific recipe.  It is also simple bookkeeping to keep track of the agent locations, velocities, and orientations with respect to the scanning beam boresite at every pulse.  While conceptually easy, the process is computationally expensive, requiring large arrays to store the pertinent information.  



\section{Examples}
\label{sim:examples}
\subsection{Agent-based Behavioral Model}
\label{sim:abm}
The first step before any radar simulation can be performed is producing a set of LaGrangian agents that will act as the scatterers and move through the radar domain.  The behavioral simulator is based on the `boid' model described by \cite{Reynolds1987,Reynolds1999}, and determines an agent's flight velocity by the position of the surrounding individuals.  From this method, a relatively simple set of behavioral attributes, or `rules', defines the decision-making process of each agent as it moves through space and time.  The use of a rule-based technique enables dynamical interactions to emerge across a wide range of spatial and temporal scales.  Three major rules control the behavioral dynamics of the group: collision avoidance, group cohesion, and flock alignment.  To emulate the behavior of a desired organism, the parameters controlling these characteristics must be tuned.  These parameters include the average agent velocity, minimum spacing between individuals, angular field of view for group cohesion ($\varphi_{co}$), radius of influence for group cohesion ($r_{co}$), angular field of view for flock alignment ($\varphi_{al}$), radius of influence for flock alignment ($r_{al}$), initial position, and initial velocity of each agent.  Figure \ref{sim:agnts}a shows an example of these regions for the current agent (black triangle) with respect to the surrounding agents (gray triangles).  The agent symbols are aligned along their current velocity vectors.  Based on these inputs, several velocity vectors may be calculated.  A cohesion vector ($v_{co}$) points toward the center of mass of the individuals within the cohesion field of view and radius of influence of the active agent, and tends to draw the agent into the group.  A separation vector ($v_{sp}$) points away from the agent's nearest neighbor and promotes collision avoidance.  An alignment vector ($v_{al}$) is the averaged velocity vectors of individuals within an agent's alignment field of view and radius of influence.  An additional vector can limit acceleration by being set to the velocity of the previous time step ($v_{(t-1)}$).  Figure \ref{sim:agnts}b shows these resulting vector components for the agents and regions in Fig. \ref{sim:agnts}a.  The final velocity vector for each agent, $v(t)$, is the weighted sum of these individual components,
%
\begin{equation} 
v(t) = \frac{w_1 v_{co} +  w_2 v_{sp} +  w_3 v_{al} +  w_4 v_{(t-1)}} {w_1+w_2+w_3+w_4}.
 \label{eq:velWt}
\end{equation}
%
\begin{figure*}[t]
\centering
\includegraphics[width=6in]
  {figures/sim_agnts.pdf}
  \caption{Geometry of agent decision-making process. (a) The agent's regions of influence for cohesion and alignment. (b) The resulting motion vectors for cohesion, alignment, separation, and momentum.} 
\label{sim:agnts}
\end{figure*}
The weights, $w_n$, must be chosen to best represent the desired behavior of the group.  For example, a migratory flock of birds would require a large weight on the alignment component to ensure the flock travels in the same direction.  Conversely, simulating bees swarming around a hive would require very little directional consistency, and therefore a low weight would be placed on alignment.  Additional rules may be developed \emph{ad hoc} to evoke more specialized behaviors.  As the algorithm iterates over time, the calculation in (\ref{eq:velWt}) is performed for each agent.  The final result is an array of the three-dimensional coordinates of the location and velocity of each agent for every time step, with a temporal resolution on the order of one second.  The resulting group behavior exhibits dynamic interactions across a range of scales as agents adapt their motion to deal with their surroundings.


\begin{figure*}[h!]
\centering
\includegraphics[width=6in]
  {figures/sim_steps.pdf}
  \caption{A simplified simulation case illustrating transformation from agent-based model to radar power product. (a) A 2D agent based model in which triangles indicate the agent position and direction of motion. A sample resolution volume is highlighted in red.  (b) the co-polar horizontal time-series components for the two agents in the highlighted resolution volume. (c) The co-polar time-series for the highlighted resolution volume created from the coherent addition of the agent time-series contributions, and the corresponding Doppler spectrum. (d) The resulting PPI of log power for this sweep.} 
\label{sim:steps}
\end{figure*}
\clearpage

\subsection{Test Cases}
\label{sim:cases}
The first application of the biological model is meant to illustrate the process of radar simulation in the simplest form.  A number of flying agents are defined to move at approximately 10 m s$^{-1}$  and are confined to a 2D plane in space.  As the agents move across the plane, they are scanned by a radar with a beam width of $11.25^\circ$ and elevation angle of $0^\circ$ located in the same plane (Fig. \ref{sim:steps}a).  The isosceles triangles in Fig. \ref{sim:steps}a denote the positions of agents, with the acute angle pointing in the velocity direction.  The resolution volume highlighted in red in Fig. \ref{sim:steps}a contains two agents moving away from the radar.  Figure \ref{sim:steps}b shows the time-series data for the horizontal polarization for each agent over the 340 samples in the azimuth.  Coherently summing the time-series for each agent yields the time-series for the resolution volume (Fig. \ref{sim:steps}c, top) which can be decomposed into the Doppler spectrum (Fig. \ref{sim:steps}c, bottom).  The Doppler spectrum indicates two modes closely spaced near $+10$ m s$^{-1}$, indicating the radial velocity of the agents.  Following spectral processing, the corresponding radar products can be displayed (e.g., echo power, Fig. \ref{sim:steps}d).  Even the most complicated simulations follow this basic process, although the number of agents or polarization contributions can be larger.


The second test extends the simulation to three dimensions---both in agent motions and scanning pattern.  In this case, a population of swarming agents are confined to an imaginary 1-km cube that is located several kilometers to the southeast of a simulated NEXRAD site (Fig. \ref{sim:dumbDemo}a).  As the agents swarm throughout the domain, the NEXRAD scans in VCP 32, resulting in the corresponding reflectivity factor PPIs that update at approximately 10-minute resolution (Fig. \ref{sim:dumbDemo}b).  The comparison of agent positions to reflectivity PPI demonstrate the effects of inter-scan agent motions; the positions of the agents were moving throughout the 10-minute scan.


\begin{figure*}[t]
\centering
\includegraphics[width=6in]
  {figures/sim_dumbDemo.pdf}
  \caption{Screenshot of a simple test case of 3D swarming agents confined to an imaginary cube and the resulting radar reflectivity factor.} 
\label{sim:dumbDemo}
\end{figure*}


\subsection{Simulating Bat Emergences}
\label{sim:batEmer}
To demonstrate the utility of this method in a realistic application, the nocturnal dispersion of a colony of Brazilian free-tailed bats from a cave roost was simulated, and this emergence event was virtually sampled by an S-band weather radar. The inputs to the biological behavioral model were based on the BATOIDS agent-based model, in which the primary behavioral rules include velocity matching and collision avoidance \citep{Hallam++2006}. The total population of the colony was set to one hundred thousand agents. While this value is representative of many maternity roosts of free-tailed bats in central Texas, some colonies may have populations on the order of one million \citep{Betke++2008}. The model was initialized with all bats within the cave (below ground level) and they were allowed to exit the cave at a set rate through a finite-sized hole at the surface. This produced the columnar group formation characteristic of these emergences \citep{Wilkins1989}. Additional rules were set to mimic the transition from emergence to dispersion flight modes at a set altitude. The biological behavioral simulation presented in this example spans a full emergence event, lasting one hour. Fig. \ref{sim:emerge} shows four snapshots from the first twenty minutes of this simulation, with time increasing from left to right. The top row of panels shows a projection of the location of all 100,000 agents onto a horizontal plane. The bottom row shows the projection of the agents onto a vertically oriented plane running west to east. Initially, the bats exit the cave and gain altitude in a dense, columnar formation, reducing the individual risk of predation. The risk of encountering aerial predators decreases away from the cave mouth, allowing the bats to transition to a dispersion flight mode. As the emergent column continues the transition, the horizontal projection begins to develop a hole in the center of the group as bats diverge from the cave location. At the final stages of the emergence, the bats descend to the height at which they forage for insects \citep{Wilkins1989}. While this visualization only captures the overall group motion, it is important to note that complex behavioral dynamics are also occurring on much smaller scales as individuals interact to avoid collisions and modify their position within the group.

\begin{figure*}[t]
\centering
\includegraphics[width=6in]
  {figures/stepa1.pdf}
  \caption{Biological behavioral model snapshots from 5, 10, 15, and 20 minutes into the simulation (left to right). The top row shows the horizontal projection of agent locations onto the ground. The bottom row shows the projection of agent locations onto a vertically-oriented plane running west to east. Each point represents the location of a single agent (i.e., bat).} 
\label{sim:emerge}
\end{figure*}

A simple scattering model was chosen for this example to aid in the interpretation of the radar moments. Each agent is represented as a liquid water sphere having a mass of ten grams, typical of an adult Brazilian free-tailed bat. The RCS is calculated using the full Mie solution \citep{Mie1908}, yielding a $\sigma_{b,i}$ of 10.59 cm$^2$. The use of a constant sized, spherical model avoids the resonant and angular effects characteristic of real-life bioscatter, and allows a straightforward understanding of the radar simulation output. The radar system configuration was set to emulate a standard S-band, NEXRAD weather surveillance radar running in VCP 32, the standard clear-air mode \citep{Serafin+2000}. The location of this simulated radar with respect to the cave mouth was chosen to mimic the location of an actual NEXRAD site in Del Rio, Texas (KDFX) which overlooks a known cave of Brazilian free-tailed bats (Frio Cave). Because KDFX regularly detects the nightly emergence from Frio Cave, an abundance of real-life data is available for comparisons with the simulation output. Fig. \ref{sim:frio} shows the 0.5$^\circ$ plan position indicator (PPI) of reflectivity, velocity, and spectrum width from KDFX on May 16, 2011 for 1:02:20 UTC (left) and 1:12:27 UTC (right) as the Frio Cave colony emerges after sunset. The output from the biological behavioral model was ingested into the radar driver for the specified radar configuration, and the time-series signals were synthesized as described. The resulting baseband signals were processed using the auto-covariance pulse-pair processor to calculate the radar moments \citep{DoviakZrnic}. The radar reflectivity factor (Z) was calibrated and range corrected by a factor of $r^2$. Fig. \ref{sim:simfrio} shows the resulting 0.5$^\circ$ PPIs from the simulated NEXRAD from approximately 30 minutes and 40 minutes into the simulation, representing two consecutive volume scans from VCP 32. In addition, the instantaneous positions of the sampled agents are shown, providing the �ground truth� reference (Fig. \ref{sim:simfrio}, bottom row). The spatial domain and color scales for Figs. \ref{sim:frio} and \ref{sim:simfrio} are similar, providing direct comparison.


By comparing the simulation results to the KDFX data, it is possible to reach various conclusions about the characteristics of actual bats with respect to the behavioral simulations. For example, the radar reflectivity factor in the simulation is found to be consistently higher than what is observed on KDFX. There are at least two explanations which could lead to this discrepancy. It is possible that the number concentrations in the simulation were too high as the bats dispersed, suggesting that they may travel in less dense groups than were simulated. Alternatively, it may be that the RCS used in the simulation was too high. Considering that the RCS was obtained using the equi-mass, liquid water sphere assumption, it could be concluded that a portion of a bat�s mass is less reflective than water. It is also possible that the spherical geometry enhances the backscatter as compared to the complex shape of an actual bat. Similar comparisons can be made between velocity and spectrum width estimates from the two data sets. For example, the velocity data resulting from the simulation study are generally larger than those observed on KDFX. It could either be that the simulated agents are flying too quickly, or that real bats have more directional variation, leading to lower measured velocity magnitudes. Because the simulated spectrum width agrees with the KDFX observations, it is likely that the directional variation introduced in the behavioral simulation is correct, and indeed the agent velocities are simply too high. Taken collectively, results from the simulator when used in conjunction with actual radar observations of birds and bats could improve our ability to use radar data to functionally estimate the distribution of aerial individuals across the observable landscape and to serve as empirical estimates of behavioral models and utilization distributions of assemblages of aerial species. Utilization distributions are bivariate probability distribution functions that describe the locations of animals (x, y) during their normal activities of foraging, mating, and caring for their young \citep{VanWinkle1975}. Typically utilization distributions are constructed statistically based on the observed behaviors of one of more selected individuals over time. Consequently, the calculations may require extensive data sets of animal tracking records and the resulting modeled distribution functions are not available to researchers in real time. If properly diagnosed, radar observations can provide direct empirical estimates of utilization distributions of foraging groups of animals, such as bats, and be used to refine agent-based models of their behavior. Consider, for example, the simulated results shown in Fig. \ref{sim:simfrio}. The images of radar reflectivity factor shown in the upper panels can ostensibly be used as a proxy for the density and distribution of the simulated bats represented in the lower panels. The simulator will serve to enable needed comparisons of radar output with actual underlying collections of aerial individuals to determine the extent to which utilization distribution functions can be extracted from such observations. Being able to directly estimate utilization distributions in quasi real time would provide a powerful means of exploring the plasticity in individual positioning within foraging groups and determine whether landscape heterogeneity influences group cohesion during foraging. Such information is needed in order to fulfill the demands of ecological forecasting in an era of rapid climate and environmental change \citep{Clark++2001}.

\begin{figure*}[t]
\centering
\includegraphics[width=6in]
  {figures/stepa2.pdf}
  \caption{The emergence of a Brazilian free-tailed bat colony, Frio Cave, as seen by the Del Rio, Texas NEXRAD (KDFX) at 1:02:20 UTC (left column) and 1:12:27 UTC (right column). Range rings are at 10 km intervals, centered on the radar location. Moments include radar reflectivity factor (top), radial velocity (middle), and spectrum width (bottom).} 
\label{sim:frio}
\end{figure*}

\begin{figure*}[t]
\centering
\includegraphics[width=6in]
 {figures/stepa3.pdf}
  \caption{The emergence of a Brazilian free-tailed bat colony synthesized by the radar simulator to emulate KDFX in VCP32 at 30 minutes (left column) and 40 minutes (right column) into the biological simulation. Range rings are at 10 km intervals, centered on the radar location. Moments include radar reflectivity factor (top), radial velocity (top middle), and spectrum width (bottom middle). The corresponding �ground truth� agent locations are shown as black points (bottom).} 
\label{sim:simfrio}
\end{figure*}


Moreover, the simulator can be used to assess the impacts of detectability issues related to a particular species in conjunction with radar coverage, sampling effects, and other factors. As an example, the emergence presented in Fig. \ref{sim:emerge} was used to simulate signals from a rapid scanning mobile radar (Fig. \ref{sim:mobile}).  Because the underlying biological model is identical, direct comparisons can be made between these data and those from the NEXRAD simulation.  



\begin{figure*}[t]
\centering
\includegraphics[width=6in]
 {figures/sim_mobile.pdf}
  \caption{The emergence of a Brazilian free-tailed bat colony synthesized by the radar simulator to emulate RaXpol at 10 minutes (left column) and 30 minutes (right column) into the biological simulation. } 
\label{sim:mobile}
\end{figure*}




%\begin{figure*}[t]
%\centering
%\includegraphics[width=6in]
%  {figures/sim_batPPI.pdf}
%  \caption{. } 
%\label{sim:batPPI}
%\end{figure*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Radar Validation using Acoustic Localization}
\label{local}
\section{Background}
\label{local:back}
Measurements of abundance, distribution, and composition of nocturnal migrants are valuable for both basic science and applied problems in ornithology. Of the existing techniques for monitoring avian life in the airspace, few can singlehandedly provide information on all three of these metrics simultaneously. In some applications, however, acoustics can provide an exception to this general limitation.  This ability can enable validation of polarimetric measurements, or provide collocated, species-specific observations that can form the basis of polarimetric training data for machine learning applications.

Acoustic observations of nocturnal flight calls have long been a source of information on the presence and identity of birds in the airspace \citep{Libby1899, Farnsworth2005}. The development of amplification and recording devices propelled acoustic methods into regular use in avian field studies \citep{Graber+1959}. Some applications include the use of acoustic proxies for abundance \citep{Farnsworth++2004}, as well as regionally distributed recording stations for broad scale distribution studies \citep{Evans+1999}. Recent advances in wireless electronics and digital recording have resulted in sophisticated audio processing techniques \citep{Blumstein++2011}, including automated call detection \citep{Raven2011}, recognition \citep{Kogan+1998, Cortopassi+2000, Baker+2003}, and localization.

Acoustic localization (sometimes referred to as �triangulation�) is the process of identifying the source location of sounds using recordings from multiple time-synchronized microphones \citep{Blumstein++2011}. Bioacoustic localization of birds has been developed theoretically \citep{Magyar++1978, Spiesberger2001, Spiesberger2005} and demonstrated in several field trials. While the utility of these techniques has been proven, the existing applications have been limited exclusively to terrestrial or near-terrestrial environments.  For example, \cite{Magyar++1978} retrieved the ground position of a Bobwhite Quail (\emph{Colinus virginianus}).  \cite{McGregor++1997} focused on retrievals of several synthesized and natural calls, all located between 1 and 1.5 meters above ground level.  The Red-winged Blackbirds (\emph{Agelaius phoeniceus}) studied by \cite{Spiesberger1999} were located in small trees below 3 meters, and similarly, the duetting calls broadcasted by \cite{Mennill++2006} represented birds perched 1 meter above ground level. \cite{Collier++2010} chose the Mexican Antthrush (\emph{Formicarius moniliger}) as a model species for their localization studies for its ground-based vocalizations that facilitate two-dimensional localizations. The demonstration by \cite{Mennill++2012} broadcasted all calls from 1.5 AGL, and both \cite{Wang++2005} and \cite{Wilson++2013} describe all calls and retrievals as occurring in the same horizontal plane.  Most recently, \cite{Frommolt+2014} performed a multi-year study on the Eurasian Bittern (\emph{Botaurus stellaris}) using all two-dimensional retrievals.  Overall, no published studies have localized birdcalls above 13.5 meters \citep{Ali++2009}, or in flight within the airspace. Extension of bioacoustic localization to three-dimensional space provides surveillance of calling nocturnal migrants that can be collocated with polarimetric radar measurements. With this application in mind, an acoustic microphone array has been designed and constructed at the University of Oklahoma as a flight call localization proof-of-concept testbed. 

\section{Methods and Materials}
\label{local:method}
\subsection{Computational Methods using Time Difference of Arrival}
\label{local:method:TDOA}
Several techniques exist for extracting sound sources from multiple recordings \citep{Blumstein++2011}.  Many of these techniques have been developed for applications ranging from acoustic aircraft surveillance \citep{Blumrich+2000} to enemy gunshot positioning \citep{Ferguson++2002}.  The Time Difference of Arrival (TDOA) method has been successfully transitioned to a number of biological applications, including monitoring marine and terrestrial wildlife \citep{Magyar++1978, Clark+2000, Collier++2010}.  The fundamental TDOA technique was developed for radio navigation in the early 1970�s \citep{VanEtten1970, Schmidt1972}, and has been subsequently applied to several bioacoustic software packages [Raven Pro (Cornel Lab of Ornithology, Ithaca, New York, USA); Avisoft-SASLab Pro (Avisoft Bioacoustics, Berlin, Germany); SIGNAL (Engineering Design, Belmont, Massachusetts, USA); ArrayGUI (J. Burt, Seattle, Washington, USA); Sound Finder \citep{Wilson++2013}].  Although these software packages exist, we chose to write our own TDOA implementation in MATLAB (2010b) to provide easier access and modification to the underlying algorithm. The TDOA workflow is as follows:
\begin{enumerate}
\item Record six synchronized channels of audio from the acoustic array (detailed in the following section).
\item Use a single band-limited energy detector in Raven Pro 1.4 \citep{Raven2011} to identify potential calls within the audio file.
\item Manually screen the automated selections to eliminate false alarms and ensure that each call is detected on at least five of the six channels.
\item Calculate the temporal cross-correlation of the audio waveforms from each channel to obtain the arrival time lags \citep{Spiesberger+1990}. 
\item Calculate the sound source location from the five time lags using the set of equations presented by \citep{Spiesberger2001}.
\item When detections occur in all six channels, perform Step 5 on each of the possible subsets of five channels and take the median retrieval location. 
\end{enumerate}
The results that are presented through the duration of this chapter were obtained using this workflow, but it is anticipated that similar results would be obtained using any of the listed software packages that support 3D TDOA retrievals (e.g., Sound Finder).




\subsection{An Array Design Facilitating Three-dimensional Retrievals}
\label{local:method:array}
The basic hardware requirement for 3D TDOA localization is a distributed network of five or more time synchronized recording devices \citep{Spiesberger2001}; however, it is the placement of these microphones combined with the recorder sample rate that determines whether practical 3D localization can be achieved. To demonstrate this dependence, consider a vertical tower with a microphone (M1) located at the base, and a second microphone (M2) located 154 cm directly above M1 (Fig.~\ref{local:err}). We will call this separation distance between the microphones $d$.  Both microphones are synchronized and recording at a typical rate of 22050 samples per second (i.e., $\tau = 22050$ Hz), and the atmospheric speed of sound, $v$, is 340 m~s$^{-1}$.  In this case, the distance corresponding to each recorded sample is $\delta = v/\tau = 15.4$ cm.  
\begin{figure}[t]
\centering
\includegraphics[width=4in]
  {figures/local_err.pdf}
  \caption{Schematic of all possible hyperboloids for $\ell_{max}=10$.  The red circle indicates the sound source location of a calling bird.} 
\label{local:err}
\end{figure}
When a flying bird calls (Fig.~\ref{local:err}, red circle), the sound will eventually arrive at both microphones and the offset number of samples, or lag, between the arrivals can be computed. The maximum possible lag, $\ell_{max}$, will occur when the call is located on the line emanating out of the top of the tower (Fig.~\ref{local:err}, green line), and is equal to the maximum lag samples that fit between M1 and M2.  That is, $\ell_{max} = d/\delta = 10$ samples.  Of course, the minimum possible lag is zero, and will occur when the call is located on the plane passing between the two microphones (Fig.~\ref{local:err}, blue line).  As a result, there are only 11 possible lags that can occur for a bird calling: $\ell = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10$.  Each of lags 1 through 9 creates a unique hyperboloid passing between the vertical line and the horizontal plane (Fig.~\ref{local:err}, black curves).  The more black curves that exist, the greater number of possible localization outcomes, and therefore, higher retrieval accuracy.  The number of black curves will always equal ($\ell_{max}$-1), so there are only two ways to increase accuracy: increase the audio sampling rate or increase the distance between microphones.  In this example, no hyperboloid passes through the location of the bird, and so the retrieval must choose one of its neighboring hyperboloids.  This necessary deviation from the true bird location results in retrieval error.  By doubling the distance between microphones ($d = 308$ cm), $\ell_{max}$ will increase to 20 samples and the number of hyperboloid solutions will double, effectively placing an additional black curve between each existing one, and decreasing the error in the bird location solution. 


In short, accurate 3D localization requires sufficient microphone height diversity in the array layout. Many microphone arrays are distributed with all microphones at the same or similar height, resulting in extreme retrieval uncertainty in height.  With this in mind, we achieved vertical microphone separations using three 30-foot tall towers. Each tower was constructed from three connected 10-foot segments of black iron pipe using standard pipe couplings, and was held upright by several guy wires (Fig.~\ref{local:setup}a).  Microphones were secured to the top and bottom of each tower, with towers arranged in an equilateral triangle with vertices 20 m apart (Fig.~\ref{local:setup}b).  Rather than placing the lower microphones directly on the ground, they were secured approximately 1.5 m high on the tower to inhibit gnawing or nesting rodents, as well as mitigate noise from insects on the ground.  In this configuration, we were able to raise the towers with only two people; however, we strongly suggest a minimum of three to avoid perilous results.
%
\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/local_setup.pdf}
  \caption{Acoustic array setup and components. (a) Photo of array deployment beside radars and lidar at the ARM SGP site in spring 2014. (b) Schematic of array layout. (c) Close-up on one microphone enclosure with protective cloth cover. (g) Inside of microphone enclosure revealing foam baffling surrounding flowerpot microphone. (e) Central enclosures holding amplifiers (bottom) and laptop (top). (d) Inside of amplification enclosure. (f) Laptop for data acquisition and storage.} 
\label{local:setup}
\end{figure}
%

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/local_kite.pdf}
  \caption{Validation using helikite and test samples. (a) helikite with attached GPS recorder, mp3 player, and speakers. (b-l) Spectrograms for test sample recordings: (b) Black-throated Blue Warbler, (c) Dickcissel, (d) Indigo Bunting, (e) Ovenbird, (f) Summer Tanager, (g) Swainson's Thrush, (h) Vesper Sparrow, (i) Wood Thrush, (j) Yellow-billed Cuckoo, (k) Yellow Warbler, and (l) synthetic signal.} 
\label{local:kite}
\end{figure}



\section{Test using Kite-lofted Speaker and GPS}
\label{local:kite}
With the array constructed, our first test was to characterize the errors in call localizations aloft.  This verification required some method of generating calls at known locations high above ground level.  Our solution was to lift a small speaker, mp3 player, and GPS unit using a helium balloon-kite hybrid (helikite, Fig.~\ref{local:kite}).  A loop of ten varying flight call recordings and one synthetic tone sweep was broadcasted throughout different positions in the airspace by raising, lowering, and walking with the helikite tether line (Fig.~\ref{local:kite}).  The ten flight calls from \cite{Evans+2002} were chosen to cover a wide range of frequencies, durations, and bandwidths, and are detailed in Figure~\ref{local:kite}.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/local_scatter.pdf}
  \caption{Comparison of localization results to GPS measurements.  (left) Longitudinal retrieval comparisons.  (center) Latitudinal retrieval comparisons.  (right) Altitudinal retrieval comparisons in height above ground level.  The solid line denotes the one-to-one boundary.  The region bounded by the dashed lines indicates the reported measurement uncertainty of the GPS unit (�5 m).} 
\label{local:scatter}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/local_box.pdf}
  \caption{Call-specific errors. (left) Longitudinal retrieval errors.  (center) Latitudinal retrieval errors.  (right) Altitudinal retrieval errors.  Red dots indicate outliers. The region bounded by the dashed lines indicates the reported measurement uncertainty of the GPS unit (�5 m).  Four-letter alpha codes correspond to calls listed in Fig.~\ref{local:kite}.} 
\label{local:box}
\end{figure}

During the course of the experiment, the helikite attained a maximum altitude of 140 m above ground level, and maximum total distance of 175 m from the array center.  Unfortunately, due to the small, lightweight speaker, we could only consistently detect calls up to approximately 50 meters above ground level before the signal extinguished into the ambient noise. The localization algorithm was run on all recorded calls having a signal-to-noise ratio greater than 10 decibels, and compared to the GPS measurements (Fig.~\ref{local:scatter}).  Because the GPS unit reports with 5-meter accuracy, localization results within 5 meters of the GPS are considered perfect retrievals.  Additionally, these retrieval errors are summarized in terms of call-specific variations (Fig.~\ref{local:box}).

Comparison with the GPS reveals high localization retrieval accuracy for all of the detected calls (Fig. \ref{local:scatter}).  Similar studies have noted that artificially broadcasted calls can be much lower in amplitude than those emitted from actual birds \citep{Ali++2009}, and our need to keep the speaker light exacerbated this effect.  Nonetheless, we believe that the errors recorded within our height range are still characteristic of potential retrievals higher aloft provided similar atmospheric conditions.  Admittedly, the only way to prove this accuracy at higher altitudes would require a more powerful speaker that can replicate true call amplitudes and a much larger balloon.

Considering species-specific errors, it is clear that some calls yielded consistently poor retrievals (Fig. \ref{local:box}).  A dominant factor in these poor retrievals is the frequency of the underlying call.  The atmosphere acts as a low-pass acoustic filter, and so low frequency calls should attenuate the least along their path \citep{Horton++2015att}. As a result, we would expect that low frequency calls should be the most detectable; however, ambient noise is also attenuated at higher frequencies, resulting in greater noise amplitudes at lower frequencies. The Yellow-billed Cuckoo call resides in this elevated noise region (see Fig. \ref{local:kite}j), resulting in the worst retrieval errors (Fig. \ref{local:box}, YBCU).  Conversely, the impulse-like call of the Black-throated Blue Warbler (Fig. \ref{local:kite}b) is high enough in frequency to avoid the elevated noise levels, yielding exceptionally good retrievals (Fig. \ref{local:box}, BTBW). 

Several additions to the retrieval technique can enhance the overall method.  One addition would be the use of acoustic self-surveys as by \cite{Collier++2010}.  By periodically transmitting an acoustic impulse from a known location, the exact microphone positions can be regularly surveyed to yield better retrievals.  This process would be especially beneficial in long-term field deployments when microphone locations may slowly change in time.  For example, as guy-wires gradually stretch and are retightened, towers can lean slightly off vertical, resulting in horizontal changes in microphone locations---especially at the top of the tower.  Regular acoustic self-surveys can mitigate this effect.  

At higher altitudes, the propagation of bird calls will have a greater atmospheric dependence \citep{Horton++2015att}.  Factors influencing propagation of calls include the variable speed of sound in regions of vertical temperature gradients, and call drift from winds.  Generally these local meteorological measurements will not be available, motivating retrieval techniques that can account for these effects.  Spiesberger (1999 and 2005) demonstrates two methods that can solve for atmospheric conditions as well as call sources.

%%%%%%%%%%%%%%%  END NEEDS EDITING %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Field Demonstration at the ARM SGP Site}
\label{local:demo}
Acoustic flight call recordings can be compared to other remote sensing measurements such as radar or thermal images to better characterize animals in the airspace \citep{Larkin++2002, Farnsworth++2004, Horton++2015}.  It is generally the case, however, that a recorded call cannot be directly attributed to a specific animal in other observations.  For example, a flight call may be recorded while several birds are observed flying overhead, but it is usually unclear which bird uttered the call.  Localization of the calls can solve this problem by providing the source position of the sound.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/local_killdeer.pdf}
  \caption{A 40-sec audio spectrogram of a Killdeer (Charadrius vociferus, pictured) calling continuously while flying past the acoustic recording array location on 11 April 2014.  Brighter shades denote louder acoustic amplitudes.  The zoomed region spans approximately seven seconds, and shows six successive Killdeer calls.} 
\label{local:kill}
\end{figure}

As a first attempt to connect flight call localization retrievals to radar observations, we organized a field campaign in spring of 2014 at the Atmospheric Radiation Measurement (ARM) program�s Southern Great Plains (SGP) climate research facility in Billings, Oklahoma \citep{Stepanian+2014NAME}.  Several vertically pointing remote sensors are operated continuously at the ARM SGP site, including two polarimetric cloud radars, a wind profiler, and a Doppler lidar \citep{ARM2013}.  From 20 March through 20 May 2014, the acoustic array was deployed alongside the ARM SGP remote sensors, collecting collocated measurements (Fig. \ref{local:setup}).  We are currently in the process of analyzing the resulting dataset, but some preliminary results have been obtained that demonstrate the value of making such measurements.  


\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/local_killdeer6chan.pdf}
  \caption{Sample 1.6-sec audio spectrograms from the six recording channels of two Killdeer calls, illustrating call lags. Brighter shades denote louder acoustic amplitudes. Vertical red lines reference time of first call detection (in channels 5 and 6). } 
\label{local:kill6}
\end{figure}

The SGP site is home to a large number of Killdeer (\emph{Charadrius vociferus})---a medium-sized plover---that produce distinct and long-lasting strings of calls as they move around the site during the day (Fig. \ref{local:kill}).  While these are not nocturnal flight calls, they do serve as good validation sets for ensuring localization retrievals are producing physically realistic results, i.e., smooth, spatially continuous flight tracks.  Figure \ref{local:kill} shows a single channel for one case in which a Killdeer flew over the array while continuously calling.  The calls are faint at the start of the segment; they increase in amplitude, and then decrease again.  This pattern suggests that the bird was moving toward and then away from the array (i.e., flying over or past).  These Killdeer flyovers were quite common through the daylight portions of the experiment.  Figure \ref{local:kill6} shows a short clip of two consecutive Killdeer calls as recorded on the six microphone channels.  The birdcall wave front initially arrived at microphones 5 and 6 (times denoted by red lines) followed by the other channels. Qualitatively, this means that the Killdeer is located closest to Tower \#1, and closer to Tower \#3 than Tower \#2. Applying the localization algorithm to the six channels, each call can be localized to a source in 3D space.  Because calls occur at a sub-second frequency, their retrieval locations provide the continuous flight track of the Killdeer over the site (Fig. \ref{local:track}).  Naturally, these Killdeer are not engaging in migratory movements, but the ability to form coherent flight tracks is a good quality check for the retrieval technique.

\begin{figure}[t]
\centering
\includegraphics[width=6in]
  {figures/local_trackPrelim.pdf}
  \caption{Killdeer flight track over the ARM site. Each yellow marker indicates a call localization result.} 
\label{local:track}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Conclusions and Future Work}
\label{conc}
Radar aeroecology has steadily progressed since its inception in the 1950's, and additional advancements in radar technology, networks, modeling, and validation continue to drive the science further.  As human impacts steadily increase within the airspace and across the landscape, conflicts with wildlife will motivate new solutions for mitigating these negative interactions.  Radar will continue to play a valuable role in observing and quantifying life in the atmosphere. Radar meteorology has laid a foundation that is convenient for building biological parallels; however, blindly adopting this paradigm can often lead to limitations in our interpretations of radar observations.  As the application of weather radars in the field of aeroecology continues to progress, new challenges arise in the effort to extract additional biologically relevant information from radar data. The increasing interest in biological quantification and polarimetry is necessitating a deeper understanding of the relationships between the behavior and physical characteristics of individual organisms and the resulting radar moments. In some respects, two paths are forming toward very different goals in radar aeroecology.

The first goal is the mechanistic understanding of how collections of organisms and their behaviors create unique manifestations in polarimetric variables.  Ultimately, this seems like the ideal achievement in radar aeroecology because it would provide biological interpretations for radar measurements.  Two aspects of this process will need to be considered in addressing the challenges associated with achieving this goal.  First, interpreting polarimetric signatures of biological scatterers will not be achieved through inspecting radar data.  That is, real data has too many unknown variables to yield any rigorous, quantitative `truths'.  Our understanding of these data will only be possible through rigorous laboratory measurements (as in Ch. 4), electromagnetic models, and ensemble radar simulations (as in Ch. 7).  Most common electromagnetic models (e.g., T-matrix) are not sophisticated enough to reproduce bird-like scattering signatures, and many alternative modeling packages have never been applied to biological scatterers.  Developing adequate models will not be trivial, and even once sufficient techniques exist it will take considerable effort to model a range of taxa.  Similar challenges will be faced in measurement and simulation applications, requiring continuous, cohesive efforts over many years.  To assume that we are close to cracking biological polarimetry would be far from true.

While rigorous understanding and interpretation of biological polarimetry at the radar volume level may be decades away, that is not to say that polarimetry is not currently useful.  In fact, many topics can be addressed with even basic understanding of biological polarimetry.  Such an application in deducing flight orientations was presented in Chapter 6, and is currently being used to explore questions in migratory theory.  While precise mechanistic descriptions of polarimetric signatures might still elude us, general patterns and trends can still have great value.  With this utility in mind, the second goal of radar aeroecology is the immediate application of radar to address basic and applied problems in ecology.  Many of the motivating statements used to support efforts in radar aeroecology use appeals to conservation---many animals are experiencing dire declines from a multitude of causes, and some form of guidance is needed to inform policy and decision-making.  In this light, the purpose of aeroecology is to perform ecological triage, focusing on near-term solutions for environmental impacts and increasing knowledge of animal ecology.  This immediate application of radar aeroecology may require additional information to help compensate for our lack of understanding of biological polarimetry.  These complementary data sources might include acoustic recordings (as in Ch. 8), or citizen science reports---both of which can supplement polarimetric measurements.  Many of these near-term tasks will provide additional support to basic science in ornithology, entomology, and chiropterology. 

With these goals in mind, polarimetric radio measurements of a Brazilian free-tailed bat and Brown-headed cowbird were performed to begin to quantify the scattering characteristics of common aerial animals.  The resulting diagrams of radar cross sections versus polar look-angle were revealing; Not only does RCS vary greatly in azimuth, but the differential RCS can span over 20 decibels depending on the look-angle.  While these results were limited to X-band, it will be these types of data that ultimately lead to accurate interpretations of NEXRAD products.  Future studies will need to consider a larger range of scatterers, from insects to large birds.  Additionally, it will be critical to determine to what extent dead specimens accurately emulate the radio scattering characteristics of living subjects.  For their convenient size range and invasive status, some promising focal subjects for these studies may be House Sparrows (\emph{Passer domesticus}), European Starlings (\emph{Sturnus vulgaris}), and Rock Doves (\emph{Columba livia}).  Larger birds such as the American Crow (\emph{Corvus brachyrhynchos}), ducks, and geese will also need to be characterized.  Generally it is difficult and time consuming to make these measurements. Future efforts in radio modeling may be validated with a limited set of measurements, providing an alternative method for further electromagnetic characterization.     

Typical manifestations of biological scatterers on radar were reviewed.  The cases highlighted included diurnal echoes in the convective boundary layer, evening emergences of Brazilian free-tailed bats, morning roost exoduses of Purple Martins, and seasonal migrations of passerines.  In each case, the physical attributes and behavioral drivers of the specific organisms yield signatures that can be identified regularly in NEXRAD data.  Many nuances of these signatures can be attributed to the STAR operation of NEXRAD systems, and can vary substantially among radar sites.  Much work remains in isolating these radar-specific artifacts and distilling the radar products to their purely-biological information.  Nonetheless, polarimetric signatures in NEXRAD products (namely, $\rho_{HV}$) have been used to deduce migratory flight orientation as a function of altitude, providing an approach to analyze migratory theory in an empirical framework.  Current efforts are underway in applying this technique to address topics in avian flight in windy conditions with respect to drift versus compensation theory.  

An agent-based, time-series radar simulation suite has been developed to explore and develop techniques for the biological application of weather radars. This includes a method for emulating the behavior of specific biological agents, defining the scattering characteristics of each individual, and calculating the baseband radar return for a specified system and scanning strategy. By splitting the routines into three separate modules, each part can be varied independently to test the response to modifications to the individual modules. The emergence of a colony of Brazilian free-tailed bats was simulated and scanned by a virtual radar system, demonstrating the utility of this method for the synthesis of realistic time series signals. The simulated images were compared and contrasted with real NEXRAD data, demonstrating one of the uses of the simulation suite. While this example provides a sample of the capabilities of the simulator, additional work is needed to accurately model the behavior and scattering characteristics of bats and other organisms, thereby providing more realistic outputs. Once a robust simulator has been completed, it can be used to investigate and verify different RCS and behavioral model input parameters. That is, the input parameters can be adjusted until the simulated radar fields satisfactorily replicate actual observations. Benefits of this approach include: i) simulated movement tracks can be compared to actual foraging tracks for validation, ii) aggregations of simulated movement tracks, generated from multi-agent modeling, can be used to predict aggregate behaviors of foragers under different environmental conditions, and iii) a better understanding of bioscatter for more precise interpretation of new and existing radar observations.


Finally, a method of obtaining collocated observations of calling birds using passive acoustic localization was described.  Even without radar, the ability to connect flight calls with their location in the airspace adds value to bioacoustic recordings.  These data can provide species-specific altitudinal distributions of migrants, and their transitions within and across nights and seasons.  Measurements of altitudinal preferences during migration can be compared to meteorological conditions to deduce the decision-making processes of animals on the move.  Furthermore, our ability to provide an exact position of calling birds can allow better risk assessments.  For example, we can determine whether calling birds are flying above, below, or within the height of wind turbine rotors.  Coupled with polarimetric radar, these calls can provide organismal compositions of the airspace within radar sampling volumes, lending guidance to polarimetric interpretation, or providing a type of ground-truth.

Radar polarimetry---especially in biological applications---is still a young area of research, with much work still to be done.  As scientists and engineers continue to probe the polarimetric products, future developments in our understanding will enhance the use of polarimetric information as a tool in the study of airborne organisms. 


%\references{agsm}{myreferences}
\references{amermeteorsoc}{bibliography/master}
%comment out the next four lines if you don't have an appendix
%\appendix
%\chapter{Some Bullshit}

\end{document}