[1]Arunkumar Govindarajan

[2]Subhankar Chattoraj [1]Aarthi Scans & Labs, Chennai, Tamil Nadu, India
2]Qure.ai, Mumbai, Maharashtra, India

# Identification of Hemorrhage and Infarct Lesions on Brain CT Images using
Deep Learning

*<EMAIL_ADDRESS>Arjun Agarwal<EMAIL_ADDRESS><EMAIL_ADDRESS>Dennis Robert<EMAIL_ADDRESS>Satish Golla<EMAIL_ADDRESS>Ujjwal Upadhyay<EMAIL_ADDRESS>Swetha Tanamala<EMAIL_ADDRESS>Aarthi Govindarajan<EMAIL_ADDRESS>*  [

###### Abstract

Head Non-contrast computed tomography (NCCT) scan remain the preferred primary
imaging modality due to their widespread availability and speed. However, the
current standard for manual annotations of abnormal brain tissue on head-NCCT
scans involves significant disadvantages like lack of cutoff standardization
and degeneration identification. The recent advancement of deep learning-based
computer-aided diagnostic (CAD) models in the multidisciplinary domain has
created vast opportunities in neurological medical imaging. Significant
literature has been published earlier in the automated identification of brain
tissue on different imaging modalities. However, determining Intracranial
hemorrhage (ICH) and infarct can be challenging due to image texture, volume
size, and scan quality variability. This retrospective validation study
evaluated a DL-based algorithm identifying Intracranial hemorrhage (ICH) and
infarct from head-NCCT scans. The head-NCCT scans dataset was collected
consecutively from multiple diagnostic imaging centers across India. The study
exhibits the potential and limitations of such DL-based software for
introduction in routine workflow in extensive healthcare facilities.

###### keywords:

Computer-aided-diagnostic solution, Non-contrast CT, Intracranial hemorrhage,
Infarcts, Deep learning, Clinical evaluation

##  1 Introduction

In Cognitive Neuroscience, Neuropsychological investigation of stroke patients
is widely utilized in advancing our knowledge of brain functions. The
considerable insight into the relation of the brain function to its anatomy
has been determined via correlation analysis between physical brain damage and
impaired behavior [1][2][3]. The stroke topology can be broadly classified
into two types: $1)$ Intracranial hemorrhage (ICH), the rupture blood vessel
within the brain which causes bleeding. The common factors related to the
cause of ICH are advanced age, heavy alcohol usage, and high blood pressure
(hypertension) [4]. As per some recent studies, although ICH accounts for
10–15% of all stroke-related deaths, over the last thirty years, the mortality
and morbidity have not changed, particularly in developing countries [5]. $2)$
Ischemic stroke or infarct, is interruption of blood flow due to blood clot.
Infarct is generally caused by the buildup of plaques (atherosclerosis) over
time in the arteries. Globally, over 13.7 million individuals have a stroke
each year, of which approximately 70%, i.e., 9.5 million, are infarct [6].
Presently, mapping of the stroke lesion is regularly done using Computed
tomography (CT) and magnetic resonance imaging (MRI). The MR (T1- weighted and
T2- weighted) anatomical images are acquired as a part of routine practice for
stroke patients. In stroke suspected patients with negative CT scans, MRI can
also be performed. After the first few hours of onset, the ischemic stroke can
be identified using the MRI. Additionally, the differentiation of irreparably
damaged brain tissue and the tissue at risk due to infraction can be diagnosed
using the MRI. However, CT is the preferred imaging modality over MRI in acute
stroke care units and clinical trials due to the reduced exclusion criteria
compared to MRI, affordability, speed, and accessibility [7]. In CT,
hemorrhage is percieved as the bright region (hyper-dense) exhibiting sharp
contrast and infarct as dark region (hypo-dense) depending on the time
progressed after the onset.

The manual annotations of abnormal brain tissue by trained neuroradiologists
is currently the present standard method for lesion identification [8].
However, the manual annotations of abnormal brain tissue have many
disadvantages [9]. $1)$ Lack of cutoff standardization:, There is no standard
protocol for explicit cutoff, particularly around the ventricles, to
differentiate lesioned and non-lesioned tissues; as a result, this approach
produces large variability and lacks reproducibility across operators. $2)$
Degeneration identification: The stroke-induced degeneration occurring in
chronic stroke patients outside the lesion is not captured in the standard
manual annotations process, even though a significant clinical impact on
patients is caused due to the stroke-induced degeneration. The recent
advancement of deep learning based computer aided diagnostic (CAD) models in
medical imaging and signal processing can significantly assist in overcoming
the existing challenges [10][11][12][13][14]. In addition, the manual editing
combined with an automated detection solution of hypo- or hyper-dense regions
that remains under operator supervision and can assist in overcoming the
present challenges [15]. More recently, a study using large CT datasets to
remove the inter-subject variability in brain lesion characterization using an
automated approach was proposed [16]. Several state-of-the-art algorithms have
been proposed for lesion segmentation in MR images over the past few years,
but very few have been developed to address stroke lesions on CT scans. Most
of the earlier work published to validate automated solutions was directed
toward identifying ICH. As the ICH appears bright in CT scans, developing an
automated solution based on supervised or unsupervised learning algorithm or
extracting morphological features from labeled images to differentiate between
true lesioned and non-lesioned tissues is less challenging [17] [18]. Infarct
identification, on the other hand, is a less popular problem statement among
researchers compared to ICH detection due to its challenging nature. To
address this issue very recently, a rule-based approach based on seeded
region-growing algorithms was proposed via extracting hand-crafted features
such as relative position for an axis of symmetry, texture, and brightness
[19]. However, the primary disadvantage of this study is that the seeded
region-growing algorithms may not be able to define the boundaries of the
stroke region distinctively.

In this study, we have evaluated an Artificial Intelligence (AI) based
automated CAD algorithm based on deep learning, capable of identifying ICH and
infarct on Head-Non-contrast Computed Tomography (Head-NCCT) scan. The
solution has been earlier validated on detecting ICH on Head-NCCT scan images
[14]. The Institutional Review Board (IRB) has approved the proposed
retrospective study. We demonstrated the effectiveness and validity of the
automated CAD solution in detecting ICH infarct and quantifying infarct on
Head-NCCT scan. Our proposed validation will provide a rapid and efficient
tool for both research and clinical application. It will assist in the broader
adaptation of automated CAD solutions at extensive clinical facilities.

##  2 Material and Methods

The study was a HIPAA-compliant retrospective study with Institutional Review
Approval (IRB) from Royal Pune Independent Ethics Committee (RPIEC) (IRB No.
RPIEC240123). Informed consent was obtained from all participants. All methods
were carried out in accordance with relevant guidelines and regulations.

The primary objective was to evaluate the commercially available deep
learning-based algorithm qER (Qure.ai Technologies, Mumbai, India) in terms of
Area Under the Receiver Operating Characteristics Curve (AUC) in triaging
Head-NCCT scan in detection and quantification of infarcts. It was estimated
that a minimum sample of 418 Head-NCCT scans (167 Head-NCCT scans image with
radiologist-confirmed infarcts, 251 Head-NCCT scans images without infarcts,
2:3 ratio) would provide a minimum of 80% power to estimate an anticipated AUC
of 80% with 7% precision assuming a Type I error rate of 5% [20][21]. The
Head-NCCT scans, and their signed-off original radiological report performed
from 01-September-2021 to 31-August-2022 were acquired from diagnostic imaging
centers across India. A total of 1878 Head-NCCT scan were collected. The
original radiological report of these scans was subjected to a manual review
by a clinical data abstractor to classify the scans into infarct, and non-
infarct reported scans based on the original radiological report. A stratified
random sample of 500 Head-NCCT scans stratified by the presence and absence of
infarct (based on the original radiological reports) were then selected for
independent ground truthing by a radiologist with more than fourteen years of
experience. The inclusion criteria were Head-NCCT scans with soft
reconstruction kernel covering the complete brain, slice thickness $\leq 6mm$.
The exclusion criteria were Head-NCCT scans with obvious postoperative defects
or from patients who had previously undergone brain surgery, Head-NCCT scans
with artifacts such as burr holes, shunts or clips, Head-NCCT scans containing
metal artifacts, excessive motion artifacts, Head-NCCT scans containing
missing and improperly ordered slices. The ground truther radiologist had
access to the original head NCCT scan image but was blinded to the original
radiology report. The ground truther reviewed all the Head-NCCT scans and
provided segmentation boundaries for infarcts and intracranial hemorrhages.
The ground truther radiologist also provided a binary response for the
presence or absence of cranial fracture, midline shift, and mass effect. The
ground truth output was the reference standard for all downstream statistical
analyses, not the original radiological report.

The sensitivity and specificity were estimated based on a default device
threshold (available from the manufacturer based on internal testing), and the
optimum threshold was based on Youden’s index. The 95% confidence intervals
for sensitivity and specificity are reported based on exact method [22]. AUC
and 95% confidence interval (CI) was estimated based on the empirical method
and De Long methodology, respectively [23]. The segmentation provided by the
ground truther radiologist was utilized for the quantification analysis of the
error in the predicted infarct volume by the DL-based algorithm. Absolute
errors in infarct volume estimation in milliliter (mL), and summary statistics
of absolute errors were reported. The statistical analyses were performed
using RStudio (RStudio version 2022.07.1, R version 4.2.1) and Python version
3.9.7.

Table 1: Distribution of Head-NCCT Scan Images used for the Analysis
Stratified by Findings as Determined by Ground Truth.

Infarct (%) | Non-Infarct (%)
---|---
N (%) | Other subgroups | n # | N (%) | Other subgroups | n #
| No other target abnormality* | 170 |  | No other target abnormality* | 212
187 | ICH | 7 | 241 | ICH | 14
(43.7%) | Cranial Fracture | 5 | (56.3%) | Cranial Fracture | 18
| Midline Shift | 2 |  | Midline Shift | 2
| Mass Effect | 7 |  | Mass Effect | 8

  * *

