Terahertz hyperspectral image restoration method and system based on low-rank joint sparsity

Abstract

The invention discloses a terahertz hyperspectral image recovery method and system based on low-rank joint sparsity, and mainly solves the problems that in the prior art, a terahertz hyperspectral image recovery method adopts a single primary function for representation, detail loss and artifacts are easily caused, and the imaging quality is poor in a low sampling rate or low signal-to-noise ratio environment. According to the implementation scheme, the method comprises the following steps: sparsely acquiring a terahertz hyperspectral image of an image domain to obtain actual observation data; defining a sparse forward transformation operator and a sparse inverse transformation operator by using multiple groups of wavelet bases and Dirac bases; the method comprises the following steps: constructing a terahertz hyperspectral image recovery model based on low-rank joint sparse constraint by using actual observation data and a sparse forward transformation operator, and solving the image recovery model by using a primal-dual splitting method to obtain a recovered terahertz hyperspectral image. The terahertz hyperspectral image recovery method realizes high-resolution terahertz hyperspectral image recovery in an extremely low sampling rate or low signal-to-noise ratio environment, and can be used for medical auxiliary diagnosis, security check screening and drug detection.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, and specifically relates to a terahertz hyperspectral image restoration method, which can be used for medical auxiliary diagnosis, security screening and drug detection. Background Technology

[0002] Terahertz imaging technology, due to its unique fingerprint spectral recognition capabilities and penetration into non-polar materials, holds significant value in fields such as biomedicine, security inspection, and drug analysis. Traditional terahertz imaging methods acquire images through mechanical scanning, which involves moving the object under test using a high-precision stage or shifting the focus of the terahertz beam using optical lenses to traverse every spatial position within the imaging area. The system can only acquire the complete time-domain waveform data of one pixel at a time. This operating mode results in extremely low data acquisition efficiency; acquiring a single high-resolution hyperspectral image often takes several hours or even longer, failing to meet real-time requirements and severely limiting the practical applications of terahertz hyperspectral imaging.

[0003] In their paper "Time-domain terahertz compressive imaging" (OpticsExpress, 2020, 28(3): 3795-3802), Zanotto et al. proposed a single-pixel compressed imaging method based on a time-domain terahertz spectral system. This method utilizes a digital micromirror device (DMD) to spatially encode the beam to modulate the terahertz wavefront. The modulated terahertz signal is collected by a single-point detector, and a multi-dimensional terahertz image is reconstructed using a compressed sensing algorithm. However, this hardware-encoded compressed sensing method relies on spatial modulation devices and cannot be directly applied to existing terahertz time-domain spectral systems with scanning image modes. Furthermore, it leads to blurred image details and texture loss under low sampling rates or low signal-to-noise ratios, requiring a sampling rate of at least 40% to ensure image quality.

[0004] Patent document CN201610829573.4 discloses a terahertz time-domain spectral sparse imaging method. This method utilizes an exponential function to perform a nonlinear transformation on wavelet coefficients to enhance their sparsity, and employs a split-iteration algorithm to decompose the optimization cost function into multiple sub-problems for rapid solution, thereby obtaining the contour features of the imaged object at a low sampling rate of 10%. However, because this method primarily selects the peak values ​​of the terahertz time-domain waveform to construct a two-dimensional imaging model, it loses the frequency domain information of the terahertz wave, thus it cannot be directly applied to the sparse restoration of terahertz hyperspectral images. Furthermore, since this method relies only on a single time-invariant wavelet basis, its reconstruction accuracy and detail capability are low when facing targets of different scales and shapes in complex scenes.

[0005] Patent document CN201610805487.X discloses a hyperspectral image restoration method based on sparse and low-rank matrix approximation. This method divides multi-band hyperspectral data into several small three-dimensional data blocks, then assembles them into two-dimensional data matrices. A weighted Schatten-p paradigm low-rank matrix approximation model is constructed, and the extended Lagrange multiplier method is used for iterative solution to remove mixed noise and restore the multi-band hyperspectral image. However, this method requires iterative singular value decomposition (SVD) of a large number of matrix blocks during the solution process, resulting in extremely high computational complexity and excessively long image restoration time, making it difficult to meet the real-time processing requirements of rapid imaging. Summary of the Invention

