# Key-Frame Detection and Super-Resolution of Hyperspectral Video via Sparse-Based Cumulative Tensor Factorization

Authors: Ruofei Zhou; Gang Wang; Bo Li; Jinlong Wang; Tianzhu Liu; Chungang Liu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Thanks to the rapid development of hyperspectral sensors, hyperspectral videos (HSV) can now be collected with high temporal and spectral resolutions and utilized to handle invisible dynamic monitoring missions, such as chemical gas plume tracking. However, using such sequential large-scale data effectively is challenged, because the direct process of these data requires huge demands in terms of computational loads and memory. This paper presents a key-frame and target-detecting algorithm based on cumulative tensor CANDECOMP/PARAFAC (CP) factorization (CTCF) to select the frames where the target shows up, and a novel super-resolution (SR) method using sparse-based tensor Tucker factorization (STTF) is used to improve the spatial resolution. In the CTCF method, the HSV sequence is seen as cumulative tensors and the correlation of adjacent frames is exploited by applying CP tensor approximation. In the proposed STTF-based SR method, we consider the HSV frame as a third-order tensor; then, HSV frame super-resolution problem is transformed into estimations of the dictionaries along three dimensions and estimation of the core tensor. In order to promote sparse core tensors, a regularizer is incorporated to model the high spatial-spectral correlations. The estimations of the core tensor and the dictionaries along three dimensions are formulated as sparse-based Tucker factorizations of each HSV frame. Experimental results on real HSV data set demonstrate the superiority of the proposed CTCF and STTF algorithms over the comparative state-of-the-art target detection and SR approaches.


---

## Body

## 1. Introduction

Hyperspectral imaging has been one of the most popular research fields due to its ability of identifying the materials from very high spectral resolution and coverage. In the last decade, researchers focused on the processing and application of hyperspectral image (HSI), such as denoising [1, 2], feature extraction [3, 4], classification [5–11], detection [12–14], and super-resolution (fusion) [15–18]. In this section, researching of the latter two fields which are related to this paper will be briefly introduced.Basically, target detection is a kind of binary classifier with the purpose of labeling every image pixel as a target or background. In HSIs, pixels with a significantly different spectral signature from their neighboring background pixels are defined as spectral anomalies. Anomaly detectors are statistical or pattern recognition methods used to detect distinct pixels that differ from the background. It is worth mentioning that, in spectral anomaly detection approaches [19–22], such as Reed-Xiaoli (RX) algorithm [23], no prior information of the target spectral signature is assumed or used. However, we focus on the detection of invisible gas plumes in this paper, and the prior knowledge of the desired targets spectral characteristics is assumed to be known. In such cases, signature-based target detection algorithms are presented instead of anomaly detection. In these algorithms, the spectral characteristics of the target can be represented by a target subspace or a single target spectrum [24]. Likewise, the characteristics of background can be statistically expressed by a Gaussian distribution or a subspace defining the local or whole background statistics. As for this category, the matched subspace detector (MSD) method [25] is one of the most typical algorithms. In the MSD, the target pixel vectors are represented by a linear combination of the target spectral signature and the background spectral signature, which stand for the subspace target spectra and the subspace background spectra, respectively. Then, the generalized likelihood ratio test (GLRT) is applied, using projection matrices associated with the background subspace and the target-and-background subspace. At last, the comparison between the output of GLRT and a preset threshold makes a final decision about whether the target is absent or present. From pixel level to subpixel level, a single pixel may contain several distinct pure materials (endmembers), also known as the mixed pixel. The presence of mixed pixels is a tough problem caused by the low spatial resolution of HSIs. Accordingly, some unmixing approaches [26–28] have been designed to compute fractional abundance of endmembers. In [29], a hyperspectral unmixing approach based on constrained matrix factorization (CMF) was proposed. Unlike conventional methods, each column vector of endmember matrix is represented as a nonnegative linear combination of pixel spectra. After endmember matrix and the corresponding fractional abundance matrix are obtained by solving optimization problems, abundance map of the target endmember shows the detection result.As mentioned before, the HSIs often suffered from low spatial resolution. To acquire an HSI, the number of sun photons in each spectral band has to be greater than a minimum value, and the number of spectral bands is so huge in an HSI that the spatial resolution has to be sacrificed. Therefore, super-resolution (SR) techniques have aroused great interest in the last decade. Generally, the SR methods of HSI can be classified into four categories: Bayesian [30], component analysis [31], deep learning [32], and sparse representation. Due to the limited length of this paper, we focus on the introduction of sparse-based algorithms. In such HSI super-resolution schemes, images are expressed by dictionaries and corresponding sparse coefficients. On the basis of the spatial-spectral sparsity in the HSIs, the dictionaries and sparse coefficients are estimated jointly [33]. Huang et al. [34] introduced a fusion method of multispectral images (MSIs) with different spectral and spatial resolutions based on sparse matrix factorization. Akhtar et al. [35] presented an MSI-HSI fusion approach using sparse coding and Bayesian dictionary learning. Moreover, some algorithms based on matrix factorization [36–38] or unmixing [39] can also be regarded as the sparse representation schemes because the source images are decomposed into some basis and the corresponding coefficients. Yokoya et al. proposed a couple nonnegative matrix factorization (CNMF) [40] algorithm, where the unmixing techniques are employed to yield the endmember matrices and the high-resolution (HR) abundance matrices of HSI. In [41], Lanaras et al. suggested a joint scheme to solve the spectral unmixing problems. In [42], Zhang et al. fused the low-resolution (LR) HSI and HR-MSI based on the group spectral embedding and low-rank factorization.However, the matrix factorization based schemes cannot fully exploit the spatial-spectral correlations of the HSIs. It is believed that considering HSIs as tensors is better because an HSI can be naturally expressed as a third-order tensor. In this paper, a detection algorithm based on cumulative tensor CP factorization (CTCF) is proposed. The sequential HSV data is expressed as a four-dimensional (4D) cumulative tensor; factor matrices are obtained by decomposing original 4D tensor using CP factorization. When a new frame presents and is added to the time dimension of the original tensor, this 4D cumulative tensor is updated together with the factor matrices. Consequently, a CP tensor approximation of the new frame is computed by updated factor matrices and the fitness between the new frame and the approximation is calculated. After comparing the fitness to a preset threshold, we can make the decision that whether the new frame continues to be used to update the cumulative tensor or the new frame is the key-frame where the target presents. CTCF-based method exploits not only the spatial-spectral correlations of the HSIs by applying tensor model, but also the temporal correlation between adjacent frames of the HSV.On the other hand, tensor-based analysis has also been widely used in HSI super-resolution [43–45]. To the best of our knowledge, most of the SR algorithms enhance spatial resolution by fusing high-resolution MSI (HR-MSI) and low-resolution HSI (LR-HSI) from the same scene. Unfortunately, it is less practical in real application. In some situations, LR-HSI is the only data we have rather than both. In this paper, we suggest an SR algorithm using sparse-based tensor Tucker factorization (STTF). Inspired by the Tucker factorization and its related works, the HSV frames are represented as third-order tensors, which are approximated by the multiplication of the dictionaries along three dimensions (i.e., the dictionaries of the height mode, the width mode, and the spectral mode: they are named “three modes dictionaries” for short in the rest of this paper) and a core tensor. Then, the problem of SR is transformed into the estimations of the three modes dictionaries and estimation of the core tensor. Specifically, the spatial information is represented by the height mode dictionary and the width mode dictionary, the spectral information is represented by the spectral mode dictionary, and the correlations of the three modes dictionaries are modeled by the core tensor. HSIs are generally self-similar so that a sparse prior can be imposed on the core tensor; then, the estimations of the core tensor and three modes dictionaries are formulated as the STTF of the LR and HR HSV frames. In the iteration of STTF, core tensor and dictionaries are all updated and accurate estimates are yielded when convergence is achieved.The remainder of this paper is organized as follows.Section 2 presents the materials and methods, including the basic notations and preliminaries of tensor and tensor factorization, the proposed CTCF approach for key-frame detection, and the proposed STTF method for key-frame super-resolution problem. In Section 3, experimental results on real HSV and the discussions are given. The paper is summarized in Section 4 with ideas for future work along the path presented here.


## 2. Materials and Methods

### 2.1. Tensor Notations and Preliminaries

#### 2.1.1. Tensor Notations

In this paper, vectors are denoted by boldface lowercase lettersa,b,c,⋯, matrices are denoted by boldface capital letters A,B,C,⋯, and tensors are denoted by bold Euler script letters A,B,C,⋯. Generally, a tensor is a kind of multidimensional array, denoted by A∈ℝI1×I2×⋯×IN. Here, tensor A is an Nth-order tensor and In1≤n≤N is the dimension of the nth mode. Obviously, vectors are first-order tensors and matrices are second-order tensors. We use Ai1,⋯,in−1,:,in+1,⋯,iN to denote the mode-n fiber, which are vectors yielded from tensor A by changing index in with other indexes fixed. The mode-n unfolding matrix of tensor A is generated by placing all the mode-n fibers in a matrix as columns, denoted by An∈ℝIn×I1,⋯,In−1In+1,⋯,IN.An important calculation between a tensor and a matrix is then-mode product, which is defined as(1)F=A×nB,where B∈ℝJn×In and F∈ℝI1×⋯×In−1×Jn×In+1×⋯×IN. The elements of A are denoted by ai1i2,⋯,iN, so the elements of F are computed by(2)fi1,⋯,in−1jnin+1,⋯,iN=∑inai1,⋯,in−1inin+1,⋯,iNbjnin.Given the definition ofn-mode product, we can obtain(3)F=A×nB⇔Fn=BAn.For continuous multiplication of a tensor and matrices in distinct modes, the result is not affected by the multiplication order, described by(4)A×nB×mC=A×mC×nBm≠n.If the modes are equivalent, equation (4) is transformed into(5)A×nB×nC=A×nCB.Suppose thatEn∈ℝJn×In1≤n≤N is a collection of matrices; we define tensor G∈ℝJ1×J2×⋯×JN as(6)G=A×1E1×2E2⋯×NEN.The matricization form of equation (6) is presented by(7)g=EN⊗EN−1⊗,⋯,⊗E1a,where g=vecG∈ℝJ (J=∏n=1NJn) and a=vecA∈ℝI (I=∏n=1NIn) are vectors yielded by arranging the mode-1 fibers of the tensors G and A. The Kronecker product is denoted by symbol “⊗.”Moreover, given the tensorA, A0 represents the ℓ0-norm which equals the number of nonzero elements of A, A1=∑i1,i2,⋯,iNai1i2,⋯,iNdenotes the ℓ1-norm, and AF=∑i1,i2,⋯,iNai1i2,⋯,iN2 denotes the Frobenius norm.The definition of rank-one tensor is introduced at last. TheNth-order tensor A is rank-one if it can be written as the outer product of N vectors, i.e., A=a1∘a2∘⋯∘aN. The symbol “∘” denotes the vector outer product [46].


#### 2.1.2. Tensor Factorizations

