Acoustic Lucky Imaging for microphone phased arrays

. Introduction Measuring the strength and location of sound sources is crucial to the development of silent technologies. Acoustic imaging echniques allow for visualization and quantification of sound sources in wind tunnel experiments. An array of microphones samples he sound waves emitted from the sound sources, similar to a lens collecting light waves. The microphone signals are focused igitally to obtain the acoustic image. This focusing method is knows as the beamforming method [1,2]. A requisite for beamforming s sufficient coherence between microphone signals. Considering one sound source, the signals measured by different microphones hould be identical up to a time shift and amplitude scaling. However, when the acoustic signal propagates through a shear layer the urbulence continuously alters the acoustic propagation. This leads to random distortion of the microphones signals and consequently he coherence between two microphone signals is reduced. Coherence loss between microphone signals is undesired because it reduces the peak value of the Sound Pressure Level (SPL) n the acoustic image [3]. This hampers the accurate evaluation of the acoustic source strength. Furthermore, the resolution of the coustic image is reduced, making it difficult to separate individual sound sources. Oerlemans and Sijtsma [4] investigated acoustic easurements performed on a scaled aircraft model placed in the open jet of a large industrial wind tunnel. The acoustic waves raveled through a thick shear layer to the microphone array placed outside the jet. For frequencies ranging from 8 kHz to 12 kHz difference of 10 – 15 dB was observed between the peak SPL in the acoustic image and the measurement of a single microphone ∗ Corresponding author at: University of Twente, Department of Thermal and Fluid Engineering, Engineering Fluid Dynamics research group, the Netherlands. E-mail address: j.biesheuvel@utwente.nl (J. Biesheuvel). vailable online 19 October 2022 022-460X/© 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license http://creativecommons.org/licenses/by/4.0/). ttps://doi.org/10.1016/j.jsv.2022.117357 eceived 5 April 2022; Received in revised form 10 August 2022; Accepted 3 October 2022 Journal of Sound and Vibration 544 (2023) 117357 J. Biesheuvel et al. f f a t c a t m p b d r m c t s h h a a t c s a e m t l w A u w d s t f o a l t l d l t a d s m i t s w a b b p e at 50 m⋅s−1 and an estimated shear layer thickness of ∼1 m. The results also indicate that spatially integrated source SPL compared avorably with the SPL measured by a single microphone, the difference being less than 2 dB. This suggests that the largest error is ound in the peak level and not in the integrated level. There are several ways to increase the signal coherence and consequently the coustic image quality. The signal coherence is improved by removing the shear layer, e.g., using a closed test section, or a hybrid est section [5]. However, this is not always a possible or desirable solution. Improved post-processing techniques, able to correct for oherence loss, were therefore investigated. Several techniques have been developed for use in aeroacoustic measurements, nautical pplications, or astronomy. In the field of astronomy images of stars or planets are blurred by atmospheric turbulence influencing he propagation of light, analogue to the blurring of acoustic images by shear layer turbulence. Research regarding the loss of coherence between two microphone signals is closely related to the distortion of the signals easured by the individual microphones, commonly used to measure the acoustic spectrum. The turbulence induced fluctuations in ropagation time cause a tonal peak in the spectrum to spread over a wider frequency range. This phenomenon is known as spectral roadening, or more colloquially as ‘‘hay-stacking’’. A fundamental theory regarding the scattering of sound due to turbulence was erived by Howe [6]. This theory included the absorption, scattering and spectral broadening of the acoustic spectrum. The results egarding spectral broadening and scattering agreed with predictions based on geometrical acoustics. Aside from the spectrum easured by a single microphone additional information regarding the phase of the sound waves is obtained from the crossorrelation between two microphone signals. Yang [7] used the theory of geometrical acoustics combined with a path integral o compute the coherence between two rays (signals) moving through a medium disturbed by linear waves. In a more recent tudy McAlpine and Tester [8] derived an analytical model to describe the two-point correlation for sound propagation through omogeneous axisymmetric turbulence. This model showed better agreement with measurements compared to models based on omogeneous isotropic turbulence. The signal measured by the microphone equals the signal emitted by the sound source shifted by the acoustic delay time. This coustic time delay can be averaged to obtain a time-invariant delay caused by the finite propagation speed of the sound and dvection by the mean flow. The time-variant delay, i.e., the deviation from the average, is due to the turbulent shear layer. When he shear layer is thin compared to the wavelength the time-invariant delay time is predicted by the method of Amiet [9,10]. In this ase the influence of the turbulence can be neglected. However, when the shear layer is not thin compared to the wavelength and the hear layer is strongly turbulent the amplitude of the time-variant delay increases. Freund and Fleischman [11] compared ray-traced coustics with and without turbulence. Their most notable finding was that the evolving turbulence allowed the sound rays to more asily escape the potential core when traveling upstream. Koop et al. [12] investigated the time fluctuations as experienced by the icrophone array. They performed numerical as well as experimental studies to predict and measure the phase variations due to urbulence interacting with the sound waves. The use of a known source allowed the time fluctuations to be measured. Coherence oss in acoustic measurements was also observed by Dougherty [13] in the form of an apparent increase of the sound power level hen switching from a larger diameter array optimized for low frequencies to a small diameter array designed for higher frequencies. model for wave-coherence based on the work of Tatarskii [14] was used to predict the coherence loss. The model relies on the se of a statistical description of the turbulent fluctuations to assess the decorrelation of the sound waves. A structure function as used to determine the correlation between the turbulence and the fluctuations in acoustic propagation time. Blacodon [15] erived a correction method using a known reference source to derive a deconvolution operation to remove the smearing of the ources in the acoustic map. The author compared the integrated spectrum of an open test section measurement, computed using he deconvolution operation, with the spectrum of a single microphone measurement in a closed test section. The results compared avorably up to 10 kHz. Tatarskii’s theory was also used by Ernst et al. [16] to predict the coherence loss due to shear layers in pen jet wind tunnels. With the model of Tatarskii Ernst et al. predicted the magnitude of the coherence loss. They conclude that n impractical measurement time is necessary for statistical convergence of the coherence estimator in the case of severe coherence oss. Pires et al. [17] used Tatarskii’s model to study the coherence loss caused by the boundary layer of a closed test section wind unnel. Both Ernst et al. and Pires et al. performed flow measurements and compared the experimentally obtained acoustic coherence oss to predictions based on experimentally measured structure functions. The model captured well the coherence loss for larger istances between microphone pairs. This corresponds to the limit where the turbulence encountered over the acoustic paths is no onger correlated. Biesheuvel et al. [18] used Tatarskii’s model to assess the effect of coherence loss on beamforming results due o thick boundary layers (∼0.2 m) as encountered in industrial scale wind tunnels. For frequencies between 20 kHz and 50 kHz decrease in beamforming SPL was predicted ranging from 2 dB to 5 dB. In experiments with models the frequencies of interest uring the measurement scale inversely proportional to the model scale. Considering a 1:10 model scale, the frequencies on the full cale range from 2 kHz to 5 kHz, which is well within the audible range. A simple method to increase overall microphone coherence is to discard the incoherent microphone signals. In general, icrophone weighting methods are known as shading techniques [19]. In practice the microphones near the outer rim are discarded n the processing, effectively scaling down the microphone array dimension. Microphones that are spatially close are more likely o measure similar wavefront distortions and retain their signal coherence. Severe loss of signal coherence may even limit the tatistical convergence of the test results [16]. Shading based on the coherence of microphone pairs was used by Amaral et al. [20], ho defined the weighting factor of a microphone as the mean coherence between the microphone and all other microphones of the rray. This weighting improved the results of a deconvolution technique [21]. Bahr and Lockard [22] developed shading weights ased on the distance between the microphones and the frequency. They note a visual improvement for beamforming maps, judged y the absence of blurred sources. In general shading can recover some of the lost resolution by effectively removing microphone airs that, due to the coherence loss, have a cross-correlation that is not converged. The downside of applying shading is that the 2 ffective diameter of the microphone array becomes smaller, which also decreases the acoustic image resolution. Koop et al. [12] Journal of Sound and Vibration 544 (2023) 117357 J. Biesheuvel et al. h k t r b m t i m s r t a t r d t p r s o r i ‘ s o S o o T d h o f t s d r Fig. 1. Overview of the Acoustic Lucky Imaging methodology. int at the possibility to use a known source as a ‘‘guide star’’. This was done by Sijtsma [23] who used the presence of a few nown sources to estimate the wavefront disturbance and correct the phase of the microphone signals. A similar algorithm is used o correct images of stars using known stars or lasers, hence the technique is named ‘‘guide star’’. The propagation time from the est of the wind tunnel model could be estimated through interpolation between speaker locations. This method improved the eamforming results near the reference speakers. The technique required a large computational effort. A time domain beamforming ethod is presented by Cousson et al. [24]. Although, their intention was to present a method to identify moving acoustic sources he method implicitly incorporates a non-stationary delay time. When the acoustic source is allowed to move over time, this also mplies that the acoustic propagation is time-dependent. A localization method for acoustic sources in random moving media by eans of beamforming was presented by Gay et al. [25]. The authors used Kirchhoff migration (back propagation of the measured ignals) and coherent interferometry (back propagation of the signal cross-correlation [26]) to localize acoustic sources in a jet. The esults were favorable, given that the turbulence induced perturbations were not too large. In the field of nautical research acoustic measurements are degraded by ocean turbulence. Unknown fluctuations in propagation ime degrade results obtained for passive ranging. Ge and Kirsteins [27,28] developed a method to classify acoustic snapshots s either coherent or incoherent. Before processing, the data was ranked and the highly coherent snapshots were accumulated, aking advantage of the stochastic nature of the distortions. Ge and Kirsteins then demonstrated a reduction of smearing effects in ange-bearing plots. Looking further, into the field of astronomy, similar issues are encountered. Loss of resolution occurs when temperature and ensity variations due to atmospheric turbulence alter the local refraction coefficient of the air, resulting in a distorted image of he stars. One algorithm used to correct for this is called the ‘‘Luck exposure’’ [29] or ‘‘Lucky Imaging’’ [30] algorithm. For a articular section of the sky a large number of images are made using a very short exposure time. Only the images with good esolution are further processed. In this manner, the algorithm exploits the stochastic nature of the turbulence, which implies that ome of the images acquired are less disturbed by the turbulence than others. When the observation time is long enough a subset f ‘‘good quality’’ images with high resolution is selected for further processing. Dantowitz [31] used the method to create highesolution images of the planet Mercury. Law [32] also used to method to increase the resolution of an image of the stars, the ncreased resolution made more faint objects visible. Baldwin [33] showed an improvement in image quality by retaining only the ‘best’’ 1% exposures. In a later study Baldwin [34] performed numerical simulations using the random phase screen method. In this tudy a Kolmogorov spectrum was used to simulate the turbulence. The results compared favorably with results obtained from an bservatory. The most notable efforts to explicitly correct acoustic beamforming maps for turbulence induced fluctuations have been made by ijtsma, Blacodon and Gay et al. [15,23,25]. In addition to these efforts we propose a novel method to correct acoustic images based n the ‘‘Lucky Imaging’’ technique. This technique offers various advantages: (1) it does not require the use of ‘‘guiding’’ sources n the model or reference acoustic sources before or after the test. (2) The method is relatively simple and easy to implement. (3) he ‘‘Lucky imaging’’ methodology relies on filtering and not on the ensemble average of the microphone cross-correlations. This ifference is mainly relevant when performing acoustic tests in large industrial wind tunnels which have a thick shear layer or at igh frequencies. The large coherence loss may then lead to excessive measurement times. In this paper the Lucky Imaging methodology is developed to be used in acoustic wind tunnel tests [35]. Fig. 1 shows a schematic verview of the (Acoustic) Lucky Imaging method. The method consists of three stages. In the first stage acoustic images are obtained rom measurement intervals short enough to ‘‘freeze’’ the turbulence induced distortions. In the second stage those acoustic images hat are severely distorted are removed from the processing. The third, and final, stage combines all the ‘‘good’’ images into a ingle acoustic image. To support the correction methodology a model for coherence loss is derived. This model is based on the iscretization of the wavefront distortions. The model is used to show why and when Acoustic Lucky Imaging improves beamforming esults. The paper is structured into the following (independent) sections: • Acoustic Imaging with very short measurement intervals. • Breaking up the array into sub-arrays. • Modeling the phase distortions. • Short-interval images under different conditions. • Statistics: number of short-interval images of ‘‘good’’ quality. • Experimental results: acoustic speaker & aircraft model in a 8 × 6 m2 open jet wind tunnel.


