A Novel Hyperspectral Image Simulation Method Based on Nonnegative Matrix Factorization

: Hyperspectral (HS) images can provide abundant and ﬁne spectral information on land surface. However, their applications may be limited by their narrow bandwidth and small coverage area. In this paper, we propose an HS image simulation method based on nonnegative matrix factorization (NMF), which aims at generating HS images using existing multispectral (MS) data. Our main novelty is proposing a spectral transformation matrix and new simulation method. First, we develop a spectral transformation matrix that transforms HS endmembers into MS endmembers. Second, we utilize an iteration scheme to optimize the HS and MS endmembers. The test MS image is then factorized by the MS endmembers to obtain the abundance matrix. The result image is constructed by multiplying the abundance matrix by the HS endmembers. Experiments prove that our method provides high spectral quality by combining prior spectral endmembers. The iteration schemes reduce the simulation error and improve the accuracy of the results. In comparative trials, the spectral angle, RMSE, and correlation coe ﬃ cient of our method are 5.986, 284.6, and 0.905, respectively. Thus, our method outperforms other simulation methods.


Introduction
Hyperspectral (HS) remote sensing technology can acquire a wide-wavelength-range and fine spectral information of an observation area through numerous spectral channels. The abundant spectral information enables widespread applications, such as soil contamination [1], geological mineral mapping [2], fire detection [3], and vegetation monitoring [4]. However, the coverage area and observation time restrict wider applications of HS images [5]. For example, the well-known HS sensor Hyperion has a swath width of 7.5 km, which is considerably narrower than that of other satellite images [6]. The observed objects in a scene are significantly decreased. Meanwhile, the visit cycle of Hyperion is 200 days, which limits the application of Hyperion data [7]. For cloudy areas, the data collection will be further difficult because of the long revisiting period [8,9]. Therefore, HS images are still difficult to obtain from large-area and more frequent coverages. In comparison with HS images, multispectral (MS) data have the advantages of economic accessibility, frequent observation, and global coverage [10]. Furthermore, MS images can be provided by numerous satellites, such as Sentinel-2 [11], Proba-V [12], ALOS [13], and Landsat series [14], which have accumulated massive and long-term historical record of land surface data. Hence, a simulation method must be designed to generate HS images on the basis of the existing MS data of the study area. With these simulated HS images, the application fields of the HS technology can be extended in land cover mapping [15], soil monitoring [16], and agriculture management [17]. which can transform the HS endmembers into MS ones. Then, we develop the new simulation method based on the coupled NMF [32] to factorize the HS and MS images. The final simulated HS image is produced by combining the abundance matrix of the MS image and the HS endmember matrix. This paper is organized as follows. Section 2 initially introduces the experimental datasets and then discusses the novel spectral transformation matrix and proposed HS image simulation method. Section 3 presents the experiments and results on four datasets and reports the comparisons between the proposed and other methods. Section 4 provides additional discussions about our approach and a detailed parameter analysis. Finally, Section 5 draws the conclusions.

Datasets
In this paper, four datasets from Earth Observation-1 (EO-1) and Huangjing-1A (HJ-1A) satellites are selected as experimental data. Each dataset contains HS and MS images of the same area. The first dataset covers Hongkong City, with dense buildings, abundant vegetation, and sea. The second dataset covers Longyan City in Fujian Province, with urban, mining, and vegetated areas. The third dataset presents a suburban area of Wuhan in the middle of China, with agriculture regions, and a part of the Yangtze River. The fourth dataset locates near Baotou, with agriculture and mountain areas in the north of China. The training and test images of Dataset 4 are not from the same scene, which could provide cross validation of the simulation methods. Figure 1 shows the geographic locations of the images in Dataset 4. Table 1 presents further detailed information about the test datasets. The test datasets can objectively validate the performance of the proposed method in different environments. the proposed and other methods. Section 4 provides additional discussions about our approach and a detailed parameter analysis. Finally, Section 5 draws the conclusions.

