## Abstract

For a number of maritime tasks there is a short time period, typically only a few tens of seconds, where a critical event occurs that defines a limiting wave height for the whole operation. Examples are the recovery of fixed and rotary winged aircraft, cargo transfers, final pipe mating in fluid transfer operations, and launch/recovery of small craft. The recovery of a 30-t rescue submersible onto a mother ship in the North Atlantic Treaty Organization (NATO) Submarine Rescue System is a prime example. In such applications short-term deterministic sea wave prediction (DSWP) can play a vital role in extending the sea states under which the system can be safely deployed. DSWP also has great potential in conducting experimental sea wave research at full scale. This report explores the feasibility of using data from an experimental wave profiling radar in achieving DSWP. The report includes theory, simulation, and field testing. Two forms of DSWP are employed: a fixed point system based upon a restricted set of wave directions from which some success is obtained and the other a fully two-dimensional technique that requires further development. The main finding is that using wave profiling radar for DSWP offers promise but requires improvements both to the spatial reliability and the resolution of the wave profiling radar and to the temporal resolution of its sweep before the technique can be considered to be viable as a usable tool.

## 1. Introduction

The aim of this study is to explore the feasibility of using a wave profiling radar system called the Wave and Surface Current Monitoring System II (WaMoS II; Nieto Borge et al. 2004) as a data source for multidirectional extensions of linear deterministic sea wave prediction (DSWP) (Morris et al. 1992, 1998; Edgar et al. 2000; Belmont et al. 2003, 2006; Abusedra and Belmont 2011). The paper covers theory, potential sources of error and their mitigation, simulation testing, and results from field work.

DSWP uses measurements of the past motion of the sea’s surface to predict the actual profile of the sea surface for a short period into the future. In contrast to the mature discipline of the statistical prediction of waves (Pierson et al. 1955; Kinsman 1984; Tucker and Pitt 2001), this is a relatively new area with a very modest-sized literature (Morris et al. 1992, 1998; Prislin et al. 1997; Zhang et al. 1999; Belmont et al. 2003, 2006; Wu et al. 2000; Edgar et al. 2000; Janssen et al. 2002; Wu 2004; Blondel et al. 2008; Naaijen and Huijsmams 2008; Naaijen et al. 2009; Abusedra and Belmont 2011) that is only beginning to make the move from a research topic into practical applications.

It has been known for some time that a number of maritime operations can benefit from short-term deterministic knowledge of the sea surface shape. These range from various vessel-based applications (Belmont et al. 1995) to improvements in the performance of wave energy converters (Falnes 2002; Belmont 2009, 2010). For many of these activities there is a short time period, typically only a few tens of seconds, where a critical event occurs that defines a limiting wave height for safely conducting the whole operation. Examples of vessel applications are the recovery of fixed and rotary winged aircraft, cargo transfers, final pipe mating in fluid transfer operations, and launch/retrieval of small craft. The recovery of a 30-t rescue submersible onto a mother ship in the North Atlantic Treaty Organization (NATO) Submarine Rescue System is a prime example. In such applications, short-term DSWP can play a vital role in extending the sea states under which the system can be safely deployed (B. Ferrier et al. 2012, meeting presentation). In addition to practical applications, DSWP provides a low-cost approach to performing experimental research on sea wave dynamics at full scale.

Unlike using DSWP for research work, almost all of its practical maritime applications require real-time prediction. The maximum predict-ahead time available is set by the propagation time for waves to travel from the region where they are measured to the prediction site. Given that each set of wave measurements represents only a very partial window of the sea surface, it is necessary to treat each batch of wave data used to build a prediction model as independent. This means that all computations needed to create a prediction model must be completed in times much shorter than the prediction time because all such calculations deduct directly from that prediction time. Moving vessels operating independently must typically use some form of remote sensing for wave measurements, such as the experimental WaMoS II (OceanWaves GMBH 2013) wave profiling radar system, which is based upon work by Nieto Borge (1998), Nieto Borge et al. (1999, 2004), Alpers and Hasselmann (1982), Ziemer and Günther (1994), Plant and Zurk (1997), and Hessner et al. (2002), or shallow-angle wave profiling lidar (Belmont et al. 2007). Given that the measurement horizons for these techniques are of kilometer scale, this produces maximum predict-ahead times of the order of 1 min, meaning that all data processing and prediction model building must be completed in a few seconds. Thus, while work has been undertaken on nonlinear DSWP (Prislin et al. 1997; Zhang et al. 1999; Wu et al. 2000; Wu 2004; Blondel et al. 2008), the need for short computational time scales means that generating nonlinear multidirection prediction models is unrealistic without supercomputer resources. So, for most practical maritime operations, real-time DSWP is restricted to using linear sea models (Morris et al. 1992, 1998; Edgar et al. 2000; Belmont et al. 2003, 2006; Abusedra and Belmont 2011).

