Low Bit Rate SAR Image Compression Based on Sparse Representation

Synthetic aperture radar (SAR) is an active remote sensing tool operating in the microwave range of the electromagnetic spectrum. It uses the motion of the radar transmitter to synthesize an antenna aperture much larger than the actual antenna aperture in order to yield high spatial resolution radar images (Curlander & McDonough, 1991). It has been applied to military survey, terrain mapping, and other fields for its characteristics of working in all weather during day and night.


Introduction
Synthetic aperture radar (SAR) is an active remote sensing tool operating in the microwave range of the electromagnetic spectrum.It uses the motion of the radar transmitter to synthesize an antenna aperture much larger than the actual antenna aperture in order to yield high spatial resolution radar images (Curlander & McDonough, 1991).It has been applied to military survey, terrain mapping, and other fields for its characteristics of working in all weather during day and night.
In the last few years, high quality images of the Earth produced by SAR systems, carried on a variety of airborne and space borne platforms, have become increasingly available.With the increasing resolution of SAR images, it is of great interest to find efficient ways to store the high volume of SAR image data at real time or to compress SAR images with higher compression performances for limited bandwidth of communication channel.
There are some special characteristics of SAR images so different from incoherent, optical images, to which compression algorithms are commonly applied, that sensitively affect the design of an image compression algorithm: First of all the speckle phenomena, which results from the coherent radiation and processing and severely degrades the quality of SAR images.When illuminated by the SAR, each target contributes backscatter energy which, along with phase and power changes, is then coherently summed for all scatterers.This summation can be either high or low, depending on constructive or destructive interference.Second the very high dynamic range of SAR images also attributable to the coherent nature of the imaging process.Within a resolution cell of an image, the transduced image domain value is related to the radar cross section per unit area of the corresponding patch of illuminated terrain (Eichel & Ives, 1999) This specific cross section can vary over a considerable range.Most natural terrain, being rough relative to the wavelengths employed, exhibit relatively low values of this parameter, in the vicinity of -15 dBsm/m 2 ; while flat, smooth surfaces such as lake exhibit even lower values.On the other hand, manmade objects, especially of conducting materials with large flat surfaces and right angles, can have specific cross sections of +60 dBsm/m 2 and higher.
These differences mean that encoding/decoding algorithms designed for optical data may not be optimized or even appropriate for SAR data.
The most widely used compression techniques are based on Block Adaptive Quantization (BAQ), due to its simplicity for coding and decoding.The algorithm is based on the consideration that the SAR (complex) signal has commonly a Gaussian, with real and imaginary parts mutually independent and practically uncorrelated between adjacent pixels.It divides the data in blocks, and for each block computes the standard deviation  , in order to determine the optimum quantizer, which adapts to the changing levels of the signal (Kwok & Jhonson, 1989).A non-uniform quantizer (or Lloyd-Max quantizer) that minimizes the Mean Squared Error (MSE) for a given number of quantization levels (Goyal et al., 1998) is commonly used.The minimum block size is selected in order to guarantee a Gaussian statistic within a block; the maximum block size is limited by the fact that the signal power should be approximately constant in the block.
Transform coding algorithms have been also applied on SAR intensity images (Dony & Haykin, 1997, Eichel & Ives, 1999, Baxter, 1999) and on SAR raw data (Pascazio & Schirinzi, 2003).They are based on the decomposition of the signal to be encoded in an orthonormal basis.Then, each decomposition coefficient is approximated by a quantized variable.The role of the signal decomposition is to decorrelate the signal and to make the subsequent quantization process easier.The coding performance depends on the choice of the basis.The best basis compacts the image energy into the fewest coefficients.The small number of significant coefficients in the transformed domain results in a sparse representation, that can be coded with fewer bits, due to its low entropy.An entropy coder can be then used as the last step in the compression scheme.In (Dony & Haykin, 1997) a method which combines Karhunen-Loeve transform and Vector Quantisation is proposed, while in (Eichel & Ives, 1999) simply 2-D Fourier transformation is used, in (Baxter, 1999) a compression system based on the Gabor transform is adopted, in (Pascazio & Schirinzi, 2003) a transform coding compression method using wavelet basis is applied to SAR raw data.Wavelets have also been applied on SAR images by (Zeng & Cumming, 2001, Xingsong et al. 2004), in the first case it has been proposed a treestructured wavelet transform, while in the second case a compression scheme which combines the wavelet packet transform, the quadtree classification and universal trelliscoded quantization has been adopted.
With the aim of reducing the number of significant representation coefficients, and obtaining a sparse representation, overcomplete dictionaries, or frames, have been recently proposed (Goyal et al., 1998).A frame is a set of column vectors, just as a transform, but with a larger number of vectors than the number of elements in each vector.The representation of the observed data in terms of overcomplete basis in not unique, then a constraint have to be enforced to recover uniqueness.To achieve a sparse representation the introduced constraint can be the minimization of the significant coefficients by using an  1 norm based penalty.It can be shown that in this case the problem to be solved is a linear programming problem, that can be viewed as a Maximum a Posteriori (MAP) estimation problem with a Laplacian prior distribution assumption (Hyvarinen et al., 1999).
A possible choice of the basis is the overcomplete Independent Component Analysis (ICA) basis (Hyvarinen et al., 1999, Algra, 2000), allowing to model the data as a mixture of non Gaussian and "almost" statistically independent sources, so that the representation coefficients, due to their scarce correlation, can be efficiently coded using a scalar quantizer.
In this paper the performance of a compression method based on an overcomplete ICA representation, coupled with the use of an entropy constrained scalar quantizer (Pascazio & Schirinzi, 2003), optimized for the Laplacian statistics of the ICA coefficients, and using a proper bit allocation strategy, first proposed in (Budillon et al., 2005) is analyzed in details and proved on different set of real data, obtained with ERS1, COSMO SkyMed and TerraSAR-X sensors.