Datasets
In this paper, four datasets from Earth Observation-1 (EO-1) and Huangjing-1A (HJ-1A) satellites are selected as experimental data. Each dataset contains HS and MS images of the same area. The first dataset covers Hongkong City, with dense buildings, abundant vegetation, and sea. The second dataset covers Longyan City in Fujian Province, with urban, mining, and vegetated areas. The third dataset presents a suburban area of Wuhan in the middle of China, with agriculture regions, and a part of the Yangtze River. The fourth dataset locates near Baotou, with agriculture and mountain areas in the north of China. The training and test images of Dataset 4 are not from the same scene, which could provide cross validation of the simulation methods. Figure 1 shows the geographic locations of the images in Dataset 4. Table 1 presents further detailed information about the test datasets. The test datasets can objectively validate the performance of the proposed method in different environments.    In Datasets 1, 3 and 4, the Hyperion data provide 242 bands with the wavelength ranges from 0.356 µm to 2.577 µm [7]. However, some bands produce a considerable amount of noise because of the malfunction in the Hyperion sensor. We manually select 154 bands (10th−57th, 80th−120th, 134th−164th, 187th−220th bands) with good quality to compose the test HS image [33,34]. The ALI data covers the wavelength range from 0.433 µm to 2.35 µm by 9 bands. All the bands of ALI data are used in our experiments. In Dataset 2, the original HSI data contains 115 bands. The wavelength ranges from 0.45 µm to 0.95 µm. We also selected 92 bands (21st−112th bands) with little noise as our experiment HS data. The CCD data covers a 0.43-0.9 µm wavelength range. All the 4 bands in CCD data are contained in the test MS image. Before experiments, all the test images are preprocessed by the commercial software ENVI 5.3. The preprocessing steps include radiation calibration, atmospheric correction and image registration. More detailed information about the spectral wavelengths of the sensors are presented in Figure 2. In Datasets 1, 3 and 4, the Hyperion data provide 242 bands with the wavelength ranges from 0.356 μm to 2.577 μm [7]. However, some bands produce a considerable amount of noise because of the malfunction in the Hyperion sensor. We manually select 154 bands (10 th -57 th , 80 th -120 th , 134 th -164 th , 187 th -220 th bands) with good quality to compose the test HS image [33,34]. The ALI data covers the wavelength range from 0.433 μm to 2.35 μm by 9 bands. All the bands of ALI data are used in our experiments. In Dataset 2, the original HSI data contains 115 bands. The wavelength ranges from 0.45 μm to 0.95 μm. We also selected 92 bands (21 st -112 th bands) with little noise as our experiment HS data. The CCD data covers a 0.43-0.9 μm wavelength range. All the 4 bands in CCD data are contained in the test MS image. Before experiments, all the test images are preprocessed by the commercial software ENVI 5.3. The preprocessing steps include radiation calibration, atmospheric correction and image registration. More detailed information about the spectral wavelengths of the sensors are presented in Figure 2.  Figure 3 presents an overview of our proposed method. We prepare MS and overlapped HS images. A small subarea in the MS image and the corresponding area in the HS image are extracted as the training MS and HS images, respectively. The rest of the MS image is the test MS image, which is used to simulate the HS image. In our method, the simulated images are produced by combining an HS endmember matrix and an abundance matrix. First, we extract HS endmembers from the training HS image and obtain corresponding MS endmembers through spectral transformation matrix. The spectral transformation matrix is estimated using the proposed method in Section 2.2.2. Second, we use the iteration computations of coupled NMF to update the MS and HS endmembers alternately. The iteration scheme improves the accuracy of the HS and MS endmembers in the factorization of the HS and MS images. Afterward, we utilize the updated MS endmembers to factorize the test MS image through a new round of iterations based on multiplicative update rules. The abundance matrix of the test MS image is then achieved. The simulated HS image is produced by multiplying the achieved abundance matrix by the updated HS endmember matrix.

Methods
In this part, Section 2.2.1 reviews the nonnegative matrix factorization of images. The estimation of the spectral transformation matrix is presented in Section 2.2.2. The iteration schemes of image simulation are introduced in detail in Section 2.2.3.  Figure 3 presents an overview of our proposed method. We prepare MS and overlapped HS images. A small subarea in the MS image and the corresponding area in the HS image are extracted as the training MS and HS images, respectively. The rest of the MS image is the test MS image, which is used to simulate the HS image. In our method, the simulated images are produced by combining an HS endmember matrix and an abundance matrix. First, we extract HS endmembers from the training HS image and obtain corresponding MS endmembers through spectral transformation matrix. The spectral transformation matrix is estimated using the proposed method in Section 2.2.2. Second, we use the iteration computations of coupled NMF to update the MS and HS endmembers alternately. The iteration scheme improves the accuracy of the HS and MS endmembers in the factorization of the HS and MS images. Afterward, we utilize the updated MS endmembers to factorize the test MS image through a new round of iterations based on multiplicative update rules. The abundance matrix of the test MS image is then achieved. The simulated HS image is produced by multiplying the achieved abundance matrix by the updated HS endmember matrix.

Methods
In this part, Section 2.2.1 reviews the nonnegative matrix factorization of images. The estimation of the spectral transformation matrix is presented in Section 2.2.2. The iteration schemes of image simulation are introduced in detail in Section 2.2.3.

Nonnegative Matrix Factorization of Images
Each pixel in the remote sensing image covers an area on the ground. The covered area often contains multiple classes of materials due to the complex distribution of ground objects. Each material has a typical spectrum. Therefore, the spectral information in a pixel can be divided into several types of basic pure spectrums. The linear spectral unmixing model has been widely used to represent the relationship between pixel and basic spectrums [35]. This model assumes that each pixel can be approximated to a weighted sum of several pure spectrums of different materials. In many previous studies, the pure spectrums and their corresponding weight coefficients are known as spectral endmembers and abundance [35,36]. The mathematical expression of the linear spectral unmixing model is as follows: where is the pixel in the remote sensing image, is the spectral endmembers, is the number of endmembers, represents the weight coefficient of each endmember, and is the residual. The linear spectral unmixing model has advantages of intuitive interpretation, sufficient accuracy, and simplicity [37].
By using the linear spectral unmixing model, each pixel in images can be decomposed into a group of endmembers and their abundance. Equation (1) is rewritten in matrix form as where = , … , ∈ ℝ × is a pixel with bands. ∈ ℝ × is the endmember matrix with endmembers, in which each column is an endmember spectrum. = , … , ∈ ℝ × is the abundance matrix, where denotes the abundance value of the ith endmember. = , … , ∈ ℝ × is the residual matrix. Generally, the abundance is interpreted as the area proportions of materials in a pixel.
The physical meaning indicates that the abundance matrix should be constrained by two conditions: 1) all elements in matrix are nonnegative and 2) the sum of any row in is one. Remote sensing image is composed of numerous pixels. To factorize an image, we initially arrange the image in a new matrix with the size of λ × .
denotes the number of pixels in the image. Subsequently, the image is transformed from a 3D matrix to a 2D one. Each column represents a pixel

