A spectral imaging method based on high-order mathematical modeling and fitting calibration

Abstract

The present invention relates to a spectral imaging method based on high-order mathematical modeling and fitting calibration, belonging to the technical field of spectral imaging. An equivalent correction code is set in the established high-order mathematical model, and while implementing the encoding operation, it characterizes the response non-uniformity. Before and after the equivalent correction code of the established high-order mathematical model, a group or multiple groups of two-dimensional convolution kernels are set to comprehensively characterize non-ideal factors such as aberration, pixel mismatch, high-order dispersion, and assembly error, improving the coincidence degree between the mathematical model and the actual projection measurement process. The calibration process of the high-order mathematical model parameters is constructed as an inverse optimization problem, training sample data is collected, and an optimization algorithm is used to solve the constructed optimization problem to complete the calibration of the high-order mathematical model parameters, improving the practicability of the calibration method. After the spectral data cube is modulated by the coded aperture and the dispersion prism, it is imaged onto the detector. Finally, based on the established high-order mathematical model, a reconstruction algorithm is used to reconstruct the spectral image.

Description

Technical Field

[0001] The present invention relates to a spectral imaging method based on high - order mathematical modeling and fitting calibration, belonging to the technical field of spectral imaging. Background Technique

[0002] Spectral imaging technology can obtain spatial images of a target scene in different bands, has the ability to identify substances and components non - destructively, and has wide applications in fields such as biochemical analysis, remote sensing monitoring, and precision agriculture. A spectral image is three - dimensional data, also known as a "spectral data cube", and the smallest data unit therein is called a "voxel". Common image detectors (abbreviation: "detector") usually have a planar structure and can only detect two - dimensional data each time they image. Therefore, the core problem of spectral imaging lies in how to use a low - dimensional detector to obtain a three - dimensional spectral image.

[0003] In 2008, based on the theory of compressive sensing, researchers proposed the coded aperture snapshot spectral imaging (CASSI) technology. The CASSI technology performs spatial encoding and dispersion translation on the spectral data cube, detects two - dimensional projection measurement images in which all spectral bands are aliased, and by establishing a mathematical model for the above - mentioned projection measurement process and using a reconstruction algorithm, can reconstruct the spectral image to be detected. The CASSI technology can reconstruct a high - dimensional spectral image from low - dimensional projection data, is expected to solve the problem of mutual restriction of multiple detection performances in traditional spectral imaging technology, and realize the simultaneous improvement of spatial resolution, spectral resolution, time resolution, and detection signal - to - noise ratio.

[0004] The mathematical model of CASSI is used to describe the entire - link imaging process from the target spectral data cube to the detector projection measurement image. The parameters in the mathematical model need to be calibrated according to the actual CASSI system. Since the mathematical model is the basis for the reconstruction algorithm to reconstruct the spectral image, the accuracy of the model has a crucial impact on the imaging accuracy of the CASSI technology. The CASSI technology initially adopted an ideal mathematical model, which assumes that the imaging process of the CASSI system is an ideal point - to - point imaging, and each discrete band with a finite bandwidth in the spectral data cube is approximated as an ideal monochromatic band. In 2013, Arguello et al. studied the limitations of the ideal mathematical model. In general detection scenarios, the optical wavelength is continuously distributed, and the dispersion translation amount caused by the prism changes continuously with the optical wavelength, rather than being misaligned in units of voxels. The dispersion phenomenon exists not only between the bands of the spectral data cube but also within the bands. This non - ideal factor is called "higher - order dispersion". In response, this group of researchers proposed a higher - order dispersion mathematical model, which can characterize the higher - order dispersion effect and improve the accuracy of the CASSI mathematical model. However, the higher - order dispersion model needs to use a method of scanning band - by - band and field - by - field to calibrate the higher - order dispersion parameters, with complex operations and poor practicability.

[0005] In 2015, Galvis et al. pointed out in their research that if the coded aperture pixels and the detector pixels cannot be accurately aligned (referred to as "pixel mismatch"), it will also cause the ideal mathematical model to deviate from the actual detection process, thereby reducing the imaging accuracy of the CASSI system. In response, the researchers in this group used the actual photographed coded pattern to replace the original coded pattern and proposed a mathematical model based on the actual photographed code, which can characterize the pixel mismatch effect. However, the mathematical model based on the actual photographed code needs to re-photograph the actual photographed code and calibrate the pixel mismatch parameters every time a new coded pattern is used, reducing the practicality of the CASSI technology.