Introduction
Measuring the strength and location of sound sources is crucial to the development of silent technologies. Acoustic imaging techniques allow for visualization and quantification of sound sources in wind tunnel experiments. An array of microphones samples the sound waves emitted from the sound sources, similar to a lens collecting light waves. The microphone signals are focused digitally to obtain the acoustic image. This focusing method is knows as the beamforming method [1,2]. A requisite for beamforming is sufficient coherence between microphone signals. Considering one sound source, the signals measured by different microphones should be identical up to a time shift and amplitude scaling. However, when the acoustic signal propagates through a shear layer the turbulence continuously alters the acoustic propagation. This leads to random distortion of the microphones signals and consequently the coherence between two microphone signals is reduced.
Coherence loss between microphone signals is undesired because it reduces the peak value of the Sound Pressure Level (SPL) in the acoustic image [3]. This hampers the accurate evaluation of the acoustic source strength. Furthermore, the resolution of the acoustic image is reduced, making it difficult to separate individual sound sources. Oerlemans and Sijtsma [4] investigated acoustic measurements performed on a scaled aircraft model placed in the open jet of a large industrial wind tunnel. The acoustic waves traveled through a thick shear layer to the microphone array placed outside the jet. For frequencies ranging from 8 kHz to 12 kHz a difference of 10 -15 dB was observed between the peak SPL in the acoustic image and the measurement of a single microphone at 50 m⋅s −1 and an estimated shear layer thickness of ∼1 m. The results also indicate that spatially integrated source SPL compared favorably with the SPL measured by a single microphone, the difference being less than 2 dB. This suggests that the largest error is found in the peak level and not in the integrated level. There are several ways to increase the signal coherence and consequently the acoustic image quality. The signal coherence is improved by removing the shear layer, e.g., using a closed test section, or a hybrid test section [5]. However, this is not always a possible or desirable solution. Improved post-processing techniques, able to correct for coherence loss, were therefore investigated. Several techniques have been developed for use in aeroacoustic measurements, nautical applications, or astronomy. In the field of astronomy images of stars or planets are blurred by atmospheric turbulence influencing the propagation of light, analogue to the blurring of acoustic images by shear layer turbulence.
Research regarding the loss of coherence between two microphone signals is closely related to the distortion of the signals measured by the individual microphones, commonly used to measure the acoustic spectrum. The turbulence induced fluctuations in propagation time cause a tonal peak in the spectrum to spread over a wider frequency range. This phenomenon is known as spectral broadening, or more colloquially as ''hay-stacking''. A fundamental theory regarding the scattering of sound due to turbulence was derived by Howe [6]. This theory included the absorption, scattering and spectral broadening of the acoustic spectrum. The results regarding spectral broadening and scattering agreed with predictions based on geometrical acoustics. Aside from the spectrum measured by a single microphone additional information regarding the phase of the sound waves is obtained from the crosscorrelation between two microphone signals. Yang [7] used the theory of geometrical acoustics combined with a path integral to compute the coherence between two rays (signals) moving through a medium disturbed by linear waves. In a more recent study McAlpine and Tester [8] derived an analytical model to describe the two-point correlation for sound propagation through homogeneous axisymmetric turbulence. This model showed better agreement with measurements compared to models based on homogeneous isotropic turbulence.
The signal measured by the microphone equals the signal emitted by the sound source shifted by the acoustic delay time. This acoustic time delay can be averaged to obtain a time-invariant delay caused by the finite propagation speed of the sound and advection by the mean flow. The time-variant delay, i.e., the deviation from the average, is due to the turbulent shear layer. When the shear layer is thin compared to the wavelength the time-invariant delay time is predicted by the method of Amiet [9,10]. In this case the influence of the turbulence can be neglected. However, when the shear layer is not thin compared to the wavelength and the shear layer is strongly turbulent the amplitude of the time-variant delay increases. Freund and Fleischman [11] compared ray-traced acoustics with and without turbulence. Their most notable finding was that the evolving turbulence allowed the sound rays to more easily escape the potential core when traveling upstream. Koop et al. [12] investigated the time fluctuations as experienced by the microphone array. They performed numerical as well as experimental studies to predict and measure the phase variations due to turbulence interacting with the sound waves. The use of a known source allowed the time fluctuations to be measured. Coherence loss in acoustic measurements was also observed by Dougherty [13] in the form of an apparent increase of the sound power level when switching from a larger diameter array optimized for low frequencies to a small diameter array designed for higher frequencies.
A model for wave-coherence based on the work of Tatarskii [14] was used to predict the coherence loss. The model relies on the use of a statistical description of the turbulent fluctuations to assess the decorrelation of the sound waves. A structure function was used to determine the correlation between the turbulence and the fluctuations in acoustic propagation time. Blacodon [15] derived a correction method using a known reference source to derive a deconvolution operation to remove the smearing of the sources in the acoustic map. The author compared the integrated spectrum of an open test section measurement, computed using the deconvolution operation, with the spectrum of a single microphone measurement in a closed test section. The results compared favorably up to 10 kHz. Tatarskii's theory was also used by Ernst et al. [16] to predict the coherence loss due to shear layers in open jet wind tunnels. With the model of Tatarskii Ernst et al. predicted the magnitude of the coherence loss. They conclude that an impractical measurement time is necessary for statistical convergence of the coherence estimator in the case of severe coherence loss. Pires et al. [17] used Tatarskii's model to study the coherence loss caused by the boundary layer of a closed test section wind tunnel. Both Ernst et al. and Pires et al. performed flow measurements and compared the experimentally obtained acoustic coherence loss to predictions based on experimentally measured structure functions. The model captured well the coherence loss for larger distances between microphone pairs. This corresponds to the limit where the turbulence encountered over the acoustic paths is no longer correlated. Biesheuvel et al. [18] used Tatarskii's model to assess the effect of coherence loss on beamforming results due to thick boundary layers (∼0.2 m) as encountered in industrial scale wind tunnels. For frequencies between 20 kHz and 50 kHz a decrease in beamforming SPL was predicted ranging from 2 dB to 5 dB. In experiments with models the frequencies of interest during the measurement scale inversely proportional to the model scale. Considering a 1:10 model scale, the frequencies on the full scale range from 2 kHz to 5 kHz, which is well within the audible range.
A simple method to increase overall microphone coherence is to discard the incoherent microphone signals. In general, microphone weighting methods are known as shading techniques [19]. In practice the microphones near the outer rim are discarded in the processing, effectively scaling down the microphone array dimension. Microphones that are spatially close are more likely to measure similar wavefront distortions and retain their signal coherence. Severe loss of signal coherence may even limit the statistical convergence of the test results [16]. Shading based on the coherence of microphone pairs was used by Amaral et al. [20], who defined the weighting factor of a microphone as the mean coherence between the microphone and all other microphones of the array. This weighting improved the results of a deconvolution technique [21]. Bahr and Lockard [22] developed shading weights based on the distance between the microphones and the frequency. They note a visual improvement for beamforming maps, judged by the absence of blurred sources. In general shading can recover some of the lost resolution by effectively removing microphone pairs that, due to the coherence loss, have a cross-correlation that is not converged. The downside of applying shading is that the effective diameter of the microphone array becomes smaller, which also decreases the acoustic image resolution. Koop et al. [12] J. Biesheuvel et al. hint at the possibility to use a known source as a ''guide star''. This was done by Sijtsma [23] who used the presence of a few known sources to estimate the wavefront disturbance and correct the phase of the microphone signals. A similar algorithm is used to correct images of stars using known stars or lasers, hence the technique is named ''guide star''. The propagation time from the rest of the wind tunnel model could be estimated through interpolation between speaker locations. This method improved the beamforming results near the reference speakers. The technique required a large computational effort. A time domain beamforming method is presented by Cousson et al. [24]. Although, their intention was to present a method to identify moving acoustic sources the method implicitly incorporates a non-stationary delay time. When the acoustic source is allowed to move over time, this also implies that the acoustic propagation is time-dependent. A localization method for acoustic sources in random moving media by means of beamforming was presented by Gay et al. [25]. The authors used Kirchhoff migration (back propagation of the measured signals) and coherent interferometry (back propagation of the signal cross-correlation [26]) to localize acoustic sources in a jet. The results were favorable, given that the turbulence induced perturbations were not too large.
In the field of nautical research acoustic measurements are degraded by ocean turbulence. Unknown fluctuations in propagation time degrade results obtained for passive ranging. Ge and Kirsteins [27,28] developed a method to classify acoustic snapshots as either coherent or incoherent. Before processing, the data was ranked and the highly coherent snapshots were accumulated, taking advantage of the stochastic nature of the distortions. Ge and Kirsteins then demonstrated a reduction of smearing effects in range-bearing plots.
Looking further, into the field of astronomy, similar issues are encountered. Loss of resolution occurs when temperature and density variations due to atmospheric turbulence alter the local refraction coefficient of the air, resulting in a distorted image of the stars. One algorithm used to correct for this is called the ''Luck exposure'' [29] or ''Lucky Imaging'' [30] algorithm. For a particular section of the sky a large number of images are made using a very short exposure time. Only the images with good resolution are further processed. In this manner, the algorithm exploits the stochastic nature of the turbulence, which implies that some of the images acquired are less disturbed by the turbulence than others. When the observation time is long enough a subset of ''good quality'' images with high resolution is selected for further processing. Dantowitz [31] used the method to create highresolution images of the planet Mercury. Law [32] also used to method to increase the resolution of an image of the stars, the increased resolution made more faint objects visible. Baldwin [33] showed an improvement in image quality by retaining only the ''best'' 1% exposures. In a later study Baldwin [34] performed numerical simulations using the random phase screen method. In this study a Kolmogorov spectrum was used to simulate the turbulence. The results compared favorably with results obtained from an observatory.
The most notable efforts to explicitly correct acoustic beamforming maps for turbulence induced fluctuations have been made by Sijtsma, Blacodon and Gay et al. [15,23,25]. In addition to these efforts we propose a novel method to correct acoustic images based on the ''Lucky Imaging'' technique. This technique offers various advantages: (1) it does not require the use of ''guiding'' sources on the model or reference acoustic sources before or after the test. (2) The method is relatively simple and easy to implement. (3) The ''Lucky imaging'' methodology relies on filtering and not on the ensemble average of the microphone cross-correlations. This difference is mainly relevant when performing acoustic tests in large industrial wind tunnels which have a thick shear layer or at high frequencies. The large coherence loss may then lead to excessive measurement times.
In this paper the Lucky Imaging methodology is developed to be used in acoustic wind tunnel tests [35]. Fig. 1 shows a schematic overview of the (Acoustic) Lucky Imaging method. The method consists of three stages. In the first stage acoustic images are obtained from measurement intervals short enough to ''freeze'' the turbulence induced distortions. In the second stage those acoustic images that are severely distorted are removed from the processing. The third, and final, stage combines all the ''good'' images into a single acoustic image. To support the correction methodology a model for coherence loss is derived. This model is based on the discretization of the wavefront distortions. The model is used to show why and when Acoustic Lucky Imaging improves beamforming results. The paper is structured into the following (independent) sections: • Acoustic Imaging with very short measurement intervals.
• Breaking up the array into sub-arrays.
• Modeling the phase distortions.