Nonnegative Matrix Factorization of Images
Each pixel in the remote sensing image covers an area on the ground. The covered area often contains multiple classes of materials due to the complex distribution of ground objects. Each material has a typical spectrum. Therefore, the spectral information in a pixel can be divided into several types of basic pure spectrums. The linear spectral unmixing model has been widely used to represent the relationship between pixel and basic spectrums [35]. This model assumes that each pixel can be approximated to a weighted sum of several pure spectrums of different materials. In many previous studies, the pure spectrums and their corresponding weight coefficients are known as spectral endmembers and abundance [35,36]. The mathematical expression of the linear spectral unmixing model is as follows: p = c 1 e 1 + c 2 e 2 + · · · + c n e n + r where p is the pixel in the remote sensing image, e i is the spectral endmembers, n is the number of endmembers, c i represents the weight coefficient of each endmember, and r is the residual.
The linear spectral unmixing model has advantages of intuitive interpretation, sufficient accuracy, and simplicity [37]. By using the linear spectral unmixing model, each pixel in images can be decomposed into a group of endmembers and their abundance. Equation (1) is rewritten in matrix form as where P = (p 1 , . . . , p λ ) T ∈ R λ×1 is a pixel with λ bands. E ∈ R λ×n is the endmember matrix with n endmembers, in which each column is an endmember spectrum. C = (c 1 , . . . , c n ) T ∈ R n×1 is the abundance matrix, where c i denotes the abundance value of the ith endmember. N = (r 1 , . . . , r λ ) T ∈ R λ×1 is the residual matrix. Generally, the abundance is interpreted as the area proportions of materials in a pixel. The physical meaning indicates that the abundance matrix should be constrained by two conditions: (1) all elements in matrix C are nonnegative and (2) the sum of any row in C is one.
Remote sensing image is composed of numerous pixels. To factorize an image, we initially arrange the image in a new matrix with the size of λ × L. L denotes the number of pixels in the image. Subsequently, the image is transformed from a 3D matrix to a 2D one. Each column represents a pixel with λ channels. The abundance of the pixel is a column vector. To represent the abundance of all pixels, we increase the number of columns in the abundance matrix. Then, the image matrix can be reconstructed by multiplying an endmember matrix with a 2D abundance matrix where u ∈ R λ×L denotes the rearranged image, E ∈ R λ×n is the endmember matrix, and C ∈ R n×L represents the abundance matrix of the image. The residual N ∈ R λ×L is ignored.

Estimation of Spectral Transformation Matrix
In the HS image simulation problem, the simulated HS image is constructed on the basis of a test MS image of the same area. The distributions of ground objects are the same in the HS and MS images. Therefore, each pixel in HS image should contain similar materials and proportions to the corresponding pixel in the MS image. If endmembers of the same group of materials are used to unmix the HS and MS images, then the abundance matrix of the HS image should be equal to that of the MS image. We assume that the simulated HS image can be obtained by multiplying the HS endmember matrix by the MS abundance matrix. On the basis of this hypothesis, we must build the transformation relation between the endmember matrixes of the HS and MS images.
HS images can provide wide and continuous ranges of spectral information. As one of the well-known HS sensors, Hyperion covers a 356-2577 nm wavelength range [38]. Meanwhile, the high spectral resolution (10 nm) makes Hyperion image contain 242 bands [38]. The wavelength ranges of MS images are more discrete and considerably narrower than those of HS images. For example, the MS sensor Advanced Land Imagery (ALI) provides discrete wavelength coverage in nine bands. Figure 2a presents an intuitive comparison between wavelength ranges of Hyperion and ALI data.
As presented in Figure 2, the wavelength range of the bands in MS image is divided into many small spectral ranges in HS image. The HS bands can be seen as fine divisions of the MS channels. We assume that an MS band can be obtained by linearly combining the HS bands in the same wavelength range with the target MS band. For example, the 8th band of the ALI image can be reconstructed using the 141th-160th bands of the Hyperion image ( Figure 2a). We do not utilize the other bands in HS image, outside the range of 141th-160th, in this reconstruction. The reason is that the wavelength range of the other bands is excluded in the range of the target MS band. The spectral features acquired by the other HS bands may not be contained in this target MS band. Therefore, only HS bands in the same wavelength range are selected to produce the corresponding MS channel. The assumption is formulized as where M i is the produced MS band, H j . . . H j+t are the HS channels in the given wavelength ranges, w i1 . . . w it denote the weights of each channel, and r i is the residual. Then, we need to estimate the weights of the HS channels. The multiple linear regression model is selected in our approach. We initially use the training HS and MS images to develop the multiple linear regression model. The relation function is constructed as where u mi ∈ R 1×L is the training MS channel rearranged in a 1 × L matrix, u hj ∈ R t×L is the training HS bands rearranged in a t × L matrix, W i ∈ R 1×t is the weight matrix that need to be estimated, and N i is the residual. The weight matrix can be estimated using the linear regression model: Following this strategy, the weight matrix of each MS band can be obtained. We integrate the weight matrixes of all MS channels as an overall matrix.
where W ∈ R m×h and m is the number of bands in the MS image. In our method, W is the spectral transformation matrix, which can realize the spectral degradation of the global HS image. Finally, we can use the spectral transformation matrix to build the image relation function: where u m ∈ R m×L and u h ∈ R h×L are the rearranged training MS and HS images, respectively. h denotes the number of channels in the HS image, and the residual is ignored.
In addition, it is worth noting that not all the MS channels are included in the wavelength range of the simulated HS image. For instance, the HS data acquired by Huanjing-1A satellite (HJ-1A) covers the wavelength from 450 nm to 950 nm [39]. When we use Landsat MS images to simulate the HJ-1A HS data, the shortwave infrared bands of Landsat are out of the wavelength range. These MS channels are inappropriate to be reconstructed by HS bands with different wavelengths. Thus, our method disregards such MS channels in simulating the HJ-1A HS data. If the Landsat MS image is used to simulate the Hyperion data, then the shortwave infrared bands are applicable. In our method, the spectral relations are primarily built between the MS and HS bands with the same wavelength.

