Imaging 3 D Sea Surfaces from 3 D Dual-Sensor Towed Streamer Data

3D realistic sea surface imaging from 3D dual-sensor towed streamer data is presented. The technique is based on separating data acquired by collocated dual-sensors into up-going and down-going wavefields. Subsequently, these wavefields are extrapolated upwards in order to image the sea surface. This approach has previously been demonstrated using 2D data examples. Here, the focus is on 3D data. Controlled 3D data based on the Kirchhoff-Helmholtz algorithm is generated, and the 3D sea surface imaging technique is applied. For coarsely spaced streamers from 3D field data, the technique is applied streamerwise (i.e., 2D wavefield separation, extrapolation, and imaging). In the latter case, the resulting sea surface profiles corresponding to each time frame are interpolated to demonstrate that themain sea surface characteristics are preserved, and artefacts due to 2Dprocessing of 3D data are mainly limited to areas corresponding to large angles of incidence. Time-varying sea surfaces from two different 3D field data are imaged.The data examples were acquired under different weather conditions. The imaged sea surfaces show realistic wave heights, and their spectra suggest plausible speeds and directions.


Introduction
In marine seismic acquisition, rough sea surfaces cause amplitude and phase perturbations in the acquired seismic data [1].However, the corrections applied to the data during processing are often inadequate.In addition to the timevarying nature of the sea surface, seismic streamer depth may also vary with time.This is more visible in time-lapse seismic imaging [1] where successive seismic images of producing field aid geophysicists in identifying bypassed oil, but the effects of time-varying sea surface prove to be a challenge when these images are matched.This is because the temporal distortions introduced in the data by a time-varying rough sea (and fluctuating streamer depth) are neglected.
Remote sensing can provide detailed (i.e., smaller wavelengths of the sea surface are imaged) digital elevation models of a sea surface using across-track interferometric synthetic aperture radar (InSAR) [2] (e.g., Schulz-Stellenfleth et al. 2001).However, these elevation models are distorted (depending on the amplitude of the ocean swell) because of the continuous motion of the sea wave.In addition, timevarying sea surface information is not available from satellites because of data size.Moreover, the acquired data may not be ideal for correcting sea surface distortion of marine seismic data.On the other hand, dual-sensor technology can provide time-varying sea surfaces which are smoothed because of a relative low frequency band of the seismic reflection response.
Laws and Kragh [3] propose that specialized deconvolution methods can be used to remove such distortions if timevarying sea surface elevation measurements are available.They derive sea surface elevations and vertical streamer displacements from time-varying pressure measurements using a specialized hydrophone setup.However, this method is dependent on the precision of the pressure sensors at very low frequencies (between 0.05 Hz and 0.3 Hz), which is well below normal seismic bandwidth.A method of imaging 2D sea surfaces from 2D seismic data acquired by collocated dual-sensor systems (comprising pressure and particle velocity sensors) has been described by Orji et al. [4][5][6].This paper investigates the robustness of the proposed imaging International Journal of Geophysics technique to reconstruct 3D time-varying sea surfaces from 3D towed dual-sensor streamer data.
Many authors discuss scattering of wavefields from rough surfaces based on the Kirchhoff methods (i.e., the Kirchhoff-Helmholtz integral where the Kirchhoff approximation is used in solving the normal derivative of the wavefield at the free surface).Among these authors we cite [7][8][9][10][11] among others.Given the frequency band of typical marine seismic sources and the degree of sea surface roughness, the Kirchhoff method seems to be a valid approach when computing scattered seismic data for marine seismic cases.Orji et al., [5,6] presented a forward modeling approach based on the 2D time-domain Kirchhoff-Helmholtz integral, which can handle 2D time-varying rough sea surfaces, and used it to validate the sea surface imaging technique.Here, 3D extension of the forward modeling tool is used in generating controlled data.The effort in obtaining 3D imaged sea surfaces is based on the fact that most field data nowadays are 3D.Given the constraints in 3D data acquisition and processing (e.g., crossline sampling), the imaging technique needs to be verified by 3D in order to demonstrate various practical challenges encountered in imaging 3D sea surfaces from coarsely spaced streamers.For spectral analysis of the imaged sea surface, data acquired from a larger area may give better information about the sea surface.Thus, in this work we compare 2D and 3D processings and also carried out spectral analysis where we show that limited sea surface area poses a problem.
Up-going and down-going pressure wavefields are computed for 3D realistic sea surface conditions.The sea surface conditions are computed based on wind wave models such as the Pierson-Moskowitz spectrum [12] and Hasselmann's directivity function [1,13,14].More specifically, an up-going pressure wavefield is generated by placing one or several secondary sources below the receivers assuming a homogenous (water) half space.The dow-ngoing wavefield is then formed when these upward travelling waves are scattered by the rough sea surface.To simulate a dual-sensor configuration, vertical particle velocity fields are also computed.Wavefield separation in 3D (e.g., [15][16][17][18][19][20]) is then applied to the modeled dual-sensor data to recover the up-going and down-going pressure fields.Subsequently, the separated wavefields are extrapolated upwards in small steps, and an imaging condition [5,6,17,21,22] is applied at every step in order to recover the modeled 3D sea surface.The depth positions of the maximum image values define the sea surface.Finally, by using a sliding window through the data and repeating the imaging technique for each window, time variable changes of the sea surface can be accounted for [4][5][6].
In a typical 3D seismic acquisition, the wavefield is much more densely sampled in inline direction in comparison to cross-line direction and at such full 3D processing is not feasible.Alternatively, a 2D processing approach which implies streamerwise imaging of the sea surface is pursued.This imaging approach followed by interpolation is shown to give promising results after proper postimaging filtering.This is demonstrated by comparing the spectra of the imaged sea surfaces.Finally, based on the strategy developed using controlled data, two 3D dual-sensor field data examples were analyzed.The first field data example employed one source, while the second example employed dual-simultaneous sources towed at different depths.The former was acquired from deep water under moderate weather conditions, whereas the latter was acquired from relatively shallower water under marginal weather conditions.For the first time, the feasibility of imaging 3D time-varying sea surfaces from 3D dual-sensor streamer data is demonstrated.

