Change detection in SAR time-series based on the coefficient of variation
CChange detection in SAR time-series based on the coefficient ofvariation
Elise Colin-Koeniguer , ‡ , Fabrice Janez , ‡ and Jean-Marie Nicolas Onera, Universit´e Paris Saclay; [email protected] Telecom ParisTech
Abstract
This paper discusses of change detection in SAR time-series. Several criteria based on thetemporal coefficient of variation are proposed, dedicated to different kind of changes. Firstly,the coefficient of variation is used alone to detect any changes; some statistical propertieshighlights its pertinence. Secondly, two other criteria based on ratios of coefficients ofvariations are proposed to detect long events such as construction test sites or point-eventsuch as vehicles. These detection methods are evaluated both on theoretical statisticalsimulation and on real data for different types of scenes and sensors (Sentinel-1, TerraSAR-X, UAVSAR). In particular, a quantitative evaluation is performed with a comparison ofour solutions with a reference one.
Keywords:
Multitemporal; change detection; time series; SAR; coefficient of variation
1. Introduction
Since the launch of the Copernicus Programme, data becomes available on a full, openand free-of-charge basis; thus many new services can be developed whether for environ-mental, civil, industrial, defense or surveillance applications. Easy access to time-series isparticularly interesting for change detection scenarios: • Using a higher number of images enables to improve the robustness of change detection • Temporal information can help discriminate different classes of change and thus con-tribute to the characterization of changes. Information about a change can also beenriched with its duration. This information can be particularly interesting to controlthe activity rate over a given area.Although the use of time series is particularly interesting in this context, relatively littlework can be found in the literature because this context of access to data is relatively recent,especially for SAR images. However, this data is particularly valuable in this context:unlike optical data, the SAR image is almost insensitive to cloud cover and independentof illumination conditions. Change detection performance is known to be very interesting.Also, this article focuses on the detection of changes in the SAR time series.For this purpose, we can consider two levels of products:
Preprint submitted to Elsevier April 26, 2019 a r X i v : . [ phy s i c s . d a t a - a n ] A p r The first type of product is a visualization product: it can be considered as a soft levelof decision. • The second type of product is an automatic detection binary mask.Classically for visualization, RGB products can represent changes occurring between twoor three dates in a colored composition. When the time series becomes more significant,it is then possible to use these colored visualizations as frames of a video. This idea isproposed in [5] where each frame contains the RGB composition for three successive dates.An alternative approach defines each frame as an RGB composition for a triplet of imagesconstituted by the first reference image and two successive images. This strategy makes itpossible to distinguish between one-point changes or more extended changes.Very few works propose to synthesize the temporal information in a single image prod-uct. The work in [1] proposes a representation approach including interferometric coherenceinformation. However, although this composition relies on the entire time series, it is morefor a classification product than a change detection product.Recently, we proposed a method called REACTIV (Rapid and EAsy Change detectionon Time-series using the coefficient of Variation). To our knowledge, this is the only visu-alization method in the literature dedicated to change. The approach is based on the HSVcolor space, where the critical parameter is the temporal coefficient of variation, which isused to encode the color saturation. Indeed, this parameter has been proven to be particu-larly sensitive to changes. By this way, the brighter the color, the more likely the change is.The assigned color or Hue canal corresponds to the most probable date of the change.If this visualization has been proved to be useful to visualize changes in a given areaquickly, it is natural to look at the automatic detection methods. [8] gives a clear overviewof the multiple change-points detection in multivariate time series. Methods are separatedinto two classes: sequential or online methods, and retrospective or offline analysis. The firstones are dedicated to the processing of data as they are acquired [11], to detect a changepoint as soon as possible after it occurs. The second ones consider the complete data set atonce and look back in time to recognize where a change occurred.In this paper, we are interested in the offline case: we want to detect in a SAR timeseries all changes that occurred in a given area and a given period. A first strategy consistsin considering binary-temporal change detection tests between all possible pairs of datesexhaustively. Resulting tests are represented by a change detection matrix (CDM) [10] Thisapproach is intuitively compelling because it considers the whole time series and the possi-bility of change characterization [19]. In the same vein, multiplying inter-image comparisonsalso makes it possible to provide greater robustness by examining the consistency there isbetween the different tests [3]. However, it is a very combinatorial problem.Therefore, another alternative is to consider the detection of different ruptures. Theomnibus method proposes to test the hypothesis that all polarimetric signals belong to asingle statistical population. This global test can be factored into a binary set of tests,which test iteratively that the element t has the same realization polarimetric matrix thanthe first element of the series. As soon as the test is invalidated, then the test is rejected.It is then possible to consider the following sub-series to search again for another possible2upture. In practice, that means that a subset of binary tests is conducted among all possiblepairs corresponding to the first line in the change detection matrix. This approach has beenapplied on data from different sensors by using a statistical test on polarimetric Wishartdistribution in [16], [4], [12], [17].Finally, the last category of approaches focuses on how to define a criterion to judge thehomogeneity of a whole statistical population. It the hypothesis is rejected, then it leadsto detect a potential statistical rupture without trying to date it precisely. In this case, anentire group is processed without considering bi-date tests.The MIMOSA method studies the difference in statistical properties of geometric orarithmetic means in the presence of change. Always in this category, the coefficient ofvariation is another potentially advantageous candidate because of its simplicity, and its re-markable statistical properties to detect a change, already illustrated as part of the efficiencyof REACTIV visualization method.Note that to our knowledge, there is no published work on the implementation of a cri-terion about the homogeneity of the temporal profile by including the phase information.Coherent change detection his is indeed an even more ambitious challenge: the interfero-metric phase can be affected by changes in the position of objects between two acquisitionsat the wavelength scale. For this reason, in this article, even if we plan to work with imagesin interferometric conditions, we will restrict ourselves to parameters addressing only theamplitudes.The coefficient of variation is, therefore, an ideal criterion for detecting a change whateverits category. The statistical properties of this criterion are at the heart of the strategies forsetting the detection functions. They are therefore studied and presented in section 2.In section 3, we propose, based on these properties, other criteria more specific to distinc-tive categories of changes; including profiles exhibiting a target appearing only on a singledate, typically, a vehicle; and profiles with high returns of longer duration, typically corre-sponding to the activity of a building site. Section 4 evaluates the three previous detectionscenarios. After the presentation of the evaluation protocol, we give the performances of themethods.
2. Theoretical considerations for the coefficient of variation
The coefficient of variation (CV), also known as relative standard deviation, is mathemat-ically defined in probability theory and statistics by σ/µ , where σ is the standard deviationof the signal and µ the mean value. Therefore, it is considered as a normalized measurementof the dispersion of a probability distribution.Because the availability of large time series of SAR images is quite recent, studies ontemporal statistical SAR distributions are still marginal. However, these studies are essentialfor setting up detection methods based on this coefficient of variation.In the specific area of SAR images, the coefficient of variation has already been consid-ered for DInSAR (Differential Interferometry SAR) applications. This technique uses theinformation from selected points having high backscattering and stable through time, that3 igure 1: Three generic case in temporal profiles with increasing coefficient of variation: the deterministicstable class, the natural non deterministic speckle class, and changes we call Permanent Scatterer. In this framework, the coefficient of variation is one candidatefor the detection of these particularly stable scatterers.In the framework of REACTIV, first theoretical studies have shown that the coefficientof variation is also interesting to detect changes in speckle areas.It can, therefore, be considered that the coefficient of variation has different statisticalproperties for at least three categories of temporal profiles: that of a permanent scatterer;that of a natural area of stable speckle, not necessarily correlated in time, but which isstationary; and finally, that of a non-stationary area that we will generally be interpretedas a change. These three general cases are represented in Fig. 1.The first two of these categories concern cases without any changes. On both these cases,it is possible to derive useful mathematical properties about the coefficient of variations.Also, we discuss these properties in the following paragraphs. In order to identify stationary targets, at least two strategies can be chosen. The firstone is to use the interferometric coherence level. The main difficulty relies on the spatialestimation process that implies a spatial resolution loss. The second one relies on theanalysis of the amplitude values A k along the time axis by the so-called dispersion index,that is precisely the coefficient of variation.The amplitude values follow the Rice distribution: it is the superposition of the backscat-tering of a deterministic target with constant values of phase and amplitude, with an additivefully developed speckle. Let µ c be the backscattering parameter of the target, then the Ricelaw can be expressed as: 4 C [ µ c , µ ]( x ) = 2 xµ e ( x µc µ ) I ( 2 µ c xµ )where I is modified Bessel function of the first kind of order 0.Thus, the shape of the Rice distribution depends on the ratio between the amplitude ofthe deterministic component and the amplitude of the speckle component.The coefficient of variation, defined as the ratio of the standard deviation and the average,can be expressed in terms of the statistical moments of any distribution following: γ = σµ = (cid:112) m − m m (1)It is possible to literally express m , m and m in terms of the law parameters. One ofthe most elegant methods is to go through the characteristic function of the second-orderexpressed through Mellin transformations. It makes it possible to obtain the analytic ex-pression of the first moments of the law, by involving the confluent hypergeometric function,also called the Kummer function. All mathematical demonstrations are described in [14].The expressions of the moments thus obtained are not trivial. However, they enable toobtain an analytical expression of the coefficient of variation, either in terms of Kummerfunctions, or in terms of Bessel functions: γ = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) e µ cµ F (2; 1; µ c µ ) π ( F ( ; 1; µ c µ )) − γ = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) e µ cµ (1 + µ c µ ) π ((1 + µ c µ ) I ( µ c µ ) + µ c µ I ( µ c µ )) − λ = µ c µ ,i.e. on how the deterministic target emerges from the local noise.Moreover, it is possible to find an approximation for the asymptotic cases: • When µ c tends towards infinity, that means a target whose value is infinitely greaterthan that of the speckle, we find a Normal law N ( µ c , µ √ ), and the coefficient of vari-ation tends towards 0. Moreover, we can find the asymptotic approximation γ = √ λ • When µ c or λ tends towards 0, we find the Rayleigh law which followed the fullydeveloped speckle. Then the coefficient of variation will correspond to the case offully developed speckle γ = 0 . . . It is possible to find the asymptotic relation γ = 0 . − . λ for small λ , satisfactory for λ < . g ( m , m ) by performing a first-order limited expansion of thefunction around the values m , and m . . This method applied to the coefficient of variationexpressed in terms of m and m in Eq. 1 leads to find the following formula, as alreadyproposed in [13]: var ( γ ) = 14 N m − m m + m m − m m m m ( m − m ) , (4)where N is the number of images in the sequence.It is of no use to express this variance for the coefficient of variation of Rice’s law:the even moments involving exponential functions and the odd moments of the modifiedBessel functions of the first kind, no simplification can be expected. However, it has beendemonstrated that the moment of order n is proportional to µ n , and we can deduce fromthis relation that the variance of the estimator of the coefficient of variation depends onlyon the parameter λ and N , and can be written: var ( γ ) = 1 √ N f ( λ ) (5)In all cases, the coefficient of variation associated with a Rice law is statistically lowerthan that of the corresponding Rayleigh law, and this, especially since the ratio µ c /µ islarge. A speckle area is typically a forest or area of bare soil. In SAR images, it is commonly ac-cepted that the amplitude of a texture-free speckle distribution follows a Rayleigh-Nakagamilaw [7]: RN [ µ, L ]( u ) = 2 µ √ L Γ( L ) ( √ Luµ ) L − e − ( √ Lµ u ) (6)where µ is the shape parameter and L the Equivalent Number of Looks (ENL).Let us note that the case of the fully developed speckle corresponds to the case L = 1,and to Rayleigh’s law. Rayleigh Nakagami law generalizes the framework of the studyto multilooked speckle, as in the case of Sentinel 1 GRD data. Note that Rice’s law onlyaddresses the fully developed speckle as observed in the SLC data. The case of a deterministictarget in multilooked data has not been studied. Its formalism would indeed be even morecomplicated: it would be necessary to study the spatial averaging of an area containinga deterministic component in a pixel and surrounding pure speckle. This study is out ofthe scope of this paper. However, it is useful to show, as part of the speckle, how the L parameter can influence the behavior of the coefficient of variation.Generally, the parameters of this law are estimated spatially from a homogeneous area.In our approach, we are interested in temporal statistics. We then assume that, for a pixelbelonging to a Rayleigh Nakagami law speckle area, and having not perturbed by a change,the various realizations of amplitudes over time of this pixel also follow a Rayleigh Nakagami6aw. In a way, we consider that the hypothesis of spatial ergodicity can be transformed intoa hypothesis of temporal ergodicity. In this case, we will speak in the following of stablespeckle .From Rayleigh Nakagami law, it is possible to derive the expression of empirical moments.We have in particular [13]: m = µ Γ( L + ) √ L Γ( L ) and m = µ , (7)which allows to find the following expression of the coefficient of variation, defined as theratio of the standard deviation and the average: γ = σµ = (cid:112) m − m m = (cid:115) Γ( L )Γ( L + 1)Γ( L + 1 / − the coefficient of variation willhave the same value for all stable speckle zones, whatever the average amplitudeof this speckle. In the particular case L = 1 we find γ = 0 . L : varγ ) = 14 N L Γ( L ) (4 L Γ( L ) − L Γ( L + ) − Γ( L + ) )Γ( L + ) ( L Γ( L ) − Γ( L + ) ) (9) In the particular case L = 1, we have var ( γ ) = 0 . /N , σ ( γ ) = 0 . / √ N respec-tively for the variance and the standard deviation.The application of a multilook on a SAR image has the effect of modifying the L param-eter. The Sentinel-1 GRD data are delivered with a displayed equivalent number of looksof L = 4 .
9, calculated for the theoretical value of the middle of swath and in the middleof orbit [18]. For a speckle area, this value of L gives γ = 0 . γ ) = 0 . /N and σ ( γ ) = 0 . / √ N .This formula shows a second attractive property on this parameter: the standarddeviation of the coefficient of variation decreases with the number of images N in √ N .We can thus imagine that, on temporally stable speckle areas, the coefficient of variationwill be constant, with a variance decreasing with the number of images in the sequence.A last significant result is that the coefficient of variation, calculated on amplitudesaccording to a Rayleigh Nakagami law, seems to follow a Rayleigh Nakagami law again.This property is still a postulate, verified for many simulations and statistical analyses. We have theoretically analyzed the coefficient of variation (CV) for the temporal am-plitude profile of a pixel, modeled as the superposition of a fully developed speckle and adeterministic component. 7he following properties have been demonstrated:- The higher the deterministic component, the lower the Coefficient of Variation. It tendsto 0 as the relative power between the deterministic component and the speckle increases.It depends only on this ratio.- When the deterministic component tends towards 0, then CV tends towards 0.52. Weare then in the presence of a pure decorrelated speckle. In this case, one finds empiricallythat the distribution of the CV corresponds to a Rayleigh law.- In all cases, the variance of the estimator is proportional to 1 /N where N is the numberof images in the time-series.
3. Change detection criteria
Thanks to the properties highlighted previously, we propose several decision criteria forthe change detection, based empirically on the coefficient of variation.
Let S be a sequence of N co-registered SAR images acquired at time { t , t , . . . , t N } .In this paper, a purely temporal approach is proposed. That means that we develop criteriathat are defined for each pixel independently, by considering only its temporal profile: x ( t ) = { A ( t ) , t = t , ..., t N } where A ( t ) is the amplitude signal for this pixel in the imageacquired at time t .From this profile, a first intention is to determine if there is a change, of any naturewhatsoever, that means with at least one break. A second intention is to detect a specifickind of change. However, we do not determine the different break dates, unlike intended in[4]. Therefore, we will propose now the coefficient of variation as generic change detection,and other criteria for specific change detection.Change detection can be formulated as a classical detection problem using a hypothesistesting, with two hypothesis, H =”no change” against H =”a change”. Then the math-ematical expression of the test corresponds simply to the thresholding of the criterion mapfollowing: f ( x ( t )) H ≶ H λ (10)with f the criteria function and λ the threshold value.We describe three different functions in the next subsections for generic and specificchanges. In section 2, we have illustrated that the coefficient of variation is even lower when thesignal is stable. Our first generic change detector is simply this criterion. In this case,function f is defined by: f ( x ( t )) = γ ( x ( t )) = σµ (11)8ith µ = 1 N N (cid:88) t =1 x ( t ) , σ = (cid:118)(cid:117)(cid:117)(cid:116) N N (cid:88) t =1 x ( t ) (12) We call a point-event or one-time event, a time profile for which a high-backscatteringdeterministic target is seen at a single date in the profile. In practice, this will correspondto the case of detection of a vehicle. For example, it is the case for boats on the ocean as theprobability of having different boats in the same place at different moments is small, evenin a long sequence. It is also the case for vehicles in the desert.The empirical following criterion is proposed. It is defined as the ratio of the coefficientsof variation, the first one computed on the profile without the minimal amplitude value, andthe other one computed for the profile without its maximal value.The expression is: f ( x ( t )) = γ ( x ( t ) t ∈{ ...N }\{ t max } ) γ ( x ( t ) t ∈{ ...N }\{ t min } ) (13)where t min and t max are the indices on which the signal is minimal or maximal respectively.Note that a sequential version of this criterion is directly obtained as a special case: it isto consider the ratio of coefficient of variation for the profiles respectively without the firstand the last date obtained. This makes it possible to obtain a detector of an event arrivingfor the last date of acquisition. In this case, the threshold will be applied to: f (cid:48) ( x ( t )) = γ ( x ( t ) t ∈{ ...N } ) γ ( x ( t ) t ∈ ...N − ) (14)We also propose to compare the f criterion with the one computed only on the amplitudemeans: f ( x ( t )) = µ ( x ( t ) t ∈{ ...N }\{ t max } ) µ ( x ( t ) t ∈{ ...N }\{ t min } ) (15) Among the changes that do not belong to the previous category, we will find othercategories, such as chaotic changes in radiometry, or the presence of deterministic targetsover more extended periods.The first ones are observed in practice in all agricultural areas. The coefficient of variationis sensitive to these variations related to seasonal changes, moisture changes, or plant growth.These variations can be considered as randomly spread over time. Because changes relatedto agricultural surfaces are not of interest in some applications, we are more interested inthe second category of changes. Therefore, they correspond to changes related to objects9uch as infrastructures construction or destruction. The typical profiles generally includesteps, that is to say, several dates contiguous with a strong deterministic signal.We then propose empirical criteria dedicated to this case. The idea is to split the sequenceinto two parts, to compute the ratio between the coefficient of variation on both sides. Then,we take advantage of the temporal redundancy by scanning the different possible cut-offpoints and by averaging the results thus obtained.We expect then to have two typical behaviors: • If the temporal evolution is simply a radiometric hazard such as those encounteredwith agricultural surfaces, then the coefficients of variation between the different sub-parts of the profile remain statistically equivalent. It is then expected to obtain a ratioclose to 1. • On the other hand, in the presence of a change as infrastructure construction or aparticular object, several cuts of the profile will lead to finding coefficients of variationsignificantly different. The ratios should thus have values deviating from 1.Let p be the variable corresponding to the break date of the sequence. We obtain two sub-parts of the temporal profile noted x p ( t ) and x p +1 ( t ). Then the calculation of the coefficientof variation from these profiles is respectively noted γ p and γ p +1 . To have a significantmeaning for the ratio measure, we must respect the statistical constraint for calculating thecoefficient of variation, which is to have a sufficient number of images. If we note M theminimum number of images, then p must vary between M and ( N − M ).Then, the function f for the hypothesis test is defined by : f ( x ( t )) = 1 − N − M + 1 N − M (cid:88) p = M min( γ p γ p +1 , γ p +1 γ p ) (16)In this expression, we have a min operator to have a resulted value between 0 and 1.Moreover, we take one minus the obtained value in order to have low values in case of nochange and close to 1 in case of change.Once again, we propose to consider the same idea based on mean amplitudes instead ofthe coefficient of variation: f ( x ( t )) = 1 − N − M + 1 N − M (cid:88) p = M min( µ p µ p +1 , µ p +1 µ p ) (17)
4. Analysis of performance thanks to statistical simulations
To evaluate the detection performance based on the previous criteria, we rely on manysimulations. These consider different kinds of change. Therefore, we are now analyzingdifferent kinds of profiles. 10 igure 2: Examples of simulated profiles with Rayleigh and one-point event, for two contrast ratios
We consider here the profiles containing, for the most part, uncorrelated speckle, andfor which we have at a given date, the addition of a deterministic target. In this context,we make the number of images of the temporal stack N vary, and the ratio between theamplitude of the target µ c and the mean amplitude of the speckle µ . Note that µ is notprecisely equal to the µ parameter, but is deduced by Eq. 7. Given the particular dynamicsof the radar, we vary this contrast linearly in decibels. Two examples of such profiles aregiven in Fig. 2.In order to estimate the performance of a given criterion, we simulate a large number ofprofiles. For each category of changes, we simulate K = 10 reference profiles without anychange and K profiles with a change. The investigated criterion is calculated for these twopopulations. This computation enables to achieve K estimation of this criterion, with andwithout change. From these K realizations, the empirical distributions of this parameterwith and without change are estimated. It is thus classical to compute a Probability ofDetection/False Alarm curve for this parameter and this scenario by varying the thresholdof the studied criterion.For the one-point event, we have compared the three criteria f (coefficient of variation), f (ratio of the coefficient of variation, without min and max value) and f (ratio of amplitudemeans), by fixing PFA to 0.1%. Fig. 3 shows a colored representation of the Probability ofdetection for the three criteria: Red channel corresponds to the f criterion, Green channelto f and Blue channel for f . Black color corresponds to a zero PD for all criteria whereasWhite color means PD=1 for all criteria. • For PFA=0.1%, a target can be detected as soon as its contrast is more than 8db.Such a profile example is given in the figure 2. • for N <
20, the coefficient of variation f and the intensity mean ration f should bepreferred. 11 igure 3: Probability of detection of f , f , f , for False Alarm 0.1 %, for profiles containing a target at asingle time, for varying number of images and varying contrasts • for N >
20, ratio criteria f and f have better performances than f for contrasts lyingbetween 8 and 13dB; for contrasts beyond 13dB, all criteria have same performances. We consider here profiles mixing two different populations of speckle, with a speckle oflevel given in proportion P , and a speckle of higher level, for the rest of the profile (proportion1- P ). This type of profile can correspond, for proportions close to 1/2, to profiles of typesof agricultural plots, for which different backscattering levels may appear. Note that when P tends to 0 or 1, we find ourselves in cases of profiles closer those corresponding to one-offevents.In this case, we have again set the number of images to N = 100, we then varied P , thepercentage of dates at which speckle 1 is present relative to speckle 2 and the ratio betweenthe two average amplitudes of speckle.In the figure 4, two examples of profiles are given by setting P = 50, for two contrastvalues between speckle: 15dB and 8dB. These examples illustrate the difficulty of detectingwhether a population corresponds in practice to a single Rayleigh population or a mixtureof two populations.The figure 5 shows the sensitivity of our different f , f and f criteria to this type ofmixture profiles. Only the f criterion is totally insensitive to this mixture, unless it occursonly on one date, which brings us back to the previous case of the point-event target.The coefficient of variation ( f ) remains very sensitive to a mixture of this type. Fornatural areas with radiometric evolutions, it will lead to change detection. The f criterionremains a little less sensitive than f but will still detect such events as a change in the caseof mixing over low durations. 12 igure 4: Examples of simulated profiles with mixture of two Rayleigh distributions, for two contrast ratiosFigure 5: Probability of detection of f , f , f , for mixture of two Rayleigh speckle, for False Alarm 0.1 %,for varying percentage of duration ratios and varying contrasts. igure 6: Examples of simulated profiles with Rayleigh and one step event, for two contrast ratios We now consider profiles containing speckle, on which a deterministic type target issuperimposed for a date proportion equal to P . In this case, we have set the number ofimages to N = 100; we make P vary, the percentage of dates at which the deterministiccomponent is present, and the ratio between the magnitude of the target and the amplitudeof the speckle. This type of profile typically corresponds to a site for which there is thepresence of a deterministic signal on several dates. It is typical of building construction forexample.The figure 6 illustrates two examples of a profile with a deterministic component over25% of the observation duration, for two contrast values.In a first simulation illustrated in figure 7, probability of detection is represented for thethree criteria f , f and f for varying proportions P and varying contrasts. These criteriaare independent of the order of the values of the profile. Results are very similar to the onespresented in figure 5. • The criterion f is the only one which is specific to a very time-limited change. • The criterion f can detect step with duration up to 20% of total duration. It is ofless good specificity. • The general criterion f detect only step with duration less than the half observationperiod. For longer events, its failure is probably because the deterministic componentof the step becomes predominant in the profile, which makes the coefficient of variationdecrease again.In another simulation illustrated in figure 8, we consider different fixed values of theproportion P of points that constitutes a step in the profile, and we study the impact of the14 igure 7: Probability of detection of CV and ”one-point” detectors, for False Alarm 0.1 % corresponding toa step signal for varying length of times and contrasts date of the beginning of this step and the contrast. Probability of detection is representedfor the three criteria f , f and f .These simulations lead us to the following observations: • Sensitivity to the relative duration of the step ( P ) – the f criterion only manages to detect steps of durations lower than half of theduration of observation. – the criterion f can detect only steps with durations greater than half the durationof observation. Beyond a certain duration (P¿ 80%), it is the only criterion whichmanages to detect the step – the criterion f is the most versatile compared to the duration of the step. • Sensitivity to the date of the beginning of the step – only the criterion f does not depend, by construction, on the position of thestep. – the criteria f and f depend on this position and detect it all the more difficultsince it is centered temporally, and also for particularly brief events ( P ¡25%) orparticularly long events ( P ¿75%). – for high contrast, more than 10 dB , there is no more impact on performances15 igure 8: Probability of detection of f , f and f for False Alarm 0.1 %, for steps with 6 different lengthsof time, for varying step position and contrasts. . Performance Evaluation on real SAR data Now that our criteria have been evaluated through statistical simulations, we are con-sidering real data validation cases.
The first site of interest is chosen on the Saclay area (near Paris, France) between June2015 and June 2017. Indeed, this site includes many construction sites related to the devel-opment of the University Paris-Saclay. Construction dates are available thanks to the localauthorities.The considered area is about 15 km x 12 km. A precise Ground Truth Database has beenestablished, by using two geospatial vector databases: a topographic database produced byIGN (French National Geographic Agency), and OSM database, corresponding respectivelyat dates before and after the observed time range. We kept several standard semantic classesto help interpret the results: water surface, crops, forests, roads, and buildings. Both vectorfiles have been converted into raster images by using the GDAL rasterize utility, on a 10m x10 m common georeferencing grid, and have been subtracted from one another. Then, all thechanges found were manually validated or invalidated using past versions of High-Resolutionoptical images on the Google Earth timeline. The resulting ground truth in Fig. 9 showsdifferent classes. Among them, parking areas, building construction or destruction, andground changes are assigned to a generic change class used in the performance evaluation.The numbers of change for each class are indicated in table 1.Building 1 Parking area 2 Ground area140 111 42
Table 1: Number of changes for each category
This site is observed by the Sentinel 1 sensor. Change detection methods have beenapplied on SLC (Single Look Complex) images and not on GRD (Ground Range Detected)ones. We have made this choice for several reasons: • Statistical properties of speckle are preserved in the whole SLC image and do notrequire to take into account the effect of multilooking. • full polarimetric information is preserved • in future works, coherent methods could be applied and compared in the presentresults.We have then computed all the criteria maps in the SLC reference grid, before exportingthem to the Lambert Conformal Conic projection coordinate system of the Ground Truth.In order to achieve this, the dense non-rigid transformation between SLC and GRD images17 igure 9: Ground truth established on a test site around Saclay, France igure 10: REACTIV method applied to an area of the Grizzly Bay, San Francisco (Left) and correspondingmanual data entry for Ground Truth (Right), 600x600 pixels has been computed thanks to the GeFolki algorithm [2]. Then this transformation can beapplied to any product computed from the SLC time-series.In order to take into consideration also one-event targets, a second validation scenariohas been considered, with 68 UAVSAR images in the Grizzly Bay, around San Francisco.Several boats are available in the images. They belong to the Fifth Reserve Fleet whichis docked off the coast of Benicia. The ground truth has been established by manuallysegmenting the vessels on each image of the stack. For some pixels, an overlap of two boatspresent at two different dates may occur. In the following, we will estimate the detectionperformance only on the parts outside these recoveries, to correspond precisely to the ”singleevent” case.Note that in this data set, we do not strictly have a Ground Truth since we do not useinformation external to our images; this manual data entry is intended to observe statisticalbehaviors. The Saclay area contains a significant proportion of cultivated fields and vegetation, inaddition to urban areas. The ground truth mainly involves changes in the building elements,some elements of ground modifications and the parking areas. In this context, we will onlyevaluate the generic and specific detection of step signals. Besides, we compare our resultswith those obtained with two other reference methods: • the reference method of the literature on the detection of a change in SAR time seriesby the so-called ”omnibus” method [15], • a state of the art method of change detection by deep learning [6], applied to Sentinel2 optical images on the same area. This method has been tested only on buildingchanges for which it has been trained. 19he first method consists of a statistical homogeneity test based on a likelihood-ratiotest. Authors propose two versions of this method. A first version considers the test of globalhomogeneity of the complete temporal profile. It is possible to demonstrate that in this case,the test leads to threshold a criterion proportional to the ratio between geometric mean andarithmetic mean of the intensity signal. When polarimetric information is available, thedeterminant of the polarimetric covariance matrix replaces the single-channel intensity.In a second version of the code, several homogeneity tests are carried out iterativelyon subsets of the profile: the intensity considered at date j is compared with the intensitypopulation estimated for all previous dates without breaks. This second version of the codeis more comprehensive as it finds the location of all probable breaks. However, it is muchmore expensive in computing time. In this paper, the first version will be called omnibus method and the second version of the full omnibus methodRegarding the deep learning method, the comparison must be made with much morecaution: the detection of change has been done in a bi-date way, between the start date andthe end date. We took care to select images without cloud cover. However, it is likely thatmany changes over a limited period within the range cannot be seen in this context. Onthe other hand, the learning has been done mainly on changes relating to new or missingbuildings, and therefore, only the performances related to this type of change are consideredlater. This comparison aims to demonstrate the difficulty of the task of detection of a changefrom optical images, including with the most advanced algorithms today in deep learning,compared to that of detection from radar images.We obtain each point of the resulting ROC curves by fixing a threshold on one criterion.For a given threshold, the Probability of False Alarm (PFA) is given by the percentage ofpixels detected that are not lying in the change class of the ground truth. The Probabilityof Detection (PD) is the percentage of change class objects that are detected; we consideran object as detected as soon as a minimum of 5% pixels of the entire object has beendetected. We justify this choice by the very sparse nature of a radar signal. For a givenbuilding footprint, only a few points are bright points. Thus, we can not hope to detect abreak on all the pixels of the considered footprint. Lastly, the Precision is the percentageof pixels declared as detected that corresponds effectively to a change in the Ground Truth.We obtain two types of curves, the (PD, PFA) curve and the (Precision, Recall) curve withthe Recall equal to Probability of Detection (PD).Figure 11 shows the performance for the detection of changes in all classes. Overall, wespecify above all that we focus our analysis on low false alarm rates. Indeed, we believethat few operational scenarios can accept 5% False Alarm, mainly when they are denselydistributed among the whole image. With this point of view, we can not expect a detectionrate higher than 70%. For general changes, the cumulative ratio of CV ( f ) and the fullomnibus method obtain the best performances. Nevertheless, when we restrict the changedetection to the building class, the figure 12 shows that several methods achieve a muchbetter result, such as 70% for less than 1% as False Alarm Rate. This analysis indicatesthat the main difficulty relies on park areas and ground changes, that cannot be robustlydetected on these Sentinel 1 images.For change detection restricted to buildings in figure 12, except the deep learning method,20 igure 11: Generic Change Detection (Saclay site, France), PD/PFA on the left, Precision/Recall on therightFigure 12: Case of construction site detection (Saclay site, France) - PD/PFA performances / Preci-sion/Recall performances all methods obtain similar performances for low False Alarm Rates, and two methods arebetter for higher false alarm rates, that are the cumulative ratio of the coefficient of variation f and the cumulative ratio of mean amplitude f .When qualitatively analyzing the criteria maps, we can see that most false alarms aredue to high sensitivity to changes in an agricultural area, except for f . Depending onweather conditions and crop growth, it is likely that the electromagnetic mechanisms andtherefore the backscattering levels vary significantly from one date to another. In this case,simulations have proven that cumulative ratio methods f are less sensitive to this kind ofchanges than f and f . On the other hand, the main problem with f highlighted by ourprevious simulations is the difficulty in detecting events shorter than half the duration ofobservation. Thus, the f and f criteria have complementary behaviors for comparableperformance.Finally, we note the performance significantly worse with the method of deep learning,21 igure 13: case of construction site detection of size ¿ 20px (Saclay site, France) PD/PFA performances -Precision/recall highlighting the difficulty to detect changes in optics.In a second step, the performances have been calculated by removing any zone of size lessthan 20 pixels from the ground truth and the change detection results. This restriction affectsmainly the performance on building detection, as shown in figure 13 where the improvementis significant, mainly for two methods, the cumulative ratio of CV ( f ) and the full omnibusmethod.In summary, two methods show excellent performance for change detection on largeenough buildings: the cumulative ratio of CV ( f ) and the full omnibus method. Never-theless, remember that the full omnibus method is time-consuming (50 minutes) while theother method is high-speed (4 minutes) for the change detection on radar series composedof 2 (polar) x 64 images of size (1133 x 3205).CV or Omnibus Ratio of CV or Amplitude Omnibus (full)6 seconds 4 minutes 50 minutes Table 2: Time calculation on Intel(R) Core(TM) i7-4710HQ CPU @ 2.50GHz
On this second set of data, we address a very different case from the previous one: thatof one-point events, with an example of boat detection.The performance analysis conducted as previously at the object level would have nomeaning here: firstly because some criterion thresholds allow us to obtain a Probability ofdetection equal to 100% for a null probability of false alarm; secondly, because we do nothave Ground Truth external to the data.For these reasons, we analyze performances on this real data set, statistically at the pixelscale. 22 igure 14: Performances of change detection for ships at the pixel level, on full polarimetric UAVSAR dataset
We have selected all the pixels belonging to one single boat without superposition, thoserepresented in gray in the figure 10. The histograms of all these pixels were calculated forthe three criteria f f and f , and compared to those of a sea area. The ROC curvescorresponding to these comparisons of distributions are presented in figure 14.The f and f criteria are the ones that lead to the best performances. Be careful,note that the y-axis starts at 80% detection. These two specific criteria are statisticallybetter than the generic criteria for the coefficient of variation and the polarimetric omnibuscriterion computed on the diagonal covariance matrix.
6. Conclusion
In this paper, change detection in SAR time-series has been studied. We have proposedseveral criteria based on the temporal coefficient of variation.This parameter has indeed revealed useful statistical properties: its average values andits variances do not depend on the scale parameter of the speckle. In the case of the presenceof a deterministic component, it depends only on the contrast between this and the clutter.Statistics properties of this criterion have made it a key parameter for deploying robust andextremely fast change detection strategies.A first one, the coefficient of variation itself, is dedicated to detection for any change.Also, we proposed to use it as a generic change detector. Then, since the notion of changecorresponds to several types of events, we have proposed other detectors based on ratios, inorder to adapt to either short period event or one-point events such as a vehicle, or moreextended events such as construction sites.These criteria were first analyzed using statistical simulations. These simulations provedthe great genericity of the coefficient of variation to detect any change, but also the specificityof the other criteria to address specific changes. Then these criteria were evaluated on realdata. 23or step type changes, the specific criteria lead to excellent performance, at least equalto that of the state-of-the-art method ”full omnibus” which is much more expensive incomputing time. Performances are also much better than those that would be obtained usingoptical data on a similar scenario. Moreover, the two proposed criteria are complementary:one is likely to respond better to short-steps, in return is more sensitive to changes in allagricultural areas.For ”one-event” changes, the validation has been done on a manual statistical analysisof real high-resolution data containing a large number of ships.Here again, the two specific criteria lead to the best performances. The simulations tendto say that the criterion based on the CV ratio is even more specific to these short periodevents. This conclusion allows us to consider using this parameter when detecting a newevent with each new acquisition.In the future, future work will focus on proposing a fusion of the various criteria devel-oped. Ideally, a complete schema could go as far as categorizing change profiles. In parallel,we will analyze the contribution of polarimetric information.
Aknowledments we gratefully acknowledges Behnaz Pirzamanbein for her kindly help on the full omnibusmethod, and her team for sharing the code. We are very grateful for Rodrigo Dault, Ph.D.student in Onera, for having tested his algorithm on the Sentinel 2 images over Saclay. Wefinally thanks for the long-term discussions at the origin of this work our colleagues BeatricePinel-Puyssegur and Jean-Michel Lagrange from CEA-DAM.UAVSAR data are courtesy NASA/JPL-Caltech, BD TOPO courtesy IGN (France),Sentinel 1 data courtesy ESA.The APC was funded by Onera, as part of the MEDUSA Research project.The formal analysis in this paper has been conducted by E.K. and J-M. N. Conceptual-ization of this paper has been proposed by E.K.; the project administration and writingo-riginal draft preparation has been proposed by F.J. The methodology, softwares, validation,writingreview and editing have been performed by F.J. and E.K.24 bbreviations
The following abbreviations are used in this manuscript:CDM Change Detection MatrixCV Coefficient of VariationDInSAR Differential Interferometry SARENL Equivalent Number of LooksGRD Ground Radar DetectedHSV Hue Saturation ValueIGN Institut Gographique National (French National Geographic Agency)MIMOSA Method for generalized Means Ordered Series AnalysisONERA Office National D’Etudes et de Recherches Arospatiales (The French Aerospace Lab)OSM Open Street MapPD Probability of DetectionPFA Probability of False AlarmREACTIV Rapid and EAsy Change detection on Time-series using coefficient of VariationRGB Red Green BlueROC Receiver Operating CharacteristicSAR Synthetic Aperture RadarSLC Single Look Complex
References [1] Donato Amitrano, Gerardo Di Martino, Antonio Iodice, Daniele Riccio, and Giuseppe Ruello. Rgb sarproduct exploiting multitemporal: general processing and applications. In
Analysis of MultitemporalRemote Sensing Images (MultiTemp), 2017 9th International Workshop on the , pages 1–4. IEEE, 2017.[2] Guillaume Brigot, Elise Colin-Koeniguer, Aur´elien Plyer, and Fabrice Janez. Adaptation and evaluationof an optical flow method applied to coregistration of forest remote sensing images.
IEEE Journal ofSelected Topics in Applied Earth Observations and Remote Sensing , 9(7):2923–2939, 2016.[3] L Bruzzone, M Bertoluzza, and Francesca Bovolo. A new paradigm for the exploitation of the semanticcontent of large archives of satellite remote sensing images. In
Conference on Big Data from Space ,pages 114–117. JRC, 2017.[4] Knut Conradsen, Allan Aasbjerg Nielsen, and Henning Skriver. Determining the points of change intime series of polarimetric sar data.
IEEE Transactions on Geoscience and Remote Sensing , 54(5):3007–3024, 2016.[5] Shiyong Cui.
Spatial and temporal SAR image information mining . PhD thesis, University of Siegen,2014.[6] Rodrigo Caye Daudt, Bertr Le Saux, and Alexandre Boulch. Fully convolutional siamese networksfor change detection. In , pages4063–4067. IEEE, 2018.[7] J.W. Goodman. Some fundamental properties of speckle.
JOSA , 66(11):1145–1150, 1976.[8] Flore Harl´e.
D´etection de ruptures multiples dans des s´eries temporelles multivari´ees: application `al’inf´erence de r´eseaux de d´ependance . PhD thesis, Grenoble Alpes, 2016.[9] M. Kendall and A. Stuart. The advanced theory of statistics. Vol. 1: Distribution theory.
London:Griffin, 1977, 4th ed. , 1977.[10] Thu Trang Lˆe.
Extraction d’informations de changement `a partir des s´eries temporelles d’images radar`a synth`ese d’ouverture . PhD thesis, Grenoble Alpes, 2015.
11] Gr´egoire Mercier. Progressive change detection in time series of sar images. In
Geoscience and RemoteSensing Symposium (IGARSS), 2010 IEEE International , pages 3086–3089. IEEE, 2010.[12] Javier Muro, Morton Canty, Knut Conradsen, Christian H¨uttich, Allan Aasbjerg Nielsen, HenningSkriver, Florian Remy, Adrian Strauch, Frank Thonfeld, and Gunter Menz. Short-term change detectionin wetlands using sentinel-1 time series.
Remote Sensing , 8(10):795, 2016.[13] JM Nicolas. Application de la transform´ee de Mellin: ´etude des lois statistiques de l’imagerie coh´erente.
Rapport de recherche, 2006D010 , 2006.[14] JM Nicolas. Un nouveau formalisme pour la loi de Rice, 2018. https://perso.telecom-paristech.fr/nicolas/jmnicolas rice 2018D003.pdf.[15] Allan A Nielsen, Knut Conradsen, and Henning Skriver. Change detection in a time series of polari-metric sar data by an omnibus test statistic and its factorization (conference presentation). In
Imageand Signal Processing for Remote Sensing XXII , volume 10004, page 1000411. International Society forOptics and Photonics, 2016.[16] Allan A Nielsen, Knut Conradsen, Henning Skriver, and Morton J Canty. Change detection in a seriesof sentinel-1 sar data. In
Analysis of Multitemporal Remote Sensing Images (MultiTemp), 2017 9thInternational Workshop on the , pages 1–3. IEEE, 2017.[17] Joshua Rutkowski, Morton J Canty, and Allan A Nielsen. Site monitoring with sentinel-1 dual polar-ization sar imagery using google earth engine.
Journal of Nuclear Materials Management , 46(3):48–59,2018.[18] Sentinel 1 Team. Sentinel 1 User Handbook, 2013. https://sentinel.esa.int/documents/247904/685163/Sentinel-1 User Handbook.[19] Xin Su, Charles-Alban Deledalle, Florence Tupin, and Hong Sun. Norcama: Change analysis in sartime series by likelihood ratio change matrix clustering.
ISPRS Journal of Photogrammetry and RemoteSensing , 101:247–261, 2015., 101:247–261, 2015.