Target abnormalities were custom defined and consisted of infarct,
intracranial hemorrhage, cranial fracture, midline shift and mass effect.

  * #

Denotes the absolute number of scans in the subgroups, and the numbers will
not add up to N because one scan may contain multiple target abnormalities.

##  3 Experimental Results

###  3.1 Identification of ICH and Infarct

The ground truthing was completed for 428, while 22 Head-NCCT scan were
excluded due to the inclusion and exclusion criteria mentioned in section 2. A
total of 187 Head-NCCT scan confirmed (based on ground truth) the presence,
while 241 Head-NCCT scan confirmed the absence of any infarcts. This
distribution of scans with and without infarcts met the minimum sample size
requirements described earlier in 2. In addition, 21 scans with intracranial
hemorrhages (ICH) and 23 scans with cranial fractures were present in the
sample. A total of 212 (49.5%) of the 428 Head-NCCT scans did not contain any
infarcts, intracranial hemorrhages, cranial fracture, midline shift, or mass
effect. The distribution of the Head-NCCT scans is shown in Table. 1.

It can be observed from Table. 2 that the DL-based algorithm achieved an AUC
of 86.8% (95% CI: 83.4 - 90.2) in detecting scans with the presence of
infarcts while the sensitivity and specificity were estimated to be 66.8% (95%
CI: 59.6-73.5)and 86.7% (95% CI: 81.8-90.7) respectively at the default
threshold. The optimum operating threshold was determined using Youden’s
index. At this optimum threshold, it was observed that the sensitivity of the
DL-based algorithm improved to 80.2% (95% CI: 73.8 - 85.7) without substantial
reduction in specificity 80.1% (95% CI: 74.5 - 84.9). For ICH, an AUC of 94.8%
(95% CI: 87.4 - 100) was achieved. There was no change in sensitivity compared
to the default and optimum threshold, while the specificity increased by 3%
using the optimum threshold. In contrast, the sensitivity of cranial fracture
compared to the default and optimum threshold, an enhancement of 15.8% was
observed while the specificity decreased by 2.7%. In Fig. 1, the AUC-ROC plot
for Cranial Fracture, ICH, and Infarct is given.

Table 2: Performance Evaluation of DL-based algorithm at Default and Optimum
Threshold; AUC: area under the curve; SN: sensitivity; SP: specificity, TP:
true positive, TN: true negative, AP1: all positive and AN1: all negative

| Target
---
| AUC (95% CI)
---
Default threshold | Optimum threshold
Abnormality |  |  | SN (95% CI)
---
(TP/AP)
| SP (95% CI)
---
(TN/AN)
| SN (95% CI)
---
(TP/AP)
| SP (95% CI)
---
(TN/AN)
Infarct |  | 86.8%
---
(83.4-90.2)
| 66.8%
---
(59.6-73.5)
(125/187)
| 86.7%
---
(81.8-90.7)
(209/241)
| 80.2%
---
(73.8-85.7)
(150/187)
| 80.1%
---
(74.5-84.9)
(193/241)
| Intracranial
---
hemorrhage
| 90.5%
---
(69.6 – 98.8)
| 90.5%
---
(69.6-98.8)
(19/21)
| 93.1%
---
(90.2 – 95.4)
(379/407)
| 90.5%
---
(69.6-98.8)
(19/21)
| 96.1%
---
(93.7 – 97.7)
(391/407)
| Cranial
---
fracture
| 95.2%
---
(90.1 – 100)
| 78.3%
---
(56.3 – 92.5)
(18/23)
| 96.8%
---
(94.6 – 98.3)
(392/405)
| 91.3%
---
(72.0 – 98.9)
(21/23)
| 94.1%
---
(91.3 – 96.2)
(381/405)

  * 1

$AP=TP+FN$, $AP=TN+FP$.

  * *

Default threshold: A threshold determined from internal testing results during
model training.

Figure 1: AUC-ROC plot: (i) Cranial Fracture (ii) ICH and (iii) Infarct.

###  3.2 Quantification of Infarct Volume

The DL-based algorithm for identifying infarcts produces the infarct volume in
mL. A total of 150 true positive scans for which both DL-based algorithms
predicted volume and ground truth volume were available for this analysis. The
reference standard was radiologist annotations done for each Head-NCCT scan
images.

Table 3: Performance Evaluation of DL-based algorithm in Quantifying Infarct; MSE: Mean absolute error; STD: standard deviation. Infarct Volume | NI |  | MAE
---
| Median
---
| STD
---
Percentile
|  |  |  |  |  | 10%
---
| 90%
---
0 - 5 mL | 108 | 3.200 | 0.158 | 1.121 | 0.032 | 2.746
$>$ 5 mL | 42 | 8.557 | 25.514 | 34.751 | 5.657 | 80.236
Overall | 150 | 4.700 | 0.600 | 24.397 | 0.042 | 42.080
Figure 2: 1) Scatter plot for infarct volume quantification with volume
distribution of 0-5 mL and $>$ 5 mL. 2) Bland-Altman plot for infarct volume
quantification

The mean absolute error (MAE) was 4.7 mL for overall scans. Based on ground
truth volume, the scans were further divided into two categories - scans with
0 - 5 mL and $>$ 5 mL infarcts volume, respectively. It can be observed from
Table. 3 that the MAE for 0 - 5 mL and $>$ 5 mL scans were found to be 3.2 mL
and 8.557 mL, respectively. In Fig. 2 from the scatter plot of infarct volumes
(1), it can be observed that with an increase in infarct volume, there is a
positive correlation between DL-based algorithm volume and ground-truth
annotated volume. The Bland-Altman plots showing good agreement between the
ground truther annotation and predicted volume by the DL-based algorithm are
shown in Fig. 2 (2).

###  3.3 Visual Explanations of DL-based Algorithm

The experimental findings depict that the evaluated DL-based algorithm
achieved superior performance represented in Table. 2 and 3. In most DL-based
models the rationale behind the prediction is not reveled explicitly. Since
these DL black box models can not be decomposed into intuitive and
comprehensive modules, these models are hard to interpret. Consequently, the
end-users develop skepticism and find the model difficult to trust. The
emergence of explainable artificial intelligence (XAI) is an essential aspect
of model transparency and the social right to explain DL inferences [24],[25].
XAI encompasses a better understanding of incoherent output, isolates failure
modes, and builds trust in intelligent systems for effective incorporation
into our everyday lives [26]. The present evaluated DL-based algorithm outputs
a boundary across the infarcts which revels rationale behind the superior
performance. In Fig. 3 is can be observed that for both small and large
infarcts volume on Head-NCCT scan, the model predicted boundary clearly
overlaps with the ground truther boundary.

Figure 3: Representative two head NCCT scan images with infarct. (a) The head
NCCT scan with presence of infarct (b) The ground truther annotation around
the region of infarct with red boundary (c) DL-based algorithm volume
quantification (blue boundary) overlapping with ground truther annotation (red
boundary).

##  4 Discussion

This retrospective study evaluated a deep learning algorithm for detecting
infarcts in Head-NCCT scans. The algorithm had a good AUC of about 86% in
detecting infarcts. After adjusting for thresholds, a balanced sensitivity of
80.2% and specificity of 80.1% was estimated to detect infarcts. The
algorithm’s sensitivity in detecting infarcts in scans with no other target
abnormalities was found to be 80% (136 correctly detected out of 170). It did
not differ from the overall sensitivity at optimum sensitivity. This states
the robustness of the DL-based algorithm to identify infarcts with negligible
drop in sensitivity with presence of other abnormalities. Additionally, it is
to be noted that the sensitivity of Head-NCCT scans in detecting infarcts is
generally considered low, especially in the case of hyperacute and acute
ischemic strokes. In one study, the sensitivity of detecting acute ischemic
stroke on head NCCT scans ranged from 57% to 71% with considerable inter-
reader variability [27][28]. Additionally, we evaluated the performance to
detect ICH and cranial fracture, and both had excellent AUC. However, the
interpretation is limited by low sample sizes for these two abnormalities. Our
results also show that threshold adjustments might be needed before using such
algorithms routinely for clinical decision support.

Deep learning or big data are often called ”black box” and represent
substantial obstacles in introducing intuitive and comprehensive modules into
actual clinical practice; these models are challenging to interpret. However,
the DL-based method validated in this study provides a post-hoc attention tool
for the clinician to identify the lesion visually. In addition, the DL-based
algorithm validated in this study encompasses a better understanding of
incoherent output, isolates failure modes, and builds trust in intelligent
systems for effective incorporation into routine clinical practice. Moreover,
the proposed validation of the DL-based algorithm will be beneficial in the
resource constraint areas with a limited number of radiologists or with only
access to teleradiology facilities.

Our study has limitations. First, the differentiation of infarct into acute
and chronic infarct was not analyzed. Second, the ground truthing for the head
NCCT scans images with the presence of infarcts was done by a single
radiologist. Thirdly, there were not enough scans for the ICH and cranial
fracture to estimate performance metrics with sufficient precision.

##  5 Conclusion

The present study evaluated a DL-based algorithm to determine the presence and
absence of ICH and infarcts on head-NCCT scans. The DL-based algorithm
demonstrated high detection performance rate in identifying infarcts, ICH, and
cranial fracture. Additionally, the DL-based algorithm exhibits a positive
correlation between DL-based algorithm volume and ground-truth annotated
volume. The study demonstrated the performance of ICH detection and infarcts
detection and quantification to indicate the feasibility of introduction of
such DL-algorithms in routine workflow in extensive healthcare facilities.

##  6 Data Availability

The datasets used or analyzed during the current study are available from the
corresponding author on reasonable request.

