# DATA FUSION OF SATELLITE IMAGERY FOR GENERATION OF DAILY CLOUD-FREE IMAGES
AT HIGH RESOLUTION LEVEL

Natalya Ivanchuk<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Computer Sciences and Applied Mathematics, National University
of Water and Environmental Engineering, Soborna str., 11, Rivne, Ukraine,
33028; EOS Data Analytics Ukraine, Desyatynny lane, 5, 01001 Kyiv, Ukraine
Peter Kogut<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Differential
Equations, Oles Honchar Dnipro National University, Gagarin av., 72, 49010
Dnipro, Ukraine; EOS Data Analytics Ukraine, Gagarin av., 103a, Dnipro,
Ukraine Petro Martyniuk<EMAIL_ADDRESS><EMAIL_ADDRESS>) Institute of Automation, Cybernetics and Computer Engineering, National
University of Water and Environmental Engineering, Soborna str., 11, Rivne,
Ukraine, 33028; EOS Data Analytics Ukraine, Gagarin av., 103a, Dnipro, Ukraine

###### Abstract

In this paper we discuss a new variational approach to the Date Fusion problem
of multi-spectral satellite images from Sentinel-2 and MODIS that have been
captureed at different resolution level and, arguably, on different days. The
crucial point of our approach that the MODIS image is cloud-free whereas the
images from Sentinel-2 can be corrupted by clouds or noise.

###### keywords:

Data fusion, variational approach, Sentinel-2, Moderate Resolution Imaging
Spectroradiometer (MODIS), time series, image restoration, data assimilation,
constrained minimization problems.

###### AMS:

94A08, 49Q20, 49K20, 49J45.

##  1 Introduction

Following [23], the Image Fusion is a process of combining the relevant
information from a set of images of the same scene into a single image and the
resultant fused image must be more informative and complete than any of the
input images. At the same time, when we deal with the data fusion problem for
satellite images, such images, as a rule, are multi sensor, multi-modal,
multi-focus and multi temporal. Moreover, the data fusion problem is often
exacerbated by cloud contamination. In some cloudy areas, researchers are
fortunate to get 2–3 cloud-free satellite scenes per year, what is
insufficient for many applications that require dense temporal information,
such as crop condition monitoring and phenology studies [19, 39]. In view of
this, we can indicate the following general requirements for the satellite
image fusion process:

  * (i)

The fused image should preserve all relevant information from the input
images;

  * (ii)

The image fusion should not introduce artifacts which can lead to wrong
inferences.

In spite of the fact that the first requirement (item (i)) sounds rather
vague, we give a precise treatment for it in Section 5, making use of a
collection of special constrained minimization problems (see (37)–(38)). As
for the second item, it is important to emphasize that we are mainly
interesting by satellite images that can be useful from agricultural point of
view (land cover change mapping, crop condition monitoring, yield estimation,
and many others). Because of this an important option in the image data fusion
is to preserve the precise geo-location of the existing crop fields and avoid
an appearance of the so-called false contours and pseudo-boundaries on a given
territory.

In this paper we mainly focus on the image fusion problem coming from two
satellites — Satnitel-2 and Moderate Resolution Imaging Spectroradiometer
(MODIS). Since each band (spectral channel) in Sentinel images has $10$, $20$,
or $60$ meters in pixel size, it gives an ideal spatial resolution for
vegetation mapping at the field scale. Moreover, taking into account that
Sentinel-2 has $3$–$5$ revisit cycle over the same territory, it makes its
usage for studying global biophysical processes, which allows to evolve
rapidly during the growing season, essentially important. The unique problem
that drastically restricts its practical implementation, is the fact that the
satellite images, as a rule, are often contaminated by clouds, shadows, dust,
and other atmospheric artifacts.

One of possible solutions for practical applications is to make use of
frequent coarse-resolution data of the MODIS. Taking into account that the
MODIS data can be delivered with the daily repeat cycle and 500-m surface
reflectance, the core idea is to use the Sentinel and cloud-free MODIS data to
generate synthetic ’daily’ surface reflectance products at Sentinel spatial
resolution [31].

The problem we consider in this paper can be briefly described as follows. We
have a collection of multi-band images
$\left\\{S_{1},S_{2},\dots,S_{N}:G_{H}\to\mathbb{R}^{m}\right\\}$ from
Sentinel-2 that were captured at some time instances
$\left\\{t_{1},t_{2},\dots,t_{N}\right\\}$, respectively, and we have a MODIS
image $M:G_{L}\to\mathbb{R}^{n}$ from some day $t_{M}$. It is assumed that all
of these images are well co-registered with respect to the unique geographic
location. We also suppose that the MODIS image is cloud-free and the day
$t_{M}$ may does not coincide with any of time instances
$\left\\{t_{1},t_{2},\dots,t_{N}\right\\}$. Meanwhile, the Sentinel images
$\left\\{S_{2},\dots,S_{N}:G_{H}\to\mathbb{R}^{m}\right\\}$ can be corrupted
by clouds. The main question is how to generate a new synthetic ’daily’ multi-
band image of the same territory from the day $t_{M}$ at the Sentinel-2
spatial resolution $G_{H}$, utilizing for that the above mentioned data.

In principle, this problem is not new in the literature (see, for instance,
[19, 20, 24, 34, 39]). For nowadays the spatial and temporal adaptive
reflectance fusion model (STARFM) is one of the most popular model where the
idea to generate a new synthetic ’daily’ satellite images at high resolution
level has been realized. Basing on a deterministic weighting function computed
by spectral similarity, temporal difference, and spatial distance, this model
(as many other generalizations) allows to predict daily surface reflectance at
high spatial resolution and MODIS temporal frequency. However, its performance
essentially depends on the characteristic patch size of the landscape and
degrades somewhat when used on extremely heterogeneous fine-grained landscapes
[19].

As it was mentioned in [5], the majority image interpolation techniques are
essentially based on certain assumptions regarding the data and the
corresponding imaging devices. First, it is often assumed that the image of
high spatial resolution is a linear combination of the spectral channels with
known weights [33]. Second, the loss of resolution is usually modeled as a
linear operator which consists of a subsampled convolution with known kernel
(point spread function). While both assumptions may be justified in some
applications, it may be difficult to measure or estimate the weights and the
convolution kernel in a practical situation.

Instead of this, we mainly focus on the variational approach to the satellite
image data fusion. In order to eliminate the above mentioned restrictions, we
formulate the data fusion problem as the two-level optimization problem. At
the first level, following a simple iterative procedure, we generate the so-
called structural prototype for a synthetic Sentinel image from the given day
$t_{M}$. The main characteristic feature of this prototype is the fact that,
it must have a similar geometrical structure (namely, precise location of
contours and field boundaries) to the nearest in time ’visible’ Sentinel
images, albeit they may have rather different intensities in all bands. Since
the revisit cycle of Sentinel-2 is $3$–$5$ days, such prototype can be easily
generated (for the details, we refer to Section 4). In fact, we consider the
above mentioned structural prototype as a reasonable input data for ’daily’
prediction problem that we formulate in the form of a special constrained
minimization problem, where the cost functional has a nonstandard growth and
the edge information for restoration of MODIS cloud-free images at the
Sentinel resolution is accumulated both in the variable exponent of
nonlinearity and in the directional image gradients which we derive from the
predicted structural prototype.

Our approach is based on the variational model in Sobolev-Orlicz space with a
non-standard growth condition of the objective functional and on the
assumption that, to a large extent, the image topology in each spectral
channel is similar to the topographic map of its structural prototype. It is
worth to emphasize that this model is considerably different from the
variational model for P+XS image fusion that was proposed in [2]. We discuss
the well thoroughness of the above approach and consistency of the
corresponding variational problem, and we show that this problem admits a
unique solution. We also derive some optimality conditions and supply our
approach by results of numerical simulations with the real satellite images.

The paper is organized as follows. Section 2 contains some preliminaries,
auxiliary results, and a non-formal statement of the data fusion problem. In
Section 3 we begin with some key assumptions and after we give a precise
statement of the satellite image data fusion in the form of two-level
constrained optimization problem with a nonstandard growth energy functional.
We show that, in principle, we can distinguish three different statements of
the data fusion problem. Namely, it is the so-called restoration problem (when
the main point is to restore the information in the cloud-corrupted zone for
Sentinel images), the interpolation problem (when the day $t_{M}$ is
intermediate for some time instances of
$\left\\{t_{1},t_{2},\dots,t_{N}\right\\}$), and the so-called extrapolation
problem (when $t_{M}>t_{N}$). We also discuss the specific of each of these
problems and their rigorous mathematical description. Section 4 is devoted to
the study of a particular model for the prediction of structural prototypes.
We also illustrate this approach by some numerical simulations.

Consistency issues of the proposed minimization problems, optimality
conditions, and their substantiation are studied in Sections 5 and 6. For
illustration of this approach, we give in Section 7 some results of numerical
experiments with real satellite images. The experiments undertaken in this
study confirmed the efficacy of the proposed method and revealed that it can
acquire plausible visual performance and satisfactory quantitative accuracy
for agro-scenes with rather complicated texture of background surface. In
Appendix we give the main auxiliary results concerning the Orlicz spaces and
the Sobolev-Orlicz spaces with variable exponent.

##  2 Non-Formal Statement of the Problem

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded connected open set with a
sufficiently smooth boundary $\partial\Omega$ and nonzero Lebesgue measure. In
majority cases $\Omega$ can be interpreted as a rectangle domain. Let $G_{H}$
and $G_{L}$ be two sample grids on $\Omega$ such that
$G_{H}=\widehat{G}_{H}\cap\Omega$ and $G_{H}=\widehat{G}_{H}\cap\Omega$, where

$\displaystyle\widehat{G}_{H}$
$\displaystyle=\left\\{(x_{i},y_{j})\left|\begin{array}[]{c}x_{1}=x_{H},\
x_{i}=x_{1}+\Delta_{H,x}(i-1),\ i=1,\dots,N_{x},\\\ y_{1}=y_{H},\
y_{j}=y_{1}+\Delta_{H,y}(j-1),\ j=1,\dots,N_{y},\end{array}\right.\right\\},$
$\displaystyle\widehat{G}_{L}$
$\displaystyle=\left\\{(x_{i},y_{j})\left|\begin{array}[]{c}x_{1}=x_{L},\
x_{i}=x_{1}+\Delta_{L,x}(i-1),\ i=1,\dots,M_{x},\\\ y_{1}=y_{L},\
y_{j}=y_{1}+\Delta_{L,y}(j-1),\ j=1,\dots,M_{y},\end{array}\right.\right\\},$

with some fixed points $(x_{H},y_{H})$ and $(x_{L},y_{L})$. Here, it is
assumed that $N_{x}>>M_{x}$ and $N_{y}>>M_{y}$.

Let $[0,T]$ be a given time interval. Normally, by $T$ we mean a number of
days. Let $t_{M}$ and $\left\\{t_{k}\right\\}_{k=1}^{N}$ be moments in time
(particular days) such that $0\leq t_{1}<t_{2}<\dots<t_{N}\leq T$ and
$t_{1}<t_{M}<T$. Let
$\left\\{S_{1},S_{2},\dots,S_{N}:G_{H}\to\mathbb{R}^{m}\right\\}$ be a
collection of multispectral images of some territory, delivered from
Sentinel-2, that were taken at time instances $t_{1},t_{2},\dots,t_{N}$,
respectively. Hereinaftre, $m=13$ and it stands for the number of spectral
channels in images from Sentinel-2. Let $M:G_{L}\to\mathbb{R}^{n}$, with
$n=6$, be a MODIS image of the same territory and this image has been captured
at time $t=t_{M}$. It is assumed that:

  1. (i)

The Sentinel-2 images $S_{k}:G_{H}\to\mathbb{R}^{m}$, $k=2,\dots,N$ can be
corrupted by some noise, clouds and blur. However, the first one
$S_{1}:G_{H}\to\mathbb{R}^{m}$ is a cloud-free image;

  2. (ii)

For further convenience we divide the set of all bands for Sentinel images
onto two parts $J_{1}$ and $J_{2}$ with $\mathrm{dim}(J_{1})=6$ and
$\mathrm{dim}(J_{2})=7$, where

Table 1: Spectral channels of $J_{1}$ Band | Resolution | Central Wavelenth | Description
---|---|---|---
$B_{2}$ |  $10$ m |  $490$ m | Blue
$B_{3}$ |  $10$ m |  $560$ m | Green
$B_{4}$ |  $10$ m |  $665$ m | Red
$B_{8a}$ |  $20$ m |  $865$ m | Visible and Near Infrared (VNIR)
$B_{11}$ |  $20$ m |  $1610$ m | Short Wave Infrared (SWIR)
$B_{12}$ |  $20$ m |  $2190$ m | Short Wave Infrared (SWIR)
Table 2: Spectral channels of $J_{2}$ Band | Resolution | Central Wavelenth | Description
---|---|---|---
$B_{1}$ |  $60$ m |  $443$ m | Ultra Blue (Coastal and Aerosol)
$B_{5}$ |  $20$ m |  $705$ m | Visible and Near Infrared (VNIR)
$B_{6}$ |  $20$ m |  $740$ m | Visible and Near Infrared (VNIR)
$B_{7}$ |  $20$ m |  $783$ m | Visible and Near Infrared (VNIR)
$B_{8}$ |  $10$ m |  $842$ m | Visible and Near Infrared (VNIR)
$B_{9}$ |  $60$ m |  $940$ m | Short Wave Infrared (SWIR)
$B_{10}$ |  $60$ m |  $1375$ m | Short Wave Infrared (SWIR)
  3. (iii)

Each spectral channel of the MODIS image
$M=\left[M_{1},M_{2},\dots,M_{6}\right]:G_{L}\to\mathbb{R}^{6}$ has the
similar spectral characteristics to the corresponding channel of $J_{1}$-group
$\left\\{B_{2},B_{3},B_{4},B_{8a},B_{11},B_{12}\right\\}$, respectively;

  4. (iv)

The principle point is that the MODIS image $M:G_{L}\to\mathbb{R}^{6}$ is
visually sufficiently clear and does not corrupted by clouds or its damage
zone can be neglected;

  5. (v)

The MODIS image $M:G_{L}\to\mathbb{R}^{6}$ and the images
$\left\\{S_{1},S_{2},\dots,S_{N}:G_{H}\to\mathbb{R}^{m}\right\\}$ from
Sentinel-2 are rigidly co-registered. This means that the MODIS image after
arguably some affine transformation and each Sentinel images after the
resampling to the grid with low resolution $G_{L}$, could be successfully
matched according to the unique geographic location.

In practice, the co-registration procedure can be realized using, for
instance, the open-source LSReg v2.0.2 software [37, 40] that has been used in
a number of recent studies [18, 38], or the rigid co-registration approach
that has been recently developed in [27, 28]. However, in both cases, in order
to find an appropriate affine transformation, we propose to apply this
procedure not to the original images, but rather to the contour’s map of their
spectral energies $Y_{M}:G_{L}\to\mathbb{R}$ and
$Y_{S_{j}}:G_{H}\to\mathbb{R}$, where the last ones should be previously
resampled to the grid of the low resolution $G_{L}$. Here,

