Imaging through deep turbulence using single-shot digital holography data

Cameron J. Radosevich; Casey J. Pellizzari; Samuel Horst; Mark F. Spencer

doi:10.1364/OE.395674

1. Introduction

Despite multiple decades of outstanding research on deep turbulence, it remains an unsolved problem. This statement is in reference to one’s ability to resolve distant objects in the presence of deep-turbulence conditions. In general, deep-turbulence conditions arise from aberrations and the resulting phase errors that distribute themselves along the propagation path. Given the coherent illumination of optically rough extended objects, phenomena such as speckle and scintillation result, in addition to blurring and anisoplanatism, making for a multifaceted problem. This multifaceted problem has multiple empirically determined limitations and demonstrating performance in the presence of these limitations is the purpose of this work.

In this paper, we use an experimental setup consisting of phase plates to create distributed-volume phase errors. These phase plates allow us to replicate deep-turbulence conditions with path-integrated Kolmogorov statistics. In particular, we formulate four imaging scenarios with weak-to-deep turbulence strengths in a scaled-laboratory environment for the purpose of demonstrating performance in the presence of blurring, anisoplanatism, and scintillation. We parameterize the strength of these phenomena using the Fried parameter, the isoplanatic angle, and the Rytov number, respectively. In turn, this calibrated experimental setup allows us to validate the performance of our novel approach to imaging through deep turbulence.

This aforementioned approach starts with a digital-holography receiver [1–4]. After collecting single-shot digital holography data (i.e., one speckle measurement from the coherent illumination of an optically rough extended object), we then use a recently published algorithm referred to as multi-plane iterative reconstruction (MIR) [5–7]. In essence, the MIR algorithm allows us to sense and correct for the distributed-volume phase errors caused by deep-turbulence conditions.

Similar to model-based iterative reconstruction (MBIR) [8–10], the MIR algorithm uses a Bayesian framework. Such a framework allows us to estimate the focused, speckle-free image, $r$, in addition to the distributed-volume phase errors, $\phi = [\phi _1,\phi _2,\ldots ,\phi _{K}]$, where $K$ denotes the $k^{th}$ plane. This last point distinguishes the MIR algorithm from MBIR, which only senses and corrects for phase errors in the pupil plane of the digital-holography receiver. In particular, the MIR algorithm uses expectation maximization to jointly estimate $r$ and $\phi =[\phi _1,\phi _2,\ldots ,\phi _{K}]$. In its current form, the MIR algorithm also uses a Q-generalized Gaussian Markov random field to enforce correlation between neighboring pixels. This approach works well for simple objects (e.g., binary resolution charts), and ongoing efforts in denoising using a plug-in-play framework is extending the MIR algorithm to work well for complex objects [11–13].

With the details of the MIR algorithm in mind, this paper presents results from two MIR algorithm configurations, in addition to the pupil-plane-only corrections of MBIR. Similar to the approach taken in [5], the first configuration, MIR-A, assumes that we have a priori knowledge of the placement of the phase plates, so that we sense and correct in the exact locations of the phase errors. New to this paper, the second configuration, MIR-B, then assumes that we have no a priori knowledge of the placement of the phase plates, so that we sense and correct in two fixed planes for all phase-error combinations. Overall, the results show that both configurations outperform MBIR and that MIR-B performs comparably to MIR-A. Such results are promising for tactical applications, where one might not have a priori knowledge of the strength of the blurring, anisoplanatism, and scintillation caused by the deep-turbulence conditions.

Before moving on to the next section, the reader should note that this paper significantly builds upon the experimental results found in [5]. A shortcoming of the experimental results found in [5] is that the authors did not use a calibrated experimental setup with known turbulence strengths. Specifically, they used CD jewel cases with unknown path-integrated statistics. This shortcoming makes the experimental results presented here that much more meaningful, as our calibrated experimental setup allows us to validate the performance of the MIR algorithm in the presence of multiple empirically determined limitations, demonstrating that we have a novel approach to imaging through deep turbulence.

In what follows, we review the details of our experimental setup in Sec. 2. Then in Sec. 3., we review the details of the MIR algorithm. Section 4. provides the results and discussion for this paper, whereas Sec. 5. provides the conclusion.

2. Experimental setup

To create a calibrated experimental setup, we wanted to replicate deep-turbulence conditions with path-integrated Kolmogorov statistics. To do so, we distributed phase plates along the propagation path between our extended object and the lens to our digital-holography receiver. We then used our digital-holography receiver to collect single-shot digital holography data. Overall, we collected single-shot digital holography data with four distinct turbulence strengths. This data allowed us to validate the performance of the MIR algorithm. Thus, in what follows, we provide the details associated with our optical layout, weak-to-deep turbulence strengths, and path-integrated Kolmogorov statistics.

2.1 Optical layout

