Passive Imaging with Ambient Noise Under Wave Speed Mismatch: Mathematical Analysis and Wave Speed Estimation

P assiv e Imaging with Am bien t Noise Under W a v e Sp eed Mismatc h: Mathematical Analysis and W a v e Sp eed Estimation Zetao F ei ∗ and Josselin Garnier † F ebruary 18, 2026 Abstract It is known that wa ves generated by am bient noise sources and recorded b y passiv e receiv ers can b e used to image the reflectivities of an unknown medium. Ho wev er, reconstructing the reflectivit y of the medium from partial b oundary measuremen ts remains a challenging problem, particularly when the bac kground wa ve speed is unknown. In this paper, w e in vestigate passiv e correlation-based imaging in the da ylight configuration, where uncontrolled noise sources illuminate the medium and only am bient fields are recorded by a sensor arra y . W e first analyze dayligh t migration for a p oint reflector em b edded in a homogeneous background. By in tro ducing a searching w av e speed into the migration functional, w e deriv e an explicit characterization of the deterministic shift and defocusing effects induced by wa v e-sp eed mismatc h. W e sho w that the maximum of the env elop e of the resulting functional pro vides a reliable estimator of the true w av e speed. W e then extend the analysis to a random medium with correlation length smaller than the wa v elength. Lev eraging the shift formula obtained in the homogeneous case, w e introduce a virtual guide star that remains fixed under migration with different searc hing sp eeds. This prop erty enables an effective wa ve-speed estimation strategy based on spatial a veraging around the virtual guide star. F or both homogeneous and random media, we establish resolution analyses for the prop osed w av e-speed estimators. Numerical experiments are conducted to v alidate the theoretical result. 1 In tro duction P assive correlation–based imaging reconstructs information ab out a medium or subsurface using only wa ve- fields transmitted by opp ortunistic or ambien t noise sources. This tec hnique w as first demonstrated exp eri- men tally in [1] and has since gained widespread use in v arious fields, including helioseismology [2], v olcano monitoring [3], and reservoir monitoring [4]. F or a broad interdisciplinary review of correlation-based meth- o ds, we refer the reader to [5]; an extensive ov erview of the mathematical foundations can b e found in the monograph [6]. It is now well understoo d that key information—such as trav el times and Green functions—is enco ded in the cross-correlation of recorded signals [7]–[10]. This observ ation underlies the success of imaging reflectors using migration of the empirical cross-correlation. In [11], the authors developed a unified mathematical framew ork for analyzing passiv e imaging across diverse configurations. The uniqueness of passiv e imaging in verse problems has also b een established in both the time and frequency domains; see, for instance, [12], [13]. In this pap er, w e fo cus on passive imaging in the daylight c onfigur ation , where a sensor arra y is lo cated b et w een randomly distributed noise sources and the medium of interest. This setting is common in the literature on virtual-source imaging, where the receiv er acts as a virtual source. The scattered w av es recorded at the arra y are numerically bac k-propagated as if emitted from the receiver, allo wing wa vefield reconstruction and improv ed fo cusing in complex media. This metho d is widely used in seismic exploration [14]–[16], and ∗ CMAP , Ecole Polytec hnique, Institut Polytechnique de Paris, 91120 Palaiseau, F rance (zetao.fei@polytechnique.edu). † CMAP , CNRS, Ecole Polytec hnique, Institut Polytec hnique de Paris, 91120 P alaiseau, F rance (jos- selin.garnier@polytechnique.edu). This work was partially supp orted by Agence Nationale de la Recherc he (Grant No. ANR-23-CE30-0021). 1 its mathematical foundations are discussed in [17]. W e first in vestigate the case of a p oin t-like reflector em b edded in a homogeneous bac kground, a geometry closely related to classical reflector-imaging problems. Cross-correlation based method in this setting has been v alidated in physical experiments [18] and widely applied in seismology , e.g., [19]–[21]. A detailed resolution analysis of this approach w as provided in [22], [23]. W e then extend the prop osed metho d to the a second scenario in whic h the medium is random, with a correlation length muc h smaller than the wa velength. W e emphasize that the setup considered in this paper differs from typical virtual-source imaging: in our case, there is no controlled source arra y , and the imaging functional in dayligh t geometry takes a different form. Nev ertheless, the underlying idea remains similar— reco vering information ab out the medium of interest from the scattered field generated under uncon trolled (da ylight) illumination. A ma jor fo cus of this w ork concerns the estimation of the bac kground w av e speed from partial passiv e b oundary measurements. Accurate kno wledge of the bac kground w av e sp eed is crucial in passiv e imaging, b oth for properly fo cusing the image and for interpreting medium prop erties. Eik onal-based metho ds are widely used for wa v e-sp eed estimation in ambien t noise tomograph y [24], [25]; see also the review [26]. Ho wev er, their implemen tation in a da ylight configuration is challenging. Optimization-based approaches suc h as F ull W av eform Inv ersion (FWI) [27], [28] offer an alternative strategy by minimizing the discrepancy b et w een recorded and simulated data, and recent efforts ha ve extended such methods to quantitativ e passive imaging [29]. Nonetheless, these techniques remain computationally demanding and suffer from in trinsic non- con vexit y . In the literature, cross-correlation–based metho ds ha ve b een applied to image media with three- dimensional rapid fluctuations and a one-dimensional slo wly v arying bac kground in the parab olic regime [30]. Ho wev er, to our knowledge, few studies ha ve considered estimating the effectiv e w av e sp eed in a random medium under passiv e, da ylight-t yp e illumination. In the con text of active imaging, recen t developmen ts in ab erration correction and w av e-sp eed estimation using reflection-matrix analysis ha ve introduced techniques that do not require a strong guide star (i.e. a p oin t-lik e strong reflector). These methods exploits the recorded backscattered field by generating a syn thetic fo cus that acts as a virtual guide star [31]–[33]. This idea motiv ates the in tro duction of a virtual guide star in our passive imaging setting. A mathematical framew ork for in terpreting such exp erimen ts was recently prop osed in [34], and a related concept has b een explored in passiv e imaging in [35]. The contributions of this pap er are summarized below. First, in the da ylight configuration with a p oint- lik e reflector in a homogeneous background, we introduce a searc hing wa ve sp eed c s in to the traditional da ylight imaging functional and deriv e a quantitativ e characterization of the effects of w av e-sp eed mismatc h. Our analysis shows how an incorrect wa ve sp eed leads to shifts and blurring in the reconstructed reflector lo cation. W e further demonstrate that the maximum of the en velope of the generalized imaging functional yields a robust estimator for the bac kground wa ve speed. Second, w e extend this metho dology to the random- medium setting, where no strong guide star is present. By migrating the empirical cross-correlation using differen t searc hing wa ve speeds, we construct a virtual guide star whose ph ysical location is stabilized by the analytical shifting form ula established in the homogeneous case (see (3.8)). The background wa ve sp eed is then estimated from spatial av erages of the imaging functional around this virtual guide star. F or both settings, w e pro vide a comprehensive theoretical analysis, including resolution properties, and v alidate the prop osed estimators through n umerical exp erimen ts. The remainder of the pap er is organized as follows. Section 2 in tro duces the mathematical framew ork and reviews k ey prop erties of cross-correlation in passiv e imaging. Section 3 analyzes the generalized dayligh t migration imaging functional in the presence of a p oin t-like reflector and presents a w av e-sp eed estimator. Section 4 extends the method to random media and introduces the virtual guide-star–based estimator. Numerical experiments are given in Section 5, and concluding remarks app ear in Section 6. 2 Preliminaries In this section, we presen t the mathematical mo dels and assumptions underlying the problem studied in this pap er. Section 2.1 provides a detailed description of the problem setup. In Section 2.2, w e briefly review key prop erties of cross-correlation, a fundamental to ol in passive imaging. F or a more comprehensive discussion, w e refer the reader to [6]. 2 2.1 Problem Outline W e consider the scalar wa vefield u gov erned by the inhomogeneous acoustic w av e equation 1 c 2 ( x ) ∂ 2 u ∂ t 2 − ∆ x u = n ( t, x ) , (2.1) where c ( x ) denotes the spatially v arying wa ve sp eed and n ( t, x ) is a random noise field. The noise is assumed to b e zero-mean and stationary in time, with frequency band B cen tered at ω c and bandwidth B . Its second-order statistics are sp ecified b y E [ n ( t 1 , y 1 ) n ( t 2 , y 2 )] = F ε ( t 2 − t 1 ) Γ( y 1 , y 2 ) , (2.2) where F ε enco des the temp oral correlation and satisfies F ε (0) = 1. Its F ourier transform, b F ε ( ω ) = Z R F ε ( t ) e iω t dt, is real-v alued, even, nonnegativ e, and prop ortional to the source p ow er spectral densit y . W e assume that the decoherence time of the sources is short relativ e to the t ypical trav el times. Introducing a small scale parameter ε ≪ 1, we write F ε ( t ) = F  t ε  , (2.3) so that the central frequency scales as ω c = ω 0 /ε , and the bandwidth scales as B = B H /ε . Spatially , the sources are taken to b e delta-correlated: Γ( y 1 , y 2 ) = K ( y 1 ) δ ( y 1 − y 2 ) , where K describ es their spatial distribution. Medium model. The w av e sp eed is describ ed as a small p erturbation of a homogeneous bac kground: 1 c 2 ( x ) = 1 c 2 0 (1 + ρ ( x )) . W e consider t wo types of p erturbations: • Poin t-like reflector. A compact reflector lo cated near z r is modeled by ρ ( x ) = σ r V r ( x ) , V r ( x ) = 1 Ω r ( x − z r ) , where Ω r has characteristic size ℓ r . The p oin t-reflector regime corresp onds to ℓ r ≪ λ , and the reflector is assumed w eak, σ r ≪ 1. Neither the reflector position z r nor the bac kground v elo cit y c 0 is assumed to be known. • Random medium. In this regime, the p erturbations are supported in a region D and mo deled b y a rapidly oscillatory random field: ρ ( x ) = ρ µ ( x ) = µ  x ℓ c  , x ∈ D , where ℓ c ≪ λ is the correlation length. The random field µ is assumed to b e stationary , mean-zero, and has co v ariance function Σ( x ) = Co v ( µ ( x ) , µ (0)) , deca ying at the rate | Σ( x ) | ≲ (1 + | x | ) − 4 − η for some η > 0. The fluctuations are assumed weak, | Σ(0) | ≪ 1. While the mean-zero assumption is restrictive theoretically , it is not an obstacle in practice b ecause the b oundary reflections due to a nonzero mean can b e mitigated or are negligible in the imaging geometry presented b elow. 3 Data acquisition geometry . W e work in a passiv e daylight configuration in which an array A records the signals u ( t, x j ), j = 1 , . . . , N , and is p ositioned betw een the random noise sources (i.e. the support of K ) and the medium of interest (i.e. z r or D ). The short decoherence time implies a high central frequency . The arra y radius a is assumed larger than the w a velength λ , but smaller than the distance to the reflector | z r | . Under the scaling (2.3), we write ω c = ω 0 ε , λ = ε 2 π c 0 ω 0 , a = ε 1 / 2 a 0 , and parameterize the arra y as A = a A 0 for some reference geometry A 0 . Here, w e assume that the typical tra vel distance is muc h larger than the w av elength with dist( z r , A ) ∼ 1 or dist( D , A ) ∼ 1. The primary goal of this pap er is to prop ose estimators for the background wa ve sp eed c 0 for b oth configurations and to lo calize the p oin t-lik e reflector in the first configuration. Figure 1 illustrates the da ylight configuration for the point-reflector and random-medium setups. A z r (a) Dayligh t configuration for point-lik e reflector. A D ℓ c (b) Dayligh t configuration for random medium. Figure 1 : Dayligh t configuration for passiv e imaging of tw o differen t types of medium. 2.2 Cross-correlation and imaging functional W e start the introduction to the cross-correlation with the representation of the wa ve function (2.1). By the form ulation provided previously , the solution to the w av e equation has the follo wing integral representation: u ( t, x ) = Z Z dsd y n ( t − s, y ) G ( s, x , y ) . (2.4) The function G is the time-dep enden t causal Green’s function, i.e., the fundamental solution to the equation ( 1 c 2 ( x ) ∂ 2 G ∂ t 2 − ∆ x G = δ ( t ) δ ( x − y ) G ( t, x , y ) = 0 , ∀ t < 0 . (2.5) W e also in tro duce the time-harmonic Green’s function b G 0 for the background homogeneous medium, i.e., the solution to the Helmholtz equation    ∆ x b G 0 + ω 2 c 2 0 b G 0 = − δ ( x − y ) , lim | x |→∞ | x |  x | x | · ∇ x − i ω c 0  b G 0 ( ω , x , y ) = 0 . (2.6) W e consider the empirical cross-correlation ov er a time interv al (0 , T ) of the signal recorded at x 1 , x 2 ∈ A defined b y C T ( τ , x 1 , x 2 ) = 1 T Z T 0 u ( t, x 1 ) u ( t + τ , x 2 ) dt. (2.7) The statistical stabilit y of the ab o ve quantit y is w ell established, as sho wn in the following proposition [6]: Prop osition 2.1. (1) The exp e ctation of the empiric al cr oss-c orr elation C T (with r esp e ct to the statistic al law of the sour c es) is indep endent of T : E C T ( τ , x 1 , x 2 ) = C (1) ( τ , x 1 , x 2 ) . (2.8) 4 The the or etic al (or statistic al) cr oss-c orr elation C (1) is given by C (1) ( τ , x 1 , x 2 ) = 1 2 π Z Z d y dω b F ( ω ) K ( y ) b G ( ω , x 1 , y ) b G ( ω , x 2 , y ) e − iω τ , (2.9) wher e b G ( ω , x , y ) is the time-harmonic Gr e en ’s function, that is, the F ourier tr ansform of G ( t, x , y ) , and the overline stands for c omplex c onjugation. (2) The empiric al cr oss-c orr elation C T is a self-aver aging quantity: C T ( τ , x 1 , x 2 ) T →∞ − → C (1) ( τ , x 1 , x 2 ) , (2.10) in pr ob ability with r esp e ct to the statistic al law of the sour c es. With the recorded signals { u ( t, x j ) } t ∈ R , j =1 ,...,N , we compute the corresp onding cross-correlations. The in verse problem considered in this pap er is then reformulated as extracting information about the medium from the data set { C T } . T o this end, we in tro duce the gener alize d daylight migr ation imaging functional I ( z s , c s ) = N X j,l =1 C T ( T s ( z s , x l ) + T s ( z s , x j ) , x j , x l ) , (2.11) where z s denotes the search point, c s is the sound sp eed used in migration, and T s ( z s , x ) is the generalized tra vel time defined by T s ( z s , x ) := | z s − x | c s . F urthermore, w e denote the ratio of c s and c 0 b y v = c s /c 0 . When c s = c 0 (i.e. v = 1), T s coincides with the tra vel time in a homogeneous medium, T 0 ( z s , x ) := | z s − x | c 0 , and the imaging functional reduces to the classical dayligh t migration functional. In the dense-arra y limit, the imaging functional can be approximated by I ( z s , c s ) = N 2 |A| 2 Z Z A 2 d x 1 d x 2 C T ( T s ( z s , x 2 ) + T s ( z s , x 1 ) , x 1 , x 2 ) . (2.12) The construction of the imaging functional motiv ates the analysis of the theoretical cross-correlation C (1) . Under the assumption that the p erturbation of the background medium is weak, we may inv ok e the Born appro ximation for the Green’s function, b G ( ω , x , y ) ≃ b G 0 ( ω , x , y ) + ω 2 c 2 0 Z d z e ρ ( z ) b G 0 ( ω , x , z ) b G 0 ( ω , z , y ) , (2.13) where e ρ ( z ) = ( σ r ℓ 3 r δ ( z − z r ) , point-lik e reflector , ρ µ ( z ) , random medium . Com bining (2.9), (2.13), and the scaling regime introduced in the previous section, we obtain the decom- p osition C (1) ( τ , x 1 , x 2 ) = C (1) 0 ( τ , x 1 , x 2 ) + C (1) i ( τ , x 1 , x 2 ) + C (1) ii ( τ , x 1 , x 2 ) , (2.14) where C (1) 0 ( τ , x 1 , x 2 ) = 1 2 π Z Z d y dω K ( y ) b F ( ω ) b G 0  ω ε , x 1 , y  b G 0  ω ε , x 2 , y  e − i ωτ ε , C (1) i ( τ , x 1 , x 2 ) = 1 2 π c 2 0 ε 2 Z d z e ρ ( z ) Z Z d y dω ω 2 K ( y ) b F ( ω ) b G 0  ω ε , x 2 , z  b G 0  ω ε , z , y  b G 0  ω ε , x 1 , y  e − i ωτ ε , C (1) ii ( τ , x 1 , x 2 ) = 1 2 π c 2 0 ε 2 Z d z e ρ ( z ) Z Z d y dω ω 2 K ( y ) b F ( ω ) b G 0  ω ε , x 1 , z  b G 0  ω ε , z , y  b G 0  ω ε , x 2 , y  e − i ωτ ε . 5 The term C (1) 0 corresp onds to the con tribution of the direct w a ves, and it inv olves the rapid phase ω ε  | x 2 − y | c 0 − | x 1 − y | c 0 − τ  . F or fixed x 1 , x 2 ∈ A , this con tribution is supported on the time in terv al [ − ∆ τ 0 , ∆ τ 0 ], where ∆ τ 0 := sup y ∈ supp( K )   | y − x 1 | − | y − x 2 |   c 0 ≤ | x 1 − x 2 | c 0 ≤ diam( A ) c 0 , whic h is of order O ( ε 1 2 ) in the regime under consideration. In con trast, the scattered-w av e contributions C (1) i and C (1) ii in volv e the rapid phases ω ε  | x 2 − z | c 0 + | z − y | c 0 − | x 1 − y | c 0 − τ  and ω ε  | x 2 − y | c 0 − | z − y | c 0 − | x 1 − z | c 0 − τ  , resp ectiv ely . Recall that w e consider the da ylight configuration in this pap er, these terms are supp orted on time in terv als satisfying | τ | ≍ 2 dist( z r , A ) c 0 or | τ | ≍ 2 dist( D , A ) c 0 , b oth of which are quantities of order one. Consequen tly , within the considered regime, the direct and scattered con tributions to the cross-correlation are well separated in time. Therefore, the comp onen t C (1) 0 can be effectively remov ed from the recorded data by time-windowing. In the remainder of the analysis, w e fo cus exclusively on the scattered-w av e con tributions. With a slight abuse of notation, we therefore redefine C (1) ( τ , x 1 , x 2 ) = C (1) i ( τ , x 1 , x 2 ) + C (1) ii ( τ , x 1 , x 2 ) . (2.15) 3 Guide star assisted w a v e sp eed estimation In this section, we analyze the prop erties of the imaging functional I ( z s , c s ) defined in (2.11) when a point- lik e reflector is present in the medium. The analysis quantitativ ely c haracterizes the influence of w a ve sp eed mismatc h on the p oin t spread function, including the phenomenon of fo cusing p oin t shift. Based on these insigh ts, we prop ose metho ds for estimating the background wa ve sp eed c 0 and localizing the p oin t-like reflector. 3.1 Analysis of the imaging functional T o present the properties of the imaging functional, we first in tro duce a co ordinate frame ( b e 1 , b e 2 , b e 3 ) satis- fying • the sensor arra y is cen tered around 0 and A ⊂ span { b e 1 , b e 2 } , • the reflector z r ∈ span { b e 1 , b e 3 } with co ordinates z r = ( z r, 1 , 0 , z r, 3 ). Recall the regime we adopt for the sensor arra y , for x 1 ∈ A , w e write x 1 = ε 1 2 f x 1 = ε 1 2 ( g x 1 , 1 , g x 1 , 2 , 0). W e also in tro duce the following orthonormal frame { b f 1 , b f 2 , b f 3 } : b f 1 = 1 | z r | ( z r, 3 b e 1 − z r, 1 b e 3 ) , b f 2 = b e 2 , b f 3 = 1 | z r | ( z r, 1 b e 1 + z r, 3 b e 3 ) . (3.1) Differen t from the traditional da ylight migration imaging functional, in I ( z s , c s ), the sp eed of sound used for backpropagation c s can hav e a mismatch with c 0 . This mismatch can lead to a shift in the fo cusing p oin t of the image deriv ed from I ( z s , c s ). A similar phenomenon is observed in the active imaging configuration, 6 see [34]. W e denote the focusing point in the image as z f ( v ), and z f when there is no am biguity . W e define the map ϕ v : R 3 → R 3 as ϕ v ( z r ) =  v 2 z ⊥ r , v q ( z q r ) 2 + (1 − v 2 ) | z ⊥ r | 2  . (3.2) When the p oin t-like reflector is p ositioned at z r , w e sho w the fo cusing point is located at z f = ϕ v ( z r )m with v = c s /c 0 under certain condition in Theorem 3.1. Before diving in to the mathematical calculation, we first in tro duce an orthonormal frame { b g 1 , b g 2 , b g 3 } based on z f as follo ws: b g 1 = 1 | z f | ( z f , 3 b e 1 − z f , 1 b e 3 ) , b g 2 = b e 2 , b g 3 = 1 | z f | ( z f , 1 b e 1 + z f , 3 b e 3 ) . (3.3) W e illustrate the three co ordinate systems in Figure 2. A z r b e 1 b e 3 0 b f 1 b f 3 z f b g 1 b g 3 Figure 2 : The sensor arra y and three reference frames ( b e 1 , b e 3 ), ( b f 1 , b f 3 ), and ( b g 1 , b g 3 ). The unit vector b e 2 = b f 2 = b g 2 are orthogonal to the plane of the figure. In the case that c s = c 0 , it is known that the resolution of the image dep ends on the tilt of the arra y relativ e to b f 3 , see [6]. In the analysis, w e observ e that the imaging function also dep ends on the tilt of the arra y relative to b g 3 . Therefore, for the conv enience of presentation, w e introduce the cosine of the tilt angles: α r = D b f 3 , b e 3 E = z r, 3 | z r | , α f = ⟨ b g 3 , b e 3 ⟩ = z f , 3 | z f | , (3.4) where z f is giv en b y W e define the normalized point spread function as G α r ,α f ( ξ 1 , ξ 2 , ξ 3 , ξ 4 ) = 1 |A 0 | Z A 0 du 1 du 2 exp  − i ( α f u 1 ξ 1 + u 2 ξ 2 ) − i u 2 1 α 2 r + u 2 2 2 ξ 3 − i u 2 1 + u 2 2 2 ξ 4  . (3.5) In the theoretical results of the imaging functional, { ξ j } corresp ond to different physical meanings. ξ 1 and ξ 2 is related to the transv erse direction, ξ 3 corresp onds to the range direction and the influence of the w av e sp eed mismatch is reflected through ξ 4 . W e observe that the function G α r ,α f is a generalization of the point spread function G α r defined in Section 6, [6]. The tw o functions hav e similar qualitative b ehaviors when ξ 3 = 0 (or ξ 4 = 0). When ξ 3 , ξ 4  = 0, the behavior of the function G α r ,α f is complicated, which leads to the complicated coupling phenomenon of the range comp onent and the mismatch of the wa ve sp eed in the p oin t spread function, w hic h is discussed in 3.2. Figure 3 illustrates the b eha vior of G α r ,α f through an example. In the follo wing, we show the asymptotic expression for the generalized imaging function under the regime w e consider. Through stationary phase analysis, w e quan tify the dep endence of the imaging functional on the follo wing factors: • the normalized arra y radius a 0 , • the distance of the reflector to the cen ter of the sensor array | z r | , • the cen tral frequency ω c = ω 0 ε , • the effectiv e bandwidth B H ε , 7 -10 -5 0 5 10 1 -10 -5 0 5 10 3 -10 -5 0 5 10 1 -10 -5 0 5 10 -10 -5 0 5 10 1 -10 -5 0 5 10 0.2 0.4 0.6 0.8 Figure 3 : W e pic k the example that z r = ( − 5 , 0 , 50), v = 1 . 1 and calculate α r and α f according to (3.4) with z f = ϕ v ( z r ). W e plot the function G α r ,α f ( ξ 1 , 0 , ξ 3 , 0) showing its real part, imaginary part, and mo dulus from left to right. • the c hosen sp eed of sound for migration c s . Theorem 3.1. L et us assume that b F ( ω ) has the form b F ( ω ) = 1 B H  b F 0  ω 0 − ω B H  + b F 0  ω 0 + ω B H  , (3.6) wher e ω 0 , B H > 0 . The gener alize d daylight imaging functional is non-ne gligible (i.e. it pr esents a wel l- define d p e ak) when the se ar ching wave sp e e d c s satisfies v = c s c 0 ≥ 1 −  z r, 3 | z ⊥ r |  2 . (3.7) F or a given c s that satisfying (3.7), as a function of z s , the imaging functional I ( z s , c s ) define d by (2.11) has the fo cusing p oint z f = ϕ v ( z r ) =  v 2 z ⊥ r , v q ( z q r ) 2 + (1 − v 2 ) | z ⊥ r | 2  . (3.8) A t the se ar ching p oint z s : z s = z f + r s = z f + ε 1 2 v 2 η 1 b g 1 + ε 1 2 v 2 η 2 b g 2 + εv η 3 b g 3 , (3.9) the imaging functional has the asymptotic form (as ε → 0 ): I ( z s , c s ) ≈ σ r ℓ 3 r N 2 2 6 π 3 c 0 ε K ( 0 , z r ) | z r | 2 · P ( r s , c s ) , (3.10) wher e K ( x , z ) = Z ∞ 0 K  x + x − z | x − z | · l  dl, (3.11) and P ( r s , c s ) is the p oint spr e ad function define d by P ( r s , c s ) = Re  Z B dω  − iω b F ( ω )  exp  − i ω c 0  2 η 3 + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |  × G 2 α r ,α f − a 0 ω c 0 | z r | η 1 , − a 0 ω c 0 | z r | η 2 , 0 , − a 2 0 ω c 0 | z r |  c 0 c s  2 − 1 !!) . (3.12) 8 If additional ly B H ≪ ω 0 , then the p oint spr e ad function has the form P ( r s , c s ) ≃ 4 π ω 0 Re  − i exp  − i ω 0 c 0  2 η 3 + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |  F 0  − B H c 0  2 η 3 + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |  × G 2 α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , 0 , − a 2 0 ω 0 c 0 | z r |  c 0 c s  2 − 1 !!) , (3.13) whose slow ly varying envelop e is given by P E ( r s , c s ) ≃ 2 π ω 0     F 0  − B H c 0  2 η 3 + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |      ×      G α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , 0 , − a 2 0 ω 0 c 0 | z r |  c 0 c s  2 − 1 !!      2 . (3.14) Pr o of. See App endix A.1. Here, we emphasize that the expression abov e is v alid when the bandwidth of the noise sources is not v ery small (compared to the order O ( ε )). When the bandwidth is of order O ( ε ), the expression of the imaging functional is sligh tly differen t from the previous one. Prop osition 3.2. L et us assume that b F ( ω ) has the form b F ( ω ) = 1 εB 0  b F 0  ω 0 − ω εB 0  + b F 0  ω 0 + ω εB 0  , (3.15) wher e ω 0 , B 0 > 0 . Under the c ondition (3.7), at the se ar ching p oint z s : z s = z f + r s = z f + ε 1 2 v 2 η 1 b g 1 + ε 1 2 v 2 η 2 b g 2 + v η 3 b g 3 , (3.16) such that | r s | ≪ | z r | , the imaging functional has the asymptotic form (as ε → 0 ): I ( z s , c s ) ≈ σ r ℓ 3 r N 2 2 6 π 3 c 0 ε K ( 0 , z r ) | z r | 2 · e P ( r s , c s ) , (3.17) wher e e P ( r s , c s ) is the p oint spr e ad function define d by e P ( r s , c s ) =4 π ω 0 Re  − i exp  − i ω 0 c 0  2 η 3 ε + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |  F 0  − 2 B 0 η 3 c 0  ×G 2 α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , − a 2 0 ω 0 c 0 | z r | 2 η 3 , − a 2 0 ω 0 c 0 | z r |  c 0 c s  2 − 1 !!) . (3.18) If B 0 is lar ger than the critic al value B c := ω 0 2 a 2 0 | z r | 2 , (3.19) then the p oint spr e ad function has the form e P ( r s , c s ) =4 π ω 0 Re  − i exp  − i ω 0 c 0  2 η 3 ε + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |  F 0  − 2 B 0 η 3 c 0  ×G 2 α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , 0 , − a 2 0 ω 0 c 0 | z r |  c 0 c s  2 − 1 !!) , (3.20) 9 whose slow ly varying envelop e is given by P E ( r s , c s ) = 2 π ω 0     F 0  − 2 B 0 η 3 c 0      ·      G α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , 0 , − a 2 0 ω 0 c 0 | z r |  c 0 c s  2 − 1 !!      2 . (3.21) If B 0 is smal ler than the critic al value B c , then the p oint spr e ad function has the form e P ( r s , c s ) =4 π ω 0 Re  − i exp  − i ω 0 c 0  2 η 3 ε + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |  F 0 (0) ×G 2 α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , − a 2 0 ω 0 c 0 | z r | 2 η 3 , − a 2 0 ω 0 c 0 | z r |  c 0 c s  2 − 1 !!) , (3.22) whose slow ly varying envelop e is given by P E ( r s , c s ) = 2 π ω 0 | F 0 (0) | ·      G α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , − a 2 0 ω 0 c 0 | z r | 2 η 3 , − a 2 0 ω 0 c 0 | z r |  c 0 c s  2 − 1 !!      2 . (3.23) Pr o of. See App endix A.2. Here, we discuss and interpret the results derived in the theorem and proposition abov e. F or c s = c 0 , we reco ver the w ell-kno wn result for the resolution analysis of the da ylight passiv e imaging configuration. W e refer to a more detailed in tro duction to [6]. F or completness, w e summarize the results as follows: • the frames { b f 1 , b f 2 , b f 3 } and { b g 1 , b g 2 , b g 3 } coincide, and the image is cen tered at z r , • the cross-range resolution is: λ | z r | aα r in the b f 1 direction and λ | z r | a in the b f 2 direction, • the range resolution is εc 0 2 B H in the broadband case ( B H ε ≫ B c ) and λ | z r | 2 a 2 in the narro wband case ( B H ε ≪ B c ). The effect caused by the mismatch of c s and c 0 can be summarized in the follo wing. • When the p oin t-like reflector has mo derate p erp endicular distance to the b e 3 axis, i.e, the condition 3.7 is satisfied, the image is not negligible. • In the broadband case, the image fo cuses at the p oin t z f = ϕ v ( z r ). The vectors b g 1 and b g 2 define the transv erse directions of the image, and b g 3 defines the range direction. F urthermore, the resolution is the same as describ ed in the case c s = c 0 . • In the narrowband case, the mismatch couples with the range direction behavior, which causes further shift of the exact fo cusing p oint of the image in the range direction. W e illustrate this phenomenon by plotting the function |G α r ,α f | 2 for three differen t c hoices of c s in Figure 4. • The mismatc h causes the reduction of the p eak v alue at the fo cusing p oin t for all cases except for the narro wband case with α r = 1. This fact pav es the wa y for the estimation of c 0 in the general case as discussed in the next section. • The narrowband case with α r = 1 is a sp ecial case where the range comp onen t couples with the mismatc h in the follo wing w ay: P E ( r s , c s ) ∼      G α r ,α f − a 0 ω 0 c 0 | z r | η 1 , − a 0 ω 0 c 0 | z r | η 2 , − a 2 0 ω 0 c 0 | z r | 2 η 3 + | z r |  c 0 c s  2 − 1 !! , 0 !      2 . This means that along the line defined b y η 1 = η 2 = 0 and η 3 + | z r |  ( c 0 /c s ) 2 − 1  = 0, the p oint spread function remains at its maximum. Consequently , in this particular case, the p eak of the p oin t spread function env elop e cannot b e used to estimate c 0 . Therefore, this scenario is unfav orable for sim ultaneously estimating the wa ve sp eed and imaging the unknown reflector. How ever, if either the bac kground sp eed or the reflector p osition is kno wn, it is still p ossible to estimate the remaining unkno wn parameter. 10 -5 0 5 1 -30 -20 -10 0 10 20 3 -5 0 5 1 -30 -20 -10 0 10 20 -5 0 5 1 -30 -20 -10 0 10 20 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Figure 4 : W e pick the example that a 0 = 20, ω 0 = 4, c 0 = 1, and z r = (10 , 0 , 20 √ 2). W e plot the profile of     G α r ,α f  − a 0 ω 0 c 0 | z r | η 1 , 0 , − a 2 0 ω 0 c 0 | z r | 2 η 3 , − a 2 0 ω 0 c 0 | z r |   c 0 c s  2 − 1      2 for three differen t choices of c s : c s = 0 . 8 (left), c s (middle) and c s = 1 . 2 (righ t). The peaks of the three plots takes v alue 0 . 87, 1 and 0 . 96 resp ectiv ely . 3.2 Estimator for the w av e sp eed Previously , w e ha ve derived the asymptotic expression of the generalized da yligh t migration imaging func- tional when there is a mismatc h b et ween the chosen wa ve sp eed for bac kpropagation and the true wa v e sp eed in the medium. Based on that, we propose the estimators for the p osition of the reflector and the w av e sp eed c 0 except for the narrowband case with α r = 1. W e start with the following corollary . Corollary 3.3. The function P E define d in (3.14), (3.23) and (3.21) achieves its glob al maximum at r s = 0 , c s = c 0 . Pr o of. W e only need to notice that the function |G α r ,α f ( ξ 1 , ξ 2 , ξ 3 , ξ 4 ) | achiev es its global maximum at ξ j = 0, j = 1 , 2 , 3 , 4, and the function F 0 ( ω ) ac hiev es its maxim um at ω = 0. Based on the corollary ab o ve, we prop ose the follo wing estimator: Estimator. ( b z r , b c ) = argmax z s ,c s En velope {I ( z s , c s ) } . (3.24) T o implement the proposed estimator, w e first note that the generalized dayligh t migration imaging functional do es not directly provide its env elop e. T o extract the env elop e, we apply the Hilb ert transform to I ( z s , c s ). Specifically , we discretize the search domain in spherical co ordinates. F or eac h fixed pair of angular v ariables, the Hilb ert transform is p erformed along the radial direction. This pro cess is then rep eated to encompass the en tire searc h domain. The corresp onding pseudo-code is pro vided b elo w. The resolution of the estimator for the reflector is discussed in Section 3.1. Similarly , based on the theoretical results, we observe that the relative resolution for the wa ve sp eed estimation is giv en by λ | z r | a 2 for b oth narrowband and broadband cases. W e discuss this result in the follo wing. • Compared with the relativ e cross-range resolution of reflector lo calization, given b y λ | z r | a , the relativ e resolution for w a ve speed estimation differs b y an additional factor of | z r | /a and the cosine of the tilt angle. • In the narro wband case, the relativ e resolution for w av e speed estimation coincides with the relativ e range resolution for reflector lo calization. In con trast, in the broadband case, the relativ e resolution for w av e sp eed estimation is p o orer than the relativ e range resolution for reflector lo calization. 11 Algorithm 1: Reflector lo cation and wa ve speed estimation Input: Search domain Ω ⊂ R 3 , a searc hing grid of c s , measuremen t data u ( t, x r ). Output: The estimators ( b z r , b c ) 1 Calculate the cross-correlation; 2 Discretize the searc h domain Ω in spherical co ordinates { ( r j , θ k , φ l ) } ; 3 foreach fixe d se ar ching wave sp e e d c s do 4 foreac h fixe d angular p air ( θ k , φ l ) do 5 P erform the Hilb ert transform of I (( r j , θ k , φ l ) , c s ) along the radial v ariable r ; 6 Store the resulting env elop e along the radial direction; 7 Assemble the env elop es obtained for all ( θ k , φ l ) to co v er the full searc h domain; 8 Calculate the estimators using (3.24); 9 return The estimators ( b z r , b c ); • The relativ e resolution scales in versely with the F resnel num b er, defined as F num ∼ a 2 λ | z r | . When F num ≫ 1, the resolution is high, consistent with the geometric optics approximation adopted in the analysis. When F num ≲ 1, diffraction effects b ecome significant, the geometric optics appro ximation breaks do wn, and the resolution deteriorates accordingly . 4 Effectiv e w a v e sp eed estimation through virtual guide star In this section, we dev elop a metho d for estimating the effective w a ve sp eed in random media, extending the guide-star strategy introduced earlier. Since a random medium do es not con tain a ph ysical guide star, w e instead in tro duce a virtual guide star at a lo cation of our o wn c hoice within the medium. This virtual guide star is not a physical reflector; rather, it is realized indirectly through a family of searching p oin ts in the imaging functional. The theoretical principles underlying this construction are presen ted in Section 4.1, and the corresponding practical estimator is describ ed in Section 4.2. 4.1 Probing the sp ec kle In practice, our goal is to prob e a fixed physical lo cation in the medium, denoted b y z g and referred to as the virtual guide star . The k ey observ ation is that a mismatch betw een the true propagation sp eed and the bac kpropagation sp eed c s causes a reflector lo cated at z r ∈ D to fo cus at a shifted lo cation z f = ϕ v ( z r ) in the imaging domain. F or any fixed v , the mapping ϕ v : R 3 → R 3 , defined in (3.2), is bijectiv e. T o prob e the medium appropriately , w e pro ceed as follo ws. F or a giv en backpropagation speed c s , we prescrib e a trav el time t g and define the corresp onding searc hing p oin t in the imaging functional as z s ( c s , t g ) = (0 , 0 , c s · t g ) . With this c hoice, the associated ph ysical lo cation of the virtual guide star is given by z g = ϕ − 1 v  z s ( c s , t g )  = (0 , 0 , c 0 · t g ) , whic h is indep enden t of c s . Therefore, although the searching point z s v aries with the assumed bac kpropa- gation speed, it alw ays corresponds to the same physical lo cation z g in the medium. In principle, the virtual guide star can b e placed an ywhere within the domain D by a similar procedure describ ed abov e. F or simplicity , we restrict our atten tion to lo cations along the b e 3 -axis. This is perfectly 12 A b e 1 b e 3 0 z s z g D ℓ c Figure 5 : Illustration for the effectiv e wa ve sp eed estimation. F or an y fixed c s , we pick the searc hing point at z s , and the virtual guide star is p ositioned at z g = ϕ − 1 v ( z s ) with v = c s /c 0 . legitimate as random scatterers are distributed throughout the medium. The corresp onding geometric con- figuration is illustrated in Figure 5. The following theorem c haracterizes the b eha vior of the imaging functional at the searching p oints asso- ciated with the virtual guide star and the packpropagation wa ve sp eed. Theorem 4.1. L et us assume that b F ( ω ) has the form b F ( ω ) = 1 B H  b F 0  ω 0 − ω B H  + b F 0  ω 0 + ω B H  , (4.1) wher e ω 0 , B H > 0 . F or any fixe d c s and a pr escrib e d tr avel time t g , cho osing the se ar ch p oint z s ( c s , t g ) = (0 , 0 , c s · t g ) induc es a virtual guide star lo c ate d at z g = ϕ − 1 v  z s ( c s , t g )  = (0 , 0 , c 0 · t g ) , which is indep endent of c s . Denote the ve ctor η = ( η 1 , η 2 , η 3 ) , and r ( η ) = ( ε 1 2 η 1 , ε 1 2 η 2 , εη 3 ) , then the imaging functional has the asymptotic form (as ε → 0 ): I ( z s ( c s , t g ) , c s ) ≈ N 2 ε 2 6 π 3 c 0 K ( 0 , z g ) | z g | 2 · Z R 3 ρ µ ( z g ( t g ) + r ( η )) P µ ( t g , η , c s ) d η , (4.2) wher e P µ ( t g , η , c s ) is the p oint spr e ad function define d by P µ ( t g , η , c s ) = Re  Z B dω  − iω b F ( ω )  exp  i ω c 0  2 η 3 + ( η 2 1 + η 2 2 ) | z g |  ×G 2 1 , 1 a 0 ω c 0 | z g | η 1 , a 0 ω c 0 | z g | η 2 , 0 , − a 2 0 ω c 0 | z g |  c 0 c s  2 − 1 !!) . (4.3) If additional ly B H ≪ ω 0 , then the p oint spr e ad function has the form P µ ( t g , η , c s ) ≃ 4 π ω 0 Re  − i exp  i ω 0 c 0  2 η 3 + ( η 2 1 + η 2 2 ) | z g |  F 0  B H c 0  2 η 3 + ( η 2 1 + η 2 2 ) | z g |  × G 2 1 , 1 a 0 ω 0 c 0 | z g | η 1 , a 0 ω 0 c 0 | z g | η 2 , 0 , − a 2 0 ω 0 c 0 | z g |  c 0 c s  2 − 1 !!) . (4.4) Pr o of. See App endix B.1. Remark 4.2. In this se ction, we only c onsider the c ase wher e the b andwidth is not very smal l c omp ar e d to the or der of the wavelength, O ( ε ) . We exp e ct that the similar metho dolo gy c an b e extende d to the narr owb and c ase as wel l. 13 W e observ e that the mismatc h betw een c s and c 0 pla ys a role in the profile of the p oin t spread function, whic h is a random quantit y that dep ends on the realization of the random medium ρ µ . The following prop osition offers insight that w e can extract the information of c 0 from the second-order momen t of the imaging functional. Prop osition 4.3. (se c ond or der moment) F ol lowing the assumptions in The or em 4.1, we have E |I ( z s ( c s , t g ) , c s ) | 2 ∼ ℓ 3 c ∥ Σ ∥ L 1 ( R 3 ) · Z R 3 |P µ ( t g , η , c s ) | 2 d η . (4.5) If B H ≪ ω 0 , we have E |I ( z s ( c s , t g ) , c s ) | 2 ∼ Z R 3     F 0  B H c 0  2 η 3 + η 2 1 + η 2 2 | z g |      2 ×      G 1 , 1 a 0 ω 0 c 0 | z g | η 1 , a 0 ω 0 c 0 | z g | η 2 , 0 , − a 2 0 ω 0 c 0 | z g |  c 0 c s  2 − 1 !!      4 dη 1 dη 2 dη 3 . (4.6) Pr o of. See App endix B.2. Although the second order moment provides the information of c 0 , we should b e aw are that we access to only one realization of the random medium in practice. The follo wing t w o propositions sho w that by c hanging the searching p oin t with distance of the order of the wa velength, the imaging functional is stationary , which links the spatial av erage to the ensemble a v erage. Prop osition 4.4. (L o c al ly stationary) L et ∆ z = (∆ z ⊥ , ∆ z ∥ ) , and p (∆ z , c s ) =  ∆ z ⊥ v 2 , ∆ z ∥ v  , we have I ( z s ( c s , t g ) + ε ∆ z , c s ) ≃ N 2 ε 2 6 π 3 c 0 K ( 0 , z g ) | z g | 2 · Z R 3 ρ µ ( z g + ε p (∆ z , c s ) + r ( η )) P µ ( t g , η , c s ) d η , (4.7) In p articular, I is stationary in the le ading or der. The pro of of the prop osition is almost the same calculation as the one in Theorem 4.1. No w, we are able to link the spatial and ensem ble av erages by showing the ergo dicit y . Prop osition 4.5. (sp atial aver age to ensemble aver age) L et S ( l ) = [ − l 2 , l 2 ] 3 , we have 1 |S ( l ) | Z S ( l ) |I ( z s ( c s , t g ) + ε ∆ z , c s ) | 2 d ∆ z l →∞ − → E |I ( z s ( c s , t g ) , c s ) | 2 a.s. (4.8) Pr o of. See App endix B.3. In summary , by calculating the spatial av erage of |I ( z s ( c s , t g ) + ε ∆ z , c s ) | 2 for ∆ z ∈ S ( l ) with sufficiently large l , we can appro ximate the ensemble av erage E |I ( z s ( c s , t g ) , c s ) | 2 , whic h contains the information of c 0 . Hence, the spatial av erage can be used to estimate the effective wa ve sp eed. 4.2 Effectiv e w a v e sp eed estimation In this section, we discuss the pro cedure of the effectiv e w av e sp eed estimation. W e start with the following c haracterization of the second order moment for the regime B H ≪ ω 0 . Corollary 4.6. F or E |I ( z s ( c s , t g ) , c s ) | 2 given in (4.6), we have c 0 = arg max c s E |I ( z s ( c s , t g ) , c s ) | 2 Pr o of. See App endix B.4. Then, it is straigh tforward to prop ose the following estimator for the effectiv e wa ve speed based on the spatial a verage of I ( z s ( c s , t g ) + ∆ z , c s ) as follo ws: 14 Estimator. F or an appr opriate choic e of l > 0 , b c eff = argmax c s 1 |S ( l ) | Z S ( l ) |I ( z s ( c s , t g ) + ε ∆ z , c s ) | 2 d ∆ z . (4.9) W e summarize the whole pro cedure of estimating the effective wa ve sp eed in the following pseudo-co de. Algorithm 2: Effective wa v e sp eed estimation Input: A searching grid of c s , a grid for spatial a verage { ∆ z j } N j =1 , an appropriate choice of t g , t ypical wa ve length in the medium ε , measuremen t data u ( t, x r ). Output: The estimator b c eff 1 Calculate the cross-correlation; 2 foreach fixe d se ar ching wave sp e e d c s do 3 Calculate the imaging p oin t z s ( c s , t g ) = (0 , 0 , c s t g ); 4 Calculate the quan tit y I = 1 N P j |I ( z s ( c s , t g ) + ε ∆ z j ) | 2 ; 5 Calculate the estimators using (4.9); 6 return The estimator b c eff ; Finally , w e notice that the formula (4.6) indicates that the realtive resolution of the effective wa v e speed estimator is giv en b y λ | z g | a 2 . 5 Numerical exp eriments In this section, we presen t the numerical exp erimen ts carried out for the estimators introduced in sections 3- 4. All sim ulations are p erformed in MA TLAB on a mac hine equipp ed with an M3 8-core CPU. F or reference, the computation time for each numerical exp erimen t is rep orted. W e note that all algorithms can b e further accelerated using high-performance GPUs, as the co de is primarily based on matrix and tensor op erations. W e assume that the av ailable data is the theoretical cross-correlation C (1) giv en in (2.9). Owing to the statistical stability of the cross-correlation, one may alternativ ely start from the raw sensor-arra y recordings, pro vided that the recording time window is sufficiently large. 5.1 Homogeneous medium with a p oin t reflector W e consider a homogeneous background medium with wa ve speed c 0 = 1. The distribution of the noise sources is described b y K ( x ) = exp  −  ( x 1 + 50) 2 500 + ( x 2 + 50) 2 500 + ( x 3 − 42 . 5) 2 10  . The p o wer sp ectral densit y is assumed to b e b F ( ω ) = ω 2 e − ω 2 . The sensor array lies in the plane span { e 1 , e 2 } with co ordinates (5 j, 5 k , 0) for j, k = − 2 , − 1 , . . . , 2. The p oin t reflector is p ositioned at ( − 5 , 0 , 50) with σ r l 3 r = 0 . 01. The numerical configuration is illustrated in Figure 6. F or generating the syn thetic wa ve data, w e employ the domain [ − 50 , 50] × [ − 50 , 50] × [ − 50 , 35] in the { b e 1 , b e 2 , b e 3 } coordinates, discretized in to a 101 × 101 × 16 grid for K ( x ). The frequency range is discretized as 0 : 5 × 10 − 3 : 4, and the time-lag v ariable τ is sampled on the grid 40 : 10 − 2 : 300. The data-generation stage tak es appro ximately 6.4 min utes. 5.1.1 Localization and w av e sp eed estimation W e c ho ose the imaging window to be a square of side length 30 in the plane span { e 1 , e 3 } , cen tered at the reflector, with a grid step size of 0 . 5. Figure 7 displa ys the numerical results of the generalized dayligh t migration imaging functional. W e present the cases c s = 0 . 8, c s = 1, and c s = 1 . 2. The black diamond indicates the theoretical fo cusing lo cation given by (3.8). The total computation time for pro ducing the 15 A z r e 3 e 1 − 50 0 50 − 50 0 50 Figure 6 : The sensor array is p ositioned at the cen ter of the plane span { e 1 , e 2 } , the reflector is p ositioned at ( − 5 , 0 , 50), and the circles on the left represen t the source region which is describ ed b y K ( x ). Figure 7 : Numerical simulation for different choices of c s . The black diamond represents the theoretically predicted fo cus point of the image, z f . Sp ecifically , the diamond in the middle graph shows the exact p osition of the reflector. three images is ab out 3 seconds. The sim ulation agrees w ell with the theoretical prediction, and the shift in the focusing p oin t as c s v aries is clearly observ able. Next, we conduct numerical exp erimen ts to estimate the wa ve sp eed. T o this end, we consider the profile of the function ϕ 1 ( c s ) := max z s En velope {I num ( z s , c s ) } , (5.1) obtained from n umerical sim ulations. Here, the imaging functional I num is defined as I num ( z s , c s ) := N 2 |A| 2 Z Z A 2 d x 1 d x 2 C (1)  T s ( z s , x 2 ) + T s ( z s , x 1 ) , x 1 , x 2  , (5.2) where C (1) is giv en in (2.9). The search domain for z s is discretized in spherical co ordinates, with radius ranging from 30 to 70, p olar angle from 0 to 12 ◦ , and azimuthal angle from − 8 ◦ to 8 ◦ , using a 300 × 80 × 80 grid. The wa ve speed is sampled from 0 . 7 to 1 . 3 with step size 0 . 02. The en tire computation tak es appro ximately 3 min utes. F or comparison, w e also compute the theoretical prediction ϕ 2 ( c s ) := En velope {P (0 , c s ) } ≃      Z ∞ 0 dω  − iω b F ( ω )  G 2 α r ,α f 0 , 0 , 0 , − a 2 0 ω c 0 | z r |  c 0 c s  2 − 1 !!      . Figure 8 displays b oth ϕ 1 ( c s ) and ϕ 2 ( c s ), eac h appropriately normalized. The theoretical and n umerical curv es agree w ell, and the w av e sp eed can b e accurately reco vered by maximizing ϕ 1 ( c s ). 16 0.7 0.8 0.9 1 1.1 1.2 1.3 0.4 0.6 0.8 1 theoretical result numerical simulation Figure 8 : Simulation and theoretical results for the w a ve s peed estimation Error Mean V ariance RMSE E r,est 3 . 795 × 10 − 2 1 . 517 × 10 − 4 3 . 988 × 10 − 2 E r,true − 9 . 776 × 10 − 4 3 . 262 × 10 − 12 9 . 776 × 10 − 4 E c,est 2 . 366 × 10 − 3 1 . 050 × 10 − 6 2 . 576 × 10 − 3 E c,true 2 . 241 × 10 − 3 7 . 607 × 10 − 7 2 . 403 × 10 − 3 T able 1: Mean, v ariance, and root mean square error (RMSE) of relative range and cross-range lo calization errors o ver 100 Monte Carlo trials. All errors are normalized b y the true target range. 5.1.2 Mon te Carlo exp erimen t W e conduct Mon te Carlo sim ulations to assess the statistical stability of the proposed w av e-sp eed estimator in a homogeneous background medium. W e consider 100 independent realizations of the noisy sources. F or eac h realization, 500 sampling points are dra wn according to the density function K ( x ), and the corresp onding cross-correlations are computed. The discretization grid for the wa ve data generalization is the same as that used in the previous exp erimen ts. F or wa ve-speed estimation, the searc hing grid for c s is chosen as 0 . 9:0 . 005:1 . 1, and the searc h domain is the same as that used in the previous exp erimen ts except for the refinement of the spatial grid in the region of in terest to 500 × 100 × 100. The statistical results are describ ed in the following. The estimated mean wa ve sp eed is 1 . 0379, the sample v ariance is 1 . 4151 × 10 − 4 , and the root mean square error (RMSE) is 0 . 0397. The total computational time for the Mon te Carlo sim ulations is appro ximately 20 hours. W e also ev aluate the relative lo calization errors for the reflector using the estimated wa ve sp eed and the true wa ve sp eed from each Monte Carlo trial. W e denote the estimator using estimated wa v e sp eed as b z r,est , and the one using the true w av e sp eed as b z r,true . The relative range and cross-range errors are defined as E r,est = ⟨ b z r,est , z r ⟩ | z r | 2 , E c,est = p | b z r,est − z r | 2 − E r,est | z r | , E r,true = ⟨ b z r,true , z r ⟩ | z r | 2 , E c,true = p | b z r,true − z r | 2 − E r,true | z r | . W e displa y the mean, sample v ariance, and RMSE of these errors in T able 1. W e observ e that the localization error in the range direction is impacted b y the uncertain ty in the estimated wa v e sp eed. The lo calization error in the cross-range direction, ho wev er, is not sensitive to the estimated w av e sp eed. 17 5.2 Random medium with ℓ c ≪ λ In this section,w e conduct numerical experiments to v alidate the metho dology for estimating the effectiv e w av e sp eed in a random medium, w e plot the profile of E |I ( z s ( c s , t g ) , c s ) | 2 , where the p oint spread function without the narrowband approximation is given by (4.6). The parameters are set as a 0 = 10, c 0 = 1, and t g = 50. W e c ho ose F 0 ( t ) = e − t 2 2 σ 2 cos(2 π ν 0 t ) , ν 0 = 10 , σ = 0 . 3 , for which the corresp onding bandwidth B H is close to 1. The resulting theoretical profile is sho wn in Figure 9, whose peak is attained at c s = c 0 = 1. 0.7 0.8 0.9 1 1.1 1.2 1.3 0 0.2 0.4 0.6 0.8 1 Figure 9 : Plot of the profile of E |I ( z s ( c s , t g ) , c s ) | 2 giv en in (4.6). The p eak is achiev ed at c s = c 0 = 1. 6 Conclusion In this paper, w e in vestigated passiv e imaging in the dayligh t configuration and developed a mathematical framew ork to quantitativ ely analyze the impact of wa ve-speed mismatc h in backpropagation-based reflector imaging. When the searching sp eed matches the true background sp eed, our results recov er the classical resolution analysis for p oin t-reflector imaging. In the presence of mismatch, we deriv ed explicit form ulas c haracterizing b oth the spatial shift of the reflector and the defo cusing of the asso ciated point spread function. This analysis naturally leads to a no vel estimator for the bac kground w av e speed. W e further extended the framework to random media and introduced the concept of a virtual guide star, whic h enables an effective and robust w a ve-speed estimator under dayligh t illumination. F or clarit y of exp osition and without loss of generality , our analysis fo cused on the case where the virtual guide star lies on the axis p erp endicular to the sensor arra y . W e expect that the theoretical framew ork dev elop ed here pro vides a foundation for further adv ances in w av e-sp eed estimation within the broader field of passive imaging. 18 App endix A Pro ofs of the results in Section 3 App endix A.1 Pro of of Theorem 3.1 T o b egin with, we observe that C (1) ( τ , x 1 , x 2 ) = C (1) ( − τ , x 2 , x 1 ) , w e write I ( z s , c s ) ≈ I i + ( z s , c s ) + I ii + ( z s , c s ) + I i − ( z s , c s ) + I ii − ( z s , c s ) , where I i ± ( z s , c s ) = N 2 2 |A| 2 Z Z A 2 d x 1 d x 2 C (1) i ( ± ( T s ( z s , x 2 ) + T s ( z s , x 1 )) , x 1 , x 2 ) , and I ii ± ( z s , c s ) = N 2 2 |A| 2 Z Z A 2 d x 1 d x 2 C (1) ii ( ± ( T s ( z s , x 2 ) + T s ( z s , x 1 )) , x 1 , x 2 ) . W e first focus on the calculation of I i + . W e ha ve C (1) i ( τ , x 1 , x 2 ) = σ r ℓ 3 r 2 7 π 4 c 2 0 ε 2 Z Z d y dω ω 2 b F ( ω ) K ( y ) | x 2 − z r || z r − y || x 1 − y | e i Φ i ( ω, y ) ε , where Φ i ( ω , y ) = ω ( T ( x 2 , y ) + T ( z r , y ) − T ( x 1 , y ) − τ ) . Therefore, w e deriv e that I i + ( z s , c s ) = σ r ℓ 3 r 2 8 π 4 c 2 0 ε 2 N 2 |A| 2 Z Z A 2 dσ ( x 1 ) σ ( x 2 ) Z Z d y dω ω 2 b F ( ω ) K ( y ) | x 2 − z r || z r − y || x 1 − y | e i Φ i + ( ω, y , x 1 , x 2 ) ε , (A.1) where Φ i + ( ω , y , x 1 , x 2 , z s ) = ω ( T ( x 2 , z r ) + T ( z r , y ) − T ( x 1 , y ) − T s ( x 1 , z s ) − T s ( x 2 , z s )) The analysis of the imaging function highly relies on the analysis of the term Φ i + , in the following, we divide the calculation of the (generalized) tra vel times into tw o parts. First, w e consider T ( x 2 , z r ) − T s ( x 2 , z s ). Let the searc hing point b e z s = z f + ε 1 2 r , where z f is the focusing p oin t to be determined. By T aylor expansion, we derive | z r − x 2 | = | z r | − ε 1 2 ⟨ z r , f x 2 ⟩ | z r | + O ( ε ) , | z s − x 2 | = | z s | − ε 1 2 ⟨ z s , f x 2 ⟩ | z s | + O ( ε ) = | z f | + ε 1 2 ⟨ z f , r ⟩ | z f | − ε 1 2 ⟨ z r , f x 2 ⟩ | z r | + O ( ε ) . T o av oid rapid phase terms in the integral (A.1), we derive • F or O (1) terms: | z r | c 0 = | z f | c s . • F or O ( ε 1 2 ) terms: ⟨ z f , r ⟩−⟨ z f , f x 2 ⟩ c s | z f | + ⟨ z r , f x 2 ⟩ c 0 | z r | = ⟨ z f , r ⟩ v 2 | z r | c 0 + ⟨ z r − z f v 2 , f x 2 ⟩ c 0 | z r | ∼ O ( ε 1 2 ) holds for all r ∈ R 3 , and f x 2 ∈ R 2 . 19 The t wo conditions give the follo wing equations for z f : | z f | = v | z r | , z ⊥ f = v 2 z ⊥ r , whose solution pro vides the p osition of z f : z f =  v 2 z ⊥ r , v q ( z q r ) 2 + (1 − v 2 ) | z ⊥ r | 2  , pro vided by v ≥ 1 − ( z q r ) 2 | z ⊥ r | 2 . In the follo wing calculation, we denote z f = ( z f , 1 , 0 , z f , 3 ) for simplicit y . Mean while, the tw o conditions ab o ve lead us to the following correct scale to proceed the analysis to the profile of imaging function. W e consider z s = z f + r s , where r s = ε 1 2 v 2 η 1 b g 1 + ε 1 2 v 2 η 2 b g 2 + εv η 3 b g 3 . W e now pro ceed to the detailed calculation of T ( x 2 , z r ) − T s ( x 2 , z s ). By T aylor’s expansion, we get | z r − x 2 | ≃ | z r | − ε 1 2 ⟨ z r , f x 2 ⟩ | z r | + ε | z r | 2 | f x 2 | 2 − ⟨ z r , f x 2 ⟩ 2 2 | z r | 3 c 0 , (A.2) and | z s − x 2 | ≃ | z s | − ε 1 2 ⟨ z s , f x 2 ⟩ | z s | + ε | z s | 2 | f x 2 | 2 − ⟨ z s , f x 2 ⟩ 2 2 | z s | 3 c 0 = v | z r | − ε 1 2 v ⟨ z r , f x 2 ⟩ | z r | + ε 2 v | z r |  − 2 v 2 ( η 1 ⟨ f x 2 , b g 1 ⟩ + η 2 ⟨ f x 2 , b g 2 ⟩ ) + | f x 2 | 2 − v 2 D f x 2 , b f 3 E + 2 v 2 | z r | η 3 + v 4 ( η 2 1 + η 2 2 )  (A.3) Hence, w e deriv e that T ( x 2 , z r ) − T s ( x 2 , z s ) ≃ − ε η 3 c 0 + ε | f x 2 | 2 2 | z r | c 0  1 − 1 v 2  + ε | z r | c 0  η 1 ⟨ f x 2 , b g 1 ⟩ + η 2 ⟨ f x 2 , b g 2 ⟩ − v 2 2 ( η 2 1 + η 2 2 )  = − ε η 3 c 0 + ε | f x 2 | 2 2 | z r | c 0  1 − 1 v 2  + ε | z r | c 0  g x 2 , 1 z f , 3 | z f | η 1 + g x 2 , 2 η 2 − v 2 2 ( η 2 1 + η 2 2 )  (A.4) No w, we fo cus on the term T ( z r , y ) − T ( x 1 , y ) − T s ( x 1 , z s ). F or the leading order, we ha v e T ( z r , y ) − T ( x 1 , y ) − T s ( x 1 , z s ) ≃ | z r − y | c 0 − | y | c 0 − | z f | c s = | z r − y | − | y | − | z r | c 0 . T o av oid the quic k phase term in the integral (A.1), we parametrize y as follo ws: y = | z r |  − s 3 b f 3 + ε 1 2 s 1 b f 1 + ε 1 2 s 2 b f 2  . The calculation of | z s − x 1 | is similar to (A.3), and w e calculate that | z r − y | ≃ | z r |  (1 + s 3 ) + ε s 2 1 + s 2 2 2(1 + s 3 )  , (A.5) 20 and | x 1 − y | ≃ | z r | s 3 + ε 1 2 D f x 1 , b f 3 E + ε    − D f x 1 , b f 3 E 2 2 | z r | s 3 + 1 2 | z r | s 3  | f x 1 | 2 + | z r | 2 ( s 2 1 + s 2 2 ) − 2 | z r |  s 1 D f x 1 , b f 1 E ) + s 2 D f x 1 , b f 2 E    (A.6) Then, w e deriv e that T ( z r , y ) − T ( x 1 , y ) − T s ( x 1 , z s ) ≃ − ε η 3 c 0 − ε v 2 ( η 2 1 + η 2 2 ) 2 | z r | c 0 − ε | z r | c 0 s 2 1 + s 2 2 2 s 3 (1 + s 3 ) + ε ⟨ f x 1 , z r ⟩ 2 2 | z r | 3 c 0  1 + 1 s 3  − ε | f x 1 | 2 2 | z r | c 0  1 s 3 + 1 v 2  + ε s 1 D f x 1 , b f 1 E ) + s 2 D f x 1 , b f 2 E c 0 s 3 + ε η 1 ⟨ f x 1 , b g 1 ⟩ + η 2 ⟨ f x 1 , b g 2 ⟩ | z r | c 0 = − ε η 3 c 0 − ε v 2 ( η 2 1 + η 2 2 ) 2 | z r | c 0 − ε | z r | c 0 s 2 1 + s 2 2 2 s 3 (1 + s 3 ) + ε s 1 c 0 s 3 z r, 3 | z r | g x 1 , 1 + ε s 2 c 0 s 3 g x 1 , 2 + ε g x 1 , 1 z f , 3 c 0 | z r || z f | η 1 + ε g x 1 , 2 | z r | c 0 η 2 + ε g x 1 , 1 2 z 2 r, 1 2 c 0 | z r | 3  1 + 1 s 3  − ε | f x 1 | 2 2 c 0 | z r |  1 s 3 + 1 v 2  (A.7) W e compute the in tegrals in s 1 and s 2 Z Z ds 1 ds 2 exp  i ω ε ( T ( z r , y ) − T ( x 1 , y ) − T s ( x 1 , z s ))  ≃ − 2 π ic 0 s 3 (1 + s 3 ) ω | z r | exp  i ω c 0  | f x 1 | 2 2 | z r |  1 − 1 v 2  − η 3 − v 2 ( η 2 1 + η 2 2 ) 2 | z r | c 0 + g x 1 , 1 z f , 3 | z f || z r | η 1 + g x 1 , 2 | z r | η 2  . It giv es Z Z ds 1 ds 2 Φ i + ( ω , y , x 1 , x 2 , z s ) ≃ − 2 π ic 0 s 3 (1 + s 3 ) ω | z r | × exp  i ω c 0  | f x 1 | 2 + | f x 2 | 2 2 | z r |  1 − 1 v 2  − 2 η 3 − v 2 ( η 2 1 + η 2 2 ) | z r | + z f , 3 η 1 | z f || z r | ( g x 1 , 1 + g x 2 , 1 ) + η 2 | z r | ( g x 1 , 2 + g x 2 , 2 )  . Therefore, w e obtain I i + ( z s , c s ) = σ r ℓ 3 r N 2 K ( 0 , z r ) 2 6 π 3 c 0 | z r | 2 ε Z B dω  − iω b F ( ω )  exp  − i ω c 0  2 η 3 + c 2 s ( η 2 1 + η 2 2 ) c 2 0 | z r |  × G 2 α r ,α f − a 0 ω c 0 | z r | η 1 , − a 0 ω c 0 | z r | η 2 , 0 , − a 2 0 ω c 0 | z r |  c 0 c s  2 − 1 !! Similar calculation can be done for I ii + , I i − , and I ii − . Only I ii − pro vides non-v anishing con tribution. Hence, the expression in (3.10) is deriv ed. App endix A.2 Pro of of Prop osition 3.2 W e revise the parametrization of the searching p oin t as z s = z f + ε 1 2 v 2 η 1 b g 1 + ε 1 2 v 2 η 2 b g 2 + v η 3 b g 3 . Compared to the previous one, the scaling in direction b g 3 is c hanged. W e conduct a similar calculation and get 21 | z s − x 2 | ≃ | z s | − ε 1 2 ⟨ z s , f x 2 ⟩ | z s | + ε | z s | 2 | f x 2 | 2 − ⟨ z s , f x 2 ⟩ 2 2 | z s | 3 c 0 = | z f | + v η 3 − ε 1 2 ⟨ f x 2 , b g 3 ⟩ + ε v 4 ( η 2 1 + η 2 2 ) 2( | z f | + v η 3 ) − v 2 ( η 1 ⟨ f x 2 , b g 1 ⟩ + η 2 ⟨ f x 2 , b g 2 ⟩ ) | z f | + v η 3 + | f x 2 | 2 − ⟨ f x 2 , b g 3 ⟩ 2 2( | z f | + v η 3 ) ! (A.8) Then, w e ha ve T ( x 2 , z r ) − T s ( x 2 , z s ) ≃ − η 3 c 0 + ε c 0    | f x 2 | 2 2 | z r | − | f x 2 | 2 2 v ( | z f + v η 3 | ) + ⟨ f x 2 , b g 3 ⟩ 2 2 v ( | z f + v η 3 | ) − D f x 2 , b f 3 E 2 2 | z r | + v ( η 1 ⟨ f x 2 , b g 1 ⟩ + η 2 ⟨ f x 2 , b g 2 ⟩ ) | z f | + v η 3 − v 3 ( η 2 1 + η 2 2 ) 2( | z f | + v η 3 )  = − η 3 c 0 + ε c 0 | f x 2 | 2 2  1 | z r | − 1 v 2 ( | z r | + η 3 )  + g x 2 , 1 2 z 2 r, 1 2  1 | z r | + η 3 − 1 | z r |  g x 2 , 1 | z r | + η 3 z f , 3 | z f | η 1 + g x 2 , 2 | z r | + η 3 η 2 − v 2 ( η 2 1 + η 2 2 ) 2( | z r | + η 3 )  , (A.9) and T ( z r , y ) − T ( x 1 , y ) − T s ( x 1 , z s ) ≃ − η 3 c 0 + ε c 0 − ( s 2 1 + s 2 2 ) | z r | 2 s 3 (1 + s 3 ) + g x 1 , 1 2 z 2 r, 1 2  1 s 3 | z r | + 1 | z r | + η 3  − | f x 1 | 2 2  1 | z r | s 3 + 1 v 2 ( | z r | + η 3 )  + s 1 s 3 z r, 3 z r g x 1 , 1 + s 2 s 3 g x 1 , 2 − v 2 ( η 2 1 + η 2 2 ) 2( | z r | + η 3 ) + z f , 3 | z f | g x 1 , 1 η 1 | z r | + η 3 + g x 1 , 2 η 2 | z r | + η 3 (A.10) Computing the in tegral in s 1 and s 2 , w e deriv e Z Z ds 1 ds 2 exp  i ω ε ( T ( z r , y ) − T ( x 1 , y ) − T s ( x 1 , z s ))  ≃ − 2 π ic 0 s 3 (1 + s 3 ) ω | z r | exp  i ω c 0  − 2 η 3 ε − v 2 ( η 2 1 + η 2 2 ) | z r | + η 3 + φ ( f x 1 ) + φ ( f x 2 )  , where, for x = ( x 1 , x 2 ) ∈ R 2 , φ ( x ) is defined by φ ( x ) = z f , 3 | z f | x 1 η 1 | z r | + η 3 + x 2 η 2 | z r | + η 3 + x 2 1 2 z 2 r, 3 | z r | 2 + x 2 2 2 ! η 3 | z r | ( | z r | + η 3 ) + | x | 2 2 1 − 1 v 2 | z r | + η 3 . (A.11) Therefore, w e ha ve I i + ( z s , c s ) = σ r ℓ 3 r N 2 K ( 0 , z r ) 2 6 π 3 c 0 | z r | 2 ε Z B dω  − iω b F ( ω )  exp  − i ω c 0  2 η 3 ε + c 2 s ( η 2 1 + η 2 2 ) c 2 0 ( | z r | + η 3 )  × G 2 α r ,α f − a 0 ω c 0 ( | z r | + η 3 ) η 1 , − a 0 ω c 0 ( | z r | + η 3 ) η 2 , − a 2 0 ω c 0 | z r | ( | z r | + η 3 ) η 3 , − a 2 0 ω c 0 ( | z r | + η 3 )  c 0 c s  2 − 1 !! Similar calculation can be done for I ii + , I i − , and I ii − . Only I ii − pro vides non-v anishing con tribution. Hence, w e hav e I ( z s , c s ) ≈ σ r ℓ 3 r N 2 K ( 0 , z r ) 2 6 π 3 c 0 | z r | 2 ε Re  Z B dω  − iω b F ( ω )  exp  − i ω c 0  2 η 3 ε + c 2 s ( η 2 1 + η 2 2 ) c 2 0 ( | z r | + η 3 )  ×G 2 α r ,α f − a 0 ω c 0 ( | z r | + η 3 ) η 1 , − a 0 ω c 0 ( | z r | + η 3 ) η 2 , − a 2 0 ω c 0 | z r | ( | z r | + η 3 ) η 3 , − a 2 0 ω c 0 ( | z r | + η 3 )  c 0 c s  2 − 1 !!) (A.12) 22 F urther, by the assumption that the bandwidth is of the order O ( ε ): b F ( ω ) = 1 εB H  b F 0  ω 0 − ω εB H  + b F 0  ω 0 + ω εB H  , w e derive the following expression I ( z s , c s ) ≈ σ r ℓ 3 r N 2 K ( 0 , z r ) ω 0 ε 2 4 π 2 c 0 | z r | 2 Re  − i exp  − i ω c 0  2 η 3 ε + c 2 s ( η 2 1 + η 2 2 ) c 2 0 ( | z r | + η 3 )  F 0  − 2 B 0 η 3 c 0  ×G 2 α r ,α f − a 0 ω c 0 ( | z r | + η 3 ) η 1 , − a 0 ω c 0 ( | z r | + η 3 ) η 2 , − a 2 0 ω c 0 | z r | ( | z r | + η 3 ) η 3 , − a 2 0 ω c 0 ( | z r | + η 3 )  c 0 c s  2 − 1 !!) (A.13) The condition that | r s | ≪ | z r | giv es the expression (3.18). When B 0 ≪ B c , we get the expression (3.22). When B 0 ≫ B c , w e get the expression (3.20). App endix B Pro ofs of the results in Section 4 App endix B.1 Pro of of Theorem 4.1 Similar to the pro of of Theorem 3.1, we observe that C (1) ( τ , x 1 , x 2 ) = C (1) ( − τ , x 2 , x 1 ) , w e write I ( z s , c s ) ≈ I i + ( z s , c s ) + I ii + ( z s , c s ) + I i − ( z s , c s ) + I ii − ( z s , c s ) , where I i ± ( z s , c s ) = N 2 2 |A| 2 Z Z A 2 d x 1 d x 2 C (1) i ( ± ( T s ( z s , x 2 ) + T s ( z s , x 1 )) , x 1 , x 2 ) , and I ii ± ( z s , c s ) = N 2 2 |A| 2 Z Z A 2 d x 1 d x 2 C (1) ii ( ± ( T s ( z s , x 2 ) + T s ( z s , x 1 )) , x 1 , x 2 ) . W e first focus on the calculation of I i + . W e ha ve C (1) i ( τ , x 1 , x 2 ) = 1 2 7 π 4 c 2 0 ε 2 Z D d z ρ µ ( z ) Z Z d y dω ω 2 b F ( ω ) K ( y ) | x 2 − z || z − y || x 1 − y | e i Φ i ( ω, y , z ) ε , where Φ i ( ω , z , y ) = ω ( T ( x 2 , y ) + T ( z , y ) − T ( x 1 , y ) − τ ) . Therefore, w e deriv e that I i + ( z s , c s ) = 1 2 8 π 4 c 2 0 ε 2 N 2 |A| 2 Z Z A 2 dσ ( x 1 ) σ ( x 2 ) Z D ρ µ ( z ) d z Z Z d y dω ω 2 b F ( ω ) K ( y ) e i Φ i + ( ω, y , x 1 , x 2 , z s , z ) ε | x 2 − z r || z r − y || x 1 − y | , (B.1) where Φ i + ( ω , y , x 1 , x 2 , z s , z ) = ω ( T ( x 2 , z ) + T ( z , y ) − T ( x 1 , y ) − T s ( x 1 , z s ) − T s ( x 2 , z s )) The analysis of the imaging function highly relies on the analysis of the term Φ i + , in the following, we divide the calculation of the (generalized) tra vel times into tw o parts. 23 First, w e consider T ( x 2 , z ) − T s ( x 2 , z s ): By T aylor expansion, w e hav e | x 2 − z s | ≃ | z s | + | f x 2 | 2 2 | z s | ε, (B.2) and | x 2 − z | ≃ | z g | + ε 2 | z g |  η 2 1 + η 2 2 + 2 η 3 | z g | + | f x 2 | 2 − 2( η 1 ⟨ b e 1 , f x 2 ⟩ + η 2 ⟨ c e 2 , f x 2 ⟩ )  (B.3) Then, w e ha ve T ( x 2 , z ) − T s ( x 2 , z s ) ≃ ε c 0  η 2 1 + η 2 2 2 | z g | + η 3 + | f x 2 | 2 2 | z g |  1 − 1 v 2  − 1 | z g | ( η 1 ⟨ b e 1 , f x 2 ⟩ + η 2 ⟨ c e 2 , f x 2 ⟩ ))  (B.4) T o calculate the term T ( z , y ) − T ( x 1 , y ) − T s ( x 1 , z s ), w e in tro duce the parametrization of y as: y = | z g |  − s 3 b e 3 + ε 1 2 s 1 b e 1 + ε 1 2 s 2 b e 2  . By T aylor expansion, w e hav e | x 1 − z s | ≃ | z s | + | f x 1 | 2 2 | z s | ε, (B.5) | z − y | ≃ (1 + s 3 ) | z g | + ε 2(1 + s 3 ) | z g |  2 η 3 (1 + s 3 | z g | ) + ( η 1 − z g s 1 ) 2 + ( η 2 − z g s 2 ) 2  , (B.6) | x 1 − y | ≃ | z g | s 3 + ε  | f x 1 | 2 2 | z g | s 3 + ( s 2 1 + s 2 2 ) | z g | 2 s 3 − s 1 ⟨ f x 1 , b e 1 ⟩ + s 2 ⟨ f x 1 , b e 2 ⟩ s 3  . (B.7) Hence, w e ha ve T ( z , y ) − T ( x 1 , y ) − T s ( x 1 , z s ) ≃ ε c 0  − | f x 1 | 2 2 | z g |  1 v 2 + 1 s 3  − ( s 2 1 + s 2 2 ) | z g | 2 s 3 + s 1 ⟨ f x 1 , b e 1 ⟩ + s 2 ⟨ f x 1 , b e 2 ⟩ s 3 + η 3 + 1 2(1 + s 3 ) | z g |  | z g | 2 ( s 2 1 + s 2 2 ) − 2 | z g | ( s 1 η 1 + s 2 η 2 ) + ( η 2 1 + η 2 2 )   . (B.8) Computing the in tegral in s 1 and s 2 , w e deriv e Z Z ds 1 ds 2 exp  i ω ε ( T ( z , y ) − T ( x 1 , y ) − T s ( x 1 , z s ))  ≃ − 2 π ic 0 s 3 (1 + s 3 ) ω | z g | exp  i ω c 0  2 η 3 + η 2 1 + η 2 2 | z g | + ψ ( f x 1 ) + ψ ( f x 2 )  , where, for x = ( x 1 , x 2 ) ∈ R 2 , φ ( x ) is defined by ψ ( x ) = | x | 2 2 | z g |  1 − 1 v 2  − η 1 x 1 + η 2 x 2 | z g | . (B.9) Therefore, w e obtain I i + ( z s , c s ) = N 2 K ( 0 , z g ) ε 2 6 π 3 c 0 | z g | 2 Z ρ µ ( z g ( t g ) + r ( η )) Z B dω  − iω b F ( ω )  exp  i ω c 0  2 η 3 + ( η 2 1 + η 2 2 ) | z g |  × G 2 1 , 1 a 0 ω c 0 | z g | η 1 , a 0 ω c 0 | z g | η 2 , 0 , − a 2 0 ω c 0 | z g |  c 0 c s  2 − 1 !! d η . Similar calculation can be done for I ii + , I i − , and I ii − . Only I ii − pro vides non-v anishing con tribution. Hence, w e hav e I ( z s , c s ) = N 2 K ( 0 , z g ) ε 2 6 π 3 c 0 | z g | 2 Z ρ µ ( z g ( t g ) + r ( η )) Re  Z B dω  − iω b F ( ω )  exp  i ω c 0  2 η 3 + ( η 2 1 + η 2 2 ) | z g |  ×G 2 1 , 1 a 0 ω c 0 | z g | η 1 , a 0 ω c 0 | z g | η 2 , 0 , − a 2 0 ω c 0 | z g |  c 0 c s  2 − 1 !!) d η . 24 App endix B.2 Pro of of Prop osition 4.3 W e observe that E |I ( z s , c s ) | 2 ∼ Z D Z D E ( ρ µ ( z ) ρ µ ( z ′ )) · P ( z , z s , c s ) P ( z ′ , z s , c s ) d z d z ′ = Z D Z D Σ  z − z ′ ℓ c  · P ( z , z s , c s ) P ( z ′ , z s , c s ) d z d z ′ = ℓ 3 c · Z D Z D Σ ( x ) · P ( z ′ + ℓ c x , z s , c s ) P ( z ′ , z s , c s ) d x d z ′ ≃ l 3 ∥ Σ ∥ L 1 ( D ) · Z R 3 |P ( z , z s , c s ) | 2 d z , where w e use Lemma C.1 in the last inequalit y . T o calculate the profile of the p oin t spread function for B H ≪ ω 0 , w e first denote ϕ ( η ) = 1 c 0  2 η 3 + η 2 1 + η 2 | z g |  , G ( η ) = G 1 , 1 a 0 ω 0 c 0 | z g | η 1 , a 0 ω 0 c 0 | z g | η 2 , 0 , − a 2 0 ω 0 c 0 | z g |  c 0 c s  2 − 1 !! , and f ( η ) = − i exp ( iω 0 ϕ ( η )) F 0 ( B H ϕ ( η )) · G 2 ( η ) Then, w e notice that Z R 3 |P µ ( t g , η , c s ) | 2 d η = 1 2 Z R 3 | f ( η ) | 2 d η + 1 2 Re  Z R 3 f 2 ( η ) d η  ≜ I 1 + I 2 . (B.10) F or I 1 , w e ha ve I 1 = 1 2 Z R 3 | F 0 ( B H ϕ ( η )) | 2 |G ( η ) | 4 d η . F or I 2 , w e notice that ∂ ∂ η 3 ϕ ( η ) = 2 c 0 . (B.11) F or any fixed ( η 1 , η 2 ), w e let ζ = B H ϕ ( η ) and implement the change of v ariable. W e calculate that η 3 := Ψ( η 1 , η 2 , ζ ) = c 0 2  ζ B H − η 2 1 + η 2 2 | z g |  , (B.12) and the Jacobian is ∂ η 3 ∂ ζ = c 0 2 B H . Then, w e deriv e that Z R 3 f 2 ( η ) d η = − c 0 2 B H Z R 2 dη 1 dη 2 Z R e i 2 ω 0 B H ζ F 2 H, 0 ( ζ ) G 4 ( η 1 , η 2 , ζ ) dζ Due to the existence of the quic k phase term, the integral is negligible. Therefore, w e ha ve Z R 3 |P µ ( t g , η , c s ) | 2 d η ≃ 1 2 Z R 3     F 0  B H c 0  2 η 3 + η 2 1 + η 2 2 | z g |      2      G 1 , 1 a 0 ω 0 c 0 | z g | η 1 , a 0 ω 0 c 0 | z g | η 2 , 0 , − a 2 0 ω 0 c 0 | z g |  c 0 c s  2 − 1 !!      4 dη 1 dη 2 dη 3 . 25 App endix B.3 Pro of of Prop osition 4.5 Here, w e provide a pro of for the case ρ µ is a stationary Gaussian pro cess. F or general stationary pro cess, the proof can b e done by a v ariance calculation. W e notice that the cov ariance function of the stationary Gaussian pro cess ρ µ is L 1 , then ρ µ is ergodic. Therefore |I ( z s , c s ) | 2 is also ergo dic as a measurable function of ρ µ . Combined with the stationarit y shown in the previous prop osition, b y Birkhoff ’s Theorem, w e deriv e the result. App endix B.4 Pro of of Corollary 4.6 Here, w e offer a qualitativ e argument. W e notice that for fixed c s , the essen tial support G 1 , 1 concen trates around η = 0, and also the essential supp ort of F 0 concen trate around 0. The quantit y E |I ( z s , c s ) | 2 b eha v es qualitativ ely as the function     G 1 , 1  0 , 0 , 0 , − a 2 0 ω 0 c 0 | z g |   c 0 c s  2 − 1      4 up to some multiplicit y constant. Then, it is clear that the maximum is achiev ed at c s = c 0 . App endix C Auxiliary lemmas Lemma C.1. F or any deterministic function f ∈ C 1 ( D ) , g ∈ C 0 ( D ) , Z D × D Co v ( n ε ( x ) , n ε ( y )) f ( x ) g ( y ) d x d y =  Z R d ε d Σ( z ) d z  · Z D f ( x ) g ( x ) d x + O  ε d +1  Pr o of. By definition and b y a c hange of v ariable Z D × D Co v ( n ε ( x ) , n ε ( y )) f ( x ) g ( y ) d x d y = Z D × D Σ  x − y ε  f ( x ) g ( y ) d x d y = ε d Z D Z D ε Σ( z ) f ( ε z + y ) d z ! ¯ g ( y ) d y Therefore,      Z D × D Co v ( n ε ( x ) , n ε ( y )) f ( x ) g ( y ) d x d y − ε d Z D ε Σ( z ) d z ! Z D f ( x ) g ( x ) d x      ≲ ε d Z D ε  Z D | g ( y ) || f ( ε z + y ) − f ( y ) | d y  | Σ( z ) | d z ≲ ε d +1 Z D ε | Σ( z ) || z | d z !  Z D |∇ f ( y ) || g ( y ) | d y  (C.1) Since Σ has the deca ying rate | Σ( x ) | ≲ (1 + | x | ) − 4 − η , for some η > 0, the right-hand side of the difference con verges, which leads to the result. Lemma C.2. We pr esent an inte gr al which is fr e quently use d in the c alculation of the p oint spr e ad function: Z R e i ( − Ax 2 + B x ) dx = r π A e − i π 4 e i B 2 4 A . (C.2) Pr o of. Z R e i ( − Ax 2 + B x ) dx = Z R e − iA ( x − B 2 A ) 2 · e i B 2 4 A dx = r π A e − i π 4 e i B 2 4 A , (C.3) where w e use the iden tit y that Z R e − i s 2 2 ds = √ 2 π e − i π 4 . 26 References [1] R. L. W eav er and O. I. Lobkis, “Ultrasonics without a source: Thermal fluctuation correlations at mhz frequencies,” Physic al R eview L etters , v ol. 87, no. 13, p. 134 301, 2001. [2] T. L. Duv all Jr, S. Jeffferies, J. Harv ey, and M. Pomeran tz, “Time–distance helioseismology,” Natur e , v ol. 362, no. 6419, pp. 430–432, 1993. [3] K. G. Sabra, P . Roux, P . Gerstoft, W. Kuperman, and M. C. F ehler, “Extracting coherent co da arriv als from cross-correlations of long p eriod seismic wa ves during the mount st. helens 2004 eruption,” Ge ophysic al r ese ar ch letters , vol. 33, no. 6, 2006. [4] A. Curtis, P . Gerstoft, H. Sato, R. Snieder, and K. W ap enaar, “Seismic in terferometry—turning noise in to signal,” The L e ading Edge , vol. 25, no. 9, pp. 1082–1092, 2006. [5] E. Larose, L. Margerin, A. Dero de, et al. , “Correlation of random wa v efields: An interdisciplinary review,” Ge ophysics , v ol. 71, no. 4, SI11–SI21, 2006. [6] J. Garnier and G. P apanicolaou, Passive Imaging with A mbient Noise . Cambridge Univ ersity Press, 2016. [7] O. I. Lobkis and R. L. W eav er, “On the emergence of the green’s function in the correlations of a diffuse field,” The Journal of the A c oustic al So ciety of A meric a , v ol. 110, no. 6, pp. 3011–3017, 2001. [8] R. Snieder, “Extracting the green’s function from the correlation of co da wa ves: A deriv ation based on stationary phase,” Phys. R ev. E , vol. 69, p. 046 610, 4 Apr. 2004. doi : 10.1103/PhysRevE.69.046610 . [Online]. Av ailable: https://link.aps.org/doi/10.1103/PhysRevE.69.046610 . [9] Y. C. De V erdi ` ere, “Semiclassical analysis and passiv e imaging,” Nonline arity , v ol. 22, no. 6, R45, 2009. [10] C. Bardos, J. Garnier, and G. Papanicolaou, “Iden tification of green’s functions singularities b y cross correlation of noisy signals,” Inverse Pr oblems , vol. 24, no. 1, p. 015 011, Jan. 2008. doi : 10.1088/0266- 5611/24/1/015011 . [Online]. Av ailable: https://doi.org/10.1088/0266- 5611/24/1/015011 . [11] J. Garnier and G. P apanicolaou, “Passiv e sensor imaging using cross correlations of noisy signals in a scattering medium,” SIAM Journal on Imaging Scienc es , vol. 2, no. 2, pp. 396–437, 2009. doi : 10 . 1137 / 080723454 . eprin t: https : / / doi . org / 10 . 1137 / 080723454 . [Online]. Av ailable: https : //doi.org/10.1137/080723454 . [12] T. Helin, M. Lassas, L. Oksanen, and T. Saksala, “Correlation based passive imaging with a white noise source,” Journal de Math´ ematiques Pur es et Appliqu´ ees , vol. 116, pp. 132–160, 2018, issn : 0021- 7824. doi : https : / / doi . org / 10 . 1016 / j . matpur . 2018 . 05 . 001 . [Online]. Av ailable: https : //www.sciencedirect.com/science/article/pii/S0021782418300710 . [13] A. Agaltsov, T. Hohage, and R. Novik ov, “Global uniqueness in a passive inv erse problem of helioseis- mology ,” Inverse Pr oblems , vol. 36, no. 5, p. 055 004, 2020. [14] A. Bakulin and R. Calvert, “The virtual source metho d: Theory and case study,” Ge ophysics , v ol. 71, no. 4, SI139–SI150, 2006. [15] K. W ap enaar and J. Thorbeck e, “Virtual sources and their responses, part i: Time-reversal acoustics and seismic in terferometry,” Ge ophysic al Pr osp e cting , v ol. 65, no. 6, pp. 1411–1429, 2017. [16] K. W ap enaar, J. Thorb ec ke, J. v an der Neut, E. Slob, and R. Snieder, “Review paper: Virtual sources and their resp onses, part ii: Data-driven single-sided fo cusing,” Ge ophysic al Pr osp e cting , vol. 65, no. 6, pp. 1430–1451, 2017. doi : https : / / doi . org / 10 . 1111 / 1365 - 2478 . 12495 . eprin t: https : / / onlinelibrary . wiley . com / doi / pdf / 10 . 1111 / 1365 - 2478 . 12495 . [Online]. Av ailable: https : //onlinelibrary.wiley.com/doi/abs/10.1111/1365- 2478.12495 . [17] J. Garnier and G. Papanicolaou, “Role of scattering in virtual source array imaging,” SIAM Journal on Imaging Scienc es , v ol. 7, no. 2, pp. 1210–1236, 2014. [18] P . Gou ´ edard, L. Stehly, F. Brenguier, et al. , “Cross-correlation of random fields: Mathematical ap- proac h and applications,” Ge ophysic al Pr osp e cting , vol. 56, no. 3, pp. 375–393, 2008. doi : https : / / doi . org / 10 . 1111 / j . 1365 - 2478 . 2007 . 00684 . x . eprint: https : / / onlinelibrary . wiley . com / doi / pdf / 10 . 1111 / j . 1365 - 2478 . 2007 . 00684 . x . [Online]. Av ailable: https : / / onlinelibrary . wiley.com/doi/abs/10.1111/j.1365- 2478.2007.00684.x . 27 [19] U. Harmank ay a, A. Kaslilar, J. Thorb ec ke, K. W ap enaar, and D. Draganov, “Lo cating near-surface scatterers using non-ph ysical scattered wa ves resulting from seismic in terferometry,” Journal of Applie d Ge ophysics , vol. 91, pp. 66–81, 2013. [20] A. Kaslilar, U. Harmank ay a, K. W ap enaar, and D. Draganov, “Estimating the location of a tunnel using correlation and in version of ra yleigh wa ve scattering,” Ge ophysic al R ese ar ch L etters , vol. 40, no. 23, pp. 6084–6088, 2013. [21] L. A. Konstantaki, D. Dragano v, T. Heimo v aara, and R. Ghose, “Imaging scatterers in landfills using seismic in terferometry,” Ge ophysics , vol. 78, no. 6, EN107–EN116, 2013. [22] J. Garnier and G. Papanicolaou, “Resolution analysis for imaging with noise,” Inverse Pr oblems , vol. 26, no. 7, p. 074 001, 2010. [23] J. Garnier and G. P apanicolaou, “Correlation-based virtual source imaging in strongly scattering ran- dom media,” Inverse Pr oblems , vol. 28, no. 7, p. 075 002, Jun. 2012. doi : 10.1088/0266- 5611/ 28 / 7/ 075002 . [Online]. Av ailable: https://doi.org/10.1088/0266- 5611/28/7/075002 . [24] N. M. Shapiro, M. Campillo, L. Stehly, and M. H. Ritzw oller, “High-resolution surface-w av e tomogra- ph y from am bient seismic noise,” Scienc e , v ol. 307, no. 5715, pp. 1615–1618, 2005. [25] F.-C. Lin, M. H. Ritzwoller, and R. Snieder, “Eikonal tomograph y: Surface wa v e tomography b y phase fron t tracking across a regional broad-band seismic arra y,” Ge ophysic al Journal International , vol. 177, no. 3, pp. 1091–1110, 2009. [26] M. H. Ritzw oller, F.-C. Lin, and W. Shen, “Am bient noise tomograph y with a large seismic array,” Comptes R endus Ge oscienc e , vol. 343, no. 8-9, pp. 558–570, 2011. [27] K. Sager, L. Ermert, C. Bo ehm, and A. Fic htner, “T o wards full w av eform ambien t noise in version,” Ge ophysic al Journal International , v ol. 212, no. 1, pp. 566–590, 2018. [28] A. Fich tner and V. C. Tsai, “Theoretical foundations of noise in terferometry,” in Seismic A mbient Noise , N. Nak ata, L. Gualtieri, and A. Fich tner, Eds., Cambridge Universit y Press, 2019, ch. 4, pp. 109– 143. [29] B. M ¨ uller, T. Hohage, D. F ournier, and L. Gizon, “Quan titative passive imaging by iterative hologra- ph y: The example of helios eismic holography ,” Inverse Pr oblems , v ol. 40, no. 4, p. 045 016, 2024. [30] J. Garnier and K. Sølna, “Bac kground v elo cit y estimation with cross correlations of incoherent wa v es in the parab olic scaling,” Inverse Pr oblems , v ol. 25, no. 4, p. 045 005, F eb. 2009. doi : 10.1088/0266- 5611/25/4/045005 . [Online]. Av ailable: https://doi.org/10.1088/0266- 5611/25/4/045005 . [31] W. Lambert, J. Robin, L. A. Cobus, M. Fink, and A. Aubry, “Ultrasound matrix imaging—part i: The fo cused reflection matrix, the f-factor and the role of m ultiple scattering,” IEEE T r ansactions on Me dic al Imaging , vol. 41, no. 12, pp. 3907–3920, 2022. [32] W. Lam b ert, L. A. Cobus, J. Robin, M. Fink, and A. Aubry, “Ultrasound matrix imaging—part ii: The distortion matrix for ab erration correction ov er multiple isoplanatic patches,” IEEE T r ansactions on Me dic al Imaging , vol. 41, no. 12, pp. 3921–3938, 2022. doi : 10.1109/TMI.2022.3199483 . [33] F. Bureau, E. Giraudat, A. L. Ber, et al. , “Reflection matrix imaging for w av e velocity tomograph y,” arXiv pr eprint arXiv:2409.13901 , 2024. [34] J. Garnier, L. Giov angigli, Q. Go epfert, and P . Millien, “Probing the sp eckle to estimate the effectiv e sp eed of sound, a first step to wards quantitativ e ultrasound imaging,” Inverse Pr oblems and Imaging , 2025, issn : 1930-8337. doi : 10.3934 /ipi .2026001 . [Online]. Av ailable: https :// www. aimsciences. org/article/id/692581ef1511553245b010e2 . [35] E. Giraudat, A. Burtin, A. Le Ber, M. Fink, J.-C. Komorowski, and A. Aubry, “Matrix imaging as a to ol for high-resolution monitoring of deep v olcanic plumbing systems with seismic noise,” Communic ations Earth & Envir onment , v ol. 5, no. 1, p. 509, 2024. 28