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 [Reference Hnatushenko, Kogut and Uvarow21] 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 [Reference Ballester, Caselles, Igual, Verdera and Rougé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 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}

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 [Reference Bungert, Coomes, Ehrhardt, Rasch, Reisenhofer and Schönlieb8, Reference Bungert and Ehrhardt9]), 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 [Reference Bungert, Coomes, Ehrhardt, Rasch, Reisenhofer and Schönlieb8, Reference Bungert and Ehrhardt9], 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 [Reference Blomgren, Chan, Mulet and Wong7] 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 [Reference Antil and Rautenberg4, Reference Chen, Levine and Rao11], 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 [Reference Blomgren6, Reference Blomgren, Chan, Mulet and Wong7, Reference Bungert, Coomes, Ehrhardt, Rasch, Reisenhofer and Schönlieb8, Reference Chen, Levine and Rao11] for some specific details).

Following in many aspects the recent studies [Reference D’Apice, De Maio and Kogut15, Reference Kogut24], we give the following existence result.

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 [Reference De Maio, Kogut and Manzo16, Reference Kogut, Kupenko and Manzo28, Reference Manzo33] 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 [Reference Kogut25, Reference Kogut26, Reference Kogut27, Reference Kogut and Leugering30].

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}

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.

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 [Reference Evans and Gariepy18]). 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 [Reference Lions32, 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.