RESTORATION TECHNIQUE FOR PLEIADES-HR PANCHROMATIC IMAGES

Pleiades-HR is a high resolution remote sensing syst em developed by the French Space Agency (CNES) that was launched on the 17 of December 2011 from Kourou Space Centre, French G uyana. Like others high resolution optical satellites, it acquires both panchromatic images, with 70cm spatia l resolution, and lower resolution multispectral images with 2.8m spatial r esolution. Pleiades-HR is an optimized system, which means that the Modulation Transfer Function has a low value a t Nyquist frequency, in order to reduce both the telescope diameter and ali asing effects. Shannon sampling condition is thus m et at first order, which also makes classical ground processing, such as image ma tching or resampling, more justified for a mathemat ical point of view. Raw images are thus blurry which implies a deconvolu tion stage that restores sharpness but also increas es the noise level in the high frequency domain. A denoising step, based upon wave let packet coefficients thresholding/shrinkage tech nique, allows controlling the final noise level. Each of these methods includ es numerous parameters that have to be assessed dur ing the inflight commissioning period: deconvolution filter that dep ends on MTF assessment, instrumental noise model, n oise level target for denoised images, wavelet packet decomposition level . This paper aims to precisely describe the deconvolu tion/denoising algorithms and how their main parame ters have been set up during the inflight commissioning stage. Special attention will be given to structured noise induced by Pleia d s-HR on board wavelet-based compression algorithm


INTRODUCTION
Pleiades-HR is a high resolution optical remote sensing system developed by the French Space Agency (CNES), consisting in a twosome satellite constellation.The first one was successfully launched on the 17 th of December 2011 from Kourou Space Centre, French Guyana.Since it was launched, Pleiades 1A high resolution optical satellite has been thoroughly tested and validated during the commissioning phase led by CNES ( [1]).
Images are simultaneously acquired in both Panchromatic (PA) and multi-spectral (Blue, Green, Red and Near Infra Red), with spatial resolution respectively equal to 70cm and 2.8m for nadir viewing conditions.Panchromatic acquisition is based upon Time Delay Integration (TDI) detectors with a maximum of 20 integration stages in order to comply with Signal to Noise requirements, while multispectral mode relies on a classical CCD detector.Pleiades-HR is an optimized system, which means that the panchromatic Modulation Transfer (MTF) function has a relatively low value at Nyquist frequency, in order to reduce both the telescope diameter and aliasing effects.Shannon sampling condition is thus met at first order, which also makes classical ground processing, such as image matching or resampling, more justified for a mathematical point of view.Raw images are thus blurry which implies a deconvolution stage that restores panchromatic sharpness but also increases the noise level in the high frequency domain.A denoising step, based upon wavelet packet coefficients thresholding technique, allows controlling the final noise level.Following chapters detailed this restoration technique and how restoration parameters have been set up, with a particular focus on compression noise.