Acoustic imaging with distorted wavefronts
Acoustic beamforming is used to show the strength and location of sound sources, and is a method that numerically generates acoustic images from the measured microphone signals. As the microphones have different locations on the array the sound waves reach microphones at different times. The beamformer uses the time differences to localize the sound source. To estimate the time difference it is convenient to transform the microphone signal into the Fourier domain and work with the phase difference instead. This allows focusing on specific frequencies of the sound source. Mathematically the beamformed map is generated by evaluating the integrals of Eq. (1).
Each point ′ , ′ in the source distribution 0 is mapped to a point , in the image by a transfer-function . The capital , 0 , and are all complex numbers and therefore have both a phase and a magnitude. The integral of all the measured points on the array yields the response of the total array. The time integral indicates the finite measurement time, where 0 is the start of the measurement and the duration. A typical acoustic measurement time ranges from a few seconds to a few minutes. The source here is assumed stationary. However, the transfer-function is assumed non-stationary as this includes the effect of turbulence on the propagation from the source, through the turbulent shear layer, to the microphones. Due to the turbulence of the shear layer the wavefronts are disturbed. The transfer-function thus inherently depends on the turbulence of the shear layer and is time dependent. For a very short measurement time , in the order of milliseconds, the distortion of the wavefront is approximately constant. Eq. (1) is then approximated as: The acoustic image is still dependent on 0 in the transfer-function, because the turbulence velocities are different for each measurement. Consequently, the amount of image distortion also varies with 0 . A propagating wavefront reaching the array is focused as shown schematically in Fig. 2(a). For undisturbed waves the beamforming algorithm will show a virtual sound source at the location of the real sound source. The turbulence induced velocity fluctuations disturb the wavefront. When the distortion is linear the wavefront remains planar and the virtual source location only appears shifted, see Fig. 2(b). The disturbances can also be non-linear. However, locally the wavefront is still approximately planar and the virtual source is split into multiple virtual sources as shown in Fig. 2(c). The acoustic image of a point source is fully described by the transfer-function . A schematic representation of the time dependent transfer-function is shown in Fig. 3(a). Because the transfer-function is changing over time the average transfer-function, see Fig. 3(b), is of lower resolution than the theoretical transfer-function. Lucky Imaging methods aim to use filtering and shifting of the instantaneous transfer-functions to produce a sharp average transfer-function.
The specific implementation of the (Acoustic) Lucky Imaging method in this paper is discussed in order to create a context for the rest of the paper. The principle is to compute acoustic images from short data blocks. Typically this results in about 5,000 to 50,000 acoustic images. The images are then filtered to remove images that are severely distorted, this stage typically rejects 50% to 99% of the images, depending on the severity of the distortions. The images that have passed are then further processed to obtain a final acoustic image. The algorithm implemented for this work is summarized as: • Splitting of the acoustic data signals into blocks. The time duration of a single block should be sufficiently short such that Eq. (2) approximates Eq. (1). • Perform a Fast Fourier Transform (FFT) on a single data block for all microphone channels.   [36].
• Select which frequencies to include in the processing.
• For the selected frequencies compute an acoustic image for all data blocks. These images are referred to as the ''short-interval images''. • Using the short-interval images compute the average image for each frequency. This is used as the reference image.
• Compute the normalized (spatial) cross-correlation between the short-interval images and the reference image.
• Check that the normalized cross-correlation is larger than a chosen threshold for a minimum of frequencies, e.g., 8 out of 10 selected frequencies should have a normalized cross-correlation larger than 0.8. • Remove the data block from the processing if the requirement is not met.
• Using the normalized cross-correlation, compute the average displacement for a single block using all frequencies. Spatially shift the images corresponding to the block using the computed average displacement. • Compute the average of all the shifted images on a per frequency basis. This is the final acoustic image.
A supporting theory is developed in order to obtain a relation between the improvement of the processed acoustic image and the distortion of the wavefront. This theory is not necessary to implement the steps above which work with conventional microphone arrays and conventional beamforming methods. Another key issue is the number of high quality images that are obtainable in a reasonable measurement time. The next sections elaborate on these issues, after which the experimentally obtained results of the above algorithm are presented.