Methodology
2.1.Time-Varying Sea Surface Imaging.Computing scattered data from a time-varying 3D sea surface is achieved by employing the 3D Kirchhoff-Helmholtz forward modeling tool discussed in Appendix B. In order to recover the sea surface information from the recorded scattered data, the backward approach could be used.However, our interest is to further validate the sea surface imaging technology proposed by [3][4][5] which is applicable to field data.It is reintroduced here in 3D in order to bring the present work into context.The technique involves wavefield separation, wavefield extrapolation, and an adequate imaging condition.

Wavefield Separation and Extrapolation.
Collocated dual-sensor systems simultaneously measure pressure wavefield and vertical component of the particle velocity wavefields.This measurement can be separated into the component up-going and down-going pressure or particle velocity wavefields if properly combined in data processing [4][5][6][15][16][17].In a homogenous medium, assuming plane wave representation, the total pressure wavefield P is recorded by hydrophones: where Ũ and D are, respectively, the up-going and downgoing pressure wavefields in frequency-wavenumber domain,   and   are, respectively, the inline and cross-line spatial wavenumbers and obey the dispersion relation ),   is the vertical wavenumber,  is angular frequency, and tilde designates the frequency-wavenumber transformed wavefields [16,18,23].Exploiting now the relationship between vertical particle velocity and pressure gradient, the velocity wavefield recorded by motion sensors is then in frequency-wavenumber domain [4,16,17]: where  is the density of the medium (water in this case).
Since the positive -axis is pointing downwards, the upgoing and down-going wavefields are both recorded with negative polarities for the velocity sensor but with, respectively, positive and negative polarities for the pressure sensor.This is because of the negative reflection coefficient at the sea surface and the fact that particle velocity field is a vector quantity, while pressure field is a scalar quantity.Thus, the up-going pressure wavefield is obtained by rearranging (2) and subtracting from (1): whereas the down-going pressure wavefield is obtained by adding the two equations: The up-going and down-going wavefield components of the vertical particle velocity wavefields can also be similarly obtained [4].In a selected window of traces, the separated upgoing and down-going pressure or vertical particle velocity wavefields are then moved (extrapolated) to the sea surface in small steps upwards towards the sea surface.The decomposed wavefields are extrapolated towards the sea surface, by applying proper phase shifts [4,17,24,25] where the up-going wavefield is forward propagated in time, while the down-going wavefield is backward propagated in time.Consequently, the difference in their travel times continues reducing until it becomes zero at the sea surface location where the arrival times and the amplitudes of the two wavefields are identical (for a special case of a flat sea surface), except for that they are 180 ∘ out of phase.

