Reconstruction of the Frequency-Wavenumber Spectrum of Water Waves With an Airborne Acoustic Doppler Array for Noncontact River Monitoring

This work presents a novel method to reconstruct the frequency-wavenumber spectrum of water waves based on the complex acoustic Doppler spectra of scattered sound measured with an array of microphones. The reconstruction is based on a first-order small-roughness-amplitude expansion of the acoustic wave scattering equation, which is discretized and inverted by means of a singular value decomposition. An analogy of this approach to the first-order Bragg scattering problem is demonstrated by means of a stationary phase expansion. The approach enables the reconstruction of the dispersion relation of water waves when the ratio between roughness height and acoustic wavelength is less than 0.1, and when the surface wavelength is larger than 1/2 of the acoustic wavelength. The method is validated against synthetic data and data from laboratory and field experiments, to demonstrate its applicability to 2-D and 3-D complex patterns of water waves, and specifically to the surface deformations that arise naturally in a turbulent open-channel flow. Fitting the reconstructed data with the analytical dispersion relation enables the noncontact estimate of the underlying flow velocity for hydraulic conditions where the coexistence of different types of turbulence-forced and freely propagating water waves would limit the accuracy of standard noncontact Doppler velocimetry approaches, paving the way for robust and accurate noncontact river monitoring using acoustics.


I. INTRODUCTION
T HE properties of water waves measured by means of scat- tered electromagnetic or acoustic signals have been used to monitor currents [1], bathymetry [2], and significant wave height and period [3] in the oceans, and surface flow velocities in rivers [4], [5].In most cases, the dominant mechanism that influences the scattered acoustic or electromagnetic signals is Bragg scattering [6].The scattered signal is selectively sensitive to water waves with the Bragg wavenumber κ B = k 0 [cos(φ i ) + cos(φ r )], where k 0 = ω 0 /c 0 is the wavenumber modulus of the incident signal, ω 0 is the corresponding angular frequency, c 0 is the speed of light or sound (as appropriate), and φ i and φ r are the elevation angles of the wavenumber vectors of the incident and scattered waves relative to the horizontal plane.For a monostatic setup where the signal is emitted and received at the same location, κ B = 2 k 0 cos(φ i ).The frequency of the sound or electromagnetic wave reflected from the rough moving water surface changes by an amount ω proportional to the wave speed (Doppler shift).Then, the Doppler spectrum of the scattered signal at the frequency ω 0 + ω is proportional to the amplitude of the water waves with wavenumber κ B , where ω = κ B c and c is the phase speed of the Bragg waves on the rough surface.Measurements of the Doppler spectra provide an estimate of the speed of water waves, which can be used to infer the velocity of the underlying flow [7].
Applications of Doppler spectra analysis for the noncontact measurement of flow velocities in rivers or open-channel flows employed both radar [4], [5], [8], [9], [10] and acoustic [11] signals.These applications can suffer from large uncertainties because of the complex dynamics of the water waves when they interact with a turbulent flow.Turbulent forcing can modify the speed of water waves [12] or generate additional wave patterns [13] making the identification of the Doppler velocity and their link with the underlying depth-averaged flow velocity nontrivial [4], [11].Measurements with a bistatic geometry (separate source and receiver) allow focusing on longer and more predictable waves while increasing the signal-to-noise ratio, but at the expense of the ease of data interpretation [8], [14].The main obstacle to the interpretation of the Doppler 1558-0644 © 2024 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
spectra for river monitoring is the fact that waves with different wavelengths, speeds, and/or directions of propagation can produce a peak at the same Doppler frequency [11].The frequency-wavenumber spectrum of the water surface fluctuations allows the separation of the contributions from the different types of waves (gravity-capillary waves and turbulence-forced patterns) that appear in a turbulent flow, according to their speed of propagation and dispersion relation [15], [16], [17], [18].The analysis of such spectra enables a more accurate and robust estimate of the time-averaged surface flow velocity [3], and also of the water depth [19].However, the calculation of frequency-wavenumber spectra requires measurements of the surface deformations in time and in space which is challenging.Numerous methods have been developed for the reconstruction of the ocean wave spectrum based on scattered radar signals, although their application to river flows is not straightforward.There is a direct relationship between the acoustic/electromagnetic wave Doppler spectrum and the water wave spectrum, which can be revealed by expanding the scattering equations with respect to a small ratio of surface amplitude over signal wavelength [20].Voronovich and Zavorotny [21] proposed a method to reconstruct the spatial spectrum of water waves based on the first-order expansion of the Doppler HF/VHF radar signal emitted and recorded by a moving antenna on an airplane, using a synthetic aperture approach.The detection of the Bragg peaks in the Doppler spectrum is facilitated in the case of radar scattering in the ocean, where measurements are usually performed in the Fraunhofer zone, i.e., for small values of the ratio between the effective size of the illuminated area and measurement range [22], and where the angle of incidence of the signal and the Bragg wavenumber are well defined.This facilitates the estimation of the wave speed but provides little information about the spectrum of scales other than the ones identified by the Bragg wavenumber.In turn, information about the broader water wave spectrum can be extracted from the continuous second-order terms of the backscattering Doppler spectra in the neighborhood of the Bragg peaks [23], [24].Such an approach has been applied mainly to sea surface scattering of high-frequency radar, and was later improved and extended to multiple radar systems to address an ambiguity in the wave travel direction [25], [26].These methods require prior knowledge of the dispersion relation (and therefore of the speed) of water waves to relate the temporal information in the Doppler spectra with the spatial spectrum of the rough surface, therefore, they are not applicable in river flows where the dispersion relation depends on the (unknown) flow velocity.Alternatively, the continuous spectrum of relatively long water waves can be estimated from the modulation of shorter Bragg waves, which results in a pattern of spatial and/or temporal variations of the Doppler velocity [27].The approach employs the linear water wave theory to express the modulation as a function of wave frequency.It assumes the velocity of the flow to be constant, and requires a capability to clearly identify the frequency of the Bragg peaks with a high spatial and/or temporal resolution.
More recent approaches attempted to directly reconstruct the instantaneous shape of the water surface based on simultaneous measurements of the scattered electromagnetic [28], [29] or acoustic [30], [31], [32], [33], [34], [35] signal at different receivers along an array.These methods, derived for static surfaces, neglect the effects of the surface dynamics on the scattered field (e.g., Doppler shift), which can limit the reliability of the reconstruction over longer periods of time.
Here, we propose a method for reconstructing the frequency-wavenumber spectrum of water waves from the first-order Doppler spectra of scattered acoustic signals without any prior knowledge of the surface dynamics.The method produces an estimate of the dispersion relation of the water waves that enables the characterization of the surface and flow dynamics even in turbulence-dominated environments such as rivers and partially filled pipes.Combining a single omnidirectional narrow-band acoustic source with an array of receivers arranged within a region of specular reflection, the method allows the characterization of the behavior in space and time of small-amplitude (relative to the acoustic wavelength) water waves.The method is derived in Section II and validated both numerically (Section III) and experimentally (Section IV) in conditions that are representative of a small river or open channel.The potential application for the noninvasive monitoring of river flow is demonstrated by estimating the time-averaged flow velocity based on the reconstructed spectra in Section IV-B.

