In this paperwe discuss a newvariational 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 assumption to our approach is that the MODIS image has to be cloud-free whereas the images from Sentinel-2 can be corrupted by clouds or noise. We formulate the data fusion problem as the two-level optimization problem. We discuss the well thoroughness and consistency of the proposed variational models. We also derive some optimality conditions and supply our approach by results of numerical simulations with the real satellite images.
It is well-known that 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. 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 (22)). 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 — Setnitel-2 and Moderate Resolution Imaging Spectroradiometer (MODIS) [1]. 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 images from Sentinel-2, 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 MODIS data to generate synthetic ’daily’ surface reflectance products at Sentinel spatial resolution [2].
The problem we consider in this paper can be briefly described as follows. We have a collection of multi-band images { , ,… , : → ℝ } from Sentinel-2 that were captured at some time instances { , ,… , }, respectively, and we have a MODIS image : → ℝ from some day . 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 may does not coincide with any of time instances { , ,… , }. Meanwhile, the Sentinel images { , ,… , : → ℝ } 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 at the Sentinel-2 spatial resolution , utilizing for that the above mentioned data.
In principle, this problem is not new in the literature [3, 4, 5]. 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 [6, 7]. However, its performance essentially depends on the characteristic patch size of the landscape and degrades somewhat when used on extremely heterogeneous fine-grained landscapes [3].
Instead of this, we mainly focus on the variational approach to the satellite image data fusion (see [8, 9, 10]). 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 . 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 2–3 days, such prototype can be easily generated. We consider the above mentioned structural prototype as a reasonable input data for ’daily’ prediction problem that we formulate 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. It is worth to emphasize that our model is considerably different from the variational model for P+XS image fusion that was proposed in [11].
Let Ω ⊂ ℝ be a bounded connected open set with a sufficiently smooth boundary Ω and nonzero Lebesgue measure. In majority cases Ω can be interpreted as a rectangle domain. Let and be two sample grids on Ω .
Let [0, ] be a given time interval. Normally, by we mean a number of days. Let and { } be moments in time (particular days) such that 0 ≤ < < ⋯ < ≤ and < < . Let { , ,… , : → ℝ } be a collection of multispectral images of some territory, delivered from Sentinel-2, that were taken at time instances , ,… , , respectively. Let : → ℝ be a MODIS image of the same territory and this image has been captured at time = . It is assumed that: 1. The Sentinel-2 images { } can be corrupted by some noise, clouds and blur, whereas is a cloud-free image;
2. We divide the set of all bands for Sentinel images onto two parts and such that each spectral channel of the MODIS image has the similar spectral characteristics to some channel of -group;
3. The MODIS image : → ℝ is not corrupted by clouds or its damage zone can be neglected;
4. The MODIS image and the images { , ,… , } from Sentinel-2 are rigidly co-registered [12, 13]. This means that the MODIS image after arguably some affine transformation and each Sentinel images after the resampling to the grid with low resolution , could be successfully matched according to the unique geographic location [14, 15]. 2.1.
where ∥ ∥ ( ;ℝ )= ∫ | ( )| for 1 ≤
Given a real Banach space , we will denote by [0, ] into . We recall that, for 1 ≤ < ∞ , :[0, ]→ such that [16, 17] (Ω ;ℝ ) =
:Ω → ℝ / ≤ + ∞ , the space ∶ ∥ ∥ ( ;ℝ )< + ∞ ,
(Ω ;ℝ ) is defined by < + ∞ . ([0, ]; ) the space of all continuous functions from (0, ; ) is the space of all measurable functions ∥ ∥ ( , ; )= ∫ ∥ ( ) ∥ < ∞ , while (0, ; ) is the space of measurable functions such that ∥ ∥ ( , ; )= sup ∈[ , ] ∥ ( ) ∥ < ∞ .For more detailed presentation of the theory of these spaces and the Sobolev spaces with variable exponents, we refer to [18, 19, 20, 21, 22, 23]. 2.2.
Following the main principle of the Mathematical Morphology, a scalar image :Ω → ℝ is a representative of an equivalence class of images obtained from via a contrast change, i.e., = ( ), where is a continuous strictly increasing function. Under this assumption, a scalar image can be characterized by its upper level sets ( ) = { ∈ Ω ∶ ( ) ≥ }. Moreover, each image can be recovered from its level sets by the reconstruction formula ( ) = sup{ ∶ ∈ ( )}. Thus, according to the Mathematical Morphology Doctrine, we can suppose that the reliable geometric information about a scalar image is contained in those level sets.
In order to describe the level sets by their boundaries, ( ), we assume that ∈ , (Ω ), where , (Ω ) stands for the standard Sobolev space of all functions ∈ (Ω ) with respect to the norm ∥ ∥ , ( )= ∥ ∥ ( ) + ∥ ∇ ∥ ( ) . Then at almost all points of almost all level sets of ∈ , (Ω ) we may define a unit normal vector ( ) [24]. This vector field formally satisfies the following relations ( ,∇ ) = |∇ | and | |≤ 1 a.e.in. In the sequel, we will refer to the vector field as the vector field of unit normals to the topographic map of a function . So, we can associate with the geometry of the scalar image [25, 26].
Remark 2.1 In practice, a better choice for some small value of > 0, where ( , , ) is a solution (for small value of ( , ) would be to compute it as the ration ∇ ( ,⋅) for
|∇ ( ,⋅)| > 0) of the following initial-boundary value problem with 1 -Laplace operator in the principle part
∇ = div |∇ | , ∈ (0,+ ∞ ), ( , ) ∈ Ω , (0, , ) = ( , , )
( , ), ( , ) ∈ Ω ,
Definition 2.1 We say that a function :(0, )× → ℝ if it is defined by the rule :(0, )×
→ ℝ is the texture index of a given optical ( , ):= 1 + ∫ |(∇ ∗ ( ,⋅))( )| , for all ( , ) ∈ , where denotes zero extension of from = (0, )× Ω to ℝ , stands for the two-dimentional Gaussian of width (standard deviation) > 0 :[0,∞ ) → (0,∞ ) is the edge-stopping function which we take in the form of the Cauchy law ( ) = with > 0 small enough, and ℎ > 0 is a small positive value. ∥ ∥ ([ , ]; ( )) such that
Since ∈ (ℝ ) , it follows from (
Lemma 2.1 Let ∈ ([0, ]; ( )) be a measurable function extended by zero outside of . Let be the corresponding texture index. Then there exists a constant > 0 depending on , , and := 1 + ≤ ( , ) ≤ ∈ , ( := 2, ∀ ( , ) ∈ ), , where = ℎ ℎ+ ∥ ∥ ( ) |Ω |∥ ∥ ( , ; ( )) .
Let { , ,… , : → ℝ } be a collection of multispectral images of some territory from Sentinel-2 that were taken at time instances { , ,… , }⊂ [0, ], respectively. Let { , ,… , }⊂ 2 be a collection of damage regions for the corresponding Sentinel-images. So, in fact, we deal with the set of images { : \ → ℝ }, = 1,..., . Let : → ℝ be a MODIS image of the same territory and this image has been captured at time = ∈ ( , ).
We begin with the following assumption: = ∅ and the damage zones for the rest images from Sentinel-2 are such that each , = 2,… , , is a measurable closed subset of Ω with property ℒ ( ) ≤ 0.6 ℒ (Ω ), where ℒ ( ) stands for the 2-D Lebesgue measure of .
We say that multi-band images
cloud-corrupted ones { :
\
→ ℝ
if they are defined as follows:
are structural prototypes of the corresponding
(
, = [ [ , ] \ := ℒ ( \ )∑ ∈ \ , ( ).
Remark 3.1 As follows from the rule (
Since the MODIS image has been taken at a time instance
cases:
∈ ( , ), we can have three possible
(А1) there exists an index ∗ ∈ {1,..., } such that = ∗;
(А2) there exists an index ∗ ∈ {1,..., − 1} such that ∗ <
(А3) < < .
<
∗ ;
In view of this, we will distinguish three different statements of the data fusion problem:
(Restoration Problem) The problem (A1) consists in restoration of the damaged multi-band optical
image ∗: \ ∗ → ℝ using result of its fusion with the cloud-free MODIS image : → ℝ of the
same territory. It means that, we have to create a new image ∗ : → ℝ , which would be well defined
on the entire grid , such that
∗ ( ) =
∗( ), ∀
= ( , ) ∈
\ ∗,
(
The precise description of the class of admissible (or feasible) images ℐ will be given in the next section.
(Interpolation Problem) The problem (A2) consists in generation of a new multi-band optical image : → ℝ at the Sentinel-level of resolution using result of the fusion of cloud-free MODIS image : → ℝ with the predicted structural prototype : → ℝ from the given day . In fact, in this case we deal with the two-level problem. At the first level, having the collection of structural prototypes : → ℝ which is associated with the time instances { , ,… , }⊂ [0, ], we create a new ’intermediate’ image : → ℝ that can be considered as daily prediction of the topographical map of a given territory from the day . Then, at the second level, we realize the fusion procedure of this predicted image with the cloud-free MODIS image : → ℝ of the same territory.
(Extrapolation Problem) The problem (A3) consists in generation of a new multi-band optical image : → ℝ using result of the data assimilation from the cloud-free MODIS image : → ℝ into the structural prototype : → ℝ of the Sentinel-image : \ → ℝ .
Let ∗, and ∗ , be structural prototypes of the corresponding images from given days ∗ and ∗ . Since ∗, and ∗ , are well co-registered images, it is reasonable to assume that they have the similar geometric structure albeit they may have very different intensities.
The main question we are going to discuss in this section is: how to correctly define the ’intermediate’ image : → ℝ that can be considered as daily prediction of the topographical map of a given territory from the day . With that in mind, for each spectral channel ∈ {1,2,… , }, we make use of the following model − div|∇ | ( , ) ∇ =
in ( ∗, ∗ )× Ω , = 0 on ( ∗, ∗ )× Ω , ( ∗,⋅)= ∗, (⋅) in Ω , where ( , ) stands for the texture index of the scalar image (see Definition 2.1), and ∈ ( ∗, ∗ ; (Ω )) is an unknown source term that should be defined in the way to guarantee the fulfillment (with some accuracy) of the relation
( ∗ ,⋅)≈ ∗ , (⋅) in Ω .
Here, in (
∗ ∫ ∗ ∫
− holds true for any function + (|∇ | ∇ ,∇ )
∗ = ∫ ∗
∫ ∈ Φ, where Φ = + ∫ ∈ ∗, |
∗ ([ ∗, ∗ ]× Ω ) ∶ | ∗ = 0 .
(
Utilizing the perturbation technique and the classical fixed point theorem of Schauder, it has been recently proven the following existence result.
Remark 4.1 The main characteristic feature of the proposed initial-boundary value problem (IBVP) is
the fact that the exponent depend not only on ( , ) but also on ( , ). It is well-known that the
variable character of exponent causes a gap between the monotonicity and coercivity conditions.
Because of this gap, equations of the type (
∗, (⋅) and
∗ , (⋅).
We note that the distributed control in the right hand side of (
Definition 4.1 We say that, for given
to the problem (
∈ ( ) and ∈ [ ∗, ∗ ], ∗, ∈ ( ), a function is a weak solution ( ∗, ∗ ; (Ω )), ( ,⋅)∈
∗ , (Ω ) , ∫ ∗ ∫ |∇ | ( , ) < + ∞ ,
Theorem 4.1 [28] Let ∈ ( ) and ∗, ∈ ( ) be given distributions. Then initial-boundary
value problem (
∈ ( ∗, ∗ ; (Ω )), ∈ , (( ∗, ∗ )× Ω ), ∈ ( ∗, ∗ ; (Ω )),
where the exponent is given by the rule
(
In order to satisfy the condition (
∗, , div(|∇ | ∇ )| = [div |∇ ∗, | ∗, ∇ ∗, + div |∇ ∗ , | ∗ , ∇ ∗ , ].
It is clear that a minimum point ∗ ∈ (Ω ) to unconstrained minimization problem (
Setting = ∗ in (
This observation leads us to the following conclusion: the ’intermediate’ image defined as follow: : → ℝ can be ( ) = ∗( , ), ∀ ∈ .
Coming back to the principle cases (A1)–(A3), that have been described in Section 3, we can suppose
that a structural prototype from the given day is well defined. As it was pointed out in Section 3,
this prototype coincides either with one of the images in cases (A1) and (A3), or it should be defined
using the solutions of the problem (
Let ∈ {1,..., } be a fixed index value (the number of spectral channel). Let :Ω → ℝ be the
texture index of the -th band for the structural prototype : → ℝ . Let ∈ (
We define the linear operator ( ),∇
, ( ) ℝ , :ℝ → ℝ , ∇ := ∇ −
= |∇ as follows ,∇
, = ∫ \ ∗ | ∗, ( )| − ∫ \ ∗ | , ( )| is a solutions of the following constrained minimization problem ℱ ( ) = ∫ ( ) | , ∇ ( )| ( ) + ∫ ∇ ( )− ∇ , ( ) (22) (23) + ∫ Π ( ∗ − , ) → i∈nf, Ξ = ∈ , (⋅)(Ω ) ∶ 0 ≤ ( ) ≤ a.e.in Ω stands for the set of feasible solutions, and , (⋅)(Ω ) denotes the Sobolev space with variable exponent. As for the constants , their choice depends on the format of signed integer numbers in which the corresponding intensities , ( ) are represented. In particular, it can be = 2 − 1, = 2 − 1, and so on.
Problem, if it is given by the rule : → ℝ , with ∗ < < ∗ , is a solution of the Interpolation , ∀ ∈ = ∫ | , ( )| − ∫ | , ( )| is a solutions of the constrained minimization problem (22).
and where and
Problem, if it is given by the rule : → ℝ , with < <
, is a solution of the Extrapolation
, ( ) , ∀ ∈ , is a solutions of the constrained minimization problem (22), and ∀ ∈ , ∀ ∈ ,
Let us briefly discuss the relevance of the proposed minimization problem. The first term in (22) can be considered as the regularization in the Sobolev-Orlicz space , (⋅)(Ω ) [29]. As for the second term in (22), 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 : → ℝ . We interpret this closedness in its simplified form, namely, in the sense of -norm of the difference of the corresponding gradients [6, 30, 31]. The last term in (22) represents an -distortion between a -th spectral channel in the MODIS image and the corresponding channel of the retrieved image which is resampled to the grid of low resolution
Existence Result and Optimality Conditions for Constrained Minimization Problem ( )
Theorem 6.1 Let : Sentinel-2. Then for any given a unique solution ∈ .
→ ℝ
∈ ,
be a given structural prototype for unknown image from
> 0, > 0, and ∈ (
Theorem 6.2 Let : → ℝ be a given structural prototype for unknown image from Sentinel-2. Let : → ℝ be a given MODIS image. Let stands for the texture index of the -th band for the predicted structural prototype . Then the unique minimizer ∈ to the minimization problem ℱ ( ) is characterized by the following variational inequality ∈ ∫ | , ∇
( )| ( ) + ∫ Π + ∫ ∗ ∗
, ∇ ∇
( ), , ∇ ( )− , ∇ ( )− ∇ , ( ),∇ ( ) ∗ − , ≥ 0, ∀ ( ) ∈ Ξ .
Remark 6.1 In practical implementation, it is reasonable to define an optimal solution ∈ using a ’gradient descent’ strategy. Indeed, following the standard procedure and starting from the initial image , , we can pass to the following initial value problem for the quasi-linear parabolic equations with Nuemann boundary conditions
− div | , ∇ = − div | , ∇ ( )| ( ) ( )| ( ) , ∇
( ) , ∇ + div∇ ( )− ∇ ( )| ( ) , ( ) − ,
Π = 0 on ( ), ∗ ∗ ∗ − | , ∇ , ∇ Ω , 0 ≤ ( ) ≤ a.a.in Ω ,
(0, ) = , ( ), ∀ ∈ Ω . (24) In principle, instead of the initial condition (24) we may consider other image that can be generated from , and the bicubic interpolation of the MODIS band , onto the entire domain Ω following one of well-known simple data fusion methods [32].
In order to illustrate the proposed approach for the restoration of satellite multi-spectral images we have used two Sentinel-2 images (610 × 699 in pixels) over the Dnipro Airport area, Ukraine (see Fig. 1) that were captured at different time instances. What is the most important, it was a growing season when the global biophysical processes are rapid enough. This region represents a typical agricultural area with medium sides fields of various shapes. Figure
2: MODIS image with the date of generation 2019/07/01
We also have a cloud-free MODIS image (67 × 77 in pixels) from 2019/07/01 (see Fig. 2). We solve the interpolation problem (A2). As a result, a new image from the date 2019/07/01 at the Sentinel-level of resolution is depicted in Fig. 3. These numerical simulations have been supplied by a detailed analysis of the obtained results using the special validation metrics. 8. References [19] D. Cruz-Uribe, A. Fiorenza, Variable Lebesgue Spaces: Foundations and Harmonic Analysis,
[20] L. Diening, P. Harjulehto, P. Hästö, M. R ̇ uˆ ziˆ cka, Lebesgue and Sobolev Spaces with Variable
[21] L. Nirenberg, Topics in Nonlinear Analysis, Lecture Notes, New York University, 1974. [22] V. Zhikov, Solvability of the three-dimensional thermistor problem, Proceedings of the Steklov
[23] V. Zhikov, On variational problems and nonlinear elliptic equations with nonstandard growth
conditions, Journal of Mathematical Sciences 173 (
Transactions on Image Processing 10 (
Access 8 (2020) 222944–222955. [31] P. Blomgren, T. Chan, P. Mulet, C. Wong, Total variation image restoration: Numerical methods and extensions, Proceedings of the 1997 IEEE International Conference on Image Processing III (1997) 384–387. [32] K. Rani, R. Sharma, Study of different image fusion algorithm, International Journal of Emerging