Imaging Condition.
The positions of coincidence of the extrapolated up-going and down-going wavefields are determined by applying an adequate imaging condition at each extrapolation depth.This follows from the fact that conventional seismic imaging is based on the assumption that wavefields originating from a seismic source propagate to and interact with a discontinuity (i.e., the sea surface in this case) as an incident wavefield (or the up-going wavefield) before returning to the sensor(s) as a reflected seismic wavefield (or the down-going wavefield) [25,26].Therefore, the two wavefields kinematically coincide at the discontinuity.Thus, at every extrapolation depth, the classical / imaging condition could be applied [17].However, this imaging condition becomes unstable for small values of the up-going wavefield.
On the one hand, least-square imaging condition involving a division of the cross-correlation of the up-going and down-going wavefields by the autocorrelation of the up-going wavefields at each extrapolation depth [4,5,17,21,22] has been shown to be robust [21]: where   and   , respectively, denote the extrapolated up-going and down-going wavefields in frequency-space domain,  is the current extrapolation depth,  and  define the lateral coordinates of each receiver, and bar indicates conjugation.The coordinates of the maximum image point values define the sea surface position at that point.The extrapolation depth associated with the maximum image value (, , ) for a given sensor corresponds to the sea surface elevation () at the sensor's spatial location (, ).A plot of these extrapolation depths against these spatial locations gives a time stationary image of the sea surface corresponding to the selected time window.

Processing Strategy for 3D Data
In order to establish a robust processing strategy of 3D field data, we first analyze controlled scattered data associated with synthetic sea surfaces.Thus, a gridded sea surface model was simulated based on (A.1) in Appendix A. See the following references on simulation of realistic sea surfaces [1, 5, 10, 12-14, 27, 28]. Figure 1(a) shows frozen time frame of the sea surface for a wind speed of 17 m/s blowing at an angle of  = 45 ∘ , whereas Figure 1(b) shows the associated ensemble averaged spectrum.An ensemble averaged spectrum was used to represent a stochastic process of ocean waves; nevertheless, the truncated sea surface area employed in the synthetic study is chosen in order to properly investigate the effects of truncation in field data.By exploiting deep water dispersion relation, the sea surface input parameters such as its speed and direction can be recovered from the spectrum.Speed of approximately 19 m/s and direction of 45 ∘ are obtained from the spectral peak by employing (A.5) and (A.8), respectively.These parameters correlate well with the input values.A water-filled half space with a free surface boundary was assumed.Receiver grids containing 16384 channels and separated by 6.25 m in the inline and cross-line directions were placed at a depth of 50 m below the sea surface.The upward travelling wavefield was simulated by assuming a "secondary" source at a depth of 700 m below the middle receiver.The source pulse was a Ricker wavelet with a centre frequency of 90 Hz and temporal sampling interval of 2 ms.In the simulations, the -direction was taken as the inline direction, while the -direction was taken as the cross-line direction.
The sampling of the sea surface requires at least five discretized points per acoustic wavelength of the wavefield impinging the sea surface.In addition, to minimize boundary effects, the sea surface area has to be extended beyond the region of interest.Thus, a sea surface covering an area of about 1.1 km 2 was employed in computing the data.The sea surface is gridded with a grid size of 3 × 3 m 2 implying total scattering points on the surface of 116281 for one time instant of the time-varying sea surface.Due to the large number of receivers (16384) only one source was used to generate the upward travelling waves (i.e., static sea surface).The computational time would otherwise become prohibitive.In [5] 2D time-varying sea surfaces were modeled based on the 2D Kirchhoff-Helmholtz formulation.The time-varying changes of the sea surfaces were obtained by applying the imaging technique in a sliding window mode.Our goal here is to analyze 3D sea surfaces and further validate our imaging technique.Figure 2 shows a time slice taken at 0.5 s of the generated total pressure wavefield for streamers at 20 m depth.Observe that the up-going waves (the outer circle) did not interact with the sea surface and only show the primary wavefield with the associated geometrical spreading effects.However, the down-going waves (the inner circle) represent the scattered wavefield with the trailing coda from the complex sea surface.Notice also that the up-going and down-going wavefields are about 180 ∘ out of phase.The data is generated in a homogenous medium, thus the circular event.Before applying the imaging technique as discussed earlier, the vertical velocity wavefields were obtained from the generated pressure wavefields by use of the constitutive equations.Then wavefield separation was performed to obtain the separated up-going and down-going pressure or velocity wavefields.Though this is not needed for synthetic data, the above procedure was followed in order to study how possible errors resulting from the wavefield separation may distort the synthetic data.Notice the good match between the modeled and imaged sea surfaces.However, due to the limited frequency band employed in the modelling (i.e., in data generation), the imaged sea surface is a filtered version of the true sea surface model (i.e., high frequency variations of the sea surface were not resolved).Our goal here is to demonstrate that the imaging technique is able to resolve the predominant structures of a given sea surface that affect seismic wavefields of a dualsensor type of acquisition.Spectral analysis of the truncated sea surface model shown in Figure 3(a) was carried out, and, from the peak of this spectrum, the recalculated sea surface wind speed and direction are 22 m/s and 63 ∘ , respectively.The estimated speed of the main-wind-driven events obtained from the spectrum of the truncated sea surface model is close to the true value (the true wind speed obtained using (A.7) is 19 m/s).Also, the dominant wave direction is not that far from the chosen wind direction of 45 ∘ .Since these estimates are calculated from a rather limited sea surface area, the slight differences between the true sea surface wave speed (and direction) and that obtained from the truncated sea surface can be attributed solely to the truncation error.Figure 3(c) shows the spectrum of the imaged sea surface.The peak of the spectrum also corresponds to sea surface wave speed and direction of, respectively, 22 m/s and 63 ∘ .Hence the ability of the imaging technique to image 3D realistic sea surfaces is demonstrated as the peaks of the spectra of the truncated, modeled, and imaged sea surfaces show.the requirements of the processing and the freedom allowed by the survey design.When streamers are coarsely spaced, artefacts resulting from aliasing put a limit to the bin sizes (in the cross-line direction).Thus, during sea surface imaging, a practical way of processing 3D narrow azimuth streamer data is to apply the imaging technique in a streamerwise manner (i.e., by applying wavefield separation, extrapolation, and imaging condition to one streamer at a time).The imaged sea surface profiles along each streamer corresponding to the same observation time are juxtaposed and interpolated (using, e.g., bicubic interpolation) to give a 3D sea surface.By moving the imaging window down the data and repeating this procedure, sea surface images for different observation times are obtained.Finally, by concatenating the sea surfaces, a time-varying sea surface can be obtained.Figure 4(a) shows the same imaged sea surface as in Figure 3(b), but this time the imaging technique has been applied streamerwise.Observe that the central parts of the sea surface being well illuminated are properly resolved in Figure 4(a) (since the "secondary" source generating the upward travelling waves is placed in the middle of the image).The differences in image quality are also due to the fact that the 2D extrapolation of the wavefields is only accurate at the vertical plane containing the "secondary" sources.Further away from this vertical plane, edge effects combined with wavefield separation and extrapolation errors distort the image at its upper and lower borders (see Figure 4(a)).Nevertheless, by comparing Figures 3(a) and 4(a), one can see that many of the main features of the modeled sea surface are still present.

