Lightweight starshade position sensing with convolutional neural networks and simulation-based inference

Starshades have the potential to discover and characterize the atmospheres of Earth-like exoplanets in the habitable zone of nearby stars[1, 2, 3]. Their ability to achieve high contrast while maintaining high optical throughput and broad wavelength coverage make them the most promising technology to produce the first spectrum of an exo-Earth atmosphere. Starshades are currently in the technology development phase, with a targeted effort to advance critical starshade technologies (such as formation sensing and control) to Technology Readiness Level (TRL) 5[4]. A number of mission concepts have been proposed, including: the Starshade Rendezvous Mission[2] (SRM), which pairs a 26 meter diameter starshade with the 2.4 meter diameter Nancy Grace Roman Space Telescope (NGRST; formerly known as WFIRST) separated by 27,000 km; and the Habitable Exoplanet Explorer[3] (HabEx), which pairs a 52 meter diameter starshade with a 4 meter diameter telescope separated by 72,000 km. In this work, we use SRM as a case study to demonstrate our methodology, but it is applicable to a starshade mission of any size.

To maintain high contrast during observations, the telescope must remain in the deep shadow cast by the starshade; typically 2 meters in diameter larger than the size of the primary mirror. The size of the shadow can be designed to be arbitrarily large, but at the expense of a larger starshade farther away, which takes more fuel and time to switch between targets. Therefore, improved formation flying allows for a more efficient starshade architecture. It has long been recognized that the most difficult aspect to starshade formation flying is accurately sensing the relative position between the starshade and telescope over thousands of kilometers separation. Knowledge of the lateral position of the starshade relative to the telescope’s line of sight to the target star should be accurate to tens of centimeters in order to maintain position inside the typical $\pm 1$ meter tolerance. For instance, both SRM and HabEx set a lateral position sensing requirement of 30 cm $(3\sigma)$. Provided accurate sensing information is available, controlling the starshade to maintain $\pm 1$ meter alignment is relatively straightforward: the starshade observatory operates at the Sun-Earth Lagrange point (L2) where the dominant disturbances, gravity gradients, and solar radiation pressure are relatively benign and easily handled with linear control[5, 6, 7].

Sensing the starshade’s relative position to centimeter-level accuracy over 27,000 km is made possible by exploiting the nature of diffraction. At wavelengths outside of the starshade’s designed bandpass (“out-of-band light”), its suppression performance rapidly degrades and at the center of the shadow, a bright signal of diffracted starlight reemerges as the spot of Arago (see Figure 1). The position of this diffraction peak corresponds exactly with the position of the starshade. By imaging the diffraction pattern at the entrance pupil of the telescope, the starshade’s lateral position can thus be determined to high accuracy[8].

The formation sensing problem is now a matter of quickly, accurately, and reliably extracting the position of the spot of Arago from pupil plane images. Previous methods¹¹1Neural networks were also used in a very recent work[10]. that have demonstrated sufficient accuracy in this task include an image matching algorithm[11] and a non-linear least squares (NLLS) fit to a Bessel function model[7]. The image matching algorithm is slow to compute, which could limit high cadence updates to the spacecraft control. It also requires a large image library, which would occupy the limited memory of the spacecraft’s computer, and is limited by the finite resolution of the image library. The fit to the Bessel model is computationally intensive for state-of-the-art flight computers and requires an accurate initial guess to converge to the correct solution[12].

To overcome the limitations of these previous methods, we propose a convolutional neural network (CNN) to extract the starshade position from pupil plane images. In this paper we demonstrate through simulation and experiment that the CNN achieves accurate and robust position sensing with substantially reduced computing power and memory compared to previous methods. We then make use of the CNN to produce a summary statistic for a computationally efficient simulation-based inference (SBI) technique that yields an improved estimate of the starshade position and its uncertainty. We validate this method experimentally in the Princeton Starshade Testbed[9], where 10^-10 contrast has been demonstrated on 1/1000${}^{\text{th}}$ scale starshades at a flight-like Fresnel number.

The remainder of this paper is organized as follows. In Section 2 we pose the position sensing problem and describe the CNN architecture and training strategy. We show results on simulated and experimental data in Section 3 and Section 4, respectively. We further augment the CNN approach with a SBI technique in Section 5, and discuss our findings and future directions of this work in Section 6.

