Sparse Bayesian remote sensing image reconstruction method and system based on Laplace level prior and GAMP

Abstract

The invention discloses a sparse Bayesian remote sensing image reconstruction method and system based on Laplace level prior and GAMP, and the method comprises the following steps: S1, selecting a remote sensing image data set, setting image pixels, converting a remote sensing image into an amplitude diagram, and carrying out the sparse processing of the amplitude diagram; s2, establishing a sparse Bayesian model; s3, calculating a likelihood function, Laplacian hierarchical prior distribution of sparse signals and a marginal likelihood function; s4, performing Bayesian inference according to a maximum likelihood method, and calculating posterior distribution of sparse signals; s5, obtaining a mean value and a variance of a sparse signal and an updating formula of a hyper-parameter through an expectation maximization algorithm; s6, introducing a GAMP algorithm to update the mean value and the variance; s7, adding a damping factor to the GAMP structure; s8, judging whether an iteration termination condition is met or not, if so, outputting posteriori estimation of the sparse signal, and otherwise, executing S5; and S9, performing inverse wavelet transform on the output sparse signal to obtain a reconstructed image. The method is low in reconstruction error and high in reconstruction precision and efficiency.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing image reconstruction, and in particular relates to an efficient sparse Bayesian remote sensing image reconstruction method and system based on Laplace hierarchical prior and GAMP. Background Art

[0002] With the continuous development of remote sensing technology, high-resolution remote sensing images are increasingly used in environmental monitoring, urban planning, military reconnaissance, transportation and other fields. Remote sensing technology refers to the technology of obtaining information on the earth's surface from a long distance through various sensors. With the advancement of sensor technology, satellite technology and computer technology, the acquisition of remote sensing data has become more efficient and accurate. With the continuous improvement of the degree of digitization, the demand for high-performance remote sensing image data processing technology is increasing. For these high-dimensional data analysis problems containing a large number of features, algorithms such as dimensionality reduction, denoising, compressed sensing (CS), and sparse representation (SR) are currently used to preprocess and optimize the data and remove redundant information from the feature data. In this way, the problem being processed is simplified without losing feature information or affecting performance, so as to improve data processing efficiency.

[0003] Compressed sensing technology reconstructs sparse signals from a small number of non-adaptive linear measurements. A major problem in compressed sensing is sparse regularization, which can be solved by methods such as Bregman splitting, orthogonal matching pursuit (OMP), Bayesian estimation, iterative threshold or compressed sampling matching pursuit (CoSaMP) algorithm. However, the computational complexity of remote sensing image reconstruction technology based on compressed sensing is high at this stage, involving iterative operations and high-order matrix operations, which may not meet the needs of fast response when processing large-scale data. And because compressed sensing technology relies on the sparsity of the signal, if the signal itself is not sparse or the sparsity is not obvious, the reconstruction quality may be affected. In addition, the presence of noise may also cause the reconstruction result to be unstable. In contrast, sparse representation SR technology obtains highly sparse signals by using as few non-zero coefficients as possible to represent the main information of the signal, providing good feature data for further image processing. As the inverse problem of SR, sparse reconstruction has received widespread attention in recent years. Sparse Bayesian Learning (SBL) is one of the most commonly used methods in sparse reconstruction. Based on the Bayesian idea, it performs probability modeling on the unknown parameters in the model, and then obtains a sparse solution through Bayesian reasoning. And the estimated value of the unknown parameters is solved through the expectation maximization (EM) algorithm. Compared with other sparse reconstruction models, the SBL algorithm has more inference advantages.

[0004] However, the current remote sensing image reconstruction algorithm based on the SBL model has the following defects: the EM algorithm is used to iteratively update the parameters, which involves the inversion of high-dimensional matrices and has high computational complexity; when processing large-scale data or images with complex backgrounds, the iterative convergence speed will slow down. The Generalized Approximate Message Passing (GAMP) algorithm proposed in recent years can effectively approximate the mean and variance required to be solved in the E step of the EM algorithm in the SBL model by transforming the estimation problem, using Gaussian approximation and the central limit theorem, and effectively avoids the inversion calculation of the matrix in the solution of the model, greatly improving the reconstruction efficiency. However, in the SBL model based on GAMP, multiple random variables may present a thick-tailed distribution, and the prior distribution of the Gaussian mixture function in GAMP cannot be well fitted. Therefore, it may lead to local optimization during iteration. This means that more time is needed to approximate the hyperparameters in the M step of the EM algorithm. The Laplace distribution has good robustness to thick-tailed distributions. It can easily extract samples far from the mean (abnormal samples). Therefore, the Laplace prior is more suitable for the reconstruction of natural images. Moreover, the Laplace hierarchical prior (LHP) can improve the sparsity of the unknown signal and obtain a lower reconstruction error. Summary of the invention