$\displaystyle Y_{M}(z)$
$\displaystyle:=\alpha_{1}M_{1}(z)+\alpha_{2}M_{2}(z)+\alpha_{3}M_{3}(z),\quad\forall\,z=(x,y)\in
G_{L},$ $\displaystyle Y_{S_{i}}(z)$
$\displaystyle:=\alpha_{1}S_{i,1}(z)+\alpha_{2}S_{i,2}(z)+\alpha_{3}S_{i,3}(z),\quad\forall\,z=(x,y)\in
G_{H}$

with $\alpha_{1}=0.114$, $\alpha_{2}=0.587$, and $\alpha_{3}=0.299$.

######  Remark 2.1.

Let us mention that in the case of digital images, the only accessible
information is a sampled and quantized version of $I:\Omega\to\mathbb{R}^{m}$,
$I(x_{i},y_{j})$, where $\left\\{(x_{i},y_{j})\in\Omega\right\\}$ is a set of
discrete points and for each spectral channel $k=1,\dots,m$,
$I_{k}(x_{i},y_{j})$ belongs in fact to a discrete set of values,
$0,1,\dots,255$ in many cases. Due to the Shannon’s theory, it is plausible to
assume that $I_{k}$ is recoverable at any point $(x,y)\in\Omega$ from the
samples $I_{k}(x_{i},y_{j})$. So, we may assume that the image $I$ is known in
a continuous domain, up to the quantization noise (see [2]). However, in
practice, such reconstruction is not a trivial problem.

###  2.1 Functional Spaces

Let us recall some useful notations. For vectors $\xi\in\mathbb{R}^{2}$ and
$\eta\in\mathbb{R}^{2}$, $\left(\xi,\eta\right)=\xi^{t}\eta$ denotes the
standard vector inner product in $\mathbb{R}^{2}$, where t stands for the
transpose operator. The norm $|\xi|$ is the Euclidean norm given by
$|\xi|=\sqrt{(\xi,\xi)}$. Let $\Omega\subset\mathbb{R}^{2}$ be a bounded open
set with a Lipschitz boundary $\partial\Omega$. For any subset
$E\subset\Omega$ we denote by $|E|$ its $2$-dimensional Lebesgue measure
$\mathcal{L}^{2}(E)$. Let $\overline{E}$ denote the closure of $E$, and
$\partial E$ stands for its boundary. We define the characteristic function
$\chi_{E}$ of $E$ by

$\chi_{E}(x):=\left\\{\begin{array}[]{ll}1,&\ \text{for }\ x\in E,\\\ 0,&\
\text{otherwise}.\end{array}\right.$

Let $X$ denote a real Banach space with norm $\|\cdot\|_{X}$, and let
$X^{\prime}$ be its dual. Let $\left<\cdot,\cdot\right>_{X^{\prime};X}$ be the
duality form on $X^{\prime}\times X$. By $\rightharpoonup$ and
$\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}$ we denote the weak and
weak∗ convergence in normed spaces, respectively.

For given $1\leq p\leq+\infty$, the space $L^{p}(\Omega;\mathbb{R}^{2})$ is
defined by

$L^{p}(\Omega;\mathbb{R}^{2})=\left\\{f:\Omega\rightarrow\mathbb{R}^{2}\ :\
\|f\|_{L^{p}(\Omega;\mathbb{R}^{2})}<+\infty\right\\},$

where $\|f\|_{L^{p}(\Omega;\mathbb{R}^{2})}=\left(\int_{\Omega}\lvert
f(x)\rvert^{p}\,dx\right)^{1/p}$ for $1\leq p<+\infty$. The inner product of
two functions $f$ and $g$ in $L^{p}(\Omega;\mathbb{R}^{2})$ with
$p\in[1,\infty)$ is given by

$\left(f,g\right)_{L^{p}(\Omega;\mathbb{R}^{2})}=\int_{\Omega}\left(f(x),g(x)\right)\,dx=\int_{\Omega}\sum_{k=1}^{2}f_{k}(x)g_{k}(x)\,dx.$

We denote by $C_{c}^{\infty}(\mathbb{R}^{2})$ a locally convex space of all
infinitely differentiable functions with compact support in $\mathbb{R}^{2}$.
We recall here some functional spaces that will be used throughout this paper.
We define the Banach space $H^{1}(\Omega)$ as the closure of
$C^{\infty}_{c}(\mathbb{R}^{2})$ with respect to the norm

$\|y\|_{H^{1}(\Omega)}=\left(\int_{\Omega}\left(y^{2}+\lvert\nabla
y\rvert^{2}\right)\,dx\right)^{1/2}.$

We denote by $\left(H^{1}(\Omega)\right)^{\prime}$ the dual space of
$H^{1}(\Omega)$.

Given a real Banach space $X$, we will denote by $C([0,T];X)$ the space of all
continuous functions from $[0,T]$ into $X$. We recall that a function
$u:[0,T]\to X$ is said to be Lebesgue measurable if there exists a sequence
$\left\\{u_{k}\right\\}_{k\in\mathbb{N}}$ of step functions (i.e.,
$u_{k}=\sum_{j=1}^{n_{k}}a_{j}^{k}\chi_{A_{j}^{k}}$ for a finite number
$n_{k}$ of Borel subsets $A_{j}^{k}\subset[0,T]$ and with $a_{j}^{k}\in X$)
converging to $u$ almost everywhere with respect to the Lebesgue measure in
$[0,T]$.

Then for $1\leq p<\infty$, $L^{p}(0,T;X)$ is the space of all measurable
functions $u:[0,T]\to X$ such that

$\|u\|_{L^{p}(0,T;X)}=\left(\int_{0}^{T}\|u(t)\|^{p}_{X}\,dt\right)^{\frac{1}{p}}<\infty,$

while $L^{\infty}(0,T;X)$ is the space of measurable functions such that

$\|u\|_{L^{\infty}(0,T;X)}=\sup_{t\in[0,T]}\|u(t)\|_{X}<\infty.$

The full presentation of this topic can be found in [16].

Let us recall that, for $1\leq p\leq\infty$, $L^{p}(0,T;X)$ is a Banach space.
Moreover, if $X$ is separable and $1\leq p<\infty$, then the dual space of
$L^{p}(0,T;X)$ can be identified with $L^{p^{\prime}}(0,T;X^{\prime})$.

For our purpose $X$ will mainly be either the Lebesgue space $L^{p}(\Omega)$
or $L^{p}(\Omega;\mathbb{R}^{2})$ or the Sobolev space $W^{1,p}(\Omega)$ with
$1\leq p<\infty$. Since, in this case, $X$ is separable, we have that
$L^{p}(0,T;L^{p}(\Omega))=L^{p}(Q_{T})$ is the ordinary Lebesgue space defined
in $Q_{T}=(0,T)\times\Omega$. As for the space
$L^{p}(0,T;W^{1,\alpha}(\Omega))$ with $1\leq\alpha,p<+\infty$, it consists of
all functions $u:[0,T]\times\Omega\to\mathbb{R}$ such that $u$ and $|\nabla
u|$ belongs to $L^{p}(0,T;L^{\alpha}(\Omega))$. Moreover,

$\left(\int_{0}^{T}\left[\left(\int_{\Omega}|u|^{\alpha}\,dx\right)^{\frac{p}{\alpha}}+\left(\int_{\Omega}|\nabla
u|^{\alpha}\,dx\right)^{\frac{p}{\alpha}}\right]\,dt\right)^{\frac{1}{p}}$

defines the norm in $L^{p}(0,T;W^{1,\alpha}(\Omega))$.

###  2.2 Topographic Maps and Geometry of Satillite Multispectral Images

Following the main principle of the Mathematical Morphology, a scalar image
$u:\Omega\rightarrow\mathbb{R}$ is a representative of an equivalence class of
images $v$ obtained from $u$ via a contrast change, i.e., $v=F(u)$, where $F$
is a continuous strictly increasing function. Under this assumption, a scalar
image can be characterized by its upper (or lower) level sets
$Z_{\lambda}(u)=\left\\{x\in\Omega\ :\ u(x)\geq\lambda\right\\}$ (resp.
$Z^{\prime}_{\lambda}(u)=\left\\{x\in\Omega\ :\ u(x)\leq\lambda\right\\}$).
Moreover, each image can be recovered from its level sets by the
reconstruction formula $u(x)=\sup\left\\{\lambda\ :\ x\in
Z_{\lambda}(u)\right\\}$. Thus, according to the Mathematical Morphology
Doctrine, the reliable information in the image contains in the level sets,
independently of their actual levels (see [7] for the details). So, we can
suppose that the entire geometric information about a scalar image is
contained in those level sets.

In order to describe the level sets by their boundaries, $\partial
Z_{\lambda}(u)$, we assume that $u\in W^{1,1}(\Omega)$, where
$W^{1,1}(\Omega)$ stands for the standard Sobolev space of all functions $u\in
L^{1}(\Omega)$ with respect to the norm

$\|u\|_{W^{1,1}(\Omega)}=\|u\|_{L^{1}(\Omega)}+\|\nabla
u\|_{L^{1}(\Omega)^{2}},$

where the distributional gradient $\nabla u=\left(\frac{\partial u}{\partial
x_{1}},\frac{\partial u}{\partial x_{2}}\right)$ is represented as follows

$\int_{\Omega}u\frac{\partial\phi}{\partial
x_{i}}\,dx=-\int_{\Omega}\phi\frac{\partial u}{\partial
x_{i}}\,dx,\quad\forall\,\phi\in C^{\infty}_{0}(\Omega),\ i=1,2.$

It was proven in [1] that if $u\in W^{1,1}(\Omega)$ then its upper level sets
$Z_{\lambda}(u)$ are sets of finite perimeter. So, the boundaries $\partial
Z_{\lambda}(u)$ of level sets can be described by a countable family of Jordan
curves with finite length, i.e., by continuous maps from the circle into the
plane $\mathbb{R}^{2}$ without crossing points. As a result, at almost all
points of almost all level sets of $u\in W^{1,1}(\Omega)$ we may define a unit
normal vector $\theta(x)$. This vector field formally satisfies the following
relations

$\left(\theta,\nabla u\right)=|\nabla u|\quad\text{and}\quad|\theta|\leq 1\
\text{a.e. in $\Omega$}.$

In the sequel, we will refer to the vector field $\theta$ as the vector field
of unit normals to the topographic map of a function $u$. So, we can associate
$\theta$ with the geometry of the scalar image $u$.

In the case of multi-band satellite images $I:\Omega\to\mathbb{R}^{m}$, we
will impose further the following assumption: $I\in
W^{1,1}(\Omega;\mathbb{R}^{m})$ and each spectral channel of a given image $I$
has the same geometry. We refer to [8] for the experimental discussion.

######  Remark 2.2.

In practice, at the discrete level, the vector field $\theta(x,y)$ can be
defined by the rule $\theta(x_{i},y_{j})=\frac{\nabla u(x_{i},y_{j})}{|\nabla
u(x_{i},y_{j})|}$ when $\nabla u(x_{i},y_{j})\neq 0$, and $\theta=0$ when
$\nabla u(x_{i},y_{j})=0$. However, as was mentioned in [2], a better choice
for $\theta(x,y)$ would be to compute it as the ration $\frac{\nabla
U(t,\cdot)}{|\nabla U(t,\cdot)|}$ for some small value of $t>0$, where
$U(t,x,y)$ is a solution of the following initial-boundary value problem with
$1D$-Laplace operator in the principle part

$\displaystyle\frac{\partial U}{\partial t}=\mathrm{div}\,\left(\frac{\nabla
U}{|\nabla U|}\right),\quad t\in(0,+\infty),\ (x,y)\in\Omega,$ (1)
$\displaystyle U(0,x,y)=u(x,y),\quad(x,y)\in\Omega,$ (2)
$\displaystyle\frac{\partial U(0,x,y)}{\partial\nu}=0,\quad t\in(0,+\infty),\
(x,y)\in\partial\Omega.$ (3)

As a result, for any $t>0$, there can be found a vector field

$\xi\in L^{\infty}(\Omega;\mathbb{R}^{2})\ \text{ with }\
\|\xi(t)\|_{L^{\infty}(\Omega;\mathbb{R}^{2})}\leq 1$

such that

$\left(\xi(t),\nabla U(t,\cdot)\right)=|\nabla U(t,\cdot)|\ \text{in }\
\Omega,\quad\xi(t)\cdot\nu=0\ \text{on }\ \partial\Omega,$ (4)

and $U_{t}(t,x,y)=\mathrm{div}\,\xi(t,x,y)$ in the sense of distributions on
$\Omega$ for a.a. $t>0$.

We notice that in the framework of this procedure, for small value of $t>0$,
we do not distort the geometry of the function $u(x,y)$ in an essential way.
Moreover, it can be shown that this regularization of the vector field
$\theta(x,y)=\frac{\nabla U(x,y)}{|\nabla U(x,y)|}$ satisfies condition
$\mathrm{div}\,\theta\in L^{2}(\Omega)$.

###  2.3 Texture Index of a Gray-Scale Image

Let $u\in C([0,T];L^{2}(\Omega))$ be a given function. For each $t\in[0,T]$,
we associate the real-valued mapping $u(t,\cdot):\Omega\mapsto\mathbb{R}$ with
a gray-scale image, and the mapping $u:(0,T)\times\Omega\rightarrow\mathbb{R}$
with an optical flow. A widely-used way to smooth $u(t,\cdot)$ is by
calculating the convolution

$u_{\sigma}(t,x):=\left(G_{\sigma}\ast\widetilde{u}(t,\cdot)\right)(x)=\int_{\mathbb{R}^{2}}G_{\sigma}(x-y)\widetilde{u}(t,y)\,dy,$

where $\widetilde{u}$ denotes zero extension of $u$ from
$Q_{T}=(0,T)\times\Omega$ to $\mathbb{R}^{3}$, and $G_{\sigma}$ stands for the
two-dimentional Gaussian of width (standard deviation) $\sigma>0$:

$G_{\sigma}(x)=\frac{1}{2\pi\sigma^{2}}\exp\left(-\frac{|x|^{2}}{2\sigma^{2}}\right).$

######  Definition 1.

We say that a function $p_{u}:(0,T)\times\Omega\to\mathbb{R}$ is the texture
index of a given optical flow $u:(0,T)\times\Omega\rightarrow\mathbb{R}$ if it
is defined by the rule

$p_{u}(t,x):=1+g\left(\frac{1}{h}\int_{t-h}^{t}\left|\left(\nabla
G_{\sigma}\ast\widetilde{u}(\tau,\cdot)\right)(x)\right|^{2}\,d\tau\right),\
\forall\,(t,x)\in Q_{T},$ (5)

where $g{:}[0,\infty)\rightarrow(0,\infty)$ is the edge-stopping function
which we take in the form of the Cauchy law $g(s)=\frac{a}{a+s}$ with $a>0$
small enough, and $h>0$ is a small positive value.