Overcomplete ICA
Independent Component Analysis (ICA) (Hyvarinen et al., 1999) has been proposed as a statistical generative model that allows to represent the observed data as a linear transformation of variables that are non Gaussian and mutually independent.The model is the following: where is the random vector representing the observed data, is the random vector of the independent components, A is an unknown constant matrix, called the mixing matrix or basis matrix.
The overcomplete ICA paradigm (Hyvarinen et al., 1999) assumes n>m.This means we have a larger number of independent components and we can more easily adapt to signal statistics.Since the matrix A is not invertible, even if it is known, the estimation of the independent components is an undetermined problem that does not admit a unique solution.Then a constraint on the statistical distribution of the ICA coefficients is introduced to solve the problem.
It can be shown (Hyvarinen et al., 1999) that assuming a Laplacian distribution, the optimal estimation of the coefficients ˆi s leads to the minimization of their 1  -norm with the constraint  x As The choice of the Laplacian distribution is convenient with respect to the compression application since permits to have a sparse representation with a small number of non zero coefficient ˆi s .
For the estimation of the basis matrix we adopted a modification of FastICA proposed in (Hyvarinen et al., 1999), that searches for "quasi-orthogonal" basis.

Entropy constrained scalar quantizer and bit allocation
In this section we analyze the performance of a scalar quantizer, optimized for Laplaciandistributed coefficients.In particular, we consider an entropy constrained scalar quantizer, defined as the quantizer minimizing the quadratic distortion for a given value of the entropy of the quantized coefficients (Jayant & Noll, 1984).
It is already known that under the high-resolution quantization hypothesis, i. e. if the number of quantization levels is sufficiently large a uniform quantizer is optimal for a large class of probability distributions [9].However, at low bit rates, a uniform quantizer is not optimal.Among the non-uniform quantizers, we consider the threshold quantizer (or quasiuniform quantizer), which assigns zero to the coefficients whose amplitude lies inside a proper interval [ , ] T T  , and uniformly quantizes the others with a quantization step All the samples whose absolute value exceeds a certain saturation amplitude K are quantized to the highest or lowest quantization level (see Fig. 1) depending on their sign.The saturation factor K can thus be defined as the ratio between the quantizer saturation value and the standard deviation  of the coefficient.Note that this quantizer is symmetric and works with an even number of quantization levels.If the decomposition basis is chosen so that many coefficients are close to zero and few of them have a large amplitude, the threshold quantizer tends to an optimal entropy constrained quantizer (Mallat & Falzon, 1998).The larger the number of coefficients close to zero, the closer the threshold quantizer is to the optimal one.At low bit rates, the decomposition coefficients are coarsely quantized, and many are set to zero.Thus it is convenient to scan the coefficients frame in a predefined order and store the position of zero versus non-zero quantized coefficients in a binary significance map.In the same scanning order, the amplitudes of the non-zero quantized coefficients are also entropy encoded with a Huffman or an arithmetic coding and transmitted together with the coded map.

