Development of a Radar/SAR Assimilation System for Internal Wave Prediction
Dr. Christopher Wackerman 1200 Joe Hall Drive, P.O. Box 990, Ypsilanti MI 48197 phone: (734) 4805413
fax: (734) 4805367 email: chris.wackerman@gdais.com Award Number: N0001405C0204
LONGTERM GOALS
Oceanic internal waves, particularly large nonlinear ones, can have a significant impact on ship and
submarine operations when they move through a region due to the surface currents and buoyancy issues such waves
induce. Thus the Navy has a need for a predictive system that can tell a ship or submarine what the future
internal wave effects will be in their region.
The long termgoal of this project is to provide a component of such a predictive system that will take remote
sensing imagery of a region and estimate the current internal wave characteristics within that region. It is then
anticipated that these internal wave characteristics will be handedoff to an internal wave propagation model that
will generate the future characteristics of the internal waves.
OBJECTIVES
This program will focus on utilizing synthetic aperture radar (SAR) or realaperture radar (RAR) imagery to
characterize internal waves. A SAR or RAR system will only image surface effects on the ocean (such microwave
sensors do not penetrate into ocean any significant depth), so internal wave characteristic will be derived from
their surface manifestations. The most significant surface effect will be a modulation of the surface current due
to the passage of the internal wave. Thus the program will focus on estimating surface currents from radar cross
section imagery (from either SAR or RAR systems) and using those to estimate internal wave characteristics. The
program has four objectives.
(1) Inversion of radar cross section signatures to surface current gradients and internal wave
characteristics. This is the development of a tool to take the radar cross section image and generate estimates of
surface current gradients (2) Automated internal wave detection. This will develop a tool to automatically
locate internal wave signatures in radar cross section imagery and estimate their propagation characteristics from
multiple imagery.
(3) Development of the endtoend system. The two components developed above will be put together into a single,
automated, system and validated at a sponsordetermine location with in situ observations of the internal waves
for comparison.
(4) Documentation. The final objective will be to provide documentation for the final system and its validation.
APPROACH
(1) Inversion of radar cross section signatures to surface current gradients and internal wave
characteristics.
This effort will be to build the inversion scheme to estimate surface currents and model parameters from
RCS signatures. We anticipate that this will be an iterative approach, where we start with initial guesses for the
current gradients and model parameters, and then successively improve them until the simulated RCS signatures from
the forward model matches the observed RCS signature. The key to such a scheme is the estimation of the gradient
of the error (where the error is between the simulated and actual RCS observations) with respect to the parameters
being estimated (i.e. current gradient field and model parameters). We will investigate a number of ways to
generate this error gradient. One will be to assume a functional form for the surface currents induced by the
internal wave, then numerically calculate the error gradient via small perturbations of the parameters for the
surface current function as well as the forward model. This will be effective as long as the number of parameters
is small. We will also investigate removing the need for a functional form of the current, and directly estimate
the twodimensional current gradient field. This will require estimating the error gradient via the inversion of
the forward model, since the number of parameters will be too large for numerical techniques. The model inversion
will be performed via adjoint techniques, for which the models are run backwards with errors as “inputs” and thus
gradients as “outputs” to the adjoint or inverse model.
In addition, we will examine exploitation of the surface wave signatures to determine surface current maps.
High resolution (36m) satellite SAR imagery could be available that could provide some amount of dwell time over
the scene, thus allowing the estimate of wave motion in addition to wave length. This would generate both period
and wave number, which could then be inverted to derive the surface current that the wave is feeling (assuming
that the wave is not feeling the bottom). If the dwell time is not sufficient to estimate wave period, we will
also investigate the use of the socalled wave ray equations, that determine changes in wave length and direction
as a wave propagates over the current field, to invert the observed changes in the wave field into surface
currents. Note that this approach assumes that the “same” wave field exists throughout the image.
(2) Automated internal wave detection
Some work in this area has been going on in Europe under the Marine SAR Analyses and Interpretation System
(MARSAIS) program using wavelet analysis. Under our existing NOAA/NESDIS programs we have been collaborating with
the MARSAIS staff on wind and ship algorithms; we anticipate continuing this collaboration to include internal
wave algorithms. We will first determine what the state of the art is in the MARSAIS project, and then use this as
a base for our development. The biggest problem for an operational system will be the rejection of false alarms;
something that the MARSAIS project has not focused on. We anticipate that this will require spatially locating
proposed internal wave packets on a map, and then utilizing their relative relationship with each other (i.e. true
internal wave packets will be along some line with some generally known spacing) as well as their relationship to
the coastline (i.e. true internal waves will propagate in known areas and directions). This will help eliminate
errant signatures from surfactants or wind rows that may locally have similar characteristics to internal wave
signatures, but globally are not in the appropriate locations.
Under this task we will also address the multiple image problem. That is, how to determine the same
internal wave packet in two images of the same location taken at different times. We will derive metrics to
characterize the internal wave packets from each image (number of waves, length, curvature, strength, etc.) as
well as specify their spatial location. We will then perform the “matching” algorithm on the packet
characterizations to find matches that are consistent with known internal wave propagation. This will generate a
master list of timeseries properties from a given internal wave packet.
The final step in this task will be to take the master list of internal wave packet properties and derive
local oceanic environmental conditions. This part of the system will allow the user to specify sets of assumptions
(i.e. a twolayer model, constant depth of the upper layer, known bathymetry, etc.) and will then derive the
oceanic conditions consistent with the internal wave packet characterizations and these assumptions. We will
utilize the same techniques as have been discussed in the literature; thus this step will be dominantly an
implementation of known techniques.
(3) Development of the endtoend system.
The task will combine the two previously developed tools and put them together into a single, automated,
tool. This approach will be straightforward connection of the two tools. Then the system will be validated on some
sponsordetermined location.
(4) Documentation
Documentation for the tool and the validation process will be generated
WORK COMPLETED
We are currently ending the third year of this program. The system to invert radar cross section
modulations (Task 1 above) is completed and has been tested on a limited data set (see below). We anticipate that
additional testing will be performed when additional ground truth data becomes available. We are in the process of
developing the automated internal wave detection codes (Task 2 above) which will be finalized at the end of the
third year; there are no results to report yet. The development of the endtoend system (Task 3 above) will be
performed in the fourth year of the program. The original program plan had the fifth year to perform documentation
(Task 4 above).
RESULTS
The inversion tool for estimating surface currents from radar cross section modulations (RCSM) has been
completed and tested on a limited set of data. Figure 1 shows a flowchart of the tool which consists of two parts.
The top section of Figure 1 is a forward model that generates simulated RCSM from an assumed surface current
field. The bottom section of Figure 1 is the inverse model that takes as input the error between the simulated
RCSM and the actual RCSM extracted from the SAR image, then determines the change in surface current required to
decrease this error. The bottom section is done by estimating the gradient of the error with respect to the
surface current, then moving the surface current values in the conjugate gradient direction based on the
gradients. Thus the whole system is a conjugate gradient search for the surface current field that minimizes the
error between the simulated and actual RCSM. The real challenge in developing the tool is in estimating the error
gradient; i.e. the lower, right box in Figure 1.
The forward model in Figure 1 first solves the wave action balance equation to determine the spatial
changes in surface wave spectra as they propagate across the surface current field. These wave spectra are then
put into a radar cross section (RCS) model that generates the normalized RCSM to compare to the SAR imagery. The
wave action balance equation (see, for example, Lyzenga and Bennett (1988) for a full description), written in a
simplified form for clarity that assumes that the surface currents are only in one direction and only changing in
the xdirection, is where A is the action defined as the
Figure 1: Flowchart for the inversion tool.
wave spectrum divided by its frequency, the left hand side represents the transport of action in both space
and spectral domains, and FS represents all of the processes that are putting energy into the waves or taking
energy out of the waves. Action is conserved as the waves propagate through the currents, thus we use Eq. (1) to
solve for the action and then multiply by the frequency of each spectral component to generate the wave spectrum.
Initially we examined existing tools for solving Eq. (1) that used differential methods, but we eventually
developed a ray trace approach that uses a fifthorder RungeKutta method for solution. The ray trace approach
solves for the action at a give spatial location and for a given spectral location by first determining where in
the spatial domain this spectral component came from and then determining how its action changes as it propagated
to its current position. Thus Eq. (1) is first solved backwards in time with the right hand side, F_{S},
set to zero. This results in solving the socalled raytrace equations that determine the changes in (k_{x},
k_{y}) for a given wave as it moves through the spatially changing currents. This backward solution is
stopped when we are well outside of the region of changing currents. At this point we assume that the wave energy
is strictly determine by the model for the equilibrium wind wave spectrum. We then set the action using
the
equilibrium spectrum, then solve Eq. (1) forward in time with the correct sources and sinks of energy in F_{S}
until we get back to the original spatial location (and time value). We now have the action value at the original
(x,y) and (k_{x}, k_{y}) locations.
We chose a raytrace RungeKutta solver for two main reasons. First, it allowed the user to specify the
output locations (in both space and spectral domains) of the wave spectra without any connection to the grid used
to solve Eq. (1). This saved us a considerable amount to computation time, since we could specify only the
spectral locations required to generate the RCSM to compare to the SAR imagery, versus most differential solvers
that require calculating the full spectrum at each grid point needed to solve the differential equation. Second,
the RungeKutta approach allowed us to automatically determine the time spacing, and to change it dynamically, to
ensure accuracy in the solution. There was a disadvantage however; we only have a single spectral component
available as we solve the equation along the ray. That means that we have difficulty implementing equations for FS
that require the entire spectrum to be available for each spatial location.
A significant question in using Eq. (1) is the functional form for F_{S}. There is a large amount
of discussion in the literature as to the right form and parameter values. We implemented a perturbation model of
the form
where Ao is the equilibrium action, β is the wind growth term (that determines how much energy is put into
the wave system from the wind), ν is the surface viscosity (which takes energy out of the wave system), and the
last term in Eq. (2) lumps together all the nonlinear sinks (such as wave breaking) of wave energy. Note that Eq.
(2) is defined so that F_{S}(A_{o}) = 0, which is why it is a perturbation model. Different
researchers have used different values of n (the nonlinear power for the action) from 2 to 6 depending upon the
reference. In addition, there is often a spectral energy transport term added on to the expression in Eq. (2), but
that is exactly the term that a raytrace solution has difficulty implementing. There is also no real agreement on
what the wind growth term, β, should be and values in the literature can vary by a factor of 5. Finally, there is
no real agreement on what the equilibrium spectrum should be (and thus the equilibrium action, A_{o}) for
wind waves. So one of the first things we needed to do was determine the effect of the various function forms for
FS (in terms of Eq. (2) the effect of different values of n), the different models for β, and the different models
for the equilibrium wave spectra. We show this below.
The final part of the forward model is the generation of the RCSM from the wave spectra. This utilizes a
model that had been developed under previous programs which contains three scattering contributions to the RCSM.
The first is a tilted Bragg twoscale model that incorporates the tilting of the smallscale Bragg waves by the
largescale waves and the hydrodynamic modulations of the Bragg wave amplitudes by the surface currents induced by
the orbital velocities of the largescale waves. The second is a physical optics term for specular scattering from
ocean surfaces that are "glinting" to the radar. See Wackerman et al., (2002) for details on these terms. The
third term is for breaking water regions where the fraction of breaking is determined from the distribution of
downward acceleration within a ocean surface cell and the RCS for the breaking region is determine from in situ
observations that have been done on breaking water (Ericson et al., 1999; Haller and Lyzenga, 2003). All of these
terms can be generated from the wave spectrum, thus each location where a twodimensional wave spectrum has been
generated via the solution of Eq. (1) we then generate a single values for the RCSM at that location. Doing this
at a series of spatial locations we can generate a plot of RCSM values across the surface current field.
As mentioned above, there are a number of possible choices for the wind growth term in Eq. (2), a number of
choices for the value of n, and a number of choices for the equilibrium action Ao. We performed a sensitivity
study to determine how the final values of RCSM vary based on which option we used. We programmed in all eight of
the major function forms for β as well as the different values for n. Figure 2 shows a plot of the eight different
β models for a 1m wave as a function of wind speed. Note that in Figure 2 one can see the range of a factor of 5
between the values. The largest values come from a model by Plant and Wright, and the smallest comes from a model
by Kukula and Hara. We used those two β models in a forward model run where n=2 in Eq. (2) and the surface current
field was assumed to be four internal waves in a packet with peak surface currents of 0.6 m/s. We modeled the
functional form of each internal wave current field as a sech2. We used an equilibrium spectrum model that came
from Lyzenga and Bennett (1988). The result is shown in Figure 3 where we have plotted the value of radar cross
section in dB versus the normalized modulation in radar cross section. Note that the choice of the wind growth
model does not have a significant effect on the result. We then modified the value of n in Eq. (2) using the Plant
and Wright wind growth model and the Lyzenga and Bennett equilibrium spectrum. The results are shown in Figure 4.
Note that the assumption of n does have a noticeable effect on the results; larger values of n generate smaller
perturbations in RCS. Finally we looked at the effect of the equilibrium model. Figure 5 shows RCS plots for the
Lyzenga and Bennet model, a model referred to as CMOD4 that was developed to generate Cband and Xband RCS
observations of the ocean (Wackerman et al., 2002) and the standard PiersonMoskowitz spectrum. In this case we
plot the RCS modulations, since the absolute value of RCS varies significantly for the ambient region under these
different spectra. Note that here also there is a noticeable difference, although the CMOD4 (which is derived from
the PiersonMoskowitz) and the PiersonMoskowitz are similar. Thus the sensitivity study indicated that the choice
of wind growth was not significant, the choice of equilibrium spectrum was significant, but if we stayed with
"similar" spectra models the details would not matter, and the choice of n did matter with larger values of n
generating smaller RCS modulations.
Figure 2: Comparison of models for the wind growth term.
Figure 3. Results of the forward model for the smallest and largest wind growth models.
Note that there is no significant effect.
Figure 4: Results of the forward model for different functional forms of the righthand side, Fs.
Note that this is an impact with larger n generating smaller modulations.
Figure 5: Results of the forward model for different equilibrium spectrum models.
The true test is whether the forward model can predict observed RCS. To perform this test we require in
situ observations of ocean surface currents coincident with SAR imagery of the region. This was accomplished
during the Shallow Water '06 experiment off the New Jersey coast. Surface currents were measured with a towed
ADCP, a shipborne radar (the WaMOS system) and a set of ASIS buoys. To date we have only had access to the towed
ADCP data during one month of the experiment. In that data set the towed ADCP observed an internal wave packet
just before two SAR sensors collected imagery over the region. Thus we have one data set of in situ surface
currents with SAR imagery. Figure 6 shows surface currents across the internal wave packet interpolated from the
ADCP measurements (which did not measure currents all the way to the surface) via a linear interpolation scheme.
One can observe in Figure 6 currents that appear to come from a four wave internal wave packet with peaks of
around 0.6 m/s (which actually drove our simulated waves in the sensitivity study above). Figure 7 shows the RCS
modulations derived from the SAR image of the wave packet (taken from the approximate locations that the ship
moved through) compared to the forward model prediction (which is actually the PiersonMoskowitz spectrum results
from Figure 5). Note that the model gets the peaks of the RCS oscillations relatively well; the data goes from 0.2
to 0.3 and the model goes from a little less than 0.2 to a little larger. The troughs are not modeled well, but
this is due mostly to a noise floor problem in the actual SAR imagery which clips the modulations. This can easily
be incorporated into the model. These results are for the PiersonMoskowitz equilibrium spectrum, a value of 2 for
n, and the Plant & Wright beta model.
Figure 6: Surface currents interpolated from the towed ADCP during Shallow Water 06.
The four waves in the middles are across an internal wave packet.
Shallow Water 06 experiment to perform more tests and to determine the robustness of our parameter choices.
The second component of the system is the inversion model; the lower half of Figure 1. We have derived a
general expression for the gradient in the inversion
where
is the error we are trying to minimize,
Figure 7: Results of the forward model (black) using the measure surface currents in Figure 6 and
compared to the SAR image modulations (red) collected almost simultaneously. Note the
good agreement in the modulations peaks. The troughs disagree due to a thermal noise
floor in the actual SAR imagery.
and f is the perturbation in action so that f = A/Ao1. We have written Eq. (3) with respect to the current
gradient rather than the current since the RCS modulations which go into the error term in Eq. (4) are
proportional to the current gradient under many conditions.
If we put into Eq. (3) our forward model for RCS discussed above (the twoscale tilted Bragg model), we get
where s is a vector of surface tilts, G is the geometric factor for Bragg
scattering off a surface with tilts s , P( s )
is the probability of having such a tilt, and S_{f} is the wave spectrum which in Eq. (7) is evaluated at
the Bragg wave number k_{s}. The wave spectrum can be written in term of the action perturbation as
So far we have only used the RCS model in our derivation. Now we would need to use the actual wave action
balance equation to derive the expression for the derivative of the action perturbation, f, with respect to the
current gradient that is the last term in Eq. (8). The general solution requires an inversion of the full wave
action balance equation, however if we make the assumption that the relaxation rate for the Bragg waves is large
(often referred to as the lossy limit) we get a much simpler result (Lyzenga, 1991)
were the delta function is a result of the large relaxation rate assumption (in other words there is no
memory in the solution so the RCS modulations as location x only depend on other quantities at the same location),
p(k) is the falloff in k of the equilibrium spectrum, φ is the look direction of the radar and βr is the
relaxation rate which, using Eq. (2) can be written as
Substituting Eq. (9) into Eq. (8), then that into Eq. (7) generates an analytical expression for the
gradient which can be solved by simply running the RCS model once. Although it is not a general result, the
assumption that creates Eq. (9) is rather limiting, the hope was that it would be good enough to get close to the
solution. Note that we are using the assumption that generates Eq. (9) strictly for estimating the gradient, not
for the forward model.
Figure 8 shows a result of using this inverse model to iteratively reconstruct a forward model simulation.
We used a single internal wave, modeled with a sech^{2} function, that had a peak current of 0.4 m/s.
Figure 8 shows five iterations of the full system in Figure 1 where, after the fifth one, we have essentially
reproduced the RCS modulation plot. Note that this is with the simplified gradient discussed above with the full
forward model. Figure 9 shows how the current gradient changes with iteration and shows that we have essentially
estimated the current surface current gradient after the fifth iteration.
Thus the full system appears to be operating reasonably. As mentioned above we are still waiting some
additional surface current ground truth to perform additional tests.
Figure 8: Result of the full inversion system using a single internal wave as input. Note that after five
iterations we have essentially reconstructed the RCS modulations.
Figure 9: Results of the full inversion system for current gradient. Nota that we essentially reconstruct
the gradient after five iterations.
REFERENCES
Ericson, E.A., D.R. Lyzenga, D.T. Walker, "Radar backscatter from stationary breaking waves," J. Geophy.
Res., vol. 104, 2967929695, 1999.
Haller, M.C., D.R. Lyzenga, "Comparison of radar and video observations of shallow water breaking waves,"
IEEE Trans. Geosc. Remote Sens., vol. 41, 832844, 2003.
Lyzenga, D.R., J.R. Bennett," Fullspectrum modeling of synthetic aperture radar internal wave signatures,"
J. Geophy. Res., vol. 93, 1234512354, 1988.
Lyzenga, D.R., "Interaction of short surface and electromagnetic waves with ocean fronts," J. Geophys.
res., vol. 96, 1076510772, 1991.
Wackerman, C.C., P. ClementeColon, W.G. Pichel, X. Li, "A twoscale model to predict Cband VV and HH
normalized radar cross section values over the ocean." Can. J. Remote Sensing, vol. 28, 367384, 2002.
IMPACT/APPLICATIONS
If successful, the resulting system will be one component of an operational Navy tool to allow prediction
of future internal wave activity in a region so that the Navy vessels can maneuver appropriately.
RELATED PROJECTS
There are no ongoing related projects that are closely identified with this project.
Here is the original document:
Development of a Radar/SAR Assimilation
System for Internal Wave Prediction