Figure 1(a) shows a schematic of the source board. Here, we split the light from a 700 mW continuous-wave, master-oscillator laser with a wavelength of 1064 nm (JDSU NPRO Series 126) into two single-mode, polarization-maintaining fibers through a Faraday isolator, lens, half-wave ($\lambda /2$) plate, polarizing beam splitter (PBS) cube, variable neutral density (VND) filters, and microscope objectives on 5-axis fiber couplers. These fibers formed the signal fiber and the local oscillator (LO) fiber, respectively.

Fig. 1. Overview of the experimental setup. Here, (a) describes the source board, which couples light into signal and LO fibers; (b) describes the entire experimental setup; and (c) shows the actual experimental setup with the extended object highlighted in orange.

Download Full Size | PDF

Figure 1(b) then shows a schematic of the entire experimental setup. Here, we directed the light from the signal fiber to a lens and a mirror to create a diverging beam that coherently illuminated the extended object [i.e., the chrome-on-glass 1951 US Air Force Bar Chart backed by Labsphere Spectralon highlighted in orange in Fig. 1(c)]. Note that by translating the Labsphere Spectralon relative to the stationary chrome-on-glass 1951 US Air Force Bar Chart, we achieved independent speckle realizations from measurement to measurement. Also note that we placed this extended object 2.50 m away from the lens of our digital-holography receiver, which was in the off-axis image plane recording geometry (IPRG) [1,3,14,15]. This lens had a focal length of 100 mm.

In accordance with the off-axis IPRG, we imaged the signal light scattered from the optically rough extended object onto the focal plane array (FPA) of our camera (Allied Vision Guppy PRO f-125 FireWire). For all the measurements, we used a 768x768 pixel region of interest on the FPA, which had a 3.75 $\mu$m pixel pitch. We placed the tip of the LO fiber next to the lens of our digital-holography receiver, so that it produced diverging-tilted reference light that interfered with our converging-aberrated signal light. Also in accordance with the off-axis IPRG, we measured the resulting interference pattern in an image plane, creating single-shot digital holography data.

By placing phase plates along the propagation path between the extended object and the lens of the digital-holography receiver, we replicated deep-turbulence conditions of varying turbulence strengths. Figure 1(c) shows the actual experimental setup with four phase plates, enabling an imaging scenario with deep-turbulence conditions. The next subsection provides further details on the weak-to-deep turbulence strengths used in this paper.

2.2 Weak-to-deep turbulence strengths

In order to achieve weak-to-deep deep turbulence strengths, we used four phase plates manufactured by Lexitek, Inc. [16]. Note that Lexitek, Inc. used random phase screens (generated from a Kolmogorov power-spectral density similar to the approach used by Schmidt [17]) to manufacture these phase plates with prescribed Fried parameter values of 0.55 mm, 0.85 mm, 1.20 mm, and 2.35 mm. Also note that using digital-holographic detection, we determined the percentage errors between the prescribed and measured Fried parameter values to within $\pm 20\%$ [18]. By distributing the phase plates along the propagation path, we were then able to achieve path-integrated values for the spherical-wave Fried parameter, $r_0$, using the following discrete summation [17]:

(1)$$r_0 = \left[ \sum _{i=1}^{K}r_{0_i}^{-\frac{5}{3}} \left( \frac{z_{i}}{\Delta z} \right) ^{\frac{5}{3}} \right] ^{-3/5},$$

where $r_{0_i}$ is the prescribed Fried parameter of an individual phase plate, $z_i$ is the location of the phase plate relative to the extended object, and $\Delta z$ is the object distance (i.e., the distance between the extended object and the lens of the digital-holography receiver). In addition, we were able to achieve path-integrated values for the isoplanatic angle, $\theta _0$, using the following discrete summation [17]:

(2)$$\begin{aligned} \theta_{0}= \left[ \frac{2.91}{0.423}\sum _{i=1}^{K}r_{0_i}^{-\frac{5}{3}} \left( \Delta z -z_i \right) ^{\frac{5}{3}} \right] ^{-3/5} , \end{aligned}$$

and the spherical-wave Rytov number, $\sigma _{\chi }^{2}$, using the following discrete summation [17]:

(3)$$\begin{aligned} \sigma_\chi^{2}=1.33k^{-5/6} \Delta z^{\frac{5}{6}} \sum _{i=1}^{K}r_{0_i}^{-\frac{5}{3}} \left( \frac{z_{i}}{\Delta z} \right) ^{\frac{5}{6}} \left( 1-\frac{z_{i}}{\Delta z} \right) ^{\frac{5}{6}}, \end{aligned}$$

where $k=2\pi /\lambda$ is the angular wavenumber and $\lambda$ is the wavelength.

Examining Eqs. (1)–(3), a discrete solution space exists for the spherical-wave Fried parameter, $r_{0}$, the isoplanatic angle, $\theta _{0}$, and the spherical-wave Rytov number, $\sigma _{\chi }^{2}$, given the phase plates with prescribed Fried parameter values, $r_{0_i}$, and a fixed object distance, $\Delta z$. Thus, we defined "targeted" values for weak, moderate (isoplanatic), moderate (anisoplanatic), and deep turbulence strengths. We then implemented a Matlab-based search algorithm to find the locations of each phase plate in order to obtain the closest physical match or "achieved" values. Table 1 shows these targeted values versus the achieved values.