This limitation to linear models may seem very restrictive because all wave systems are at least weakly nonlinear. However, in most practical vessel-based operations the use of DSWP is in so-called quiescent period prediction (QPP) (B. Ferrier et al. 2012, meeting presentation) QPP exploits the well-known phenomenon that, under most sea conditions, intervals of large waves alternate with smaller ones; the aim is then to predict the occurrence of the more quiescent intervals. In such cases it is only necessary to determine that the wave height is less than some value; the precise amplitude and profile are not significant and thus predicting the detailed wave shape is not critical. Only in highly cluttered, fully two-dimensional, highly nonlinear seas driven by local wind waves would this not be sufficient. For such seas to be large enough to inhibit most maritime tasks of interest, the prevailing wind conditions would be so severe as to cause the operations to be suspended. Such wave-limited, vessel-based applications are normally constrained by large swells created by a very limited number of remote large storm systems, with modest directional spreading, and are generally accepted to be well described by linear wave models (Kinsman 1984; Tucker and Pitt 2001).

Clearly an affordable, convenient remote sensing system is a key requirement in DSWP. Satellite techniques at present cannot provide either the required spatial resolution or the coverage. Shallow angle wave profiling lidar (Belmont et al. 2007) is very much a research technique, and wave buoys are of no value in moving vessel applications. This leaves wave radar as the remaining candidate. Traditionally, wave radars return wave statistics; however, a commercial product called WaMoS II (Nieto Borge et al. 2004) is being further developed to provide deterministic wave data over a two-dimensional region using standard vessel navigation radars. The aim of the present study is to explore the feasibility of multidirectional extensions of linear DSWP (Morris et al. 1992, 1998; Edgar et al. 2000; Belmont et al. 2003, 2006; Abusedra and Belmont 2011) using this system as a data source.

This paper describes two approaches to linear DSWP using wave profiling radar. The first extends the so-called *fixed point* method introduced for one-dimensional seas by Morris et al. (1992, 1998) to multidirectional seas, using wave height time series data obtained at a modest number of fixed locations that are equivalent to the outputs from an array of heave-only wave sensor buoys. The second method involves a fully two-dimensional sea model that uses all of the available wave data obtained from an area scan of the sea surface.

In section 2 we introduce the theory behind fixed point linear DSWP. There are a number of potential sources of error in this model; these are described in section 3, along with methods for mitigating against them, or at least measuring their effect. In section 4 we describe the results from simulations using synthetic data and in section 5 real results from a sea trial. The fully two-dimensional DSWP model is described in section 6, with some results from both synthetic and real data. Section 7 concludes.

## 2. Theory: MFP linear DSWP

As stated the multiple fixed point (MFP) method of linear DSWP uses wave height time series data obtained at a modest number of fixed locations that are equivalent to the outputs from an array of heave-only wave sensor buoys. Each wave time history (one from each sensor location) makes it possible to model a separate, polychromatic wave system propagating in a different direction. This approach assumes that the wave system is created by a modest number of remote storm systems, each with a relatively small angular spread in wave directions. Given that the wave radar system makes data available ideally as a sequence of “snapshots” over a finely spatially sampled area, this means that only a small fraction of the available wave data is used in this technique.

### MFP DSWP for multidirectional seas

We employ the standard linear oceanographic wave model (Kinsman 1984; Tucker and Pitt 2001): the wave height, *h*(*x*, *y*, *t*), at the spatial coordinates *x*, *y*, and temporal coordinate *t* is given by

where *N* is the number of frequencies employed, *R* is the number of significant storm directions, *A*(*ω*_{n}, *θ*_{r}) is the directional magnitude spectrum, *θ*_{r} is the propagation direction of an individual wave component, *ω*_{n} is the angular frequency, *k*_{n} is its wave number, and Φ(*ω*_{n}, *θ*_{r}) is the phase. In this case *R* is modest, typically *R* < 10. For convenience Eq. (1) is recast into complex form as

