METHOD AND DEVICE FOR GENERATING SATELLITE REMOTE-SENSING IMAGE WITH HIGH SPATIAL, TEMPORAL AND SPECTRAL RESOLUTIONS

Abstract

A method and device for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions is provided. The method includes: fusing a multispectral image with high temporal and spatial resolutions but a low spectral resolution and a hyperspectral image with a high spectral resolution but low temporal and spatial resolutions to generate a fused image with high spatial, temporal and spectral resolutions. Compared with a current spatial-temporal-spectral integrated fusion method, the proposed method is more in line with resolution characteristics of current mainstream satellite-borne remote-sensing data and can achieve high-fidelity fusion performance.

Description

TECHNICAL FIELD

The present disclosure relates to the field of optical remote-sensing images, and more specifically, to a method and device for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions.

BACKGROUND

Remote sensing has characteristics of real-time, accurate, and continuous acquisition of large-scale earth surface information, and has enormous application value. It has been widely used in various fields, such as natural resource monitoring, environmental monitoring, mineral identification, precision agriculture, and other fields. However, due to limitations of imaging sensor technologies and costs, a satellite sensor has to make a trade-off between spatial, temporal and spectral resolutions, and it is very difficult to obtain data with high spatial, temporal and spectral resolutions. This limits further application of current remote-sensing data in many fields. As a result, there is more available remote-sensing data, but a relative proportion of data that can be truly used is still low.

In current satellite remote-sensing imaging, a hyperspectral image has a high spectral resolution, but has lower spatial and temporal resolutions than a multispectral image. However, a spectral resolution of the multispectral image is relatively low, and there are generally less than ten bands. China's hyperspectral satellite ZY-1 02D was launched in 2020, and its hyperspectral sensor can capture 166 spectral bands with a spatial resolution of 30 meters, and its revisit period is 55 days. A low spatial resolution (usually 30 m) and a long revisit period (longer than 15 days) make it difficult to finely monitor a rapidly changing scenario based on existing hyperspectral satellite data. Compared with the hyperspectral sensor, a multispectral sensor typically has higher spatial and temporal resolutions. For example, the twin satellites Sentinel-2 a/b launched in 2015 and 2017 can access the earth surface with a 5-day revisit period, providing 4-band data with a spatial resolution of 10 meters. A short revisit period and a fine spatial resolution of the multispectral satellite make it possible to perform rapid monitoring with a long timing sequence, but a small quantity of bands make it difficult to accurately identify a surface feature. Remote-sensing image fusion, as an information processing technology, can systematically synthesize two or more remote-sensing images with complementary information in a same region to generate an image that provides more visual perception and computer processing information than a single image that can be obtained. The remote-sensing image fusion has become one of main means to improve a spatial/temporal/spectral resolution of a single-source satellite. However, current remote-sensing image fusion methods mainly focus on fusion of two of three resolution attributes of a remote-sensing sensor (spatial-spectral fusion, and temporal-spatial fusion). The only few spatial-temporal-spectral fusion methods are mostly used to fuse hyperspectral images having a high temporal resolution, without considering a resolution characteristic of a current satellite-borne remote-sensing image, that is, a characteristic that the spatial and temporal resolutions of the hyperspectral image are lower than those of the multispectral image.

Therefore, based on resolution characteristics of the current remote-sensing data, it is of great significance for the further application of remote sensing data to carry out the research on a multi-source remote-sensing spatial-temporal-spectral collaborative fusion method, which integrates complementary advantages of spatial, temporal and spectral resolutions of a multi-source remote-sensing image and generates image data with high spatial, temporal and spectral resolutions through fusion.

SUMMARY

In order to overcome shortcomings that a current remote-sensing image fusion method insufficiently considers resolution characteristics of current remote-sensing data and cannot fully leverage a collaborative advantage of current multi-source remote-sensing data to generate data with high spatial, temporal and spectral resolutions, the present disclosure provides a method and device for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions to obtain the data with high spatial, temporal and spectral resolutions.

According to a first aspect, a method for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions is provided, including:

• • step 1: obtaining a hyperspectral image H₁ and a multispectral image M₁ at a time point T₁, and a multispectral image M₂ at a time point T₂, and performing preprocessing to obtain an upsampled hyperspectral image Ĥ₁;
  • step 2: based on a local linear constraint, searching for an adjacent similar full-band block M₁(Ω_i,j,k) in the multispectral image M₁ to represent an image block centered at (i, j) in the multispectral image M₂, and obtaining a representation weight function of the similar full-band block; and sharing the representation weight function with the upsampled hyperspectral image Ĥ₁ to obtain an upsampled hyperspectral image Ĥ₂^*∈ at a prediction time point;
  • step 3: obtaining a time variation image represented as T=Ĥ₂^*−Ĥ₁ for an initialized target fused image;
  • step 4: decomposing a target image X₂ at the time point T₂ into an endmember matrix E₂ and an abundance matrix A₂ based on a spectral linear mixed model, where X₂=E₂A₂, and the target image X₂ is an image with high spatial, temporal and spectral resolutions;
  • step 5: establishing observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions;
  • step 6: introducing a generalized linear mixed model, and representing the endmember matrix E₂ as E₂=PeE₁, where the symbol e represents a Hadamard product of element-by-element multiplication, P represents a time variation matrix of the target image, and E₁ represents an endmember matrix at the time point T₁;
  • step 7: obtaining the abundance matrix A₂, the time variation matrix P and a time variation image T of the target image; and
  • step 8: obtaining the target image finally by multiplying the endmember matrix E₁ at the time point T₁ by the time variation matrix P at the time point T₂ element by element and then by the abundance matrix A₂ with high resolution: X=(P^eE₁)A₂.

Preferably, the step 1 includes:

• • step 1.1: obtaining the hyperspectral image H₁∈ and the multispectral image M₁∈ at the time point T₁, and the multispectral image M₂∈ at the time point T₂, where l and L represent quantities of bands, w and W represent quantities of pixels, l<L, and w<W;
  • step 1.2: normalizing the hyperspectral image H₁ and the multispectral images M₁ and M₂, and upsampling a normalized hyperspectral image H₁ to a spatial size of the multispectral image to obtain the Ĥ₁∈; and
  • step 1.3: folding the multispectral images M₁ and M₂ and the upsampled hyperspectral image Ĥ₁ into a three-dimensional form along a spectral dimension, and splitting the three-dimensional form into √{square root over (c)}×√{square root over (c)}×l full-band blocks with both a width and a height being √{square root over (c)} and the spectral dimension being l.

Preferably, in the step 5, the observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions are represented as follows:

M_i=RX_i=RE_iA_i

H_i=X_iFS=E_iA_iFS

where i∈{1,2}; X_i∈, H_i∈, and M_i∈ respectively represent an image with high spatial, temporal and spectral resolutions, a hyperspectral image, and a multispectral image at a time point T_i; E_i∈z,186 and A_i∈ respectively represent an endmember matrix and an abundance matrix of the image with high spatial, temporal and spectral resolutions at the time point T_i, Q represents a quantity of endmembers, and <<L; and R∈ represents a spectral downsampling matrix, and F∈ and S∈ respectively represent a spatial blur matrix and a spatial downsampling matrix.