CANDECOMP/PARAFAC (CP) factorization decomposes a tensor into a sum of component rank-one tensors [47]. For example, given a third-order tensor X∈ℝI×J×K, we may formulate it as(8)X≈∑r=1Rar∘br∘cr,where R is a positive integer and ar∈ℝI, br∈ℝJ, and cr∈ℝK (r=1,2,⋯,R). The element of tensor X can be computed by(9)xijk≈∑r=1Rairbjrckr,i=1,⋯,I,j=1,⋯,J,k=1,⋯,K.CP factorization is illustrated in Figure1.Figure 1
CP factorization of a third-order tensor.The factorization result can be expressed by factor matrices of three dimensions. Factor matrices refer to the combination of the vectors from the rank-one components; i.e.,(10)A=a1,a2,⋯,aR,B=b1,b2,⋯,bR,C=c1,c2,⋯,cR.Following [48], the CP model can be concisely represented as(11)X≈A,B,C≡∑r=1Rar∘br∘cr.On the basis of factor matrices, the mode-n unfolding matrices Xn (n=1,2,3) of X can be represented as(12)X1≈AC⊙BT,X2≈BC⊙AT,X3≈CB⊙AT,where the symbol “⊙” denotes the Khatri-Rao product [49]. In this way, loss functions can be modeled as the approximation of the mode-n unfolding matrices; then the factor matrices of CP factorization can be obtained by solving the corresponding optimization problem.Tucker factorization is another popular tensor decomposing approach [50]. It decomposes a tensor into a core tensor multiplied by a matrix along each mode. Thus, in the same case as above where X∈ℝI×J×K, the factorization can be described as(13)X≈Z×1A×2B×3C=∑p=1P∑q=1Q∑r=1Rzpqrap∘bq∘cr,where A∈ℝI×P, B∈ℝJ×Q, and C∈ℝK×R are factor matrices which can be regarded as the principal components in each mode. Therefore, Tucker factorization is a form of higher-order principal component analysis (PCA). Tensor Z∈ℝP×Q×R is the core tensor and its elements stand for the correlation level between the different components. Similar to (11), the Tucker model can be concisely represented by X≈Z;A,B,C. Elementwise equation (13) can be represented as(14)xijk≈∑p=1P∑q=1Q∑r=1Rzpqraipbjqckr,i=1,⋯,I,j=1,⋯,J,k=1,⋯,K.The Tucker factorization is illustrated in Figure2.Figure 2
Tucker factorization of a third-order tensor.


### 2.2. The Proposed CTCF-Based Detection Method

In this subsection, the optimization problem of updating factor matrix is presented, followed with the proposed cumulative tensor CP factorization (CTCF) of third-order tensors. It is then extended toNth-order tensors. The CTCF-based detection method is described in the end of this subsection with its flowchart shown in Figure 3.Figure 3
Flowchart of CTCF-based detection method.


#### 2.2.1. CP Tensor Approximation by Factor Matrices

Similar to equation (12), the mode-n unfolding matrix of X∈ℝI1×I2×⋯×IN can be approximated by factor matrices; i.e.,(15)Xn≈AnAN⊙⋯⊙An+1⊙An−1⊙⋯⊙A1T=An⊙i≠nNAiT,where the factor matrices A1,⋯,AN are obtained by CP factorization. The corresponding loss function is(16)ℒ=12Xn−An⊙i≠nNAiT2.The Alternating Least Squares (ALS) algorithm is often applied to obtain factor matrices by solving the following optimization problem:(17)An=argminAn12Xn−An⊙i≠nNAiT2.When the tensor updates, the new tensor can be computed by the updated factor matrices which are given by equation (17).


#### 2.2.2. CTCF of Third-Order Tensor

Generally, an image is a second-order tensor; then sequential images form a third-order tensor, i.e., a video, adding a temporal dimension on two spatial dimensions. When a new video frame presents and is added to the time dimension of the original tensor, it is defined as a three-dimensional (3D) cumulative tensor. With the number of new frames increasing, the 3D cumulative tensor updates frame by frame.In conventional CP tensor approximation, whenever a new frame of image is added in the time dimension, ALS algorithm needs to be reused to approximate the new cumulative tensor, which is a time consuming process. In addition, the temporal correlation between neighboring frames is not exploited in the decomposition of the cumulative tensor. This paper proposes CTCF to update the CP factorization of original cumulative tensor, obtain the updated factor matrices, and approximate the new frame.Given an original 3D cumulative tensorXori∈ℝI×J×Tori, the result of CP factorization is denoted by Xori≈Aori,Bori,Cori. When a new tensor Xnew∈ℝI×J×Tnew is added in the time dimension, the updated cumulative tensor is X∈ℝI×J×Tori+Tnew, of which the CP factorization appears as X≈A,B,C. We focus on obtaining A, B, and C by updating Aori, Bori, and Cori.The updating process is operated in an alternating way. Firstly, temporal dimensional factor matrixC is computed while factor matrices A and B are fixed; i.e.,(18)C=argminC12X3−CB⊙AT2=argminC12Xori3Xnew3−C1C2B⊙AT2=argminC12Xori3−C1B⊙ATXnew3−C2B⊙AT2,where C is divided into two terms. For A and B are fixed as Aori and Bori, the first row of (18) will be minimized if C1=Cori. To minimize the second row, according to (12), the optimal solution of C2 is Xnew3B⊙AT†, where the symbol “†” denotes Moore–Penrose pseudoinverse of the matrix [51]. So, C can be updated by adding Cnew which is represented by(19)C=CoriCnew=CoriXnew3B⊙AT†.Secondly, factor matrixA is computed while factor matrices B and C are fixed. Similar to 16, the loss function of estimating A is written as(20)ℒ=12X1−AC⊙BT2.Deriveℒ with respect to A; then, we have(21)∂ℒ∂A=X1C⊙B−AC⊙BTC⊙B.To simplifyequation (21), denote P=X1C⊙B and Q=C⊙BTC⊙B; thus, when ∂ℒ/∂A=0, we have A=PQ−1. According to [47], Q can be rewritten as(22)Q=C⊙BTC⊙B=CTC⊗BTB.For computingP, we also divide X1 and C into two terms; i.e.,(23)P=X1C⊙B=Xori1,Xnew1CoriCnew⊙B=Xori1,Xnew1Cori⊙BCnew⊙B=Xori1Cori⊙B+Xnew1Cnew⊙B.SinceB are fixed as Bori, the first term of equation (23) contains only the information of original tensor, which can be expressed by(24)Pori=Xori1Cori⊙B,so,equation (23) is rewritten as(25)P=Pori+Xnew1Cnew⊙B,Hence, P can be updated from Pori using mode-1 unfolding matrix of Xnew and factor matrix mentioned above Cnew. Generally, P is initialized by Xτ∈ℝI×J×τ, which is a small front part of Xori, and updated iteratively by (25). Analogously, the update process of Q can be represented by(26)Q=Qori+Cnew⊙BTCnew⊙B=Qori+CnewTCnew⊗BTB.The update ofA may be summarized as(27)P←P+Xnew1Cnew⊙B,Q←Q+CnewTCnew⊗BTB,A←PQ−1.Finally, the update of factor matrixB may likewise be expressed by(28)U←U+Xnew2Cnew⊙A,V←V+CnewTCnew⊗ATA,B←UV−1,where U=X2C⊙A and V=CTC⊗ATA.To make the process clearer, the proposed CTCF of third-order tensor is summarized by Algorithm1.Algorithm 1: CTCF of third-order tensor.
Input: original 3D cumulative tensor Xori∈ℝI×J×Tori≈Aori,Bori,Cori new tensor Xnew∈ℝI×J×TnewStep 1: new tensor is added in the time dimension andX∈ℝI×J×Tori+Tnew is obtainedStep 2: decomposeX by CP factorization X≈A,B,CStep 3: updateC by (19), with A and B are fixedStep 4: updateA by (27), with B and C are fixedStep 5: updateB by (28), with A and C are fixedStep 6: estimateX by updated A, B and COutput: approximation of updated cumulative tensor X


#### 2.2.3. CTCF ofNth-Order Tensor

On the basis ofSection 2.2.2, we try to extend CTCF to higher-order tensors. Suppose an N-dimensional cumulative tensor Xori∈ℝI1×I2×⋯×IN−1×Tori where the last dimension is temporal dimension. The CP factorization of Xori is represented as Xori≈Aori1,⋯,AoriN−1,AoriN. When a new tensor Xori∈ℝI1×I2×⋯×IN−1×Tnew is added in the time dimension, the updated cumulative tensor is X∈ℝI1×I2×⋯×IN−1×Tori+Tnew, of which the CP factorization is denoted by X≈A1,⋯,AN−1,AN.Similar toSection 2.2.2, temporal dimensional factor matrix AN is firstly updated with other N−1 matrices fixed. Like 17, the optimization problem of estimating AN is formulated by(29)AN=argminAN12XN−AN⊙N−1AiT2.We also separate original part from new added part; i.e.,(30)AN=AoriNAnewN=AoriNXnew⊙N−1AiT†,i=1,2,⋯,N−1.The original part is minimized by fixing the firstN−1 factor matrix and the new part is updated by Xnew.The updates of nontemporal dimensional factor matricesAn (n∈1,N−1) may refer to the ones of factor matrices A and B in Section 2.2.2. The loss function ℒn of estimating An is the same as 16. Let ∂ℒn/∂An=0 and introduce matrices Pn and Qn; the update of An may be summarized as(31)Pn←Pn+XnewnAnewN⊙Kn,Qn←Qn+AnewNTAnewN⊗Rn,An←PnQn−1,where Kn=⊙i≠nN−1Ai and Rn=⊗i≠nN−1AiTAi.


#### 2.2.4. CTCF-Based Detection Method

In HSV, the sequential data is expressed as a 4D cumulative tensor; the temporal dimension increases with new frames are added in. Whenever a new frame presents, the results of original cumulative tensor CP factorization are updated to obtain the factor matrices of the new cumulative tensor, and the CP tensor approximation of the newly added frame is obtained at the same time. If the target is absent, the CP tensor approximation will lead to a small error, since the background information is similar between adjacent frames. On the contrary, if the error is large, the target is likely to present. We define the fitness between the new frame and its approximation in34. If the fitness is smaller than the threshold, the target is supposed to appear in the new frame. Otherwise, the new frame is added in the temporal dimension and used to update original cumulative tensor.The original 4D cumulative tensor is denoted byXori∈ℝI1×I2×I3×n; n denotes the frame number of initial video. The factor matrices of four dimensions are represented as(32)Xori≈Aori1,Aori2,Aori3,Aori4,where Aori1∈ℝI1×R, Aori2∈ℝI2×R, Aori3∈ℝI3×R, and Aori4∈ℝn×R and R denotes the number of component rank-one tensors in CP factorization. When a new frame Xnew∈ℝI1×I2×I3 is added in the temporal dimension of original 4D cumulative tensor, the 4D cumulative tensor is updated and denoted by X∈ℝI1×I2×I3×n+1. The factor matrices of X are expressed by(33)X≈A1,A2,A3,A4,where A1∈ℝI1×R, A2∈ℝI2×R, A3∈ℝI3×R, and A4∈ℝn+1×R. Based on Section 2.2.3, we estimate A1∼A4 and obtain the approximation of X and X˜new, where X˜new≈A1,A2,A3,A4n+1,:. Actually, it is the specific case when N=4.We define thefitness (Xnew, X˜new) as(34)fitnessXnew,X˜new=1−Xnew−X˜newFXnewF.If the target does not appear, the approximation error is small and the result of fitness is large. Given a preset thresholdη, when fitnessXnew,X˜new>η, i.e., the fitness is larger than η, we decide that the target is absent. Then, the nontarget frame is added in temporal dimension and the updated 4D cumulative tensor becomes the new original 4D cumulative tensor, which can be expressed as(35)Xori⇐X,Aorin⇐An,n∈1,4.If the target appears, the approximation error is large and the fitness is smaller thanη. The residual of Xnew and X˜new is the approximation of the target tensor; i.e.,(36)Xtar≈Xnew−X˜new.The target of each frame will be shown in 2D form by taking the maximum value of every spectrum. In this way, the proposed CTCF-based detection method can extract not only the key-frames where the target presents, but also the approximate region of target in every key-frame. The flowchart of the proposed method is shown in Figure3. In Section 3, experiments on real HSV data are conducted and the proposed method is compared with some representative techniques.


### 2.3. The Proposed STTF-Based Super-Resolution Method

InSection 2.2, we present an approach to detect the frames where the target appears in HSV and the approximate region of the target. However, as discussed in Section 1, there has to be a tradeoff between spectral resolution and the spatial resolution in HSI imaging systems [52]. The spatial resolution is always low since high spectral resolution is required in HSIs and HSV. So, we are interested in improving the spatial resolution of targets after the detecting process. Instead of fusing HR-MSI and LR-HSI, we try to handle the target SR problem by what we have got, which is more practical in real cases.


