# Opening the Black-Box: A Systematic Review on Explainable AI in Remote
Sensing

Adrian Höhl1∗, Ivica Obadic1,2∗, Miguel-Ángel Fernández-Torres3, Hiba Najjar4,
Dario Oliveira5,
Zeynep Akata6,7, Andreas Dengel4 and Xiao Xiang Zhu1,2 1Data Science in Earth
Observation, Technical University of Munich (TUM)  2Munich Center for Machine
Learning, 80333 Munich, Germany  3Image Processing Laboratory (IPL),
Universitat de València (UV)  4University of Kaiserslautern-Landau and German
Research Center for Artificial Intelligence (DFKI)  5School of Applied
Mathematics, Getulio Vargas Foundation  6Institute for Explainable Machine
Learning at Helmholtz Munich  7Interpretable and Reliable Machine Learning,
Technical University of Munich (TUM)  ∗Shared first authorship

###### Abstract

In recent years, black-box machine learning approaches have become a dominant
modeling paradigm for knowledge extraction in Remote Sensing. Despite the
potential benefits of uncovering the inner workings of these models with
explainable AI, a comprehensive overview summarizing the used explainable AI
methods and their objectives, findings, and challenges in Remote Sensing
applications is still missing. In this paper, we address this issue by
performing a systematic review to identify the key trends of how explainable
AI is used in Remote Sensing and shed light on novel explainable AI approaches
and emerging directions that tackle specific Remote Sensing challenges. We
also reveal the common patterns of explanation interpretation, discuss the
extracted scientific insights in Remote Sensing, and reflect on the approaches
used for explainable AI methods evaluation. Our review provides a complete
summary of the state-of-the-art in the field. Further, we give a detailed
outlook on the challenges and promising research directions, representing a
basis for novel methodological development and a useful starting point for new
researchers in the field of explainable AI in Remote Sensing.

###### Index Terms:

earth observation, explainable AI (xAI), explainability, interpretable ML
(IML), interpretability, remote sensing

##  I Introduction

Figure 1: The number of publications of ML and xAI in EO, obtained by using
the search query described in Section III-B and Appendix
LABEL:apx:search_queries, differ by a factor of $\approx 70$. (*Calculated
amount of publications given the first ten months of the year and assuming a
linear trend in 2023.)

Machine Learning (ML) methods have shown outstanding performance in numerous
Earth Observation (EO) tasks [1, 2], but mostly they are complex and lack the
interpretability and explanation of their decisions. In EO applications,
understanding the model’s functioning and visualizing the interpretations for
analysis is crucial [3], as it allows practitioners to gain scientific
insights, discover biases, trustworthiness and fairness for policy decisions,
and to debug and improve a model. Currently, the European Parliament is
adopting an Artificial Intelligence (AI) Act to ensure that the methods
developed and used in Europe align with fundamental rights and values such as
safety, privacy, transparency, explicability, as well as social and
environmental wellbeing [4]. It is anticipated that other governments
worldwide will implement similar regulations [5]. Many applications in EO
could potentially violate these values when data and AI are employed for
analysis and decision-making. Nevertheless, explainable AI (xAI) can
contribute to aligning these practices with rights and laws. Hence, xAI
emerges as a promising research direction to tackle the above-mentioned
scientific and regulatory challenges with observational data [6].
Despite these potential benefits, currently, there is a gap between the usage
of ML methods in EO and the works that aim to reveal the workings of these
models. This gap is illustrated in Figure 1, where the blue curve shows that
the number of ML papers in Remote Sensing (RS) has drastically increased in
the past few years. Although the works dealing with xAI in EO shown by the
green curve have also increased rapidly, there is still a gap by a factor of
$\approx 70$ compared to the number of papers of ML in RS. This increasing
number of xAI in RS papers motivates us to summarize the existing work in the
field and provide an overview to RS practitioners about the recent
developments, which might lead to narrowing this gap and making xAI approaches
more common in the field for RS.

xAI methods are typically designed to work on natural images. However, the RS
imagery has different properties than natural images [7]. First, images are
captured from above. This perspective comes with unique scales, resolutions,
and shadows. For instance, a RS image can cover whole landscapes, thousands of
square kilometers, while natural images can only cover a tiny fraction of it
[8]. Second, RS captures images in other electromagnetic spectra, apart from
the usual RGB channels. From hyperspectral over Synthetic Aperture Radar (SAR)
to Light Detection and Ranging (LiDAR), RS does cover a wide range of
reflectance data. Third, usual Red-Green-Blue (RGB) cameras are primarily
passive, while RS can be active, which changes properties like the radar
shadows, foreshortening, layover, elevation displacement, and speckle effects
[9]. Besides the different image properties, the tasks differ as well.
Computer Vision (CV) primarily addresses dynamic scenarios where objects can
be spatially covered by others. In contrast, the observed processes in RS
happen on different spectral and spatiotemporal scales, long-term and short-
term, and the systems modeled are very complex and diverse. Often, the
observed systems are not fully understood and can not be observed completely
or only indirectly. For instance, weather events, hydrology, hazards,
ecosystems, and urban dynamics. Although there are numerous reviews for xAI in
the literature [10, 11, 12], they typically do not reflect on the works
specific for RS nor reveal how the existing xAI approaches tackle the above
challenges related to remote sensing data. Therefore, a review of xAI tailored
to the field of RS is necessary to reveal the key trends, common objectives,
challenges, and latest developments.

While current reviews of xAI in EO are focused on social, regulatory, and
stakeholder perspectives [13, 14], specific subtopics in RS [15], or do not
provide a broad literature database [16, 15, 17], this paper targets the
applications and approaches from xAI in RS and follows a systematic approach
to gather a comprehensive literature database. To this end, we conduct a
systematic literature search for xAI in RS in three commonly used literature
databases in RS. From this literature database, we identify the relevant
papers and provide an overview of typical usages, objectives, and new
approaches in the field. Additionally, we propose a categorization of the
employed xAI methods, which allows to get a detailed overview and
understanding of the xAI taxonomy and techniques. We found that xAI methods
for RS are most often used for landcover monitoring, agriculture monitoring,
and natural hazard monitoring. They cover 40% of all 207 papers, with 34, 26,
and 26 papers, respectively. Overall, SHapley Additive exPlanations (SHAP) is
the most frequent method, used in over 38% of all publications. However, we
observe distinct variations of xAI methods in the EO tasks. For instance,
besides local approximation methods (29%), backpropagation methods (23%),
particularly Class Activation Mapping (CAM) methods (18%), are very prominent
in landcover monitoring. Conversely, behind local approximation (22%),
perturbation methods (22%) are widely utilized in interpreting agricultural
models. It should be noted, however, that most explanation outcomes are not
evaluated quantitatively and over 90% of the authors only provide anecdotal
evidence. Our work also includes a table with concise information about each
paper included in our study. Furthermore, we discuss the alignment of the
usage of xAI in RS with standard practices in xAI and the evaluation of xAI in
RS. As such, this review aims to assist users in the field of EO for the usage
of xAI. Finally, the combination of xAI with related fields, the unique
properties of RS data, and the interpretability of Deep Learning (DL) and
their implication into the lack of labels in RS are the challenges and
limitations we identified while conducting this review.

In summary, this review provides the following main contributions:

  * •

First, we present an overview and categorization of current xAI methods.

  * •

Second, we summarize the state-of-the-art (SOTA) xAI approaches in RS through
the analysis of a comprehensive literature database.

  * •

Third, we identify objectives and practices for evaluating xAI methods in RS.

  * •

Finally, we discuss challenges, limitations, and future directions for xAI in
RS.

The remainder of this review is structured as follows. We discuss the related
work and draw the distinction between existing xAI reviews and ours in Section
II. The methodological approach for our systematic review and our research
questions are introduced in Section III. Then, in Section IV, the taxonomy of
xAI is clarified, a categorization of the methods is developed, and the most
common xAI methods are explained. The results of this review are presented in
Section V. They show the current state of xAI in the RS literature, trends,
and novel approaches for RS tasks. Last but not least, Section
LABEL:sec:discussion discusses these results and how they align with the
recommended practices in xAI. Furthermore, connections to current challenges
and limitations are drawn. In the end, we summarize the main takeaways and
conclusions from our review and present an outlook.

##  II Related Work

Numerous resources for either the field of xAI [18, 10, 19] or ML for EO
applications [1, 20] are available in the literature. In contrast, to the best
of our knowledge, in the overlapping area of these two fields, there are only
two reviews [13, 17]. [13] aims to summarize the existing works of xAI in EO
and addresses the xAI usage from a regulatory and societal perspective,
discussing the requirements and type of xAI that is needed in EO from policy,
regulation, and politics [13]. However, only a small literature database is
used in this study and the author relies on a high-level categorization of the
xAI methods. As such, this work does not provide a comprehensive overview of
current xAI approaches in EO. Further, the work of [17] categorizes the
identified works in xAI according to the general challenges in the bio- and
geosciences [17]. Their categorization of xAI properties and their emphasis on
considering expert knowledge constitute two highlights of this review.
Furthermore, the presented challenges are still faced by researchers today.
Yet, they do not use a broad literature database to extensively cover the xAI
approaches used in EO. There also exist xAI reviews for certain EO tasks or
perspectives [14, 15, 16]. Still, they do not comprehensively summarize the
work done in the broad research field of RS in EO because they focus only on
specific parts of the literature. [14] conduct a review on the stakeholders
and goals within human-centered xAI applied in RS [14]. Their findings
indicate an underrepresentation of non-developer stakeholders in this area.
[15] [15] review DL methods and investigate to which degree the methods can
explain human wealth or poverty from satellite imagery. [16] [16] focus on xAI
in conjunction with Deep Neural Networks (DNNs) that incorporate geographic
structures and knowledge. They discuss three challenges when applying xAI to
geo-referenced data: challenges from xAI, geospatial AI, and geosocial
applications. Based on a short use case on land use classification and relying
on SHAP explanations, they show that the geometry, topology, scale, and
localization are of great importance. Finally, a group at Colorado State
University published a survey on their work using xAI for climate and weather
forecasting [21].

To overcome these shortcomings, we approach xAI in EO systematically and
provide an extensive literature database, resulting in a comprehensive summary
of the current literature. This work is characterized by our interest in the
application and usage of xAI, putting a special emphasis on RS sensors, while
others focus on the application of geological features, natural sciences, or
social implications. Not only do we present an overview of the current
challenges in the field, but we also highlight the state-of-the-art methods
for tackling these challenges. We believe this could provide valuable insight
into current limitations faced in the field. Compared to the existing
literature, we look at the topic from a technical perspective without
reflecting on regulatory or ethical implications originating from integrating
xAI into the field of EO. Because of the contrasting approaches relative to
xAI, this review excludes works from related and overlapping domains, such as
physics-aware ML, uncertainty quantification, or causal inference. Instead, we
refer the reader to overviews in their fields: [22], [23], and [24, 25],
respectively.

##  III Research method

In this review, we aim for transparency and reproducibility of our work by
following the PRISMA scheme [26]. We apply the appropriate PRISMA requirements
for our field to provide a transparent, complete, and trustworthy review.
Since literature search engines, such as Google Scholar, do not provide the
same search results for all users (the results depend on the geographic
location and time) [27], we avoided the use of such search engines and relied
on databases where the results can be reproduced. In addition to the PRISMA
guidelines, we provide the full search results and the exclusion criteria in
the supplementary material to further improve transparency and
reproducibility.

###  III-A Research Objective

The final goal of this survey is to summarize and compare the work conducted
in the field of xAI and EO with a focus on RS. Furthermore, we want to show
where and how the xAI methods are adopted, why they are used, and which
problems arise when applying them. In particular, we will attempt to answer
the following research questions:

  * •

RQ1 - Which explainable AI approaches have been used and which methods have
been developed in the literature for EO tasks?

  * •

RQ2 - How are xAI explanations analyzed, interpreted, and evaluated?

  * •

RQ3 - What are the objectives and findings of using xAI in RS?

  * •

RQ4 - How do the utilized xAI approaches in RS align with the recommended
practices in the field of xAI?

  * •

RQ5 - What are the limitations, challenges, and new developments of xAI in RS?

###  III-B Search procedure

The search query consists of two major parts: keywords related to xAI and
keywords related to EO. All the nested keywords are connected via an
$\operatorname{OR}$ operator, while the two parts are connected via the
$\operatorname{AND}$ operator. Due to the interchangeably used taxonomy in
both fields, we added additional keywords to the generally known terms to
receive as many relevant papers as possible, attempting not to excessively
increase the false positive rate. For example, we included specific types of
remote sensing sensors and used wildcards to cover different ways authors
might refer to xAI.

⬇

$\Bigl{[}$Earth observation OR remote sensing OR earth science OR
$\bigl{(}$(satellite OR aerial OR airborne OR spaceborne OR radar) AND (image
OR data)$\bigl{)}$ OR LiDAR OR SAR OR UAV OR Sentinel OR Landsat OR MODIS OR
gaofen OR ceres$\Bigl{]}$

AND

$\Bigl{[}$xai OR $\bigl{(}$(interpret* OR explain*) AND (deep learning OR
machine learning OR artificial intelligence OR dl OR ml OR ai OR
model)$\bigl{)}$$\Bigl{]}$

This general search query was adapted to the different search filters in three
databases: Scopus111https://www.scopus.com/,
Springer222https://link.springer.com/, and
IEEE333https://ieeexplore.ieee.org/. Appendix LABEL:apx:search_queries lists
the exact search queries for each database.

Figure 2: Flow Diagram of PRISMA. Figure adapted from [26].

The field of xAI started to gain momentum from 2014 [28]. However, we found
that search results of xAI papers in RS before 2017 are rarely relevant,
supporting that the methods have only been used in the last few years in this
area. Therefore, we searched journal and conference papers between 01.01.2017
and 31.10.2023, which is the last time the search was executed in the
databases. Hence, our review comprises more than the last six years. We
considered papers discussing RS in EO, namely RS sensors mounted to aerial
vehicles and satellites. We included all the available data from these
sensors, as well as high-level and compound products, such as Digital
Elevation Models (DEMs) and the ERA5 dataset [29]. By contrast, the filtering
or exclusion criteria we followed can be summarized as follows:

  * •

EC1 - Publication unrelated to explainable AI

  * •

EC2 - Publication unrelated to Remote Sensing

  * •