II. DERIVATION A. Three-Dimensional Surface
Fig. 1 illustrates schematically the problem of acoustic scattering by a dynamically rough water surface.The surface roughness is defined as function ζ expressing a perturbation of the plane x − y, z = ζ (x, y), with z pointing upward.A harmonic signal with angular frequency ω 0 and wavenumber k 0 is emitted by an acoustic source located at S = (x s , y s , z s ) in the far-field with respect to the scattering surface.The signal scattered by the surface is recorded by N M receivers with co-ordinates M m = (x m , y m , z m ), m = 1, . . ., N M that are also set up in the far-field.The directivity of the source is modeled by the function D(φ s (x, y, z)), where φ s is the angle between the normal to the acoustic wavefront and the axis of the speaker, identified by the unit vector n s .For an acoustic wave with wavenumber vector k 0 , k 0 • n s = k 0 cos(φ s ).Here, D(φ s (x, y, z)) is approximated by the directivity pattern of a piston of radius a 0 with an infinite baffle as [36] where J 1 (•) is the Bessel function of the first kind and first order.It is assumed that.
1) The maximum surface elevation is small relative to the distance between the source and the surface, i.e., max{ζ (x, y)} ≪ |z − S|, and to the distance between the surface and the receivers, i.e., max{ζ (x, y)} ≪ |M − z|.
2) The maximum surface slope is small relative to the angle of incidence of the signal with respect to the horizontal plane, i.e., max{∇ζ (x, y)} ≪ tan(φ i ), where φ i = − sin −1 (e z • k 0 /k 0 ).
3) The Kirchhoff approximation [37] is valid, i.e., 2k 0 r c sin 3 φ i > 1, where r c is the curvature radius of the scattering surface.4) The deformation of the surface shape in time is slow relative to the speed of sound.Condition 2 is necessary for the linearization of the scattering equation with an acoustically rigid boundary and ensures that shadowing can be neglected [37].
Under these conditions, the scattered acoustic pressure at the location M m can be approximated by the Kirchhoff integral equation [22] p(M m , t) = exp(−iω 0 t) H 0 (M m , x, y) where s are the squared distances to the source and receivers from a point with coordinates (x, y, 0) and is the projection of the rough surface on the x y-plane.The kernel H 0 (M m , x, y) is a function of the geometry of the array that is independent of the shape of the surface and of time t Assuming a small relative surface deformation, k 0 ζ ≪ 1, and expanding (2) in a Taylor series, one obtains where It is further assumed that the scattering surface can be described as the inverse Fourier transform of its complex frequency-wavenumber spectrum ψ(κ x , κ y , ω), e.g., κ y , ω e i(κ x x+κ y y−ωt) dκ x dκ y dω (7) where ψ(−κ x , −κ y , −ω) = ψ * (κ x , κ y , ω) to ensure that ζ is real, ψ * is the complex conjugate of ψ, and κ = (κ x , κ y ) is the wavenumber of the rough surface.The Fourier transform in time of the scattered signal p(M m , t) up to the second order yields where is the spatial Fourier transform of the kernel H n .
The zeroth-order term, δ(ω − ω 0 ) Ĥ0 (M m , 0, 0), in (8) contributes only to the signal at the carrier frequency ω 0 .The integral that includes the first term Ĥ1 (M m , −κ x , −κ y ) can be approximated numerically: discretizing the wavenumber and frequency axes with regular grids, . ., N y /2; and ω = n ω, n = −N ω /2, . . ., N ω /2; and assuming a surface spectrum with an upper cutoff to satisfy the Nyquist-Shannon sampling theorem (ψ(κ, Equation ( 10) can be rewritten as a matrix product Pm,n = p(M m , n ω) is the complex discrete Fourier spectrum coefficient in the spectrum of the scattered signal corresponding to the frequency n ω and recorded by the mth microphone. (2), . . ., (N y ) ] T , with elements (b)  n,a = ψ(a k x , b k y , n ω − ω 0 ), is the complex discrete frequency-wavenumber spectrum of the scattering surface at the Doppler frequency n ω − ω 0 , reshaped as a 2-D matrix.
The frequency-wavenumber spectrum of the surface can be estimated by inverting (11).Note that the choice of the wavenumber discretization is arbitrary (despite some requirements that will be highlighted in Section III).However, in typical applications, the number of discrete wavenumbers N κ exceeds the number of physical receivers N m .In this case, the inversion can be achieved numerically, for example, by means of a singular value decomposition of the transfer matrix Ĥ through regularization.Here, we employed the Tikhonov regularization, and the regularization parameter was calculated based on the generalized cross-validation method.Krynkin et al. [32] applied a similar approach to reconstruct the shape of a surface directly, but neglecting the Doppler shift.
Once the complex frequency-wavenumber spectrum is estimated, it is possible to reconstruct either the frequency spectrum at the location (x, y) or even the instantaneous shape of the water surface by means of (7).

B. Two-Dimensional Surface
Assuming no variations of the surface elevation along the y-axis (ζ = ζ (x, t) with ∂ζ /∂ y = 0) and source and receivers situated on the x z-plane, the 3-D Kirchhoff integral equation ( 2) can be simplified as a single integral along the projection x of the rough surface on the x-axis by means of a stationary phase expansion along y [38], yielding (13) where . Equations ( 5)-( 12) remain valid after replacing ψ(κ x , κ y , ω) with

C. Bragg Scattering Interpretation
To clarify the range of applicability of the method, it is helpful to examine the integrals in ( 9) and ( 15) by means of a stationary phase expansion.The case of a 2-D surface [see (15)] is considered for simplicity, but the results are valid in a more general 3-D case.The superscript (2D) and subscript x in κ x are dropped for brevity.It is assumed that the scattering surface area x is large compared to the acoustic wavelength λ 0 and that the directivity D(φ s ) varies slowly along x.The phase of the integrand is The gradient of the phase equals zero at the stationary phase point, x 0 (M m , κ), which is given implicitly by The first term on the right-hand side of ( 17) identifies the point of specular reflection, which is a stationary phase point when κ = 0.It is noted that in a monostatic configuration (M = S, N m = 1) (17) yields 2k 0 (x 0 − x s )/r s = 2k 0 cos(φ i ) = κ, which enables to identify κ as the Bragg wavenumber.In the more general bistatic configuration, the expansion of ( 15) at x 0 (M m , κ) yields Equation ( 18) links the Fourier transform Ĥn of the kernel at the wavenumber κ with the location x 0 on the surface that satisfies the Bragg resonance condition for κ. (17) has solutions only for |κ| ≤ 2k 0 .Therefore, the Doppler spectrum as a first-order approximation is governed by roughness wavelengths larger than λ 0 /2.Let us consider a sinusoidal scattering surface with the wavenumber κ 0 , moving at the speed c, with the following monochromatic frequency-wavenumber spectrum where ψ 0 is a constant.Substituting into (8) yields a result first obtained by Barrick [39] in the context of radar scattering.The zeroth-order spectrum at the carrier frequency ω 0 is Ĥ0 (M m , 0), which is due to reflections at the specular point.At first order, the Doppler spectrum has two peaks, at the Doppler frequencies ω − ω 0 = ±κ 0 c.These are linked to scattering at two stationary phase points, x 0 (M m , −κ 0 ) and x 0 (M m , κ 0 ), respectively.Assuming (without any loss of generality) that the surface moves toward the positive x-direction, and that x m > x s , then the peak at the positive Doppler frequency, κ 0 c, is due to scattering from x 0 (M m , −κ 0 ) < x 0 (M m , 0) which is to the left of the specular point (closer to the source-see Fig. 1), and the peak at the negative Doppler frequency −κ 0 c, is due to scattering from Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.x 0 (M m , κ) is the location of the stationary phase point identified from the Bragg resonance condition for the mth receiver and for the wavenumber κ [see (17)].
x 0 (M m , 0) identifies the specular point for the mth receiver.The two areas shaded in pink and in blue identify respectively the portions of rough surface where Bragg scattering by waves with the wavenumbers −κ (delimited by dashed lines) and κ (delimited by dashed-dotted lines) is predominant.The two areas overlap in the area shaded in purple.
x 0 (M m , κ 0 ) > x 0 (M m , 0) (closer to the receiver).The spatial dependence of the Doppler spectrum is highlighted in Fig. 1, where the areas shaded in pink and in blue represent the areas of the surface where Bragg scattering is sensitive to waves with the wavenumbers −κ (producing the Doppler spectrum at frequency κc) and κ (producing the Doppler spectrum at frequency −κc), respectively.These areas overlap only in a small portion of the surface.Second-order terms in (20) contribute to the Doppler spectrum at the carrier frequency ω 0 , and at the second-order harmonics ω − ω 0 = ±2κ 0 c of the Bragg frequencies.Including three terms in the expansion would lead to additional peaks at the frequencies ω − ω 0 = ±κ 0 c ± κ 0 c ± κ 0 c, and so on.Higher order terms of the small perturbation expansion yield spurious harmonics in the estimated surface spectrum.For a broadband roughness spectrum, the harmonics are replaced by the autoconvolution of the roughness frequency spectrum, yielding an additional broadband contribution to the Doppler spectrum which is not easily separated from the first-order terms, unless the measurements are performed in the Fraunhofer zone.In this case, the Bragg wavenumber is uniquely defined by the geometry, and second-order Doppler inversion [23], [24] becomes a more favorable approach.

D. Scattering From Water Waves
The spectrum of water waves can be described as ψ(κ, ω) = (κ)δ(ω − κc), where c is the phase velocity defined by the gravity-waves dispersion relation where g is the acceleration due to gravity, B = ρgγ −1 κ −2 is the Bond number, which indicates the ratio of gravity and surface tension effects, γ is the surface tension coefficient, ρ is the water density, and d is the water depth.
In rivers and open-channel flows, the mean flow modifies the frequency of gravity-capillary waves as follows [16]: where U is the flow velocity (here assumed to be constant in space and in time).In addition to gravity waves, other types of surface deformations such as turbulent boils and vortex dimples, (see Muraro et al. [13]) are induced by turbulent structures.The frequency of these so-called "forced" waves [12] is equal to [16] The coexistence of multiple types of waves with different speeds and directions of propagation makes the wavenumber-frequency relationship nonunique, causing multiple peaks in the Doppler spectra of backscattered acoustic [11] or radar [4] signals, and resulting in the ambiguity of flow velocity estimations.
On the other hand, the knowledge of the frequencywavenumber spectrum of the surface and its comparison with the theoretical dispersion relations [see (22) and (23)] would allow a more robust estimation of the surface flow velocity.The approach has been applied extensively in oceanography [40], [41], and was recently extended to the noncontact monitoring of the velocity [42] and water depth [19] of river flows through the analysis of videos of their water surface.Specifically, Dolcetti et al. [19] employed an optimization algorithm to identify the flow parameters (velocity and water depth) that provide the best fit between the measured frequency-wavenumber spectrum and a synthetic spectrum that satisfies the theoretical dispersion relations, (22) and (23).The applicability of this latter approach to the spectra reconstructed from the scattered acoustic signal is tested in Section IV-B, to estimate the flow velocity.

III. NUMERICAL VALIDATION
The numerical tests were aimed at characterizing the performance of the method in ideal conditions.For the numerical validation of the method, random surface realizations with different wave amplitudes and with different spectra were constructed by means of Fourier synthesis [37].The initial surface shape was evolved in time according to (21) (i.e., without mean flow), assuming infinitely large depth  [11], [14].The surface parameters used for the simulations in the present work ensured the validity of the small-roughness-amplitude approximation.
The Doppler spectrum calculated for the sinusoidal surface in Fig. 2(a) has a series of peaks at the harmonics of the frequency κc.These represent second-and higher order terms like those defined in (20).Since positive and negative Doppler frequencies relate to scattering at different locations on the surface (x 0 (M m , −κ 0 ) and x 0 (M m , κ 0 ), see (20) and Fig. 1), the Doppler spectrum is not symmetric with respect to ω −ω 0 = 0 [e.g., Fig. 2(d)].In Fig. 2(b) and (e), the crosses and the dashed line indicate the theoretical dispersion relation of the water waves for the surfaces with the sinusoidal and broadband roughness, respectively.The main ridges of the reconstructed frequency-wavenumber spectrum follow these lines closely, indicating a successful reconstruction.The wavenumber resolution of the reconstructed frequency-wavenumber spectrum is not governed by the aperture of the array, but by the distance between the first and last stationary phase points.Since this distance is different for positive and negative frequencies, the resolution is also different and generally higher at larger wavenumbers [e.g., Fig. 2

(b) and (e)].
The dashed-dotted lines in Fig. 2(c) and (f) indicate the position of the stationary phase points x 0 (M 1 , κ(ω)) and x 0 (M 30 , κ(ω)) calculated for the first and last receivers as a function of κ(ω) = ω/c, i.e., for the wavenumber that corresponds to the frequency ω according to the dispersion relation (21).The frequency spectra of the surface roughness were independent of x, whereas the reconstructed spectra appear to decay rapidly to zero outside of the region delimited by the first and last stationary phase points.This can be explained by ( 18) and ( 20), since only the portion of the surface where a stationary phase point can exist for the wavenumber κ contributes to the Doppler spectrum at the frequency −κc.Therefore, only this portion of the surface can be reconstructed based on the Doppler spectrum at such a frequency.Averaging the reconstructed surface frequency spectrum E(x, ω) along x must be done carefully because of the spurious dependence of E(x, ω) on x.Here, the averaging was performed by means of a frequency-dependent rescaling factor w(ω) = /[x 0 (M 30 , κ(ω)) − x 0 (M 1 , κ(ω))], i.e., Ē(ω) = w(ω) i E(x i , ω)/N x .Fig. 3 shows the results for synthetic data calculated with λ 0 = 23 mm and with surface amplitude varying between λ 0 /160 and λ 0 /10.Fig. 3(c) and (f) shows the average (and rescaled) reconstructed frequency spectrum of the scattering surface.As the relative amplitude of the water waves increased, the amplitude of the Doppler spectra in Fig. 3(a) and (d) also increased, as expected.An increase in the number of harmonics and a significant broadening of the Doppler spectra was also observed for the sinusoidal and broadband surfaces, respectively.The reconstructed frequency-wavenumber spectra in Fig. 3(b) and (e) relate to cases with λ 0 / max(|ζ |) = 10 [Fig.3(b)] and λ 0 /rms(ζ ) = 10 [Fig.3(e)].They both show some spurious ridges/harmonics at the frequencies ω = ±nκc, n = 1, 2, 3.These spurious harmonics cause an overestimation of the reconstructed frequency spectra amplitude, as seen in Fig. 3(c) and (f).
From the Bragg scattering interpretation in Section II-D and from the results of the numerical validation study, it is possible to draw the following preliminary conclusions regarding the range of applicability of the proposed methodology.
1) According to (9), the Doppler spectrum contribution at the first order at the Doppler frequency ω − ω 0 is related to the spectrum of the wave with the wavenumber κ that satisfies the Bragg resonance condition.From the point of view of the reconstruction of the surface spectrum, the smallest surface wavelength that can be reconstructed with an omnidirectional source is equal to |κ| < 2k 0 .
2) The contribution to the first-order Doppler spectrum of a surface wave with wavenumber κ is mediated by the value of the kernel H 1 at the location x 0 .Hence, the value of the reconstructed spectrum at each surface wavenumber κ is representative of the behavior of the surface in a limited region, whose size and location depend on κ and on the geometry of the acoustic setup [see Fig.

IV. EXPERIMENTAL VALIDATION A. Paddle-Generated Waves in a Wave Tank
A first series of experiments was performed in a 12 m long and 0.5 m wide rectangular wave tank in the LVV facility at the University of Sheffield.The tank was filled with water up to a depth of 1 m.The surface waves were generated by means of a purposely made CNC-controlled vertically oscillating wave-maker, in the absence of a unidirectional flow.The experimental setup was identical to that of the numerical simulations and photographs are shown in Fig. 4. The acoustic setup consisted of 30 1/4" microphones (G.R.A.S. 40PH) and loudspeaker (Visaton G25FFL).The incident acoustic signal had a frequency of 14 kHz with the wavelength of 24.3 mm.The wave-maker was installed at x = −756 mm behind the source and generated water waves propagating toward the positive x-direction.Time series of the surface elevation were measured at seven locations before and after the measurement area with a set of calibrated wave gauges.The amplitude of the water waves decreased exponentially with the distance from the wave maker.The wave amplitude at the array location was  The results of the analysis applied to experimental data are shown in Fig. 5.Despite the presence of some noise [see Fig. 5(b) and (e)] that was absent in the synthetic data sets, the spectrum reconstruction results based on experimental data are consistent with those observed in the numerical simulations.The dispersion relation of the water waves is easily identified by the maxima of the reconstructed frequency-wavenumber spectrum [Fig.5(b) and (e)].These relations follow those predicted by the theory [see (21)] over a significant range of scales.Compared to the frequency spectra of the water waves measured directly by wave gauges, the reconstructed frequency spectra in Fig. 5(c) and (f) reveal some influence of higher order effects (spurious harmonic peaks), which are consistent with the relatively large surface wave amplitude (λ 0 /max(|ζ |) = 18 and λ 0 /rms(ζ ) = 17, respectively).

B. Water Surface Fluctuations Induced by a Turbulent Flow
Turbulent open-channel flows display strongly 3-D water surface deformations.Therefore, the analysis was performed in 3-D with measurements obtained with a 2-D microphone arrangement on a plane.The acoustic device used in these experiments was a Simcenter Sound Camera (Siemens Digital Industries Software, Belgium), which is an array of MEMS microphones distributed around a video camera, and designed primarily for sound source localization methods, both nearand far-field.The Sound Camera was operated in two different configurations: a configuration with 49 MEMS microphones arranged along nine spiral-shaped rays with an aperture of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The laboratory experiments were performed in a recirculating hydraulic flume at the University of Sheffield.The dimensions of the flume were 15 m long and 0.506 m wide.The flume bed had a slope of 0.002 and was covered with a layer of spheres with a diameter of 25 mm.Four steady flow rates were tested: 5.3, 10.3, 13.2, and 20.0 l•s −1 (flow conditions Lab 1, Lab 2, Lab 3, and Lab 4, respectively).The water surface velocity was measured in each condition by timing neutrally buoyant particles across a known distance.The surface flow velocity varied between 0.30 and 0.59 m s −1 , while the uniform water depth varied between 42.2 and 115.1 mm across all the flow conditions (see Table I).The uncertainty of the surface velocity measurements, calculated as half of the maximum difference between all measurements in each flow condition, was below 5%.The Froude number of the flow (F = U/(gd) 1/2 ) varied between 0.47 and 0.56 (see Table I).The sound camera (in the base configuration) and source were installed 9.5 m downstream of the inlet.The center of the source was 250 mm away from the center of the camera and they were elevated 340 and 320 mm from the bed of the flume, respectively.The sound camera and source had the same inclination angle of approximately 60 • , with the sound camera facing downstream and the source facing upstream.vector components and frequency).Fig. 7(a)-(d) shows the spectral cross section at κ y = 0 and Fig. 7(e)-(h) shows the spectral cross section at κ x = 0, where κ x is in the streamwise direction (along the source axis plane) and κ y is in the transverse direction.The spectra in the κ y − ω plane are nearly symmetrical around κ y = 0, while those in the κ x − ω plane are strongly tilted upward toward positive κ x because of the advection of the waves by the mean flow [see (22) and ( 23)].Three spectra ridges are easily identifiable in Fig. 7(e)-(h).They are representative of the dispersion relation of gravity-capillary waves [upper and lower ridge, (22)] and of turbulence-forced surface deformations [central ridge, (23)], as demonstrated by the close proximity with the theoretical relations also plotted in Fig. 7.
Fig. 8 shows the normalized frequency-wavenumber spectra of the water surface reconstructed from the analysis of the scattered acoustic signals measured in the field.Compared to the laboratory data, the central ridge corresponding to turbulence-forced surface deformations is nearly absent in Fig. 8(a)-(c).The balance between gravity-capillary waves and turbulent-forced deformations at the free surface of a turbulent flow is the result of the equilibrium between the turbulent forcing and the stabilizing effect of gravity and surface tension, determined by the characteristic velocity and length scales of the flow [13], [44].While the Lab Conditions showed a combination of both types of waves which is indicative of a regime of weak turbulence, the lower Froude number in the Field Conditions is representative of a regime of gravity-dominated turbulence, where waves are expected to be predominant [44].The transverse spectra in Fig. 8(d)-(f) are not symmetric demonstrating the capability of the method to account for the velocity components in the xand y-directions.The more ragged ridges in the streamwise spectra for Field Condition 3 [Fig.8(c)] are likely to be an artifact of the normalization or higher order effects due to larger surface deformations for this flow condition.
Fig. 9 shows the reconstructed frequency spectra of the surface fluctuations in the laboratory [Fig.9(a)] and field [Fig.9(b)] experiments.These are compared with the spectra of the data recorded with wave gauges for the laboratory conditions in Fig. 9(a).The increase in the spectral energy from Lab Condition 1 to Lab Condition 4 is captured well by the reconstruction, although the energy at frequencies below 4 rad s −1 is largely overestimated.This is due to leakage of the zeroth-order term of the Doppler spectrum at the Doppler frequency χ − ω 0 ≈ 0. Similar behavior of the reconstructed frequency spectra can be seen in Fig. 9(b) for the field dataset, although the lack of direct measurements of the surface fluctuations in the field does not allow an evaluation of the accuracy of the reconstruction for this dataset.The relatively small slope of the reconstructed spectrum for Field Condition 3 suggests a noisier and less accurate reconstruction, which is confirmed by the more ragged appearance of the reconstructed frequency-wavenumber spectrum in Fig. 9(c) and (f), and which could be due to the larger amplitude of the surface fluctuations in this flow condition.
The reconstructed frequency-wavenumber spectra in the laboratory and in the field were fit with the theoretical dispersion relations [see (22) and ( 23)] using the MATLAB code kOmega [45], based on the work of Dolcetti et al. [19] to estimate the time-averaged flow velocity U = |U|.The measured and estimated velocities are reported in Table I.The difference  between the measured and estimated values is smaller than 10% for the laboratory data and 14% for the field data, which is comparable with a typical uncertainty in the direct flow measurements.

