New approach to measure the on-orbit point spread function for spaceborne imagers

Lingling Guo; Zepeng Wu; Liguo Zhang; Jianyue Ren

doi:10.1117/1.OE.52.3.033602

4 March 2013 New approach to measure the on-orbit point spread function for spaceborne imagers

Lingling Guo, Zepeng Wu, Liguo Zhang, Jianyue Ren

Author Affiliations +

Optical Engineering, Vol. 52, Issue 3, 033602 (March 2013). https://doi.org/10.1117/1.OE.52.3.033602

Abstract

A new approach to estimate the point spread function (PSF) of remotely sensed images is proposed here based on multiple natural point-like sources, i.e., subimages of the observed image. First the conditional subimages are extracted, and then a blind deconvolution technique is used to derive the PSF from these subimages. For a sampled imaging system with signal-to-noise ratio <20 dB , this method can provide a relatively accurate PSF result. Moreover, the estimated PSF can be applied to image restoration using Wiener filtering or other nonblind restoration methods.

1. Introduction

Spatial characteristics of remote sensing images are generally described by point spread function (PSF) for sampled imaging systems. PSF describes the blur produced in the image by a point in the scene. Unfortunately, a sampled-data system is not shift-invariant on a microscale. Park et al.¹ pointed out that the system PSF had a spurious location-dependent phase dependence, which made the PSF measurement very complicated. Nonetheless, a sampled imaging system can be considered quasilinear over a restricted operating region² and described by a degradation model.

Although the PSFs of spaceborne imagers are measured in laboratory before launch, they may change owing to vibration during launch or aging of materials over time. For this reason, on-orbit PSF measurement is necessary to monitor the actual performance of spaceborne imagers and significant for further image restoration.

The majority of existing PSF estimation methods are developed based on the theoretical method proposed by Smith,³ who aimed at calculating modulation transfer function (MTF, the magnitude of the Fourier transformation of PSF) from its knife-edge response. These methods are frequently used for the reason that knife-edge or step profiles can be easily found in observed images, such as rooflines, farmlands, roadways, and also tarps laid on the ground ahead of time.⁴ In astronomical applications, reasonably bright, isolated stars can serve as point sources to derive the PSF.⁵ Similarly, for earth-observing satellites, PSF can also be obtained by imaging a point source lying on the ground.⁶ In this case, parametric models or multiple point sources are needed to rebuild the PSF. In addition, recently, methods using nonspecific views have been paid more attention. Delvit et al.⁷ proposed a univariate approach to assess MTF using nonspecific scenes.

In this paper, we present a new approach to measure the on-orbit PSF using multiple opportunistic point-like spots. Unlike the point source methods that involve single or multiple artificial spotlights,⁶ this method makes use of cross-relation among subimages so that we do not have to lay artificial targets ahead of time or put effort into collecting information about the real feature in each subimage. First of all, more than one point-like subimage in certain conditions is extracted from the observed remote sensing image. Once the subimages are detected, the PSF can be obtained using blind deconvolution. As a result, we get a two-dimensional (2-D) PSF directly. In the presence of noise, the diversity gain resulted from inherent diversity of natural subimages provides significant improvement in the accuracy of PSF measurement.

The rest of the paper is organized as follows. Section 2 introduces the detailed scheme of our approach. In Sec. 3, we examine the proposed method for its ability to assess PSF and the impact of additive noise on the accuracy. In Sec. 4, we illustrate that this method can be successfully applied to spaceborne imagers for the purpose of image restoration. The last section is for the conclusion and final remarks.

2. Methodology

During the acquisition process, optical imaging systems causes image degradations due to cumulative effects of the instrumental optics (diffraction, aberrations, focusing error) and satellite movement. Typically, the degradation can be modeled by a linear system characterized by its PSF.²^,⁸ The relation between observed image $z$ and the original scene $u$ is described by

Eq. (1)

z = h * u + η,

where

h