[0006] In addition to high-order dispersion and pixel mismatch, there are also non-ideal factors such as response non-uniformity, aberration, and assembly error in the CASSI system. The existing mathematical models cannot effectively characterize these factors. Moreover, the existing mathematical models only consider individual non-ideal factors and cannot comprehensively describe the combined influence of the above-mentioned multiple non-ideal factors on the projection measurement image. In addition, the existing method for calibrating the parameters of the mathematical model requires complex operations and poor practicality, restricting the application fields of the CASSI technology.

[0007] In summary, the mathematical modeling method and the method for calibrating the parameters of the mathematical model in the CASSI technology need to be further improved. Summary of the Invention

[0008] The technical problem to be solved by the present invention is to overcome the deficiencies of the prior art and propose a spectral imaging method based on high-order mathematical modeling and fitting calibration.

[0009] The technical solution of the present invention is as follows:

[0010] A spectral imaging method based on high-order mathematical modeling and fitting calibration. In the high-order mathematical model established in this method, an equivalent correction code is set to characterize the response non-uniformity while performing the coding operation. Before and after the equivalent correction code in the established high-order mathematical model, a set or multiple sets of two-dimensional convolution kernels are set to comprehensively characterize non-ideal factors such as aberration, pixel mismatch, high-order dispersion, and assembly error, improving the degree of coincidence between the mathematical model and the actual projection measurement process. Moreover, the calibration process of the high-order mathematical model parameters is constructed as an inverse optimization problem, training sample data is collected, and an optimization algorithm is used to solve the constructed optimization problem to complete the calibration of the high-order mathematical model parameters, improving the practicality of the calibration method. After the spectral data cube is modulated by the coded aperture and the dispersion prism, it is imaged onto the detector. Finally, based on the established high-order mathematical model, a reconstruction algorithm is used to reconstruct the spectral image. The steps of this method include:

[0011] The first step is to irradiate the spectral data cube to be measured with a light source The spectral data cube X is imaged onto the coded aperture through an imaging lens, encoded by the coding pattern, and represented by the matrix indicating the coding pattern loaded on the coded aperture during the s-th projection measurement. The encoded spectral data cube X is dispersed in the spatial horizontal or vertical direction by a dispersion prism, and the dispersed spectral data cube X is imaged onto an array detector through a relay lens. The array detector integrates the spectral data cube X along the spectral dimension to obtain the two-dimensional projection measurement image in the s-th projection measurement where V = (N + L - 1);

[0012] Second, establish a high-order mathematical model for the projection measurement process in the first step;

[0013] Third, collect the spectral data cube and projection measurement images of the calibration target as training sample data, then construct an optimization problem, solve the optimization problem through a numerical algorithm, fit the training sample data, and calibrate the parameters of the high-order mathematical model established in the second step;

[0014] Fourth, perform vectorized scanning on the two-dimensional projection measurement image Y s to convert the two-dimensional projection measurement image Y s into a column vector Perform vectorized scanning on the spectral data cube X to be measured, and convert the three-dimensional tensor X into a column vector Represent the projection matrix corresponding to the s-th projection measurement process by the matrix According to the high-order mathematical model calibrated in the third step, establish the matrix high-order mathematical model y s = H s x;