## References

  *   * Fellows et al. [2005] Fellows, L.K., Heberlein, A.S., Morales, D.A., Shivde, G., Waller, S., Wu, D.H.: Method matters: An empirical study of impact in cognitive neuroscience. Journal of Cognitive Neuroscience 17, 850–858 (2005)
  * Geva et al. [2011] Geva, S., Jones, P.S., Crinion, J.T., Price, C.J., Baron, J.C., Warburton, E.A.: The neural correlates of inner speech defined by voxel-based lesion–symptom mapping. Brain 134, 3071–3082 (2011)
  * Utz et al. [2012] Utz, S., Humphreys, G.W., Chechlacz, M.: Parietal substrates for dimensional effects in visual search: evidence from lesion-symptom mapping. Brain : a journal of neurology 136 Pt 3, 751–60 (2012)
  * Haji and Naval [2020] Haji, S., Naval, N.S.: Management of intracerebral hemorrhage. Evidence-Based Critical Care (2020)
  * Rymer [2011] Rymer, M.M.: Hemorrhagic stroke: intracerebral hemorrhage. Missouri medicine 108 1, 50–4 (2011)
  * Phipps and Cronin [2020] Phipps, M.S., Cronin, C.: Management of acute ischemic stroke. BMJ 368 (2020)
  * Rorden et al. [2012] Rorden, C., Bonilha, L., Fridriksson, J., Bender, B., Karnath, H.-O.: Age-specific ct and mri templates for spatial normalization. NeuroImage 61, 957–965 (2012)
  * Fiez et al. [2000] Fiez, J.A., Damasio, H., Grabowski, T.J.: Lesion segmentation and manual warping to a reference brain: Intra‐ and interobserver reliability. Human Brain Mapping 9 (2000)
  * Ashton et al. [2003] Ashton, E.A., Takahashi, C., Berg, M.J., Goodman, A., Totterman, S.M.S., Ekholm, S.: Accuracy and reproducibility of manual and semiautomated quantification of ms lesions by mri. Journal of Magnetic Resonance Imaging 17 (2003)
  * Pratiher and Chattoraj [2018] Pratiher, S., Chattoraj, S.: Diving deep onto discriminative ensemble of histological hashing & class-specific manifold learning for multi-class breast carcinoma taxonomy. ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 1025–1029 (2018)
  * Nawn et al. [2020] Nawn, D., Pratiher, S., Chattoraj, S., Chakraborty, D., Pal, M., Paul, R.R., Dutta, S., Chatterjee, J.: Multifractal alterations in oral sub-epithelial connective tissue during progression of pre-cancer and cancer. IEEE Journal of Biomedical and Health Informatics 25, 152–162 (2020)
  * Govindarajan et al. [2022] Govindarajan, A., Govindarajan, A., Tanamala, S., Chattoraj, S., Reddy, B., Agrawal, R., Iyer, D., Srivastava, A., Kumar, P., Putha, P.: Role of an automated deep learning algorithm for reliable screening of abnormality in chest radiographs: A prospective multicenter quality improvement study. Diagnostics 12 (2022)
  * Pratiher et al. [2018] Pratiher, S., Chattoraj, S., Mukherjee, R.: Stationplot: A new non-stationarity quantification tool for detection of epileptic seizures. 2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP), 499–503 (2018)
  * Chilamkurthy et al. [2018] Chilamkurthy, S., Ghosh, R., Tanamala, S., Biviji, M., Campeau, N.G., Venugopal, V.K., Mahajan, V., Rao, P., Warier, P.: Deep learning algorithms for detection of critical findings in head ct scans: a retrospective study. The Lancet 392, 2388–2396 (2018)
  * Wilke et al. [2011] Wilke, M., Haan, B., Juenger, H., Karnath, H.-O.: Manual, semi-automated, and automated delineation of chronic brain lesions: A comparison of methods. NeuroImage 56, 2038–2046 (2011)
  * Rekik et al. [2012] Rekik, I., Allassonnière, S., Carpenter, T., Wardlaw, J.M.: Medical image analysis methods in mr/ct-imaged acute-subacute ischemic stroke lesion: Segmentation, prediction and insights into dynamic evolution simulation models. a critical appraisal. NeuroImage : Clinical 1, 164–178 (2012)
  * Chan [2007] Chan, T.: Computer aided detection of small acute intracranial hemorrhage on computer tomography of brain. Computerized medical imaging and graphics : the official journal of the Computerized Medical Imaging Society 31 4-5, 285–98 (2007)
  * Liu et al. [2008] Liu, R., Tan, C.L., Leong, T.-Y., Lee, C.K., Pang, B.C., Lim, C.C.T., Tian, Q., Tang, S., Zhang, Z.: Hemorrhage slices detection in brain ct images. 2008 19th International Conference on Pattern Recognition, 1–4 (2008)
  * Matesin et al. [2001] Matesin, M., Lonvaric, S., Petravic, D.: A rule-based approach to stroke lesion analysis from ct brain images. ISPA 2001. Proceedings of the 2nd International Symposium on Image and Signal Processing and Analysis. In conjunction with 23rd International Conference on Information Technology Interfaces (IEEE Cat., 219–223 (2001)
  * Obuchowski [1994] Obuchowski, N.A.: Computing sample size for receiver operating characteristic studies. Investigative Radiology 29, 238–243 (1994)
  * Zhou et al. [2002] Zhou, X.-H., Obuchowski, N.A., McClish, D.K.: Statistical methods in diagnostic medicine. (2002)
  * Clopper and Pearson [1934] Clopper, C.J., Pearson, E.S.: The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika 26, 404–413 (1934)
  * DeLong et al. [1988] DeLong, E.R., DeLong, D.M., Clarke‐Pearson, D.L.: Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach. Biometrics 44 3, 837–45 (1988)
  * Došilović et al. [2018] Došilović, F.K., Brčić, M., Hlupić, N.: Explainable artificial intelligence: A survey. In: 2018 41st International Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO), pp. 0210–0215 (2018). IEEE
  * Adadi and Berrada [2018] Adadi, A., Berrada, M.: Peeking inside the black-box: A survey on explainable artificial intelligence (xai). IEEE Access 6, 52138–52160 (2018)
  * Gunning [2017] Gunning, D.: Explainable artificial intelligence (xai). Defense Advanced Research Projects Agency (DARPA), nd Web 2(2) (2017)
  * Broderick et al. [2007] Broderick, J.P., Connolly, S.E., Feldmann, E., Hanley, D.F., Kase, C.S., Krieger, D., Mayberg, M.R., Morgenstern, L.B., Ogilvy, C.S., Vespa, P.M., Zuccarello, M.: Guidelines for the management of spontaneous intracerebral hemorrhage in adults. 2007 update. guideline from the american heart association/american stroke association stroke council, high blood pressure research council, and the quality of care and outcomes in research interdisciplinary working gro. Stroke 38, 2001–2023 (2007)
  * Srinivasan et al. [2006] Srinivasan, A., Goyal, M., Azri, F.A., Lum, C.: State-of-the-art imaging of acute stroke. Radiographics : a review publication of the Radiological Society of North America, Inc 26 Suppl 1, 75–95 (2006)


# Multitemporal SAR images change detection and visualization using RABASAR
and simplified GLR

Weiying Zhao, Henri Maître, Jean-Marie Nicolas and Florence Tupin
LTCI, Télécom ParisTech
Université Paris-Saclay
75013 Paris, France
&Charles-Alban Deledalle
MB, CNRS, Bordeaux INP
University of Bordeaux
F-33405 Talence, France
&Loïc Denis
Laboratoire Hubert Curien UMR 5516, UJM-Saint-Etienne, CNRS
Institut d Optique Graduate School, University of Lyon
F-42023 SAINT-ETIENNE, France

###### Abstract

Understanding the state of changed areas requires that precise information be
given about the changes. Thus, detecting different kinds of changes is
important for land surface monitoring. SAR sensors are ideal to fulfil this
task, because of their all-time and all-weather capabilities, with good
accuracy of the acquisition geometry and without effects of atmospheric
constituents for amplitude data. In this study, we propose a simplified
generalized likelihood ratio ($S_{GLR}$) method assuming that corresponding
temporal pixels have the same equivalent number of looks (ENL). Thanks to the
denoised data provided by a ratio-based multitemporal SAR image denoising
method (RABASAR), we successfully applied this similarity test approach to
compute the change areas. A new change magnitude index method and an improved
spectral clustering-based change classification method are also developed. In
addition, we apply the simplified generalized likelihood ratio to detect the
maximum change magnitude time, and the change starting and ending times. Then,
we propose to use an adaptation of the REACTIV method to visualize the
detection results vividly. The effectiveness of the proposed methods is
demonstrated through the processing of simulated and SAR images, and the
comparison with classical techniques. In particular, numerical experiments
proved that the developed method has good performances in detecting farmland
area changes, building area changes, harbour area changes and flooding area
changes.

_K_ eywords Multitemporal $\cdot$ synthetic aperture radar $\cdot$ change
detection $\cdot$ clustering $\cdot$ visualization

##  1 Introduction

Timely, accurate and continuous monitoring of land cover and land use changes
is important for land resource management. According to [1], land cover
changes can be classified into several categories: changed from one land cover
class to another, change of shape, shrink or transform, change of position,
fragment or merge of adjacent regions. Based on the change reason or change
type, changes can be classified into: short-term change (synoptic weather
events), cyclic change (seasonal phenology), directional change (urban
development), multidirectional change (deforestation & regeneration), event
change (catastrophic fires). Changes happen spatially and temporally. More and
more available multispectral, multitemporal, multisensor satellite data
enhance the capability of detecting, identifying, mapping and monitoring these
changes [2, 3].