#### 2.3.1. Problem Formulation

In this subsection, HSIs are represented as 3D tensors with three indexes (H,W,S), which stand for the height, width, and spectral modes. X∈ℝH×W×S denotes the HR-HSI and the LR-HSI is denoted by Y∈ℝh×w×S, where W>w and H>h. The goal is to estimate X from Y.There are two significant characteristics of HR-HSIs [53]: the first one is that spectral vectors can be well approximated in low dimensional subspaces, and the second one is that HSIs are spatially self-similar. This means that sparsity exists in both spectral and spatial dimensions. Inspired by sparse representation [54], the low dimensionality in spectral domain gives the possibility to form a spectral mode dictionary S with few nonzero atoms; the self-similarities in spatial domain guarantee the sparse representations of the height and width modes with spatial dictionaries H and W. In this way, the conventional Tucker factorization is transformed into the multiplication of the core tensor and three modes dictionaries. The factorization is illustrated in Figure 4. The HR-HSI is represented as(37)X=Z×1H×2W×3S,where H∈ℝH×zh, W∈ℝW×zw, and S∈ℝS×zs. The variables zh, zw, and zs denote the atoms (i.e., the number of columns) of H, W, and S, respectively. The core tensor Z contains the coefficients of X over three modes dictionaries. We can see that 37 incorporates the information of separated modes into a unified framework.Figure 4
Sparse-based tensor Tucker factorization of an HSI.The LR key-frame of HSVY can be seen as the spatially downsampled version of HR-HSI X, which is written as(38)Y=X×1D1×2D2,where D1∈ℝh×H and D2∈ℝw×W are downsampling matrices of the height and width modes. Substituting 37 into (38), Y is represented by(39)Y=Z×1D1H×2D2W×3S=Z×1H∗×2W∗×3S,where H∗=D1H∈ℝh×zh and W∗=D2W∈ℝw×zw denotes the downsampled dictionary of height and width modes. To recover X, we focus on estimating the dictionaries H, W, and S and the core tensor Z.


#### 2.3.2. The Proposed STTF-Based SR Algorithm

SinceY is a downsampled version, recovering X from Y is a typical inverse problem, which is badly ill-posed. So, some prior knowledge of X is needed to regularize the super-resolution problem. In HSI processing, the spectral sparsity is a widespread regularizer applied to solve varieties of ill-posed problems [55–58]. In such regularization, spectral vectors are linearly combined by a small quantity of different spectral signatures. However, these schemes only take advantage of the sparsity of the spectral domain. In the proposed algorithm, taking into account the HSI self-similarity, sparsity regularization is extended to the spatial domain by exploiting the sparse-based tensor Tucker factorization (STTF). In STTF, the HR-HSI performs a united sparse representation of the core tensor and three modes dictionaries.On the basis of equation (39), the HSV frame super-resolution is formulated as a constrained least-squares optimization problem:(40)minH,W,S,ZY−Z×1H∗×2W∗×3SF2,s.t.Z0≤θ,where ⋅F represents the Frobenius norm and θ denotes the number of maximum nonzero elements of Z. Because of the ℓ0-norm constraint, equation (40) is nonconvex. To make the optimization processable, the ℓ0-norm is replaced by the ℓ1-norm and 40 is transformed into an unconstrained version:(41)minH,W,S,ZY−Z×1H∗×2W∗×3SF2+λZ1,where λ is the parameter of sparse regularizer. Equation (41) is also nonconvex, and the solutions of H, W, and S and Z are not unique. Nonetheless, if we focus on only one variable with other variables fixed, the objective function in equation (41) is convex. Inspired by [59, 60], equation (41) can be solved by proximal alternating optimization scheme, which is guaranteed to reach convergence in a particular situation. Concretely, H, W, S, and Z are updated iteratively by(42)H=argminHfH,W,S,Z+αH−HpreF2,W=argminWfH,W,S,Z+αW−WpreF2,S=argminSfH,W,S,Z+αS−SpreF2,Z=argminZfH,W,S,Z+αZ−ZpreF2,where ⋅pre denotes the previous estimation in the last iteration and α denotes a positive number. Equation (41) defines the object function fH,W,S,Z. The optimizations of H, W, S, and Z will be presented detailedly in the appendix. The conjugate gradient (CG) method [61] and the alternating direction method of multipliers (ADMM) [62] will be used in the optimizations.


#### 2.3.3. Initialization of the Proposed Method

Since the optimization problem in (41) is nonconvex, the solution would result in poor local minima if we set the initialization carelessly. In this paper, we initialize the spatial dictionaries H∗ and W∗ from Y1 and Y2 dictionary-updates-cycles KSVD (DUC-KSVD) [63]; this method can promote sparse representations. Then, initialization of spectral dictionary S is accomplished by simplex identification split augmented Lagrangian (SISAL) algorithm [64]; this approach can efficiently identify a minimum unit that contains the spectral vectors.The proposed STTF-based SR algorithm is summarized in Algorithm2.Algorithm 2: STTF-based SR method.
Input: LR-HSI YInitializeS with SISALInitializeH∗ and W∗ with DUC-KSVDInitializeZ with (39)while no convergence doStep 1: updateH by solving (A.3) with CGH∗=D1H,Hpre=HStep 2: updateW by solving (A.6) with CGW∗=D2W,Wpre=WStep 3: updateS by solving (A.9) with CGSpre=SStep 4: updateZ by solving (A.15) with CGZpre=Zend whileEstimateX by (37)Output: HR-HSI X


## 2.1. Tensor Notations and Preliminaries

### 2.1.1. Tensor Notations

In this paper, vectors are denoted by boldface lowercase lettersa,b,c,⋯, matrices are denoted by boldface capital letters A,B,C,⋯, and tensors are denoted by bold Euler script letters A,B,C,⋯. Generally, a tensor is a kind of multidimensional array, denoted by A∈ℝI1×I2×⋯×IN. Here, tensor A is an Nth-order tensor and In1≤n≤N is the dimension of the nth mode. Obviously, vectors are first-order tensors and matrices are second-order tensors. We use Ai1,⋯,in−1,:,in+1,⋯,iN to denote the mode-n fiber, which are vectors yielded from tensor A by changing index in with other indexes fixed. The mode-n unfolding matrix of tensor A is generated by placing all the mode-n fibers in a matrix as columns, denoted by An∈ℝIn×I1,⋯,In−1In+1,⋯,IN.An important calculation between a tensor and a matrix is then-mode product, which is defined as(1)F=A×nB,where B∈ℝJn×In and F∈ℝI1×⋯×In−1×Jn×In+1×⋯×IN. The elements of A are denoted by ai1i2,⋯,iN, so the elements of F are computed by(2)fi1,⋯,in−1jnin+1,⋯,iN=∑inai1,⋯,in−1inin+1,⋯,iNbjnin.Given the definition ofn-mode product, we can obtain(3)F=A×nB⇔Fn=BAn.For continuous multiplication of a tensor and matrices in distinct modes, the result is not affected by the multiplication order, described by(4)A×nB×mC=A×mC×nBm≠n.If the modes are equivalent, equation (4) is transformed into(5)A×nB×nC=A×nCB.Suppose thatEn∈ℝJn×In1≤n≤N is a collection of matrices; we define tensor G∈ℝJ1×J2×⋯×JN as(6)G=A×1E1×2E2⋯×NEN.The matricization form of equation (6) is presented by(7)g=EN⊗EN−1⊗,⋯,⊗E1a,where g=vecG∈ℝJ (J=∏n=1NJn) and a=vecA∈ℝI (I=∏n=1NIn) are vectors yielded by arranging the mode-1 fibers of the tensors G and A. The Kronecker product is denoted by symbol “⊗.”Moreover, given the tensorA, A0 represents the ℓ0-norm which equals the number of nonzero elements of A, A1=∑i1,i2,⋯,iNai1i2,⋯,iNdenotes the ℓ1-norm, and AF=∑i1,i2,⋯,iNai1i2,⋯,iN2 denotes the Frobenius norm.The definition of rank-one tensor is introduced at last. TheNth-order tensor A is rank-one if it can be written as the outer product of N vectors, i.e., A=a1∘a2∘⋯∘aN. The symbol “∘” denotes the vector outer product [46].


### 2.1.2. Tensor Factorizations

CANDECOMP/PARAFAC (CP) factorization decomposes a tensor into a sum of component rank-one tensors [47]. For example, given a third-order tensor X∈ℝI×J×K, we may formulate it as(8)X≈∑r=1Rar∘br∘cr,where R is a positive integer and ar∈ℝI, br∈ℝJ, and cr∈ℝK (r=1,2,⋯,R). The element of tensor X can be computed by(9)xijk≈∑r=1Rairbjrckr,i=1,⋯,I,j=1,⋯,J,k=1,⋯,K.CP factorization is illustrated in Figure1.Figure 1
CP factorization of a third-order tensor.The factorization result can be expressed by factor matrices of three dimensions. Factor matrices refer to the combination of the vectors from the rank-one components; i.e.,(10)A=a1,a2,⋯,aR,B=b1,b2,⋯,bR,C=c1,c2,⋯,cR.Following [48], the CP model can be concisely represented as(11)X≈A,B,C≡∑r=1Rar∘br∘cr.On the basis of factor matrices, the mode-n unfolding matrices Xn (n=1,2,3) of X can be represented as(12)X1≈AC⊙BT,X2≈BC⊙AT,X3≈CB⊙AT,where the symbol “⊙” denotes the Khatri-Rao product [49]. In this way, loss functions can be modeled as the approximation of the mode-n unfolding matrices; then the factor matrices of CP factorization can be obtained by solving the corresponding optimization problem.Tucker factorization is another popular tensor decomposing approach [50]. It decomposes a tensor into a core tensor multiplied by a matrix along each mode. Thus, in the same case as above where X∈ℝI×J×K, the factorization can be described as(13)X≈Z×1A×2B×3C=∑p=1P∑q=1Q∑r=1Rzpqrap∘bq∘cr,where A∈ℝI×P, B∈ℝJ×Q, and C∈ℝK×R are factor matrices which can be regarded as the principal components in each mode. Therefore, Tucker factorization is a form of higher-order principal component analysis (PCA). Tensor Z∈ℝP×Q×R is the core tensor and its elements stand for the correlation level between the different components. Similar to (11), the Tucker model can be concisely represented by X≈Z;A,B,C. Elementwise equation (13) can be represented as(14)xijk≈∑p=1P∑q=1Q∑r=1Rzpqraipbjqckr,i=1,⋯,I,j=1,⋯,J,k=1,⋯,K.The Tucker factorization is illustrated in Figure2.Figure 2
Tucker factorization of a third-order tensor.


## 2.1.1. Tensor Notations