2D or Streamerwise
Ignoring the cross-line angles in the streamerwise processing of 3D data introduces strong "DC errors (i.e., high amplitude values at wavenumbers close to zero)" which renders spectral analyses impossible (if not corrected).These DC values are associated with areas of high cross-line angles which are incorrectly extrapolated and tend to image continuously shallower with increasing cross-line angles (e.g., Figure 4(a)).As a pragmatic correction attempt, we subtracted the mean value of each imaged sea surface profile corresponding to a given time before juxtaposition and interpolation.Figure 4(b) shows the imaged sea surface obtained after the mean values of the streamerwise processed data have been removed.Figure 4(c) shows the spectrum of the imaged sea surface with DC corrections.After such corrections, it can be seen that the spatial wavenumber values corresponding to the peak are close to those obtained from 3D imaging (see Figures 3(b), 3(c), 4(b), and 4(c)).Nevertheless, truncation artefacts are still present in both spectra.These are seen as large amplitude values especially close to the zero wavenumber limits and are caused by the limited sea surface area.

Effects of Sparsely Spaced Streamers.
The latest advancements within 3D streamer technology allow seismic vessels to tow up to 20 streamers.However, the streamers are still sparsely spaced due to the complex relationship between load, speed, the limited frequency-band of the source, and other survey and navigational parameters associated with the particular seismic vessel and acquisition area.Even if all these parameters are balanced such that denser streamer spacing is permitted, maintaining densely spaced streamers can be challenging.Possible future streamer spacing is not believed to be less than 25 m.In order to investigate the effect(s) of feasible streamer spacing in the imaged sea surface, the streamer spacing used in the synthetic data was gradually increased from 25 m up to 100 m, and the imaging technique was applied streamerwise in each case.As before, the mean of each imaged sea surface profile was removed before interpolation.Figure 5(a) shows the imaged sea surface using the same survey configuration as before but using a streamer spacing of 25 m (i.e., keeping every 4th streamer of the original dense survey used in the modeling).Observe that the imaged sea surface in Figure 5(a) retains most features of Figure 4(b), albeit at a lower resolution.Using a streamer spacing of 50 m (i.e., keeping every 8th streamer) showed that the main features are still present but at a lower resolution when compared to Figure 5(a).Finally, a spacing of 100 m was used (i.e., keeping every 16th streamer), and the imaged sea surface shown in Figure 5(b) is similar to applying a harsh low pass filter to the imaged sea surface of Figure 4(b).Nevertheless, the main features are still present as indicated by the spectrum shown in Figure 5(c) (compare with Figure 4(c)).Thus, in terms of imaging, the effect of sparse streamer spacing is a reduced resolution as can be observed by comparing