Optical remote sensing images with high spatial and spectral resolutions can
be easily acquired and have been widely used for land cover monitoring [4].
The passive acquisition model of optical sensors and the use of near-visible
light wavelengths, sun illumination (or thermal radiation) and cloud-free
weather requirements heavily limit its wide application.

However, SAR imagery has the main advantage of being an all-time and all-
weather sensor, with good accuracy for the acquisition geometry and with few
atmospheric effects for amplitude data. It has been widely used for
environmental change detection, urban area change detection and disaster
monitoring. The past few years have seen the launch of numerous synthetic
aperture radar (SAR) sensors and new sensors will also be launched soon, such
as TSX-NG, Cosmo-SkyMed second generation, Radarsat constellation and NISAR.
The increasing availability of SAR data allows the high accuracy of change
detection, such as abrupt (step) changes, seasonal changes and longer-term
developments. Unlike optical images, SAR images are seriously affected by
speckle noise, which makes the SAR image change detection much harder. SAR
image change detection can be applied on image pairs or multitemporal series.
Different SAR image characteristics can be used for change detection. The
current change detection methods are mainly based on likelihood ratio [5, 6,
7, 8], coherence [9], image texture and structure analysis [10], and deep
learning based methods [11, 12]. Due to the multiplicative noise in coherent
SAR images, the likelihood ratio test is popularly used for change detection.
A number of studies have applied to SAR data change detection methods based on
mean-ratio [13], image ratio [14, 15], log-ratio operator [16, 17]. The true
distribution of these ratio images depends on the relative change of the SAR
reflectivities [18]. These methods are easily applied and the associated
threshold can be calculated automatically.

Visualizing the changes is also important. Changes usually represent
transitions that occur between states [19]. Using a colourful image can ease
results interpretation highlighting changed areas. In the change detection
field, a number of studies have used different colours to show different
changes. Su et al. [6] associate different colours to highlight different
change types of time series, Mou et al. [20] propose using different colours
to represent different change phenomenons, and Domínguez et al. [21] use
different fusion strategies to illustrate the real changes and false alarms.
With short time series, Nielsen et al. [22] propose to use RGB colours to
represent different change times and use black colour to represent the
unchanged areas. In addition, Amitrano et al. [23] proposed using new bi-
temporal and multitemporal RGB combination frameworks to illustrate temporal
SAR images. The effectiveness of these two RGB visualization approaches has
been verified by change detection and classification, respectively. However,
none of them associates the colours with the times of change in a long time
series. Recently, Rapid and EAsy Change detection in radar TIme-series by
Variation coefficient (REACTIV) method has been proposed by Koeniguer et al.
[24]. It is a simple and highly efficient time series change detection and
visualization algorithm. It is based on HSV space and exploits only time
domain estimates without any spatial estimation. The colour saturation is
coded by the temporal coefficient of variation. However, the detection results
are corrupted by speckle noise. Even using some state-of-the-art denoising
methods, the bias estimation in vegetation areas still prohibits REACTIV to
provide accurate performance. In addition, the colour in REACTIV results only
represents the appearing date of maximum values.

The inherent speckle which is attached to any coherent imaging system affects
the analysis and interpretation of synthetic aperture radar (SAR) images.
Thanks to the recently proposed Ratio-Based Multitemporal SAR Images Denoising
(RABASAR) method [25], we can significantly suppress the negative effect of
SAR speckle. In this paper, we derive a simplified generalized log-likelihood
ratio criterion ($S_{GLR}$) based on gamma distribution using RABASAR
denoising data. Based on the obtained similarity function, strategies to apply
this similarity criterion to image pair change detection, cumulative change
detection and change classification are introduced. In addition, we adapt
REACTIV method to integrate RABASAR denoising results. Then, we apply the
simplified GLR function in this framework to detect different change times of
interest, such as change starting and ending times under predefined threshold,
or maximum change magnitude appearing time.

##  2 Image pair change area detection and change magnitude visualization

In this section, we derive a generalized log-likelihood ratio criterion based
on gamma distribution using denoised data. Based on the obtained similarity
functions, strategies to apply this similarity criterion to image pair change
detection, cumulative change detection, change classification and change time
detection are introduced.

###  2.1 Change area detection

Under Goodman’s hypothesis [26], the fully developed intensity speckle follows
a Gamma distribution $\mathcal{G}[u,L]$ depending on the number of looks $L$
and the mean reflectivity $u$ of the scene:

$\mathcal{G}[u,L](y)=\frac{L}{u\Gamma(L)}\bigg{(}\frac{Ly}{u}\bigg{)}^{L-1}e^{-\frac{Ly}{u}}$
(1)

To compare the similarity of two gamma distributed variables
$(y_{t},y_{t^{\prime}})$, a log-likelihood ratio test can be used. Speckle is
an inherent problem for SAR image interpretation which brings many drawbacks
for traditional SAR image change detection methods. When dealing with denoised
images $\hat{u}_{t}$ and $\hat{u}_{t^{\prime}}$, with associated ENL $L_{t}$
and $L_{t^{\prime}}$, we have:

$S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})=L_{t}\log\frac{L_{t}\hat{u}_{t}+L_{t^{\prime}}\hat{u}_{t^{\prime}}}{\hat{u}_{t}(L_{t}+L_{t^{\prime}})}+L_{t^{\prime}}\log\frac{L_{t}\hat{u}_{t}+L_{t^{\prime}}\hat{u}_{t^{\prime}}}{\hat{u}_{t^{\prime}}(L_{t}+L_{t^{\prime}})}$
(2)

where $\hat{u}$ is the estimated reflectivity, $t$ represents the time index
in the time series.

Unlike $CGLRT$ and $AGLRT$ methods [6], we fully trust the denoising results
and do not take the noisy data into account any more. In practice, we directly
calculate the similarity of the multi-looked SAR data with the use of their
corresponding ENL. Suppose the corresponding pixels in despeckling images have
the same ENL $L_{t}=L_{t^{\prime}}=\hat{L}$, the simplified GLR method [27,
28] turns out to be:

$S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})=2\hat{L}\log{\bigg{(}\sqrt{\frac{\hat{u}_{t}}{\hat{u}_{t^{\prime}}}}+\sqrt{\frac{\hat{u}_{t^{\prime}}}{\hat{u}_{t}}}}\bigg{)}-2\hat{L}\log{2}$
(3)

Defining a global threshold is a simple and widely used approach to
distinguish changes from unchanged points. To detect the changed areas, we
used a thresholding function:

$\varphi[S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})]=\left\\{\begin{array}[]{cc}1,&\mathrm{if\
}S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})\geq\tau\\\ 0,&\mathrm{otherwise}\\\
\end{array}\right.$ (4)

The threshold definition methods are under the same no-change hypothesis
framework that can be used for this method. As introduced by [29], when the
sample size approaches infinity, the log-likelihood statistic model
asymptotically converges towards chi-squared distributed probability under the
null hypothesis. Thus, we can use the chi-square cumulative function to
estimate the change probability of $S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})$
with:

$\rho=1-\frac{1}{4\hat{L}}$ (5)
$\omega_{2}=-\frac{1}{4}(1-\frac{1}{\rho})^{2}$ (6) $P\\{2\rho
S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})\leq\delta\\}\simeq
P\\{\chi^{2}(1)\leq\delta\\}+\omega_{2}[P\\{\chi^{2}(5)\leq\delta\\}-P\\{\chi^{2}(1)\leq\delta\\}]$
(7)

where $\delta$ is the statistical significance. For the detailed derivation of
this probability calculation method, we recommend referring to [7, 30]. Unlike
Conradsen’s PolSAR change detection analysis, we mainly pay attention to
spatially adaptive denoised single-channel SAR images with the same estimated
ENL. To robustly estimate $\hat{L}$ in the denoised image, the log-cumulant
method is used [31]. In the following sections, we mainly use this way to
define the threshold.

###  2.2  Change magnitude index for visualization

To distinguish appearing from disappearing changes, we used a signum function
$\texttt{sign}(x)$ to convert $S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})$ to
positive and negative values. In this case, if we set the image acquired at
time $t$ as the reference image, the positive and negative values correspond
to the increase and decrease of the object backscattering values.

$\texttt{sign}(x)\begin{cases}-1&\texttt{if }x<0\\\ 0&\texttt{if }x=0\\\
1&\texttt{if }x>0\end{cases}$ (8)

$x=\log\bigg{(}\sqrt{\frac{\hat{u}_{t^{\prime}}}{\hat{u}_{t}}}\bigg{)}$ (9)

To clearly illustrate the temporal changes, the similarity ratio is normalized
and transformed to values within the range [0, 255].

$S^{conv}_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})=\left\\{\begin{array}[]{cc}255,&\mathrm{if\
}\frac{2(S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})-\alpha_{1})}{\alpha_{2}-\alpha_{1}}\geq
2\\\
127\frac{2(S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})-\alpha_{1})}{\alpha_{2}-\alpha_{1}}+1,&\mathrm{otherwise}\\\
\end{array}\right.$ (10)

where $\alpha_{1}$ and $\alpha_{2}$ represent the minimum and maximum values
in the temporal dissimilarities $S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})$.
To suppress the outliers, we empirically set $\alpha_{1}=-2$ and
$\alpha_{2}=2$. The value range will be converted to [-255, 255] by
multiplying $S^{conv}_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})$ with
$\texttt{sign}(x)$. The rainbow index colour is used to represent different
change magnitudes (appearing, disappearing and slow changes).

In practice, we can arbitrarily combine the reference and slave image through
RGB composition (R: slave image, G: reference image, B: slave image). Although
different colours could indicate the increase and decrease of backscattering
values, the illustration performance is not as good as the aforementioned
strategy.

###  2.3 Time series change type classification

During the time series acquisition, changes may occur multiple times and with
different magnitudes. To detect the change types, we propose an improved
change classification method inspired by NORCAMA method [6] and spectral
clustering method [32, 33]. In practice, the change types are transferred into
a partitioning problem and detected using spectral clustering.

