A very serious obstacle to utilization of optical remote sensing satellite images is a risk of cloud and cloud shadow distortion issue (referred to as cloud contamination, hereafter). It has been reported in [1] that over 50% of all the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument aboard the Terra and Aqua satellites are covered by clouds or cloud-contaminated globally. Moreover, it is a typical situation when the measure of degradation of optical images is such that we can not rely even on the brightness inside of the damaged regions. As a result, there is a risk of information loss, some subdomains of such images become absolutely invisible. So, there is a great deal of missing information in optical remote sensing data, and a huge gap still exists between the satellite data we acquire and the data we require. Therefore, the reconstruction of missing information in remote sensing data becomes an active research field.

Many solutions have been developed to remove the clouds from multispectral images. (for the technical review, we refer to [2]). Formally, the present traditional algorithms can be primarily classified into four categories: 1) spatial-based methods, without any other auxiliary information source; 2) spectral-based methods, which extract the complementary information from other spectra; 3) multitemporal-based methods, which extract the complementary information from other data acquired at the same position and at different time periods; and 4) hybrid methods, which extract the complementary information by a combination of the three previous approaches. In parallel to the traditional approach, data-driven machine learning algorithms are actively developing since 2014 [3]. However, the necessity of large datasets and volatility to errors in input data limits its performance.

Our main effort in this research is to develop a new variational model for the exact restoration of the damaged multi-band optical satellite images that will meet demands from the agro application, that is, it must be applicable for large areas in different climate zones, and preserves the crop fields borders within damaged regions. In some sense, this model combines the spacialbased method with the multitemporal one. Therefore, in contrast to the standard variationalbased methods that are often optimased for a specific region or significantly blur textures, (see, for instance, [4] and the references therein), we focus on the global texture reconstruction inside damage regions. With that in mind we assume that the texture of a corrupted image can be predicted through a number of past cloud-free images of the same region from the time series. In order to describe the texture of background surface in the damage region, we follow the paper [5], where the authors experimentally checked the hypothesis that the essential geometric contents of a color image is contained in the level lines of the corresponding total spectral energy of such image.

We also pay much attention to the faithfulness of the reconstruction problem in the framework of the proposed model and supply this approach by the rigorous mathematical substantiation. 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.

Let Ω ⊂ R 2 be a bounded image domain with a Lipschitz boundary 𝜕Ω. With each particular image 𝑢 ⃗ = [𝑢 1 (𝑥), 𝑢 2 (𝑥), . . . , 𝑢 𝑀 (𝑥)] 𝑡 : Ω → R 𝑀 , where each coordinate represents the intensity of the corresponding spectral channel, we associate the panchromatic image 𝑢 (the so-called total spectral energy of 𝑢 ⃗ )

Here, 𝛼 1 , . . . , 𝛼 𝑀 are some weight coefficients.

Let 𝐷 ⊂ Ω be a Borel set with non empty interior and sufficiently regular boundary and such that |Ω ∖ 𝐷| > 0. We call 𝐷 the damage region. Let 𝑢 ⃗ 0 ∈ 𝐿 2 (Ω ∖ 𝐷; R 𝑀 ) be an image of interest which is assumed to be corrupted by clouds, and 𝐷 is the zone of missing information.

As it was mentioned before, we deal with the case where we have no information about the original image 𝑢 ⃗ 0 inside 𝐷. Instead of this, we assume that the texture of background surface in the damage region 𝐷 can be predicted with some accuracy by a number of past cloud-free images of the same region from the time series of satellite images. Unlike the wellknown 'chronochrome method' [6] which essentially assumes that the background in 𝐷 is stationary in wide sense, we admit that the image time series follows smooth variations over land (background), the time-series data are strictly chronological, and display regular fluctuations.