where *C*(*ω*_{n}, *θ*_{r}) are the complex Fourier coefficients. The sgn(*n*) function, which returns the sign of *n*, is required because the deep-water form of the dispersion relationship typically used means that the real and imaginary parts of the right-hand side of Eq. (2) would otherwise not satisfy the necessary symmetry relationships needed to ensure that *h*(*x*, *y*, *t*) is a real function.

Using sea surfaces obtained from overlapping wave radar scans mean that at each measurement location, designated by *x*_{i}, *y*_{i}, a time series of wave height values, *h*(*x*_{i}, *y*_{i}, *t*_{l}), where 1 ≤ *i* ≤ *R* and 0 ≤ *l* ≤ *N* − 1, is available. This allows a discrete Fourier transform (DFT) of Eq. (2) to be taken at each measurement location, giving a system of *R* complex frequency domain equations of the form

where *H*(*x*_{i}, *y*_{i}, *m*) is the DFT of the set *h*(*x*_{i}, *y*_{i}, *t*_{l}). This system is conveniently written in matrix form as

where is the *R* × 1 spectral data vector, [Coeff] is an *R* × *R* coefficient matrix whose elements are the terms exp {*j*[sgn(*m*)*k*_{m} cos(*θ*_{r})*x*_{i} + sgn(*m*)*k*_{m} sin(*θ*_{r})*y*_{i}]}, and is the *R* × 1 vector of unknown Fourier coefficients. Inverting the matrix Eq. (4), either explicitly or via a least squares fitting process, yields estimates of the vector of Fourier coefficients that when substituted into Eq. (1) produce a prediction model for the wave height *h*(*x*, *y*, *t*) at the required spatial and temporal prediction coordinates of the prediction site.

This assumes that the wave radar only delivers wave elevation at each point; hence, *R* data locations are equivalent to *R* heave-only wave buoys. If, however, the north–south, east–west motions can be extracted from the radar (e.g., via wave-induced surface velocities), then any two of the three equivalent buoy motions can be used rather than merely heave. This allows *R* wave measurement locations to model 2*R* wave directions. The *R* new equations are simply the appropriate projections of *π*/2 phase-shifted versions of the heave equations (one-quarter of a particle rotation later than heave) of the form

The reason for 2*R* directions and not 3*R* is because any two motions of the three (heave, north–south, east–west) are sufficient to define the circular water particle motion for each wave component; hence, only two of these datasets are linearly independent.

#### Wave directions

The above-mentioned process requires estimates of the set of wave directions, *θ*_{r}, where 1 ≤ *r* ≤ *R*. Given the physical assumption that the sea of interest is created by a modest number of remote storm systems, it is reasonable to assume that these directions remain unchanged over periods of at least several minutes. This is in contrast to the *C*(*ω*_{m}, *θ*_{r}) Fourier coefficients which, for reasons given earlier, must be updated with every new wave dataset. Consequently, standard statistical directional spectral estimation methods can be employed to estimate the *θ*_{r}. In real-time DSWP, this statistical estimation is run as a background task, updated on a time scale of approximately 10 min. A particularly convenient approach is an extension to multiple data gathering sites of the original Longuet-Higgins et al. (1963) cross-spectral method, using Borgman (1979). Because of problems such as the generation of nonphysical negative power density spectra, caused by Gibbs phenomena, such methods are normally superseded by various forms of maximum likelihood estimators. However, as only the locations of the peaks in the directional spectra are needed for DSWP, the nonphysicality is not critical. Alternatively, for the reasons given below, it is computationally realistic to employ a fast optimization process to estimate the wave directions:

The accuracy of the predictions can be measured at the prediction site “after the event.”

The term

*R*is small,*R*≤ 10.The

*θ*_{r}changes slowly and hence angle estimation is a background task repeated at the fastest rate of every 10 min. Furthermore, rapid convergence will be helped by using very good initial guesses based upon the previous values.

## 3. Potential sources of error and their mitigation

Linear DSWP can only approximately predict the sea surface profile, and clearly it is necessary to assess the level of confidence that can be placed in predictions and where possible how to mitigate the errors, or at least measure their effects.

### Potential errors

The major factors affecting errors in linear DSWP are as follows: (i) the sea is only ever approximately linear; (ii) the dataset used for estimation is finite; (iii) the sea is not periodic, whereas the Fourier series sea wave model constituting Eq. (1) is strictly periodic; and (iv) the accuracy of estimating the Fourier coefficients in the matrix Eq. (4) is affected by the extent to which the system is mathematically well conditioned.

#### Addressing errors