EC3 - Review/survey/short conference paper

  * •

EC4 - Publication published before 2017

  * •

EC5 - Publication not written in English

The covered work discusses the interpretability of the method or the
explanation results. Hence, it is not enough to only use a model which is
interpretable by itself, like a decision tree or linear regression (EC1). We
explicitly excluded in-situ measurements and pictures actively taken in the
surrounding environment, which are non-RS products, as well as papers out of
the scope of RS and xAI (EC2). Furthermore, review and survey papers were
excluded (EC3) since they would be included in the related work section. Also,
short conference papers that do not exceed a minimum of 5 pages are excluded
(EC3) because they usually contain preliminary or incomplete results. The
results of the search were filtered in three steps: (i) removing duplicates,
conference abstracts, and reports, (ii) screening through the abstract, and
(iii) screening through the full text. At least one author read each abstract
and paper; whenever there was doubt, the other authors read the corresponding
paper as well. In total, our search results in 1075 papers, merging the
different sources and removing duplicates, left us with 964 papers. In the
first shallow abstract screening, we removed 607 papers since they were
unrelated to our review topic. Three papers had to be excluded due to a lack
of access rights. After the full-text screening, we were left with 147 papers.
These papers were accompanied by 60 papers we had in our libraries. The
procedure is summarized in Figure 2. Reasons for removing a specific paper can
be found in the supplementary material.

##  IV Explanation methods in ML

This chapter provides a general overview of xAI. We first present the taxonomy
used to describe common distinctions between explanation methods. Then, we
introduce our categorization of xAI methods. Subsequently, we describe in
detail a selection of the methods commonly used in the field of RS. Finally,
we give an overview of the different metrics proposed in the literature for
evaluating these methods and present the main objectives of using xAI.

There exist several terms for explainable AI, such as interpretable Machine
Learning (IML) and interpretable AI. These terms often refer to the same
concept: the explanation or interpretation of AI models [30]. Hence, we will
use them interchangeably in this paper.
In the literature, three common distinctions exist when categorizing xAI
methods: _ante-hoc_ vs. _post-hoc_ , _model-agnostic_ vs. _model-specific_ ,
and _local_ vs. _global_ [12, 31, 32, 18]. The _ante-hoc_ and _post-hoc_
taxonomy refer to the stage where the explanation is generated. A xAI method
that provides interpretations within or simultaneously with the training
process is called ante-hoc. In contrast, a post-hoc method explains the model
after the training phase using a separate algorithm that could even be
distinct from the one followed by the model. _Model-agnostic_ methods have the
ability to generate explanations for any model. The method does not access the
model’s internal state or parameters, and the explanations are created by
analyzing the changes in the model’s output when modifying its inputs. _Model-
specific_ methods are exclusively designed for specific architectures and
typically have access to the model’s inner workings. _Local_ methods explain
individual instances and the model’s behavior at a particular sample. In
contrast, _global_ methods explain the model’s behavior on the entire dataset.
In practice, local explanations can be leveraged to achieve global
explanations. Through aggregation over a set of input instances chosen to
represent the dataset, local explanations can provide insights into the
general behavior of the inspected model. The aggregation mechanism should be
carefully defined since straightforward aggregation on some xAI methods might
lead to erroneous results. Some researchers further investigated the question
of finding meaningful aggregation rules of local explanations [33, 34]. The
methods presented in this review are categorized into the three universal
categories introduced above, as shown in Figure 3.

###  IV-A xAI methods categorization

In this section, we introduce a categorization of the most important xAI
methods in the literature. We build on the categorical foundation of [10] and
further adjust these categories to capture a larger set of xAI methods. The
categories are structured in a tree-like design. The tree has two internal
layers, which describe a hierarchy of primary and secondary categories and a
leaf layer which indicates an individual or group of specific methods. In the
following, a high-level overview of the categories of xAI methods is explained
following the structure of Figure 3. Our four primary categories are _feature
attribution_ , _distillation_ , _intrinsic explanations_ , and _contrastive
examples_. _Feature attribution_ methods highlight the input features,
significantly influencing the output. Alternatively, _distillation_ builds a
new interpretable model from the behavior of the complex model. _Intrinsic_
methods focus on making the model itself or its components inherently
interpretable. Lastly, _contrastive examples_ concentrate on showing simulated
or real examples and allow an explanation by comparing them. Each of these
categories splits into 2 to 3 subcategories.

For the sake of completeness, it should be noted that _feature selection_
methods are not considered in our categorization and are listed separately in
the results. Even though they are related to feature interpretation, they
serve different purposes. Feature selection can be defined as a strategy to
reduce the dimensionality of input space to improve the model performance and
reduce its computational cost [35]. While feature selection can constitute the
first step in a ML pipeline, feature interpretation is usually the last step
and typically involves more advanced techniques than only looking at the
predictive performance. Therefore, feature selection and xAI can be
complementary. On the one hand, feature selection reduces the input space that
needs to be interpreted. On the other hand, xAI can provide more qualitative
insights for selection, such as uncovering a bias introduced by a feature.

[xAI, for tree= child anchor=west, parent anchor=east, grow’=east, minimum
size=.65cm, draw, rounded corners, anchor=west, edge path= [edge] (.child
anchor) -— +(-5pt,0) – +(-5pt,0) —- (!u.parent anchor)edge label; , where
level=0text width=2cm, where level=1text width=3.4cm, where level=2text
width=4cm, where level=3text width=7cm, [Feature attribution,
[Backpropagation, [Activation Maximization [36], tikz=[draw=black, diamond,
fill=blue!30, inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black,
rectangle, fill=blue!30, inner sep=2.5pt, ] at([xshift=-15pt].north
east);[draw=black, circle, fill=black, inner sep=2pt, ]
at([xshift=-25pt].north east);] [Gradient [37], tikz=[draw=black, diamond,
fill=black, inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black,
rectangle, fill=blue!30, inner sep=2.5pt, ] at([xshift=-15pt].north
east);[draw=black, circle, fill=black, inner sep=2pt, ]
at([xshift=-25pt].north east);] [Integrated Gradients [38], tikz=[draw=black,
diamond, fill=black, inner sep=1.5pt] at([xshift=-5pt].north
east);[draw=black, rectangle, fill=blue!30, inner sep=2.5pt, ]
at([xshift=-15pt].north east);[draw=black, circle, fill=black, inner sep=2pt,
] at([xshift=-25pt].north east);] [Deconvolution [39], tikz=[draw=black,
diamond, fill=black, inner sep=1.5pt] at([xshift=-5pt].north
east);[draw=black, rectangle, fill=blue!30, inner sep=2.5pt, ]
at([xshift=-15pt].north east);[draw=black, circle, fill=black, inner sep=2pt,
] at([xshift=-25pt].north east);] [Class Activation Mapping (CAM) [40]
and variants (e.g. Grad-CAM [41]), tikz=[draw=black, diamond, fill=black,
inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black, rectangle,
fill=blue!30, inner sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black,
circle, fill=black, inner sep=2pt, ] at([xshift=-25pt].north east);] [Layer-
wise Relevance Propagation (LRP) [42], tikz=[draw=black, diamond, fill=black,
inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black, rectangle,
fill=blue!30, inner sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black,
circle, fill=black, inner sep=2pt, ] at([xshift=-25pt].north east);] ]
[Perturbation, [Occlusion Sensitivity [43], tikz=[draw=black, diamond,
fill=black, inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black,
rectangle, fill=black, inner sep=2.5pt] at([xshift=-15pt].north
east);[draw=black, circle, fill=black, inner sep=2pt, ]
at([xshift=-25pt].north east);] [Partial Dependence Plot (PDP) [44],
tikz=[draw=black, diamond, fill=blue!30, inner sep=1.5pt]
at([xshift=-5pt].north east);[draw=black, rectangle, fill=black, inner
sep=2.5pt] at([xshift=-15pt].north east);[draw=black, circle, fill=black,
inner sep=2pt, ] at([xshift=-25pt].north east);] [Accumulated Local Effects
(ALE) [45], tikz=[draw=black, diamond, fill=blue!30, inner sep=1.5pt]
at([xshift=-5pt].north east);[draw=black, rectangle, fill=black, inner
sep=2.5pt] at([xshift=-15pt].north east);[draw=black, circle, fill=black,
inner sep=2pt, ] at([xshift=-25pt].north east);] ] ] [Distillation [Local
Approximation [Local Interpretable Model-agnostic Explanation (LIME) [46],
tikz=[draw=black, diamond, fill=black, inner sep=1.5pt] at([xshift=-5pt].north
east);[draw=black, rectangle, fill=black, inner sep=2.5pt]
at([xshift=-15pt].north east);[draw=black, circle, fill=black, inner sep=2pt,
] at([xshift=-25pt].north east);] [SHapley Additive exPlanations (SHAP) [47],
tikz=[draw=black, diamond, fill=black, inner sep=1.5pt] at([xshift=-5pt].north
east);[draw=black, rectangle, fill=black, inner sep=2.5pt]
at([xshift=-15pt].north east);[draw=black, circle, fill=black, inner sep=2pt,
] at([xshift=-25pt].north east);] ] [Model Translation [Rule [48, 49], Tree
[50, 51, 52], Graph [53, 54], tikz=[draw=black, diamond, fill=blue!30, inner
sep=1.5pt] at([xshift=-5pt].north east);[draw=black, rectangle, fill=black,
inner sep=2.5pt] at([xshift=-15pt].north east);[draw=black, circle,
fill=black, inner sep=2pt, ] at([xshift=-25pt].north east);] ] ] [Intrinsic
[Interpretable by Design [Decision Rule and Decision Tree [55],
tikz=[draw=black, diamond, fill=blue!30, inner sep=1.5pt]
at([xshift=-5pt].north east);[draw=black, rectangle, fill=blue!30, inner
sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black, circle, fill=blue!30,
inner sep=2pt, ] at([xshift=-25pt].north east);] [Linear Regression [55],
Generalized Linear Model (GLM) [56],
Generalized Additive Model (GAM) [57], tikz=[draw=black, diamond,
fill=blue!30, inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black,
rectangle, fill=blue!30, inner sep=2.5pt, ] at([xshift=-15pt].north
east);[draw=black, circle, fill=blue!30, inner sep=2pt, ]
at([xshift=-25pt].north east);] [Latent Dirichlet Allocation (LDA) [58],
tikz=[draw=black, diamond, fill=black, inner sep=1.5pt] at([xshift=-5pt].north
east);[draw=black, rectangle, fill=blue!30, inner sep=2.5pt, ]
at([xshift=-15pt].north east);[draw=black, circle, fill=blue!30, inner
sep=2pt, ] at([xshift=-25pt].north east);] ] [Embedding Space [Attention
Mechanism [59, 60], tikz=[draw=black, diamond, fill=black, inner sep=1.5pt]
at([xshift=-5pt].north east);[draw=black, rectangle, fill=blue!30, inner
sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black, circle, fill=blue!30,
inner sep=2pt, ] at([xshift=-25pt].north east);] [Activation Assessment,
tikz=[draw=black, diamond, fill=blue!30, inner sep=1.5pt]
at([xshift=-5pt].north east);[draw=black, rectangle, fill=blue!30, inner
sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black, circle, fill=black,
inner sep=2pt, ] at([xshift=-25pt].north east);] [Concept Discovery (e.g. TCAV
[61], ACE [62]), tikz=[draw=black, diamond, fill=blue!30, inner sep=1.5pt]
at([xshift=-5pt].north east);[draw=black, rectangle, fill=blue!30, inner
sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black, circle, fill=black,
inner sep=2pt, ] at([xshift=-25pt].north east);] ] [Joint Training
[Explanation Association [63, 64], tikz=[draw=black, diamond, fill=black,
inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black, rectangle,
fill=blue!30, inner sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black,
circle, fill=blue!30, inner sep=2pt, ] at([xshift=-25pt].north east);]
[Prototype Learning [65, 66, 67], tikz=[draw=black, diamond, fill=black, inner
sep=1.5pt] at([xshift=-5pt].north east);[draw=black, rectangle, fill=blue!30,
inner sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black, circle,
fill=blue!30, inner sep=2pt, ] at([xshift=-25pt].north east);] [Model
Association (e.g. text explanations [68]), tikz=[draw=black, diamond,
fill=black, inner sep=1.5pt] at([xshift=-5pt].north east);[draw=black,
rectangle, fill=blue!30, inner sep=2.5pt, ] at([xshift=-15pt].north
east);[draw=black, circle, fill=blue!30, inner sep=2pt, ]
at([xshift=-25pt].north east);] ] ] [Contrastive Examples [Counterfactuals
[[69] [69], tikz=[draw=black, diamond, fill=black, inner sep=1.5pt]
at([xshift=-5pt].north east);[draw=black, rectangle, fill=blue!30, inner
sep=2.5pt, ] at([xshift=-15pt].north east);[draw=black, circle, fill=black,
inner sep=2pt, ] at([xshift=-25pt].north east);] ] [Example-based [[70] [70],
tikz=[draw=black, diamond, fill=black, inner sep=1.5pt] at([xshift=-5pt].north
east);[draw=black, rectangle, fill=blue!30, inner sep=2.5pt, ]
at([xshift=-15pt].north east);[draw=black, circle, fill=black, inner sep=2pt,
] at([xshift=-25pt].north east);] ] ] ] [below=1ex,draw,rectangle, rounded
corners, inner sep=5pt] at (current bounding box.south) post-hocante-
hoclocalglobalmodel-agnosticmodel-specific ; (BB)[use as bounding box, inner
sep=0, outer sep=0, fit to=tree]; (BB.north
west)rectangle([yshift=-3ex]BB.south east);

Figure 3: Categorization of xAI methods based on [10] [10]

####  IV-A1 Feature attribution

Feature attribution methods rely on the trained ML model to estimate the
importance of the input features. Depending on whether the explanations are
generated by inspecting the model internals or by analyzing the changes in the
model’s output after modifying the input features, the methods in this
category are further split into backpropagation and perturbation methods. The
output of these methods is commonly a saliency plot which determines the
contribution of the input features to the model prediction. In the case of
imagery inputs, the output of these methods can be visualized as a heatmap,
also called a saliency map, which highlights regions relevant to the model
prediction.

##### Backpropagation