[0005] In view of the above situation, the present invention proposes an efficient sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP, called LHP-GAMP-SBL method, which can effectively improve the reconstruction accuracy of remote sensing images. The present invention integrates the GAMP framework and Laplace hierarchical prior (LHP) into the SBL framework. First, the GAMP algorithm is introduced on the basis of SBL to replace the E step of the classical EM algorithm, avoiding the matrix inversion of solving the mean of the unknown signal, and greatly reducing the computational complexity of the model. The Laplace hierarchical prior is used to implement the M step in SBL, which has good robustness for the fat-tail phenomenon in natural signals, and the introduction of the damping factor enhances the convergence of the model. The integration of the two structures further deepens the hierarchical structure of the model, enhances the representation ability of the model, and improves the reconstruction efficiency and accuracy. The present invention improves the image reconstruction accuracy and convergence speed by using the GAMP algorithm and the Laplace hierarchical prior. Compared with the existing image reconstruction technology, the image reconstruction using the present invention has significant improvements in reconstruction accuracy and efficiency, and can well meet the application requirements in actual scenes.

[0006] The present invention adopts the following technical scheme:

[0007] An efficient sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP includes the following steps:

[0008] S1, obtain remote sensing image data set, set the pixel size of all remote sensing images, convert remote sensing images into amplitude maps and perform sparse preprocessing on each amplitude map;

[0009] S2, establish a sparse Bayesian model;

[0010] S3, calculate the likelihood function, the Laplace stratified prior distribution of the sparse signal x, and the marginal likelihood function;

[0011] S4, calculating the posterior distribution of the sparse signal x by Bayes' theorem according to the likelihood function, Laplace stratified prior distribution and marginal likelihood function calculated in step S3;

[0012] S5. Obtain the posterior mean and posterior variance of the sparse signal and the update formula of the hyperparameters through the expectation maximization (EM) algorithm;

[0013] S6, introduce the GAMP algorithm to replace the E step of the expectation-maximization algorithm to update the mean and variance in a low-complexity way;

[0014] S7, adding damping factor to GAMP structure to enhance convergence;

[0015] S8, judging whether the iteration termination condition is met, if so, outputting the estimation of the sparse signal; if not, executing step S5 until the iteration termination condition is met;

[0016] S9. Perform inverse wavelet transform on the output sparse signal x to obtain a reconstructed remote sensing image.

[0017] The reconstruction performance can be evaluated by four indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM) and peak signal-to-noise ratio (PSNR).

[0018] Preferably, in step S1, all original images are adjusted to 256×256 pixels, and multi-band remote sensing images are converted into grayscale images and processed into .mat format files for storing feature data of the images.

[0019] Preferably, in step S1, the multi-band image is converted into an amplitude map of a remote sensing image using principal component analysis (PCA); finally, each amplitude map is subjected to sparse preprocessing by adding a "symmlet,8" wavelet transform with a coarsest scale of 3 and a finest scale of 7, and the image is reconstructed using a multi-scale compressed sensing (CS) scheme. The wavelet coefficients to be estimated are set to: n = 16320, and the number of compressed samples is set to: m = 10608.

[0020] Preferably, in step S2, a sparse Bayesian (SBL) mathematical model is established: y = Dx + e. is the observation vector, is a sparse signal, is a known dictionary matrix, is the noise vector. According to Bayes’ theorem, the estimate of x is obtained by maximizing the posterior distribution p(x,ω,λ,θ|y):

[0021]

[0022] Among them, ω, λ, and θ represent hyperparameters, and p(x, ω, λ, θ|y) is a conditional probability distribution, which is represented as the posterior distribution in the SBL model.

[0023] Preferably, in step S3, the likelihood function p(y|x,θ) is calculated according to the posterior distribution expression in step S2:

[0024] Assume that the noise e has a mean of 0 and a variance of θ -1 Gaussian distribution of:

[0025]

[0026] The hyperparameter θ follows the Gamma prior:

[0027]

[0028] Among them, shape parameter a>0, scale parameter b>0, Γ(a θ ) is the Gamma function, and the expectation of θ is expressed as:

[0029] Preferably, step S3 calculates the Laplace hierarchical prior distribution p(x|λ) of the sparse signal x according to the posterior distribution expression in step S2:

[0030] Using a hierarchical model prior, the first layer:

[0031]

[0032] Where ω=(ω1,ω2,...,ω N ),ω i is control x i The sparsity hyperparameter, i = 1, 2, ..., N; the second layer hyperparameter ω i Modeling:

[0033]

[0034] Among them, ω i ≥0, λ≥0. Therefore, the Laplace stratified prior distribution of the sparse signal x is:

[0035]

[0036] Among them, the parameter λ obeys the Gamma super prior, that is:

[0037] Preferably, in step S3, according to the posterior distribution expression in step S2, the marginal likelihood function p(y) is calculated:

[0038]

[0039] Where C = (θ -1 I+DΛ -1 D T ), I represents the identity matrix, is a diagonal matrix.