Discretization of the wavefront and the response of a continuous array
The wavefront is subdivided in sections that are sized according to the linearly disturbed wavefront segments. In this section the response of a continuous acoustic beamformer due to the subdivided wavefront is computed. Fig. 4 shows a schematic representation of a typical acoustic setup. The acoustic array, with diameter , is at a distance from the acoustic source, located on the -axis. The two-dimensional array or ''lens'' is parameterized by the coordinates and and the focus-plane, i.e., the beamformed map, is parameterized by the coordinates and . The source acts as a perfect monopole. The beamformer only corrects for time delays, thus pressure corrections are omitted. Furthermore, it is assumed that all phases have been corrected for mean flow advection, and therefore the wave theory valid in quiescent air is used. With these assumptions, the incoming wave at the acoustic-lens and the focus-plane is written as: where is the pressure at the array, ′ is the array pressure propagated to the focal plane, is time, the wave number, the angular frequency, and the is a distance function that computes the distance between two points. The quantity of interest is the difference ( ; ; ; ) − ( ; ) which is approximated using a Taylor expansion, resulting in: This approximation is better known as the Fresnel approximation and the derivation is found in Appendix A. The approximation is valid when: J. Biesheuvel et al. To reduce the symbols involved in the mathematical expressions the scaled wavenumber̂= ∕ is used in the rest of the paper. Integrating the contributions ′ over the array surface gives the total beamformed response, i.e., where the right hand side integral can be recognized as the Fraunhofer diffraction equation. The array is now discretized by subdividing the array into sub-arrays. The shape of the sub-array is a hexagon as shown in Fig. 5 (similar to the James Webb telescope). The hexagon is parameterized by two length scales and . The length scales model the size of the distortions, which are caused by the turbulent shear layer. Using two length scales allows the length scale in streamwise direction to be different from the length scale in spanwise direction. For mathematical convenience a slope is defined, which is entirely defined by the length scales. For an isogonal hexagon is equal to 1 2 √ 3 , from which follows that for arbitrary and the expression for is equal to 1 2 ∕ since scaling is a linear operation. Physically this assumes that both length scales of the distortions are equal. The center of the sub-array is at ( ; ). The contribution of a single hexagon to the total beamformed response is found by integrating Eq. (7) over the hexagon. The integration is performed in Appendix B, the result is the following expression: Thus the response of a sub-array is equal to some constant , that incorporates the exact shape, with a phase shift, that depends on the location in the array. The total beamformed response is equal to the sum of all sub-arrays .
The constant has been moved out of the sum and is the number of sub-arrays. The first term produces a phase shift depending on the location in the focal plane. Since a phase shift does not change the magnitude it is not of importance for a beamformer. The second term is governed by the shape of the sub-array. The last ''resolution'' term is shown schematically in Fig. 6(a). Near the origin this term is effectively a sum of cosines, and the coordinates and produce different wavelengths. The only shared maximum of these waves is at zero and consequently, this produces a sharp peak at the origin. The last term effectively states that the sum of the sub-arrays achieves a higher resolution than a single element. The subdivision of the array is done numerically using an infinite tiling while discarding sub-arrays outside the array. The remaining sub-arrays all have unique center coordinates and .

