Beamforming in multipath environment using the stable components of wave field
aa r X i v : . [ phy s i c s . c l a ss - ph ] J u l Beamforming in multipath environment using thestable components of wave field
A.L. Virovlyansky
Institute of Applied Physics, Russian Academy of Science, 46 Ul’yanov Street, Nizhny Novgorod, Russia,[email protected]
Abstract
The paper describes the beamforming procedures in an acoustic waveguide based onrepresenting the field on the antenna as a superposition of several stable componentsformed by narrow beams of rays [A.L. Virovlyansky, J. Acoust. Soc. Am. , 1180–1189 (2017)]. A modification of the matched field processing method is proposed, basedon the transition from comparing the measured and calculated fields on the antennato comparing their stable components. The modified approach becomes less sensitiveto the inevitable inaccuracies of the environmental model. In the case of a pulsedsource, the stable components carry signals whose arrival times can be taken as inputparameters in solving the inverse problems. The use of the stable components as theinitial fields on the aperture of the emitting antenna makes it possible to excite narrowcontinuous wave beams propagating along given ray paths. Introduction
The sound field at the antenna aperture in a multipath environment is a superposition ofseveral waves coming from different directions [1]. Each wave is formed by a beam of rayshitting the aperture which we will call an eigenbeam. If the eigenbeam is so narrow thatits width is small compared to the spatial scale of a weak sound speed perturbation, then,passing through the perturbation, all the rays forming this eigenbeam acquire approximatelythe same phase increments φ . In Ref. [2], the field component formed by such a beam iscalled stable. If the wave field is excited by a tonal source, then the stable components inthe perturbed and unperturbed waveguide differ by only a constant phase factor exp( iφ ) .In the case of transient wave field the perturbation causes only an additional time delayof the stable component as a whole. This paper considers the use of vertical antennas forreceiving and emitting sound waves in an underwater waveguide and proposes beamformingprocedures based on the use of the stable field components.An analog of conventional plane-wave beamforming in an inhomogeneous medium isthe matched field processing (MFP) [3, 4]. This method is based on comparing the vector v , whose elements are the complex amplitudes of the signals measured by the antennaelements, with the vector u representing the theoretical estimate of v calculated using theavailable environmental model. When solving the problem of source localization [5,6] and/orreconstruction of environmental parameters [7], the measured vector v is compared withthe vectors u calculated for different source positions and/or different values of unknownenvironmental parameters. The desired estimates are given by the values of the sourcecoordinates and medium parameters corresponding to the maximum of the scalar product of u and v . However, this approach is effective only in the case of a fairly accurate environmentalmodel. Otherwise, that is, under conditions of uncertain environment, it may turn out thatthe maximum of the scalar product corresponds to values of unknown parameters that differsignificantly from the true ones.A number of approaches have been developed for working in uncertain environments.One of them is based on reducing the sensitivity of MFP through the use of a multiplyconstrained beamformer [8]. This suggests that opening up the search window in one ormore of these parameters would make the beamformer more tolerant of uncertainty in theother parameters [6, 9, 10]. Another well-known approach is based on solving the problem oflocalization by incorporating environmental variability a priori [11–13]. The current state ofresearch related to the use of MFP in ocean acoustics is presented in the review [14].This paper discusses an alternative approach, which we call the generalized matched field2rocessing (GMFP). Its idea is to move from comparing vectors v and u to comparing theirstable components. This relaxes the accuracy requirements for the environment model. Theuse of GMFP is illustrated by a numerical example.In the case of a pulsed source, a procedure has been proposed for isolating signals carriedby stable components, that is, arriving at the antenna through individual eigenbeams. Thearrival times of these signals can be used as input parameters in solving inverse problems.It is also shown that the use of stable components in combination with the phase con-jugation method [15, 16], allows one to create the field distributions on the aperture of theradiating antenna for emitting narrow wave beams, propagating along given ray paths.The paper is organized as follows. The representation of the field at the aperture ofthe receiving antenna as a superposition of stable components is introduced in Sec. 2. Amethod for isolating the stable components from the total field is outlined in Sec. 3. Section4 discusses a generalization of the MFP method. Section 5 describes a procedure of isolatingthe stable components from the sound field excited by a broadband source. The use of avertical emitting antenna for exciting a narrow wave beam propagating along a given raypath is discussed in Sec. 6. The results of the work are summarized in Sec. 7. To illustrate the general statements discussed in this and subsequent sections, we will use anidealized 2D model of an underwater acoustic waveguide in the deep sea with coordinates ( r, z ) , where r is the distance and z is the depth. The z-axis is directed downward, the watersurface is in the plane z = 0 , and the bottom in the plane z = h .The vertical antenna is located along the straight line r = 0 and covers the depth interval z up ≤ z ≤ z do . We will consider the sound speed field in the form c ( r, z ) = c b ( z ) + δc ( r, z ) , where c b ( z ) is the background sound speed profile and δc ( r, z ) is the sound speed fluctuation.As an example, in this work, we use a waveguide of depth h = 3 km with the profile c b ( z ) shown in the left panel of Fig. 1. The refractive index is ν ( r, z ) = c /c ( r, z ) , where c = 1.5km/s is the reference sound speed.In what follows we consider the sound field excited by the source set at the point ( r , z ) ,where r = 30 km and z = 0.7 km. The field is recorded by a vertical antenna with a lengthof L = 0.25 km. The coordinates of its end points are z up = 0.75 km and z do = 1 km. It isassumed that the bottom is strongly absorbing and therefore the contributions of the bottom3 .46 1.475 c (km/s) z ( k m ) r (km)1 2 3antenna (r ,z ) Figure 1: Left panel. Background sound speed profile c b ( z ) . Right panel. Beams of ray(eigenbeams) hitting the antenna aperture. Next to each eigenbeam, its number is indicated.reflected waves are negligible. We will consider the CW fields at the carrier frequency f =500 Hz. When simulating pulsed signals, this will be the center frequency.The launch angles χ of rays propagating without reflections from the absorbing bottomsatisfy the condition | χ | < χ max , where χ max is the critical angle [5]. It is assumed thatthere are N eigenbeams arriving at the antenna. The launch angles of the rays forming the n -th eigenbeam fill the subinterval χ ′ ,n < χ < χ ′′ ,n of the angular interval ( − χ max , χ max ) .The subintervals corresponding to different n do not overlap. If the antenna length tendsto zero, then eigenbeams turn into ordinary eigenrays, that is, into rays arriving at a givenpoint.In our example, χ max = 12 . ◦ and there are three eigenbeams ( N = 3) shown in the rightpanel of Fig. 1. Launch angles of eigenbeams 1, 2, and 3 lie in the intervals ( − . ◦ , − ◦ ) , (1 . ◦ , . ◦ ) , and (9 . ◦ , . ◦ ) , respectively.The rays belonging to the n -th eigenbeam in the unperturbed waveguide form the com-ponent of the total field, which we denote by u n ( z ) . The total field on the antenna is asuperposition of these components u ( z ) = N X n =1 u n ( z ) . (1)If the eigenbeam is narrow enough, then it forms a component, which in Ref. [2] is called4table. Let us dwell on this issue.In the geometric optics approximation, the contribution to the total field from the raywith the launch angle χ is A ( χ ) exp [ ikS ( χ )] , where A ( χ ) and S ( χ ) are the amplitudeand eikonal of the ray at the observation range, respectively, and k = 2 πf /c is the referencewavenumber [5, 17]. On short paths, a variation of the ray trajectory and amplitude in thepresence of a weak perturbation δc can be neglected [17, 18]. In this case, the influence ofthe perturbation is taken into account by replacing the eikonal S ( χ ) by S ( χ ) + δS ( χ ) ,where δS ( χ ) = − c Z Γ δc ds (2)is the eikonal increment, Γ is the unperturbed ray path, ds is the arc length. In a deepocean at frequencies of order 100 Hz this approximation is applicable at ranges up to a fewhundred kilometers [18].If the vertical spread of the n -th eigenbeam does not exceed the vertical scale of the per-turbation δc , all rays from this eigenbeam intersect approximately the same inhomogeneitiesand acquire approximately the same eikonal increments δS ( χ ) . Then the phase incrementsof these rays, φ ( χ ) = kδS ( χ ) , are close to the same value, which we denote φ n . If thiscondition is satisfied, we call the component u n ( z ) stable. The influence of the perturbationis manifested only in the multiplication of u n ( z ) by the phase factor exp ( iφ n ) independentof z . The values of | φ n | are not necessarily small. They can significantly exceed π .If the eigenbeam is not narrow enough, it can be divided into several narrower eigenbeams.In the limiting case, when the antenna covers the entire cross-section of the waveguide,the entire interval of launch angles ( − χ max , χ max ) can be divided into small subintervals sothat each of them corresponds to a stable component. This situation was considered inRefs. [19, 20].If all components u n ( z ) are stable, then in the presence of the perturbation, the field u ( z ) at the antenna becomes v ( z ) = N X n =1 γ n u n ( z ) , (3)where γ n = exp( iφ n ) . Since different eigenbeams intersect different inhomogeneities, therandom values of φ n for different n are independent.As a quantitative characteristic of the spread of φ n corresponding to different rays forming5he n -th eigenbeam, we take ∆ φ n = k χ n Z χ ′′ ,n χ ′ ,n (cid:10) q ( χ ) (cid:11) dχ ! / , (4)where q ( χ ) = δS ( χ ) − δS ( ¯ χ n ) , ¯ χ n = ( χ ′′ ,n + χ ′ ,n ) / is the launch angle of the central ray, ∆ χ n = χ ′′ ,n − χ ′ ,n , the angular brackets denote averaging over the random realizations of δc .Strictly speaking, Eq. (3) is valid if ∆ φ n ≪ π. (5)In this paper, we model the perturbation δc ( r, z ) by a Gaussian random function withzero mean ( h δc i = 0 ) and the correlation function h δc ( r, z ) δc ( r ′ z ′ ) i = δc rms × exp − π ( r − r ′ ) l r − π ( z − z ′ ) l z ! , where δc rms is the rms amplitude of the sound speed fluctuations, l r and l z are the horizontaland vertical correlation scales, respectively. This simplest model significantly differs frommore realistic models used in ocean acoustics to describe the sound speed fluctuations in thedeep sea [18, 21]. Nevertheless, as shown in Ref. [19], it is suitable for the general analysis ofstable components. In what follows we take δc rms = 0.25 m/s, l r = 5 km, l z = 0.5 km.Using this model for the 1st, 2nd, and 3rd eigenbeams shown in the right panel of Fig. 1,we find ∆ φ = 0.73, ∆ φ = 1.83, and ∆ φ = 0.69, respectively. Thus, in our example, thecondition (5) is not met. However, we assume that if the weaker condition ∆ φ n < π (6)is satisfied, the field v ( z ) can be approximately represented as a linear combination of stablecomponents u n ( z ) with weight factors γ n different from exp ( iφ n ) .Representation of the field on the antenna as a superposition of stable components un-derlies the beamforming procedures discussed below. In order to use this representationone should calculate the functions u n ( z ) . In free space, where there is only one eigenbeam( N = 1 ) the solution is obvious. If the antenna is located far from the source, the eigenbeamis formed by almost parallel straight lines and u ( z ) is a fragment of a plane wave. However,in a refractive medium, the situation becomes more complicated: a beam of rays describes awave whose phase front curvature depends on the range and, generally, varies from zero toinfinity. In such a medium, there are caustics in the vicinity of which the ray approximationfails [5, 17, 22]. 6ince the stable components are introduced using the ray-based representation of thewave field, we have no rigorous definition of u n ( z ) . Two heuristic methods for isolatingstable components from the total field u ( z ) are proposed in Ref. [2]. Both methods comefrom the relationship between the ray-based and wave-based descriptions of sound fields andgive similar results. In this paper, we will use one of these methods outlined in the nextsection. This section presents the method of isolating the stable components using the coherent stateexpansion developed in quantum mechanics [23–25].
To describe the ray trajectories, we apply the Hamiltonian formalism [26, 27]. In the scopeof this formalism, the trajectory at a range r is defined by its depth z and momentum p = ν ( r, z ) sin χ , where χ is the grazing angle at the point ( r, z ) . The functions p ( r, p , z ) and z ( r, p , z ) , representing the solutions of the Hamilton ray equations with the initialconditions p = p and z = z at r = r , determine the ray path in the phase space ( r, p, z ) .In the case of a source located at a point ( r , z ) , all the rays start from this pointwith different launch angles χ and, accordingly, with different initial momenta and p = ν ( r , z ) sin χ . In the phase plane (momentum P , depth Z ) at the observation distance r ,the arrival of one ray is represented by a point. These points form a curve representing aLagrange manifold [22]. We will call this curve a geometric ray line or simply a ray line. Itis defined parametrically by the equations P = p ( r, p , z ) and Z = z ( r, p , z ) with fixed r and z .In Fig. 2, the solid line shows the ray line at a distance r = 0 in our example introducedin Sec. 2. The horizontal dashed lines indicate the depths of the antenna endpoints. Thesestraight lines “cut out” the segments of the ray line shown by bold curves. Each segmentrepresents the arrivals of the rays forming one of the three eigenbeams shown in Fig. 1. Thenumbers of eigenbeams are indicated next to the corresponding segments.To isolate the contributions of waves arriving in a neighborhood of the depth Z undergrazing angles close to χ = arcsin( P/ν ( r, Z )) , following [2, 28], we use the coherent state7 P Z ( k m ) Figure 2: Geometric ray line (solid curve) and fuzzy ray line (gray area) at the observationdistance r = 0 km. Dashed lines show the horizons of the antenna endpoints. Bold segmentsof the geometric ray line represent ray arrivals that hit the antenna and form the eigenbeams.Fuzzy segments corresponding to individual eigenbeams are highlighted in dark gray.8xpansion. The coherent state associated with the point µ = ( P, Z ) of the phase plane isgiven by the function [24, 25, 29] Y µ ( z ) = 1 √ ∆ z exp " ikP ( z − Z ) − π ( z − Z ) z , (7)where ∆ z is the vertical scale.Although the coherent states are not orthogonal, they form a complete system of functionsand an arbitrary function u ( z ) can be represented as an expansion [24, 25] u ( z ) = λ − Z dµ a µ Y µ ( z ) , (8)where λ = 2 π/k is the wavelength, a µ = Z dz u ( z ) Y ∗ µ ( z ) (9)and superscript asterisk denotes complex conjugation. The integration with respect to µ formally goes over the entire phase plane and the integration with respect to z goes over theentire vertical axis. However, from (7) it is clear that the main contribution to the integral(9) comes from the interval Z ± ∆ z .The closeness of the coherent states associated with the points of the phase plane µ =( P, Z ) and µ = ( P , Z ) , can be quantitatively characterized by their squared scalar product (cid:12)(cid:12)(cid:12)(cid:12)Z dz Y µ ( z ) Y ∗ µ ( z ) (cid:12)(cid:12)(cid:12)(cid:12) = e − d ( µ,µ ) , (10)where d ( µ, µ ) = π ( P − P ) ∆ p + π ( Z − Z ) ∆ z , (11) ∆ p = λ/ (2∆ z ) . We will interpret the quantity d ( µ, µ ) as a dimensionless distance betweenthe points µ and µ . Coherent states are close if this distance is small compared with unity.The distance from the point µ to a curve in the phase plane (for example, to the ray line orto its segment) is the distance from µ to the nearest point of the curve.In quantum mechanics, the wave function Y µ ( z ) defines a state with a minimum uncer-tainty [24, 25, 29]. In acoustics Y µ ( z ) can be interpreted as the cross section of a Gaussianwave beam arriving at the grazing angle χ = arcsin ( P/ν ) in a neighborhood of the depth Z . According to (9) – (11), the coherent state amplitude, a µ , is formed by contributions ofwaves arriving in the depth interval Z ± ∆ z / under grazing angles corresponding to themomentum interval P ± ∆ p / . 9quations (9) – (11) suggest that the distribution of the squared coherent state amplitude | a µ | in the phase plane is localized mainly in the region, formed by points located at distances d < (12)from the ray line. We will call this region a fuzzy ray line and denote it by the symbol σ .In Fig. 2, the area occupied by the fuzzy ray line is highlighted in gray. Here the coherentstate expansion is performed with ∆ z = 90 m. With this choice of ∆ z , the area of the fuzzyray line, and hence its average width, in our example takes the minimum value. Consider a segment of the geometrical ray line presenting the arrivals of the rays formingthe n -th eigenbeam. It is natural to assume that the contribution of these rays to the totalfield is represented by a superposition of coherent states associated with points of the P − Z plane located at distances d < from the segment. We will call the area occupied by thesepoints a fuzzy segment and denote it by σ n . In Fig. 2, the fuzzy segments are highlighted indark gray.According to this assumption, the contribution of the n -th eigenbeam to the total field is u n ( z ) = λ − Z σ n dµ a µ Y µ ( z ) . (13)Using Eq. (9), this expression can be rewritten as u n ( z ) = Z dz ′ Ξ n ( z, z ′ ) u ( z ′ ) , (14)where Ξ n ( z, z ′ ) = λ − Z σ n dµ Y µ ( z ) Y ∗ µ ( z ′ ) . (15)Equation (14) explicitly defines the procedure for isolating the desired component u n ( z ) fromthe total field u ( z ) calculated or measured using a sufficiently long antenna. The components u n ( z ) found in this way are determined in the entire cross-section of the waveguide.Note that locally horizontal portions of the ray line correspond to caustics [22]. In Fig. 2it is seen that such a portion is present on the bold segment formed by the rays forming the2nd eigenbeam. In the right panel of Fig. 1 we see that the antenna crosses a caustic formedby rays from this eigenbeam. Therefore, the field component formed by this eigenbeamcannot be described in the geometrical optics approximation. However, for the applicabilityof the method outlined here, it is not an obstacle. The total field u ( z ) at the observation10 e ( u ) R e ( u ) z (km) R e ( u ) bac Figure 3: The real parts of the functions u ( z ) (a), u ( z ) (b), and u ( z ) (c) on the antennaaperture.range present in Eq. (14) should be obtained by a full-wave calculation. In this paper, weapply the method of wide-angle parabolic equation [5]. The ray tracing is used only to findthe fuzzy segment σ n .Figure 3 shows the real parts of the field components u ( z ) , u ( z ) , and u ( z ) at theantenna aperture considered in our example. The classical MFP method [3–5] is based on the assumption that the available environmentalmodel is accurate enough to correctly calculate the field on the antenna for the known sourceposition. In this section, a generalized version of this method is introduced, based on theweaker assumption that the field on the antenna can be represented as Eq. (3), that is, asthe sum of the known stable component u n ( z ) with unknown weights γ n .11 .1 Projection of the measured field onto stable components We assume that the antenna is an array consisting of a large number of elements denselyfilling the aperture (e.g., ten elements per wavelength). As in the Introduction, we willrepresent the calculated and measured fields by the vectors u and v , respectively, whoseelements are the complex amplitudes of the signals at the antenna elements. The proximityof vectors u and v is quantitatively characterized by their normalized squared scalar product K = | u + v | | u | | v | = v + P v | v | , (16)where P = uu + | u | (17)is the projection matrix. The coefficient (16) is proportional to the square of the projectionof the measured vector v onto the calculated vector u .If u is considered as a function of unknown parameters θ that specify the coordinatesof the source and/or parameters of the medium, then the similarity coefficient K is alsoa function of θ . The maximum of this function corresponds to the actual source positionand/or the correct value of the medium parameters. This can be used in solving inverseproblems. However, as indicated above, this is true only when the model of the medium issufficiently accurate.Our idea is to relax the requirements for the accuracy of the environmental model byusing the field representation (3). In matrix notation, this expression takes the form v = N X n =1 γ n u n , (18)where u n are vectors representing the stable field components at the antenna aperture. Sucha representation of the measured vector v adequately describes the situation arising frommultipath sound propagation.Note that u and v are elements of the vector space Ω formed by vectors of size N a × ,where N a is the number of antenna elements. The subspace of Ω , whose basis is given by thevectors u n , n = 1 , . . . , N , we denote by Ω ′ . Let us introduce the matrix W = [ u , . . . , u N ] ,whose columns are the vectors u n , and use its singular value decomposition [30] W = N X n =1 α n ξξξ n ηηη + n , α n are the singular numbers, ξξξ n and ηηη n are the singular vectors. The projection matrix P = N X n =1 ξξξ n ξξξ + n (19)projects any vector from Ω onto Ω ′ . The new similarity coefficient, characterizing the prox-imity of u and v , we define as K = v + Pvv + v . (20)It differs from (16) by replacing the projection matrix P by P .If the vectors u and v correspond to the same source position and v can be representedin the form (18), then Pv = v and K = 1 . The value of K averaged over the ensemble ofsound speed fluctuations can be easily estimated under the assumption that for n = m thevectors u n and u m are almost orthogonal, that is, (cid:12)(cid:12) u + n u m (cid:12)(cid:12) ≪ | u n | | u m | . (21)This situation is typical of an antenna whose length L is large compared to the wavelength λ . If the vectors u n represent the arrivals of quasi-plane waves, then (21) is satisfied underthe condition λ/L ≪ χ max . We also assume that the phase increments φ n , n = 1 , . . . , N , areindependent random variables uniformly distributed in the interval (0 , π ) . Then h K i = P Nn =1 | u n | (cid:16)P Nn =1 | u n | (cid:17) . It follows that h K i is always in the range from /N (if all | u n | are equal) to 1. For individualrealizations of perturbation, the value of K can be less than /N .Note that if the condition (21) is satisfied and ξξξ n ≃ u n / | u n | , then v + Pv ≃ N X n =1 (cid:12)(cid:12) u + n v n (cid:12)(cid:12) , where v n = γ n u n . This relation allows us to interpret the transition from K to K , asthe transition from comparing the measured and calculated fields to comparing their stablecomponents.Let us turn to the example described in Sec. 2. Figure 4 shows the results of calculatingthe coefficients K (circles) and K (asterisks) for 40 realizations of the perturbation δc . Thecalculations were performed for the positions of the source and antenna shown in Fig. 1. Thefact that the coefficient K in Fig. 4 is less than unity indicates that the field representation13
10 20 30 40
Realization S i m il a r i t y c oe ff i c i en t K K Figure 4: Similarity coefficients K and K calculated for 40 realizations of random per-turbation δc ( r, z ) . Both coefficients are found for the perturbed and unperturbed fieldscorresponding to the same source position ( r , z ) .(18) is not entirely accurate and the vector v is only approximately described by a superpo-sition of stable components. Nevertheless, for each realization of δc , the value of K exceeds K . This is consistent with our expectation that the coefficient K is less sensitive to soundspeed fluctuations than K .The similarity coefficient K can be used to solve the same inverse problems as K . TheMFP method, modified by replacing K with K , we will call generalized matched fieldprocessing (GMFP). In the absence of multipath, N = 1 and GMFP reduces to MFP. It isnatural to expect that the GMFP method is more robust and less sensitive to inaccuraciesin the environmental model. In the next section, this will be demonstrated by a numericalexample. Consider the use of MFT and GMFT to estimate the source coordinates in our waveguidemodel. It is assumed that the antenna receives the signals of the tonal source placed atthe point ( r , z ) indicated in Fig. 1. The vector v , calculated for a realization of δc ,simulates the measured field. It is compared with the fields on the antenna in the unperturbedwaveguide ( δc = 0 ) excited by sources located at different points ( r s , z s ) . The vectors ofsignal amplitudes u , stable components u n , and matrices P and P were computed for a setof source positions ( r s , z s ) . Then, using Eqs. (16) and (20), similarity coefficients K and K were found for each position. The arguments of the uncertainty functions K ( r s , z s ) and K ( r s , z s ) , corresponding to the main maxima of these functions, are taken as the estimates14 .40.60.81 z s ( k m )
24 26 28 30 32 34 36 r s (km) z s ( k m ) K , KK , MFPK, GMFP Figure 5: The uncertainty functions K ( r s , z s ) (upper panel) and K ( r s , z s ) (lower panel) inthe unperturbed waveguide ( δc = 0 ).of the actual source coordinates.Figure 5 shows the uncertainty functions in the unperturbed waveguide ( δc = 0 ). As itshould be, both functions take maximum values at the point with coordinates z s = z = 0 . km and r s = r = 30 km. The maximum of the function K ( r s , z s ) is so sharp that it ishardly distinguishable in the figure. The function K ( r s , z s ) has a much wider maximum.This indicates a weak sensitivity of the coefficient K not only to the sound speed variationsbut also to the variations in the source position.Similar uncertainty functions calculated for δc = 0 , are shown in Fig. 6. In the presenceof fluctuations, the main maximum of K ( r s , z s ) splits into many small local maxima. Thefunction K ( r s , z s ) remains smooth and takes the largest values (exceeding the values of K ( r s , z s ) ) in the vicinity of the point ( r , z ) .The uncertainty functions K ( r s , z s ) and K ( r s , z s ) were calculated for 40 realizations ofthe perturbation δc . The functions averaged over all the realizations are shown in Fig. 7.Both averaged functions have wide maxima centered approximately at the point ( z , r ) .However, the values of h K i in the vicinity of this point are two times higher than the valuesof h K i . This is consistent with the data presented in Fig. 4, which shows the values of K and K at the point ( r , z ) for the same 40 realizations of δc . In this paper, we do not15 .40.60.81 z s ( k m )
24 26 28 30 32 34 36 r s (km) z s ( k m ) K , KK, GMFPK , MFP Figure 6: The same as in Fig. 5, but in the presence of a realization of δc ( r, z ) . z s ( k m )
24 26 28 30 32 34 36 r s (km) z s ( k m ) K , KK , MFPK, GMFP Figure 7: The same as in Figs. 5 and 6, but the uncertainty functions are averaged over 40realizations of δc . 16ake into account the external noise. However, the results presented in Fig. 7 suggest thatthe transition from MFP to GMFP can enhance the signal-to-noise ratio at the beamformeroutput. Let us proceed to the analysis of the field excited by a pulsed point source and consider theisolation of sound pulses arriving at the antenna via individual eigenbeams. In this case, thefields on the antenna in a perturbed and unperturbed waveguide can be represented as V ( z, t ) = Z df v ( z, f ) s ( f ) e − πift (22)and U ( z, t ) = Z df u ( z, f ) s ( f ) e − πift , (23)respectively, where t is the time, s ( f ) is the spectrum of the emitted pulse. Here we explicitlyindicate the argument f of functions v ( z, f ) and u ( z, f ) , which is omitted in the othersections of this paper. Since the ray trajectories are frequency independent, we will assumethat the eigenbeams are the same at all frequencies.In principle, the conditions (5) or (6) may not be satisfied for all the frequencies from thesource bandwidth. Therefore, at some frequencies, it may be necessary to split eigenbeamsinto narrower ones, as indicated in Sec. 2. We do not consider such a situation, assumingthat the bandwidth is not very large.Functions u n ( z, f ) representing the contributions of individual eigenbeams can be foundat any frequency f using Eq. (14). To isolate the contribution of the n -th eigenbeam tothe total pulsed field, we will proceed as follows. Assuming that functions u n ( z, f ) withdifferent n are almost orthogonal, that is, if the condition (21) is satisfied at each frequency,we calculate the projections of the field v ( z, f ) on the function u n ( z ) g n ( f ) = R z do z up dz u ∗ n ( z, f ) v ( z, f ) (cid:16)R z do z up dz | u n ( z, f ) | (cid:17) / . (24)We define the contribution of the n -th eigenbeam, as V n ( z, t ) = Z df g n ( f ) u n ( z, f ) e − πift . (25)17 .30.60.9 z ( k m ) z ( k m ) z ( k m ) t (s) z ( k m ) |V(z,t)|ab |V (z,t)|c |V (z,t)|d |V (z,t)| Figure 8: Amplitudes of the received signals in the plane ’arrival time – depth’, ( t, z ) . Thesolid line, the same on all panels, represents the timefront, that is, the distribution of rayarrivals in the ( t, z ) plane. Dashed lines show the horizons of the antenna endpoints. (a) Thetotal field at the observation distance r = 0 in the unperturbed waveguide. (b,c,d) Pulsedfield components formed by the 1-st, 2-nd, and 3-rd eigenbeams, respectively, in the presenceof a realization of δc ( r, z ) . 18he results of applying this procedure in our example are presented in Fig. 8. Thecalculations were performed for an emitted signal with the spectrum s ( f ) = exp (cid:2) − πτ ( f − f ) (cid:3) , where f = 500 Hz, τ = 0.01 s. Figure 8a shows the distribution of the total field amplitudein the unperturbed waveguide ( δc = 0 ) in the plane "time - depth" ( t, z ) . Figures 8b, c,and d present similar distributions of the field components formed by the 1st, 2nd, and3rd eigenbeam, respectively, in the presence of perturbation δc . The solid line, the sameon all panels, shows the so-called timefront, that is, the distribution of ray arrivals in the ( t, z ) plane at the observation distance r = 0. Dashed lines show the depths of the antennaendpoints. The timefront segments between the dashed lines depict the arrivals of the raysbelonging to the eigenbeams. Thus, we see that the described procedure allows one to isolatethe field components coming through different eigenbeams.The pulse G n ( t ) = Z df g n ( f ) e − πift can be interpreted as the total signal arriving at the antenna via the n -th eigenbeam. Notethe interesting fact that in the unperturbed waveguide, the arrival times of all pulses G n ( t ) are the same. Indeed, if v ( z, f ) = u ( z, f ) and the condition (21) is satisfied at all frequencies,then g n ( f ) ≃ Z z do z up dz | u n ( z, f ) | ! / and all G n ( t ) take maximum values at the same moment t = 0. If the fields U ( z, t ) and V ( z, t ) are excited at different times t u and t v , respectively, then all pulses G n ( t ) will takemaximum values at time t = t v − t u .The pulses G n ( t ) calculated in the unperturbed (thin lines) and perturbed (thick lines)waveguide are shown in Fig. 9. Figure 9a, b, and c show the arrivals of pulses representingthe contributions of the 1st, 2nd, and 3rd eigenbeams, respectively. In accordance withthe above, in an unperturbed waveguide, all pulses G n ( t ) arrive at the same time. In thepresence of the sound speed perturbation δc , the arrival time of the n -th pulse changes by δt n = δS n /c , where δS n is the eikonal increment, which is approximately the same for all raysbelonging to the eigenbeam. These delays are approximately equal to those delays in the raytravel times, which are used as input parameters in the ocean acoustic tomography [31, 32].19 | G ( t ) | | G ( t ) | -20 -10 0 10 20 t (ms) | G ( t ) | abc Figure 9: Sound pulses arriving at the antenna via the 1-st (a), 2-nd (b), and 3-rd (c)eigenbeams in an unperturbed waveguide (thin solid lines) and in the presence of a realizationof δc ( r, z ) (thick solid lines). 20
10 20 30 r (km) z ( k m ) Figure 10: Sound field focused at the point ( r , z ) using the phase conjugation. Dashed linesshow the central rays of the eigenbeams. Next to each line, the number of the correspondingeigenbeam is indicated. So far, we have been considering a receiving antenna. Now we turn to the emitting an-tenna and show how to use a stable component for exciting a narrow continuous wave beampropagating along a ray path.We will apply the well-known method of phase conjugation, which is used to focus thefield at a given point [15, 16]. Let us denote by u ( z ) the complex amplitude of the signalemitted by the antenna element at a depth z . In order to focus the field at a given point ( r , z ) , one should take u ( z ) = u ∗ ( z ) , where u ( z ) is the field at the antenna aperturereceived from the source placed at ( r , z ) .Consider an example of such focusing. The antenna shown in the right panel of Fig.1, now will be considered as the emitting one, and the source position ( r , z ) will be con-sidered as the focus point. Figure 10 presents the sound field excited in the unperturbedwaveguide by this antenna with the initial field u ( z ) = u ∗ ( z ) B ( z ) , where B ( z ) is theweight factor that ensures smooth decay of the sound field to the antenna endpoints. Weuse B ( z ) = exp (cid:0) − π ( z − z c ) /L (cid:1) , where z c = ( z do + z up ) / is the depth of the antennacenter. The excited field represents a superposition of three wave beams propagating alongthe eigenbeams. The dashed lines show the trajectories of the central rays of the eigenbeamsfrom Fig. 1.Our idea is to use as the initial fields the components u ∗ n ( z ) corresponding to individual21 z ( k m ) r (km) a bc Figure 11: The wave beams propagating along the paths shown by the dashed lines. Thepaths represent the central rays of the 1-st (a), 2-nd (b), and 3-rd (c) eigenbeams.eigenbeams. The antenna with an initial field u ∗ n ( z ) should emit a wave beam propagatingalong the n -th eigenbeam. Figures 11a, b, and c show the wave fields excited by our antennawith u ( z ) equal to u ∗ ( z ) B ( z ) , u ∗ ( z ) B ( z ) , and u ∗ ( z ) B ( z ) , respectively, where u ( z ) , u ( z ) , and u ( z ) are the field components whose real parts are shown in Fig. 3. Thesebeams are calculated in the unperturbed waveguide. In the presence of the sound speedfluctuations, they change insignificantly (not shown).The simulation results presented in Fig. 11 show that the proposed method makes itpossible to excite narrow wave beams associated with individual eigenbeams, and therebyuse the antenna in the inhomogeneous medium as an acoustic searchlight. The problem of beamforming in an inhomogeneous environment in this paper is addressedusing the notion of the stable components of the sound field introduced in Ref. [2]. Theassumptions underlying this notion do not yet have a rigorous justification and are basedonly on simple estimates derived in the geometrical optics approximation. Nevertheless, the22umerical simulation confirms that the components, called stable, are indeed less sensitiveto the sound speed variations than the total wave field [19, 20, 28].In Sec. 4 it is shown that the use of the field representation as a superposition of thestable components makes it possible to modify the traditional MFP method and make it lesssensitive to the inevitable inaccuracies of the environmental model. This approach can beapplied to develop robust methods of source localization and reconstruction of environmentalparameters. In Refs. [28] and [20], the stable components have already been used for solvingsimilar problems. However, the approaches considered in these works imply the coherentstate expansion of not only the calculated fields but the measured field as well. This requiresa long antenna, whose size significantly exceeds the scale of the coherent state ∆ z . Incontrast, in the present paper, only the calculated fields are decomposed into the coherentstates. The vectors u n are determined by fragments of the calculated stable componentswithin the depth interval overlapped by the antenna. However, in this case, there are stillrestrictions on the antenna length. On a very short antenna, the vectors u n with different n become indistinguishable. To fulfill the condition (21), the antenna aperture should be largeenough.The isolation of stable components from the total field excited by a pulsed source isconsidered in Sec. 5. Here, the pulses G n ( t ) are defined, which are interpreted as signalsarriving at the antenna through individual eigenbeams. In the presence of perturbation,their arrival times acquire different delays that carry information about the inhomogeneitiesthrough which the eigenbeams pass.The main attention in this paper is given to the beamforming for a receiving antenna.But in Sec. 6 it is demonstrated numerically that the use of the stable components to formthe amplitude-phase distributions on the antenna elements allows one to emit narrow wavebeams propagating along given ray paths.The author is grateful to Dr. L. Ya. Lyubavin and Dr. A. Yu. Kazarova for valuablediscussions. The research was carried out within the state assignment of IAP RAS (Project0035-2019-0019). References [1] H. L. V. Trees,
Optimum Array Processing: Detection Estimation and Modulation The-ory Part IV (NY:John Wiley and Sons, New York, 2002).232] A. Virovlyansky, “Stable components of sound fields in the ocean,” J. Acoust. Soc. Am. (2), 1180–1189 (2017).[3] M. Hinich, “Maximum-likelihood signal processing for a vertical array,” J. Acoust. Soc.Am. (2), 499–453 (1973).[4] H. Bucker, “Use of calculated sound fields and matched-field detection to locate soundsources in shallow water,” J. Acoust. Soc. Am. (2), 368–373 (1976).[5] F. Jensen, W. Kuperman, M. Porter, and H. Schmidt, Computational Ocean Acoustics (Springer, New York, 2011).[6] A. Baggeroer, W. Kuperman, and P. Mikhalevsky, “An overview of matched field meth-ods in ocean acoustics,” J. Ocean. Eng. (4), 401–424 (1993).[7] A. Tolstoy, “Review of matched field processing for environmental inverse problems,”International journal of modern physics C (4), 691–708 (1992).[8] H. Schmidt, A. Baggeroer, W. Kuperman, and E. Scheer, “Environmentally tolerantbeamforming for high-resolution matched field processing: Deterministic mismatch,” J.Acoust. Soc. Am. (4), 1851–1862 (1990).[9] J. Krolik, “Matched-field minimum variance beamforming in a random ocean channel,”J. Acoust. Soc. Am. (3), 1408–1419 (1992).[10] G. Byun, F. Akins, K. Gemba, H. Song, and W. A. Kuperman, “Multiple constraintmatched field processing tolerant to array tilt mismatch,” J. Acoust. Soc. Am. (2),1231–1238 (2020).[11] A. Richardson and L. Nolte, “A posteriori probability source localization in an uncertainsound speed, deep ocean environment,” J. Acoust. Soc. Am. (5), 2280–2284 (1991).[12] S. Dosso and M. Wilmut, “Comparison of focalization and marginalization for bayesiantracking in an uncertain ocean environment,” J. Acoust. Soc. Am. (2), 717–722(2008).[13] S. Dosso and J. Dettmer, “Bayesian matched-field geoacoustic inversion,” Inverse Prob-lems , 055009 (23pp) (2011).[14] A. Sazontov and A. Malekhanov, “Matched field signal processing in underwater soundchannels (review),” Acoust. Phys. (2), 213–230 (2015).2415] D. Jackson and D. Dowling, “Phase-conjugation in underwater acoustics,” J. Acoust.Soc. Am. , 171–181 (1991).[16] M. Fink and C. Prada, “Acoustic time-reversal mirrors,” Inverse Problems , R1–R38(2001).[17] L. Brekhovskikh and Y. Lysanov, Fundamentals of Ocean Acoustics (Springer-Verlag,New York, 2003).[18] S. Flatte, R. Dashen, W. Munk, K. Watson, and F. Zakhariasen,
Sound transmissionthrough a fluctuating ocean (Cambringe U.P., London, 1979).[19] A. Virovlyansky and Y. Makarova, “On spatial structure of the wave field in a verticalsection of a deep water acoustic waveguide,” EPL , 54004 (2018).[20] A. Virovlyansky, A. Kazarova, and L. Lyubavin, “Matched field processing in phasespace,” J. Ocean Eng. (Early Access) (2019) doi: 10.1109/JOE.2019.292765.[21] J. Colosi,
Sound propagation through the stochastic ocean (Cambridge University Press,New York, 2016).[22] M. Alonso, “Rays and waves,” in
Phase-space optics. , edited by M. Testorf, B. Hennely,and J.Ojeda-Castaneda (McGraw-Hill, New York, 2010), Chap. 8, pp. 237–277.[23] R. Glauber,
Quantum Theory of Optical Coherence. Selected Papers and Lectures (Wiley-VCH, Weinheim, 2007).[24] J. Klauder and E. Sudarshan,
Fundamentals of quantum optics (W.A. Benjamin, NewYork, 1968).[25] W. Schleich,
Quantum Optics in Phase Space (Wiley-VCH, Berlin, 2001).[26] A. Virovlyansky, A. Kazarova, and L. Lyubavin, “Statistical description of chaotic raysin a deep water acoustic waveguide,” J. Acoust. Soc. Am. (5), 2542–2552 (2007).[27] D. Makarov, S. Prants, A. Virovlyansky, and G. Zaslavsky,
Ray and wave chaos in oceanacoustics (Word Scientific, New Jersey, 2010).[28] A. Virovlyansky, “Matched shadow processing,” JASA Express Letters. J. Acoust. Soc.Am. (1), EL136–EL142 (2017).[29] L. Landau and E. Lifshitz,
Quantum mechanics (Pergamon Press, Oxford, 1977).2530] G. H. Golub and C. V. Loan,
Matrix computations (The John Hopkins University Press,Baltimore, 1989).[31] W. Munk and C. Wunsch, “Ocean acoustic tomography: A scheme for large scale mon-itoring,” Deep-Sea Res. , 123–161 (1979).[32] B. Howe, P. Worcester, and R. Spindel, “Ocean acoustic tomography: mesoscale veloc-ity,” J. Geophys. Res.92