[0040] Preferably, in step S4, according to Bayes' theorem, the convolution result of two Gaussian functions is still a Gaussian function, so it can be calculated that the posterior density p(x, ω, λ, θ|y) of x also obeys the Gaussian distribution:

[0041]

[0042] Preferably, in step S5, the mean μ of the sparse signal x is obtained by using the expectation maximization algorithm according to the posterior distribution p(x, ω, λ, θ|y) obtained in step S4. x and variance∑ x And the update formula of the hyperparameters (ω, λ, θ) is as follows:

[0043] μ x =θ∑ x D T y

[0044] ∑ x =(θD T D+Λ) -1

[0045]

[0046] The hyperparameter γ is updated as follows:

[0047]

[0048] in, is the digamma function, right The derivative of . When the hyperparameter γ is estimated by the above formula, the reconstruction performance is slightly worse than when γ is zero, so γ is fixed to 0.

[0049] Preferably, in step S6, a GAMP algorithm is introduced to replace the E step of the expectation-maximization algorithm to update the mean and variance in a low-complexity manner.

[0050] For a given prior distribution information p(x|ω) and marginal likelihood function p(y), GAMP constructs an iterative expression through quadratic approximation and Taylor series expansion; through the sum-product algorithm (SPA), the input function And the output function It can be defined as:

[0051]

[0052] Among them, m∈M represents the mth atom of the input function, n∈N represents the nth atom of the output function, and the intermediate variable and denote the approximate estimate of the sparse signal x and the noise-free model k=Dx, respectively. and Corresponding to the variance of x and k respectively. Therefore, the update of the input function and the output function can be obtained as:

[0053]

[0054] The iterative update of the mean and variance is expressed as:

[0055] Preferably, in step S7, a damping factor is added to the GAMP structure to enhance the convergence of the model, and the damping factor is calculated by the following formula:

[0056]

[0057] Among them, Ω max represents the maximum singular value of dictionary D; Σ Ω represents the sum of the singular values ​​of D; α and β are constant terms; damping factor and is the root of the above formula, and Used to control Used to control x.

[0058] Preferably, step S8 determines whether x satisfies the iterative convergence condition. If so, output an estimate of x: x = μ x t ; If not, execute step S5. The iterative convergence condition is expressed as:

[0059]

[0060] Preferably, in step S9, the output x is subjected to an inverse wavelet transform to obtain a reconstructed image. The calculation formula for the reconstruction error is: in, represents the reconstructed image, and S represents the real image.

[0061] Finally, the reconstruction performance can be evaluated by four indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM) and peak signal-to-noise ratio (PSNR).

[0062] The present invention also discloses a sparse Bayesian remote sensing image reconstruction system based on Laplace hierarchical prior and GAMP, which includes the following modules based on the above method:

[0063] Remote sensing image acquisition and preprocessing module: obtain remote sensing image data sets, set the pixels of remote sensing images, convert remote sensing images into amplitude maps and perform sparse processing on the amplitude maps;

[0064] Sparse Bayesian model building module: build sparse Bayesian model;

[0065] Computation module: calculate likelihood function, Laplace stratified prior distribution of sparse signals, marginal likelihood function;

[0066] Posterior distribution calculation module: Based on the likelihood function, Laplace stratified prior distribution and marginal likelihood function obtained by the calculation module, the posterior distribution of the sparse signal is calculated by Bayesian theorem;

[0067] Update module: obtain the posterior mean and posterior variance of the sparse signal and the update formula of the hyperparameters through the expectation maximization algorithm;

[0068] Mean and variance update module: introduce the GAMP algorithm to update the mean and variance;

[0069] Damping factor adding module: adding damping factor to GAMP structure;

[0070] Iteration module: determines whether the iteration termination condition is met. If so, it outputs the estimation of the sparse signal.

[0071] Inverse wavelet transform module: Perform inverse wavelet transform on the output sparse signal to obtain the reconstructed remote sensing image.

[0072] The beneficial effects of the present invention are:

[0073] The present invention introduces a generalized approximate message passing algorithm (GAMP) to approximate the mean and variance required to be solved in the E step of the SBL model. In solving the model, the inversion calculation of the matrix is ​​effectively avoided, and the reconstruction efficiency is greatly improved. The present invention improves the sparsity of the model through the Laplace hierarchical prior, and can well handle the thick tail phenomenon. The convergence of the model is improved by adopting an appropriate damping coefficient. The integration of the two structures further deepens the hierarchical structure of the model, improves the robustness of the model and the reconstruction accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0075] Figure 1 It is a model structure diagram of a preferred embodiment of the present invention, which is an efficient sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP.

[0076] Figure 2 It is an algorithm flow chart of an efficient sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP in a preferred embodiment of the present invention.

[0077] Figure 3 It is a reconstructed image sample of the remote sensing image data set adopted by the present invention.

[0078] in, Figure 3 (a), 3(b) and 3(c) represent two example images from the DIOR, DOTA and UCM datasets, respectively. From left to right, they represent the original image and the reconstruction effects of the six comparison methods. The values ​​in brackets represent the reconstruction error (ERROR).