Table 1. Summary of the weak-to-deep turbulence strengths used in this paper.

View Table | View all tables in this article

Notice that in Table 1, we parameterize the turbulence strengths in terms of $D/r_0$ and $\theta _0/(\lambda /D)$. These normalizations allow us to account for the aperture diameter, $D=11.1$ mm, of our digital-holography receiver. In general, by normalizing the aperture diameter by the spherical-wave Fried parameter (i.e., $D/r_0$), one can gauge the strength of the blurring. If $D/r_0>4$, then one obtains turbulence-limited images that require higher-order correction in order to track features on the extended object. The reader should note that this regime is an empirically determined limitation and helps to motivate the results presented in this paper.

Similarly, by normalizing the isoplanatic angle by the diffraction-limited half angle [i.e., ${{\tilde {\theta }}_{0}}=\theta _0/(\lambda /D)$], one can gauge the strength of the anisoplanatism. If ${{\tilde {\theta }}_{0}}>10$, then the imaging system is generally isoplanatic, such that the distributed-volume phase errors provide a linear, shift-invariant imaging system. In this regime, we find that the imaging system nominally has one point-spread function (PSF) across its full field of view (FOV). As ${{\tilde {\theta }}_{0}}\to 1$, however, the imaging system becomes increasingly anisoplanatic, such that the distributed-volume phase errors provide a linear, shift-varying imaging system. Here, we find that the imaging system has multiple PSFs across its full FOV. In the limit that ${{\tilde {\theta }}_{0}}\to 1$, it becomes increasingly more difficult to form a focused image, even with higher-order corrections in multiple planes. The reader should also note that these regimes are empirically determined limitations that help to motivate the results presented in this paper.

Lastly, in Table 1, the spherical-wave Rytov number $\sigma _\chi ^{2}$ gives us a gauge for the amount of scintillation. In this paper, the term "scintillation" refers to the constructive and destructive interference that results from distributed-volume phase errors and should not be confused with the term "speckle," which refers to the constructive and destructive interference that results from rough-surface scattering. Recall that the log-amplitude variance associated with a propagated point-source beacon, in practice, increasingly saturates when $0.1<\sigma _\chi ^{2}<1$ [19,20]. Generally, when $\sigma _\chi ^{2}>0.1$ and $\sigma _\chi ^{2}\to 0.5$, one is transitioning from the weak to moderately saturated scintillation regime, and when $\sigma _\chi ^{2}>0.5$ and $\sigma _\chi ^{2}\to 1$, one is transitioning from the moderately to fully saturated scintillation regime. In both regimes, the branch-point density tends to grow linearly without bound as a function of $\sigma _\chi ^{2}$ [19,20]. These regimes, yet again, refer to empirically determined limitations that help motivate the results presented in this paper.

Before moving on to the next subsection, it is important to note that Appendix A tabulates the locations of the phases plates needed to achieve the weak-to-deep turbulence strengths reported in Table 1. Recall that we determined the "targeted" strengths using the empirically determined limitations outlined above. Also recall that we calculated the "achieved" strengths using our Matlab-based search algorithm, which makes use of Eqs. (1)–(3).

2.3 Path-integrated Kolmogorov statistics

To verify that our phase-plate locations achieved the parameters reported in Table 1, we modified the experimental setup shown in Fig. 1. To this end, we placed the tip of the signal fiber in the object plane and centered it on axis. Note that the signal fiber approximated a point-source beacon. Also note that we imaged the light from this point-source beacon onto the FPA of our camera, interfering the converging-aberrated signal light with the diverging-tilted reference light. In accordance with the off-axis IPRG, we measured the resulting interference pattern in an image plane, creating multiple digital holograms for each measurement.

In total, we measured 12 unique digital holograms for the four turbulence strengths given Table 1, each with different phase-plate realizations. We achieved different phase-plate realizations by rotating the phase plates by 30$^{\circ }$, so that we exposed the light from the point-source beacon to different annular sections of the phase plates. To obtain estimates of the complex-optical field, we demodulated the digital holograms using the approach outlined in Appendix B. In accordance with the off-axis IPRG, these complex-optical field estimates were in the pupil plane of our digital holography receiver; thus, in the analysis that follows, we refer to these estimates as ${{\hat {U}}_{P}}(\mathbf {r})$.

To determine path-integrated values for the measured spherical-wave Fried parameter, $r_0'$, we first calculated the mutual coherence function, $\Gamma (\textbf {r}_1,\textbf {r}_2)$, from the complex-optical field estimates. In particular,

(4)$$\begin{aligned} \Gamma(\textbf{r}_1,\textbf{r}_2) = \langle {\hat{U}}_{P}(\textbf{r}_1){\hat{U}}^{*}_{P}(\textbf{r}_2) \rangle, \end{aligned}$$