is the system PSF and

η

is an additive noise. About the noise

η

, we know only some statistics, such as its mean and variance. As we can see, the observed image

z

is the only known variable in Eq. (1). It is an ill-posed problem to solve and find the solution of

h

. Besides, the current trend is to make spaceborne cameras with higher resolution and wider field of view (FOV) so that the image

z

is considerably large in size. This paper is built on the idea of reducing computational cost by using subimages of

z

. If

P

different subimages are available, then they satisfy the following equation:

Eq. (2)

z_{p} = h * u_{p} + η_{p}, p = 1, 2, \dots P,

where

z_{p}

is the

p

th subimage from the observed image

z

of uniform size

m_{z} \times n_{z}

,

u_{p}

is the original object corresponding to

z_{p}

, and

η_{p}

is independent noise. As pointed out by Holst,² the system can be considered globally shift-invariant so that for each subimage,

h

is considered unchanged. When higher-order aberrations such as comatic aberration cannot be ignored in the system, the assumption of shift-invariance is no longer valid and thus by imposing such an assumption, an averaged PSF is eventually obtained.

2.1.

Blind Deconvolution Method

To solve Eq. (2), i.e., determine the PSF evaluation $\hat{h}$ , we adopt a classic approach of minimizing a regularized energy function. The energy function takes the following form:

Eq. (3)

E (h, u_{1}, u_{2}, \dots, u_{P}) = \frac{1}{2} \sum_{p = 1}^{P} {‖ h * u_{p} - z_{p} ‖}^{2} + λ_{u} \sum_{p = 1}^{P} Q (u_{p}) + λ_{h} Q (h) + μ R (u_{1}, u_{2}, \dots, u_{P}) .

The first term in $E (h, u_{1}, u_{2}, \dots, u_{P})$ measures the fidelity to the data. The remaining three are regularization terms with positive weighting scalars $λ_{h}$ , $λ_{u}$ , and $μ$ . Regularization $Q (u_{p})$ [and $Q (h)$ ] is a smooth term of the total variation form

Eq. (4)

Q (u_{p}) = \int | {\nabla u}_{p} | d x,

where

| \nabla u_{p} |

is the size of subimage gradient.

x = (x, y)

denotes location in the domain of image-intensity functions.

Q (\cdot)

can also have other forms, such as Tichonov regularization and hypersurface minimal function.⁸

Based on the commutativity of convolution operator, we observe that

Eq. (5)

u_{i} * z_{j} = u_{j} * z_{i} for 1 \leq i, j \leq P

in a noiseless case. [Equation (5) holds only in the case that the convolution in Eq. (2) is full. This discrepancy will be dealt with in Sec. 2.2.] Therefore the regularization term

R (u_{1}, u_{2}, \dots, u_{P})

is defined as

Eq. (6)

R (u_{1}, u_{2}, \dots, u_{P}) = \frac{1}{2} \sum_{1 \leq i < j \leq P} {‖ z_{i} * u_{j} - z_{j} * u_{i} ‖}^{2} .

As a functional of several variables, $E (h, u_{1}, u_{2}, \dots, u_{P})$ is not convex everywhere and allows infinite solution. The mean-value ambiguity is removed by imposing a constraint on $h$ . It takes the following form:

Eq. (7)

\int h d x = 1 .

On the other hand, an appropriate setting of the PSF support $m_{h} \times n_{h}$ can alleviate the shift ambiguity. For fixed $h$ or $u_{p}$ ( $p = 1, 2, \dots, P$ ), $E (h, u_{1}, u_{2}, \dots, u_{P})$ is convex. To find a minimizer of the energy function in Eq. (3), we perform alternating minimizations (AM) of $E$ over $h$ and $u_{p}$ . It alternates between the following two steps:⁸

Eq. (8)

{u_{1}^{(n)}, u_{2}^{(n)}, \dots, u_{P}^{(n)}} = argmin E (h^{(n - 1)}, u_{1}, u_{2}, \dots, u_{P})