[0079] Figure 4 It is a comparison chart of the reconstruction performance of the remote sensing image data set adopted by the present invention.

[0080] in, Figure 4 (a), 4(b), 4(c), and 4(d) respectively show the comparison graphs of reconstruction error (ERROR), average running time (TIME), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM) of 20 randomly selected images; Figure 4 (e) and 4(f) respectively represent the reconstruction error (ERROR) and the average running time (TIME) at each iteration.

[0081] Figure 5 It is a block diagram of an efficient sparse Bayesian remote sensing image reconstruction system based on Laplace hierarchical prior and GAMP in a preferred embodiment of the present invention. DETAILED DESCRIPTION

[0082] The following describes the implementation of the present invention through specific embodiments, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific implementations, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention.

[0083] The preferred embodiment of the present invention discloses an efficient sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP, referring to Figure 1 As shown in the figure, it is the model structure diagram of this method. Figure 2 , which is the algorithm flow chart of this method. Figure 3 , which is a reconstructed image sample of the remote sensing image data set used in the present invention. All images are processed according to the process of the method of this embodiment to illustrate in detail the effect of the method of the present invention on improving the accuracy of remote sensing image reconstruction.

[0084] like Figure 1-3 As shown, this embodiment provides an efficient sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP, comprising the following steps:

[0085] S1. Obtain remote sensing image datasets; set the pixel size of all remote sensing images to 256×256, convert multi-band remote sensing images into grayscale images and process them into .mat format files to save image feature data; use principal component analysis (PCA) to convert multi-band images into amplitude maps of remote sensing images; finally, add wavelet transform to each amplitude map for sparse preprocessing, and use multi-scale compressed sensing (CS) scheme for image reconstruction.

[0086] Compared with general image reconstruction technology, remote sensing image reconstruction is different from other image reconstruction technologies in the preprocessing process. In order to highlight different types of ground object information in remote sensing images and reduce the influence of factors such as atmosphere and illumination, it is necessary to extract the intensity information of multi-band remote sensing images and generate amplitude maps of remote sensing images to facilitate subsequent image reconstruction. Therefore, this embodiment mainly generates amplitude maps of remote sensing images through the PCA method.

[0087] PCA is a linear transformation method that is simple to calculate, easy to implement, and does not require complex parameter adjustments. PCA can compress the information in multi-band remote sensing images into fewer principal components, thereby reducing the dimension of the data, helping to reduce storage requirements and computational complexity. As the number of principal components increases, the proportion of noise contained usually increases. Therefore, by selecting the first few principal components, noise can be effectively suppressed, the quality of the image can be improved, and subsequent classification tasks can be facilitated.

[0088] S2. Establish a sparse Bayesian (SBL) mathematical model: y=Dx+e.

[0089] in, is the measurement vector, is a sparse signal, is a known dictionary matrix, is the noise vector. According to Bayes’ theorem, the estimate of x is obtained by maximizing the posterior distribution p(x,ω,λ,θ|y):

[0090]

[0091] S3. Calculate the likelihood function according to the posterior distribution expression in step S2:

[0092]

[0093] Among them, the noise e has a mean of 0 and a variance of θ -1 The Gaussian distribution of . The hyperparameter θ follows the Gamma prior:

[0094]

[0095] Among them, shape parameter a>0, scale parameter b>0, Γ(a θ ) is the Gamma function, and the expectation of θ is expressed as:

[0096] According to the posterior distribution expression in step S2, calculate the prior distribution:

[0097]

[0098] in, ω=(ω1,ω2,...,ω N ),ω i is control x i The sparsity hyperparameter, i = 1, 2, ..., N. For the hyperparameter ω i Modeling:

[0099]

[0100] Among them, ω i ≥0, λ≥0. The parameter λ obeys the Gamma hyper-prior:

[0101] According to the posterior distribution expression in step S2, the marginal likelihood function is calculated:

[0102]

[0103] Where C = (θ-1 I+DΛ -1 D T ), I represents the identity matrix, is a diagonal matrix.

[0104] S4. According to Bayes' theorem, we can calculate that the posterior density of x, p(x, ω, λ, θ|y), also obeys the Gaussian distribution:

[0105]

[0106] S5. According to the posterior distribution p(x,ω,λ,θ|y) obtained in step S4, the mean μ of the sparse signal x is obtained by the expectation maximization algorithm. x and variance Σ x And the update formula of the hyperparameters (ω, λ, θ) is as follows:

[0107] μ x =θΣ x D T y

[0108] ∑ x =(θD T D+Λ) -1

[0109]

[0110]

[0111] The hyperparameter γ is updated as follows:

[0112]

[0113] in, is the digamma function, right The derivative of .

[0114] S6. The GAMP algorithm is introduced to replace the E step of the expectation-maximization algorithm to update the mean and variance in a low-complexity way.

