In this work, vectors are represented with bold lowercase letters and the matrices, with bold uppercase characters. The expression means the space of real matrices of size 𝑁_𝑟 𝖷 𝑁_𝑐. In the vector column , the element placed in the 𝑛𝑟 th position is denoted by . In the matrices, the element placed in the 𝑛𝑟 th row and 𝑛𝑐 th column is denoted by. Additionally, the vector located in the 𝑛_𝑟 th row of matrix 𝑿 is denoted by and, in the same way, for the 𝑛_𝑐 th column, be, represent by . Let the hyperspectral matrix, where 𝑁_𝑝 is the number of pixels; N _b , the number of bands; and , the 𝑛𝑝 -th spectral vector. Due to the fact that the adjacent bands of spectral images are highly correlated, the columns of 𝐅 lie in a subspace with dimension 𝑁𝑠, such that 𝑁𝑠 ≪ 𝑁𝑏 [17]. Then, using the matrix form, the spectral image 𝐅 can be written as the combination of a low-rank parametrization and additive noise see (1)

where is the true signal of numerical rank 𝑁_𝑠 and is the additive noise matrix that is assumed to be distributed as the matrix normal distribution as (2)

with noise covariance . Traditionally, the noise in HSIs is modelled as additive Gaussian and i.i.d. with covariance [18].

Recently, noise has been assumed to be non-i.i.d., having a diagonal noise covariance matrix as ^{[17] [19]}. In order to reconvert the non-i.i.d. model into an i.i.d. one, Equation (1) is rewritten as (3).

Therefore, noise matrix is assumed to be i.i.d. in accordance with the matrix normal distribution with noise covariance . Specifically, if represents the so-called Frobenius norm of the matrix 𝐗, then, the solution to the constrained optimization problem by (4).

Where is the truncated SVD of When is estimated from, the true signal is . Since there is a correlation between bands, the spectral correlation can be removed successfully using multiple linear regression by estimating the diagonal matrix 𝚺 [17].

In accordance with the LMM described in [20], it is assumed that the spectral scene, for a given pixel, is generated by a linear combination of a small number of deterministic spectral signatures known as endmembers. The basic vectors may be interpretable as the mixing matrix, containing the spectral signatures of the endmembers, and as the abundance vector of the spectral vector From (1), the matrix form of LMM is given by (5).

Where is the matrix of abundances and is the matrix of endmembers. Due to physical considerations, the abundance coefficients should satisfy the non-negativity constraints (ANC). When the hyperspectral dataset is denoised, the proposed method only requires the initialization of the endmember matrix 𝚬. The pure pixel-based algorithms assume the presence, in the data, of at least one pure pixel per endmember ^[11]. The Vertex Component Analysis (VCA) is a pure pixel-based algorithm that iteratively projects data onto a direction orthogonal to the subspace spanned by the endmembers already determined ^[21]. However, the pure pixel assumption is a strong requirement that may not hold in many datasets. For that purpose, the Simplex Identification via Split Augmented Lagrangian (SISAL) method ^[22] is a state-of-the-art method that does not require the presence of pure pixels in the image.

The unmixing procedure can be formulated as a constrainednorm optimization problem, in which the observations are affected by Gaussian noise ^[23]. Furthermore, the abundance sum-to-one constraint (ASC) is defined as; then, the optimization problem is formulated as (6). This problem is solved using the SUnSAL TV algorithm ^[24] . Such formulation looks for the sparsest non-negative solution 𝚨 which explains data 𝚾, given endmembers 𝚬. However, traditional LMM approaches fail to fully account for the spectral variability associated with spatial changes. In this work, the LMM has been extended into an augmented linear mixing model by introducing a weighting vector introducing the weighinto (6), the bands that contribute the least to erroneous fraction estimates are prioritized, and the effect of unfavorable bands can be mitigated.

Figure 1 provides an overview of the weighting operation. Note that, when the signatures are normalized by the weight factor, they are very similar. The plot shows that the influence of spectral variability weight on original signatures is more beneficial for the unmixing problem. For that reason, the weighting factor is introduced in (6) to reduce erroneous abundances. The estimation becomes (7).

Fig. 1 Source: Authors’ own work.

Despite the fact that this regularizer is not convex, the corresponding function also has a simple closed form: the so-called hard-threshold function, where hard (·, a) denotes the component-wise application of the function . In this work, the SUNSAL TV algorithm was modified [24]. In (2), the soft operator is changed by a weighted hard thresholding operator. Additionally, the TV regularizer term is introduced to impose spatial consistency in hyperspectral sparse unmixing solutions.

Where is a vector extension of the non-isotropic TV and are regularization parameters.

Where is the number of non-zero entries in and is the weighted regularization parameter,commonly called the “norm” (although it is not a norm).

The proposal is summarized in Algorithm 1.

Algorithm 1: Hyperspectral Unmixing

// Input parameters //

1: Input:,

// Denoising procedure//

// (Section 2.2) //

2:𝚺→noise_estimation(F) // Alg. 1 [17] //

3:

4:

5:

// Endmember extraction //

// (Section 2.3) //

6: // Algorithm 1[22] //

// Abundance estimation //

// (Section 2.4) //

7:

8:

9:𝚨→ modified_sunsal_tv

// Algorithm 2 [24] //

// Output parameters //

10: Output:

This paper uses a dataset of an oil palm plantation acquired through hyperspectral airborne remote sensing in Colombia. The dataset was recorded and supplied by the Colombian company Quimbaya Aerial Services. These images were acquired in the visible and near infrared (VNIR) range using a HySpex VNIR-1600 hyperspectral camera (Norsk Elektro Optikk AS, Norway) over oil palm crops in the southeastern region of Colombia in 2017.