and

Eq. (9)

h^{(n)} = argmin E (h, u_{1}^{(n)}, u_{2}^{(n)}, \dots, u_{P}^{(n)})

for

n \geq 1

. Generally the AM method used here is a variation of the steepest-descent algorithm. It starts with descending in the subimage space to reach a minimum and then advances in the PSF space in the direction orthogonal to the previous one. The alternating minimization scheme repeats until the result trends to stable. More precisely, we use the preconditioned conjugate gradient method (PCG function in Matlab) to solve the unconstrained minimization problem [Eq. (8)] and fmincon function to solve Eq. (9). Appropriate initial setting of

h

may help the iterative algorithm to converge to the global minimum quickly.

2.2.

Data Preprocessing

If the blind deconvolution method described above is applied, we first need to extract the subimages from the captured image. Because of the cropping at the boundaries of subimages, the convolution here is a partial convolution. However, Eq. (5) holds only in the case that the subimages $z_{p}$ are of full size, i.e., the convolutions in Eq. (1) are full convolutions. The discrepancy between full and partial convolution operators needs to be eliminated ahead of time.

For notational simplicity, we use one-dimensional signals to describe the difference between full and partial convolutions. For the same convolution result $d = {[d_{0}, d_{1}, \dots, d_{M + N - 2}]}^{T}$ and same convolution kernel $b = {[b_{0}, b_{1}, \dots, b_{M - 1}]}^{T}$ , full convolution takes the form

Eq. (10)

[\begin{matrix} d_{0} \\ d_{1} \\ ⋮ \\ d_{M + N - 2} \end{matrix}] = {(a * b)}_{full} = [\begin{matrix} b_{0} \\ b_{1} & b_{0} \\ ⋮ & ⋮ & ⋱ \\ b_{M - 1} & b_{M - 2} & \dots & b_{0} \\ b_{M - 1} & ⋮ & b_{1} \\ ⋱ & ⋮ \\ b_{M - 1} \end{matrix}] [\begin{matrix} a_{0} \\ a_{1} \\ ⋮ \\ a_{N - 1} \end{matrix}],

while partial convolution takes the form

Eq. (11)

[\begin{matrix} d_{0} \\ d_{1} \\ ⋮ \\ d_{M + N - 2} \end{matrix}] = {(a * b)}_{partial} = [\begin{matrix} b_{M - 1} & b_{M - 2} & \dots & b_{0} \\ b_{M - 1} & b_{M - 2} & \dots & b_{0} \\ ⋱ & ⋮ & ⋮ & ⋱ \\ b_{M - 1} & b_{M - 2} & \dots & b_{0} \end{matrix}] [\begin{matrix} \begin{matrix} a_{- M + 1} \\ ⋮ \\ a_{- 1} \\ a_{0} \end{matrix} \\ a_{1} \\ ⋮ \\ \begin{matrix} a_{N - 1} \\ a_{N} \\ ⋮ \\ a_{M + N - 2} \end{matrix} \end{matrix}] .

Obviously, they involve different lengths of the input signal $a$ . Moreover, partial and full convolutions are the same if input signal $a$ satisfies the following condition:

Eq. (12)

{[\begin{matrix} a_{- M + 1} & \dots & a_{- 1} & a_{0} & a_{1} & \dots & a_{N - 1} & a_{N} & \dots & a_{M + N - 2} \end{matrix}]}^{T} = {[\begin{matrix} 0 & \dots & 0 & a_{0} & a_{1} & \dots & a_{N - 1} & 0 & \dots & 0 \end{matrix}]}^{T} .