HS Image Simulation and Iterative Calculation Scheme
On the basis of the NMF of images, our approach simulates HS images by multiplying an HS endmember matrix and the abundance matrix of the test MS image.
where H ∈ R h×L is the simulated HS image, E H ∈ R h×n is the HS endmember matrix, and C M ∈ R n×L is the abundance matrix. To achieve an ideal result, we need to minimize the residual error N. We build a cost function F . E H and C M , which can minimize the cost function, are our ideal results. However, the simulated image H is unknown and unavailable in solving this function. Then, we use training HS and MS images to compute the accurate E H and C M .
In our approach, the training HS image is already prepared in the spectral transformation matrix estimation. We assume that the endmembers, which can properly unmix the training HS data, are also suitable for the global simulated HS image. Therefore, the required endmembers can minimize the residuals in the factorizations of the training HS and MS images. We express these residual minimization as cost functions , where u h ∈ R h×L and u m ∈ R m×L are the overlapped training HS and MS images. As u h and u m are factorized by endmembers from the same materials, their abundance matrix should be equal. Meanwhile, the HS and MS endmembers are from the same group of materials. The MS endmembers can be considered as degraded HS endmembers where W is the spectral transformation matrix estimated in Section 2.2.2. To obtain the endmembers that minimize the residuals, we use multiplicative update rules in our simulation approach. The multiplicative update rules, which were proposed by Lee et al., can converge to the local optimal solution under the nonnegativity constraints of two factorized matrixes [29,40]. Yokoya et al. developed a coupled NMF method, which achieves the accurate endmember and abundance matrixes of the test images, on the basis of the multiplicative update rules [32]. Initially, we use vertex component analysis to extract HS endmembers from the training HS image [41]. The number of endmembers is manually set as n. The extracted HS endmembers are arranged as the endmember matrix E h of the training HS image. The abundance matrix C h of the training HS image is initialized as 1/n. Subsequently, we use two rounds of iterative computations to optimize the matrixes of training HS and MS images. In the first round, we fix E h . C h is updated using the iterative computation in Equation (11) until the convergence. In the implementation, to avoid excessive computation, the iterative times of Equation (11) is set to no larger than I 1 , namely, the inner iteration times. The updated C h and initial E h are then alternately optimized by Equations (11) and (12). The iterative calculation is terminated when E h and C h stabilize. In the implementation, the iterative times of Equations (11) and (12) are also limited to no larger than I 1 . In the second round, we use iterative computations to achieve E m and C m . The MS endmember matrix E m is initialized by degrading E h using Equation (10). The initial abundance matrix C m of the training MS image is equal to C h . C m is updated by Equation (13) with fixed E m . The updated C m and initial E m are then alternately optimized by Equations (13) and (14). Both iterative computations are ended with the convergence. In the implementation, the maximum iteration times in these two computations are set to I 1 . Afterward, we repeat the above two rounds of computations as global optimization to obtain the optimal E h and E m . In the first round of repetitions, initial C h is set equal to the updated C m of Equation (13) in the last iteration. E h is updated by Equation (11) with fixed C h . The following alternate updates of E h and C h are the same as above. The maximum repeat times of the two rounds are set as I 2 , namely, the outer iteration times.

Results
In this section, we initially introduce our evaluation protocol of the image simulation. Then, the experimental results and discussions are provided. We compare the spatial and spectral qualities of our method, the PHITA [5], Liu's method [19], and Chen's method [18].