Let {𝑢 ⃗ 𝑡−1 , . . . , 𝑢 ⃗ 𝑡−𝑛 } be a given collection of past cloud-free images of the same region, where we set 𝑢 ⃗ 𝑡 = 𝑢 ⃗ 0 . We suppose that each cloud-free image of this time series should be well co-registered with 𝑢 ⃗ 0 ∈ 𝐿 2 (Ω ∖ 𝐷; R 𝑀 ) in Ω ∖ 𝐷 [7]. With each particular image 𝑢 ⃗ 𝑡−𝑘 in this series, we associate its total spectral energy 𝑢 𝑡−𝑘 using the standard rule (1). So, each element of the time series {𝑢 𝑡 , 𝑢 𝑡−1 , . . . , 𝑢 𝑡−𝑛 } is well-defined in Ω.

Let 𝑢 * be a predicted total spectral energy of 𝑢 ⃗ 0 in the damage region 𝐷. This prediction can be done following the regularized regression model and the available information in the time series {𝑢 𝑡 , 𝑢 𝑡−1 , . . . , 𝑢 𝑡−𝑛 } (for the details we refer to [8]).

where 𝜆 > 0 is the regularization parameter, 𝛽 1 > 0 and 𝛽 2 > 0 are the parameters that control the importance of the prediction and estimation terms, respectively. Seeing that the prediction and estimation errors in (2) are normalized by the volume of samples contributing to each term, we can constrain the values of 𝛽 1 = 𝛽 ∈ [0, 1] and 𝛽 2 = 1 − 𝛽 to control their relevance with a single parameter 𝛽.

As a result, setting 𝑤 0 = Argmin ℒ(𝑤), the total spectral energy 𝑢 * in the damage region 𝐷 can be estimated as follows 𝑢 * = ̂︀ 𝑢 in 𝐷, where

In order to reconstruct the texture (or geometry) of 𝑢 ⃗ 0 in the damage region 𝐷, we assume that predicted total energy 𝑢 * is a function of bounded variation, i.e. 𝑢 * ∈ 𝐵𝑉 (𝐷), and all spectral channels of the damaged image should share the geometry of the panchromatic image 𝑢 * ∈ 𝐿 2 (𝐷) in 𝐷. Hence, at most all points of almost all level sets of 𝑢 * ∈𝐵𝑉 (𝐷) we can define a normal vector 𝜃(𝑥), i.e., it formally satisfies (𝜃, 𝑢 * ) = |∇𝑢 * | and |𝜃| ≤ 1 a.e. in 𝐷.

In view of the risk of cloud distortion, the problem is to reconstruct the original multi-band image 𝑢 ⃗ 0 in the damage region 𝐷 using the knowledge of its texture (geometry) on the subset 𝐷 together with the exact information about this image in Ω ∖ 𝐷 (the undamaged region). We say that a function 𝑢 ⃗ = [𝑢 1 , 𝑢 2 , . . . , 𝑢 𝑀 ] 𝑡 : Ω → R 𝑀 is the result of restoration of a cloud corrupted image 𝑢 ⃗ 0 : 𝐷 → R 𝑀 if for given regularization parameters 𝜇>0, 𝛼>1, and 𝜆 𝑗 >0, 𝑗 = 1, 2, each spectral component 𝑢 𝑖 is the solution of the following constrained minimization problem with the nonstandard growth energy functional

• Ξ stands for the set of feasible solutions to the problem (3) which we define as follows

Here, 𝑊 By default we assume that the functions 𝑢 0,𝑖 and 𝑢 * are extended by zero outside of domains Ω ∖ 𝐷 and 𝐷, respectively.

The proposed model (3) provides a completely new approach to restoration of non-smooth multi-spectral images ⃗ 𝑢 0 with the gap in damage region. The main characteristic feature of this model is that we involve into consideration the energy functional with the nonstandard growth condition. The variable character of the exponent 𝑝(𝑥) provides more flexibility in terms of regularity for the recovered images. Since the first term in (3) is the regularization and the second one is the so-called data fidelity, it is worth to emphasize the role of the rest terms in (3). Taking into account the fact that the magnitude F(𝑣(𝑥)) is close to one at those points, where the spectral energy 𝑣 is slowly varying, and this value is close to zero at the edges of 𝑣, it follows that the edge information in the non-damage zone for the reconstruction is derived from the given image 𝑢 ⃗ 0 . So, in view of the estimate