The limitation is not on ${[a_{0}, a_{1}, \dots, a_{N - 1}]}^{T}$ , but on the elements outside the subscript range of $0 \sim N - 1$ . Similar to Eq. (12), the difference between 2-D partial and convolutions can be eliminated by demanding the pixels outside the border to be zeros. Unfortunately, pixels in natural scenes and captured images are seldom zeros. (For 10-bit remotely sensed data [range: $0 \sim 1023$ ], the minimal value can be about 200.) A propitious thing we learn from experience is that it is much easier to find point-like spots with uniform background in remote sensing images, which have the form

Eq. (13)

u_{p}^{full} = [\begin{matrix} c_{p}^{D} & c_{p}^{V} & c_{p}^{D} \\ c_{p}^{H} & u_{p} & c_{p}^{H} \\ c_{p}^{D} & c_{p}^{V} & c_{p}^{D} \end{matrix}] = [\begin{matrix} 0 & 0 & 0 \\ 0 & u_{p} - c_{p}^{M} & 0 \\ 0 & 0 & 0 \end{matrix}] + [\begin{matrix} c_{p} & \dots & c_{p} \\ ⋮ & ⋱ & ⋮ \\ c_{p} & \dots & c_{p} \end{matrix}],

where

c_{p}^{V}

,

c_{p}^{H}

,

c_{p}^{D}

, and

c_{p}^{M}

are block matrices consisting of same constants (

c_{p}

) of size

(m_{h} - 1) \times n_{z}

,

m_{z} \times (n_{h} - 1)

,

(m_{h} - 1) \times (n_{h} - 1)

and

m_{z} \times n_{z}

, respectively.

c_{p}

is the local uniform background intensity. Considering Eqs. (7) and (13), we have

Eq. (14)

(z_{p} - c_{p}^{M}) = h * (u_{p} - c_{p}^{M})

in a noiseless case. To apply the blind deconvolution method, therefore, we only need to extract subimages with flat borders and then replace

z_{p}

with

(z_{p} - c_{p}^{M})

in Eq. (2). In the presence of noise, the observed intensity of flat areas may be subject to tiny fluctuations. So when searching for subimages with local uniform background, it is reasonable to allow small variation of intensity in the surrounding area. In summary, the blind deconvolution method proceeds with the following preprocessing steps: (1) Search for conditional subimages in observed remote sensing image; (2) Record the original data

{\tilde{z}}_{p} (p = 1, 2, \dots P)

, each of size

m_{z} \times n_{z}

; (3) For each subimage, compute the local background intensity

c_{p} = \min ({\tilde{z}}_{p})

; (4) Update the subimage data

z_{p} = {\tilde{z}}_{p} - c_{p}

, (

p = 1, 2, \dots P

).

3. Simulation

The simulation demonstrates the capability of our method to compute PSF accurately from multiple degraded images. First, six different images were blurred with the same $5 \times 5 pixel$ PSF (see Fig. 1). The blurred images are directly below their corresponding original images in Fig. 1. The boundaries of the input images were extended by padding with a constant value so that the boundary condition mentioned in Sec. 2.2 could be satisfied. Initial parameters were set as follows: $h^{(0)} = delta function$ ; $m_{h} \times n_{h} = 5 \times 5$ ; $λ_{h} = 1000$ ; $λ_{u} = 0.01$ ; $μ = 10$ .

Fig. 1

(a) Six original images, each of size $5 \times 5 pixel$ . (b) A $5 \times 5 pixel$ point spread function (PSF). (c) Six degraded images, each of size $9 \times 9 pixel$ .

To quantitatively measure the accuracy of estimated PSFs, we use mean squared errors (MSE) here since the real PSF is actually known. The estimation results when $P$ different degraded images were involved are listed in Table 1. As we can see, the MSEs are $< 1 %$ when utilizing more than three degraded images. A similar simulated experiment was conducted again, but Gaussian noise was added to the blurred images in Fig. 1(c). The estimation results are listed in Table 2. For heavy noisy data ( $SNR = 10$ , 20 dB), estimation accuracy is significantly hindered by the noise, especially when using just two images from Fig. 1(c). However, practical spaceborne imagers usually provide data with higher SNRs. In addition, this simulation demonstrates the importance of the total number of subimages, i.e., $P$ , to provide better estimation accuracy.