[0006] The purpose of this invention is to address the shortcomings of the prior art by proposing a terahertz hyperspectral image restoration method and system based on low-rank joint sparsity, so as to achieve efficient and high-resolution terahertz hyperspectral image restoration in environments with extremely low sampling rates and low signal-to-noise ratios. This effectively avoids detail loss and artifacts when processing targets with complex textures, multiple scales, and multiple morphologies, and significantly improves the universality and robustness of image restoration.

[0007] The technical approach to achieving the objective of this invention is to construct a terahertz hyperspectral image restoration model based on low-rank joint sparse constraints, and to introduce a linear extrapolation primal-dual splitting method for solving the model. This enables efficient and high-resolution terahertz hyperspectral image restoration under extremely low sampling rates and low signal-to-noise ratios. Furthermore, by constructing a sparse transform operator with a joint dictionary definition that includes multiple sets of wavelet bases and Dirac bases, the loss of details and artifacts are effectively avoided when processing targets with complex textures, multiple scales, and multiple morphologies, thereby enhancing the universality and robustness of this invention.

[0008] Based on the above ideas, the technical solution of the present invention includes:

[0009] 1. A terahertz hyperspectral image restoration method based on low-rank joint sparsity, characterized in that it includes:

[0010] (1) Acquire terahertz hyperspectral images in the sparse image domain to obtain actual observation data. ;

[0011] (2) Using multiple sets of wavelet bases and Dirac bases, define a sparse positive transform operator. and sparse inverse transform operator ;

[0012] (3) Utilizing actual observation data Sparse positive transformation operator A terahertz hyperspectral image restoration model based on low-rank joint sparse constraints is constructed:

[0013] ;

[0014] ;

[0015] in, For the terahertz hyperspectral image to be recovered, For sparse sampling mask, For the regularization parameters of the kernel norm, for The regularization parameter of the norm term, For error tolerance, Represents the nuclear norm. express Norm, express Norm, It represents the Hadamardi (or Hadama) stack;

[0016] (4) Solve the above-mentioned terahertz hyperspectral image restoration model based on low-rank joint sparse constraints to obtain the restored terahertz hyperspectral image.

[0017] 2. A terahertz hyperspectral image restoration system based on low-rank joint sparsity, characterized in that it comprises:

[0018] The image domain sparse acquisition module is used to acquire actual observation data. ;

[0019] The sparse transform operator module is used to define sparse positive transform operators using multiple sets of wavelet bases and Dirac bases. and sparse inverse transform operator ;

[0020] The model building module is used to build a terahertz hyperspectral image restoration model based on low-rank joint sparse constraints.

[0021] The model solving module is used to solve the terahertz hyperspectral image reconstruction model to obtain the reconstructed terahertz hyperspectral image. .

[0022] Furthermore, the model building module includes:

[0023] The data fidelity constraint submodule is used to constrain data fidelity residuals. Norm limit within error tolerance Inside, ensure image recovery With observation data As close as possible;

[0024] The nonnegativity constraint submodule is used for reconstructing the image. Apply nonnegativity constraints to satisfy Physical constraints on images in the real number domain;

[0025] The image restoration model construction submodule is used to solve terahertz hyperspectral images. The sparse recovery problem is formulated as an optimization problem of low rank and joint sparsity.

[0026] Compared with the prior art, the present invention has the following advantages:

[0027] 1. By utilizing multiple sets of wavelet bases and Dirac bases to define sparse transform operators, this invention can overcome the loss of details and artifacts caused by single basis function representation, thereby improving the accuracy and robustness of image restoration under complex textures and multi-scale targets.

[0028] 2. This invention utilizes the low rank of terahertz hyperspectral images in the spectral dimension and the joint sparsity in the transform domain to construct a terahertz hyperspectral image restoration model based on low-rank joint sparsity constraints. This model can not only be directly applied to existing terahertz time-domain spectral systems with scanning image modes, but also achieve high-resolution terahertz hyperspectral image restoration under extremely low sampling rates and low signal-to-noise ratio environments.

[0029] 3. This invention employs a linear extrapolation primal-dual splitting method to solve the terahertz hyperspectral image restoration model based on low-rank joint sparse constraints. By introducing dual variables, it effectively separates each regularization term and obtains explicit analytical solutions to subproblems. Therefore, it eliminates the need for internal iteration to solve subproblems, requiring only one iteration layer. This simplifies the computation process, significantly shortens the restoration time, and improves the applicability of terahertz hyperspectral imaging in rapid screening and real-time processing scenarios. Attached Figure Description