the third term is also fidelity term which forces the texture (or topological map) of minimizer 𝑢 in domain Ω ∖ 𝐷 to stay close to the texture of a given spectral energy 𝑢 ⃗ 0,𝑖 . As for the last term in (3), we notice that since 𝜃 ∈ 𝐿 ∞ (𝐷, R 2 ) is a vector field with indicated properties, it follows that 𝜃(𝑥) has the direction of the normal to the level lines of 𝑢 * . Therefore, the counterclockwise rotation of angle 𝜋/2, denoted by 𝜃 ⊥ , represents the tangent vector to the level lines of 𝑢 * . In this case, if the spectral channels 𝑢 𝑖 share the geometry of the panchromatic image 𝑢 * , we have

Therefore, we impose them in the energy functional 𝐽 𝑖 in the form of the last term.

In practice, at the discrete level, 𝜃 can be defined by the relation

However, a better choice would be to compute 𝜃(𝑥 1 , 𝑥 2 ) by first regularizing ̂︀ 𝑢 by the equation

coupled with the initial and Neumann boundary conditions

Then, for any 𝑡 ≥ 0, there is a vector field 𝜉(𝑡) ∈ 𝐿 ∞ (Ω) with ‖𝜉(𝑡)‖ 𝐿 ∞ (Ω) ≤ 1 such that (see [9,10] for the details)

and

As a result, in order to characterize the texture of the cloud contaminated image ⃗ 𝑢 0 in the damage region 𝐷, we may take 𝜃 = 𝜉(𝑡) for some small value of 𝑡 > 0. As was mentioned in [9], following this way, we do not not distort the geometry of 𝑢 ︀ in an essential way.

In this section we show that constrained minimization problem (3) is consistent and admits at least one solution (𝑢 𝑟𝑒𝑐 𝑖 , 𝑝 𝑟𝑒𝑐 𝑖 ) ∈ Ξ, where 𝑝 𝑟𝑒𝑐 𝑖 (𝑥) = F(𝑢 𝑟𝑒𝑐 𝑖 (𝑥)) in Ω. We note that because of the specific form of the energy functional 𝐽 𝑖 (𝑣, 𝑝) in (3), the standard approaches are no longer applicable in its study, especially with respect to the existence of minimizers and their basic properties. It makes the minimization problem (3) rather challenging.

We begin with some auxiliary results which will play a crucial role in the sequel. Lemma 4.1. Let {𝑣 𝑘 } 𝑘∈N be a sequence of measurable non-negative functions 𝑣 𝑘 : Ω → [𝛾 0 , ∞) such that {𝑣 𝑘 } 𝑘∈N are uniformly bounded in 𝐿 1 (Ω) and 𝑣 𝑘 (𝑥) → 𝑣(𝑥) almost everywhere in Ω for some 𝑣 ∈ 𝐿 1 (Ω).

where

.

Proof. Since {𝑣 𝑘 } 𝑘∈N is the bounded sequence in 𝐿 1 (Ω), by smoothness of the Gaussian filter kernel 𝐺 𝜎 , it follows that

, ∀ 𝑥 ∈ Ω.

Then 𝐿 1 -boundedness of {𝑣 𝑘 } 𝑘∈N guarantees the existence of a positive value 𝛿 ∈ (0, 1) such that estimate (4) holds true for all 𝑘 ∈ N. Moreover, as follows from the estimate

and smoothness of the function ∇𝐺 𝜎 (•), there exists a positive constant 𝐶 𝐺 > 0 independent of 𝑘 such that

Setting

we see that