MTF
From a signal processing point of view, the acquisition system may be considered as a low-pass filter characterized by its MTF, Fourier transform of the Point Spread Function.MTF is a bidimensional function of (f x , f y ) frequency variables, respectively along the image horizontal and vertical image directions, which tells how high spatial frequencies are weakened by the acquisition chain.In order to minimize aliasing effects during sampling, MTF must have a low value for frequencies lying outside the Fourier domain recovered by Pleiades panchromatic sampling ([-fe x /2;fe x /2]x([-fe y /2;fe y /2], where fe x and fe y stands for sampling frequency along horizontal and vertical , which means relatively low values near Nyquist frequencies.This is also a guarantee for system optimization : MTF high value at Nyquist means that some high frequency details, captured by the instrument, have been lost because of an unsuitable sampling scheme. As depicted in Figure 2, Panchromatic MTF value at Nyquist frequencies is close to 0.16 along both f rows and f columns axes.These estimations are issued from star acquisition thanks to Pleiades pointing capability ( [2]).Stars may be considered as punctual objects, therefore ideally suited for PSF and MTF computation.

Deconvolution
Deconvolution takes place just after ground decompression on each image generated by the 5 TDI panchromatic arrays.Radiometric correction, including both offset and inter detector equalization ([3]), has already been applied on board, in order to minimize compression artefacts and overlapping Inter-Array-Zone (IAZ) have already been rebuilt.The main point is that deconvolution occurs before any resampling, so that deconvolution filter does not depend on the geometric acquisition conditions.Deconvolution filter is defined in the Fourier domain as , where MTF target is the MTF goal after deconvolution.Ideally, the target should be a box function equal to 1 but it would induce high values of D for high frequencies, where MTF is low.The deconvolution filter impulse response would then have an oscillating shape and could induce ringing artefacts near large radiometric transitions.Specular reflections occurences increase with spatial resolution due to highly reflective buildings, car coating or windshields and ringing artefacts may be particularly disturbing, taking into account the 12 bits dynamic range for pixel coding.
Several deconvolution filters were tested and the final choice, depicted in Figure 3, was a compromise between ringing artefacts and final sharpness.

Signal to noise
Noise comes from several sources and may be classified in two main categories: white noise and structured noise White noise spectral density is independent from spatial frequencies before deconvolution.Its standard deviation may be modelled as Assessing the (A,B) parameters set is a major goal of the inflight commissioning since SNR is a key image quality specifications ([3)]).Several assessment techniques have been implemented during the commissioning period, using dedicated acquisitions (slow motion acquisitions, uniform landscapes) or standard acquisitions thanks to an innovative method named BETSI.BETSI takes advantage from the MTF low pass effect, considering that noise dominates signal within some high frequency subset of the Fourier plane.Constant noise variance A may be easily computed from sky acquisition, thanks to Pleiades steering ability.
These techniques finally led to the parameters listed in Figure 5.  On the contrary, structured noise spectral density is non uniform and may be induced by residuals after radiometric corrections or on-board compression algorithm.

Denoising
Denoising technique implemented in Pleiades ground segment deals only with white noise.
Deconvolution transforms white noise into colored noise : since the deconvolution filter magnifies the high spatial frequency content, noise spectral density is higher near Nyquist frequencies.
The objective of the denoising stage is twofold : .identification of areas to be denoised .local transformation of the colored noise into a white noise with specified standard deviation.
Since noise within a deconvolved image depends both on spatial location (x,y) and frequency location (f x ,f y ), denoising technique relies upon a spatial-frequency transform : the wavelet packet decomposition.
Wavelet packet transform may be seen as a filter bank (cf Figure 6) : the initial image is decomposed in several sub images, each of them representing the initial one in a frequency subband with a degraded spatial resolution.Denoising technique consists in reducing wavelet packet amplitude when it is below a threshold that corresponds to what would be observed if there was only noise.Such a threshold may be easily computed by applying the instrumental noise model to the local signal value issued from the lowest frequency subband and multiplying the result by a deconvolution amplication coefficient that depends on the frequency range of the considered subband.In order to make the denoising stage more efficient and to benefit from both spatial and frequency accuracy of each decomposition level, these operations are not realized for one depth D, but rather spred among the whole set of recomposition levels.
Denoising efficiency depends on two main parameters : . The local SNR threshold under which a wavelet coefficient is considered as highly corrupted by instrumental noise and thus to be shrunk .The final noise standard deviation goal.A small target value could induce classical wavelet artefacts such as excessive blur or "butterfly" wavelet artefacts in low radiance areas.This value was tuned to σ noise target =4 digital count, which corresponds to the noise level before deconvolution for signal value S=279 digital count, or S=24Wm -2 Str -1 µm -1 in radiance units.
The denoising stage is time consuming because of the wavelet packet analysis/synthesis complexity, particularly when using spline filters with 27 tap coefficients, in order to be accurate enough in the frequency domain.

Compression algorithm
Onboard compression is mandatory is order to comply with on board storage and data link telemetry constraints.Pleiades onboard compression algorithm is based upon a three level wavelet transform with 9-7 classical Daubechies filters and a bit plane encoding technique [5].It produces an embedded binary bit stream where each block of 16 rows x 1496columns of wavelet coefficients is exactly coded with R*16*1496 bits, where R is the programmed bit rate, obtained through binary stream truncation.This set of wavelet coefficients roughly corresponds to a 16 rows x 1496 columns image subset : Pleiades compression is a bit fixed rate algorithm and the bit rate, expressed in bits per pixel (bpp) units, is thus the main parameter to be tuned during the commissioning period.On board compression applies on the whole set of each spectral band.A discrete set of 6 predefined bit rates may be used : 2bpp, 2.22 bpp, 2.5 bpp, 2.86 bpp, 3.33 bpp and 4 bpp.

Compression noise
Compression noise is somewhat peculiar since it depends both on local image entropy and on the decorrelation transform it is based upon.With a wavelet transform, excessive compression induces blur over low variance areas because nearly all the small amplitude wavelet coefficients are set to 0 by bit stream truncation that leaves apart the least significant bit planes.First coefficients set to zero belongs to the highest frequency subband because of the decrease of coefficient amplitude with frequency, due both to natural landcape behaviour and MTF low pass filtering effect."Butterfy" effects occurs when quite a few coefficients remains : it has the same shape as Daubechies wavelet impulse response.Linear transients may also suffer from aliasing effects.These artefacts are magnified by the deconvolution cannot be dealt with by the Pleiades wavelet packet denoising stage.The only way is thus to find the minimum bit rate that comply with image quality requirement of the final restored image.

Compression bit rate tuning
Due to the embedded bit stream property, images compressed with the higher bit rate (4 bpp) can be decompressed on the ground for every available programming bit rates.Comparing these images to the fully decoded bit stream (4bpp) allows to measure the qualitative and quantitative impact of compression, at any stage of ground processing : decoding, deconvolution, denoising or pan-sharpening ouputs.
Pleiades compression algorithmcommissioning involved a large panel of 24 panchromatic/multispectral acquisitions focusing on large cities acquired with maximum 4bpp rate, with high entropy value, supposed to be worst cases for compression.Comparison criteria included both quantitative and qualitative criteria.
Quantitative measurements were based upon the image difference between the full decoded image and the partially decoded one.Comparison is realized per strips generated by each of the 5 CCDs covering Pleiades swath, both in multispectral and panchromatic modes (cf Fig 7a and 7b).Quantitative measurements include classical global criteria such as Root Mean Square Error (RMSE), Peak to Signal to Noise Ratio (PSNR), but also local criteria measured on small size blocks, such as local signal variance divided by local compression variance error.Local criteria also prove useful for visual investigation.However, qualitative analysis is unavoidable, because structured compression artefacts may not be revealed by any quantitative criteria, neither global nor local.
While quantitative analysis was only done on raw images obtained just after the decoding stage, visual investigation rather involves higher level processed images : restored panchromatic and pan-sharpened multispectral products.
Compression artefacts are more likely to appear in low variance areas included in complex block of lines, which constitutes the bit allocation unit of Pleiades compression algorithm : implicit quantization is suited to the block mean complexity, but not to these areas.
The occurrence probability of such a phenomenon increases with the spatial resolution because of the high percentage of shadowy areas in urban cities.
It should be noticed that completely masking wavelet compression artefacts is tricky : a completely uniform zone, included in a complex block of lines, will be corrupted by instrumental white noise which is likely to contain some tiny "butterfly" artefacts after compression.Using local dynamic adaptation will reveal these artefacts, even if there was no information loss.Finally, quantitative and visual analysis led to slightly increase the bit rate of both multispectral and panchromatic respectively up to 2.86 bpp and 3.33 bpp while preflight foreseen values were 2.5 bpp and 2.86 bpp.

PAN-SHARPENING
Pan-sharpening ground processing allows to generate a multispectral image with the panchromatic spatial resolution.
Such an image cannot be directly acquired for two main reasons .multispectral SNR would be too low because multispectral bandwidths are very narrow compared to panchromatic one and, secondly .the bit budget would explode and would not comply with on board storage and telemetry constraints.Such a technique proves very efficient when panchromatic and multispectral images are natively registered but some problems arise when resampling is needed because multispectral images have a very high MTF level , above 0.3 at Nyquist, which means aliasing artefacts during step 1.Moreover, "soft" panchromatic is basically alias-free since it is a low pass filtered version of the original panchromatic, which is already not so far from meeting Shannon sampling requirement.Therefore, "soft" panchromatic does no longer looks like the zoomed in multispectral aliased images.Consequently, irisation effects may be quite noticeable around high radiometric transitions, as shown in Figure 12a.
A temporary solution was found by tuning some processing chain parameters (restoration and interpolation filters), that minimizes these artefacts without any change in the algorithm.
Moreover, the final solution has already been mocked up Figure 12b clearly points out that this improvment of the current pan-sharpening algorithm completely solves these aliasing problems.This significant evolution-should be deployed within Pleiades ground segments by the end of the year 2012.
Figure 12 pan-sharpened Pleiades images for initial (left) and improved (right) merging algorithms.

CONCLUSIONS
With Pleiades-HR successful in flight commissioning, instrument/sampling optimization ideas, issued from SPOT5 Supermode experience, and applied since the very early stages of Pleiades project, proved to be the right option.
MTF value at Nyquist has to be kept low enough in order to prevent aliasing effects and to minimize the telescope diameter, and thus the system overall cost.With its 0.16 MTF value at Nyquist, Pleiades panchromatic images may be even be considered a bit aliased and its resolution could be improved with a finer sampling grid.
Pleiades incredible agility proved useful to tune restoration parameters, with MTF assessment from star acquisition available just a few days after launch, or instrumental noise model computed thanks to the "slow motion "mode.
Aliasing is definitely an important matter : the multispectral aliasing artefacts have led to significantly modify the initial Pleiades-HR pan-sharpening algorithm.
Finally, Pleiades innovative wavelet compression algorithm gave satisfying results and its embedded bit stream proved useful to adjust the operational bit rate.
Pleiades-1A inflight commissioning period was an exciting period and the whole team is looking forward for next PHR-1B launch.

Figure
Figure 4a : before deconvolution

,
where S(x,y) is the local pixel value for location (x,y) and A, B two constant parameters.A represents the variance of constant noise gathering quantization noise, amplification chain noise and darkness noise.BS refers to the Poisson noise variance.Such a model implies that SNR evolves linearly for low radiances, when Poisson noise is neglectible versus darkness noise, and as ) , ( y x S for high radiance values.

Figure 6
Figure 6 Wavelet packet decomposition A decomposition depth D generates 2 D x2 D subbands.Frequency and spatial accuracy of each subband evolves oppositely with D, according to the Heisenberg principle.The lowest frequency subband may be seen as a 2 D fold zoomed out version of the full resolution image.Each wavelet packet belonging to a given subband may thus be associated to a local signal value in the lowest frequency subband.D=5 was the value adopted for Pleiades wavelet packet analysis.

Figure 7 a
Figure 7 a : Sanaa Pleiades image (5 strips A1 to A5) Fig 8 point out the compression blurring effect on high frequency details for bit rate values less than 2.86 .