where here, the angle brackets denote a correlation with respect to the complex-optical field estimates at two points in space. Normalizing $\Gamma (\textbf {r}_1,\textbf {r}_2)$, we obtain the coherence factor, $\mu (\textbf {r}_1,\textbf {r}_2)$, such that

(5)$$\begin{aligned} \mu ({{\mathbf{r}}_{1}},{{\mathbf{r}}_{2}})=\frac{|\Gamma ({{\mathbf{r}}_{1}},{{\mathbf{r}}_{2}})|}{\sqrt{\left| \Gamma ({{\mathbf{r}}_{1}},{{\mathbf{r}}_{1}}) \right|\left| \Gamma ({{\mathbf{r}}_{2}},{{\mathbf{r}}_{2}}) \right|}}, \end{aligned}$$

which we can rewrite in terms of the wave-structure function, $D(\textbf {r}_1,\textbf {r}_2)$, viz.

(6)$$\begin{aligned} \mu ({{\mathbf{r}}_{1}},{{\mathbf{r}}_{2}})=\textrm{exp}\left[-\frac{1}{2}D({{\mathbf{r}}_{1}},{{\mathbf{r}}_{2}}) \right]. \end{aligned}$$

According to the Rytov approximation [17], atmospheric turbulence has a wave-structure function of the following form:

(7)$$\begin{aligned} D(\Delta r)=6.88{\left( \frac{\Delta r}{{{r}_{0}'}} \right)}^{5/3}, \end{aligned}$$

where $\Delta r=\left | {{\mathbf {r}}_{2}}-{{\mathbf {r}}_{1}} \right |$. Equations (4)–(7) then yield a coherence factor, where

(8)$$\begin{aligned} \mu (\Delta r)=\textrm{exp}\left[ -3.44{{\left( \frac{\Delta r}{{{r}_{0}'}} \right)}^{5/3}} \right]. \end{aligned}$$

If we apply a least-squares fit to Eq. (8) using Matlab, then we can determine path-integrated values for $r_0'$. Figure 2 provides an example for the imaging scenario with weak-turbulence conditions (cf. Table 1).

Fig. 2. An example of how we determined the path-integrated values for the measured spherical-wave Fried parameter, $r_0'$. Here, (a) shows the average coherence factor for 12 phase-plate realizations for the imaging scenario with weak-turbulence conditions (cf. Table 1), and (b) shows the associated azimuthal average from the measured data (blue curve), in addition to the theoretical coherence factor (orange curve) and the least-squares fit (yellow curve).

Download Full Size | PDF

Table 2 presents (1) the path-integrated values for the theoretical spherical-wave Fried parameter, $r_0$ (obtained from our Matlab-based search algorithm), (2) the path-integrated values for the measured spherical-wave Fried parameter, $r_0'$ [obtained from our least-squares fits (cf. Fig. 2)], and (3) the associated percent errors. Since all the percent errors are <10%, we are confident that we accurately replicated path-integrated Kolmogorov statistics using our experimental setup.

Table 2. Path-integrated values for the theoretical and measured spherical-wave Fried parameters, $r_0$ and $r_0'$, respectively, and the associated percent errors. For completeness, we also present values for $D/r_0$ and $D/r_0'$, where $D$ is the aperture diameter. This normalization provides slightly higher percent errors.

View Table | View all tables in this article

3. MIR Algorithm

The MIR algorithm developed in [5–7] estimates the focused, speckle-free image, $r$, and the distributed-volume phase errors, $\phi = [\phi _1,\phi _2,\ldots ,\phi _{K}]$, from single-shot digital holography data. For a digital holography receiver using the off-axis IPRG [1,3,14,15], we model the complex data as

(9)$$\begin{aligned} y = & A_{\phi} g+w , \end{aligned}$$

where $g \in \mathbb {C}^{M}$ is the complex reflection coefficient of the object, $w \in \mathbb {C}^{M}$ is additive Gaussian measurement noise, and $A_{\phi } \in \mathbb {C}^{M\times M}$ is a linear transformation that accounts for the propagation geometry and is dependent on $\phi$ [5–7]. It is common to model $g$ as a complex Gaussian random variable with variance $E_g[|g|^{2} | r]=r$, where $E_g[\cdot |\cdot ]$ indicates the conditional expectation with respect to $g$ and $r$ is the real-valued reflectance of the scene [21].

Note that $r$ is the quantity we typically observe in incoherent images. Simply inverting the data to estimate $g$ results in images that are blurry, noisy, and speckled. Therefore, it is common to estimate $r$ by averaging estimates of $|g|^{2}$ from many independent speckle realizations and performing phase corrections as an independent step. Alternatively, MIR uses a Bayesian framework to jointly compute the maximum a posterior estimates of both $r$ and $\phi$ directly from a single realization of $y$. We obtain our estimates with an iterative optimization given by

(10)$$(\hat{r},\hat{\phi}) = \underset{r, \bar{\phi} } {\mathrm{argmin}} \left \{ - \textrm{E}_g \left[ \log p(y,g~ | r,\bar{\phi}) | y, r', \bar{\phi}' \right] - \log p \left( r \right) - \log p \left( \bar{\phi} \right)\right \},$$

where $\bar {\phi }$ is a low-resolution version of $\phi$, $p(\cdot )$ and $p(\cdot |\cdot )$ represent conventional and conditional probability distributions, respectively, and the primes, $r'$ and $\bar {\phi }'$, indicate the current estimate of the respective variables. Since the actual log-likelihood function is not easily evaluated, we use the expectation maximization to introduce $g$ into Eq. (10), simplifying its functional form. Then, we marginalize over $g$ using the expectation operator. We obtain a full resolution estimate of the phase errors according to $\phi = P(\bar {\phi })$, where $P(\cdot )$ is an 2D bi-linear interpolation function. Finally, we use Markov random fields for our prior models, $p(r)$ and $p(\bar {\phi })$. Full details on the models used in Eq. (10) and the optimization update steps are provided in [5–7].

The cost function that we minimize in Eq. (10) is nonconvex and therefore we usually converge to a local minimum that is not global. To ensure that we find a good local minima, in the sense that it produces a focused image, we iteratively change our initial conditions. At the start of the algorithm, we initialize the phase-error function with zeros. Next, we run the algorithm for a fixed number of steps to improve our estimates. We then restart the algorithm with the improved estimate of $\phi$ but reinitialize the reflectance according to $r=A_{\phi }^{H}y$. This bootstrapping processes is repeated a set number of times. We have found this initialization process to work well in the presence of very strong phase errors [5–7,9,10].

In this paper, we added a new feature to the MIR bootstrapping process. Similar to Thurman [22], we begin the bootstrapping process with $\bar {\phi }$ represented using a low-resolution grid. This low-resolution representation reduces the number of unknowns early in the estimation process and reduces variation in our cost function. After a set number of bootstrapping iterations, we increase the resolution of $\bar {\phi }$ by a factor of $2\times$ in each dimension. We repeat this process until we reach our final high-resolution representation of $\bar {\phi }$. As an example, in this paper, we began the bootstrapping process with $\bar {\phi }\in \mathbb {R}^{M/16^{2}}$, where $M=288\times 288$ was the total number of pixels, and we iteratively increased the resolution by $2\times$ in each dimension until $\bar {\phi }\in \mathbb {R}^{M/2^{2}}$. To avoid noise in the high-resolution estimates of $\bar {\phi }$, as we increase the resolution, we also decreased the variation, $\sigma _{\bar {\phi }}^{2}$, that we allow in our model of the phase errors, $p(\bar {\phi })$. We found this variable-resolution approach to produce superior results when compared to the MIR algorithm using a fixed resolution for $\phi$.

We also compared the MIR algorithm against MBIR [8–10]. In practice, MBIR is the precursor to the MIR algorithm; however, this approach was only designed to estimate isoplanatic phase errors that exist near the pupil plane of the digital-holography receiver. As a result, MBIR fails when the phase errors become anisoplanatic. Thus, we include MBIR reconstructions in our results to help visually highlight the degree of anisoplanatism of our data and to show the utility of the MIR algorithm. All MBIR parameters were set in accordance with [10].

For the MIR algorithm, we followed the parameter selection used in [5] with the following exceptions. We set the number of bootstrapping iterations to $N_L=30$ and the number of expectation maximization iterations per bootstrap iteration to $N_K=100$. During the final bootstrap iteration, we set the image regularization parameter to $\gamma = 2.2$ and set $T=0.01$, where $T$ controls how similar two neighboring pixels must be in value to have a strong correlation in $p(r)$.

We began the MIR algorithm with $n_b=16$ and $\sigma _{\bar {\phi }}=2$ rad. Here, $n_b$ is the down sampling factor used to represent $\phi$ on a low-resolution grid and $\sigma _{\bar {\phi }}$ is the standard deviation of $\bar {\phi }$ in our model, $p(\phi )$. Over the course of the bootstrapping process, we progressively reduced $n_b$ and $\sigma _{\bar {\phi }}$ by a factor of two until we ended up with $n_b=2$ and $\sigma _{\bar {\phi }}=0.25$ rad for the final estimate. We implemented these changes when the total bootstrapping process was $10\%$, $20\%$, $40\%$, and $60\%$, complete.

As an extension of the experimental results presented in [5], we tested two configurations of the MIR algorithm in this paper. For the first configuration, MIR-A, we estimated $\phi$ at the exact locations of the phase screens used in each experimental setup. In the second configuration, MIR-B, we only estimated $\phi$ at two planes that we fixed for all data sets at $z=1.80$ m and $z=2.40$ m. These values represent $72\%$ and $96\%$ of the total propagation distance and were heuristically chosen because they produced good results for all of the data sets. While the MIR-A configuration represents an idealized approach in which we have a priori knowledge of the distribution of the phase errors along the propagation path, the MIR-B configuration represents a more generalized approach that might be used in tactical applications.

4. Results and discussion

Figure 3 shows the resulting reconstructions with the strength of the turbulence increasing from left to right, beginning with the no-turbulence case. On the the top row, we show the original back-projected image, $r=A_{0}^{H}y$, where the subscript $0$ implies that $\phi =0$ in our linear transformation, $A_\phi$, that accounts for the propagation geometry. As such, this top row represents the original images obtained directly from the data with no additional processing. In addition to measurement noise and speckle, these images increasingly contain more and more blurring, anisoplanatism, and scintillation from left to right. The second, third, and fourth rows then show the MBIR, MIR-A, and MIR-B reconstructions, respectively.

Fig. 3. Reconstructions for the weak-to-deep turbulence strengths used in this paper (cf. Table 1) with the no-turbulence images in the far-left column. The top row shows the original image, which in addition to measurement noise and speckle, increasingly contains more and more blurring, anisoplanatism, and scintillation from left to right. Respectively, the remaining rows show the corresponding MBIR, MIR-A, and MIR-B reconstructions. These results show the utility of the MIR algorithm over MBIR in the presence of increasingly anisoplanatic conditions.

Download Full Size | PDF

The MBIR reconstructions show some noticeable residual blurring for the imaging scenario where $[D/r_0,~\tilde {\theta }_0,~\sigma _{\chi }^{2}]=[3.30,~10.3,~0.04]$ and significant residual blurring for the stronger-turbulence cases. From left to right, the residual blurring highlights the effects of increasingly anisoplanatic conditions and the inability to correct for distributed-volume phase errors when only estimating for $\phi$ in the pupil plane of the digital-holography receiver. Alternatively, both MIR algorithm configurations correct for the distributed-volume phase errors by estimating for $\phi$ at multiple planes along the propagation path. This multi-plane approach allows both MIR algorithm configurations to produce focused, speckle-free images in increasingly anisoplanatic conditions.

With the above information in mind, the MIR reconstructions work up to the imaging scenario where $[D/r_0,~\tilde {\theta }_0,~\sigma _{\chi }^{2}]=[10.2,~3.41,~0.27]$. While there is still some noticeable residual blurring in this case, the MIR-A and MIR-B reconstructions are significantly improved compared to the original image and the MBIR reconstructions. The MIR reconstructions unfortunately fail for the imaging scenario where $[D/r_0,~\tilde {\theta }_0,~\sigma _{\chi }^{2}]=[19.7,~1.53,~0.92]$; however, we can still observe some detail in the largest set of bars for both algorithm configurations. In this case, ${{\tilde {\theta }}_{0}}=\theta _0/(\lambda /D) \approx 1$ and MIR reconstructions succumb to the the effects of anisoplantism, since every diffraction-limited half angle within the full FOV of the digital-holography receiver is starting to experience a completely different PSF. This fundamental limitation with respect to anisoplanatism, coupled with the effects of blurring and scintillation, tends to hinder the effectiveness of the MIR reconstructions the most when in the presence of deep-turbulence conditions.

An interesting result shown in Fig. 3 is that the generalized MIR algorithm configuration, MIR-B, performs just as well as MIR-A. This is encouraging since it indicates that we can achieve high-quality estimates with no a priori knowledge of the distribution of the phase errors along the propagation path and with only two planes of correction.

5. Conclusion

In this paper, we used an experimental setup consisting of phase plates and a digital-holography receiver to validate the performance of a novel approach to imaging through deep turbulence. This approach makes use of single-shot digital holography data (i.e., one speckle measurement from the coherent illumination of an optically rough extended object), in addition to the MIR algorithm, which works by sensing and correcting for the distributed-volume phase errors caused by deep-turbulence conditions. As such, we first showed that our distributed-volume phase errors, created using the phase plates, followed path-integrated Kolmogorov statistics for weak-to-deep turbulence strengths. We then presented results from two MIR algorithm configurations: (A) where we had a priori knowledge of the placement of the phase plates, so that we corrected in the exact locations of the phase errors, and (B) where we did not have a priori knowledge of the placement of the phase plates, so that we corrected in two planes that we fix for all phase-error combinations. Given the weak-to-deep turbulence strengths created in our calibrated experimental setup, the results showed that the two MIR algorithm configurations performed comparably for the four imaging scenarios tested. Such results are promising for tactical applications, where one might not have a priori knowledge of the deep-turbulence conditions.

Appendix A

This Appendix tabulates the locations of the phase plates used to create the weak-to-deep turbulence strengths discussed in Section 2.2. In Table 3, the first column provides each phase plate’s prescribed Fried parameter $r_{0_i}$. The remaining columns then show the location of each plate for each imaging scenario. Note that we measured the distances with respect to the extended object, which we placed placed 2.50 m away from the lens of our digital holography receiver. Also note that only the imaging scenario with deep-turbulence conditions required the use of all four phase plates.

Table 3. Summary of the phase-plate locations needed for the weak-to-deep turbulence strengths used in this paper.

View Table | View all tables in this article

Recall that we used a Matlab-based search algorithm to obtain the phase-plate locations provided in Table 3. Referencing Eqs. (1)–(3), the free variable in the experimental setup is $z_i$, which is the location of the phase plates relative to the extended object. Therefore, we implemented a brute-force calculation for the following parameters: $D/r_0$, $\theta _0/(\lambda /D)$, and $\sigma _\chi ^{2}$ (cf. Table 1), varying $z_i$ for each individual phase plate along the propagation path, as well as for all possible combinations of the phase plates. For each turbulence strength listed in Table 1, we then chose the phase-plate locations that minimized the percent error between our targeted values verses our achieved values.

Appendix B

This Appendix outlines the procedures used to demodulate the digital holograms referenced in Section 2.3. Specifically, by averaging over 200 frames per phase-plate realization, we significantly improved the signal-to-noise ratio (SNR) of our complex-optical field estimates by an average factor of 8.4 (+ 9.0 dB). Figure 4 (top left) shows a frame-averaged digital hologram for a single phase realization for the imaging scenario with weak-turbulence conditions (cf. Table 1).

Fig. 4. An example of how we demodulated our digital holograms in order to obtain an estimate of the complex-optical field.

Download Full Size | PDF

We also collected 200 signal-only frames and 200 reference-only frames. Figure 4 (top left) shows the frame-averaged reference and signal irradiance patterns for a single, phase-plate realization, for the imaging scenario with weak-turbulence conditions given in Table 1. It is important to note that we subtracted these frame-averaged signal and reference irradiance patterns from the frame-averaged digital holograms to further improve the SNR by an average factor of 22 (+ 13.4 dB).

To obtain an estimate of the complex-optical field, we computed a 2D inverse fast Fourier transform (ifft) using Matlab from the frame-averaged data. Figure 4 shows the irradiance pattern in the Fourier plane (bottom right) and the associated amplitude and phase of the circularly filtered complex-optical field estimate (bottom left). In accordance with the off-axis IPRG, this complex-optical field estimate was in the pupil plane of our digital holography receiver.

Acknowledgments

The authors of this paper would like to thank the Joint Directed Energy Transition Office for sponsoring this research, and D. E. Thornton for many insightful discussions regarding the results presented within.

Disclosures

The authors declare no conflicts of interest.

References

1. M. F. Spencer, R. A. Raynor, M. T. Banet, and D. K. Marker, “Deep-turbulence wavefront sensing using digital-holographic detection in the off-axis image plane recording geometry,” Opt. Eng. 56(3), 031213 (2016). [CrossRef]

2. M. T. Banet, M. F. Spencer, and R. A. Raynor, “Digital-holographic detection in the off-axis pupil plane recording geometry for deep-turbulence wavefront sensing,” Appl. Opt. 57(3), 465–475 (2018). [CrossRef]

3. M. F. Spencer, “Spatial heterodyne,” Encyclopedia of Modern Optics, 2nd Edition (Elsevier, 2018).

4. D. E. Thornton, M. F. Spencer, and G. P. Perram, “Deep-turbulence wavefront sensing using digital holography in the on-axis phase shifting recording geometry with comparisons to the self-referencing interferometer,” Appl. Opt. 58(5), A179–A189 (2019). [CrossRef]

5. C. J. Pellizzari, M. F. Spencer, and C. A. Bouman, “Imaging through distributed-volume aberrations using single-shot digital holography,” J. Opt. Soc. Am. A 36(2), A20–A33 (2019). [CrossRef]

6. C. J. Pellizzari, M. F. Spencer, and C. A. Bouman Jr, “Single shot imaging using digital holography receiver,” (2019). US Patent 10,416,609.

7. C. J. Pellizzari, M. F. Spencer, and C. A. Bouman Jr, “Method of single shot imaging for correcting phase errors,” (2020). US Patent 10,591,871.

8. C. A. Bouman, “Model based image processing,” https://engineering.purdue.edu/bouman/publications/pdf/MBIP-book.pdf (2013).

9. C. J. Pellizzari, M. F. Spencer, and C. A. Bouman, “Phase-error estimation and image reconstruction from digital-holography data using a bayesian framework,” J. Opt. Soc. Am. A 34(9), 1659–1669 (2017). [CrossRef]

10. C. J. Pellizzari, M. T. Banet, M. F. Spencer, and C. A. Bouman, “Demonstration of single-shot digital holography using a bayesian framework,” J. Opt. Soc. Am. A 35(1), 103–107 (2018). [CrossRef]

11. G. T. Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman, “Plug-and-play unplugged: Optimization-free reconstruction using consensus equilibrium,” SIAM J. Imaging Sci. 11(3), 2001–2020 (2018). [CrossRef]

12. C. J. Pellizzari, M. F. Spencer, and C. A. Bouman, “Optically coherent image reconstruction in the presence of phase errors using advanced-prior models,” Proc. SPIE 10650, 106500B (2018). [CrossRef]

13. C. P. Pellizzari, M. F. Spencer, and C. A. Bouman, “Coherent plug-and-play: Digital holographic imaging through atmospheric turbulence using model-based iterative reconstruction and convolutional neural networks,” Submitted to IEEE Transactions on Computational Imaging (2020).

14. D. E. Thornton, M. F. Spencer, C. A. Rice, and G. P. Perram, “Digital holography efficiency measurements with excess noise,” Appl. Opt. 58(34), G19–G30 (2019). [CrossRef]

15. D. E. Thornton, D. Mao, M. F. Spencer, C. A. Rice, and G. P. Perram, “Digital holography experiments with degraded temporal coherence,” Opt. Eng. 59(10), 1 (2020). [CrossRef]

16. S. M. Ebstein, “Pseudo-random phase plates,” Proc. SPIE 4493, 150–155 (2002). [CrossRef]

17. J. D. Schmidt, Numerical Simulation of Optical Wave Propagation, with examples in Matlab (SPIE, 2010).

18. S. Horst, C. J. Radosevich, C. J. Pellizzari, and M. F. Spencer, “Measuring the fried parameter of transmissive phase screens using digital-holographic detection,” Proc. SPIE 11135, 111350D (2019). [CrossRef]

19. M. Spencer, “Wave-optics investigation of turbulence thermal blooming interaction: I. using steady-state simulations,” Opt. Eng. 59(8), 081804 (2020). [CrossRef]

20. M. F. Spencer, “Wave-optics investigation of turbulence thermal blooming interaction: Ii. using time-dependent simulations,” Opt. Eng. 59(8), 081805 (2020). [CrossRef]

21. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts & Company, 2006).

22. S. T. Thurman, “Phase-error correction in digital holography using single-shot data,” J. Opt. Soc. Am. A 36(12), D47–D61 (2019). [CrossRef]

Turbulence strengths	$D / r_{0}$ (targeted / achieved)	$θ_{0} / (λ / D)$ (targeted / achieved)	$σ_{χ}^{2}$ (targeted / achieved)
Weak	5.00 / 3.30	10.0 / 10.3	0.10 / 0.04
Moderate (isoplanatic)	10.0 / 9.82	10.0 / 10.0	0.10 / 0.08
Moderate (anisoplanatic)	10.0 / 10.2	5.00 / 3.41	0.50 / 0.27
Deep	20.0 / 19.7	1.00 / 1.53	1.00 / 0.92

	$r_{0}$ (mm) theoretical	$r_{0}^{'}$ (mm) measured	percent error	$D / r_{0}$ theoretical	$D / r_{0}^{'}$ measured	percent error
Weak	3.36	3.16	−6.17%	3.30	3.52	6.58%
Moderate (isoplanatic)	1.13	1.10	−2.96%	9.82	10.12	3.05%
Moderate (anisoplanatic)	1.04	1.07	3.11%	10.17	10.34	−3.16%
Deep	0.56	0.51	−9.77%	19.6	21.8	10.8%

$r_{0_{i}}$	Weak	Moderate (isoplanatic)	Moderate (anisoplanatic)	Deep
0.55 mm	Not used	Not used	Not used	1.69 m
0.85 mm	Not used	Not used	1.84 m	1.75 m
1.2 mm	Not used	2.37 m	Not used	1.62 m
2.35 mm	1.75 m	1.80 m	1.80 m	1.54 m

Turbulence strengths	$D / r_{0}$ (targeted / achieved)	$θ_{0} / (λ / D)$ (targeted / achieved)	$σ_{χ}^{2}$ (targeted / achieved)
Weak	5.00 / 3.30	10.0 / 10.3	0.10 / 0.04
Moderate (isoplanatic)	10.0 / 9.82	10.0 / 10.0	0.10 / 0.08
Moderate (anisoplanatic)	10.0 / 10.2	5.00 / 3.41	0.50 / 0.27
Deep	20.0 / 19.7	1.00 / 1.53	1.00 / 0.92

	$r_{0}$ (mm) theoretical	$r_{0}^{'}$ (mm) measured	percent error	$D / r_{0}$ theoretical	$D / r_{0}^{'}$ measured	percent error
Weak	3.36	3.16	−6.17%	3.30	3.52	6.58%
Moderate (isoplanatic)	1.13	1.10	−2.96%	9.82	10.12	3.05%
Moderate (anisoplanatic)	1.04	1.07	3.11%	10.17	10.34	−3.16%
Deep	0.56	0.51	−9.77%	19.6	21.8	10.8%

Imaging through deep turbulence using single-shot digital holography data

Abstract

1. Introduction

2. Experimental setup

2.1 Optical layout

2.2 Weak-to-deep turbulence strengths

2.3 Path-integrated Kolmogorov statistics

3. MIR Algorithm

4. Results and discussion

5. Conclusion

Appendix A

Appendix B

Acknowledgments

Disclosures

References

Cited By

Figures (4)

Tables (3)

Equations (10)

Optics Express