Variational model with nonstandard growth conditions for restoration of satellite optical images using synthetic aperture radar

Abstract

In this paper, the problem of restoration of cloud contaminated optical images is studied in the case when we have no information about brightness of such images in the damage region. We propose a new variational approach for exact restoration of optical multi-band images utilising Synthetic Aperture Radar (EOS – Spatial Data Analytics, GIS Software, Satellite Imagery – is a cloud-based platform to derive remote sensing data and analyse satellite imagery for business and science purposes) images of the same regions. We prove existence of solutions, propose an alternating minimisation method for computing them, prove convergence of this method to weak solutions of the original problem and derive optimality conditions.

1 Introduction

It is well-known that optical satellite multi-band images have a high resolution and can be easily captured by low-cost cameras. However, they are often corrupted because of poor weather conditions, such as rain, clouds, fog and dust conditions. Moreover, it is a typical situation that the measure of degradation of optical images is such that it cannot even rely one brightness value being available inside the damaged regions. As a result, some subdomains of such images become absolutely invisible. Thus, in spite of the fact that in the literature there are many approaches to the reconstruction of images when information of the colours is not everywhere available (see, for instance, [23, 31, 37, 38, 40]), the traditional approaches to the exact restoration of damaged optical images are no longer applicable in this case and it makes this problem challenging.

However, in contrast to optical observations, radar images do not depend on reflected sunlight and they can be used at night and under poor weather conditions. In the vegetation case, instead of giving an indication on biophysical processes in the plant, the radar backscatter rather contains information on the structure and moisture content of vegetation and the underlying soil. Therefore, the fusion of Synthetic Aperture Radar (SAR) and optical images is very important for classification of land cover [34] and estimation of soil moisture to remove vegetation cover effects from radar backscattering coefficient [19, 20, 36]. At the same time, because of the distinct natures of SAR and optical images, there exist a huge radiometric difference between optical and synthetic aperture of radar images.

It is well-known that different wavelengths encode different object properties, therefore leading to significant intensity differences between SAR and optical satellite images for the same object. Because of this, it would be naive to suppose that the intensities in all spectral channels of cloud contaminated optical images can be successfully recovered inside the damaged regions at high level of accuracy through the corresponding intensities of the SAR images of this region. At the same time, when dealing with optical and SAR images of some agricultural areas of various shapes or/and farmland, it can be observed that the textures of such images have many common features (see, for instance, [22]). Mainly inspired by this observation, we propose a new variational model for exact restoration of the damaged multi-band optical satellite images using results of their co-registration with SAR images of the same regions.

Let $\Omega\subset\mathbb{R}^2$ be a bounded image domain with Lipschitz boundary $\partial\Omega$ , and let $D\subset\Omega$ be a Borel set with non-empty interior and sufficiently regular boundary and such that $|\Omega\setminus D|>0$ . We call D the damage region of a given multi-band image $\vec{u}_0=\left[u_{1,0}, u_{2,0},\dots,u_{M,0}\right]^t\in L^2(\Omega\setminus D;\mathbb{R}^M)$ where the optical image $\vec{u}_0$ is corrupted by clouds. As it was mentioned before, we deal with the case where we have no information about the original image $\vec{u}_0$ inside D. Instead of this, we assume that a SAR image $u_{SAR}:\Omega\to\mathbb{R}$ of the same region is given, and this image is well co-registered with $\vec{u}_0$ in $\Omega\setminus D$ . The problem is to reconstruct the intensities $u_{i,0}(x)$ , $i=1,\dots,M$ , of the original multi-band image $\vec{u}_0$ through a variational model starting from the knowledge of SAR image on the subset D (the damaged region) together with the exact information of $\vec{u}_0$ on $\Omega\setminus D$ (the undamaged region).

Mostly motivated by the recent studies in this field [4] (see also [6, 7, 8, 9, 11] for comparison), we propose a new strategy for the restoration problem of cloud contaminated multi-band optical images through their fusion with the corresponding SAR images of the same territory. The variational approach we consider is inspired, in some sense, by the recent paper [4], where the authors consider the minimisation of the following energy functional

(1.1) \begin{equation}u\mapsto \frac{1}{2}\|\left(-\Delta \right)^\frac{s}{2} u\|^2_{L^2(\Omega)}+\frac{\alpha}{2}\|\left(-\Delta \right)^{-\frac{\beta}{2}}(u-g)\|^2_{L^2(\Omega)}\end{equation}

with $0<s<1$ and $\beta\in[0,1]$ , where $\left(-\Delta \right)^s$ denotes the fractional power of the Laplacian with zero Neumann boundary conditions. Then, the first necessary and sufficient optimality condition determines the unique minimiser u via

\begin{gather*}\left(-\Delta \right)^s u +\alpha \left(-\Delta \right)^{-\beta}(u-g)=0\quad \text{in }\ \Omega,\\\partial_\nu u=0\quad\text{on }\ \partial\Omega,\end{gather*}

which is a linear elliptic partial differential equation that can be efficiently solved using, for instance, the Fourier spectral method [3] or the Stinga-Torrea extension [35]. Since, from the reconstruction point of view, it is desirable that the regularity of the solution to (1.1) is low in places in $\Omega$ where edges or discontinuities are present in u, and that is high in places where u is smooth or contains homogeneous features, it is of interest to consider (1.1) where $s :\Omega \to [0, 1]$ is not a constant. So, the choice of parameter s has a direct influence on the global regularity of the solution to problem (1.1). However, the definition of the fractional Laplacian in terms of the Caffarelli–Silvestre [10] or the Stinga–Torrea [39] extension is well-known only for a constant $s\in(0,1)$ , whereas such a result remains open when $s(x)\in (0,1)$ for $x\in\Omega$ . So, there is no obvious way how to correctly define the operator $\left(-\Delta \right)^{s(x)}$ and for nowadays it looks as an open question. For the details, we refer to [4].

In contrast to the standard setting of restoring colour damaged images where the starting point is either the knowledge of the grey level of the original colour image $\vec{u}_0$ on a given open subset D of $\Omega$ (the damaged region) together with the exact information of $\vec{u}_0$ on $\Omega\setminus D$ (the undamaged region) or the grey level information in the damage region $D\subset \Omega$ is modelled as a nonlinear distortion of the colours, in this paper we deal with the case where we do not have any information about $\vec{u}_0$ inside D but instead we assume that a SAR image $u_{SAR}:\Omega\to\mathbb{R}$ of the same region is given. So, the purpose of this paper is to study the faithfulness of the reconstruction following the proposed variational model and supply this approach by the rigorous mathematical substantiation.

The paper is structured as follows. Section 2 contains some preliminaries and auxiliary results. For other details related to Orlicz and variable exponent Sobolev spaces, we refer to Appendices A–C. In Section 3, we begin with some key assumptions and after that we give a precise statement of the restoration problem in the form of some constrained minimisation problems with nonstandard growth energy functionals. We discuss the consistency issues of the proposed minimisation problems and the existence of the corresponding minimisers in Section 4. To prove the existence result, mainly because of the specific form of the objective functional, we follow the technique that was recently developed by the authors in [14, 15, 24].

In spite of the fact that the proposed minimisation problems are well-posed and possess good approximating properties, their practical implementation to the restoration of real satellite images is a tricky matter because of the non-convexity and complicated structure of the corresponding optimality conditions. Therefore, our main intention in Section 5 is to present ‘an approximation approach’, which is based on the concept of relaxation of extremal problems and their variational convergence. With that in mind, we propose to pass to some relaxation of the original problem using a special iterative algorithm. We show that at each step of the iteration procedure, we obtain a strictly convex optimisation problem which admits a unique solution. It is established that the so-defined sequence of approximations is precompact in some Hausdorff topology, and each convergent subsequence leads to a weak solution of the original problem. Optimality conditions for the relaxed version of a given restoration problem and their substantiation are studied in Section 6. The experiments undertaken in this study (see Section 7) 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.

2 Preliminaries and some auxiliary results

We begin with some notation. For vectors $\xi\in\mathbb{R}^N$ and $\eta\in\mathbb{R}^N$ , $\left(\xi,\eta\right)=\xi^t\eta$ denotes the standard vector inner product in $\mathbb{R}^N$ , where $^t$ denotes 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