Imaging the telescope’s entrance pupil in out-of-band light shows the starshade’s diffraction pattern and provides a bright signal (out-of-band light is only suppressed by a factor of 10³, compared to 10¹⁰ with in-band light) indicating the starshade position. As an example, Figure 1 shows a simulated pupil image of NGRST including obstructions from the secondary mirror and supports.

Telescopes equipped with a coronagraph typically have a camera that takes images of the pupil plane to allow low-order wavefront sensing for the coronagraph. Using pupil images therefore provides positional information without the need for additional hardware.

We first investigate the performance of the CNN to estimate the starshade position. We run a baseline simulation with 96$\times$96 pixels sampling the telescope’s entrance aperture and a Peak SNR of 5, which corresponds to a 1 second exposure time of a 4${}^{\text{th}}$ magnitude star with NGRST. We simulate 10,000 images with random starshade positions distributed over a 3 meter square. Figure 3 shows the distribution of the $(x,y)$ and magnitude $R$ position errors for the baseline simulation. The mean and standard deviation of the magnitude of the error are 6.3 cm and 4.8 cm, respectively; and 99.7% of the errors are less than 25 cm. Figure 4 plots these errors as a function of the starshade’s position. The error increases when the spot of Arago is either obscured by the secondary mirror in the center or off the edge of the aperture. The latter situation is responsible for most of the tail of the error distribution. We thus consider 6 cm to be a conservative estimate because with an active control system the starshade should never stray past $\pm 1$ meter offset during observations. If we only consider times when the offset is less than 1 meter, the mean error drops to 2.8 cm, and 99.7% are less than 10.5 cm.

It is promising to see the good performance outside of the telescope’s aperture, which suggests this method could be extended farther outside the $\pm 1$ meter control box and could be used earlier in the transition between acquisition and observation guiding modes[12]. We save the investigation into how far this method can be extended for a later date.

To examine the robustness of the CNN solution over the width of the shadow, we run a large number of simulations with the starshade positioned on a 20 x 20 grid over a 3 meter width. Each point on the grid, 5,000 simulations are run with the Peak SNR randomly pulled from a uniform distribution. Figure 5 shows the result of this Monte Carlo simulation: the ellipses show the $3\sigma$ covariance of predicted positions centered on the mean prediction; the grid of true starshade positions are marked with black plus signs. Inside the telescope aperture, the covariance is on the order of a few centimeters and radially correlated. The left panel of Figure 5 is the moderate SNR case with the Peak SNR drawn from 5 - 15; the right panel is the high SNR case with Peak SNR drawn from 15 - 25. There is a slight bias in the predicted positions towards the edge of the aperture, but generally this is smaller than the $3\sigma$ level. Outside the aperture, there are larger covariances and biases, with a strong inward bias at the corners. The higher Peak SNR reduces the covariances, but has little to no effect on the biases.

Figure 6 shows the distribution of position errors as a function of Peak SNR. For each Peak SNR, 1000 images were simulated with randomly distributed positions; two curves are shown which limit the testset to offsets of a certain radius. We find that above a Peak SNR of 5, there is little change in the position error. Limiting to offsets of 1 meter results in a 50% decrease in position error.

Finally, we examine the trend in performance with the number of pixels sampling the pupil. Figure 7 shows the distribution of position errors for CNNs trained on images with sizes 48$\times$48, 64$\times$64, and 96$\times$96. We find that the position error scales approximately linearly with the number of pixels sampling the pupil (in one dimension); the mean position error for 48, 64, and 96 pixels is 12.4 cm, 10.5 cm, and 6.3 cm, respectively. As the CNN is inherently lightweight and computationally fast, the cost of a larger image size is minimal compared to the image matching and model fitting methods, and the gain in performance is worth the cost. While the choice of pupil sampling size will ultimately depend on the optical layout of the pupil imaging system and its detector electronics, the CNN benefits from a larger image size.