Since $G_{\sigma}\in C^{\infty}(\mathbb{R}^{2})$, it follows from (5) and
absolute continuity of the Lebesgue integral that $1<p_{u}(t,x)\leq 2$ in
$Q_{T}$ and $p_{u}\in C^{1}([0,T];C^{\infty}(\mathbb{R}^{2}))$ even if $u$ is
just an absolutely integrable function in $Q_{T}$. Moreover, for each
$t\in[0,T]$, $p_{u}(t,x)\approx 1$ in those places of $\Omega$ where some
edges or discontinuities are present in the image $u(t,\cdot)$, and
$p_{u}(t,x)\approx 2$ in places where $u(t,\cdot)$ is smooth or contains
homogeneous features. In view of this, $p_{u}(t,x)$ can be interpreted as a
characteristic of the sparse texture of the function $u$ that can change with
time. The following result plays a crucial role in the sequel (for the proof
we refer to [30, Lemma 2.1]).

######  Lemma 2.

Let $u\in C([0,T];L^{2}(\Omega))$ be a measurable function extended by zero
outside of $Q_{T}$. Let

$p_{u}=1+g\left(\frac{1}{h}\int_{t-h}^{t}\left|\left(\nabla
G_{\sigma}\ast\widetilde{u}(\tau,\cdot)\right)\right|^{2}\,d\tau\right)$

be the corresponding texture index. Then there exists a constant $C>0$
depending on $\Omega$, $G$, and $\|u\|_{C([0,T];L^{2}(\Omega))}$ such that

$\displaystyle\alpha:=1+\delta\leq
p_{u}(t,x)\leq\beta:=2,\quad\forall\,(t,x)\in Q_{T},$ (6) $\displaystyle
p_{u}\in C^{0,1}(Q_{T}),\ |p_{u}(t,x)-p_{u}(s,y)|\leq
C\left(|x-y|+|t-s|\right),\ \forall\,(t,x),(s,y)\in\overline{Q_{T}},$ (7)

where

$\displaystyle\delta$
$\displaystyle={ah}\left[{ah+\|G_{\sigma}\|^{2}_{C^{1}(\overline{\Omega-\Omega})}|\Omega|\|u\|^{2}_{L^{2}(0,T;L^{2}(\Omega))}}\right]^{-1},$
(8) $\displaystyle\|G_{\sigma}\|_{C^{1}(\overline{\Omega-\Omega})}$
$\displaystyle:=\max\limits_{z=x-y\atop
x\in\overline{\Omega},y\in\overline{\Omega}}\Big{[}|G_{\sigma}(z)|+|\nabla
G_{\sigma}(z)|\Big{]}=\frac{e^{-1}}{\left(\sqrt{2\pi}\sigma\right)^{2}}\left[1+\frac{1}{\sigma^{2}}\mathrm{diam}\,\Omega\right].$
(9)

##  3 Data Fusion Problem. Main Requirements to the Formal Statement

Let $\left\\{S_{1},S_{2},\dots,S_{N}:G_{H}\to\mathbb{R}^{m}\right\\}$, with
$m=13$, be a collection of multispectral images of some territory from
Sentinel-2 that were taken at time instances
$\left\\{t_{1},t_{2},\dots,t_{N}\right\\}\subset[0,T]$, respectively. We admit
that these images can be corrupted because of poor weather conditions, such as
rain, clouds, fog, and dust conditions. Typically, the measure of degradation
of optical satellite images can be such that we cannot even rely on some
reliability of pixel’s values inside the damaged regions for each of spectral
channels. As a result, some subdomains of such images become absolutely
invisible. Let $\left\\{D_{1},D_{2},\dots,D_{N}\right\\}\subset 2^{\Omega}$ be
a collection of damage regions for the corresponding Sentinel-images. So, in
fact, we deal with the set of images

$\left\\{S_{i}:G_{H}\setminus D_{i}\to\mathbb{R}^{m},\ i=1,\dots,N\right\\}.$

Let $M:G_{L}\to\mathbb{R}^{n}$, with $n=6$, be a MODIS image of the same
territory and this image has been captured at time $t=t_{M}\in(t_{1},T)$.

Before proceeding further, we begin with the following assumption:

  1. (a)

$D_{1}=\emptyset$ and the damage zones for the rest images from Sentinel-2 are
such that each $D_{i}$, $i=2,\dots,N$, is a measurable closed subset of
$\Omega$ with property $\mathcal{L}^{2}(D_{i})\leq
0.6\,\mathcal{L}^{2}(\Omega)$, where $\mathcal{L}^{2}(D_{i})$ stands for the
$2$-D Lebesgue measure of $D_{i}$;

  2. (b)

The MODIS image $M:G_{L}\to\mathbb{R}^{n}$ is assumed to be cloud-free;

  3. (c)

The images $M:G_{L}\to\mathbb{R}^{6}$ and
$\left\\{S_{1},S_{2},\dots,S_{N}:G_{H}\to\mathbb{R}^{m}\right\\}$ are rigidly
co-registered. This means that the MODIS image after arguably some affine
transformation and each Sentinel images after the resampling to the grid with
low resolution $G_{L}$, could be successfully matched according to the unique
geographic location;

  4. (d)

There exists an impulse response $\mathcal{K}$ such that, for any multi-band
image with high resolution $I:G_{H}\rightarrow\mathbb{R}^{m}$, its resampling
to the grid with low resolution $G_{L}$ can be expressed as follows

$I(x_{i},y_{j})=\left[\mathcal{K}\ast
I\right](x_{i},y_{j}),\quad\forall\,i=1,\dots,M_{x},\
\forall\,j=1,\dots,M_{y},$

where $\mathcal{K}\ast I$ stands for the convolution operator. For instance,
setting $\mathcal{K}=[k_{p,q}]_{p,q=1,\dots,K}$, we have

$\left[\mathcal{K}\ast
I\right](x_{i},y_{j})=\sum_{p=1}^{K}\sum_{q=1}^{K}k_{p,q}I(x_{i-p+1},y_{j-q+1})$

provided $I(x,y)=0$ if $(x,y)\not\in\Omega$. In majority of cases, it is
enough to set up $k_{p,q}=\frac{1}{K^{2}}$, $\forall\,p,q=1,\dots,K$, with an
appropriate choice of $K\in\mathbb{N}$. However, if we deal with satellite
images containing some agricultural areas with medium sides fields of various
shapes, then the more efficient way for the choice of kernel
$\mathcal{K}=[k_{p,q}]_{p,q=1,\dots,K}$ is to define it using the weight
coefficients of the Lanczos interpolation filters.

For our further analysis, we make use of the following notion. We say that the
multi-band images
$\left\\{\widehat{S}_{i}:G_{H}\rightarrow\mathbb{R}^{m}\right\\}_{i=1}^{N}$
are structural prototypes of the corresponding cloud-corrupted ones
$\left\\{S_{i}:G_{H}\setminus D_{i}\to\mathbb{R}^{m}\right\\}_{i=1}^{N}$ if
they are defined as follows:

$\begin{split}\widehat{S}_{1,k}(z)&=S_{1,k}(z),\\\
\widehat{S}_{i,k}(z)&=\left\\{\begin{array}[]{ll}S_{i,k}(z),&z\in
G_{H}\setminus D_{i},\\\ \gamma_{i,k}\widehat{S}_{i-1,k}(z),&z\in G_{H}\cap
D_{i},\end{array}\right\\},\ i=2,\dots,N,\ k=1,\dots,m,\end{split}$ (10)

where

$\gamma_{i,k}=\frac{\Big{[}\chi_{\Omega\setminus
D_{i}}S_{i,k}\Big{]}\cdot\Big{[}\chi_{\Omega\setminus
D_{i}}\widehat{S}_{i-1,k}\Big{]}}{\big{\|}\chi_{\Omega\setminus
D_{i}}\widehat{S}_{i-1,k}\big{\|}^{2}_{\mathcal{L}(\mathbb{R}^{2},\mathbb{R}^{2})}},$

$A\cdot B$ stands for the scalar product of two matrices $A$ and $B$, and
$\|\cdot\|_{\mathcal{L}(\mathbb{R}^{2},\mathbb{R}^{2})}$ denotes the Euclidean
norm of a matrix. Moreover, each structural prototype
$\widehat{S}_{i}:G_{H}\rightarrow\mathbb{R}^{m}$ is rigidly related to the
corresponding day $t_{i}$ when the image $S_{i}:G_{H}\setminus
D_{i}\to\mathbb{R}^{m}$ had been captured.

######  Remark 3.1.

As follows from the rule (10), this iterative procedure should be applied to
each spectral channel of all multi-band images from Sentinel-2. Since the
revisit time for Sentinel-2 is $3$–$5$ days and the collection of images
$\left\\{S_{i}:G_{H}\setminus D_{i}\to\mathbb{R}^{m}\right\\}_{i=1}^{N}$ is
rigidly co-registered, it follows from (10) that the structural prototypes
$\left\\{\widehat{S}_{i}:G_{H}\rightarrow\mathbb{R}^{m}\right\\}_{i=1}^{N}$
are also well co-registered and they have the similar topographic maps with
respect to their precise space location, albeit some false contours can appear
along the boundaries of the damage zones $D_{i}$. In fact, in order to avoid
the appearance of the false contours, the weight coefficients $\gamma_{i,k}$
have been introduced.

Since the MODIS image has been captured at a time instance
$t_{M}\in(t_{1},T)$, we can have three possible cases:

(A1)

there exists an index $i^{\ast}\in\left\\{1,2,\dots,N\right\\}$ such that
$t_{M}=t_{i^{\ast}}$;

(A2)

there exists an index $i^{\ast}\in\left\\{1,2,\dots,N-1\right\\}$ such that
$t_{i^{\ast}}<t_{M}<t_{i^{\ast}+1}$;

(A3)

$t_{N}<t_{M}<T$.

In view of this, we will distinguish three different statements of the data
fusion problem:

Case (A1)

(Restoration Problem) The problem (A1) consists in restoration of the damaged
multi-band optical image $S_{i^{\ast}}:G_{H}\setminus
D_{i^{\ast}}\to\mathbb{R}^{m}$ using result of its fusion with the cloud-free
MODIS image $M:G_{L}\to\mathbb{R}^{6}$ of the same territory. It means that,
we have to create a new image $S^{rest}_{i^{\ast}}:G_{H}\to\mathbb{R}^{m}$,
which would be well defined on the entire grid $G_{H}$, such that

$\displaystyle S^{rest}_{i^{\ast}}(z)=S_{i^{\ast}}(z),\quad\forall\,z=(x,y)\in
G_{H}\setminus D_{i^{\ast}},$ (11) $\displaystyle\sum\limits_{z\in G_{L}\cap
D_{i^{\ast}}}\Big{(}\left(\mathcal{K}\ast
S^{rest}_{i,k}\right)(z)-M_{k}(z)\Big{)}^{2}=\inf\limits_{I\in\mathcal{I}}\sum\limits_{z\in
G_{L}\cap D_{i^{\ast}}}\Big{(}\left(\mathcal{K}\ast
I\right)(z)-M_{k}(z)\Big{)}^{2},\atop\forall\,k\in J_{1},$ (12) $\displaystyle
S^{rest}_{i^{\ast},k}(z)=\widehat{S}_{i^{\ast},k}(z),\quad\forall\,z\in
G_{H},\ \forall\,k\in J_{2}.$ (13)

The precise description of the class of admissible (or feasible) images
$\mathcal{I}$ will be given in the next section.

Case (A2)

(Interpolation Problem) The problem (A2) consists in generation of a new
multi-band optical image $S^{int}_{t_{M}}:G_{H}\to\mathbb{R}^{m}$ at the
Sentinel-level of resolution using result of the fusion of cloud-free MODIS
image $M:G_{L}\to\mathbb{R}^{6}$ with the predicted structural prototype
$\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ from the given day
$t_{M}$. In fact, in this case we deal with the two-level problem. At the
first level, having the collection of structural prototypes
$\left\\{\widehat{S}_{i}:G_{H}\rightarrow\mathbb{R}^{m}\right\\}_{i=1}^{N}$
which is associated with the time instances
$\left\\{t_{1},t_{2},\dots,t_{N}\right\\}\subset[0,T]$, we create a new
’intermediate’ image $\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ that
can be considered as daily prediction of the topographical map of a given
territory from the day $t_{M}$. Then, at the second level, we realize the
fusion procedure of this predicted image with the cloud-free MODIS image
$M:G_{L}\to\mathbb{R}^{6}$ of the same territory. As a result, we have to
create a new image $S^{interp}_{t_{M}}:G_{H}\to\mathbb{R}^{m}$ with
properties:

$\displaystyle\sum\limits_{z\in G_{L}}\Big{(}\left(\mathcal{K}\ast
S^{int}_{t_{M},k}\right)(z)-M_{k}(z)\Big{)}^{2}=\inf\limits_{I\in\mathcal{I}}\sum\limits_{z\in
G_{L}}\Big{(}\left(\mathcal{K}\ast
I\right)(z)-M_{k}(z)\Big{)}^{2},\atop\forall\,k\in J_{1},$ (14) $\displaystyle
S^{int}_{t_{M},k}(z)=\widehat{S}_{t_{M},k}(z),\quad\forall\,z\in G_{H},\
\forall\,k\in J_{2}.$ (15)

Case (A3)

(Extrapolation Problem) The problem (A3) consists in generation of a new
multi-band optical image $S^{ext}_{t_{M}}:G_{H}\to\mathbb{R}^{m}$ using result
of the data assimilation from the cloud-free MODIS image
$M:G_{L}\to\mathbb{R}^{6}$ into the structural prototype
$\widehat{S}_{N}:G_{H}\rightarrow\mathbb{R}^{m}$ of the Sentinel-image
$S_{N}:G_{H}\setminus D_{N}\to\mathbb{R}^{m}$. Here, it is assumed that the
level sets of the given territory (or topographical map) for each Sentinel
spectral channel from the day $t_{M}$ have the same geo-location as they have
in $\widehat{S}_{N}:G_{H}\rightarrow\mathbb{R}^{m}$. So, we can set
$\widehat{S}_{t_{M}}=\widehat{S}_{N}$. Thus, in the framework of this problem,
we have to retrieve a new image $S^{ext}_{t_{M}}:G_{H}\to\mathbb{R}^{m}$,
which would be well defined on the entire grid $G_{H}$, such that

$\displaystyle\sum\limits_{z\in G_{L}}\Big{(}\left(\mathcal{K}\ast
S^{ext}_{t_{M},k}\right)(z)-M_{k}(z)\Big{)}^{2}=\inf\limits_{I\in\mathcal{I}}\sum\limits_{z\in
G_{L}}\Big{(}\left(\mathcal{K}\ast
I\right)(z)-M_{k}(z)\Big{)}^{2},\atop\forall\,k\in J_{1},$ (16) $\displaystyle
S^{ext}_{t_{M},k}(z)=\widehat{S}_{t_{M},k}(z),\quad\forall\,z\in G_{H},\
\forall\,k\in J_{2}.$ (17)

To provide the detailed analysis of the above mentioned problems, we begin
with some auxiliaries.

##  4 The Model for Prediction of Structural Prototypes

