Joint 3D Localization and Classification of Space Debris using a Multispectral Rotating Point Spread Function

TL;DR

This paper presents a novel three-stage method for joint 3D localization and material classification of space debris using multispectral rotating point spread functions, improving accuracy over single-band approaches.

Contribution

The authors introduce a new three-stage approach combining optimization, spectral estimation, and classification for space debris analysis using multispectral RPSF data.

Findings

Enhanced localization accuracy with multispectral data

Improved material classification through spectral analysis

Method outperforms single-band techniques in numerical tests

Abstract

We consider the problem of joint three-dimensional (3D) localization and material classification of unresolved space debris using a multispectral rotating point spread function (RPSF). The use of RPSF allows one to estimate the 3D locations of point sources from their rotated images acquired by a single 2D sensor array, since the amount of rotation of each source image about its x, y location depends on its axial distance z. Using multi-spectral images, with one RPSF per spectral band, we are able not only to localize the 3D positions of the space debris but also classify their material composition. We propose a three-stage method for achieving joint localization and classification. In Stage 1, we adopt an optimization scheme for localization in which the spectral signature of each material is assumed to be uniform, which significantly improves efficiency and yields better localization results than possible with a single spectral band. In Stage 2, we estimate the spectral signature and refine the localization result via an alternating approach. We process classification in the final stage. Both Poisson noise and Gaussian noise models are considered, and the implementation of each is discussed. Numerical tests using multispectral data from NASA show the efficiency of our three-stage approach and illustrate the improvement of point source localization and spectral classification from using multiple bands over a single band.

Tables

Table 1. Table 1: Comparisons of 3D localization and classification by single-band and multispectral images.

Localization	No. bands		Localization		Classification
method	Stage 1	Stage 2	Recall	Precision	OA	kappa
Single	1	-	87.47%	39.22%	-	-
	1	5	86.00%	83.95%	88.09%	84.81%
Multiple	4	-	95.20%	63.78%	-	-
	4	5	94.80%	95.34%	93.28%	91.57%

Table 2. Table 2: Comparisons using class accuracy and a confusion matrix.

Localization method	Confusion matrix
Single	96.00	0.00	1.66	1.66	1.33
	0.00	88.00	0.00	0.00	0.67
	0.00	0.00	95.67	0.00	4.00
	2.00	24.67	3.33	74.33	6.33
	4.67	8.00	2.00	7.00	82.33
Multiple	96.00	0.00	0.66	3.00	0.66
	0.00	98.67	0.00	1.00	0.67
	0.00	0.00	93.33	0.00	0.66
	1.33	7.00	4.67	87.67	2.67
	4.00	4.67	0.00	6.33	86.67

Table 3. Table 3: Comparisons of 3D localization and classification accuracies with changing threshold.

$\gamma$	Localization		Classification
	Recall	Precision	OA	kappa
0.1	94.80%	90.45%	93.28%	91.57%
0.2	94.80%	95.34%	93.28%	91.57%
0.3	91.33%	96.61%	91.55%	89.17%
0.4	82.00%	97.06%	89.93%	86.17%
0.5	72.53%	97.82%	85.08%	79.48%

Equations

$$G(\mathbf{s})=\mathcal{N}\left(\sum _{j=1}^M{\mathcal{H}}_j(\mathbf{s}-{\mathbf{s}}_j)f_j+b\right),$$

$$\mathbf{s}=\frac{R}{\lambda z_I}\mathbf{r}$$

$${\mathcal{H}}_j(\mathbf{s})=\frac{1}{\pi }{\left|\int P(\mathbf{u})\exp \left[\iota (2\pi \mathbf{u}\cdot \mathbf{s}+{\zeta }_j{u}^2-\psi (\mathbf{u}))\right]d\mathbf{u}\right|}^2,$$

$${\zeta }_j=\frac{\pi (l_0-z_j)R^2}{\lambda l_0{z}_j}.$$

$$\psi (\mathbf{u})=l{\varphi }_{\mathbf{u}},\text{for }\sqrt{\frac{l-1}{L}}\le u\le \sqrt{\frac{l}{L}},l=1,\cdot \cdot \cdot ,L,$$

$$z_j=\frac{\pi l_0{R}^2}{\lambda l_0{\zeta }_j+\pi R^2}.$$

$$\left(-\frac{z_j{x}_j}{z_I},-\frac{z_j{y}_j}{z_I},z_j\right),$$

$$\{(x_j,y_j,{\zeta }_j,{\mathbf{f}}_j),j=1,2,\dots ,M\},$$

$$G^{(i)}({\mathbf{s}}^{(i)})=\mathcal{N}\left(\sum _{j=1}^M{\mathcal{H}}_j^{(i)}({\mathbf{s}}^{(i)}-{\mathbf{s}}_j^{(i)})f_j^{(i)}+b\right),$$

$${\mathbf{s}}^{(i)}=\frac{R}{{\lambda }^{(i)}{z}_I}\mathbf{s},{\mathbf{s}}_j^{(i)}=\frac{R}{{\lambda }^{(i)}{z}_I}{\mathbf{r}}_j,\mathrm{and}{\zeta }_j^{(i)}=\frac{\pi (l_0-z_j)R^2}{{\lambda }^{(i)}{l}_0{z}_j}$$

$$G=\mathcal{N}\left(\mathcal{T}(\mathcal{A}\ast \mathcal{X})+b1\right),$$

$$D_s(\mathcal{X},\mathcal{A},G)=⟨1,\mathcal{T}(\mathcal{A}\ast \mathcal{X})-G\log (\mathcal{T}(\mathcal{A}\ast \mathcal{X})+b1)⟩,$$

$$\underset{\mathcal{X}\ge 0}{\min }\left\{D_s(\mathcal{X},\mathcal{A},G)+\mu \sum _{p,q,r=1}^{m,n,d}\frac{|{\mathcal{X}}_{pqr}|}{a+|{\mathcal{X}}_{pqr}|}\right\},$$