Dual-Sensors and Conventional
Source.The first field data analyzed was acquired offshore Bahia, Brazil, by PGS using a 3D setup of dual-sensor streamers.The dual-sensor setup comprises collocated pressure and vertical particle velocity sensors that allow the up-going and down-going pressure or velocity wavefields to be separated during data processing.A total of 10 streamers separated by 100 m were towed at a depth of 15 m.Each streamer contains 648 collocated sensors separated by 12.5 m.The source is a conventional source consisting of arrays of air guns towed at a depth of 7 m. Figure 6 shows a shot record of the measured total pressure wavefield for the outer (streamer 1).The events at the upper part of the figure are primary events, while the weaker events at the lower part are the first water bottom multiples.Between these two dominant zones of events, primary arrivals with poor signal-to-noise ratio (S/N) can be seen.The time of the first arrivals shows that these data were acquired in deep water.After wavefield separation and extrapolation to the data acquisition datum level, a plot of the total pressure wavefield obtained by summing the up-going and the downgoing wavefields shows residue amplitudes suggesting that the sea surface has varied spatially and temporally during the data acquisition.
The data was processed in the streamerwise manner because of the limited number of streamers used in the acquisition and the characteristics of the source signature associated with the survey.The length of the imaging window was 1 s, and it was shifted 10 ms each time in order to recover the time-varying changes of the sea surface.Sea surface profiles corresponding to the same time were interpolated to obtain time-varying 3D sea surfaces.Streamer depth fluctuations were corrected for [5], and DC errors were removed by subtracting the mean value of each sea surface profile before interpolation.Figures 7(a) and 7(b) show a time frame of the imaged sea surface based on, respectively, the primary events and the deeper multiples of the first shot record.Observe that the sea surface shows spatial variations in both and -directions with the main direction of the wind-driven events appearing to be parallel to the direction.
From visual inspection, the statistical properties of the imaged sea surfaces appear to be similar.The changes between the sea surfaces of Figures 7(a) and 7(b) appear to be generally smooth as expected since the time difference (about 2.5 s) between the images is small compared to the main period of the sea surface (about 14 s obtained from spectral analysis).Time frames of the sea surfaces from the second shot imaged based on the primary events and the first multiples show similar characteristics as Figures 7(a) and 7(b).In particular, the difference between the sea surface from the primary events of the second shot and that from the first multiples of the first shot is small implying a smoothly varying sea surface as expected.Considering the difference in time between the two consecutive shots (about 12 s) one can conclude that the imaged sea surfaces (i.e., Figures 7(a) and 7(b)) are acceptable since the sea surface general characteristics appear to be retained (by examining different time frames of the imaged sea surfaces).However, due to windowing, the limited frequency band of the source, and the coarse spacing of the streamers, the imaged sea surfaces represent only filtered versions of the true sea state.
Spectral analysis of the imaged sea surfaces was carried out by averaging over all the available spectra for the first shot, and this was compared with that from the second shot (further analysis could be carried out by taking more shots).Since the imaging technique is based on the strength of the reflected data present in the imaging window, this implies that the area of low reflection amplitude between the primary and the multiple events was avoided (see Figure 6) in obtaining the time-varying sea surface.An ensemble averaged spectrum was then obtained by summing over the spectra of the time-varying sea surface.The number of time windows over which the ensemble average spectrum was taken is limited (337 time windows), and truncation errors falling close to zero wavenumbers in the spectrum were muted.Figures 8(a) and 8(b) show the averaged spectrum of the imaged sea surfaces corresponding to, respectively, the first and second shots.Observe the good correlation between these spectra.This is expected since the shots are consecutive.The peak wavenumbers are the same in both plots and correspond to a sea surface wave speed of about 21 m/s with the main-wind-driven events moving along the -direction (i.e., along the streamers).
The dominant wave height of the imaged sea surface is about 2.5 m.The observer's log file from this survey reported an average sea state of 4-5 on the Beaufort Sea State Scale.Based on this scale, a sea state of 4-5 is moderate with a prevailing wind speed of about 6-11 m/s and a subjective wave height of about 1-3 m.The log file also shows that data acquisition was stopped when the subjective wave heights reached about 3.5 m.Visual inspection of the imaged sea surfaces shows that the wave heights obtained from the imaged sea surfaces are reasonable when compared to the Beaufort Sea State Scale.The sail line direction at the acquisition time of these shots was SSE, and the wind direction changed from E to SSE at the same time, as reported by the observer's log file.This implies that the main wave fronts of the sea surface waves are approximately perpendicular to the streamers.Based on this fact we may conclude that the estimated directions of the main-wind-driven events given in Figures 8(a) and 8(b) are plausible.
However, the sea surface wave speed obtained from the spectral analysis appears to deviate from the reported results (i.e., wind speed obtained from the observer's log file).This can be due to the following reasons: (i) the 3D field data were processed and imaged streamerwise, (ii) heavy interpolation of the imaged sea surface in the cross-line direction (in order to increase the accuracy of the Fourier transform), (iii) truncation of the imaged sea surface (i.e., the sea surface area used in the spectral analysis not being wide enough to contain the dominant wavelengths and wave heights a sufficient number of times in order for the sea surface to be considered statistically stationary), and (iv) uncertainty in the observer's log.

Dual-Sensors and
Dual-Simultaneous Sources.Imaged sea surfaces obtained in the previous section show poor resolution when events from areas between the primary events and the first multiples were used.This is because this zone is characterized by poor S/N (see Figure 6).This demonstrates that the ability of the imaging technique to properly resolve sea surfaces from a given data set is hinged on the strength of the reflected events present in the data.To further examine the method, a second field data set acquired by PGS was imaged.This survey used dual-simultaneous sources and was acquired offshore the North Sea under marginal weather conditions.The 3D survey setup consisted of 6 dual-sensor streamers separated by 100 m with each streamer containing 480 collocated dual-sensors separated by 12.5 m and towed at 20 m depth in order to reduce noise.The source consists of primary and secondary air gun arrays towed, respectively, at 5 m and 9 m depths and fired almost simultaneously (the timing difference between the two sources is 0.795 s).This enables more reflected events to be recorded within a given imaging window.
Figure 9(a) shows a shot record of the recorded total pressure wavefield from the outermost streamer.When compared to the measured total pressure wavefields in Figure 6, one can see easily that the field data from Brazil was acquired in a relatively deeper sea with a harder sea floor (presence of strong multiples in Figure 6), while the field data from the North Sea was acquired from a relatively shallower sea with a softer sea floor (absence of strong multiples in Figure 9(a)).However, the second train of events with strong amplitudes seen below the first ensemble of events in Figure 9(a) comes from the second source.The visible vertical stripes seen across the data suggest that the data was acquired under noisy conditions.As before, streamerwise processing of the data and imaging was applied.
Figure 9(b) shows the imaged sea surface taken around the primary events.Visual inspection of the imaged sea surface suggests that the sea surface condition at the time of acquisition was marginally rough (considering the dominant wave heights and wavelengths).The imaging technique was able to resolve the time-varying sea surface, despite the noise contamination, due to the dual-simultaneous sources used in the acquisition.However, the imaging window had to be restricted to the region in the data with high density of reflected events.Thus, by considering Figure 9(a) it follows that well-resolved sea surfaces cannot be obtained beyond 3 s from this particular example of dual-simultaneous sources.The ensemble averaged spectrum of the imaged sea surfaces was obtained as before.The spectrum was computed from consecutive time frames of the imaged sea surfaces obtained from data 0.3 s-2.5 s.The spectral peak indicates that the dominating wave direction is parallel with the streamers (i.e., in the -direction) with a prevailing wind speed of about 25 m/s.Such a value is considered a marginal sea surface wave speed according to the Beaufort Sea State Scale.The observer's log file reported that the wind was blowing westwards with a speed of about 14 m/s and also showed that the sail line direction is WSW.However, as discussed earlier the spectral peak value is likely to be contaminated by truncation.

Conclusions
By Applying wavefield separation to collocated dual-sensor towed streamer data followed by upwards extrapolating of the resulting up-going and down-going wavefields and then applying an adequate imaging condition at every extrapolation depth, sea surface variations are imaged.Repeating the process in a sliding window gives sea surface image at a different time instant.Using existing 3D sea surface modelling formulation and combining it with the 3D extension of the time-domain Kirchhoff-Helmholtz integral including the Kirchhoff approximation, we were able to model acoustic wavefields scattering from realistic time-varying sea surfaces in an efficient manner.Spectral analysis of the sea surfaces shows that the peak wavenumbers of the spectra are related to the wind speeds and wind directions.These data enabled us to investigate the performance of the sea surface imaging technique in case of realistic 3D sea conditions.The robustness of the imaging algorithm in 3D was demonstrated employing synthetic data generated using the realistic Pierson-Moskowitz sea surfaces for different feasible 3D streamer configurations.The effect of applying the imaging technique in 2D was shown to be poor resolution of the imaged sea surfaces as the angle of incidence increases.However, a significant improvement could be obtained if the mean values of each imaged sea surface profile were subtracted.Truncation artefacts introduced in the spectra of the imaged sea surfaces due to a limited aperture were reduced by using ensemble averaged spectra.Sea surface wave speeds and directions recovered from such spectra of the imaged sea surfaces correlate well with the true values.Nevertheless, these values are still affected by errors due to truncation and 2D processing.
Two different field data sets from different 3D seismic surveys and different sea surface conditions were processed.The first field data was acquired in a deep water area with hard sea floor using collocated dual-sensor streamers and conventional sources under moderate sea surface conditions.The second field data was acquired, in relatively shallower water with a softer sea floor, using dual-simultaneous sources towed at different depths under rough weather conditions.Time-varying sea surface images obtained from the data confirm that the ability of the imaging technique to resolve sea surfaces is hinged on the strength of the reflected events present in the sliding imaging window (or alternatively to the S/N ratio in the window).In addition, given the typical source frequency band characteristics and sparse spacing of the streamers in the cross-line direction, the imaged sea surfaces represent filtered versions of the true sea conditions.
We have demonstrated sea surface imaging from dualsensor data by modeling, imaging, and analyzing realistic sea surfaces in a systematic manner using 3D acquisition geometries.Ongoing research effort is geared towards improving the imaging algorithm to be able to resolve sea surfaces in areas where the data show poor S/N ratio.The imaged sea surfaces could be "already made" input to many seismic processing schemes.For example, based on the imaged sea surface information, static corrections could be now applied during data processing in a proper manner.In addition, streamer depth bias can also be corrected by carrying out spectral analysis of the imaged sea surface.Access to the imaged time-varying sea surfaces could provide insight about energy distribution on the sea surface during seismic data acquisition.This could help in making strategic decisions when the acquired data are processed.In time-lapse seismic, where slight changes in the reservoir are pertinent information, imaged sea surfaces could be an indispensable tool.Another possible application of the results is in the statistical analyses of the reflection coefficient estimates obtained by imaging the sea surface which is a prioritized ongoing research topic.

Appendices A. Realistic Time-Varying Sea Surfaces
We outline the general method of computing realistic timevarying sea surfaces similar to that used by [1].Spatiotemporally varying sea surface height function (  ,   , ) is given as [1,27]  (  ,   , ) where (  ,   ) is isotropic of the sea surface combined with complex symmetric random numbers and some scaling factors (see [5,10]).  = 2/  and   = 2/  are, respectively, the and -components of the absolute wavenumbers   = √ 2  +  2  , while   and   are the and -components of the sea surface area.The wind wave model implied in the term (  ,   ) in (A.1) is taken to be the Pierson-Moskowitz spectrum (in this work).The spectrum in 3D is [12,28]  (A.4) The normalization constant where Γ is the gamma function.In (A.3)   is the frequency which relates to the wavenumber   through deep water dispersion relation given as where   is the peak frequency.This implies that each space harmonic component of the sea surface is made to move with a definite phase velocity (i.e., sea surface waves are dispersive).The sea surface moves forward (if Ω is negative) or backwards (if Ω is positive) as discussed by [27].
The angle between the wind vector and the wave propagation direction measured from north in clockwise direction (i.e., from   -axis to   -axis) is given as Thus, by changing  the main-wind-driven events of the sea surface can be steered in a preferred direction.Since ocean waves are dispersive, the spectrum is such that waves of longer wavelengths travel faster relative to waves of shorter wavelengths, and the peak (i.e., main-wind-driven events) occurs for a wave whose speed is comparable to the prevailing wind speed and is given as [28] (Ω peak ) The direction of the main-wind-driven events is recovered from the spectrum of the sea surface by applying ) .
(A.8) Thus, by performing spectral analysis of the sea surfaces (i.e., the Fourier transform of the sea surfaces), the spectral peaks can be used to approximately estimate the speeds (using (A.5)) and directions (using (A.8)).