Due to the iterative procedure (10), we can define the so-called structural
prototypes
$\left\\{\widehat{S}_{i}:G_{H}\rightarrow\mathbb{R}^{m}\right\\}_{i=1}^{N}$
for each cloud-corrupted Sentinel image $\left\\{S_{i}:G_{H}\setminus
D_{i}\to\mathbb{R}^{m}\right\\}_{i=1}^{N}$. Let $i^{\ast}$ be an integer such
that $t_{i^{\ast}}<t_{M}<t_{i^{\ast}+1}$. Let $\widehat{S}_{i^{\ast},j}$ and
$\widehat{S}_{i^{\ast}+1,j}$ be structural prototypes of the corresponding
images from given days $t_{i^{\ast}}$ and $t_{i^{\ast}+1}$. Since
$\widehat{S}_{i^{\ast},j}$ and $\widehat{S}_{i^{\ast}+1,j}$ are well co-
registered images, it is reasonable to assume that they have the similar
geometric structure albeit they may have rather different intensities.

The main question we are going to discuss in this section is: how to correctly
define the ’intermediate’ image
$\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ that can be considered as
daily prediction of the topographical map of a given territory from the day
$t_{M}$. With that in mind, for each spectral channel
$j\in\left\\{1,2,\dots,m\right\\}$, we make use of the following model

$\displaystyle\displaystyle\frac{\partial u}{\partial
t}-\mathop{\mathrm{div}}\left(|\nabla u|^{p_{u}(t,x)-2}\nabla
u\right)=v\quad\text{in }\ (t_{i^{\ast}},t_{i^{\ast}+1})\times\Omega,$ (18)
$\displaystyle\partial_{\nu}u=0\quad\text{on }\
(t_{i^{\ast}},t_{i^{\ast}+1})\times\partial\Omega,$ (19) $\displaystyle
u(t_{i^{\ast}},\cdot)=\widehat{S}_{i^{\ast},j}(\cdot)\quad\text{in }\ \Omega,$
(20)

where $p_{u}(t,x)$ stands for the texture index of the scalar image $u$ (see
Definition 1), and $v\in L^{2}(t_{i^{\ast}},t_{i^{\ast}+1};L^{2}(\Omega))$ is
an unknown source term that has to be defined in the way to guarantee the
fulfillment (with some accuracy) of the relation

$u(t_{i^{\ast}+1},\cdot)\approx\widehat{S}_{i^{\ast}+1,j}(\cdot)\quad\text{ in
}\ \Omega.$ (21)

Here, we assume that the images $\widehat{S}_{i^{\ast},j}$ and
$\widehat{S}_{i^{\ast}+1,j}$ in (20) and (21) are well defined into the entire
domain $\Omega$.

######  Remark 4.1.

The main characteristic feature of the proposed initial-boundary value problem
(IBVP) is the fact that the exponent $p_{u}$ depend not only on $(t,x)$ but
also on a solution $u(t,x)$ of this problem. It is well-known that the
variable character of the exponent $p_{u}$ causes a gap between the
monotonicity and coercivity conditions. Because of this gap, equations of the
type (18) can be termed as equations with nonstandard growth conditions. So,
in fact, we deal with the Cauchy-Neumann IBVP for a parabolic equation of
$p_{u}=p(t,x,u)$-Laplacian type with variable exponent of nonlinearity. It was
recently shown that the model (18)–(20) naturally appears as the Euler-
Lagrange equation in the problem of restoration of cloud contaminated
satellite optical images [14, 25]. Moreover, the above mentioned problem can
be considered as a model for the deblurring and denoising of multi-spectral
images. In particular, this model has been proposed in [15, 31] in order to
avoid the blurring of edges and other localization problems presented by
linear diffusion models in images processing. We also refer to [29], where the
authors study some optimal control problems associated with a special case of
the model (18)–(20) and show that the given class of optimal control problems
is well posed.

Before proceeding further, we note that the distributed control $v$ in the
right hand side of (18) describes the fictitious sources or sinks of the
intensity $u$ that may have a tendency to change at most pixels even for co-
registered structural prototypes $\widehat{S}_{i^{\ast},j}(\cdot)$ and
${S}_{i^{\ast}+1,j}(\cdot)$. As for the Neumann boundary condition
$\partial_{\nu}u=0$ on $\partial\Omega$, this condition corresponds to the
reflection of the image across the boundary and has the advantage of not
imposing any value on the boundary and not creating ’edges’ on it. So, it is
very natural conditions if we assume that the boundary of the image is an
arbitrary cutoff of a larger scene in view.

In order to characterize the solvability issues of the IBVP (18)–(20), we
adopt the following concept.

######  Definition 3.

We say that, for given $v\in L^{2}(\Omega)$ and $\widehat{S}_{i^{\ast},j}\in
L^{2}(\Omega)$, a function $u$ is a weak solution to the problem (18)–(20) if

$u\in L^{2}(t_{i^{\ast}},t_{i^{\ast}+1};L^{2}(\Omega)),\ u(t,\cdot)\in
W^{1,1}(\Omega)\ \text{ for a.a.
$t\in[t_{i^{\ast}},t_{i^{\ast}+1}]$},\atop\displaystyle\int_{t_{i^{\ast}}}^{t_{i^{\ast}+1}}\int_{\Omega}|\nabla
u|^{p_{u}(t,x)}\,dxdt<+\infty,$ (22)

and the integral identity

$\int_{t_{i^{\ast}}}^{t_{i^{\ast}+1}}\int_{\Omega}\left(-u\frac{\partial\varphi}{\partial
t}+\left(|\nabla u|^{p_{u}}\nabla u,\nabla\varphi\right)\right)\,dxdt\\\
=\int_{t_{i^{\ast}}}^{t_{i^{\ast}+1}}\int_{\Omega}v\varphi\,dxdt+\int_{\Omega}\widehat{S}_{i^{\ast},j}\varphi|_{t=t_{i^{\ast}}}\,dx$
(23)

holds true for any function $\varphi\in\Phi$, where $\Phi=\left\\{\varphi\in
C^{\infty}([t_{i^{\ast}},t_{i^{\ast}+1}]\times\overline{\Omega})\
:\left.\varphi\right|_{t=t_{i^{\ast}+1}}=0\right\\}$.

The following result highlights the way in what sense the weak solution takes
the initial value $u(t_{i^{\ast}},\cdot)=\widehat{S}_{i^{\ast},j}(\cdot)$.

######  Proposition 4 ([30]).

Let $v\in H^{1}(\Omega)$ and $\widehat{S}_{i^{\ast},j}\in L^{2}(\Omega)$ be
given distributions. Let $u$ be a weak solution to the problem (18)–(20) in
the sense of Definition 3. Then, for any $\eta\in
C^{\infty}(\overline{\Omega})$, the scalar function
$h(t)=\displaystyle\int_{\Omega}u(t,x)\eta(x)\,dx$ belongs to
$W^{1,1}(t_{i^{\ast}},t_{i^{\ast}+1})$ and
$h(0)=\displaystyle\int_{\Omega}\widehat{S}_{i^{\ast},j}(x)\eta(x)\,dx$.

Utilizing the perturbation technique and the classical fixed point theorem of
Schauder [35], it has been recently proven the following existence result.

######  Theorem 5 ([30]).

Let $v\in L^{2}(\Omega)$ and $\widehat{S}_{i^{\ast},j}\in L^{2}(\Omega)$ be
given distributions. Then the initial-boundary value problem (18)–(20) admits
at least one weak solution $u=u(t,x)$ with the following higher inegrability
properties

$u\in L^{\infty}(t_{i^{\ast}},t_{i^{\ast}+1};L^{2}(\Omega)),\ u\in
W^{1,\alpha}(\left(t_{i^{\ast}},t_{i^{\ast}+1}\right)\times\Omega),\atop\displaystyle
u\in L^{2\alpha}\left(t_{i^{\ast}},t_{i^{\ast}+1};L^{2\alpha}(\Omega)\right),$
(24)

where the exponent $\alpha$ is given by the rule

$\displaystyle\alpha$
$\displaystyle={ah}\left[{ah+\|G_{\sigma}\|^{2}_{C^{1}(\overline{\Omega-\Omega})}|\Omega|\left(\|v\|^{2}_{L^{2}(Q_{T})}+2\|\widehat{S}_{i^{\ast},j}\|^{2}_{L^{2}(\Omega)}\right)}\right]^{-1}$

In order to satisfy the condition (21) and define an appropriate source term
$v=v(t,x)$, we utilize some issues coming from the well-known method of Horn
and Schunck [21] that has been developed in order to compute optical flow
velocity from spatiotemporal derivatives of image intensity. Following this
approach, we define the function $v^{\ast}$ as a solution of the problem

$\int_{\Omega}\left(\left.\frac{\partial Y}{\partial
t}\right|_{t=(t_{i^{\ast}}+t_{i^{\ast}+1})/2}-\left.\mathop{\mathrm{div}}\left(|\nabla
Y|^{p_{Y}}\nabla
Y\right)\right|_{t=(t_{i^{\ast}}+t_{i^{\ast}+1})/2}-v\right)^{2}\,dx\\\
+\lambda_{1}^{2}\int_{\Omega}|\nabla v|^{2}\,dx\rightarrow\inf_{v\in
H^{1}(\Omega)},$ (25)

where $\lambda_{1}>0$ is a tuning parameter (for numerical simulations we take
$\lambda_{1}=0.5$), and the spatiotemporal derivatives are computed by the
rules

$\displaystyle\left.\frac{\partial Y}{\partial
t}\right|_{t=(t_{i^{\ast}}+t_{i^{\ast}+1})/2}$
$\displaystyle=\frac{1}{t_{i^{\ast}+1}-t_{i^{\ast}}}\left(\widehat{S}_{i^{\ast}+1,j}-\widehat{S}_{i^{\ast},j}\right),$
$\displaystyle\left.\mathop{\mathrm{div}}\left(|\nabla Y|^{p_{Y}}\nabla
Y\right)\right|_{t=(t_{i^{\ast}}+t_{i^{\ast}+1})/2}$
$\displaystyle=\frac{1}{2}\Big{[}\mathop{\mathrm{div}}\left(|\nabla\widehat{S}_{i^{\ast},j}|^{p_{\widehat{S}_{i^{\ast},j}}}\nabla\widehat{S}_{i^{\ast},j}\right)$
$\displaystyle\quad+\mathop{\mathrm{div}}\left(|\nabla\widehat{S}_{i^{\ast}+1,j}|^{p_{\widehat{S}_{i^{\ast}+1,j}}}\nabla\widehat{S}_{i^{\ast}+1,j}\right)\Big{]}.$

It is clear that a minimum point $v^{\ast}\in H^{1}(\Omega)$ to unconstrained
minimization problem (25) is unique and satisfies necessarily the Euler-
Lagrange equation

$\displaystyle\lambda_{1}^{2}\Delta v^{\ast}+\left(\left.\frac{\partial
Y}{\partial
t}\right|_{t=(t_{i^{\ast}}+t_{i^{\ast}+1})/2}-\left.\mathop{\mathrm{div}}\left(|\nabla
Y|^{p_{Y}}\nabla
Y\right)\right|_{t=(t_{i^{\ast}}+t_{i^{\ast}+1})/2}-v^{\ast}\right)=0$ (26)

with the Nuemann boundary condition $\partial_{\nu}v^{\ast}=0$ on
$\partial\Omega$.

Setting $v=v^{\ast}$ in (18), we can define a function
$u^{\ast}=u^{\ast}(t,x)$ as the weak solution of the IBVP (18)–(20). Numerical
experiments show that, following this way, we obtain a function $u^{\ast}$
with properties (22) and (24) such that

$u^{\ast}(t_{i^{\ast}},x)=\widehat{S}_{i^{\ast},j}(x)\quad\text{and}\quad
u^{\ast}(t_{i^{\ast}+1},x)\approx\widehat{S}_{i^{\ast}+1,j}(x)\quad\text{in }\
\Omega,$

where the peak signal-to-noise ratio (PSNR) between images
$u^{\ast}(t_{i^{\ast}+1},x)$ and $\widehat{S}_{i^{\ast}+1,j}(x)$ is
sufficiently large, $PSNR>46$.

This observation leads us to the following conclusion: the ’intermediate’
image $\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ can be defined as
follow:

$\widehat{S}_{t_{M}}(x)=u^{\ast}(t_{M},x),\quad\forall\,x\in G_{H}.$ (27)

To illustrate how the proposed model (18)–(20) works, we consider as an input
data two images $S_{1}$, $S_{2}$ of some region that represent a typical
agricultural area in Australia the resolution $20m/pixel$. These images have
been delivered from Sentinel-2 and captured at the time instances
$t_{1}=^{\prime}July,08^{\prime}$ and $t_{2}=^{\prime}August,25^{\prime}$,
respectively. Both of these images are cloud-free (so, we can set
$\widehat{S}_{1}=S_{1}$ and $\widehat{S}_{2}=S_{2}$) and their spectral
energies ($Y_{S_{1}}$ and $Y_{S_{2}}$) are depicted in Figure 1.

Fig. 1:  The spectral energies of the Sentinel images $S_{1}$ and $S_{2}$ with
resolution $20m/pixel$. The real size of each image is $1000\times 1000$
pixels.

Setting $p_{u}$ in the form (5) with $h=0.1$, $a=0.01$, $\sigma=1$, and
defining the function $v$ as a solution of the Neumann boundary value problem
(26), we pass to the numerical solution of the IBVP (18)–(20). To this end, we
use an implicit discretization in time of equation (18), and after we apply
the conjugate gradient method. As a result, we define a function
$u^{\ast}(t,x)$ such that $u^{\ast}(t_{1},x)=S_{1}(x)$, $\forall\,x\in G_{H}$,
and the peak signal-to-noise ratio between $u^{(}t_{2},\cdot)$ and
$S_{2}(\cdot)$ is equal to $36.41$. It means that we can guarantee the
equality $u^{\ast}(t_{2},\cdot)\approx S_{2}(\cdot)$ at the high level of
accuracy. Let $t_{3}=^{\prime}July,18^{\prime}$ and
$t_{4}=^{\prime}August,10^{\prime}$ be some intermediate time point on the
interval $[t_{1},t_{2}]$. In Figure 2 the screenshots $u^{\ast}(t_{3},\cdot)$
and $u^{\ast}(t_{4},\cdot)$ of the solution to the problem (18)–(20) are
shown. Thus, the function $u^{\ast}$ can be considered as an acceptable
approximation to the evolution of the spectral energy of Sentinel images over
the time interval $[t_{1},t_{2}]$.

Fig. 2:  The screenshots $u^{(}t_{3},\cdot)$ and $u^{\ast}(t_{4},\cdot)$ of
the solution to the problem (18)–(20) taken at the time instances $t_{3}$ and
$t_{4}$, respectively.

##  5 Variational Statements of the Data Fusion Problems (A1)–(A3)

Coming back to the principle cases (A1)–(A3), that have been described in
Section 3, we can suppose that a structural prototype
$\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ from the given day
$t_{M}$ is well defined. As it was emphasized in Section 3, this prototype
coincides either with one of the images
$\left\\{\widehat{S}_{i}:G_{H}\rightarrow\mathbb{R}^{m}\right\\}_{i=1}^{N}$ in
cases (A1) and (A3), or it is defined using the solutions of the problem
(18)–(20), (26) for each $j=1,\dots,m$ in the case (A2)-problem (see the rule
(27)). For further convenience, we assume that
$\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ is zero-extended outside
of $\Omega$.