Evaluation and Comparative Methods
The quality assessment of the simulated HS images is important in the experiments. To use a real HS image as reference, we evaluate the simulated results by referring to the evaluation protocol of the image pansharpening [44]. Figure 4 shows the evaluation flow. Initially, we prepare MS and HS images in the same location. The MS image is resampled to the spatial resolution of the HS image. Then, we extract the overlapped areas in the two images as experiment data. To train the proposed method, we select a subarea in the MS image and the corresponding part in the HS image as the training MS and HS images, respectively. The remaining MS image is used as the test MS image to simulate the HS image. For methods without training data, the HS image is directly produced using the test MS image. Afterward, the remaining HS image is regarded as the reference with ideal quality. The HS images simulated by different approaches are compared with the reference image to assess their qualities objectively.  The quality assessment of the simulated HS images primarily focuses on the visual appearance and quantitative evaluation. The visual appearance provides a general idea of the results. The quantitative evaluation provides fine and objective comparisons. The evaluation indices in our experiments include spectral angle mapper (SAM) [45], RMSE [46], relative dimensionless global error in synthesis (ERGAS) [47], correlation coefficient (CC) [44], universal image quality index (UIQI) [48], adaptive cosine/coherent estimator (ACE) [49,50].

Experimental Results
The proposed HS image simulation method is compared with PHITA [5], Liu's method [19], and Chen's method [18]. All the simulation methods are realized by programing in MATLAB. In the experiments, we test the four methods on each dataset. Figures 5-8 present the result images of the four datasets, respectively. As the test datasets have long and narrow shapes, the images presented in these figures are parts of the results. Table 2 provides the color qualities of the result images. Table  3 show the quantitative evaluations and the statistics of the four datasets, including means and standard deviations (Stds). In this table, the method with the best quantitative performance is represented by a bold font. The key parameters of our method in the experiments are selected by sensitivity analysis, as presented in Section 4.2.

Spatial Performance
We initially compare the simulated images with the reference HS images on the basis of their visual appearance. Figure 5, 7, and 8 present the natural color images of Datasets 1, 3, and 4 (R: 29th band, G: 20th band, B: 12th band), respectively. Figure 6 shows the false color images of Dataset 2 (R: 100th band, G: 75th band, B: 55th band). Each simulation method provides a practicable result. Chen's method simulates images with larger distortions than the other methods. In Figure 5(c), the colors of The quality assessment of the simulated HS images primarily focuses on the visual appearance and quantitative evaluation. The visual appearance provides a general idea of the results. The quantitative evaluation provides fine and objective comparisons. The evaluation indices in our experiments include spectral angle mapper (SAM) [45], RMSE [46], relative dimensionless global error in synthesis (ERGAS) [47], correlation coefficient (CC) [44], universal image quality index (UIQI) [48], adaptive cosine/coherent estimator (ACE) [49,50].

Experimental Results
The proposed HS image simulation method is compared with PHITA [5], Liu's method [19], and Chen's method [18]. All the simulation methods are realized by programing in MATLAB. In the experiments, we test the four methods on each dataset. Figures 5-8 present the result images of the four datasets, respectively. As the test datasets have long and narrow shapes, the images presented in these figures are parts of the results. Table 2 provides the color qualities of the result images. Table 3 show the quantitative evaluations and the statistics of the four datasets, including means and standard deviations (Stds). In this table, the method with the best quantitative performance is represented by a bold font. The key parameters of our method in the experiments are selected by sensitivity analysis, as presented in Section 4.2. Sensors 2018, 18, x FOR PEER REVIEW 11 of 22   The spatial performance of the simulated images is then quantitatively evaluated referring to the real HS image. The main spatial quality evaluation indexes are CC and UIQI in our experiments. Chen's method achieves the lowest CCs and UIQIs in Table 3. The results are consistent with the The spatial performance of the simulated images is then quantitatively evaluated referring to the real HS image. The main spatial quality evaluation indexes are CC and UIQI in our experiments. Chen's method achieves the lowest CCs and UIQIs in Table 3. The results are consistent with the

