Virtual Head Waves in Ocean Ambient Noise: Theory and Modeling

The Green’s function retrieval in media with horizontal boundaries usually only considers the extraction of direct and reflected waves but ignores the virtual head waves, which have been observed experimentally from ocean ambient noise and used to invert for geometric and environmental parameters. This paper derives the extraction of virtual head waves from ocean ambient noise using a vertically spaced sensor pair in a Pekeris waveguide. Ocean ambient noise in the water column is a superposition of direct, reflected, and head waves. The virtual head waves are produced by the cross-correlations between head waves and either reflected waves or other head waves. The locations of sources that contribute to the virtual head waves are derived based on the method of stationary phase. It is the integration over time of contributions from these sources that makes the virtual head waves observable. The estimation of seabed sound speed with virtual head waves using a vertical line array is also demonstrated. The slope of the virtual head waves is different from that of direct and reflected waves in the virtual source gather; it is therefore possible to constructively stack the virtual head waves. The predictions are verified with simulations. VC 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1121/10.0002926 (Received 10 September 2020; revised 19 November 2020; accepted 26 November 2020; published online 18 December 2020) [Editor: Nicholas P. Chotiros] Pages: 3836–3848


I. INTRODUCTION
The Green's function estimate between two receivers in both open and closed systems can be extracted by crosscorrelating the field measured from sources that surround the receivers; this is referred to as Green's function retrieval 1,2 or seismic interferometry. [3][4][5][6]32 Experimental evidence of this result has been presented in helioseismology, 7 ultrasonics, 8 underwater acoustics, 1,9,10 and seismology. 11 In an environment with boundaries, the emergence of a Green's function estimate from correlation of ambient noise recordings was demonstrated theoretically. Lobkis and Weaver 8 demonstrated this principle using a normal-mode formulation in a reverberant ultrasonic cavity. However, the equilibration of normal modes is not a necessary condition. By cross-correlating multiple scattered waves in a freespace medium with embedded scatterers uniformly distributed, Snieder 12 proposed a stationary phase derivation to explain the extraction of the ballistic Green's function in seismology. The method of stationary phase was often combined with ray theory to explain the extraction of Green's function estimates in a medium with horizontal reflectors. For example, in a homogeneous medium with one horizontal reflector and without a free surface, Snieder et al. 13 showed the correlation of waves recorded by two receivers correctly yielded the Green's function (direct plus single-reflected wave) between the two receivers. In an ocean waveguide bounded above by a pressure release surface and below by a seabed with attenuation, Sabra et al. 9 formulated the time domain Green's function for time-averaged surface generated ambient noise cross-correlation, where the sources were modeled as point sources evenly distributed on a horizontal plane at a constant depth. Brooks and Gerstoft 14 described the relation between the stacked cross-correlations from a line of vertical sources, located in the same vertical plane as two receivers, and the Green's function between the receivers.
Recently, the virtual head waves have been observed from ocean surface generated noise both in simulation and experiment using vertical [15][16][17] and horizontal 15,18 arrays. The processing used is a generalization of the passive fathometer [19][20][21][22][23] and produces cross-beam correlations, [15][16][17] and it is equivalent to seismic interferometry techniques for delay and sum beamforming (but is not for adaptive beamforming). 15 The virtual head waves have the same phase speed as the real, acoustic head waves, [24][25][26][27] but the travel time is offset due to the cross-correlation method that the paths in common disappear and only the difference remains, thus the term virtual. The virtual head waves, also called spurious multiples, 13,14,28,29 are non-physical energy in Green's function estimates. 30 However, they are useful since a) Electronic mail: jie_li@sjtu.edu.cn, ORCID: 0000-0003-3066-1590. b) ORCID: 0000-0002-0471-062X. c) ORCID: 0000-0002-3487-6110. the travel time and angle of arrival of the virtual head waves can be used to invert for geometric and environmental properties. 15,17 Most research on Green's function retrieval in media with horizontal boundaries only mentions the extraction of direct and reflected waves between receivers 9,13,14 but not head waves or virtual head waves. The goal of this paper is to investigate the extraction of virtual head waves from ocean ambient noise in the simple case of a Pekeris waveguide with the method of stationary phase and ray theory. The differences between this work and literature mentioned above 9,13,14 are as follows: (1) besides direct and reflected waves, head waves are also considered between the noise sources and receivers and (2) instead of two receivers that are arbitrarily positioned in the x-z plane, the focus is on two vertically placed receivers in the water column, and it is extended later to the case of a vertical line array. Under these conditions, nine terms can be obtained by crosscorrelating ocean surface noise recorded at two receivers. The method of stationary phase is used to analyze each term, and it is proved that the virtual head waves are produced by three terms.
In the following, Sec. II proposes a simple model of the sea surface generated ambient noise and explains the extraction of virtual head waves for a vertically spaced sensor pair. Section III presents the details of seabed sound speed estimation using the extracted virtual head waves from a vertical line array. Section IV demonstrates the theory above with simulations. Section V contains the summary and conclusion.