For multidirectional seas a particularly important nonlinear effect is the degree of energy transfer between modes (during wave propagation after measurements) stemming from directional wave–wave resonance (Phillips 1960; Longuet-Higgins and Phillips 1962). The inherent restriction of practical DSWP to wave propagation distances of the order of 1 km naturally limits the effects of this source of error.

The departure from linearity can be measured by using standard bispectral methods (Kim and Powers 1979). This does not remove the effects of nonlinearity but does provide a metric for assessing the confidence that can be placed in predictions.

The effect of using a finite dataset involves the well-known phenomenon of “windowing leakage errors” in the frequency domain. In signal processing applications, these errors are often mitigated by using smooth window envelopes such as Hanning functions. Unfortunately, these introduce global distortion. A more suitable technique is so-called “end matching” (Morris et al. 1992, 1998; Abusedra and Belmont 2011), which searches for the most periodic subsets within the wave data used. This approach not only minimizes windowing errors (Brigham 1988) but also addresses the inappropriate periodicity problem, which is actually closely linked.

Of all the sources of error, the most significant is concerned with the degree of mathematical ill conditioning of Eq. (4). Given an uncertainty, , in the vector (e.g., experimental measurement errors), it is well known that any form of inversion produces an uncertainty, , in , which satisfies the inequality

where cond([Coeff]) is the condition number of the coefficient matrix [Coeff] and the denotes any suitable norm. Clearly for cond([Coeff]) > 1, any errors associated with are inflated.

In general the condition number of a system rises with its rank and unfortunately [Coeff] is closely related to the class of Vandemonde matrices (Golub and Van Loan 1996) that are well known to be especially poorly conditioned. The strong effect of rank in this case is why it is necessary to restrict *R* to modest values.

Now the elements exp{*j*[sgn(*m*)*k*_{m} cos(*θ*_{r})*x*_{i} + sgn(*m*)*k*_{m} sin(*θ*_{r})*y*_{i}]} of [Coeff] depend upon the wave directions and the locations of the measurement points relative to the prediction site. Thus, for a given set of sea conditions, there will be a best set of measurement locations with respect to system conditioning and as with the wave directions these can be estimated via an optimization routine run as a background processing task.

## 4. Simulation testing

Comprehensive simulation testing of DSWP implementation algorithms would require exploration of the directional characteristics, number and location of the measurement locations, data window lengths, prediction bandwidth, etc. Even presenting the results of such an undertaking is totally unrealistic and thus only individual illustrations can be shown.

### a. Tests solely on the prediction process

Figure 1 shows results designed to test only the prediction technique, and hence the mean wave directions for each wave system were taken as known quantities. The goal was to assess the goodness of fit of the predicted heave values at the prediction site estimated over an interval up to 30 s into the future. The results were obtained under “blind trial” conditions in which a third party generated the wave data from the sea model generator described.

In total 82 different sea scenarios were tested, each of which was provided with 3600 observations, sampled at 1 Hz. For each scenario, 1000 s of data were used to predict at 1 Hz up to 30 s into the future. A sliding window approach was used, moving along the dataset by one observation each time, resulting in 2571 sets of 30-s predictions for each scenario. These predictions were measured against the actual values using Pearson’s linear correlation (Pearson 1920). The boxes in Fig. 1 represent the 25th–75th percentile of the resulting 2571 correlation values for each scenario; the center of the box marks the median and the whiskers the full range of values.

The sea models used were produced over a range of sea states from various standard forms (Tucker and Pitt 2001). To avoid any individual wave systems being too long crested their directional magnitude spectra were very finely resolved, with 360 directions modeled for each storm system over a typical angular spread of 20°. The phases for each separate direction in the spectra were uniformly randomly distributed over the range −*π* to +*π*.

Three measurement points were used in the prediction process, specified over the half plane known to contain the wave directions at distances greater than 200 m from the prediction site. This allowed only three wave directions to be modeled in contrast to the 360 used in the sea model.

Given the coarse nature of the directional modeling by the DSWP as compared to the finely resolved sea model and the deliberately nonoptimum distribution of measurement locations, the overall performance of the prediction process is deemed satisfactory, especially recalling the requirements of quiescent period prediction. What is clear is that the distribution of goodness of fit is bimodal; the correlations for the various sea conditions are either good or poor. This reinforces the point that in practical DSWP applications, it is important to provide feedback from the actual predicted quantity of interest, typically wave-induced vessel motion, for under certain sea conditions the DSWP process fails.