T -T T+
The total number of bits necessary to encode a data frame is given by the number of bits necessary to encode the significance map, plus the number of bits necessary to encode the significant coefficients.Of course, the total bit rates decreases as the coefficients probability density function (pdf) becomes more peaked, so that the number of significant coefficients decreases and the significance map becomes more correlated (Mallat & Falzon, 1998).
The performance of the quantizer considered depends on three parameters: the threshold value T, the number of the quantization levels L, and the saturation factor K. These parameters must be set in such a way as to minimize distortion for an assigned rate value.The optimization of the threshold quantizer performance with respect to parameters T, L and K can be performed once for a unit variance Laplacian variable.Its performance for Laplacian signals with variance 2  can be simply inferred by that found for an unit variance Laplacian signal, by simply multiplying the threshold value by  , and the obtained distortion by the variance 2  .
A desired bit rate can be achieved for different values of the coder parameters.The optimal coder parameters ˆˆ, , T L K for an assigned bit rate can be found by minimizing the distortion, that for Laplacian distribution can be expressed in an analytical form, following the method presented in (Pascazio & Schirinzi, 2003), for the Gaussian case.
The minimum distortion-rate curve  

ˆˆˆ( , , , ) D R D R T L K
 for a unit variance Laplacian distribution obtained using the threshold entropy constrained quantizer, is shown in Fig. 2 (solid line), where it is compared with the curves obtained with an entropy constrained uniform quantizer (dotted line) (Algra, 2000).Fig. 2. Distortion-rate function for an optimally threshold quantized Laplacian signal of unit variance (solid line), compared with that obtained with an optimal uniform quantizer (dots).
We can note that the optimal threshold quantizer performs better than a uniform one at very low bit rates (lower than 1.8), and reduces to a uniform one for higher bit rates.This behavior can be conveniently exploited when a proper bit allocation is adopted, since in this case it can happen very often that very low rates must be assigned to certain blocks, even if the total bit budget is not very small.If we apply the same threshold quantizer to the entire frame of the coefficients, we obtain the performance shown in Fig. 2. Different results can be obtained by using different quantizer for the different blocks in which the coefficient frame is decomposed.In particular, we can use different quantizers and associated coders optimized for the statistic of each block.A bit allocation algorithm can then be used to distribute bits among the blocks.
A procedure that can be used for optimally distributing the assigned number of bits among the different blocks is described in (Pascazio & Schirinzi, 2003).It imposes an equal average distortion per block and assigns more bits to the blocks with a larger variance respect to those assigned to the blocks with a lower variance.The optimal number of bits for each block is determined by exploiting the minimum distortion curve of Fig. 2 (Pascazio & Schirinzi, 2003).