Phase distortions
Assuming the sub-arrays are of a size over which the wavefront distortion can be assumed constant, the beamform equation, Eq. (10), is adapted by adding phase shifts to the sum of sub-arrays: where is the delay time distortion. Without the delay time distortion the sum of exponentials produces a sharp peak as illustrated in Fig. 6(a). Fig. 6(b) shows the same sum with the delay time distortions that cause misalignment of the waves. The different wavelengths are due to the ( + ) term in the exponential. This reduces the ability of the beamformer to locate the sound sources correctly. In most beamforming applications the final quantity of interest is the source power, i.e., the square of the amplitude of the sound source. This quantity is computed from the beamformed response by multiplying with the complex conjugate: Which due to the symmetry can be written in terms of sines and cosines, and split into a response due to the auto-powers and cross-powers: This expression computes the source power using the beamforming method. The assumption is that the time fluctuations are constant, similar to Eq. (2). In general, beamforming results are averaged quantities rather than instantaneous quantities, e.g., the mean source power. A generic beamforming equation for the instantaneous source power can be written as: where is a complex valued steering vector which acts as a transfer-or Green's function and is the cross-spectral matrix, containing the Fourier transform of the microphone cross-correlations. The cross-spectral matrix is computed from the measured microphone signals as = * . The average of the beamformed source power follows as: J. Biesheuvel et al. where the latter equality holds if the steering vectors are not depending on time. In other words one obtains the same result either by computing the average of the instantaneous images or beamform using the average cross-spectral matrix. The latter approach is considered as the standard. To show the effect of the time fluctuations on the average source power, the expected value of Eq. (14) is computed: where is a probability density function (PDF) for the time fluctuations . For simplicity a Gaussian PDF is used. Thus is equal to: where is the variance of the differences calculated directly from the variance of the fluctuations as = √ 2 . As found in Gradshteyn [37] (section 17.23, entry 13), the expected value of the sine and cosine parts are conveniently found as: With = ( √ 2 ) −1 and = this results in: The expected value of the beamformed power becomes: ] ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ cross-powers shown in Fig. 7(d). The contribution of the auto-powers is shown in Fig. 7(c). The auto-power term behaves as if an array with the size of the sub-array had been used, and thus features poor resolution due to the Rayleigh criterion. Fig. 7(b) shows the response only due to the cross-powers, this result illustrates why diagonal removal improves the image quality. The cross-powers, containing the phase information, cause the focussing effect. The cross-powers are however, subject to the time fluctuations, i.e., the exponential term, and are suppressed if variances are large. The optimal response, with = 0, is shown in Fig. 7(a). Eq. (21) also provides an explanation why integrated beamforming levels often match the spectra obtained from single microphone measurements. In case of severe coherence loss, and if the object is directly in front of the array, the response of the beamformer corresponds to the averaged auto-powers. However, no directional information is present in the integrated levels and thus individual sources cannot be distinguished.