Let $j\in\left\\{1,2,\dots,m\right\\}$ be a fixed index value (the number of
spectral channel). Let $q_{j}:\Omega\rightarrow\mathbb{R}$ be the texture
index of the $j$-th band for the structural prototype
$\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$, i.e.

$q_{j}(x):=1+g\left(\left|\left(\nabla
G_{\sigma}\ast\widehat{S}_{t_{M},j}\right)(x)\right|^{2}\right),\
\forall\,x\in\Omega,$ (28)

where $g{:}[0,\infty)\rightarrow(0,\infty)$ is the edge-stopping function
which we take in the form of the Cauchy law $g(s)=\frac{a}{a+s}$ with $a>0$
small enough. Let $\eta\in(0,1)$ be a given threshold. Let
$\theta_{j}=\left[\theta_{j,1}\atop\theta_{j,2}\right]\in
L^{\infty}(\Omega,\mathbb{R}^{2})$ be a vector field such that

$|\theta_{j}(x)|_{\mathbb{R}^{2}}\leq 1\ \text{ and }\
\left(\theta_{j}(x),\nabla\widehat{S}_{t_{M},j}(x)\right)_{\mathbb{R}^{2}}=|\nabla\widehat{S}_{t_{M},j}(x)|_{\mathbb{R}^{2}}\quad\text{a.e.
in $\Omega$.}$

As it was mentioned in Subsection 2.2, for each spectral channel
$j\in\left\\{1,\dots,m\right\\}$ this vector-field can be defined by the rule
$\theta_{j}(x)=\frac{\nabla U_{j}(t,x)}{|\nabla U_{j}(t,x)|}$ with $t>0$ small
enough, where $U_{j}(t,x)$ is a solution the following initial-boundary value
problem

$\displaystyle\frac{\partial U}{\partial t}=\mathrm{div}\,\left(\frac{\nabla
U}{|\nabla U|+{\varepsilon}}\right),\quad t\in(0,+\infty),\ x\in\Omega,$ (29)
$\displaystyle U(0,x)=\widehat{S}_{j,t_{M}}(x),\quad x\in\Omega,$ (30)
$\displaystyle\frac{\partial U(0,x)}{\partial\nu}=0,\quad t\in(0,+\infty),\
x\in\partial\Omega$ (31)

with a relaxed version of the $1D$-Laplace operator in the principle part of
(29). Here, ${\varepsilon}>0$ is a sufficiently small positive value.

Taking into account the definition of the Directional Total Variation (see
[5]), we define a linear operator
$R_{j,\eta}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ as follows

$R_{j,\eta}\nabla v:=\nabla v-\eta^{2}\left(\theta_{j},\nabla
v\right)_{\mathbb{R}^{2}}\theta_{j},\quad\forall\,v\in W^{1,1}(\Omega).$ (32)

It is clear that $R_{j,\eta}\nabla v$ reduces to $(1-\eta^{2})\nabla v$ in
those regions where the gradient $\nabla v$ is co-linear to $\theta_{j}$, and
to $\nabla v$, where $\nabla v$ is orthogonal to $\theta_{j}$. It is important
to emphasize that this operator does not enforce gradients in the direction
$\theta_{j}$.

Let $\delta_{(x_{i},y_{j})}$ be the Dirac’s delta at the point
$(x_{i},y_{j})$. Then $\Pi_{S_{L}}=\sum_{(x_{i},y_{j})\in
S_{L}}\delta_{(x_{i},y_{j})}$ stands for the Dirac’s comb defined onto the
sample grid $G_{L}$.

We are now in a position to give a precise meaning of the solutions to
Problems (A1)–(A3). We say that:

  1. (A1)

A multi-band image $S^{rest}_{i^{\ast}}:G_{H}\to\mathbb{R}^{m}$, where
$t_{i^{\ast}}=t_{M}$, is a solution of the Restoration Problem, if it is given
by the rule

$\displaystyle S^{rest}_{i^{\ast},j}(x)$
$\displaystyle=\left\\{\begin{array}[]{ll}S_{i^{\ast},j}(x)&,\ \forall\,x\in
G_{H}\setminus D_{i^{\ast}},\\\ \beta_{j}u^{0}_{j}(x)&,\ \forall\,x\in
G_{H}\cap D_{i^{\ast}},\end{array}\right.\ \forall\,j\in J_{1},$ (35)
$\displaystyle S^{rest}_{i^{\ast},j}(x)$
$\displaystyle=\widehat{S}_{t_{M},j}(x),\quad\forall\,x\in G_{H},\
\forall\,j\in J_{2}.$ (36)

Here, $\beta_{j}$ is the weight coefficient and we define it as follows

$\beta_{j}=\left[\int_{\Omega\setminus
D_{i^{\ast}}}S_{i^{\ast},j}(x)u^{0}_{j}(x)\,dx\right]/\left[\int_{\Omega\setminus
D_{i^{\ast}}}|u^{0}_{j}(x)|^{2}\,dx\right]$

and $u^{0}_{j}$ is a solutions of the following constrained minimization
problem

$\left(\mathcal{P}\right)\qquad\qquad\mathcal{F}_{j}\left(u^{0}_{j}\right)=\inf_{u\in\Xi_{j}}\mathcal{F}_{j}(u),$
(37)

where

$\mathcal{F}_{j}(u):=\int_{\Omega}\frac{1}{q_{j}(x)}|R_{j,\eta}\nabla
u(x)|^{q_{j}(x)}\,dx+\frac{\mu}{2}\int_{\Omega}\left|\nabla
u(x)-\nabla\widehat{S}_{t_{M},j}(x)\right|^{2}\,dx\\\
+\frac{\gamma}{2}\int_{\Omega\setminus
D_{i^{\ast}}}\left|u(x)-S_{t_{M},j}(x)\right|^{2}\,dx+\frac{\vartheta}{2}\int_{\Omega}\Pi_{S_{L}}\Big{(}\left[\mathcal{K}\ast
u-M_{t_{M},j}\right]\Big{)}^{2}\,dx,$ (38)

$\Xi_{j}=\left\\{u\in W^{1,q_{j}(\cdot)}(\Omega)\ :\ 0\leq u(x)\leq C_{j}\
\text{a.e. in }\ \Omega\right\\}$ stands for the set of feasible solutions,
$\mu>0$, $\gamma>0$, $\vartheta>0$ are some weight coefficients, and
$W^{1,q_{j}(\cdot)}(\Omega)$ denotes the Sobolev space with variable exponent
(for the details we refer to Appendix C). As for the constants $C_{j}$, their
choice depends on the format of signed integer numbers in which the
corresponding intensities $S_{i,j}(x)$ are represented. In particular, it can
be $C_{j}=2^{8}-1$, $C_{j}=2^{16}-1$, and so on.

  2. (A2)

A multi-band image $S^{int}_{t_{M}}:G_{H}\to\mathbb{R}^{m}$, with
$t_{i^{\ast}}<t_{M}<t_{i^{\ast}+1}$, is a solution of the Interpolation
Problem, if it is given by the rule

$S^{int}_{t_{M},j}(x)=\left\\{\begin{array}[]{ll}\beta_{j}u^{0}_{j}(x)&,\
\forall\,j\in J_{1},\\\ \widehat{S}_{t_{M},j}(x)&,\ \forall\,j\in
J_{2},\end{array}\right.\ \forall\,x\in G_{H},$ (39)

where

$\beta_{j}=\left[\int_{\Omega}\widehat{S}_{t_{M},j}(x)u^{0}_{j}(x)\,dx\right]/\left[\int_{\Omega}|u^{0}_{j}(x)|^{2}\,dx\right]$
(40)

and $u^{0}_{j}$ is a solutions of the constrained minimization problem
(37)–(38) with $D_{i^{\ast}}=\Omega$.

  3. (A3)

A multi-band image $S^{ext}_{t_{M}}:G_{H}\to\mathbb{R}^{m}$, with
$t_{N}<t_{M}<T$, is a solution of the Extrapolation Problem, if it is given by
the rule

$S^{ext}_{t_{M},j}(x)=\left\\{\begin{array}[]{ll}\beta_{j}u^{0}_{j}(x)&,\
\forall\,j\in J_{1},\\\ \widehat{S}_{t_{M},j}(x)&,\ \forall\,j\in
J_{2},\end{array}\right.\ \forall\,x\in G_{H},$ (41)

where $u^{0}_{j}$ is a solutions of the constrained minimization problem
(37)–(38) with $D_{i^{\ast}}=\Omega$, and $\beta_{j}$ is defined as in (40).

Let us briefly discuss the relevance of the proposed minimization problem
$\left(\mathcal{P}\right)$. We begin with the motivation to the choice of the
energy functional in the form (37)–(38).

The first term in (38) can be considered as the regularization in the Sobolev-
Orlicz space $W^{1,q_{j}(\cdot)}(\Omega)$ because, for each spectral channel,
we have

$(1-\eta^{2})|\nabla u|\leq|R_{j,\eta}\nabla u|\leq|\nabla u|\quad\text{in
$\Omega$}$ (42)

with a given threshold $\eta\in(0,1)$. Hence,

$\displaystyle\int_{\Omega}\frac{1}{q_{j}(x)}|R_{j,\eta}\nabla
u(x)|^{q_{j}(x)}\,dx\geq(1-\eta^{2})^{2}\int_{\Omega}|\nabla
u(x)|^{q_{j}(x)}\,dx,\quad\forall\,u\in W^{1,q_{j}(\cdot)}(\Omega)$ (43)

and, therefore, if

$u\in\Xi_{j}\subset W^{1,q_{j}(\cdot)}(\Omega)\cap
L^{\infty}(\Omega)\quad\text{and}\quad\mathcal{F}_{j}(u)<+\infty,$ (44)

then

$\displaystyle\|u\|^{\alpha}_{W^{1,q_{j}(\cdot)}}$
$\displaystyle=\left(\|u\|_{L^{q_{j}(\cdot)}(\Omega)}+\|\nabla
u\|_{L^{q_{j}(\cdot)}(\Omega;\mathbb{R}^{2})}\right)^{\alpha}$
$\displaystyle\leq C\left(\|u\|^{\alpha}_{L^{q_{j}(\cdot)}(\Omega)}+\|\nabla
u\|^{\alpha}_{L^{q_{j}(\cdot)}(\Omega;\mathbb{R}^{2})}\right)$
$\displaystyle\stackrel{{\scriptstyle\text{by
\eqref{A1.2}}}}{{\leq}}C\left(\int_{\Omega}|u(x)|^{q_{j}(x)}\,dx+\int_{\Omega}|\nabla
u(x)|^{q_{j}(x)}\,dx+2\right)$ $\displaystyle\stackrel{{\scriptstyle\text{by
\eqref{5.8a}}}}{{\leq}}C\left(|\Omega|C_{j}^{2}+\int_{\Omega}|\nabla
u(x)|^{q_{j}(x)}\,dx+2\right)$ $\displaystyle\stackrel{{\scriptstyle\text{by
\eqref{5.8}}}}{{\leq}}C\left(|\Omega|C_{j}^{2}+2+\frac{1}{(1-\eta^{2})^{2}}\int_{\Omega}\frac{1}{q_{j}(x)}|R_{j,\eta}\nabla
u(x)|^{q_{j}(x)}\,dx\right)$ $\displaystyle\stackrel{{\scriptstyle\text{by
\eqref{5.5ab}}}}{{\leq}}C\left(|\Omega|C_{j}^{2}+2+\frac{1}{(1-\eta^{2})^{2}}\mathcal{F}_{j}(u)\right)<+\infty.$
(45)

On the other side, this term plays the role of a spatial data fidelity.
Indeed, what we are going to achieve in this interpolation problem, it is to
preserve the following property for the retrieved images at the Sentinel-2
resolution level: the geometry of each spectral channel of the retrieved image
has to be as close as possible to the geometry of the predicted structural
prototype $\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ that we obtain
either as a solution of the problem (18)–(20), (26), or as a result of the
iterative procedure (10). Formally, it means that relations

$\left(\theta_{j}^{\perp},\nabla u\right)=0\quad\text{a.e. in $\Omega$},\
\forall\,j\in J_{1},$ (46)

have to be satisfied. Hence, the magnitude
$\int_{\Omega}\Big{|}\left(\theta_{j}^{\perp},\nabla u\right)\Big{|}\,dx$ must
be small enough for each spectral channel, where $\theta_{j}$ stands for the
vector field of unit normals to the topographic map of the predicted band
$\widehat{S}_{t_{M},j}:G_{H}\rightarrow\mathbb{R}$. In order to achieve this
property, we observe that the expression $R_{j,\eta}\nabla u$ can be reduced
to $(1-\eta^{2})\nabla u$ in those places of $\Omega$ where $\nabla u$ is co-
linear to the unit normal $\theta_{j}$, and to $\nabla u$ if $\nabla u$ is
orthogonal to $\theta_{j}$.

Thus, gradients of the intensities $u$ that are aligned/co-linear to
$\theta_{j}$ are favored as long as $|\theta_{j}|>0$. Moreover, this property
is enforced by the special choice of the exponent $q_{j}(x)$. Since
$q_{j}(x)\approx 1$ in places in $\Omega$ where edges or discontinuities are
present in the predicted band $\widehat{S}_{t_{M},j}$, and $q_{j}(x)\approx 2$
in places where $\widehat{S}_{t_{M},j}$ is smooth or contains homogeneous
features, the main benefit of the energy functional (38) is the manner in
which it accommodates the local image information. For the places where the
gradient of $\widehat{S}_{t_{M},j}$ is sufficiently large (i.e. likely edges),
we deal with the so-called directional TV-based diffusion [5, 6], whereas in
the places where the gradient of $\widehat{S}_{t_{M},j}$ is close to zero
(i.e. homogeneous regions), the model becomes isotropic. Specifically, the
type of anisotropy at these ambiguous regions varies according to the strength
of the gradient. Apparently, the idea to involve the norm of
$W^{1,q_{j}(\cdot)}(\Omega)$ with a variable exponent $q_{j}(x)$ was firstly
proposed in [4] in order to reduce the staircasing effect in the TV image
restoration problem.

As for the second term in (38), it reflects the fact that the topographic map
of the retrieved image should be as close as possible to the topographic map
of predicted structural prototype
$\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$. We interpret this
closedness in its simplified form, namely, in the sense of $L^{2}$-norm of the
difference of the corresponding gradients.

It remains to say a few words about the last term in (38). Basically, this
term represents an $L^{2}$-distortion between a $j$-th spectral channel in the
MODIS image $M=\left[M_{1},M_{2},\dots,M_{6}\right]:G_{L}\to\mathbb{R}^{6}$
and the corresponding channel of the retrieved image $u^{0}_{j}$ which is
resampled to the grid of low resolution $G_{L}$.

##  6 Existence Result and Optimality Conditions for Constrained Minimization
Problem $\left(\mathcal{P}_{k}\right)$

Our main intention in this section is to show that, for each $j\in J_{1}$,
constrained minimization problem (37)–(38) is consistent and admits at least
one solution. Because of the specific form of the energy functional
$\mathcal{F}_{j}(u)$, the minimization problem (37)–(38) is rather challenging
and we refer to [3, 4, 5, 9] for some specific details.

Following in many aspects the recent studies [11, 26] (see also [12, 14, 13,
22]), we can give the following existence result.

######  Theorem 6.

Let $\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ be a given structural
prototype for unknown image ${S}_{t_{M}}$ from Sentinel-2. Then for any given
$j\in J_{1}$, $\mu>0$, $\gamma>0$, $\vartheta>0$, and $\eta\in(0,1)$, the
minimization problem (37)–(38) admits a unique solution $u_{j}^{0}\in\Xi_{j}$.

In order to derive some optimality conditions to the problem (37)–(38) and
characterize its solution $u_{j}^{0}\in W^{1,q_{j}(\cdot)}(\Omega)$, we show
that the cost functional $\mathcal{F}_{j}:\Xi_{j}\rightarrow\mathbb{R}$ is
Gâteaux differentiable. To this end, we note that, for arbitrary $v\in
W^{1,q_{j}(\cdot)}(\Omega)$, the following assertion

$\frac{|R_{j,\eta}\nabla u^{0}_{j}(x)+tR_{\eta}\nabla
v(x)|^{q_{j}(x)}-|R_{j,\eta}\nabla u^{0}_{j}(x)|^{q_{j}(x)}}{q_{j}(x)t}\\\
\rightarrow\left(|R_{j,\eta}\nabla u^{0}_{j}(x)|^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j}(x),R_{j,\eta}\nabla v(x)\right)\ \text{ as }\ t\to 0$

holds almost everywhere in $\Omega$. Indeed, by convexity,

$|\xi|^{q_{j}(x)}-|\eta|^{q_{j}(x)}\leq
2q_{j}(x)\left(|\xi|^{q_{j}(x)-1}+|\eta|^{q_{j}(x)-1}\right)|\xi-\eta|,$

it follows that

$\left|\frac{|R_{j,\eta}\nabla u^{0}_{j}(x)+tR_{j,\eta}\nabla
v(x)|^{q_{j}(x)}-|R_{j,\eta}\nabla
u^{0}_{j}(x)|^{q_{j}(x)}}{q_{j}(x)t}\right|\\\ \leq 2\left(|R_{j,\eta}\nabla
u^{0}_{j}(x)+tR_{j,\eta}\nabla v(x)|^{q_{j}(x)-1}+|R_{j,\eta}\nabla
u^{0}_{j}(x)|^{q_{j}(x)-1}\right)|R_{j,\eta}\nabla v(x)|\\\
\leq\mathrm{const}\,\left(|R_{j,\eta}\nabla
u^{0}_{j}(x)|^{q_{j}(x)-1}+|R_{j,\eta}\nabla
v(x)|^{q_{j}(x)-1}\right)|R_{j,\eta}\nabla v(x)|.$ (47)

Taking into account that

$\int_{\Omega}|R_{j,\eta}\nabla u^{0}_{j}(x)|^{q_{j}(x)-1}|R_{j,\eta}\nabla
v(x)|\,dx\\\ \stackrel{{\scriptstyle\text{by
\eqref{AKog_1.2.1}}}}{{\leq}}2\|R_{j,\eta}\nabla
u^{0}_{j}(x)|^{q_{j}(x)-1}\|_{L^{(q_{j})^{\prime}(\cdot)}(\Omega)}\|R_{j,\eta}\nabla
v(x)|\|_{L^{q_{j}(\cdot)}(\Omega)}\\\ \leq 2\|\nabla
u^{0}_{j}(x)|^{q_{j}(x)-1}\|_{L^{(q_{j})^{\prime}(\cdot)}(\Omega,\mathbb{R}^{2})}\|\nabla
v(x)\|_{L^{q_{j}(\cdot)}(\Omega,\mathbb{R}^{2})},$

and

$\displaystyle\int_{\Omega}|R_{j,\eta}\nabla
v(x)|^{q_{j}(x)}\,dx\stackrel{{\scriptstyle\text{by
\eqref{5.8.0}}}}{{\leq}}\int_{\Omega}|\nabla v(x)|^{q_{j}(x)}\,dx$
$\displaystyle\leq\|\nabla
v\|^{2}_{L^{q_{j}(\cdot)}(\Omega,\mathbb{R}^{2})}+1,$

we see that the right hand side of inequality (47) is an $L^{1}(\Omega)$
function. Therefore,

$\int_{\Omega}\frac{|R_{j,\eta}\nabla u^{0}_{j}(x)+tR_{j,\eta}\nabla
v(x)|^{q_{j}(x)}-|R_{j,\eta}\nabla u^{0}_{j}(x)|^{q_{j}(x)}}{q_{j}(x)t}\,dx\\\
\rightarrow\int_{\Omega}\left(|R_{j,\eta}\nabla
u^{0}_{j}(x)|^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x),R_{j,\eta}\nabla
v(x)\right)\,dx\ \text{ as }\ t\to 0$

by the Lebesgue dominated convergence theorem.

Since the cost functional $\mathcal{F}_{j}:\Xi_{j}\rightarrow\mathbb{R}$ can
be cast in the form

$P_{k}(u)=\mathcal{A}_{j}(u)+\mu\mathcal{B}_{j}(u)+\gamma\mathcal{C}_{j}(u)+\vartheta\mathcal{D}_{j}(u),$

where

$\displaystyle\mathcal{A}_{j}(u)$
$\displaystyle=\int_{\Omega}\frac{1}{q_{j}(x)}|R_{j,\eta}\nabla
u(x)|^{q_{j}(x)}\,dx,$ $\displaystyle\mathcal{B}_{j}(u)$
$\displaystyle=\frac{1}{2}\int_{\Omega}\left|\nabla
u(x)-\nabla\widehat{S}_{t_{M},j}(x)\right|^{2}\,dx,$
$\displaystyle\mathcal{C}_{j}(u)$
$\displaystyle=\frac{1}{2}\int_{\Omega\setminus
D_{i^{\ast}}}\left|u(x)-S_{t_{M},j}(x)\right|^{2}\,dx,$
$\displaystyle\mathcal{D}_{j}(u)$
$\displaystyle=\frac{1}{2}\int_{\Omega}\Pi_{S_{L}}\Big{(}\left[\mathcal{K}\ast
u-M_{t_{M},j}\right]\Big{)}^{2}\,dx,$

we deduce that

$\mathcal{A}_{j}^{\prime}(u^{0}_{j})[v]=\int_{\Omega}\left(|R_{j,\eta}\nabla
u^{0}_{j}(x)|^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x),R_{j,\eta}\nabla
v(x)\right)_{\mathbb{R}^{2}}\,dx,$ (48)