####  2.3.1 Change Criterion Matrix (CCM)

By making use of the eigenvalues of the similarity matrix of the data,
spectral clustering techniques perform dimensionality reduction before
clustering in fewer dimensions. It has been successfully used to cluster the
temporal pixels based on their similarity symmetric matrix [33]. Given a time
series $\\{\hat{u}_{1},\hat{u}_{2}\cdots\hat{u}_{M}\\}$, the affinity matrix
is defined as a symmetric matrix $A(s)$, with elements
$S(\hat{u}_{t},\hat{u}_{t^{\prime}})$ representing the change criterion
between different data points.

$A(s)=\begin{bmatrix}S(\hat{u}_{1},\hat{u}_{1})&S(\hat{u}_{1},\hat{u}_{2})&S(\hat{u}_{1},\hat{u}_{3})&\dots&S(\hat{u}_{1},\hat{u}_{t^{\prime}})\\\
S(\hat{u}_{2},\hat{u}_{1})&S(\hat{u}_{2},\hat{u}_{2})&S(\hat{u}_{2},\hat{u}_{3})&\dots&S(\hat{u}_{2},\hat{u}_{t^{\prime}})\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
S(\hat{u}_{t},\hat{u}_{1})&S(\hat{u}_{t},\hat{u}_{2})&S(\hat{u}_{t},\hat{u}_{3})&\dots&S(\hat{u}_{M},\hat{u}_{M})\end{bmatrix}$
(11)

where $A(s)$ is a symmetric change criterion matrix with size $M\times M$, $s$
is the location in one image, $t$ and ${t^{\prime}}$ are time index with
$1\leq t\leq M$ and $1\leq{t^{\prime}}\leq M$.

To avoid the overlapping of different clusters, Xin $et$ $al.$ [6] proposed to
binarize the change criterion matrix $A(s)$. The binary process can tighten
the clusters, but it will force the clustering results seriously depending on
the used thresholds. In practice, we use a binarized change criterion matrix.
In addition, the $k$-nearest neighbours algorithm could be used to classify
this change criterion matrix as well. To suppress the temporal variance caused
by the residual speckle, we can apply the exponentially weighted moving
average [34] to the time series.

####  2.3.2 Clustering by spectral clustering method

Based on the acquired change criterion matrix, the Laplacian matrix $A^{L}(s)$
is computed by:

$A^{L}(s)=D(s)-A(s)$ (12)

$D(s)=\begin{bmatrix}\sum S(\hat{u}_{1},\hat{u}_{t^{\prime}})&0&\dots&0\\\
0&\sum S(\hat{u}_{2},\hat{u}_{t^{\prime}})&\dots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\dots&\sum
S(\hat{u}_{t},\hat{u}_{t^{\prime}})\end{bmatrix}$ (13)

$\sum
S(\hat{u}_{t},\hat{u}_{t})=\sum_{{t^{\prime}}=1}^{M}S(\hat{u}_{t},\hat{u}_{t^{\prime}})$
(14)

The former steps (Eq.(11)$\sim$(14)) are the same as NORCAMA [6] except the
similarity calculation. Then, the Laplacian matrix $A^{L}(s)$ is normalized
using:

$A^{L}_{norm}={D}^{-1/2}A^{L}{D}^{-1/2}$ (15)

The eigenvalues $\lambda$ are computed through:

$A^{L}_{norm}V=\lambda V$ (16)

where $V$ is the eigenvector. After sorting the acquired eigenvalues
$\\{\lambda_{1},\lambda_{2},\cdots,\lambda_{M}\\}$ in ascending order, the
clustering number $k$ is calculated using the eigengap heuristic method [35]:

$k=\operatorname*{arg\,max}_{1\leq t<M}(\lambda_{t+1}-\lambda_{t})$ (17)

To reduce the data dimension, only the eigenvectors $\upsilon_{t}$ ( $M\times
1$ column vector) which correspond to the $k$ largest eigenvalues of
$L^{norm}$ are used with $U=[\upsilon_{1},\upsilon_{2}\cdots\upsilon_{k}]$. To
obtain the unit norm, we re-normalize the matrix rows. Finally, the $k$-means
method is used to cluster each row $u_{t}$ in $U$ and the cluster labels
$\\{l_{1},l_{2},...l_{M}\\}$ are assigned to each cluster element with $1\leq
l_{t}\leq k$.

The aforementioned method is similar to the normalized cut method [32, 6].
However, they normalize the rows of $A(s)$ to sum to 1 and use its
eigenvectors instead of the normalized Laplacian matrix $A^{L}(s)$ calculated
using equation (15). In addition, they do not re-normalize the rows of $U$ to
unit length [33].

####  2.3.3 Change type recognition

Based on the number of clusters $k$ and cluster labels
{$l_{1},l_{2},...l_{M}$} acquired by the $k$-means algorithm, the change type
of the time series points can be recognized [6] according to Table 1.

Table 1: Label of different change types Classes | Types | $k$ | Label series {$l_{1},l_{2},...l_{M}$}
---|---|---|---
1 | Unchanged | 1 | 1, 1, …1
2 | Step | 2 | 1, 1, …1, 2, 2, …2
3 | Impulse | 2 | 1, 1, …1, 2, 2, …2, 1, 1, …1
4 | Cycle | 2 | 1, …1, 2, …2, 1, …1, 2,…2,…
5 | Complex |  $\geq$3 | 1, 1, …, 2, 2…, 3, 3…4, 4…

##  3 Change times of interest detection and visualization with extended
REACTIV

In this section, our objective is to adapt the REACTIV method to integrate
RABASAR denoising results. Then, using the proposed simplified GLR function to
detect the change time of interest.

The principle of REACTIV is to exploit the HSV colour space and a temporal
stack of SAR images. The hue channel H represents the time, the saturation
channel S corresponds to the temporal coefficient of variation, and the value
V corresponds to the maximum radar intensity of the temporal series in each
pixel [36].

###  3.1 Times of interest (Hue)

As introduced in [36], one can associate a colour with a particular time
according to the change. During the procedure, different change types can be
considered, such as abrupt change, seasonal change, deforestation &
regeneration, etc. The REACTIV visualization method chooses to highlight the
appearing time of the maximum value. Although the REACTIV visualization method
can well associate the maximum value appearing time with the colour, the first
and last dates have very similar colours because the HSV colour palette is
continuous and loops on itself. Thus, we propose to highlight the interested
time using the normalized time in the interval:

$f_{t}=\frac{5}{6}\times\frac{t-t_{1}}{t_{2}-t_{1}}$ (18)

where $t_{1}$ and $t_{2}$ are the first and last image acquisition time in the
time series, and $t$ is the time of interest. ${5}/{6}$ is used to suppress
the time interval, so as to avoid using the starting color category (loop)
again.

With the time series
$\\{\hat{u}_{1}(s),\hat{u}_{2}(s),\cdots,\hat{u}_{M}(s)\\}$, we can use
$S_{GLR}$ to detect times of interest:

  * •

Start changing time

When detecting the start changing time in the time series, points similarities
$S_{GLR}(\hat{u}_{1}(s),\hat{u}_{t}(s))$ are calculated with the reference to
the value of the first date $\hat{u}_{1}(s)$. After transforming the
similarity to change probability $P\\{2\rho
S_{GLR}(\hat{u}_{1},\hat{u}_{t^{\prime}})\\}$, we can decide whether there is
a change or not based on a predefined change probability $\tau$ (such as
$99\%$).

$T_{start}=\begin{cases}t&\texttt{if }P\\{2\rho
S_{GLR}(\hat{u}_{1},\hat{u}_{t^{\prime}})\\}>\tau\\\
0&\texttt{else}\end{cases}$ (19)

where $t$ is the start changing time with $1<t\leq M$. It corresponds to the
first time that the change probability is larger than the threshold $\tau$.

  * •

Maximum changing time

Generally, the abrupt changes are associated with large change magnitude [37].
The $S_{GLR}(\hat{u}_{t}(s),\hat{u}_{t^{\prime}}(s))$ is supposed to be the
maximum change. In this case, $t$ and $t^{\prime}$ are set as the adjacent
times. $t^{\prime}$ is the maximum change time.

$T_{max}=\begin{cases}t^{\prime}&\text{if
}S_{GLR}(\hat{u}_{t}(s),\hat{u}_{t^{\prime}}(s))=\tau_{max}\\\
0&\text{else}\end{cases}$ (20)

where $\tau_{max}$ equals to the maximum dissimilarity
$\tau_{max}=\operatorname*{arg\,max}\limits_{t^{\prime}}S_{GLR}(\hat{u}_{t}(s),\hat{u}_{t^{\prime}}(s))$.

  * •

Stop changing time

For the detection of stop changing time, the last changed point in the time
series is used as the reference point. In practice, the calculation is carried
out in reverse order.

$T_{stop}=\begin{cases}t&\text{if
}P\\{2\rho(S_{GLR}(\hat{u}_{t}(s),\hat{u}_{M}(s)))\\}>\tau\\\
0&\text{else}\end{cases}$ (21)

where $t$ corresponds to the dates going from $M,M-1,\cdots$ to 1.

###  3.2 Saturation (S)

This colour component defines whether there are changes or not in the time
series. Unlike popularly used change detection methods which mainly pay
attention to the intensity value changes, the REACTIV method uses the dynamics
of the coefficient of variation. Based on Rayleigh Nakagami’s distribution
[38], it is possible to derive the empirical moments’ expression of pure
speckle with [39]:

$m_{1}=u_{A}\frac{\Gamma(L+\frac{1}{2})}{\sqrt{L}\Gamma(L)}$ (22)
$m_{2}=u_{A}^{2}$ (23)

Based on the ratio of standard deviation and the amplitude average, the
coefficient of variation can be calculated through:

$\gamma=\frac{\sigma}{u_{A}}=\frac{\sqrt{m_{2}-m_{1}^{2}}}{m_{1}}=\sqrt{\frac{\Gamma(L)\Gamma(L+1)}{\Gamma(L+1/2)^{2}}-1}$
(24)

To know the behaviour of this parameter, the variance of this estimator can be
calculated according to [39, 40]:

$\mathrm{Var}(\gamma)=\frac{1}{4M}\frac{4m_{2}^{3}-m_{2}^{2}m_{1}^{2}+m_{1}^{2}m_{4}-4m_{1}m_{2}m_{3}}{m_{1}^{4}(m_{2}-m_{1}^{2})}$
(25)

where $M$ is the number of samples to compute the estimation which corresponds
here to the temporal images considered. In addition, based on the order 1 to 4
moments of the Nakagami distribution and $L$, the variance of this estimator
can be written:

$\mathrm{Var}(\gamma)=\frac{1}{4M}\frac{L\Gamma(L)^{4}(4L^{2}\Gamma(L)^{2}-4L\Gamma(L+\frac{1}{2})^{2}-\Gamma(L+\frac{1}{2})^{2})}{\Gamma(L+\frac{1}{2})^{4}(L\Gamma(L)^{2}-\Gamma(L+\frac{1}{2})^{2})}$
(26)

Usually, the number of available images on the same area varies. In order to
overcome the dependency of the coefficient of variation on the number of
images, it is possible to normalize the distribution with the theoretical mean
and standard deviation. Koeniguer et al. [36] propose to use the following
empirical normalization:

$\gamma\leftarrow\frac{\gamma-\mathbb{E}[\gamma]}{10\sigma(\gamma)}+0.25$ (27)

with the theoretical mean and standard deviation values for a "stable" speckle
and $L=4.9$ for Sentinel-1 GRD data. This empirical normalization aims to
reduce the saturation values of the stable zones around a low saturation to
0.25 and to spread the changes on higher saturation.

To speed up the processing, we will use a coefficient of variation method to
detect time series changes. Then, detect the different change times of
interest.

###  3.3 Value (V)

The REACTIV visualization method uses the maximum amplitude value of each time
series as the value channel. The hue colour component is computed using
equation (18) with time $t$ corresponding to the maximum amplitude value
appearing time. This choice is particularly suitable for an abrupt event (such
as the presence of a boat).

Although using maximum values can highlight abruptly appearing objects, the
acquired results seem too noisy. Apart from this choice, we can use the
denoising results or the temporal average image to have a clear vision of the
ground. We only use the maximum time series values in the following
experiments.

##  4 Experimental results and discussion

To illustrate and compare the proposed methods with state-of-the-art change
detection methods, simulated SAR images and SAR images are tested in this
section. All the data are despeckled by RABASAR before change analysis: image
pair change detection in section 4.2, continuous change monitoring in section
4.3, change classification experiments in section 4.4 and change time
detection in section 4.5.

###  4.1 Experimental data introduction

To test the improved algorithms, we prepared four kinds of data. In this
paper, we mainly consider SAR images acquired through the same orbit with
similar incidence angles. During the RABASAR denoising, we mainly use
arithmetic mean (AM), denoised arithmetic mean (DAM) and denoised binary
weighted arithmetic mean (DBWAM) of the time series.

  * •

Simulated SAR data

Many SAR image simulations are based on reflectivity maps obtained from
optical images. However, real SAR images exhibit strong and persistent
scatterers, especially in urban areas which can hardly be simulated using
optical images. Therefore, we propose to use the arithmetic mean image of long
time series of SAR images, considered as a noise-free image (a reflectivity
map $u$) to create realistic simulations of SAR images. With the acquired
arithmetic mean image $y^{AM}$, we add different kinds of object changes
according to their values in real SAR images. Then, we simulate the temporal
images by multiplying $y^{AM}$ with different simulated Gamma distribution
noise $v_{t}$.

  * •

TerreSAR-X data

The TerraSAR-X images are acquired over the harbour area of $Sendai$. 9
temporal well-registered SAR images are used for the preparation of denoised
data, both for 2SPPB [28] and for RABASAR. Only the two images which were
acquired on 06/05/2011 and 08/06/2011 are used for the change detection.

  * •

Sentinel-1 single look SAR data

The used Sentinel-1 IW VV single look SAR data are acquired over the Saclay
area, South of Paris. All the images are registered using geometric-based
subpixel images registration method [41]. This area mainly contains farmland
area, forest area, building area, etc. Saclay area has been chosen since it
has been chosen to receive the future scientific area of the University of
Paris-Saclay. Starting in 2010, constructions and public works were decided to
convert agricultural terrains into research and education buildings, mostly 2
to 5 storey compact and geometrical structures made of concrete, steel and
glass. Many plots have been barred from vegetation, and the excavated heavy
plant machinery and trucks have been parked in some places. All these elements
greatly influence the SAR reflectivity.

  * •

Sentinel-1 GRD data

The Sentinel-1 GRD data are acquired over Saddle Dam D, Southern Laos. All the
images are coregistered using geometric-based registration method111SNAP:
http://step.esa.int/main/toolboxes/snap/. All the historical Sentinel-1 CRD VV
descending images acquired over this area are used to compute the temporal
mean image, which is used by the RABASAR method for image denoising. Google
Engine Engine is used during the image preparation.

###  4.2 Image pair change detection

To evaluate the change detection performances and validate the effectiveness
of $S_{GLR}$ method, the image pair change detection results are compared with
Conradsen’s method [7, 8] and $CGLRT$ method [6].

|
---|---
(a) | (b)
|
(c) | (d)
Figure 1: False positive vs true positive curves comparison based on simulated
SAR images. (c) ROC curve results using simulated Sentinel-1 data (a), (d) ROC
curve results using simulated TerraSAR-X data (b). Multilooked data with
$3\times 3$ window size, 2SPPB and RABASAR provided data are used for the
comparison of Conradsen’s method, $CGLRT$ method and $S_{GLR}$ method,
respectively. Different colours represent different object changes:
green=farmland, yellow=forest, red=appearing, blue=appearing then
disappearing, cyan=disappearing.

During the acquisition of SAR time series, changes may happen for different
kinds of objects, such as buildings, farmland and forest, etc. Generally,
different kinds of objects have different change magnitudes. To
comprehensively and quantitatively evaluate the performances of different
methods, we processed the simulated SAR images which have different kinds of
object changes.

Because of the small change magnitude in SAR intensity images, identifying
forest area changes is much harder than farmland and building changes. It is
obvious that $S_{GLR}$ method can obtain the best detection results. In
addition, $S_{GLR}$ method can obtain better results with RABASAR-DAM provided
data. Compared with other methods, the $CGLRT$ method is good at detecting
building area changes which have large change magnitude. This characteristic
causes $CGLRT$ method to provide worse results when the false positives are
larger than $1\%$. According to the results shown in Figure 1, Conradsen’s
method always have more false positive because of using multilooked data.

|
---|---
(a) Noisy image acquired on 06/05/2011 | (b) Noisy image acquired on 08/06/2011
|
(c) Optical image acquired on 06/04/2011 | (d) RGB combination of the noisy images
|
(e) MIMOSA with false alarm rate equal to $1\%$ | (f) $CGLRT$ with 2SPPB data
|
(g) $S_{GLR}$ with RABASAR-DAM data | (h) $S_{GLR}$ with RABASAR-DBWAM data
Figure 2:  $Sendai$ SAR image pair change detection comparison. Pink
represents disappearing areas, while green represents appearing areas. The
comparison areas are indexed using red circles with associated numbers.

To fairly compare with the state-of-the-art change detection methods, we
processed popularly used TerraSAR-X images acquired over $Sendai$. Although
MIMOSA [42] detects all the changes directly using the noisy data, there are
too much wrong detections in the unchanged areas. The detections are seriously
influenced by the noise. $CGLRT$ can provide good results, but the changed
area boundaries are blurred compared with $S_{GLR}$. In addition, $CGLRT$
provides some wrong results in the water area.

24/12/2014-05/05/2015 |  |  |
---|---|---|---
24/12/2014-05/05/2015 |  |  |
24/12/2014-05/05/2015 |  |  |
24/12/2014-05/05/2015 |  |  |
| (a) | (b) | (c)
| decrease  increase
Figure 3: Continuous change monitoring results. (a) the reference image
divided by the other images, (b) changes detection results, (c) changes with
change magnitude weights. The thresholds were chosen for column (b) with a
false alarm rate equal to 0.54% and 1% for column (c). The time intervals of
the changes are shown on the left. Different colours represent the decrease
and increase in the magnitude of the backscattering values. The temporal
images are denoised using RABASAR-DAM.

$S_{GLR}$ method obtains the best results both with RABASAR-DAM and RABASAR-
DBWAM. They all precisely detect the 5 typical changed areas (Figure 2 (c)).
However, there is some noise inside the changed areas when using RABASAR-
DBWAM. $S_{GLR}$ method obtained better results for the objective change
detection with RABASAR-DAM.

###  4.3 Continuous change monitoring

Continuous change monitoring is a good way to track the development of object
changes. In this section, 5 Sentinel-1 images are processed using $S_{GLR}$
method, which provides satisfying multitemporal change detection results.
Yellow lines are used to highlight the farmland area boundaries, as shown in
the first column of Figure 3. Compared with the farmland backscattering values
in the reference image, the others seem to have smaller values. Thus, we
create the background image by dividing the reference image with other images,
so as to highlight the changed farmland areas.

There are valleys in this area, and even the farmland areas are not flat. The
threshold is defined empirically, according to the detecting change types or
changed areas. For example, big threshold values lead to the detection of high
change magnitude areas. With the false alarm rate equal to 0.54% (computed
using the denoised simulated SAR images without change), we can detect the
farmland area changes (Figure 3).