Preferably, the step 7 includes:

• • step 7.1: separately introducing a regularization term of the abundance matrix A₂, a regularization term of the time variation matrix P and a regularization term of the time variation image T of the target image to obtain an energy function for a fusion problem;
  • step 7.2: splitting the energy function into three sub-functions corresponding to variables of the abundance matrix A₂, the time variation matrix P and the time variation image T, introducing a plurality of auxiliary variables for each sub-function to split the sub-function into sub-problems, and constructing an augmented Lagrangian function for each sub-problem; and
  • step 7.3: using an alternating direction method of multipliers to perform optimal solving on the augmented Lagrangian function to obtain the abundance matrix A₂, the time variation matrix P and the time variation image T of the target image.

Preferably, in the step 7.1, the energy function is represented as follows:

$\arg \underset{A_2,P,T}{\min }\frac{1}{2}{M_2-R\left({\mathrm{PeE}}_1\right)A_2}_2^2+\frac{1}{2}{D-R\left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{1}{2}{\mathrm{RT}-D}_2^2+\frac{1}{2}{T-\left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)}_2^2+\Phi \left(A_2\right)+\Phi \left(P\right)+{\rm Y}\left(T\right)s.t.D=M_2-M_1$

where D represents a differential image of the multispectral images M₂ and M₁. Preferably, the step 7.1 includes:

• • step 7.1.1: introducing the regularization term of the abundance matrix A₂ of the target image, where the regularization term is a spatial smoothing constraint and represented as follows:

$\Phi \left(A_2\right)=\frac{{\lambda }_A}{2}\left({H_h{A}_2}_{2,1}+{H_v{A}_2}_{2,1}\right)$

where H_h and H_v represent linear operators respectively used to calculate horizontal and vertical gradients between components of adjacent pixels in a two-dimensional signal, ∥·∥_2,1 represents a mixed L_2,1 norm, and λ_A represents a regularization parameter of the abundance matrix A₂;

• • step 7.1.2: introducing the regularization term of the time variation matrix P, where the regularization term is a spectral smoothing constraint and represented as follows:

$\Psi \left(P\right)=\frac{{\lambda }_1}{2}{H_{\ell }P}_F^2+\frac{{\lambda }_2}{2}{P-1_L{1}_Q^T}_F^2$

where H₁ represents a differential operator in a spectral dimension direction, 1_L1_Q^T∈ represents a matrix with all elements in L rows and Q columns being 1, and λ₁ and λ₂ represent regularization parameters of the time variation matrix P; and

• • step 7.1.3: introducing the regularization term of the time variation image T, where the regularization term is a low-rank constraint and represented as follows:

${\rm Y}\left(T\right)={\lambda }_T{𝒯}_{\mathrm{TNN}}$