Acoustic Lucky Imaging & short-interval acoustic images
Acoustic Lucky Imaging operates on short-interval acoustic images. Therefore, the nature of the short-interval images very much influences the results. The three operations that Acoustic Lucky Imaging applies are filtering, spatially shifting, and averaging.  During filtering the severely distorted images are removed, e.g., in this section the filtering criterion is based on the maximum response of the short-interval image which must be greater than a threshold. The images that have passed the criterion are then spatially shifted according to an alignment criteria, e.g., align the maximum response. Finally, the images are averaged to obtain a final result. In this section the influence of the magnitude of the delay time distortions and the size of the wavefront distortions is discussed. Depending on the wind tunnel test environment four different regimes of distortion are identified. The regime determines how much the filtering and shifting will improve the results. Fig. 8 schematically shows the relation between the test environment and the wavefront distortion. If the turbulence structures are small the disturbed wavefront is approximated by more piecewise constant segments than if the turbulence structures are large. The amplitude of the time delay distortions is assumed to increase with the wind tunnel velocity, hence higher wind tunnel velocities are expected to increase the wavefront distortion. The response of a beamformer in case D, a low amplitude linear distortion, is shown in Fig. 9(a). The beamforming images obtained from very short-intervals of measurement data are schematically shown in Fig. 9(b). The response of the average beamformer can be improved by discarding instantaneous images with a maximum response below a threshold. The remaining images feature a high maximum response, but due to the linear disturbance the maximum response is at the wrong location. This is corrected by shifting the image such that the maximum response it at the origin again. To obtain statistical convergence the selected and shifted images are averaged to obtain a final beamformed image, see Fig. 9(a).
If the wind tunnel velocity is increased the linear wavefront distortions increase in amplitude, this is case B in Fig. 8. The principle of image selection and subsequent shifting of the maximum response can still increase the beamformed image quality. However, due to the increased amplitude of the wavefront distortions large shifts are to be expected as shown in Fig. 10(b). Furthermore, the response of the instantaneous images decreases when the source is substantially shifted. Thus even with selection and shifting a loss of peak power is inevitable. Still a significant improvement with respect to conventional beamforming is expected since the large shifts in source location rapidly decrease the average peak as shown in Fig. 10(a).
In case C the low amplitude wavefront distortion is non-linear. This results in instantaneous images with a maximum response near the origin. The wavefront distortion will decrease the peak response. In this case shifting the images is not expected to increase J. Biesheuvel et al. As well can the turbulence be strong or weak, which is proportional to the wind tunnel speed. the final image quality. However, as shown in Fig. 11(a), discarding images that are too distorted still improves final results. Since the majority of samples is rejected this does require longer measurements times. Finally, Fig. 12(a) shows case A where the wavefront is severely distorted in a non-linear manner. The instantaneous images, as shown in Fig. 12(b), are so severely distorted that all images are discarded. In this regime the method is not expected to function. This section showed how the filtering and shifting operations of Acoustic Lucky Imaging improves beamforming results. The improvement is obtained by first removing the data that is severely distorted. Subsequently, the selected data is aligned and averaged to obtain a statistical significant result. In the next section a methodology is derived to estimate the number of required instantaneous images. For explanatory reasons we have adapted the simple quality criterion based on the maximum response. In practice this is not a good approach. Improved alternatives are discussed in Appendix C.

Success rate of Acoustic Lucky Imaging
Acoustic Lucky Imaging discards measurement data that is severely distorted. The useful data should be obtainable within a reasonable measurement time for the Acoustic Lucky Imaging method to be applicable. Ideally an a-priori estimate provides the percentage of expected high-quality images so that the required measurement time can be calculated. In this section a simple estimation method is developed.
In astronomy, Fried [38] derived an expression for the probability of obtaining usable data, however, this does not easily translate to Acoustic Lucky Imaging. An alternative method by Biesheuvel et al. [18] solved this issue by fitting an approximate PDF using J. Biesheuvel et al.
It is assumed that the maximum response is near the origin and thus and are small. Furthermore, the normalized wave number should also be small. With these assumptions the expression for the source power reduces to: The maximum response with no distortions is equal to one, therefore, the value of must equal 1∕( −1) 2 . We now require that each cosine term contributes at least to the total sum. Requiring this for terms containing the 1th microphone, yields the inequality: J. Biesheuvel et al. where is the index of any other microphone. Rewriting to obtain an inequality for leads to: where is used as a shorthand for the cosine term. The phase differences between other pairs of microphones follows from their definition, i.e., If the phase differences between the 1th microphone and all other microphones are smaller than ∕2, the phase differences between all other possible pairs of microphones must be smaller than . The probability that the phase differences between the 1th microphone and all other microphones is smaller than ∕2 is given by: The PDF of follows from Eq. (18). The cumulative PDF is found by integrating Eq. (18), and equals: where is the variance of the phase differences. This results in the probability that all shifts are within being equal to: For a fixed threshold the probability for a given number of accepted images becomes smaller with increasing frequency. Furthermore, non-linear distortion of the wavefront must be modeled by a higher number of sub-arrays , leading to a lower probability of finding good images. Similarly, the large time fluctuations increase the variance and decrease the probability. In practice the values for and are not always known. However, a speaker representing a monopole source could be used to compute the cumulative PDF and subsequently, the and are computed with a fitting procedure.

Numerical example
In this section a numerical example is provided to showcase the Acoustic Lucky Imaging method and how to use Eq. (30) to make estimates. A set of 15,000 short-interval beamforming images is generated numerically, see Fig. 13 For the Acoustic Lucky Imaging methodology the threshold was set to = 0.71. Thus the images with a maximum value above 0.71 are selected. The instantaneous images, one of which is shown in Fig. 13(a), feature an off center maximum response. This significantly decreases the average response plotted in Fig. 13(b). Discarding instantaneous acoustic images with a maximum below J. Biesheuvel et al. the threshold results in Fig. 13(c). A significant improvement can already be obtained by just removing severely distorted images. The result can however, be improved more by also aligning the maximum responses of the selected instantaneous acoustic images, the result is a response that is almost restored as shown in Fig. 13(d). In this example of 15,000 instantaneous acoustic images 312 passed the threshold criteria, which is 2.1%. The estimate based on Eq. (30) is 2.2% as shown in Fig. 14(a). The distribution of the maximum responses of all 15,000 acoustic images is shown in Fig. 14(b), this allows the computation of an approximate cumulative PDF. From the cumulative PDF and the threshold the number of good images can also be estimated. This can be compared to the analytical estimate based upon Eq. (30), denoted by the ×, both show indeed a similar probability.

Application to a large open jet wind tunnel experiment
To demonstrate the applicability of the method to experimental data, measurements were performed in the acoustic wind tunnel (LLF) of the German-Dutch wind tunnels [39] with an open jet test section. A calibrated speaker source was placed inside the potential core of the wind tunnel using a sting, see Appendix F. On the floor an acoustic array, with 140 microphones, was mounted on a traverse system. The traverse allowed to change the relative angle of the array with respect to the source and change the acoustic path length through the turbulence. Turbulence induced time delay variations were measured using the procedure discussed in Appendix D. The time delay variations allow to cross-validate the results of the Acoustic Lucky Imaging method. The wind tunnel velocity was varied from 34 m⋅s −1 to 68 m⋅s −1 . The wind tunnel nozzle dimensions are 8 m × 6 m. The sample time was set to 15.26 μs (2 16 = 65, 536 Hz). A high pass filter of 500 Hz was applied. The measurement time was equal to 45 s. Two signals were used as input to the speaker. The first was a saw-tooth signal with a fundamental of 2 kHz. The second signal was a white noise signal, band filtered between 5 kHz and 20 kHz.
The block size during post-processing was 256 samples which makes the short-interval equal to 3.91 ms. With this data 10,000 short-interval images were generated. Combining this short-interval with the average flow velocity the turbulent structures will J. Biesheuvel et al.   have moved 66 mm to 134 mm within the interval. The selection criterion was based upon a normalized image cross-correlation. The shifts were also obtained using the same normalized image cross-correlation. The input parameters are listed in Tables 1-2.
The value of indicates a threshold value of the normalized cross-correlation between the short-interval images and the reference image. If a short-interval image correlates less than this value it is rejected. The selection can be performed independently for multiple frequencies. Subsequently, a second requirement is that multiple frequencies have a normalized cross-correlation value greater than . The second number in the row labeled A indicates how many frequencies are analyzed (6 or 10). However, for each frequency the normalized image cross-correlation may be higher or lower than the threshold. The first number indicates the minimum number of images that must have a higher normalized cross-correlation, e.g., 8 out of 10. For more details on the subject of selection criteria see Appendix C.

