Wishart-based adaptive temporal filtering of polarimetric SAR imagery

: Temporal ﬁltering for speckle reduction of polarimetric SARimages is described. The method is based on a sequential complex Wishart-based change detection algorithm which is applied to polarized SAR imagery, including the dual-polarization intensity data archived on the Google Earth Engine (GEE). Software for convenient application and analysis is presented. Results compare favorably with, and improve upon, standard spatial and temporal ﬁlters for speckle reduction.


Introduction
A characteristic of SAR images, when compared to their visual/infrared counterparts, is the disconcerting "speckle" effect which sometimes makes visual interpretation very difficult. Speckle gives the appearance of random noise, but it is actually a deterministic consequence of the coherent nature of the radar signal. Although disadvantageous in some contexts, coherence effects in the received radar signal constitute the basis of SAR interferometry and its many applications. However, one is often concerned primarily with SAR intensity (as opposed to amplitude and phase) observations, for example, for target detection, image classification, or forest monitoring, where speckle is indeed a nuisance.
Many adaptive spatial filtering methods have been developed in the past, some of the most popular of which are described in [1][2][3][4][5][6][7][8]. In addition, temporal [9] and hybrid temporal [10] methods have been proposed for time series of SAR imagery. In this paper, we describe a straightforward adaptive temporal filter, which is essentially a by-product of a sequential change detection algorithm for times series of polarimetric SAR images [11][12][13]. In order to determine the so-called Löwner order (direction) of significant changes in the algorithm, it is necessary to keep a running average image whose pixels are updated whenever a change occurs. At the conclusion of this process, the running average is the temporally filtered version of the last image in the series. Pixels which have not changed, or which change only early in the sequence, are well averaged and exhibit little speckle. Those which change later will be less well despeckled. Since there is no convolution involved with a spatial kernel, there is no degradation of spatial resolution whatsoever. We refer to this procedure in the following as an adaptive temporal spatial filter (ATSF) and the despeckled result as the ATSF image. The method is especially useful if one is also performing a change detection study over many time points, in which case, the despeckling is a (free) by-product. The filter is unique, resolution preserving, and may be of interest for purposes other than simple visualization (speckle reduction), such as preprocessing for feature detection or for classification and segmentation by means of machine learning methods, which are insensitive to the equivalent number of looks (ENL), e.g., convolutional neural networks.

Theory
Complex Wishart-based change detection algorithms for bitemporal and multitemporal SAR images have been described in a series of publications [14][15][16][17][18], and applications to long time series of Sentinel-1 imagery archived on the Google Earth Engine (GEE) [19] are given in [13]. To briefly summarize the methodology, given a sequence of k SAR images in polarimetric matrix form, all characterized by the same ENL, the statistical null hypothesis Σ 1 = Σ 2 = · · · = Σ k , where Σ j is the covariance matrix characterizing the complex Wishart probability distribution of an observed pixel in the jth image, is tested against all alternatives. This is referred to as the Q test in the cited papers. If this test is accepted, we are done: for this pixel there is no change over the time span t 1 to t k . If the test is rejected, the hypotheses are tested (called the R j tests [15]). In the event of rejection of the null hypothesis at any time point (significant change), the tests are restarted from that point for the image pixel concerned. See [15] for details. Subsequently, a direction measure for the observed significant changes is determined as follows: The covariance matrix prior to change is averaged over all preceding no-change images in the sequence and subtracted from the covariance matrix immediately after significant change. The Löwner order, i.e., the positive-, negative-, or in-definiteness of the difference matrix is then taken to be a measure of direction [11]. The running average image needed for the Löwner order constitutes an adaptive temporal filter on the last image in the sequence. Pixels which do not exhibit change over the full sequence are maximally averaged, and changed pixels are averaged over the period following the most recent change. In order to reduce the speckle for regions of very recent change, all averages for which the number of contributing pixels is below a threshold may be replaced by the output of an adaptive spatial filter applied to the original image. In our work, we chose the refined Lee spatial filter with 7 × 7 directional kernel [2,20], but any appropriate spatial filter can be used. It should be remarked that the only other free parameter involved in determining the ATSF is the (common) significance level chosen for the Q and R j tests. Of course, as the algorithm involves statistical testing, false negatives may result in changes being overlooked, depending on the significance level chosen.