#### Overlapping estimates (forward marching prediction)

Given a continuous stream of wave data from each measurement site, it is possible to update the DSWP prediction model on a regular basis, up to a maximum rate of once each time step. This produces a common forward marching prediction window within which multiple estimates of the prediction at a given future instant are available. The possibility of such multiple estimates allows for considerable increases in the confidence level in the predictions. An illustration of this for a well-predicted scenario is shown in Fig. 2 in which the common interval is 30 samples, and hence 30 separate estimates are available at each predicted instant. In the figure the black line is the actual heave, each of the predictions is represented by a gray line, and the forward marching process moves forward a total of 200 samples.

### b. Condition number effects

Taken together, cond([Coeff]) and are influenced by all the relevant physical parameters defining both the sea conditions and the DSWP process. Thus, again an exhaustive exploration of condition number behavior over the whole parameter space is unrealistic. Furthermore, inequality [Eq. (6)] provides an upper bound and not an actual value for error. What can be said anecdotally, based upon very partial parameter explorations, is that under the general conditions of the simulations, large condition numbers—for example, cond([Coeff]) > 100—ensure that predictions were very poor, while conditions producing substantially lower values did not necessarily guarantee optimum prediction results.

As with the wave directions, the best set of measurement locations can be determined by fast optimization run as a background task, again using “after the event” measured prediction errors as the cost function. To illustrate this process, a search based on a genetic algorithm (GA) (Goldberg 1989) was undertaken to identify the best measurement positions. In the simulation the wave system was set up to be centered upon a north-northwesterly direction, and hence the search was restricted to the region shown. Given the wave nature of the problem, it was expected that there would be many near-optimal solutions and this was found in practice. Five almost equivalent clusters of locations are illustrated in Fig. 3.

### c. Measures of nonlinearity

As discussed in section 1, practical constraints mean that only linear DSWP can be used for real-time DSWP; thus, there are almost certain to be situations when the linear DWSP predictions will be very poor. It is thus necessary to provide a measure of confidence in the prediction accuracy and hence indicate when the predictions can be safely used in advising operational decisions.

The natural approach to this is to simply compare the predicted wave-induced vessel motion (using DSWP as input to a vessel motion model) with those measured and assume that the results are representative. The problem with solely relying on this approach is that it is an “after the fact” test and that there may be conditions when one good prediction does not absolutely guarantee that the next prediction will be equally good. These circumstances might occur when there is a significant level of nonlinearity present, where in addition to the obvious direct nonlinear effects even the zeroth-order linear aspects of the wave system can become highly nonstationary due to phenomena such as directional wave–wave resonance (Phillips 1960; Longuet-Higgins and Phillips 1962). This strongly influences the estimates of wave direction that are assumed to remain unchanged for significant periods of time.

Thus, in addition to the “after the fact” measures, what is also required is a metric that measures the departure of the present sea wave conditions from the linear behavior assumed by linear DSWP. In the present context linearity is equivalent to the wave system statistics being Gaussian and nonlinearity induces higher-order statistics. A well-established technique that tests equivalently for departure from Gaussian and for nonlinearity is the bispectral method (Kim and Powers 1979). This is based upon a third-order correlation function and allows the derivation of a measure that can be shown to be zero for linear/Gaussian processes and increases with departure from these conditions. Figure 4 shows a frequency domain bicoherence plot (Kim and Powers 1979) for a sequence of real sea wave data as estimated by the WaMoS II wave radar system. The grayscale starts at white for linear/Gaussian with increasing lightness indicating the presence of increasing levels of nonlinearity. The two axes are frequencies, and the location of the lighter pixels indicate the spectral interaction regions. In effect such a plot shows the presence of frequency domain mixing. Such two-dimensional plots are very data intensive, so in order to extract a simple single-value metric for nonlinearity, the magnitude of the plots are averaged over the frequency domain.

Given such an integrated bispectral metric (IBSM), if the time evolution of this number correlates significantly with the prediction quality, it is very probable that nonlinearity is significantly affecting the prediction process and hence linear DSWP should not be relied upon.

## 5. Sea trial

A sea trial was conducted as part of the NATO Submarine Rescue System (NSRS) program in order to begin the process of assessing the suitability of both a DSWP system and of deterministic wave profiling radar (as a data input to the DSWP process). The eventual goal was to produce a quiescent period prediction system aimed at extending the sea-state operational envelope of NSRS. The wave profiling equipment used was the WaMoS II, which is signal processing/software technology, undergoing commercial development, that employs data from standard marine radars (Dittmer 1995) to estimate deterministic wave profiles as well as standard sea wave statistics. A fundamental challenge inherent in the trial was that very limited independent validation data were available for the wave radar, and hence considerable care was needed in interpreting the outcomes of the trial.