In this paper, vectors are denoted by boldface lowercase lettersa,b,c,⋯, matrices are denoted by boldface capital letters A,B,C,⋯, and tensors are denoted by bold Euler script letters A,B,C,⋯. Generally, a tensor is a kind of multidimensional array, denoted by A∈ℝI1×I2×⋯×IN. Here, tensor A is an Nth-order tensor and In1≤n≤N is the dimension of the nth mode. Obviously, vectors are first-order tensors and matrices are second-order tensors. We use Ai1,⋯,in−1,:,in+1,⋯,iN to denote the mode-n fiber, which are vectors yielded from tensor A by changing index in with other indexes fixed. The mode-n unfolding matrix of tensor A is generated by placing all the mode-n fibers in a matrix as columns, denoted by An∈ℝIn×I1,⋯,In−1In+1,⋯,IN.An important calculation between a tensor and a matrix is then-mode product, which is defined as(1)F=A×nB,where B∈ℝJn×In and F∈ℝI1×⋯×In−1×Jn×In+1×⋯×IN. The elements of A are denoted by ai1i2,⋯,iN, so the elements of F are computed by(2)fi1,⋯,in−1jnin+1,⋯,iN=∑inai1,⋯,in−1inin+1,⋯,iNbjnin.Given the definition ofn-mode product, we can obtain(3)F=A×nB⇔Fn=BAn.For continuous multiplication of a tensor and matrices in distinct modes, the result is not affected by the multiplication order, described by(4)A×nB×mC=A×mC×nBm≠n.If the modes are equivalent, equation (4) is transformed into(5)A×nB×nC=A×nCB.Suppose thatEn∈ℝJn×In1≤n≤N is a collection of matrices; we define tensor G∈ℝJ1×J2×⋯×JN as(6)G=A×1E1×2E2⋯×NEN.The matricization form of equation (6) is presented by(7)g=EN⊗EN−1⊗,⋯,⊗E1a,where g=vecG∈ℝJ (J=∏n=1NJn) and a=vecA∈ℝI (I=∏n=1NIn) are vectors yielded by arranging the mode-1 fibers of the tensors G and A. The Kronecker product is denoted by symbol “⊗.”Moreover, given the tensorA, A0 represents the ℓ0-norm which equals the number of nonzero elements of A, A1=∑i1,i2,⋯,iNai1i2,⋯,iNdenotes the ℓ1-norm, and AF=∑i1,i2,⋯,iNai1i2,⋯,iN2 denotes the Frobenius norm.The definition of rank-one tensor is introduced at last. TheNth-order tensor A is rank-one if it can be written as the outer product of N vectors, i.e., A=a1∘a2∘⋯∘aN. The symbol “∘” denotes the vector outer product [46].


## 2.1.2. Tensor Factorizations

CANDECOMP/PARAFAC (CP) factorization decomposes a tensor into a sum of component rank-one tensors [47]. For example, given a third-order tensor X∈ℝI×J×K, we may formulate it as(8)X≈∑r=1Rar∘br∘cr,where R is a positive integer and ar∈ℝI, br∈ℝJ, and cr∈ℝK (r=1,2,⋯,R). The element of tensor X can be computed by(9)xijk≈∑r=1Rairbjrckr,i=1,⋯,I,j=1,⋯,J,k=1,⋯,K.CP factorization is illustrated in Figure1.Figure 1
CP factorization of a third-order tensor.The factorization result can be expressed by factor matrices of three dimensions. Factor matrices refer to the combination of the vectors from the rank-one components; i.e.,(10)A=a1,a2,⋯,aR,B=b1,b2,⋯,bR,C=c1,c2,⋯,cR.Following [48], the CP model can be concisely represented as(11)X≈A,B,C≡∑r=1Rar∘br∘cr.On the basis of factor matrices, the mode-n unfolding matrices Xn (n=1,2,3) of X can be represented as(12)X1≈AC⊙BT,X2≈BC⊙AT,X3≈CB⊙AT,where the symbol “⊙” denotes the Khatri-Rao product [49]. In this way, loss functions can be modeled as the approximation of the mode-n unfolding matrices; then the factor matrices of CP factorization can be obtained by solving the corresponding optimization problem.Tucker factorization is another popular tensor decomposing approach [50]. It decomposes a tensor into a core tensor multiplied by a matrix along each mode. Thus, in the same case as above where X∈ℝI×J×K, the factorization can be described as(13)X≈Z×1A×2B×3C=∑p=1P∑q=1Q∑r=1Rzpqrap∘bq∘cr,where A∈ℝI×P, B∈ℝJ×Q, and C∈ℝK×R are factor matrices which can be regarded as the principal components in each mode. Therefore, Tucker factorization is a form of higher-order principal component analysis (PCA). Tensor Z∈ℝP×Q×R is the core tensor and its elements stand for the correlation level between the different components. Similar to (11), the Tucker model can be concisely represented by X≈Z;A,B,C. Elementwise equation (13) can be represented as(14)xijk≈∑p=1P∑q=1Q∑r=1Rzpqraipbjqckr,i=1,⋯,I,j=1,⋯,J,k=1,⋯,K.The Tucker factorization is illustrated in Figure2.Figure 2
Tucker factorization of a third-order tensor.


## 2.2. The Proposed CTCF-Based Detection Method

In this subsection, the optimization problem of updating factor matrix is presented, followed with the proposed cumulative tensor CP factorization (CTCF) of third-order tensors. It is then extended toNth-order tensors. The CTCF-based detection method is described in the end of this subsection with its flowchart shown in Figure 3.Figure 3
Flowchart of CTCF-based detection method.


### 2.2.1. CP Tensor Approximation by Factor Matrices

Similar to equation (12), the mode-n unfolding matrix of X∈ℝI1×I2×⋯×IN can be approximated by factor matrices; i.e.,(15)Xn≈AnAN⊙⋯⊙An+1⊙An−1⊙⋯⊙A1T=An⊙i≠nNAiT,where the factor matrices A1,⋯,AN are obtained by CP factorization. The corresponding loss function is(16)ℒ=12Xn−An⊙i≠nNAiT2.The Alternating Least Squares (ALS) algorithm is often applied to obtain factor matrices by solving the following optimization problem:(17)An=argminAn12Xn−An⊙i≠nNAiT2.When the tensor updates, the new tensor can be computed by the updated factor matrices which are given by equation (17).


### 2.2.2. CTCF of Third-Order Tensor

Generally, an image is a second-order tensor; then sequential images form a third-order tensor, i.e., a video, adding a temporal dimension on two spatial dimensions. When a new video frame presents and is added to the time dimension of the original tensor, it is defined as a three-dimensional (3D) cumulative tensor. With the number of new frames increasing, the 3D cumulative tensor updates frame by frame.In conventional CP tensor approximation, whenever a new frame of image is added in the time dimension, ALS algorithm needs to be reused to approximate the new cumulative tensor, which is a time consuming process. In addition, the temporal correlation between neighboring frames is not exploited in the decomposition of the cumulative tensor. This paper proposes CTCF to update the CP factorization of original cumulative tensor, obtain the updated factor matrices, and approximate the new frame.Given an original 3D cumulative tensorXori∈ℝI×J×Tori, the result of CP factorization is denoted by Xori≈Aori,Bori,Cori. When a new tensor Xnew∈ℝI×J×Tnew is added in the time dimension, the updated cumulative tensor is X∈ℝI×J×Tori+Tnew, of which the CP factorization appears as X≈A,B,C. We focus on obtaining A, B, and C by updating Aori, Bori, and Cori.The updating process is operated in an alternating way. Firstly, temporal dimensional factor matrixC is computed while factor matrices A and B are fixed; i.e.,(18)C=argminC12X3−CB⊙AT2=argminC12Xori3Xnew3−C1C2B⊙AT2=argminC12Xori3−C1B⊙ATXnew3−C2B⊙AT2,where C is divided into two terms. For A and B are fixed as Aori and Bori, the first row of (18) will be minimized if C1=Cori. To minimize the second row, according to (12), the optimal solution of C2 is Xnew3B⊙AT†, where the symbol “†” denotes Moore–Penrose pseudoinverse of the matrix [51]. So, C can be updated by adding Cnew which is represented by(19)C=CoriCnew=CoriXnew3B⊙AT†.Secondly, factor matrixA is computed while factor matrices B and C are fixed. Similar to 16, the loss function of estimating A is written as(20)ℒ=12X1−AC⊙BT2.Deriveℒ with respect to A; then, we have(21)∂ℒ∂A=X1C⊙B−AC⊙BTC⊙B.To simplifyequation (21), denote P=X1C⊙B and Q=C⊙BTC⊙B; thus, when ∂ℒ/∂A=0, we have A=PQ−1. According to [47], Q can be rewritten as(22)Q=C⊙BTC⊙B=CTC⊗BTB.For computingP, we also divide X1 and C into two terms; i.e.,(23)P=X1C⊙B=Xori1,Xnew1CoriCnew⊙B=Xori1,Xnew1Cori⊙BCnew⊙B=Xori1Cori⊙B+Xnew1Cnew⊙B.SinceB are fixed as Bori, the first term of equation (23) contains only the information of original tensor, which can be expressed by(24)Pori=Xori1Cori⊙B,so,equation (23) is rewritten as(25)P=Pori+Xnew1Cnew⊙B,Hence, P can be updated from Pori using mode-1 unfolding matrix of Xnew and factor matrix mentioned above Cnew. Generally, P is initialized by Xτ∈ℝI×J×τ, which is a small front part of Xori, and updated iteratively by (25). Analogously, the update process of Q can be represented by(26)Q=Qori+Cnew⊙BTCnew⊙B=Qori+CnewTCnew⊗BTB.The update ofA may be summarized as(27)P←P+Xnew1Cnew⊙B,Q←Q+CnewTCnew⊗BTB,A←PQ−1.Finally, the update of factor matrixB may likewise be expressed by(28)U←U+Xnew2Cnew⊙A,V←V+CnewTCnew⊗ATA,B←UV−1,where U=X2C⊙A and V=CTC⊗ATA.To make the process clearer, the proposed CTCF of third-order tensor is summarized by Algorithm1.Algorithm 1: CTCF of third-order tensor.
Input: original 3D cumulative tensor Xori∈ℝI×J×Tori≈Aori,Bori,Cori new tensor Xnew∈ℝI×J×TnewStep 1: new tensor is added in the time dimension andX∈ℝI×J×Tori+Tnew is obtainedStep 2: decomposeX by CP factorization X≈A,B,CStep 3: updateC by (19), with A and B are fixedStep 4: updateA by (27), with B and C are fixedStep 5: updateB by (28), with A and C are fixedStep 6: estimateX by updated A, B and COutput: approximation of updated cumulative tensor X


### 2.2.3. CTCF ofNth-Order Tensor

On the basis ofSection 2.2.2, we try to extend CTCF to higher-order tensors. Suppose an N-dimensional cumulative tensor Xori∈ℝI1×I2×⋯×IN−1×Tori where the last dimension is temporal dimension. The CP factorization of Xori is represented as Xori≈Aori1,⋯,AoriN−1,AoriN. When a new tensor Xori∈ℝI1×I2×⋯×IN−1×Tnew is added in the time dimension, the updated cumulative tensor is X∈ℝI1×I2×⋯×IN−1×Tori+Tnew, of which the CP factorization is denoted by X≈A1,⋯,AN−1,AN.Similar toSection 2.2.2, temporal dimensional factor matrix AN is firstly updated with other N−1 matrices fixed. Like 17, the optimization problem of estimating AN is formulated by(29)AN=argminAN12XN−AN⊙N−1AiT2.We also separate original part from new added part; i.e.,(30)AN=AoriNAnewN=AoriNXnew⊙N−1AiT†,i=1,2,⋯,N−1.The original part is minimized by fixing the firstN−1 factor matrix and the new part is updated by Xnew.The updates of nontemporal dimensional factor matricesAn (n∈1,N−1) may refer to the ones of factor matrices A and B in Section 2.2.2. The loss function ℒn of estimating An is the same as 16. Let ∂ℒn/∂An=0 and introduce matrices Pn and Qn; the update of An may be summarized as(31)Pn←Pn+XnewnAnewN⊙Kn,Qn←Qn+AnewNTAnewN⊗Rn,An←PnQn−1,where Kn=⊙i≠nN−1Ai and Rn=⊗i≠nN−1AiTAi.


### 2.2.4. CTCF-Based Detection Method