where represents a third-order tensor form of the T; ∥·∥_TNN represents a tensor nuclear norm, with ∥∥_TNN=Σ_i=1³Σ_jσ_j(T⁽ⁱ⁾ where σ_j(T⁽ⁱ⁾) represents a j^th singular value of T⁽ⁱ⁾, and T is obtained by performing Fourier transform on the T; and λ_T represents a regularization parameter of the time variation image T.

Preferably, in the step 7.2, a sub-problem of the abundance matrix A₂ of the target image is as follows:

$A_2=\arg \underset{A_2}{\min }\frac{1}{2}{M_2-R\left(\mathrm{Pe}{E}_1\right)A_2}_2^2+\frac{1}{2}{D-R\left(\left(\mathrm{Pe}{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{\eta }{2}{T-\left(\left(\mathrm{Pe}{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{{\lambda }_A}{2}\left({H_h{A}_2}_{2,1}+{H_v{A}_2}_{2,1}\right)$

where η represents a term parameter of the time variation image and is used to control a contribution of the time variation image to a model; and

• • auxiliary variables B₁=A₂, B₂=H_h(A₁), and B₃=H_v(A₁) are introduced, u=vec(A₂) is set, a symbol vec(·) represents a vectorization operation, and an augmented Lagrangian function for the sub-problem of the abundance matrix A₂ is given as follows:

$ℒ\left(u,B_1,B_2,B_3,V_1,V_2,V_3\right)=\frac{1}{2}{\mathrm{vec}\left(M_2\right)-{\mathrm{rep}}_M\left(R\left(\mathrm{Pe}{E}_1\right)u\right)}_2^2+\frac{1}{2}{\mathrm{vec}\left(D\right)-{\mathrm{rep}}_M\left(R\left(\left(\mathrm{Pe}{E}_1\right)u-\right)\right)+\mathrm{vec}\left(E_1{A}_1\right)}_2^2+\frac{\eta }{2}{\mathrm{vec}\left(T\right)-{\mathrm{rep}}_M\left(\left(\mathrm{Pe}{E}_1\right)u\right)+\mathrm{vec}\left(E_1{A}_1\right)}_2^2+{\lambda }_A\left({B_2}^2+{B_3}^2\right)+\frac{\rho }{2}\left({\mathrm{vec}\left(B_1\right)-u+V_1}_2^2+{\mathrm{vec}\left(B_2\right)-H_h\mathrm{vec}\left(B_1\right)+V_2}_2^2+{\mathrm{vec}\left(B_3\right)-H_v\mathrm{vec}\left(B_1\right)+V_3}_2^2\right)$

where a matrix V_i, i=1,2,3 represents a dual variable of a Lagrangian function; ρ represents a step size, where ρ>0; and rep_M(V) represents a block diagonal matrix of a matrix V and means copying the matrix V along a diagonal for M times.

Preferably, in the step 7.2, a sub-problem of the time variation matrix P of the target image is as follows:

$P=\arg \underset{P}{\min }\frac{1}{2}{M_2-R\left(\mathrm{Pe}{E}_1\right)A_2}_2^2+\frac{1}{2}{D-R\left(\left(\mathrm{Pe}{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{1}{2}{T-\left(\left(\mathrm{Pe}{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{{\lambda }_1}{2}{P-1_L{1}_Q}_2^2+\frac{{\lambda }_2}{2}{H_{l}P}_2^2$

an auxiliary variable B=PeE₁ is introduced, and an augmented Lagrangian function for the sub-problem of the time variation matrix P is as follows:

$ℒ\left(B,P,V\right)=\frac{1}{2}{M_2-RB{A}_2}_2^2+\frac{1}{2}{D-R\left(B{A}_2-E_1{A}_1\right)}_2^2+\frac{\eta }{2}{T-\left(B{A}_2-E_1{A}_1\right)}_2^2+\frac{{\lambda }_1}{2}{P-1_L{1}_Q}_2^2+\frac{{\lambda }_2}{2}{H_{l}P}_2^2+\frac{\rho }{2}\left({B-\mathrm{Pe}{E}_1+V}_2^2\right)$

a sub-problem of the time variation image T of the target image is as follows:

$T=\arg \underset{T}{\min }\frac{1}{2}{\mathrm{RT}-D}_2^2+\frac{1}{2}{T-\left(\left(\mathrm{Pe}{E}_1\right)A_2-E_1{A}_1\right)}_2^2+{\lambda }_T{𝒯}_{\mathrm{TNN}}$

an auxiliary variable = is introduced, and an augmented Lagrangian function for the sub-problem of the time variation image T is as follows: ?

$T=\arg \underset{T}{\min }\frac{1}{2}{\mathrm{RT}-D}_2^2+\frac{1}{2}{T-\left(\left(\mathrm{Pe}{E}_1\right)A_2-E_1{A}_1\right)}_2^2+{\lambda }_T{ℬ}_{\mathrm{TNN}}.$

Preferably, in the step 2, the image block centered at the (i, j) in the multispectral image M₂ is represented using the adjacent similar full-band block M₁(Ω_i,j,k) in the multispectral image

$M_1:M_2\left(i,j\right)=\sum _{{\Omega }_{i,j,k}\in {\Omega }_{i,j}}{\omega }_{i,j,k}{M}_1\left({\Omega }_{i,j,k}\right),$

where ω_i,j,k represents a weight of a k^th image block in the Ω_i,j,k, and the ω_i,j,k is obtained based on the local linear constraint:

$\min {M_2\left(i,j\right)-M_1\left({\Omega }_{i,j}\right){\omega }_{i,j}}^2+\lambda {d_{i,j}\odot {\omega }_{i,j}}^{2}s.t.1^T{\omega }_{i,j}=1$

where ω_i,j represents a weight corresponding to the image block M₂(_i,j), ∥d_i,jeω_i,j ∥² represents a regularization term of the local linear constraint, d_i,j represents a Euclidean distance between a similar pixel and a center of M₁(_i,j), and λ represents a regularization term parameter; and a solution to an optimization problem of the local linear constraint is as follows:

${\omega }_{i,j}=\left(\left(C_{i,j}+\lambda \mathrm{diag}\left(d_{i,j}\right)\right)\backslash 1\right){\left(1^T\left(C_{i,j}+\lambda \mathrm{diag}\left(d_{i,j}\right)\right)\backslash 1\right)}^{-1}$

where C_i,j=(M₁(Ω_i,j)−M₂(i,j)^T)(M₁(Ω_i,j)−1M₂(i,j)^T)^T represents a covariance matrix of the image block, diag (d_i,j) represents a diagonal matrix with d_i,j being a diagonal element; afterwards, the weight ω_i,j obtained from the multispectral image is shared with the upsampled hyperspectral image Ĥ₁ to generate the upsampled hyperspectral image

${\hat{H}}_2^{*}=\sum _{{\Omega }_{i,j,k}\in {\Omega }_{i,j}}{\omega }_{i,j,k}{\hat{H}}_1\left({\Omega }_{i,j,k}\right)$

at the prediction time point.

According to a second aspect, a device for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions is provided, which is configured to execute the method for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions according to any one first aspect, and includes:

• • a first obtaining module configured to obtain a hyperspectral image H₁ and a multispectral image M₁ at a time point T₂, and a multispectral image M₂ at a time point T₂, and perform preprocessing to obtain an upsampled hyperspectral image Ĥ₁;
  • a search module configured to: based on a local linear constraint, search for an adjacent similar full-band block M₁(Ω_i,j,k) in the multispectral image M₁ to represent an image block centered at (i,j) in the multispectral image M₂, and obtain a representation weight function of the similar full-band block; and share the representation weight function with the upsampled hyperspectral image Ĥ₁ to obtain an upsampled hyperspectral image Ĥ₂^*∈ at a prediction time point;
  • a second obtaining module configured to obtain a time variation image represented as T=Ĥ₂^*−Ĥ₁ for an initialized target fused image;
  • a decomposition module configured to decompose a target image X₂ at the time point T₂into an endmember matrix E₂ and an abundance matrix A₂ based on a spectral linear mixed model, where X₂=E₂A₂, and the target image X₂ is an image with high spatial, temporal and spectral resolutions;
  • an establishment module configured to establish observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions;
  • a representation module configured to introduce a generalized linear mixed model, and represent the endmember matrix E₂ as E₂=PeE₁, where the symbol e represents a Hadamard product of element-by-element multiplication, P represents a time variation matrix of the target image, and E₁ represents an endmember matrix at the time point T_i;
  • a third obtaining module configured to obtain the abundance matrix A₂, the time variation matrix P and a time variation image T of the target image; and
  • an obtaining module configured to obtain the target image finally by multiplying the endmember matrix E₁ at the time point T_i by the time variation matrix P at the time point T₂element by element and then by the abundance matrix A₂ with high resolution: X=(P^eE₁)A₂.

The present disclosure has following beneficial effects:

1. The present disclosure mainly focuses on fusing a hyperspectral remote-sensing image with low temporal and spatial resolutions and a multispectral remote-sensing image with high temporal and spatial resolutions. Compared with existing remote-sensing fusion methods, the present disclosure can fuse attributes of spatial, temporal and spectral resolutions at a time to generate a fused image with high spatial, temporal and spectral resolutions. In addition, a fusion task targeted by the present disclosure is more in line with resolution characteristics of current satellite-borne remote-sensing data.

2. A spatial-temporal-spectral fusion method established in the present disclosure provides a generalized fusion framework, including a time variation model and a resolution enhancement reconstruction model for a time variation image. The time variation model is used to depict and model a time variation between images, and the resolution enhancement reconstruction model accurately estimates and obtains the time variation by reconstructing a resolution-enhanced time variation image model. The time variation model and the resolution enhancement reconstruction model for the time variation image are coupled to achieve high-fidelity spatial-temporal-spectral fusion performance.

3. Based on a local linear constraint, the present disclosure represents a multispectral image M₂ at prediction time T₂ using a multispectral image M₁ observed at observation time T₁. After the image is divided into blocks, a similar block M₁(Ω_i,j,k) is searched for around a full-band block of M₁(i,j) at a previous time point to represent a target block M₂(i,j), so as to obtain a weight coefficient. The learned weight coefficient is shared with an upsampled hyperspectral image Ĥ₁ observed at the T₁ (observation time) to generate an upsampled hyperspectral image Ĥ₂^* at the prediction time. This fully utilizes a local self-similarity in a remote-sensing image to learn a time variation feature.

4. The present disclosure establishes a time variation image model. An initialized time variation image generated by Ĥ₂ ^* and Ĥ₁ includes a time variation characteristic of a target image between two time points, providing accurate time variation information and hyperspectral resolution information for subsequent estimation of the target image.

5. The present disclosure depicts the time variation based on a generalized linear mixed model, and then coupled with the resolution enhancement reconstruction model for the time variation image to propose a spatial-temporal-spectral fusion model. The estimation of the target image is transformed into estimation of an abundance matrix A₂, a time variation matrix P and a time variation image T of the target image. The time variation matrix P in the generalized linear mixed model can capture the time variation of the target image, and the estimation of the time variation image T can serve as a constraint for obtaining the time variation.

6. The present disclosure introduces a spatial smoothing prior for abundance of the target image, a spectral smoothing prior for the time variation, and a low-rank prior for the time variation image model. The introduction of prior information enables the model to obtain an accurate solution, achieving a high-fidelity fused image X.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram of generating an initialized time variation image based on a local linear constraint according to the present disclosure; and

FIG. 2 is a schematic diagram of a spatial-temporal-spectral integrated fusion model according to the present disclosure.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The present disclosure will be further described below with reference to embodiments. The following description of the embodiments is only for helping to understand the present disclosure. It should be noted that, several improvements and modifications may be made by a person of ordinary skill in the art without departing from the principle of the present disclosure, and these improvements and modifications should also fall within the protection scope of the present disclosure.

Embodiment 1

A spatial-temporal-spectral fusion problem addressed by the present disclosure is to predict target image X₂∈ with high spatial, temporal and spectral resolutions at prediction time T₂ based on multispectral image M₁∈ (with l bands and W pixels) and hyperspectral image H₁∈ (with L bands and w pixels) that are observed at observation time T₁, and multispectral image M₂∈ observed at the T2, where l<L, and w<W. A method provided in the present disclosure includes the following steps:

Step 1: Hyperspectral image H₁ and multispectral image M₁ at time point T₁, and multispectral image M₂ at time point T₂ are obtained and preprocessed to obtain upsampled hyperspectral image Ĥ₁.

The step 1 includes the following sub-steps:

Step 1.1: The hyperspectral image H₁∈ and the multispectral image M₁∈ at the time point T₁, and the multispectral image M₂∈ at the time point T₂ are obtained, where l and L represent quantities of bands, w and W represent quantities of pixels, l<L, and w<W.

Step 1.2: The hyperspectral image H₁ and the multispectral images M₁ and M₂ are normalized, and normalized hyperspectral image H₁ is upsampled to a spatial size of the multispectral image to obtain the Ĥ₁ ∈, where an upsampling method for the hyperspectral image H₁ may be bilinear interpolation.

Step 1.3: The multispectral images M₁ and M₂ and the upsampled hyperspectral image Ĥ₁ are folded into a three-dimensional form along a spectral dimension, and the three-dimensional form is split into √{square root over (c)}×√{square root over (c)}×l full-band blocks.

Step 2: Based on a local linear constraint, adjacent similar full-band block M₁(Ω_i,j,k) in the multispectral image M₁ is searched for to represent an image block centered at (i, j) in the multispectral image M₂, and a representation weight function of the similar full-band block is obtained; and the representation weight function is shared with the upsampled hyperspectral image Ĥ₁ to obtain upsampled hyperspectral image Ĥ₂^*∈ at a prediction time point.

In the step 2, the image block centered at the (i,j) in the multispectral image M₂ is represented using the adjacent similar full-band block M₁(Ω_i,j,k) in the multispectral image M₁:

$M_2\left(i,j\right)=\sum _{{\Omega }_{i,j,k}\in {\Omega }_{i,j}}{\omega }_{i,j,k}{M}_1\left({\Omega }_{i,j,k}\right),$

where ω_i,j,k represents a weight of a k^th image block in the Ω_i,j,k, and the ω_i,j,k is obtained based on the local linear constraint:

$\min {M_2\left(i,j\right)-M_1\left({\Omega }_{i,j}\right){\omega }_{i,j}}^2+\lambda {d_{i,j}e{\omega }_{i,j}}^{2}s.t.1^T{\omega }_{i,j}=1$

where ω_i,j represents a weight corresponding to the image block M₂(i,j), ∥d_i,jeω_i,j∥² represents a regularization term of the local linear constraint, d_i,j represents a Euclidean distance between a similar pixel and a center of M₁(i,j), and λ represents a regularization term parameter; and a solution to an optimization problem of the local linear constraint is as follows:

${\omega }_{i,j}=\left(\left(C_{i,j}+\lambda \mathrm{diag}\left(d_{i,j}\right)\right)\backslash 1\right){\left(1^T\left(C_{i,j}+\lambda \mathrm{diag}\left(d_{i,j}\right)\right)\backslash 1\right)}^{-1}$

where C_i,j=(M₁(Ω_i,j)−M₂(i,j)^T)(M₁(Ω_i,j)−1M₂(i,j)^T)^T represents a covariance matrix of the image block, diag(d_i,j) represents a diagonal matrix with d_i,j being a diagonal element; afterwards, the weight ω_i,j obtained from the multispectral image is shared with the upsampled hyperspectral image Ĥ₁ to generate the upsampled hyperspectral image

${\hat{H}}_2^{*}=\sum _{{\Omega }_{i,j,k}\in {\Omega }_{i,j}}{\omega }_{i,j,k}{\hat{H}}_1\left({\Omega }_{i,j,k}\right)$

at the prediction time point.

Step 3: A time variation image represented as T=Ĥ₂^*−Ĥ₁ is obtained for an initialized target fused image.

Step 4: Target image X₂ at the time point T₂ is decomposed into endmember matrix E₂ with retained spectral information and abundance matrix A₂ with retained spatial information of the target image based on a spectral linear mixed model, where X₂=E₂A₂, and the target image X₂ is an image with high spatial, temporal and spectral resolutions.

Step 5: Observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions are established.

In the step 5, the observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions are represented as follows:

M_i=RX_i=RE_iA_i

H_i=X_iFS=E_iA_iFS

where i∈{1,2}; X_i∈, H_i∈, and M_i∈ respectively represent an image with high spatial, temporal and spectral resolutions, a hyperspectral image, and a multispectral image at a time point T_i; E_i∈ and A_i∈ respectively represent an endmember matrix and an abundance matrix of the image with high spatial, temporal and spectral resolutions at the time point T_i, Q represents a quantity of endmembers, and <<L; and R∈ represents a spectral downsampling matrix, and F∈ and S∈ respectively represent a spatial blur matrix and a spatial downsampling matrix.

Step 6: A generalized linear mixed model is introduced, and the endmember matrix E₂ is represented as E₂=PeE₁, where the symbol e represents a Hadamard product of element-by-element multiplication, P represents a time variation matrix of the target image, and E₁ represents an endmember matrix at the time point T₁.

In the step 6, after the generalized linear mixed model is introduced, an estimation problem of the target image X₂ is transformed into estimation of the abundance matrix A₂ of the target image and the time variation matrix P of the target image.

Step 7: The abundance matrix A₂, the time variation matrix P, and time variation image T of the target image are obtained.

In the step 7, initialized T is obtained based on T⁽⁰⁾=Ĥ₂^*−Ĥ₁. Initialized E₁ is extracted from the hyperspectral image H₁ by using an endmember extraction algorithm. The time variation matrix P is initialized as P⁽⁰⁾=1_L1_Q^T. Initialized abundance A₂ of the target image is obtained by classifying the multispectral image M₂ or by spectral unmixing.

The step 7 includes the following sub-steps:

Step 7.1: A regularization term of the abundance matrix A₂, a regularization term of the time variation matrix P and a regularization term of the time variation image T of the target image are separately introduced to obtain an energy function for a fusion problem.

In the step 7.1, the energy function is represented as follows:

$\arg \underset{A_2,P,T}{\min }\frac{1}{2}{M_2-R\left(Pe{E}_1\right)A_2}_2^2+\frac{1}{2}{D-R\left(\left(Pe{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{1}{2}{\mathrm{RT}-D}_2^2+\frac{1}{2}{T-\left(\left(Pe{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\Phi \left(A_2\right)+\Psi \left(P\right)+{\rm Y}\left(T\right)$ $s.t.D=M_2-M_1$

The step 7.1 includes the following sub-steps:

Step 7.1.1: The regularization term of the abundance matrix A₂ of the target image is introduced, where the regularization term is a spatial smoothing constraint and represented as follows:

$\Phi \left(A_2\right)=\frac{{\lambda }_A}{2}\left({H_h{A}_2}_{2,1}+{H_v{A}_2}_{2,1}\right)$

where H_h and H_v represent linear operators respectively used to calculate horizontal and vertical gradients between components of adjacent pixels in a two-dimensional signal, ∥·∥_2,1 represents a mixed L_2,1 norm, and λ_A represents a regularization parameter of the abundance matrix A₂.

Step 7.1.2: The regularization term of the time variation matrix P is introduced, where the regularization term is a spectral smoothing constraint and represented as follows:

$\Psi \left(P\right)=\frac{{\lambda }_1}{2}{H_{\ell }P}_F^2+\frac{{\lambda }_2}{2}{P-1_L{1}_Q^T}_F^2$

where H₁ represents a differential operator in a spectral dimension direction, and λ₁ and λ₂ represent regularization parameters of the time variation matrix P.

Step 7.1.3: The regularization term of the time variation image T is introduced, where the regularization term is a low-rank constraint and represented as follows:

${\rm Y}\left(T\right)={\lambda }_T{𝒯}_{TNN}$

where represents a third-order tensor form of the T; ∥·∥_TNN represents a tensor nuclear norm, with ∥∥_TNN=Σ_i=1³Σ_jσ_j(T⁽ⁱ⁾), where σ_j(T⁽ⁱ⁾) represents a j^th singular value of T, and T is obtained by performing Fourier transform on the T; and λ_T represents a regularization parameter of the time variation image T.

In the step 7.1, the present disclosure couples a generalized spectral linear mixed model and a time variation image reconstruction model, and introduces the abundance matrix A₂, the time variation matrix P of the target image, and a prior term of reconstructed time variation image T to construct a spatial-temporal-spectral fusion model.

Step 7.2: The energy function is split into three sub-functions corresponding to variables of the abundance matrix A₂, the time variation matrix P and the time variation image T, a plurality of auxiliary variables are introduced for each sub-function to split the sub-function into sub-problems, and an augmented Lagrangian function is constructed for each sub-problem.

In the step 7.2, a splitting strategy used in the present disclosure is to separately solve one of the variables while fixing the other two variables, and use an alternating least squares method to minimize a loss function corresponding to each variable, in order to obtain a local stationary point.

In the step 7.2, a sub-problem of the abundance matrix A₂ of the target image is as follows:

$A_2=\arg \underset{A_2}{\min }\frac{1}{2}{M_2-R\left(Pe{E}_1\right)A_2}_2^2+\frac{1}{2}{D-R\left(\left(Pe{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{\eta }{2}{T-\left(\left(Pe{E}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{{\lambda }_A}{2}\left({H_h{A}_2}_{2,1}+{H_v{A}_2}_{2,1}\right)$

auxiliary variables B₁=A₂, B₂=H_h(A₁), and B₃=H_v(A₁) are introduced, u=vec(A₂) is set, symbol vec(·) represents a vectorization operation, and an augmented Lagrangian function for the sub-problem of the abundance matrix A₂ is given as follows:

$ℒ\left(u,B_1,B_2,B_3,V_1,V_2,V_3\right)=\frac{1}{2}{\mathrm{ve}c\left(M_2\right)-{\mathrm{rep}}_M\left(R\left(Pe{E}_1\right)u\right)}_2^2+\frac{1}{2}{\mathrm{ve}c\left(D\right)-{\mathrm{rep}}_M\left(R\left(\left(Pe{E}_1\right)u-\right)\right)+vec\left(E_1{A}_1\right)}_2^2+\frac{\eta }{2}{\mathrm{ve}c\left(T\right)-{\mathrm{rep}}_M\left(\left(Pe{E}_1\right)u\right)+vec\left(E_1{A}_1\right)}_2^2+{\lambda }_A\left({B_2}^2+{B_3}^2\right)+\frac{\rho }{2}\left({\mathrm{ve}c\left(B_1\right)-u+V_1}_2^2+{\mathrm{ve}c\left(B_2\right)-H_{h}vec\left(B_1\right)+V_2}_2^2+{\mathrm{ve}c\left(B_3\right)-H_{v}vec\left(B_1\right)+V_3}_2^2\right)$

where matrix V_i, i=1,2,3 represents a dual variable of a Lagrangian function; ρ represents a step size, where ρ>0; and rep_M(V) represents a block diagonal matrix of matrix V and means copying the matrix V along a diagonal for M times.

The present disclosure splits a solution of a Lagrangian function of the abundance matrix A₂ into sub-problems about variables u, B₁,B₂,B₃,V₁,V₂,V₃.

The sub-problem of the variable u of the abundance matrix A₂ is as follows:

$u\in \arg \underset{u}{\min }ℒ\left(u,B_1,B_2,B_3,V_1,V_2,V_3\right)=\arg \underset{u}{\min }\frac{1}{2}{\mathrm{ve}c\left(M_2\right)-{\mathrm{rep}}_M\left(R\left(Pe{E}_1\right)u\right)}_2^2+\frac{1}{2}{\mathrm{ve}c\left(D\right)-{\mathrm{rep}}_M\left(R\left(\left(Pe{E}_1\right)u\right)\right)+vec\left(E_1{A}_1\right)}_2^2+\frac{\eta }{2}{\mathrm{ve}c\left(T\right)-{\mathrm{rep}}_M\left(\left(Pe{E}_1\right)u\right)+vec\left(E_1{A}_1\right)}_2^2+\frac{\rho }{2}{\mathrm{ve}c\left(B_1\right)-u+V_1}_2^2$

A solution to the above equation is: u*=Ω₁⁻¹Ω₂, where Ω₁ and Ω₂ are given according to the following equations:

${\Omega }_1=2{\left(R\left(Pe{E}_1\right)\right)}^{T}R\left(Pe{E}_1\right)+{\eta \left(Pe{E}_1\right)}^T\left(Pe{E}_1\right)+\rho I$ ${\Omega }_2={\left(R\left(Pe{E}_1\right)\right)}^T\left(M_2+D+H_1\right)+{\left(Pe{E}_1\right)}^T\left(T+E_1{A}_1\right)+\rho \left(B_1+V_1\right)$

The sub-problem of the variable B₁ of the abundance matrix A₂ is as follows:

$B_1\in \arg \underset{B_1}{\min }ℒ\left(u,B_1,B_2,B_3,V_1,V_2,V_3\right)=\arg \underset{B_1}{\min }\frac{\rho }{2}\left({\mathrm{ve}c\left(B_1\right)-u+V_1}_2^2+{\mathrm{ve}c\left(B_2\right)-H_{h}vec\left(B_1\right)+V_2}_2^2+{\mathrm{ve}c\left(B_3\right)-H_{v}vec\left(B_1\right)+V_3}_2^2\right)$

A solution to the above equation is:

$vec\left(B_1\right)={\left(I+H_h^T{H}_h+H_v^T{H}_v\right)}^{-1}\left(u-V_1+H_h^{T}vec\left(B_2\right)+H_v^T{V}_2+H_v^{T}vec\left(B_3\right)+H_v^T{V}_3\right)$

The symbol vec(·) in the above equation represents an operator that converts a matrix into a vector.

The sub-problem of the variable B₂ of the abundance matrix A₂ is as follows:

$B_2\in \arg \underset{B_2}{\min }ℒ\left(u,B_1,B_2,B_3,V_1,V_2,V_3\right)=\arg \underset{B_2}{\min }\frac{{\lambda }_A}{2}{H_h{A}_2}_{2,1}+{\mathrm{ve}c\left(B_2\right)-H_{h}vec\left(B_1\right)+V_2}_2^2$

A solution to the above equation is equivalent to a proximal operator of an L₂ norm. This solution can be obtained by using a soft-threshold operator, which is expressed as follows:

$B_2^{*}={\mathrm{soft}}_{{\lambda }_A/\rho }\left(H_{h}vec\left(B_1\right)+V_2\right)$

In the above equation, the soft-threshold operator is defined as follows:

${\mathrm{soft}}_a\left(b\right)=\{\begin{array}{cc}\max \left(1-\frac{a}{b},0\right)b,& b\ne 0\\ 0,& b=0\end{array}$

The sub-problem of the variable B₃ of the abundance matrix A₂ is as follows:

$B_3\in \arg \underset{B_3}{\min }ℒ\left(u,B_1,B_2,B_3,V_1,V_2,V_3\right)=\arg \underset{B_3}{\min }\frac{{\lambda }_A}{2}{H_h{A}_2}_{2,1}+{\mathrm{ve}c\left(B_2\right)-\left(B_1\right)+V_2}_2^2$

Similar to an optimization strategy of the sub-problem of the B₂, a solution is obtained by a soft threshold:

$B_3^{*}={\mathrm{soft}}_{{\lambda }_A/\rho }\left(H_v\mathrm{vec}\left(B_1\right)+V_3\right)$

An update rule for the dual variable of the variable of the abundance matrix A₂ is as follows:

$V_1\leftarrow V_1+\mathrm{vec}\left(B_1\right)-u$ $V_2\leftarrow V_2+vec\left(B_2\right)-H_{h}vec\left(B_1\right)$ $V_3\leftarrow V_3+vec\left(B_3\right)-H_{v}vec\left(B_1\right)$

A sub-problem of the time variation matrix P of the target image is as follows:

$P=\arg \underset{P}{\min }\frac{1}{2}{M_2-R\left({\mathrm{PeE}}_1\right)A_2}_2^2+\frac{1}{2}{D-R\left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{1}{2}{T-\left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{{\lambda }_1}{2}{P-1_L{1}_Q}_2^2+\frac{{\lambda }_2}{2}{H_{l}P}_2^2$

Auxiliary variable B=PeE₁ is introduced, and an augmented Lagrangian function for the sub-problem of the time variation matrix P is as follows:

$ℒ\left(B,P,V\right)=\frac{1}{2}{M_2-{\mathrm{RBA}}_2}_2^2+\frac{1}{2}{D-R\left({\mathrm{BA}}_2-E_1{A}_1\right)}_2^2+\frac{\eta }{2}{T-\left({\mathrm{BA}}_2-E_1{A}_1\right)}_2^2+\frac{{\lambda }_1}{2}{P-1_L{1}_Q}_2^2+\frac{{\lambda }_2}{2}{H_{l}P}_2^2+\frac{\rho }{2}\left({B-{\mathrm{PeE}}_1+V}_2^2\right)$

The present disclosure splits a solution of the augmented Lagrangian function of the time transformation P into sub-problems about variables B, P, V.

The sub-problem of the variable B of the time transformation P is as follows:

$B\in \arg \underset{B}{\min }ℒ\left(B,P,V\right)=\arg \underset{B}{\min }\frac{1}{2}{M_2-{\mathrm{RBA}}_2}_2^2+\frac{1}{2}{D-R\left({\mathrm{BA}}_2-Y_1\right)}_2^2+\frac{1}{2}{T-\left({\mathrm{BA}}_2-Y_1\right)}_2^2+\frac{\rho }{2}\left({B-{\mathrm{PeE}}_1+V}_2^2\right)$

A solution to the above equation is:

$\left(2R^{T}R+1\right){\mathrm{BA}}_2{A}_2^T+\rho B=R^T\left(M_2+D+{\mathrm{RE}}_1{A}_1\right)A_2^T+\left(T+E_1{A}_1\right)A_2^T+\rho \left({\mathrm{PeE}}_1-V\right)$

The above equation forms a Silvestre equation, which can be solved using command X=dlyap(A,B,C) in matlab.

The sub-problem of the variable P of the time transformation P is as follows:

$P\in \arg \underset{P}{\min }ℒ\left(B,P,V\right)=\arg \underset{P}{\min }\frac{{\lambda }_1}{2}{P-1_L{1}_Q}_2^2+\frac{{\lambda }_2}{2}{H_{l}P}_2^2+\frac{\rho }{2}\left({B-{\mathrm{PeE}}_1+V}_2^2\right)$

The above equation is rewritten into a form of a vector for solving:

$P^{*}=\mathrm{invec}({\lambda }_{1}1+$ $\rho \mathrm{vec}\left(\left(E_1\mathrm{eB}+E_1\mathrm{eV}\right)\right)/\left({\lambda }_{1}I+{\lambda }_2{\mathrm{rep}}_Q\left(H_L^T{H}_L\right)+\rho \mathrm{diag}\left(\mathrm{vec}\left(E_1{\mathrm{eE}}_1\right)\right)\right)$

An update rule for a dual variable of the time transformation P is as follows:

$V\leftarrow V+B-{\mathrm{PeE}}_1$

A sub-problem of the time variation image T of the target image is as follows:

$T=\arg \underset{T}{\min }\frac{1}{2}{\text{?}-D}_2^2+\frac{1}{2}{T-\left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)}_2^2+{\lambda }_T{𝒯}_{\mathrm{TNN}}$ $\text{?}\text{indicates text missing or illegible when filed}$

Auxiliary variable = is introduced, and an augmented Lagrangian function for the sub-problem of the time variation image T is as follows:

$T=\arg \underset{T}{\min }\frac{1}{2}{\mathrm{RT}-D}_2^2+\frac{1}{2}{T-\left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)}_2^2+{\lambda }_T{ℬ}_{\mathrm{TNN}}$

The present disclosure splits a solution of a Lagrangian function of the time variation image T into sub-problems about variables , , , where the , , respectively represent third-order tensor forms of matrices T, B, V.

The sub-problem of the variable of the time variation image T is as follows:

$\begin{array}{c}𝒯\in \arg \underset{ℬ}{\min }ℒ\left(𝒯,ℬ,𝒱\right)\\ =\arg \underset{𝒯}{\min }\frac{1}{2}{\mathrm{RT}-D}_2^2+\frac{1}{2}{T-\left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)}_2^2+\frac{\rho }{2}\left({ℬ-𝒯+𝒱}_2^2\right)\end{array}$

The , , are expanded into matrices, and ∂/∂T is set to be equal to 0 to obtain an equation for solving the T:

$T^{*}={\left(\frac{1}{1+\rho }{R}^T\mathrm{RT}+1\right)}^{-1}\left(\frac{1}{1+\rho }\left(R^{T}D+\eta \left(\left({\mathrm{PeE}}_1\right)A_2-E_1{A}_1\right)+\rho \left(B+V\right)\right)\right)$

The sub-problem of the variable of the time variation image T is as follows:

$ℬ\in \arg \underset{ℬ}{\min }ℒ\left(𝒯,ℬ,𝒱\right)=\arg \underset{ℬ}{\min }\frac{\rho }{2}{ℬ-𝒯+𝒱}_2^2+{\lambda }_T{ℬ}_{\mathrm{TNN}}$

A solution to the above equation is:

$B^{*}=\sum _{i=1}^3{\mathrm{fold}}_i\left(S_{\frac{\rho }{{\lambda }_T},\epsilon }\left(B+V\right)\right),\mathrm{where}{S}_{\frac{\rho }{{\lambda }_T},\epsilon }(·)=U\left(\Sigma -\frac{\rho }{{\lambda }_T}\mathrm{diag}\left(\epsilon \right)\right)V^T$

represents a singular value contraction operator, and UΣV^T represents a singular value decomposition process.

An update rule for a dual variable of the time transformation T is as follows:

$𝒱\leftarrow 𝒱+ℬ-𝒯$

Step 7.3: An alternating direction method of multipliers is used to perform optimal solving on the augmented Lagrangian function to obtain the abundance matrix A₂, the time variation matrix P and the time variation image T of the target image.

Step 8: The target image is finally obtained by multiplying the endmember matrix E₁ at the time point T₁ by the time variation matrix P at the time point T₂ element by element and then by the abundance matrix A₂ with high resolution: X=(PeE₁) A₂.

Embodiment 2

A device for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions includes:

• • a first obtaining module configured to obtain hyperspectral image H₁ and multispectral image M₁ at time point T₁, and multispectral image M₂ at time point T₂, and perform preprocessing to obtain upsampled hyperspectral image Ĥ₁;
  • a search module configured to: based on a local linear constraint, search for adjacent similar full-band block M₁(Ω_i,j,k) in the multispectral image M₁ to represent an image block centered at (i, j) in the multispectral image M₂, and obtain a representation weight function of the similar full-band block; and share the representation weight function with the upsampled hyperspectral image Ĥ₁ to obtain upsampled hyperspectral image Ĥ₂^*∈ at a prediction time point;
  • a second obtaining module configured to obtain a time variation image represented as T=Ĥ₂^*−Ĥ₁ for an initialized target fused image;
  • a decomposition module configured to decompose target image X₂ at the time point T₂ into endmember matrix E₂ and abundance matrix A₂ based on a spectral linear mixed model, where X₂=E₂A₂, and the target image X₂ is an image with high spatial, temporal and spectral resolutions;
  • an establishment module configured to establish observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions;
  • a representation module configured to introduce a generalized linear mixed model, and represent the endmember matrix E₂ as E₂=PeE₁, where the symbol ⊙ represents a Hadamard product of element-by-element multiplication, P represents a time variation matrix of the target image, and E₁ represents an endmember matrix at the time point T₁;
  • a third obtaining module configured to obtain the abundance matrix A₂, the time variation matrix P, and time variation image T of the target image; and
  • an obtaining module configured to obtain the target image finally by multiplying the endmember matrix E₁ at the time point T₁ by the time variation matrix P at the time point T₂ element by element and then by the abundance matrix A₂ with high resolution: X=(PeE₁)A₂.

In conclusion, the present disclosure proposes a time variation model for depicting and modeling a time variation characteristic of an image, and a resolution enhancement reconstruction model that is of a time variation image and achieves a fidelity optimization solving of a time variation. The time variation model and the resolution enhancement reconstruction model of the time variation image are coupled to establish a generalized spatial-temporal-spectral fusion framework. The time variation model models and depicts time variation characteristics between images, while the resolution enhancement reconstruction model of the time variation image obtains a fidelity constraint on the time variation. This method can achieve spatial-temporal-spectral fusion of a hyperspectral image with low temporal and spatial resolutions and a multispectral image with high temporal and spatial resolutions to generate a fused image with high spatial and spectral fidelity and high spatial, temporal and spectral resolutions.

Claims

1. A method for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions, comprising: step 1: obtaining a hyperspectral image H1 and a multispectral image M1 at a time point T1, and a multispectral image M2 at a time point T2, and performing preprocessing to obtain an upsampled hyperspectral image Ĥ1;step 2: based on a local linear constraint, searching for an adjacent similar full-band block M1(Ωi,j,k) in the multispectral image M1 to represent an image block centered at (i,j) in the multispectral image M2, and obtaining a representation weight function of the adjacent similar full-band block; and sharing the representation weight function with the upsampled hyperspectral image Ĥ1 to obtain an upsampled hyperspectral image Ĥ2*∈ at a prediction time point;step 3: obtaining a time variation image represented as T=Ĥ2*−Ĥ1 for an initialized target fused image;step 4: decomposing a target image X2 at the time point T2 into an endmember matrix E2 and an abundance matrix A2 based on a spectral linear mixed model, wherein X2=E2A2, and the target image X2 is an image with high spatial, temporal and spectral resolutions;step 5: establishing observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions;step 6: introducing a generalized linear mixed model, and representing the endmember matrix E2as E2=PeE1, wherein a symbol e represents a Hadamard product of element-by-element multiplication, P represents a time variation matrix of the target image, and E1 represents an endmember matrix at the time point T1;step 7: obtaining the abundance matrix A2, the time variation matrix P and a time variation image T of the target image; andstep 8: obtaining the target image finally by multiplying the endmember matrix E1 at the time point T1 by the time variation matrix P at the time point T2 element by element and then by the abundance matrix A2 with high resolution: X=(PeE1)A2.

2. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 1, wherein the step 1 comprises: step 1.1: obtaining the hyperspectral image Hi∈ and the multispectral image M1∈ at the time point T1, and the multispectral image M2∈ at the time point T2, wherein l and L represent quantities of bands, w and W represent quantities of pixels, l<L, and w<W;step 1.2: normalizing the hyperspectral image H1 and the multispectral images M1 and M2, and upsampling a normalized hyperspectral image H1 to a spatial size of the multispectral image to obtain the Ĥ1∈; andstep 1.3: folding the multispectral images M1 and M2 and the upsampled hyperspectral image Ĥ1 into a three-dimensional form along a spectral dimension, and splitting the three-dimensional form into full-band blocks with both a width and a height being √{square root over (c)} and the spectral dimension being l.

3. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 2, wherein in the step 5, the observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions are represented as follows: Mi=RXi=REiAi Hi=XiFS=EiAiFSwherein i∈{1,2}; Xi∈, Hi∈, and Mi∈ respectively represent an image with high spatial, temporal and spectral resolutions, a hyperspectral image, and a multispectral image at a time point Ti; Ei∈ and Ai∈ respectively represent an endmember matrix and an abundance matrix of the image with high spatial, temporal and spectral resolutions at the time point Ti, and Q represents a quantity of endmembers; and R∈ represents a spectral downsampling matrix, and F∈ and S∈ respectively represent a spatial blur matrix and a spatial downsampling matrix.

4. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 3, wherein the step 7 comprises: step 7.1: separately introducing a regularization term of the abundance matrix A2, a regularization term of the time variation matrix P and a regularization term of the time variation image T of the target image to obtain an energy function for a fusion problem;step 7.2: splitting the energy function into three sub-functions corresponding to variables of the abundance matrix A2, the time variation matrix P and the time variation image T, introducing a plurality of auxiliary variables for each sub-function to split the sub-function into sub-problems, and constructing an augmented Lagrangian function for each sub-problem; andstep 7.3: using an alternating direction method of multipliers to perform optimal solving on the augmented Lagrangian function to obtain the abundance matrix A2, the time variation matrix P and the time variation image T of the target image.

5. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 4, wherein in the step 7.1, the energy function is represented as follows:

6. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 5, wherein the step 7.1 comprises: step 7.1.1: introducing the regularization term of the abundance matrix A2 of the target image, wherein the regularization term of the abundance matrix A2 of the target image is a spatial smoothing constraint and represented as follows:

7. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 6, wherein in the step 7.2, a sub-problem of the abundance matrix A2 of the target image is as follows:

8. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 7, wherein in the step 7.2, a sub-problem of the time variation matrix P of the target image is as follows:

9. The method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 8, wherein in the step 2, the image block centered at the (i,j) in the multispectral image M2 is represented using the adjacent similar full-band block M1(Ωi,j,k) in the multispectral image

10. A device for generating a satellite remote-sensing image with high spatial, temporal and spectral resolutions, configured to execute the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 1, and comprising: a first obtaining module configured to obtain the hyperspectral image H1 and the multispectral image M1 at the time point T1, and the multispectral image M2 at the time point T2, and perform preprocessing to obtain the upsampled hyperspectral image Ĥ1;a search module configured to: based on the local linear constraint, search for the adjacent similar full-band block M1(Ωi,j,k) in the multispectral image M1 to represent the image block centered at (i,j) in the multispectral image M2, and obtain the representation weight function of the adjacent similar full-band block; and share the representation weight function with the upsampled hyperspectral image Ĥ1 to obtain the upsampled hyperspectral image Ĥ2*∈ at the prediction time point;a second obtaining module configured to obtain the time variation image represented as T=Ĥ2*−Ĥ1 for the initialized target fused image;a decomposition module configured to decompose the target image X2 at the time point T2 into the endmember matrix E2 and the abundance matrix A2 based on the spectral linear mixed model, wherein X2=E2A2, and the target image X2 is the image with high spatial, temporal and spectral resolutions;an establishment module configured to establish the observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions;a representation module configured to introduce the generalized linear mixed model, and represent the endmember matrix E2 as E2=PeE1, wherein the symbol e represents the Hadamard product of element-by-element multiplication, P represents the time variation matrix of the target image, and E1 represents the endmember matrix at the time point T1;a third obtaining module configured to obtain the abundance matrix A2, the time variation matrix P and the time variation image T of the target image; andan obtaining module configured to obtain the target image finally by multiplying the endmember matrix E1 at the time point Ti by the time variation matrix P at the time point T2 element by element and then by the abundance matrix A2 with high resolution: X=(PeE1)A2.

11. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 10, wherein the step 1 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions comprises: step 1.1: obtaining the hyperspectral image Hi∈ and the multispectral image M1∈ at the time point T1, and the multispectral image M2∈ at the time point T2, wherein l and L represent quantities of bands, w and W represent quantities of pixels, l<L, and w<W;step 1.2: normalizing the hyperspectral image H1 and the multispectral images M1 and M2, and upsampling a normalized hyperspectral image H1 to a spatial size of the multispectral image to obtain the Ĥ1∈; andstep 1.3: folding the multispectral images M1 and M2 and the upsampled hyperspectral image Ĥ1 into a three-dimensional form along a spectral dimension, and splitting the three-dimensional form into full-band blocks with both a width and a height being √{square root over (c)} and the spectral dimension being l.

12. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 11, wherein in the step 5 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions, the observation models of the multispectral image, the hyperspectral image, and the image with high spatial, temporal and spectral resolutions are represented as follows: Mi=RXi=REiAi Hi=XiFS=EiAiFSwherein i∈{1,2}; Xi∈, Hi∈, and Mi∈ respectively represent an image with high spatial, temporal and spectral resolutions, a hyperspectral image, and a multispectral image at a time point Ti; Ei∈ and Ai∈ respectively represent an endmember matrix and an abundance matrix of the image with high spatial, temporal and spectral resolutions at the time point Ti, and Q represents a quantity of endmembers; and R∈ represents a spectral downsampling matrix, and F∈ and S∈ respectively represent a spatial blur matrix and a spatial downsampling matrix.

13. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 12, wherein the step 7 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions comprises: step 7.1: separately introducing a regularization term of the abundance matrix A2, a regularization term of the time variation matrix P and a regularization term of the time variation image T of the target image to obtain an energy function for a fusion problem;step 7.2: splitting the energy function into three sub-functions corresponding to variables of the abundance matrix A2, the time variation matrix P and the time variation image T, introducing a plurality of auxiliary variables for each sub-function to split the sub-function into sub-problems, and constructing an augmented Lagrangian function for each sub-problem; andstep 7.3: using an alternating direction method of multipliers to perform optimal solving on the augmented Lagrangian function to obtain the abundance matrix A2, the time variation matrix P and the time variation image T of the target image.

14. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 13, wherein in the step 7.1 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions, the energy function is represented as follows:

15. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 14, wherein the step 7.1 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions comprises: step 7.1.1: introducing the regularization term of the abundance matrix A2 of the target image, wherein the regularization term of the abundance matrix A2 of the target image is a spatial smoothing constraint and represented as follows:

16. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 15, wherein in the step 7.2 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions, a sub-problem of the abundance matrix A2 of the target image is as follows:

17. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 16, wherein in the step 7.2 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions, a sub-problem of the time variation matrix P of the target image is as follows:

18. The device for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions according to claim 17, wherein in the step 2 of the method for generating the satellite remote-sensing image with high spatial, temporal and spectral resolutions, the image block centered at the (i,j) in the multispectral image M2 is represented using the adjacent similar full-band block M1(Ωi,j,k) in the multispectral image