The positive values (red) indicate the increase of backscattering values
according to the reference data, and the negative values (blue) represent the
decrease of the backscattering values. The appearing or disappearing buildings
always have a large change magnitude. Seasonal changed areas (like farmland
areas, and some kinds of forest areas) have dynamic changes during the
timeline. With the false alarm rate equal to 0.54% (Figure 3), the proposed
method detects 83.23% of the appearing and disappearing buildings.

Since farmland areas and building areas have different mean intensity values,
we could suppress the detection of farmland area changes by adding this
information as a weight, with
$\exp({\sqrt{\hat{u}_{t}}+\sqrt{\hat{u}_{t^{\prime}}}}/{2})S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})$.
The exponential function is used to enlarge the backscattering value
differences between different objects.

In addition, these two areas usually have different change magnitudes. For
example, with weights calculated using the log version distance of the
correponding amplitude values:

$S^{D}_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})=\log(|\sqrt{\hat{u}_{t}}-\sqrt{\hat{u}_{t^{\prime}}}|)S_{GLR}(\hat{u}_{t},\hat{u}_{t^{\prime}})$
(28)

we could acquire new change detection results (third column of Figure 3).
However, after multiplying the weights, the threshold has to be defined
empirically.

###  4.4 Change classification

In this section, the proposed change classification method is compared with
NORCAMA [6] with the use of multitemporal Sentinel-1 images. This process can
distinguish farmland area changes from building area changes, which have
seasonal and non-seasonal changes, respectively.

####  4.4.1 Change classification with Sentinel-1 data

With high-frequency acquisition data, the whole duration of the changed
buildings will be monitored. During the construction of the building, its
backscattering values may keep changing which leads to complex change
monitoring results. To suppress the complex or high-frequency changes in
construction areas and keep the cycle changes of farmland areas, we under
sampled the frequency of the time series. The acquisition time of the
Sentinel-1 time series are 24/12/2014, 05/05/2015, 25/11/2015, 05/04/2016,
02/10/2016 and 18/01/2017.

However, all the Sentinel-1 IW VV polarization SAR images acquired through 110
orbit from 24/12/2014 to 18/01/2017 are used during the speckle reduction
process, so as to acquire better denoising results. RABASAR-DAM denoised
images and all change classification results are illustrated in Figure 4.
Since the actual backscattering values of the farmland areas are controlled by
the surface roughness and soil moisture, we could observe the weak backscatter
fields in SAR images acquired in spring. This phenomenon also reflects the
seasonal changes in the time series.

Compared to NORCAMA provided results, $S_{GLR}$ based change classification
method provides much better results. All the detected changed areas are well
corresponding to the previous cumulative change detection results demonstrated
in Figure 3. Visually, $S_{GLR}$ gives better detection results (Figure 4(h))
when using RABASAR-DAM provided data. There are fewer isolated points in the
detection results and the changed farmland areas are very smooth. The changed
types are similar to that using RABASAR-DBWAM provided data.

Global threshold is used for $S_{GLR}$ based change type detection method, so
as to speed up the processing. The denoising results of RABASAR seem good when
using denoised binary weighted arithmetic mean image, but the change type
detection results are not better than that using the denoised arithmetic mean
image. When using RABASAR provided data, NORCAMA method can obtain much better
results than using 2SPPB provided data.

|  |
---|---|---
(a) 24/12/2014 | (b) 05/05/2015 | (c) 25/11/2015
|  |
(d) 05/04/2016 | (e) 02/10/2016 | (f) 18/01/2017
|  |
(g) | (h) | (i)
Figure 4: Sentinel-1 time series change classification. (a-f) Sentinel-1 images, (g) NORCAMA, (h) $S_{GLR}$ with RABASAR-DAM provided data, (i) $S_{GLR}$ with RABASAR-DBWAM provided data. 6 images are used for the change type detection. The change type results are: white: no change, red: step change, green: impulse change, blue: cycle change and cyan: complex change. |
---|---
(a) Maximum value time | (b) Ground truth map
24/12/2014  18/01/2017
Figure 5: Maximum amplitude value time provided by REACTIV method. 10 Sentinel-1 images are used for the comparison. Ground truth map is prepared according to the arithmetic mean image, with different colours representing different objects. We mainly pay attention to the changed building areas (blue) and farmland areas (yellow, green). |
---|---
(a) Maximum value time | (b) Maximum change time
|
(c) Start changing time | (d) Stop changing time
24/12/2014  18/01/2017
Figure 6: Different change time detection comparison with improved REACTIV method. 10 Sentinel-1 images are used. The time series data are denoised by RABASAR-DAM with 69 time series images. A threshold is set at $99\%$ on the change probability.  |  |  |  |  |
---|---|---|---|---
06/05/2018 | 18/05/2018 | 30/05/2018 | 11/06/2018 | 23/06/2018
|  |  |  |
05/07/2018 | 17/07/2018 | 29/07/2018 | 10/08/2018 |
|
(a) Maximum value time | (b) Maximum change time
|
(c) Start changing time | (d) Stop changing time
2018/05/06  2018/08/10
Figure 7: Flooding area change time detection comparison with 9 Sentinel-1 GRD
data. 9 noisy images are shown above. The images are acquired over Xe-Pian Xe-
Namnoy dam in the southeastern province of Attapeu in Laos. 87 temporal
Sentinel-1 GRD images are used for the preparation of the arithmetic mean
image. All the test images are provided by RABASAR-AM.

The same crop with similar growth periods monitored in different years with
the same SAR sensor, they have similar backscattering values. It seems that
using fewer images is not enough to detect all kinds of vegetation interannual
variations. It is recommended to use $S_{GLR}$ change type detection method
with RABASAR-DAM provided data.

###  4.5 Change time detection and visualization with extended REACTIV method

In this section, the time series generated from Sentinel-1 single-look complex
images and Sentinel-1 GRD images are separately used to illustrate the
capability of the proposed method.

####  4.5.1 Farmland area and building area monitoring

To make a comprehensive evaluation of the method, we processed the temporal
Sentinel-1 SAR images with the original REACTIV method. The detection results
are shown in Figure 5(a). Although the changing area can be roughly detected
by the test statistics with original SAR data, the results are seriously
affected by the noise. The previously prepared ground truth map is shown in
Figure 5(b) to help the interpretation and evaluation of change detection
results.

With the RABASAR provided data, different change time detection strategies
described in Section 3.1 are utilized to process the Sentinel-1 time series.
The spatial adaptive denoising leads to spatially variable ENL. During the
change time detection test, we suppose they have the same ENL so as to speed
up the process. In addition, the red colour at the end of the colour bar is
removed so as to avoid the mix of red colours. All the time are compressed to
the set of colour index bar.

Compared with the original method, the improved REACTIV method can obtain much
better change detection results (Figure 6(a)). According to the ground truth,
the improved REACTIV method can detect the changed building areas and farmland
areas. With the default parameter, it can not detect the changes in forest
areas which have low change magnitude in the temporal amplitude SAR images.

With reference to different object changes, the maximum change magnitude time
and start changing time are the same for disappearing building areas. The
detection results in the appearing building areas are much more complex. Based
on the detection results, we can even distinguish different kinds of
farmlands. Since most farmland areas are not totally flat, this may cause
different parts of the same farmland to reach their maximum value at different
times. Therefore, the detection results in these areas are not smooth enough.
Generally, construction areas are larger than the under-construction building,
and this area may keep changing before the building is completed. All these
phenomena lead to complex change shapes in the detection results.

####  4.5.2 Monitoring abrupt floods in Southern Laos

Change detection is a significant application of remote sensing technology. In
this section, we try to apply the improved method to flooding area monitoring.
There exist a large number of Sentinel-1 GRD data over the test site, all the
images which have similar acquisition geometry are used to prepare the super-
image.

Sentinel-1 GRD images have the ability to monitor large-area changes. We only
test the improved method over the water storage area. The amount of flooding
water can be estimated according to the changes in the water area and local
digital elevation model. As shown in Figure 7, most of the areas are changed
between the image pairs which were acquired on 17/07/2018 and 29/07/2018. The
black areas surrounded by blue areas are the final water area during the image
acquisition period. The detection results have high similarity with the
results provided by ESA222ESA:
https://www.esa.int/Our_Activities/Observing_the_Earth/Copernicus/Sentinel-1/Sentinel-1_maps_flash_floods_in_Laos.

##  5 Conclusion

In this paper, we proposed a simple $S_{GLR}$ based similarity test which
could be applied and benefit to any denoised SAR images. The simple $S_{GLR}$
similarity is based on gamma distribution and used for the calculation of the
criteria map. Based on the prefiltered data, this method has been used for
image pair change detection, continuous change monitoring and change
classification. In particular, we mainly used RABASAR provided data in this
paper. The processing results of simulated and real SAR images show that
$S_{GLR}$ based change detection method provided good results both in the
processing of image pairs and temporal images. $S_{GLR}$ method gives better
change classification results compared to NORCAMA method. Using RABASAR-DAM
provided data, $S_{GLR}$ acquired much better change classification results,
with smooth changed areas and less noisy points.

In addition, we used the $S_{GLR}$ function and RABASAR denoising data to
improve REACTIV method. Based on the detection areas acquired by REACTIV
method (dynamics of time series coefficient of variation), we associated the
colours with different kinds of change times and changed the background with a
denoised image or arithmetic mean image. By only using part of the hue colour
channel, we successfully avoid the mixture of the red colour index. The
results obtained by the improved method provided useful information and allow
extended interpretation. The change time detection is much more effective for
homogeneous area changes and for abrupt changes, which is suitable for
monitoring farmland areas, flooding areas and some human activities (harbour
activities, urbanization and airport dynamics). However, the method has less
capability to detect seasonal changes in forest areas.