A. Ambient noise cross-correlation function
The model geometry is depicted in Fig. 1. The water column is bounded above by a pressure-release surface and below by a semi-infinite bottom layer. The density and sound speed of the water and the bottom are given by q 1 , v 1 and q 2 , v 2 , respectively. Consider an infinite plane parallel to the surface and located below the surface at depth z s . 31 In this plane, let each noise source strength be Sðr s ; z s Þ, r s ¼ ðr s cos u s ; r s sin u s Þ; and u s 2 ½0; 2p. The two receivers V 1 and V 2 are on the z-axis with coordinates ðr 1 ; z 1 Þ and ðr 2 ; z 2 Þ, r 1 ¼ r 2 ¼ ð0; 0Þ; 0 < z 1 z 2 < Z. The field then becomes independent of azimuthal angle u. The pressure from noise source S to receiver V 1 is Sðr s ; z s ; xÞGðr 1 ; z 1 ; r s ; z s ; xÞd 2 r s ; (1) where Gðr 1 ; z 1 ; r s ; z s ; xÞ is the Green's function between S and V 1 . In this section, frequency dependence (x) is suppressed since the derivation is in the frequency domain. Based on ray theory, the full Green's function consists of three terms, where G D , G R , and G H are direct, reflected, and head waves between S and V 1 . Assuming that sources are close to the surface (z s % 0), there are only down-going waves from the source, but both down-going (h 2 ð0; p=2) and up-going (h 2 ½Àp=2; 0Þ) waves at receivers. The grazing angle h shown in Fig. 1 is down-going. Since G R and G H have upand down-going contributions they can be expanded to, where A D is the amplitude of the direct wave down-going from S and down-going to V 1 ; A R 6 , and A H 6 are those of the reflected and head waves down-going from S and up-(À) or down-going (þ) to V 1 ; z 1s ¼ ½z 1 ; z s is the depths of receiver and source, c wb ¼ ½Z; q 1 ; q 2 ; v 1 ; v 2 contains parameters for the water column and sea bottom, h c ¼ arccosðv 1 =v 2 Þ is the critical grazing angle. Due to the ocean surface and bottom, the reflected and head waves can bounce several times in the waveguide thus have various ray paths from S to V 1 .
Defining m 2 N 0 as the bounce number from the bottom interface at z ¼ Z that occur between S and V 1 ; G R 6 , and G H 6 are expressed as the superposition of wavefields from each ray path. For the reflected waves, m 2 ½1; 1, while for the head waves, m 2 ½0; l 6 , where m ¼ 0 corresponds to the case when the horizontal range between S and V 1 is less than the critical offset defined as, X 6 1 ðmÞ ¼ ð2mZ À z s 6z 1 Þcot h c , 17 with A H 6 ðm ¼ 0Þ ¼ 0. This leads to the largest number of head waves bounces that are up-or down-going from S to V 1 as FIG. 1. (Color online) Definition of coordinate system and geometric quantities. The noise source S is located on a plane at (r s cos u s ; r s sin u s ; z s ), and receivers V 1 and V 2 are located at ð0; 0; z 1 Þ and ð0; 0; z 2 Þ, respectively, within a waveguide of depth Z.
Note, this paper is aimed at discussing the location of sources that contribute to virtual head waves, therefore the amplitudes are not expressed accurately. The expressions of A D , A R , and A H can be found in other references. 9,14,26 The cross-correlation between noise sources at receivers V 1 and V 2 is where Eq. (1) is used. Inserting Eq. (2) into Eq. (4), is the crosscorrelation between the three types of waves and For uncorrelated noise, assume that hSðr s ; z s ÞS Ã ðr s 0 ; z s Þi ¼ Qðr s Þdðr s À r s 0 Þ, where h i is the ensemble average, Qðr s Þ the noise power spectrum density, and d the Dirac delta function. For uniformly distributed noise, Qðr s Þ ¼ Q is a constant. In the high-frequency regime, the spatial integration in Eq. (6) over the distribution of noise sources can be estimated using a stationary phase approximation where the phase of C ab is / ab ¼ xðt a À t b Þ. The spatial integration is estimated by finding the stationary phase points, where / ab has vanishing derivatives, d/ ab =dr s ¼ 0. Since the two receivers are on the z-axis, u s 2 ½0; 2p.
For a ¼ H or b ¼ H, the horizontal range between the noise source S and receivers V 1 and V 2 ; jr s À r 1 j ¼ jr s À r 2 j ¼ r s , should be greater than the critical offset X 6 1 ðmÞ ¼ ð2mZ À z s 6z 1 Þcot h c or X 6 2 ðnÞ ¼ ð2nZ À z s 6z 2 Þcot h c , 17 to produce head waves between them, where n 2 N þ is the number of bounces from the bottom interface at depth Z that occur between S and V 2 .