for each $v\in W^{1,q_{j}(\cdot)}(\Omega)$.

As for the rest terms $\mu\mathcal{B}_{j}(u)$, $\gamma\mathcal{C}_{j}(u)$, and
$\lambda\mathcal{D}_{j}(u)$ in the cost functional
$\mathcal{F}_{j}:\Xi_{j}\rightarrow\mathbb{R}$, utilizing the similar
arguments, we have the following representation for their Gâteaux derivatives.

######  Proposition 7.

For a given MODIS image $M:G_{L}\to\mathbb{R}^{6}$, the functionals
$\mathcal{B}_{j},\mathcal{C}_{j},\mathcal{D}_{j}:L^{2}(\Omega)\rightarrow\mathbb{R}$
are convex and Gâteaux differentiable in $L^{2}(\Omega)$ with

$\displaystyle\mathcal{B}_{k}^{\prime}(u^{0}_{j})[v]$
$\displaystyle=\int_{\Omega}\left(\nabla
u^{0}_{j}(x)-\nabla\widehat{S}_{t_{M},j}(x),\nabla v(x)\right)\,dx,$ (49)
$\displaystyle\mathcal{C}_{k}^{\prime}(u^{0}_{j})[v]$
$\displaystyle=\int_{\Omega\setminus
D_{i^{\ast}}}\left(u^{0}_{j}(x)-S_{t_{M},j}(x)\right)v(x)\,dx,$ (50)
$\displaystyle\mathcal{D}_{k}^{\prime}(u^{0}_{j})[v]$
$\displaystyle=\int_{\Omega}\Pi_{L}\left(\left[\mathcal{K}\ast
u^{0}_{j}\right]-M_{t_{M},j}\right)\left[\mathcal{K}\ast v\right]\,dx$
$\displaystyle=\int_{\Omega}\Pi_{L}\left[\mathcal{K}^{\ast}\ast\left(\left[\mathcal{K}\ast
u^{0}_{j}\right]-M_{t_{M},j}\right)\right]v\,dx,$ (51)

for all $v\in W^{1,q_{j}(\cdot)}(\Omega)$.

Thus, in order to derive some optimality conditions for the minimizer
$u^{0}_{j}\in W^{1,q_{j}(\cdot)}(\Omega)$ to the problem
$\inf\limits_{u\in\Xi_{j}}\mathcal{F}_{j}(u)$, we note that $\Xi_{j}$ is a
nonempty convex subset of $W^{1,q_{j}(\cdot)}(\Omega)\cap L^{\infty}(\Omega)$
and the objective functional $\mathcal{F}_{j}:\Xi_{j}\rightarrow\mathbb{R}$ is
strictly convex. Hence, the well known classical result (see [32, Theorem
1.1.3]) leads us to the following conclusion.

######  Theorem 8.

Let $\widehat{S}_{t_{M}}:G_{H}\rightarrow\mathbb{R}^{m}$ be a given structural
prototype for unknown image ${S}_{t_{M}}$ from Sentinel-2. Let
$M:G_{L}\to\mathbb{R}^{6}$ be a given MODIS image. Let $q_{j}$ stands for the
texture index of the $j$-th band for the predicted structural prototype
$\widehat{S}_{t_{M}}$ (see (28)). Then the unique minimizer
$u^{0}_{j}\in\Xi_{j}$ to the minimization problem
$\inf\limits_{u\in\Xi_{j}}\mathcal{F}_{j}(u)$ is characterized by the
following variational inequality

$\int_{\Omega}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x),R_{j,\eta}\nabla
v(x)-R_{j,\eta}\nabla u^{0}_{j}(x)\right)\,dx\\\ +\mu\int_{\Omega}\left(\nabla
u^{0}_{j}(x)-\nabla\widehat{S}_{t_{M},j}(x),\nabla v(x)\right)\,dx\\\
+\gamma\int_{\Omega\setminus
D_{i^{\ast}}}\left(u^{0}_{j}(x)-S_{t_{M},j}(x)\right)v(x)\,dx\\\
+\vartheta\int_{\Omega}\Pi_{L}\left[\mathcal{K}^{\ast}\ast\left(\left[\mathcal{K}\ast
u^{0}_{j}\right]-M_{t_{M},j}\right)\right]v\,dx\geq 0,\
\quad\forall\,v\in\Xi_{j}.$ (52)

######  Remark 6.1.

In practical implementation, it is reasonable to define an optimal solution
$u^{0}_{j}\in\Xi_{j}$ using a ’gradient descent’ strategy. Indeed, observing
that

$\int_{\Omega}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x),\nabla
v(x)\right)\,dx\\\
=-\int_{\Omega}\mathop{\mathrm{div}}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x)\right)v\,dx\\\
+\int_{\partial\Omega}\left(\Big{|}R_{j,\eta}\nabla
I_{k,i}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j}(x),\nu\right)v\,d\mathcal{H}^{1},$

and

$\int_{\Omega}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j}(x),\left(\theta,\nabla v\right)\theta\right)\,dx\\\
=-\int_{\Omega}\mathop{\mathrm{div}}\left(\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j}(x),\theta\right)\theta\right)v\,dx\\\
+\int_{\partial\Omega}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j}(x),\theta\right)\left(\theta,\nu\right)v\,d\mathcal{H}^{1},$

we see that

$\int_{\Omega}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x),R_{j,\eta}\nabla
v(x)\right)\,dx\\\ =\int_{\Omega}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x),\nabla
v-\eta^{2}\left(\theta,\nabla v\right)\theta\right)\,dx\\\
=-\int_{\Omega}\mathop{\mathrm{div}}\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x)\right)v\,dx\\\
+\eta^{2}\int_{\Omega}\mathop{\mathrm{div}}\left(\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j}(x),\theta\right)\theta\right)v\,dx$

provided

$\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j},\nu\right)=0\quad\text{on}\ \partial\Omega.$

Thus, following the standard procedure and starting from the initial image
$\widehat{S}_{t_{M},j}$, we can pass to the following initial value problem
for the quasi-linear parabolic equations with Nuemann boundary conditions

$\displaystyle\frac{\partial u^{0}_{j}}{\partial t}-$
$\displaystyle\mathop{\mathrm{div}}\left(|R_{j,\eta}\nabla
u^{0}_{j}(x)|^{q_{j}(x)-2}R_{j,\eta}\nabla u^{0}_{j}(x)\right)$
$\displaystyle=-\eta^{2}\mathop{\mathrm{div}}\left(\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j}(x),\theta\right)\theta\right)$
$\displaystyle\quad+\mu\mathop{\mathrm{div}}\left(\nabla
u^{0}_{j}(x)-\nabla\widehat{S}_{t_{M},j}(x)\right)-\gamma\left(u^{0}_{j}(x)-{S}_{t_{M},j}(x)\right)$
$\displaystyle\quad-\vartheta\Pi_{L}\left[\mathcal{K}^{\ast}\ast\left(\left[\mathcal{K}\ast
u^{0}_{j}\right]-M_{t_{M},j}\right)\right],$ (53)
$\displaystyle\left(\Big{|}R_{j,\eta}\nabla
u^{0}_{j}(x)\Big{|}^{q_{j}(x)-2}R_{j,\eta}\nabla
u^{0}_{j},\nu\right)=0\quad\text{on}\ \partial\Omega,$ (54) $\displaystyle 0$
$\displaystyle\leq u^{0}_{j}(x)\leq C_{j}\ \text{a.a. in}\ \Omega,$ (55)
$\displaystyle
u^{0}_{j}(0,x)=\widehat{S}_{t_{M},j}(x),\quad\forall\,x\in\Omega.$ (56)

In principle, instead of the initial condition (56) we may consider other
image that can be generated from $\widehat{S}_{t_{M},j}$ and the bicubic
interpolation of the MODIS band $M_{t_{M},j}$ onto the entire domain $\Omega$.
For instance, it can be one of well-known simple data fusion methods (for the
details, we refer to [36]).

##  7 Numerical Experiments