The standard operating conditions for the NSRS deployments are that the mother ship, typically a vessel of approximately 3000–8000 tonnes equipped to launch and recover the 30-t rescue submersible, steams at extremely slow speed in the direction of the prevailing sea wave system. This allowed multiple scans of the wave profiling radar, ahead of the vessel, to be overlaid producing a 1.2 km × 1.2 km region common to a large number of the scans. Within this common region, fixed spatial points could be selected at which wave height time series could be obtained to serve as inputs for the MFP form of DSWP. The wave height data at an additional point, “down wave,” within the common region was intended to check the accuracy of the predictions. The standard statistical processing mode of WaMoS II also provided one possible source of the directional wave spectral data required for assigning the wave directions needed by the prediction model. In this feasibility study, the DSWP system was not intended to operate in real time, with all calculations being performed ashore after completion of the trial. This made it possible to corroborate the WaMoS II wave directions with those produced by the hind casting MetOcean wave model operated by the Met Office (Stretch 2012).

### Experimental process

#### 1) Radar estimation of sea surface elevation

The following provides a brief summary of the wave profile estimation process. Interested readers are referred to the cited literature.

It is known that under various conditions, signatures of the sea surface are visible in the near range (<3 n mi) of marine radar images. These signatures are known as sea clutter because they are undesirable for navigation purposes and generally suppressed by filter algorithms; the longer waves become visible in the radar images because they modulate the sea clutter signals (Hessner et al. 1999; Nieto Borge 1997; Seeman et al. 1997). This modulation is a nonlinear process and so the sea clutter radar image intensities do not have a one-to-one mapping with the sea surface elevation.

The standard WaMoS II wave analysis is derived from a Cartesian subset of the radar images, typically 0.25–2.25 km^{2}. The digitized electromagnetic (EM) intensity, that is, “sea clutter,” Ξ(*x*, *y*, *t*), is Fourier transformed in space and time to yield the three-dimensional image spectrum *χ*(*k*_{x}, *k*_{y}, *ω*):