In HSV, the sequential data is expressed as a 4D cumulative tensor; the temporal dimension increases with new frames are added in. Whenever a new frame presents, the results of original cumulative tensor CP factorization are updated to obtain the factor matrices of the new cumulative tensor, and the CP tensor approximation of the newly added frame is obtained at the same time. If the target is absent, the CP tensor approximation will lead to a small error, since the background information is similar between adjacent frames. On the contrary, if the error is large, the target is likely to present. We define the fitness between the new frame and its approximation in34. If the fitness is smaller than the threshold, the target is supposed to appear in the new frame. Otherwise, the new frame is added in the temporal dimension and used to update original cumulative tensor.The original 4D cumulative tensor is denoted byXori∈ℝI1×I2×I3×n; n denotes the frame number of initial video. The factor matrices of four dimensions are represented as(32)Xori≈Aori1,Aori2,Aori3,Aori4,where Aori1∈ℝI1×R, Aori2∈ℝI2×R, Aori3∈ℝI3×R, and Aori4∈ℝn×R and R denotes the number of component rank-one tensors in CP factorization. When a new frame Xnew∈ℝI1×I2×I3 is added in the temporal dimension of original 4D cumulative tensor, the 4D cumulative tensor is updated and denoted by X∈ℝI1×I2×I3×n+1. The factor matrices of X are expressed by(33)X≈A1,A2,A3,A4,where A1∈ℝI1×R, A2∈ℝI2×R, A3∈ℝI3×R, and A4∈ℝn+1×R. Based on Section 2.2.3, we estimate A1∼A4 and obtain the approximation of X and X˜new, where X˜new≈A1,A2,A3,A4n+1,:. Actually, it is the specific case when N=4.We define thefitness (Xnew, X˜new) as(34)fitnessXnew,X˜new=1−Xnew−X˜newFXnewF.If the target does not appear, the approximation error is small and the result of fitness is large. Given a preset thresholdη, when fitnessXnew,X˜new>η, i.e., the fitness is larger than η, we decide that the target is absent. Then, the nontarget frame is added in temporal dimension and the updated 4D cumulative tensor becomes the new original 4D cumulative tensor, which can be expressed as(35)Xori⇐X,Aorin⇐An,n∈1,4.If the target appears, the approximation error is large and the fitness is smaller thanη. The residual of Xnew and X˜new is the approximation of the target tensor; i.e.,(36)Xtar≈Xnew−X˜new.The target of each frame will be shown in 2D form by taking the maximum value of every spectrum. In this way, the proposed CTCF-based detection method can extract not only the key-frames where the target presents, but also the approximate region of target in every key-frame. The flowchart of the proposed method is shown in Figure3. In Section 3, experiments on real HSV data are conducted and the proposed method is compared with some representative techniques.


## 2.2.1. CP Tensor Approximation by Factor Matrices

Similar to equation (12), the mode-n unfolding matrix of X∈ℝI1×I2×⋯×IN can be approximated by factor matrices; i.e.,(15)Xn≈AnAN⊙⋯⊙An+1⊙An−1⊙⋯⊙A1T=An⊙i≠nNAiT,where the factor matrices A1,⋯,AN are obtained by CP factorization. The corresponding loss function is(16)ℒ=12Xn−An⊙i≠nNAiT2.The Alternating Least Squares (ALS) algorithm is often applied to obtain factor matrices by solving the following optimization problem:(17)An=argminAn12Xn−An⊙i≠nNAiT2.When the tensor updates, the new tensor can be computed by the updated factor matrices which are given by equation (17).


## 2.2.2. CTCF of Third-Order Tensor

Generally, an image is a second-order tensor; then sequential images form a third-order tensor, i.e., a video, adding a temporal dimension on two spatial dimensions. When a new video frame presents and is added to the time dimension of the original tensor, it is defined as a three-dimensional (3D) cumulative tensor. With the number of new frames increasing, the 3D cumulative tensor updates frame by frame.In conventional CP tensor approximation, whenever a new frame of image is added in the time dimension, ALS algorithm needs to be reused to approximate the new cumulative tensor, which is a time consuming process. In addition, the temporal correlation between neighboring frames is not exploited in the decomposition of the cumulative tensor. This paper proposes CTCF to update the CP factorization of original cumulative tensor, obtain the updated factor matrices, and approximate the new frame.Given an original 3D cumulative tensorXori∈ℝI×J×Tori, the result of CP factorization is denoted by Xori≈Aori,Bori,Cori. When a new tensor Xnew∈ℝI×J×Tnew is added in the time dimension, the updated cumulative tensor is X∈ℝI×J×Tori+Tnew, of which the CP factorization appears as X≈A,B,C. We focus on obtaining A, B, and C by updating Aori, Bori, and Cori.The updating process is operated in an alternating way. Firstly, temporal dimensional factor matrixC is computed while factor matrices A and B are fixed; i.e.,(18)C=argminC12X3−CB⊙AT2=argminC12Xori3Xnew3−C1C2B⊙AT2=argminC12Xori3−C1B⊙ATXnew3−C2B⊙AT2,where C is divided into two terms. For A and B are fixed as Aori and Bori, the first row of (18) will be minimized if C1=Cori. To minimize the second row, according to (12), the optimal solution of C2 is Xnew3B⊙AT†, where the symbol “†” denotes Moore–Penrose pseudoinverse of the matrix [51]. So, C can be updated by adding Cnew which is represented by(19)C=CoriCnew=CoriXnew3B⊙AT†.Secondly, factor matrixA is computed while factor matrices B and C are fixed. Similar to 16, the loss function of estimating A is written as(20)ℒ=12X1−AC⊙BT2.Deriveℒ with respect to A; then, we have(21)∂ℒ∂A=X1C⊙B−AC⊙BTC⊙B.To simplifyequation (21), denote P=X1C⊙B and Q=C⊙BTC⊙B; thus, when ∂ℒ/∂A=0, we have A=PQ−1. According to [47], Q can be rewritten as(22)Q=C⊙BTC⊙B=CTC⊗BTB.For computingP, we also divide X1 and C into two terms; i.e.,(23)P=X1C⊙B=Xori1,Xnew1CoriCnew⊙B=Xori1,Xnew1Cori⊙BCnew⊙B=Xori1Cori⊙B+Xnew1Cnew⊙B.SinceB are fixed as Bori, the first term of equation (23) contains only the information of original tensor, which can be expressed by(24)Pori=Xori1Cori⊙B,so,equation (23) is rewritten as(25)P=Pori+Xnew1Cnew⊙B,Hence, P can be updated from Pori using mode-1 unfolding matrix of Xnew and factor matrix mentioned above Cnew. Generally, P is initialized by Xτ∈ℝI×J×τ, which is a small front part of Xori, and updated iteratively by (25). Analogously, the update process of Q can be represented by(26)Q=Qori+Cnew⊙BTCnew⊙B=Qori+CnewTCnew⊗BTB.The update ofA may be summarized as(27)P←P+Xnew1Cnew⊙B,Q←Q+CnewTCnew⊗BTB,A←PQ−1.Finally, the update of factor matrixB may likewise be expressed by(28)U←U+Xnew2Cnew⊙A,V←V+CnewTCnew⊗ATA,B←UV−1,where U=X2C⊙A and V=CTC⊗ATA.To make the process clearer, the proposed CTCF of third-order tensor is summarized by Algorithm1.Algorithm 1: CTCF of third-order tensor.
Input: original 3D cumulative tensor Xori∈ℝI×J×Tori≈Aori,Bori,Cori new tensor Xnew∈ℝI×J×TnewStep 1: new tensor is added in the time dimension andX∈ℝI×J×Tori+Tnew is obtainedStep 2: decomposeX by CP factorization X≈A,B,CStep 3: updateC by (19), with A and B are fixedStep 4: updateA by (27), with B and C are fixedStep 5: updateB by (28), with A and C are fixedStep 6: estimateX by updated A, B and COutput: approximation of updated cumulative tensor X


## 2.2.3. CTCF ofNth-Order Tensor

On the basis ofSection 2.2.2, we try to extend CTCF to higher-order tensors. Suppose an N-dimensional cumulative tensor Xori∈ℝI1×I2×⋯×IN−1×Tori where the last dimension is temporal dimension. The CP factorization of Xori is represented as Xori≈Aori1,⋯,AoriN−1,AoriN. When a new tensor Xori∈ℝI1×I2×⋯×IN−1×Tnew is added in the time dimension, the updated cumulative tensor is X∈ℝI1×I2×⋯×IN−1×Tori+Tnew, of which the CP factorization is denoted by X≈A1,⋯,AN−1,AN.Similar toSection 2.2.2, temporal dimensional factor matrix AN is firstly updated with other N−1 matrices fixed. Like 17, the optimization problem of estimating AN is formulated by(29)AN=argminAN12XN−AN⊙N−1AiT2.We also separate original part from new added part; i.e.,(30)AN=AoriNAnewN=AoriNXnew⊙N−1AiT†,i=1,2,⋯,N−1.The original part is minimized by fixing the firstN−1 factor matrix and the new part is updated by Xnew.The updates of nontemporal dimensional factor matricesAn (n∈1,N−1) may refer to the ones of factor matrices A and B in Section 2.2.2. The loss function ℒn of estimating An is the same as 16. Let ∂ℒn/∂An=0 and introduce matrices Pn and Qn; the update of An may be summarized as(31)Pn←Pn+XnewnAnewN⊙Kn,Qn←Qn+AnewNTAnewN⊗Rn,An←PnQn−1,where Kn=⊙i≠nN−1Ai and Rn=⊗i≠nN−1AiTAi.


## 2.2.4. CTCF-Based Detection Method

In HSV, the sequential data is expressed as a 4D cumulative tensor; the temporal dimension increases with new frames are added in. Whenever a new frame presents, the results of original cumulative tensor CP factorization are updated to obtain the factor matrices of the new cumulative tensor, and the CP tensor approximation of the newly added frame is obtained at the same time. If the target is absent, the CP tensor approximation will lead to a small error, since the background information is similar between adjacent frames. On the contrary, if the error is large, the target is likely to present. We define the fitness between the new frame and its approximation in34. If the fitness is smaller than the threshold, the target is supposed to appear in the new frame. Otherwise, the new frame is added in the temporal dimension and used to update original cumulative tensor.The original 4D cumulative tensor is denoted byXori∈ℝI1×I2×I3×n; n denotes the frame number of initial video. The factor matrices of four dimensions are represented as(32)Xori≈Aori1,Aori2,Aori3,Aori4,where Aori1∈ℝI1×R, Aori2∈ℝI2×R, Aori3∈ℝI3×R, and Aori4∈ℝn×R and R denotes the number of component rank-one tensors in CP factorization. When a new frame Xnew∈ℝI1×I2×I3 is added in the temporal dimension of original 4D cumulative tensor, the 4D cumulative tensor is updated and denoted by X∈ℝI1×I2×I3×n+1. The factor matrices of X are expressed by(33)X≈A1,A2,A3,A4,where A1∈ℝI1×R, A2∈ℝI2×R, A3∈ℝI3×R, and A4∈ℝn+1×R. Based on Section 2.2.3, we estimate A1∼A4 and obtain the approximation of X and X˜new, where X˜new≈A1,A2,A3,A4n+1,:. Actually, it is the specific case when N=4.We define thefitness (Xnew, X˜new) as(34)fitnessXnew,X˜new=1−Xnew−X˜newFXnewF.If the target does not appear, the approximation error is small and the result of fitness is large. Given a preset thresholdη, when fitnessXnew,X˜new>η, i.e., the fitness is larger than η, we decide that the target is absent. Then, the nontarget frame is added in temporal dimension and the updated 4D cumulative tensor becomes the new original 4D cumulative tensor, which can be expressed as(35)Xori⇐X,Aorin⇐An,n∈1,4.If the target appears, the approximation error is large and the fitness is smaller thanη. The residual of Xnew and X˜new is the approximation of the target tensor; i.e.,(36)Xtar≈Xnew−X˜new.The target of each frame will be shown in 2D form by taking the maximum value of every spectrum. In this way, the proposed CTCF-based detection method can extract not only the key-frames where the target presents, but also the approximate region of target in every key-frame. The flowchart of the proposed method is shown in Figure3. In Section 3, experiments on real HSV data are conducted and the proposed method is compared with some representative techniques.


## 2.3. The Proposed STTF-Based Super-Resolution Method

InSection 2.2, we present an approach to detect the frames where the target appears in HSV and the approximate region of the target. However, as discussed in Section 1, there has to be a tradeoff between spectral resolution and the spatial resolution in HSI imaging systems [52]. The spatial resolution is always low since high spectral resolution is required in HSIs and HSV. So, we are interested in improving the spatial resolution of targets after the detecting process. Instead of fusing HR-MSI and LR-HSI, we try to handle the target SR problem by what we have got, which is more practical in real cases.