Spatial Performance
We initially compare the simulated images with the reference HS images on the basis of their visual appearance. Figures 5, 7 and 8 present the natural color images of Datasets 1, 3, and 4 (R: 29th band, G: 20th band, B: 12th band), respectively. Figure 6 shows the false color images of Dataset 2 (R: 100th band, G: 75th band, B: 55th band). Each simulation method provides a practicable result. Chen's method simulates images with larger distortions than the other methods. In Figure 5c, the colors of vegetation and water are considerably different from the reference HS image. In Figure 6c, the spatial details, such as urban area and river, are not efficiently preserved. In Figures 7c and 8c, the farmland is falsely reconstructed with distorted textures and coclors. UPDM provides better visual performance than Chen's method. The images (d) in Figures 5-8 have similar boundaries and textures as the reference HS image. The dense buildings in Figure 5 and the river in Figure 6 are preserved in the results. However, UPDM still produces slight color changes. The colors of vegetation, water, and soil are slightly brighter than those in the reference, as shown in Figures 5, 6 and 8. PHITA and the proposed method provide the two best results. The topographic features and textures of these two results are nearly the and 0.742 (Table 3), respectively, which are higher than those of Chen's method. The best two results in the comparisons are produced by the proposed method and PHITA. The proposed method achieves the highest CC and UIQI, as shown in Table 3. The highest mean CC (0.905) and UIQI (0.842) in Table 3 reveal that our method generates the best quality spatial information in the four datasets. This quantitative evaluation can also be confirmed by visual appearance. Figure 7b presents similar same as those of the references. The proposed method achieves higher visual quality in Datasets 2 with more realistic colors in water and vegetation than PHITA. We also compute the CIEDE2000 of the RGB bands in each dataset ( Table 2). The CIEDE2000 of each method conforms to their visual appearance. The CIEDE2000 of the proposed method is least in most cases (Dataset 1, 2, and 4). PHITA achieves comparable performance. UPDM and Chen's method rank behind the proposed method and PHITA. In addition, the simulated images provide more clear textures than the reference in Figures 5 and 6. The noises and vertical stripes of the simulated images are also removed in the figures.
The spatial performance of the simulated images is then quantitatively evaluated referring to the real HS image. The main spatial quality evaluation indexes are CC and UIQI in our experiments. Chen's method achieves the lowest CCs and UIQIs in Table 3. The results are consistent with the geometry and texture information in the mountain area ( Figure 5c) and river (Figure 6c). The mean CC and UIQI of Chen' method are 0.683 and 0.484 (Table 3), respectively, which are worse than those of other methods. UPDM provides more spatial details and less color distortion in the result images, especially in Figures 6d, 7d and 8d, than Chen's method. The mean CC and UIQI of UPDM are 0.896 edges and colors of the farmland. In Figure 6b, the textures of mountain and urban areas are even clearer than those of the reference. In Figure 8b, most of the spatial information is well generated. PHITA provides comparable results in the figures. In quantitative assessment, the mean CC (0.901) and UIQI (0.841) of PHITA are lower than 0.905 and 0.842, which are achieved by our method.
We further subtract the CCs of UPDM's results from those of the proposed method's and PHITA's results to show the superiority of our method further clearly. As the CCs of Chen's method are considerably less than those of other methods, we do not present its results. Figure 9 shows the subtraction. A positive value means that the CC performance is better than that of UPDM's result, whereas a negative value is the opposite. Figure 10 shows the RMSE of each band of PHITA's and our method's results. The Figures 9 and 10 confirm that our method has higher qualities than other methods in many bands. geometry and texture information in the mountain area ( Figure 5c) and river (Figure 6c). The mean CC and UIQI of Chen' method are 0.683 and 0.484 (Table 3), respectively, which are worse than those of other methods. UPDM provides more spatial details and less color distortion in the result images, especially in Figures 6d, 7d, and 8d, than Chen's method. The mean CC and UIQI of UPDM are 0.896 edges and colors of the farmland. In Figure 6b, the textures of mountain and urban areas are even clearer than those of the reference. In Figure 8b, most of the spatial information is well generated. PHITA provides comparable results in the figures. In quantitative assessment, the mean CC (0.901) and UIQI (0.841) of PHITA are lower than 0.905 and 0.842, which are achieved by our method. We further subtract the CCs of UPDM's results from those of the proposed method's and PHITA's results to show the superiority of our method further clearly. As the CCs of Chen's method are considerably less than those of other methods, we do not present its results. Figure 9 shows the subtraction. A positive value means that the CC performance is better than that of UPDM's result, whereas a negative value is the opposite. Figure 10 shows the RMSE of each band of PHITA's and our method's results. The Figure 9 and 10 confirm that our method has higher qualities than other methods in many bands.