\begin{equation*}\chi_E(x):=\left\{\begin{array}{ll}1,&\ \text{for }\ x\in E,\\ \\[-7pt] 0,&\ \text{otherwise}.\end{array}\right.\end{equation*}

Let $D\subset\Omega$ be a Borel set with non-empty interior and sufficiently regular boundary and such that $|\Omega\setminus D|>0$ . Let $\vec{u}_0=\left[u_{1,0}, u_{2,0},\dots,u_{M,0}\right]^t\in L^2(\Omega\setminus D;\mathbb{R}^M)$ be an image of interest, where each coordinate represents the intensity of the corresponding spectral channel. We associate with $\vec{u}_0$ the panchromatic image u (the so-called total spectral energy of $\vec{u}$ )

(2.1) \begin{equation}u(x)=\alpha_1 u_1(x)+\dots \alpha_M u_M(x).\end{equation}

In particular, if we take into account in (2.1) just the weight coefficients for the RGB-channels with

\begin{equation*}\alpha_R=0.299,\quad \alpha_G=0.587,\quad \alpha_B=0.114,\end{equation*}

then u can be interpreted as the luma component of $\vec{u}$ and it represents the perceptual brightness of the multispectral image $\vec{u}:\Omega\setminus D\rightarrow \mathbb{R}^M$ .

We assume that the multi-band image $\vec{u}_0$ is corrupted by clouds, and D is the zone of missing information. We call D the damage region. Since we have no information about the original image $u_0$ inside D, we assume that a SAR image $u_{SAR}:\Omega\to\mathbb{R}$ of the same region is given, and this image is well co-registered with $u_0\in L^2(\Omega\setminus D;\mathbb{R}^M)$ in $\Omega\setminus D$ . This means that there exists an affine transformation $\mathcal{F}: \mathbb{R}^2\rightarrow \mathbb{R}^2$ of the form

(2.2) \begin{equation}\mathcal{F}(x)=Bx+a,\quad\forall\,x\in \mathbb{R}^2,\end{equation}

where

\begin{equation*}a=\left[\begin{array}{c}a_1\\ \\[-7pt] a_2\end{array}\right]\quad\text{ and }\quad =\left[\begin{array}{c@{\quad}c}b_{11} & b_{12}\\ \\[-7pt] b_{21} & b_{22}\end{array}\right]\end{equation*}

such that the SAR image after the affine transformation $u_{SAR}\left(\mathcal{F}^{-1}(\cdot)\right)$ and panchromatic image $u(\cdot)$ could be successfully matched in $\Omega\setminus D$ (see [21] for the details).

We assume that the SAR image $u_{SAR}:\Omega\to\mathbb{R}$ is a function of bounded variation, that is, $u_{SAR}\in BV(\Omega)$ . Then, $u_{SAR}\in L^2(\Omega)$ and almost all level sets $\left\{x\in\Omega: u_{SAR}(x)\ge \lambda\right\}$ are sets of finite perimeter. Hence, at almost all points of almost all level sets of $u_{SAR}\in BV(\Omega)$ we can define a normal vector $\theta(x)$ . This vector field of normals $\theta(x)$ can be also defined as the Radon-Nikodym derivative of the measure $\nabla u_{SAR}$ with respect to $|\nabla u_{SAR}|$ , that is, it formally satisfies the following relations

\begin{equation*}\left(\theta, \nabla u_{SAR}\right) = |\nabla u_{SAR}|\quad\text{and}\quad|\theta|\le 1 \ \text{a.e. in $\Omega$}.\end{equation*}

In the sequel, we will refer to the vector field $\theta$ as the vector field of unit normals to the topographic map of the image $u_{SAR}$ .

Remark 2.1 In practice, at the discrete level, $\theta(x,y)$ can be defined by the rule $\theta(x_i,y_j)=\frac{\nabla u_{SAR}(x_i,y_j)}{|\nabla u_{SAR}(x_i,y_j)|} $ when $\nabla u_{SAR}(x_i,y_j)\ne 0$ , and $\theta=0$ when $\nabla u_{SAR}(x_i,y_j)=0$ . However, as was mentioned in [5], a better choice for $\theta(x,y)$ would be to compute it as $\theta=\xi(t)$ for some small value of $t > 0$ , where $\xi(t)=\frac{\nabla U(t,\cdot)}{|\nabla U(t,\cdot)|}$ and U(t,x,y) is a solution of the following initial-boundary value problem with 1-Laplacian in the right hand side

(2.3) $$\matrix{ {{{\partial U} \over {\partial t}} = {\rm{div}}{\mkern 1mu} \left( {{{\nabla U} \over {|\nabla U|}}} \right),\quad t \in (0, + \infty ),\;(x,y) \in \Omega ,} \hfill\cr}$$

(2.4) \begin{align} &U(0,x,y)= u_{SAR}(x,y),\quad (x,y)\in\Omega,\end{align}

(2.5) \begin{align} &\frac{\partial U(0,x,y)}{\partial\nu}=0,\quad t\in(0,+\infty),\ (x,y)\in\partial\Omega. \end{align}

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

\begin{equation*} \xi\in L^\infty(\Omega;\mathbb{R}^2)\ \text{ with }\ \|\xi(t)\|_{L^\infty(\Omega;\mathbb{R}^2)}\le 1 \end{equation*}

such that

(2.6) \begin{equation} \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, \end{equation}

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

We notice that following this procedure, for small value of $t > 0$ , we do not distort the geometry of the function $u_{SAR}(x,y)$ in an essential way. Moreover, it can be shown that this regularisation of the vector field $\theta$ satisfies condition ${div}\,\theta \in L^2(\Omega)$ .

Since the main problem we are going to consider in this article is to reconstruct the original multi-band image $\vec{u}_0\in L^2(\Omega\setminus D;\mathbb{R}^M)$ in the damage region D using the knowledge of geometrical texture of the SAR image $u_{SAR}\in BV(\Omega)$ on the subset D together with the exact information about this image in $\Omega\setminus D$ (the undamaged region), we associate with each spectral channel of the restored optical image $\vec{u}=\left[u_1, u_2,\dots,u_M\right]^t:\Omega\rightarrow \mathbb{R}^M$ the so-called texture index $p_i:\Omega\rightarrow\mathbb{R}$ following the rule

(2.7) \begin{equation}p_i(x):=\mathfrak{F}(u_i(x))=1+g\left(\left|\left(\nabla G_\sigma\ast u_i\right) (x)\right|\right),\quad\forall\,x\in\Omega,\ \forall\,i=1,\dots,M,\end{equation}

where $g{:}[0,\infty)\rightarrow (0,\infty)$ is an edge-stopping function which we take in the form of the Cauchy law $g(t)=\frac{1}{1+(t/a)^2}$ with an appropriate $a>0$ , $\left(\nabla G_{\sigma}\ast u_i\right) (x)$ denotes the convolution of function $u_i$ with the two-dimensional Gaussian filter kernel $G_\sigma$ ,

(2.8) \begin{gather}G_{\sigma}(x)=\frac{1}{2\pi \sigma^2}e^{-\frac{|x|^2}{2\sigma^2}},\quad x\in\mathbb{R}^2,\end{gather}

(2.9) \begin{align}\left(\nabla G_{\sigma}\ast u_i\right)(x):=\int_{\Omega} \nabla G_{\sigma}(x-y) u_i(y)\,dy,\quad\forall\,x\in\Omega.\\ \nonumber\end{align}

Here, the parameter $\sigma{>}0$ determines the spatial size of the image details which are removed by this 2D filter. By default, we assume that the functions $u_i$ are extended by zero outside of domain $\Omega$ .

As follows from (2.7), the magnitude $g\left(\left|\left(\nabla G_\sigma\ast u_i\right) (x)\right|\right)$ is close to one at those points, where the spectral intensity $u_i$ is slowly varying, and this value is close to zero at the edges of $u_i$ . In view of this, the function $p_i(x)$ can be interpreted as a characteristic of texture of the image $\vec{u}$ in its i-th spectral channel. The following result plays a crucial role in the sequel.

Lemma 2.1 Let $\left\{v_k\right\}_{k\in\mathbb{N}}\subset L^\infty(\Omega)$ be a sequence of measurable non-negative functions such that $v_k \rightarrow v $ weakly- $\ast$ in $L^\infty(\Omega)$ for some $v\in L^\infty(\Omega)$ , and each element of this sequence is extended by zero outside of $\Omega$ . Let $\left\{p_k=1+g\left(\left|\left(\nabla G_\sigma\ast v_k\right)\right|\right)\right\}_{k\in\mathbb{N}}$ be the corresponding sequence of texture indices. Then,

(2.10) \begin{align} \notag p_k(\cdot)\rightarrow p(\cdot)=1+g\left(\left|\left(\nabla G_\sigma\ast v\right)(\cdot)\right|\right)\quad\text{uniformly in $\overline{\Omega}$ as $k\to\infty$},\\[3pt] \alpha:=1+\delta\le p_k(x) \le \beta:=2,\quad \forall\,x\in\Omega,\ \forall\,k\in\mathbb{N}, \end{align}

(2.11) \begin{align} \text{with }\quad\delta={a^2}\left[{a^2+\|G_{\sigma}\|^2_{C^1(\overline{\Omega-\Omega})} \sup_{k\in \mathbb{N}}\|v_k\|^2_{L^1(\Omega)}}\right]^{-1}. \end{align}

Proof. In view of the initial assumptions, the sequence $\left\{v_k\right\}_{k\in\mathbb{N}}$ is uniformly bounded in $L^1(\Omega)$ . Hence, by smoothness of the Gaussian filter kernel $G_\sigma$ , it follows that

\begin{align*} \left|\left(\nabla G_\sigma\ast v_k\right) (x)\right|&\le \int_{\Omega} \left|\nabla G_{\sigma}(x-y)\right| v_k(y)\,dy\le \|G_{\sigma}\|_{C^1(\overline{\Omega-\Omega})} \|v_k\|_{L^1(\Omega)},\\[3pt] p_k(x)&=1+\frac{a^2}{a^2+(\left|\left(\nabla G_\sigma\ast v_k\right) (x)\right|)^2}\\[3pt] &\ge 1+\frac{a^2}{a^2+\|G_{\sigma}\|^2_{C^1(\overline{\Omega-\Omega})} \|v_k\|^2_{L^1(\Omega)}},\quad\forall\,x\in\Omega, \end{align*}

where $\|G_{\sigma}\|_{C^1(\overline{\Omega-\Omega})}=\max\limits_{z=x-y\atop x\in\overline{\Omega}, y\in\overline{\Omega}} \left[|G_\sigma(z)|+|\nabla G_\sigma(z)|\right]$ . Then, $L^1$ -boundedness of $\left\{v_k\right\}_{k\in\mathbb{N}}$ guarantees the existence of a positive value $\delta\in (0,1)$ (see (2.11)) such that estimate (2.10) holds true for all $k\in\mathbb{N}$ .

Moreover, as follows from the relations

(2.12) \begin{align} \notag \left|p_k(x)-p_k(y)\right|&\le a^2\left| \frac{\left|\left(\nabla G_\sigma\ast v_k\right) (x)\right|^2-\left|\left(\nabla G_\sigma\ast v_k\right) (y)\right|^2}{\left(a^2+\left|\left(\nabla G_\sigma\ast v_k\right) (x)\right|^2\right)\left(a^2+\left|\left(\nabla G_\sigma\ast v_k\right) (y)\right|^2\right)}\right|\\[3pt] &\le \frac{2\|G_{\sigma}\|_{C^1(\overline{\Omega-\Omega})} \|v_k\|_{L^1(\Omega)}}{a^2} \left|\left|\left(\nabla G_\sigma\ast v_k\right) (x)\right|-\left|\left(\nabla G_\sigma\ast v_k\right) (y)\right|\right| \end{align}

\begin{align} \notag &\le \frac{2\|G_{\sigma}\|_{C^1(\overline{\Omega-\Omega})} \gamma_1^2|\Omega|}{a^2} \int_{\Omega}\left|\nabla G_\sigma(x-z)-\nabla G_\sigma(y-z)\right|\,dz,\\[3pt] \notag &\qquad\forall\,x,y\in\Omega\quad\text{with }\ \gamma_1=\sup_{k\in \mathbb{N}} \|v_k\|_{L^\infty(\Omega)}, \end{align}

and smoothness of the function $\nabla G_\sigma(\cdot)$ , there exists a positive constant $C_G>0$ independent of k such that

\begin{equation*} \left|p_k(x)-p_k(y)\right|\le \frac{2\|G_{\sigma}\|_{C^1(\overline{\Omega-\Omega})} \gamma_1^2|\Omega|C_G}{a^2}|x-y|, \ \forall\,x,y\in\Omega. \end{equation*}

Setting

(2.13) \begin{equation} C:=\frac{2\|G_{\sigma}\|_{C^1(\overline{\Omega-\Omega})} \gamma_1^2|\Omega|C_G}{a^2}, \end{equation}

we see that

\begin{equation*} \left\{p_k(\cdot)\right\}\subset \mathfrak{S}=\left\{h\in C^{0,1}(\Omega)\ \left| \begin{array}{c} |h(x)-h(y)|\le C|x-y|,\ \forall\,x,y,\in \Omega,\\ \\[-7pt] 1<\alpha\le h(\cdot)\le\beta\ \text{in}\ \overline{\Omega}. \end{array} \right. \right\} \end{equation*}

Since $\max_{x\in\overline{\Omega}}|p_{k}(x)|\le\beta$ and each element of the sequence $\left\{p_k\right\}_{k\in\mathbb{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 $\left\{p_{k}\right\}_{k\in \mathbb{N}}$ is relatively compact with respect to the strong topology of $C(\overline{\Omega})$ . Taking into account the estimate (2.12) and the fact that the set $\mathfrak{S}$ is closed with respect to the uniform convergence and

\begin{equation*} \left(\nabla G_\sigma\ast v_k\right)(x)\rightarrow \left(\nabla G_\sigma\ast v\right)(x)\ \text{ as $k\to\infty$},\ \forall\, x\in \Omega \end{equation*}

by definition of the weak- $\ast$ convergence in $L^\infty(\Omega)$ , we deduce: $p_{k}(\cdot)\rightarrow p(\cdot)$ uniformly in $\overline{\Omega}$ as $k\to\infty$ , where $p(x)=1+g\left(\left|\left(\nabla G_\sigma\ast v\right)(x)\right|\right)$ in $\Omega$ . The proof is complete.

3 Problem statement

The problem of restoration of cloud contaminated optical images is to reconstruct the original multi-band image $\vec{u}_0$ in the damage region D using the exact information about this image in $\Omega\setminus D$ (the undamaged region) and the texture (geometry) of the corresponding SAR image on the subset D. We begin with the following key assumptions.

Assumption 1. The cloud contaminated optical image $\vec{u}_0\in L^2(\Omega\setminus D;\mathbb{R}^M)$ and the corresponding SAR image $u_{SAR}\in BV(\Omega)$ are rigidly co-registered.

Assumption 2. The intensities $u_i$ of all spectral channels for the retrieved image $\vec{u}$ are subjected to the constraints $\gamma_{0,i}\le u_i(x)\le \gamma_{1,i}$ a.a. in $\Omega$ , where

(3.1) \begin{equation} \gamma_{0,i}=\inf_{x\in\Omega\setminus D} u_{i,0},\qquad \gamma_{1,i}=\sup_{x\in\Omega\setminus D} u_{i,0}. \end{equation}

Assumption 3. The topographic maps for all spectral channels $u_i$ , $i\in\left\{1,\dots,M\right\}$ of the retrieved image $\vec{u}=\left[u_1, u_2,\dots,u_M\right]^t:\Omega\rightarrow \mathbb{R}^M$ have a similar geometrical structure to the topographic map of the SAR image $u_{SAR}\in BV(\Omega)$ albeit these images can have very different intensities.

As follows from Assumption 3, all spectral channels of the retrieved image should share the geometry of the panchromatic SAR image $u_{SAR}$ in D. It means that, due to the property $u_{SAR}\in BV(\Omega)$ , for almost all points of almost all level sets of $u_{SAR}$ we can define a vector field $\theta\in L^\infty(\Omega,\mathbb{R}^2)$ such that $\theta(x)$ has the direction of the normal to the level lines of $u_{SAR}$ (see Remark 2.1 for the details). Therefore, the counterclockwise rotation of angle $\pi/2$ , denoted by $\theta^\perp$ , represents the tangent vector to the level lines of $u_{SAR}$ . In this case, if the spectral channels of $\vec{u}$ share the geometry of the panchromatic image $u_{SAR}$ , we have

(3.2) \begin{equation}\left(\theta^\perp,\nabla u_i\right)=0\quad \text{a.e. in {D}}, \ \forall\, i\in \left\{1,\dots,M\right\}.\end{equation}

We say that a function $\vec{u}=\left[u_1, u_2,\dots,u_M\right]^t:\Omega\rightarrow \mathbb{R}^M$ is the result of restoration of a cloud-corrupted image $\vec{u}_0: D\rightarrow\mathbb{R}^M$ if for given regularisation parameters $\mu{>}0$ , $\alpha{>}1$ , and $\lambda{>}0$ , each spectral component $u_i$ is the solution of the following constrained minimisation problem with the nonstandard growth energy functional

(3.3) \begin{align}\notag\left(\mathcal{P}_i\right)\quad J_i(v,p)&:=\int_{\Omega}\frac{1}{p(x)}|R_\eta\nabla v(x)|^{p(x)}\,dx+\frac{\mu}{\alpha}\int_{\Omega\setminus D}\left|{v(x)}-{u}_{0,i}(x)\right|^{\alpha}\,dx\\[3pt] &\quad+\lambda\int_{D}\Big|\left(\theta^{\perp},\nabla v\right)\Big|^\alpha\,dx\longrightarrow \inf_{(v, p)\in \Xi_i},\end{align}

where

• • $\Xi_i$ stands for the set of feasible solutions to the problem (3.3) which we define as follows

  \begin{equation*} \Xi_i=\left\{(v,p)\ \left| \begin{array}{l} v\in W^{1,p(\cdot)}(\Omega),\ p\in C(\overline{\Omega}),\\ \\[-7pt] 1\le \gamma_{0,i}\le v(x)\le \gamma_{1,i}\ \text{a.a. in}\ \Omega,\\ \\[-7pt] p(x)=\mathfrak{F}(v(x))\ \text{ in $\Omega$.} \end{array} \right. \right\} \end{equation*}
  Here, $W^{1,p(\cdot)}(\Omega)$ is the Sobolev-Orlicz space (for the details we refer to Appendix B),
  \begin{equation*} \mathfrak{F}(v(x))=1+g\left(\left|\left(\nabla G_\sigma\ast v\right) (x)\right|\right), \end{equation*}
  and $g{:}[0,\infty)\rightarrow (0,\infty)$ is the edge-stopping function which we take it in the form $g(t)=\frac{1}{1+(t/a)^2}$ ;

• • $\theta\in L^\infty(D,\mathbb{R}^2)$ is a vector field such that

  \begin{equation*} |\theta(x)|\le 1\ \text{ and }\ \left(\theta(x),\nabla u_{SAR}(x)\right)_{\mathbb{R}^2}=| \nabla u_{SAR}(x)| \quad\text{a.e. in {D};} \end{equation*}

• • $R_\eta\nabla v:=\nabla v-\eta^2\chi_D \left(\theta,\nabla v\right) \theta$ stands for the so-called directional gradient (see [8, 9]), and $\eta\in (0,1)$ is a given threshold;

• • $\left(G_{\sigma}\ast v\right) (x)$ denotes the convolution of function v with the two-dimensional Gaussian filter kernel $G_\sigma$ .

Let us give a short motivation to the choice of the energy functional in the form (3.3). Since

(3.4) \begin{equation}R_\eta\nabla v=\nabla v\quad\text{in $\Omega\setminus D$ and }\ (1-\eta^2)|\nabla v|\le |R_\eta\nabla v|\le |\nabla v|\quad\text{in}\ {{D}}\end{equation}

with a given $\eta\in (0,1)$ , it follows that this term can be considered as the regularisation in the Sobolev-Orlicz space $W^{1,p(\cdot)}(\Omega)$ . On the other side, this term plays the role of a spacial data fidelity. Indeed, the main goal, we are going to follow in the restoration problem, is to preserve the following property: the geometry of each spectral channels in the damage zone of the retrieved image should be as close as possible to the geometry of the SAR image. Formally, it means that relations (3.2) have to be satisfied. Hence, the magnitude $\int_{D}\Big|\left(\theta^{\perp},\nabla v\right)\Big|^\alpha\,dx$ must be small enough for each spectral channel, where $\theta$ stands for the vector field of unit normals to the topographic map of the SAR image $u_{SAR}$ . In fact, we impose it in the energy functional (3.3) in the form of the last term. However, in order to enforce this property, we observe that the expression $R_\eta\nabla v$ can be reduced to $(1-\eta^2)\nabla v$ in those places of the damage region D where v is collinear to the unit normal $\theta$ and to $\nabla v$ if $\nabla v$ is orthogonal to $\theta$ . Thus, gradients of the spectral intensities that are aligned/collinear $\theta$ are favoured as long as $|\theta|>0$ . Moreover, this property is enforced by the exponent p(x). In fact, the third term basically has a similar effect as the first one, so, in some practical implementations it can be omitted.

Since $p(x)\approx 1$ in places in $\Omega$ where edges or discontinuities are present in the spectral channel v, and $p(x)\approx 2$ in places where v(x) is smooth or contains homogeneous features, the main benefit of the model (3.3) is the manner in which it accommodates the local image information. For the places where the spectral gradient is sufficiently large (i.e. likely edges), we deal with the so-called directional TV-based diffusion [8, 9], whereas the places where the gradient 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. This enables the model to have a much lower dependence on the approximation schemes for the variable exponent p(x) and other thresholds. Apparently, the idea to involve a fractional norm of $W^{1,p(\cdot)}(\Omega)$ with a variable exponent p(x) was firstly proposed in [7] in order to reduce the staircasing effect in the TV image restoration problem.

As for the second term in (3.3), it is the so-called data fidelity and it forces the minimiser v in domain $\Omega\setminus D$ to stay as close as possible to the spectral intensity of the cloud corrupted image $\vec{u}_0: D\rightarrow\mathbb{R}^M$ .

Thus, the proposed model (3.3) provides a completely new approach to restoration of non-smooth multi-band optical images $\vec{u}_0$ with the gap in damage region. The main characteristic feature of this model is that we involve into consideration the energy functional with nonstandard growth where the edge information for restoration of $\vec{u}_0\in L^2(\Omega\setminus D;\mathbb{R}^M)$ in D is accumulated in the variable exponent p(x) and in the directional gradients $R_\eta\nabla$ which we derive from the SAR data. However, in contrast to [4, 11], we do not apply any selection procedure for approximation of the variable exponent p(x) in (3.3).

We are now in a position to define what we mean by the solution of the restoration problem. Taking into account Assumptions (1)–(3) that we imposed above and the structure of the energy functionals $\mathcal{J}_i:\Xi_i\rightarrow \mathbb{R}$ , we say that a multi-band image $\vec{u}^{\,0}=\left[u^0_1, u^0_2,\dots,u^0_M\right]^t:\Omega\rightarrow \mathbb{R}^M$ is the result of restoration of the cloud contaminated optical image $\vec{u}_0\in L^2(\Omega\setminus D;\mathbb{R}^M)$ if, for each index $i\in\left\{1,\dots,M\right\}$ , the pair $\left(u_i^0, p_i^0\right)$ , where $p_i^0=1+g\left(\left|\left(\nabla G_\sigma\ast u_i^0\right)\right|\right)$ in $\Omega$ , is a solution of the constrained minimisation problem (3.3), that is,

\begin{equation*}\left(u_i^0, p_i^0\right)\in\Xi_i\ \text{ and }\ \mathcal{J}_i\left(u_i^0, p_i^0\right)=\inf_{\left(v,p\right)\in\Xi_i} \mathcal{J}_i\left(v,p\right).\end{equation*}

4 Existence results

Our main intention in this section is to show that constrained minimisation problem (3.3) is consistent and admits at least one solution. Because of the specific form of the energy functional $J_i(v,p)$ , the minimisation problem (3.3) is rather challenging (we can refer to [6, 7, 8, 11] for some specific details).

Following in many aspects the recent studies [15, 24], we give the following existence result.

Theorem 4.1 For each $i=1,\dots,M$ and given $\mu{>}0$ , $\lambda{>}0$ , $\eta\in (0,1)$ , $\theta\in L^\infty(D,\mathbb{R}^2)$ , and $\vec{u}_0\in L^2(\Omega\setminus D;\mathbb{R}^M)$ , the minimisation problem (3.3) admits at least one solution $(u^{rec}_i,p^{rec}_i)\in \Xi_i$ .

Proof. Since $\Xi_i\ne\emptyset$ and $0\le J_i(v,p)<+\infty$ for all $(v,p)\in \Xi_i$ , it follows that there exists a non-negative value $\zeta\ge 0$ such that $\zeta=\inf\limits_{(v,p)\in \Xi_i}J_i(v,p)$ . Let $\left\{(v_k,p_k)\right\}_{k\in \mathbb{N}}\subset \Xi_i$ be a minimising sequence to the problem (3.3), that is,

\begin{equation*} (v_k, p_k)\in \Xi_i,\ p_k(x)=1+g\left(\left|\left(\nabla G_\sigma\ast v_k\right) (x)\right|\right)\ \text{ in $\Omega$}\ \forall\,k\in \mathbb{N},\ \text{ and } \lim\limits_{k\to\infty} J_i\left(v_k,p_k\right)=\zeta. \end{equation*}

So, without lost of generality, we can suppose that $J_i\left(v_k,p_k\right)\le \zeta+1$ for all $k\in \mathbb{N}$ . From this and the initial assumptions, we deduce

(4.1) \begin{gather} \notag \int_{\Omega} |v_k(x)|^{\alpha}\,dx\le \int_{\Omega} \gamma_{1,i}^{\alpha}\,dx\le \gamma_{1,i}^{\alpha}|\Omega|,\quad\forall\,k\in\mathbb{N},\\[3pt] \int_{\Omega}|R_\eta\nabla v_k(x)|^{p_k(x)}\,dx\le 2\int_{\Omega}\frac{1}{p_k(x)}|R_\eta\nabla v_k(x)|^{p_k(x)}\,dx<2(\zeta+1), \quad\forall\,k\in\mathbb{N}, \end{gather}

where $\sup_{k\in \mathbb{N}}\left[ \sup_{x\in\Omega} p_k(x)\right]\le 2$ .

Utilising the fact that $v_k(x)\le \gamma_{1,i}$ for almost all $x\in\Omega$ , we infer the following estimate

\begin{align*} \|v_k\|_{L^1(\Omega)}&\le \gamma_{1,i}|\Omega|,\quad\forall\,k\in\mathbb{N}. \end{align*}

Then arguing as in Lemma 2.1, it can be shown that $p_k\in C^{0,1}(\overline{\Omega})$ and

(4.2) \begin{align} \alpha:=1+\delta\le p_k(x) \le \beta:=2,\quad \forall\,x\in\Omega,\quad \forall\,k\in\mathbb{N}, \end{align}

(4.3) \begin{align} \text{with }\quad \delta=\frac{a^2}{a^2+\|G_{\sigma}\|^2_{C^1(\overline{\Omega-\Omega})} \gamma_1^2|\Omega|^2}. \end{align}

Taking this fact into account, we deduce from (4.1), (4.2) and (B5) that

\begin{align*} \notag \|v_k\|_{W^{1,\alpha}(\Omega)}&=\left( \int_\Omega \left[|v_k(x)|^{\alpha}+ |\nabla v_k(x)|^{\alpha}\right]\, dx\right)^{1/\alpha}\\[1pt] \notag &\le \left(1+|\Omega|\right)^{1/\alpha}\left(\int_{\Omega}\left[|v_k(x)|^{p_k(x)}+|\nabla v_k(x)|^{p_k(x)}\right]\,dx+2\right)^{1/\alpha}\\[1pt] &\stackrel{\text{by (3.4)}}{\le}\left(\frac{1+|\Omega|}{(1-\eta^2)^2}\right)^{1/\alpha} \left(\int_{\Omega}\left[|v_k(x)|^{p_k(x)}+|R_\eta\nabla v_k(x)|^{p_k(x)}\right]\,dx+2\right)^{1/\alpha}\\[1pt] &\stackrel{\text{by (4.1)}}{\le} \left(\frac{1+|\Omega|}{(1-\eta^2)^2}\right)^{1/\alpha}\left(\gamma_1^{2}|\Omega|+2\zeta+4\right)^{1/\alpha} \end{align*}

uniformly with respect to $k\in\mathbb{N}$ . Therefore, there exists a subsequence of $\left\{v_{k}\right\}_{k\in \mathbb{N}}$ , still denoted by the same index, and a function $u^{rec}_i\in W^{1,\alpha}(\Omega)$ such that

(4.4) \begin{gather} \notag v_{k}\rightarrow u^{rec}_i\ \text{strongly in $L^q(\Omega)$ for all $q\in [1,\alpha^\ast)$},\\[3pt] v_{k}\rightharpoonup u^{rec}_i\ \text{weakly in $W^{1,\alpha}(\Omega)$ as $k\to\infty$}, \end{gather}

where, by Sobolev embedding theorem, $\alpha^\ast=\frac{2\alpha}{2-\alpha}=\frac{2+2\delta}{1-\delta}>2$ .

Moreover, passing to a subsequence if necessary, we have (see Proposition B.2 and Lemma 2.1):

(4.5) \begin{gather} \nonumber\\[-30pt] v_{k}(x)\rightarrow u^{rec}_i(x)\ \text{ a.e. in $\Omega$. } \\[3pt] \notag v_{k}\rightharpoonup u^{rec}_i\ \text{ weakly in $L^{p_k(\cdot)}(\Omega)$},\\[3pt] \notag \nabla v_{k}\rightharpoonup \nabla u^{rec}_i\ \text{ weakly in $L^{p_k(\cdot)}(\Omega;\mathbb{R}^2)$},\\[3pt] \notag p_k(\cdot)\rightarrow p^{rec}_i(\cdot)=1+g\left(\left|\left(\nabla G_\sigma\ast u^{rec}_i\right)(\cdot)\right|\right)\quad\text{uniformly in $\overline{\Omega}$ as $k\to\infty$}, \end{gather}

where $u^{rec}_i\in W^{1,p^{rec}(\cdot)}(\Omega)$ with $p^{rec}(x)=1+g\left(\left|\left(\nabla G_\sigma\ast u^{rec}_i\right) (x)\right|\right)$ in $\Omega$ .

Since $\gamma_{0,i}\le v_k(x)\le \gamma_{1,i}$ a.a. in $\Omega$ for all $k\in\mathbb{N}$ , it follows from (4.5) that the limit function $u^{rec}_i$ is also subjected to the same restriction. Thus, $\left(u^{rec}_i,p^{rec}_i\right)$ is a feasible solution to minimisation problem (3.3).

Let us show that $(u^{rec}_i, p^{rec}_i)$ is a minimiser of this problem. With that in mind, we note that in view of the obvious inequality

\begin{equation*} \left|{v_k(x)}-u_{0,i}(x)\right|^{\alpha}\le 2^{\alpha-1}\left(\left|{v_k(x)}\right|^{\alpha}+\left|u_{0,i}(x)\right|^{\alpha}\right) \end{equation*}

and the fact that $u_{0,i}\in L^2(\Omega\setminus D)$ , we have: the sequence $\left\{v_k(x)-u_{0,i}(x)\right\}_{k\in \mathbb{N}}$ is bounded in $L^{\alpha}(\Omega\setminus D)$ and converges weakly in $L^{\alpha}(\Omega\setminus D)$ to $u^{rec}_i-u_{0,i}$ . Hence, by Proposition B.2 (see (B17)), $u^{rec}_i-u_{0,i}\in L^{p^{rec}(\cdot)}(\Omega\setminus D)$ and

(4.6) \begin{equation} \liminf_{k\to\infty} \int_{\Omega\setminus D} |v_k(x)-u_{0,i}(x)|^{\alpha}\,dx\ge \int_{\Omega\setminus D} |u^{rec}_i(x)-u_{0,i}(x)|^{\alpha}\,dx. \end{equation}

As for the last terms in (3.3), in view of the weak convergence $v_{k}\rightarrow u^{rec}_i$ in $W^{1,\alpha}(\Omega)$ (see (4.4)), we have

(4.7) \begin{align} \liminf_{k\to\infty} \int_{D}\Big|\left(\theta^{\perp},\nabla v_k\right)\Big|^\alpha\,dx\ge\int_{D}\Big|\left(\theta^{\perp},\nabla u^{rec}_i\right)\Big|^\alpha\,dx\,dx. \end{align}

It remains to notice that due to the properties (4.1), (4.4) and the fact that $\theta\in L^\infty(D,\mathbb{R}^2)$ , the sequence $\left\{|R_\eta\nabla v_k|\in L^{p_k(\cdot)}(\Omega)\right\}_{k\in \mathbb{N}}$ is bounded and weakly convergence to $|R_\eta\nabla u^{rec}_i|$ in $L^\alpha(\Omega)$ . Hence, by Proposition B.2, the following lower semicontinuous property

(4.8) \begin{equation} \liminf_{k\to\infty} \int_{\Omega} \frac{1}{p_k(x)}|R_\eta\nabla v_k(x)|^{p_k(x)}\,dx\ge \int_{\Omega} \frac{1}{p(x)}|R_\eta\nabla u^{rec}_i(x)|^{p(x)}\,dx \end{equation}

holds true.

As a result, utilising relations (4.6), (4.7) and (4.8), we finally obtain

\begin{align*} \zeta=\inf\limits_{(v,p)\in \Xi_i}J_i(v,p)=\lim_{k\to\infty} J_i\left(v_k,p_k\right)=\liminf_{k\to\infty}J_i(v_{k}, p_k)\ge J_i(u^{rec}_i, p^{rec}_i). \end{align*}

Thus, $(u^{rec}_i,p^{rec}_i)$ is a minimiser to the problem (3.3), whereas its uniqueness remains as an open question.

5 Iterative algorithm based on the variational convergence of extremal problems

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

\begin{gather*}1\le \gamma_{0,i}\le v(x)\le \gamma_{1,i}\ \text{a.a. in}\ \Omega,\quad p(x)=1+g\left(\left|\left(\nabla G_\sigma\ast v\right) (x)\right|\right)\end{gather*}

that we impose on the set of feasible solutions $\Xi_i$ . This motivates us to pass to some relaxation (we refer to [16, 28, 33] where the similar approach has been utilised). 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 [25, 26, 27, 30].

At the first step, we set up

(5.1) \begin{align}p_0(x)=\left\{\begin{array}{l}1+g\left(\left|\left(\nabla G_\sigma\ast u_{0,i}\right) (x)\right|\right), \ \text{ if }\ x{\in} \Omega\setminus D,\\ \\[-7pt] 1+g\left(\left|\left(\nabla G_\sigma\ast u_{SAR}\right) (x)\right|\right), \ \text{ if }\ x{\in} D,\end{array}\right\}\end{align}

(5.2) \begin{align}u^{0}=\mathop{{Argmin}}_{v\in \mathcal{B}_{i,p_0(\cdot)}} J_i(v,p_0(\cdot)).\end{align}

Here, $\mathcal{B}_{i,p(\cdot)}=\big\{v\in W^{1,{p(\cdot)}}(\Omega): 1{\le} \gamma_{0,i}\le v(x)\le \gamma_{1,i}$ a.a. in $ \Omega\big\}$ .

Then, for each $k\ge 1$ , we set

(5.3) \begin{gather}p_k(x)=1+g\left(\left|\left(\nabla G_\sigma\ast u^{k-1}\right)(x)\right|\right),\quad\forall\,x\in\Omega,\quad u^{k}=\mathop{Argmin}_{v\in \mathcal{B}_{i,p_k(\cdot)}} J_i(v,p_k(\cdot)).\end{gather}

Before proceeding further, we set

(5.4) \begin{gather}\mathfrak{S}=\left\{h\in C(\overline{\Omega})\ \left|\begin{array}{c}|h(x)-h(y)|\le C|x-y|,\ \forall\,x,y\in \Omega,\\ \\[-8pt] \alpha:=1+\delta\le h(x) \le \beta:=2,\quad \forall\,x\in\Omega,\end{array}\right. \right\}\end{gather}

where $C>0$ and $\delta>0$ are defined by (2.13) and (4.3), respectively.

Arguing as in the proof of Theorem 4.1 and using the convexity arguments, it can be shown that, for each $p(\cdot)\in \mathfrak{S}$ , there exists a unique element $u^{0,p(\cdot)}_i\in \mathcal{B}_{i,p(\cdot)}$ such that $u^{0,p(\cdot)}_i=\mathop{Argmin}_{v\in \mathcal{B}_{i,p(\cdot)}} J_i(v,p(\cdot))$ . Moreover, it can be shown that, for given $i=1,\dots,M$ , $\mu>0$ , $\lambda{>}0$ , $\eta\in (0,1)$ , $u_{SAR}\in BV(\Omega)$ , and $\vec{u}_0\in L^2(\Omega\setminus D,\mathbb{R}^M)$ , the sequence $\left\{ u^{k}\in W^{1,{p_k(\cdot)}}(\Omega)\right\}_{k\in\mathbb{N}}$ is compact with respect to the weak topology of $W^{1,\alpha}(\Omega)$ , whereas the exponents $\left\{ p_k\right\}_{k\in\mathbb{N}}$ are compact with respect to the strong topology of $C(\overline{\Omega})$ .

We say that a pair $(\widehat{u}_i,\widehat{p}_i)$ is a weak solution to the original problem (3.3) if

(5.5) \begin{equation}\begin{array}{c}\widehat{u}_i=\mathop{Argmin}\limits_{v\in \mathcal{B}_{\widehat{p}_i(\cdot)}} J_i(v,\widehat{p}_i(\cdot)),\ \widehat{u}_i\in \mathcal{B}_{i,\widehat{p}_i(\cdot)},\\ \\[-7pt] \widehat{p}_i(x)=1+g\left(\left|\left(\nabla G_\sigma\ast \widehat{u}_i\right)(x)\right|\right),\ \forall\,x\in\Omega.\end{array}\end{equation}

Remark 5.1 The relation between a weak solution and a solution to the problem (3.3) is very intricate. Since the uniqueness of solutions to (3.3) is a rather questionable option, it follows that, in principle, these definitions can lead to the different pairs in $\Xi_i$ . As immediately follows from (5.5), a weak solution is a merely feasible one to the original problem. However, if the problem (3.3) admits a unique solution $(u_i^0,p_i^0)\in \Xi_i$ , then (5.5) implies that this pair is a weak solution. On the other hand, if $(u_i^\ast,p_i^\ast)\in \Xi_i$ is some solution to (3.3), then setting in (5.1) $p_0(x)=p_i^\ast(x)$ for all $x\in\Omega$ , we immediately arrive at the conclusion: $(u^k,p_k)=(u_i^\ast,p_i^\ast)$ for all $k\in\mathbb{N}$ and, therefore, $(u_i^\ast,p_i^\ast)$ is a weak solution to the above problem.

Our main goal in this section is to establish the existence of a weak solution to the original problem (3.3) and show that it can be attained by some iterative algorithm.

We begin with some technical results.

Lemma 5.1 For given $\mu{>}0$ , $\lambda{>}0$ , $\eta\in (0,1)$ , $u_{SAR}{\in} BV(\Omega)$ and $\vec{u}_0{\in} L^2(\Omega{\setminus} D,\mathbb{R}^M)$ , the sequence $\left\{ u^{k}\in W^{1,{p_k(\cdot)}}(\Omega)\right\}_{k\in\mathbb{N}}$ is compact with respect to the weak topology of $W^{1,\alpha}(\Omega)$ .

Proof. To begin with, let us show that the sequence of minimisers $\left\{ u^{k} \right\}_{k\in\mathbb{N}}$ is bounded in the sense of condition (B10). Let $\widehat{u}\in C^1(\overline{\Omega})$ be an arbitrary function such that $\gamma_{0,i}\le \widehat{u}(x)\le \gamma_{1,i}$ in $\Omega$ . Since

(5.6) \begin{equation} J_i(u^k,p_k(\cdot)) \le J_i(\widehat{u},p_k(\cdot)),\quad\forall\,k=0,1,2,\dots \end{equation}

and

\begin{align*} \int_{\Omega}\frac{1}{p_k(x)}&|R_\eta\nabla \widehat{u}(x)|^{p_k(x)}\,dx \le \int_{\Omega}\frac{1}{p_k(x)}|\nabla \widehat{u}(x)|^{p_k(x)}\,dx\\[3pt] &\le \int_{\Omega}\frac{1}{p_k(x)}\left(1+\|\widehat{u}\|_{C^1(\overline{\Omega})}\right)^{p_k(x)}\,dx\le \frac{|\Omega|}{\alpha}\left(1+\|\widehat{u}\|_{C^1(\overline{\Omega})}\right)^2,\\[3pt] \int_{\Omega\setminus D}\frac{1}{p_k(x)}&\left|{\widehat{u}(x)}-u_{0,i}(x)\right|^{\alpha}\,dx\le 2\int_{\Omega\setminus D}\left[\left(1+\|\widehat{u}\|_{C^1(\overline{\Omega})}\right)^{\alpha}+ \left|u_{0,i}(x)\right|^{\alpha}\right]\,dx\\[3pt] &\le 2|\Omega\setminus D|\left(1+\|\widehat{u}\|_{C^1(\overline{\Omega})}\right)^\alpha+2|\Omega|^\frac{2-\alpha}{2}\|u_{0,i}\|^\alpha_{L^{2}(\Omega)}, \end{align*}

it follows that

\begin{equation*} \sup_{k\in \mathbb{N}}J_i(u^k,p_k(\cdot))\le \sup_{k\in \mathbb{N}}J_i(\widehat{u},p_k(\cdot))\le \widehat{C} \end{equation*}

with some constant $\widehat{C}>0$ .

From this and estimate (3.4), we deduce

(5.7) \begin{align} \int_{\Omega} |u^{k}(x)|^{\alpha}\,dx\le \gamma_{1,i}^{\alpha}|\Omega|,\quad\forall\,k\in\mathbb{N}, \end{align}

(5.8) \begin{align} \int_{\Omega}|\nabla u^{k}(x)|^{p_k(x)}\,dx\le 2\int_{\Omega}\frac{1}{p_k(x)}|\nabla u^{k}(x)|^{p_k(x)}\,dx<\frac{2}{(1-\eta^2)^2}\widehat{C}, \quad\forall\,k\in\mathbb{N}. \end{align}

Then estimate (B5) implies that the sequence $\left\{ u^{k} \right\}_{k\in\mathbb{N}}$ is bounded in $W^{1,\alpha}(\Omega)$ . So, its weak compactness is a direct consequence of the reflexivity of $W^{1,\alpha}(\Omega)$ .

We notice that boundedness of the sequence $\left\{ u^{k} \right\}_{k\in\mathbb{N}}$ in $W^{1,{\alpha}}(\Omega)$ and compactness of the embedding $W^{1,\alpha}(\Omega)\hookrightarrow L^q(\Omega)$ for $q\in \left[1,\frac{2\alpha}{2-\alpha}\right)$ imply the existence of element $u^\ast\in W^{1,\alpha}(\Omega)$ such that, up to a subsequence,

(5.9) \begin{align} u^{k}(x)\rightarrow u^\ast(x)\ \text{a.e. in }\ \Omega,\end{align}

(5.10) \begin{align}u^{k}\rightharpoonup u^\ast\ \text{ in }\ L^\alpha(\Omega),\ \text{ and }\ \nabla u^{k}\rightharpoonup \nabla u^\ast\ \text{ in }\ L^\alpha(\Omega;\mathbb{R}^2).\end{align}

Then using (5.9) and passing to the limit in two-side inequality $\gamma_{0,i}\le u^k(x)\le \gamma_{1,i}$ , we obtain

(5.11) \begin{equation}\gamma_{0,i}\le u^\ast(x)\le \gamma_{1,i}\quad\text{for a.a. }\ x\in\Omega.\end{equation}

Utilising this fact together with the pointwise convergence (5.9), by the Lebesgue dominated convergence theorem, we have

(5.12) \begin{align}\lim_{k\to\infty} \mathfrak{F}(u^{k}(x))&=1+\frac{a^2}{a^2+(\lim\limits_{k\to\infty}\left|\left(\nabla G_\sigma\ast u^{k}\right) (x)\right|)^2} \nonumber\\ & =1+\frac{a^2}{a^2+\Big(\left|\left(\nabla G_\sigma\ast \lim\limits_{k\to\infty}u^{k}\right) (x)\right|\Big)^2}=\mathfrak{F}(u^\ast(x)),\quad\forall\,x\in\Omega.\end{align}

Since, by Arzelà–Ascoli theorem, the set $\left\{p_k=1+g\left(\left|\left(\nabla G_\sigma\ast u^{k-1}\right)(x)\right|\right)\right\}_{k\in\mathbb{N}}$ is compact with respect to the strong topology of $C(\overline{\Omega})$ , it follows from (5.12) (see also the proof of Lemma 2.1) that

(5.13) \begin{equation}p_k\rightarrow p^\ast=\mathfrak{F}(u^\ast(x))\ \text{ strongly in $C(\overline{\Omega})$ as $k\to\infty$}, \ \text{ and $p^\ast\in\mathfrak{S}$}.\end{equation}

Then properties (5.9)–(5.13) and Proposition B.2 imply:

(5.14) \begin{equation}u^\ast\in \mathcal{B}_{i,p^\ast(\cdot)}=\left\{u\in W^{1,{p^\ast(\cdot)}}(\Omega)\ :\ 1\le \gamma_{0,i}\le u(x)\le \gamma_{1,i}\ \text{a.a. in}\ \Omega\right\}.\end{equation}

Thus, $(u^\ast,p^\ast)\in \mathcal{B}_{i,p^\ast(\cdot)}\times \mathfrak{S}$ is a cluster point of iterative procedure (5.3) with respect to the convergence (5.9)–(5.10), (5.13).

The main result of this section can be stated as follows:

Theorem 5.1 Let $\mu{>}0$ , $\lambda{>}0$ , $\eta\in (0,1)$ , $u_{SAR}{\in} BV(\Omega)$ and $\vec{u}_0{\in} L^2(\Omega{\setminus} D,\mathbb{R}^M)$ be given data. Then, for each $i\in\left\{1,\dots,M\right\}$ , the sequence of approximated solutions $\left\{(u^k,p_k)\right\}_{k\in\mathbb{N}}$ possesses the following asymptotic properties:

(5.15) \begin{align} u^{k}(x)\rightarrow \widetilde{u}_i(x)\ \text{a.e. in }\ \Omega, \end{align}

(5.16) \begin{align} u^{k}\rightharpoonup \widetilde{u}_i\ \text{ in }\ L^\alpha(\Omega),\ \text{ and }\ \nabla u^{k}\rightharpoonup \nabla \widetilde{u}_i\ \text{ in }\ L^\alpha(\Omega;\mathbb{R}^2),\end{align}

(5.17) \begin{align} p_k\rightarrow \widetilde{p}_i=\mathfrak{F}(\widetilde{u}_i)\ \text{ strongly in $C(\overline{\Omega})$ as $k\to\infty$}, \end{align}

where $(\widetilde{u}_i,\widetilde{p}_i)$ is a weak solution to the original problem (3.3), that is,

\begin{gather*} \widetilde{u}_i\in \mathcal{B}_{i,\widetilde{p}_i(\cdot)},\quad \widetilde{u}_i=\mathop{Argmin}_{v\in \mathcal{B}_{i,\widetilde{p}_i(\cdot)}} J_i(v,\widetilde{p}_i(\cdot)), \end{gather*}

and, in addition, the following variational property holds true

(5.18) \begin{align} \notag \lim_{k\to\infty}J_i(u^{k},p_k(\cdot))&= \lim_{k\to\infty}\left[\inf_{v\in \mathcal{B}_{i,p_k(\cdot)}}J_i(v,p_k(\cdot))\right]\\[4pt] &=\inf_{v\in \mathcal{B}_{i,\widetilde{p}_i(\cdot)}}J_i(v,\widetilde{p}_i(\cdot))=J_i(\widetilde{u}_i,\widetilde{p}_i(\cdot)). \end{align}

Proof. Let $(\widetilde{u}_i,\widetilde{p}_i)$ be a cluster point of the sequence $\left\{(u^k,p_k)\right\}_{k\in\mathbb{N}}$ with respect to the convergence (5.15)–(5.17). The existence of such pair immediately follows from Lemma 5.1 and reasoning given above. Let’s assume the converse – namely, there is a function $z\in \mathcal{B}_{i,\widetilde{p}_i(\cdot)}$ such that

(5.19) \begin{equation} J_i(z,\widetilde{p}_i(\cdot))=\inf_{v\in \mathcal{B}_{i,\widetilde{p}_i(\cdot)}} J_i(v,\widetilde{p}_i(\cdot))<J_i(\widetilde{u}_i,\widetilde{p}_i(\cdot)). \end{equation}

Using the procedure of the direct smoothing, we set

\begin{gather*} u_{\varepsilon}(x)=\frac{1}{{\varepsilon}^2}\int_{\mathbb{R}^2} K\left(\frac{x-s}{{\varepsilon}}\right)\widetilde{z}(s)\,ds, \end{gather*}

where ${\varepsilon} >0$ is a small parameter, K is a positive compactly supported smooth function with properties

\begin{equation*} K\in C_0^\infty(\mathbb{R}^2),\ \int_{\mathbb{R}^2}K(x)\,dx=1,\ \text{ and }\ K(x)=K(-x), \end{equation*}

and $\ \widetilde{z}\ $ is zero extension of z outside of $\Omega$ .

Since $z\in W^{1,\widetilde{p}(\cdot)}(\Omega)$ and $\widetilde{p}(x)\ge \alpha=1+\delta$ in $\Omega$ , it follows that $z\in W^{1,\alpha}(\Omega)$ . Then,

(5.20) \begin{gather} \notag \text{$u_{\varepsilon}\in C^\infty_0(\mathbb{R}^2)$ for each ${\varepsilon}>0$},\\[3pt] u_{\varepsilon}\rightarrow z\ \text{ in }\ L^\alpha(\Omega),\quad \nabla u_{\varepsilon}\rightarrow \nabla z\ \text{ in }\ L^\alpha(\Omega;\mathbb{R}^2) \end{gather}

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

(5.21) \begin{equation} u_{\varepsilon}(x)\rightarrow z(x)\ \text{ a.e. in $\Omega$.} \end{equation}

Moreover, taking into account the estimates

\begin{align*} u_{\varepsilon}(x)&=\int_{\mathbb{R}^2} K\left(y\right)\widetilde{z}(x-{\varepsilon} y)\,dy\le \gamma_{1,i}\int_{\mathbb{R}^2} K\left(y\right)\,dy=\gamma_{1,i},\\[3pt] u_{\varepsilon}(x)&\ge\int_{y\in {\varepsilon}^{-1}(x-\Omega)} K\left(y\right)\widetilde{z}(x-{\varepsilon} y)\,dy\ge \gamma_{0,i}\int_{y\in {\varepsilon}^{-1}(x-\Omega)} K\left(y\right)\,dy\ge\gamma_{0,i}, \end{align*}

we see that each element $u_{\varepsilon}$ is subjected to the pointwise constraints

\begin{gather*} \gamma_{0,i}\le u_{\varepsilon}(x)\le \gamma_{1,i}\ \text{a.a. in}\ \Omega,\quad \forall\,{\varepsilon}>0. \end{gather*}

Since, for each ${\varepsilon}>0$ , $u_{\varepsilon}\in W^{1,p_k(\cdot)}(\Omega)$ for all $k\in\mathbb{N}$ , it follows that $u_{\varepsilon}\in \mathcal{B}_{i,p_k(\cdot)}$ , that is, each element of the sequence $\left\{u_{\varepsilon}\right\}_{{\varepsilon}>0}$ is a feasible solution to all approximating problems $\big\langle\inf_{v\in \mathcal{B}_{i,p_k(\cdot)}} J_i(v,p_k(\cdot)) \big\rangle$ . Hence,

(5.22) \begin{align} J_i(u^{k},p_k(\cdot))\le J_i(u_{\varepsilon},p_k(\cdot)),\quad\forall\,{\varepsilon}>0,\ \forall\,k=0,1,\dots \end{align}

Further we notice that

(5.23) \begin{equation} \liminf_{k\to\infty} J_i(u^{k},p_k(\cdot))\ge J_i(\widetilde{u}_i,\widetilde{p}_i(\cdot)) \end{equation}

by Proposition B.2 and Fatou’s lemma, and

(5.24) \begin{align} \notag \lim_{k\to\infty} J_i(u_{\varepsilon},p_k(\cdot)) &= \lim_{k\to\infty} \int_{\Omega}\frac{1}{p_k(x)}|R_\eta\nabla u_{\varepsilon}(x)|^{p_k(x)}\,dx\\[3pt] &\quad +\frac{\mu}{\alpha}\int_{\Omega\setminus D}\left|{u_{\varepsilon}(x)}-{u}_{0,i}(x)\right|^{\alpha}\,dx +\lambda\int_{D}\Big|\left(\theta^{\perp},\nabla u_{\varepsilon}\right)\Big|^{\alpha}\,dx. \end{align}

Since

\begin{gather*} \frac{1}{p_k(x)}|R_\eta\nabla u_{\varepsilon}(x)|^{p_k(x)}\rightarrow \frac{1}{\widetilde{p}_i(x)}|R_\eta\nabla u_{\varepsilon}(x)|^{\widetilde{p}_i(x)}\quad\text{uniformly in $\Omega$ as $k\to\infty$}, \end{gather*}

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

(5.25) \begin{equation} \lim_{k\to\infty} J_i(u_{\varepsilon},p_k(\cdot))= J_i(u_{\varepsilon},\widetilde{p}_i(\cdot)),\quad\forall{\varepsilon}>0. \end{equation}

As a result, passing to the limit in (5.22) and utilising properties (5.23)–(5.25), we obtain

(5.26) \begin{align} \notag J_i(\widetilde{u}_i,\widetilde{p}_i(\cdot))&\le J_i(u_{\varepsilon},\widetilde{p}_i(\cdot)) = \int_{\Omega}\frac{1}{\widetilde{p}_i(x)}|R_\eta\nabla u_{\varepsilon}(x)|^{\widetilde{p}_i(x)}\,dx\\[3pt] &\quad +\frac{\mu}{\alpha}\int_{\Omega\setminus D}\left|{u_{\varepsilon}(x)}-{u}_{0,i}(x)\right|^{\alpha}\,dx +\lambda_2\int_{D}\Big|\left(\theta^{\perp},\nabla u_{\varepsilon}\right)\Big|^{\alpha}\,dx, \end{align}

for all ${\varepsilon}>0$ . Taking into account the pointwise convergence (see (5.21) and property (5.20))

\begin{gather*} |R_\eta\nabla u_{\varepsilon}(x)|^{\widetilde{p}_i(x)}\rightarrow |R_\eta\nabla z(x)|^{\widetilde{p}_i(x)},\\[3pt] \left|{u_{\varepsilon}(x)}-{u}_{0,i}(x)\right|^{\alpha}\rightarrow \left|{z(x)}-{u}_{0,i}(x)\right|^{\alpha},\\[3pt] \Big|\left(\theta^{\perp},\nabla u_{\varepsilon}\right)\Big|^{\alpha}\rightarrow \Big|\left(\theta^{\perp},\nabla z\right)\Big|^{\alpha} \end{gather*}

as ${\varepsilon}\to 0$ , and the fact that, for ${\varepsilon}$ small enough,

\begin{align*} |R_\eta\nabla u_{\varepsilon}(x)|^{\widetilde{p}_i(x)}&\le \left(1+\left|\nabla {z(\cdot)}\right|\right)^{\widetilde{p}_i(\cdot)}\in L^1(\Omega)\quad\text{a.e. in $\Omega$},\\[3pt] \left|{u_{\varepsilon}(\cdot)}-{u}_{0,i}(\cdot)\right|^{\alpha}&\le \left[2\left(1+\left|z(\cdot)\right|\right)^{\alpha}+2\left(1+\left|{u}_{0,i}(\cdot)\right|\right)^{2}\right]\in L^1(\Omega)\quad\text{a.e. in $\Omega\setminus D$},\\[3pt] \Big|\left(\theta^{\perp},\nabla u_{\varepsilon}\right)\Big|^{\alpha}&\le \|\theta\|_{L^\infty(D,\mathbb{R}^2)} \left(1+\left|\nabla z(\cdot)\right|\right)^{\widetilde{p}_i(\cdot)}\in L^1(\Omega)\quad\text{a.e. in $\Omega$}, \end{align*}

we can pass to the limit in (5.26) as ${\varepsilon}\to 0$ by the Lebesgue dominated convergence theorem. This yields

\begin{equation*} J_i(\widetilde{u}_i,\widetilde{p}_i(\cdot))\le\lim_{{\varepsilon}\to 0} J_i(u_{\varepsilon},\widetilde{p}_i(\cdot))=J_i(z,\widetilde{p}_i(\cdot)). \end{equation*}

Combining this inequality with (5.26) and (5.19), we finally get

\begin{equation*} J_i(z,\widetilde{p}_i(\cdot))=\inf_{v\in \mathcal{B}_{i,\widetilde{p}_i(\cdot)}} J_i(v,\widetilde{p}_i(\cdot)) <J_i(\widetilde{u}_i,\widetilde{p}_i(\cdot))\le J_i(z,\widetilde{p}_i(\cdot)), \end{equation*}

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

(5.27) \begin{equation} J_i(\widetilde{u}_i,\widetilde{p}_i(\cdot))=\inf_{v\in \mathcal{B}_{i,\widetilde{p}_i(\cdot)}} J_i(v,\widetilde{p}_i(\cdot)) \end{equation}

and, therefore, $(\widetilde{u}_i,\widetilde{p}_i)$ is a weak solution to the original problem (3.3). As for the variational property (5.18), it is a direct consequence of (5.27) and (5.25).

6 Optimality conditions

To characterise the solution $u^{0,p(\cdot)}\in \mathcal{B}_{i,p(\cdot)}$ of the approximating optimisation problem $\langle\inf_{v\in \mathcal{B}_{i,p(\cdot)}}J_i(v,p(\cdot))\rangle$ , we check that the functional $J_i(v,p(\cdot))$ is Gâteaux differentiable with respect to v, that is,

(6.1) \begin{align}&\lim_{t\to 0} \frac{J_i(u^{0,p(\cdot)}+tv,p(\cdot))-J_i(u^{0,p(\cdot)},p(\cdot))}{t} \nonumber\\[3pt]&\quad =\int_{\Omega}\left(|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)-2}R_\eta\nabla u^{0,p(\cdot)}(x),R_\eta\nabla v(x)\right)\,dx \nonumber\\[3pt]&\qquad +\mu\int_{\Omega\setminus D}\left|u^{0,p(\cdot)}(x)-{u}_{0,i}(x)\right|^{\alpha-2}u^{0,p(\cdot)}(x)v(x)\,dx \nonumber\\[3pt]&\qquad +\alpha\lambda\int_{\Omega}\Big|\left(\theta^{\perp},\nabla u^{0,p(\cdot)}\right)\Big|^{\alpha-1} \left(\theta^{\perp},\nabla v\right)\,dx,\quad\forall\,v\in W^{1,{p(\cdot)}}(\Omega).\end{align}

To this end, we note that

\begin{align*}&\frac{|R_\eta\nabla u^{0,p(\cdot)}(x)+tR_\eta\nabla v(x)|^{p(x)}-|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)}}{p(x) t}\\[3pt] &\qquad\qquad\qquad\qquad \rightarrow\left(|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)-2}R_\eta\nabla u^{0,p(\cdot)}(x),R_\eta\nabla v(x)\right)\ \text{ as }\ t\to 0\end{align*}

