CN107784278B - Using structured prior constraints to improve sparse image reconstruction accuracy and reduce complexity - Google Patents

Abstract

Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained, belong to Remote Sensing Image Compression technical field, it is characterized in that, it is a kind of method for introducing image information on the basis of more observation vector models to improve precision reduction complexity, successively the following steps are included: input observing matrix, the matrix of more observation vector compositions, the wavelet transformation number of plies, every layer of wavelet coefficient number；Initialize system model parameter, image wavelet coefficient, maximum number of iterations, precision property index；Image wavelet coefficient is reconstructed with the gibbs sampler based on Markov chain monte carlo method, then wavelet inverse transformation is carried out to matrix of wavelet coefficients after reconstruct, to obtain reconstructed image.The present invention can apply required precision property index to practical reconstructed image, complexity be reduced using more observation vector models, while adjusting reconstruction accuracy further through the method for the adjustment wavelet transformation number of plies, so that the quality of reconstructed image meets performance indicator requirement.

Description

Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained

Technical field

The invention belongs to field of signal processing, can be applied to the sparse reconstruct of image information, in particular to are seen based on more Survey the low complex degree structuring reconstructing method of vector model.

Background technique

According to compressed sensing (Compressive Sensing, CS) theory, if signal meets centainly in a certain transform domain Sparsity can then sample signal using the rate far below nyquist sampling rate, and realize to the accurate of signal Reconstruct.Currently, compressed sensing technology has application, such as Sub-nyquist sampling system, compression imaging in many engineering fields System, compression sensing network etc..Particularly, compressed sensing technology can also be applied on sparse image reconstruction, and Through being quite widely applied in key areas such as medical imagings.

However, existing compressed sensing algorithm is mostly based on single observation vector (SMV) when solving the problems, such as image reconstruction Entire sparse image is considered as one-dimensional vector and handled by model；The problems of the processing mode is: when our needs When high-definition picture is reconstructed, the one-dimensional vector dimension being drawn by large-scale image is quite high, therefore relates in algorithm And the observing matrix and covariance matrix arrived is in large scale, so that sizable storage pressure can be brought to existing computing platform With calculating pressure.For example, if opening the image of 1024 × 1024 dimensions using the compressed sensing algorithm process one based on SMV, Dimension will be more than 1.04 × 10 after the image array being related in algorithm transforms into a column vector⁶Dimension, stores such scale Covariance matrix need be more than 64GB computing platform memory；Even if computing platform memory is sufficiently large, matrix storage can satisfy Demand simultaneously makes these methods can be with successful operation, but the computing resource for the calculating time and occupancy thus expended is also to be difficult to hold It receives.Therefore, in high-definition picture reconstruction, do not have just with the compressed sensing algorithm based on SMV model existing Sincere justice.

For this purpose, the present invention fully considers the characteristics of image sparse reconstruct, problem is characterized again using MMV model. MMV model handles each column vector of original signal as independent subsignal, and in array signal processing, cognition The fields such as radio have been achieved for corresponding research achievement, but are not applied to the sparse reconstruction of image.Therefore, such as More observation vector models are introduced large-scale image reconstruction to fruit by us, and each column vector of image array are regarded as It is an observation vector in MMV model, then such model characteristic manner does not further relate to convert image to high-dimensional single The operation of vector, corresponding problem scale will be reduced significantly, so that being solved using the principle of compressed sensing big The reconstruct of scale image becomes the engineering problem that can be realized in most of computing platforms.But if by existing solution The MMV restructing algorithm of other field problem directly applies to when solving the problems, such as image reconstruction, then having ignored image information is had Architectural characteristic, to cause the serious loss and the deterioration of reconstruction accuracy of information.

It therefore, is the reconstruction accuracy deterioration problem being likely to occur when further overcoming directly using MMV model, institute of the present invention Reconstructing method is stated on the basis redescribed using MMV model to problem, introduces in image information and potentially ties Structureization is prior-constrained, be embodied as of both structuring prior information: on the one hand, image it is dilute on wavelet transformed domain Sparse coefficient has the probability dependency of quaternary tree shape structure, on the other hand, between the column vector of image array and dependent, and It is that there is certain spatial coherence.

In conclusion the present invention is by introducing high-definition picture reconstruction for MMV model, can be realized makes its storage Complexity and computation complexity decline to a great extent；Further, prior-constrained by introducing the structuring of image information, it ensure that The precision of image reconstruction under MMV model.

Summary of the invention

(1) technical problems to be solved

The storage complexity faced present invention aim to address the sparse reconstruct of high-definition picture and time are complicated Spend the limited problem of higher and reconstruction accuracy.

(2) technical solution

In order to solve the above-mentioned technical problem, the present invention proposes a kind of structure based on more observation vectors characterization of low complex degree Change reconstructing method.The purpose of the present invention is what is be achieved through the following technical solutions: input observing matrix, more observation vector compositions Matrix, the wavelet transformation number of plies, every layer of wavelet coefficient number；Initialize system model parameter, image wavelet coefficient, greatest iteration time Number, precision property index；Weight is carried out to image wavelet coefficient with the gibbs sampler based on Markov chain monte carlo method Structure, then wavelet inverse transformation is carried out to matrix of wavelet coefficients after reconstruct, so that reconstructed image is obtained,

1, sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained, which is characterized in that be one Kind is the remote sensing images X of N × N=1024 × 1024 to resolution ratio₀After being redescribed using more observation vector models, then benefit Use X₀Sparse coefficient on wavelet transformed domain has the probability dependency and more observation vector matrixes of quaternary tree shape structure The structuring of spatial coherence these two aspects between column vector it is prior-constrained come reconstruct resolution ratio be N × N=1024 × 1024 remote sensing images X^*Method, successively realize according to the following steps in a computer:

Step (1) inputs following parameter to computer:

The remote sensing images matrix X of N × N=1024 × 1024₀, abbreviation matrix X₀,

It is generated with the tool box MATLAB and is based on the matrix X₀Quaternary tree shape structure 2-d discrete wavelet transformation of coefficient Matrix Θ, resolution ratio are N × N=1024 × 1024, and the number of plies S of two-dimensional discrete wavelet conversion is based on the matrix X₀Point Resolution is N × N=1024 × 1024, the maximum value S of the number of plies_max=8, level number s=1,2 ..., 8, quaternary tree shape structure The number of each layer of wavelet coefficient is (NC) in wavelet coefficient_s, root node is in first layer S₁, appoint and take S=4,

With the tool box MATLAB generate M × N random observation matrix, M≤N, observing matrix abbreviation matrix Φ, similarly hereinafter,

The matrix Y of more observation vectors composition based on the matrix Θ M × N generated, this is to be with the matrix Φ For the matrix X under conditions of observing matrix₀Observation, i.e. Y=Φ Θ,

The remote sensing images X obtained after input reconstruct^*Two image property indexs: the setting value of Y-PSNR PSNR, The setting value of mean square error MSE,

Step (2) initializes Modulus Model parameter, comprising:

1. the regularized covariance Matrix C of the wavelet coefficient of different level numbers₁, C₂..., C₄, when initialization is all unit square Battle array, abbreviation C_s={ C₁, C₂..., C₄,

2. controlling each C_sThe parameter lambda of value, enables λ be initialized as 2,

3. the precision parameter for the Gaussian Profile that observation noise is obeyed is the α reciprocal of the variance of observation noise_n, at the beginning of enabling it Beginning turns to α_n=10,

4. the precision parameter α for the Gaussian Profile that the wavelet coefficient of different level numbers is obeyed₁, α₂..., α₄, it is each layer small echo The inverse of variance between coefficient column vector, when initialization, are equal to 1,

5. setting: control the probability that each layer wavelet coefficient takes zero:

For first layer S₁: control matrix of wavelet coefficients Θ⁽¹⁾In elementValue take zero probability to be initialized asAnd enable W₁=0.3, correspondingly,The probability of negated zero is 1-W₁, wherein symbol Θ⁽¹⁾、In on be designated as the level number of each layer matrix of wavelet coefficients Θ, subscript (ij) is sequence of each element W in row, column Number, similarly hereinafter,

For the second layer to the 4th layer of S₂, S₃, S₄:

If matrix of wavelet coefficients Θ^(s), the element of 2≤s≤4The coefficient value of its father's node under quad-tree structureThen controlValue takes zero probability2 s≤4 <, and initialize W0_s=0.3, symbol W0_sIn Symbol " 0 " indicate father's nodeCoefficient value be zero, symbol " π " indicate father's node, similarly hereinafter,

If matrix of wavelet coefficients Θ^(s), the element of 2≤s≤4The coefficient value of its father's node under quad-tree structureThen controlValue takes zero probability2 s≤4 <, and initialize W1_s=0.7, symbol W1_sIn Symbol " 1 " indicate father's nodeCoefficient value be 1,

6. system hyper parameter a₀,b₀,c₀,d₀,e₀,f₀, the hyper parameter is that the parameter in a kind of pair of step (2) is obeyed The parameter obtained after prior probability distribution parametrization, in which:

a₀,b₀,c₀,d₀In initialization all equal to 10^-6, it is constant for controlling the parameter of probability-distribution function,

Symbol e₀It include: (ec)₁, and (ecz)_s, 2≤s≤4:

(ec)₁For controlling W₁The probability density function obeyed, value are as follows:

(ec)₁=0.9 × (NC)₁,

(ecz)₂, (ecz)₃, (ecz)₄, it is respectively used to W0 when 2≤s of control≤4₂, W0₃, W0₄The probability density obeyed point Cloth function, value are as follows:

(ecz)₂=0.1 × (NC)₂, (ecz)₃=0.1 × (NC)₃, (ecz)₄=0.1 × (NC)₄, (eco)₂, (eco)₃, (eco)₄It is respectively used to W1 when 2≤s of control≤4₂, W1₃, W1₄The probability density function obeyed, value are as follows:

(eco)₂=0.5 × (NC)₂, (eco)₃=0.5 × (NC)₃, (eco)₄=0.5 × (NC)₄,

Symbol f₀It include: (fc)₁, and (fcz)_s, 2≤s≤4:

(fc)₁For controlling W₁The probability density function obeyed, value are as follows:

(fc)₁=0.1 × (NC)₁,

(fcz)₂, (fcz)₃, (fcz)₄It is respectively used to W0 when 2≤s of control≤4₂, W0₃, W0₄The probability density obeyed point Cloth function, value are as follows:

(fcz)₂=0.9 × (NC)₂, (fcz)₃=0.9 × (NC)₃, (fcz)₄=0.9 × (NC)₄, (fco)₂, (fco)₃, (fco)₄It is respectively used to W1 when 2≤s of control≤4₂, W1₃, W1₄The probability density function obeyed, value are as follows:

(fco)₂=0.5 × (NC)₂, (fco)₃=0.5 × (NC)₃, (fco)₄=0.5 × (NC)₄,

7. the parameter t, T of initialization control iterative process_max, in which:

T is the serial number of the number of iterations, T_maxFor maximum number of iterations, thus the setting value of the serial number of the number of iterations be t=1, 2,...,T_max}；