Since max 𝑥∈Ω |𝑝 𝑘 (𝑥)| ≤ 𝛽 and each element of the sequence {𝑝 𝑘 } 𝑘∈N has the same modulus of continuity, it follows that this sequence is uniformly bounded and equi-continuous. Hence, by Arzelà-Ascoli Theorem the sequence {𝑝 𝑘 } 𝑘∈N is relatively compact with respect to the strong topology of 𝐶(Ω). Taking into account that the set S is closed with respect to the uniform convergence and 𝑣 𝑘 (𝑥) → 𝑣(𝑥) almost everywhere in Ω, we deduce:

The proof is complete. Following in many aspects the resent studies [11,12], we give the following existence result. Theorem 4.2. For each 𝑖 = 1, . . . , 𝑀 and given 𝜇>0, 𝜆 1 >0, 𝜆 2 >0, 𝜃 ∈ 𝐿 ∞ (𝐷, R 2 ), and 𝑢 0,𝑖 ∈ 𝐿 2 (Ω ∖ 𝐷), the minimization problem (3) admits at least one solution (𝑢 𝑟𝑒𝑐 𝑖 , 𝑝 𝑟𝑒𝑐 𝑖 ) ∈ Ξ. Proof. Since Ξ ̸ = ∅ and 0 ≤ 𝐽 𝑖 (𝑣, 𝑝) < +∞ for all (𝑣, 𝑝) ∈ Ξ, it follows that there exists a non-negative value 𝜁 ≥ 0 such that 𝜁 = inf (𝑣,𝑝)∈Ξ 𝐽 𝑖 (𝑣, 𝑝). Let {(𝑣 𝑘 , 𝑝 𝑘 )} 𝑘∈N ⊂ Ξ be a minimizing sequence to the problem (3), i.e.

So, without lost of generality, we can suppose that 𝐽 𝑖 (𝑣 𝑘 , 𝑝 𝑘 ) ≤ 𝜁 + 1 for all 𝑘 ∈ N. From this and the initial assumptions, we deduce

where

Utilizing the fact that 𝑣 𝑘 (𝑥) ≤ 𝛾 1 for almost all 𝑥 ∈ Ω, we infer the following estimate

Then arguing as in Lemma 4.1 it can be shown that 𝑝 𝑘 ∈ 𝐶 0,1 (Ω) and

with

Taking this fact into account, we deduce from ( 6), (7), and (??) that

where

in Ω for all 𝑘 ∈ N, it follows from (9) that the limit function 𝑢 𝑟𝑒𝑐 𝑖 is also subjected the same restriction. Thus, 𝑢 𝑟𝑒𝑐 𝑖 is a feasible solution to minimization problem (3).

Let us show that (𝑢 𝑟𝑒𝑐 𝑖 , 𝑝 𝑟𝑒𝑐 𝑖 ) is a minimizer of this problem. With that in mind we note that in view of the obvious inequality

and the fact that 𝑢 0,𝑖 ∈ 𝐿 2 (Ω ∖ 𝐷), we have: the sequence {𝑣 𝑘 (𝑥) − 𝑢 0,𝑖 (𝑥)} 𝑘∈N is bounded in 𝐿 𝛼 (Ω ∖ 𝐷) and converges weakly in 𝐿 𝛼 (Ω ∖ 𝐷) to 𝑢 𝑟𝑒𝑐 𝑖 − 𝑢 0,𝑖 . Hence, by Proposition A.3 (see

As for the rest terms in (3), in view of the strong convergence

Hence,

As a result, utilizing relations ( 10), ( 11), (12), and the lower semicontinuity property (32), we finally obtain

Thus, (𝑢 𝑟𝑒𝑐 𝑖 , 𝑝 𝑟𝑒𝑐 𝑖 ) is a minimizer to the problem (3), whereas its uniqueness remains as an open question.

It is clear that because of the nonstandard energy functional and its non-convexity, constrained minimization problem (3) is not trivial in its practical implementation. The main difficulty in its study comes from the state constraints

that we impose on the set of feasible solutions Ξ. This motivates us to pass to some relaxation. In view of this, we propose the following iteration procedure which is based on the concept of relaxation of extremal problems and their variational convergence [13,14,15,16]. At the first step we set up

Then, for each 𝑘 ≥ 1, we set

Here,

where 𝐶 > 0 and 𝛿 > 0 are defined by ( 5) and ( 8), respectively.

Arguing as in the proof of Theorem 1 and using the convexity arguments, it can be shown that, for each 𝑝(•) ∈ S, there exists a unique element 𝑢

). Moreover, it can be shown that, for given 𝑖 = 1, . . . , 𝑀 , 𝜇 > 0, 𝜆 1 >0, 𝜆 2 > 0, 𝑢 * ∈ 𝑊 1,𝛼 (𝐷), and 𝑢 ⃗ 0 ∈ 𝐿 2 (Ω∖𝐷, R 𝑀 ), the sequence {︀ 𝑢 𝑘 ∈ 𝑊 1,𝑝 𝑘 (•) (Ω) }︀ 𝑘∈N is compact with respect to the weak topology of 𝑊 1,𝛼 (Ω), whereas the exponents {𝑝 𝑘 } 𝑘∈N are compact with respect to the strong topology of 𝐶(Ω).