[0115] For a given prior distribution information p(x|ω) and marginal likelihood function p(y), GAMP constructs an iterative expression through quadratic approximation and Taylor series expansion. During the message propagation process, the function is input through the sum-product algorithm (SPA) And the output function It can be defined as:

[0116]

[0117] Among them, m∈M represents the mth atom of the input function, n∈N represents the nth atom of the output function, and the intermediate variable and They represent the approximate estimates of the sparse signal x and the noise-free model k=Dx, and their variances are and Therefore, the update of the input function and the output function can be obtained as:

[0118]

[0119] The iterative update of the mean and variance is expressed as:

[0120] S7. Adding a damping factor to the GAMP structure enhances the convergence of the model. The formula for solving the damping factor is as follows:

[0121]

[0122] The GAMP structure represents a model calculated by the GAMP algorithm. The GAMP structure is introduced into the SBL model to replace the E step of the EM algorithm to update the posterior mean and posterior covariance.

[0123] In the above formula, Ω max represents the maximum singular value of dictionary D; Σ Ω represents the sum of the singular values ​​of D; α and β are constant terms; damping factor and is the root of the above formula, and Used to control Used to control x.

[0124] S8. Determine whether x satisfies the iterative convergence condition. If so, output the estimate of x: x = μ x t ; If not, execute step S5. The iteration termination condition is:

[0125]

[0126] S9. Perform inverse wavelet transform on the output x to obtain a reconstructed image. The calculation formula for the reconstruction error is: in, represents the reconstructed image, and S represents the real image.

[0127] Finally, the reconstruction performance is evaluated by four indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM) and peak signal-to-noise ratio (PSNR).

[0128] The following describes in detail each step of the above embodiment in conjunction with specific data.

[0129] In step S1, the remote sensing data sets obtained are specifically: DIOR optical remote sensing image data set, DOTA optical remote sensing image data set and UC Merced Land Use (UCM) remote sensing image data set. 50 images (a total of 150 images) are randomly selected from each data set for the experiment. In order to highlight the different types of ground object information in the remote sensing image and reduce the influence of factors such as atmosphere and illumination, the intensity information of each band is extracted in the experiment to generate the amplitude map of the remote sensing image, which is convenient for subsequent image reconstruction.

[0130] In step S1, the pixel size of all remote sensing images is set to 256×256, 50 images are randomly selected from each data set and numbered in sequence. The multi-band remote sensing images are converted into grayscale images and read into .mat files to save the feature data of the images.

[0131] Step S1, first use the principal component analysis (PCA) method to convert the multi-band image into a set of unrelated components, where the first principal component usually contains the most important information in the image, that is, the amplitude map of the remote sensing image. Then add the "symmlet,8" wavelet transform with the coarsest scale of 3 and the finest scale of 7 to each image for sparse preprocessing, and use the multi-scale compressed sensing (CS) scheme for image reconstruction. The number of compressed samples is set to 10672, and the band estimated wavelet coefficients are set to 16320.

[0132] Step S2, based on the sparse representation theory (SR) and Bayesian theorem, establish a sparse Bayesian (SBL) mathematical model: y = Dx + e. Where: is the measurement vector, is a sparse signal, is a known dictionary matrix, is the noise vector. According to Bayes’ theorem, the estimate of x is obtained by maximizing the posterior distribution p(x,ω,λ,θ|y).

[0133]

[0134] Step S3, calculate the likelihood function p(y|x,θ) according to the posterior distribution expression in step S2. Assume that the noise e has a mean of 0 and a variance of θ -1 Gaussian distribution of:

[0135]

[0136] The hyperparameter θ follows the Gamma prior:

[0137]

[0138] Among them, shape parameter a>0, scale parameter b>0, is the Gamma function, and the expectation of θ is expressed as:

[0139] According to the posterior distribution expression in step S2, calculate the Laplace hierarchical prior distribution p(x|λ) of the sparse signal x. Assume that each element in the sparse coefficient x is independent and identically distributed. Add a prior density to x to achieve sparsity. According to the SBL theory, a hierarchical model prior is used, the first layer:

[0140]

[0141] Where ω=(ω1,ω2,...,ω N ), is a diagonal matrix, ω i is control x i The sparsity hyperparameter, i = 1, 2, ..., N; the second layer hyperparameter ω i Modeling:

[0142]

[0143] Among them, ω i ≥0, λ≥0. Therefore, the Laplace stratified prior distribution of the sparse signal x is obtained as:

[0144]

[0145] Among them, the parameter λ obeys the Gamma hyper-prior:

[0146] According to the posterior distribution expression in step S2, the marginal likelihood function p(y) is calculated. The marginal likelihood function is calculated by the total probability formula:

[0147]

[0148] Where C = (θ -1 I+DΛ -1 D T ), I represents the identity matrix.

[0149] Step S4, according to Bayes’ theorem, the convolution result of two Gaussian functions is still a Gaussian function, so it can be calculated that the posterior density of x p(x, ω, λ, θ|y) also obeys the Gaussian distribution:

[0150]