[0030] Figure 1 This is a flowchart illustrating the terahertz hyperspectral image restoration method provided in Embodiment 1 of the present invention.

[0031] Figure 2 This is a block diagram of the terahertz hyperspectral image restoration system module provided in Embodiment 2 of the present invention;

[0032] Figure 3 These are simulation diagrams of the recovery results of this invention at different sampling rates;

[0033] Figure 4 These are simulation diagrams of the recovery results of this invention under different signal-to-noise ratios;

[0034] Figure 5 The image restoration simulation diagrams of the present invention, the traditional linear interpolation algorithm, and the orthogonal matching pursuit (OMP) algorithm at different sampling rates at 0.43 THz are shown.

[0035] Figure 6The image restoration simulation diagrams at 0.43THz show the results of this invention, traditional linear interpolation algorithms, and orthogonal matching pursuit (OMP) algorithms with different signal-to-noise ratios. Detailed Implementation

[0036] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention and not all embodiments. Based on the embodiments of the present invention, other embodiments obtained by those skilled in the art without creative effort should all fall within the protection scope of the present invention.

[0037] It should be noted that the step numbers in the specification and claims of this invention are only for the purpose of clearly describing the embodiments of this invention and facilitating understanding, and their order is not limited.

[0038] Example 1: Terahertz hyperspectral image restoration method based on low-rank joint sparsity.

[0039] Reference Figure 1 The implementation steps for this example include the following:

[0040] Step 1. Obtain actual observation data .

[0041] By using terahertz time-domain spectrometers to sparsely acquire terahertz hyperspectral image data, actual observation data can be obtained. :

[0042] ;

[0043] in, For the hyperspectral image to be recovered, For sparse sampling mask, Based on actual observation data, It is additive Gaussian noise. It represents the Hadamah accumulation. This represents the number of pixels in each frame of a hyperspectral image. Indicates the number of frames.

[0044] Step 2. Define the sparse positive transformation operator and sparse inverse transform operator .

[0045] Since terahertz images are sparse in wavelet bases and terahertz hyperspectral images exhibit joint sparsity in the wavelet domain, the structural similarity of real signals across different spectra can be used to preserve effective information, remove independently distributed noise, and improve image restoration accuracy. However, because a single wavelet basis is limited by its fixed vanishing moment and support length, the restoration results become unstable when dealing with multi-scale and multi-morphological targets in complex scenes, making it difficult to meet the sparsity requirements of different structures. Therefore, this example uses multiple sets of wavelet bases and Dirac bases to define a sparse positive transform operator. and sparse inverse transform operator , means as follows:

[0046] ,

[0047] ,

[0048] in, Indicates the first A positive transformation basis, Indicates the first An inverse transform basis, From 1 to , The number of bases.

[0049] Step 3. Utilize actual observation data Sparse positive transformation operator A terahertz hyperspectral image restoration model based on low-rank joint sparse constraints was constructed.

[0050] 3.1) To ensure image recovery With observation data Get as close as possible to the data residuals to ensure data fidelity. of Norm, limited to the error tolerance Inside, it is represented as follows:

[0051] ;

[0052] 3.2) In order to meet the requirements of the terahertz hyperspectral image to be recovered Due to the physical limitations of images in the real number domain, Apply nonnegativity constraints, i.e. ;

[0053] 3.3) Under the constraints of steps 3.1) and 3.2), the solution will be... The sparse recovery problem is formulated as a low-rank and jointly sparse optimization problem, and the nuclear norm is used to characterize the low-rank property of terahertz hyperspectral images in the spectral dimension. The norm characterizes the joint sparsity of terahertz hyperspectral images in the transform domain, leading to a terahertz hyperspectral image reconstruction model based on low-rank joint sparsity constraints:

[0054] ,

[0055] ,

[0056] in, For the terahertz hyperspectral image to be recovered, For sparse sampling mask, For the regularization parameters of the kernel norm, for The regularization parameter of the norm term, For error tolerance, Represents the nuclear norm. express Norm, express Norm, It represents the Hadamardi (or Hadama) stack.

[0057] Step 4. Solve the above terahertz hyperspectral image restoration model based on low-rank joint sparse constraints to obtain the restored terahertz hyperspectral image.