These methods leverage the inherent structure of DNNs and estimate the
relevant features by propagating the output values at the top or intermediate
layers of the network to the input features. In the majority, these methods
compute gradients to this end. The _Deconvolution_ method, also known as
Deconvolutional Neural Network (NN) [39], is designed to reverse the
convolutional operations from Convolutional Neural Networks (CNNs).
Reconstructing the input space from the feature maps of the CNN allows
visualizing which information was learned and how the input is transformed
across different network layers. Similarly, _Layer-wise Relevance Propagation
(LRP)_ [42] calculates relevance scores for individual input features through
layerwise backpropagating the neuron’s activations from the output, utilizing
specialized propagation rules. The scores indicate the significance of the
connection between input and output. Various adaptations have been proposed
which apply different propagation rules based on the design of the networks
[71, 72, 73]. In contrast, the _Gradient_ or Saliency method [37] uses the
partial derivative with respect to the input to create the attribution maps.
Rather than computing the gradient once, _Integrated Gradients_ [38]
calculates the integral of gradients with respect to the input features along
an interpolation or path defined between a baseline input and the instance to
be explained. _CAM_ [40] visually explains CNNs through attribution heatmaps
by introducing a global pooling layer right before the top fully connected
one. Using the weights of the latter layer for a particular class, the heatmap
is generated by computing the weighted average of the activation maps in the
last convolutional layer before being upsampled to match the size of the input
tensor for explainability purposes. One extension of CAM is Gradient-weighted
Class Activation Mapping (Grad-CAM) [41], which replaces these weights by the
gradient of the output with respect to the last convolutional layer, thus
removing the original requirement of a final global pooling layer.

##### Perturbation

The perturbation methods assess feature importance by measuring the
sensitivity of model predictions to changes in the input features. These
methods are distinguished by how the features are perturbed. Among others,
perturbations include blurring, averaging, shuffling, or adding noise. For
example, the _Occlusion_ method [43] tries to remove a feature by occluding
the features with a neutral value. _Permutation Feature Importance (PFI)_ [74]
permutes the features along their dimension, destroying the original
relationship between input and output values. The _Partial Dependence Plot
(PDP)_ [44] method is designed to show the average influence of a single input
feature on the decision while marginalizing the remaining features, which are
fixed. Therefore, it assumes feature independence. A similar approach called
_Accumulated Local Effects (ALE)_ [45] can handle correlated features by
averaging over the conditional distribution.

####  IV-A2 Model distillation

Model distillation methods approximate the predictive behavior of a complex
model by training a simpler surrogate model that is usually interpretable-by-
design. By replicating the predictions of the complex model, the surrogate
model offers hypotheses about the relevant features and the correlations
learned by the complex model without providing further insights into its
internal decision mechanism. Distillation approaches are categorized into a)
local approximation methods, which train the surrogate model in a small
neighborhood around an individual local example, and b) model translation
methods, which replicate the behavior of the complex model over the entire
dataset.

##### Local approximation

Approaches in this category focus on explaining individual predictions of the
complex model by inspecting a small neighborhood around the instances to be
explained. In contrast to the backpropagation and perturbation methods, which
operate on the raw input features, the local approximation approaches
transform the input features into a simplified representation space, such as
superpixels for imagery inputs. A prominent approach in this category is
_Local Interpretable Model-agnostic Explanation (LIME)_ [46]. It creates a new
dataset in the neighborhood of the target instance by perturbing its
simplified representation. Next, an interpretable surrogate model is trained
to approximate the predictions of the complex model on this newly created
dataset. Hence, the explanation for the complex model is distilled to the
interpretation of the surrogate model. A similar strategy is employed by the
_SHAP_ framework [47]. Concretely, [47] [47] introduce Kernel SHAP, which
utilizes the LIME’s framework under specific constraints to obtain the feature
importance by approximating their Shapley values, a method grounded in game
theory for estimating the player’s contribution in cooperative games [75].
While Kernel SHAP is a model-agnostic approach, Shapley values can also be
approximated with model-specific approaches such as Deep SHAP [47] for neural
networks and Tree SHAP [76] which enables fast approximation of these values
for tree-based models.

##### Model translation

All the model’s decisions on the entire dataset are approximated into a simple
surrogate model. Typically, the interpretable-by-design methods summarized in
the next section are used as surrogate models, such as rule-based [77, 78],
tree-based [50, 79] or graph-based [53, 80].

####  IV-A3 Intrinsic

ML models, which provide an explanation by themselves based on their
structure, components, parameters, or outputs. Alternatively, a human-
interpretable explanation can be obtained by visualizing them.

##### Interpretable-by-Design

These methods are interpretable by humans because of their simplicity in
design, architecture, and decision process. _Decision Rules_ are hierarchical
IF-THEN statements, assessing conditions and determining a decision. Fuzzy
rules [81] are designed to address uncertainty and imprecision, which are
frequently encountered in nature and classical precise rules struggle to
represent, e.g., by partial membership of classes (fuzzy sets) [82]. Since the
proximity of these rules to natural language, they are interpretable by
humans. The _Decision Tree_ [83] greedily learns decision rules. Their
internal structure is a binary tree, where each internal node is a condition
and each leaf is a decision.

_Generalized Additive Models (GAMs)_ [57] are statistical modeling techniques
that approximate data using an additive function:
$g(E(y))=f_{1}(X_{1})+f_{2}(X_{2})+...+f_{i}(X_{i})+\beta_{0}+\epsilon$, where
$f_{i}$ are smooth, non-linear functions transforming the features $X_{i}$,
$\beta_{0}$ the intercept coefficient, and $\epsilon$ represents the error
term. On the other hand, _Generalized Linear Models (GLMs)_ [56] consider a
linear relationship defined by means of a specific distribution. Each $f_{i}$
becomes a coefficient $\beta_{i}$ and the function $g$ now computes a weighted
sum of features, which allows representing the mean of various exponential-
family distributions. A simple example of this model type is _Linear
Regression (LR)_ , which assumes a Gaussian data distribution where the
function $g$ is the identity. In this context, the model explanation can be
obtained by examining the coefficients.

A well-known approach for generative probabilistic topic modeling is _Latent
Dirichlet Allocation (LDA)_ [58]. A dataset is assumed to be organized in
corpora or collections, and each collection contains discrete units, such as
documents comprised of words. A distribution of these units characterizes both
the collections and the topics. By analyzing the proportions of the units in
the collections, LDA estimates the underlying topics. The decisions of the
Latent Dirichlet Allocation (LDA) model and the identified topics can be
interpreted by examining the predicted proportions for each collection.

##### Embedding space

These approaches process the activations in the latent space of DNNs to
interpret its workings. The _Attention_ Mechanism [59, 60] creates high-level
feature representations by using the attention weights to model the
dependencies between the different elements in the input. Hence, visualizing
the attention weights is a common procedure to assess the relevant features
for the model decisions. _Activation Assessment_ analyzes the activations in
the latent space of a NN based on projection techniques. Commonly used are
dimensionality reduction approaches [84] (e.g., t-distributed Stochastic
Neighbor Embedding (t-SNE) [85] and Uniform Manifold Approximation and
Projection (UMAP) [86]) or neuron receptive fields in CNNs [87]. _Concept-
based Explanations_ summarize the activations in the latent space in terms of
interpretable, high-level concepts. The concepts represent visual patterns in
computer vision, and they are either user-defined in Testing with Concept
Activation Vectors (TCAV) [61] or learned in an unsupervised manner [62].
These approaches provide global explanations by quantifying the concept
relevance per class.

##### Joint training

This category provides ante-hoc explanations by introducing an additional
learning task to the model. This task is jointly optimized with the original
learning objective and is used as an explanation. The methods in this family
typically differ based on the explanation representation and how the
additional task is integrated with the original model.

_Explanation association_ methods impose the inference of a black-box model to
rely upon human-interpretable concepts. A prominent approach in this category
is the concept bottleneck models that represent the concepts as neurons in the
latent space of the model. They introduce an additional task, first predicting
the interpretable concepts in an intermediate layer of a DL model. Then, the
model predictions are derived based on these concepts. During training, a
regularization term is added to the loss function that enforces alignment of
the latent space according to the interpretable concepts [63]. To estimate the
concept importance directly, [64] predict the final output by linearly
combining the concept activation maps [64].

_Prototype learning_ approaches aim to identify a set of representative
examples (prototypes) from the dataset and enable an interpretable decision
mechanism by decomposing the model predictions based on the instance’s
similarity with the learned prototypes [65]. Thus, visualizing the prototypes
enables global model interpretability, while the similarity with an input
instance offers local model explanations. One popular approach for prototype
learning on image classification tasks is the ProtoPNet architecture
introduced by [66] [66]. The prototypes represent image parts and are encoded
as convolutional filters in a prototype layer of the proposed network. Their
weights are optimized with the supervised learning loss of the network and
additional constraints that ensure both the clustering of the prototypes
according to their class and the separability from the other classes. One
extension of this approach is the Neural Prototype Trees [67], which organizes
the prototypes as nodes in a binary decision tree. Each node computes the
similarity of the corresponding prototype with an instance. These similarities
are used to route the instance towards the leaves of the tree containing the
class predictions.

In contrast to the previous approaches, which encode the explanation within
the model that performs inference for the original learning task, the _model
association_ methods introduce an external model that generates explanations.
These approaches are often utilized to provide textual explanations for CV
tasks. An example of such an approach is presented by [68], who derive text
explanations for self-driving cars based on jointly training a vehicle
controller and a textual explanation generator [68]. The vehicle controller is
a CNN model that recognizes the car’s movements with spatial attention maps.
Next, the explanation generator, which is a Long Short-Term Memory neural
network (LSTM) model, processes the context vectors and the spatial attention
maps from the controller to produce the text explanation.

####  IV-A4 Contrastive Examples

Methods within this category provide alternative examples to an input instance
and allow obtaining an explanation by comparing them. Usually, examples that
are close to each other in the input space yet lead to a different outcome
than the original input instance are shown.

##### Counterfactuals

This explanation type aims to discover the smallest change required for an
instance to achieve a predefined prediction. Essentially, they answer the
question, ”Why does it yield output X rather than output Y?” which is very
close to human reasoning [88]. They have a close proximity to adversarial
examples, although their objectives differ significantly. Adversarial examples
usually want to achieve a confident prediction with a minimal perturbed
instance, whose change should remain imperceptible for humans. Conversely,
counterfactuals aim to provide a diverse set of examples and should allow
representing the decision boundary of the model. [69] [69] introduce an
optimization problem whose primary objective is to find a counterfactual,
denoted as $x^{\prime}$, which is as close as possible to the original input
$x$. As such, the distance function $d(\cdot,\cdot)$ between the
counterfactual and the original input should be minimized. The corresponding
loss function is as follows:

$L(\lambda,x_{i},x^{\prime},y^{\prime})=\lambda(f_{w}(x^{\prime})-y^{\prime})^{2}+d(x_{i},x^{\prime}),$
(1)

where $y^{\prime}$ represents the desired output, $f_{w}$ is the model with
fixed weights, and $\lambda$ serves as a regularization parameter to balance
the contribution of the proximity of the original output to the desired output
$y^{\prime}$ with respect to the similarity between counterfactual and input.

##### Example-based explanations

Unlike counterfactuals, which can generate artificial instances, example-based
explanations usually present existing ”historical” training instances and
showcase similar instances to the input under consideration [84]. The user can
connect, correlate, and reason based on the analogies. The explanatory
approach aligns with case-based interpretable-by-design model explanations,
e.g., k-Nearest-Neighbors (kNNs). For example, [70] [70] train a skip-gram
model and evaluate the model using the nearest neighbors determined by the
distances in the embedded space. Furthermore, they illustrate that the
acquired word representations have a linear relationship, allowing for the
computation of analogies through vector addition.

###  IV-B Common Explanation methods in RS

In this section, we provide an in-depth explanation of several popular xAI
methods in the field of RS.

####  IV-B1 Class Activation Mapping (CAM) and Grad-CAM

[40] introduce CAM [40] to obtain visual heatmaps for class discriminative
localization in CNNs. This technique aims to identify important regions within
an image that influence the model’s decision toward each of the $C$ classes.
The approach requires the presence of a global average pooling layer between
the last convolutional layer and the top fully connected layer of the CNN
architecture. To calculate the class activation maps, the weights $w_{c,k}$ at
the fully connected layer associated with a particular class $c$ are used to
estimate the importance of each feature or activation map $A_{k}$ at the input
of the global average pooling layer. Finally, the saliency map highlighting
the discriminative regions for the class $c$ is computed as the weighted sum
of these $K$ activation maps with the following equation:

$\operatorname{S}^{CAM}_{c}=\sum_{k=1}^{K}w_{c,k}A_{k}$ (2)

[41] generalize CAM to a broader range of CNN architectures by introducing
Grad-CAM [41]. In contrast to CAM, which makes use of the weights
corresponding to class $c$, Grad-CAM utilizes the average gradient of the
logit $y_{c}$ for each of the $K$ feature maps in the last convolutional layer
to assess their importance for class $c$. Additionally, a ReLU function is
applied, propagating only the positive values. Therefore Eq. 2 becomes:

$\operatorname{S}^{Grad-
CAM}_{c}=\operatorname{ReLU}\bigg{(}\sum_{k=1}^{K}\Big{(}\frac{1}{N}\sum_{i}\sum_{j}\frac{\partial
y_{c}}{\partial A_{i,j}}\Big{)}A_{k}\bigg{)},$ (3)

where $(i,j)$ stand for the spatial coordinates of each of the $N$ locations
in every activation map $A_{k}$. Because the convolutional layers of the
networks usually reduce the input size, the activation maps need to be
upsampled to the original input size. However, this upsampling process can
lead to imprecisions and other backpropagation methods produce more fine-
grained heatmaps, which directly attribute to each pixel, such as Layer-wise
Relevance Propagation (LRP) or Integrated Gradients (IG). Nevertheless, recent
method revisions, such as LayerCAM [89], have aimed to produce more detailed
heatmaps.

####  IV-B2 Occlusion senstivity

[43] propose a perturbation-based method [43] to visualize the importance of
different image regions. Their approach slides a patch over the image and
observes the sensitivity of the model’s prediction. Different values for the
patch, its size, and sampling techniques can be considered. The assigned
importance is directly proportional to the drop in performance of the model
after occluding the patch. Consequently, this method is model-agnostic and can
be applied to any kind of architecture.

####  IV-B3 Local Interpretable Model-agnostic Explanation (LIME)