Spectral Performance
We present the SAM distribution images of the Datasets 1-3 in Figures 11-13 to visualize the spectral qualities. The pixels with low SAM errors are represented by blue, whereas those with high SAM errors are represented by red. Chen's method has the most yellow and red pixels in Figures 11-13, thereby showing the largest spectral distortion in the comparisons. UPDM produces pixels with large SAM error in global images. For example, Figure 11c consists of light blue and red, indicating that the spectral information of water is not accurately constructed. The images resulted from PHITA have more dark blue pixels and less red pixels, especially in vegetation and water, than those from the other methods. As shown in Figure 11d, the upper mountain area covered by vegetation has minimal SAM error. Therefore, the spectral distortions generated by PHITA are considerably less than those of Chen's method and UPDM. The proposed method generates a large number of dark blue pixels and little red pixels. Our method constructs more accurate spectrums in vegetation area and water than PHITA. Figures 11a and 12a contain more dark blue area than the other images. In Figures 11a and 13a, the sea and river contain less red and yellow pixels than those of PHITA's results. Meanwhile, the proposed method produces the least red pixels in Figure 11, which consists of multiple materials and abundant texture features. Thus, our method can reduce spectral distortions in complex areas and has significant advantage in generating high-spectral-quality results.
In quantitative evaluation, the main spectral quality measurements include SAM, ERGAS, RMSE and ACE. Chen's method provides the poorest spectral quality performance in four datasets (Table 3). UPDM reconstructs a spectrum in each pixel by combining standard spectrums in the library. The mean SAM, ERGAS, RMSE, and ACE of this method are 10.92, 57.005, 735.1, and 0.08 (Table 3), respectively, which are all better than those of Chen's method. PHITA obtains good spectral qualities in experiments. In Datasets 1, 3, and 4, the spectral indexes of PHITA are all better than those of Chen's method and UPDM. In Dataset 2, the SAM and ERGAS of PHITA are 3.196 and 13.310, respectively, which are the best in the comparison. The best overall spectral quality on the four datasets is achieved by the proposed method. The mean SAM (5.986), ERGAS (16.817), RMSE (284.6), and ACE (0.165) of our method rank first in Table 3. The best ACE performances are achieved by our method, which demonstrates that our method could preserve subpixel signatures better. The spectral quality of our method is slightly worse than that of PHITA in Dataset 2. The reason is that some of the spectrums in the simulated area are different from the spectrums of the same classes in the training area. The spectral endmembers are extracted from the training area. When we use endmembers to reconstruct these special spectrums, large errors may be obtained. For instance, we present the vegetation and water spectrums of the training and simulated areas in Figure 14. The spectral features of water in training and simulated areas are considerably different, such as the 1st-25th and 40th-60th bands. In Figure 14c, the reflectance values of water in the simulated area are considerably less than those in the training area. Therefore, these special spectrums are constructed with large errors, which reduce the global spectral quality on Dataset 2. In Dataset 4, the spectral qualities of our method and PHITA are a little worse than the other three datasets. The reason is that when the distance between the training and test areas is large, the difference between the spectral information in the two images will increases. Thus, error in the simulation image may be increased. However, if the distance is not too large, the error can be accepted.

Overall Quality
Chen's method loses considerable spatial details and produces large spectral distortions in the results. UPDM ranks third among the compared methods. This method constructs detailed textures and geometric information but still leads to spectral distortions. PHITA provides good spectral and spatial qualities but is worse than our method in quantitative evaluations. The proposed method achieves the best performance in the global evaluation ( Table 3). The advantages of the proposed method are presented as follows: 1. The proposed method builds the spectral relation between the MS and HS bands in the same wavelength ranges. The simulated images can be generated by using prior HS endmembers extracted from training HS images. In this manner, the relations between the bands of

Overall Quality
Chen's method loses considerable spatial details and produces large spectral distortions in the results. UPDM ranks third among the compared methods. This method constructs detailed textures and geometric information but still leads to spectral distortions. PHITA provides good spectral and spatial qualities but is worse than our method in quantitative evaluations. The proposed method achieves the best performance in the global evaluation ( Table 3). The advantages of the proposed method are presented as follows: 1. The proposed method builds the spectral relation between the MS and HS bands in the same wavelength ranges. The simulated images can be generated by using prior HS endmembers extracted from training HS images. In this manner, the relations between the bands of

Overall Quality
Chen's method loses considerable spatial details and produces large spectral distortions in the results. UPDM ranks third among the compared methods. This method constructs detailed textures and geometric information but still leads to spectral distortions. PHITA provides good spectral and spatial qualities but is worse than our method in quantitative evaluations. The proposed method achieves the best performance in the global evaluation ( Table 3). The advantages of the proposed method are presented as follows: 1.
The proposed method builds the spectral relation between the MS and HS bands in the same wavelength ranges. The simulated images can be generated by using prior HS endmembers extracted from training HS images. In this manner, the relations between the bands of spectral endmembers are preserved, and fine spectral features are achieved in the simulated images.

2.
The proposed method reconstructs the simulated images by combining the spectrums of related materials, pixel by pixel. Following this strategy, the simulation of each pixel is independent, which can improve the spectral quality of the areas with complex materials and objects.

3.
We utilize iterative schemes to optimize the endmembers and abundance matrixes of the images. This method reduces the residual error in unmixing and reconstruction, thereby further improving the global quality of the results. Our method achieves the best performance in Table 3.
All simulation methods are realized using MATLAB. A laptop computer with a 10 G memory, Intel Core i5 CPU, and Windows 10 system is used for the experiments. The processing time of the proposed method, UPDM, Chen's method and PHITA are presented in Table 4. Our method consumes more time than UPDM (52.55.64s), PHITA (18.23s) and Chen's method (28.84s). However, this limitation can be acceptable in many application fields. The proposed HS image simulation method can provide results with higher spatial and spectral qualities than Chen's method, UPDM, and PHITA with acceptable efficiency.

Discussion
In this section, we discuss the features of our approach. Afterward, a detailed sensitivity analysis is conducted to select the optimal parameters of the four datasets.