We say that a pair (̂︀ 𝑢 𝑖 , ̂︀ 𝑝) is a weak solution to the original problem (3) if

Our main result can be stated as follows: Theorem 5.

and, in addition, the following variational property holds true

Proof. Let's assume the converse -namely, there is a function

Using the procedure of the direct smoothing, we set

where 𝜀 > 0 is a small parameter, 𝐾 is a positive compactly supported smooth function with properties

and

by the classical properties of smoothing operators (see [17]). From this we deduce that

Moreover, taking into account the estimates

we see that each element 𝑢 𝜀 is subjected to the pointwise constraints

Since, for each 𝜀 > 0, 𝑢 𝜀 ∈ 𝑊 1,𝑝 𝑘 (•) (Ω) for all 𝑘 ∈ N, it follows that 𝑢 𝜀 ∈ ℬ 𝑝 𝑘 (•) , i.e., each element of the sequence {𝑢 𝜀 } 𝜀>0 is a feasible solution to all approximating problems

Further we notice that lim inf

by Proposition A.3 and Fatou's lemma, and

Since

it follows from the Lebesgue dominated convergence theorem and (20) that

As a result, passing to the limit in (18) and utilizing properties ( 19)-( 21), we obtain

for all 𝜀 > 0. Taking into account the pointwise convergence (see (17) and property ( 16))

as 𝜀 → 0, and the fact that, for 𝜀 small enough,

we can pass to the limit in (22) as 𝜀 → 0 by the Lebesgue dominated convergence theorem. This yields

Combining this inequality with ( 22) and ( 15), we finally get

that leads us into conflict with the initial assumption. Thus,

and, therefore, (̃︀ 𝑢, ̃︀ 𝑝) is a weak solution to the original problem (3). As for the variational property (14), it is a direct consequence of (23) and (21).

To characterize the solution 𝑢 0,𝑝(•) ∈ ℬ 𝑝(•) of the approximating optimization problem

for all 𝑣 ∈ 𝑊 1,𝑝(•) (Ω), where

To this end, we note that

almost everywhere in Ω. Since, by convexity,

Taking into account that Utilizing the similar arguments to the rest terms in (3), we deduce that the representation (24) for the Gâteaux differential of 𝐽 𝑖 (•, 𝑝(•)) at the point 𝑢 0,𝑝(•) ∈ ℬ 𝑝(•) is valid.

Thus, in order to derive some optimality conditions for the minimizing element 𝑢 0,𝑝(•) ∈ ℬ 𝑝(•) to the problem inf