where F(·) is the Fourier transform operator, **k** = (*k*_{x}, *k*_{y})^{T} is the wave vector, |**k**| = *k* = 2*π*/*λ* is the wavenumber, with *λ* as the wavelength, and *ω* = 2*πT* is the angular frequency with the period *T*. The energy related to the surface waves is localized in the image spectrum following the dispersion relation for linear gravity waves:

where *g* is the acceleration of gravity, *H* is the local water depth, and **U** = (*U*_{x}, *U*_{y})^{T} is the surface current. By determining the current and filtering *χ*(*k*_{x}, *k*_{y}, *ω*), the wave image spectrum is obtained:

As the digitized EM intensity is not linearly related to sea surface elevation, a modulation transfer function, MTF(·), (Plant 1989; Ziemer and Günther 1994; Nieto Borge et al. 2004) is required to derive the wave spectrum. By integrating over all positive frequencies the directional wavenumber spectrum is obtained:

The directional wavenumber spectrum is then inverse Fourier transformed to yield the sea surface elevation:

#### 2) Experimental data

Figure 5 shows the arrangement of the vessel and the common measured region, approximately 1.2 km × 1.2 km, used for the experiment. The measurement locations used were jointly optimized for the condition number of the coefficient matrix and for prediction quality. The three black spots mark the chosen measurement locations, while the small black square is the prediction site. For these tests, 200 observations were used to predict 20 observations ahead, with the observations spaced at 2.52-s intervals.

#### 3) Results

The black line in Fig. 6 is a plot of the time evolution of the prediction correlation coefficient at the prediction site. It clearly shows the presence of quasi-periodic variations from reasonable prediction quality (certainly good enough for QPP) to antiphase predictions. Figure 7 provides an illustration of wave profile fits for a short period around test 30.

#### 4) Consideration of errors

Given that the measurement locations had been optimized as described above, it is unlikely that condition number was the cause of the cyclic errors. To determine the effect of nonlinearity, the time evolution of the IBSM measure was plotted against prediction quality, as illustrated in Fig. 6; the gray lines are the bicoherence measure at the prediction site (solid line) and at each of the measurement sites (dashed and dotted lines). The bicoherence was measured over the 200 observation training interval, which was split up into four segments of 50 observations each in order to perform the calculation. No evident relationship was found, so nonlinearity was ruled out as a case of the quasi-periodic error. There are no reasons to expect that windowing or discretization effects would induced the observed long-period quasi periodicity. The only remaining source was some form of time-dependent resolution error in the radar-estimated wave profiles that was further supported by anecdotal reports of the need to locally spatially realign wave height estimates in previous wave radar trials.

## 6. Two-dimensional linear DSWP (TWD)

As stated in section 4, the general findings from the field work using WaMoS II data as input to MFP indicate that there are almost certainly some spatial referencing issues with the wave profile estimation that require resolving. Furthermore, MFP inherently involves optimizing the measurement locations used. Thus, providing the remote storm-generated wave systems have reasonably narrow angular spreads, this approach is likely to lead to best-case results. In contrast the fully two-dimensional approach to linear DSWP requires the use of all the available data across the measurement region and in no sense is it optimized with respect to data locations. Thus, it would be expected that under the present circumstances, the two-dimensional prediction might perform significantly worse than MFP. This was found to be the case for the field data and so the work on the two-dimensional reported here has been restricted to merely demonstrating the basic approach; no efforts were made to refine the technique (as had been undertaken in the MFP case).

### a. The two-dimensional prediction process

Given an appropriately sampled spatial array of data, it is possible to employ a two-dimensional Fourier transform approach to linear DSWP. The starting point is a continuum version of Eq. (2), that is,

Again, the sgn(*ω*) functions return the sign of *ω* and are included to force the appropriate symmetries that guarantee that *h*(*x*, *y*, *t*) remains a real function. By employing a “Fourier like” orthogonal integral transform, Eq. (12) can be inverted to produce

Substituting for *C*(*ω*, *θ*) into Eq. (12) provides a two-dimensional frequency-domain-based prediction model. For sampled data the above equation discretizes in an obvious manner to

where *P*.*δx* and *Q*.*δy* satisfy both the spatial form of the Nyquist sampling theorem and the domain size requirements for the lowest value of wavelength employed. Thus, a direct approach to the two-dimensional linear DSWP involves the following steps:

Acquire the sampled spatial data,

*h*(*p.δx*,*q.δy*,*t*_{0}), over the region 0 ≤*x*≤*P.δx*, 0 ≤*y*≤*Q.δy*at a given time*t*_{0}, and transform the spatial data according to Eq. (13), which in practice can be recast into a standard base two fast Fourier transform.Substitute the resulting estimate for the complex vector

*C*(*ω*_{n},*m.δθ*) into Eq. (2), which can be used to predict the time wave elevation at the desired location and time.

### b. Space domain wave filters

By substituting for *C*(*ω*, *θ*) from Eq. (13) into Eq. (12), it can be shown that the prediction problem can be recast into a spatial convolution:

where ϒ(*x*, *y*, *t*, *g*) is the spatial impulse response function:

Alternatively, the problem can be recast into the wavenumber domain, producing this time a spatial impulse response function of the form

This is an infinite impulse response wave-filter approach to prediction as opposed to the frequency domain phase-shifting approach employed so far. The filter technique offers some advantages and some disadvantages. The advantages are the avoidance of explicit periodicity and that asymptotic analytic forms can be developed for the impulse responses to reduce computational costs. However, two new problems raised by the methodology are (i) the impulse responses are noncausal and (ii) they have unbounded derivatives in the limit of large *x*, *y*, *t*. Resolving these two issues requires conjugate domain truncation, which is explored in a one-dimensional treatment of wave filtering (Belmont et al. 2006).

### c. Simulation testing of two-dimensional DSWP

Simulation of the two-dimensional frequency domain prediction method is illustrated for the case of a single-frequency sinusoid with a 10-s period propagating at an angle of *π*/4 rad. The results are shown in Fig. 8. The modest errors are primarily a consequence of a combination of discretization and finite data windowing effects.

### d. Sea trial using two-dimensional DSWP

Using the prediction code employed in the simulation shown in section 6c, the sea model data were replaced with WaMoS II wave profile estimates from the sea trial to estimate wave height time series at a single point in the measured region. The same dataset was used as in the MFP method, although clearly unlike in MFP, all the spatial data (for each radar scan) were employed in the two-dimensional predictions. The results, presented in Fig. 9, are very poor compared to those in Fig. 6, produced using MFP. Figure 10 illustrates wave profile fits for the best case obtained. However, for the seas observed the ratio of spatial resolution to wavelength was much larger than expected under NSRS operating conditions when significantly better performance was expected.