Result Analysis
As shown in Section 3.2, the proposed method achieves better performance than Chen's method, UPDM, and PHITA in spatial and spectral qualities. The main features of our approach focus on using the prior MS and HS endmembers and independently constructing pixels and the iteration schemes based on multiplicative update rules.
The proposed method uses the endmembers extracted from training images in HS image simulation. The pixels are produced by the prior spectral information of the similar materials. Following this strategy, all bands in a pixel are constructed simultaneously by the same endmembers and coefficients. Therefore, the fine spectral features of the materials can be efficiently restored. PHITA utilizes different MS bands and linear coefficients to produce each HS band independently. The relations between the simulated bands may be not well considered, which easily leads to spectral distortions in the result. In addition, the proposed method achieves better performance in high spatial resolution image. The reason may be the images with high spatial resolution contain more pure pixels, which can be simulated easily.
For the areas with complex and multi-class objects, the HS image simulation may suffer from large errors because of the various spectral features and high heterogeneity. Our method produces HS images pixel by pixel. Thus, each pixel is generated independently and is slightly influenced by the neighboring pixels and environment. The pixels with special spectral information or textures can be efficiently processed through our method. Chen's method provides the results by classifying the test MS image and replacing the MS information. In this method, the pixels in complex areas are easily misclassified. Consequently, the spectral and spatial qualities of the results will be significantly reduced.
The iteration schemes of our method can effectively optimize the HS and MS endmembers. After alternately updating, the spectral endmembers factorize the images with less residuals, indicating that the endmembers can compose most of the spectrums in the images. Thus, the spectral endmembers become further suitable for the training and test images after our iteration schemes. Then, the accuracy of the results can be improved. UPDM simulated images are also based on spectral unmixing. However, the standard spectrums are often different from the actual spectrums in the study area, which easily leads to large residuals in pixel unmixing. The large residuals of UPDM reduce the quality of the result.

Sensitivity Analysis
In the proposed method, the accuracy of the simulated images is significantly influenced by the parameters, including the inner iteration times I 1 , the outer iteration times I 2 , and the number of spectral endmembers n. To find appropriate parameters of the proposed method, we initially provide a detailed parameter analysis on Dataset 1. Then, the parameters of other datasets are selected by fine-tuning those of Dataset 1. In addition, finding the optimal parameters needs testing all the parameters by lots of experiments. The number of experiments in this section could not ensure that the selected parameters are optimal. The parameters we provided is just suitable for the test datasets.

Iteration Times
The optimal endmember and abundance matrixes are obtained through iterative calculation. When the iteration times are increasing, the factorized matrixes of images will be gradually optimized. Therefore, we must find suitable times for obtaining optimal results with low computational cost. The iterative times in our method include I 1 and I 2 . I 1 is the iteration time of each equation. I 2 is the repeat time of the two rounds of computations. Figure 15 presents the quantitative results with different iteration times. When I 1 is larger than 250, the evaluation indexes vibrate near stable values. When I 2 reaches 5, the evaluation indexes are not continuously improved. Therefore, to achieve the optimal result with less calculation amount, we set I 1 = 250 and I 2 = 5 in the experiment on Dataset 1.  We then extend the obtained iteration parameters to the other datasets by fine tuning. In Dataset 2, the iteration times of our method are set to I 1 = 150 and I 2 = 5. The parameters on Dataset 3 and 4 are both I 1 = 250 and I 2 = 5.

Number of Endmembers
In our approach, a large number of endmembers can provide abundant spectral information, which helps improve the accuracy of the simulated images. However, massive endmembers will not contribute to the quality of the results. Meanwhile, numerous endmembers may increase the computation complexity of our method. Thus, we performed detailed experiments to find the suitable number of endmembers. Figure 16 shows the results. If the number of endmembers is small, then the result HS images have large errors, such as n = 5. With the increase in n, the evaluation indexes are gradually improved. When n > 40, the result qualities are not continuously increased. Therefore, we set the number of endmembers to 40 in the experiments on Dataset 1.   To find the suitable numbers of endmembers on Datasets 2 and 3, we conduct a detailed analysis. After parameter tuning, the suitable numbers of endmembers on Datasets 2, 3, and 4 are determined as 30, 40, 40, respectively. When an unknown dataset is processed by our method, the initial number of the endmembers is suggested to be about 40. Then, this number can be fine-tuned referring to the result quality and computing time to obtain an appropriate value.

Conclusions
In this paper, we propose a novel HS image simulation method based on NMF. Our main contributions are developing the spectral transformation matrix between the HS and MS endmembers and designing the simulation method for constructing the HS image. The spectral transformation matrix is estimated between each MS band and the HS bands in the same wavelength range by using the linear regression model. This matrix helps obtain the corresponding MS and HS endmembers. The proposed method uses iterative computations based on NMF to optimize the extracted endmembers, which are used to simulate the final HS image.
Experiments are performed on four datasets from EO-1 and HJ-1A satellites. Results illustrate that the proposed method can provide simulated images with good spectral and spatial qualities. In comparison with UPDM, PHITA, and Chen's method, the proposed method has advantages in producing clear textures and reducing spectral distortions in the areas with complex and multi-class objects. The four evaluation indexes, namely, SAM, ERGAS, RMSE, CC, UIQI and ACE, in Table 3 are 5.986, 16.817, 284.6, 0.905, 0.842, and 0.165, respectively, which are better than those of other methods.
However, in our approach, if the spectrums in the simulated HS image are considerably different from the prior endmember of the same material in the training images, then these spectrums are difficult to produce precisely. This case may lead to large errors in the results. In addition, the proposed method consumes more time than other methods. In the future, we will focus on correcting the endmembers referring to the spatial location and corresponding MS information in the simulation of each pixel. The acceleration of our method will also be considered in our future work.