VARIATIONAL APPROACH FOR THE RECONSTRUCTION OF DAMAGED OPTICAL SATELLITE IMAGES THROUGH THEIR CO-REGISTRATION WITH SYNTHETIC APERTURE RADAR

In this paper the problem of reconstruction of damaged multi-band optical images is studied in the case where we have no information about brightness of such images in the damage region. Mostly motivated by the crop field monitoring problem, we propose a new variational approach for exact reconstruction of damaged multi-band images using results of their co-registration with Synthetic Aperture Radar (SAR) images of the same regions. We discuss the consistency of the proposed problem, give the scheme for its regularization, derive the corresponding optimality system, and describe in detail the algorithm for the practical implementation of the reconstruction procedure.


Introduction
It is well-known that visible red, green, and blue bands and also near-infrared (NIR) and SWIR regions of the electromagnetic spectrum of optical satellite images have been used successfully to monitor crop cover, crop health, soil moisture, nitrogen stress, and crop yield (see, for instance, [14,37,43]). In view of this, qualitative analysis of vegetation and detection of changes in vegetation patterns are the keys to natural resource assessment and monitoring. Thus, it comes as no surprise that the detection and quantitative assessment of crop cover and green vegetation biomass is one of the major applications of remote sensing for environmental resource management and decision making.
However, in spite of the fact that optical satellite multi-band images have a high resolution and can be easily captured by low-cost cameras, they are often corrupted because of the poor weather conditions, such as rain, clouds, fog, and dust conditions. Moreover, it is a typical situation when the measure of degradation of optical images is such that we can not rely even on the brightness recovery of the damaged regions. As a result, some subdomains of such images become invisible or their coverage by the spectral vegetation indices is rather far from to be reliable and consistent. However, in contrast to the optical observation, the 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 SAR and optical images is very important for classification of land cover [42] and estimation of soil moisture to remove vegetation cover effects from radar backscattering coefficient [19,20,47]. At the same time, because of the distinct natures of SAR and optical images, there exist huge radiometric and geometric differences between optical and synthetic aperture of radar images. As a result, their structure and texture are drastically different. Because of this, it would be naive to suppose that the gray level of the original color image u 0 or its brightness on the damaged region D can be successfully recovered at high level of accuracy through SAR images of this region. Thus, in spite of the fact that in the literature there are many approaches to the reconstruction of an image when information of the colors is not everywhere available (see, for instance, [26,35,48,49,51]), the traditional approaches to the exact reconstruction of damaged color images are no longer applicable in this case and it makes this problem challenging.
The aim of this paper is propose and study the variational model for exact reconstruction of damaged multi-band optical satellite images using results of their co-registration with SAR images of the same regions. The variational approach we consider is inspired, in some sense, by the famous ROF model for denoising, introduced by Rudin, Osher and Fatemi in the context of grey level functions (see [45]): to minimize BV (Ω) u → |Du|(Ω) + λ u − u 0 2 L 2 (Ω), (1.1) where Ω ⊂ R 2 denotes the image domain, u 0 ∈ L 2 (Ω) is the given image, and λ > 0 is a tunning parameter. In order to be able to reconstruct edges in the image, it is represented by a function of bounded variation u ∈ BV (Ω). So, the first term in (1.1) is the total variation |Du|(Ω) which has a regularizing effect but at the same time allows for discontinuities which may represent adges in the image. The second term is the fidelity term which measures the distance to the given image. Often, some weaker norms such as the H −1 -norm are considered to define the latter term. However, the non-differentiability of the total variation is challenging from a computational point of view. Moreover, in some papers (see, for instance, [3,23]), it was shown that natural images are incompletely represented by BV (Ω) ∩ L 2 (Ω) functions. In fact, it can be shown that BV (Ω) ∩ L ∞ (Ω) is contained in the fractional Sobolev space H s (Ω) for 0 < s < 1/2 [52]. As a result, they propose to replace the total variation term in (1.1) by a squared fractional Sobolev norm. In other words, they involve in minimization the following energy with 0 < s < 1 and β ∈ [0, 1], where (−∆) s denotes the fractional power of the Laplacian with zero Neumann boundary conditions. Then the first necessary and sufficient optimality condition determines the unique minimizer u via 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 [44]. Since, from the reconstruction point of view, it is desirable that the regularity of the solution to (1.2) is low in places in Ω where edges or discontinuities are present in u true, and that is high in places where u true is smooth or contains homogeneous features, it is of interest to consider (1.2) where s : Ω → [0, 1] is not a constant. For the details we refer to the recent paper [4].
Later on it was shown that the ROF model can be extended to various image processing problems, one of which is T V inpainting method [11]. Since the colorization task can also be understood as inpainting the colors (as Sapiro's insight [48]), the T V minimizing approach has been widely used for different problems related with colorization of black and white images in computer graphics, and also with reconstruction of damaged color images [18].
However, in contrast to the standard setting of reconstruction problem for color damaged images where the starting point is either the knowledge of the gray level of the original color image u 0 on a given open subset D of Ω (the damaged region) together with the exact information of u 0 on Ω \ D (the undamaged region) or the grey level information in the damage region D ⊂ Ω is modeled as a nonlinear distortion of the colors, in this paper we deal with the case where we do not have any information about u 0 inside D but instead we assume that a SAR image u SAR : Ω → R of the same region is given.
When dealing with multi-band optical satellite images, a color of each pixel can be identified with a vector ξ = (ξ 1 , . . . , ξ M ) t ∈ R M , where components ξ i corresponds to the different channels which are crucial for calculation of the major vegetation indices that have a wide application in many agricultural monitoring services. For instance, These indexes are well established proxies for the crop conditions and give us early insights into how well the crops are doing and if they are in need of water or nutrients [56]. So, when we speak about multi-band optical satellite images, we suppose that they have at least red, green, blue, NIR, and SWIR channels. Let M be the number of channels that we are allowed to use.
There are two preferred ways to represent multi-band images mathematically. The first one is the model, where an image u ∈ BV (Ω, R M ) is represented via its M channels u = (u 1 , . . . , u M ). The other way to represent an image is called Chromaticity/Brightness (or CB-model), where a multi-band image u ∈ BV (Ω, R M ) is decomposed into two components: its chromaticity C := u/|u| = (u 1 /|u|, . . . , u M /|u|) and its brightness B := |u| = u 2 1 + · · · + u 2 M . It is wellknown that treating brightness separately from the chromaticity can give more flexibility in detail recovery (see, for instance [17,27]). In particular, in [10], the authors showed that CB model gives better color control and detail recovery for color image denoising compared to different color settings. In this paper, we use CB multi-band model for the reconstruction of optical satellite images.
The paper is organized as follows. After recalling basic notions and background in Section 2, we give in Section 3 the precise setting of the reconstruction problem for multi-band damaged optical images. We do it in several steps. First we discuss the procedure of denoising and selective smoothing for the SAR data of the same region assuming that these data are well co-registered with the original optical image. Moreover, the main focus of this procedure is to preserve the edges inside the damage region D after transformations related to denosing and selective smoothing of SAR images. At the second step, we set the reconstruction problem of the brightness B 0 = |u 0 | ∈ BV (Ω \ D) of a multi-band optical image in the damage domain D ⊂ Ω. The novelty of model that we propose, is that the edge information for the brightness reconstruction is derived from the SAR data. In some sense, this model is related to piecewise smooth Mumford-Shah segmentation [39], since our model yields smooth spectral bands for each homogeneous region, while the edge information is enforced by the special weight function. We also present in this section a variational model for recovery of the chromaticity data C in the damage domain D ⊂ Ω provided the chromaticity components C 0 for the original multi-band image u 0 are well defined in Ω \ D.
With that in mind we introduce a special constrained minimization problem where the minimization of the cost functional with respect to the chromaticity C is proposed to affect by the brightness data B 0 in Ω \ D and by the smoothed SAR data u * SAR in the damage region D. As a result, we expect that the diffusion of chromaticity is inhibited across the edges of B 0 in Ω \ D and the edges of u * SAR in D, yielding a sharp transition in the function C in the rest regions.
After the chromaticity C rec and the brightness B rec are recovered by solving of the corresponding minimization problems, the fully multi-band image u in Ω can be restored by the rule u = C rec B rec . Thus, we give the way for the reconstruction of the major vegetation indices in the damage region D. We notice that the proposed reconstruction scheme is rather flexible. Each spectral band diffuses in D as long as it meets the boundary of region enforced by the SAR information. Moreover, by using the chromaticity model, natural spectral band blending is possible with respect to geodesic direction in chromaticity space S M −1 . In a homogeneous region, if different bands are given, this model will naturally diffuse the spectral band by diffusing the vector values in S M −1 in geodesic direction between bands.
Our main intention in Section 4 is to show that constrained minimization problems for the chromaticity and the brightness recovery, are consistent and the corresponding sets of their solutions are nonempty. Since the chromaticity recovery problem is not trivial in its practical implementation (because of the nonconvex state constraint C(x) ∈ S M −1 ), we discuss its relaxation and asymptotic properties of the sequence of minimizers for penalized minimization problems in Section 5.
Optimality conditions for the penalized chromaticity recovery problem and the brightness reconstruction problem are studied in Sections 6-7. The last section is devoted to the short description of the crucial steps of alternative minimization method that we suggest for the numerical simulations of the spectral indices by the damaged multi-band satellite optical images.

Notation and Preliminaries
We begin with some notation. For vectors ξ ∈ R N and η ∈ R N , (ξ, η) = ξ t η denotes the standard vector inner product in R N , where t denotes the transpose operator. The norm |ξ| is the Euclidean norm given by |ξ| = (ξ, ξ).
Let Ω ⊂ R 2 be a bounded open set with a Lipschitz boundary ∂Ω. For any subset D ⊂ Ω we denote by |D| its 2-dimensional Lebesgue measure L 2 (D). For a subset D ⊆ Ω let D denote its closure and ∂D its boundary. We define the characteristic function χ D of D by For a function u we denote by u| D its restriction to the set D ⊆ Ω, and by u ∂D its trace on ∂D. Let C ∞ 0 (Ω) be the infinitely differentiable functions with compact support in Ω. For a Banach space X its dual is X * and ·, · X * ;X is the duality form on X * × X. By and * we denote the weak and weak * convergence in normed spaces, respectively. For given M ∈ N and 1 p +∞, the space L p (Ω; R M ) is defined by The inner product of two functions f and g in L p (Ω; Let D (Ω) be the dual of the space C ∞ 0 (Ω), i.e. D (Ω) is the space of distributions in Ω. By H 1 0 (Ω) we denote the closure of C ∞ 0 (Ω)-functions with respect to the norm Then H 1 0 (Ω) is a Banach space and the norm in H 1 0 (Ω) can be defined by We denote the dual of H 1 0 (Ω) by H −1 (Ω). Then (see [46, p.401]), H −1 (Ω) is isometrically isomorphic to the Hilbert space of all distributions F ∈ D (Ω) satisfying It can be shown that the standard norm in H −1 (Ω) is equivalent to the following one (for the details we refer to [24]) BV (Ω; R M ), endowed with the norm , and the supremum is taken on the set of functions Remark 2.1. For a Borel set B ⊂ Ω and an arbitrary function u ∈ BV (Ω; R M ), the mapping B → |Du|(B) is a Radon measure that is lower semi-continuous with respect to the L 1 -convergence of sets.
We recall that a sequence ) and, therefore, the notatioń Ω φ dDf k should be interpreted as followŝ converges strongly to some f in L 1 (Ω; R M ) and sup k∈N´Ω d|Df k | < +∞, then (see, for instance, [1] and [5]) So, a simple criterion for weak * convergence can be states as follows (see [1, p.125]): and converges to u strongly in L 1 (Ω; R M ).
The following embedding results for BV -function plays a crucial role for variational problems that we study in this paper (see [5, p.378]).
We also recall the Poincare-Wirtinger inequality: in two dimensional case, there exists a constant C P W such that, for any u ∈ BV (Ω; R M ), we have By analogy with the theory of Sobolev spaces, the notion of trace operator can be extended for BV -functions. Namely, the following result is well-known [ • the Green's formula holds: ∀ ϕ ∈ C 1 (Ω; R 2 ), where ν(x) is the outer unit normal at H 1 -almost all x in ∂Ω.
A special case of functions of bounded variation are those that are characteristic functions of sets of finite perimeter and were introduced by R. Caccioppoli in [8].
Let χ E be its characteristic function. We say that E is a set with finite perimeter in Ω if χ E ∈ BV (Ω). This means that the distributional gradient Dχ E is a vectorvalued measure with finite total variation. The total variation |Dχ E |(Ω) is called the perimeter of E in Ω, i.e., P (E, Ω) = |Dχ E |(Ω) and, therefore, It is known that for every set E with finite perimeter in Ω the following generalized Gauss-Green formula holdŝ where ν E is the inner unit normal to E.
Since the sets with finite perimeter are not smooth in general, the correct way to represent the measure Dχ E is to introduce the so-called reduced boundary ∂ * E.
In this way, for every set E of finite perimeter we have The following properties are well-known: (c) E → P (E, Ω) is lower semicontinuous with respect to the convergence in measure in Ω; (d) For any two sets E 1 and E 2 with finite perimeters in Ω, the relation holds, with equality holding if only the distance between these sets in Euclidean space R 2 is non-zero.
For further information concerning functions of bounded variation, we refer to [1,5].
We recall also some auxiliary results concerning the vector fields z ∈ L ∞ (Ω; R 2 ) whose divergence in the sense of distribution are L 2 (Ω)-functions. We denote by L ∞ 2,div (Ω; R 2 ) the space of all such vector fields z. Then the trace of the normal component of the field z ∈ L ∞ 2,div (Ω; R 2 ) on ∂Ω can be defined as distribution Tr(z, ∂Ω) in the sense of Anzellotti (see [2]). In particular, having assumed that the original domain Ω ⊂ R 2 is of class C 1 , the trace of the normal component of z on ∂Ω is the distribution defined as follows For example, if z is a piecewise continuous vector field that can be extended continuously in Ω, then (see [12, p. 22]) Utilizing the property ∂Ω = ∂ * Ω, we have the following result (the so-called Gauss-Green formula) (see [12,Theorem 5.1] and [28, Proposition 6.12] for the details).
Lemma 2.1. For any u ∈ BV (Ω) and z ∈ L ∞ 2,div (Ω; R 2 ) there exists a Radone measure on Ω denoted by (z, Du) and a function Tr u ∂Ω ∈ L 1 (∂Ω, dH 1 ) stands for the trace of u ∈ BV (Ω) on ∂Ω. The measure (z, Du) and | (z, Du) | are absolutely continuous with respect to |Du| and, for any open Ω ⊂ Ω, ∀ ϕ ∈ C ∞ 0 ( Ω), and for all Borel sets Ω ⊂ Ω, we have Moreover, it turns out that in this case the following estimate holds true

Setting of the Problem
Let Ω ⊂ R 2 be a bounded image domain with Lipschitz boundary ∂Ω, and let D ⊂ Ω be a Borel set with non empty interior and sufficiently regular boundary and such that |Ω \ D| > 0. We call D the damage region of a given multi-band image u 0 ∈ BV (Ω \ D; R M ). As it was mentioned in Introduction, we deal with the case where we have no information about the original image u 0 inside D. So, the brightness B 0 = |u 0 | and the chromaticity components C 0 = u 0 /B 0 are only defined in Ω\D. Instead of this, we assume that a SAR image u SAR : Ω → R of the same region is given, and this image is well co-registered with u 0 ∈ BV (Ω\D; R M ) in Ω \ D. By default, we assume that all functions C 0 , B 0 , and u SAR take values in the set of strictly positive real numbers R + . The problem is to reconstruct the original multi-band image u in the damaged domain D ⊂ Ω.
We will do it in several steps. First, we realize the denoising and selective smoothing procedure for the SAR data u SAR in order to preserve the edges inside the damage region D. With that in mind, we propose to make use of Perona-Malik equation [41] in its combination with the image empirical mode decomposition method (IEMD) [38].

Denoising and selective smoothing procedure for the SAR data
Formally, this procedure can be described as follows.
A2. Perform the procedure of selective smoothing and denoising for the image J n : Ω → R. With that in mind, we define a smooth version V n of J n as a solution of the following initial-boundary value problem where ∂Vn ∂ν is the normal derivative of V n at the boundary of Ω, g : [0, ∞) → (0, ∞) is a continuous monotone decreasing function such that g(0) = 1 and g(t) > 0 for all t > 0 with lim t→+∞ g(t) = 0, and G σ * I determines the convolution of function I with the Gaussian kernel G σ : Here, σ determines the width of the Gaussian kernel and we will refer to σ as the inner scale.
A3. Find the IEMD for Z n := V n + R n , i.e., represent Z n in the form of its decomposition of its intrinsic oscillatory modes (known as IMFs) and the last residue, A4. Set up n := n + 1, R n := R L(n),n , and J n = L(n) k=1 C k,n (x, y); A5. If n N 0 then go to step A2. Otherwise the procedure should be stopped and image u * SAR := Z n declared as a denoised and selective smoothed version of the SAR image u SAR .
We notice that the time variable t in the evolution equation (3.1) corresponds to a spatial scale analogous to σ, and by the regularity result in [41], we have V n ∈ C 1 loc (0, ∞; H 1 (Ω)) for each n ∈ N. Typical examples of the edge function g are

Reconstruction of the brightness
At this stage, we deal with the reconstruction problem of the brightness B 0 = |u 0 | ∈ BV (Ω\D) in the damage domain D ⊂ Ω. In order to recover the brightness data everywhere in D, we propose to solve the following constrained minimization problem (see [50] for comparison) where χ Ω\D (x) = 1 for x ∈ Ω \ D Here, from physical point of view, we suppose that a feasible brightness should be non-negative. As for the regularizing term R(B), its role is to fill in the brightness content into the damage domain D, e.g., by diffusion and/or transport. After pioneering works of Masnou and Morel [40] and Bertalmio et. al [7], the typical choice of regularizing term is to use the total variation, i.e. R(B) = |DB|(Ω) =´Ω d|DB|. Other examples to be mentioned for R(B) are the active contour model based on Mumford and Shah's segmentation [54], the inpainting scheme based on the Mumford-Shah-Euler image model [16], and the Euler elastica model, where with positive weights a and b. As follows from the coarea formula, we have the following properties where Γ s = {x ∈ Ω : B(x) = s} is the level line of the brightness for the grayvalue s. This circumstance motivate us to specify the cost functional J(B) in (3.7) as follows where λ > 0 is a weight coefficient, g : [0, ∞) → (0, ∞) represents an edge function with properties described in A2 (see Subsection 3.1), u * SAR is the denoised and selective smoothed version of the SAR image u SAR , * denotes the convolution operator, and G σ stands for the Gaussian kernel (3.4).
As follows from (3.7), the minimization of the functional (3.9) with respect to the brightness B is now affected not only by the brightness data B 0 in Ω \ D and but also by the smoothed SAR data u * SAR in the damage region D. The first term in the functional J is the functional weighted by the linear combination of the edge functions g (see [27,57] for its utilization in the color image inpainting problems). In view of the properties of the edge function g, we see that the value of χ Ω\D g(|∇G σ * B 0 |)+χ D g(|∇u * SAR |) is close to one in regions of Ω\D where the original brightness B 0 is slowly varying, and in regions of D where the smoothed SAR image u * SAR is also slowly varying. At the same time, this value is small at the edges of brightness and the SAR image, if both σ and the constant a in (3.6) are small enough. Hence, the first term of J acts as a regularization functional such that the diffusion of brightness is inhibited across the edges of B 0 in Ω \ D and across the edges of u * SAR in D, yielding a sharp transition in the function B. At the same time, the second term in the cost functional J requires the BV -function B to be close to the brightness data B 0 in Ω \ D.
Thus, the novelty of model (3.7), (3.9) is that the edge information for the brightness reconstruction is derived from the SAR data. In some sense, this model is related to piecewise smooth Mumford-Shah segmentation [39], since our model yields smooth spectral bands for each homogeneous region, while edge information is enforced by the weight function χ Ω\D g(|∇G σ * B 0 |) + χ D g(|∇u * SAR |).

Recovery of chromaticity via SAR-weighted harmonic map
In this subsection we present a variational model for recovery of the chromaticity data C in the damage domain D ⊂ Ω provided the chromaticity components C 0 for the original multi-band image u 0 are only defined in Ω \ D.
We consider the following cost functional for reconstruction of the chromaticity data C in D: where α > 0 is a weight coefficient, and the rest counterparts are as in the cost functional (3.9). We associate with (3.10) the corresponding constrained minimization problem So, by analogy with the previous subsection, the minimization of the functional (3.10) with respect to the chromaticity C is also proposed to affect by the brightness data B 0 in Ω \ D and by the smoothed SAR data u * SAR in the damage region D. As a result, we expect that the diffusion of chromaticity is inhibited across the edges of B 0 in Ω \ D and the edges of u * SAR in D, yielding a sharp transition in the function C.
Thus, the novelty of model (3.11) is that the edge information for chromaticity recovery is derived from the brightness and SAR data, and what is more important, this problem can be solved separately from the brightness reconstruction one.

Final recovery of multi-band optical image
After the chromaticity C rec and the brightness B rec are recovered by solving minimization problems (P 1 ) and (P 2 ), respectively, the fully multi-band image u in Ω will be given by u = C rec B rec . As a result, we give the way for the reconstruction of the major vegetation indices in the damage region D. In order to weight up all pros and cons of this approach, we notice that the proposed reconstruction scheme is rather flexible. Each spectral band diffuses in D as long as it meets the boundary of region enforced by the SAR information. Moreover, by using the chromaticity model, natural spectral band blending is possible with respect to geodesic direction in chromaticity space S M −1 . In a homogeneous region, if different bands are given, this model will naturally diffuse the spectral band by diffusing the vector values in S M −1 in geodesic direction between bands.

Existence Results
Our main interest in this section is to show that constrained minimization problems (3.7) and (3.11) are consistent and the corresponding sets of their solutions are nonempty. We begin with the following result (for comparison, we refer to [9,36,50]). Proof. First we notice that minimization problem (3.7)-(3.9) is consistent, that is, J(B) < +∞ for any feasible B ∈ Ξ. Indeed, since ∂Ω is a Lipschitz domain in R 2 , it follows from Poincare inequality (2.4) that BV (Ω) space is continuously embedded into L 2 (Ω). On the other hand, by Sobolev embedding theorem we have a compact injection H 1 0 (Ω) → L 2 (Ω). Hence, by the duality arguments and Riesz representation theorem, we deduce that L 2 (Ω) = L 2 (Ω) * → H 1 0 (Ω) * = H −1 (Ω) with a compact embedding as well. Thus, to finalize the remark about consistency, it is enough to observe that B 0 ∈ L 2 (Ω) and g(|∇G σ * B 0 |) together with g(|∇G σ * u * SAR |) are continuous functions. Hence, and therefore, J(B) < +∞ provided B ∈ Ξ. As follows from (3.9), the infimum in (3.7) is finite. Let {B k } k∈N be a minimizing sequence for (3.7), i.e. lim Let us show that, in fact, the minimizing sequence {B k } k∈N is bounded in BV (Ω). Indeed, by smoothness of |∇G σ * B 0 | and |∇G σ * u * SAR |, and positiveness of B 0 and u * SAR , we deduce the existence of a constant γ > 0 such that Hence,ˆΩ d|DB k | by (4.14) Then, by Poincaré-Wirtinger inequality (2.4), there exists a constant C P W > 0 depending only on Ω such that Let us show that, in fact, we can indicate a constant C > 0 such that Arguing as in [36], we set w k = 1 |Ω|´Ω B k dx χ Ω and v k = B k − w k . Then it is clear that w k , v k ∈ BV (Ω) for all k ∈ N. Moreover, by compactness of the embedding Since it follows from (4.5) and the inequality with some constant C > 0 independent of k. This implies that Hence, Thus, due to compactness of the embedding L 2 (Ω) → H −1 (Ω), there exists a constant C * > 0 depending on Ω such that From this and (4.3), we finally deduce that uniformly with respect to k ∈ N. Therefore, there exists a subsequence of {B k } k∈N , still denoted by the same index, and a function B rec ∈ BV (Ω) such that Moreover, passing to a subsequence if necessary, we have: in Ω, and B k B rec weakly in L 2 (Ω). (4.10) in Ω for all k ∈ N, it follows from (4.10) 1 that the limit function B rec is also subjected the same restriction. Thus, B rec is a feasible solution of reconstruction problem (3.7)-(3.9). Taking into account the weak convergence in L 2 (Ω) (see (4.10) 2 ), by compactness of the embedding L 2 (Ω) → H −1 (Ω), it implies that B k → B rec strongly in H −1 (Ω).
(4.11) Therefore, utilizing properties (4.9) 2 and (4.11), we get As a result, it follows from the above consideration that Since the set of feasible solutions Ξ is convex and closed in BV (Ω), by Mazur's theorem we have that this set is sequentially closed with respect to the weak * convergence in BV (Ω). Thus, B rec ∈ Ξ and B rec is a minimizer for constrained minimization problem (3.7)-(3.9). It remains to show that B rec is a unique minimizer for this problem. Indeed, let us assume the converse. Let B ∈ Ξ and D ∈ Ξ be two minimizers for the problem (3.7)-(3.9). Then by the strict convexity of norm · H −1 (Ω) and convexity of the set of feasible solutions Ξ, we have which brings us into conflict with the initial assumptions. Thus, B rec is a unique minimizer to the problem (3.7)-(3.9). The proof is complete.
We proceed further with the study of the SAR-weighted chromaticity recovery problem (3.11). Our aim is to show that this problem has a solution. The proof relies on the following Poincaré type of inequality (see [13,Lemma 4.1]). (4.13) We are now in a position to prove the existence result concerning the constrained minimization problem (3.11) with the convex functional F : BV (Ω; R M ) → R + and the non-convex sphere constraint C(x) ∈ S M −1 for a.a. x ∈ Ω. Since the indicated type of constraints is not trivial to satisfy, in practice the problem (P 2 ) requires a relaxation in passing from the non-convex set |C| = 1 to the convex unit ball |C| 1 with the corresponding penalization of this constraint into the minimizing functional. This issue will be considered in details in the next section. Let Ω ⊂ R 2 be a bounded image domain with Lipschitz boundary ∂Ω, and let D ⊂ Ω be a damage region such that |Ω \ D| > 0 and int D = ∅. Then, for any given B 0 ∈ BV (Ω \ D) and u SAR ∈ L ∞ (Ω) such that B 0 > 0 and u * SAR > 0 a.e. in Ω, the minimization problem (P 2 ) has a solution, i.e. there exists a function C rec ∈ BV (Ω; S M −1 ) such that Proof. Since 0 F (C) < +∞ for all C ∈ BV (Ω; S M −1 ), it follows that there exists a non-negative value µ 0 such that µ = inf Let {C k } k∈N ⊂ BV (Ω; S M −1 ) be a minimizing sequence to the problem (3.11), i.e. C k ∈ BV (Ω; S M −1 ), ∀ k ∈ N, and lim k→∞ F (C k ) = µ. So, we can legitimately suppose that F (C k ) µ + 1 for all k ∈ N. Then By smoothness of |∇G σ * B 0 | and |∇G σ * u * SAR |, and positiveness of B 0 and u * SAR , we deduce the existence of a positive constant γ such that χ Ω\D g(|∇G σ * B 0 |) + χ D g(|∇G σ * u * SAR |) > γ in Ω. Hence, Ω d|DC k | by (4.14) < γ −1 (µ + 1).

(4.16)
Taking into account the estimatê it follows from Lemma 4.1 that So, there exists a function of chromaticity C rec ∈ BV (Ω; R M ) such that, up to a (not relabeled) subsequence, C k → C rec strongly in L 1 (Ω; R M ) and C k (x) → C rec (x) almost everywhere in Ω. Since C k (x) ∈ S M −1 for a.a. x ∈ Ω, it follows from the pointwise convergence property that C rec ∈ BV (Ω; S M −1 ). It remains to take into account estimate (4.16) and the lower semi-continuity property (2.3) in order to deduce |DC rec |(Ω) lim inf k→∞ |DC k |(Ω). (4.17) Utilizing again the pointwise convergence C k (x) → C rec (x) almost everywhere in Ω and Fatou's lemma, we get Hence, by (4.17) and (4.18), we finally obtain so that, C rec ∈ BV (Ω; S M −1 ) is a solution of the minimization problem (3.11).

On Relaxation of the Chromaticity Recovery Problem
As it was mentioned in the previous sections, the constrained minimization problem (P 2 ) is not trivial in its practical implementation because of the nonconvex state constraint C(x) ∈ S M −1 . It is worth to notice that in recent years, minimization problems over manifold-valued data has been attracted many interest (see, for instance, [6,55,57] and the other references therein). Many interesting and well-suited approaches for non-convex optimization problems were proposed including augmented Lagrangian methods, penalty methods, alternative direction minimization, and others. In this section, by analogy with the recent results developed in [6] (see also [29][30][31][32]34]), we make use of the relaxation approach passing from the non-convex constrained set |C(x)| = 1 to the convex unit ball |C(x)| 1 in R M with further penalization of this constraint in the corresponding minimization cost functional. With that in mind, for any real number ε > 0 and any given B 0 ∈ BV (Ω\D) and u SAR ∈ L ∞ (Ω), we consider the convex functional and the corresponding minimization problem with a convex constraint Hereinafter, we assume that the parameter ε varies within a strictly decreasing sequence of positive real numbers which converge to 0. When we write ε > 0, we consider only the elements of this sequence. We begin with the following existence result for the penalized variational problem (P ε ).
Proof. The proof of an existence result is similar to that of Theorem 4.1. For the uniqueness of such solution it is enough to apply the convexity arguments. Therefore, we sketch only the main points. Due to the standard properties of the convolution operator, we see that |∇G σ * B 0 | and |∇G σ * u * SAR | are bounded functions in the closure of Ω \ D and D, respectively. Hence, by the definition of the edge function g, it follows that there exists a positive constant γ such that As a result, we deduce: the cost functional F ε is coercive on the set of feasible solutions Λ = C ∈ BV (Ω; R M ) : |C(x)| 1 for a.a. x ∈ Ω and lower semi-continuous with respect to the weak * convergence in BV (Ω; R M ). Then a solution of the variational problem (P ε ) exists for any ε > 0.
Let C rec ε be a solution of (P ε ) for a given ε > 0. Let us show that the sequence {C rec ε } ε>0 is bounded in BV (Ω; R M ). To this end, we set C is a constant vector field in R M such that | C| = 1 at each point x ∈ R M . Then it is clear that C ∈ Λ. Therefore, Since the right-hand side of (5.3) does not depend on ε, it follows that where µ : Taking into account that sup where the constant C > 0 comes from the Poincaré inequality (4.13).
We proceed further with the study of asymptotic properties of the sequence of minimizers {C rec ε } ε>0 as ε → 0. For the technique and more details, we refer to [33].
Theorem 5.1. Let {C rec ε } ε>0 be a sequence of minimizers to the penalized minimization problems (P ε ). Then this sequence is compact with respect to the weak * convergence in BV (Ω; R M ), and each its weakly * cluster point is a solution to the chromaticity recovery problem (3.11).
Proof. Since the set Λ is sequentially closed with respect to the weak * convergence in BV (Ω; R M ), and the uniformly bounded sets in BV -norm are relatively compact in L 1 (Ω), it follows from (5.6) that there exists a function C 0 ∈ Λ ⊂ BV (Ω; R M ) and a subsequence C rec ε k k∈N such that C rec ε k → C 0 strongly in L 1 (Ω; R M ) and almost everywhere in Ω as k → ∞. From this we deduce that Since, for any ε > 0, we have Utilizing this fact together with (5.7), we obtain: |C rec ε k | → 1 strongly in L 2 (Ω), and |C 0 (x)| = 1 for a.a. x ∈ Ω. Hence, C 0 ∈ BV (Ω; S M −1 ).
It remains to show that the function C 0 solves the chromaticity recovery problem (P 2 ). Indeed, for an arbitrary function C ∈ BV (Ω; S M −1 ), we obviously have |C|(x) = 1 a.e. in Ω, and, therefore, Then passing in the right-hand side of (5.9) to the limit as k → ∞ and taking into account the Fatou's lemma and the lower semi-continuity property of the total variation (as we did in the proof of Theorem 4.1), we have Hence, and this implies that F C 0 = inf C∈BV (Ω;S M −1 ) F (C).

Optimality System for the Penalized Chromaticity Recovery Problem
The main object of our consideration in this section is the constrained minimization problem (5.2). Following in general aspects the technique of Temam for the problem of minimal surfaces [53] and duality results from [15], we derive in this section the necessary optimality conditions in order to characterize the solution C rec ∈ BV (Ω; S M −1 ) the problem (P ε ). To begin with, we reformulate problem (5.2) as an equivalent problem on the space X := L 2 (Ω; R M ). For this reason we define the following functional on X: where K (B 0 , u SAR ) is given by (4.12).
Notice that all functional are well-defined on X. In view of this, we extend the energy functional F ε to the entire set L 2 (Ω; R M ) by the rule Then it is clear that a minimizer C rec ∈ Ξ of (5.2) is also a minimizer of the modified problem because the inclusion C ∈ Λ is equivalent to the consistency condition of the problem (6.6), i.e.
C ∈ Λ ⇔ F ε (C) < +∞ and G(C) 0. (6.7) Our next step is to derive some necessary condition for minimizers of a constrained minimization problem (6.6). A fundamental difficulty that typically appears in this case is the lack in differentiability of the energy functional (6.5). Since the functionals F ε : L 2 (Ω; R M ) → R and G : L 2 (Ω; R M ) → R are convex, a necessary condition for a minimizer C rec ε ∈ Λ should employ the convex subdifferential ∂ F ε (C rec ε ) which is some subsets of the dual space L 2 (Ω; R M ) * = L 2 (Ω; R M ).
Before proceeding further, we recall the definition of the subdifferential ∂F (u) of a convex proper functional F : X → R ∪ +∞ at some element u ∈ X. If X = L 2 (Ω), then X * = X and, for given u ∈ X, an element ξ ∈ X belongs to ∂F (u) if and only if, ∀ v ∈ X, Thus, ξ ∈ ∂F (u) if u is a minimizer on X of the following variational problem We begin with the following technical results (for the proof and their substantiation we refer to [24, Section 3]). Proposition 6.1. Let B 0 ∈ BV (Ω\D), C 0 ∈ BV (Ω; S M −1 ), and u SAR ∈ L ∞ (Ω) be given functions. Then the functionals E 2 , E 2 : L 2 (Ω; R M ) → R are convex and Gâteaux differentiable in L 2 (Ω; R M ) with , ∀ H ∈ L 2 (Ω; R M ). (6.10) Our next step is to define the structure of subdifferential for the functional E 1 given by the rule (6.1). To do so, we make use of the following result (for the proof we refer to [24,Proposition 4.4]).
As immediately follows from (6.11), a vector field z can be formally identified with the quotient DC |DC| provided |DC| is nonzero and well defined at a given point x ∈ Ω. However, due to the Azellotti's theory of pairing, the correct interpretation of the quotient DC |DC| can be made through the equality (z, DC) = |DC|, where the field z ∈ L ∞ 2,div (Ω; R M ×2 ) is such that z L ∞ (Ω;R M ×2 ) 1 (see [12] for the details).
We are now in a position to derive the necessary conditions for a unique minimizer C rec ε ∈ Λ ⊂ BV (Ω; R M ) of constrained minimization problem (5.2). Since the functionals J i : L 2 (Ω; R M ) → R and G : L 2 (Ω; R M ) → R are convex, necessary conditions should employ the convex subdifferentials ∂ J i (C rec ε ) and ∂G (C rec ε ). As a result, utilizing Proposition 6.3 in [28], we arrive at the following result.
of convex functionals saying that we deduce the existence of elements Hence, in view of Propositions 6.1 and 6.2, we arrive at the announced relations (6.13)-(6.16).

Optimality Conditions for Brightness Reconstruction Problem
The main object of our consideration in this section is the following constrained minimization problem where the set of admissible solutions Ξ is defined in (3.8), and the weight multiplayer K (B 0 , u SAR ) is given by (4.12). By analogy with the previous section, we define the following functionals +∞, B ∈ L 2 (Ω) \ BV (Ω). Notice that all functional are well-defined on L 2 (Ω). In view of this, we extend the energy functional J to the entire set L 2 (Ω) by the rule Then it is clear that a minimizer B rec ∈ Ξ of (7.1) is also a minimizer of the modified problem because the inclusion B ∈ Ξ is equivalent to the consistency condition of the problem (7.1), i.e.

Schemes for Numerical Simulations
In this section we shortly describe the crucial steps of alternative minimization method that we suggest for the numerical simulations of the spectral indices by the damaged multi-band satellite optical images. Since the reconstruction problem is split up onto two independent constrained minimization problems (P 1 ) and (P ε ), we begin with the chromaticity recovery problem (3.11).
Due the operator splitting technique (see, for instance, [21,22]), we pass to the following regularized version of the problem (5.2) Find C 0 ε ∈ BV (Ω; R M ) and V 0 ε ∈ L 2 (Ω; R M ) such that F ε C 0 , V 0 = inf C,V F ε (C, V ) subject to V = C, |V | = 1, Following method of multipliers and gradients that was proposer in [25], we associate with constrained minimization problem (8.1) the following augmented Lagrangian functional Here, λ 1 = λ 1 (x) ∈ R M and λ 2 = λ 2 (x) ∈ R are the Lagrange multipliers associated to the constraints V = C and |V | = 1, respectively, and r 1 > 0 is the penalty parameter corresponding to the Lagrange multiplier λ 1 . Let at some iteration k ∈ N, C k , V k , λ k 1 , and λ k 2 be given. Then the idea is to exploit the alternating minimization approach to find the saddle points of the Lagrangian L (C, V , λ 1 , λ 2 ) iteratively. With that in mind, we propose the following iteration scheme (see [57,Section 2] for comparison) where τ 2ˆΩ C − C k 2 dx, with τ > 0, plays the role of the proximal term. So, the original chromaticity recovery problem can be reduced to the sequences of two unconstrained minimization problems and two multiplier updates. Moreover, their solutions can be obtained separately for each spectral channel. For the details of this procedure, its substantiation, and practical implementation, we refer to the recent paper [57]. As for the numerical scheme for the brightness reconstruction problem (7.1), we refere to [50], where some generalization of Chambolle's algorithm to the case of an H −1 -constrained minimization of the total variation in the case of T V -H −1 inpainting problem has been proposed (see [50,Section 2.3]).