Time fluctuations
The method proposed by Koop and Sijtsma, see also Appendix D, is used to determine the shear layer induced time fluctuations. The method requires that the carrier frequency is clearly measurable and that the spectral broadening effects are not overlapping. Fig. 15(a) shows the Welch power estimation of the reference signal with no flow and the power estimation with the wind tunnel at the maximum velocity to show that these requirements have been met. Furthermore, a measurement without any signal is plotted to show the background noise of the wind tunnel. The measured time fluctuations are summarized in Fig. 15(b). A wrapped normal distribution has been fitted to compute the variance of the fluctuations. The dependence of the variance of the fluctuations on the wind tunnel velocity is shown in Fig. 16(a). To obtain some redundancy the time fluctuations obtained from two carrier frequencies have been analyzed, both yielding similar results, indicating the proper functioning of the method. Fig. 16(b) shows an estimate of the sub-array size. These estimates have been obtained by fitting the PDF of measured normalized cross-correlations to a Monte-Carlo simulation based on Eq. (13). This is however, a rough estimate due to the discrete nature of the theory which only allows for an integer number of sub-arrays regions. The sub-array sizes have been averaged over different frequencies. Normalizing the sub-array size with the array diameter yields a value of 0.2 m, with an average ≈ 0.04 ms, the measurement appears to be similar to case C from Fig. 8.

Source power levels
This section presents the results obtained with Acoustic Lucky Imaging. Fig. 17(a) shows the maximum response for Acoustic Lucky Imaging and conventional beamforming as function of frequency. The benchmark is the acoustic spectrum as measured by an individual microphone. The auto-powers are used as a comparison since they are not effected by the coherence loss. The wind tunnel velocity was 34 m⋅s −1 . The Acoustic Lucky Imaging is able to correct 2-−3 dB in the frequency range shown. Above 10 kHz the coherence loss becomes dominant and the improvement of the acoustic maps decreases. Fig. 17(b) shows the results as function of the wind tunnel velocity at 9.1 kHz. The increase of turbulence with increasing wind tunnel velocity causes the beamforming process to show lower maxima. The Acoustic Lucky Imaging methodology is able to partially recover this. Examples of the final results are shown in Appendix E, see Fig. E.25, as well as some short-interval maps in Fig. E.26. The improvement of resolution is measured using the main lobe width. The main lobe width, based on the 3 dB down criterion, decreased with a factor 1.28-1.50. The computation time for processing the noise source was 6 minutes for 10 frequencies and a grid of 128 × 128 points, using a computer equipped with an Intel ® Core™i7-6700HQ processor and a NVIDIA Quadro M1000M graphical processor.

Parameter sensitivity
The parameters in Tables 1-2 are chosen such that a minimum number of images remained after the selection procedure. In this section we would like to expand on these choices. Fig. 18(a) shows the PDF for the normalized cross-correlation values for different J. Biesheuvel et al.  frequencies. Most of these distributions show a sharp peak exactly at the point where the threshold values would be, i.e., in the range 0.8 to 1.0. This means that increasing the threshold by a small amount can cause the amount of rejected images to decrease relatively fast. The image quality can be computed independently for a number of frequencies. For the noise signal the computation was performed for 10 frequencies. If the wavefront distortion is indeed small the normalized image cross-correlation should be high for most frequencies. Fig. 18(b) shows the number of passed images as function of the minimum number of frequencies with a normalized cross-correlation greater than the threshold, e.g., 6000 images have a normalized cross-correlation higher than 0.83 for at least 6 frequencies. This criterion can be used in reverse, e.g., further processing requires at least 500 images, then an image should pass for at least 9 frequencies.

Acoustic Lucky Imaging on aircraft model: showcase
The Acoustic Lucky Imaging method is used on acoustic data obtained in a test with a scaled aircraft model, to showcase the method in an industrial setting. The wind tunnel velocity equaled 62 m⋅s −1 . The frequency of interest is 2.56 kHz. The result obtained with conventional beamforming is shown in Fig. 19(a). The image generated by the Acoustic Lucky Imaging method is shown in Fig. 19(b). The peak SPL has been increased by 6 dB. Also we observe that the greatest difference is near the dominant sound source at (x = −0.8 m; = 1.1 m). The difference in improvement is probably due to the non-linearity of the wavefront distortion. The smaller the region-of-interest the more linear the distortion. The spatial shifts for the strong source are probably not the same for J. Biesheuvel et al. Fig. 19. Acoustic measurements of a scaled aircraft model. the weaker source at the end of the wing, which therefore, is not increased in the same proportion. To process the acoustic map of the whole aircraft model Acoustic Lucky Imaging must be employed on multiple sub-regions. This would require a method to stitch processed maps together. This issue is similarly encountered in the ''guide star'' method, where the delay time distortions corrections are only valid near the reference sources.

Conclusion
Acoustic Lucky Imaging can be used to improve the resolution of acoustic images obtained using beamforming when the acoustic waves have passed through turbulence. Beamforming images based on very short measurement intervals proved useful for further processing. Subsequently, the images that were severely distorted by the shear layer turbulence are discarded from further processing. Linear wavefront distortions are corrected for in a post-processing step by spatially shifting the acoustic image. The methodology is supported by a mathematical model. This model shows the requirements for the methodology to be successfully applied. The model also allows to predict the number of images rejected.
An experiment with a loudspeaker in a large industrial wind tunnel confirmed the applicability of the proposed methodology. The wind tunnel featured an open jet test section and has nozzle dimensions of 8 m × 6 m. The methodology was able to increase the source map resolution. Furthermore, an increase in source power up to 3 dB was observed with respect to conventional beamforming, at a frequency of 8 kHz, a wind tunnel speed of 34 m⋅s −1 , and a shear layer thickness of ∼1 m. Compared to the microphone autospectrum a difference remained present, this is in accordance with the derived theory, where Acoustic Lucky Imaging can only partially recover the undisturbed source map. Another wind tunnel test using an aircraft model showed an improvement of 6 dB.
Further research will focus on extending the basic Acoustic Lucky Imaging framework. The three stages (image generation, image grading, and image combination) could all be replaced by more elaborate schemes. A further improvement would allow for a finer frequency resolution.

Funding sources
This research was developed in the context of the Silent Approach project receiving partial funding from the Top Technology Twente program in the framework of the TKI-HTSM roadmap Aeronautics.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The authors do not have permission to share data.

Appendix A. Fresnel approximation
The mathematical model presented in Eq. (3)-(4), in Section 3, is simplified by approximating the distance function . First, the distance is written in a simple parametric form: with dummy parameters and . Subsequently,̂is approximated with a multivariate Taylor expansion up to second order, i.e., This yields as the approximation: In optical terms we are performing a Fresnel approximation. The approximation is valid if the ignored third order terms make a neglectable contribution to the phase in the exponential. This is true if the value is much smaller the period 2 , i.e.,