almost everywhere in $\Omega$ . Since, by convexity,

\begin{equation*}|\xi|^p-|\eta|^p\le 2p\left(|\xi|^{p-1}+|\eta|^{p-1}\right)|\xi-\eta|,\end{equation*}

it follows that

(6.2) \begin{align}&\left|\frac{|R_\eta\nabla u^{0,p(\cdot)}(x)+tR_\eta\nabla v(x)|^{p(x)}-|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)}}{p(x) t}\right|\nonumber\\[3pt]&\quad \le 2\left(|R_\eta\nabla u^{0,p(\cdot)}(x)+tR_\eta\nabla v(x)|^{p(x)-1}+|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)-1}\right)|R_\eta\nabla v(x)| \nonumber\\[3pt] &\quad \le {const}\, \left(|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)-1}+|R_\eta\nabla v(x)|^{p(x)-1}\right)|R_\eta\nabla v(x)|.\end{align}

Taking into account that

\begin{align*} &\int_\Omega |R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)-1}|R_\eta\nabla v(x)|\,dx\\[3pt] & \quad \stackrel{\text{by (B.18)}}{\le}2\|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)-1}\|_{L^{p^\prime(\cdot)}(\Omega)}\|R_\eta\nabla v(x)|\|_{L^{p(\cdot)}(\Omega)}\\[3pt] & \quad \le 2\|\nabla u^{0,p(\cdot)}(x)|^{p(x)-1}\|_{L^{p^\prime(\cdot)}(\Omega,\mathbb{R}^2)}\|\nabla v(x)\|_{L^{p(\cdot)}(\Omega,\mathbb{R}^2)},\end{align*}