[0151] Step S5, according to the posterior distribution p(x,ω,λ,θ|y) obtained in step S4, the mean μ of the sparse signal x is obtained by the expectation maximization algorithm x and variance∑ x And the update formula of the hyperparameters (ω, λ, θ) is as follows:

[0152] μ x =θΣ x D T y

[0153] ∑ x =(θD T D+Λ) -1

[0154]

[0155] The hyperparameter γ is updated as follows:

[0156]

[0157] in, is the digamma function, right The derivative of ; When the hyperparameter γ is estimated by the above formula, the reconstruction performance is slightly worse than when γ is zero, so γ is fixed to 0.

[0158] Step S6, introducing the GAMP algorithm to replace the E step of the expectation-maximization algorithm, updates the mean and variance in a low-complexity manner.

[0159] For a given prior distribution information p(x|ω) and marginal likelihood function p(y), GAMP constructs an iterative expression through quadratic approximation and Taylor series expansion. The LHP-GAMP-SBL model structure is as follows: Figure 1 As shown. Figure 1 In the message propagation process, the function is input through the sum-product algorithm (SPA) And the output function It can be defined as:

[0160]

[0161] Among them, m∈M represents the mth atom of the input function, n∈N represents the nth atom of the output function, and the intermediate variable and denote the approximate estimate of the sparse signal x and the noise-free model k=Dx, respectively. and Corresponding to the variance of x and k respectively. Therefore, according to step S3, the above formula can be simplified to obtain the update of input function and output function:

[0162]

[0163] The iterative update of the mean and variance is expressed as:

[0164] Step S7, adding a damping factor to the GAMP structure to enhance the convergence of the model, the damping factor solution formula is as follows:

[0165]

[0166] Among them, Ω max represents the maximum singular value of dictionary D, Σ Ω represents the sum of the singular values ​​of D, and α and β are constant terms. Based on the above theory, it is further assumed that in the process of message propagation, there are intermediate variables propagated and variance Damping factor and is the root of the above formula, and Used to control Used to control the sparse signal x.

[0167] Step S8, determining the mean μ in step S6 x ←G x Whether the iterative convergence condition is met, if so, output the estimate of x: x = μ x t ; If not, execute step S5. The iterative convergence condition is expressed as:

[0168]

[0169] Step S9, perform inverse wavelet transform on the output x obtained in step S8 to obtain a reconstructed image. The calculation formula for the reconstruction error is: in represents the reconstructed image, and S represents the real image.

[0170] Finally, the reconstruction performance is evaluated by four indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM) and peak signal-to-noise ratio (PSNR).

[0171] Table 1 is a comparison table of the computational complexity of the E-step and M-step of different reconstruction models;

[0172] Table 1

[0173]

[0174] Table 2 is a comparison table of reconstruction performance of two indicators, reconstruction error (ERROR) and average running time (TIME) for all images;

[0175] Table 2

[0176]

[0177] * A(DIOR Dataset)B(DOTA Dataset)C(UCM Dataset)

[0178] Figure 4 This is a comparison chart of the reconstruction performance of 20 randomly selected sample images.

[0179] Table 1 summarizes the computational complexity of the E-step and M-step of different image reconstruction methods. It can be seen that the LHP-GAMP-SBL method significantly reduces the computational complexity by using GAMP to estimate the mean and variance in the E-step and using Laplace hierarchical priors to update hyperparameters in the M-step.

[0180] Table 2 summarizes the two reconstruction performance indicators of ERROR and TIME on the three datasets. Each value in Table 2 represents the average of the results of 150 images. It can be seen from Table 2 that the LHP-GAMP-SBL method has better reconstruction accuracy than other methods. In terms of time, the LHP-GAMP-SBL method has a greater advantage than other GAMP-based methods. By using the Laplace hierarchical prior to estimate the hyperparameters in the M step, the image can be reconstructed with fewer iterations, which greatly improves the reconstruction efficiency.

[0181] Figure 4 The ERROR, TIME, PSNR and SSIM results of 20 randomly selected example images are shown in Figure 2. Figure 4 (a)-4(d) shows more intuitively that the LHP-GAMP-SBL method uses GAMP to estimate the mean and variance in the E step and uses LHP to update the hyperparameters in the M step, achieving better reconstruction efficiency. It also obtains higher PSNR and SSIM. The algorithm with a higher PSNR value has a smaller error and a higher SSIM. Figure 4 (e) and 4(f) show that, compared with the VGAMP-SBL and GGAMP-SBL algorithms, the LHP-GAMP-SBL method requires about 50 iterations to complete convergence, and the error is very small. Compared with other methods, the LHP-GAMP-SBL method of the present invention takes less running time per iteration, and the use of the Laplace hierarchical prior leads to better sparsity of the model.