B. The Kirchhoff-Helmholtz Integral for Time-Varying Boundaries
Scattered wavefields are computed from the simulated timevarying sea surface height functions employing the timedomain Kirchhoff-Helmholtz integral.Orji et al. [4] computed scattered pressure wavefields from a time-varying 2D realistic sea surface employing the 2D time-domain Kirchhoff-Helmholtz forward modeling formulation where the source and the receiver are taken to be at the far field of the sea surface.Here, 3D extension of the forward modeling formulation is outlined.In order to establish the notation being used here, the following quantities are introduced with reference to Figure 10.(, ) represents the surface height function; [  ,   , (  ,   )] defines the position of a scattering point on the surface;  ⇀   defines a vector from the origin to the running scattering point;  ⇀  defines a vector from the origin to a fixed receiver position; ⃗   −  ⇀  is a vector from a given receiver position to the running scattering point;  ⇀   is a vector from a source position to the receiver;  ⇀   defines a vector from the origin to the fixed source position;  ⇀  defines the vector from the fixed source position to the running scattering point; the unit vectors n and ρ, respectively, denote the normal to the surface and the unit vector direction of the incident field at [  ,   , (  ,   )]; the obliquity factor is given by cos  = n ⋅ ρ ≡ (  ,   ) (see Figure 10).The up-going and the corresponding down-going pressure wavefields caused by scattering from the sea surface are computed from the 3D Kirchhoff-Helmholtz integral with the Kirchhoff approximation included (i.e., valid when the sea surface is considered locally planar on the scale of the acoustic wavelength): ).The validity of the Kirchhoff-Helmholtz integral with the Kirchhoff approximation included has been demonstrated by many authors (e.g., [5,7,10,11]).