[46] proposed LIME [46], a model-agnostic method that approximates the
behavior of a complex model locally, in the neighborhood of a target instance.
Concretely, to explain the prediction of a complex model $f$ for a target
instance $x$, LIME performs the following main steps: 1) a dataset is created
around the neighborhood of $x$ by randomly performing perturbations on it
(e.g., adding noise, hiding or blurring parts of the input, etc.), 2) an
interpretable by-design surrogate model $g$ is trained on this dataset and 3)
the internals of $g$ are inspected to provide an explanation. It is important
to note that LIME generates explanations on a simplified representation space
that is interpretable to humans. For example, when the input is an image, the
simplified representation can correspond to a binary vector, indicating the
presence of superpixels decomposing the image. Therefore, in the first step,
LIME creates a dataset by perturbing the simplified representation of $x$.
This dataset is labeled according to the predictions of the complex model $f$
on the perturbed instances (the perturbed instances are reverse-transformed
into the original input representation before feeding it into the complex
model). In the second step, the surrogate model $g$ is trained on this dataset
by weighting the perturbed samples based on their similarity with $x$.
Finally, the internals of $g$ are inspected to explain the prediction of $f$
for $x$. For instance, in case $g$ is a linear regression model, its
coefficients can be used to assess the feature importance. In contrast, if $g$
is a decision tree, inspecting its rules can serve to explain the predictions
of the complex model $f$.

####  IV-B4 SHapley Additive exPlanations (SHAP)

The Shapley values [75] constitute an approach from cooperative game theory
used to estimate the importance of the input features for the prediction of a
ML model $f$ on a local instance, as well as their average marginal
contribution across all possible coalitions of features. Concretely, given an
instance $x$, the Shapley value $\phi_{i}$ for feature $i$ is computed as
follows:

$\phi_{i}(f,x)=\sum_{S\subseteq
F\backslash\\{i\\}}\frac{\left|S\right|!\left(F-\left|S\right|-1\right)!}{F!}\left[f_{x}\left(S\cup\\{i\\}\right)-f_{x}\left(S\right)\right],$
(4)

where $F$ is the set of all input features, $S$ is a coalition of features,
and $f_{x}\left(S\right)$ is the marginalized prediction over the features not
included in $S$ while the features in the coalition $S$ take the values of the
instance $x$. The Shapley values are considered to fairly assign the
contribution of the input features to the model prediction as they satisfy the
following properties: 1) the sum of the feature contributions adds up to the
difference between the model prediction for the instance $x$ and the average
model prediction on the dataset (efficiency), 2) the contributions of two
features are the same if they equally contribute to all coalitions (symmetry),
3) a zero contribution is assigned to the features that do not change the
model prediction (dummy), and 4) the Shapley value of a feature for an
ensemble of ML models can be computed by aggregating the individual Shapley
values across the models in the ensemble (additivity) [19].

The analytical solution for the Shapley values can be a computationally
expensive operation for models trained on more than a few features, as Eq. 4
requires iterating over all possible feature coalitions and computing the
marginal contribution in each coalition. Therefore, in practice, the Shapley
values are estimated with approximation techniques. [47] introduce the SHAP
framework [47] as a unified approach for model interpretability based on the
family of additive feature attribution methods. This family represents the
explanation through the coefficients of a linear model and has a similar set
of desired explanation properties to the ones of the Shapley values. Although
other xAI methods like LIME, LRP, and Deep Learning Important FeaTures
(DeepLIFT) [90] can be represented in the form of additive feature attribution
methods, the only additive method that satisfies the desired explanation
properties is the one having the Shapley values as its linear coefficients
(referred as SHAP values). This formulation enables the approximation of the
Shapley values with a model-agnostic approach based on LIME. Concretely, the
authors introduce Kernel SHAP, an approach that estimates the Shapley values
by constraining the LIME method to rely on a linear model as a surrogate model
and to use a specific similarity function for the weighting of the perturbed
instances. To leverage the internal structure of the ML models for a fast
approximation of the Shapley values, the authors also propose model-specific
approaches, namely Deep SHAP for DNNs and Tree SHAP [76] for tree-based
models.

####  IV-B5 Attention mechanism

Attention mechanisms represent an integral component of neural networks that
mimics cognitive attention [91]. Although initially proposed for natural
language processing tasks like machine translation and sentiment analysis [59,
60], they have also been leveraged recently in CV applications [92], as well
as for graph-structured data [93]. These modules induce attention weights on
the input features, determining how a neural network combines them to produce
a high-level feature representation. Thus, visualizing the attention weights
is a well-known approach to understanding the relevant features for the model
predictions and assessing the interaction of the input features in the context
of the learning task [94, 95].
An attention mechanism is usually defined for a query $\bf q$ and matrices of
key and value pairs $\mathbf{K}=[\mathbf{k}_{1},...,\mathbf{k}_{L}]$ and
$\mathbf{V}=[\mathbf{v}_{1},...,\mathbf{v}_{L}]$, respectively. It outputs a
high-level representation $\bf c$ depending on the values and based on the
alignment of the keys $\bf K$ with the query $\bf q$. The alignment function
is specified according to the attention weights and can be computed using
various functions proposed in the literature [96]. One of the most widely used
mechanisms is the _scaled dot-product attention_ introduced by [60] in [60],
where attention weights $\alpha$ are computed as follows:

$\alpha=\operatorname{softmax}(\frac{\mathbf{q}\mathbf{K}^{T}}{\sqrt{d_{k}}}),$
(5)

being $d_{k}$ the embedding dimension of the keys. Then, the high-level value
representation $\bf v$ is computed as a linear combination of the attention
weights and the value vectors $\bf c=\alpha\bf V$.

###  IV-C Evaluation of xAI Methods

Evaluating explanation quality and its trustworthiness is an essential
methodological challenge in xAI, which has received considerable attention in
recent years. The existing evaluation strategies can be categorized into (i)
functional approaches based on quantitative metrics and (ii) user studies
[97].

####  IV-C1 Functional evaluation metrics

These metrics assess the explanation quality by quantitatively describing to
which extent an explanation satisfies a certain set of desired properties. A
commonly evaluated explanation property is _faithfulness_ (also called
correctness) [98], which asserts how close an explanation method approximates
the actual model workings. Various functional metrics are proposed to evaluate
explanation faithfulness. For example, metrics based on randomization tests
are introduced in [99] to evaluate the explanation sensitivity to
randomization in model weights and label permutation. The results reveal that
most of the evaluated backpropagation methods do not pass these tests. Another
common approach to evaluate explanation faithfulness is based on the
perturbation of the input features. For instance, [100] measures the changes
in the model output after perturbation of the supposedly important features,
as estimated by the explanation method. Next, metrics based on perturbation
are also used to evaluate other explanation properties like _robustness_ (also
referred to as explanation sensitivity) which inspect the impact of small
perturbations of the input features on the resulting explanation [101, 102].
Explanations with low sensitivity are preferred as this indicates that the
explanation is robust to minor variations in the input. However, it is worth
noting that the perturbation-based metrics might result in examples with
different distribution than the instances used for model training, which
questions whether the drop in model performance can be attributed to the
distribution shift or to the perturbation of the important features. To
address this issue, [103] [103] propose model retraining on a modified dataset
where a fraction of the most important features identified by the xAI method
is perturbed. Further, an improved evaluation strategy that avoids the need
for model retraining based on information theory is presented by [104] [104].
Additionally, other evaluation approaches also include properties like
_localization_ , which determines whether the explanation identifies the
ground-truth region of interest, or _complexity_ , with the aim of measuring
the sparsity of an explanation [101]. For a detailed categorization of the
different explanation quality properties and a summary of the approaches used
in the literature to functionally evaluate them, the reader is referred to
[98].

####  IV-C2 User studies

Conversely, experiments in which humans evaluate the quality of explanations
can be conducted. User studies in [97] are further categorized into an
application-grounded evaluation and human-grounded evaluation, depending on
the evaluation task, the type of participants, and the considered explanation
quality criteria.

On the one hand, _Application-grounded_ evaluation studies typically involve
domain experts who evaluate the explanation in the context of the learning
task. For instance, the authors in [105] perform a case study where physicians
assess the predictions of a CNN model that classifies electrocardiogram
heartbeats. This study measures physicians’ agreement with the CNN model when
they are presented with (i) explanations computed with the LIME method and
(ii) nearest-neighbour examples obtained from the latent space of a trained
CNN model. The results indicate higher agreement for the nearest-neighbor
example explanations, which simultaneously enable the physicians to relate
model outcomes to clinically relevant concepts in contrast to the LIME feature
importance visualization.

On the other hand, _human-grounded_ evaluation consists of user studies where
the participants are typically non-domain experts that evaluate more general
notions of explanation quality. For example, [106] [106] evaluate whether
explanations computed with the LRP method help lay users recruited via an
online crowdsourcing platform to understand the decisions of CNN model for
image classification. The study’s conclusions identify a significant
statistical result, which shows that users predict the model decisions with
higher accuracy when presented with LRP explanations. A detailed survey on the
user studies for xAI methods evaluation is presented in [107].

###  IV-D xAI Objectives

In this study, we group the objectives of utilizing xAI in RS according to the
following four reasons for using xAI, as defined in [30]: (1) explain to
justify, (2) explain to control, (3) explain to discover, and (4) explain to
improve. _Explain to justify_ is motivated by the need to explain individual
outcomes, which ensures that the ML systems comply with legislations, such as
enabling users the ”right to explanations”. Furthermore, the explanations can
enable a detailed understanding of the workings of the ML model. Hence, the
_explain to control_ objective is relevant for assessing model trustworthiness
and can help to identify potential errors, biases, and flaws of the ML model.
These insights can be used to _discover_ scientific knowledge and new insights
about the underlying process that is modeled with the ML system or to further
_improve_ the existing model. [108] classify the improvement techniques based
on the xAI insights into (1) augmenting the input data, (2) augmenting the
intermediate features, (3) augmenting the loss function, (4) augmenting the
gradient, and (5) augmenting the ML model [108]. For the sake of completeness,
we consider adapting existing xAI methods as another improvement strategy.

##  V XAI in RS

Here we summarize the research on xAI in RS and answer the research questions
RQ1, RQ2, and RQ3. First, we bring attention to the common practices and
highlight new approaches (RQ1). Next, we summarize the methodologies to
understand and evaluate the model explanations (RQ2). Further, we outline the
common research questions across the different RS tasks that practitioners aim
to answer with xAI (RQ3).

###  V-A RQ1: Usage/Applications of xAI in Remote Sensing

Table LABEL:tab:summary contains the full list of publications included in our
study. The papers are arranged by considering groups for the different EO
tasks (see Appendix LABEL:sec:glossary) and following the xAI categorization
shown in Figure 3. xAI methods, objectives, and evaluation types are also
listed in the table.

Figure 4: Number of publications per different xAI methods, grouped according
to our categorization in Figure 3.

####  V-A1 Application of existing xAI methods to RS

Figure 4 shows the number of papers using a xAI method grouped by the
categories introduced, while Figure LABEL:fig:xai_models_categories_methods
(see Appendix LABEL:sec:add_eo_task_plots) illustrates all the combinations of
model and methods in the literature. Local approximation methods are the most
frequently used in over 94 publications. Through their popularity, they are
used on the most diverse set of EO tasks and models. In most cases, they
interpret tree-based models (i.e., Random Forest (RF) and tree ensembles),
followed by CNNs and Multilayer Perceptrons (MLPs). Also, with over 65%, most
of these publications rely solely on local approximation methods without
evaluating other methods. Backpropagation methods follow with a small gap of
72 papers, leveraging mostly CAM variants. While a lot of CNN architectures
are interpreted, most time series models, like the LSTM, are also located
here. In contrast, most papers using transformers are among the 35
publications that leverage embedding space interpretation techniques, which
was to be expected since attention is already the centerpiece of the
architecture. Many publications in this category also assess the models’
feature space. Followed by a large gap to the SHAP and LIME methods, 56
publications use various perturbation methods. The number of publications is
fairly well proportioned between the methods, but PFI is the most widely used,
followed by occlusion and PDP. 15 publications leverage a diverse set of
interpretable-by-design models. Although Generalized Linear Models (GLMs),
particular LR, are the most common, newer models like the Explainable Boosting
Machines (EBM) are gaining recognition. Overall, only a few papers employ
joint training, like prototypes or explanation associations. Even less popular
are model translations, counterfactuals, and example-based explanations.

The Mean Decrease in Impurity (MDI) or Gini importance is often used in
feature selection for global importance measurements and can be easily
obtained for tree-based methods [109, 110, 111, 112, 113, 114, 115, 116, 117,
118, 119, 120, 121, 122, 123, 124, 125]. This shows that xAI and feature
selection have fluent boundaries, as the Gini index is used for feature
selection, i.e., when deciding about the depth in the tree with the purity of
a split, but also allows for the interpretation of the decision process of the
model. Nevertheless, we will refer to them as feature selection methods.

##### Approaches in common EO tasks

(a) Landcover mapping.

(b) Agricultural monitoring.

(c) Natural hazard monitoring.

Figure 5: The number of times models, xAI categories, and xAI methods are used
for the most common EO tasks.

To provide a deeper insight, we will now discuss in detail the three most
representative EO tasks according to the number of papers, which are landcover
mapping, agriculture monitoring, and natural hazard monitoring. The graphs in
Figure 5 show these tasks and the frequency of the combination in the usage of
ML models and xAI methods. Figure LABEL:fig:xai_models_categories_methods in
Appendix LABEL:sec:add_eo_task_plots provides a detailed overview of the usage
patterns for all the EO tasks considered in this study.

_Landcover monitoring_
Due to the well-established datasets, the most prominent EO task is often
employed for evaluating or developing new xAI methods. Further, a low amount
of expert knowledge is needed to elucidate the model outputs for this task
compared to others, such as atmosphere monitoring or ecosystem interactions.
The majority of the 36 landcover mapping papers present CNN architectures,
leading to the favored use of backpropagation methods, especially from the CAM
family. For example, [126] leverage Grad-CAM to identify the reasons for
misclassified images [126], while [127] assess transfer learning and the
interpretation resilience with CAM [127]. The numerous publications proposing
novel or evaluating methods can be found in the corresponding Sections V-A2
and V-B2, respectively. Albeit with lower frequency than the CAM approaches,
the workings of the CNN models are also interpreted based on the local
approximation approaches. For example, SHAP is utilized to reveal the salient
pixels and the global importance of the spectral bands in [128], LIME is used
to improve model performance on the misclassified examples in [129], and both
approaches are utilized in [130] to validate the model predictions.