V. CONCLUSION
A novel method for reconstructing the frequencywavenumber spectrum of a dynamic water surface based on the Doppler spectra of the scattered acoustic field recorded with an array of microphones has been proposed.The method does not require any prior knowledge of the free surface dynamics to produce an estimate of the dispersion relation of the water waves, enabling the characterization of the surface dynamics in turbulence-dominated environments such as rivers or partially filled pipes.The spectrum reconstruction is based on a small roughness-amplitude expansion that can be interpreted in terms of first-order Bragg scattering theory.The limits of applicability of the method have been investigated with a numerical model.The shortest roughness scale that can be reconstructed is limited to half of the acoustic wavelength.The spectral resolution of the reconstructed spectra is a function of the surface wavelength and is inversely proportional to the maximum distance between the stationary phase points for all microphones in the array.
The method has been validated in the laboratory and in the field for 2-D and 3-D surface roughness patterns, demonstrating its capability to detect the presence of different types of waves (dispersive gravity-capillary waves and nondispersive turbulence-forced surface fluctuations) coexisting simultaneously on the surface.Fitting the reconstructed frequency-wavenumber spectra with theoretical dispersion relations for both types of surface deformations using an existing algorithm [45] has allowed for the estimation of the flow velocity with an error of the order of 10%-14%.By accounting explicitly for the complexity of the water surface dynamics in turbulent flows and by being based on a bistatic geometry to provide a high signal-to-noise ratio, the proposed method is a robust and compact alternative for existing noncontact river monitoring approaches.It could be rapidly deployed on unmanned aerial vehicles [46] to monitor floods and representative sections of river reaches.The approach could be extended to electromagnetic waves, complementing existing techniques for the reconstruction of