Materials and Methods
Software for the sequential SAR change detection algorithm and for the adaptive temporal filter described in the preceding section is provided in Matlab (https://people.compute.dtu.dk/~{}alan/ software.html) and in a Python implementation (https://github.com/mortcanty/EESARDocker) written specifically against the GEE API. The latter runs in a Docker container serving Jupyter notebooks to the user's local browser, also see [13,20]. The precompiled Docker container can be pulled and run directly from Docker Hub and is platform independent. The filtered image is exported to the user's Google Drive. Postprocessing Python routines for thresholding/substitution with the Lee or Gamma MAP filters are included in the container. Alternatively, refined Lee filter thresholding/substitution can be run directly within the GEE application. Long time series of archived Sentinel-1 images suitable for adaptive filtering are easily generated with the GEE, but time series from other platforms, such as RADARSAT-2, TerraSAR-X, etc., can be uploaded to the GEE platform and processed in the same way.

Results and Discussion
A Sentinel-1 image of an open-cast coal mine in Germany from 27 September 2018 is shown in Figure 1, together with the ATSF. A time series of 29 acquisitions terminating with the 27 September image was analyzed with the sequential omnibus algorithm, generating the ATSF as a side effect. As can be seen, the speckle is virtually absent from the filtered image, and land-cover boundaries and features are more easily discerned. In Figure 2, the ATSF is compared with refined Lee [2], Gamma MAP [3], and Frost [4] spatial adaptive filters, as well as with the empirical mode decomposition (EMD) temporal filter [9] and an average over the full time series. The time average is essentially the Quegan filter [21] for uncorrelated images with the same numbers of looks and equal means. The three spatially filtered images were produced with the Sentinel 1 Toolbox using default parameters. The EMD filtered image was calculated using the full 29-image time series in decibel (log) scale augmented with five additional images following the target date. The additional images were included in order to avoid edge effects in the empirical mode decomposition, as suggested in [9]. The Python bindings to the C-package libeemd [22] were used to program the method.
Only the Quegan filter is comparable in visual quality to the ATSF image, but being a simple time average, the earth moving machines at their locations on 27 September are no longer visible. This is also the case for the EMD filter, which otherwise has a good visual appearance. In this sense, the ATSF is, like the spatial filters but unlike the Quegan and EMD temporal filters, a despeckled version of the September 27 target image. Figure 3 shows a subset of the reference image from 27 September with the Jülich Research Center and environs. The center itself and the forest surrounding it exhibit very little speckle after filtering (top right). However, speckle is still present in some of the nearby agricultural fields due to recent changes (harvesting). These areas also appear darker because of the lower radar cross-section in the harvested fields. The bottom left figure shows, per-pixel, the number of images used in forming the ATSF as a gray scale image, from black = 1 through to white = 29. The thresholded ATSF image (see Section 2) is shown bottom right for a threshold of 7. The speckle in the recently changed, harvested agricultural fields is now reduced by the Lee adaptive filter.  Figure 4 shows on the right, the ATSF result for a time series of 12 quad-pol RADARSAT-2 images from 2009 and 2010 over an agricultural/mixed forest area near Bonn, Germany. The ENL of the images in the series was 12, so that speckle is in any case reduced in the original image (left). However, the ATSF still improves the visual impression considerably.
Another example of the application of the ATSF filter is given in Figure 5. With its large scale structures (roads, runways, etc.), the image sequence will be used below to illustrate edge preservation/enhancement.  Evaluating and comparing speckle filters is particularly difficult because, apart from in simulations, there is no reference other than the original unfiltered image itself. Table 1 compares the mean bias criterion (MB) as a fidelity measure comparing various filtered images with the original image as suggested in [10], namely where µ 1 and µ 2 are the mean intensities of the original and filtered images, respectively. Four adaptive spatial filters were applied using the ESA Sentinel-1 Toolbox with default parameters along with the EMD and Quegan filters. Larger values indicate better agreement with the original image. The Gamma MAP and Frost filters perform best with this metric on the chosen datasets, but the ATSF achieves acceptable results compared to the spatial filters (the EMD filter was not applied to the RadarSat-2 imagery due to the relatively short length of the time series). The MB metric is marginal, that is, bands are compared separately. On the other hand, the average difference in spectral angle (ADSA) (see, e.g., [23]) allows for the comparison of the complete polarimetric images: where x 1 and x 2 are observation vectors for the original and filtered images. Typical values are shown in Table 2. Here, the adaptive spatial filters slightly outperform the ATSF and EMD filters. Table 2. Comparison of mean spectral angle differences in degrees, best results in boldface. Not surprisingly, the ATSF outperforms the adaptive spatial filters at edge preservation and enhancement, since there is no spatial averaging taking place whatsoever. This can be seen clearly in Figure 6 where the Canny edge detector is applied to several spatial filters, to the EMD filter, and to the ATSF of the Frankfurt airport image in Figure 5. The EMD filter also performs better than the spatial filters but not as well as the ATSF.

Scene (Platform) Lee [2] Gamma [3] Frost [4] IDAN [5] EMD [9] Quegan [21] ATSF
Another means to compare speckle filters is on the basis of the ENL characterizing an image. The ENL describes the degree of averaging applied to the SAR measurements during data formation and postprocessing. The ENL of a speckle filter result is essentially the number of independent looks which would be needed to obtain the same degree of smoothing as that achieved by the filter itself, so that the range of ENL values within an image is an indirect measure of speckle reduction. Figure 7 shows histograms of the ENL values. These were estimated with a multivariate method described in [24]. A robust ENL estimate is extracted from the distribution of small sample estimates collected over the whole scene. Histograms of these small sample estimates collected within a 7 × 7 window passed over the scene are shown in the figure, in this case, for a heavily forested area in Germany exhibiting well-developed speckle. The upper-left histogram is for the target image without filtering and has a mode at about five equivalent looks, corresponding to the 5 × 1 (range × azimuth) multilooking carried out in the original processing. The upper-right hand histogram is for the ATSF using a time series of 29 images ending with the target image. It approximately compares to the refined Lee spatial filter with a 5 × 5 kernel (lower left). Finally, the ratio of the original image intensities to those of the filtered image should ideally, because of the multiplicative noise assumption, be featureless noise [9]. This is largely the case for the refined Lee and Gamma MAP spatial filters (Figure 8). The ratios for the EMD and ATSF temporal filters show some land-use features. Note that the earth moving machines are absent in the ATSF ratio image, as they should be. A summarized comparison of the various speckle filters examined is given in Table 3.

Conclusions
The adaptive temporal speckle filter (ATSF) described here is a side effect of time series analysis with the sequential omnibus algorithm and is seen to be comparable to a conventional spatial filter with respect to the amount of smoothing achieved and fidelity to the original image, while at the same time exhibiting a much better performance at edge preservation and enhancement. It presupposes a sufficiently long and contiguous time series (in the order of 10 or more images) preceding the acquisition date of the target image, as well as a predominance of unchanged pixels over the series. The Python software described in Section 3 exploits the huge Sentinel-1 archive on the GEE to make the filter application simple and fast. The effect of recent changes in the time series can be mitigated by using a spatial filter for the corresponding pixels. In the worst case of predominating recent changes, the ATSF (in our software implementation) reduces to the refined Lee spatial filter.