The second outstanding ML model for landcover monitoring is RF and exclusively
model-agnostic xAI methods are used for its interpretation. [112] explain the
correlations regarding land use for tea cultivation through LIME [112]. The
most frequently used approach, SHAP, is leveraged in [131], while PFI is
explored in [110] and [132] to rank reflectance data from various satellites
and DEM information. The DEM information seems to be valuable since both [131]
and [132] identify slope as the most important feature, while elevation is
ranked the highest in [110].

Unique to this task, the less common counterfactual and example-based methods
are applied to CNNs. [133] utilize a Generative Adversarial Network (GAN) to
generate counterfactual time series for the Normalized Difference Vegetation
Index (NDVI) of different landcover classes. The generator introduces
perturbations to instances and attempts to change the prediction of the pre-
trained classifier. Meanwhile, a discriminator ensures the quality of the
newly generated instance [133]. Moreover, [134] show examples of the dataset
belonging to the same class by computing the similarity in the latent space of
the last layer. These examples can be used to identify Out-Of-Distribution
(OOD) instances or cases where the classifier lacks generalization [134].

_Agricultural monitoring_
The primary tasks in this group are crop yield prediction and crop type
classification. Two additional tasks are tackled: irrigation scheme
classification [135] and lodging detection [136]. The methodologies employed
revolve around model-agnostic and backpropagation approaches, namely
perturbation, local approximation techniques, and CAM methods. These methods
are frequently integrated with gradient boosting and RF tree-based models.
Gradient boosting models are consitently interpreted with SHAP [137, 136, 138,
139, 140]. For example, [136] combine gradient boosting with SHAP to identify
maize lodging from Unmanned Aerial Vehicles (UAV) images [136]. The key
features identified are plant height computed from a DEM and Digital Surface
Model (DSM) and textural features from the Gray-Level Co-Occurrence Matrix.

Many different xAI strategies have been considered for crop yield prediction.
[138] assess the feature importance of weather and spectral bands across three
different crop phenology times on administrative boundary levels by
summarizing the SHAP values [138], while [137] examine the relevance of soil
features and DEM [137]. Both studies emphasize the importance of specific time
steps and that the spectral bands provide valuable information. [139]
investigate soil features, DEM, and spectral reflectance with a finer spatial
resolution at subfield yields. Interestingly, no spatially consistent limiting
factors are identified for an entire field, indicating the potential of
countermeasures on a subfield level [139].

Other studies utilize RF and perturbation methods, specifically PFI [141] and
ALE or PDP [142]. The feature importance of soil, weather, and spectral data
on yield potential is analyzed in [141], while PDP and ALE are used to assess
the interactions between management, weather, and soil data in [142]. Their
findings suggest that residue management or rate application decisions can
significantly influence crop yield.

Time series data is frequently employed for crop mapping, leading to the
utilization of Recurrent Neural Networks (RNNs) [143, 144, 145, 121] and
transformers [146, 147, 148, 121] architectures. Given LSTM models for crop
yield prediction, IG and SHAP are used to attribute soil moisture, weather,
and reflectance data in [144] and [145]. Both studies reach the same
conclusion: high temperatures during the growth season have a negative impact
on crop yield. LSTM with attention and transformers allow to differentiate
between corn and soybeans in [121]. Both models are interpreted not only
through their corresponding attention mechanisms but also considering
gradients and activation projection with t-SNE. Overall, both models agree in
their attribution, emphasizing the middle of the year as an important period
when corn starts to silk and soybeans begin to bloom. [149] apply CNNs to time
series for predicting crop yield from vegetation indices and weather data.
Attributions and scenario analysis for different weather conditions are
provided with a method derived from CAM and modified for regression, named
Regression Activation Mapping (RAM) [149].

Crop classification with prototypes through a CNN encoder is proposed by [150]
[150]. The LLP-Co method uses a priori proportions of the classes to match the
instance proportions assigned to the prototypes. A distinctive approach based
on a variational adversarial network for crop yield and irrigation scheme
classification is introduced by [135] [135]. To learn more meaningful latent
representations, the discriminator of the GAN architecture also needs to
classify the latent representations from the encoder into the correct classes.
Additionally, LRP is used to attribute the inputs. [117] use a model
translation method to identify RF subtrees, common rules within the RF model,
and their associated error rates [117]. Other works in agriculture monitoring
encompass Gaussian Processs (GPs) [151, 152] and EBM [153]. Notably, the
methods for agricultural monitoring are evaluated solely with anecdotal
evidence. No study employs quantitative metrics, and only [145] carry out a
small user study.

_Natural hazard monitoring_
The third main task consists of 26 papers on landslide susceptibility, fire
and flood monitoring, as well as geological hazard monitoring. The typical
features used as input for landslide susceptibility assessment are DEM
variables, landcover and vegetation (e.g., NDVI) information, and weather
variables such as temperature or precipitation. Often, they are supplemented
by human factors like the distance to roads or hydrological properties (e.g.,
drainage density, soil type). Usually, these models are interpreted with SHAP:
RFs [154, 155], Support Vector Machines (SVMs) [155], CNNs [156], Feed Forward
Neural Networks (FNNs) [156, 157, 156], and Gradient Boosting (GB) [158, 159].
Additionally, LIME, PFI or PDP are applied in the same fashion in [160, 156].
Other approaches include interpretable models [161, 162, 163].

In summary, rainfall [162, 164, 156], slope [161, 158], elevation, aspect
[164], curvature [154], distance to road [158, 160], and NDVI [154, 160] are
among the most important features for this task. Moreover, higher NDVI
suggests a decrease in landslide probability [160]. In contrast, mines can
increase the probability [162].

The flood mapping approaches rely on similar features as landslide
susceptibility (being rainfall a permanent component) and leverage mostly
tree-based models or CNNs with SHAP [165, 166, 167], PDP [167] or prototypes
[168]. The elevation and slope have the biggest influence on floods worldwide
and in Turkey, according to [165, 166], while rainfall, road density, and
building density have more influence on urban floods following [167]. [168]
detect floods by performing adaptive k-means clustering of the image pixels in
the latent space of a U-Net encoder [168]. This approach associates the
prototypes with the cluster centers and enables the interpretation of the
model decisions in terms of linguistic _IF … THEN_ rules.

Model-agnostic methods (SHAP, PFI, PDP) are mostly applied to explain fire
susceptibility prediction models (RFs and FNNs, see [169, 170]), which
typically rely on landcover, weather, and DEM data, as the former tasks. [170]
shows that a climate fire index (created other work with DNN), NDVI, and slope
are critical indicators for Italy, while [169] finds that humidity, wind
speed, and rainfall most are most important factors for Australia. Further,
[171] reveals that fires are more severe in areas with higher elevation and in
those where the dominant vegetation types are shrubs and open forests.

Various works explore geological hazards. For example, [172] uses CNN to
classify various disaster events (e.g., building damage, fire) from aerial
imagery and reveals the salient regions with a weighted combination of LIME
and SHAP attributions. Further, [173] utilizes Grad-CAM to identify the
mistakes of the CNN model for predicting volcano deformation patterns and uses
t-SNE to evaluate the differences between the latent space representations of
real and simulated data. Regarding other geological hazards, [174] explores
the tunnel geothermal disaster susceptibility based on land surface
temperature, river density, and other geological factors. They predict the
susceptibility based on an ensemble of machine learning models and evaluate
the factor importance with PFI, PDP, and LIME. Their analysis reveals that
land surface temperature, fault density, earthquake peak acceleration, and
river density are among the critical indicators. In another study, [175]
explores earthquake probability prediction with eXtreme Gradient Boosting
(XGBoost) and explains the model predictions with SHAP.

##### Distinct approaches across EO tasks

Besides the above-described common practices of applying xAI in RS, we also
identified distinct modeling approaches applied to specific EO tasks. For
instance, intrinsic CNN architectures based on the xAI category explanation
association are introduced for socioeconomic indicators estimation. The
semantic bottleneck models are particularly used to decompose the model
prediction into a linear combination of human-understandable concepts, as
described in Section IV-A. Consequently, these models explain the complex
relationship between the remote sensing scene and the socioeconomic indicators
through a set of interpretable concepts. In [176], a semantic bottleneck model
aims to predict environment scenicness which is shown to be an important proxy
for socioeconomic indicators, such as quality of life or health. The proposed
model is jointly optimized to first predict intermediate concepts
corresponding to landcover classes, followed by a linear layer that estimates
the scenicness based on the predicted landcover classes. This architecture
allows the authors to directly relate the scenicness prediction to landcover
classes which are human-understandable concepts. Similarly, for the problem of
estimating the living quality from aerial images, a semantic bottleneck model
is proposed in [177] with the interpretable concepts being different
liveability dimensions like population statistics, building quality, physical
environment, safety, and access to amenities.

Although interpretable-by-design models provide lower performances than DL
architectures, they are leveraged in diverse contexts. Linear regression is
used to support more complex methods and evaluates the linear trend for a
feature [142, 178]. The application of GLM is observed for small sample sizes
[179] and Generalized Additive Models (GAMs) are used for larger datasets
[180, 153, 161, 162]. [180] reduce the feature space with Principal Component
Analysis (PCA) in [180], while advanced versions based on EBMs are introduced
in [153, 161]. EBMs use pairwise feature interactions within the GAM and
consider gradient boosting to train each feature function consecutively. [162]
follows this idea and uses FNNs as function approximations [162]. A decision
tree incorporating linear regression models at the terminal leaves is employed
in [181] to gain scientific insights into the partitioning of precipitation
into evapotranspiration and runoff. [182] apply LDA to SAR images. Their bag-
of-words approach uses superpixels as words to do landcover mapping [182,
183]. Last but not least, [152] analyze the weights of the GP to find
anomalous samples [152].

_Fuzzy logic-based models_ are another approach that is mainly used to
evaluate the trustworthiness of ML models. An Ordered Weighted Averaging (OWA)
fusion function is presented in [184] for burned area mapping, which allows
controlling if the fusion results are affected by more false positives than
false negatives, and vice versa. At the same time, it foresees if there are
only a few highly or many low relevant factors when providing a particular
output. The outputs of Fuzzy Logic Systems (FLSs) for tree monitoring [185]
can also be easily validated. Last but not least, measure-, integral- and
data-centric indices based on the Choquet integral (an aggregation function
defined with respect to the fuzzy measure) are introduced in [186, 187, 188,
189] intending to develop more understandable ensembles for landcover mapping.

xAI methods are not only used to gain insights about the model behavior. Some
works explore their abilities to _classify or enhance training_ [190, 191].
[191] [191] leverage Grad-CAM to create masks that occlude features the
network has emphasized, encouraging the network to exploit other features
[191]. The output of three NNs is merged through attention in another study
[164]. The outputs and the attributions of the DeepLIFT backpropagation method
are the key values of the attention layer. In contrast, [190] [190] directly
classify objects on top of the CAM attribution map. Hence, their weakly
supervised method does not need bounding boxes, making the image labels
sufficient. They propose three distinct methods to segment the attribution
maps into bounding boxes, each performing differently depending on the number,
heterogeneity, and complexity of the objects. Furthermore, the same object
detection strategy is used by [192] to evaluate the xAI methods for their
localization abilities.

####  V-A2 Adapted xAI approaches

Figure 6: Number of papers adapting existing xAI methods, grouped per year.

The xAI methods commonly used in RS are originally designed to work on
classical computer vision datasets, such as ImageNet [193]. These are mostly
composed of natural images which, as described in Section I, significantly
differ from remote sensing acquisitions. This raises the question of whether
the utilized xAI methods fit remote sensing data well. In this respect, we
identified several works that propose new approaches considering remote
sensing data properties to produce better explainability insights. This is
typically performed by adapting the existing xAI methods or by proposing new
DL architectures tailored to RS applications.

Recently, as illustrated in Figure 6, there has been an increase in the
tailoring of xAI methods for RS properties. 50% of the publications have been
published in the last year, and no novel method was identified before 2021. A
strong focus is given to modifying the CAM and Grad-CAM methods for different
RS applications and domains. For example, [194] exploit that the target
objects in SAR images occupy only a small portion of the image to propose a
new CAM method which, instead of upsampling the feature map of the
convolutional layer to the input image, downsamples the input image to the
feature map of the last convolutional layer [194]. This operation results in
saliency maps that localize precisely the targets in SAR images compared to
the Grad-CAM method. Additionally, a CAM method able to produce much more
fine-grained saliency maps than the prior CAM methods is introduced by [192]
in [192]. Similar to Layer-CAM [89], they use shallow layers to get more fine-
grained results but also rely on scores, following the idea of Score-CAM
[195], which are not as noisy as gradients. Another attempt to improve current
CAM methods was proposed by [196], who utilize the attribution maps from all
network layers and decrease their number through only retaining maps which
minimize the information loss according to the Kullback-Leiber (KL)-divergence
[196]. Additionally, local attribution maps are generated by masking the image
and weighting the maps by the corresponding bounding box and their prediction
confidence. These local maps need to be smoothed with a Gaussian kernel to
avoid sharp boundaries in the resulting CAM. This novel approach, called
Crown-CAM, is evaluated on a localization metric and outperforms (augmented)
Score-CAM and Eigen-CAM on a tree crown localization task. CAM variants for
hyperspectral images are developed in [197]. The saliency map is now a 3D
volume instead of a 2D image, and each voxel attributes the different channels
in depth, which provides pixel-wise and spectral-cumulative attributions.
Other Grad-CAM adaptions proposed are a median pooling [198] and a pixel-wise
[199] variant.

Regarding deep learning approaches, a model prototype approach for RS is
proposed in [200], where the ProtoPNet architecture [66] is adapted to also
consider the location of the features. The network is iteratively trained in 3
stages. Firstly, the encoder and prototype layers are trained to produce the
prototypes. Secondly, the prototypes are replaced by the nearest prototype of
the corresponding class. Lastly, only the output layer weights are trained to
produce the final prediction. In contrast to ProtoPNet, the prototype
similarity is scaled with a location value learned by the network. This
acknowledges the location of the prototypes in the image and makes them
location-aware. Another approach is presented in [201], where a reconstruction
objective is added to the loss function to enable the Grad-CAM ++ method to
more accurately localize multiple target objects within an aerial image scene.