Table 1

MSEs of estimated PSFs in the noise-free case.

	$P = 2$	$P = 3$	$P = 4$	$P = 5$	$P = 6$
MSE (%)	2.2076	1.2887	0.6298	0.5853	0.4722

Table 2

MSEs of estimated PSFs in presence of different levels of noise.

	$P = 2$				$P = 6$
SNR (dB)	10	20	30	40	10	20	30	40
MSE (%)	14.3232	12.5709	5.2611	2.5556	9.2063	6.6830	1.9772	0.6929

4. Application

4.1.

On-Orbit Evaluation of PSF

In this section, the proposed method is tested using data in the panchromatic band of the China YG-8 spaceborne camera. As discussed earlier, the conditional subimages were extracted from the captured image, as shown in Fig. 2. Based on our experiments, we find that the needed subimages are mostly located in lakes, seas, and broad roads.

Fig. 2

(a) A decreased fraction of the original remote sensing image where the central locations of subimages are marked with squares. (b) Eight $9 \times 9 pixel$ subimages. (The subimages are enlarged and enhanced to make them easier to see.)

To achieve more accurate estimation of PSF, we used all eight subimages in Fig. 2(b) to carry out the blind deconvolution method. Initial parameters were set the same as in simulation (Sec. 3) except for PSF size $m_{h} \times n_{h} = 7 \times 7$ . The 2-D PSF result is shown in Fig. 3.

Fig. 3

Normalized 2-D PSF.

In addition, to verify the correction of the estimated results, we compared the results of our method with slanted-edge method.⁹ The knife-edge method can only obtain line spread function (LSF) and 1-D MTF directly. We accordingly computed the integral of 2-D PSF as well as the MTF curves along both directions. The results are shown in Figs. 4 and 5, respectively.

Fig. 4

Comparison of line spread functions (LSFs) of our method and slanted-edge method.

Fig. 5

Comparison of modulation transfer functions (MTFs) of our method and slanted-edge method.

4.2.

Image Restoration

To improve the quality of observed image, the evaluated PSF is applied to restore the degraded image. Because the PSF is a known variable in Eq. (1) now, we can adopt nonblind restoration methods like Wiener filtering. If the estimated PSF is inaccurate, applying it to Wiener filtering will introduce severe ringing effects and cause deterioration of image quality. Two truncated images were restored using the estimated PSF, as shown in Fig. 6. Apparently, the restored images have better visual effects, and no obvious ringing is exhibited.

Fig. 6

Examples of degraded images [(a) and (c)] and restored images [(b) and (d)].

5. Conclusion

Our proposed method makes use of multiple opportunistic point-like sources to measure on-orbit PSF. The main advantage of this new method is the independence of artificial ground targets. The simulation results have shown that it can provide a relatively accurate result in presence of less heavy noise ( $SNR > 20 dB$ ). Moreover, the PSF of a satellite panchromatic band is estimated, and the derived LSF is also properly close to the result of the slanted-edge method. In addition, the image quality has been improved by using Wiener filtering and the estimated PSF. Because the isolated point-like sources in a relatively uniform background can be found in parking lots, roadways, seas, rivers, lakes, desert scenes, farm ground, etc., the proposed method can be successfully applied to spaceborne sensors as long as these objects are available.

References

1.

S. K. ParkR. SchowengerdtM.-A. Kaczynski, “Modulation-transfer-function analysis for sampled image systems,” Appl. Opt., 23 (15), 2572 –2582 (1984). http://dx.doi.org/10.1364/AO.23.002572 APOPAI 0003-6935 Google Scholar

2.

G. C. Holst, “Imaging system fundamentals,” Opt. Eng., 50 (5), 052601 (2011). http://dx.doi.org/10.1117/1.3570681 OPEGAR 0091-3286 Google Scholar