[0182] The results of the above embodiments confirm the rationality and effectiveness of the Laplace hierarchical prior and damped GAMP algorithm used in steps S3-S7. The GAMP algorithm avoids matrix inversion and adopts a suitable damping coefficient to improve the convergence of the algorithm. The Laplace hierarchical prior is used to iteratively update the hyperparameters, which improves the sparsity and reconstruction efficiency of the model. The integration of the two further deepens the hierarchical structure of the model and enhances the representation ability of the model, and its efficiency and performance are significant.

[0183] like Figure 5 As shown, this embodiment discloses a sparse Bayesian remote sensing image reconstruction system based on Laplace hierarchical prior and GAMP, based on the above method embodiment, including the following modules:

[0184] Remote sensing image acquisition and preprocessing module: obtain remote sensing image data sets, set the pixels of remote sensing images, convert remote sensing images into amplitude maps and perform sparse processing on the amplitude maps;

[0185] Sparse Bayesian model building module: build sparse Bayesian model;

[0186] Computation module: calculate likelihood function, Laplace stratified prior distribution of sparse signals, marginal likelihood function;

[0187] Posterior distribution calculation module: Based on the likelihood function, Laplace stratified prior distribution and marginal likelihood function obtained by the calculation module, the posterior distribution of the sparse signal is calculated by Bayesian theorem;

[0188] Update module: obtain the posterior mean and posterior variance of the sparse signal and the update formula of the hyperparameters through the expectation maximization algorithm;

[0189] Mean and variance update module: introduce the GAMP algorithm to update the mean and variance;

[0190] Damping factor adding module: adding damping factor to GAMP structure;

[0191] Iteration module: determines whether the iteration termination condition is met. If so, it outputs the estimation of the sparse signal.

[0192] Inverse wavelet transform module: Perform inverse wavelet transform on the output sparse signal to obtain the reconstructed remote sensing image.

[0193] For other contents of this embodiment, please refer to the above method embodiment.

[0194] In summary, the present invention discloses an efficient sparse Bayesian remote sensing image reconstruction method and system based on Laplace hierarchical prior and GAMP. The preferred embodiment of the present invention selects a public remote sensing image data set, adjusts all original images to 256×256 pixels, converts multi-band remote sensing images into grayscale images and processes them into .mat format files to save the feature data of the image. The image is preprocessed by sparseness using wavelet transform. The present invention first introduces the GAMP algorithm to replace the E step of the classic EM algorithm on the basis of SBL, avoids the matrix inversion of solving the mean of the unknown signal, and greatly reduces the computational complexity of the model. The Laplace hierarchical prior is used to implement the M step in SBL, which has good robustness for the thick tail phenomenon in natural signals, and introduces the damping factor to enhance the convergence of the model. The integration of the two structures further deepens the hierarchical structure of the model, enhances the representation ability of the model, and improves the reconstruction efficiency and accuracy. The present invention improves the image reconstruction accuracy and convergence speed by using the GAMP algorithm and the Laplace hierarchical prior.

[0195] Compared with the remote sensing image reconstruction method in the prior art, the efficient sparse Bayesian remote sensing image reconstruction method and system based on Laplace hierarchical prior and GAMP of the present invention have faster reconstruction time and better reconstruction accuracy and reconstruction efficiency.

[0196] The embodiments described above are merely descriptions of preferred implementations of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope of the present invention.

Claims

1. Sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP, Its characteristics are The steps include: S1, obtain a remote sensing image data set, set the pixels of the remote sensing image, convert the remote sensing image into an amplitude map and perform sparse processing on the amplitude map; S2, establish a sparse Bayesian model; S3, calculate the likelihood function, the Laplace stratified prior distribution of the sparse signal, and the marginal likelihood function; S4, calculating the posterior distribution of the sparse signal by Bayes' theorem according to the likelihood function, Laplace stratified prior distribution and marginal likelihood function calculated in step S3; S5. Obtain the posterior mean and posterior variance of the sparse signal and the update formula of the hyperparameters through the expectation maximization algorithm; S6, introduce the GAMP algorithm to update the mean and variance; S7, adding damping factor to the GAMP structure; S8, judging whether the iteration termination condition is satisfied, if satisfied, outputting the estimation of the sparse signal; if not satisfied, executing step S5; S9. Perform inverse wavelet transform on the output sparse signal to obtain a reconstructed remote sensing image.

2. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 1, Its characteristics are: In step S2, a sparse Bayesian model is established, as follows: y=Dx+e in, is the observation vector, R represents the set of real numbers, is a sparse signal, is a known dictionary matrix, M represents the dimension of the observation vector y, N represents the dimension of the sparse signal x, is the noise vector; represents a set of real numbers, M represents the dimension of the observation vector y, and N represents the dimension of the sparse signal x. According to Bayes' theorem, the estimate of x is obtained by maximizing the posterior distribution p(x,ω,λ,θ|y): Among them, ω, λ, θ represent hyperparameters, and p(x, ω, λ, θ|y) is a conditional probability distribution.

3. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 2, Its characteristics are: In step S3, the likelihood function p(y|x,θ) is calculated according to the posterior distribution expression in step S2: Assume that the noise e has a mean of 0 and a variance of θ -1 Gaussian distribution of: in, Represents the Gaussian distribution symbol; The hyperparameter θ follows the Gamma prior: Among them, shape parameter a>0, scale parameter b>0; Γ(a θ ) is the Gamma function; the expectation of θ is expressed as: According to the posterior distribution expression in step S2, calculate the sparse Laplace stratified prior distribution p(x|λ): Using a hierarchical prior model, the first layer: Where ω=(ω 1 ,ω 2 ,...,ω N ),ω i is control x i The sparsity hyperparameter, i = 1, 2, ..., N; the second layer hyperparameter ω i Modeling: Among them, ω i ≥0, λ≥0; therefore, the Laplace stratified prior distribution of the sparse signal x is: Among them, the parameter λ obeys the Gamma hyper-prior: According to the posterior distribution expression in step S2, calculate the marginal likelihood function p(y): Where C = (θ -1 I+DΛ -1 D T ), I represents the identity matrix, is a diagonal matrix.

4. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 3, Its characteristics are: In step S4, according to Bayes' theorem, the posterior distribution of x p(x, ω, λ, θ|y) is calculated to also obey the Gaussian distribution:

5. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 4, Its characteristics are: In step S5, the mean μ of the sparse signal x is obtained by the expectation maximization algorithm. x and variance∑ x And the update formula of the hyperparameters (ω, λ, θ) is as follows: μ x =θ∑ x D T y S x =(θD T D+L) -1 The hyperparameter γ is updated as follows: in, is the digamma function, right The derivative of .

6. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 5, Its characteristics are: In step S6, the GAMP algorithm is introduced to replace the E step of the expectation-maximization algorithm to update the mean and variance: For a given prior distribution information p(x|ω) and marginal likelihood function p(y), GAMP constructs an iterative expression through quadratic approximation and Taylor series expansion; through the sum-product algorithm, the input function And the output function Defined as: Where m∈M represents the mth atom of the input function, n∈N represents the nth atom of the output function, and k m represents the mth atom of the noise-free model k = Dx, y m represents the mth atom of the observation vector y, Indicates k m The mean of Indicates k m The variance of x n represents the nth atom of the sparse signal x, ω n represents the nth atom of the hyperparameter ω, Represents x n The mean of Represents x n Variance of; intermediate variable and denote the approximate estimate of the sparse signal x and the noise-free model k=Dx, respectively. and Corresponding to the variance of x and k respectively; therefore, the update of the input function and the output function is: The iterative update of the mean and variance is expressed as: μ x ←G x , 7. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 6, Its characteristics are: In step S7, a damping factor is added to the GAMP structure to enhance the convergence of the model. The damping factor is calculated as follows: Among them, η represents the damping factor to be solved, const represents a constant, and Ω max represents the maximum singular value of dictionary D; Σ Ω represents the sum of the singular values ​​of D; α and β are constant terms; damping factor and is the root of the above formula, and Used to control Used to control x.

8. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 7, Its characteristics are: The iteration termination condition of step S8 is: Among them, μ x t represents the mean of the current iteration, μ x t-1 represents the mean of the previous iteration, ε em Indicates the iteration termination parameter; Determine whether the iterative convergence condition is met. If so, output the estimate of x: x = μ x t ; If not satisfied, execute step S5.

9. The sparse Bayesian remote sensing image reconstruction method based on Laplace hierarchical prior and GAMP as claimed in claim 8, Its characteristics are: In step S9, the output sparse signal x is subjected to inverse wavelet transform to obtain a reconstructed image, and the calculation formula of the reconstruction error is: Where F represents the F norm, represents the reconstructed remote sensing image, and S represents the real remote sensing image.

10. A sparse Bayesian remote sensing image reconstruction system based on Laplace hierarchical prior and GAMP, based on the method described in any one of claims 1 to 9, Its characteristics are Includes the following modules: Remote sensing image acquisition and preprocessing module: obtain remote sensing image data sets, set the pixels of remote sensing images, convert remote sensing images into amplitude maps and perform sparse processing on the amplitude maps; Sparse Bayesian model building module: build sparse Bayesian model; Computation module: calculate likelihood function, Laplace stratified prior distribution of sparse signals, marginal likelihood function; Posterior distribution calculation module: Based on the likelihood function, Laplace stratified prior distribution and marginal likelihood function obtained by the calculation module, the posterior distribution of the sparse signal is calculated by Bayesian theorem; Update module: obtain the posterior mean and posterior variance of the sparse signal and the update formula of the hyperparameters through the expectation maximization algorithm; Mean and variance update module: introduce the GAMP algorithm to update the mean and variance; Damping factor adding module: adding damping factor to GAMP structure; Iteration module: determines whether the iteration termination condition is met. If so, it outputs the estimation of the sparse signal. Inverse wavelet transform module: Perform inverse wavelet transform on the output sparse signal to obtain the reconstructed remote sensing image.

Images