In order to illustrate the proposed approach for the restoration of satellite
multi-spectral images we have used a series of Sentinel-2 images ($725\times
600$ in pixels) over the South Dakota area (USA) with resolution $10m/pixel$
(see Fig. 3 that were captured at different time instances in period from July
10 to July 15, 2021, when the global biophysical processes are rapid enough.
So, it was a period of an active vegetation growing, and what is more
important, each of these images have been captured at a cloud-free day. As
follows from Fig. 3, this region represents a typical agricultural area with
medium sides fields of various shapes. Since each of these images contains
four bands — $B_{4}$ (red), $B_{3}$ (green), $B_{2}$ (blue), and $B_{8a}$
(near infrared), we denote these images as
$S_{1},S_{2},S_{3}:G_{H}\to\mathbb{R}^{4}$, respectively. We also have a
cloud-free MODIS image $M:G_{L}\to\mathbb{R}^{6}$ ($29\times 24$ in pixels)
from 2021/08/13 with resolution $250m/pixel$ (see Fig. 4).

Fig. 3: Images from Sentinel-2. Date of generation: (left up) 2021/08/10,
(right up) 2021/08/13, , (center down) 2021/08/15 Fig. 4: MODIS image with the
date of generation 2021/08/13

To emulate the interpolation problem (A2), we assume that the image $S_{2}$ is
cloud-corrupted. With that in mind, we create some artificial clouds on it
(see the left picture in Fig. 5) and consider a new image
$\mathcal{S}_{2}:G_{H}\setminus D\to\mathbb{R}^{4}$ as an image with damage
zone $D$ in all bands. So, the problem consists in generation of a new multi-
band optical image $\mathcal{S}^{int}_{2}:G_{H}\to\mathbb{R}^{4}$ at the
Sentinel-level of resolution using result of the fusion of cloud-free MODIS
image $M:G_{L}\to\mathbb{R}^{6}$ with the predicted structural prototype
$\widehat{\mathcal{S}}_{2}:G_{H}\rightarrow\mathbb{R}^{m}$ from the given day
2021/08/13.

Fig. 5:  (left) Image with damage, (right) Its structural prototype

At the first step, following the procedure described in Section 4 and
utilizing for that the two cloud-free Sentinel images $S_{1}$ and $S_{3}$, we
create for each band its daily prediction of the topographical map from the
day 2021/08/13 (see the right picture in Fig. 5). After that, we realize the
fusion procedure of this predicted image with the cloud-free MODIS image
$M:G_{L}\to\mathbb{R}^{6}$ of the same territory. In all numerical
simulations, we set $h=0.1$, $a=0.01$, $\sigma=1$, ${\varepsilon}=0.001$,
$\eta=0.95$, $\mu=2.5$, $\vartheta=1$, $\gamma=0$. As a result, a new image
from the date 2021/08/13 at the Sentinel-level of resolution has been
generated and it is depicted in Fig. 6.

Fig. 6:  (left) Original image from Sentinel from 2021/08/13, (right) Result
of its restoration

To evaluate the given interpolation result, we make use of the following
validation metrics:

RMSE | $\mathrm{RMSE}\left(S_{2},\mathcal{S}^{int}_{2}\right)=\frac{1}{|G_{H}|}\sum_{x\in G_{H}}\left(S_{2}(x)-\mathcal{S}^{int}_{2}(x)\right)^{2}$
---|---
Corr | $\mathrm{Corr}\left(S_{2},\mathcal{S}^{int}_{2}\right)=\frac{\mathrm{cov}\,\left(S_{2},\mathcal{S}^{int}_{2}\right)}{\sqrt{\mathrm{cov}\,\left(S_{2},S_{2}\right)\mathrm{cov}\,\left(\mathcal{S}^{int}_{2},\mathcal{S}^{int}_{2}\right)}}$
CorrLaplace | $\mathrm{CorrLaplace}\left(S_{2},\mathcal{S}^{int}_{2}\right)=\mathrm{Corr}\left(\Delta S_{2},\Delta\mathcal{S}^{int}_{2}\right)$
SSIM | The structural similarity index combining local image structure,
| luminance, and contrast into a single local quality score
HaarPSI | The Haar wavelet-based perceptual similarity index
| that aims to correctly assess the perceptual similarity
| between two images with respect to a human viewer

where the above abbreviations stand for: RMSE (The Root-Mean-Square
Deviation), Corr (The Correlation Coefficient), CorrLaplace (The Correlation
of Laplace Feature).

Table 3 contains the comparison results of the original image $S_{2}$ with its
restored variant $\mathcal{S}^{int}_{2}$ for each band separately, whereas in
Table 4 we present the similarity results for the normalized difference
vegetation indices (NDVI) which can be calculated from these images as
follows:

$NDVI(S_{2})=\frac{B_{8a}(S_{2})-B_{4}(S_{2})}{B_{8a}(S_{2})+B_{4}(S_{2})},\qquad
NDVI(\mathcal{S}^{int}_{2})=\frac{B_{8a}(\mathcal{S}^{int}_{2})-B_{4}(\mathcal{S}^{int}_{2})}{B_{8a}(\mathcal{S}^{int}_{2})+B_{4}(\mathcal{S}^{int}_{2})}.$

Table 3:  Similarity characteristics between $S_{2}$ and its restored variant $\mathcal{S}^{int}_{2}$ Metrics | $B_{2}$ | $B_{3}$ | $B_{4}$ | $B_{8a}$
---|---|---|---|---
RMSE | 1669.9299 | 6339.3894 | 5587.9002 | 143613.5512
Corr | 0.9720 | 0.9272 | 0.9769 | 0.9220
CorrLaplace | 0.9490 | 0.8623 | 0.9229 | 0.3107
SSIM | 0.9428 | 0.9094 | 0.9313 | 0.8665
HaarPSI | 0.8008 | 0.6994 | 0.7907 | 0.6498
Table 4:  Similarity characteristics between NDVIs for $S_{2}$ and its restored variant $\mathcal{S}^{int}_{2}$ Metrics | Value
---|---
$\mathrm{RMSE}\left(\mathrm{NDVI}(S_{2}),\mathrm{NDVI}(\mathcal{S}^{int}_{2})\right)$ | 0.0002
$\mathrm{SSIM}\left(\mathrm{NDVI}(S_{2}),\mathrm{NDVI}(\mathcal{S}^{int}_{2})\right)$ | 0.9999
$\mathrm{HaarPSI}\left(\mathrm{NDVI}(S_{2}),\mathrm{NDVI}(\mathcal{S}^{int}_{2})\right)$ | 0.9999

As for the choice of numerical schemes for the problems (26) and (53)–(56),
their consistency and substantiation, and also other scenario for simulations
including the restoration and prediction problems, these issues will be a
subject of a forthcoming paper. We will show that the proposed approach is
appealing for automated processing of large data sets.

##  8 Conclusion

We propose a novel model for the restoration of satellite multi-spectral
images. This model is based on the solutions of special variational problems
with nonstandard growth objective functional. Because of the risk of
information loss in optical images, we do not impose any information about
such images inside the damage region, but instead we assume that the texture
of these images can be predicted through a number of past cloud-free images of
the same region from the time series. So, the characteristic feature of
variational problems, which we formulate for each spectral channel separately,
is the structure of their objective functionals. On the one hand, we involve
into consideration the energy functionals with the nonstandard growth $p(x)$,
where the variable exponent $p(x)$ is unknown a priori and it directly depends
on the texture of an image that we are going to restore. On the other hand,
the texture of an image $\vec{u}$, we are going to restore, can have rather
rich structure in the damage region $D$. In order to identify it, we push
forward the following hypothesis: the geometry of each spectral channels of a
cloud corrupted image in the damage region is topologically close to the
geometry of the total spectral energy that can be predicted with some accuracy
by a number of past cloud-free images of the same region. As a result, we
impose this requirement in each objective functional in the form of a special
fidelity term. In order to study the consistency of the proposed collection of
non-convex minimization problems, we develop a special technique and supply
this approach by the rigorous mathematical substantiation.

##  Appendix A On Orlicz Spaces

Let $p(\cdot)$ be a measurable exponent function on $\Omega$ such that
$1<\alpha\leq p(x)\leq\beta<\infty$ a.e. in $\Omega$, where $\alpha$ and
$\beta$ are given constants. Let
$p^{\prime}(\cdot)=\frac{p(\cdot)}{p(\cdot)-1}$ be the corresponding conjugate
exponent. It is clear that

$1\leq\underbrace{\frac{\beta}{\beta-1}}_{\beta^{\prime}}\leq
p^{\prime}(x)\leq\underbrace{\frac{\alpha}{\alpha-1}}_{\alpha^{\prime}}\ \
a.e.\ in\ \Omega,$

where $\beta^{\prime}$ and $\alpha^{\prime}$ stand for the conjugates of
constant exponents. Denote by $L^{p(\cdot)}(\Omega)$ the set of all measurable
functions $f(x)$ on $\Omega$ such that
$\int_{\Omega}|f(x)|^{p(x)}\,dx<\infty$. Then $L^{p(\cdot)}(\Omega)$ is a
reflexive separable Banach space with respect to the Luxemburg norm (see [10,
17] for the details)

$\|f\|_{L^{p(\cdot)}(\Omega)}=\inf\left\\{\lambda>0\ :\
\rho_{p}(\lambda^{-1}f)\leq 1\right\\},\ $ (57)

where $\rho_{p}(f):=\int_{\Omega}|f(x)|^{p(x)}\,dx$.

It is well-known that $L^{p(\cdot)}(\Omega)$ is reflexive provided $\alpha>1$,
and its dual is $L^{p^{\prime}(\cdot)}(\Omega)$, that is, any continuous
functional $F=F(f)$ on $L^{p(\cdot)}(\Omega)$ has the form (see [42, Lemma
13.2])

$F(f)=\int_{\Omega}fg\,dx,\quad\ with\ \ g\in L^{p^{\prime}(\cdot)}(\Omega).$

As for the infimum in (57), we have the following result.

######  Proposition 9.

The infimum in (57) is attained if $\rho_{p}(f)>0$. Moreover

$if\ \lambda_{\ast}:=\|f\|_{L^{p(\cdot)}(\Omega)}>0,\ \ then\ \
\rho_{p}(\lambda_{\ast}^{-1}f)=1.$ (58)

Taking this result and condition $1\leq\alpha\leq p(x)\leq\beta$ into account,
we see that

$\displaystyle\frac{1}{\lambda^{\beta}_{\ast}}\int_{\Omega}\left|f(x)\right|^{p(x)}\,dx$
$\displaystyle\leq\int_{\Omega}\left|\frac{f(x)}{\lambda_{\ast}}\right|^{p(x)}\,dx\leq\frac{1}{\lambda^{\alpha}_{\ast}}\int_{\Omega}\left|f(x)\right|^{p(x)}\,dx,$
$\displaystyle\frac{1}{\lambda^{\beta}_{\ast}}\int_{\Omega}\left|f(x)\right|^{p(x)}\,dx$
$\displaystyle\leq
1\leq\frac{1}{\lambda^{\alpha}_{\ast}}\int_{\Omega}\left|f(x)\right|^{p(x)}\,dx.$

Hence, (see [10, 17, 41] for the details)

$\displaystyle\|f\|^{\alpha}_{L^{p(\cdot)}(\Omega)}$
$\displaystyle\leq\int_{\Omega}|f(x)|^{p(x)}\,dx\leq\|f\|^{\beta}_{L^{p(\cdot)}(\Omega)},\
\ if\ \|f\|_{L^{p(\cdot)}(\Omega)}>1,$
$\displaystyle\|f\|^{\beta}_{L^{p(\cdot)}(\Omega)}$
$\displaystyle\leq\int_{\Omega}|f(x)|^{p(x)}\,dx\leq\|f\|^{\alpha}_{L^{p(\cdot)}(\Omega)},\
\ if\ \|f\|_{L^{p(\cdot)}(\Omega)}<1,$ (59)

and, therefore,

$\displaystyle\|f\|^{\alpha}_{L^{p(\cdot)}(\Omega)}-1$
$\displaystyle\leq\int_{\Omega}|f(x)|^{p(x)}\,dx\leq\|f\|^{\beta}_{L^{p(\cdot)}(\Omega)}+1,\quad\forall\,f\in
L^{p(\cdot)}(\Omega),$ (60) $\displaystyle\|f\|_{L^{p(\cdot)}(\Omega)}$
$\displaystyle=\int_{\Omega}|f(x)|^{p(x)}\,dx,\ \ if\ \
\|f\|_{L^{p(\cdot)}(\Omega)}=1.$ (61)

The following estimates are well-known (see, for instance, [10, 17, 41]): if
$f\in L^{p(\cdot)}(\Omega)$ then

$\displaystyle\|f\|_{L^{\alpha}(\Omega)}$
$\displaystyle\leq\left(1+|\Omega|\right)^{1/\alpha}\|f\|_{L^{p(\cdot)}(\Omega)},$
(62) $\displaystyle\|f\|_{L^{p(\cdot)}(\Omega)}$
$\displaystyle\leq\left(1+|\Omega|\right)^{1/\beta^{\prime}}\|f\|_{L^{\beta}(\Omega)},\quad\beta^{\prime}=\frac{\beta}{\beta-1},\quad\forall\,f\in
L^{\beta}(\Omega).$ (63)

Let $\left\\{p_{k}\right\\}_{k\in\mathbb{N}}\subset
C^{0,\delta}(\overline{\Omega})$, with some $\delta\in(0,1]$, be a given
sequence of exponents. Hereinafter in this subsection we assume that

$\displaystyle p,p_{k}\in C^{0,\delta}(\overline{\Omega})\ \ for\
k=1,2,\dots,\ \mbox{and}\atop p_{k}(\cdot)\rightarrow p(\cdot)\
\mbox{uniformly in}\ \overline{\Omega}\ \mbox{as}\ k\to\infty.$ (64)

We associate with this sequence the following collection $\left\\{f_{k}\in
L^{p_{k}(\cdot)}(\Omega)\right\\}_{k\in\mathbb{N}}$. The characteristic
feature of this set of functions is that each element $f_{k}$ lives in the
corresponding Orlicz space $L^{p_{k}(\cdot)}(\Omega)$. We say that the
sequence $\left\\{f_{k}\in L^{p_{k}(\cdot)}(\Omega)\right\\}_{k\in\mathbb{N}}$
is bounded if

$\limsup_{k\to\infty}\int_{\Omega}|f_{k}(x)|^{p_{k}(x)}\,dx<+\infty.$ (65)

######  Definition 10.

A bounded sequence $\left\\{f_{k}\in
L^{p_{k}(\cdot)}(\Omega)\right\\}_{k\in\mathbb{N}}$ is weakly convergent in
the variable Orlicz space $L^{p_{k}(\cdot)}(\Omega)$ to a function $f\in
L^{p(\cdot)}(\Omega)$, where $p\in C^{0,\delta}(\overline{\Omega})$ is the
limit of $\left\\{p_{k}\right\\}_{k\in\mathbb{N}}\subset
C^{0,\delta}(\overline{\Omega})$ in the uniform topology of
$C(\overline{\Omega})$, if

$\lim_{k\to\infty}\int_{\Omega}f_{k}\varphi\,dx=\int_{\Omega}f\varphi\,dx,\quad\forall\,\varphi\in
C^{\infty}_{0}(\mathbb{R}^{N}).$ (66)

We make use of the following result (we refer to [42, Lemma 13.3] for
comparison) concerning the lower semicontinuity property of the variable
$L^{p_{k}(\cdot)}$-norm with respect to the weak convergence in
$L^{p_{k}(\cdot)}(\Omega)$.

######  Proposition 11.

If a bounded sequence $\left\\{f_{k}\in
L^{p_{k}(\cdot)}(\Omega)\right\\}_{k\in\mathbb{N}}$ converges weakly in
$L^{\alpha}(\Omega)$ to $f$ for some $\alpha>1$, then $f\in
L^{p(\cdot)}(\Omega)$, $f_{k}\rightharpoonup f$ in variable
$L^{p_{k}(\cdot)}(\Omega)$, and

$\liminf_{k\to\infty}\int_{\Omega}|f_{k}(x)|^{p_{k}(x)}\,dx\geq\int_{\Omega}|f(x)|^{p(x)}\,dx.$
(67)

######  Remark A.1.

Arguing in a similar manner and using, instead of (LABEL:A1.5.2), the estimate

$\liminf_{k\to\infty}\int_{\Omega}\frac{1}{p_{k}(x)}|f_{k}(x)|^{p_{k}(x)}\,dx\\\
\geq\int_{\Omega}f(x)\varphi(x)\,dx-\int_{\Omega}\frac{1}{p_{k}^{\prime}(x)}|\varphi(x)|^{p^{\prime}(x)}\,dx,$

it can be shown that the lower semicontinuity property (67) can be generalized
as follows

$\liminf_{k\to\infty}\int_{\Omega}\frac{1}{p_{k}(x)}|f_{k}(x)|^{p_{k}(x)}\,dx\geq\int_{\Omega}\frac{1}{p(x)}|f(x)|^{p(x)}\,dx.$
(68)

The following result can be viewed as an analogous of the Hölder inequality in
Lebesgue spaces with variable exponents (for the details we refer to [10,
17]).

######  Proposition 12.

If $f\in L^{p(\cdot)}(\Omega)^{N}$ and $g\in
L^{p^{\prime}(\cdot)}(\Omega)^{N}$, then $\left(f,g\right)\in L^{1}(\Omega)$
and

$\int_{\Omega}\left(f,g\right)\,dx\leq
2\|f\|_{L^{p(\cdot)}(\Omega)^{N}}\|g\|_{L^{p^{\prime}(\cdot)}(\Omega)^{N}}.$
(69)

##  Appendix B Sobolev Spaces with Variable Exponent

We recall here well-known facts concerning the Sobolev spaces with variable
exponent. Let $p(\cdot)$ be a measurable exponent function on $\Omega$ such
that $1<\alpha\leq p(x)\leq\beta<\infty$ a.e. in $\Omega$, where $\alpha$ and
$\beta$ are given constants. We associate with it the so-called Sobolev-Orlicz
space

$W^{1,p(\cdot)}(\Omega):=\left\\{u\in
W^{1,1}(\Omega):\int_{\Omega}\left[|u(x)|^{p(x)}+|\nabla
u(x)|^{p(x)}\right]\,dx<\infty\right\\}$ (70)

and equip it with the norm
$\|u\|_{W^{1,p(\cdot)}(\Omega)}=\|u\|_{L^{p(\cdot)}(\Omega)}+\|\nabla
u\|_{L^{p(\cdot)}(\Omega;\mathbb{R}^{N})}$.

It is well-known that, in general, unlike classical Sobolev spaces, smooth
functions are not necessarily dense in $W=W^{1,p(\cdot)}_{0}(\Omega)$. Hence,
with variable exponent $p=p(x)$ ($1<\alpha\leq p\leq\beta$) we can associate
another Sobolev space,

$H=H^{1,p(\cdot)}(\Omega)\ \mbox{as the closure of the set}\
C^{\infty}(\overline{\Omega})\ \mbox{in}\ W^{1,p(\cdot)}(\Omega)\mbox{-norm}.$

Since the identity $W=H$ is not always valid, it makes sense to say that an
exponent $p(x)$ is regular if $C^{\infty}(\overline{\Omega})$ is dense in
$W^{1,p(\cdot)}(\Omega)$.

The following result reveals the important property that guarantees the
regularity of exponent $p(x)$.

######  Proposition 13.

Assume that there exists $\delta\in(0,1]$ such that $p\in
C^{0,\delta}(\overline{\Omega})$. Then the set $C^{\infty}(\overline{\Omega})$
is dense in $W^{1,p(\cdot)}(\Omega)$, and, therefore, $W=H$.

###### Proof.

Let $p\in C^{0,\delta}(\overline{\Omega})$ be a given exponent. Since

$\lim\limits_{t\to 0}|t|^{\delta}\log(|t|)=0\quad\ with\ \delta\in(0,1],$ (71)

it follows from the Hölder continuity of $p(\cdot)$ that

$|p(x)-p(y)|\leq
C|x-y|^{\delta}\leq\left[\sup_{x,y\in\Omega}|x-y|^{\delta}\log(|x-y|^{-1})\right]\omega(|x-y|),\quad\forall\,x,y\in\Omega,$

where $\omega(t)=C/\log(|t|^{-1})$, and $C>0$ is some positive constant.

Then property (71) implies that $p(\cdot)$ is a log-Hölder continuous
function. So, to deduce the density of $C^{\infty}(\overline{\Omega})$ in
$W^{1,p(\cdot)}(\Omega)$ it is enough to refer to Theorem 13.10 in [42]. ∎

## References

  * [1] L. Ambrosio, V. Caselles, S. Masnou and J. M. Morel, _The connected components of sets of finite perimeter_ , European Journal of Math., 3 (2001), 39–92.
  * [2] C. Ballester, V. Caselles, L. Igual, J. Verdera, B. Rougé, A Variational Model for P+XS Image Fusion, International Journal of Computer Vision, 69 (2006), 43–58.
  * [3] P. Blomgren Total variation methods for restoration of vector valued images, Ph.D. Thesis, (1998), 384–387.
  * [4] P. Blomgren, T.F. Chan, P. Mulet, C. Wong, Total variation image restoration: Numerical methods and extensions, Proceedings of the 1997 IEEE International Conference on Image Processing, (1997), III:384–387.
  * [5] L. Bungert, D.A. Coomes, M.J. Ehrhardt, J. Rasch, R. Reisenhofer, R. & C.-B. Schönlieb, Blind image fusion for hyperspectral imaging with the directional total variation, Inverse Problems, 34 (4) (2018), Article 044003.
  * [6] L. Bungert, M.J. Ehrhardt Robust Image Reconstruction with Misaligned Structural Information, IEEE Access, 8 (2020), 222944–222955.
  * [7] V. Caselles, B. Coll, J.M. Morel, Topographic maps and local contrast changes in natural images, IEEE Transactions on Image Processing, 10(8) (1999), 5–27.
  * [8] V. Caselles, B. Coll, J.M. Morel, Geometry and color in natural images, J. Math. Imaging and Vision, 16 (2002), 89–107.
  * [9] Y. Chen, S. Levine, M. Rao, Variable exponent, linear growth functionals in image restoration, SIAM Journal Appl. Math. 66(4), (2006), 1383–1406.
  * [10] D.V. Cruz-Uribe, A. Fiorenza, Variable Lebesgue Spaces: Foundations and Harmonic Analysis, Birkhäuser, New York, 2013.
  * [11] C. D’Apice, U. De Maio, P.I. Kogut, An indirect approach to the existence of quasi-optimal controls in coefficients for multi-dimensional thermistor problem, in “Contemporary Approaches and Methods in Fundamental Mathematics and Mechanics”, Editors: Sadovnichiy, Victor A., Zgurovsky, Michael (Eds.). Springer. Chapter 24, (2020), 489–522.
  * [12] C. D’Apice, P.I. Kogut, R. Manzo, On Coupled Two-Level Variational Problem in Sobolev-Orlicz Space, Differential and Integral Equations, (2022), (in press).
  * [13] C. D’Apice, P.I. Kogut, O. Kupenko, R. Manzo, On Variational Problem with Nonstandard Growth Functional and Its Applications to Image Processing, Journal of Mathematical Imaging and Vision, accepted 1 Nov. 2022; https://doi.org/10.1007/s10851-022-01131-w. Publishes online.
  * [14] C. D’Apice, P.I. Kogut, R. Manzo, M.V. Uvarov, Variational Model with Nonstandard Growth Conditions for Restoration of Satellite Optical Images Using Synthetic Aperture Radar, Europian Journal of Applies Math., Published online by Cambridge University Press: 11 March 2022, https://doi.org/10.1017/S0956792522000031.
  * [15] C. D’Apice, P.I. Kogut, R. Manzo, M.V. Uvarov, On Variational Problem with Nonstandard Growth Conditions and Its Applications to Image Processing, Proceeding of the 19th International Conference of Numerical Analysis and Applied Mathematics, ICNAAM 2021, 20?26 September 2021, Location: Rhodes, Greece.
  * [16] R. Dautray, J.L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, Vol.5, Springer-Verlag, Berlin Heidelberg, 1985.
  * [17] L. Diening, P. Harjulehto, P. Hästö, M. Ru̇ẑiĉka, Lebesgue and Sobolev Spaces with Variable Exponents, Springer, New York, 2011.
  * [18] D. Frantz, _Landsat+ Sentinel-2 analysis ready data and beyond_ , Remote Sens., 11, 1124, 2019.
  * [19] F. Gao, J. Masek, M. Schwaller, F. Hall, _On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance_ , IEEE Tran. Geosci. Remote Sens., 44(8) (2006), 2207–2218.
  * [20] T. Hilker, M. A. Wulder, N. C. Coops, J. Linke, G. McDermid, J. G. Masek, F. Gao, J. C. White, _A new data fusion model for high spatial- and temporal-resolution mapping of forest disturbance based on Landsat and MODIS_ , Remote Sens. Environ., 113(8) (2009), 1613–1627.
  * [21] B.K. Horn, B.G. Schunck, Determining optical flow, Artificial Intelligence, 17 (1981), 185–203.
  * [22] T. Horsin, P. Kogut, Optimal $L^{2}$-control problem in coefficients for a linear elliptic equation. I. Existence result, Mathematical Control and Related Fields, 5 (1) (2015), 73–96.
  * [23] M.V. Joshi, K.P. Upla, Multi-resolution Image Fusion in Remote Sensing, Cambridge University Press, Cambridge, 2019.
  * [24] J. Ju, D. P. Roy, _The availability of cloud-free Landsat ETM+ date over the conterminous united states and globally_ , Remote Sens. Environ., 112(3) (2008), 1196–1211.
  * [25] P. Khanenko, P. Kogut, M. Uvarov, On Variational Problem with Nonstandard Growth Conditions for the Restoration of Clouds Corrupted Satellite Images, CEUR Workshop Proceedings, the 2nd International Workshop on Computational and Information Technologies for Risk-Informed Systems, CITRisk-2021, September 16-17, 2021, Kherson, Ukraine, Volume 3101, 6–25, 2021.
  * [26] P.I. Kogut, _On optimal and quasi-optimal controls in coefficients for multi-dimensional thermistor problem with mixed Dirichlet-Neumann boundary conditions_ , Control and Cybernetics, 48(1) (2019), 31–68.
  * [27] V.V. Hnatushenko, P.I. Kogut, M.V. Uvarov, _On flexible co-registration of optical and SAR satellite images_ , in ”Lecture Notes in ”Computational Intelligence and Decision Making” (series ’Advances in Intelligent Systems and Computing’), Springer, 2021, 515–534.
  * [28] V.V. Hnatushenko, P.I. Kogut, M.V. Uvarov, _Variational approach for rigid co-registration of optical/SAR satellite images in agricultural areas_ , Journal of Computational and Applied Mathematics, 400 (2022), Id 113742.
  * [29] P. Kogut, Ya. Kohut, R. Manzo, _Fictitious Controls and Approximation of an Optimal Control Problem for Perona-Malik Equation_ , Journal of Optimization, Differential Equations and Their Applications (JODEA), 30 (1) (2022), 42–70.
  * [30] P. Kogut, Ya. Kohut, N. Parfinovych, Solvability Issues for Some Noncoercive and Nonmonotone Parabolic Equations Arising in the Image Denoising Problems, Journal of Optimization, Differential Equations and Their Applications (JODEA), 30 (2) (2022), 42–70.
  * [31] P.I. Kogut, O.P. Kupenko, N.V. Uvarov, _On increasing of resolution of satellite images via their fusion with imagery at higher resolution_ , J. of Optimization, Differential Equations and Their Applications (JODEA), 29(1) (2021), 54–78.
  * [32] J.-L. Lions, Optimal Control of Systems Governed by Partial Differential Equations. Springer, Berlin, 1971.
  * [33] L. Loncan, L.B. De Almeida, J.V. Bioucas-Dias, X. Briottet, J. Chanussot, N. Dobigeon, S. Fabre, W. Liao, G.A. Licciardi, M. Simoes, J.Y. Tourneret, M.A. Veganzones, G. Vivone, Q. Wei, N. Yokoya, Hyperspectral pansharpening: a review, IEEE Geoscience and Remote Sensing Magazine, 3 (3) (2015), 27–46.
  * [34] J. G. Masek, E. F. Vermote, N. E. Saleous, R. Wolfe, F. G. Hall, F. Huemmrich, F. Gao, J. Kutler, T. K. Lim, _A Landsat surface reflectance data set for North American, 1990–2000_ , IEEE Geosci. Remote Sens. Lett., 3(1), (2006), 69–72.
  * [35] L. Nirenberg, _Topics in Nonlinear Analysis_ , Lecture Notes, New York University, New York, 1974.
  * [36] K. Rani, R. Sharma, _Study of Different Image fusion Algorithm_ , International Journal of Emerging Technology and Advanced Engineering, 3 (5) 2013, 288–290.
  * [37] D. P. Roy, J. Li, H. K. Zhang, L. Yan, _Best practices for the reprojection and resampling of Sentinel-2 Multi Spectral Instrument Level 1C data_ , Remote Sens. Lett., 7 (2016), 1023–1032.
  * [38] D. P. Roy, H. Huang, L. Boschetti, L. Giglio, H. K. Zhang, J. Li, _Landsat-8 and Sentinel-2 burned area mapping — a combined sensor multi-temporal change detection approach_ , Remote Sens. Environ., 231, 111254, 2019.
  * [39] P. Wang, F. Gao, J. G. Masek, _Operational data fusion framework for building frequent Landsat-like imagery_ , IEEE Transactions on Geoscience and Remote Sensing, 52(11), (2014), 7353–7365.
  * [40] L. Yan, D. P. Roy, H. Zhang, J. Li, H. Huang, _An automated approach for sub-pixel registration of Landsat-8 Operational Land Imager (OLI) and Sentinel-2 Multi Spectral Instrument (MSI) imagery_ , Remote Sens., 8, 520, 2016.
  * [41] V.V. Zhikov, Solvability of the three-dimensional thermistor problem, Proceedings of the Steklov Institute of Mathematics, 281 (2008), 98–111.
  * [42] V.V. Zhikov, On variational problems and nonlinear elliptic equations with nonstandard growth conditions, Journal of Mathematical Sciences, 173(5) (2011), 463–570.