and $\int_\Omega |\nabla v(x)|^{p(x)}\,dx\stackrel{\text{by (B5)}}{\le} \|\nabla v\|^2_{L^{p(\cdot)}(\Omega,\mathbb{R}^2)}+1,$ we see that the right hand side of inequality (6.2) is an $L^1(\Omega)$ function. Therefore,

\begin{align*}&\int_\Omega\frac{|R_\eta\nabla u^{0,p(\cdot)}(x)+tR_\eta\nabla v(x)|^{p(x)}-|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)}}{p(x) t}\,dx\\[3pt] &\quad \rightarrow\int_\Omega\left(|R_\eta\nabla u^{0,p(\cdot)}(x)|^{p(x)-2}R_\eta\nabla u^{0,p(\cdot)}(x),R_\eta\nabla v(x)\right)\,dx\ \text{ as }\ t\to 0\end{align*}

by the Lebesgue dominated convergence theorem.

Utilising similar arguments to the rest terms in (3.3), we deduce that the representation (6.1) for the Gâteaux differential of $J_i(\cdot,p(\cdot))$ at the point $u^{0,p(\cdot)}\in \mathcal{B}_{i,p(\cdot)}$ is valid.

Thus, in order to derive some optimality conditions for the minimising element $u^{0,p(\cdot)}\in \mathcal{B}_{i,p(\cdot)}$ to the problem $\inf\limits_{v\in \mathcal{B}_{p(\cdot)}} J_i(v,p(\cdot))$ , we note that $\mathcal{B}_{i,p(\cdot)}$ is a non-empty convex subset of $W^{1,{p(\cdot)}}(\Omega)$ and the objective functional $J_i(\cdot,p(\cdot)): \mathcal{B}_{i,p(\cdot)}\rightarrow\mathbb{R}$ is strictly convex. Hence, the well-known classical result (see [32, Theorem 1.1.3]) and representation (6.1) lead us to the following conclusion.

