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) 480-5413 fax: (734) 480-5367 email: firstname.lastname@example.org
Award Number: N00014-05-C-0204
Oceanic internal waves, particularly large non-linear 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 term-goal 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 handed-off to an internal wave propagation model that will generate the future characteristics of the internal waves.
This program will focus on utilizing synthetic aperture radar (SAR) or real-aperture 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 end-to-end system. The two components developed above will be put together into a single, automated, system and validated at a sponsor-determine 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.
(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 two-dimensional 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 (3-6m) 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 so-called 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 time-series 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 two-layer 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 end-to-end 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 sponsor-determined location.
Documentation for the tool and the validation process will be generated
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 end-to-end 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).
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 x-direction, 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 fifth-order Runge-Kutta 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, FS, set to zero. This results in solving the so-called ray-trace equations that determine the changes in (kx, ky) 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 FS until we get back to the original spatial location (and time value). We now have the action value at the original (x,y) and (kx, ky) locations.
We chose a ray-trace Runge-Kutta 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 Runge-Kutta 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 FS. 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 non-linear sinks (such as wave breaking) of wave energy. Note that Eq. (2) is defined so that FS(Ao) = 0, which is why it is a perturbation model. Different researchers have used different values of n (the non-linear 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 ray-trace 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, Ao) 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 two-scale model that incorporates the tilting of the small-scale Bragg waves by the large-scale waves and the hydrodynamic modulations of the Bragg wave amplitudes by the surface currents induced by the orbital velocities of the large-scale 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 two-dimensional 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 C-band and X-band RCS observations of the ocean (Wackerman et al., 2002) and the standard Pierson-Moskowitz 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 Pierson-Moskowitz) and the Pierson-Moskowitz 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 right-hand 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 Pierson-Moskowitz 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 Pierson-Moskowitz 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
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/Ao-1. 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 two-scale 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 Sf is the wave spectrum which in Eq. (7) is evaluated at the Bragg wave number ks. 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 fall-off 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 re-construct a forward model simulation. We used a single internal wave, modeled with a sech2 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 re-constructed the RCS modulations.
Figure 9: Results of the full inversion system for current gradient. Nota that we essentially reconstruct the gradient after five iterations.
Ericson, E.A., D.R. Lyzenga, D.T. Walker, "Radar backscatter from stationary breaking waves," J. Geophy. Res., vol. 104, 29679-29695, 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, 832-844, 2003.
Lyzenga, D.R., J.R. Bennett," Full-spectrum modeling of synthetic aperture radar internal wave signatures," J. Geophy. Res., vol. 93, 12345-12354, 1988.
Lyzenga, D.R., "Interaction of short surface and electromagnetic waves with ocean fronts," J. Geophys. res., vol. 96, 10765-10772, 1991.
Wackerman, C.C., P. Clemente-Colon, W.G. Pichel, X. Li, "A two-scale model to predict C-band VV and HH normalized radar cross section values over the ocean." Can. J. Remote Sensing, vol. 28, 367384, 2002.
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.
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