Step (3) successively calculate by following procedure the iterative process of matrix of wavelet coefficients Θ:

Step (3.1) successively carries out first time iteration according to the following steps, enables t=1,

Step (3.1.1) sets various system model parameters in constraint condition and step (2):

Constraint condition are as follows: the remote sensing images matrix X of N × N=1024 × 1024₀, the wavelet systems of N × N=1024 × 1024 Matrix number Θ is quad-tree structure, S_max=8, appoint and takes S=4,

When initial, the matrix of wavelet coefficients Θ of each layer is full null matrix,

The initial value according to system model parameter and wavelet coefficient is successively calculated as follows in step (3.1.2), calculates intermediate Variable

Wherein:

||·||₂Indicate the 2- norm of amount of orientation or matrix,

The element Θ arranged for the i-th row of Θ matrix, jth_ijThe precision parameter for the probability distribution obeyed,

The element Θ arranged for the i-th row of Θ matrix, jth_ijThe mean of a probability distribution parameter obeyed,

The element Θ arranged for the i-th row of Θ matrix, jth_ijThe adjustment parameter of zero probability is taken,

α_sAnd α_nSystem model parameter value when expression is initial,System model when expression is initial The intermediate variable that parameter is calculated,

Φ_-iIndicate to remove the matrix of other all column compositions except its i-th column in Φ,

Φ_·iIndicate the vector of all elements composition of the i-th column of Φ matrix,Indicate vector Φ_·iTransposition,

Θ_(-i)jIt indicates in vector theta_·jThe vector of the middle other all elements composition removed except its i-th of element,

Step (3.1.3), according to intermediate variable obtained in step (3.1.2)Calculate wavelet coefficient Θ institute Posterior probability distribution p (the Θ of obedience_ij):

Wherein: δ₀The uni-impulse function at origin is represented,

Posterior probability distribution p (the Θ obeyed according to wavelet coefficient Θ_ij), sampling obtains the i-th row of matrix Θ, jth column Element Θ_ijValue to get arrive matrix Θ sampled value,

Step (3.1.4), according to the sampling of matrix Θ obtained in the initial value of system model parameter and step (3.1.3) Value calculates system model parameter is obeyed in first time iteration posterior probability Density Distribution and Matrix C_sFirst time change For result:

Wherein:

The value set of footmark s is { 1 ..., 4 }, it is noted that only in the S of input >=2, and the number of plies s currently considered When greater than 1, just there is W0_sAnd W1_s,

The Frobenius norm of representing matrix,

||·||₀The 0- norm of representing matrix (or vector),

Θ^(s)It is that s layers of wavelet conversion coefficient combines the matrix to be formed according to natural order,

R_sIt is Θ^sLine number,Indicate Θ^sR row,

Θ_s,kIndicate s layers of k-th of wavelet coefficient, Θ_π(s,k)It is Θ_s,kFather's node coefficient value,

P (x)=Gamma (a, b) indicates that stochastic variable x obeys the gamma that parameter is (a, b) and is distributed, and expanded form isWherein Γ () indicates gamma function, by taking a is independent variable as an example, the expanded form of gamma function For

P (x)=Beta (a, b) indicates that stochastic variable x obeys the beta that parameter is (a, b) and is distributed, and expanded form is

Indicating indicative function, i.e., functional value takes 1 when the expression formula in bracket is true, otherwise functional value takes 0,

According to above system model parameter α_n、α_s、W₁、W0_sAnd W1_sThe posterior probability Density Distribution obeyed, to system mould Shape parameter is sampled, and updates corresponding system model parameter value, for iterative process next time,

Step (3.1.5), the sampled value of matrix Θ obtained in storing step (3.1.3), first time iterative process are formal Terminate,

Step (3.2) carries out second to T according to the description of step (3.1.1)~step (3.1.5)_maxSecondary iteration Calculating,