[0058] Existing methods for solving terahertz hyperspectral image restoration models include the Alternating Direction Multiplier Method (ADMM), the Fast Iterative Shrinking Thresholding Algorithm (FISTA), the Proximal Alternating Linearization Algorithm (PALM), and the Primitive-Dual Split-Solve Method (PDS). Based on the performance of each algorithm, this example adopts, but is not limited to, the Primitive-Dual Split-Solve Method (PDS) for solution. Its specific implementation includes the following:

[0059] 4.1) Introducing characteristic functions transforms the constrained optimization problem into the following unconstrained optimization problem:

[0060] ,

[0061] in, Indicated by observation data Centered on the ball, with error tolerance Indicator function of a set of spheres with radius , Represents the characteristic function of the nonnegative real number field;

[0062] 4.2) Solve the unconstrained optimization problem in step 4.1) to obtain the recovered terahertz hyperspectral image. :

[0063] 4.2.1) Initialize all parameters:

[0064] Let the regularization parameter of the nuclear norm be... , Norm term regularization parameter ,

[0065] Let the original variable iteration step size be... The iteration step size of the dual variable , , ,

[0066] Set the initial values ​​of the hyperspectral image to be recovered. Initial value of extrapolated variable ,

[0067] Let the initial values ​​of the dual variables of the nuclear norm term be... , Initial values ​​of the dual variable of the norm term The initial value of the dual variable of the data fidelity term ,

[0068] Set error tolerance Outer layer residual tolerance The relative change limit for stopping iterations Maximum number of iterations ;

[0069] 4.2.2) Utilizing the dual variable of the nuclear norm term , Dual variable of norm term The dual variable of the data fidelity term Original variables Calculate the original variable after the current update. :

[0070] ,

[0071] in, This represents the sparse inverse transform operator. This represents the transpose of the sparse sampling mask; It is a projection operator for the characteristic function of the non-negative real number field, which means projecting the input variable onto the non-negative real number field.

[0072] 4.2.3) Utilizing the updated original variables Apply linear extrapolation to calculate and update the extrapolated variables for the current iteration. :

[0073] ;

[0074] 4.2.4) Utilizing updated extrapolated variables Update the dual variable of the nuclear norm term. :

[0075] ,

[0076] in, The iteration step size represents the dual variable of the nuclear norm term; This represents the singular value soft thresholding operation, which is used to preserve input matrices larger than or equal to 0. Set the singular values ​​of the singular values ​​to zero, and set the rest of the singular values ​​to zero.

[0077] 4.2.5) Utilizing updated extrapolated variables ,renew Dual variable of norm term :

[0078] ,

[0079] in, This represents the sparse positive transformation operator. express The iteration step size of the dual variable of the norm term.

[0080] yes The proximal operator of the norm, representing the expression for each row of the input matrix. The norm is used for soft threshold shrinkage, which is used to preserve... norm greater than Set the first row to zero and the rest of the rows to zero.

[0081] 4.2.6) Utilizing updated extrapolated variables Update the dual variable of the data fidelity term. :

[0082] ,

[0083] in, The iteration step size represents the dual variable of the data fidelity term. Indicates a sparse sampling mask. The projection operation, representing the data fidelity term, is used to project input variables onto a matrix. Centered on the ball, with error tolerance On a sphere with radius [radius];

[0084] 4.2.7) Utilize the original variables from before the current iteration update and the original variables after the current iteration update Calculate the original variables relative change :

[0085] ;

[0086] 4.2.8) For the data fidelity residuals in step 3.1), of Norms and relative changes in step 4.2.7) Make a judgment:

[0087] If relative change is satisfied Or data fidelity residuals norm If any of the following conditions are met, and the number of iterations is satisfied, then... If the current iteration number is incremented by 1, proceed to the next iteration, and return to step 4.2.2), update the original variable for the next iteration;

[0088] If the conditions are not met, the algorithm iteration ends, and the recovered terahertz hyperspectral image is obtained. .

[0089] Example 2: Terahertz hyperspectral image restoration system based on low-rank joint sparsity.

[0090] Reference Figure 2 This example includes: an image domain sparse acquisition module 1, a sparse transformation operator module 2, a model building module 3, and a model solving module 4. The model building module 3 includes a data fidelity constraint submodule 31, a nonnegativity constraint submodule 32, and an image restoration model building submodule 33. The model solving module 4 includes an unconstrained form transformation submodule 41, a parameter initialization submodule 42, an original variable update submodule 43, an extrapolated variable update submodule 44, and a nuclear norm term dual variable update submodule 45. Norm term dual variable update submodule 46, data fidelity term dual variable update submodule 47, original variable The relative change calculation submodule 48 and the iterative convergence judgment submodule 49.