Figure 2 shows the RGB map of the study area, which is located in the Department of Caquetá, Colombia. The study area is a 130-hectare oil palm plantation and the HSIs were acquired through nine flights lines at an altitude of 430 m. Specifically, the dataset considered in this work has a spatial resolution of 0.6 x 0.6 m, and a spectral resolution of 3.7 nm in the spectral range of 415 - 990 nm. The dimensions of the hypercube are spectral bands. Fig. 2 shows three subsections in the scene, and each subsection includes various land cover types, e.g., oil palm plantation areas, other vegetation, bare land, buildings, and roads, etc. The first scene exhibits three materials (Oil Palm, Oil Palm Disease, and Grass); the second scene, four materials (Soil, Tree, Palm, and Grass); and the third scene, five materials (Roof, Soil, Grass, and Palm).

Fig. 2 Source: Authors’ own work.

This section explores the effect of selecting the parameter 𝜆_𝑇𝑉 and the number of features 𝑁_𝑠 on the performance of the proposed method for oil palm identification. In the inference experiments, such as unmixing, classification, and detection, the estimators have been improved when denoising techniques are applied as a preprocessing step. Fig. 3(a) is the original noisy spectral band 160, while Figs. 3(b) and 3(c) show the clean spectral band 160 when the rank of the matrix 𝚾 is 𝑁_𝑠 = 10 and 𝑁_𝑠 = 20, respectively. As can be seen in Fig. 3(d) (diagonal covariance matrix), the high noise is present in the 160 band, which can be observed in Fig. 3(e) that the pixel in the clean dataset is smooth with respect to the original pixel in the noisy dataset. The estimation of the number of endmembers is a critical stage in the unmixing. In this work, we consider two of the most widely used approaches: virtual dimensionality (VD) ^[25], which allows some flexibility in the estimation by having an additional input parameter that enables the control of the sensitivity of the method, and hyperspectral signal identification with minimum error (HySime)^[17], which offers advanced features when modeling the noise present in the HSIs. The average number of endmembers estimated by the HySime algorithm for the three scenes was 9 endmembers. However, based on knowledge of the actual terrain and each scene used in this work, the endmembers with the greatest abundance and interest in this study are shown (diseased oil palm, healthy oil palm, grass, soil, and roof).

Fig. 3 Source: Authors’ own work.

The SISAL regularization parameters were set to 𝜆 = 10, and the remaining parameters were set to 𝜏 = 1 and 𝜇 = 10^-4. Fig. 4 shows the endmember signatures estimated by using SISAL (first column) and VCA (second column) for each scene. Note that the results of endmembers estimated by VCA in the three scenarios are very similar to those of SISAL. However, when abundances are estimated, the difference between VCA and SISAL is noticeable. Figures 5 and 6 show the abundances using the proposed method for the first and second scene, respectively, when the endmembers matrix is initialized using VCA (first row) and SISAL (second row). In optimal conditions, the algorithm without the pure pixel assumption offers a better performance. In the abundance results in the second row in Fig. 6, the red points indicate the oil palms detected by the proposed method and the blue spaces, other types of objects. Note that the proposed method achieves great oil palm identification results using SISAL, as opposed to when the endmembers are initialized with VCA. Furthermore, Fig. 7 shows the effect of the parameter 𝜆_𝑇𝑉 on palm identification; the proposed method is applied to the dataset with three different values of 𝜆_𝑇𝑉 = 0.025, 0.01, 0.001 and 𝑁_𝑠 = 10. It can be seen that using 𝜆_𝑇𝑉 = 0.01 yields better results than other cases.

Fig. 4 Source: Authors’ own work.

Fig. 5 Source: Authors’ own work.

Fig. 6 Source: Authors’ own work.

Fig. 7 Source: Authors’ own work.

In this section, the proposed framework is compared with two state-of-art unmixing methods: SUnSAL ^[23] and SUnSAL TV ^[24]. Fig. 8 illustrates the performance of the proposed computational framework with the third scene. Comparisons with pure pixel estimation methods were conducted, and the obtained results shown the effectiveness of the proposed method to map oil palm plantations. It can be seen that this area is not homogeneous. Namely, it contains objects that are not oil palm trees, i.e., buildings, other trees, and shadows.

Fig. 8 Source: Authors’ own work.

Moreover, the SUnSAL and SUnSAL TV methods could merge some trees together due to the presence of spectral variability or when the distance between them is short. In that case, the regions obtained as palms and non-palms using the proposed method are quite evident.

In addition to that, an oil palm detection procedure was applied for evaluation purposes. After the abundance maps of the oil palm plantation were estimated, oil palm detection was performed using the Circle Hough Transform (CHT) algorithm ^[26], the Producer Accuracy (Pacc) and User Accuracy (Uacc) by (8).

Where 𝑁_{𝑝𝑎𝑙𝑚} is the total number of palm trees in the image; 𝑇_𝑝, the true positive detection indicating the number of detected palms trees; and F _p , the false negative detection indicating the number of palm trees that were not detected (missed). Fig. 9 shows some palm and no-palm templates highlighted with circles bythe CHT algorithm. In the first image, the white circles represent the true location of oil palms. In the second and third images, the white circles represent T _p , and the red circles represent F _p , which are estimated from abundance map by the proposed method and SUnSAL TV, respectively. For example, the oil palms missing in the detection are marked with red circles. Oil palm trees can be missing from the detection because of the presence of other objects or tree crowns close together. Table 1 illustrates the detection results using different measurements. Note that the results suggest that the proposed method is the most promising strategy as preprocessing step before oil palm mapping.

Fig. 9 Source: Authors’ own work.

Table 1

Source: Authors’ own work.