### 2.3.1. Problem Formulation

In this subsection, HSIs are represented as 3D tensors with three indexes (H,W,S), which stand for the height, width, and spectral modes. X∈ℝH×W×S denotes the HR-HSI and the LR-HSI is denoted by Y∈ℝh×w×S, where W>w and H>h. The goal is to estimate X from Y.There are two significant characteristics of HR-HSIs [53]: the first one is that spectral vectors can be well approximated in low dimensional subspaces, and the second one is that HSIs are spatially self-similar. This means that sparsity exists in both spectral and spatial dimensions. Inspired by sparse representation [54], the low dimensionality in spectral domain gives the possibility to form a spectral mode dictionary S with few nonzero atoms; the self-similarities in spatial domain guarantee the sparse representations of the height and width modes with spatial dictionaries H and W. In this way, the conventional Tucker factorization is transformed into the multiplication of the core tensor and three modes dictionaries. The factorization is illustrated in Figure 4. The HR-HSI is represented as(37)X=Z×1H×2W×3S,where H∈ℝH×zh, W∈ℝW×zw, and S∈ℝS×zs. The variables zh, zw, and zs denote the atoms (i.e., the number of columns) of H, W, and S, respectively. The core tensor Z contains the coefficients of X over three modes dictionaries. We can see that 37 incorporates the information of separated modes into a unified framework.Figure 4
Sparse-based tensor Tucker factorization of an HSI.The LR key-frame of HSVY can be seen as the spatially downsampled version of HR-HSI X, which is written as(38)Y=X×1D1×2D2,where D1∈ℝh×H and D2∈ℝw×W are downsampling matrices of the height and width modes. Substituting 37 into (38), Y is represented by(39)Y=Z×1D1H×2D2W×3S=Z×1H∗×2W∗×3S,where H∗=D1H∈ℝh×zh and W∗=D2W∈ℝw×zw denotes the downsampled dictionary of height and width modes. To recover X, we focus on estimating the dictionaries H, W, and S and the core tensor Z.


### 2.3.2. The Proposed STTF-Based SR Algorithm

SinceY is a downsampled version, recovering X from Y is a typical inverse problem, which is badly ill-posed. So, some prior knowledge of X is needed to regularize the super-resolution problem. In HSI processing, the spectral sparsity is a widespread regularizer applied to solve varieties of ill-posed problems [55–58]. In such regularization, spectral vectors are linearly combined by a small quantity of different spectral signatures. However, these schemes only take advantage of the sparsity of the spectral domain. In the proposed algorithm, taking into account the HSI self-similarity, sparsity regularization is extended to the spatial domain by exploiting the sparse-based tensor Tucker factorization (STTF). In STTF, the HR-HSI performs a united sparse representation of the core tensor and three modes dictionaries.On the basis of equation (39), the HSV frame super-resolution is formulated as a constrained least-squares optimization problem:(40)minH,W,S,ZY−Z×1H∗×2W∗×3SF2,s.t.Z0≤θ,where ⋅F represents the Frobenius norm and θ denotes the number of maximum nonzero elements of Z. Because of the ℓ0-norm constraint, equation (40) is nonconvex. To make the optimization processable, the ℓ0-norm is replaced by the ℓ1-norm and 40 is transformed into an unconstrained version:(41)minH,W,S,ZY−Z×1H∗×2W∗×3SF2+λZ1,where λ is the parameter of sparse regularizer. Equation (41) is also nonconvex, and the solutions of H, W, and S and Z are not unique. Nonetheless, if we focus on only one variable with other variables fixed, the objective function in equation (41) is convex. Inspired by [59, 60], equation (41) can be solved by proximal alternating optimization scheme, which is guaranteed to reach convergence in a particular situation. Concretely, H, W, S, and Z are updated iteratively by(42)H=argminHfH,W,S,Z+αH−HpreF2,W=argminWfH,W,S,Z+αW−WpreF2,S=argminSfH,W,S,Z+αS−SpreF2,Z=argminZfH,W,S,Z+αZ−ZpreF2,where ⋅pre denotes the previous estimation in the last iteration and α denotes a positive number. Equation (41) defines the object function fH,W,S,Z. The optimizations of H, W, S, and Z will be presented detailedly in the appendix. The conjugate gradient (CG) method [61] and the alternating direction method of multipliers (ADMM) [62] will be used in the optimizations.


### 2.3.3. Initialization of the Proposed Method

Since the optimization problem in (41) is nonconvex, the solution would result in poor local minima if we set the initialization carelessly. In this paper, we initialize the spatial dictionaries H∗ and W∗ from Y1 and Y2 dictionary-updates-cycles KSVD (DUC-KSVD) [63]; this method can promote sparse representations. Then, initialization of spectral dictionary S is accomplished by simplex identification split augmented Lagrangian (SISAL) algorithm [64]; this approach can efficiently identify a minimum unit that contains the spectral vectors.The proposed STTF-based SR algorithm is summarized in Algorithm2.Algorithm 2: STTF-based SR method.
Input: LR-HSI YInitializeS with SISALInitializeH∗ and W∗ with DUC-KSVDInitializeZ with (39)while no convergence doStep 1: updateH by solving (A.3) with CGH∗=D1H,Hpre=HStep 2: updateW by solving (A.6) with CGW∗=D2W,Wpre=WStep 3: updateS by solving (A.9) with CGSpre=SStep 4: updateZ by solving (A.15) with CGZpre=Zend whileEstimateX by (37)Output: HR-HSI X


## 2.3.1. Problem Formulation

In this subsection, HSIs are represented as 3D tensors with three indexes (H,W,S), which stand for the height, width, and spectral modes. X∈ℝH×W×S denotes the HR-HSI and the LR-HSI is denoted by Y∈ℝh×w×S, where W>w and H>h. The goal is to estimate X from Y.There are two significant characteristics of HR-HSIs [53]: the first one is that spectral vectors can be well approximated in low dimensional subspaces, and the second one is that HSIs are spatially self-similar. This means that sparsity exists in both spectral and spatial dimensions. Inspired by sparse representation [54], the low dimensionality in spectral domain gives the possibility to form a spectral mode dictionary S with few nonzero atoms; the self-similarities in spatial domain guarantee the sparse representations of the height and width modes with spatial dictionaries H and W. In this way, the conventional Tucker factorization is transformed into the multiplication of the core tensor and three modes dictionaries. The factorization is illustrated in Figure 4. The HR-HSI is represented as(37)X=Z×1H×2W×3S,where H∈ℝH×zh, W∈ℝW×zw, and S∈ℝS×zs. The variables zh, zw, and zs denote the atoms (i.e., the number of columns) of H, W, and S, respectively. The core tensor Z contains the coefficients of X over three modes dictionaries. We can see that 37 incorporates the information of separated modes into a unified framework.Figure 4
Sparse-based tensor Tucker factorization of an HSI.The LR key-frame of HSVY can be seen as the spatially downsampled version of HR-HSI X, which is written as(38)Y=X×1D1×2D2,where D1∈ℝh×H and D2∈ℝw×W are downsampling matrices of the height and width modes. Substituting 37 into (38), Y is represented by(39)Y=Z×1D1H×2D2W×3S=Z×1H∗×2W∗×3S,where H∗=D1H∈ℝh×zh and W∗=D2W∈ℝw×zw denotes the downsampled dictionary of height and width modes. To recover X, we focus on estimating the dictionaries H, W, and S and the core tensor Z.


## 2.3.2. The Proposed STTF-Based SR Algorithm

SinceY is a downsampled version, recovering X from Y is a typical inverse problem, which is badly ill-posed. So, some prior knowledge of X is needed to regularize the super-resolution problem. In HSI processing, the spectral sparsity is a widespread regularizer applied to solve varieties of ill-posed problems [55–58]. In such regularization, spectral vectors are linearly combined by a small quantity of different spectral signatures. However, these schemes only take advantage of the sparsity of the spectral domain. In the proposed algorithm, taking into account the HSI self-similarity, sparsity regularization is extended to the spatial domain by exploiting the sparse-based tensor Tucker factorization (STTF). In STTF, the HR-HSI performs a united sparse representation of the core tensor and three modes dictionaries.On the basis of equation (39), the HSV frame super-resolution is formulated as a constrained least-squares optimization problem:(40)minH,W,S,ZY−Z×1H∗×2W∗×3SF2,s.t.Z0≤θ,where ⋅F represents the Frobenius norm and θ denotes the number of maximum nonzero elements of Z. Because of the ℓ0-norm constraint, equation (40) is nonconvex. To make the optimization processable, the ℓ0-norm is replaced by the ℓ1-norm and 40 is transformed into an unconstrained version:(41)minH,W,S,ZY−Z×1H∗×2W∗×3SF2+λZ1,where λ is the parameter of sparse regularizer. Equation (41) is also nonconvex, and the solutions of H, W, and S and Z are not unique. Nonetheless, if we focus on only one variable with other variables fixed, the objective function in equation (41) is convex. Inspired by [59, 60], equation (41) can be solved by proximal alternating optimization scheme, which is guaranteed to reach convergence in a particular situation. Concretely, H, W, S, and Z are updated iteratively by(42)H=argminHfH,W,S,Z+αH−HpreF2,W=argminWfH,W,S,Z+αW−WpreF2,S=argminSfH,W,S,Z+αS−SpreF2,Z=argminZfH,W,S,Z+αZ−ZpreF2,where ⋅pre denotes the previous estimation in the last iteration and α denotes a positive number. Equation (41) defines the object function fH,W,S,Z. The optimizations of H, W, S, and Z will be presented detailedly in the appendix. The conjugate gradient (CG) method [61] and the alternating direction method of multipliers (ADMM) [62] will be used in the optimizations.


## 2.3.3. Initialization of the Proposed Method

Since the optimization problem in (41) is nonconvex, the solution would result in poor local minima if we set the initialization carelessly. In this paper, we initialize the spatial dictionaries H∗ and W∗ from Y1 and Y2 dictionary-updates-cycles KSVD (DUC-KSVD) [63]; this method can promote sparse representations. Then, initialization of spectral dictionary S is accomplished by simplex identification split augmented Lagrangian (SISAL) algorithm [64]; this approach can efficiently identify a minimum unit that contains the spectral vectors.The proposed STTF-based SR algorithm is summarized in Algorithm2.Algorithm 2: STTF-based SR method.
Input: LR-HSI YInitializeS with SISALInitializeH∗ and W∗ with DUC-KSVDInitializeZ with (39)while no convergence doStep 1: updateH by solving (A.3) with CGH∗=D1H,Hpre=HStep 2: updateW by solving (A.6) with CGW∗=D2W,Wpre=WStep 3: updateS by solving (A.9) with CGSpre=SStep 4: updateZ by solving (A.15) with CGZpre=Zend whileEstimateX by (37)Output: HR-HSI X


## 3. Results and Discussion

### 3.1. Experimental Data Set

To highlight the advantages of HSIs, we choose invisible gas plume to be the target. The proposed algorithms can be extended to other types of data reasonably. In this section, the HSV data set is acquired by the infrared imaging spectrometer “HyperCam-LW.” Sulfur hexafluoride (SF6) is chosen to be the target, since it is a kind of odorless and colorless gas plume with a distinct absorption peak in LWIR range. The HSV data set consists of 60 infrared hyperspectral frames with the size of 128×320×127. The imaging interval is 4.8 s, and the wavelength of the data ranges from 7.8 μm to 11.8 μm.In SR method, only the middle128×128 pixels are used in the experiment (specifically, column 71 to column 198) for reasons connected with the algorithm process. And we remove the spectral band 41–127 because of water vapor absorption and extremely low SNR. At last, the size of input LR-HSI is 128×128×40.


### 3.2. Compared Methods

For CTCF-based detection method, we compare it with two representative methods: MSD (matched subspace detector) [25] and CMF (constrained matrix factorization) [29]. For STTF-based SR method, we compare it with three state-of-the-art algorithms: bicubic interpolation, sparse representation-based SR method [54], and sequence information-based SR method [65].