[0091] The working principle of the entire system is as follows:

[0092] The image domain sparse acquisition module 1 is used to acquire actual observation data. and the actual observation data Transmitted to model building module 3;

[0093] The sparse transform operator module 2 is used to define sparse positive transform operators using multiple sets of wavelet bases and Dirac bases. and sparse inverse transform operator The sparse forward transform operator and the sparse inverse transform operator are then transmitted to the model building module 3.

[0094] The model building module 3 is used to build the model based on the actual observation data obtained by the image domain sparse acquisition module 1. The sparse positive transform operator defined in sparse transform operator module 2 constructs a terahertz hyperspectral image restoration model based on low-rank joint sparse constraints. The data fidelity constraint submodule 31 is used to convert the data fidelity residuals... Norm limit within error tolerance Inside, ensure image recovery The actual observation data obtained in the image domain sparse acquisition module 1 As close as possible, and transmit this data fidelity constraint to the image restoration model construction submodule 33, to obtain the data fidelity residuals. The norm is passed to the iterative convergence judgment submodule 49; the nonnegativity constraint submodule 32 is used to process the recovered image. Apply nonnegativity constraints to satisfy The physical constraints of the real-domain image are applied, and this non-negativity constraint is transmitted to the image restoration model construction submodule 33. Under the data fidelity constraints in the data fidelity constraint submodule 31 and the non-negativity constraints in the non-negativity constraint submodule 32, the image restoration model construction submodule 33 will solve for the terahertz hyperspectral image. The sparse restoration problem is formulated as a low-rank and jointly sparse optimization problem; finally, the image restoration model is transmitted to the model solving module 4.

[0095] The model solving module 4 is used to solve the image restoration model constructed by the model building module 3 to obtain the restored terahertz hyperspectral image. ,in:

[0096] The unconstrained form transformation submodule 41 is used to transform the constrained optimization problem in the model building module 3 into an unconstrained optimization problem by introducing an indicator function, thereby obtaining an unconstrained optimization model, and then transmitting the unconstrained optimization model to the parameter initialization submodule 42.

[0097] The parameter initialization submodule 42 is used to initialize various parameters based on the unconstrained form optimization model obtained in the unconstrained form transformation submodule 41, and then transmit them to the original variable update submodule 43, the extrapolation variable update submodule 44, and the nuclear norm dual variable update submodule 45, respectively. Norm term dual variable update submodule 46, data fidelity term dual variable update submodule 47, original variable The relative change calculation submodule 48 and the iterative convergence judgment submodule 49;

[0098] The original variable update submodule 43 is used to utilize the dual variable of the nuclear norm term. , Dual variable of norm term The dual variable of the data fidelity term Original variables Calculate the original variables after the current iteration update. and will Transmitted to extrapolation variable update submodule 44 and original variable The relative change calculation submodule 48;

[0099] Extrapolation variable update submodule 44 is used to update the original variables after iterative updates in submodule 43 using the original variables. Apply linear extrapolation to calculate the extrapolated variables after the current iteration update. and will Transmitted to the kernel norm dual variable update submodule 45 Norm term dual variable update submodule 46 and data fidelity term dual variable update submodule 47;

[0100] Kernel norm dual variable update submodule 45 is used to update the extrapolated variables iteratively updated in submodule 44 using extrapolated variables. Update the dual variable of the nuclear norm term. ;

[0101] Norm term dual variable update submodule 46 is used to update the extrapolated variables iteratively updated in submodule 44 using extrapolated variables. ,renew Dual variable of norm term ;

[0102] Data fidelity term dual variable update submodule 47 is used to update the extrapolated variables iteratively updated in extrapolated variable update submodule 44. Update the dual variable of the data fidelity term. ;

[0103] primitive variables The relative change calculation submodule 48 is used to utilize the original variables before the current iteration update. The original variables updated iteratively in the original variable update submodule 43 Calculate the original variables relative change and this relative change Transmitted to the iterative convergence judgment submodule 49;

[0104] Iterative convergence judgment submodule 49 utilizes the data fidelity residuals from data fidelity constraint submodule 31. Norm, primitive variables The relative change calculation submodule 48 Error tolerance in parameter initialization submodule 42 Outer layer residual tolerance The relative change limit for stopping iterations and maximum number of iterations To determine whether the termination condition of the iteration is met, if not, increment the current iteration number by 1 and return to the original variable update submodule 43 to continue the iteration; if the condition is met, stop the iteration and obtain the recovered terahertz hyperspectral image.