Future work will take into account the object attributes of the changed areas,
so as to acquire better analysis results. To precisely identify the seasonal
change of vegetation areas, we will pay attention to the multitemporal
coherence maps.

## References

  * [1] S. Khorram.  Accuracy assessment of remote sensing-derived change detection.  Asprs Publications, 1999.
  * [2] P. Coppin, I. Jonckheere, K. Nackaerts, B. Muys, and E. Lambin.  Review article digital change detection methods in ecosystem monitoring: a review.  International journal of remote sensing, 25(9):1565–1596, 2004\.
  * [3] T. Habib, J. Chanussot, J. Inglada, and G. Mercier.  Abrupt change detection on multitemporal remote sensing images: a statistical overview of methodologies applied on real cases.  In Geoscience and Remote Sensing Symposium, 2007. IGARSS 2007. IEEE International, pages 2593–2596. IEEE, 2007.
  * [4] L. Bruzzone and F. Bovolo.  A novel framework for the design of change-detection systems for very-high-resolution remote sensing images.  Proceedings of the IEEE, 101(3):609–630, 2013.
  * [5] P. Lombardo and T. Pellizzeri.  Maximum likelihood signal processing techniques to detect a step pattern of change in multitemporal SAR images.  IEEE Transactions on Geoscience and Remote Sensing, 40(4):853–870, 2002.
  * [6] X. Su, C.A. Deledalle, F. Tupin, and H. Sun.  NORCAMA: Change analysis in SAR time series by likelihood ratio change matrix clustering.  ISPRS Journal of Photogrammetry and Remote Sensing, 101:247–261, 2015.
  * [7] K. Conradsen, A.A. Nielsen, J. Schou, and H. Skriver.  A test statistic in the complex Wishart distribution and its application to change detection in polarimetric SAR data.  IEEE Transactions on Geoscience and Remote Sensing, 41(1):4–19, 2003.
  * [8] K. Conradsen, A.A. Nielsen, and H. Skriver.  Determining the points of change in time series of polarimetric SAR data.  IEEE Transactions on Geoscience and Remote Sensing, 54(5):3007–3024, 2016.
  * [9] L. Pulvirenti, M. Chini, N. Pierdicca, and G. Boni.  Use of SAR data for detecting floodwater in urban and agricultural areas: The role of the interferometric coherence.  IEEE Transactions on Geoscience and Remote Sensing, 54(3):1532–1544, 2016.
  * [10] M.T. Pham, G. Mercier, and J. Michel.  Change detection between SAR images using a pointwise approach and graph theory.  IEEE Transactions on Geoscience and Remote Sensing, 54(4):2020–2032, 2016.
  * [11] J. Zhao, M. Gong, J. Liu, and L. Jiao.  Deep learning to classify difference image for image change detection.  In Neural Networks (IJCNN), 2014 International Joint Conference on, pages 411–417. IEEE, 2014.
  * [12] M. Gong, J. Zhao, J. Liu, Q. Miao, and L. Jiao.  Change detection in synthetic aperture radar images based on deep neural networks.  IEEE transactions on neural networks and learning systems, 27(1):125–138, 2016.
  * [13] R. Touzi, A. Lopes, and P. Bousquet.  A statistical and geometrical edge detector for SAR images.  IEEE Transactions on geoscience and remote sensing, 26(6):764–773, 1988.
  * [14] C. Oliver and S. Quegan.  Understanding synthetic aperture radar images.  SciTech Publishing, 2004.
  * [15] V. Akbari, S.N. Anfinsen, A.P. Doulgeris, T. Eltoft, G. Moser, and S.B. Serpico.  Polarimetric SAR Change Detection With the Complex Hotelling-Lawley Trace Statistic.  IEEE Transactions on Geoscience and Remote Sensing, 54(7):3953–3966, 2016.
  * [16] Y. Bazi, L. Bruzzone, and F. Melgani.  An unsupervied approach based on the generalized Gaussian model to automatic change detection in multitemporal SAR images.  IEEE Transactions on Geoscience and Remote Sensing, 43(4):2972–2982, 2005.
  * [17] F. Bovolo and L. Bruzzone.  A detail-preserving scale-driven approach to change detection in multitemporal SAR images.  IEEE Transactions on Geoscience and Remote Sensing, 43(12):2963–2972, 2005.
  * [18] F. Bovolo and L. Bruzzone.  The time variable in data fusion: A change detection perspective.  IEEE Geoscience and Remote Sensing Magazine, 3(3):8–26, 2015.
  * [19] S. Aminikhanghahi and D.J. Cook.  A survey of methods for time series change point detection.  Knowledge and information systems, 51(2):339–367, 2017.
  * [20] L. Mou, L. Bruzzone, and X.X. Zhu.  Learning spectral-spatial-temporal features via a recurrent convolutional neural network for change detection in multispectral imagery.  arXiv preprint arXiv:1803.02642, 2018.
  * [21] E.M. Domínguez, E. Meier, D. Small, M.E. Schaepman, L. Bruzzone, and D. Henke.  A multisquint framework for change detection in high-resolution multitemporal SAR images.  IEEE Transactions on Geoscience and Remote Sensing, 56(6):3611–3623, 2018.
  * [22] A.A. Nielsen, K. Conradsen, and H. Skriver.  Change detection in a short time sequence of polarimetric C-band SAR data.  In ESA Living Planet Symposium 2016. European Space Agency, ESA, 2016.
  * [23] D. Amitrano, F. Cecinati, G. Di Martino, A. Iodice, P.P. Mathieu, D. Riccio, and G. Ruello.  An end-user-oriented framework for rgb representation of multitemporal sar images and visual data mining.  In Image and Signal Processing for Remote Sensing XXII, volume 10004, page 100040Y. International Society for Optics and Photonics, 2016.
  * [24] E. Koeniguer, A. Boulch, P. Trouvé, and F. Janez.  Colored visualization of multitemporal SAR data for change detection : issues and methods.  In Synthetic Aperture Radar (EUSAR), 20018 12th European Conference on. VDE, 2018.
  * [25] Weiying Zhao, Charles-Alban Deledalle, Loïc Denis, Henri Maître, Jean-Marie Nicolas, and Florence Tupin.  Ratio-based multitemporal sar images denoising: Rabasar.  IEEE Transactions on Geoscience and Remote Sensing, 57(6):3552 – 3565, 2019.
  * [26] J.W. Goodman.  Speckle phenomena in optics: theory and applications.  Roberts and Company Publishers, 2007.
  * [27] C.A. Deledalle, L. Denis, and F. Tupin.  Iterative weighted maximum likelihood denoising with probabilistic patch-based weights.  IEEE Transactions on Image Processing, 18(12):2661–2672, 2009.
  * [28] X. Su, C.A. Deledalle, F. Tupin, and H. Sun.  Two-step multitemporal nonlocal means for synthetic aperture radar images.  IEEE Transactions on Geoscience and Remote Sensing, 52(10):6181–6196, 2014.
  * [29] S.S. Wilks.  The large-sample distribution of the likelihood ratio for testing composite hypotheses.  The Annals of Mathematical Statistics, 9(1):60–62, 1938.
  * [30] T.W. Anderson and E.U. Mathématicien.  An introduction to multivariate statistical analysis, volume 2.  Wiley New York, 1958.
  * [31] C. Tison, J-M. Nicolas, F. Tupin, and H. Maître.  A new statistical model for Markovian classification of urban areas in high-resolution SAR images.  IEEE transactions on geoscience and remote sensing, 42(10):2046–2057, 2004.
  * [32] J. Shi and J. Malik.  Normalized cuts and image segmentation.  IEEE Transactions on pattern analysis and machine intelligence, 22(8):888–905, 2000.
  * [33] A.Y. Ng, M.I. Jordan, and Y. Weiss.  On spectral clustering: Analysis and an algorithm.  In Advances in neural information processing systems, pages 849–856, 2002.
  * [34] E. Schubert, M. Weiler, and H.P. Kriegel.  Signitrend: scalable detection of emerging topics in textual streams by hashed significance thresholds.  In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 871–880. ACM, 2014.
  * [35] U. Von Luxburg.  A tutorial on spectral clustering.  Statistics and computing, 17(4):395–416, 2007.
  * [36] E.C. Koeniguer, J-M. Nicolas, B. Pinel-Puyssegur, J-M. Lagrange, and F. Janez.  Visualisation des changements sur séries temporelles radar: méthode REACTIV évaluée à l’échelle mondiale sous Google Earth Engine.  Conférence Française de Photogrammétrie et de Télédétection (CFPT), 2018.
  * [37] Michèle Basseville, Igor V Nikiforov, et al.  Detection of abrupt changes: theory and application, volume 104\.  Prentice Hall Englewood Cliffs, 1993.
  * [38] J.W. Goodman.  Some fundamental properties of speckle.  Journal of the Optical Society of America, 66(11):1145–1150, 1976\.
  * [39] J-M. Nicolas.  Application de la transformée de Mellin: étude des lois statistiques de l’imagerie cohérente.  Rapport de recherche, 2006D010, 2006.
  * [40] M. Kendall and A. Stuart.  The advanced theory of statistics. Vol. 1: Distribution theory.  London: Griffin, 1977, 4th ed., 1977.
  * [41] J-M. Nicolas, E. Trouve, R. Fallourd, F. Vernier, F. Tupin, O. Harant, M. Gay, and L. Moreau.  A first comparison of Cosmo-Skymed and TerraSAR-X data over Chamonix Mont-Blanc test-site.  In Geoscience and Remote Sensing Symposium (IGARSS), pages 5586–5589, July 2012.
  * [42] G. Quin, B. Pinel-Puyssegur, J-M. Nicolas, and P. Loreaux.  MIMOSA: An automatic change detection method for SAR time series.  IEEE Transactions on Geoscience and Remote Sensing, 52(9):5349–5363, 2014.