Finally, we also identified one approach addressing the model agnostic
methods. Specifically, [202] adapts the model-agnostic PFI method to
incorporate spatial distances [202]. Analogous to the PFI method, the features
are permuted, and the mean decrease in predictive accuracy is assessed.
Notably, features are permuted across various predefined distances, revealing
the spatial importance or sensitivity of the model.

###  V-B RQ2: Interpretation and evaluation of xAI explanations

####  V-B1 Understanding and validating explanations

The properties of remotely sensed data mentioned in Section I can hinder the
intuitive understanding of the semantics of the objects or the individual
pixels in the remote sensing scene. Therefore, an obstacle that arises when
applying xAI in RS is explanation interpretation, as the relevant features
often do not have a straightforward interpretation. We identify that this
challenge is frequently tackled by transforming the raw features into
interpretable indices used for model training [203, 144] or by associating
domain knowledge with the explanation at the post-hoc stage [146, 121].

##### The creation of interpretable input spaces

The data preprocessing in ML usually generates meaningful features for models.
Ensuring human-understandable features is essential for comprehending input-
output relationships or gaining knowledge of the model. For instance,
understanding xAI attributions of SAR images can pose challenges. However,
[199] address this by transforming the feature space into human-interpretable
factors using a U-Net architecture [199]. They derive three interpretable
variables from the VH polarization backscatter coefficients and the VV
polarization interferometric coherence of the Sentinel-1 images, providing
insights into the temporal variance and the temporal minimum. While the
temporal variance changes between different crops and landforms over time, the
temporal minimum is specific for flooded rice fields due to their proximity to
water. This facilitates the understanding of the attribution of the applied
Grad-CAM method.
A well-known preprocessing technique is the creation of spectral indices, like
the NDVI vegetation index. This approach leverages the known characteristics
of the various spectral channels and facilitates a more accessible
representation and interpretation of the features. Our study observes a
significant adoption of these indices.
Encoding the feature space through dimensionality reduction approaches is
commonly believed to decrease the interpretability. However, [132]
demonstrates that a helpful interpretation from a complex, correlated feature
space can still be maintained [132]. The author employs structured PCA for
feature space reduction and RF for classification. The xAI methods ALE, PDP,
SHAP, and PFI are evaluated both in the original and the transformed feature
space. The results indicate that the behavior of the features can be
identified in the principal components and allow the extraction of the
relationship between the main feature groups.

##### Explanation interpretation with domain knowledge

When raw inputs are used for model training, domain knowledge is often
utilized to reveal the semantics of the relevant features. It is usually
derived from already established indicators for the task under study, by
utilizing external geolocated data sources, or based on expert knowledge.

One task where existing indicators are commonly used is agriculture
monitoring, where insights about crop phenology are utilized to uncover the
semantics of the relevant features. Concretely, [146] associate the NDVI index
with the time points that are highly attended by the transformer encoder model
to reveal the key phenological events for crop discrimination [146].
Similarly, agronomical knowledge is considered in [121, 118] to reveal the
phenological stages by which the crops are classified. Further, the SHAP
values of the input features used in [138] for crop yield estimation are
aggregated according to the observation acquisition date into the stages of
crop development and crop harvest. This aggregation reveals the criticality of
the observations close to the harvest stage for the yield prediction. Similar
efforts for explanation interpretation can also be found for other tasks, such
as for predicting landslide susceptibility. For instance, [159] [159] relate
the spatial heterogeneity of the SHAP values to the natural characteristics
and human activities for the following factors: lithology, slope, elevation,
rainfall, and NDVI. They argue the differences in factor contributions can be
attributed to local regional characteristics such as topography, geology, or
vegetation.

Regarding geolocated data sources, we identified that landcover class labels
are associated with the input observations to interpret the explanation. [204]
identify the relationship between urban topology and the average household
income [204]. For this purpose, the Grad-CAM attributions of the image pixels
are related to their landcover classes. Commercial/residential units are
characterized by low income, while natural areas link to higher income.

In certain works, expert knowledge is used to validate or interpret the
findings. For example, the most important drivers of landslides noted by the
field investigation reports are compared to the features sorted by explanation
magnitudes in [164]. LDA is leveraged in [182] for unsupervised sea ice
classification, and closely related classes are identified by KL divergence.
The LDA derived topics and probabilities, together with the interclass
distances and segmented images, enable experts to assess the physical
relationship between these classes. For example, water bodies, melted snow,
and water currents have a similar topic distribution and a substantial
physical similarity: liquid water. In [205], a comprehensive approach is
employed for integrating expert knowledge to interpret and guide the process
of finding and validating dwelling styles and their evolution within ethnic
communities. Initially, the building footprint is extracted from satellite
imagery, and with the help of experts, the styles are classified into
different types. Subsequently, xAI is applied, leveraging XGBoost and SHAP to
determine the importance, with the experts inferring the semantic meaning. In
the last step, the experts guide the clustering of the proximity of the styles
and the analysis of the geographical distribution. The results reveal the
emergence of mixed ethical styles inheriting from the three traditional
styles, which can be correlated to migration records.

####  V-B2 Evaluation of xAI methods for RS

Figure 7: The number of times different evaluation types are considered in the
literature, based on their corresponding number of papers.

As indicated in Section IV-C, xAI evaluation poses an open challenge. Figure 7
illustrates that the majority of the literature relies on anecdotal evidence,
often involving the visualization of arbitrarily or cherry-picked examples. In
addition to this informal evaluation, some authors [200, 206, 207] evaluated
their methods on straightforward toy tasks where humans can easily identify
the sought relationships.

[145] exclusively conduct a user survey to assess (1) the importance of the DL
features by experts and (2) judge the importance by the xAI method [145]. Five
crop modeling experts assigned important scores to the features, subsequently
compared to the post-hoc SHAP importance scores. Afterward, the experts
categorized the model explanations into four categories (strong) agree,
(strong) disagree, and should provide a justification. Overall, it is
demonstrated that experts can understand the model explanations, and the
explanations enable the experts to get insights into the models. However, the
task remains challenging and has the potential for misconceptions about the
model behavior.

Regarding quantitative evaluation, we identified 16 studies testing the xAI
methods for RS problems with functional metrics. These metrics mainly asses
the explanation quality properties described in Section IV-C and the
localization ability of the xAI methods. The backpropagation methods are most
frequently evaluated with a particular focus on the CAM approaches.

The explanation quality properties are evaluated in [208, 197, 198, 209, 194,
210, 127]. [208] [208] evaluate Saliency, Input*Gradient (I*G), IG,
GuidedBackprop, (Guided) Grad-CAM, Occlusion, DeepLift, LIME, and Smooth
Gradient (SmoothGrad) on landcover mapping tasks. Utilized metrics are max-
sensitivity, file size, computation time, and Most Relevant First (MoRef).
Max-sensitivity measures the reliability (maximum change in explanation) when
the input is slightly perturbed, the file size is used as a proxy for
explanation sparsity, and MoRef measures how fast the classification accuracy
declines when removing the most relevant explanations. The results indicate no
obvious choice for this task. While Occlusion, Grad-CAM, and LIME were the
most reliable according to the max-sensitivity metric, they lack high-
resolution explanations. Grad-CAM also emerged as the most computationally
efficient choice. Evaluating faithfulness with MoRef is also conducted in
[197] where the performance of different CAM methods is compared. Also, a
class sensitivity metric is employed to measure the correlation between the
attributions of different target classes. Additionally, [198] leverages
median-pooling Grad-CAM, Grad-CAM, Grad-CAM ++, and SmoothGrad-CAM ++ on the
metric drop/increase in confidence when occluding important regions from the
image. Similar metrics are also evaluated in [194], where two different
occlusion tests are employed to benchmark various CAM approaches with proposed
Self-Matching CAM. While the methods perform similarly when the most
influential pixels are perturbed, only Self-Matching CAM observes a minor drop
in the prediction difference when only the most salient pixels are not
occluded, demonstrating that its explanation focuses on the target object in
the image Further, [209] develop attention NN and compare it to attention
networks from the literature, CAM, Grad-CAM and Layer-CAM. The metrics max-
sensitivity and average % drop/increase in confidence are utilized. In a
different study, Low faithfulness of the Grad-CAM explanations based on
similar metrics is demonstrated in [210] for the task of real estate
appraisal. This study also conducts model and data randomization tests to find
that Grad-CAM is sensitive to changes in network weights and label
randomization. A similar evaluation for sensitivity to changes in the network
weights is conducted in [127], where the goal is to assess how damaging the
network weights affect the inference of the CNN models. To evaluate this, the
authors compare the CAM maps of the original model against the CAM maps of the
model with disabled weights. Their results indicate that the ResNet model
exhibits the highest resilience to network changes as the CAM map model with
disabled weights matches most closely the CAM map of the original ResNet
model.

The localization ability is evaluated in [211, 192, 196] For instance, [211]
[211] use various CAM approaches to evaluate the ResNet model on the tasks of
dealing with large variance problems and its localization ability.
Particularly, segmentation maps were derived from CAM, Grad-CAM, Grad-CAM ++,
SmoothGrad-CAM ++, and Score-CAM approaches by thresholding their
attributions. Grad-CAM showed superior performance in localization accuracy
and the ability to extract complex features in images with large variance. In
[192], the localization capabilities of the new CSG-CAM method are compared to
CAM methods, Grad-CAM, Grad-CAM ++ and Score-CAM. To this end, a weakly
supervised segmentation task driven by the xAI attributions maps is
considered. If the attributions succeed a threshold for a certain class, the
class is assumed to be contained in the scene. Similarly, [196] evaluates the
localization ability of different CAM approaches for tree crown detection.

When it comes to the evaluation of other xAI methods, despite their high
usage, the SHAP explanations are quantitatively evaluated only in [212] where
a FNN trained using the features with the highest SHAP values slightly
outperforms the model trained on the complete set of features. Further, the
attention weights are evaluated in [146] by inspecting drops in the accuracy
for crop mapping when the transformer model is trained on a subset of dates
with the highest attention values. The results verify that attention weights
select the key dates for crop discrimination as training the model with only
the top 15 dates is sufficient to approximate the accuracy of the model
trained on the complete dataset. Finally, [133] [133] employ distinct metrics
for their counterfactual generation. These metrics aim to ensure a certain
quality of the counterfactuals. They used proximity (evaluating closeness to
the original input instance, measured by $l_{2}$ distance), compactness
(ensuring a small number of perturbations across time steps), stability
(measuring the consistency for comparable input samples) and plausibility
(measuring adherence to the same data distribution).

###  V-C RQ3: xAI Objectives and Findings in RS

(a) Number of times xAI objectives are considered in the publications.

(b) Number of times xAI objectives are considered in the publications per EO
task.

Figure 8: The frequency of the objectives for using xAI in RS identified in
the analyzed works, categorized according to the scheme introduced in Section
IV-D. Papers can appear in multiple categories due to ambiguous meanings or
multiple motivations and objectives.

In this section, we analyze the motivation for applying xAI in RS according to
the common objectives specified in Section IV-D. The frequency of the
objectives for using xAI displayed in Figure 8a reflects a similar trend as in
[13]. Namely, the objective to control is the most commonly found in 127
works. It is followed by a large gap in the objective to discover insights,
which is met in 51 works. Next, follow the objectives to improve and to
justify which are found in 27 and 5 works, respectively. Moreover, Figure 8b
shows that the objective to discover has a unique distribution across the EO
tasks: it frequently occurs in studies monitoring the atmosphere, vegetation,
and the human environment interaction. In contrast to other objectives, which
are often identified in the studies related to the three main EO tasks
described in Section V-A.

The observed distribution of the EO tasks with the objective _to discover_
indicates that these studies focus on knowledge extraction for more recent EO
tasks, which are used to deal with monitoring and understanding the driving
factors behind extreme events related to climate changes or natural hazards.
These studies are of utmost importance as their insights can enable the
application of early preventive measures for disaster management. For
instance, the problem of uncovering the key drivers for wildfire
susceptibility is tackled in [213]. Applying SHAP and PDP on the trained ML
model reveals that soil moisture, humidity, temperature variables, wind speed,
and NDVI are among the most important factors associated with wildfires.
Regarding monitoring natural hazards, [214] uses SHAP to identify that
volcanic deposits, terrain properties, and vegetation types are strongly
linked to vegetation vulnerability after volcanic eruptions. Concretely,
increased vegetation vulnerability is associated with higher lapilli
accumulations, being crops and forests the most and the least susceptible
vegetation types, respectively [214].
Extracting novel scientific insights can also be valuable to gain new
knowledge about ambiguously defined concepts. One such concept is wilderness.
Although monitoring wilderness areas can be an important indicator of
sustainable development, there is no uniform definition in the literature for
wilderness, which is usually described with several philosophical reflections.
[215] aim to infer wilderness characteristics from Sentinel-2 images. By
analyzing occlusion sensitivity maps, they reveal that wilderness is
characterized by large areas containing natural undisrupted soils, in contrast
to anthropogenic areas that have specific edge shapes and lie close to
impervious structures [215]. Besides discovering new knowledge, scientific
insights can also be used to identify the reasons for inaccuracies in Earth
system models. [216] provide insights into a lightning model’s structural
deficits by predicting its error with a gradient boosting algorithm and
interpreting it with SHAP. The error is computed as the difference between the
output of the model and satellite observational data. Their analysis reveals
potential deficits in high convective precipitation and landcover
heterogeneities [216].

On the other hand, the studies with the objective _to control_ are mainly
comparing the inference mechanisms of various established ML models used for
EO. Consequently, they are mostly conducted in landcover mapping. One example
of such study is presented in [217], where SHAP feature importance rankings
are compared across several ML models trained on Sentinel-2 images,
topography, phenology, and texture features. This study illustrates that the
used models consistently rank the relevance of the red edge bands for
predicting the agriculture class and the chlorophyll-sensitive bands for
identifying deciduous trees. Conversely, inconsistencies among the models are
found when analyzing the importance of the texture features.