### 3.3. Qualitative and Quantitative Metrics

For detection methods, receiver operating characteristic (ROC) curves [66] are used to evaluate the performance. Generally, a detector outperforms another one if the area under its ROC curve is larger [67]. As suggested in [68], the area under the ROC curve (AUC) is also calculated as a measure of performance of these detection methods. Usually, a better detector gets a higher AUC value.For SR algorithms, since we directly process the LR-HSI, there is no original HR-HSI (i.e., the ground truth) for reference. Thus, some popular quantitative metrics are not available, such as RMSE (root-mean-square error) [69], PSNR (peak signal to noise ratio), and SAM (spectral angle mapper). In this section, entropy and average gradient are introduced to evaluate the performance of SR methods.


#### 3.3.1. Entropy

Super-resolution aims to introduce more useful information into images, so we may measure the performance of SR methods by calculating the contained information in the experimental results. The entropy is indicated as(43)E=−∑i=1nPilog2Pi.The probability of a pixeli in the image is denoted by Pi and n denotes the grey value range 0∼255. The larger the entropy value of the image, the richer the information contained in the image.


#### 3.3.2. Average Gradient

Another assessment to measure the performance of super-resolution is the change of the amount of detailed information in the image. We may evaluate the experimental results by average gradient, since it can reflect the ability of expressing the details and measuring the clarity of the image. The gradient increases if the greyscale level rate in one direction of the image varies quickly. The average gradient is formulated as(44)G¯=1m×n∑i=1m∑j=1n∂fi,j/∂xi2+∂fi,j/∂yi2/2,where m and n denote the height and width of the image, respectively; fi,j denotes the greyscale value of pixel i,j in the image. The larger the average gradient value of the image is, the clearer the image will be.Besides, the visual quality of output images is an important qualitative metric.


### 3.4. Parameters Setting

In MSD, we pick 463 spectrums of gas target and 846 spectrums of background from the 12th frame of HSV to build up the training set. The size of the target subspace and background space is127×112 and 127×115, respectively. In CMF, the number of endmembers is 3, the sparsity of factor matrices is 2, and number of iteration is 3. In the proposed CTCF-based method, the original cumulative tensor is obtained by ALS, the tensor rank is 3, the maximum iteration number is 100, and the reconstruction error is 10−8; in update stage, the threshold of fitness is 0.9. In the proposed STTF-based SR method, the number of iterations is 5; the parameter α is the weight in (42) and we set α=10−3; parameter λ controls the sparsity of Z; we set λ=10−5; parameter μ is set by μ=10−2; the size of Z is set by zh=240, zw=240, and zs=12. The parameters above are decided after sufficient number of experiments to make a balance between efficiency and stability.


### 3.5. Experimental Results and Discussion

In this subsection, we show the experimental results of the various methods for detection and super-resolution.After processing the HSV by the proposed CTCF-based method, we compute the values of Frobenius norm of each frame, which are presented in Figure5. It is obvious that target gas appears in the 12th frame and disappears in the 51st frame. Figure 6 compares the ROC curves of test methods on four frames in detail, and Figure 7 illustrates the general trends of ROC curves of MSD, CMF, and CTCF, respectively. As can be seen from Figures 6 and 7, the proposed CTCF-based detection algorithm outperforms the other two methods. The AUC values of three approaches are shown in Table 1. In each row, the bold value represents the highest AUC value. Although the AUC values of CMF in some frames are better, we can see that the AUC values of CMF in some other frames are very low (less than 0.98). On the contrast, all the results of CTCF lie in the range of 0.98 to 1. From the average value and the variance (the bold value represents the highest value), we can conclude that the proposed method is superior and more stable. The graphical results are illustrated in Figure 8.Figure 5
The value of Frobenius norm of each frame.Figure 6
The comparison of ROC curves of three detection methods: (a) frame 13; (b) frame 23; (c) frame 34; (d) frame 44.
(a)(b)(c)(d)Figure 7
The ROC curves of 39 key-frames by three detection methods: (a) MSD; (b) CMF; (c) CTCF.
(a)(b)(c)Table 1
Detection quantitative results (AUC value) of the test methods on key-frames.
FrameMSD [25]CMF [29]CTCF120.96550.99930.9980130.84620.99950.9980140.88780.99940.9981150.81890.99650.9981160.87340.99950.9977170.93480.99870.9946180.57920.74770.9915190.78940.99880.9958200.93360.89910.9934210.82220.99800.9966220.90010.99150.9969230.83880.99860.9945240.89140.99900.9983250.91690.99890.9961260.92540.99470.9989270.87220.99740.9951280.85030.99780.9930290.94900.98920.9974300.90110.93410.9885310.91570.98670.9965320.88810.85820.9811330.93450.97710.9933340.90070.99770.9922350.92730.99520.9933360.93490.99890.9950370.95280.99860.9984380.92990.98750.9981390.88380.99790.9962400.92950.99760.9966410.91650.99390.9976420.96650.99860.9988430.96600.97630.9969440.90830.99950.9979450.90460.99770.9935460.91560.99640.9918470.92250.99790.9950480.86230.99800.9967490.88940.99920.9980500.86400.99690.9947Average0.89260.98170.9954Variance0.4407 × 10−20.2290 × 10−20.1134 × 10−6Figure 8
The AUC values 39 key-frames by three detection methods.The target of each key-frame is shown in 2D form (grey image) by taking the maximum value of every spectrum. To save the length of the paper, we choose 8 frames to show the comparison of three detectors, which are shown in Figure9. The first row to the eighth row present the detection result of the chosen frame, of which the frame number is 15, 18, 22, 28, 31, 39, 48, and 50. The higher the greyscale of the pixel in the image is, the closer it is to the target. It is apparent that our method extracts more accurate targets.Figure 9
The 2D form of the detection results by three detection methods.Table2 shows the entropy and average gradient of the key-frames by four SR algorithms. Since sequence-based method needs 5 LR frames to form 1 HR frame, the compared frame number is changed from range 12∼50 to range 14∼48. In each row, the bold values represent the highest entropy value and the highest average gradient value. From Table 2, we can conclude that firstly, although interpolation can add more information in the frame, the details of the target are lost; secondly, sparse representation SR and sequence information SR have almost the same entropy, but the latter approach offers more details because in the method the HR dictionary is formed by several LR dictionaries; finally, the proposed STTF-based SR method outperforms the other three methods in both metrics.Table 2
SR quantitative results (entropy and average gradient) of the test methods on key-frames.
MethodsLR frameBicubic interpolationSparse representation-based SR [54]Sequence information-based SR [65]STTF-based SRFrameEntropyAverage gradientEntropyAverage gradientEntropyAverage gradientEntropyAverage gradientEntropyAverage gradient145.16030.00765.37440.00525.49960.00785.42590.00905.60980.0121154.70860.00775.14070.00635.23980.00895.21840.01035.36780.0139165.55210.00845.80130.00605.87650.00855.82030.01035.97720.0135175.52930.00865.70540.00565.79180.00815.61250.00945.88310.0129184.29890.00724.83390.00624.97940.00885.04230.01065.11080.0138194.48430.00735.00450.00635.13270.00895.18310.01065.26440.0140205.14420.00755.43070.00605.50390.00865.39870.00995.61220.0137214.88210.00715.22340.00605.32640.00865.25780.01005.44910.0137224.34720.00674.89290.00654.98580.00904.94090.01025.12610.0141234.61270.00675.04300.00615.15340.00865.08060.00985.28230.0135244.41890.00644.86880.00604.98500.00854.89160.00965.11680.0133254.32730.00664.90910.00665.00610.00914.96070.01035.14940.0143264.13940.00644.75890.00664.84380.00904.80780.01034.99190.0141274.11270.00654.69950.00664.78780.00914.77130.01044.93660.0142283.96570.00614.65760.00664.74670.00914.75070.01074.89020.0142294.18850.00634.75450.00644.86110.00884.84810.01015.00660.0137303.96720.00634.63450.00674.72530.00924.72410.01064.87520.0144313.94400.00614.61350.00654.71310.00904.73490.01044.86540.0141323.86610.00604.57990.00644.69140.00884.69710.01014.84970.0138334.04790.00604.71000.00644.81260.00884.81140.01004.96310.0137344.16910.00604.78240.00664.86210.00884.84620.01025.00720.0136354.09330.00624.71690.00674.80100.00914.82450.01084.95150.0143363.91570.00634.59950.00674.68810.00924.67120.01034.85080.0142373.78100.00594.50280.00644.60880.00894.59840.01004.76660.0138383.88140.00614.54790.00654.64830.00904.63950.01014.80500.0140394.31680.00604.83970.00614.94060.00844.91350.00995.07920.0130403.93330.00614.65970.00674.73800.00914.72090.01044.89060.0142414.20090.00634.78970.00664.87110.00894.83460.01025.01380.0138424.10830.00634.75140.00674.83620.00914.83980.01034.98360.0142434.04850.00634.68270.00674.76020.00914.71170.01014.91090.0142444.05210.00624.65100.00634.75870.00874.72730.00974.91650.0136454.34420.00604.90110.00615.00790.00854.93800.00975.14130.0134464.00060.00614.59130.00624.70800.00874.64670.00974.85870.0136474.47490.00625.01090.00605.10150.00845.00070.00955.23020.0131484.06850.00654.68430.00654.78590.00904.76070.01014.93840.0141Avg.4.31670.00664.86710.00634.96510.00884.93290.01015.10490.0138Figure10 presents the visual quality of the results obtained by four test methods. We choose the 16th, 21st, 34th, and 47th frames as a representative. The smaller one with size of 128×128 is the LR 2D-form frame. The bigger ones with size of 256×256 are the SR results of different algorithms. As can be seen from Figure 10, the proposed approach yields clearer outputs with sharper edges and more textures. A drawback is the “checkerboard artifacts,” which may be caused by the deconvolution operations in the method. We desired to fix it in our future work.Figure 10
The SR results of four test methods, from left to right: original LR frame, bicubic interpolation, sparse representation-based SR method [54], sequence information-based SR method [65], the proposed STTF-based SR method: (a) frame 16; (b) frame 21; (c) frame 34; (d) frame 47.
(a)(b)(c)(d)


## 3.1. Experimental Data Set

To highlight the advantages of HSIs, we choose invisible gas plume to be the target. The proposed algorithms can be extended to other types of data reasonably. In this section, the HSV data set is acquired by the infrared imaging spectrometer “HyperCam-LW.” Sulfur hexafluoride (SF6) is chosen to be the target, since it is a kind of odorless and colorless gas plume with a distinct absorption peak in LWIR range. The HSV data set consists of 60 infrared hyperspectral frames with the size of 128×320×127. The imaging interval is 4.8 s, and the wavelength of the data ranges from 7.8 μm to 11.8 μm.In SR method, only the middle128×128 pixels are used in the experiment (specifically, column 71 to column 198) for reasons connected with the algorithm process. And we remove the spectral band 41–127 because of water vapor absorption and extremely low SNR. At last, the size of input LR-HSI is 128×128×40.


## 3.2. Compared Methods

For CTCF-based detection method, we compare it with two representative methods: MSD (matched subspace detector) [25] and CMF (constrained matrix factorization) [29]. For STTF-based SR method, we compare it with three state-of-the-art algorithms: bicubic interpolation, sparse representation-based SR method [54], and sequence information-based SR method [65].


## 3.3. Qualitative and Quantitative Metrics

For detection methods, receiver operating characteristic (ROC) curves [66] are used to evaluate the performance. Generally, a detector outperforms another one if the area under its ROC curve is larger [67]. As suggested in [68], the area under the ROC curve (AUC) is also calculated as a measure of performance of these detection methods. Usually, a better detector gets a higher AUC value.For SR algorithms, since we directly process the LR-HSI, there is no original HR-HSI (i.e., the ground truth) for reference. Thus, some popular quantitative metrics are not available, such as RMSE (root-mean-square error) [69], PSNR (peak signal to noise ratio), and SAM (spectral angle mapper). In this section, entropy and average gradient are introduced to evaluate the performance of SR methods.