3.

P. L. Smith, “New technique for estimating the MTF of an imaging system from its edge response,” Appl. Opt., 11 (6), 1424 –1425 (1972). http://dx.doi.org/10.1364/AO.11.001424 APOPAI 0003-6935 Google Scholar

4.

P. W. Nugentet al., “Measuring the modulation transfer function of an imaging spectrometer with rooflines of opportunity,” Opt. Eng., 49 (10), 103201 (2010). http://dx.doi.org/10.1117/1.3497051 OPEGAR 0091-3286 Google Scholar

5.

R. Luptonet al., “The SDSS imaging pipelines,” in Astronomical Data Analysis Software and System X, in ASP Conf. Ser., 0101420 (2001). Google Scholar

6.

D. LegerJ. DuffautF. Robinet, “MTF measurement using spotlight,” in Proc. Geosci. Remote Sens. Symp., 2010 –2012 (1994). Google Scholar

7.

J.-M. Delvitet al., “Modulation transfer function estimation from nonspecific images,” Opt. Eng., 43 (6), 1355 –1365 (2004). http://dx.doi.org/10.1117/1.1724838 OPEGAR 0091-3286 Google Scholar

8.

F. ŠroubekP. Milanfar, “Robust multichannel blind deconvolution via fast alternating minimization,” IEEE Trans. Image Process., 21 (4), 1687 –1700 (2012). http://dx.doi.org/10.1109/TIP.2011.2175740 IIPRE4 1057-7149 Google Scholar

9.

J. T. OlsonR. L. EspinolaE. L. Jacobs, “Comparison of tilted slit and tilted edge superresolution modulation transfer function techniques,” Opt. Eng., 46 (1), 016403 (2007). http://dx.doi.org/10.1117/1.2430503 OPEGAR 0091-3286 Google Scholar

Biography

Lingling Guo obtained her bachelor’s degree in electronic engineering from University of Science and Technology of China in 2011. Then she attended Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences, for her MSc degree. She has been involved in the field of quality assessment and image restoration of Earth observation satellites, especially in the measurement of modulation transfer function of satellite imaging systems and the image restoration and enhancement based on the estimated PSF.

Zepeng Wu obtained his bachelor’s degree in electronic engineering & information science from University of Science and Technology of China in 2011. He attended Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences, for his MSc degree after that. He has been working on image processing and analysis since 2009. His current research involves system design for image processing based on FPGA and real-time signal processing algorithm design.

Liguo Zhang is a researcher at Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences, and has been involved in the research of space optical remote sensors since 2000. He is in charge of the thermal control system design and image quality enhancement for optical remote sensing images.

Jianyue Ren is a researcher at Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences. He has been engaged in research of space optical remote sensors since 1996 and has been in charge of many major state scientific research projects. His work won first prize of the 2007 army scientific and technological progress and second prize of the 2008 national science and technology invention. He and his team have done a significant job to the development of Earth observation satellites.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Lingling Guo, Zepeng Wu, Liguo Zhang, and Jianyue Ren "New approach to measure the on-orbit point spread function for spaceborne imagers," Optical Engineering 52(3), 033602 (4 March 2013). https://doi.org/10.1117/1.OE.52.3.033602

Published: 4 March 2013

Access the abstract

JOURNAL ARTICLE
7 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 2 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Point spread functions

Convolution

Imaging systems

Signal to noise ratio

Deconvolution

Modulation transfer functions

Image restoration

1.

Introduction

2.

Methodology

Eq. (1)

Eq. (2)

2.1.

Blind Deconvolution Method

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

2.2.

Data Preprocessing

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

3.

Simulation

Fig. 1

Table 1

Table 2

4.

Application

4.1.

On-Orbit Evaluation of PSF

Fig. 2

Fig. 3

Fig. 4

Fig. 5

4.2.

Image Restoration

Fig. 6

5.

Conclusion

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years