[0015] Fifth, change the coding pattern E s , repeat the content described in the first step, perform a total of S projection measurements, and represent all projection measurement images by the vector y = [y 1T , y 2T , …, y ST T where · T represents the transpose operation. Substitute the changed coding patterns into the high-order mathematical model calibrated in the third step in turn, repeat the content described in the fourth step, and establish the matrix high-order mathematical model y = Hx corresponding to a total of S projection measurements, where H = [H 1T , H 2T , …, H ST T ;

[0016] ​​Step 6: Since there is a high correlation between adjacent bands and adjacent spatial positions in the spectral data cube, it is assumed that x is sparsely represented on the sparse basis Ψ, and θ represents the sparse coefficient. According to the principle of compressive sensing, y = Hx is represented by the underdetermined equation y = HΨθ;

[0017] Step 7: Use a reconstruction algorithm to solve the underdetermined equation y = HΨθ established in Step 6 to recover the spectral data cube X.

[0018] In the first step described above, the spectral data cube to be measured is represented by a three-dimensional tensor where M, N, and L represent the sizes of the spectral data cube in the vertical spatial dimension, horizontal spatial dimension, and spectral dimension, respectively;

[0019] In the second step described above, the method for establishing a high-order mathematical model is as follows:

[0020] Step 201: Based on the two-dimensional coding pattern Construct the three-dimensional equivalent correction coding in the s-th projection measurement according to the formula where the three-dimensional tensor represents the response non-uniformity, the subscript · represents the two-dimensional matrix corresponding to the k-th band in the three-dimensional tensor, ⊙ represents the Hadamard product operation, and the equivalent correction coding C k represents the response non-uniformity while performing the coding operation. If there is no obvious response non-uniformity during the projection measurement process, then R is a tensor with all element values being 1, and the equivalent correction coding degenerates into the original coding pattern at this time; s In the meantime of implementing the coding operation, it represents the response non-uniformity. If there is no obvious response non-uniformity during the projection measurement process, then R is a tensor with all element values being 1, and the equivalent correction coding degenerates into the original coding pattern at this time;

[0021] Step 202: Divide the spatial region and band groups in the spectral data cube. Let U, V, and W represent the total number of divisions of the spectral data cube in the vertical spatial dimension, horizontal spatial dimension, and spectral dimension, respectively. Let X u，v，w represent the subset of the spectral data cube obtained by the division. Then X = {X u，v，w}, where u, v, and w represent the indices of the vertical spatial region, horizontal spatial region, and band group, respectively. Each group of index values points to a subset of the data cube, u ∈ {1, 2,..., U}, v ∈ {1, 2,..., V}, w ∈ {1, 2,..., W}. If each band is separately divided into a band group, then this subset is a matrix; otherwise, this subset is a three-dimensional tensor;

[0022] Step 203: Set one or more groups of two-dimensional convolution kernels before and after the equivalent correction coding. Let B represent the number of convolution kernel groups, let W b represent one of the groups of convolution kernels, b ∈ {1, 2,..., B}, and let represent a single two-dimensional convolution kernel in each group of convolution kernels. Then For the data cube subset X u，v，w The corresponding data is subjected to three-dimensional convolution operation, and the convolution kernel comprehensively characterizes various non-ideal factors such as the aberration of the imaging lens group, the aberration of the relay lens group, the pixel mismatch between the coded aperture and the detector, the high-order dispersion of the prism, and the assembly error of the system;

[0023] Step 204, combining the dispersion and integration effects in the projection measurement process, constructing a high-order mathematical model containing equivalent correction codes and one or more groups of two-dimensional convolution kernels. Taking three groups of convolution kernels as an example, the high-order mathematical model is constructed as follows: in It represents the three-dimensional convolution operation of dividing space into regions and bands into small groups. The data corresponding to all data cube subsets are subjected to three-dimensional convolution operation with the two-dimensional convolution kernels of the same index according to the indexes of u, v, and w. If there is a convolution kernel group consisting of only unit impulse functions, the high-order mathematical model of three groups of convolution kernels degenerates into a high-order mathematical model of two groups or one group of convolution kernels. If step 202 divides each band of the data cube into a band group separately, the three-dimensional convolution operation at this time degenerates into a two-dimensional convolution operation. disp(·) represents the dispersion operation. Taking the matrix corresponding to the first band of the three-dimensional tensor as the reference, the matrices corresponding to the remaining bands are respectively translated by one unit in the positive direction of the horizontal or vertical dimension of the space relative to the matrix of the previous band. ∑ k Represents the accumulation operation along the spectral dimension, accumulating the spectral data at each spatial position of the three-dimensional tensor, and the result of the operation is a matrix;

[0024] In the third step, the specific calibration method of the high-order mathematical model parameters is:

[0025] Step 301: Establish a low-order mathematical model Z corresponding to the p-th projection measurement p =∑ k disp(X⊙E p ), where Z p is the projection measurement image calculated by the low-order mathematical model, for the matrix Z p Perform vectorized scanning, convert Zp into column vector zp, use Gp to represent the projection matrix corresponding to the p-th low-order projection measurement process, then the matrix low-order mathematical model of the total P projection measurements is z = Gx = GΨθ, where z = [z 1T , z 2T , …, z PT ] T , G=[G 1T , G 2T , …, G PT ] T , Ψ represents a sparse basis;

[0026] Step 302: Use a standard whiteboard as the sample target, and change the coding pattern E p , and under the illumination of the light source used during imaging, capture P projection measurement images For perform vectorized scanning to obtain the corresponding column vectors According to the formula and Use to represent all P projection measurement images captured, and construct an optimization problem where τ represents the weight factor, ||·|| t represents the calculation of the l t norm;

[0027] Step 303: Use a compressed sensing reconstruction algorithm to solve for r, and arrange r into a three-dimensional tensor Calibrate to obtain the response non-uniformity parameter;

[0028] Step 304: On the basis of the known response non-uniformity parameter R, use an object with certain texture features as the sample target, and change the coding pattern E q , and under the illumination of the light source used during imaging, capture Q projection measurement images If the spectral data cube of the sample target is unknown, use other spectral imaging systems to obtain the spectral data cube of the sample target If the spectral data cube of the sample target is known, use it directly, and then construct an optimization problem where Y q is the projection measurement image calculated by the high-order mathematical model in the q-th projection measurement, and C q is the equivalent correction code in the q-th projection measurement;

[0029] Step 305: Use an optimization algorithm to solve for W 1 , W 2 , …, W B , and calibrate to obtain a set or multiple sets of two-dimensional convolution kernel parameters set

[0030] Beneficial effects

[0031] (1) The high-order mathematical modeling method of the present invention sets an equivalent correction code and a set or multiple sets of two-dimensional convolution kernels, comprehensively characterizing various non-ideal factors such as response non-uniformity, aberration, pixel mismatch, high-order dispersion, and assembly error.

[0032] (2) The high-order mathematical model established by the high-order mathematical modeling method of the present invention has a higher degree of coincidence with the actual projection measurement process compared to the existing mathematical models, thereby improving the reconstruction accuracy of spectral images.

[0033] (3) The fitting and calibration method of the high-order mathematical model of the present invention constructs the calibration process as an inverse optimization problem, collects training sample data, uses an optimization algorithm to solve the constructed optimization problem, and completes the calibration of the parameters of the high-order mathematical model.

[0034] (4) The fitting and calibration method of the high-order mathematical model of the present invention avoids complex scanning operations compared with the existing calibration methods of mathematical models, and does not require re-calibration every time a new coding pattern is used, thus having higher practicability.

[0035] (5) In the method of the present invention, after the spectral data cube is modulated by the coded aperture and the dispersion prism, it is integrated on the detector to form a projection measurement image. Then, based on the established high-order mathematical model and the detected projection measurement image, a reconstruction algorithm is used to reconstruct the spectral image, which has a high detection efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 is a schematic flow chart of the method of the present invention;

[0037] Figure 2 is a structural diagram of the CASSI system applied in the present invention;

[0038] Figure 3 is a comparison of the reconstruction results of the high-order mathematical model and the fitting calibration method, the mathematical model based on real-shot coding and its calibration method under the same number of projection measurements and the same reconstruction algorithm used;

[0039] Figure 4 is a sampling comparison of the reconstructed spectral images of two mathematical models and their calibration methods. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0040] The present invention will be further described in detail below with reference to the accompanying drawings.

[0041] The object of the present invention is to provide a high-order mathematical modeling and fitting calibration method. The established high-order mathematical model sets an equivalent correction code, which characterizes the response non-uniformity while implementing the coding operation. Before and after the equivalent correction code of the established high-order mathematical model, a group or multiple groups of two-dimensional convolution kernels are set to comprehensively characterize non-ideal factors such as aberration, pixel mismatch, high-order dispersion, and assembly error. And, the calibration process of the high-order mathematical model parameters is constructed as an inverse optimization problem, training sample data is collected, and an optimization algorithm is used to solve the constructed optimization problem to complete the calibration of the high-order mathematical model parameters. The high-order mathematical modeling method improves the accuracy of the mathematical model, and the fitting calibration method improves the practicability of the calibration method.

[0042] As Figure 1 shown, a spectral imaging method based on high-order mathematical modeling and fitting calibration has the following specific steps:

[0043] Step 101: Use a three-dimensional tensor to represent the spectral data cube to be measured, where M, N, and L represent the sizes of the spectral data cube in the vertical spatial dimension, horizontal spatial dimension, and spectral dimension, respectively.

[0044] Step 102: As shown in the CASSI system structure diagram, use a light source to irradiate the spectral data cube X to be measured. The data cube is imaged onto the coded aperture through an imaging lens group and is encoded by the coded pattern. Use the matrix Figure 2 to represent the coded pattern loaded on the coded aperture during the s-th projection measurement. The encoded data cube is dispersed in the horizontal or vertical spatial direction by a dispersion prism. The dispersed data cube is imaged onto an array detector through a relay lens group. The detector integrates the modulated data cube along the spectral dimension to obtain the two-dimensional projection measurement image in the s-th projection measurement where V = (N + L - 1).

[0045] Step 103: Establish a high-order mathematical model for the projection measurement process in Step 102.

[0046] The specific construction method of the high-order mathematical model in Step 103 of the present invention is as follows:

[0047] Step 201: Based on the two-dimensional coded pattern construct the three-dimensional equivalent correction code in the s-th projection measurement according to the formula where the three-dimensional tensor represents the response non-uniformity. The subscript · represents the two-dimensional matrix corresponding to the k-th band in the three-dimensional tensor. ⊙ represents the Hadamard product operation. The equivalent correction code C k represents the response non-uniformity while performing the encoding operation. If there is no obvious response non-uniformity during the projection measurement process, then R is a tensor with all element values being 1, and the equivalent correction code degenerates into the original coded pattern at this time. s

[0048] Step 202: Divide the spatial region and band groups in the spectral data cube. Use U, V, and W to represent the total number of divisions of the spectral data cube in the vertical spatial dimension, horizontal spatial dimension, and spectral dimension, respectively. Use X u，v，w to represent the subset of the spectral data cube obtained by the division. Then X = {X u，v，w ​​}, where u, v, and w respectively represent the indices for the vertical region of space, the horizontal region of space, and the band group. Each set of index values points to a subset of the data cube, with u ∈ {1, 2,..., U}, v ∈ {1, 2,..., V}, and w ∈ {1, 2,..., W}. If each band is separately divided into a band group, then this subset is a matrix; otherwise, this subset is a three-dimensional tensor.

[0049] Step 203: Set one or more sets of two-dimensional convolutional kernels before and after the equivalent correction encoding. Let B represent the number of sets of convolutional kernels, and let W b represent one set of convolutional kernels among them, where b ∈ {1, 2,..., B}, and let represent a single two-dimensional convolutional kernel in each set of convolutional kernels. Then is used to perform three-dimensional convolutional operations on the data corresponding to the data cube subset X u，v，w respectively. The convolutional kernels comprehensively represent various non-ideal factors such as the aberration of the imaging lens group, the aberration of the relay lens group, the pixel mismatch between the coded aperture and the detector, the high-order dispersion of the prism, and the assembly error of the system.

[0050] Step 204: Combine the dispersion and integration effects in the projection measurement process to construct a high-order mathematical model containing the equivalent correction encoding and one or more sets of two-dimensional convolutional kernels. Taking the case of three sets of convolutional kernels as an example, the high-order mathematical model is constructed as where represents the three-dimensional convolutional operation for spatial sub-regions and band sub-groups. For the data corresponding to all data cube subsets, according to the indices of u, v, and w, perform three-dimensional convolutional operations with the two-dimensional convolutional kernels of the same index respectively. If there is a set of convolutional kernels composed only of unit impulse functions, the high-order mathematical model of the three sets of convolutional kernels degenerates into the high-order mathematical model of two or one set of convolutional kernels. If in Step 202 each band of the data cube is separately divided into a band group, the three-dimensional convolutional operation at this time degenerates into a two-dimensional convolutional operation. disp(·) represents the dispersion operation. Taking the matrix corresponding to the first band of the three-dimensional tensor as a reference, shift the matrices corresponding to the remaining bands one unit in the positive direction of the spatial horizontal or vertical dimension relative to the matrix of its previous band. ∑ k · represents the cumulative operation along the spectral dimension, and accumulates the spectral data at each spatial position of the three-dimensional tensor. The operation result is a matrix.

[0051] Step 104: Collect the spectral data cube and projection measurement images of the sample target as training sample data, then construct an optimization problem, solve the optimization problem through a numerical algorithm, fit the training sample data, and calibrate the parameters of the high-order mathematical model.

[0052] The specific method for calibrating the parameters of the high-order mathematical model in Step 104 is as follows:

[0053] Step 301: Establish a low-order mathematical model \(Z\) corresponding to the \(p\)th projection measurement p =\(\sum\) k disp(X⊙E p ), where \(Z\) p is the projection measurement image calculated by the low-order mathematical model. Vectorize and scan the matrix \(Z\) p and convert \(Z\) p into a column vector \(z\) p . Use \(G\) p to represent the projection matrix corresponding to the \(p\)th low-order projection measurement process. Then the matrix low-order mathematical model for a total of \(P\) projection measurements is \(z = Gx = G\Psi\theta\), where \(z = [z 1T , z 2T , …, z PT T , \(G = [G 1T , G 2T , …, G PT T , and \(\Psi\) represents the sparse basis.

[0054] Step 302: Use a standard whiteboard as the sample target. Change the encoding pattern \(E_p\) and take \(P\) projection measurement images under the illumination of the light source used during imaging Perform vectorized scanning on to obtain the corresponding column vector According to the formula and Use to represent all \(P\) projection measurement images obtained by shooting, and construct an optimization problem where \(\tau\) represents the weight factor, and \(\|\cdot\|\) t represents the calculation of the \(l\) t norm.

[0055] Step 303: Use the compressive sensing reconstruction algorithm to solve for \(r\) and arrange \(r\) into a three-dimensional tensor Calibrate to obtain the response non-uniformity parameter.

[0056] Step 304: On the basis of the known response non-uniformity parameter \(R\), use a target with certain texture features as the sample target. Change the encoding pattern \(E_q\) and take \(Q\) projection measurement images under the illumination of the light source used during imaging If the spectral data cube of the sample target is unknown, use other spectral imaging systems to obtain the spectral data cube of the sample target If the spectral data cube of the sample target is known, use it directly, and then construct an optimization problem where \(Y q is the projection measurement image calculated by the high-order mathematical model in the \(q\)th projection measurement, and \(C q ​​is the equivalent correction coding in the q-th projection measurement.

[0057] Step 305: Solve for W using an optimization algorithm 1 , W 2 , …, W B , and calibrate to obtain a set or multiple sets of two-dimensional convolution kernel parameters that are set.

[0058] Step 105: Perform a vectorized scan on the projection measurement image Ys, and convert the matrix Y s into a column vector Perform a vectorized scan on the data cube X, and convert the three-dimensional tensor X into a column vector Use the matrix to represent the projection matrix corresponding to the projection measurement process. According to the high-order mathematical model calibrated in step 104, establish the matrix high-order mathematical model for the s-th projection measurement process as y s = H s x.

[0059] Step 106: Change the encoding pattern E s , repeat the content described in step 102, perform a total of S projection measurements. All projection measurement images can be represented by the vector y = [y 1T , y 2T , …, y ST T , where · T represents the transpose operation. Substitute the changed encoding patterns into the high-order mathematical model calibrated in step 104 in sequence, and repeat the content described in step 105. The matrix high-order mathematical model for a total of S projection measurements is y = Hx, where H = [H 1T , H 2T , …, H ST T .

[0060] Step 107: Since there is a high correlation between adjacent bands and adjacent spatial positions in the spectral data cube, therefore, assume that x is sparsely represented on the sparse basis Ψ, θ represents the sparse coefficient. According to the principle of compressive sensing, y = Hx can be represented by the underdetermined equation y = HΨθ.

[0061] Step 108: Use a reconstruction algorithm to solve the underdetermined equation y = HΨθ established in step 107 and recover the spectral data cube X.

[0062] Embodiment of the present invention:

[0063] Figure 3 It is a comparison of the reconstruction results of the high-order mathematical model and the fitting calibration method, the mathematical model based on actual shooting coding and its calibration method under the same number of projection measurements and the same reconstruction algorithm used. Figure 3 ​​(a) is the ground truth of the spectral image, Figure 3 (b) is the spectral image reconstructed using the mathematical model based on real-shot encoding and its calibration method, Figure 3 (c) is the spectral image reconstructed using the high-order mathematical model and the fitting calibration method. The central wavelength of the band is marked at the lower left corner of each band image in each group of spectral images, and the peak signal-to-noise ratio (PSNR) is marked at the upper right corner. Compared with the mathematical model based on real-shot encoding, the high-order mathematical model has higher reconstruction accuracy.

[0064] Figure 4 It is a sampling comparison of the reconstructed spectral images of the two mathematical models and their calibration methods. Figure 4 (a) is the grayscale image of the imaging target object, on which three sampling positions of the spectral curve are marked, Figure 4 (b)-(d) show the spectral curves detected by different systems at the three sampling positions. The triangular curve represents the ground truth of the spectral curve, the circular curve represents the spectral curve reconstructed using the mathematical model based on real-shot encoding and its calibration method, and the square curve represents the spectral curve reconstructed using the high-order mathematical model and the fitting calibration method. For the reconstruction result corresponding to the mathematical model based on real-shot encoding, the characteristics of the target spectral curve are not reconstructed at all. The reconstruction result corresponding to the high-order mathematical model is basically consistent with the reference curve.

[0065] Although the specific embodiments of the present invention have been described in conjunction with the accompanying drawings, for those skilled in the art, without departing from the principle of the present invention, several modifications, substitutions and improvements can still be made, which should also be regarded as falling within the protection scope of the present invention.

Claims

1. A spectral imaging method based on high-order mathematical modeling and fitting calibration, characterized in that The steps of the method include: In the first step, a light source is used to irradiate the spectral data cube to be measured. The spectral data cube is successively imaged by an imaging mirror, encoded by an encoding aperture, dispersed by a dispersion prism, imaged again by a relay mirror, and integrated by an array detector to complete the projection measurement process, and a two-dimensional projection measurement image in the projection measurement is obtained; In the second step, a high-order mathematical model is established for the projection measurement process in the first step; In the third step, the spectral data cube and the projection measurement image of the calibration target are collected as training sample data, an optimization problem is constructed, and then the parameters of the high-order mathematical model established in the second step are calibrated; In the fourth step, the two-dimensional projection measurement image obtained in the first step is vectorized and scanned to be converted into a column vector, and the spectral data cube to be measured is vectorized and scanned to be converted into a column vector. According to the high-order mathematical model calibrated in the third step, a matrix high-order mathematical model corresponding to one projection measurement process is established; In the fifth step, the encoding pattern is changed, the content described in the first step is repeated, multiple projection measurements are performed, the changed encoding patterns are successively substituted into the high-order mathematical model calibrated in the third step, and the content described in the fourth step is repeated to establish a matrix high-order mathematical model describing all projection measurement processes; In the sixth step, a reconstruction algorithm is used to restore the spectral data cube according to the projection measurement data obtained in the fifth step and the established matrix high-order mathematical model; The high-order mathematical modeling refers to establishing a high-order mathematical model. In the second step, the method for establishing the high-order mathematical model is: Step 201, based on the two-dimensional coding pattern According to the formula Constructing the three-dimensional equivalent correction code in the s-th projection measurement The three-dimensional tensor represents the response inhomogeneity, subscript · k represents the two-dimensional matrix corresponding to the kth band in the three-dimensional tensor, ⊙ represents the Hadamard product operation, and the equivalent correction code C s While implementing the coding operation, the response inhomogeneity is characterized. If there is no response inhomogeneity in the projection measurement process, R is a tensor whose element values ​​are all 1, and the equivalent correction code at this time degenerates into the original coding pattern; Step 202, divide the spatial region and band groups in the spectral data cube. Use U, V, and W to represent the total number of divisions of the spectral data cube in the vertical spatial dimension, horizontal spatial dimension, and spectral dimension respectively. Use X u,v,w to represent the subset of the spectral data cube obtained by the division. Then X = {X u,v,w}, where u, v, and w represent the indices of the vertical spatial region, horizontal spatial region, and band group respectively. Each group of index values points to a subset of the data cube, u ∈ {1, 2, …, U}, v ∈ {1, 2, …, V}, w ∈ {1, 2, …, W}. If each band is separately divided into a band group, then this subset is a matrix; otherwise, this subset is a three-dimensional tensor. Step 203, before and after the equivalent correction coding, set one or more groups of two-dimensional convolutional kernels. Let B denote the number of groups of convolutional kernels, and let W b represent one group of convolutional kernels among them, where b ∈ {1, 2, …, B}, and let represent a single two-dimensional convolutional kernel in each group of convolutional kernels. Then is used to perform three-dimensional convolutional operations on the data corresponding to the data cube subset X u,v,w respectively; Step 204, construct a high-order mathematical model containing an equivalent correction code and one or more groups of two-dimensional convolution kernels.

2. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 1, characterized in that: In the first step described above, the spectral data cube to be measured is irradiated with a light source The spectral data cube X is imaged onto the coded aperture through an imaging mirror assembly, encoded by the coded pattern, and the matrix represents the coded pattern loaded on the coded aperture during the s-th projection measurement. The encoded spectral data cube X is dispersed in the horizontal or vertical direction in space by a dispersive prism. After dispersion, the spectral data cube X is imaged onto an array detector through a relay mirror assembly. The array detector integrates the spectral data cube X along the spectral dimension to obtain a two-dimensional projection measurement image in the s-th projection measurement where V = (N + L - 1).

3. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 1 or 2, characterized in that: In the third step, the spectral data cube and the projection measurement image of the calibration target are collected as training sample data, then an optimization problem is constructed, the optimization problem is solved by a numerical algorithm, the training sample data is fitted, and the parameters of the high-order mathematical model established in the second step are calibrated.

4. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 3, characterized in that: In the fourth step described above, for the two-dimensional projection measurement image Y obtained in the first step s perform vectorized scanning on the two-dimensional projection measurement image Y s and convert it into a column vector Perform vectorized scanning on the spectral data cube X to be measured, and convert the three-dimensional tensor X into a column vector Use the matrix to represent the projection matrix corresponding to the s-th projection measurement process. According to the high-order mathematical model calibrated in the third step, establish the matrix high-order mathematical model y corresponding to the s-th projection measurement process s = H s s.

5. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 4, characterized in that: In the fifth step described above, change the encoding pattern E s , repeat the content described in the first step, perform a total of S projection measurements, and represent all projection measurement images with the vector y = [y 1T , y 2T , …, y ST T . Here, T represents the transpose operation. Substitute the changed encoding patterns into the high-order mathematical model calibrated in the third step in sequence, repeat the content described in the fourth step, and establish the matrix high-order mathematical model y = Hx corresponding to a total of S projection measurements, where H = [H 1T , H 2T , …, H ST T .​​ 6. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 5, characterized in that: Assume that x is sparsely represented on the sparse basis Ψ, θ represents the sparse coefficient, and according to the principle of compressive sensing, y = Hx is represented by the underdetermined equation y = HΨθ.

7. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 1, characterized in that: In the first step described above, the spectral data cube to be measured is represented by a three-dimensional tensor where M, N, and L represent the sizes of the spectral data cube in the vertical dimension, horizontal dimension, and spectral dimension of space, respectively.

8. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 1, characterized in that: In the said step 204, the high-order mathematical models of the three groups of convolution kernels are constructed as where represents a three-dimensional convolution operation of spatial sub-regions and band sub-groups. The data corresponding to all data cube subsets are respectively subjected to a three-dimensional convolution operation with two-dimensional convolution kernels of the same index according to the indices of u, v, and w. If there is a convolution kernel group composed only of unit impulse functions, the high-order mathematical models of the three groups of convolution kernels degenerate into the high-order mathematical models of two groups or one group of convolution kernels. If step 202 separately divides each band of the data cube into a band sub-group, the three-dimensional convolution operation at this time degenerates into a two-dimensional convolution operation. disp(·) represents a dispersion operation. Based on the matrix corresponding to the first band of the three-dimensional tensor, the matrices corresponding to the remaining bands are respectively translated one unit in the positive direction of the spatial horizontal or vertical dimension relative to the matrix of its previous band. ∑ k · represents an accumulation operation along the spectral dimension, accumulating the spectral data at each spatial position of the three-dimensional tensor. The operation result is a matrix.

9. The spectral imaging method based on high-order mathematical modeling and fitting calibration according to claim 1, characterized in that: In the third step, the specific calibration method for the high-order mathematical model parameters is: Step 301, establish a low-order mathematical model \(Z\) corresponding to the \(p\)th projection measurement p =\(\sum\) k \(disp(X\odot E\) p ), where \(Z\) p is the projection measurement image calculated by the low-order mathematical model. Perform vectorized scanning on the matrix \(Z\) p and convert \(Z\) p into a column vector \(z\) p . Let \(G\) p represent the projection matrix corresponding to the \(p\)th low-order projection measurement process. Then the matrix low-order mathematical model for a total of \(P\) projection measurements is \(z = Gx = G\varPsi\theta\), where \(z=[z\) 1T ,z\) 2T ,\(\cdots,z\) PT T , \(G=[G\) 1T ,G\) 2T ,\(\cdots,G\) PT T , and \(\varPsi\) represents the sparse basis;​​ Step 302: Use a standard whiteboard as the sample target object and change the encoded pattern E p , and under the illumination of the light source used during imaging, capture P projection measurement images Perform vectorized scanning to obtain the corresponding column vectors According to the formula and Use to represent all P projection measurement images captured, and construct an optimization problem where τ represents the weight factor, and ‖·‖ t represents the calculation of the l t norm; Step 303: Solve for r using a compressive sensing reconstruction algorithm and arrange r into a three-dimensional tensor Calibrate to obtain the response non-uniformity parameter; Step 304, on the basis of the known response non-uniformity parameter R, use the object with texture features as the sample object and change the encoding pattern E q , under the illumination of the light source used during imaging, capture Q projection measurement images If the spectral data cube of the sample object is unknown, use other spectral imaging systems to obtain the spectral data cube of the sample object If the spectral data cube of the sample object is known, directly use it and then construct an optimization problem where Y q is the projection measurement image calculated by the high-order mathematical model in the q-th projection measurement, and C q is the equivalent correction code in the q-th projection measurement; Step 305, solve for W using an optimization algorithm 1 , W 2 , …, W B , and calibrate to obtain the set of one or more sets of two-dimensional convolution kernel parameters set.

Images