The Princeton Starshade Testbed[9] was built for optical model validation experiments at high contrast levels with a flight-like Fresnel number, i.e. the parameter governing optical diffraction, but at only $1/1000$ of the scale of a real starshade. Properties of the lab and two flight configurations are listed in Table 1. The sub-scale starshades are lithographically etched from a silicon wafer and characterized well enough to predict the starshade performance at 10^-9 contrast. As such, the optical model should be very accurate in predicting the out-of-band diffracted light at the 10^-3 contrast level.

The testbed provides the unique opportunity to demonstrate that our approach can be trained on simulated data and then accurately be applied to experimental images. There are, however, a number of factors unique to the testbed that make the experimental results a conservative estimate of the on-orbit performance. First, the testbed operates in the atmosphere, and changes in the air temperature and density result in motion of the diffraction pattern and variation in the optical power incident on the starshade. Second, a low level background is observed from light scattering off aerosols as it propagates the 80 meters from laser to telescope. Third, the pupil imaging optics are off-the-shelf and not calibrated to the extent that the flight system will be, so any model mismatch seen in the lab is dominated by the pupil imaging optics and not the starshade. Lastly, the flight missions propose using an electron-multiplying CCD, while the lab detector (properties listed in Table 2) is operated with a conventional amplifier and thus experiences a larger contribution from read-out noise.

To correct for atmosphere-induced motions in the diffraction pattern and variations in the brightness, the true position and brightness of each experimental image is solved for via the NLLS algorithm [7]. Each image is normalized by the NLLS-solved brightness, and the NLLS-solved position is used as the true position when determining the accuracy of the CNN-DELFI position.

Figure 8 shows the distribution of the X, Y, and magnitude (R) position errors for 23,700 images with Peak SNR = 5 - 8. The mean and standard deviation of the magnitude of the error are 6.5 cm and 4.0 cm, respectively, and 99.7% of the errors are less than 25 cm. Figure 9 plots these errors as a function of the starshade’s position. The error follows the same structure as that of the simulation data, with the worst performance in the center and towards the edges, when the spot of Arago is blocked by the secondary mirror and moves off the edge of the aperture.

Sampled over the same region, the positional accuracy for experimental images is worse by about a factor of 2 compared to simulated images, likely due to the aforementioned issues with testing in the atmosphere, but the performance still exceeds our goal of 30 cm.

The CNN produces an accurate point-estimate of the starshade position. Although it is entirely usable on its own, as we have demonstrated in the previous section, it lacks uncertainty estimates. These would be beneficial for the control system during formation flying because the UKF could then appropriately reduce the importance of the CNN predictions in regions of large uncertainty. We therefore couple the CNN with a subsequent calibration procedure, which provides a quantification of its associated uncertainty and yields an improved position estimate.

In keeping with the spirit of this work, we choose the lightweight method DELFI[21]. In short, it amounts to fitting a Gaussian mixture model (GMM) to pairs of true $(x,y)$ and CNN-estimated $(x^{\prime},y^{\prime})$ positions obtained from the same simulations we introduced in Section 2. With a pretrained CNN, this first step is done once and the resulting GMM is stored. At test time, the GMM is evaluated at the position $(x^{\prime},y^{\prime})$ as reported by the CNN, which yields a probability distribution $p(x,y\mid x^{\prime},y^{\prime})$. The mode of this distribution is the most likely true position, while the standard deviation provides an estimate of the uncertainty. Because of the mathematical structure of the GMM, the calculation at test time can be carried out analytically and require only matrix and vector operations in $\mathbb{R}^{2}$. We explain the approach in detail below.

The results of Section 3 and Section 4 show that the CNN alone can easily meet the position estimation requirement of 30 cm set in mission designs[2, 3]. For the entirety of the observation phase, the control law should not allow the misalignment to exceed 1 meter, therefore we focus on results where the starshade position is less than 1 meter. When doing so, the 99.7 percentile position error for the baseline simulation is 10.5 cm, and for the experimental results is 21.5 cm.³³3The experimental results should be viewed as a conservative estimate of the performance as the mismatch between experimental and simulated training data is dominated by the effects of the atmosphere and uncharacterized pupil imaging optics. For an actual mission, atmospheric issues will not be present in space, and the imaging system will be better characterized on the ground and after launch for an accurate model of the instrument response. These results meet the 30 cm requirement that defines the formation sensing milestone for TRL 5[4]. As we demonstrate in Section 5, augmenting the CNN with our simulation-based calibration method DELFI further improves the estimation accuracy and pushed out typical errors of less than 10 cm to well beyond the telescope pupil.