Numerical results
To test the performance of the proposed method we considered different SAR images obtained with different SAR sensors, using ERS1, COSMO SkyMed and TerraSAR-X data.We wanted to test the performance also on different kinds of areas that due to different backscatter characteristics may have different local statistical characteristics: we used three kinds of images that cover agricultural, suburban, and countryside areas.
Each single-look intensity image has been subdivided in data frames in the azimuth and range directions, respectively.Note that the SAR image pixels are floating point valued with a dynamic range of about 50 dB.Moreover, the SAR images are affected by the presence of speckle, typical of images generated by coherent systems, that has to be preserved to keep the information contained in the image.
Each frame has been subdivided in blocks of 8x8 pixels.The overcomplete ICA basis have been computed using the algorithms presented in [6], and starting from a set of 8x8 training vectors.The set size has to be larger than m.The used value for the ratio m/n is about 0.7.
Then, the ICA coefficients of each frame have been computed using Eq.(2).The ICA coefficients have, then, been quantized using the optimal threshold quantizer of Fig. 2, that, besides exhibiting a better performance, has the advantage of allowing any fractional bit rate value.For the bit budget distribution among the different coefficients vectors, the bit allocation procedure referred in Section III has been adopted.To sum up we followed for each image the scheme reported in Fig. 3 Different quality parameters can be considered to evaluate the performance of the compression method.In particular, one of the most meaningful parameters is the signal-tonoise ratio (SNR), or its reciprocal, the normalized average distortion D computed on the SAR images obtained after the processing of the compressed data.We have chosen to evaluate the normalized distortion in order to compare the obtained results with the distortion-rate curve reported in Fig. 2. We reported also the distortion evaluated on the ICA coefficients.First a single look ERS-1 intensity image, relative to Fleevoland, in The Netherlands, was considered.The original image is shown in Fig. 4.
The (equalized) intensity of a frame is shown in Fig. 5.The obtained statistical distribution of each ICA coefficient, as expected, has a Laplacian behaviour as shown in Fig. 6.The ICA basis are reported in Fig. 7.
The average distortion obtained for different bit rates are presented in Table 1.Note that the bit rate values represent the average number of bits per pixel of the image frame.Being the dimension of the ICA basis larger than that of the observation domain, the average number of bits per coefficient is smaller (it is scaled by the factor m/n).
Quantized coefficients have then been used to reconstruct the corresponding image using Eq. ( 1).The image frame obtained with an average rate per sample R=2, is shown in Fig. 8.
Secondly we considered a single look COSMO SkyMed of Naples surroundings, in Italy, shown in Fig. 9.
The (equalized) intensity of a frame is shown in Fig. 10.The obtained statistical distribution of each coefficient, as expected, has a Laplacian behaviour as shown in Fig. 11.The ICA basis are reported in Fig. 12.
The average distortion obtained for different bit rates are presented in Table 2.
The image frame obtained with an average rate per sample R=2, is shown in Fig. 13.
Thirdly we considered a single look TerraSAR-X of Frankfurt surrounding, in Germany, shown in Fig. 14.
The (equalized) intensity of a frame is shown in Fig. 15.The obtained statistical distribution of each coefficient, as expected, has a Laplacian behaviour as shown in Fig. 16.The ICA basis are reported in Fig. 17.
The average distortion obtained for different bit rates are presented in Table 3.
The image frame obtained with an average rate per sample R=2, is shown in Fig. 18.
It can be noted that in all cases the average image distortions reported in Table 1,2 and 3 are below the value of -11 dB obtained for rate 2 in the curve reported in Fig. 2 for an optimally threshold quantized Laplacian signal of unit variance.
Moreover in all the cases there is no visual appreciable degradation in the reconstructed images using the quantized coefficients of the overcomplete ICA basis.-11.4 -16.4 Table 3. Rate-Distortion values for the TerraSAR-X image.

Conclusion
In this paper the performance of a compression method based on an overcomplete ICA representation, coupled with the use of an entropy constrained scalar quantizer, optimized for the Laplacian statistics of the ICA coefficients, and using a proper bit allocation strategy, has been analyzed in details and proved on different set of real data, obtained with ERS1, COSMO SkyMed and TerraSAR-X sensors.The best performances in terms of ratedistortions are obtained on the TerraSAR-X data frame, since it is relative to an enough uniform reflectivity area, while the worse on the COSMO SkyMed image frame, where more details were present.In any cases the image average distortions are below the one obtained with an optimally threshold quantized Laplacian signal of unit variance.The ICA coefficients exhibits the statistical behaviour forced by the fast ICA algorithm, in fact the probability empirical distributions are in all cases a good approximation of a Laplacian distribution.This behaviour allows to use a threshold quantizer optimized for this particular statistical distribution, allowing to discard many coefficients and to use a high bit rate only for few significant coefficients in order to keep low bit rates.

Table 1 .
Rate-Distortion values for the SAR ERS-1 image.