Theorem 6.1. Let $p_k(\cdot)\in \mathfrak{S}$ be an exponent given by the iterative rule (5.3). Then, the unique minimiser $u^{k}\in \mathcal{B}_{i,p_k(\cdot)}$ to the approximating problem $\inf\limits_{v\in \mathcal{B}_{i,p_k(\cdot)}}J_i(v,p_k(\cdot))$ is characterised by

(6.3) \begin{align} & \int_{\Omega}\left(\Big|R_\eta\nabla u^{k}(x)\Big|^{p_k(x)-2}R_\eta\nabla u^{k}(x),R_\eta\nabla v(x)-R_\eta\nabla u^{k}(x)\right)\,dx \nonumber\\[3pt] &\quad\qquad\qquad\quad +\mu\int_{\Omega\setminus D}\Big|u^{k}(x)-{u}_{0,i}(x)\Big|^{\alpha-2}u^{k}(x)\left(v(x)-u^{k}(x)\right)\,dx \nonumber \\[3pt] &\quad\qquad +\alpha\lambda\int_{\Omega}\Big|\left(\theta^{\perp},\nabla u^{0,p(\cdot)}\right)\Big|^{\alpha-1} \left(\theta^{\perp},\nabla v-\nabla u^{k}\right)\,dx\ge 0,\quad\forall\,v\in \mathcal{B}_{i,p_k(\cdot)}. \end{align}