### 3.3.1. Entropy

Super-resolution aims to introduce more useful information into images, so we may measure the performance of SR methods by calculating the contained information in the experimental results. The entropy is indicated as(43)E=−∑i=1nPilog2Pi.The probability of a pixeli in the image is denoted by Pi and n denotes the grey value range 0∼255. The larger the entropy value of the image, the richer the information contained in the image.


### 3.3.2. Average Gradient

Another assessment to measure the performance of super-resolution is the change of the amount of detailed information in the image. We may evaluate the experimental results by average gradient, since it can reflect the ability of expressing the details and measuring the clarity of the image. The gradient increases if the greyscale level rate in one direction of the image varies quickly. The average gradient is formulated as(44)G¯=1m×n∑i=1m∑j=1n∂fi,j/∂xi2+∂fi,j/∂yi2/2,where m and n denote the height and width of the image, respectively; fi,j denotes the greyscale value of pixel i,j in the image. The larger the average gradient value of the image is, the clearer the image will be.Besides, the visual quality of output images is an important qualitative metric.


## 3.3.1. Entropy

Super-resolution aims to introduce more useful information into images, so we may measure the performance of SR methods by calculating the contained information in the experimental results. The entropy is indicated as(43)E=−∑i=1nPilog2Pi.The probability of a pixeli in the image is denoted by Pi and n denotes the grey value range 0∼255. The larger the entropy value of the image, the richer the information contained in the image.


## 3.3.2. Average Gradient

Another assessment to measure the performance of super-resolution is the change of the amount of detailed information in the image. We may evaluate the experimental results by average gradient, since it can reflect the ability of expressing the details and measuring the clarity of the image. The gradient increases if the greyscale level rate in one direction of the image varies quickly. The average gradient is formulated as(44)G¯=1m×n∑i=1m∑j=1n∂fi,j/∂xi2+∂fi,j/∂yi2/2,where m and n denote the height and width of the image, respectively; fi,j denotes the greyscale value of pixel i,j in the image. The larger the average gradient value of the image is, the clearer the image will be.Besides, the visual quality of output images is an important qualitative metric.


## 3.4. Parameters Setting

In MSD, we pick 463 spectrums of gas target and 846 spectrums of background from the 12th frame of HSV to build up the training set. The size of the target subspace and background space is127×112 and 127×115, respectively. In CMF, the number of endmembers is 3, the sparsity of factor matrices is 2, and number of iteration is 3. In the proposed CTCF-based method, the original cumulative tensor is obtained by ALS, the tensor rank is 3, the maximum iteration number is 100, and the reconstruction error is 10−8; in update stage, the threshold of fitness is 0.9. In the proposed STTF-based SR method, the number of iterations is 5; the parameter α is the weight in (42) and we set α=10−3; parameter λ controls the sparsity of Z; we set λ=10−5; parameter μ is set by μ=10−2; the size of Z is set by zh=240, zw=240, and zs=12. The parameters above are decided after sufficient number of experiments to make a balance between efficiency and stability.


## 3.5. Experimental Results and Discussion

In this subsection, we show the experimental results of the various methods for detection and super-resolution.After processing the HSV by the proposed CTCF-based method, we compute the values of Frobenius norm of each frame, which are presented in Figure5. It is obvious that target gas appears in the 12th frame and disappears in the 51st frame. Figure 6 compares the ROC curves of test methods on four frames in detail, and Figure 7 illustrates the general trends of ROC curves of MSD, CMF, and CTCF, respectively. As can be seen from Figures 6 and 7, the proposed CTCF-based detection algorithm outperforms the other two methods. The AUC values of three approaches are shown in Table 1. In each row, the bold value represents the highest AUC value. Although the AUC values of CMF in some frames are better, we can see that the AUC values of CMF in some other frames are very low (less than 0.98). On the contrast, all the results of CTCF lie in the range of 0.98 to 1. From the average value and the variance (the bold value represents the highest value), we can conclude that the proposed method is superior and more stable. The graphical results are illustrated in Figure 8.Figure 5
The value of Frobenius norm of each frame.Figure 6
The comparison of ROC curves of three detection methods: (a) frame 13; (b) frame 23; (c) frame 34; (d) frame 44.
(a)(b)(c)(d)Figure 7
The ROC curves of 39 key-frames by three detection methods: (a) MSD; (b) CMF; (c) CTCF.
(a)(b)(c)Table 1
Detection quantitative results (AUC value) of the test methods on key-frames.
FrameMSD [25]CMF [29]CTCF120.96550.99930.9980130.84620.99950.9980140.88780.99940.9981150.81890.99650.9981160.87340.99950.9977170.93480.99870.9946180.57920.74770.9915190.78940.99880.9958200.93360.89910.9934210.82220.99800.9966220.90010.99150.9969230.83880.99860.9945240.89140.99900.9983250.91690.99890.9961260.92540.99470.9989270.87220.99740.9951280.85030.99780.9930290.94900.98920.9974300.90110.93410.9885310.91570.98670.9965320.88810.85820.9811330.93450.97710.9933340.90070.99770.9922350.92730.99520.9933360.93490.99890.9950370.95280.99860.9984380.92990.98750.9981390.88380.99790.9962400.92950.99760.9966410.91650.99390.9976420.96650.99860.9988430.96600.97630.9969440.90830.99950.9979450.90460.99770.9935460.91560.99640.9918470.92250.99790.9950480.86230.99800.9967490.88940.99920.9980500.86400.99690.9947Average0.89260.98170.9954Variance0.4407 × 10−20.2290 × 10−20.1134 × 10−6Figure 8
The AUC values 39 key-frames by three detection methods.The target of each key-frame is shown in 2D form (grey image) by taking the maximum value of every spectrum. To save the length of the paper, we choose 8 frames to show the comparison of three detectors, which are shown in Figure9. The first row to the eighth row present the detection result of the chosen frame, of which the frame number is 15, 18, 22, 28, 31, 39, 48, and 50. The higher the greyscale of the pixel in the image is, the closer it is to the target. It is apparent that our method extracts more accurate targets.Figure 9
The 2D form of the detection results by three detection methods.Table2 shows the entropy and average gradient of the key-frames by four SR algorithms. Since sequence-based method needs 5 LR frames to form 1 HR frame, the compared frame number is changed from range 12∼50 to range 14∼48. In each row, the bold values represent the highest entropy value and the highest average gradient value. From Table 2, we can conclude that firstly, although interpolation can add more information in the frame, the details of the target are lost; secondly, sparse representation SR and sequence information SR have almost the same entropy, but the latter approach offers more details because in the method the HR dictionary is formed by several LR dictionaries; finally, the proposed STTF-based SR method outperforms the other three methods in both metrics.Table 2
SR quantitative results (entropy and average gradient) of the test methods on key-frames.
MethodsLR frameBicubic interpolationSparse representation-based SR [54]Sequence information-based SR [65]STTF-based SRFrameEntropyAverage gradientEntropyAverage gradientEntropyAverage gradientEntropyAverage gradientEntropyAverage gradient145.16030.00765.37440.00525.49960.00785.42590.00905.60980.0121154.70860.00775.14070.00635.23980.00895.21840.01035.36780.0139165.55210.00845.80130.00605.87650.00855.82030.01035.97720.0135175.52930.00865.70540.00565.79180.00815.61250.00945.88310.0129184.29890.00724.83390.00624.97940.00885.04230.01065.11080.0138194.48430.00735.00450.00635.13270.00895.18310.01065.26440.0140205.14420.00755.43070.00605.50390.00865.39870.00995.61220.0137214.88210.00715.22340.00605.32640.00865.25780.01005.44910.0137224.34720.00674.89290.00654.98580.00904.94090.01025.12610.0141234.61270.00675.04300.00615.15340.00865.08060.00985.28230.0135244.41890.00644.86880.00604.98500.00854.89160.00965.11680.0133254.32730.00664.90910.00665.00610.00914.96070.01035.14940.0143264.13940.00644.75890.00664.84380.00904.80780.01034.99190.0141274.11270.00654.69950.00664.78780.00914.77130.01044.93660.0142283.96570.00614.65760.00664.74670.00914.75070.01074.89020.0142294.18850.00634.75450.00644.86110.00884.84810.01015.00660.0137303.96720.00634.63450.00674.72530.00924.72410.01064.87520.0144313.94400.00614.61350.00654.71310.00904.73490.01044.86540.0141323.86610.00604.57990.00644.69140.00884.69710.01014.84970.0138334.04790.00604.71000.00644.81260.00884.81140.01004.96310.0137344.16910.00604.78240.00664.86210.00884.84620.01025.00720.0136354.09330.00624.71690.00674.80100.00914.82450.01084.95150.0143363.91570.00634.59950.00674.68810.00924.67120.01034.85080.0142373.78100.00594.50280.00644.60880.00894.59840.01004.76660.0138383.88140.00614.54790.00654.64830.00904.63950.01014.80500.0140394.31680.00604.83970.00614.94060.00844.91350.00995.07920.0130403.93330.00614.65970.00674.73800.00914.72090.01044.89060.0142414.20090.00634.78970.00664.87110.00894.83460.01025.01380.0138424.10830.00634.75140.00674.83620.00914.83980.01034.98360.0142434.04850.00634.68270.00674.76020.00914.71170.01014.91090.0142444.05210.00624.65100.00634.75870.00874.72730.00974.91650.0136454.34420.00604.90110.00615.00790.00854.93800.00975.14130.0134464.00060.00614.59130.00624.70800.00874.64670.00974.85870.0136474.47490.00625.01090.00605.10150.00845.00070.00955.23020.0131484.06850.00654.68430.00654.78590.00904.76070.01014.93840.0141Avg.4.31670.00664.86710.00634.96510.00884.93290.01015.10490.0138Figure10 presents the visual quality of the results obtained by four test methods. We choose the 16th, 21st, 34th, and 47th frames as a representative. The smaller one with size of 128×128 is the LR 2D-form frame. The bigger ones with size of 256×256 are the SR results of different algorithms. As can be seen from Figure 10, the proposed approach yields clearer outputs with sharper edges and more textures. A drawback is the “checkerboard artifacts,” which may be caused by the deconvolution operations in the method. We desired to fix it in our future work.Figure 10
The SR results of four test methods, from left to right: original LR frame, bicubic interpolation, sparse representation-based SR method [54], sequence information-based SR method [65], the proposed STTF-based SR method: (a) frame 16; (b) frame 21; (c) frame 34; (d) frame 47.
(a)(b)(c)(d)


## 4. Conclusions

In this paper, aiming at hyperspectral video, we propose a novel key-frame and target detection method based on cumulative tensor CP factorization, termed as CTCF, and a super-resolution algorithm based on sparse-based tensor Tucker factorization, called STTF. Unlike conventional matrix factorization based methods, CTCF considers hyperspectral video (HSV) as 4D cumulative tensor and approximates new added frames by updating factor matrices. To break the limit of conventional methods and make super-resolution (SR) more practical, STTF exploits the sparsity of HSV frames and factorizes them as a sparse core tensor multiplied by three modes dictionaries. In this way, spatial resolution of LR-HSI is enhanced directly without HR samples. The experimental results systematically prove that the proposed CTCF and STTF methods outperform other state-of-the-art algorithms.In the future works, we focus on tensor factorization based target tracking methods which are able to extract target region more accurately and clearly. For super-resolution, we aim at exploiting nonlocal similarities in tensor factorization framework, which has been widely used in inverse problems. Besides target tracking and super-resolution, regions of interest (ROI) approaches will be investigated, in order to make HSV target recognition more efficient and full featured. Inspired by [70] and other related works, we believe that the researches of chemical gas detecting methods will benefit the agricultural application of HSI/HSV. These studies will be of great significance in internet of things (IoT), smart agriculture, pollution monitoring, etc.



---

*Source: 9548749-2020-07-14.xml*