Appendix B. Beamforming integrals
The contribution of a single hexagon to the total beamformed response is found by integrating Eq. (7) over the hexagon, which is conveniently expressed in terms of as: . (B.1) The first integrals are integrated as follows. After which a second arrangement follows, that takes into account the complex conjugates. Noting that the cosine function is even, i.e., cos( ) = cos(− ), further simplification is possible.

Appendix C. Acoustic Lucky Imaging implementation
A general Acoustic Lucky Imaging algorithm implementation consists of three stages, see Fig. 1, repeated below in Fig. C.20. In this appendix we discuss the concrete implementation details for each stage.

C.1. Stage I: short-interval map generation
Acoustic Lucky Imaging works by sampling acoustic arrays for a very short time and subsequently generating beamforming images. The short interval should be chosen as to temporally freeze the turbulence, or more specifically the time variations as measured on the array. Beamforming however, includes a Discrete Fourier Transform (DFT) of the acoustic data. This requires a minimum number of waves to be captured within the measurement time in order for the DFT to have statistical meaning. Furthermore, shorter measurement times translate to a smaller data block size, which in turn yields large frequency bins in the DFT domain. This can cause background noise to suppress the sound source. Thus the block size has to be as short as possible, while still featuring minimal signal-to-noise ratio.

C.1.1. Stochastic noise sources-block overlap
The Acoustic Lucky Imaging methodology requires the use of very short block sizes. For periodic sound sources this does not require extra processing. However, a situation shown in Fig. C.21 may occur for noise sources that are measured on large arrays or under large aperture angles. The data blocks traveling along the acoustic rays arrive at the microphones. However, the difference in arrival time is larger than the block size. This means that during post-processing two blocks of uncorrelated noise are used, leading to spurious results. This is easily prevented by reading the data blocks offset. The offset is computed using the center source location and converted for each data channel into an integer offset using the sample time . After computing the DFT of the signal the data is shifted back in the DFT domain by the complex rotation . Both expressions are given as: where the ⌊… ⌉ denotes the rounding operation.

C.1.3. Speeding up the process for a single data block
If only a single block of data is used it is not necessary to compute the cross-spectral matrix in order to compute the beamforming image. The source power is conventionally computed as: where is the source power, the cross-spectral matrix, and are steering vectors. The cross-spectral matrix follows from the complex pressure signals as = * . However, if is based on a single block of data is more efficiently computed by first beamforming with the acoustic pressure signal and then squaring the result to obtain the source power, i.e., where is a normalization factor. Similarly, beamforming without the auto-powers can be done more efficiently in a two step process, i.e., An added benefit is that these formulations allow easier parallelization of the computational code.

C.2. Stage II: evaluation & discarding
A critical point of Acoustic Lucky Imaging is the definition of the image quality or ''goodness''. A simple criterion is the maximum value found in the beamformed map. The rationale behind this approach is that severe coherence loss in a short-interval image leads to a spread of the available source power across the image, thereby, decreasing the maximum value. Thus short-intervals with little coherence loss should have a high maximum, and therefore be selected. However, if a source is distributed this may not be a valid criterion. Furthermore, acoustic sources may vary in strength and thus the cross-correlation value varies accordingly. This can especially be the case for a non-stationary source such as (aeroacoustic) noise.

C.2.1. Distributed sources & frequency ranges
To allow for the processing of general shaped source distributions the short-interval images are cross-correlated with a reference image. This reference image may be the average of all non-shifted short-interval images. If the strength of the distribution is unknown, with the location presumed by the model geometry, e.g., a slat or a bracket, a generated image could also be used as reference for the cross-correlation. The cross-correlation can be normalized to put emphasis on the similarity than on the strength of the source. If a sound source spans across multiple frequencies in the acoustic spectrum, it is possible to cross-reference the image quality across all the measured frequencies. If the image quality is sufficient across most of the span this gives more confidence than if the data block only resulted in a good images for one particular frequency. For example the source is visible at 1, 3, 4 and 5 kHz and the image quality is found to be 0.9, 0.7, 0.9, and 0.9 respectively. The fact that 3 out of 4 frequencies agree on the high image quality gives more confidence than if only 1 out of 4 had a high image quality. This can be used as a second selection criteria. Mathematically the proposed method for qualification is summarized as follows: compute the (spatial) cross-correlation for each frequency bin using a Fourier transform: where is the cross-correlation, the short-interval image, and pad indicates the zero padding of the images. This computation is performed for each block of data independently. Subsequently, the cross-correlation is normalized with the auto-correlation: (C.6) The selection passes if the image qualities: summed together: are larger than a set criterion .

C.2.2. Non-stationary sources
If the strength of a noise source is varying with time, instantaneous images with a low response may still be of high quality. However, a bias towards short-intervals with high source power is found when images are selected by their maximum response or their cross-correlation with the reference image. To mitigate this issue a normalized cross-correlation provides an alternative. With the normalization only the relative source strength is compared and not the absolute source strength.

C.3. Stage III: shifting & combining
For a single source the approximate shift is obtained from the maximum source location. However, even with two point sources this type of shift determination is invalid. Any of the two sources can be the maximum response at any time. Thus the algorithm would place the source on the wrong location. A more robust method therefore is to use a cross-correlation to determine the shortinterval shift. For efficient cross-correlation a two-dimensional DFT is often used, which requires a regular two-dimensional grid to be used. The difference between two grid points determines the shifting resolution, which should be less than the beam width of the array, i.e., grid < 1. 22 2 , (C. 9) where is the distance between source and array, the speed of sound, the array radius, and the frequency [1].

Appendix D. Measurement of time travel variations
This section expands the method to compute the time variations proposed by Koop [12] and Sijtsma [23]. Assuming a sinusoidal input signal with Fourier transform and a measured signal with Fourier transform defined as: where is a window function,̄is a carrier angular frequency, is time, ′ are time variations caused by the acoustic waves traveling through the shear layer (dependent on time) and 0 is the angular frequency of the Fourier transform. For simplicity a generic hat window, also known as rectangular window, is chosen as with width . The Fourier transform of the input signal can now be evaluated as: Since the signal is only defined at the carrier frequency the frequency of interest is 0 . However, this requires resolving the singularity, which is done by taking the limit 0 →̄. ] e ī′ d (D.7) The limit is rewritten in terms of =̄− 0 . If is sufficiently small such that the exponent is almost constant in the short-interval corresponding to , the time variation follows as the argument of the exponential. The integral effectively works as a low-pass filter, thus the value of must be as small as allowed by the SNR such that a sufficiently high frequency response can be achieved. The main criteria here is that the dominant time variations, i.e., the ones of highest amplitude, are captured in the spectra. The SNR criterium is essential since smaller values of will increase the bin size allowing more noise to be present in the bin containing the carrier frequency. Since the power of the carrier remains constant a large bin size will capture more noise into the bin.

Appendix E. Beamforming maps
The results obtained using conventional beamforming and Acoustic Lucky Imaging are shown in Figs. E.25 and E.26.

Appendix F. Industrial scale wind tunnel: measurement setup
The experimental setup is shown in Fig. F.27.