Step (3.3), to the serial number of all the number of iterations, t={ 1 ..., T_max, it calculates storage in each iteration and obtains All matrix Θ arithmetic mean of instantaneous value, obtained average value is denoted as Θ^*, to average value Θ^*Carry out S layers of 2-d wavelet inversion It changes, obtains X^*,

Step (3.4), counts reconstructed image X^*Relative to original image X₀Error extension, PSNR and MSE, if this Two indexes meet the performance indicator requirement of input, then algorithm stops, the wavelet coefficient number of plies otherwise inputted in set-up procedure (1) S, and S:=S+1 is enabled, go to step (2) again, again operation initialization and iterative process, until reconstructed image is opposite Meet performance indicator requirement in the error extension of original image.

Beneficial effect

The present invention gives a kind of sparse reconstructing methods of the high-definition picture of low complex degree.The reconstructing method is based on MMV model is designed；To further increase reconstruction accuracy, the structuring that this method introduces image is prior-constrained.This method Have the advantage that (1) makes the storage complexity of problem by O (N by introducing MMV model⁴) it is reduced to O (N²), it calculates complicated Degree is by O (M²N²) it is reduced to O (MN), to provide feasible scheme for high-definition picture reconstruction；(2) pass through introducing The constraint of structuring prior probability distribution, improves the accuracy of image reconstruction under MMV model；(3) it utilizes and is derived in this method The Gibbs sampling method arrived, can conclude that obtain image coefficient and model parameter under the structuring MMV model it is theoretical most Excellent solution.

Detailed description of the invention

Fig. 1 is the flow chart that the present invention realizes sparse image reconstruction.

Fig. 2 features the quaternary tree shape probability dependency structure of image multilayer wavelet transform, and wherein s is that multi-level Wavelet Transform becomes The number of plies changed.

Fig. 3 illustrates the stratification Bayesian model based on structuring prior information that the present invention establishes.Wherein a₀,b₀, c₀,d₀,e₀,f₀The hyper parameter of picture engraving coefficient and model parameter prior distribution, Θ, Y be respectively image wavelet coefficient and Based on the sparse observing matrix that more observation vector models obtain, α_s,C_s,ω_sIt is the parameter for portraying the probability density characteristics of Θ, α_n It is the parameter for portraying the probability density characteristics of Y.

Fig. 4 is that image reconstructing method proposed by the present invention (cM-SSBL) and the existing sparse reconstruct based on MMV model are calculated Contrast on effect of the method (OMP-MMV, M-FOCUSS and MSBL) in the sparse reconstruction of processing large-scale image.Horizontal axis is normalizing Change sample rate M/N, the longitudinal axis is Relative reconstruction error；It indicates to use OMP-MMV algorithm,Using M-FOCUSS Algorithm,It indicates to use MSBL algorithm,It indicates to use cM-SSBL algorithm.

Fig. 5 is that image reconstructing method proposed by the present invention (cM-SSBL) and the existing sparse reconstruct based on MMV model are calculated Contrast on effect of the method (OMP-MMV, M-FOCUSS and MSBL) in the sparse reconstruction of processing large-scale image.Horizontal axis is normalizing Change sample rate M/N, the longitudinal axis is Y-PSNR (PSNR)；It indicates to use OMP-MMV algorithm,Using M- FOCUSS algorithm,It indicates to use MSBL algorithm,It indicates to use cM-SSBL algorithm.

Fig. 6 is to use resolution ratio for 1024 × 1024 standard testing image " cameraman ", is in normalization sample rate In the case where 0.4, image is carried out to OMP-MMV, M-FOCUSS, MSBL and tetra- kinds of methods of cM-SSBL proposed by the invention The visual effect display diagram of reconstitution experiments.Fig. 6 (a) is the reconstruction result of OMP-MMV algorithm and corresponding PSNR, Fig. 6 (b) are M- The reconstruction result of FOCUSS algorithm and corresponding PSNR, Fig. 6 (c) are the reconstruction result and corresponding PSNR, Fig. 6 of MSBL algorithm (d) be cM-SSBL algorithm proposed by the present invention reconstruction result and corresponding PSNR.

Fig. 7 is the reference curve during parameter preferably relevant to iteration control, and horizontal axis is Y-PSNR, the longitudinal axis The number of iterations, in order to show the influence of the variation of the number of iterations to reconstruct error performance, Fig. 7 intercepted preferred process carry out to Finally, Y-PSNR situation of change when the number of iterations is relatively fewer.

Specific embodiment

With reference to the accompanying drawings and examples, specific embodiments of the present invention will be described in further detail.Implement below Example is not limited the scope of the invention for illustrating the present invention.

As shown in Figure 1, the purpose of the present invention is be achieved through the following technical solutions: input observing matrix, observe more to Measure the matrix of composition, the wavelet transformation number of plies, every layer of wavelet coefficient number；Initialization system model parameter, image wavelet coefficient, most Big the number of iterations, precision property index；With the gibbs sampler based on Markov chain monte carlo method to image wavelet system Number is reconstructed, then carries out wavelet inverse transformation to matrix of wavelet coefficients after reconstruct, to obtain reconstructed image.The present invention can be to reality Border reconstructed image applies required precision property index, reduces complexity using more observation vector models, while further through adaptive The method adjustment the number of iterations answered, so that the precision of reconstructed image meets performance indicator requirement.

Now for carrying out sparse reconstruct to standard testing image X of the width dimension for N × N (N=1024),

Sparse observing matrix is determined first, and the observation of M × N (M=512) is generated using MATLAB generating random number tool box Matrix Φ,

According to the basic principle that image sparse characterizes, for original remote sensing images under conditions of matrix Φ is characterization substrate X₀It is observed, obtains the matrix Y of M × N of observation vector composition,

The number of plies S=4 of two-dimensional discrete wavelet conversion is set, then the number of every layer of wavelet coefficient also determines therewith: in N= In the case where 1024, S=4, the 0th layer of wavelet coefficient number is (NC)₀=3 × 64²=12288, the 1st layer of wavelet coefficient number be (NC)₁=3 × 128²=49152, the 2nd layer of wavelet coefficient number is (NC)₂=3 × 256²=196608, the 3rd layer of wavelet coefficient Number is (NC)₃=3 × 512²=786432；

(this sentences MSE=16.3 setting performance indicator to setting reconstructed image performance indicator PSNR=36, and effect is of equal value )；

Next, the sparse image reconstruction algorithm based on structuring prior information proposed according to the present invention, successively carries out Following solution procedure:

The wavelet coefficient Θ for initializing N × N enables Θ be equal to full null matrix；

Initialize system model parameter, comprising:

The regularized covariance Matrix C of the wavelet coefficient of different level numbers₀,C₁,...,C_S-1, unit matrix is all when initial, Abbreviation C_S,C_S={ C₀,C₁,...,C_S-1}；

Control C_SThe parameter lambda of size, enables λ=2；

The precision parameter for the Gaussian Profile that observation noise is obeyed is the inverse of the variance of observation noise, is denoted as α_n, α_n= 10；

The precision parameter α for the Gaussian Profile that the wavelet coefficient of different level numbers is obeyed₀,α₁,...,α_S-1, it is s layers respectively The inverse of covariance between wavelet coefficient column vector, when initialization, are equal to 1；

W₁Control first layer wavelet coefficient takes zero probability, and sequence s, s=2 ..., S, W0 are designated as under₂..., W0_SAnd W1₂..., W1_SThe probability that wavelet coefficient takes zero is controlled when being all s > 1,

For scale be N × N wavelet coefficient Θ the i-th row corresponding at different sequence s, jth column member The setting value of element, i=1,2, i .N, j=1,2, j., N:

W₁Θ when for first layer s=1_ij|_S=1WhenValue,

W0_sFor s > 1 and Θ_π(ij)When=0Value,Wherein, Θ_π(ij)It is the Θ in the case where Θ is quad-tree structure_ijFather's node coefficient value, similarly hereinafter, s value range be { 2 ..., S }, i= 1,2, i .N, j=1,2, j., N,

W1_sFor s > 1 and Θ π_(ij)When=1Value,S value range For { 2 ..., S }, i=1,2, i .N, j=1,2, j., N,

Initialization system hyper parameter a₀,b₀,c₀,d₀,e₀,f₀, in which:

a₀,b₀,c₀,d₀In initialization all equal to 10-⁶,

Symbol e₀It include: (ec)₁, and (ecz)₂..., (ecz)_s(eco)₂..., (eco)_s:

(ec)₁For controlling W₁The probability density function obeyed, value are as follows:

(ec)₁=0.9 × (NC)₁,

(ecz)₂..., (ecz)_SFor controlling W0₂..., W0_SThe probability density function obeyed, value are as follows:

(ecz)₂=0.1 × (NC)₂..., (ecz)_S=0.1 × (NC)_S, (eco)₂..., (eco)_SFor controlling W1₂..., W1_SThe probability density function obeyed, value are as follows:

(eco)₂=0.5 × (NC)₂..., (eco)_S=0.5 × (NC)_S,

Symbol f₀It include: (fc)₁, and (fcz)₂..., (fcz)_S(fco)₂..., (fco)_S:

(fc)₁For controlling W₁The probability density function obeyed, value are as follows:

(fc)₁=0.1 × (NC)₁,

(fcz)₂..., (fcz)_SFor controlling W0₂,...,W0_S-1The probability density function obeyed, value is such as Under:

(fcz)₂=0.9 × (NC)₂..., (fcz)_S=0.9 × (NC)_S,

(fco)₂..., (fco)_SFor controlling W1₂,...,W1_S-1The probability density function obeyed, value is such as Under:

(fco)₂=0.5 × (NC)₂..., (fco)_s=0.5 × (NC)_s；

The parameter t, T of initialization control iterative process_max, in which:

T is the serial number of the number of iterations,

Set maximum number of iterations T_max=55.

Successively calculate by following procedure the iterative process of wavelet coefficient Θ,

Step (1) determines the initial value of system model parameter according to current wavelet coefficient number of plies S；

Step (2), T_maxFor the good constant of current setting, the periodicity at most iterated to calculate is represented, On be designated as the number of iterations,

Step (3) successively realizes iterative calculation for the first time, serial number t=t as follows⁽¹⁾=1,

Initial value of the step (3.1) using the initial value of each system model parameter and wavelet coefficient as first time iteration,

The initial value according to system model parameter and wavelet coefficient is successively calculated as follows in step (3.2), calculates intermediate become Amount

Step (3.3) calculates the Posterior probability distribution p (Θ that wavelet coefficient Θ is obeyed_ij):

Step (3.4), the Posterior probability distribution p (Θ obeyed according to wavelet coefficient Θ_ij), sampling obtains the of matrix Θ I row, jth column element Θ_ijValue to get the sampled value for arriving matrix Θ, according to the initial value of system model parameter and matrix Θ Sampled value calculates system model parameter is obeyed in first time iteration posterior probability Density Distribution and Matrix C_sFirst Secondary iteration result,

Step (4) carries out second to T according to the description of step (3.1)~step (3.4)_maxThe calculating of secondary iteration,

Step (5), to all Θ^(t),Calculate all Θ^(t)Arithmetic mean of instantaneous value, what is obtained is flat Mean value is denoted as Θ^*, to average value Θ^*The wavelet inverse transformation for carrying out S layers, obtains X^*,

Step (6), counts reconstructed image X^*PSNR, do not meet the performance indicator of PSNR=36 yet at this time, then by T_max Increase 10, at this time T_max=65, and wavelet coefficient number of plies S is increased by 1, go to step (1) again, again operation iterative calculation Process, PSNR > 36 of this obtained reconstructed image then stop operation, export reconstructed image.

On the basis of above-mentioned algorithm executes step, the prior-constrained MMV of structuring is introduced to of the present invention here Restructing algorithm is illustrated compared to the reconstruction property gain that existing MMV restructing algorithm can be provided.We obtain reconstruct The Relative reconstruction error and Y-PSNR of wavelet coefficient investigate existing 3 kinds of MMV restructing algorithms as evaluation criterion respectively The standard that OMP-MMV, M-FOCUSS, MSBL and cM-SSBL algorithm proposed by the present invention are 1024 × 1024 to resolution ratio is surveyed Attempt the effect being reconstructed as " cameraman ".

Fig. 4 illustrates Relative reconstruction error with the situation of change of normalization sample rate M/N, it can be seen from the figure that returning One changes sample rate in the variation range from 0.25 to 0.55, cM-SSBL algorithm and OMP-MMV, M-FOCUSS, and MSBL algorithm is compared It is obviously improved in terms of error performance.By taking experimental result when normalizing sample rate and being 0.4 as an example, cM-SSBL algorithm Relative reconstruction error is 0.03 or so, and the Relative reconstruction error of OMP-MMV, M-FOCUSS, MSBL algorithm is respectively 0.075, 0.06,0.05 or so, it is 0.4 in normalization sample rate in other words, test image is that resolution ratio is 1024 × 1024 In the case where " cameraman ", cM-SSBL algorithm is compared with OMP-MMV algorithm, precision improvement 60% or so, with M-FOCUSS Algorithm is compared, precision improvement 50% or so, and compared with MSBL algorithm, precision improvement 40% or so.

In the case that Fig. 5 is illustrated using Y-PSNR as evaluation criterion, OMP-MMV, M-FOCUSS, MSBL, and The effect that " cameraman " test image that resolution ratio is 1024 × 1024 is reconstructed in cM-SSBL algorithm proposed by the present invention Fruit, it can be seen from the figure that in normalization sample rate in the variation range from 0.25 to 0.55, cM-SSBL algorithm and OMP- MMV, M-FOCUSS, MSBL algorithm are equally obviously improved compared to the aspect of performance in Y-PSNR.

In order to compare OMP-MMV, M-FOCUSS, MSBL and cM-SSBL algorithm proposed by the present invention from visual effect Performance, we use resolution ratio to carry out for 1024 × 1024 standard testing image " cameraman " to above-mentioned three kinds of methods Image reconstruction test, normalization sample rate are 0.4.From fig. 6, it can be seen that the quality reconstruction of cM-SSBL algorithm will be substantially better than OMP-MMV, M-FOCUSS and MSBL algorithm.Meanwhile evaluated from objective PSNR angle, cM-SSBL method compared to OMP-MMV, M-FOCUSS, MSBL algorithm also provide 5 gains for arriving 8dB.This explanation, introducing structuring proposed by the present invention are first It is effective for testing the idea of constraint.

In order to examine the robustness of cM-SSBL algorithm, we have chosen different test images to OMP-MMV, M- FOCUSS, MSBL and cM-SSBL algorithm are tested, and are counted Relative reconstruction error and Y-PSNR respectively and are carried out pair According to experimental result is as shown in table 1.Table 1 be using different test images (resolution ratio of all test images is 1024 × 1024, normalization sample rate image reconstructing method (cM-SSBL) proposed by the present invention and existing is based on MMV mould when being 0.4) Property of the sparse restructing algorithm (OMP-MMV, M-FOCUSS and MSBL) of type when handling the sparse reconstruction of large-scale image It can comparison.

Algorithms of different reconstruction property comparison of the table 1 for 1024 × 1024 dimension standard testing images

Table 1 illustrates the test result of 4 kinds of standard testing images, and 4 kinds of images are " Lena " respectively, " Boats ", " Barbara ", and " Goldhill ", all test image resolution ratio are 1024 × 1024, and normalization sample rate is 0.4.From Table 1 can visually see, cM-SSBL algorithm under different test images with OMP-MMV, M-FOCUSS and MSBL algorithm It compares, can be obviously improved reconstruction property.

From figure 7 it can be seen that the number of iterations is more, the Y-PSNR of reconstruct is higher, further considers timeliness Can, we preferably go out T in embodiments_max=65 as the iteration total degree in emulation example.

A specific embodiment of the invention is described in conjunction with attached drawing above, but these explanations cannot be understood to limit The scope of the present invention, protection scope of the present invention are limited by appended claims, any in the claims in the present invention base Change on plinth is all protection scope of the present invention.

Claims (1)

1. The method of improving sparse image reconstruction accuracy and reducing complexity with structured prior constraints is characterized in that it is a method to re-describe a remote sensing image X ₀ with a resolution of N×N=1024×1024 using a multi-observation vector model Then, using the structural prior constraints that the sparse coefficients of X _o in the wavelet transform domain have the probability dependence of the quad-tree structure and the spatial dependence of the multi-observation vector matrix between the column vectors to reconstruct The method of obtaining a remote sensing image X ^* with a resolution of N×N=1024×1024 is implemented in the computer according to the following steps: Step (1), input the following parameters to the computer: N×N=1024×1024 remote sensing image matrix X ₀ , referred to as matrix X ₀ , The MATLAB toolbox is used to generate a two-dimensional discrete wavelet coefficient transform matrix Θ based on the quad-tree structure of the matrix X _o , and its resolution is N×N=1024×1024, the number of layers S of the two-dimensional discrete wavelet transform, based on The resolution of the matrix X ₀ is N×N=1024×1024, the maximum value of the number of layers S _max = 8, the layer number s = 1.2, . . . , 8, among the wavelet coefficients of the quadtree structure The number of wavelet coefficients in each layer is (NC) _s , the root node is in the first layer S ₁ , and S=4 is chosen arbitrarily, The M×N random observation matrix generated by the MATLAB toolbox, M≤N, the observation matrix is referred to as matrix Φ, the same below, The matrix Y composed of M×N multi-observation vectors generated based on the matrix Θ is the observed value of the matrix X ₀ under the condition that the matrix Φ is the observation matrix, that is, Y=ΦΘ, Enter the two image performance indicators of the reconstructed remote sensing image X ^* : the set value of the peak signal-to-noise ratio PSNR, the set value of the mean square error MSE, Step (2), initialize the coefficient model parameters, including: _{①Regularized covariance matrices} C ₁ , _{C 2} , _. }, ②Control the parameter λ of each C _s value, let λ be initialized to 2, ③ The accuracy parameter of the Gaussian distribution obeyed by the observation noise is the reciprocal α _n of the variance of the observation noise, which is initialized as α _n =10, ④ The precision parameters α ₁ , α ₂ , ..., α ₄ of the Gaussian distribution obeyed by the wavelet coefficients of different layer numbers are the reciprocal of the variance between the column vectors of the wavelet coefficients of each layer, which are all equal to 1 during initialization, ⑤Setting: Control the probability that the wavelet coefficients of each layer take zero: For the first layer S ₁ : control the elements in the wavelet coefficient matrix Θ ⁽¹⁾ The probability that the value of is zero is initialized as And let W ₁ =0.3, correspondingly, The probability of taking a non-zero value is 1-W ₁ , where the symbols Θ ⁽¹⁾ , The superscript in is the layer number of the wavelet coefficient matrix Θ of each layer, and the subscript (ij) is the serial number of each element W in the row and column, the same below, For the second to fourth layers S ₂ , S ₃ , S ₄ : If the wavelet coefficient matrix Θ ^(s) , the elements of 2≤s≤4 The coefficient value of its parent node in the quadtree structure then control the probability that the value will be zero 2<s≤4, and initialize W0 _s = 0.3, the symbol "0" in the symbol W0 _s represents the parent node The coefficient value of is zero, the symbol "π" represents the parent node, the same below, If the wavelet coefficient matrix Θ ^(s) , the elements of 2≤s≤4 The coefficient value of its parent node in the quadtree structure then control the probability that the value will be zero 2<s≤4, and initialize W1 _s = 0.7, the symbol "1" in the symbol W1 _s represents the parent node The coefficient value of 1 is 1, ⑥System hyperparameters a ₀ , b ₀ , c ₀ , d ₀ , e ₀ , f ₀ , the hyperparameters are parameters obtained by parameterizing the prior probability distribution obeyed by the parameters in step (2). ,in: a ₀ , b ₀ , c ₀ , d ₀ are all equal to 10 ^-6 when initialized, and are used to control the parameters of the probability distribution function, which are constants, The symbol e ₀ includes: (ec) ₁ , and (ecz) _s , 2≤s≤4: (ec) ₁ is used to control the probability density distribution function obeyed by W ₁ , and its values are as follows: (ec) ₁ =0.9×(NC) ₁ , (ecz) ₂ , (ecz) ₃ , (ecz) ₄ are respectively used to control the probability density distribution function that W0 ₂ , W0 ₃ , and W0 ₄ obey when 2≤s≤4, and the values are as follows: (ecz) ₂ =0.1×(NC) ₂ , (ecz) ₃ =0.1×(NC) ₃ , (ecz) ₄ =0.1×(NC) ₄ , (eco) ₂ , (eco) ₃ , (eco) ₄ are respectively used to control the probability density distribution function that W1 ₂ , W1 ₃ , and W1 ₄ obey when 2≤s≤4, and the values are as follows: (eco) ₂ =0.5×(NC) ₂ , (eco) ₃ =0.5×(NC) ₃ , (eco) ₄ =0.5×(NC) ₄ , The symbol f ₀ includes: (fc) ₁ , and (fcz) _s , 2≤s≤4: (fc) ₁ is used to control the probability density distribution function obeyed by W ₁ , and its values are as follows: (fc) ₁ =0.1×(NC) ₁ , (fcz) ₂ , (fcz) ₃ , (fcz) ₄ are respectively used to control the probability density distribution functions that W0 ₂ , W0 ₃ , and W0 ₄ obey when 2≤s≤4, and the values are as follows: (fcz) ₂ =0.9×(NC) ₂ , (fcz) ₃ =0.9×(NC) ₃ , (fcz) ₄ =0.9×(NC) ₄ , (fco) ₂ , (fco) ₃ , (fco) ₄ are respectively used to control the probability density distribution functions that W1 ₂ , W1 ₃ , and W1 ₄ obey when 2≤s≤4. The values are as follows: (fco) ₂ =0.5×(NC) ₂ , (fco) ₃ =0.5×(NC) ₃ , (fco) ₄ =0.5×(NC) ₄ , ⑦Initialize the parameters t, _Tmax that control the iterative process, where: t is the number of iterations, and _Tmax is the maximum number of iterations, so the set value of the number of iterations is t={1, 2, ..., T _max }; Step (3), carry out the iterative process of calculating wavelet coefficient matrix Θ successively according to the following process: Step (3.1), perform the first iteration according to the following steps in turn, let t=1, Step (3.1.1), set constraints and various system model parameters in step (2): Constraints are: N×N=1024×1024 remote sensing image matrix X ₀ , N×N=1024×1024 wavelet coefficient matrix Θ is a quadtree structure, S _max =8, arbitrarily take S=4, Initially, the wavelet coefficient matrix Θ of each layer is an all-zero matrix, Step (3.1.2), calculate the intermediate variables according to the initial value of the system model parameters and wavelet coefficients according to the following formula in: ||·|| ₂ represents the 2-norm of the vector or matrix, is the precision parameter of the probability distribution that the element Θij of the i-th row and j-th column of the Θ _matrix obeys, is the mean parameter of the probability distribution to which the element Θ _ij of the i-th row and j-th column of the Θ matrix obeys, is the adjustment parameter of the probability that the element Θ _ij of the i-th row and the j-th column of the Θ matrix takes zero, α _s and α _n represent the initial system model parameter values, represents an intermediate variable calculated from the initial system model parameters, Φ _-i represents a matrix composed of all columns except the i-th column in Φ, Φ _.i represents a vector composed of all elements of the i-th column of the Φ matrix, represents the transpose of the vector _Φ.i , Θ _{(-i)j represents} a vector composed of all elements except the i-th element in the vector Θ.j, Step (3.1.3), according to the intermediate variable obtained in step (3.1.2) Compute the posterior probability distribution p(Θ _ij ) obeyed by the wavelet coefficients Θ: where: δ ₀ represents the unit shock function at the origin, According to the posterior probability distribution p(Θ _ij ) obeyed by the wavelet coefficients Θ, the value of the element Θ _ij in the i-th row and the j-th column of the matrix Θ is obtained by sampling, that is, the sampling value of the matrix Θ is obtained, Step (3.1.4), according to the initial value of the system model parameters and the sampling value of the matrix Θ obtained in step (3.1.3), calculate the posterior probability density distribution obeyed by the system model parameters in the first iteration, and the matrix The result of the first iteration of C _s : p(W ₁ )=Beta((ec) ₁ +||Θ ⁽¹⁾ || ₀ , (fc) ₁ +N ₁ -||Θ ⁽¹⁾ || ₀ ), in: The value set of the subscript s is {1,...,4}. Note that W0 _s and W1 _s exist only when the input S ≥ 2 and the currently considered layer number s is greater than 1. represents the Frobenius norm of the matrix, represents the 0-norm of a matrix (or vector), Θ ^(s) is a matrix formed by combining the wavelet transform coefficients of the sth layer in natural order, R _s is the number of rows in Θ ^s , represents the rth row of Θ ^s , Θ _s,k represents the k-th wavelet coefficient of the s-th layer, Θ _π(s,k) is the coefficient value of the parent node of Θ _s,k , p(x)=Gamma(a,b) means that the random variable x obeys the gamma distribution with parameters (a,b), and its expansion form is where Γ( ) represents the gamma function. Taking a as the independent variable as an example, the expanded form of the gamma function is: p(x)=Beta(a,b) means that the random variable x obeys the beta distribution with parameters (a,b), and its expansion form is Π(·) represents the indicative function, that is, when the expression in parentheses is true, the function value takes 1, otherwise the function value takes 0, According to the posterior probability density distribution obeyed by the above system model parameters α _n , α _s , W ₁ , W0 _s and W1 _s , the system model parameters are sampled, and the corresponding system model parameter values are updated for the next time. Iteration process, Step (3.1.5), store the sampling value of the matrix Θ obtained in step (3.1.3), the first iteration process officially ends, In step (3.2), according to the description of steps (3.1.1) to (3.1.5), the calculation of the second to T _max iterations is performed, Step (3.3), for the serial numbers of all iterations, t={1,...,T _max }, calculate the arithmetic mean of all matrices Θ stored in each iteration, and the obtained mean is denoted as Θ ^* , Perform the inverse two-dimensional wavelet transform of the S layer on the mean Θ ^* to obtain X ^* , Step (3.4), count the error index, PSNR and MSE of the reconstructed image X ^* relative to the original image X ₀ , if these two indices meet the performance requirements of the input, the algorithm stops, otherwise adjust the input in step (1). and let S:=S+1, jump to step (2) again, run the initialization and iterative calculation process again, until the error index of the reconstructed image relative to the original image meets the performance requirements.