[0105] It should be noted that the above functional modules can be implemented, in whole or in part, through software, hardware, firmware, or any combination thereof. When implemented in software, they can be implemented, in whole or in part, as a program instruction product. A program instruction product includes one or a set of program instructions. When the program instructions are loaded and executed on a computer, the described process or function is generated, in whole or in part. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The program instructions can be stored in a computer-readable and writable storage medium, or transferred from one computer's readable and writable storage medium to another.

[0106] In this embodiment, the direct coupling or communication connection between the modules can be achieved through indirect coupling or communication connection via interfaces, devices, or modules. The functional modules and sub-modules in this embodiment can dynamically reside within a single processing unit, or each module can exist physically independently, or two or more modules can dynamically reside within a single processing unit. When these dynamic components are implemented as software functional modules and sold or used as independent products, they can also be stored in a computer-readable and writable storage medium. This storage medium can be a memory, disk, or optical disc, etc.

[0107] The effectiveness of this invention can be further demonstrated through the following simulation.

[0108] 1. Simulation experimental conditions.

[0109] The software platform for the simulation experiment is Windows 11 operating system and MATLAB R2024b.

[0110] The experimental data used in the simulation experiment consisted of measured data from plastic objects embedded with the letters "HT," with an object size of approximately 40mm x 40mm. A sample set was compiled using images in the frequency band from 0.43THz to 1.18THz, with each image measuring 256 x 256 pixels.

[0111] The control group data in the simulation experiment is the complete terahertz hyperspectral image data obtained by scanning the plastic embedded with the letter HT point by point using a terahertz time-domain spectrometer. The original image of the control group can be obtained by directly visualizing the data using MATLAB R2024b.

[0112] 2. Simulation experiment content and result analysis.

[0113] Simulation Experiment 1: Using the method of this invention, image restoration was performed on the above experimental data at different sampling rates and four frequency points: 0.43THz, 0.68THz, 0.93THz, and 1.18THz. The observed data images, restored images, and original images under the corresponding conditions were obtained, as shown below. Figure 3 As shown, where:

[0114] Figure 3 (a) The result of image restoration at the four frequency points with a 10% sampling rate;

[0115] Figure 3 (b) The results of image restoration at the four frequency points with a sampling rate of 20%;

[0116] Figure 3 (c) The result of image restoration at the four frequency points with a sampling rate of 30%;

[0117] Figure 3 (d) shows the results of image restoration at the four frequency points with a sampling rate of 50%.

[0118] Depend on Figure 3 It is evident that the method of this invention can achieve high-quality image restoration even at a sampling rate of 10%.

[0119] Simulation Experiment 2: Using the method of this invention, image restoration was performed on the above experimental data under different signal-to-noise ratios and four frequency points (0.43THz, 0.68THz, 0.93THz, and 1.18THz). The observed data images, restored images, and original images under the corresponding conditions were obtained, as shown below. Figure 4 As shown, where:

[0120] Figure 4 (a) is the result of image restoration at the four frequency points with a signal-to-noise ratio of 10;

[0121] Figure 4 (b) The result of image restoration at the four frequency points with a signal-to-noise ratio of 20;

[0122] Figure 4 (c) The result of image restoration at the four frequency points with a signal-to-noise ratio of 30;

[0123] Figure 4 (d) shows the results of image restoration at the four frequency points with a signal-to-noise ratio of 50.

[0124] Depend on Figure 4As can be seen, the method of the present invention can recover the image well when the signal-to-noise ratio is 20, indicating that the method of the present invention has a good ability to suppress noise. In terahertz experimental systems, the signal-to-noise ratio is generally higher than 40 dB, so the method of the present invention has good applicability.

[0125] Simulation Experiment 3 involves using the present invention and existing linear interpolation algorithms Lerp and OMP (Orthogonal Matching Pursuit) to perform image restoration on the above experimental data at 0.43 THz under sampling rates of 0.1, 0.3, and 0.5, respectively. The results are as follows: Figure 5 As shown.