𝐽 𝑖 (𝑣, 𝑝(•)), we note that ℬ 𝑝(•) is a nonempty convex subset of 𝑊 1,𝑝(•) (Ω) and the objective functional 𝐽 𝑖 (•, 𝑝(•)) : ℬ 𝑝(•) → R is strictly convex. Hence, the well known classical result (see [18,Theorem 1.1.3]) and representation (24) lead us to the following conclusion. Theorem 6.1. Let 𝑝 𝑘 (•) ∈ S be an exponent given by the iterative rule (13). Then the unique minimizer 𝑢 𝑘 ∈ ℬ 𝑝 𝑘 (•) to the approximating problem inf 𝑣∈ℬ 𝑝 𝑘 (•) 𝐽 𝑖 (𝑣, 𝑝 𝑘 (•)) is characterized by

∇𝑢 𝑘 (𝑥), ∇𝑣(𝑥) − ∇𝑢 𝑘 (𝑥)

In order to illustrate the proposed algorithm for the restoration of satellite multi-spectral images we have provided some numerical experiments. As input data we have used a series of Sentinel-2 L2A images over the Dnipro Airport area, Ukraine (see Fig. 1, 2). This region represents a typical agricultural area with medium sides fields of various shapes. As a final result, we obtain in Fig. 3. Comparing the restored image and the contaminated one we could see that the texture of original image is well preserved. However, overall colors of different fields are shifted due to colorization part of algorithm. This problem has to be addressed in the following research.

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 (see [19] for the details), 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 𝑝(𝑥), where the variable exponent 𝑝(𝑥) 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 𝑢 ⃗ , we are going to restore, can have rather rich structure in the damage region 𝐷. 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.

Let 𝑝(•) be a measurable exponent function on Ω such that 1 ≤ 𝛼 ≤ 𝑝(𝑥) ≤ 𝛽 < ∞ a.e. in Ω, where 𝛼 and 𝛽 are given constants. Let 𝑝 ′ (•) = 𝑝(•) 𝑝(•)−1 be the corresponding conjugate exponent. It is clear that

where 𝛽 ′ and 𝛼 ′ stand for the conjugates of constant exponents. Denote by 𝐿 𝑝(•) (Ω) the set of all measurable functions 𝑓 (𝑥) on Ω such that ∫︀ Ω |𝑓 (𝑥)| 𝑝(𝑥) 𝑑𝑥 < ∞. Then 𝐿 𝑝(•) (Ω) is a reflexive separable Banach space with respect to the Luxemburg norm (see [20] for the details)

where 𝜌 𝑝 (𝑓 ) := ∫︀ Ω |𝑓 (𝑥)| 𝑝(𝑥) 𝑑𝑥. It is well-known that 𝐿 𝑝(•) (Ω) is reflexive provided 𝛼 > 1, and its dual is 𝐿 𝑝 ′ (•) (Ω), that is, any continuous functional 𝐹 = 𝐹 (𝑓 ) on 𝐿 𝑝(•) (Ω) has the form (see [21,Lemma 13.2])

As for the infimum in (26), we have the following result. Hence, (see [20] for the details) Let {𝑝 𝑘 } 𝑘∈N ⊂ 𝐶 0,𝛿 (Ω), with some 𝛿 ∈ (0, 1], be a given sequence of exponents. Hereinafter in this subsection we assume that 1 ≤ 𝛼 ≤ 𝑝 𝑘 (𝑥) ≤ 𝛽 < ∞ a.e. in Ω for 𝑘 = 1, 2, . . . , and 𝑝 𝑘 (•) → 𝑝(•) in 𝐶(Ω) as 𝑘 → ∞.

( (•) (Ω) to a function 𝑓 ∈ 𝐿 𝑝(•) (Ω), where 𝑝 ∈ 𝐶 0,𝛿 (Ω) is the limit of {𝑝 𝑘 } 𝑘∈N ⊂ 𝐶 0,𝛿 (Ω) in the uniform topology of 𝐶(Ω), if

For our further analysis, we need some auxiliary results (we refer to [21,Lemma 13.3] for comparison). We begin with the lower semicontinuity property of the variable 𝐿 𝑝 𝑘 (•) -norm with respect to the weak convergence in 𝐿 𝑝 𝑘 (•) (Ω).