Figure 1 :
Figure 1: (a) A time frame of the time-varying directional Pierson-Moskowitz sea surface for a wind blowing with a speed of 17 m/s at an angle of 45 ∘ between the positive and -directions.(b) The ensemble averaged spectrum corresponding to the sea surface conditions shown in (a).The peak of the spectrum is indicated.
Figure 3(a)  shows a surface plot of a selected part of the modeled sea surface from Figure1(a) to be imaged (the area covers approximately the

Figure 2 :
Figure 2: Time slice of the total pressure wavefield taken at time 0.5 s for hydrophones placed at 20 m depth.

Figure 3 :
Figure 3: (a) Surface plot of the selected part of the modeled sea surface from Figure 1(a) that is to be imaged, (b) the imaged sea surface, and (c) the spectral plot of the imaged sea surface shown in (b).The indicated peak of the spectrum corresponds to these spatial wavenumber values:   = 0.01704 and   = 0.00852.The imaging technique was applied in 3D.

InternationalFigure 4 :
Figure 4: (a) Streamerwise imaged sea surface covering the same area as in Figure 3(a), (b) the same sea surface as in (a) but with the mean of each sea surface profile subtracted, and (c) the spectral plot of the imaged sea surface shown in (b).The indicated peak of the spectrum in (c) corresponds to   = 0.01704 and   = 0.00852.The imaging technique was applied in 2D (i.e., streamerwise with each streamer taken in the inline direction).

Figure 5 :
Figure 5: Streamerwise imaged sea surface as shown in Figure 4, but now the spacings between the streamers are, respectively, (a) 25 m, (b) 100 m, and (c) the ensemble averaged spectrum of the imaged sea surface in case of 100 m spacing.The indicated peak of the spectrum corresponds to   = 0.01704 and   = 0.00852.

Figure 6 :
Figure 6: Total pressure wavefield at 15 m depth for the outermost streamer.

Figure 7 :
Figure 7: Frozen time frame of the imaged time-varying sea surface based on (a) the primary events of the first shot and (b) the first multiples of the first shot.

Figure 8 :
Figure 8: (a) Ensemble averaged spectrum based on all the imaged sea surfaces (i.e., from the primary events to the first multiple events after DC correction) of the first shot and (b) the ensemble averaged spectrum from the second shot.The indicated peaks of the spectra each correspond to   = 0.02104 and   = 0.

Figure 9 :
Figure 9: (a) A shot record of the measured total pressure wavefield at 20 m depth for the outermost streamer and (b) the imaged sea surface taken from the area around the primary event.

[
Figure 10: A sketch showing the coordinates of the source , the receiver , and the running scattering point [  ,   , (  ,   )] on the sea surface (, ).