7 Numerical experiments

In order to illustrate the proposed algorithm for the restoration of the cloud-corrupted satellite multispectral optical images, we have provided some numerical experiments. As input data, we have used a series of optical images and a radar image that were delivered from two twin satellites, Sentinel-2A and Sentinel-2B. Each of the optical image is corrupted by clouds with a different level of degradation (see Figures 1, 2, 3, 4). The corresponding SAR image is shown in Figure 5.

Figure 1. Optical image from 2021/040.

Figure 2. Optical image from 2021/032.

Figure 3. Optical image from 2021/041.

Figure 4. Optical image from 2021/0406.

Figure 5. The SAR image of the same territory.

The results of restoration after 5 iterations of the proposed algorithm are shown in Figures 6, 7, 8, 9 for $\eta=0.8$ , $\mu=10$ , and $\lambda=20$ . As for the exponent $\alpha$ , we define it following the rule (2.10) with $\delta={a^2}\left[a^2+\|G_{\sigma}\|^2_{C^1(\overline{\Omega-\Omega})} 255^2|\Omega|^2\right]^{-1}$ and $a=0.01$ .

Figure 6. Result of its restoration.

Figure 7. Result of its restoration.

Figure 8. Result of its restoration.

Figure 9. Result of its restoration.

Comparing the restored images and the contaminated ones, we could see that the texture of original images is well preserved. However, some details in the damage zones are blurred. This problem will be addressed in the following research.