Our results compare favorably to alternative approaches with image matching[11] and NLLS[7] algorithms, all of which have been shown to satisfy accuracy requirements. The most significant advantages of our approach is that it is lightweight in memory and fast to compute. The CNN trained on images of 96$\times$96 pixels has $4\times 10^{5}$ parameters and takes up 1.6 MB in disk space. A CNN trained on 32$\times$32 images only needs 92 kB of disk space, 500$\times$ less than the 46 MB required for the image library (32$\times$32) of the matching algorithm[11]. A DELFI model with $K=250$ Gaussian components also only occupies 42 kB. The CNN needs 5.3 MFLOPs (million floating point operations) to process a 96$\times$96 image; 0.36 MFLOPs to process a 32$\times$32 image. The computational cost of DELFI is entirely negligible. For the image library method, Ref. 11 reports 70 MFLOPs to search the whole library, but suggest that a reduced grid search around an initial guess could reduce that to 0.5 MFLOPs. The NLLS algorithm involves computing Equation 6 and its analytical Jacobian at each pixel once per iteration; with a good initial guess, the Levenberg-Marquardt solution typically converges within 15 iterations. The NLLS computation time is dominated by the calculation of the Bessel function and its derivative ($J_{0}$, $J_{1}$), which we estimate to require $3BN$ FLOPs each ($B$ is the number of terms kept in the Bessel function series expansion and $N$ is the total number of pixels). Therefore, the NLLS takes 14 MFLOPs to process a 96$\times$96 image; 1.6 MFLOPs to process a 32$\times$32 image.

Our method thus substantially reduces storage space and computation time over the image matching algorithm. It also produces a continuous position estimate, whereas the image matching results are discrete due to the finite offset sampling in the image library. Our method also outperforms the NLLS algorithm, though to a lesser degree, and provides the following extra benefits: 1) unlike NLLS, it does not require a valid initial guess to properly converge, and 2) the Bessel function model used by the NNLS is only valid out to an offset of 1.5 meters. Future work can investigate whether additional training of the CNN model and DELFI renders it effective for larger offsets.

One important aspect of practical applicability we did not address in this work is related to variations in the stellar spectrum. Due to the monochromatic laser in the experimental testbed, this study was limited to light of a single wavelength and therefore does not capture the effect of the target star’s spectral type $t$. While we expect that CNN estimates to only weakly depend on $t$,⁴⁴4Over the 100 nm wide guiding bandpass, the diffraction pattern is well-approximated by the Bessel function model: $$I(s)=J^{2}_{0}\left(\frac{2\pi Rs}{\lambda z}\right)\,,$$ (6) where $I$ is the intensity as a function radial coordinate $s$, $R$ is the starshade radius, $z$ is the telescope–starshade separation, and $\lambda$ is the wavelength of light. Inverting Equation 6, the width of the spot of Arago scales linearly with wavelength (as $\lambda z/R$). To first order, and due to the stellar flux increasing with wavelength on the Wien’s tail for solar-type stars, the diffraction pattern will be dominated by the spot of Arago broadened by the longest wavelength in the bandpass. As the CNN is trained on units of suppression and images are normalized by the unblocked star’s expected flux, the CNN is expected to be weakly dependent on spectral type $t$.future works should be able to improve the accuracy for broadband images by training the CNN on image with a range of spectral types. For predictions it can then either implicitly marginalize over an unknown spectral type or, even better, condition its estimates on the known spectral type $t$ of the target star. The latter option can be achieved by adding the spectral type as an extra variable to the last fully connected layer of the CNN, at a minor increase of required memory. Our calibration method DELFI can also take the spectral type into account by evaluating $p(x,y\mid x^{\prime},y^{\prime},t)$, i.e. by modestly increasing the task of density estimation from a four- to a five-dimensional space.

We therefore conclude that position sensing with CNN and DELFI is a viable, accurate, and computationally efficient option for formation flying of future Starshade missions. To help with the adoption of our method, we publicly release the code at https://github.com/astro-data-lab/starshade-xy.