[0126] Depend on Figure 5 As can be seen, linear interpolation can reconstruct the contour of an image, but it will result in lines in the middle and some image loss. The higher the sampling data rate, the higher the image quality. The orthogonal matching pursuit algorithm only has a good restoration effect when the sampling data rate is greater than or equal to 50%, while the method of this invention can restore the image with high quality even at a lower sampling data rate, such as 10%, verifying that the method of this invention can restore the image with high quality even at extremely low sampling rates.

[0127] Simulation Experiment 4: Using the present invention and existing linear interpolation algorithms Lerp and OMP, image restoration was performed on the above experimental data at 0.43 THz under signal-to-noise ratios of 10, 20, 30, 40, and 50. The results are as follows. Figure 6 As shown.

[0128] Depend on Figure 6 As can be seen, linear interpolation still results in a linear pattern even at a signal-to-noise ratio (SNR) as high as 50 dB, leading to poor image restoration quality. The orthogonal matching pursuit algorithm only achieves good restoration results when the SNR is greater than or equal to 50 dB, while the method of this invention can restore images well even at an SNR greater than or equal to 20 dB, verifying that the method of this invention has excellent noise resistance and can restore images with high quality even at low SNR.

[0129] The simulation results above show that the present invention can achieve high-resolution terahertz hyperspectral image recovery under extremely low sampling rate and low signal-to-noise ratio conditions.

Claims

1. A terahertz hyperspectral image restoration method based on low-rank joint sparsity, characterized in that, include: (1) Acquire terahertz hyperspectral images in the sparse image domain to obtain actual observation data. ; (2) Using multiple sets of wavelet bases and Dirac bases, define a sparse positive transform operator. and sparse inverse transform operator ; (3) Utilizing actual observation data Sparse positive transformation operator A terahertz hyperspectral image restoration model based on low-rank joint sparse constraints is constructed: ； ； in, For the terahertz hyperspectral image to be recovered, For sparse sampling mask, For the regularization parameters of the kernel norm, for The regularization parameter of the norm term, For error tolerance, Represents the nuclear norm. express Norm, express Norm, It represents the Hadamardi (or Hadama) stack; (4) Solve the above-mentioned terahertz hyperspectral image restoration model based on low-rank joint sparse constraints to obtain the restored terahertz hyperspectral image.

2. The method according to claim 1, characterized in that, The actual observation data obtained in (1) , means as follows: ； in, For the hyperspectral image to be recovered, For sparse sampling mask, Based on actual observation data, It is additive Gaussian noise. This represents the number of pixels in each frame of a hyperspectral image. Indicates the number of frames.

3. The method according to claim 1, characterized in that, In (2), multiple sets of wavelet bases and Dirac bases are used to define a sparse positive transform operator. and sparse inverse transform operator It is represented as follows: ； ； in, The number of bases, Indicates the first A positive transformation basis, Indicates the first A 5 inverse transform basis.

4. The method according to claim 1, characterized in that, The terahertz hyperspectral image restoration model constructed in (3) based on low-rank joint sparse constraints includes the following implementation: (3a) To ensure image recovery With observation data As close as possible to the residual Norm, limited to the error tolerance within, that is ; (3b) In order to satisfy Due to the physical limitations of images in the real number domain, Apply nonnegativity constraints, i.e. ; (3c) Under the constraints of (3a) and (3b), the solution will be... The sparse restoration problem is formulated as a low-rank and jointly sparse optimization problem, resulting in a terahertz hyperspectral image restoration model based on low-rank joint sparse constraints: ； 。 5. The method according to claim 1, characterized in that, The solution to the terahertz hyperspectral image restoration model based on low-rank joint sparse constraints in (4) includes: (4a) Introducing the characteristic function transforms the above constrained optimization problem into the following unconstrained optimization problem: ， in, Indicated by observation data Centered on the ball, with error tolerance Indicator function of a set of spheres with radius , Represents the characteristic function of the nonnegative real number field; (4b) Solve the unconstrained optimization problem in (4a) to obtain the recovered terahertz hyperspectral image. .