Fig. 1 .
Fig. 1.Schematic representation of the geometry of the problem.S and M m represent the position of the acoustic source and of the mth receiver, respectively.x0 (M m , κ) is the location of the stationary phase point identified from the Bragg resonance condition for the mth receiver and for the wavenumber κ [see(17)].x0 (M m , 0) identifies the specular point for the mth receiver.The two areas shaded in pink and in blue identify respectively the portions of rough surface where Bragg scattering by waves with the wavenumbers −κ (delimited by dashed lines) and κ (delimited by dashed-dotted lines) is predominant.The two areas overlap in the area shaded in purple.

Fig. 2 .
Fig. 2. Examples of spectrum reconstruction based on synthetic data for: (a)-(c) sinusoidal surface with λ 0 /max(|ζ |) = 80 and (d)-(f) broadband surface roughness λ 0 /std(ζ ) = 80.(a) and (d) Doppler spectrum at the first receiver, p(M 1 , χ ).(b) and (e) Reconstructed frequency wavenumber spectrum ψ(κ, ω).(c) and (f) Reconstructed frequency spectrum E(x, ω).The crosses in (b) indicate the actual wavenumber and frequency of the scattering surface.The dashed line in (e) indicates the dispersion relation of gravity waves.The dashed lines in (c) and (f) indicate the location of the stationary phase point x 0 (M m , κ(ω)) calculated for the first and last receiver of the array.
2(c) and (f)].3) The spectral resolution in wavenumber space varies with the wavenumber κ, and it is proportional to 2π/max(|x 0 (M n , κ)−x 0 (M n ′ , κ)|), i.e., it depends on the maximum distance between the stationary phase points obtained with (17) for all microphones in the array, as a function of κ [see Fig. 2(b) and (e)].