Agricultural monitoring is the second most common task where this objective is
considered. These studies typically focus on assessing the reliability of the
proposed models by quantifying the relevance of the multitemporal information.
For instance, [147] use the gradient method to measure the temporal importance
assigned by various DL models for crop classification. Their analysis reveals
that the transformer and the LSTM approaches ignore the observations obscured
by clouds and focus on a relatively small number of observations when compared
to the CNN models [147]. Following a similar approach, [121] evaluate the
generalization capabilities of these models when inference is performed in
different years than the ones used for model training [121]. Their experiments
indicate that the LSTM model adapts better to changes in crop phenology
induced by late plantation compared to the transformer model.

The objective to control also supports different types of studies that
anticipate the model decisions in scenarios that can occur in practical
applications. For instance, [218] investigate the impact of exposing a DNN,
initially trained on cloud-free images, to cloudy images [218]. They apply
Grad-CAM to identify crucial regions for the classifier in both image types.
Their findings reveal different factors why the network misclassifies the OOD
examples, including the coverage of structures through cloud cover and
shadows, as well as the homogeneity or heterogeneity induced by different
cloud types. Forecasting weather extremes constitutes another application
where verifying that ML models generalize well enough to predict unlikely and
rare events is essential. [219] provide a showcase for extreme rainfall
downscaling [219]. By interpreting a CNN with Grad-CAM they show that the
model is able to learn the most relevant meteorological features.

The studies with the objective _to justify_ usually discuss why the insights
about the model workings are relevant for the policymakers. For instance,
[143] argue that the interpretability of the DL models in agricultural
applications is critical to ensure fair payouts to the farmers according to
the EU Common Agricultural Policy (CAP) [143]. By applying a perturbation
approach, the authors find that the summer acquisitions and the red and near-
infrared Sentinel-2 spectral bands carry essential information for the
decisions of the used RNN model. Further, the human footprint index, which
represents the human pressure on the landscape and can be a valuable metric
for environmental assessments, is predicted from Landsat imagery in [220]. LRP
is then leveraged to visually highlight the relevant features in the images.
Besides policy-makers, the approach presented in [134] can also assist
individual users during the production phase of a ML model, as it justifies
the model’s validity for inference by providing example-based explanations. If
the explanation example does not fit the input instance, the model is
considered unreliable for the RS image classification task.

Lastly, a large part of the works with the objective _to improve_ focus on
adapting the existing xAI methods to RS tasks and are already described in
Section V-A2. Concerning the other techniques for model improvement based on
xAI insights, we identified that data augmentation and model augmentation are
commonly performed.

Data augmentation is applied by [221], who simulate synthetic data for
training a CNN model for the detection of volcanic deformations [173]. By
conducting an explainability analysis with Grad-CAM on real data, the authors
identify that the model wrongly predicts volcanic deformations on, e.g., salt
lakes and slope-induced signals, which are patterns that were not considered
in the simulated data. These insights are used to improve the prediction
performance by fine-tuning the last layer of the CNN model on a hybrid
synthetic-real dataset that accounts for these patterns. An iterative model
improvement for satellite onboard processing through a weakly supervised
human-in-the-loop is proposed by [209] in [209]. In their scenario, the
satellite has a convolutional attention network for object classification
onboard. Initially, the uncertain explanations and images are identified by
assessing the similarity of the attention maps across the attention blocks. An
inconsistency metric is introduced to measure the similarity of the
attribution maps emphasizing commonly highlighted regions. If a sample
performs badly on that metric, experts on the ground refine it by labeling the
incorrect pixels in the attention map. In the last step, a local model is
retrained with the refined data, and the onboard model gets updated. A last
example of data augmentation is the work by [191] [191] where, for every
image, a mask based on the Grad-CAM output is used during training to occlude
those regions for which the model provides the highest activation, attempting
to encourage the network to explore other features in the image.

On the other hand, model augmentation is performed in [222], where xAI is
initially used to compare the workings of the Vision Transformer (ViT) and CNN
models for monocular height estimation. Upsampling the feature maps of the
penultimate layer to the input image unveils that even though the ViT learns
more disentangled representations than the CNN, a neuron in the penultimate
layer of the ViT still encodes semantics of multiple objects. This is resolved
by proposing a new transformer model that learns distinct latent
representations per semantic class.