6. The method according to claim 5, characterized in that, The solution to the unconstrained optimization problem in (4b) is implemented as follows: (4b1) Initialize all parameters: Kernel norm regularization parameter , Norm term regularization parameter , Original variable iteration step size The iteration step size of the dual variable , , , Initial values ​​of the hyperspectral image to be recovered Initial value of extrapolated variable , Initial values ​​of the dual variables of the nuclear norm term , Initial values ​​of the dual variable of the norm term The initial value of the dual variable of the data fidelity term , Error tolerance Outer layer residual tolerance The relative change limit for stopping iterations Maximum number of iterations ; (4b2) Using the dual variable of the nuclear norm term , Dual variable of norm term The dual variable of the data fidelity term , original variables Calculate the original variable after the current update. : ， in, This represents the sparse inverse transform operator. This represents the transpose of the sparse sampling mask. It is a projection operator for the characteristic function of the non-negative real number field, which means projecting the input variable onto the non-negative real number field. (4b3) Using the updated original variables Apply linear extrapolation to calculate and update the extrapolated variables for the current iteration. : ； (4b4) Using updated extrapolated variables Update the dual variable of the nuclear norm term. : ， in, The iteration step size represents the dual variable of the nuclear norm term; This represents the singular value soft thresholding operation, which is used to preserve input matrices larger than or equal to 0. Set the singular values ​​of the singular values ​​to zero, and set the rest of the singular values ​​to zero. (4b5) Using updated extrapolated variables ,renew Dual variable of norm term : ， in, This represents the sparse positive transformation operator. express The iteration step size of the dual variable of the norm term. yes The proximal operator of the norm, representing the expression for each row of the input matrix. The norm is used for soft threshold shrinkage, which is used to preserve... norm greater than Set the first row to zero and the rest of the rows to zero. (4b6) Using updated extrapolated variables Update the dual variable of the data fidelity term. : ， in, The iteration step size represents the dual variable of the data fidelity term. Indicates a sparse sampling mask. The projection operation, representing the data fidelity term, is used to project input variables onto a matrix. Centered on the ball, with error tolerance On a sphere with radius [radius]; (4b7) Use the original variables from before the current iteration update and the original variables after the current iteration update Calculate the original variables relative change : ， (4b8) Determine whether the relative change condition is satisfied. Or data fidelity residuals If any of the following conditions are met, and the number of iterations is satisfied, then... : If satisfied, increment the current iteration number by 1 and proceed to the next iteration, return to step (4b2), and update the original variable for the next iteration; If the conditions are not met, the algorithm iteration ends, and the recovered terahertz hyperspectral image is obtained. .

7. A terahertz hyperspectral image restoration system based on low-rank joint sparsity, characterized in that, include: The image domain sparse acquisition module is used to acquire actual observation data. ; The sparse transform operator module is used to define sparse positive transform operators using multiple sets of wavelet bases and Dirac bases. and sparse inverse transform operator ; The model building module is used to build a terahertz hyperspectral image restoration model based on low-rank joint sparse constraints. The model solving module is used to solve the terahertz hyperspectral image reconstruction model to obtain the reconstructed terahertz hyperspectral image. .

8. The system according to claim 7, characterized in that, The model building module includes: The data fidelity constraint submodule is used to limit the data fidelity items within a specific error range to ensure image recovery. With observation data As close as possible; The nonnegativity constraint submodule is used for reconstructing the image. Apply nonnegativity constraints to satisfy Physical constraints on images in the real number domain; The image restoration model construction submodule is used to solve terahertz hyperspectral images. The sparse recovery problem is formulated as an optimization problem of low rank and joint sparsity.

9. The system according to claim 7, characterized in that, The model solving module includes: The Unconstrained Form Transformation submodule is used to transform constrained optimization problems into unconstrained optimization problems by introducing characteristic functions; The parameter initialization submodule is used to initialize the various parameters in the algorithm. The original variable update submodule is used to utilize the dual variable of the nuclear norm term. , Dual variable of norm term The dual variable of the data fidelity term , original variables To obtain the updated original variables ; The extrapolation variable update submodule is used to update the original variable. Apply linear extrapolation to calculate the extrapolated variables after the current iteration update. ; The nuclear norm dual variable update submodule is used to utilize the updated extrapolated variables. Update the dual variable of the nuclear norm term. ; The norm term dual variable update submodule is used to utilize the updated extrapolated variables. ,renew Dual variable of norm term ; The data fidelity term dual variable update submodule is used to utilize the updated extrapolated variables. Update the dual variable of the data fidelity term. ; primitive variables The relative change calculation submodule is used to utilize the original variables before the current iteration update. and the original variables after the current iteration update Calculate the original variables relative change ; The iterative convergence judgment submodule is used to determine whether the termination condition of the iteration is met. If it is met, the iteration stops and the recovered terahertz hyperspectral image is obtained; if it is not met, the current iteration number is incremented by 1 and the original variable update submodule is returned to continue the iteration.

Images