Fig. 4 .
Fig. 4. Photographs of the experimental setup with the wave tank in the LVV facility at the University of Sheffield.

Fig. 6 .
Fig. 6.Photographs of the measurement setup for the laboratory (a) and field (b) experiments with turbulence-induced water surface fluctuations.The flow direction is from left to right.
The acoustic signal radiated by the source had a frequency of 19 kHz.The instantaneous surface elevation was measured by means of three wave gauges positioned at a distance of 10 m from the flume inlet.The measured rms average of the surface fluctuations varied between 0.09 mm (λ 0 /rms(ζ ) = 199, condition Lab 1) and 0.25 mm (λ 0 /rms(ζ ) = 72, condition Lab 4).Fig. 6(a) shows the photograph of the experimental setup in the flume.Field measurements were performed on two separate dates in the River Loxley (Sheffield, U.K.) at three different locations within the same reach.The River Loxley is a river draining a catchment of 43.5 km 2 .The mean annual discharge recorded at Damflask Reservoir (approximately 1.5 km upstream of the measurement reach) is 0.563 m 3 s −1 , with 5% percentile of 1.512 m 3 s −1 and 95% percentile flow of 0.116 m 3 s −1 [43].The water depth and water surface velocity at each measurement site were: 350 mm, 0.36 m s −1 (Field Condition 1); 500 mm, 0.41 m s −1 (Field Condition 2); and 550 mm, 0.52 m s −1 (Field Condition 3).The estimated uncertainty of the surface velocity measurements between 7% for condition Field Condition 1 and 14% for Field Conditions 2 and 3.The Froude number of the flow varied between 0.19 and 0.22.The speaker and sound camera were held at a height between 500 and 700 mm above the water surface by means of a cantilever beam installed on two tripods [Fig.6(b)], with a distance of 700 to 750 mm between the speaker and the center of the array.The incident acoustic signal had a frequency of 17.5 kHz.The water surface fluctuations could not be measured directly in the field.Their rms average was estimated by integrating the frequency power-spectral-density spectra estimated acoustically, E(x, y, ω).The estimated rms(ζ ) varied between 0.66 mm (λ 0 /rms(ζ ) = 29, Field Condition 1 and Field Condition 2), and 4.63 mm (λ 0 /rms(ζ ) = 4.2, Field Condition 3), although the latter could have been overestimated with the proposed method.Fig. 7 shows the normalized frequency-wavenumber spectra of the water surface reconstructed from the acoustic signals recorded in the laboratory flume.The normalization was performed by dividing the spectra by their maximum at each frequency to compensate for the rapid decay at high frequencies and to facilitate the visualization over a wider range of scales.The reconstructed spectra are 3-D (two wavenumber

Fig. 9 .
Fig. 9. Frequency spectrum, E(ω), reconstructed from the experimental data.(a) From laboratory data.(b) From field data.(dashed lines): frequency spectrum of the surface fluctuations measured with wave gauges.