## 7. Conclusions concerning the sea trials results

The general finding from the sea trials is that using the multiple fixed point (MFP) method, there are significant intervals where the prediction is good enough for many practical marine tasks, such as quiescent period prediction for the NSRS application. However, these alternate, on a reasonable predictable quasi-periodic basis, with intervals of increasingly poor prediction.

We also find that even under the observed conditions, which are at the threshold of viability for the wave profiling radar system, for substantial time periods the wave prediction technology is sufficiently successful to already be of practical value in certain quiescent period prediction applications. The sponsor of this work, the U.K. Ministry of Defence, was sufficiently encouraged by the results presented here to support further work on this technology with the practical goal of developing an operational sea-going system.

Consideration of the various sources of error suggests that there are problems with the accuracy of spatial locations of the wave height values estimated by the wave radar; resolving this will clearly involve further development of this technology. Despite this drawback the relatively smooth time dependence of the prediction quality indicates that even at the present state of development, a combination of the MFP method of DSWP using WaMoS II wave data can provide valuable additional guidance in wave-limited marine operations.

## Acknowledgments

The authors acknowledge funding from the European Union FP7 and U.K. Ministry of Defence NSRS projects.

## REFERENCES

*Offshore Focus,*June 1995 edition.

*J. Mar. Eng. Technol.,*

**2003,**21–26.

*Proceedings of the 23rd International Workshop on Water Waves and Floating Bodies,*H. S. Choi and Y. Kim, Eds., IWWWFB, 13–16. [Available online at http://www.iwwwfb.org/Abstracts/iwwwfb23/iwwwfb23_04.pdf.]

*Ocean Wave Climate,*M. D. Earle and A. Malahoff, Eds., Marine Science Series, Vol. 8, Plenum Press, 269–300.

*The Fast Fourier Transform and Its Applications.*Prentice Hall, 448 pp.

*Ocean Waves and Oscillating Systems: Linear Interactions Including Wave-Energy Extraction.*Cambridge University Press, 286 pp.

*Genetic Algorithms in Search, Optimization and Machine Learning.*Addison-Wesley Longman Publishing Co., 432 pp.

*Matrix Computations.*3rd ed. Johns Hopkins University Press, 728 pp.

*IGARSS ’99 Proceedings,*Vol. 1, IEEE, 500–502, doi:.

*Ocean Wave Measurement and Analysis,*B. L. Edge and J. M. Hemsley, Eds, ASCE, 221–230, doi:.

*Ocean Wave Measurement and Analysis,*B. L. Edge and J. M. Hemsley, Eds, ASCE, 377–387, doi:.

*Wind Waves: Their Generation and Propagation on the Ocean Surface.*Dover Publications, 676 pp.

*Ocean Wave Spectra,*Prentice-Hall, 111–136.

*Manoeuvring and Control of Marine Craft: Proceedings of the Second International Conference,*P. A. Wilson, Ed., WIT Press,

*Ocean Engineering: Offshore Renewable Energy,*Vol. 4,

*Proceedings of the ASME 2008 27th International Conference on Offshore Mechanics and Arctic Engineering,*ASME, OMAE2008-57804, 607–614, doi:.

*Ocean Engineering; Ocean Renewable Energy; Ocean Space Utilization, Parts A and B,*Vol. 4,

*Proceedings of the ASME 2009 28th International Conference on Offshore and Arctic Engineering,*ASME, OMA2009-79366, 243–255, doi:.

*Radar Scattering from Modulated Wind Waves: Proceedings of the Workshop on Modulation of Short Wind Waves in the Gravity-Capillary Range by Non-Uniform Currents,*G. Komen and W. Oost, Eds., Kluwer Academic Publishers, 155–172.

*Oceans ’97 MTS/IEEE: Conference Proceedings,*Vol. 2, MTS/IEEE,

*Waves in Ocean Engineering.*Elsevier Ocean Engineering Series, Vol. 5, Elsevier, 548 pp.

*Proceedings of the 15th International Workshop on Water Waves and Floating Bodies,*IWWWFB, 191–194.

*Second Int. Conf. on Air–Sea Interaction and on Meteorology and Oceanography of the Coastal Zone,*Lisbon, Portugal, Amer. Meteor. Soc., 117–118.