B. Evaluation of cross-correlation functions
Equation (5) shows that, by decomposing the noise field on each receiver into direct, reflected, and head waves, nine terms are obtained after the cross-correlation processing. In fact, since the reflected and head waves have both up-and down-going propagation to the receiver ( , the nine terms can be further separated into 25 terms. Specifically, the 25 terms are the result of expanding C DR ; C DH ; C RD , and C HD , where each consists of two terms (e.g., , and expanding C RR ; C RH ; C HR , and C HH , where each consists of four terms (e.g., The cross-correlation terms C DD , C DR , C RD , and C RR have been analyzed in an ocean waveguide for two arbitrarily spaced receivers in the x-z domain. 9,14 It was shown that the stationary source points at the surface are aligned with the line joining the two receivers, or joining images of two receivers. Here, a special case of two vertically spaced receivers on the z-axis is considered. The stationary sources for these four terms are therefore located right above the two receivers (r s ¼ 0). After cross-correlating noise measured on these two receivers, the direct and surface or bottom reflected waves between two receivers are obtained. In the following, only the results of these four terms are shown, the analysis will be ignored. However, the cross-correlations between head waves and direct, reflected, and head waves, C DH , C HD , C RH , C HR , and C HH , are discussed in detail in Secs. II B 2-II B 5).
Based on Eqs. (3) and (5), the cross-correlations C DD , C DR , C RD , and C RR are expressed as C DD % Qa DD e Àixðz 1 Àz 2 Þ=v 1 ; where D m ¼ m À n is the difference of seabed bounces of rays from noise source to two receivers.

C DH and C HD
Equation (6) represents the cross-correlation between direct and head waves when n¼0 a DH 6 ðn; r s Þe Ài/ DH 6 ðn;r s Þ dr s ; (8) where 2p is due to the integration of u s over 0 to 2p; a DH 6 ðn; r s Þ ¼ A D ðr s ; z 1s ÞA Ã H 6 ðn; r s ; z 2s ; c wb Þ, To find the location of stationary sources of 6 is the horizontal range of the stationary point. However, r DH 6 does not satisfy the condition for head waves as the horizontal distance between source and receiver (jr s À r 2 j ¼ r DH 6 ) is less than the critical range [X 6 2 ðnÞ], thus there are no head waves between two receivers, C DH ¼ 0. Similarly, C HD ¼ 0.

C RH
Equation (6) represents the cross-correlation between reflected and head waves when a R 6 H 6 ðm;n;r s Þe Ài/ R 6 H 6 ðm;n;r s Þ dr s ; where a R 6 H 6 ðm;n; r s Þ ¼ A R 6 ðm;r s ;z 1s ; c wb ÞA H 6 ðn;r s ;z 2s ; c wb Þ, Let d/ R 6 H 6 ðm; n; r s Þ=dr s ¼ 0, the horizontal range for the stationary point is r s ¼ r R 6 H 6 ðmÞ ¼ X 6 1 ðmÞ. Here, C R 6 H À satisfies jr s À r 2 j > X 6 2 ðnÞ when D m ! 0, while the condition is met for C R 6 H þ when D m ! 1. Therefore, C R 6 H 6 has contributions from discrete sources on concentric circles with coordinates S R 6 H 6 ðmÞ ¼ S R 6 H 6 ðr R 6 H 6 ðmÞ cos u s ; r R 6 H 6 ðmÞ sin u s ; z s Þ; u s 2 ½0; 2p. In detail, we have The location of stationary sources at the ocean surface is depicted in Fig. 2(a). For each cross-correlation term Fig. 2(a). The corresponding ray geometry for sources in panel (a) in the x-z domain are shown in panel (b). The expression of C RH is simplified as The expressions of a R À H þ ðD m Þ; a R þ H À ðD m Þ, and a R þ H þ ðD m Þ are similar to a R À H À ðD m Þ, and not shown here. From Eq. (13), by cross-correlating up-and down-going reflected and head waves, one can obtain four types of virtual head waves with different amplitudes 2pQa R 6 H 6 ðD m Þ and travel times, The virtual head wave travel times t R 6 H À ðD m ¼ 0Þ and t R 6 H þ ðD m ¼ 1Þ are equal to the travel time difference of sound waves on path I and path II (t I À t II ) in Fig. 2(b).
The virtual head waves have the same phase speed as the head waves, v 1 = sin h c , but the travel time is offset due to the cross-correlation processing, thus the term virtual. 17 From Eq. (15), the virtual head wave travel times are a function of D m . Similarly, the stationary source ranges r R 6 H 6 ðmÞ are also dependent on D m by using m ¼ D m þ n, Therefore, the virtual head wave travel times t R 6 H 6 ðD m Þ and stationary source ranges r R 6 H 6 ðD m ; nÞ are related through D m [see Fig. 3

C HR
Equation (6) represents the cross-correlation between head waves and reflected waves when a ¼ H; b ¼ R. As with C RH , the horizontal range of the stationary point is r s ¼ r H 6 R 6 ðnÞ ¼ X 6 2 ðnÞ by letting d/ H 6 R 6 ðm; n; r s Þ=dr s ¼ 0. For C H 6 R À , it satisfies jr s À r 1 j > X 6 1 ðmÞ when D m À1, while for C H 6 R þ , the condition is met when D m 0. Therefore, C H 6 R 6 has contributions from sources on concentric circles with coordinates S H 6 R 6 ðnÞ ¼ S H 6 R 6 ðr H 6 R 6 ðnÞ cos u s ; r H 6 R 6 ðnÞ sin u s ; z s Þ; u s 2 ½0; 2p. In detail, we have The distribution of stationary noise sources and the corresponding ray geometry are shown in Fig. 4. The expression of C HR is where j ¼ À1 for a H 6 R À and j ¼ 0 for a H 6 R þ . Similar to C RH , there are four types of virtual head waves with amplitudes 2pQa H 6 R 6 ðD m Þ and travel times, The virtual head wave travel times t H 6 R À ðD m ¼ À1Þ and t R 6 H þ ðD m ¼ 0Þ are equal to the travel time difference of sound waves on path II and path I (t II À t I ) in Fig. 4(b). The expressions for a H 6 R 6 , not shown here, are similar to a R À H À . The stationary source ranges r H 6 R 6 ðnÞ are dependent on D m by using n ¼ m À D m , Figure 3(b) shows the relation between virtual head wave travel times t H 6 R 6 ðD m Þ and stationary source ranges r H 6 R 6 ðD m ; mÞ. As predicted by Eq. (20), D m > À1 does not exist for C H 6 R À , and D m > 0 does not exist for C H 6 R þ . The virtual head waves produced by C H 6 R 6 at a fixed D m have contributions from multiple stationary sources [r H 6 R 6 ðD m ¼ constant; mÞ; m 2 N þ ] simultaneously. For example, for C H À R À , the virtual head waves at D m ¼ À1 (À0:05 s) have contributions from stationary sources at two discrete ranges where a H 6 H 6 ðm;n;r s Þ ¼ A H 6 ðm;r s ;z 1s ;c wb ÞA H 6 ðn;r s ;z 2s ;c wb Þ, Note, the sources in the inner diameter of the annuli, X 6 1 ðmÞ or X 6 2 ðnÞ, do not contribute to C HH . Besides, the outer diameter of each annulus is not infinite, but related to attenuation and not discussed here. The sources on the three annuli with smallest inner diameter S H À H À ðm ¼ 1Þ; S H þ H À ðm ¼ 1Þ, and S H 6 H þ ðn ¼ 1Þ are shown in Fig. 5(a), and the coordinates in the brackets are ignored. The corresponding ray geometry for panel (a) is shown in panel (b). Sources at the critical offsets X À 1 ðm ¼ 1Þ; X þ 1 ðm ¼ 1Þ; X þ 2 ðn ¼ 1Þ, and X þ 2 ðn ¼ 1Þ (hollow star) are shown in the top-to-bottom panels. The expression for C HH is where a H À H À ðD m Þ ¼ b 1 > X À 1 ðmÞ; b 2 > X À 2 ðm À D m Þ. As with C RH and C HR , there are four types of virtual head waves with amplitudes 2pQa H 6 H 6 ðD m Þ and travel times, The virtual head wave travel times t H 6 H À ðD m ¼ 0Þ are equal to t I À t II in the top two panels of Fig. 5(b), while t H 6 H þ ðD m ¼ 0Þ are equal to t II À t I in the bottom two panels. The source ranges r H 6 H 6 ðmÞ and r H 6 H 6 ðnÞ are expressed in terms of D m by using m ¼ D m þ n and n ¼ m À D m ,  Fig. 2(b), but for C HH . Sources at the critical offsets X À 1 ðm ¼ 1Þ; X þ 1 ðm ¼ 1Þ; X þ 2 ðn ¼ 1Þ, and X þ 2 ðn ¼ 1Þ (hollow star) are shown from the top to bottom panels.
a 66 ðD m Þe Àixð2D m Z6z 1 7z 2 Þ sin h c =v 1 ; From Eq. (28), each type of virtual head waves are produced by C RH , C HR , and C HH simultaneously. Therefore, the virtual head wave travel times t 66 ðD m Þ ¼ ft  Fig. 6(b). All the sources greater than or equal to the critical offsets X 6 1 ðD m þ 1Þ and X 6 2 ð1 À D m Þ [on the right hand side of Eq. (30)] have contributions to the virtual head waves. Although the virtual head waves produced by a single noise source are weak, they are observable by averaging over time of multiple noise sources.

The sum of nine terms
Based on Eqs. (7) and (28), the cross-correlation result for a vertical sensor pair is the sum of nine terms, a 66 ðD m Þe Àixð2D m Z6z 1 7z 2 Þ sin h c =v 1 : (31) The first four terms on the right hand side of Eq. (31) are from the cross-correlation of direct and reflected waves (C DD ; C DR ; C RD , and C RR ), producing physical waves including direct, surface reflected, bottom reflected, and surface-bottom reflected waves between V 1 and V 2 . The last term is from the cross-correlation between reflected and head waves (C RH , C HR ), and between head and head waves (C HH ), producing virtual head waves that do not propagate between two receivers.

III. SEABED SOUND SPEED ESTIMATION WITH VIRTUAL HEAD WAVES
It has been shown that the time domain cross-correlation function can be extracted through beamforming or seismic interferometry, and the seabed sound speed can be estimated by stacking the virtual head waves in the time domain crosscorrelation function. 15 In this section, we will show theoretically that the time domain cross-correlation function of vertical sensor pairs consists of direct, reflected, and virtual head waves based on the analysis in Sec. III. The slope of the direct and reflected waves is different from that of virtual head waves, therefore it is possible to stack the virtual head waves while other waves are added destructively.
Assuming a vertical array with N R hydrophones, V j and V k are two array elements, j; k 2 ½1; …; N R . Then, the time domain cross-correlation function c jk ðsÞ, also called virtual source gather with the virtual source at V k , 15 is obtained by transforming Eq. (31) to the time domain, while the subscripts 1 and 2 are replaced by j and k, c jk ðsÞ¼F À1 Cðr k ;z k ;r j ;z j ; x À Á a 66 ðD m Þs sÀð2D m Z6z j 7z k Þsinh c =v 1 À Á : Note, cðs; s v ; À+z jk Þ and cðs; s v ; +z jk Þ contain the same information, 17 thus only the former is summed in Eq. (38). Therefore, one obtains The virtual head waves at different D m in Eq. (39) can also be delayed and summed,

IV. SIMULATION
The model geometry is shown in Fig. 7. The full wavefield from each source (stars) to all receivers (triangles) was simulated using wavenumber integration. The full wavefield from the active sources (solid stars) and noise sources (hollow stars) was obtained from the ocean acoustic and seismic exploration synthesis (OASES) package. The noise sources have random phase and placed near the surface, approximating the noise field generated by breaking waves. Two types of sources and a vertical array form two acquisition geometries, which are, vertical array passive noise (VP) and vertical array active source (VA). It has previously been shown that both geometries are able to retrieve virtual head waves, 15 but the processing techniques are different.
According to Eqs. (31) and (32), the passive correlation c jk ðsÞ is defined as where j; k 2 ½1; …; 200. However, for the active correlation, the processing involves measuring the acoustic field at V 1 through V 200 due to the controlled, impulsive point sources near the surface. This is followed by cross-correlating over sources N S (1 N S 2401) and summing to create a virtual source at V k , 15 c A jk ðsÞ ¼ X N s l¼1 h l F À1 Pðr j ; z j ; r l ; z s ; xÞP Ã ðr k ; z k ; r l ; z s ; xÞ À Á : Here, the passive case demonstrates the theory in Sec. III that five types of waves including direct, surface reflected, bottom reflected, surface-bottom reflected waves between two receivers and virtual head waves can be extracted at the same time from noise cross-correlations. By stacking the virtual head waves, it is possible to estimate the seabed sound speed. Active sources at different ranges are simulated to validate the relation between virtual head wave travel times and location of sources that contribute to the virtual head waves. For both simulations, the frequencies are computed every 0.25 Hz from 400 to 1000 Hz. Figure 8 shows the envelope of virtual source gather c jk ðsÞ, k ¼ 1, 67, 133, and 199. That is, this plot is similar to what would be observed if the receiver (V k ) were a source and the time-domain response is plotted for all j receivers , the virtual head waves are therefore observed at In the following, active sources at different ranges are summed to validate the relation between virtual head wave travel times [t 66 ðD m Þ] and source ranges [r 66 ðD m ; mÞ and r 66 ðD m ; nÞ]. Noise sources are not considered because these are already summed in Pðr j ; z j ; xÞ and Pðr k ; z k ; xÞ before the correlation. Figure 10 shows the virtual head waves at different vertically spaced receiver pairs from sources at different ranges. The y axis of panels (a)-(d) represents sources from ranges [600: 1: SR] m (corresponding to N s ¼ SR-599 in Eq. (42), where SR ¼ [700: 50: 3000]). The reason for not using sources with ranges less than 600 m is to avoid the appearance of direct, surface reflected, bottom reflected, and surface-bottom reflected waves.
In panel (a), virtual head waves at different D m are observed, however, for short distances ([600, 1500] m), only virtual head waves at D m ¼ 0 are observed. With increasing range ([1500, 2500] m), virtual head waves at D m ¼ À1 and D m ¼ 1 begin to appear. As range increases to 2500 m, virtual head waves at D m ¼ À2 and D m ¼ 2 are weakly observed due to attenuation. Similar phenomena are observed in panels (b)-(d), but the virtual head waves shift in time and range due to different receiver depth. As explained in Sec. III, the virtual head waves are produced by three terms, C RH , C HR , and C HH , and they have contributions from all the sources greater than or equal to the critical offsets X 6 j ðD m þ nÞ and X 6 k ðm À D m Þ, see Fig. 6(b). Here, only sources at the critical offsets X 6 j ðD m þ 1Þ and X 6 k ð1 À D m Þ are plotted to simplify. The theoretical predictions in each panel match well with the simulation result.

V. DISCUSSION AND CONCLUSION
This study derives the extraction of virtual head waves from sea surface generated ambient noise in a Pekeris waveguide using a vertically spaced sensor pair. Based on ray theory, it is assumed that noise recorded at each receiver consists of direct, reflected, and head waves. After crosscorrelating noise measured on two receivers, nine terms are obtained.
Analyzing these nine terms with the method of stationary phase, it is shown that four types of physical waves (direct, surface reflected, bottom reflected, and surfacebottom reflected waves) between two receivers and virtual head waves can be extracted from ocean ambient noise. The virtual head waves are produced from cross-correlations of reflected waves and head waves, and head waves and head waves. They have contributions from sources located on an annulus, where the inner radius is the critical offset, and the external radius is the farthest source-receiver horizontal distance that makes the head waves detectable. These sources' contributions are integrated over time, making it possible to observe the virtual head waves.
The slope of the virtual head waves is different from that of physical waves in the virtual source gather, making it possible to constructively stack the virtual head waves, while the other waves are summed destructively through this processing. The virtual head waves are therefore observed at the seabed sound speed.
The simulation with noise sources confirms the theoretically predicted five types of waves and the estimation of seabed sound speed by stacking the virtual head waves. The controlled, active sources at different ranges are simulated to validate the relation between the virtual head wave travel times and ranges of sources that contribute to the virtual head waves.