TABLE I: A complete list of all relevant papers in this review, aggregated by EO Task and xAI Category. (†xy indicates a new method which was derived from xy; a full list of all used acronyms can be found in Appendix LABEL:sec:glossary.) EO Task |  xAI Category |  Paper, xAI Methods |  Model |  Evaluation |  |  |  Objective |  |  |
---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  Toy Task |  Anecdotal |  Quantitative |  Control |  Improve |  Discover |  Justify
Agricultural Monitoring |  Backpropagation |  [199] PWGrad-CAM†Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [149] RAM |  CNN |  |  ✓ |  |  |  |  |
|  Backpropagation, Embedding Space |  [147] Gradient, Attention, Activation Assessment |  Transformer, CNN, ConvLSTM, LSTM |  |  ✓ |  |  ✓ |  |  |
|  Backpropagation, Embedding Space, Feature Selection |  [121] MDI, Gradient, Attention, Activation Assessment |  RF, Transformer, aLSTM |  |  ✓ |  |  ✓ |  |  |
|  Backpropagation, Joint Training |  [135] Explanation Association, LRP |  GAN |  |  ✓ |  |  ✓ |  ✓ |  |
|  Backpropagation, Local Approximation |  [144] SHAP, IG |  LSTM |  |  ✓ |  |  ✓ |  |  |
|  |  [145] SHAP, IG |  LSTM |  |  ✓ |  |  ✓ |  |  |
|  Embedding Space |  [148] Attention |  Transformer |  |  ✓ |  |  ✓ |  |  |
|  Feature Selection |  [118] MDI, GFFS |  RF |  |  ✓ |  |  ✓ |  |  |
|  Interpretable by Design |  [153] EBM |  EBM |  |  ✓ |  |  ✓ |  |  |
|  Interpretable by Design, Perturbation |  [152] GP, Occlusion |  GP |  |  ✓ |  |  ✓ |  |  |
|  |  [142] GLM, PFI, ALE, PDP |  GLM, RF |  |  ✓ |  |  |  |  |
|  Joint Training |  [150] Prototype |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Local Approximation |  [137] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [136] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [138] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [223] SHAP |  RF |  |  ✓ |  |  |  |  |
|  |  [139] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [140] SHAP |  GLM, SVM, FNN, RF, kNN, GB |  |  ✓ |  |  ✓ |  |  |
|  Model Translation, Feature Selection |  [117] MDI, Rule Extraction |  RF, DT |  |  ✓ |  |  ✓ |  |  |
|  Perturbation |  [202] Spatial Variable Importance Profiles†PFI |  LDA, kNN LDA, RF |  |  ✓ |  |  |  ✓ |  |
|  |  [141] PFI |  RF |  |  ✓ |  |  ✓ |  |  |
|  |  [151] |  GP |  |  ✓ |  |  ✓ |  |  |
|  Perturbation, Embedding Space |  [143] Occlusion, Activation Assessment |  LSTM |  |  ✓ |  |  |  |  |  ✓
|  |  [146] Occlusion, Attention |  Transformer |  |  |  ✓ |  ✓ |  |  |
Atmosphere Monitoring |  Backpropagation |  [224] IG |  NN Ensemble |  |  ✓ |  |  |  |  |
|  Backpropagation, Local Approximation, Perturbation |  [225] SHAP, XRAI, Occlusion |  CNN |  |  ✓ |  |  |  |  |
|  Embedding Space |  [226] Activation Assessment |  FNN |  |  ✓ |  |  |  |  |
|  |  [227] Attention |  FNN |  |  ✓ |  |  |  |  |
|  Embedding Space, Joint Training, Feature Selection |  [122] Model Association, MDI, Activation Assessment |  RF, NN+ML Ensemble, FNN |  |  ✓ |  |  ✓ |  |  |
|  Feature Selection |  [123] MDI |  Tree Ensemble |  |  ✓ |  |  ✓ |  |  |
|  |  [124] MDI |  Tree Ensemble |  |  ✓ |  |  ✓ |  |  |
|  |  [125] MDI |  Tree Ensemble |  |  ✓ |  |  ✓ |  |  |
|  Interpretable by Design |  [228] GLM†GLM |  GLM |  |  ✓ |  |  |  ✓ |  |
|  Local Approximation |  [229] SHAP |  FNN |  |  ✓ |  |  |  |  |
|  |  [230] SHAP |  GLM, GB |  |  ✓ |  |  |  |  |
|  |  [231] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [232] SHAP |  GB |  |  ✓ |  |  |  |  |
|  Local Approximation, Feature Selection |  [115] MDI, SHAP |  RF |  |  ✓ |  |  ✓ |  |  |
|  Perturbation |  [233] PFI |  aCNN |  |  ✓ |  |  ✓ |  |  |
|  |  [234] PFI, ALE |  RF |  |  ✓ |  |  ✓ |  ✓ |  |
|  Perturbation, Feature Selection |  [120] MDI, PDP |  RF |  |  ✓ |  |  |  |  |
|  Perturbation, Local Approximation |  [235] SHAP, LIME, Ceteris Paribus Profiles, PFI |  RF |  |  ✓ |  |  |  |  |
|  |  [236] SHAP, PDP |  GB |  |  ✓ |  |  ✓ |  |  |
Building Mapping |  Backpropagation |  [237] Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [238] A*G, DeepLift |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [239] Grad-CAM |  CoAtNet |  |  ✓ |  |  ✓ |  |  |
|  |  [211] CAM, Score-CAM, Grad-CAM++, SmoothGrad-CAM++, Grad-CAM |  CNN |  |  |  ✓ |  ✓ |  |  |
|  Backpropagation, Embedding Space |  [240] Attention, Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Backpropagation, Embedding Space, Perturbation |  [222] IG, Occlusion, Activation Assessment |  Transformer, CNN |  |  ✓ |  |  ✓ |  ✓ |  |
|  Local Approximation |  [205] SHAP |  GB |  |  ✓ |  |  |  |  |
|  Local Approximation, Feature Selection |  [114] MDI, SHAP |  RF, FNN |  |  ✓ |  |  |  ✓ |  |
Ecosystem Interactions |  Interpretable by Design |  [180] GAM |  GAM |  |  ✓ |  |  |  |  |
|  Interpretable by Design, Local Approximation |  [178] GLM, SHAP |  GLM, GB |  |  ✓ |  |  |  |  |
|  Local Approximation |  [241] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [242] SHAP |  GB |  |  ✓ |  |  |  |  |
|  Perturbation, Local Approximation |  [243] SHAP, PFI |  GB |  |  ✓ |  |  ✓ |  |  |
Human Environment Interaction |  Backpropagation |  [204] GuidedGrad-CAM |  CNN |  |  ✓ |  |  |  |  |
|  |  [220] LRP |  CNN |  |  ✓ |  |  ✓ |  |  |  ✓
|  Backpropagation, Embedding Space, Perturbation |  [210] Occlusion, Grad-CAM, Activation Assessment |  GLM, NN Ensemble |  |  |  ✓ |  |  |  |  ✓
|  Interpretable by Design |  [244] GLM |  GLM |  |  ✓ |  |  ✓ |  |  |
|  Joint Training |  [177] Explanation Association |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [176] Explanation Association |  CNN |  |  ✓ |  |  |  |  |
|  Local Approximation |  [245] SHAP |  RF |  |  ✓ |  |  |  |  |
|  |  [246] SHAP |  GB |  |  ✓ |  |  |  |  |
|  |  [247] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [248] SHAP |  GB |  |  ✓ |  |  |  |  |
|  |  [249] SHAP, LIME |  RF |  |  ✓ |  |  |  |  |
|  Perturbation |  [250] ALE |  RF |  |  ✓ |  |  |  |  |
|  |  [251] PFI |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Perturbation, Embedding Space |  [215] Occlusion, Activation Assessment |  CNN |  |  ✓ |  |  |  |  |
|  Perturbation, Local Approximation |  [252] SHAP, LIME, Occlusion |  FNN |  |  ✓ |  |  |  |  |
Hydrology Monitoring |  Backpropagation |  [206] Gradient |  GP |  ✓ |  ✓ |  |  |  ✓ |  |
|  |  [253] Gradient |  GP |  |  ✓ |  |  ✓ |  |  |
|  |  [254] Grad-CAM |  CNN |  |  ✓ |  ✓ |  ✓ |  ✓ |  |
|  |  [255] Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [256] Gradient |  CNN |  |  ✓ |  |  |  |  |
|  |  [257] Grad-CAM |  CNN |  |  ✓ |  |  |  ✓ |  |
|  Backpropagation, Local Approximation |  [258] SHAP, IG |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Backpropagation, Perturbation |  [259] PDP, IG, DeepLift, Expected Gradients |  LSTM |  |  ✓ |  |  ✓ |  |  |
|  Embedding Space |  [260] Attention |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [207] smlp†Attention |  FNN |  ✓ |  ✓ |  |  |  ✓ |  |
|  Interpretable by Design, Local Approximation, Perturbation |  [181] Cubist, ALE, LIME, PFI |  Cubist, GB |  |  ✓ |  |  |  |  |
|  Local Approximation |  [261] SHAP |  FNN |  |  ✓ |  |  ✓ |  |  |
|  |  [262] SHAP |  RF, GB, DT, GB |  |  ✓ |  |  ✓ |  |  |
Landcover Mapping |  Backpropagation |  [197] 3DGrad-CAM†Grad-CAM |  CNN |  |  ✓ |  ✓ |  ✓ |  ✓ |  |
|  |  [218] Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [192] CSG-CAM†Grad-CAM, Score-CAM, Grad-CAM++, Grad-CAM |  CNN |  |  ✓ |  ✓ |  |  ✓ |  |
|  |  [201] ERC-CAM†CAM |  CNN |  |  ✓ |  |  |  ✓ |  |
|  |  [263] XRAI |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [264] Grad-CAM |  CNN |  |  ✓ |  ✓ |  |  ✓ |  |
|  |  [126] Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [198] MPGrad-CAM†Grad-CAM, Grad-CAM++, SmoothGrad-CAM++, Grad-CAM |  CNN |  |  ✓ |  ✓ |  |  ✓ |  |
|  |  [127] CAM |  CNN |  |  ✓ |  ✓ |  ✓ |  |  |
|  Backpropagation, Embedding Space |  [209] CAM, Attention, Grad-CAM, Layer-CAM, Attention, Attention |  CNN |  |  |  ✓ |  |  ✓ |  |
|  Backpropagation, Local Approximation, Perturbation |  [208] IG, DeepLift, Gradient, GuidedBackprop, LIME, GuidedGrad-CAM, Grad-CAM, I*G, Occlusion |  CNN |  |  ✓ |  ✓ |  ✓ |  |  |
|  Counterfactuals |  [133] GAN |  CNN |  |  |  ✓ |  ✓ |  |  |
|  Embedding Space |  [265] Activation Assessment |  CNN, CapsuleNet |  |  ✓ |  |  ✓ |  |  |
|  |  [266] Activation Assessment |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [267] Attention |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [268] Attention, Activation Assessment |  CNN |  |  ✓ |  |  |  ✓ |  |
|  Embedding Space, Local Approximation |  [189] SHAP, Activation Assessment |  NN Ensemble |  |  ✓ |  |  ✓ |  |  |
|  |  [129] LIME, Activation Assessment |  CNN |  |  ✓ |  |  |  ✓ |  |
|  |  [187] SHAP, Activation Assessment |  NN Ensemble |  |  ✓ |  |  ✓ |  ✓ |  |
|  |  [188] SHAP, Activation Assessment |  NN Ensemble |  |  ✓ |  |  ✓ |  |  |
|  Example-based |  [134] WIK |  CNN |  |  ✓ |  |  ✓ |  ✓ |  |  ✓
|  Interpretable by Design |  [269] Rules |  NN+Tree Ensemble |  |  ✓ |  |  ✓ |  |  |
|  |  [183] LDA |  LDA |  |  ✓ |  |  ✓ |  |  |
|  Interpretable by Design, Embedding Space |  [270] LDA, Activation Assessment |  NN Ensemble |  |  ✓ |  |  |  ✓ |  |
|  Joint Training |  [271] Prototype |  CNN |  |  ✓ |  |  |  ✓ |  |
|  Local Approximation |  [131] SHAP |  RF |  |  ✓ |  |  ✓ |  |  |
|  |  [130] SHAP, LIME |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [217] SHAP |  SVM, RF, GB, CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [128] SHAP |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [272] LIME |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Local Approximation, Feature Selection |  [111] MDI, SHAP |  RF, CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [112] MDI, LIME |  RF |  |  ✓ |  |  ✓ |  |  |
|  Perturbation, Local Approximation |  [132] SHAP, ALE, PDP, PFI |  RF |  |  ✓ |  |  |  ✓ |  |
|  Perturbation, Local Approximation, Feature Selection |  [110] MDI, SHAP, PFI |  RF |  |  ✓ |  |  ✓ |  |  |
Natural Hazard Monitoring |  Backpropagation, Embedding Space |  [164] DeepLift, Attention |  NN Ensemble |  |  ✓ |  |  |  ✓ |  |
|  Backpropagation, Embedding Space, Perturbation |  [221] Occlusion, Grad-CAM, Activation Assessment |  CNN |  |  ✓ |  |  ✓ |  ✓ |  |
|  Backpropagation, Local Approximation, Perturbation |  [213] SHAP, IG, PDP |  LSTM |  |  ✓ |  |  |  |  |
|  Embedding Space, Local Approximation |  [162] SHAP, Activation Assessment |  GAMI-Net |  |  ✓ |  |  ✓ |  |  |
|  Interpretable by Design |  [161] EBM |  EBM |  |  ✓ |  |  ✓ |  |  |
|  |  [184] Fuzzy Rules |  Rules |  |  ✓ |  |  ✓ |  |  |  ✓
|  |  [163] SNN†FNN |  GAM |  |  ✓ |  |  |  ✓ |  |
|  Joint Training |  [168] Prototype |  CNN |  |  ✓ |  |  |  ✓ |  |
|  Local Approximation |  [169] SHAP |  FNN |  |  ✓ |  |  |  |  |
|  |  [155] SHAP |  SVM, RF |  |  ✓ |  |  ✓ |  |  |
|  |  [165] SHAP |  RF, GB |  |  ✓ |  |  ✓ |  |  |
|  |  [157] SHAP |  FNN |  |  ✓ |  |  ✓ |  |  |
|  |  [273] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [158] SHAP |  GB |  |  ✓ |  |  |  ✓ |  |
|  |  [175] SHAP |  NN+Tree Ensemble |  |  ✓ |  |  ✓ |  |  |
|  |  [166] SHAP |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [172] SHAP, LIME |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [154] SHAP |  RF |  |  ✓ |  |  ✓ |  |  |
|  |  [159] SHAP |  GB |  |  ✓ |  |  |  |  |
|  Perturbation |  [171] PFI, PDP |  RF |  |  ✓ |  |  |  |  |
|  Perturbation, Feature Selection |  [119] MDI, PFI |  GLM, SVM, FNN, RF, kNN, GLM, Tree Ensemble, NB |  |  ✓ |  |  ✓ |  |  |
|  Perturbation, Local Approximation |  [156] SHAP, PDP, PFI |  CNN, NN Ensemble, FNN |  |  ✓ |  |  ✓ |  |  |
|  |  [174] LIME, Occlusion, PDP |  Tree Ensemble |  |  ✓ |  |  ✓ |  |  |
|  |  [170] SHAP, PFI |  RF |  |  ✓ |  |  ✓ |  |  |
|  |  [160] LIME, PDP |  RF |  |  ✓ |  |  ✓ |  |  |
|  |  [167] SHAP, PDP |  GB |  |  ✓ |  |  ✓ |  |  |
Soil Monitoring |  Backpropagation, Perturbation |  [274] SquareGrad, VarGrad, Gradient, PFI, SmoothGrad, I*G, IG |  ConvLSTM |  |  ✓ |  |  ✓ |  |  |
|  Local Approximation |  [275] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [276] SHAP |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [277] SHAP |  SVM |  |  ✓ |  |  ✓ |  |  |
|  Perturbation |  [278] PFI |  RF |  |  ✓ |  |  ✓ |  |  |
Surface Temperature Prediction |  Backpropagation |  [279] LRP |  FNN |  |  ✓ |  |  ✓ |  |  |
|  Local Approximation |  [280] SHAP |  GB |  |  ✓ |  |  |  |  |
|  |  [281] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [282] SHAP |  GB |  |  ✓ |  |  |  |  |
|  |  [283] SHAP |  GB |  |  ✓ |  |  |  |  |
|  Local Approximation, Feature Selection |  [115] MDI, SHAP |  RF |  |  ✓ |  |  ✓ |  |  |
Target Mapping |  Backpropagation |  [194] Self-Matching-CAM†Grad-CAM |  CNN |  |  |  ✓ |  |  ✓ |  |
|  |  [284] CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [190] CAM |  CNN |  |  |  |  |  ✓ |  |
|  |  [285] IG, Score-CAM |  NN Ensemble |  |  ✓ |  |  ✓ |  |  |
|  |  [286] Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Backpropagation, Embedding Space |  [287] CAM, Activation Assessment |  BagNet |  |  ✓ |  |  ✓ |  |  |
|  Backpropagation, Embedding Space, Perturbation |  [191] Attention, Empirical Receptive Field, Grad-CAM, Activation Assessment |  CNN |  |  ✓ |  |  ✓ |  ✓ |  |
|  Backpropagation, Local Approximation |  [288] GuidedBackprop, SHAP |  CNN |  |  ✓ |  |  |  ✓ |  |
|  |  [289] CAM, LIME |  CNN |  |  ✓ |  |  |  ✓ |  |
|  Embedding Space |  [290] Attention |  CNN |  |  |  |  |  ✓ |  |
|  |  [291] Activation Assessment |  VAE |  |  ✓ |  |  ✓ |  |  |
|  |  [292] Activation Assessment |  BagNet |  |  ✓ |  |  ✓ |  |  |
|  |  [293] Attention |  aCNN |  |  ✓ |  |  ✓ |  |  |
|  Local Approximation |  [294] SHAP |  CNN |  |  ✓ |  |  |  ✓ |  |
|  |  [295] LIME |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [296] LIME |  CNN |  |  ✓ |  |  |  ✓ |  |
|  Perturbation |  [297] Occlusion |  NN Ensemble |  |  ✓ |  |  ✓ |  |  |
Vegetation Monitoring |  Backpropagation |  [196] Eigen-CAM, AugScore-CAM, Crown-CAM†Grad-CAM, Score-CAM |  CNN |  |  |  ✓ |  |  ✓ |  |
|  |  [298] GuidedGrad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Backpropagation, Local Approximation |  [299] Deconvolution, Gradient, SHAP, GuidedGrad-CAM |  CNN |  |  |  ✓ |  |  ✓ |  |
|  Interpretable by Design |  [179] GLM |  GLM |  |  ✓ |  |  |  |  |
|  |  [185] Fuzzy Rules |  FLS |  |  ✓ |  |  ✓ |  |  |
|  Joint Training |  [300] Explanation Association |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Local Approximation |  [212] SHAP |  FNN |  |  |  ✓ |  ✓ |  ✓ |  |
|  |  [203] SHAP |  SVM, RF, GB |  |  ✓ |  |  ✓ |  |  |
|  |  [301] SHAP |  Extra Tree |  |  ✓ |  |  ✓ |  |  |
|  |  [302] SHAP |  RF |  |  ✓ |  |  |  |  |
|  Local Approximation, Feature Selection |  [109] MDI, LIME |  GB |  |  ✓ |  |  ✓ |  |  |
|  Model Translation |  [303] Rule Extraction |  RF |  |  ✓ |  |  |  |  |
|  |  [304] Rule Extraction |  RF |  |  ✓ |  |  |  ✓ |  |
|  Perturbation |  [305] PFI, PDP |  RF |  |  ✓ |  |  |  |  |
|  Perturbation, Local Approximation |  [214] SHAP, PFI |  GB |  |  ✓ |  |  |  |  |
|  |  [306] SHAP, PFI |  RF |  |  ✓ |  |  |  |  |
Weather Climate Prediction |  Backpropagation |  [307] LRP |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [308] LRP |  FNN |  |  ✓ |  |  ✓ |  |  |
|  |  [309] LRP |  FNN |  |  ✓ |  |  |  |  |
|  |  [219] Grad-CAM |  CNN |  |  ✓ |  |  ✓ |  |  |
|  |  [310] GuidedBackprop |  CNN |  |  ✓ |  |  |  ✓ |  |
|  Backpropagation, Feature Selection |  [113] MDI, Grad-CAM |  RF, CNN |  |  ✓ |  |  ✓ |  |  |
|  Feature Selection |  [311] TreeInterpreter |  RF |  |  ✓ |  |  |  |  |
|  Joint Training |  [200] Prototype†ProtoPNet |  CNN |  ✓ |  ✓ |  |  |  ✓ |  |
|  Local Approximation |  [312] SHAP |  GB |  |  ✓ |  |  ✓ |  |  |
|  |  [216] SHAP |  GB |  |  ✓ |  |  |  ✓ |  |
|  Perturbation, Local Approximation |  [313] SHAP, PFI |  GLM, RF, SVM, FNN |  |  ✓ |  |  ✓ |  |  |
Other |  Backpropagation |  [314] Gradient |  FNN |  |  ✓ |  |  ✓ |  |  |
|  Feature Selection |  [116] MDI |  RF |  |  ✓ |  |  ✓ |  |  |
|  Interpretable by Design |  [182] LDA |  LDA |  |  ✓ |  |  |  |  |
|  Interpretable by Design, Embedding Space |  [270] LDA, Activation Assessment |  NN Ensemble |  |  ✓ |  |  |  ✓ |  |
|  Joint Training |  [315] Explanation Association |  CNN |  |  ✓ |  |  ✓ |  |  |
|  Local Approximation |  [316] SHAP |  RF |  |  ✓ |  |  ✓ |  |  |
|  |  [317] SHAP |  RF |  |  ✓ |  |  ✓ |  |  |
|  |  [318] SHAP |  GB |  |  ✓ |  |  |  |  |
|  Perturbation |  [319] PFI, PDP |  RF |  |  ✓ |  |  |  |  |

##  VI Discussion

In this section, we address research questions RQ4 and RQ5 and discuss the
usage of xAI methods in RS, the adaptation of xAI methods to RS problems, as
well as the evaluation of such methods (RQ4). Furthermore, we highlight the
challenges, limitations, and emerging research directions in the field (RQ5).

xAI in RS is still a very young and dynamic field. Figure
LABEL:fig:yearly_eo_xaif should give a broad overview of the development and
evolution of the field in the last years and provide the reader with the
context for the usage as well as the faced challenges of xAI in RS. Figure
LABEL:fig:yearly_eo_xaifa illustrates varying trends in the EO tasks, with
vegetation, atmosphere, or natural hazard monitoring recently getting more
attention. Meanwhile, agricultural monitoring and target mapping show a
constant trend, while landcover mapping publications are decreasing. This
suggests a varying time lag before new approaches from ML are widely applied
in RS. Indeed, landcover mapping is one of the most established EO tasks, with
a variety of benchmark datasets and models. While many approaches have already
been applied to landcover mapping, other tasks are currently being explored,
creating a trend in these areas. Similarly, Figure LABEL:fig:yearly_eo_xaifb
reveals increasing trends for local approximation and perturbation methods
while backpropagation methods stagnate. This stagnation may be attributed to
an interconnection between EO tasks. Initial assessments of xAI often involve
simple and operable methods, specifically, local approximation and
perturbation methods, which are model-agnostic. Hence, reflecting the high
usage of xAI in RS while these xAI methods are extensively explored.

###  VI-A RQ4: Recommended practices in xAI for EO

####  VI-A1 The usage of xAI

In order to get meaningful and reliable explanations, xAI methods should be
only applied to ML models that have achieved a good generalization [320].
Furthermore, the robustness of ML models and xAI methods should be quantified,
as it is common for them to disagree [321]. To maintain consistent outcomes,
it is advisable to consider a set of initializations of xAI methods and
different ML models, as it is done when benchmarking the performance of a ML
model. For instance, in the works identified in our review, different model
types [165, 203, 156, 254, 155, 121, 259], model seeds or configurations [145,
259], and xAI method seeds [233] were used. Furthermore, the outcomes of
different xAI methods can be compared, as it is done in [155, 181, 192, 238,
259, 209, 282, 196, 144, 142, 117, 198, 211, 178]. For example, [259] apply
IG, Expected Gradients (EG), and DeepLift to exhaustively interpret different
RNN-based architectures (RNN, LSTM, GRU), also varying their numbers of
neurons [259]. Other related works are [178], where [178] compare the results