$$\{\begin{array}{ll}{w}_{pqr}^k& =\frac{a\mu }{{\left(a+{\hat{\mathcal{X}}}_{pqr}^k\right)}^2},\forall p,q,r;\\ {\hat{\mathcal{X}}}^k& =\arg \underset{\mathcal{X}\ge 0}{\min }\left\{D_s(\mathcal{X},\mathcal{A},G)+\sum _{p,q,r=1}^{m,n,d}{w}_{pqr}^k|{\mathcal{X}}_{pqr}|\right\}.\end{array}$$

$$\begin{array}{rl}& ⟨1,\mathcal{T}{\mathcal{U}}_0-G\log (\mathcal{T}{\mathcal{U}}_0+b1)⟩+\sum _{p,q,r=1}^{m,n,d}{w}_{pqr}^k|({\mathcal{U}}_1{)}_{pqr}|\\ +& \frac{{\beta }_0}{2}\Vert {\mathcal{U}}_0-\mathcal{A}\ast \mathcal{X}-{\eta }_0{\Vert }^2+\frac{{\beta }_1}{2}\Vert {\mathcal{U}}_1-\mathcal{X}-{\eta }_1{\Vert }^2+I_{+}({\mathcal{U}}_1),\end{array}$$

$$I_{+}(\mathcal{X})=\{\begin{array}{ll}0,& \mathcal{X}\ge 0,\\ \infty ,& \text{otherwise}.\end{array}$$

$${\mathcal{U}}_0^{t+1}$$

$${\mathcal{U}}_1^{t+1}$$

$${\mathcal{X}}^{t+1}$$

$${\eta }_0^{t+1}$$

$${\eta }_1^{t+1}$$

$${\left({\mathcal{U}}_0^{t+1}\right)}_{pqr}=\{\begin{array}{ll}\frac{-({\xi }_1^t{)}_{pqr}+\sqrt{({\xi }_1^t{)}_{pqr}^2+4{\beta }_0{G}_{pq}}}{2{\beta }_0},& \text{if }r=d,\\ ({\xi }_0^t{)}_{pqr},& \text{otherwise}.\end{array}$$

$${\left({\mathcal{U}}_1^{t+1}\right)}_{pqr}=\max \left\{{\mathcal{X}}^t+{\eta }_1^t-w_{pqr}^k/{\beta }_1,0\right\}.$$

$${\mathcal{F}}^{-1}\left\{{\Omega }_s\cdot \left(\overline{\mathcal{F}\{\mathcal{A}\}}\cdot \mathcal{F}\{{\mathcal{U}}_0^{t+1}-{\eta }_0^t\}+\frac{{\beta }_1}{{\beta }_0}\mathcal{F}\{{\mathcal{U}}_1^{t+1}-{\eta }_1^t\}\right)\right\},$$

$$G^{(i)}\approx \mathcal{N}\left(\mathcal{T}({\mathcal{A}}^{(i)}\ast {\mathcal{X}}^{(i)})+b1\right),i=1,2,\dots ,K,$$

$$G^{(i)}\approx \mathcal{N}\left(\mathcal{T}({\mathcal{A}}^{(i)}\ast \mathcal{X})+b1\right).$$

$$\underset{\mathcal{X}\ge 0}{\min }\left\{\sum _{i=1}^K{D}_s(\mathcal{X},{\mathcal{A}}^{(i)},G^{(i)})+\mu \sum _{p,q,r=1}^{m,n,d}\frac{|{\mathcal{X}}_{pqr}|}{a+|{\mathcal{X}}_{pqr}|}\right\},$$

$$\{\begin{array}{ll}{w}_{pqr}^k& =\frac{a\mu }{{\left(a+{\hat{\mathcal{X}}}_{pqr}^k\right)}^2},\forall p,q,r;\\ {\hat{\mathcal{X}}}^k& =\arg \underset{\mathcal{X}\ge 0}{\min }\left\{\sum _{i=1}^K{D}_s(\mathcal{X},{\mathcal{A}}^{(i)},G^{(i)})+\sum _{p,q,r=1}^{m,n,d}{w}_{pqr}^k|{\mathcal{X}}_{pqr}|\right\}.\end{array}$$

$$\begin{array}{rl}& \sum _{i=1}^K⟨1,\mathcal{T}{\mathcal{U}}_0^{(i)}-G^{(i)}\log (\mathcal{T}{\mathcal{U}}_0^{(i)}+b1)⟩+\sum _{p,q,r=1}^{m,n,d}{w}_{pqr}^k|({\mathcal{U}}_1{)}_{pqr}|\\ +& \frac{{\beta }_0}{2}\sum _{i=1}^K\Vert {\mathcal{U}}_0-{\mathcal{A}}^{(i)}\ast \mathcal{X}-{\eta }_0{\Vert }^2+\frac{{\beta }_1}{2}\Vert {\mathcal{U}}_1-\mathcal{X}-{\eta }_1{\Vert }^2+I_{+}({\mathcal{U}}_1).\end{array}$$

$$\{\begin{array}{ll}{\mathcal{U}}_0^{t+1}& =\underset{{\mathcal{U}}_0}{\mathrm{argmin}}{\mathcal{L}}_m({\mathcal{U}}_0,{\mathcal{U}}_1^t,{\mathcal{X}}^t,{\eta }_0^t,{\eta }_1^t)\\ {\mathcal{U}}_1^{t+1}& =\underset{{\mathcal{U}}_1\ge 0}{\mathrm{argmin}}{\mathcal{L}}_m({\mathcal{U}}_0^{t+1},{\mathcal{U}}_1,{\mathcal{X}}^t,{\eta }_0^t,{\eta }_1^t)\\ {\mathcal{X}}^{t+1}& =\underset{\mathcal{X}}{\mathrm{argmin}}{\mathcal{L}}_m({\mathcal{U}}_0^{t+1},{\mathcal{U}}_1^t,\mathcal{X},{\eta }_0^t,{\eta }_1^t)\\ ({\eta }_0^{(i)}{)}^{t+1}& =({\eta }_0^{(i)}{)}^t-\rho (({\mathcal{U}}_0^{(i)}{)}^{t+1}-{\mathcal{A}}^{(i)}\ast {\mathcal{X}}^{t+1})\\ & i=1,2,\dots ,K\\ {\eta }_1^{t+1}& ={\eta }_1^t-\rho ({\mathcal{U}}_1^{t+1}-{\mathcal{X}}^{t+1}),\end{array}$$

$$\begin{array}{r}{\mathcal{F}}^{-1}\{{\Omega }_m(\sum _{i=1}^K\overline{\mathcal{F}\{{\mathcal{A}}^{(i)}\}}\cdot \mathcal{F}\left\{{\left({\mathcal{U}}_0^{(i)}\right)}^{t+1}-{\left({\eta }_0^{(i)}\right)}^t\right\}\\ +\frac{{\beta }_1}{{\beta }_0}\mathcal{F}\{{\mathcal{U}}_1^{t+1}-{\eta }_1^t\})\},\end{array}$$

$$(\mathcal{U}^{(i)})^{t+1}_{pqr}=\begin{cases}\frac{-\Big{(}\xi^{(i)}_{1}\Big{)}^{t}_{pqr}+\sqrt{\left\{\Big{(}\xi^{(i)}_{1}\Big{)}^{t}_{pqr}\right\}^{2}+4\beta_{0}G_{pq}^{(i)}}}{2\beta_{0}},&\text{if }r=d,\\ \Big{(}\xi^{(i)}_{0}\Big{)}^{t}_{pqr},&\text{otherwise}.\end{cases}$$

Full text

1 Introduction

Image data taken with a specially designed point spread function (PSF) that encodes, via a simple rotation [1, 2, 3, 4], the source distance can be employed to acquire a three dimensional (3D) field of unresolved sources like space debris [5, 6]. By employing a spiral phase mask with the phase winding number changing in regular integer steps from one annular zone to the next and the zone radii having a square root dependence on their index [2], one can create a PSF that rotates uniformly with changing source distance while largely maintaining its shape and size. Specifically, the off-center, shape-preserving PSF is rotated relative to a reference PSF corresponding to the point source being at the plane of best focus, with the rotation angle being approximately proportional to the source misfocus from the plane of best focus.

Hyperspectral and multispectral imaging techniques constitute another important topic in optics and remote sensing problems. Multiple images are collected in contiguous spectral bands for a wide range of applications from space situational awareness [7, 8, 9, 10], agriculture [11, 12] to military detection [13, 14]. One can classify the material composition of an object using detailed spectral information.

Combining the techniques of PSF engineering and multispectral imaging is a promising direction. In [15], a joint deblurring and sparse unmixing is shown to efficiently analyze the image objects using hyperspectral image data. In [7], the authors consider the different blurring effects in different spectral wavelengths and address the space unmixing problem with estimated PSF in each band. Another direction is PSF engineering with multispectral/multicolor imaging. [16] extends 3D localization using double helix PSF in multispectral imaging problem but its considerations are limited to detecting point sources with the spectral intensity mainly distributed in two bands, i.e., red (440 nm) and green (550 nm). This modification can be performed by adding a spectrally-sensitive, phase-modulating element to the imaging path which is positioned in a plane conjugate to the back focal plane of the imaging device. Fluorescent labels encode point sources into red and green classes for biological contextual information. Furthermore, via chromatic dispersion, the multiple colors are encoded in the PSF [17, 18]. Then the spectral information can be detected even in a single 2D snapshot. However, the spectral energy of each targeted object is required to be peaked in a characteristic narrow band, which is not true for space debris, a fact that precludes their classification in this way. For a description of the NASA/DOD space object material spectral signature library we are using see [19].

We discuss here the problem of 3D localization and classification of closely spaced point sources from simulated noisy image data obtained by using Prasad’s rotating-PSF imager [2] in multiple bands. For space debris, this would require an active laser illumination system and collection of light scattered by the debris in multiple narrow bands. The localization problem is discretized on a cubical lattice where the coordinates and values of its nonzero entries represent the 3D locations and fluxes of the sources, respectively. The flux value for one point source in the $i$-th band is the $i$-th entry in its corresponding spectral signature. Finding the locations and spectral signatures of a few point sources on a large lattice is a large-scale sparse 3D inverse problem. Using both Poisson and Gaussian statistical noise models, we describe the results of simulations and propose a three-stage method to extract both the 3D location coordinates and spectral signatures of individual debris particles from rotating-PSF imagery, thus achieving both 3D localization and classification. The Poisson noise case involves an EMCCD sensor operating in the photon-counting (PC) regime [20] while the Gaussian noise case involves conventional CCD sensors operating at high per-mode photon numbers and large read-out noise. We show different implementations for these two noise models and compare our results to those obtained when using only a single band. The joint approach obtains more accurate results than one that first localizes the sources with single-band data and then classifies them with multispectral image data. We use the terms “band”, “channel”, and “spectral wavelength” interchangeably, depending on the context.

The rest of the paper is organized as follows. Section 2 describes the physics model for the single-lobe rotating PSF for both single and multiple bands. In Section 3, we review a non-convex optimization method to solve the point source localization problem for a single band image using the Poisson noise model. In Section 4, we introduce multispectral rotating PSF imaging and propose a three-stage method, which we analyze for the Poisson noise case. We adapt this method For the Gaussian noise case in Section 5. We discuss numerical experiments, including comparisons with the single band model, in Section 6, and present some concluding remarks in Section 7.

2 Physics Models for RPSFs at Single and Multiple Spectral Wavelengths

In this section, we review the single-lobe rotating PSF forward model [2] for single and multiple bands. If the number of point sources is $M$, then the 2D observed image $G$ in a single band has the form,

[TABLE]

where $\mathcal{N}$ is the noise operator and $b$ is the uniform background value. Here $\mathcal{H}_{j}(\mathbf{s}-\mathbf{s}_{j})$ is the rotating PSF for the $j$-th point source of flux $f_{j}$. The depth information $z_{j}$ in the 3D position coordinates $(x_{j},y_{j},z_{j})$ is encoded in $\mathcal{H}_{j}$, and $\mathbf{s}=(x,y)$ is a scaled version of the image-plane position vector, $\mathbf{r}$, namely

[TABLE]

where $z_{I}$ is the distance between the image plane and the lens. $\mathbf{s}_{j}=(x_{j},y_{j})$ is the similarly scaled version of the corresponding 2D image-plane location, $\mathbf{r}_{j}$, of the $j$th source. The imaging wavelength is denoted by $\lambda$ and the radius of the pupil is $R$. Following the Fourier optics model, the incoherent PSF for a clear aperture containing a phase mask with optical phase retardation, $\psi$, is given by

[TABLE]

where $\iota=\sqrt{-1}$ and the depth information is encoded in the defocus parameter

[TABLE]

Here, $l_{0}$ denotes the distance between the lens and the best focus point in the object space. The indicator function for the pupil of radius $R$ is represented as $P(\mathbf{u})$. The normalized pupil-plane position vector $\boldsymbol{u}$ is given by dividing the actual pupil-plane position vector, $\boldsymbol{\rho}$, by the pupil radius, $\mathbf{u}=\frac{\boldsymbol{\rho}}{R}$. For the single-lobe rotating PSF, $\psi(\mathbf{u})$ is chosen to be the spiral phase profile defined as

[TABLE]

in which $L$ is the number of concentric annular zones in the phase mask and $\mathbf{u}$ is in the polar coordinate $(u,\phi_{\mathbf{u}})$. In the following, we directly focus on estimation of the location in the $x$-$y$ plane, $(x_{j},y_{j})$, as well as the defocus parameter $\zeta_{j}$. From (4), we directly obtain the depth as

[TABLE]

Since the $(x,y)$ location of the $j$-th point source in the object space is related to its $(x,y)$ location in the image plane by the simple magnification factor $-z_{j}/z_{I}$, its 3D location in the object space is determined in terms of its image plane coordinates as

[TABLE]

in which $z_{j}$ is calculated from the source defocus parameter, $\zeta_{j}$, using (5).

In the multispectral problem, we have a sequence of 2D observed images $\{G^{(i)},\quad i=1,2,\dots,K\}$ corresponding to $K$ spectra. We need to solve for the 3D locations of the point sources and their flux values for each spectral wavelength, namely

[TABLE]

where $\mathbf{f}_{j}=(f^{(1)}_{j},f^{(2)}_{j},\dots,f^{(K)}_{j})$ and $\zeta_{j}$ is determined from Eq. (4) in which $\lambda$ is set equal to 400 nm. Each point source has a value for its flux in each spectral band, and the vector of those values defines its spectral signature. Thus the forward model is

[TABLE]

where $\mathcal{H}^{(i)}_{j}$ is defined in terms of $\zeta_{j}^{(i)}$. Here the symbols $\mathbf{s}^{(i)}$, $\mathbf{s}^{(i)}_{j}$, and $\zeta_{j}^{(i)}$ generalize the corresponding single-wavelength symbols to multiple wavelengths, namely

[TABLE]

and the spatially uniform background flux $b$ is assumed for simplicity to be the same for each spectral wavelength.

Note that the depth information $z$ is transformed into a defocus parameter called $\zeta$. For a given defocused source, since $\boldsymbol{s}$ and $\zeta$ both scale inversely with the wavelength, its PSF at a longer wavelength is less rotated but spatially more extended than that at a shorter wavelength. Thus, for example, for the RGB (red green blue) wavelenths of $\lambda=660,\,550,$ and 440 nm, if we denote the $\zeta$ parameter of the source for the blue wavelength as $\zeta_{B}$, then the same parameter for the green and red wavelengths will be equal to $\zeta_{G}=\frac{4}{5}\zeta_{B}$ and $\zeta_{R}=\frac{2}{3}\zeta_{B}$, respectively, corresponding to progressively smaller rotations from the blue to the red wavelengths. Correspondingly, the spatial footprints of the PSF blur functions, on the other hand, get increasingly extended, a fact that we incorporate in our simulations, in effect, by changing the physical dimensions of the image plane in proportion to the wavelength.

3 Numerical Localization Schemes for a Single Wavelength

Here, we review the forward model for the single wavelength problem covered in [21, 22]. In order to estimate the 3D locations of the point sources, we assume their distribution is approximated by a discrete lattice $\mathcal{X}\in\mathds{R}^{m\times n\times d}$. The indices of the nonzero entries of $\mathcal{X}$ are the 3-dimensional locations of the point sources and the values at these entries correspond to the fluxes, i.e., the energy emitted by the illuminated point source. The 2D observed image $G\in\mathds{R}^{m\times n}$ can be represented as

[TABLE]

where $b$ is background signal value, $\mathbb{1}$ is a matrix of $1$s of size the same as the size of $G$, and $\mathcal{N}$ is the noise operator. Here $\mathcal{A}\ast\mathcal{X}$ is the convolution of $\mathcal{X}$ with the 3D PSF $\mathcal{A}$, in which the latter is a 3D tensor that comprises a sequence of images with respect to point sources at different depths. The $j$-th slice is the image corresponding to a point source at the origin in the $(x,y)$ plane and at the defocus parameter $\zeta_{j}$ which encodes the depth information. The dictionary $\mathcal{A}$ is constructed by sampling depths at regular intervals in the range, $\zeta_{j}\in[-\pi L,\ \pi L]$ over which the PSF performs one complete rotation about the geometric image center before it largely break apart. Here $\mathcal{T}$ is an operator for extracting the last slice of the cube $\mathcal{A}\ast\mathcal{X}$ since the observed information is a 2D snapshot.

Here we consider $\mathcal{N}$ to be the Poisson noise operator, which is data-dependent. Determining the source locations and fluxes is a large-scale sparse 3D inverse problem. In [22], Kullback-Leibler (KL) divergence [23] is used as the data-fitting term,

[TABLE]

to which is added a nonconvex regularization term to enforce sparsity, and the minimization problem for 3D localization for a single wavelength amounts to

[TABLE]

where $a$ is fixed and determines the degree of non-convexity and $\mu$ is the regularization parameter. The minimization problem (9) is solved by the iteratively reweighted $\ell_{1}$ algorithm [24] with the iterative scheme being

[TABLE]

In the $\mathcal{X}$-subproblem of (10), the alternating direction method of multipliers (ADMM) [25] is used by introducing two auxiliary variables, namely $\mathcal{U}_{0}$ and $\mathcal{U}_{1}$. Then the augmented Lagrangian function $\mathcal{L}_{s}(\mathcal{U}_{0},\mathcal{U}_{1},\mathcal{X},\eta_{0},\eta_{1})$ is

[TABLE]

where $\eta_{0},\eta_{1}\in\mathds{R}^{m\times n\times d}$ are the Lagrange multipliers and $\beta_{0},\beta_{1}>0$ and $\|\mathcal{X}\|$ is equal to the $\ell_{2}$ norm of the vectorized $\mathcal{X}$. Here, $I_{+}(\mathcal{X})$ enforces the nonnegative constraint,

[TABLE]

Therefore, the iterative scheme for the $\mathcal{X}$-subproblem of (10) is

[TABLE]

where the superscript $t$ indexes the inner loop iterations, as opposed to the superscript $k$ used for outer iterations in (10). Here $\rho\in\left(0,\frac{1+\sqrt{5}}{2}\right)$ is the dual steplength. We choose $\rho=1.618$ in our numerical tests [22]. By using [22, Proposition 4.1], we can solve the $\mathcal{U}_{0}$-subproblem (11a) by

[TABLE]

where $\xi_{1}^{t}=1-\beta_{0}b-\beta_{0}\xi_{0}^{t},$ and $\xi_{0}^{t}=\mathcal{A}\ast\mathcal{X}^{t}+\eta_{0}^{t}.$ The $\mathcal{U}_{1}$-subproblem (11b) is solved by soft-thresholding under nonnegative constraint. So the closed-form solution is given by

[TABLE]

The $\mathcal{X}$-subproblem (11c) is a least squares problem. We turn the convolution into componentwise multiplication by using the Fourier transform. Then the closed-form solution for $\mathcal{X}^{t+1}$ reads as

[TABLE]

where $\Omega_{s}=\left(|\mathcal{F}\{\mathcal{A}\}|^{2}+\frac{\beta_{1}}{\beta_{0}}\right)^{-1}.$ Here $\left|\mathcal{X}\right|^{2}$ and $\overline{\mathcal{X}}$ are the componentwise operations of the square of the absolute value of $\mathcal{X}$, and complex conjugate of $\mathcal{X}$, respectively. By assuming the boundary condition to be periodic, we use the 3D fast Fourier transformation (3D FFT) to compute (18) efficiently. For real data, the point sources may be not on a grid, which means the discrete model may not be accurate. In order to avoid missing point sources, the regularization parameter $\mu$ is kept small, which can potentially lead to over-fitting. Our optimization solution generally contains tightly clustered point sources, so we need to regard any such cluster of point sources as a single point source. The same phenomenon has been observed in [26, 21, 27]. We apply a post-processing approach following [27]. The method is based on the well-defined tolerance distance for recognizing tightly clustered neighbors as single point sources. We compute the centroid of each such cluster, which we regard as the location of a single point source.

4 Numerical Scheme for Joint Localization and Classification using Multiple Bands

We now propose a three-stage method for multispectral rPSF imaging. Our aim is to localize in 3D and to classify the materials of unresolved space debris (i.e., point sources) based on multispectral information.

Denote the observed image for the $i$-th spectral band as $G^{(i)}$. We need to restore the 3D tensor $\mathcal{X}^{(i)}$ from the observed image. Note that the positions of the non-zero entries of $\{\mathcal{X}^{(i)}\}$ corresponding to the source locations are the same for all $i=1,\dots,K.$ The values of nonzero entries in each spectral band correspond to the fluxes of the point sources in that band. Since the PSF differs across spectral wavelengths, our dictionaries corresponding to different spectral bands are also different. The discretized forward model becomes

[TABLE]

where $\mathcal{A}^{(i)}$ is the 3D tensor corresponding to the $i$-th spectral band. In this section, we consider the Poisson noise case for imaging space debris [20]. The Gaussian noise case will be introduced in Section 5.

4.1 Stage 1: 3D Localization Using Multiple Bands

In the first stage, we consider 3D localization via all the spectral images. One approach is to restore $\mathcal{X}^{(i)}$ by each $G^{(i)}$, separately, based on the algorithm in Section 3. However, the separate single-wavelength problems are coupled by the fact that the nonzero pattern must be consistent across $\{\mathcal{X}^{(i)}\}$, as it corresponds to the true physical locations of the point sources. Enforcing this constraint while solving all single-wavelength problems would be computationally expensive.

Another reasonable approach is to perform the localization using a single wavelength. Shorter wavelengths lead to smaller side length of the PSF image, so choosing the image corresponding to the smallest available wavelength might be a good choice. However, this approach ignores the observed information from all other spectral bands, which we will see leads to suboptimal results.

Here, we try to utilize all the information of spectral images while maintaining computational efficiency. We make the simplifying assumption that for each point source the intensities in the various spectral bands are all equivalent. That is to say, we do the 3D localization based on assuming a uniform spectral signature for each point source. We correct this assumption in Stage 2, where we more accurately estimate the intensities of the point sources in the different spectral bands. Therefore, in Stage 1, we assume that $\mathcal{X}^{(i)}$ are the same for all $i$, and the discretized forward model becomes

[TABLE]

Our numerical scheme attempts to minimize the sum of the errors across bands, and we adapt the KL-NC model in [22] from single wavelength to multispectral images. The optimization (9) becomes

[TABLE]

where $K$ is the number of bands. Similarly, we use the iteratively reweighted $\ell_{1}$ algorithm [24] to solve (15):

[TABLE]

Here the subproblem of (16) is solved by ADMM. We introduce $(K+1)$ auxiliary variables, namely $\mathcal{U}_{0}^{(1)},\mathcal{U}_{0}^{(2)},\dots,\mathcal{U}_{0}^{(K)}$ and $\mathcal{U}_{1},$ and we have $(K+1)$ Lagrange multipliers. Then the augmented Lagrangian function $\mathcal{L}_{m}(\mathcal{U}_{0},\mathcal{U}_{1},\mathcal{X},\eta_{0},\eta_{1})$ is

[TABLE]

The iterative scheme is

[TABLE]

The $\mathcal{U}_{1}$-subproblem is the same as (11b) in the single band image, so it has the same closed-form solution (12).

The $\mathcal{X}$-subproblem is a least squares problem. Compared to (11c), the first part of objective function is the summation. Similarly, we rewrite the convolution into componentwise multiplication by using the Fourier transform. The closed-form solution of $\mathcal{X}^{t+1}$ reads as

[TABLE]

where $\Omega_{m}=\left(\sum\limits_{i=1}^{K}\left|\mathcal{F}\{\mathcal{A}^{(i)}\}\right|^{2}+\frac{\beta_{1}}{\beta_{0}}\right)^{-1}$. As in the single band case, we use the 3D FFT to compute (18) efficiently.

The solution of the $\mathcal{U}_{0}$-subproblem is the same as (11a) for each $\mathcal{U}^{(i)}_{0}$. Therefore, we need Proposition 1 in [28]. Then the closed-form solution of $\mathcal{U}_{0}^{(i)}$ is

[TABLE]

where \Big{(}\xi^{(i)}_{1}\Big{)}^{t}=1-\beta_{0}b-\beta_{0}\Big{(}\xi^{(i)}_{0}\Big{)}^{t}, and \Big{(}\xi^{(i)}_{0}\Big{)}^{t}=\mathcal{A}^{(i)}\ast\mathcal{X}^{t}+(\eta_{0}^{(i)})^{t}.

Note that $\{\mathcal{U}_{0}^{(i)}\}$ and $\eta_{0}^{(i)}$ are a sequence of 3D tensors. In our implementation, we store them in two 4D tensors and use component-wise operators for 4D tensors in the updating scheme for the $\mathcal{U}_{0}$ and $\eta_{0}$ subproblems. In addition, both $\{\eta_{0}^{(i)}\}$ and $\eta_{1}$ can be stored in the Fourier domain to reduce per-iteration computations.

Remarks: Similar to the single band case, we need to do post-processing to remove some clustered false positives. We compute the centroid for each clustered point source and then represent the clustered set by the centroid [22]. We subsequently perform false positive removal, as described in the following section.

4.2 Stage 2: Alternating approach for estimating spectral signatures and removing false positives

In the 1st stage, we obtain the locations of the point sources

[TABLE]

where $M$ is the number of point sources. In the 2nd stage, our goal is to more accurately estimate the spectral signature of each point source. Here the vectorized PSF image corresponding to the $j$-th point source and the $i$-th band is denoted as $\mathbf{h}_{j}^{(i)}$, which is normalized, i.e., $\mathbf{1}^{T}\mathbf{h}_{j}^{(i)}=1$, where $\mathbf{1}$ is the vector of all 1’s. We let $L$ be the total number of pixels in the vectorized data array, so $L=mn.$ Then we define a system PSF matrix $H^{(i)}$ for the $i$-th band with

[TABLE]

The vectorized observed image in the $i$-th band is denoted by

[TABLE]

where $\mathbf{f}^{(i)}\in\mathds{R}^{M}$ is the flux vector for all estimated point sources in the $i$-th band. In second stage, we follow [22] to estimate the flux information in each band via the iterative algorithm:

[TABLE]

in which $\mathbf{f}^{(i)}_{G}=\left(H^{(i)}\right)^{+}(\mathbf{g}^{(i)}-b\mathbf{1})$, $\left(H^{(i)}\right)^{+}$ being the pseudo inverse of $H^{(i)}$, is the solution corresponding to the Gaussian noise model and

[TABLE]

where $\mathbf{e}_{p}$ is the $p$-th canonical basis unit vector.

The iterative scheme can be implemented as a fixed point iteration away from the starting Gaussian noise solution $\mathbf{v}_{G}$ to the final Poisson noise solution, a process that shows the connection between these two noise models.

In multispectral images, this iterative scheme is processed for each $\mathbf{f}^{(i)}$ separately. Then the estimated spectral signature for the $j$-th point source is

[TABLE]

Unlike the single band image case, we next utilize the estimated spectral information to refine the localization and spectrum estimation. We do these processes in an alternating fashion for further removal of false positives in the estimated point sources. In [22], we observe that the flux value in false positives is usually lower than the value of most true positives. However, it is impractical to set up the threshold with a single band as it requires some information on the flux value in ground truth. In multispectral images, with more information on fluxes in each band, we can estimate the threshold.

Here we propose two criteria for identifying the false positives. The first is the positivity of the fluxes. Note that the iteration scheme (19) does not require a positivity constraint. In the numerical tests, we observe negative flux in some bands for some of the false positives. This can be explained as the effect of false positives compensating for an overestimate of some true positive flux values. The second criterion involves the summation of the spectral values in each point source. We know the flux varies for different materials in each band and we may get a large flux value for false positives when considering only one band. However, because the PSF images are different in different bands, it is very seldom that we have large flux values in every band for false positives. Thus, the probability of having a large flux value in every band for false positives is small. In other words, multispectral imaging data contain more information to detect the false positives. Moreover, the ground-truth spectrum is normalized for each point source. The summation of the spectra of ground truths is in a certain range, while we expect, as we confirmed in our numerical tests, false positives fail to reach this range of energy. In fact, we observed false positives typically to have much lower summation values

Here, we can regard the point source a as false positive if its corresponding estimated spectral signature $\mathbf{f}$ satisfies

[TABLE]

where the threshold

[TABLE]

and $\mathbf{f}^{\ast}_{k}$ denotes the true spectral signature for the $k$-th material. The value $\min\limits_{k\neq l}\frac{\mathbf{1}^{T}\mathbf{f}^{\ast}_{k}}{\mathbf{1}^{T}\mathbf{f}^{\ast}_{l}}$ is computed from the spectral library [36] and the choice of $\gamma$ is robust in not being tightly linked to the ground truth; see also Section 6.3. The algorithm for our 2nd stage is summarized as Algorithm 1.

4.3 Stage 3: Multispectral Classification of Point Sources

In the last stage, we classify the materials of space debris based on the estimated spectra in the 2nd stage. There are many methods for classification or clustering developed for hyperspectral images. See e.g., [29, 30, 31], which determine the underlying materials in each pixel. Here, for the multispectral rotating PSF image, we detect the material for each point source instead of for each pixel in the observed image. The measured spectrum at each point source is assumed to be a linear combination of spectral signatures (called endmembers), so we may formulate the linear spectral mixture model as follows:

[TABLE]

where the estimated spectra is

[TABLE]

Here, $M\in\mathds{R}^{N\times K}$ is a spectral library containing spectral signatures of $N$ endmembers with $K$ spectral bands and $E$ is the computational error from the former two stages. Here we need to solve $X\in\mathds{R}^{M\times N}$ whose entries are the coefficients for the linear combination of different endmembers. Here the summation of entries in each row of $X$ is one. We assume each space debris contains only a few kinds of materials, so $X$ is a sparse and nonnegative matrix. We formulate the classification problem as a least-squares minimization problem. In order to determine the matrix $X$, we add the abundance nonnegativity constraint $X\geq 0$ (i.e., each entry is nonnegative) and the abundance sum-to-one constraint

[TABLE]

Therefore, the minimization for spectral unmixing classification is

[TABLE]

Note that we do not need to involve a total variation (TV) term for requiring the sparsity of the gradient in $X$ as in [7, 32, 33, 15]. This is because the objects we classify are not pixels in the image but point sources of space debris which are independent of one another. The minimization is solved by an Interior Point Least Squares solver [34, 35]. The proposed algorithm for all the stages is summarized as Algorithm 2.

5 Gaussian Noise Case

In this subsection, we briefly discuss cases where the point source images are contaminated by additive Gaussian noise, rather than signal-dependent Poisson noise. The noisy forward model now becomes:

[TABLE]

where $N^{(i)}$ is additive Gaussian noise with zero mean and standard deviation equal to 10% of the highest pixel value in the $i$-th band of the original image $\mathcal{T}(\mathcal{A}^{(i)}\ast\mathcal{X}^{(i)})$.

In Stage 1, we use the $\ell_{1}$ norm as the regularization term. Then the optimization model is

[TABLE]

Because the model is convex, we use ADMM directly without any outer iteration. Upon introducing $(K+1)$ auxiliary variables $\mathcal{U}_{1}$ and $\mathcal{U}_{0}^{(i)}\ i=1,\dots,K$, we obtain the following augmented Lagrangian function $\mathcal{L}_{g}(\mathcal{U}_{0},\mathcal{U}_{1},\mathcal{X})$:

[TABLE]

The iterative scheme is

[TABLE]

In Stage 2, we still use an alternating scheme to estimate the spectral signatures and then refine the 3D localizations. However, the inner loop only minimizes a least squares error, which yields the estimated flux in the $i$-th band, $\mathbf{f}^{(i)}$ as

[TABLE]

The classification algorithm for the Gaussian noise case in Stage 3 is the same as the one for the Poisson noise case, as it is independent of the noise model.

6 Numerical Experiments

Here we test our proposed numerical algorithms for the localization and classification of point sources using multispectral rotating PSFs and compare the results with those obtained using the single band approach. All the numerical experiments are conducted on a standard desktop with Intel i7-6700, 3.4GHz CPU and $\mathrm{MATLAB\ 9.2\ (R2017a)}.$

We use real spectral data provided by Kira Abercombie [36]. This test data corresponds to materials which make up various satellites, rocket bodies, etc., from the NASA JSC Space Object Materials Spectral Database. Some of the data was collected from lab tests, and other data was from flown satellites imaged from ground-based telescopes.

Similar data was also used in [7, 10, 15]. We assume here that the space debris in question consist of parts of man-made orbiting objects with each debris being unresolved and thus treatable as a point source by our imaging system. The signatures cover a band of spectra from 400nm to 2500nm over 100 evenly distributed sampling points. Each material type shows a different spectrum based on its composition. Using low-resolution reflectance spectroscopy and comparing the absorption features and overall shape of the spectra, it is possible to determine the material type of man-made orbiting objects.

We assume that each point source contains only one material and test our ideas on 3D localization and classification. We consider five different kinds of materials. The spectral signatures and the names as well as the corresponding false colors for the materials are shown in Figure 1.

For the Poisson noise case, only the bands with wavelengths between 400 nm and 1000 nm are considered. In the Gaussian noise case, we only consider the higher wavelength bands between 1500nm and 2500nm, over which the photon number per mode is likely to be large compared to 1 and the statistics of the intrinsic photon number more likely to be Gaussian, much like the read noise, which would combine to yield Gaussian statistics for the observed counts at the sensor.

We present results of detailed numerical experiments only for the case of Poisson noise. In Section 6.6, we merely summarize our results for Gaussian noise. Here all 2D simulated observed images are described by 96-by-96 matrices. The number of zones in the spiral phase mask that generates rotating PSF imagery is 7, and the aperture-plane side length is 4 for the Poisson noise case and 2 for the Gaussian noise case. Based on [22], both the 3D dictionary and our discretized 3D space contain 21 slices in the axial direction, with the corresponding values of the defocus parameter, $\zeta$, distributed uniformly over the range, $[-21,21]$. Since the ground truth point sources may not be on grid points, we consider a point source to be recovered if its location is estimated correctly to within one pixel in the axial direction and two pixels in the transverse directions in our comparisons with ground truth [22].

The fidelity of localization is assessed in terms of the recall rate, defined as the number of correctly identified true point sources divided by the number of true point sources, and the precision rate, defined as the number of correctly identified true positive point sources divided by the number of all point sources obtained by the algorithm, which includes the false positives; see [37]. To quantitatively evaluate the performance of classification, we use the following two widely-used metrics: overall accuracy (OA) - the percentage of correctly classified point sources among the identified true positive point sources and kappa coefficient (kappa) - the percentage of correctly classified pixels corrected by the number of agreements that would be expected purely by chance [38].

The original images and observed images for 5 bands are shown in Figure 3. The corresponding center wavelengths in these five bands are 400.0nm, 548.5nm, 697.0nm, 845.5nm, and 993.9nm, respectively. The number of photons emitted by each point source in a particular band is taken to be 2000 times the normalized spectral power of that band. The corresponding sequence of the 15 simulated space debris, with 3 different debris for each of the 5 materials, is shown in Figure 2. Each space debris is represented as a point source with a different spectral trace representing a different material. In Figure 3, we plot the corresponding original images in each band for this 15 point sources case, as well as the noisy observed images. Some of the PSF images overlap and we cannot recognize their rotation. The angle of rotation changes with changing wavelength. With longer wavelengths, the PSF images spread out due to greater diffraction. All these factors make the 3D localization of the point sources very challenging.

6.1 Comparisons with the Single Band Approach

We compare the use of single band images here by replacing the algorithm in Stage 1 by the single band approach for localization that we reviewed in Section 3. That is to say, we estimate the 3D locations of point sources via the method in [22] for a single frame $G^{(1)}$ corresponding to the shortest wavelength and then we do the same processing in Stage 2 and 3. Note that there is no chance to do the classification without multispectral information for these tests. The single band algorithm [22] can thus achieve only localization. So, the single band result in Figure 5 is set up to use a single band in Stage 1 with our proposed stopping criterion (see Section 6.2 for the details) and use multiband information in Stages 2 and 3.

6.1.1 Sample Trial

In Figures 4 and 5 we show an example of the 15 source case. Here we use information in 5 bands for Stage 2 for both cases and 4 bands in Stage 1 for multispectral images. Figure 4 shows the localization results for single and multispectral bands. Not only do we recover more estimated point sources using multispectral images but we also obtain more identified true positives compared to the single band case. There are 14 identified true positives and 10 false positive point sources in the single band image while multispectral images identifies 15 true positives with 11 false positives. The recall rate is 100% and precision rate is 57.69% for the multispectral images while it is 93.33% and 58.33% for single band images. Figure 5 shows both classification and localization results for these two cases. We also see that the single band image has more false positives than our multiple bands case even though we have run the same process Stage 2 to further remove false positives. The recall rate is 93.33% and precision rate is 87.50% for single band images, which means there are 16 estimated point sources including 14 identified true positives. In multispectral images, both recall and precision rates are 100%. As for classification, we misclassify two point sources for single band case, while only one point is misclassified in multispectral case. Here the overall average is 85.71% and kappa is 82.28% for single band images, while it is 93.33% and 91.67% for the multispectral images. To further explore the performance of classification, we plot the estimated spectra for both cases in Figure 6. We observe that their performances are similar except for the fourth material. The estimated spectra in the single band image is further away from the ground truth. Here we misclassify two point sources out of three, while for our multispectral one we only misclassify one. Note that we plot all the estimated spectra corresponding only to identified true positives with respect to location. For example, in the 2nd material the single band image only has two estimated spectra.

6.1.2 Repeated Trials

We also test for localization and classification using 50 random tests for each case and compute the average value of the performance metrics; see Table 1. Here we use the number of bands in Stage 1 and Stage 2 to represent the different conditions we test. A missing number in Stage 2 means that we skip Stage 2 and show the localization result from Stage 1. In Table 1, the multispectral images perform much better than the single band image in both Stages 1 and 2. We also plot the confusion matrix in Table 2. Neither case performs well in the last two classes. In the single band image, 24.67% point sources in class 4 are wrongly estimated as in class 2. This is because the spectral signature is not estimated well. We know the false positives affect the estimation of flux very seriously and the classification results relies on the performance of localization.

6.2 Effects of the Stopping Criterion

In this subsection, we propose a stopping criterion based on our observations and empirically show its efficiency. In the numerical test, we observe the relative change of $\mathcal{U}_{1}$ fails to be a stopping criterion. In the field of inverse problems, some works have focused on the behavior of data-fitting term to adjust the parameter in the optimization [39]. One of the more famous methods is called the discrepancy principle [40, 41, 42] which tunes the regularization parameter until the data-fitting term is close to the noise level. It motivates us to consider the behavior of data-fitting term for our problem.

Figure 7(a) shows the value of the data-fitting term in each iteration. Here we set the maximum number of outer iterations as 2 and maximum number of inner iterations as 400. The curve of the data fitting term is very flat and then significantly drops down for some iterations before becoming flat again. Note that the data-fitting term in the Poisson noise case does not reach the noise level nor do we add any prior information on the noise level. Therefore instead of the data-fitting term, we turn to its relative change between the previous iteration and the current one, i.e., $\frac{|D_{m}(\mathcal{U}_{1}^{t})-D_{m}(\mathcal{U}_{1}^{t-1})|}{|D_{m}(\mathcal{U}_{1}^{t})|},$ where

[TABLE]

The corresponding plot is shown in Figure 7(b). Even though the curve fluctuates a lot, it drops down to a sufficiently low level after some iterations. Based on this observation, we set up the stopping criterion as

[TABLE]

We generate 50 random tests for the case of 15 point sources and run our algorithm with different values of $\epsilon$. In Figure 8, we see that the recall, precision and OA are similar with different values of $\epsilon$, but the computational time is significantly reduced by setting up $\epsilon$ to be not too small. In the following tests, we set $\epsilon=10^{-5}$.

6.3 Effects of Thresholds in Stage 2

In this subsection, we discuss the behaviors of algorithm when changing the threshold in Stage 2, i.e., using different values of $\gamma$. In the space object spectral data [19] that we use for testing purposes in 5 bands (see the blue curves in Figure 6), the minimal summation of spectra over the maximal summation is 0.425, i.e.,

[TABLE]

According to (20), we need to choose $\gamma$ to be smaller than 0.425. This is to avoid the possible over-estimation in spectra. Here we choose different values of $\gamma$ between 0.1 and 0.5 and evaluate the accuracy in both localization and classification. In Table 3, we observe that the performance is similar when $\gamma<0.4$ while the recall rate decrease greatly when $\gamma\geq 0.4$. These illustrate the robustness of algorithm when $\gamma$ satisfies (20). In the numerical tests, we set $\gamma=0.2$ for both Poisson and Gaussian noise cases.

6.4 Effects of the Number of Bands Used

In the first two stages described earlier, we used only a few bands. We now consider a general multi-spectral problem and ask the question — how many bands do we need to use in stages 1 and 2 to produce acceptable localization and classification results?

We still consider 5 different materials and 15 point sources. Each material is assigned to 3 point sources, as before. In the first test, we fix the number of bands in Stage 2 as 5 and try various numbers of bands in Stage 1. The 5 bands are consistent with those in the previous subsection. When the number of bands in Stage 1 is 1, we use Band 1 for Stage 1, similarly, we use both Band 1 and Band 2 when the number of bands in Stage 1 is 2. Here, we randomly generate 50 tests and compute the average of our performance metrics — recall rate, precision rate, overall accuracy (oa), and kappa as well as the computational times in Stage 1, which are shown in Figure 9. We observe that both localization and classification results get better as the number of bands increases. The improvement is especially significant when increasing the number of bands from 1 to 2. The curves gradually become flat, especially after 4 bands. The computational time in Stage 1 increases almost linearly. All these tests use our proposed stopping criterion. The computational time is around 17 seconds, for even a single band image.

In the second test, we fix the number of bands in Stage 1 as 4, and try different numbers of bands in Stage 2 from 2 to 8. The corresponding eight wavelengths are 400.00nm, 484.84nm, 569.69nm 654.54nm, 739.39nm, 824.24nm, 909.09nm, and 993.93nm, respectively. The fixed four bands in Stage 2 use the first four wavelength’s information. The average recall rate, precision rate, overall accuracy (oa), and kappa as well as computational time in Stage 2 using 50 random tests are shown in Figure 10. We observe that using at least 4 bands leads to acceptable results. Using two bands in Stage 2 results in a value of kappa less than 50%. In addition, with more spectral information, the precision rate increases as we remove more false positives, but the recall rate decreases. This means we have also removed some true positives. The computational time in Stage 2 increases linearly but it is much less than in Stage 1, when we compare Figure 9(b) with Figure 10 (b).

6.5 Effects of Regularization

In this subsection, we discuss effect of regularization in the optimization model used for Stage 1. Since we propose Stage 2 to remove false positives, it is worth our attention to consider whether the convex regularization is good enough for classification. Here we test the $\ell_{1}$ regularization term and compare it with the results of our non-convex approach and similarly use, as before, a set of 50 random tests for each case. The minimization model for the $\ell_{1}$ regularization in the multispectral image is

[TABLE]

The problem is the same as the $\mathcal{X}$-subproblem in (16), except that weights $w^{k}_{pqr}$ have all been set equal to a single parameter, $\mu$. Therefore, we use the same algorithm to solve (25). In Figure 11, localization and classification results are shown. Both models have very similar performance. However, the convex model only takes 10.77 seconds while the nonconvex one needs 16.72 seconds, on average. The performance also illustrates the efficiency of the three-stage method. The second stage improves the localization and loosens the requirement of sparsity in the first stage.

6.6 Gaussian Noise Case

In the last subsection, we summarize our results for the Gaussian noise case. Figure 12 shows a sequence of 2D original images as well as the observed images in 5 bands corresponding to the wavelengths, 1530.0nm, 1736.4nm, 1969.7nm, 2203.0nm, and 2436.4nm. Since the wavelengths corresponding to the selected bands are different from the ones in the Poisson noise case, the effect of diffraction is larger for the Gaussian noise case. To keep the image blurring comparable across the two different wavelength ranges corresponding to the two different noise models, we assume the aperture diameter to be twice as large in the Gaussian noise case. Correspondingly, the side lengths are only half as great as the ones in the Poisson noise case. In Figure 12, the diffraction effect is serious even in band 1. Comparing the results to those in Figure 3 for the Poisson noise case, we see that the background noise is much stronger since Gaussian noise is additive and independent of the signal.

We do the same experiments with different numbers of bands in Stages 1 and 2, and plot the performance in Figure 13. The corresponding wavelengths in Figure 13(a) are the five ones as in above discussion. The eight wavelengths for Figure 13 are 1503.0nm, 1630.3nm, 1757.6nm, 1884.8nm, 2012.1nm, 2139.4nm, 2266.7nm, and 2393.9nm, respectively. In Figure 13(a), we observe that the performance metrics are lower than for the Poisson noise case, but the trend is the same — there is significant improvement from single band to multispectral images. In Figure 13(b), we observe that the classification results do not improve much with an increase in the number of bands used only when that number becomes larger than 6. This is a basis for the poor performance in Figure 13(a). In the Gaussian noise case, we need to use more spectral images in Stage 2.

7 Conclusions and Future Work

In this paper we have proposed a three-stage method for employing multispectral RPSF images to localize and classify space debris. With image data from multiple bands, not only the localization results based on the use of a single band can be improved, but also classification can be better addressed for identifying the material components. We have extended the single band’s nonconvex optimization model and adapted it for the multiple band case in the first stage of the proposed three-stage approach. In Stage 2, we have estimated the spectral signatures using each band separately. An alternative scheme was proposed which can further remove the false positives and achieve better estimation for both locations and spectral signatures. Stage 3 was devoted to classification from the estimated spectral signatures. This method can be adapted to different noise models, as we have showed here for Poisson and Gaussian models, and different regularization terms. The numerical results were tested for the case of a pure material for each space debris component. We illustrated the efficiency and stability in the numerical results and illustrated that much better accuracy can be obtained by using multiple bands. In the future, we will study the more involved joint localization-spectral unmixing problem in which each piece of space debris may contain multiple material components. In addition, spectral image based recovery of extended space debris of different shapes, sizes, and material compositions will be of interest as a future direction for our work reported here.

Funding

The research was supported by the US Air Force Office of Scientific Research under grant FA9550-15-1-0286.