8 Conclusion

This paper proposes a novel restoration model for the satellite optical images that is based on their co-registration with SAR images of the same region and the solution of a special variational problem with nonstandard growth objective functional. The characteristic feature of the restoration problem is the fact that the optical images are supposed to be corrupted and there is a subdomain (the so-called damage region) where we do not have any information about the brightness of such images. Instead, we propose to make use of the structure of the SAR images of the same regions. The novelty of our approach is that we involve into consideration the objective functional with nonstandard growth p(x), where the variable exponent is unknown a priori and it directly depends on the texture of an image that has to be restored. In order to study the consistency of the proposed non-convex minimisation problem, we develop a special technique and supply this approach by the rigorous mathematical substantiation.

However, this problem is not trivial from its practical implementation and has limited performance from numerical point of view. With that in mind, we propose an alternating minimisation procedure which is based on the concept of relaxation of extremal problems with the state constraints and their variational convergence. We show that at each step of the iteration procedure, we obtain a strictly convex optimisation problem which admits a unique solution. It is established that the so-defined sequence of approximations is precompact and each convergent subsequence leads to the so-called weak solutions of the original problem. So, following this way some weak solutions of the restoration problem can be attained in the suitable topology by a special iterative algorithm.

As for the numerical scheme for practical implementation of the proposed procedure, we use the straightforward finite difference method (see [12] for some details) and it will be described in the forthcoming paper. Besides of that, in the second part of this paper we mainly focus on the solvability issues of variational inequality (6.3), its qualitative analysis and the study of its asymptotic behaviour as k tends to $\infty$ .